Eigenvalue–eigenfunction analysis of infinitely fast reactions and micromixing regimes in regular and chaotic bounded flows

Chemical Engineering Science 59 (2004) 2125 – 2144

www.elsevier.com/locate/ces

Eigenvalue–eigenfunction analysis of in!nitely fast reactions and micromixing regimes in regular and chaotic bounded %ows Stefano Cerbelli, Valerio Vitacolonna, Alessandra Adrover, Massimiliano Giona∗ Dipartimento di Ingegneria Chimica, Facolta di Ingegneria, Centro Interuniversitario sui Sistemi Disordinati, e sui Frattali nell’Ingegneria Chimica, Universita di Roma “La Sapienza”, via Eudossiana 18, 00184 Roma, Italy Received 28 June 2003; received in revised form 13 January 2004; accepted 11 February 2004

Abstract We analyze the dynamics of a single irreversible reaction A + B → Products, occurring in a bounded incompressible %ow. Within the limits of in!nitely fast kinetics, the system is reduced to an advection–di7usion equation for the scalar , representing the di7erence between the reactant concentrations. By the linearity of the governing PDE, the system evolution is determined by the properties of the eigenvalue–eigenfunction spectrum associated with the advection–di7usion operator. In particular, the dependence of the dominant eigenvalue —yielding the time-scale controlling the asymptotic reactant decay—as a function of the molecular di7usivity, D, for di7erent stirring protocols is analyzed. We !nd  ∼ D , where the exponent  ∈ [0; 1] depends upon the kinematic features of the stirring %ow. When the kinematics is regular within most of the %ow domain (e.g. two-dimensional autonomous %ows or time-periodic protocols possessing large quasiperiodic islands) a purely di7usive scaling,  = 1 settles as D → 0. The singular scaling  = 0 is found in the case of globally chaotic kinematics, whereas mixed regimes, 0 ¡  ¡ 1, occur in %ows that are characterized by the coexistence of quasiperiodic and chaotic behavior. The analysis of spectral properties of the advection–di7usion operator provides a new classi!cation of micromixing regimes, and new mixing indices for quantifying homogenization performances in the presence of di7usion. ? 2004 Elsevier Ltd. All rights reserved. Keywords: Fluid dynamics; Micromixing; Reacting %ows

1. Introduction and motivation Mixing-controlled reactions have been widely investigated and are of great importance for the optimization of industrial processes (Baldyga and Bourne, 1999). Contributions to this subject include theoretical, numerical and experimental works which have focused on a variety of physical systems, ranging from the simplest example of a single irreversible reaction A + B → Product (Sokolov and Blumen, 1991; Muzzio and Ottino, 1989), to more complex schemes of parallel competitive/consecutive reactions (see e.g. Zalc and Muzzio, 1999 and references therein). In the case of a single irreversible reaction occurring in a closed %ow, imperfect mixing re%ects into a slower rate of reaction (apparent rate) than that predicted on a chemical kinetics basis, whereas the equilibrium state is represented by a uniform product distribution. Moreover, in multiple reaction schemes, mixing can also in%uence the overall amount ∗ Corresponding author. Tel.: +39-06-445-85-892; fax: +39-06-445-85-339. E-mail address: [email protected] (M. Giona).

0009-2509/$ - see front matter ? 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2004.02.005

of products at equilibrium, even if all of the reactions are irreversible. Experimental evidence of mixing-sensitive yield is reported in e.g. Baldyga and Bourne (1990) and Nienow et al. (1992) for competitive-parallel reactions, and in e.g. Bourne and Kozicki (1977) and Garcia-Rosas and Petrozzi (1990) for competitive-consecutive reactions. Yield sensitivity is particularly appealing when an experimental validation of mixing models is sought, and has motivated continuous e7ort in mixing modeling that goes beyond the practical relevance of these types of systems (Baldyga and Bourne, 1999). Despite the manifold reasons that justify interest in multiple reactions, there are many fundamental questions concerning the single reaction scheme that are yet to be answered. The ultimate goal would be to “design” a stirring !eld (i.e. to specify proper geometric and operating conditions for the reactor and for the stirring equipment) so as to make the cooperative action of convection di7usion as effective as possible in order to maximize the overall reaction rate. Where the characteristic time of reaction is much smaller than those of convection and di7usion, the overall rate does

2126

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

not depend upon the chemical kinetics (Chella and Ottino, 1984). In this case, usually referred to as in3nitely fast reaction, the evolution of the system is described by a single linear advection–di7usion equation governing transport of the scalar !eld , representing the di7erence between the concentration of the two reacting species. This situation represents a sort of “ePciency limit”, meaning that the overall rate of reaction of a !nitely fast kinetics will always be lower than that predicted by using the in!nitely fast-reaction approximation. At the same time, it provides a useful paradigm to single out the impact of mixing on reaction dynamics, thus allowing us to compare the performance of di7erent stirring protocols on a quantitative basis, and to determine the fundamental qualitative features characterizing homogenization dynamics in complex %ows. Many studies focusing on in!nitely fast reactions have approached the problem within a one-dimensional framework by analyzing the pattern dynamics and apparent rate in the case where the two reactants are initially distributed in an alternate array of striations, also referred to as lamellae (see e.g. Sokolov and Blumen, 1991 and therein cited references). The results of these works were extended to include the simultaneous action of convection by making use of the so-called warped-time transformation (Ranz, 1979; Chella and Ottino, 1984; Sokolov and Blumen, 1991), which models the stretching and folding processes induced by the stirring %ow through a one-dimensional e7ective velocity, ve (x; t). The main limitation of this approach is essentially that it does not include the fundamental properties of advection associated with the folding of material interfaces within the mixing space. For this reason, more recent studies have considered reaction dynamics in a higher-dimensional setting, by performing direct numerical simulations of the corresponding advection–di7usion–reaction equations. Numerical simulations of advection–di7usion–reaction kinetics in cellular and two-dimensional laminar chaotic %ows have been performed by Reigada et al. (1996, 1997) and by Muzzio and Liu (1996). Reigada et al. focus mainly on the scaling of the reactant quantities with time in order to determine the controlling regime from characteristic exponents, whereas Muzzio and Liu analyze the e7ects of the %ow protocol on the overall product formation and on the spatial heterogeneity of the product distribution. In previous articles (Giona et al., 2002a, b), we analyzed the evolution of an in!nitely fast reaction occurring in a two-dimensional %ow in the presence of di7usion by undertaking a geometric approach. In particular, the length of the interface separating the segregated reactants (reaction interface) was tracked for a !xed initial distribution of reactants and for di7erent mixing protocols. The information about the interface length was exploited in the framework of the warped-time approach to derive a one-dimensional model for predicting the overall reaction decay. Stemming from the analogy between homogenization of an advecting di7using scalar and in!nitely fast reactions, work relevant to the subject also includes contributions

from the !elds of applied mathematics and physics. In particular, several authors analyzed the asymptotic regimes of an advecting–di7using scalar in steady unbounded %ows of either cellular (i.e. spatially periodic) (Rosenbluth et al., 1987; Fannjiang and Papanicolaou, 1994), or random types (Castiglione et al., 1999; Fannjiang and Papanicolaou, 1997). The asymptotic regime was quanti!ed by the e5ective di5usivity, De , e.g. characterizing the long-time behavior of the mean square displacement r 2 (t) of an ensemble of di7using particles advected by the stirring %ow. The aim of this article is to address thoroughly the phenomenological picture associated with the interplay of advection and di7usion in closed incompressible laminar %ows, providing a functional–theoretical justi!cation for it. For this reason, we considered di7erent stirring protocols which cover the typical situations occurring in Lagrangian kinematics. The evolution of in!nitely fast reactions is thus framed in a broader perspective, where attention is focused on how the asymptotic behavior of the reactant decay depends upon the interaction between advection and di7usion in the limit of very small molecular di7usivity. The main central di7erence with respect to previous work focusing on e7ective di7usivity is that the physical problem examined here naturally leads to consider the advection–diffusion equation as de!ned in a bounded domain. Thus, the dynamics of any physical quantity (for example, the decay rate or the species distribution within the mixing space) can be completely characterized through the eigenvalue–eigenfunction spectrum associated with the linear dissipative advection–di7usion operator (or with the stroboscopic operator obtained from the advection–di7usion equation in the case of time-periodic %ows). In particular, the eigenvalue spectrum de!nes the hierarchy of decay exponents (i.e. of time scales) that are to be expected for a !xed initial distribution of reactants. Given that the asymptotic decay rate for a generic initial distribution is expressed by the dominant eigenvalue of the spectrum, it is fundamental to determine how the dominant eigenvalue, say , depends upon the system di7usivity D, for di7erent stirring protocols. E7ective convective mixing is expected to result in an enhanced regime,  ∼ D , where  ¡ 1 ( = 1 corresponds to the same behavior that one would observe in a still medium, i.e. without convection). Throughout this article, we consider simple deterministic %ow protocols giving rise to complex and chaotic Lagrangian kinematics. This corresponds to the analysis of the e7ect of a mean velocity !eld on dispersion and homogenization. In the case of turbulent %ow conditions, for which the velocity %uctuations are described by means of stochastic models, the spectral results obtained in this article should be complemented with statistical mechanical approaches (Majda and Kramer, 1999; Baldyga and Bourne, 1999). The paper is organized as follows. Section 2 brie%y reviews the in!nitely fast reaction approximation and the class of %ow systems analyzed throughout the article. The

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

