2.1 Problem description

In the laboratory, shock tubes are most often used to generate large relative velocities between the liquid droplet and surrounding gas flow. Normal shock waves, in and of themselves, have little effect on the droplet. Instead, they serve as a reliable and repeatable way to create a high-speed flow around the droplet, which is responsible for the droplet’s subsequent deformation and disintegration (Hanson et al. Reference Hanson, Domich and Adams1963; Ranger & Nicholls Reference Ranger and Nicholls1968; Joseph et al. Reference Joseph, Belanger and Beavers1999). In figure 1, a schematic is shown of the initial condition and computational grid. The shock wave of strength $M_{s}=1.47$ (modelled as a step discontinuity in fluid properties) is travelling in air towards the water droplet, which is initialized as a spherical interface with diameter $D_{0}$ on the cylindrical grid. The drop and the ambient air downstream of the shock are initialized at rest. The entire 3D flow field is simulated with no enforced symmetries. Non-reflective boundary conditions (NRBC), applied at all computational domain boundaries, approximately extend the surrounding air to infinity. The simulation is performed on the computational domain $\unicode[STIX]{x1D6FA}=[-7D_{0},15.5D_{0}]\times [0,6D_{0}]\times [0,2\unicode[STIX]{x03C0}]$ with a spatial resolution of $(N_{z},N_{r},N_{\unicode[STIX]{x1D703}})=(800,600,320)$ grid cells. The grid is stretched towards the axial boundaries using a hyperbolic tangent function. The most refined portion of the grid is located near the initial location of the droplet and in the region of the near-field wake. In this region, the nominal axial and radial grid resolution corresponds to 100 cells per original droplet diameter. The azimuthal resolution is chosen such that the cells near the spherical droplet interface are close to regular. Previous testing of grid resolution sensitivities in 2D (Meng & Colonius Reference Meng and Colonius2015; Meng Reference Meng2016) suggests the present spatial resolution captures the salient flow features without being computationally cumbersome. The simulation is performed with a constant Courant–Friedrichs–Lewy (CFL) number of 0.2. Since the initial condition is axisymmetric, the gas is initially seeded with small random radial and azimuthal velocity perturbations, the largest of which have approximate magnitudes of $O(10^{-4}u_{s})$ , where $u_{s}$ is the post-shock gas velocity. These white-noise-type perturbations are generated via the Fortran compiler’s intrinsic random number generator, and are applied to both the pre-shock and post-shock gas in the initial condition.

Figure 1.

Schematic of the 3D initial condition and non-uniform computational domain at 1:10 of the actual resolution. The grid extends radially outward for $6D_{0}$ , and the axial extents are $-7D_{0}\leqslant z\leqslant 15.5D_{0}$ .

2.2 Governing equations

We model the flow with the compressible multicomponent Euler equations. In addition to being compressible, each fluid is considered immiscible and does not undergo phase change. In the absence of mass transfer and surface tension, material interfaces are simply advected by the local flow velocity. For a system of two fluids, this gives rise to a five-equation model, introduced in its inviscid form by Allaire, Clerc & Kokh (Reference Allaire, Clerc and Kokh2002). Implications of the choice to neglect viscous and capillary effects are addressed in detail in § 3.2. The five-equation model, (2.1), consists of individual continuity equations for each of the fluids, (2.1a ) and (2.1b ), mixture momentum, (2.1c ), and energy, (2.1d ), equations, and a transport equation for the gas volume fraction, (2.1e ):

where $\unicode[STIX]{x1D70C}$ is the density, $\unicode[STIX]{x1D6FC}$ is the volume fraction, $\boldsymbol{u}$ is the velocity vector, $p$ is the pressure, $E$ is the total energy defined as $E=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D700}+(1/2)\unicode[STIX]{x1D70C}\Vert \boldsymbol{u}\Vert ^{2}$ , $\unicode[STIX]{x1D700}$ is the specific internal energy, and the subscripted $g,l$ represent, respectively, the gas and liquid fluids. Each equation in the five-equation model, (2.1), is evolved independently of all others. While alternate models for multicomponent flows exist within the literature (e.g. Kapila et al. Reference Kapila, Menikoff, Bdzil, Son and Stewart2001; Murrone & Guillard Reference Murrone and Guillard2005; Saurel, Petitpas & Berry Reference Saurel, Petitpas and Berry2009; Pelanti & Shyue Reference Pelanti and Shyue2014), the present model is chosen for its simplicity and desirable conservation properties. The simplified interface model theoretically ensures volume fraction positivity, and does not require regularization to fully specify shock jump conditions. Though it is not capable of physically handling mixture regions, the model is sufficient for the immiscible fluids assumed in the multicomponent flows of interest.

