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

HAL
open science

# A phenomenological analysis of droplet shock-induced cavitation using a multiphase modelling approach

Luc Biasiori-Poulanges, Kevin Schmidmayer

► To cite this version:

Luc Biasiori-Poulanges, Kevin Schmidmayer. A phenomenological analysis of droplet shock-induced cavitation using a multiphase modelling approach. Physics of Fluids, 2023, 35, pp.013312. 10.1063/5.0127105. hal-03894523

HAL Id: hal-03894523

https://hal.science/hal-03894523

Submitted on 12 Dec 2022

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
# A phenomenological analysis of droplet shock-induced cavitation using a multiphase modelling approach

L. Biasiori-Poulanges¹ and K. Schmidmayer²

¹) ETH Zurich, Department of Mechanical and Process Engineering, Institute of Fluid Dynamics, Sonneggstrasse 3, Zurich 8092, Switzerland.

²) INRIA Bordeaux Sud-Ouest, project-team CAGIRE, Université de Pau et des Pays de l'Adour, E2S UPPA, Laboratory of Mathematics and Applied Mathematics (LMAP), Pau, France.³

(*Electronic mail: kevin.schmidmayer@inria.fr)

(*Electronic mail: lbiasiori@ethz.ch)

(Dated: 2 December 2022)

Investigations of shock-induced cavitation within a droplet is highly challenged by the multiphase nature of the mechanisms involved. Within the context of heterogeneous nucleation, we introduce a thermodynamically well-posed multiphase numerical model accounting for phase compression and expansion, which relies on a finite pressure-relaxation rate formulation. We simulate (i) the spherical collapse of a bubble in a free field, (ii) the interaction of a cylindrical water droplet with a planar shock wave, and (iii) the high-speed impact of a gelatin droplet onto a solid surface. The determination of the finite pressure-relaxation rate is done by comparing the numerical results with the Keller-Miksis model, and the corresponding experiments of Sembian *et al*. and Field, Dear, and Ogren, respectively. For the latter two, the pressure-relaxation rate is found to be $\mu = 3.5$ and $\mu = 0.5$, respectively. Upon validation of the determined pressure-relaxation rate, we run parametric simulations to elucidate the critical Mach number from which cavitation is likely to occur. Complementing simulations with a geometrical acoustic model, we provide a phenomenological description of the shock-induced cavitation within a droplet, as well as a discussion on the bubble-cloud growth effect on the droplet flow field. The usual prediction of the bubble cloud center, given in the literature, is eventually modified to account for the expansion wave magnitude.

## I. INTRODUCTION

The interaction of a liquid droplet with a shock wave, or the impact of the droplet on a solid substrate, results in the transmission of a confined shock to the droplet. Because of the large acoustic impedance ratio between the surrounding gas and the liquid, and so a poor gas-to-liquid energy transfer³, the droplet interface acts as a nearly perfect mirror which traps the transmitted wave energy within the droplet. The confined shock therefore experiences near total reflections, as well as focusing which results in the amplification of the shock local interaction with the liquid⁴. The reflected wave of the transmitted shock being an expansion wave, some regions of the droplet are thus exposed to a tensile force which, under some conditions, may generate hydrodynamic cavitation (i.e., isothermal inertia-driven phase change). Liquid rupture arises when subject to tension exceeding a threshold value which depends on the nature of the liquid and its purity. Pure liquids cavitate when the random thermal motions of molecules cause microscopic voids⁵. This process is usually referred to as homogeneous cavitation. When liquids are not pure, i.e. contain pre-existing nuclei/impurities, cavitation results from the expansion of submicroscopic gas pockets trapped on particles present in the liquid. The process of bubble formation by this mechanism is referred to in the literature as heterogeneous nucleation.

Shock-induced cavitation within a droplet, upon impact or interaction with a shock wave, occurs in wide range of appli-

cations, as a desired or adverse effect, ranging from raindrop impact on aircraft⁶, to combustion and detonation of multiphase mixtures⁷, through ink-jet printing or liquid jet-based physical cleaning⁸,⁹, to name but a few. The comprehension of the bubble dynamics within the droplet is thus of major importance to evaluate the erosion efficiency of the bubble-compounded droplet, related to the collapse and jetting processes of the cavitation bubbles. Given that the presence of cavities inside the droplet alters its interfacial dynamics¹⁰,¹¹, under some conditions, changes in the fragmentation process are to be expected.

The experimental characterization of shock-induced cavitation within a droplet is particularly challenging¹¹. By reducing the droplet to a water column, Sembian *et al*. and Field, Dear, and Ogren however successfully imaged the growth of a bubble cloud during the interaction of a cylindrical droplet with a planar shock wave, and the high-speed impact of a cylindrical droplet with a solid substrate, respectively. To overcome the experimental limitations and address the shock-induced cavitation within a droplet under near-reality conditions, previous attempts to explicit the conditions for the bubbles to grow mostly relied on numerical simulations. As a first approach, past numerical studies used numerical models not accounting for phase change or phase expansion¹,⁴,¹²⁻¹⁴. The occurrence and intensity of the cavitation were evaluated by probing the pressure field and comparing the low-pressure region magnitude to the cavitation threshold given by the classical nucleation theory for homogeneous cavitation (134 MPa at 300 K)¹⁵, or the Blake threshold pressure¹⁶. Recently, Kyriazis, Koukouvinis, and Gavaises simulated the experiment of Field, Dear, and Ogren, i.e. high-speed droplet impact,

<footnote>a) https://kevinschmidmayer.github.io/</footnote>
using a thermodynamically well-posed model incorporating phase change. They successfully demonstrated that such models are well adapted to simulate the growth of bubbles and to examine its effect on the droplet dynamics. However, a direct comparison of the numerical results with the experimental observations revealed that the numerical model significantly overestimates the size of the bubble cloud. This is because of the thermodynamic equilibrium assumption, corresponding to an instantaneous equilibrium of pressures, temperatures and velocities, in other words an analogous to infinite relaxation rates for the pressures, temperatures and velocities, which enable the instantaneous expansion of the gas phase when subjected to an tensile wave.

In this work, we introduce a multiphase numerical model, in velocity equilibrium, with a finite pressure-relaxation rate, $\mu$, to address the over-expansion of the gas phase as previously reported for an infinite $\mu$. The finite pressure-relaxation rate is defined on the $[0, \infty]$ range. Shock-induced cavitation primary resulting from heterogeneous cavitation, we do not account for phase change. The droplet containing pre-existing nuclei is 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. We first simulate the spherical collapse of an air bubble in a free field, over the $\mu$ range, and compare the results with the solution of the Keller-Miksis equation to eliminate $\mu$ values that do not agree with the theoretical bubble behavior. Secondly, we simulate the experiment of Sembian et al. with a Mach 2.4 shock wave for which, cavitation bubbles have been imaged. After investigating the influence of $\mu$ on the shock-induced cavitation, we calibrate the finite pressure-relaxation rate against the experimental image, and eventually validate the calibrated value by computing the experiment with a Mach 1.75 shock wave. A phenomenological analysis of the shock-induced cavitation is eventually proposed by interpreting the droplet internal wave pattern, computed from the numerics, and complemented with the geometrical acoustic model of Biasiori-Poulanges and El-Rabii. We finally examine the sensitivity of $\mu$ on the material properties of the liquid mixture constituting the droplet, by simulating the high-speed droplet impact experiments of Field, Dear, and Ogren.

## II. PROBLEM DESCRIPTION

The interaction of a confined fluid volume with a shock wave is known to generate a complex time dependent wave pattern. Accounting for the compression and expansion effects in a two-phase liquid-gas droplet, the canonical wave structure is modified. Based on the droplet internal wave structure, this section first gives a phenomenological description of the shock-induced cavitation within a liquid droplet. It also introduces recent works on the analytical description of the wave pattern, which has been interpreted using the classical ray-tracing approach to geometrical acoustics. Note that, in this section, the description is based on the interaction of a shock wave with a droplet, but the phenomenology is also

valid for the high-speed droplet impact.

### A. Phenomenology

