2.1. Problem description

The physical configuration depicted in figure 1 illustrates a shock wave with the shock Mach number $M_s$ moving to the right and engaging with a two-dimensional (cylindrical) droplet of diameter $D_0$. Due to challenges in maintaining a perfectly spherical droplet during experiments, geometric factors could influence the shock–droplet interaction (Xiang & Wang Reference Xiang and Wang2017). However, numerical simulations can mitigate this influence. For computational efficiency, this study focuses solely on the top portion of the physical field in a two-dimensional context, delineated by the blue dashed line within a computational domain of $14D_0 \times 6D_0$ with 800 grid cells per diameter. A grid-independent test is detailed in Appendix A. The reflective boundary condition (RBC) is applied to the lower boundary to account for symmetry, while the non-reflective boundary condition (NRBC) is implemented on the upper, left and right boundaries due to the small scale of the droplet relative to the entire physical spatial scale. Detailed boundary condition definitions are available in Thompson (Reference Thompson1987, Reference Thompson1990). Initial conditions will introduce a starting error, manifesting as a sound wave opposing the direction of airflow due to the shock wave's motion, as elaborated by LeVeque (Reference LeVeque2002). Despite this minor error, the negligible intensity of the sound wave compared to the shock wave ensures that it does not impact the droplet's aerodynamic deformation or downstream effects.

Figure 1. Schematic diagram of two-dimensional shock–droplet interaction. The blue dashed line denotes the computational domain.

Both the pre-shock air and liquid droplet are initially set at atmospheric pressure. The shock wave propagates at the $M_s$ shock Mach number, inducing post-shock air characterized by increased pressure and density. Table 1 presents the pressure and density ratios ($p_{post}/p_{pre}$ and ${\rho _{post}}/{\rho _{pre}}$, respectively) between post-shock air and pre-shock air, alongside the gas–liquid wave velocity ratio ($n$) at varied shock Mach numbers for water and $n$-hexane phases. The gas–liquid wave velocity ratio is derived as the ratio of shock wave velocities in the gas and liquid phases, defined by $n={u_{g}}/{u_{l}}$, where $u_{g}$ and $u_{l}$ represent the shock wave velocities in the gas and liquid phases. Following the Boyd & Jarrahbashi (Reference Boyd and Jarrahbashi2021) methodology, cases with $n>1$ are deemed high gas–liquid wave velocity ratios whereby incident shock waves surpass internal transmitted waves in speed. Conversely, instances of $n<1$ constitute low gas–liquid wave velocity ratios, indicating shock waves trailing internal transmitted waves in velocity.

Table 1. Simulated cases and relevant shock parameters.

The Reynolds number ($Re$) and the Weber number ($We$) indicate the significance of viscosity and surface tension, respectively. The definitions in the Meng (Reference Meng2016) research of shock–droplet interaction are

(2.1)$$\begin{gather} Re=\frac{\rho_{post}\,u_{post}\,D_0 }{\mu }, \end{gather}$$

(2.2)$$\begin{gather}We=\frac{\rho_{post}\,u^{2}_{post}\,D_0 }{\sigma }, \end{gather}$$

where $\rho _{post}$, $u_{post}$, $D_0$, $\sigma$ and $\mu$ are the post-shock air density, post-shock air velocity, initial droplet diameter, surface tension coefficient and dynamic viscosity coefficient of the post-shock air, respectively. The ranges considered for $Re$ and $We$ in this study span from $2.01\times 10^{3}$ to $1.19\times 10^5$, and from $9.42 \times 10^{3}$ to $4.10 \times 10^{6}$, respectively. Notably, previous research by Meng (Reference Meng2016) and Kaiser et al. (Reference Kaiser, Winter, Adami and Adams2020) indicates that the influence of surface tension and viscous forces can be negligibly small compared to inertial forces, especially during the early stages of shock–droplet interactions.

2.2. Numerical method

This research utilizes the open-source Multi-component Flow Code (MFC), a high-order, multi-component, multi-phase and multi-scale compressible flow solver developed by Bryngelson et al. (Reference Bryngelson, Schmidmayer, Coralic, Meng, Maeda and Colonius2021) for the simulations. The governing equations are represented as

(2.3)$$\begin{gather} \frac{\partial{\alpha_1\rho_1}} {\partial{t}} +\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha_1\rho_1{\boldsymbol{u}})=0, \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial{\alpha_2\rho_2}} {\partial{t}} +\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha_2\rho_2{\boldsymbol{u}})=0, \end{gather}$$

