Problem with Parabolic Linear PDE-crank-nicolson
10 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
zanubah adillah rahmah
on 3 Jun 2023
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
0 Comments
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!