The simulation of material interfaces is made possible by the volume-of-fluid method, which belongs to the broader class of interface-capturing schemes. One common characteristic of these schemes is the relaxation of the natural sharpness of material discontinuities. That is, the interfaces are allowed to numerically diffuse, resulting in an interface region of small, but finite, thickness. Within this interface region, a non-physical fluid mixture exists, which is appropriately treated using mixture rules that are defined in § 2.2.1. For numerical stability purposes, material interfaces are not initialized as sharp discontinuities, but are smeared over a few grid cells. Based on previous results (Johnsen Reference Johnsen2007; Coralic Reference Coralic2015) and sensitivity testing by the authors for the specific problem of aerobreakup, this initial artificial smearing over a few cells is known to have negligible impact on the computational results.

3.1 Spatial and temporal discretizations

The five-equation model, (2.1), is spatially discretized on a 3D cylindrical grid in the following form:

where $\boldsymbol{q}$ is the vector of conservative variables, $\boldsymbol{f}(\boldsymbol{q}),\boldsymbol{g}(\boldsymbol{q}),$ and $\boldsymbol{h}(\boldsymbol{q})$ are flux vectors in each of the coordinate directions, and $\boldsymbol{s}(\boldsymbol{q})$ is a source term vector (see details in appendix A). In cylindrical coordinates, the divergence operator on an arbitrary vector $\boldsymbol{v}=(v_{z},v_{r},v_{\unicode[STIX]{x1D703}})^{\text{T}}$ is

Following Johnsen (Reference Johnsen2007), all $v_{r}/r$ terms are treated as geometrical source terms in $\boldsymbol{s}(\boldsymbol{q})$ . By discretizing the governing equations in this manner, the source terms account for all cylindrical geometry effects. The axis singularity is treated following Mohseni & Colonius (Reference Mohseni and Colonius2000) with the definition of a new radial coordinate that spans both positive and negative radius, i.e.

Additional details of the numerical method, as well as its verification and validation, can be found in Meng (Reference Meng2016). Within the finite-volume formulation, the reconstruction of the state variables at the cell boundaries is performed using a third-order weighted essentially non-oscillatory (WENO) scheme developed for Cartesian coordinates. Using the Cartesian WENO weights and polynomials, formal order of accuracy is retained in the axial and radial coordinates, while not guaranteed in the azimuthal coordinate. Convergence studies of the numerical method, however, show that second-order accuracy is retained (Meng Reference Meng2016). Finally, the system of equations is temporally integrated using a third-order total variation diminishing Runge–Kutta scheme.

3.2 Modelling approximations and implications

The governing equations of § 2.2 are notably missing physical models for molecular viscosity and surface tension. While the algorithms required to implement viscous and capillary models in the governing equations are known (Meng Reference Meng2016) (the five-equation model of Allaire et al. (Reference Allaire, Clerc and Kokh2002) was subsequently extended to include viscous effects by Perigaud & Saurel (Reference Perigaud and Saurel2005), and capillary effects can be modelled by including additional source terms), the computational resources required for such a direct numerical simulation are well beyond what is currently feasible. The present spatial resolution of $(N_{z},N_{r},N_{\unicode[STIX]{x1D703}})=(800,600,320)$ is already a grid with $1.536\times 10^{8}$ cells. For the high Reynolds numbers associated with aerobreakup in the SIE regime, the grid resolution requirements to reach grid convergence would be extraordinary. As an example, the axisymmetric work of Chang et al. (Reference Chang, Deng and Theofanous2013) employed an eight-level adaptive mesh refinement (AMR) scheme with an equivalent resolution of 800 points per diameter. In order to sufficiently resolve the viscous boundary layer, the AMR mesh had to be augmented with an additional body-fitted structured conformal mesh that had a spatial resolution of 4000 points per diameter. Taking even the coarsest resolution requirement of 800 cells per diameter would be an increase by a factor of 8 from the present resolution. Additionally, while it is more difficult to estimate a precise grid requirement for capillary effects, a conservative estimate of a factor of 10 in each coordinate direction would result in a requirement of $O(10^{14})$ grid points to reach ab initio physical fidelity. Even with AMR savings, we suspect this problem to be intractable. Therefore, while the physics can be (and has been) implemented with relative ease, the numerical viscosity would dominate these physical effects unless the above spatial resolutions were somehow achieved.

Given that a compromise must be made, we consider the ‘inviscid’ case without surface tension. Firstly, it is generally understood that so-called ‘inviscid’ simulations using shock- and interface-capturing schemes, i.e. those performed without explicit molecular viscosity modelling, inherently include numerical viscosity whose magnitude is dependent on the spatial resolution. This compromise thus precludes providing evidence of convergence, because such convergence will not exist when the flow becomes unstable and turbulent – the smallest length scales will be at the scale of the grid. It is clear that at the present spatial resolution, we are unable to capture fine-scale instabilities that might arise on the droplet’s surface such as the KH waves experimentally observed by Theofanous et al. (Reference Theofanous, Mitkin, Ng, Chang, Deng and Sushchikh2012) (recall that Chang et al. (Reference Chang, Deng and Theofanous2013) used a resolution of 4000 nodes per diameter to capture the viscous boundary layer). It should be noted that an inability to capture these fine-scale instabilities is solely a consequence of the employed spatial resolution. Interfacial instabilities (at least KH, RT, Richtmyer–Meshkov, and capillary) are driven by shear, which arises even in ‘inviscid’ simulations. For a solid surface, the no-slip condition, which of course requires the existence of viscosity, gives rise to shear, which in turn can drive both viscous (e.g. Tollmien–Schlichting) and inviscid instabilities. In the interface between two immiscible fluids, the conservation of momentum requires equality of the tangential viscous stresses, which in turn leads to a continuity of the velocity. In the true inviscid limit, the fluids would slip, but in the presence of any viscosity, and, indeed, any artificial viscosity, the velocity is again continuous. This leads to shear and, in turn, interfacial instabilities. With the exception of capillary instabilities, the five-equation model can support all the other mentioned surface instabilities, though they are not always adequately resolved.