characterization of di7erent %ow protocols is based upon their kinematic performance, i.e. the structure of invariant sets associated with the advection equation x˙ = v(x; t). Section 3 discusses the functional properties stemming from the linear dissipative structure of the advection–diffusion operator. Section 4 succinctly addresses the numerical issues in approaching the advection–di7usion equation, and the full eigenvalue spectrum of its Floquet operator in time-periodic velocity !elds. Section 5 analyzes the spectral (eigenvalue–eigenfunction) structure of the advection–di7usion operator for several stirring protocols. The asymptotic scaling regimes associated with di7erent stirring protocols, ranging from completely regular (autonomous) %ows, to a globally chaotic conditions, are discussed in detail. The dominant eigenvalue, yielding the asymptotic decay exponent, is computed via direct numerical simulation of the governing PDE through the scaling of the overall mass of reactants. The spatial properties of the eigenfunctions are also addressed. Finally, Section 6 frames the results obtained in Section 5 on the occurrence of di7erent mixing regimes within the more general discussion of micro- and macromixing in complex %ows. 2. Statement of the problem and ow systems 2.1. Reacting systems and mathematical formulation Throughout this article, we consider an elementary bimolecular reaction A + B → Product that occurs in a incompressible medium %owing with velocity v(x; t) in a closed region M. The stirring !eld v(x; t) is incompressible, i.e. ∇ · v = 0, and independent of the mixture composition (i.e. it is not modi!ed by the reaction extent). Furthermore, we make the following additional assumptions: (i) all of the physical properties of the mixture (density, viscosity, etc.) are independent of the composition; (ii) the reactants A and B are characterized by identical di7usivities within the mixture, DA = DB = D = constant. The last assumption is crucial in the mathematical treatment of in!nitely fast reactions, albeit in some practical cases may be oversimpli!ed. By setting  = CA − CB , (where CA , CB are dimensionless molar concentrations of the two reactants, rescaled by the same reference concentration) the material balance equation reduces to the linear advection–di7usion equation for the concentration di7erence , which in dimensionless form reads as @ 1 + v · ∇ = R; @t Pe

(1)

where the Peclet number, Pe, is the dimensionless group expressing the ratio of the characteristic time for di7usion to that of advection. Applications usually show large values of Peclet, in the range from Pe = 102 for laminar %ames to

2127

Pe = 1010 or higher in polymerizations (Muzzio and Liu, 1996). The di7usionless limit is hence given by Pe → ∞. Within the in!nitely fast reaction approximation, reactants that are initially segregated remain segregated at all times, and the evolution of the system is completely speci!ed once the solution (x; t) of Eq. (1) is known, since CA (x; t) = ((x; t) + |(x; t)|)=2, CB (x; t) = (−(x; t) + |(x; t)|)=2. The zero-level set of the -function, 0 (t) = {x | (x; t) = 0} possesses therefore a direct physical meaning, that is, it identi!es the reaction interface separating the species A and B at time t, see e.g. Giona et al. (2002b). Without loss of generality, we consider throughout this article stoichiometric loading conditions, which correspond for the reaction model considered to equal initial masses of the two reactants, and imply that the initial condition (x; 0) for (x; t) possesses zero mean. Under stochiometric loading, both reactants persist asymptotically, and the reaction interface is a non-empty set of points for any t ¿ 0. 2.2. Flow systems As a model of stirring %ow, we use a simple yet widely investigated one-parameter class of two-dimensional %ows, namely the sine %ow (SF) system (Franjione and Ottino, 1992; Liu et al., 1994), which provides a complete spectrum of qualitative kinematic behaviors observed in physically realizable (Ottino, 1989) as well as in industrial %ows (Zalc et al., 2001; Harvey and Rogers, 1996). The SF system is de!ned as the periodic sequence of two steady sinusoidal %ows, v1 (x) = (v1; x ; v1; y ) = (sin(2y); 0), and v2 (x)=(0; sin(2x)), acting alternately for a time T=2 on points of the two-dimensional torus T2 (or equivalently, on the points x=(x; y) of the unit square I2 =(0; 1)×(0; 1) with opposite edges identi!ed). In particular, we refer to such protocols as time-periodic sine %ow (TPSF), as opposed to the stationary %ow de!ned by the velocity !eld v1 (x) acting for an in!nite time autonomous sine %ow—ASF. This latter case can be regarded as a simple paradigm of a kinematically regular %ow, in which passively advected particles move along torus meridians (i.e. circumferences at constant y), each with its own revolution period,  = 1=|sin(2y)|. As far as time-periodic protocols are concerned, their kinematic analysis is best performed in a time-discrete framework, i.e. by considering the PoincarUe map, P(x), obtained by integrating the advection equation x˙ = v(x; t)

(2)

over the %ow period T . Therefore, the current positions, xn = x(t)|t=nT , of a non-di7using passively advected scalar, initially located at x0 , at multiple integers nT of the %ow period T , are given by the autonomous system xn = P(xn−1 ) = Pn (x0 ), where Pn denotes the nth-fold recursive application of the PoincarUe map P. Since in TPSF the stirring protocol is generated by periodically blinking each T=2 time units two steady %ows v1 (x)

2128

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

Fig. 1. PoincarUe sections of the TPSF system: (A) T = 0:4; (B) T = 0:56; (C) T = 0:8; (D) T = 1:18.

and v2 (x), the PoincarUe map P can be obtained explicitly. It results
 Px (x; y) P(x) = Py (x; y)
 x + (T=2) sin(2y) = mod: 1; y + (T=2) sin[2x + (T=2) sin(2y)] (3) where the “mod. 1” condition accounts for the spatial periodicity of the %ow domain. The fact that the %ow is divergence-free (∇ · v = 0) re%ects into the area-preserving property of P. Qualitatively di7erent behaviors of the trajectories associated with the PoincarUe map stem from its nonlinear dependence on the spatial coordinates. Speci!cally, di7erent kinematic conditions may occur, corresponding to point particle trajectories living (i) on a !nite set of points (periodic orbits), (ii) on a collection of closed one-dimensional loops, or (iii) densely !lling a two-dimensional subregion of the %ow domain (chaotic region). As far as kinematic mixing is concerned, the presence of closed invariant curves centered around periodic elliptic points are associated with low mixing ePciency, as regions enclosed by these loops— referred to as islands—do not exchange material with the

surrounding %uid. Therefore, optimal convective mixing requires the chaotic region, say Cc , to be as large as possible (i.e. coinciding with the entire %ow domain for an ideal globally chaotic protocol). Figs. 1(A)–(D) show the phase space structure of the PoincarUe section for four stirring protocols corresponding to T = 0:4, 0.56, 0.8 and 1.18, respectively. These four protocols show the occurrence of a central X-shaped chaotic region (Fig. 1(A)), the structure and the width of which change with the %ow period T . As T increases, the extension of the chaotic region becomes broader, and tends to invade the whole mixing space. This occurs, e.g. at T = 1:6 (not shown), a case in which quasiperiodic islands are no longer detectable. This situation is typically referred to as “global chaos”, though the wording “nearly globally chaotic condition” would be more appropriate. Values of T in the range [0:4; 1:6] give rise to di7erent con!gurations of the quasiperiodic islands. For T = 0:4, two main quasiperiodic islands are present in the mixing space. The protocol corresponding to T = 0:56 is characterized by two shrunk main islands surrounded by an archipelago of newborn islets. At T = 0:8, the two main regular islands break up into four separate smaller islands. At T = 1:16, the size of these four islands is signi!cantly smaller than the size of the main chaotic region. The values of the Lyapunov exponent per unit time Lyap (obtained by dividing the Lyapunov exponent

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144 Table 1 Lyapunov exponent per unit time Lyap vs. the %ow period T for the TPSF system Period T

Lyapunov exponent Lyap

0.4 0.56 0.8 1.18 1.6 2.0

0.435 0.584 0.758 1.165 1.205 1.170

of the PoincarUe map P by the %ow period T ) associated with the trajectories belonging to the main chaotic region are reviewed in Table 1.

where ∇ 2L2 is the square integral norm of the gradient ∇. Eq. (6) states the dissipative nature of the advection –di7usion equation for any Pe ¿ 0, since it indicates that the L2 -norm of any solution of Eq. (4) with zero mean is strictly decreasing with time. This property of  2L2 stems from the fact that any function , di7erent from the zero function, possesses a non-vanishing ∇ L2 , since constant functions do not belong to H . Moreover, Eq. (6) states that the monotonically decreasing evolution of the L2 -norm of the solutions of the advection–di7usion equation is controlled by the L2 -norm of its gradient. For this reason, it is useful to express the evolution equations for ∇ 2L2 . By setting g = ∇, it follows: d g 2L2 = −2 dt

 I2

(∇v)T : gg dx

2  2  − ∇gi · ∇gi dx; Pe I2

3. Advection–di usion equation: spectral properties

2129

(7)

i=1

