78 views (last 30 days)
Raj001 on 8 Jul 2018
Commented: Torsten on 5 Jul 2022
I have solved the advection-diffusion-reaction (del_C/del_t)=D[del^2)_C/(del_x)^2]+U[del_C/del_x]+kC equation numerically using Matlab. I have used Crank-Nicolson method to solve the problem. Is the scheme choose is perfect for better stability?
I am facing some problems here. While using the value of U, the results seems not perfect. Solution provides result for small value of U. What is the reason?
For very small value of U, the solution at some grid points is crossing below at boundary value. what is the reason?
My MATLAB code is here.
clc
close all
clear window
% 1-D Unsteady state convection diffusion Reaction problem in cartesian co-ordinate
%Governing Equation (del_C/del_t)=D(del^2_C/del_x^2)+U(del_C/del_x)+kC
numx = 70; %number of grid points in x direction
numt = 10000; %number of time steps
L = 0.1; %Thickness of the plate (m)
dx = L/(numx - 1); %grid size in x direction
dt = 0.1; %grid size in time direction
D = 1.0e-06; %species diffusivity (m2/s)
C0 = 1.0; %species concentration at plate(kg/m3)
Cinf = 0.1; %species concentration at the surface/bulk(kg/m3)
x=0:dx:L; %vector of x positions
U=0.1; %Magnitude of the convective velosity
k=-1.67e-3; %Reaction rate constant k (1/s)
r=D*dt/(2*dx*dx);
s=U*dt/(4*dx);
Peclet=U*dx/dt; %Stability: Peclet<2
%Initial condition at time t=0 for 0<x<L
for i=1:numx
C(i,1)=C0;
end
%Boundary condition at x=0 for time t>0
for j=1:numt
C(1,j)=C0;
end
%Boundary condition at x=L for time t>0
for j=1:numt
C(numx,j)=Cinf;
end
t(1)=0;
for j=1:numt-1
t(j+1)=t(j)+dt;
end
%Using Crank-Nicolson scheme and converting the system of linear equation
....having numx-2 equations and numx unknowns into matrix form.
% The matrix that provide solution at each grid is:....
...ML*C(i,j+1)+r(j+1)=MR*C(j)-r(j)
%Tridiagonal matrix at Left side in the form of square matrix for time j+1
aal(1:numx-3)=(-r-s);
bbl(1:numx-2)=1.0+2.0*r;
ccl(1:numx-3)=(-r+s);
MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);
%Tridiagonal matrix at Right side in the form of square matrix for time j
aar(1:numx-3)=(r+s);
bbr(1:numx-2)=1.0-(2.0*r)+(k*dt);
ccr(1:numx-3)=(r-s);
MMr=diag(bbr,0)+diag(aar,-1)+diag(ccr,1);
%To convert the left side matrix into square form using the boundary
....conditions (For known values at C(1,j+1) & C(numx,j+1)
rl=zeros(numx-2,1);
for j=1:numt-1
rl(1,1)=(-r-s)*C(1,j+1);
rl(numx-2,1)=(-r+s)*C(numx,j+1);
end
rr=zeros(numx-2,1);
for j=1:numt-1
rr(1,1)=(r+s)*C(1,j);
rr(numx-2,1)=(r-s)*C(numx,j);
end
%Solving matrix to find C(i,j) at each grid points:....
.....C(i,j+1)=inv(ML)*[(MR*C(i,j))+r(j)-r(j+1)]
for j=1:numt-1
uu=C(2:numx-1,j);
C(2:numx-1,j+1)=inv(MMl)*((MMr*uu)+rr-rl);
end
%Plotting graph x/L vs C/C0 for different time
figure(1);
title('Concentration profile in x direction at different time')
hold on;
plot(x/L,C(:,1)/C0,x/L,C(:,500)/C0,x/L,C(:,1000)/C0, x/L,C(:,1500)/C0, x/L,C(:,2000)/C0,x/L,C(:,4000)/C0,x/L,C(:,5000)/C0, 'linewidth',1.5)
xlabel('x/L');
ylabel('C(x,t)/C0');
grid on
Torsten on 5 Jul 2022
Either you program your Crank-Nicolson or you use MATLAB's pdepe.

Hassan Abdullah Saleem on 12 Feb 2020
When I tried your code, it works when velocity (U) is a small value U < 0.01 ( example U = 0.0001)
Also it would be more important if you also calculate courant number in addetion to peclet number, since both in most cases must satisfiy a numerical condition.
• One other thing I think you might look for is to make sure to calculate your advection term using forward-upwind scheme and avoid central differences for the advection term.
• I have a question for you why your right hand side is a square materix? Can you please send me the full mathematical formula for your problem and the FD fomulation?.
• One other thing, I believe it would be more acurate if you use nx1 vector for the right hand side, try to use different formulation such as the general formulation that involve alpha (between 0 to 1), you can find it in many text book such as Applied Conatmenant Transport: by Zheng and Bennet
• The general formulation that has alpha with
• when alpha = 0 --> explicit
• if alpha = 1 -->implicit
• and if alpha = 0.5 --> it is Crank Nicolson

Hellow, I'm tryin to solve a problem : using discontinuous Galerkin finite elements method (DGFEM) for solving steady-state diffusion-onvection-reaction equations. the main programme in some research paper begin with
clear all
clc
% Generate the mesh
% Nodes
Nodes = [ 0 , 0 ; 0.5 , 0 ; 1 , 0 ; 0 , 0.5 ; 0.5 , 0.5 ; 1 , 0.5 ; 0 , 1 ; 0.5 , 1 ; 1 , 1 ] ;
% El ement s
Elements = [ 4 , 1 , 5 ; 1 , 2 , 5 ; 5 , 2 , 6 ; 2 , 3 , 6 ; 7 , 4 , 8 ; 4 , 5 , 8 ; 8 , 5 , 9 ; 5 , 6 , 9 ] ;
% Dirichlet bdryedge s
Dirichlet = [ 1 , 2 ; 2 , 3 ; 1 , 4 ; 3 , 6 ; 4 , 7 ; 6 , 9 ; 7 , 8 ; 8 , 9 ] ;
% Neumann bdryedge s
Neumann = [ ];
% Initial mesh struct
mesh=getmesh(Nodes,Elements,Dirichlet,Neumann,4,5) ;
------------------------------------------------------
matlab show us : ??? Undefined function or method 'getmesh' for input arguments of type 'double'.
---------------------------------------------------------------
the problem we have not the code of getmesh function