Direct numerical simulation of waterethanol flows in a Tmixer
Abstract
The efficient mixing of fluids is key in many applications, such as chemical reactions and nanoparticle precipitation. Detailed experimental measurements of the mixing dynamics are however difficult to obtain, and so predictive numerical tools are helpful in designing and optimizing many processes. If two different fluids are considered, the viscosity and density of the mixture depend often nonlinearly on the composition, which makes the modeling of the mixing process particularly challenging. Hence waterwater mixtures in simple geometries such as Tmixers have been intensively investigated, but little is known about the dynamics of more complex mixtures, especially in the turbulent regime. We here present a numerical method allowing the accurate simulation of twofluid mixtures. Using a highperformance implementation of this method we perform direct numerical simulations resolving the spatial and temporal dynamics of waterethanol flows for Reynolds numbers from 100 to 2000. The flows states encountered during turbulence transition and their mixing properties are discussed in detail and compared to waterwater mixtures.
keywords:
Direct numerical simulation, miscible fluids, waterethanol, turbulent flows, Tmixer, mixing, transition to turbulence1 Introduction
The bioavailability of pharmaceutical drugs is the degree to which a drug becomes available to the target tissue after administration. It is associated with the aqueous dissolution rate of the drug and determined by the size of its constituent particles. Nearly 70% of the drugs developed in recent years show poor water solubility and hence manufacturing processes to produce monodisperse particles of the smallest possible size have been the focus of much research. These processes can be either topdown or bottomup and usually deal with the production of nanocrystals containing mainly active pharmaceutical ingredient. Liquid antisolvent (LAS) precipitation is a promising, yet poorly understood, bottomup process. Thorat and Dalvi (2012) recently provided a review of experimental and numerical approaches to LAS precipitation, highlighting its paramount importance for the pharmaceutical industry. Its advantages are simple upscaling, economical production and the generation of crystals in the nanometer range. In LAS precipitation the drug is initially dissolved in an organic solvent such as ethanol, which is then mixed with an antisolvent as water. In the mixing zone a thermodynamic driving force (supersaturation) arises and the drug precipitates into nanoparticles. Their size is determined by the supersaturation, which in turn strongly depends on the mixing efficiency. Consequently, a detailed understanding of the fluid dynamics of such mixing phenomena and its connection to particle formation is required to control and optimize the outcome of LAS.
(a)  (b) 
Despite the growing attention that LAS precipitation has attracted in the last decade, few works have focused on the accurate modelling of the mixing of the solvent (alcohol) and the antisolvent (water), see Orsi et al. (2013); Wang et al. (2012). Waterwater (WW) mixtures in laminar and slightly transient regimes in Tshaped micromixers have been intensively investigated Hoffmann et al. (2006); Bothe et al. (2006, 2011); Fani et al. (2013, 2014); Orsi et al. (2013). Fig. 1(a) shows a schematic of a simple Tmicromixer as usually used to investigate the physics of mixing processes. The most important parameter of the system is the Reynolds numbers , where is the mean inlet velocity, the diameter of the inlet and the kinematic viscosity of the fluid. The Schmidt number , where is the diffusion coefficient, quantifies the rate of momentum to mass transfer is typically large in fluids. The flow patterns and corresponding mixing behaviour have been thoroughly characterised as a function of , and an evident improvement in mixing efficiency has been observed for increasing the Reynolds number Hoffmann et al. (2006); Bothe et al. (2006, 2011); Fani et al. (2013, 2014); Orsi et al. (2013). In addition, direct numerical simulations (DNS) of turbulent mixing, complemented with the simulation of appropriate population balance equations, have succeeded in predicting the experimental outcome of the precipitation of inorganic compounds such as barium sulphate Gradl et al. (2006); Schwertfirm et al. (2007).
At very low the fluids streams flow side by side along the main channel. In the case that the concentration is a passive scalar, as in a waterwater system, the velocity field and concentration field are topdown and leftright symmetric. Here mixing occurs only by diffusion at the straight contact plane and completely mixed state is reached after the distance . However, is typically very small and thus in applications is prohibitively long. The simulations of Bothe et al. (2006, 2011) and experiments of Hoffmann et al. (2006) showed that by increasing vortex pairs arise at the junction, see the top panels in Fig. 2. However, in this socalled stratified regime the vortices are topdown and leftright symmetric. Hence a straight contact plane between the fluids persists and the mixing remains purely diffusive. By further increasing the Reynolds number they identified an additional steady regime, in which the leftright and topdown symmetries are broken and the resulting vortices intertwine the two inlet streams as shown in Fig. 2. In this engulfment regime the contact surface between the fluids becomes convoluted and henceforth the mixing efficiency is enhanced. Note that the engulfment flow remains symmetric with respect to a 180 rotation. Fani et al. (2013) analysed the instability mechanism and Thomas and Ameel (2010) experimentally reported that at even higher Reynolds numbers the flow becomes periodic and afterwards chaotic. A stability analysis of these secondary transitions was carried out by Fani et al. (2014), who precisely determined the Reynolds numbers at which they occur. Both Thomas and Ameel (2010) as well as Fani et al. (2014) discussed qualitatively the enhancement of mixing because of transient flow.
(a)  (b)  (c) 

