T-shaped micromixers aligned in a row: characterization of the engulfment regime

Accurate control of mixing between two fluids is a fundamental aspect in many applications and generally implies the use of small devices operating at low velocities. This is often achieved using micromixers which, due to the combination of small dimensions and low velocities, work in the laminar regime and generally process very limited flow rates. The flow rates can be increased using more mixers at the same time. In this respect it is appealing to use several micromixers connected in a unique device. In this paper we propose and characterize a simple strategy to connect several T-mixers together by aligning them in a row and feeding them by shared inlet channels. Since engulfment is the preferred flow regime for mixing, we investigate the proposed devices specifically focusing on the onset and on the properties of the engulfment regime. This investigation is carried out by a combined use of numerical simulation and linear stability analysis. The results reported here show that the proposed strategy can lead to compact devices in which the single mixers manifest engulfment, although some important differences in comparison with isolated T-mixers may exist depending on the spacing between the inlet/outlet channels in the device.


Introduction
Accurate control of mixing between two different fluids, which is a fundamental aspect in many applications, can be achieved by processing small quantities of fluid and by controlling residence time in a dedicated mixing device. Satisfying these requirements leads to devices which usually work in the laminar flow regime. Moreover, most of the mixing devices are passive, i.e. do not require active control nor an external power source to achieve the design performances (see, for instance, Ref. [18] for a review). The resulting problem, i.e. designing devices which work in the laminar regime are passive and lead to a controlled mixing between the two fluids entering the device, is a very active research area in the scientific community. Many devices having the above-mentioned characteristics have been proposed in the literature, and several of them are mentioned in the following review papers: [6,18,20,21,28]. Some recent promising configurations can also be mentioned here (see for instance Refs. [4,9,25,26,33]) so as to indicate the continuous research efforts in this direction. Among the large variety of micromixers which can be found in the literature, T-mixers are widely used and geometrically simple. They are made by three rectilinear conduits, two of which are parallel and aligned, and the third one is orthogonal to the first two, so that the flow enters the device through the two parallel conduits and outflows from the orthogonal one. The simplicity of their geometry makes T-mixers very attractive. Moreover, besides their geometrical simplicity, they also have good mixing properties (see for instance [3,5,11,22]). Concerning mixing, a particularly appealing flow regime that can be found in T-mixers is denoted as engulfment. In this steady regime, which can be observed when the Reynolds number Re is larger than a critical value (Re cr ), the flow is asymmetric although the geometry and the incoming flow are the colours of the streamlines indicate the inlet surface from which they are computed characterized by two mutually orthogonal reflection symmetries. As a result of the flow asymmetry, convection substantially enhances mixing, which is highly controllable due to the steadiness of the flow. Conversely, below the critical Reynolds number for the onset of engulfment, the flow shares the same reflectional symmetries of the geometry (provided that the inflow from the external environment is also symmetric) and settles in a regime denoted as vortex regime. In the vortex regime the two flows entering the mixer come into contact with the outlet channel through a planar material surface, and they mix mainly due to diffusion, which implies a generally low mixing degree. A pictorial example of the substantial difference between the vortex and the engulfment regimes illustrated above is given in Fig. 1, in which two families of streamlines are reported, starting respectively from the two inlets of a T-mixer, both in vortex (left column) and engulfment regimes (right column). For the reasons highlighted above, the ideal operative regime for T-mixers is the engulfment, and several works in the literature are focused on investigating the operative conditions for its onset (see, for instance, the review in Ref. [6]). In particular, onset of engulfment depends on the channel geometry (see Refs. [1,30,31,35,37]), fluid properties (see Refs. [15,16,22,23,29,30,34,36]), and external conditions (see Ref. [17]). The onset of engulfment has been also investigated by direct/adjoint stability analysis in Ref. [12], showing that it is related to a supercritical pitchfork bifurcation of the vortex regime.
Considering the properties of the working fluids that are typically employed in applications and the range of Reynolds numbers to operate in the engulfment regime, it results that the processed mass flow rates in a single T-mixer are generally small due to a combination of small dimensions and low operating velocities. Increasing the processed flow rates can be particularly important, as, for instance, for the miniaturization of industrial processes. Several strategies have been investigated in the literature in order to scale the flow rates of microreactors for industrial-scale applications, and most of them are reviewed in Ref. [10]. The two main strategies adopted in this respect are numbering-up, i.e. the use of parallel arrangements of identical microreactors, and sizing-up, i.e. a size increase of some among the geometrical characteristics of a single microreactor. However, for fine control of mixing the numbering-up approach is preferred. The major challenge in numbering-up strategies is to grant a uniform flow rate to the parallel microreactors. For cases in which the compounds are of high value it is feasible to equip each microreactor with a dedicated flow distributor. In other cases a shared distribution unit is employed, and several efforts are dedicated to the design of the unit (see for instance [10,32]). In this work, in the spirit of a numbering-up strategy, we focus on mixing devices which are themselves clusters of T-mixers, connected in parallel while preserving individual properties of the T-mixers composing the cluster, so as to increase almost linearly the processed flow rates using one dedicated distributor. Moreover, the use of several independent clusters, each one collecting more interconnected mixers, can be an appealing strategy for process parallelization which may show practical advantages over the use of many independent T-mixers.
The purpose of the present paper is to investigate a strategy for assembling more T-mixers in a unique device by connecting them in parallel and to characterize the behaviour of the resulting device focusing on the onset and on the properties of the engulfment regime, especially in comparison with the behaviour of a single T-mixer. The ideal design objective is to obtain a compact device while maintaining the key properties of engulfment. We propose here to align the T-mixers in a row and to feed them by inlet channels so that each inlet channel is shared by two neighbouring mixers. All conduits are rectilinear, with a rectangular cross section, and intersect forming right angles. As a result, the final devices are an alternate sequential distribution of T-joints [2] and T-mixers. T-joints are geometrically identical to T-mixers but in that case the flow enters from the orthogonal channel and exits from the two aligned channels. It is shown in Ref. [2] that the flow in T-joints is characterized by a couple of counter-rotating vortices that can be found in each one of the two outlet channels. These vortices form at the intersection between the channels and can undergo breakdown for sufficiently large values of the Reynolds number. Differently, in T-mixers the vortex and engulfment regimes are dominated by a system of vortices forming at the confluence between the two incoming flows in proximity of the outlet channel, as detailed for instance in Ref. [12]. In the devices proposed here we show that both the vortical structures typical of T-joints and T-mixers are present, and the resulting flow depends on their mutual interaction. The degree of interaction depends on the spacing between the inlet/outlet channels in the device. In particular, it is shown that when a large spacing is adopted, the two systems of vortices interact marginally, and engulfment is similar to the one for isolated T-mixers. Conversely, the interaction degree increases as the inlet/outlet channels are positioned progressively closer, and, in that case, the flow scenario can be substantially different from isolated mixers. The flow in the devices is characterized here in the case of a single working fluid, and the investigation is carried out by a combined use of direct numerical simulation (DNS) and linear stability analysis.
The paper is organized as follows: the considered geometries are defined in Sect. 2. The results, obtained numerically using the tools described in Sect. 3, are discussed in Sect. 4. Finally, conclusions are reported in Sect. 5.

