Nelson Poumaëre

Benoît Pier

Florence Raynal

Residence time distributions in unstable channel flow

I. INTRODUCTION

Stirring in laminar flows is an essential step in all mixing processes [START_REF] Villermaux | Mixing versus stirring[END_REF]. Whilst mixing is generally characterized through the decay of concentration variance, measuring stirring efficiency requires additional tools in the absence of molecular diffusion. This can be done for instance using the stretching rate of fluid elements; in the case of flows subject to chaotic advection [START_REF] Ottino | Mixing, chaotic advection, and turbulence[END_REF][START_REF] Aref | Frontiers of chaotic advection[END_REF], Lyapunov exponents and Poincaré sections are more commonly employed. In this article we rather focus on residence time distributions (RTDs), as first introduced by Danckwerts [START_REF] Danckwerts | Continuous flow systems. distribution of residence times[END_REF] in his seminal paper to characterize transport properties of open flows, i.e. flows with clearly defined input and output. The applications are numerous and diverse, for example oil transport in pipelines (Hull & Kent [5]), contaminant dispersion in ground watersheds (Haitjema [6]), or microfluidic devices (Bošković & Loebbecke [START_REF] Bošković | Modelling of the residence time distribution in micromixers[END_REF], Trachsel et al. [START_REF] Trachsel | Measurement of residence time distribution in microfluidic systems[END_REF]). Note that RTD can also be used to obtain qualitative information about mixing by distinguishing between residence times that are greater or smaller than the typical time needed for diffusion to be effective, as explained by Gouillart et al. [START_REF] Gouillart | Measures of mixing quality in open flows with chaotic advection[END_REF].

As schematized in figure 1, for a given fluid system, the RTD provides information about the probability for a fluid particle to cross the system in a certain lapse of time τ = t outt in . As shown by Spalding [START_REF] Spalding | A note on mean residence-times in steady flows of arbitrary complexity[END_REF] and Danckwerts [START_REF] Danckwerts | Local residence-times in continuous-flow systems[END_REF] in steady conditions, the mean residence time τ m is given by

τ m = τmax τmin τ f (τ ) dτ = V Q , (1)
where f (τ ) is the RTD, V the volume in the system between inlet and outlet, and Q the volumetric flow rate. The minimum and maximum residence times τ min and τ max respectively correspond to the fastest and slowest trajectories attainable in the system. Ideally, a mixer should mimic as closely as possible the plug-flow model, where all particles cross the system in the same time τ = τ min = τ max , and for which the RTD is essentially a Dirac pulse. In practice however, the maximum residence time τ max is generally infinite, which implies that the RTD f (τ ) in equation ( 1) must decrease to zero at large times. Whilst the decay may be exponential in a pure extensional flow (Raynal & Carrière [12], Hervieu et al. [START_REF] Hervieu | Residence times of polydisperse dilute suspensions in sheared and extensional flows[END_REF]), the situation is very different when considering flows with walls, since the decay is thus algebraic, with a τ -3 tail [START_REF] Raynal | The distribution of "time of flight" in three dimensional stationary chaotic advection[END_REF], also called heavy tail. This tail invalidates the use of a second-order moment of the distribution f , and thus a standard deviation does not exist. As a measure for the width of the distribution, we will therefore use the mean absolute deviation (Poumaëre et al. [START_REF] Poumaëre | Residence time distributions for in-line chaotic mixers[END_REF]). When a flow is subject to hydrodynamic instabilities, the Eulerian flow structure can be substantially modified. Whether this modification can profoundly impact the mixing, transport and more especially residence time properties is an important question. In the case where instabilities lead to turbulent flow, these properties are indeed deeply impacted as was already emphasized by Hull & Kent [START_REF] Hull | Radioactive tracers to mark interfaces and measure intermixing in pipelines[END_REF], who showed that RTDs in pipelines are more stretched in the laminar case than in the turbulent configuration. They thus recommended that the flow is maintained in a turbulent regime so as to limit mixing between two consecutive pockets of petroleum.

However hydrodynamic instabilities do not always result in fully developed turbulence. In closed flows for instance, the instabilities may lead to a time-periodic flow that exhibits chaotic advection; see for instance Oteski et al. [START_REF] Oteski | Lagrangian chaos in confined two-dimensional oscillatory convection[END_REF] who studied a two-dimensional confined convection flow, or the works by Meunier [START_REF] Meunier | Geoinspired soft mixers[END_REF], along with Lefranc et al. [START_REF] Lefranc | Gas injection into a tilted rotating cylinder[END_REF], taking advantage of triadic resonances in a tilted tank of water to efficiently mix with limited stress, or also Qu et al. [START_REF] Qu | Chaotic mixing in an acoustically driven cavity flow[END_REF] who showed that acoustic streaming could enable global chaos in a parallelepipedic cavity using an oscillatory instability of the flow arising at moderate power of the source.

The case of mixing by instabilities in open flows is even more complex. Lou et al. [START_REF] Lou | Computational fluid dynamics simulations of unsteady mixing in spacer-filled direct contact membrane distillation channels[END_REF] studied how laminar vortex shedding can affect temperature and concentration fields in the vicinity of direct contact membrane distillation channels. Venditti et al. [START_REF] Venditti | Impact of transversal vortices on the performance of opentubular liquid chromatography[END_REF] showed how the presence of transversal vortices greatly improved the performances of open-tubular liquid chromatography by containing axial dispersion. The case of the Couette-Taylor-Poiseuille flow, or spiral Poiseuille flow (SPF) was studied in terms of residence time distributions. Kataoka et al. [START_REF] Kataoka | Ideal plug-flow properties of taylor vortex flow[END_REF][START_REF] Kataoka | Heat/mass transfer in taylor vortex flow with constant axial flow rates[END_REF] noticed that when the vortex flow regime was maintained in the presence of an axial flow rate, the SPF displayed plug-flow properties, with all particles inside one vortex having the same residence time. When axial flow was increased however, the "doubly periodic" regime -wavy vortices -led to profound broadening of the RTD, due to enhanced intermixing between neighbouring vortices.

In the present work, we characterize another type of perturbed open flow in terms of residence time distributions: the plane Poiseuille flow (PPF), in which a fully developed, nonlinear Tollmien-Schlichting wave propagates at constant speed.

Even though hydrodynamic stability is not the main subject of the present article, a quick review of the main results about the linear and nonlinear stability theory for the plane Poiseuille flow is required here.

Among the numerous authors who have concerned themselves with the linear stability of the PPF (for example Lin [START_REF] Lin | On the stability of two-dimensional parallel flows, Part III.-Stability in a viscous fluid[END_REF], Thomas [START_REF] Thomas | The stability of plane Poiseuille flow[END_REF] and Grosch & Salwen [START_REF] Grosch | The stability of steady and time-dependant plane Poiseuille flow[END_REF]), it is Orszag [START_REF] Orszag | Accurate solution of the Orr-Sommerfeld stability equation[END_REF] who eventually established the critical Reynolds number to its now classically known value of Re c = 7696 (with the definition adopted in the present article). Making use of a Chebyshev expansion and a spectral method, he numerically solved the Orr-Sommerfeld equation for the normal velocity component ṽ of an infinitesimal perturbation written in complex form as

