<figure><img src="image_1.png" /></figure>

HAL
open science

# Geometry effects on the droplet shock-induced cavitation

Kevin Schmidmayer, Luc Biasiori-Poulanges

► To cite this version:

Kevin Schmidmayer, Luc Biasiori-Poulanges. Geometry effects on the droplet shock-induced cavitation. Physics of Fluids, 2023, 35 (6), pp.063315. 10.1063/5.0151404. hal-04100852

HAL Id: hal-04100852

https://hal.science/hal-04100852

Submitted on 18 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geometry effects on the droplet shock-induced cavitation

Kevin Schmidmayer<fnref n="1" />*, Luc Biasiori-Poulanges<fnref n="2" />

**Abstract**

Assessment of geometry effects affecting shock-induced cavitation within a droplet is investigated for the first time. To do this, we use a thermodynamically well-posed multiphase numerical model accounting for phase compression and expansion, which relies on a finite pressure-relaxation rate formulation and which allow for heterogeneous nucleation. These geometry effects include the shape of the transmitted wave front, which is related to the shock speed to droplet sound speed ratio, and the droplet geometry (cylindrical versus spherical). Phenomenological differences between the column and the droplet configurations are presented. In addition, the critical Mach number for cavitation appearance is determined for both cases: between $M = 1.8$ and $M = 2$ for the column, and between $M = 2$ and $M = 2.2$ for the droplet. Based on the transmitted wavefront geometry, with Mach number varying from 1.6 to 6, two cavitation regimes have been identified and the transition has been characterised: an exponentially ($M < 4.38$) and a linearly ($M > 4.38$) increasing bubble-cloud volume. On more applied aspects, we also investigate the influence of the bubble cloud on the interface disruption and compare the results against the pure liquid droplet test case. A parallel with the technique of effervescent atomization is eventually presented.

**Keywords**

Multiphase — Compressible flow — Cavitation — Shock — Droplet — Atomization

## 1. Introduction

Liquid droplets are well known to experience aerodynamic fragmentation upon interaction with a plane shock wave [1, 2]. More recently, it has been shown that in the very early stages of the shock-droplet interaction – long before the droplet starts to deform and break up – the growth of cavitation bubbles inside the droplet may occur [3, 4, 5, 6]. Such a process is likely to occur in a large spectrum of applications, as a desired or adverse effect, ranging from raindrop impact on aircraft [7], to combustion and detonation of multiphase mixtures [8], through ink-jet printing or liquid jet-based physical cleaning [9, 10], to name but a few. Because the experimental characterisation of the shock-induced cavitation within a droplet is particularly challenging [11], recent research efforts have focused on the development of models and numerical methods to account for phase change or gas phase growth within the liquid phase. Using a model incorporating phase change, Kyriazis et

al. [4] recently simulated the high-speed droplet impact on a solid substrate as experimented by Field et al. [12]. Laying the groundwork for the numerical simulation of shock-induced droplet cavitation, authors successfully showed the ability of such models to simulate the growth of bubbles. Comparing the numerical results with the experimental observations however displayed a large overestimation in the size of the bubble cloud. This is a direct consequence of the thermodynamic-equilibrium assumption which corresponds to an instantaneous equilibrium of pressures, temperatures, velocities and chemical potentials. Indeed, this approach, analogous to infinite relaxation rates for the pressures, temperatures, velocities and chemical potentials, enables the instantaneous expansion of the gas phase when subjected to a tensile wave [13, 14]. Very recently, Forehand et al. [15] modelled and studied shock-induced cavitation in cylindrical and spherical droplets and were able to well capture the initial wave dynamics. However, their investigation

<footnote n="1">¹Inria center at the university of Bordeaux, project-team CAGIRE, Université de Pau et des Pays de l'Adour, E2S UPPA, Laboratory of Mathematics and Applied Mathematics (LMAP), Pau, France</footnote>

<footnote n="2">²ETH Zurich, Department of Mechanical and Process Engineering, Institute of Fluid Dynamics, Sonneggstrasse 3, Zurich 8092, Switzerland</footnote>

<footnote>*Corresponding author: kevin.schmidmayer@inria.fr</footnote>
also did not align with the experiments conducted by Sembian et al. [3] on cylindrical droplets regarding the activity of the bubble cloud. As a result, the authors' conclusions only suggest that, like cylindrical droplets and like previous studies suggested, cavitation can occur for spherical droplets.

Within the context of heterogeneous nucleation, we recently introduced, assessed and validated a thermodynamically well-posed multiphase numerical model accounting for phase compression and expansion, which relies on a finite pressure-relaxation rate formulation [6]. Upon validation, we exploited the model to describe the phenomenology of the shock-induced cavitation at relatively low shock-wave Mach number ($1 < M < 3$) and for a cylindrical droplet, for which experimental results have already been reported [3, 12]. The bubble-cloud activity was for the first time successfully captured. Accordingly, to our knowledge, there are currently no experimental or numerical studies available in the literature that examine the impact of cavitation events on interface disruption and atomization for cylindrical or spherical droplets, regardless of whether they present low or high Mach numbers.

Consequently, we extend herein our previous work through parametric simulations aiming to evaluate for the first time the geometry effects on the droplet shock-induced cavitation process. These geometry effects include the shape of the transmitted wave front, which is related to the ratio of the shock speed to the droplet sound speed, and the droplet geometry (cylindrical versus spherical). Based on the transmitted wavefront geometry, two cavitation regimes have been identified and the transition has been characterised. On more applied aspects, we also investigate the influence of the bubble cloud<fnref n="1" /> on the interface disruption and compare the results against the pure liquid droplet test case. A parallel with the technique of effervescent atomization is eventually presented. Note that this work only considers heterogeneous cavitation, i.e., no phase change. The droplet initially contains pre-existing nuclei modelled as a liquid-gas mixture. Considering the difference in the acoustic impedance between both phases, such a modelling enables to simulate each phase response, within the mixture, to compression and expansion effects.

## 2. Problem description

### 2.1 Phenomenology

