On aerobreakup in the stagnation region of high-Mach-number flow over a bluff body

Abstract

In this paper, aerobreakup in the stagnation region of high-Mach-number flow over a bluff body is studied with experiments and computations. Water drops of diameter 0.51–2.30 mm were acoustically levitated at sea level along the flight path of a rectangular $100\ {\rm mm} \times 150\ {\rm mm}$ rail-gun launched projectile. This enabled the study of aerobreakup at high Mach (3.03–5.12), post-shock Mach (1.5–1.9), Weber $(5 \times 10^4\unicode{x2013}4 \times 10^5)$ and Reynolds $(6 \times 10^4\unicode{x2013}3 \times 10^5)$ numbers. High-speed backlit shadowgraphy is used to record the flow structure. Computations are made for two cases, and it was found that the drop behaviour is not dominated by viscous or surface-tension effects and can be adequately captured by treating the gas as calorically perfect with the ratio of specific heats set to 1.3 to account for thermochemical effects. To assess drop surface stability at early breakup times, results from Newton's inclination method are used to determine the flow along the drop surface and input to a linear-stability analysis; from this, it was found that viscosity and surface tension can be neglected. Moreover, the acceleration term dominates the shear term at the stagnation point, a point accentuated as a drop flattens; this relation inverts closer to the drop equator. Linear-stability analysis was insufficient for modelling late-stage aerobreakup because the predicted wavelengths were too small and the expected aerobreakup times were non-physically short. To address this discrepancy, a nonlinear instability model with constant-rate growth is used that treats the accelerated drop surface as analogous to bubbles rising through a liquid; agreement with computations is good.

1. Introduction

High-speed-vehicle design requires the consideration of many fundamental fluid-mechanic phenomena (Leyva 2017). Impacts with drops or small particles at high velocity can degrade aerospace materials, reducing transparency and other desirable properties (Lapp, Stutzman & Wahl 1954; Jenkins 1966; Reinecke 1974; Schmitt 1975; Smith 1976; Callens & Lawrence 1981; Gohardani 2011; Moylan, Landrum & Russell 2013). The sphericity and size of a drop can change the degree of damage upon impact (Adler 1975; Adler & Hooker 1978; Adler 1982, 1987, 1995a,b, 1999; Lee et al. 2012). Additionally, the flow structure about a high-speed vehicle, namely the leading shock wave, can disrupt a drop, potentially altering the deleterious effects of impact (Engel 1958; Nicholson 1967; Ranger & Nicholls 1969, 1972; Reinecke & McKay 1969; Waldman & Reinecke 1971; Simpkins & Bales 1972; Waldman, Reinecke & Glenn 1972; Barber et al. 1975; Reinecke & Waldman 1975; Barber 1976). Numerous researchers have devised drop breakup models, such as the Taylor analogy breakup (TAB) model from O'Rourke & Amsden (1987). Researchers have added to these breakup models (Beale & Reitz 1999; Schmehl 2002; Stefanitsis et al. 2019), and, more recently, have worked towards integrating these models in high-speed, multidimensional flow (Hess et al. 2021; Aguilera et al. 2023; Daniel et al. 2023).

The first phase of this problem is processing of the drop by a bow shock, the strength of which is determined by the Mach number of the projectile (here, $M_p>3$). The drop begins the aerobreakup process post-shock in a non-constant flow, juxtaposed to the near-constant flow in drop demise experiments in a shock tube. Then, the altered drop impacts a surface at hypervelocity, which results in a complex multiphase flow field. Major dimensionless groups that govern aerobreakup and drop impact that have been identified include the Weber number, the Reynolds number, the Ohnesorge number, the Mach number and the Knudsen number. The Weber number is the ratio of the aerodynamic force to the restoring capillary force on the drop and can be written as

(1.1)\begin{equation} We_A = \frac{\rho_g U_A^{2} d_{d}}{\sigma}.\end{equation}

Here, $\rho _g$, $U_A$, $d_d$ and $\sigma$ are, respectively, the local gas density, the velocity difference between the local gas and the drop, the drop diameter and the surface tension for the drop–gas interface. The drop impact and aerobreakup research communities often use different values for density and velocity in this calculation. The Reynolds number is the ratio of the inertial to the viscous forces and can be written as

(1.2)\begin{equation} Re_A = \frac{\rho_g U_A d_{d}}{\mu_g}.\end{equation}

Here, $\mu _g$ is the dynamic viscosity of the local gas. The Ohnesorge number relates viscous forces to the inertial and surface-tension forces as

(1.3)\begin{equation} Oh = \frac{\mu_{d}}{\sqrt{\rho_{d} d_{d} \sigma}} ,\end{equation}

noting the impact and aerobreakup communities write it the same way and that, here, $\rho _d$ and $\mu _d$ refer to the liquid drop density and viscosity. The Mach number is velocity relative to the local sound speed as

(1.4)\begin{equation} M_A = \frac{U_A}{a_g} = \frac{U_A}{\sqrt{\gamma_g R_g T_g}},\end{equation}

where $a_g$, $\gamma _g$, $R_g$ and $T_g$ are the sound speed, ratio of specific heats, the gas constant and the temperature of the local gas, respectively.

The time, $t$, since interaction with the bow shock of the projectile is typically non-dimensionalized by the so-called ‘characteristic breakup time’ (Ranger & Nicholls 1969; Reinecke & McKay 1969) or ‘Rayleigh time’ (Hébert et al. 2020),

(1.5)\begin{equation} t_{C2} = \sqrt{\frac{\rho_d}{\rho_{g2}}} \frac{ d_{d} }{ U_{2LF} },\end{equation}

as $\tau = t/t_{C2}$. Here, the conditions immediately behind the bow shock of the projectile are chosen to define the characteristic time so that comparisons are consistent with the literature.

The first part of the problem of high-speed-vehicle/rain encounter is drop aerobreakup, which will be the focus of this paper. Work on the impact portion of this problem will be left to a follow-on paper. Aerobreakup is often characterized by different breakup regimes corresponding to increasing aerodynamic force on the liquid drop (Hinze 1949a,b, 1955; Lane 1951; Pilch & Erdman 1987; Hsiang & Faeth 1995). As the ratio of aerodynamic force to restoring capillary force on the drop increases (non-dimensionalized as the Weber number, $We_A$), the drop is said to pass through different breakup regimes. At low Weber numbers ($We_A < 100$), the drop is said to pass from a vibrational breakup regime through the bag breakup regime to the bag-and-stamen breakup regime and then the multimode breakup regime (Pilch & Erdman 1987).

Above $We_A \approx 350$, there exists controversy in the literature on what the breakup mechanism is. Additionally, there is uncertainty on what value of $We_A$ for which there exists, if any, a high-Weber-number regime, whereby the Weber number no longer delineates breakup mechanisms. Pilch & Erdman (1987) propose that, at higher Weber numbers, the drop enters a sheet-stripping regime ($100 \lessapprox We_{A} \lessapprox 350$), a wave-crest stripping regime ($We_{A} \gtrapprox 350$) and then a regime they term ‘catastrophic’ also for $We_{A} \gtrapprox 350$. The ‘catastrophic’ regime is characterized by relatively long-wavelength surface instabilities that grow in amplitude until they are comparable to the drop diameter, at which point the drop breaks into smaller droplets. This can occur before wave-crest stripping can completely reduce the drop mass.

Theofanous, Li & Dinh (2004), Theofanous et al. (2007, 2012) and Theofanous (2011) propose a different understanding of drop demise. In their model, Rayleigh–Taylor (RT) surface instability dominates initially at higher ($We_{A} \approx 1\times 10^2$) Weber numbers. In the RT-dominated regime (dubbed RT piercing, or RTP), surface waves are comparable in amplitude to the drop diameter. Starting from $We_A \approx 1\times 10^2\unicode{x2013}1\times 10^3$, and as $We_A$ becomes large, the regime transitions to shear-induced entrainment (SIE), where mass is stripped from the drop. Several researchers have performed research at Weber numbers in the regime where RTP is expected to transition to SIE. Wierzba & Takayama (1988) experimentally studied aerobreakup at $We_A = 600\unicode{x2013}7600$ with holography and presented a four-stage stripping-type mechanism. Chou, Hsiang & Faeth (1997) studied experimentally the temporal properties of drop breakup in the shear breakup for $We_A=125\unicode{x2013}375$. Dorschner et al. (2020) studied numerically and experimentally the formation and recurrent shedding of ligaments in drop aerobreakup in this transitional regime. The introduction of that paper gives a good review of the different breakup regimes and asserts that the ligament formation process in the vicinity of the RTP–SIE transition and beyond (i.e. $We_A >1\times 10^2$) is a subject of current investigation (Jalaal & Mehravaran 2014; Meng & Colonius 2015, 2018; Jain et al. 2019; Wang et al. 2020). Theofanous & Li (2008) used laser-induced fluorescence to image the aerobreakup of drops in Mach 3 flow at Weber numbers up to $5.4\times 10^3$. The liquid was tributylphosphate (TBP), which those researchers call ‘water-like’, given it is Newtonian, has a higher viscosity than water (4 mPa s) and lower surface tension (27.3 mN m$^{-1}$) at room temperature. Those researchers assert that previous research performed with shadowgraphy ‘allowed misinterpretations that have led to inappropriate conceptualizations (and theory) of the physics that govern breakup at high Weber numbers’. Continuing, ‘for $We_A>1\times 10^3$, breakup occurs by SIE, which is a regime dominated by shear (and damping of RT waves due to shear-induced straining) and one that is mechanistically different from (the previously hypothesized) stripping’. The researchers also discuss potential phase change and surface-tension reduction in high-temperature experiments. Theofanous et al. (2012) studied water drops at shock Mach numbers $M_s=1.1\unicode{x2013}2.1$ and $We_A=60-3\times 10^4$. They also studied TBP at $M_s=1.02\unicode{x2013}2.5$ and $We_A=12\unicode{x2013}1.5\times 10^5$. The results of that work lead the authors to assert that SIE is the ‘second and terminal criticality’. Theofanous (2011) provides a complete review of aerobreakup, discussing the relevant non-dimensional numbers pertinent to Newtonian and viscoelastic fluids. In that review, § 4.2 titled ‘Shear-Induced Entrainment is the Terminal Regime’ discusses how the ‘so-called catastrophic regime, thought to exist for Weber numbers greater than $\approx 1\times 10^3$, was an artifact of shadowgraphs’. They then suggest that, for increasing Weber number, SIE is increasingly favoured over RTP.

