L Biasiori-Poulanges

K Schmidmayer
email: kevin.schmidmayer@inria.fr

Droplet shock-induced cavitation

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 µ = 3.5 and µ = 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 3 , 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 4 . 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 [START_REF] Balibar | Metastable liquids[END_REF] . 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 applia) https://kevinschmidmayer.github.io/ cations, as a desired or adverse effect, ranging from raindrop impact on aircraft [START_REF] Ando | Effects of polydispersity in bubbly flows[END_REF] , to combustion and detonation of multiphase mixtures [START_REF] Meng | Numerical simulations of the early stages of high-speed droplet breakup[END_REF] , through ink-jet printing or liquid jet-based physical cleaning [START_REF] Okorn-Schmidt | Particle cleaning technologies to meet advanced semiconductor device process requirements[END_REF][START_REF] Tatekura | Observation of water-droplet impacts with velocities of o (10 m/s) and subsequent flow field[END_REF] , 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 bubblecompounded 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 [START_REF] Avila | Fragmentation of acoustically levitating droplets by laser-induced cavitation bubbles[END_REF][START_REF] Biasiori-Poulanges | Multimodal imaging for intradroplet gas-cavity observation during droplet fragmentation[END_REF] , 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 [START_REF] Biasiori-Poulanges | Multimodal imaging for intradroplet gas-cavity observation during droplet fragmentation[END_REF] . 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 shockinduced 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 1,4, [START_REF] Haller | Computational study of high-speed liquid droplet impact[END_REF][START_REF] Kondo | One-way-coupling simulation of cavitation accompanied by high-speed droplet impact[END_REF][START_REF] Biasiori-Poulanges | Dynamics of a shocked bubble-encapsulated droplet[END_REF] . 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) [START_REF] Fisher | The fracture of liquids[END_REF] , or the Blake threshold pressure [START_REF] Apfel | Acoustic cavitation prediction[END_REF] . Recently, Kyriazis, Koukouvinis, and Gavaises simulated the experiment of Field, Dear, and Ogren, i.e. high-speed droplet impact, 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, µ, to address the over-expansion of the gas phase as previously reported for an infinite µ. The finite pressure-relaxation rate is defined on the ]0, ∞] range. Shock-induced cavitation primary resulting from heterogeneous cavitation, we do not account for phase change. The droplet containing preexisting 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 µ range, and compare the results with the solution of the Keller-Miksis equation to eliminate µ 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 µ on the shock-induced cavitation, we calibrate the finite pressurerelaxation 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 µ on the material properties of the liquid mixture constituting the droplet, by simulating the highspeed 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-toair 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 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.

and its multiple internal reflections, which read

x M =[c l t -nR d (1 -cos α) -2(k -1)R d cos θ ] cos(γ k -θ ) -R d cos(γ k -2θ ), (1a)
y M =[c l t -nR d (1 -cos α) -2(k -1)R d cos θ ] sin(γ k -θ ) -R d sin(γ k -2θ ), (1b)
where k = 1 corresponds to the transmitted shock and k = 2, 3, ... to the successive internal reflections, n is the water-toair sound speed ratio c l /c g and γ k = 2kθα -(k -1) π. R d is the droplet radius. We denote F k the wavefront associated to the k-th reflection. The incident and refraction angles, α and θ , are related by the fundamental law of refraction, sin θ = n sin α.

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 wavefront, which are given by

x caustic =R d f (α) cos γ k + R d [ f (α) -1] cos(γ k -2θ ), (2a) y caustic = R d f (α) sin γ k + R d [ f (α) -1] sin(γ k -2θ ), (2b)
where The red-to-blue colored segments are the singly-reflected wavefronts as time proceeds.

