# Impact of inlet gas turbulence on the formation, development, and breakup of
interfacial waves in a two-phase mixing layer

D. Jiang  Y. Ling1<EMAIL_ADDRESS>1Department of Mechanical
Engineering, Baylor University, Waco, TX 76798, USA

###### Abstract

Understanding the development and breakup of interfacial waves in a two-phase
mixing layer between the gas and liquid streams is paramount to atomization.
Due to the velocity difference between the two streams, the shear on the
interface triggers a longitudinal instability, which develops to interfacial
waves that propagate downstream. As the interfacial waves grow spatially,
transverse modulations arise, turning the interfacial waves from quasi-2D to
fully 3D. The inlet gas turbulence intensity has a strong impact on the
interfacial instability. Therefore, parametric direct numerical simulations
are performed in the present study to systematically investigate the effect of
the inlet gas turbulence on the formation, development, and breakup of the
interfacial waves. The open-source multiphase flow solver, PARIS, is used for
the simulations and the mass-momentum consistent volume-of-fluid method is
used to capture the sharp gas-liquid interfaces. Two computational domain
widths are considered and the wide domain will allow a detailed study of the
transverse development of the interfacial waves. The dominant frequency and
spatial growth rate of the longitudinal instability are found to increase with
the inlet gas turbulence intensity. The dominant transverse wavenumber,
determined by the Rayleigh-Taylor instability, scales with the longitudinal
frequency, so it also increases with the inlet gas turbulence intensity. The
holes formed in the liquid sheet is important to the disintegration of the
interfacial waves. The holes formation is influenced by the inlet gas
turbulence. As a result, the sheet breakup dynamics and the statistics of the
droplets formed also change accordingly.

###### keywords:

Interfacial waves; Atomization; Direct numerical simulation; Two-phase mixing
layer

##  1 Introduction

When a fast gas stream interacts with a parallel co-flowing liquid stream of
lower velocity, the gas-liquid interface is unstable and interfacial waves
will form and develop. As the interfacial waves are advected downstream, the
waves grow in amplitude, deform, and eventually break into small droplets. The
droplets are dispersed by the gas stream, forming a two-phase mixing layer
between the gas and liquid streams. The two-phase mixing layer is at the heart
of numerous twin-fluid atomization processes, such as air-assisted and air-
blast atomizations (Lefebvre, 1980, 1988).

The formation and early development of the interfacial waves are mainly
controlled by the longitudinal shear-induced Kelvin-Helmholtz-like
instability. The most unstable mode in the longitudinal instability dictates
the frequency and the wavelengths of the interfacial waves. Linear stability
analysis has been carried out to predict the most-unstable longitudinal wave
frequency $f$ and the wavelength $\lambda$. While inviscid analysis (Raynal,
1997; Marmottant & Villermaux, 2004; Matas et al., 2011) yielded reasonable
predictions for the frequency but under-predicted the spatial growth rate,
viscous temporal analysis (Boeck & Zaleski, 2005) well predicted the growth
rate but overestimated frequency. It was later shown by Otto et al. (2013) and
Fuster et al. (2013) that viscous spatial-temporal stability analysis is
required to well predict both the frequency and the growth rate. Due to the
velocity deficit induced by the wake of the separator plate, the instability
can transfer from convective to absolute instability regimes. When the
instability is absolute, a dominant unstable mode will arise. Matas (2015)
further showed that the confinement effect induced by the stream thickness
needs to be taken into account to yield a good prediction of the most unstable
mode. When the longitudinal wavelength $\lambda$ is significantly lower than
the stream thickness, then it scales with the gas stream boundary layer
thickness $\delta$ at the inlet and the absolute instability is mainly
controlled by the surface tension mechanism. In contrast, if the most-unstable
wavelength is comparable or higher than the stream thickness, then the
inviscid confinement absolute instability will become dominant (Matas et al.,
2018). In such a case, the wave propagation speed was found to agree well with
the Dimotakis speed (Dimotakis, 1986),
$U_{D}=(\sqrt{\rho_{l}}U_{l}+\sqrt{\rho_{g}}U_{g})/(\sqrt{\rho_{l}}+\sqrt{\rho_{g}})$,
where $\rho_{l},U_{l}$ and $\rho_{g},U_{g}$ represent the densities and
velocities for the gas and liquid streams, respectively. The most unstable
wavelength and frequency can thus be related to each other by the Dimotakis
speed as $\lambda=U_{D}/f$.

Conventionally, the stability analysis and numerical studies of the
longitudinal shear-induced instability assume both the gas and liquid streams
are laminar when they meet (Boeck & Zaleski, 2005; Otto et al., 2013; Fuster
et al., 2013; Agbaglah et al., 2017; Ling et al., 2017, 2019). Nevertheless,
turbulent fluctuations may exist in the gas stream in experiment due to the
high Reynolds numbers. As indicated by Matas et al. (2015) that, this may be a
potential reason for the discrepancies between different experiments. The
impact of the gas inlet turbulence on the longitudinal instability has been
recently investigated through experiment by Matas et al. (2015) and later
using direct numerical simulation by Jiang & Ling (2020). Both the frequency
and the spatial growth rate were observed to increase with the inlet gas
turbulence intensity $I$, when $I$ is over a threshold (Matas et al., 2015;
Jiang & Ling, 2020). Attempts have been made to incorporate the effect of
inlet gas turbulence in the linear viscous spatial-temporal stability analysis
based on turbulent viscosity models. The modified stability analysis
reasonably capture the trend, i.e., the frequency increases with $I$, but
underestimates the values (Jiang & Ling, 2020). An accurate linear stability
theory that captures the most unstable modes for a turbulent gas stream
remains to be developed.

Transverse modulations on the longitudinal waves develop when the waves grow
and propagate downstream. The Rayleigh-Taylor (RT) instability has been shown
to be the primary driving mechanism for the transverse instability in a
cylindrical coaxial configuration (Varga et al., 2003; Marmottant &
Villermaux, 2004). The azimuthal wavelength was estimated based on inviscid RT
instability on a planar surface with undulations from the longitudinal
instability. When the interface accelerates toward liquid or decelerates
toward gas, the interface is unstable and the most-unstable wavelength is
dictated by the surface tension and the interface acceleration. In the work of
Marmottant & Villermaux (2004), the maximum acceleration is estimated based on
experimental correlation and is a function of $\text{We}_{\delta}$, a Weber
number based on $\delta$. Numerical studies by Jarrahbashi & Sirignano (2014)
showed that both the RT instability (baroclinic effects) and strain-vorticity
interaction contribute to the transverse instability development. The latter
is more important when the gas-to-liquid density ratio is high. The sequential
development of the transverse variation of the interfacial wave is influenced
by the interaction between the interfacial waves and the gas stream. As the
interfacial waves grow and invade into the gas stream, the accelerated gas
flow above the wave crest induces a Bernoulli depression, which enhances the
growth of the wave (Hoepffner et al., 2011). Since the transverse instability
is closely tied to the longitudinal counterpart and the latter is in turn
influenced by the inlet gas turbulence, it is expected that the inlet gas
turbulence will also play a significant role in the transverse instability and
the subsequent transverse development of the interfacial waves. A detailed
analysis on the effect of inlet gas turbulence on transverse instability
features such as the dominant transverse wavenumber remains absent in the
literature.

There exist multiple pathways for the 3D interfacial waves to break into
filaments and droplets: the fingering and the hole-in-sheet modes.
Visualization of these two breakup modes have been presented in high-fidelity
numerical simulations (Jarrahbashi et al., 2016; Ling et al., 2017; Zandian et
al., 2018). The fingering modes typically occur when the disintegration of the
interfacial wave is relatively mild. In such a case, the Rayleigh-Plateau (RP)
instability gets a chance to develop at the Taylor-Culick rims on the edge of
the liquid sheets extended from the waves (Roisman, 2010; Agbaglah et al.,
2013). The RP instability results in liquid fingers which are approximately
aligned with the streamwise direction. The number of fingers formed is related
to the dominant transverse wavenumber (Marmottant & Villermaux, 2004). The
liquid fingers will continue to break into droplets.

The formation of droplets due to pinching of a filament is by itself a
complicated subject (Eggers, 1993; Ambravaneswaran et al., 2002; Castrejón-
Pita et al., 2015; Zhang et al., 2019). In general, the primary droplets
formed are similar to the local diameter of the filament, while the secondary
satellite droplets can be much smaller. Since the liquid fingers typically
exhibit irregular shapes, the breakup dynamics is more complex than the
classic Rayleigh breakup of a liquid cylinder and generally leads to a
distribution of droplets of different sizes (Villermaux et al., 2004; Ling et
al., 2017). The size distribution of droplets formed is essential to spray
applications and different distribution models have been proposed, , the log-
normal, exponential, Poisson, Weibull-Rosin-Rammler, Pareto, and gamma
distributions (Villermaux et al., 2004; Herrmann, 2011; Marty, 2015; Ling et
al., 2017; Kooij et al., 2018; Balachandar et al., 2020). These distribution
functions agree with different sets of experimental or simulation data in some
extent. Whether there exists a universal distribution of droplet size remains
an unresolved question.

When long liquid sheets extend from the interfacial wave crest and has a
strong interaction with the gas stream, the disintegration of the interfacial
wave is generally more violent and tends to follow the hole-in-sheet mode
(Ling et al., 2017). The hole expansion speed follows the Taylor-Culick
velocity (Opfer et al., 2014; Marston et al., 2016; Ling et al., 2017;
Agbaglah, 2021). The holes grow and merge, and eventually lead to a violent
rupture of the sheet, producing separate ligament, droplets of different sizes
and orientations, and fingers that remain attached to the liquid sheet. The
formation of a hole in a liquid sheet is due to the pinching of the two
surfaces of the liquid sheet, similar to the pinching of a filament in drop
formation. The pinching of surfaces is induced by the disjoining pressure when
the distance between the two surfaces is sufficiently small (O(10 nm)). In
numerical simulations, the cell size is usually much larger, so the disjoining
pressure is typically not included in the physical model. As a result, the
minimum cell size serves as the numerical cut-off length scale to pinch a
liquid sheet and to form holes. The absence of disjoining pressure seems to
have little effect on the surface pinching, since the process for low-
viscosity liquids is mainly dictated by the fluid inertia, similar to the
droplet formation due to the pinching of filament (Zhang et al., 2019).
Nevertheless, a careful grid-refinement study is still required to verify the
simulation results, in particular for the statistics of the droplets formed
(Ling et al., 2017). Former studies on the interfacial waves breakup assume
that the inlet gas stream is laminar. The effect of the inlet gas turbulence
on the sheet breakup dynamics and the droplets statistics remains unclear.

The goal of the present study is to investigate the fate of the interfacial
waves in a two-phase mixing layer between parallel liquid and gas streams
through direct numerical simulations. As an extension of our former studies
(Ling et al., 2017, 2019; Jiang & Ling, 2020), the present study is focused on
the transverse interfacial instability and the impact of inlet gas turbulence
on the formation, development, and breakup of the three-dimensional
interfacial waves. To allow a detailed investigation of the transverse
features of interfacial waves, we have used a computational domain that is
three times as wide as that in our former studies. The key questions we aim to
address include:

  1. 1.

What are the physical mechanisms that drive the transverse development of the
interfacial waves? Can one predict the transverse wavenumber based on the
stability theory?

  2. 2.

How will the inlet gas turbulence influence the longitudinal and transverse
interfacial instability and the development of the 3D interfacial waves?

  3. 3.

What are the pathways for the interfacial waves to disintegrate into filaments
and droplets? What is the effect of the inlet gas turbulence on the
interfacial wave breakup dynamics and the statistics of the droplets formed?

The rest of the paper will be organized as follows. The problem description
and the simulation approaches, including the governing equations, the
numerical methods, and the simulation setup, will be presented in section 2.
The simulation results will be shown in section 3. The longitudinal and
transverse instabilities of the interfacial waves, the development and breakup
of the interfacial waves, and the droplet statistics will be discussed in
sequence. The key conclusions of the present study will be summarized in
section 4.

