- Email: [email protected]
Modeling of axially symmetric flow reactors
MODELING National Measurement
OF AXIALLY
SYMMETRIC
FLOW
REACTORS
R. L. BROWN Laboratory, Center for Chemical Physics, Chemical Kinetics Division, National Bureau of Standards, Washington, DC 20234, U.S.A. (Receioed
15
April 1983)
Abstract-A method of calculating the velocity profiles and species concentrations in axially symmetric flow reactors is presented. The method is illustrated with a reactor consisting of three concentric tubes arranged so that reactants flow through the inner tube and inner annulus to a mixing region and the products flow out through the outer annulus. A single bimolecular reaction is used in the example. The method involves two steps; calculation of the velocity field from a solution of the Navier-Stokes equations, followed by solution of the species conservation equations. The technique provides a way of analyzing reactors with variable cross sections and sampling ports near mixing regions. The method used for solving the equations of motion is one originally proposed by Thorn (1933). It has been used subsequently by other workers (Lester, 1961; Mills, 1968). From the Navier-Stokes equations in terms of velocities one generates two second order coupled partial differential equations, one for the vorticity 4, and the other for the Stokes stream function I,$. These equations can be solved by standard finite difference methods. By then differentiating the stream function one can get the velocity field for subsequent use in the species continuity equations. For the analysis, incompressible flow was assumed.
INTRODUCHON
The flow tube reactor has been one of the more powerful tools for generating gas phase chemical kinetic data (Howard, 1979). It is normally used under laminar flow conditions at pressures sufficiently low to allow diffusional dispersion to minimize radial concentration gradients. At higher pressures, or in the presence of wall reactions, corrections for radial gradients must be applied to the measurements (Brown, 1978; Kaufman, 1961; Ogren, 1975; Porier & Carr, 1971). In addition, this reactor is normally operated under conditions of a fully developed parabolic velocity profile. Thus, measurements in the region of inlet tubes or in reactors with variable cross sections are avoided. It is, however, relatively easy to mathematically model more complex reactors. Furthermore, the increasing use of low cost high,speed microcomputers by many laboratories makes such modeling an attractive way to extend the range of applicability of flow reactors. The purpose of this work is to show how one can model laminar flow reactors having axial symmetry. As an example of the method a reactor consisting of three concentric tubes has been considered. The inner tube and inner annulus contain reactants in flowing carrier gases which mix and impinge on a wall oriented at right angles to their flow. This reverses their direction and forces them to flow out through the outer annulus. This particular arrangement was chosen to simulate the conditions near a mass spectrometer sampling orifice located very close to the point of mixing of the reactants. The problem is handled by splitting it into two parts. To do this, it is necessary to assume that the reactants are sufficiently buffered with inert carrier gas so that any heat released or absorbed by the reaction does not alter the temperature of the mixture enough to significantly affect the flow fields in the reactor. Then it is possible to first solve the Navier-Stokes equations of motion to get the flow pattern. Following this, the species continuity equations involving diffusion and chemical reaction can be solved using the previously calculated velocity profiles. If the gas flows are kept fixed, these same profiles can be used repeatedly to calculate concentration profiles for entirely different chemical systems. Thus, the equations of motion can be solved for a limited set of flow conditions and the velocity profiles stored for future use.
CALCULATION
OF FLOW FIELD
Derivabn of the difference equations In a cylindrical coordinate system having axial symmetry, the equations of motion for an incompressible fluid with constant viscosity are (Bird et a[., 1960)
au,
u, -
a~,
a.
+ v, -
8,
= -;$+$[;-&!$)+g]
(lb)
where V, and U, are the components of velocity in the radial and axial directions, respectively; p is the density, p is the pressure, and ,Uis the viscosity. There is also an equation for the conservation of mass which must be satisfied; this is
To set up a dimensionless introduce the variables u = V,/%>
w = V,l%,
R=r/a,
z=
z/u,
set of equations,
we
p = P l(P%% Re = av,,p 1~
(3)
where u0 is a reference velocity and a is a reference distance. These will be chosen later. 139
140
R. L. BROWN
Substituting these dimensionless variables into Eqs (la), (lb) and (2), we get the following set of equations,
Equations (8) and (9) can be solved by finite difference methods. To convert these equations into finite difference form, three point centered difference approximations ,for the derivatives were used with a square grid having a dimensionless spacing H. If G represents either q or I& then the formulas used were aG -z
G,,, +
ilR
aw
waR+waz=-z
a~
aw
a2G @= aG Fz” d2G @=
(5) We next introduce a dimensionless stream function $J by the relations
aw
au
“=E-x
~=‘f!i R
az'
vorticity
W=
9 and
(6)
_ltP.
RdR
I - G,,,- I 2H
’
Gi,j+~ - 2G, + G,>j_, ff= ’ (10)
G,+cij--t-u 2H
’
G,+,j-2Gij+Gi_,,j
HZ
.
These are approximations to the derivatives at the grid point (i, j), where the index i refers to points along the Z axis, and j refers to points along the R axis. Substituting these formulas into Eqs. (8) and (9) and solving for GLiJand v,~ gives the formulas
(7)
Substitution of Eq. (7) into Eq. (5) shows that the continuity equation is automatically satisfied. To derive an equation for the stream function, we take the derivative of U with respect to Z and that of W with respect to R in Eq. (7). Adding these and inserting Eq. (6) we get
(8) -
An equation for the vorticity can be gotten by taking the derivative of Eq. (4a) with respect to 2, subtracting the derivative of Eq. (4b) with respect to R, and inserting the vorticity wherever the expression of Eq. (6) appears. The resulting equation is
Mi,,
+ I -@t+i.j+I)(tli+l.j-
Except for a difference indices, these equations by Mills (1968).
vi-
,Jl/W4)}
in the identification of the are the same as those given
Description of reactor model The geometry of the sample reactor is shown in Fig. 1. This reactor consists of three concentric tubes with walls of finite thickness. Carrier gas containing reactant B flows through tube 1 toward the vertical wall at the origin. This tube ends at a distance z, to the right of this wall. We call the part of tube 1 to the
wa -
Fig. 1. Geometry
of the axially
(13
symmetric flow reactor used in the sample
calculation.
Modeling
of axially symmetric
right of z,, region B. Reactant A, contained in the same type of carrier gas, flows in the same direction between the outer wall of tube 1 and the inner wall of tube 2. This second tube ends at a distance z2 from the origin. The annular region between tubes 1 and 2 to the right of Z, we denote as region A. These reactants mix by diffusion in the region to the left of Z, and flow away from the origin through the annular space between the outer wall of tube 2 and the inner wall of tube 3. The portion of this space to the right of z2 will be called region C. The flows in these three regions will have parabolic velocity profiles provided one is sufficiently far away from the origin. There, the flow will have only an axial component. These limiting profiles for the three regions are (Bird er al., 1960),
141
Row reactors
sionless relation (kX5 + &I
Re@
between
the Reynolds
1 + k,, ReP)
= (1 + k&Re(C)
(20)
Let us now write dimensionless versions of Eqs. (13HlS) for the limiting velocity profiles. The reference length and velocity indicated in Eqs. (3) are taken to be a5 and {wc), respectively. We get )
x[l -(z>‘-_fAln(~)], k25SR5k35(22) WB=
WB/{ WC}= - 2
z
(ks, - k,,)
+-($1, (14)
numbers:
W, = wc,/(wc}
= 2
(23)
OSRSk,5 [I - RZ -f&
In(R)],
kd55 R 5 1. (24)
In these equations, (w”). (ws>, and (wc> are the average axial velocities, and the a, are the tube radii; the other quantities are given by k,i = q/a,
(161
1 +!c;+“&
g,=
next define
ReV
Reynolds
I=
numbers
- (% - a&
ReVI
= - alp (w,>/P
Re(C)
= G+ - a&
(17b) for
these re-
{w,>/P (18)
(wc>i~.
-w(u3+aJRe(A)
mB = - apa, Re (B) mc = np(a,
(191
+ KJ Re (C)
For conservation of mass, these three equations must add to zero. (Note that {w”} and (wa> are negative in this system.) Summing Eqs. (19) gives a dimen-
s
RWdR+O
(25)
where the integration is over the range of R values covered by a particular region, and 8 is the constant of integration. Integrating Eqs. (22)-(24) gives the limiting stream functions. lird = - $AR’[l
Consider the total mass flow rate of gas m, through any one of these limiting regions. This is given by the integral m=2~Spwrdr which is equal to 2np f wr dr since p is assumed to be constant. Because the average velocity here is 27~j wr dr /A, where A is the cross sectional area of the region, the mass flow rate can be written in terms of the average velocity as m = p (w )A. The mass flow rates can then be expressed in terms of the Reynolds numbers through use of Eq. (18). Thus, for the three limiting regions, the mass flow rates are m,=
$=-
(17a)
& = ( 1 - k$)/ln (kg) We gions;
Boundary conditions on 9 and q From the above limiting velocities far upstream and far downstream, we can derive expressions for the limiting values of the stream function. These can be calculated by integrating the expression for W in Eq. (7). We get
- #lU2
l+k,= -:BRV Gc = -
-&On
(R/k,,)
-+(R/k,$]
R*[l - iR* -j&In
R -i)J/g,,
-
+ ea + Br
$1 + 8, 126) (27) (28)
where A
=
_
2Re(AI(1 - k44 g2dWC)h - kz5)
and B = -
2Re(B)(k,,
- k,,)/Re(C)
To obtain expressions for the integration constants, we make the following observations. At any wall in the reactor the normal velocity component must vanish. Likewise the tangential component must also vanish if there is to be a no-slip condition at the walls. From Eq. (7), this means that the gradient of + must be zero at the walls. Thus, 1(1is constant along a wall. In addition, from symmetry, the radial velocity U must vanish at R = 0. Thus, from Eq. (71, + must be a constant all along the reactor axis. Since the axis hits the vertical wall at the end of tube 3, IL
R. L. BROWN
142
at the center must be the same as at the inner wall of tube 3 if we are to have continuity in # at the origin. If we arbitrarily set $ equal to zero along the axis, we get from Eq. (27) that 8, = 0. Then we can immediately get 0, from the condition I,!J,=(%= 1) = 0. The constant 0, can be obtained from the condition GA(R = k,,) = @_XR = ke5) which must be satisfied if II, is to be constant along the inner and outer surfaces of tube 2. The results are Be = (1 +J43/(2g‘S)
8, = ~~-Ak:,(l+hd
-
k:,(k:, +f45M2gM)+ 0,.
+fa/R)
‘I~ = 2BR /k f5 ‘~c =
2(2R + f.slR h’g=w
%j~t., = 3(+_, + I J -
Jl,,,Y~* - fR,vr+ IJ
(36)
For the vertical walls at the ends of tubes 1 and 2, we have %,rli., = 3(+,-i,,
- i&)/H2
- +R,?, - I.,
(37)
(29)
Having obtained values for all the mesh points on the boundaries, we can now show how Eqs. (11) and (12) were solved.
(30)
Solution
These expressions for the limiting stream functions allow us to calculate all the boundary conditions on JI. Thus, far upstream in regions A and B, Eqs. (26) and (27) aive 11 as a function of R. Far downstream, this is &en by Eq. (28). The constancy of I& at the walls and reactor axis give all the other required boundary values. The limiting expressions for the vorticity are obtained from the velocities WA, W, and WC by differentiation. We get rtA = A(2%lk:,
gives
(31) (32) (33)
These can be used to calculate the boundaries for g, far upstream and far downstream. Because U and dW/dR are zero along the axis, we have from Eq. (6) that ?I is also zero along this boundary. At the remaining boundaries, rf is not known but must be caluclated from the flow field in the region near the boundary. Lester (1961) has derived a finite difference formula for q on a boundary wall inclined at an angle cp to the z axis. If rl, is the value at the wall and n, the value at a point one mesh interval H into the flow field in a direction normal to the wall, then
x [2 + H cos cp/R, - +H’ cos2 cp/Rw2] - ‘. (34) The angle cp is measured in a counter clockwise direction and when cp = 0, the normal is taken to be in the positive R direction. Thus, for the horizontal walls corresponding to the inner surfaces of the three tubes we have rp = 180” since their normals into the fluid point in the negative R directions. For q on the outer walls of tubes 1 and 2, we have cp = 0. Thus, if qij is one of these horizontal walls, we have from Eq. (34)
where the upper signs correspond to the outer walls. For the vertical wall at the end of the reactor we have ‘p = - 90”. Therefore, for qi,j on this wall, Eq. (34)
of finite
dzxerence
equations
for
flow
field
The following values were used for the reactor parameters; u, = 0.5 cm, a2 = 0.6 cm, uj = 1.2 cm, a, = 1.3 cm and u5 = 1.7 cm. The wall thickness of the tubes was chosen to be equal to one mesh width. This produced 18 mesh points along the r axis. For cases where the wall thickness is not a multiple of the mesh width, modified finite difference formulas for the boundaries would have to be derived. We did not consider this problem. The dimensionless mesh width was H = h/a,, = 1/17. A small section of the mesh is shown in Fig. 1. The distances z, and rS of the ends of tubes 1 and 2 from the origin were taken to be 16h and 7h, respectively. The Reynolds numbers for the flow in regions A and % were both set equal to 1. Because of the relation between the Reynolds numbers, Eq. (20), this choice gives Re(C) = 0.759. For these conditions, the upstream and downstream limiting values for the axial velocities, stream functions, and vorticities are shown in Fig. 2 as a function of R. At the start of the iteration process, the limiting profiles were assigned to both @ and 9 at the upstream and downstream boundaries. The distances at which these boundaries are placed can be varied to see if their positions have any effect on the results in the mixing region. For the present conditions, distances of 46h from the origin were found to be satisfactory. At higher flows, the downstream distance for the limiting boundary would probably have to be increased. For the stream function, the limiting profile values at R = 0 and at the tube walls were continued around all of the walls. The limiting vorticity values at the walls were continued only to the ends of the tubes leaving a discontinuity there; as can be seen from Fig. 2, in the final solution, the vorticity along a wall must change continuously from the limiting value on one side to that on the other. For the vertical wall at the origin, the vorticity was initially given the average of its limiting values at R = 0 and 1. For the first run, all the internal mesh points for IJ and tf were given the value zero. In subsequent runs at different Re values, previously calculated results were used for the initial conditions. The iterative routine described by Mills (1968) was used. For this, one performs a complete iteration of the II, Eeld with Eq. (11) followed by a complete iteration of the q field with Eq. (12). In between these iterations the vorticity boundary values are calculated from Eqs. (35), (36) or (37). The gelds were scanned in the following way: starting at the upstream boundaries in regions A and B, at a fixed value of Z, R was scanned from R = 0 to the inner wall of tube 2; Z was then decreased to the next downstream mesh value
Modeling of axially symmetric fiow reactors
Fig. 2. The limiting profiles for the velocity, stream function, and vorticity in the upstream B and the downstream region C: Re(A) = Re(B) = 1, Re(C) = 0.759. and this R scan was repeated. This process was continued until Z, was reached. For the region bounded by q, a.,, and the end wall of the reactor, the scan was 6rst made for Z = 0 to Z, with R fixed at U. This 2 scan was repeated with R increased one mesh point until R reached the position of the outer wall of tube 2. Then the scan was again made over R, from the outer wall of tube 2 to the inner wall of tube 3, starting at Z = N and increasing 2 until the Iimiting downstream position in region C was reached. The comers at the ends of tubes 1 and 2 were treated according to the way suggested by Thorn (1934) as discussed by Lester (1961) and Mills (1968).
143
r+ons
_,gand
The stabiiity and convergence of the method can be estimated (Bird et ok, 1960) from the condition W s 4@/Re. Here, the Re refers to the largest Reynolds number in the system. There are reasons for limiting the method to cases where Re 5 50, and Milk (1968) gives a discussion of these. At pressures below 10 or 20 torr, this range of Reynolds numbers covers the flow conditions most often used in these reactors. For example, in a nitrogen carrier at 5 torr, the average flow velocity in tube 1 is 47 Re(A) cm/s. The results for the sample reactor are shown in Fig. 3. This is a contour diagram for the stream function JI. The direction of the flow is indicated on each
... ... ... ... -.. ... ... .-- ..* ... ... -.. ... ._. .._... .. .. . ... ... ._. .._ ... ... ... _.. ... .-. .._ ... ... ... _.. _..... .. Fig. 3. Stream line contours:
&(A
) = Re(B) = 1.
0
R. L.
144
stream line. As mentioned earlier the velocity field is obtained from the gradient of $. Normally, one can say that the closer spaced are the stream function contours the higher the fluid velocity. This is deceptive, however, for cylindrical systems because the gradient of $ must be weighted by l/R to get the velocity field. Thus, the average upstream velocity in the annulus between tubes 1 and 2 appears from the stream line density to be much greater than that in tube 1, whereas, for this case, it happends to be equal. Having shown how to calculate the velocity field, we can now consider how to calculate the species concentrations in a reacting system. DETERMINATION M
OF
A CHEMICALLY
CONCENTRATION REACTING
PROFILES SYSTEM
Let c” be the concentration of some species in the reactor. Its conservation equation is (Bird et al., 1960)
where dA is the diffusion coefficient of the species in the carrier gas and So is its production rate through chemical reactions. There will be one of these equations for each different reacting species. Introduction of the dimensionless variables of EQ. {3) gives
where D, = dJ(av,), S, = asA/(cOvO) and C, = cl/C’,. The quantity c, is some reference concentration. Since U and W are to be obtained from the previous calculation, we must here use the same reference length and velocity, namely a5 and (wc). We assume that S, is proportional to C, and write S, = PACA, where P, is some function of the species concentrations possibly including that of A. Using the derivative approximations given by Eq. (10) in Eq. (39) we get the difference formula
BROWN
gradient approximation can be used. A five point difference formula was used to calculate the concentration at a boundary. This was obtained by setting a five point formula for the derivative equal to zero. If C, is the concentration on the boundary and C,+i are the values at successive mesh points along a line perpendicular to the boundary and into the gas, then for a zero derivative at the boundary, we used the formula C, = 12(4Cb+,
- 3Cb+z + 4C,+,/3 - c&,/4)/25.
A three point formllla probably would have been satisfactory for these calculations but was not tested. Instead of the zero gradient downstream approximation another approach might be to set the axial gradient equal to S,/ W which is a reasonable approximation when S, is small. Then one could consider reactions which’had not gone as far toward completion at the downstream boundary. For the iterations, the field scanning was done in the same manner as for the $ and r] fields. To show the kind of results that can be obtained, a single bimolecular reaction was examined. One reactant entered from region A, and the other from region B. The rate function PA in this case is c,akCJu,, where k is the bimolecular rate constant for the reaction. The rate function P, is cgzkC,fv,. The flow field shown in Fig. 3 was used. Initial dimensionless concentrations were taken to be C, = 0.5 and C, = 1.7. These values were chosen so that the molar flow rates into the reactor of the two species would be the same. Dimensionless diffusion coefficients were given 7 the values D, = 0.093 and D, = 0.052. The value of the dimensionless rate constant c,ak/v, was taken to be 1.3. The results for the sohttion of Eq. (38) are shown in Figs. 4 and 5 as contour diagrams for the concentrations C,, and C,. These indicate clearly the combined effects of reaction and diffusion. If the reaction rate were zero, the limiting downstream concentrations in this case would be 0.386. Both species spread out considerably in the region to the left of the end of tube I. At these flows, there is, however, very little back diffusion into the adjacent tube by either species. The effect of the parabolic velocity profile is evident in the downstream regions in the form of higher concentrations in the middle of the flow stream than near the walls. CONCLUSIONS
For simplicity the species subscript has been dropped in this equation. For boundary conditions, we have a zero concentration gradient at the walls and perpendicular to the Z axis at R = 0. Far upstream in regions A and B the concentrations are fixed at their initial values, while far downstream we assume that the reactions have proceeded so far toward completion that a zerotThese values for the dimensionless coefficients are applicable to the diKusion in a reactor of the present dimensions of N atoms and NO in N2 at pressures around 5 torr.
The aim of this work has been to show in some detail one method of modeling an axially symmetric flow reactor containing a relatively simple chemical system. The specific example chosen for modeling has shown the type of results obtainable. These consist of concentration distributions for the individual chemical species as a function of carrier flow rates, species diffusion coefficients, and chemical rate constants. Since these calculations give the concentrations throughout the reactor, one can relate the concentration measured at any particular point to the parameters of the system. For the reactor studied here, the most obvious sampling point would be at the origin. A variety of experiments can be envisioned with this reactor. For example, one could measure diffusion coefficients by putting only one species at a time through the system, or by using non-reactive
Modeling
of axially symmetric flow reactors
Fig. 4. Concentration
contours for species
145
A.
-------
Fig. 5. Concentration species. The measured concentrations can be studied as a Function of overall flow rate or relative flow rates. Initial concentration ratios can be varied. The most time consuming part of the calculation is likely to be the solution of the species conservation equations. This method was not examined for systems having complex chemistry. In particular, systems having widely different rates, i.e. stiffness, may not be easily handled by the method. It would appear to be most suitable for systems containing four or five species interacting via reactions which do not proceed at widely different rates. The particular reactor considered here should be especially suitable for studying systems having secondary reactions and for using mass spectrometric detection. The flow field calculation should present few problems. As already pointed out, the same field can be used for different chemical systems. In fact, if the relative flow rates of the carrier in tubes 1 and 2 are kept the same and only the overall flow changed, it should be possible to simply scale the field resulting
contours for species B. from a single calculation. For example, the arrangement of the stream lines for Re(A) and Re(B) equal to 15 was very similar to that shown in Fig. 3 for a value of 1. Scaling and interpolation among a stored set of selected profiles could make the flow calculation a one time affair for a particular reactor.
REFERENCES
Bird, R. B., Stewart, W. E. & Lightfoot, E. N. (1960) in Transport Phenomena, New York, Wiley. Brown, R. L. (1978). J. &s. Not. Bur. Stand. (US) 83, 1. Howard, C. J. (1979) J. Phy. Chem. 83, 3. Kaufman, F. ( 1961), Prog. React. Kiner. I, 3. Lester, W. G. S. (19611, Rep. Memo. Aeronaur. Rex Com-
mun. 3240. Mills, R. D. (1968), J. Me&. Ennng Sci. 10, 133. Ogren, P. J. (19751, J. Phys. Chem. 79, 1749. Porier, R. V., Carr, R., Jr. (1971), J. Phys. Chem. 75, 1593. Thorn, A. (1933), Proc. Roy. SOC. 141A, 651. Thorn, A. (1934), Rep. Memo. Aeronaur. Res. Commun. 1604.
OF AXIALLY
SYMMETRIC
FLOW
REACTORS
R. L. BROWN Laboratory, Center for Chemical Physics, Chemical Kinetics Division, National Bureau of Standards, Washington, DC 20234, U.S.A. (Receioed
15
April 1983)
Abstract-A method of calculating the velocity profiles and species concentrations in axially symmetric flow reactors is presented. The method is illustrated with a reactor consisting of three concentric tubes arranged so that reactants flow through the inner tube and inner annulus to a mixing region and the products flow out through the outer annulus. A single bimolecular reaction is used in the example. The method involves two steps; calculation of the velocity field from a solution of the Navier-Stokes equations, followed by solution of the species conservation equations. The technique provides a way of analyzing reactors with variable cross sections and sampling ports near mixing regions. The method used for solving the equations of motion is one originally proposed by Thorn (1933). It has been used subsequently by other workers (Lester, 1961; Mills, 1968). From the Navier-Stokes equations in terms of velocities one generates two second order coupled partial differential equations, one for the vorticity 4, and the other for the Stokes stream function I,$. These equations can be solved by standard finite difference methods. By then differentiating the stream function one can get the velocity field for subsequent use in the species continuity equations. For the analysis, incompressible flow was assumed.
INTRODUCHON
The flow tube reactor has been one of the more powerful tools for generating gas phase chemical kinetic data (Howard, 1979). It is normally used under laminar flow conditions at pressures sufficiently low to allow diffusional dispersion to minimize radial concentration gradients. At higher pressures, or in the presence of wall reactions, corrections for radial gradients must be applied to the measurements (Brown, 1978; Kaufman, 1961; Ogren, 1975; Porier & Carr, 1971). In addition, this reactor is normally operated under conditions of a fully developed parabolic velocity profile. Thus, measurements in the region of inlet tubes or in reactors with variable cross sections are avoided. It is, however, relatively easy to mathematically model more complex reactors. Furthermore, the increasing use of low cost high,speed microcomputers by many laboratories makes such modeling an attractive way to extend the range of applicability of flow reactors. The purpose of this work is to show how one can model laminar flow reactors having axial symmetry. As an example of the method a reactor consisting of three concentric tubes has been considered. The inner tube and inner annulus contain reactants in flowing carrier gases which mix and impinge on a wall oriented at right angles to their flow. This reverses their direction and forces them to flow out through the outer annulus. This particular arrangement was chosen to simulate the conditions near a mass spectrometer sampling orifice located very close to the point of mixing of the reactants. The problem is handled by splitting it into two parts. To do this, it is necessary to assume that the reactants are sufficiently buffered with inert carrier gas so that any heat released or absorbed by the reaction does not alter the temperature of the mixture enough to significantly affect the flow fields in the reactor. Then it is possible to first solve the Navier-Stokes equations of motion to get the flow pattern. Following this, the species continuity equations involving diffusion and chemical reaction can be solved using the previously calculated velocity profiles. If the gas flows are kept fixed, these same profiles can be used repeatedly to calculate concentration profiles for entirely different chemical systems. Thus, the equations of motion can be solved for a limited set of flow conditions and the velocity profiles stored for future use.
CALCULATION
OF FLOW FIELD
Derivabn of the difference equations In a cylindrical coordinate system having axial symmetry, the equations of motion for an incompressible fluid with constant viscosity are (Bird et a[., 1960)
au,
u, -
a~,
a.
+ v, -
8,
= -;$+$[;-&!$)+g]
(lb)
where V, and U, are the components of velocity in the radial and axial directions, respectively; p is the density, p is the pressure, and ,Uis the viscosity. There is also an equation for the conservation of mass which must be satisfied; this is
To set up a dimensionless introduce the variables u = V,/%>
w = V,l%,
R=r/a,
z=
z/u,
set of equations,
we
p = P l(P%% Re = av,,p 1~
(3)
where u0 is a reference velocity and a is a reference distance. These will be chosen later. 139
140
R. L. BROWN
Substituting these dimensionless variables into Eqs (la), (lb) and (2), we get the following set of equations,
Equations (8) and (9) can be solved by finite difference methods. To convert these equations into finite difference form, three point centered difference approximations ,for the derivatives were used with a square grid having a dimensionless spacing H. If G represents either q or I& then the formulas used were aG -z
G,,, +
ilR
aw
waR+waz=-z
a~
aw
a2G @= aG Fz” d2G @=
(5) We next introduce a dimensionless stream function $J by the relations
aw
au
“=E-x
~=‘f!i R
az'
vorticity
W=
9 and
(6)
_ltP.
RdR
I - G,,,- I 2H
’
Gi,j+~ - 2G, + G,>j_, ff= ’ (10)
G,+cij--t-u 2H
’
G,+,j-2Gij+Gi_,,j
HZ
.
These are approximations to the derivatives at the grid point (i, j), where the index i refers to points along the Z axis, and j refers to points along the R axis. Substituting these formulas into Eqs. (8) and (9) and solving for GLiJand v,~ gives the formulas
(7)
Substitution of Eq. (7) into Eq. (5) shows that the continuity equation is automatically satisfied. To derive an equation for the stream function, we take the derivative of U with respect to Z and that of W with respect to R in Eq. (7). Adding these and inserting Eq. (6) we get
(8) -
An equation for the vorticity can be gotten by taking the derivative of Eq. (4a) with respect to 2, subtracting the derivative of Eq. (4b) with respect to R, and inserting the vorticity wherever the expression of Eq. (6) appears. The resulting equation is
Mi,,
+ I -@t+i.j+I)(tli+l.j-
Except for a difference indices, these equations by Mills (1968).
vi-
,Jl/W4)}
in the identification of the are the same as those given
Description of reactor model The geometry of the sample reactor is shown in Fig. 1. This reactor consists of three concentric tubes with walls of finite thickness. Carrier gas containing reactant B flows through tube 1 toward the vertical wall at the origin. This tube ends at a distance z, to the right of this wall. We call the part of tube 1 to the
wa -
Fig. 1. Geometry
of the axially
(13
symmetric flow reactor used in the sample
calculation.
Modeling
of axially symmetric
right of z,, region B. Reactant A, contained in the same type of carrier gas, flows in the same direction between the outer wall of tube 1 and the inner wall of tube 2. This second tube ends at a distance z2 from the origin. The annular region between tubes 1 and 2 to the right of Z, we denote as region A. These reactants mix by diffusion in the region to the left of Z, and flow away from the origin through the annular space between the outer wall of tube 2 and the inner wall of tube 3. The portion of this space to the right of z2 will be called region C. The flows in these three regions will have parabolic velocity profiles provided one is sufficiently far away from the origin. There, the flow will have only an axial component. These limiting profiles for the three regions are (Bird er al., 1960),
141
Row reactors
sionless relation (kX5 + &I
Re@
between
the Reynolds
1 + k,, ReP)
= (1 + k&Re(C)
(20)
Let us now write dimensionless versions of Eqs. (13HlS) for the limiting velocity profiles. The reference length and velocity indicated in Eqs. (3) are taken to be a5 and {wc), respectively. We get )
x[l -(z>‘-_fAln(~)], k25SR5k35(22) WB=
WB/{ WC}= - 2
z
(ks, - k,,)
+-($1, (14)
numbers:
W, = wc,/(wc}
= 2
(23)
OSRSk,5 [I - RZ -f&
In(R)],
kd55 R 5 1. (24)
In these equations, (w”). (ws>, and (wc> are the average axial velocities, and the a, are the tube radii; the other quantities are given by k,i = q/a,
(161
1 +!c;+“&
g,=
next define
ReV
Reynolds
I=
numbers
- (% - a&
ReVI
= - alp (w,>/P
Re(C)
= G+ - a&
(17b) for
these re-
{w,>/P (18)
(wc>i~.
-w(u3+aJRe(A)
mB = - apa, Re (B) mc = np(a,
(191
+ KJ Re (C)
For conservation of mass, these three equations must add to zero. (Note that {w”} and (wa> are negative in this system.) Summing Eqs. (19) gives a dimen-
s
RWdR+O
(25)
where the integration is over the range of R values covered by a particular region, and 8 is the constant of integration. Integrating Eqs. (22)-(24) gives the limiting stream functions. lird = - $AR’[l
Consider the total mass flow rate of gas m, through any one of these limiting regions. This is given by the integral m=2~Spwrdr which is equal to 2np f wr dr since p is assumed to be constant. Because the average velocity here is 27~j wr dr /A, where A is the cross sectional area of the region, the mass flow rate can be written in terms of the average velocity as m = p (w )A. The mass flow rates can then be expressed in terms of the Reynolds numbers through use of Eq. (18). Thus, for the three limiting regions, the mass flow rates are m,=
$=-
(17a)
& = ( 1 - k$)/ln (kg) We gions;
Boundary conditions on 9 and q From the above limiting velocities far upstream and far downstream, we can derive expressions for the limiting values of the stream function. These can be calculated by integrating the expression for W in Eq. (7). We get
- #lU2
l+k,= -:BRV Gc = -
-&On
(R/k,,)
-+(R/k,$]
R*[l - iR* -j&In
R -i)J/g,,
-
+ ea + Br
$1 + 8, 126) (27) (28)
where A
=
_
2Re(AI(1 - k44 g2dWC)h - kz5)
and B = -
2Re(B)(k,,
- k,,)/Re(C)
To obtain expressions for the integration constants, we make the following observations. At any wall in the reactor the normal velocity component must vanish. Likewise the tangential component must also vanish if there is to be a no-slip condition at the walls. From Eq. (7), this means that the gradient of + must be zero at the walls. Thus, 1(1is constant along a wall. In addition, from symmetry, the radial velocity U must vanish at R = 0. Thus, from Eq. (71, + must be a constant all along the reactor axis. Since the axis hits the vertical wall at the end of tube 3, IL
R. L. BROWN
142
at the center must be the same as at the inner wall of tube 3 if we are to have continuity in # at the origin. If we arbitrarily set $ equal to zero along the axis, we get from Eq. (27) that 8, = 0. Then we can immediately get 0, from the condition I,!J,=(%= 1) = 0. The constant 0, can be obtained from the condition GA(R = k,,) = @_XR = ke5) which must be satisfied if II, is to be constant along the inner and outer surfaces of tube 2. The results are Be = (1 +J43/(2g‘S)
8, = ~~-Ak:,(l+hd
-
k:,(k:, +f45M2gM)+ 0,.
+fa/R)
‘I~ = 2BR /k f5 ‘~c =
2(2R + f.slR h’g=w
%j~t., = 3(+_, + I J -
Jl,,,Y~* - fR,vr+ IJ
(36)
For the vertical walls at the ends of tubes 1 and 2, we have %,rli., = 3(+,-i,,
- i&)/H2
- +R,?, - I.,
(37)
(29)
Having obtained values for all the mesh points on the boundaries, we can now show how Eqs. (11) and (12) were solved.
(30)
Solution
These expressions for the limiting stream functions allow us to calculate all the boundary conditions on JI. Thus, far upstream in regions A and B, Eqs. (26) and (27) aive 11 as a function of R. Far downstream, this is &en by Eq. (28). The constancy of I& at the walls and reactor axis give all the other required boundary values. The limiting expressions for the vorticity are obtained from the velocities WA, W, and WC by differentiation. We get rtA = A(2%lk:,
gives
(31) (32) (33)
These can be used to calculate the boundaries for g, far upstream and far downstream. Because U and dW/dR are zero along the axis, we have from Eq. (6) that ?I is also zero along this boundary. At the remaining boundaries, rf is not known but must be caluclated from the flow field in the region near the boundary. Lester (1961) has derived a finite difference formula for q on a boundary wall inclined at an angle cp to the z axis. If rl, is the value at the wall and n, the value at a point one mesh interval H into the flow field in a direction normal to the wall, then
x [2 + H cos cp/R, - +H’ cos2 cp/Rw2] - ‘. (34) The angle cp is measured in a counter clockwise direction and when cp = 0, the normal is taken to be in the positive R direction. Thus, for the horizontal walls corresponding to the inner surfaces of the three tubes we have rp = 180” since their normals into the fluid point in the negative R directions. For q on the outer walls of tubes 1 and 2, we have cp = 0. Thus, if qij is one of these horizontal walls, we have from Eq. (34)
where the upper signs correspond to the outer walls. For the vertical wall at the end of the reactor we have ‘p = - 90”. Therefore, for qi,j on this wall, Eq. (34)
of finite
dzxerence
equations
for
flow
field
The following values were used for the reactor parameters; u, = 0.5 cm, a2 = 0.6 cm, uj = 1.2 cm, a, = 1.3 cm and u5 = 1.7 cm. The wall thickness of the tubes was chosen to be equal to one mesh width. This produced 18 mesh points along the r axis. For cases where the wall thickness is not a multiple of the mesh width, modified finite difference formulas for the boundaries would have to be derived. We did not consider this problem. The dimensionless mesh width was H = h/a,, = 1/17. A small section of the mesh is shown in Fig. 1. The distances z, and rS of the ends of tubes 1 and 2 from the origin were taken to be 16h and 7h, respectively. The Reynolds numbers for the flow in regions A and % were both set equal to 1. Because of the relation between the Reynolds numbers, Eq. (20), this choice gives Re(C) = 0.759. For these conditions, the upstream and downstream limiting values for the axial velocities, stream functions, and vorticities are shown in Fig. 2 as a function of R. At the start of the iteration process, the limiting profiles were assigned to both @ and 9 at the upstream and downstream boundaries. The distances at which these boundaries are placed can be varied to see if their positions have any effect on the results in the mixing region. For the present conditions, distances of 46h from the origin were found to be satisfactory. At higher flows, the downstream distance for the limiting boundary would probably have to be increased. For the stream function, the limiting profile values at R = 0 and at the tube walls were continued around all of the walls. The limiting vorticity values at the walls were continued only to the ends of the tubes leaving a discontinuity there; as can be seen from Fig. 2, in the final solution, the vorticity along a wall must change continuously from the limiting value on one side to that on the other. For the vertical wall at the origin, the vorticity was initially given the average of its limiting values at R = 0 and 1. For the first run, all the internal mesh points for IJ and tf were given the value zero. In subsequent runs at different Re values, previously calculated results were used for the initial conditions. The iterative routine described by Mills (1968) was used. For this, one performs a complete iteration of the II, Eeld with Eq. (11) followed by a complete iteration of the q field with Eq. (12). In between these iterations the vorticity boundary values are calculated from Eqs. (35), (36) or (37). The gelds were scanned in the following way: starting at the upstream boundaries in regions A and B, at a fixed value of Z, R was scanned from R = 0 to the inner wall of tube 2; Z was then decreased to the next downstream mesh value
Modeling of axially symmetric fiow reactors
Fig. 2. The limiting profiles for the velocity, stream function, and vorticity in the upstream B and the downstream region C: Re(A) = Re(B) = 1, Re(C) = 0.759. and this R scan was repeated. This process was continued until Z, was reached. For the region bounded by q, a.,, and the end wall of the reactor, the scan was 6rst made for Z = 0 to Z, with R fixed at U. This 2 scan was repeated with R increased one mesh point until R reached the position of the outer wall of tube 2. Then the scan was again made over R, from the outer wall of tube 2 to the inner wall of tube 3, starting at Z = N and increasing 2 until the Iimiting downstream position in region C was reached. The comers at the ends of tubes 1 and 2 were treated according to the way suggested by Thorn (1934) as discussed by Lester (1961) and Mills (1968).
143
r+ons
_,gand
The stabiiity and convergence of the method can be estimated (Bird et ok, 1960) from the condition W s 4@/Re. Here, the Re refers to the largest Reynolds number in the system. There are reasons for limiting the method to cases where Re 5 50, and Milk (1968) gives a discussion of these. At pressures below 10 or 20 torr, this range of Reynolds numbers covers the flow conditions most often used in these reactors. For example, in a nitrogen carrier at 5 torr, the average flow velocity in tube 1 is 47 Re(A) cm/s. The results for the sample reactor are shown in Fig. 3. This is a contour diagram for the stream function JI. The direction of the flow is indicated on each
... ... ... ... -.. ... ... .-- ..* ... ... -.. ... ._. .._... .. .. . ... ... ._. .._ ... ... ... _.. ... .-. .._ ... ... ... _.. _..... .. Fig. 3. Stream line contours:
&(A
) = Re(B) = 1.
0
R. L.
144
stream line. As mentioned earlier the velocity field is obtained from the gradient of $. Normally, one can say that the closer spaced are the stream function contours the higher the fluid velocity. This is deceptive, however, for cylindrical systems because the gradient of $ must be weighted by l/R to get the velocity field. Thus, the average upstream velocity in the annulus between tubes 1 and 2 appears from the stream line density to be much greater than that in tube 1, whereas, for this case, it happends to be equal. Having shown how to calculate the velocity field, we can now consider how to calculate the species concentrations in a reacting system. DETERMINATION M
OF
A CHEMICALLY
CONCENTRATION REACTING
PROFILES SYSTEM
Let c” be the concentration of some species in the reactor. Its conservation equation is (Bird et al., 1960)
where dA is the diffusion coefficient of the species in the carrier gas and So is its production rate through chemical reactions. There will be one of these equations for each different reacting species. Introduction of the dimensionless variables of EQ. {3) gives
where D, = dJ(av,), S, = asA/(cOvO) and C, = cl/C’,. The quantity c, is some reference concentration. Since U and W are to be obtained from the previous calculation, we must here use the same reference length and velocity, namely a5 and (wc). We assume that S, is proportional to C, and write S, = PACA, where P, is some function of the species concentrations possibly including that of A. Using the derivative approximations given by Eq. (10) in Eq. (39) we get the difference formula
BROWN
gradient approximation can be used. A five point difference formula was used to calculate the concentration at a boundary. This was obtained by setting a five point formula for the derivative equal to zero. If C, is the concentration on the boundary and C,+i are the values at successive mesh points along a line perpendicular to the boundary and into the gas, then for a zero derivative at the boundary, we used the formula C, = 12(4Cb+,
- 3Cb+z + 4C,+,/3 - c&,/4)/25.
A three point formllla probably would have been satisfactory for these calculations but was not tested. Instead of the zero gradient downstream approximation another approach might be to set the axial gradient equal to S,/ W which is a reasonable approximation when S, is small. Then one could consider reactions which’had not gone as far toward completion at the downstream boundary. For the iterations, the field scanning was done in the same manner as for the $ and r] fields. To show the kind of results that can be obtained, a single bimolecular reaction was examined. One reactant entered from region A, and the other from region B. The rate function PA in this case is c,akCJu,, where k is the bimolecular rate constant for the reaction. The rate function P, is cgzkC,fv,. The flow field shown in Fig. 3 was used. Initial dimensionless concentrations were taken to be C, = 0.5 and C, = 1.7. These values were chosen so that the molar flow rates into the reactor of the two species would be the same. Dimensionless diffusion coefficients were given 7 the values D, = 0.093 and D, = 0.052. The value of the dimensionless rate constant c,ak/v, was taken to be 1.3. The results for the sohttion of Eq. (38) are shown in Figs. 4 and 5 as contour diagrams for the concentrations C,, and C,. These indicate clearly the combined effects of reaction and diffusion. If the reaction rate were zero, the limiting downstream concentrations in this case would be 0.386. Both species spread out considerably in the region to the left of the end of tube I. At these flows, there is, however, very little back diffusion into the adjacent tube by either species. The effect of the parabolic velocity profile is evident in the downstream regions in the form of higher concentrations in the middle of the flow stream than near the walls. CONCLUSIONS
For simplicity the species subscript has been dropped in this equation. For boundary conditions, we have a zero concentration gradient at the walls and perpendicular to the Z axis at R = 0. Far upstream in regions A and B the concentrations are fixed at their initial values, while far downstream we assume that the reactions have proceeded so far toward completion that a zerotThese values for the dimensionless coefficients are applicable to the diKusion in a reactor of the present dimensions of N atoms and NO in N2 at pressures around 5 torr.
The aim of this work has been to show in some detail one method of modeling an axially symmetric flow reactor containing a relatively simple chemical system. The specific example chosen for modeling has shown the type of results obtainable. These consist of concentration distributions for the individual chemical species as a function of carrier flow rates, species diffusion coefficients, and chemical rate constants. Since these calculations give the concentrations throughout the reactor, one can relate the concentration measured at any particular point to the parameters of the system. For the reactor studied here, the most obvious sampling point would be at the origin. A variety of experiments can be envisioned with this reactor. For example, one could measure diffusion coefficients by putting only one species at a time through the system, or by using non-reactive
Modeling
of axially symmetric flow reactors
Fig. 4. Concentration
contours for species
145
A.
-------
Fig. 5. Concentration species. The measured concentrations can be studied as a Function of overall flow rate or relative flow rates. Initial concentration ratios can be varied. The most time consuming part of the calculation is likely to be the solution of the species conservation equations. This method was not examined for systems having complex chemistry. In particular, systems having widely different rates, i.e. stiffness, may not be easily handled by the method. It would appear to be most suitable for systems containing four or five species interacting via reactions which do not proceed at widely different rates. The particular reactor considered here should be especially suitable for studying systems having secondary reactions and for using mass spectrometric detection. The flow field calculation should present few problems. As already pointed out, the same field can be used for different chemical systems. In fact, if the relative flow rates of the carrier in tubes 1 and 2 are kept the same and only the overall flow changed, it should be possible to simply scale the field resulting
contours for species B. from a single calculation. For example, the arrangement of the stream lines for Re(A) and Re(B) equal to 15 was very similar to that shown in Fig. 3 for a value of 1. Scaling and interpolation among a stored set of selected profiles could make the flow calculation a one time affair for a particular reactor.
REFERENCES
Bird, R. B., Stewart, W. E. & Lightfoot, E. N. (1960) in Transport Phenomena, New York, Wiley. Brown, R. L. (1978). J. &s. Not. Bur. Stand. (US) 83, 1. Howard, C. J. (1979) J. Phy. Chem. 83, 3. Kaufman, F. ( 1961), Prog. React. Kiner. I, 3. Lester, W. G. S. (19611, Rep. Memo. Aeronaur. Rex Com-
mun. 3240. Mills, R. D. (1968), J. Me&. Ennng Sci. 10, 133. Ogren, P. J. (19751, J. Phys. Chem. 79, 1749. Porier, R. V., Carr, R., Jr. (1971), J. Phys. Chem. 75, 1593. Thorn, A. (1933), Proc. Roy. SOC. 141A, 651. Thorn, A. (1934), Rep. Memo. Aeronaur. Res. Commun. 1604.