The first investigation of the mixing of water and ethanol (WE) was performed by Wang et al. (2012), who conducted LIF measurements in a Tmixer. By studying the LAS precipitation of curcumin they found a strong correlation of the mean particle size and the mixing efficiency, whereby a higher mixing quality resulted in a smaller mean particle size. Those findings motivated Orsi et al. (2013) to precisely study the mixing of water with ethanol in the steady flow regime () with numerical simulation. They ascertained that flow regimes analogous to the waterwater system are observed, however the transitions are hampered by the increase of viscosity as the two fluids mix. The strong dependence of the viscosity on the concentration of the mixture is shown in Fig. 1(b).
Schwertfirm et al. (2007) and Gradl et al. (2006) demonstrated that DNS combined with appropriate models of the precipitation kinetics can quantitatively predict the precipitation of inorganic compounds. These elucidations outline the motivation of this work, namely to develop a code for DNS of miscible fluids with varying physical properties at operating Reynolds numbers. Such a code would be the first step toward a predictive tool for the outcome of LAS precipitation. Note however that real operating conditions for LAS precipitation are typically beyond , which has not been reached with DNS for fluids with constant physical properties, let alone complex mixtures such as WE. The need for the DNS arises from the naturally small molecular diffusion compared to the momentum diffusion, characterized by the Schmidt number , and the strong nonlinear change of the viscosity with concentration. Under these circumstances a direct application of LES and RANS additionally would bring along with the turbulent diffusion an unquantifiable parameter to this already challenging problem. Our work may serve as a benchmark for the development of industry accessible tools such as RANS and LES.
2 Modeling assumptions and governing equations
We here consider miscible fluids with physical properties, such as density and dynamic viscosity that depend on the normalized concentration of the mixture, . Fitted functions based on experimental data Dizechi and Marschall (1982) give smooth dependencies of the density and the dynamic viscosity as a function of the composition, see Fig. 1(b). The transport of species (without chemical reactions) is modeled here with a convectiondiffusion equation, which in conservative form reads
(1) 
where is the fluid velocity. The mixing process is considered to be isothermal, as justified in detail by Orsi et al. (2013). Furthermore, the volume excess naturally occurring while mixing ethanol and water is neglected. Although the density depends implicitly on time , we assume that the total mass is conserved. This is enforced by setting
(2) 
which is known as the quasiincompressibility assumption. A detailed justification of the validity of this hypothesis can be found in Joseph (1990); Joseph et al. (1996). The fluid motion is governed by the Navier–Stokes equations, which in conservative form read
(3) 
The dynamics of the fluid motion is governed by the Reynolds number, defined as , where , , are the mean velocity in the inlet, diameter of the inlet, and kinematic viscosity of water, respectively. Note that as the surface area of the mixing channel is equal to the sum of the inlet channel, the mean velocity of inlet and outlet coincides. The diffusion coefficient is fixed to , which yields a Schmidt number for water and for ethanol. The Peclet number is . The diffusion coefficient of waterethanol mixtures depends on the concentration Legros et al. (2015), and this dependence should be considered in subgrid models of diffusion controlled micromixing. Here we are interested in the macroscopic mixing, which is controlled by convection already at moderate Reynolds numbers (). We hence assume that the dependence of the diffusion coefficient on should have a minor impact on the macroscopic properties of mixing and so we neglect it. However, we note that it would be straightforward to incorporate a concentrationdependent diffusion coefficient in our model and numerical code.
2.1 Boundary conditions
The velocity profile at the inlet is taken as fully developed laminar flow in a square duct and can be obtained by solving the twodimensional Poisson equation . The right hand side with pressure gradient is then scaled to the corresponding Reynolds number. The concentration is set to (water) and (ethanol) at the inflow planes. The noslip boundary condition at the wall is employed for the velocity field, whereas a zeroflux boundary condition is used for the concentration. At the outlet plane a mass preserving convective boundary condition is used, defined as
where the advective velocity is set to be the mean velocity . Furthermore a zero flux boundary condition at the outlet plane is chosen for the concentration. A convective boundary condition was also tested and found not to influence the concentration field, while hindering convergence of our iterative solvers.
2.2 Diagnostic quantities
We quantify the mixing efficiency with the degree of mixing
(4) 
where is the maximum of the variance achieved by completely segregated streams, and is the mean deviation of the concentration profile as
(5) 
where is the mean concentration of a cross section. The energy input is here quantified with the pressure loss
(6) 
where and are the surfaceaveraged pressure at the outlet and inlet cross section, respectively. Using the pressure loss the Darcy friction factor is calculated as
(7) 
where is the total length of the Tmixer and is the hydraulic diameter of the mixing outlet channel. Here is chosen to be the mean density of the mixture.
We characterise the relative importance of the acting transport mechanisms by considering three key time scales: the advective time scale , the viscous time scale and the diffusive time scale . Note that these can differ by several orders of magnitude depending on the flow state.
3 Numerical discretization
Direct Numerical Simulations of turbulent flows must in principle resolve all flow scales in order to produce physically accurate results. In our case the largest eddies span the whole channel, whereas the smallest eddies are set by the Kolmogorov scale . In practice a grid spacing of a few is required Pope (2000); Moin et al. (1998), which results in very large numbers of grid points even at moderate Reynolds numbers. In addition, a statistical characterisation of the computed flow regimes implies that DNS need to span long time scales, in wallbounded shear flows several advective time units (). Hence DNS are computationally very expensive: the required computing power scales as under the assumption that algorithms scale optimally with problem size Pope (2000); Moin et al. (1998). This makes DNS prohibitive at large and requires accurate numerical methods and their efficient implementation even for moderate as in LAS precipitation.
3.1 Spatial discretization
We discretize the Navier–Stokes equations (1)–(2) using a order central finitevolume scheme in a collocated grid arrangement based on the work of Peric et al. Peric et al. (1988). For a detailed introduction to the finite volume method we refer the reader to the classic books of Ferziger and Peric (2012) and Versteeg and Malalasekera (2007). The code used here (FASTEST3D) has been developed at the Institute of Fluid Mechanics at the FriedrichAlexanderUniversität ErlangenNürnberg. It can solve for laminar and turbulent flows in complex geometries with DNS, LES and RANS. The code can be used also to solve fluidstructure Breuer et al. (2012) and fluidacoustic interaction Schäfer et al. (2010), as well as chemical reactions Enger et al. (2012). FASTEST3D is parallelized with a hybrid OpenMPMPI implementation Scheit et al. (2013) and has been optimized for both vector computers with high workload per process and for massively parallel systems with relatively low workload per process but great communication overhead.
The spatial discretization of the convection equation (3) is more challenging because of the large Schmidt number in our problem (). Let us consider a WW mixture with a mean velocity of , corresponding to in the inlet channels. At the flow is stationary but in the engulfment regime, which is characterized by two counterrotating vortices originating at the junction due to a KevinHelmholtz instability Fani et al. (2013). Fig. 2 shows that these vortices fold and stretch the contact plane between the two incoming fluid streams, resulting in sharp gradients in the concentration field . The discretization of sharp solutions by central discretization schemes lead to wiggles, overshoots and numerical oscillations whereas pure upwind discretization schemes promote mixing by introducing numerical diffusion Godunov (1959). Note that yields , and in order to get smooth numerical approximations with a second order central scheme prohibitive grid spacings would be required. To alleviate this problem the total variation diminishing (TVD) criterion was developed to ensure monotonicity of the bounded variables in time and space Harten (1983); Yee (1987); Waterson and Deconinck (2007). Additionally, limiter functions are used with the help of locally calculated successive gradients, to add numerical diffusion only where sharp gradients occur and to preserve boundedness. The fundamental finite volume discretization of those can be found in Versteeg and Malalasekera (2007).
Here we do not intend to present a new mathematical approach but rather show that TVD schemes perform remarkably well and are compulsory to maintain physical validity. We performed a grid convergence study for the WW mixture at in the engulfment regime shown. The relative error of the concentration field
(8) 
shown in Fig. 3(a), reflects an almost order accuracy. It is worth noting that sharp gradients are at about 45 degrees with respect to the control volume surfaces of our structured cartesian grid, corresponding to a worst case scenario.
3.2 Temporal discretization
DNS are often performed using semiimplicit schemes for advancing the Navier–Stokes equations in time Moin et al. (1998). We here employ an explicit method. Let us consider a WW mixture for ( in the inlet channel), which corresponds to a chaotic flow state. In this case the flow field is well resolved with a grid spacing of . In order that explicit or semiimplicit temporal schemes be stable the Courant–Friedrichs–Lewy (CFL) condition must be satisfied, which requires in this case a timestep size of about for a CFL constant of . The molecular diffusion and viscous transport restrictions on stability for an explicit scheme are and . We conclude that for both fully explicit and semiimplicit schemes are stable provided that the CFL condition is satisfied and this remains so as the Reynolds number further increases. Hence we decided to perform our DNS with fully explicit scheme because of the much lower computing cost involved.
We chose a lowstorage RungeKutta scheme decoupling equations (1) and (2) from (3). The right hand sides of (1) and (3) depend only on the spatial discretization, and for simplicity, these will be referred to as and , respectively. The timestepping algorithm is given below. Here the subscripts , refer to the central and surrounding center points of the control volumes respectively, as commonly used for finite volume discretizations Peric et al. (1988); Versteeg and Malalasekera (2007), and is the volume of the cell:

A predictor step that solves the momentum equation with a 3step RungeKutta time scheme
(9) where is not divergence free.

Accounting for equation (2) we formulate a corrector step for the pressure, leading to a Poisson type equation
(10) 
Solving the linear system for

Eq. 1 is solved with a 3step RungeKutta time scheme and the new velocity field
(13) (14) (15) 
update viscosity and density with
3.3 Grid study
We determined the required spatial resolution by two criteria. Firstly, we enforce that the difference in the pressure drop (6) computed with the chosen grid and a coarser grid is less than 1%. Secondly, in the mixing zone, where strong vortical structures arise, the grid spacing should be in the order of a few Kolmogorov length scales . In fact, we find that for the turbulent flow states suffices. The settings for the simulations in this paper are summarized in Table 1.
mixture  , , in [mm]  #grid points  in [s]  

220–650  WW  , ,  3.5–1  30–50  
300650  WE  , ,  2–1  30–50  
1100  WW,WE  , ,  0.4  4060  
2000  WW,WE  , ,  0.08  4060 
3.4 Code performance
(a)  (b) 

Simulating turbulent flow at requires the numerical code to scale in a supercomputer. The performance of our code is demonstrated in Fig. 3 with a speedup calculation for two sample setups grid64 (7 Mio. control volumes) and grid128 (55 Mio. control volumes). The speedup is calculated by , where is the mean execution time for one time step and is the mean execution time for the smallest number of cores. Ten sample time steps for a fixed number of pressure corrections is used to calculate the mean execution time. Notably, the large setup loses only 5% performance up to 3000 processors. A plain MPIparallelization strategy is adopted for the implementation of the code. The scaling test has been realized on SuperMUC at the Leibniz Supercomputing Centre (LRZ) in Germany.
4 Benchmark
We have tested our code for a wide range of Reynolds numbers and different geometries, including classical test cases for DNS such as channel flow Moin and Kim (1982). In this section we present two test cases of mixing of miscible fluids in a Tjunction. We provide a quantitative comparison against published results and discuss the physics of mixing as a reference for the cases at higher Reynolds number presented in §5.
4.1 Mixing of waterwater in a timeperiodic regime
(a)  (b)  (c) 

