2.1 Governing equations

The governing equations are the continuity, momentum and LS/VoF equations. Since both the gas and the liquid are incompressible, the continuity equation is

(2.1) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the velocity vector. The momentum equation, including the viscous diffusion and surface tension forces and neglecting the gravitational force, is

(2.2) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}\boldsymbol{u})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u}\boldsymbol{u})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(2\unicode[STIX]{x1D707}\unicode[STIX]{x1D63F})-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FF}(d)\boldsymbol{n}, & & \displaystyle\end{eqnarray}$$

where $p$ , $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D707}$ are the pressure, density and dynamic viscosity of the fluid, respectively. $\unicode[STIX]{x1D63F}$ is the rate of deformation tensor,

(2.3) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D63F}={\textstyle \frac{1}{2}}\left[(\unicode[STIX]{x1D735}\boldsymbol{u})+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}\right]. & & \displaystyle\end{eqnarray}$$

The last term in (2.2) is the surface tension force per unit volume $\boldsymbol{F}=-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FF}(d)\boldsymbol{n}$ ; where $\unicode[STIX]{x1D70E}$ is the surface tension coefficient, $\unicode[STIX]{x1D705}$ is the surface curvature, $\unicode[STIX]{x1D6FF}(d)$ is the Dirac delta function, $d$ is the distance from the interface, and $\boldsymbol{n}$ is the unit vector normal to the liquid–gas interface.

The LS function $\unicode[STIX]{x1D719}$ is defined as a smooth distance function, which enables evaluation of the density, viscosity and surface tension at any distance from the interface in either gas or liquid zones. This method was developed by Osher and his coworkers (Zhao et al. Reference Zhao, Chan, Merriman and Osher1996; Sussman et al. Reference Sussman, Fatemi, Smereka and Osher1998; Osher & Fedkiw Reference Osher and Fedkiw2001). In this algorithm, the interface $\unicode[STIX]{x1D6E4}$ is the zero level set of $\unicode[STIX]{x1D719}$ ,

(2.4) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6E4}=\left\{\boldsymbol{x}|\unicode[STIX]{x1D719}(\boldsymbol{x},t)=0\right\}. & & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D719}<0$ in the liquid region and $\unicode[STIX]{x1D719}>0$ in the gas region are taken. $\boldsymbol{u}$ is continuous across the interface. Since the interface moves with the fluid particles, the evolution of $\unicode[STIX]{x1D719}$ is then given by

(2.5) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0, & & \displaystyle\end{eqnarray}$$

which is called the LS equation. If the initial distribution of the level set is a signed distance function, after a finite time of being convected by a non-uniform velocity field, it will not remain a distance function. Therefore, the LS function should be re-initialized in such a way that it will be a distance function, without changing the zero level set (position of the interface). This is achieved by solving the following differential equation (Sussman et al. Reference Sussman, Fatemi, Smereka and Osher1998);

(2.6) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x2202}d}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=\text{sign}(\unicode[STIX]{x1D719})(1-|\unicode[STIX]{x1D735}d|), & & \displaystyle\end{eqnarray}$$

with the initial condition

(2.7) $$\begin{eqnarray}\displaystyle \displaystyle d(\boldsymbol{x},0)=\unicode[STIX]{x1D719}_{0}(\boldsymbol{x}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}$ is a pseudotime. The steady solutions of (2.6) are distance functions. Furthermore, since $\text{sign}(0)=0$ , then $d(\boldsymbol{x},t)$ has the same zero level set as $\unicode[STIX]{x1D719}(\boldsymbol{x})$ .

Since the density and viscosity are constant in each fluid, they take on two different values depending on the sign of $\unicode[STIX]{x1D719}$ ; we can write

(2.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D70C}_{l}+(\unicode[STIX]{x1D70C}_{g}-\unicode[STIX]{x1D70C}_{l})\text{H}(\unicode[STIX]{x1D719}), & & \displaystyle\end{eqnarray}$$

and

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D707}_{l}+(\unicode[STIX]{x1D707}_{g}-\unicode[STIX]{x1D707}_{l})\text{H}(\unicode[STIX]{x1D719}), & & \displaystyle\end{eqnarray}$$

where $\text{H}(\unicode[STIX]{x1D719})$ is a smoothed Heaviside function, and the subscripts $g$ and $l$ refer to gas and liquid, respectively. Using these expressions, the governing equation for the fluid velocity $\boldsymbol{u}$ (2.2) can be written as a single equation containing both liquid and gas properties.

At low density ratios, we use a transport equation similar to the level set (2.5) for the volume fraction $f$ , also called the VoF variable, in order to describe the temporal and spatial evolution of the two-phase flow (Hirt & Nichols Reference Hirt and Nichols1981);

(2.10) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f=0, & & \displaystyle\end{eqnarray}$$

where the VoF variable $f$ represents the volume of (liquid phase) fluid fraction as follows:

(2.11) $$\begin{eqnarray}\displaystyle f(t)=\left\{\begin{array}{@{}cl@{}}0, & \text{outside of liquid phase}\\ 0<f<1, & \text{at the interface}\\ 1, & \text{inside the liquid phase}.\end{array}\right. & & \displaystyle\end{eqnarray}$$

The normal direction of the fluid interface is found where the value of $f$ changes most rapidly. With this method, the free surface is not defined sharply; instead, it is distributed over the height of the cell. Since fluid properties are required every time step in order to solve the Navier–Stokes equations, the density and viscosity change continuously based on (2.8) and (2.9). However, in these equations, the argument of the Heaviside function is replaced by $f$ , in the VoF method.

The formulation for the fully conservative momentum convection and volume fraction transport, the momentum diffusion, and the surface tension are treated explicitly. To ensure a sharp interface of all flow discontinuities and to suppress numerical dissipation of the liquid phase, the interface is reconstructed at each time step by the PLIC (piecewise linear interface calculation) method proposed by Rider & Kothe (Reference Rider and Kothe1998). The liquid phase is transported on the basis of its reconstructed distribution. The capillary effects in the momentum equations are represented by a capillary tensor as introduced by Scardovelli & Zaleski (Reference Scardovelli and Zaleski1999).

2.2 Flow configuration

The computational domain, shown in figure 3, consists of a cube, which is discretized into uniform-sized cells. The liquid segment, which is a sheet of thickness $h$ ( $h=50~\unicode[STIX]{x03BC}\text{m}$ for the thin sheet and $200~\unicode[STIX]{x03BC}\text{m}$ for the thick sheet in this study), is centred at the middle of the box. The domain size in terms of the sheet thickness $h$ is $16h\times 10h\times 10h$ , in the $x$ , $y$ and $z$ directions, respectively, for the thin sheet, and $4h\times 4h\times 8h$ for the thick sheet. The liquid segment is surrounded by gas zones on top and bottom. The gas moves in the positive $x$ - (streamwise) direction, with a constant velocity ( $U=100~\text{m}~\text{s}^{-1}$ ) at the top and bottom boundaries, and its velocity diminishes to the interface velocity with a boundary layer thickness obtained from 2D full-jet simulations. In the liquid, the velocity decays to zero at the centre of the sheet with a hyperbolic tangent profile. For more detailed description of the initial conditions, see figure 12 of Zandian, Sirignano & Hussain (Reference Zandian, Sirignano and Hussain2016).

Figure 3.

The computational domain with the initial liquid and gas zones.

Our study involves a temporal computational analysis with a relative velocity between the two phases. Due to friction, the relative velocity decreases with time. Furthermore, the domain is several wavelengths long in the streamwise direction so that some spatial development occurs. In order to reduce the dependence of the results on details of the boundary conditions, specific configurations (e.g. air-assist or air-blast atomization) are avoided. However, calculations are made with the critical non-dimensional parameters in the ranges of practical interest.

The liquid–gas interface is initially perturbed symmetrically on both sides with a sinusoidal profile of predefined wavelength and amplitude obtained from the 2D full-jet simulations (see Zandian et al. Reference Zandian, Sirignano and Hussain2016). Two analyses without forced or initial surface perturbations – the full-jet 2D simulations (figure 11 of Zandian et al. Reference Zandian, Sirignano and Hussain2016) and the initially non-perturbed 3D simulations (figure 13 of Zandian et al. Reference Zandian, Sirignano and Hussain2016, and figure 51 of the current article) – show KH wavelengths in the moderate range of $80$ – $125~\unicode[STIX]{x03BC}\text{m}$ over a wide range of $Re_{l}$ , $We_{g}$ and $\hat{\unicode[STIX]{x1D70C}}$ studied here. In order to expedite the appearance and growth of the KH waves, initial perturbations with a wavelength of $100~\unicode[STIX]{x03BC}\text{m}$ are imposed on the interface, with a small amplitude of $4~\unicode[STIX]{x03BC}\text{m}$ . The amplitude of the perturbations is small enough so that any other subharmonic wavelengths would have a chance to form and grow. As shown by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) and Zandian et al. (Reference Zandian, Sirignano and Hussain2017), at higher $Re_{l}$ and lower $\hat{\unicode[STIX]{x1D70C}}$ , a variety of smaller and larger wavelengths also form; at very low $We_{g}$ , the waves also merge to create longer wavelengths. Both streamwise ( $x$ -direction) and spanwise ( $y$ -direction) perturbations are considered in this study. The amplitude of the spanwise perturbations is reduced to $2.5~\unicode[STIX]{x03BC}\text{m}$ to let the other subharmonic waves appear naturally. At higher $Re_{l}$ , smaller spanwise waves (formed first as corrugations) appear superimposed on the initial perturbations (Zandian et al. Reference Zandian, Sirignano and Hussain2017). Similarly, the waves merge to create larger waves (lobes) at lower $We_{g}$ . Periodic boundary conditions for all components of velocity as well as the LS/VoF variable are imposed on the four sides of the computational domain; i.e. the $x$ - and $y$ -planes.

The most important dimensionless groups in this study are the Reynolds number ( $Re$ ), the Weber number ( $We$ ), the gas-to-liquid density ratio ( $\hat{\unicode[STIX]{x1D70C}}$ ) and viscosity ratio ( $\hat{\unicode[STIX]{x1D707}}$ ), as defined below. The initial KH wavelength-to-sheet-thickness ratio ( $\unicode[STIX]{x1D6EC}$ ) is also an important parameter that defines the relative length of the initial perturbations;

(2.12a-e ) $$\begin{eqnarray}\displaystyle Re=\frac{\unicode[STIX]{x1D70C}_{l}Uh}{\unicode[STIX]{x1D707}_{l}},\quad We=\frac{\unicode[STIX]{x1D70C}_{l}U^{2}h}{\unicode[STIX]{x1D70E}},\quad \hat{\unicode[STIX]{x1D70C}}=\frac{\unicode[STIX]{x1D70C}_{g}}{\unicode[STIX]{x1D70C}_{l}},\quad \hat{\unicode[STIX]{x1D707}}=\frac{\unicode[STIX]{x1D707}_{g}}{\unicode[STIX]{x1D707}_{l}},\quad \unicode[STIX]{x1D6EC}=\frac{\unicode[STIX]{x1D706}}{h}. & & \displaystyle\end{eqnarray}$$

The sheet thickness $h$ is considered as the characteristic length, and the relative gas–liquid velocity $U$ as the characteristic velocity. Theoretically, if the flow field is infinite in the streamwise direction (as in our study), a Galilean transformation shows that only the relative velocity between the two streams is consequential. Spatially developing flow fields, however, are at most semi-infinite so that both velocities at the flow-domain entry and their ratio (or their momentum ratio) are important. A wide range of $Re$ and $We$ at high and low density ratios is covered in this research. The main six cases studied in this article are presented in table 1. The cases are called in a ‘Dnx’ format, where ‘n’ presents the Domain number (1, 2 or 3), and ‘x’ presents the range of the case (‘a’ for high density ratio, and ‘b’ for low density ratios). Two cases – one at high density ratio and fairly high $We_{g}$ , and one at low $\hat{\unicode[STIX]{x1D70C}}$ and low $We_{g}$ – are studied in each domain to clearly show the effects of density ratio on the vortex dynamics of each atomization process. The ‘b’ cases are more practical and in the ranges usually seen in most atomization applications. The range of parameters for other cases presented in this article is indicated in their corresponding figure captions.

Table 1.

The main cases studied with their dimensionless parameter values.

2.3 Numerical methods

Direct numerical simulation (DNS) is done by using an unsteady 3D finite-volume solver to solve the Navier–Stokes equations for the planar incompressible liquid sheet segment (initially stagnant), which is subject to instabilities due to a gas stream that flows past it on both sides. The LS and VoF methods are used for liquid–gas interface tracking. The temporal deformation of the liquid–gas interface is predicted, resulting in 3D instabilities, that can lead to ligament formation and sheet breakup. The original code containing the LS subroutines was developed and applied by Dabiri, Sirignano & Joseph (Reference Dabiri, Sirignano and Joseph2007), and later by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014), Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016), and Zandian et al. (Reference Zandian, Sirignano and Hussain2016). Because of the weakness of the LS method in mass conservation at low density ratios, the VoF method is used at low gas densities.

A uniform staggered grid is used with a mesh size of $2.5~\unicode[STIX]{x03BC}\text{m}$ and a time step of 5 ns – a finer grid resolution of $1.25~\unicode[STIX]{x03BC}\text{m}$ is used for the case with higher $We_{g}$ ( $We_{g}=115\,000$ and $Re_{l}=320$ ) or higher $Re_{l}$ ( $Re_{l}=5000$ and $We_{g}=3000$ , $7250$ , $20\,000$ , and $36\,000$ ). The third-order-accurate QUICK scheme is used for spatial discretization and the Crank–Nicolson scheme for time marching. The velocity–pressure coupling is established using the SIMPLE algorithm. The initial velocity profile and some other initial estimations, e.g. initial surface perturbation wavelength, are obtained from a 2D planar full-jet simulation with the same schemes and methodology as those of Zandian et al. (Reference Zandian, Sirignano and Hussain2016).

The grid independency tests were performed previously by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014), Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016), and Zandian et al. (Reference Zandian, Sirignano and Hussain2016). They showed that the errors in the size of the ligaments, penetration length of the liquid jet, and the magnitude of the velocity computed using different mesh resolutions were within an acceptable range. The effects of the mesh size, the thickness of the fuzzy zone between the two phases, where properties have large gradients to approximate the discontinuities, and mass conservation of the LS method have been previously addressed by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014). The effects of mesh resolution on the most important flow parameters, e.g. surface structures, velocity and vorticity profiles, are studied in detail in § 2.4. The domain-size independency was also checked in both streamwise and spanwise directions to make sure that the resolved wavelengths were not affected by the domain length or width. The normal dimension of the domain was chosen such that the top and bottom boundaries remain far from the interface at all times, so that the surface deformation is not directly affected by the boundary conditions. In addition, accuracy tests and validation with experiments and other numerical approaches were performed previously (Dabiri et al. Reference Dabiri, Sirignano and Joseph2007; Jarrahbashi & Sirignano Reference Jarrahbashi and Sirignano2014), and will not be repeated here.

2.4 Effects of mesh resolution

The computational grid should be able to resolve well the boundary layer of the gas and liquid near the interface in order to capture the frequency and growth of the interfacial instabilities accurately (Fuster et al. Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013). More importantly, the requirement for sufficient numerical resolution to compute the formation of holes and corrugations on the lobes indeed is a stricter restriction on the mesh size in the current study. To assess the effects of grid resolution on the numerical results of the DNS simulation of the atomization, three grid resolutions are considered in this section: M1, M2 and M3, with $\unicode[STIX]{x1D6E5}=5~\unicode[STIX]{x03BC}\text{m}$ , $2.5~\unicode[STIX]{x03BC}\text{m}$ and $1.25~\unicode[STIX]{x03BC}\text{m}$ , respectively.

