P G Correa
email: pcorrea@exa.unicen.edu.ar

J M Gomba
email: jgomba@exa.unicen.edu.ar

J R Mac Intyre

S Ubal

C A Perazzo Imettyb

J P Hulin

H Auradou

Influence of aspect ratio on vortex formation in X-junctions: direct numerical simulations and eigenmode decomposition

We study numerically the appearance and number of axial vortices in the outlets of X-shaped junctions of two perpendicular channels of rectangular sections with facing inlets. We explore the effect of the aspect ratio of the cross section, AR, on the number of vortices created at the center of the junction. Direct numerical simulations (DNS) performed for different values of the Reynolds number Re and AR demonstrate that vortices with their axis parallel to the outlets, referred to as axial vortices, appear above critical Reynolds numbers Re c . As AR increases from 1 to 11, the number of vortices observed increases from 1 to 4, independently of Re. For AR = 1, the single axial vortex induces an interpenetration of the inlet fluids in the whole section; instead, for larger AR's for which more vortices appear, the two inlet fluids remain largely segregated in bands, except close to the vortices. The linear stability analysis demonstrates that only one leading eigenmode is unstable for a given set of values of AR and Re. This mode provides a simplified model of the flow field, reproducing its key features such as the number of vortices and their distance. Its determination with this method requires a much smaller computational load than the DNS. This approach is shown to allow one to determine quickly and precisely the critical Reynolds number Re c and the sensitivity function S which characterizes the influence of variations of the base flow on the unstable one.

I. INTRODUCTION

One of the challenges of the lab-on-a-chip technology is to develop a micro reactor and analytical equipment that operates at the scale of a few micrometers and can be used to mix [START_REF] Capretto | Micromixing within microfluidic devices[END_REF][START_REF] Lee | Passive mixers in microfluidic systems: A review[END_REF] , separate [START_REF] Gossett | Label-free cell separation and sorting in microfluidic systems[END_REF][START_REF] Jackson | Advances in microfluidic cell separation and manipulation[END_REF][START_REF] Sajeesh | Particle separation and sorting in microfluidic devices: a review[END_REF] , trap [START_REF] Nilsson | Review of cell and particle trapping in microfluidic systems[END_REF] and transport chemicals or small particles [START_REF] Amini | Intrinsic particle-induced lateral transport in microchannels[END_REF] . To allow all those operations, a perfect control of the fluids and of their flow is required. Knowing and understanding how vortices are created is, for instance, crucial to combine fluids and to ensure a good mixing between them in order to enhance mass or heat transfer [START_REF] Abed | Heat transfer enhancement a cross-slot micro-geometry[END_REF] . But, due to the small dimensions of those systems and because of the viscosity of the fluids, vortices cannot be obtained by making the fluid turbulent. Yet, even in relatively viscous flows and in geometries of small size, vortices can be created by carving cavities [START_REF] Howell | A microfluidic mixer with grooves placed on the top and bottom of the channel[END_REF] or by deviating and bending the flow, thereby playing with the channels geometry [START_REF] Sudarsan | Multivortex micromixing[END_REF][START_REF] Sudarsan | Fluid mixing in planar spiral microchannels[END_REF] . This can be achieved by injecting fluids in X-and T-junctions in opposing channels: this provides hope for the development of passive mixers [START_REF] Kumar | Single-phase fluid flow and mixing in microchannels[END_REF] , using only static parts to promote mixing without the need of any external energy supply or movable elements [START_REF] Sudarsan | Multivortex micromixing[END_REF] . Another important practical issue is the active control of the flow within such junctions and, particularly, of the development of instabilities [START_REF] Aryshev | Stability of colliding flows[END_REF] .

Originally, X-junctions were developed because their geometry allows the generation of a stagnation point away from channel walls. In the vicinity of this particular point, the flow is purely planar, elongational, and is free of shear and vorticity along the symmetry axes [START_REF] Haward | Elastic instabilities in planar elongational flow of monodisperse polymer solutions[END_REF] . These properties have been put at work in many areas of research, including studies of polymer macromolecules dynamics, such as DNA [START_REF] Perkins | Single polymer dynamics in an elongational flow[END_REF][START_REF] Xu | Exploring both sequence detection and restriction endonuclease cleavage kinetics by recognition site via single-molecule microfluidic trapping[END_REF] or fluid rheology [START_REF] Rocha | On extensibility effects in the cross-slot flow bifurcation[END_REF][START_REF] Haward | Stagnation point flow of wormlike micellar solutions in a microfluidic cross-slot device: Effects of surfactant concentration and ionic environment[END_REF][START_REF] Arratia | Elastic instabilities of polymer solutions in cross-channel flow[END_REF][START_REF] Haward | Microfluidic extensional rheometry using stagnation point flow[END_REF][START_REF] Burshtein | Inertioelastic Flow Instability at a Stagnation Point[END_REF] and for imposing controlled deformations to cells, vesicles or droplets [START_REF] Kantsler | Critical dynamics of vesicle stretching transition in elongational flow[END_REF] .

The flow field is quite simple at a low velocity, but increasing it leads to a more complex vorticity field in the intersection and in the outlets. For a Newtonian fluid of constant viscosity ν, the onset of the different regimes is determined by the Reynolds number Re = U W /ν, where U is a characteristic velocity and W is the dimensional width of the channels. For an X-junction with channels crossing at an angle α = 90 o , the flow along the outlet channels at low Re's is symmetric with respect to the plane y = 0 (see Figs. 1(b) and 1(c)): the two injected fluids remain segregated. At higher Re's, an axial vortex appears at the intersection of the channels and extends towards each outlet. In channels with a square cross-section [START_REF] Abed | Heat transfer enhancement a cross-slot micro-geometry[END_REF][START_REF] Haward | Microfluidic extensional rheometry using stagnation point flow[END_REF][START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF] , the transition between the two regimes occurs for Re c ∼ 40 , while for circular channels [START_REF] Correa | Three-dimensional flow structures in X-shaped junctions: Effect of the Reynolds number and crossing angle[END_REF] it is slightly higher with Re c ∼ 48.

In these square cross-section junctions, the redirection of the flow around the corners at the junction of the streams leads to the formation of small Dean vortices at the intersection of the channels [START_REF] Dean | Fluid Motion in a Curved Channel[END_REF][START_REF] Zhang | Experimental investigation of three-dimensional flow regimes in a cross-shaped reactor[END_REF] , even for Re < Re c . This secondary flow consists of a double pair of counter rotating vortices positioned symmetrically on the four corners of the outlet channels, but out of the z-axis (see Fig. 1 in Ref. 27). Both fluids remain segregated by the plane y = 0 despite the presence of these small structures. The intensity of these vortices increases with the flow rate and, at the critical Reynolds number Re c , two opposite vortices, out of the four, intensify and the symmetry is broken. These latter vortices finally merge into a single, steady, stream-wise, vortex centered on the z-axis. This vortex was shown recently to be strong enough to deform the membrane of living cells that in turn become porous to nano materials [START_REF] Kang | Intracellular nanomaterial delivery via spiral hydroporation[END_REF] which opens the possibility to use such flow for hydroporation. The creation of "hot spots" by the vortices where chemical reactions are enhanced were also found to be a possible method to study chemical reactions [START_REF] Damian | From flow focusing to vortex formation in crossing microchannels[END_REF][START_REF] Borgogna | Inertia-driven enhancement of mixing efficiency in microfluidic cross-junctions: a combined Eulerian/Lagrangian approach[END_REF][START_REF] Moloudi | Inertial particle focusing dynamics in a trapezoidal straight microchannel: application to particle filtration[END_REF][START_REF] Lee | Three-Dimensional Vortex-Induced Reaction Hot Spots at Flow Intersections[END_REF] . If the velocity is further increased, the flow becomes unsteady at Re ∼ 100 and periodic oscillations are observed [START_REF] Zhang | Experimental investigation of three-dimensional flow regimes in a cross-shaped reactor[END_REF] .

