Foundations of laminar chaotic mixing and spectral theory of linear operators

Chemical Engineering Science 61 (2006) 2754 – 2761 www.elsevier.com/locate/ces

Foundations of laminar chaotic mixing and spectral theory of linear operators S. Cerbellia , A. Adrovera , F. Cretab , M. Gionaa,∗ a Dipartimento di Ingegneria Chimica, Universitá di Roma “La Sapienza”, via Eudossiana 18, 00184 Roma, Italy b Dipartimento di Meccanica e Aeronautica, Università di Roma “La Sapienza”, via Eudossiana 18, 00184 Roma, Italy

Received 10 June 2005; received in revised form 3 October 2005; accepted 4 October 2005 Available online 4 January 2006

Abstract This article addresses the potentiality and the range of application of the spectral theory of mixing as a simple and objective way to define in a rigorous way mixing performance and stirring efficiency. Attention is mainly oriented to the class of laminar flows which is the typical condition occurring in micromixers and flow microdevices. The application of the theory is highlighted by considering the cavity flow as a case study. 䉷 2005 Elsevier Ltd. All rights reserved. Keywords: Mixing; Advection–diffusion equations; Laminar flow; Microdevices; Transport processes; Diffusion

1. Introduction Mixing theory has attracted the interest of fluid dynamicists and engineers for its relevance in the understanding some of the fundamental problems involving fluid flows, and for its practical impact in connection with chemical, pharmaceutical and other manufacturing industrial processes (Baldyga and Bourne, 1999). The first quantitative approach to mixing dates back to Danckwerts (1952, 1958), who introduced global indices of the “degree of mixedness” such as the intensity of segregation, which essentially corresponds to the concentration variance of an advecting–diffusing dye carried by the flow, and the linear scale of segregation, yielding the average “diameter” of the partially mixed structures. Owing to its conceptual simplicity and experimental feasibility, this approach is still largely used to quantify mixing, especially when complex flows (e.g. flows in industrial equipment) are to be dealt with (Müller et al., 2004). About two decades ago, a new approach to characterize mixing in flow systems spread out in the fluid dynamics community (Aref, 1984; Ottino, 1989). This approach stemmed from the observation that even simple, large-scale velocity fields can generate mixing structures at arbitrarily small scale. This idea constituted a major breakthrough for people involved in fluid mixing. In fact, until then, the wordings “mixing flow” and

∗ Corresponding author. Tel.: +11390644585609; fax: +11390644585451.

E-mail address: [email protected] (M. Giona). 0009-2509/$ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2005.10.062

“turbulent flow” were used almost interchangeably. Methods and quantities derived from dynamical system theory were directly applied to simple flows, laboratory equipment, stirred tanks, static mixers, and constituted what is now referred to as the Lagrangian, or kinematic theory of laminar mixing (Ottino, 1989; Beigie et al., 1994). The kinematic approach has contributed to a better understanding of existing mixing devices as well as to the rational design of new ones (Harvey and Rogers, 1996; Zalc et al., 2001). In the last few years, the development of technology for designing flow circuits at micro lengthscales (typically ranging form 1 to 100
m), have generated renewed interest in the foundations of laminar chaotic mixing, since the flow regime in microfluidic devices is typically laminar. Microflow technology opens up new perspectives in approaching a rational design of industrial operations (Hessel et al., 2004), in the design of new micromixers (Hessel et al., 2005; Aubin et al., 2005), as well as in designing processes involving biological molecules, e.g. those associated with DNA manipulations (Tay, 2002). This novel and promising field forces the research community involved in laminar mixing theory to a critical analysis of the state of the art of the theory, with the scope of identifying the fundamental issues that should be addressed in the next few years and, hopefully, of giving a definite answer to them. In this perspective, i.e., in the light of concrete application to microdevices, even the kinematic approach suffers some conceptual drawbacks. The most important is that it does not account for diffusion, which is, ultimately, the only mechanisms that

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

ensures the approach toward a spatially homogeneous mixture. Also, even though the kinematic approach defines objective indices assessing the quality of the mixing process within each of the invariant sets associated with the flow, it does not provide any information about how this region-dependent indices should be combined together to construct a global measure of mixing in the flow domain as a whole. This last issue is clearly connected with the first in that the only flux exchange of the scalar entity (be it a chemical species or energy) that is being transported between kinematically invariant sets relies exclusively on the presence of the diffusion mechanism. In this article, we discuss the fundamental physical and mathematical aspects associated with laminar mixing and propose a unified approach referred to as the spectral theory of mixing (Cerbelli et al., 2003, 2004). The aim of this article is to show that the spectral approach provides an efficient and complete global characterization of mixing, and can be applied to general flow devices of practical interest. Throughout this article, we consider the twodimensional cavity flow as a prototypical model for describing and highlighting the theoretical and computational methods. The discussion of the computational issues in the case of three-dimensional flows is briefly addressed in the concluding section. 2. The advection–diffusion equation and flow systems Let us consider a closed domain
representing the mixing space, and the velocity field v(x, t) defined on it. The mixing space is filled with an incompressible fluid, so that ∇ ·v(x, t)=0 for any t. The physical origin of v(x, t) is immaterial in the present analysis. Let (x, t) be a scalar field representing the concentration of the transported entity (which may represent the concentration of some dye or temperature). The basic phenomenology of mixing is described by means of the advection–diffusion equation for the concentration field (x, t), which, in dimensionless form reads: j 1 = −v · ∇ + ∇ 2 , jt Pe