As discussed in the previous section, the hypothesis of in!nitely fast reaction allows to recast the system of advection–di7usion–reaction into a pure advection–di7usion problem Eq. (1) with no loss of information. Owing to the linear structure of the latter equation a formal analysis can be pursued, by framing Eq. (1) within the theory of in!nite-dimensional dynamical systems (Temam, 1997; Robinson, 2001). This section skips the functional analytical subtleties, and focuses on the spectral (eigenvalue/eigenfunction) properties which admit direct practical implications. In the dynamical system framework, the advection–di7usion equation can be envisioned as an evolution equation @ = L[]; @t

(4)

where the di7erential operator L, L[; t] = −v(x; t) · ∇ +

1 R Pe

(5)

acts on elements of a suitable functional space H . (In the case of the SF, the space H is the space of square summable and periodic functions on the unit square I2 , possessing zero mean. This space, usually indicated as L˙2per (I2 ), is a Hilbert space, possessing the norm f L2 , where f 2L2 =  |f(x)|2 dx:) Existence and uniqueness of the solutions I2 of Eq. (4) can be found elsewhere (Temam, 1997; Sell and You, 2000). Several important properties of Eq. (4) can be derived from the analysis of the evolution of the norm of the solutions. Let us consider a real-valued solution . By multiplying Eq. (4) by  and integrating over I2 , it follows after some algebra: d  2L2 2 = − ∇ 2L2 ; dt Pe

(6)

where the diadic term at the rhs of Eq. (7) is given by (∇v)T : gg =

2  @vi gi gj : @xj

(8)

i; j=1

Eq. (7) indicates that the evolution of the gradient norm is a consequence of two con%icting actions: the tendency to increase it due to the stretching action of the velocity !eld (!rst term at the rhs of Eq. (7)), and the smoothening e7ect of di7usion (second term at the rhs of Eq. (7)). Eq. (4) is a linear functional equation in a Hilbert space, and its fundamental properties stem from the spectral (eigenvalue) structure of the evolution operator L. In the analysis of the spectral properties, a distinction should be made between the autonomous and the time-periodic case. Let us !rstly consider the autonomous case, that is, the case where the velocity !eld v(x) does not depend explicitly on time. Under fairly general conditions on v(x) (e.g. it is suPcient to assume continuity), the spectrum of the advection–di7usion operator L possesses the following properties: (i) the spectrum reduces to a point spectrum formed by the eigenvalues , de!ned by the equation − v(x) · ∇ (x) +

1 R (x) =  (x); Pe

(9)

(ii) the eigenvalue spectrum is formed by a countable set of eigenvalues {n }∞ n=1 , such that the sequence 1=n admits zero as the unique accumulation point; (iii) the system of generalized eigenfunctions { n (x)}∞ n=1 , associated with the eigenvalues n forms a complete basis for the Hilbert space H . Proof of these results can be found in e.g. Hess (1986) and Faierman (1995).

2130

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

Before commenting on these results and their practical applications, let us formulate the spectral problem for time-periodic velocity !elds. In the time-periodic case, due to the explicit dependence of L upon the time variable t, it is meaningless to search for the eigenvalues of L. Nevertheless, by enforcing the time-periodicity of the stirring !eld, this problem can be overcome by de!ning a Floquet operator (also referred to as PoincarUe operator) F associated with the evolution of the advection–di7usion equation sampled at periodic instants of time equal to the period of the velocity !eld. Therefore, F maps the di7erence function (x; t = nT ), into the function (x; t = (n + 1)T ) after one period of motion. It can be shown that, if the velocity !eld is a bounded function of the position, i.e. there exists a positive constant Vm such that |v(x; t)| 6 Vm , for any x and t, the Floquet operator F is compact (Giona et al., 2004). This property admits a wealth of important implications (Reed and Simon, 1980). For one thing, it implies that eigenvalues $n of F, de!ned by F[ n (x)] = $n

n (x)

(10)

are either in !nite or countably in!nite number, and that the spectrum admits at the most zero as the unique accumulation point. Moreover, the compactness of F implies that F represents the limit operator of a sequence of !nite-rank operators Fn that can be obtained by projecting F into a nested sequence of !nite-dimensional subspaces. The convergence of Fn towards F is strong, in the operator norm (Reed and Simon, 1980). For any n, the reduced Floquet operator admits complete rank, meaning that the generalized eigenfunctions of Fn span the !nite-dimensional subspace (Giona et al., 2004). In the case of piecewise-steady time-periodic %ows (such as the TPSF system), the operator L[; t] is reduced to two distinct autonomous operators 1 Li = −vi (x) · ∇ + R; i = 1; 2; (11) Pe which act during the !rst and second half period of motion, respectively. For this class of %ows, the Floquet operator takes on the formal expression L2 T=2

F[] = e

◦e

L1 T=2

;

(12)

where “◦” indicates composition. It should be noted that when the %ow is autonomous, the Floquet operator can be de!ned by considering an arbitrary period T = ¿ 0, and the completeness of the system of eigenfunctions is inherited by the corresponding property of L. Like the advection di7usion operator, the Floquet operator is linear, and therefore all its properties depend upon its eigenvalue–eigenfunction spectrum. Let us denote the (countable) set of eigenvalues and eigenfunctions with ∞ {$k }∞ k=1 and { k }k=1 , respectively. Given that the advection–di7usion operator is not self-adjoint, the eigenvalues might be complex numbers. One important property that follows directly from the presence of the di7usive term R=Pe is that the modulus

|$n | of any eigenvalue is strictly less than unity for any !nite Pe, regardless of the structure of the convective !eld v(x; t). The eigenvalue $0 = 1 that corresponds to the uniform eigenfunction (x; y) = constant is not considered in the present discussion as the stoichiometric loading condi tion I2 (x; 0) dx = 0 implies that the di7erence function  possesses zero projection upon this eigenspace. The strict inequality |$k | ¡ 1 can be derived as a straightforward consequence of Eq. (6), which expresses the dissipative nature of the advection–di7usion equation. Indeed, let us assume that an eigenfunction (x) with |$| ¿ 1 exists, and let us take (x; t =0)= (x). Given that (x) is an eigenfunction, it would result (x; nT ) L2 = |$|n (x) L2 , which would imply that the norm of  diverges to in!nity, in contradiction with Eq. (6). Since the spectrum of F admits zero as the only possible limit point, its eigenvalues can be ordered in a non-increasing sequence such that 1 ¿ |$1 | ¿ |$2 | ¿ · · · ¿ |$k | ¿ · · ·. The largest eigenvalue $1 is the referred to as the dominant one. Since the system of generalized eigenfunctions forms a basis for H , any initial condition (x; 0)can be expanded ∞ in terms of eigenfunctions as (x; 0) = k=1 k (0) k (x). For the sake of notational simplicity, let us assume that the algebraic multiplicity of any eigenvalue equals its geometric multiplicity (i.e. the operator F attains a diagonal form when represented in the eigenfunction basis). This assumption can be easily removed without altering the main qualitative results derived below. It follows that the evolution of , at multiple integers of T reads as (x; nT ) =

∞ 

k (0)$kn

k (x):

(13)

k=1

For generic initial conditions, such that 1 (0) = 0, the asymptotic dynamics is governed by the !rst (dominant) eigenvalue $1 . Two cases may arise, depending on whether $1 is real or complex. If $1 is real, the scalar !eld (x; t) at large times is given by (x; nT ) ∼ 1 (0)$1n

1 (x):

(14)

As a consequence, the stroboscopic sampling of the normalized function norm (x; nT ) =

(x; nT ) $1n

(15)

asymptotically converges towards a time-invariant function, which is the dominant eigenfunction. √ Where $1 =|$|1 e2i!1 , i= −1, is complex (and therefore R I 1 (x) = 1 (x; y) + i 1 (x) is also complex), we have (x; nT ) ∼ 2|$1 |n {[A1 − [A1

I 1 (x)

R 1 (x)

+ B1

− B1

R 1 (x)] cos(2!1 n)

R 1 (x)] sin(2!1 n)}:

(16)

The asymptotic behavior of the normalized function norm (x; nT ) = (x; nT )=(2|$1 |n ) depends on !1 . If !1

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

is rational, say !1 = p=q with p and q integers, then norm (x; y; nT ) is a periodic function of n of period q as n → ∞. In the case where !1 is irrational, quasiperiodic oscillations arise, and the normalized di7erence function norm (x; nT ) will never converge towards a strictly time-invariant function. From Eqs. (14) and (16), it follows that the exponential asymptotic decay of a generic initial condition in the time-continuous frame is given by =−

log(|$1 |) : T

(17)

Thus, the !rst and foremost issue in quantifying the interplay between advection and di7usion in closed systems is to determine how the dominant exponent  and the dominant eigenfunction(s) associated with it depend on the velocity !eld v(x; t) and on the Pe number. More speci!cally, it is natural to ask whether stirring protocols possessing di7erent kinematic structures lead to qualitatively di7erent dependencies of the dominant exponent as a function of Pe. This is thoroughly addressed in Section 5, by enforcing several numerical approaches which are brie%y described in the next section. Henceforth, we will use the wording “eigenvalue” or “exponent” as synonymous to indicate the quantity  de!ned by Eq. (17) starting from the eigenvalue $1 of the Floquet operator. Observe that the quantity  in Eq. (17) is always positive.

