Direct numerical simulation of water-ethanol flows in a T-mixer

Direct numerical simulation of water-ethanol flows in a T-mixer


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 water-water mixtures in simple geometries such as T-mixers 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 two-fluid mixtures. Using a high-performance implementation of this method we perform direct numerical simulations resolving the spatial and temporal dynamics of water-ethanol 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 water-water mixtures.

Direct numerical simulation, miscible fluids, water-ethanol, turbulent flows, T-mixer, mixing, transition to turbulence

1 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 top-down or bottom-up and usually deal with the production of nanocrystals containing mainly active pharmaceutical ingredient. Liquid antisolvent (LAS) precipitation is a promising, yet poorly understood, bottom-up 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 up-scaling, 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 anti-solvent 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)
Figure 1: (a) T-shaped mixer consisting of two quadratic inlet channels and a rectangular mixing outlet channel. The width and height of the inlet channel is mm and the origin of the coordinate system is located at the middle of the mixing and inlet channel. The dark-shaded cross section at is used later to characterise the properties of the mixing process. (b) Kinematic viscosity (y-axis) and density (alt y-axis) as a function of the concentration for a water-ethanol mixtures. Here and correspond to the pure water and ethanol phases, respectively. The solid lines indicate functions fitted to the experimental data (dots and squares), see Dizechi and Marschall (1982).

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 anti-solvent (water), see Orsi et al. (2013); Wang et al. (2012). Water-water (WW) mixtures in laminar and slightly transient regimes in T-shaped 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 T-micromixer 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 water-water system, the velocity field and concentration field are top-down and left-right 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 so-called stratified regime the vortices are top-down and left-right 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 left-right and top-down 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)
Figure 2: Colormap of the concentration at three cross-sections along the main channel for WW in the steady vortex regime ( and , top row) and steady engulfment regime ( and , bottom row). The arrows show the cross-sectional velocity field. At , the spatially-averaged mixing degree (see eq. 4) is ), respectively, while at the mixing degree is .

The first investigation of the mixing of water and ethanol (WE) was performed by Wang et al. (2012), who conducted -LIF measurements in a T-mixer. 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 water-water 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 convection-diffusion equation, which in conservative form reads


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


which is known as the quasi-incompressibility 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


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 water-ethanol mixtures depends on the concentration Legros et al. (2015), and this dependence should be considered in sub-grid 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 concentration-dependent 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 two-dimensional 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 no-slip boundary condition at the wall is employed for the velocity field, whereas a zero-flux 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


where is the maximum of the variance achieved by completely segregated streams, and is the mean deviation of the concentration profile as


where is the mean concentration of a cross section. The energy input is here quantified with the pressure loss


where and are the surface-averaged pressure at the outlet and inlet cross section, respectively. Using the pressure loss the Darcy friction factor is calculated as


where is the total length of the T-mixer 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 wall-bounded 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 finite-volume 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 (FASTEST-3D) has been developed at the Institute of Fluid Mechanics at the Friedrich-Alexander-Universität Erlangen-Nü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 fluid-structure Breuer et al. (2012) and fluid-acoustic interaction Schäfer et al. (2010), as well as chemical reactions Enger et al. (2012). FASTEST-3D is parallelized with a hybrid OpenMP-MPI 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 counter-rotating vortices originating at the junction due to a Kevin-Helmholtz 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


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 semi-implicit 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 semi-implicit temporal schemes be stable the Courant–Friedrichs–Lewy (CFL) condition must be satisfied, which requires in this case a time-step 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 semi-implicit 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 low-storage Runge-Kutta 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 time-stepping 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:

  1. A predictor step that solves the momentum equation with a 3-step Runge-Kutta time scheme


    where is not divergence free.

  2. Accounting for equation (2) we formulate a corrector step for the pressure, leading to a Poisson type equation

  3. Solving the linear system for

  4. We obtain a new pressure field which complies the enforced Eq. 2 by updating the velocity field


    and the pressure field


    The corrector steps 2-4 are repeated until Eq. 2 is fulfilled within machine precision and .

    • In the first correction loop

    • If the stop criterion is not fulfilled and

  5. Eq. 1 is solved with a 3-step Runge-Kutta time scheme and the new velocity field

  6. 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
