Large-Eddy Simulation of Droplet Deformation and Fragmentation Under Shock Wave Impact

Abstract

This study employs the large-eddy simulation (LES) approach, together with the hybrid volume of fluid—discrete phase model, to examine the deformation and breakup of a water droplet impacted by a traveling shock wave. The research investigates the influence of Weber number on transient deformation and breakup characteristics. Particular focus is given to the detailed analysis of sub-droplet-size distributions, which are frequently overlooked in existing studies, providing a novel insight into droplet fragmentation dynamics. The predicted deformation and breakup patterns of droplets in the shear breakup regime align well with experimental data, validating the computational approach. Notably, LES is able to reproduce the underlying physical mechanisms, highlighting the significant role of recirculation zones and the progression of Kelvin–Helmholtz instabilities in droplet breakup. Additionally, it is shown that higher Mach numbers significantly amplify both cross-stream and streamwise deformations, leading to earlier breakup at higher airflow pressures. Increasing the Weber number from 205 to 7000 results in 25% reduction in the average size of the sub-droplets, indicating the strong influence of aerodynamic forces on droplet fragmentation. This comprehensive analysis, while aligning with experimental observations, also provides new insights into the complex dynamics of droplet breakup under post-shock conditions, highlighting the robustness and applicability of the proposed hybrid Eulerian–Lagrangian formulation for such advanced applications in fluid engineering.

1. Introduction

Droplet aerobreakup, also known as secondary atomization, refers to the disintegration of liquid droplets induced by high-speed gas flows. This phenomenon poses a complex challenge in multiphase flow investigations owing to its highly transient characteristics spanning both spatial and temporal dimensions. Atomization driven by aerodynamic forces is critically important in various engineering and industrial applications, including liquid jet atomization, agricultural spraying, energy systems, and nuclear fusion processes [1,2,3,4]. Also, the breakup of liquid droplets in high-pressure airflows is particularly important during compression strokes in internal combustion engines and compressors in gas turbines [5]. Additionally, aerodynamic breakup triggered by passing shock waves is a fundamental aerospace engineering concern, with advanced applications ranging from fuel injection in rocket engines to erosion damage in supersonic flights [6]. Understanding droplet deformation and breakup is especially important in applications where increased airflow pressure adds complexity to droplet behavior. Thus, gaining a deeper understanding of droplet deformation and breakup is crucial under both atmospheric and higher-pressure conditions.

In the last decades, extensive experimental and numerical research has been carried out to investigate droplet aerobreakup. Various classifications of droplet breakup regimes are documented in the literature, defined by critical dimensionless parameters. The latter include the Weber number $\left(We\right)$, which characterizes the predominance of aerodynamic (inertial) forces over surface tension (restorative) forces in free-stream conditions, and the Ohnesorge number $\left(Oh\right)$, which accounts for the influence of liquid viscosity on breakup dynamics. These two parameters are defined as follows:

$$We={\displaystyle \frac{{\rho }_g{u}_g^2{D}_0}{\sigma }},$$

$$Oh={\displaystyle \frac{{\mu }_l}{\sqrt{{\rho }_l\sigma D_0}}},$$

depending on factors such as initial droplet diameter ($D_0$), surface tension coefficient ($\sigma$), post-shock gas density (${\rho }_g$) and velocity ($u_g$), liquid density (${\rho }_l$), and dynamic viscosity (${\mu }_l$). For instance, Joseph et al. [7] documented the breakup of droplets through the bag and stamen process at elevated Weber numbers for both viscoelastic and Newtonian fluids. Importantly, as $Oh$ increases, a delay in the deformation and breakup process is observed, with the viscosity resulting in having minimal influence for $Oh<0.1$ [8].

Several comprehensive reviews, such as those by Pilch and Erdman [9] and Theofanous [10], have individuated various modes of droplet fragmentation, along with their associated criteria for regime transitions. Traditional works have identified five aerobreakup modes, characterized by distinct ranges at an increasing Weber number, namely, bag breakup, bag-stamen/plume breakup, multibag breakup, sheet thinning/stripping breakup, and catastrophic breakup [11,12,13]. However, some breakup investigations conducted for air-jet flows, at moderately low $We$ values, have revealed multimode structural formations due to Rayleigh–Taylor instability (RTI). Furthermore, the transition from other modes to shear breakup regime, distinguished by equatorial liquid sheet formation, was demonstrated to be influenced by viscosity and the density ratios between the two phases [8]. Therefore, also based on these studies, various breakup modes have been reclassified taking into account the hydrodynamic instabilities responsible for the physical phenomenon [10]. Specifically, for low Weber numbers ($We<{10}^2$), the previously identified bag, bag-stamen, and multibag breakup modes have been consolidated into the single Rayleigh–Taylor piercing (RTP) mode, attributed to RTI. Moreover, at higher Weber numbers ($We>{10}^3$), the previously identified sheet thinning/stripping and catastrophic modes have been combined into the shear-induced entrainment (SIE) regime, associated with Kelvin–Helmholtz (KH) instability.

