Study on the Interface Instability of a Shock Wave–Sub-Millimeter Liquid Droplet Interface and a Numerical Investigation of Its Breakup

Abstract

This study investigated the influence of instability on the interaction between sub-millimeter liquid droplets and shock waves. Experiments were conducted using 0.42 mm diameter droplets with varying shock wave Mach numbers. The investigation quantified the effects of Weber numbers and initial diameters on the development of Rayleigh–Taylor and Kelvin–Helmholtz instabilities at the shock wave–sub-millimeter liquid droplet interface. Three-dimensional numerical simulations were performed to investigate the deformation and breakup behaviors of sub-millimeter liquid droplets under the impact of a shock wave with a Mach number of 2.12. The post-shock gas flow environment in this condition was in a supersonic state. The simulations utilized the volume-of-fluid method to model the gas–liquid interface, employed unsteady Reynolds-averaged Navier–Stokes methods to simulate turbulence, and incorporated grid gradient adaptive technology to enhance computational efficiency. The results revealed that by increasing the Weber number or decreasing the initial diameter, both the growth rate and the wavenumber extremum of the Rayleigh–Taylor and Kelvin–Helmholtz instability waves increased. The variation in the K–H instability’s growth rate extremum increasing Weber number surpassed that of the R–T’s instability. This indicated that both the R–T and K–H waves on sub-millimeter liquid droplets tended to exhibit increased growth rates and reduced scales. Moreover, as the Weber number increased, the K–H instability became dominant in the aerodynamic fragmentation. The numerical simulations showed good qualitative agreement with the experimental data, affirming the viability of numerical methods for addressing such challenges. The evolution of the sub-millimeter liquid droplets was marked by two primary stages, flattening and shear stripping, signifying that the K–H instability-driven SIE mechanism governed the aerodynamic breakup in the supersonic post-shock airflow.

1. Introduction

The research into the interaction between shock waves generated by hypersonic airflow and liquid droplets is involved in various fields, including liquid atomization, multiphase combustion, and discrete multiphase flow. Meanwhile, such issues also have practical applications in daily life and engineering. For example, during high-speed aircraft flights through thunderstorms, investigations on the erosive effects of raindrops on the aircrafts’ surface coatings have shown reduced aircraft corrosion [1]. In natural disaster prevention, the dispersion of inert liquid substances under shock wave influences has ensured their wide distribution, effectively preventing fires and their spread. Moreover, due to the complex shock waves present in a combustion chamber when a hypersonic aircraft is in operation, the fragmentation and atomization of liquid droplets after interactions with shock waves affect the subsequent combustion stages and determine an engine’s performance [2]. Consequently, the investigation of liquid droplet deformation and fragmentation under shock wave influences holds substantial academic and engineering significance.

The deformation and aerodynamic breakup of liquid droplets under shock wave influence represents a classic phenomenon in multiphase fluid dynamics. Initially, research was primarily based on shock tube experiments and focused on the deformation and breakup of Newtonian fluid droplets in low-velocity flow conditions. The goal was to characterize various droplet breakup mechanisms, measure breakup times, and quantitatively assess the impacts of a shock wave’s intensity and viscosity ratio, among other parameters [3,4,5,6,7,8,9,10]. Theoretically, the Weber number (We), which represents the ratio of disruptive aerodynamic and restorative surface tensions ($W_e={\rho }_g{u}_g{}^2{d}_o/\sigma$), and the Ohnesorge number (Oh), which indicates the ratio of the liquid viscosity and surface tension ($O_h={\mu }_l/\sqrt{\sigma d_o{\rho }_l}$), serve as pivotal parameters for describing the aerodynamic breakup of liquid droplets. It is commonly accepted that when Oh < 0.1, the influence of liquid viscosity can be neglected so that the aerodynamic breakup process of liquid droplets is practically controlled by the Weber number [6]. Pilch and Erdman [7] classified five modes of liquid droplet aerodynamic breakup based on Weber numbers, which included vibrational breakup, bag breakup, multiple breakups, sheet-stripping breakup, and catastrophic breakup. Subsequently, as experimental techniques advanced, Theofanous [11,12,13] utilized laser-induced fluorescence (LIF) to address the challenge of liquid droplet obscuration during observations. Their research revealed that at a high Weber number, the aerodynamic breakup of liquid droplets results from shear-induced interface detachment, thus refuting the conventional catastrophic breakup mechanism. Based on these findings they introduced two novel droplet breakup mechanisms: RTP (Rayleigh–Taylor piercing) and SIE (shear-induced entrainment). While experimental research on the interactions between shock waves and liquid droplets has continued to advance, it is noteworthy that liquid droplets often have relatively large initial diameters. For instance, Poplavski [14] examined the breakup process of 2.7 mm water droplets under shock waves with Mach numbers ranging from 1.109 to 1.34. However, there remains a lack of systematic research on the aerodynamic breakup of smaller-scale liquid droplets, such as sub-millimeter-sized droplets, under the impact of a shock wave. In recent years, Shi et al. [15] shifted their research focus to explore the impact of shock waves on sub-millimeter liquid droplets. They devised a spatiotemporal evolution model based on high-speed schlieren imaging for various modes of sub-millimeter liquid droplet fragmentation. More recently, Zhang et al. [16] employed a combination of dual-pulse laser holography and an experimental apparatus with shock wave-driven droplets, and they acquired high-resolution experimental image data at the micron-nanosecond level.

