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 plot Figure 6 of Ref.2 and compare it to Figure 2 of Ref. 1. Both Figures give the dimensionless bulk concentration versus the Damköhler group I when we consider a laminar flow reactor and a non-Newtonian fluid (the Prandtl-Eyring fluid).
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 a Prandtl-Eyring fluid
V[r_]:=1-Cosh[br]/Cosh[b]
In[]:=
BETA=1/w;w=1-2Tanh[b]/b+2/b^2-2/b^2/Cosh[b]
Out[]=
1+--
2
2
b
2Sech[b]
2
b
2Tanh[b]
b
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
In[]:=
BC1=Table[Ca[[1+(j-1)(Np+1)]]==1,{j,1,Np+1}];
In[]:=
BC2=Table[L1r[[i+Np(Np+1)]].Ca==0,{i,2,Np}];
In[]:=
BC3=Table[L1r[[i]].Ca==0,{i,2,Np}];
In[]:=
BC4=Table[L1z[[Np+1+(j-1)(Np+1)]].Ca==0,{j,1,Np+1}];
Unknown concentrations
In[]:=
unknowns=Join[Flatten[Table[{{Ca[[i]],0.25}},{i,1,(Np+1)(Np+1)}],1]];
In[]:=
Length[unknowns]
Out[]=
1681
System of equations to be solved using FindRoot
In[]:=
sys=Flatten[Join[{BC1,BC2,BC3,BC4,Flatten[Table[{Eq1[i+(j-1)(Np+1)]},{i,2,Np},{j,2,Np}],1]}],1];
In[]:=
Length[sys]
Out[]=
1681
Parameters KAPPA is 0.01 - same as Fig. 2 of Ref. 1
In[]:=
KAPPA=0.01;
Plot for b=1.0
b=1.0;
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[]:=
MCUBC[z_?NumericQ]:=NIntegrate[rV[r]f04[r,z],{r,0,1}]/NIntegrate[rV[r],{r,0,1}];
In[]:=
p1=LogPlot[MCUBC[z],{z,0,L},PlotRange->{{0,7},{0.01,1}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Red,Dashed}];
Plot for b=3.0
b=3.0;
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[]:=
MCUBC[z_?NumericQ]:=NIntegrate[rV[r]f02[r,z],{r,0,1}]/NIntegrate[rV[r],{r,0,1}];
In[]:=
p2=LogPlot[MCUBC[z],{z,0,L},PlotRange->{{0,7},{0.01,1}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Green,Dashed}];
Plot for b=5.0
b=5.0;
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[]:=
MCUBC[z_?NumericQ]:=NIntegrate[rV[r]f08[r,z],{r,0,1}]/NIntegrate[rV[r],{r,0,1}];
In[]:=
p3=LogPlot[MCUBC[z],{z,0,L},PlotRange->{{0,7},{0.01,1}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Blue,Dashed}];
Plot for b=10.0
b=10.0;
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[]:=
MCUBC[z_?NumericQ]:=NIntegrate[rV[r]f06[r,z],{r,0,1}]/NIntegrate[rV[r],{r,0,1}];
In[]:=
p4=LogPlot[MCUBC[z],{z,0,L},PlotRange->{{0,7},{0.01,1}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Magenta,Dashed}];
Plot for b=1000.0
b=1000;
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[]:=
MCUBC[z_?NumericQ]:=NIntegrate[rV[r]f0[r,z],{r,0,1}]/NIntegrate[rV[r],{r,0,1}];
In[]:=
p5=LogPlot[MCUBC[z],{z,0,L},PlotRange->{{0,7},{0.01,1}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Cyan,Dashed}];
Plot for b=0.001
b=0.001;
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];f1=Interpolation[a,InterpolationOrder3];
In[]:=
MCUBC[z_?NumericQ]:=NIntegrate[rV[r]f1[r,z],{r,0,1}]/NIntegrate[rV[r],{r,0,1}];
In[]:=
p6=LogPlot[MCUBC[z],{z,0,L},PlotRange->{{0,7},{0.01,1}},Frame->True,GridLines->Automatic,AspectRatio->1,PlotStyle->{Thick,Darker@Orange,Dashed}];
In[]:=
Showp1,p2,p3,p4,p5,p6,ImageSize600{1,1},FrameLabel->Style["Domköhler Group I",16],Style" bulk (or cup mixing) concentration",16
Out[]=
Superposing our results and Fig. 2 of Ref. 1
Superposing our results and Fig. 2 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 Prandtl–Eyring fluid
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/3274265
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/3274265