Quantitative Analysis of Mixing Structures in Chaotic Flows Generated

5722

J. Phys. Chem. A 2002, 106, 5722-5736

Quantitative Analysis of Mixing Structures in Chaotic Flows Generated by Infinitely Fast Reactions in the Presence of Diffusion M. Giona,* S. Cerbelli, and A. Adrover Dipartimento di Ingegneria Chimica, Centro InteruniVersitario sui Sistemi Disordinati, e sui Frattali nell’Ingegneria Chimica, UniVersita` di Roma “La Sapienza”, Via Eudossiana 18, 00184 Roma, Italy ReceiVed: October 10, 2001; In Final Form: January 29, 2002

The interplay between diffusive and convective mixing processes may have a strong impact upon apparent reaction rates. This paper analyzes the interaction of convection and diffusion mechanisms by considering an infinitely fast irreversible reaction A + B f products, occurring in a two-dimensional chaotic flow. Attention is focused on the geometric properties of mixing patterns and on the overall reactant consumption. We show that the length of the reaction interface undergoes a transition from a kinematics-dominated exponential growth to a persistent oscillatory regime. This regime results from two competing mechanisms, namely, recursive stretching and folding of the interface caused by chaotic advection and merging of contiguous striations patterns owed to diffusive transport. In the case of globally chaotic flows, a singular transition is observed in the scaling of the dominant eigenvalue with the Peclet number. The geometric information arising from the analysis of the reaction interface is also exploited for deriving a simple one-dimensional model that predicts the apparent rates over a wide range of Peclet numbers.

1. Introduction The possibility of obtaining high reaction yields in due time in large scale volumes relies heavily upon the efficiency by which the reacting species are brought into contact with one another.1 A particularly critical case is represented by highly viscous reactants that are mechanically stirred within industrial equipment. In this case, the stirring process is designed to accomplish a 2-fold task: (a) reduce the segregation length down to the diffusive lengthscale and (b) generate interface between the reactants by stretching and folding the reaction zone into a foliated continuous lamellar structure. If the characteristic time of reaction is short compared to those of convection and diffusion (mixing-controlled reaction), then the reaction zone is reduced to a two-dimensional interface that acts as a boundary between the segregated reactants,2,3 (see section 3). Besides its unquestionable practical relevance, this physical setting deserves particular attention as it bridges the notion of mixing structure from a purely kinematic (i.e., convective) frame to systems with molecular diffusion. The nonlinear dependence of nearly all of the (bounded) velocity fields W(x, t) of practical interest upon the spatial coordinates makes closed-form solutions of the advectiondiffusion equation, with or without chemical reactions, generally unattainable. For this reason, research efforts have been directed toward the characterization of the mixing action due to convection alone (i.e. in the absence of diffusion), assuming that the effects owed to diffusion could be somehow superimposed a posteriori to the purely kinematic picture. (A conceptually analogous approach has also been applied to the advectiondiffusion of a vector quantity (i.e., the magnetic field) in the so-called fast dynamo problem, by introducing the concept of pulsed system.4) * To whom correspondence should be addressed.

In the diffusionless case, point tracers move according to the kinematic equation

x3 ) v(x, t)

(1)

where x is the tracer position. The solutions of eq 1 (i.e., the point tracer trajectories) may exhibit remarkably complicated features even when the flow field takes on very simple functional forms. This phenomenon, generically referred to as Lagrangian chaos, has been investigated in a wealth of systems, encompassing prototypical models as well as industrially relevant flows, to define what the essential features that make a flow field an efficient mixer are. (see, e.g., refs 5-7 and references therein). Motivated by the broader perspective of understanding the interplay between convection-diffusion and chemical reactions, more recent work has focused on the mechanisms of deformation of material interfaces as opposed to single particle trajectories. It has been found that material interfaces, advected by laminar chaotic flows, possess global invariant properties that can be directly connected to the existence of invariant manifolds within the chaotic region (a concise review of the main results concerning chaotic advection is reported in section 3). In parallel with this approach, several articles have focused on the analysis of simplified premixed systems (lamellar systems) undergoing finite-rate and infinitely fast reactions.8-10 These works are explicitly or implicitly based on a time-splitting between the action of convection and that of diffusion, in the meaning that mechanical agitation was accounted for only in defining the initial conditions (i.e., the spatial distribution of a one-dimensional array of lamellae of various thicknesses). Starting from this “frozen” picture, the modification induced by diffusion on the segregation patterns, within the otherwise still medium, were theoretically and numerically investigated in order to derive quantitative information on the overall reaction rate and product formation. The results obtained in the above

10.1021/jp013781e CCC: $22.00 © 2002 American Chemical Society Published on Web 05/16/2002

Analysis of Mixing Structures in Chaotic Flows framework were extended to include the simultaneous action of convection by using the concept of “warped time”,2,3 that is, by assuming that the medium is subjected to a convective field equivalent to the flow along the stable direction of a hyperbolic stagnation point. This assumption is motivated by the fact that, in a chaotic flow, lamellae are actually shrunk at an exponential rate to maintain incompressibility, whereas they are locally stretched along the unstable directions of the flow (see section 3). In point of fact, the complexity of convective mixing manifests itself through global features, namely, the folding of lamellae, which brings into close contact portions of the fluid that were not contiguous in the one-dimensional approximation of the mixing space within the lamellar approach. This mechanisms contrasts the action of lamellar merging by molecular diffusion, by introducing continuous generation of newborn lamellae, which enter the one-dimensional space from the “orthogonal direction”, i.e., from the dilating direction of the flow. A direct approach to simultaneously convecting-diffusing systems (without chemical reaction) in the context of twodimensional periodically forced chaotic flows has been undertaken in refs 11-13. In particular, Rom-Kedar and Poje11 focus on the impact of the frequency of velocity oscillations upon the maximum mass flux and find that a “Lagrangian steady state” is established as a result of the combined action of advection and diffusion. Tang and Bozer12 reformulate the advection-diffusion equation in a Lagrangian framework. Adrover et al.13 analyze the influence of chaotic advection upon the diffusive transport of a tracer concentrated along a one-dimensional subset of the mixing space. By tracking the structure of the concentration contours, two global quantities related to the structure of the concentration field are defined, namely, the diffusional thickness and the area of diffusional influence. It is shown that these quantities undergo nonmonotonic behavior with a transition from a kinematics-dominated template to a complex oscillatory state. Numerical simulations of advection-diffusion-reaction kinetics in cellular and laminar chaotic flows have been performed by Reidaga et al.,14,15,16 by Liu and Muzzio,17 and by Zalc and Muzzio,18 for finitely fast bimolecular reactions (parallel and consecutive, respectively). Reigada et al. focus mainly on the scaling of the reactant quantities with time in order to determine the controlling regime from characteristic scaling exponents, whereas Liu and Muzzio and Zalc and Muzzio analyze the effects of the flow protocol on the overall product formation and on the spatial heterogeneity of the product distribution. By enforcing the analogy between heat and mass transfer, it is worth mentioning the results by Sawyers et al.,19,20 by Mokrani et al.,21 and by Raynal and Gence22 on heat transfer in laminar chaotic flows, which focus on the effects of chaotic stirring on the performances of heat exchangers. From this concise overview of the state of the art, three main questions arise: (a) how does the presence of molecular diffusion modify the geometry of partially mixed structures caused by passive advection, (b) do invariant geometrical patterns exist, arising from the interplay between chaotic advection and diffusion, and finally (c) how can the geometrical characterization derived from points a and b be applied for the prediction of the overall conversion in the presence of chemical reactions? Throughout this article, we attempt to give an answer to these fundamental issues by considering the dynamics of reaction

J. Phys. Chem. A, Vol. 106, No. 23, 2002 5723 interfaces in the case of infinitely fast (instantaneous) reactions. The application of Gale¨rkin (spectral) methods provides the most convenient computational strategy for approaching the resulting balance equations, in that they allow to bypass the spurious numerical diffusion associated finite difference/volume/element simulations.13 The main goals of this article can be summarized as follows: (i) illustrate qualitatively and quantitatively the geometric features of the partially mixed structures generated by an infinitely fast reaction in two-dimensional chaotic flows, (ii) show the occurrence of persistent/invariant patterns in the dynamics of reaction interfaces, intrinsically controlled by the presence of diffusion, the existence of which can be proved from the analysis of the corresponding advection/diffusion equation regarded as a dynamical system evolving in an abstract functional space, (iii) provide numerical evidence that, however small the diffusivity, the structure of invariant patterns in the presence of diffusion differs from that associated with pure advection (the existence of this singular behavior demands that mixing indices, aimed at establishing the quality of the stirring protocol in systems with molecular diffusion, be derived from the analysis of the advection-diffusion equation), and (iv) use the information arising form the geometric analysis to derive a simple one-dimensional model for predicting overall reaction rates and product formation. The article is organized as follows. Section 2 formulates the mathematical setting of the problem. Section 3 reviews succinctly the geometric properties of interface dynamics in chaotic flows in the diffusionless setting. Section 4 defines the flow systems and the numerical techniques used to approach instantaneous reactions and to track the reaction interface. Results on the time evolution of the reaction interface length L(t) and on the spatial interface structure are analyzed in section 5 for a wide range of Peclet numbers and for different mixing protocols. This section also exploits the functional setting of the advection-diffusion equation in order to explain the phenomenology observed numerically. Section 6 discusses the singular limit of the advection/diffusion equation for vanishing diffusivities in the case of globally chaotic flows and addresses thoroughly the asymptotic decay of reactant quantities. Finally, in section 7, the knowledge of overall interface length is used for deriving a one-dimensional model that predicts the overall reaction rate over a wide range of the Peclet number. 2. Statement of the Problem Throughout this article, we consider a bimolecular reaction A + B f product that occurs in a incompressible medium flowing with velocity v(x, t) in a bounded region M. Flow incompressibility dictates the velocity field be solenoidal, i.e., ∇‚v ) 0. We make the following assumptions about the system: (i) the reaction is irreversible, i.e., it proceeds until complete consumption of the limiting reactant, (ii) the physical properties of the mixture (density, viscosity, etc.) do not depend on the mixture composition, and (iii) reactants A and B are characterized by identical diffusivities within the mixture, DA ) DB ) D ) constant. Under these conditions, the mass balance equations for the reactants read

∂CA + v‚∇CA ) D∆CA - kCACB ∂t

(2)

∂CB + v‚∇CB ) D∆CB - kCACB ∂t

(3)

5724 J. Phys. Chem. A, Vol. 106, No. 23, 2002

Giona et al.

where CA and CB are the molar concentrations of the reacting species, ∆ denotes the Laplacian operator, and k is the rate constant of the reaction. Equations 2 and 3 are equipped with the initial conditions CA(x, t)|t)0 ) CA,0(x) and CB(x, t)|t)0 ) CB,0(x). Equations 2 and 3 can be made dimensionless by rescaling concentration, length, and velocity with some reference values Cref, L, and Vref. This induces a rescaling of the time scale by a factor L/Vref. Henceforth, we consider exclusively the dimensionless formulation associated with eqs 2 and 3, and we use the same symbols CA, CB, x, v, and t to indicate the corresponding dimensionless quantities. By setting φ ) CA - CB, and by subtracting eq 3 from eq 2, one obtains a linear advection-diffusion equation for the concentration difference φ:

∂φ 1 + v‚∇φ ) ∆φ ∂t Pe

(4)

where the Peclet number, Pe ) VrefL/D, is the dimensionless group expressing the ratio of the characteristic time for diffusion to that of advection. The diffusionless limit is hence given by Pe f ∞. When the characteristic time of reaction is much shorter than those of convection and diffusion (infinitely fast reaction, k f ∞), the concentration product CACB vanishes within all of the flow domain except on the interface between A and B, that is, on the reaction interface.2 Thus, reactants remain segregated at all times, and the evolution of the system expressed by eqs 2 and 3 is completely specified once the solution φ(x, t) of eq 4 is known, because

CA(x, t) )

φ(x, t) + |φ(x, t)| , 2 CB(x, t) )

-φ(x, t) + |φ(x, t)| (5) 2

The zero-level set of the φ function

γ0(t) ) {x| φ(x, t) ) 0}

(6)

identifies the reaction interface separating the species A and B at time t. The asymptotic persistence of the reaction interface depends on the initial loading conditions, that is, on the global reactant quantities

mR(0) )

∫MCR,0(x) dx,

R ) A, B

(7)

initially present in the system. If Φ0 ) mA(0) - mB(0) * 0, one of the two reactants (B if Φ0 > 0, A otherwise) disappears in finite time. Conversely, in the case where Φ0 ) 0 (stoichiometric loading conditions), neither of the two reactants is completely consumed in finite time, and the reaction interface γ0(t) is a nonempty set of points for all positive values of t. Unless otherwise specified, in the remainder of this paper, we will focus exclusively on stoichiometric loading, to be able to perform an asymptotic analysis of reaction interface and segregation patterns. In a broader context, the analysis of the level set γ(t), where

)

∫Mφ(x, 0)

1 mis(M)

(8)

where mis(M ) denotes the measure (volume, area) of the mixing space, provides an intrinsic geometric characterization of the

advection/diffusion equation eq 4, independently of the physical meaning of eq 4 within the context of infinitely fast irreversible reactions. 3. Diffusionless Setting: Chaotic Advection of Material Interfaces In the limit of vanishing diffusivity, eq 4 becomes

∂φ Dφ + v‚∇φ ) )0 ∂t Dt

(9)

where D/Dt denotes the material derivative, i.e., the derivative taken along the trajectory of a fluid particle that moves with velocity v(x, t). Assume the two reactants A and B are initially segregated within the mixing space M into two disjoint regions A, B with A ∪ B ) M. Let their concentration be uniform within each region, i.e., CA(x, 0) ) CA,0 if x ∈A, CA(x, 0) ) 0 elsewhere, and CB(x, 0) ) CB,0 if x ∈B, CB(x, 0) ) 0 elsewhere. The difference φ ) CA - CB takes values φ ) CA,0 if x ∈A and φ ) -CB,0 if x ∈B, and it is discontinuous with a jump CA,0 - CB,0 at the boundary ∂A ) ∂B ) γ0(0) between the two regions, with γ0(0) being a curve or a closed surface, depending on whether the flow is two- or three-dimensional. Because the solution of eq 9 is constant along the flow determined by v(x, t), the knowledge of the interface γ0(t) determined by advecting the points x ∈γ0(0) through eq 1 specifies completely the structure of the concentration field at any given time t. (Clearly, being D ) 0, the mass of both reactants is conserved, as they cannot be transported to the interface to react with one another.) As pointed out in the Introduction, the qualitative and statistical features of the solutions of eq 1 are nontrivial even when the structural form of the velocity field is simple. The intense research activity in the field of chaotic advection has produced a wealth of results and observations, the majority of which applies to laminar flows, that constitute the simplest physical frame amenable to direct investigation. In particular, theoretical,5,7,26 computational,27,28 and experimental6,29 studies focusing on two-dimensional time-periodic flows proved that, in the limit of vanishing diffusivity, best mixing performance is achieved when the stirring field is chaotic (i.e., when there exists a subset of positive measure within the mixing space characterized by a positive Lyapunov exponent). From the geometric viewpoint, this condition ensures that a generic segregated region be recursively stretched and folded by the stirring flow toward a continuous yet recursively nested filamented structure that invades densely all of the chaotic subregion of the mixing space.23-25 The combined action of stretching and folding causes a sustained exponential growth of the intermaterial contact area, along with exponential shrinking of lengths in the transverse direction, with the shrinking being a consequence of flow incompressibility. Most importantly, interface dynamics possesses invariant properties in the meaning that the geometric structure attained by the interface at integer multiples of the period of the flow becomes independent of both time and initial condition, as the interface is progressively transformed into one of the unstable leaves associated with the Poincare` map of the flow. An example of these properties is given in Figure 3 parts I and II and Figure 4 parts I and II for the Sine Flow. Details about the flow system an the initial interface structure are given in section 4. For the time being, we only want to point out how the interface quickly develops into a convoluted curve that fills densely the chaotic region of