More closely related to this paper is research at high Weber ($We_{A} > 1\times 10^4$) and Mach number ($M_s>3$). Ranger & Nicholls (1969, 1972) studied drop demise in a shock tube at shock Mach number 1.5–3.5, in atmospheric air, yielding Weber numbers above $1\times 10^4$. Notably, those researchers devised the ‘characteristic breakup time’ and proposed a detailed model of boundary-layer stripping, noting that ‘[t]he shearing action exerted by the high-speed flow causes a boundary layer to be formed in the surface of the liquid and the stripping away of this layer accounts for the breakup’. This is likely among the first times SIE was recognized as a drop-demise mechanism.

Reinecke & McKay (1969) studied water drop breakup behind Mach 3–12 shocks at relatively high pressure, 90–760 torr, yielding results with Weber numbers greater than $1\times 10^4$. They found that ‘at shock speeds from Mach 3 to Mach 6 the drops are progressively stripped by the shearing of the surrounding air$\ldots$ Above Mach 6, where the drops experience accelerations in excess of $1\times 10^6$ g, an apparent surface instability results in the destruction of the drop before breakup is accomplished by the stripping action’.

Krauss (1970) studied water drop displacement and flattening behind Mach 1.6–3 shocks in a shock tube. A stream of liquid drops of diameter 0.5 to 3.2 mm was injected into a test section at 1 atmosphere, yielding aerobreakup Weber numbers in the range $1.3\times 10^3$ to $1.2\times 10^5$. Jaffe (1975) performed similar work, studying drop displacement and flattening behind Mach 7.5–13.9 shocks in the NASA EAST shock tube with initial pressure of 0.13–0.50 atm, yielding Weber numbers greater than $1\times 10^5$.

In their review of the literature, Pilch, Erdman & Reynolds (1981) referred to Ranger & Nicholls (1969) and Ranger & Nicholls (1972) and their analysis of boundary-layer stripping. He computes the time required for the boundary-layer stripping process to reduce the mass of a liquid drop to zero and compares that calculated time with the measured drop-demise time found by Engel (1958) and Ranger & Nicholls (1969). He concludes that the required time for drop demise due to boundary-layer stripping alone is up to an order of magnitude greater than the observed demise time. He concludes that boundary-layer stripping cannot then be the sole mechanism of drop aerobreakup, and that a breakup mechanism driven by surface instability must also be present.

Boiko, Papyrin & Poplavskii (1987) studied the interaction of drops at an initial pressure of 20–100 kPa with a shock wave at Mach numbers $M_s = 2\unicode{x2013}6$ and Weber numbers in the range $1\times 10^4\unicode{x2013}1.3\times 10^5$. Those researchers tracked the drop movement. They also assert that a ‘major role in disintegration of low viscosity liquids [referring to water] in a shock wave with $M_s = 2\unicode{x2013}4$ is played by the mechanism of removal of a surface layer of liquid, which acts after an induction period $t_i$ and continues through the complete decay of the droplet’.

Joseph, Belanger & Beavers (1999) studied drop demise in a shock tube at shock Mach number 2–4 with pre-shocked air pressures of 50–58 kPa, yielding Weber numbers in the range $4\times 10^4\unicode{x2013}1.6\times 10^5$ for drops of water in addition to several other liquids. ‘The thesis of [that] paper is that breakup at high accelerations, corresponding to high Weber numbers, is controlled at early times by Rayleigh–Taylor instabilities’. Joseph, Beavers & Funada (2002) continued that work for Newtonian and viscoelastic fluids. Joseph et al. (1999, 2002) both provide a detailed stability analysis to support their assertions based on experimental data.

Hébert et al. (2020) studied aerobreakup experimentally and numerically at Mach 4.4 and Weber numbers above $1\times 10^5$, which overlaps the parameter space of the work presented in this paper. In that work, they report on the drop displacement and flow morphology during three breakup steps delineated by non-dimensional breakup time. The drops were suspended at the crossing of two 25 $\mathrm {\mu }$m copper wires by capillary force. Virot et al. (2023) used a similar strategy to study aerobreakup at Mach numbers 4.3 and 10.6 for long non-dimensional breakup time at Weber numbers $5\times 10^4\unicode{x2013}1.1\times 10^5$. The authors make comments on heat transfer and the mist formation at the differing Mach numbers.

Salauddin et al. (2023) images the breakup of liquid rocket propellant drops at $M_s=5.1$ and $We_A=1.2\times 10^5$ in a detonation facility. They argue that the ‘transition from a shock wave to a detonation suppresses the deformation of the droplet and augments small-scale breakup’. They go on to discuss that, at high Weber and Mach numbers, ligaments are generated and finally, a ‘jellyfish’ structure is formed during the drop breakup process.

This review of the current literature identifies a dearth of high-quality aerobreakup data at high Mach, Weber and Reynolds numbers. Moreover, much of the research in the literature was performed with a shock tube, where there is a well-defined, impulsive start to the flow. In the present work, new experimental data and computations are presented on aerobreakup at high Mach and Weber numbers in the stagnation region of high-speed flow over a bluff body.

Additionally, a long-standing discrepancy in the literature is identified where some researchers assert that linear instability growth results in rapid and sudden drop breakup when the instability amplitude becomes comparable to the drop thickness (Reinecke & McKay 1969; Joseph et al. 1999). Others, including Theofanous & Li (2008) and Theofanous (2011), assert that this is incorrect and that drop breakup is governed exclusively by mass stripping. The present work develops a model for surface instabilities based on linear-stability theory from Funada & Joseph (2001) and shows a mismatch between the predicted amplitude and wavelength from that model and those observed in either shadowgraph images or numerical simulations. A model for nonlinear growth from Pilch et al. (1981) was used for large instability amplitudes and further developed; agreement between this model and computations is good.

2. Facility and experimental set-up

The experiments were performed at the Naval Surface Warfare Center Dahlgren Division's electromagnetic launcher (EML) located at the Hypersonic Integration and Test (HIT) facility at the Potomac River Test Range located in Dahlgren, Virginia, USA (Tadjdeh 2017). The EML digitally prescribes projectile launch acceleration profiles so it provides a ground-test capability with a broad range of precisely controlled velocities for weather encounters and terminal impacts with a launch cadence of multiple shots per day (O'Rourke 2017). For this testing campaign, eight bore-rider projectiles were fired at velocities between 1000 and 1800 m s$^{-1}$ into a sea-level atmosphere. The flat-faced, bore-rider projectiles were approximately rectangular and 100 mm wide $\times$ 150 mm tall. Water drops of diameter from 0.51–2.30 mm were suspended along the projectile's line of flight using acoustic levitators (Marzo, Barnes & Drinkwater 2017; Sheldon 2019), noting that the use of acoustic levitation in shock tubes has been done before by, for example, Hanson, Domich & Adams (1963) and Dorschner et al. (2020). Acoustic levitators were chosen for this application instead of suspending the drops from wire, as in Hébert et al. (2020), or releasing the drops from a dripper, as in many other works. This is because levitators made it easier to hold drops in the small field of view in the camera's relatively narrow depth of field. The levitators did not observably disturb the drop surfaces; pressure waves from levitator transmitters are not present in shadowgraphs and there is no observable distortion of the projectile bow shock or on the drop surface before the projectile arrives.

In figure 1, three frames of high-speed video looking uprange toward the rail-gun muzzle are presented along with a schematic of the experimental set-up. A coordinate system is placed in each to orient the reader. The projectile can be seen passing through the acoustic levitator frame before, during and after interaction with the drop. In figure 1-top, before the interaction, a green laser beam that acts as a ‘tripwire’ to initiate the high-speed cameras can be seen as a semitransparent green line. In figure 1-middle, during the interaction, a high-speed laser light source, seen as a semitransparent red line, is initiated and synchronized to a high-speed camera. Finally, in figure 1-bottom, after the interaction, the high-speed projectile is seen having passed cleanly through the levitator frame.