4. Numerical approaches The numerical methods applied in this article are oriented towards two main goals, namely: (i) an ePcient and accurate direct numerical simulation of the advection–di7usion equation, and (ii) the estimate of the full eigenvalue spectrum. The !rst problem can be tackled by means of the Galerkin expansion. This can be conveniently implemented by expanding the di7erence function (x; t) with respect to the complex exponential basis {+r; s (x; y) = e2i(rx+sy) }∞ r; s=−∞ ,  (x; t) = r; s (t)e2i(rx+sy) ; (18) r; s

which forms the eigenbasis of the Laplacian operator. For a thorough analysis of the Galerkin method and its advantages in solving the advection–di7usion and the advection– di7usion–reaction kinetics in the presence of chaotic %ows see (Adrover et al., 2002). In the case of the TPSF system, the evolution equation for the Fourier coePcients r; s (t) reads as dr; s (t) 42 (r 2 + s2 ) =− r; s (t) dt Pe −r(r; s−1 (t) − r; s+1 (t))

(19)

2131

for nT=2 6 t 6 (2n + 1)T=2, n = 0; 1; : : :, and 42 (r 2 + s2 ) dr; s (t) =− r; s (t) dt Pe − s(r−1; s (t) − r+1; s (t))

(20)

for (2n + 1)T=2 6 t 6 (n + 1)T . The dissipative nature of the advection–di7usion operator ensures that a !nite number of Fourier modes, r; s ∈ [ − N; N ], yields an accurate approximate solution. Thus, by truncating the expansion to the !rst (2N + 1) modes for each coordinate axis, the evolution equation for the coePcients can be rearranged in the form d
(t) = A(t) ·
(t); (21) dt where
is a coePcient vector of (2N +1)2 entries, and A(t) is a piecewise-constant time-periodic (2N + 1)2 × (2N + 1)2 sparse matrix. In fact, the number of ODE to be solved is (2N + 1)2 − 1, since 0; 0 = 0 identically, due to the zero-mean condition. The choice of N , i.e. of the proper Galerkin truncation, depends signi!cantly upon the Peclet value (see Giona et al., 2002a). We use up to N = 200 for simulating the advection–di7usion equation for high Peclet number (Pe = 105 ). The computation of the eigenvalues spectrum of the Floquet operator F is obtained by solving the eigenvalue problem for the [(2N + 1)2 − 1] × [(2N + 1)2 − 1] matrix F de!ned by F = eA2 T=2 eA1 T=2 ;

(22)

A1 and A2 being two constant real matrices corresponding to the action of the !rst and of the second velocity !eld vi (x), i=1; 2. The matrix F represents the !nite-rank Floquet operator Fn with respect to the complex exponential basis. In order to !nd the eigenvalues and the eigenvectors of the matrix F (i.e. of the truncation of the Floquet operator), a QZ algorithm is used. The QZ algorithm is a variant of the well-known QR algorithm, originally devised for solving the generalized eigenvalue problem Fw = Bw (Golub and Van Loan, 1983), where F and B are two square matrices. The QZ algorithm can therefore be used to solve the classical eigenvalue problem Fw = w, when B is set equal to the identity matrix. The choice of the algorithm was restricted by the fact that the matrices dealt with are not symmetric. Furthermore, the whole spectrum needed to be obtained and not just a few selected eigenvalues. The QZ algorithm, (for B = identity matrix) consists of the reduction of the matrix F by Jordan blocks. More speci!cally, we used the Fortran subroutine rgg.f, which belongs to the public-domain package eispack and is available at the URL site www.gams.nist.gov. Figs. 2 and 3 show the spectrum of the Floquet operator associated with the TPSF system for several values of the Peclet number for T = 0:4 (Fig. 2) and for T = 1:18 (Fig. 3). It can be observed (see e.g. Figs. 2(C) and (D)) that as the number of modes increases, the eigenvalue spectrum converges towards an invariant spectral pattern

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

1

1

0.5

0.5 Im(µ)

Im(µ)

2132

0 -0.5

-0.5

-1

-1 -1

-0.5

(A)

0 Re(µ)

0.5

1

-1

-0.5

0 Re(µ)

0.5

1

-1

-0.5

0 Re(µ)

0.5

1

(B)

1

1

0.5

0.5 Im(µ)

Im(µ)

0

0 -0.5

0 -0.5

-1

-1 -1

-0.5

(C)

0 Re(µ)

0.5

1 (D)

1

1

0.5

0.5 Im(µ)

Im(µ)

Fig. 2. Eigenvalue spectra of the Floquet operator associated with the TPSF system at T = 0:4: (A) Pe = 100, N = 20; (B) Pe = 500, N = 20; (C) Pe = 1000, N = 20; (D) Pe = 1000, N = 30.

0

0 -0.5

-0.5

-1

-1 -1

-0.5

0.5

1

1

1

0.5

0.5

0

-0.5

0 Re(µ)

0.5

1

-1

-0.5

0

0.5

1

0 -0.5

-0.5

-1

-1 -1 (C)

-1 (B)

Im(µ)

Im(µ)

(A)

0 Re(µ)

-0.5

0 Re(µ)

0.5

1 (D)

Re(µ)

Fig. 3. Eigenvalue spectra of the Floquet operator associated with the TPSF system at T = 1:18: (A) Pe = 500, N = 30; (B) Pe = 1000, N = 30; (C) Pe = 2000, N = 30; (D) Pe = 5000, N = 30.

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

(depicted as the point set in the Re($)–Im($) plane). Given that as N increases new eigenvalues must appear, “spectral invariance” means that the newborn eigenvalues accumulate towards the origin of the complex plane, which represents the only possible limit point of the spectrum. For an accurate estimate of the dominant eigenvalue and of the corresponding eigenfunction we use the classical power method (Golub and Van Loan, 1983), which applies solely when the dominant eigenvalue is real. The power method is a direct numerical approach, based on the iterative application of the recursion F
n
n+1 = |F
n |

|λn|

100

1

0.01

(23)

starting from a generic initial condition, until convergence is reached.

2133

C branch

0

500

S branch

1000

1500

2000

n Fig. 4. Spectral branches of the Floquet operator associated with the TPSF system at T = 0:56. n = −log(|$n |)=T vs. n (n represents a labelling for ordering the eigenvalues in a increasing sequence for each of the two branches). The arrows indicate increasing values of the Peclet number Pe = 100; 200; 500; 1000.

5. Spectral structure and scaling properties The linear nature of the advection–di7usion problem and the analysis developed in Section 4, indicates that the key issue in understanding the interplay between advection and di7usion lies in the scaling of the eigenvalue spectrum as a function of the Peclet number, with particular attention on the behavior of the dominant eigenvalue in the limit for Pe → ∞. This section thoroughly addresses this issue by considering the TPSF system for several values of the period T in order to cover, in a fairly complete way, the di7erent kinematic conditions, corresponding to non-chaotic Lagrangian advection, partially chaotic advection, and global Lagrangian chaos. The cases of autonomous (two-dimensional %ows) and globally chaotic conditions (e.g. the TPSF system at T =1:6) have been analyzed elsewhere (Giona et al., 2002a; Cerbelli et al., 2003), and therefore are succinctly reviewed. All the other cases, which correspond to kinematic regimes displaying the coexistence of a chaotic region and of a manifold of quasiperiodic islands of di7erent relative sizes are addressed in this article for the !rst time. Before speci!cally discussing the spectral properties, let us mention an important symmetry of the TPSF system. Eqs. (19) and (20) shows that the dynamics of the Fourier coePcients { r; s = r;Rs + i r;I s } is characterized by a real-valued coePcient matrix. This implies that the real { r;Rs } and the imaginary { r;I s } part of the Fourier coePcients are decoupled. More speci!cally, any initial pro!le (x; t = 0), for which the Fourier representation is purely real generates, as a solution of the advection–di7usion equation, a pro!le which possesses this property for any time t ¿ 0. A complementary behavior is observed for those initial pro!les which admit a zero real part of the Fourier coePcients, and a non-vanishing imaginary part. As a result of this decoupling between the real and the imaginary part of the Fourier expansion, two spectral

branches may be identi!ed, referred to as the C (for cosine) and S (for sine) branch, respectively. The C branch is formed by eigenvalues whose eigenfunctions possess vanishing imaginary part ( r;I s = 0, for all r; s). As a consequence, the eigenfunctions of this branch can be expressed as a cosine series (C)

(x) = 2

∞ 

R r; s

cos[2(rx + sy)]:

(24)

r; s¿0

Conversely, the S branch is characterized by the eigenfunction property r;Rs = 0 for all r; s, so that a generic eigenfunction of this branch reads as (S)

(x) = 2

∞ 

I r; s

sin[2(rx + sy)]:

(25)

r; s¿0