The phenomenology of the shock-induced cavitation within a liquid droplet is sketched in Fig. 1, where the wave pattern inside the droplet is drawn as time proceeds. The time origin, $t=0$, corresponds to the instant at which the shock wave interacts with the droplet. This interaction results in the transmission of a shock to the droplet [Fig. 1(a)], while part of the incident shock is diffracted around the droplet. The transmitted shock is a compression wave that spherically propagates in the stream direction. When the transmitted shock meets the droplet boundary, and as a consequence of the large water-to-air acoustic impedance ratio, the transmitted shock reflects at the interface as an expansion wave [Fig. 1(a)], thereby forming low pressure regions in the internal flow field. At the early stage, the acoustic ray theory has shown this reflection to be a two-segment wavefront. On reaching the downstream droplet surface, the transmitted shock is completely reflected back [Fig. 1(b)]. The only expansion wave remains and propagates upstream by converging and amplifying due to the droplet curvature. The low-pressure region generated by the expansion wave thus locally exposes the liquid to a pulling force which, under some conditions, results in the cavitation and growth of bubbles [Fig. 1(a-d)]. Once the convergence of the expansion wave is completed, it diverges by shaping an horseshoe [Fig. 1(d)]. Before this transition, portion of the expansion wave crosses at the droplet axis where the bubbles eventually meet, and form a single bubble cloud. The bubble cloud collapses over time [Fig. 1(d-f)], while the expansion wave continues to propagate upstream. When the collapse is completed, a spherical shock wave (CiS) originating from the cloud center is emitted [Fig. 1(f)]. Upon reaching the droplet interface, the CiS similarly reflects as an expansion wave [Fig. 1(g)] which, under some conditions, may also result in the cavitation and growth of bubbles. The wave pattern drawn in Fig. 1(f) and Fig. 1(g) corresponds to some of the successive reflections of the transmitted shock, i.e., the reflections of the TS reflection. Note that when the expansion wave reflects at the droplet interface, it transforms into a compression wave.

### B. Analytical wave description

The wave pattern inside a spherical water droplet impacted with a planar shock wave has recently been extensively analyzed. Within the context of ray theory, Biasiori-Poulanges and El-Rabii described the time-dependent shape of the internal wavefront whose dominant feature has been shown to be the existence of cusp singularities, and examined in details the focusing of the singly-reflected wavefront. Authors also derived the parametric equations for the transmitted wavefront
<figure><img src="image_3.png" /><figcaption>FIG. 1: (a-g) Sketch of the internal wave pattern and phenomenology of the shock-induced cavitation within a liquid droplet. Time is indicated with a red-to-blue colormap, with t = 0 the instant at which the incident shock reaches the droplet. Not all internal reflections are drawn for the sake of clarity and educational purposes.</figcaption></figure>

and its multiple internal reflections, which read

$$x_M = [c_l t - n R_d (1 - \cos \alpha) - 2(k-1)R_d \cos \theta] \cos(\gamma_k - \theta) - R_d \cos(\gamma_k - 2\theta), \quad (1a)$$

$$y_M = [c_l t - n R_d (1 - \cos \alpha) - 2(k-1)R_d \cos \theta] \sin(\gamma_k - \theta) - R_d \sin(\gamma_k - 2\theta), \quad (1b)$$

where $k = 1$ corresponds to the transmitted shock and $k = 2, 3, ...$ to the successive internal reflections, $n$ is the water-to-air sound speed ratio $c_l/c_g$ and $\gamma_k = 2k\theta - \alpha - (k-1)\pi$. $R_d$ is the droplet radius. We denote $\mathcal{F}_k$ the wavefront associated to the $k$-th reflection. The incident and refraction angles, $\alpha$ and $\theta$, are related by the fundamental law of refraction, $\sin \theta = n \sin \alpha$.

As visible in Fig. 1, internal reflections are two-segment fronts which exhibit a singular point. These points trace out surfaces that are the caustics of the associated $k$-th reflection, which also have a singular point (Fig. 2). From the singularities of the energy flux density of the refracted wave, Biasiori-Poulanges and El-Rabii derived the parametric equations of the caustic surface associated with the $k$-th reflected wave-front, which are given by

$$x_{\text{caustic}} = R_d f(\alpha) \cos \gamma_k + R_d [f(\alpha) - 1] \cos(\gamma_k - 2\theta), \quad (2a)$$

$$y_{\text{caustic}} = R_d f(\alpha) \sin \gamma_k + R_d [f(\alpha) - 1] \sin(\gamma_k - 2\theta), \quad (2b)$$

where

$$f(\alpha) = \frac{1}{2} \frac{2n^2(k-1)\sin 2\alpha - \sin 2\theta}{n^2(2k-1)\sin 2\alpha - \sin 2\theta}. \quad (3)$$

<figure><img src="image_4.png" /><figcaption>FIG. 2: Caustic traced out by the singular point ($\mathcal{S}$) of the singly-reflected wavefront. The caustic exhibits a cusp ($\mathcal{C}$). The red-to-blue colored segments are the singly-reflected wavefronts as time proceeds.</figcaption></figure>

The cuspidal point of the caustic is located on the droplet axis $y_{\mathcal{C}} = 0$, while the $x$-coordinate has been found to be

$$x_{\mathcal{C}} = \frac{(-1)^k n}{(2k-1)n-1} R_d. \quad (4)$$
<table><caption>TABLE I: Weber and Reynolds numbers associated to the experiments of Sembian *et al.* and Field, Dear, and Ogren.</caption><thead><tr><th>Configuration</th><th>We</th><th>Re</th></tr></thead><tbody><tr><td>Shock-droplet<sup>1</sup></td><td>~ 10<sup>3</sup></td><td>~ 10<sup>6</sup></td></tr><tr><td>High-speed impact<sup>2</sup></td><td>~ 10<sup>6</sup></td><td>~ 10<sup>6</sup></td></tr></tbody></table>

### C. 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 (5)$$

In the configuration of the shock-droplet interaction (resp. high-speed droplet impact), $U_s$ is the incident shock wave velocity (resp. impact velocity), $c$ is the gas sound speed in the pre-shocked state (resp. sound speed in the liquid), $\rho$ is the density of the post-shocked gas (resp. droplet density), $U$ is the post-shocked gas velocity (resp. impact velocity), $\mu$ is the dynamic viscosity of the gas (resp. dynamic viscosity of the liquid), $\sigma$ is the surface tension coefficient and $d_0$ is the diameter of the cylindrical droplet.

Table I reports high values of $We$ and $Re$ indicating that, in both experiments, the inertial forces dominate the flow over the surface tension and the viscous forces, respectively.

In addition to the shock-droplet interaction and high-speed droplet impact configurations, we herein also simulate the spherical collapse of an air bubble in a free field. Viscous and capillary effects are trivially shown to be also negligible by computing the Rayleigh-Plesset equation. In this work, in- viscid flows are therefore modelled and capillary effects are not accounted for.

The phenomenology of the shock-induced cavitation within a droplet is described using 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$ 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}, \quad (6)$$

where $(\tilde{\gamma})$ denotes a non-dimensional quantity.

### III. NUMERICAL MODELLING

We use herein a slightly modified version of the modelling proposed by Schmidmayer *et al.*¹⁸ 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) and is detailed in Section III B.

### A. 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{align}
\frac{\partial \alpha_k}{\partial t} + \mathbf{u} \cdot \nabla \alpha_k &= \delta p_k, \nonumber \\
\frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}) &= 0, \tag{7} \\
\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I}) &= \mathbf{0}, \nonumber \\
\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, \nonumber
\end{align}$$

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, \quad (8)$$

while the mixture total energy is

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

where $e$ is the mixture specific internal energy

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

In (10), $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}. \quad (11)$$

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

$$p_g = \rho_g(\gamma_g - 1)(e_g - e_{g,\text{ref}}), \quad (12)$$

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

$$p_l = \rho_l(\gamma_l - 1)(e_l - e_{l,\text{ref}}) - \gamma_l \pi_\infty, \quad (13)$$

where $\gamma$, $e_{ref}$ and $\pi_\infty$ are model parameters¹⁹. The interfacial pressure is defined as

$$p_I = \frac{\sum_k^N (p_k \sum_{j \neq k}^N z_j)}{(N-1) \sum_k^N z_k}, \quad (14)$$

where $z_k = \rho_k c_k$ and $c_k$ are the acoustic impedance and speed of sound of the phase $k$, respectively.

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

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

We note that (15) 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, \quad (16) $$

which is found to be in agreement with previously reported expression$^{20}$.

We also recall that the model is in velocity equilibrium, 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.

## B. Expression of $\delta p_k$

For the pressure-relaxation terms between the phases, $\delta p_k$ reads under its general form

$$ \delta p_k = \sum_{j \neq k}^{N} \mu_{k,j} (p_k - p_j), \quad (17) $$

where the relaxation coefficients $\mu_{k,j}$, related to the $k-j$ interactions ($j$ are components different from $k$) and appearing in the original form of the complete disequilibrium model$^{21}$, can be expressed under different forms.

In most if not all the literature, e.g. Schmidmayer et al. $^{18}$, Saurel, Petitpas, and Berry$^{20}$, Baer and Nunziato$^{22}$, Saurel and Abgrall$^{23}$, the relaxation coefficients are taken as unique and constant for all interactions, i.e. $\mu_{k,j} = \mu$.