Geometry of the devices
We consider here several T-mixers which are connected in parallel in a single device. We focus the attention on planar and homogeneous assemblies, i.e. on a distribution of mixers in the device made by the repetition of one periodic pattern. Furthermore, we impose that inlet and outlet channels are parallel and all conduits intersect at right angles (see Fig. 2). Each resulting device is thus made by a long rectilinear conduit on which inlet and outlet channels are alternately connected in an equi-spaced distribution. The inlet/outlet channels can be placed on opposite sides of the conduit or on the same side. The shape of the channels' cross section and their spacing are free parameters of the design. Considering the direction of the flow in the channels, the proposed devices can be interpreted as an alternate distribution of T-joints and T-mixers. The geometries investigated here are sketched schematically in Fig. 2. In particular, two combinations of cross-section shapes are considered, i.e. one in which all channels of the device have the same square-shaped cross section (with reference to Fig. 2 configurations A and C) and one in which inlet/outlet channels are rectangular (aspect ratio 1:2) and the connecting conduits are square-shaped (configurations B and D). Moreover, two arrangements of inlet/outlet channels are considered, i.e. placed on the same side (configurations C and D) or on opposite sides (configurations A and B) of the conduits connecting inlet and outlet channels. Globally, four configurations are thus considered (A, B, C, and D in Fig. 2) as a result of all the possible combinations of the options listed above. For each one of these four combinations, three different spacings of the inlet/outlet channels are investigated, whose geometrical details are summarized in Table 1. With reference to Fig. 2 we indicate with W and W g the width of the outlet and connecting channels, respectively, while both conduits have a depth equal to H z (which is the extension of the device in the z-direction). Moreover, the shortest distance between two contiguous inlet/outlet channels is indicated with L W , so that the distance between their axes is equal to L W + W . We underline that the devices considered here are indefinitely long, i.e. they comprise an infinite number of inlet/outlet channels which are alternately distributed along the x-direction. This is representative of the case in which many T-mixers are connected together in a unique device. However, the end effects, which are always present when a finite number of mixers is considered, are also investigated in Sect. 4.5.
In the following discussion all quantities are normalized using (i) the hydraulic diameter of each outlet channel, D h (D h = 4 A/P, A and P being the area and the perimeter of the cross section of each outlet channel, respectively), and (ii) the average bulk velocity in each outlet channel, U b , defined as: where I is the global mass flow rate entering the device, A is the cross section area of each outlet channel, and N o is the number of outlet channels in the device. In the case of an infinitely long device the same equation can be applied for one periodic pattern composing the device. Consequently, the Reynolds number Re adopted here is defined as follows: where ν is the kinematic viscosity of the fluid. The maximum value of the Reynolds number considered in this work is equal to 500, but generally lower values are of interest, roughly in the range Re ≤ 350. Moreover, the working flow is assumed to be incompressible.