The effects of surface tension in the aerobreakup problem become increasingly important with time. In particular, capillary instabilities are of primary importance when it comes to the disintegration of the thin sheets of liquid that are stripped from the droplet’s periphery (discussed in § 4). In the absence of surface tension modelling, there does not exist a numerical mechanism that approximates capillary effects, i.e. there is no capillary counterpart to numerical viscosity. Without such a mechanism, we do not capture the ultimate capillary-driven breakup mechanism of a droplet undergoing aerobreakup. Instead, the actual mechanism of breakup in the results is modelled in terms of numerical diffusion effects, i.e. breakup occurs as a consequence of numerical diffusion and finite spatial resolution. Given this modelling approximation, we expect the numerical results to have increasingly large uncertainties as time progresses. Nevertheless, the agreement we find with experimental results suggests that the important mechanisms are captured; a definitive quantification of the error awaits future work.

While traditional grid convergence cannot be shown for the present numerical results, we can estimate an approximate, effective Reynolds number associated with the numerical viscosity for the given grid. A series of tests using 2D viscous simulations (Meng & Colonius Reference Meng and Colonius2015; Meng Reference Meng2016) indicate that the effective Reynolds number is no less than 500 at the spatial resolution of our 3D simulation. While this estimate is crude, and much lower than the corresponding experiments, it provides some assurance that we are in the inertia-dominated regime, where we might expect viscosity to play a small role at the scales of interest, at least at early times before the flow is fully turbulent and the droplet is disintegrated.

3.3 Droplet diagnostics

As the droplet undergoes aerobreakup, its centre-of-mass properties are of interest. Taking advantage of the type of quantitative analysis allowed by simulations, integral expressions have been derived (Meng & Colonius Reference Meng and Colonius2015) for the droplet’s centre-of-mass velocity and acceleration that minimize unnecessary noise that would be introduced by differentiating position data:

where the integrated volume is that of the entire computational domain, $\unicode[STIX]{x1D6FA}$ . The liquid partial density, $\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}$ , is then the parameter that restricts the integration to cells with non-zero liquid volume fractions. Additional details of the derivations of (3.4)–(3.6) can be found in Meng (Reference Meng2016). It should be noted that (3.5) and (3.6) are valid as long as the liquid mass in the domain, $m_{l}=\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\,\text{d}V$ , remains a constant; once liquid mass flux through the domain boundaries is non-zero, we terminate their calculation.

Given that our volume-of-fluid numerical method involves a diffuse interface that is smeared across a few grid cells, isopleths and isosurfaces, representing the coherent droplet body in the following numerical results, are shown for multiple values of the liquid volume fraction, $\unicode[STIX]{x1D6FC}_{l}$ . The impact of various choices of $\unicode[STIX]{x1D6FC}_{l}$ on the interpretation of the numerical results is explicitly discussed in the respective results sections. Fortunately, the ambiguity in the exact interface location does not adversely impact the conclusions that are drawn from our results.

4.1 Centre-of-mass properties

The droplet’s centre-of-mass drift, velocity, and acceleration in the streamwise (axial) direction are plotted in figure 4. While the drift curve appears to be roughly parabolic, the subsequent velocity and acceleration curves reveal that a constant acceleration assumption would be erroneous. Focusing on the acceleration curve, the initial spike in acceleration is the passage of the shock wave over the droplet. The maximum acceleration occurs when the shock reflection on the droplet’s surface transitions from a regular reflection to a Mach reflection. This is further discussed in § 4.3. Following this initial spike, we note the existence of a brief period immediately following the passage of the shock when the droplet is subject to constant acceleration. During this time period, the droplet is adjusting to the step change in ambient flow conditions, and is still well-approximated as a rigid sphere. This delay, which is not observed in 2D simulations (Meng & Colonius Reference Meng and Colonius2015) and is, perhaps, related to the flow-relieving effect of the third dimension, ceases when the droplet begins to pancake and its drag properties substantially change. Secondly, a dominant low-frequency oscillation is observed in the droplet’s acceleration curve. In figure 5, we replot the acceleration curve using standard convective time units for flow past a rigid sphere, $tu_{s}/D_{0}$ , to check if this low frequency corresponds with wake instability. The vertical gridlines are spaced $5D_{0}/u_{s}$ apart to coincide with the expected period of a wake instability (vortex shedding) which, for flow over a rigid sphere, occurs with Strouhal number $St=fD_{0}/u_{s}=0.2$ . The largest oscillations in the acceleration are roughly commensurate with this frequency, but owing to the relatively short data available, it is difficult to draw a firm conclusion.