(a)  (b) 

Thomas and Ameel (2010) and Fani et al. (2014) showed experimentally and numerically that as the Reynolds number increases the steady engulfment regime shown in Fig. 2 is superseded by an unsteady periodically oscillating regime. The intrinsic KevinHelmholtz vortex structures persist, but, in comparison to the steady state, the vortex core moves with time as shown in the two snapshots of Fig. 4(a). Note also that the vortices are more intense and convoluted, and hence the spatiallyaveraged mixing degree at and has more than doubled from to . The timeperiodic nature of this regime can be seen in Fig. 5(a), which shows timeseries of for the three crosssections in Fig. 4. The frequency of the oscillation , expressed with the dimensionless Strouhal number , is shown in Fig. 5(b) as a function of Reynolds number. Here the onset of unsteadiness occurs at and our simulations are in very good agreement with those of Fani et al. (2014) in spite of the different widthtoheight aspectratio of the inlets (0.75 in Fani et al. (2014) and 1 here), outlet (1.5 in Fani et al. (2014) and 2 here) and length of the outlet channel (25 in Fani et al. (2014) and 12.5 here). Surprisingly, the frequency of the natural flow pulsation is rather insensitive to the details of the geometry. Note that at there is a sharp change in the oscillation frequency, which was attributed in Thomas and Ameel (2010); Fani et al. (2014) to a transition to a different more symmetric timeperiodic flow regime. As this is also observed here, it appears that this is also a robust feature of WW flows at Tjunctions.
4.2 Waterethanol mixing in the steady engulfment regime
We validated our code for WE mixtures by performing simulations of the steady engulfment regime at and . The crosssectional profiles of the concentration shown in Fig. 6 are in excellent agreement with those of Orsi et al. (2013) (see their Figure 6c–d). Interestingly, the emergence of the vortices near the junction as seen in Fig. 6(a) is comparable to that of the WW case as shown in Fig. 2(a). Orsi et al. (2013) found that in comparison to WW, in WE the transition to the engulfment regime is delayed to larger and that the flow remains steady up to larger . This is because of the stabilizing effect resulting from the increased viscosity of the mixture (see Fig. 1b). Thus, a comparison of the flow features at (WW) and (WE) is appropriate. In particular the mixing degree at is found to be in the same order as for WW and for WE, respectively. Note however that in WE the nonlinear dependence of the viscosity on the concentration makes the vortices asymmetric. In particular, the different strength of the vortices at their origin () generates a completely asymmetric flow pattern along the mixing channel.
(a)  (b)  (c) 

5 Transition to turbulence and multiplicity of flow states
(a)  (b) 