Figure 1. Three frames looking uprange toward the rail-gun muzzle. The approximately rectangular (100 mm wide $\times$ 150 mm tall) projectile is seen traversing the flight path through the levitator. The three frames depict periods before the imaging system is triggered by the ‘beam trip’, during shadowgraph recording, and after the projectile has traversed the field of view. Cartoons of drop and laser beams are not to scale. For additional context, a movie of this projectile flying through the levitator is provided as supplementary movie 7 available at https://doi.org/10.1017/jfm.2024.1092.

Back-lit shadowgraphy and front-lit imaging were used to image the drop (Hauser et al. 1936; Settles 2001; Danehy et al. 2018). For the back-lit shadowgraphy, the drop–projectile interaction is illuminated from behind by a laser light source synchronized to a Kirana high-speed camera operated at 2–4 million frames per second (f.p.s.) with 100 ns exposure. The lens used on the Kirana was a AF-S NIKKOR 200 mm f/2G ED-VR-II positioned approximately 2.8 m from the drop behind 50 mm of bulletproof glass. To achieve approximately 20 pixels mm$^{-1}$ resolution, one $2\times$ teleconverter, one $1.4\times$ teleconverter and 92 mm of lens tube were placed between the lens and the camera body. The experiment was designed such that an object moving at 500 m s$^{-1}$ would be imaged with one pixel of image blur. A bandpass filter at the laser wavelength was placed in the lens tube to attenuate muzzle flash to mitigate oversaturation. For the front-lit imaging, a Phantom TMX 7510 camera was run at 300k f.p.s., with a 190 ns exposure time; the light source was the rail-gun muzzle flash. A Sigma 150–600 mm f/5–6.3 was used with one $2\times$ and one $1.4\times$ teleconverter to provide 11.5 pixels mm$^{-1}$ resolution. The front-lit imaging did not yield useful data for this experimental campaign. The projectile interrupts a 532 nm laser beam, which is pointed at a DET36A2 biased photodetector with an FL532 bandpass filter. When that beam is interrupted, the falling edge signal from the photodetector triggers an oscilloscope's output, which triggers a pulse generator, in turn triggering the high-speed cameras.

The ambient and launch conditions of the six shots are reported in table 1. Here, $d_{d}$, $T_{1}$, $P_{1}$ and $\phi _{R}$ are the drop diameter, the ambient temperature and pressure and the relative humidity, respectively. Also, $U_{P}$, $M_{P}$, $\varDelta _{P}$, $E$ and $h_{T}$ are the projectile velocity, projectile Mach number, shock standoff distance, the programmed launch energy of the EML and total enthalpy of the flow, respectively.

Table 1. Dahlgren shot conditions.

3. Stagnation point flow model

In this section, a simple model is devised to find the flow parameters (e.g. Weber or Reynolds number) that define the problem of aerobreakup along the stagnation streamline of a bluff body in high-speed flow. This model enables the quick and clear definition of the problem and parameter space so that the relevancy of phenomena and appropriateness of assumptions may be assessed (e.g. calorically perfect gas, vibrational relaxation, thermochemistry, humidity, etc.). Here, it is assumed the drop is processed by the projectile bow shock and then by the gas along the stagnation streamline of the projectile; also, the drop does not disturb the mean flow. Further, a one-dimensional, reacting gas dynamics model, where the mass flux linearly decreases from the point behind the bow shock to zero at the corresponding stagnation point, is assumed to adequately capture the major characteristics of the flow, given the drops are much smaller than the projectile size ($>50\,:\,1$). The different states of the flow are marked in figure 2. This stagnation point problem is solved twice for each case, first for the projectile along $x_p$ from state 2 to state $sp$, and then for the drop itself along $x_d$ from state 3 to state $sd$. The calculations to predict the flow along the stagnation streamline were performed with Cantera (Goodwin 2003; Goodwin et al. 2022) using the shock and detonation toolbox (SDT) (Browne et al. 2018) in MATLAB. Technical details of the calculations and some nomenclature are in Appendix A.

Figure 2. Definition of the different states. Here, 1, 2 and $sp$ mark the ambient state (table 1), the post-bow shock state and the stagnation point on the projectile, respectively. The parameters $\varDelta _p$ and $x_p$ are the projectile bow shock standoff distance and the coordinate along the stagnation streamline from 2 to $sp$, respectively; 3 and $sd$ are the post-bow shock state of the drop and the stagnation state of the drop, respectively. Finally, $\varDelta _d$ and $x_d$ are the drop bow shock standoff distance and the coordinate along the stagnation streamline from 3 to $sd$, respectively.

Relevant post-shock quantities are non-dimensionalized and presented in figure 3. Here, a linearly decreasing mass flux is shown, per the model, along with a decreasing velocity in the projectile frame, $U_{PF}/U_2$. The gas velocity that the drop experiences is that in the laboratory frame; that is, $U_A$ from (1.1) is $U_{LF}/Up$. That the velocity experienced by the drop increases from the post-shock state to the velocity of the projectile is an illustration of the importance of considering the reference frame. The post-shock Mach number is supersonic in the laboratory frame, increasing from 1.9 to over 2.2. The thermodynamic parameters do not change much in the stagnation region of the projectile, namely the temperature ($<2$ %), density ($<8$ %) and ratio of specific heats ($<1$ %). This is in part because the enthalpy of these shots is not high enough to induce substantive chemical reaction. Using the definition of the Knudsen number from Loth (2008), it was determined that all experiments were executed in the continuum regime. Importantly, the ratio of specific heats is observed to be constant and $\gamma \approx 1.3$. The aerobreakup Weber number, $We_A$, and Reynolds number, $Re_A$, increase by ${\approx }50\,\%$ and ${\approx }25\,\%$, respectively. Although these increases are appreciable, they most likely do not alter the drop-breakup regime. Additionally, $We_A = \rho _g U_{LF}^{2} d_{d} / \sigma$, so it essentially parameterizes the aerodynamic pressure on the drop; that is, the aerodynamic pressure increases ${\approx }50\,\%$.

Figure 3. Non-dimensional quantities for the solution along the stagnation streamline for shot 3. Here, lines are from the thermicity model and the dots represent the simulations detailed in § 4.

A similar solution procedure is used to calculate the conditions behind the bow shock of the drop, a compression from state 3 to a hypothetical stagnation point on the drop, state $sp$ (figure 2). The initial conditions are derived from the jump conditions; however, state 1 is replaced with values computed along the stagnation streamline, and state 2 becomes state 3. The case where the drop is half-way along the stagnation streamline, $x_p = \varDelta _p/2$, is highlighted. The solution is presented in figure 4 with non-dimensional quantities and mole fractions. Some chemical non-equilibrium is apparent, as up to 4 % NO is being produced at a hypothetical stagnation point on the drop face. The ratio of specific heats is nearly constant, $\gamma \approx 1.28$.

Figure 4. Non-dimensional quantities (a) and mole fractions (b) along the stagnation streamline from behind the bow shock of the drop to the drop face. Here, lines are from the thermicity model and the dots represent the simulations detailed in § 4.

Dimensional and non-dimensional conditions at the post-shock states 2 and 3 are summarized in tables 2 and 3, respectively. Aerobreakup Weber/Reynolds numbers at state 3 are reported because these are the conditions to which the drop is exposed. This is analogous to boundary-layer studies on high-speed bodies referencing the Mach/Reynolds number at the boundary-layer edge, not the free-stream Mach/Reynolds number. Because the jump condition for mass dictates that the product of $\rho U$ is constant and the dynamic viscosity will go up due to higher temperature in state 3, $Re_{A2}$ is similar to $Re_{A3}$. The Weber numbers corresponding to state 3 are approximately 30 %–60 % of that in state 2. As such, it is asserted that the non-dimensional numbers at state 2 offer a sufficient point of comparison with other aerobreakup research at high Weber and Reynolds numbers because the values do not change enough as the drop traverses the stagnation streamline to change the aerobreakup regime. Additionally, the gas in the stagnation region of the projectile at these conditions can be approximated as a calorically perfect gas with a ratio of specific heats of $\gamma \approx 1.3$. This assumption was used to complete the detailed simulations in § 4.

Table 2. Dimensional conditions in states 1, 2 and 3.

Table 3. Non-dimensional parameters for the current test campaign.

The drop is subject to supersonic flow in the laboratory frame for all cases presented in this work. In figure 5, the shock Mach number, $M_s$, is plotted against the Mach number of the drop in the laboratory frame for a perfect gas and a gas in thermochemical equilibrium. For M$_s\approx 2$, the drop will be in transonic flow in the laboratory frame. Much research in this area is performed in facilities with M$_s\approx 2$. To highlight the effect on the flow field, (14) from Loth (2008) is used to plot the drag coefficient for a sphere using $M_{2LF}$ as the input. All of the conditions in this work are at projectile Mach numbers of at least $M_P=M_s=3$, so the Mach number in the laboratory frame for this work is not transonic, as $M_{2LF}>1.4$.