Figure 4 shows the mass conservation error in different mesh resolutions. The vertical axis measures the instantaneous volume of the liquid phase $V$ normalized by the initial volume of the liquid in the domain $V_{0}$ . The coarse grid M1 clearly has the largest mass loss, with an error of approximately 0.5 % near the end of the computation. The M2 and M3 grids, however, both have fair mass conservation, with low errors of $0.15$  % and ${<}0.1\,\%$ , respectively. The results shown in figure 4 are for a case with $\hat{\unicode[STIX]{x1D70C}}=0.05$ . At higher density ratios, the mass conservation is significantly better.

Figure 4.

Mass conservation check for different mesh resolutions.

Figure 5.

Close-up view of the liquid surface structures formed at the wave crest for different mesh resolutions. Hole formation in Domain II, for $Re_{l}=2500$ and $We_{g}=7250$ at $90~\unicode[STIX]{x03BC}\text{s}$ (a–c); corrugation formation in Domain III, for $Re_{l}=5000$ and $We_{g}=7250$ (Case D3a) at $90~\unicode[STIX]{x03BC}\text{s}$ (d–f); droplet formation in Domain II, for $Re_{l}=2500$ and $We_{g}=7250$ at $95~\unicode[STIX]{x03BC}\text{s}$ (g–h). Comparison of the relative sizes of holes, ligaments and droplets with the grid size for $Re_{l}=5000$ and $We_{g}=36\,000$ at $13~\unicode[STIX]{x03BC}\text{s}$ (i). $\hat{\unicode[STIX]{x1D70C}}=0.1$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ and $\unicode[STIX]{x1D6EC}=0.5$ for all images.

A close-up view of the different liquid surface structures found in different atomization domains are shown in figure 5 for different mesh resolutions of the same cases and at the same time steps. Figure 5(a–c) shows the effects of mesh resolution on depicting the hole formation in Domain II. Physically, holes are formed on the liquid lobes only when the lobe thickness is very small ( $O(10)~\text{nm}$ ) and the disjoining pressure becomes active (Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017). Here in the simulations, holes appear when the thickness of the liquid lobe (sheet) decreases to approximately the cell size $\unicode[STIX]{x1D6E5}$ . Since mechanisms of sheet rupture, such as disjoining pressure, are absent in the present study, and because of this numerical cutoff criterion, the initial perforation of the liquid sheets is highly dependent on the mesh resolution. As can be seen in figure 5(a–c), the M1 grid does not have enough resolution to capture the hole formation accurately. The holes form much sooner and expand much faster in the M1 grid compared to the other two finer grids. As seen in figure 5(c), the large perforation seen on the lobe is a result of merging of several smaller holes on the liquid sheet, which can be captured only by a high-resolution grid such as M3. The entire size of the hole is fairly similar for the M2 and M3 grids. The lobe rim is much thicker for finer resolutions. Consequently, a further increase of mesh resolution will only delay the pinch-off point but will not affect the ligaments formed from the expansion of the holes. Hole formation and rim dynamics are well captured by the M2 and M3 grids.

At high $Re_{l}$ , in Domain III of figure 1, corrugations form on the lobe rim. There are typically three or four such small-scale corrugations (per wavelength) on the lobe rims, each having a size (thickness) of approximately $15$ – $25~\unicode[STIX]{x03BC}\text{m}$ . Clearly, the M1 grid does not have enough resolution to capture these corrugations, and only the larger ligaments are resolved (figure 5 d). The M2 and M3 grids both resolve the corrugations correctly with very similar sizes and structures (figure 5 e,f). Figure 5(g,h) compares the droplet formation and the cascade of structures for the M2 and M3 grids. The general aspects of the process, e.g. size and location of the lobes, ligaments, and holes and thickness of the rims, are fairly similar for both grids. However, as described before, because of the numerical cutoff criterion, the holes form and expand slightly sooner ( ${<}1~\unicode[STIX]{x03BC}\text{s}$ ) for the M2 grid. The size of the resolved droplets and ligaments are also slightly smaller in the M2 grid. This is similar to what was seen by Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017) and Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) in their studies of the grid resolution. It is well known that the size of corrugations and the resulting ligaments and droplets decrease as $Re_{l}$ and $We_{g}$ increase. In the current study, the M3 grid is used only for those higher $Re_{l}$ and $We_{g}$ , and for the rest of the cases the M2 grid is adequate to capture the physics of the cascade process. As shown in figure 5(i) for a high $Re_{l}=5000$ and high $We_{g}=36\,000$ , the M3 grid resolution is fine enough to capture the smallest holes, ligaments and droplets. The radii of curvature of the corrugations and holes and ligaments can go down to approximately $6~\unicode[STIX]{x03BC}\text{m}$ , which can be resolved well only by the M3 grid with $\unicode[STIX]{x1D6E5}=1.25~\unicode[STIX]{x03BC}\text{m}$ .

Figure 6.

The streamwise velocity profile (a,c) and the spanwise vorticity profile (b,d) near the wave crest for Case D3a with different mesh resolutions at an early time $t=10~\unicode[STIX]{x03BC}\text{s}$ (a,b), and at a later time $t=70~\unicode[STIX]{x03BC}\text{s}$ (c,d). The gas–liquid interface locations are denoted by the dashed lines.

Since the major part of this study is concerned with the vortex dynamics of the cascade process, it is important that we resolve the velocity and vorticity in the gas–liquid mixing layer correctly. A comparison of the streamwise velocity and spanwise vorticity profiles near the wave crest at an early time ( $t=10~\unicode[STIX]{x03BC}\text{s}$ ) and a later time ( $t=70~\unicode[STIX]{x03BC}\text{s}$ ) are shown in figure 6. These profiles are obtained at the front-most tip of one of the lobe crests at two distinct times. All three grids resolve the gas and liquid boundary layers correctly at early times, with the M1 profile slightly shifted towards the sheet centre near the interface (figure 6 a). The gas vorticity layer $\unicode[STIX]{x1D6FF}_{g}$ imposed in this simulation is initially approximately $200~\unicode[STIX]{x03BC}\text{m}$ , and the velocity gradient, hence the vorticity magnitude near the interface increases with time (see figure 6 c). At later times (figure 6 c,d), the M1 grid fails to capture the velocity profile correctly and is also unable to resolve the small-scale velocity fluctuations and the small vortices near the interface. The M2 and M3 grids, however, resolve the velocity fluctuations similarly with very little difference, and vividly portray the location and amplitude of those fluctuations. The liquid sheet bulk velocity increases in the meantime.

In terms of calculating the vorticity, figure 6(b) shows that the magnitude of the spanwise vorticity obtained from M1 grid is significantly lower than that of M2 and M3 even at early times. The location of the vorticity peak is also closer to the sheet centre and inside the liquid phase because the velocity gradient and the vorticity magnitude cannot be predicted correctly with the M1 grid. At later times (figure 6 d), this difference between the vorticity magnitudes of M1 and M2/M3 becomes even more noticeable, where the M1 grid cannot resolve the small vortices. The M2 and M3 grids have much better consistency in resolving the number of peaks, their locations and their magnitudes; however, the M3 grid is required for capturing all the vortices correctly at such high $Re_{l}$ values. The spanwise vorticity is not equally distributed between the two phases, and it mainly sits on the gas side since the vorticity thickness (and the boundary layer thickness) is larger in the gas phase compared to the liquid phase. The inclination of the vortices towards the gas phase is more pronounced at lower density ratios (Hoepffner, Blumenthal & Zaleski Reference Hoepffner, Blumenthal and Zaleski2011).

Experiments and simulations show that as KH waves amplify, the convective velocity of the waves becomes very close to the Dimotakis speed (Dimotakis Reference Dimotakis1986) defined as $U_{D}=(U_{l}+\sqrt{\hat{\unicode[STIX]{x1D70C}}}U_{g})/(1+\sqrt{\hat{\unicode[STIX]{x1D70C}}})$ . Considering the case shown in figure 6(a) with $\hat{\unicode[STIX]{x1D70C}}=0.5$ and $U_{l}=20~\text{m}~\text{s}^{-1}$ , the Dimotakis speed gives $U_{D}\approx 53~\text{m}~\text{s}^{-1}$ . The interface velocity (at the base of the KH wave) in our simulation (figure 6 a) indicates a value of $U_{int}\approx 52~\text{m}~\text{s}^{-1}$ , which agrees well with the expected Dimotakis wave speed.

2.5 Data analysis

Our goal is to study the vorticity dynamics as well as the liquid surface dynamics in order to understand breakup mechanisms at different flow conditions. To this end, $\unicode[STIX]{x1D706}_{2}$ contours at different cross-sections of the domain are analysed in time. Here, we briefly review the definition of the $\unicode[STIX]{x1D706}_{2}$ method.

An objective definition of a vortex should permit the use of vortex dynamics concepts to identify coherent structures (CS), to explain formation and evolutionary dynamics of CS, and to explore the role of CS in turbulence phenomena. Jeong & Hussain (Reference Jeong and Hussain1995) define a vortex core as a connected region with two negative eigenvalues of $\unicode[STIX]{x1D64E}^{2}+\unicode[STIX]{x1D734}^{2}$ ; where, $\unicode[STIX]{x1D64E}$ and $\unicode[STIX]{x1D734}$ are the symmetric and antisymmetric components of $\unicode[STIX]{x1D735}u$ ; i.e.  $\unicode[STIX]{x1D61A}_{ij}=(u_{i,j}+u_{j,i})/2$ and $\unicode[STIX]{x1D6FA}_{ij}=(u_{i,j}-u_{j,i})/2$ . If $\unicode[STIX]{x1D706}_{1}$ , $\unicode[STIX]{x1D706}_{2}$ and $\unicode[STIX]{x1D706}_{3}$ are the eigenvalues such that $\unicode[STIX]{x1D706}_{1}\geqslant \unicode[STIX]{x1D706}_{2}\geqslant \unicode[STIX]{x1D706}_{3}$ , this definition is equivalent to the requirement that $\unicode[STIX]{x1D706}_{2}<0$ within the vortex core, since $\unicode[STIX]{x1D706}_{3}$ is always negative because the sum of the normal viscous stresses is zero. This definition is proven to meet the requirements for existence of a vortex core in different flow conditions (Jeong & Hussain Reference Jeong and Hussain1995), while the vortex identification by the Q-definition (Kolář Reference Kolář2007) may be incorrect when vortices are subjected to a strong external strain (Jeong & Hussain Reference Jeong and Hussain1995), as in our study.

3.1 Vortex dynamics and surface dynamics

To interpret the surface deformation via vorticity concepts and to more clearly delineate the complex 3D flow physics, the dual approaches of vortex dynamics and surface wave dynamics are employed. The language of wave dynamics focuses on crests, troughs and lobes, while the vortex dynamics terminology refers to vortex rollers (eddies), braids and hairpins. The instability starts from the initially symmetric KH waves. The crest of the KH waves, when amplified and slightly stretched in the flow direction, is referred to as ‘KH crest’. As the 3D character of instability develops, ‘KH crests’ divide into distinct ‘lobes’. The forward-most tip of the lobe is referred to as ‘spanwise crest’, and the region between two adjacent crests is called ‘spanwise trough’. In terms of the vorticity dynamics, the ‘KH crests’ and ‘KH troughs’ are equivalent to ‘vortex rollers’ and ‘braids’, respectively; see figure 7.

Figure 7.

Terminology used for vorticity and liquid–gas interface deformations from side view (a) and top view (b).

The spanwise vorticity rolls up as KH rollers. The braid connects two adjacent rollers separated in $x$ , which stretch the fluid in between. There are two modes of unstable waves, corresponding to the two surface waves oscillating exactly in or out of phase, commonly referred to as the sinuous (antisymmetric) and varicose (symmetric) modes, respectively; see figure 7.

3.2 Hole and bridge formations (LoHBrLiD mechanism)

Figure 8.

Liquid surface deformation in the LoHBrLiD mechanism; $Re_{l}=320$ , $We_{g}=115\,000$ ( $Oh_{m}=1.06$ ),  $\hat{\unicode[STIX]{x1D70C}}=0.5$ and $\hat{\unicode[STIX]{x1D707}}=0.0022$ , at $t=18~\unicode[STIX]{x03BC}\text{s}$ (a), $22~\unicode[STIX]{x03BC}\text{s}$ (b), $26~\unicode[STIX]{x03BC}\text{s}$ (c), $28~\unicode[STIX]{x03BC}\text{s}$ (d), $30~\unicode[STIX]{x03BC}\text{s}$ (e) and $32~\unicode[STIX]{x03BC}\text{s}$ (f).

The $LoHBrLiD$ mechanism occurs at medium $Re_{l}$ and high $We_{g}$ ; see figure 1. This process is shown in figure 8. The lobes form and thin on the primary KH wave crests. The middle section of the lobes (the braid), where the highest strain occurs, thins faster and thus perforates, creating a hole and a bridge on the lobe rim. Bridges become thinner as the holes expand. Finally, the bridges break and create one or two ligaments depending on the breakup location. The ligaments stretch and eventually break into droplets under capillary action.

At lower $We_{g}$ , the surface tension force resists the formation of holes. In this range, usually one hole forms on the liquid lobe and stretches into a large one. At higher $We_{g}$ , lobes thin much more easily since the resistance due to the surface tension forces is not large enough to stabilize the growth of the instabilities. In this range usually several small-scale holes appear at different locations on the lobe, as seen in the close-up view of figure 9(a). As these perforations grow, several of these holes merge to create larger holes (figure 9 b), and a thick liquid bridge is created on the lobe rim (figure 9 c). As the lobes continue to stretch and the holes continue to expand, the bridges become thinner (figure 9 d) and finally break to create several ligaments (figure 9 e). The mechanism of hole formation, to be discussed in this section, is believed to be the same at high and low $We_{g}$ values and density ratios; the only difference is in the time, location and size of the holes.

Figure 9.

Close-up view of the hole formation and expansion on a single liquid lobe; $Re_{l}=5000$ , $We_{g}=36\,000$ , $\hat{\unicode[STIX]{x1D70C}}=0.1$ and $\hat{\unicode[STIX]{x1D707}}=0.0066$ . $t=10~\unicode[STIX]{x03BC}\text{s}$ (a), $11~\unicode[STIX]{x03BC}\text{s}$ (b), $12~\unicode[STIX]{x03BC}\text{s}$ (c), $13~\unicode[STIX]{x03BC}\text{s}$ (d) and $14~\unicode[STIX]{x03BC}\text{s}$ (e).

Case D2a (table 1), which falls in the $LoHBrLiD$ class, is chosen in this section to study and explain the vortex dynamics in the hole-formation mechanism. At the end of this section, Case D2b is studied to clarify the effects of lowering density ratio. Zandian et al. (Reference Zandian, Sirignano and Hussain2016) and Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) related the hole formation to the overlapping of the hairpin vortices that form on the braids. Their finding is confirmed here and more details in the vortex overlapping process are revealed. In this section, $\unicode[STIX]{x1D706}_{2}$ contours are elucidated in two cross-sections in $x$ – $z$ planes – one passing through the spanwise crest and the other through the trough – along with the instantaneous liquid–gas interface (red lines) position at different times.

Figure 10.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest (left) and the top view of the liquid surface (right), at $t=6~\unicode[STIX]{x03BC}\text{s}$ (a), $10~\unicode[STIX]{x03BC}\text{s}$ (b) and $14~\unicode[STIX]{x03BC}\text{s}$ (c) of Case D2a.

