Problem with Parabolic Linear PDE-crank-nicolson
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
zanubah adillah rahmah
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Boundary Conditions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!