Numerical simulations
Numerical simulations of the time-dependent incompressible flow in the considered devices have been carried out using the code Nek5000 [14], which is a spectral-element code widely used in the scientific community. The code employs a Legendre-based spectral-element formulation on cubic reference elements. Time advancing is carried out using a second-order backward differentiation formula for the diffusive terms while a second-order explicit extrapolation formula is applied for the convective terms. A splitting method based on a pressure correction step is adopted so as to advance the equations in time [13,24]. Spurious pressure modes have been eliminated using a P N − P N −2 formulation, i.e. using P N elements for velocity and P N −2 for pressure. Details on the numerical methods employed in Nek5000 can be found in Ref. [8].
For what concerns the computational domain and the related boundary conditions, in all cases only the periodic pattern composing the device has been simulated. Specifically, this pattern is represented in blue colour in Fig. 2 for each of the investigated configurations. Moreover, periodic boundary conditions are imposed on the appropriate boundaries that are clear from Fig. 2. No-slip boundary conditions are imposed on all the solid walls of the device. On the inflow boundaries a fully developed laminar flow is assumed, which depends on the shape of the inlet cross section (see, for instance, Ref. [12]). On the outlet boundaries, stress-free conditions are applied. The lengths of inlet (L i ) and outlet (L o ) channels in the computational domains used here are such that L i ≥ 7 D h and L o ≥ 15 D h , which is appropriate for studying engulfment in T-mixers, as demonstrated in Ref. [12].
The free numerical parameters of the discretization have been fixed on the basis of previous simulations of T-mixers carried out by the author and documented in the literature (see for instance [12]) and on the basis of a dedicated grid-convergence study (not reported here for the sake of brevity). In particular, we have fixed N = 6 and we have used a quasi-uniform grid (not uniform but with very gentle refinements when approaching corners) where the size of the elements is approximately equal to 0.12 W on each section of the channels. The grid is stretched along the axes of inlet/outlet conduits so that it becomes coarser towards the inlet/outlet boundaries. Globally, depending on the geometry, a variable number of elements has been used which varies between 9 · 10 3 and 1.6 · 10 4 . An example of grid used is reported in Fig. 3, where we represent a few detailed views of the grid used for case CA_2 by plotting the collocation points in case N = 9. Simulations are advanced in time keeping a constant CFL number equal to 0.4. The resulting numerical discretization has been

Stability analysis
Linear stability analysis is used here to estimate accurately and at an acceptable computational cost the critical value Re cr for the onset of engulfment. To this purpose we have used the linearized version of Nek5000 as a time-stepper for advancing the disturbance equations. The baseflow is computed by the nonlinear solver, and, when Re > Re cr , it is stabilized by taking advantage of the geometric symmetries of the problem. The linearized code is run with the same numerical discretization of the nonlinear simulations. The time-stepper has been used together with a power method for estimating the dominant eigenvalue of the linearized Navier-Stokes operator. We focus on the flow instability leading to engulfment, which has been shown in Ref. [12] to be related to a supercritical pitchfork bifurcation and is thus associated with a real-valued unstable eigenvalue.

A note on grid convergence
Even if grid convergence is not shown in detail for the sake of brevity, one result is given here so as to provide a quantitative indication on the level of grid convergence of the presented results. The test concerns the estimation of Re cr for one configuration, namely CA_2. The estimation of Re cr comprises two steps: (i) the evaluation of the baseflow and (ii) the stability analysis carried out on the baseflow itself. Thus, a grid convergence indication on the value of Re cr validates the spatial resolution for both the numerical simulation and the stability analysis at the same time. The grids used for all the simulations documented here employ N = 6, as stated above, which leads to a discretization for CA_2 involving about 4 · 10 6 degrees of freedom. With this particular grid we estimated a value of Re cr for CA_2 equal to Re cr 321.2. In order to check grid convergence, we tested an N-refinement using N = 9, thus leading to a discretization made of about 1.7 · 10 7 degrees of freedom. With this new grid, involving about 4 times the degrees of freedom of the previous one, we estimated Re cr 321.5. The negligible discrepancy between the two results versus the large disparity in the number of degrees of freedom justified the use of the coarser grid for all the other cases described in the following.