The shock-induced cavitation within a liquid droplet is initiated with the interaction of the shock wave with the droplet at time $t = 0$. Upon interaction, the shock is transmitted to the droplet, while part of the incident shock is diffracted around the droplet. It results in a compression wave propagating within the droplet in the stream direction. We refer to this wave as the transmitted wavefront or transmitted sound wave, denoted TSW. As a consequence of the large water-to-air acoustic impedance ratio, the TSW reflects at the droplet interface as a converging expansion wave. This first internal reflection of the TSW is denoted TSWr. As it propagates within the droplet, the TSWr generates low pressure regions in the internal flow field which, under some conditions, result in the cavitation and growth of bubbles forming a cloud. This bubble cloud eventually collapses and generates a spherical shock wave (CiS, for collapse-induced shock) originating from the cloud center. Upon reaching the droplet interface, the CiS similarly reflects as an expansion wave which, under some conditions, may also result in the cavitation and growth of bubbles. Note that when the expansion wave reflects at the droplet interface, it transforms into a compression wave. In this work, we focus on the effects of the TSW geometry on the shock-induced cavitation phenomenology. The parametric equations of the TSW read [5]

$$ 
\begin{aligned}
x_M &= [c_l t - n R_d (1 - \cos \alpha)] \cos(\theta - \alpha) - R_d \cos(\alpha), \\
y_M &= [c_l t - n R_d (1 - \cos \alpha)] \sin(\theta - \alpha) + R_d \sin(\alpha),
\end{aligned} 
\quad (1) 
$$

where $n$ is the water-to-air sound speed ratio $c_l/c_g$ and $R_d$ is the droplet radius. The incident and refraction angles, $\alpha$ and $\theta$, are related by the fundamental law of refraction, $\sin \theta = n \sin \alpha$. More details can be found in [5, 6].

### 2.2 Problem dimensions

The Mach number $M$ of the shock wave, the Weber number $We$, and the Reynolds number $Re$ are defined as

$$ M = \frac{U_s}{c}, \quad We = \frac{\rho U^2 d_0}{\sigma} \quad \text{and} \quad Re = \frac{\rho U d_0}{\mu}. \quad (2) $$

$U_s$ is the incident shock wave velocity, $c$ is the gas sound speed in the pre-shocked state, $\rho$ is the density of the post-shocked gas, $U$ is the post-shocked gas velocity, $\mu$ is the dynamic viscosity of the gas, $\sigma$ is the surface tension coefficient and $d_0$ is the diameter of the droplet.

<footnote n="1">¹Note that in this work, the “shock-induced bubble cloud” terminology stands for the liquid-gas mixture.</footnote>
Herein, we report high values of We, ranging from ~10⁴ to ~10⁶, and Re, from ~10⁶ to ~10⁷, for Mach numbers from 1.6 to 6. This indicates that the inertial forces dominate the flow over the surface tension and the viscous forces, respectively. These effects are therefore neglected in our modelling.

In order to facilitate the description of the shock-induced cavitation within a droplet, we use dimensionless parameters. Unless otherwise specified, non-dimensionalization of the space and time variables, L and T, is done using the initial droplet diameter $d_0 = 22$ mm and the sound speed in water $c_l$

$$ \tilde{L} = \frac{L}{d_0} \quad \text{and} \quad \tilde{T} = T \frac{c_l}{d_0}, \qquad (3) $$

where (·) denotes a non-dimensional quantity. Note that the time shown in the figures starts from the first interaction of the shock wave with the interface. We also define the dimensionless volume

$$ V^* = \frac{V_c}{V_d}, \qquad (4) $$

where $V_c = \sum_i \alpha_{g,i}V_i$ and $V_d$ are the volumes of the bubble cloud and of the droplet, respectively. $\alpha_g$ is the volume fraction of gas within the droplet and $V_i$ is the volume of the i-th cell. Cylindrical and spherical droplet volumes are identical in the beginning of the simulation to allow for better comparison. Since the initial radius is the same in both configurations, the depth of the cylindrical droplet is adapted in order to match the volume of the spherical droplet.

### 3. Numerical modelling

We use herein the modelling proposed by Biasiori-Poulanges & Schmidmayer [6] which is a slightly modified version of the modelling proposed by Schmidmayer et al. [16] to simulate the compression and expansion of each phase within the liquid-gas mixture, while ignoring phase change. The modification is only related to the form of the pressure-relaxation terms (right-hand side). The details of this modelling are provided below for the self-consistency of the paper.

#### 3.1 Governing equations

The thermodynamically well-posed, pressure- and temperature-disequilibrium, multi-component flow model conserves

mass, momentum and total energy. It reads for N phases

$$ \begin{aligned} \frac{\partial \alpha_k}{\partial t} + \mathbf{u} \cdot \nabla \alpha_k &= \delta p_k, \\ \frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}) &= 0, \\ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I}) &= \mathbf{0}, \\ \frac{\partial \alpha_k \rho_k e_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k e_k \mathbf{u}) + \alpha_k p_k \nabla \cdot \mathbf{u} &= -p_I \delta p_k, \end{aligned} \qquad (5) $$

where $\alpha_k$, $\rho_k$, $p_k$ and $e_k$ are the volume fraction, density, pressure and internal energy of each phase, respectively, and for which k indicates the phase index. The mixture density and pressure are

$$ \rho = \sum_{k=1}^{N} \alpha_k \rho_k \quad \text{and} \quad p = \sum_{k=1}^{N} \alpha_k p_k, \qquad (6) $$

while the mixture total energy is

$$ E = e + \frac{1}{2} \|\mathbf{u}\|^2, \qquad (7) $$

where e is the mixture specific internal energy

$$ e = \sum_{k=1}^{N} Y_k e_k (\rho_k, p_k). \qquad (8) $$

In (8), $e_k(\rho_k, p_k)$ is defined via an equation of state (EOS) and $Y_k$ are the mass fractions

$$ Y_k = \frac{\alpha_k \rho_k}{\rho}. \qquad (9) $$

Herein, we consider two-phase mixtures of gas (g) and liquid (l), in which the gas is modelled by the ideal-gas EOS