Figure 5. Shock number, $M_s$, vs post-shock Mach number in the laboratory frame, $M_{2LF}$, for a perfect gas with $\gamma =1.4$ and in chemical equilibrium. Ratio of specific heats as computed by Cantera, $\gamma$, and $C_D$ as computed with (14) of Loth (2008) for $M_{2LF}$.

The calorically perfect gas assumption is sometimes not appropriately applied in the high-speed aerobreakup literature. The calculations done by Theofanous & Li (2008) in their tables I and II while discussing the work of Reinecke & McKay (1969) do not account for thermochemistry. Those calculations are performed assuming a calorically perfect gas with the ratio of specific heats set to 1.4 and predict temperatures that are too high, and the resulting assertions that are made may not be correct. Accounting for thermochemistry in those calculations would reduce the post-shock temperature considerably (Thompson 1988; Anderson 2011). It is imperative that the often-used calorically perfect gas assumption with the ratio of specific heats set to 1.4 be evaluated for shock Mach numbers exceeding 3, with or without a projectile. The post-shock thermodynamic parameters ($P, T, \rho$) can be incorrectly predicted, which can change basic flow features, such as the shock standoff distance of the drop.

4. Numerical simulation set-up

This computational study utilizes the Allaire, Clerc & Kokh (2002) viscous five-equation model, which encompasses two equations dedicated to mass conservation of the liquid and gas phases, momentum equations for each spatial direction, an energy equation and a non-conservative volume-fraction advection equation. The complete governing equations can be found in Appendix B.

The Allaire model has been used in previous aerobreakup and shock–drop interaction works (Meng & Colonius 2018) due to its ability to suppress pressure oscillations at material interfaces (Abgrall 1996). The model was designed to characterize immiscible fluids, and so, the interfacial region is subjected to artificial mixing due to resolution constraints and missing surface tension. Nevertheless, comparison with experiments, at least for the high-Mach-number conditions considered in this work, seems to build confidence in the ability of the model to capture the major instability mechanisms in the front face, as well as the evolution of the mist boundaries (see figure 18 in § 8.3).

The governing equations are solved with a fifth-order weighted compact nonlinear scheme adapted from Wong et al. (2021) which employs higher-order flux and variable reconstruction, as well as a positivity-preserving limiting procedure switching to a first-order Harten–Lax–van Leer contact Riemann solver to ensure robustness. The viscous terms are discretized using regular centred second-order accurate conservative finite difference discretization. For time integration, a third-order strong stability-preserving Runge–Kutta method is employed. The numerical simulation approach was previously validated in several works (Viqueira-Moreira et al. 2022; Stoffel et al. 2023; Viqueira-Moreira et al. 2023a,b).

For the current simulations, the projectile with a rectangular cross-section is immersed in the flow domain using an established immersed boundary method (Brehm, Hader & Fasel 2015; Brehm, Barad & Kiris 2019; Boustani et al. 2021). Firstly, a perfect gas base flow is computed to establish a flow field around the projectile. Secondly, a drop is inserted in front of the shock front convecting with the mean flow. For computational efficiency, adaptive mesh refinement (AMR) is used, tracking the pressure gradient to follow the shock front and the water volume fraction to track the drop. A total of 12 refinement levels are used to be able to provide sufficient grid resolution around the drop to maintain at least 150 grid points per diameter. A schematic of the simulation set-up is shown in figure 6(a), the flow field around the projectile in figure 6(b) and an overview of the AMR grid topology with colour contours of $u$-velocity contours in figure 6(c).

Figure 6. Schematic showing adaptive mesh refinement tracking of the drop.

5. Overview of results and discussion

5.1. Overview of experimental results

The high-speed movies from each of the cases are included as supplementary material and are numbered according to the shot they depict. Supplementary movie 1 shows the first successful experiment, illustrating that aerobreakup and impact can be visualized in the unique environment of an outdoor rail-gun facility. Supplementary movie 2 shows two drops in the same stagnation region. Here, the importance of the characteristic breakup time, $t_{C2}$, is seen as the smaller of the two drops breaks up appreciably earlier and does not impact the projectile. Supplementary movies 3–5 were recorded with nominally matching conditions, enabling the study of different sized drops (thus, different Weber and Reynolds numbers) at projectile Mach number $M_P=5 +0.12/-0.08$. Supplementary movie 6 was recorded with $M_P=3.03$ to overlap with previous shock tube aerobreakup studies and complete the parameter space as best possible at the time. Two drops are seen in the stagnation region, the smaller of which appears to break up before impact. Finally, video of the projectile flying through the levitator is provided as supplementary movie 7.

Hereafter, the results from the drop 3 experiment (supplementary movie 3) and the drop 4 experiment (supplementary movie 4), which are at nominally matching conditions, are compared with the results from the drop 6B experiment (supplementary movie 6). These cases are compared because they highlight distinct differences in the flow topology observed at noticeably different Mach and Weber numbers. The naming convention of each drop follows tables 1 and 3. Four sample shadowgraph images from the drop 3 experiment are shown in figure 7 to give the reader an idea about the finer flow features that could be captured in the experiments.

Figure 7. Four shadowgraphs of the high-speed movie for drop 3 showing the shock structure about the drop. Short red ticks denote millimetre increments. Long red ticks denote centimetre increments. Full movie included as supplementary movie 3.

5.2. Overview of computational results and discussion

Figure 8 gives an overview of the aerobreakup process as observed in the simulations. Drops 3 (Mach 5 case) and 6B (Mach 3 case) display similar overall dynamics, but some key differences are observed and will be analysed hereafter. As the drop enters the shock layer (shock structure is discussed in detail in § 6), the aerobreakup process begins, causing a substantial change to the drop shape. Different colour schemes were chosen for velocity contours in figure 8 for inside, outside and the wake of the drop, as the velocity scales in these regions of the flow domain are widely different. Firstly, the drop experiences a strong, non-spatially uniform deceleration due to compressibility effects. The front of the drop experiences much larger acceleration than the back, causing flattening. The flattening of the back of the drop shape has been discussed in prior work. For example, Meng & Colonius (2018) note that high pressures at the forward and rear stagnation points contribute to the flattening. In the present simulations, the pressure at the front face is the main driver of this ‘pancaking’, with much lower pressures at the rear. Negative pressures can be observed inside the drop, meaning it undergoes internal tension (Azouzi et al. 2012; Caupin et al. 2012). This is modelled through the stiffened gas equation of state (see Appendix B). Although cavitation effects could play a role in the flow field evolution, they are disregarded in this study. Additionally, the computations predict that the aspect ratio of the liquid mass of the drop changes appreciably with time. That is, mass is shed from the windward drop face while the aerodynamic forces flatten the drop, accentuating this increase in aspect ratio. Figure 8 shows that the liquid drop has an aspect ratio of $b/c\approx 1, 2, 4$ (spanwise vertical, $y$ to streamwise, $z$) at times $\tau \approx 0.25, 0.50, 0.75$, respectively. Drop displacement and flattening are discussed in detail in § 7. The wake dynamics is distinct in figure 8, with the wake of drop 6B becoming turbulent sooner than that of drop 3. Instability at the windward gas/liquid interface is readily apparent in figure 8. The presence of the RT and Kelvin–Helmholtz (KH) instabilities can be seen in the 0.1-volume-fraction contour lines where RT seems to dominate the front centre region of the drop, and the KH instability becomes more apparent closer towards the equator. Once the drops have substantially flattened, strong bursting events can be observed, e.g. for drop 3 at $\tau =1.0$. Both drops are close to catastrophic breakup at $\tau =1.0$, displaying drop surface waviness amplitudes reaching the thickness of the drop. With increasing simulation time, the dominant wavelength of the deformed gas–liquid interface is increasing; e.g. compare snapshots $\tau =0.25$ and $\tau =1.0$. These characteristics will be investigated further with linear- and nonlinear-stability analyses in § 8.

Figure 8. Velocity and pressure contours in centre cut plane for drops 3 (a) and 6B (b) at four different times. Note different colour maps correspond to different ranges for inside the drop, outside the drop and in the wake.

Before discussing further details of the aerobreakup process observed in the simulations, it should be noted that the simulations neglect some of the multiphase physics, such as surface tension and the possible formation of a layer of mist or dispersed phase surrounding the drop which may alter the simulated breakup process. Furthermore, the grid resolution with 150 grid points per drop diameter is not sufficient to accurately resolve the most amplified linear-instability mode with an estimated wavelength of ${\approx }30\ \mathrm {\mu }$m, as will be discussed in § 8.2. The good agreement between the simulations and the experiments suggests, however, that the dominant physical mechanisms are properly captured.

6. Shock structure

In figure 7, the projectile is coming from the left and the projectile bow shock is seen to process the drop. Then, the shock structure develops around the drop. In contrast to a solid sphere, there is no separation shock at the drop equator due to the change in the boundary condition. The wake shock does appear and its two termini are the base of the turbulent multiphase wake, and the bow shock; this bell shape of the wake shock is consistent among all cases.

For both the experiments and simulations, time is zero when the bow shock first intercepts the drop. In figures 9 and 10, the experimental snapshots for drops 3 and 6B are compared with the simulations, respectively. The angle of the main shock is steeper in the free stream for drop 3 in comparison with drop 6B due to the larger Mach number.

