# Residence time distributions in unstable channel flow

Nelson Poumaëre, Benoît Pier, Florence Raynal

► To cite this version:

Nelson Poumaëre, Benoît Pier, Florence Raynal. Residence time distributions in unstable channel flow. Physical Review Fluids, 2024, 9 (10), pp.104501. 10.1103/PhysRevFluids.9.104501. hal-04296888v2

HAL Id: hal-04296888

https://hal.science/hal-04296888v2

Submitted on 18 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Residence time distributions in unstable channel flow

Nelson Poumaëre, Benoît Pier, and Florence Raynal

Laboratoire de Mécanique des Fluides et d'Acoustique, CNRS, École centrale de Lyon,
INSA Lyon, Université Claude Bernard Lyon 1, 69134 Écully Cédex, France

(Dated: October 18, 2024)

The transport properties of the plane Poiseuille flow in which a two-dimensional, nonlinearly saturated Tollmien-Schlichting wave is propagating are studied in terms of residence time distributions (RTDs). First, a method for computing RTDs in any type of open flows is developed, making use of a single trajectory over a long period of time, with a controlled level of diffusion. With this method, RTDs of this perturbed flow are computed, along with a quantitative measure of their dispersion through the mean absolute deviation. Depending on the travel distance, RTDs display two kinds of pattern. For short travel distances, a pattern of peaks and valleys is observed for long residence times, originating in regions of negative streamwise velocity produced by the wave. For longer travel distances, a large probability peak is observed at $t = \tau_{wave}$, the time needed for the wave to travel one section downstream. This peak is attributed to the cat's eyes pattern characteristic of this type of traveling wave. It is shown that the increased dispersion of the RTD is mainly due to the nonlinear correction of the mean velocity profile.

## I. INTRODUCTION

Stirring in laminar flows is an essential step in all mixing processes [1]. 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 [2, 3], 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 [4] 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 [7], Trachsel et al. [8]). 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. [9].

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 $\tau = t_{out} - t_{in}$. As shown by Spalding [10] and

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic of an open flow with fluid trajectories, and corresponding RTD, with $\tau_m$ the mean residence time, and $\tau_{min}$ and $\tau_{max}$ respectively the minimum and maximum residence times.</figcaption></figure>

Danckwerts [11] in steady conditions, the mean residence time $\tau_m$ is given by

$$ \tau_m = \int_{\tau_{min}}^{\tau_{max}} \tau f(\tau) d\tau = \frac{V}{Q}, \qquad (1) $$

where $f(\tau)$ 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 $\tau_{min}$ and $\tau_{max}$ respectively correspond to the fastest and slowest trajectories attainable in the system. Note that equation 1 fails for flows with closed fixed recirculations, as can be seen for instance behind an obstacle or a step. Indeed, without diffusion, those recirculations would not be reached by particles seeded at the entrance of the system. Therefore in this latter case, the volume $V$ should be replaced by a smaller volume $V^*$ for which the volume occupied by the recirculations has been removed.
Ideally, a mixer should mimic as closely as possible the plug-flow model, where all particles cross the system in the same time $\tau = \tau_{min} = \tau_{max}$, and for which the RTD is essentially a Dirac pulse. In practice however, the maximum residence time $\tau_{max}$ is generally infinite, which implies that the RTD $f(\tau)$ 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. [13]), the situation is very different when considering flows with walls, since the decay is thus algebraic, with a $\tau^{-3}$ tail [12], 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. [14]).

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 [5], 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. [15] who studied a two-dimensional confined convection flow, or the works by Meunier [16], along with Lefranc et al. [17], taking advantage of triadic resonances in a tilted tank of water to efficiently mix with limited stress, or also Qu et al. [18] 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. [19] studied how laminar vortex shedding can affect temperature and concentration fields in the vicinity of direct contact membrane distillation channels. Venditti et al. [20] 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. [21, 22] 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 [23], Thomas [24] and Grosch & Salwen [25]), it is Orszag [26] 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 $\tilde{v}$ of an infinitesimal perturbation written in complex form as

$$ \mathbf{u}(x, y, t) = \tilde{\mathbf{u}}(y) \exp [i(\alpha x - \omega t)] + c.c., \quad (2a) $$

$$ p(x, y, t) = \tilde{p}(y) \exp [i(\alpha x - \omega t)] + c.c., \quad (2b) $$

where $\alpha$ denotes the real wavenumber of the perturbation, $\omega = \omega_r + i\omega_i$ its complex frequency, and c.c. the complex conjugate. For a given setting of Re and $\alpha$, the linear eigenvalue problem arising from the Orr-Sommerfeld equation is solved, yielding a discrete spectrum of eigenvalues $\omega$ and the associated eigenfunctions $\tilde{u}(y)$ and $\tilde{p}(y)$. If the growth rate $\omega_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 $\alpha_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 Tollmien-Schlichting (TS) wave. Note that the TS mode exhibits sinuous symmetry. In the $(Re, \alpha)$-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 \approx 1500)$ the critical Reynolds number (Davies & White [27]). 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

$$ \mathbf{u}(x, y, t) = \sum_{n} \mathbf{u}^{(n)}(y, t) \exp [in\alpha x]. \qquad (3) $$

Early works (Meksyn & Stuart [29], Stuart [30]) 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 [31] and Watson [32] 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. [33] or Jiménez [34] to study different solution branches bifurcating from basic plane Poiseuille flow. Ehrenstein & Koch [35], using Fourier expansion and pseudo-arclength continuation, thoroughly computed two- and 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 the system is ergodic, and 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 \times L$ with $M \gg 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 [12, 36]. 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 $\tau_m$ defined in equation (1). 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 [37], and more recently Dentz et al. [38] and Puyguiraud et al. [39]. 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. [40]). 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 [10], 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,

$$ \partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} = \nu \Delta \mathbf{u} - \frac{1}{\rho} \nabla p, \quad (4a) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (4b) $$