f (α) = 1 2 2n 2 (k -1) sin 2α -sin 2θ n 2 (2k -1) sin 2α -sin 2θ . ( 3
The cuspidal point of the caustic is located on the droplet axis y C = 0, while the x-coordinate has been found to be

x C = (-1) k n (2k -1)n -1 R d . (4)

C. Problem dimensions

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

M = U s c , We = ρU 2 d 0 σ and Re = ρUd 0 µ . (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), ρ is the density of the post-shocked gas (resp. droplet density), U is the post-shocked gas velocity (resp. impact velocity), µ is the dynamic viscosity of the gas (resp. dynamic viscosity of the liquid), σ 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, inviscid 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

L = L d 0 and T = T c l d 0 , (6)
where (• ) denotes a non-dimensional quantity.

III. NUMERICAL MODELLING

We use herein a slightly modified version of the modelling proposed by Schmidmayer et al. [START_REF] Schmidmayer | Modelling interactions between waves and diffused interfaces[END_REF] 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 (righthand 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

∂ α k ∂t + u • ∇α k = δ p k , ∂ α k ρ k ∂t + ∇ • (α k ρ k u) = 0, ∂ ρu ∂t + ∇ • (ρu ⊗ u + pI) = 0, ∂ α k ρ k e k ∂t + ∇ • (α k ρ k e k u) + α k p k ∇ • u = -p I δ p k , (7)
where α k , ρ 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

ρ = N ∑ k=1 α k ρ k and p = N ∑ k=1 α k p k , (8)
while the mixture total energy is

E = e + 1 2 u 2 , ( 9
)
where e is the mixture specific internal energy

e = N ∑ k=1 Y k e k (ρ k , p k ) . (10)
In (10), e k (ρ k , p k ) is defined via an equation of state (EOS) and Y k are the mass fractions

Y k = α k ρ k ρ . (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 = ρ g (γ g -1)(e g -e g,ref ), (12)
and the liquid is modeled by the stiffened-gas (SG) EOS

p l = ρ l (γ l -1)(e l -e l,ref ) -γ l π ∞ , (13)
where γ, e ref and π ∞ are model parameters [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF] . The interfacial pressure is defined as

p I = ∑ N k p k ∑ N j =k z j (N -1) ∑ N k z k , (14)
where z k = ρ 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

∂ ρE ∂t + ∇ • [(ρE + p) u] = 0. ( 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 = N ∑ k=1 Y k c 2 k , (16)
which is found to be in agreement with previously reported expression [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] .

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 ± c, where u is the velocity in the x-direction.

B. Expression of δ p k

For the pressure-relaxation terms between the phases, δ p k reads under its general form

δ p k = N ∑ j =k µ k, j (p k -p j ) , (17)
where the relaxation coefficients µ 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 [START_REF] Drew | Theory of Multicomponent fluids[END_REF] , can be expressed under different forms.

In most if not all the literature, e.g. Schmidmayer et al. [START_REF] Schmidmayer | Modelling interactions between waves and diffused interfaces[END_REF] , Saurel, Petitpas, and Berry 20 , Baer and Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] , Saurel and Abgrall [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] , the relaxation coefficients are taken as unique and constant for all interactions, i.e. µ k, j = µ.

Herein, we propose to use a different approach and to express them as µ k, j = α k α j µ, where µ is a constant parameter. This leads to

δ p k = µα k N ∑ j =k α j (p k -p j ) . (18)
First, this expression is consistent for N phases, meaning ∑ k δ 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 volumefraction configuration within the computational cell. Note that µ is a finite parameter which can be selected in the ]0, ∞] 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 finitevolume 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 [START_REF] Schmidmayer | Modelling interactions between waves and diffused interfaces[END_REF] .

The relaxation terms (system of ordinary differential equations) are integrated with a first-order, explicit, Euler scheme with time-step subdivisions [START_REF] Schmidmayer | Modelling interactions between waves and diffused interfaces[END_REF] . 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 [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] 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 (µ = ∞), the model converges to the mechanical-equilibrium model of Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF] and the effective mixture speed of sound matches Wood's