Herein, we propose to use a different approach and to express them as $\mu_{k,j} = \alpha_k \alpha_j \mu$, where $\mu$ is a constant parameter. This leads to

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

First, this expression is consistent for $N$ phases, meaning $\sum_k \delta p_k = 0$, and it does not alter the model properties (first and second law of thermodynamics and hyperbolicity). Second, the combination of the volume fractions allows specific behaviours:

* a dilute phase takes time to reach equilibrium with the carrier fluid,

* two phases with approximately the same volume fraction (close to 0.5), e.g. interfaces, or bubbles or droplets of approximately the size of the computational cells, quickly reach equilibrium, and

* two dilute phases within a carrier fluid hardly communicate and therefore take a significantly long time to reach equilibrium through their own interactions.

Hence the local relaxation rate adapts over time to the volume-fraction configuration within the computational cell. Note that $\mu$ is a finite parameter which can be selected in the $[0, \infty]$ range. However, for a given mixture and flow regime, only one value within this range accurately reproduces the physics. This value change from one configuration to another and must be determined by comparison with appropriated experimental data.

## C. Numerical method

We numerically solve Eq. (7) 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 relation (15). 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$^{18}$.

The relaxation terms (system of ordinary differential equations) are integrated with a first-order, explicit, Euler scheme with time-step subdivisions$^{18}$. 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$^{20}$ to guarantee a unique pressure and better estimate the solution. This also assures a faster computation.

As a side note, after applying an infinite pressure relaxation ($\mu = \infty$), the model converges to the mechanical-equilibrium model of Kapila et al. $^{24}$ and the effective mixture speed of sound matches Wood's

$$ \frac{1}{\rho c_w^2} = \sum_{k=1}^{N} \frac{\alpha_k}{\rho_k c_k^2}. \quad (19) $$

A second-order-accurate MUSCL scheme with two-step time integration is used$^{25}$, where the first step is a predictor step for the second and the usual piece-wise linear MUSCL reconstruction$^{26}$ with the monotonized central (MC)$^{27}$ 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$^{28}$. 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)} > \epsilon, \quad (20) $$
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 $\epsilon = 0.02$. The above refinement criterion is tested for density, velocity, pressure and volume fraction and refines the cell if the criterion is fulfilled for any of these variables. In addition, neighboring cells of refined cells are also refined to prevent oscillations as well as loss of precision.

This modelling is implemented in ECOGEN²⁵, which has been validated, verified and tested for finite-relaxation rate in various setups such as gas bubble dynamics problems, including 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.,¹⁸,²⁹⁻³³).

## D. Computational setup

The determination of the pressure-relaxation rate $\mu$ is a two-step approach. The first step consists in reducing the $\mu$-range by simulating the spherical collapse of a bubble in a free field and determined $\mu$ values that agree with the predicted bubble dynamics given by the Keller-Miksis equation³⁴. The second step consists in calibrating $\mu$ against the experiment of Sembian et al., where a Mach 2.4 planar shock wave interacts with a cylindrical water droplet. In this experiment, the growth of a bubble cloud has been imaged. The calibrated $\mu$ is eventually validated against a second experiment of Sembian et al. with $M = 1.75$, for which no bubble cloud has been recorded. This determination procedure is done for a given fluid, that is water droplet. To evaluate the sensitivity of $\mu$ on the material properties of the mixture constituting the droplet, the experiment of Field, Dear, and Ogren, also showing shock-induced cavitation, is simulated. It consists in a spherical gelatin droplet which is impacted by a metallic slider at 110 m/s.

### 1. Spherical bubble collapse in a free field

This test case aims to present the behavior of the relaxation rate during the spherical collapse of a bubble in a free field [Fig. 3(a)]. To reduce the computational cost, a one-dimensional (1D) domain of 3 mm long is used with spherical axi-symmetry to mimic a three-dimensional bubble³¹. The domain consists in high-pressure water at $p_\infty = 50$ atm with density $\rho_l = 1000$ kg/m³. An air bubble is located at the origin of the computational domain. The initial bubble radius is $R_0 = 0.1$ mm. The bubble pressure and density are $p_b = 3550$ Pa (i.e., vapor pressure) and $\rho_b = 0.027$ kg/m³, respectively. Initial velocities are nulls. One could note that this configuration enforces an initial interface disequilibrium. A non-reflecting boundary condition is used at the far field limit while a symmetry boundary condition is used at the origin of the domain. The mesh contains 150 cells from 0 to 0.3 mm, which corresponds to 100 cells per bubble diameter, and then

the grid is stretched non-uniformly to accommodate the large computational domain.

<figure><img src="image_5.png" /><figcaption>FIG. 3: (a) Volumetric representation of the spherical collapse of a bubble in a free field. (b,c) Computational setups corresponding to the experiments of Sembian et al. and Field, Dear, and Ogren, respectively.</figcaption></figure>

### 2. Cylindrical droplet interaction with a planar shock wave

The two-dimensional (2D) computational setup, corresponding to Sembian et al. experiments, is shown in Fig. 3(b), where the x-axis is the axis of symmetry on which the center of the cylindrical droplet of radius $R_d$ is located. Simulations are performed in a $[12R_d \times 6R_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 three grid levels and adapted to follow the flow discontinuities is used. The AMR level is selected based on the analysis of the grid sensitivity (see section III E). 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³ and a
1 atm pressure. The water has a density of 1028 kg/m³ and is modelled using the SG EOS¹⁹,ᵃ with $\gamma = 2.35$ and $\pi_{\infty} = 10^9$. The initial air volume fraction in water is $10^{-6}$. This corresponds to the preexisting 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).

### 3. High-speed cylindrical droplet impact

The computational setup corresponding to Field, Dear, and Ogren experiments is shown in Fig. 3(c). The spherical droplet is modelled using a 2D formulation, where the y-axis is the axis of symmetry on which the center of the droplet of radius $R_d$ is located. Simulations are performed in a $[6R_d \times 6R_d]$ square computational domain. A wall boundary condition and a symmetry boundary condition are applied to the bottom and left sides of the computational domain, respectively. Non-reflective boundary conditions are imposed to the two remaining boundaries. The droplet moves downward with a velocity $u_i = 110$ m/s. The initial droplet is resolved by 100 cells per diameter and the three-grid level AMR is used. In the experiments, the droplet is made of a 12 wt% gelatin. It is modelled using the SG EOS¹³,³⁵ with $\gamma = 6.72$ and $\pi_{\infty} = 3.70 \times 10^8$. Similarly to the previous computational setup, the initial air volume fraction within the droplet is $10^{-6}$. As a simplification of the experiments, the present setup does not account for the material properties of the metallic slider, used in the experiments, and simulates a droplet impacting a wall at velocity $110$ m/s¹⁷.

### E. Grid convergence for finite pressure relaxation

To consider the spatial convergence of the numerical method, a grid resolution study is performed by simulating the interaction of the cylindrical water droplet with a Mach 2.4 plane shock wave at four different resolutions (Fig. 4). Keeping constant the initial mesh size to 100 cells per diameter, four AMR grid levels are used which eventually result in 100, 200, 400 and 800 cells per diameter. For this study, the pressure-relaxation rate is chosen equal to $\mu = 10$. As time proceeds, the grid sensitivity is examined by comparing the growth and collapse of the shock-induced bubble cloud, within the droplet. Note that in this work, the “shock-induced bubble cloud” is a liquid–gas mixture. Figure 4 plots the equivalent radius of the gas phase within the droplet, denoted $R$ and defined as

$$R = \left( \frac{1}{\pi} \int_{V_d} \alpha_g dV \right)^{1/2}, \qquad (21)$$

where $V_d$ is the droplet volume.

<figure><img src="image_6.png" /><figcaption>FIG. 4: Growth and collapse of the gas phase within the droplet at four grid resolutions. The pressure-relaxation rate is $\mu = 10$. The droplet is initially resolved by 100 cells per diameter.</figcaption></figure>

One can observe convergence of the bubble-cloud radius. We consider the solution with three levels of refinement to be sufficiently converged for the purpose of the paper.

The rate of convergence is also presented in Figure 5 in terms of the discrete $L_2$-error $\epsilon$ as

$$\epsilon = \frac{1}{N_t} \sum_{i=0}^{N_t} \frac{\|R(t_i) - R_{\text{ref}}(t_i)\|}{R_{\text{ref}}(t_i)}, \qquad (22)$$

where $N_t$ is the number of time steps in the temporal window $\tilde{t} \in [0, 3.25]$, $R(t_i)$ is the bubble radius at time $t_i$ of our simulations and $R_{\text{ref}}$ is the reference solution, here chosen as the solution of the simulation AMR level 3. We see that the method converges at 1.7 order, matching the expected rate for AMR simulations of flows mainly governed by tension waves and exhibiting shocks and interfaces.