Figure 10(a) shows the formation of lobes and hairpin vortices that occur on the braid on both top and bottom sheet surfaces, at very early time $t=6~\unicode[STIX]{x03BC}\text{s}$ . The vortex structure is symmetric at this time and consists of a large vortex just downstream of the wave, which is hereafter called the ‘KH vortex’ (shown by the white arrows), and hairpin vortices on the braid, between two adjacent KH rollers. The location of the vortices relative to the interface are the same on the spanwise crest and trough cross-sections at this time; i.e. the vortices on the spanwise trough cross-section (not shown here) are slightly upstream of those on the spanwise crest; see the top views of figure 10. The vortices have some undulations in $x$ and manifest a hairpin structure similar to what was observed by Bernal & Roshko (Reference Bernal and Roshko1986). Lasheras & Choi (Reference Lasheras and Choi1988) attributed the appearance of three-dimensionality to the stretching along the principal direction of the positive strain. The maximum amplification of the vortex lines occurs near the braid region, where positive strain is the maximum. Therefore, they called these hairpin vortices ‘strain-oriented vortex tubes’ and described a mechanism for evolution of the 3D instabilities, where vortices enhanced by stretching are pulled more strongly into the streamwise direction until a series of hairpins extend from the underside of one roller to the top of its neighbour. A similar mechanism is seen here for the two-phase flow.

As shown in figure 10(b), the upstream hairpins are pulled downstream by the KH vortex on the outer side of the roller, while the downstream hairpins are pulled upstream on the inner side of the KH roller (easily understood from a frame fixed with the KH vortex). The direction of the streamwise hairpin stretch is determined by the global induction of the KH roller, as indicated by the black arrows in figure 10(b). Martin & Meiburg (Reference Martin and Meiburg1991) explained that the sign of the streamwise vorticity component is determined by the competition between global and local inductions; i.e. between the overall effect of the vorticity field and the locally self-induced velocity of a vortex tube. While the global induction effect continues to determine the sign of the streamwise braid vorticity, the direction of the streamwise vorticity component in the crest region is determined by the local induction of the vortex tube (Martin & Meiburg Reference Martin and Meiburg1991).

Figure 11.

3D Schematics showing the overlapping of the two hairpin vortices – one from the lobe crest (outer black tube, pointing downstream), and the other from the braid (inner red tube, pointing upstream) (a); $A$ is the plane in which (b) is drawn; cross-sectional view of the $A$ -plane, showing the thinning of the lobe sheet due to the combined induction of the two oppositely oriented overlapping hairpins (b). The vortex schematics are periodic in $x$ - and $y$ -directions.

The motion of the upstream and downstream hairpins on the outer and inner sides of the KH roller – due to the roller’s induced motion – causes these two hairpins to align spanwise and overlap – one layer locating on the outer surface of the lobe, i.e. on the streamwise wave crest, and the other layer on the inner side of the lobe. This hairpin overlapping has been observed and explained in both mixing layers (see figures 3 through 6 of Comte et al. Reference Comte, Lesieur and Lamballais1992) and in two-phase round jets (see figure 9 of Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016). Comte et al. (Reference Comte, Lesieur and Lamballais1992) discuss helical pairing of the oppositely oriented hairpins and the formation of diamond-shaped vortex-lattice structures. As shown in the illustrative sketches of figure 11, and also described by Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016), the liquid sheet between a pair of these overlapping, oppositely pointed hairpins becomes thinner.

Figure 12.

$\unicode[STIX]{x1D706}_{2}$ isosurfaces in a top close-up view of a liquid lobe in Case D2a (a), and a 3D view of the same snapshot (b), at $t=16~\unicode[STIX]{x03BC}\text{s}$ . Top view of the box drawn in image (a), showing the overlapping of the red and green isosurfaces without the blockage of the grey KH isosurface and blue liquid surface (c). The isosurface values are: $\unicode[STIX]{x1D706}_{2}=-2\times 10^{10}~\text{s}^{-2}$ (grey), $-3\times 10^{10}~\text{s}^{-2}$ (green) and $-3\times 10^{9}~\text{s}^{-2}$ (red). The grey isosurface is made transparent in subfigures (a) and (b) to display the inner hairpin underneath it. The hatched area denotes the approximate hairpin overlapping region, where the hole would appear on the lobe later.

Figure 11(a) schematically depicts two overlapping hairpin vortices in the liquid–gas interface region – one outer vortex originating from the lobe crest and stretching downstream (the slender black tube), and the other inner vortex originating from the braid and stretching upstream (the red tube). The KH vortex is shown by the thicker grey tube in this figure. Figure 11(b) shows a cross-sectional view of the vortex structure along with the lobe located between the two hairpin vortices on the $A$ -plane of figure 11(a). The induced velocity of the two oppositely oriented overlapping hairpin vortices (see the qualitative streamlines shown by the black and red arrows in figure 11 b) pushes the top surface of the lobe downward and the bottom surface upward, causing the lobe to become thinner in the middle and thus become vulnerable to puncture at that region. This occurs at $t=16~\unicode[STIX]{x03BC}\text{s}$ , on the top surface of the liquid sheet, as shown in figure 12. Three $\unicode[STIX]{x1D706}_{2}$ isosurfaces are shown in this figure: the grey isosurface denotes the KH vortex (the grey tube in figure 11); the green isosurface denotes the outer hairpin (the black hairpin in figure 11); and the red isosurface denotes the inner hairpin (the red hairpin in figure 11). The $\unicode[STIX]{x1D706}_{2}$ isosurface tracks the strength of the vorticity but circulation on its surface is not constant (as it is not a vortex surface); thus, it slightly differs from the surface of a vortex tube. The $\unicode[STIX]{x1D706}_{2}$ magnitude of each isosurface is denoted in the figure caption. The KH vortex in figure 12 is the strongest (i.e. greatest circulation) and, at its core, has larger $\unicode[STIX]{x1D706}_{2}$ values than the hairpins; however, the outer isosurface, with a lower $\unicode[STIX]{x1D706}_{2}$ value, is depicted to give a better indication of its size. The hairpin overlapping region is denoted by a hatched area in figure 12(b). The overlapping of the outer and inner hairpins can be clearly observed in figure 12(c), where the grey KH vortex isosurface and the blue liquid surface have been removed from the box drawn in figure 12(a). The hatched area shows the zone where the lobe thinning occurs and the first hole forms a few microseconds later; see figures 14 and 15(a).

Figure 13.

Schematic showing the transition of a symmetric thin sheet (a) to antisymmetric mode (c). The dashed lines denote the sheet centre plane. The black and dark grey circles denote the two adjacent KH vortices on the top surface, and the light grey and white circles are two adjacent KH vortices on the bottom surface.

The sketches in figure 13 show the transition of a thin sheet from a symmetric mode to an antisymmetric mode with the qualitative location of the KH vortices at the two liquid surfaces. In the beginning, the top and bottom sheet surfaces as well as their KH vortices are symmetric with respect to the centre plane, as shown in figure 13(a). All plane jets in reality exit with symmetric perturbations because of the long-wavelength perturbations due to the upstream chamber’s Helmholtz resonance modes or driving compressor blade wakes, but the antisymmetric mode has a much higher growth rate and thus eventually dominates. The true physical explanation of this transition is intriguing, but remains somewhat elusive, and deserves careful examination and explanation. Ashurst & Meiburg (Reference Ashurst and Meiburg1988) showed that the loss of symmetry in the two-vortex-layer calculation is related to the enhancement and reduction of relative streamwise displacements for filaments in the stronger layer by the addition of the second layer of vorticity.

Due to this transition towards antisymmetry, the KH rollers on one of the layers (the bottom layer in our case) are not able to stretch the bottom lobes and cause the hairpins overlapping at early times; thus, lobe formation and thinning are delayed on the bottom surface. This corresponds to the instant in figure 13(b) and the period $18$ – $34~\unicode[STIX]{x03BC}\text{s}$ in our simulation, shown in figures 14 through 19. As the sheet gets thicker, i.e. the two vorticity layers get more independent, the transition towards antisymmetry is delayed, and both sides have time to roll up the surface waves and stretch them into lobes. For the rest of this section, our focus will be only on the top interface, where lobe stretching and hole formation are evident at early times.

Figure 14.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest (a) and the spanwise trough (b), with the top view of the liquid surface on the right, at $t=18~\unicode[STIX]{x03BC}\text{s}$ of Case D2a.

As figure 14 shows, the KH vortex starts to move away from the KH wave at $18~\unicode[STIX]{x03BC}\text{s}$ . However, it has had enough time to roll the lobe and make the overlapping occur on the top sheet surface. Notice that the KH vortex on the trough cross-section (figure 14 b) is still closer to its original location compared to the same vortex on the crest cross-section (figure 14 a). As described earlier, this results in further stretching of the hairpins near the sides of the lobe compared to its crest; hence, holes occur sooner near the lobe sides than at the lobe crest. This motivates the idea that the KH vortex location and strength are important in formation or inhibition of holes at early times. The KH vortex transition and structure is directly related to the liquid viscosity, and thus the Reynolds number. At high $Re_{l}$ , even though vorticity diffusion is slower, the fluids motion is less constrained by the viscous forces; thus, the KH vortices advect away from the interface much more easily, and do not acquire enough time to roll up the lobes and stretch the hairpins. Hole formation also depends on the surface tension. The inertia should have enough strength to overcome the surface tension forces in order to perforate the lobe. In conclusion, the $LoHBrLiD$ mechanism becomes less probable as $We_{g}$ decreases or $Re_{l}$ becomes very high or very low; see figure 1.

Figure 15.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest (left) and the top view of the liquid surface (right), at $t=20~\unicode[STIX]{x03BC}\text{s}$ (a) and $24~\unicode[STIX]{x03BC}\text{s}$ (b) of Case D2a.

$\unicode[STIX]{x1D706}_{2}$ contours of figure 15 clearly show the overlapping of the outer and inner hairpins at the location of the holes. The holes expand and create bridges at $t=20~\unicode[STIX]{x03BC}\text{s}$ (figure 15 a). Meanwhile, the KH vortices advect both downstream and away from the interface. Later at $24~\unicode[STIX]{x03BC}\text{s}$ (figure 15 b), the bridges break and create spanwise ligaments. If the bridge breakup were to have happened at the spanwise crest instead of its trough, the resulting two ligaments would have been streamwise-oriented, as in the sketch of figure 2.

Figure 16.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest (a) and the spanwise trough (b), with the top view of the liquid surface on the right, at $t=28~\unicode[STIX]{x03BC}\text{s}$ of Case D2a.

The broken bridges undergo capillary instability and break further into droplets and smaller ligaments, as shown in figures 16 and 19. Meanwhile, the hairpin vortices that are overlapping near the spanwise trough, thin the bridges, as has been indicated by the red arrows in figure 16(b). This thinning along with the stretching of the bridges by the induced velocity of the KH roller deforms the bridge and makes it thinner and on the verge of further breakup. This completes the $LoHBrLiD$ breakup mechanism at primary atomization. The liquid interface at this time has become completely antisymmetric (figure 16) following the KH vortices that became antisymmetric earlier.

The role of vortices in the $LoHBrLiD$ breakup mechanism is summarized schematically in figure 17. This figure shows the liquid surface and also the qualitative location of the nearby vortices at four consecutive times. At an early time $t_{1}$ , hairpin vortices form on the braids due to the strain caused by the neighbouring KH vortices. The hairpins closer to the KH wave crest – shown by black lines – are stretched downstream, and the hairpins near the KH wave trough – shown by red lines – are stretched upstream by the induced motion of the KH rollers. The hairpin parts that are stretched downstream are rolled over the KH vortex tube, and are denoted by solid lines, while the hairpin parts that are stretched upstream are pulled under the KH vortex tube, and are denoted by dashed lines.

Figure 17.

Schematics of the $LoHBrLiD$ process at four consecutive times. The liquid–gas interface is shown in blue, and the KH vortex by black tubes. The red and black lines denote the inner and outer hairpin vortices near the KH wave trough and crest, respectively. The solid or dashed lines denote where the hairpins are stretched upstream and inward, or downstream and outward, respectively.

Later at $t_{2}$ , the KH vortices deflect more under the induction of the hairpin vortices. The KH vortex stretches the lobes over itself and covers underneath the lobe surface at later times (see times $t_{3}$ and $t_{4}$ in figure 17). The black and red hairpins that roll over and under the KH vortex, respectively, overlap later at the centre of the lobe as well as the two sides of the lobe; i.e. at the spanwise troughs. These overlapping regions have been denoted in figure 17 at $t_{3}$ . The lobe sheet fills the vertical gap between these overlapping hairpins, causing the lobe to become thinner and thus vulnerable to puncture at the overlapping regions, as was described in figure 11.

Whether the liquid sheet subject to these conditions punctures or not depends on other flow conditions, particularly the surface tension. At high $We_{g}$ (high $Oh_{m}$ ), the hole formation prevails and the lobes perforate at the predicted locations, as shown at $t_{4}$ in figure 17. As the overlapping hairpin filaments continue to stretch, the holes also stretch and expand, creating even larger holes and thinner bridges. If $We_{g}$ is not large enough, the liquid lobe in the overlapping region can recover instead. In this case, hole formation is inhibited and the lobes stretch directly into ligaments via $LoLiD$ or $LoCLiD$ mechanisms, as discussed later.

Figure 18.

$\unicode[STIX]{x1D706}_{2}$ isosurface in a top close-up view of a liquid lobe in Domain II at low density ratio (Case D2b) at $90~\unicode[STIX]{x03BC}\text{s}$ (a); the solid black line shows the lobe front edge location and the hatched area shows the location of hairpin overlapping. The isosurfaces represent: the KH vortex with $\unicode[STIX]{x1D706}_{2}=-4\times 10^{9}~\text{s}^{-2}$ (grey), the outer crest hairpin with $\unicode[STIX]{x1D706}_{2}=-6\times 10^{9}~\text{s}^{-2}$ (green) and the inner trough hairpin with $\unicode[STIX]{x1D706}_{2}=-3\times 10^{9}~\text{s}^{-2}$ (red). Lobe surface showing the hole formation from a top view at $90~\unicode[STIX]{x03BC}\text{s}$ (b), $92~\unicode[STIX]{x03BC}\text{s}$ (c) and $94~\unicode[STIX]{x03BC}\text{s}$ (d) of Case D2b.

The mechanism of hole formation at low density ratios is similar to the high $\hat{\unicode[STIX]{x1D70C}}$ described in figure 17, with minor differences. Figure 18 shows the vortex structures and the lobe deformation at $90$ – $94~\unicode[STIX]{x03BC}\text{s}$ for $\hat{\unicode[STIX]{x1D70C}}=0.05$ (Case D2b). Three $\unicode[STIX]{x1D706}_{2}$ isosurfaces corresponding to the KH vortex (grey), the outer hairpin (green) and the inner hairpin (red) are shown in figure 18(a) at $90~\unicode[STIX]{x03BC}\text{s}$ . The liquid surface has been removed from this image to reveal all the vortex structures, but the black solid line indicates the location of liquid lobe front edge for comparison. The liquid surface at the same time is shown in figure 18(b). The outer and inner hairpins at this low $\hat{\unicode[STIX]{x1D70C}}$ stretch and wrap around the KH vortex in opposite directions, and overlap on top and bottom of the lobe. The difference here is that the hairpins are more stretched and the overlapping region is more elongated in the streamwise direction, as shown by the hatched area in figure 18(a,b). The holes that are formed near the hatched area at later times (figure 18 c,d) are also more elongated in the streamwise direction compared to the high $\hat{\unicode[STIX]{x1D70C}}$ case shown in figure 15, where the holes were more spanwise oriented. This finding agrees with the predictions of hole location by Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) (see their figure 20). There when the gas density increased, the edge curvature decreased and the locations of the holes in neighbouring lobes became closer to each other. This caused the neighbouring holes to merge and expand the hole in the azimuthal (spanwise) direction.