Fig. 4 shows the two spectral branches for the SF system at T = 0:56, for several values of the Peclet number. In this !gure, the eigenvalues of each branch are ordered in an increasing way. 5.1. The autonomous case The autonomous case is interesting for highlighting how di7erent scaling regimes coexist within one and the same spectrum. These regimes correspond to qualitatively di7erent interactions between advection and di7usion. The autonomous sine %ow is described by Eq. (19). The dynamics of the Fourier modes decouples into a countable family of independent subsystems parametrized with respect to the integer r. Two di7erent cases arise, depending on whether r = 0 or = 0. In the !rst case, Eq. (19) is reduced to a pure di7usion equation, and admits a family of eigenvalues given by {−42 s2 =Pe} s = 1; 2; : : :, corresponding to the eigenfunctions of the Laplacian operator e2isy . This family of eigenfunction/eigenvalues form the

2134

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

1

1

b

c

Λ

λn

b

a

a

0.01

0.1

0.0001

100

1000

10000

100000

Fig. 5. Dominant eigenvalues vs. the Peclet number for the di7usive branch (◦) and the CoED branch (•) for the autonomous sine %ow. The solid lines are the scalings  ∼ Pe−1 (curve a), and  ∼ Pe−1=2 (curve b).

di5usive branch, and are characterized by a dominant di5usive exponent di7 = 42 =Pe. In the case r = 0, the coePcient matrix of each of the subsystems attains a tridiagonal form, and an e7ective coupling between advection and di7usion occurs. Spectral results indicate that the eigenvalues are complex conjugate and proportional to the reciprocal of the square root of the Peclet number. This means that there exists an eigenvalue branch the eigenvalues of which scale with Pe in a qualitatively di7erent way with respect to those belonging to the di7usive branch, the dependence on Pe being expressed by 1 ∼ √ : (26) Pe This branch is referred to as the convection-enhanced di5usive branch (CoED) which was coined by Papanicolau and coworkers in their analysis of asymptotic properties of two-dimensional cellular %ows in the presence of di7usion in open systems (Fannjiang and Papanicolaou, 1994). Fig. 5 shows the scaling of the dominant eigenvalues of these two branches. Throughout this section n indicates the norm of the eigenvalues of the Floquet operator in the time-continuous frame, n = −log(|$n |)=T . 5.2. The globally chaotic case The globally chaotic case has been addressed in detail elsewhere (Cerbelli et al., 2003). In this subsection, the results are brie%y reviewed, and new and more complete data are presented. Fig. 6 shows the dominant eigenvalues of the two spectral branches C and S. For large Peclet numbers, the dominant eigenvalue which is real and belongs to the S branch, appears to become independent of Pe, i.e. for Pe → ∞:

1000

10000

100000

Pe

Pe

 ∼ constant

100

(27)

This new scaling, which is di7erent from CoED, has been named chaos-enhanced di5usion (ChED), and a

Fig. 6. Dominant eigenvalues for the two spectral branches C (◦) and S of the TPSF system at T =1:6 vs. the Peclet number. Empty (◦) and !lled (•) circles refer to the dominant eigenvalue of the S branch (curve a) estimated by means of the QZ algorithm and the power method, respectively. Triangles ( and curve c) refer to the !rst complex eigenvalue of the C branch, and boxes ( and curve b) to the !rst real eigenvalue of the C branch. Solid lines are interpolating curves, depicted for enhancing visualization.

scaling model justifying it has been proposed, starting from the two fundamental equations (6) and (7), see Cerbelli et al. (2003). It is interesting to observe that not only the dominant eigenvalue, but also the second eigenvalue, third and so forth, ful!ll the same scaling law expressed by Eq. (27). For this reason, ChED regime may be referred to as a homogeneous coupling between convection and di7usion (see further Section 6). It is worth mentioning that results in line with the relationship expressed by Eq. (27) have been obtained by other researches (see e.g. Sukhatme and Pierrehumbert, 2002), for simpli!ed pulsed-system models, which are obtained by splitting the interplay between advection and di7usion in two separate steps: a !rst step in which advection is active and no di7usion occurs; a second step without advection, in which di7usion smoothens the pro!les. However, the use of pulsed systems for interpreting the interplay between advection and di7usion in real physical system is questionable, and is still matter of discussion and debate. Besides these authors (Giona et al., 2002a), Eq. (27) was proposed in Toussaint et al. (1995, 2000) for the case of a three-dimensional autonomous %ow giving rise to global Lagrangian chaos. Speci!cally, the connection between the eigenvalue scaling and the properties of the energy spectrum at intermediate wavenumbers (referred to as the Batchelor scaling) is addressed in Toussaint et al. (2000). Spectral data depicted in Fig. 6 indicate the occurrence of a crossing between the C and the S branch. For Pe values less than about 2000, the dominant eigenvalue belongs to the C branch and is a complex conjugate, while for larger Pe the !rst eigenvalue of the S branch (which is real) becomes the slowest. Fig. 7(A) shows the contour plot of the normalized dominant eigenfunction for Pe = 105 , and Fig. 7(B) the contour plot of the norm of its gradient. It can be observed that the

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

2135

Fig. 7. Dominant eigenfunction for the TPSF system at T = 1:6 and Pe = 105 : (A) contour plot of the normalized eigenfunction; (B) contour plot of its gradient.

1

S a

0.01

}

C

0.0001

5.3. Time-periodic cases with imperfect mixing The case of velocity !elds giving rise to imperfect mixing, that is, the coexistence of chaotic regions and quasiperiodic islands, is a completely unexplored !eld as far as a spectral characterization of the advection di7usion equation is concerned. Next, we discuss the results of spectral analysis. Let us start with the case T = 0:4, which in the di7usionless setting yields large recirculation islands separated by a narrow X-shaped chaotic region (Fig. 1(A)). The dominant eigenvalues of the two (C- and S-) spectral branches are characterized by two di7erent scaling properties (Fig. 8). The dominant eigenvalue of the C-branch is real and scales proportionally to the inverse of the Peclet number (which corresponds to a di7usive regime). This corresponds to the fact that advection does not signi!cantly in%uence the dominant mode, which is localized within the quasiperiodic islands. Fig. 9 shows the spatial structure of the dominant eigenfunction at Pe = 105 (N = 150), and Fig. 10 its contour plot. The dominant eigenvalue of the S branch is complex conjugate, and it scales with Pe as a CoED mode (i.e. proportional to Pe−1=2 ). A slightly di7erent spectral structure occurs for Pe=0:56, which is characterized by two main quasiperiodic islands and a rich structure of smaller ones embedded into the main chaotic region. Fig. 11 shows the dominant eigenvalues of the C and S branch. In this case, the dominant eigenvalues

b

λn

spatial region characterized by high values of the gradient tends to invade the whole mixing region, and that the contour plot (or the reaction interface associated with the dominant eigenfunction) closely resembles the structure of the leaves of the unstable foliation (Giona and Adrover, 1998). The data depicted in Fig. 7 are particularly neat, owing to the fact that the number of modes used to obtain the pro!les is 160 000 (corresponding to N = 200 modes in each spatial direction).

100

1000

10000

100000

Pe Fig. 8. Dominant eigenvalues of the C and S branch of the spectrum vs. the Peclet number for the TPSF system at T = 0:4. Empty circles (◦) represent the eigenvalue of the C branch: the !rst, second and third eigenvalues (see the direction of the arrow) are depicted. Solid circles (•) indicate the dominant eigenvalue of the S branch.√Curves (a) and (b) represent the scalings n = 42 =Pe, and n = A= Pe (A = 14:9), respectively.

of the two branches scale in the same way as a function of the Peclet number:  ∼ Pe−

(28)

with a scaling exponent  = 0:865 unambiguously di7erent from 1. The dominant eigenvalue belonging to the C branch is real and the structure of the dominant eigenfunctions is depicted in Fig. 12(A)–(C) for several values of the Peclet number. The occurrence of a scaling exponent  between 12 and 1 is a clear indication of the strong in%uence of the regions of poor mixing within which the dominant eigenfunction is essentially localized, but also of the interaction of chaotic advection, which modi!es the scaling with respect to a purely di7usive regime, as observed, e.g. in the case T = 0:4. This phenomenon can be better appreciated by the analysis of the gradient of the dominant eigenfunction for high Peclet number (see Fig. 13), which is characterized by

2136

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

1 λn

2

b

0.1

ψ

S 1

0.01

0

0.001

a

100

1000

} 10000

C 100000

Pe -1

-2 0

0 0.25

0.25

y

0.5

0.5

0.75

0.75 1

x

1

Fig. 9. Dominant eigenfunction pro!le for the TPSF system at T = 0:4, Pe = 105 (N = 150).

Y

1

-1.7 -1.2 -0.7 -0.2 0.2 0.7 1.2 1.7

0.5

0 0

0.5 X

Fig. 11. Dominant eigenvalues of the C and S branch of the spectrum vs. the Peclet number for the TPSF system at T = 0:56 . Empty circles (◦) represent the dominant eigenvalues of the C branch. The second ( ) and third () eigenvalues are also depicted (see the direction of the arrow). Solid circles (•) indicate the dominant eigenvalue of the S branch. Curves (a) and (b) represent the scaling law n ∼ Pe− with  = 0:865.

1

Fig. 10. Contour plot of the dominant eigenfunction the TPSF system at T = 0:4, Pe = 105 (N = 150).

