Diffusion Advection Reaction Equation

43 次查看（过去 30 天）
Raj001
Raj001 2018-7-8
评论： Torsten 2022-7-5
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
  3 个评论
显示 1更早的评论
Bright Takyi
Bright Takyi 2022-7-5
Please is it possible to solve ADRE equation with crank-nicolson finite difference method with MATLAB Pdepe? Please I will like to have a document that will guide me if it possible. Thank you
Torsten
Torsten 2022-7-5
Either you program your Crank-Nicolson or you use MATLAB's pdepe.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Hassan Abdullah Saleem
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
azertazert azertazertazertazert
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

类别

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by