Eigenmode analysis of advective-diffusive mixing

(1)

Eigenmode analysis of advective-diffusive mixing

Citation for published version (APA):

Gorodetskyi, O. (2014). Eigenmode analysis of advective-diffusive mixing. Technische Universiteit Eindhoven. https://doi.org/10.6100/IR761739

DOI:

10.6100/IR761739

Document status and date: Published: 01/01/2014 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

(2)

Eigenmode analysis of advective-diffusive

mixing

(3)

A catalogue record is available from the Eindhoven University of Technology Library. ISBN: 978-90-386-3519-4

Subject heading: Fluid Mixing, Mapping Method, Molecular Diffusion, Spectral Anal-ysis, Numerical Simulations.

This thesis was prepared with the LA_{TEX 2ε documentation system.} Cover design: Adobe InDesign, O.Gorodetskyi.

(4)

Eigenmode analysis of advective-diffusive

mixing

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Eindhoven, op gezag van de rector magnificus, prof.dr.ir. C.J. van Duijn, voor een

commissie aangewezen door het College voor Promoties in het openbaar te verdedigen op donderdag 16 januari 2014 om 16.00 uur

door

Oleksandr Gorodetskyi

(5)

voorzitter: prof.dr L.P.H. de Goey

1e promotor: prof.dr.ir. P.D. Anderson

copromotor: dr.ir. M.F.M. Speetjens

leden: prof.dr. M. Giona (Universit`a di Roma)

prof.dr. R. Mauri (Universit`a di Pisa) prof.dr.ir. F.N. van de Vosse

(6)

Summary

Mixing of fluids is highly important and widely spread phenomenon, ubiquitously occurring both in natural and artificially created processes. Fluid mixing impacts on a variety of different fields, covering science, industry, everyday life, etc. In many cases, mixing is associated with turbulent motion, for which turbulent fluctuations enhance fluid mixing. However, mixing of highly viscous fluids, usually associated with polymers and pharmaceutical processing, or microfluidic mixing, observed in biotechnology, analytical chemistry, micro-reactors, represent another class of prob-lems where mixing occurs in purely laminar conditions. Laminar chaotic mixing is a powerful alternative to turbulent flows. The general mechanisms that make laminar mixing effective are the repetitive stretching and folding of material elements leading to a reduction of the characteristic length scales, and to an increase of the interma-terial contact area. This leads to the formation of zones with small-scale structures, in which molecular diffusion becomes highly efficient and drives the homogenization process. Direct analysis of the advection-diffusion problem is computationally expen-sive, the development of efficient numerical algorithms for simulating chaotic mixing, especially in the presence of diffusion, is a central issue in fluid mixing theory. The analysis of fluid mixing related to process industry applications, where the fi-nal homogenized mixture is the outcome of mixing units that operates either using stirring parts or using complex static geometrical structures (static mixers), provokes the highlighted interest. The mapping method, based on the distribution matrices represented on a finite grid, offers an efficient and computationally feasible way to investigate the kinematics of mixing for generic mixing protocols in time and spatial periodic flows. Classical mapping approach utilizes coarse-grained discretization of the action of a convective mixing field onto a discretized concentration. The method describes the pure kinematics, thus neglecting the influence of molecular diffusion. However, molecular diffusion and its interplay with advection is often essential to the transport properties. Present thesis proposes an extension of the mapping matrix approach to include the effect of molecular diffusion. The method chosen is based

(11)

on the stochastic formulation of advection-diffusion kinematics (Langevin equations) and on the analogy between stochastic differential equations driven by vector-valued Wiener processes and advection-diffusion dynamics. The stochastic approach is of straightforward implementation and can be regarded as the natural extension of the classical mapping approach. This greatly expands the area of application of the map-ping technique, since diffusion is often significant.

Extension of the mapping approach allows in a quantitative new way to interpret the purely advective mapping method. Coarse-graining acts essentially as an additional superimposed diffusional contribution, and purely advective analysis, performed using the advective mapping matrix, provides a reliable approximation for the advection-diffusion operator, essentially as a consequence of numerical advection-diffusion. The impact of numerical diffusion is associated with some well defined value of the Péclet number, namely effective Péclet number, that is a function of the linear lattice size of the computational grid. The quantitative applicability of the purely advective mapping analysis in modeling of the advection-diffusion processes has been verified comparing spectral data derived from purely advective mapping matrix and Fourier-series expan-sion of the continuum advection-diffuexpan-sion operator. Comparative analysis is applied to prototypical Time Periodic Sine Flow (TPSF) and more complex and realistic Partitioned Pipe Mixer (PPM). Numerical results have shown that the dominant eigenvalues deriving from purely advective mapping simulations and corresponding eigenfunctions provide an accurate estimate of the corresponding quantities associ-ated with the continuum advection-diffusion operator at a value of the Péclet number equal to the effective Péclet number. In this way, it is possible to simulate mixing process at extremely high values of Péclet number (P e ≃ 107_{− 10}9_{), just increasing} grid resolution and performing purely advective (non-stochastic) kinematics. This is a major practical result compared to the existing numerical methods (finite elements, fi-nite volumes or Galerkin expansions), the range of validity of which in mixing analysis over realistic computational structures hardly exceed P e = 106_{. Thereby, mapping} method covers the whole range of Péclet numbers, at lower P e diffusive mapping method is applied, and at high values of Péclet numbers numerical diffusion associ-ated with effective Péclet number is used to simulate advective-diffusive mixing. Diffusive mapping method was implemented to investigate mixing for variety of proto-typical 2D-3D mixing devices at macro and micro levels. Analysis of TPSF and PPM has shown an excellent agreement with results obtained by analysis of advection-diffusion operator using Galerkin expansions, thereby proving validity of the dif-fusive mapping technique. Using extendend mapping approach, advective-difdif-fusive mixing was studied for microflows. The analysis is adopted to three-dimensional Staggered Herringbone Micro-mixer (SHM). Simulations with the diffusive mapping matrix achieve a close agreement with experimental observations from the literature. Finally, method was applied to study mixing performance for industrial-orientated realistic 3D mixer

s that possess complex inner geometries. Mapping method is appliable to mixers with complex geometries, simply adapting reflecting boundary conditions on the solid walls of the mixer. As prototypes were used well-known Twin-Screw Extruder and recently introduced Reacom. The flow field is induced by the rotation of the moving inner

(12)

vii

parts and computed by XFEM. Starting from simple 2D flows and extending study to a complex realistic 3D systems, we have shown that mapping method is an effi-cient numerical tool for understanding the interplay between advection and diffusion, reliably approximating realistic mixing process.

Spectral (eigenvalue/eigenfunction) analysis has found a wide application in investi-gation of distributive mixing. Governing advection-diffusion equation is linear and thus amenables characterization in terms of the eigenvalues and eigenvectors of the evolution operator. The physics of the homogenization process can be easily in-terpreted both as regards of homogenization time- and length-scales, and the spa-tial patterns associated with parspa-tially mixed structures, by considering the eigen-values/eigenfunctions of the mapping matrix. Complex dynamics of the flows may assume existence of invariant sets that can be explored by considering the Poincaré sections of the Lagrangian dynamics associated with the kinematics of fluid elements. However, Poincaré sections do not take into account influence of the molecular diffu-sion. Spectral analysis of the diffusive mapping matrix is an efficient tool for inves-tigating the fine structures associated with advective-diffusive transport. In point of fact, from the structure of the dominant eigenfunctions and of the higher-order modes, coupled with the estimate of their mutual spectral gaps, it is possible to achieve a very detailed qualitative understanding of the dynamics of a mixing process and of the geometry of the resulting partially mixed structures. This qualitative picture, deriving from the analysis of the eigenfunction profiles, is in some sense more general, and physically more significant, than the inspection of the pure diffusionless Poincaré section.