(2.5)$$\begin{gather}\frac{\partial{\rho{\boldsymbol{u}}}} {\partial{t}} +\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho{\boldsymbol{u}}\times {\boldsymbol{u}}+p{\boldsymbol{I}})=0, \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial{E}} {\partial{t}} +\boldsymbol{\nabla}\boldsymbol{\cdot}((E+p){\boldsymbol{u}})=0, \end{gather}$$

(2.7)$$\begin{gather}\frac{\partial{\alpha _1}} {\partial{t}} +{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha _1=0, \end{gather}$$

where $\rho$ is the density, $\alpha$ is the volume fraction, $\boldsymbol {u}$ is the velocity vector, $p$ is the pressure, $\boldsymbol {I}$ is the identity matrix, and $E$ is the total energy defined as $E=\sum _{i=1}^{N_i}\alpha _i\rho _ie_i+{ \| \boldsymbol {u} \|}^2/2$, with $N_i$ the number of involved fluids, which is two in this study. The subscripted variables represent different fluids. It should be noted that there is an expansion term $K\,\boldsymbol {\nabla } \boldsymbol {\cdot } u$ in (2.3) and (2.4) for the model of Kapila et al. (Reference Kapila, Menikoff, Bdzil, Son and Stewart2000). The effect of this expansion term is provided in Appendix B, and it shows that the term has a negligible effect on the result.

The stiffened gas equation of state (SGEOS) is employed to close the five-equation model, defined as

(2.8)\begin{equation} p_i=({\gamma }_i -1){\rho}_ie_i-\gamma_i{{\rm \pi} }_{\infty,i}, \end{equation}

where $\gamma$ is the specific heat ratio, and ${\rm \pi} _{\infty }$ is the liquid stiffness (if the component is gas, then the value equals 0). The SGEOS applies to both components, including gas and liquid. For liquid, the above two parameters can be fitted by the experimental data, and the specific method can be found in Johnsen (Reference Johnsen2008). The SGEOS parameters of water can be found in Meng & Colonius (Reference Meng and Colonius2018), while the data for $n$-hexane are obtained by fitting the experimental data of Marsh (Reference Marsh1980). The obtained parameters are listed in table 2. A comprehensive selection of SGEOS parameters for liquids is presented, and the original experimental data of these parameters are assessed in Appendix C.

Table 2. Fluid properties at normal temperature and pressure and the SGEOS parameters.

Finally, the entire governing equations are closed with a series of mixture relationships:

(2.9)$$\begin{gather} 1=\alpha_1+\alpha_2, \end{gather}$$

(2.10)$$\begin{gather}\rho=\alpha_1\rho_1+\alpha_2\rho_2, \end{gather}$$

(2.11)$$\begin{gather}{\rho}e={\alpha_1\rho_1}e_1+{\alpha_2\rho_2}e_2. \end{gather}$$

The governing equations are solved using the finite-volume method and shock-capturing schemes, coupled with the Harten–Lax–van Leer–contact Riemann solver. Spatial discretization employs a fifth-order weighted essentially non-oscillatory scheme for flux reconstruction, while time integration utilizes a third-order total variation diminishing Runge–Kutta scheme. Further insights into the numerical methodologies can be found in the works of Johnsen & Colonius (Reference Johnsen and Colonius2009) and Bryngelson et al. (Reference Bryngelson, Schmidmayer, Coralic, Meng, Maeda and Colonius2021).

2.3. Validation

The two-dimensional shock–droplet interaction with $M_s=2.4$ is first simulated with the same conditions as in Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). Comparisons between the present numerical simulation results and experiments of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016) are shown in the left-hand sides of each of figures 2(a–f). The top image in each panel shows the numerical schlieren, and the bottom shows the experimental shadowgraph image of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). The instantaneous pressure contours are also shown on the right-hand sides of each of figures 2(a–f). The non-dimensional time is defined as $t^{\ast }=c_{l}t/D_0$, where $c_{l}$ is the sound speed in the liquid phase, $t$ is the physical time, and $D_0$ is the initial diameter of the droplet.

Figure 2. Comparison between the simulation results in the top part of each panel and the experimental shadowgraphs at bottom left. The pressure contours are also shown, at bottom right.