u(x, y, t) = ũ(y) exp [i(αx -ωt)] + c.c., (2a)
p(x, y, t) = p(y) exp [i(αx -ωt)] + c.c., (2b)
where α denotes the real wavenumber of the perturbation, ω = ω r + iω i its complex frequency, and c.c. the complex conjugate. For a given setting of Re and α, the linear eigenvalue problem arising from the Orr-Sommerfeld equation is solved, yielding a discrete spectrum of eigenvalues ω and the associated eigenfunctions ũ(y) and p(y). If the growth rate ω i of a given mode is positive, it is exponentially amplified, otherwise it is neutrally stable, or decays asymptotically. Above the threshold value of Re c , to which corresponds the critical wavenumber α c , a continuous band of unstable wavenumbers exists, and for each of these unstable wavenumbers, the eigenmode with the largest positive growth rate is referred to as the Tollmien-Schlichting (TS) wave. Note that the TS mode exhibits sinuous symmetry. In the (Re, α)-plane, the locus of neutrally stable configurations is called the neutral curve, and is reproduced (numerical methods in the next section) in figure 3(a), along with contours of constant growth rates. Interest in the nonlinear stability properties of the plane Poiseuille flow (and shear flows in general) has developed in the wake of advancements in the linear theory, mainly to explain why turbulence is observed in experiments way below (Re ≃ 1 500) the critical Reynolds number (Davies & White [START_REF] Davies | An experimental study of the flow of water in pipes of rectangular section[END_REF]). Nonlinear effects manifest themselves in essentially three ways (Stuart [28]): distortion of the mean flow, modification of the fundamental from linearized theory, and generation of higher harmonics. Essentially all works relied on some kind of spatial Fourier expansion of the total flow

u(x, y, t) = n u (n) (y, t) exp [inαx]. (3)
Early works (Meksyn & Stuart [29], Stuart [START_REF] Stuart | On the non-linear mechanics of hydrodynamic stability[END_REF]) made use of the mean-field approximation, where this expansion is truncated up to |n| < 2, and only the interaction between the mean flow and the fundamental is considered. By truncating the expansion up to the second harmonic (|n| < 3), and assuming the first harmonic to be the normalized eigenmode of linear theory multiplied by an amplitude function of time, Stuart [START_REF] Stuart | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows part 1. the basic behaviour in plane Poiseuille flow[END_REF] and Watson [START_REF] Watson | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows part 2. the development of a solution for plane poiseuille flow and for plane couette flow[END_REF] were able to study the subcritical and supercritical behavior of the flow in the vicinity of critical conditions, notably establishing the existence of equilibrium solutions consisting of steadily traveling two-dimensional finite-amplitude waves. This kind of expansion was used in consecutive works like those of Zahn et al. [START_REF] Zahn | Nonlinear cellular motions in Poiseuille channel flow[END_REF] or Jiménez [START_REF] Jiménez | Bifurcations and bursting in two-dimensional Poiseuille flow[END_REF] to study different solution branches bifurcating from basic plane Poiseuille flow. Ehrenstein & Koch [START_REF] Ehrenstein | Three-dimensional wavelike equilibrium states in plane Poiseuille flow[END_REF], using Fourier expansion and pseudo-arclength continuation, thoroughly computed twoand three-dimensional nonlinear traveling-wave solutions, under both sub-and supercritical conditions.

The present article is concerned with this type of solutions, but considers only those obtained under supercritical conditions, i.e. for linearly unstable base flows.

There are different ways to calculate numerically a RTD in an open flow. In the case of in-line chaotic mixers made of a succession of identical elements of length L each, chaos is global so that a given particle visits the whole periodic domain; a natural method is thus to follow a unique particle over a long lapse of time, long enough for it to travel a distance M × L with M ≫ 1, and to monitor the time needed to cross each portion. The RTD is thus computed as the histogram of the resulting M residence time values [START_REF] Raynal | The distribution of "time of flight" in three dimensional stationary chaotic advection[END_REF][START_REF] Artuso | Dynamical and transport properties in a family of intermittent areapreserving maps[END_REF]. Indeed, the particle naturally visits regions of high velocity more often than regions of low velocity, and the residence time values are already properly weighted; this also implies that the mean value of the residence times in the sample is the Eulerian mean time τ m defined in equation [START_REF] Villermaux | Mixing versus stirring[END_REF]. This method can be qualified as a Lagrangian process, and is available for stationary flows; it is also very cheap, as the different residence time values are obtained simultaneously with the points when calculating the Poincaré section. When the flow is not chaotic and stationary, a simple technique for the calculation of RTD is to seed homogeneously a large number of particles at the entrance of the section, follow them on the longitudinal distance L and note the residence times; the RTD is then the pdf of the residence times weighted by the longitudinal velocity at entrance. This method of weighted residence times is best known in porous media studies, in the context of the Continuous Time Random Walk (CTRW) method, see for example the early work of Saffman [START_REF] Saffman | A theory of dispersion in a porous medium[END_REF], and more recently Dentz et al. [START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF] and Puyguiraud et al. [START_REF] Puyguiraud | Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media[END_REF]. It can however become difficult to implement when the flow-field is time-dependent, since the velocity at entrance is both position and time dependent. The flow considered here is not only time-dependent, but also non chaotic as it is two-dimensional and can be made stationary by a Galilean change of frame (see Cox et al. [START_REF] Cox | Chaotic advection of irrotational flows and of waves in fluids[END_REF]). We propose here to generalize the Lagrangian method by allowing the particle to diffuse and thus "sample" the entire cross section of the flow. This was first proposed by Spalding for stationary flows [START_REF] Spalding | A note on mean residence-times in steady flows of arbitrary complexity[END_REF], but never used in practice. The idea is here to have a diffusion large enough so as to favour the sampling, but at the same time small enough so as not to alter much the RTDs obtained.

The paper is organized as follows. In section II, the geometry studied and general equations are recalled, followed by a quick review of the linear and nonlinear stability of the plane Poiseuille flow. Then we explain the method for computing the nonlinear Tollmien-Schlichting waves, and introduce a reference case typical of the situation encountered in the problem. Section III explains the diffusive method for the calculation of the RTDs. The procedure is then validated on the parabolic plane Poiseuille flow. Section IV first focuses on the reference case, and explains the links between the properties of the saturated flow and the corresponding RTD. From what was learned with the Lagrangian method, we then propose an Eulerian method, without diffusion, for the calculation of the RTD. Finally we vary the physical parameters of the problem, and give the conclusions and perspectives in the last section.

II. PROBLEM FORMULATION AND METHODS

A. Geometry and governing equations

Throughout this investigation, the flow between two infinite parallel plates is considered, governed by the incompressible Navier-Stokes equations,

∂ t u + (u • ∇)u = ν∆u - 1 ρ ∇p, (4a)
∇ • u = 0, (4b)
where ρ and ν respectively denote the constant density and kinematic viscosity of the fluid. The entire study takes place within a two-dimensional framework, assuming neither spanwise variations nor spanwise velocity component. Adopting Cartesian coordinates, with x denoting the streamwise and y the wall-normal coordinates, the velocity u and pressure p then depend on position x = xe x + ye y and time t. The channel is of width h and the fluid put into motion by a constant streamwise pressure gradient G. No-slip boundary conditions prevail along the walls at y = ±h/2 (see figure 2). The steady base Poiseuille flow solution follows the classical parabolic profile

U (y) e x = 3 2 U m 1 - 4y 2 h 2 e x , (5)
with U m = h 2 12ν G the mean velocity. Throughout this article, the definition of the Reynolds number is based on the channel width h and the mean velocity U m :

Re = hU m ν . (6)
Note that the flow dynamics is entirely controlled by the choice of the Reynolds number. For deriving numerical values, without loss of generality, we specify two dimensional parameters h = 1 and U m = 1.

B. Linear and nonlinear stability of the plane Poiseuille flow

In order to investigate the stability of the Poiseuille profile, the total instantaneous flow fields are separated into basic and perturbation quantites as