Analysis of Mixing Structures in Chaotic Flows

J. Phys. Chem. A, Vol. 106, No. 23, 2002 5725

Figure 2. Poincare´ section for the Prototypical Cavity Flow at T ) 0.6. The figure shows a segment of 103 iterates of a hundred particles that were initially uniformly distributed along the line y ) 1/2. The phase space landscape shows a main chaotic region intertwined with quasiperiodic islands.

Figure 1. Poincare´ section of the Sine Flow system with (a) T ) 0.4, (b) 0.8, and (c) 1.6. It can be observed how the size of quasiperiodic islands is progressively shrunk as the flow period increases.

the mixing space while maintaining the overall time-invariant shape. As a macroscopic outcome of this growth template, the overall length of the curve increases exponentially in time. Another interesting aspect that emerges from the analysis of chaotic advection is that all of the qualitative features described above are generic in the meaning that they are shared by a large class of nonlinear incompressible flows. From the point of view of mixing efficiency, the enhancement of the overall reaction rate caused by chaotic advection is accomplished through two different mechanisms, namely, the amplification of concentration gradients by the shrinking of lengths in the stable direction of the chaotic region and the stretching of the interface along the unstable directions. On the other hand, to have a nonvanishing rate of reactants consumption, a positive molecular diffusivity must be present in the system, which modifies drastically the mixing patterns as it will be discussed in section 5.

