ABSTRACT(original article): The laminar flow reactor (LFR) is one of the most comprehensive problems in chemical reaction engineering, as its modeling involves mass, heat, and momentum conservation equations coupled with chemical reaction rate equations. It is relatively easy to grasp its basic operation, but the solution of the problem is far from trivial. Although there are analytical solutions that simplify the problem, students and tutors must use efficient numerical strategies to appropriately solve some LFR problems. In the present investigation, we solve several different cases of two-dimensional cylindrical LFR, beginning with the isothermal case as a benchmark. Calculations were performed using the Chebyshev orthogonal collocation technique by custom scripts in Scilab, Mathematica©, and Matlab®, and were compared with solutions available from the ECRE version of COMSOL®, which was used as the reference, as well as available analytical solutions in the literature. After analyzing the case of the isothermal LFR with Newtonian fluid, we explored non-Newtonian fluids, including Carreau and Bingham fluids, whose LFR results are not available in the preceding literature. For Newtonian fluids, besides (a) the isothermal case, we also explore three other non-isothermal design cases: (b) cooling jacket with a fixed heat exchange coefficient at the wall, (c) adiabatic operation, (d) non-isothermal LFR with isothermal wall. Through a performance criterion, the different operation models are compared and we show that non-isothermal design cases perform better than the isothermal case. Computation times for different scenarios are quite short, and taking 25 nodal points or more suffices for an accurate and timely solution of the more complex problem, Case (b). The numerical modeling approach is useful from a pedagogical standpoint, as one can compare numerical results with classical assumptions, and progress from a more restrictive conceptual model (segregated flow) to an array of different kinds of operation with the LFR, in which one needs to consider diffusion and temperature distribution. CITATION (original article): A.-H. Soares, L. de Sousa Santos, H. Binous, A. A. Shaikh, and A. Bellagi, Exploring tubular steady-state laminar flow reactors with orthogonal collocation, Comput. Appl. Eng. Educ. (2024), e22795. DOI: https://doi.org/10.1002/cae.22795
We reproduce Figure 4 of Ref. 1 using Chebyshev orthogonal collocation.
Number of collocation points is Np+1
Np=40;
Differentiation matrix
Ds=Table[d[i,j],{i,0,Np},{j,0,Np}];d[0,0]=(2Np^2+1)/6//N;d[Np,Np]=-(2Np^2+1)/6//N;d[i_,i_]:=-y[i]/(2(1-y[i]^2));y[i_]:=Cos[iPi/Np]//N;d[j_,k_]:=If[j!=k,c[j](-1)^(j+k)/(c[k](y[j]-y[k]))];c[i_]:=If[i==0||i==Np,2,1];Dsq=Ds.Ds;
Reactor dimensions
R=1;L=10.;
Partial derivative operators
L2r=(2/R)^2KroneckerProduct[Dsq,IdentityMatrix[Np+1]];L1r=(2/R)KroneckerProduct[Ds,IdentityMatrix[Np+1]];L1z=(2/L)KroneckerProduct[IdentityMatrix[Np+1],Ds];
Concentrations
Ca=Table[ca[i+(j-1)(Np+1)],{j,1,Np+1},{i,1,Np+1}]//Flatten;
Collocation point positions:
Radial positions
In[]:=
r=Table[R/2(1+y[i]),{i,Np,0,-1}]
Out[]=
{0.,0.00154133,0.00615583,0.013815,0.0244717,0.0380602,0.0544967,0.0736799,0.0954915,0.119797,0.146447,0.175276,0.206107,0.238751,0.273005,0.308658,0.345492,0.383277,0.421783,0.46077,0.5,0.53923,0.578217,0.616723,0.654508,0.691342,0.726995,0.761249,0.793893,0.824724,0.853553,0.880203,0.904508,0.92632,0.945503,0.96194,0.975528,0.986185,0.993844,0.998459,1.}
Axial positions
In[]:=
z=Table[L/2(1+y[i]),{i,Np,0,-1}]
Out[]=
{0.,0.0154133,0.0615583,0.13815,0.244717,0.380602,0.544967,0.736799,0.954915,1.19797,1.46447,1.75276,2.06107,2.38751,2.73005,3.08658,3.45492,3.83277,4.21783,4.6077,5.,5.3923,5.78217,6.16723,6.54508,6.91342,7.26995,7.61249,7.93893,8.24724,8.53553,8.80203,9.04508,9.2632,9.45503,9.6194,9.75528,9.86185,9.93844,9.98459,10.}
Axial velocity for power-law fluids
V[r_]:=1-(r/R)^((n+1)/n)
Governing equations
Eq1[i_]:=KAPPA(-1/r[[Quotient[i-1,Np+1]+1]]L1r[[i]].Ca+L2r[[i]].Ca)+BETAV[r[[Quotient[i-1,Np+1]+1]]]L1z[[i]].Ca-ra[i]0;
Rate equation - 1st order
ra[i_]:=Ca[[i]]
Boundary conditions
BC1=Table[Ca[[1+(j-1)(Np+1)]]==1,{j,1,Np+1}];BC2=Table[L1r[[i+Np(Np+1)]].Ca==0,{i,2,Np}];BC3=Table[L1r[[i]].Ca==0,{i,2,Np}];BC4=Table[L1z[[Np+1+(j-1)(Np+1)]].Ca==0,{j,1,Np+1}];
Unknown concentrations
unknowns=Join[Flatten[Table[{{Ca[[i]],0.25}},{i,1,(Np+1)(Np+1)}],1]];
System of equations to be solved using FindRoot
sys=Flatten[Join[{BC1,BC2,BC3,BC4,Flatten[Table[{Eq1[i+(j-1)(Np+1)]},{i,2,Np},{j,2,Np}],1]}],1];
Check : the number of equations matches the number of unknowns
In[]:=
Length[unknowns]
Out[]=
1681
In[]:=
Length[sys]
Out[]=
1681
Values for KAPPA and BETA
In[]:=
KAPPA=0.01;
Expression of BETA is specific to power-law fluids
BETA=(3n+1)/(n+1);
In[]:=
n=0.4;
In[]:=
sol=FindRoot[sys,unknowns];
In[]:=
a=Flatten[Table[{r[[j]],z[[i]],Ca[[i+(j-1)(Np+1)]]/.sol},{j,1,Np+1},{i,1,Np+1}],1];f04=Interpolation[a,InterpolationOrder3];
In[]:=
p1=Plot[f04[r,5],{r,0,1},PlotRange->{{0,1},{0.0,0.08}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Red,Dashed}];
In[]:=
n=0.2;
In[]:=
sol=FindRoot[sys,unknowns];
In[]:=
a=Flatten[Table[{r[[j]],z[[i]],Ca[[i+(j-1)(Np+1)]]/.sol},{j,1,Np+1},{i,1,Np+1}],1];f02=Interpolation[a,InterpolationOrder3];
In[]:=
p2=Plot[f02[r,5],{r,0,1},PlotRange->{{0,1},{0.0,0.08}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Green,Dashed}];
In[]:=
n=0.8;
In[]:=
sol=FindRoot[sys,unknowns];
In[]:=
a=Flatten[Table[{r[[j]],z[[i]],Ca[[i+(j-1)(Np+1)]]/.sol},{j,1,Np+1},{i,1,Np+1}],1];f08=Interpolation[a,InterpolationOrder3];
In[]:=
p3=Plot[f08[r,5],{r,0,1},PlotRange->{{0,1},{0.0,0.08}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Blue,Dashed}];
In[]:=
n=0.6;
In[]:=
sol=FindRoot[sys,unknowns];
In[]:=
a=Flatten[Table[{r[[j]],z[[i]],Ca[[i+(j-1)(Np+1)]]/.sol},{j,1,Np+1},{i,1,Np+1}],1];f06=Interpolation[a,InterpolationOrder3];
In[]:=
p4=Plot[f06[r,5],{r,0,1},PlotRange->{{0,1},{0.0,0.08}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Magenta,Dashed}];
In[]:=
n=0.0001;(*ornapproximately0*)
In[]:=
sol=Quiet@FindRoot[sys,unknowns];
In[]:=
a=Flatten[Table[{r[[j]],z[[i]],Ca[[i+(j-1)(Np+1)]]/.sol},{j,1,Np+1},{i,1,Np+1}],1];f0=Interpolation[a,InterpolationOrder3];
In[]:=
p5=Plot[f0[r,5],{r,0,1},PlotRange->{{0,1},{0.0,0.08}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Cyan,Dashed}];
In[]:=
n=1.;(*Newtonianfluidcase*)
In[]:=
sol=FindRoot[sys,unknowns];
In[]:=
a=Flatten[Table[{r[[j]],z[[i]],Ca[[i+(j-1)(Np+1)]]/.sol},{j,1,Np+1},{i,1,Np+1}],1];f1=Interpolation[a,InterpolationOrder3];
In[]:=
p6=Plot[f1[r,5],{r,0,1},PlotRange->{{0,1},{0.0,0.08}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Orange,Dashed}];
Our computed figure for a dimensionless axial position ζ=5 and for κ=0.01
Our computed figure for a dimensionless axial position ζ=5 and for κ=0.01
In[]:=
Show[p1,p2,p3,p4,p5,p6,ImageSize550{1,1},FrameLabel->{Style["Dimensionless Radius ",16],Style[" Dimensionless Concentration ",16]}]
Out[]=
Superposing our results and Fig. 4 of Ref. 1
Superposing our results and Fig. 4 of Ref. 1
References
References
◼
[1] R. V. Homsy and R. D. Strohman, Diffusion and chemical reaction in a tubular reactor with non-Newtonian laminar flow, AIChE J. 17 (1971), 215–219. https://doi.org/10.1002/aic.690170142
◼
[2] A.-H. Soares, L. de Sousa Santos, H. Binous, A. A. Shaikh, and A. Bellagi, Exploring tubular steady-state laminar flow reactors with orthogonal collocation, Comput. Appl. Eng. Educ. (2024), e22795. https://doi.org/10.1002/cae.22795
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Homsy's isothermal laminar flow reactor with power-law fluids
by André von-Held Soares, Lizandro de Sousa Santos, Housam Binous, Abdullah A. Shaikh and Ahmed Bellagi
Wolfram Community, STAFF PICKS, September 14, 2024
https://community.wolfram.com/groups/-/m/t/3273698
by André von-Held Soares, Lizandro de Sousa Santos, Housam Binous, Abdullah A. Shaikh and Ahmed Bellagi
Wolfram Community, STAFF PICKS, September 14, 2024
https://community.wolfram.com/groups/-/m/t/3273698