The breakup mechanism involving hole and bridge formation and the role of vortical structures in this regard were discussed in this section via consecutive tracking of $\unicode[STIX]{x1D706}_{2}$ contours near the interface. Vortex dynamics is able to explain the breakup process at high $We_{g}$ and medium $Re_{l}$ . In the remainder of this section, $\unicode[STIX]{x1D706}_{2}$ contours along with the interface location are shown at a few later time steps to track the deformation of the ligaments and also to see the evolution of the surface waves after the KH vortices have left the interface.

Figure 19.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest (left) and the top view of the liquid surface (right), at $t=30~\unicode[STIX]{x03BC}\text{s}$ (a), and $34~\unicode[STIX]{x03BC}\text{s}$ (b) of Case D2a.

Figure 19 shows $\unicode[STIX]{x1D706}_{2}$ contours of the liquid jet in the period $30$ – $34~\unicode[STIX]{x03BC}\text{s}$ . Even though the original KH vortex (denoted by the white arrow) is now advected downstream and away from the waves, the vorticity contours on the braid stretch and collect at the crest of the new KH waves and create new KH vortices (figure 19 a). These vortices roll a new KH wave and stretch new sets of lobes similar to the original ones; see figure 19(b). The direction of the induced motion by these new vortices are indicated by the black arrows in figure 19(b). The direction in which the fluids are swirling around the KH vortices are also denoted by the blue curly arrows. At $t=34~\unicode[STIX]{x03BC}\text{s}$ the sheet and the vortices are antisymmetric with half a wavelength phase difference between the top and bottom surfaces in the streamwise direction. The vortices on the bottom side are now nearly stationary with respect to the interface; hence, they have enough time to roll the waves and form stretched lobes on the bottom surface (figure 19 b). As the new lobes get stretched, new pairs of hairpins form on their braids, and the whole $LoHBrLiD$ process repeats, creating new holes, bridges, ligaments and droplets. Meanwhile, the formerly broken bridges undergo capillary instability and thin at the necks and break into droplets, as shown in the top view of figure 19(b).

Overlapping of the oppositely oriented hairpins is the cause of lobe perforation, as described in this section; however, the locations of the perforations and the direction and rate of their growth depend on the overlapping location with respect to the lobe sheet. What was presented in this section was an ideally simple case to explain the process. In reality, the hairpin structures can be more irregular than what was presented here, which would shift the overlapping zones and hence change the location of the first perforations. More hole-formation examples are presented by Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) and Zandian et al. (Reference Zandian, Sirignano and Hussain2016). Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) also showed the direction in which the holes merge and expand at high and low density ratios.

3.3 Corrugation formation at high $Re_{l}$ ( $LoCLiD$ mechanism)

The $LoCLiD$ mechanism occurs at high $Re_{l}$ and low $We_{g}$ (low $Oh_{m}$ ), as indicated in figure 1. This process is shown in figure 20. The lobes form similar to the previous case, but do not stretch as much. Corrugations form on the lobes’ front edge and stretch to create ligaments. Multiple ligaments are formed per lobe, typically shorter and thinner compared to the ligaments seen in the Domain II for $LoHBrLiD$ . Eventually, the ligaments detach from the liquid jet and break up into droplets by capillary action. These droplets are consequently smaller than the ones formed in the $LoHBrLiD$ mechanism. Sequential evolution of $\unicode[STIX]{x1D706}_{2}$ contours on $y$ -plane cross-sections are used in this section to delineate the physics of this process, and more particularly describe why lobes do not perforate but get corrugated at higher $Re_{l}$ and lower $We_{g}$ . Case D3a (see table 1) is analysed in this section; at the end of this section, the effects of lower density ratio are described by analysing Case D3b.

Figure 20.

Liquid surface deformation following the $LoCLiD$ mechanism in Case D3a, at $t=44~\unicode[STIX]{x03BC}\text{s}$ (a), $48~\unicode[STIX]{x03BC}\text{s}$ (b), $50~\unicode[STIX]{x03BC}\text{s}$ (c),  $52~\unicode[STIX]{x03BC}\text{s}$ (d), $56~\unicode[STIX]{x03BC}\text{s}$ (e) and $60~\unicode[STIX]{x03BC}\text{s}$ (f).

Figure 21.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest cross-section, at $t=6~\unicode[STIX]{x03BC}\text{s}$ (a) and $10~\unicode[STIX]{x03BC}\text{s}$ (b) of Case D3a.

Figure 21 shows $\unicode[STIX]{x1D706}_{2}$ contours on a $y$ -plane passing through the spanwise crest at $6~\unicode[STIX]{x03BC}\text{s}$ and $10~\unicode[STIX]{x03BC}\text{s}$ . Since the vortex structure passing through the spanwise trough has similar features, albeit with some phase difference with respect to the vortices at the spanwise crest, it will not be shown here. In the beginning, the vorticity field contains a series of hairpin filaments on the braid and a stronger KH vortex just downstream of the KH wave crest (indicated by the white arrows). The hairpins are stretched by the neighbouring rollers as discussed before. As also described by Martin & Meiburg (Reference Martin and Meiburg1991), since stretching is the most intense in the downstream half of the braid region, hairpin vortices acquire a larger streamwise component in the upstream neighbourhood of a vortex roller (i.e. on the lobe crest) than its downstream side (i.e. on the trough). Hence, the hairpin filaments that are closer to the streamwise wave crest (denoted by the grey arrows) are stronger than the upstream braid hairpins. Ashurst & Meiburg (Reference Ashurst and Meiburg1988) also found that, while in one vortex layer the filaments at the centre between two spanwise rollers experience the most stretching, if a second vortex layer is added, the filaments with the strongest spanwise modulation will be located closer to the downstream roller; i.e. three-dimensionality occurs first in the downstream half of the braid region. The vorticity field as well as the sheet itself are initially symmetric with respect to the centre plane. The hairpin filaments reside near the interface, slightly inclined towards the gas zone.

Figure 22.

Schematic of how KH vortex splits into a slow-moving vortex (denoted by $s$ ) in the liquid and a fast-moving vortex (denoted by $f$ ) in the gas zone, at four consecutive times. The qualitative velocity magnitudes are denoted by the straight arrows in the gas and the liquid. The vortices are nearly uniform in the $y$ direction, normal to the paper.

Figure 23.

$\unicode[STIX]{x1D706}_{2}$ isosurfaces in a close-up view of a liquid lobe in Case D3a, from a top view (a) and a 3D view (b), at $t=12~\unicode[STIX]{x03BC}\text{s}$ . The isosurface values are: $\unicode[STIX]{x1D706}_{2}=-8\times 10^{10}~\text{s}^{-2}$ (grey), $-4\times 10^{10}~\text{s}^{-2}$ (green) and $-10^{10}~\text{s}^{-2}$ (red). The liquid surface is shown in blue.

At high $Re_{l}$ , inertia dominates the viscous forces and the fluid particles are more free to move around as the viscous forces do not produce enough resistance against their motion – especially in the gas zone. The higher velocity of the gas layer compared to the liquid layer causes the KH roller to split into two vortices at $10~\unicode[STIX]{x03BC}\text{s}$ (see figure 21 b); this process is shown schematically in figure 22 and using the $\unicode[STIX]{x1D706}_{2}$ isosurfaces in figure 23. The outer part of the KH vortex (indicated by the white arrow in figure 21 b), which resides in the fast-moving gas layer, separates from the part that is inside the surface of the lower speed liquid (indicated by the black arrow in figure 21 b). The splitting of the KH vortices starts at $10~\unicode[STIX]{x03BC}\text{s}$ (figure 21 b) and continues until $12~\unicode[STIX]{x03BC}\text{s}$ (figure 23). The grey, green and red $\unicode[STIX]{x1D706}_{2}$ isosurfaces in figure 23 denote the KH vortex, the outer (crest) hairpin and the inner (braid) hairpin, respectively (the magnitude of each $\unicode[STIX]{x1D706}_{2}$ isosurface is indicated in the caption). The split slow- and fast-moving vortices are denoted by ‘ $s$ ’ and ‘ $f$ ’ in figure 22, respectively.

Figure 24.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest cross-section, at $t=14~\unicode[STIX]{x03BC}\text{s}$ (a), $16~\unicode[STIX]{x03BC}\text{s}$ (b) and $18~\unicode[STIX]{x03BC}\text{s}$ (c) of Case D3a.

As demonstrated in figures 22 and 24(a), part of the KH vortex that resides in the faster-moving gas, near the liquid interface, advects downstream with the gas, while the slow-moving vortex (indicated by the black arrow in figure 24 a) advects more slowly with the interface velocity, remaining stationary relative to the liquid surface. The two vortices are completely separated at $t_{4}$ in figure 22, corresponding to $t=16~\unicode[STIX]{x03BC}\text{s}$ in our simulation (see figure 24 b). This vortex separation has two significant consequences. (i) The slow-moving vortex (black vortex) downstream of the KH wave is not strong enough to curl the KH wave and stretch the lobe downstream over itself. Consequently, the outer hairpins do not overlap in $x$ with the inner trough hairpins, as in the Domain II for $LoHBrLiD$ ; hence, the hole formation is inhibited at early times. (ii) The fast-moving vortex (the white vortex) gets closer to the downstream outer hairpin as it advects away from the upstream hairpins, residing in the vicinity of the trough; see figure 24(b). This successive variation in the distance between the hairpins and the fast-moving vortex induces a fluid motion that stretches the hairpins in the opposite directions, resulting in a less-orderly hairpin structure with more undulations, as shown schematically in figure 25.

Figure 25.

Schematic of the hairpin undulation formation under the influence of the fast-moving KH vortex (denoted by 1 $f$ ) at four consecutive times $t_{1}$ – $t_{4}$ , from a frame of reference moving with the slow-moving KH vortices (denoted by 1 $s$ and 2 $s$ ).

After splitting of the KH vortex into fast- and slow-moving vortices ( $t_{2}$ in figure 25), the fast-moving vortex (1 $f$ ) advects downstream over the outer hairpin ( $t_{3}$ in figure 25). As the faster-moving vortex (1 $f$ ) passes over the hairpin, its induced motion pulls the crest of the hairpin in the upstream direction under vortex 1 $f$ , causing the hairpin to undergo an undulation with smaller local wavelength ( ${\approx}45~\unicode[STIX]{x03BC}\text{m}$ ). Meanwhile, the tip of the hairpin (now having two crests) is stretched downstream over the neighbouring slow-moving vortex (2 $s$ ), and the trough of the hairpin is stretched upstream under the upstream slow-moving vortex (1 $s$ ) ( $t_{3}$ in figure 25). The streamwise stretching on different parts of the hairpins is indicated by the black curly arrows. As the faster-moving vortex (1 $f$ ) moves further downstream, the newly formed undulation wraps around it and stretches downstream, as shown at $t_{4}$ in figure 25. Therefore, another turn is created on the hairpin vortex and the local wavelength of the undulations decreases ( $15$ – $35~\unicode[STIX]{x03BC}\text{m}$ ). As the fast-moving vortex keeps moving downstream, the hairpin corrugations stretch with it. It will be shown below that this hairpin structure results in less stretched (slower stretching) and more corrugated lobes. The direction of the streamwise stretch on the hairpin filaments – based on the velocity field induced by the two split rollers – is denoted by the red arrows in figure 24(c). The sheet is still symmetric at this moment, while the vortices have lost their symmetry with respect to the centre plane. The fast-moving vortex (1 $f$ ) finally reaches the downstream slower-moving vortex (2 $s$ ) and combines with it to create a stronger KH vortex, which stretches the now corrugated hairpin in the downstream direction.

Figure 26.

Schematic of the $LoCLiD$ process at four consecutive times. The liquid surface is shown in blue, and the KH vortex is denoted by the black tube. The red and black lines denote the inner and outer hairpin vortices near the trough and crest of the KH wave, respectively. The solid and dashed lines indicate the parts of the hairpins that are stretched upstream and inward, or downstream and outward, respectively.

To better understand the consequence of hairpin vortex structure on the liquid surface deformation at high $Re_{l}$ and low $We_{g}$ , the evolution of the vortices and the interface in the $LoCLiD$ process (Domain III) are schematically depicted in figure 26 at four consecutive instances. At an early time $t_{1}$ , the braid regions connecting the emerging KH rollers become progressively more depleted of vorticity. The spanwise vortices on the braid deflect due to the induced motion of the neighbouring KH rollers in both the upstream and downstream directions – creating the hairpin vortex structures with a spanwise size equal to the spanwise perturbation wavelength ( $100~\unicode[STIX]{x03BC}\text{m}$ ). Two hairpins are formed on the braid – one located near the lobe crest, called the outer hairpin and denoted by the black line (corresponding to the green isosurface in figure 23), and the other slightly upstream near the trough, called the inner hairpin and denoted by the red line (corresponding to the red isosurface in figure 23). The deflected hairpin filaments form the lobes as they are stretched by the KH roller. So far, the process is similar to the $LoHBrLiD$ process.

As discussed earlier and demonstrated in figure 25, the hairpin vortices get corrugated under the influence of the fast-moving vortex that hovers over the interface after splitting of the KH vortex. This corresponds to $t_{2}$ and $t_{3}$ in figure 26. Only the deformation of the outer hairpin (black hairpin) was shown in figure 25; the inner hairpin (red hairpin) undergoes the same process, but in the opposite direction – stretching downstream and over the fast-moving vortex. The upstream turns and bends on the hairpins prevents further downstream stretching of the lobes. Consequently, the lobes are less stretched and more blunt compared to the lobes in Domain II (compare $t_{2}$ in figure 26 with $t_{2}$ in figure 17). The fast-moving vortex moves closer to the interface to merge with the downstream slow-moving vortex and later passes through the lobe and reaches the slower vortex below the lobe. In this process, the faster vortex carries along the corrugated outer hairpin that is wrapping around it. Consequently, the corrugated outer hairpins penetrate through the lobe and locate underneath the lobe in contrast to Domain II, where the outer hairpins were located on top of the lobe. Since this process occurs quickly, the corrugated hairpins do not have enough time to completely wrap around the faster vortex and create a combined vortex; therefore, after the merging of the two split KH vortex counterparts, only two pairs of undulated hairpins – which are separate from the KH vortex – remain on top and bottom of the reformed KH vortex, right below the lobe. The liquid surface approximately follows the hairpin structures with some delay at this high- $Re_{l}$ range – as the vortex lines are nearly material lines. Because of these shorter hairpin wavelengths, corrugations with length scales comparable to the local hairpin wavelengths ( $15$ – $25~\unicode[STIX]{x03BC}\text{m}$ ) form on the front-most edge of the lobes, as shown in figure 26 at $t_{3}$ . Both experimental observations (Lasheras & Choi Reference Lasheras and Choi1988; Liepmann & Gharib Reference Liepmann and Gharib1992) and numerical simulations (Brancher et al. Reference Brancher, Chomaz and Huerre1994; Danaila et al. Reference Danaila, Dušek and Anselmet1997) for homogeneous jets show that the size of vortex pairs (lobes in two-phase flows) decreases with increasing Reynolds number.

Upon creation of a stronger KH vortex downstream of the KH waves at $t_{4}$ – after merging of the fast- and slow-moving vortices – the new KH roller, which is now located under the lobe (now a thicker tube) is strong enough to stretch the hairpins and the corrugations. The two hairpin layers overlap as illustrated in figure 26 at $t_{4}$ . In this figure, the dashed lines represent the hairpins stretching upstream and on the streamwise trough (i.e. lower surface of the gas tongue), while the solid lines denote parts of the hairpins that are stretched downstream below the lobe (i.e. upper surface of the gas tongue).

Figure 27.

3D Schematics showing the overlapping of the outer hairpin (black slender tube) and the inner hairpin (red tube) resulting in formation of lobe corrugations (a) – $A$ is the plane in which (b) is drawn; cross-sectional view of the $A$ -plane showing the corrugation formation and thinning of the gas tongue due to the combined induction of the overlapping hairpins (b). The vortex schematics are periodic in $x$ - and $y$ -directions.