Figure 3. Comparison of the kinematic (Pe ) ∞) and reaction interfaces at Pe ) 104 in the Sine-Flow at T ) 1.6 (the globally chaotic protocol) (a-d) mixing patterns (white and gray) and reaction interface (black line) at times T, 2T, 3T, and 5T. (I and II) Kinematic interface at the end of the first (I) and second (II) period.

4. Flow Systems and Numerical Techniques As a model system, we consider a well-know paradigm of chaotic behavior, the Sine Flow system (SF), obtained by blinking every T/2 time units the two steady fields v1 ) (sin(2πy), 0) and v2 ) (0, sin(2πx)).30,24 This model flow is defined


Giona et al. where i ) (-1)1/2. By substituting this expression into eq 4, the following infinite-dimensional linear tri-diagonal system of ODEs for the coefficients during the first half-period T/2 (motion along the x-axis) is obtained:

φ˙ h,k ) -

Figure 4. Comparison of the kinematic (Pe ) ∞) and reaction interfaces at Pe ) 104 in the Sine-Flow for a mixing protocol with large quasiperiodic islands (T ) 0.4). (a-d): mixing patterns (white and gray) and reaction interface (black line) at times 20T, 30T, 50T, and 60T. (I-II): Kinematic interface at time 20T (I) and 30T (II).

on the two-dimensional torus (i.e., on the unit square M ) I2 equipped with periodic boundary conditions). While the geometric and statistical properties of this type of flow have been shown to be similar to those of physically realizable systems, the absence of boundaries in the mixing space simplifies remarkably the numerical solution of eq 4, allowing us to perform accurate simulations at feasible computational expenses. Figure 1 depicts the Poincare` section (i.e., the superimposed plot of the particle trajectories at multiple integers of the flow period, starting from different initial positions) associated with the kinematic equation of motion eq 1 for three values of the period. The case T ) 0.4 (Figure 1a) is characterized by two large islands of quasiperiodicity, which surround the chaotic region, identified by the X-shaped random cloud of points. At T ) 0.8 (Figure 1b), the spatial extension of quasiperiodic islands is significantly reduced, and for T ) 1.6 (Figure 1c), the chaotic region appears to invade the entire mixing space. The most convenient approach for solving numerically eq 4 is the Gale¨rkin (spectral) expansion of the function φ, because it is intrinsically free of numerical diffusion, which is instead unavoidably associated with finite differences (FD), finite volume (FV), and even finite elements (FE) methods.31 The spurious effects of numerical diffusion can be overcome by using a spectral expansion, i.e., an integral approximation based on a generalized Fourier series expansion with respect to the complete system of eigenfunctions of the Laplacian operator. Owing to the periodic structure of the torus and of the velocity field defining the SF system, this representation coincides with the classical Fourier series expansion: ∞

φ(x, t) )

∑ φh,k(t) exp[2πi(hx + ky)] h,k)-

(10)

4π2 2 (h + k2)φh,k - πh(φh,k-1 - φh,k+1) (11) Pe

where φ˙ h,k ) dφh,k/dt. An analogous equation is obtained for the second half-period by interchanging the indices h and k. Equation 11 was solved numerically for 102 e Pe e 105 with h, k ∈ - N, ..., N, where N (the number of modes is (2N + 1)2) varied from 80 (for the lowest Pe) to 300 (for Pe ) 105). The number of modes was chosen so as to ensure N independence of the solution in the norm L2. For the time-integration, we used both an implicit second-order scheme and a fourth-order Runge-Kutta algorithm. The two methods yield basically identical results. To determine both the spatial patterns, and the overall interface length, the values of φ(x, t) on a square mesh of 512 × 512 nodes were computed by means of a standard fast-Fourier transform (FFT) routine. The use of the FFT algorithm is essential in that a direct computation the nodal values of φ(x, t) through eq 10 would result in exceedingly long CPU time. Once the function φ(x, t) was computed on a square grid of points, its nodal values were linearly interpolated, and the reaction interface γ0(t) was determined as the intersection of the globally continuous piecewise triangular surface with the plane φ(x, t) ) 0. The mesh utilized proved fine enough to resolve the details of mixing patterns for all of the finite Pe values considered. To verify the generality of the results derived, we also consider the case of a wall-bounded two-dimensional flow, the prototypical cavity flow (PCF). While the PCF closely mimics the structure of the Stokes flow within a rectangular cavity with moving walls (cavity flow), it allows us to avoid the interpolation problems that are associated with the numerical solution of the Stokes equation in a rectangular cavity. The PCF stems from the following definition of the flow stream function:

ψ(x, y) ) ψo sin2

() ( ) πx πy2 sin 2 Lx Ly

(12)

The resulting velocity field is given by

Vx ) -

Vy )

() ( ) () () ( )

2πψoy 2 πx ∂ψ πy2 )sin cos ∂y Lx Ly2 Ly2

∂ψ 2πψo πx πx πy2 sin cos sin 2 ) ∂x Lx Lx Lx Ly

(13)

The stream function ψ(x, y) defined by eq 12 satisfies the boundary conditions for a cavity flow system in which the upper wall parallel to the x axis is moving and the other walls are static: The velocity field associated with the motion of the

Vy|x)0 ) Vy|x)Lx ) Vy|y)0 ) Vy|y)Ly ) 0 Vx|x)0 ) Vx|x)Lx ) Vx|y)0 ) 0, Vy)Ly )

( )

2πψo 2 πx sin Ly Lx

(14)

bottom wall can be readily obtained from eq 12 by enforcing the symmetries Vxbw(x, y) ) -Vxtw(x, Ly - y) and Vybw(x, y) ) Vytw(x, Ly - y), where vtw and vbw are the velocity fields



associated with the motion of the top and bottom walls, respectively. We consider the case ψo ) 1, Lx ) Ly ) 1, so that the period T is the parameter characterizing the stirring protocol. Figure 2 shows the Poincare` section of the PCF flow at T ) 0.6. This value of the period shows the occurrence of a main central chaotic region surrounded by smaller invariant chaotic sets, intertwined with islands of quasiperiodic motion. The difference field in the square cavity flow (Lx ) Ly ) 1) satisfies the zero-flux boundary condition at the cavity walls

∂φ ∂φ ) 0, )0 | | ∂x x)0,Lx ∂y y)0,Ly

(15)

The system of eigenfunctions of the Laplacian operator equipped with these boundary conditions {Rn,m(x, y)} is given by

Rn,m(x, y) ) cos(nπx) cos(mπy), n,m ) 0, 1, ... (16) Likewise for the SF system, the advection/diffusion equation can be reduced to a system of ordinary differential equations by enforcing a Gale¨rkin expansion with respect to the orthogonal system {Rn,m(x,y)}. Details are not reported here, because the resulting equations are lengthy and their derivation is redundant for the main goals of this article. 5. Reaction Interface Dynamics As an initial condition to investigate the evolution of interface dynamics, we chose

φ(x, 0) ) 2 - 4η(x - 1/2)

(17)