Initially, in figure 2(a), the shock wave is observed attached to the left pole of the droplet. On the right-hand side of the shock wave, a stationary droplet and pre-shock air are visible, with post-shock air located on the left-hand side. As the shock wave progresses to the right, it generates a reflected wave moving in the opposite direction outside the droplet, creating a transmitted wave moving inward in the same direction within the droplet (as depicted in figure 2b). This transmitted wave inside the droplet continues its rightward trajectory (figure 2c), with a Mach stem emerging at the intersection of the incident shock wave and the reflected wave outside the droplet. Some of the transmitted wave near the upper and lower boundaries of the droplet is reflected back due to the concave shape of the droplet, leading to a pressure decrease. Eventually, the transmitted wave reaches the right-hand boundary of the droplet and reflects as an expansion wave (figure 2d). In figure 2(e), the reflected expansion wave converges at a point, creating a small area with highly negative pressure on the droplet's centreline. The expansion wave continues its propagation until its strength diminishes, as seen in figure 2(f). Overall, the simulation results align closely with the experimental observations, which demonstrates the accuracy of the present method.

4.1. The original theoretical model

The theoretical analysis proposed in Cervenỳ (Reference Cervenỳ2001) and subsequently utilized by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) and Xu et al. (Reference Xu, Fan, Wu, Wen and Wang2023) employs the ray-tracing method to elucidate the wave transmission within the droplet. Figure 6 illustrates the schematic diagram of the initial theoretical model based on the ray-tracing method. This model operates under four fundamental assumptions (Biasiori-Poulanges & El-Rabii Reference Biasiori-Poulanges and El-Rabii2021; Xu et al. Reference Xu, Fan, Wu, Wen and Wang2023):

Assumption 1 As shock waves engage with a droplet, the droplet's corresponding location initiates a disturbance. The wavelet will spread at the speed of sound with the disturbance, and the effect of the wavelet on the wavefront can be simplified as rays according to Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) and Xu et al. (Reference Xu, Fan, Wu, Wen and Wang2023). The location of the disturbance $\beta$ and the direction of ray propagation $\theta$ should follow

(4.1)\begin{equation} \sin\theta=\sin\beta/n. \end{equation}

Figure 6. Schematic diagram of the original theoretical model (refer to Biasiori-Poulanges & El-Rabii Reference Biasiori-Poulanges and El-Rabii2021).

Assumption 2 Only rays with a generation speed exceeding the propagation speed will impact the envelope of the transmitted wave within the droplet:

(4.2)\begin{equation} \beta<\beta_{c}={\rm arcsin}(n), \end{equation}

where $\beta _{c}$ is the critical angle used to distinguish the areas that can affect the envelope. It should be noted that the critical angle differs from the limiting angle. The limiting angle is defined as the angle from the shock regular reflection to the occurrence of Mach stem according to Vijayashankar, Kutler & Anderson (Reference Vijayashankar, Kutler and Anderson1976). The limiting angle obtained in our numerical simulation results is approximately 48$^\circ$, similar to the 46$^\circ$ given by Geva, Ram & Sadot (Reference Geva, Ram and Sadot2018). The critical angle mentioned in the present work is a variable that varies with the gas–liquid wave velocity ratio.

Assumption 3 The positions where the internal rays reach the droplet boundary first and second, denoted as $Q_1$ and $Q_2$, are

(4.3)$$\begin{gather} \lambda_1=2\theta-\beta, \end{gather}$$

(4.4)$$\begin{gather}\lambda_2={\rm \pi}+4\theta-\beta, \end{gather}$$

where $\theta$ represents the initial direction of the rays upon entering the droplet, also signifying the angle at each reflection.

Assumption 4 The correlation between the trajectory of rays and time has been studied by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) and Xu et al. (Reference Xu, Fan, Wu, Wen and Wang2023):

(4.5)$$\begin{gather} t=t_{S}(\beta)+2(j-1)\,\frac{r}{u_{l}}\cos\theta+\frac{l_{Q_{j-1}M}}{u_{l}}, \end{gather}$$

(4.6)$$\begin{gather}t_{S}(\beta)=\frac{r}{u_{l}}\,(1-\cos\beta), \end{gather}$$

where $j$ represents the number of reflections within the liquid droplet, $l_{Q_{j-1}M}$ indicates the distance between the point where the liquid droplet reflects $j-1$ times, $Q_{j-1}$, and the current location, $M$. Equation (4.5) features three components: $t_{S}(\beta )$ signifies the time necessary for the incident shock wave to reach a location $\beta$ on the droplet boundary; the second component involves the time for ray propagation along the path $Q_{j-2}Q_{j-1}$; and the third component reflects the time for ray propagation from the initial contact with the droplet boundary to the current point $M$. By determining and connecting the positions $M$ of all rays at a chosen time $t$, the envelope of the wave is established. Furthermore, the disparity in expressions for $u_l$ between (4.5) and (4.6) is expounded upon in Appendix D.

4.2. The limitation of the original theoretical model

