Eigenmode analysis of advective–diffusive transport by the compact mapping method

European Journal of Mechanics B/Fluids 49 (2015) 1–11

Contents lists available at ScienceDirect

European Journal of Mechanics B/Fluids journal homepage: www.elsevier.com/locate/ejmflu

Eigenmode analysis of advective–diffusive transport by the compact mapping method O. Gorodetskyi, M.F.M. Speetjens, P.D. Anderson ∗ Department of Mechanical Engineering, Eindhoven University of Technology, The Netherlands

article

info

Article history: Received 7 June 2013 Received in revised form 5 June 2014 Accepted 30 June 2014 Available online 15 July 2014 Keywords: Mapping method Eigenmode decomposition Diffusion Mixing

abstract The present study concerns an efficient spectral analysis of advective–diffusive transport in periodic flows by the way of a compact version of the diffusive mapping method. Key to the compact approach is the representation of the scalar evolution by only a small subset of the eigenmodes of the mapping matrix, and capturing the relevant features of the transient towards the homogeneous state. This has been demonstrated for purely advective transport in an earlier study by Gorodetskyi et al. (2012). Here this ansatz is extended to advective–diffusive transport and more complex 3D flow fields, motivated primarily by the importance of molecular diffusion in many mixing processes. The study exposed an even greater potential for such transport problems due to the progressive widening of the spectral gaps in the eigenvalue spectrum of the mapping matrix with increasing diffusion. This facilitates substantially larger reductions of the eigenmode basis compared to the purely advective limit for a given approximation tolerance. The compact diffusive mapping method is demonstrated for a representative three-dimensional prototype micro-mixer. This revealed a reliable prediction of (transient) scalar evolutions and mixing patterns with reductions of the eigenmode basis by up to a factor 2000. The accurate estimation of the truncation error from the eigenvalue spectrum enables the systematic determination of the spectral cut-off for a desired degree of approximation. The validity and universality of the presumed correlation between spectral cut-off and truncation error have been established. This has the important practical consequence that the cut-off can a priori be chosen such that the truncation error remains within a preset tolerance. This offers a way to systematically (and reliably) employ the compact mapping method for an in-depth analysis of advective–diffusive transport. © 2014 Elsevier Masson SAS. All rights reserved.

1. Introduction Fluid mixing is a widespread and complex phenomenon observed in a great variety of natural and man-made processes at length scales ranging from microns to hundreds of kilometers. Fluid mixing can be categorized according to the flow and mixing conditions or the fluid properties. Many industrial applications concern mixing under purely laminar flow conditions. This regime is characterized by a low Reynolds number Re = UL/ν , with U and L typical velocity and length scales, respectively, and ν the kinematic viscosity, and can occur basically due to either highly viscous fluids or very small system dimensions. Examples of the former include chemical engineering, cosmetic industry, food processing, polymer and pharmaceutical processing [1,2]; the latter encompasses the highly-relevant field of micro-fluidics, which is at the

∗

Corresponding author. E-mail address: [email protected] (P.D. Anderson).

http://dx.doi.org/10.1016/j.euromechflu.2014.06.007 0997-7546/© 2014 Elsevier Masson SAS. All rights reserved.

heart of many emerging technologies. Here Reynolds numbers often are in the order of unity or less [3–6]. Efficient mixing is the ultimate goal in the overwhelming majority of micro-fluidic systems. However, turbulence, the prime generator of mixing in macroscopic flows, usually is unattainable in micro-scale systems. Numerous studies [7–12] have demonstrated that in laminar flows efficient mixing can be accomplished by an alternative mechanism: chaotic advection. This is the process of exponential stretching and folding of material fluid regions, which enhances mixing by creating an as large as possible ‘‘working area’’ for molecular diffusion to exchange scalar quantities (e.g. chemical species, reactants, nutrients) between fluid parcels. Hence, mixing in laminar flows is an interplay between advection and diffusion. The mixing efficiency is determined strongly by the stretching rate of material lines and surfaces by the advective transport [13–15]. Stretching and folding by chaotic advection rapidly creates thin fluid filaments with extremely long interfaces that promote fast homogenization of the mixture. However, despite the essential part of molecular diffusion in the mixing

2

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

Fig. 1. Computation of the coefficients Φij of the mapping matrix 8. Initially (panel (a)) a cell Ωj is uniformly filled by Mj particles that are tracked over one period (panel (b)). The ratio of the number of markers received by the recipient cell Ωi to the initial number of markers in Ωj determines Φij (here Φij = 3/16).