Figure 9. Comparison of shadowgraphs (above) vs numerical schlieren (below) for drop 3. Length scale represents the initial drop diameter.

Figure 10. Comparison of shadowgraphs (above) vs numerical schlieren (below) for drop 6B. Length scale represents the initial drop diameter.

The simulation results are in good agreement with the experiments and properly capture the position and shape of the main shock in front of the drop. The recompression shock is significantly harder to model as it is affected by some of the wake dynamics behind the drop. Overall, good agreement in shock position and shape has been obtained, providing confidence in the simulation results, and the dominant features can be captured in the simulations.

The initial transient in the shock dynamics around the drop is similar to what is observed for a rigid sphere, except for some behaviour in the wake region. Numerical schlieren images with highlighted flow features are provided in figure 11. The simulation results are qualitatively similar to what has been observed for the shock–cylinder interaction problems in the literature (Bryson & Gross 1961). First, the shock intercepts the front face of the drop and is reflected; this shock forms the bow shock over the drop. Following the terminology in Bryson & Gross (1961), a curved ‘Mach shock’ and contact discontinuity form, establishing a triple point on the incident shock. Note this triple point is three-dimensional and revolves around the centre axis. Again, similar to Bryson & Gross (1961), a second ‘Mach shock’ and contact discontinuity form, coming to a triple point on the initial ‘Mach shock’.

Figure 11. Simulation results at drop 6B conditions. Comparison of the initial transient for shock–drop (top) and shock–sphere interactions (bottom). Ref. S = reflected shock, Rec. S. = recompression shock, S.S. = separation shock, C.D. 1 = contact discontinuity, M.S. = Mach shock, T.P. = triple point, V.P. = vortical protuberance. Terminology from Bryson & Gross (1961).

The separation shock does not appear on the drop; this is because of the partial-slip boundary condition on the drop surface. There is an additional feature that is not observed in shock–cylinder interactions: a protuberance that arises on the back of the drop and initially has a smaller radius. This feature comprises vortical structures, which rotate in the opposite direction to the upcoming first Mach shock that generates them after colliding in the wake region.

7. Drop displacement and flattening, and comparison with existing models

Several researchers have characterized the displacement and flattening of drops in high-speed flow (Engel 1955, 1958; Reinecke & McKay 1969; Hébert et al. 2020; Wang et al. 2020). Drop deformation is primarily driven by the difference in pressure between the gas at the stagnation point on the windward face of the drop and the decreasing pressure of the gas as it traverses the drop surface to the equator (Engel 1958). In this section, a comparison between experimental observations and computations is made for drop diameter and flattening; these data are compared with existing predictive models in the literature. Moreover, the computations are analysed to investigate the drop deformation mechanics using cut planes showing pressure and velocity contours.

As in Hébert et al. (2020), the displacement of the mist leading edge (MLE) from its initial position was used as a proxy for the displacement of the drop centre of mass, which cannot be seen in shadowgraphs. The complete data set and the process by which the location of the MLE was determined are discussed in Appendix C. Those data indicate an average experimental value and standard deviation for a drag coefficient of $C_D=2.32\pm 0.76$. This is nominally consistent with Hébert et al. (2020) ($C_D=2.15$) and Pilch et al. (1981) ($C_D=2.52$), bringing confidence to the current results. A subset of experimental and computational displacement data for drops 3, 4 and 6B are presented, along with error bars, in figure 12(a). Drops 3 and 4 were chosen because those data were taken at nominally the same Mach 5 condition, so repeatability could be assessed. Drop 6B was a Mach 3 shot, which represents the largest available contrast in Mach/Weber number in this dataset (see table 3). The computational predictions fall within the error bars for drop MLE displacement.

Figure 12. (a) Non-dimensional displacement of MLE and (b) non-dimensional diameter for drops 3, 4 and 6B, comparing simulation output with experimental data.

Experimental and computational drop flattening data for drops 3, 4 and 6B are compared in figure 12(b) with three flattening models found in the literature. The Reinecke & Waldman (1975) model is a piecewise function where drop diameter is a function of the maximum pressure coefficient on the drop. The TAB model from O'Rourke & Amsden (1987) approximates the water drop as a constant-mass spring–damper system, where the restoring force is provided by drop surface tension while the drop viscosity dissipates energy. Inputs to this model are the Weber and Ohnesorge numbers, along with a constant that we have set as $C_2=2$, which is nearly the same value for the constant as in Daniel et al. (2023). The nonlinear TAB (NTAB) model from Schmehl (2002) introduces additional terms to account for the deformation of the liquid drop toward an ellipsoidal shape. Inputs to this model are the Weber and Reynolds numbers, and a constant which, as in the TAB model, we have set as $C_2 = 2$. Both the TAB and NTAB models include $1/We_{A2}$ or $1/Re_{A2}$ terms, which decrease to very small values in this study; consequently, the growth rates predicted by each model for drops 3, 4 and 6B are almost indistinguishable when plotted, and only the model results for drop 6B are plotted in figure 12(b). Drops 3 and 4 exhibit similar flattening trends, but differ from one another enough that their error bars barely overlap. This indicates that there is a statistical spread in the flattening data. The root-mean-square difference between the models and experiment is approximately 5 %–10 %. A so-called ‘burst’ is highlighted in figure 12(b) as the drop diameter temporarily changed substantially when a chunk of the drop was cleaved off and entrained into the wake.

In contrast to the Reineke correlation, and the TAB/NTAB deformation models, Su et al. (1996) present a drop-breakup model, the basis of which is derived from KH and RT linear-stability analysis. This KH–RT hybrid model was designed to analyse the breakup of liquid sprays at velocities of ${\approx }100\ {\rm m}\ {\rm s}^{-1}$ in a near-standard-temperature-and-pressure environment with $We_A<1000$. The times calculated for liquid drop breakup for cases 1–6B by this hybrid model are ${\approx }1\ \mathrm {\mu }$s for the RT mode and ${\approx }500\ \mathrm {\mu }$s for the KH mode. These calculated breakup times differ greatly from the observed drop-breakup times in this work and from the Rayleigh breakup time $t_{C2}$ computed in prior high-speed aerobreakup research. It is determined that the hybrid KH–RT model presented by Su et al. (1996) is not applicable to these flow conditions. This motivates further discussion of the linear- and nonlinear-instability regimes, as is done in § 8.

8. Instability mechanisms and models

Linear- and nonlinear-stability analyses are used to attempt to explain the developing structure at the liquid–air interface with the ultimate aim of eventually predicting drop demise. Computational results at early breakup times motivate a linear-stability analysis. A disconnect between the structures observed in both the computations and experiments and those predicted by the linear-stability analysis motivates a nonlinear-stability analysis. A comparison of the present work with that in the literature is made.

8.1. Computational results for early breakup time

The simulated drop exhibits ripples across its surface right after its interaction with the projectile bow shock (figure 13). The drop starts forming distinct annular structures and those structures keep growing throughout the test time. The shape of the back of the drop takes the characteristic ‘cupcake’ shape seen in lower-Mach-number experiments and simulations (Meng & Colonius 2018), but the front face is different, as it is not smooth, which is the result of higher shear in the present experiments.

Figure 13. Three-dimensional water volume fraction contours of 0.1 for drops 3 and 6B at early stages.

8.2. Simple linear-stability model for early breakup time

A simple treatment of the stability problem at the drop/gas interface is presented in this section to estimate the importance of shear, acceleration, viscosity, surface tension and geometry (aspect ratio) for early breakup times. The drop is sketched in figure 14(a). It is assumed that the front face takes the form of an ellipse as it deforms. Figure 14(b) shows a schematic of the classical stability treatment of an accelerated interface from Funada & Joseph (2001) of the Kelvin (1871) and Helmholtz (1868) instability. Here, it is assumed that the initial stability problem can be treated as a locally accelerated interface, $\zeta (x',t) = \zeta _0 \exp ( {\rm i} \kappa x' + s t )$, and the stability analysis is locally valid as a function of $\phi$ from the stagnation point of the drop. Funada & Joseph (2001) show, just after their (2.18), that the growth rate is

(8.1)\begin{align} s &={-} {\rm i} \kappa \frac{\rho_g u_g + \rho_d u_d} { \rho_g + \rho_d } \underbrace{- \frac{\kappa^2 (\mu_g + \mu_d)}{\rho_g + \rho_d } }_{Viscosity} \pm\left[ \underbrace{\frac{\kappa^2 \rho_g \rho_d \left(u_g - u_d \right)^2 }{(\rho_g + \rho_d)^2}}_{Shear} \underbrace{- \frac{\kappa a_{{\perp}} (\rho_g - \rho_d)}{\rho_g + \rho_d}}_{Acceleration} \right. \nonumber\\ & \quad \left. \underbrace{- \frac{\kappa^3 \sigma}{\rho_g + \rho_d} }_{{Surface\ tension}} \underbrace{+ \frac{\kappa^4 (\mu_g + \mu_d)^2}{(\rho_g + \rho_d)^2} + 2{\rm i} \kappa^3 \frac{(\rho_g \mu_d - \rho_d \mu_g) (u_g-u_d)}{(\rho_g + \rho_d)^2} }_{Viscosity} \right]^{1/2}. \end{align}