The numerical schlieren images under varying shock Mach numbers are displayed in figure 7, illustrating the progression of the Mach stem as the incident shock wave traverses from left to right. The interaction with the droplet surface triggers disturbance at specific points, distinguished by a blue line (disturbed) versus a red line (undisturbed). These disturbances extend to form the shock envelope (wavefront) indicated as the transmitted wave in figure 7. In figure 7(a), due to a low gas–liquid wave velocity ratio, the critical angle, as per (4.2) in Assumption 2, governing the transmitted wave's formation is relatively minimal. Here, a distinct separation exists between the Mach stem or the incident shock wave and the transmitted wave, indicating that disturbances beyond the critical angle do not significantly impact the wavefront's construction. Consequently, at low gas–liquid wave velocity ratios, whether it be the incident shock wave or the Mach stem dictating disturbance onset, the internal wave structure within the droplet remains largely consistent.

Figure 7. Comparison of the numerical schlieren between low and high velocity ratios with $n$-hexane as liquid phase. The blue and red lines denote the interfacial regions with disturbance by waves and without disturbance, respectively. (a) $Ms\ 2.4,\, n = 0.755$, $t^{*} = 0.848$, (b) $Ms\ 3.0,\, n = 0.943$, $t^{*} = 0.810$, (c) $Ms\ 4.8,\, n = 1.509$, $t^{*} = 0.504$, (d) $Ms\ 4.8,\, n = 1.509$, $t^{*} = 0.756$.

With an escalation in the gas–liquid wave velocity ratio, the concavity of the transmitted wave intensifies gradually, and the gap between the incident shock wave and the Mach stem, reaching the droplet, widens. The Mach stem assumes a pivotal role in triggering disturbances. In figure 7(d), the incident shock wave has moved to the right of the right pole of the droplet, underscoring the Mach stem's continued connection to the droplet, with select points yet to experience any disturbance. The original theoretical model utilizing solely the incident shock wave as the disturbance initiator is incongruent with scenarios involving high gas–liquid wave velocity ratios. Consequently, the factor governing disturbance onset ought to transition from the incident shock wave to the Mach stem.

4.3. Improvement to the original theoretical model

The subsequent discussion will elaborate on the replacement of the controlling factor for the onset of disturbance from the incident shock wave with the Mach stem. The core concept involves substituting the first term of (4.5) in Assumption 4 as delineated in § 4.1, denoted as $t_{S}(\beta )$, representing the time necessary for the incident shock wave to reach point $\beta$ on the droplet, with $t_{M}(\beta )$, denoting the time required for the Mach stem to reach the same point on the droplet. Insight into the Mach stem's behaviour can be gleaned through post-processing of numerical simulation results, enabling the fitting of the Mach stem trajectory $\beta (t_{M})$ and the derivation of the inverse function to determine $t_{M}(\beta )$.

Observing figure 8, the trajectories of the Mach stem closely align with temporal variations at various Mach numbers. In figure 8(a), the spatial positioning of the Mach stem trajectory within the Cartesian coordinate framework is presented. The spatial coordinates are non-dimensionalized by the initial droplet diameter $D_0$, so the shape is roughly like a circular shape enveloping the liquid droplet. The start of observation of the Mach stem is represented by the leftmost point, while detachment occurs as it reaches the droplet's right extremity, illustrated as the rightmost point in figure 8(b).

Figure 8. Trajectory of Mach stem with different $M_s$, in (a) Cartesian coordinate system, (b) polar coordinate system. The black line is the fitting expression for the improved model. The red line is the shock trajectory used in the original model for comparison here.

The Mach stem's trajectory is estimated by multiplying speed by time. In numerical simulations, time $t$ equates to the product of the time step $N$ and the time step size ${\rm d}t$, leading to the representation of the Mach stem's trajectory as $\mathscr {L}=uN\,{\rm d}t$. The time step size ${\rm d}t$ in this analysis is inferred to be inversely related to the shock Mach number $M_s$, denoted as ${\rm d}t \approx O({M_s}^{-1})$, with the derivation

(4.7)\begin{equation} {\rm d}t=\frac{\mathit{CFL}\times{\rm d}\kern 0.06em x}{\sqrt{{1.4p_{post}}/{\rho_{air}}}}, \end{equation}

where $CFL$, ${{\rm d}\kern 0.06em x}$, $p_{post}$ and $\rho _{air}$ are the CFL number, grid spatial resolution, pressure of post-shock air and density of pre-shock air, respectively. They are all constants except $p_{post}$, which has expression