In fact, the above-mentioned hydrodynamic instabilities are crucial for determining droplet morphology during the aerodynamic atomization process. RTI typically drives breakup modes involving bag formation, with various types of bag breakups observed. For instance, recent studies by Theofanous and co-workers [10,14] have utilized finite-thickness RTI analysis to identify the conditions leading to initial criticality for the bag formation process. In contrast, shear-induced droplet breakup is driven by KH instability on the liquid droplet surface [14,15]. Several studies have explored the complex process of liquid droplet breakup using the shock tube technique. For instance, Sembian et al. [16] combined experimental and numerical methods to investigate the dynamics of reflected, transmitted, and diffracted waves during the interaction of an incident shock wave with a cylindrical water column [17]. Their findings indicated that, at higher shock Mach numbers, the focusing of expansion waves within the droplet could also lead to cavitation. Moreover, a number of numerical simulations have demonstrated the formation of a recirculating flow in the equatorial region, together with an upstream jet flow in the droplet’s wake, which promotes droplet deformation and breakup [18,19,20].

The rapid aerodynamic fragmentation of liquid droplets, which is particularly pronounced at high Weber numbers, poses significant challenges for numerical simulations due to its complex multiscale nature and short characteristic time scales. Achieving high spatio-temporal resolutions would be essential to observe the transformation of millimeter-sized primary droplets into micron-sized secondary droplets within these brief time intervals. While increasing computational resources can enhance numerical simulations, comprehensive studies that effectively align with experimental data across a range of parameters, and fully capture the evolving dynamics of this phenomenon, are still lacking. Indeed, the limitations and challenges associated with current numerical methodologies for simulating droplet aerobreakup highlight a significant gap in the literature regarding a detailed analysis of dispersed mist development under the SIE regime. This gap primarily arises from the prevalence of multiscale structures and the substantial computational costs required for thorough analyses. Despite extensive research, the exploration of these fundamental mechanisms remains limited, and comprehensive three-dimensional numerical investigations under high-pressure conditions are relatively rare. Direct numerical simulation (DNS), with a full resolution of the phase interface, would be the only universal tool. However, due to the associated immense computational cost, DNS can be applied only to a very limited class of simple model problems, while being unaffordable for applications of engineering interest. Moreover, some limitations of the commonly employed volume of fluid (VOF) method are well-documented, particularly in resolving fine-scale features and accurately capturing the interface dynamics of sub-droplets. These limitations include interface smearing, resolution dependency, and a restricted particle tracking capability.

In this complex scenario, the main aim of this study was to test a hybrid Eulerian–Lagrangian approach, combining VOF with the discrete phase model (DPM), e.g., [21]. This method exploits the VOF technique’s strength in capturing initial breakup dynamics, together with the DPM’s ability to accurately track individual sub-droplets. The proposed computational model makes use of the large-eddy simulation (LES) approach, where the large turbulent structures are resolved, while modeling the effect of unresolved small turbulent eddies. The LES-based hybrid VOF-DPM method allows us to explore the multiscale characteristics of droplet aerobreakup at a feasible computational cost. In fact, the present method effectively predicts both the mechanisms of early-stage breakup and the subsequent dispersion of fragments generated by shock waves at different Mach numbers, as demonstrated by comparisons against experimental observations. Due to the DPM approach, a key contribution of this study is the comprehensive quantitative analysis of droplet fragmentation dynamics and mist formation.

The rest of the manuscript is structured as follows. Section 2 presents the Eulerian–Lagrangian approach that is proposed, while the corresponding applied CFD model is introduced in Section 3. Numerical results for two different flow configurations at different Weber numbers are presented and discussed in Section 4. Finally, Section 5 summarizes the key conclusions of this research.

2. Hybrid Eulerian–Lagrangian Formulation

In this section, the three building blocks of the proposed hybrid Eulerian–Lagrangian approach are briefly presented. The interested reader is referred to the mentioned references for further details.

2.1. VOF Method

The VOF method [22] represents the most widely utilized technique among algorithms relying on continuous volume markers. Its popularity derives from the simple implementation and high efficiency, while it has demonstrated excellent performance in calculating interfacial flows. The core concept involves treating the liquid and gas as a single two-phase medium, where a special marker function $F(x,y,z,t)$ defines the space–time distribution of the phases within the computational domain. The volume fraction of the liquid phase within a computational cell is characterized as follows: $F=0$ indicates a cell without liquid, $F=1$ denotes a cell entirely filled with liquid, and $0<F<1$ represents a cell containing the air–liquid interface. The movement of the interface is tracked by solving the following transport equation for the liquid volume fraction, as the interface moves along with the liquid:

$${\displaystyle \frac{\partial F}{\partial t}}+\nabla ·\left(F\mathbf{V}\right)=0,$$

(1)

where $\mathbf{V}$ is the velocity vector of the two-phase medium. The latter is determined by solving a system of compressible flow governing equations, including the mass conservation equation

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathbf{V}\right)=0,$$