u tot (x, y, t) = U (y) e x + u pert (x, y, t), (7)
p tot (x, y, t) = -Gx + p pert (x, y, t), (8)
whether the perturbation quantities are of small amplitude (for linear stability analyses) or not (for computing the saturated fully developed nonlinear flow solutions). Replacing these expansions in the Navier-Stokes equations, the equations governing the evolution of the perturbation fields u pert and p pert are obtained as

∂ t u pert = -L(u pert ) -N (u pert ) + ν∆u pert - 1 ρ ∇p pert , (9a)
∇ • u pert = 0, (9b)
where L(u pert ) = U ′ v pert e x + U ∂ x u pert is the linear term accounting for the advection by the base flow and N (u pert ) = (u pert • ∇)u pert is the nonlinear self-advection term. For the rest of this paper, we will only consider the perturbation velocity and pressure fields u pert and p pert and, for simplicity, omit their subscript. When considering small-amplitude perturbations, the quadratic term N (u) is neglected and the flow fields are sought in normal-mode form as

u(x, y, t) = ũ(y) exp [i(αx -ωt)], (10a)
p(x, y, t) = p(y) exp [i(αx -ωt)], (10b)
where α is the streamwise wavenumber and ω the frequency. Substitution of this expansion into the linearized perturbation equations results in a classical eigenvalue problem, the solution of which yields the frequency spectrum in the complex ω-plane for a given wavenumber α (Schmid & Henningson [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]). A positive imaginary part ω i of the frequency associated by the linear dispersion relation to a given wavenumber α signals a linear instability. An initial spatially periodic perturbation of wavenumber α is then exponentially amplified in time with growth rate ω i so as to reach finite amplitude levels. The dynamics of the fully developed perturbation is then governed by the complete equations (9a,9b) since the nonlinear terms can no longer be neglected. After expanding the spatially periodic finite-amplitude perturbation as

u(x, y, t) = n u (n) (y, t) exp [inαx], (11a)
p(x, y, t) = n p (n) (y, t) exp [inαx], (11b)
and substituting in the nonlinear evolution equations (9a)(9b) leads to a system of coupled partial differential equations governing the temporal evolution of the spatial Fourier components u (n) (y, t) and p (n) (y, t). This is equivalent to considering the dynamics in a finite domain of length 2π/α with periodic boundary conditions in the streamwise x-direction.

C. Numerical computation of the Tollmien-Schlichting waves

The numerical implementation closely follows the methods developed by Pier & Schmid [START_REF] Pier | Linear and nonlinear dynamics of pulsatile channel flow[END_REF] in the context of pulsatile channel flow. The spatial discretization involves a Chebyshev collocation method, using N c Gauss-Lobatto collocation points to cover the channel diameter -h/2 ≤ y ≤ h/2. Thus, the differential operator ∂ y is recast as a matrix operating on the array consisting of the values at the collocation points.

After this spatial discretization, the linear eigenvalue problem governing the eigenfunctions ũ(y) and p(y) (10a),(10b) is recast as an algebraic eigenproblem, the solution of which then yields the linear dispersion relation.