(4.8)\begin{equation} \frac{p_{post}}{p_{pre}}=\frac{2\gamma_{air}{(M_{s}^{2}-1)}}{\gamma_{air}+1}+1, \end{equation}

where $p_{pre}$ and $\gamma _{air}$ are constant pressure of pre-shock air, and specific heat ratio of air in table 2, respectively. Now $p_{post}$ can be taken into ${\rm d}t$ to get

(4.9)\begin{equation} {\rm d}t=\frac{1.810\times10^{{-}7}}{\sqrt{2.8M_s^2-0.4}}. \end{equation}

Due to the phenomenon that the Mach stem always stays close to the incident shock wave that has velocity $M_s$, the motion speed of the Mach stem is believed to be proportional to the shock Mach number, which is $u \approx O({M_s})$. Therefore, in the present numerical simulation, the influence of the shock Mach number on the motion trajectory of the Mach stem can be eliminated due to the multiplication of the time step and the motion speed, resulting in the phenomenon that the motion trajectory of the Mach stem at different shock Mach numbers is approximately at the same spatial position at the same time step. Therefore, only a set of Mach stem is needed to fit the $\beta (t_M)$ function. The non-dimensional time $T_M$ can be expressed as

(4.10)\begin{equation} T_{M}=\frac{t_{M}}{{\rm d}t}. \end{equation}

By fitting the Mach stem trajectory in figure 8(b), we obtain the scaling law for all Mach numbers as

(4.11)\begin{equation} \beta=f(T_{M})={-}8\times10^{{-}5}(T_{M})^{2}+0.2337T_{M}+30.531. \end{equation}

After obtaining the motion trajectory of the Mach stem $\beta =f(T_{M})$, take its inverse function

(4.12)\begin{equation} T_{M}=f^{{-}1}(\beta)=1458.75-\sqrt{2.51\times10^{6}-1.25\times10^{4}\beta}, \end{equation}

and take (4.9) and (4.10) into (4.12) to get

(4.13)\begin{equation} t_{M}=\frac{2.64\times10^{{-}4}-\sqrt{8.223\times10^{{-}8}-4.095\times10^{{-}10}\beta}}{\sqrt{2.8M_s^2-0.4}}. \end{equation}

Replace the first term of (4.5) in Assumption 4 to obtain

(4.14)\begin{align} t&=\frac{2.64\times10^{{-}4}-\sqrt{8.223\times10^{{-}8}-4.095\times10^{{-}10}\beta}}{\sqrt{2.8M_s^2-0.4}}\nonumber\\ &\quad +2(k-1)\,\frac{r}{u_{l}}\cos\theta+\frac{l_{Q_{K-1}M}}{u_{l}}. \end{align}

By applying (4.14) to determine the ending rays for wavefront formation, the consideration of the Mach stem's impact is included. Appendix E addresses the implications of employing various trajectory fitting functions for the Mach stem.

Figure 8(b) contrasts the refined expression (shown as the black line) with the initial expression (depicted in red, derived from (4.6)). The discrepancy is minimal until $\beta$ approaches approximately 110$^\circ$, indicating that disregarding the Mach stem's effect has negligible consequences on theoretical analyses when the angle defined by (4.2) in Assumption 2 is below 110$^\circ$, particularly at low wave velocity ratios. Subsequently, the disparity between the two models widens notably due to curvature effects. The enhanced model aligns more closely with real-world scenarios compared to the original version, a topic elaborated on in the next subsection.

5.1. Difference between the NPP and focus points

Figure 5 illustrates that the NPP point's locations indicated by numerical simulations align with the theoretical results from Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) under low gas–liquid wave velocity ratios. However, a substantial variance emerges at higher gas–liquid wave velocity ratios. To analyse the reasons for this distinction, an improved theoretical model is employed to deduce the focus points at varying gas–liquid wave velocity ratios. This study defines the focus point obtained through theoretical analysis as the instant when the quantity of reflected waves traversing the central axis peaks within a specific range. At this moment, the location of the central axis wavefront is termed the focus point. As depicted in figure 11, the focus point is identified as the leftmost position of the wavefront (denoted by the blue line) within a range selected from one-fortieth to one-twentieth of the diameter, i.e. $R = D_0/a$, $a\in [20,40]$.

Figure 11. The definition of focus in the selected region in the theoretical model.