(1)

where Pe = VL/D is the Peclet number (V , L and D being a characteristic velocity, lengthscale and the diffusivity, respectively), and t and v are the dimensionless time and velocity field, respectively. In its ultimate essence, finding optimal mixing conditions expresses into maximizing the interplay between advection and diffusion so as to obtain a homogeneous concentration field  in the shortest possible time. Since we are considering a closed domain for which v·n|j
= 0 at the boundary j
of the domain
(n is the normal unit vector to j
), the boundary conditions for  are of Neumann type: ∇ · n|j
= 0, and Eq. (1) is equipped with the initial conditions (x, t)|t=0 = 0 (x). Throughout this article, we consider the two-dimensional cavity flow under creeping flow conditions, defined on the twodimensional unit square
= {(x, y) | 0
x, y
1} as a proto-

2755

typical model flow. The cavity flow stems from a streamfunction (x, y) (v1 = j/jy, v2 = −j/jx, v = (v1 , v2 )), which is the solution of the biharmonic equation ∇ 4  = 0, equipped with vanishing boundary conditions on  and j/jn all over the boundary but on y = 1, for which v1 |y=1 = j/jy|y=1 = 1. For the scope of this article, it is sufficient to consider an approximation for the stream-function of the cavity flow developed by Chella and Ottino (1985), which derives from a first-order expansion of the solution of the biharmonic equation obtained by applying the method of weighted residuals. Within this approximation the streamfunction (x, y) reads as (x, y) = Cx 2 (1 − x)2 k0 (y),

(2)

where C = 21, k0 (y) = A2 cosh(1 y) sin(2 y) + sinh(1 y) [B1 sin(2 y) + B2 cos(d2 y)], where 1 = 4.1503, 2 = 2.28582, and the coefficients A2 , B1 , and B2 are given by A2 =

1 sinh(1 ) sin(2 ) (2 sinh(1 ))2 − (1 sin(2 ))2

B1 = − B2 = −

,

1 cosh(1 ) sin(2 ) − 2 sinh(1 ) cos(2 ) (2 sinh(1 ))2 − (1 sin(2 ))2 A 2 2 . 1

, (3)

Within this approximation, the velocity at the boundary is given by v1 |y=1 = Cx 2 (1 − x)2 . The flow associated with the streamfunction Eq. (2) is the autonomous cavity flow, and cannot exhibit Lagrangian chaos since it is two-dimensional and time-independent. In order to improve mixing, it is natural to consider timeperiodic protocols giving rise to Lagrangian chaos. This is one of the simple “rules of thumb” representing the legacy of the kinematic theory of mixing. Let us consider the time-periodic (TP) cavity flow obtained by blinking alternately two steady flows associated with the streamfunction (1) (x, y) = (x, y) (giving rise to the velocity field v(1) ) for the first half period Tp /2, and the streamfunction (2) (x, y) = −(x, 1 − y) (velocity field v(2) ) for the second half period. The TP cavity flow defined above originates from the instantaneous switch of the upper (y = 1) and lower (y = 0) walls motion, with a wall ve(1) (2) locity oriented in positive x-direction, i.e., v1 |y=1 =v1 |y=0 = x 1 (1 − x)2 . The assumption of instantaneous switching between the two steady flows represent a fluid dynamic approximation, which is accurate at low Reynolds numbers Re < 1 (Liu et al., 1994). Fig. 1(A)–(D) show the Poincaré section (i.e., the global representation of the orbit dynamics x(n+1) = P (x(n) ), where x(n) = x(nT p ), n = 0, 1, . . .), of the TP cavity flow for several values of the flow period Tp . For lower values of Tp , the mixing space is almost entirely covered with periodic and quasiperiodic orbits (see Fig. 1(A) at Tp = 1). As the flow period increases, regions of chaotic motion appear, and progressively invade larger and larger portions of the mixing space (see Fig. 1 from (B) to (D)). For a generic value of the period Tp , the phase portrait of the kinematic equations of motion is characterized by the

2756

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

