필터 지우기
필터 지우기

Problem with Parabolic Linear PDE-crank-nicolson

조회 수: 3 (최근 30일)
doliph
doliph 2013년 11월 10일
답변: zanubah adillah rahmah 2023년 6월 3일
Hi!
I am using crank nicolson method to implicitly solve a mass diffusion equation. Where a gas concentration above a 10cm column of water is held at C=.01mol/m3. The Diffusion constant is D=2*10^-9m2/s. Boundary conditions are assumed to be C=0 @ L=10cm. I am getting a resultant plot and solution but it doesn't necessarily make sense with what you would expect the results to be.
Any suggestions would be greatly appreciated!
function Crank
% Parameters needed to solve the equation via Crank-Nicholson method
timesteps = 100; % Number of time steps
L =10.; % characteristic length cm
dt = .01;
n = 50.; % Number of space steps
dz = L/n;
diffco = (2*10^-9)*10000*3600.*24.; % diffusion coefficient in days
xi = diffco*dt/(dz*dz); % Parameter of the method
% Initial Concentration
for i = 1:n+1
z(i) =(i-1)*dz;
u(i,1) =.01;
end
% Concentration at the boundary (C=0)
for k=1:timesteps+1
u(1,k) = 0;
u(n+1,k) = 0.;
time(k) = (k-1)*dt;
end
% Defining the Matrices M_right and M_left in the method
aL(1:n-2)=-xi;
bL(1:n-1)=2.+2.*xi;
cL(1:n-2)=-xi;
MMl=diag(bL,0)+diag(aL,-1)+diag(cL,1);
aR(1:n-2)=xi;
bR(1:n-1)=2.-2.*xi;
cR(1:n-2)=xi;
MMr=diag(bR,0)+diag(aR,-1)+diag(cR,1);
% Implementation of the Crank-Nicholson method
for k=2:timesteps % Time Loop
uu=u(2:n,k-1);
u(2:n,k)=inv(MMl)*MMr*uu;
end
% Graphical representation of the temperature at different selected times
figure(1)
plot(z,u)
title('Concentration-Crank-Nicholson method')
xlabel('Z')
ylabel('C')
%u'
figure(2)
mesh(z,time,u')
title('Concentration within the Crank-Nicholson method')
xlabel('Z')
ylabel('Concentration')
end

답변 (1개)

zanubah adillah rahmah
zanubah adillah rahmah 2023년 6월 3일
method crank nicolson delT/delt = k*del^2T/delx^2
Heat conduction in an aluminum rod long flat,
stem length L=10cm, delta x = 1 cm
time step, delta t = 0.1s
coefficient diffusion thermal, k = 0.835 cm^2/s
boundary conditions: T(x=0,t)= 100 celcius and T(x=20,t)=50 celcius
initial value: T(x,t=0)= 0 celcius

카테고리

Help CenterFile Exchange에서 Boundary Conditions에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by