The temporal evolution of the number of reflected waves (the time scale is equivalent to (4.10)) is portrayed in figure 12, with the dashed line indicating the instant when the NPP emerges, as determined by numerical simulations. For shock Mach number 2.4, regardless of the range size $a$, the focusing time – marked by the maximum number of reflected waves – aligns closely with the NPP moment, suggesting congruence between the NPP and focus points in scenarios with lower wave velocity ratios. This finding supports the inference that the convergence of expansion waves can lead to pressure reduction under lower gas–liquid wave velocity ratios. Conversely, with shock wave Mach number 3.6, discrepancies between the focusing and NPP times grow; better alignment is observed with smaller ranges. At shock wave Mach numbers 5.4 and 8.4, significant deviations between focusing and NPP times occur across all selected range sizes, indicating non-coincidence of the NPP and focus points for elevated gas–liquid wave velocity ratios. Notably, at $M_{s}=8.4$, the ray promptly focuses near the right-hand wall of the droplet post-reflection, inducing a sudden shift in the focus degree ($k_{inc}/K$) from 0, as illustrated in figure 12(d). This observation demonstrates that the pressure decline at the focus point instigated by expansion waves fails to drive the pressure to its minimum, i.e. the NPP value.

Figure 12. Comparison of the NPP and focus points for three selected shock Mach numbers with $n$-hexane as the liquid. The black dashed line is the instant when the NPP happens, $k_{inc}$ is the number of reflected waves in the selected region, and $K$ is the total number of reflected waves. Here, (a) $M_{s}=2.4$, (b) $M_{s}=3.6$, (c) $M_{s}=5.4$, (d) $M_{s}=8.4$.

Figure 13 contrasts the focus point locations derived from the improved theoretical model with the NPP point location from numerical simulations, the theoretical predictions by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) and Xu et al. (Reference Xu, Fan, Wu, Wen and Wang2023), and the experimental findings of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). The horizontal axis replaces the shock Mach number $M_s$ with the gas–liquid wave velocity ratio $n$. It unveils that our improved theoretical forecast of the focus point's position aligns closely with the theoretical prediction by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) at the low $n$ region.

Figure 13. Comparison of the focus point, the NPP point with water and $n$-hexane as liquid phase, the theoretical results by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021), and the experimental results by Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016).

Our focus locations, along with data from Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) and Biasiori-Poulanges & Schmidmayer (Reference Biasiori-Poulanges and Schmidmayer2023), exhibit a consistent downward trend with increasing gas–liquid wave velocity ratios. Observably, under lower gas–liquid wave velocity ratios, the NPP point lies between our focus spot and the Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) estimate, emphasizing its dependence solely on the gas–liquid wave velocity ratio rather than the liquid phase (water or $n$-hexane) selected. As the gas–liquid wave velocity ratio escalates, disparities between focus point and NPP point's location emerge. High gas–liquid wave velocity ratios prompt numerical simulations to revert to a constant NPP point's location, aligning with the observations of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016).

In conclusion, the convergence of expansion waves contributes to pressure reduction at lower gas–liquid wave velocity ratios. However, at elevated gas–liquid wave velocity ratios, the pressure decrease arising from the expansion wave's focus point fails to reach the minimum necessary to designate this spot as an NPP point. Understanding the mechanism behind the generation of minimum pressure at NPP points entails comprehensive scrutiny of pressure-altering factors beyond solely expansion waves at the NPP and focus points.

5.2. Factors of the pressure variation

Initially, the pressure variation within the droplet during the early stages of shock–droplet interaction is categorized into four phases: the shock wave effect, relaxation effect, fluctuation effect and expansion wave effect, as illustrated in figure 14. As the transmitted wave traverses a point inside the droplet, for example, slightly to the left of the origin at $-D_{0}/8$, the pressure sharply escalates at that specific location, as observed in the studies by Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016), Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021), Xu et al. (Reference Xu, Fan, Wu, Wen and Wang2023) and Xiang & Wang (Reference Xiang and Wang2017). Subsequently, owing to uneven pressure distribution within the droplet, this point gradually aligns with the neighbouring pressure, characterizing what we term the relaxation effect in this investigation.

Figure 14. The typical variation of pressure and its classification at point $(-D_{0}/8, 0)$ inside the droplet.

Demonstrated in figure 15, the fluctuation effect arises from sustained high pressure on the droplet's upstream side, exemplified by $M_s=2.4$ and 4.8 for low and high gas–liquid wave velocity ratios, respectively. Following the interaction of the external incident shock wave with the droplet, a region of heightened pressure forms diametrically opposite to the incident direction, inducing periodic pressure surges downstream, particularly at the droplet's centre. This fluctuation phenomenon is also noted in the experimental studies of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). In the research by Schmidmayer & Biasiori-Poulanges (Reference Schmidmayer and Biasiori-Poulanges2023), the existence of the fluctuation effect is observable irrespective of whether the droplets are cylindrical or spherical.