Considering that the deformation and breakup of liquid droplets are associated with fluid dynamic interface instability, further research has been initiated. Commonly encountered fluid dynamics interface instabilities include Richtmeyer–Meshkov (R–M), Rayleigh–Taylor (R–T), and Kelvin–Helmholtz (K–H) instabilities. Xiang [17] suggested that at high Mach numbers, the rapid growth of high-wavenumber disturbances during the breakup of Newtonian liquid droplets results from weak viscous dissipation and surface tension. The Rayleigh–Meshkov instability can be ignored due to the significant acoustic impedance mismatch at the interface, with the primary mechanisms affecting droplet breakup being the Rayleigh–Taylor and Kelvin–Helmholtz instabilities. Jalaal and Mehravaran [18] attributed the growth and breakup of liquid sheets sandwiched between shearing airflows during droplet fragmentation to the development of Rayleigh–Taylor instability waves on the sheet. Additionally, Theofanous [19,20] conducted an experimental verification of Rayleigh–Taylor piercing during the droplet breakup process and quantified critical physical parameters related to shear-induced entrainment. Their findings highlight the predominant role of instability wave penetration in droplet breakup, particularly in systems with a high gas–liquid density ratio and high surface tension. Their research primarily concentrates on investigating the influence of viscosity and density ratios on the development of instability. However, a quantitative analysis of instability development at the liquid droplet interface under other parameters is conspicuously absent.

Advances in computational capabilities and technology have transformed computational fluid dynamics (CFD) into a valuable complement to experimental and analytical methods. Compared to experiments, CFD numerical simulations excel at capturing the intricate details of flow fields and gas–liquid interfaces. However, the high cost of three-dimensional numerical calculations has led to the predominance of simulations conducted in two dimensions [21,22] or in a two-dimensional axisymmetric setting [23]. Meng et al. [24] performed numerical simulations analyzing liquid droplet characteristic breakup and edge tip generation. They also investigated the SIE mechanism, elucidating the shear, peel, and fracture processes of liquid droplets and the flow field evolution. Luo [25] employed the volume-of-fluid (VOF) method and turbulence models for numerical simulations, revealing that the Weber number enhanced liquid droplet deformation and breakup. Shi et al. [26] employed the VOF multiphase flow model and the k-ε turbulence model; they analyzed, using two-dimensional numerical simulations, the impact of different Weber numbers on the deformation and evolution processes of sub-millimeter liquid droplets. Poplavski et al. [14] adapted dynamic grid technology and studied the structure of the flow near and in the wake of a drop, the features of the flow around a drop, the type of the shape evolution, and the character of the mass entrainment, based on the use of the volume-of-fluid (VOF) method to resolve the phase interface. Rossano et al. [27] creatively applied a hybrid approach using VOF and Lagrangian computational fluid dynamics methods to investigate liquid droplet fragmentation. Their research accurately predicted the deformation and fragmentation progression of liquid droplets, as well as the subsequent formation of liquid sprays, replicating essential features of droplet deformation and fragmentation under shear-induced entrainment conditions. Nevertheless, there is currently limited research on three-dimensional numerical simulations of the aerodynamic breakup of sub-millimeter-sized droplets.

While many studies have investigated droplet breakup with shock wave interactions through experiments and numerical simulations, there remains a notable gap in the understanding of the dynamic mechanisms governing sub-millimeter droplets generated after secondary fragmentation. Particularly, there is an absence of quantitative analysis concerning instability development at the sub-millimeter droplet–shock wave interface, as well as a dearth of accurate three-dimensional numerical simulations regarding the evolution of sub-millimeter droplets when the airflow environment is in a supersonic state after a shock wave interaction. Therefore, this study first conducted experiments on the aerodynamic fragmentation of sub-millimeter liquid droplets with an initial diameter of 0.42 mm under various shock wave impacts, based on the experimental designed by Zhang [16]. Quantitatively, we investigated the impact of sub-millimeter droplet initial diameter and Weber number on the development of the R–T and K–H instabilities following the extraction of spatiotemporal parameters for the fragmentation evolution of liquid droplets. Subsequently, a three-dimensional numerical simulation study was conducted to investigate the characteristics of the aerodynamic breakup of sub-millimeter liquid droplets in a supersonic gas velocity environment behind a Mach 2.12 shock. The VOF method was employed to capture the deformation of the liquid droplet interface and the evolution of the flow field. Moreover, adaptive grid refinement technology was employed to improve the accuracy and efficiency of the simulation. The numerical simulation result exhibited good qualitative agreement with the experimental data. These findings serve as a reference for further investigations into the development of interface instability and the dynamic mechanisms governing sub-millimeter liquid droplets generated after secondary fragmentation due to shock wave actions.

2. Instability of Shock Wave–Sub-Millimeter Liquid Droplet Interfaces

2.1. Experimental Parameters

Liquid droplet breakup in a shock tube experiment is predominantly caused by the substantial relative velocity produced around the droplet following the passage of a positive shock wave [22]. To investigate the influence of the initial droplet’s diameter on the instability of the shock wave–sub-millimeter droplet interface, sub-millimeter droplets with initial diameters of 0.42 mm were subjected to varying Mach numbers. The shock Mach number, defined as the ratio of the shock wave velocity to the local speed of sound, is expressed as follows:

$$M_a=\frac{u_s}{c},$$

(1)

where $u_s$ is the shock wave speed and $c$ is the local speed of sound.

Our research employed the experimental method developed by Zhang [16], which integrated dual-pulse laser holographic testing technology with a shock tube device. Comprehensive details regarding the experimental procedures and principles can be found in reference [16]. To quantitatively study the instability development of the initial sub-millimeter droplet’s diameter and Weber number on the shock wave–droplet interface, the sub-millimeter droplet size in the experiment was 0.42 mm. This comparative analysis utilized experimental data from Shi [15] (sub-millimeter) and Zhang [16] (0.24 mm). In addition, this experiment also investigated the evolution of sub-millimeter liquid droplet deformation and fragmentation under shock wave conditions with a Ma value of 2.12, where the airflow environment post-shock wave achieved supersonic velocities. The experiment employed two working fluids—air and water droplets—and the aerodynamic experimental parameters are presented in

Table 1. Experimental operating conditions.

Ma	p2/ (Pa)	ug/(m/s)	ρg/(kg/m3)	T2/(K)	We
1.09	123,701	48.83	1.38	315.52	19
1.13	138,337	69.34	1.47	323.12	42
1.27	174,391	136.58	1.76	349.44	196
1.49	247,241	231.75	2.22	391.82	713
2.12	514,480	466.5	3.4	533.64	4272

In the table, p₂, u_g, ρ_g, and T₂ represent the post-shock airflow pressure, velocity, density, and temperature, respectively.

. Further, the airflow and water droplet properties that were imposed as initial conditions are summarized in

Table 2. Pre-shock air and water droplet properties.

Properties	Symbol	Value
Water surface tension	σ	0.0728 N/m
Water density	ρl	998 kg/m3
Water viscosity	µl	1.003 × 10−3 Pa/s
Air temperature	T1	298 K
Air viscosity	µg	1.79 × 10−5 Pa/s
Air pressure	p1	101.325 kpa

.

2.2. The Influence of Weber Number and Initial Diameter on the R–T Instability

Under low Weber number conditions, initial perturbations manifested on the surfaces of the sub-millimeter liquid droplets prior to their fragmentation. These perturbations primarily underwent small-scale development, rendering their observation challenging. Nevertheless, with the ongoing interaction between the high-speed airflow and the droplet, the R–T instability underwent progressive evolution, culminating in the aggregation of large-scale R–T waves. The fundamental reason for the formation of one or more bag-like structures during bag breakup (We = 19) and multiple breakup (We = 42) was the continued piercing of the R–T instability waves, as depicted in Figure 1.

Following the passage of a shock wave over a sub-millimeter liquid droplet, the droplet was immediately exposed to high-speed post-shock airflow. Assuming exponential growth ($e^{nt}$) in the surface amplitude of the droplet, the dispersion equation for the R–T instability waves can be expressed as follows [28,29,30,31]:

$$n_{RT}=-k_{RT}{}^2\frac{{\mu }_l+{\mu }_g}{{\rho }_l+{\rho }_g}+{(k_{RT}\frac{{\rho }_l-{\rho }_g}{{\rho }_l+{\rho }_g}a-k_{RT}{}^3\frac{\sigma }{{\rho }_l+{\rho }_g}+k_{RT}{}^4{(\frac{{\mu }_l+{\mu }_g}{{\rho }_l+{\rho }_g})}^2)}^{0.5},$$

(2)

where $n_{RT}$ is the growth rate of the R–T instability waves, $k_{RT}$ is the wavenumber of the R–T instability waves, ${\mu }_l$ and ${\mu }_g$, respectively, represent the liquid viscosity and gas viscosity, and a is the acceleration of the liquid droplet’s center of mass, which is expressed as follows:

$$a=\frac{\mathrm{\Delta }{l}_2-\mathrm{\Delta }{l}_1}{\mathrm{\Delta }{t}^2},$$

(3)

where $\mathrm{\Delta }{l}_1{}_{}$ and $\mathrm{\Delta }{l}_2$ represent the displacements of the sub-millimeter liquid droplets within the time intervals shown in Figure 2.

Figure 3 illustrates how the growth rate of a sub-millimeter liquid droplet (with an initial diameter of 0.42 mm) in R–T instability waves changed as a function of the wavenumber. Initially, the growth rate of R–T instability waves increased with the wavenumber, followed by an overall decrease. This occurred due to the aerodynamic breakup of the sub-millimeter liquid droplets, resulting in the continuous detachment of sub-droplets from the main droplet. This resulted in droplet sizes that were insufficient to support the growth of the R–T instability waves. As the Weber number parameter increased, the extremum shifted toward the upper right. This indicated that the surface R–T instability waves of the sub-millimeter liquid droplets tended to develop in the direction of the increasing wavenumbers, decreasing scales, and increasing growth rates. Variations in growth rates among the different wavelengths would become prominent. At this stage, the substantial growth rate increase in the R–T instability waves accelerated the R–T instability development, rendering the sub-millimeter liquid droplets more prone to fragmentation.