## IV. RESULTS AND DISCUSSION

In this section, we discuss the results from the three-step procedure for the determination of the pressure-relaxation rate.

### A. Spherical bubble collapse in a free field

The finite pressure-relaxation-based method allows for an infinite range of pressure-relaxation rate. As a first approach, it is instructive to evaluate the sensitivity of $\mu$ by simulating the spherical collapse of a bubble in free field [Fig. 3(a)] for various orders of magnitude of $\mu$, and compare the bubble response to the analytical solution of the Keller–Miksis equation³⁴; the compressible form of the Rayleigh–Plesset equa-

<footnote>ᵃ The water density is calculated to agree with the sound speed in water calculated from the experimental observations of Sembian *et al.* (≈ 1512 m/s), when using the Eq. 13.</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 5: Convergence order of the numerical method. The discrete $L_2$-error, $\epsilon$, representative of the convergence rate is given by Eq. 22.</figcaption></figure>

tion. Assuming a spherical collapse, the Keller–Miksis equation is based on an asymptotic expansion in Mach number of the Bernoulli equation. The use of the Keller–Miksis equation in the present work is predicated on larger measured relative errors than errors related to the asymptotic expansion and the inherent presumption of sphericity. This assumption is born out by the results displayed in Fig. 6, which discloses the comparison between numerical simulations and the Keller-Miksis solution by plotting the radial bubble-wall evolution. In agreement with the 1D formulation, the effective bubble radius, $R_b$, is defined as

$$R_b = \sum_{i=1}^{N} \alpha_{g,i} V_{c,i}, \qquad (23)$$

where $N$ is the number of grid cells, and $\alpha_{g,i}$ and $V_{c,i}$ are the gas volume fraction and the volume of the $i$-th cell, respectively. The dimensionless time is given by the ratio of the dimensional time $t$ with the Rayleigh collapse time,

$$t_r = 0.915 R_0 \sqrt{\frac{\rho_l}{p_\infty}}, \qquad (24)$$

which is the nominal total collapse time, i.e. the time required for the bubble to complete its collapse³⁶. Note that the solutions are only displayed until $t = 1.05t_c$, right after the minimum bubble radius is reached, since the subsequent rebounds for large pressure ratios for the Keller–Miksis equation are well-known to be physically inaccurate³⁷. When comparing the Keller-Miksis equation with our simulations, our results are expected to converge towards the analytical solution up to the first rebound because of the reducing diffusion. Beyond the collapse, the solution of the Keller-Miksis equation is here not accurate as the Mach of the interface is high at the time of the collapse, and so the compressibility effects are important. These conditions appear to be out of the validity domain on which the Keller-Miksis equation has been derived (i.e. low Mach number). The convergence should thus not occur.

<table><caption>TABLE II: Relative error, $\epsilon_r$, between the numerics and the Keller–Miksis solution on the minimum bubble radius.</caption><thead><tr><th>μ</th><th>0.01</th><th>0.1</th><th>1</th><th>10</th><th>∞</th></tr></thead><tbody><tr><th>ε<sub>r</sub> (%)</th><td>735</td><td>275</td><td>10</td><td>19</td><td>18</td></tr></tbody></table>

Figure 6 and Table II show that pressure-relaxation rates $1 \lesssim \mu \lesssim \infty$ are in satisfying agreement with theory, while underestimating the minimum radius up to 19%. The best agreement with theory is given for $\mu = 1$, where the relative error $\epsilon_r$ on the minimum bubble radius, with respect to the Keller-Miksis solution, is 10%. Note that $\min(R_b/R_0)|_{\mu=10} \approx \min(R_b/R_0)|_{\mu=\infty}$. One should also note that decreasing errors is expected for refined grids and for three-dimensional simulations, for which the non-conservative terms of the axisymmetry are absent. Major discrepancies are however reported for $\mu \lesssim 0.1$, which significantly overestimate the minimum bubble radius.

From this first sensitivity analysis on $\mu$, the range of pressure-relaxation rate can be restricted to $[0.1, \infty]$. Note that this range can be extended when refining the mesh size, so that $\mu \in [a, \infty]$ with $a \to 0$ for an infinitely small mesh size. The $\mu$-range $[0.1, \infty]$ is here valid for a reasonable resolution of 100 cells/diameter. In the following, we calibrate $\mu$ upon experimental observations which is therefore expected to be in the $\mu$-domain here determined.

## B. Cylindrical droplet interaction with a planar shock wave

In this section, we simulate the experiments of Sembian *et al.* to first investigate the influence of $\mu$ on the droplet internal flow field, and then to calibrate and validate the corresponding pressure-relaxation rate. Upon validation, parametric simulations are eventually run to evaluate the critical Mach number from which, bubbles are likely to grow.

### 1. Influence of the pressure-relaxation rate

The effect of the pressure-relaxation rate is investigated by comparing numerical results at four different rates, 1, 10, 100 and $\infty$, and by using the computational setup dedicated to the shock-droplet interaction [Fig. 3(b)]. Figure 7 shows that larger values of $R$ are reached for higher values of $\mu$. The maximum radius asymptotically increases with $\mu$, so that $\max(R/R_0)|_{\mu=100} \approx 2\max(R/R_0)|_{\mu=10}$ and $R(\mu = 100) \approx R(\mu = \infty)$. The simulation thus converges to an infinite pressure-relaxation rate modelling as $\mu$ increases. Conversely, smaller $R$ are reported for lower values of $\mu$. Numerical results computed for $\mu = 1$ seem to approximate simulations that do not account for cavitation.

The image sequences of Fig. 8 show the internal flow field of the droplet and illuminate the various discontinuities using a numerical schlieren (i.e., exponential of the negative, normalized density gradient³²,³⁸). Image sequences are dis-
<figure><img src="image_8.png" /><figcaption>FIG. 6: (a) Radial bubble-wall evolution for a spherical bubble collapse with $p_\infty/p_b = 1427$ and $N = 100$ cells, with (b) a magnified view in the rebound region.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7: Effect of the pressure-relaxation rate on the growth and collapse of the gas phase within the cylindrical droplet.</figcaption></figure>

played for simulations with $\mu$ in the [1, 10, 100] range. For $\mu = 1$, the numerical schlieren images are overlaid with the predicted wavefronts computed from Eq. 1 (solid red lines), as well as the caustic associated to the singly-reflected wavefront given by Eq. 2 (dashed red line). This superposition reveals an excellent agreement between the simulation and the theory. The parametric equations of the transmitted shock and the k-th reflections being derived by assuming they propagate in an homogeneous phase, this agreement indicates the absence of bubble growth within the droplet. This is consistent with Fig. 7, where the dimensionless ratio $R/R_0$ does not increase as time proceeds. Note that both the theoretical and numerical singular points (S) well trace out the caustic associated to the first reflection ($k=2$).

Now considering the simulation results for $\mu = 10$ and

$\mu = 100$, Fig. 8 shows a very similar phenomenology. It however displays strong discrepancies with the simulation with $\mu = 1$, that the comparison of the theoretical wave pattern with the numerical schlieren, along with Fig. 7, evidences to be related to the growth of the gas phase within the droplet. First, the transmitted shock propagates downstream. Upon interaction with the droplet boundary, it reflects as an expansion wave. As previously described, this wavefront consists in two segment, $\mathcal{F}_{2,a}$ and $\mathcal{F}_{2,b}$ [see Fig. 8 at $\tilde{t} = 0.87$ for $\mu = 1$ and $\mu = 10$, and the schematics in Fig. 8(a)]. This expansion wave generates a low pressure region which, under some conditions, is likely to expand the gas phase within the droplet. When the expansion wave is strong enough to generate a bubble cloud, it modifies the canonical wave pattern observed in the simulation with $\mu = 1$. In this configuration, the forehead segment $\mathcal{F}_{2,a}$ initiates the growth of the gas phase which scatters the incoming subsequent segment $\mathcal{F}_{2,b}$ [see Fig. 8 from $\tilde{t} = 0.87$ for $\mu = 100$ and from $\tilde{t} = 1.21$ for $\mu = 10$, and the schematics in Fig. 8(a-b)]. As a result of this scattering process, only $\mathcal{F}_{2,a}$ remains and continues to propagate while expanding the gas phase within the droplet [see Fig. 8 at $\tilde{t} = 1.21$ for $\mu = 10$ and $\mu = 100$, and the schematics in Fig. 8(c)]. In this context and assuming an infinite pressure-relaxation rate, the volume of the bubble cloud corresponds to volume through which, segment $\mathcal{F}_{2,a}$ propagates. The complete singly-reflected wavefront ($\mathcal{F}_{2,a}$ and $\mathcal{F}_{2,b}$) is known to be a spherically converging wavefront whose intensity is amplified due to the focusing. The maximum intensity is reached at the cuspidal point of the caustic (C). When the wavefront meets C, it then spherically diverges by propagating upstream (Fig. 8 at $\tilde{t} = 1.56 - 1.90$ for $\mu = 1$). However, due to the scattering of $\mathcal{F}_{2,b}$, the energy collected at C decreases as $\mu$ increases. This is qualitatively visible in Fig. 8 at time $\tilde{t} = 1.56$, where a slight diverging wave is observed for $\mu = 10$ (light blue box), and no discontinuity is detected for $\mu = 100$. Note that the schlieren in the light blue box has been post-processed with a contrast en-
<figure><img src="image_10.png" /><figcaption>FIG. 8: Influence of $\mu$ on the internal flow field of a 22 mm-diameter cylindrical droplet interacting with a planar shock wave propagating at $M = 2.4$. (a-c) Schematics of the bubble formation and the scattering process.</figcaption></figure>

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

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

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