The mapping matrix reveals detailed properties of the distributive mixing and all relevant information is stored in its eigenspectrum. Any vector that describes a dis-tribution of concentration can be expanded in the complete system of linearly inde-pendent eigenvectors of the mapping matrix. The rapid decay of the contribution of each mode in the eigenmode decomposition allows for a truncation of the eigenmode expansion from the whole spectrum to only the dominant eigenmodes, characterized by a decay rate significantly lower than the duration of the mixing process. This truncated decomposition adequately represents the distribution of concentration in-side the flow domain already after a low number of periods, while contributions of all non-dominant eigenmodes rapidly become insignificant. The truncation is deter-mined independently of the initial distribution of concentration and based on the decay rates of the eigenmodes. The compact mapping method is demonstrated for TPSF and SHM. It reveals a reliable prediction of (transient) scalar evolutions and mixing patterns with reductions of the eigenmode basis by up to a factor 2000. The thesis has resulted in the following publications:

• O. Gorodetskyi, M. Giona, P. D. Anderson, Spectral analysis of mixing in chaotic flows via the mapping matrix formalism - Inclusion of molecular dif-fusion and quantitative eigenvalue estimate in the purely convective limit, Phys. Fluids, 24, 073603, (2012)

(13)

study transport and chaotic mixing for extremely large P´eclet values, Europhys. Lett., 97, 14002, (2012)

• 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), 053602, (2012)

• O. Gorodetskyi, M. Giona, M. F. M. Speetjens, P. D. Anderson, Analysis of the advection-diffusion mixing by the mapping method formalism in 3D open-flow devices, AIChE J., (accepted), (2013)

• O. Gorodetskyi, M. F. M. Speetjens, P. D. Anderson, Simulation and eigenmode analysis of advective-diffusive transport in micro-mixers by the diffusive mapping method, Chem. Eng. Sci., (accepted), (2013)

• O. Gorodetskyi, M. F. M. Speetjens, P. D. Anderson, Eigenmode analysis of the advective-diffusive mixing by means of mapping matrix, Eur. J. Mech. B/Fluids, (submitted), (2014)

• O. Gorodetskyi, I.Vivat, M. F. M. Speetjens, P. D. Anderson, Advective-diffusive transport inside complex mixing devices, Chem. Eng. Sci., (submitted) , (2014) The author contributed to other publications:

• P. E. Neerincx, O. Gorodetskyi, H. E. H. Meijer, One step creation of fractal hierarchical structures, Science, (in preparation), (2014)

• I. Vivat, O. Gorodetskyi, P. D. Anderson, Distributive mixing analysis inside complex geometries. Co-rotating mixers., Macromol. Theory Simul., (in prepa-ration), (2014)

(14)

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. The transfer of chemical substances within oceans and in the atmosphere [1] is controlled by mixing and dispersion. Mixing plays also an important role in the maintenance of flora and fauna [2]. The prediction and analysis of mixing and transport patterns are also relevant in geophysical, oceanographic, and atmospheric flows. Further examples include population genetics, distribution of nutrients, oil spills and the prediction of the weather [1, 3–7].

Fluid mixing can be categorized according to the flow and mixing conditions or the fluid properties. In many cases, mixing is associated with turbulent motion, for which turbulent fluctuations enhance fluid mixing [8]. However, many industrial applications concern mixing under purely laminar flow conditions. This regime is characterized by a low Reynolds number Re = U L/ν, 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 for-mer include chemical engineering, cosmetic industry, food processing, polyfor-mer and pharmaceutical processing [9–11]; the latter encompasses the highly-relevant field of micro-fluidics, which is at the heart of many emerging technologies. Here Reynolds numbers often are in the order of unity or less [12–15]. In compounding of various substances it is necessary to obtain fluid mixtures with specified homogeneity prop-erties, such as necessary for solutions, emulsions, suspensions. In such problems, the efficiency of mixing is determined by the repeated stretching and folding of material lines and surfaces due to advective transport of fluid volumes [16, 17]. The stretching and folding mechanism is enhanced if the flow induces Lagrangian chaos in the kine-matics of fluid particles [17, 18]. The mixing efficiency is determined strongly by the stretching rate of material lines and surfaces by the advective transport[16, 17, 19]. As thoroughly addressed by a huge Literature in the last thirty years [18, 20–27], lam-inar chaotic mixing provides a powerful alternative to turbulent flows. Stretching and

(15)

folding cause rapid reduction of the characteristic size of zones filled by uniform com-ponent, thus increasing the intermaterial contact area [28, 29]. Stretching and folding by chaotic advection leads to the formation of zones with small-scale structures, in which molecular diffusion becomes highly efficient and promotes fast homogenization of the mixture. However, despite the essential part of molecular diffusion in the mixing process, it has received far less attention than advection.

Focusing attention on industrial mixers, typically the flow device consists in a set of periodic units, and the same stretching and folding protocol of material elements re-peats in every unit. Direct analysis of the advection-diffusion problem for such devices is computationally expensive, especially when the geometry of the mixer is complex, as for static mixers, mixers with rotating inner parts, etc. There is a significant moti-vation for developing efficient numerical algorithms for investigating chaotic mixing, especially in the presence of diffusion, that can be applied to complex flows, and generic mixing devices. For discussed above mixers (i.e. with laminar flows possess-ing either spatial or time periodicity), Spencer and Wiley [30] proposed a method that describes transport of a conservative quantity from one state to another by means of a discretized mapping matrix Φ. The mapping method [31, 32] offers an efficient way to investigate the kinematics of mixing for generic mixing protocols inside periodic flow systems. A variety of two and three-dimensional prototypes of flows and mixers has been studied by using the mapping method: simple prototypical flows as time peri-odic sine flow [33], cavity flow [34, 35], journal bearing flow [36], Kenics mixer [37, 38] and other industrial inline mixers [39–41], the Rotated Arc Mixer (RAM) [42], va-riety of micromixers [43–45], granular flows [46], various extruders [47, 48], etc. In these works, the evolution of a scalar field is studied by the iterative application of the the mapping matrix to the vector describing the concentration distribution after each period of mixing. From the mathematical point of view, the mapping matrix can be viewed as the coarse-grained representation of the Frobenius-Perron operator [49] associated with the kinematics of fluid elements, i.e., a modification of the Ulam’s approximation scheme[50] for discretizing the Frobenius-Perron operator associated with the kinematics[51, 52]. By applying appropriate mixing measures to these discretized concentration fields, mixing processes could be understood and stirring protocols optimized.

Microfluidics

Microfluidic mixing and scalar transport in flow domains of size O(100µm) is key to a wide range of technologies and applications as biotechnology, analytical chemistry, micro-reactors, micro-total-analysis systems and labs-on-a-chip [53–56]. Design and optimization of efficient micromixers is essential to technological advancement in these fields. This has triggered development of a great variety of micro-mixers during the last decade, encompassing both miniaturized versions of large-scale industrial static mixers [57] and dedicated systems like the three-dimensional (3D) serpentine mixer [58], the barrier-embedded mixer [59], the serpentine laminating micro-mixer [60], T-junction and multilaminated/elongational-flow micro-micro-mixers [61], the staggered herringbone mixer [62], etc. Refer to[63–66] for exhaustive reviews on

(16)

micro-3