We investigated the transition to turbulence of WE mixtures by initializing our simulations from the steady engulfment regime at and then stepwise changing to a prescribed value. We found that for the flow starts pulsating periodically after a long transient of about . The resulting timeseries of the pressure loss necessary to drive the flow at the left (x) and right (+x) inlet are shown in Figure 7(a) for . As ethanol is more viscous than water, a higher pressure difference is necessary to drive ethanol at the same volumetric flow rate as water. Figure 7(b) shows the time evolution of the spatiallyaveraged mixing degree at three selected crosssections. As in the WW case shown in Fig. 5, a quiescent phase with relatively poor mixing is followed by a more dynamical phase characterized by a sharp peak in the mixing degree. By averaging these time series over a period we find that the mixing degree at is more than twice that of the steady engulfment regime at . Overall, our results are in qualitative agreement with the experiments of Wang et al. (2012), who reported a strong increase in the mixing efficiency at for WE mixtures. Note however that their experiments were carried out in a Tjunction with the same dimensions for the inlet and outlet channels, which prevents a direct comparison of our simulations and their experiments.
(a)  (b)  (c) 

Overall the WE timeperiodic engulfment regime is qualitatively similar to that in WW reported by Fani et al. (2013) and described in §4 for our geometry. However, the asymmetry in the viscosity as a function of concentration leads to different shear stresses at the contact plane between the two fluids and hence the arising KelvinHelmholtz vortices shown in Fig. 8(a) are slightly asymmetric. The growth and the intense convolution of the vortices substantially enhances mixing along the main channel, as can be seen in the snapshots of Fig. 8(b)–(c). This increase in mixing carries an increase of the mixture’s viscosity and hence towards the outlet vortices are damped and the flow starts to relaminarize. This process can be seen in the timeaveraged streamwise velocity at the same cross sections (bottom row of Fig. 8). Near the junction () an almost rotationally symmetric flow profile, which is a landmark of the engulfment state, is observed and progressing towards the outlet the irregular pattern of slowly evolves into a quasiparabolic flow profile at .
(a)  (b)  (c) 

(a)  (b) 

As the Reynolds number is further increased the flow becomes first quasiperiodic and finally chaotic. The top panel of Fig. 9(a) shows that in average an approximate rotational symmetry is retained near the junction, and hence we term this as the chaotic engulfment regime. The onset of chaotic mixing continues to enhances the performance of the Tmixer, until at the mixing degree suddenly drops by a factor of three (see Fig. 10a). This dramatic change in the mixing behavior is caused by the transition to a new topdown symmetric timeperiodic flow state (see the bottom row of Fig. 9). We investigated this surprising phenomenon by subsequently reducing the Reynolds number and found a hysteresis region in which the chaotic engulfment regime and the symmetric regime coexist. Note in fact the two flow states shown in Fig. 9 were both obtained at . Interestingly, Figure 10(b) shows that the pressure drop to drive the flow at a constant mass flux at a given number is higher for the engulfment regime than for the symmetric regime.
In the symmetric regime the water and ethanol streams flow in a nearly parallel fashion, thus strongly hindering mixing. Although there are vortices deforming the interface between the fluids locally, the intertwining of the two streams characteristic of the engulfment regime is entirely missing. By subsequently increasing , the symmetric flow suffers a transition to turbulence. As a result the mixing quality starts to slowly recover, however, it is only at that the same mixing degree as at in the engulfment regime is recovered. We note that the transition to turbulence and evolution of the mixing as a function of Reynolds number are similar for WW (see Fig. 10a) and hence not discussed here in detail.
6 Turbulent regime
(a) WE  (b) WW 

(a)  (b)  (c) 

Fig. 11 shows the time evolution of the crosssectionally averaged mixing degree at in the turbulent regime. Here the flow at the inlets is still assumed laminar because flows in straight square ducts do not sustain a turbulent state for Uhlmann et al. (2007). WW and WE have nearly the same mixing efficiency, but in WW the mixing degree fluctuates more strongly in time, as illustrated by the error bars of Fig. 10(a). The main flow features of WE are shown in Fig. 12. Near the junction the water and ethanol streams collide close to the center plane and the fluids are pressed towards the bottom and top wall as shown in the transversal velocities pointing in the direction. There, small irregular vortices arise and are transported along the mixing channel as they grow. Large fluctuations in the time evolution of show the raise of the ensuing large scale turbulent motions, whose strong transversal velocities promote mixing as depicted in the top and middle panels of Fig. 12(b)–(c). For WE the local increase of viscosity damps these fluctuations earlier than for WW. Despite the spatiotemporal chaotic nature of the flow, some poorly mixed patches of fluid reach the outlet.
At the flow is turbulent and hence the snapshots of the concentration in the top and middle rows of Fig. 12 do not display symmetric flow structures. However, the topdown symmetry is preserved in average as shown in the bottom row of Fig. 12. Note that the average is not leftright symmetric because of the different shear stress at the two sides of the contact plane between ethanol and water. The average pattern exhibits a pair of two counterrotating vortices, which are in fact quite similar to those of the steady vortex regime at low Reynolds number Fani et al. (2013); Thomas and Ameel (2010); Orsi et al. (2013), see the top panels of Fig. 2. At , the mean vortical structures are streaked with smallscale motions, while, close to the outlet at , the flow has nearly recovered to a quasiparabolic profile. The flow features of WW are very similar and hence not shown here. The main difference is that the timeaveraged velocity profile is leftright symmetric because in contrast to WE, the viscosity has no dependence on the concentration.
7 Statistical analysis of mixing properties
(a)  (b)  