Figure 4.

Droplet streamwise centre-of-mass (a) drift, (b) velocity, and (c) acceleration.

Figure 5.

Droplet (a) acceleration in convective time units and (b) unsteady drag coefficient for a range of $D_{d}$ (dependent upon choice of $\unicode[STIX]{x1D6FC}_{l}$ ).

Finally, figure 5(b) shows the droplet’s unsteady drag coefficient, $C_{D}$ . The frontal area used in the non-dimensionalization assumes a circular upstream projected area based on the droplet’s deformed diameter, $D_{d}$ . From figure 2, this assumption is seen to be a reasonable approximation for the times shown in figure 5. $D_{d}$ is taken as the maximal spatial extents in the $\tilde{r}$ -coordinate, which depends on the choice of $\unicode[STIX]{x1D6FC}_{l}$ . Thus, in figure 5, we show bounds for the $C_{D}$ spanning $0.25\leqslant \unicode[STIX]{x1D6FC}_{l}\leqslant 0.99$ . For $t^{\ast }<0.3$ , the choice of $\unicode[STIX]{x1D6FC}_{l}$ has little impact on the $C_{D}$ . During the period of constant acceleration, and before the droplet has had sufficient time to significantly deform, we observe that the drag coefficient for a rigid sphere, $C_{D}=0.5$ , is approximately recovered. As the droplet begins pancaking, the drag coefficient transitions to be comparable to that for a flat disk, which has $C_{D}\approx 1$ . For $t^{\ast }>0.3$ , $C_{D}$ depends more strongly on choice of $\unicode[STIX]{x1D6FC}_{l}$ -value, representing the uncertainty in the deformed diameter. Unsteady effects begin to dominate and the drag coefficient subsequently exhibits fluctuations about an increasing average value. While it is impossible to state the true value of $C_{D}$ , the oscillatory behaviour (present also in the acceleration time history) is believed to be physical, and the plot is thought to reasonably bound the physical values at late times.

4.2 Comparison with experimental visualizations

We now compare the simulated breakup with the experimental visualizations of Theofanous et al. (Reference Theofanous, Mitkin, Ng, Chang, Deng and Sushchikh2012). They studied the breakup of water droplets (among other liquids) in a helium shock tube. Given the uncertainty in the exact interface location and the numerical diffusion of the material discontinuity, we compare the experimental images to both the $\unicode[STIX]{x1D6FC}_{l}=0.01$ isosurface and sliced isopleths for various $\unicode[STIX]{x1D6FC}_{l}$ -values. Despite a mismatch in flow conditions and a dearth of timing data, our numerical results show good qualitative agreement with the experimental visualizations of the SIE phenomenology. The small value of the $\unicode[STIX]{x1D6FC}_{l}$ -isosurface is believed to be a fair comparison for visualization purposes, as the experimental images are also obscured by the fine mist that is generated during breakup. Additionally, comparisons with the sliced isopleths confirm a qualitative agreement largely independent of the numerics.

Figure 6.