mixing technology.

Numerical simulations of fluid flow and scalar transport are an indispensable element in the design, analysis and optimization of micromixers. Many devices operate on the inline-mixing principle: stirring of a continuous axial throughflow by systematic reorientation and periodic repetition of the transverse flow forcing [67]. Aim of the reorientation is achievement of chaotic advection in transverse direction and, in con-sequence, efficient mixing. Thereby, the mapping method is an efficient and flexible computational method for this broad class of micro-mixers. Mapping method has been applied to analyze fluid mixing inside a variety of microfluidic devices such as barrier embedded micromixer (BEM) [68], staggered herringbone mixer (SHM) [44], serpentine mixer [43, 69], slanted ridge mixer (SRM) [45].

Molecular diffusion

There is a major conceptual limitation on the original formulation of the mapping method discussed above, that becomes physically relevant especially when applied to practical industrial devices. It describes the pure kinematics, thus neglecting the influence of molecular diffusion. However, molecular diffusion and its interplay with advection is often essential to the transport properties. Present thesis extends the advective formulation of the mapping matrix approach by including the effects of molecular diffusion in the construction of the matrix. The method is based on the stochastic formulation of advection-diffusion kinematics and on the analogy between stochastic differential equations and advection-diffusion dynamics [49, 70]. The inclu-sion of diffuinclu-sion in the classical kinematic formulation of the mapping method gives a valuable contribution in the investigation of mixing problems, providing a closer cor-respondence between simulations and experiments i.e. the “real” advective-diffusive mixing, where diffusion is the ultimate driving force to homogenization. This allows a much more accurate prediction of the outcome of a mixing process, hereby reducing time and energy costs when real industrial devices have to be analyzed.

Spectral analysis

In the last decades spectral analysis has become popular for the investigation of scalar transport, as the governing advection-diffusion equation is linear and thus amenable to characterization in terms of the eigenvalues and eigenvectors of the evolution operator[71]. Spectral analysis provides the proper formal setting to ex-plain the asymptotic exponential decay of the scalar field (controlled by the eigen-values of the transport operators), and the occurrence of repeated mixing patterns (controlled by the corresponding eigenfunctions), the intensity of which decay expo-nentially as a tracer field approaches the completely mixed state. Such repeating patterns, observed in chaotic flows, were called strange eigenmodes by Pierrehumbert [72]. Strange eigenmodes were observed experimentally by Rothstein [73], and fur-ther analyzed by Liu and Haller [74]. Both Pierrehumbert [72] and Rothstein [73] supported the observation that strange eigenmodes appear also in aperiodic flows as well. Toussaint et al. [75] were among the first that considered spectral properties

(17)

of the advection-diffusion operator in a three-dimensional chaotic flow. Gleeson [76] applied spectral methods to study simple chaotic flows in an annular geometry. The group of Giona at “La Sapienza” performed detailed analysis of the strange eigen-modes of the advection-diffusion operator in two-dimensional laminar chaotic mixing flows [23, 77, 78]. Simple non-chaotic autonomous flows were thoroughly investigated [79, 80], using a Schrödinger-like formulation of the advection-diffusion equation. The above-mentioned studies investigate the advective-diffusive transport in terms of eigenmode decompositions of the advection-diffusion operator. The connection between the spectral properties of the mapping matrix and the advection-diffusion operator has been revealed in Singh et al. [33] in the purely advective limit of the scalar transport. This paper is mainly focused on the eigenfunctions of the purely advective mapping matrix and on their relation to the structure of the quasi-periodic regions of the Poincaré maps associated with the flow. The inclusion of the diffusion into the description of the mixing process gives qualitatively new results [81–84]. Structure of the dominant eigenfunctions and of the higher-order modes allows to achieve a very detailed qualitative understanding of the dynamics of the advective-diffusive mixing process and to investigate influence of the diffusion on the geometry of the resulting partially mixed structures. In some sense, such analysis is more general and physically more significant, than the inspection of the pure diffusionless “frozen” picture arising from the analysis of the kinematic Poincaré sections.

Numerical diffusion

Usually effect of numerical diffusion arises as undesirable numerical error, induced by discrete representation of the concentration field on the finite grid. Present thesis provides a new view on the effect of numerical diffusion and shows that the purely advective mapping matrix can be interpreted in a new quantitative way. Coarse-graining acts essentially as an additional superimposed diffusional contribution, and purely advective mapping method provides a reliable approximation for the advection-diffusion operator, essentially as a consequence of numerical advection-diffusion induced by discretization. Impact of the numerical diffusion in the mixing process is associated at some well defined value of the Péclet number, namely effective Péclet number, that is a function of the computational grid resolution. It allows simulating of the advective-diffusive processes at large Péclet values (P e ≥ 107_{) by simply varying the} size of the computational grid, that is unaffordable by existing numerical methods (finite elements, finite volumes or Galerkin expansions, etc.).

Compact mapping method

Spectral analysis of the mapping matrix reveals detailed properties of the distributive mixing as all relevant information about mixing properties is stored in its eigenmodes. Any vector that describes a distribution of concentration can be expanded in the com-plete system of linearly independent eigenvectors of the mapping matrix. Eigenmode decompositions of the mapping matrix was presented by Singh et al.[33] at the purely advective limit. Due to the rapid decay of the contribution of the major part of

(18)

5

the eigenspectrum, especially in presence of diffusion, eigenmode decomposition al-lows truncation of the eigenmode expansion from the whole spectrum to only the dominant eigenmodes and representation in the compact form. Compact approach adequately represents of the scalar evolution by only a small subset of the eigenmodes already after a low number of periods.

Thesis overview

The thesis is organized as follows. Chapter 2 overviews classical mapping approach and describes extension of the method to the advective-diffusive case. Spectral prop-erties of the mapping matrix, reliability of the mapping approach, introduction and estimation of the numerical diffusion are discussed in Chapter 3 for the case of simple 2D Time Periodic Sine Flow (TPSF). Chapter 4 extends application of the diffusive mapping method to realistic 3D open flows, using the Partitioned Pipe Mixer (PPM) as a prototypical system. Advective-diffusive transport in microchannels is studied in Chapter 5. As an example of the micromixer Staggered Herringbone Mixer (SHM) is considered. Implementation of the compact mapping method to simple 2D TPSF in the purely advective limit is discussed in Chapter 6. Compact approach is extended to advective-diffusive case for 3D SHM in Chapter 7. Chapter 8 investigates mixing inside realistic industrial mixers. Conclusions and recommendations are drawn in the last Chapter 9.

(19)

(20)

2 _{Inclusion of molecular}

diffu-sion in the mapping method

1

Abstract

Present Chapter extends the original mapping matrix formalism to include the effects of molecular diffusion in the analysis of mixing processes. The approach followed is Lagrangian, by considering the stochastic formulation of advection-diffusion processes via the Langevin equation for passive fluid particle motion and on the analogy be-tween stochastic differential equations and advection-diffusion dynamics. Proposed technique provides straightforward extension of the mapping method for simulating advective-diffusive transport problems.

2.1 Introduction