where $\rho$ and $\nu$ 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 $\mathbf{u}$ and pressure $p$ then depend on position $\mathbf{x} = x\mathbf{e}_x + y\mathbf{e}_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 = \pm h/2$ (see figure 2).

<figure><img src="image_3.png" /><figcaption>FIG. 2. Channel geometry and parabolic plane Poiseuille flow profile.</figcaption></figure>

The steady base Poiseuille flow solution follows the classical parabolic profile

$$ U(y) \mathbf{e}_x = \frac{3}{2} U_m \left( 1 - \frac{4y^2}{h^2} \right) \mathbf{e}_x, \quad (5) $$

with $U_m = \frac{h^2}{12\nu}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$:

$$ \mathrm{Re} = \frac{hU_m}{\nu}. \quad (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 quantities as

$$ \mathbf{u}_{\mathrm{tot}}(x, y, t) = U(y) \mathbf{e}_x + \mathbf{u}_{\mathrm{pert}}(x, y, t), \quad (7) $$

$$ p_{\mathrm{tot}}(x, y, t) = -Gx + p_{\mathrm{pert}}(x, y, t), \quad (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 $\mathbf{u}_{\mathrm{pert}}$ and $p_{\mathrm{pert}}$ are obtained as

$$ \partial_t \mathbf{u}_{\mathrm{pert}} = -\mathcal{L}(\mathbf{u}_{\mathrm{pert}}) - \mathcal{N}(\mathbf{u}_{\mathrm{pert}}) + \nu \Delta \mathbf{u}_{\mathrm{pert}} - \frac{1}{\rho} \nabla p_{\mathrm{pert}}, \quad (9a) $$

$$ \nabla \cdot \mathbf{u}_{\mathrm{pert}} = 0, \quad (9b) $$
where $\mathcal{L}(\mathbf{u}_{\text{pert}}) = U'v_{\text{pert}}\mathbf{e}_x + U\partial_x\mathbf{u}_{\text{pert}}$ is the linear term accounting for the advection by the base flow and $\mathcal{N}(\mathbf{u}_{\text{pert}}) = (\mathbf{u}_{\text{pert}} \cdot \nabla)\mathbf{u}_{\text{pert}}$ is the nonlinear self-advection term. For the rest of this paper, we will only consider the perturbation velocity and pressure fields $\mathbf{u}_{\text{pert}}$ and $p_{\text{pert}}$ and, for simplicity, omit their subscript.

When considering small-amplitude perturbations, the quadratic term $\mathcal{N}(\mathbf{u})$ is neglected and the flow fields are sought in normal-mode form as

$$ \mathbf{u}(x, y, t) = \tilde{\mathbf{u}}(y) \exp [i(\alpha x - \omega t)], \quad (10a) $$

$$ p(x, y, t) = \tilde{p}(y) \exp [i(\alpha x - \omega t)], \quad (10b) $$

where $\alpha$ is the streamwise wavenumber and $\omega$ 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 $\omega$-plane for a given wavenumber $\alpha$ (Schmid & Henningson [41]).

A positive imaginary part $\omega_i$ of the frequency associated by the linear dispersion relation to a given wavenumber $\alpha$ signals a linear instability. An initial spatially periodic perturbation of wavenumber $\alpha$ is then exponentially amplified in time with growth rate $\omega_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

$$ \mathbf{u}(x, y, t) = \sum_{n} \mathbf{u}^{(n)}(y, t) \exp [in\alpha x], \quad (11a) $$

$$ p(x, y, t) = \sum_{n} p^{(n)}(y, t) \exp [in\alpha x], \quad (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 $\mathbf{u}^{(n)}(y, t)$ and $p^{(n)}(y, t)$. This is equivalent to considering the dynamics in a finite domain of length $2\pi/\alpha$ 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 [42] 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 \le y \le h/2$. Thus, the differential operator $\partial_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 $\tilde{\mathbf{u}}(y)$ and $\tilde{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 $\mathbf{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 [43], improved by Raspo et al. [44]. More details of the general method can be found in [42, 45]. The numerical implementation in C++ is based on the ‘home-spun’ PackstaB library [46]. The required numerical effort is reduced by explicitly taking into account symmetry properties prevailing for sinusoidal 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 $(\text{Re}, \alpha)$-plane, by solving the algebraic eigenvalue problem and monitoring the most unstable eigenvalue of the spectrum obtained for each choice of $\text{Re}$ and $\alpha$. Figure 3(a) shows isolines of constant growth rate $\omega_i$ along with the neutral curve. Note that the most unstable eigenmode is always of sinusoidal symmetry. The associated phase velocity $\omega_r/\alpha$ is plotted as dashed lines in Figure 3(b). The region inside the neutral curve $\omega_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 instabilities, that would lead to
<figure><img src="image_4.png" /><figcaption>FIG. 3. Numerical parameters for both figures: $N_c = 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.</figcaption></figure>

more complex dynamics. Thus, in this finite-amplitude regime, the temporal evolution of the different Fourier components in (11a)(11b), is of the form $\mathbf{u}^{(n)}(y, t) = \tilde{\mathbf{u}}^{(n)}(y) \exp(-in\omega^{nl}t)$, for some real-valued nonlinear frequency $\omega^{nl}$. By monitoring the different components during the temporal integration and computing the frequency $\omega^{nl}$, it is thus verified that, after convergence, the perturbations may indeed be written as a steadily traveling nonlinear wave

$$ \mathbf{u}(x, y, t) = \sum_{|n|<N_h} \tilde{\mathbf{u}}^{(n)}(y) \exp [in(\alpha x - \omega^{nl}t)]. \quad (12) $$

This nonlinear Tollmien-Schlichting wave travels with a phase velocity $c^{nl} = \omega^{nl}/\alpha$; 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 $\lambda = 2\pi/\alpha$, 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 $\phi$ with

$$ \phi = \alpha x - \omega^{nl} t \quad \mod 2\pi \text{ with } 0 \le \phi \le 2\pi. \quad (13) $$

Similarly, we define the periodic longitudinal position $x'$ in the moving reference frame as

$$ x' = \phi/\alpha = x - c^{nl} t \quad \mod \lambda \text{ with } 0 \le x' \le \lambda. \quad (14) $$

**D. A reference case: $\alpha = 2$ and $\text{Re} = 10000$**

As an example, the result of a simulation for $\text{Re} = 10000$ and $\alpha = 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 $\omega = 0.74967 + 0.00670i$. For this case we have $\lambda = \lambda_{\text{ref}} = 2\pi/\alpha = \pi$, the nonlinear frequency is $\omega^{nl} = 1.17$, yielding a nonlinear wavespeed $c^{nl} = 0.587$ and a time period $T = 2\pi/\omega^{nl} = 5.36$. The longitudinal extent for the contour map is the spatial pseudo-period $\lambda = \pi$. 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
<figure><img src="image_5.png" /><figcaption>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 = 10000$ and $\alpha = 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: $N_c = 65$ and $N_h = 5$.</figcaption></figure>

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.

In figure 5, we represent snapshots of the trajectories of 2500 particles over one period of the flow. As was described by Weiss [47] 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 shown in figure 4 b. In particular, particles traveling slowly near the separatrix are easily distinguishable at time $t = T$.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Snapshots of the positions of 2500 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 = 10000$ and $\alpha = 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.</figcaption></figure>
### III. COMPUTATION OF RESIDENCE TIME DISTRIBUTIONS

In this section we develop a general method for the calculation of RTDs, that could in principle be available for any type of in-line flow or mixer where no fixed recirculation is present. Note that in our case the recirculations shown in figure 4 b are not fixed since they move a the speed of the wave. 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; we allow the particle to diffuse so as to sample the entire cross-section of the flow and render the system ergodic. The advection of a single diffusive particle can be modeled with a generalized Langevin equation

$$ \frac{d\mathbf{x}}{dt} = \mathbf{u}(\mathbf{x}, t) + \zeta(t), \qquad (15) $$

where $\zeta(t)$, a Gaussian decorrelated process such that $\langle \zeta_i(t)\zeta_j(t') \rangle = 2D\delta_{ij}\delta(t-t')$, with D the diffusion coefficient, $\delta_{ij}$ the Krönecker symbol and $\delta$ the Dirac delta function, models the molecular diffusion [48]. 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$,

$$ \mathrm{Pe} = \frac{hU_m}{D}. \qquad (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 $\phi_{j-1} = \alpha \times (j-1)L - \omega^{nl}t_{j-1}$; then we calculate the residence time $\tau_j = t_j - t_{j-1}$ in this element as the time it takes to reach the exit plane $x = jL$ 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 = \sum_{i=1}^{j} \tau_i$. In practice, we follow the particle on a very long time so as to obtain a large sample of M data $(\tau_j, y_j, \phi_j)_{1 \le j \le M}$. The process is illustrated in figure 6. The RTD can thus be obtained as the histogram of the resulting M residence time values $\tau_j$; the

<figure><img src="image_7.png" /><figcaption>FIG. 6. Schematic of the partitioning of the flow in sections of length L, along with a particle trajectory and retained data $(\tau_j, y_j, \phi_j)$.</figcaption></figure>

other parameters $(y_j, \phi_j)$ will be used for other purposes, see section IV. In practice, we store these data for $L = \lambda$. The case $L = n \times \lambda$ ($2 \le n \le 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 $(\tau_j)_{1 \le j \le M}$ of residence times is measured through the mean absolute deviation

$$ \sigma_1 = \frac{1}{M} \sum_{j=1}^{M} \left| \frac{\tau_j}{\tau_m} - 1 \right|, \qquad (17) $$

where $\tau_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_{\text{conv}} = L/U_m$. This means that (i) the characteristic longitudinal diffusive time $t_{\text{diff},x} = L^2/D$ must be large compared
to $t_{\text{conv}}$, i.e. $t_{\text{diff},x} \gg t_{\text{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_{\text{diff},y} = h^2/D$ must be large compared to $t_{\text{conv}}$, i.e. $t_{\text{diff},y} \gg t_{\text{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 \times L$ in the $x$-direction, associated to a typical lapse of time $M \times t_{\text{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_{\text{diff},y} \ll Mt_{\text{conv}}$. Conditions (i)—(iii) can be rewritten as:

$$Pe \gg h/L, \qquad (18)$$

$$Pe \gg L/h, \qquad (19)$$

$$M \gg Pe h/L. \qquad (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 (19) along with condition (20) is sufficient in practice. Note finally that condition (20) can be rewritten $L_{\text{tot}} \gg Pe h$, where $L_{\text{tot}} = M \times 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 $(10L, M/10)$, since both cases correspond to the same total longitudinal length traveled $L_{\text{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. [49] or Mezić et al. [50] 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 $\lambda = 2\pi/\alpha$, a natural possibility is to use $L = n \times \lambda$, 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 (5). 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 \to \infty$)

$$f_{\text{Poi}}(\tau) = \frac{1}{3} \frac{\tau_m^2}{\tau^3} \left(1 - \frac{2\tau_m}{3\tau}\right)^{-1/2}, \qquad (21)$$

with the corresponding analytical value of the mean absolute deviation, see appendix A,

$$\sigma_{1,\text{Poi}} = 2/(3\sqrt{3}) \approx 0.385. \qquad (22)$$

Note that, unlike the perturbed flows considered in the present study, the parabolic flow field has no natural length; thus equations (21, 22) 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 IID with wavenumber $\alpha = 2$, and equivalently $\lambda_{\text{ref}} = \pi$.

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 $\sigma_1$: a total longitudinal distance $L_{\text{tot}} = M \times \pi$ with $M = 10^8$ is fairly sufficient in practice. Then we compute the resulting RTDs (normalized by their corresponding mean residence time $\tau_m$) for the different Péclet numbers Pe and travel distances $L$, and calculate the associated mean absolute deviations $\sigma_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 (19): 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 longer possible. This results in a cut-off of the algebraic $\tau^{-3}$ tail obtained without diffusion [14]. When increasing the Péclet number Pe for the fixed travel distance $L = \lambda$, diffusion effects are less important, and the cut-off is shifted toward higher residence times as observed in figure 7 (a). Reversely, when
<figure><img src="image_8.png" /><figcaption>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 = \lambda_{\text{ref}} = \pi$, except for the dashed line for which $L = 10\lambda_{\text{ref}} = 10\pi$. 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 = \pi h$) and ($Pe = 10^6, L = 10\pi h$), which correspond to the same value 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 = \pi h$) and ($Pe = 10^6, L = 10\pi h$), which correspond to the same value of $Pe L/h$ and for which the curves have a similar behavior at short times.</figcaption></figure>

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 = \pi$ and $L = 10\pi$. Note also that the tails of the RTDs of the cases ($Pe = 10^4, L = \pi$) and ($Pe = 10^5, L = 10\pi$) are almost identical. Overall, when the quantity $Pe h/L$ increases, the cut-off is shifted toward higher residence times, as expressed in condition (19).

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, $\tau_{\min} = 2/3 \times \tau_m$, see equation (21). When $Pe$ is increased with a fixed travel distance, diffusion effects are less important, and the RTDs become more peaked around $\tau_{\min}$, as shown in the figure for $L = \lambda$. 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 = \lambda$ and $L = 10\lambda$. Note finally that the shape of the RTDs at small time is very similar for the two cases ($Pe = 10^5, L = 10\pi$) and ($Pe = 10^6, L = \pi$). Overall, the RTDs are more peaked on $\tau_{\min}$ when increasing the quantity $PeL/h$, as expressed in condition (18).

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 $\sigma_{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\pi$ with $1 \le n \le 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 $Pe h/L$ corresponding to the parameter involved in the most restrictive condition (equation 19). For small values of $Pe h/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 (20) 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$, ($\tau_i$)$_{1 \le i \le M}$, with $M_{\max} = 10^8$. The results are shown in figure 8 (b), where the computed values of $\sigma_1$ are plotted as a function of the non dimensional parameter $ML/hPe$ involved in equation (20). Note that we recover that for small values of $Pe h/L$,
<figure><img src="image_9.png" /><figcaption>FIG. 8. (a) Evolution of the mean absolute deviation with 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.</figcaption></figure>

as discussed before, the curves do not converge toward the correct theoretical value. This is the case for instance
for the parameters ($Pe = 10^4$, $L = 10\pi$), whilst the theoretical value is reached for the same Péclet number and
a smaller travel distance, i.e. for the couple ($Pe = 10^4$, $L = \pi$). However, although they do not converge toward
the same value, the two curves behave similarly with the parameter $ML/hPe$; indeed, for a given value of the
parameter, as explained before, the total distance $L_{tot} = M \times L$ is the same (with a sample 10 times smaller
for $L = 10\pi$ than for $L = \pi 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 = \pi$)
and ($Pe = 10^5$, $L = 10\pi$), 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 $ML/hPe \ge 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 \ge 300 \times Peh/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 \approx 300 Peh/L = 3 \times 10^7 h/L$.
For instance, for the reference case for which $\lambda = \pi$ (with $h = 1$), we need a sample of size $M \approx 10^7$.

IV. RESULTS

In this section, we first consider the reference case, the nonlinear TS wave introduced in section IID, 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 $\tau_{wave} = L/c^{nl}$ it takes for one wave train to travel the distance $L$ as the wave residence time; similarly, $\tau_{min}$ is the minimum of the residence time for the distance $L$. We introduce the nondimensional wave residence time $\hat{\tau}_{wave} = \tau_{wave}/\tau_m = U_m/c^{nl}$ and also $\hat{\tau}_{min} = \tau_{min}/\tau_m = U_m/U_{max}$. Note that both $\hat{\tau}_{wave}$ and $\hat{\tau}_{min}$ are independent on the travel distance $L$.
A. The reference case $\alpha = 2$ and Re = 10 000

We first begin with the reference case already described in paragraph II D, for which $\alpha = 2$, Re = 10 000, and $\lambda = 2\pi/\alpha = \pi$. 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 \times 10^6$ sections of length $L = \lambda$, resulting in a sample $(\tau_i, y_i, \phi_i)$ of size $M$.

<figure><img src="image_10.png" /><figcaption>FIG. 9. Normalized RTDs of the reference flow (solid line) in lin-lin scale; left: $L = \lambda$; right: $L = 10\lambda$. The sample corresponds to a single trajectory of total longitudinal length $L_{tot} = 12 \times 10^6 \times \lambda$. 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 $\hat{\tau}_{wave}$ is more pronounced for $L = 10\lambda$.</figcaption></figure>

Figure 9 shows the RTDs obtained for the travel distances $L = \lambda$ and $L = 10\lambda$. For the purpose of comparison, the RTD of the parabolic Poiseuille flow, obeying to equation (21), 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 $\hat{\tau}_{wave}$. The peak is much higher and more marked for the larger travel distance $L = 10\lambda$.

To explain this behavior, the samples are visualized in the $(\phi, y, \tau)$ 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 $\hat{\tau}_{wave}$ (green color). This behavior is more pronounced as $L$ is increased. Indeed, focus on the lower recirculation region in figure 10(a) ($L = \lambda$): 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 $\tau_{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\lambda, n \ge 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 $\tau_{wave} = n\lambda/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 $\lambda$. Therefore, the residence time in the laboratory frame is less than that of the central point, such that $(n-1) \times \lambda/c^{nl} \le \tau \le n \times \lambda/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 $\lambda$ in the moving frame, and the residence time is higher and verifies $n \times \lambda/c^{nl} \le \tau \le (n+1) \times \lambda/c^{nl}$. Overall, all points in the recirculation have a residence time such that $(n-1) \times \lambda/c^{nl} \le \tau \le (n+1) \times \lambda/c^{nl}$. Finally, the nondimensional residence time $\hat{\tau}$ of a point in the
<figure><img src="image_11.png" /><figcaption>FIG. 10. Colour maps of residence time values in the $(\phi, y)$ plane for the reference case $\alpha = 2$ and $\text{Re} = 10000$, built as follows: for a given trajectory, the position of the particle entering the $j$-th element at $(j-1) \times L$ is located by its coordinates $(\phi_{j-1}, y_{j-1})$ and associated to the residence time $\tau_j$, as defined in figure 6. The colour maps are obtained as a 2D histogram $(250 \times 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 $\tau_m(L) = L/U_m$ 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.</figcaption></figure>

recirculation is such that

$$ (1 - 1/n) \hat{\tau}_{\text{wave}} \le \hat{\tau} \le (1 + 1/n) \hat{\tau}_{\text{wave}}, \qquad (23) $$

which approach $\hat{\tau}_{\text{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 = \lambda$ to $L = 10\lambda$ in log-log scale: although we globally recover the $t^{-3}$ behavior and its cut-off due to diffusion [14, 38, 39], 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. Note finally that this does not break the ergodicity necessary to the validity of the Lagrangian method: these regions of the flow are nevertheless visited, but these events do not correspond to a crossing of a line $x = n \times L$, and therefore are only forbidden from the Poincaré section point of view inherent to residence times calculated with a Lagrangian method.

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
<figure><img src="image_12.png" /><figcaption>FIG. 11. Normalized RTDs of the perturbed flow for 4 travel distances from $L = \lambda$ to $L = 10\lambda$, in log-log scale, for the reference case, $\alpha = 2$ and Re = 10000. All RTDs are plotted as histograms constituted of 500 logarithmically distributed bins. The sample corresponds to a single trajectory of total longitudinal length $L_{tot} = 1.2 \times 10^7 \times \lambda$. The dotted line represents the analytical RTD of the theoretical plane Poiseuille flow.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. Blow-up of figure 4b near the bottom wall (streamlines of the reference case, $\alpha = 2$ and Re = 10000, in the referential of the wave). 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' = \phi/\alpha, y$) plane when entering the $j$th element at $x = (j - 1) \times \lambda$ (figure 6) and can be seen as a Poincaré section of the diffusing trajectory.</figcaption></figure>

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 $\lambda/c^{nl}$ in the moving frame. 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.

Let us now focus on figure 12: a particle starting outside the forbidden region is always allowed to travel a distance equal to $n \times \lambda$ 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 $\tau \approx m \times \lambda/c^{nl}$, where $m$ is an integer which is at least such that $m \ge 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) \times \lambda$, which makes this distance in the moving frame less probable. This corresponds to a rare travel time $\tau \approx (m + 1/2) \times \lambda/c^{nl}$.

Let us come back now to figure 11: in the case $L = \lambda$, the mean travel time is $\tau_m = \lambda/c^{nl}$. Therefore the most probable nondimensional times are

$$ \hat{\tau}_{\text{peak}} \approx m \times \hat{\tau}_{\text{wave}}; \qquad (24) $$

Similarly, the positions of less probable nondimensional times (valleys) are given by

$$ \hat{\tau}_{\text{valleys}} \approx \left(m + \frac{1}{2}\right) \times \hat{\tau}_{\text{wave}}. \qquad (25) $$

This is indeed what is observed in figure 11a. When $L = n \times \lambda$ with $n \ge 2$, as in figures 11(b–d), the mean traveling time is $\tau_m = n \times \lambda/c^{nl}$, and the peaks are located at $\hat{\tau}_{\text{peak}} \approx m/n \times \hat{\tau}_{\text{wave}}$, with the minima at $\hat{\tau}_{\text{valleys}} \approx (m + 1/2)/n \times \hat{\tau}_{\text{wave}}$. This explains why two events (peaks or valleys) get mechanically closer to each other for increasing $L$, as shown in the figure.

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.

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 $\sigma_{1,Poi} \approx 0.385$ for the basic plane Poiseuille flow, mean absolute deviations of $\sigma_1 = 0.431$ for $L = \lambda$ and $\sigma_1 = 0.424$ for $L = 10\lambda$ 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 $\hat{\tau} = \hat{\tau}_{\text{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 $\sigma_1$ decreases from $\sigma_1 = 0.431$ for $L = \lambda$ to $\sigma_1 = 0.424$ for $L = 10\lambda$, which corresponds to a 1.6% decrease only. This is because the main peak is located at $\tau = \tau_{\text{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 $\phi_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 \to \infty$ 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 \le i \le M$), the initial position is characterized by $x_i = 0$, and an initial position $y_i$ and phase $\phi_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 $\tau_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 [14, 37, 38]. 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.

<figure><img src="image_14.png" /><figcaption>FIG. 13. RTDs of the reference flow ($\alpha = 2$, Re = 10000), for the Eulerian method (dark colours) and for the Lagrangian method with diffusion (Pe = $10^5$, clear colours), for (a) $L = \lambda$ and (b) $L = 10\lambda$. The diffusive RTDs are the same as in figures 11 (a) and (d). The dotted line is the RTD of the Poiseuille flow.</figcaption></figure>

The resulting RTDs are shown for the reference case in figure 13, with (a) $L = \lambda$ and (b) $L = 10\lambda$. For sake of comparison the diffusive RTDs obtained by the Lagrangian method with Pe = $10^5$ (figures 11 (a) and (d)) are also plotted.

In the case $L = \lambda$ (figure 13a), 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 $\tau^{-3}$ tail is recovered, although modulated by the *peak and valley* pattern already addressed. In the case $L = 10\lambda$ (figure 13b), 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 $\tau_{\text{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 $\tau_{\text{wave}}$ as $L$ is increased, (equation 23, see also figure 10), the peak at $\tau = \tau_{\text{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 $\sigma_1$ within a 1% difference only, both for $L = \lambda$ and $L = 10\lambda$. This is because the most probable values (short times and peak at $\tau = \tau_{\text{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 Re and $\alpha$

Up to now we have only considered the so-called reference case, $\alpha = 2$ and Re = 10000. 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 $\alpha$ 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 $\alpha$ 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.

<figure><img src="image_15.png" /><figcaption>FIG. 14. RTDs for $L = 10\lambda$ resulting from nonlinear Tollmien-Schlichting waves for (a) fixed Reynolds number $Re = 10000$ and varying wavenumber, from $\alpha_{min} = 1.76$ to $\alpha_{max} = 2.19$ in steps of 0.04; (b) fixed $\alpha = 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.</figcaption></figure>

We first address figure 14(a), for which the Reynolds number is $Re = 10000$, and $\alpha$ varies from $\alpha_{min} \simeq 1.76$ to $\alpha_{max} \simeq 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 $\alpha$ is increased. Indeed, the peak occurs at $\hat{\tau}_{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 \le x' \le \lambda; -h/2 \le y \le h/2]$ (figure 4b), all the more so as the wave length $\lambda$ decreases when increasing $\alpha$. Thinking in Lagrangian terms, this means that a given particle will be more likely trapped in a recirculation region, associated with a residence time $\tau \simeq \tau_{wave}$, which explains the higher peaks when $\alpha$ is increased.

Those changes are still present but much less pronounced in the case when the wavenumber $\alpha$ is kept constant, and the Reynolds number is varied, as shown in figure 14b, where $\alpha = 2$ and the Reynolds number is varied between 8000 and 12000: 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 $\alpha$ 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 = \lambda$ and $L = 10\lambda$.

In terms of width of distributions, a map of values of the mean absolute deviation $\sigma_1$, calculated for $L = \lambda$ 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
<figure><img src="image_16.png" /><figcaption>FIG. 15. (a) Levels of mean absolute deviation $\sigma_1$, from 0.385 (plane Poiseuille flow value) to 0.445 in steps of 0.002, for $L = \lambda$. (b) Levels of mean flow correction at the center of the channel, from 0.0175 to 0.0675 in steps of 0.005.</figcaption></figure>

is surprising from what was previously shown: because the peak at $\tau = \tau_{\text{wave}}$ increases and gets closer to $\tau_m$ as $\alpha$ or Re increase, one could imagine that $\sigma_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 $\tau_{\min}$ decreases, and curves are very similar in the lower part of the instability region, when the peak at $\tau = \tau_{\text{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.

## Acknowledgments

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 all present versions of this manuscript up to this Author Accepted Manuscript.

## 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

$$U_{\text{Poi}}(y) = \frac{3}{2} U_m \left(1 - \frac{4y^2}{h^2}\right). \qquad (A1)$$

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

$$g_{\text{Poi}}(y) = \frac{3}{h} \left(1 - \frac{4y^2}{h^2}\right). \qquad (A2)$$

Let $f_{\text{Poi}}(\tau)$ 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 $\tau$ and $\tau + d\tau$ is equal to the probability for the particle to begin the crossing between ordinates $y$ and $y + dy$, that is:

$$f_{\text{Poi}}(\tau) \, d\tau = g_{\text{Poi}}(y) \, dy. \qquad (A3)$$

We denote by $\tau_m$ the mean time to travel the distance $L$; we have:

$$\tau = \frac{L}{U_{\text{Poi}}(y)} = \frac{2}{3} \frac{\tau_m}{1 - 4y^2/h^2}, \qquad (A4)$$

which implies:

$$y = \frac{h}{2} \sqrt{1 - \frac{2 \tau_m}{3 \tau}}. \qquad (A5)$$

This can be differentiated into

$$\frac{dy}{d\tau} = \frac{h}{6} \frac{\tau_m}{\tau^2} \left(1 - \frac{2\tau_m}{3\tau}\right)^{-1/2}. \qquad (A6)$$

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

$$f_{\text{Poi}}(\tau) = \frac{1}{3} \frac{\tau_m^2}{\tau^3} \left(1 - \frac{2\tau_m}{3\tau}\right)^{-1/2}, \qquad (A7)$$

Note that despite being singular in the limit $\tau \to \tau_{\min} = \frac{2}{3}\tau_m$, the following classical properties hold:

$$\int_{\tau_{\min}}^{+\infty} f_{\text{Poi}}(\tau) \, d\tau = 1, \quad \text{and} \quad \int_{\tau_{\min}}^{+\infty} \tau f_{\text{Poi}}(\tau) \, d\tau = \tau_m. \qquad (A8)$$
We can now calculate the mean absolute deviation of this distribution:

$$ \sigma_{1,Poi} = \int_{\tau_{min}}^{+\infty} \left| \frac{\tau}{\tau_m} - 1 \right| f_{Poi}(\tau) d\tau = \int_{\tau_{min}}^{\tau_m} \left( 1 - \frac{\tau}{\tau_m} \right) f_{Poi}(\tau) d\tau + \int_{\tau_m}^{+\infty} \left( \frac{\tau}{\tau_m} - 1 \right) f_{Poi}(\tau) d\tau. \quad (A9) $$

First, by setting $u = \tau/\tau_m$, we obtain

$$ \sigma_{1,Poi} = \int_{2/3}^{1} \frac{1-u}{3u^3} \left(1 - \frac{2}{3u}\right)^{-1/2} du + \int_{1}^{+\infty} \frac{u-1}{3u^3} \left(1 - \frac{2}{3u}\right)^{-1/2} du, \quad (A10) $$

$$ = \int_{2/3}^{1} \left( \frac{1}{3u} - \frac{1}{3} \right) \left( 1 - \frac{2}{3u} \right)^{-1/2} \frac{du}{u^2} + \int_{1}^{+\infty} \left( \frac{1}{3} - \frac{1}{3u} \right) \left( 1 - \frac{2}{3u} \right)^{-1/2} \frac{du}{u^2}. \quad (A11) $$

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

$$ \sigma_{1,Poi} = \int_{0}^{1/\sqrt{3}} \frac{1-3x^2}{2} dx + \int_{1/\sqrt{3}}^{1} \frac{3x^2-1}{2} dx \quad (A12) $$

$$ = \left[ \frac{x(1-x^2)}{2} \right]_0^{1/\sqrt{3}} + \left[ \frac{x(x^2-1)}{2} \right]_{1/\sqrt{3}}^1, \quad (A13) $$

and finally:

$$ \sigma_{1,Poi} = 2/(3\sqrt{3}). \quad (A14) $$

## Appendix B: PDFs of position and phase

Figure 16 shows the pdfs of (a) transversal positions $y_i$ and (b) phases $\phi_i$ corresponding to the sample obtained by the Lagrangian method in the reference case ($\alpha = 2$, Re = 10 000) with a Péclet number $Pe = 10^5$. The positions follow a distribution based on the mean velocity profile and the phases conform to an almost uniform distribution.

<figure><img src="image_17.png" /><figcaption>FIG. 16. PDF of the (a) transversal positions $y_i$ and (b) phases $\phi_i$ corresponding to the sample studied in section IV A.</figcaption></figure>

## Appendix C: RTD of low energy wave

Figure 17 shows the RTDs for $L = \lambda$ and $L = 10\lambda$ in a case of very low wave energy, when there are no negative longitudinal velocities: the patterns of peaks and valleys are no longer present at large residence times.
<figure><img src="image_18.png" /><figcaption>FIG. 17. RTDs of the perturbed flow for (a) $L = \lambda$ and (b) $L = 10\lambda$ normalized by the mean residence time $\tau_m = L/U_m$. The physical parameters are here $Re = 12000$ and $\alpha = 1.6475$, the resulting nonlinear wavespeed is $c^{nl} = 0.331$.</figcaption></figure>

[1] E. Villermaux, "Mixing versus stirring," *Annual Review of Fluid Mechanics*, vol. 51, pp. 245-273, 2019.

[2] J. M. Ottino, "Mixing, chaotic advection, and turbulence," *Annual Review of Fluid Mechanics*, vol. 22, no. 1, pp. 207-254, 1990.

[3] H. Aref, J. R. Blake, M. Budišić, S. S. S. Cardoso, J. H. E. Cartwright, H. J. H. Clercx, K. El Omari, U. Feudel, R. Golestanian, E. Gouillart, G. F. van Heijst, T. S. Krasnopolskaya, Y. Le Guer, R. S. MacKay, V. V. Meleshko, G. Metcalfe, I. Mezić, A. P. S. de Moura, O. Piro, M. F. M. Speetjens, R. Sturman, J.-L. Thiffeault, and I. Tuval, "Frontiers of chaotic advection," *Rev. Mod. Phys.*, vol. 89, p. 025007, Jun 2017.

[4] P. V. Danckwerts, "Continuous flow systems. distribution of residence times," *Chem. Eng. Sc.*, vol. 2, pp. 1-13, 1953.

[5] D. E. Hull and J. W. Kent, "Radioactive tracers to mark interfaces and measure intermixing in pipelines," *Ind. & Eng. Chem.*, vol. 44, pp. 2745-2750, 1952.

[6] H. M. Haitjema, "On the residence time distribution in idealized groundwatersheds," *J. Hydrol.*, vol. 172, pp. 127-146, 1995.

[7] D. Bošković and S. Loebbecke, "Modelling of the residence time distribution in micromixers," *Chem. Eng. J.*, vol. 135S, pp. S138-S146, 2008.

[8] F. Trachsel, A. Günther, S. Khan, and K. F. Jensen, "Measurement of residence time distribution in microfluidic systems," *Chem. Eng. Sci.*, vol. 60, pp. 5729-5737, 2005.

[9] E. Gouillart, O. Dauchot, and J.-L. Thiffeault, "Measures of mixing quality in open flows with chaotic advection," *Physics of fluids*, vol. 23, no. 1, 2011.

[10] D. B. Spalding, "A note on mean residence-times in steady flows of arbitrary complexity," *Chem. Eng. Sc.*, vol. 9, pp. 74-77, 1958.

[11] P. V. Danckwerts, "Local residence-times in continuous-flow systems," *Chem. Eng. Sc.*, vol. 9, pp. 78-79, 1958.

[12] F. Raynal and P. Carrière, "The distribution of "time of flight" in three dimensional stationary chaotic advection," *Phys. Fluids*, vol. 27, p. 043601, 2015.

[13] B. Hervieu, A. Besq, and J. Angilella, "Residence times of polydisperse dilute suspensions in sheared and extensional flows," *Physics of Fluids*, vol. 34, no. 5, 2022.

[14] N. Poumaère, B. Pier, and F. Raynal, "Residence time distributions for in-line chaotic mixers," *Phys. Rev. E*, vol. 106, p. 015107, 2022.

[15] L. Oteski, Y. Duguet, and L. Pastur, "Lagrangian chaos in confined two-dimensional oscillatory convection," *Journal of Fluid Mechanics*, vol. 759, pp. 489-519, 2014.

[16] P. Meunier, "Geoinspired soft mixers," *J. Fluid Mech.*, vol. 903, p. A15, 2020.

[17] C. Lefranc, O. Detournay, and P. Meunier, "Gas injection into a tilted rotating cylinder," *Phys. Fluids*, vol. 35, p. 066602, 2023.

[18] J. Qu, D. Henry, S. Miralles, V. Botton, and F. Raynal, "Chaotic mixing in an acoustically driven cavity flow," *Physical Review Fluids*, vol. 7, no. 6, p. 064501, 2022.

[19] J. Lou, J. Johnston, T. Y. Cath, D. Martinand, and N. Tilton, "Computational fluid dynamics simulations of unsteady mixing in spacer-filled direct contact membrane distillation channels," *J. Membr. sci.*, vol. 622, p. 118931,
2021.

[20] C. Venditti, V. Biagioni, A. Adrover, and S. Cerbelli, "Impact of transversal vortices on the performance of open-tubular liquid chromatography," *Journal of Chromatography A*, vol. 1685, p. 463623, 2022.

[21] K. Kataoka, H. Doi, T. Hongo, and M. Futagawa, "Ideal plug-flow properties of taylor vortex flow," *J. Chem. Eng. Japan*, vol. 8, pp. 472-476, 1975.

[22] K. Kataoka, H. Doi, and T. Komai, "Heat/mass transfer in taylor vortex flow with constant axial flow rates," *Int. J. Heat Mass Transf.*, vol. 20, pp. 57-63, 1977.

[23] C. C. Lin, "On the stability of two-dimensional parallel flows, Part III.—Stability in a viscous fluid," *Q. Appl. Math.*, vol. 3, pp. 277-301, 1946.

[24] L. H. Thomas, "The stability of plane Poiseuille flow," *Phys. Rev.*, vol. 91, pp. 780-783, 1953.

[25] C. E. Grosch and H. Salwen, "The stability of steady and time-dependent plane Poiseuille flow," *J. Fluid Mech.*, vol. 34, pp. 169-187, 1968.

[26] S. A. Orszag, "Accurate solution of the Orr-Sommerfeld stability equation," *J. Fluid Mech.*, vol. 50, pp. 689-703, 1971.

[27] S. J. Davies and C. M. White, "An experimental study of the flow of water in pipes of rectangular section," *Proc. R. Soc. Lond. A*, vol. 119, pp. 92-107, 1928.

[28] J. T. Stuart, "Nonlinear stability theory," *Annu. Rev. Fluid Mech.*, vol. 3, pp. 347-370, 1971.

[29] D. Meksyn and J. T. Stuart, "Stability of viscous motion between parallel planes for finite disturbances," *Proc. R. Soc. Lond. A*, vol. 208, pp. 517-526, 1951.

[30] J. T. Stuart, "On the non-linear mechanics of hydrodynamic stability," *J. Fluid Mech.*, vol. 4, pp. 1-21, 1958.

[31] J. T. Stuart, "On the non-linear mechanics of wave disturbances in stable and unstable parallel flows part 1. the basic behaviour in plane Poiseuille flow," *J. Fluid Mech.*, vol. 9, pp. 353-370, 1960.

[32] J. 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," *J. Fluid Mech.*, vol. 9, pp. 371-389, 1960.

[33] J.-P. Zahn, J. Toomre, E. A. Spiegel, and D. O. Gough, "Nonlinear cellular motions in Poiseuille channel flow," *J. Fluid Mech.*, vol. 64, pp. 319-346, 1974.

[34] J. Jiménez, "Bifurcations and bursting in two-dimensional Poiseuille flow," *Phys. Fluids*, vol. 30, p. 3644, 1987.

[35] U. Ehrenstein and W. Koch, "Three-dimensional wavelike equilibrium states in plane Poiseuille flow," *J. Fluid Mech.*, vol. 228, pp. 111-148, 1991.

[36] R. Artuso, L. Cavallasca, and G. Cristadoro, "Dynamical and transport properties in a family of intermittent area-preserving maps," *Physical Review E*, vol. 77, no. 4, p. 046206, 2008.

[37] P. G. Saffman, "A theory of dispersion in a porous medium," *J. Fluid Mech.*, vol. 6, pp. 321-349, 1959.

[38] M. Dentz, M. Icardi, and J. J. Hidalgo, "Mechanisms of dispersion in a porous medium," *J. Fluid Mech.*, vol. 841, pp. 851-882, 2018.

[39] A. Puyguiraud, P. Gouze, and M. Dentz, "Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media," *Phys. Rev. Lett.*, vol. 126, 2021.

[40] S. M. Cox, P. G. Drazin, S. C. Ryrie, and K. Slater, "Chaotic advection of irrotational flows and of waves in fluids," *J. Fluid Mech.*, vol. 214, pp. 517-534, 1990.

[41] P. J. Schmid and D. S. Henningson, *Stability and Transition in Shear Flows*. Springer, New York, 2001.

[42] B. Pier and P. Schmid, "Linear and nonlinear dynamics of pulsatile channel flow," *J. Fluid Mech.*, vol. 815, pp. 435-480, 2017.

[43] K. Goda, "A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows," *J. Comp. Phys.*, vol. 30, pp. 76-95, 1979.

[44] I. Raspo, S. Hugues, E. Serre, A. Randriamampianina, and P. Bontoux, "A spectral projection method for the simulation of complex three-dimensional rotating flows," *Computers and Fluids*, vol. 31, pp. 745-767, 2002.

[45] B. Pier, *Dynamique des écoulements ouverts : instabilités et transition, courbure et rotation*. Habilitation à diriger des recherches. École centrale de Lyon & Université de Lyon, 2015.

[46] B. Pier, "PackstaB," Mar. 2023.

[47] J. B. Weiss, "Transport and mixing in traveling waves," *Phys. Fluids A*, vol. 3, pp. 1379-1384, 1991.

[48] H. Aref and S. W. Jones, "Enhanced separation of diffusing particles by chaotic advection," *Physics of Fluids A: Fluid Dynamics*, vol. 1, no. 3, pp. 470-474, 1989.

[49] D. V. Khakhar, J. G. Franjione, and J. M. Ottino, "A case study of chaotic mixing in deterministic flows: The partitioned-pipe mixer," *Chem. Eng. Sci.*, vol. 42, pp. 2909-2926, 1987.

[50] I. Mezić, S. Wiggins, and D. Betz, "Residence-time distributions for chaotic flows in pipes," *Chaos*, vol. 9, pp. 173-182, 1999.