Flow peculiarities in the considered devices: a preliminary discussion
At difference with an isolated T-mixer, the proposed devices are an alternate distribution of interconnected T-joints and T-mixers. The common arms of T-joints and T-mixers are denoted here as connecting channels. T-joints are the portions of the device where inlet and connecting channels intersect, while T-mixers are the portions where the connecting and outlet channels intersect. In proximity of the T-joints, a system of vortices forms, which is new with respect to the flow scenario in isolated T-mixers. This system of vortices, denoted in the following as inlet vortices, is studied in detail in [2]. In particular, in each connecting channel, inlet vortices are a couple of intense counter-rotating vortices whose axis is parallel to that of the connecting channels themselves, i.e. the x-axis in our specific case. The main effect of the inlet vortices is to redistribute the axial momentum in the connecting channels by means of the lift-up mechanism. More specifically, their rotation sign leads to a region of higher axial velocity near the wall of the connecting channels that is opposite to the side on which inlet channels are attached. This is illustrated in several Figures in the following discussion; for the moment see Fig. 4 as concerns the axial velocity distribution and Fig. 5a-c for the identification and position of inlet vortices.
The intensity of the inlet vortices progressively decreases along the connecting channels due to diffusion and annihilation with the vorticity of opposite sign generated on the solid walls. It is thus reasonable to expect that, for sufficiently long connecting channels, i.e. for a sufficient distance between inlet and outlet channels, inlet vortices become weak, the axial velocity in the connecting channels approach that of a fully developed channel flow, and, consequently, the device behaves as a set of independent T-mixers. Conversely, if the connecting channels are not long enough, the intensity of the input vortices is still large when they enter the T-mixers. In this case, the flow scenario in the outlet channels of the device can be substantially different from that of an isolated T-mixer for two main reasons.
Firstly, the velocity distribution of the flow entering each outlet channel can be substantially different from that of an isolated mixer because inlet vortices generate a significant redistribution of axial momentum in the connecting channels, as mentioned above. This is illustrated quantitatively in Fig. 4. As shown in the Figure, which refers to configuration A, the secondary velocity induced by input vortices, whose clear trace can be observed in Fig. 4c, d, can substantially affect the momentum distribution in proximity of the outlet channels. Conversely, for L w = 3, Fig. 4b shows that the effect of inlet vortices is almost vanishing and the resulting velocity distribution is close to that observed in the isolated equivalent T-mixer, Fig. 4a.
Secondly, input vortices can interact directly with another system of vortices, denoted here as confluence vortices. Confluence vortices, which are typical of isolated T-mixers, form in the recirculation region that can be found at the confluence of the two incoming flows in proximity of each outlet channel. We refer again to Fig. 5 for the identification and shape of the confluence vortices in one particular device. It is shown in Ref. [12] that confluence vortices play a fundamental role for the onset of engulfment in isolated T-mixers. For this reason, interaction between confluence and inlet vortices is expected to substantially alter the onset and the characteristics of engulfment in the devices.
All the aspects that have been qualitatively introduced here so as to preliminarily clarify the peculiarities of the problem at issue will be discussed in quantitative terms in the next Sections.

Critical Reynolds numbers for the onset of the engulfment regime
We identify here the critical Reynolds numbers (Re cr ) for the onset of the engulfment regime. In each configuration, the value of Re cr has been estimated by a dedicated linear stability analysis. More specifically, symmetric baseflows (vortex regime) and their stability properties have been evaluated for a set of Reynolds numbers across the stability limit, and, successively, the value of Re cr has been estimated by interpolating the computed growth rates associated with the dominant eigenvalues versus the corresponding value of Re. We remind that Re cr is associated with a null growth rate of the leading eigenvalue, which is thus marginally stable. The results are reported for all the considered configurations in Table 2. In the same Table we report the values of Re cr , both recomputed here and reported in the literature, for isolated T-mixers having the same shape of those composing the devices. In this way it is possible to appreciate the variation in terms of Re cr observed when a given T-mixer is collocated in the device. The value of Re cr for isolated T-mixers has been recomputed using the same numerical discretization adopted for the corresponding device for two reasons: (i) to provide a reference term for comparison obtained with the same numerics, and (ii) to validate the numerical discretization used here against the literature. Concerning this second aspect, Table 2 shows that the results provided by our analysis are in reasonable agreement with the literature. As a first comment, it can be noted that Re cr always increases in the considered configurations in comparison with the corresponding value for isolated T-mixers, with the only exception of case CD_1, where Re cr slightly decreases, but relative variations are very small (of the order of 3%). Moreover, the variation of Re cr is generally weaker when the spacing between the inlet and outlet channels, L w in Fig. 2, is the largest one, while they become progressively more important as L w is decreased. This is in agreement with the qualitative discussion in Sect. 4.1. In particular, for sub-configurations 1 (CA_1, CB_1, CC_1 and CD_1), i.e. those associated with the longest value of L w (L w = 3 W ), variations of Re cr in comparison with the isolated case are modest: 5.2% for CA_1, 3.9% for CC_1 and −3.1% for CD_1. Only for case C B_1 the variation is more significant, specifically of the order of 18.1%. Note that sub-configurations 1 lead to quite compact devices, as shown for instance in Fig. 2e, f where CA_1 and CB_1 are plotted in scale; thus, the fact that Re cr is modestly affected is a positive aspect from a practical viewpoint if a compact parallelization of T-mixers is searched.
As concerns the effect of L w on Re cr , in all cases Re cr significantly increases as L w is decreased. For instance, when L w = W (which corresponds to very compact devices) variations of Re cr with respect to the  isolated T-mixer are: 26.4% for CA_2, 42.1% for CB_2, and 48.2% for CD_2. Instead, for case CC_2 the flow is fully stabilized on the vortex regime, and engulfment has not been detected for Re ≤ 400. The same behaviour, i.e. the absence of the engulfment and a consequent large stability of the vortex regime, has been observed for all the configurations in case of an extremely compact design, i.e. when L w = 0.5 W (sub-configurations 3). The only exception is case CC_3, for which we found that the vortex regime undergoes a bifurcation very soon, precisely at Re cr = 77.6 (we remind that Re cr 254.0 for the isolated mixer). However, the flow regime observed in this specific case, when Re > 77.6, is very different from a typical engulfment regime and does not lead to a significant improvement in mixing although reflectional symmetries are broken. This regime is described in detail in Sect. 4.3.3.