In Figure 4, a comparison is presented between this study and the experiments conducted by Shi [15] and Zhang [16]. It shows how R–T instability wave growth rates vary with the wavenumbers for sub-millimeter liquid droplets of varying initial diameters under identical Weber number conditions. Evidently, decreasing the initial diameters of sub-millimeter liquid droplets leads to gradual increases in both the peak growth rates and wavenumbers of the R–T instability waves. This implies that, under a constant Weber number, reducing the droplet diameter causes the R–T instability waves on the droplet’s surface to evolve toward smaller scales and higher growth rates. The analysis indicated that the growth of R–T waves was primarily linked to acceleration, where a higher acceleration led to faster R–T instability development. Under a constant Weber number, reducing the initial diameters of sub-millimeter liquid droplets increased the airspeeds around the droplets, consequently leading to greater acceleration for smaller droplets compared to larger ones in the same time frame. Consequently, experimental observations revealed that sub-millimeter liquid droplets with smaller initial diameters displayed reduced-scale surface-piercing behaviors and underwent more rapid development. Furthermore, it became apparent that with increasing Weber numbers, the disparities in the peak growth rates and maximum wavenumbers of the R–T instability waves between sub-millimeter liquid droplets with different initial diameters gradually diminished. This suggested that when the Weber number exceeded a specific threshold, variations in the initial diameters of sub-millimeter liquid droplets had a relatively minor influence on the progression of R–T instabilities.

2.3. The Influence of Weber Number and Initial Diameter on K–H Instability

When air flows over the surface of a sub-millimeter liquid droplet, a tangential velocity difference between the gas–liquid phases induces shear forces that generate Kelvin–Helmholtz instability waves at the interface. These waves continue to develop over time. Figure 5 depicts the evolution of K–H instability waves. With an increasing Weber number, R–T instability is suppressed and the K–H instability progressively intensifies. Since K–H instability primarily affects the equatorial plane of a sub-millimeter liquid droplet and exerts a more pronounced influence compared to R–T instability, the process involved droplet shear and stripping, followed by rapid downstream travel. This was distinct from the bag and multiple breakup mechanisms observed at low Weber numbers.

The growth process of a K–H instability wave directly affects the aerodynamic breakup of a liquid droplet [20]. While investigating the theory of interface instability, Wang [28] derived the following dispersion equation for planar K–H instability:

$$n_{KH}=-k_{KH}{}^2\frac{{\mu }_l+{\mu }_g}{{\rho }_l+{\rho }_g}+{(k_{KH}\frac{{\rho }_l-{\rho }_g}{{\rho }_l+{\rho }_g}a+\frac{k_{KH}{}^2{\rho }_l{\rho }_g\left(u_{gs}-u_{ls}\right)}{{({\rho }_l+{\rho }_g)}^2}-k_{KH}{}^3\frac{\sigma }{{\rho }_l+{\rho }_g}+k_{KH}{}^4{(\frac{{\mu }_l+{\mu }_g}{{\rho }_l+{\rho }_g})}^2)}^{0.5},$$

(4)

where $n_{KH}$ is the growth rate of the K–H instability waves, $k_{KH}$ is the wavenumber of the K–H instability waves, and $u_{gs}$ and $u_{ls}$ are the surface air velocity and fluid motion velocity of the droplet, respectively.

Figure 6 illustrates the relationship between the growth rate of K–H instability waves in 0.42 mm liquid droplets and wavenumber variations. Similar to the R–T instability waves, the growth rate of K–H instability waves initially increased and then decreased. Notably, the extremum of the K–H instability wave growth rate and the wavenumber variations were significantly more pronounced than they were in the R–T instability shown in Figure 3 (reaching a value of ${10}^6$). Moreover, as the Weber number increased, both the extremum of the K–H instability wave growth rate and the wavenumber steadily rose. This indicated that at low Weber numbers, K–H instability exhibited minimal growth, whereas under high Weber number conditions, the influence of K–H instability on sub-millimeter liquid droplet breakup surpassed that of R–T instability. In such cases, K–H instability predominated in sub-millimeter liquid droplet breakup.

Figure 7 compares the growth rates of the K–H instability waves with the wavenumbers for sub-millimeter liquid droplets with various initial diameters but under constant Weber number conditions. As the initial diameters of the sub-millimeter liquid droplets decreased, the K–H instability wave curve showed similarities to the R–T instability wave curve. Both the peak growth rates and the corresponding wavenumbers increased gradually. This trend implied that, under constant Weber conditions, decreasing the droplet diameter led to the development of K–H waves on the droplets’ surfaces at smaller scales and higher growth rates.

This corresponded to the experimental observation that K–H instability waves on sub-millimeter liquid droplet surfaces had smaller sizes and led to faster droplet breakup. K–H instability wave growth is primarily associated with the velocity difference between the gas and liquid phases. With an increasing velocity difference, K–H instability develops rapidly. When Weber numbers are equal, the velocity difference between the gas and liquid phases increases as the initial diameters of sub-millimeter liquid droplets decrease. Consequently, rapid K–H instability development occurs on the smaller droplets’ surfaces. Moreover, as Weber numbers increase, unlike R–T instability waves, the disparity between the maximum growth rate and the maximum wavenumber of the K–H instability waves in sub-millimeter droplets of varying diameters also expands. The analysis indicated that increasing Weber numbers and reducing the initial diameters of the sub-millimeter liquid droplets amplified the velocity difference between the gas and liquid phases, thereby fostering K–H instability development on the droplets’ surfaces and gradually increasing the effects on the aerodynamic breakup of droplets.