$$ p_g = \rho_g(\gamma_g - 1)e_g, \qquad (10) $$

and the liquid is modelled by the stiffened-gas (SG) EOS

$$ p_l = \rho_l(\gamma_l - 1)e_l - \gamma_l \pi_{\infty,l}, \qquad (11) $$

where $\gamma$ and $\pi_{\infty}$ are model parameters [17]. Herein, similarly to [6], we use $\gamma_g = 1.4$, $\gamma_l = 2.35$ and $\pi_{\infty,l} = 10^9$. The interfacial pressure is defined as

$$ p_I = \frac{\sum_k^N \left( p_k \sum_{j \neq k}^N z_j \right)}{(N-1) \sum_k^N z_k}, \qquad (12) $$
where $z_k = \rho_k c_k$ and $c_k$ are the acoustic impedance and speed of sound of the phase $k$, respectively.

For the pressure-relaxation terms between the phases, following Biasiori-Poulanges & Schmidmayer [6], $\delta p_k$ reads

$$ \delta p_k = \mu \alpha_k \sum_{j \neq k}^{N} \alpha_j (p_k - p_j). \qquad (13) $$

Note that $\mu$ is a finite constant parameter which can be selected in the $[0, \infty]$ range. However, for a given mixture and flow regime, only one value within this range accu- rately reproduces the physics. This value changes from one configuration to another and must be determined by comparison with appropriated experimental data.

Since pressures are in disequilibrium here, the to- tal energy equation of the mixture is replaced by the internal-energy equation for each phase. Nevertheless, the conservation of the mixture total energy can be writ- ten in its usual form

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot [(\rho E + p) \mathbf{u}] = 0. \qquad (14) $$

We note that (14) is redundant when the internal energy equations are also computed. However, in practice, we include it in our computations to ensure that the total energy is numerically conserved, and thus preserve a correct treatment of shock waves.

Based on the hyperbolic study, the mixture speed of sound, also called frozen speed of sound, is derived as

$$ c^2 = \sum_{k=1}^{N} Y_k c_k^2, \qquad (15) $$

which is found to be in agreement with previously re- ported expression [18].

We also recall that the model is in velocity equilib- rium, respects the second law of thermodynamics and is hyperbolic with eigenvalues either equal to $u$ or $u \pm c$, where $u$ is the velocity in the $x$-direction.

## 3.2 Numerical method

We numerically solve Eq. (5) using a splitting procedure between the left-hand-side terms associated with the flow and the right-hand-side terms associated with our relaxation procedure.

The left-hand-side terms are solved by an explicit finite-volume Godunov scheme where, to ensure the conservation of total energy, a procedure correcting the non-conservative terms of the internal-energy equations

is required and it uses the mixture total-energy rela- tion (14). The method corrects the total energy before the relaxation procedure, during the flux computation of the hyperbolic step, and therefore allows finite or infinite relaxations [16].

The relaxation terms (system of ordinary differential equations) are integrated with a first-order, explicit, Euler scheme with time-step subdivisions [16]. The number of subdivisions is adapted at each time step to verify the volume-fraction and pressure constraints. During this procedure, if the pressures are completely relaxed, i.e. a unique pressure for all phases, we terminate the Euler scheme and we perform from the initial state an infinite- relaxation procedure [18] to guarantee a unique pressure and better estimate the solution. This also assures a faster computation.

A second-order-accurate MUSCL scheme with two- step time integration is used [19], where the first step is a predictor step for the second and the usual piece-wise linear MUSCL reconstruction [20] with the monotonized central (MC) [21] slope limiter is used for the primitive variables.

In order to resolve the wide range of spatial and temporal scales of wavefronts and interfaces, an adaptive mesh refinement technique is employed [22]. The cell $i$ is refined when the following criterion is fulfilled

$$ \frac{|X_{Nb(i,j)} - X_i|}{\min(X_{Nb(i,j)} - X_i)} > \varepsilon, \qquad (16) $$

where $X$ is a given flow variable. The criterion is tested for all neighboring cells, denoted by the subscript $Nb(i, j)$, where the $j$-th cell is the corresponding neighbor of the $i$-th cell. The threshold is conservatively set to $\varepsilon = 0.02$. The above refinement criterion is tested for density, ve- locity, pressure and volume fraction and refines the cell if the criterion is fulfilled for any of these variables. In ad- dition, neighboring cells of refined cells are also refined to prevent oscillations as well as loss of precision.

This modelling is implemented in ECOGEN [19], which has been validated, verified and tested for finite- relaxation rate in various setups such as droplet shock- induced cavitation, gas bubble dynamics problems, in- cluding free-space and near-wall bubble collapses, and liquid-gas shock tubes. Using infinite-relaxation rate, it has also been validated for surface-tension problems as well as column and droplet breakup due to high-speed flow (see, e.g., [6, 16, 23] and references therein).
## 3.3 Computational setup

The calibration of the pressure-relaxation rate, $\mu$, was done by Biasiori-Poulanges & Schmidmayer [6] against the experiment of Sembian *et al.* [3], where a Mach 2.4 planar shock wave interacts with a cylindrical water droplet (column) of 22 mm in diameter. In this experiment, the growth of a bubble cloud has been imaged. The calibrated pressure-relaxation rate was found to be $\mu = 3.5$ for an initial air volume fraction in water of $\alpha_g = 10^{-6}$. This corresponds to the pre-existing nuclei in non-purified water. We recall that considering the difference in the acoustic impedance between both phases, the modelling enables to simulate each phase response, within the mixture, to compression and expansion effects, i.e. heterogeneous cavitation (without phase change). $\mu$ was again validated against a second experiment of Sembian *et al.* [3] with $M = 1.75$, for which no bubble cloud has been recorded. Note that the results presented herein are based on the assumption that this calibration on the cylindrical droplet of Sembian *et al.* also works for spherical droplets with the same initial purity (nucleus concentration). Regardless, the phenomenology would be qualitatively unaffected as it is mostly governed by the internal wave dynamics and geometry effects.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Computational setup corresponding to a planar shock wave interacting with a water droplet. The four blue segments aligned on the droplet vertical axis show the locations of the probes used to plot the pressure profiles in Fig. 3.</figcaption></figure>