process, it has received far less attention than advection. Notable exceptions are the few studies like by Raynal et al. [16], Christov et al. [17], Schlick et al. [18]. This underexposed role of diffusion motivates this present study. Mixing studies on the basis of the spectrum of the linear advection–diffusion operator have proven efficient and successful. These methods in essence describe the evolution and decay rate in terms of the eigenmodes of this operator [19–21]. The eigenmodes either decay exponentially or persist indefinitely, causing initial fields to always follow the same scenario: progression towards an asymptotic steady via a transient that quickly becomes dominated by the slowest-decaying eigenmode. The latter is named ‘‘dominant’’ (or ‘‘strange’’) eigenmode and is the principal subject of investigation in many studies [22–32]. Systems with moderate to high diffusion typically have one distinct dominant mode. However, in the advective limit (i.e. for very weak diffusion), usually subsets of modes emerge with comparably slow decay rates and, inherently, similar impact on the transient. These are termed ‘‘transient eigenmodes’’ hereafter [33]. Spectral analyses, though a powerful tool, are difficult to perform for realistic flows in complex geometries. Hence, studies to date primarily concern highly-idealized systems [25,19,26]. An efficient and versatile numerical strategy is an absolute necessity for making this approach accessible to realistic configurations. The mapping method offers a solution as a well-established simulation technique for complex (industrial) mixing problems in timeperiodic and spatially-periodic laminar flows [34–41]. The link between the spectral properties of the mapping matrix and the transport operator in the purely advective limit has been explored by Singh et al. [42] and Speetjens et al. [43]. This has recently been extended to advective–diffusive processes by the so-called diffusive mapping matrix [44–46]. Gorodetskyi et al. [33] have demonstrated that the scalar evolution in the purely advective limit can, at least for practical analyses, be sufficiently accurately reproduced by a subset of transient eigenmodes that is substantially smaller than the whole eigenspace. The present study extends the compact mapping method to the case of advective–diffusive transport (denoted ‘‘diffusive compact mapping method’’). Diffusion drastically reduces the subset of eigenmodes in the compact method for a given truncation error and thus yields an enormous saving in eigenspectrum computations compared to the purely advective case. (Computation times for subsets of M eigenmodes typically increase exponentially: t ∼ expM .) This is of great practical importance, since mapping matrices normally are huge, sparse and non-symmetrical, rendering evaluation of (subsets of) the eigenspectrum extremely costly if not impossible altogether. Thus the diffusive compact mapping method, by capturing the relevant physics in manageable subsets of eigenmodes, is particularly attractive for spectral decomposition and analysis of mixing processes. Mixing simulation using an eigenmode representation is in itself computationally more expensive than directly using the original mapping matrix due to the considerable effort involved with, first, the spectral decomposition of the matrix [42] and, second, evaluation of the expansion coefficients [33] for given initial conditions. The diffusive compact method offers tremendous reduction [45,46] compared to the purely advective compact method

yet remains relatively costly. Hence, computational considerations alone are insufficient to justify this approach. The benefit of the compact method must be sought at a more fundamental level. The essential transport characteristics of a system are (at least for practical purposes) namely embedded in a subset M ≪ N of the N eigenmodes, which can be exploited in two ways. First one, for deepening physical understanding of fluid dynamics and mixing process by examining the spatio-temporal behavior of the retained modes. Second, for the optimization and control of the mixing process. Representation by a compact eigenmode basis namely enables inversion of the problem with sufficient accuracy and acceptable cost. This facilitates the systematic determination of the initial conditions required to attain a desired final state and may thus potentially revolutionize design procedures for mixing processes and equipment. (Design and optimization happens to date mainly on an empirical basis by experimentally determining which initial conditions give the ‘‘best’’ final state.) The present study takes a first step in this direction by investigating (i) the properties of the compact eigenmode basis specifically in the presence of diffusion and (ii) the manner in which this connects initial and final states. The compact diffusive mapping method is demonstrated by a representative three-dimensional (3D) prototype micro-mixer: the staggered herringbone micromixer (SHM) [47,48]. Diffusion is especially important in such micro-fluidic devices due to the smaller Péclet numbers Pe = UL/D, with U , L as before and D the diffusivity, compared to macroscopic systems. The discussion is organized as follows. Section 2 briefly reviews the (diffusive) mapping method. The SHM is introduced in Section 3. The general spectral properties of the diffusive mapping matrix associated with the SHM and its reduction to a compact form are discussed in Section 4. The performance of the compact mapping method is investigated in Section 5. Conclusions are drawn in Section 6. 2. The mapping method for advective–diffusive mixing The mapping method is based on the tracking of sets of particles for specific intervals in time or space (period Tp ). To this end the flow domain is subdivided into small cells {Ωk }k=1,N that are uniformly filled with particles. The mapping matrix 8 captures the transport in terms of the redistribution of the particles over the cells within one period. This approach is shown schematically in Fig. 1. The original mapping method concerns purely advective mixing and has recently been extended to advective–diffusive transport by Gorodetskyi et al. [44]. The diffusive method is based on the stochastic formulation of advective–diffusive transport. Here diffusion is incorporated into the Lagrangian equations of motion for fluid parcels via a stochastic term that represents the random Brownian motion due to molecular effects. This is described by the non-dimensional stochastic differential equation,

 dx(t ) = v(x(t ), t ) dt +

2 Pe

dξ(t ),

(1)

with Pe the Péclet number introduced before and dξ(t ) infinitesimal increments of a d-dimensional Wiener process (with d the dimension of the flow domain) [44]. Eq. (1) can be solved numerically

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

b

3

c

a

Fig. 2. Schematic of the staggered herringbone micro-mixer (SHM): (a) configuration and dimensions; (b) alternating groove patterns on bottom wall for one mixing unit; (c) cross-sectional streamfunctions of the transverse flow for each of the groove patterns.

a

b

c

Fig. 3. Cross-sectional advection patterns at the inlet of each mixing unit (visualized by Poincaré sections) as a function of the non-dimensional separatrix position: (a) Π = 0.2/2.35; (b) Π = 0.65/2.35; (c) Π = 1/3.

by an Euler–Langevin explicit scheme [49], according to

 xn+1 = xn + vn h +

2 Pe

h1/2 rn ,

(2)