The illustrative sketches of figure 27 show the corrugated hairpin structures and their position with respect to the lobe in the $LoCLiD$ process (corresponding to $t_{4}$ in figure 26). The black tube represents an outer hairpin on the lobe crest and the red tube represents an inner hairpin on the streamwise trough (originating from the braid). The two hairpin layers – after getting corrugated due to the induced motion of the two halves of the split KH vortex (see $t_{2}$ and $t_{3}$ in figure 25) – overlap under the lobe, as shown in figure 27(a). Following the fluid motion induced by the vortex pairs, the lobe front edge gets corrugated with comparable length scales to the hairpin undulation wavelength; see figure 27(a). The layer of the outer hairpin (black tube) is located underneath the lobe (on top of the gas layer), and the layer of the inner hairpin (red tube) is located on top of the interface at the trough; see figure 27(b). The induced flow creates undulations on both the bottom surface of the lobe and the trough surface. The combined induction of the two oppositely oriented overlapping hairpin layers thins the gas layer that fills the vertical gap underneath the lobe; i.e. the lobe collapses on the liquid sheet as the bottom surface of the lobe descends and the trough surface ascends; see the qualitative streamlines in figure 27(b).

As the counter-rotating pairs of hairpins stretch under the induction of the KH vortex, the corrugations on the lobes stretch with them and form thin ligaments (see $t_{4}$ in figure 26). The ligaments stretch downstream and break up into droplets as they undergo capillary instabilities (the ligament breakup mechanism is discussed in § 3.4). In the meantime, the eddies cascade into smaller vortical structures as transition to turbulence occurs. The smaller vortices are moved towards the gas by the induction of the larger eddies. The outward movement of the vortices spreads the droplets in the normal direction, helps the expansion of the spray angle, and enhances the two-phase mixing.

The process of lobe-sheet collision after the merging of the two split KH vortices is shown in the side views of the liquid sheet along with $\unicode[STIX]{x1D706}_{2}$ contours near the interface in figure 28. The antisymmetric vortices, which now are strong enough to stretch the lobes, create antisymmetric KH waves, at $44~\unicode[STIX]{x03BC}\text{s}$ (figure 28 a). The hairpin vortices on the KH wave crest and trough stretch in the opposite directions following the induced fluid motion of the KH roller. Following the direction of the swirl, the hairpin projections on the trough are stretched upstream from the inner side of the roller, while the hairpin projections on the KH crest are stretched downstream from the outer side of the roller. The direction of streamwise stretch on the fluid elements is shown by the red arrows in figure 28.

Figure 28.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest cross-section, at $t=44~\unicode[STIX]{x03BC}\text{s}$ (a), $46~\unicode[STIX]{x03BC}\text{s}$ (b) and $48~\unicode[STIX]{x03BC}\text{s}$ (c) of Case D3a.

Recall that at lower $Re_{l}$ the two overlapping hairpins lie above and below the lobe and cause its thinning. At higher $Re_{l}$ , on the other hand, the outer and inner hairpins are inside the liquid (very close to the surface), and overlap on outer and inner sides of the gas layer that penetrates under the lobe. To understand this better, consider the opposite perspective; i.e. the layer of gas that fills the vertical gap underneath the KH waves can be considered as an upstream-pointed gas tongue, surrounded by the liquid; see figure 28(b). The combined induction of the overlapping hairpins thins the tongue similarly to the process described in figure 11. The thinning regions are indicated by the two thin black arrows in figure 28(c). Since it would be confusing to describe a hole formation in the tongue, we could interpret this phenomenon as the collapse of the lobe on the jet, encapsulating the gas in between.

Figure 29.

Top view of a liquid lobe in Case D3a at $t=50~\unicode[STIX]{x03BC}\text{s}$ (left) and the $\unicode[STIX]{x1D706}_{2}$ contours on three $yz$ -cross-sections (right); $A$ – $A$ passing upstream of the lobe front, $B$ – $B$ passing through the corrugations, and $C$ – $C$ passing downstream of the lobe.

Figure 29 shows the top view of a corrugated lobe and the $\unicode[STIX]{x1D706}_{2}$ contours at three spanwise ( $y$ – $z$ ) planes passing though the lobe body (plane $A$ – $A$ ), through lobe corrugations (plane $B$ – $B$ ), and through the braid (plane $C$ – $C$ ), at $50~\unicode[STIX]{x03BC}\text{s}$ . The $\unicode[STIX]{x1D706}_{2}$ contours clearly show the vortical structures suggested by the schematics of figure 27. Upstream of the lobe front, at plane $A$ – $A$ , one downstream stretching vortex pair is seen just underneath the lobe – indicated by the black curly arrows – and three upstream stretching counter-rotating vortex pairs are seen just above the trough surface – indicated by the red curly arrows. At this cross-section, the red hairpin shows the three vortex pairs shown in figure 27(b), but the black hairpin still has the main two legs and does not manifest undulations in this plane. That is, the undulations have not stretched enough to reach the $A$ – $A$ plane yet. The effects of these vortex pairs on the liquid surface is evident in this figure. Moving slightly downstream, in plane $B$ – $B$ , both black and red hairpins manifest three vortex pairs, located under the lobe and above the trough, respectively. The $\unicode[STIX]{x1D706}_{2}$ contours at this cross-section agree well with the scenario depicted in figure 27(b) on plane $A$ . However, the vortex pairs are not nicely aligned in $y$ , as shown in figure 27(b). The distance between the two legs of each vortex pair has been denoted on the figure. As indicated before, the vortex undulations with wavelengths in the range $15$ – $35~\unicode[STIX]{x03BC}\text{m}$ are seen in this picture. The creation of bumps and craters on the liquid surface induced by these vortex pairs is also seen here. Further downstream, on plane $C$ – $C$ , only the two main legs of the red hairpin and parts of the KH vortex are seen. This is what would be expected based on the 3D schematic of figure 27(a). The distance between the two legs of this hairpin vortex is $70~\unicode[STIX]{x03BC}\text{m}$ at this cross-section. The wavelength of the undulations and the scale of surface corrugations presumably depend on other parameters such as surface tension and liquid viscosity; however, analysis of the effects of these parameters is beyond the scope of this article and would be an interesting subject for a prospective study.

Figure 30.

$\unicode[STIX]{x1D706}_{2}$ isosurface in a top close-up view of a liquid lobe in Domain III at low density ratio (Case D3b) at $64~\unicode[STIX]{x03BC}\text{s}$ (a); the solid black line shows the lobe front edge. The isosurfaces represent: the fast- and slow-moving KH vortices with $\unicode[STIX]{x1D706}_{2}=-7\times 10^{10}~\text{s}^{-2}$ (grey), the outer crest hairpin with $\unicode[STIX]{x1D706}_{2}=-4\times 10^{10}~\text{s}^{-2}$ (green) and the inner trough hairpin with $\unicode[STIX]{x1D706}_{2}=-3\times 10^{10}~\text{s}^{-2}$ (red). Lobe surface showing the corrugation formation from a top view at $64~\unicode[STIX]{x03BC}\text{s}$ (b), $66~\unicode[STIX]{x03BC}\text{s}$ (c) and $68~\unicode[STIX]{x03BC}\text{s}$ (d) of Case D3b.

The mechanism of corrugation formation at lower density ratios is similar to what was shown above in this section. However, the hairpins are not as easily bent (corrugated more slowly) as in high density ratios and, consequently, the number of corrugations are less. Figure 30 shows the vortex structures along with the images of a lobe from a top view in the $64$ – $68~\unicode[STIX]{x03BC}\text{s}$ period of Case D3b (see table 1). The vortex structures in figure 30(a) correspond to the liquid lobe in figure 30(b). The liquid surface has been removed from figure 30(a) to show all the vortices under the liquid surface; to facilitate the comparison between the vortex structures and surface structures, the lobe front edge is denoted in this figure by the solid black line. The outer hairpin (green isosurface) and the fast-moving KH vortex still lie above the lobe surface at $64~\unicode[STIX]{x03BC}\text{s}$ , while the inner hairpin (red isosurface) and the slow-moving KH vortex lie underneath the lobe and closely downstream of the lobe front. Following the illustrations in figure 25 (at $t_{3}$ and $t_{4}$ ), the tip of the outer hairpin stretches upstream while the inner hairpin stretches downstream under the influence of the fast-moving vortex. This process creates the hairpin undulations and increases the number of counter-rotating vortex pairs along the lobe span. However, comparison of figures 30(a) and 29 reveals that the hairpin corrugation appears much more quickly at high gas density, as three pairs of counter-rotating hairpins are seen at an earlier time in figure 29. At low gas density, the hairpins are not easily deflected and bent; the reason is conjectured to be the lower gas inertia. While the hairpins are in the gas, the lower inertia of the gas slows down the process of hairpin deflection. By the time the outer hairpin and the fast-moving vortex move into the liquid, the hairpin deformation hastens and another turn forms at its trough (at the centre of figure 30 a), where the green hairpin stretches downstream, over the fast-moving vortex. Thereafter, the number of undulations increases and the corrugations stretch. The relatively higher inertia of the neighbouring liquid retards the rate of growth of the instability. That is, the reaction of the hairpins to the KH vortex is slowed.

Figure 30 confirms that the lobe deformation is in fact an outcome of the KH and hairpin vortex distortions, not the other way around. The vortex and liquid structures at $64~\unicode[STIX]{x03BC}\text{s}$ , shown in figure 30(a,b), prove that the vortex deformation precedes the surface deformation. At this time, the hairpins have been already distorted and manifest two pairs of counter-rotating streamwise vortices near the lobe front, while the lobe front itself does not show any corrugation. At later times (figure 30 c,d), two corrugations start to form following the structure of the undulated hairpins. This transforms the lobe from a singular protrusion to two smaller ones growing out of the corrugated rim. While the lobe rim stretches under the influence of the KH vortex (after the two split counterparts merge), it also retracts towards the centre by the capillary forces. This brings the hairpins and the corrugations closer to the centre of the lobe, as seen in figure 30(d). This retraction is faster at lower $We_{g}$ ; compare figures 30 and 29.

3.4 Lobe and ligament stretching at low $Re_{l}$ ( $LoLiD$ mechanism)

At low $Re_{l}$ and low $We_{g}$ , the surface tension force resists perforation. The liquid viscosity is also fairly high and damps the small-scale corrugations on the lobe front edge. Consequently, the entire lobe stretches slowly into a thick and long ligament, which eventually breaks into large droplets. This terminates the $LoLiD$ breakup mechanism, which prevails in Domain I, shown in figure 1. The $LoLiD$ process is shown step by step in figure 31. The case shown in this figure (Case D1a) is taken for vortex analysis in this section. $\unicode[STIX]{x1D706}_{2}$ contours on $y$ -planes and $\unicode[STIX]{x1D706}_{2}$ isosurfaces in the interface vicinity are used in this section to explain the vortex dynamics of the $LoLiD$ process.

Figure 31.

Liquid surface deformation following the $LoLiD$ mechanism in Case D1a, at $t=26~\unicode[STIX]{x03BC}\text{s}$ (a), $36~\unicode[STIX]{x03BC}\text{s}$ (b), $40~\unicode[STIX]{x03BC}\text{s}$ (c), $44~\unicode[STIX]{x03BC}\text{s}$ (d), $46~\unicode[STIX]{x03BC}\text{s}$ (e), $48~\unicode[STIX]{x03BC}\text{s}$ (f) and $52~\unicode[STIX]{x03BC}\text{s}$ (g).

Figure 32.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest cross-section, at $t=6~\unicode[STIX]{x03BC}\text{s}$ (a), $10~\unicode[STIX]{x03BC}\text{s}$ (b) and $14~\unicode[STIX]{x03BC}\text{s}$ (c) of Case D1a.

As shown in figure 32(a), the process starts with a large KH roller downstream of the KH wave (denoted by the white arrow), hairpin filaments on the braid, and a much stronger hairpin filament on the KH wave crest (denoted by the dark grey arrow). Because of the high gas viscosity, the KH vortex diffuses much faster into the gas and at $6~\unicode[STIX]{x03BC}\text{s}$ the vortex core is entirely in the gas zone. Since the gas phase has higher velocity compared to the liquid phase, the gas distant from the liquid surface advects faster with respect to the interface. Hence, a few microseconds later, say at $10~\unicode[STIX]{x03BC}\text{s}$ (figure 32 b), the entire KH roller advects downstream with respect to the interface as its core moves farther from the interface via diffusion. The KH roller gains speed as it moves away from the interface. In contrast, the hairpin filaments that are closer to the liquid surface remain almost stationary with respect to the interface; see figure 32(c). The KH vortex gets larger with time as it diffuses; compare figure 32(a,c).

The KH roller reaches the neighbouring downstream crest hairpin at $16~\unicode[STIX]{x03BC}\text{s}$ (figure 33 a). Meanwhile, the braid hairpins overlap with the crest hairpins, constraining the lobe sheet in between. So far, the vortex dynamics manifest the conditions required for hole formation, i.e. overlapping of two oppositely oriented hairpins on top and bottom of the lobe, as well as for corrugation formation; i.e. constant pull in opposite directions induced by a moving vortex hovering over the hairpins. However, none of these structures are seen on the lobe at this moment (see figure 31 a). The inhibition of hole formation is due to the high surface tension. In such a low $We_{g}$ , the inertia and viscous forces are not strong enough to overcome the surface tension force to stretch and thin the lobe; hence, the lobe perforation is inhibited. Also, because of the high liquid viscosity at such a low $Re_{l}$ , the liquid surface deformation is much slower. Hence, the corrugation formation on the lobe edge does not occur as quickly as for the higher $Re_{l}$ . Small-scale corrugations are damped by both the high viscous and high surface tension forces. Thus, the lobe slowly stretches into a thick ligament. Moreover, the KH roller is also farther away from the interface in this case compared to the $LoCLiD$ process (compare figures 33 a and 24 c), which means that it has a much weaker influence on the hairpin filaments.

Figure 33.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest cross-section, at $t=16~\unicode[STIX]{x03BC}\text{s}$ (a), $20~\unicode[STIX]{x03BC}\text{s}$ (b) and $22~\unicode[STIX]{x03BC}\text{s}$ (c) of Case D1a.

Figure 34.

Schematic of the recirculation zone downstream of the KH wave. The hollow circles denote the KH vortex and the black circles denote the newly formed vortex in the recirculation zone. The induction of each vortex on the neighbouring vortices is indicated by the straight arrows with the same shading as the source vortex. The curly arrows qualitatively indicate the streamlines near the interface.

The tip of the KH wave acts like a backward-facing step for the gas flow and a temporary recirculation zone forms immediately downstream of the wave, as shown in figure 34. This is very similar to the scenario seen at low density ratios by Hoepffner et al. (Reference Hoepffner, Blumenthal and Zaleski2011). They also showed that the head of the wave appears to the gas stream as a fixed obstacle, with the ensuing vortex shedding. The high strain rate at the braid depletes the vorticity upstream of the braid and collects the vorticity into a new vortex just downstream of the KH wave, in the recirculation zone at $20~\unicode[STIX]{x03BC}\text{s}$ . This new vortex is indicated by the thick black arrow in figure 33(b) and the black circles in figure 34. When the KH roller (vortex 1) passes over the KH wave front (see $t_{1}$ in figure 34), the mutual induction of this roller and the new vortex (vortex 3) causes the KH roller to advect farther downstream, further stretching the KH vortex. The local induction of each vortex on the neighbouring vortices is indicated by the straight arrows of the same shading as the source vortex in figure 34. When the KH vortex (vortex 1) passes over the downstream wave at $22~\unicode[STIX]{x03BC}\text{s}$ (figure 33 c), the mutual induction of the vortices pushes away the KH vortex, advecting it in the $z$ -direction; see $t_{2}$ in figure 34. The stronger induction by vortex 4 pushes vortex 1 in the direction normal to the line connecting the centres of the two vortices (in the direction of the black arrow). Vortex 3 also induces a flow that would bring vortex 1 closer to the interface; however, this induction is much weaker as vortex 3 is more distant from vortex 1 compared to vortex 4. The new vortex wraps under and around the KH roller at its crest and later over the trough of the adjacent downstream KH roller, as shown in figure 35. Thus, this new vortex (the red tube in figure 35) transforms into a hairpin vortex that is stretched between two adjacent KH vortices and wraps around them at their crests and troughs.