Fig. 1. Poincaré sections of TP cavity flows. (A) Tp = 1. (B) Tp = 1.5. (C) Tp = 2. (D) Tp = 4.

coexistence of regions of chaotic motion and regions of periodic/quasiperiodic behavior. Even at parameter values ensuring satisfactory “mixing performances” according to the Lagrangian paradigm, the chaotic region is massive but not ubiquitous. The occurrence of Lagrangian chaos intermingled with regions of regular motion, (this heterogeneous situation is commonly encountered in physically realizable laboratory flows or flows inside industrial equipment) restricts and limits the use of the classical quantifiers of Lagrangian chaos, such as the maximum Lyapunov exponent or related stretching-based quantities, as global performance indexes for quantifying mixing performances. In fact, the dynamics of stretch factors undergoes qualitatively different evolutions, dependently on the target invariant region considered. In order to clarify this issue, consider e.g. the Poincaré section for Tp =2 Fig. 1(C), in which a main central chaotic region C coexists with two smaller archipelagos of quasiperiodic islands. A maximum Lyapunov exponent can be defined within each invariant set (i.e., within a set A ⊆
such that P (A)=A), such as the chaotic region C depicted in Fig. 1(C). The Lyapunov exponent LLyap (C) can be computed by considering the stretching rates of vectors w = (w1 , w2 )T evolving according to the tangent dynamics dwi /dt = (jvi (x(t))jxj )wj (here x1 = x, x2 =y) coupled with the kinematic equations of motion (Beigie et al., 1994). Let w(n) (x0 ) = w(nT p ; x0 , w0 ) be the vector advected through the tangent dynamics after n periods starting from x = x0 and from an initial unit vector w0 (throughout this Section we always assume that the initial vector possesses unit norm). The Lyapunov exponent LLyap (C), restricted to C, is defined by the scaling of the quantity log |w(n) |C ∼ nLLyap (C) for large n, where the average ·C is either a spatial average performed over C or a time-average over an ergodic trajectory belonging to C. For the TP cavity flow at Tp = 2, the value of the Lyapunov exponent referred to the main chaotic region is LLyap (C) = 0.242. However, it can be readily observed that the Lyapunov exponent LLyap does not give any global indication of the stretching rate throughout the whole mixing space, since it ignores, by definition, the behavior within the periodic/quasiperiodic islands. Indeed, the only case where LLyap (C) referred to a chaotic region yields a global estimate of the stretching rate throughout the mixing space occurs if the chaotic region

possesses full measure (this situation being usually referred to as global chaotic conditions). A possible way for overcoming this problem, thus introducing global Lagrangian quantities, is to define a sort of “Lyapunov exponent” referred to the whole mixing space, and not to a specific chaotic region. This quantity, which is not a Lyapunov exponent according to the classical definition given in dynamical system theory (Eckmann and Ruelle, 1985), can be referred to as the global Lyapunov-type number NLyap , and is defined via the scaling relation of the stretching rates averaged over the whole mixing space,
log |w(n) |
= log |w(n) (x0 )|w0 dx0 ∼ nN Lyap , (4)


which holds true at large n. The inner average ·w0 is performed with respect to the initial orientations of the unit vector w0 . For example, for the TP cavity flow at Tp =2, NLyap =0.214, which is, as expected, smaller than LLyap (C) due to the contribution of quasiperiodic islands. This is because islands are characterized by vanishing Lyapunov exponent as the stretch factor within these invariant sets increases only linearly with time. We will see below (Section 4), through a concrete example that even this global quantity is insufficient and inadequate for describing mixing performances in an objective way. 3. Spectral theory of mixing In order to overcome the shortcomings of the Lagrangian theory associated with the global definition of mixing performances (“the mixing index”), it is convenient to tackle directly the problem starting from the advection–diffusion equation. Without loss of generality, we can always consider that the initial concentration field 0 (x) possesses zero mean, i.e.,
0 (x) dx = 0. This is because, the set of constant functions on
is an invariant linear subspace for Eq. (1), stemming this property from mass conservation in closed domains. A distinction should be drawn between perfectly mixed conditions as a status of a fluid continuum, and good mixedness as a dynamic process. In a static perspective, perfectly mixed conditions correspond to the occurrence of a constant concentration field throughout the mixing space. Once this condition is achieved in a closed system, neither diffusion nor advection can perturb or modify it. Although perfectly mixed conditions

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