The dynamics of the destabilization and merging of the Dean vortices was studied experimentally using time-resolved flow velocimetry by Burshtein et al [START_REF] Burshtein | Controlled symmetry breaking and vortex dynamics in intersecting flows[END_REF] . They also investigated the influence of the aspect ratio AR = H/ W , where W and H are respectively the width and height of the channels, on the dynamics of the formation of the central vortex. They confirmed that the geometry influences the nature of the transition. For wide channels corresponding to AR < 0.5, the flow recovers its symmetry at similar Re c values for both decreasing or increasing inlet velocity variation ramps. For AR ≥ 0.5, the transition changes and becomes subcritical. As a result, when the flow velocity decreases, the flow configuration reverts from a single vortex centered in the outlet channel to a Dean-like one with symmetric vortices at a Reynolds number Re * c , lower than Re c . This confirmed previous results obtained numerically by Haward et al 23 : these authors had shown that the symmetry-breaking flow bifurcation which is supercritical (non symmetric) for wide channels becomes subcritical (pitchfork bifurcation) for deep ones. These former studies demonstrate therefore that the geometry of the channel section influences the critical Reynolds number for the transition between the different regimes. Similar observations were made in Y or T junctions [START_REF] Chan | Microscopic investigation of vortex breakdown in a dividing T-junction flow[END_REF] and in X-junctions with varying crossing angles [START_REF] Correa | Three-dimensional flow structures in X-shaped junctions: Effect of the Reynolds number and crossing angle[END_REF] .

All these studies dealt with the dynamics of vortex formation and of the steady "engulfment" regime at relatively low values of AR (0.4 < AR < 3.8). Moreover, in this range, there is only one axial vortex inside the flow section. Higher aspect ratios (AR 4) are however also of interest, for instance for rheology, because high AR's ensure that the strain rate is approximately uniform throughout the height which makes rheological measurements easier. To our knowledge, larger aspect ratios were only used in the experiments of Kalashnikov and Tsiklauri [START_REF] Kalashnikov | Ordered three-dimensional structures resulting from instability of two-dimensional flow in crossed channels[END_REF] . For very large values AR = 262 and AR = 32, these authors observed a periodic array of vortices stacked over the full height of the channels. When increasing the flow, this structure appeared above respective critical Reynolds numbers Re c = 55 and Re c ∼ 43. The transition is, again, subcritical so that, when the flow is lowered, the vortices disappear for Re c = 38.5 and Re c ∼ 30, respectively. Qualitative visualizations suggested that such devices may mix large quantities of fluids at a relatively low Reynolds number.

Recently, the so-called "structural sensitivity" formalism has been developed to predict the effect of perturbations brought to a flow in localized regions [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF] : it allows detecting the locations where an external actuation either triggers an instability, or delays it. Giannetti & Luchini [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF] and Marquet et al. [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF] applied this approach to a two dimensional analysis of the wake generation downstream of a stationary cylinder. Although the magnitude of linear instability modes increases continuously with the distance, the region where a feedback force produces the largest change of the leading eigenvalues (i.e. where the structural sensitivity is largest), is close to the obstacle and is the best location to control the wake generation. Using the same approach, Lashgari et al. [START_REF] Lashgari | The planar X-junction flow: Stability analysis and control[END_REF] studied the flow instabilities in an X-junction (three inlets, one outlet), for which the base flow is two dimensional and the perturbations are decomposed into biglobal modes [START_REF] Theofilis | Global Linear Instability[END_REF] . In this case, the sensitivity is highest at the edge of recirculation bubbles at the corners of the junction. Chen et al. [START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF] performed a global sensitivity analysis to study flow through a T-shaped pipe bifurcation. They observed recirculation zones resembling the traditional bubble-type breakdown. These regions are highly sensitive to localized feedback forces and, in contrast with observations on three inlet X-junctions [START_REF] Lashgari | The planar X-junction flow: Stability analysis and control[END_REF] , the flow separation at the corners of the T does not display a clear-cut sensitivity.

The present paper is devoted to the numerical study of the transitions between the different regimes discussed above for flow in X-junctions of channels of rectangular cross-sections. More specifically, we report direct numerical simulation (DNS) for AR values between 1 and 11 corresponding to the appearance of 1 to 4 vortices at the junction intersection and along the height (z-axis); note that these computations become more demanding as AR increases. For each value of AR, the influence of the Reynolds number on the structure of the flow is studied at Reynolds numbers below 100 corresponding to stationary flow regimes and the possible subcritical nature of the instability leading to the appearance of the vortices is investigated. We analyze in particular whether the number of vortices in the height only depends on AR or whether Re also has an influence. A three dimensional global stability analysis has then been performed to model the instability leading to the formation of the vortices: we shall compare its predictions for the critical values of AR and Re, and the number and spacing of the vortices to those of the DNS. The sensitivity of the flow to the application of feedback forces is finally considered, which may be of interest for flow control applications.

II. PROBLEM STATEMENT AND FORMULATION

We consider the flow of two facing incompressible Newtonian liquids within an X-junction with an angle 90 o between the inlet and outlet channels. The device consists of two perpendicular channels of length 2L and identical rectangular cross-sections, with the inlets aligned with the y-axis and the outlets parallel to the x-axis, as shown in Fig. 1(a). The dimensional width and height of the channels are W and H, respectively, AR = H/ W being the aspect ratio. The origin O of the coordinate system is located at the center of symmetry of the junction. In the following, all lengths are normalized by the width W , and the velocities by the mean velocity at the entrances, U in . With these scales, the normalized height and width of the channels are named respectively H(= AR) and W (= 1), respectively.

The flow within the junction is governed by the time-dependent three-dimensional incompressible Navier-Stokes equations:

∂u ∂t + (u • ∇)u = -∇p + 1 Re ∇ 2 u, (1a)
∇ • u = 0, (1b)
where u ≡ (u, v, w) and p are the normalized velocity and pressure, respectively. The Reynolds number is Re = U in W /ν, where ν is the viscosity of the liquid. We impose the same parabolic inflow conditions at both inlets (Poiseuille solution for rectangular pipes), a stress-free outflow for the outlets, and u = 0 at the walls [START_REF] Correa | Three-dimensional flow structures in X-shaped junctions: Effect of the Reynolds number and crossing angle[END_REF] .

A. Linear stability analysis