The two-dimensional (2D) computational setup, corresponding to a planar shock wave interacting with a water droplet, is shown in Fig. 1, where the x-axis is the axis of symmetry on which the center of the droplet of radius $R_d = 11$ mm is located. Cylindrical axi-symmetry is used to model the spherical droplet while it is disabled to model the cylindrical droplet. Simulations are performed in a $[6R_d \times 3R_d]$ rectangular computational domain. A symmetric boundary condition is applied

to the bottom side of the computational domain, and non-reflective boundary conditions are imposed to the remaining boundaries. The droplet is initially located at the center, and is assumed to be in mechanical equilibrium with the surrounding air. The initial droplet is resolved by 100 cells per diameter. Adaptive mesh refinement (AMR) composed out of two grid levels, leading to 400 cells per diameter, and adapted to follow the flow discontinuities is used. The AMR level is selected based on the analysis of the grid sensitivity [6]. The shock wave is initialized inside the domain, and travels from left to right in air at atmospheric conditions. For the incident shock Mach number $M$, the initial flow field is determined from the Rankine-Hugoniot jump relations using a downstream density of 1.204 kg/m$^3$ and a 1 atm pressure. The water has a density of 1028 kg/m$^3$. This was calculated to agree with the sound speed in water calculated from the experimental observations of Sembian *et al.* [3] ($\approx 1512$ m/s), when using the Eq. 11.

In addition, note that a 3D simulation has been carried out on a coarser mesh to assess the quality of the 2D axi-symmetric approximation and only marginal differences were observed, arising from differences in numerical dissipation. Indeed, differences should emerge for longer observation times where surface tension and viscosity play a role. Hence, for computational reasons, 2D axi-symmetric simulations are only performed herein to represent a spherical droplet.

Further, note that the computation of the bubble cloud volume in the axi-symmetric configuration takes into account the volume of each cell, $V_i$, with a cylindrical revolution.

# 4. Results and Discussion

## 4.1 Column versus droplet

### 4.1.1 Phenomenological differences for Mach 2.4

In order to assess the phenomenological differences between the cylindrical and the spherical droplet, we take the case where the incident shock wave travels at Mach 2.4.

The first phenomenological difference concerns the wave dynamics. The geometrical difference influences the shape and amplitude of the reflected shock wave from the impingement of the incident shock wave on the droplet, and, more importantly, of the transmitted sound wave (TSW) within the droplet. Fig. 2 depicts these shape differences with the help of numerical schlieren (i.e., the exponential of the negative, normalized density
Droplet shock-induced cavitation

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of the wave dynamics before complete convergence of the rarefaction waves. Incident shock wave propagates at $M = 2.4$. Schlieren is presented for both the column and the droplet.</figcaption></figure>

gradient [24, 25]). We observe that the reflected shock wave is delayed in the spherical droplet case. This behavior is explained by the fact that the reflected shock wave travels spherically (i.e. with two curvature angles), while it is cylindrically reflected (i.e. with one curvature angle) in the case of the column. Regarding the TSW, almost no influence of the curvature angles is observed by only looking at the schlieren images. Indeed, the positions of the waves are matching between the two cases. However, pressure profiles presented in Fig. 3 show that the TSW, which is a compression wave, has a weaker amplitude in the case of the droplet. This indicates that the additional curvature angle induces a weaker energy transmission of the incident shock wave towards the interior of the droplet. We also note that the initial compression jump of the TSW is followed by a slowly decreasing pressure. Then the TSWr, which are rarefaction waves, focus and generate a sudden drop in pressure. This pressure drop is observed for the column in Fig. 3 but not for the droplet. This difference is explained by a peculiar phenomenon happening in the case of the droplet. Fig. 4 presents pressure fields in two configurations, pure water and water with initial nuclei, for both the column and the droplet. In the case of pure water, cavitation is absent. Hence, the peculiar event is easier to elucidate. During the focus of the TSWr ($\tilde{t} = 1.22$ and $1.25$), we observe higher amplitudes of negative pressures for the droplet. Although the amplitude of the TSW is weaker for the droplet, this phenomenon is expected since the additional curvature angle induces a stronger focus of the TSWr. At $\tilde{t} = 1.29$, the peculiar phenomenon appears for the droplet case. Instead of having primarily one low-pressure point (tension) as for the column, we ob-

serve two low-pressure points with short but sufficient distance from one another to generate compression in between. This compression then travels and counter-balances the rarefaction wave propagating in direction of the center of the droplet ($\tilde{t} = 1.32$). For this reason, in opposition of the column case, no sudden pressure drop is observed at the center of the droplet but rather a sudden compression ($\tilde{t} = 1.56$). Note that the compression is also observed after the rarefaction wave for the column, although with lower pressures. When nuclei are initially present, this phenomenon is weakened since a significant part of the energy is absorbed by the gas phase to generate the bubble cloud.

The second phenomenological difference concerns the bubble cloud resulting from these wave dynamics. In fact, in the case of the spherical droplet, one could think that this more pronounced curvature would induce a stronger focus of the TSWr and therefore a greater cavitation phenomenon. However, contrary to appearances, although the focus is stronger, it appears that the combination, of (i) the lower transmitted energy within the droplet with (ii) the peculiar wave-dynamics phenomenon explained above, counterbalances this focus. This results, as shown in Fig. 5, in a smaller volume of bubble cloud for the droplet compared to the column. Notably, the maximum bubble-cloud volume for the droplet is ≈ 23% of the one for the column. Further, note that the bubble clouds of both configurations start to grow at the same time and with very similar growth rates. However, the growth phase is longer for the column case.

Fig. 6 presents at longer time instants the internal structure in regard of the wave dynamics and of the bubble-cloud activity, through the numerical schlieren
<figure><img src="image_5.png" /><figcaption>Figure 3. Pressure profiles over time for four probes placed on the r axis from the droplet center, as pictured by the blue segments in Fig. 1. The results for the column and the droplet are presented for an incident shock wave propagating at $M = 2.4$.</figcaption></figure>