hancement algorithm to highlight the wavefront. In the case of $\mu = 10$, the diverging wavefront propagating upstream weakens as travelling and reflects again at the droplet boundary (Fig. 8 at $\tilde{t} = 1.90$). The resulting multiple wavefronts are fi-

nally invisible on the numerical schlieren and other computed flow fields. From time 1.56 to 2.56 (resp. 2.93), the simulation with $\mu = 10$ (resp. $\mu = 100$) does not exhibit any flow discontinuity related to the internal wave reflections. Between
time 2.56 and 2.76 (resp. 2.93 and 3.28), the bubble cloud col-
lapse has been completed in the case $\mu = 10$ (resp. $\mu = 100$)
and the collapse-induced shock (CiS) has been emitted (Fig. 8
at $\tilde{t} = 2.76$ for $\mu = 10$). Upon collapse, the bubble cloud does
not experience an additional growth phase, which allows the
reflection of the CiS at the rear side of the droplet to propagate
upstream without being scattered (Fig. 8 at $\tilde{t} = 2.93$ for $\mu = 10$
and $\tilde{t} = 3.28$ for $\mu = 100$). The CiS is a compression wave,
which implies that its reflection is an expansion wave. This
expansion is strong enough to generate a new bubble cloud
as clearly visible with the third discontinuities at time 3.28
for $\mu = 100$. This last phenomenon is consistent with Figs. 4
and 7, where a rebound in the time-dependent $R/R_0$-ratio is
observed between $\tilde{t} = 3.00$ and $\tilde{t} = 3.25$.

## 2. Pressure-relaxation rate calibration with experiment

The pressure field within the droplet has been previously
shown to be dependent on the pressure-relaxation rate (Sec-
tion IV B 1). We now determine the pressure-relaxation rate
by calibrating the simulation against the experimental data of
Sembian et al. We first compare and analyze three different
orders of magnitude of $\mu$. For $\mu$ equals to 1, 10 and 100,
Fig. 9 displays a comparison of the numerical solutions with
an experimental schlieren photograph at same initial condi-
tions and time. The numerics plot the colored volume fraction
of air (yellow-to-black colormap) overlaid with the grayscale
schlieren. Note that the numerical schlieren can only con-
tour the interface between the gas and the liquid phase, while
an experimental schlieren photograph, sensitive to the first
derivative in density, images the line-path integrated volume
of the cloud (dark area).

Test case $\mu = 1$ shows a very good agreement in both
shape and location of the internal wavefronts, but does not dis-
close the existence of a bubble cloud as evidenced by the vol-
ume fraction mapping and the lack of collapse-induced shock
(CiS). This suggests that $\mu = 1$ underestimates the experi-
mental relaxation rate. The internal wave structure reported
for $\mu = 10$ does not agree neither with the experiments nor
the theory, but however reveals a CiS generated by the previ-
ous growth and collapse of a bubble cloud. The absence of
the wavefront as predicted, and seen for $\mu = 1$, results from
the early development of the cloud that scattered the segment
$\mathcal{F}_{2,b}$. Discrepancies are also observed by comparing the cen-
ter of the experimental bubble cloud and the origin of the
CiS. These observations indicate that $\mu = 10$ overestimates
the experimental relaxation rate. The pressure-relaxation rate
$\mu = 100$ is eventually simulated. As expected, it also over-
estimates the experimental relaxation rate. The internal wave
structure does not agree with the experiments due to the obvi-
ous growth of the gas phase which has not yet completed its
collapse. Finally, comparing simulation results for $\mu$ equals to
1, 10 and 100 indicates that the experimental relaxation rate
should be between 1 and 10.

To identify the pressure-relaxation rate exhibiting the better
agreement with the experiments of Sembian et al., we com-
pared the center of the bubble cloud and the collapse time, de-

noted $x_c$ and $t_c$ respectively, between the experiments and the
numerical simulations. The experimental center of the cloud
is determined by detecting the cloud contour on the lower
halves of Fig. 9 using an edge detection algorithm, and com-
puting the center-of-mass. The collapse time is determined
based on the image sequence available in Figure 8 in Sembian
et al. from which we assume the collapse to occur between
frame (e) and (f). Using the position of the internal and ex-
ternal wavefronts (see Fig. 11), and knowing the size of the
droplet as well as the shock wave Mach number, the dimen-
sional time of frames (e) and (f) has been determined. We
denote $\Delta t_c$ the time interval between frames (e) and (f).

Figure 10 shows the functional dependency of $x_c$ and $t_c$
on the pressure-relaxation rate which is estimated using non-
linear least square fits of the form $a\mu^b + c$. The two plots
exhibit an asymptotic behavior as $\mu \to \infty$, which is consistent
with the analysis of Fig. 9.

Figure 12 plots the $(x_c, t_c)$-coordinates for various $\mu$ in the
$[3, \infty]$ range. The dark solid line corresponds to the experi-
mental $x_c$-coordinate, and the gray filled area to $\Delta t_c$. It appears
that the dimensionless cloud center coordinate $\tilde{x}_c$ linearly de-
pends on the collapse time $\tilde{t}_c$. The linear interpolation of the
$(x_c, t_c)$-coordinates intersects the experimental $x_c$ coordinate
on $\mu = 3.5$, and within the $\Delta t_c$ window.

Figure 13 shows the numerical results for $\mu = 3.5$ com-
pared with the experimental schlieren image. The upper half
of Fig. 13(a) plots the volume fraction of air with a numerical
schlieren, and the upper half of Fig. 13(b) maps the mixture
pressure field overlaid with the numerical schlieren. Fig. 13(a)
reports an excellent agreement both on the internal wave struc-
ture (dashed red line), and the size and location of the bubble
cloud (solid red line). Within the region $\mathcal{R}$ delimited by the
dashed red line, one can note the growth of a second bub-
ble cloud. It results from the interaction of the 3-rd reflected
transmitted shock ($\mathcal{F}_3$, red dashed line) with the CiS reflec-
tion. It is evidenced by Fig. 13(b) which shows two regions
$\mathcal{R}_1$ and $\mathcal{R}_2$ with opposite curvatures and propagation direc-
tions. Region $\mathcal{R}_1$ is driven by the downstream propagation of
the $\mathcal{F}_3$ wavefront, while region $\mathcal{R}_2$ is induced by the upstream
propagation of the CiS reflection. The interplay between $\mathcal{F}_3$
and the CiS reflection results in a low pressure region at the
origin of the growth of a second bubble cloud. A closer look
at the experimental schlieren image (see the magnified view
in Fig. 13(b)) shows discontinuities, highlighted with dashed-
dotted white lines, which align with the upstream contour of
region $\mathcal{R}_2$ (white dashed line). We infer these discontinuities
to be the CiS reflection. In light of the comparison displays
in Fig. 13, the numerical simulation very well agrees with the
experiment of Sembian et al.$^b$, thus validating the pressure-
relaxation rate $\mu = 3.5$.

A descriptive comparison of the numerical simulations with
the experimental observations for $\mu \in [1, 3.5, 10]$, and differ-
ent times, is given in Appendix A.