For nonlinear temporal evolution problems, the streamwise Fourier expansions (11a)(11b) for u(x, y, t) and p(x, y, t) are truncated at |n| < N h . Then, the Chebyshev collocation method in y transforms the governing partial differential equations (9a),(9b) into an algebraic system of coupled ordinary differential equations in time. Temporal integration is carried out using a predictor-corrector fractional-step method derived from the work of Goda [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows[END_REF], improved by Raspo et al. [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF]. More details of the general method can be found in [START_REF] Pier | Linear and nonlinear dynamics of pulsatile channel flow[END_REF][START_REF] Pier | Dynamique des écoulements ouverts : instabilités et transition, courbure et rotation[END_REF]. The numerical implementation in C++ is based on the 'home-spun' PackstaB library [START_REF] Pier | PackstaB[END_REF]. The required numerical effort is reduced by explicitily taking into account symmetry properties prevailing for sinuous or varicose perturbations: the different flow components are then either symmetric or antisymmetric in y. Thus only half of the collocation points need to be considered, using the associated discretised differential operators of corresponding symmetry.

The essential properties of linear and nonlinear Tollmien-Schlichting waves are illustrated in figure 3. First, maps of the linear growth rate and phase velocity are computed over the (Re, α)-plane, by solving the algebraic eigenvalue problem and monitoring the most unstable eigenvalue of the spectrum obtained for each choice of Re and α. Figure 3(a) shows isolines of constant growth rate ω i along with the neutral curve. Note that the most unstable eigenmode is always of sinuous symmetry. The associated phase velocity ω r /α is plotted as dashed lines in Figure 3(b). The region inside the neutral curve ω i = 0 in Figure 3(a), corresponds to linearly unstable configurations, for which temporal integration of the nonlinear evolution problem leads to fully developed perturbations that saturate at finite amplitude, after an initial phase of exponential growth. For the parameter ranges considered in the present investigation, the perturbation energy is found to converge towards a constant level, and the saturated finite-amplitude perturbations propagate as a whole at a constant velocity without modulations. This behavior suggests the absence of secondary instablities, that would lead to more complex dynamics. Thus, in this finite-amplitude regime, the temporal evolution of the different Fourier components in (11a)(11b), is of the form u (n) (y, t) = ũ(n) (y) exp(-inω nl t), for some real-valued nonlinear frequency ω nl . By monitoring the different components during the temporal integration and computing the frequency ω nl , it is thus verified that, after convergence, the perturbations may indeed be written as a steadily traveling nonlinear wave

u(x, y, t) = |n|<N h ũ(n) (y) exp [in(αx -ω nl t)]. (12)
This nonlinear Tollmien-Schlichting wave travels with a phase velocity c nl = ω nl /α; in a reference frame comoving with the wave at speed c nl , the perturbation fields are thus steady in time and spatially periodic with period λ = 2π/α, the wavelength of the perturbation. The solid lines in Figure 3(b) show contours of the nonlinear phase velocity c nl , computed for the parameter region corresponding to linear instability. Note that from equation ( 12), the velocity field is only a function of wall-normal position y and phase ϕ with

ϕ = αx -ω nl t mod 2π with 0 ≤ ϕ ≤ 2π. (13)
Similarly, we define the periodic longitudinal position x ′ in the moving reference frame as

x ′ = ϕ/α = x -c nl t mod λ with 0 ≤ x ′ ≤ λ. (14)
D. A reference case: α = 2 and Re = 10 000

As an example, the result of a simulation for Re = 10 000 and α = 2 can be seen on figure 4, which displays the corrected mean profile and a contour map of the streamfunction of the total flow in the frame comoving with the wave, frame in which the flow is stationary. This configuration is linearly unstable with ω = 0.74967 + 0.00670i. For this case we have λ = λ ref = 2π/α = π, the nonlinear frequency is ω nl = 1.17, yielding a nonlinear wavespeed c nl = 0.587 and a time period T = 2π/ω nl = 5.36. The longitudinal extent for the contour map is the spatial pseudo-period λ = π. The distorted mean flow is compared to the base flow on the left, an important feature of this distortion being that the maximal value of the mean flow is greater than in the case of the base flow. On the right the streamlines of the total flow in the frame comoving with the wave are visualized. The sinuous structure of the wave is clearly recognized here in the form of the Kelvin-Stuart cat's eyes pattern, a series of cellular recirculation zones traveling downstream of the flow at a velocity equal to the wave speed. These waves also exhibit small near-wall regions where the longitudinal velocity, measured in the laboratory frame, is negative. Basically, if a particle is close to the wall and happens to be just below (or above) a cat's eye, it is traveling upstream because of the vortex motion.

-h/2 0 h/2 0 U m 2U m (
In figure 5, we represent snapshots of the trajectories of 2 500 particles over one period of the flow. As was described by Weiss [START_REF] Weiss | Transport and mixing in traveling waves[END_REF] for the general case of a traveling wave, the flow is essentially separated into three distinct regions: the core region where particles travel faster than the wave, the near wall region where they lag behind it, and the trapping regions where particles are bound to travel around the center of the cat's eye cells. In particular, particles traveling slowly near the separatrix are easily distinguishable at time t = T . 5. Snapshots of the positions of 2 500 particles, initially uniformly placed in the channel over two wavelengths of the TS wave, for different instants between t = 0 and t = T . The parameters are those of the reference case, Re = 10 000 and α = 2. In order to visualize the velocity of the wave, the blue and red vertical lines drawn at each side of the initial blob are artificially moved at the wave speed.

t = 0 t=T/5 t=2T/5 t=3T/5 t=4T/5 t=T 0 ⩽ x ⩽ 6λ FIG.

III. COMPUTATION OF RESIDENCE TIME DISTRIBUTIONS

In this section we develop a general method for the calculation of RTDs, that could in principal be available for any type of in-line flow or mixer. The method is validated afterwards in the case of the non perturbed parabolic Poiseuille flow.

A. The Lagrangian method

As explained in the introduction, we propose here to compute the residence time distributions on a given longitudinal distance L using a single long trajectory, and allow the particle to diffuse so as to sample the entire cross-section of the flow. The advection of a single diffusive particle can be modeled with a generalized Langevin equation

dx dt = u(x, t) + ζ(t), (15)
where ζ(t), a Gaussian decorrelated process such that ⟨ζ i (t)ζ j (t ′ )⟩ = 2Dδ ij , with D the diffusion coefficient and δ ij the Krönecker symbol, models the molecular diffusion [START_REF] Aref | Enhanced separation of diffusing particles by chaotic advection[END_REF]. Although diffusion is artificial in the method (it just allows to visit the whole width), we introduce the Péclet number based on the mean velocity

U m , Pe = hU m D . (16)
We then proceed as follows: we solve equation ( 15) for a single diffusive particle. When the particle enters the jth element at plane x = (j -1) L and time t j-1 , we note its vertical position y j-1 and the corresponding phase ϕ j-1 = α × (j -1) Lω nl t j-1 ; then we calculate the residence time τ j = t jt j-1 in this element as the time it takes to reach the exit plane x = j L once entered. Note that, if the particle enters the first element at time t 0 = 0, the total time at the outlet of the jth element is thus t = j i=1 τ i . In practice, we follow the particle on a very long time so as to obtain a large sample of M data (τ j , y j , ϕ j ) 1≤j≤M . The process is illustrated in figure 6. The RTD can thus be obtained as the histogram of the resulting M residence time values τ j ; the FIG. 6. Schematic of the partitioning of the flow in sections of length L, along with a particle trajectory and retained data (τj, yj, ϕj).

other parameters (y j , ϕ j ) will be used for other purposes, see section IV. In practice, we store these data for L = λ. The case L = n × λ (2 ≤ n ≤ 10) is obtained by summing the residence times on n consecutive elements, resulting in a smaller sample of size M/n and a mean residence time n times larger than for n = 1. As explained in the introduction, because of the heavy tail of the distribution, the dispersion of the sample (τ j ) 1⩽j⩽M of residence times is measured through the mean absolute deviation

σ 1 = 1 M M j=1 τ j τ m -1 , (17)
where τ m = L/U m is the mean residence time on the distance L.

Since the only purpose of diffusion here is to enable the computation of the residence time distribution of the underlying flow, care must be taken when choosing the Péclet number. Indeed, diffusion must allow the particle to sample the whole flow, without affecting the RTD.

First of all diffusion must not statistically modify the time needed to travel the distance L, typically t conv = L/U m . This means that (i) the characteristic longitudinal diffusive time t diff,x = L 2 /D must be large compared to t conv , i.e. t diff,x ≫ t conv . Since the traveling time depends on the streamline, this also implies that the particle must statistically not be allowed to change streamline when traveling the distance L; this involves that (ii) the characteristic transversal time t diff,y = h 2 /D must be large compared to t conv , i.e. t diff,y ≫ t conv . Finally, the trajectory must be long enough for the statistics to be converged. Denoting by M the size of the sample, the particle has traveled a total length M × L in the x-direction, associated to a typical lapse of time M × t conv . The statistics are converged if the particle has statistically visited all the streamlines of the flow. As the particle changes streamline thanks to diffusion, this implies that (iii) t diff,y ≪ M t conv . Conditions (i)-(iii) can be rewritten as:

Pe ≫ h/L , ( 18
) Pe ≫ L/h , (19)
M ≫ Pe h/L . (20)
Note that conditions ( 18) and ( 19) are very similar. Because in most cases the length L considered is such that L > h, condition ( 18) is less restrictive, and condition [START_REF] Lou | Computational fluid dynamics simulations of unsteady mixing in spacer-filled direct contact membrane distillation channels[END_REF] along with condition ( 20) is sufficient in practice. Note finally that condition (20) can be rewritten L tot ≫ Pe h, where L tot = M × L is the total longitudinal distance traveled by the particle. Therefore, this last criterion should be fulfilled equally for a couple of parameters (L, M ) or (10 L, M/10), since both cases correspond to the same total longitudinal length traveled L tot .

B. Choice of the travel distance L

In spatially periodic flows, a natural choice for the travel distance L on which the RTD is computed could be the period or a multiple number of periods, see for instance Khakhar et al. [START_REF] Khakhar | A case study of chaotic mixing in deterministic flows: The partitioned-pipe mixer[END_REF] or Mezić et al. [START_REF] Mezić | Residence-time distributions for chaotic flows in pipes[END_REF] in the partitioned-pipe mixer. In the configuration at hand, the flow is only pseudo-periodic and there is no true spatial periodicity. However, as the flow is periodic in the frame comoving with the wave with period λ = 2π/α, a natural possibility is to use L = n × λ, with n reasonably chosen between 1 and 10 [14].

C. Validation of the method

We propose to validate those conditions on the unperturbed plane Poiseuille flow, with parabolic profile of equation [START_REF] Hull | Radioactive tracers to mark interfaces and measure intermixing in pipelines[END_REF]. In such a flow a non diffusing tracer goes straight without changing streamline; it is therefore a drastic candidate to test the method. Moreover, this flow is simple enough so as to obtain an analytical expression for the RTD (valid in the limit Pe → ∞)

f Poi (τ ) = 1 3 τ 2 m τ 3 1 - 2τ m 3τ -1/2 , (21)
with the corresponding analytical value of the mean absolute deviation, see appendix A,

σ 1,Poi = 2/(3 √ 3) ≃ 0.385. (22)
Note that, unlike the perturbed flows considered in the present study, the parabolic flow-field has no natural length; thus equations [START_REF] Kataoka | Ideal plug-flow properties of taylor vortex flow[END_REF][START_REF] Kataoka | Heat/mass transfer in taylor vortex flow with constant axial flow rates[END_REF] are valid whatever the choice of L. Because we need a travel distance for the validation, we choose L as the wavelength of the reference perturbed case depicted in paragraph II D with wavenumber α = 2, and equivalently λ ref = π.

In order to test which Péclet numbers Pe and sample size M may fulfil conditions (18-20), we proceed as follows: for a given Péclet number, we compute a single very long trajectory, long enough to yield a converged value for σ 1 : a total longitudinal distance L tot = M × π with M = 10 8 is fairly sufficient in practice. Then we compute the resulting RTDs (normalized by their corresponding mean residence time τ m ) for the different Péclet numbers Pe and travel distances L, and calculate the associated mean absolute deviations σ 1 .

Figure 7 (a) displays the RTDs in log-log scale, so as to focus on the long residence times behavior. The figure clearly illustrates criterion (ii) associated with condition [START_REF] Lou | Computational fluid dynamics simulations of unsteady mixing in spacer-filled direct contact membrane distillation channels[END_REF]: in the presence of diffusion, a particle located close to the wall (where the velocity is weak) cannot remain forever on its streamline, and arbitrary long residence times are no more possible. This results in a cut-off in the algebraic τ -3 tail obtained without diffusion [START_REF] Poumaëre | Residence time distributions for in-line chaotic mixers[END_REF]. When increasing the Péclet number Pe for the fixed travel distance L = λ, diffusion effects are less important, and the cut-off is shifted toward higher residence times as observed in figure 7 (a). Reversely, when increasing the travel distance L for a given Péclet number, a particle is more likely to change streamline so that diffusion effects are more pronounced, and the cut-off is shifted toward smaller residence times, as illustrated in the case Pe = 10 5 for L = π and L = 10 π. Note also that the tails of the RTDs of the cases (Pe = 10 4 , L = π) and (Pe = 10 5 , L = 10 π) are almost identical. Overall, when the quantity Peh/L increases, the cut-off is shifted toward higher residence times, as expressed in condition [START_REF] Lou | Computational fluid dynamics simulations of unsteady mixing in spacer-filled direct contact membrane distillation channels[END_REF]. Figure 7 (b) displays the RTDs in lin-lin scale, so as to focus on the small residence times behavior. The figure now clearly illustrates criterion (i) associated with the less restrictive condition (equation 18): diffusion can modify the time needed to travel the distance L, leading to residence times that can be significantly smaller than the minimum theoretical time without diffusion, τ min = 2/3 × τ m , see equation [START_REF] Kataoka | Ideal plug-flow properties of taylor vortex flow[END_REF]. When Pe is increased with a fixed travel distance, diffusion effects are less important, and the RTDs become more peaked around τ min , as shown in the figure for L = λ. Identically, when the travel distance L is increased at fixed Péclet number, diffusion is more likely to move the particle equally forward and backward on the larger distance, and the RTD is more peaked, as illustrated in the case Pe = 10 5 for L = λ and L = 10λ. Note finally that the shape of the RTDs at small time is very similar for the two cases (Pe = 10 5 , L = 10 π) and (Pe = 10 6 , L = π). Overall, the RTDs are more peaked on τ min when increasing the quantity PeL/h, as expressed in condition [START_REF] Qu | Chaotic mixing in an acoustically driven cavity flow[END_REF].

In order to obtain a quantitative criterion, we now turn to the computation of the mean absolute deviation. The results are compared to the theoretical mean absolute deviation σ 1,Poi (equation 22); we consider that a given point is acceptable if it lies within 1% of this theoretical value.

Figure 8 (a) shows the mean absolute deviations obtained for the different trajectories corresponding to different Péclet numbers, for varying values of the travel distance, L = n π with 1 ≤ n ≤ 10. For each point the size of the sample is M/n, that is, larger than 10 7 (since M = 10 8 ). The data are plotted as a function of the quantity Peh/L corresponding to the parameter involved in the most restrictive condition (equation 19). For small values of Peh/L, the values obtained are underestimated, which is explained by the fact that the cut-off in figure 7 (a) is too much shifted toward small residence times. As seen in the figure, all points corresponding to Pe = 10 5 or higher are acceptable.

Condition [START_REF] Venditti | Impact of transversal vortices on the performance of opentubular liquid chromatography[END_REF] is related to the size of the sample M : in order to optimize the method, we must evaluate the minimum size M for which the mean absolute deviation is sufficiently converged. In practice, we use the long trajectories calculated previously and consider subsets of increasing size M , (τ i ) 1⩽i⩽M , with M max = 10 8 . The results are shown in figure 8 (b), where the computed values of σ 1 are plotted as a function of the non dimensional parameter M L/hPe involved in equation [START_REF] Venditti | Impact of transversal vortices on the performance of opentubular liquid chromatography[END_REF]. Note that we recover that for small values of Pe h/L, as discussed before, the curves do not converge toward the correct theoretical value. This is the case for instance for the parameters (P e = 10 4 , L = 10 π), whilst the theoretical value is reached for the same Péclet number and a smaller travel distance, i.e. for the couple (P e = 10 4 , L = π). However, although they do not converge toward the same value, the two curves behave similarly with the parameter M L/hPe; indeed, for a given value of the parameter, as explained before, the total distance L tot = M × L is the same (with a sample 10 times smaller for L = 10 π than for L = π h). Therefore the convergence is identical since the two curves are calculated with the same part of the trajectory. This is even more striking for the two sets of parameters (Pe = 10 5 , L = π) and (Pe = 10 5 , L = 10 π), where the two curves converge toward the same parameter and are nearly identical.

Finally, in all cases, the convergence toward the final value is effective for M L/hPe ≥ 300.

To sum up, the Péclet number must be equal to 10 5 or higher, and the size M of the sample must satisfy M ≥ 300 × Pe h/L. Similar thresholds were obtained when considering different trajectories (different initial points, and of course different aleatory diffusion increments). In order to have a reasonable sample size, we therefore choose the minimum acceptable Péclet number, Pe = 10 5 , for the computation of residence times hereafter. The size of the sample is thus chosen as the minimum acceptable size, M ≈ 300 Pe h/L = 3 × 10 7 h/L. For instance, for the reference case for which λ = π (with h = 1), we need a sample of size M ≈ 10 7 .

IV. RESULTS

In this section, we first consider the reference case, the nonlinear TS wave introduced in section II D, for which we study the influence of the travel distance L on the RTD. From this example we show how to calculate RTDs properly with an Eulerian method (without diffusion) in the case of TS waves. Then we vary the wavenumber and Reynolds number.

As explained in the previous section, the Péclet number chosen for the calculation of the RTD is Pe = 10 5 , and the minimum sample size is M = 10 7 ; numerical parameters for the computation of the wave are N cheby = 65 and N h = 5.

In the following, we refer to the time τ wave = L/c nl it takes for one wave train to travel the distance L as the wave residence time; similarly, τ min is the minimum of the residence time for the distance L. We introduce the nondimensional wave residence time τwave = τ wave /τ m = U m /c nl and also τmin = τ min /τ m = U m /U max . Note that both τwave and τmin are independent on the travel distance L.

A. The reference case α = 2 and Re = 10 000

We first begin with the reference case already described in paragraph II D, for which α = 2, Re = 10 000, and λ = 2π/α = π. For this case we compute the RTD, using a single diffusive tracer, following the method outlined in section III; a single trajectory is computed throughout M = 12 × 10 6 sections of length L = λ, resulting in a sample (τ i , y i , ϕ i ) of size M . Figure 9 shows the RTDs obtained for the travel distances L = λ and L = 10 λ. For the purpose of comparison, the RTD of the parabolic Poiseuille flow, obeying to equation [START_REF] Kataoka | Ideal plug-flow properties of taylor vortex flow[END_REF], is plotted with a dotted line.

First we note that the minimum residence time is significantly smaller for this reference case than for the Poiseuille flow. Indeed, the maximum mean velocity at the center of the flow is greater than that of the Poiseuille flow, as previously seen in figure 4(a).

The second characteristic is the existence of a peak located at the nondimensional wave residence time τwave . The peak is much higher and more marked for the larger travel distance L = 10λ.

To explain this behavior, the samples are visualized in the (ϕ, y, τ ) plane in figure 10 for increasing travel distances, with points colorized according to the corresponding residence time. The figure essentially shows that a point starting in a recirculation region is associated with a nondimensional residence time close to τwave (green color). This behavior is more pronounced as L is increased. Indeed, focus on the lower recirculation region in figure 10(a) (L = λ): the point at the center of the recirculation (which is fixed in the moving referential) travels exactly at the velocity of the wave in the laboratory frame, and its residence time is exactly equal to τ wave . Consider now a point initially slightly above the center: the gyratory motion sends the particle to the low-velocity, near-wall region, and the resulting residence time is higher. Conversely, a point initially just below the center is entrained toward the central jet, where the velocities are higher, resulting in a lower residence time. When L is increased (L = nλ, n ≥ 2), points that are in the recirculation region turn around the center in the wave referential. If they make an exact number of turns, their residence time is exactly equal to τ wave = nλ/c nl ; if they do not, their residence time may be slightly higher or lower. Consider finally the extreme case of a point located on one of the separatrices: if located on the upper separatrix, it moves to the right toward the next saddle point on the right, and the distance covered in the moving frame is at most equal to λ. Therefore, the residence time in the laboratory frame is less than that of the central point, such that (n -1) × λ/c nl ≤ τ ≤ n × λ/c nl . On the contrary, a point situated on the lower separatrix moves to the left in the wave referential, with a maximal distance covered λ in the moving frame, and the residence time is higher and verifies n × λ/c nl ≤ τ ≤ (n + 1) × λ/c nl . Overall, all points in the recirculation have a residence time such that (n -1) × λ/c nl ≤ τ ≤ (n + 1) × λ/c nl . Finally, the nondimensional residence time τ of a point in the recirculation is such that

(1 -1/n) τwave ≤ τ ≤ (1 + 1/n) τwave , (23)
which approach τwave as n increases.

We now focus on long residence times, corresponding to particles located near the walls, where the longitudinal velocity is weak. Figure 11 shows the RTD obtained for different travel distances ranging from L = λ to L = 10λ in log-log scale: although we globally recover the t -3 behavior and its cut-off due to diffusion [START_REF] Poumaëre | Residence time distributions for in-line chaotic mixers[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF][START_REF] Puyguiraud | Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media[END_REF], some residence times seem underrepresented, so that a peculiar pattern of peaks and valleys is observed for long residence times. The reason for these less probable residence times is of Eulerian nature. Indeed, the total flow field exhibits some small regions delimited in blue in figure 4, near the upper and lower walls and right at the vertical of the recirculation regions, where the longitudinal velocity is negative in the laboratory frame: in these regions a particle actually travels upstream. This implies that a particle cannot start nor end its trajectory on a section of length L in these particular regions (one cannot cross the start or finish line moving backwards). This is particularly visible in figure 10, which represent a 2D histogram of points entering in an element: the regions delimited in blue in figure 4(b) roughly correspond to those empty of points in figure 10. A blow-up of the bottom near-wall region is shown in figure 12. As can be seen on the figure, the forbidden region is rather larger than the zone delimited in blue. This can be explained as follows: when a particle is in a zone of negative longitudinal velocity, it goes backwards. Therefore it must go forward for some time before crossing the finish line. Therefore the forbidden region extends to the left of the blue zone. More generally, when approaching the blue zone, the velocity decreases and approaches zero; therefore a particle is very unlikely to cross the finish line in the vicinity of the blue zone, and the density of points decreases as approaching this particular region.

In order to understand why those forbidden regions lead to less probable residence times, let us consider a particle near the bottom wall that has just entered the j-th element: because it is located near the wall, its longitudinal velocity is very weak (and sometimes negative when crossing the region delimited in blue in figure 4 b). Therefore it travels a long time before leaving the j-th element of length L. The trajectory of such a particle in the moving frame is depicted in figure 12: in this referential the particle moves very rapidly toward the left, at a velocity close to c nl . When the particle leaves on the left in the periodic frame of figure 12 it reenters on the right, on the same trajectory if diffusion is negligible; the time taken to cross once from left to right is roughly λ/c nl . Because it travels a long time, it moves over a long distance in the referential of the wave, and reenters many times in the domain. In blue, region of negative longitudinal velocities in the laboratory frame. In red, trajectory of a particle near the bottom wall in the referential of the wave. The points correspond to the different positions of the diffusing particle in the (x ′ = ϕ/α, y) plane when entering the jth element at x = (j -1) × λ (figure 6) and can be seen as a Poincaré section of the diffusing trajectory.

Let us now focus on figure 12: a particle starting outside the forbidden region is always allowed to travel a distance equal to n × λ in the moving frame (corresponding to the distance L in the laboratory referential), because then the finish line is also outside the forbidden region. This very probable configuration corresponds to a travel time τ ≈ m × λ/c nl , where m is an integer which is at least such that m ≥ 2 (because of the long times involved in this near-wall region). On the other hand, because of the large forbidden and rarefied region, it is more likely to end in the forbidden region when traveling a distance equal to (m + 1/2) × λ, which makes this distance in the moving frame less probable. This corresponds to a rare travel time τ ≈ (m + 1/2) × λ/c nl .

Let us come back now to figure 11: in the case L = λ, the mean travel time is τ m = λ/c nl . Therefore the most probable nondimensional times are τpeak ≈ m × τwave ;

Similarly, the positions of less probable nondimensional times (valleys) are given by Note also that those patterns become less pronounced when L is increased (see also section IV B later for the non diffusive case). There are two reasons for this: the first one is that diffusion effects are more important for larger travel distances; the second one is that the events (peaks or valleys) get mechanically closer to each other for increasing L, and sampling very rare events close to each other tend to blur the results.

τvalleys ≈ m + 1 2 × τwave . (25)
As the perturbed Poiseuille flow is expected to mix better than the basic Poiseuille flow, we would also expect the width of the RTD (measured using the mean absolute deviation) to be smaller in the perturbed case. Surprisingly, the width is larger in the perturbed configuration for the reference case.

Indeed, starting with σ 1,Poi ≃ 0.385 for the basic plane Poiseuille flow, mean absolute deviations of σ 1 = 0.431 for L = λ and σ 1 = 0.424 for L = 10λ are attained, corresponding to a significant increase of 12% and 10% respectively. This is due to the increase in maximum velocity at the center of the channel due to the nonlinear correction of the mean velocity profile. Indeed, increasing the maximum velocity while keeping the mean velocity constant automatically produces residence time values that are not only farther away from the mean residence time, but also highly probable (since associated with a high velocity region). The main peak located at τ = τwave leads to the same kind of effect, although less prominent given the much smaller height of the peak as seen in lin-lin scale in figure 9.

Finally, note that the width of the distributions is only slightly modified when increasing the traveling length L: the mean absolute deviation σ 1 decreases from σ 1 = 0.431 for L = λ to σ 1 = 0.424 for L = 10λ, which corresponds to a 1.6% decrease only. This is because the main peak is located at τ = τ min and is related to the maximal velocity at the center, see discussion later at the end of the section.

B. RTDs in the non diffusive limit

From what was learned before, is it possible to calculate the RTD by seeding particles at entrance of a section? Indeed, when computing the pdfs of the positions y i and the phases ϕ i corresponding to the residence time samples previously obtained in the reference case, it appears that the positions follow a distribution based on the mean velocity profile and that the phases obey an almost uniform distribution (see appendix B).

Therefore we propose to build a RTD in the limit Pe → ∞ following an Eulerian method, by seeding a large number M of non diffusing tracers at the entrance of a section of length L, as follows:

1. for each tracer i (1 ≤ i ≤ M ), the initial position is characterized by x i = 0, and an initial position y i and phase ϕ i chosen randomly, both following a uniform distribution;

2. the trajectory of each tracer is calculated until reaching the exit (x i = L);

3. the corresponding residence time τ i is weighted by the initial longitudinal velocity at entrance. This weighting reflects the fact that regions of high longitudinal velocity are more frequently crossed by traveling particles [START_REF] Poumaëre | Residence time distributions for in-line chaotic mixers[END_REF][START_REF] Saffman | A theory of dispersion in a porous medium[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF]. If the longitudinal velocity is negative, which corresponds to a point going backwards and initially located in one of the blue zones of figure 4(b), the point is omitted.

Note that we do not omit here all the points in the forbidden region (which is larger than the blue zone); however, the transverse velocities are very small in the whole forbidden region, which gives a negligible weight to all the forbidden points. Therefore we expect the resulting RTDs to be very close to those calculated with diffusion. In the case L = λ (figure 13 a), the RTDs are almost identical, except for the very long times where the artificial cut-off due to diffusion is no more present with the Eulerian method. Note indeed that with this latter method, the algebraic τ -3 tail is recovered, although modulated by the peak and valley pattern already addressed. In the case L = 10λ (figure 13 b), although the RTDs are identical for the most probable values (above 10 -1 ), the effect of diffusion is much more visible, not only for the very long times, but also in the vicinity of the main peak at τ wave . Indeed, without diffusion, a particle remains forever in the recirculation region where it was initially located. Since the residence time of all points in a recirculation region approaches τ wave as L is increased, (equation 23, see also figure 10), the peak at τ = τ wave gets narrower, which can be associated with unattainable residence time values; those zero probability regions are erased by diffusion, which allows particles to escape the recirculations on long distances, yielding intermediate residence values. Finally, at long residence times, the algebraic tail is also recovered (with a less visible pattern) with the Eulerian non diffusive method.

Overall, in terms of dispersion of the RTDs, the two methods lead to very similar values, with a mean absolute dispersion σ 1 within a 1% difference only, both for L = λ and L = 10 λ. This is because the most probable values (short times and peak at τ = τ wave ) are identical with both methods. Indeed, the RTDs obtained in figure 13 are almost identical when plotted in lin-lin scale (not shown here).

C. Variation of physical parameters

Up to now we have only considered the so-called reference case, α = 2 and Re = 10 000. We found that in that case, the TS waves tend to increase the width of the RTDs compared to the unperturbed Poiseuille flow, and also uncovered the existence of a peculiar pattern at long residence times. We now vary the wavenumber α and Reynolds number Re in order to explore the whole range of parameters of instability.

We begin with two sets of numerical experiments using the Lagrangian method: in the first one, the Reynolds number is kept fixed and α is varied; reversely, in the second one, the wavenumber is kept fixed, and the Reynolds number is varied. The results are displayed in figure 14.

We first address figure 14 to α max ≃ 2.18; those values correspond to the lower and upper limits of the neutral curve at this Reynolds number. What is most noticeable is the shift of the position of the peak corresponding to the TS waves towards lower residence time values when α is increased. Indeed, the peak occurs at τwave = U m /c nl , with U m = 1 for the whole study, and the nonlinear phase speed c nl is an increasing function of the wave number, as can be derived from figure 3(b).

Another characteristic of importance is the amplitude of the main peak, which increases with the wavenumber. This is a direct consequence of the fact that the total energy of the perturbation is also an increasing function of the wavenumber: when the total energy increases, the propagating vortices also proportionally grow in size in the periodic region [0 ≤ x ′ ≤ λ; -h/2 ≤ y ≤ h/2)] (figure 4 b), all the more so as the wave length λ decreases when increasing α. Thinking in Lagrangian terms, this means that a given particle will be more likely trapped in a recirculation region, associated with a residence time τ ≃ τ wave , which explains the higher peaks when α is increased.

Those changes are still present but much less pronounced in the case when the wavenumber α is kept constant, and the Reynolds number is varied, as shown in figure 14 b, where α = 2 and the Reynolds number is varied between 8 000 and 12 000: this is because the characteristics of the TS waves change much less in that case, in terms of nonlinear phase speed c nl and total energy.

Finally, note that the peculiar pattern of peaks and valleys observed in the reference case at small residence lengths L (and less visible when L is increased, as explained before) is in practice nearly always observed whatever the parameters α and Re taken inside the neutral curve. Only at very low wave energy, slightly above the bottom neutral curve, when no negative longitudinal velocities exist, can the pattern disappear. An example of the RTDs found in such a rare case is shown in appendix C for L = λ and L = 10λ.

In terms of width of distributions, a map of values of the mean absolute deviation σ 1 , calculated for L = λ with the Eulerian method (so as to compare with the values already obtained with the Lagrangian method and validate the method with more cases) is shown in figure 15(a).

As can be seen, the width of the distribution always increases when going away from the neutral curve. This is surprising from what was previously shown: because the peak at τ = τ wave increases and gets closer to τ m as α or Re increase, one could imagine that σ 1 would consequently decrease. However, the highest peak in the distributions is the one at the minimum residence time, corresponding to the maximum velocity at the center. A map of the mean flow correction at the center is given in figure 15(b): as can be seen, the velocity at the center increases away from the neutral curve, so that the minimum time τ min decreases, and curves are very similar in the lower part of the instability region, when the peak at τ = τ wave is small, while the effect of this peak is more visible in the upper part.

V. SUMMARY AND CONCLUSION

In this paper, the flow resulting from the nonlinear saturation of a linearly unstable Tollmien-Schlichting wave in the plane Poiseuille flow has been studied in terms of residence time distributions. A method for computing the RTD of this non chaotic flow, making use of diffusion as a way for a single particle to visit the entire cross-section of the channel has been outlined and validated on the simple case of the basic parabolic plane Poiseuille flow.

Since the choice of the travel distance was arbitrary, its impact on the RTD of the flow has first been assessed, assessment which at the same time brought to light the main features of the RTD and their links with the physical properties of the wave, which can be summarized as follows:

(i) For relatively small travel distances, the existence of regions of the flow where the longitudinal velocity is negative leads to a peak-valley pattern situated on the RTD's tail. This pattern indicates that some long residence times are much less likely to be encountered than others, entirely because of the steadily propagating vortices.

(ii) For longer travel distances, a main peak appears on the RTD, precisely located on the wave residence time. This peak is the consequence of particles being trapped and entrained in the propagating vortices, which travel downstream at constant speed.

(iii) No matter what the travel distance is, the distortion of the mean velocity profile decreases the minimum residence time (keeping the same mean residence time) compared to the basic plane Poiseuille flow.

This distortion of the mean velocity profile has been shown to be the main reason for increased dispersion of the RTD. Indeed, by varying the wavenumber of the disturbance, at constant Reynolds, and computing the dispersion of the RTD resulting from the corresponding saturated nonlinear Tollmien-Schlichting wave, good agreement has been found between the dispersion of the RTD -quantified by the mean absolute deviation -and the maximum velocity attained at the center of the channel. It is thus interesting to note that despite having better mixing capacities than the basic plane Poiseuille flow, saturated PPF suffers the counterpart of having increased dispersion of residence time values, which in some applications is of critical importance.

The Lagrangian method introduced in this article for the calculation of RTDs could easily be applied to any kind of open flow. For example, intermediate equilibrium states of the channel flow (3D disturbances, modulated Tollmien-Schlichting waves, etc.) could be considered up until complete destabilisation to turbulent state, in order to observe how the dispersion of residence time values evolves from the laminar value to the presumably very low dispersion attained with the turbulent velocity field.

FIG. 1 .

 1 FIG. 1. Schematic of an open flow with fluid trajectories, and corresponding RTD, with τm the mean residence time, and τmin and τmax respectively the minimum and maximum residence times.

FIG. 2 .

 2 FIG. 2. Channel geometry and parabolic plane Poiseuille flow profile.

FIG. 3 .

 3 FIG. 3. Numerical parameters for both figures: Nc = 65 and N h = 5 for the nonlinear computations. (a) Growth rate contours of the most unstable eigenmode of the linear eigenvalue problem. Inside the neutral curve (black solid line), growth rates are positive. Contour values shown are -0.12 to 0.0 in steps of 0.02, and 0.0 to 0.01 in steps of 0.002. (b) Phase speed contours for the linear (dashed lines) and nonlinear (solid lines) Tollmien-Schlichting waves. Contour values shown are 0.26 to 0.50 for linear phase speeds and 0.32 to 0.64 for nonlinear phase speeds, both in steps of 0.02.

FIG. 4 .

 4 FIG. 4. (a): Parabolic Poiseuille flow (dotted line) and distorted mean flow (solid line); (b): streamlines of the total flow for the saturated nonlinear TS wave at Re = 10 000 and α = 2 in the frame comoving with the wave (therefore in the (x ′ , y) plane). The small near-wall regions inside the blue curves display negative longitudinal velocities measured in the laboratory frame. Numerical parameters for the computation of nonlinear TS wave are: Nc = 65 and N h = 5.

Pe = 10 3

 3 FIG. 7. RTDs for the plane Poiseuille flow for different Péclet numbers. The size of the samples for each Péclet number is M = 10 8 . The travel distance is L = λ ref = π, except for the dashed line for which L = 10 λ ref = 10 π. The black solid line represents the analytical RTD in the non diffusive limit. (a) log-log scale; (b) lin-lin scale, with a focus on the lowest residence time values, indicated schematically on panel (a). Note that the tails of cases (Pe = 10 5 , L = π h) and (Pe = 10 6 , L = 10 πh), which correspond to the same of Pe h/L, are almost superimposed on panel (a); a similar remark can be made on panel (b) for cases (Pe = 10 6 , L = π h) and (Pe = 10 6 , L = 10 πh), which correspond to the same value of Pe L/h and for which the curves have a similar behavior at short times.

Pe = 10 3 ,

 3 FIG. 8. (a) Evolution of the mean absolute deviation whit the travel distance in the plane Poiseuille flow, for different Péclet numbers. (b) Evolution of the mean absolute deviation along the trajectory, for different Péclet numbers and travel distances. These diagrams correspond to the RTDs of figure 7.

FIG. 9 .

 9 FIG. 9. Normalized RTDs of the reference flow (solid line) in lin-lin scale; left: L = λ; right: L = 10λ. The sample corresponds to a single trajectory of total longitudinal length Ltot = 12. 10 6 ×λ. The dotted line represents the theoretical RTD of the basic Poiseuille flow (equation 21): the minimum residence time is clearly smaller in the reference case. The peak corresponding to the nondimensional wave residence time τwave is more pronounced for L = 10λ.

FIG. 10 .

 10 FIG.10. Colour maps of residence time values in the (ϕ, y) plane for the reference case α = 2 and Re = 10 000, built as follows: for a given trajectory, the position of the particle entering the j-th element at (j -1) × L is located by its coordinates (ϕj-1, yj-1) and associated to the residence time τj, as defined in figure6. The colour maps are obtained as a 2D histogram (250 × 250 uniform bins) of all the points of the long trajectory used in figure 11, weighted by the residence time values, and normalized by the mean residence time τm(L) = L/Um in each case. The colour scale is identical for all figures; the colour bar on the right is indicative and not to scale. The white regions correspond to a forbidden zone empty of points. The streamlines of figure 4(b) are reproduced on the last panel to emphasize the connection between the different regions of the flow and the resulting residence time properties.

FIG. 11 .FIG. 12 .

 1112 FIG. 11. Normalized RTDs of the perturbed flow for 4 travel distances from L = λ to L = 10λ, in log-log scale, for the reference case, α = 2 and Re = 10 000. All RTDs are plotted as histograms constituted of 500 logarithmically distributed bins. The sample corresponds to a single trajectory of total longitudinal length Ltot = 1.2 × 10 7 × λ. The dotted line represents the analytical RTD of the theoretical plane Poiseuille flow.

  This is indeed what is observed in figure11 a. When L = n × λ with n ≥ 2, as in figures 11(b-d), the mean traveling time is τ m = n × λ/c nl , and the peaks are located at τpeak ≈ m/n × τwave , with the minima at τvalleys ≈ (m + 1/2) /n × τwave . This explains why two events (peaks or valleys) get mechanically closer to each other for increasing L, as shown in the figure.

FIG. 13 .

 13 FIG. 13. RTDs of the reference flow (α = Re = 10 000), for the Eulerian method (dark colours) and for the Lagrangian method with diffusion (Pe = 10 5 , clear colours), for (a) L = λ and (b) L = 10λ. The diffusive RTDs are the same as in figures 11 (a) and (d). The dotted line is the RTD of the Poiseuille flow.

FIG. 14 .

 14 FIG. 14.RTDs for L = 10λ resulting from nonlinear Tollmien-Schlichting waves for (a) fixed Reynolds number Re = 10 000 and varying wavenumber, from αmin = 1.76 to αmax = 2.19 in steps of 0.04; (b) fixed α = 2, and varying Reynolds number, from Re = 8000 to Re = 11500 in steps of 500. Main plots are in log-log scale, while subplots focusing on the peaks are in lin-lin scale. At the top right of each plot are schematized the ranges considered throughout the neutral curve.

  FIG. 15. (a) Levels of mean absolute deviation σ1, from 0.385 (plane Poiseuille flow value) to 0.445 in steps of 0.002, for L = λ. (b) Levels of mean flow correction at the center of the channel, from 0.0175 to 0.0675 in steps of 0.005.

FIG. 17 .

 17 FIG. 17. RTDs of the perturbed flow for (a) L = 1λ and (b) L = 10λ normalized by the mean residence time τm = L/Um. The physical parameters are here Re = 12 000 and α = 1.6475, the resulting nonlinear wavespeed is c nl = 0.331.

Acknowlegments

Support from the Pôle de modélisation et de calcul en sciences de l'ingénieur et de l'information (PMCS2I) of École centrale de Lyon is gratefully acknowledged for providing access to high-performance computing resources. For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present manuscript and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

Appendix A: Calculation of the RTD and its mean absolute deviation for the plane Poiseuille flow

The probability for a fluid particle to cross a given cross-section of the channel at this distance must be proportional to the velocity field at this point. If we take y as the distance from the centre of the channel, the velocity profile of the plane Poiseuille flow is

Normalizing the velocity field in order to have an integral unity across half the channel width, we obtain the probability density

Let f Poi (τ ) be the probability density of the time needed to cross an element of length L. Thus the probability for a fluid particle to travel a length L in a time between τ and τ + dτ is equal to the probability for the particle to begin the crossing between ordinates y and y + dy, that is:

We denote by τ m the mean time to travel the distance L; we have:

which implies:

This can be differentiated into

Combining (A2), (A3), and (A6), we finally obtain:

Note that despite being singular in the limit τ → τ min = 2 3 τ m , the following classical properties hold: We can now calculate the mean absolute deviation of this distribution:

First, by setting u = τ /τ m , we obtain

Then, by setting x 2 = 1 -2/3u, so that xdx = 1/3du/u 2 and 1/3u = (1x 2 )/2, we obtain

and finally:

Appendix B: Pdfs of position and phase