Comparisons with experimental visualizations of the aerobreakup of a water droplet from figure 33 of Theofanous et al. (Reference Theofanous, Mitkin, Ng, Chang, Deng and Sushchikh2012). Flow is from right to left. Each frame consists of an experimental image reproduced from the video stored online (https://doi.org/10.1063/1.3680867.12) (left), the numerical $\unicode[STIX]{x1D6FC}_{l}=0.01$ isosurface (centre), and sliced isopleths of $\unicode[STIX]{x1D6FC}_{l}$ (right) corresponding to the legend in figure 3. The timing information is for the numerical results; no experimental timing information is available.

The post-shock flow in the experiment has a Mach number of 0.32, which corresponds to experimental Reynolds and Weber numbers of $Re=2.2\times 10^{4}$ and $We=780$ . The large Reynolds and Weber numbers allow us to make qualitative comparisons with our numerical results. This comparison is shown in figure 6. First, we observe the same initial deformation of the droplet into a muffin-like shape. The upstream side of the droplet remains spherical, while the downstream side is flattened into a planar surface. What appears to be a thin liquid sheet coming from the spherical lip is visible beginning at $t^{\ast }=0.162$ . This sheet quickly disintegrates into a mist that obscures the coherent part of the droplet. The coherent droplet body is continually flattened in the streamwise direction, and liquid material is constantly stripped off near its equator. At late times in the SIE process, the liquid fragments form a cylindrical curtain around a cavity behind the coherent droplet. In figure 7, we compare experimental and numerical images of this liquid sheet. As discussed in § 3.2, we do not capture the ultimate capillary-driven breakup of the liquid sheet. Though these experimental comparisons are not without uncertainty, we find that our numerical results are, by and large, not far off from the experimental visualizations of the SIE regime.

Figure 7.

Comparison of experimental and numerical liquid ‘curtains.’ The experimental image taken from figure 17(b) of Theofanous et al. (Reference Theofanous, Mitkin, Ng, Chang, Deng and Sushchikh2012) (a) is for Mach 0.59 flow at $t^{\ast }=0.81$ . The numerical isosurface of $\unicode[STIX]{x1D6FC}_{l}=0.01$ (b) and sliced isopleths of $\unicode[STIX]{x1D6FC}_{l}$ (c) are for Mach 0.5775 flow at $t^{\ast }=0.808$ . Isopleth colours correspond to the legend in figure 3.

4.3 Flow field in the gas

From their experimental results, Liu & Reitz (Reference Liu and Reitz1997) proposed a rough classification of the breakup process into two stages. They experimentally observed the first stage as a pressure-driven shape change of the droplet. This deformation stage, shared amongst multiple breakup regimes, is characterized by a flattening of the droplet in the streamwise direction. This pancaking is a consequence of the non-uniform pressure distribution around the droplet surface. High pressures at the forward and rear stagnation points, as well as low pressures at the equator due to the acceleration of the gas, both contribute to the flattening of the droplet. The second stage, as described by Liu & Reitz (Reference Liu and Reitz1997), is characterized by droplet disintegration and is when the phenomenology diverges for the various breakup regimes. For the stripping regime, they observed the edges of the droplet being drawn out into a thin liquid sheet by drag forces, and the subsequent breakup of the sheet into fine ligaments.

Contrary to Liu & Reitz’s description of the breakup process as two consecutive stages, our numerical results suggest that the phenomenology of stripping may be better described as the simultaneous flattening and stripping of liquid material. In fact, not only are the flattening and the disintegration processes occurring concurrently throughout the breakup process, they are also intricately connected by the dynamic behaviour of the surrounding gas flow.

Figure 8.

Filled contour slices of velocity magnitude, $\Vert \boldsymbol{u}\Vert ^{\ast }$ , and pressure, $p^{\ast }$ , and isopleths of the numerical schlieren function. Flow is from bottom left to top right. Offset 2D slices are taken at $x^{\ast }=0,y^{\ast }=0$ .

Figure 9.

Isosurfaces of positive (red) and negative (blue) azimuthal vorticity, $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}$ . Flow is from top left to bottom right. Offset 2D slices are taken at $x^{\ast }=0,y^{\ast }=0,z^{\ast }=1$ .

In order to elucidate why the liquid behaves as it does, we look to the behaviour of the surrounding gas flow as visualized in figures 8 and 9 (images are ordered top to bottom, left to right). Figure 8 shows 2D slices taken through the centre of the droplet, $x^{\ast },y^{\ast }=0$ , that are offset for unobstructed viewing of both planes. The vertical plots are coloured by velocity magnitude normalized by the post-shock gas velocity, $\Vert \boldsymbol{u}\Vert ^{\ast }$ , while the horizontal plots are coloured by pressure, $p^{\ast }$ . Isopleths of the numerical schlieren function, $\unicode[STIX]{x1D711}$ , reveal the intricate and dynamic flow structures that develop in the wake. Following Quirk & Karni (Reference Quirk and Karni1996), the numerical schlieren function is computed as the exponential of the negative, normalized density gradient

where $\unicode[STIX]{x1D6FD}$ is a scaling parameter that allows simultaneous visualization of waves in both fluids. Following Johnsen (Reference Johnsen2007), $\unicode[STIX]{x1D6FD}_{air}=40$ and $\unicode[STIX]{x1D6FD}_{water}=400$ . The three-dimensionality of the aerobreakup process is well captured in figure 9, which plots various isosurfaces of azimuthal vorticity, $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}$ . Since the innermost regions of the flow are obstructed from view, slices are again taken through the centre of the drop and offset to both sides. The transverse slice, offset upstream of the droplet, is taken at $z^{\ast }=1$ .

The incident and reflected shock waves, as well as the secondary wave system generated by the convergence of Mach stems at the rear stagnation point, are visible in the first few snapshots of figure 8. The peak drag of the droplet is marked by the transition of the reflected shock from a regular reflection to a Mach reflection at some inclination angle preceding the droplet equator. This phenomenon has been studied in the literature for rigid cylinders and spheres (Takayama & Itoh Reference Takayama and Itoh1986; Tanno et al. Reference Tanno, Itoh, Saito, Abe and Takayama2003), and is also observed in the case of 2D aerobreakup (Meng & Colonius Reference Meng and Colonius2015). Promptly after the passage of the incident shock, the flow is accelerated to approximately $1.5u_{s}$ at the droplet’s equator (pink colouring is visible at $t^{\ast }=0.044$ in figure 8), which is the value expected from potential flow theory for flow past a rigid sphere. In the early stages of droplet deformation, $t^{\ast }\leqslant 0.435$ , the transmitted wave that propagates into the liquid from the incident shock wave is observed to bounce back and forth within the droplet. From our numerical results and the aerobreakup literature, this transmitted wave inside the liquid droplet is not thought to play a significant role in the aerobreakup process.