<footnote>$^b$ Note that Fig. 11 is plotted for $\mu = 3.5$.</footnote>
<figure><img src="image_14.png" /><figcaption>FIG. 9: Comparison of the internal droplet structure, at $\tilde{t} = 2.72$, between (top) numerical simulations at different $\mu$ and (bottom) experiments of Sembian *et al.* The upper-halves display the volume fraction of air (yellow-to-black colormap) overlaid with numerical schlieren images (white). The lower-halves disclose experimental schlieren images which, for $\mu = 1$, is compared with the theoretical predictions given by Eqs. 1 for $k \in [1, 7]$. Reproduced from S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis, "Plane shock wave interaction with a cylindrical water column", Phys. Fluids, 28, 056102, 2016, with the permission of AIP Publishing.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 10: Dependence of $\tilde{x}_c$ and $\tilde{t}_c$ coordinates on $\mu$. Both fits are non-linear least square fits of the form $a\mu^b + c$ with (a) $a = -0.15$, $b = -0.90$ and $c = 0.36$, and (b) $a = -2.61$, $b = -0.73$ and $c = 3.10$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 11: Interaction of a planar shock wave at $M = 2.4$ with a 22-mm diameter cylindrical water droplet. The upper half is a numerical schlieren and the lower-half is a schlieren photograph from Sembian *et al.* Reproduced from S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis, "Plane shock wave interaction with a cylindrical water column", Phys. Fluids, 28, 056102, 2016, with the permission of AIP Publishing. The agreement in the internal and external wave locations enables to calibrate the dimensional time of the experiments. Both the numerics and the experiment image the irregular Mach reflection over the droplet, which consists in the reflected shock wave, the Mach stem, the slip surface and the incident shock wave (not visible here) all connected by the triple point. They also image the third internal reflection, $\mathcal{F}_3$, of the transmitted shock.</figcaption></figure>

### 3. Validation and critical Mach number

Sembian *et al.* carried out experiments at $M = 2.4$ and $M = 1.75$. They observed the growth of a bubble cloud for

$M = 2.4$, but they did not capture it for $M = 1.75$. We thus simulate the experiments at $M = 1.75$ with $\mu = 3.5$ to assess the determined value of the pressure-relaxation rate. Figure 14 plots the dimensionless equivalent radius of the
<figure><img src="image_17.png" /><figcaption>FIG. 12: Identification of the pressure-relaxation rate corresponding to the experiment of Sembian *et al.* and based on the time and x-location of the cloud center.</figcaption></figure>

cloud for various shock wave numbers. It shows that the gas phase inside the droplet does not expand for $M = 1.75$, which agrees with the experimental observations ($\max(\tilde{R}) \sim 10^{-10}$). Running the simulation for $M = 1.75$, but with $\mu = \infty$, however shows the significant expansion of the gas phase, which is completely off the experimental observations. Hence, we validate the value $\mu = 3.5$. Note that similar results than for $(M, \mu) = (1.75, \infty)$ would also be expected for thermodynamically-consistent, mechanical-equilibrium or thermodynamical-equilibrium models. In addition, we also run parametric simulations for $M \in [1.75, 2.4]$ to evaluate the critical Mach number, $M_c$, from which cavitation occurs. From Fig. 14, we estimate $M_c$ to be in the [1.9,2.0] range where $\max(\tilde{R})$ varies from $\sim 10^{-9}$ to $\sim 10^{-8}$. For $M = 1.9$ (resp. $M = 2.0$), the minimum pressure behind the expansion wave $\mathcal{F}_2$, measured from its total reflection to complete focusing, ranges from $-0.37$ MPa to $-1.9$ MPa (resp. from $-0.5$ MPa to $-2.4$ MPa). This is consistent with the pressure thresholds for the heterogeneous cavitation reported in the literature and found to be between $-0.1$ MPa and $-1$ MPa³⁹.

### C. High-speed droplet impact

We previously determined and validated the pressure-relaxation rate for a water droplet impacted by a planar shock wave. To evaluate the dependence of $\mu$ on the material properties, we now simulate the high-speed impact of a cylindrical gelatin droplet on a solid substrate [Fig. 3(c)] and compare the results to the experiment of Field, Dear, and Ogren.

After parametric investigation on the pressure-relaxation rate and calibration on the experimental observations, the simulation which agrees best with the experiment uses $\mu = 0.5$. A comparison between the numerics and the experiment is shown in Fig. 15. The left-hand side of Fig. 15(a) plots the volume fraction of air, superposed with a numerical schlieren.

<figure><img src="image_18.png" /><figcaption>FIG. 13: Simulation results for $\mu = 3.5$ compared with the experimental schlieren image of Sembian *et al.* Reproduced from S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis, "Plane shock wave interaction with a cylindrical water column", Phys. Fluids, 28, 056102, 2016, with the permission of AIP Publishing. The upper half shows (a) the volume fraction of air with numerical schlieren and (b) the mixture pressure field superposed with the numerical schlieren.</figcaption></figure>

The only numerical schlieren is plotted in Fig. 15(b). Figure 15(c) sheds light on common features reported on both numerical and experimental images, used to validate the simulation. The dashed gray line corresponds to the contour of the volume fraction of air as drawn in Fig. 15(a) with the dashed white line. The red dashed lines are discontinuities extracted from the photograph, while the white dashed lines have been plotted from the numerical discontinuities. While discrepancies in the shape and the location of all these dis-
<figure><img src="image_19.png" /><figcaption>FIG. 14: Parametric simulations of the shock-droplet interaction for various shock wave Mach numbers ranging from 1.75 to 2.4. Labels indicate the (M, µ) combination used in the simulations.</figcaption></figure>

continuities are observed, numerics and experiments present a very similar pattern. Among others, a very good agreement is reported on the shape and location of the denser region of the bubble cloud, i.e. the centered dark region on the schlieren images. Finally, we plotted the contour of the bubble cloud computed from the numerical simulation [white dashed line in Fig. 15(d)]. Remember that the numerical schlieren only contours the cloud. The palm-like shape of the bubble cloud relatively well agrees with the experimental observation. The location and the curvature of the stem perfectly match the darker region in the experiments, while the top of the palm connects the streamers experimentally observed.

## V. LOCATION OF THE BUBBLE CLOUD

The location of the caustic's cuspidal point, $(x_\mathcal{C}, y_\mathcal{C})$ has been proven to be the highest density of shock-induced cavitation⁴⁰. This is because the density of rays tangent to the caustic, which gives a relative measure of the focusing strength over the caustic, is maximum at the intersection of the caustic and the droplet axis ($y = 0$)⁴. Consequently, previous research efforts on shock-induced cavitation within a cylindrical droplet assumed the cavitation bubble cloud to appear at the focus of the reflected wavefront $(x_\mathcal{C}, y_\mathcal{C})$²,¹³. However, a close examination of the bubble cloud center as seen in the experimental observations shows discrepancies between the theoretical $x_\mathcal{C}$ value and experimental measurements, errors of 42 ± 10% and 23 ± 5% for the shock-droplet interaction and the high-speed droplet impact, respectively. As clearly visible in Fig. 8, the center of the cloud on the x-axis does not agree with the cuspidal point of the caustic ($\mathcal{C}$). In the following, we denote $x_{c,n}$ the location of the cloud determined from the numerical simulations. The top graph in Fig. 16 plots

the $x_{c,n}/x_\mathcal{C}$ ratio over the pressure-relaxation rate $\mu$. It is obvious that the theoretical value $x_\mathcal{C}$ does not agree with the numerical value $x_{c,n}$. Note that, in the absence of the bubble cloud growth ($\mu = 1$), a very good match has nevertheless been reported between $x_\mathcal{C}$ and the focal point of the reflected wave simulated. As $\mu$ decreases, the $x_{c,n}$ converges towards $x_\mathcal{C}$, while an asymptotic behavior of the relative error between $x_{c,n}$ and $x_\mathcal{C}$ is reported when $\mu \to \infty$.

Three regimes of droplet shock-induced cavitation can be defined: (i) no-cavitation regime, (ii) convergence-driven cavitation regime, and (iii) the immediate cavitation regime. When the transmitted shock is not strong enough, the droplet never experiences cavitation and bubble growth. This is the no-cavitation regime. We note $\mathcal{I}_1$ the intensity of the transmitted shock in this regime. In the other extreme, the strongest transmitted shock instantaneously results in cavitation upon the reflection of the transmitted shock as an expansion wave. This corresponds to the immediate cavitation. We note $\mathcal{I}_2$ the intensity of the transmitted shock in this regime. However, for transmitted shock with intensity $\mathcal{I}_1 < \mathcal{I} < \mathcal{I}_2$, cavitation does not immediately occur upon reflection, but during the convergence of the expansion wave which induces its amplification. This corresponds to the transitional convergence-driven cavitation regime.

We denote, $x_t$, the location of the expansion wave $\mathcal{I}_2$, during its convergence, where the critical amplitude for the expansion wave to sufficiently expand the gas is reached. Locating the droplet center at the domain origin so that $(x, y) = (0, 0)$, the $x_t$-coordinate is bounded as $x_\mathcal{C} < x_t < R_d$. Note that $x_t = R_d$ is the condition for the immediate cavitation regime. The diameter of the bubble cloud measured on the x-axis is then given by $x_t - x_\mathcal{C}$. Assuming a volumetric collapse of the cloud, this implies that the cloud center $\Delta_x$ is given by