and the volume fraction of gas, respectively. We observe that the TSWr are more mitigated for the droplet, and the dynamics and position of the cloud is captured. For the droplet, the bubble cloud is no longer present from the image at $\tilde{t} = 1.48$. One can also note that the shock wave which could be emitted by the collapse of the bubble cloud (CiS) is not perceptible in the case of the droplet, whereas it is observed, although weak, for the column.

Consequently, these phenomenological differences, here depicted for a Mach number of 2.4, indicate that cavitation is in general less likely to occur in the spherical droplet case. The aim of the next section is therefore to confirm it and to evaluate the critical Mach number for shock-induced cavitation.

### 4.1.2 Critical Mach for cavitation appearance

In order to assess the critical Mach number presenting the first signs of cavitation appearance, it is necessary to define an indicator offering more information that a simple threshold of volume of gas phase. In consequence, we consider the relative variation of volume between two consecutive simulations:

$$ \delta V_{\max}^{*} = \frac{V_{\max,p}^{*} - V_{\max,p-1}^{*}}{V_{\max,p-1}^{*}}, \qquad (17) $$

where $V_{\max}^{*}$ is the maximum bubble-cloud volume reached during the simulation. Subscript $p$ denotes the evaluated Mach-number point and $p-1$ is therefore the preceding point (e.g., $p$ for the simulation at $M = 1.8$ and $p-1$ for the simulation at $M = 1.6$).

Fig. 7 presents $V_{\max}^{*}$ and $\delta V_{\max}^{*}$ from $M = 1.6$ to 3. Results from the simulations (markers) are completed by splines. The first observation is that the cloud volume is significantly larger for all Mach numbers in the case of the column compared to the droplet, but for below $M = 1.8$ where no cavitation is detected in both configurations. From $\delta V_{\max}^{*}$, we note for the column (droplet) that the volume is increased by 64% and 264% (48% and 179%) at $M = 1.8$ and $M = 2$ ($M = 2$ and $M = 2.2$), respectively. This indicates that cavitation is very slightly starting at $M = 1.8$ for the column ($M = 2$ for the droplet). Since it started from a very small gas volume composed of only nuclei, we do not consider it sufficient to call it a cavitation phenomenon. However, based on the sudden and significant increase in volume between $M = 1.8$ and $M = 2$ for the column ($M = 2$ and $M = 2.2$ for the droplet), we can consider cavitation to first appear within this range, which is represented with a shaded area for each configuration.

Moreover, vertical lines represent the first observations of collapse-induced shock waves (CiS). One can note that the CiS are observed for larger Mach numbers than for the first cavitation observations: $M = 2.3$ and $M = 2.5$ for the column and droplet, respectively. In fact, the smaller bubble-cloud volume combined with the milder pressure field around the cloud didn't allow for a shock wave to be emitted at the end of the collapse of the cloud. In addition, as expected, a CiS is observed for a lower Mach number for the column compared to the droplet. However, one could have expected similar cloud volume for the first CiS apparition between both configurations, but they present $V^{*} \approx 3.2 \times 10^{-5}$ and $V^{*} \approx 1.6 \times 10^{-5}$, respectively. This is explained by the fact that the cloud covers the entire thickness of the column and cylindrically collapses, whereas it presents higher concentration of bubbles/gas content and spherically collapses for the droplet.

In conclusion, we can confirm the indication of the previous section, i.e. cavitation is less likely to occur in the droplet case.

## 4.2 High Mach numbers

The influence of the Mach number on the internal pressure field is twofold: (i) the transmitted energy increases
Droplet shock-induced cavitation

<figure><img src="image_6.png" /></figure>

**Figure 4.** Comparison of the pressure fields within the column and the droplet. Results are presented for simulations with pure water and for water with initial nuclei for an incident shock wave propagating at $M = 2.4$. Note that for visualization purposes, the colorbar does not start at the minimum observed pressure, therefore, small regions of low pressure might be saturated in dark blue.

<figure><img src="image_7.png" /></figure>

**Figure 5.** Evolution of the bubble-cloud volume for the column and the droplet. The incident shock wave propagates at $M = 2.4$.

with $M$, and (ii) the geometry of the transmitted front changes from concave to convex. In the context of shock-induced cavitation, the two effects combine and affect the growth of the gas phase in the droplet. As shown in the bottom graph of Fig. 8, $V_{\max}^*$ is constant and equals to zero for $1.0 < M \lesssim 2.2$. This is due to incident shock waves not strong enough to generate cavitation upon reflection and amplification. At high Mach number, $M \gtrsim 4.5$, $V_{\max}^*$ exhibits a linear behavior, while for $2.2 \lesssim M \lesssim 4.5$ stands out as a transitional region. This latter region is a direct consequence of the change in the geometry of the transmitted wavefront (see top

graph, Fig. 8) from concave to convex with an infinite radius of curvature, $\mathcal{R} \to \infty$, appearing for $M = 4.38$. This critical Mach number corresponds to a shock speed equals to the sound speed in water, $n = c_l/u_s = 1$. Consequently, the fundamental law of refraction reduces to $\sin \theta = \sin \alpha$, where $\alpha$ and $\theta$ are the incident and refraction angles, respectively. According to the geometrical ray acoustics, the transmitted rays to the droplet are then aligned with the incident rays and propagate at the same speed which results in a plane wavefront. For $n > 1$ the acoustic rays diverge and propagate faster than the shock wave outside, which draw a spherically diverging (concave) shock. Conversely, $n < 1$ results in a spherically converging (convex) shock. The shape of the TSW as a function of the shock-wave Mach number is shown in Fig. 9. The linear trend reported at high Mach numbers are consistent with the small variation in $\mathcal{R}$ (see the dotted-dashed blue line, bottom graph of Fig. 8). The dimensionless volume $V_{\max}^*$ depends on the pressure induced by the expansion wave which corresponds to the first reflection of the transmitted wave (TSWr). The intensity of the expansion wave is related to its focusing, which depends on the shape of the TSW and the geometry of the reflector, here the droplet back face. For two TSW with relatively same $\mathcal{R}$, the amplification rate during the expansion wave focusing is nearly equal since the shape of the TSW and the geometry of the reflector are both conserved. Consequently, the only transmitted energy acts on $V_{\max}^*$, which approximately linearly increases over $M \gtrsim 6$. For large variation in $\mathcal{R}$, that is
Droplet shock-induced cavitation

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of the internal structure between the column and the droplet until the bubble-cloud collapse. The incident shock wave propagates at *M* = 2.4. The numerical schlieren is displayed and overlaid with the volume fraction of gas (white-to-blue colormap), with an opacity function to render translucent surfaces.</figcaption></figure>

