Numerical Simulation of Shock Wave Propagation Through Multiple Raindrops

1. Introduction

Shock wave interacting with liquid droplets arise in a wide variety of environments, including transonic/supersonic flight through precipitation, detonation-droplet interaction, and high-speed propulsion systems. In particular, aircraft operating in rainy weather inevitably encounter raindrops and shock wave interactions.

Numerous works have been carried out to study the interactions between shock waves and droplets [1,2,3,4]. Experiments using shock tubes and high-speed imaging have revealed that the breakup behavior strongly depends on the Mach number and the pressure loading across the interface [5,6,7]. Despite these advances, experimental approaches remain limited in their ability to capture the detailed wave structure and droplets interaction. The characteristic timescales are extremely short, and the spatial scales are often very small. As a result, most experimental measurements focus on macroscopic quantities such as droplet shape and breakup patterns [3,8], while detailed flow information such as pressure distribution, wave propagation and local Mach number remains difficult to obtain. Consequently, numerical simulations have become an essential tool to investigate the interactions between shock waves and droplets.

Numerical simulations of compressible multiphase flows have advanced significantly over the past years [9,10,11,12,13,14,15,16,17,18,19]. Existing computational methods can be classified into two categories: diffuse-interface method and sharp-interface method. Diffuse-interface method represents the gas-liquid interface as a mixture region with finite thickness over several computational cells [20]. This method is robust and widely used in compressible multiphase flow simulations. However, the smeared representation of the interface often result in loss of geometric accuracy, particularly in flows involving strong shock waves and large density ratios. Under such conditions, accurately capturing interface shape becomes challenging.

Sharp-interface method treats the interface as a discontinuity and explicitly reconstruct its geometry [21]. Approaches such as level-set [22,23,24], coupled volume-of-fluid (VoF) [25,26,27] and other reconstruction schemes [28] have demonstrated improved accuracy in representing interface dynamics. However, these methods were originally designed for incompressible or weakly compressible flows [21,29,30,31,32,33,34,35], often relying on pressure-based strategy [36,37,38,39,40]. When applied to strongly compressible flows involving shock waves, the complex coupling procedures and iterative pressure correction steps will compromise computational efficiency and stability. As a result, accurately resolving the interaction between shock waves and deformable liquid droplets remains challenging for existing numerical methods.

Another limitation of current studies is that most investigations consider the interaction between a shock wave and a single droplet, whereas realistic rain environments contain large numbers of raindrops. In such cases, raindrops are not isolated but exist in clusters or arrays, leading to strong aerodynamic interactions. When a shock wave propagates through this field, complex wave patterns may develop due to shock diffraction, reflection, and interference. Although several studies have investigated droplet clusters or clouds [41,42,43], the detailed mechanisms governing shock wave and multiple raindrops interaction remain rarely studied.

Motivated by these challenges, the present study performs a numerical study of the interaction between a shock wave and multiple raindrops. Density based strategy is chosen due to its superiority over pressure based method for compressible problems. Piecewise linear reconstruction (PLIC) method is used to reconstruct the geometry of interfaces. To combine these two methods, a space-time flux-integral approach is designed. First, the numerical method is validated using canonical liquid-gas shock wave problems and shock wave interface interactions. Second, the interactions between a shock wave and a single droplet and two droplets examined, respectively. Finally, shock wave and multiple raindrops interaction is studied. Through detailed flow field analysis, the present work aims to provide new insights into the physics of shock wave and multiple raindrops interaction.

2. Methodology

2.1. Governing Equations