Analysis of the flow in the considered configurations
We now investigate the flow in the considered configurations. The discussion is carried out separately for each configuration.

Configuration A
In configuration A all channels have the same cross section, and inlet and outlet channels are connected on opposite sides of the connecting channels. In the vortex regime the flow coming from each inlet channel distributes equally between the two connecting channels departing from it. Consequently, the bulk velocity is equal in inlet and outlet channels, and it is the double of that measured in the connecting channels. As a result, as concerns the mixers of the device, this is an accelerating configuration, i.e. the flow accelerates when entering the outlet channels. In Fig. 5 we show a 3D view of one periodic pattern of the device, together with the vortical structures identified by the λ 2 criterion [19] and some representative sections of the flow. Inlet and confluence vortices are explicitly labelled in Fig. 5c, and the secondary flow that they induce is represented by in-plane velocity vectors in sections S1 (subfigure (a)) and S2 (subfigure (b)). As concerns inlet vortices, their intensity in the formation region (see Fig. 5a) is large, as highlighted by the vorticity contour maps, and progressively decreases along the connecting channel so that, as shown in Fig. 5c, they are almost vanishing near the entrance of outlet channel.
Since there is not a direct interaction between inlet and confluence vortices, the engulfment in CA_1 is similar to the case of isolated T-mixers and is solely related to the confluence vortices and to the separation region where they take origin. This is shown in Fig. 6, where it is clear that engulfment originates from a tilting of the separation region at the confluence between the two incoming flows (see Fig. 6a) and results in an asymmetric distribution of vorticity among the confluence vortices so that, already at a short distance from the beginning of the outlet channel, the secondary flow is dominated by two co-rotating vortices (see Fig. 6c) leading to a significant increase of mixing between the two incoming streams.
When configuration CA_2 is considered, in which L w = W , the flow scenario is substantially different. Focusing firstly on the vortex regime, Fig. 7b shows that now the intensities of inlet and confluence vortices at the entrance of the outlet channel are comparable. Moreover, as shown in Fig. 4c, the distribution of velocity at the entrance of the outlet channel is substantially affected by the inlet vortices. Figure 7b also shows that a couple of secondary weaker vortices forms in the connecting channel and merges with the confluence vortices. As a result, eight vortices enter the outlet channel as shown in Fig. 7c; the four main vortices are the inlet vortices, while at the centre of the section other four weaker vortices are present, and they are the result of the union between the confluence and the secondary vortices. Further downstream in the outlet channel only the four inlet vortices remain as they are more intense, while the confluence vortices annihilate with the vorticity of opposite sign of the adjacent inlet vortices, as shown in Fig. 7d. As a result Fig. 7d, which at first glance seems similar to Fig. 5e, is very different indeed since the vortices dominating the secondary flow are almost in the same position in the two Figures but they rotate in opposite directions. This happens because in CA_1 (Fig. 5e) these vortices are the confluence vortices while in CA_2 (Fig. 7d) they are the inlet vortices, so that in the two cases vortices have a different origin and opposite sign.
As a result of the dominance of inlet vortices, engulfment in CA_2 is also substantially different in comparison with an isolated T-mixer, as shown in Fig. 8. In particular, Fig. 8a shows that confluence vortices and the related tilting of the separation region in x − z planes is still the driving mechanism for the onset of the engulfment. This is clear in Fig. 8c, where again it is shown that eight vortices enter the outlet channel (as in Fig. 7c). However, due to the symmetry breaking, the four inlet vortices are almost symmetrical, similarly to the vortex regime (see Fig. 7c), while the four confluence vortices are markedly asymmetric, with two corotating cores which are as intense as the inlet vortices and other two co-rotating vortices which are very weak and hardly visible in the Figure. As a net result, further downstream (see Fig. 8d), the vorticity redistributes so that only four vortices remain, made by two couples of counter-rotating vortices, which, however, are in asymmetric configuration with respect to the x − y and x − z planes. They mainly collect the vorticity of the inlet vortices, but they are in an asymmetric configuration due to the velocity induced in the first part of the outlet channel by the asymmetric confluence vortices. As a global effect, since the symmetry is broken, convection plays a role in the mixing process between the two streams entering the outlet channel.
Finally, when L w is further reduced, i.e. for case CA_3 where L w = 0.5 W , the scenario is similar to that observed in the vortex regime for CA_2 in Fig. 7, with the difference that inlet vortices are now even more intense at the confluence region. As a result, the role of the confluence vortices is further reduced, and the vortex regime remains stable, at least up to Re = 500, which is the maximum value simulated here.

