Ind. Eng. Chem. Res. 2005, 44, 5031-5035
5031
On Axial Diffusion in Laminar-Flow Reactors E. Bruce Nauman* and Ashish Nigam The Isermann Department of Chemical and Biological Engineering, Rensselaer Polytechnic Institute, Troy, New York 12180
Although unimportant in conventional laminar-flow reactors, axial diffusion can be the dominant mechanism for mixing in micrometer-scale devices. The low Reynolds numbers characteristic of these devices allow numerical solution of the convective diffusion equation, coupled when necessary with the equations of motion. The method of false transients provides a simple way to obtain such numerical solutions. For the case of reaction within a confined zone, there is an internal optimum with respect to diffusion. Unlike the lumped-parameter axial dispersion model, the conversion and residence time distribution depend on the boundary conditions. The physically realizable open-open boundary conditions give higher conversions and more uniform residence times. Micrometer-scale devices will have such high levels of diffusion that mixing between separately fed, compatible components will become trivially easy. Reactor performance will approximate that of a CSTR with attendant low yields and poor selectivities. Introduction This paper is focused on numerical solutions to the convective diffusion equation and on the circumstances in which various terms in that equation can become important. Specifically considered are flows between flat plates, in ducts, or in tubes. Suppose there is a single reactive species that obeys Fick’s law of diffusion and has a constant molecular diffusivity, D. Then, the governing equation, written in rectangular coordinates, is
(
)
∂c ∂c ∂c ∂c ∂2c ∂2c ∂2c + ux + uy + uz ) D 2 + 2 + 2 + Rc ∂t ∂x ∂y ∂z ∂x ∂y ∂z (1) Suppose there is some predominant flow direction, z, which we will term the axial direction, and that uX and uy are the relatively minor cross-channel components that are induced by changes in viscosity as the reaction proceeds or by the presence of mixing elements (i.e., motionless mixers, static mixers). Cleland and Wilhelm1 gave numerical solutions for a first-order reaction in a circular tube with a parabolic velocity profile and with radial diffusion but no axial diffusion. Specifically, they solved
uz
(
)
∂c ∂c 1 ∂c ∂2c + Rc (2) ) 2u j (1 - r2/R2) ) D + ∂z ∂z r ∂r ∂r2
When the reaction is first-order, Rc ) -kc, eq 2 can be transformed into the Graetz equation for which an analytical solution is known. For practical purposes, however, numerical solutions are preferred. Equation 2 has usually been solved using the method of lines with the lines running in the axial direction and with fully explicit finite differences in the radial direction. The key to this solution technique is the absence of the second derivative in the z direction that arises when axial diffusion is considered. When there is only radial diffusion, eq 2 can be approximated as a set of first* To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: 518-276-6726. Fax: 518-382-2912.
Figure 1. First-order reaction with kth ) 1 in a laminar-flow reactor with a parabolic velocity profile, as calculated from eq 2.
order ordinary differential equations in the independent variable z that have known initial values, c(r,0) ) cin(r). Inclusion of the second derivative converts the problem to a numerically awkward set of two-point boundary value problems with Danckwerts-type boundary conditions. Figure 1 illustrates the solution to eq 2 for the case of laminar flow between flat plates. The distance between the plates is 2Y. Results are similar for flow in rectangular ducts and circular tubes.2 As shown, the performance of the laminar-flow reactor asymptotically approaches that of a piston-flow reactor as the diffusivity increases. However, upon reflection, it is apparent that increasing the diffusion group, Dth/Y2, to arbitrarily high values will drive the reactor toward complete mixing [i.e., toward CSTR (continuous stirred tank reactor) behavior] if the aspect ratio, L/Y, is fixed. The results in Figure 1 apply only in the limit of high aspect ratio. Essentially all studies to date have ignored axial diffusion. There are two reasons for this: Axial diffusion is, in fact, unimportant in conventional laminar-flow reactors such as those used for commercial polymerizations, and inclusion of axial diffusion prevents numerical solution by the established mythology. The diffusion group, Dth/R2 or Dth/Y2, is typically less than 0.003, so that even radial diffusion can be ignored.3 Microfluidic devices fabricated by soft lithography have a size scale of Y ≈ 100 µm, so that Dth/Y2 ≈ 0.1 for fast reactions in water-like media. As indicated in Figure 1, the effects
10.1021/ie049677b CCC: $30.25 © 2005 American Chemical Society Published on Web 08/24/2004
5032
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Figure 2. Combined effects of radial and axial diffusion on the conversion of a first-order reaction with kth ) 1.
of even cross-channel diffusion remain rather small. However, values of Dth/Y2 of 100 or more can be expected as device sizes shrink toward 1 µm. Such large values of the diffusion group lead to very different performance than predicted from eq 2. Specifically, the cross-channel gradients in concentration will have been eliminated, and even the axial gradient will have lessened. Equation 1 can be solved using the method of lines in the time direction and fully explicit finite differences in the spatial directions. Integrating to long times gives the steady-state solution. The approach is called the method of false transients, although the transients in the present case could be achieved physically. Details on the solution technique are postponed until after the results are discussed. First-Order Reactions Figure 2 shows results, computed by the method of false transients, for the flat-plate version of eq 2, but with inclusion of the axial diffusion term. The equation being solved is
1.5u j (1 - y2/Y2)
(
)
∂2c ∂2c ∂c ) D 2 + 2 - kc ∂z ∂y ∂z
(3)
The boundary conditions are those for an open system with indefinitely long entrance and exit regions in which there is convection and diffusion but no reaction
c(x,y,-∞) ) cin
(4)
dc(x,y,∞) )0 dz
(5)
The striking thing about Figure 2 is the internal optimum. For a fixed aspect ratio, L/Y, there is a value for the diffusion group, Dth/Y2, that gives the best performance. Values for Dth/Y2 that are much lower than the optimum give the performance of a laminar-flow reactor without diffusion (LFR). Values much higher than the optimum give the performance of a CSTR, although this fact is not clearly shown in the figure because of numerical ill-conditioning of the solution technique at very high values of the diffusion parameter. Intermediate values can give conversions superior to either the LFR or CSTR limits, with values near piston flow for long reactors. Figure 3 shows the location of the interior optimum. From a design perspective, the value for Dth/Y2 and kth will be set by the size scale of the fabrication technique, by the desired extent of reaction, and by the physical
Figure 3. Location of optimum yield for first-order reactions in a laminar-flow flat-plate reactor with diffusion in the cross-channel and axial directions.
Figure 4. Effect of combined cross-channel and axial diffusion on a laminar-flow reactor.
properties of the reaction mixture. The optimal aspect ratio can then be determined from Figure 3. The apparent linearity of the results suggests that extrapolations to large values of Dth/Y2 are possible. It is clear that, with Dth/Y2 ≈ 100, reactors with very high aspect ratios, L/Y ≈ 1000, will be needed achieve performance within the range conventionally expected from a laminar-flow reactor, i.e., reaction yields bounded by those for an LFR and a piston-flow reactor. Such long reactors might not be feasible due to pressure drop. Although the Reynolds number is less than 1, the pressure drop will be approximately 108 Pa/m at a flow rate of 10-2 m/s for a channel with Y ≈ 1 µm. An acceptable pressure drop4 of about 10 kPa restricts L/Y to about 100 so that axial diffusion will be important. Short reactors with L/Y < 10, might be suitable for radiation-induced reactions,5 will approximate the behavior of a CSTR, with attendant lower yields and selectivities. (We ignore autocatalytic reactions for which a CSTR can perform better than a piston-flow reactor.) Extrapolation of the results in Figure 3 to low values of Dth/Y2 indicates that essentially any level of molecular diffusion decreases the performance of very short reactors, so that an interior optimum with respect to Dth/Y2 no longer exists. This is intuitively reasonable given that the beneficial effects of cross-channel diffusion are outweighed by the detrimental effects of axial diffusion when L/Y is small. Figure 4 shows the effect of the aspect ratio on the yield of a first-order reaction with kth ) 1 and Dth/Y2 ) 1. This represents about a factor of 3 reduction from the current scale of soft lithography. It is apparent that axial diffusion will shortly become important in the design of microfluidic devices and that the current difficulty in achieving adequate cross-channel mixing (e.g., between initially unmixed reactants) will cease to be a problem.
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5033
Figure 5. Residence time density functions for values of Dth/Y2.
Residence Time Distributions In closed systems, the residence time distribution can be measured using inert tracers. For example, the output response of the system to a unit impulse of tracer (i.e., a Dirac input) gives the residence time density or differential distribution function, f(t). Such measurements are theoretically impossible in an open system when there is significant axial diffusion because time spent outside the system boundaries is included in the inert tracer results. As one consequence, the apparent value for ht as determined from inert tracer experiments will be higher than the actual value calculated using ht ) L/u j . A reasonable approximation of the extent of overestimation can be obtained from
ht apparent Dth Y 2 )1+2 2 ht actual Y L
()
(6)
This result was obtained for the axial dispersion model.6 Its use here for reactors with high values of Dth/Y2 is justified on the grounds that diffusion will have substantially eliminated cross-channel concentration gradients. For the case of Dth/Y2 ) 3 and L/Y ) 3, inert tracer experiments would overestimate ht by 67%. If the residence time distribution is desired, for example, to calculate bounds on the conversion of reactions other than first order, it can be calculated from solutions of eq 3 for various values of the rate constant, k. The residence time density function is then determined by inverse Laplace transformation
cout ) cin
∫0∞e-ktf(t) dt
(7)
Figure 5 shows results for a reactor with L/Y ) 3. For this case, Dth/Y2 ) 3 is near the optimum and gives the highest conversion. Dth/Y2 ) 9 approximates the residence time distribution of a CSTR, whereas Dth/Y2 ) 1 begins to show the slowly convergent tail that corresponds to laminar flow without diffusion. Open versus Closed Systems The assumption of an open system deserves comment. In the lumped-parameter axial dispersion model, open and closed systems have the same residence time distribution within the reaction zone and will have the same reactor performance.6 This is no longer true for the distributed-parameter models of eqs 1-3. Closed boundary conditionssthese being generalizations of the
Figure 6. Comparison of conversions for open-closed versus open-open conditions.
Danckwerts’ boundary conditions that conserve molecular fluxsare
uz(x,y)cin ) uz(x,y) c(x,y,0) - D dc(x,y,L) )0 dz
(∂z∂c)
z)0
(8) (9)
In the axial dispersion model, the concentrations at location z ) 0 and z ) L turn out to be identical for open and closed systems. In the distributed-parameter case, concentrations vary in the cross-channel directions, and cross-channel diffusion will exist outside the reaction zone in an open system. For example, c(x,y,L) varies from point to point on the exit plane and has different values under eq 9 than it would have if downstream diffusion were allowed to occur in the indefinitely long exit zone of eq 5. Figure 6 illustrates the difference between open and closed boundary conditions at the reactor exit. The inlet boundary conditions used in Figure 6 were open (eq 4). The difference is appreciable for high values of Dth/Y2. Diffusion in the exit zone increases conversion by promoting a more uniform concentration profile and thus a more uniform residence time distribution within the reaction zone. Second-Order Reactions We have so far considered only the case where Rc ) 0 outside the prescribed reaction zone, 0 < z < L. This situation is physically realizable for radiation-induced reactions5 and might be a reasonable approximation when heat can be applied to a limited portion of a flow channel or when the channel walls are catalytic within the reaction zone. However, the majority of reactions anticipated for microfluidic devices require mixing of initially unmixed components. Here, open boundary conditions are suitable to model the separate entrance channels for the components of a second-order reaction. The geometry is shown in Figure 7. Versions of eq 1 are written for the individual components. Reaction is possible everywhere, including the separate entrance channels, and will begin well upstream of the point
5034
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Numerical Solution Techniques
Figure 7. Geometry of an open system suitable for a second-order reaction.
Figure 8. Mixing-cup average concentration profiles for a secondorder reaction in the flat-plate geometry of Figure 6.
Figure 9. Outlet concentration from the two-channel reactor as a function of the diffusion group.
where the two channels join if Dth/Y2 is large.5 This is illustrated in Figure 8 for an elementary second-order reaction of the form A + B f product with dimensionless rate constant ainkth ) 1, where a ) [A]. The results apply to the symmetrical case of equal concentrations, channel heights, and flow rates. Note that the reaction continues indefinitely downstream. Thus, the magnitude of ht is somewhat arbitrary. It corresponds to the value ht ) L/u j , where L is the, again somewhat arbitrary, reactor outlet. Figure 9 shows the concentration at a fixed downstream location, the “outlet”, as a function of Dth/Y2. Three regions can be identified: For low values of Dth/Y2, the reaction is limited by mixing of the separately fed components; this is radial or cross-channel mixing. At intermediate values of Dth/Y2, mixing in the crosschannel direction is essentially complete, and the reaction yield is limited by the value of ht. At high values of Dth/Y2, the effective volume of the reaction zone upstream of the outlet is expanded by axial diffusion. The phenomenon is similar to that discussed in connection with eq 6. With an open boundary condition in the downstream direction, the potential volume of the reactor is unbounded, and the reaction goes to completion at high values of Dth/Y2.
Second-order finite-difference approximations are used for all spatial derivatives in eq 1. This gives a set ODEs, one for each spatial point, that are coupled through dependence on concentrations at nearestneighbor points. Zero-flux boundary conditions with second-order convergence are applied at all solid boundaries (or at lines of symmetry). Inlet and outlet boundary conditions of the open type are numerically approximated by using inlet and outlet zones that are arbitrarily long and are governed by eq 1 but with Rc ) 0. All crosschannel concentrations at the first and second points of the inlet zone are set equal to cin, and the concentrations at the last point of the outlet zone are set equal to those at the previous point. This forces zero slope at the beginning of the inlet zone and at the end of the outlet zone. The adequacy of the lengths assumed for the inlet and outlets regions is checked by testing that the derivatives are indeed small at the two ends and that the computed average outlet concentration is insensitive to changes in the lengths. The initial condition is simply that c ) cin everywhere. Euler’s method is preferred to integrate the set of time-dependent ODEs because only the long-time solution is sought. The solution technique must be stable, but computational accuracy is unnecessary at small times. An outlet boundary condition of the closed type, eq 9, is implemented by forcing zero slope at z ) L. It is much more difficult to implement a closed boundary condition (eq 8) at the reactor inlet because the numerical solution requires guessing a function, c(x,y,0). However, when Dth/Y2 is small, open and closed systems will have c ≈ cin at every point on the inlet plane to the reactor. The zero-slope condition at the outlet will also be satisfied. Thus, the open and closed systems give identical results for diffusion-free laminar flow. Identical results are also obtained in the high conversion limit, but appreciable differences exist at intermediate values of Dth/Y2. As a practical consideration, a closed system would be difficult to design in a microfluidic device. As suggested in the Introduction, changes in the axial velocity profile might arise because of changes in viscosity upon reaction or as a result of the installation of mixing elements. When Dth/Y2 is large, cross-channel viscosity differences will vanish, and there is little need for mixing elements. However, changes in the axial velocity profile can easily be incorporated into the numerical scheme used to calculate concentration profiles. For tubular reactors, there is an established way to calculate uz(r) when viscosity is a function of strong function of r but a slowly varying function of z (e.g., see Chapter 8 of Nauman2), and this methodology is easily extended to the flat-plate geometry. For rectangular ducts, an analytical solution for the axial velocity profile is awkward even for isoviscometric flow. However, it is easily calculated by the method of false transients. The equation to be solved is
F
duz dP )+ ∇µ∇uz dt dz
(10)
subject to no-slip boundary conditions at the walls. A suitable initial condition is that u(0,x,y,z) ) u j . The average velocity, u j , is proportional to the axial pressure drop, but the shape of the velocity profile is independent of the pressure drop. When viscosity is variable, the pressure drop will also vary as a function of z. In
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5035
numerical calculations, the pressure drop is adjusted as necessary to maintain a constant mass flow down the channel. Equation 6 assumes that the cross-channel velocity components are negligible. This assumption is reasonable at low Reynolds numbers. However, it has been shown7 that, although negligible for purposes of calculating uz, the cross-channel components can have a significant, cumulative effect on reaction yields. For flow in tubes and between flat plates, there is only one velocity component, and it can be recovered from the continuity equation. For flow in rectangular ducts, there are two cross-channel components, and no simple way has been developed to calculate them. Conclusions Microfluidic devices are becoming important in analytical and combinatorial chemistry, as well as for the production of designer molecules. This paper has explored previously ignored effects of axial diffusion in laminar-flow reactors. It is predicted that these effects will soon become important, when the length-scale of microfluidic devices is reduced by a factor of about 3 from what is currently possible with soft lithography. Indeed, if a microfluidic device were designed for a gasphase reaction, axial diffusion would be important at the current length scale. When microfluidic devices are actually at a micrometer length scale, the effects of axial diffusion will become manifestly important to the point of eliminating current design concepts for liquid-phase laminar-flow reactors. Fortunately, the low Reynolds numbers characteristic of microfluidic devices allow numerical solution of the convective diffusion equation, coupled when necessary with the equations of motion. The use of method of false transients provides a simple and easy way to obtain such numerical solutions. Notation a ) concentration of reactant A ain ) inlet concentration of reactant A aout ) mixing-cup average outlet concentration of reactant A A ) reactant A
bin ) inlet concentration of reactant B B ) reactant B c ) concentration of reactant cin ) inlet concentration of reactant cout ) mixing-cup average outlet concentration D ) diffusivity of reactant f(t) ) differential distribution function for residence times k ) reaction rate constant P ) pressure r ) radial coordinate R ) radius of the reactor RC ) rate of reaction t ) time ht ) mean residence time u ) velocity u j ) average axial velocity x ) x coordinate, cross-channel direction Y ) half-distance between parallel plates y ) y coordinate, coordinate in the cross-channel direction z ) z coordinate F ) mass density µ ) viscosity
Literature Cited (1) Cleland, F. A.; Wilhelm, R. H. Diffusion and reaction in viscous-flow tubular reactor. AIChE J. 1956, 2, 489-497. (2) Nauman, E. B. Chemical Reactor Design Optimization and Scale-up; McGraw-Hill: New York, 2002; pp 284, 197-201. (3) Merill, L. S., Jr.; Hamrin, C. E., Jr. Conversion and Temperature Profiles for Complex Reactions in Laminar and Plug Flow. AIChE J. 1970, 6, 194-198. (4) Papautsky, I.; Brazzle, J.; Ameel, T.; Frazier, A. B. Laminar fluid behavior in microchannels using micropolar fluid theory. Sens. Actuators A: Phys. 1999, 73, 101-108. (5) Nauman, E. B.; Nigam, A. Mixing in Sub-micron Ducts. Chem. Eng. Technol. 2004, 27, 293-296. (6) Nauman, E. B. Invited Review: Residence Time Distributions and Micromixing. Chem. Eng. Commun. 1981, 8, 53-131. (7) McLaughlin, H. S.; Mallikarjun, R.; Nauman, E. B. The Effect of Radial Velocities on Laminar Flow, Tubular Reactor Models. AIChE J. 1986, 32, 419-425.
Received for review April 21, 2004 Revised manuscript received July 9, 2004 Accepted July 14, 2004 IE049677B
