Write script matlab for solving the convection-diffusion equation of a compound in a 3D annular reactor (diffusion in two directions with convection)

hello,

I would like help to write a matlab for solving the convection-diffusion equation of a compound in a 3D annular reactor (diffusion in two directions r and teta with convective transport in the main direction, see figure below) with a reaction on the surface of inner cylinder. The elements of the problem are as follows:

a 3D model is proposed. This model is based on a convection / radial diffusion principle in the annular space of the reactor. The mass balance in the gas phase can be written (eq. (1)):

(1)

0 < z < L

Inlet concentration : C0=4.5 10-5 mol.m-3

Where r is radial coordinate (m), z is axial coordinate, u is flow velocity (m s-1) and C, is the concentration of NCl3 (mol m-3).

The model takes the velocity profile developed in the annular space into account, as expressed by the following relationship:

(2)

With

(3)

z=1m

Um=0.5m s-1

ri=0.04m, the internal radius

re=0.05m, the external radius of reactor (m).

Boundary conditions are:

• at the inner surface of the external cylinder, r=re:

(4)

• at the catalyst surface, r=ri,

(5)

Where rp is the photocatalytic rate (mol m-2 s-1).

(6)

I=9.8 W m-2, n =1, kp = 2.6x10-6 mol s-1 W-1and Kp = 3x103m3 mol-1.

and for and

Thank you very much for your help

Answer 1

[Torsten]

Since you don't seem to have experience with discretization and the method of lines, I'd strongly recommend a ready-to-use 3d solver for partial differential equations like the MATLAB PDE toolbox or COMSOL Multiphysics.

[Fabien]

Thank's a lot for your answer!

Regards

[Torsten]

I forgot to say that you should check whether the solvers can work in cylindrical coordinates. But since the diffusion coefficients in r- and theta- direction seem to be the same in your case, you could equally well work in cartesian coordinates, I guess.

[Fabien]

Dear Torsten,

Thank you for your comment. I agree with you, but in my case, there is a surface reaction (rp) with irradiance term that depends on theta!

Regards

F

[Torsten]

I agree with you, but in my case, there is a surface reaction (rp) with irradiance term that depends on theta!

This does not matter - you can work in cartesian coordinates. If necessary, use "cart2pol":

https://uk.mathworks.com/help/matlab/ref/cart2pol.html

Even for the diffusion coefficients: if they were given in (r,theta,z)-direction and different in each direction, you could have transformed them to diffusion coefficients in (x,y,z)-direction, I guess.

[Fabien]

Dear Torsten,

Thank you for you answer, i wrote the following script from your advices but i did't achieve a solution with a error message such as:

Error using pde.internal.FESolverUtilities/solveStationaryNonlinear

Step size too small.

Error in pde.PDEModel/solvepde (line 79)

u = self.solveStationaryNonlinear(coefstruct, u0);

Error in Reacteur_annulaire_solaire (line 34)

results = solvepde(model);

Could you help me please

Kind regards

Fabien

%%%%%%%%%%%%%
%Convection diffusion reaction model by finite elements
close all;
clear all;
global cc0 Da ks ki Ka zlist U
ri=0.045; %inner radius of annular space m
re=0.05; %external radius of annular space m
s= pi*(re^2-ri^2);%annular area m2
Da =9.4e-6;% diffusion coefficient m2/s
qair=8;% inlet flow rate of the reactor L/min
%**********************Case 1 - Q=8L/min - C0=1.13.x10-4 mol/m3 - 9.2 W/m2********************
U=(qair/60/1000)/s;%inlet air velocity of the reactor 
a=2.6e-6; % kinetic constant LH model  mol.m2/m3/s/W 
ki=9.2;% irradiance W/m2 
ks=ki*a;%calcul of kinetic constant with irradiance
Ka=3e4;% adsorption constant LH model 
cc0=1.13e-4; %inlet concentration  
model = createpde;
importGeometry(model,"Reacteur2.stl");% import geometry
applyBoundaryCondition(model,"dirichlet",Face=1,u=cc0);% F1= inlet annular face of the reactor
applyBoundaryCondition(model,"neumann",Face=2,g=@gcoef);% F2= inner cylinder surface with catalytic reaction surface 
applyBoundaryCondition(model,"neumann",Face=3,g=0);% F3= external cylinder surface  
%applyBoundaryCondition(model,"dirichlet",Face=4,u=50);% F4= outlet annular face of the reactor
zlist = -500:10:500; % step of z axis
nr=length(zlist);
specifyCoefficients(model,m=0,d=0,c=Da,a=0,f=@fcoef);
setInitialConditions(model,U,0);
pdegplot(model.Geometry, ...
    FaceLabels="on", ...
    FaceAlpha=0.5)
generateMesh(model);%,'Hmax', h/2);
results = solvepde(model);%,zlist);
u = results.NodalSolution;
msh = results.Mesh;
figure
pdeplot3D(model,ColorMapData=u)
sol = results.NodalSolution;
figure (3)
A= findElements(msh,"region",Cell=1);
pdemesh(msh,ElementLabels="on")
%%Boundary Condition (at the inner cylinder surface)
function g = gcoef(location,state)
global ks Ka
g  = -ks*Ka*state.u/(1+Ka*state.u);
end
%  2d member of the convection/diffusion equation
function f = fcoef(location,state)
global U  
f=U*state.uz;
end 

[Torsten]

The z-coordinate in your model takes the role of the t-coordinate in the PDE Toolbox. So I guess you have to set d = U and f = 0 in specifyCoefficients ? And perform a "time-dependent" (z-dependent) instead of a stationary simulation ? And change your geometry from 3d to 2d - thus only do a simulation in a circular ring ?

[Fabien]

Dear Torsten,

Thank you very much for your help, I have to apply your suggestions right now.

Kind regards.

F