Reaction enhancement of point sources due to vortex stirring

PHYSICAL REVIEW E 74, 016307 共2006兲 Reaction enhancement of point sources due to vortex stirring 1 John P. Crimaldi,1,* Jillian R. Hartford,1 and Jeffrey B. Weiss2 Civil and Environmental Engineering, University of Colorado, Boulder, Colorado 80309-0428, USA Department of Atmospheric and Oceanic Sciences, University of Colorado, Boulder, Colorado 80309-0311, USA 共Received 17 November 2005; published 31 July 2006兲 2 We investigate a class of reactive advection-diffusion problems motivated by an ecological mixing process. We use analytical and numerical methods to determine reaction rates between two initially distinct scalar point masses that are separated from one another by a third 共nonreactive兲 scalar. The scalars are stirred by a single two-dimensional vortex in a variety of geometrical configurations. We show that the aggregate second-order reaction rate in the low-concentration limit is enhanced by the instantaneous stirring processes, relative to the rate predicted by an equivalent eddy diffusivity. The peak reaction rate grows as P1/3, and the time to reach the peak decreases as P−2/3, where P is the Péclet number. The results of this study have important implications not only for ecological modeling, but for the general understanding of turbulent reactive flows. DOI: 10.1103/PhysRevE.74.016307 PACS number共s兲: 47.27.tb, 47.15.ki, 47.70.Fw Chemical mixing and reaction rates in turbulent flows play a crucial role in many fields including ecology 关1兴, atmospheric chemistry 关2,3兴, combustion 关4兴, and microfluidics 关5兴. The dual nonlinearities of fluid advection and chemical reaction rates combine to make this a particularly challenging problem. Most studies up to this point consider two scalars that are already in contact with one other in the initial condition. In one common scenario, the first scalar forms an island that is surrounded by the second scalar 关3兴; in another, each of the two scalars separately occupies one-half of the domain 关6兴. Results from this class of problems indicate that advective transport can enhance reaction rates between the two scalars, principally through filamentation that increases contact area between the neighboring scalars. In the present study, we investigate a different class of problems that has received little attention in the literature; namely, we consider two reactive scalars that are initially distinct, separated from one another by a third, nonreactive, scalar. This class of problems is inherently different because filamentation of the two reactive scalars by the advective field is not automatically accompanied by contact between those two species 共since the third scalar can act as a spatial barrier兲. These different problems can have different scaling behavior: the separation between the two scalars introduces a length scale, while the case of scalars filling the two halves of the domain has no intrinsic length scale. An ecological process called broadcast spawning motivated our study, and the process serves as an excellent example of this class of problems. Broadcast spawning is a fertilization strategy used by various benthic invertebrates 共e.g., sea urchins, anemones, corals兲 whereby male and female adults release sperm and egg into the surrounding flow 关7兴. The sperm and egg 共the two reactive scalars兲 are initially separated by ambient water 共the third, nonreactive, scalar兲. For effective fertilization, the turbulent stirring must not just produce filaments of the reactive scalars, it must also bring these filaments 共from the initially distinct scalars兲 into close contact. *Author to whom correspondence should be addressed. Electronic address: [email protected] 1539-3755/2006/74共1兲/016307共4兲 The aggregate effect of instantaneous turbulent stirring and mixing is often modeled using an eddy viscosity 共for momentum兲 or an eddy diffusivity 共for scalars兲. While this method is sometimes satisfactory, it is intrinsically a timeaveraged approach that can fail spectacularly. For example, in the broadcast spawning problem, field measurements of fertilization rates are rarely less than 5%, and are often as high as 90% 关8–10兴. However, numerical modeling based on a turbulent eddy diffusivity predicts fertilization rates of only 0.01–1 %, due to the strong time-average turbulent dilution of the gametes 关11兴. Many studies 关3,12兴 have shown that the instantaneous details of advective transport 共not captured by the eddy diffusivity approach兲 can enhance reactions rates for scalars that are initially in contact. We hypothesize that similar reaction enhancements are produced in problems like broadcast spawning, even though the reactive scalars are initially distinct. The results of this study have important implications not just for ecological modeling, but for the general understanding of turbulent reactive flows. By way of introduction, consider the turbulent transport of two scalars A and B. Each local scalar concentration CA and CB can be decomposed into mean and fluctuating components C共t兲 = C̄ + c⬘共t兲. For simplicity, we assume that CA and CB react at the second-order rate kCACB, where k is the dimensional reaction rate constant. At any location, the timeaveraged value of the reaction rate is then kCA共t兲CB共t兲 = kCACB + kcA⬘ 共t兲cB⬘ 共t兲, 共1兲 where the first term on the right-hand side is simply the product of the time-averaged concentrations, and the second term is the correlation between the fluctuating concentrations. By definition, models based on an effective turbulent diffusivity incorporate only the first term. The question then becomes, is the second term nonzero, and why? In this paper we examine the ability of structured stirring processes to impose scalar correlations on initially distinct point scalar masses. The basic mechanism by which advection can increase the reaction rate is simple, in principle. For incompressible flows, turbulent stirring can only stretch and fold scalar filaments, with no direct mechanism for decreasing the filament 016307-1 ©2006 The American Physical Society PHYSICAL REVIEW E 74, 016307 共2006兲 CRIMALDI, HARTFORD, AND WEISS FIG. 1. Geometry for two-scalar vortex stirring problem. concentration. Only molecular diffusion intersperses fluid between adjacent filaments to produce scalar dilution. Thus, the parametrization of advective stirring by an eddy diffusion is intrinsically flawed since it incorrectly produces a systematic decrease in concentrations and thus a decrease in reaction rates. If, on the other hand, advection brings two scalar filaments together before molecular diffusion decreases their concentration, then the reaction rate will be enhanced. Oceanic flows tend to be populated by coherent vortices with long lifetimes. We thus consider the problem of two reactive point-source scalars stirred by a vortical flow u␾共r兲 共see Fig. 1兲. The point masses are placed at distances RA and RB from the vortex center, with initial angular and linear separations ␾ and L. Both scalars have mass M and molecular diffusivity D. We consider only the low-concentration limit where the scalar concentrations are small and the second-order reactive removal of scalar mass is negligible. The total dimensional reaction rate ⍜̂共t兲 corresponds simply to the scalar overlap integrated over the domain: ⍜̂共t兲 = 冕冕 kCACBdA. 共2兲 To establish a baseline case that excludes instantaneous stirring processes, we first consider the case of pure diffusional spreading 共u␾ = 0兲 due to molecular diffusion. We work in nondimensional variables where we scale time by a diffusive time td = L2 / 8D, distances by L, and the total reaction rate by kM 2 / ␲eL2. The resulting expression for ⍜ in the low-concentration limit is 关13兴 ⍜= 冉 冊 1 e exp − . t t 共3兲 With this scaling, ⍜ reaches its peak value of unity at t = 1. Equation 共3兲 is shown as the dashed line in Fig. 2. For small times, ⍜ is effectively zero due to the initially distinct scalars, and, at large times, ⍜ asymptotes to zero due to diffusional dilution of CA and CB. It is only at intermediate time scales that diffusion has produced enough spreading for interactions between the scalars, but not so much spreading that they become overly dilute. To consider the effect of structured stirring processes, we now add an ideal point-vortex velocity field u␾共r兲 = ⌫ / 共2␲r兲, where ⌫ is the vortex circulation. The advective time scale based on one revolution at r = R ⬅ 共RA + RB兲 / 2 is ta = 2␲R2 / ⌫. This leads to a Péclet number FIG. 2. Numerical simulations of ⍜ vs t for P = 0 , 10, 100. The P = 0 results match Eq. 共3兲, which is shown with a dashed line. 冉 冊 td L ⌫ = ta 4␲D 2R P= 2 共4兲 , where P = 0 corresponds to the pure diffusion case and P becomes large as advection dominates. We first determine ⍜ for the simplest case RA = RB = L / 2 and ␾ = ␲, such that P = ⌫ / 4␲D. In the limit of large P, ⍜ can be calculated analytically. Due to the radial shear of the vortex, particles are rapidly dispersed in the circumferential direction by shear dispersion 关14兴. For large P, the radial diffusion is small, the nonlinear dependence of the velocity on radius is well approximated by its linear shear, and the polar geometry 共r , ␾兲 can be approximated by a periodic Cartesian geometry x = 共x , y兲. The velocity is then u = 共 −␴y , 0兲, and the shear is ␴ = ⌫ / ␲R2 = 16PD / L2. Thus, in the large-P limit, the problem depends only on the local shear, independent of other aspects of the vortex profile. In a nonperiodic Cartesian geometry an initial ␦-function point source at the origin evolves under a linear shear into a Gaussian concentration C共x , t兲 = N exp共xA−1x / 2兲, where A is the covariance matrix, A= 冉 P2t3/3 − Pt2/4 − Pt2/4 t/4 冊 共5兲 , and N is a normalization factor. Note that the P2t3 shear dispersion in the x direction dominates over the linear-intime dispersion in the y direction. For the periodic Cartesian geometry intrinsic to vortex stirring, the periodic images of C must be summed. The sum in the resulting integral, Eq. 共2兲, gives the Jacobi elliptic theta function ␽2共z , q兲 关15兴. The reaction rate for different Péclet numbers is self-similar for large P, and the explicit dependence on P can be eliminated by using the similarity variables ⍜s = ⍜P−1/3, and ts = tP2/3, giving ⍜s共ts兲 = 冑3e ts2 冉 冋 册冊 ␽2 0,exp − 3␲2 ts3 . 共6兲 Equation 共6兲 is shown as the dashed line in Fig. 3. The peak reaction rate grows as P1/3, while the time to reach the peak shrinks as P−2/3, indicating that vortex stirring can produce significant amplification and acceleration of the reaction over pure diffusion. A similar analysis for two scalars filling two halves of an infinite domain gives, due to the lack of a length 016307-2 PHYSICAL REVIEW E 74, 016307 共2006兲 REACTION ENHANCEMENT OF POINT SOURCES DUE¼ FIG. 3. Simulation results for large Péclet numbers, plotted in self-similar variables. As P increases, the simulation results approach Eq. 共6兲, shown with the dashed line. scale, a similarity solution where different values of shear and diffusion give the same solution but at a different time, which is qualitatively different from the scaling with P seen here for two point sources. Turbulent stirring is often parametrized by eddy diffusion. If stirring does act diffusively, then one can choose an eddy diffusivity De so the diffusive time scale for the eddy diffusivity matches the advective time scale, De = PD. The reaction rate for a process governed by such an eddy diffusion is the diffusive reaction rate Eq. 共3兲, with time scaled by P: ⍜共t兲 = 冉 冊 1 e exp − . Pt Pt as P increases are evident. The P = 0 simulation matches the pure-diffusion analytical solution 关Eq. 共3兲兴, shown by the dashed line. Simulation results for higher Péclet numbers are shown in Fig. 3, this time plotted in self-similar coordinates. As P increases, the simulation results approach the highPéclet analytical solution 关Eq. 共6兲兴, shown with the dashed line. In order to test the effects of the initial geometric placement of scalars A and B on the reaction rate, we now vary the initial geometry of the scalars relative to the vortex. The relative radii of the initial scalar positions is described by an eccentricity 共7兲 The peak reaction rate is still unity, but occurs at an accelerated time which scales as P−1. Eddy diffusion can thus speed up reactions, but since it necessarily spreads scalars as they diffuse, it cannot enhance the reaction rate over the values attained under molecular diffusion. This illustrates the fundamental failure of eddy diffusion models of turbulent stirring of reactive scalars. With analytical ⍜ expressions now defined for P = 0 关Eq. 共3兲兴, and for large P 关Eq. 共6兲兴, we now turn to a numerical modeling approach to calculate ⍜ for a range of intermediate P values. The model uses a Lagrangian particle-tracking approach, with a random-walk model for molecular diffusion. The two-dimensional position of a particle xn is determined from its previous position xn−1 using the updating relationship 关16兴 xn = xn−1 + u共xn−1兲⌬t + Z冑2D⌬t, FIG. 4. Simulations of ⍜ for eccentric 共E ⫽ 0兲 initial conditions at P = 100 and ␾ = ␲. 共8兲 where u is the local velocity vector, and Z is a twodimensional random vector with mean magnitude zero and variance one. The particle-tracking method does not exhibit any numerical dispersion in the classical sense 关16兴, and the computational effort is proportional to the number of particles and is therefore concentrated into areas with the highest concentrations 关17兴. In the simulation, we vary P, but hold kM / D 共a Damkoler number兲 constant at unity. Simulation results for ⍜ vs t for P = 0, 10, and 100 are shown in Fig. 2, again for the simplest case RA = RB = L / 2 and ␾ = ␲. The enhancement and acceleration of the reaction rate E= RA − RB , RA + RB 共9兲 where E = 0 is a symmetric initial condition, and E = 1 is the limiting case where one scalar mass coincides with the vortex center. Simulation results for ⍜ for six values of E at P = 100 and ␾ = ␲ are shown in Fig. 4. As E increases from zero, the peak value of ⍜ first increases 共as seen in the E = 0.15 curve兲 and then decreases. As E approaches unity, the ⍜ curve reverts to the pure diffusion solution. For intermediate values of E, two local ⍜ maxima occur. The first maximum corresponds to an asymmetric distribution of scalar overlap 关the integrand of Eq. 共2兲兴 as the leading edge of one of the scalar filaments overtakes the trailing edge of the second 共see Fig. 5兲. The second maximum corresponds to a largely symmetric overlap distribution as radial diffusion dominates the scalar dispersion. The nondimensional reaction rate ⍜ is insensitive to changes in the initial angle ␾ for the pure diffusion case, but it becomes increasingly sensitive as P increases. The variation of the peak value of ⍜ as a function of ␾ is shown for three values of P in Fig. 6共a兲. Finally, we consider the significance of the chosen vortex velocity distribution u␾共r兲 on ⍜. An interesting comparison case is the distribution given by u␾共r兲 = 共⌫ / 2␲r兲兵1 − exp关 −共2r / B兲2兴其. Often called an Oseen vortex, this distribution behaves like the ideal point vortex when 2r / B Ⰷ 1, and like solid-body rotation when 2r / B Ⰶ 1. Figure 6共b兲 is a plot of the peak value of ⍜ in an Oseen vortex for P = 1000 and ␾ = ␲ as a function of initial scalar separation L / B. When the scalars are placed outside the rotational core of the vortex 016307-3 PHYSICAL REVIEW E 74, 016307 共2006兲 CRIMALDI, HARTFORD, AND WEISS FIG. 5. 共Color online兲 Snapshots of a simulation for P = 1000, E = 0.16, and ␾ = ␲ showing 共a兲 scalar A 关red 共light gray兲兴 and B 关blue 共dark gray兲兴 particle locations, 共b兲 distribution of the scalar overlap 关the integrand of Eq. 共2兲兴, and 共c兲 the time 共shown by the vertical line兲 of the snapshots. 共L / B Ⰷ 1兲, ⍜ has the same enhancement and acceleration as the ideal vortex result. However, when the scalars are placed inside the core, where there is no straining 共and hence no shear dispersion兲, ⍜ reverts to the unenhanced pure diffusion result. The results presented in this paper demonstrate that instantaneous stirring processes have a profound impact on reaction rates between two initially distinct point scalars. Stirring imparts spatial correlations from the velocity field onto the initially uncorrelated scalar field. This structured stirring results in enhanced reaction rates that are not predicted by eddy diffusivity models, even when the average scalar variance growth is modeled correctly. The reaction enhancement is highest at intermediate time scales, when scalar filaments are stirred into close proximity but are not yet strongly diffused. Any eddy diffusivity approach 共regardless of the diffusivity magnitude兲 will produce the same reaction rates as the P = 0 pure diffusion case shown in Fig. 2, with no reaction enhancement. Thus, the results indicate a fundamental 关1兴 A. Okubo and S. A. Levin, Diffusion and Ecological Problems: Modern Perspectives, 2nd ed. 共Springer-Verlag, New York, 2001兲. 关2兴 S. F. Maria, L. M. Russell, M. K. Gilles, and S. C. B. Myneni, Science 306, 1921 共2004兲. 关3兴 S. Edouard, B. Legras, and V. Zeitlin, J. Geophys. Res., 关Atmos.兴 101, 16771 共1996兲. 关4兴 N. Peters, Turbulent Combustion 共Cambridge University Press, Cambridge, U.K., 2000兲. 关5兴 V. Hessel, H. Lowe, and F. Schonfeld, Chem. Eng. Sci. 60, 2479 共2005兲. 关6兴 F. J. Muzzio and M. Liu, Chem. Eng. J. 64, 117 共1996兲. 关7兴 D. Levitan, in Ecology of Marine Invertebrate Larvae, edited by L. McEdward 共CRC Press, Boca Raton, FL, 1995兲, pp. 123–156. FIG. 6. Sensitivity of the peak value of ⍜ to 共a兲 variations in the initial angle ␾ for the ideal vortex for P = 10, 100, and 1000, and 共b兲 variations of the initial scalar separation L / B in an Oseen vortex for ␾ = ␲ and P = 1000. The dashed lines indicated the peak ⍜ values for an ideal point vortex at P = 0 and 1000. flaw in relying on time-averaged techniques such as eddy diffusion when modeling these types of reactions. Although we chose an extremely simplistic velocity field 共to facilitate both analytical and numerical approaches兲, the results suggest similar implications for more generalized flowfields. If turbulence is conceptualized as a collection of interactive vortices, then, at local scales, the instantaneous flow fields resemble those used in this study. In this case, two reactive scalars in the vicinity of a vortex are likely to experience reaction enhancement and acceleration analogous to that presented in the study. The relative insensitivity of the results to the initial geometry of the scalars facilitates this extension. Indeed, preliminary simulations of twodimensional flows with multiple, interactive vortices exhibit enhanced reaction rates relative to that predicted by an eddy diffusivity approach. Detailed studies are needed to assess possible extensions of these results in two- and threedimensional flows. The authors would like to thank M. Denny for ecological insights to this problem. This work is supported by the National Science Foundation under CAREER Grant No. 0348855. 关8兴 关9兴 关10兴 关11兴 关12兴 关13兴 关14兴 关15兴 J. E. Eckman, J. Exp. Mar. Biol. Ecol. 200, 207 共1996兲. J. Pennington, Biol. Bull. 169, 417 共1985兲. P. Yund, Trends Ecol. Evol. 15, 10 共2000兲. M. W. Denny and M. F. Shibata, Am. Nat. 134, 859 共1989兲. J. M. Ottino, Chem. Eng. Sci. 49, 4005 共1994兲. J. P. Crimaldi and H. S. Browning, J. Mar. Syst. 49, 3 共2004兲. G. Taylor, Proc. R. Soc. London, Ser. A 219, 186 共1953兲. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, edited by M. Abramowitz and I. A. Stegun 共Dover Publications, New York, 1974兲. 关16兴 W. Kinzelbach, in Groundwater Flow and Quality Modeling, edited by E. Custodio 共Reidel, Rotterdam, 1988兲, pp. 227-245. 关17兴 K. Dimou and E. Adams, Estuarine Coastal Shelf Sci. 37, 99 共1993兲. 016307-4