where η(x) is the unit step function (η(x) ) 1 for x > 0, and η(x) ) 0 for x > 0), corresponding to a unit initial mass of the reactants that are spatially organized into two vertical stripes each of which occupies half of the unit square I 2. For the SineFlow system, the boundary between reactants is given by the two segments x ) 0 and x ) 1/2 (0 < y < 1) as the two vertical boundaries of the unit square x ) 0 and x ) 1 are identified on the torus. The kinematic interface in the limit Pe ) ∞ is obtained by advecting the points of the segment through eq 1. Figure 3a-d shows the segregation patterns (white and gray) and the reaction interface γ (black line) in the SF for a mixing protocol specified by a period T ) 1.6 at Pe ) 104 at times nT, with n ) 1, 2, 3, and 5, respectively. In the same figure, the structure of the kinematic interface (Pe ) ∞) is shown at the end of the first (I) and second (II) periods. It is worth noting that the value T ) 1.6 of the switching period yields a nearly globally chaotic protocol (compare with Figure 1). In this situation, the stroboscopic evolution of the kinematic interface undergoes invariant space-filling exponential growth, as can be observed from Figure 3 parts I and II. Snapshots of the kinematic interface at later times (not shown here for brevity) yield what is essentially the same structure supplemented with increasingly fine detail. Comparison of Figure 3 parts a and I shows that at the end of the first period the reaction interface is indistinguishable from the corresponding kinematic (Pe ) ∞) interface. The situation changes drastically at the end of the second period, as can be observed by comparing Figure 3 parts b and II. Diffusion swiftly merged and erased many of the fine scale structures, causing an enormous reduction of the overall interface length. The patterns corresponding to the third and fifth period (Figure 3 parts c and d) show the attainment of a persistent oscillatory behavior in the geometry of the structures. These oscillatory

patterns can be seen qualitatively as the resultant of a dynamical equilibrium between two competing mechanisms, namely the convection-driven generation of interface and the merging of neighboring patterns caused by molecular diffusion. Although these mechanisms have been hypothesized since the early studies focusing on the interaction between convection and diffusion, direct quantitative investigation of this phenomenon with respect to interface evolution in two-dimensional time-dependent flows has never been undertaken, to the best of our knowledge. An important question arises as to the role of chaos in determining the dynamics of reaction interfaces. To explore this issue, we analyzed the case of a stirring protocol that possesses large islands, as can be obtained, for example, by setting the periodicity of the sine flow system to the value T ) 0.4 (Figure 1a). Figure 4 parts I and II shows the snapshots of the kinematic (Pe ) ∞) interface at times nT, with n ) 20 and 30, respectively. In this case, the material interface undergoes altogether different stretching processes inside and outside the chaotic region, with the overall rate of growth being exponential in the chaotic region (X-shaped area in the figures), and linear within the islands. The corresponding reaction interfaces at Pe ) 104 (Figure 4 parts a and b) coincide with the kinematic template within the islands, whereas the fine structure of the diffusionless limit inside the chaotic region is evidently blurred into one large lamella at the times considered. Snapshots of the reactive patterns at later times (n ) 50, Figure 4c, and n ) 60, Figure 4d) again point to a persistent oscillatory evolution in the mixing patterns. An overall quantitative description of interface dynamics can be obtained by tracking the length of the reaction interface L(t) vs time. We consider a wide range of Pe ) 102 ÷ 105 for the two mixing protocols (T ) 0.4 and 1.6). Figure 5a shows the results for the nearly globally chaotic case T ) 1.6, at Pe ) 103, 104, and 105 (continuous lines), together with the growth of the length of the kinematic interface (line with points). For each of the Pe considered, we can unambiguously identify a crossover time t* ) t*(Pe) (and a corresponding breakup length L*) beyond which the dynamics of reaction interfaces departs irreversibly from exponential kinematic growth and settles into a bounded oscillating pattern around a characteristic average length depending on both Pe and the mixing protocol. The case T ) 0.4 (Figure 5b) displays more complex features. Here it is possible to identify two separate crossover times t/C and t/P corresponding to the breakup of the kinematic interface inside and outside the chaotic region, respectively. For t/C < t < t/P, the interface undergoes an oscillating restructuring driven by repeated merging events inside the chaotic region, whereas the overall trend of monotonic growth is driven by the fraction of interface that falls within the islands. If we target the Pe ) 104 case (second curve from below), the breakup time t/P is of the order of t/P ) 11.6, which corresponds to approximately n ) 28 periods, whereas t/C ) 2.1. Comparison of the kinematic and reaction interfaces at n ) 30 (t ) 12; Figure 3 parts b II) supports the observation that significant merging of structures within the islands begins only at times larger than t/P ) 11.6. We note that the high-frequency fluctuations (in the time-scale of a half-period) clearly detectable for T ) 0.4, and especially at high Peclet numbers, are not spurious consequences of numerics but rather derive from merging and restructuring events between parts of the interface that lie “at the border” between the quasiperiodic and chaotic regions.


