Problem with Parabolic Linear PDE-crank-nicolson
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
zanubah adillah rahmah
le 3 Juin 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 commentaires
Voir également
Catégories
En savoir plus sur Boundary Conditions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!