Apart from the fluid dynamic applications, the analysis of chaotic mixing has deep and diversified implications in many other fields of science. Laminar chaotic mixing is a powerful alternative to turbulent flows in polymer processing and in the presence of highly viscous fluids [10], and it represents a condition frequently encountered in microfluidic theory and practice [14]. We are focusing attention on periodic flows, where the same stretching and folding protocol of material repeats after specified period. For such flows (i.e. for flows possessing either spatial or time periodicity), Spencer and Wiley [30] proposed a method that describes transport of a conservative quantity from one state to another by means of a discretized mapping represented, on a finite mesh, by a matrix. The mapping method [31, 32], based on the distri-bution matrices, offers an efficient way to investigate the kinematics of mixing for 1 _{Reproduced from: 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, 073603, (2012)

(21)

generic mixing protocols in time and spatial laminar periodic flows. A variety of two and three-dimensional prototypes of flows and mixers has been studied by using the mapping method [48], including the Kenics and other static mixers [37, 39, 40], the Rotated Arc Mixer (RAM) [42], micromixers [44], various extruders [47], etc. The evolution of a scalar field is studied by the iterative application of the the mapping matrix to the vector describing the initial/inlet concentration distribution. From the mathematical point of view, the mapping matrix can be viewed as the coarse-grained representation of the Frobenius-Perron operator [49] associated with the kinematics of fluid elements, i.e., a modification of the Ulam’s approximation scheme[50] for discretizing the Frobenius-Perron operator associated with the kinematics[51, 52]. There is a major conceptual limitation of the original mapping method described above, that becomes physically relevant especially when applied to practical industrial devices. It describes the pure kinematics, thus neglecting the influence of molecular diffusion. The development of efficient numerical algorithms for simulating chaotic mixing, especially in the presence of diffusion is a central issue in fluid mixing theory.

2.2 Statement of the problem

This Section defines the problem setting and provides a concise overview of the tech-nical tools applied throughout the Chapter.

Let c(x, t) be a scalar field associated with a physical transported entity C (a dif-fusing solute, temperature, etc.) defined in a mixing domain M, and advected by a deterministic laminar incompressible flow V(x, t), ∇ · V = 0. Let D be the diffusiv-ity of C in the fluid continuum. The evolution of C within M is described by the advection-diffusion equation

∂c

∂t = −∇ · (V c) + D∇

2_{c ,} _(2.1)

for its concentration field c. Let cref be a reference value for c, L0 a characteristic length-scale of the mixing domain M, and V0 a characteristic velocity scale. Intro-ducing u(x, t) = c(x, t)/cref, v(x, t) = V(x, t)/V0, and t V0/L0 7→ t, x/L0 7→ x, the non-dimensional form of Equation (2.1) reads as

∂u

∂t = −∇ · (v u) + 1 P e∇

2_{u ,} _(2.2)

where P e = V0L0/D is the P´eclet number. For notational simplicity we still use the symbols t and x to indicate the non-dimensional time and position vector, as well as the nabla operator ∇ entering Equation (2.2) refers to the non-dimensional coordinates. Henceforth, unless otherwise specified, solely non-dimensional variables are used.

2.3 Bounded and unbounded mixing systems

The mixing domain M can be bounded or unbounded. In closed bounded sys-tems, mixing quantification can be obtained by considering the norm decay of the

(22)

2.4. Original mapping method 9

field u(x, t) as a function of time t. Since in a closed and bounded system the normal component of the flux vanishes at the boundary of M, the mean value u(t) = (1/meas(M))RMu(x, t) dx where meas(M) is the Lebesgue measure of M -is conserved, i.e., u(t) = u(t = 0) = u0 for any t > 0. It is convenient to define a new field φ(x, t) = u(x, t) − u0 possessing zero mean. The field φ(x, t) satisfies the same advection-diffusion equation defining the evolution of u(x, t),

∂φ

∂t = −∇ · (v φ) + 1 P e∇

2_{φ .} _(2.3)

In open mixing systems, such as a mixing tube, static mixers, etc., the mixing domain M can be represented as the Cartesian product M = (0, 1) × Σ⊥, where Σ⊥ is the transverse cross section. In many applications of open flow devices steady-state mixing properties are particularly relevant.

Let x = (x⊥, z), where x⊥and z are respectively the transverse and axial coordinates, and v = (v⊥, w) be the analogous decomposition for the non-dimensional velocity field. For high P´eclet numbers, if the characteristic transverse length is much smaller than the axial length of the device, axial diffusion is negligible. Therefore, at steady state Equation (2.2) simplifies as:

w∂u ∂z = −v⊥· ∇⊥u + 1 P e∇ 2 ⊥u , (2.4)

where ∇⊥ is the transverse nabla operator and incompressibility has been enforced. Let us assume w = w(x⊥) to be independent of the axial coordinate (as for Stokes flows in straight channels), and that w(x⊥) ≥ 0. Under these assumptions, Equation (2.4) is similar to the time-evolution of a scalar field in a two-dimensional mixing domain (in this analogy represented by the cross section Σ⊥), where the time vari-able is replaced by the axial coordinate weighted with respect to the axial velocity component.

Let φ = u−uw, where uw(z) = (1/meas(Σ⊥))R_Σ

⊥w(x⊥) u(x⊥, z) dx⊥is the weighted

average of u with respect to the axial coordinate. Since the flux vanishes at the boundary of the cross section, from Equation (2.4) it follows that uw(z) = uw,0 = constant. Therefore, φ(x, z) satisfies Equation (2.4).

2.4 Original mapping method

The mapping method is based on a suggestion due to Spencer and Wiley [30] that describes the transport of a conservative quantity from one state to another by means of a discretized mapping represented, on a finite grid, by a matrix. In their 1951-paper, Spencer and Wiley clearly state their idea: “...Now define a distribution matrix [Dij], where Dij is the fraction of material originating in the i-th cell which is to be forced in the j-th cell after the mixing operation.”

Several decades after this suggestion, a wealth of articles applied to different pro-totypical, industrial, and micro-mixers have shown that the mapping approach (i.e. the coarse-grained discretization of the action of a convective mixing field onto a discretized concentration field) provides an efficient computational tool for obtaining

(23)

an accurate characterization of mixing properties in the purely kinematic case at a feasible computational cost [32, 44, 85]. The entries of this matrix contain the frac-tions of material from one part of the domain which is transferred to various parts of the domain when specific flow is applied. Figure 2.1 depicts how the entries of the mapping matrix are approximated. The flow domain is divided into a number of cells

Figure 2.1: Illustration of the computation of the entries Φijof the mapping matrix Φ. The cell Ωj at t = t0 is covered with a number of markers that are tracked during flow in ∆t to arrive at the final cross section t = t0+ ∆t.The ratio of the number of markers received by the recipient cell Ωi to the initial number of markers in Ωj is determined. In this case Φij= 3/16.

N which are numbered sequentially and N2

p markers are uniformly placed inside each cell and tracked in time/space (period of the flow protocol) from t = t0to t = t0+ △t. If the number of markers in the donor cell Ωj is Mj= Np2 at t = t0 and the number of markers found after tracking in the recipient cell Ωiis Mij at t = t0+ ∆t, then the entry Φij of the mapping matrix Φ is calculated [85] as

Φij = Mij

Mj

. (2.5)

Once the mapping matrix is evaluated over one period, the evolution of the concen-tration field at discrete time/space intervals multiple of the period can be estimated as

Cn= (Φ(Φ(...(Φ

| {z }

n times

C0)...))). (2.6)

where C0is a N -dimensional vector representing the initial distribution of concentra-tion inside the flow domain. The vector Cn provides the coarse-grained description of the volume fraction (dimensionless concentration) of the transported entity in the moving continuum, and its i-th entry corresponds to the concentration locally aver-aged at the cell Ωi. Constructed mapping matrix concerns purely advective transport and does not take into account last term in the Equation (2.2) that corresponds to dif-fusion. The main motivation is inclusion of the molecular diffusion in the formulation of the mapping approach.

(24)

2.5. Diffusive mapping matrix 11

2.5 Diffusive mapping matrix

Consider a d-dimensional flow (d = 2, 3) associated with a time-periodic non-dimensional velocity field v(x, t + Tp) = v(x, t), defined on a bounded domain M.

The Lagrangian motion of a passive fluid particle in a d-dimensional flow domain M, initially located at x0∈ M, is affected by two factors: (i) the advective contribution associated with a given stirring mechanics described by a velocity field v(x, t), and the influence of random fluctuations at a molecular level, and described by the non-dimensional stochastic differential equation[49, 70, 86]

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

2

P edξ(t) , (2.7)

where dξ(t) are the infinitesimal increments of a d-dimensional Wiener process, i.e., of a vector-valued stochastic process ξ(t) = (ξ1(t), . . . , ξd(t)) in which ξh(t) and ξk(t) are independent stochastic processes for h 6= k, possessing uncorrelated increments distributed in a normal way, i.e. hξh(t)i = 0, and h(ξh(t + τ ) − ξh(t))2i = τ, for any h = 1, . . . , d and t, τ > 0, and such that ξ(0) = 0. Equation (2.7) is equipped with the initial condition x(0) = x0. Let P (x, t/x0) be the probability density function (pdf) associated with the stochastic process (2.7). For notational simplicity we indicate this quantity as P (x, t), omitting the reference on the initial condition x0. The pdf P (x, t) satisfies the stochastic differential equation[49],

∂P (x, t)

∂t = −∇ · (v(x, t) P (x, t)) + 1 P e∇

2_{P (x, t) .} _(2.8)

Equation (2.8) corresponds to the classical non-dimensional advection-diffusion equa-tion (2.2), where P e is the P´eclet number.

The stochastic differential equation (2.7) can be solved numerically by means of an Euler-Langevin explicit scheme [49]. Let h be the time step, xn= x(nh), and vn(x) = v(x, nh). The Euler-Langevin solver for Equation (2.7) is simply given by

xn+1= xn+ vnh + r 2 P eh 1/2_r n , (2.9)

where rn = (r1,n, . . . , rd,n) is the realization of a d-dimensional normally distributed independent random variable. This means that the d random variables r1,n, . . . , rd,n, are independent of each other and distributed according to a Gaussian pdf with zero mean and unit variance.

For any finite value of P e, the mapping matrix obtained by means of the method de-scribed above will be indicated as the diffusive mapping matrix. The matrix obtained in the limit case of P e → ∞, i.e., when no stochastic contribution is superimposed to deterministic advection induced by the flow, will be referred to as the purely advective mapping matrix.

2.6 Conclusions

By making use of the stochastic formulation of advection-diffusion processes the map-ping matrix can be constructed by following exactly the same computational strategy

(25)

adopted in the purely advective case. Specifically:

• The mixing domain M is discretized with a structured grid of N = Nd c boxes (Nc× Nc for two-dimensional simulations).

• Within each box, say the i-th, Nd

p particles are considered, uniformly distributed at internal points. Let x(i,j) _{be the position vector of the j-th particle j =} 1, . . . , Nd

p.

• For each x(i,j)_{, let x}(i,j)

0 = Φ−1(x(i,j)) the effective initial position, where Φ−1 is the inverse of the kinematic Poincar´e map.

• The Langevin equation (2.7) is solved according to Equation (2.9) for a time interval Tp, equal to the flow period, starting from the initial condition x(i,j)0 . Let Xf = {x(i,j)f }

Nd p

j=1 the ensemble of the final positions for the (Nc× Np)d particles considered.

• The mapping matrix Φ is evaluated in the usual way i.e. by means of Equa-tion (2.5), where Mij is the number of particles in the ensemble Xf that fall within the i-th box of the discretization, given that their initial position is within the j-th box, and Mj is the number of initial positions falling within the j-th box.

(26)

3 _{Spectral analysis of mixing in}

chaotic flows

2

Abstract

Inclusion of diffusional effects in the mapping matrix formalism permits to frame the spectral properties of mapping matrices in the purely advective limit in a quantitative way. Specifically, the effects of coarse-graining can be quantified by means of an effective P´eclet number that scales as the second power of the linear lattice size. This simple result is sufficient to predict the scaling exponents characterizing the behavior of the eigenvalue spectrum of the advection-diffusion operator in chaotic flows as a function of the P´eclet number, exclusively from purely kinematic data, by varying the grid resolution. Simple but representative model systems is considered under a wealth of different kinematic conditions - from the presence of large quasi-periodic islands intertwined by chaotic regions, to almost globally chaotic conditions, to flows possessing “sticky islands” - providing a fairly comprehensive characterization of the different numerical phenomenologies that may occur in the quantitative analysis of mapping matrices, and ultimately of chaotic mixing processes.

3.1 Introduction

In a physical perspective, fluid mixing is the result of the interaction of a stirring flow with diffusion. Therefore, it is physically relevant to frame spectral analysis in this perspective, by considering the interaction between an advecting field and molecular 2_{Reproduced from: 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, 073603, (2012)

O. Gorodetskyi, M. Giona, P. D. Anderson, Exploiting numerical diffusion to study transport and chaotic mixing for extremely large P´eclet values, Europhys. Lett., 97, 14002, (2012)

(27)

diffusion, i.e. by analyzing the eigenvalue/eigenfunction properties of the advection-diffusion operator. Spectral analysis provides the proper formal setting to explain the asymptotic exponential decay of mixing norms (controlled by the eigenvalues of the transport operators), and the occurrence of repeated mixing patterns (controlled by the corresponding eigenfunctions) the intensity of which decay exponentially as a tracer field approaches the completely mixed state. Such repeating patterns, observed in chaotic flows, were called strange eigenmodes by Pierrehumbert [72].

The advective-diffusive transport in terms of eigenmode decompositions of the advection-diffusion operator has been widely investigated in the literature [23, 75–78]. The connection between the spectral properties of the mapping matrix and the advection-diffusion operator has been revealed in Singh et al.[33] in the purely advective limit of the scalar transport. In this paper, Singh et al. apply the mapping matrix formalism to address the spectral properties of purely advected scalar fields, by focusing mainly on the eigenfunctions of the purely advective mapping matrix and on their symmetries that are related to the structure of the quasi-periodic regions of the Poincar´e maps associated with the flow.

The present Chapter is essentially aimed at a quantitative analysis of the eigenvalues of mapping matrices in order to relate the spectral structure of the mapping matrix to the corresponding properties of the advection-diffusion operator. For this reason, we consider a simple flow systems, namely the time-periodic sine flow, as a paradigmatic example of a bounded chaotic flow. Essentially, the scope of the Chapter is to ana-lyze the numerical diffusion effects induced by the coarse-graining in the discretized representation of the mapping matrix in the purely advective case, where solely the action of the flow field is considered, and in the general case where it superimposes to molecular diffusion.

As a relevant byproduct of this analysis we show that even when the purely advective case is considered, numerical diffusion can be exploited to mimic molecular diffusion, and a strict quantitative correspondence between the spectral structure of the purely advective mapping matrix and that of the continuum advection-diffusion operator can be defined. More precisely, for any coarse-graining resolution, i.e., for any size of the discretized grid defining a given mapping-matrix approximation, an effective Péclet number can be defined (that is associated with numerical diffusion), so that the mapping matrix at this level of resolution can be viewed as a quantitatively reliable approximation of the continuum advection-diffusion operator for a value of the Péclet number equal to the effective Péclet number of the grid resolution (see Section 3.6 for details).

3.2 Problem definition

Consider a closed mixing system in the presence of a time-periodic velocity field v(x, t + Tp) = v(x, t), where Tp > 0 is the flow period. In order to obtain an autonomous (time-independent) representation of the concentration field dynamics, it is convenient to sample its evolution at time instants multiple of the flow period. Therefore, the stroboscopic operator P : L2_{(M) → L}2_{(M) can be defined (from} Equation (2.2)) in such way that it transforms the scalar field un(x) = u(x, nTp)

(28)

3.3. Model system: time-periodic sine flow (TPSF) 15

(solution of Equation (2.2)) at time t = nTpinto the field un+1(x) = u(x, (n + 1)Tp) after one flow period,

un+1(x) = P[un(x)] . (3.1)

The operator P is essentially the Floquet operator associated with an advection-diffusion equation [78].

For any P e < ∞, the operator P is compact, and therefore it admits a countable spectrum of eigenvalues {µn}∞n=1, lying in the unit circle, that can be ordered in a non-increasing way with respect to the moduli, |µ1| > |µ2| ≥ |µ3| ≥ · · · |µn| ≥ · · · . Let {ψn}∞n=1be the corresponding eigenfunctions. Due to conservation of the transported entity, the first eigenvalue equals 1, i.e. |µ1| = 1, and the corresponding eigenfunction is uniform, ψ1 = const. Starting from µ2, all the other eigenvalues admit norms strictly less than 1. Therefore, for generic initial conditions (i.e. apart from a set of initial conditions of zero measure), the norm decay of scalar fields transported according to Equation (2.2) is controlled by the second eigenvalue µ2 (referred to, for this reason, as the dominant eigenvalue). Specifically, since φ0(x) = u0(x) − u0 possesses zero mean, then

||φn||L2 ∼ |µ|n= e−nTpΛ, (3.2)

where the exponent Λ will be referred to as the dominant decay exponent of the time-periodic advection-diffusion operator,

Λ = −log |µ_T 2| p

. (3.3)

Let {λn} be the eigenvalues of mapping matrix Φ, that lie within the unit circle. Similarly, eigenvalues can be ordered in a decreasing way with respect to their moduli. For a closed mixing system, the first eigenvalue equals |λ1| = 1, the second eigenvalue λ2 is less than one, and, within the mapping matrix formalism is related to the dominant decay exponent Λmp, via the relation

Λmp= −log |λ2| Tp

, (3.4)

where Tp is the flow period. Henceforth, we indicate with the subscript “mp” the dominant decay exponent obtained from the mapping matrices, and with Λ the cor-responding exponent obtained from the spectral analysis of the continuous advection-diffusion operator, obtained e.g. by Fourier methods as in Cerbelli et al.[78]

3.3 Model system: time-periodic sine flow (TPSF)

As a benchmark model we consider the time-periodic sine flow (TPSF) [87, 88]. TPSF is defined as periodic sequences of two steady sinusoidal flows, which act during the half of the period Tpinside the unit bi-periodic square Ω = (0, 1) × (0, 1),

v1= (vx, vy) = (sin (2πy), 0) for 0 ≤ t < Tp 2 , v2= (vx, vy) = (0, sin (2πx)) for

Tp

2 ≤ t < Tp.

(29)

The qualitative kinematic properties of TPSF protocols (obtained by varying the pe-riod Tp) are analogous to those observed in physically realizable flow systems. There-fore, the use of TPSF allows to gather valuable information on fluid mixing under a wealth of different kinematic conditions at moderate computational costs. TPSF admits an analytical expression for fluid particle orbits. Its Poincar´e map is given by:

Ψ(xi, yi) =

xi = xi−1+ Tpsin (2πyi−1)/2 mod. 1,

yi= yi−1+ Tpsin (2πxi)/2 mod. 1. (3.6)

(a)

(b)

(c)

(d)

Figure 3.1: Poincar´e sections of the TPSF for different periods Tp. Panel (a) corresponds to Tp= 0.56; panel (b) to Tp= 0.8; panel (c) to Tp= 1.18; panel (d) to Tp= 1.6 Floquet operator P obtained by Fourier-series expansion and its spectral properties has been studied by Cerbelli et al. [78], and we use the data reported in this article (together with new data obtained using the same Fourier expansion) to compare

(30)

3.4. Constructing mapping matrix 17

the spectral properties of the mapping matrix with those of the advection-diffusion operator.

The flow structure and corresponding mixing properties can be made visible by Poincar´e sections [17]. This is demonstrated in Figure 3.1 for the four stirring pro-tocols corresponding to Tp = 0.56, Tp = 0.8, Tp = 1.18 and Tp = 1.6, exposing the characteristic flow structure of 2D mixing flows: islands and chaotic seas. The pro-gression with growing Tp reveals a gradual diminution of islands in favor of a chaotic sea until a state of chaos sets at Tp= 1.6.

3.4 Constructing mapping matrix

3.4.1 Numerical issues

For fixed grid resolution Nc(the effects of the coarse-grained resolution are thoroughly addressed in Section 3.6.1), the construction of the mapping matrix depends on Np and, in the diffusive case, of the time step h.

A preliminary analysis of the effect of the time step h indicates that the choice h = 10−2 _{provides a reliable trade-off between accuracy and computational speed in} the simulation of the Langevin equation through the Euler-Langevin algorithm (2.9). The other numerical parameter is the number of particles Nd

p (Np2in two-dimensional simulations) that are placed within each element of the discretization grid for obtain-ing the mappobtain-ing matrix as discussed in Section 2.4.

Consider first the purely advective case. Figure 3.2 depicts the behavior of the domi-nant decay exponent (the superscript ‘c’ indicates the purely advective case) obtained for the purely kinematic mapping matrix for several typical values of the flow period, and grid resolutions.

101 102 103 10−3 10−2 10−1 100 Np Λ c m p

Figure 3.2: Dominant decay exponent Λc

mp of the purely advective mapping matrix as a

function of the specific particle number Np. Symbols ‘’ refer to Tp = 0.8, Nc = 60, ‘’ to Tp = 0.8 and Nc = 100, ‘◦’ to Tp = 1.6 and Nc = 60, ‘•’ to

Tp = 1.6 and Nc = 100, symbols (▽, N, △) to Tp = 2 at Nc = 60, 80, 100,

(31)

The typical situation, occurring for flow protocols the Poincar´e maps of which possess either large islands of quasi-periodicity (Tp = 0.8 corresponding to symbols ‘’ and ‘’) or seemingly globally chaotic conditions (Tp= 1.6 corresponding to symbols ‘◦’ and ‘•’), is that is sufficient to take Np= 20 (i.e. 400 particles per discretization box) to achieve an accurate estimate of the mapping matrix. However, situations occur in which this may not be sufficient. This is for example the case of the TPSF at Tp= 2.0 (symbols ‘▽’, ‘N’, ‘△’), for which at least Np= 100 is required for convergence. This case corresponds to a Poincar´e map possessing small islands of quasi-periodicity or parabolic periodic points, characterized by long-range correlations [89, 90] in the particle kinematics.

An analogous phenomenon occur when diffusion is accounted for via the stochastic approach described in Section 2.5. Figure 3.3 depicts the behavior of Λmp at Tp= 2, which is the worst scenario among those considered as it regards the influence of Npof the convergence of the mapping-matrix construction algorithm. As can be observed,

101 102 103 10−1 100 Np Λm p

Figure 3.3: Dominant decay exponent Λmp of the diffusive mapping matrix vs Npat Tp= 2 and at different values of the P´eclet number. Markers ‘’ refer to P e = 102_{, ‘△’} to P e = 103_{, ‘♦’ to P e = 10}4_{, ‘◦’ to P e = 10}5_.

Np= 50 is the relevant number of particle per discretization box that are required to reliable and accurate construction of the mapping matrix in the presence of diffusion. 3.4.2 Statistical analysis of the mapping matrix

In the presence of diffusion, the mapping matrix obtained via the stochastic Langevin equation (2.7) is a stochastic process itself. This Section investigates its statistical properties as it regards its spectral structure with a focus on its second eigenvalue or equivalently on the associated dominant decay exponent Λmp. For the sake of simplicity, we consider in this Section two-dimensional systems such as the TPSF. From the previous Section, the construction of the mapping matrix depends essentially on the two integers: N2

c, expressing the total number of boxes in the discretization of the mixing domain, and N2

p corresponding to the total number of particles per box. As Npincreases, a statistical averaging over a larger ensemble of fluid particle

(32)

trajec-3.4. Constructing mapping matrix 19

tories is performed, the effects of stochastic fluctuations on the spectral properties of Φbecome less pronounced as Np increases, see previous Section 3.4.1. We performed the statistical analysis by considering an ensemble of Nr = 103 realizations of the diffusive mapping matrix. Let hΛmpi and σΛ2 be the mean and the variance of the dominant decay exponent Λmp, respectively, obtained with respect to the ensemble of Nrrealizations of Φ: hΛmpi = Nr−1

PNr

α=1Λmp,α, σ2Λ= Nr−1Pα=1Nr (Λmp,α− hΛmpi)2, where Λmp,αis the dominant decay exponent of the α-th realization. Results are lim-ited to Nc = 150 as for higher grid resolutions the computational cost in evaluating the spectral structure of Nr matrices becomes prohibitive.

Figure 3.4 depicts the behavior of hΛmpi vs Np for three values of the flow period at Nc = 50, P e = 103 (panel (a)) and Nc = 150, P e = 104 (panel (b)). The horizontal lines correspond to the values of Λ(P e) obtained by means of Fourier analysis of the advection-diffusion equation. The values of the P´eclet number (P e = 103 _and P e = 104_{) are chosen so that the effect of numerical diffusion for this grid resolution} (Nc= 50 and Nc= 150) is negligible, see Section 3.6.4.

0 20 40 60 80 100 0.1 0.2 0.3 0.4 0.5 < Λm p > Np 0 20 40 60 80 100 0.1 0.2 0.3 0.4 < Λm p > Np (a) (b)

Figure 3.4: Dominant decay exponent hΛmpi vs the number of particles per unit interval Np

for the TPSF at Nc = 50 and P e = 103 (panel (a)); Nc = 150 and P e = 104

(panel (b)). Symbols ‘△’ refers to Tp= 0.8; ‘’ to Tp= 1.18; ‘♦’ to Tp = 1.6. The horizontal lines represent the corresponding values of Λ obtained by means of Fourier analysis.

Starting from Np ≥ 10 the average value of Λmp practically coincides with that of Λ(P e).

More interesting is the scaling of σΛ (depicted in Figure 3.5) for Tp= 0.8, 1.18, 1.6 at Nc= 50 for P e = 103and Nc= 150 for P e = 104, as can be observed

σΛ∼ Np−1. (3.7)

In order to provide a direct quantification of the statistical error induced by stochas-ticity in the construction of Φ, it is convenient to consider the normalized squared root of the variance σ∗

Λ, namely, σ∗Λ=

σΛ hΛmpi

(33)

101 102 10−5 10−4 10−3 10−2 A B σΛ Np 10 1 102 10−6 10−4 10−2 A B σΛ Np (a) (b)

Figure 3.5: Variance σΛ vs Npfor the TPSF at Tp= 0.8 (symbols ‘△’); Tp= 1.18 (symbols ‘’) Tp= 1.6 (symbols ‘♦’) at Nc= 50 for P e = 103 (panel (a)) and Nc= 150 for P e = 104 _{(panel (b)). Lines (A) and (B) represent the scaling σ}

Λ∼ Np−1.

which is depicted in Figure 3.6 for several values of the flow period and of the P´eclet number. 100 101 102 10−5 10−4 10−3 10−2 A B

σ

∗ Λ

N

p

Figure 3.6: Normalized variance σ∗

Λ vs Np for several TPSF protocols: Tp = 0.8, Nc= 50 and P e = 102 _{(symbols ‘◦’); T}

p= 1.18, Nc= 80 and P e = 102 (symbols ‘+’);

Tp = 1.6, Nc = 80 and P e = 103 (symbols ‘•’); Tp = 1.18, Nc = 150 and

P e = 104_{(symbols ‘△’); T}

p= 1.6, Nc= 150 and P e = 104 (symbols ‘’). Lines (A) and (B) represent the scaling σΛ∗∼ N

−1 p .

From Equation (3.7) it follows that σ∗

Λ ∼ Np−1 and this power-law scaling sets in starting from low values of Np ≥ 5. Regarding the values attained by σ∗Λ, starting from Np= 10, σΛ∗ < 10−3 and for Np= 50 its numerical value smaller than 10−4 for all considered flow protocols. This means that, for Np> 10 the error in the estimate of Λmp is order of 0.1-1 percent of the value of the dominant decay exponent. This

(34)

3.4. Constructing mapping matrix 21

indicates that, at least for the TPSF, from a single realization of the diffusive mapping matrix, the statistical error in the estimate of Λmpcan be kept well below 1 % provided that Np is taken large enough Np≥ 10. Throughout this Chapter all the simulations of the TPSF refers to Np = 100, and this is the reason why in numerical analysis of Λmp the value of the dominant decay exponent has been estimated from a single realization of the mapping matrix without reporting the error bars.

Figure 3.7 panel (a) shows the normalized variance σ∗Λ for TPSF at Tp = 0.8 and P e = 103_{for several values of N}

c (Nc= 50, 80, 100, 150). Let s∗Λ be the prefactors of the scalings of σ∗

Λas a function of Np(Np≥ 5), i.e., σ∗Λ= s∗Λ/Np. The information on the statistical properties of σ∗

Λ can be visualized in a compact way by depicting the behavior of the prefactor s∗

Λ as a function of the linear grid resolution Nc (see panel (b)). As can be observed, s∗ Λ scales as a function of Nc as s∗Λ∼ Nc−1. 100 101 102 10−4 10−3 10−2 A σ ∗ Λ Np 10 2 10−3 10−2 10−1 s ∗ Λ Nc (a) (b)

Figure 3.7: Panel (a): Normalized variance σ∗

Λ for Tp= 0.8 and P e = 103 for several values of Nc: Nc = 50 (symbol ‘◦’); Nc = 80 (symbol ‘+’); Nc = 100 (symbol ‘△’);

Nc = 150 (symbol ‘’). Lines (A) represents the scaling σ∗Λ ∼ Np−1 Panel

(b): Prefactor s∗

Λ of the normalized variance vs Nc for the TPSF at Tp= 0.8,

P e = 103. The solid line represents the scaling s∗

Λ∼ Nc−1.

In point of fact, the two scalings σ∗

Λ∼ Nc−1and s∗Λ∼ Np−1can be unified into a single expression by observing that the total number of fluid particles is Ntot= (Nc× Np)2, and therefore express the property that σ∗

Λis inversely proportional to the square root of the total particle number, i.e.,

σ∗Λ∼ 1 √ Ntot = 1 Nc× Np . (3.9)

This result can be predicted from the theory of fluctuations as a function of the particle number of elementary statistical mechanics[91, 92]. The result expressed by Equation (3.9), that combines the effect of discretization (Nc) and that of the statistical accuracy (Np) is depicted in Figure 3.8 for several TPSF protocols and several values of Nc at a fixed value of the P´eclet number. Independently of the flow period and of Nc, all the data collapse for Nc× Np≥ 102 into a single master curve corresponding to the scaling (3.9).

(35)

102 103 104 10−5 10−4 10−3 10−2 σ ∗ Λ Np× Nc

Figure 3.8: Invariant rescaling of σ∗

Λvs Nc× Npfor several TPSF protocols at a fixed values of P e = 103_{: ‘◦’ T}

p = 0.8, Nc = 50; ‘♦’ Tp = 1.18, Nc = 50; ‘’ Tp = 1.18,

Nc = 80; ‘•’ Tp = 1.6, Nc = 100; ‘△’ Tp = 0.8, Nc = 150; ‘∗’ Tp = 1.18,

Nc= 150. The solid line represents the scaling σΛ∗ = a/(Nc× Np) with a = 5.

3.5 Observations on spectral properties

This Section addresses the higher-order spectral structure of the diffusive mapping matrix and the generic behavior of Λmp as a function of the P´eclet number.

A comprehensive qualitative picture of the spectrum of the diffusive mapping matrix is reported in Figure 3.9. Panels (a)-(c) of this figure show the eigenvalue spectrum of the diffusive mapping matrix for different values of the P´eclet number, for the TPSF at Tp= 0.8, obtained by considering Np= 50 and Nc= 100. As can be observed, the eigenvalue spectrum converges (for fixed Nc) towards an invariant structure as P e → ∞ which corresponds to the eigenvalue spectrum of the purely advective mapping matrix (see panel (d)).

Once the numerical properties of the diffusive mapping matrix have been character-ized, consider the behavior of the dominant decay exponent as a function of the P´eclet number. Figure 3.10 depicts the typical behavior of Λmp vs the P´eclet number for the diffusive mapping matrix.

Specifically, the case of Tp = 0.8 and Nc = 50 is considered, and compared with the results for Λ obtained by considering a Fourier-series expansion of the advection-diffusion operator. Up to a given critical P´eclet number, that depends on the dis-cretization, i.e., on Nc (this critical P e value is order of 103 for the case depicted in Figure 3.10), the values of Λmp practically coincides with those obtained by cal-culating the dominant decay exponent of the Floquet operator associated with the advection-diffusion equation (2.2) (line ‘a’). As the P´eclet number increases, Λmp saturates towards a constant value Λc

mp as a consequence of numerical diffusion in-duced by discretization. The value Λc

mpcorresponds to the dominant decay exponent obtained from the analysis of the purely advective mapping matrix.

Two main conclusions can be drawn from the analysis of the data depicted in Figure 3.10. (i) Up to a given critical P´eclet number, the eigenvalues of the diffusive

(36)

map-3.5. Observations on spectral properties 23 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Im( λ ) Re(λ) (a) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Im( λ ) Re(λ) (b) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Im( λ ) Re(λ) (c) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Im( λ ) Re(λ) (d)

Figure 3.9: Panels (a)-(c): Eigenvalue spectrum of the diffusive mapping matrix (Nc= 100) of the sine flow at Tp= 0.8 for several values of the P´eclet number: (a) P e = 104, (b) P e = 105, (c) P e = 106. Panel (d) depicts the eigenvalue spectrum of the purely advective mapping matrix.

101 100 10-1 10-2 10-3 1010 108 106 104 102 Λmp , Λ Pe a b

Figure 3.10: Dominant decay exponent of the mapping matrix Λmp(dots •) vs P e at Tp=

0.8, Nc= 50. Line ‘a’ represents the scaling of the dominant decay exponent Λ of the advection-diffusion equation, dotted horizontal line ‘b’ is the limit value Λc

mp of Λmpin the purely advective case.

ping matrix (obtained by applying the simple stochastic approach) provide a reliable and accurate approximation for the eigenvalues of the propagator of the

(37)

advection-diffusion equation; (ii) coarse-graining introduces an additional diffusive contribution (numerical diffusion) that superimposes to the diffusive term associated with random fluctuations. Therefore, from the analysis of Λc

mp it would be possible to infer quan-titative information on the effect of numerical diffusion influencing the properties of the mapping matrix in the purely advective setting (diffusionless limit).

3.6 Numerical diffusion

3.6.1 Introduction to numerical diffusion

The latter observation indicates that the mapping matrix for a given value of Nc (i.e. for a given coarse-grained resolution) can be regarded as a coarse-grained approxima-tion of an advecapproxima-tion-diffusion process,

∂c

∂t = −∇ · (V c) + Dtot∇

2_{c ,} _(3.10)

associated with a cumulative diffusion coefficient Dtot expressing the superposition of the molecular diffusivity D and of the numerical diffusivity Dnum induced by coarse graining,

Dtot= D + Dnum, (3.11)

where Dnum= Dnum(Nc) depends on the grid resolution. It is important to observe that Equation (3.10) is just a leading-order approximation of the effect of coarse-graining as numerical diffusion may possess highly non-trivial correlation properties that contribute in Equation (3.10) as higher order spatial derivatives of the concen-tration field. Therefore, Equation (3.10) can be regarded as a first starting point to tackle and quantify the action of numerical diffusion on the spectral properties of the matrix Φ. In dimensionless term, the resulting P´eclet number P etotis therefore given by 1 P etot = 1 P e+ 1 P eeff , (3.12)

where P eeff = V0L0/Dnum is the effective P´eclet number associated with numerical diffusion, which equals the value of the P´eclet number pertaining to the purely ad-vective mapping matrix for a given grid resolution Nc. These two issues are carefully addressed in the remainder of this Chapter by considering a detailed analysis of the TPSF for different, and characteristic values, of the flow period Tp.

Equations (3.11)-(3.12) represent a leading-order description of the impact of coarse-graining (numerical diffusion) in the construction of the mapping matrices on a discrete grid. Although several studies address the effect of numerical diffusion on advection-diffusion algorithms[93, 94], mostly with reference to finite-difference schemes, very little is known on the influence of numerical diffusion on spectral proper-ties of advection-diffusion operators, especially for generic mixing (partially or globally chaotic) flows. Consistently with the main motivation of the present analysis, which

Eigenmode analysis of advective-diffusive mixing

Eigenmode analysis of advective-diffusive mixing

Eigenmode analysis of advective-diffusive

mixing

Eigenmode analysis of advective-diffusive

mixing

PROEFSCHRIFT

Oleksandr Gorodetskyi

Contents

Summary

1

Introduction

Microfluidics

Molecular diffusion

Spectral analysis

Numerical diffusion

Compact mapping method

Thesis overview

2

Inclusion of molecular

diffu-sion in the mapping method

Abstract

2.1 Introduction

2.2 Statement of the problem

2.3 Bounded and unbounded mixing systems

2.4 Original mapping method

2.5 Diffusive mapping matrix

2.6 Conclusions

3

Spectral analysis of mixing in

chaotic flows

Abstract

3.1 Introduction

3.2 Problem definition

3.3 Model system: time-periodic sine flow (TPSF)

(a)

(b)

(c)

(d)

3.4 Constructing mapping matrix

σ

N

3.5 Observations on spectral properties

3.6 Numerical diffusion

_Introduction

_{Inclusion of molecular}

_{Spectral analysis of mixing in}