Figure 15. Fluctuation effect in (a) low wave velocity ratio ($M_s=2.4$) and (b) high wave velocity ratio ($M_s=4.$). The numerical schlieren is shown in the upper part of the image, and the pressure contour is shown in the lower part.

Nevertheless, the fluctuation effect is relatively minor in contrast to the impact of the shock wave. Subsequently, reaching this point, the expansion wave triggers a substantial pressure decrease. What follows is the natural evolution, encompassing fluctuation pressure influences, eventually stabilizing the pressure from its negative peak value. This marks the completion of a cycle. When the expansion wave strikes the wall edge and reflects, the subsequent traversal through that point initializes the ensuing cycle.

Figure 16 provides a depiction of the pressure distribution over time at various points along the central axis at different shock Mach numbers. Each selected point experiences the four delineated processes before the second convergence of the wave system: shock effect, fluctuation effect, relaxation effect and expansion wave effect. Notably, in figure 16(a), as the point nears the droplet's downstream side, the pressure increment from the shock wave's impact diminishes. Conversely, with an escalation in the shock Mach number, as demonstrated in figure 16(d), when the point approaches the droplet's downstream region, the pressure surge induced by the shock wave intensifies considerably.

Figure 16. Temporal pressure variation at different points ($x=-D_{0}/4$, $-D_{0}/8$, 0, $D_{0}/8$, $D_{0}/4$, NPP point, focus point) in the $n$-hexane droplet at different cases, with (a) $M_s=2.4$, (b) $M_s=4.2$, (c) $M_s=4.8$, (d) $M_s=5.4$.

5.2.1. Shock effect

The shock wave effect phenomenon can be elucidated through theoretical analysis models to derive the internal wave system structure corresponding to the scenarios depicted in figures 16(a,d). In figure 17, two observation points are selected, $-D_{0}/4$ and $D_{0}/4$ from the droplet's origin, to exemplify the comparison between shock Mach numbers 2.4 and 5.4. The internal wave structure is illustrated when the ray from the droplet's left pole reaches the observation point. For a lower gas–liquid wave velocity ratio (figures 17a,b), the rays emitted from the disturbed part of the droplet, impacting the envelope surface of the internal waves (now a transmitted wave), exhibit an outward propagation angle.

Figure 17. Theoretical analysis of shock effect at the points (a,c) $x=-D_{0}/4$, (b,d) $x=D_{0}/4$, with shock Mach numbers (a,b) $M_s=2.4$, (c,d) $M_s=5.4$.

Consequently, the envelope surface shows a convex shape, termed diverging cases by Boyd & Jarrahbashi (Reference Boyd and Jarrahbashi2021). Since the selected points in figure 16 lie on the central axis, moving downstream results in the surrounding rays propagating outwards, gradually moving away from the central axis, reducing the impact from the rays. Hence in instances of lower gas–liquid wave velocity ratios, as depicted in figure 16(a), the pressure increase caused by the shock wave effect decreases as the point moves closer to the downstream wall of the droplet.

Conversely, when considering higher gas–liquid wave velocity ratios (figures 17c,d), the inward propagation angle of emitted radiation from the disturbance creates a concave shape on the envelope surface, termed converging cases by Boyd & Jarrahbashi (Reference Boyd and Jarrahbashi2021). As the observation point shifts to the right, the influence of each ray intensifies. Consequently, under a high gas–liquid wave velocity ratio illustrated in figure 16(d), points near the downstream wall experience a more pronounced pressure increase due to the shock wave effect.

5.3. Decreasing pressure caused by the reflected expansion wave

In the preceding subsection, it was observed that the focus of the reflected expansion wave does not inevitably result in the formation of NPP points in high gas–liquid wave velocity ratios, thus suggesting that it may not be the cause of cavitation. Furthermore, it is proposed that during the initial stage of shock–droplet interaction, the pressure variation within the droplet can be delineated into four phases: shock wave effect, fluctuation effect, relaxation effect and expansion wave effect.

What follows first defines the pressure decrease contribution caused by expansion waves in the pressure variation curve, $\Delta P_{exp}$, and compares it with the number of reflected waves obtained from the improved theoretical model. This comparison is employed to illustrate quantitatively that the pressure reduction attributed to reflected expansion waves constitutes just one component of the overall pressure decline. Noteworthy is the fact that focusing represents only the maximum pressure drop caused by the reflected expansion wave at that point, and it may not necessarily lead to the NPP considering the pressure variation with the four parts of the factors.