Figure 35.

3D Schematics showing the vortex structures in the $LoLiD$ Domain I (a) – $A$ is the plane in which (b) is drawn; cross-sectional view of the $A$ -plane, showing the spanwise squeezing of the lobe by induced flow of the hairpin vortices (b). The vortex schematics are periodic in $x$ - and $y$ -directions.

Figure 36.

$\unicode[STIX]{x1D706}_{2}$ contours on the spanwise crest cross-section at $t=26~\unicode[STIX]{x03BC}\text{s}$ (a) and $30~\unicode[STIX]{x03BC}\text{s}$ (b) of Case D1a.

Schematics of the vortex structures corresponding to $t=26~\unicode[STIX]{x03BC}\text{s}$ (shown in figure 36 a) are illustrated in figure 35. The KH vortex has a larger undulation in this domain compared to the other two domains and is also farther away from the interface in the gas zone (compare figure 35 with figures 11 and 27). Two pairs of counter-rotating hairpins – one on the lobe crest (the black tube), and the other on the trough (the red tube) – stretch and wrap around the KH vortex (figure 35 a). Notice that the red hairpin in figure 35 is actually the new vortex that was formed at the wave front end and was indicated by the black arrows in figures 33 and 36. These hairpins are periodic in both spanwise ( $y$ ) and streamwise ( $x$ ) directions; i.e. the tubes that emerge at the bottom corner or the left side of figure 35(a) re-enter from the top corner or the right side of the sketch, respectively. As shown in the cross-sectional view of the $A$ -plane passing through the lobe in figure 35(b), both the black and the red hairpins are located slightly above and on both sides of the lobe at this moment. While the flow induced by the KH vortex creates a streamwise flow on top and bottom of the lobe, the gas flow induced by these two counter-rotating hairpins (shown by the curly arrows is figure 35 b) generates a spanwise flow towards the lobe midplane. Consequently, the lobe is both squeezed in the spanwise direction – via the induced flow of the hairpins – and stretched in the streamwise direction by the induced flow of the KH vortex. The gas flow induced by the hairpins also lifts the lobe in the $z$ -direction; see figure 35(b).

Figure 37.

Velocity fluctuation vector field with respect to the average velocity of the vortex cores near the interface, and $\unicode[STIX]{x1D714}_{y}$ contours on the spanwise crest cross-section at $t=36~\unicode[STIX]{x03BC}\text{s}$ of Case D1a. The thick arrows refer to the vortices shown in figure 36(b). The star symbol marks the saddle point, and the liquid is shown in light blue.

While the KH vortex keeps diffusing due to the high gas viscosity, the dark grey and black vortices wrap around it and are stretched with it as it moves downstream (figure 36 b). Meanwhile, another vortex forms downstream of the wave crest at $30~\unicode[STIX]{x03BC}\text{s}$ , indicated by the light grey arrow in figure 36(b). This vortex advects downstream closely hovering over the interface at $36~\unicode[STIX]{x03BC}\text{s}$ , while new vortices keep forming downstream of the KH wave crest at the recirculation zone to replace it; see figure 37. Hoepffner et al. (Reference Hoepffner, Blumenthal and Zaleski2011) were the first to identify this vortical mechanism at low density ratios. Similar to figure 37, they showed that two vortices exist in the gas – one is at the instant of leaving the shelter of the wave (the dark grey curly arrow in figure 37), and the second further downstream is the result of the previous shedding event (the light grey curly arrow in figure 37). The induced motion of these vortices entrains the gas and draws it under the protruding lobe (follow the velocity vectors in the gas zone near the interface in figure 37). As is well known, the streamwise strain rate is highest at the saddle (Hussain Reference Hussain1986; Lasheras & Choi Reference Lasheras and Choi1988). The location of the saddle point is marked by a star symbol in figure 37. The vortices that reach the saddle point stretch in the streamwise direction and transform into hairpin vortices with streamwise legs. The location of the black and dark grey hairpins that are wrapping around the KH vortex are denoted by straight arrows of the corresponding vortex colour. Since these vortices are mostly streamwise at this cross-section, they do not have $\unicode[STIX]{x1D714}_{y}$ contours. The streamwise elongation of the vortices progressively aligns the two counter-rotating legs of the hairpins in the streamwise direction. The self-induction of the counter-rotating legs of the KH vortex and the hairpin vortices, moves the vortex tubes in the normal direction away from the interface, as shown schematically in figure 38.

Figure 38.

Schematic showing the self-induction of the streamwise counter-rotating hairpin legs (a), and in a $yz$ -plane crossing the legs (b). The grey arrows in (b) show the direction of the self-induction.

Figure 39.

$\unicode[STIX]{x1D706}_{2}$ isosurface in a top close-up view of a liquid lobe in Domain I (Case D1a) (left) and the liquid surface from a top view (right), at $t=6~\unicode[STIX]{x03BC}\text{s}$ (a), $16~\unicode[STIX]{x03BC}\text{s}$ (b) and $26~\unicode[STIX]{x03BC}\text{s}$ (c). The solid black line in the $\unicode[STIX]{x1D706}_{2}$ isosurface images shows the location of the lobe front edge. The isosurfaces represent: the KH vortex with $\unicode[STIX]{x1D706}_{2}=O(-10^{10})$ (grey), the crest hairpin with $\unicode[STIX]{x1D706}_{2}=O(-10^{11})$ (green) and the trough hairpin with $\unicode[STIX]{x1D706}_{2}=O(-10^{11})$ (red).

The temporal evolution of the vortex structures in the $LoLiD$ process is illustrated in figure 39 from a top view of the liquid surface. The liquid surface is illustrated on the right panel of this figure, and the $\unicode[STIX]{x1D706}_{2}$ isosurfaces (vortex filaments) are depicted on the left. In the beginning of the process (figure 39 a), the KH vortex (grey tube) and the crest hairpin filament (green tube) are completely spanwise oriented (in the $y$ -direction). The braid hairpins (thinner green tubes) are also primarily spanwise with some undulations and are half a wavelength ( $180^{\circ }$ ) out of phase with respect to the crest hairpins in $y$ . This is consistent with the findings of Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) for non-homogeneous round jets, and Danaila et al. (Reference Danaila, Dušek and Anselmet1997) for like-density jets. At $16~\unicode[STIX]{x03BC}\text{s}$ , the KH vortices stretch more; the hairpins on the braid also collect into the crest hairpin and manifest higher undulations (green isosurface in figure 39 b). The trough vortex (red isosurface), which was first seen in figure 33(b) (indicated by black arrow) also manifests at this time. Both the green and the red vortex filaments slide underneath the KH roller at the spanwise wave crest (lobe crest), and wrap over it at the spanwise wave trough (lobe sides). The wave front has been marked with a solid black line in figure 39 to help find the location of the lobe crest and trough with respect to the vortices. The KH roller in this process is more stretched than the other two mechanisms (compare figure 39 with figures 17 and 26).

The vortices become streamwise near the braid in both spanwise crest and trough at $26~\unicode[STIX]{x03BC}\text{s}$ (figure 39 c). Except for the tip of the KH roller, which is still spanwise, the rest of the vortex structures are nearly streamwise-oriented. Both the green and the red hairpins wrap around the KH vortex on top of each other (see figures 35 a and 39 c). As the vortices get stretched in the streamwise direction, the lobes follow their shape and get thinner in the spanwise direction because of the induced gas flow by the hairpin legs; compare the shape of the lobes with that of the KH vortex just downstream of the lobe in figure 39(c) to observe their similarity. A neck starts to form on the lobe edge right between the two hairpins on the two sides of the lobe. This completely follows the mechanism introduced in figure 35.

Figure 40.

Temporal evolution of the vortices indicated by the $\unicode[STIX]{x1D706}_{2}$ isosurface coloured by the streamwise velocity contours, at $t=6~\unicode[STIX]{x03BC}\text{s}$ (a), $16~\unicode[STIX]{x03BC}\text{s}$ (b), $26~\unicode[STIX]{x03BC}\text{s}$ (c), $40~\unicode[STIX]{x03BC}\text{s}$ (d), $48~\unicode[STIX]{x03BC}\text{s}$ (e) and $52~\unicode[STIX]{x03BC}\text{s}$ (f) of Case D1a. The liquid surface is shown in blue. The arrows refer to the vortices denoted in figures 32–36.

Figure 40 shows the cascade and deformation of the vortex filaments in the period $6~\unicode[STIX]{x03BC}\text{s}$ – $52~\unicode[STIX]{x03BC}\text{s}$ . The $\unicode[STIX]{x1D706}_{2}$ isosurfaces are coloured by the streamwise velocity contours here. The vortices start from a spanwise orientation and gradually turn streamwise. The vortex structures cascade from thick and uniform vortices to thin and chaotic structures as they stretch. The gradual departure of the vortices away from the liquid surface is also shown in this figure, where the KH roller, indicated by the white arrow, starts on the liquid surface at $6~\unicode[STIX]{x03BC}\text{s}$ (figure 40 a) and advects away from the interface as it gets streamwise due to the self-induction process described in figure 38.

At $40~\unicode[STIX]{x03BC}\text{s}$ (figure 40 d), almost all of the vortex structures have become streamwise except for a short section at the tip of the KH roller. This confirms that the vortices become streamwise-oriented as they get close to the saddle point. Meanwhile, the vortices cascade into smaller structures due to turbulence; compare figures 40(c) and 40(d). As the legs of the hairpins get closer to each other, the mutual induction of the legs lifts the hairpins in the normal direction (figure 38), and locates the hairpins above the lobes. The lobes get thinner as the streamwise counter-rotating legs induce a gas flow in the spanwise direction towards the lobe midplane, squeezing the lobe and transforming it into a ligament. Ligament creation is strongly correlated with its local velocity field, and is induced by the local shear, as indicated by Shinjo & Umemura (Reference Shinjo and Umemura2010). Nearby vortices determine the ligament formation direction. Spanwise vortices form ligaments normal to the injection direction (i.e. spanwise) and streamwise vortices form ligaments parallel to the injection direction (i.e. streamwise). The influence of the vorticity field on the ligament orientation is consistent with the findings of Shinjo & Umemura (Reference Shinjo and Umemura2010). Following the streamwise-oriented hairpin vortices in the $LoLiD$ process, mostly streamwise ligaments form in Domain I.

Figure 41.

Pressure contours on the surface of a ligament at $t=42~\unicode[STIX]{x03BC}\text{s}$ (a), $44~\unicode[STIX]{x03BC}\text{s}$ (b), $45~\unicode[STIX]{x03BC}\text{s}$ (c) and $46~\unicode[STIX]{x03BC}\text{s}$ (d) of Case D1a, showing the ligament pinch-off following the short-wave breakup mode.

Velocity contours in figure 40(e) show that the streamwise velocity of the vortices increases at higher $z$ -levels. This can be clearly seen from the gradient of colours from green ( $50~\text{m}~\text{s}^{-1}$ ) – near the liquid surface – to red ( $90~\text{m}~\text{s}^{-1}$ ) – at the tip of the KH vortex. The liquid surface also experiences a similar velocity gradient, where the velocity at the tip of the lobe is higher compared to its root; i.e. where the lobe connects to the liquid sheet. This velocity gradient manifests how the lobe elongates under the streamwise strain and forms the ligament. The ligaments finally pinch off and create droplets, as shown in figure 40(f). The ligament pinch-off follows the short-wave breakup mode introduced by Shinjo & Umemura (Reference Shinjo and Umemura2010). The ligament acts as a very small round liquid jet emanating from the sheet surface. The ligament tip pressure is high in the beginning and contracts due to surface tension and pushes the inner part along the ligament axis (figure 41 a). This motion emanates compression waves in the upstream direction and a neck forms. The tip bulb grows as it absorbs the liquid from the upstream liquid sheet by contraction (figure 41 b). As the tip bulb size grows, its inner pressure drops. Consequently, the tip bulb sucks the liquid further from the neck, as shown in figure 41(b) – causing the neck to become narrower. When the neck becomes thin enough and its pressure high enough, the circumferential surface tension cuts the neck and a droplet pinches off (figure 41 c). This process reiterates from the beginning, and the next droplets pinch off identically. The droplets fly away from the interface under the vortex induction and gain a higher velocity than the liquid sheet; therefore, the droplets advect downstream with respect to the jet (figure 40 f).

Figure 42.

$\unicode[STIX]{x1D706}_{2}$ isosurface in a top close-up view of a liquid lobe in Domain I at low density ratio (Case D1b) at $16~\unicode[STIX]{x03BC}\text{s}$ (a); the solid black line shows the lobe front edge. The isosurfaces represent: the KH vortex with $\unicode[STIX]{x1D706}_{2}=-10^{11}~\text{s}^{-2}$ (grey), the outer crest hairpin with $\unicode[STIX]{x1D706}_{2}=-2\times 10^{11}~\text{s}^{-2}$ (green) and the inner trough hairpin with $\unicode[STIX]{x1D706}_{2}=-3\times 10^{11}~\text{s}^{-2}$ (red). Lobe surface showing the corrugation formation from a top view at $16~\unicode[STIX]{x03BC}\text{s}$ (b), $18~\unicode[STIX]{x03BC}\text{s}$ (c) and a 3D view at $22~\unicode[STIX]{x03BC}\text{s}$ (d) of Case D1b.

The vortex structures and lobe deformation at low density ratio (Case D1b) are shown in figure 42. The lattice made by the KH vortex and the hairpin vortices is very similar to what was proposed in figure 35 for high density ratios. The only difference being that as gas density decreases, the KH vortices are less bound to the liquid and depart from the interface much more easily and meanwhile stretch in both streamwise and normal $z$ -directions (figure 42 a). All the vortices in figure 42(a) exist in the gas zone and above the liquid interface. Following this vortex deformation, the lobe is squeezed and thinned in the spanwise direction, while being stretched and lifted by the KH vortex normally outward. This causes the ligaments to protrude more and quicker in the normal direction (figure 42 d) compared to higher $\hat{\unicode[STIX]{x1D70C}}$ . This figure also proves that the vortex deformation leads the surface deformation, and not the other way around.

As will be discussed in § 3.5, when density ratio gets higher, i.e. higher $We_{g}$ , the vortices are closer to the interface and are less stretched, and more spanwise- rather than streamwise-oriented. Hence, the gas flow that is induced by the KH vortex shears the lobe from its top and bottom sides. In this scenario, the lobe thins in the normal direction instead of the spanwise direction, as shown schematically in figure 43(a). The blue curly arrows illustrate the qualitative streamlines in this situation, and the red straight arrows denote the direction of lobe squeeze. The shear due to the gas flow induced by the spanwise vortices stretches the lobe in the streamwise direction and thins the lobe and makes it more vulnerable to puncture. This is consistent with the map in figure 1; as $We_{g}$ increases, the breakup mechanism shifts from $LoLiD$ to $LoHBrLiD$ .

Figure 43.

Schematic showing the lobe thinning in the normal direction in the $LoHBrLiD$ process from a $xz$ -plane (a), and lobe thinning in the spanwise direction in the $LoLiD$ process from a $xy$ -plane (b).