Configuration B
We now focus on configuration B which, at difference with configuration A, has inlet and outlet channels with rectangular cross sections (aspect ratio 1 : 2), whose area is the double of that of the connecting channels. Thus, in the vortex regime the bulk velocity is constant in all the conduits of the device, and the configuration is thus non-accelerating.
The main difference between configurations A and B concerns the evolution of the inlet vortices. More specifically, in configuration B the intensity of the inlet vortices is initially lower than in A due to the differences in the geometry. From a quantitative viewpoint, the intensity of inlet vortices on the symmetry plane of the inlet channel in configuration B is approximately equal to 0.4 times that in configuration A. As a consequence, the intensity of inlet vortices at the entrance of the outlet channel is lower for a fixed value of L w for configurations B in comparison with A. At the same time, the qualitative behaviour of the flow as L w is reduced is analogous to that discussed for configuration A. In particular, as L w is reduced there is a progressively stronger interaction between inlet and confluence vortices. When input vortices become dominant the vortex regime is stabilized. As concerns the stabilized vortex regime in configuration CB_3, the flow scenario is very similar to configuration CA_3. Also in configuration CB_3 the flow in the outlet channel is dominated by eight vortices, four inlet and four confluence vortices, even if their intensity and position is different in comparison with CA_3. This is shown, for instance, in Fig. 10 where vorticity and in-plane velocity fields at the same y-section (y = 0.55 D h ) are reported together for CA_3 (left Subfigure) and CB_3 (right Subfigure). Concerning the engulfment regime and the interaction between inlet and confluence vortices, in configuration CB_2 confluence vortices are still dominant, consistently with the fact that inlet vortices are weaker in configuration B than in configuration A. As an example, we report in Fig. 9 some representative views of the engulfment regime for CB_2. As can be noted in the Figure, tilting of the confluence vortices is evident as for the equivalent isolated T-mixer (not reported here for the sake of brevity). Moreover, secondary flow in the outlet channel is dominated by two co-rotating confluence vortices. Figure 9 clearly shows that in configuration CB_2 confluence vortices are still dominant in the flow dynamics, i.e. both for the onset of the engulfment and for the secondary velocity field in the outlet channel (Fig. 10).

Configuration C
Configuration C has the same geometric characteristics of configuration A with the difference that in C inlet and outlet channels are positioned on the same side of the connecting channels. The main effect of this peculiar arrangement is that inlet vortices lie close to the wall of the connecting channel that is opposite to the side on which outlet channels are attached. This fact has two consequences. Firstly, the secondary motion induced by inlet vortices tends to accelerate the axial velocity in the region where the confluence vortices originate. This aspect is highlighted in Fig. 11, where the normal vorticity and the in-plane velocity fields are reported for CA_2 and CC_2 on a section of the connecting channels. The two sections in Fig. 11 have been arranged so that in both cases the outlet channels are positioned on the bottom side of the Figure. We remind that, according to previous investigations in the literature and, in particular, to Ref. [12], the region close to the wall which is opposite to the outlet channel contains the core of the engulfment instability, and, thus, modifications of the flow field in that region are expected to have a substantial effect on the onset of engulfment. Thus, for a fixed value of L w , it is expected that configuration C manifests a larger effect on engulfment than observed in configuration A. This is confirmed by comparing the values of Re cr for the two types of devices reported in Table 2. Indeed, as can be noticed from Table 1, already for CC_2 (L w = W ) engulfment is fully suppressed. Moreover, a completely new flow regime is observed further reducing L w , specifically in configuration CC_3, which manifests at Re cr 77.6. Before entering in a detailed analysis of this new regime, we would like to underline a second major difference between configurations A and C. More specifically, besides their position, also the sign of inlet vortices is opposite in comparison with those of configuration A, as shown in Fig. 11 by the sense of rotation indicated by the in-plane velocity vectors. We thus deduce that in configuration C, at difference with configurations A, the couples of inlet and confluence vortices which approach at the beginning of the outlet channel have the same sign (in configuration A they were of opposite sign). Thus, they merge gaining intensity as they enter the outlet channel. This peculiarity is present and even more evident in configuration D, and it will be illustrated in the discussion of that case in Section 4.3.4 (for the moment see Fig. 13).
Returning now to the peculiar flow configuration found in CC_3 when Re > 77.6, we report in Fig. 12 some representative views of the flow. As can be evinced by the Figure, when Re > 77.6 the flow sets in a new regime in which most of the fluid coming from one inlet channel enters only one among the two nearest outlet channels. In this way, each outlet channel receives fluid mainly from only one inlet channel. This aspect is clearly illustrated by the in-plane velocity vectors on the symmetry plane z = H z /2, reported in Fig. 12a. Due to this peculiar flow configuration, there is not any more a separation region generated by two opposite streams in proximity of each outlet channel. Consequently, confluence vortices, which are generated in that region, are not present at all in this flow regime. Thus, only inlet vortices are present, and this is clear from the vortical structures identified by the λ 2 criterion reported in Fig. 12b. In particular, in that Figure it is shown that each outlet channel receives the couple of counter-rotating inlet vortices generated in only one among the two adjacent inlet channels. The resulting secondary flow in the outlet channel is thus strongly asymmetric, as shown in the two representative sections in Fig. 12c, d, where the in-plane vectors highlight the presence of the two inlet vortices. In this case secondary motion is strong and asymmetric enough to induce a visible redistribution of the axial velocity in the first part of the outlet channel, leading to a higher velocity near the wall that is opposite to the inlet channel where the input vortices originate. However, since the flow in the outlet channel comes mainly from only one among the two adjacent inlet channels, although secondary motions in the outlet channel are intense and strongly asymmetric, mixing remains very poor.