$$ \Delta_x = x_\mathcal{C} + \frac{x_t - x_\mathcal{C}}{2} \qquad (25) $$

The bottom graph in Fig. 16 plots the $x_{c,n}/\Delta_x$ ratio over the pressure-relaxation rate $\mu$, and shows a very good agreement between the bubble cloud center measured on the numerical simulation and the computed $\Delta_x$. Note that the second term in the right-hand side of the Eq. 25, $x_t - x_\mathcal{C}/2$, is a corrective term to Obreschkow et al.'s theory to estimate the location of the cavitation region in the convergence-driven cavitation regime. Eq. 25 agrees to locate the highest density of shock-induced cavitation at the caustic's cuspidal point ($\lim_{x_t \to x_\mathcal{C}} \Delta_x = x_\mathcal{C}$). In the immediate cavitation regime, $\Delta_x = x_\mathcal{C} + (R_d - x_\mathcal{C})/2$. In agreement with the numerical simulations and the experimental observations, and using Eq. 4 for $k = 2$, the location of the bubble cloud center during droplet shock-induced cavitation thus obeys

$$ \frac{n}{3n-1} \le \frac{\Delta_x}{R_d} \le \frac{1}{2} \left( 1 + \frac{n}{3n-1} \right) \qquad (26) $$

Note that the analysis on the cloud center, here shown for the configuration of Sembian et al., is also valid for the high-speed droplet impact configuration.
<figure><img src="image_20.png" /><figcaption>FIG. 15: Comparison of the numerical simulation (left) with the experiments of Field, Dear, and Ogren² (right). Reproduced from J. E. Field, J. P. Dear, and J. E. Ogren, "The effects of target compliance on liquid drop impact", J. Applied Phys., 65, 533–540, 1989, with the permission of AIP Publishing. The numerical image shows (a) the colored volume fraction of air overlaid with a schlieren image. (b) the numerical schlieren image alone. (c) is a magnified view of the bubble cloud observed on (b). (d) is a second magnified view of the bubble cloud showing, in white dashed line, the bubble cloud contour computed from the numerics.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIG. 16: Comparison between the bubble cloud center measured from the numerical simulations and theoretical predictions.</figcaption></figure>

## VI. CONCLUSION

In this paper, we introduce a multiphase numerical model, in velocity equilibrium, using a finite pressure-relaxation rate. In the context of heterogeneous cavitation, we demonstrate the finite formulation to be more suitable and effective to simulate shock-induced cavitation. Based on the shock-droplet interaction experiments of Sembian *et al.*, where cavitation bubbles within the droplet have been reported, we calibrate and validate the pressure-relaxation rate, which is found to be $\mu = 3.5$. A parametric investigation on $\mu$ shows the effect of the bubble cloud of the internal wave structure which, for a significant growth, scatters the incoming wavefronts resulting in their annihilation. For $\mu > 3.5$, the size, the location and the collapse time of the bubble cloud is overestimated as $\mu$ increases, and conversely, these parameters are underestimated when $\mu$, below 3.5, decreases. We eventually determined the critical shock Mach number $M_c$ from which, shock-induced cavitation is possible, to be $1.9 < M_c < 2.0$. This is consistent with Sembian *et al.* who observed cavitation at $M = 2.4$, while no bubbles has been imaged at $M = 1.75$.

Complementing the phenomenological analysis, based on the numerical results, with ray theory from geometrical acoustics, we discussed the theoretical location of the bubble cloud center. Although usually approximated to be located at the focal point of the singly reflected wavefront, we actually show
that it depends on the magnitude of the continuously ampli-
fying expansion wave, so that $\Delta_x \ge x_{\mathcal{E}}|_{k=2}$. The correction of
the $\Delta_x$ is done by adding the $(x_t - x_{\mathcal{E}})/2$ term to the $\mathcal{F}_2$'s focal
point location $x_{\mathcal{E}}$. Future work should address the analytical
determination of the $x_t$ location, which is here numerically de-
termined.

To evaluate the sensitivity of the pressure-relaxation rate on
the material properties where the cavitation occurs, we then
simulated the experiment of Field, Dear, and Ogren consisting
in the impact of a solid surface onto a gelatin droplet. Field,
Dear, and Ogren observed the cavitation of bubbles under the
internal reflection of the water hammer shock. Our calibration
procedure identified $\mu = 0.5$ to best match the experimental
observations. The droplet dependency on $\mu$ is consistent with
the previous observations in the simulation of Sembian et al.'s
configuration. Simulating the spherical collapse of a bubble in
a free field and comparing the bubble dynamics to the analyt-
ical solution of the Keller-Miksis equation, the $[0.1, \infty]$ range
has been found to very well agree with theory. Encouragingly,
the two pressure-relaxation rates determined fits in this $\mu$ in-
terval. Note that, in this work, the calibrated values of $\mu$ are
only valid for the two configurations we have been concerned
with and for an initial gas volume fraction $\alpha_g = 10^{-6}$. Ef-
fects of $\alpha_g$ on the phenomenology of the shock-induced cav-
itation require additional research efforts based on comple-
mentary experiments, which are out of the scope of this paper.
However, to shed light on the possible influence of $\alpha_g$ on the
physics and the calibrated $\mu$, additional simulations have been
run and discussed in Appendix B.

Future works would be to integrate phase change in the
numerical model which would require to calibrate again the
pressure-relaxation rate along with the chemical-potential re-
laxation rate. A closer examination on the dependence of $\mu$
on the material properties should also be considered. In this
work, we infer the difference in the $\mu$ values between the two
configurations investigated to be related to the changes in the
properties of the material constituting the droplet, i.e. water
versus gelatin. The present modelling and existing experimen-
tal datasets do not allow for the identification of the governing
parameters, while the viscoelasticity, the spatial arrangement
of molecules, the concentration and size of nuclei, or the equa-
tions of state should be considered. Additional experiments,
varying the material properties and the regimes, should be
conducted to independently investigate the influence of these
parameters on the cavitation dynamics.

ACKNOWLEDGMENTS

Authors thank Fabien Petitpas and Nicolas Favrie for their
participation in the determination of the $\mu_{k,j}$ formulation. Au-
thors acknowledge the financial support from the ETH Zurich
Postdoctoral Fellowship program.

CONFLICT OF INTEREST

The authors have no conflicts to disclose.

DATA AVAILABILITY

The data that support the findings of this study are available
from the corresponding author upon reasonable request.

Appendix A: Time sequences of shock-droplet interaction

Figure 17 shows image sequences for the interaction of a 22-mm-diameter cylindrical droplet with a planar shock wave propagating at Mach 2.4. The upper halves are the colored volume fraction of air (yellow-to-black colormap) overlaid with a numerical schlieren in grayscale. The lower-halves are the experimental schlieren photographs from Sembian et al. The comparison is disclosed for $\mu$ equals 1, 3.5 and 10. On frames (a) and (b), the three simulations are in excellent agreement with the experiments. The transmitted shock and the expansion wave both agree on time and space. The simulations do not exhibit the expansion of the gas phase within the droplet, which is in accordance with the absence of bubble cloud on the photographs. On frame (c), the three numerical wave patterns are still in line with the experiments. However, discrepancies on the volume fraction of air arise for $\mu = 10$ (red box), which shows the growth of the gas phase downstream the expansion wave, while no bubble cloud is imaged in the experiment. On frame (d), only $\mu = 1$ and $\mu = 3.5$ well simulate the diverging expansion wave, which is invisible in the simulations with $\mu = 10$. Simulations with $\mu = 1$ and $\mu = 3.5$ also show an horseshoe-like shape bubble cloud behind the expansion wave. It is difficult to discuss the validity of this cloud against the experiment, as the schlieren photograph exhibits an intricate structure inside the horseshoe shape which complicates the analysis of the image. However, it is obvious that the size of the cloud simulated with $\mu = 10$ is significantly overestimated. On frame (e), experiments show a bubble cloud (contoured with the red dashed line). In the simulation with $\mu = 1$, no cloud is visible. In addition, the intensity of the numerical schlieren is the same for the various internal reflections, while the $\mathcal{F}_2$ wavefront is significantly stronger on the experimental schlieren. Note that these observations holds for frame (f). The wave pattern and the bubble cloud observed in the simulation with $\mu = 3.5$ are however in a very good agreement, which is also true for frame (f). Although the simulation with $\mu = 10$ significantly fails to reproduce the wave pattern, we note that a relatively good match is reported on the bubble cloud. However, on frame (f) and still for $\mu = 10$, the CiS reveals a late collapse of the bubble when comparing with the discontinuities seen in the experiment [see the magnified view in Fig. 13(b)].