for $3 \lesssim M \lesssim 5$, the amplification rate strongly varies as the geometrical configuration for the TSWr significantly changes<fnref n="2" />.

Fig. 10 shows the variation of the time-dependent gas volume within the droplet. It reports both an increase in $V^*$ over *M* as well as a phase shift in the temporal location of the maximum gas volume experienced by the droplet. Both the increase in the volume and the phase shift are related to the *M*-dependent transmitted energy, as discussed previously with Fig. 8. Dissociating the role of the TSW geometry on the phase shift from the contribution of the energy transmitted is not straightforward. However, we here assume the TSW shape to not drive the phase shift as the time required to complete the focusing is fixed by the non-deviated ray propagating along the droplet axis.

## 4.3 Bubble-cloud influence on the droplet interface

In order to study the potential influence of the bubble cloud on the droplet interface and its consequence on the atomization process, we focus on the highest Mach number: 6. Indeed, we observed the largest cloud volume for this Mach number and it is therefore expected to observe its greatest influence on the droplet dynamics, at least at the back of the droplet where the cloud is located.

To assess this influence, we compare in Fig. 11 the results given by a simulation with a droplet of water containing initial nuclei, to the results given by a simulation

involving pure water (no bubble-cloud appearance). At $\tilde{t} = 0.79$, results are almost identical, in particular in regard of the position and strength of the waves. Later, at $\tilde{t} = 1.74$, we observe the bubble cloud for the simulation containing initial nuclei. Note that the TSW reflections are completely absorbed by the cloud during its formation. Whereas waves continue to propagate for the pure-water case with potential perturbation of the interface at the front of the droplet. The cloud then collapses until $\tilde{t} \approx 3.37$ where it reaches its minimum volume and emits a shock wave (CiS). At this time, we observe a shift of the position of the back of the droplet between the two cases. In fact, the cloud during its collapse pulls the back of the droplet. During the rebound phase, this shift is reduced until almost no difference is observed and it remains as it is a *posteriori*. Indeed, only a small activity of the cloud is afterwards noticed. This leads to the assessment of the global disruption of the interface caused by the cloud. At $\tilde{t} = 5.09$, from Fig. 11 with the schlieren or from Fig. 12 where interface contours are superposed, only marginal differences of interface dynamics are observed. In other words, the bubble cloud, although significantly impacting the wave activity, did not impact the overall deformation and disruption of the interface known to be mainly governed by the interface instabilities such as Kelvin–Helmholtz or Rayleigh–Taylor instabilities [23]. This indicates that the cloud activity, for this unique impulse scenario (a unique incident shock wave), does not play a role within the atomization process. Note that similar conclusions were observed for $M = 3$, therefore they are not presented here. In addition, although using a different numerical modelling, this result is consistent with the previous work of Nykteri & Gavaises [26] mentioning that the fragmentation of a 1.9 mm droplet exposed to a $M = 2.64$ shock wave is not altered by the presence of the bubble cloud.

<footnote n="2">²Note that the present numerical model does not account for the ionization of air. At low hypersonic regimes, the air temperature increases as the molecular bonds of the air molecule increase their vibration. In the high hypersonic regimes, the high temperature of air ($> 2000^{\circ}\text{C}$) results in the atom dissociation of the oxygen molecules. During the ionization process, atoms lose electrons to form a plasma which could strongly affect the aerodynamic and hydrodynamic processes. The critical Mach number for which ionization should occur is not clear, but most probably around 8.</footnote>
Droplet shock-induced cavitation

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Evolution and variation of ($\delta$) of the maximum bubble-cloud volume, $V_{\max}^*$, for a shorter range of Mach number. The shaded areas represent Mach-number ranges where we observed the first signs of cavitation appearance. Dash-dotted vertical lines represent the first observation of a collapse-induced shock wave (CiS). The coloring of the shaded areas and of the vertical lines is consistent with the column and droplet configurations.</figcaption></figure>

When relating this result to effervescent atomization, this seems counter intuitive. Effervescent atomization is a method of atomization that involves inserting a small amount of gas into the liquid before it is atomized. This technique leads to significant improvements in performance in terms of smaller drop sizes and/or lower injection pressures [27]. As a result, one could have guessed that the presence of the bubble cloud would also have led to improved atomization. Thereby, in order to better mimic such technique, we propose an additional test case where a bubble is initially placed within the droplet. We assume that this bubble could have been voluntarily placed there or simply created during a previous pulse

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The top graph shows the variation of the transmitted wavefront curvature over $M$, where $\mathcal{R}$ is the radius of curvature of the wavefront and $\mathcal{R}^* = \mathcal{R}/R_d$. The dashed red line locates the critical Mach number $M = 4.38$ from which the wavefront is plan. For lower $M$, the transmitted wave is spherically diverging and for higher $M$ it is a spherically converging. The bottom graph shows the effect of the Mach number on $V_{\max}^*$.</figcaption></figure>

(incident shock wave). Hence, we choose the volume of the bubble to equal the maximum volume of the bubble cloud of the previous simulation with nuclei. At $M = 6$, $V_{\max}^* \approx 3.865 \times 10^{-3}$, this leads to a bubble of initial radius $R_b \approx 1.726$ mm, to compare with the 11 mm-radius droplet. We place this bubble approximately at the location of the bubble-cloud collapse point, i.e. 8.5 mm from the center of the droplet, on the right.