Configuration D
Configuration D has the same characteristics of configuration B but input/output channels are installed on the same side of the connecting channels. As for the comparison between configuration A and B, also in case D the inlet vortices are markedly less intense then in C due to the different shapes of the channels. As a result, interaction between inlet and confluence vortices, which leads to modifications of the flow regimes in the outlet channels with respect to the isolated case, are weaker than in configuration C for a fixed value of L w . This can be observed by the behaviour of Re cr vs L w in Table 2. In particular, configuration D manifests engulfment also for CD_2, while this is fully suppressed in CD_3. Moreover, at difference with configuration C, no other regimes are observed except for the vortex and the engulfment ones.
As concerns the concordant rotation sense of inlet and confluence vortices and their constructive sum in the confluence region between inlet and outlet channels, which was also mentioned for configuration C, we show an example in Fig. 13. Figure 13 shows some representative views of the vortex regime in CD_2. It is clear from Fig. 13a that the couples of inlet and confluence vortices which merge have the same rotation sense. The net result in the outlet channel, depicted in Fig. 13c, is that only four vortices are observed which are the result of the constructive union between inlet and confluence vortices.

Assessment of mixing efficiency
In this Section the mixing capabilities of some of the proposed configurations are investigated. In particular, cases with square-shaped inlet/outlet channels have been considered. In order to carry out a meaningful comparison among different configurations and against the case of an isolated mixer, we investigate the mixing efficiency by fixing the flow Reynolds number to Re = 350. At this Reynolds number engulfment is present in all configurations with square inlet/outlet channels that may experience this flow regime, i.e. configurations CA_1, CA_2, and CC_1. In order to have a reference for comparison, mixing efficiency has been investigated also for the equivalent isolated T-mixer. Mixing efficiency of each mixer has been measured by simulating the dynamics of a passive scalar injected from one inlet of the mixer. This has been done using the same code employed for the fluid simulations, but using a dedicated grid which is more refined that those used for the flow, uniform, and built ad hoc for the simulations of the scalar dynamics. The numerical spatio-temporal resolution used for simulating the dynamics of the passive scalar is the same for all the four configurations considered, so that differences observed cannot be ascribed to differences in the numerics but are entirely due to the differences in the flow velocities. In simulating the dynamics of the passive scalar it is assumed that its non-dimensional diffusivity is equal to 1.0 · 10 −3 , which is a reasonable value to highlight the differences between the different configurations within the first 15 hydraulic diameters of the outlet channels. The scalar distribution is characterized by the scalar concentration Φ, which ranges from 0 (no scalar) to 1 (100% scalar presence). Once the steady distribution is reached, following the literature [7], efficiency δ m (y 0 ) is computed versus the distance y 0 from the beginning of the outlet channel, where δ m is defined as follows: where and Quantity v(x, y, z) in the above definitions indicates the distribution of the velocity y-component in the outlet channel.
A further characterization of the role of convection for mixing has been carried out by following a large set of Lagrangian particles which are injected from one specific inlet channel. This test has been carried out for each configuration. As a result, by plotting the particle trajectories it is possible to highlight the sharp boundary that would confine a passive scalar introduced by the same inlet channel in case of null diffusivity. This picture helps in understanding and interpreting the role of convection on mixing, especially in comparison with the effective scalar distribution obtained for a fixed value of diffusivity (here equal to 10 −3 , as already stated above). The mixing efficiency δ m computed in outlet channels of the considered configurations is reported in Fig. 14. As shown in the Figure, mixing capabilities of the configurations proposed here are always higher than that of an isolated T-mixer at a sufficient distance from the beginning of the outlet channel, i.e. roughly for y 0 > 5. Conversely, in the first part of the outlet channel (y 0 < 5) the proposed configurations perform worse than the isolated mixer. Actually, the region y 0 < 5 is not important for the mixers since it is only the very first part of the outlet channel, which is usually much longer in real applications. Nevertheless, it is interesting to highlight this aspect in commenting the results reported in Fig. 14. The reason for this behaviour is that the natural vortices in the isolated case are stronger, and they form closer to the top wall of the mixers, as highlighted in the previous Sections, this contributing to mixing substantially before entering the outlet channel. However, for y 0 > 5, in the isolated case the secondary flow is dominated only by the natural vortices which lead to the accumulation of the scalar on a rather confined region of the channel, thus limiting its diffusion. Conversely, the other mixer configurations, especially CA_2, have a more complex set of vortices contributing to the secondary flow, which in all cases distributes the scalar in different patterns which are advantageous for mixing. These aspects are highlighted in Figs. 15 and 16, showing the distribution of the passive scalar at sections y 0 = 2 and y 0 = 10, respectively. In each of the two Figures we report on the top row the distribution of the scalar concentration and on the bottom the Lagrangian particles entering from two contiguous inlet channels coloured in blue and red, respectively, which highlight the interface that the passive scalar would have if its diffusivity would be null. The combined analysis of the two rows of Figures helps in understanding the role of convection in the mixing process.
Concluding, the configurations proposed here generally perform better than the case of an isolated T-mixer as concerns the mixing efficiency along the outlet channels. Among the different arrangements, the one with inlet/outlet channels positioned on opposite sides of the connecting conduits (CA_1 and CA_2) performs better. Moreover, mixing efficiency is increased as the distance between inlet and outlet channels is reduced, provided that engulfment may be still experienced by the device. As concerns the pressure drop across the proposed devices, this has been evaluated considering a tract of inlet channel whose length is 7 D h and a tract of outlet channel which is 14 D h long. The device with the largest pressure drop is configuration CA_1, for which ΔP/(ρ U 2 b ) 3.45. Configurations CA_2 and CC_1 have lower pressure drops, equal to 0.92 and 0.97 times the one experienced for CA_1, respectively. Finally, although these results have been obtained for the mixers with square-shaped inlet/outlet channels, we believe that they are indicative also for configurations with rectangular cross sections.