300-650 WE , , 2–1 30–50
1100 WW,WE , , 0.4 40-60
2000 WW,WE , , 0.08 40-60
Table 1: Overview of simulation settings. Grid spacings are given for the mixing channel. In the junction an equidistant grid spacing is chosen having the finer shown resolution. Time steps are selected such that more than 80% of the simulation time the CFL condition is satisfied with everywhere in the domain. Depending on the branch and initial guess it takes up to 60 advective time units to reach statistical steady state (last column).

3.4 Code performance

(a) (b)
Figure 3: (a) Relative error of the concentration field (8) as function of the grid spacing . The TVD scheme with a van Leer limiter function exhibits spatial convergence of approximately seocod order . (b) Speed-up (strong scaling) over cores for two different grids sizes. Lines show perfect scaling for the corresponding setup.

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 MPI-parallelization 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 T-junction. 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 water-water in a time-periodic regime

(a) (b) (c)
Figure 4: Snapshots of the colormap of the concentration for WW in the time-periodic engulfment regime () at (top row) and (bottom row) corresponding to the two different phases of mixing as it can be seen in Fig. 5. The arrows show the cross-sectional velocity field.
(a) (b)
Figure 5: (a) Time-evolution of the mixing degree at three selected cross sections for in the time-periodic engulfment regime. The time-averaged mixing degree is , respectively. (b) Strouhal number as a function of Reynolds number (red circles, the line is to guide the eyes). The data from the simulations of Fani et al. (2014) is shown as black triangles. Because of the slightly different geometry their data have been scaled here with the inlet hydraulic diameter to allow for comparison.

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 Kevin-Helmholtz 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 spatially-averaged mixing degree at and has more than doubled from to . The time-periodic nature of this regime can be seen in Fig. 5(a), which shows time-series of for the three cross-sections 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 width-to-height aspect-ratio 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 time-periodic flow regime. As this is also observed here, it appears that this is also a robust feature of WW flows at T-junctions.

4.2 Water-ethanol mixing in the steady engulfment regime

We validated our code for WE mixtures by performing simulations of the steady engulfment regime at and . The cross-sectional 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)
Figure 6: Colormaps of the concentration for WE in the steady engulfment regime at (top row) and (bottom row). The arrows show the cross-sectional velocity field.

5 Transition to turbulence and multiplicity of flow states

(a) (b)
Figure 7: Time-periodic engulfment regime for WE at . (a) Time evolution of surface averaged pressure loss computed for water (; ) and ethanol (; ) inlets. (b) Time evolution of the mixing degree at three selected cross sections (see legend). The time-averaged mixing degree is , respectively.

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 time-series 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 spatially-averaged mixing degree at three selected cross-sections. 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 T-junction with the same dimensions for the inlet and outlet channels, which prevents a direct comparison of our simulations and their experiments.

(a) (b) (c)
Figure 8: Time-periodic engulfment regime for WE at . The top and middle row show snapshots of the concentration (colormap) and cross-sectional velocity field (arrows) at and , respectively. The bottom row shows the time-averaged streamwise (colormap) and cross-sectional velocities (arrows).

Overall the WE time-periodic 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 Kelvin-Helmholtz 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 time-averaged 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 quasi-parabolic flow profile at .

(a) (b) (c)
Figure 9: Colormap of the time-averaged streamwise velocity for WE at in the chaotic engulfment regime (top row) and time-periodic symmetric regime (bottom row). The arrows show the time-averaged cross-sectional velocity field.
(a) (b)
Figure 10: (a) Mixing degree evaluated at cross section as a function of Reynolds number. (b) Darcy friction factor as a function of Reynolds number in a log-log scale. Error bars illustrate the maximum and minimum values of and for an unsteady flow, respectively.