Giona et al. breakup of the reaction interface, which coincides with the kinematic interface up to that time, occurs when the diffusional lengthscale ldiff ) (2t/Pe)1/2 is of the same order of magnitude as the average lamellar thickness 1/Lkin(t), i.e., ldiff(t*) ) RL-1kin(t*). In the last relationship, the prefactor R ∼ O(1) accounts for the heterogeneity in the spatial behavior of the short-time Lyapunov exponents, which determines the fine structure of the local striation thickness, and can be estimated from a single simulation experiment at low Peclet number, e.g., Pe ) 102. The length of the kinematic interface is given by Lkin(t) ) L0 exp(htopt/T), L0 ) 2, where htop is the topological entropy of the system estimated from its Poincare` map (htop ) 2.33 for T ) 1.6). t* is therefore implicitly expressed by the equation Pe ) 8R2t* exp(2θt*/T), with R ) 0.62 for T ) 1.6. The agreement of this simple model (solid line in Figure 8) with the simulations (triangular symbols) is satisfactory in view of the fact that interface breakup is an extremely complex process, in which the spatial heterogeneity of the local stretching field plays a significant role. All of these qualitative observations find a rigorous explanation in the analysis of eq 4 within the context of the theory of infinite dimensional dynamical systems32 and follow from spectral theory of linear (albeit not self-adjoint) operators in Hilbert spaces.33 Let us consider the case of the SF system defined on the unit square with periodic boundary conditions. Equation 4 is an evolution equation generated by the advection/ diffusion operator

L[φ;t] ) - v(x, t)‚∇φ + Figure 5. Length of the reaction interface L(t) vs time for the Sine Flow at T ) 1.6 (a) and 0.4 (b). Continuous lines: Pe ) 103, 104, and 105. Dotted line: kinematic interface.

Qualitatively analogous results are observed for the PCF (Figure 6). Again, at short times, kinematic (Figure 6I) and reaction (Figure 6a) interface are substantially identical, whereas at later times, the fine-scale structure of the kinematic interface (Figure 6 parts II and III) finds no correspondence in the reaction interface (Figure 6 parts b and c), because molecular diffusion cannot support the presence of such small length scales. Figure 7 depicts the behavior of the length of the reaction interface for the PCF flow at T ) 0.6 and Pe ) 103. As already observed for the SF system, the length of the reaction interface follows initially the kinematic template and settles asymptotically onto a regime of persistent oscillations. With our attention turned back to the SF system, it is worth noting that the breakup time t* does not correspond to the extinction of the reaction. Indeed, by targeting the T ) 1.6 case, one obtains mA(t*) ) 0.40, 0.56, and 0.66 for Pe ) 103, 104, and 105, respectively. This implies that for high Peclet numbers the major contribution to reactant consumption occurs in the mixed regime (t > t*) rather than in the kinematics-controlled regime (t < t*). Moreover, in the mixed regime, the partially mixed structures attain an almost constant average thickness approximately equal to 1/L*. Figure 8 shows the scaling of the breakup length L* vs the Peclet number for both T ) 0.4 and 1.6. In the first case, both L(t/C) and L(t/P) are shown. In all of the cases examined, L*(Pe) ∼ Peν over three decades, where ν ) 1/8 for L(t/C), ν ) 1/3 for L(t/P), and ν ) 0.4 for L* at T ) 1.6. With our attention focused on the nearly globally chaotic case (which is not further complicated by the occurrence of merging events in the chaotic and quasiperiodic regions), the behavior of L* vs Pe can be explained by means of elementary scaling arguments. The first

1 ∆φ Pe

(18)

in a subspace [this subspace is the Sobole¨v space H˙ 1per(I 2) of square integrable functions on I 2 ) [0, 1] × [0, 1] possessing zero mean square integrable first order partial derivatives (gradient)] of the functional space of square summable functions on the unit torus L2per(I 2), the norm of which is given by

||φ(t)|| ) [

∫I φ 2(x, t) dx]1/2, 2

φ ∈ L2per(I 2)

(19)

In the case of piecewise-steady time-periodic flows (such as the SF system), the operator L [φ;t] reduces to two distinct autonomous operators

Li ) - vi(x)‚∇ +

1 ∆, i ) 1, 2 Pe

(20)

where i ) 1, 2 corresponds to the first and second half period of motion. Associated with eq 4 in the presence of time-periodic velocity fields with instantaneous switching, a Poincare´ operator in the space of square summable functions P, L2per(I 2) f L2per(I 2), can be defined as

P [φ] ) expL2T/2°expL1T/2 φ

(21)

where ° indicates composition. The Poincare´ operator maps the difference function φ(x, t ) nT) into the function φ(x, t ) (n + 1)T) after one period of motion. The Poincare´ operator is linear, and therefore, its asymptotic properties depend on its eigenvalue-eigenfunction structure. A first property follows from the dissipative nature of the advection-diffusion equation eq 4. After elementary manipulations (by multiplying eq 4 with φ(x, t), and by integrating over



Figure 6. Comparison of the kinematic (Pe ) ∞) and reaction interfaces at Pe ) 103 for the prototypical cavity flow at T ) 0.6. (a-c) Mixing patterns (white and gray), and reaction interface (black line) at the end of the first, second, and third period, respectively. (I-III) Structure of the kinematic interface at the corresponding time instants.

I 2), it is straightforward to show that for any t > 0

d||φ||2 2 ) - ||∇φ||2 dt Pe where

||∇φ||2 )

() ∂φ

∑i ∫ ∂x

(22)

2

dx

(23)

i

Consequently, all of the eigenvalues {µi} associated with the Poincare´ operator possess a modulus strictly less than 1 for any Pe > 0.

The asymptotic regime is controlled by the fundamental mode (or modes) corresponding to the eigenfunction (or to the eigenfunctions) associated with the eigenvalue possessing lowest modulus, say µ1. Two cases may arise, depending on whether µ1 is real or complex. If the dominant eigenvalue µ1 is real, the stroboscopic sampling of the function φ(x, t) behaves as

φ(x, nT) = (-sign(µ1))na0e-nlog(1/|µ1|)ψ1(x)

(24)

where ψ1(x) is the dominant eigenfunction of the Poincare´ operator associated with µ1, a0 is a constant related to the


Figure 7. Length of the kinematic (line a) interface and of the reaction interface (line b) vs time for the PCF flow at T ) 0.6 and Pe ) 103.

Figure 8. Scaling of breakup length L*(t/C) vs Pe for the SF. (4), T ) 1.6; (°), T ) 0.4. The continuous line represents the prediction of L* ) L*(Pe) derived from scaling arguments (see main text) for the case T ) 1.6. The boxes (0) represent L*(t/P) at T ) 0.4.

projection of the initial condition onto the subspace spanned by ψ1(x), and sign(‚) is the sign function (sign(µ1) ) 1 if µ1 > 0 and sign(µ1) ) - 1 if µ1 < 0). In this case, the reaction interface sampled at t ) nT coincides with the zero level set of the eigenfunction ψ1(x). The stroboscopic sampling of the length of the reaction interface saturates toward a constant value (given by the overall length of the zero level set of the real eigenfunction ψ1(x)), and the decay of the L2 norm of the difference function is exponential. This is for example the case of the SF system at T ) 1.6, i.e., in the globally chaotic case. Conversely, if the dominant eigenvalue is complex, µ1 ) |µ1|eiω1, ω1 * 0,π (and consequently there exists also its complex conjugate |µ1|e-iω1), the dominant eigenspace is two-dimensional and is spanned by the real and imaginary parts of the complex eigenfunction ψ1(x) ) ψr1(x) + iψi1(x) associated with µ1. In this case, ψ(x, nT) scales asymptotically as

φ(x, nT) = (-sign(µ1))ne-nlog(1/|µ1|){[a0ψr1(x) - b0ψi1(x)] cos(nω1) - [a0ψi1(x) + b0ψr1(x)]sin(nω1)} (25) where a0 and b0 are the coefficients of the projection of the initial profile onto the subspace spanned by ψr1(x) and ψi1(x). The stroboscopic sampling of the length of the reaction interface displays periodic or quasiperiodic oscillations depending on whether ω1 is rational or irrational. An example of this behavior is given by SF protocols that possess islands of quasiperiodic motion, e.g., T ) 0.4 and 0.8.

Giona et al. To conclude the analysis of the geometry of reaction interface, let us briefly address the connection between finitely and infinitely fast reactions. The understanding of reaction interface dynamics for instantaneous reactions constitutes an important preliminary step to approach the finite-rate reaction case, which is further complicated by the presence of the nonlinear reaction term, such as - kCACB in the bimolecular case. In this case, the advection-diffusion-reaction equation is characterized by two dimensionless parameters, namely, the Peclet number Pe already introduced and a dimensionless rate constant κ ) kCrefL/ Vref ) (kCrefL2/D)/Pe. Intuitively, one expects that in cases where the reaction is not too slow (i.e., κ . 1) the structure of the reaction zone, identified through the product kCACB, must have a close connection with that of the infinitely fast case. As a quantitative confirmation of this observation, Figure 9 shows the comparison between the infinitely fast reaction interface (Figure 9a-c) at Pe ) 5 × 103, κ ) 10Pe in the SF system, and the contour plot of the normalized product CACB/maxx∈M(CACB). It can be noticed that the “hot spots” of the reaction, i.e., the spatial locations where the rate of product generation is the highest, are centered around the reaction interface corresponding to the limit of infinitely fast reaction, which acts as a sort of backbone around which the kinetic process is organized. This geometric observation is further supported by the analysis of overall reactant consumption (Figure 10). Because κ ) 10Pe corresponds to a relatively fast reaction, the overall reactant decay at short and intermediate time scales follows the trend characterizing the instantaneous case. Of course, at very large t, the decay of the finitely fast bimolecular case follows the scaling mA(t) ∼ t-1, whereas instantaneous reactions decay exponentially mA(t) ∼ exp(-λt). Nevertheless, at short/intermediate time scales, i.e., in the convection controlled regime, the geometric information related to the evolution of the reaction interface may prove useful as qualitative and quantitative tool for the interpretation of the dynamics of finitely fast chemical reactions. The scaling theory of finitely fast reactions has been considered in detail by Sokolov and Blumen.9 6. Asymptotic Scaling in Globally Chaotic Flows This section addresses the scaling properties of reactant consumption vs time in the asymptotic regime. The presence of a diffusive contribution, no matter how large Pe, can be viewed as a singular perturbation of the advection equation, which switches from a first order (hyperbolic) equation, the solution of which is controlled by its characteristic lines, to a second order (parabolic) equation, possessing dissipative behavior. This is a generic statement, valid for arbitrary velocity fields, yet it tells us very little about the coupling between diffusion and advection. In point of fact, the interplay between molecular diffusion and convection in a globally chaotic timeperiodic flow seems to yield a new singular phenomenon which has no counterpart for two-dimensional steady flows and, to the best of our knowledge, has not yet been reported in the literature. Specifically, we provide numerical evidence for a singular behavior, analogous to a phase transition, of the exponent controlling the asymptotic rate of reactant decay as a function of the Peclet number. We support this observation by means of scaling arguments. Indeed, from the functional setting outlined at the end of section 5, it follows that for square-integrable and bounded velocity fields, and for any value of 1/Pe > 0, reactant decay (or, equivalently, the L2 norm of φ) undergoes an asymptotic exponential decay:

||φ(t)|| ∼ mA(t) ∼ e-λ(Pe)t

(26)



Figure 9. Comparison between the reaction interface for an infinitely fast reaction, and the reaction zone at Pe ) 5 × 103 and κ ) 10Pe (SF system at T ) 1.6). (a-c): infinitely fast reaction interface at the end of the first- second- and third period. (d-f): contour plot of the product CACB/maxx∈M(CACB) at the same time instants.

where λ(Pe) ) log(1/|µ1|)/T is related to the absolute value (or to the norm, depending on whether µ1 is real or complex) of the dominant eigenvalue of the Poincare` operator, P, defined by eq 21. This result is in apparent contrast with the superexponential (exponential of an exponential, EE) scaling law determined by Tang and Boozer.12 The paradox arises as a direct