(2)

and the momentum balance equation

$${\displaystyle \frac{\partial \left(\rho \mathbf{V}\right)}{\partial t}}+\nabla ·\left(\rho \mathbf{V}\mathbf{V}\right)=-\nabla p+\nabla ·{\sigma }^v+{\mathbf{F}}^S.$$

(3)

Here, $\rho$ is the density of the two-phase medium, p is the static pressure, ${\sigma }^v$ is the viscous stress tensor, and ${\mathbf{F}}^S$ is the volumetric force vector field. The components of the viscous stress tensor are expressed as follows:

$${\sigma }_{ij}^v=\mu \left({\displaystyle \frac{\partial V_i}{\partial x_j}}+{\displaystyle \frac{\partial V_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\mu {\displaystyle \frac{\partial V_k}{\partial x_k}}{\delta }_{ij},$$

(4)

where $\mu$ represents the dynamic viscosity of the two-phase medium. The local density and molecular viscosity of the mixture fluid are determined using the volume fraction of the liquid in the cell, according to the mixture rule, namely:

$$\rho =F{\rho }_l+(1-F){\rho }_g,$$

(5)

$$\mu =F{\mu }_l+(1-F){\mu }_g,$$

(6)

where ${\rho }_l$ and ${\mu }_1$ are the density and viscosity of the liquid, and ${\rho }_g$ and ${\mu }_g$ are the density and viscosity of the gas, respectively.

When analyzing fluid flows with a gas–liquid interface, special consideration has to be given to the phenomenon of surface tension. In fact, studying flows affected by surface tension forces represents a complex and distinct task. One of the benefits of the VOF method is in the ability to facilitate the relatively straightforward inclusion of surface tension effects. The typical continuous surface force (CSF) algorithm proposed by Brackbill et al. [23] is employed here to simulate surface tension within the current VOF framework. This algorithm introduces the additional volumetric force ${\mathbf{F}}^S$ into the momentum equation, which is determined by the following relation:

$${\mathbf{F}}^S=\sigma \kappa \nabla F,$$

(7)

where $\sigma$ is the air–liquid surface tension coefficient, and $\kappa$ is the curvature of the free surface, defined as the divergence of the normal unit vector, namely:

$$\kappa =\nabla ·{\displaystyle \frac{\mathbf{n}}{\left|\mathbf{n}\right|}}.$$

(8)

Finally, the energy conservation equation reads

$${\displaystyle \frac{\partial \left(\rho E\right)}{\partial t}}+\nabla ·\left(\rho \mathbf{V}(E+{\displaystyle \frac{p}{\rho }})\right)=\nabla ·(\lambda \nabla T),$$

(9)

where E represents the total energy, which is the sum of the internal and kinetic energies. To account for the thermal effects, the thermal conductivity and heat capacity of the two-phase medium are defined using the mixing rule:

$$\lambda =F{\lambda }_l+(1-F){\lambda }_g,C_p=F{C}_{pl}+(1-F)C_{pg},$$

(10)

where ${\lambda }_l$ and ${\lambda }_g$ are the thermal conductivities of the liquid and gas phases, respectively, and $C_{pl}$ and $C_{pg}$ are the corresponding heat capacities.

2.2. LES Approach

As far as fluid turbulence is concerned, in this work, the LES approach is followed, where part of the turbulent scales are directly resolved, while modeling the effect of the residual part. Specifically, the spatially filtered Navier–Stokes equations are numerically solved [24], which allows us to accurately describe the space–time evolution of the mixture flow with affordable computational complexity. In contrast, the Reynolds-averaged Navier–Stokes (RANS) models provide only averaged flow characteristics and do not adequately capture the development of instabilities in droplet–shock wave interactions, e.g., [25]. Note that classical numerical estimates for this particular two-phase flow problem indicate that DNS calculations would require computational grids with hundreds of millions of cells for typical configurations, while LES represents a viable numerical technique for aerobreakup studies, e.g., [26,27].

The system of compressible LES governing equations can be expressed as follows:

$${\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial x_i}}=0,$$

(11)

$${\displaystyle \frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho {\overline{u}}_i{\overline{u}}_j\right)=-{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+{\displaystyle \frac{\partial {\overline{\sigma }}_{ij}}{\partial x_j}}+{\displaystyle \frac{\partial {\tau }_{ij}}{\partial x_j}}+F_i^S,$$

(12)

where ${\overline{\sigma }}_{ij}$ is the filtered viscous stress tensor, which is defined in terms of resolved velocity gradients by

$${\overline{\sigma }}_{ij}=\mu \left({\displaystyle \frac{\partial {\overline{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\overline{u}}_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\mu {\displaystyle \frac{\partial {\overline{u}}_k}{\partial x_k}}{\delta }_{ij}.$$

(13)

The tensor components ${\tau }_{ij}$, which are referred to as subgrid scale (SGS) stresses, represent the effect of unresolved turbulent fluctuations upon the dynamics of resolved turbulent eddies. Analogous to RANS models, these unknown quantities can be determined using the Boussinesq approximation:

$${\tau }_{ij}-{\displaystyle \frac{1}{3}}{\tau }_{kk}{\delta }_{ij}=-2{\mu }_t{\overline{S}}_{ij},$$

(14)

where ${\overline{S}}_{ij}$ is the resolved strain-rate tensor, namely,

$${\overline{S}}_{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial {\overline{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\overline{u}}_j}{\partial x_i}}\right).$$

(15)

Following the classical Smagorinsky model [28], the SGS viscosity parameter is provided as ${\mu }_t=\rho L_S^2\left|\overline{S}\right|$, where $L_S$ denotes the mixing length of the (unresolved) subgrid scales, namely:

$$L_S=min(\kappa d,C_S{V}_{\mathrm{cell}}^{1/3}),$$

(16)

where $\kappa$ is the von Karman constant, d the distance to the nearest wall, $V_{\mathrm{cell}}$ the volume of the computational cell, and $C_S$ the Smagorinsky coefficient. In this study, employing the dynamic formulation of this model, the latter parameter is not constant but evaluated based on the localized instantaneous flow conditions, improving the accuracy of the LES solution [29,30]. The dynamic model equations, which are not reported here for brevity, can be found, for instance, in [31], where an assessment of various SGS models, in a similar DPM-LES framework, was performed for a different fluid engineering application, namely, water-reactive fire suppression.

2.3. DPM Method

The present research employs the hybrid VOF-DPM approach, with VOF capturing the dynamics and deformation of primary droplets and DPM tracking the formation and trajectories of small liquid fragments. The method detects the formation of small sub-droplets by monitoring the breakup of large droplets. This is achieved through the VOF-LES method, which tracks the interface between the water bodies and surrounding air. Practically, when the primary droplets break up, the VOF technique identifies the regions where new sub-droplets are formed. These regions are then treated as discrete particles by the DPM methodology. This way, the DPM algorithm can track the trajectories of the above water fragments by solving the corresponding Lagrangian equations of motion, where the momentum balance considers the various forces acting on the particles, including aerodynamic loads and weight. By solving these equations, the method calculates the position and velocity of each identified sub-droplet at each discrete time step, allowing for a detailed analysis of their movement and interaction within the airflow. Importantly, the DPM algorithm also accounts for the effects of resolved turbulence and recirculation zones, which influence the trajectories of sub-droplets, their dispersion, and further breakup. As the method continuously updates the position and velocity of each sub-droplet, throughout the simulation, the dynamic tracking provides real-time insights into the dispersion patterns and interactions of the particles with the surrounding airflow, while characterizing their statistical distribution, in terms of both size and velocity.

3. Computational Model

In this section, the details of the overall applied CFD model for the droplet breakup simulation are provided, including the flow geometry and the main numerical settings.

3.1. Flow Geometry

The compressible flow governing equations were solved within a rectangular computational domain defined as $\mathrm{\Omega }=[-80D_0,120D_0]\times [-30D_0,30D_0]\times [-30D_0,30D_0]$, where $D_0$ represents the initial diameter of the spherical droplet. The domain size was determined based on prior studies to minimize boundary condition effects while maintaining affordable computational complexity. Practically, a sufficiently large computational domain was chosen to ensure that no significant reflections occur at the boundaries (sides, inlet, and outlet) within the time interval of interest.

Differently from analogous aerobreakup studies, e.g., [16,27], fully three-dimensional simulations were performed without imposing any symmetries or simplifications to the computed flow field, which allows us to accurately capture the non-symmetric complex behavior of the deforming liquid droplet surface and the surrounding airflow. The simplified flow geometry is illustrated in Figure 1, along with a schematic of the physical model, including the traveling normal shock front, before the impact on the droplet. The reference coordinate system $(x,y,z)\equiv (x_1,x_2,x_3)$ is oriented with the first axis in the stream-wise direction, with the origin being located at the leading edge of the spherical water droplet in its initial position. Regarding boundary conditions, pressure inlet and pressure outlet conditions were imposed in the streamwise direction, while the four lateral sides were treated as symmetry boundaries, thus limiting artificial reflections of waves or disturbances, which could otherwise affect the simulation accuracy.

In contrast to previous studies [21,32], the flow within the shock tube was not explicitly simulated. Rather, the discontinuous airflow conditions across the moving shock front were directly imposed. This simplified simulation setup expedited the calculations, where the reduced computational complexity allowed us to allocate more resources to resolving the droplet dynamics, thus focusing on the primary research interest. However, it should be mentioned that such a procedure can lead to losing detailed information about the ambient airflow, while not exactly reproducing the initial conditions for the breakup process, which potentially implies discrepancies between simulation results and experimental observations [25,33].

The simulations began with the shock front located one diameter away from the droplet’s leading edge ($x/D_0=-1$), thus separating the domain into post-shock and pre-shock sections. It is important to note that the shock front, and consequently the post-shock airflow, moves in the positive x-axis direction.

3.2. Case Studies

The simulations of the water droplet aerobreakup were carried out for two distinct Weber numbers: 205 and 7000. Following previous studies [21,34], the physical properties used in the calculations were prescribed as follows: a water density of 998 ${\mathrm{kg}/\mathrm{m}}^3$, a water viscosity of $1.003\times {10}^{-3}$ $\mathrm{kg}/(\mathrm{m}·\mathrm{s})$, and a surface tension coefficient of 0.073 N/m.

The corresponding Mach numbers of the propagating shock waves, the post-shock flow conditions, and the Weber and Reynolds numbers calculated using post-shock air properties are presented in

Table 1. Shock wave and post-shock conditions for the two different aerobreakup configurations.

Case	${\mathit{D}}_0$ (mm)	${\mathit{M}}_{\mathit{s}}$	${\mathit{u}}_{\mathit{g}}$ (m/s)	${\mathit{\rho }}_{\mathit{g}}$ (kg/m3)	$\mathit{Re}$	$\mathit{We}$
1	2.4	1.12	65	1.44	$1.2\times {10}^4$	205
2	4.8	1.47	226	2.17	$1.07\times {10}^5$	$7\times {10}^3$

. These quantities can be used to non-dimensionalize the various flow variables, following the usual methodology, e.g., [18], facilitating a direct comparison across different flow cases. In particular, the non-dimensional time variable is defined as follows:

$$t^{\ast }={\displaystyle \frac{u_{\mathrm{g}}}{D_0\sqrt{\epsilon }}}(t-t_0),$$

(17)

where $\epsilon$ stands for the ratio between liquid and gas densities, with $t_0$ representing the time instant of the shock impact on the droplet.

Previous similar research works based on no-model calculations [35,36] employed a computational grid with a spatial resolution exceeding 1000 cells per initial diameter. This resolution was necessary to accurately resolve the viscous boundary layer and predict the occurrence of KH instabilities. However, these studies were limited by their substantial computational demands, focusing primarily on demonstrating early-stage instabilities occurring on the cohesive surface of droplets. These studies did not extend the analysis to investigate the subsequent stages involving fragmentation and mist formation, which constitute the principal focus of the present research. Here, due to the adopted LES approach, which utilizes an SGS model to approximate the effects of unresolved turbulent eddies, a spatial resolution of 100 cells per original diameter was employed in the flow region of primary interest. This resolution was also chosen based on similar works, where it was proven to be effective in accurately capturing the overall deformation of the droplet surface while providing an approximate estimation of KH hydrodynamic instabilities [37,38]. According to previous research findings at comparable Mach numbers, e.g., [27], the current flow configurations are characterized by a Kolmogorov length scale of approximately $0.5\mathsf{\mu }\mathrm{m}$, for which the DNS approach would be practically unaffordable.

3.3. Simulation Setup

The commercial CFD solver ANSYS Fluent 23R2 was employed for practical calculations, where the VOF model is available with the pressure-based solver, with the two-phase flow governing equations being numerically solved using the finite volume method, e.g., [39]. The rectilinear mesh was made of cubic cells, with an automatic mesh refinement at the air–water interface.

Regarding spatial discretization, a second-order central difference scheme was used to approximate the convective terms in the filtered momentum equation. The volume fraction equation was solved using an explicit method employing the geo-reconstruct discretization scheme [36,40]. An implicit first-order method was utilized to approximate the unsteady terms in the LES governing equations. Diffusion fluxes and source terms were approximated with a second-order accuracy. Also, second-order upwind discretization was used for the continuity and energy equations. The coupling between the velocity and pressure fields was established using the SIMPLEC (semi-implicit method for pressure linked equations-consistent) procedure, which successfully addresses the challenges associated with resolving the transient phase interface, e.g., [25]. The resulting system of discretized equations was solved iteratively using a multigrid solver. Staggered grids were employed alongside the PRESTO! discretization scheme for pressure, which was chosen for its ability to yield more accurate results by mitigating interpolation errors, while avoiding assumptions about pressure gradients at boundaries. This scheme is particularly effective for problems involving strong body forces, such as surface tension, and high-density ratios [26].

Regarding temporal integration, the varying time step size was controlled using a specified maximum value for the Courant–Friedrichs–Lewy (CFL) parameter, which was 0.25. This low value was also selected based on previous similar studies, e.g., [41], to ensure numerical stability and accuracy, practically balancing computational efficiency with the need to accurately capture transient phenomena. Correspondingly, typical time step values ranged from $1\times {10}^{-8}$ to $4\times {10}^{-8}$ s depending on the flow conditions.

4. Results and Discussion

In this section, the results of the different simulations are presented and discussed. Comprehensive qualitative and quantitative analyses of droplet drift, deformation, and breakup are conducted while making comparisons with reference data from the relevant literature.

4.1. Shock–Droplet Interaction

When the shock front encounters the droplet surface, a range of complex flow structures is produced, including reflected, transmitted, and diffracted waves, profoundly altering the local flow conditions around and inside the droplet. While the complex internal and external flow fields are challenging to capture experimentally, recent advances in numerical methods and availability of suitable computational power allow for their observation in numerical simulations. Indeed, as demonstrated by the following results, the proposed LES model is able to accurately reproduce the shock–droplet interaction mechanism and subsequent breakup.

The compressible flow dynamics associated with the computed shock–droplet interaction is illustrated in Figure 2 for the configuration at $M_s=1.47$. The type of reflected wave (either shock or expansion) actually depends on the relative acoustic impedance of the media on both sides of the boundary [42]. In this case, since $Z_{\mathrm{water}}>Z_{\mathrm{air}}$, both reflected and transmitted waves result in being shocks (panel (a)). The transmitted wave creates a high-pressure region behind it, driving fluid motion towards the droplet’s downstream side [19]. Given the different impedance values, pressure waves travel faster in water than in air, causing the transmitted wave to detach from the incident wave, and first reach the leeward side of the droplet [16,20]. The diffraction of the Mach stem occurs after it traverses the droplet’s equatorial plane, contributing to the delay of the Mach wave near its contact point with the droplet’s surface. This Mach wave forms axisymmetrically around the droplet and converges near the rear stagnation point, intensifying the shock and creating a localized high-pressure region (panel (b)). The internal flow remains directed towards the leeward side of the droplet until the shock converges at the rear stagnation point. The high-pressure forces induce localized deformation of the droplet surface, leading to an upstream flow within the droplet. The two opposing flows interact at a saddle point inside the droplet, redirecting the flows towards the droplet’s equator. Additionally, as the external airflow evolves around the droplet’s periphery, a toroidal (primary) vortex forms near the flow separation point, similar to the flow around a solid sphere (panel (c)). Over time, this vortex structure migrates downstream, creating a recirculation zone and generating a jet flow at the rear stagnation region (panel (d)). Hence, the interaction between the incident shock wave and the droplet results in a combined effect where the recirculated jet flow, along with the expansion of the high-pressure region, causes the leeward side of the droplet to flatten, while taking a cupcake-like shape. Overall, the observed wave dynamics during the interaction process align with previous experimental studies [16,43], as well as various numerical simulations of shock-induced vaporization of droplets [19,20].

4.2. Verification and Validation

To validate the present LES model, the simulation results are compared with experimental data by Poplavski et al. [26], experimental images provided by Theofanous et al. [14], and numerical results by Meng and Colonius [41]. Figure 3 shows some side views of the deforming droplet in the meridian plane, plotted against non-dimensional time. An analysis of snapshots within the time frame $0.074\le t^{\ast }\le 0.49$ indicates that the numerical results for the configuration at $M_s=1.47$ (as detailed in Table 1), correspond well to the reference experiment examining the breakup of a water droplet at $We=780$ for a post-shock flow Mach number of 0.32 [14]. Generally, the qualitative agreement with experimental photographs is evident, both in terms of droplet shape and deformation dynamics. However, at the last time point, a slight difference can be observed, where the numerical solutions show a cupcake-like shape, while the experimental shape looks like an inverted bowl, which is likely due to the different Weber number.

Beyond the previous qualitative analysis, a key quantitative indicator of droplet deformation is represented by the time increasing cross-sectional extent of the droplet, which also affects the droplet dynamics. Figure 4 illustrates the relative growth of the droplet’s midsection, represented as $D_{cro}/D_0$, which is the ratio of the current transverse droplet dimension to its initial size, for $M_s=1.12$ (see Table 1). The figure indicates a strong correlation between the average deformation rate obtained from an experiment [26] and the current LES data. In fact, despite minor discrepancies, the change in the droplet’s cross-stream deformation generally aligns acceptably well with the experimental findings in terms of both trend and magnitude. Though the reference data seem slightly underpredicted, the proposed VOF-LES method is demonstrated to effectively capture the droplet deformation and initial breakup characteristics, even considering that the experimental uncertainty is not known.

Basically, the primary objective of numerical simulations focusing on the early stages of shock–droplet interactions is to accurately predict the motion and deformation of the water body, which is critical for understanding the subsequent breakup phase. The droplet’s kinematic characteristics are usually assessed by analyzing the streamwise position and velocity of its center-of-mass (CM). These time-dependent parameters are defined, according to the adopted VOF formulation, as follows:

$$x_{CM}\left(t\right)={\displaystyle \frac{{\int }_{\mathrm{\Omega }}F{\rho }_{l}xd\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}F{\rho }_{l}d\mathrm{\Omega }}},$$

(18)

$$u_{CM}\left(t\right)={\displaystyle \frac{{\int }_{\mathrm{\Omega }}F{\rho }_{l}ud\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}F{\rho }_{l}d\mathrm{\Omega }}}.$$

(19)

The above integrals, while spanning the entire computational domain $\mathrm{\Omega }$, exclusively account for contributions from space regions where the liquid phase is present, as determined by a non-zero volume fraction F. The displacement and velocity parameters are normalized here as $\mathrm{\Delta }{x}^{\ast }={\displaystyle \frac{x_{CM}\left(t\right)-x_{CM}\left(0\right)}{D_0}}$ and $u^{\ast }={\displaystyle \frac{u_{CM}\left(t\right)}{u_g}}$, respectively. Figure 5 shows the time histories of these two parameters at $M_s=1.47$. The present LES results demonstrate a good correlation with the reference DNS data [41]. However, while these transient parameters are captured with a satisfactory accuracy, a slight underestimation is observed at later time instants. Apparently, the error corresponding to the approximation of unresolved flow scales in the LES modeling approach increases with time.

4.3. Droplet Breakup Dynamics

The flow patterns for the two test cases under study are analyzed in Figure 6, showing the velocity magnitude corresponding to the (x, y) plane at various time instants. For the first configuration at $We=205$, the breakup regime falls into the transition zone between RTI- and shear stripping-based mechanisms. This region is characterized by the simultaneous development of multiple breakup mechanisms, leading to what is referred to as a mixed or multimode regime [44]. At low $We$ numbers, where the droplet experiences minimal deformation, the surrounding air flow resembles that corresponding to a solid sphere at Reynolds numbers between ${10}^3$ and ${10}^4$, characterized by a vortex separation near the midsection, as seen at $t^{\ast }=0.03$. Seemingly, at later time instants ($t^{\ast }=0.06$), the wake flow is characterized by two vortical zones, which appear enlarged at $t^{\ast }=0.14$, confirming previous results by Boiko and Poplavski [45]. These characteristic structures remain discernible until approximately $t^{\ast }=0.33$, when a toroidal vortex and reverse flow develop, accompanied by the emergence of two annular waves on the droplet’s surface. Then, the wake flow loses its axial symmetry and breaks down into smaller, non-stationary structures ($t^{\ast }=0.54$). These structures form a large stagnant region with a low average velocity, similar to the velocity in the flow deceleration zone in front of the droplet’s windward surface. This flow results in the surface flattening and causes the air to spread radially across the rear surface. This process leads to the formation of a liquid disk whose edge facilitates the transfer of droplet fragmentation products into the aerodynamic wake. Overall, this type of droplet breakup actually corresponds to the so-called “sheet stripping” mechanism.

The above breakup dynamics are characteristic across the range of $We$ numbers investigated. The analysis of the air velocity field reveals flow characteristics that align with experimental observations. Specifically, following [45], the external flow around the droplet is inherently unsteady, yet it comprises some stable gas-dynamic structures, each influencing the droplet in a unique manner. Initially, a prominent wave known as the generatrix bend forms. This wave emerges when the droplet is still spherical, at the separation point, where the velocity reaches its maximum. Beyond this point, a recirculation zone with a counter-current flow along the droplet’s surface becomes evident. Furthermore, in the wake, a secondary stable structure forms as a recirculation zone with an intense reverse flow along the longitudinal axis. This flow causes surface flattening, leading the air to spread radially across the rear surface. In this region, a second break in the generatrix occurs, which is linked to the radial spreading of liquid water. This process leads to the formation of a liquid disk whose edge facilitates the transfer of droplet fragmentation products into the aerodynamic wake. Additionally, the air flow separates near the second wave, radiating outward across the rear surface. Part of this flow re-enters the initial separation zone along the droplet’s surface, while the remainder transitions into a third stable structure, represented by an isolated toroidal vortex that partially separates the previous two. Overall, this type of droplet breakup corresponds to the so-called “sheet stripping” mechanism.

Furthermore, the interaction between shock front and water droplet induces a range of intricate flow phenomena within the airflow field, including the formation and development of important vortical structures. Figure 7 shows the vorticity field for $M_s=1.47$ at the same time instants of the previous figure. As the flow field evolution at $M_s=1.12$ results in being substantially similar, the higher Mach number case highlights the investigated effects more clearly. By analyzing the droplet morphology at early times (panel (a)), it is observed that a recirculation zone on the windward side begins to form on the droplet’s surface. As time advances, the recirculation zone expands and slowly transitions from the center towards the equator, coinciding with the complete passage of the shock wave through the droplet (panel (b)). Importantly, the larger the windward recirculation zone, the greater the shearing force exerted by the airflow. The expanding windward recirculation zone leads to the formation of KH waves on the windward surface, similar to those described by Jalaal and Mehravaran [46]. At $t^{\ast }=0.14$, flow separation leads to the formation of the wake recirculation zone near the leeward stagnation point (panel (c)). As flow separation intensifies, the wake recirculation zone becomes more pronounced over time (panel (d)), following the shock–droplet interaction, eventually resulting in a flattened shape of the droplet on the leeward side (panel (e)). Basically, over time, the KH waves continue to grow, merge, and evolve into sheets or ligaments near the equator. The stretching of these liquid structures, along with the stripping of KH waves, contributes to the overall breakup process. The interaction between the sheets and wake recirculation zones on the leeward side of the droplet produces child droplets, primarily stripped from the radial edges of the parent droplet.

4.4. Sub-Droplets Distribution

Sub-droplets are formed when liquid sheets, ligaments, and lips detach from the parent droplet [9]. Practically, the initial breakup event occurs when the first ligament or sub-droplet separates from the droplet under shear breakup conditions. Due to the DPM methodology, in the framework of the employed CFD solver, the progression from initiation to the completion of the ligament fragmentation process can be tracked while counting the number of sub-droplets that are formed. Here, differently from previous research employing the unsteady RANS approach [21,25], the LES formulation allows for a more detailed analysis of sub-droplet size distributions while investigating the impact of varying the $We$ parameter.

Figure 8 demonstrates the generation of sub-droplets in terms of the percentage distribution of diameters at various Mach and Weber numbers, making a comparison against experimental data by Sharma et al. [47]. The distribution of droplet sizes tends to shift towards smaller values as $We$ increases, which represents a trend well documented in the literature. In fact, at higher $We$ values, the proportion of smaller droplets increases, whereas at lower $We$ values, larger droplets make up a greater portion of the total liquid mass. Note that sub-droplets generated by the direct stripping of KH waves result in being generally smaller than those formed from sheets, ligaments, and lips detaching from the edges and leeward side of the parent droplet. Initially, aerodynamic forces can effectively break apart larger sheets and ligaments into multiple larger sub-droplets. Over time, these larger fragments may further disintegrate into smaller sub-droplets. Notably, the present trend in particle size distribution aligns with the experimental data, confirming the accuracy of the LES approach. It is important to highlight that the reference experimental distribution actually corresponded to the temporal range $0.33<t^{\ast }<0.67$, whereas the preset numerical data corresponded to a given time instant: $t^{\ast }$ = 0.55. This difference represents a key factor contributing to the observed discrepancy in the diameter sizes. Nevertheless, the overall predicted log-normal distribution results in a consistency between the two different approaches, supporting the validity of the scale-resolving DPM-LES approach.

Furthermore, Figure 9 displays the relationship between particle velocities and associated diameters. Apparently, higher Weber and Mach numbers correspond to increased particle velocities. To examine the sub-droplets distribution, the Sauter mean diameter (SMD) can be computed. This latter parameter is defined as follows:

$$D_{32}={\displaystyle \frac{{\sum }_i{n}_i{d}_i^3}{{\sum }_i{n}_i{d}_i^2}},$$

(20)

where $n_i$ represents the number of droplets with a (discretized) diameter $d_i$. Practically, this expression calculates the mean diameter by weighting the volume and surface area of the tracked droplets. SMD is highly effective for characterizing the distribution of spray droplets, particularly in evaporation processes. The SMD values calculated throughout the current atomization process are reported in

Table 2. Overall SMD value as a function of Weber number.

We	SMD [$\mathsf{\mu }$m]
205	61.9
7000	15.0

, indicating a significant reduction in SMD as the Weber number increases.

5. Conclusions

This study investigated the aerodynamic breakup of water droplets induced by air shock waves traveling at different Mach numbers. The numerical simulations were carried out using the LES approach combined with a hybrid VOF-DPM technique. The obtained results confirm previous theoretical and experimental findings regarding the physical mechanisms driving droplet fragmentation.

In fact, the LES method, balancing computational efficiency with solution accuracy, effectively captures flow instabilities and the main turbulent flow structures during the shock–droplet interaction, offering key insights into water body deformation and breakup under high-speed airflow. At lower Weber numbers, the breakup transitioned between RTI and stripping regimes, indicating a multimodal regime. The flow around the droplet was similar to that of a sphere at comparable Reynolds number, with vortex shedding near the midsection. A toroidal vortex and reverse flow initially developed, generating annular waves that influenced the subsequent fragmentation. The sheet stripping mechanism was consistently observed across simulations, with KH waves becoming particularly pronounced at higher Mach number. In fact, two distinct stripping mechanisms affect the droplet surface evolution: the primary mechanism generates a main stream of fragments from the droplet’s equator to the flattened rear, while the secondary mechanism operates at the front. Also, larger fragments form in regions of intense stripping, then evolving into a fine mist. The droplet undergoes continuous stretching and compression, with more pronounced deformation at higher Mach numbers. The sub-droplet size distribution shifts toward smaller values with increasing Weber number, confirming the model’s accuracy compared to experimental observations.

Overall, the present hybrid Eulerian–Lagrangian approach proved effective in capturing the complexities of droplet aerobreakup, as demonstrated by both qualitative and quantitative analyses. Compared to previous RANS-based studies, the present LES method offers a more robust and detailed understanding of the complex dynamics involved in droplet breakup under post-shock conditions. However, future research employing different SGS models could provide further assessments of the proposed formulation, where the combination of LES with DPM and VOF techniques appears very promising. Potentially, the method may lead to more accurate predictive models for advanced engineering applications such as fuel atomization and spray dynamics.