As the Reynolds number is further increased the flow becomes first quasi-periodic 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 T-mixer, 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 top-down symmetric time-periodic 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
Figure 11: Time evolution of mixing degree at for WE (a) and WW (b) at .
(a) (b) (c)
Figure 12: Turbulent flow of WE at . The top and middle row show snapshots of the concentration (colormap) and cross-sectional velocity field (arrows) at and , respectively. The bottom row shows the time-averaged streamwise (colormap) and cross-sectional (arrows) velocities.

Fig. 11 shows the time evolution of the cross-sectionally 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 spatio-temporal 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 top-down symmetry is preserved in average as shown in the bottom row of Fig. 12. Note that the average is not left-right 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 counter-rotating 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 small-scale motions, while, close to the outlet at , the flow has nearly recovered to a quasi-parabolic profile. The flow features of WW are very similar and hence not shown here. The main difference is that the time-averaged velocity profile is left-right symmetric because in contrast to WE, the viscosity has no dependence on the concentration.

7 Statistical analysis of mixing properties

(a) (b)
Figure 13: Stream-wise evolution of the mixing for in the chaotic engulfment regime at . (a) Time-averaged probability density function of the concentration at selected cross sections in the mixing channel. (b) Snapshots of the concentration at mid-height () at time (left) and (right).
(a) (b)
Figure 14: (a) Time-averaged probability density function of concentration for WE at and selected cross sections in the mixing channel. (b) Snapshots of colormap of concentration of WE for cross section at time (left) and (right) at corresponding to the time in Fig. 12.

In the previous sections the mixing efficiency was evaluated with the mixing degree defined in equation (4), which entails a spatial-averaging over a channel cross-section. 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 cross-sections 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 Kelvin-Helmholtz 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 plateau-like 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 top-down symmetric and nearly left-right 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 large-scale 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 spatio-temporal 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 low-storage Runge-Kutta method with decoupled transport equations and allows to simulate complex mixtures at the same computational cost as the transport of a passive scalar. Appropriate second-order 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 water-water and water-ethanol mixtures for Reynolds numbers up to at a T-mixer.

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 water-water case, this steady stratified vortex flow possess all the symmetries of the system, whereas in the water-ethanol 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 Kelvin-Helmholtz vortices at the junction, which break the left-right and up-down 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 time-periodic and later chaotic and hence mixing continues to improve. Note that in average the rotationally symmetric large-scale 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 time-periodic is top-down symmetric. In addition, near the junction it is nearly left-right 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 water-water 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 water-ethanol mixture being only approximately symmetric in average.