consequence of the fact that the authors consider a discrete model (the Arnold’s cat map, ref 34) in place of the continuous advection-diffusion equation. On the basis of this assumption, they derive a purely diffusive equation characterized by tensor diffusivity whose component along the stable direction increases exponentially in time, while being independent of the spatial


Giona et al.

Figure 10. Reactant concentration decay for a finitely fast reaction (κ ) 10Pe, Pe ) 5 × 103) for the SF system at T ) 1.6 (line a). Line b is the decay for the corresponding infinitely fast reaction under the same conditions.

Figure 11. Reactant concentration decay for the autonomous SF system v(x) ) v1(x). The arrow indicates increasing values of Pe ) 102, 103, and 104.

position. An analogous assumption is discussed in ref 9 in the context of one-dimensional approximate lamellar systems, leading to the same (EE) law. In the case 1/Pe ) 0, from eq 22, it follows ||φ(t)|| ) constant, and therefore

λ(∞) ) λ(Pe)|1/Pe)0 ) 0

(27)

We are interested in the scaling behavior of the decay exponent λ(Pe) vs Pe at high Pe, because the decay exponent can be viewed as the order parameter of the singular transition from pure convection to convection-diffusion phenomena. To this purpose, it is useful to compare the scaling behavior of λ(Pe) for several typical situations. In the case of pure diffusion on the torus, it results

λ(Pe) )

4π2 Pe

(28)

(The meaning of the Peclet number in the absence of convection might seem ambiguous because Vref is undefined. However, the torus topology allows us to bypass this ambiguity by considering a constant vector field of constant orientation and magnitude Vref. Independently of the orientation and Vref, the advectiondiffusion equation is transformed into a diffusion equation with constant diffusivity when written in the convected coordinates. Therefore, the physical meaning of eq 28 is that the decay exponent is proportional to the diffusion coefficient. We prefer to keep the Pe notation even for this case, as we are interested in comparing the diffusional scaling eq 28 with the dependence of λ(Pe) vs Pe in other regimes.) According to eq 28, the order parameter in pure diffusion converges smoothly to λ(∞) ) 0 with a Pe-1 behavior. Figure 11 shows the decay of mA(t) in a two-dimensional autonomous flow, over a broad range of Pe. The velocity field is given by v(x) ) v1(x), i.e., the velocity field in the first half period of the time-periodic SF system. We refer to this case as the autonomous SF system. The initial conditions are given by eq 17 as in section 5. The autonomous SF system possesses intrinsic symmetries, because the velocity field has solely a nonvanishing component directed along the x direction. For any initial distribution of the reactants that does not depend on the x coordinate, the advective contribution vanishes, and the balance equation reduces to a purely diffusive motion. For generic initial conditions violating this elementary symmetry of the system, an effective interplay between advection and

Figure 12. Reactant concentration decay for the SF system at T ) 1.6 for Pe ) 103, 2 × 103, 5 × 103, 104, 4 × 104.

diffusion occurs. The asymptotic slopes depicted in Figure 13 line b indicate that λ(Pe) converges to zero as

λ(Pe) ∼

A , Pe f ∞ xPe

(29)

where A is a positive constant. (The numerical results for the dominant eigenvalue obtained from the asymptotic exponential decay of reactant are confirmed by the direct estimate of the eigenvalue spectrum. For the Autonomous SF, the coefficient matrix reduces to a block structure, in which each block is tridiagonal, see eq 11. In point of fact, the use of exponential mass-loss scaling with time is a classical, robust, and widely applied method to estimate the dominant eigenvalue in linear transport theory.37) The exponent 1/2 in the scaling of the decay exponent with Pe is the signature of the convection-enhanced regime typical of two-dimensional autonomous flows. Equation 29 can be viewed as the consequence of the scaling properties of the effective diffusion coefficient, justified by Childress and Soward35 by means of boundary-layer arguments and proved mathematically by Fannjiang and Papanicolau36 through variational methods. Let us now consider the case of a globally chaotic flow, such as the SF system at T ) 1.6. Figure 12 shows the reactant decay in the range [102, 105]. As can be observed, the slope of these curves in a normal-log plot saturates toward a constant value as depicted in Figure 13 curve c. This means that

λ(Pe) = λ∞ ) constant * 0, Pe f ∞

(30)


J. Phys. Chem. A, Vol. 106, No. 23, 2002 5733 the quantity Peλ(Pe). Therefore, the ratio

||∇φ|| ||g|| ) = s-1(Pe) ||φ|| ||φ||

Figure 13. λ(Pe) vs Pe. Squares: SF at T ) 1.6. Circles: SF at T ) 0.8. The continuous lines represent the scaling laws λ(Pe) ) R/Pe (where R is a constant), and the scaling λ(Pe) ) 4π2/Pe that corresponds to a purely diffusive process.

(33)

becomes asymptotically a function of the Peclet number, and the quantity s(Pe) can be physically interpreted as the width of the boundary layer on which the concentration gradients are localized. Equation 33 can be viewed as a definition of s-1(Pe) in terms of λ(Pe), s-1(Pe) ) [Pe λ(Pe)]1/2. The issue is therefore reduced to obtaining an expression for the L2 norm of the gradient ||g||. This can be achieved by manipulating eq 4 (essentially by taking the gradient, performing the dot product of the resulting equation with ∇φ, and finally integrating over the space) so as to obtain the evolution equation for ||g||2:

d||g||2

) -2

dt

∫I

2

(∇v): g g dx -

2

||∇gi||2 ∑ Pe i

(34)

where ||∇gi|| is the norm of gradient of the ith component of the vector g and

-

∂Vi

∫I (∇v): g g dx ) - ∑∫I ∂x gigj dx 2

2

i,j

(35)

j

The latter term is the rate of stretching of the square norm of a vector field, which for a globally chaotic flow, is proportional to ||g||2

-

Figure 14. Time behavior of ||g(t)||2/||φ(t)||2 for the SF system at T ) 1.6 for several values of the Peclet number. (a) Pe ) 100, (b) 1000, (c) 10000. The straight lines parallel to the x axis are the corresponding values of the quantity Peλ(Pe).

The scaling exponent of the order parameter with Pe is therefore equal to 0, and eq 30 implies the occurrence of a discontinuity for Pe f ∞

lim λ(Pe) ) λ∞ * λ(∞) ) 0

Pef∞

(31)

This discontinuity can be viewed as a kind of a phase transition, which finds no correspondence in the case of two-dimensional autonomous velocity fields. The comparison of the scaling laws expressed by eqs 30 and 29 points to a definitely more pronounced interplay between advection and diffusion in the globally chaotic case than in a regular (non chaotic) flow. Intuitively, this can be explained as the balancing between competing mechanisms, namely, of the tendency of chaotic advection to increase exponentially the length of material interfaces and the action of diffusion that levels concentration gradients out. Next, we justify eq 30 through scaling arguments. Because d||φ(t)||2/dt ∼ -2λ(Pe)||φ(t)||2, asymptotically, from eq 22 it follows for large t

λ(Pe) =

1 ||g||2 Pe ||φ||2

(32)

where we defined g ) ∇φ. To provide an example confirming eq 32, Figure 14 shows the time behavior of the quantity ||g||2/ ||φ||2, for several values of the Peclet number, compared with

∫I (∇v): g g dx = K||g||2 2

(36)

where K is a positive constant that does not depend on Pe. The conjecture eq 36 stems from the properties of globally chaotic flows. Equation 36 implies that the rate of stretching of the norm of a vector because of the effects of advection

∆(t) ) -

1 ||g||

∫I (∇v): g g dx 2