Fig. 13 presents a comparison for different time instants between a droplet with and without an initial bubble inside the droplet. The incident shock wave propagates once again at $M = 6$ and the water contains initial nuclei. At $\tilde{t} = 0.88$, we observe especially TSW reflection (rarefaction wave) from the impingement of TSW (compression wave) on the bubble. This TSW–bubble interaction induces the collapse of the bubble, the later occurring until $\tilde{t} \approx 3.2$. At $\tilde{t} = 1.31$, the bubble cloud is formed in both configurations. In the case with the bubble, the cloud surrounds the bubble with a small
Droplet shock-induced cavitation

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Shape of the transmitted wavefront (TSW) as a function of the shock-wave Mach number when the on-axis transmitted ray reaches the droplet center.</figcaption></figure>

layer of almost pure water separating the bubble and the cloud. Note that the outer shape of the cloud is similar in both cases. At $\tilde{t} = 3.2$, the cloud and the bubble have collapsed and CiS are observed. In addition, the formation of a thin jet is taking place at the back of the droplet. Indeed, similarly to bubble collapse near rigid or free surfaces or bubble collapse induced by a shock wave [28, 29, 30], the non-spherical collapse of the present bubble induces a jet, here directed toward the back of the droplet. This jet is more prominent at $\tilde{t} = 5.09$. Whereas, also here, the cloud activity starts to fade. Fig. 14 shows the contours at the same latest time and allows a better observation of the jet but also of the well-known shape of a bubble jetting. Note that we measure an average speed of the jet of $\approx 107$ m/s over a duration of $\tilde{t} \approx 2.06$ (30 µs), with jet speed starting approximately 40% above this value and slowly decreasing over time. In comparison to the case without the initial bubble, we note a significant difference of the droplet contour at the back of the droplet where the jet takes place, but almost none elsewhere. This perturbation of the interface from the back is expected to play a key role in the atomization process and will certainly increase its speed. However, the complete assessment for longer times is out of the scope of the paper and implies to take into account viscous and surface tension effects, while also computing 3D simulations since axisymmetry breaks for longer times [23, 31]. Furthermore,

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Growth of the gas phase within the droplet as time proceeds and as a function of the Mach number.</figcaption></figure>

this problem is a toy problem and it would make more sense to invest time and resources on real-world applications where multiple pulses are encountered and may result in enhanced atomization due to these cavitation events. In many applications, successive shocks may be transmitted to the droplet. It typically happens when a shock wave is reflected at one or multiple walls, under the droplet exposure to a pulsed source, or when the droplet-shock interaction is followed by the droplet high-speed impact. This latter example occurs during supersonic flights, when raindrops interact with shock waves (e.g., the detached bow shock) before impacting the aircraft structure and cause rain erosion damage [8].

One should note that we undertook as well a simulation where the bubble was placed in the center of the droplet, resulting in marginal differences with respect to interface disruption. This indicates that the change in water volume to gas volume ($\approx 0.3865\%$) is not the main reason interface disruption differs. Hence, the position of the cloud is important, although this position is, most of the time, not practically controlled and is directly linked to the shape and size of the droplet.

5. Conclusion

Shock-induced cavitation within a column and within a droplet have been presented with phenomenological differences arising, which constitutes the first studied geometry effect. We observed that the energy transmitted within the droplet is smaller than within the column. This, combined with a peculiar wave dynamics involving a compression wave to appear between two low-pressure points (tension), indicates that cavitation is less likely
Droplet shock-induced cavitation

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison over time between a pure-water droplet and a droplet of water containing initial nuclei. Incident shock wave propagates at *M* = 6. Schlieren is presented.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Interface-contour comparison at $\bar{t} = 5.09$ between a pure-water droplet and a droplet of water containing initial nuclei. Incident shock wave propagates at *M* = 6. Contour is taken for $\alpha_l = 0.5$.</figcaption></figure>

to occur for a droplet. The latter was then confirmed over a large Mach-number range. In addition, the critical Mach number for cavitation appearance was found to be between $M = 1.8$ and $M = 2$ for the column, and between $M = 2$ and $M = 2.2$ for the droplet.

The second geometry effect is linked to the shape of the transmitted front which changes from concave to convex when increasing the Mach number. The transition between the two appears for $M = 4.38$, corresponding to a shock speed equaling the sound speed in water. This gives rise to two regimes of cavitation: an exponentially ($M < 4.38$) and a linearly ($M > 4.38$) increasing bubble-cloud volume.

Finally, the study of the bubble-cloud influence on the droplet interface concluded that, counter intuitively, the cloud has only marginal effect on the disruption of the interface and therefore on the atomization process when the droplet is subjected to a unique pulse (unique shock wave). However, exercising with a toy model, where there initially is a bubble within the droplet, strongly

suggests that the presence of cavitation events within the droplet will enhance the atomization process when the droplet is subjected to multiple pulses. This result draws parallels with the technique of effervescent atomization. Ultimately, if time and resources are invested for real-world applications where multiple pulses are encountered, three-dimensional, viscous and surface tension effects should certainly be added to the modelling to assess interface disruption at longer times and therefore to assess atomization performance.

**Acknowledgment**

Authors acknowledge the financial support from the ETH Zurich Postdoctoral Fellowship program.

**Declaration of competing interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Author Contributions**

**Kevin Schmidmayer:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal). **Luc Biasiori-Poulanges:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal).

**Data availability**

The data that support the findings of this study are available from the corresponding author upon reasonable request.
Droplet shock-induced cavitation

<figure><img src="image_15.png" /></figure>

**Figure 13.** Comparison over time between a droplet with and without an initial bubble inside the droplet. Incident shock wave propagates at *M* = 6 and the water contains initial nuclei. Schlieren is presented.

<figure><img src="image_16.png" /></figure>

**Figure 14.** Interface-contour comparison at $\bar{t} = 5.09$ between a droplet with and without an initial bubble inside the droplet. Incident shock wave propagates at *M* = 6 and the water contains initial nuclei. Contour is taken for $\alpha_l = 0.5$.

References