(a)  (b)  

In the previous sections the mixing efficiency was evaluated with the mixing degree defined in equation (4), which entails a spatialaveraging over a channel crosssection. Hence cannot be used to infer the spatial homogeneity of the mixture, which is important for chemical reactions and precipitation processes Schwertfirm et al. (2007); Gradl et al. (2006); Wang et al. (2012). For the example of LAS precipitation, the driving force is supersaturation and this can vary by orders of magnitude depending on the local value of the concentration. In this section we briefly analyze the spatial homogeneity of the WE mixing process for the chaotic engulfment regime at and the (symmetric) turbulent regime at . Note that both cases exhibit a similar mixing degree at the outlet .
Fig. 13(a) shows probability distribution functions (PDF) of the concentration at several crosssections along the main channel for at . These PDFs were constructed by sampling a snapshot of a cross section every and then taking their time average. Fig. 13(b) shows two snapshots of the concentration at further illustrating the spatial evolution of the mixing process along the main channel. Between and the water and ethanol streams swap sides in the mixing channel, which is induced by the the growth of the KelvinHelmholtz vortices arising at the junction (see the top panel of Fig. 9a). This is a salient feature of the engulfment regime and accounts for its higher mixing degree, when compared to the symmetric regime at the same Reynolds number (see the lower panel of Fig. 9a). The intertwining of the fluid streams close to the junction result in a swift change from a bimodal to a plateaulike PDF. This can be interpreted following Schwertfirm et al. (2007), who found that the streamwise evolution of the PDF allows to determine the transport mechanisms dominating in the flow. In particular, a fast transition from a bimodal to plateau shape of the PDF, as seen in Fig. 13(a), implies that convective mixing is the dominant process, whereas an alteration in the plateau shape around points at mixing induced by molecular diffusion. Convective mixing is here ascribed to the stretching and folding of the contact plane between the fluids by vortices, which results in an increase of the flux of molecular diffusion. Toward the outlet the flow progressively relaminarises and molecular diffusion becomes the dominant process. In spite of the relatively high mixing degree at the outlet (), there are many parcels of fluid leaving the channel practically unmixed, as illustrated by the dark blue and red regions close to the outlet in Fig. 13(b).Thus it appears that the vortex roll up in the engulfment regime results in rapid mixing but mainly around the vortex core. The vortices arising at the junction strongly rotate as the grow along the mixing channel but fail to involve the total mass of fluid because of their gradual viscous decay in the streamwise direction.
It is useful to compare these results to the situation at , for which the mean flow structures are topdown symmetric and nearly leftright symmetric, as shown in Fig. 12. Despite the turbulent nature of the flow, Figure 14(b) shows that the incoming water and ethanol streams flow nearly parallel to each other. Here convective mixing occurs solely because of fluctuating vortices acting across the interface, whereas the mean vortical pattern does barely contribute to mixing. This is in sharp contrast to the engulfment regime, for which most of the convective mixing is induced by the largescale mean vortices swapping the streams immediately after the junction. This qualitative difference in the mixing process can be seen in the PDFs of Fig. 14(a), which remain largely bimodal. The peaks on the water () and ethanol () sides shrink steadily while moving progressively toward a completely mixes state (). Beyond the fluctuating vortices lose gradually their intensity, albeit they are stronger than at because of the lesser effect of viscosity, and there is a transition toward diffusion dominated mixing.
8 Conclusion
Efficiently mixing fluids is key to many applications in chemical engineering. Hence it is desirable to develop simulation tools for detailed and accurate predictions of mixing processes, for example of their spatiotemporal homogeneity. In the case of two different miscible fluids, such as water and ethanol in LAS precipitation, the viscosity and density of the mixture depend strongly on the concentration and so the mixing process is particularly challenging to simulate. We designed and implemented an accurate numerical scheme to compute temporally and spatially resolved velocity and concentration fields in the turbulent regime. This is necessary for many applications, as in LAS precipitation, where particle formation depends on the local concentration distribution and its evolution along particle trajectories. Our scheme is based on a fully explicit lowstorage RungeKutta method with decoupled transport equations and allows to simulate complex mixtures at the same computational cost as the transport of a passive scalar. Appropriate secondorder spatial discretizations with flux limiters were adopted to handle large Schmidt numbers, here , as usually occurring for miscible fluids. We demonstrated the potential of our methods by simulating waterwater and waterethanol mixtures for Reynolds numbers up to at a Tmixer.
At low Reynolds numbers the flow is stratified and the mixing is purely diffusive. Here the two fluid streams flow parallel to each other, despite the presence of vortices confined to each of the streams. For the waterwater case, this steady stratified vortex flow possess all the symmetries of the system, whereas in the waterethanol case these symmetries are only approximate because of the complex dependence of the viscosity and density on the composition. As increases the steady engulfment regime emerges. This is characterized by strong KelvinHelmholtz vortices at the junction, which break the leftright and updown symmetries of the stratified flow. Nevertheless the flow remains symmetric with respect to a 180 rotation (again for WE the symmetry is only approximate). This qualitative change in the vortical pattern results in the convolution of the two streams along the main channel, which greatly enlarges the contact surface between them and thus enhances mixing considerably. As further increases the flow becomes timeperiodic and later chaotic and hence mixing continues to improve. Note that in average the rotationally symmetric largescale vortices typical of the engulfment regime are preserved. Surprisingly, further increasing the Reynolds number results in a sudden reduction of the mixing efficiency, occurring for WW at and WE at . This is caused by a transition to a distinct flow regime, which despite being timeperiodic is topdown symmetric. In addition, near the junction it is nearly leftright symmetric and so the two fluid streams flow parallel to each other as in the stratified vortex flow at low Reynolds numbers. Here the mean vortical pattern contributes little to mixing and it is only through temporal fluctuations, i.e. vortices which act across the interface between the two fluids, that convective mixing occurs. As increases, the fluctuations become larger but it is only at that the same mixing degree as for the chaotic engulfment at is reached. However, it is worth noting that the mixing by temporal fluctuations at is spatially more homogenous than at , as in the latter many parcels of fluids leave the channel being practically unmixed.
An interesting feature of this transition between the two flow types is the existence of a hysteresis region. Hence depending on the initial conditions or mode of operation, one state or the other can be realized, which actually results in a factor of three different mixing degree. Although a transition to a new more symmetric flow state had been previously identified by DNS (Fani et al., 2014) as well as experiments (Thomas and Ameel, 2010) in the waterwater case, the finding of a hysteresis region and analysis of the mixing dynamics had not been reported. In agreement with previous simulations ((Orsi et al., 2013)), we found that the transition sequence occurs at larger for WE than for WW, which implies that at the same WW is better mixed. For the symmetric flow, however, this difference rapidly diminishes as increases, which can be attributed to the increasing importance of turbulent fluctuations and also to waterethanol mixture being only approximately symmetric in average.
Our numerical approach, and potential improvements such as the multipleresolution strategy of OstillaMonico et al. (2015), opens up the possibility to explore higher Reynolds number corresponding to real operating conditions with DNS, which will be pursued in the future. Our results contribute to understanding experimental measurements by providing detailed information of the velocity and concentration fields, and may be used as a reference to develop and test LES and RANS models. Finally, tracking the concentration and turbulent properties along flow trajectories will allow to deploy micromixing models and population balance equation solvers to predict the particle formation of LAS precipitation. This is part of our ongoing work.
The authors would like to acknowledge the funding of the Deutsche Forschungsgemeinschaft (DFG) through the Cluster of Excellence Engineering of Advanced Materials (EAM) and Bayer Technology Services GmbH (BTS). Computing time from the Regionales Rechenzentrum Erlangen (RRZE) is acknowledged.
Footnotes
 journal: Journal of Chemical Engineering Journal