Our numerical approach, and potential improvements such as the multiple-resolution strategy of Ostilla-Monico 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.


  1. journal: Journal of Chemical Engineering Journal


  1. 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.
  2. M. Dizechi, E. Marschall, Viscosity of some binary and ternary liquid mixtures, Journal of Chemical and Engineering Data 27 (1982) 358–363.
  3. G. Orsi, M. Roudgar, E. Brunazzi, C. Galletti, R. Mauri, Water–ethanol mixing in t-shaped microdevices, Chemical Engineering Science 95 (2013) 174 – 183.
  4. W. Wang, S. Zhao, T. Shao, Y. Jin, Y. Cheng, Visualization of micro-scale mixing in miscible liquids using μ-lif technique and drug nano-particle preparation in t-shaped micro-channels, Chemical Engineering Journal 192 (2012) 252 – 261.
  5. M. Hoffmann, M. Schlüter, N. Räbiger, Experimental investigation of liquid–liquid mixing in t-shaped micro-mixers using -lif and -piv, Chemical Engineering Science 61 (2006) 2968 – 2976. Fluid Mixing {VIII} International ConferenceFluid Mixing {VIII} International Conference.
  6. D. Bothe, C. Stemich, H.-J. Warnecke, Fluid mixing in a t-shaped micro-mixer, Chemical Engineering Science 61 (2006) 2950 – 2958. Fluid Mixing {VIII} International ConferenceFluid Mixing {VIII} International Conference.
  7. D. Bothe, A. Lojewski, H.-J. Warnecke, Fully resolved numerical simulation of reactive mixing in a t-shaped micromixer using parabolized species equations, Chemical Engineering Science 66 (2011) 6424 – 6440. Novel Gas Conversion Symposium- Lyon 2010, C1-C4 Catalytic Processes for the Production of Chemicals and Fuels.
  8. A. Fani, S. Camarri, M. V. Salvetti, Investigation of the steady engulfment regime in a three-dimensional t-mixer, Physics of Fluids 25 (2013).
  9. A. Fani, S. Camarri, M. V. Salvetti, Unsteady asymmetric engulfment regime in a t-mixer, Physics of Fluids 26 (2014).
  10. J. Gradl, H.-C. Schwarzer, F. Schwertfirm, M. Manhart, W. Peukert, Precipitation of nanoparticles in a t-mixer: 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.
  11. 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.
  12. S. Thomas, T. A. Ameel, An experimental investigation of moderate reynolds number flow in a t-channel, Experiments in fluids 49 (2010) 1231–1245.
  13. D. Joseph, Fluid dynamics of two miscible liquids with diffusion and gradient stresses, European journal of mechanics. B, Fluids 9 (1990) 565–596.
  14. D. D. Joseph, A. Huang, H. Hu, Non-solenoidal velocity effects and korteweg stresses in simple mixtures of incompressible liquids, Physica D: Nonlinear Phenomena 97 (1996) 104 – 125.
  15. 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.
  16. S. Pope, Turbulent Flows, Cambridge University Press, 2000. URL:
  17. P. Moin, , K. Mahesh, Direct numerical simulation: A tool in turbulence research, Annual Review of Fluid Mechanics 30 (1998) 539–578.
  18. M. Peric, R. Kessler, G. Scheuerer, Comparison of finite-volume numerical methods with staggered and colocated grids, Computers and Fluids 16 (1988) 389 – 403.
  19. J. H. Ferziger, M. Peric, Computational methods for fluid dynamics, Springer Science & Business Media, 2012.
  20. H. K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, Pearson Education, 2007.
  21. M. Breuer, G. De Nayer, M. Münsch, T. Gallinger, R. Wüchner, Fluid–structure interaction using a partitioned semi-implicit predictor–corrector coupling scheme for the application of large-eddy simulation, Journal of Fluids and Structures 29 (2012) 107–130.
  22. F. Schäfer, S. Müller, T. Uffinger, S. Becker, J. Grabinger, M. Kaltenbacher, Fluid-structure-acoustic interaction of the flow past a thin flexible structure, AIAA journal 48 (2010) 738–748.
  23. S. Enger, F. Schäfer, M. Breuer, F. Durst, Czochralski melt flow on high-performance, 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.
  24. C. Scheit, G. Hager, J. Treibig, S. Becker, G. Wellein, Optimization of FASTEST-3D for modern multicore systems, CoRR abs/1303.4538 (2013).
  25. S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Matematicheskii Sbornik 89 (1959) 271–306.
  26. A. Harten, High resolution schemes for hyperbolic conservation laws, Journal of computational physics 49 (1983) 357–393.
  27. H. Yee, Construction of explicit and implicit symmetric tvd schemes and their applications, Journal of Computational Physics 68 (1987) 151–179.
  28. N. Waterson, H. Deconinck, Design principles for bounded higher-order 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.
  29. P. Moin, J. Kim, Numerical investigation of turbulent channel flow, Journal of fluid mechanics 118 (1982) 341–377.
  30. M. Uhlmann, A. Pinelli, G. Kawahara, A. Sekimoto, Marginally turbulent flow in a square duct., Journal of fluid mechanics 588 (2007) 153–162.
  31. R. Ostilla-Monico, Y. Yang, E. van der Poel, D. Lohse, R. Verzicco, A multiple-resolution strategy for direct numerical simulation of scalar turbulence, Journal of Computational Physics 301 (2015) 308 – 321.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description