(37)

is independent of the Peclet number. Indeed, in the diffusionless setting (1/Pe ) 0), the square norm of a vector field grows exponentially as a function of time, ||g(t)||2 = ||g(0)||2e2θt, where θ is the topological entropy of the flow.38,39 Consequently, ∆(t) is a function of time, the mean value of which

∆I(t) )

∫0t∆(τ) dτ

1 t

(38)

converges to the topological entropy θ for t f ∞. Thus, the estimate eq 36 implies that the quantity ∆(t), where the vector field g is just the gradient of φ, would behave as in the diffusionless limit for large Pe in the presence of a globally chaotic flow. To support this conjecture, Figure 15A shows the behavior (Within the numerical approach based on Gale¨rkin expansion, the integral quantities entering eq 37 associated with field gradients can be directly expressed in terms of the spectral coefficients φnm of the field φ, thus avoiding discretization in performing the derivatives.) of ∆(t) in the case of the SF at T ) 1.6 for two values of Pe spanning two decades, and Figure 15B shows the averaged quantity ∆I(t). As conjectured, ∆(t) is on average independent of Pe, and ∆I(t) converges for any Pe large enough toward a constant value independent of Pe, close to the topological entropy of the flow. Because d||g||2/dt = 2λ(Pe)||g||2, it results ||∇gi||2 = As2-(Pe)||gi||2 ∼ As2-(Pe)||g||2/


Giona et al. The nonchaotic case gives rise to a typical convection-enhanced diffusion characterized by a modified dependence (ζ ) 1/2) of the dominant eigenvalue on the Peclet number, when compared to a purely diffusive case. For a globally chaotic flow, the competition between advection and diffusion provides a further improvement (which indeed is the essence of mixing): the decay exponent λ(Pe) is asymptotically independent of Pe. From eqs 32 and 33 this further implies that

s(Pe) ∼ Pe-1/2

(42)

i.e., the width of the boundary layer scales as the reciprocal of the square root of the Peclet number. A similar scaling result was obtained by Klapper40 by using shadowing techniques on Wiener trajectories simulating the advection-diffusion of a scalar. The asymptotic properties of flow systems giving rise to the occurrence of regions of chaotic and quasiperiodic motion (e.g., the SF system at T ) 0.4 or 0.8) deserve particular attention and will be treated in full detail elsewhere. 7. Prediction of Overall Reactant Decay

Figure 15. (A) Time behavior of ∆(t) defined by eq 37 for the SF at T ) 1.6 for two values of the Peclet number Pe ) 102 (solid line) and 104 (dashed line). (B) time behavior of ∆I(t) defined by eq 38: Pe ) 102 (solid line) and 104 (dashed line). The dotted line represents the topological entropy θ ) 1.425 ( 0.03 for the SF at T ) 1.6.

2, where A is a positive constant independent of Pe. This scaling expression follows straightforwardly from eq 33 applied to φ ) gi, i ) 1, 2, and from the fact that ||g1|| ∼ ||g2|| ∼ ||g||, because the two components of the vector g scale with time t in the same way. By combining this result with eq 32, it follows

λ(Pe) ∼ -K +

As-2(Pe) 2Pe

(39)

From eq 33, λ(Pe) ) s-2(Pe)/Pe, and therefore

λ(Pe) )

2K ) constant A-2

(40)

which corresponds to the scaling observed in the numerical simulations. (To sum up, the scaling argument proposed is derived through simple algebraic steps starting from the exact eq 34 and on the basis of the scaling assumption eq 36 for the rate of growth of the norm of a generic vector in a globally chaotic flow. In the development of this scaling analysis, the hypothesis of global chaos flow is explicitly used in writing eq 36. For a flow possessing large islands of quasiperiodicity, eq 36 is no longer true, because the reaction interface (around which the highest concentration gradients are localized) may be located with the islands of quasiperiodic motion.) To sum up, the order parameter λ(Pe) in advection-diffusion phenomena scales with Pe as a power law λ(Pe) ) Pe-ζ for Pe f ∞, and the value of the scaling exponent ζ is indicative of the transport conditions, i.e., of the interplay between advection and diffusion:

{

1 pure diffusion ζ ) 1/2 two-dimensional autonomous flows 0 two-dimensional globally chaotic flows

}

(41)

This section provides an illustration of how the knowledge of reaction interface dynamics in chaotic flows can be used to predict in an approximate way the overall reaction rate and reactant decay of infinitely fast reactions. To this end, in line with the well-established lamellar approach, we model the interplay between advection and diffusion through a onedimensional system composed by two lamellae of the same thickness xj(t)/2: 2 ∂ψ(x, t) x dxj(t) ∂ψ(x, t) 1 ∂ ψ(x, t) + ) , ∂t ∂x Pe ∂x2 xj(t) dt

x ∈ (0,xj(t)) (43) equipped with periodic boundary conditions and initial condition ψ(x, 0) ) 2 - 4η(x - 1/2). The periodic boundary conditions are in this case straightforwardly inherited by the spatial periodicity of the flow domain. The linear velocity profile appearing in the convective term models the flow along the stable direction of a hyperbolic stagnation point within a onedimensional space (0,xj(t)) that is globally shrinking at a rate of dxj(t)/dt. The implementation of this approach rests upon the knowledge of the time-behavior of dxj(t)/dt. Explicit expressions for reactant decay were previously obtained by assuming a constant rate of shrinking R, i.e., d log xj(t)/dt ) R ) constant.9 With this assumption, in particular conditions of symmetry for the initial condition (e.g., with alternated lamellae all of the same thickness), it is possible to obtain a superexponential EE decay. Geometric information about reaction interface enters directly the definition of xj(t): we assume the simplest condition that is consistent with flow incompressibility, i.e.

xj(t) )

L0 L(t)

(44)

where L0 and L(t) are respectively the initial and the current overall length of the reaction interface. We observe that the geometric counterpart of a EE decay is represented by a monotonic unbounded exponential growth of the reaction interface, which is in contrast with what obtained in section 5. It is important to stress that this is a classical “engineering” model, and the aim of this section is just to show that the geometric information deriving from the knowlegde of the reaction interface length with time can be fruitfully used to



Figure 16. Overall reactant decay for the protocols T ) 1.6 (circles) and 0.4 (boxes). The continuous lines represent the predictions obtained through eq 45.

predict the order of magnitude of reactant decay. This approach is essentially valid for short/intermediate time scales and should not be confused with the asymptotic scaling analysis for the dominant eigenvalue developed in section 6. By enforcing the symmetries, the total amount of reactant A can be expressed as mA(t) ) ∫1/2 j(t) ξ, t) dξ, and the 0 ψ(x solution of eq 43 yields

4

mA(t) )

π

∞

∑ 2 n)1

1 - (-1)n n2

exp

(

-4π2n2 Pe

∫0t

dτ jx2(τ)

)

(45)

Figure 16 shows the comparison of eq 45 with simulations for Pe ) 104 at T ) 0.4 and 1.6. Analogous results are obtained for values of Pe in the range considered above. The level of agreement is satisfactory and, in some cases (e.g., T ) 0.4), extremely good. This indicates that the dynamics of the reaction interface provides the fundamental information required for the prediction of the evolution of infinitely fast reactions in chaotic flows, at least for short/intermediate time scales. This observation applies in a broader sense to finite rate reactions, where the reaction interface associated with an infinitely fast reaction provides the backbone around which the kinetic process is organized. Indeed, the extension of eqs 43 and 45 to finitely fast reactions is not particularly difficult. 8. Conclusions The main goal of this work was to extend the geometric approach to mixing from a purely kinematic context to systems subjected to molecular diffusion. As a natural physical framework for pursuing this extension, we considered a twodimensional stirred system in the presence of two segregated reactants undergoing infinitely fast reaction. The fact that the characteristic time of reaction is arbitrarily smaller (for instantaneous reactions the characteristic reaction time is zero) than those of convection and diffusion, implies that the species remain segregated at all times, thus allowing us to extend the notion of partially mixed structures and intermaterial contact area even in the presence of diffusion. The intermaterial contact surface coincides, in this case, with the reaction interface. The main goals of this paper are (i) to provide a qualitative understanding of the geometry of reaction interfaces, (ii) to provide a simple functional analytic background on the time behavior of the reaction interface supporting the numerical findings, and (iii) to show the occurrence of a new phenomenol-