As discussed earlier, the non-uniform pressure distribution around the droplet is the principal mechanism driving droplet deformation. This is clearly seen in the filled pressure contours of figure 8. Unlike the case of steady separated flow past a rigid sphere, the pressure at the rear stagnation point remains, during the initial flattening, larger than the pressure at the droplet’s equator.

The equatorial recirculation region, visible in figure 9, is formed by the interaction of two opposite-sign azimuthal vorticity streams generated by baroclinicity. The formation of this equatorial recirculation region has previously been observed to also occur for 2D aerobreakup (Meng & Colonius Reference Meng and Colonius2015). Its location coincides with the sides of the muffin-shaped droplet behind the lip of the spherical upstream droplet surface. This equatorial recirculation region thus serves as a possible explanation for the muffin-like shape of the deformed droplet, and, for the duration of its existence, is at least partially responsible for the liquid sheets that are drawn out from both the spherical lip and the planar back of the droplet.

Other notable flow features visible in figure 8 include a wake recirculation region and two concentric shear layers at the droplet equator. Behind the droplet, the wake recirculation region is quickly established, as evidenced by the departure of the flow reversal region (i.e. the white patch in the wake directly behind the droplet where the streamwise velocity changes direction) from the rear stagnation point at $t^{\ast }=0.099$ in figure 8 (at $t^{\ast }=0.126$ in figure 9, the wake recirculation region already exists). The flow reversal patch also serves to demarcate the end of an upstream jet that is created in the wake. The wake recirculation region is bounded by shear layers that are subject to KH instability, and the unsteady vortex shedding drives the development of a complicated wake. KH roll-up and subsequent vortex shedding from the shear layers can be seen particularly clearly from the numerical schlieren isopleths at $t^{\ast }=0.435,0.544,0.781$ in figure 8.

For clarity, it would be ideal if the aforementioned flow phenomena could be discussed separately as distinct flow features. However, in reality, they are so interconnected that such a disjointed discussion would be both incomplete and overly simplistic. In the remainder of this section, we attempt a comprehensive examination and discussion of these flow phenomena.

The formation of the wake recirculation region initiates a strongly coupled, self-sustaining set of flow phenomena that evolve with ever-increasing complexity. The wake recirculation region, created by the stream of negative azimuthal vorticity from the upstream side of the droplet, remains in the near-field wake region. It is perpetually sustained by the same vorticity stream, and entrains the surrounding fluid, which is pulled into an upstream jet that impinges on the rear stagnation point of the droplet.

This upstream jet, driven by the recirculation region, preserves high pressure at the rear stagnation point that contributes to both the pancaking of the coherent droplet, as well as the generation of positive vorticity along the back of the droplet. The positive vorticity is transported towards the equator by the recirculation region, and interacts with the negative vorticity stream to create, at early times, the equatorial recirculation region. Some positive vorticity is also transported downstream, forming a parallel stream inside the negative vorticity stream coming from the upstream side of the droplet. The shear flow that is created when the upstream jet impinges on the back of the droplet may also contribute to the liquid sheet that arises from the downstream side.

As the equatorial liquid sheet is blown downstream by inertial forces from the surrounding gas flow, shear layers are formed on both sides of the sheet. The gas that is accelerating around the droplet creates shear on the exterior, while the wake recirculation region, that exists inside the cavity enveloped by the sheet, creates shear along the interior. Both the interior and exterior shear layers are visible in figure 8. As the liquid sheet flaps, generating longitudinal ripples, the shear layers, which are subject to KH instability, periodically shed vortices that are either entrained by the wake recirculation region, or are convected downstream. Entrained vortices (of both signs) by the wake recirculation region result in the upstream jet being characterized by concentric layers of alternating vorticity sign (visible from $t^{\ast }=0.535{-}0.781$ in figure 9). Entrainment of shed vortices is also associated with a temporary increase in upstream jet velocity that results in a cyclic pumping of fluid onto the back side of the droplet. Downstream-convected vortices, that are not entrained, quickly lose their initial axisymmetry due to the instability of vortex rings, and subsequently develop into fully 3D flow features.

In time, as the liquid sheet is drawn downstream and the coherent droplet diameter expands laterally, the flow within the enveloped cavity at the back of the droplet, encompassing the wake recirculation region and the upstream jet, correspondingly grows in size and complexity. Loss of axisymmetry is observed to first occur in the core region of wake, i.e. small $r$ . However, before it radially expands to encompass the entire wake region, another instability is observed to emerge along the still-axisymmetric positive vorticity sheet (located just inside the equatorial liquid sheet). From $t^{\ast }=0.754{-}0.781$ in figure 9, we see what appear to be RT fingers or mushroom-like features propagating inwards towards the core that cause azimuthal rippling in the previously axisymmetric vorticity sheet. This instability, further discussed in § 5, generates the transverse azimuthal modulations observable on the liquid sheet visible for $t^{\ast }\geqslant 0.808$ in figure 2. As the entire wake region devolves into chaotic, turbulent-like flow, the general coherence of the aforementioned phenomena is lost, as seen for $t^{\ast }\geqslant 1.054$ in figure 9. At these late times, the coherent droplet body presents an essentially blunt body to the oncoming free-stream flow such that the highest pressures are found on the upstream side of the flattened disk-like droplet.