To perform a standard global linear stability analysis [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Fani | Investigation of the steady engulfment regime in a three-dimensional T-mixer[END_REF] , the variables in Eq.( 1) are written as the sum of a steady base flow {u 0 , p 0 } and an unsteady small perturbation {u ′ , p ′ }. The base state shares the same initial and boundary conditions as {u, p}. In our case, it is a symmetrical solution of Eq.( 1) in which the liquid coming from each inlet splits equally between both outlets, with the streamlines completely segregated by the plane y = 0, as shown in Figures 1(b) and 1(c) [START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF][START_REF] Correa | Three-dimensional flow structures in X-shaped junctions: Effect of the Reynolds number and crossing angle[END_REF][START_REF] Cachile | Stokes flow paths separation and recirculation cells in X-junctions of varying angle[END_REF] .

The perturbations are decomposed into global modes, i.e. {u ′ , p ′ }={ û, p}(x, y, z) exp (λ t), where λ = σ + i ω is a complex eigenvalue. The real part σ is the growth rate and the imaginary part ω is the oscillation angular frequency of the perturbation. If the growth rate is positive for, at least, one eigenvalue, the flow is linearly unstable, otherwise the perturbation decays to zero. This linearization of the flow around the base flow and the subsequent eigenmode decomposition, result in the following direct eigenvalue equation for the perturbations [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF] λ û + (u

0 • ∇) û + ( û • ∇) u 0 = 1 Re ∇ 2 û -∇p , (2a)
∇ • û = 0. ( 2b
)
Since {u, p} and {u 0 , p 0 } share the same boundary conditions, the perturbations satisfy homogeneous conditions at all boundaries. The components of each velocity field are u 0 := (u 0 , v 0 , w 0 ) and û := (û, v, ŵ).

B. Adjoint problem and structural sensitivity

The adjoint of a linear operator is a useful concept in functional analysis that has been widely applied to problems in turbulence control, receptivity, and transition, and it has recently been used for the analysis of flow within micro-junctions [START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF][START_REF] Fani | Investigation of the steady engulfment regime in a three-dimensional T-mixer[END_REF] . Following Chomaz [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF] , and in order to evaluate the sensitivity of the solutions of eq. ( 1), we apply some well known concepts related to the adjoint problem. In the next lines, we summarize the main equations that we need to solve in order to compute the sensitivity function. A more detailed derivation of this theory can be found in refs. [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF][START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF][START_REF] Hill | A theoretical approach for analyzing the restabilization of wakes[END_REF] . Hill [START_REF] Hill | A theoretical approach for analyzing the restabilization of wakes[END_REF] and Giannetti & Luchini [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF] developed the theory of structural analysis, and showed that the adjoints fields {u + , p + } = { û+ , p+ }(x, y, z) exp (-λt) associated to the global mode { û, p, λ} satisfy the eigenvalue problem:

-λ * û+ -(∇u 0 ) • û+ + u 0 • ∇ û+ = ∇p + - 1 Re ∇ 2 û+ , (3a)
∇ • û+ = 0 . ( 3b
)
where we use the notation ((∇a) • b) i := j b j ∂a j /∂x i 36 . The boundary conditions for the adjoint modes are p+ n = (u 0 • n) û+ + Re -1 (n • ∇ û+ ) at the outlets [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Fani | Investigation of the steady engulfment regime in a three-dimensional T-mixer[END_REF] , and û+ = 0 for the rest of the boundaries.

Solving equations ( 2) and (3) allows one to detect the core region of the instability (for some geometries the region is named wavemaker [START_REF] Lashgari | The planar X-junction flow: Stability analysis and control[END_REF][START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF][START_REF] Juniper | The structural sensitivity of open shear flows calculated with a local stability analysis[END_REF][START_REF] Qadri | Structural sensitivity of spiral vortex breakdown[END_REF] ), i.e. the spots where a local feedback force results in the largest drift of the most "dangerous" eigenvalue [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF] . This force could be, for example, the result of the action of an actuator that reacts to the local velocity of the flow at the point where the actuator is operating. Giannetti & Luchini find that the maximum change in the dominant eigenvalue is induced at the location of the greatest overlap of the direct and adjoint modes, and define the structural sensitivity function as [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF][START_REF] Qadri | Structural sensitivity of spiral vortex breakdown[END_REF] :

S(x, y, z) = û û+ û, û+ , (4)
where a, b ≡ V a * • b dV , the asterisk denoting the conjugate of a complex quantity [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Fani | Investigation of the steady engulfment regime in a three-dimensional T-mixer[END_REF] .

The regions of the flow where S is large are the sites where a feedback force will produce the strongest disturbance on the flow, and it identifies the region where the triggering of the instability occurs [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF] . Notice that the location of large values of the direct eigenfunction does not necessarily play a special role in determining the spectrum of a stability equation unless the adjoint eigenfunction is also large at the same spot. This fact is also helpful for the numerical simulations, because it gives a criterion for detecting the region where the mesh should be denser in order to capture accurately the global modes.

C. Numerical implementation and validation

Equations ( 1)-( 3) were discretized and solved numerically by a finite element method. Polynomial shape functions P2 and P1 were used for the discretization of the velocity and pressure, respectively. The time dependent Navier-Stokes equations (1) were solved by means of a backward differentiation scheme with adaptive time stepping. The steady-state base flow configuration {u 0 , p 0 }(x, y, z), was solved through an iterative method, the Generalized Minimal RESiduals (GMRES), preconditioned using a standard multigrid algorithm.

A convergence criterion of 10 -3 is used for the relative error defined by a weighted euclidean norm for two successive iteration steps (see Correa et al. [START_REF] Correa | Three-dimensional flow structures in X-shaped junctions: Effect of the Reynolds number and crossing angle[END_REF] ). In the stability problem given in equations ( 2), the eigenvalues were computed employing a variant of the implicitly restarted Arnoldi method in the ARPACK routine [START_REF] Lehoucq | ARPACK Users' Guide[END_REF] . The junction domain was meshed with tetrahedral elements for the time dependent problem and hexahedral elements for the base flow and eigenvalue problems, with a higher concentration near the walls and in the crossing region of the X-junction. We carried out convergence studies to estimate the minimal number of grid elements needed to obtain accurate results.

The accuracy of our numerical procedure is established by a grid convergence study over the leading eigenvalue (see Table I). Direct and adjoint eigenvalues, Eqs. ( 2) and (3), were obtained for four meshes with AR = 1 and L = 7. The differences among the eigenvalues is less that 1.3% for all the cases considered, so we used M3 in the rest of the study in order to achieve an accurate spatial description of the corresponding eigenmodes. The minimum size of the domain was determined by analyzing the influence of the lengths of the channels on the value of λ 38 . According to between reliable results and computational cost to capture fully the dynamics of the flow and its instability.

We also validated the Direct Numerical Simulations (DNS) codes and the stability analysis method by considering the flow in a three- They reported a segregated flow regime for Re < 175 and a first flow bifurcation occurring at Re = 175. For larger Re's, the symmetry of the flow is broken in an engulfment regime that mixes both incoming streams. We reproduced their results by means of the finite element method. The accuracy of the computation was controlled by a convergence study on the mesh parameters. A reliable result for the leading eigenvalues is obtained in a mesh domain of 23460 elements, finding σ = -1.542 × 10 -2 . This value differs by only 1.4% from the eigenvalue reported in ref. [START_REF] Fani | Investigation of the steady engulfment regime in a three-dimensional T-mixer[END_REF] . Moreover, the shapes of the direct and adjoint modes associated with the leading eigenvalue match very well.

III. RESULTS

A. DNS investigation of the flow structure

We analyzed first the flow at low Re values, for which the flow remains stable. For instance, for Re = 22 (see Fig. 1(b) and 1(c)), the inlet flows are split into two equal parts after they meet in the crossing region and are in contact inside the outlet channels only in the plane y = 0. This behavior is observed at small Reynolds numbers, independent of the value of the aspect ratio AR > 1. Figure 1(c) shows that this flow displays a noticeable extensional structure in the z direction. For Re > Re c , the segregated flow configuration becomes unstable and its symmetry is broken by the appearance of vortices with their axis parallel to the x direction. Their formation is triggered by inertial effects which let the incoming liquids cross the plane y = 0 (see below), resulting in a swirling motion toward the outlets. The number of layers is directly related to the number n of vortices. As shown above, in the plane y = 0, each vortex is located between stripes of two different fluids. Since there are fluid layers between both the upper and lower walls and the nearest vortices, the total number of layers is n + 1. We note that the layered distribution of the two fluids is the same in both exits, with the stripes of a given fluid at the same location. After a transient flow settling phase, the flow field and the fluid distribution become time-independent in the range of Reynolds numbers investigated. Quantitatively, we observed one single vortex for AR 3.8, in agreement with Refs. [START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF][START_REF] Zhang | Trapping region of impinging jets in a cross-shaped channel[END_REF] , two vortices for 3.8 AR 8.5 and three for 8.5

AR < 11. For AR = 11, which represents the highest demand for computationalpower which we can deal with, four vortices appear.

Figures 3(a)-3(c) display the streamlines associated with the vortices in the middle plane x = 0 of the junction, respectively for 1, 2 and 3 vortices. For AR = 1, Fig. 3(a), the upward deflection of the blue fluid (on the figure) combined with the downward deflection of the orange one serves as a source for a counter-clockwise rotation which propagates thereafter in the z direction. For AR = 6, Fig. 3(b), the orange flow is deflected toward the upper and lower walls of the channel while the blue one is focused toward the middle. One has this time two sources of local rotation in opposite directions which also propagate along z, creating two vortices of opposite circulations parallel to the direction of the outlet. The vortices result therefore from an inertia driven distortion of the initial separation plane y = 0 of the two opposite flows encountering in the junction. If AR increases further, more distortions appear, leading to more vortices as shown in Fig. 3(c).

In Fig. 3(a), one notices that, since the two fluids are identical, a configuration in which the orange fluid is deflected upward and the blue fluid downward is equally possible, leading to a clockwise rotating vortex. In the same way, in Fig. 3(b), the orange fluid might as well advance farther in the center and the blue fluid near the walls. This, too, would result in a reversal of the two vortices with an orange fluid layer in the center of the section and blue layers near the ends. Similarly, for other values of AR, one may swap the locations of the two fluids in the outlet channels provided the rotation of the vortices is also reversed. Experimentally, both configurations should be observed with equal probabilities.

The effect of the transit of the liquid through the output channels is demonstrated by comparing Figs. 3(a,b,c) (x = 0) and 3(d,e,f) (x = 7). For all AR's, the vortical flow is limited to the region close to the center of the section at x = 0. Nevertheless, the aspect ratio affects the distribution of the two fluids at the outlets. For AR = 1, Fig. 3(d), the blueorange streamlines define an helical-shaped interface that is longer than that observed in the layered fluid distribution for AR = 6 and AR = 10 shown in Figs. 3(e) and 3(f), respectively. The vorticity distribution at the center differs therefore strongly from that at the outlet. We identify the 3D structure of the vortices by Q-contours [START_REF] Zhang | Experimental investigation of three-dimensional flow regimes in a cross-shaped reactor[END_REF][START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF][START_REF] Jeong | On the identification of a vortex[END_REF][START_REF] Chakraborty | On the relationships between local vortex identification schemes[END_REF] . The criterion Q is defined as

Q = 1 2 ( Ω 2 F -S 2 F ),
where Ω = 1 2 ∇u -(∇u) ⊤ and S = 1 2 ∇u + (∇u) ⊤ are the rotation and the strain-rate tensors, respectively, and ( . F ) stands for the Frobenius norm. In Figures 3(g)-3(i), Q isosurfaces are represented in red for AR = 1, 6 and 10; in these three cases, the vortices develop only in the region where the channels meet and hardly penetrate into the outlet channels. As a result, close to the exits, the axial velocity disappears and the parabolic profile and its vorticity distribution are recovered. The physical mechanism of the development of the vortices is the same for all AR analyzed; it is illustrated in Fig. 4 for AR = 6 by a time sequence of maps of the axial (x-component) vorticity and streamlines. Initially, in Fig. 4(a), the vorticity is non zero only close to the upper and lower walls and the flow lines are straight and horizontal. In Fig. 4(b), after Re has increased with time, the vorticity has diffused from the top and bottom walls towards z = 0, especially in the region close to the interface. In Fig. 4(c), four Dean cells start to develop on both walls of the intersection with an inversion of the sign of the vorticity close to these walls. In Figs. 4(d) and 4(e), the largest vorticity keeps developing close to the walls and the flow symmetry is lost. Finally the steady state corresponding to Fig. 3(b) above is reached in Fig. 4(f). The key influence on the triggering of the instability of the formation of the Dean vortices and the diffusion of the x-component of the vorticity was first explained for AR = 1 by Haward et al. [START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF] and Burshtein et al [START_REF] Burshtein | Inertioelastic Flow Instability at a Stagnation Point[END_REF] . The Multimedia view attached to Fig. 4 displays this same transition from zero to n vortices respectively for the streamlines and the vorticity for AR = 1, 6, and 10. We will show in Section III B that the regions with the maximum value for the sensitivity S develop four local x-vorticity maxima (close to the upper and lower walls) for z ± 1.6 (in the present case). Then, the actuation on these regions may hamper or boost vorticity generation. A first important result of the simulations is that the number n of vortices observed depends only on the aspect ratio AR and not on Re provided Re > Re c (Re c depends on AR). The variation of n with AR is shown in Fig. 5; one observes a stepwise increase with no overlap between the different steps. In the same way, we performed these simulations at different Re's (≥ Re c ) for AR = 6 (resp. 10) corresponding to n = 2 (resp. 3) vortices. We also included one case with AR = 11, for which one observes 4 vortices; further increments of AR are beyond our current computing capabilities. Figure 7(a) displays in the two cases the distances between the different vortices as a function of Re: their relative variation is less than 5%. The onset of the instability is analyzed in Fig. 8. To characterize it, we selected v max,z as the order parameter and explored its variation with Re for three different AR's [START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF] . These curves were obtained by means of numerical simulations using increasing and decreasing ramps of Re with steps ∆Re = 1 in the vicinity of the critical value Re c . The segregated stable flow is characterized by v max,z = 0, which means that the interface between the two fluids remains in the plane y = 0. The sudden increase of v max,z indicates the onset of the instability. As can be seen from the three examples in Fig. 8, there is an hysteresis for all values of the number n of vortices. The critical Re's for the increasing ramps are Re c = 44, 37 and 41, and for the decreasing ramps Re * c = 41, 30 and 36, for AR = 1, 6 and 10, respectively. The fact that Re * c < Re c for all the explored values of AR implies therefore that the instability leading to the appearance of vortices is subcritical in all these cases.

B. Linear stability analysis

In this section we investigate the global stability of the flow and compare the results, such as the values of Re c , the number of vortices and their spacing, to those reported in section III A) for the same aspect ratios AR. The global stability of the flow is analyzed by the computation of the eigenvalues and eigenmodes of the direct problem (Eqs. 2). The natural choice of the base flow required in the calculation is the steady-state symmetric segregated flow field. In order to obtain the latter for Re higher than Re c , we compute the solution in half of the domain and impose a symmetry boundary condition in the plane y = 0. The solution is then mirrored to the other half of the domain.

A discrete set of eigenvalues is shown in Fig. 9 for stable (a) and unstable (b) flows with AR = 1. The excellent agreement between the direct and adjoint spectra demonstrates the reliability of the numerical procedure. The eigenvalues λ = σ + iω are distributed symmetrically with respect to the real axis. In all the simulations, the eigenvalue with the largest growth rate σ corresponds to a steady-state perturbation (ω = 0). This is also the case for the two other eigenvalues with σ closest to 0, labeled as 2 and 3 in Fig. 9. The change of the growth rate σ of the leading eigenvalue from negative to positive shown in Fig. 10 allows us to determine precisely the critical Reynolds number, which is found to be Re c = 43.5 ± 0.15. The three-dimensional global mode associated to the leading eigenvalue (labeled 1 in Fig. 9) is shown in Fig. 11 for Re = 50. The eigenmodes are spatially located mainly in the outlet pipes. The v and ŵ components are composed of two main lobes, both symmetrical with respect to the plane x = 0 and with a maximum of intensity in the intersection (see arrows in Fig. 12). On the other hand, û consists of four lobes, anti-symmetrical with respect to the plane x = 0, with their highest and lowest intensities outside the junction, at a distance x = ±1.7 (see Fig. 12). Yellow corresponds to the value -10 -6 and green to 10 -6 . The velocity field of the leading eigenmode is displayed in Fig. 12 in four sections of an outlet branch. Similar to the corresponding DNS, one observes in these (y, z) cut planes a single vortical motion around the junction center. Although all perturbation components are of the same order of magnitude, v is approximately 2.5 times higher than û and ŵ. The maximum of v occurs at (x, y, z) = (0, 0, ±0.25) where (û, ŵ) = (0, 0). The component ŵ reaches its maximum value in the plane z = 0 near the lateral walls, at the beginning of the outlet branches (x = ±0.5).

The stability has also been studied for the junctions with AR = 6 and 10 and the results compared to those of the DNS simulations. Like for AR = 1, the leading eigenvalues are real and their variation with Re is shown in Fig. 13. The critical Reynolds numbers are respectively Re c = 39 ± 0.15 and 42.4 ± 0.1 for AR = 6 and 10.

Figure 14 displays isosurfaces of the velocity perturbation components of the mode associated with the leading eigenvalue for AR and Re values corresponding to 2, Figs. 14(a,b,c), and 3, Figs. 14(d,e,f), vortices. In these cases, like for AR = 1, the values of v and ŵ are largest inside the intersection and the maximum for û lies in the outlet branches for û. Also, û (1.5 × 10 -6 ).

remains antisymmetric with respect to the plane x = 0 while v and ŵ are still symmetric.

There are however significant changes compared to AR = 1. For û, the group of four lobes in each outlet channel is split, for AR = 6, into two distant groups with a new pair of lobes at half height while, for AR = 10, an additional pair appears. For v, there are three lobes for AR = 6, instead of two for AR = 1, and they are less elongated along x; for AR = 10, four elongated lobes are obtained. For ŵ, there are, for AR = 6, four pairs of medium size, mildly elongated lobes and two small instead of two large, very elongated ones; for AR = 10, there are two more lobes and they are more elongated. Regarding the symmetry of the lobes with respect to the plane z = 0, it is determined by the even (n = 2 for AR = 6) or odd (n = 1 (resp. 3) for AR = 1 (resp. 10)) values of the number of vortices. In Fig. 15, the three velocity perturbation components of the corresponding leading eigenmode for each AR are represented at different distances x along the outlet. For AR = 6, Fig. 15(a), the field (v, ŵ) displays a double vortex structure in the x = 0 plane as in the DNS velocity field of Fig. 3(b). In the central region, the perturbation points directly toward the left inlet, extending slightly beyond the section shown in Fig. 15, while the vortices are distributed in the upper and lower portion of the plane. For AR = 10, Fig. 15(b), one observes, as expected, 3 vortices like in Fig. 3(c). For both AR's, the v and ŵ modes decay slowly along the outlet branch of the junction and û changes sign at the plane y = 0 and at the height z at which a vortex center is located. For AR = 6, the absolute value of û increases from the center and is maximal at approximately x = 1 while, for AR = 10, it is highest at approximately x = 1.5. For AR = 6 and 10, one finds, like for AR = 1, that the three eigenvalues closer to σ = 0 correspond to eigenmodes with one, two or three vortices (see Fig. 9); however, only one of these three eigenvalues may become positive when Re increases and corresponds to number n of vortices depending only of AR. In Fig. 16, we compare the variations of Re c with AR (1 < AR < 10) obtained from the linear stability analysis (LSA) using Eq. 2 to those determined from DNS simulations. The values of Re c found by the stability analysis are in good agreement with those of the DNS obtained by using an increasing ramp for Re and display, like them, an initial linear increase of Re c for AR 2, leveling off for AR ≥ 5. This suggests that, at large values of AR, Re c increases only slowly with H; ref. 34 reports indeed a critical value Re c = 55 at AR = 262 for increasing flow rates and Re c = 43 for decreasing ones. The minimum of Re c is found at approximately AR = 2, so that, for Re < 23, the flow is completely segregated within the range of AR values analyzed. The results for both increasing and decreasing Re ramps compare well with those of Haward et al. [START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF] , who reported an increase of Re c for decreasing AR when AR 2. The numbers of vortices as a function of AR from the linear analysis have been superimposed in Fig. 5 onto the equivalent data points from the DNS. The results are perfectly compatible and the transition from n to n + 1 vortices can be determined more precisely by means of the linear analysis due to the reduced computation time.

Let us compare now the spatial structure of these modes to the single and multiple vortex instabilities of the nonlinear DNS results. As shown above, the leading mode that dictates the shape of the base flow instability clearly captures the number of vortices predicted by the DNS. The distances between adjacent vortices arising from the two methods have been superimposed in Fig. 7: the values obtained from the linear analysis are slightly higher (5%) that those from the DNS. One notes however that, while accurate Re c values may be expected from the linear analysis, the spacing between vortices is approximated since, in the DNS, the instability is already developed and the flow structure may be influenced by nonlinear terms.

Moreover, the null frequency component ω of the leading eigenvalues is compatible with the steady state of the flow obtained for Re > Re c by means of the direct numerical simulations. The linear stability analysis indicates indeed that the instability of the flow occurs through a pitchfork bifurcation. Additionally, the global mode structures suggest in all cases that the effect of inertia is higher in the intersection of the channels, as shown by the large values of the component v. The accommodation of the vortex flow in the downstream direction leads to a strong perturbation of the velocity component u 0 along the outlet branches. Finally, we observe that, when Re increases, there is a spatial elongation of all eigenmodes toward the outlets.

C. Sensitivity function

We investigate now the spatial variations of the sensitivity function S defined by Eq. 4 for flows in which axial vortices are present. The occurrence and location of regions of large sensitivity to local feedback forces is indeed closely related to the global mode dynamics [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF] . For instance, if variations of the flow are induced in regions of low S values, this influences very little the leading eigenvalue. Any strategy aimed at controlling the instability must therefore be applied to a region of the flow where S is large. are not shown for concision but will be briefly mentioned below. The values of S for the three isosurfaces shown are S = 0.25, 0.5, and 0.75 and the outer (inner) shells correspond to the smaller (higher) sensitivities.

For AR = 1, Fig. 17(a), S is largest inside two symmetrical regions of the crossing zone of the junction elongated towards the inlets. More precisely, the maximum of S is located on the z-axis at z/H ≈ ±0.15. The inspection of the data shows that these two points are closer to the center of the junction than the locations of the maxima of w 0 and v. We also observe that the region where the sensitivity is highest differs from the location of the vortex (z = 0). This latter feature resembles the result of Chen et al. [START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF] for T-junctions, where the sensitivity is highest in lobes located in the exterior region of the vortices. However, in the geometry of this latter work, the S lobes are elongated in the direction of the outlets instead of the inlets, like in the present work. Interestingly, the minimum of S is located at the center of the X-junction where the vortex is generated.

For AR = 6, Fig. 17(b), the sensitivity map displays one more lobe in the intersection of the channels than seen above for AR = 1 (3 instead of 2). S reaches its maximum values on the z-axis but, here, there is a local maximum at z = 0 instead of a minimum as for AR = 1. The absolute maximum for S occurs at two points close the upper and lower walls for both AR = 6 and 10 (in this latter case, S displays four local maxima). This feature may be interesting for control applications using an external actuator because the most sensitive (target flow) region is close to the boundaries and not immersed in the bulk of the fluid. Interestingly, as AR increases, the location of the maximum value for S is closer to that of the maximum of w 0 , although the latter is not involved in the definition of the sensitivity. On the other hand, for all AR ′ s, the minimum of S is located at the points where the vortices appear. Following Ref. [42], this suggests that, in order to control the instability within X-junctions, perturbations must not be applied at the centers of the vortices. Despite the differences between the distributions of S for different AR ′ s, an important feature is that, in all cases, the sensitivity S is largest in the region of the interface where the two fluids first meet. For a T-shaped junction, Fani et al. [START_REF] Fani | Investigation of the steady engulfment regime in a three-dimensional T-mixer[END_REF] also report a maximum of the sensitivity in the crossing region. There are indeed two facing inlets both in the X-and T-junctions; however, the values of Re c are lower in X-than in T-junctions.

In order to understand better the instability, we compare now the spatial distributions of the sensitivity S and of the time derivative of the kinetic energy per unit volume Ė ≡ u ′ • ∂u ′ /∂t for AR = 1, Figs. 18(a,b,c), and AR = 6, Figs. 18(d,e,f). Larger values of Ė are concentrated in two (AR = 1) and three (AR = 6) main regions (Figs. 18(b) and 18(e)) with the maxima located very close to those of the sensitivity (Figs. 18(a) and 18(d)). This emphasizes the importance of these regions for the development of the instability. Moreover, the dominant contribution to Ė corresponds to the exchange of energy between the perturbation and the base flow represented by the term v v ∂v 0 /∂y [START_REF] Lashgari | The planar X-junction flow: Stability analysis and control[END_REF][START_REF] Schmid | Stability and Transition in Shear Flows[END_REF] , as shown in Figs. 18(c) and 18(f); this is due to the strong v-component and the large gradient of the y-component v 0 of the base flow, as shown in section III B. This analysis of the energy exchange confirms therefore the important contribution of the interaction between the two facing flows in the inlets of the junction to the triggering of the instability.

IV. CONCLUSIONS

In this paper, we have studied numerically the flow structure and the instabilities creating axial vortices in X-junctions of perpendicular channels of rectangular cross sections with aspect ratios AR. Previous works had only dealt either with low values of AR (one vortex) [START_REF] Haward | Tricritical spiral vortex instability in cross-slot flow[END_REF][START_REF] Zhang | Experimental investigation of three-dimensional flow regimes in a cross-shaped reactor[END_REF][START_REF] Zhang | Trapping region of impinging jets in a cross-shaped channel[END_REF] or large ones (many vortices) [START_REF] Kalashnikov | Effect of polymer additives on ordered threedimensional structures arising in cross-slot flow[END_REF] . The present study has been focused, instead, on a transition range: 1 < AR < 11 for which 1 to 4 vortices are observed.

We first used 3D DNS simulations to determine the global structure of the flow field as a function of AR and of the Reynolds number Re = U in W /ν. Up to a critical Reynolds Re c , one has, for all aspect ratios, segregated outflows of the two fluids in the two outlets, each on a side of the mid-plane y = 0. Above Re c , steady vortex structures appear at the intersection of the junction and induce some local mixing of the fluids. Although the geometry of the domains are different, these steady vortex structures are reminiscent of the one reported by Kerr & Dold, who analyzed the stability of an stagnation point flow within an infinite domain [START_REF] Kerr | Periodic steady vortices in a stagnation-point flow[END_REF] .

For 1 ≤ AR ≤ 3.8 and Re ≥ Re c , a single vortex with the axis parallel to each outlet develops, in good agreement with the results of ref. 23. This feature has also been observed at the intersection of circular tubes in the range of crossing angles: 68 • ≤ α ≤ 90 •24 .

For AR ≥ 3.8, more vortices stacked along the z-axis appear for Re > Re c and their number n increases steadily with AR. In the studied cases, the number of vortices only depends on AR and not on Re(> Re c ). Also, the instability leading to the appearance of the vortices is always subcritical irrespective of AR: the vortices appear and disappear at different thresholds Re c (respectively Re * c ) when Re follows an increasing (resp. decreasing) ramp. Re * c is always smaller than Re c and both numbers vary with n but retain similar orders of magnitude.

Compared to the case n = 1, the flow structures for n = 2, 3, and 4 display an important difference: for n = 1, the vortex interlaces the streamlines of the two fluids across the whole section. For n > 1, one observes instead, in the outlets, n + 1 alternate stripes of the two pure fluids separated by zones close to the vortex centers where their streamlines are interlaced: mixing due to the vortices is therefore less thorough. Another important feature is that the normalized velocity profiles v(z/H)/v max,z corresponding to different AR values collapse precisely. As a result, the distance d v between adjacent vortices is proportional to H and increases therefore with AR for a given aperture W : the velocity field corresponding to each vortex is then more and more elongated until a new vortex appears.

The DNS simulations are heavily time consuming which makes difficult, for instance, the precise determination of the threshold of the instability. In order to obtain such information and understand better the dynamics of the system, we performed a global linear stability and sensitivity analysis in which the steady segregated flow is used as the base state. The critical values Re c for the transition from zero to n vortices obtained in this way agree well with those determined from the DNS by increasing Re. For Re > Re c , only one eigenvalue is both positive and real, which agrees with the idea that both zero and n vortices configurations are steady, as observed in the DNS simulations. For a given AR value, the corresponding eigenmode has the same number of vortices as that determined by the DNS, and the locations of these vortices predicted by both methods agree well.

The analysis of the leading eigenmodes also provides interesting information on the perturbation fields of the instability. They do not reach their highest values at the center of the intersection (origin of coordinates) but on the z-axis and close to the top and bottom walls (see Figs. 15 and16). The sensitivity study, for which the adjoint modes must be considered, shows that, in these spots, the receptivity to feedback local forces is highest: thus, these are the regions to actuate in order to control the instability. This result may be interpreted in terms of the kinetic energy variation with time, which reaches its maximum values close to the regions where S is maximum; this variation is mostly due to the transfer of momentum from the base state to the v component (this agreement between the locations of the maxima of S and the variation of the kinetic energy has been also reported for X-junctions but with three inlets and one outlet [START_REF] Lashgari | The planar X-junction flow: Stability analysis and control[END_REF] ). The sensitivity analysis also shows that the core of the instability is outside, and not inside, the vortex structures, as also occurs in T-junctions [START_REF] Chen | Vortex breakdown, linear global instability and sensitivity of pipe bifurcation flows[END_REF] .

As indicated above, the total computing time is significantly smaller for the global stability and sensitivity analysis than for for the DNS simulations. Typically, the CPU time to solve the non-linear problem is about eight times the required to the computation of five modes in the linear analysis (direct and adjoint problems). The linear analysis is, therefore, a robust alternative to study and predict the flow structure, and, also, for the detection of the most sensitive regions of the flow, which is a key asset for flow control strategies. However, computing the eigenmodes for large AR values requires a large amount of memory for the meshing process. This has limited up to now the values of AR which we have reached and we are currently working to overcome this limitation.

We left for future work the evaluation of the efficiency of active perturbations in the selection of one of the two steady flow configurations. This control would be achieved by, for example, suction/injection of fluid through the upper and lower walls at the intersection. The technique should be the appropriate for our geometry, because the maximum of the sensitivity is measured close to the walls where it is easier to apply an active external perturbation. This control strategy was considered by Lashgari [START_REF] Lashgari | The planar X-junction flow: Stability analysis and control[END_REF] , also for X-junctions but with three inlets and one outlet, and it was found to be successful and might be applicable in future work to our flow configuration. Another important issue to be consider in future works is the mixing efficiency of the junction, which is out of the scope of this article. As an illustration, we present a preliminary evaluation in the appendix.

As mentioned in the Introduction, junctions are potentially of interest for fluid mixing applications [START_REF] Dreher | Characterization of laminar transient flow regimes and mixing in T-shaped micromixers[END_REF][START_REF] Lobasov | Analyzing mixing quality in a T-shaped micromixer for different fluids properties through numerical simulation[END_REF] . We examined therefore briefly the mixing performance of the present Xjunctions by adding to our direct numerical simulations a transport equation and assuming a uniform concentration in each inlet (the dimensionless concentration is unity in one inlet and zero in the other). In a preliminary evaluation, we characterized mixing by the quality index M [START_REF] Engler | Numerical and experimental investigations on liquid mixing in static micromixers[END_REF][START_REF] Dreher | Characterization of laminar transient flow regimes and mixing in T-shaped micromixers[END_REF] defined by: M = 1 -σ 2 /σ 2 0 ; σ 2 is the variance of the concentration field in a section of the outlet channels, and σ 2 0 corresponds to the maximum variance in this section. For a perfectly segregated flow M is equal to zero, and for a completely mixed one M = 1.

A qualitative trend of mixing performance, for the range of AR values considered here, is shown in Fig. 19: it displays the variation of M with AR at the cross section x = 7 of one outlet channel, far from the center of the junctions. The index M decreases with AR for a fixed number n of vortices: its value is then highest when AR is near the lower limit of a range corresponding to a given n value and lowest near the upper limit. This is likely due to the fact that the height along z of a band of pure injected liquid located between two vortices (Fig. 2) is smaller near the lower limits mentioned above and larger near the upper ones. The local maximum of M decreases with n (mixing is most efficient for AR = 1 in the range considered here) while the local minimum varies less. Summarizing, when a single junction is used, increasing AR does not enhance mixing in itself, but increases the number of alternate streams of the two fluids.

The development of new microfluidic and 3D printing techniques opens the possibility to combine junctions [START_REF] Haward | Microfluidic analog of an opposed-jets device[END_REF] and to build complex 3D structures with improved mixing characteristics [START_REF] Bhargava | Discrete elements for 3D microfluidics[END_REF][START_REF] Lai | A Rubik's microfluidic cube[END_REF] . The formation of vortices and their localization in the network will be key factors of such improvements. In future studies, it will be interesting to consider as a first step the mixing properties of two junctions placed one behind the other.

Figure 1 :

 1 Figure 1: (a) X-junction of channels of rectangular cross-sections with normalized width and height W and H, respectively. The inlets are perpendicular to the outlets. (b) Perspective view of streamlines corresponding to the segregated base flow {u 0 , p 0 } discussed in Sec. II A. Streamlines corresponding to the fluid entering each inlet are labeled with a different color. (c) Streamlines of base flow {u 0 , p 0 } in the plane x = 0. The streamlines shown in graphs (b) and (c) correspond to Re = 22 and AR = 1.8

  dimensional T-shaped channel with two inlets and one outlet. This flow configuration was first studied experimentally and numerically by Engler et al. (2004), and Soleymani et al. (2008), respectively 47,48 . Fani et al. (2013) studied the linear stability of the T-junction flow by means of a spectral element method 40 .

Figure 2 :

 2 Figure 2: Streamlines obtained by DNS simulations. (a) AR = 1 and Re = 50, (b) AR = 6 and Re = 42, and (c) AR = 10 and Re = 44. The blue and orange colors identify liquids coming from each inlet. Left: 3D perspective views; Right: plane cut showing representative streamlines.

  Figure 2 shows the streamlines in X-junctions for AR = 1, 6 and 10 and Re > Re c . A key result is that, for AR > 3.8 and Re > Re c , more than one vortex appear in the flow. In this case, the incoming streams are sorted vertically in the outlets into alternating layers of the two pure fluids separated by zones where their streamlines are interlaced (these zones are located at z-coordinates similar to those of the vortex centers, as shown on the right sides of Figs. 2(b) and 2(c)). In these two latter cases, the total height of these zones of interlacement represent a smaller fraction of H than for AR = 1, as shown in Fig. 2(a) (right).

Figure 3 :

 3 Figure 3: (a,b,c) Streamlines in plane x = 0 for same DNS as in Fig. 2 with (a) AR = 1, Re = 50, (b) AR = 6, Re = 42, (c) AR = 10, Re = 44. Vertical lines: location of outlets walls. (d,e,f) Flow lines distribution for the same experiments in section x = 7 of outlet channel. Streamlines intersect the section at a right angle and appear as dots. (g,h,i) 3D views of Q isosurfaces (in red) for the same three experiments. In blue: streamlines from one entrance. (Multimedia view): development with time of streamlines in the plane x = 0 for these experiments. Re increases at first linearly with time from 1 (t = 0) to the Re values corresponding to cases (a,b,c), reached respectively for t = 100, 50 and 50, and remaining constant afterwards.

Figure 4 :

 4 Figure 4: Time sequence of maps of the x vorticity component (color scale) in plane x = 0 for the same experiment as in Fig. 3(b) with AR = 6 and Re increasing linearly with time from Re = 1 (t = 0) to Re = 42 (t = 50) and remaining constant thereafter. Black lines: streamlines. (a) t = 0, (b) t = 20, (c) t = 60 (insets: local inversion of the sign of the vorticity), (d) t = 80, (e) t = 94, and (f) t = 130. (Multimedia view): development with time of the x-vorticity in the plane x = 0 for the same three experiments as in Fig. 3.

Figure 5 :

 5 Figure 5: Number n of vortices as function of aspect ratio AR. Results from DNS simulations: ; results from linear stability analysis (LSA): (see Sec. III B for explanations).

Figure 6

 6 Figure 6 displays the variations of the y-component of the velocity along the z-axis (x = 0, y = 0) for different AR values, for which respectively n = 2 (Fig.6(a)) and n = 3 (Fig.6(b)) vortices are observed. The very good collapse of the normalized velocity profiles (v max,z is the maximum absolute value of v along the z-axis) velocity profiles in each graph implies that the normalized distance d v /H between the centers of two adjacent vortices is nearly constant with AR: one has d v /H ≃ 0.5 for n = 2 and d v /H ≃ 0.33 for n = 3. Therefore, the simulations suggest that the vortices are equidistant along the z axis; the distance d v /H between two adjacent vortex-centers is ≃ 1/n, and the distance of the top and bottom vortices to the adjacent walls is ≃ 1/2n.

Figure 6 :

 6 Figure 6: Variation of the normalized y-velocity component v/v max,z with z/H along the z-axis (x = 0, y = 0) for several aspect ratios AR and Re = 44. Profiles correspond to: (a) n = 2, and (b) n = 3. Boldface numbers: labels for vortex centers; d ij v /H = normalized distance between the centers of vortices i and j .

Figure 7 :

 7 Figure 7: (a)Normalized distances d v /H between vortices as a function of Re > Re c . Circles correspond to DNS simulations (squares for linear stability analysis, Sec. III B). For AR = 6, d v = d 12 v : ⊕ (⊞) and for AR = 10, d v = (d 12 v +d 23 v )/2: ⊗ (⊠), d v = d 13 v : • ( ). (b) Variation with AR of the mean of the dimensionless distances d i,i+1 v

  Figure 7(b) shows the variation with AR of the spacing d v of adjacent vortices (averaged over several pairs of vortices when n ≥ 2). Note that the variations of d v with AR are consistent with the estimation d v ≃ AR/n (dashed lines in Fig.7(b)) based on the order of magnitude: d v /H ≃ 1/n mentioned above.

Figure 8 :

 8 Figure 8: Maximum v max,z , along the z-axis, of the velocity component v as a function of Re. In (a), (b) and (c), the instability corresponds to the appearance of one (AR = 1), two (AR = 6) and three (AR = 10) vortices, respectively.

Figure 9 :

 9 Figure 9: Eigenvalue spectrum for the direct ( ) and adjoint (+) problems for AR = 1 and (a) Re = 40 < Re c and (b) Re = 50 > Re c . The horizontal and vertical axis correspond respectively to the real (σ) and imaginary (ω) parts of the eigenvalue. Red numbers above the leading eigenvalues are the number of vortices in the corresponding eigenmodes.

Re

2 AR = 1 Figure 10 :

 2110 Figure 10: Real part of the leading eigenvalue as function of Re for AR = 1. Vertical dashed line corresponds to the critical Reynolds number: Re c = 43.5.

Figure 11 :

 11 Figure 11: Isosurfaces for AR = 1 and Re = 50 of the û components: (a) û, (b) v, and (c) ŵ.

Figure 12 :

 12 Figure 12: Perturbation solution in yz-planes for AR = 1 and Re = 50 at several distances x along an outlet branch (from left to right: x = 0, 0.5, 1.5 and 3). The length of the vectors is given by the values of the components v and ŵ, and colors code the component û.

Figure 13 :

 13 Figure 13: Real part of leading eigenvalues versus Re for (a) AR = 6, and (b) AR = 10. Critical Reynolds numbers are respectively Re c = 39 and Re c = 42.4.

Figure 14 :

 14 Figure 14: Isosurfaces of the velocity perturbation components. Top : AR = 6, Re = 42. (a) û, (b) v, and (c) ŵ. Yellow (green) corresponds to the value -3 × 10 -6 (3 × 10 -6 ). Bottom: AR = 10, Re = 44. (d) û, (e) v, and (f) ŵ. Yellow (green) corresponds to the value -1.5 × 10 -6

Figure 15 :

 15 Figure 15: Map of the velocity perturbation components (û, v, ŵ) of the leading eigenmode in different sections x = cst. for (a) AR = 6 and Re = 42; (b) AR = 10 and Re = 44. Black arrows correspond to (v, ŵ), and the color map codes the values of û.

Figure 16 :

 16 Figure 16: Critical Reynolds numbers Re c as a function of AR: open (resp. filled) symbols for increasing (resp. decreasing) Re's. DNS: , (•); Linear stability analysis (LSA) method: ; Haward et al. 23 (inset of Fig. 2(i)): △,( ).

Figure 17 :

 17 Figure 17: Maps of isosurfaces of the sensitivity function S for (a) AR = 1, Re = 50 and (b) AR = 6, Re = 42. Isosurfaces shown correspond to S = 0.25 (yellow), S = 0.5 (orange) and S = 0.75 (red).

Figure 17

 17 Figure17displays isosurfaces corresponding to several normalized values of the sensitivity S for the cases AR = 1 (a) and AR = 6 (b) studied above; the results obtained for AR = 10 are not shown for concision but will be briefly mentioned below. The values of S for the three isosurfaces shown are S = 0.25, 0.5, and 0.75 and the outer (inner) shells correspond to the smaller (higher) sensitivities.

Figure 18 :

 18 Figure 18: Comparison between S and Ė at x = 0 for AR = 1 (a,b,c) and AR = 6 (d,e,f). (a,d): Sensitivity S, (b,e): time derivative Ė of total energy, (c,f): component v v ∂v 0 /∂y of Ė.

Table I :

 I Mesh convergence for the leading eigenvalues of the direct and adjoint problems with AR = 1, L = 7 and Re = 50.

	Mesh Grid elements λ D	λ A
	M1	10368	0.08021 0.07979
	M2	16072	0.07886 0.07988
	M3	22680	0.07828 0.07829
	M4	28800	0.07835 0.07802
		L λ D
		4 0.08102
		6 0.07887
		7 0.07889

Table II, L = 6 and L = 7 result in accurate calculations of λ and the corresponding eigenmode. Solving with L = 8 shows negligible effects on the leading eigenvalue and the corresponding eigenmode. Choosing L = 7 is a good compromise

Table II: Variation of the leading eigenvalue for different channel lengths, L, for AR = 1 and Re = 50.

ACKNOWLEDGEMENTS

This work is supported by public grants overseen by the Secretaría de Políticas Universitarias through the program 'Fortalecimiento de la Ciencia y Tecnología en las Universidades Nacionales' (reference: 03-PIO-44H), and by the French National Research Agency (ANR) through the 'Laboratoire d'Excellence Physics Atom Light Mater' (LabEx PALM) as part of the 'Investissements d'Avenir' program (ANR-10-LABX-0039). We also acknowledge support by the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) and the Franco-Argentinian International Research Project CNRS-CONICET: 'Ingénieries Vertes par la Mécanique des Fluides'. (IRP-IVMF). Data points from reference 23 have been included in Fig. 16 under the rules of the Creative Commons Attribution 3.0 License.

DATA AVAILABILITY.

The data that support the findings of this study are available within the article (and its supplementary material). Re c

APPENDIX: MIXING EFFICIENCY AND ASPECT RATIO