ogy characterizing the interplay between diffusion and advection in globally chaotic flows, expressed by the scaling of the dominant eigenvalue controlling reactant decay as a function of the Peclet number. The evolution of the geometric structures generated by the stirring protocol are characterized by the occurrence of two regimes, namely, a kinematics-dominated growth, during which the reaction interface behaves like a material line passively advected by the flow, and a mixed regime resulting in persistent oscillations. The latter regime results from the intertwined action of stretching and folding of the mixing patterns because of chaotic advection and the merging of contiguous lamellae by the action of diffusion. To analyze the specific role of chaos in determining the fate of segregation patterns, we considered mixing protocols characterized by large, medium-sized, and negligible islands. In the case of large islands, the different time scales associated with convective mixing resulted into two growth regimes of the reaction interface. These phenomenological observations were made quantitative by analyzing eq 4 as a dynamical system evolving in an infinite dimensional functional space. In particular, the establishment of a persistent oscillatory regime for the geometry of partially mixed structures is a consequence of the dissipative nature of the evolution equation. By exploiting flow incompressibility, the information about the time behavior of the length of the reaction interface allowed us to derive a simple one-dimensional model that predicts the overall reaction rate and product generation over a wide range of Pe numbers with acceptable accuracy, at least for short and intermediate time scales. The asymptotic properties of the solution of the advectiondiffusion equation are particularly interesting. We provided numerical evidence of, and scaling arguments on, the occurrence of a scaling regime typical of globally chaotic flows, characterized by a constant value of the decay exponent independent of the Peclet number for high Pe. The convergence of the dominant eigenvalue λ(Pe) toward a constant value for high Pe should not be confused with the behavior of the conversion-time curves. The conversion-time curves do not collapse onto each other, but rather shif parallel to each other maintaining for large Pe the same asymptotic slope (as can be observed from Figure 12). Indeed. in the singular limit 1/Pe ) 0, no reaction occurs, and mA(t) ) mA(0) (see eq 31). In other words, the asymptotic behavior is reached after a transient that lasts longer and longer as Pe increases. At Pe ) ∞, the transient is infinitely long, and the conversion-time curve is just a constant corresponding to the initial reactant mass. The analysis of the asymptotic properties for flows which are not globally chaotic, and are characterized (in the diffusionless limit) by the simultaneous presence of invariant regions of quasiperiodic and chaotic motion, is more subtle. Chaotic and quasiperiodic regions behave as two different “phases”, with the mass transfer that occurs across the separatrices being controlled by diffusion. As a consequence, multiple crossovers in the scaling of the order parameter λ(Pe) as a function of Pe may occur. The phenomenology is further complicated by its significant dependence on the initial conditions. The specific treatment of this case will be developed in detail elsewhere. As a final remark, it is important to stress that the results obtained for two-dimensional flows can be directly extended to three-dimensional systems. The equation for the evolution of the norms of φ and g ) ∇φ, eqs 22 and 34, hold in any dimensions and so do the results deduced from these equations.

5736 J. Phys. Chem. A, Vol. 106, No. 23, 2002 When moving from two to three dimensions, the main difference is that three-dimensional autonomous velocity fields may give rise to chaotic behavior. References and Notes (1) Baldyga, J.; Bourne, J. R. Turbulent Mixing and Chemical Reactions; J. Wiley: New York, 1999. (2) Sokolov, I. M.; Blumen, A. Phys. ReV. A 1991, 43, 2714. (3) Chella, R.; Ottino, J. M. Chem. Eng. Sci. 1983, 39, 551. (4) Childress, S.; Gilbert, A. D. Stretch, Twist, Fold: The Fast Dynamo; Springer-Verlag: Berlin, 1995. (5) Aref, H. J. Fluid Mech. 1984, 143, 1. (6) Ottino, J. M. The kinematics of mixing: stretching, chaos and transport; Cambridge University Press: Cambridge, U.K., 1989. (7) Beigie, D.; Leonard, A.; Wiggins, S. Chaos Sol., & Fract. 1994, 4, 749. (8) Muzzio, F. J.; Ottino, J. M. Phys. ReV. Lett. 1989, 63, 47. (9) Sokolov, I. M.; Blumen, A. Int. J. Mod. Phys. B 1991, 20, 3127. (10) Muzzio, F. J.; Ottino, J. M. Phys. ReV. A 1989, 40, 7182. (11) Rom-Kedar, V.; Poje, A. C. Phys. Fluids 1999, 11, 2044. (12) Tang, X. Z.; Boozer, A. H. Physica D 1996, 95, 283. (13) Adrover, A.; Cerbelli, S.; Giona, M. J. Phys. Chem. A 2001, 105, 4908. (14) Reigada, R.; Sagues, F.; Sokolov, I. M.; Sancho, J. M.; Blumen, A. Phys. ReV. Lett. 1997, 78, 741. (15) Reigada, R.; Sagues, F.; Sokolov, I. M.; Sancho, J. M.; Blumen, A. J. Chem. Phys. 1996, 105, 10925. (16) Reigada, R.; Sagues, S.; Sokolov, I. M.; Sancho, J. M.; Blumen, A. J. Chem. Phys. 1997, 107, 843. (17) Muzzio F. J.; Liu, M. Chem. Eng. J. 1996, 64, 117. (18) Zalc, J. M.; Muzzio, F. J. Chem. Eng. Sci. 1999, 54, 1053. (19) Sawyers, D. R.; Sen, M.; Chang, H. C. Chem. Eng. J. 1996, 64, 129.

Giona et al. (20) Sawyers, D. R.; Sen, M.; Chang, H. C. Int. J. Heat Mass Transfer 1998, 41, 3559. (21) Mokrani, A.; Castellain, C.; Peerhossaini, H. Int. J. Heat Mass Transfer 1997, 40, 3089. (22) Raynal F.; Gence, J. N. Int. J. Heat Mass Transfer 1997, 40, 3267. (23) Giona, M.; Adrover, A.; Muzzio, F. J.; Cerbelli, S.; Alvarez, M. M. Physica D 1999, 132, 298. (24) Alvarez, M. M.; Muzzio, F. J.; Cerbelli, S.; Adrover, A.; Giona, M. Phys. ReV. Lett. 1998, 81, 3395. (25) Giona, M.; Adrover, A. Phys. ReV. Lett. 1998, 81, 3864. (26) D’Alessandro, D.; Daleh, M.; Mezic, I. IEEE Trans. Automat. Control 1999, 44, 1852. (27) Aref, H.; Balachandar, S. Phys. Fluid 1986, 29, 3515. (28) Muzzio, F. J.; Swanson, P. D.; Ottino, J. M. Phys. Fluids A 1991, 3, 822. (29) Sadhan, J. C.; Metcalfe, G.; Ottino, J. M. J. Fluid Mech. 1994, 269, 199. (30) Liu, M.; Muzzio, F. J.; Peskin, R. L. Chaos, Sol., & Fract. 1994, 4, 869. (31) Ferziger, J. H.; Peric, M. Computational Methods for Fluid Dynamics; Springer-Verlag: Berlin, 1996. (32) Temam, R. Infinite-dimensional dynamical systems in mechanics and physics; Springer-Verlag: Berlin, 1997. (33) Dunford, N.; Schwartz, J. T. Linear Operators; J. Wiley: New York, 1988. (34) Arnold, V. I.; Avez, A. Ergodic problems of classical mechanics; Benjamin: New York, 1968. (35) Childress, S.; Soward, M. J. Fluid Mech. 1989, 205, 99. (36) Fannjiang, A.; Papanicolau, G. SIAM J. Appl. Math. 2001, 62, 129. (37) Crank, J. The mathematics of diffusion; Clarendon Press: Oxford, U.K., 1975. (38) Eckmann, J.-P.; Ruelle, D. ReV. Mod. Phys. 1985, 57, 617. (39) Newhouse, S.; Pignataro T. J. Stat. Phys. 1993, 72, 1331. (40) Klapper, I. Phys. Fluids A 1992, 4, 861.

Quantitative Analysis of Mixing Structures in Chaotic Flows Generated

Recommend Documents