##  2 Simulation Methods

###  2.1 Problem description

The two-phase mixing layer to be considered in the present study is
illustrated in figure 1. Two parallel planar gas and liquid streams enter the
domain from the left. Two thin solid plates are placed near the inlet to
separate the two streams, mimicking the injector nozzle. The inclusion of the
separator plate has been shown to be important to the interfacial instability
(Otto et al., 2013). The two streams meet at the end of the lower separator
plate. The liquid inflow is laminar, while turbulent velocity fluctuations of
different intensity levels are present at the gas inlet. The mean gas velocity
is significantly higher than that for the liquid. As a result, the gas-liquid
interface is unstable, and wavy structures develop on the interface and
propagate downstream. For convenience of discussions, we refer to the $x$,
$y$, and $z$ directions as the longitudinal, vertical, and transverse
directions, respectively. The interfacial waves, when they are just formed,
are approximately 2D and longitudinal. Yet as they grow and propagate
downstream, transverse modulations arise and the waves evolve to be fully 3D.
The interfacial waves will interact with the gas stream as the amplitudes
become finite. The eventual breakups of the interfacial waves produce small
ligaments and droplets that are mixed with the gas stream, forming a two-phase
mixing layer. After an initial transition period for the two streams to
progressively enter the domain, the wave formation and the two-phase turbulent
flow reaches a statistically stationary state (Ling et al., 2017, 2019).

Figure 1: Simulation setup for the interfacial waves between parallel liquid
and gas streams.

###  2.2 Governing equations

The two-phase interfacial flows are governed by the incompressible Naviers-
Stokes equation with surface tension. The one-fluid approach is employed,
where the gas and liquid phases are treated as one fluid with material
properties change abruptly across the interface. The momentum and continuity
equations are expressed as

$\displaystyle\rho\left(\frac{\partial u_{i}}{\partial t}+u_{i}\frac{\partial
u_{j}}{\partial x_{j}}\right)$ $\displaystyle=-\frac{\partial p}{\partial
x_{i}}+\frac{\partial}{\partial x_{j}}\left[\mu\left(\frac{\partial
u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial
x_{i}}\right)\right]+\sigma\kappa\delta_{s}n_{i}\,,$ (1)
$\displaystyle\frac{\partial u_{i}}{\partial x_{i}}$ $\displaystyle=0\,.$ (2)

where $\rho,u_{i},p,\mu$ represent density, velocity, pressure and viscosity,
respectively. The last term on the right hand side of the momentum equation
represents the surface-tension, which is a singular force localized on the
sharp interface using the Dirac distribution function $\delta_{s}$. The
surface tension coefficient $\sigma$ is taken to be constant, and $\kappa$ and
$n_{i}$ represent the curvature and normal vector of the interface.

The gas and liquid phases are distinguished by a characteristic function
$\chi$, while $\chi=1$ and $0$ represent the liquid and gas phases,
respectively. The interface evolution can be captured by solving the advection
equation of $\chi$.

$\displaystyle\frac{\partial\chi}{\partial
t}+u_{i}\frac{\partial\chi}{\partial x_{i}}$ $\displaystyle=0\,.$ (3)

The mean value of $\chi$ in a computational cell is defined as

$\displaystyle f=\frac{1}{\Delta\Omega}\int_{\Omega}\chi dV\,.$ (4)

which also represents the volume fraction of liquid ($\chi=1$) in a cell.
Correspondingly, the gas volume fraction in a cell is $\hat{f}=1-f$. The fluid
properties in interfacial cells with $0<f<1$ are calculated based on
arithmetic mean

$\displaystyle\rho=f\rho_{l}+\hat{f}\rho_{g}\,,\quad\mu=f\mu_{l}+\hat{f}\mu_{g}\,.$
(5)

where the subscripts $l$ and $g$ represent variables corresponding to the
liquid and gas phases, respectively.

###  2.3 Numerical methods

The governing equations are solved by the finite volume method on a staggered
grid. The advection equation, Eq. (3), is solved using a geometric volume-of-
fluid (VOF) method. The interface normal is computed following the mixed
Young’s-centred method of Aulisa et al. (2007). The Lagrangian-explicit scheme
of Li (1995) is used for the VOF advection (Scardovelli & Zaleski, 2003). The
convection term in the momentum equation, Eq.(1), is discretized consistently
with the VOF method (Arrufat et al., 2020), and this mass-momentum consistence
has been shown to be crucial in capturing interfacial dynamics when large
velocity and density contrasts are present at the interface (Rudman, 1998;
Ling et al., 2017; Vaudor et al., 2017; Zhang et al., 2020; Arrufat et al.,
2020). The incompressibility condition is incorporated using the projection
method (Chorin, 1968). The pressure Poisson equation is solved using PFMG
multigrid solver in the HYPRE library. The viscous term is discretized
explicitly using the second-order centered difference scheme. The interface
curvature is calculated using the height-function method of Popinet (2009) and
the balanced continuous-surface-force method is used to discretize the surface
tension term (Renardy & Renardy, 2002; Francois et al., 2006; Popinet, 2009).
The time integration is done by a second-order predictor-corrector method. To
capture the dynamics of under-resolved droplets less erroneously than by just
quasi-fragment VOF patches, droplets of size smaller than about two cells are
converted into Lagrangian point-particles and are traced under the one-way
coupling approximation, following the approach of Ling et al. (2015).

The aforementioned numerical methods have been implemented in the open-source
solver, _PARIS-Simulator_. Detailed implementations and validation of the code
can be found in previous studies (Tryggvason et al., 2011; Ling et al., 2015,
2017, 2019; Arrufat et al., 2020; Aniszewski et al., 2021).

###  2.4 Simulation setup

The computational domain is a cuboid. The dimensions in $x$ and $y$ directions
for the two domains are the same, i.e., $L_{x}=16H$ and $L_{y}=8H$, where $H$
is the height of the liquid stream at the inlet. Periodic boundary conditions
are applied to the front and back of the domain. A free boundary condition is
applied at the top, so gas is allowed to freely flow through the boundary. The
velocity outflow boundary condition is imposed at the right surface. The
bottom surface is treated as a slip wall to avoid the modeling challenge for
moving contact lines on a no-slip surface. The contact angle is specified as
90 degrees on the bottom surface. As a result, the bottom is equivalent to a
symmetric boundary and the present simulation setup can be viewed as a
symmetric model for the airblast atomization configurations with gas streams
on both sides of the liquid stream (Chaussonnet et al., 2020), which could
capture the interfacial stability when the symmetric mode is dominant. A full
simulation with both gas streams will be required to resolve the asymmetric
instability mode and breakup dynamics (Delon et al., 2018).

The length and height of the domain and boundary conditions used have been
examined in previous studies, which are shown to be sufficient to resolve
present problem (Ling et al., 2019). Two different domain widths, $L_{z}=2H$
and $6H$, are considered. Some of the simulation results for the narrow domain
have been shown in our previous study (Jiang & Ling, 2020). While the narrow
domain is useful in capturing the longitudinal wave formation, the wide domain
($L_{z}=6H$) is required to investigate the transverse instability and the
transverse development of the interfacial waves.

The gas stream has a similar thickness as the liquid stream, i.e.,
$H-\eta_{y}$, where $\eta_{y}$ is the thickness of the thin separator plate.
The liquid and gas properties are similar to those of water and pressurized
air, following the previous studies (Ling et al., 2017, 2019; Jiang & Ling,
2020), see table 1.

The mean flow at the inlet is horizontal, so the $y$\- and $z$-components of
the mean velocity are zero, i.e., $\overline{v}_{0}=\overline{w}_{0}=0$. The
$x$-component of the mean velocity at the inlet is expressed as