3. Numerical Simulation Study of Sub-Millimeter Liquid Droplets

Although there have been numerous numerical simulations investigating shock wave interactions with droplets [23,24,32], research on droplets with initial diameters at or below the sub-millimeter scale in a supersonic post-shock wave airflow environment is lacking. Consequently, it is imperative to conduct numerical simulations to analyze the aerodynamic breakup of individual sub-millimeter liquid droplets when exposed to high Mach number shock waves. We conducted simulations using the parameters detailed in Table 1.

3.1. Computational Analysis

A shock wave itself exerts limited influence on the breakup of sub-millimeter liquid droplets. Instead, the primary factor is the significant relative velocity generated in the flow field surrounding a droplet after a shock wave’s passage, initiating a sequence of deformations and fragmentations [33]. In practice, a droplet experiences acceleration under the influence of air that is thousands of times greater than gravitational acceleration. Consequently, gravity’s effect is usually neglected in the description of droplet evolution and quantitative parameter analysis [34]. When dimensionally analyzing droplet evolution, we can assess interactions between various parameters. The dimensionless windward displacement, denoted by $X$, is defined as the ratio of the distance moved along the airflow from the beginning of the windward to the initial droplet diameter ($S/d_0$). The dimensionless transverse deformation width, denoted as $d_c$, signifies the ratio of the width of the transverse deformation along the vertical airflow direction during the droplet’s evolution to the initial diameter ($\left(D-d_0\right)/d_0$). Furthermore, the dimensionless axial deformation width, denoted as $w$, is defined as the ratio of the axial deformation width in the direction of the airflow propagation to the initial diameter ($\left(d_l-d_0\right)/d_0$), as illustrated in Figure 8. To eliminate the influence of airflow-related parameters, the following formulas were used to nondimensionalize the fragmentation time:

$$T=u_g\frac{t}{d_o\sqrt{{\rho }_l/{\rho }_g}}$$

(5)

where t is the dimensional time, $d_o$ is the initial diameter of the droplet, and ${\rho }_l$ and ${\rho }_g$ represent the densities of the liquid phase and the gas phase, respectively. The moment when the shock wave contacts the droplet is defined as time zero.

In the simplified representation shown in Figure 9, the computational domain model depicts a left-to-right propagating shock wave. A sub-millimeter liquid droplet is symmetrically positioned on the central axis. The shock tube’s cross-sectional area is 60 mm × 60 mm, and the experimental section has a length of 500 mm. This study exclusively examined the physical attributes of the flow field around sub-millimeter liquid droplets, eliminating the need for a numerical simulation of the entire shock tube structure. The simplified calculation spatial domain size was chosen such that $-20<x/d_0<30$, $-10<y/d_0<10$, and $-10<z/d_0<10$. This approach not only mitigated the influence of reflected waves but also conserved computational resources, to some extent. Given the high flow field symmetry, a one-quarter model was selected for the numerical simulation. “Symmetry” boundary conditions were imposed on the two central planes containing sub-millimeter liquid droplets. At the left boundary of the computational domain, pressure inlet boundary condition was applied and at the right boundary, a pressure outlet was assigned. Other boundary conditions, which included the top and back boundaries of the computational domain, were set to symmetry. The flow field used structured grids, and the sub-millimeter liquid droplet section employed a grid gradient adaptation to improve the computational efficiency. Dynamic grid adaptation based on phase fraction gradient was implemented inANSYS Fluent V20.2. This involved adjusting the mesh dynamically by refining regions with high gradient values and coarsening regions with low gradient values. To achieve this, the refinement thresholds were set to refine the grids in areas with large gradients, while coarsening thresholds were set to revert the mesh in areas with small gradients. For unsteady flow environments, it was essential to control the entire computational time by standardizing the gradient values. In the scaling method, the “Scale by Global Maximum” option was applied. Finally, for the phase fraction gradient adaptation, specific refinement and coarsening values were provided, with a refinement value of 0.06, a coarsening value of 0.05, and a maximum refinement level that was set to 4. Grid refinement was automatically conducted during the calculations based on the liquid volume fraction gradient to achieve an accurate solution domain. Furthermore, under these thermodynamic conditions, the environmental pressure in the flow field exceeded the saturated vapor pressure of the water. Consequently, during the sub-millimeter liquid droplet deformation and breakup, no phase change occurred, and it could be assumed that the gas–liquid two-phase was immiscible.

In this research, the multiphase flow model selected was the VOF (volume-of-fluid) model due to its excellent mass conservation properties and effective interface capturing capabilities. It can accurately describe mass transfer in the fluid and interface diffusion [14,32,35]. This method is based on the idea of treating liquid and gas as a single two-component medium. Within the computational domain, the spatial distribution of the liquid phase was characterized by a step function, denoted as F (x, y, z, t), with values that ranged between 0 and 1. If the entire computational grid was occupied by the liquid phase, F (x, y, z, t) was equal to one, if it was entirely occupied by the gas phase, F (x, y, z, t) was equal to zero. In cases where both the liquid and gas phases coexisted within the computational grid, 0 < F (x, y, z, t) < 1. As the value of F (x, y, z, t) changed over time, solving the transport equation for the volume fraction of the liquid phase within the grid allowed for the determination of the distribution of the two-phase system. This, in turn, enabled the tracking of the evolution of the liquid phase interface.