are the ultimate goal of any mixing equipment, the characterization of these conditions is trivial and moreover it is uninteresting in mixing theory, which is focused in describing and optimizing the evolution of a fluid mixture as a dynamic process towards good mixedness. Good mixedness corresponds to the situation for which the relaxation towards perfectly mixed conditions proceeds at a sufficiently fast rate. In fact, optimization of mixing in a closed domain implies finding the stirring protocol that yields the minimum possible time to reach prescribed homogenization conditions. Elementary definitions from functional analysis are necessary in order to frame the problem in a setting that can be approached by exploiting the theory of linear operators. Let us define with L2 (
) the space of square summable functions in
and with L˙ 2 (
) ⊂ L2 (
) the subspace of L2 (
) of zero-mean functions. The square integral norm (t) of  and of its gradient ∇(t) are defined as

2 (x, t) dx, ∇2 (t) = ∇ · ∇ dx. (5) 2 (t) =





Let g(x, t) = ∇(x, t). Simple analytical manipulations of Eq. (1) yield the following evolution equations for the square integral norm of  and of its gradient (Cerbelli et al., 2004): d2 2 = − ∇2 , dt Pe
dg2 jvi 2  = −2 gi gj dx − ∇gi 2 . dt Pe
jxj

(6) (7)

i

It follows from Eq. (6) that the higher the L2 -norm of the gradient the faster the relaxation towards perfectly mixed conditions. More precisely, the relevant quantity for describing and optimizing mixing performance is the averaged ratio
1 t ∇2 () (t; Pe) = d, (8) t 0 Pe2 () since (t) = (0) exp(−t(t; Pe)). The latter observation can be exploited for constructing a useful definition of mixing as a dynamic process. Indeed, from Eqs. (1), (6)–(8) mixing can be viewed as the interplay between advection and diffusion in a continuum, the efficiency of which depends on the capability of the stirring protocol to contrast the smoothening action of diffusion in order to maintain the highest admissible gradients in the concentration field. In other words, the efficiency of a mixing protocol relies in the ability of the stirring action of pushing the concentration field towards the highest possible spatial frequencies, so that for any fixed time instant the quantity (t; Pe) is maximal. These observations explain why a spectral approach (i.e., the representation of the concentration field in terms of a functional basis of L2 (
), ordered with respect to the norm of the gradient) is the most natural setting for characterizing mixing performance and for addressing optimal mixing problems. The most natural functional basis fulfilling the requirements stated above is the eigenfunction system of the Laplacian

2757

operator defined in
and equipped with the Neumann boundary condition, which forms a basis for the functional space L2 (
). For the unit square domain
considered in this article, this eigensystem is given by B = {m,n (x) = cos(mx) cos(ny)}∞ n,m=0 , and the concentration field can  be expanded with respect to this basis as (x, t) = ∞ m,n=0 m,n (t)m,n (x). The coefficients m,n (t) = 1/( m n )
(x, t)m,n (x) dx are the spectral coefficients of  (here 0 = 1, and q = 21 for q > 0). Since B is ordered with respect to the norm of the gradient (specifically, ∇m,n 2 = m n 2 (m2 + n2 )), the spectral coefficients m,n provide a convenient representation of the spatial information on the concentration field  for addressing mixing problems. The spectral decomposition permits to express Eq. (1) in an alternative way. By enforcing a classical Galerkin expansion, the following system of equations for the Fourier coefficients is obtained dm,n 1  2 (m2 + n2 ) = A(m,n),(p,q) p,q − m,n , dt m n p,q Pe (9) where A(m,n),(p,q) = −





m,n (x)[v(x, t) · ∇p,q (x)] dx.

(10)

Since the advection–diffusion operator is compact (Reed and Simon, 1980), the number of spectral coefficients may be truncated to a finite number N = (M + 1)2 , ensuring that the truncated system yields arbitrarily accurate approximations of the original dynamics. Of course, the higher the value of Pe, the greater the number of modes N that should be considered to achieve a given accuracy. Throughout this article, we set M =30 up to M = 60, depending on the value of the Peclet number (which corresponds to N = 961 up to N = 3721 modes). By truncating Eq. (9), the advection–diffusion equation becomes a system of linear differential equations. Let c = {m,n }M n,m=0 the vector of the spectral coefficients. Eq. (1) can be expressed as dc 1 = Ac − Dc, dt Pe

(11)

where A is the matrix whose entries are A(m,n),(p,q) accounting for the advection term, and D is a diagonal matrix with entries D(m,n),(m,n) = 2 (m2 + n2 ), representing diffusion. Next, let us apply the spectral analysis to quantify the mixing action of different flow protocols. First, let us consider an autonomous flow. The matrix A entering Eq. (11) is a constant-coefficient matrix, and therefore Eq. (11) is a system of linear differential equations with constant coefficients. The decay of the concentration field is controlled by the eigenvalue of the matrix A − D/Pe possessing the largest nonvanishing real part (all the eigenvalues possess nonpositive real part because of Eq. (6), and mass conservation implies that the coefficient matrix possesses the eigenvalue = 0 associated with the constant function). This eigenvalue is referred to as the dominant eigenvalue and its real part is the dominant scaling

