1 Introduction
Mixing layers formed between parallel gas and liquid streams are commonly seen in nature and industrial applications, e.g., breaking ocean waves and injection of liquid fuels into engines. Typically a velocity difference exists between the two streams, which triggers a shear instability on the gas–liquid interface. The interfacial instability grows and eventually causes the bulk liquid to break into small droplets, forming a two-phase mixing layer between the two streams. At the continuum level, the gas and liquid streams are immiscible, so the ‘mixing’ layer here indeed refers to a layer consisting of a mixture of gas and a dispersion of droplets generated from the bulk liquid disintegration. The process where the bulk liquid stream breaks into a large number of small droplets is often referred to as ‘atomization’ and the resulting gas–droplets mixture as a spray. Since the breakup of the liquid stream can be significantly enhanced by the parallel fast gas stream, this co-flowing configuration (also known as air-blast atomization) is widely used in fuel injectors (Lefebvre & McDonell Reference Lefebvre and McDonell2017).
1.1 Problem description
In the present study we focus on modelling the wall-bounded two-phase mixing layer experiment by Matas, Marty & Cartellier (Reference Matas, Marty and Cartellier2011), which is illustrated in figure 1 by a snapshot of present simulation results (details of simulation are to be presented later). The grey surface is the liquid–gas interface and the background is the
$z$
-component (in spanwise direction) of vorticity. The parallel gas and liquid streams, separated by a small separator plate, enter the domain from the left. The thicknesses of the two streams at the inlet are the same. A mixing layer is formed between the gas stream and the stagnant gas, as indicated by purple dashed lines. Similarly, a two-phase mixing layer is formed between the gas and liquid streams, between the two green dashed lines. The two-phase mixing layer near the inlet is simply a downstream propagating interfacial wave, strictly speaking there is no ‘mixing’. The mixing between the two immiscible phases only occur further downstream when the liquid stream breaks up, forming a mixture of gas and droplets.
A two-phase mixing layer between parallel gas and liquid streams. The grey surface is the liquid–gas interface and the background is the
$z$
-component of vorticity.
As indicated in this figure, the two-phase mixing layer is a phenomenon of enormous complexity that involves interfacial instability, two-phase turbulence and topological changes owing to liquid breakups occurring at a wide variety of spatial and temporal scales. Owing to its important application to fuel injection, the problem has attracted increasing attention in recent years and extensive theoretical, experimental and numerical investigations have be performed. Some of the important previous works in the literature are discussed in the next sections.
1.2 Shear-induced interfacial instability
The destabilization of the liquid stream in the present problem is initiated by the instability at the gas–liquid interface, which is in turn induced by the velocity difference between parallel streams and the resulting shear on the interface near the nozzle (Taylor Reference Taylor and Batchelor1963; Renardy Reference Renardy1985; Rangel & Sirignano Reference Rangel and Sirignano1988; Lasheras, Villermaux & Hopfinger Reference Lasheras, Villermaux and Hopfinger1998; Lasheras & Hopfinger Reference Lasheras and Hopfinger2000; Matas et al. Reference Matas, Marty and Cartellier2011; Matas Reference Matas2015). Different mechanisms that drive the interfacial instability, including the classic Kelvin–Helmholtz instability (Helmholtz Reference Helmholtz1868; Thomson Reference Thomson1871), the instability owing to viscosity contrast (Yih Reference Yih1967), the inviscid Rayleigh (inflection-point) and viscous Tollmien–Schlichting mechanisms, have been addressed in previous works (Ozgen, Degrez & Sarma Reference Ozgen, Degrez and Sarma1998; Otto, Rossi & Boeck Reference Otto, Rossi and Boeck2013). Stability analysis at the interface between two immiscible fluids of different densities and viscosities has been conducted for both planar (Renardy Reference Renardy1985; Rangel & Sirignano Reference Rangel and Sirignano1988, Reference Rangel and Sirignano1991; Matas et al. Reference Matas, Marty and Cartellier2011) and cylindrical geometries (Raynal Reference Raynal1997; Lasheras et al. Reference Lasheras, Villermaux and Hopfinger1998; Marmottant & Villermaux Reference Marmottant and Villermaux2004).
Conventionally, the development of the interfacial wave is investigated through a linear analysis of small perturbations of the two-dimensional base flow (with no variation in the transverse direction) and studies focus on predicting the most unstable wavelength and frequency. The linear instability studies that yield theoretical prediction of the most unstable frequency were first carried out assuming inviscid flows (Marmottant & Villermaux Reference Marmottant and Villermaux2004; Eggers & Villermaux Reference Eggers and Villermaux2008; Matas et al. Reference Matas, Marty and Cartellier2011), and extensions to viscous regime have been made in recent years (Boeck & Zaleski Reference Boeck and Zaleski2005; Sahu et al. Reference Sahu, Valluri, Spelt and Matar2007; Fuster et al. Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013; Otto et al. Reference Otto, Rossi and Boeck2013; O’Naraigh, Spelt & Shaw Reference O’Naraigh, Spelt and Shaw2013; O’Naraigh et al. Reference O’Naraigh, Valluri, Scott, Bethune and Spelt2014; Matas et al. Reference Matas, Marty, Dem and Cartellier2015). Research efforts have also been made to investigate the effect of confinement on the instability of mixing layers (Juniper & Candel Reference Juniper and Candel2003; Juniper Reference Juniper2006; Juniper, Tammisola & Lundell Reference Juniper, Tammisola and Lundell2011; Matas Reference Matas2015).
While the inviscid stability theory has been shown to well predict the scaling relation of the most unstable wavelength and frequencies (Rangel & Sirignano Reference Rangel and Sirignano1988, Reference Rangel and Sirignano1991; Raynal Reference Raynal1997; Marmottant & Villermaux Reference Marmottant and Villermaux2004; Eggers & Villermaux Reference Eggers and Villermaux2008), the viscous stability analysis is required in general to yield accurate prediction of the magnitudes of the most unstable frequency and wavelength (Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013; Otto et al.
Reference Otto, Rossi and Boeck2013; O’Naraigh et al.
Reference O’Naraigh, Spelt and Shaw2013; Matas et al.
Reference Matas, Marty, Dem and Cartellier2015). The linear stability analysis in both inviscid and viscous regime confirms that the vorticity layer thickness of the gas stream at the inlet, denoted by
$\unicode[STIX]{x1D6FF}_{g}$
, is the characteristic length scale that controls the selection of the most unstable wavelength.
It is clearly shown by previous works (Raynal Reference Raynal1997; Hoepffner, Blumenthal & Zaleski Reference Hoepffner, Blumenthal and Zaleski2011; Ling et al.
Reference Ling, Fuster, Zaleski and Tryggvason2017) that the propagation speed of the interfacial wave is well predicted by the Dimotakis speed. The Dimotakis speed
$U_{D}$
is defined as (Dimotakis Reference Dimotakis1986)
$$\begin{eqnarray}\displaystyle U_{D}=\frac{\sqrt{\unicode[STIX]{x1D70C}_{l}}U_{l}+\sqrt{\unicode[STIX]{x1D70C}_{g}}U_{g}}{\sqrt{\unicode[STIX]{x1D70C}_{l}}+\sqrt{\unicode[STIX]{x1D70C}_{g}}}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}_{l}$
and
$\unicode[STIX]{x1D70C}_{g}$
are the liquid and gas densities and the subscripts
$g$
and
$l$
represent the gas and liquid phases, respectively. The velocities of the liquid and gas stream at the inlet are denoted by
$U_{l}$
and
$U_{g}$
, and
$U_{D}$
is obtained through a phenomenological approach, assuming the gas and liquid dynamic pressures are in a balance in the reference frame moving with the wave speed. With the Dimotakis speed, the frequency and wavelength of the most unstable wave can then be related to each other as
$f=U_{D}/\unicode[STIX]{x1D706}$
.
Recently, through viscous spatial–temporal analysis, Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) and Otto et al. (Reference Otto, Rossi and Boeck2013) showed that the interfacial instability can be absolute or convective, depending mainly on the dynamic pressure ratios between the two phases,
$M$
, defined as
$$\begin{eqnarray}\displaystyle M=\frac{\unicode[STIX]{x1D70C}_{g}U_{g}^{2}}{\unicode[STIX]{x1D70C}_{l}U_{l}^{2}}. & & \displaystyle\end{eqnarray}$$
When
$M$
is large, the interfacial instability is absolute and the wave frequency predicted by stability analysis was found to agree well with experiments and simulations (Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013; Agbaglah, Chiodi & Desjardins Reference Agbaglah, Chiodi and Desjardins2017).
Studies of interfacial instability eventually aim to shed light on understanding the behaviour of liquid breakups occurring further downstream. Nevertheless, connections between the upstream interfacial instability and the downstream turbulence and spray characteristics have not been investigated thoroughly in previous studies. A possible reason is that a three-dimensional simulation that can resolve both the interfacial instability and the resulting turbulent spray is too expensive. In this study we only consider one specific case of two-phase mixing layer which clearly lies in the absolute instability regime. As a result, there is a dominant interfacial stability frequency and we will compare the numerical results with the spatial viscous stability theory of Otto et al. (Reference Otto, Rossi and Boeck2013).
1.3 Two-phase turbulent coherent structures
As the interfacial wave is formed, turbulent coherent structures simultaneously appear in the gas stream near the interface (Bernal & Roshko Reference Bernal and Roshko1986). Owing to the significant difference in velocities and viscosities between the gas and liquid streams, the gas–liquid interface acts like a deforming wavy ‘wall’ to the gas stream. The resulting turbulent vortical structures are similar to those in boundary layers (Wu & Moin Reference Wu and Moin2009; Jodai & Elsinga Reference Jodai and Elsinga2016). The growing interfacial waves significantly perturb the gas stream and play a significant role in the transition to turbulence. When the amplitude of the interfacial wave is large compared with the thickness of the gas stream, it appears as an obstacle to the gas flow, causing the latter to separate downstream of the wave. As discussed above, the frequency of wave formation will correspond to the fastest growing mode if the instability is absolute. Therefore, the turbulence production will also be related to the wave frequency. Not only can the interfacial wave development modulate the gas flow, the vortices in the gas flow also influence the wave evolution and the subsequent breakup (Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016). The liquid stream eventually disintegrates into a large number of droplets with a wide range of sizes. These droplets are dispersed in the turbulent flow. Secondary breakup or coalescence may also occur. The present study aims to provide rigorous statistics of multiphase turbulence in the mixing layer.
1.4 Direct numerical simulation of chaotic liquid breakups and topology changes
Thanks to the rapid development of computer power and numerical methodology, direct numerical simulations (DNS) of atomization have become viable in the past decade and recent simulations have provided high-resolution details of atomization, including interfacial instability development, interaction between the interfacial wave and the turbulent gas stream, and formation of liquid sheets, ligaments and droplets (Ménard, Tanguy & Berlemont Reference Ménard, Tanguy and Berlemont2007; Shinjo & Umemura Reference Shinjo and Umemura2010; Rana & Herrmann Reference Rana and Herrmann2011; Le Chenadec & Pitsch Reference Le Chenadec and Pitsch2013; Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016; Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017; Agbaglah et al. Reference Agbaglah, Chiodi and Desjardins2017). In particular, different droplets formation mechanisms have been observed. When the interfacial waves roll up and develop into liquid sheets, Taylor–Culick rims form at the edges of liquid sheets. Rayleigh–Taylor (RT) or Rayleigh–Plateau (RP) instabilities in the transverse direction then develop at the rims, generating liquid fingers and filaments (Marmottant & Villermaux Reference Marmottant and Villermaux2004; Roisman, Horvat & Tropea Reference Roisman, Horvat and Tropea2006; Agbaglah, Josserand & Zaleski Reference Agbaglah, Josserand and Zaleski2013). These filaments finally break into a distribution of small droplets. In addition to this well-known finger mechanism, simulation results also reveal a less-established mechanism, i.e., holes form in liquid sheets and the spontaneous expansion of these holes causes liquid sheets to rupture violently, producing numerous filaments and droplets of different sizes (Shinjo & Umemura Reference Shinjo and Umemura2010; Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016; Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017; Zandian, Sirignano & Hussain Reference Zandian, Sirignano and Hussain2017, Reference Zandian, Sirignano and Hussain2018). The hole-induced breakup of a thin liquid sheet is also observed in the bag breakup of a drop in secondary atomization (Opfer et al. Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014) and splashes (Marston et al. Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016). The mechanisms that cause sheet deformation and hole formation have been investigated recently via vortex dynamics by Zandian et al. (Reference Zandian, Sirignano and Hussain2018). In current interface-resolved simulations, disjoining pressure is generally ignored as the affordable minimum mesh size is still far larger than the sheet thickness where molecular forces are active in collapsing a liquid sheet. The holes observed in the simulations are thus an outcome of the numerical cut-off length scale, i.e., the cell size (typically in micrometres or sub-micrometres). Nevertheless, recent experiments in splash and secondary breakup interestingly show that holes indeed form in liquid sheets when the thickness is around micrometres (Opfer et al. Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014; Marston et al. Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016). The reasons for holes arising in a thicker sheet are not fully understood, but experiments seem to indicate that the holes observed in atomization simulations are not far from what is observed in reality.
1.5 Modelling of turbulent atomization
An important future direction of atomization simulations is the development of sub-grid models like large-eddy simulation (LES) for turbulent single-phase flow (Pope Reference Pope2000). Spatial scales involved in atomization processes, varying from the size of the injector to the diameter of the smallest droplet, can easily go beyond three or four orders of magnitudes. If one has to fully resolve all the scales to guarantee reasonably accurate macro-scale features, the impact of numerical simulations to practical atomization applications will be limited by their extreme costs.
Attempts to combining interface-capturing schemes and Lagrangian point-particle models have been proposed in recent years (Herrmann Reference Herrmann2010; Tomar et al. Reference Tomar, Fuster, Zaleski and Popinet2010; Ling, Zaleski & Scardovelli Reference Ling, Zaleski and Scardovelli2015; Zuzio, Estivalezes & DiPierro Reference Zuzio, Estivalezes and Dipierro2017). In these combined approaches, the interfaces of the small droplets are not resolved as for the macro-scale interfaces; instead, the droplets are treated as point masses. Since the droplet-scale flows are not resolved, closure models of the force and heat transfer between the droplets and the surrounding flow are needed. As the Weber numbers of these droplets/bubbles are typically small, they are not much different from solid particles. Thus, modelling efforts on force and heat transfer for dispersed multiphase flows or particle-laden flows are directly applicable (Magnaudet & Eames Reference Magnaudet and Eames2000; Balachandar & Eaton Reference Balachandar and Eaton2010; Ling, Parmar & Balachandar Reference Ling, Parmar and Balachandar2013; Ling, Balachandar & Parmar Reference Ling, Balachandar and Parmar2016). However, the above modelling efforts have not yet been able to resolve the fundamental challenge of sub-grid modelling of atomization, i.e., how to accurately represent the under-resolved formation of sub-scale droplets. It is expected that statistics of droplets from different formation mechanisms will vary significantly. The size distribution of droplets generated in ligament breakup owing to RP instability will, for example, be different from that for droplets produced in a secondary breakup.
In the literature, there are also simulations which combine interface-capturing methods and LES filtering to the turbulent gas flows (Labourasse et al. Reference Labourasse, Lacanette, Toutant, Lubin, Vincent, Lebaigue, Caltagirone and Sagaut2007; Larocque et al. Reference Larocque, Vincent, Lacanette, Lubin and Caltagirone2010; Lakehal, Labois & Narayanan Reference Lakehal, Labois and Narayanan2012; Aniszewski Reference Aniszewski2016). These modelling efforts are mainly focused on the sub-scale surface tension effect since the small-scale variation of curvature is under-resolved. The robustness and accuracy of these models in capturing flows with significant topological changes are still to be explored. We believe that a viable sub-grid modelling approach will need to be event based. In other words, the model has to be able to identify the droplet formation event and the corresponding breakup mechanism based on topological configurations of the macro-scale liquid structures. To develop sub-grid models like this, rigorous data of the droplets statistics covering a sufficiently large number of events has to be collected from fully resolved simulations.
1.6 Effect of mesh resolution
While simulations of bubbles and drops retaining their identities have been shown to produce fully converged solutions (Lu & Tryggvason Reference Lu and Tryggvason2013; Dodd & Ferrante Reference Dodd and Ferrante2016), liquid–gas multiphase flow simulations where the topology changes through breakup and coalescence generally result in spontaneously generated small-scale features that are difficult to resolve. This is particularly true for almost all simulations of large-scale atomization (Shinjo & Umemura Reference Shinjo and Umemura2010; Le Chenadec & Pitsch Reference Le Chenadec and Pitsch2013; Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016; Agbaglah et al. Reference Agbaglah, Chiodi and Desjardins2017; Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017). The general consensus among researchers has been that while the small-scale physics are under-resolved, the large-scale flow remains correct. Since small droplets and filaments contain little mass, leaving them unresolved should have only minor impact on the overall results. We have recently started to examine this assumption in more detail, by extensive grid refinement studies varying from 8 million to 4 billion cells (number of cells to resolve the initial liquid stream thickness varying from 32 to 256) (Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017). While the results show that some of the large-scale statistics converge, considerable sensitivity on the resolution has also been observed, such as for the droplet size distribution. In particular, small-scale instabilities can generate drops larger than the most-unstable wave length and the error resulting from not resolving the smallest scales fully thus manifests itself at much larger scales.
1.7 Goals of this study
The purpose of the present study is to answer the following important questions for simulations of spray formation in a two-phase mixing layer between parallel gas and liquid streams.
(i) What are the Kolmogorov and Hinze scales in the present two-phase mixing layer and do they effectively represent the flow physics in atomization?
(ii) What is the mesh requirement to fully resolve turbulent atomization?
(iii) Will the large-scale multiphase turbulence statistics be affected if the small scale versions are under-resolved?
(iv) How does the interfacial instability affect the multiphase turbulence development?
Particular attention will be focused on obtaining the statistics of multiphase turbulence and on the impact of the upstream interfacial instability on the turbulence. As an extension to our previous work (Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017), the simulation for the most refined mesh (M3) has been run for about twice as long, so that the statistically converged multiphase turbulence statistics, in particular those of higher order, can be obtained.
2 Methodology
2.1 Governing equations
The one-fluid approach is employed to resolve the two-phase flow, where the liquid and gas phases are treated as one fluid with material properties (such as density and viscosity) that change abruptly across the interface. The incompressible two-phase flows are governed by the Navier–Stokes equations with surface tension,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}t}+u_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right)=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left[\unicode[STIX]{x1D707}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)\right]+f_{s,i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}$
and
$\unicode[STIX]{x1D707}$
are the fluid density and viscosity,
$u$
and
$p$
the velocity and pressure fields. The surface tension term is expressed as
$$\begin{eqnarray}\displaystyle f_{s,i}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FF}_{s}n_{i}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70E}$
is the surface tension coefficient (assumed to be constant here), and
$\unicode[STIX]{x1D705}$
and
$n_{i}$
are the local curvature and unit normal of the interface. The surface tension is a singular term, with a Dirac distribution function
$\unicode[STIX]{x1D6FF}_{s}$
localized on the interface.
The volume fraction
$c$
is introduced to distinguish between the two different phases. Here,
$c=1$
in computational cells with only the liquid phase, and its time evolution satisfies the advection equation (Hirt & Nichols Reference Hirt and Nichols1981)
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+u_{i}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x_{i}}=0. & & \displaystyle\end{eqnarray}$$
The fluid density and viscosity are calculated based on the arithmetic mean as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}=c\unicode[STIX]{x1D70C}_{l}+(1-c)\unicode[STIX]{x1D70C}_{g}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}=c\unicode[STIX]{x1D707}_{l}+(1-c)\unicode[STIX]{x1D707}_{g}. & \displaystyle\end{eqnarray}$$
Detailed discussion about using arithmetic or harmonic means for viscosity has been given by Boeck et al. (Reference Boeck, Li, López-Pagés, Yecko and Zaleski2007). It is shown that both viscosity methods yield similar results for sufficiently high mesh resolution.
2.2 Numerical methods
The governing equations are solved by the open source code PARIS-Simulator. The details of the numerical methods implemented in PARIS-Simulator can be found in previous works (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011; Ling et al. Reference Ling, Zaleski and Scardovelli2015; Bnà et al. Reference Bnà, Manservisi, Scardovelli, Yecko and Zaleski2016; Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017) and the code webpage. (The PARIS-Simulator Code, available from http://www.ida.upmc.fr/∼zaleski/paris.) Only the numerical aspects that are relevant to the present study are summarized here.
The Navier–Stokes equations (2.1)–(2.2), are solved by the finite volume method on a staggered grid. The fields are discretized using a fixed regular cubic grid (with cell size
$\unicode[STIX]{x1D6E5}$
) and we use a projection method for the time stepping to incorporate the incompressibility condition (Chorin Reference Chorin1968). The temporal integration is done by a second-order predictor–corrector method. The interface is tracked using a volume-of-fluid (VOF) method with the mixed Young’s-centred implementation of Aulisa et al. (Reference Aulisa, Manservisi, Scardovelli and Zaleski2007) to determine the normal vector and the Lagrangian-explicit scheme of Li (Reference Li1995) for the VOF advection (Scardovelli & Zaleski Reference Scardovelli and Zaleski2003). The advection of momentum near the interface is implemented in a manner consistent with the VOF advection, similar to the methods of Rudman (Reference Rudman1998) and Vaudor et al. (Reference Vaudor, Ménard, Aniszewski, Doring and Berlemont2017). The superbee limiter is applied in the flux calculation (Roe Reference Roe1986). The viscous term is treated explicitly with a second-order centered difference scheme. Curvature is computed using the height-function method of Popinet (Reference Popinet2009). Surface tension is computed from the curvature by a balanced continuous-surface-force method (Renardy & Renardy Reference Renardy and Renardy2002; Francois et al.
Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006; Popinet Reference Popinet2009). To capture the dynamics of under-resolved droplets less erroneously than by just quasi-fragment VOF patches, droplets of size smaller than four cells are converted into Lagrangian point-particles and are traced under the one-way coupling approximation, following the approach of Ling et al. (Reference Ling, Zaleski and Scardovelli2015).
2.3 Simulation set-up
2.3.1 Computational domain
As shown in figure 1 the computational domain is a rectangular cuboid. The domain is initially filled with stationary gas (at
$t=0$
) and then liquid and gas streams progressively enter it. The
$x$
-coordinate is aligned with the stream velocity, whereas
$y$
and
$z$
are along the height and width of the stream. The thicknesses of the liquid and gas streams at the inlet are represented by
$H_{l}$
and
$H_{g}$
, respectively. Here,
$H_{l}$
is chosen to be the characteristic length scale. Then the length (
$x$
), height (
$y$
) and width (
$z$
) of the domain are taken to be
$L_{x}=16H_{l}$
,
$L_{y}=8H_{l}$
and
$L_{z}=2H_{l}$
, respectively. The thickness and the length of the separator plate are denoted as
$l_{y}$
and
$l_{x}$
. The separator plate is included to mimic the effect of the fuel injection nozzle and the need for such a plate to accurately capture interfacial instability and wave breakups has been addressed by Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) and will also be discussed later.
To reduce the computational cost, a relatively small domain width
$L_{z}$
is used, compared with
$L_{x}$
and
$L_{y}$
. The characteristic length scale for the interfacial instability development is the vorticity layer thickness
$\unicode[STIX]{x1D6FF}_{g}$
. The current domain width is significantly larger than
$\unicode[STIX]{x1D6FF}_{g}$
, i.e.,
$L_{z}/\unicode[STIX]{x1D6FF}_{g}=16$
, and therefore is sufficient to capture the development of interfacial stability and wave formation. When the transverse instability develops at the rim further downstream, the domain width used here may not be sufficient to resolve the large wavelengths. The effect of the domain width
$L_{z}$
to the simulation results are discussed in the appendix C, in which we have shown results with a domain four times wider than the present one (namely
$L_{z}=8H_{l}$
). The results of the present and the wider domains for both low- and high-order two-phase turbulence statistics (mean velocity and dissipation) agree with the results of the present domain in general, suggesting that the important conclusions made in the present study remain valid. The discrepancy mainly lies at the unbroken liquid stream near the bottom of the domain, which indicates the constraint of the domain width indeed influences the transverse instability development and interfacial wave breakup downstream in some extent. A high-resolution simulation using a wider domain is computationally expensive and will be relegated to future work.
2.3.2 Boundary conditions
An inflow boundary condition is applied to the left of the domain (
$x=0$
), with the velocity specified as
$$\begin{eqnarray}\displaystyle u_{x=0}=\left\{\begin{array}{@{}ll@{}}\displaystyle U_{l}\,\text{erf}\,\frac{(H_{l}-y)}{\unicode[STIX]{x1D6FF}_{l}}, & 0\leqslant y<H_{l},\\ 0, & H_{l}\leqslant y<H_{l}+l_{y},\\ \displaystyle U_{g}\,\text{erf}\,\frac{[y-(H_{l}+l_{y})]}{\unicode[STIX]{x1D6FF}_{g}}\,\text{erf}\,\frac{[(H_{l}+l_{y}+H_{g})-y]}{\unicode[STIX]{x1D6FF}_{g}}, & H_{l}+l_{y}\leqslant y<H_{l}+l_{y}+H_{g},\\ 0, & \text{otherwise}.\end{array}\right. & & \displaystyle\end{eqnarray}$$
The separator plate is located at
$H_{l}\leqslant y<H_{l}+l_{y}$
. The error function, defined as
$$\begin{eqnarray}\displaystyle \text{erf}(y)=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{y}\exp (-\unicode[STIX]{x1D712}^{2})\,\text{d}\unicode[STIX]{x1D712}, & & \displaystyle\end{eqnarray}$$
is known to be the exact solution of the first Stokes problem and is employed to represent the vorticity layers on the top and bottom boundaries of the gas stream and the top of the liquid stream, following the previous works (Otto et al.
Reference Otto, Rossi and Boeck2013; Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013). The thickness of the vorticity layers at the top and bottom boundaries of the gas stream is denoted by
$\unicode[STIX]{x1D6FF}_{g}$
and that for the vorticity (boundary) layer at the top of the liquid stream by
$\unicode[STIX]{x1D6FF}_{l}$
. (There is no vorticity layer at the bottom of the liquid stream since the domain bottom is considered to be a slip wall.) For the velocity profile defined here, the displacement boundary layer thickness is
$\unicode[STIX]{x1D6FF}/\sqrt{\unicode[STIX]{x03C0}}$
and the boundary layer thickness corresponding to
$u=0.99U_{g}$
is about
$2\unicode[STIX]{x1D6FF}$
(Ling et al.
Reference Ling, Fuster, Zaleski and Tryggvason2017). The volume fraction function at the inlet is specified as
$$\begin{eqnarray}\displaystyle c_{x=0}=\left\{\begin{array}{@{}ll@{}}1, & 0\leqslant y<H_{l},\\ 0, & \text{otherwise}.\end{array}\right. & & \displaystyle\end{eqnarray}$$
The bottom of the domain (
$y=0$
) is taken to be a slip wall and periodic boundary conditions are used at the back and front boundaries (
$z=0$
and
$z=L_{z}$
).
To minimize the effect of the finite size of the domain, additional attention is required for the boundary conditions at the top (
$y=L_{y}$
) and the right of the domain (
$x=L_{x}$
). In general, there are two options for the top boundary: (1) symmetric boundary (or slip wall) (Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013; Agbaglah et al.
Reference Agbaglah, Chiodi and Desjardins2017); and (2) free boundary that allows the gas to freely enter or leave the boundary (Taub et al.
Reference Taub, Lee, Balachandar and Sherif2013; Ling et al.
Reference Ling, Zaleski and Scardovelli2015, Reference Ling, Fuster, Zaleski and Tryggvason2017; Almagro, Garcia-Villalba & Flores Reference Almagro, Garcia-Villalba and Flores2017). If the former condition is used, a recirculating flow will form on top of the parallel streams (Agbaglah et al.
Reference Agbaglah, Chiodi and Desjardins2017). The recirculation is less favorable since obviously it may influence the physics of interest, such as carrying coherent structures downstream back to the inlet, unless the domain is so large that the effect of the recirculation becomes negligibly weak (Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013). Owing to high computational cost, a relatively small domain is used in the present study, although
$L_{y}$
and
$L_{x}$
are already 8 and 16 times of the initial liquid stream thickness and are large enough to capture the physics near the parallel streams. Therefore, the free boundary conditions is chosen for the top boundary in the present set-up to minimize the effect of recirculation. Since we have used the free boundary condition on the top, the outlet condition on the right surface of the domain requires the convective velocity to be specified. (If a pressure outflow boundary condition is invoked, then the flow is under constrained and may exit at the top boundary, breaking the parallelism of the two streams.) The outflow velocity profile imposed at the right of the domain will affect the mean flow. To mimic the development of the gas stream, we specify the outflow velocity based on the average velocity of a planar turbulent jet (Pope Reference Pope2000),
$$\begin{eqnarray}\displaystyle u_{x=L_{x}}=\left\{\begin{array}{@{}ll@{}}U_{c}\,\text{sech}^{2}(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D709}), & \text{if }y>H_{l}+1/2H_{g},\\ U_{c}, & \text{otherwise},\end{array}\right. & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D709}=y/y_{1/2}$
, and
$y_{1/2}$
is the half-width of the turbulent jet. The convective velocity
$U_{c}$
is determined by mass balance so that the flux into the domain given in (2.7) is equal to that leaving the domain. The variation
$y_{1/2}$
in
$x$
is found to be linear, namely,
$\text{d}y_{1/2}/\text{d}x=S$
. The parameters
$S$
and
$\unicode[STIX]{x1D6FC}$
are constant, the values of which are given as 0.10 and 0.88, respectively (Pope Reference Pope2000). A Neumann boundary condition
$\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}x=0$
is applied for the volume fraction function at the right boundary.
It is noted that the outflow velocity profile used here does not represent the exact condition at the outlet of the domain, since the flow is turbulent and time dependent. Equation (2.10) is thus only an approximation to the mean flow at the outlet. It is expected that the overall mean flow can be affected to a certain extent by the outflow velocity profile, in particular, a small region near the outlet will be influenced. Nevertheless, it is shown from simulation results that the overall boundary conditions applied here can effectively minimize the recirculation on the top of the parallel streams (see figure 1) and also convect the vortices and droplets out of the domain.
To thoroughly examine the effect of the present boundary conditions on the simulation results, we have performed simulations with a larger domain (
$L_{x}$
and
$L_{y}$
are 1.5 times those of the current set-up) on a coarser mesh. The details of the tests for different domain sizes are given in the appendix C. The results show that the key conclusions made in the present study are not influenced by the boundary conditions and the domain size.
Physical parameters.
2.3.3 Physical parameters
The material properties of the two fluids (
$\unicode[STIX]{x1D70C}_{l}$
,
$\unicode[STIX]{x1D70C}_{g}$
,
$\unicode[STIX]{x1D707}_{l}$
,
$\unicode[STIX]{x1D707}_{g}$
,
$\unicode[STIX]{x1D70E}$
) and the injection conditions of the two streams (
$U_{l}$
,
$U_{g}$
,
$H_{l}$
,
$H_{g}$
,
$\unicode[STIX]{x1D6FF}_{l}$
,
$\unicode[STIX]{x1D6FF}_{g}$
), values given in table 1, fully characterize the resulting multiphase flow. To simplify the analysis, we take
$H_{l}=H_{g}+l_{y}$
and
$\unicode[STIX]{x1D6FF}_{l}=\unicode[STIX]{x1D6FF}_{g}$
. Following the Buckingham
$\unicode[STIX]{x03C0}$
theorem, the dimensional physical parameters in table 1 are expressed in dimensionless form as shown in table 2.
The liquid properties used here are the same as those of water. The gas is similar but not identical to pressurized air. Instead of using exact air properties in experiments (Matas et al.
Reference Matas, Marty and Cartellier2011; Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013), we consider a case of moderate density ratio (Ling et al.
Reference Ling, Fuster, Zaleski and Tryggvason2017) that is less expensive for numerical simulation. (As will be shown later, even for this ‘easier’ case, we barely reach fully resolved results, and thus a DNS at the exact experiment condition will be exceedingly expensive for currently available computer power.) Therefore, the fluid properties and injection conditions here are not chosen to match any realistic fuel injection condition. A larger gas density is adopted here so that the liquid-to-gas density ratio is equal to 20. The gas viscosity is chosen here so that kinematic viscosity for the two phases are the same. The dynamic pressure ratio
$M$
has been shown to be the primary parameter determining the macroscale behaviour of a two-phase mixing layer (Lasheras & Hopfinger Reference Lasheras and Hopfinger2000) and whether the interfacial instability is absolute or convective (Otto et al.
Reference Otto, Rossi and Boeck2013; Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013). To place the interfacial instability in the absolute instability regime, a large gas-to-liquid dynamic pressure ratio is needed, and in the present study
$M$
is taken to be 20. Since the liquid-to-gas density ratio
$r$
and viscosity ratio
$m$
are all equal to 20, we referred to this case of atomization in a two-phase mixing layer as the ‘A20’ case.
Key dimensionless parameters.
The gas vorticity layer thickness
$\unicode[STIX]{x1D6FF}_{g}$
is the characteristic length scale for the interfacial instability (Eggers & Villermaux Reference Eggers and Villermaux2008; Matas et al.
Reference Matas, Marty and Cartellier2011). The vorticity layer thickness
$\unicode[STIX]{x1D6FF}_{g}$
varies with gas properties and injection conditions, and thus a precise value of
$\unicode[STIX]{x1D6FF}_{g}$
is generally unknown a priori. In the experiment by Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013), the injected air is at standard condition and an empirical correlation of
$\unicode[STIX]{x1D6FF}_{g}$
was given as a function of the Reynolds number of the gas stream,
$Re_{g,H}$
, which is defined as
$$\begin{eqnarray}\displaystyle Re_{g,H}=\frac{\unicode[STIX]{x1D70C}_{g}U_{g}H_{g}}{\unicode[STIX]{x1D707}_{g}}. & & \displaystyle\end{eqnarray}$$
For the present simulation, the gas properties and injection condition are different; therefore, the empirical correlation of Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) is not applicable. The vorticity layer thickness in the present simulation is an independent parameter and the value used in the present stimulation is chosen as
$\unicode[STIX]{x1D6FF}_{g}/H_{l}=1/8$
(or
$\unicode[STIX]{x1D6FF}_{g}/l_{y}=4$
).
The effect of
$\unicode[STIX]{x1D6FF}_{g}$
(and
$\unicode[STIX]{x1D6FF}_{l}$
) on the development of interfacial instability has been discussed extensively by Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) through 2D simulations. To investigate the influence of
$\unicode[STIX]{x1D6FF}_{g}$
on the interfacial wave breakup and the multiphase turbulence, simulations have to be extended to three dimensions. A parametric study of
$\unicode[STIX]{x1D6FF}_{g}$
with 3D simulations is interesting yet out of the scope of the present work.
The separator plate thickness
$l_{y}$
can have a significant impact on the interfacial instability if it is larger than or comparable with
$\unicode[STIX]{x1D6FF}_{g}$
(Fuster et al.
Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013). When
$l_{y}/\unicode[STIX]{x1D6FF}_{g}$
is sufficiently small, then the effect of
$l_{y}/\unicode[STIX]{x1D6FF}_{g}$
becomes negligible. Here, we chose
$l_{y}/\unicode[STIX]{x1D6FF}_{g}=1/4$
, which is significantly smaller than the threshold value of
$l_{y}/\unicode[STIX]{x1D6FF}_{g}=1$
given by Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) and thus the specific value of
$l_{y}$
is immaterial to the results presented here.
The Reynolds and Weber numbers of the gas stream based on the vorticity layer thickness at the inlet
$\unicode[STIX]{x1D6FF}_{g}$
, namely,
$$\begin{eqnarray}\displaystyle & \displaystyle Re_{g,\unicode[STIX]{x1D6FF}}=\frac{\unicode[STIX]{x1D70C}_{g}U_{g}\unicode[STIX]{x1D6FF}_{g}}{\unicode[STIX]{x1D707}_{g}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle We_{g,\unicode[STIX]{x1D6FF}}=\frac{\unicode[STIX]{x1D70C}_{g}U_{g}^{2}\unicode[STIX]{x1D6FF}_{g}}{\unicode[STIX]{x1D70E}}, & \displaystyle\end{eqnarray}$$
are also key dimensionless parameters for the interfacial instability (Otto et al. Reference Otto, Rossi and Boeck2013).
Summary of simulation runs.
2.3.4 Mesh resolution and time step
The fields are discretized using a fixed regular cubic grid (with cell size
$\unicode[STIX]{x1D6E5}$
). Simulations are performed on four meshes referred to as M0, M1, M2 and M3, so that M
$n$
has
$H_{l}/\unicode[STIX]{x1D6E5}=32\times 2^{n}$
points in the liquid stream layer
$H_{l}$
; see table 3. The time step in the simulation for each mesh is computed based on time step restrictions for the convection term (the Courant–Friedrichs–Lewy (CFL) condition), the diffusion term and the surface tension term,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}t\leqslant \unicode[STIX]{x0394}t_{conv}=\frac{\unicode[STIX]{x1D703}\unicode[STIX]{x1D6E5}}{u_{max}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}t\leqslant \unicode[STIX]{x0394}t_{visc}=\frac{\unicode[STIX]{x1D6E5}^{2}}{6\unicode[STIX]{x1D708}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}t\leqslant \unicode[STIX]{x0394}t_{surf}=\sqrt{\frac{(\unicode[STIX]{x1D70C}_{g}+\unicode[STIX]{x1D70C}_{l})\unicode[STIX]{x1D6E5}^{3}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}}}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D703}$
is the CFL number.
For the M3 mesh,
$\unicode[STIX]{x1D6E5}=3.15\times 10^{-6}$
m,
$\unicode[STIX]{x0394}t_{conv}=1.27\times 10^{-7}$
s (assuming
$u_{max}=U_{g}=10~\text{m}~\text{s}^{-1}$
and
$\unicode[STIX]{x1D703}=0.4$
),
$\unicode[STIX]{x0394}t_{visc}=1.63\times 10^{-6}$
, and
$\unicode[STIX]{x0394}t_{surf}=4.52\times 10^{-7}~\text{s}$
. As can be seen here the convection time step restriction is the most demanding and as a result dictates the time step in the simulation. The small time step given by the CFL condition is due to the large gas injection velocity. In the simulation, the CFL number
$\unicode[STIX]{x1D703}$
is taken to be 0.4 in general. To confirm the simulation results are independent of time step, smaller
$\unicode[STIX]{x1D703}$
such as 0.2 has also been used and it is confirmed that the time step is sufficiently small.
The domain is initially filled with stationary gas (at
$t=0$
) and then liquid and gas streams progressively enter it. It takes a time period of about
$tU_{g}/H_{l}\approx 200$
for the flow to reach a statistically steady state. The transient process has been shown in previous works (Ling et al.
Reference Ling, Fuster, Zaleski and Tryggvason2017). For the M0, M1 and M2 meshes, the simulations all start from
$t=0$
and end at about
$tU_{g}/H_{l}=1000$
, 880 and 650, respectively. For the M3 mesh, the simulation was performed using
$4.29\times 10^{9}$
cells and 16 384 processors. Owing to the high computational cost for the M3 simulation, the simulation starts from a checkpoint of the M2 simulation at about
$tU_{g}/H_{l}=200$
, and is continued only up to about
$tU_{g}/H_{l}=450$
.
The M3 simulations are split into multiple runs, which are conducted on the supercomputers CINECA-FERMI in Italy, LRZ-superMUC in Germany and TGCC-CURIE in France. The M0, M1 and M2 simulations are all performed on TGCC-CURIE. The total simulation time for all four meshes is over
$15\times 10^{6}$
CPU core-hours. The results presented correspond to the M3 mesh, unless stated otherwise.
The results of grid and statistical convergence studies, namely evaluating the effects of the mesh resolution and the averaging time on the present results, are shown in the appendices A and B.
3 Results
3.1 General behaviour
Development of interfacial waves and coherent vortical structures. (See the supplementary movie available at https://doi.org/10.1017/jfm.2018.825.)
When the two streams are injected into the domain, both of them are laminar and the gas–liquid interface is perfectly flat. As the two streams meet at the downstream end of the separator plate, the velocity difference between the two streams introduces a shear on the interface, which then triggers a Kelvin–Helmholtz instability. As a response of this shear-induced instability, an interfacial wave is formed as shown in the right column of figure 2. The shape of the wave at early stage is mainly influenced by the density ratio, as explained by Hoepffner et al. (Reference Hoepffner, Blumenthal and Zaleski2011). The wave propagates downstream with the Dimotakis velocity (1.1), which is in between the gas and liquid injection velocities (see the right column of figure 2). This is consistent with experiments by Raynal (Reference Raynal1997), Hoepffner et al. (Reference Hoepffner, Blumenthal and Zaleski2011) and Jerome et al. (Reference Jerome, Marty, Matas, Zaleski and Hoepffner2013). The amplitude of the wave grows in time. At a certain stage the wave amplitude becomes comparable with the gas stream thickness, and the interaction between the interfacial wave and the gas stream becomes strong. The interaction causes the liquid sheet pulled from the wave crest to roll and to flap and eventually the liquid sheet breaks violently.
A snapshot of the vortical structures near the interface. Vortices are visualized by the
$\unicode[STIX]{x1D706}_{2}$
criterion. The colour on the
$\unicode[STIX]{x1D706}_{2}$
iso-surfaces represents the
$z$
-component of vorticity.
At the same time, instability also develops at the gas stream vorticity layer near the interface due to the shear. Owing to the lower velocity of the liquid stream, the gas–liquid interface is seen as a deformable and wavy wall by the gas stream. The evolution of vortical structures near the interface is visualized by the
$\unicode[STIX]{x1D706}_{2}$
criterion (Jeong & Hussain Reference Jeong and Hussain1995). To distinguish the vortex rotation direction, the
$\unicode[STIX]{x1D706}_{2}$
iso-surface is coloured by the
$z$
-component of the vorticity. As a result, the red and blue vortices are aligned with the
$z$
direction and rotate counter-clockwise and clockwise, respectively. On the other hand, vortices with green colour are aligned with the stream direction. A 3D snapshot of the vortical structure is shown in figure 3. It can be observed that the vortical structures upstream of the interfacial wave are quite similar to those in a turbulent boundary layer (Wu & Moin Reference Wu and Moin2009). The laminar vorticity layer transitions to turbulence and hairpin vortices are clearly seen near the transition region. As the amplitude of the interfacial wave becomes large and acts as an obstacle to the gas flow, the flow separates at the downstream face of the wave, forming a turbulent wake. As a result, the interfacial wave is immersed in these complex turbulent vortices and the stretching and breaking of the wave take place in a fully 3D chaotic manner.
3.2 Statistics of multiphase turbulence
3.2.1 Reynolds averaging and the mean flow
The mean flow for the present problem is 2D (
$x$
–
$y$
), so averaging of quantities obtained from the DNS is conducted both temporally and spatially in the
$z$
direction. The time and spatial (in
$z$
direction) averaging operator
$\overline{()}$
is defined as
$$\begin{eqnarray}\displaystyle \overline{u}(x,y)\equiv \frac{1}{t_{1}-t_{0}}\frac{1}{L_{z}}\int _{t_{0}}^{t_{1}}\int _{0}^{L_{z}}u(x,y,z,t)\,\text{d}z\,\text{d}t, & & \displaystyle\end{eqnarray}$$
where
$t_{0}$
and
$t_{1}$
are the starting and ending time for averaging. In the present study,
$t_{0}U_{g}/H_{l}=200$
when the statistically steady state is reached, and
$t_{1}$
is the end time of the computation. The mean quantities are time-independent if
$t_{1}$
is sufficiently large.
The average liquid volume fraction and the streamwise velocity are shown in figure 4. The contour line in figure 4 corresponds to
$\overline{c}=0.5$
, which can be viewed as the ‘average’ boundary of the unbroken liquid stream. The streamwise evolutions of the profiles of
$\overline{c}$
and
$\overline{u}$
are shown in figure 5.
Average (a) liquid volume fraction and (b)
$u$
-velocity. The black curve corresponds to
$\overline{c}=0.5$
. The orange rectangle near the inlet represents the separator plate.
Mean flow profiles for (a)
$\overline{c}$
and (b)
$\overline{u}$
at different streamwise locations.
The fluctuation deviating from the average quantity is given as
$$\begin{eqnarray}\displaystyle u^{\prime }=u-\overline{u}, & & \displaystyle\end{eqnarray}$$
which is denoted by a prime
$^{\prime }$
. Conventionally, the Reynolds stress tensor divided by density (or often simply referred to as Reynolds stress (Pope Reference Pope2000)) is expressed as the velocity covariances,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{RA}/\unicode[STIX]{x1D70C}=\overline{u_{i}^{\prime }u_{j}^{\prime }}, & & \displaystyle\end{eqnarray}$$
where the subscript
$RA$
represent Reynolds averaging. The results for
$\overline{u_{i}^{\prime }u_{j}^{\prime }}$
are shown in figure 6. It can be observed that the maximum magnitude of
$\overline{u^{\prime }u^{\prime }}/U_{g}^{2}$
is about 0.12, which is much larger than those of other components.
It should be noted that, the Reynolds stress tensor expression given in (3.3) is strictly valid only for single-phase incompressible flows. For the present problem that involves two fluids of different density, the Reynolds stress tensor based on Favre averaging will better characterize the turbulent two-phase flows, discussed in the following section.
Reynolds stresses (divided by density) for Reynolds averaging. The black curve corresponds to
$\overline{c}=0.5$
. The orange rectangle near the inlet represents the separator plate.
3.2.2 Favre averaging and averaged momentum equation
For the present problem, the density at a given location may exhibit temporal fluctuations, although the density in each phase remains constant. The density fluctuations are due to the unsteady motion of the gas–liquid interface and, thus, are generally strong near the gas–liquid interface. As a result, the average density
$\overline{\unicode[STIX]{x1D70C}}$
varies spatially. For turbulent flows with variable density, such as compressible turbulent flows (Huang, Coleman & Bradshaw Reference Huang, Coleman and Bradshaw1995), the Favre-averaging (density-weighted-averaging) technique is commonly employed to develop the averaged equations. The Favre averaging or decomposition have also been applied to gas–liquid flow for turbulence statistics analysis and model development (Vallet, Burluka & Borghi Reference Vallet, Burluka and Borghi2001; Demoulin et al.
Reference Demoulin, Beau, Blokkeel, Mura and Borghi2007; Mortazavi et al.
Reference Mortazavi, Le Chenadec, Moin and Mani2016).
The Favre-averaging operator
$\tilde{()}$
is defined as
$$\begin{eqnarray}\displaystyle \tilde{u} =\overline{\unicode[STIX]{x1D70C}u}/\overline{\unicode[STIX]{x1D70C}}, & & \displaystyle\end{eqnarray}$$
and the fluctuation away from the Favre-averaged quantity can be expressed as
$$\begin{eqnarray}\displaystyle u^{\prime \prime }=u-\tilde{u} , & & \displaystyle\end{eqnarray}$$
which is denoted by a double prime
$^{\prime \prime }$
. It can be easily shown that
$\overline{\tilde{u} }=\tilde{u}$
and
$\overline{u^{\prime \prime }}=\overline{u}-\tilde{u}$
.
The two-dimensional averaged momentum equation can be written as
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{i}\tilde{u} _{j}}{\unicode[STIX]{x2202}x_{j}}=-\frac{\unicode[STIX]{x2202}\overline{p}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\overline{\left[\unicode[STIX]{x1D707}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)\right]}+\overline{f_{s,i}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70F}_{ij}$
is the Reynolds stress tensor
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{ij}=-\overline{\unicode[STIX]{x1D70C}}\widetilde{u_{i}^{\prime \prime }u_{j}^{\prime \prime }}. & & \displaystyle\end{eqnarray}$$
In two-phase flows with a sharp interface, the viscosity can be expressed in terms of the Heaviside function as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{l}H+\unicode[STIX]{x1D707}_{g}(1-H), & & \displaystyle\end{eqnarray}$$
where
$H=1$
and 0 in liquid and gas, respectively. The volume fraction function
$c$
in VOF methods is identical to
$H$
in cells with either liquid or gas; while for cells with an interface,
$c$
is the integral of
$H$
divided by the cell volume. The viscosity computed by (2.6) is exact in cells with only liquid or gas and is a numerical approximation in cells with an interface. As a result, the average viscosity
$\overline{\unicode[STIX]{x1D707}}$
can be related to
$\overline{c}$
as
$$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D707}}=(\unicode[STIX]{x1D707}_{l}-\unicode[STIX]{x1D707}_{g})\overline{c}+\unicode[STIX]{x1D707}_{g}, & & \displaystyle\end{eqnarray}$$
and it can be easily shown that
$$\begin{eqnarray}\displaystyle \frac{\overline{\unicode[STIX]{x1D707}u}}{\overline{\unicode[STIX]{x1D707}}}=\frac{(\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g}-1)\overline{cu}+\overline{u}}{(\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g}-1)\overline{c}+1}. & & \displaystyle\end{eqnarray}$$
Similarly, for the average density, we have
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D70C}}=(\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{g})\overline{c}+\unicode[STIX]{x1D70C}_{g}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\overline{\unicode[STIX]{x1D70C}u}}{\overline{\unicode[STIX]{x1D70C}}}=\frac{(\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}-1)\overline{cu}+\overline{u}}{(\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}-1)\overline{c}+1}. & \displaystyle\end{eqnarray}$$
In the present study,
$\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}=\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g}$
, thus
$$\begin{eqnarray}\displaystyle \frac{\overline{\unicode[STIX]{x1D707}u}}{\overline{\unicode[STIX]{x1D707}}}=\frac{\overline{\unicode[STIX]{x1D70C}u}}{\overline{\unicode[STIX]{x1D70C}}}\equiv \tilde{u} , & & \displaystyle\end{eqnarray}$$
and (3.6) can be simplified as
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{i}\tilde{u} _{j}}{\unicode[STIX]{x2202}x_{j}}=-\frac{\unicode[STIX]{x2202}\overline{p}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left[\overline{\unicode[STIX]{x1D707}}\left(\widetilde{\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}}+\widetilde{\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}}\right)\right]+\overline{f}_{s,i}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}. & & \displaystyle\end{eqnarray}$$
The Reynolds stress tensor can be computed by
$$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D70C}}\widetilde{u_{i}^{\prime \prime }u_{j}^{\prime \prime }}=\overline{\unicode[STIX]{x1D70C}u_{i}u_{j}}-\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{i}\tilde{u} _{j}, & & \displaystyle\end{eqnarray}$$
which involves third-order statistics.
It can be observed from the comparison between figures 6 and 7 that, all the four components of Reynolds stress tensors from Favre averaging,
$\widetilde{u_{i}^{\prime \prime }u_{j}^{\prime \prime }}$
, are quite similar to those from Reynolds averaging,
$\overline{u_{i}^{\prime }u_{j}^{\prime }}$
. The discrepancy between the Reynolds- and Favre-averaged quantities is mainly located in the gas–liquid mixing layer (
$y/H_{l}\sim 1$
), particularly in the region where the interfacial waves form and grow (
$4<x/H_{l}<8$
). Above the contour line of
$\overline{c}=0.5$
, the magnitudes of
$\widetilde{u_{i}^{\prime \prime }u_{j}^{\prime \prime }}$
are shown to be much lower than those of
$\overline{u_{i}^{\prime }u_{j}^{\prime }}$
. The unsteady motion of the gas–liquid interface and the substantial difference in turbulence intensity in the gas and liquid streams (the liquid stream remains laminar) have a strong impact on Reynolds stresses near the interface. As in the present problem, the liquid density is significantly larger than the gas density, the mass-weighted (Favre-)averaged properties (such as Reynolds stresses) are weighted toward the liquid properties. Since the velocity fluctuations in the liquid stream are much weaker than those in the gas flow, the magnitudes of
$\widetilde{u_{i}^{\prime \prime }u_{j}^{\prime \prime }}$
become smaller than
$\overline{u_{i}^{\prime }u_{j}^{\prime }}$
in the gas–liquid mixing layer. As there is no density fluctuations in the gas–gas mixing layer in general (except when the interfacial wave occasionally invades into the gas–gas mixing layer), the difference between the results from Reynolds and Favre averaging is generally small.
As shown in (3.6) that the Favre-averaging technique allows one to write the mean flow momentum equation without including density fluctuations. This is an important useful feature for two-phase turbulence modelling as already shown by Vallet et al. (Reference Vallet, Burluka and Borghi2001). A more detailed analysis and modelling of Favre-averaged Reynolds stress tensor are of interest, yet which is out of the scope of the present work.
Reynolds stresses (divided by density) for Favre averaging. The black curve corresponds to
$\overline{c}=0.5$
. The orange rectangle near the inlet represents the separator plate.
3.2.3 Turbulent kinetic energy budget
The equation for the kinetic energy of the instantaneous flow can be obtained by multiplying the momentum equation (2.1) by
$u_{i}$
, giving
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}\frac{1}{2}u_{i}u_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}u_{j}\frac{1}{2}u_{i}u_{i}}{\unicode[STIX]{x2202}x_{j}} & = & \displaystyle -\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}u_{i}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left[\unicode[STIX]{x1D707}u_{i}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)\right]\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D707}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+f_{s,i}u_{i},\end{eqnarray}$$
where
$u_{i}u_{i}/2$
is the kinetic energy per unit mass. (Hereafter, we simply refer kinetic energy per unit mass as ‘kinetic energy’ unless otherwise specified.)
Similarly, we can obtain the equation for the kinetic energy of the mean flow by multiplying the mean momentum equation (3.6) by
$\tilde{u} _{i}$
,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{j}\frac{1}{2}\tilde{u} _{i}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}} & = & \displaystyle -\frac{\unicode[STIX]{x2202}\overline{p}}{\unicode[STIX]{x2202}x_{i}}\tilde{u} _{i}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left[\overline{\unicode[STIX]{x1D707}}\tilde{u} _{i}\left(\frac{\unicode[STIX]{x2202}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\tilde{u} _{j}}{\unicode[STIX]{x2202}x_{i}}\right)\right]\nonumber\\ \displaystyle & & \displaystyle -\,\overline{\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x2202}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\tilde{u} _{j}}{\unicode[STIX]{x2202}x_{i}}\right)\frac{\unicode[STIX]{x2202}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}+\overline{f}_{s,i}\tilde{u} _{i}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}-\unicode[STIX]{x1D70F}_{ij}\frac{\unicode[STIX]{x2202}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}.\end{eqnarray}$$
The turbulent kinetic energy (TKE),
$k$
, is defined as
$$\begin{eqnarray}\displaystyle k=\widetilde{u_{i}^{\prime \prime }u_{i}^{\prime \prime }}/2 & & \displaystyle\end{eqnarray}$$
and is equal to the trace of the tensor
$\widetilde{u_{i}^{\prime \prime }u_{j}^{\prime \prime }}$
, the components of which are already shown in figure 7.
The equation for TKE can be obtained by subtracting (3.17) from the averaged (3.16),
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle -\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{j}k}{\unicode[STIX]{x2202}x_{j}}-\overline{\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x_{i}}u_{i}^{\prime \prime }}-\frac{\unicode[STIX]{x2202}\frac{1}{2}\overline{\unicode[STIX]{x1D70C}}\widetilde{(u_{j}^{\prime \prime }u_{i}^{\prime \prime }u_{i}^{\prime \prime })}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left[\overline{\unicode[STIX]{x1D707}}\widetilde{u_{i}^{\prime \prime }\left(\frac{\unicode[STIX]{x2202}u_{i}^{\prime \prime }}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}^{\prime \prime }}{\unicode[STIX]{x2202}x_{i}}\right)}\right]\nonumber\\ \displaystyle & & \displaystyle -\,\overline{\unicode[STIX]{x1D707}}\widetilde{\left(\frac{\unicode[STIX]{x2202}u_{i}^{\prime \prime }}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}^{\prime \prime }}{\unicode[STIX]{x2202}x_{i}}\right)\frac{\unicode[STIX]{x2202}u_{i}^{\prime \prime }}{\unicode[STIX]{x2202}x_{j}}}-\unicode[STIX]{x1D70F}_{ij}\frac{\unicode[STIX]{x2202}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}+\overline{f_{s,i}^{\prime }u_{i}^{\prime \prime }}.\end{eqnarray}$$
Similar TKE equations have also been presented by Vallet et al. (Reference Vallet, Burluka and Borghi2001) and Mortazavi et al. (Reference Mortazavi, Le Chenadec, Moin and Mani2016). The terms on the right-hand side are advection, pressure diffusion, turbulent diffusion, viscous diffusion, dissipation, production and surface-tension-induced diffusion, respectively. The profiles of these terms at different streamwise locations are shown in figure 8. The magnitudes of all the terms generally decrease when the sampling location moves downstream. The downstream results are more noisy (which may be due to the fact that the averaging time is still not long enough), but their contribution to the overall turbulence statistics is relatively small.
The TKE budget terms can be further averaged over the domain height
$L_{y}$
to obtain a 1D distribution of TKE budget along the streamwise direction as shown in figure 9. Owing to the existence of a large number of droplets, the term due to surface tension is generally very noisy, in particular in the downstream region where interfacial waves break into droplets, see figure 9(b). The pressure-diffusion term in (3.19) includes two contributions: the pressure fluctuations due to turbulent motion and those due to surface tension at the interface. Similar to the surface-tension term in TKE budget, the pressure diffusion also exhibits significant fluctuations (similar magnitude but with an opposite sign) downstream. Note that these fluctuations are mainly induced by the Laplace pressure at droplet interfaces instead of turbulence. To identify the pressure diffusion of TKE only related to turbulent flow motion, we plot the pressure diffusion without the contribution of surface tension (namely taking away the Laplace pressure from the total pressure), as shown in figure 9(a), the profile of which is seen to be much smoother. Furthermore, since the pressure fluctuation due to the contribution of Laplace pressure is present as long as there are interfaces. Even in the region with smooth interfaces without droplets (such as, e.g.,
$3<x/H_{l}<6$
), the reduction of the magnitude of the pressure diffusion by removing the contribution of Laplace pressure is also profound.
As the TKE budget terms involve high-order statistics, it is more difficult to obtain converged solutions. Results for the 1D TKE budget obtained from different meshes are shown figure 9(c). The magnitudes of the advection and the production terms are generally significantly larger than the other three terms, so three separate figures (see figure 9
d–f) are plotted to show the effect of mesh resolution on the dissipation, the pressure diffusion and the turbulent diffusion terms. It can be observed that the dissipation and the pressure terms are more sensitive to the mesh resolution than other terms. In particular, minimum dissipation decreases from about
$-$
0.0005 to
$-$
0.0016 when the mesh is refined from M0 to M2. The M2 and M3 results for the dissipation agree quite well. Similar observation can be made for the pressure term. While the cell size decreases from M0 to M2, the pressure diffusion increases substantially. The results of the M2 and M3 are similar, although due to the noise in the M3 results the agreement is not as good as for the dissipation. The generally good agreement between the M2 and M3 results of high-order turbulence statistics indicates that the M3 mesh is adequate to resolve the turbulence in the present problem. (Further evidence for this conclusion is to be given later based on the enstrophy calculation and the estimated Kolmogorov scale.)
TKE budget at different streamwise locations. The TKE terms in (3.19) are normalized by
$\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l}$
.
One-dimensional TKE budget along the streamwise direction. The TKE terms in (3.19) are normalized by
$\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l}$
and averaged over domain height
$L_{y}$
. (a,b) TKE budget with and without the Laplace pressure contribution subtracted from the pressure diffusion term, respectively. (c) TKE budget computed with different mesh resolutions. (d–f) Closeups of the dissipation, pressure, and turbulent diffusion terms shown in (c).
It can be observed from figure 8(a) that, near the inlet at
$x/H_{l}=4$
, TKE production in the gas–liquid mixing layer is much stronger than the gas–gas counterpart. This is consistent with previous observations in figure 2 that the two-phase mixing layer is more unstable and transits to turbulence earlier. When moving downstream, the interfacial wave grows and vortices generated on the wave interact with the gas–gas mixing layer, accelerating its transition to turbulence. At
$x/H_{l}=6$
, the gas–gas mixing layer produces TKE comparable with the gas–liquid counterpart. The TKE for both mixing layers are diffused effectively and at
$x/H_{l}=8$
, the two mixing layers merge together although the results are somehow noisy. A smoother representation of the TKE budget at the downstream region would require a long averaging time, which in turn can be achieved by running the simulation for a much longer time. A longer simulation with M3 mesh is beyond the current resources available to us and will be relegated to future works.
3.2.4 Turbulence dissipation
The distribution of the turbulence dissipation, denoted as
$\unicode[STIX]{x1D716}$
, for the different meshes is shown in figure 10. The results for different mesh resolutions are plotted with the same legend. (Further results of grid refinement studies can be found in the appendix A, see figure 19
d.) It is clear that a fine mesh is required to capture the dissipation. While the M0 and M1 meshes underpredict turbulence dissipation, the M2 and M3 meshes yield similar results. The turbulence dissipation is generally located at the gas–liquid mixing layer. In the region of where the dissipation is larger (i.e.,
$4<x/H_{l}<6$
), there remains a small discrepancy between the M2 and M3 results. More results for grid-refinement studies are shown in figure 19 in appendix A. The difference between the M2 and M3 results are obviously much smaller than that between the M1 and M2 results. Therefore, although it would require a mesh finer than M3 (such as M4) to fully confirm grid independence, we believe the M3 results of dissipation presented in figure 10(a) are not far from the grid-converged solution.
Distribution of TKE dissipation for different mesh resolutions. The black curve corresponds to
$\overline{c}=0.5$
. The orange rectangle near the inlet represents the separator plate.
Profiles of TKE dissipation along the lines (a)
$y/H_{l}=1$
and (b)
$x/H_{l}=4$
, 6, 8 and 10.
The profiles of
$\unicode[STIX]{x1D716}$
along the lines
$y/H_{l}=1$
and
$x/H_{l}=4$
, 6 and 8 are plotted in figure 11. As shown in figure 11(a) the magnitude of turbulence dissipation starts to increase at about
$x/H_{l}=2$
where the interfacial wave starts to develop and the laminar vorticity layer transits to turbulence. The dissipation grows along
$x$
as turbulence develops and reaches a maximum of about
$\unicode[STIX]{x1D716}_{max}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})=-0.01$
at about
$x/H_{l}=5.5$
. After that,
$\unicode[STIX]{x1D716}$
decreases gradually. Near the outlet
$x/H_{l}=14.5$
,
$\unicode[STIX]{x1D716}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})=-0.002$
. From figure 11(b) it is seen that the distribution of
$\unicode[STIX]{x1D716}$
is initially similar to a Gaussian profile and symmetric about the line
$y/H_{l}=1$
. As the gas–liquid mixing layer develops, the profile of
$\unicode[STIX]{x1D716}$
expands in the
$y$
direction and loses its symmetry owing to the influence of the bottom wall. The bottom boundary of the non-zero
$\unicode[STIX]{x1D716}$
region is aligned with the contour line
$\overline{c}=0.5$
. The top boundary, e.g. defined as
$\unicode[STIX]{x1D716}=20\,\%\unicode[STIX]{x1D716}_{max}$
, is about a straight line with the slope
$\text{d}y/\text{d}x=0.25$
. This expansion with a constant slope ends at about
$x/H_{l}=8$
, where the interfacial wave amplitude becomes comparable with the stream thickness and the two mixing layers merge. Then the distribution of
$\unicode[STIX]{x1D716}$
becomes more uniform within the two merged mixing layers (
$0<y/H_{l}<2$
) as shown in figures 10(a) and 11(b).
When details of the turbulent flow are unknown, a simple estimate of
$\unicode[STIX]{x1D716}$
is often made based on the integral velocity
$U_{0}$
and length scale
$l_{0}$
as
$$\begin{eqnarray}\displaystyle \frac{|\unicode[STIX]{x1D716}|}{\unicode[STIX]{x1D70C}_{g}}\approx \frac{U_{0}^{3}}{l_{0}}. & & \displaystyle\end{eqnarray}$$
If we take
$U_{0}=U_{D}$
and
$l_{0}=H_{l}$
, then
$|\unicode[STIX]{x1D716}|/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})\approx U_{D}^{3}/U_{g}^{3}$
. For the current problem with larger
$M$
and
$\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}$
, the Dimotakis speed can be approximated as
$U_{D}\approx U_{g}\sqrt{(\unicode[STIX]{x1D70C}_{g}/\unicode[STIX]{x1D70C}_{l})}$
, then
$|\unicode[STIX]{x1D716}|/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})\approx U_{D}^{3}/U_{g}^{3}\approx (\unicode[STIX]{x1D70C}_{g}/\unicode[STIX]{x1D70C}_{l})^{3/2}=0.011$
, which is close to the maximum magnitude of dissipation obtained in simulation (see figure 11). It can be also proved that, if the gas inflow velocity
$U_{g}$
is used as
$U_{0}$
, the
$\unicode[STIX]{x1D716}$
will be significantly overestimated if
$H_{l}$
remains to be used as the length scale. This seems to indicate that the interfacial wave advection speed
$U_{D}$
is better than the inflow gas stream velocity
$U_{g}$
in characterizing the integral scale of the turbulent flow motion.
3.2.5 Estimates of Kolmogorov and Hinze scales
With
$\unicode[STIX]{x1D716}$
obtained above, we can estimate the Kolmogorov length scale in the gas–liquid mixing layer. The expression of the Kolmogorov length scale is given as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}=\left(\frac{\unicode[STIX]{x1D708}_{g}^{3}}{\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70C}_{g}}\right)^{1/4}. & & \displaystyle\end{eqnarray}$$
It is shown in figure 11(a) that the maximum value of
$\unicode[STIX]{x1D716}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})$
is about 0.01. Then the corresponding Kolmogorov length scale is about
$3.0~\unicode[STIX]{x03BC}\text{m}$
(
$\unicode[STIX]{x1D702}/H_{l}=0.0038$
). In the downstream region of the gas–liquid mixing layer,
$\unicode[STIX]{x1D716}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})$
decreases to about 0.002, for which
$\unicode[STIX]{x1D702}\approx 4.5~\unicode[STIX]{x03BC}\text{m}$
(
$\unicode[STIX]{x1D702}/H_{l}=0.0056$
). According to the DNS resolution criterion given by Pope (Reference Pope2000), the smallest turbulent scales will be well resolved if
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x1D702}}\lesssim 2.1. & & \displaystyle\end{eqnarray}$$
The cell size for the M3 mesh,
$\unicode[STIX]{x1D6E5}_{M3}=3.125~\unicode[STIX]{x03BC}\text{m}$
, clearly satisfies the criterion. Even the M2 mesh cell size is close to the required resolution. This is consistent with the observation that the M2 and M3 meshes yield similar results for dissipation (see figure 19
d) and confirms that the M3 mesh is adequate to provide a resolved simulation of the present problem and the multiphase turbulence statistics presented above are grid-independent.
Based on a scaling argument focusing on the balance between the inertia force due to turbulent motion and the surface tension, the maximum stable droplet diameter for the droplet size was proposed by Kolmogorov (Reference Kolmogorov1949) and Hinze (Reference Hinze1955) as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{H}=C\left(\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}_{g}}\right)^{3/5}\left(\frac{\unicode[STIX]{x1D716}}{\unicode[STIX]{x1D70C}_{g}}\right)^{-2/5}, & & \displaystyle\end{eqnarray}$$
where
$C\approx 0.725$
is a constant and
$\unicode[STIX]{x1D702}_{H}$
is often referred to as the Hinze scale. For droplets/bubbles larger than
$\unicode[STIX]{x1D702}_{H}$
, the surface tension will not be sufficient to balance the dynamic pressure fluctuations and these droplet/bubbles will break into smaller ones. Therefore, the Hinze scale indicates the smallest droplet size which can exist in a turbulent flow.
Similar to the Kolmogorov scale, we can also estimate the Hinze scale in the present problem with the turbulence dissipation obtained in simulation. For
$\unicode[STIX]{x1D716}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})=0.01$
,
$\unicode[STIX]{x1D702}_{H}\approx 264~\unicode[STIX]{x03BC}\text{m}$
; for
$\unicode[STIX]{x1D716}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})=0.002$
at the downstream mixing layer,
$\unicode[STIX]{x1D702}_{H}\approx 502~\unicode[STIX]{x03BC}\text{m}$
.
The size distribution of droplets has been shown in our previous studies (Ling et al.
Reference Ling, Fuster, Zaleski and Tryggvason2017) and it was found that the majority of droplets generated in the mixing layer are significantly smaller than
$\unicode[STIX]{x1D702}_{H}$
obtained above. (The measured mean volume-based diameter is about
$50~\unicode[STIX]{x03BC}\text{m}$
, see figure 13(c) in the work by Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017).) Therefore, the Hinze scale does not well represent the size of droplets formed in the present problem.
The maximum stable droplet diameter from the Kolmogorov–Hinze theory assumes that the breakup of a large droplet is mainly dictated by the turbulent velocity fluctuation over a length comparable with the droplet diameter. Therefore, it is a good estimate of the droplet size when turbulence is responsible for breaking bulk liquids into small droplets. The disagreement between the Hinze scale and the droplet size in the present problem seems to indicate that, although the breakups of liquid sheets and ligaments are surrounded by turbulent vortices, the turbulent velocity fluctuations are not the dominant breakup mechanism and do not dictate the size of the droplets formed. During the disintegration of a ligament or a liquid sheet, droplets much smaller than the smallest wave length are observed. The satellite droplets generated from a ligament breakup, for example, are significantly smaller than the main droplets which are of the scale of the ligament diameter. If the generated small droplets were contained in a finite region and would coalesce, the size may experience a reverse cascade back to the Hinze scale. Nevertheless, in the present problem, the droplets are rapidly convected and dispersed downstream and coalescence is rarely observed. Therefore, the Hinze scale is not a relevant length scale in describing the smallest droplet size formed in the two-phase mixing layer.
3.2.6 Energy spectra
The temporal velocity spectra at different streamwise locations of the gas–liquid mixing layer are shown in figure 12. The purpose of showing the velocity spectra is to examine whether the inertial subrange can be identified, which in turn can show whether the turbulence at different streamwise locations is in equilibrium. The spectra of the
$u^{\prime }$
and
$v^{\prime }$
fluctuations do not show a clear range with the
$-$
5/3 slope. The
$-$
5/3 slope can be better discerned from the spectrum of
$w^{\prime }$
, namely the velocity fluctuations in the homogeneous direction. The difficulty in identifying the inertial subrange is most likely due to the moderate Reynolds number in the present problem, for which the width of the inertial range is not very large. Further than that, the results shown in figure 12 are from temporal data for a given spatial location and are, thus, quite noisy, making the identification of the inertial subrange somewhat tentative.
Velocity spectra at different locations.
To better show the velocity spectra, we plot the spatial velocity spectra in the
$z$
direction at different streamwise locations in figure 13. To reduce the noise, the spectra are averaged in time. Then the inertial subrange with a
$-5/3$
power law is more clearly revealed in all three velocity components between
$k/H_{l}=0$
and 1. The lower bound wavenumber is dictated by the domain width
$L_{z}/H_{l}=2$
. The inertial subrange seems to be more clear and wider when the sampling location moves downstream. This seems to indicate that the turbulence at the upstream location (
$x/H_{l}=4$
) is out of equilibrium, which in turn is due to the strong wave–turbulence interaction in the upstream region as shown in figure 2. After the wave breaks downstream, the turbulence equilibrates and the inertial subrange is better established.
Velocity spectra at different streamwise locations.
3.3 Interfacial instability regime and dominant frequency
Following the analysis of Otto et al. (Reference Otto, Rossi and Boeck2013), we solve the Orr–Sommerfeld equations to investigate the viscous spatio-temporal instability of the two-phase mixing layer. The details of the approach can be found in their paper (Otto et al. Reference Otto, Rossi and Boeck2013) and thus are not repeated. Here only the essential steps are briefly summarized for clarity.
The base flow is 2D and the streamwise velocity profile is taken as
$$\begin{eqnarray}\displaystyle & \displaystyle u_{g}(y-H_{l})=-U_{l}\,\text{erf}\left(\frac{y-H_{l}}{\unicode[STIX]{x1D6FF}_{l}}\right)+U_{s}\left[1+\text{erf}\left(\frac{y-H_{l}}{\unicode[STIX]{x1D6FF}_{d}}\right)\right],\quad y<H_{l}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{l}(y-H_{l})=-U_{g}\,\text{erf}\left(\frac{y-H_{l}}{\unicode[STIX]{x1D6FF}_{g}}\right)+U_{s}\left[1-\text{erf}\left(\frac{y-H_{l}}{\unicode[STIX]{x1D6FF}_{d}}\right)\right],\quad y>H_{l}, & \displaystyle\end{eqnarray}$$
where
$U_{s}$
is the interface velocity obtained from continuity of shear stresses across the interface as
$$\begin{eqnarray}\displaystyle \frac{U_{s}}{U_{a}}=\frac{\unicode[STIX]{x1D6FF}_{l}/\unicode[STIX]{x1D6FF}_{g}(1+{\mathcal{M}})+m(1-{\mathcal{M}})}{1+m}\frac{\unicode[STIX]{x1D6FF}_{d}}{\unicode[STIX]{x1D6FF}_{g}}, & & \displaystyle\end{eqnarray}$$
where
${\mathcal{M}}=(U_{g}-U_{l})/(U_{g}+U_{l})$
,
$m=\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g}$
and
$U_{a}=(U_{g}+U_{l})/2$
. The parameter
$\unicode[STIX]{x1D6FF}_{d}$
is an adjusting parameter to mimic the velocity deficit behind the separator plate. Numerical simulations and experimental data reported by Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) have confirmed that this velocity deficit is important to capture the correct transition from convective to absolute instability.
The model base flow profiles used in stability analysis for different
$\unicode[STIX]{x1D6FF}_{d}$
are shown in figure 14, which also includes the mean streamwise velocity profile at
$x/H_{l}=0.75$
in the present simulation for comparison.
Spatio-temporal viscous stability: (a) models for initial streamwise velocity profile and (b)
$\unicode[STIX]{x1D714}_{i}$
contours to show the pinching point.
The perturbation about the base flow is given in the form of streamfunction
$\unicode[STIX]{x1D713}$
and takes the form of normal modes
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{g,l}(x,y,t)=\unicode[STIX]{x1D719}(y)_{g,l}\exp (\text{i}\unicode[STIX]{x1D6FC}x-\unicode[STIX]{x1D714}t). & & \displaystyle\end{eqnarray}$$
Then the Orr–Sommerfeld equations for the
$\unicode[STIX]{x1D719}_{g}(y)$
and
$\unicode[STIX]{x1D719}_{l}(y)$
are expressed as
$$\begin{eqnarray}\displaystyle & \displaystyle \left[(-\text{i}\check{\unicode[STIX]{x1D714}}+\text{i}\check{\unicode[STIX]{x1D6FC}}{\check{U}}_{l})\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\check{y}^{2}}-\check{\unicode[STIX]{x1D6FC}}^{2}\right)-\frac{m}{r\,Re_{a,\unicode[STIX]{x1D6FF}}}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\check{y}^{2}}-\check{\unicode[STIX]{x1D6FC}}^{2}\right)^{2}-\text{i}\check{\unicode[STIX]{x1D6FC}}{\check{U}}_{l}^{\prime \prime }\right]\check{\unicode[STIX]{x1D719}}_{l}=0,\quad y<H_{l},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[(-\text{i}\check{\unicode[STIX]{x1D714}}+\text{i}\check{\unicode[STIX]{x1D6FC}}{\check{U}}_{g})\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\check{y}^{2}}-\check{\unicode[STIX]{x1D6FC}}^{2}\right)-\frac{1}{Re_{a,\unicode[STIX]{x1D6FF}}}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\check{y}^{2}}-\check{\unicode[STIX]{x1D6FC}}^{2}\right)^{2}-\text{i}\check{\unicode[STIX]{x1D6FC}}{\check{U}}_{g}^{\prime \prime }\right]\check{\unicode[STIX]{x1D719}}_{g}=0,\quad y>H_{l},\qquad & \displaystyle\end{eqnarray}$$
where
$\check{()}$
denote non-dimensional variables with
$U_{a}$
and
$\unicode[STIX]{x1D6FF}_{g}$
as typical velocity and length scales.
The branches of the imaginary part of the complex growth rate
$\unicode[STIX]{x1D714}_{i}$
are shown in the complex spatial wavenumber plane (
$\unicode[STIX]{x1D6FC}_{r}$
–
$\unicode[STIX]{x1D6FC}_{i}$
plane) in figure 14(b). The method of Bers (Reference Bers1983) is used to determine the transition from convective to absolute instability. It is observed that the two branches for
$\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D6FF}_{g}/U_{a}=0.04$
reconnect at a saddle point (
$\unicode[STIX]{x1D6FC}_{r}\unicode[STIX]{x1D6FF}_{g}=0.5$
and
$\unicode[STIX]{x1D6FC}_{i}\unicode[STIX]{x1D6FF}_{g}=-0.3$
). The value of
$\unicode[STIX]{x1D714}_{r}$
corresponding to the saddle point is the dominant angular frequency emerging in the flow field. The dominant frequency,
$f_{0}$
, can be calculated as
$f_{0}=\unicode[STIX]{x1D714}_{r}/2\unicode[STIX]{x03C0}$
. The values obtained from the theory range from
$\unicode[STIX]{x1D714}_{r}\unicode[STIX]{x1D6FF}_{g}/U_{a}\approx 0.09$
to 0.1 (
$f_{0}H_{l}/U_{g}\approx 0.06$
to 0.067) for
$\unicode[STIX]{x1D6FF}_{d}/\unicode[STIX]{x1D6FF}_{g}=0.1$
to 0.5. Similarly, the value of
$-\unicode[STIX]{x1D6FC}_{i}\unicode[STIX]{x1D6FF}_{g}\approx 0.3$
is the dominant spatial growth rate. The theoretical predictions of dominant spatial growth rate and frequency are compared with simulation results in figure 15.
In the simulations, the mixing layer thickness here is estimated as
$H_{l}-y_{c}$
, where
$y_{c}$
is the mean interfacial height corresponding to
$\overline{c}=0.5$
. The spatial growth of the mixing layer thickness is shown in figure 15(a) for different meshes. It should be noted that along the streamwise direction,
$H_{l}-y_{c}$
is first negative and then becomes positive (see figure 4
a). Since
$H_{l}-y_{c}$
is plotted on a logarithmic scale in figure 15(a), the results for negative
$H_{l}-y_{c}$
(for
$(x-l_{x})/\unicode[STIX]{x1D6FF}_{g}\lesssim 9$
) will not be shown. Strictly speaking,
$H_{l}-y_{c}$
serves as a good measure of the mixing layer thickness only when the value is not too small, i.e.,
$\log [(H_{l}-y_{c})/\unicode[STIX]{x1D6FF}_{g}]>-2$
. The rapid growth in the region
$8\lesssim (x-l_{x})/\unicode[STIX]{x1D6FF}_{g}\lesssim 11$
is due to the artefact that
$H_{l}-y_{c}$
transits from negative to positive values.
It is observed that all the M1 to M3 curves superpose for
$(x-l_{x})/\unicode[STIX]{x1D6FF}_{g}\gtrsim 12$
. An exponential spatial growth can be identified in the region
$11<(x-l_{x})/\unicode[STIX]{x1D6FF}_{g}<16$
. The exponential growth rates for the M1 to M3 meshes are quite similar. These rates agree very well with the theoretical predicted value
$-\unicode[STIX]{x1D6FC}_{i}=0.3/\unicode[STIX]{x1D6FF}_{g}$
. The exponential growth region for the M2 mesh is wider than that for the M3 mesh, which may be due to the fact that the M3 simulation has only been run for a relatively short time and there remain spatial fluctuations in
$y_{c}$
(see figure 4). After that region, a more gradual nonlinear growth is seen. The linear stability theory is valid only when the amplitude of the interfacial wave and the perturbation caused by the wave to the gas stream remain small. As the interfacial wave grows and propagates downstream, the nonlinear effect will eventually become important (such as for
$\log [(H_{l}-y_{c})/\unicode[STIX]{x1D6FF}_{g}]\gtrsim -0.5$
) and the simulation results will deviate from the linear theory. The fact that the exponential growth appears only in an intermediate region has also been observed in experiments (Matas et al.
Reference Matas, Marty and Cartellier2011) and simulations (Agbaglah et al.
Reference Agbaglah, Chiodi and Desjardins2017).
(a) Spatial growth of the mixing layer thickness measured in simulations, compared with the exponential growth predicted using the maximum spatial growth rate. The dominant frequency indicated by the numerical results of (b) the spectra of the liquid and gas enstrophy over the whole domain and (c) the spectra of the interfacial height at
$x/H_{l}=0.75$
, compared with the predicted value by stability analysis.
In the simulations, the dominant frequency for the wave formation is measured through the spectra of the interfacial height, as shown in figure 15(c). Here the interfacial height is measured near the nozzle exit (at
$x/H_{l}=0.75$
) when the wave amplitude remains small, and
$A$
denotes the Fourier transform of the interfacial height. The dominant frequency is found to be about
$f_{0}H_{l}/U_{g}\approx 0.05$
, which is quite close to the theoretical prediction (
$f_{0}H_{l}/U_{g}\approx 0.06$
to 0.067). The discrepancy between the stability theory and simulation results may be due to the fact that the effect of finite stream thickness has not been considered in the stability analysis, see Matas (Reference Matas2015). Since the interfacial wave propagates with
$U_{D}$
, the wavelength can be estimated as
$\unicode[STIX]{x1D706}=U_{D}/f_{0}\approx 4.5H_{l}$
. The wavelength estimated here is consistent with the observation in figure 2. The most unstable wavelength
$\unicode[STIX]{x1D706}$
is significantly larger than the liquid stream thickness at the inlet
$H_{l}$
.
The dominant frequency (or period) is also observed in the temporal evolution and spectra of the liquid and gas enstrophy, see figures 16 and 15(b). The liquid and gas enstrophy are computed by integration over the whole domain as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{l}=\frac{1}{2}\int _{0}^{L_{x}}\int _{0}^{L_{y}}\int _{0}^{L_{z}}c\unicode[STIX]{x1D701}_{i}\unicode[STIX]{x1D701}_{i}\,\text{d}x\,\text{d}y\,\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{g}=\frac{1}{2}\int _{0}^{L_{x}}\int _{0}^{L_{y}}\int _{0}^{L_{z}}(1-c)\unicode[STIX]{x1D701}_{i}\unicode[STIX]{x1D701}_{i}\,\text{d}x\,\text{d}y\,\text{d}z, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D701}_{i}$
is the vorticity. The dominant wave period is about
$27H_{l}/U_{g}$
and is the inverse of the dominant frequency obtained from figure 15, which is about
$fH_{l}/U_{g}\approx 0.04$
. The dominant frequency in enstrophy spectra is quite close to the dominant wave frequency measured in the simulations and is slightly lower than the theoretical prediction.
Figure 16(a) also shows the enstrophy evolution for different mesh resolutions. The observed oscillations in the enstrophy are due to the periodic formation and breakup of interfacial waves. Both the mean value and the oscillation amplitude increase when the mesh is refined from M0 to M3, while the dominant oscillation frequencies for different mesh resolutions are similar. The coarse meshes M0 and M1 significantly underpredict the enstrophy magnitude. The time-averaged enstrophy is plotted as a function of grid size in figure 16(b). The average enstrophy increases with number of cells used to resolve the initial liquid stream thickness, i.e.,
$H_{l}/\unicode[STIX]{x1D6E5}$
, until it saturates at the M3 mesh (
$H_{l}/\unicode[STIX]{x1D6E5}=256$
). This again indicates that the finest mesh, M3, is necessary and close to fully resolving the multiphase turbulence. We need to admit that the average enstrophy still increases about 10 % from the M2 to the M3 mesh. The discrepancy is consistent to the observation in figure 10. To fully confirm the grid-independence of the M3 results, a simulation with an even finer mesh is required, but this will be relegated to future work.
(a) Temporal evolution of gas enstrophy of the whole domain and (b) time-averaged liquid and gas enstrophy for different mesh resolutions.
The dominant frequency in the enstrophy spectra as shown in figure 15(b) is an important observation, since it clearly indicates that the turbulence production follows the same frequency as the interfacial instability. This is due to the fact that the formation and growth of the interfacial wave has a strong impact on the turbulence transition and development, see figures 2 and 3. The interfacial wave behaves as an obstacle to the gas stream and its interaction with the gas flow generates a large number of turbulent vortices both upstream and downstream of the interfacial wave. Therefore, the growth of the interfacial wave enables the kinetic energy transfer from the gas mean flow to the turbulent fluctuations. On the other hand, the flow remains laminar within the liquid stream. As a result there is a strong intermittency in the two-phase mixing layer.
To better illustrate the impact of interfacial wave on turbulence, we plot the temporal evolution of the root mean square (RMS) of velocity fluctuations along the
$z$
direction for different streamwise locations in figure 17. The averaging operator over the domain width (
$L_{z}$
) is denoted by
$\hat{()}$
, defined as
$$\begin{eqnarray}\displaystyle \hat{u} (x,y,t)\equiv \frac{1}{L_{z}}\int _{0}^{L_{z}}u(x,y,z,t)\,\text{d}z & & \displaystyle\end{eqnarray}$$
and the fluctuation away from the mean value is given as
$$\begin{eqnarray}\displaystyle u^{\ast }=u-\hat{u} . & & \displaystyle\end{eqnarray}$$
Low-frequency fluctuations can be seen in the temporal evolution of
$(\widehat{u^{\ast }u^{\ast }})^{1/2}$
measured upstream (
$x/H_{l}=4$
), see figure 17(a), which are clearly due to the passage of the interfacial wave. The average liquid volume fraction
${\hat{c}}$
is equal to zero most of the time, except when the interfacial wave passes the sampling location,
${\hat{c}}$
jumps up to about unity, (
${\hat{c}}=1$
indicates that the interfacial wave spans over
$L_{z}$
). The occurrence of spikes in
${\hat{c}}$
follows the period of wave formation and agree well with the dominant frequency predicted by instability theory. During the passage of the interfacial wave (within the spike),
$(\widehat{u^{\ast }u^{\ast }})^{1/2}$
drops to zero; but it jumps up to a large value after the wave passes. This is due to the turbulent flows developing on the upstream side of the interfacial wave, see figure 3. Then
$(\widehat{u^{\ast }u^{\ast }})^{1/2}$
will continue to decrease until the arrival of the subsequent wave.
Further downstream (
$x/H_{l}=8$
and 12), the wave breaks and the two mixing layers merge together, the effect of the interfacial instability frequency becomes less profound and the amplitude of low-frequency fluctuations in
$(\widehat{u^{\ast }u^{\ast }})^{1/2}$
becomes smaller.
Temporal evolution of
$(\widehat{u^{\ast }u^{\ast }})^{1/2}$
(purple solid lines) and
$\overline{c}$
(green dashed lines) at different streamwise locations.
The spectra of
$\widehat{u^{\ast }u^{\ast }}$
and
$\widehat{v^{\ast }v^{\ast }}$
are shown in figure 18. For
$x/H_{l}=4$
a dominant frequency is clearly seen at about
$fH_{l}/U_{g}\approx 0.03$
–0.04, which agrees reasonably well with the dominant frequency in the spectra for the interfacial height and enstrophy. When moving downstream at
$x/H_{l}=10$
, the spectrum function decreases with frequency smoothly, entering the inertial regime, but no dominant frequency is observed.
Figure 18 clearly shows that the integral time scale is dictated by the dominant frequency (the most unstable mode) in the interfacial instability. The interfacial wave development is the driving force and feed energy to the resulting turbulent flows near the interface. This is also consistent with the previous observation that using the Dimotakis speed as the integral velocity scale in (3.20) better captures the dissipation.
Finally, the spectra drops at about
$fH_{l}/U_{g}=10$
. With the dissipation measured in simulation, the Kolmogorov frequency can be estimated,
$$\begin{eqnarray}\displaystyle f_{\unicode[STIX]{x1D702}}=\left(\frac{\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70C}_{g}}{\unicode[STIX]{x1D708}_{g}}\right)^{1/2}. & & \displaystyle\end{eqnarray}$$
For
$\unicode[STIX]{x1D716}/(\unicode[STIX]{x1D70C}_{g}U_{g}^{3}/H_{l})=0.01$
,
$f_{\unicode[STIX]{x1D702}}H_{l}/U_{g}=8.94$
, which is close to the value measured above from the spectra.
Temporal evolution of RMS of velocity fluctuations at different streamwise locations.
4 Conclusions
DNS of a two-phase mixing layer between parallel gas and liquid streams have been performed. Particular attention has been focused on obtaining high-order statistics for the multiphase turbulence arising in the two-phase mixing layer. Extensive grid refinement studies have been carried out with four different meshes with the number of cells across the liquid stream thickness varying from 32 to 256. The finest mesh (M3) consists of about 4 billion cells and has been shown to be necessary and adequate to resolve the multiphase turbulence, yielding converged high-order statistics.
Owing to the presence of fluids with different densities, the averaged momentum equation and the TKE transport equation have been developed based on the Favre-averaging technique. The results for the mean flow, Reynolds stresses and TKE budget terms have been presented. The turbulence dissipation obtained has been used to estimate the Kolmogorov and Hinze scales in the present problem. The estimated Kolmogorov length scale is similar to the resolution of the finest mesh (M3) used in the present simulation, confirming that the smallest turbulent eddies are well captured. The Hinze scale is significantly larger than the typical size of droplets formed in atomization. The chaotic breakups of ligaments and sheets generate droplets that are much smaller than the most unstable wavelength and the rapid droplets dispersion leave few opportunities for coalescence, therefore, the Hinze scale does not seem to well represent droplet size in primary atomization.
Viscous stability analysis has also been performed on the present problem following the previous works of Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013) and Otto et al. (Reference Otto, Rossi and Boeck2013). The theory predicts that the instability of the present two-phase mixing layer is absolute, since branches on the complex wavenumber space reconnect at a ‘pinching’ point. The outcome of absolute instability is that a dominant frequency will arise. The predicted value agrees well with the dominant frequency in the spectra of interface motion. The dominant frequencies in interfacial instability and enstrophy are similar, which indicates that the interfacial wave development is strongly coupled with the turbulence development. Temporal evolutions of the RMS of velocity fluctuations for different streamwise locations are then presented. It is observed that near the inlet, there is a strong intermittence effect, i.e., RMS of velocity fluctuations drops to zero when a interfacial wave passes the sampling point and then rises up after the wave passage. The spectra of velocity fluctuation RMS exhibit a dominant frequency which matches well that in the spectra of interfacial motion. The most unstable mode of interfacial instability dictates the interfacial wave period and also the integral time scale for turbulence.
Finally the important simulation results for the two-phase mixing layer are summarized in table 4.
Summary of important results for the two-phase mixing layer.
Acknowledgements
This work has been supported by the Department of Mechanical Engineering at Baylor University in the United States and the ANR MODEMI project (ANR-11-MONU-0011) program in France. This work was granted access to the high-performance computing (HPC) resources of TGCC-CURIE and CINES-Occigen under the allocations t20152b7325, t20162b7760 and 2017tgcc0080, made by GENCI. HPC resources at CINECA and LRZ based in Italy and Germany have been used for the M3 mesh simulations, supported by PRACE (2014112610). This work has used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC05-00OR22725. The authors would also acknowledge the Texas Advanced Computing Center for providing HPC resources that have contributed to the simulation results. The Baylor High Performance and Research Computing Services (HPRCS) have been used to process the simulation results reported in this paper. We thank R. Scardovelli, W. Aniszewski, J. Lu, L. Malan and D. Jiang for their contribution to the development of the code PARIS-Simulator and also thank Dr T. Otto and T. Boeck for sharing their spatio-temporal stability code. We also appreciate the helpful discussions with A. Cartellier and J.-P. Matas. Finally, the simulation data are visualized by the software VisIt developed by the Lawrence Livermore National Laboratory.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2018.825.
Appendix A. Effect of mesh resolution
The results of grid convergence studies for the multiphase turbulence statistics are shown in figure 19. Four mesh resolutions are considered in the present study and the details are listed in table 3. It is seen that all four meshes used here capture the mean flow properties, including
$\overline{c}$
and
$\overline{u}$
, very well. When it comes to higher-order statistical terms, such as the Reynolds stress and the turbulence dissipation, then the M0 and M1 meshes are shown to be insufficient. It can be seen from figure 19(d) that when the cell size decreases from the M0 to M2 meshes, the magnitude of the turbulent dissipation increases significantly. In other words, the M0 and M1 meshes significantly underpredict the dissipation. The results for both Reynolds stress and turbulence dissipation for the M2 and M3 meshes generally agree very well, indicating that the M3 mesh is adequate to resolve the multiphase turbulence in the present problem. The agreement between the M2 and M3 results for the Reynolds stress component and the dissipation at
$y/H_{l}=1$
and 2 is not as good as in other regions. The discrepancy is more profound for the downstream location (
$x/H_{l}=10$
). This is mainly due to the fluctuations in the M3 results. The fluctuations in the M3 results can be further reduced by running the simulation for a longer time. Yet owing to the extreme computational cost, the longer run can only be left for future work. In spite of the noise, the current results are sufficient to provide reasonable estimate of Reynolds stress and turbulent dissipation in the two-phase mixing layer.
Grid convergence studies for multiphase turbulence statistics.
Appendix B. Effect of averaging time
The results of statistical convergence study for the multiphase turbulence statistics are shown in figure 20. Here, three different time durations are used for averaging in (3.1) where the averaging times T1, T2 and T3 are 43, 77 and 135
$H_{l}/U_{g}$
, respectively. The results clearly show that all three cases well capture the mean flow properties (
$\overline{c}$
and
$\overline{u}$
). Nevertheless, for higher-order statistics such as the Reynolds stress downstream and the turbulence dissipation, a longer averaging time is required. The simulation length T3 seems to yield converged results, although the turbulence dissipation at downstream location is still somewhat noisy. All the results for the M3 mesh presented in the results section are averaged over T3.
Statistical convergence studies for multiphase turbulence statistics.
Appendix C. Effect of the domain size
To examine the effect of the domain size, we have also considered two different domains that are larger than the present set-up: a wider domain (
$L_{x}/H_{l}=16$
,
$L_{y}/H_{l}=8$
and
$L_{z}/H_{l}=8$
) and a longer and higher domain (
$L_{x}/H_{l}=24$
,
$L_{y}/H_{l}=12$
and
$L_{z}/H_{l}=2$
). The results of the interfacial instability and the multiphase turbulence statistics are shown in figures 21 and 22, respectively. The tests are performed with a mesh resolution equivalent to M1 and the results in figure 22 are plotted with the same colour scale.
Effect of the domain width (in spanwise direction) on the spectra of (a) the enstrophy of gas in the domain and (b) the interfacial height at
$x/H_{l}=0.75$
. The domain width,
$L_{z}$
, is equal to
$2H_{l}$
and
$8H_{l}$
for the present and the wider domains, respectively. The results presented here are for the M1 mesh resolution, i.e.,
$H_{l}/\unicode[STIX]{x1D6E5}=64$
.
The development of the interfacial instability is characterized by the vorticity layer thickness (
$\unicode[STIX]{x1D6FF}_{g}$
) (Matas et al.
Reference Matas, Marty and Cartellier2011). The width of the present domain is significantly larger than
$\unicode[STIX]{x1D6FF}_{g}$
(
$H_{l}/\unicode[STIX]{x1D6FF}_{g}=16$
) and, thus, is sufficient to capture the dominant frequency arising from absolute instability. Figure 21 shows the spectra of the gas enstrophy and the interfacial height for the present and the wider domains, i.e.,
$L_{z}/H_{l}=2$
and 8. The simulation for the wider domain was run for a shorter time, so the spectra are more noisy; nevertheless, the results for the two domains generally agree well with each other. A dominant frequency is observed in both cases and is in reasonable agreement with the theoretical prediction. The dominant frequency seems to shift slightly to the right for the wider domain.
As the interfacial wave grows as it propagates downstream, transverse instability develops and the wave becomes fully 3D (Zandian et al.
Reference Zandian, Sirignano and Hussain2018). Then the domain constraint in the transverse direction will influence the transverse instability (since the long-wavelength instability will not be resolved) and later on wave breakup. From figure 22(a,b) it can be observed that although the results of
$\overline{u}$
and
$\unicode[STIX]{x1D716}$
for the two cases generally agree well, there exists a discrepancy in the contour line
$\overline{c}$
, showing that the unbroken liquid stream is shorter for the case with the wider domain. This indicates that the downstream dynamics of the two-phase mixing layer would require a wider domain to avoid any influence of boundary conditions but this is left for future research.
To avoid generating a recirculation above the parallel streams, we apply a Neumann boundary condition for the velocity on the top boundary to allow fluid to freely enter or leave the domain. Accordingly a velocity Dirichlet boundary condition is applied to the right boundary with the velocity profile specified as (2.7) and (2.10). The outflow velocity profile given is to mimic the mean flow near the outlet and doest not represent the exact time-dependent outflow condition. The results of the mean liquid volume fraction (the black lines indicated
$\bar{c}=0.5$
), the mean velocity and the turbulence dissipation for the present domain (
$L_{x}/H_{l}=16$
,
$L_{y}/H_{l}=8$
and
$L_{z}/H_{l}=2$
) are compared with those for a longer and higher domain (
$L_{x}/H_{l}=24$
,
$L_{y}/H_{l}=12$
and
$L_{z}/H_{l}=2$
). It can be seen that the results of the present domain in general agree well with those in the larger domain, except a small region near the outlet.
Therefore, the important observations made in the results section are thus confirmed to be not influenced by the domain size and the applied boundary conditions.
Comparison of results for different domain size: (a) the present domain (
$L_{x}/H_{l}=16$
,
$L_{y}/H_{l}=8$
and
$L_{z}/H_{l}=2$
), (b) the wider domain (
$L_{x}/H_{l}=16$
,
$L_{y}/H_{l}=8$
and
$L_{z}/H_{l}=8$
) and (c) the longer and higher domain (
$L_{x}/H_{l}=24$
,
$L_{y}/H_{l}=12$
and
$L_{z}/H_{l}=2$
). The results presented here are for the M1 mesh resolution, i.e.,
$H_{l}/\unicode[STIX]{x1D6E5}=64$
. The left and right columns are the mean velocity and the turbulent dissipation, respectively. The black curve corresponds to
$\overline{c}=0.5$
. The box with green dashed lines in (c) indicate the present domain size. The orange rectangle near the inlet represents the separator plate.