[1] D.R. Guildenbecher, C. López-Rivera, and P.E. Sojka. Secondary atomization. *Experiments in Fluids*, 46(3):371, 2009.

[2] T.G. Theofanous. Aerobreakup of newtonian and viscoelastic liquids. *Annu. Rev. Fluid Mech.*, 43:661–690, 2011.

[3] S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis. Plane shock wave interaction with a cylindrical water column. *Phys. Fluids*, 28(5):056102, 2016.

[4] N. Kyriazis, P. Koukouvinis, and M. Gavaises. Modelling cavitation during drop impact on solid surfaces. *Adv. Colloid Interface Sci.*, 260:46–64, 2018.

[5] L. Biasiori-Poulanges and H. El-Rabii. Shock-induced cavitation and wavefront analysis inside a water droplet. *Phys. Fluids*, 33(9):097104, 2021.

[6] L. Biasiori-Poulanges and K. Schmidmayer. A phenomenological analysis of droplet shock-induced

cavitation using a multiphase modeling approach. *Phys. Fluids*, 35(1):013312, 2023.

[7] K. Ando. *Effects of polydispersity in bubbly flows*. PhD thesis, California Institute of Technology, 2010.

[8] J. C. Meng and T. Colonius. Numerical simulations of the early stages of high-speed droplet breakup. *Shock waves*, 25(4):399–414, 2015.

[9] H. F. Okorn-Schmidt, F. Holsteyns, A. Lippert, D. Mui, M. Kawaguchi, C. Lechner, P. E. Frommhold, T. Nowak, F. Reuter, M. B. Piqué, C. Cairós, and R. Mettin. Particle cleaning technologies to meet advanced semiconductor device process requirements. *ECS J. Solid State Sci. Technology*, 3(1):N3069, 2013.

[10] Y. Tatekura, T. Fujikawa, Y. Jinbo, T. Sanada, K. Kobayashi, and M. Watanabe. Observation of water-droplet impacts with velocities of o (10 m/s) and subsequent flow field. *ECS J. Solid State Sci. Technology*, 4(9):N117, 2015.

[11] L. Biasiori-Poulanges and H. El-Rabii. Multimodal imaging for intra-droplet gas-cavity observation during droplet fragmentation. *Optics Letters*, 45(11):3091–3094, 2020.

[12] J. E. Field, J. P. Dear, and J. E. Ogren. The effects of target compliance on liquid drop impact. *J. Applied Phys.*, 65(2):533–540, 1989.

[13] R. Saurel, F. Petitpas, and R. Abgrall. Modelling phase transition in metastable liquids: application to cavitating and flashing flows. *J. Fluid Mech.*, 607:313–350, 2008.

[14] M. Pelanti and K.-M. Shyue. A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves. *J. Comp. Phys.*, 259:331–357, 2014.

[15] R. W. Forehand, K. C. Nguyen, C. J. Anderson, R. Shannon, S. M. Grace, and M. P. Kinzel. A numerical assessment of shock-droplet interaction
tion modeling including cavitation. *Phys. Fluids*, 35(2):023315, 2023.

[16] K. Schmidmayer, J. Cazé, F. Petitpas, E. Daniel, and N. Favrie. Modelling interactions between waves and diffused interfaces. *Int. J. Numer. Meth. Fl.*, 95(2):215–241, 2023.

[17] O. Le Métayer, J. Massoni, and R. Saurel. Elaborating equations of state of a liquid and its vapor for two-phase flow models. *Int. J. Therm. Sci.*, 43:265–276, 2004.

[18] R. Saurel, F. Petitpas, and R. A. Berry. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. *J. Comp. Phys.*, 228(5):1678–1712, 2009.

[19] K. Schmidmayer, F. Petitpas, S. Le Martelot, and E. Daniel. ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows. *Comp. Phys. Com.*, 251:107093, 2020.

[20] E. F. Toro. *Riemann solvers and numerical methods for fluid dynamics*. Springer Verlag, Berlin, 1997.

[21] B. Van Leer. Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow. *J. Comp. Phys.*, 23(3):263–275, 1977.

[22] K. Schmidmayer, F. Petitpas, and E. Daniel. Adaptive Mesh Refinement algorithm based on dual trees for cells and faces for multiphase compressible flows. *J. Comp. Phys.*, 388:252–278, 2019.

[23] B. Dorschner, L. Biasiori-Poulanges, K. Schmidmayer, H. El-Rabii, and T. Colonius. On the formation and recurrent shedding of ligaments in droplet aerobreakup. *J. Fluid Mech.*, 904:A20, 2020.

[24] Y. A. Pishchalnikov, W. M. Behnke-Parks, K. Schmidmayer, K. Maeda, T. Colonius, T. W. Kenny, and D. J. Laser. High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone. *J. Acoust. Soc. Am.*, 146:516–531, 2019.

[25] J. J. Quirk and S. Karni. On the dynamics of a shock–bubble interaction. *J. Fluid Mech.*, 318:129–163, 1996.

[26] G. Nykteri and M. Gavaises. Droplet aerobreakup under the shear-induced entrainment regime using a multiscale two-fluid approach. *Phys. Rev. Fluids*, 6(8):084304, 2021.

[27] S. D. Sovani, P. E. Sojka, and A. H. Lefebvre. Efficient atomization. *Prog. Energy Combustion Sc.*, 27(4):483–521, 2001.

[28] J. R. Blake and D. C. Gibson. Cavitation bubbles near boundaries. *Annu. Rev. Fluid Mech.*, 19(1):99–123, 1987.

[29] W. Lauterborn and T. Kurz. Physics of bubble oscillations. *Rep. Prog. Phys.*, 73(10):106501, 2010.

[30] O. Supponen, D. Obreschkow, M. Tinguely, P. Kolbel, N. Dorsaz, and M. Farhat. Scaling laws for jets of single cavitation bubbles. *J. Fluid Mech.*, 802:263–293, 2016.

[31] J. C. Meng and T. Colonius. Numerical simulation of the aerobreakup of a water droplet. *J. Fluid Mech.*, 835:1108–1135, 2018.