high-gradient %ame-like structures embedded and lapping the chaotic region. As a !nal remark resulting from the analysis of these two cases, it can be observed that the localization of the dominant eigenfunction for high Peclet numbers tends to mimic the structure of the PoincarUe sections very closely (see Fig. 1). This is particularly evident from the inspection of the dominant eigenfunctions at T = 0:56 for increasing values of the Peclet number. As Pe increase, the archipelago of islets surrounding the main islands that characterize the kinematic PoincarUe section emerges with an increasing level

of detail from the contour plot of the dominant eigenfunction (see e.g. Fig. 12(C) for Pe = 105 ). In order to highlight the complex role of the quasiperiodic islands in the interplay between di7usion and advection, and the role of symmetry, the case T = 0:8 is particularly interesting. Fig. 14 shows the scaling of the dominant eigenfunctions of the C and S branches: both eigenvalues are real, yet their scaling is di7erent. The smallest C eigenvalue, which is dominant in the Peclet range [102 ; 106 ], follows Eq. (28) with  = 0:745, whereas the dominant S-eigenvalue scales di7usively, i.e. it is proportional to 1=Pe. The functional decoupling between C and S branches, which is typical of the SF system, indicates, that for very large Peclet values (approximatively of order 2 × 106 ), a crossing occurs between these two branches, and the S-branch di7usive eigenvalue will eventually dominate the asymptotic scaling. The physical picture underlying this phenomenon can be gained by the analysis of the eigenfunctions associated with these two branches. Figs. 15 and 16 show the contour plots (Fig. 15) of the dominant eigenfunctions of the two branches, and the contour plots of their gradient norm (Fig. 15), respectively. Essentially, these two di7erent scaling are associated with two di7erent mixing conditions corresponding to di7erent geometric localization properties of the eigenfunctions within the quasiperiodic islands. In the case of the C branch, the two reacting species (each of which corresponds to positive and negative values of the eigenfunction, respectively) are localized within the two di7erent super-islands separated by the main chaotic region (see Fig. 15(A)). (In order to avoid ambiguity, the four main egg-like islands appearing in the PoincarUe section of the SF system at T = 0:8 can be regarded as being organized in two main super-islands, since each couple of nearby islands is “connected” to each other via a chain of smaller ones, see Fig. 1(C).) As a consequence, mass transfer is mediated by the action of advection within the chaotic region, and this transfer mechanism can be appreciated by the

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

2137

Fig. 12. Dominant eigenfunction of the TPSF system at T = 0:56: (A) Pe = 103 (N = 60); (B) Pe = 104 (N = 110); (C) Pe = 105 (N = 150).

1

0.1

λn

S C

0.01

0.001

100

1000

10000

100000

Pe

Fig. 13. Contour plot of the gradient norm of the normalized dominant eigenfunction of the TPSF system at T = 0:56 and Pe = 105 (N = 150).

Fig. 14. Dominant eigenvalues for the two spectral branches (C and S) of the TPSF system at T = 0:8 vs. the Peclet number. The solid lines represent the scalings n = A Pe− (C branch) with  = 0:745, and n = B=Pe (S branch), respectively. A and B are two positive constants.

2138

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

Fig. 15. Contour plots of the dominant eigenfunctions for the TPSF system at T = 0:8 for Pe = 4 × 104 : (A) C branch (N = 140); (B) S branch (N = 120).

Fig. 16. Contour plots of gradient modulus of the dominant eigenfunctions for the TPSF system at T = 0:8 for Pe = 4 × 104 : (A) C branch (N = 140); (B) S branch (N = 120).

inspection of the gradient norm (Fig. 16(A)), the contour plot of which shows the occurrence of typical %ame-like structures embedded within the chaotic region. This situation can be described as an inter-island transfer (i.e. mass-transfer between di7erent islands) mediated by chaotic advection. This results in a dominant exponent that possesses an intermediate scaling law between purely di7usive regime and ChED. The structure of the dominant eigenfunction of the S branch displays an altogether di7erent transfer mechanism. Both reactants are present within each super-island, as can be observed by the distribution of white and black spots in the contour plot (Fig. 15(B)). It results an intra-island transfer (i.e. a mass-transfer within the same island), completely localized within each super-island. In this case, the role of advection acting in the main chaotic region is immaterial. This phenomenon can be clearly observed by the inspection of the contour plot of the gradient norm of the S-branch-dominant eigenfunction (Fig. 16(B)), which shows that the gradient

is completely localized within the quasiperiodic islands, with negligible smearing out e7ects in the chaotic region. To conclude this phenomenological overview of the many regimes that characterize the interplay between advection and di7usion in the presence of non-global Lagrangian chaos, let us consider the case of a kinematic condition giving rise to small quasiperiodic islands, as for the TPSF system at T = 1:18 (Fig. 1(D)). In this case, four neatly separated quasiperiodic islands are present. The dominant eigenvalues of the two branches are depicted in Fig. 17, and are both real. For high Peclet numbers Pe ¿ 1000, the two dominant C and S eigenvalues attain practically the same value. Fig. 17 also shows the second eigenvalue (complex conjugate) of the S branch which follows the same scaling Eq. (28) of the dominant one, with  = 0:55. In this case, mass transfer, is essentially controlled by advection within the chaotic region, as can be observed from the contour plots of the two dominant eigenfunctions, see Figs. 18(A) and (B). The main di7erence between the eigenfunctions of

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

λn

1

0.1

0.01

100

1000

10000

100000

Pe Fig. 17. Dominant eigenvalues for the two spectral branches C (◦) and S (•) of the TPSF system at T = 1:18 vs. the Peclet number. Boxes ( ) refer to the second (complex) eigenvalue of the S branch. The solid lines represent the scaling n = A Pe− (C branch) with  = 0:55 and A a positive constant.

the two branches is related to symmetry as can be observed by their contour plots (Fig. 18). For the dominant C-branch eigenfunctions, each reactant is localized within two of the four islands, whereas the S-branch dominant eigenfunction (which is doubly degenerate) involves the interaction of only two opposite islands, the remaining couple of islands being devoid of reactants. The spectral properties described above, show that the presence of quasiperiodic islands induces a scaling of the dominant eigenvalue following Eq. (28) with an exponent  attaining values between 0.5 and 1, and one would be tempted to argue that values smaller than 0.5 may not occur (the exponent  must take on values in the interval [0; 1]). In fact, this conclusion is incorrect, as the example that follows clearly shows. Let us consider the TPSF system for T = 2:0. This is a rather peculiar situation. At !rst glance, its PoincarUe section suggests the occurrence of globally chaotic conditions. In fact, this is not true, since can be proved analytically that the PoincarUe map of the kinematics

2139

Eq. (3) possesses two elliptic !xed points, located at ( 14 ; 14 ) and ( 34 ; 34 ). A zoom-in of the PoincarUe section near the elliptic !xed point ( 14 ; 14 ) is depicted in Fig. 19. The dominant eigenvalues of the C and S branches for T =2:0 are depicted in Fig. 20, and indicate that the dominant eigenvalue is real, that it belongs to the S branch and that it satis!es Eq. (28) with =0:366. The dominant eigenfunction is depicted in Fig. 21 for Pe = 105 , and displays a strong localization around the elliptic !xed points. In conclusion, this section has described the di7erent modes of interaction between advection and di7usion that may occur in the presence of non-global chaos. In all cases, localization around the region of poor kinematic mixing, is clearly evident. At large Pe values the dominant eigenvalue is always real. It belongs to the C-branch for small values of the protocol period (T ¡ 1:18), whereas it belongs to the S-branch for larger values of T . A comprehensive analysis of the spectral phenomenology described in this section reveals that very di7erent interaction regimes between advection and di7usion can settle, depending on the symmetry and the relative spatial con!guration of quasiperiodic regions within the mixing space. A result of this manifold of regimes is re%ected in the scaling law Eq. (28), characterizing the behavior of the dominant eigenvalue vs. the Peclet number, the scaling exponent of which may attain a value ranging from 0 to 1, as summarized in Table 2. A physical discussion of these results is organically tackled in Section 6, in the light of the classical concept of micromixing.

6. Asymptotic regimes and micromixing It is important to frame the results obtained in Section 5 within the general theory of mixing in chemical processes, by making connections with the concept of micromixing. Micromixing (Bourne, 1983) involves the interaction between stirring and di7usion at the lengthscales at which

Fig. 18. Contour plots of the dominant eigenfunction of the C and S branches of the TPSF system at T =1:18 for Pe=4×104 : (A) C branch; (B) S branch.

2140

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144 Table 2 Scaling exponent  of the dominant eigenvalue as a function of Pe (for large Pe values ¿1000) vs. the period T for the TPSF system. As it regards the case T = 0:8, it is reasonable to extrapolate a di7usive dominant scaling for Pe ¿ 2 × 106 . This is motivated by the fact that the two scalings observed for the C and the S branches are associated with two di7erent transport mechanisms, and there is no physical reason that these two branches would modify their qualitative behavior for Pe → ∞

Fig. 19. Zoom-in of the PoincarUe section of the TPSF Eq. (3) at T = 2:0, close to an elliptic !xed point.

1

λn

C S

0.1

0.01

100

1000

10000

100000

Pe Fig. 20. Dominant eigenvalues of the C and S branch of the spectrum vs. the Peclet number for the TPSF system at T = 2 . Circles (◦; •) represent the eigenvalue of the S branch, squares ( ) those associated with the C branch. The curve interpolating the eigenvalues of the S branch results n = A=Pe with  = 0:366.