Our simulation calculations assumed that the liquid phase was incompressible and treated the gas phase as an ideal gas governed by the gas state equation. Gravity effects were not taken into consideration. The governing equations included the continuity equation, momentum equation, and energy equation, as follows:

$$\frac{\partial }{\partial t}({\alpha }_g{\rho }_g)+\frac{\partial }{\partial x_j}({\alpha }_g{\rho }_g{u}_{g,j})=0$$

(6)

$$\frac{\partial }{\partial t}{\alpha }_l+\frac{\partial }{\partial x_j}{\alpha }_l{u}_{l,j}=0.$$

(7)

In the above continuity equation, the subscripts $g$ and $l$ correspond to the gas phase and liquid phase, respectively, and $u$ represents velocity, $x_j$ is a Cartesian coordinate component, and ${\alpha }_l$ and ${\alpha }_g$ denote the volume fractions of the liquid phase and gas phase, respectively, which were needed to satisfy Equation (8).

$$\{\begin{array}{l}{}^{{\alpha }_g+{\alpha }_l=1}\\ \rho ={\alpha }_g{\rho }_g+{\alpha }_l{\rho }_l\end{array},$$

(8)

$$\frac{\partial }{\partial t}({\alpha }_g{\rho }_g{u}_{g,i})+\frac{\partial }{\partial x_j}\left({\alpha }_g{\rho }_g{u}_{g,i}{u}_{g,j}\right)=-{\alpha }_g\frac{\partial p}{\partial x_j}+\frac{\partial {\mathrm{\tau }}_{g,ij}}{\partial x_j}+\frac{\partial {\mathrm{\tau }}_{gs,ij}}{\partial x_j}+\frac{{\alpha }_g{\rho }_g}{{\mathrm{\tau }}_r}(u_{g,i}-u_{l,i})+{\alpha }_g{\rho }_g{g}_j-F_{g,sj},$$

(9)

$${\rho }_l\left[\frac{\partial }{\partial t}({\alpha }_l{u}_{l,i})+\frac{\partial }{\partial x_j}\left({\alpha }_l{u}_{l,i}{u}_{l,j}\right)\right]=-{\alpha }_l\frac{\partial p}{\partial x_j}+\frac{\partial {\mathrm{\tau }}_{l,ij}}{\partial x_j}+\frac{\partial {\mathrm{\tau }}_{ls,ij}}{\partial x_j}+\frac{{\alpha }_l{\rho }_l}{{\mathrm{\tau }}_r}(u_{g,i}-u_{l,i})+{\alpha }_l{\rho }_l{g}_j-F_{l,sj}.$$

(10)

In the momentum equation, ${\tau }_{m,ij}$ is the molecular stress tensor and ${\tau }_{s,ij}$ is the Reynolds stress tensor whose components have the following form:

$${\mathrm{\tau }}_{m,ij}={\mathrm{\mu }}_m\left(\frac{\partial u_{m,i}}{\partial x_j}+\frac{\partial u_{m,j}}{\partial x_i}-\frac{2}{3}{\mathrm{\delta }}_{ij}\frac{\partial u_{m,l}}{\partial x_l}\right)(m=g,l)$$

(11)

$${\mathrm{\tau }}_{ms,ij}={\mathrm{\mu }}_t\left(\frac{\partial u_{m,i}}{\partial x_j}+\frac{\partial u_{m,j}}{\partial x_i}\right)-\frac{2}{3}{\mathrm{\mu }}_t{\mathrm{\delta }}_{ij}\frac{\partial u_{m,l}}{\partial x_l}-\frac{2}{3}{\rho }_{m}k{\mathrm{\delta }}_{ij}(m=g,l),$$

(12)

where ${\tau }_r$ is the velocity relaxation time, $F_{sj}$ is the surface tension, $\mu$ is the molecular dynamic viscosity, ${\mu }_t$ is the turbulent viscosity coefficient, $m=g$ denotes the gas phase, and $m=l$ denotes the liquid phase.

$$\frac{\partial }{\partial t}(\rho E)+\frac{\partial }{\partial x_j}\left[u_j\left(\rho E+p\right)\right]=\frac{\partial }{\partial x_j}\left(k_{eff}\frac{\partial T}{\partial x_j}\right)+\frac{\partial }{\partial x_j}\left({\mathrm{\tau }}_{eff,ij}\cdot u_i\right).$$

(13)

In the above energy equation, E is the total energy, k_eff is the effective thermal transmittance, and the effective viscous stress tensor, τ_eff,ij, is given as follows:

$${\mathrm{\tau }}_{eff,ij}={\displaystyle \sum _{m=g,l}{\alpha }_m}{\rho }_m\left({\mathrm{\tau }}_{m,ij}+{\mathrm{\tau }}_{ms,ij}\right)$$

(14)

In numerical simulations, where the gas phase is considered an ideal gas, it is imperative to augment the ideal gas state equation to complete the control equations mentioned above as follows:

$$p={\rho }_{g}RT$$

(15)

where p represents the gas pressure, R is the gas constant, and T stands for the temperature.

3.2. Computational Model Establishment and Verification

The numerical calculations were conducted using the industrial simulation software ANSYS Fluent V20.2. It has been successfully applied to past aerodynamic fragmentation studies on liquid droplets under shock wave impacts [2726,3234]. Turbulent flow simulations were conducted by applying the k-w SST model in the context of unsteady Reynolds-averaged Navier–Stokes (RANS) equations. The solution method involved the use of a pressure-based coupled solver with the pressure–velocity coupling resolved through the PISO algorithm, which has been recognized for its efficiency in achieving convergence with fewer iterations. Spatial discretization of the pressure field was performed using the robust PRESTO! scheme, which is well suited for multiphase flow problems. The density, momentum, and turbulent kinetic energy were discretized with the third-order MUSCL scheme. The initial configuration of the sub-millimeter liquid droplets and the flow field was established with the Patch function. This function partitioned the computational domain into regions with uniform high- and low-pressure parameters to set the initial conditions for the propagating shock waves. The initial parameters on both sides of the shock waves were determined following the shock tube principles. Based on the flow conditions, a constant time step size of $\mathrm{\Delta }t=1\times {10}^{-8}$ s was defined for the explicit transient simulations.

Various grid densities were utilized to evaluate the grid independence in constructing the fluid domain, ensuring that the density had no impact on the computational outcomes. The non-dimensional displacement of the sub-millimeter droplets on the windward side was compared at a Mach number of 2.12 using various grid densities. CFD I utilized 1,686,528 cells, CFD II employed 2,502,892 cells, and CFD III featured 3,015,726 cells. Figure 10 demonstrates that there were no significant differences in the non-dimensional displacements of the sub-millimeter droplets over time for the three grid densities. To ensure both manageable computational demands and the complete capturing of the gas–liquid interface evolution, this study chose the grid density of CFD II for the subsequent three-dimensional numerical simulations of the sub-millimeter droplet aerodynamic breakup.

In order to validate our numerical simulation method and ensure the accuracy of the subsequent computational results, we applied this method to calculate the dimensionless transverse deformation widths for two specific parameter conditions, C2 (Ma = 1.17 and $d_0$ = 0.9) and D1 (Ma = 1.36 and $d_0$ = 0.84), as defined by Shi [26] for vertical airflow. The results are depicted in Figure 11. A comparison between these numerical simulation results and the experimental data highlighted the improved consistency when employing the numerical simulation method introduced in this paper. Nevertheless, due to the inherent uncertainties in the experimental measurements, there was no definitive value for the liquid phase volume fraction, ${\alpha }_l$, that best fit the data. This variability contributed, to some extent, to the disparity between the numerical simulation results and the experimental outcomes.

3.3. Results and Discussion

Figure 12 presents the morphological evolution of the windward and leeward surfaces at eight discrete time instances from the dimensionless times 0 to 1.497. The shapes corresponded to liquid phase isosurfaces with volume fractions of ${\alpha }_l$ = 0.99, 0.50, and 0.01 at Ma = 2.12 operating conditions. Under these operating conditions, the airflow attained a supersonic state with the post-shock Mach number $M_2=1.03$. This supersonic environment amplified the instant deformation and aerodynamic breakup of the sub-millimeter liquid droplets. During the initial stage ($T=0~0.191$), the sub-millimeter liquid droplets remained largely undeformed. Surface disturbances arising from the instability began to develop gradually. At this point, the droplets could be approximated as rigid bodies with no significant axial displacements along the airflow direction. When $T=0.322~0.518$, sub-millimeter liquid droplets experienced initial distortion and deformation. The windward side retained a relatively spherical shape, while the leeward side was rapidly compressed into a flat shape. Subsequently, when $T=0.779~1.170$, vortical structures continuously developed and expanded on the leeward sides of the sub-millimeter liquid droplets due to the airflow, resulting in the formation of noticeable cavities in the leeward regions. Concurrently, the sub-millimeter liquid droplets underwent the strong shearing effect of the K–H instability. Small-scale liquid sheets and sub-droplets continually detached from the equatorial regions and downstream of the leeward sides of the droplets. High-speed airflow accelerated and stretched them away from the parent droplets. This process, influenced by flow instability, resulted in the thinning and gradual development of ligaments from the thinner parts of the liquid sheets. Over time, these ligaments became finer, and their collisions further contributed to the droplets’ fragmentation into smaller liquid droplets. Consequently, in the later stages of fragmentation ($T=1.497$), a widespread distribution of sub-droplets became apparent in the downstream regions. Through a comprehensive examination of the sub-millimeter liquid droplet characteristics, it became evident that the deformation and breakup of the sub-millimeter liquid droplets in a supersonic airflow environment aligned with the SIE breakup mechanism as proposed by Theofanous [11].

The comparison between the numerical simulation results and the experimental visualizations at corresponding time points is presented in Figure 13. The figure illustrates a close correspondence between the deformation and breakup morphologies of the sub-millimeter liquid droplets obtained through the numerical simulations and the experimental photographs. Additionally, Figure 14 provides a parameter comparison, depicting the dimensionless transverse deformation widths of the sub-millimeter liquid droplets in the vertical airflow direction over time. The results revealed good qualitative agreement between the numerical simulations and the experimental measurements, underscoring the effectiveness of our numerical simulation method in capturing the interface evolution during the sub-millimeter liquid droplet breakup process.

Figure 15 displays the streamline, velocity contours, and pressure contours illustrating the deformation and fragmentation processes of the sub-millimeter liquid droplets subjected to shock wave conditions with an Ma value of 2.12. Small surface waves formed on the windward surfaces and propagated toward the equators from the initial moment. On the leeward surfaces of the sub-millimeter liquid droplets, airflow separation occurred, resulting in the formation of small vortexes. At this point, the airflow near the windward surfaces demonstrated high pressures and low velocities. The pressures on the windward surfaces were not uniform, and the surface waves created high-pressure stagnation points on the windward sides by impeding the airflow behind the waves. In the equatorial regions of the sub-millimeter liquid droplets, low pressures and high velocities prevailed, with the highest airflow speed observed at this position. When considering the droplets’ stagnation states, it became evident that the maximum tangential velocity differences between the gas and liquid phases occurred on the equatorial planes, resulting in the onset of Kelvin–Helmholtz instability. Subsequently, due to the significantly lower pressures in the equatorial regions compared to the pressures at the front and rear stagnation points, the sub-millimeter liquid droplets underwent flattening. This flattening reduced the curvatures on the windward sides of the droplets, leading to increased variations in the airflow velocities as they traversed the droplets. This, in turn, promoted the continuous growth of vortices on the leeward sides and propelled the propagation of the surface waves towards the equatorial planes. When examining Figure 12, it is evident that the sub-millimeter liquid droplets expanded laterally in the vertical airflow direction, while the leeward regions adopted concave shapes, causing sub-droplets to continually detach from the vicinity of the equatorial planes.

The velocity contours revealed flame-like tails extending in the airflow direction on the leeward sides of the sub-millimeter liquid droplets. In the influence of the vortices, there were observable annular low-speed regions in the middles of the velocity distributions on the leeward sides, corresponding to the outer vortex regions. Finally, the flow fields and vortex structures on the leeward sides continue to develop and became more complex. The main sizes and masses of the droplets decreased continuously due to the detachment of small droplets, while the annular low-speed regions at the tails enlarged continuously due to the expansion of the vortexes.

A quantitative investigation into the deformation process of the sub-millimeter liquid droplets was subsequently carried out. Figure 16 illustrates the dimensionless axial deformation widths aligned with the flow direction and the dimensionless transverse deformation widths aligned with vertical airflow over time. The results are presented for volume fractions of 0.01, 0.5, and 0.99. The axial widths of the sub-millimeter liquid droplets decreased approximately linearly with time in the airflow direction, while the evolution of the transverse widths in the vertical airflow direction underwent the following two stages: a slow flattening deformation stage and a strong shear stretching deformation stage, which aligned with the characteristics observed in the numerical simulations for the sub-millimeter liquid droplets’ deformation and fragmentation, as shown in Figure 12.

4. Conclusions

In this study, we conducted experimental investigations into the aerodynamic breakup of sub-millimeter liquid droplets with initial diameters of 0.42 mm under varied shock wave conditions. It quantitatively assessed the development of Rayleigh–Taylor and Kelvin–Helmholtz instabilities at different Weber numbers and initial diameters. Furthermore, three-dimensional numerical simulations, based on experimental parameters, were conducted to examine the evolution of the sub-millimeter liquid droplets and the surrounding flow fields in a supersonic post-shock gas velocity environment. The analysis centered on the deformation characteristics of these droplets and the gas–liquid interface evolution during aerodynamic breakup. The present work yielded the following conclusions:

For constant initial diameters, the growth rates of the R–T and K–H instability waves on the sub-millimeter liquid droplets’ surfaces initially increased and then decreased with varying wavenumbers. Moreover, with an increase in the Weber number, both types of instability waves exhibited inclinations towards higher growth rates and smaller scales. The differences in the growth rate extremes were notably more pronounced for the K–H instability compared to the R–T instability.

At a constant Weber number, decreasing the initial diameter of sub-millimeter liquid droplets led to gradual increases in both the maximum growth rate and the wavenumbers of the Rayleigh–Taylor and Kelvin–Helmholtz instability waves. Additionally, as the Weber number increased, the disparities between the maximum growth rates and the wavenumbers of the R–T instability waves in the sub-millimeter liquid droplets with varying initial diameters decreased, while the K–H instability waves exhibited an increasing trend.

The numerical simulations revealed that, under the supersonic gas velocities behind a Mach 2.12 shock, the robust shear forces at the gas–liquid interface induced the formation of K–H instability waves. During this phase, the sub-millimeter liquid droplets underwent aerodynamic fragmentation, adhering to the SIE mechanism and comprising the following two sequential phases: flattening and shear stripping. The flattening phase deformation arose from the pressure gradients between the high-pressure regions ahead of the sub-millimeter liquid droplets and the low-pressure zones encircling their equators. In the shear stripping phase, fragmentation occurred due to the continual development of vortex structures on the leeward sides, which resulted from the flow separations at the equatorial planes.