Diffusion Advection Reaction Equation

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]

Use "pdepe" to get rid of these technical programming problems.

[Bright Takyi]

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]

Either you program your Crank-Nicolson or you use MATLAB's pdepe.

Answer 1

[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

Answer 2

[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