1 ρc 2 w = N ∑ k=1 α k ρ k c 2 k . (19)
A second-order-accurate MUSCL scheme with two-step time integration is used [START_REF] Schmidmayer | ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows[END_REF] , where the first step is a predictor step for the second and the usual piece-wise linear MUSCL reconstruction [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] with the monotonized central (MC) [START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow[END_REF] 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 [START_REF] Schmidmayer | Adaptive Mesh Refinement algorithm based on dual trees for cells and faces for multiphase compressible flows[END_REF] . The cell i is refined when the following criterion is fulfilled

|X Nb(i, j) -X i | min(X Nb(i, j) -X i ) > ε, ( 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 ε = 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 [START_REF] Schmidmayer | ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows[END_REF] , 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., [START_REF] Schmidmayer | Modelling interactions between waves and diffused interfaces[END_REF][START_REF] Dorschner | On the formation and recurrent shedding of ligaments in droplet aerobreakup[END_REF][START_REF] Trummler | Near-surface dynamics of a gas bubble collapsing above a crevice[END_REF][START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF][START_REF] Pishchalnikov | High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone[END_REF][START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF] ).

D. Computational setup

The determination of the pressure-relaxation rate µ is a two-step approach. The first step consists in reducing the µrange by simulating the spherical collapse of a bubble in a free field and determined µ values that agree with the predicted bubble dynamics given by the Keller-Miksis equation [START_REF] Keller | Bubble oscillations of large amplitude[END_REF] . The second step consists in calibrating µ 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 µ 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 µ 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.

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 onedimensional (1D) domain of 3 mm long is used with spherical axi-symmetry to mimic a three-dimensional bubble [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF] . The domain consists in high-pressure water at p ∞ = 50 atm with density ρ l = 1000 kg/m 3 . 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 ρ b = 0.027 kg/m 3 , 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. and Field, Dear, and Ogren, respectively.

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 × 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 3 and a 1 atm pressure. The water has a density of 1028 kg/m 3 and is modelled using the SG EOS 19,a with γ = 2.35 and π ∞ = 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).

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 × 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 [START_REF] Kondo | One-way-coupling simulation of cavitation accompanied by high-speed droplet impact[END_REF][START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF] with γ = 6.72 and π ∞ = 3.70 × 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 17 .

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 µ = 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 "shockinduced 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 = 1 π V d α g dV 1/2 , (21)
a 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.

where V d is the droplet volume. AMR level 0 AMR level 1 AMR level 2 AMR level 3 FIG. 4: Growth and collapse of the gas phase within the droplet at four grid resolutions. The pressure-relaxation rate is µ = 10. The droplet is initially resolved by 100 cells per diameter.

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 ε as

ε = 1 N t N t ∑ i=0 R(t i ) -R ref (t i ) R ref (t i ) , (22)
where N t is the number of time steps in the temporal window t ∈ [0, 3.25], R(t i ) is the bubble radius at time t i of our simulations and R 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 a infinite range of pressure-relaxation rate. As a first approach, it is instructive to evaluate the sensitivity of µ by simulating the spherical collapse of a bubble in free field [Fig. 3(a)] for various orders of magnitude of µ, and compare the bubble response to the analytical solution of the Keller-Miksis equation [START_REF] Keller | Bubble oscillations of large amplitude[END_REF] ; the compressible form of the Rayleigh-Plesset equa-

10 2 2 × 10 2 3 × 10 2 4 × 10 2 N d0 10 -1 10 0 ε ∝ 1/N 1.7 d0 Numerics ∝ 1/N 2 d0
FIG. 5: Convergence order of the numerical method. The discrete L 2 -error, ε, representative of the convergence rate is given by Eq. 22.

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 borned 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 = N ∑ i=1 α g,i V c,i , (23)
where N is the number of grid cells, and α 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.915R 0 ρ l p ∞ , (24)
which is the nominal total collapse time, i.e. the time required for the bubble to complete its collapse [START_REF] Brennen | Cavitation and bubble dynamics[END_REF] . 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 [START_REF] Fuster | Liquid compressibility effects during the collapse of a single cavitating bubble[END_REF] . 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. Figure 6 and Table II show that pressure-relaxation rates 1 µ ∞ are in satisfying agreement with theory, while underestimating the minimum radius up to 19%. The best agreement with theory is given for µ = 1, where the relative error ε r on the minimum bubble radius, with respect to the Keller-Miksis solution, is 10%.

Note that min(R b /R 0 )| µ=10 ≈ min(R b /R 0 )| µ=∞ .
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 µ 0.1, which significantly overestimate the minimum bubble radius.