Period T



Observations

0.4 0.56 0.8 ” 1.18 1.6 2.0 ∞ (Aut.)

1 0.865 0.745 1 0.55 0 0.37 1

— — For Pe ¡ 106 For Pe ¿ 2 × 106 (extrapolated) — — — CoED may occur with  = 12

molecular di7usion becomes e7ective. Its study can be approached experimentally, for instance, by means of mixing-sensitive reactions (see Baldyga and Bourne, 1999; Villermaux, 1991, for a review). If micromixing is identi!ed as the coupling between the action of the mean %ow and di7usion at di7usional lengthscales, the asymptotic behavior of the advection–di7usion equation (which may represent either the outcome of a dispersion experiment or the result of an in!nitely fast reaction) is the time-domain counterpart of mixing phenomena occurring at these length-scales. This simple observation makes it possible to resolve the homogenization processes at the smallest lengthscales from time-domain analysis of simple dispersion experiments, that is, ultimately, from the spectral analysis of the advection–di7usion equation.

Fig. 21. Dominant eigenfunction of the TPSF system at T = 2:0 for Pe = 105 (N = 150): (A) contour plot of the eigenfunction; (B) contour plot of the gradient modulus of the normalized eigenfunction.

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

1 e f

Λ

In order to avoid confusion, it is important to stress that there is a fundamental conceptual di7erence between “micromixing” as a process and “micromixing models”. Micromixing, as a physical process, refers to the interaction between advection and di7usion at the lengthscales of di7usional in%uence (Danckwertz, 1957; Bourne, 1983), while micromixing models are theoretical attempts “to describe the evolution of concentration %uctuations in a spatially homogeneous reactor by following Lagrangian particles as they move through the system” (Fox, 1998). It follows that the analysis of the asymptotic properties of the response curves is by no means a micromixing model, in the meaning envisaged by Fox (1998), but rather a direct consequence of mixing, expressed mathematically through the advection– di7usion equation in higher (two- and three-) dimensional spaces. All the scaling regimes described in Section 5 which consider the dominant eigenvalue either of the advection di7usion operator (for autonomous %ows) or of the associated Floquet operator (for time-periodic %ows) should be fully ascribed to micromixing (Bourne, 1983; Fox, 1998). It is remarkable that this phenomenon involving long-term behavior (and consequently short interaction lengthscales) is in%uenced by the mean structure of the stirring %ow, and speci!cally (i) by its stretching and folding action and (ii) by the geometric structure of the PoincarUe section in the mixing space expressed by the occurrence of chaotic regions and quasiperiodic islands. In fact, the results obtained in Section 5 can be schematically summarized in a compact form, in the de!nition of di7erent classes of mixing regimes, in a way that may serve as a guideline towards a theory of micromixing for deterministic %ows. Of course, for turbulent %ows, the seemingly stochastic nature of the stirring !elds may induce other mixing regimes than those considered so far. Since the mixing regimes considered in this paper are de!ned on the basis of generic asymptotic properties (which stem from the scaling of the dominant eigenvalue Eq. (28)), they actually correspond to what is customarily referred to as micromixing. Table 2 and Fig. 22 review the results obtained in Section 5, referred to the exponent . From this table it is possible to extract a classi!cation of mixing regimes, which can be subdivided into homogeneous mixing regimes and mixed-mode regimes (Table 3). The !rst class is characterized by the occurrence of the scaling law Eq. (28) homogeneously throughout the eigenvalue spectrum. This phenomenon admits a geometric counterpart, namely that the eigenfunction are global functions, the maxima and minima of which are distributed along the whole mixing space. This spectral and spatial homogeneity is shared by the two extreme cases of interaction between advection and di7usion: purely di7usive mixing (i.e. the absence of stirring) and %ow protocols giving rise to globally chaotic kinematics. The two cases give rise to the well-de!ned values of the exponent , which

2141

d c b

0.01

a

0.0001

100

1000

10000

100000

Pe Fig. 22. Dominant eigenvalue as a function of the Peclet number for the TPSF system. Curve (a) refers to T = 0:4, (b) to T = 0:56; (c) to T = 0:8; (d) to T = 1:16; (e) to T = 1:6; (f) to T = 2:0.

Table 3 Schematic classi!cation of the micromixing regimes in advection–di7usion dynamics Regime type

Stirring conditions

Scaling exponent 

Homogeneous regimes

Pure di7usion

1

Globally chaotic kinematics

0

Mixed-mode regimes

CoED—autonomous 2-d %ows

1 2

Partially chaotic %ows

(0,1)

represent the extremes in the range of variability of , namely  = 1 in pure di7usion and  = 0 in globally chaotic advection (ChED). The latter situation may occur either in two- and three-dimensional time-periodic %ows, and even in three-dimensional autonomous stirring protocols. The other broad class of micromixing regimes may be referred to as mixed-mode regimes. This class encompasses autonomous 2-d stirring protocols, and time-periodic protocols giving rise to partially chaotic Lagrangian kinematics. The main phenomenological properties of this class of mixing regimes are the following (i) occurrence of di7erent eigenvalue branches possessing di7erent scaling laws Eq. (28), i.e. di7erent exponents ; (ii) values of  in the range (0,1); (iii) localization of the eigenfunction around the poor mixing regions (quasiperiodic islands). Therefore, both spectral and spatial heterogeneities characterize this class. In fact, the interaction between di7usion and advection can be viewed ultimately as a consequence of two competing mechanisms: the relaxation of concentration gradients due to di7usion, and the stretching of material interfaces upon which these concentration gradients are localized. The occurrence of a scaling behavior Eq. (28), with  = 1, i.e. which di7ers from purely di7usive homogenization, indicates that these two mechanisms cannot be decoupled, not even at short lengthscales. In the case of homogeneous

2142

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

mixing regimes, theory and models have been proposed for explaining and predicting the value attained by the exponent : this is rather trivial for pure di7usion, but slightly more complicated for global chaos (see e.g. Giona et al., 2002a). Essentially, these scaling theories are grounded on the application of Eqs. (6) and (7) for the time evolution of the norms of  and of its gradient, by enforcing scaling arguments which are reasonable and valid due to the spatial homogeneity of the concentration patterns (in the limit of Pe → ∞). The case of mixed-mode regimes is much more complicated, due to the heterogeneity of the interplay between advection and di7usion, which implies that the occurrence of both regions of quasiperiodic and chaotic motion within the mixing space must be accounted for. At present, (with the exception of the case of two-dimensional autonomous %ows) there are no theories or models for predicting the value of the exponent , on the basis of the pure kinematics in bounded %ows. We believe that the discovery of such theories or models is central to the theory of micromixing, regarded as the asymptotic interaction between stirring and di7usion, and moreover that the spectral analysis in deterministic %ows may provide useful hints and models for tackling the asymptotic properties of dispersion even in turbulent %ows. 7. Concluding remarks The focus of this article was to analyze the in%uence of a convective %ow upon the asymptotic properties of the advection di7usion equations. This linear, and apparently simple equation describes both homogenization dynamics (as can be experimentally measured via a tracer dispersion experiment) and the evolution of in!nitely fast reactions. Since the advection–di7usion equation is linear, all its properties can be derived from its spectral structure (i.e. from the scaling of the eigenvalues and from the spatial properties of the eigenfunctions). As a model of stirring %ow, we considered a oneparameter class of piecewise-steady time-periodic sinusoidal %ows de!ned on the two-dimensional torus (i.e. on the unit-square domain equipped with periodic boundary conditions), but the spectral results obtained are generic for more complex and realistic %ow systems (including the case of stirred tanks). We have developed a thorough investigation of a variety of di7erent stirring conditions, which has led to the identi!cation of the main qualitative properties of the interaction between advection and di7usion. In short, these properties are summarized by the scaling of the dominant eigenvalue with the Peclet number (for high Pe ¿ 1000), and by the analysis of the micromixing regimes developed in Section 6. The fundamental problem which has yet to be solved, is a theory which explains the spectral results in mixed-mode regimes, (i.e. in stirring !elds giving rise to the simultaneous presence of chaotic and quasiperiodic regions), and

thus predicts the value of the exponent  from the structure of the stirring !eld. This is certainly, a fairly diPcult and challenging goal, which is similar to the corresponding problem faced in other scienti!c !elds: e.g. quantum theory of an electron moving in a periodic potential, spectral theory of the Mathieu equations, etc. (Karpeshina, 1997). Although these two examples refer to eigenvalue problems associated with self-adjoint operators (while the advection– di7usion operator is intrinsically not self-adjoint), a sound application of perturbation theory and variational methods may serve the purpose. As a !nal comment, let us brie%y address the application of the results obtained in classical industrial applications, e.g. involving dispersion in stirred tanks. The qualitative results obtained could provide a useful paradigm for dispersion and reaction in complicated %ow patterns, such as those associated with industrial equipment. The geometric complexity, intrinsic to the structure of a stirred vessel (induced by the presence of baXes, and by the impeller geometry), severely limits the application of the Galerkin expansion to this case. Nevertheless, eigenvalue/eigenvector analysis may be still applied to discretized formulations of the advection–di7usion equation. In particular, eigenvalue –eigenfunction analysis can be performed starting from a spatial discretization of the advection di7usion equation based, e.g. on a !nite-volume formulation. Work is under way in the direction of applying this formal apparatus to such applications, by further considering physically realizable %ow systems, such as the driven cavity %ow and the %ow between concentric/eccentric cylinders which are amenable to feasible Galerkin expansion of the advection– di7usion operator, and to spectral analysis.