Here, four combinations of terms in (8.1) that might drive/damp instability growth are labelled as shear, acceleration, viscosity and surface tension. These terms are estimated as a function of $\phi$ along the front face of the drop for various flow conditions and drop geometries. This was done by considering the real parts of the growth rate, ${\rm Re} (s)$, for each term while setting the others to zero.

Figure 14. (a) Schematic of deformed drop, exaggerated for sketch. (b) Adaptation of treatment by Drazin & Reid (2004) and Funada & Joseph (2001).

First, $a_\perp = a_\perp (\phi )$ must be found, which is the acceleration normal to the interface to which a non-inertial reference frame ($x'-z'$) is fixed. The acceleration is assumed to be constant, which is supported by our data in figure 12 and data in Hébert et al. (2020). The acceleration parallel to the interface, $a_\parallel$, is assumed to negligibly affect the stability problem given that $\zeta$ is small and in the regime treated by small disturbances. That is, the component pertinent to the stability problem is only $a_\perp$. The front face of the drop is estimated to take the form of an ellipse as $F(x,y,z) = 0 = (z/c)^2 + (y/b)^2 - 1$. Then, the unit normal, $\underline {\hat {n}} = \boldsymbol {\nabla } F / |\boldsymbol {\nabla } F|$, can be used to calculate the needed quantities to approximate the acceleration as a function of $\phi$, $\tan (\phi _n) = (c/b)^2 \tan (\phi )$ and $a_\perp = |\underline {a}| \cos (\phi _n)$. The magnitude of acceleration is the drag force divided by the drop mass:

(8.2)\begin{equation} |\underline{a}| = \frac{F_d}{m_d} = \frac{3 \rho_2 C_d U_{2LF}^2}{4 \rho_d d} = \frac{3 C_d q_2}{2 \rho_d d}, \end{equation}

where $q_2=1/2 \rho_{2} U_{2LF}^2$ is the dynamic pressure based on state 2 in the laboratory frame, and $C_d$ is the drag coefficient of the drop which is estimated to be a constant $C_d\approx 2.32$ (figure 12).

Next, $u_g(\phi )$ and $\rho _g(\phi )$ must be found. A modified version of Newton's inclination method is used to find the pressure coefficient as

(8.3)\begin{equation} C_p(\phi) = \frac{p_g(\phi)-p_2}{\dfrac{1}{2} \rho_2 U_{2LF}^2} = \frac{p(\phi)-p_2}{q_2} = C_{p\,{max}}\cos^2(\phi_n).\end{equation}

Here, it is noted that the modified form of Newton's inclination method was taken from Lees (1956), who used data for the front face of a hemispherical cylinder from Stine & Wanlass (1954) to support the simple model at Mach and Reynolds numbers relevant to this work. The maximum pressure coefficient is calculated from the stagnation pressure at the drop face calculated in § 3 as $C_{p\,{max}} = (p_{sd} - p_2)/q_2$. Then, the Mach number, density and temperature of the gas at the interface may be calculated as a function of $\phi$ as

(8.4)\begin{equation} \frac{p_{sd}}{p_g(\phi)} = \left(\frac{T_{sd}}{T_g(\phi)} \right)^{{\gamma}/({\gamma-1})}= \left(\frac{\rho_{sd}}{\rho_g(\phi)} \right)^{\gamma} = \left (1 + \frac{\gamma-1}{2} M_g(\phi)^2 \right)^{{\gamma}/({\gamma-1})},\end{equation}

which yields the velocity, $u_g(\phi ) = M_g(\phi ) \sqrt {\gamma p_g(\phi )/\rho _g(\phi )}$.

The terms of (8.1) can now be estimated for an example condition, drop 3. This is done to understand which of shear, acceleration, viscosity and surface-tension terms are relevant on the drop face. In figure 15(a,c,e), the flow parameters along the front face of the drop are shown as a function of $\phi$. The solid lines are from Newton's inclination method, calculated by (8.3) and (8.4) and the dots are CHAMPS-CFD (Cartesian higher-order adaptive multiphysics solver; computational fluid dynamics) data for inviscid flow past a rigid sphere. Agreement between the inclination method and the CFD data is sufficient to justify using the former to estimate the four relevant terms in the growth rate ((8.1)). The shear (green), acceleration (red), viscosity (magenta) and surface-tension (blue) terms are shown as semitransparent surfaces in figure 15(b,d,f). Here, the negative of the surface-tension term is plotted so it may be presented in log space (the surface-tension term will always be negative and serve to dampen the instability). Also inset in black is the drop profile and aspect ratio ($b/c$). Also, the growth rate is normalized by the time scale $t_{C2}$ from (1.5).

Figure 15. Estimated flow parameters (a,c,e) and growth rate (b,d,f) for drop aspect ratios, $b/c=(1,2,4)$, respectively, for drop 3. The progression of $ab/cd/ef$ shows the effects of a drop becoming more elliptical (deforming). The lines and dots in (a,c,e) are the result of Newton's inclination method and inviscid CFD simulations for a solid sphere, respectively. The different surfaces in (b,d,f) correspond to the shear, acceleration, viscosity and surface-tension terms in (8.1). Inset in (a,c,e) is a trace of the drop shape in black.

Upon inspecting figure 15, it is apparent that viscosity and surface tension are not relevant in this simplified analysis, particularly at larger disturbance wavelength. More interestingly, as the aspect ratio, $b/c$, is increased, the terms that dominate the growth in (8.1) change. For a spherical drop, the shear term appears to dominate everywhere except near the stagnation point, where the acceleration term is most prominent. As the drop aspect ratio is increased from $b/c=1$ to $b/c=2$ and then an exaggerated $b/c=4$, the acceleration term appears to become more dominant. This is because the profile of the front face of the drop is becoming flatter; that is, the acceleration and surface-tension terms nominally stay the same but the shear term is largely reduced because $u_g(\phi )$ is smaller. This is the primary takeaway from this simplified analysis. This effect would be further exaggerated if $u_d(\phi )$ were permitted to vary from zero. The effect of $u_d(\phi )$ being non-zero is not dominant, as estimates of growth rate assuming $u_d(\phi ) \approx F u_g(\phi )$ with $F=0.25,0.5$ were performed and the results in figure 15 are not changed substantively; that is, increasing the drop aspect ratio is far more effective in reducing $u_g(\phi )$.

An additional takeaway from figure 15 is the broadband and complex nature of the most-amplified instability at early times. A range of wavelengths are amplified and the spectrum changes as the drop deforms. This analysis is a clear example of the need for high-quality images of the front face of the drop similar to those in Theofanous (2011). This simple analysis may also explain why the drop face in that work is smooth.

Treating the stability problem as a linearly perturbed interface that grows exponentially presents inconsistencies with experimental results for late breakup times. For the conditions encountered in this experimental campaign, the linear growth rate is given by (8.1) and is maximized for wavelengths ($\lambda \approx 30\ \mathrm {\mu }$m) smaller than those observed in shadowgraph images; that is, the resolution of those images (${\approx }21.5\ {\rm pixels}\ {\rm mm}^{-1} \approx 46.5\ \mathrm {\mu }{\rm m}\ {\rm pixel}^{-1}$), was not high enough to discern such small waves from image noise. Evaluating (8.1) at the stagnation point for drop 4 yields the linear amplification curve in figure 16(a). A shadowgraph image showing the incipient surface sawtooth profile of drop 4 is presented in figure 16(b). The growth rate $s$ is maximized for a surface instability of $\lambda \approx 30\ \mathrm {\mu }$m, yet the wavelength of the incipient sawtooth profile is ${\approx }0.5\unicode{x2013}1.0$ mm. This observation indicates a shortcoming in the linear-stability treatment and is consistent with observations in Reinecke & McKay (1969), Pilch & Erdman (1987), Joseph et al. (1999) and Theofanous (2011), among others.

Figure 16. (a) Linear growth rate ($s$) vs wavelength ($\lambda$) per (8.1) for drop 4 at the stagnation point showing maximum exponential growth rate at $\lambda \approx 30\ \mathrm {\mu }$m. (b) Incipient sawtooth profile of drop 4; average wavelength between ripple peaks ${\approx }\,990\ \mathrm {\mu }$m.

8.3. Experimental and computational results for late breakup time and comparison with the literature

The spatial resolution of the simulations is sufficient to capture concentric rings on the drop windward face, which the authors believe to be physical (figure 17). The exact structure, number of rings and bursts may depend on the mesh and numerical noise seeding initial instabilities; this is because of the Cartesian mesh imprint. However, a simulation sensitivity analysis was performed by increasing spatial resolution and that effort showed the main flow features are not substantively changed with increasing resolution.

Figure 17. Schematic showing the evolution of ‘bursts’ and front face instabilities for drop 3 for late breakup times.

In figure 18, numerical and experimental shadowgraphs are compared for drop 6B at later non-dimensional times. The computations are performed with no surface tension. Given the apparent qualitative agreement between experiment and computation, this brings confidence to the assumption made in the computations, at least for macroscopic characterization. The boundary layer cannot be resolved as a consequence of resolution constraints, as well as any microscopic drop generation due to smaller-scale instabilities. The serrated or ‘sawtooth’ front face seen at $\tau =1.17$, and already reported in the literature (Reinecke & McKay 1969) is present in both the experiments and simulations.