The five-equation model proposed by Grégoire Allaire et al. [44], also known as the velocity-pressure equilibrium model, assumes instantaneous equilibration of velocity and pressure among all phases within each computational cell. In the regime of strong shock wave raindrops interactions, inertial and compressibility effects dominate the flow dynamics, rendering viscous dissipation, surface tension, gravity, and thermal conduction negligible. Consiquently, these effects are omitted, and the governing equations are expressed as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial {\alpha }_1}{\partial t}}+\mathbf{u}·\nabla {\alpha }_1=0,\hfill \\ {\displaystyle \frac{\partial {\alpha }_1{\rho }_1}{\partial t}}+\nabla ·\left({\alpha }_1{\rho }_1\mathbf{u}\right)=0,\hfill \\ {\displaystyle \frac{\partial {\alpha }_2{\rho }_2}{\partial t}}+\nabla ·\left({\alpha }_2{\rho }_2\mathbf{u}\right)=0,\hfill \\ {\displaystyle \frac{\partial \rho \mathbf{u}}{\partial t}}+\nabla ·\left(\rho \mathbf{u}\otimes \mathbf{u}+p\mathbf{I}\right)=0,\hfill \\ {\displaystyle \frac{\partial \rho E}{\partial t}}+\nabla ·\left[\left(\rho E+p\right)\mathbf{u}\right]=0,\hfill \end{array}\right.$$

(1)

where the subscripts $-1$$-2$ are for phase 1 and phase 2, ${\alpha }_k$, ${\rho }_k$ ($k=1,2$) denote the volume fraction and density of each phase, $\rho$ denotes the mixture density, i.e.,

$$\rho ={\displaystyle \sum _{k=1}^2}{\alpha }_k{\rho }_k,$$

(2)

with which the mass fraction of the k-th phase is defined as:

$$Y_k={\displaystyle \frac{{\alpha }_k{\rho }_k}{\rho }},$$

(3)

E denotes the mixture total energy:

$$E=e+{\displaystyle \frac{{\mathbf{u}}^2}{2}},$$

(4)

where e is the mixture specific internal energy:

$$e={\displaystyle \sum _{k=1}^2}{Y}_k{e}_k,$$

(5)

where $e_k$ is the specific internal energy of phase k, and p denotes the mixture pressure. The frozen speed of sound is employed

$$c=\sqrt{\sum _{k=1}^n{Y}_k{c}_k^2},$$

(6)

where $c_1$ and $c_2$ are the speeds of sound for each phase, which are determined by the pressure and density according to the equation of state.

The stiffened-gas (SG) EoS is adopted for both air and water. A constant "stiffened pressure" term, $p_{\infty }$ is introduced to account for the repulsive molecular forces that dominate in dense fluids. The fundamental form of the SG EoS is given by:

$$p_k=({\gamma }_k-1){\rho }_k{e}_k-{\gamma }_k{p}_{\infty ,k},$$

(7)

where ${\gamma }_k$ is the adiabatic index and $p_{\infty ,k}$ is the stiffened parameter. For air, $\gamma =1.4$ and $p_{\infty }=0$. For dense liquids such as water, the model requires significantly different parameters to capture their low compressibility. A representative calibration for water is:

$$\gamma =4.0,p_{\infty }=6.0\times {10}^8\mathrm{Pa},$$

(8)

or alternatively,

$$\gamma =1.9276,p_{\infty }=1.1373\times {10}^9\mathrm{Pa}.$$

(9)

With the stiffened gas (SG) equation of state (EoS) and the assumption of pressure equilibrium among phases, the mixture pressure can be computed:

$$p={\displaystyle \frac{\rho e-\left({\displaystyle \frac{{\alpha }_1{\gamma }_1{p}_{\infty ,1}}{{\gamma }_1-1}}+{\displaystyle \frac{{\alpha }_2{\gamma }_2{p}_{\infty ,2}}{{\gamma }_2-1}}\right)}{{\displaystyle \frac{{\alpha }_1}{{\gamma }_1-1}}+{\displaystyle \frac{{\alpha }_2}{{\gamma }_2-1}}}}.$$

(10)

Finally, the speed of sound is given by:

$$c_k=\sqrt{{\displaystyle \frac{p+{\gamma }_k{p}_{\infty ,k}}{{\rho }_k}}}.$$

(11)

2.2. Numerical Method

The second-order finite volume method with the computation of flux integral is adopted. VoF+PLIC method, which originally designed for incompressible flows, is adopted for strong compressibility in this work. The governing equations Eq. (1) in two-dimension can be written as:

$${\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+{\displaystyle \frac{\partial {\mathbf{F}}_x}{\partial x}}+{\displaystyle \frac{\partial {\mathbf{F}}_y}{\partial y}}=\mathbf{G},$$

(12)

where

$$\mathbf{U}=\left[\begin{array}{c}\hfill {\alpha }_1{\rho }_1\\ \hfill {\alpha }_2{\rho }_2\\ \hfill \rho u\\ \hfill \rho v\\ \hfill \rho E\\ \hfill {\alpha }_1\end{array}\right],{\mathbf{F}}_x=\left[\begin{array}{c}\hfill {\alpha }_1{\rho }_{1}u\\ \hfill {\alpha }_2{\rho }_{2}u\\ \hfill \rho u^2+p\\ \hfill \rho uv\\ \hfill (\rho E+p)u\\ \hfill {\alpha }_{1}u\end{array}\right],{\mathbf{F}}_y=\left[\begin{array}{c}\hfill {\alpha }_1{\rho }_{1}v\\ \hfill {\alpha }_2{\rho }_{2}v\\ \hfill \rho uv\\ \hfill \rho v^2+p\\ \hfill (\rho E+p)v\\ \hfill {\alpha }_{1}v\end{array}\right],\mathbf{G}=\left[\begin{array}{c}\hfill 0\\ \hfill 0\\ \hfill 0\\ \hfill 0\\ \hfill 0\\ \hfill {\alpha }_1\left({\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\end{array}\right].$$

(13)

Eq. (12) is integrated in space and time as:

$$\int \int \int {\displaystyle \frac{\partial \mathbf{U}}{\partial t}}dxdydt+\int \int \int {\displaystyle \frac{\partial {\mathbf{F}}_x}{\partial x}}dxdydt+\int \int \int {\displaystyle \frac{\partial {\mathbf{F}}_y}{\partial y}}dxdydt=\int \int \int \mathbf{G}dxdydt.$$

(14)

The shape of the interface is reconstructed by using PLIC method, as shown in Figure 1, where C is the cell averaged value of $\alpha$. After the reconstruction of the interface, the space-time flux integral must be considered. Otherwise, the numerical method will blow up because of the inconsistency between VoF+PLIC and compressible solver. The flux integral is computed according to Figure 2, where $\Delta t$ is the time step, and $\Delta t_{\Sigma }$ is the time for the interface to arrive at the integral point.

2.3. Raindrop Generation

The methodology for generating statistically consistent raindrop fields consists of spectrum modeling, probabilistic diameter sampling, and spatial placement with collision control. A maritime empirical spectrum from Figure 3 is fitted with a fifth-order polynomial with respect to the diameter D:

$${log}_{10}N\left(D\right)=A_0+A_{1}D+A_2{D}^2+A_3{D}^3+A_4{D}^4+A_5{D}^5,$$

(15)

where $N\left(D\right)$ represents the number concentration function. For stochastic sampling, the spectrum is normalized to construct a probability density function,

$$\mathrm{PDF}\left(D\right)={\displaystyle \frac{N\left(D\right)}{{\int }_{D_{\min }}^{D_{\max }}N\left(D\right)dD}},$$

(16)

where $D_{\min }$ is taken as $1.0\mathrm{mm}$ and $D_{\max }$ as $5.0\mathrm{mm}$ in this work. When the distribution is truncated, the PDF is re-normalized to preserve statistical consistency.

The diameters of raindrops are sampled using Monte Carlo method. The spatial distribution of raindrops is generated using a stratified sampling approach. The raindrops are generated within the domain of $[3\mathrm{mm}\le x\le 25\mathrm{mm}]$, $[0\le y\le 25\mathrm{mm}]$. The domain is divided into quasi-uniform cells according to the number of raindrops, and one raindrop is assigned to each quasi-cell with a randomized offset from the cell center. This procedure improves spatial uniformity while maintaining randomness. Geometric overlap is prevented through a collision detection algorithm based on pairwise distance evaluation. The raindrop centers are restricted to the ’padding zone’, ensuring full containment within the domain. The generation of 20 raindrops is given in Figure 4.

3. Results and Discussion

3.1. High-Pressure Water/Low-Pressure Air Shock Tube

A shock tube with an initial interface between high-pressure liquid water and low-pressure air is set up over a $1.0m$ domain, with the interface located at $x_0=0.7m$ (illustrated in Figure 5). The initial conditions, as described in [46], are as follows:

$$\left\{\begin{array}{cccccc}\hfill liquid:& z_1=1.0-{ϵ}_{cut}\hfill & u=0\hfill & p={10}^9\left[\mathrm{Pa}\right]\hfill & \rho =1000.0[\mathrm{kg}/{\mathrm{m}}^3]\hfill & forx<x_0,\hfill \\ \hfill gas:& z_1={ϵ}_{cut}\hfill & u=0\hfill & p={10}^5\left[\mathrm{Pa}\right]\hfill & \rho =50.0[\mathrm{kg}/{\mathrm{m}}^3]\hfill & forx\ge x_0,\hfill \end{array}\right.$$

(17)

where $x_0=0.7m$, and ${ϵ}_{cut}=1.0\times {10}^{-10}$. The heat capacity ratio of air is $\gamma =1.4$, and for water $\gamma =4.4$ and $P_c=6.0\times {10}^{8}Pa$.

The solutions were generated at $t=229\mu s$ with a CFL number of $0.1$. Table 1 shows the relative error of total mass of liquid using three numerical schemes (first-order, Minmod, and van Leer) over increasing mesh resolutions $N_{cell}$. Across all resolutions and limiters, the scheme conserves the liquid mass to $0.004-0.006\%$ relative error. Among the three schemes, van Leer gives the smallest relative mass error.

Figure 6, Figure 7, and Figure 8 show the comparison for density, velocity, and pressure, respectively, with the analytical solution using the van Leer limiter and four different $N_{cell}$ numbers: (a)1000, (b)2000, (c)4000, and (d)8000. This is a very challenging problem since the pressure ratio is $10,000$, and the shock wave moving to the air phase is very weak relative to the expansion fan in the water. As can be seen from Figure 6, there is a slight undershoot at the contact interface. However, the discontinuity is nearly as sharp as the analytical jump for $N_{cell}\ge 4000$. Key shock and contact discontinuities are well captured across all resolutions. Red cross markers (liquid and gas densities) show that both phase densities remain consistent and match the analytical interface values. Increasing the mesh resolutions from $N_{cell}=1000$ to $N_{cell}=8000$, the numerical solution converges cleanly to the analytical one.

3.2. Shock-Wave in Air/Helium-Bubble Interaction

The initial conditions for the interaction between a shock wave in the air and a helium bubble are based on [37]. Originally investigated experimentally by Haas and Sturtevant [47], this setup serves as a benchmark for exploring the governing physics and validating numerical algorithms. The helium bubble, with a radius of $2.5cm$, is positioned $5.0cm$ from the initial shock wave. Upstream of the shock wave, the pressure, velocity, and temperature are:

$$\begin{array}{ccc}{p}_{II}=1.01325\times {10}^5\left[Pa\right],& u_{II}=0.0\left[m/s\right],& T_{II}=351.82\left[K\right].\end{array}$$

(18)

Downstream of the shock wave, the pressure, velocity, and temperature are

$$\begin{array}{ccc}{p}_I=1.59060\times {10}^5\left[Pa\right],& u_I=125.65\left[m/s\right],& T_I=402.67\left[K\right].\end{array}$$

(19)

The relation between pressure, temperature, and density is:

$$p={\displaystyle \frac{\gamma -1}{\gamma }}\rho C_{p}T-P_c,$$

(20)

where $C_p$ is

$$C_p=\gamma C_v.$$

(21)

Figure 9 illustrates these initial conditions. The heat capacity ratio and the specific isochoric heat capacity of air are $\gamma =1.4$, and $C_v=720[J/(kg·K)]$, respectively. For helium, these values are $\gamma =1.648$ and $C_v=2440[J/(kg·K)]$.

The computational domain spans $0.0m\le x\le 0.2m$ and $0.0m\le y\le 0.2m$, with the helium bubble centered at $(x=0.1m,y=0.1m)$. The computational domain is discretized into $800\times 800$ cells. The shock wave reaches the helium bubble at $t_0=52\mu s$. Figure 10 and Figure 10 display the Mach number contours and numerical Schileren images at $t=100\mu s$, while Figure 11 and Figure 11 show these at $t=120\mu s$. Mach number contours and Schlieren images at $t=160\mu s$ and $t=200\mu s$ are provided in Figure 12 and Figure 13, respectively. In this test case, C represents the volume fraction of helium. Comparisons between the numerical results and experimental data from [47], presented in Figure 14, Figure 15, and Figure 16, show strong agreement. The computational (on the left) and the experimental (on the right) results both illustrate the deformation of the helium bubble over time due to the shock wave interaction. The initial bubble shape (dashed line) becomes increasingly distorted (solid line) as the shock wave passes through. The numerical Schlieren results closely match the experimental observations, capturing the deformation and wave interactions effectively. This agreement confirms the model’s capability in accurately representing the key physical mechanisms and the temporal evolution of the helium bubble’s shape. As time progresses from $t-t_0=52\mu s$ to $t-t_0=102\mu s$, the bubble undergoes significant elongation and distortion. The numerical model accurately replicates these deformations and the accompanying wave structures, demonstrating strong quantitative agreement and validating the effectiveness of the numerical method in depicting shock-bubble dynamics.

3.3. Mach 1.36 Shock Wave Interacting with One Raindrop

This test case investigates the interaction between a strong planar shock wave and a liquid raindrop to assess whether surface tension effects can be considered negligible under high compressible flow conditions. The initial condition is given in Figure 17 with a coming shock wave of $Mach=1.36$. Three conditions are compared with varying the Weber number: $We=\infty$, $We=3000$, and $We=300$, allowing evaluation of whether surface tension can be neglected without significantly affecting raindrop deformation dynamics and overall flow evolution.

Figure 18 and Figure 19 compare the pressure contours and numerical Schlieren fields at early ($t=1.0\mu s$) and later ($t=10.0\mu s$) times under different Weber numbers ($We=\infty$, $We=3000$, and $We=300$). At both instants, the overall shock structure, wave patterns, and pressure peak remain nearly identical across all Weber numbers considered. In particular, the incident shock strength, its reflection at the raindrop interface, and the downstream wake topology show negligible sensitivity to surface tension. Minor differences are limited to small-scale interfacial smoothing near the raindrop boundary, which does not affect the global shock dynamics or pressure distribution. These observations indicate that, for the strong shock regime examined in this work, the characteristic inertial and compressibility effects dominate the flow evolution, while capillary forces play a secondary role. Consequently, surface tension has a negligible influence on both the early shock impact process and the subsequent wave–raindrop interaction, justifying its omission in the governing equations for the present shock-dominated flow regime.

Figure 20 display the pressure contours at $1\mu s$, $2\mu s$, and $3\mu s$, with the interface marked by a red line, while the corresponding results from [48] shown in Figure 21. Similarly, the computed pressure contours at $4\mu s$, $6\mu s$, and $8\mu s$ are presented in Figure 22, with the results from [48] displayed in Figure 23. As time progresses (from $1\mu s$ to $8\mu s$), the shock wave dynamically interacts with the droplet. From $1\mu s$ (Figure 20,Figure 21) to $3\mu s$ (Figure 20,Figure 21), the shock wave deforms as it interacts and travels through the droplet. At later times ($4\mu s$, $6\mu s$, $8\mu s$), the pressure contour indicate further deformation of the droplet. The numerical results (Figure 20, Figure 22) align well with the results from the literature (Chang and Liou [48]), demonstrating that the model’s ability to produce expected results for droplet behavior under shock wave influence. Moreover, the present method enhances the accuracy of the sharp interface representation, as evident in Figure 22 compared to Figure 23. The progression of pressure contours across the entire computational domain is shown in Figure 24, while the evolution of the Numerical Schlieren is presented in Figure 25.

Figure 26 presents the total mass error of the droplet over time during the shock-wave and droplet interaction, comparing results for two grid resolutions: $400\times 400$ cells (black line) and $1000\times 1000$ cells (red line). Both simulations show mass errors within $0.15\%$, indicating stable behavior of the droplet’s total mass.

3.4. Mach 1.73 Shock Wave Interacting with Two Raindrops

Figure 27 illustrates the initial configuration for the interaction of a planar shock wave with two water raindrops in air. A planar shock wave of $Mach=1.73$ propagates in x-direction, corresponding to a nominal supersonic flow. Two initially stationary, circular water raindrops are embedded downstream of the shock, with radius of 3.20 mm and 2.50 mm, respectively. This configuration enables the investigation of shock interaction with two raindrops of different sizes, including the effects of wave transmission, reflection, and raindrop–raindrop interaction in a shock-dominated flow regime.

Figure 28, Figure 29, Figure 30, Figure 31, Figure 32, Figure 33, Figure 34, Figure 35 and Figure 36 show the temporal evolution of Mach number, pressure contours, and numerical Schlieren. Upon impact with the upstream raindrop (Figure 28), a strong pressure augmentation forms at the windward stagnation point, accompanied by a sharp Mach number reduction. The incident shock is reflected and diffracted around the raindrop, generating curved shock fronts and lateral expansion waves along the raindrop surface. As the interaction develops (Figure 29Figure 30), the transmitted shock subsequently impinges on the downstream raindrop, leading to secondary pressure amplification and additional reflected and diffracted wave structures. At later times (Figure 31, Figure 32, Figure 33, Figure 34, Figure 35 and Figure 36), multiple shock and expansion waves interact between the two raindrops. Both raindrops undergo progressive deformation driven primarily by shock-induced pressure gradients, with the upstream raindrop experiencing stronger flattening and lateral stretching, while the downstream raindrop deforms under the combined effects of transmitted shocks and wake-induced flow structures.

3.5. Mach 1.73 Shock Wave Interacting with Twenty Raindrops

To study the interference of multiple raindrops, the two raindrops in Section 3.4 are replaced by twenty raindrops. The methodology for generating statistically consistent raindrop fields consists of spectrum modeling, probabilistic diameter sampling, and spatial placement with collision control. A maritime empirical spectrum from Figure 3 is fitted with a fifth-order polynomial with respect to the diameter D:

$${log}_{10}N\left(D\right)=3.0697+0.6732D-1.2068D^2+0.4502D^3-0.0792D^4+0.0052D^5.$$

(22)

Figure 37 and Figure 38 present the temporal evolution of the Mach number and pressure fields. At early times ($t=1.0-3.0\mu s$), the incident planar shock is progressively distorted as it encounters individual droplets. As the shock wave penetrates deeper into the droplets, multiple reflected shocks and expansion waves emerge, leading to pronounced wave curvature and strong non-uniformity in both Mach number and pressure. The interaction between neighboring raindrops results in wave interference and cumulative pressure amplification. At later times ($t=6.0-12.0\mu s$), the flow field is dominated by a complex network of interacting shocks and expansion fans. Despite the strong wave-interface interactions and significant droplet deformation, the interfaces remain sharply resolved without numerical smearing. These results highlight the importance of accounting for multi-raindrop collective effects and demonstrate the capability of the proposed sharp-interface flux-integral method to robustly capture complex shock-raindrop interactions under realistic conditions.

4. Conclusion

A density based solution strategy and geometric VoF+PLIC method are combined to study the interaction between shock wave and multiple raindrops. A space-time flux-integral method is brought about to make the density based method consistent with VoF+PLIC.

The comparison between analytical solution, experiment results and numerical results are compared for gas-liquid shock wave problem and shock wave interaction with helium bubble. The simulations of shock wave and raindrops have provided insights into wave structures and wake flow.

Future work will focus on the detailed analysis of the bow shock wave location, the displacement of individual raindrops, and cavitation phenomena within the droplets.