References
 A. A. Thorat, S. V. Dalvi, Liquid antisolvent precipitation and stabilization of nanoparticles of poorly water soluble drugs in aqueous suspensions: Recent developments and future perspective, Chemical Engineering Journal 181â182 (2012) 1 – 34.
 M. Dizechi, E. Marschall, Viscosity of some binary and ternary liquid mixtures, Journal of Chemical and Engineering Data 27 (1982) 358–363.
 G. Orsi, M. Roudgar, E. Brunazzi, C. Galletti, R. Mauri, Waterâethanol mixing in tshaped microdevices, Chemical Engineering Science 95 (2013) 174 – 183.
 W. Wang, S. Zhao, T. Shao, Y. Jin, Y. Cheng, Visualization of microscale mixing in miscible liquids using Î¼lif technique and drug nanoparticle preparation in tshaped microchannels, Chemical Engineering Journal 192 (2012) 252 – 261.
 M. Hoffmann, M. SchlÃ¼ter, N. RÃ¤biger, Experimental investigation of liquidâliquid mixing in tshaped micromixers using lif and piv, Chemical Engineering Science 61 (2006) 2968 – 2976. Fluid Mixing {VIII} International ConferenceFluid Mixing {VIII} International Conference.
 D. Bothe, C. Stemich, H.J. Warnecke, Fluid mixing in a tshaped micromixer, Chemical Engineering Science 61 (2006) 2950 – 2958. Fluid Mixing {VIII} International ConferenceFluid Mixing {VIII} International Conference.
 D. Bothe, A. Lojewski, H.J. Warnecke, Fully resolved numerical simulation of reactive mixing in a tshaped micromixer using parabolized species equations, Chemical Engineering Science 66 (2011) 6424 – 6440. Novel Gas Conversion Symposium Lyon 2010, C1C4 Catalytic Processes for the Production of Chemicals and Fuels.
 A. Fani, S. Camarri, M. V. Salvetti, Investigation of the steady engulfment regime in a threedimensional tmixer, Physics of Fluids 25 (2013).
 A. Fani, S. Camarri, M. V. Salvetti, Unsteady asymmetric engulfment regime in a tmixer, Physics of Fluids 26 (2014).
 J. Gradl, H.C. Schwarzer, F. Schwertfirm, M. Manhart, W. Peukert, Precipitation of nanoparticles in a tmixer: Coupling the particle population dynamics with hydrodynamics through direct numerical simulation, Chemical Engineering and Processing: Process Intensification 45 (2006) 908 – 916. Particulate Processes A Special Issue of Chemical Engineering and Processing.
 F. Schwertfirm, J. Gradl, H. C. Schwarzer, W. Peukert, M. Manhart, The low reynolds number turbulent flow and mixing in a confined impinging jet reactor, International Journal of Heat and Fluid Flow 28 (2007) 1429 – 1442. Revised and extended papers from the 5th conference in Turbulence, Heat and Mass TransferâThe 5th conference in Turbulence, Heat and Mass Transfer.
 S. Thomas, T. A. Ameel, An experimental investigation of moderate reynolds number flow in a tchannel, Experiments in fluids 49 (2010) 1231–1245.
 D. Joseph, Fluid dynamics of two miscible liquids with diffusion and gradient stresses, European journal of mechanics. B, Fluids 9 (1990) 565–596.
 D. D. Joseph, A. Huang, H. Hu, Nonsolenoidal velocity effects and korteweg stresses in simple mixtures of incompressible liquids, Physica D: Nonlinear Phenomena 97 (1996) 104 – 125.
 J. C. Legros, Y. Gaponenko, A. Mialdun, T. Triller, A. Hammon, C. Bauer, W. Köhler, V. Shevtsova, Investigation of fickian diffusion in the ternary mixtures of water–ethanol–triethylene glycol and its binary pairs, Physical Chemistry Chemical Physics 17 (2015) 27713–27725.
 S. Pope, Turbulent Flows, Cambridge University Press, 2000. URL: https://books.google.de/books?id=HZsTw9SMx0C.
 P. Moin, , K. Mahesh, Direct numerical simulation: A tool in turbulence research, Annual Review of Fluid Mechanics 30 (1998) 539–578.
 M. Peric, R. Kessler, G. Scheuerer, Comparison of finitevolume numerical methods with staggered and colocated grids, Computers and Fluids 16 (1988) 389 – 403.
 J. H. Ferziger, M. Peric, Computational methods for fluid dynamics, Springer Science & Business Media, 2012.
 H. K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, Pearson Education, 2007.
 M. Breuer, G. De Nayer, M. Münsch, T. Gallinger, R. Wüchner, Fluid–structure interaction using a partitioned semiimplicit predictor–corrector coupling scheme for the application of largeeddy simulation, Journal of Fluids and Structures 29 (2012) 107–130.
 F. Schäfer, S. Müller, T. Uffinger, S. Becker, J. Grabinger, M. Kaltenbacher, Fluidstructureacoustic interaction of the flow past a thin flexible structure, AIAA journal 48 (2010) 738–748.
 S. Enger, F. Schäfer, M. Breuer, F. Durst, Czochralski melt flow on highperformance, in: High Performance Scientific And Engineering Computing: Proceedings of the 3rd International FORTWIHR Conference on HPSEC, Erlangen, March 12–14, 2001, volume 21, Springer Science & Business Media, 2012, p. 201.
 C. Scheit, G. Hager, J. Treibig, S. Becker, G. Wellein, Optimization of FASTEST3D for modern multicore systems, CoRR abs/1303.4538 (2013).
 S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Matematicheskii Sbornik 89 (1959) 271–306.
 A. Harten, High resolution schemes for hyperbolic conservation laws, Journal of computational physics 49 (1983) 357–393.
 H. Yee, Construction of explicit and implicit symmetric tvd schemes and their applications, Journal of Computational Physics 68 (1987) 151–179.
 N. Waterson, H. Deconinck, Design principles for bounded higherorder convection schemes â a unified approach, Journal of Computational Physics 224 (2007) 182 – 207. Special Issue Dedicated to Professor Piet Wesseling on the occasion of his retirement from Delft University of Technology.
 P. Moin, J. Kim, Numerical investigation of turbulent channel flow, Journal of fluid mechanics 118 (1982) 341–377.
 M. Uhlmann, A. Pinelli, G. Kawahara, A. Sekimoto, Marginally turbulent flow in a square duct., Journal of fluid mechanics 588 (2007) 153–162.
 R. OstillaMonico, Y. Yang, E. van der Poel, D. Lohse, R. Verzicco, A multipleresolution strategy for direct numerical simulation of scalar turbulence, Journal of Computational Physics 301 (2015) 308 – 321.