where h is a time step and rn = (r1,n , . . . , rd,n ) is a realization of the Wiener process. The study in [44] demonstrated the close agreement between (1) and the Eulerian advection–diffusion equation. Hence, the adopted extension of the mapping method adequately incorporates diffusion in the transport simulations. The diffusive mapping matrix is constructed in essentially the same way as in the original technique by tracking particles for one period Tp via numerical integration of Eq. (1). The scalar concentration C (x, t ) is described by the vector C ∈ RN ×1 , representing a coarse-grained description of C (x, t ), where each vector component corresponds with the cell-averaged concentration. The mapping matrix enables efficient computation of the redistribution of concentration Cn after n periods following Cn = (8(8(. . . (8 C0 ) . . .),





n times

(3)

transverse flow. Here two consecutive groove patterns constitute one mixing cycle. Fig. 2(a) gives a 3D schematic representation of the SHM. The flow inside the SHM is described by the analytical Stokes flow (Re = 0) proposed by Stroock and McGraw [48]. The transverse flow due to the grooves is controlled by the location of the separatrix (position c in Fig. 2(b)) and is characterized by the non-dimensional parameter Π = c /a. An analytical solution for the transverse velocity field is obtained by the method of superposition [50,46]. The non-dimensional parameters are according to Stroock and McGraw [48] and achieve the closest agreement with experiments [47]. The pronounced elongation of the micro-channel strongly reduces the dependence on the longitudi− → nal coordinate z: u = (ux (x, y), uy (x, y), uz (x, y)). Thus the fully developed 3D flow can be represented as a superposition of transverse and axial flows. Please refer to [46] for details and explicit expressions. 3.2. Asymptotic advection patterns



for any initial state C0 . 3. Prototype micro-mixer and properties of the flow 3.1. Mixer configuration The pressure-driven static staggered herringbone micromixer (SHM) is adopted as a prototype micro-mixer. The SHM has originally been proposed in the experimental study by Stroock et al. [47] and comprises a rectangular 3D micro-channel with an off-centered herringbone pattern on the bottom wall. This pattern is created by sets of grooves that are asymmetrically arranged with respect to the centerline of the channel and induce transverse flow. The asymmetry switches every half cycle, thus resulting in a systematic alternation between two patterns and, in turn, a periodic

The model parameters have been tuned to the experimental system (Section 3.1) and remain fixed, except for the nondimensional position Π = c /a of the separatrix of the transverse flow. The advective dynamics of the flow depend significantly upon Π . This is demonstrated in Fig. 3 by way of Poincaré sections. This exposes sizeable regular zones for lower Π = c /a and their diminution with the separatrix becoming closer to the side wall. Detailed investigation of such regular zones in the purely advective limit is in [46]. Moreover, their correlation with pattern formation in the advective–diffusive transport is determined via spectral analysis of the diffusive mapping matrix. 4. Spectral properties of the diffusive mapping matrix The discussion below concentrates on the connection between the advective–diffusive transport properties and the eigenmodes

4

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

of the diffusive mapping matrix as a function of the Péclet number Pe. The behavior during the transient between initial and asymptotic states of the system is of particular interest, since many realistic systems remain in this stage for a significant portion of the duration of the entire mixing process. The study in Gorodetskyi et al. [33] demonstrated that in the purely advective limit the transient is dictated predominantly by a small subset of decaying eigenmodes (i.e. the beforementioned ‘‘transient eigenmodes’’). This notion resulted in the compact mapping matrix based upon truncation of the eigenmode expansion. Here this approach is extended to the diffusive mapping matrix. 4.1. General properties of the matrix eigenspace The diffusive mapping matrix 8 describes the redistribution of the scalar field C by advective–diffusive transport over a finite mesh (N cells) during one period Tp in a periodic flow (Section 2). The evolution over multiple periods follows from the repeated action of 8 according to (3). Important to stress is that, notwithstanding periodic flow, the evolution of the scalar field typically is not periodic, but instead undergoes a transient towards a homogeneous asymptotic state. The particle-based approach towards diffusion following (1) nonetheless admits representation of this process by way of a cyclic operation. Let {λk , νk }Nk=1 be a finite set of eigenvalues and eigenvectors of the diffusive mapping matrix 8. Boundedness of the scalar distribution restricts the eigenvalue spectrum {λk }Nk=1 to within the unit circle, i.e. |λk | ≤ 1, admitting an ordering of the eigenvalues in a descending way according to |λ1 | > |λ2 | ≥ · · · ≥ |λk | ≥ · · · ≥ |λN |. Mass conservation causes the spectrum to always possess a trivial persistent eigenmode with |λ1 | = 1 and constant eigenvector ν1 = const [42,33,43]. The eigenvalues λk for k ≥ 2 show a gradual decay (|λk | < 1), signifying an exponential decay of the associated eigenmodes and, inherently, an exponential decay of any initial scalar field towards a homogeneous state [19]. This decay rapidly becomes dictated by the slowest decaying eigenmode(s), i.e. with magnitude(s) |λk | closest to one, denoted as ‘‘dominant eigenmode’’ in the literature. The advection–diffusion operator strictly becomes a unitary operator in the purely advective limit, implying an eigenvalue spectrum restricted to the unit circle. However, numerical diffusion due to finite cell sizes causes the spectrum of the purely advective mapping matrix to nonetheless exhibit exponential decay. This numerical diffusion by coarse-graining has been studied extensively in Gorodetskyi et al. for various systems including the SHM [44–46]. Coarse-graining in fact acts as a superimposed diffusive contribution that can be represented by an effective Péclet number Peeff . Hence, the purely advective mapping matrix strictly describes advective–diffusive transport at Pe = Peeff , where the latter can be estimated numerically [44–46]. 4.2. Role of the advection characteristics The spectral properties of the mapping matrix, and consequently, the mixing quality are directly related to the flow dynamics. Key question is the role of chaotic advection in this process. Consider to this end two typical SHM protocols, viz. Π = 0.2/2.35 and Π = 1/3, which yield partially chaotic and chaotic dynamics, respectively. Case Π = 0.2/2.35 is characterized by two large regular islands (Fig. 3(a)), signifying transport barriers and thus inefficient mixing. The presence or absence of such islands has a clear impact on the spectral properties of the mapping matrix. Fig. 4 shows the dominant eigenvalues |λ2 | versus Pe for cases Π = 0.2/2.35 (markers ‘1’) and Π = 1/3 (markers ‘◦’). (The properties for Pe > 106 are determined from the purely advective mapping matrix at Pe = Peeff .) This exposes a marked increase of

Fig. 4. Modulus of the first non-trivial eigenvalue |λ2 | versus Péclet number Pe for protocol Π = 0.2/2.35 (‘1’) and Π = 1/3 (‘◦’).

|λ2 | with Pe, implying a slower transient. However, this increase is less profound for chaotic dynamics, meaning that regular zones promote this slowing down. The decay rate of the non-trivial part of the spectrum, |λk | (k ≥ 2), depends strongly on Pe, which has important ramifications for the scalar field [46]. Fig. 5 shows the first three non-trivial eigenvalues {|λ2 |, |λ3 |, |λ4 |} of the diffusive mapping matrix versus Pe for protocols Π = 0.2/2.35 (panel (a)) and Π = 1/3 (panel (b)). This reveals that the above findings on λ2 hold for the entire spectrum. Fig. 5 reveals a widening spectral gap between consecutive eigenvalues with decreasing Pe, meaning that for lower Pe (up to Pe ∼ 104 ), the evolution rapidly becomes dominated by mode k = 2. For instance, at Pe ∼ 103 , the spectral gap between k = 2 and k = 3 has already widened to |λ2 |n /|λ3 |n ∼ 102 after only n ≈ 10 periods. However, this widening substantially diminishes with increasing Pe, meaning that modes k > 2 remain an important factor in the evolution for longer periods of time. Thus for many practical systems, where processes typically have a limited duration, the influence of eigenmodes other than the dominant eigenmode k = 2 must be taken into account for a reliable spectral analysis [33]. 4.3. Compact eigenmode representation of the scalar evolution The coarse-grained description of the scalar distribution according to (3) of concentration Cn admits the eigenmode decomposition Cn =

N 

Ck0 λnk νk ,

(4)

k =1

with Ck0 the expansion coefficients, which are fully determined by the initial state, and {λk , νk }Nk=1 the set of N eigenvalue– eigenvector pairs of the mapping matrix [42,43]. However, typically only a small subset k ≤ M (M ≪ N ) of eigenmodes is relevant, comprising the slowest decaying eigenmodes (denoted as ‘‘transient eigenmodes’’ following [33]). Hence, at least for practical analyses, the scalar evolution can be described by the truncated basis k ≤ M to a sufficient degree of accuracy. This admits reduction of decomposition (4) to the compact form Cn =

M 

Ck0 λnk νk + ϵn ,

(5)

k=1

with corresponding truncation error ϵn ∼ O (|λM +1 |n ) after n periods [33]. Thus for n ≥ 1 the compact form (5) approximates the true scalar evolution with precision ϵn . The truncation error scaling

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

a

5

b

Fig. 5. Modulus of the non-trivial eigenvalues |λ2 | (‘△’), |λ3 | (‘’) and |λ4 | (‘+’) versus Pe: (a) protocol Π = 0.2/2.35; (b) protocol Π = 1/3. Table 1 Computational time for the eigenproblem versus cut-off (M /N) of the eigenspectrum and size of the mapping matrix. N ×N

2500 × 2500 10 000 × 10 000 20 000 × 20 000

M /N 1%

2%

5%

10%

20%

1 min 8 min 1 h 13 min

2 min 1 h 15 min 6 h 56 min

10 min 10 h 38 min 2 days 5 h

25 min 23 h 14 min >4 days

>2 days

1 h 5 min

as ϵn ∼ O (|λM +1 |n ) is a crucial assumption underlying the compact approach. Validity of this scaling has been established for advective transport the model flow in [33] and will be examined specifically for advective–diffusive transport in the SHM in Section 5.2. Spectral cut-off M and starting period n are interdependent and basically user-defined. Specifying M implicitly determines the period n beyond which accuracy ϵn is attained; conversely, specifying a tolerance ϵn implicitly determines the cut-off M. Once parameters (M , n) are known, the expansion coefficients Ck0 (k ≤ M ) can be evaluated by the technique proposed in [33]. First exploratory investigations of an idealized 2D flow system revealed that reduction factors of M /N = 10% are possible while maintaining an acceptable degree of approximation for (industrial) applications. Typically, the size of the mapping matrix for rather coarse meshes starts from 104 × 104 . Moreover, the matrix is sparse and unsymmetrical. Setting the cut-off M /N higher than 10%–20% is not recommended and leads to significant computational errors in solution of the eigenproblem and an exponential increase of required time. Table 1 shows the dependence of the computational time on the resolution of the mapping matrix (N × N) and cut-off of the eigenspectrum (M /N). Below, the compact eigenmode decomposition will be utilized for the investigation of advective–diffusive transport in the SHM for typical forcing protocols. Here in particular will be examined the role of Pe on performance of the compact method and reduction of computational costs.

(0 < x ≤ a) and left (−a ≤ x ≤ 0) parts are filled by a fluid with concentration C = 1 and C = 0, respectively, which is a typical inlet condition in practice. Mapping occurs between inlets of consecutive mixing segments and to this end the mixer crosssection is subdivided into equilateral cells {Ωij }, (i = 1, 235; j = 1, 100), each uniformly filled by 100 particles, yielding an eigenspace of N = 23 500 modes. This seeding density is sufficient for reliable progression of the scalar field by the mapping matrix. Advective–diffusive transport in the high-Pe regime can be reliably approximated by the purely advective mapping method by taking advantage of the numerical diffusion induced by coarsegraining (Section 4.1). The effective Péclet number Pe = Peeff can be determined from relation

5. Eigenmode analysis of advective–diffusive transport

Peeff ≃ p0 /1x2 ,

5.1. Evolution of the scalar field by the compact mapping method Investigation of the scalar evolution by means of the compact mapping method is performed for two typical situations: partially chaotic (Π = 0.2/2.35) and fully chaotic (Π = 1/3), introduced in Section 4.1. For brevity, application of the compact method to other SHM protocols (varying parameter Π ) is not considered here. However, important to note is that performance is essentially the same in all cases investigated. Considered is an initial distribution at the inlet of the micro-channel as shown in Fig. 6. Its right

Fig. 6. Initial scalar distribution at the inlet of the SHM.

(6)

with 1x the mesh size and p0 a constant prefactor [44]. For the considered protocols in the SHM this yields Peeff ≈ 1.5 × 106 (Π = 0.2/2.35) and Peeff ≈ 107 (Π = 1/3) [46]. The performance of the compact mapping method is examined via comparison with an evolution by the full mapping matrix. The error is quantified by the standard deviation

ε=

  N   (Ci − Ci∗ )2  i =1 N

,

(7)

6

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

a

b

Fig. 7. Standard deviation ε of the scalar distribution Cn∗ obtained by the compact mapping for reductions M /N = 5% (‘’) and M /N = 10% (‘△’) for the purely advective case (Pe = Peeff ): (a) protocol Π = 0.2/2.35; (b) protocol Π = 1/3.

with Ci∗ and Ci the scalar concentration in cell i according to compact and full mapping, respectively. Fig. 7 shows the standard deviation ε for Π = 0.2/2.35 (panel (a)) and Π = 1/3 (panel (b)) for truncations M /N = 5% (‘’) and M /N = 10% (‘△’). This reveals that, similar to the idealized flow in [33], the compact mapping method performs well also for the realistic case of the SHM. The approach obtains an acceptable degree of approximation of the scalar field (ε ∼ 10−2 ) for n & 10 – and thus of the outflow mixing quality of a practical SHM – with just 5%–10% of the computational effort (for given eigenmode decomposition) relative to the original approach. Important to note is that the above concerns the high-Pe regime. Introduction of significant diffusion has a considerable impact upon the spectral properties of the mapping matrix and, consequently, on the performance of the compact method. Diffusion augments the decay of the eigenvalues and, inherently, accelerates the exponential decay of the scalar field towards its homogeneous state (Fig. 5). This admits an ever stronger reduction of the eigenmode decomposition for given tolerance with decreasing Pe without compromising accuracy. Fig. 8 shows the evolution of the scalar field Cn∗ obtained by the compact mapping method for reductions M /N = 5% (center column) and M /N = 0.5% (right column) in contrast with the full mapping method (left column) at Pe = 105 . This reveals that even a reduction by a factor of 200 (right column) yields a reliable prediction of the scalar distribution after just 7 periods. Fig. 9 gives the decline of the error ε with the number of periods n for protocols Π = 0.2/2.35 (panel (a)) and Π = 1/3 (panel (b)) for various degrees of reduction M /N at Pe = 105 . This demonstrates that, depending on the flow protocol, an accuracy ε ∼ 10−2 − 10−3 is achievable within n ≈ 7–10 periods for M /N ≈ 0.5%–1%. This accuracy is comparable to the high-Pe cases in Fig. 7 for M /N = 10%. However, considerable reduction of the required sets of eigenmodes (from 10% to 0.5%–1%) significantly lowers the computational effort. The scalar evolution during the initial stages of the mixing process is for reduction M /N = 5% in fact already after n ≈ 2 periods approximated to within an error margin of ε ∼ 10−2 . This shows that a certain degree of accuracy after a given number of periods (in practice e.g. determined by the mixer design) can in principle always be achieved by a proper choice of ratio M /N. The correlation between error ε and the spectral properties of the (reduced) mapping matrix is explored further in Section 5.2. The compact mapping method has a limitation by introducing a residual error ε∞ that sets the lower bound for the maximum achievable degree of approximation. This residual error emerges in Fig. 9 as the plateau upon which ε converges for reduction M /N = 5% (‘’) after n ≈ 8 periods, implying ε∞ ∼ 10−7 as the highest

attainable precision in this case. The other cases exhibit essentially the same behavior with similar values for ε∞ and differ only in number of periods after which the lower bound is reached. Decreasing the Péclet number to Pe ∼ 104 considerably affects the correlation between accuracy ε and reduction M /N. The above trend of a declining number of relevant eigenmodes with stronger diffusion continues, allowing a further decrease in required spectral cut-off M for attaining a given accuracy. Fig. 10 clearly demonstrates this by way of the error ε versus the number of periods n for protocol Π = 0.2/2.35 (panel (a)) at Pe = 3 × 104 and protocol Π = 1/3 (panel (b)) at Pe = 104 . This exposes evolutions qualitatively similar to those shown in Fig. 9 yet at substantially lower ratios M /N, implying that comparable degrees of accuracy are reached with substantially smaller subsets of eigenmodes (M /N = 1%–2%). The residual error nonetheless is of the same magnitude: ε∞ ∼ 10−7 . Fig. 11 gives the evolution of the initial distribution following Fig. 6 obtained by the compact method for reductions M /N = 1% (center) and M /N = 0.3% (right) in comparison to that predicted by the full mapping matrix (left). This conveys essentially the same message as for Fig. 8 in that the compact method reliably captures the evolution. The distributions become indistinguishable after n = 5 periods, despite discarding respectively 99% and 99.7% of the spectrum in the compact approaches. Reductions can be even stronger upon further lowering the Péclet number. Fig. 12 shows the progression of the scalar field according to the compact method at Pe = 5 × 103 for ratios M /N = 0.5% and M /N = 0.05%, revealing that, despite using only extremely small subsets of the spectrum, departures from predictions by the full matrix quickly vanish in the course of the evolution. (Ratio M /N = 0.05% for the present mapping grid encompasses a reduced eigenspace spanning only M = 11 eigenmodes.) Fig. 13 shows the error ε for Pe = 5 × 103 as a function of the reduction M /N. For reduction M /N = 0.5% an approximation error of ε ∼ 10−2 is obtained already within n ≈ 1–2 periods; the highest accuracy ε∞ ∼ 10−7 , in turn, requires a mere n ≈ 7 periods to attain. 5.2. Scaling of the truncation error Key to reliability of the compact mapping method is a truncation error that scales as ϵn ∼ O (|λM +1 |n ) (Section 4.3). Validity of this crucial assumption has been established by Gorodetskyi et al. [33] for advective transport in a 2D model flow. Hereafter this scaling will be examined specifically for advective–diffusive transport in the SHM. Fig. 14 gives three initial conditions C0 at the inlet of the SHM typical of practical (industrial) systems that serve as test cases for

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

7

Fig. 8. Scalar evolution with period n for protocol Π = 0.2/2.35 at Pe = 105 : full mapping versus compact mapping with reductions M /N = 5% and M /N = 0.5%.

a

b

Fig. 9. Standard deviation ε of the scalar distribution Cn∗ obtained by the compact mapping for reductions M /N = 5% (‘’), M /N = 2% (‘◦’), M /N = 1% (‘’) and M /N = 0.5% (‘N’) at Pe = 105 : (a) protocol Π = 0.2/2.35; (b) protocol Π = 1/3.

8

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

a

b

Fig. 10. Standard deviation ε of the scalar distribution Cn∗ obtained by the compact mapping for reductions M /N = 2% (‘◦’), M /N = 1% (‘’), M /N = 0.5% (‘N’), M /N = 0.2% (‘•’) and M /N = 0.1% (‘’): (a) protocol Π = 0.2/2.35 at Pe = 3 × 104 ; (b) protocol Π = 1/3 at Pe = 104 .

Fig. 11. Scalar evolution with period n for protocol Π = 0.2/2.35 at Pe = 3 × 104 : full mapping versus compact mapping with reductions M /N = 1% and M /N = 0.3%.

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

9

Fig. 12. Evolution of the concentration field obtained by full and compact mapping for M /N = 0.5% and M /N = 0.05% of the eigenspectrum. Eigenspace is associated with diffusive mapping matrix for Pe = 5 × 103 in the case of SHM protocol Π = 0.2/2.35.

Fig. 13. Standard deviation ε of the scalar distribution Cn∗ obtained by the compact mapping for protocol Π = 0.2/2.35 at Pe = 5 × 103 . Reductions M /N = 0.5% (‘N’); M /N = 0.2% (‘•’); M /N = 0.1%; (‘’); M /N = 0.05% (‘+’).

a

b

evaluating the scaling of the truncation error. For brevity, results on other initial conditions C0 are not reported here. However, the analysis of different cases revealed basically the same behavior as that discussed below. The evolution of the error ε according to (7) is shown in Fig. 15 for the partially chaotic case Π = 0.2/2.35 as a function of reduction M /N (indicated by markers) and Péclet number Pe for each of the initial conditions (columns correspond with the panels in Fig. 14). Solid lines represent the estimated truncation error ϵn = |λM +1 |n . This exposes a close agreement between the computed (ε ) and estimated (ϵn ) truncations in all cases. Similar correspondence is found for other forcing protocols of the SHM. This validates the essential assumption that the truncation error scales universally with the eigenvalue λM +1 of the first discarded mode M + 1. This has the important practical consequence that the error involved with truncation can indeed be estimated a priori on the basis of the spectrum of the mapping matrix. Conversely, this offers a way to systematically (and reliably) select the appropriate spectral cut-off M for a desired degree of accuracy (Section 4.3).

c

Fig. 14. Initial scalar distributions at the inlet of the SHM.

10

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

Fig. 15. Computed truncation error of the compact mapping method versus prediction ϵn = |λM +1 |n (solid line) for protocol Π = 0.2/2.35 versus Pe and reduction M /N. Reductions M /N = 10% (‘△’); M /N = 5% (‘’); M /N = 2% (‘◦’); M /N = 1% (‘’); M /N = 0.5% (‘N’); M /N = 0.3% (‘’); M /N = 0.2% (‘•’). Columns correspond with initial conditions according the panels in Fig. 14.

6. Conclusions The present study concerns an efficient analysis of advective–diffusive transport in periodic flow systems in terms of eigenmodes of the transport operator. To this end the latter is discretized by way of a compact version of the diffusive mapping method. The diffusive mapping method and compact mapping method are two recent extensions of the conventional mapping method, and present work combines both extensions into the compact diffusive mapping method. The new approach is demonstrated and investigated for advective–diffusive transport in a representative three-dimensional (3D) prototype micro-mixer, the staggered herringbone micromixer (SHM). Key to the compact mapping method is the representation of the scalar evolution by only a small subset of the eigenmodes of the mapping matrix (termed ‘‘transient eigenmodes’’) that dominate the evolution towards the homogeneous asymptotic state. This is of particular relevance for practical systems, which typically are in transient conditions for a significant part of the process duration. To this end the spectrum is cut off such that the resulting truncation error is acceptable. This approach has already proven promising for advective transport studies and admits reductions

of the eigenmode basis by up to a factor 10 without compromising reliability. The spectral gaps in the eigenvalue spectrum of the mapping matrix widen progressively with increasing diffusion. This has the important (practical) implication that diffusion (for a given tolerance) facilitates substantially larger reductions of the eigenmode basis compared to the purely advective limit. This allows for the significant reduction of computational costs for the eigenproblem. Simulations of scalar evolutions in the SHM for various degrees of diffusion reveal in all cases a close agreement between the compact mapping method and the full method. Departures are largest during the initial stage of the transient yet quickly diminish to within acceptable tolerances after only a small number of periods. Important particularly for practical applications is that, despite possible initial differences in scalar evolutions, mixing patterns and mixing quality at the device outlet predicted by both compact and full mapping method in general are indistinguishable. However, the compact method significantly deepens the understanding of the transport phenomena and reliably captures the relevant physics through retention of only a fraction of the full set of eigenmodes (reductions up to a factor 2000). This clearly demonstrates the strong potential of the compact mapping method for an

O. Gorodetskyi et al. / European Journal of Mechanics B/Fluids 49 (2015) 1–11

in-depth (spectral) analysis of advective–diffusive mixing as well as for the systematic optimization and control of the mixing process. The validity and universality of the presumed correlation between spectral cut-off and truncation error have been established. This has the important practical consequence that the cut-off can a priori be chosen such that the truncation error remains within a preset tolerance. The diffusive mapping method is in its current form restricted to advective–diffusive transport subject to no-flux (or adiabatic) boundary conditions. Generalization to other conditions can further extend its reach to systems that involve transport across boundaries (e.g. heat transfer, osmotic processes). This may make eigenmode analyses accessible to an even wider range of practical systems. Efforts to this end are underway. References [1] D.G. Baird, D.I. Collias, Polymer Processing Principles and Design, J. Wiley Inc., New York, 1998. [2] J. Kukura, P.C. Arratia, E.S. Szalai, Understanding pharmaceutical flows, Pharmac. Techn. 26 (2002) 48–60. [3] H.A. Stone, A.D. Stroock, A. Ajdari, Engineering flows in small devices: microfluidics toward a lab-on-a-chip, Annu. Rev. Fluid Mech. 36 (2004) 381–411. [4] M.A. Stremler, F.R. Haselton, H. Aref, Designing for chaos: applications of chaotic advection at the microscale, Philos. Trans. R. Soc. Lond. 362 (2004) 1019–1036. [5] V.V. Meleshko, O.S. Galaktionov, G.W.M. Peters, H.E.H. Meijer, Threedimensional mixing in Stokes flow: the partitioned pipe mixer problem revisited, Eur. J. Mech. B Fluids 18 (1999) 783–792. [6] T.S. Krasnopolskaya, V.V. Meleshko, G.W.M. Peters, H.E.H. Meijer, Mixing in Stokes flow in an annular wedge cavity, Eur. J. Mech. B Fluids 18 (1999) 793–822. [7] H. Aref, Stirring by chaotic advection, J. Fluid Mech. 143 (1984) 1–21. [8] D.V. Khakhar, D.V. Rising, J.M. Ottino, Analysis of chaotic mixing in two model systems, J. Fluid Mech. 172 (1986) 419–451. [9] J.H.E. Cartwright, M. Feingold, O. Piro, Chaotic advection in three-dimensional unsteady incompressible laminar flow, J. Fluid Mech. 316 (1996) 259–284. [10] H. Aref, The development of chaotic advection, Phys. Fluids 14 (2002) 1315–1325. [11] I. Mezíc, Chaotic advection in bounded Navier–Stokes flows, J. Fluid Mech. 431 (2001) 347–370. [12] F. Carlsson, M. Sen, L. Lofdahl, Fluid mixing induced by vibrating walls, Eur. J. Mech. B Fluids 24 (2005) 366–378. [13] S. Childress, A.D. Gilbert, Stretch, Twist and Fold: The Fast Dynamo, SpringerVerlag, Berlin, 1995. [14] J.M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport, Cambridge University Press, Cambridge, 1989. [15] C. Castelain, A. Mokrani, Y. Le Guer, H. Peerhossaini, Experimental study of chaotic advection regime in a twisted duct flow, Eur. J. Mech. B Fluids 20 (2001) 205–232. [16] F. Raynal, J.-N. Gence, Energy saving in chaotic laminar mixing, Int. J. Heat Mass Transfer 40 (1997) 3267–3273. [17] I.C. Christov, J.M. Ottino, R.M. Lueptow, From streamline jumping to strange eigenmodes: bridging the Lagrangian and Eulerian pictures of the kinematics of mixing in granular flows, Phys. Fluids 23 (2011) 103302. [18] C.P. Schlick, I.C. Christov, P.B. Umbanhowar, J.M. Ottino, R.M. Lueptow, A mapping method for distributive mixing with diffusion: interplay between chaos and diffusion in time-periodic sine flow, Phys. Fluids 25 (2013) 052102. [19] S. Cerbelli, V. Vitacolonna, A. Adrover, M. Giona, Eigenvalue-eigenfunction analysis of infinitely fast reactions and micromixing regimes in regular and chaotic bounded flows, Chem. Eng. Sci. 59 (2004) 2125–2144. [20] O.V. Popovych, A. Pikovsky, B. Eckhardt, Abnormal mixing of passive scalars in chaotic flows, Phys. Rev. E 75 (2007) 036308. [21] G. Metcalfe, M.F.M. Speetjens, D.R. Lester, H.J.H. Clercx, Beyond passive: chaotic transport in stirred fluids, Adv. Appl. Mech. 45 (2011) 109–188. [22] R.T. Pierrehumbert, Tracer microstructure in the large-eddy dominated regime, Chaos Solitons Fractals 4 (6) (1994) 774–782.

11

[23] D. Rothstein, E. Henry, J.P. Gollub, Persistent patterns in transient chaotic fluid mixing, Nature 401 (1999) 770–772. [24] G.A. Voth, T.C. Saint, G. Dobler, J.P. Gollub, Mixing rates and symmetry breaking in two-dimensional chaotic flow, Phys. Fluids 15 (9) (2003) 2560. [25] A. Adrover, S. Cerbelli, M. Giona, A spectral approach to reaction/diffusion kinetics in chaotic flows, Comput. Chem. Eng. 26 (2002) 125–139. [26] W. Liu, G. Haller, Strange eigenmodes and decay of variance in the mixing of diffusive tracers, Physica D 188 (2004) 1–39. [27] D.R. Lester, M. Rudman, G. Metcalfe, H.M. Blackburn, Global parametric solutions of scalar transport, J. Comput. Phys. 227 (6) (2008) 3032–3057. [28] J. Sukhatme, R.T. Pierrehumbert, Decay of passive scalars under the action of single scale smooth velocity fields in bounded two-dimensional domains: from non-self-similar probability distribution functions to selfsimilar eigenmodes, Chaos Solitons Fractals 66 (2002) 056302. [29] M. Giona, S. Cerbelli, V. Vitacolonna, Universality and imaginary potentials in advection–diffusion equations in closed flows, J. Fluid Mech. 513 (2004) 221–237. [30] E. Gouillart, O. Dauchot, B. Dubrulle, S. Roux, J.-L. Thiffeault, Slow decay of concentration variance due to no-slip walls in chaotic mixing, Phys. Rev. E 78 (2008) 026211. [31] D.R. Lester, M. Rudman, G. Metcalfe, Low Reynolds number scalar transport enhancement in viscous and non-Newtonian fluids, Int. J. Heat Mass Transfer 52 (2009) 655–664. [32] D.R. Lester, M. Rudman, G. Metcalfe, M.G. Trefry, A. Ord, B. Hobbs, Scalar dispersion in a periodically reoriented potential flow: acceleration via Lagrangian chaos, Phys. Rev. E 81 (2010) 046319. [33] O. Gorodetskyi, M.F.M. Speetjens, P.D. Anderson, An efficient approach for eigenmode analysis of transient distributive mixing by the mapping method, Phys. Fluids 24 (5) (2012) 053602. [34] P.D. Anderson, H.E.H. Mejer, Chaotic mixing analyses by distribution matrices, Appl. Rheol. 10 (2000) 119–133. [35] O.S. Galaktionov, P.D. Anderson, G.W.M. Peters, H.E.H. Meijer, Mapping approach for 3D laminar mixing simulation: application to industrial flows, Internat. J. Numer. Methods Fluids 40 (1–2) (2002) 189–196. [36] M.K. Singh, P.D. Anderson, H.E.H. Meijer, Understanding and optimizing the SMX static mixer, Macromol. Rapid Comm. 30 (2009) 362–376. [37] H.E.H. Meijer, M.K. Singh, P.D. Anderson, On the performance of static mixers: a quantitative comparison, Prog. Polym. Sci. 37 (2012) 1333–1349. [38] M.K. Singh, P.D. Anderson, M.F.M. Speetjens, H.E.H. Meijer, Optimizing the rotated arc mixer, AIChE J. 54 (2008) 2809–2822. [39] M.K. Singh, T.G. Kang, H.E.H. Meijer, P.D. Anderson, The mapping method as a toolbox to analyze, design and optimize micromixers, Microfluid. Nanofluid. 5 (2008) 313–325. [40] A. Sarhangi Fard, M.A. Hulsen, H.E.H. Meijer, N.M. Famili, P.D. Anderson, Adaptive non-conformal mesh refinement and extended finite element method for viscous flow inside complex moving geometries, Internat. J. Numer. Methods Fluids 68 (8) (2012) 1031–1052. [41] T.G. Kang, M.K. Singh, T.H. Kwon, P.D. Anderson, Chaotic mixing using periodic and aperiodic sequences of mixing protocols in a micromixer, Microfluid. Nanofluid. 4 (2008) 589–599. [42] M.K. Singh, M.F.M. Speetjens, P.D. Anderson, Eigenmode analysis of scalar transport in distributive mixing, Phys. Fluids 21 (2009) 093601. [43] M.F.M. Speetjens, M. Lauret, H. Nijmeijer, P.D. Anderson, Footprints of Lagrangian flow structures in Eulerian concentration distributions in periodic mixing flows, Physica D 250 (2013) 20–33. [44] O. Gorodetskyi, M. Giona, P.D. Anderson, Spectral analysis of mixing in chaotic flows via the mapping matrix formalism—inclusion of molecular diffusion and quantitative eigenvalue estimate in the purely convective limit, Phys. Fluids 24 (2012) 073603. [45] O. Gorodetskyi, M. Giona, M.F.M. Speetjens, P.D. Anderson, Analysis of the advection–diffusion mixing by the mapping method formalism in 3D openflow devices, AIChE J. 60 (2014) 387–407. [46] O. Gorodetskyi, M.F.M. Speetjens, P.D. Anderson, Simulation and eigenmode analysis of advective–diffusive transport in micromixers by the diffusive mapping method, Chem. Eng. Sci. 107 (2014) 30–46. [47] A.D. Stroock, S.K.W. Dertinger, A. Ajdari, I. Mezić, H.A. Stone, G.M. Whitesides, Chaotic mixer for microchannels, Science 295 (2002) 647–651. [48] A.D. Stroock, G.J. McGraw, Investigation of the staggered herringbone mixer with a simple analytical model, Philos. Trans. R. Soc. Lond. 362 (2004) 971–986. [49] A. Lasota, M.C. Mackey, Chaos, Fractals, and Noise: Stochastic Aspects of Dynamics, Springer-Verlag, Berlin, 1994. [50] V.V. Meleshko, Steady Stokes flow in a rectangular cavity, Proc. R. Soc. Lond. 452 (1996) 1999–2022.