At lower density ratio and lower $We_{g}$ (Domain I), the vortices are more stretched in the streamwise direction. In this domain ( $LoLiD$ mechanism), the streamwise legs of the hairpin vortex on the two sides of the lobe induce a gas flow in the spanwise direction towards the lobe midplane, which squeezes and thins the lobe in the spanwise direction and transforms it into a thick ligament, as shown schematically from the top view of a lobe in figure 43(b). A neck forms at the location of the red arrows. The induced motion of the hairpin legs also lifts the lobe in the $z$ -direction. Eventually, the ligament pinches off at the neck and a droplet forms following the short-wave mode discussed above.

In summary, whether the lobe thins in the $z$ -direction and perforates or thins in the $y$ -direction and forms a ligament depends on the orientation of the vortices in the vicinity of the lobe. The spanwise vortices result in hole formation and spanwise bridges, while streamwise vortices result in spanwise lobe compression and streamwise ligament stretching.

3.5 Streamwise vorticity generation

The streamwise vorticity ( $\unicode[STIX]{x1D714}_{x}$ ) is crucial in initiation of the 3D instability on liquid jets. $\unicode[STIX]{x1D714}_{x}$ generation via vortex stretching and vortex tilting – i.e. strain–vorticity interactions – and baroclinic effects are studied in this section for a low density ratio of $0.05$ and a high density ratio of $0.5$ . The contributions of the different terms in the vorticity equation to $\unicode[STIX]{x1D714}_{x}$ generation are compared at two distinctly different density ratios, to understand the role of density ratio in the liquid-jet breakup. Specifically, the generation of the three-dimensionality on the liquid sheet interface is addressed in this section. Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) performed a similar analysis for the round jet at several different density ratios. They showed for a density ratio of 0.01 that the baroclinic effect, i.e. the Rayleigh–Taylor (RT) instability, is the dominant cause of the initiation of 3D structures. This is consistent with the suggestion of Marmottant & Villermaux (Reference Marmottant and Villermaux2004), who performed their experiment at lower pressures and gas density. However, for a density ratio of 0.1 or greater, Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) showed that the azimuthal tilting and radial tilting of the ring vortices are the dominant effects in streamwise vorticity generation.

The complete vorticity equation is

(3.1) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\text{D}\unicode[STIX]{x1D74E}}{\text{D}t}=(\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-\unicode[STIX]{x1D74E}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})+\unicode[STIX]{x1D735}\times \left(\frac{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}}{\unicode[STIX]{x1D70C}}\right)+\frac{1}{\unicode[STIX]{x1D70C}^{2}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\times \boldsymbol{F}_{\unicode[STIX]{x1D748}}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ and $\unicode[STIX]{x1D74E}$ are the velocity and vorticity vectors, respectively. $\unicode[STIX]{x1D749}$ is the viscous stress tensor, and $\boldsymbol{F}_{\unicode[STIX]{x1D748}}$ is the surface tension force. Since the fluids are incompressible in this study, the second term on the right-hand side is zero. A simple dimensional analysis shows that the viscous diffusion term (the third term on the right-hand side) scales as $\unicode[STIX]{x1D707}U/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6E5}^{3}$ , which for a typical case considered in our study, e.g.  $\unicode[STIX]{x1D707}=O(10^{-3})~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$ , $U=O(10)~\text{m}~\text{s}^{-1}$ , $\unicode[STIX]{x1D70C}=O(10^{3})~\text{kg}~\text{m}^{-3}$ , and a mesh size of $\unicode[STIX]{x1D6E5}=2.5~\unicode[STIX]{x03BC}\text{m}$ , gives a magnitude of ${\approx}O(10^{10})\text{s}^{-2}$ . The surface tension term (the last term) scales as $\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6E5}^{2}$ ; which for a typical case with $\unicode[STIX]{x1D70E}=O(10^{-2})~\text{N}~\text{m}^{-1}$ , and a radius of curvature of $100~\unicode[STIX]{x03BC}\text{m}$ ( $\unicode[STIX]{x1D705}=10^{4}~\text{m}^{-1}$ ), also gives a result on the order of $10^{10}~\text{s}^{-2}$ . As will be shown in this section, these two terms are two to three orders of magnitude smaller than the other terms in (3.1), and therefore have negligible contributions to vorticity generation; thus, the rate of change of $\unicode[STIX]{x1D714}_{x}$ is approximately

(3.2) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\text{D}\unicode[STIX]{x1D714}_{x}}{\text{D}t}=\unicode[STIX]{x1D714}_{x}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D714}_{y}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D714}_{z}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}+\frac{1}{\unicode[STIX]{x1D70C}^{2}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{x}$ , $\unicode[STIX]{x1D714}_{y}$ , $\unicode[STIX]{x1D714}_{z}$ and $u$ denote the streamwise, spanwise and cross-stream (normal) vorticities, and the streamwise velocity, respectively. The terms on the right-hand side denote streamwise stretching, spanwise tilting, normal tilting, baroclinic effect due to normal pressure gradient, and baroclinic effect due to spanwise pressure gradient, respectively. Density gradient is normal to the liquid interface; i.e. approximately in the $z$ direction. The spanwise density gradient, i.e.  $\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}y$ , is negligible compared to $\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}z$ , since the sheet cross-section remains fairly rectangular during early instability development. However, this term accounts for the baroclinic effect, which can deform the interface in the spanwise direction.

In the cases studied in this section, no initial perturbation is imposed on the liquid surface. Except for numerical errors which for our purpose correspond to small random physical disturbances, all terms in the $\unicode[STIX]{x1D714}_{x}$ generation (3.2) are initially zero. Namely, the first and the third terms are zero because there are no vorticity components in the $x$ and $z$ -directions initially, and the second term is zero since the streamwise velocity is uniform in the spanwise direction. The baroclinic terms are identically zero since density and pressure gradients in the $y$ -direction are initially zero.

Since the streamwise ( $\unicode[STIX]{x1D714}_{x}$ ) and normal vorticities ( $\unicode[STIX]{x1D714}_{z}$ ) cannot be generated, but can be enhanced, the main source of $\unicode[STIX]{x1D714}_{x}$ generation at early times is either the spanwise vorticity tilting or the baroclinic torque, which become non-zero as a result of small perturbations of $u$ and $p$ in the spanwise direction. This intuition is consistent with the results of Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) in round jets at early times, where for a wide range of density ratios, the baroclinic torque and the azimuthal vortex tilting terms are dominant for the first $5~\unicode[STIX]{x03BC}\text{s}$ of their computations; however, they might be overtaken later by other terms. Baroclinicity becomes more pronounced at lower density ratios, since the density gradient across the interface is higher.

In the data analysis, the gradients have been calculated and averaged over the computational interface thickness that equals three mesh points in the $z$ and tangential directions. The terms in (3.2) are:

1. (i) streamwise vortex stretching: $\unicode[STIX]{x1D714}_{x}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x)$ ;

2. (ii) spanwise vortex tilting: $\unicode[STIX]{x1D714}_{y}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y)$ ;

3. (iii) normal vortex tilting: $\unicode[STIX]{x1D714}_{z}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z)$ ;

4. (iv) baroclinic vorticity generation: $1/\unicode[STIX]{x1D70C}^{2}\left[\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}y~(\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}z)-(\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}z)~\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}y\right]$ .

Since the peak of the streamwise vorticity occurs at the wave braids (as shown by figure 49 a), the absolute value of the four above-mentioned terms are averaged only at the braid region. Both top and bottom surfaces have been considered in this measurement.

Two different density ratios have been analysed in this section. The non-dimensional characteristics of these two cases are: $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ , and $\hat{\unicode[STIX]{x1D70C}}=0.05$ and $0.5$ . The sheet thickness is $h=50~\unicode[STIX]{x03BC}\text{m}$ , and no initial perturbation is imposed on the liquid–gas interface.

Figure 44.

Contributions of streamwise vortex stretching (squares), spanwise vortex tilting (triangles), normal vortex tilting (diamonds), and baroclinicity (circles) to the generation of $\unicode[STIX]{x1D714}_{x}$ at the liquid surface for two density ratios; $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ ; $\hat{\unicode[STIX]{x1D70C}}=0.05$ (a) and $\hat{\unicode[STIX]{x1D70C}}=0.5$ (b).

Figure 44 shows the contribution of each term in the generation of $\unicode[STIX]{x1D714}_{x}$ for the high and low density ratios, in the first $20~\unicode[STIX]{x03BC}\text{s}$ . Baroclinicity (circles) is the most important factor at low density ratio, since the density gradient normal to the interface is higher, and also the local density in the gas zone near the interface is lower (see the baroclinic vorticity generation expression given above). However, at high density ratio, baroclinicity is the least significant. Baroclinicity is only slightly larger than the streamwise stretching (squares) in the beginning of the computations for high density ratio, but it is outrun by this term at approximately $13~\unicode[STIX]{x03BC}\text{s}$ and remains the lowest of all terms thence. The baroclinic vorticity generation term is an order of magnitude smaller than the vortex stretching and tilting terms at the end of the computations.

In both low and high density ratios, the spanwise and normal tilting seem to be more important than the streamwise stretching for the first 20 ms. The spanwise tilting is high because the initially spanwise vortex lines are gradually tilted in the streamwise direction. The normal tilting is also high since the velocity gradient is much higher in the normal direction ( $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z$ ) in the beginning of the computations. The streamwise stretching, however, has the highest growth rate and almost reaches the magnitude of vortex tilting at approximately $20~\unicode[STIX]{x03BC}\text{s}$ . Later on, as $\unicode[STIX]{x1D714}_{x}$ grows, the vortex stretching becomes more significant. These results are consistent with the findings of Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) for round liquid jets. Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) also concluded that, generally, the importance of baroclinic vorticity generation (RT instability) has been overemphasized in the literature, especially at very high pressures, and other important aspects of vorticity dynamics and similarities with injection into an alike fluid have been neglected. As described by them and also evident in our figures, the vortex tilting terms are the largest at early times; however, new findings show some cancellations, discussed below.

Figure 45.

Contours of streamwise vortex stretching (a), spanwise vortex tilting (b), normal vortex tilting (c) and baroclinic generation (d) on a $y$ -plane at $t=13~\unicode[STIX]{x03BC}\text{s}$ for $\hat{\unicode[STIX]{x1D70C}}=0.05$ ; $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ .

Figure 46.

Contours of streamwise vortex stretching (a), spanwise vortex tilting (b), normal vortex tilting (c) and baroclinic generation (d) on a $x$ -plane at $t=13~\unicode[STIX]{x03BC}\text{s}$ for $\hat{\unicode[STIX]{x1D70C}}=0.05$ ; $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ .

Figures 45 and 46 show the contours of the four $\unicode[STIX]{x1D714}_{x}$ -generating terms at $13~\unicode[STIX]{x03BC}\text{s}$ , on a $y$ -plane and $x$ -plane, respectively. The spanwise and normal vortex tilting terms (figures 45 b,c and 46 b,c) are stronger than the streamwise stretching, but with opposite signs. As demonstrated in figure 44, the two vortex tilting terms are also nearly equal. A closer look at their equations explains the reason:

(3.3) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D714}_{y}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y} & = & \displaystyle \frac{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{u}}{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{z}}\frac{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{u}}{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{y}}-\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y},\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D714}_{z}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z} & = & \displaystyle \frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}-\frac{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{u}}{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{y}}\frac{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{u}}{\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{z}}.\end{eqnarray}$$

The two terms in bold font are exactly the same but have opposite signs. Thus, the only difference in the magnitudes of the spanwise and normal tilting comes from the other terms (the second term in 3.3 and the first term in 3.4). However, since these terms deal with gradients of $v$ and $w$ of $O(1)$ , which are two orders of magnitude smaller than $u$ of $O(10^{2})$ , they are much smaller than the boldface terms, at the beginning of the computations. Thus, the deviation of the vortex tilting terms from the boldface terms is very small, and the two terms nearly cancel each other early on. Based on this, our earlier conclusion (and also that of Jarrahbashi & Sirignano Reference Jarrahbashi and Sirignano2014) should be modified: the spanwise and normal vortex tilting terms, even though the largest among the $\unicode[STIX]{x1D714}_{x}$ -generating terms, are not the most important in $\unicode[STIX]{x1D714}_{x}$ generation, since they nearly cancel each other. The streamwise vortex stretching and baroclinic effects (RT instability) are the most important in generation of $\unicode[STIX]{x1D714}_{x}$ , at high and low density ratios, respectively.

Figure 47.

Perturbation velocity vectors superimposed on $\unicode[STIX]{x1D714}_{y}$ contours at $t=13~\unicode[STIX]{x03BC}\text{s}$ on a $y$ -plane;   $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D70C}}=0.05$ and $\hat{\unicode[STIX]{x1D707}}=0.0066$ . Green arrows denote the maximum stretch along the diverging separatrix. The liquid–gas interface is indicated by the red line.

Figure 45(a) also confirms that the vortex stretching originates from the braids first, as the strain due to the adjacent primary vortical structures is highest at the saddle (braid) and the ribs are aligned along the diverging separatrix (Hussain Reference Hussain1986). Note that most of the coloured spots in figure 45(a) are on the braids and not the wave crests. This is consistent with the experimental observations of earlier researchers (Bernal & Roshko Reference Bernal and Roshko1986; Lasheras & Choi Reference Lasheras and Choi1988; Liepmann & Gharib Reference Liepmann and Gharib1992) for uniform-density flows. The location and direction of the stretch can also be seen in figure 47, which shows the fluctuation velocity vectors relative to the average KH vortex velocity on a blown-up section of the liquid jet at $13~\unicode[STIX]{x03BC}\text{s}$ . It is evident that the saddles with the highest strain rate are on the braids between two adjacent vortices, where the fluctuation velocity vectors depart in the opposite directions. The stretch direction at the saddle point is shown by the green arrows in figure 47. The saddle points are in the gas phase, close to the interface. The streamwise vortex stretching is highest at the saddle points (see figure 45 a), where the flow is primarily discrete ribs (Hussain Reference Hussain1986). The fluid elements are stretched along the interface, i.e. along the diverging separatrix shown by the green arrows, and compressed normal to the interface at the saddle points. The centre of the spanwise vortices (rolls) are at the crest of the interface waves, as denoted by the velocity vectors and the vorticity contours. Interestingly, the vorticity peak coincides with the interface at the crests, but not at the troughs, where the vorticity has migrated into the gas phase further away from the interface. Note that the vorticity contours look like crescents and not circular as the vorticity is not uniformly distributed as in a vortex column (like the typical Oseen vortex; see Pradeep & Hussain Reference Pradeep and Hussain2006); the non-centric vorticity distribution is due to the curved KH rollers.

Figure 48.

Contours of baroclinicity on a $y$ -plane at $t=13~\unicode[STIX]{x03BC}\text{s}$ for $\hat{\unicode[STIX]{x1D70C}}=0.5$ ; $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ .

The stretching and tilting terms are centred at the interface (see figures 45 a–c and 46 a–c), but the baroclinic torque term $1/\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \unicode[STIX]{x1D735}p$ is always larger in the gas phase (see figures 45 $d$ and 46 $d$ ). Baroclinic generation, being proportional to $1/\unicode[STIX]{x1D70C}^{2}$ , is two orders of magnitude larger in the gas phase compared to the liquid phase (for $\hat{\unicode[STIX]{x1D70C}}=0.05$ ). As density ratio increases, the difference between the gas and liquid densities decreases; thus, the contours of baroclinicity get closer to the liquid interface. This is seen in the baroclinicity contours of $\hat{\unicode[STIX]{x1D70C}}=0.5$ in figure 48. Since the density ratio is an order of magnitude higher than that of figure 45(d), the local density in the gas zone is much higher, hence the baroclinicity in the gas is lower and closer to its value in the liquid. Thus, the contours are closer to the interface (compare figures 48 and 45 $d$ ). This contributes to the $\unicode[STIX]{x1D714}_{x}$ peak being closer to the interface and growing larger compared to the low density ratios, hence creating and stretching more lobes at higher density ratios. Also, the peak of baroclinic torque is an order of magnitude smaller in figure 48 compared to figure 45(d), since the density gradient normal to the interface is lower at higher density ratios.