Figure 18. Comparison of experimental shadowgraphs (above) vs numerical shadowgraphs (below) for drop 6B at three different times.

In figure 19, the results from the present study and Hébert et al. (2020) are compared at similar conditions and times. The shock and wake structures appear to be similar, bringing confidence to the current results. More interestingly, the front surface of the drop appears to show a large-scale incipient ‘sawtooth’ or ‘fingering’ profile. Figure 20 shows results from Reinecke & McKay (1969) (top left), Hébert et al. (2020) (top right) and results from the present work (bottom). The results from the present work are for different $\tau$ but relatively similar $We_{A2}$. The non-dimensional time could potentially be explained by differing initial conditions in the three experiments. The different methods of suspending the drop could induce different initial conditions to the stability problem. However, the flow topology seems similar near the end of the breakup process. These fingering structures are similar to what is seen on the front face of the computations at later times (figure 17). Moreover, Pilch et al. (1981) and Pilch & Erdman (1987) suggest something similar in their figure 1 for ‘catastrophic breakup’. From Pilch & Erdman (1987), ‘[l]arge-amplitude, long-wavelength waves ultimately penetrate the drop creating several large fragments before wave crest stripping can significantly reduce the drop mass’. This implies a race between SIE and the larger, long-wavelength waves that penetrate the drop suggested by Joseph et al. (1999) among many others. That is, we are implying that the results from Joseph et al. (1999) and Theofanous & Li (2008) are not in conflict, as is sometimes understood, but the results are complementary. The simple linear-stability analysis done in § 8.2, and the analysis done in Joseph et al. (1999) is probably not appropriate for long non-dimensional times when these ‘fingering’ structures appear. It is most likely that a nonlinear analysis would be required. Not imaging the front face of the drop at late breakup times is a shortcoming of the experimental design and must be remedied in future work.

Figure 19. Comparison of Hébert et al. (2020) drops (above) with drops 5B and 3, noting they are recorded at similar non-dimensional (and dimensional) times and similar non-dimensional (and dimensional) conditions.

Figure 20. Comparison of Reinecke & McKay (1969) drop (upper left) with drop 6B (lower left), showing similar sawtooth profile on windward side. The Hébert et al. (2020) drop (upper right) shows incipient sawtooth formation, also visible on drop 4 (lower right); given more time, both of these drops would (most likely) develop a more pronounced sawtooth profile.

8.4. Proposed explanation of large-scale structures at late breakup times

Taking the maximum growth rate from figure 16, and assuming an initial amplitude of $\zeta _{0}=1\ \mathrm {\mu }$m, the instability is expected to grow to a size comparable to that of the drop itself in just a few microseconds; yet, no drop studied in this experimental campaign breaks up that quickly. As such, and consistent with the approach in Pilch et al. (1981) and Pilch & Erdman (1987), a ‘large-amplitude’ treatment of the stability problem is proposed to account for the fundamental discrepancy in the size of the observed structure at late breakup times.

Pilch et al. (1981), citing Davies & Taylor (1950), Garabedian (1957) and Fishburn (1974), used a model where the growth rate became a function of drop acceleration and instability wavelength after a certain instability size threshold was reached. According to this model, the growth of a surface instability on the windward face of a drop is equivalent to the penetration of liquid by an open bubble of gas, which is a common assumption in the treatment of nonlinear RT instability (Kull 1991; Zhou et al. 2021). The growth rate is constant in time and is written as

(8.5)\begin{equation} \frac{{\rm d}\zeta}{{\rm d} t}_{nonlin} = F\sqrt{(1-\epsilon )a\lambda},\end{equation}

where $a$ is the drop acceleration, $\epsilon$ is the density ratio $\rho _{g}/\rho _{d}$ and $\lambda$ is the wavelength of a surface instability. Here, $F$ is given by

(8.6)\begin{equation} F = Fr\sqrt{\left(1-4 \frac{{\rm \pi}^{2} \sigma}{(1-\epsilon)\rho_{d}a\lambda^{2}}\right)},\end{equation}

which reduces to the Froude number ($Fr$) in the case of zero surface tension. The Froude number was numerically calculated by, among others, Birkhoff (1955), Garabedian (1957) and Kull (1986) to be $\approx 0.23$ for a two-dimensional planar bubble and $\approx 0.36$ for a three-dimensional bubble. Equation (8.5) gives an upper bound for the velocity with which an open bubble penetrates a volume of liquid. Because this velocity constitutes an upper bound for the gas penetration velocity into a volume of liquid, Pilch et al. (1981) defines a transition time, $t_{t}$, between the linear (exponential) growth and the nonlinear (constant rate) growth as the time when the growth rates are equal:

(8.7)\begin{equation} \left. \frac{{\rm d}\zeta}{{\rm d} t}_{lin} \right|_{t_t} = \zeta_0 s \exp( {\rm i} \kappa x' + s t_t ) = \left. \frac{{\rm d} \zeta}{{\rm d} t}_{nonlin} \right|_{t_t} = F\sqrt{(1-\epsilon )a\lambda}.\end{equation}

Considering only the real parts of the growth rate, ignoring the effects of viscosity and restricting the analysis to the windward face of the drop, where $u_{g} \approx 0$, this relation can be solved for the transition time as

(8.8)\begin{equation} t_{t} = \frac{\ln\left( {\dfrac{F\sqrt{(1-\epsilon)a\lambda}}{s \zeta_{0}}} \right) }{s}.\end{equation}

The transition time between the linear- and nonlinear-instability growth regimes is a function of drop acceleration, $a$, surface instability wavelength, $\lambda$, the initial disturbance magnitude, $\zeta _{0}$, and the exponential instability growth rate, $s$. Pilch et al. (1981) states that there must exist a smooth transition between the exponential and constant-rate instability growth regimes. For simplicity, the surface instability growth is instead here modelled as a piecewise function. For the conditions of drop 4, the transition time $t_{t}$ was evaluated for a range of different initial disturbance magnitudes and disturbance wavelengths. Figure 21(a) shows that, for plausible initial disturbance magnitudes of the order of 1 $\mathrm {\mu }$m and surface instability wavelengths of the order of those seen in shadowgraph images, surface instability growth transitions to the nonlinear regime well before $\tau = 1$ is reached. That is, for much of the drop-breakup process, surface instability growth is nonlinear and can be modelled as in (8.5) rather than the exponential model that applies for small instability magnitudes.

Figure 21. (a) Non-dimensionalized transition time $\tau _{t} = t_{t}/t_{c2}$ for drop 4 for a range of initial disturbance magnitudes $\zeta _{0}$ and instability wavelengths $\lambda$. (b) Penetration rates for drop 4 calculated for different $\lambda$ and $\zeta _{0}$. (c) Surface instability amplitudes as a function of time for drop 4 computed for different $\lambda$ and $\zeta _0$ according to (8.9).

The instability amplitude $\zeta _{f}$ at any given time $t$ can be computed by adding the integrals of instability growth rate over time in both regimes as

(8.9)\begin{equation} \zeta_{f} = \int_{0}^{t_{t}} \zeta_{0}s\,{\rm e}^{st} \,{\rm d} t + \int_{t_{t}}^{t_{f}} F\sqrt{(1-\epsilon )a\lambda} \,{\rm d} t = \zeta_{0}\,{\rm e}^{s t_{t}} + (t_{f}-t_{t}) F\sqrt{(1-\epsilon )a\lambda}.\end{equation}

It is assumed that there is a $\zeta _{0,M}$ which is the initial disturbance necessary to maximize instability growth at the instability wavelength observed in shadowgraph photos or a numerical simulation, $\lambda _{ob}$. This is found by

(8.10)\begin{equation} \left. \frac{\partial\zeta_{f}}{\partial\lambda}\right|_{\lambda_{ob},\zeta_{0,M}} = 0. \end{equation}

Assuming surface tension is negligible, a simplified relation for $\zeta _{0,M}$ using (8.10) is

(8.11)\begin{equation} \zeta_{0,M} = \frac{Fr\sqrt{a(1-\epsilon)\lambda_{ob}}}{\exp{\dfrac{1}{2}(t_{f}\sqrt{2{\rm \pi} a (1-\epsilon)/(\lambda_{ob}(1+\epsilon))}+\ln({2{\rm \pi} a (1-\epsilon)/(\lambda_{ob}(1+\epsilon))}))}}.\end{equation}

Now, the developed model may be evaluated. For example, figure 16 shows drop 4 at time $t_{f} = 14.4\ \mathrm {\mu }$s, with $\lambda _{ob} = 990\ \mathrm {\mu }$m. The values $a$ and $\epsilon$ for this drop were computed to be $1.1\times 10^7\ {\rm m}\ {\rm s}^{-2}$ and 0.025, respectively. Here, $Fr = 0.23$ is assumed for an annular bubble (a planar bubble revolved around the drop axis). From (8.11), $\zeta _{0,M}=14.4\ \mathrm {\mu }$m for drop 4. Table 4 summarizes these results for drop 4, with those for drop 6B, based on the surface instability wavelength measured in shadowgraph images and simulations. For both drops, the transition time $t_{t}$ is well below the time of drop impact, indicating that, according to this model, much of the drop's evolution is governed by nonlinear growth.

Table 4. Observed wavelength $\lambda _{ob}$, transition time $t_t$, initial disturbance $\zeta _{0,M}$, nonlinear growth rate predicted by (8.5) for experimental and simulation data. Observed growth rate for simulations from linear fits in figure 22, noting growth rate not directly observed in experiments.

In figure 21(b), the blue curve shows the growth rate according to this model for drop 4 given $\lambda _{ob} = 990\ \mathrm {\mu }$m from the shadowgraphs and with a predicted initial disturbance amplitude of $\zeta _{0,M}=14.4\ \mathrm {\mu }$m. The red curve uses the same model but the input wavelength is not that observed from the shadowgraphs, but the wavelength that maximizes the linear growth rate in figure 16(a), which is $\lambda _{lin} = 28\ \mathrm {\mu }$m; here, the initial disturbance amplitude is predicted to be $\zeta _{0,M}=0.25$ nm. Finally, the green curve represents the result of assuming $\lambda _{lin} = 28\ \mathrm {\mu }$m with the initial disturbance amplitude estimated to be $\zeta _{0,M}=14.4\ \mathrm {\mu }$m from the prior case. Note the smaller effective growth rate and transition time for cases with $\lambda _{lin} = 28\ \mathrm {\mu }$m compared with the case with $\lambda _{ob} = 990\ \mathrm {\mu }$m. This model predicts that small-wavelength disturbances transition to the nonlinear regime relatively quickly, effectively capping the growth rate of those disturbances that would otherwise dominate. To highlight this effect, the time integral of the data in figure 21(b) yields the calculated instability magnitude, which is shown in figure 21(c).

The plausibility of the calculated initial disturbance amplitudes can be assessed. An upper bound for initial disturbance amplitude may be estimated from the experiments by inspecting the pre-shot shadowgraphs, as those data appear smooth. Given camera resolution, the initial disturbances may be no larger than 50 $\mathrm {\mu }$m. A lower bound could be estimated by considering the conditions for which the flow local to a disturbance violates the continuum assumption. That length scale is estimated to be $\approx$1–10 nm by considering the local Knudsen number to be of order unity. These upper and lower bounds bracket our predicted initial disturbance amplitudes.

The simulations can be used to assess the efficacy of this model, even if the experimental design precludes direct comparison of the growth rate. To track the mean location of the drop surface, an ellipse was fitted to the windward face as assumed to be marked by the position of the 0.01 volume fraction ($\alpha _l = 0.01$) isosurface. These data were used to compute the acceleration for each case. The amplitude of a disturbance was found by tracking the difference between adjacent ripple peaks and troughs. The adjacent peak/trough pairs were tracked backward in time via a Gaussian peakfinder (O'Haver 1997). The initial guess for each peak location is given by the prior timestep, as is done in tagging velocimetry (Segall et al. 2023). Figure 26 in Appendix D shows a snapshot of peak/trough fitting for drop 6B. Figure 22 shows the computed disturbance amplitude as a function of time for the simulations of drops 4 and 6B (the simulation of drop 3 is used for drop 4 because the flow conditions are nearly identical between the two drops). The computational data indicate a clear transition between growth regimes. The agreement between the model and the simulations during the linear stage is difficult to assess given the model assumptions and spatial resolution in the computations; the same is true of the existence and timing of transition between regimes. Finally, the model/computational agreement for the nonlinear-instability growth rate is good in each case. The value of $Fr \sqrt {\lambda a}$ ((8.5)) computed in each case matches the computations to within 25 % (for case 6B) and to within 5 % (for case 4); this was determined by comparing computed nonlinear growth rate with the slopes of the linear fits in figure 22.

Figure 22. Non-dimensionalized instability amplitude $\zeta /d_{d}$ vs non-dimensional time $t/t_{c2}$ for simulated (a) drop 4 and (b) drop 6B.

The qualitative agreement between the simple nonlinear model from Pilch et al. (1981) and the simulations brings confidence to the approach. Furthermore, the assertion that much of the aerobreakup process observed in these experiments is likely in the nonlinear regime is strengthened. Finally, the authors believe this is a plausible hypothesis to explain why the wavelengths observed in experimental shadowgraphs are far larger than those predicted by linear-stability analysis.

9. Conclusion

Experiments and simulations of drop aerobreakup at high Mach and Weber numbers were carried out. High-speed backlit shadowgraph images of drops processed by bow shock waves at Mach numbers from 3.03 to 5.12 at sea-level ambient pressures were generated, and these images were compared with numerical simulations of drop aerobreakup. The high-speed images collected in this experimental campaign represent a potentially unique data set for these flow conditions.

A one-dimensional reacting gas dynamics model was applied to predict the conditions behind the bow shock wave of the projectile and the drop. That is, a double stagnation point flow problem was solved and it was found that the gas is nearly calorically perfect with a ratio of specific heats of $\gamma =1.3$. This assumption was used to simplify the numerical simulations. Moreover, the aerobreakup Weber number, $We_A$, and Reynolds number, $Re_A$, increase along the stagnation streamline of the projectile by ${\approx }50\,\%$ and ${\approx }25\,\%$, respectively. Although these increases are appreciable, they most likely do not alter the drop-breakup regime, making this data set comparable to shock-tube aerobreakup studies. Additionally, accounting for the drop bow shock does change the Weber and Reynolds numbers, but, also most likely, not enough to change the breakup regime. The parameters $We_A$ and $Re_A$ are evaluated at state 2, behind the bow shock of the projectile, as a valid state of comparison for this work with the literature. Finally, $We_A = \rho _g U_{LF}^{2} d_{d} / \sigma$, so it essentially parameterizes the aerodynamic pressure on the drop; that is, the aerodynamic pressure increases ${\approx }50\,\%$ along the stagnation streamline of the projectile.

Experimental data for drop deformation and position were acquired and compared with simulations. Drop displacement data were used to approximate the drag coefficient for a liquid drop in aerobreakup. The approximate drag coefficient $C_{d} = 2.32$ is consistent with the range of drag coefficients found in Pilch et al. (1981) and Hébert et al. (2020), among others. Drop displacement data were used to compute the acceleration to which the drops were subjected and drop deformation data were used to compute the degree of drop flattening before impact in support of drop surface instability analysis. Numerical simulations of drop aerobreakup reach similar values of $C_{d}$ and drop acceleration to those found in the experiments.

Observations of the shock structure and wake dynamics were made. Several differences between the shock structure around a deforming liquid drop and a rigid, solid sphere were identified. No separation shock at the drop equator is observed for deforming liquid drops. A vortical protuberance in the flow field is observed in the wake of liquid drops. A Mach-number dependence on the stability of the wake potentially explains the flow topology. Higher-Mach-number flows result in less chaotic wake structures with less apparent turbulent flow. Mass shedding from the drop equator is observed in simulations and shadowgraphs of liquid drops.

Linear-stability theory results from Funada & Joseph (2001) were used to assess the dependence of growth rate on (a) location on the drop face, (b) instability wavelength, (c) relative contributions of shear, acceleration, surface tension and viscosity and (d) the shape of the drop. The calculations in this work determined that viscosity and surface tension likely do not play a significant role at these conditions. At the stagnation point, the acceleration term dominates the shear term, and this relation inverts as one moves from the stagnation point to the drop equator. Finally, drop shape has a significant role in determining which of the terms dominates the problem. That is, for flattened drops, the acceleration term's dominance tends to broaden and move towards the drop equator, minimizing the shear term's importance.

We identified a long-standing discrepancy in the literature where some researchers assert that linear-instability growth results in rapid and sudden drop breakup when instability amplitude becomes comparable to drop thickness (Reinecke & McKay 1969; Joseph et al. 1999). Others assert that no such phenomenon is present and that drop breakup is governed exclusively by mass stripping (Theofanous & Li 2008; Theofanous 2011). Surface instabilities predicted from Funada & Joseph (2001) differ in amplitude and wavelength from those observed in either shadowgraph images or numerical simulations. Linear-stability analysis was insufficient for modelling late-stage aerobreakup because the predicted wavelengths were too small and the expected aerobreakup times were non-physically short. A model for nonlinear growth at large instability amplitudes was found in the literature (Pilch et al. 1981) and further developed. This model treats the nonlinear (constant rate) surface instability growth as an analogue to bubble rise in a mass of liquid, which is common in nonlinear RT treatment (Kull 1991; Zhou et al. 2021). The computed time of transition between the linear- and nonlinear-instability growth regimes was found to be less than $\tau = 1$ for all drops studied, indicating that much of the aerobreakup process takes place in the nonlinear-instability growth regime. An analytical solution for the initial instability magnitude was derived; initial disturbances inferred from this analysis were found to be plausible for the experimental conditions ($\zeta _0\approx 1\unicode{x2013}10\ \mathrm {\mu }$m). This nonlinear-instability growth model may explain results for drop breakup in prior literature thought (or alleged) to be anomalous. That is, we are implying that the results from Reinecke & McKay (1969), Joseph et al. (1999), Theofanous & Li (2008) and Theofanous (2011) are not in conflict, as is sometimes understood, but the results are complementary.