5.1 Mechanism of instability

The upstream side of a droplet undergoing aerobreakup is susceptible to RT instability waves that arise from the acceleration of the lighter gas into the denser liquid. In the classical catastrophic breakup regime, ‘fingers of hot air’ were thought to penetrate the droplet, leading to an explosive disintegration (Joseph et al. Reference Joseph, Belanger and Beavers1999). In contrast stands the recent of work of Theofanous & Li (Reference Theofanous and Li2008) who observed that ‘[t]here are no RT waves piercing the drop…’ Indeed, no RT waves are seen on the droplet in our numerical simulation. This is shown in figure 10, where the upstream side of the drop for various isopleths of $\unicode[STIX]{x1D6FC}_{l}$ remains smooth for the entirety of the simulation. Our numerical results thus support Theofanous & Li’s claim of SIE being the terminal breakup regime. The suppression of the RT instability waves was initially explained by Theofanous & Li (Reference Theofanous and Li2008) to be a consequence of the stability of the stagnation flow (reminiscent of the lenticular shape of a gas bubble rising through liquid (Batchelor Reference Batchelor1987)). More recently, an analysis of the viscous KH instability by Theofanous et al. (Reference Theofanous, Mitkin, Ng, Chang, Deng and Sushchikh2012) found that ‘[w]ave numbers and growth factors of [KH] instability are consistently greater than those of [RT] instability by more than an order of magnitude…[T]he stretching further contributes to keeping this area mollified and free of any instability all the way to the end.’

Figure 10.

Sliced isopleths of $\unicode[STIX]{x1D6FC}_{l}$ showing droplet profiles.

In describing the instabilities that arise on the liquid sheet that lead to its disintegration, Liu & Reitz (Reference Liu and Reitz1997) proposed the ‘stretched streamwise ligament breakup’ mechanism of Stapper & Samuelsen (Reference Stapper and Samuelsen1990). This breakup mechanism is characterized by the dominant formation and growth of streamwise vortical waves on the sheet, with thin membranes formed between them. ‘Breakup occurs as the membranes are stretched thin by the rotation of the streamwise vortices and burst into small droplets. The streamwise vortical waves separate as streamwise ligaments, stretch and spin faster in the presence of the air shear, and eventually break up, contributing the larger drops to the final drop size distribution’ (Stapper & Samuelsen Reference Stapper and Samuelsen1990). From figure 11, which plots transverse slices of the liquid sheet at various streamwise locations, we see that our numerical results do not support this mechanism as the reason for sheet disintegration. Instead of a liquid sheet with variable thickness in the azimuthal coordinate, we observe relatively constant sheet thickness, and a rippling-type instability. Additionally, isosurfaces of the liquid sheet coloured by axial and radial vorticity, $\unicode[STIX]{x1D714}_{z,r}^{\ast }=\unicode[STIX]{x1D714}_{z,r}D_{0}/u_{s}$ , (as shown in figure 12) suggest that streamwise (axial) vorticity (with its smaller magnitude) does not play a dominant role in the sheet breakup. It should be noted here that while figures 11 and 12 are plotted using the $\unicode[STIX]{x1D6FC}_{l}=0.01$ isosurface, the sheet thickness is still relatively constant in the azimuthal coordinate for other values of $\unicode[STIX]{x1D6FC}_{l}$ (visible from figure 10), and the relative importance of $\unicode[STIX]{x1D714}_{z}^{\ast }$ is unchanged by the choice of $\unicode[STIX]{x1D6FC}_{l}$ . In addition to the ‘stretched streamwise ligament’ mechanism of sheet breakup, Liu & Reitz (Reference Liu and Reitz1997) also proposed another rippling mechanism based on mass conservation arguments. Though this remains a possibility, the observed phenomena are most likely the net result of several mechanisms. Given the modelling approximations and implications discussed in § 3.2, no physical length scales can be associated with the instabilities that are observed in the numerical results. Instead, we are limited to qualitatively describing the observed phenomena.

Figure 11.

Transverse slices of the liquid sheet at $t^{\ast }=0.881$ , defined by $\unicode[STIX]{x1D6FC}_{l}=0.01$ , at various streamwise locations.

Figure 12.

The liquid sheet coloured by axial and radial vorticity, $\unicode[STIX]{x1D714}_{z,r}^{\ast }$ , at $t^{\ast }=0.808$ .

Figure 13.