Figure 49.

Contours of streamwise vorticity $\unicode[STIX]{x1D714}_{x}$ (a), and spanwise vorticity $\unicode[STIX]{x1D714}_{y}$ (b) at $t=13~\unicode[STIX]{x03BC}\text{s}$ for $\hat{\unicode[STIX]{x1D70C}}=0.05$ ; $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ .

As shown in figure 45(d), the baroclinicity contours change sign (i.e. colour) continuously in $x$ . This pattern is very similar to the hairpin pattern seen in the $\unicode[STIX]{x1D714}_{x}$ contours presented by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014). The $\unicode[STIX]{x1D714}_{x}$ and $\unicode[STIX]{x1D714}_{y}$ contours for the current case are illustrated in figure 49. Comparison of the contours of $\unicode[STIX]{x1D714}_{x}$ and baroclinicity (compare figures 49 a and 45 d) shows that they follow a very similar pattern. Hence, we can conclude that baroclinicity is the most important factor in creation of the hairpin vortex structure at low density ratios. The role of the baroclinic torque in deformation of the surface waves in a stratified shear layer is studied in detail by Schowalter et al. (Reference Schowalter, Van Atta and Lasheras1994), and is consistent with our results. The vector field of figure 47 shows that $\unicode[STIX]{x1D714}_{x}$ and $\unicode[STIX]{x1D714}_{y}$ are highest at the braids and wave crests, respectively. This is also evident in the vorticity contours of figure 49.

In order to understand how fast the liquid sheet deforms and manifests 3D instabilities, the magnitudes of the vorticity components are examined through time. The 3D instabilities are directly related to the magnitudes of $\unicode[STIX]{x1D714}_{x}$ and $\unicode[STIX]{x1D714}_{z}$ against $\unicode[STIX]{x1D714}_{y}$ , which exists from the beginning when the flow is still 2D. As mentioned earlier by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014), Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016), and Zandian et al. (Reference Zandian, Sirignano and Hussain2016), the streamwise vorticity is the main cause of the 3D instabilities and interface distortion.

Figure 50.

Temporal variation of the vorticity components for $Re_{l}=2500$ , $We_{l}=14\,400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ ; $\hat{\unicode[STIX]{x1D70C}}=0.05$ (a) and $\hat{\unicode[STIX]{x1D70C}}=0.5$ (b).

The absolute value of each vorticity component, averaged over the entire liquid–gas interface, is plotted in figure 50 for low and high density ratios. The spanwise vorticity $\unicode[STIX]{x1D714}_{y}$ (KH vortex) is the only component that exists initially. In both cases, $\unicode[STIX]{x1D714}_{y}$ grows for a certain time, until approximately $15~\unicode[STIX]{x03BC}\text{s}$ , and then onward keeps more or less the same order of magnitude. The spanwise vorticity is larger for lower density ratio.

$\unicode[STIX]{x1D714}_{x}$ and $\unicode[STIX]{x1D714}_{z}$ grow very slowly up to $10~\unicode[STIX]{x03BC}\text{s}$ , after which there is a sudden increase in their growth rate. In both cases, $\unicode[STIX]{x1D714}_{x}$ and $\unicode[STIX]{x1D714}_{z}$ are of the same order of magnitude, which indicates that initially spanwise vortex filaments are lifted in the normal direction and tilted in the streamwise direction at almost equal rates. The growth rate of $\unicode[STIX]{x1D714}_{x}$ is higher at higher density ratios. $\unicode[STIX]{x1D714}_{x}$ reaches the same order of magnitude as $\unicode[STIX]{x1D714}_{y}$ at $t=20~\unicode[STIX]{x03BC}\text{s}$ , for the high density ratio. For the lower density ratio, however, the growth is slower. Hence, three-dimensionality manifests sooner at higher gas densities. For $\hat{\unicode[STIX]{x1D70C}}=0.05$ , the magnitude of $\unicode[STIX]{x1D714}_{x}$ and $\unicode[STIX]{x1D714}_{z}$ are still half of $\unicode[STIX]{x1D714}_{y}$ at $20~\unicode[STIX]{x03BC}\text{s}$ ; notice that $\unicode[STIX]{x1D714}_{y}$ always remains much higher. Thus, $\unicode[STIX]{x1D714}_{y}$ is still dominant and 2D deformations build up while the streamwise vorticity grows. $\unicode[STIX]{x1D714}_{x}$ growth has consequent impacts on the surface dynamics. In order to understand this, the liquid surface has been compared at two instances for both low and high density ratios in figure 51. The boxes in this figure show the computational domain edges. The boundary remains far away from the liquid surface; hence the results are not influenced by the domain size.

Figure 51.

Liquid–gas interface at $t=15~\unicode[STIX]{x03BC}\text{s}$ (a,b) and $18~\unicode[STIX]{x03BC}\text{s}$ (c,d) for $Re_{l}=2500$ , $We_{l}=14~400$ , $\hat{\unicode[STIX]{x1D707}}=0.0066$ ; $\hat{\unicode[STIX]{x1D70C}}=0.05$ (a,c), $\hat{\unicode[STIX]{x1D70C}}=0.5$ (b,d).

The higher $\unicode[STIX]{x1D714}_{x}$ growth rate at high density ratio causes the liquid surface to undergo 3D instabilities faster, and there are more streamwise lobes seen at high $\hat{\unicode[STIX]{x1D70C}}$ than at low $\hat{\unicode[STIX]{x1D70C}}$ . At $15~\unicode[STIX]{x03BC}\text{s}$ , the surface of the sheet with $\hat{\unicode[STIX]{x1D70C}}=0.05$ is still roughly 2D, while the high density-ratio case manifests more 3D deformations, and streamwise lobes are apparent on top of the primary KH waves. On the other hand, the higher $\unicode[STIX]{x1D714}_{y}$ compared to $\unicode[STIX]{x1D714}_{x}$ at low $\hat{\unicode[STIX]{x1D70C}}$ causes the liquid sheet to become antisymmetric much faster (compare the images in figure 51 a,b). Recall that as explained in § 3.2, transition towards antisymmetry expedites when the two vorticity layers – on top and bottom surfaces – become stronger, as $\unicode[STIX]{x1D714}_{y}$ grows.

The difference in the vorticity dynamics between the two cases has significant effects on the characteristics of the jet instabilities. As can be seen at $t=18~\unicode[STIX]{x03BC}\text{s}$ (figure 51), the low $\hat{\unicode[STIX]{x1D70C}}$ case can be characterized by roll-up of the KH waves, which creates more spanwise-aligned liquid structures and fewer stretched lobes; the entire sheet thins faster, and the liquid sheet breaks sooner. On the other hand, at high $\hat{\unicode[STIX]{x1D70C}}$ , the liquid structures orient streamwise more and manifest more lobes. The lobes are more stretched and thinned due to the larger $\unicode[STIX]{x1D714}_{x}$ of the KH vortex legs and are more prone to perforation. Hence, the hole-formation mechanism is expected to prevail over a larger area in the parameter space of $We_{l}$ versus $Re_{l}$ , at higher density ratios. This is consistent with the zones in figure 1, and is in accordance with the density-ratio effects discussed by Zandian et al. (Reference Zandian, Sirignano and Hussain2017).

The top and bottom liquid surfaces tend towards an antisymmetric mode, whether we start with a flat surface or symmetric perturbations. The antisymmetric behaviour is eventually favoured since a planar jet is more unstable to the antisymmetric mode than the symmetric mode in the parameter range of interest. The transformation towards antisymmetry occurs sooner as the density ratio is lowered.

Marmottant & Villermaux (Reference Marmottant and Villermaux2004) ‘suggest’ that, for a coaxial round jet, the transverse (azimuthal) deformation of the wave crests depends on surface tension ( $We_{g}$ ), density ratio ( $\hat{\unicode[STIX]{x1D70C}}$ ), and thickness of the vorticity layer. They propose the possibility that the transverse instability of the wave crest is the result of lobe and ligament formation due to opposing shear and surface tension. Unsteady motions at the sheet rim confer transient accelerations to the liquid perpendicular to the liquid–gas interface, which trigger a RT type of instability, producing indentations of the rim, which later result in ligaments.

In a physically different configuration for the leading rim of a transient liquid sheet, Agbaglah, Josserand & Zaleski (Reference Agbaglah, Josserand and Zaleski2013) also propose that the dynamics of a receding liquid sheet initiates a RT instability due to surface tension. They show that the growth rate of this transverse RT instability increases as the liquid sheet decelerates. Hence, they suggest that this instability is dominant at low $We_{g}$ , where the deceleration is strong enough so that the retraction forces overcome the liquid–gas inertia. Our mechanisms for the lobe and ligament formation – via vortex interactions – would not be in effect for a configuration where the oppositely oriented hairpin pairs do not exist – such as the leading rim of a liquid sheet studied by Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013). Our results here do not indicate that capillary action plays a major role in the deformation of the crest rim to create 3D structures, e.g., lobes, corrugations and ligaments. Once ligaments are formed, capillary action becomes important.

3.6 Vortex dynamics of round jets

In this section, a qualitative comparison is made between the vortex dynamics of the planar jets (studied here) and round jets (from previous studies) to show that the causes of different liquid structures, e.g. lobes, holes, corrugations and ligaments, are the same from the vortex dynamics perspective. For this purpose, the round liquid-jet computational analysis of Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) and Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) are mainly incorporated, which are most pertinent to our vortex dynamics analysis. In a few cases, numerical results of Shinjo & Umemura (Reference Shinjo and Umemura2010) are also addressed.

Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) used the vorticity dynamics to explain the formation of lobes and ligaments in their computations. They mainly used the vorticity contours and the vortex lines projection on the liquid surface in their analysis. Even though the explanation for the hole and corrugation formation was less detailed in their analysis, the vortex dynamics related to lobe formation and ligament elongation in their round jets is very similar to our findings. Similar to the planar jets, Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) and many researchers before them (Martin & Meiburg Reference Martin and Meiburg1991; Liepmann & Gharib Reference Liepmann and Gharib1992; Brancher et al. Reference Brancher, Chomaz and Huerre1994; Shinjo & Umemura Reference Shinjo and Umemura2010) observed the origination of the streamwise vorticity in the braid region of the round jet first. Later, as the vortices stretch in the streamwise direction, counter-rotating streamwise vortex pairs are formed. These are shown to be 3D hairpin vortices that wrap around the KH rollers. The streamwise vorticity in the ring originates from the upstream braid adjacent to the ring. The hairpin vortex in the ring region is $180^{\circ }$ (half wavelength) out of phase with respect to that of the braid region, similar to what we observed in our planar sheet. The lobe locations on the cone crests is correlated with locations of the hairpin vortices in the ring (crest) region.

A very similar mechanism involving hairpin vortices is proposed by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) for the formation of ligaments. They show that the streamwise vorticity projection on the ligament surface changes sign inside the ligament. Therefore, counter-rotating vortex structures exist both near the lobes in the gas phase and inside the elongated ligament. The side-jet phenomenon in homogeneous jets is similar to the formation of the ligaments in two-phase jet flow based on these hairpin vortical structures. The effects of the pressure gradient on the ligament elongation and breakup is also addressed by Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014). The absolute value of the pressure decreases towards the centre of the ligament at the smallest cross-sectional area; however, the pressure decreases radially outward at the tip of the ligament. Therefore, there is an oscillation of the pressure gradient inside a ligament that produces a waviness on its surface. According to Shinjo & Umemura (Reference Shinjo and Umemura2010), capillary waves in the short-wave mode propagate inside the ligament and decrease the diameter of the ligament locally. Surface tension pinches the ligament at minimum cross-sectional area. However, the large streamwise vorticity observed inside the ligament might produce a large angular momentum and compete with the capillary force. Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) define the competition between the vorticity effects and the capillary forces through a local Weber number, which is the ratio of the radial pressure gradient due to the transverse velocity difference caused by the streamwise vorticity inside the ligament and the radial pressure gradient due to surface tension, i.e.  $We=\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D714}_{x}^{2}R^{3}/\unicode[STIX]{x1D70E}$ ; where $R$ is the radius of the ligament. These two pressure gradients have different signs. If this $We$ is much larger than unity, the vorticity effects inside the ligament can be more important than the capillary force. In computations of Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014) $We$ is very close to unity, hence inertia can be as important as capillary forces in the ligament breakup.

Figure 52.

Isosurface of the liquid–gas interface accompanied by schematic of hairpin vortices on the ring and braid regions (a); magnified lobe with projected vortex lines, showing two families of hairpin vortices from the braid and crest and their pairing at the centre of the lobe leading to the lobe tearing and hole formation (b). $\hat{\unicode[STIX]{x1D70C}}=0.1$ , $Re_{l}=1600$ , $We_{g}=23\,000$ . The solid and dashed lines show hairpin portions stretching downstream and upstream, respectively. Gas flows from right to left. Recast from Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016).

In a later study, Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) described the mechanism of hole and bridge formation in a round jet using the vorticity dynamics. They arrived at the same conclusion of hairpin overlapping for hole formation as described herein (figure 52). However, the spanwise measure is replaced by the azimuthal angle for round liquid jets. At high $Re_{l}$ of interest, inertial effects dominate and vortex lines and material lines are almost identical. Thus, the hairpins are stretched around the braid and curled at the front of the crest following the fluid motion. In the vicinity of the lobe, the downstream stretching hairpin from the ring passes over the upstream stretching hairpin from the braid as the curling action continues, causing the phase variation sketched in figure 52(a). Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) also observed that the hairpin vortices approach each other and form a diamond-shaped region of the type shown by Comte et al. (Reference Comte, Lesieur and Lamballais1992), very similar to the vortex lattice observed experimentally for single-phase planar mixing layers (see figures 3 and 4 of Comte et al. Reference Comte, Lesieur and Lamballais1992). The hole-formation process relates to the increase in circulation originated from the hairpin vortices that envelope the lobe and make it thinner with time, as was explained earlier. Hole formation begins on a lobe sheet where the $We$ , based on lobe thickness and relative gas–liquid velocity, is too large for capillary action to be the initiating mechanism. As discussed earlier (figure 11), where two hairpin vortices with same-sign circulation overlap, the liquid sheet between them becomes thinner because their mutual induction moves the material lines closer to the same radial position and a hole can form (figure 52 b).

Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) express that, although the hole location can be predicted by hairpin overlap, other flow parameters, e.g. density ratio, viscosity ratio and surface tension, play a role in changing the flow details and the hole-formation process. For example, when the surface area of the lobe increases and its edge curvature decreases, the locations of the holes in neighbouring lobes will be closer to each other, and the holes merge to create larger holes (Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016).

The hole formation and ligament creation from the extension of the holes and tearing of the rim were observed in the computations of Shinjo & Umemura (Reference Shinjo and Umemura2010) for the same range of $Re$ and $We$ numbers as Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016). However, the lobe puncture was stated to be due to the impact of the droplets that formed earlier from the breakup of the mushroom-shaped cap on the jet core. Hence, their hole formation was inertially driven not vortically. Our temporal study produces no cap, but better represents the behaviour of spatial development after the cap has moved far downstream.

Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) relate the formation of small-scale corrugations to the smaller and less-orderly alternating streamwise vorticity regions near the lobes for lower $Re_{l}$ . However, they do not present a detailed explanation of the corrugation formation and stretching, nor do they explain the reason for such less-orderly vortical structures near the ring region. Hence, even though our corrugation formation for planar jets is also related to smaller scales of the hairpin filaments (at high $Re_{l}$ ), this mechanism cannot be directly compared with their observations.