Upon scrutinizing the pressure trend depicted in figures 14 and 16, the straightforward approach to discern whether the pressure descent stems from the relaxation effect or the expansion wave is by evaluating its rate of decline. A prominent slope in the pressure trend indicates the impact of the expansion wave, as illustrated by the first-order derivative plot of pressure against time for the case $M_s=4.2$ displayed in figure 18.

Figure 18. First derivative of pressure and the definition of decreasing pressure caused by a reflected expansion wave for $M_s=4.2$.

After being subjected to the shock wave and increasing significantly, the pressure decreases at a relatively stable speed under the influence of relaxation and fluctuation. For a more refined differentiation between the effects of expansion waves and relaxation, this study introduces the definitions of the initiation time $t_1$ and cessation time $t_3$ of expansion wave influence. Here, the initiation time marks the instance when the initial trough preceding the minimum value of the first pressure derivative ($t_2$) transpires. Conversely, the cessation time denotes the moment when the first derivative of pressure first attains zero subsequent to reaching its minimum ($t_2$). The pressure drop between $t_1$ and $t_3$ attributed to the reflection of expansion waves is denoted as $\Delta P_{exp}$.

In figure 19, the comparison between the pressure drop $\Delta P_{exp}$ from reflected expansion waves at specific points and the total pressure drop $\Delta P_{total}$ from the maximum pressure post-shock wave to the minimum pressure post-expansion wave is depicted. The ratio of reflected waves $k_{inc}/K$ is illustrated as well. Here, $k_{inc}$ represents the reflected waves at the selected point, indicating the strength of focusing, while $K$ denotes the total number of waves, including those with multiple reflections or no reflection. Notably, as the point approaches the right-hand boundary of the droplet, both the number of reflected waves ($k_{inc}$) and the pressure drop from the expansion wave ($\Delta P_{exp}$) become more pronounced.

Figure 19. Reflected rays $k_{inc}/K$ with different range $a$ (dashed lines), decreasing pressure caused by expansion wave $\Delta P_{exp}$ and total decreasing pressure $\Delta P_{total}$ (solid lines) at different cases, with (a) $M_s=2.4$, (b) $M_s=4.2$, (c) $M_s=4.8$, (d) $M_s=5.4$.

Particularly in figures 19(a,b), the variations in $k_{inc}$ and $\Delta P_{exp}$ across all locations are closely aligned, suggesting that the reflected expansion wave significantly influences $\Delta P_{exp}$. In the study by Schmidmayer & Biasiori-Poulanges (Reference Schmidmayer and Biasiori-Poulanges2023), it is posited that the focus level of the reflected expansion wave becomes fixed once the shape of the transmitted wave and reflector is established, aligning with our finding that the reduction in pressure caused by the expansion wave ($\Delta P_{exp}$) and the focus level $k_{inc}/K$ exhibit parallel trends.

Nevertheless, the total pressure drop ($\Delta P_{total}$) at the NPP point surpasses that at the focal point, notably with higher shock Mach numbers. Given the increased distance between the NPP and the droplet's right-hand boundary compared to the focus point, the NPP experiences extended durations of fluctuation and relaxation effects. This underscores the importance of considering all factors influencing pressure variations, particularly the impact of transmitted shock waves and relaxation effects when determining the NPP point's location. Consequently, the negative pressure recorded at the NPP is the most extreme within the droplet, indicating a heightened likelihood of cavitation occurrence.

In this section, focus is defined using the improved theoretical model, determining the proportion of reflected waves at points along the central axis of the droplet. Comparison between the focus from theoretical analysis and the NPP point from numerical simulation reveals a close alignment under low gas–liquid wave velocity ratios. However, at higher wave velocity ratios, the focus point shifts nearer to the droplet's downstream compared to the negative pressure point.

For the exploration of the phenomenon that the NPP value can still be reached when the number proportion of reflected waves is not the maximum, four factors are used to explain the total pressure variation in the shock–droplet interaction. Finally, by comparing the pressure decrease caused by the expansion wave ($\Delta P_{exp}$) with the total pressure drop ($\Delta P_{total}$) and the number of reflected waves ($k_{inc}$), it is found that the number of reflected waves $k_{inc}$ affects the pressure drop caused by the expansion wave $\Delta P_{exp}$, and the $\Delta P_{exp}$ at the focus is indeed the largest among all points inside the droplet. However, the $\Delta P_{total}$ of the NPP point is the highest among the points inside the droplet, and the $\Delta P_{total}$ is the consequence of considering all the four factors. Hence integrating these four pressure variation elements results in a more pronounced negative pressure at the NPP point compared to the focus point.