Filled contours and isopleths of gas partial density, $\unicode[STIX]{x1D6FC}_{g}\unicode[STIX]{x1D70C}_{g}$ , at $z^{\ast }=1$ . Darker colours correspond to larger densities (colouring and isopleth values vary between frames).

Jalaal & Mehravaran (Reference Jalaal and Mehravaran2014) proposed the RT instability as the source of the transverse azimuthal modulations. Arguing that the accelerated liquid sheet is subject to the same instability as that which forms streamwise ligaments in the case of a round liquid jet in coaxial flow, they attempted a quantitative comparison with theory, but found only marginal agreement. The general concept, though, of RT instability on the liquid sheet may, indeed, have merit, and supporting (qualitative) evidence can be found in our numerical results. As noted at the end of § 4.3, RT fingers or mushroom-like features are seen emerging along the positive vorticity sheet that lies just inside the equatorial liquid sheet. To relate these features to the RT instability, we plot filled contours and isopleths of the gas partial density for late times in figure 13. At $t^{\ast }=0.726$ , the axisymmetry of the outer regions of the wake is still preserved, as evidenced by the circular isopleths. From $t^{\ast }=0.744{-}0.790$ , however, the axisymmetry is broken by multiple fingers of denser gas propagating towards the axis into the lighter gas (darker colours correspond to larger densities). These fingers correspond exactly with the indentations generated on the liquid sheet. Not long after this loss of outer axisymmetry, the entire wake region degenerates into a chaotic, turbulent-like flow with complete loss of flow feature coherence.

Figure 14.

$L_{2}$ -norm of the (kinetic) energy (see (5.1)) for each of the $N_{\unicode[STIX]{x1D703}}/2$ modes resulting from the azimuthal Fourier decomposition of the velocity field.

Figure 15.

Isosurfaces of $\unicode[STIX]{x1D705}_{m},m=1,2,4,6,8$ at $t^{\ast }=0.808,0.935$ . Also shown is a cut isosurface of the liquid sheet visualized using $\unicode[STIX]{x1D6FC}_{l}=0.01$ . Front and side views are shown, and isosurface values change between frames.

5.2 Azimuthal Fourier decomposition

Motivated by the observed azimuthal modulations, we perform a Fourier decomposition of the velocity flow field to determine if a particular mode(s) or wavenumber(s) is associated with the loss of axisymmetry. To do this, we take a Fourier transform in the $\unicode[STIX]{x1D703}$ -coordinate to obtain the Fourier coefficients of each of the azimuthal modes, $\hat{\boldsymbol{u}}_{m}(z,r,t)$ . We then calculate an energy metric for each mode defined as

where the hat denotes the Fourier transform in the $\unicode[STIX]{x1D703}$ -coordinate. Taking an $L_{2}$ -norm of $\hat{\unicode[STIX]{x1D705}}_{m}$ over the entire computational domain, we plot the time histories of $\Vert \hat{\unicode[STIX]{x1D705}}_{m}\Vert _{2}$ for each of the $N_{\unicode[STIX]{x1D703}}/2$ modes (excluding the mean) in figure 14. Notably, the frequency response of the system shows broadband instability growth for all modes. An initial jolt applied to the system by the passage of the shock wave is mostly stabilized after approximately $t^{\ast }=0.1$ . This is followed by the exponential growth of nearly all wavenumbers. If we view the interaction of the incident shock with the droplet as an impulsive force applied to the system, the broadband instability response is not surprising. Four of the modes appear to exhibit a different type of behaviour; however, these unusual curves are actually artefacts of the random velocity perturbations seeded in the initial condition. Unfortunately, due to the unsteady, non-stationary, and nonlinear nature of the aerobreakup problem, this type of instability analysis is unable to pick out a dominant mode or wavenumber associated with the loss of axisymmetry. Despite the broadband response, visualization of the first few modes transformed back into $\unicode[STIX]{x1D703}$ -space, $\unicode[STIX]{x1D705}_{m}=u_{z,m}^{2}+u_{r,m}^{2}+u_{\unicode[STIX]{x1D703},m}^{2}$ , reveals some interesting observations about the development of the wake, and offers some general intuition about the physical spatial structure of $\unicode[STIX]{x1D705}_{m}$ . Figure 15 shows the isosurfaces of a few modes at late times in the simulation, $t^{\ast }=0.808,0.935$ , when all frequencies have saturated (see figure 14). From figure 15, we see that the structures are initially clustered around the location of the liquid sheet, and subsequently grow in size and complexity. Even at late times, the structures remain bounded by the wake region. Streaky streamwise-oriented flow structures near the wake core, visible for $\unicode[STIX]{x1D705}_{6,8}$ at $t^{\ast }=0.808$ , appear to be related to the loss of axisymmetry in the upstream jet region, while the outer structures are linked to the liquid sheet and shear layers. As this broadband instability is not directly associated with any classical hydrodynamic instability, we are unable to relate these $\unicode[STIX]{x1D705}_{m}$ -structures to other recognizable flow features.