$\overline{u}_{0}(y)=\left\\{\begin{array}[]{ll}U_{l}\,\mathrm{erf}\left[\frac{H-y}{\delta}\right],&0\leq
y<H\ \text{(liquid stream)},\\\ 0,&H\leq y<H+\eta_{y}\ \text{(separator
plate)},\\\
U_{g}\,\mathrm{erf}\left[\frac{y-(H+\eta_{y})}{\delta}\right]\mathrm{erf}\left[\frac{2H-y}{\delta}\right],&H+\eta_{y}\leq
y<2H\ \text{(gas stream)},\\\ 0,&\mathrm{else}\
\text{(wall)}.\end{array}\right.$ (6)

The velocities in the gas and liquid stream are generally uniform and are
equal to $U_{l}$ and $U_{g}$ away from the separator plates, respectively, see
figure 2(a). The error function is used to model velocity profile in the
boundary layers near the separator plates. The parameter $\delta$
characterizes the boundary layer thickness, which is taken to be the same for
both the gas and liquid streams, i.e., $\delta=H/8$. The dimensions of the two
separator plates are the same, the thickness and length of which are
$\eta_{y}=H/32$ and $\eta_{x}=H/2$, respectively. The thickness $\eta_{y}$ is
chosen to be significantly smaller than $\delta$. Based on the former study of
Fuster et al. (2013), the specific value of $\eta_{y}$ has negligible effect
on the interfacial instability. The values for the key physical parameters for
the present problem are listed in table 1.

$\rho_{l}$ | $\rho_{g}$ | $\mu_{l}$ | $\mu_{g}$ | $\sigma$ | $U_{l}$ | $U_{g}$ | $H$ | $\delta$
---|---|---|---|---|---|---|---|---
(kg/m3) | (kg/m3) | (Pa s) | (Pa s) | (N/m) | (m/s) | (m/s) | (mm) | (mm)
1000 | 50 | $\mathrm{10^{-3}}$ | 5$\times 10^{-5}$ | 0.05 | 0.5 | 10 | 0.8 | 0.1
Table 1: Physical parameters.

Pseudo turbulent velocity fluctuations are superposed on the gas velocity at
the inlet (i.e., $H+\eta_{y}\leq y<2H$) as

$u_{0}=\overline{u}_{0}+u^{\prime}(y,z,t)\,,\quad
v_{0}=v^{\prime}(y,z,t)\,,\quad w_{0}=w^{\prime}(y,z,t)\,,$ (7)

where the turbulent fluctuations $u^{\prime},v^{\prime},w^{\prime}$ are
computed using the digital filter approach of Klein et al. (2003). The filter
length is $1/4H$. Implementation details and verification for the approach to
generate turbulent velocity fluctuations have been given in previous study
(Jiang & Ling, 2020) and thus are not repeated here.

###  2.5 Data collection and processing

Time and spatial averaging are performed to process the instantaneous 3D field
of simulation data. The time averaging operator for a variable $a$ is denoted
as

$\displaystyle\overline{a}=\frac{1}{t_{1}-t_{0}}\int_{t_{0}}^{t_{1}}a(t)\,dt\,$
(8)

where $t_{0}$ and $t_{1}$ represent the beginning and the end time of the
sampling. Sampling of data will not start until the simulation reaches the
statistically steady state at about $t^{*}=200$.

The spatial averaging in the transverse $z$-direction is defined as

$\displaystyle\langle a\rangle=\frac{1}{L_{z}}\int_{0}^{L_{z}}a(z)\,dz$ (9)

The operation of double temporal and spatial averaging is denoted as
$\langle\overline{a}\rangle$.

The height functions, which are computed for the interface curvature, will
also be used to evaluate the interface location. The interfacial height, i.e.,
the $y$-coordinate of the interfacial location, is denoted by $h(t,x,z)$. The
root-mean-square of the interfacial height fluctuations,
$\overline{h^{\prime}h^{\prime}}$, represents the thickness of the two-mixing
layer or the amplitude interfacial wave, where $h^{\prime}=h-\overline{h}$.

###  2.6 Key dimensionless parameters

With $H$ and $U_{g}$ as the scaling variables, the dimensionless time,
velocity, and length are defined as $t^{*}=tU_{g}/H$, $u^{*}=u/U_{g}$,
$x^{*}=x/H$, respectively. The key parameters listed in table 1 can be
converted to the dimensionless form, see table 2. In addition to $H$, the
boundary layer thickness $\delta$ can also be important to the interfacial
instability (Otto et al., 2013). As a result, two different Reynolds numbers,
i.e., $Re_{g,H}$ and $Re_{g,\delta}$, are defined correspondingly. Due to the
high $\mathrm{Re}_{g,H}$, the gas stream will be turbulent even if the gas
inflow is laminar. The Weber number based on $\delta$ characterize the effect
of surface tension on the interfacial instability development (Otto et al.,
2013). The gas-to-liquid dynamic pressure ($M$) is important in determining
the macro-scale features, such as the breakup length (Lasheras et al., 1998),
and also the regime for the shear-induced longitudinal instability (Otto et
al. (2013); Fuster et al. (2013)). According to the previous studies (Ling et
al., 2019; Jiang & Ling, 2020), the selected dimensionless parameters place
the shear longitudinal instability in the absolute regime. There are two
different mechanisms, namely surface-tension and confinement, that can drive a
transition from convective to absolute regimes. Based on the results to be
shown later, the absolute longitudinal instability in the present problem
seems to belong to the confinement category.

###  2.7 Simulation cases

Figure 2: Temporally and spatially (in $z$-direction) averaged profiles of
(a): velocity in streamwise direction (b): normal Reynolds stress and (c):
shear Reynolds stress.

A parametric study is carried out to systematically investigate the inlet gas
turbulence intensity $I$. The simulation cases are summarized in table 3. For
the narrow domain, five different $I$ are considered (cases N0 to N4), while
for the wide domain, only three different $I$ (cases W0 to W2) are simulated
due to the higher computational costs.

$M$ | $r$ | $m$ | $\mathrm{Re}_{g,\delta}$ | $\mathrm{We}_{g,\delta}$ | $\mathrm{Re}_{g,H}$ | $I$
---|---|---|---|---|---|---
$\rho_{g}U_{g}^{2}/(\rho_{l}U_{l}^{2})$ | $\rho_{l}/\rho_{g}$ | $\mu_{l}/\mu_{g}$ | $\rho_{g}U_{g}\delta/\mu_{g}$ | $\rho_{g}U_{g}^{2}\delta/\sigma$ | $\rho_{g}U_{g}H_{g}/\mu_{g}$ | $(\sqrt{\overline{u^{\prime}u^{\prime}}})_{e}/U_{g}$
20 | 20 | 20 | 1000 | 10 | $7750$ | 0, 0.06, 0.13
Table 2: Key dimensionless parameters.

Temporally and spatially averaged profiles of the streamwise velocity
($\langle\overline{u}\rangle/U_{g}$), normal and shear Reynolds stresses
($\langle\overline{u^{\prime}u^{\prime}}\rangle/U_{g}^{2}$ and
$\langle\overline{u^{\prime}v^{\prime}}\rangle/U_{g}^{2}$) at the end of the
separator plate ($x=\eta_{x}$) for all the cases are shown in figure 2. It can
be seen that the mean velocity profiles for all the cases are almost
identical. For the cases N0 and W0, the Reynolds stresses are negligibly
small. Magnitudes of the Reynolds stresses generally increase with the inlet
gas turbulence intensity. The profiles of the normal and shear stresses are
consistent with those for planar turbulent channel flows, indicating the
turbulent fluctuations generator at the inlet has been implemented properly
(Klein et al., 2003). The normal Reynolds stress
($\langle\overline{u^{\prime}u^{\prime}}\rangle/u_{g}^{2}$) increases with
$y^{*}$ from zero at the separator plate and reaches the maximum at around the
middle of the boundary layer. Then it decreases and approaches a constant
outside of the boundary layer, the square root of the normal Reynolds
stress,i.e.,
$I=(\langle\overline{u^{\prime}u^{\prime}}\rangle_{y^{*}=1.5})^{1/2}/U_{g}$,
is used to characterize the inlet gas turbulence intensity.

The computational domains are discretized using a fixed regular cubic grid.
The cell size $\Delta=6.25$ µm ($H/\Delta=128$) is used for all the cases. The
cell size has been verified to be adequate for good estimates of high-order
two-phase turbulence statistics, such as TKE dissipation (Ling et al., 2019).
The numbers of cells for the narrow and wide domains are about 0.5 and 1.6
billions, see table 3. All cases are run for a physical time to at least
$t^{*}=450$, which is sufficient to cover the formation of about 12 to 18
waves after the statistically steady state is reached. The present simulations
are performed on the Intel Xeon Platinum 8160 (Skylake) computing nodes on the
TACC Stampede2 machine. For the narrow domain cases (N0 to N4), 22 nodes are
used, while 64 nodes are used for the wide domain cases (W0 to W2). The total
computing time used for all the cases is about 280,000 node-hours (13,400,000
core-hours).

##  3 Results and Discussion

###  3.1 General behavior

The temporal evolutions of the interfacial waves and the velocity fields for
the cases W0 and W2 are shown in figure 3. The cases W0 and W2 represent the
two distinct gas inflow conditions: laminar ($I=0$) and highly turbulent
($I=0.13$). The figures 3(a) and 3(b) in the left column show the interfacial
waves from a 3D view, with the streamwise $u$-velocity on the background. The
figures 3(c) and 3(d) in the right column are sequential snapshots of the gas-
liquid interfaces colored by the $u$-velocity from the top view.

The results for the case W0 shown in figure 3(c) represent a typical process
of interfacial wave formation and development for a planar two-phase mixing
layer (Matas et al., 2011; Agbaglah et al., 2017; Ling et al., 2017; Zandian
et al., 2018; Ling et al., 2019). When the two streams meet at the end of the
separator plate, the longitudinal instability is triggered due to the shear at
the interface. The instability develops to a 2D longitudinal interfacial wave
which propagates downstream, the wave crest of which is approximately a
straight line from the top view. Due to the transverse instability, the height
of the wave crest varies over the transverse direction. When the wave
amplitude grows, the wave interacts with the fast gas stream and develop into
a 3D wave, exhibiting a lobe shape. The liquid lobe later extends to form a
liquid sheet which bends and folds. The thickness of the sheet reduces rapidly
and unevenly due to the complex sheet deformation. Multiple holes are formed
near the edge of the sheet. The expansion and merging of the holes
disintegrate the liquid sheet violently. Multiple breakup events occur (see
snapshots between $t^{*}=335$ and 345 in figure 3(c)), until the sheet
completely breaks into droplets.

Case | $I$ |  $\Delta$ (µm) | $H/\Delta$ | $cells\\#$ | $Cores\\#$
---|---|---|---|---|---
W0 | $0.0$ | $6.25$ | $128$ | $1.6\times 10^{9}$ | $3072$
W1 | $0.06$ | $1.6\times 10^{9}$ | $3072$
W2 | $0.13$ | $1.6\times 10^{9}$ | $3072$
N0 | $0.0$ | $0.5\times 10^{9}$ | $1024$
N1 | $0.02$ | $0.5\times 10^{9}$ | $1024$
N2 | $0.039$ | $0.5\times 10^{9}$ | $1024$
N3 | $0.056$ | $0.5\times 10^{9}$ | $1024$
N4 | $0.071$ | $0.5\times 10^{9}$ | $1024$
Table 3: Summary of simulation runs. Case names starting with $W$ and $N$
represent the wide and narrow domains, respectively.  Figure 3: Snapshots for
the two-phase mixing layer for (a) the case W0 with a laminar gas inflow and
(b) the case W2 with a turbulent gas inlet. The color on the background in (a)
and (b) represents the velocity magnitude. Temporal development of the
interfacial waves for the cases W0 and W2 are shown in (c) and (d),
respectively, where the color on the interfaces represents the streamwise
velocity. The frames in (c) and (d) are moving with the waves to keep the
waves located at the center of the window. Supplementary movies are available
online.

By adding small-amplitude turbulent velocity fluctuations at the gas inlet
(case W2), several important differences are observed. First of all, the
longitudinal instability grows faster. As a result, the wave amplitude near
the inlet is much higher than that for the laminar gas inflow. Second, the
liquid sheet extended from the wave are much shorter, and experiences a weaker
interaction with the gas stream. Third, the wavenumber of the transverse
modulation of the interfacial wave increases. While there are about one to two
transverse waves shown in figure 3(c) (see $t^{*}=330$), about three to four
waves are observed in figure 3(d) (see $t^{*}=350$). Due to the increase of
wavenumber, the width of the lobes formed is reduced and the shape of the rim
formed at the edge of the liquid sheet is less regular. Fourth, the breakup
dynamics of the interfacial waves also changes significantly. Fewer holes are
seen for the case W2. The disintegration of the wave takes a different path,
i.e., forming fingers at the rim. Finally, the different breakup dynamic
impacts the statistics of the droplets generated. Significantly fewer droplets
are formed. Detailed quantitative analysis for these effects will be presented
in the following sections.

###  3.2 Longitudinal instability and wave formation

####  3.2.1 Absolute instability

The shear between the gas and liquid streams triggers a Kelvin-Helmholtz (KH)
longitudinal instability, which is the onset of the formation of interfacial
waves. The temporal evolutions of the interfacial heights are measured near
the inlet at $x^{*}=0.625$ and eleven evenly spaced locations along the
transverse direction. The temporal evolutions of the spatially-averaged
$h^{*}=h/H$ are shown in figure 4(a) for different cases of the wide and
narrow domains. At this $x$ location, the wave amplitude is small (less than
4% of $H$ for all cases), so the waves remain in the linear regime. Fourier
transform is used to cast the temporal signals to the frequency spectra, see
figure 4(b). The spectra obtained at eleven different transverse locations are
averaged. The amplitude at the dominant frequency is the maximum in the
spectra for all cases except W0, for which the dominant frequency is a local
maximum. The maximum amplitude for the case W0 is located at a very low
frequency that corresponds to the total simulation time.

Figure 4: (a) Temporal evolution of the transversely averaged interfacial
heights at $x^{*}=0.625$ for different cases. (b) Frequency spectra that
correspond to temporal signal in (a). (c) Variation of the normalized dominant
frequency for the longitudinal instability, $\omega_{L}^{*}/\omega_{L,0}^{*}$
as a function of $I$ for different cases, compared with the experimental data
of Matas et al. (2015).

The appearance of the dominant mode affirms that the cases studied are in the
absolute instability regime. The frequency of the longitudinal wave formation
is dictated by the most-unstable mode of the longitudinal instability. As
suggested by Matas (2015) and Matas et al. (2018), there are two different
mechanisms that contribute to the convective-to-absolute instability
transition: the surface tension and the confinement. Both types of absolute
instabilities have been identified in the spatial-temporal stability analysis,
as the pinching between the shear branch and the other branch controlled by
surface tension (Otto et al., 2013) or by confinement (Juniper, 2006; Healey,
2007). For both types of absolute instabilities, the imaginary part of the
frequency is positive at the pinch point. The wavenumber at the pinch point
for the confinement type is generally low, so that the most-unstable
wavelength is typically larger than the layer thickness $H$. In contrast, for
the surface-tension type, the wavenumber for the dominant mode is much higher
and the corresponding wavelength is smaller than $H$. Another fundamental
difference between the two is that, the phase velocity at the pinch point for
the confinement-triggered absolute instability follows the Dimotakis speed
$U_{D}$, while the surface-tension absolute instabilities typically exhibit
much smaller phase velocities.

The simulations for the cases with laminar gas inflows (N0 and W0) yield
similar results for the dominant frequency and wavelength,
$\omega_{L}^{*}=0.31$ and $\lambda^{*}=4.5$. Stability analysis without
confinement for the same case has been performed in previous studies (Ling et
al., 2019; Jiang & Ling, 2020), which represent the absolute instability
triggered by the surface tension mechanism. The stability analysis slightly
overestimates the dominant frequency $\omega^{*}_{\mathrm{st}}=0.40$, while
the predicted wavelength $\lambda^{*}_{\mathrm{st}}=1.6$ is significantly
lower than the simulation result, namely the surface-tension absolute
instability over-predicted the wavenumber. Furthermore, it is found in the
simulation results that, the wave propagation speed $U_{w}^{*}=0.22$, which
agrees well with the Dimotakis speed and is much higher than the wave speed
predicted by the surface-tension absolute instability,
$U^{*}_{w,\mathrm{st}}=0.10$. These discrepancies, together with the fact that
$\lambda^{*}>1$, seem to indicate the absolute shear instability observed here
belongs to the inviscid confinement mechanism, instead of the surface-tension
mechanism. Nevertheless, to fully confirm the nature of absolute instability,
stability analysis similar to Matas (2015) needs to be performed, which will
be relegated to our future work.

####  3.2.2 Dominant frequency

As shown in figure 4(b), when the gas inlet turbulence intensity $I$
increases, the dominant frequencies shift to the right for both narrow and
wide domains, and the amplitude also increases correspondingly. The normalized
frequencies for different cases, $\omega_{L}^{*}/\omega_{L,0}^{*}$, where
$\omega_{L,0}^{*}$ represents the frequency for $I=0$, are plotted in figure
4(c). The experimental data of Matas et al. (2015) are also shown for
comparison. The experiments covered a wide range of injection conditions for
different gas and liquid velocities, but the normalized frequencies
approximately collapse. Though the simulation conditions are not identical to
those for the experiments, the simulation results for the normalized frequency
agree reasonably well with the experimental results. In particular, the
general variation trend of $\omega_{L}^{*}/\omega_{L,0}^{*}$ over $I$ is well
captured, i.e., $\omega_{L}^{*}/\omega_{L,0}^{*}$ varies little for $I\lesssim
0.02$ and increases with $I$ for $I\gtrsim 0.02$. The discrepancy between the
simulation and experiment results may be due to the different density and
viscosity ratios used in the present simulation.

The increase of $\omega_{L}^{*}$ over $I$ can be explained by the viscous
spatial-temporal stability analysis with eddy viscosity model (O’Naraigh et
al., 2013; Matas et al., 2015; Jiang & Ling, 2020). The effective gas
viscosity is the sum of the molecular and eddy viscosities, which increases
with $I$ since the turbulent eddy viscosity increases with $I$. Furthermore,
the Orr-Sommerfeld system indicated that the dominant frequency increases with
the effective gas viscosity. When $I$ is small, the eddy viscosity is smaller
than the molecular counterpart and that is why $\omega_{L}^{*}$ varies little
for $I\lesssim 0.02$. More discussions on the stability analysis can be found
in our previous study (Jiang & Ling, 2020).

For similar $I$, $\omega_{L}^{*}/\omega_{L,0}^{*}$ for the wide domain is
smaller than that for the narrow domain. For example, the values of $I$ for
the cases N3 and W1 are similar, $I_{N3}=0.056$ and $I_{W1}=0.06$,
respectively, but $\omega_{L,W1}^{*}=0.38$ is about 9% lower than
$\omega_{L,N3}^{*}$. This seems to indicate that, the longitudinal instability
is more sensitive to domain width when the gas inflow is turbulent. Since the
turbulent gas flow is 3D in nature, the small domain width may have
constrained the development of the turbulence and the turbulence-interface
interaction.

####  3.2.3 Spatial growth

The spatial development of the longitudinal instability determines the
longitudinal growth of the interfacial wave amplitude. Since the interface
motion will introduce a fluctuation in the liquid volume fraction
$c^{\prime}=c-\langle\overline{c}\rangle$, the mean square of which, i.e.,
$\langle\overline{c^{\prime}c^{\prime}}\rangle$, is employed to measure the
spatial longitudinal wave amplitude, see figure 5. For a given $x^{*}$, the
vertical distance between the two contour lines of
$\langle\overline{c^{\prime}c^{\prime}}\rangle=0.02$ is defined as the
longitudinal wave amplitude $\xi_{L}^{*}$.

Figure 5: Spatial growth of the longitudinal wave amplitude (characterized by
the mean square of liquid volume fraction fluctuations
$\langle\overline{c^{\prime}c^{\prime}}\rangle$, for the cases (a) W0, (b) W1
and (c) W2. The contour line
$\langle\overline{c^{\prime}c^{\prime}}\rangle=0.02$ is used to measured the
wave amplitude $\xi_{L}^{*}$.

It can be seen that, $\xi_{L}^{*}$ grows rapidly with $x^{*}$ near the nozzle
exit and becomes more gradually downstream. From the log-linear plot shown in
figure 6(a), the spatial growth of $\xi_{L}^{*}$ near the nozzle exit is
approximately exponential, see , $0.5\lesssim x^{*}\lesssim 0.8$ for the case
W2. Nevertheless, it should be reminded that, since the longitudinal
instability here is absolute, nonlinearity will influence the properties of
the most-unstable mode and the spatial growth does not follow an exponential
function as in convective instability (Otto et al., 2013; Matas, 2015).

As the inlet gas turbulence intensity $I$ increases, the spatial growth of
$\xi_{L}^{*}$, particularly near the nozzle exit, becomes faster. To
characterize the effect of $I$, the spatial growth rate, $\alpha_{L}^{*}$,
namely the slope of the curve $\log{\xi_{L}^{*}}$-$x^{*}$, is measured at
$\log(\xi_{L}^{*})\approx-2$ for all cases. The variation of $\alpha_{L}^{*}$
over $I$ is shown in figure 6(b). Similar to the dominant frequency,
$\alpha_{L}^{*}$ also increases with $I$. If a different threshold value for
$\langle\overline{c^{\prime}c^{\prime}}\rangle$ other than 0.02 is to be used,
the values of $\xi_{L}^{*}$ would change, but those for $\alpha_{L}^{*}$ will
not be affected. For the cases W0 and N0, the variation $\xi_{L}^{*}$ over
$x^{*}$ exhibits fluctuations, which makes the measurement of $\alpha_{L}^{*}$
somewhat sensitive. As will be discussed later in section 3.3.4, these
fluctuations are induced by the spatial averaging over the transverse
direction. To reduce the influence of the fluctuations, the measurement of
$\alpha_{L}^{*}$ for the cases N0 and W0 is made for a wider range of $x^{*}$,
as indicated in figure 6(a).

Figure 6: (a) The growth of the longitudinal wave amplitude $\xi_{L}^{*}$ for
different cases. The dotted lines indicate the spatial growth rate
$\alpha^{*}_{L}$ for different cases at $\log(\xi_{L}^{*})\approx-2$. (b)
Variation of the normalized spatial growth rate $\alpha_{L}/\alpha_{L,0}$ as a
function of $I$.

It can be observed from 6(a) that the values of $\xi_{L}^{*}$ at the end of
the separator plate, $x^{*}=\eta_{x}^{*}=0.5$, are almost zero for the narrow
domain cases, but are finite for the wide domain cases. A snapshot of the
interface near the inlet for the case W0 is shown in figure 7(a). The contact
line where the three phases meet can be seen. A closeup in figure 7(c) shows
that the contact line varies in time and also in the transverse direction.
This indicates that the contact line moves on the separator front wall (facing
the streamwise direction). Accurately capturing the contact-line dynamics is
out of the scope of the present study. The boundary conditions on the solid
separator plate are no-slip wall for velocity and symmetric for the liquid
volume fraction $c$ (yielding a 90 degree contact angle). The small-amplitude
spatial variation of the contact line is due to a numerical slip, with half of
the cell size as the slip length (Snoeijer & Andreotti, 2013). For the case
N0, due to the constraint of the smaller domain width, the contact line
remains to be a straight line pinned on the lower edge of the separator front
wall, see figure 7(b). Therefore, when we compute the longitudinal wave
amplitude by averaging $\overline{c^{\prime}c^{\prime}}$ in the transverse
direction, the amplitude at $x^{*}=\eta_{x}^{*}$ for the case N0 is identical
to zero, while that for the case W0 is finite.

Figure 7: (a) Closeup of the interface near the inlet and interfaces at
$x^{*}=\eta_{x}^{*}=0.5$ for different times for the cases (a) N0 and (b) W0.
The interfaces for the case N0 collapse to a straight line.

Finally, it is worth mentioning that linear stability analysis has been
performed in previous studies to identify the most unstable mode in the
present configuration. Among the many attempts, the viscous spatial-temporal
analysis has been shown to be quite successful (Otto et al., 2013; Matas,
2015). The perturbation is introduced in the form of a 2D normal mode. The
resulting Orr-Sommerfeld equations are solved to obtain the most-unstable
frequency and growth rate, which have been shown to agree well with the
experimental and DNS results when the inlet gas is laminar (Fuster et al.,
2013). In order to predict the most unstable mode when turbulence is present
at the gas inlet, additional modeling efforts are required. Attempts have been
made by Matas et al. (2015) and Jiang & Ling (2020) by incorporating the
effect of inlet turbulence using the turbulent viscosity model. As the
turbulent viscosity increases with $I$, the Orr-Sommerfeld equations can
reproduce the trends that the frequency increases with $I$. However, the
modified stability model underestimates the values. The present simulation
results indicate that the interface exhibits 3D features right at the end of
the separator plate, the assumption of 2D mode in conventional stability
analysis may need to be relaxed to yield better prediction.

###  3.3 Transverse instability and wave development

####  3.3.1 Rayleigh-Taylor instability

Though the interfacial waves near the inlet are approximately longitudinal,
transverse modulations are also observed, see figure 3. The transverse
modulations are induced by the Rayleigh-Taylor (RT) instability. As shown in
figure 4, the longitudinal instability introduces oscillatory motion of the
interface in the vertical direction. When the interface accelerates toward the
liquid or decelerates toward the gas, the interface is unstable due to the
baroclinic effect and the RT instability is triggered, see figure 8. While the
longitudinal instability develops right at the end of the separator plate, the
transverse interfacial modulation may not grow immediately, see figure 3(a).
The reason is that the interface is stable in half time of one oscillation
cycle (when the interface decelerates toward the liquid or accelerates toward
the gas), see figure 8, so the transverse instability may first grow and then
decay.

Figure 8: Schematic for the transverse Rayleigh-Taylor (RT) instability due to
the vertical motion of interface induced by the longitudinal instability.
Figure 9: Transverse development of the Rayleigh-Taylor instability for the
case W0. (a) Top view of the interface at $t^{*}=323.75$, and the interface is
colored with $u$-velocity $(m/s)$. The vertical dashed line in (a) indicates
the wave location, and the interface profile along the dashed line is shown in
(b) on the $y$-$z$ plane to demonstrate the baroclinic effect.

The present results indicate that two conditions need to be satisfied for a
transverse modulation to grow and to transform the 2D longitudinal wave to
fully 3D. First, the longitudinal wave amplitude must be sufficiently large.
Second, the interface is decelerating toward the gas, so that the growing
interfacial wave will interact with the gas stream. A representative example
of the development of the transverse instability is shown in figure 9. The
later time evolution of this specific wave can be found in figure 3(c). At
this time ($t^{*}=323.75$) the interfacial wave is well aligned with
transverse direction if viewed from the top, with some small transverse
variations in the height of the wave crest. The interface is rising upward and
decelerating toward the gas, the misalignment between the pressure and density
gradients generates a baroclinic torque. The baroclinic torque in $x$
direction, $B_{x}=(\nabla\rho\times\nabla p)_{x}$ will cause the transverse
interfacial perturbation to grow, as indicated in figure 9(b), which shows the
interface profile on the $y$-$z$ plane. At this streamwise location
$x^{*}\approx 1.6$, the amplitude of longitudinal wave is large enough, as a
result, the wave crest will experience a significant interaction with fast gas
stream above the wave. The local acceleration of the gas flow above the wave
crest creates a Bernoulli depression, which will further pull the interface
upward (Hoepffner et al., 2011), amplifying the transverse modulation of the
interface, see $t^{*}=322.5$ to 327.5 in figure 3(c). More discussions about
the interaction between the interfacial wave and the gas stream will be
presented in section 3.4.1.

####  3.3.2 Dominant mode for transverse instability

Though the transverse modulation growth is assisted by its interaction with
the gas stream, the selection of the dominant transverse wavenumber is mainly
controlled by the RT instability. The inviscid linear temporal RT instability
with surface tension for a planar interface yields the following dispersion
relation (Chandrasekhar, 1961)

$\displaystyle\omega_{I}=\left(ak\frac{\rho_{l}-\rho_{g}}{\rho_{l}+\rho_{g}}-\frac{\sigma
k^{3}}{\rho_{l}+\rho_{g}}\right)^{1/2}\,,$ (10)

where $\omega_{I}$ is temporal growth rate, $k$ is wavenumber, and $a$ is
magnitude of the interfacial acceleration. To account for the viscous effect,
the dispersion relation becomes (Chandrasekhar, 1961; Joseph et al., 1999)

$\displaystyle\omega_{I}=-k^{2}\frac{\mu_{l}+\mu_{g}}{\rho_{l}+\rho_{g}}\pm\left(ak\frac{\rho_{l}-\rho_{g}}{\rho_{l}+\rho_{g}}-\frac{\sigma
k^{3}}{\rho_{l}+\rho_{g}}+k^{4}\left(\frac{\mu_{l}+\mu_{g}}{\rho_{l}+\rho_{g}}\right)^{2}\right)^{1/2}\,.$
(11)

Due to the low liquid viscosity in the present problem, the viscous effect on
the RT instability is very small. The difference between Eqs. (10) and (11) is
negligible, so the inviscid relation Eq. (10) will be used in the following
analysis.

Based on Eq. (10), the transverse wavenumber for the most unstable mode is

$\displaystyle k_{T}=\sqrt{\frac{a(\rho_{l}-\rho_{g})}{3\sigma}}\,,$ (12)

and the corresponding growth rate is

$\displaystyle\omega_{T}=\frac{2}{3^{3/2}}\left(\frac{a^{3}(\rho_{l}-\rho_{g})}{\sigma}\right)^{1/4}\,.$
(13)

The interfacial acceleration $a$ is induced by the longitudinal instability,
for which the characteristic length and time scales are $\delta$ and
$1/\omega_{L}$, respectively. Therefore, it can be approximated that
$a\sim\delta\omega_{L}^{2}$ and the dominant transverse wavenumber

$\displaystyle
k_{T}\approx\omega_{L}\sqrt{\frac{\delta(\rho_{l}-\rho_{g})}{3\sigma}}\,.$
(14)

As shown in section 3.2, $\omega_{L}$ increases with $I$. According to Eq.
(14), $k_{T}$ will also grow with $I$.

Using similar scaling relations, the temporal growth rate for the transverse
instability can be estimated as

$\displaystyle\omega_{T}\approx\frac{2}{3^{3/2}}\left(\frac{\delta^{3}\omega_{L}^{6}(\rho_{l}-\rho_{g})}{\sigma}\right)^{1/4}\,.$
(15)

Since $\omega_{L}\sim U_{D}/\delta$, it can be shown that

$\displaystyle{\omega_{T}}\approx{\omega_{L}}\text{We}_{L}^{1/4}\,,$ (16)

where $\text{We}_{L}$ is the Weber number for the longitudinal wave, expressed
as

$\displaystyle\text{We}_{L}=\frac{(\rho_{l}-\rho_{g})U_{D}^{2}\delta}{\sigma}\,.$
(17)

If ${\omega_{T}}$ is significantly lower than ${\omega_{L}}$, then the RT
instability will not grow fast enough to be relevant. For the present problem,
$\text{We}_{L}=9.45$, and it is estimated that ${\omega_{T}}/{\omega_{L}}\sim$
O(1). Therefore, the time scales for the longitudinal KH and the transverse RT
instabilities are comparable. In other words, the transverse RT instability
can grow in a time scale that is comparable to oscillatory interface motion
induced by the longitudinal instability. This affirms that RT instability is
responsible to the transverse development of the interfacial waves for the
present problem.

####  3.3.3 Transverse wavenumber spectra

To investigate the effect of inlet gas turbulence intensity $I$ on the
dominant transverse wavenumber, the wavenumber spectra of the transverse
interfacial modulation is examined. For a given $t^{*}$ and $x^{*}$, Fourier
transform is performed for the interfacial height $h$ along the transverse
direction. Then the wavenumber for the maximum amplitude in the spectrum,
$k_{\max}$, is measured. The contours of $k_{\max}$ for different $I$ are
shown in figure 10.

It is observed that $k_{\max}$ rises up to much larger values when
longitudinal waves pass by. This suggests that the development of the
transverse modulations are closely associated with the longitudinal waves, as
predicted by the stability analysis (Eq. (14)). As a result, the $k_{\max}$
contours also reveal important features about the formation and propagation in
the longitudinal waves. First of all, the trajectories of the longitudinal
waves appear as inclined straight lines of higher $k_{\max}$ values in the
figures. The slopes of the lines represent the wave propagation speeds in the
longitudinal direction. It is observed that the wave speeds are very similar
for all the cases and agree with the Dimotakis speed (Dimotakis, 1986),
$U_{D}$. For the present problem, $U_{D}^{*}=U_{D}/U_{g}=0.223$. As the
longitudinal wave grows in amplitude and interact with the fast gas stream,
the aerodynamic drag will cause the wave to accelerate. That is why the
longitudinal wave trajectories bend downward slightly further downstream,
which is most profound for the case W2 due to the faster spatial growth for
the longitudinal instability, see figure 10(c). The time period $\tau_{L}$ and
the wavelength $\lambda_{L}$ for the longitudinal waves can also be
identified, and they are related by $U_{D}$ as $\lambda_{L}=U_{D}\tau_{L}$. It
is also observed that $\tau_{L}$ decreases with $I$. Since $U_{D}$ remains
unchanged, $\lambda_{L}$ decreases with $I$.

Figure 10: The dominant transverse wavenumber $k_{\max}$ on the
$x^{*}$-$t^{*}$ diagrams for the cases (a) W0, (b) W1 and (c) W2. The vertical
dashed lines indicate the streamwise locations to estimate time-averaged
transverse wavenumber $\overline{k}_{T}$.

It is noted that only $k_{\max}$ on the longitudinal wave location
characterizes the transverse modulations of the interfacial waves. Between two
longitudinal waves, $k_{\max}$ represents the ripples on the interfaces, which
are less important. Therefore, we denote the $k_{\max}$ at the passages of the
longitudinal waves as the transverse wavenumber of the interfacial wave,
$k_{T}$. It is observed that $k_{T}$ for different waves are different due to
the chaotic nature of the turbulent multiphase flows. The temporal evolutions
of $k_{max}$ for different $I$ is shown in figures 11(a)-(c). For each case,
the measurement is made at $x^{*}$ locations corresponding to
$\log(\xi_{L}^{*})=-1$, i.e., $x^{*}=2.375$, 1.625, 1.25, for the cases W0,
W1, and W2, respectively, which are indicated by dashed lines in figure 10.
The RT instability predictions (Eq. (14)) are also shown in figure 11 for
comparison. It can be observed that $k_{T}$ generally increases when $I$
increases, which is consistent with the RT instability theory. The increase of
$k_{T}$ with $I$ is also consistent with the snapshots shown in figure 3.

The time average transverse wavenumber $\overline{k}_{T}$ is plotted as a
function of $I$ in figure 11(d). The measurement times for $k_{T}$ are
indicated by the vertical lines in figures 11(a)-(c), which correspond to the
longitudinal wave passage times estimated based on the dominant longitudinal
frequency $\omega^{*}_{L}$ for the given $I$. It is shown that
$\overline{k}_{T}$ increases from about 3330 to 5100 when $I$ increases from 0
to 0.13. The values of $\overline{k}_{T}$ for $I=0$ and 0.13 correspond to
about 2.5 and 3.9 waves in the domain width. The RT instability model (Eq.
(14)) well captures the increasing of trend of $\overline{k}_{T}$ over $I$.
The discrepancy between the model predictions and the simulation results can
be up to 17%, so the model predictions can only be used as approximations of
$\overline{k}_{T}$.

Figure 11: Temporal evolution of transverse wavenumber $k_{\max}$ for the
cases (a) W0, (b) W1, and (c) W2, respectively. The vertical lines indicated
the passage times of the longitudinal waves, which are estimated based on the
$\omega_{L}$for the corresponding $I$. (d) Variation of $\overline{k}_{T}$
with turbulent intensity $I$, compared with the RT predictions (Eq. (14)).

####  3.3.4 Transverse modulation amplitude

The wave amplitude for a given $(x^{*},z^{*})$ location can be characterized
by the root mean square of the interfacial height fluctuations normalized by
$H$, i.e., $\xi^{*}=(\overline{h^{\prime}h^{\prime}})^{1/2}/H$. The variations
of $\xi^{*}$ with $x^{*}$ and $z^{*}$ for different cases are shown in figure
12. The wave amplitude generally increases in the longitudinal direction,
though small variations in the transverse direction are also observed. It
should be noted that $\xi^{*}$ can be used to characterize the wave amplitude
only when the wave amplitude is small. The noises observed in figure 12(c) are
due to the rolling up and breakup of the waves. The transverse variations are
the most profound for the case W0.

Figure 12: Spatial variation of the interfacial wave amplitude
$\xi^{*}=(\overline{h^{\prime}h^{\prime}})^{1/2}/H$ for the cases (a) W0, (b)
W1, and (c) W2.

If $\xi^{*}$ is averaged in the transverse direction, we can obtain the
longitudinal wave amplitude, i.e.,
$\check{\xi}_{L}=\langle(\overline{h^{\prime}h^{\prime}})^{1/2}/H\rangle$. The
spatial variation of $\check{\xi}_{L}$ in $x^{*}$ direction is shown in figure
13(a). The results for the longitudinal wave amplitude $\xi_{L}^{*}$ defined
based on the contours of $\langle\overline{c^{\prime}c^{\prime}}\rangle=0.02$
(see figure 5(a)) are also shown for comparison. It can be observed that
though the values for $\check{\xi}_{L}$ and ${\xi}_{L}^{*}$ are different, the
ratio $\check{\xi}_{L}/{\xi}_{L}^{*}$ is approximately a constant for all
cases. Therefore, the spatial growth rates for the longitudinal instability
$\alpha_{L}^{*}$ will remain the same, no matter $\check{\xi}_{L}$ or
${\xi}_{L}^{*}$ are used. The curves of $\check{\xi}_{L}$ for the case W0
exhibits some fluctuations, which are due to stronger transverse variations of
$\xi^{*}$. The results of $\xi^{*}$ at different transverse locations are
shown in figure 14. It can be observed that the $\xi^{*}$-$x^{*}$ curves are
very different for the case W0, while the difference among different profiles
is much smaller for the cases W1 and W2. In spite of the different
$\xi^{*}$-$x^{*}$ curves at different transverse locations for the case W0,
the growth rates $\alpha_{L,W0}^{*}$ are quite similar, see figure 14(a).
Therefore, the measurement of $\alpha_{L}^{*}$ for the case W0 is not as
uncertain as it may appear in figures 6(a) and 13(a).

For a given $x^{*}$, the deviation of wave amplitude $\xi^{*}$ from the
transverse mean $\check{\xi}_{L}$ can be calculated as
$\xi^{*}-\check{\xi}_{L}=(\overline{h^{\prime}h^{\prime}})^{1/2}/H-\langle(\overline{h^{\prime}h^{\prime}})^{1/2}/H\rangle$.
Then the root mean square of the deviation,
$\check{\xi}_{T}=\langle(\xi^{*}-\check{\xi}_{L})^{2}\rangle^{1/2}$, can be
used to characterize the average amplitude of the transverse modulations,
which is plotted as a function of $x^{*}$ in figure 13(b). Near the inlet,
$\check{\xi}_{T}$ is at least an order of magnitude smaller than
$\check{\xi}_{L}$ ($\check{\xi}_{L}/\check{\xi}_{T}\gtrsim\exp(2.5)$). As the
longitudinal wave amplitude grows in $x^{*}$, the transverse modulation
amplitude also grows correspondingly. The spatial growth rates of
$\check{\xi}_{T}$ near the inlet are similar for different $I$. Yet since the
spatial growth of $\check{\xi}_{L}$ is higher for larger $I$, the wave
amplitude exceeds the linear regime at a smaller $x^{*}$. Beyond that point,
the interfacial waves start to deform and break, then $\check{\xi}_{T}$ is not
valid to characterize the transverse modulation amplitude. For the case W0,
since the longitudinal wave grows more slowly, $\check{\xi}_{T}$ has a wider
$x^{*}$ range to grow, reaching a higher value than for the cases W1 and W2.

Figure 13: Spatial variation of (a) longitudinal wave amplitude
$\check{\xi}_{L}$ and (b) transverse modulation amplitude $\check{\xi}_{T}$
along the $x^{*}$ direction for the cases W0, W1 and W2. The results for the
longitudinal wave amplitude $\xi_{L}^{*}$, which are measured based on the
liquid volume fraction fluctuations
$\langle\overline{c^{\prime}c^{\prime}}\rangle$, are also shown in (a) for
comparison. Figure 14: Spatial variation of wave amplitude $\xi^{*}$ at
different $z$ locations for the cases (a) W0, (b) W1 and (c) W2.

###  3.4 Development and disintegration of 3D interfacial waves

####  3.4.1 Wave-gas-stream interaction

When the wave amplitude grows beyond the linear regime, such as $\xi^{*}>0.1$,
the non-linear interaction between the interfacial wave and the fast gas
stream will influence the subsequent deformation of the wave and the liquid
sheet extended from the wave crest. The temporal evolutions of the interfaces,
the pressure, and the velocity magnitude for the case W0 are shown in figure
15. The first snapshots at $t^{*}=327.5$ correspond to approximately the end
of the linear regime, and the interfacial wave appears like a lobe of small
amplitude. The continuous growth of the wave introduces perpendicular exposure
of the wave crest to the fast gas stream. As a consequence of that, a high
pressure is built up on the upstream side of the wave, see figure 15(b). The
gas flow accelerates over the wave crest, introducing a low pressure above due
to the Bernoulli principle. Furthermore, the gas flow separates on the
downstream side (see figure 15(c)), which also introduces a low pressure wake
region. The pressure difference leads to an aerodynamic drag that pushes the
wave crest upward and also along the streamwise direction, further
destabilizing the interfacial wave. As a consequence of that, a liquid sheet
is extended from the wave crest, see figure 15(a) at $t^{*}=332.5$. The
Taylor-Culick rim is formed on the edge of the liquid sheet. It can be
observed the rim is not as smooth as that in the earlier time and exhibits
transverse variation with a higher wavenumber. Multiple mechanisms contribute
to the secondary transverse variations of the wave shape. At first, another
Rayleigh-Taylor instability is triggered due to the acceleration of the rim in
the longitudinal direction. Furthermore, the Rayleigh-Plateau instability of
the rim can also lead to variation in the rim diameter. Finally, the turbulent
wake of the wave and the interaction between the turbulent vortices and the
liquid sheet can also cause irregular deformation of the liquid sheet.

Figure 15: Temporal evolutions of the (a) interfaces, (b) pressure, and (c)
velocity magnitudes for the case W0, to demonstrate the interaction between
the interfacial wave and the gas stream. The interfaces in (a) are colored by
the streamwise velocity. The pressure and velocity shown in (b) and (c) are on
the $x-y$ plane through the wave. The frame is moving with the wave to keep
the wave located at the center of the window.

As the interfacial wave continues to grow, the interaction with the gas stream
becomes more intense and complex. The liquid sheet starts to bend and fold.
The thickness of the liquid sheet at the folding region becomes very small and
as a result, holes are formed in the liquid sheet. In figure 15(a), it can be
observed that at $t^{*}=337.5$ the folded segment of the liquid sheet has
completely broken, forming droplets and filaments of different sizes. A new
rim will form at the edge of the unbroken sheet that remains attached to the
wave, and another cycle of wave-gas-stream interaction will occur. Yet as the
interfacial wave is now further away from the inlet and the velocity of the
gas stream decreases along the streamwise direction due to turbulent
dissipation, the subsequent interaction between the interfacial wave and the
gas stream will be less intense.

####  3.4.2 Holes formation and sheet disintegration

The formation and development of the holes in the liquid sheet for the cases
W0 and W2 are shown in figure 16. The $w$-velocity is shown on the interfaces
and also on the cross-section $y$-$z$ planes at the hole centers. As can be
seen in figure 16(a) at $t^{*}=246.25$, the thickness of the liquid sheet is
highly uneven due to the sheet deformation. Furthermore, it can be observed
that the liquid is moving away from location of minimum thickness (see figure
16(a) at $t^{*}=246.25$). Eventually the two surfaces pinch and a hole is
formed. The pinching process observed here is reminiscent of pinching of
liquid neck in drop formation Castrejón-Pita et al. (2015); Zhang et al.
(2019). After the holes are formed, the Taylor-Culick (TC) rim develops on the
edge of the hole. Due to the capillary effect, the rim retracts radially,
making the holes expand. The speed of hole expansion is dictated by the rim
retraction velocity, namely the TC velocity,

$U_{\text{TC}}=\sqrt{\frac{2\sigma}{\rho e}}$ (18)

where $e$ is the liquid sheet thickness. The curvature of the rim has little
effect on the retraction velocity (Agbaglah, 2021). Since the thickness of the
liquid sheet is non-uniform, $e$ varies from $14.2$ to $6.5$ µm near the edge
of the hole, see figure 16(a) at $t^{*}=247.5$, for which the hole expansion
velocity predicted by theory (Eq. (18)) is $0.84-1.24$ m/s. The hole expansion
velocity measured from the simulation results is about 1.35m/s, which is in
reasonable agreement with the TC predictions. When the inlet gas is turbulent,
the gas flow around the sheet exhibits more intense fluctuations, which can be
recognized from the footprints on the interfacial velocity (see figure 16(b)).
The surrounding turbulence seems to have little effect on the hole expansion
speed, which is found to be similar to that for the case W0. This indicates
that the capillary effect still dominates the hole development. Nevertheless,
it is observed that the hole is formed closer to the rim. As a result, the
expansion of the hole will interact with the rim and cause the rim to break
(Agbaglah, 2021).

Figure 16: Formation and development of holes in liquid sheets for the cases
(a) W0 and (b) W2. The $w$-velocity is plotted on the interfaces and also on
the cross-section $y$-$z$ planes at the hole centers (indicated by the dashed
lines). The frame is moving with the hole to keep the hole located at the
center of the window.

####  3.4.3 Impact of inlet gas turbulence on sheet breakup

Though the inlet gas turbulence does not change the hole dynamics, it does
have a strong impact on the number of holes formed, and the locations where
the holes are formed. The breakups of the liquid sheet for the W0 and W2 cases
are shown in figure 17. For the case W0, the holes are formed at around
$x^{*}=4.275$, which is much further than the hole formation location
$x^{*}=3.025$ for the turbulent gas inlet. This is due to smaller growth rate
for the case W0, and as a result the interaction between the wave and the gas
stream occur further downstream. Nevertheless, the interfacial wave for the
case W0 reaches a larger amplitude and experiences a stronger interaction with
the gas stream. Furthermore, the interfacial wave is accelerated by the gas
stream to a higher velocity and is stretched to form a longer liquid sheet.
The folding of the liquid sheet enhances the reduction of sheet thickness (see
also figure 15). Multiple holes arise and expand simultaneously and the
merging of the holes leads a violent breakup of the folding segment of the
liquid sheet (Agbaglah, 2021). The rim of the liquid sheet is detached from
the wave, forming a filament aligned with the transverse direction, see
$t^{*}=238.75$ in figure 17(a). The detached filament interacts with the gas
stream and continues to break into droplets. Due to the irregular shape of the
filament, the breakup dynamics is different from the regular Rayleigh breakup
of a perturbed liquid cylinder (Villermaux et al., 2004).

For the case W2, the transverse wavenumber is higher than that for the case
W0. Therefore, the interfacial wave exhibits narrower lobes and more irregular
rims. In the snapshots shown in figure 17(b), only one hole is seen. By
examining all the time snapshots, it is affirmed that the number of holes
formed for the case W2 is significantly lower than that for the case W0. Since
only one hole is formed, instead of merging with other holes as for the case
W0, the hole expands and merges with the rim. The rim eventually breaks and
retracts to form two fingers aligned with the longitudinal direction. The
fingers are stretched by the gas stream and will break and form droplets.
Since the size and orientation of the fingers/filaments formed for the case W2
are distinct from those for the case W0, the resulting droplet statistics will
also be different.

Figure 17: Breakups of the liquid sheet due to hole formation and development
for the cases W0 and W2. The interfaces are colored with $u$-velocity.

###  3.5 Droplet statistics

To obtain statistics of droplets formed in the interfacial wave breakup, the
cells with liquid ($f>0$) that are connected together are tagged with the same
ID. Then the volume and centroid coordinates of individual liquid structures
can be computed. The details of the tagging approach can be found in previous
study (Ling et al., 2015).

####  3.5.1 Size distribution

The droplet size distribution for different cases are shown in figure 18. Only
the droplets in the region $x^{*}<14$ are considered, to exclude the effect of
the outflow boundary. Furthermore, sampling is only performed after the
droplet statistics has reached the statistically stationary state (at about
$t^{*}=200$). For all cases, droplets are collected every 50 time steps, so in
total about $N_{s}=$3000 samples have been collected. The sampling time has
been verified to be sufficiently long to yield statistically converged
results. Since the droplet number generally decreases with the droplet
diameter $d$, the size interval $\Delta d$ is increased over $d$ with a
scaling ratio 1.15, namely $\Delta d_{i+1}=1.15\Delta d_{i}$, where $i$ is the
interval number. The number of droplets falling into the size interval
centered with the diameter $d$ for all samples is denoted as $n$. The average
droplet number is then defined as $\overline{n}=n/N_{s}$.

The results for $\overline{n}$ as a function of $d$ for the cases W0, W1, and
W2 are shown in figure 18(a). It is observed that $\overline{n}$ for all $d$
generally decreases when $I$ increases. This is consistent with the
observation in figure 17, namely the interfacial wave breakup is less violent
when $I$ increases, and as a results, fewer droplets are formed. Nevertheless,
the shapes of the distribution profiles for different $I$ are quite similar,
which motivates the examination of the probability density function (PDF). The
droplet size PDF based on the droplet number is defined as

$P_{n}(d)=\frac{\overline{n}(d)}{\Delta d\sum\overline{n}(d)}\,.$ (19)

Figure 18: (a) Droplet size distribution and (b) PDF based on the time-average
droplet number for the cases W0, W1, and W2.

The droplet size PDF for different $I$ are shown in figure 18(b). The profile
of $P_{n}$ is quite complex and cannot be represented by any typical
distribution function. The recent study by Balachandar et al. (2020) indicates
that a combined log-normal and Pareto distribution is required to represent
the complex PDF for the respiratory droplets formed in a coughing event
(Duguid, 1946). While the small droplets follow the log-normal distribution,
the larger droplets follow the Pareto distribution. The drop size PDF for the
present problem for different ranges of droplet size also seem to follow
different distribution functions, i.e., the log-normal, Pareto, and
exponential distribution functions. The expressions for these three functions
are given as

$\displaystyle P_{L}(d)$
$\displaystyle=\frac{B_{L}}{{d\hat{\sigma}\sqrt{2\pi}}}\exp\left[-\frac{(\ln
d-\hat{\mu})^{2}}{2\hat{\sigma}^{2}}\right]\,,$ (20) $\displaystyle P_{P}(d)$
$\displaystyle=\frac{B_{P}}{d^{\alpha}}\,,$ (21) $\displaystyle P_{E}(d)$
$\displaystyle=B_{E}\exp(d/\lambda)\,,$ (22)

where $B_{L}$, $B_{P}$, and $B_{E}$ are the corresponding normalization
constants for the log-normal, Pareto, and exponential functions, respectively.
For the log-normal function, $\hat{\mu}$ and $\hat{\sigma}$ are the expected
value and the standard deviation of $\ln d$. While $\alpha$ is the power index
for the Pareto distribution, $\lambda$ is the characteristic length for the
exponential distribution. To compare the simulation results with the
distribution functions, we have used the data for the case W1 to fit these
distributions functions. For convenience of discussions, the droplets for
$d\lesssim 25$ µm, $25\ \lesssim d\lesssim 140$ µm, and $d\gtrsim 140$ µm are
referred to as _small_ , _medium_ and _large_ droplets, the data of which are
used to fit Eqs. (20), (21), and (22), respectively. The boundaries between
different ranges vary slightly with $I$. The values for the fitting parameters
are $B_{L}=1.14,\ \hat{\sigma}=0.35,\ \hat{\mu}=2.24$ for the log-normal
function, $B_{P}=0.099,\ \alpha=1.04$ for the Pareto function, and
$B_{E}=0.0035,\ \lambda=95.0$ for the exponential function. The units for
these parameters can be determined based on the units for $P_{n}$ and $d$,
i.e., (µm)-1 and µm, respectively. As shown in figure 18(b), the fitted
distribution functions generally agree well with the simulation results.

For the large droplets ($d\gtrsim 140$ µm), the curves of $P_{n}$ for
different $I$ approximately collapse to the same exponential function. The
exponential tail in the drop size distribution has been observed in numerous
experiments and simulations (Simmons, 1977a, b; Marmottant & Villermaux, 2004;
Ling et al., 2015). For droplets of medium sizes, $25\ \lesssim d\lesssim 140$
µm, the PDF exhibits a power-law decay, namely the Pareto distribution. The
width of the Pareto distribution region increases with the the inlet gas
turbulence intensity $I$. As a result, the power-law decay is the most
profound for the case W2. For the small droplets $d\lesssim 25$ µm, the shape
of $P_{n}$ is significantly different from that for the medium droplets. The
different PDF profiles may be due to the different droplet formation
mechanisms. While the large number of small droplets are probably generated
directly in the breakup of liquid sheets, see figure 17(a), the medium
droplets are produced in the Rayleigh-Plateau breakup of thicker ligaments and
fingers, see figure 17(b). It will be interesting for future study to identify
$P_{n}$ for different breakup events and to further reveal the physical
reasons behind the different $P_{n}$ for the small and medium droplets.
Following former studies (Marty, 2015; Ling et al., 2017; Balachandar et al.,
2020), the log-normal function has been used to fit $P_{n}$ for the small
droplets, though the narrow range of $d$ may have placed some uncertainty in
the fitting process and whether the log-normal function is indeed the best fit
for the small droplets. It should be noted that the droplets smaller than 8 µm
are not accounted here, since those droplets are smaller than about one cell
and the formation of these droplets are likely under resolved. To better
examine the distribution function for the small droplets, more data points are
required for $d<8$ µm, which would require further refined simulations. Such
simulations require computer time at least an order of magnitude higher than
that for the present simulations and will be relegated to future works.

####  3.5.2 Spatial distribution

The spatial distributions of droplet number $\overline{n}$ in the longitudinal
and vertical directions for different droplet size ranges and the cases W0,
W1, and W2 are shown in figure 19. In the longitudinal direction,
$\overline{n}$ is zero until the interfacial waves break. After that,
$\overline{n}$ increases with $x^{*}$, approximately following the hyperbolic
tangent function, see figure 19(a). For the small and medium droplets,
$\overline{n}$ reaches the maximum around $10\lesssim x^{*}\lesssim 13$. The
decrease in $\overline{n}$ for $x^{*}\gtrsim 13$ is due to the coalescence of
droplets and the merge of droplets back on the unbroken liquid layer at the
bottom (see figure 3(c) and (d)). For large droplets, the decrease of
$\overline{n}$ is less profound, and $\overline{n}$ reaches a plateau for
large $x^{*}$. When the inlet gas turbulence intensity $I$ increases, several
changes are induced. First, the onset of the increase of $\overline{n}$ occurs
at smaller $x^{*}$ due to the faster growth of the longitudinal wave and the
earlier breakup of the interfacial waves. Second, the increase of
$\overline{n}$ becomes more gradual and eventually settles on a lower plateau.
As a result, the curves for W0 and W1 cross at a critical longitudinal
location $x^{*}_{c}$. For $x^{*}<x^{*}_{c}$,
$\overline{n}_{\text{W1}}>\overline{n}_{\text{W0}}$, or vice versa. The values
of $x^{*}_{c}$ for the small and medium droplets are similar, at about
$x^{*}_{c}=7$, while that for the large droplets is larger, at around
$x^{*}_{c}=8.2$.

Figure 19: Distribution of $\overline{n}$ in (a) streamwise and (b) vertical
directions for the cases W0, W1, and W2. The droplets are collected every 50
time steps. The three columns are for the small (12.5–100 µm), medium(100– 200
µm), and large (200–300 µm) droplets, respectively.

The droplet number distributions in the vertical direction for different
droplet size ranges are shown in figure 19(d)-(f). For the small droplets, the
maxima of $\overline{n}$ are at the bottom plane of domain for all cases. The
droplet number then decreases rapidly with $y^{*}$ until $y^{*}\approx 1$,
where the interface is located at the inlet and the decreases of
$\overline{n}$ becomes more gradual. After that, $\overline{n}$ decreases
approximately exponentially to zero, see figure 19(d). The accumulation of
small droplets near the bottom domain is probably due to the fact that the
small droplets have smaller inertia and thus are easier to be carried by the
gas stream toward the bottom downstream. The medium and larger droplets have
larger inertia and tend to maintain the velocity direction when they are
formed. As a result, we see the maxima of $\overline{n}$ for the medium and
the large droplets are near $y^{*}=1$, instead. The vertical location
corresponding to the maximized $\overline{n}$ is denoted by $y^{*}_{\max}$
(see figure 19(e)). It is observed that for the medium and large droplets,
$y^{*}_{\max}$ decreases with the inlet gas turbulence intensity $I$. This is
again due to the faster growth of the longitudinal wave for the cases with
larger $I$, and as a result, the breakups of the interfacial waves occur in a
lower vertical location (see also figure 3).

##  4 Conclusions

The effect of inlet gas turbulence on the formation, development, and breakup
of the interfacial waves in a two-phase mixing layer is investigated through
direct numerical simulation. The gas and liquid properties and the mean
injection velocities are chosen to place the present case in the absolute
instability regime. Turbulent velocity fluctuations are induced at the gas
inlet and a parametric study for the inlet gas turbulence intensity $I$ is
performed. The governing equations are solved using the open-source multiphase
flow solver, _PARIS_ , in which the mass-momentum consistent volume-of-fluid
method is used to capture the sharp gas-liquid interfaces. To allow a detailed
investigation of the transverse development of the interfacial waves, two
different domain widths have been used. The wide domain consists of about 1.6
billions cells.

The temporal evolutions and the frequency spectra for the interfacial height
in the near field shows a dominant frequency, confirming that the case is in
the absolute instability regime. The dominant frequency varies little when
$I\lesssim 0.02$ but then increases over $I$ after the threshold. Similar
increasing trend over $I$ has been observed for the spatial growth rate of the
longitudinal wave amplitude. When the inlet gas is laminar, the dominant
frequencies and spatial growth rates for the narrow and wide domains are
similar. When the turbulence is present at the gas inlet, deviations between
the results for the narrow and wide domains are observed.

The longitudinal instability introduces vertical motion of the interface. When
the interface moves upward and decelerates toward the gas, the Rayleigh-Taylor
(RT) instability is triggered, introducing transverse modulations on the
interfacial waves. The transverse variation is amplified as the wave interacts
with the fast gas stream. The selection of dominant transverse wavenumber is
dictated by the RT instability and can be determined by inviscid temporal
stability theory with surface tension. The transverse wavenumber scales with
the longitudinal frequency and thus will also increase with $I$. The
simulation results for the transverse wavenumber agree well with theoretical
predictions. The ratio between the growth rate of the transverse instability
and the dominant frequency of the longitudinal instability is estimated to be
comparable, affirming the RT instability is responsible for the initiation of
the transverse modulations for the present problem.

The subsequent development of the interfacial wave is controlled by the
nonlinear interaction with the gas stream. When the interfacial wave invades
in the gas stream, an aerodynamic drag is developed, which further
destabilizes the interfacial waves and stretches the wave crest to form liquid
lobes/sheets. For the laminar gas inlet, the liquid sheet deforms and folds
when it interacts with the gas stream, leading to the formation of multiple
holes in the folding area. The expansion and merging of the holes eventually
lead to a violent disintegration of the liquid sheet, detaching filaments that
are aligned with the transverse direction. In contrast, for the turbulent gas
inlet, the lobes are narrower and the rim is less regular. Fewer holes are
formed and the expansion of the holes interacts with the rim, forming fingers
aligned with the longitudinal direction. The different breakup dynamics of the
liquid sheet results in different sizes and orientations of the filaments. The
droplets formed due to the subsequent breakup of the filaments thus exhibit
different statistics.

The numbers of droplets for all sizes reduce when $I$ increases, though the
probability density functions (PDF) appear to be similar. The droplets of
different size ranges are found to follow different distribution functions,
i.e., the log-normal, Pareto, and exponential functions for the small ($8\
\lesssim d\lesssim 25$ µm), medium ($25\ \lesssim d\lesssim 140$ µm), and
large ($d\gtrsim 140$ µm) droplets, respectively. The spatial distributions of
the droplet number for different size ranges have also been presented. The
increase of the droplet number along the longitudinal direction is similar to
the hyperbolic tangent function and the decrease of droplet number along the
vertical distance away from the mixing layer follows the exponential function.
When $I$ increases, more droplets are formed in the near field and near the
bottom of the domain.

## Acknowledgement

This research is supported by the National Science Foundation (#1942324). The
authors also acknowledge the Extreme Science and Engineering Discovery
Environment (XSEDE) for providing the computational resources that have
contributed to the research results reported in this paper. The Baylor High
Performance and Research Computing Services (HPRCS) have been used to process
the simulation data. We also thank Dr. Gretar Tryggvason, Dr. Stephane
Zaleski, and other developers of for their contributions to the _PARIS_
solver.

## Declaration of Interests

The authors report no conflict of interest.

## References

  * Agbaglah et al. (2017) Agbaglah, G., Chiodi, R. & Desjardins, O. 2017 Numerical simulation of the initial destabilization of an air-blasted liquid layer. J. Fluid Mech. 812, 1024–1038.
  * Agbaglah et al. (2013) Agbaglah, G., Josserand, C. & Zaleski, S. 2013 Longitudinal instability of a liquid rim. Phys. Fluids 25, 022103\.
  * Agbaglah (2021) Agbaglah, G. G. 2021 Breakup of thin liquid sheets through hole–hole and hole–rim merging. J. Fluid Mech. 911, A23.
  * Ambravaneswaran et al. (2002) Ambravaneswaran, B., Wilkes, E. D. & Basaran, O. A. 2002 Drop formation from a capillary tube: Comparison of one-dimensional and two-dimensional analyses and occurrence of satellite drops. Phys. Fluids 14, 2606–2621.
  * Aniszewski et al. (2021) Aniszewski, W., Arrufat, T., Crialesi-Esposito, M., Dabiri, S., Fuster, D., Ling, Y., Lu, J., Malan, L., Pal, S., Scardovelli, R., Tryggvason, G., Yecko, P. & Zaleski, S. 2021 PArallel, Robust, Interface Simulator (PARIS). Comput. Phys. Comm. 263, 107849.
  * Arrufat et al. (2020) Arrufat, T., Crialesi-Esposito, M., Fuster, D., Ling, Y., Malan, L., Pal, S., Scardovelli, R., Tryggvason, G. & Zaleski, S. 2020 A momentum-conserving, consistent, volume-of-fluid method for incompressible flow on staggered grids. Comput. Fluids 215, 104785.
  * Aulisa et al. (2007) Aulisa, E., Manservisi, S., Scardovelli, R. & Zaleski, S. 2007 Interface reconstruction with least-squares fit and split advection in three-dimensional cartesian geometry. J. Comput. Phys. 225, 2301–2319.
  * Balachandar et al. (2020) Balachandar, S., Zaleski, S., Soldati, A., Ahmadi, G. & Bourouiba, L. 2020 Host-to-host airborne transmission as a multiphase flow problem for science-based social distance guidelines. Int. J. Multiphase Flow 132, 103439.
  * Boeck & Zaleski (2005) Boeck, T. & Zaleski, S. 2005 Viscous versus inviscid instability of two-phase mixing layers with continuous velocity profile. Phys. Fluids 17, 032106.
  * Castrejón-Pita et al. (2015) Castrejón-Pita, J. R., Castrejón-Pita, A. A., Thete, S. S., Sambath, K., Hutchings, I. M., Hinch, J., Lister, J. R. & Basaran, O. A. 2015 Plethora of transitions during breakup of liquid filaments. Proc. Natl. Acad. Sci. U.S.A. 112, 4582–4587.
  * Chandrasekhar (1961) Chandrasekhar, S. 1961 Hydrodynamic and hydromagnetic stability. Oxford University Press.
  * Chaussonnet et al. (2020) Chaussonnet, G., Gepperth, S., Holz, S., Koch, R. & Bauer, H.-J. 2020 Influence of the ambient pressure on the liquid accumulation and on the primary spray in prefilming airblast atomization. Int. J. Multiphase Flow 125, 103229.
  * Chorin (1968) Chorin, A. J. 1968 Numerical solution of the Navier-Stokes equations. Math. Comput. 22, 745–762.
  * Delon et al. (2018) Delon, A., Cartellier, A. & Matas, J.-P. 2018 Flapping instability of a liquid jet. Phys. Rev. Fluids 3, 043901\.
  * Dimotakis (1986) Dimotakis, P. E. 1986 Two-dimensional shear-layer entrainment. AIAA J. 24, 1791–1796.
  * Duguid (1946) Duguid, J. P. 1946 The size and the duration of air-carriage of respiratory droplets and droplet-nuclei. Epidemiol. Infect. 44, 471–479.
  * Eggers (1993) Eggers, J. 1993 Universal pinching of 3D axisymmetric free-surface flow. Phys. Rev. Lett. 71, 3458.
  * Francois et al. (2006) Francois, M. M., Cummins, S. J., Dendy, E. D., Kothe, D. B., Sicilian, J. M. & Williams, M. W. 2006 A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. J. Comput. Phys. 213, 141–173.
  * Fuster et al. (2013) Fuster, D., Matas, J. P., Marty, S., Popinet, S., J., Hoepffner, Cartellier, A. & Zaleski, S. 2013 Instability regimes in the primary breakup region of planar coflowing sheets. J. Fluid Mech 736, 150–176.
  * Healey (2007) Healey, J. J. 2007 Enhancing the absolute instability of a boundary layer by adding a far-away plate. J. Fluid Mech. 579, 29\.
  * Herrmann (2011) Herrmann, M. 2011 The influence of density ratio on the primary atomization of a turbulent liquid jet in crossflow. Proc. Combust. Inst. 33, 2079–2088.
  * Hoepffner et al. (2011) Hoepffner, J., Blumenthal, R. & Zaleski, S. 2011 Self-similar wave produced by local perturbation of the Kelvin-Helmholtz shear-layer instability. Phys. Rev. Lett. 106, 104502.
  * Jarrahbashi & Sirignano (2014) Jarrahbashi, D. & Sirignano, W. A. 2014 Vorticity dynamics for transient high-pressure liquid injection a. Phys. Fluids 26, 73.
  * Jarrahbashi et al. (2016) Jarrahbashi, D., Sirignano, W. A., Popov, P. P. & Hussain, F. 2016 Early spray development at high gas density: hole, ligament and bridge formations. J. Fluid Mech. 792, 186–231.
  * Jiang & Ling (2020) Jiang, D. & Ling, Y. 2020 Destabilization of a planar liquid stream by a co-flowing turbulent gas stream. Int. J. Multiphase Flow 122, 103121.
  * Joseph et al. (1999) Joseph, D. D., Belanger, J. & Beavers, G. S. 1999 Breakup of a liquid drop suddenly exposed to a high-speed airstream. Int. J. Multiphase Flow 25, 1263–1303.
  * Juniper (2006) Juniper, M. P 2006 The effect of confinement on the stability of two-dimensional shear flows. J. Fluid Mech. 565, 171–195.
  * Klein et al. (2003) Klein, M., Sadiki, A. & Janicka, J. 2003 A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations. J. Comput. Phys. 186, 652–665.
  * Kooij et al. (2018) Kooij, S., Sijs, R., Denn, M. M., Villermaux, E. & Bonn, D. 2018 What determines the drop size in sprays? Phys. Rev. X 8, 031019.
  * Lasheras et al. (1998) Lasheras, J. C., Villermaux, E. & Hopfinger, E. J. 1998 Break-up and atomization of a round water jet by a high-speed annular air jet. J. Fluid Mech. 357, 351–379.
  * Lefebvre (1988) Lefebvre, A. 1988 Atomization and sprays. CRC press.
  * Lefebvre (1980) Lefebvre, A. H. 1980 Airblast atomization. Prog. Energ. Combust. Sci. 6, 233–261.
  * Li (1995) Li, J. 1995 Calcul d’interface affine par morceaux (piecewise linear interface calculation). C. R. Acad. Sci. Paris, série II b 320, 391–396.
  * Ling et al. (2019) Ling, Y., Fuster, D., Tryggvasson, G. & Zaleski, S. 2019 A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability. J. Fluid Mech. 859, 268–307.
  * Ling et al. (2017) Ling, Y., Fuster, D., Zaleski, S. & Tryggvason, G. 2017 Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup. Phys. Rev. Fluids 2, 014005.
  * Ling et al. (2015) Ling, Y., Zaleski, S. & Scardovelli, R. 2015 Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model. Int. J. Multiphase Flow 76, 122–143.
  * Marmottant & Villermaux (2004) Marmottant, P. & Villermaux, E. 2004 On spray formation. J. Fluid Mech. 498, 73–111.
  * Marston et al. (2016) Marston, J. O., Truscott, T. T., Speirs, N. B., Mansoor, M. M. & Thoroddsen, S. T. 2016 Crown sealing and buckling instability during water entry of spheres. J. Fluid Mech. 794, 506–529.
  * Marty (2015) Marty, S. 2015 Contribution a l’étude de l’atomisation assistée d’un liquide. PhD thesis, Université de Grenoble.
  * Matas (2015) Matas, J.-P. 2015 Inviscid versus viscous instability mechanism of an air–water mixing layer. J. Fluid Mech. 768, 375–387.
  * Matas et al. (2018) Matas, J.-P., Delon, A. & Cartellier, A. 2018 Shear instability of an axisymmetric air–water coaxial jet. J. Fluid Mech. 843, 575–600.
  * Matas et al. (2011) Matas, J.-P., Marty, S. & Cartellier, A. 2011 Experimental and analytical study of the shear instability of a gas-liquid mixing layer. Phys. Fluids 23, 094112.
  * Matas et al. (2015) Matas, J.-P., Marty, S., Dem, M. S. & Cartellier, A. 2015 Influence of gas turbulence on the instability of an air-water mixing layer. Phys. Rev. Lett. 115, 074501.
  * O’Naraigh et al. (2013) O’Naraigh, L., Spelt, P. D. M. & Shaw, S. J. 2013 Absolute linear instability in laminar and turbulent gas–liquid two-layer channel flow. J. Fluid Mech. 714, 58–94.
  * Opfer et al. (2014) Opfer, L., Roisman, I. V., Venzmer, J., Klostermann, M. & Tropea, C. 2014 Droplet-air collision dynamics: Evolution of the film thickness. Phys. Rev. E 89, 013023\.
  * Otto et al. (2013) Otto, T., Rossi, M. & Boeck, T. 2013 Viscous instability of a sheared liquid-gas interface: Dependence on fluid properties and basic velocity profile. Phys. Fluids 25, 032103.
  * Popinet (2009) Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228 (16), 5838–5866.
  * Raynal (1997) Raynal, L. 1997 Instabilite et entrainement a l’interface d’une couche de melange liquide-gaz. PhD thesis, Université Joseph Fourier - Grenoble I.
  * Renardy & Renardy (2002) Renardy, Y. & Renardy, M. 2002 PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method. J. Comput. Phys. 183, 400–421.
  * Roisman (2010) Roisman, I. V. 2010 On the instability of a free viscous rim. J. Fluid Mech. 661, 206–228.
  * Rudman (1998) Rudman, M. 1998 A volume-tracking method for incompressible multifluid flows with large density variations. Int. J. Numer. Meth. Fluids 28, 357–378.
  * Scardovelli & Zaleski (2003) Scardovelli, R. & Zaleski, S. 2003 Interface reconstruction with least-square fit and split eulerian–lagrangian advection. Int. J. Numer. Meth. Fluids 41 (3), 251–274.
  * Simmons (1977a) Simmons, H. C. 1977a The correlation of drop-size distributions in fuel nozzle sprays—part i: The drop-size/volume-fraction distribution. J. Engng. Power 7, 309–314.
  * Simmons (1977b) Simmons, H. C. 1977b The correlation of drop-size distributions in fuel nozzle sprays—part ii: The drop-size/number distribution. J. Engng. Power 7, 315–319.
  * Snoeijer & Andreotti (2013) Snoeijer, J. H. & Andreotti, B. 2013 Moving contact lines: scales, regimes, and dynamical transitions. Annu. Rev. Fluid Mech. 45, 269–292.
  * Tryggvason et al. (2011) Tryggvason, G., Scardovelli, R. & Zaleski, S. 2011 Direct numerical simulations of gas-liquid multiphase flows. Cambridge University Press.
  * Varga et al. (2003) Varga, C. M., Lasheras, J. C. & Hopfinger, E. J. 2003 Initial breakup of a small-diameter liquid jet by a high-speed gas stream. J. Fluid Mech. 497, 405–434.
  * Vaudor et al. (2017) Vaudor, G., Ménard, T., Aniszewski, W., Doring, M. & Berlemont, A. 2017 A consistent mass and momentum flux computation method for two phase flows. Application to atomization process. Comput. Fluids 152, 204–216.
  * Villermaux et al. (2004) Villermaux, E., Marmottant, Ph. & Duplat, J. 2004 Ligament-mediated spray formation. Phys. Rev. Lett. 92, 074501\.
  * Zandian et al. (2018) Zandian, A., Sirignano, W. A. & Hussain, F. 2018 Understanding liquid-jet atomization cascades via vortex dynamics. J. Fluid Mech. 843, 293–354.
  * Zhang et al. (2019) Zhang, B., Ling, Y., Tsai, P.-H., Wang, A.-B., Popinet, S. & Zaleski, S. 2019 Short-term oscillation and falling dynamics for a water drop dripping in quiescent air. Phys. Rev. Fluids 4, 123604.
  * Zhang et al. (2020) Zhang, B., Popinet, S. & Ling, Y. 2020 Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under non-evaporative operating conditions. Int. J. Multiphase Flow 130, 103362.