From this first sensitivity analysis on µ, the range of pressure-relaxation rate can be restricted to ]0.1, ∞]. Note that this range can be extended when refining the mesh size, so that µ ∈]a, ∞] with a → 0 for an infinitely small mesh size. The µ-range ]0.1, ∞] is here valid for a reasonable resolution of 100 cells/diameter. In the following, we calibrate µ upon experimental observations which is therefore expected to be in the µ-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 µ 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.

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 ∞, 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 µ. The maximum radius asymptotically increases with µ, so that max(R/R 0 )| µ=100 ≈ 2 max(R/R 0 )| µ=10 and R (µ = 100) ≈ R (µ = ∞). The simulation thus converges to an infinite pressure-relaxation rate modelling as µ increases. Conversely, smaller R are reported for lower values of µ. Numerical results computed for µ = 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 [START_REF] Pishchalnikov | High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone[END_REF][START_REF] Quirk | On the dynamics of a shock-bubble interaction[END_REF] ). Image sequences are dis- played for simulations with µ in the [1, 10, 100] range. For µ = 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 µ = 10 and µ = 100, Fig. 8 shows a very similar phenomenology. It however displays strong discrepancies with the simulation with µ = 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, F 2,a and F 2,b [see Fig. 8 at t = 0.87 for µ = 1 and µ = 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 µ = 1. In this configuration, the forehead segment F 2,a initiates the growth of the gas phase which scatters the incoming subsequent segment F 2,b [see Fig. 8 from t = 0.87 for µ = 100 and from t = 1.21 for µ = 10, and the schematics in Fig. 8(a-b)]. As a result of this scattering process, only F 2,a remains and continues to propagate while expanding the gas phase within the droplet [see Fig. 8 at t = 1.21 for µ = 10 and µ = 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 F 2,a propagates. The complete singly-reflected wavefront (F 2,a and 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 t = 1.56 -1.90 for µ = 1). However, due to the scattering of F 2,b , the energy collected at C decreases as µ increases. This is qualitatively visible in Fig. 8 at time t = 1.56, where a slight diverging wave is observed for µ = 10 (light blue box), and no discontinuity is detected for µ = 100. Note that the schlieren in the light blue box has been post-processed with a contrast en- hancement algorithm to highlight the wavefront. In the case of µ = 10, the diverging wavefront propagating upstream weakens as travelling and reflects again at the droplet boundary (Fig. 8 at 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 µ = 10 (resp. µ = 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 collapse has been completed in the case µ = 10 (resp. µ = 100) and the collapse-induced shock (CiS) has been emitted (Fig. 8 at t = 2.76 for µ = 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 t = 2.93 for µ = 10 and t = 3.28 for µ = 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 µ = 100. This last phenomenon is consistent with Figs. 4 and7, where a rebound in the time-dependent R/R 0 -ratio is observed between t = 3.00 and t = 3.25.

Pressure-relaxation rate calibration with experiment

The pressure field within the droplet has been previously shown to be dependent on the pressure-relaxation rate (Section 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 µ. For µ equals to 1, 10 and 100, Fig. 9 displays a comparison of the numerical solutions with an experimental schlieren photograph at same initial conditions 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 contour 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 µ = 1 shows a very good agreement in both shape and location of the internal wavefronts, but does not disclose the existence of a bubble cloud as evidenced by the volume fraction mapping and the lack of collapse-induced shock (CiS). This suggests that µ = 1 underestimates the experimental relaxation rate. The internal wave structure reported for µ = 10 does not agree neither with the experiments nor the theory, but however reveals a CiS generated by the previous growth and collapse of a bubble cloud. The absence of the wavefront as predicted, and seen for µ = 1, results from the early development of the cloud that scattered the segment F 2,b . Discrepancies are also observed by comparing the center of the experimental bubble cloud and the origin of the CiS. These observations indicate that µ = 10 overestimates the experimental relaxation rate. The pressure-relaxation rate µ = 100 is eventually simulated. As expected, it also overestimates the experimental relaxation rate. The internal wave structure does not agree with the experiments due to the obvious growth of the gas phase which has not yet completed its collapse. Finally, comparing simulation results for µ 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 compared 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 computing 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 external wavefronts (see Fig. 11), and knowing the size of the droplet as well as the shock wave Mach number, the dimensional time of frames (e) and (f) has been determined. We denote ∆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 nonlinear least square fits of the form aµ b + c. The two plots exhibit an asymptotic behavior as µ → ∞, which is consistent with the analysis of Fig. 9.

Figure 12 plots the (x c ,t c )-coordinates for various µ in the [3, ∞] range. The dark solid line corresponds to the experimental x c -coordinate, and the gray filled area to ∆t c . It appears that the dimensionless cloud center coordinate xc linearly depends on the collapse time tc . The linear interpolation of the (x c ,t c )-coordinates intersects the experimental x c coordinate on µ = 3.5, and within the ∆t c window.

Figure 13 shows the numerical results for µ = 3.5 compared 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 structure (dashed red line), and the size and location of the bubble cloud (solid red line). Within the region R delimited by the dashed red line, one can note the growth of a second bubble cloud. It results from the interaction of the 3-rd reflected transmitted shock (F 3 , red dashed line) with the CiS reflection. It is evidenced by Fig. 13(b) which shows two regions R 1 and R 2 with opposite curvatures and propagation directions. Region R 1 is driven by the downstream propagation of the F 3 wavefront, while region R 2 is induced by the upstream propagation of the CiS reflection. The interplay between 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 dasheddotted white lines, which align with the upstream contour of region 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 pressurerelaxation rate µ = 3.5.

A descriptive comparison of the numerical simulations with the experimental observations for µ ∈ [1, 3.5, 10], and different times, is given in Appendix A.

2 R R w x L X l Y 8 T R O I o M 9 8 j 2 W c O j p r o = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o M Q L 2 F X g n o M e v E Y 0 T w g W c L s p D c Z M j u 7 z M w K I e Q T v H h Q x K t f 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G t z O / 9 Y R K 8 1 g + m n G C f k Q H k o e c U W O l h z I 9 7 x V L b s W d g 6 w S L y M l y F D v F b + 6 / Z i l E U r D B N W 6 4 7 m J 8 S d U G c 4 E T g v d V G N C 2 Y g O s G O p p B F q f z I / d U r O r N I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b + 5 3 V S E 1 7 7 E y 6 T 1 K B k i 0 V h K o i J y e x v 0 u c K m R F j S y h T 3 N 5 K 2 J A q y o x N p 2 B D 8 J Z f X i X N i 4 p 3 W a n e V 0 u 1 m y y O P J z A K Z T B g y u o w R 3 U o Q E M B v A M r / D m C O f F e X c
+ F q 0 5 J 5 s 5 h j 9 w P n 8 A i 5 m N U w = = < / l a t e x i t >

q U H J F o v C V B A T k 9 n f p M 8 V M i P G l l C m u L 2 V s C F V l B m b T s G G 4 C 2 / v E q a F x X v s l K 9 r 5 Z q N 1 k c e T i B U y i D B 1 d Q g z u o Q w M Y D O A Z X u
J I x H p T k A M E 1 y x J n A Q r B N r R m Q g W D s Y 3 0 7 9 9 h P T h k f q A d K Y + Z I M F Q 8 5
J W C l x 8 p 5 j 6 s Q 0 k q / V H a r 7 g x 4 m X g 5 K a M c j X 7 p q z e I a C K Z A i q I M V 3 P j c H P i A Z O B Z s U e 4 l h M a F j M m R d S x W R z P j Z 7 O A J P r X K A I e R t q U A z 9 T f E x m R x q Q y s J 2 S w M g s e l P x P 6 cloud for various shock wave numbers. It shows that the phase inside the droplet does not expand for M = 1.75, which agrees with the experimental observations (max( R) ∼ 10 -10 ). Running the simulation for M = 1.75, but with µ = ∞, however shows the significant expansion of the gas phase, which is completely off the experimental observations. Hence, we validate the value µ = 3.5. Note that similar results than for (M, µ) = (1.75, ∞) would also be expected for thermodynamically-consistent, mechanical-equilibrium or thermodynamical-equilibrium models. In addition, we also run parametric simulations for M ∈ [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( R) varies from ∼ 10 -9 to ∼ 10 -8 . For M = 1.9 (resp. M = 2.0), the minimum pressure behind the expansion wave 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 39 .

+ b Q H j t Z 1 z F C T B F 5 4 v C R G C I 8 P R 7 P O C a U R C p J Y R q b m / F d E Q 0 o W A z K t o Q v M W X l 0 m r V v U u q x f 3 t X L 9 J o + j g I 7 R C T p D H r p C d X S H G q i J K J L o G b 2 i N 0 c 7 L 8 6 7 8 z F v X X

C. High-speed droplet impact

We previously determined and validated the pressurerelaxation rate for a water droplet impacted by a planar shock wave. To evaluate the dependence of µ 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 µ = 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. (1.75, ∞)

(1.75, 3.5)

(1.80, 3.5)

(1.90, 3.5)

(2.00, 3.5)

(2.10, 3.

(2.20, 3.5)

(2.30, 3.5)

(2.40, 3.5)

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.

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 C , y C ) has been proven to be the highest density of shock-induced cavitation [START_REF] Obreschkow | Confined shocks inside isolated liquid volumes: A new path of erosion?[END_REF] . 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) 4 . 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 C , y C ) 2,13 . However, a close examination of the bubble cloud center as seen in the experimental observations shows discrepancies between the theoretical x 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 (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 C ratio over the pressure-relaxation rate µ. It is obvious that the theoretical value x C does not agree with the numerical value x c,n . Note that, in the absence of the bubble cloud growth (µ = 1), a very good match has nevertheless been reported between x C and the focal point of the reflected wave simulated. As µ decreases, the x c,n converges towards x C , while an asymptotic behavior of the relative error between x c,n and x C is reported when µ → ∞.

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 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 I 2 the intensity of the transmitted shock in this regime. However, for transmitted shock with intensity I 1 < I < 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 convergencedriven cavitation regime.

We denote, x t , the location of the expansion wave F 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 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 tx C . Assuming a volumetric collapse of the cloud, this implies that the cloud center ∆ x is given by

∆ x = x C + x t -x C 2 (25)
The bottom graph in Fig. 16 plots the x c,n /∆ x ratio over the pressure-relaxation rate µ, and shows a very good agreement between the bubble cloud center measured on the numerical simulation and the computed ∆ x . Note that the second term in the right-hand side of the Eq. 25, x tx C /2, is a corrective term to Obreschkow et al.'s theory to estimate the location of the cavitation region in the convergencedriven cavitation regime. Eq. 25 agrees to locate the highest density of shock-induced cavitation at the caustic's cuspidal point (lim

x t →x C ∆ x = x C ). In the immediate cavitation regime, ∆ x = x C + (R d -x 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

n 3n -1 ≤ ∆ x R d ≤ 1 2 1 + n 3n -1 (26)
Note that the analysis on the cloud center, here shown for the configuration of Sembian et al., is also valid for the highspeed droplet impact configuration.

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 µ = 3.5. A parametric investigation on µ 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 µ > 3.5, the size, the location and the collapse time of the bubble cloud is overestimated as µ increases, and conversely, these parameters are underestimated when µ, 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 amplifying expansion wave, so that ∆ x ≥ x C | k=2 . The correction of the ∆ x is done by adding the (x t -x C )/2 term to the F 2 's focal point location x C . Future work should address the analytical determination of the x t location, which is here numerically determined.

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 µ = 0.5 to best match the experimental observations. The droplet dependency on µ 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 analytical solution of the Keller-Miksis equation, the ]0.1, ∞] range has been found to very well agree with theory. Encouragingly, the two pressure-relaxation rates determined fits in this µ interval. Note that, in this work, the calibrated values of µ are only valid for the two configurations we have been concerned with and for an initial gas volume fraction α g = 10 -6 . Effects of α g on the phenomenology of the shock-induced cavitation require additional research efforts based on complementary experiments, which are out of the scope of this paper. However, to shed light on the possible influence of α g on the physics and the calibrated µ, 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 relaxation rate. A closer examination on the dependence of µ on the material properties should also be considered. In this work, we infer the difference in the µ 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 experimental 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 equations 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. of Sembian et al. and Field, Dear, and Ogren, α g , is an unknown parameter. We therefore estimated α 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 α 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 α g = 10 -5 and α g = 10 -7 . Figure 18 plots the variation of the volume of the gas phase within the droplet as time proceeds for various (α g , µ) combinations. Comparing the results for the calibrated µ = 3.50 with α 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 α g = 10 -5 and α g = 10 -7 with µ = 3.50. We then calibrated again the pressure-relaxation rates for α g = 10 -5 and α g = 10 -7 and found µ = 2.50 and µ = 4.60, respectively. It appears that changing α g by one order of magnitude (lower or higher) does not change the order of magnitude of the pressure-relaxation rate.

1 S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis, "Plane shock wave interaction with a cylindrical water column," Phys. Fluids 28, 056102 (2016).

)FIG. 2 :

 2 FIG. 2: Caustic traced out by the singular point (S ) of the singly-reflected wavefront. The caustic exhibits a cusp (C ).The red-to-blue colored segments are the singly-reflected wavefronts as time proceeds.

FIG. 3 :

 3 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.

FIG. 6 :FIG. 7 :

 67 FIG. 6: (a) Radial bubble-wall evolution for a spherical bubble collapse with p ∞ /p b = 1427 and N = 100 cells, with (b) a magnified view in the rebound region.

FIG. 8 :

 8 FIG. 8: Influence of µ 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.

FIG. 9 :

 9 FIG. 9: Comparison of the internal droplet structure, at t = 2.72, between (top) numerical simulations at different µ 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 µ = 1, is compared with the theoretical predictions given by Eqs. 1 for k ∈ [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.

  t e x i t s h a 1 _ b a s e 6 4 = " K W h 0 R L J 0 b w 8 e m / x 3 P U 2 + H I l N 2 F Q = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o M Q L 2 F X g n o M e v E Y 0 T w g W c L s p D c Z M j u 7 z M w K I e Q T v H h Q x K t f 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G t z O / 9 Y R K 8 1 g + m n G C f k Q H k o e c U W O l h 3 J w 3 i u W 3 I o 7 B 1 k l X k Z K k K H e K 3 5 1 + z F L I 5 S G C a p 1 x 3 M T 4 0 + o M p w J n B a 6 q c a E s h E d Y M d S S S P U / m R + 6 p S c W a V P w l j Z k o b M 1 d 8 T E x p p P Y 4 C 2 x l R M 9 T L 3 k z 8 z + u k J r z 2 J 1 w m

  FIG. 10: Dependence of xc and tc coordinates on µ. Both fits are non-linear least square fits of the form aµ 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.

FIG. 12 :

 12 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.

FIG. 13 :

 13 FIG. 13: Simulation results for µ = 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.

FIG. 15 :FIG. 16 :

 1516 FIG. 15: Comparison of the numerical simulation (left) with the experiments of Field, Dear, and Ogren 2 (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.

FIG. 17 :

 17 FIG. 17: (a-f) Comparison of the internal droplet structure between (top) numerical simulations at different µ 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..

FIG. 18 :

 18 FIG.18: Variation of the volume of the gas phase within the droplet as time proceeds for various initial gas volume fraction α g and pressure-relaxation rate µ.

TABLE I :

 I Weber and Reynolds numbers associated to the experiments of Sembian et al. and Field, Dear, and Ogren.

	Configuration	We	Re
	Shock-droplet 1 High-speed impact 2	∼ 10 3 ∼ 10 6	∼ 10 6 ∼ 10 6

TABLE II :

 II Relative error, ε r , between the numerics and the Keller-Miksis solution on the minimum bubble radius.

	µ	0.01	0.1	1	10	∞
	ε r (%)	735	275	10	19	18

b Note that Fig.11is plotted for µ = 3.5.

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

H. Chen and S. M. Liang, "Flow visualization of shock/water column interactions," ShockWaves 17, 309-321 (2008).

 4 L. Biasiori-Poulanges and H. El-Rabii, "Shock-induced cavitation and wavefront analysis inside a water droplet," Phys. Fluids 33, 097104 (2021).

ACKNOWLEDGMENTS

Authors thank Fabien Petitpas and Nicolas Favrie for their participation in the determination of the µ k, j formulation. Authors acknowledge the financial support from the ETH Zurich Postdoctoral Fellowship program.

DATA AVAILABILITY

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

CONFLICT OF INTEREST

The authors have no conflicts to disclose.

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 µ 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 µ = 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 µ = 1 and µ = 3.5 well simulate the diverging expansion wave, which is invisible in the simulations with µ = 10. Simulations with µ = 1 and µ = 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 µ = 10 is significantly overestimated. On frame (e), experiments show a bubble cloud (contoured with the red dashed line). In the simulation with µ = 1, no cloud is visible. In addition, the intensity of the numerical schlieren is the same for the various internal reflections, while the 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 µ = 3.5 are however in a very good agreement, which is also true for frame (f). Although the simulation with µ = 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 µ = 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 collpase of the gas phase involved in the liquid-gas mixture initially constituting the droplet. This requires to set the initial gas volume fraction, α g within the droplet. In the experiments