Appendix B: Effect of the initial gas volume fraction

In this work, we are concerned with the growth and coll-
pase of the gas phase involved in the liquid-gas mixture ini-
tially constituting the droplet. This requires to set the initial
gas volume fraction, $\alpha_g$ within the droplet. In the experiments
<figure><img src="image_22.png" /><figcaption>FIG. 17: (a-f) Comparison of the internal droplet structure between (top) numerical simulations at different $\mu$ and (bottom) the experiment of Sembian *et al.* The upper halves display the volume fraction of air (yellow-to-black colormap) overlaid with numerical schlieren images (white). The lower halves disclose experimental schlieren images. Reproduced from S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis, “Plane shock wave interaction with a cylindrical water column”, Phys. Fluids, 28, 056102, 2016, with the permission of AIP Publishing. The dashed red lines on frames (e) contour the bubble cloud as indicated in Sembian *et al*.</figcaption></figure>

of Sembian *et al.* and Field, Dear, and Ogren, $\alpha_g$, is an unknown parameter. We therefore estimated $\alpha_g$ to be equal to $10^{-6}$, which is assumed to be representative of purified water. The calibrated pressure-relaxation rates for the two configurations are thus only valid for simulations with $\alpha_g = 10^{-6}$. To assess the sensitivity of the pressure-relaxation rate on the initial gas volume fraction, we run additional simulations, corresponding to the experiment of Sembian *et al.*, with $\alpha_g = 10^{-5}$ and $\alpha_g = 10^{-7}$. Figure 18 plots the variation of the volume of the gas phase within the droplet as time proceeds for various $(\alpha_g, \mu)$ combinations. Comparing the results for the calibrated $\mu = 3.50$ with $\alpha_g = 10^{-6}$, relative errors between 40% and 60% are reported on the maximal volume, which remain relatively low when compared to the errors resulting from an infinite pressure-relaxation rate (≈650%). The relative errors on the numerical location of the bubble cloud center with respect to the experiments are found to be ≈10% for both $\alpha_g = 10^{-5}$ and $\alpha_g = 10^{-7}$ with $\mu = 3.50$. We then calibrated again the pressure-relaxation rates for $\alpha_g = 10^{-5}$ and $\alpha_g = 10^{-7}$ and found $\mu = 2.50$ and $\mu = 4.60$, respectively. It appears that changing $\alpha_g$ by one order of magnitude (lower or higher) does not change the order of magnitude of the pressure-relaxation rate.

<figure><img src="image_23.png" /><figcaption>FIG. 18: Variation of the volume of the gas phase within the droplet as time proceeds for various initial gas volume fraction $\alpha_g$ and pressure-relaxation rate $\mu$.</figcaption></figure>

<footnote n="1">¹S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis, “Plane shock wave interaction with a cylindrical water column,” Phys. Fluids **28**, 056102 (2016).</footnote>

<footnote n="2">²J. E. Field, J. P. Dear, and J. E. Ogren, “The effects of target compliance on liquid drop impact,” J. Applied Phys. **65**, 533–540 (1989).</footnote>

<footnote n="3">³H. Chen and S. M. Liang, “Flow visualization of shock/water column interactions,” Shock Waves **17**, 309–321 (2008).</footnote>

<footnote n="4">⁴L. Biasiori-Poulanges and H. El-Rabii, “Shock-induced cavitation and wavefront analysis inside a water droplet,” Phys. Fluids **33**, 097104 (2021).</footnote>
⁵S. Balibar and F. Caupin, "Metastable liquids," J. Phys.: Condensed Matter **15**, S75 (2002).

⁶K. Ando, *Effects of polydispersity in bubbly flows*, Ph.D. thesis, California Institute of Technology (2010).

⁷J. C. Meng and T. Colonius, "Numerical simulations of the early stages of high-speed droplet breakup," Shock waves **25**, 399-414 (2015).

⁸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**, N3069 (2013).

⁹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**, N117 (2015).

¹⁰S. R. G. Avila and C.-D. Ohl, "Fragmentation of acoustically levitating droplets by laser-induced cavitation bubbles," J. Fluid Mech. **805**, 551-576 (2016).

¹¹L. Biasiori-Poulanges and H. El-Rabii, "Multimodal imaging for intradroplet gas-cavity observation during droplet fragmentation," Optics Letters **45**, 3091-3094 (2020).

¹²K. K. Haller, Y. Ventikos, D. Poulikakos, and P. Monkewitz, "Computational study of high-speed liquid droplet impact," J. Applied Phys. **92**, 2821-2828 (2002).

¹³T. Kondo and K. Ando, "One-way-coupling simulation of cavitation accompanied by high-speed droplet impact," Phys. Fluids **28**, 033303 (2016).

¹⁴L. Biasiori-Poulanges, G. T. Bokman, E. Baumann, and O. Supponen, "Dynamics of a shocked bubble-encapsulated droplet," Applied Phys. Letters **120**, 260601 (2022).

¹⁵J. C. Fisher, "The fracture of liquids," J. Applied Phys. **19**, 1062-1067 (1948).

¹⁶R. E. Apfel, "Acoustic cavitation prediction," J. Acoust. Soc. Am. **69**, 1624-1633 (1981).

¹⁷N. Kyriazis, P. Koukouvinis, and M. Gavaises, "Modelling cavitation during drop impact on solid surfaces," Adv. Colloid Interface Sc. **260**, 46-64 (2018).

¹⁸K. Schmidmayer, J. Cazé, F. Petitpas, E. Daniel, and N. Favrie, "Modelling interactions between waves and diffused interfaces," Int. J. Numer. Meth. Fl. x, x (2022).

¹⁹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).

²⁰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).

²¹D. A. Drew and S. L. Passman, *Theory of Multicomponent fluids*, Applied Mathematical Sciences, Vol. 135 (Springer, New York, 1998).

²²M. R. Baer and J. W. Nunziato, "A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials," Int. J. of Multiphase Flow **12**, 861-889 (1986).

²³R. Saurel and R. Abgrall, "A multiphase Godunov method for compressible multifluid and multiphase flows," J. Comp. Phys. **150**, 425-467 (1999).

²⁴A. Kapila, R. Menikoff, J. Bdzil, S. Son, and D. Stewart, "Two-phase modeling of DDT in granular materials: Reduced equations," Phys. Fluids **13**, 3002-3024 (2001).

²⁵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).

²⁶E. F. Toro, *Riemann solvers and numerical methods for fluid dynamics* (Springer Verlag, Berlin, 1997).

²⁷B. Van Leer, "Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow," J. Comp. Phys. **23**, 263-275 (1977).

²⁸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).

²⁹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).

³⁰T. Trummler, S. H. Bryngelson, K. Schmidmayer, S. J. Schmidt, T. Colonius, and N. A. Adams, "Near-surface dynamics of a gas bubble collapsing above a crevice," J. Fluid Mech. **899**, A16 (2020).

³¹K. Schmidmayer, S. H. Bryngelson, and T. Colonius, "An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics," J. Comp. Phys. **402**, 109080 (2020).

³²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).

³³K. Schmidmayer, F. Petitpas, E. Daniel, N. Favrie, and S. L. Gavrilyuk, "A model and numerical method for compressible flows with capillary effects," J. Comp. Phys. **334**, 468-496 (2017).

³⁴J. B. Keller and M. Miksis, "Bubble oscillations of large amplitude," J. Acoust. Soc. Am. **68**, 628-633 (1980).

³⁵V. Coralic and T. Colonius, "Finite-volume WENO scheme for viscous compressible multicomponent flows," J. Comp. Phys. **274**, 95-121 (2014).

³⁶C. E. Brennen, *Cavitation and bubble dynamics* (Oxford University Press, 1995).

³⁷D. Fuster, C. Dopazo, and G. Hauke, "Liquid compressibility effects during the collapse of a single cavitating bubble," J. Acoust. Soc. Am. **129**, 122-131 (2011).

³⁸J. J. Quirk and S. Karni, "On the dynamics of a shock-bubble interaction," J. Fluid Mech. **318**, 129-163 (1996).

³⁹F. Caupin and E. Herbert, "Cavitation in water: a review," Comptes Rendus Physique **7**, 1000-1017 (2006).

⁴⁰D. Obreschkow, N. Dorsaz, P. Kobel, A. de Bosset, M. Tinguely, J. Field, and M. Farhat, "Confined shocks inside isolated liquid volumes: A new path of erosion?" Phys. Fluids **23**, 101702 (2011).