2758

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

exponent − controlling the asymptotic decay of concentration starting from a generic initial condition, i.e., (t) ∼ exp(−t). In the case where the dominant eigenvalue is real, by setting ∗ (x) its associated eigenfunction, it follows from Eq. (6) that  = ∇ ∗ 2 /(Pe ∗ 2 ) which indicates that the eigenvalues of the advection–diffusion operator are related to the spectral content of the associated eigenfunction. A similar approach applies for time-periodic flows, as the TP cavity flow. The instantaneous switching between two flow protocols corresponds in Eq. (11) to the switching between the two coefficient matrices A = A(1) (corresponding to the advecting field v(1) ) and A = A(1) (associated with v(2) ). By integrating these two subsystems over a flow period, and setting c(n) = c(nT p ), it follows that (2) −D/P e)T /2 p

c(n+1) = e(A

(1) −D/P e)T /2 p

e(A

c(n) = PG c(n) ,

(12)

where PG is the Poincaré operator associated with the timeperiodic advection–diffusion equation. For time-periodic flow protocols, the slowest decay rate of the concentration field is given by the eigenvalue of PG , say ∗ , possessing the largest modulus different from = 1 ( = 1 is always an eigenvalue of PG due to mass conservation, corresponding to the constant eigenfunction). In time-periodic flow protocols, the dominant decay rate , (x, t) ∼ exp(−t), is given by  = − log | ∗ |/Tp , where | ∗ | is the modulus of ∗ . Several observations follows from the analysis developed above: 1. The basic quantities associated with the spectral theory are the eigenvalue and eigenfunction spectra of the advection–diffusion operator A − D/Pe (for autonomous flows) or of the Poincaré operator (for time-periodic flows), and specifically the dominant scaling exponent . This quantities describe global mixing properties in the entire flow domain
. 2. The spectral approach provides a way to set the optimization problem of mixing performances: given a time instant t, the optimal stirring conditions can be defined as the flow protocol maximizing the quantity (t; Pe). 3. The spectral approach provides the functional setting in which both convection and diffusion and (hypothetical) pure convection (corresponding to Pe → ∞ or, more precisely 1/Pe = 0). In the absence of diffusion, Eq. (6) becomes d2 /dt = 0, which implies that the square integral norm is constant, 2 (t) = 2 (0). In the absence of diffusion, the optimal stirring protocol maximizes the integral quan tity −
∇v : gg dx, (first term at the right-hand side of the Eq. (7)), which is equivalent to say that the gradient of the concentration function should be maximized. Moreover, the dominant eigenfunctions for Pe → ∞ provide a global template of the action of stirring throughout the whole mixing space. 4. It is important to point out that the spectral analysis described above is a theoretical framework for describing mixing processes on a global and objective way, and it is not a computational strategy. Although in solving

the advection–diffusion equation we have explicitly used a Galerkin (spectral) decomposition, this computational strategy is, in principle, uncorrelated to the spectral analysis. More precisely, the spectral quantities describing mixing, such as the eigenvalues/eigenvectors of the advection/diffusion operator, or the quantity (t; Pe) could be equally well obtained starting from any other computational method (e.g. finite difference/volume/element methods) applied to solve Eq. (1). This observation is important, since a spectral characterization can be performed on equal footing starting from CFD simulations of mixing equipments possessing more complex geometries (stirred tanks, static mixers) for which the use of the generalized Fourier expansion for solving Eq. (1) may be too cumbersome. 5. An objection against spectral analysis as compared to the Lagrangian theory of mixing is that it transforms a simple set of ordinary differential equations into a partial differential model. This could be viewed as a way to make it more complex instead of simplifying it. Indeed, it is easy to show that it is rather the opposite, essentially due to the fact that the advection–diffusion equation is a linear equation, and the results obtained out of it are of global nature.

4. Results In this Section we discuss some results deriving from the spectral analysis of the cavity flow. Due to space limitations, we pinpoint only few specific issues, leaving a more comprehensive treatment to future works. Fig. 2(A) shows the behavior of the dominant scaling exponent  for the autonomous cavity flow. A preliminary observation is that the eigenvalue spectrum can be decomposed into several distinct eigenvalue branches, each of which characterized by a different scaling behavior with the Peclet number (for Pe → ∞). For the autonomous cavity flow, two distinct branches can be observed, corresponding, respectively, to real and complex conjugated eigenvalues, the dominant exponents of which are depicted in Fig. 2(A) (dotted lines (a) and (b), respectively). At large Pe, the real branch is characterized by a diffusive scaling  ∼ 1/Pe, inversely proportional to the Peclet number, while the complex conjugated branch displays a convection-enhanced scaling  ∼ Pe−1/2 . This phenomenon is a further confirmation of the results obtained for simpler prototypical model flows by Cerbelli et al. (2004) and find a theoretical explanation in the results derived by Giona et al. (2004). It is useful to connect the scaling behavior observed in the dominant exponent  to the spatial structure of the eigenfunctions. Fig. 3(A)–(B) shows the spatial structure of the (1) eigenfunction r (x) (throughout this section, all of the eigenfunctions considered are normalized to unit L2 -norm) and (1) of the modulus of the eigenfunction cc associated with the dominant eigenvalues of the real and complex conjugated branch for Pe = 103 . The real eigenfunction (associated with the diffusive eigenvalue scaling) resembles the structure of the stream function. Conversely, the dominant eigenfunction of the

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

101

102 c

b

Λ

a

e

10-2

Λ

100

10-4 10-1 (A)

2759

d

10-1

d

101

103

c 10-3 100

105

Pe

(B)

a,b 101

102

103

104

105

Pe

Fig. 2. (A) Autonomous cavity flow. Dominant scaling exponent  vs Pe. Dotted line (a) and (•) refer to the complex conjugated branch of the spectrum, dotted line (b) and (◦) to the real branch. Lines (c) and (d) represent the diffusive scaling  ∼ Pe−1 , line (e) the convection-enhanced scaling  ∼ Pe−1/2 . (B) Dominant scaling exponent  for TP cavity flows. Lines (a)–(b) refer to TP = 1.0 and Tp = 1.5, respectively. Line (c) and (◦) to Tp = 2.0. Line (d) and (•) to Tp = 4.0. The thicker lines depict the power-law scaling Eq. (13), where the values of exponents for each flow protocol are reviewed in Table 1.

(1)

Fig. 3. Three-dimensional plots of the eigenfunctions of advection–diffusion operator for the autonomous cavity flow. (A) Dominant eigenfunction r of the (1) diffusion-controlled branch for Pe = 103 . (B) Modulus of the dominant eigenfunction cc of the convection-enhanced branch for Pe = 103 . (C) Dominant (1) (2) 4 eigenfunction r of the diffusion-controlled branch for Pe = 10 . (D) Second dominant eigenfunction r of the diffusion-controlled branch for Pe = 104 .

complex–conjugated branch displays, (as expected) higher gradients, localized near the moving wall (y = 1). As Pe increases, the dominant eigenfunctions of the diffusive branch (the real branch) converge towards an invariant shape the contour plot of which corresponds to that of the streamfunction. This is depicted in Fig. 3(C) showing the dominant real eigenfunction (1) r (x) at Pe=104 . This result is not surprising, since the family the dominant eigenfunctions of the diffusive branch are characterized by the property limP e→∞ ∇ r (x) = C = constant, i.e., their gradients converge towards a constant quantity. The higher-order eigenfunctions of the diffusive branch of the spectrum can be viewed as higher-order modulations of a fundamental shape possessing as level sets the streamlines of the stirring flow. For example, Fig. 3(D) shows the second (2) dominant eigenfunction r (x) at Pe = 104 . Let us consider some features of time-periodic protocols, for which the spectral analysis yields interesting and apparently unexpected results. Consider the TP cavity flow for Tp = 1, 1.5, 2, 4 (the Poincaré sections associated with these protocols are depicted in Fig. 1), and the time-behavior of the homogenization dynamics (represented by the graph of the function (t) vs t, where  is the solution of the advection–diffusion equation) at Pe = 103 , starting from a segregated initial condi-

√ √ tion 0 (x) = 1/ 2 for x < 21 , and 0 (x) = −1/ 2 for x > 21 . This initial profile can be viewed as the representation of two completely segregated reacting species A and B (indeed the advection–diffusion equation describe also an infinitely fast chemical reaction A + B → products, in the presence of equal diffusivities of A and B, and in this case  = cA − cB , i.e., the field  is the difference of the concentrations of the two reacting species). Fig. 4(A) shows the decay of the L2 -norm of (x, t) for several TP protocols, compared to the corresponding decay observed in the autonomous cavity flow. The behavior of the quantity ∇2 /(Pe2 ) (which is the integrand of the function (t; Pe) Eq. (8)) for the same numerical experiments is depicted in Fig. 4(B). An interesting phenomenon may be observed: time-periodic protocols (curves (b)–(d) in Fig. 4(A)) which give rise to Lagrangian chaos (at least in some regions of the mixing space) behave significantly worse, as it regards mixing performances (viewed as timescale to achieve a certain degree of homogenization) than the corresponding autonomous flow protocol (Fig. 4 line (a)). For example, the TP cavity flow protocol at Tp = 2.0, which possesses a major chaotic region, and a positive global Lyapunov-type number significantly greater than 0 (see Table 1), shows a significantly slower decay rate in the

2760

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

100

||∇ φ||2 (t)/(Pe ||φ||2)

100

||φ|| (t)

10-1

d c b

10-2 e

e 10-1 a

b-d

a 10-3

0

20

40

(A)

60

80

10-2

100

0

(B)

t

20

40

60

80

100

t

Fig. 4. TP and autonomous cavity flows at Pe = 103 . (A) (t) vs t. (B) ∇ 2 /(Pe2 ) vs t. Thicker line (a) refers to the autonomous cavity flow, line (b) to TP cavity flow at Tp = 1, line (c) to TP cavity flow at Tp = 1.5 line (d) to TP cavity flow at Tp = 2, line (e) to TP cavity flow at Tp = 4.

Table 1 Lyap , of the dominant scaling Review of the Lyapunov numbers NLyap , N exponent  of the advection–diffusion equation at Pe = 103 , and of the scaling exponent Eq. (13) for the cavity flow protocols (Aut. stands for Autonomous.) Flow protocol

NLyap

Lyap N

 (Pe = 103 )

Exp. Eq. (13)

Aut. Tp = 1.5 Tp = 2.0 Tp = 4.0

0 0.09 0.21 0.73

0 0.06 0.105 0.185

0.0595 0.0320 0.0343 0.1250

1 0.86 0.69 0.60

advection term (expressed by the matrix A in Eqs. (11)–(12) and a diffusive term, the matrix D), the spectral properties controlling the concentration decay depend in a very complicated manner on the interplay between these two terms. For instance, the convection-enhanced scaling  ∼ Pe−1/2 is a first, clear example of this complex interaction. Indeed, more complicated scaling behaviors occur for TP flow protocols. Fig. 2(B) shows the behavior of the dominant exponent  for the TP cavity flows. In general, a power-law scaling  ∼ Pe−

relaxation towards homogenization than the autonomous cavity flow. This phenomenon does not refer exclusively to the asymptotic behavior (when the L2 -norm is small), but can be significant even in the early stages of the homogenization process (i.e., for time-scales for which (t)/(0) > 5 × 10−2 ). This apparent unexpected phenomenon occurs for a broad range of Pe values ranging from 500 up to 5000. To complete the analysis, for Tp = 4.0, one may observe better mixing performance than in the autonomous case (see Fig. 4(A), line (e)). Table 1 reviews the actual values global parameters characterizing these flow protocols, namely Lyapunov type number NLyap , the rescaled Lyap = NLyap /Tp , the dominant scalLyapunov-type numberN ing exponent . Several conclusions can be drawn from the above case study. (i) As it regards homogenization performances, the occurrence of chaotic behavior (Lagrangian chaos) is not necessarily an indication of better mixing performances with respect to flow protocols displaying nonchaotic kinematics. (ii) This phenomenon is due to the complex interplay between advection and diffusion occurring within regions of regular and chaotic motion. (iii) Global kinematic quantities (e.g. the Lyapunov-type number) are unable to capture the essence of mixing process i.e., the relaxation rates towards perfectly mixed conditions, and, consequently, cannot be viewed as reliable “mixing indexes”. Although it is fully correct to view the rates controlling mixing (homogenization) as the superposition of two terms: the

(13)

is observed at large Pe values, with a value of the scaling exponent bounded between 0 and 1. The value of the exponent

is very sensitive to the flow protocol. A review of the value attained by for the TP cavity flows in reported in Table 1. To date, no theory is available for explaining and predicting the value of the exponent in the presence of a time-periodic stirring. 5. Concluding remarks In this article, we have addressed the fundamentals underlying the spectral theory of mixing, which is essentially based on the analysis of the advection–diffusion equation coupled with few elementary concepts deriving from functional analysis. We have avoided unnecessary functional-theoretical complexities, and tried to develop the topic in a form suitable to the widest possible audience. Albeit this approach may appear abstract at a first insight, one should consider that it readily yields quantitative information about the characteristic timescales to achieve homogenization (which is essentially the goal of any process engineer). Besides, the strength of the functional analysis approach is that it allows to express the mixing process in terms of a linear system of ordinary differential equations, for which a manifold of refined and accurate numerical methods exist for solving eigenvalue/eigenfunction problems for large matrices.

S. Cerbelli et al. / Chemical Engineering Science 61 (2006) 2754 – 2761

In this article, we focus on a relatively simple, physically realizable flow protocol, namely the cavity flow. The same approach can be applied to more complex flow protocols (e.g. flows in static mixers, stirred tanks, or electromagnetically driven microdevices). The only issue that matters in addressing complex mixing devices in the size (dimension) of the matrices (A and D) considered. The use of parallel numerical computation is certainly the most straightforward way to implement this approach to specific mixing devices of practical interest. Moreover, the advection–diffusion equation can be simulated within commercial CFD packages used to determine the velocity field in industrial devices, and the only additional step required for spectral analysis is the implementation of iterative methods for obtaining the dominant eigenvalue/eigenfunction. It should be observed that even in the cases where eigenvalue/eigenfunction spectra are of cumbersome computation (e.g. mixers with complex geometries), the use of the spectral analytical approach provides an objective characterization of the intrinsic timescales to homogenization (eigenvalues) and of the resulting mixing patterns (eigenfunctions). This is of direct practical use not solely for microflow devices, but also for more classical mixing units (work is in progress as it regards the application to stirred tanks). A byproduct of the spectral analysis is the formulation of an objective function, possessing a direct physical meaning, namely the function (t, Pe) defined in Eq. (8), that can be used for optimizing stirring protocols. Notation A C D LLyap (C) n

advection matrix in the spectral representation invariant chaotic region diffusion matrix in the spectral representation Lyapunov exponent referred to C iteration order in the application of time-periodic protocols NLyap global Lyapunov-type number PG Poincaré operator associated with a time-periodic advection–diffusion equation Pe Peclet number t time Tp flow period x=(x, y) position vector Greek letters

(x) (t; Pe) 

scaling exponent, see Eq. (13) eigenfunctions of the advection–diffusion operator instantaneous decay rate, see Eq. (8) dominant scaling exponent of the concentration decay

 m,n (t) (t) m,n (x) 


2761

concentration field spectral coefficients of  L2 -norm of  at time t eigenfunctions of the Laplacian operator streamfunction mixing space

References Aref, H., 1984. Stirring by chaotic advection. Journal of Fluid Mechanics 143, 1–21. Aubin, J., Fletcher, D.F., Xuereb, C., 2005. Design of micromixers using CFD modelling. Chemical Engineering Science 60, 2503–2516. Baldyga, J., Bourne, J.R., 1999. Turbulent Mixing and Chemical Reaction. Wiley, Chichester. Beigie, D., Leonard, A., Wiggins, S., 1994. Invariant manifold templates for chaotic advection. Chaos, Solitons and Fractals 4, 749–868. Chella, R., Ottino, J.M., 1985. Fluid mechanics of mixing in a singlescrew extruder. Industrial and Engineering Chemistry Fundamentals 24, 170–180. Cerbelli, S., Adrover, A., Giona, M., 2003. Enhanced diffusion regimes in bounded chaotic flows. Physics Letters A 312, 355–362. Cerbelli, S., Vitacolonna, V., Adrover, A., Giona, M., 2004. Eigenvalue–eigenfunction analysis of infinitely fast reactions and micromixing regimes in regular and chaotic bounded flows. Chemical Engineering Science 59, 2125–2144. Danckwerts, P.V., 1952. The definition and measurements of some characteristics of mixtures. Applied Science Research A 3, 279–296. Danckwerts, P.V., 1958. The effect of incomplete mixing on homogeneous reactions. Chemical Engineering Science 8, 93–99. Eckmann, J.-P., Ruelle, D., 1985. Ergodic theory of chaos and strange attractors. Reviews of Modern Physics 57, 617–655. Giona, M., Cerbelli, S., Vitacolonna, V., 2004. Universality and imaginary potentials in advection–diffusion equations in closed flows. Journal of Fluid Mechanics 513, 221–237. Harvey III, A.D., Rogers, S.E., 1996. Steady and unsteady computation of impeller-stirred reactor. A.I.Ch.E. Journal 42, 2701–2712. Hessel, V., Hardt, S., Löwe, H., 2004. Chemical Micro Process Engineering. Wiley-VCH, Weinheim. Hessel, V., Löwe, H., Schönfeld, F., 2005. Micromixers—a review on passive and active mixing principles. Chemical Engineering Science 60, 2479–2501. Liu, M., Peskin, R.L., Muzzio, F.J., Leong, C.W., 1994. Structure of the stretching field in chaotic cavity flows. A.I.Ch.E. Journal 40, 1273–1286. Müller, S.D., Mezic, I., Walther, J.H., Koumotsakos, P., 2004. Transverse momentum micromixer optimization with evolution strategies. Computers and Fluids 33, 521–531. Ottino, J.M., 1989. The Kinematics of Mixing, Stretching, Chaos and Transport. Cambridge University Press, Cambridge. Reed, M., Simon, B., 1980. Methods of modern mathematical physics: I—Functional analysis. Academic Press, New York. Tay, F.E.H., 2002. Microfluidics and BioMEMS Applications. Kluwer Academic Press, Boston. Zalc, J.M., Alvarez, M.M., Muzzio, F.J., Arick, B.E., 2001. Extensive validation of computed laminar flow in a stirred tank with three Rushton turbines. A.I.Ch.E. Journal 47, 2144–2154.