Assessment of end effects
In order to check the representativeness of the results obtained for indefinite configurations when a finite-size device is considered, one configuration has been investigated, namely CB_2, in the case in which the device is made of only five T-mixers in a series, implying five outlet and six inlet channels. The number of inlet/outlet channels in the device has been chosen as a compromise between computational costs and representativeness of the investigation. End effects, which might be significant, have been mitigated by setting the inflow mass of the two inlet channels at the extremities of the device to be 1/2 of that entering in the remaining ICs. In this way the averaged bulk velocity in the inlet and outlet channels has the same value For the considered finite-size device we have computed the value of Re cr for the onset of engulfment in the same way followed for the other configurations, i.e. by a combined use of DNS and linear stability analysis. The resulting value is Re cr 200.6. We remind that for the indefinite configuration CB_2 we found

Re cr
197.9 Thus, the discrepancy between the two values is approximately equal to 1.3% and is expected to be even lower for finite-size devices including a larger number of inlet/outlet channels.
In order to characterize the engulfment regime in the finite-size device, in Fig. 17 we report representative views of the overall flow. As illustrated in the Figure, apart for the extremities where end effects are appreciable (which, however, are limited to the sole two outlet channels at the device extremities) the engulfment regime in the remaining outlet channels is very similar to the one observed for CB_2. Moreover, the top view in Fig. 17a shows that engulfment manifests as usual by a pronounced tilting of the confluence vortices in their formation region. Vortical structures highlighted in Fig. 17b have the same typical shape observed for the corresponding indefinite configuration as in the previous Section. Globally, the illustrated flow features (compare Fig. 17 with Fig. 9), together with the estimation of Re cr described above, demonstrate that the investigation carried out simulating one minimal periodic pattern of the device is a reasonable base study to characterize finite-size devices.

Conclusions
In this work we have investigated one strategy for using multiple T-mixers which are connected in parallel inside one single device. The channels of the device are rectilinear, they intersect at right angles, and they are parallel to a common plane. We assume to have indefinitely long devices made by a periodic distributions of one fixed pattern. The resulting configurations can be seen as alternate distributions of T-joints and T-mixers. Several configurations of this kind have been studied numerically here combining different cross-sectional shapes and relative spacing of the channels, and considering two different arrangements of inlet/outlet channels. The objective of the analysis is to characterize the vortex and the engulfment regimes in such devices in comparison with the same regimes observed for isolated T-mixers. It is shown here that the flow in the considered configurations is characterized by two systems of vortices. The first one, made by what we denoted here as inlet vortices, forms at T-joints, while the second system, made by confluence vortices, is generated in the T-mixers. The resulting flow depends on the degree of interaction between these two systems of vortices. When the interaction degree is low, i.e. when inlet and outlet channels are sufficiently distant, the T-mixers of the device behave almost as independent mixers. Conversely, when inlet and outlet channels are sufficiently closed, the inlet and confluence vortices strongly interact leading to scenarios which can be substantially different from those of isolated T-mixers. In particular, as L w is decreased, the engulfment regime shows peculiar characteristics while the value of Re cr for its onset increases until the vortex regime is completely stabilized and engulfment suppressed. The mixing properties of some among the proposed devices have also been assessed, i.e. of those with square-shaped inlet/outlet channels. In this respect it has been shown that mixing efficiency obtained for a sufficiently long outlet channel is systematically higher than the one obtained by an isolated T-mixer. Moreover, efficiency increases as the spacing of the inlet/outlet channels is decreased, provided that engulfment is not suppressed. We finally show that the results obtained for indefinitely large configurations also apply when a device of finite size is considered. In the specific example provided here we consider a device with six inlet and five outlet channels.
According to the present analysis it is possible to design compact devices collecting several mixers at the same time while maintaining their properties similar to the isolated case and thus using the knowledge already available in the literature for isolated mixers. The use of such devices represents a strategy to scale the processed mass flow rate in the mixing process. Although there is an extremely wide variability in choosing a parallel arrangement of the mixers, we believe that the ones proposed here are prototypical for their simplicity. Thus, the results provided can already support the design of several devices of practical interest.