Notation Abbreviations ASF CoED ChED TPSF

autonomous sine %ow convection-enhanced di7usion chaos-enhanced di7usion time-periodic sine %ow

Latin symbols C F Fn N Pe v

molar concentration of the chemical species  matrix representation of the Floquet operator F with respect to the Fourier basis matrix representation of the truncated Floquet operator Fn truncation order of the Fourier expansion for each coordinate axis Peclet number, expressing the ratio between the characteristic time of di7usion to that of convection velocity !eld

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

x P(x) t T

position vector of a generic point of the unit square PoincarUe map expressing the new position of a massless particle at x after one %ow period time temporal %ow periodicity

Calligraphic symbols D F Fn I2 L

molecular di7usivity Floquet operator expressing the scalar evolution at multiple integers of the %ow period truncated Floquet operator unit square interval representing a global projection chart for the two-dimensional torus Advection–di7usion operator, L = −v · ∇ + (1=Pe)R

Greek letters  R  n $n  r; s n

exponent expressing the power-law scaling of the dominant eigenvalue as a function of Pe Laplacian operator dominant eigenvalue, see Eq. (17) generic eigenvalue of L generic eigenvalue of F di7erence function,  = CA − CB coePcients of the expansion of (x; t) in complex double Fourier series generic eigenfunction associated with L or with F

References Adrover, A., Cerbelli, S., Giona, M., 2002. A spectral approach to reaction –di7usion kinetics in chaotic %ows. Comput. Chem. Eng. 26, 125–139. Baldyga, J., Bourne, J.R., 1990. The e7ect of micromixing on parallel reactions. Chem. Eng. Sci. 45, 907–916. Baldyga, J., Bourne, J.R., 1999. Turbulent Mixing and Chemical Reaction. Wiley, Chichester. Bourne, J.R., 1983. Mixing on the molecular scale (micromixing). Chem. Eng. Sci. 38, 5–8. Bourne, J.R., Kozicki, F., 1977. Mixing e7ects during the bromination of 1,3,5-trimetoxybenzene. Chem. Eng. Sci. 32, 1538–1539. Castiglione, P., Mazzino, A., Muratore-Gianneschi, P., Vulpiani, A., 1999. On strong anomalous di7usion. Physica D 134, 75–93. Cerbelli, S., Adrover, A., Giona, M., 2003. Enhanced di7usion regimes in bounded chaotic %ows. Phys. Lett. A 312, 355–362. Chella, R., Ottino, J.M., 1984. Conversion and selectivity modi!cations due to mixing in unpremixed reactors. Chem. Eng. Sci. 39, 551–567. Danckwertz, P.V., 1957. Measurement of molecular homogeneity in a mixture. Chem. Eng. Sci. 7, 116–117. Faierman, M., 1995. On the spectral theory of an elliptic boundary value problem involving an inde!nite weight. In: Gohberg, I., Langer, H. (Eds.), Operator Theory and Boundary Eigenvalue Problems. BirkhYauser, Basel. Fannjiang, A., Papanicolaou, G., 1994. Convection-enhanced di7usion for periodic %ows. SIAM J. Appl. Math. 54, 333–408. Fannjiang, A., Papanicolaou, G., 1997. Convection-enhanced di7usion for random %ows. J. Stat. Phys. 88, 1033–1076.

2143

Fox, R.O., 1998. On the relationship between Lagrangian micromixing models and computational %uid dynamics. Chem. Eng. Proc. 37, 521–535. Franjione, J.G., Ottino, J.M., 1992. Symmetry concepts for geometric analysis of mixing %ows. Philos. Trans. R. Soc. London A 338, 301–323. Garcia-Rosas, J., Petrozzi, S., 1990. In%uence of mixing in the azo coupling of 1-naphtol and diatosized aniline. Chimia 44, 366–368. Giona, M., Adrover, A., 1998. Nonuniform stationary measure of the invariant unstable foliation in Hamiltonian and %uid mixing systems. Phys. Rev. Lett. 81, 3864–3867. Giona, M., Cerbelli, S., Adrover, A., 2002a. Quantitative analysis of mixing structures in chaotic %ows generated by in!nitely fast reaction in the presence of di7usion. J. Phys. Chem. 106, 5722–5736. Giona, M., Cerbelli, S., Adrover, A., 2002b. Geometry of reaction interfaces in chaotic %ows. Phys. Rev. Lett. 88, 024501. Giona, M., Cerbelli, S., Adrover, A., 2004. Homogenization exponents in regular and chaotic bounded %ows, in preparation. Golub, G.H., Van Loan, C.F., 1983. Matrix Computations. Johns Hopkins University Press, Baltimore, MD. Harvey III, A.D., Rogers, S.E., 1996. Steady and unsteady computation of impeller-stirred reactor. A.I.Ch.E. J. 42, 2701–2712. Hess, P., 1986. On the asymptotic distribution of eigenvalues of some non selfadjoint problems. Bull. London Math. Soc. 18, 181–184. Karpeshina, Y.E., 1997. Perturbation Theory for the SchrYodinger Operator with a Periodic Potential. Springer, Berlin. Liu, M., Muzzio, F.J., Peskin, R.L., 1994. Quanti!cation of mixing in aperiodic %ows. Chaos Solitons Fractals 4, 869–893. Majda, A.J., Kramer, P.R., 1999. Simpli!ed model for turbulent di7usion: theory, numerical modelling and physical phenomena. Phys. Rep. 314, 237–574. Muzzio, F.J., Liu, M., 1996. Chemical reactions in chaotic %ows. Chem. Eng. J. 64, 117–127. Muzzio, F.J., Ottino, J.M., 1989. Evolution of a lamellar system with di7usion and reaction: a scaling approach. Phys. Rev. Lett. 63, 47–50. Nienow, A.W., Drain, S.M., Boyes, A.P., Mann, R., El-Hamouz, A.M., Carpenter, K.J., 1992. A new pair of reactions to characterize imperfect macro-mixing and partial segregation in a stirred semi-batch reactor. Chem. Eng. Sci. 47, 2825–2830. Ottino, J.M., 1989. The Kinematics of Mixing. Stretching, Chaos, and Transport. Cambridge University Press, Cambridge. Ranz, W.E., 1979. Application of a stretch model to mixing, di7usion and reaction in laminar and turbulent %ows. A.I.Ch.E. J. 25, 41–47. Reed, M., Simon, B., 1980. Methods of Modern Mathematical Physics: I—Functional Analysis. Academic Press, New York. Reigada, R., Sagues, F., Sokolov, I.M., Sancho, J.M., Blumen, A., 1996. Kinetics of the A + B → O under steady and turbulent %ow. J. Chem. Phys. 105, 10925–10933. Reigada, R., Sagues, F., Sokolov, I.M., Sancho, J.M., Blumen, A., 1997. Fluctuation dominated kinetics under stirring. Phys. Rev. Lett. 78, 741–744. Robinson, J.C., 2001. In!nite-dimensional Dynamical Systems. Cambridge University Press, Cambridge. Rosenbluth, M.N., Berk, H.L., Doxas, I., Horton, W., 1987. E7ective di7usion in laminar convective %ows. Phys. Fluids 30, 2636–2647. Sell, G.R., You, Y., 2000. Dynamics of Evolutionary Equations. Springer, New York. Sokolov, I.M., Blumen, A., 1991. Mixing in reaction di7usion problems. Int. J. Mod. Phys. 20, 3127. Sukhatme, J., Pierrehumbert, R.T., 2002. Decay of passive scalars under the action of single scale smooth velocity !elds in bounded domains: from non self-similar probability distribution functions to self-similar eigenmodes. Phys. Rev. E 66, 056302. Temam, R., 1997. In!nite-dimensional Mechanical Systems in Mechanics and Physics. Springer, New York.

2144

S. Cerbelli et al. / Chemical Engineering Science 59 (2004) 2125 – 2144

Toussaint, V., Carriere, P., Raynal, F., 1995. A numerical Eulerian approach to mixing by chaotic advection. Phys. Fluids 7, 2587–2600. Toussaint, V., Carriere, P., Scott, J., Gence, J.-N., 2000. Spectral decay of a passive scalar in chaotic mixing. Phys. Fluids 12, 2834–2844. Villermaux, J., 1991. Mixing e7ects on complex chemical reactions in a stirred reactor. Rev. Chem. Eng. 7, 51–108.

Zalc, J.M., Muzzio, F.J., 1999. Parallel-competitive reactions in a two-dimensional chaotic %ow. Chem. Eng. Sci. 54, 1053–1069. Zalc, J.M., Alvarez, M.M., Muzzio, F.J., Arick, B.E., 2001. Extensive validation of computed laminar %ow in a stirred tank with three Rushton turbines. A.I.Ch.E. J. 47, 2144–2154.