D Fuster
email: fuster@dalembert.upmc.fr

A review of models for bubble clusters in cavitating flows

Keywords: Cavitation, Bubble cluster models, Single bubble dynamics

This paper reviews the various modeling strategies adopted in the literature to capture the response of bubble clusters to pressure changes. The first part is focused on the strategies adopted to model and simulate the response of individual bubbles to external pressure variations discussing the relevance of the various mechanisms triggered by the appearance and later collapse of bubbles. In the second part we review available models proposed for large scale bubbly flows used in different contexts including hydrodynamic cavitation, sound propagation, ultrasonic devices and shockwave induced cavitation processes. Finally we discuss the main challenges of cavitation models.

Introduction

Cavitation and fast bubble dynamic processes have been investigated over the last century using theoretical, experimental and numerical techniques. The first investigations, back in 1917, aimed at understanding the response of hydrodynamic cavitating flows [166] and the interaction mechanisms between the bubbles and pressure waves [146,206,167]. Since then, the prediction and measurement of high temperatures during the collapse of bubbles, the emission of shock pressure waves and the appearance of high speed jetting have helped us to unveil various mechanisms taking place in cavitating flows, being a constant source of inspiration for engineers and researchers trying to take advantage of the processes triggered by the collapse of bubbles for diverse applications.

The main difficulty when investigating bubbly flows is the existence of fast out-of-equilibrium processes involving the collective response of bubbles. Despite the development of experimental techniques, the difficulty to perform measurements of fast processes taking place at small scales together with the intrinsic difficulties to investigate multiphase flows has motivated the development of models that can eventually be used to understand and guide the optimization of the process at hand. Controlling the process of bubble inception, predicting the nucleation frequency, understanding the interactions of turbulence with heat and mass transfer in multiphase flows, and understanding the influence of collective effects during the collapse of bubble clusters are still cornerstones for the optimization and development of a wide variety of applications based on cavitation.

The development of hydrodynamic cavitation models is an active area of research motivated by the major industrial challenges when determining the role of cavitation in various applications including diessel injection nozzles [START_REF] Giannadakis | Modelling of cavitation in diesel injector nozzles[END_REF][START_REF] Andriotis | Vortex flow and cavitation in diesel injector nozzles[END_REF], high speed flows around hydrofoils [105,170,106,118,102], and hydraulic machineries [137]. Erosion [START_REF] Franc | Material and velocity effects on cavitation erosion pitting[END_REF][START_REF] Dular | Development of a cavitation erosion model[END_REF] and how the turbulence is influenced by the presence of multiple phases in cavitating flows are examples of phenomena for which there is not yet a clear modeling strategy.

Another classical problem where the development of numerical methods and models is critical to guide the further development of technologies is the optimization of high intensity ultrasounds devices widely used to enhance mixing, to fragment solid particles and to promote certain chemical reactions. Single-Bubble SonoLuminiscence [START_REF] Gaitan | Sonoluminiscence and bubble dynamics for a single, stable, cavitation bubble[END_REF] is a classical example used to illustrate the extreme conditions generated during the collapse of bubbles in high intensity ultrasound fields but the application of ultrasounds for industrial applications remain a major challenge due to reproducibility and scalability problems.

A better understanding of the interactions among collapsing bubbles and a surrounding rigid or elastic wall will find application in various fields including surface cleanning techniques [START_REF] Chahine | Modeling of surface cleaning by cavitation bubble dynamics and collapse[END_REF][START_REF] Dijkink | Measurement of cavitation induced wall shear stress[END_REF], where the collapse of bubbles can be used to clean surfaces difficult to access, and various medical applications based on the interaction between pressure waves and bubbles. Shockwave lithotripsy, where the collapse of bubbles can play a crucial role in the process of kidney stone destruction [START_REF] Bailey | Progress in lithotripsy research[END_REF][START_REF] Freund | A cumulative shear mechanism for tissue damage initiation in shock-wave lithotripsy[END_REF]110,139], High Intensity Focused Ultrasound techniques (HIFU), where the shear stresses and thermal effects are two competing mechanisms influenced by the presence of bubbles with an unclear role in the overall efficiency of medical treatment techniques, and drug delivery applications, where encapsulated bubbles need to be transported and then exploded to release drugs in specific parts of the body significantly increasing the efficiency of the treatment process, are only some examples where experimental observations do not allow to discriminate between the various mechanisms proposed to explain the effects observed. Detailed numerical studies on the interactions between pressure waves, elastic solid surfaces and bubbles are currently being used for a better understanding of the process and to define precise strategies for the further development of these technologies as well as improving our understanding of processes such as the mechanisms of head injuries induced by shock impacts [START_REF] Earls | Cavitation in medicine[END_REF].

It is also important to mention applications related to wave propagation in multiphase flows where the precise role of very small amount of bubbles in a system is not fully clarified yet. In different geophysical problems it is important to identify the influence of the presence of gas in large scale systems such as marine ecosystems (e.g. oceans), where sonnars and seismometers are usually used to inferre the properties of the sediments and crustal layers, volcanos [187,218,168] and petrophysical systems where models are crucial to interpret seismic data to detect the presence of gas in geological strata [185,198,[START_REF] Grab | The effect of boiling on seismic properties of water-saturated fractured rock[END_REF]. The interesting reflection and transmission properties of diluted bubbly systems are currently used to design new metamaterials [129] and to protect offshore structures in processes such as drilling and controlled underwater explosions [START_REF] Grandjean | Shock propagation in liquids containing bubbly clusters: a continuum approach[END_REF].

This article is a review of available models for cavitating flows and bubble dynamics. After a qualitative description of the process of cavitation, we review models used to investigate the individual response of a spherical bubble to a pressure disturbance and the information they provide about the influence of various variables on the bubble response. Then we present a survey of currently available numerical methods and models able to solve for the three-dimensional collapse of a single bubble. In the second part of the manuscript we review multi-scale models for diluted bubbly flows and the extension of these models to problems where the presence of both phases co-exist without the possibility to consider one phase as disperse. In the last section we discuss the major challenges that cavitation models face and future lines of development.

Qualitative description of relevant physical phenomena in bubble clusters

The appearance and subsequent dynamic response of a bubble cloud have associated multiple problems with very different temporal and length scales that make the overall problem of cavitation and bubble cluster dynamics difficult to model. We include below a qualitative description of a general problem where cavitation occurs.

When the liquid pressure is lowered below a certain threshold cavitation inception occurs. Once the bubble nuclei are activated the extremely fast bubble expansion velocity induces a significant agitation in the surrounding liquid. The initial bubble expansion is typically a process much faster than any characteristic diffusion time-scale and therefore it is typically assumed that neither thermal effects nor mass transfer effects play a significant role in the initial bubble growth.

Then the pressure difference between the bubble and its surrounding drops quickly and the expansion velocity decreases, while the characteristic velocity of thermal and mass diffusivity effects (proportional to D t /R 0 ) increases as a consequence of the agitation induced during the initial stages of the inception process. Thus the Peclet and Sherwood number, representing the ratio between advection and heat/mass diffusion, decrease and mass transfer processes across the interface become relevant. At this stage bubbles under expansion start filling with vapor and the gas/vapor pressure recovers up to the vapor pressure.

If the process of bubble inception is intense enough, bubble coalescence and breakup processes and turbulence interact with each other developing large scale vapor pockets where surface tension forces are not sufficient to keep bubbles spherical. In this region, mass and heat transfer processes controlling the bubble growth are significantly influenced by the flow turbulence due to the enhancement on the effective thermal and mass diffusivity and the nucleation sites activated serve as a constant source of bubbles that grow and detach in a cyclic process [START_REF] Tf Groß | Diffusion-driven nucleation from surface nuclei in hydrodynamic cavitation[END_REF].

The process finishes when the pressure recovers and the gas/vapor mixture compresses making the vapor condense. In this stage the compression velocity can be much faster than the characteristic condensation velocity and bubbles collapse violently in an adiabatic process where neither heat nor mass transfer effects play a significant role. This process is typically followed by a series of rebounds that induce significant agitation on the flow at very small scales (or the order of few microns) until the vapor fully condenses and the remaining non-condensable gas, disgregated into small bubbles after the bubble collapse, reaches equilibrium with the surrounding liquid.

The mechanisms and the response of the individual bubbles as well as their collective response is markedly different in each of the stages described above requiring ad-hoc techniques difficult to extrapolate and generalize for all the stages of the cavitation process. While models devoted to capture the individual response of bubbles in a pressure varying field are useful to isolate the main mechanisms influencing the response of a gas/vapor bubbles, large scale models are devoted to capture the large scale interaction between phases and turbulent structures created in the flow at large scales. In the next sections, we review the main models proposed in the literature used to understand the response of single bubbles and bubble clusters.

Basic equations

Before getting to the specific modeling issues of cavitating flows it is useful to recall the basic equations for a multicomponent and multiphase system from the perspective of contiuum mechanics. Of course contiuum models have limitations that will be discussed further on, but still they constitute a solid base for the development of consistent models for bubbly flows.

The continuity and momentum equation for the ith phase are

∂ρ i ∂t + ∇ • (ρ i u i ) = 0, ( 1
)
∂ρ i u i ∂t + ∇ • (ρ i u i u i ) = -∇p i + ∇ • τ i , ( 2
)
where ρ is the density, u is the fluid velocity vector, p is the fluid pressure and τ = µ ∇u + (∇u) T -2 3 µ∇ • uδ is the viscous stress tensor depending on the dynamic viscosity µ. The interface position can be defined by tracking a function I defining the position of the interface and that must obey the advection equation,

∂I ∂t + u I • ∇I = 0, (3)
where u I denotes the interface velocity.

When we consider either compressibility or thermal effects we need to add the energy equation. One possibility is to use the total energy equation

∂ρ i E i ∂t + ∇ • (ρ i E i u i ) = -∇ • (u i p i ) + ∇ • (τ i u i ) -∇ • q i , (4)
where E i = e i + 1 2 |u| 2 is the total energy and q i = -κ i ∇T i is the diffusive heat flux obtained from the Fourier's law. Note that other forms of the energy equation are alternatives to close the system of equations. For instance the internal energy equation is

ρ i c p,i DT i Dt = β T,i T i Dp i Dt + Φ v -∇ • q i , (5)
where β T is the thermal dilatation coefficient and Φ v is the viscous dissipation function. This equation can be rewritten as an equation for pressure by using the relation of density, pressure and temperature differences through the definition of the speed of sound, polytropic coefficient and thermal dilatation coefficient

dρ = γ c 2 dp -ρβdT. (6)
Combining Eqs. 1, 5, 6 the following equation for pressure is found

1 ρ i c 2 eff,i Dp i Dt = -∇ • u i + β T,i Φ v,i ρ i c p,i - β T,i ρ i c p,i ∇ • q i , (7)
where

1 ρc 2 eff,i = γi ρic 2 i - β 2
T ,i Ti ρicp,i . Note that for pure substances such as ideal gases (β T ≈ 1/T ) and liquids like water (γ ≈ 1, β T ≈ 0) the classical equation for pressure

1 ρ i c 2 i Dp i Dt = -∇ • u i (8)
is justified if we neglect viscous dissipation and thermal diffusion effects.

In those cases where we need to model the influence of various components (eg. air and water) we add the species conservation equation for the mth component

ρ m DY i,m Dt = ∇ • D M i/m ∇Y i,m (9)
where Y i,m is the mass fraction of the component in the ith phase and

D M i/m
is the mass diffusion coefficient of the component in the mixture.

Across the interface the local balances establish the relation between the liquid and gas/vapor properties at the interface

J = ρ g (u I -u g ) • n = ρ l (u I -u l ) • n (10) J u • n = p + τ • nI + σκ (11) κ∇T • n = J∆H v (12) J Y m + ρ l D M m/i ∇Y m • n = 0, ( 13
)
where • denotes the jump of a given variable through the interface region, σ is the surface tension coefficient, n is the normal vector to the interface and J is the mass transfer flux. This flux can be obtained by assuming that equilibrium conditions prevail at the interface at every instant but it is also possible to relax this assumption and to use a model providing the mass transfer flux provided the instantaneous local conditions at the interface. For example in the 101] used in the context of simulation of bubbles in [149,211,[START_REF] Hauke | Dynamics of a single cavitating and reacting bubble[END_REF], the vaporization flux is obtained as

J = α evap (p I eq -p I b,vap ) 2πr vap T int , (14)
where p I eq is the vapor equilibrium pressure computed with the interface temperature and p I b,vap is the actual vapor pressure.

The system of equations is closed by chosing an Equation Of State (EOS) for each phase. From the state equation we can obtain various quantities such as the speed of sound, the thermal dilatation coefficient, etc. It is very popular in the context of cavitating flows to use the general stiffened EOS written in the Mie-Gruneisen form (see for example [173,108])

ρ i e i = p i + Γ i Π i Γ i -1 (15)
which defines the sound speed as

c 2 i = Γ i p i + Π i ρ i . ( 16
)
but any other choice of the EOS is also possible. Understanding and modeling the response of single bubbles is crucial before building more complex models involving the interaction of various bubbles.

Rayleigh [166] and later contributions by Plesset [157] developed the simplest of the models capturing the dynamic response of a single bubble in a pressure varying field. Imposing spherical symmetry, neglecting heat and mass trasfer fluxes across the interface and considering the liquid an incompressible substance, the continuity and momentum conservation equations in the liquid (Eqs. 1-2) can be integrated from the liquid interface to infinity to obtain the Rayleigh-Plesset (RP) equation

RR + 3 2 Ṙ2 = p I l -p ∞ ρ l ( 17
)
where p ∞ is the pressure far from the bubble (assumed to be a known function of time) and p I l is connected to the bubble pressure through the Laplace equation (Eq. 11)

p I l + 2σ R + 4µ l Ṙ R = p b . (18)
When the gas content obeys the perfect gas equation, the bubble pressure is assumed to be uniform and the response of the bubble is adiabatic the bubble pressure and bubble radius are directly related through the equation,

p b = p 0 R 0 R 3γ . ( 19
)
This simple model is able to capture the radial motion of bubbles and it is still widely used nowadays to predict the temporal evolution of the radius of a gas bubble in a time-varying pressure field.

A simplified solution of equation 17 during the bubble expansion is found in the controlled inertia regime where we neglect the acceleration term and one can directly relate the interface radial velocity with the pressure difference as

Ṙ = 2 3 p l (R) -p ∞ ρ l . ( 20
)
This expression only provides a real solution when p l (R) ≥ p ∞ and it is usually used to predict the expansion rate of bubbles after inception.

Another simplified analytical solution is obtained for the so-called Rayleigh collapse problem, where a bubble with some initial pressure p 0 is suddenly exposed to a liquid at a reference pressure p ∞ . In absence of surface tension the Rayleigh collapse time is [START_REF] Brennen | Cavitation and Bubble Dynamics[END_REF] t RC = 0.915R 0 ρ l p ∞ -p 0 [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] which defines the averaged compression velocity

Ṙ = 1 0.915 p ∞ -p 0 ρ l . ( 22
)
Equations 20 and 22 indicate that bubble expansion and collapse are almost symmetric processes where the averaged bubble expansion or compression velocity is proportional to the square root of the pressure difference between the bubble and the far field pressure. As we will see later, this has been a key ingredient on the development of effective phase-change models for largescale cavitating flows.

The Rayleigh-collapse time gives a very good estimation of the collapse time experimentally observed in a variety of situations, but the instantaneous velocities that the bubble interface reaches are significantly larger than those predicted by Eq. 22. Whereas the collapse time only depends on the maximum bubble radius and the external pressure when p ∞ p 0 , the peak value of pressure determining the intensity of the collapse depends on the pressure ratio p ∞ /p 0 , which can be very large for low values of p 0 .

During the instants close to the time at which the bubble radius is minimum various effects neglected in the initial RP model play a crucial role in the conditions generated during the bubble collapse. A strong assumption during the last stages of the bubble collapse is the liquid incompressibility condition. Keller and Miksis [115], Gilmore [START_REF] Gilmore | The growth or collapse of a spherical bubble in a viscous compressible liquid[END_REF] and Tomita [193] have proposed different corrected versions of the Rayleigh-Plesset equation that accounts for compressibility effects showing that the amplitude of the rebound is significantly attenuated during the bubble collapse processes. The relevance of compressibility effects and the accuracy of the various modified equations is discussed by Prosperetti and Lezzi [165,130], and an interesting assymptotic analysis on the radiative decay of non-linear bubble oscillations at small Mach numbers has been recently presented by Smith & Wang [186]. For strong bubble collapses mechanisms such as the damping originated by emission of shock waves during the collapse are not captured by the corrections mentioned above and one has to resort to the numerical solution the full Navier-Stokes in the liquid [START_REF] Fuster | Liquid compressibility effects during the collapse of a single cavitating bubble[END_REF]122,[START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF]. These studies reveal that the errors introduced by Rayleigh-Plesset-like models is not negligible for very strong collapses, where the emission of a shock wave from the bubble seems to play a major damping role that is manifested in the differences on the amplitude of the bubble rebound amplitude predicted by RP models and Direct Numerical Simulations (DNS). Note that for very strong collapses other assumptions made in the initial model are questionable. For instance the perfect gas model is probably not sufficient to represent the bubble behavior when plasma is generated during very short instants at the bubble core and one needs to accurately represent the gas response and its properties for an accurate estimation of the conditions generated inside the bubble.

Gas/vapor bubbles

An important limitation of the original Rayleigh-Plesset model is that it neglects mass and heat transfer processes across the interface. Several theoretical and experimental studies show that heat-and mass-transfer effects have a non-negligible influence on the bubble response [START_REF] Brennen | Cavitation and Bubble Dynamics[END_REF][START_REF] Brennen | The dynamic behavior and compliance of a stream of cavitating bubbles[END_REF]128,204,[START_REF] Dular | Thermodynamic effects during growth and collapse of a single cavitation bubble[END_REF]. Models including heat and mass transfer processes are typically based on the modified version of the Rayleigh-Plesset equation

[163] R R = - 3 2 Ṙ - J ρ l 2 + R J ρ l + 2 J ρ l Ṙ - J ρ l + p b -p ∞ ρ l - 2σ ρ l R - 4µ l ρ l Ṙ R (23)
For a general case the hypothesis of an adiabatic transformation does not hold and bubble pressure and volume changes are influenced by heat and mass transfer processes across the interface. Under the assumption of uniform pressure inside the bubble it is possible to combine the continuity, momentum and the internal energy equations inside the bubble to obtain an ODE for the bubble pressure [START_REF] Commander | Linear pressure waves in bubbly liquids: Comparison between theory and experiments[END_REF]162]

ṗb = - 3 R (γ -1)q I b + γp b Ṙ -γ p b ρ I b J , (24)
where

q I b = -κ b ∂T b
∂r r=R is the heat flux across the interface from the bubble. It is easy to show that when the heat and mass transfer flux across the interface are neglected the bubble response is adiabatic and equations 19 and 24 are analogous. In the limiting case of pure vapor bubbles at low frequencies in which the bubble pressure remains constant and equal to the vapor pressure at every instant, ṗ ≈ 0, we can use Eq. 24 to obtain an explicit equation for the bubble interface velocity as a function of the heat and mass fluxes across the interface

Ṙ = - γ -1 γ q I b p b + J ρ I b . (25)
Because heat and mass fluxes are indeed related through Eq. 12 we can write the bubble expansion velocity as a function of the heat fluxes across the interface only as

Ṙ = γ -1 γ q I b p b + q I b l ∆H v ρ I b ( 26
)
which can be further simplified if we neglect temperature gradients inside the bubble to

Ṙ = κ l ∆H v ρ I b ∂T l ∂r r=R . (27)
A particular solution of this equation is obtained by Plesset [159] for vapor bubbles growing in supersaturated liquids in absence of convection

Ṙ = 3 π ρ l c p ∆T ρ v ∆H v D T l t . (28)
This equation reveals that, unlike the case of inertial bubble growth derived in Eq. 20, if the bubble pressure reaches a regime where it is able to sustain an approximately constant value close to the vapor pressure the dynamic response of the bubble is not controlled by pressure differences. Instead heat and mass transfer effects across the interface become important and parameters such as the thermal boundary layer thickness in the liquid around the bubble play a major role in the dynamic response of the bubble.

The transition between regimes where the bubble response is controlled either by pressure differences or by liquid heat diffusion processes is not clear in real cavitating flows and we may expect to find different regimes depending on the specific application. For instance, during the initial process of inception and bubble collapse the inertia controlled regime applies (Eq. 20) given that the characteristic response time of the bubble is much faster than the characterstic times of mass and heat transfer. But after inception, if the system is able to evolve to regimes where the gas/vapor pressure is approximately equal to the vapor pressure, the process of bubble growth from walls and the growth of large scale bubble clouds may be controlled by liquid heat diffusion effects (Eq. 27) where certainly turbulence will play a role in the dynamic growth of the bubbles. This is likely the case in hydrodynamic cavitation flows.

The simplest of the models able to deal with the presence of gas/vapor mixtures is the full equilibrium model, also known as cold-vapor model, which assumes that the bubble reaches equilibrium conditions with its surrounding at every instant and that the temperature remains constant. In this case, the influence of the vapor pressure on the bubble pressure is simply modeled as a function of the partial pressure of non-condensable gas, p g , and the vapor pressure at the reference temperature, p v (T 0 ),

p b (R, T ) = p g (R) + p v (T 0 ), (29)
neglecting all flux terms in the Rayleigh-Plesset equation. The full equilibrium model has strong limitations, for instance, it is unable to recover the limiting solution obtained by Plesset [159] given that it neglects the influence of any transient effect related to the mass and heat transfer process on the dynamic response of the bubbles.

Preston et al [162] have used linear solutions to analytically obtain the overall heat and mass transfer coefficient across the interface for vapor diffusion limited problems in order to estimate the heat and mass transfer fluxes needed in Eq. 24. The main advantage of this approach is that it explicitly evaluates the heat and mass transfer flux across the interface and therefore the evolution of the bubble pressure is simply given by the solution of the ODE given in Eq. 24 which is easily coupled with the solution of the Rayleigh-Plesset equation. However this model still assumes that the interface temperature is constant and that the mass flux is controlled by vapor diffusion inside the bubble, which is problematic for pure vapor bubbles where heat transfer effects in the thermal boundary layer in the liquid must be modelled.

More sofisticated models accounting for heat and mass transfer fluxes across the interface account for the radial variation of properties inside the bubble and must deal with the numerical difficulties related to the discontinuity on the velocity, pressure and the heat fluxes across the interface imposed by the local balances across an infinitely thin interface. Models initially proposed by Nigmatulin and Khabeev [149], and later developed by different authors [START_REF] Gumerov | Weakly non-linear oscillations of the radius of a vapour bubble in an acoustic field[END_REF][START_REF] Hao | The dynamics of vapor bubbles in acoustic pressure fields[END_REF]188,131,208,[START_REF] An | Model of single bubble sonoluminiscence[END_REF]211,[START_REF] Hauke | Dynamics of a single cavitating and reacting bubble[END_REF]117], include the solution of the energy equation both inside and outside the bubble and the solution of the species equation to capture the process of vapor diffusion within the gas/vapor mixture inside the bubble. Most of these works also model the kinetics of phase change through an equation that allows for an explicit calculation of the flux as a function of the difference between the equilibrium pressure at the interface conditions and the current vapor pressure (Eq. 14). This expression is useful for numerical purposes as it allows for an explicit evaluation of the flux at a given instant but it introduces some uncertainty on the values of the accomodation coefficient α evap required in the model [START_REF] Gumerov | Determination of the Accomodation Coefficient Using Vapor/Gas Bubble Dynamics in an Acoustic Field[END_REF].

The overall mass transfer flux can be seen as a combination of mechanisms in serial, each of them introducing a certain resistance on the overall mass transfer process, and therefore the total flux is given by the slowest mechanisms which usually are diffusion mechanisms such as vapor mass diffusion inside the bubble or thermal diffusion in the liquid. In these situations exact value of the accommodation coefficient is not critical to obtain physically meaningful results. In particular, in the linear oscillation regime it is possible to theoretically show that the influence of α evap is restricted to situations where the vapor content is large and the excitation frequency is very high. Otherwise heat and mass diffusion mechanisms control the overall mass flux across the interface (e.g. when the Sherwood number Sh = ωR 2 0 Dm is of order one or below [START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF]). In the case of a strong bubble collapse, Storey and Szeri [188] have shown that the presence of small amounts of gas limit the condensation rate of vapor inside the bubble by mass diffusion. This mechanism prevents the bubble to disappear and limits also the peak compression velocities that one can reach. Even when vapor mass diffusion controls the overall collapse processes numerical studies have shown that the kinetics of phase change influence the dynamics of vapor bubbles [START_REF] Gumerov | Dynamics of vapor bubbles with nonequilibrium phase transitions in isotropic acoustic fields[END_REF] and the violence of the collapse generated [START_REF] Fuster | Influence of accommodation coefficient on nonlinear bubble oscillations[END_REF]. Parametric studies based on a complete one-dimensional model accounting for mass transfer effects reveal that the dependency of the peak temperatures reached during the bubble collapse on the ambient temperature is directly linked to the mass transfer model [START_REF] Fuster | Parametric analysis for a single collapsing bubble[END_REF], revealing the intricate mechanisms between the process of phase change and the collapse process of a gas/vapor bubble.

Gas/vapor reacting bubbles

Reaction processes can become relevant during the bubble implosion when the extreme conditions reached inside the bubble trigger various chemical processes associated to the generation of radicals from gas molecules. Even from pure water vapor molecules radicals can be generated at the bubble core during the bubble implosion in an endothermic process that significantly limits the peak temperatures predicted [START_REF] Didenko | Effect of noble gases on sonoluminescence temperatures during multibubble cavitation[END_REF]. Note that these processes are significantly attenuated at the bubble interface where the interface temperature remains low and therefore radicals recombination occurs.

Particularly interesting is the investigation of bubbles containing reacting mixtures, which finds application in the promotion of chemical process by the use of cavitation [135]. Under particular conditions the mixture gas detonation leads to the appearance of shockwaves inside the bubble that break the assumption of pressure uniformity inside the bubble. Figure 1 depics the radial pressure and reaction heat profiles during the collapse of a single bubble containing an initial mixture of H 2 -O 2 -H 2 O -Ar in proportions near stoichiometric conditions where the reaction reaction heat (q r ) is represented using the derived variable, q r , q r = sign(q r ) log

10 |q r | |q r | > 1 0 |q r | ≤ 1 (30)
defined positive when the reaction is exothermic and negative when the reaction is endothermic. The simulations are obtained using the model proposed by Hauke et al [START_REF] Hauke | Dynamics of a single cavitating and reacting bubble[END_REF] and the chemical reaction mechanism proposed in [START_REF] Akbar | Detonations in H2-N2O-CH4-NH3-O2-N2 Mixtures[END_REF].

The simulations show how during the compression stage molecules contained inside the bubble are broken into radicals in an endothermic process (blue zone for t < 2.602 • 10 -5 s) when the temperature and pressure are above a given threshold. This process mainly take place at the bubble core because at the interface, the temperature remains low and no radicals are generated. In fact, the radicals which are produced inside the bubble are recombined nearby the interface when they diffuse towards it, being the reason of the exothermic zone observed near the interface in Figure 1.

As in this case, if we reach autoignition conditions at the bubble center a flame propagates towards the interface at a velocity of around 1500 m/s. As the flame approaches the interface the conditions are less extreme and the flame speed decreases until it is finally extinguished. The combustion takes place very fast and the reactives are quickly exhausted at the center whereas at the interface, where the reaction rate is negligible, the concentration remains high. The water molecules generated during the combustion process are broken in an endothermic process which again mainly takes place at the Fig. 1: Spatial and temporal representation of the reaction heat (q r ) and pressure during the implosion. The uniform pressure assumption does not hold and shockwaves induced by the detonation process are observed. bubble center. Interestingly the molecule dissociation process makes the total number of molecules at the center sudently increase creating a pressure difference between the bubble center and the interface that can lead to internal shockwaves inside the bubble that break the classical pressure uniformity assumption (see Figure 1).

Another important conclusion from these models is that the generation of radicals mainly takes place at the bubble center, being the concentration near the interface significantly reduced due to the low temperature of the interface influencing the efficiency of cavitation induced reactions for nonvolatile molecules.

Theoretical solutions for linear bubble oscillations

The linear solution of the full model equations is one of the few theoretical solutions available to better understand the response of bubbles to a pressure perturbation for a given known frequency ω. In addition they can also be used to validate numerical methods [START_REF] Bergamasco | Oscillation regimes of gas/vapor bubbles[END_REF]. Assuming solutions of the form R = R 0 (1+R exp (ıωt)) for a given external excitation p ∞ = p 0 (1+p ∞ exp (ıωt)), with both the radial and pressure perturbation being much smaller than unity, the Rayleigh-Plesset equation simplifies to

R = - 1 ω 2 0 -ω 2 + 2iδω p ∞ ρ l R 2 0 p ∞ , (31)
which establishes the relation between external pressure disturbances and bubble volume changes.

Analytical expressions for the damping factor, δ, for pure-gas bubbles have been reviewed in [START_REF] Ainslie | Review of scattering and extinction crosssections, damping factors, and resonance frequencies of a spherical gas bubble[END_REF]. The resonance frequency or Minnaert frequency [146], ω 0 , can be interpreted as the frequency at which the bubble oscillation amplitude for a given pressure perturbation is maximal.

The approximated expression of the resonance frequency,

ω 2 0 ≈ p b,0 ρ l R 2 0 (Φ), (32)
also determines the value of the amplitude of the bubble oscillation in the low frequency limit (ω ω 0 ), when Equation 31 can be rewritten as

R (ω → 0) = - p ∞ ρ l ω 2 0 R 2 0 p ∞ ≈ - p ∞ (Φ) . ( 33
)
This relation shows that volume and pressure disturbances are linked through the transfer function Φ = -p b /R , which represents the relation between the bubble volume (e.g. density) and bubble pressure fluctuations in the linear limit. It is possible to interpret this function as an effective polytropic coefficient for the bubble p b /ρ γ eff b = C, where γ eff = Φ/3. Should Φ be a real function, a barotropic Equation Of State (EOS) can be defined for the gas/vapor mixture as pressure can be directly computed from density. Otherwise both density and temperature are needed to compute the pressure inside the bubble.

Generalized solutions for gas/vapor bubbles have been discussed in [116,149,[START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF]164]. The general expression for the transfer function Φ found in [START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF] shows that although this function depends on several nondimensional quantities, the problem is reduced to the nondimensional bubble oscillation frequency (e.g. the Peclet number) and the vapor fraction Y 0 when all the physical properties of the system are fixed. Figure 2 shows the general solution of Φ as a function of the bubble Peclet number, Pe = ωR0 D T b /R0 , and the vapor content for an air/water system. The limiting solutions of the effective polytropic coefficient include the isothermal limit γ eff → 1 (for gas bubbles at low frequencies), the adiabatic limit γ eff → γ (for gas bubbles at high frequencies) and the constant vapor pressure limit for pure vapor bubbles oscillating at low frequencies (γ eff → 0) where the bubble pressure is constant if the temperature is constant, p b = p vap (T 0 ).

In order to distinguish between various oscillation regimes we can compare the full expression for Φ reported in [START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF] with that obtained with simplified models typically used to represent the bubble response such as the adiabatic gas model (Φ ad = 3γ), the full equilibrium model (Φ eq = 3(1 -Y 0 )) and models accounting for mass transfer effects limited by either vapor diffusion or thermal diffusion in the liquid boundary layer [START_REF] Bergamasco | Oscillation regimes of gas/vapor bubbles[END_REF]. Figure 3 shows the regime map of the mechanisms controlling the bubble response as a function of the vapor content and the Peclet number. The diagram is built by represeting those regions where different approximate expressions of the transfer function contain less than 10% error compared to the exact value of Φ. We [START_REF] Bergamasco | Oscillation regimes of gas/vapor bubbles[END_REF]. distinguish a high frequency regime (Pe b > 100) where neither mass nor heat have time to respond to pressure changes and both, gas and vapor respond almost adiabatically. In this regime Eq.19 is valid taking the total pressure inside the bubble, including both gas and vapor, as p 0 . This situation is likely to occur during the initial instants of the bubble expansion after inception and also during the violent collapse of bubbles when the vapor behaves almost as a non-condensable gas. At low frequencies, different regimes can be defined depending on the oscillating frequency and the vapor content. The bubble response is isothermal and the full equilibrium model is applicable for a range of vapor content that increases as we decrease the excitation frequency. Preston's model [162], that includes the full-equilibrium model as a particular solution, applies below the red line and captures the transition between the full-equilibrium model and the adiabatic model for bubbles containing a non-negligible amount of gas. As the vapor fraction approaches to one liquid heat diffusion controls the mass transfer rate across the interface and therefore the dynamic bubble response.

It is particularly interesting to investigate the transition between the full equilibrium model and the region in which heat diffusion effects in the surrounding liquid control the bubble response. Theoretically this transition depends on the critical Peclet number [START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF] Pe

c = 3 1 -Y 0 Y 0 γ -1 γ c p,b T 0 ∆H v 2 κ l κ b , (34)
where heat diffusion does control the bubble response when Pe b Pe c . This condition can be also expressed as

1 3 Y 0 1 -Y 0 γ γ -1 ρ l ρ b Σ * 1 ( 35
)
where the parameter [START_REF] Brennen | The dynamic behavior and compliance of a stream of cavitating bubbles[END_REF] to evaluate the relevance of thermal effects. This result reveals that liquid heat diffusion effects are specially relevant in systems with large enthalpy of vaporization with small liquid's diffusivities, being thermal effects always important in systems close to the saturation line (Y 0 → 1). This is consistent with the Plesset's solution obtained for pure vapor bubbles [159] where the bubble growth is given by transient heat transfer effects in the thermal boundary layer surrounding the bubble (Eq. 28). Otherwise, full-equilibrium conditions prevail at the interface and the bubble content reaches instantaneous equilibrium with its environment at every instant.

Σ * = ΣR 2 0 ω c p,b T0 √ D T L is a non-dimensional version of the Brennen's parameter Σ = ρ 2 b ∆H 2 v ρ 2 l c p,l T0 √ D T l
Also note that the imaginary part of the transfer function becomes of the order of the real part in those cases where the vapor content is large and liquid heat diffusion control the bubble response (Figure 2), therefore, the transition between the full-equilibrium model and the heat and mass diffusion regimes provide a theoretical background about the range of applicability of the barotropic approximation in the linear regime. Certainly nonlinear effects have an influence on the regime maps theoretically obtained and although preliminary numerical simulations indicate that the overall qualitative picture does not change [START_REF] Bergamasco | Oscillation regimes of gas/vapor bubbles[END_REF], more detailed analyses about regimes where either gas/vapor expansion or mass transfer processes control the dynamic response of the bubbles are required.

Direct Numerical Simulation of cavitating bubbles

While useful to discuss the main mechanisms controlling the response of bubbles to pressure changes, one-dimensional models have serious limitations. It is well known that the presence of any source of asymmetry (e.g. wall or gravity) leads to the appearance of a high speed jet during the collapse of the bubble [150,189]. The generation of a high velocity jet [127], the large viscous stresses in bubbles collapsing close to a solid wall [214] and the appearance of vortex rings responsible of a significant energy dissipation and agitation at small scales, are only few examples of mechanisms that cannot be captured by the solution of a spherically symmetric problem. Three dimensional effects are also relevant at large scales when surface tension effects are negligible and the interface is no longer spherical. Numerical methods are then a very useful and powerful tool to investigate the impact of three-dimensional effects on the dynamic response of bubbles and its consequences on the surrounding liquid and nearby walls. In this section we review the main methods developed to investigate the three-dimensional response of single bubbles.

Models

Initial models proposed for the simulation of bubbles in liquids [158,[START_REF] Blake | Growth and collapse of a vapour cavity near a free surface[END_REF][START_REF] Blake | Cavitation bubbles near boundaries[END_REF] directly solved for the advection of the interface position (Eq. 3) using the Front Tracking (or marker) method that consists in updating a set of discrete points sitting at the interface using the kinematic condition dx I dt = u I . An important limitation of this technique is that it is not trivial to deal with coallescence and breakup processes. The advantage is that like in Rayleigh-Plesset-based models, under the assumption of uniform properties inside the bubble, the equations solved are significantly simpler. As discussed above, modeling the bubble interior as a region of uniform pressure directly relating the bubble pressure and volume is a reasonable approximation in various situations where the bubble response is either fully isothermal or adiabatic. By tracking the interface position explicitely the discontinuities imposed by the local balances across the interface (Eqs. [START_REF] Andriotis | Vortex flow and cavitation in diesel injector nozzles[END_REF][START_REF] Apfel | The role of impurities in cavitation-threshold determination[END_REF] are naturally imposed. The most important methods resorting to these assumptions are Boundary Integral (BI) methods and free surface methods.

In BI methods (also known as Boundary Element Method) the liquid flow is assumed to be potential which allows directly solving for the velocity of a Lagrangian set of points following an interface without the need of discretizing the full Navier-Stokes in the bulk. The fact that the flow is potential limits the range of application of these models to Stokes and inviscid (high Reynolds) flows. Also, because the properties inside the bubble are assumed to be uniform these models cannot capture all the phenomena taking place inside the gas phase (vapor/gas diffusion, heat transfer, chemical reactions, turbulence, etc...). Firstly used to investigate the dynamic response of bubbles by Plesset & Chapman [158] and developed later by Blake and coworkers to simulate problems involving the three dimensional collapse of a single bubble [START_REF] Blake | Growth and collapse of a vapour cavity near a free surface[END_REF][START_REF] Blake | Cavitation bubbles near boundaries[END_REF], these models have been useful to investigate the bubble dynamic response near walls and the process of jet formation and its consequences at large Reynolds in incompressible liquids [START_REF] Blake | Growth and collapse of a vapour cavity near a free surface[END_REF][START_REF] Blake | Cavitation bubbles near boundaries[END_REF][START_REF] Chahine | The final stage of the collapse of a cavitation bubble near a rigid wall[END_REF]209]. Recent extensions have included liquid compressibility effects [119,201,203,202] but the fundamental limitations of these models remain.

Free surface methods generalize the models used in BI methods by taking into account viscous effects in the liquid bulk. Initially proposed by Harlow and Welch [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF] and used in the context of bubble collapse in the numerical methods proposed by [213,160,[START_REF] Delale | Direct numerical simulations of shock propagation in bubbly liquids[END_REF], these methods are based on the discretization and solution of the full Navier-Stokes equations in the liquid phase where the presence of the gas phase imposes the boundary condition at the interface through the interface balances (Eqs. [START_REF] Andriotis | Vortex flow and cavitation in diesel injector nozzles[END_REF][START_REF] Apfel | The role of impurities in cavitation-threshold determination[END_REF][START_REF] Atchley | The crevice model of bubble nucleation[END_REF][START_REF] Azouzi | A coherent picture of water at extreme negative pressure[END_REF]. Simulations reported in [160] have revealed that viscous effects become important during the process of jetting formation influencing significantly the peak velocities that can be reached. Delale et al [START_REF] Delale | Direct numerical simulations of shock propagation in bubbly liquids[END_REF][START_REF] Delale | Shock structure in bubbly liquids: comparison of direct numerical simulations and model equations[END_REF] have used these methods to investigate the structure of shock wave propagation in incompressible liquids where all the effective medium compressibility can be attributed to the presence of compressible bubbles.

In recent years, models based on Eulerian-Eulerian approaches have attracted most of the attention because they deal with the process of fragmentation/coallescence naturally without the need of highly deformed meshes that adapt to the interface contour at every instant. Two key ingredients for these methods are: (i) how to represent and advect the interface position and (ii) how to take into account the discontinuity on the field variables established by the local balances (Eqs. [START_REF] Andriotis | Vortex flow and cavitation in diesel injector nozzles[END_REF][START_REF] Apfel | The role of impurities in cavitation-threshold determination[END_REF][START_REF] Atchley | The crevice model of bubble nucleation[END_REF][START_REF] Azouzi | A coherent picture of water at extreme negative pressure[END_REF].

The interface is typically represented using the volume fraction of the reference phase α that obeys the advection equation

∂α ∂t + u I • ∇α = 0. ( 36
)
Because in the general problem where phase change occurs the velocity is not continuous through the interface it is useful to rewrite the equation above as a function of any arbitrary velocity field u as

∂α ∂t + u • ∇α = (u I -u) • nδ(x I ) ( 37
)
where we have used the relation ∇α = -nδ(x I ) to obtain a term on the right hand side of the equation that, using Eq. 10, can be related to the change on the volume fraction of the gas phase due to vaporization in those places where an interface exists. Note that most of current models presented for the simulation of single bubbles neglect mass transfer effects across the interface (large Peclet regime) and therefore we can write Eq. 37 as

∂α ∂t + ∇ • (uα) = α∇ • u. (38)
In addition to the representation of the interface position, the main differences between Eulerian-Eulerian models stem from the way the jump on properties across the interface are handled. Some models relax the constrains imposed by balances across the interface still solving for the distribution of variables inside the bubble. In the seven equation model initially proposed in [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (ddt) in reactive granular materials[END_REF], pressure temperature and velocities can be different for each phase contained in a given computational grid and additional terms are added in the equations to model the influence of the complementary phase on the evolution of momentum and energy on a given phase. Different variants of this model have been proposed based on additional constrains on the velocity, pressure and temperature differences in mixed cells [114,174,176,154] being these models an active area of research in various problems involving two phase flows, phase change, etc [START_REF] De | Homogeneous two-phase flow models and accurate steam-water table look-up method for fast transient simulations[END_REF][START_REF] Lorenzo | Hllc-type and path-conservative schemes for a single-velocity six-equation two-phase flow model: A comparative study[END_REF]. Their use for DNS of cavitating bubbles is limited given that the system of equations solved is modified to introduce relaxation parameters that only approximate the dynamics of immiscible two-phase flows. How to calibrate the relaxation time required by the model based on physical principles, the convergence properties of three dimensional methods in problems involving interfaces and how accurate these approaches represent the solution of the exact problem where the interface is infinitely thin are still open questions [175].

Another popular alternative is the homogeneous model where we write evolution equations for all the averaged quantities using the volume fraction α by multiplying each conservative equation for the ith phase by its corresponding volume fraction, summing over all the phases to obtain evolution equations for each averaged quantity. Compared to the corresponding equations for single phase flows, the conservation equations for the mixture include additional terms that account for the influence of jump conditions on the averaged fluid quantities (see for example [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] for the treatment of surface tension effects). In particular the averaged Navier-Stokes equations in presence of surface tension read

∂ρ ∂t + ∇ • (ρu) = 0, ( 39
)
∂ρu ∂t + ∇ • (ρuu) = -∇p + ∇ • τ + σκ∇α, (40)
∂ρE ∂t + ∇ • (ρEu) = -∇ • (up) + ∇ • τ u -∇ • q + u • σκ∇α, (41)
where formally we have considered that in those cells containing an interface the pressure jump due to viscous and mass transfer effects is small and therefore p 1 -p 2 ≈ σκ.

The EOS for the mixture required to close the system of equations is different depending on the model assumptions. The classical approximation for the speed of sound of a mixture

1 ρc 2 = α g ρ g c 2 g + α l ρ l c 2 l ( 42
)
is only formally valid for adiabatic systems where there is not either heat nor mass exchange across the interface (adiabatic conditions). Some other relations can be applied by assuming both phases in thermal equilibrium or in full thermodynamic equilibrium but the range of validity of these hypothesis for a general problem where dynamic effects become important remain unknown. How to define EOS for mixed cells is a major challenge for the correct simulation of cavitating flows including mass transfer effects and the interactions between interface and shockwaves.

Numerical methods

In addition to the physical equations, numerical methods have an important impact on the quality of the solutions obtained and the phenomena that can be captured. One important ingredient of numerical methods dealing with multiple phases is how to track and advect the interface position. Front Tracking methods, used in the Boundary Integral methods and free surface solvers mentioned above, are rarely developed nowadays in the context of bubble dynamics (an exception is [100]) and most of current lines of research on the numerical simulation of cavitating bubbles focus on Eulerian-Eulerian approaches.

The two most important numerical methods used to solve Eq. 38 are the Level Set (LS) and the Volume Of Fluid (VOF) method. LS methods are based on a continuous and smooth function that takes a given value at the interface location being Diffuse Interface Methods (DIM) a particular family of them [START_REF] Abgrall | Computations of compressible multifluids[END_REF]173]. These methods resort to the use of numerical schemes that introduce a certain amount of artificial diffusion relaxing the representation of the interface as a discontinuous transition between two phases. Thus, both the numerical method and the model must deal with interfaces that have a certain diffusive length. Precisely one of the main drawbacks of these algorithms is the natural trend of the method to artifically diffuse interfaces degrading the quality of the solution near the interface with respect to the description of the model problem where the interface is assumed to be sharp. Although some remedies have been proposed [151,181,[START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF], this limitation still poses problems, especially at long times, in those regions where the interface contour is not correctly resolved.

Riemann methods coupled with DIM have been sucessfully applied to investigate the bubble collapse in a compressible liquid using quasi-conservative problem proposed by Abgrall for the simulation of compressible multicomponent problems [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF]. Johnsen and Colonius [109] and Tiwari et al [190] have extended Abgrall's formulation to solve the Euler equations using shockcapturing WENO schemes, showing the capability of these kind of methods to simulate the collapse of an air bubble in water in absence of surface tension forces [110,111] as well as the collapse of bubble clusters near a wall [171,191]. More recent formulations introduce viscous effects [START_REF] Coralic | Finite-volume weno scheme for viscous compressible multicomponent flows[END_REF], mass transfer effects in high speed cavitating systems [156] and surface tension forces [START_REF] Garrick | A finite-volume hllc-based scheme for compressible interfacial flows with surface tension[END_REF]177], but in these later works the authors have not investigated the the method to simulate bubble collapse processes. The Level Set approach has been also used by Nagrath et al [148], who use a Finite Element Formulation (FEM) to report simulations of collapsing bubbles solving the full Navier-Stokes equations including surface tension forces and viscous effects. In this case, the Ghost fluid method is used to avoid the numerical problems related to the jump of properties across the interface.

Geometrical VOF methods take especial care to compute the flux of α using the geometrical reconstruction of the interface at every timestep to limit the diffusion induced by the numerical discretization of the advection equation. The price to pay is their computational cost and the difficulties to develop high order advection schemes preserving the fundamental advantages of the method (e.g. mass conservation at the discrete level in the incompressible limit). Except few exceptions [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF] these methods have not been widely developed yet in the context of cavitation. Instead algebraic VOF are becoming popular for the simulation of cavitating bubbles [145,122,121,214]. These methods are very similar to DIM and require to manipulate the fundamental transport equation for the color function to introduce terms that limit the diffusion of the interface. The impact of these artificial terms on the accuracy of the results obtained is not fully clarified yet.

VOF methods are usually combined with pressure-based solvers where density is replaced by pressure as a primitive variable that is obtained solving for Eq. 7. Some interesting properties of these methods include that: (i) it is possible to easily combine them with the VOF method, allowing for a sharp representation of the interface position; (ii) it is straightforward to implement advanced methods for surface tension extensively developed and tested in problems where the fluid can be considered as an incompressible substance; and (iii) they are not restricted by the acoustic CFL condition for stability at expenses of degrading for lengths smaller than a critical wavelength controlled by the time step and the speed of sound of the fluid [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF].

The classical version of pressure-based solvers do not guarantee exact mass, momentum and energy conservation due to discretization errors, which significantly impact the validity of the results when unresolved structures appear. In addition, because these methods do not guarantee energy conservation at the discrete level, they cannot capture the correct speed of shock Fig. 4: Numerical simulation of a bubble collapsing near a wall using the solver presented in Fuster & Popinet [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF]. Interface contour and pressure field in the middle plane.

waves, which limits the applicability of these solvers to low Mach numbers. These limitations have motivated the development of alternative formulations [207,104,[START_REF] Denner | Pressure-based algorithm for compressible interfacial flows with acoustically-conservative interface discretisation[END_REF][START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF] that keep the density and the energy as a primitive variables still solving for Eq. 7 to obtain the pressure that is used to update both momentum and the total energy after the advection step. This requires to solve an implicit equation for pressure, which significantly increases the computational cost with respect to explicit solvers. These solvers have been shown to capture the correct shock speed of supersonic flows, they are especially suitable to investigate problems where one of the phases is weakly compressible or even fully incompressible substance, and satisfy basic relations such as the balance between pressure and capillary forces in an equilibrium state [START_REF] Abu-Al-Saud | A conservative and wellbalanced surface tension model[END_REF] which is difficult to achieve in other formulations including surface tension forces. They are also very roboust and compare favorably with experimental results of the collapse of a bubble close to a wall as well as the theoretical predictions of the Keller-Miksis equation for linear, weakly non-linear and moderately intense oscillations [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF]. Figures 45show the potential of this solver to investigate the formation of the jetting phenomenon as well as the generation of vortex ring structures associated to the non-spherical bubble collapse process. Fig. 5: 2D cut of the void fraction (right) and vortex structures defined according to the λ 2 < 0 criterion (left) colored with the velocity normal to the wall after the primary collapse of a bubble near a wall depicted in Figure 4.

Modeling bubble clusters

5.1 Bubble clusters in the diluted limit

Models for diluted systems

Modeling bubble clusters is not an easy task due to the difficulties related to the representation of general multiphase flows in problems where it is not possible to capture all the scales of the flow. A subset of models assume that the bubble concentration is low, α g 1, finding equations for a system containing a known number of spherical bubbles per unit volume n and known initial bubble radius distribution where the void fraction can be obtained as α g = 4/3πnR 3 . Because in the dilute limit bubbles remain relatively far from each other it is reasonable to assume that a Rayleigh-Plesset type equation accurately reproduces the oscillatory response of bubbles at the local scale using the far field pressure obtained from the averaged equations. Note that because the averaged inter-bubble distance d scales as d ∼ R 0 α -1/3 g the threshold for the validity of the diluted limit assumption is expected to reach relatively large void fractions (α g < 0.1). In addition, bubbles remain spherical if the characteristic Weber number (We = ρ l U 2 R0 σ ) is small although DNS simulations and experimental studies have shown that the Rayleigh-Plesset model still provides a reasonable assumption of the averaged gas volume variation even in the presence of strong assymetries induced by the presence of jets which leads to the conclusion that the response of bubbles smaller than 1mm are probably well captured by Rayleigh-Plesset based models.

In the particular case of an invisicid bubbly liquid the continuity and momentum equations in the liquid (Eqs. 1-2) simplify to the Bernoulli equation

∂φ ∂t + 1 2 (∇φ) 2 + p l -p 0 ρ l,0 = 0. ( 43
)
This equation can be used to solve the problem of N oscillating bubbles in a liquid where the forcing background potential Φ ∞ is known. Because the potential is a linear function, it is possible to decompose the potentials as a sum of the background potential and the sum over the potential emitted by every single bubble in the system and to write a system of N Bernoulli equations evaluated at the interface of each individual bubble to obtain the dynamic response of all bubbles. A particular solution of this problem for the limiting case of an incompressible liquid in absence of relative motion between bubbles is

R i Ri + 3 2 Ṙ2 i = p i -p ∞ ρ l + I * , (44)
I * = - N j =i Rj R 2 j + 2 Ṙ2 j R j d ij , (45)
which recovers the traditional Rayleigh-Plesset equation when the bubble is isolated. A generalized solution including for liquid compressibility effects is derived in [START_REF] Fuster | Modeling bubble clusters in compressible liquids[END_REF].

Harkin et al [START_REF] Harkin | Coupled pulsation and translation of two gas bubbles in a liquid[END_REF] present particular solutions of this problem taking into account the relative motion between bubbles to investigate the coupled pulsation and translation of oscillating bubbles and Ilinskii et al [103] present general solutions of this approach for large bubble clusters using both Hamiltonian and Lagrangian mechanics. One interesting conclusion from this model is that if we define an effective pressure p eff = p ∞ -ρ l I * , the influence of farther away bubbles (which increase as r 2 ) collectively exert greater influence on I * than the influence of a single bubble nearby, whose effect decreases as 1/r [103]. That means that the averaged disturbance induced by all bubbles present in the system on the pressure felt by a bubble is usually more important than any pressure disturbance induced by a nearby bubble, which is a key assumption of the averaged models described below. Although these models have been sucessfully used to investigate some collective effects [START_REF] Bremond | Interaction of cavitation bubbles on a wall[END_REF]215] they are problematic and computationally intensive when N is large.

Eulerian-Eulerian models based on volume averaged equations (Eqs. [START_REF] Coralic | Finite-volume weno scheme for viscous compressible multicomponent flows[END_REF][START_REF] Coste | Sound propagation in a liquid with vapour bubbles[END_REF] were initially used in the context of wave propagation in bubbly liquids by Mallock [141], Campbell & Pitcher [START_REF] Campbell | Shock waves in a liquid containing gas bubbles[END_REF] and Van Wijngaarden [195,196].

In the dilute limit it is justified to assume that ρ ≈ ρ l (1 -α g ) and following a similar procedure to that used to derive Eq. 7 it is possible to writte an equation for the averaged pressure

1 ρ l c 2 l Dp Dt + ∇ • u = ∂α g ∂t . (46)
The system of equations is closed adding the averaged momentum equation (Eq. 40) and computing α g = 4/3πnR 3 , where the evolution of the bubble radius is given by the solution of the Rayleigh-Plesset equation. Calflish et al [START_REF] Caflisch | Effective equations for wave propagation in bubbly liquids[END_REF] used a multi-scale approach to show that this model is O(α 2 g ) accurate in problems where we can neglect advection terms in order to be consistent with the assumptions adopted.

Later works used the concept of ensemble averaging developed in the context of turbulent flows defining averaged quantities based on probabilistic functions (some examples of these formulations include [START_REF] Dopazo | On conditioned averages for intermittent turbulent flows[END_REF][START_REF] Dopazo | Intermittency in free turbulent shear flows[END_REF]). The main advantage of these methods is that, as it occurs in turbulent flows, they allow defining averaged quantities for properties defined in the bulk from an inifinite number of realizations. Zhang and Prosperetti [216,217] rigorously derive the ensemble averaged equations accounting for the relative motion between bubbles and liquid obtaining equations that are shown to be O(α g ) accurate for non-linear problems. Some other differences respect to the Van Wijngaarden's formulation is that they relax some important assumptions adopted by allowing for a relative motion between the phases and capturing the influence of fluctuations on the bubble translational velocity.

The solution of averaged equations in non-linear regimes require the use of numerical methods to solve both, the dynamic oscillation of the bubbles and the averaged equations at the large scales. Kubota [123] proposed to use this model for the prediction of cavitation activity around hydrofoils treating the liquid as an incompressible substance and assuming no slip between the two phases. Bubble dynamic effects are simply modeled by solving the Rayleigh-Plesset equation and setting p b = p v (T 0 ) constant along the simulation. Taking the averaged equations and using the Rayleigh-Plesset equation to relate the changes on the void fraction with pressure Kubota finds an equation for pressure where a source term depends on the dynamic response of the bubbles modelled using the Rayleigh-Plesset model. A similar model proposed by Preston et al [161] has shown how bubble dynamic effects are very important in non-linear regime, significantly reducing the range of applicability of barotropic models applied to the averaged quantities (e.g. density and pressure). This stresses an important point: the fact that both gas and liquid obey a barotropic relation does not imply that the averaged mixture response is barotropic too.

Various authors [START_REF] Ando | Numerical simulation of shock propagation in a polydisperse bubbly liquid[END_REF]133,134] report numerical simulations of the the ensemble averaged equations for polydisperse mixtures where Eq. 38 is particularized for the case of a system containing spherical bubbles with concentration α g ,

∂α g ∂t + ∇ • (uα g ) = 3α g R 2 Ṙ R 3 . ( 47
)
For polydisperse mixtures the method requires the solution of the RP equation for the full PDF distribution of bubble radius in each computational grid, which is an expensive computational operation for disperse bubble clouds. Alternatively Fuster & Colonius [69] propose an Eulerian-Lagrangian formulation based on the volume-averaged approach to solve the equations for the contiuum phase (e.g. the liquid) from the averaged equations for the mixture (Eqs. [START_REF] Coralic | Finite-volume weno scheme for viscous compressible multicomponent flows[END_REF][START_REF] Coste | Sound propagation in a liquid with vapour bubbles[END_REF][START_REF] Coutier-Delgosha | Numerical simulation of the unsteady behaviour of cavitating flows[END_REF]. Thus, the influence of the disperse phase is treated as source terms in the continuity, momentum, and energy equation for the liquid that for α g < 1 become [START_REF] Fuster | Modeling bubble clusters in compressible liquids[END_REF][START_REF] Colonius | Investigation of a new model for bubbly cavitating flow[END_REF]

∂ρ l ∂t + ∇ • (ρ l u) = ρ l ∂α g ∂t + u • ∇α g ( 48
)
∂ρ l u ∂t + ∇ • (ρ l uu) + ∇ • (pI) = ρ b u ∂α g ∂t + u • ∇α g ( 49
)
∂ρ l e l ∂t + ∇ • (uρ l e l ) + ∇ • (up) = ρ l e l ∂α g ∂t + u • ∇α g (50)
To compute α g , we multiply each individual bubble volume by a discrete kernel that smears its value onto the computational grid, and, at each point in the domain, sum this over all bubbles. The evolution of the volume of each individual bubble is obtained from the solution of a RP-like equation which is a generalization of the incompressible bubble-cluster model of Ilinskii et al. [103] to compressible liquids, and likewise assumes spherical bubbles and potential flow in the vicinity of the bubbles. By taking a cell based approach this model is grid dependent, much in the same way a Large Eddy Simulation of turbulence depends on the grid spacing. As the grid is refined to the size of an individual bubble, the model exactly represents the dynamics of an isolated spherical bubble. At the other extreme, the model is demonstrated to converge to the solution to ensemble-averaged models for low spatial resolution with much less computational expense compared to the corresponding ensemble-averaged approach. The interaction between diluted bubble clusters and shock waves have been investigated in recent works by Maeda et al [140] comparing favoroubly for problems related to medical applications. Models based on similar ideas have been also proposed to investigate the attenuation of waves by bubble clusters to protect offshore structures and sealife [START_REF] Grandjean | Shock propagation in liquids containing bubbly clusters: a continuum approach[END_REF] and problems related to hydrodynamic and ultrasonic cavitation [138] in regions where the gas/vapor phase can be considered disperse. One interesting observation of the results provided by these models is the prediction of significant pressure and velocity fluctuations that can be directly attributed to the heterogeneous distribution of the disperse phase in a given situation. These effects are completely neglected in ensemble averaged models that only provide information about the averaged response of the system over an infinite number of realizations.

Lagrangian-Eulerian approaches are appealing in order to account for the relative motion between the fluid and the bubbles. Indeed these models have been already used to capture the dynamics of bubble clusters where translation effects are relevant. Figure 6 shows a typical example of the response of a cylindrical bubble cloud obtained with the numerical model proposed in [192] coupled with the Epstein & Plesset model to predict the bubble volume growth [START_REF] Ps Epstein | On the Stability of Gas Bubbles in Liquid-Gas Solutions[END_REF] in order to mimick the observations of bubble clusters rising in supersaturated liquids [168]. Interestingly this very simple model is already able to capture the appearance of large scale flow structures generated by the collective motion of the bubbles showing strong coupling between the transport of gas from the liquid bulk to the bubble interior and turbulence in an auto-catalytic process that triggers the explosive growth of the bubble cluster. However these first attempts to reproduce real observations also point out the limitations of currently available models to correctly capture the correct growth of the bubble cloud and the terminal rising velocity of the cluster observed experimentally. One of the reasons is the strong simplifications adopted in the model, which neglects the influence of advection on the mass transfer flux. As we will discuss below, this is also one of the main limitation of cavitating flow models.

Linear solutions for diluted bubble clusters

The linearized version of the volume-averaged problem has been used by D 'Agostino and Brennen [44,[START_REF] Agostino | Linearized dynamics of spherical bubble clouds[END_REF] and Omta [152] to obtain solutions for a bubble cluster of finite size characterized by the cluster radius R c containing a monodisperse distribution of bubbles of size R 0 . One important result is that, unlike the case of a single bubble where a unique resonance frequency appears, bubble clusters exhibit multiple resonant frequencies given by

ω n = ω 0 1 1 + 3αg,0 (n-1/2) 2 π 2 Rc R0 2 n = 1, 2, .... ( 51
)
The highest resonant frequency corresponds to the single bubble natural frequency, while the lowest resonant frequency,

ω 1 ≈ ω 0 R0 Rc π √ 12αg,0
, is a characteristic resonance of the overall bubble cluster. These two frequencies can be used to discriminate the various frequencies appearing in a cavitating flow depending if they lay below ω 1 , above ω 0 or in-between.

Another fundamental problem widely investigated with ensemble-averaged models is the linear wave propagation in a bubbly liquid [START_REF] Campbell | Shock waves in a liquid containing gas bubbles[END_REF]195,196,[START_REF] Crespo | Sound and shock waves in liquids containing bubbles[END_REF][START_REF] Commander | Linear pressure waves in bubbly liquids: Comparison between theory and experiments[END_REF]. For a problem with a known distribution of bubble radius f (a) obeying the Rayleigh-Plesset equation, Eq. 46 takes the particular form

1 c 2 l ∂ 2 p ∂t 2 -∇ 2 p = 4πρ ∞ 0 a 2 Rf (a)da, ( 52
)
where, for a weak perturbation, R ≈ -ω 2 R 0 R . Because the radius and external pressure perturbation are connected through the linearized solution of the Rayleigh-Plesset equation (Eq. 31) it is possible to rewrite Eq. 52 as a wave equation for p,

1 c 2 eff ∂ 2 p ∂t 2 -∇ 2 p = 0 ( 53
)
where the effective sound speed, c e , is a complex quantity given by

1 c 2 eff = 1 c 2 l + 4πn ∞ 0 af (a)da ω 2 0 -ω 2 + 2iδω (54)
In most situations the effective sound speed is mainly governed by the dynamic response of the bubbles. In particular, for monodisperse bubble clouds and frequencies smaller than the bubble resonance frequency, ω ω 0 , the sound speed of the mixture simplifies to [START_REF] Didenko | Effect of noble gases on sonoluminescence temperatures during multibubble cavitation[END_REF] which is equivalent to the widely used low frequency limit [205] for pure gas bubbles in the diluted limit

c eff ≈ p 0 α g ρ l (
1 ρc 2 eff = α g ρ g c 2 g + 1 -α g ρ l c 2 l ≈ α g p 0 . ( 56
)
General solutions accounting for frequency effects have been reviewed by Commander & Prosperetti [38] finding good agreement between previously reported experimental results and theoretical predictions for the case of pure gas bubbles. This limiting solution can be used to validate numerical codes [START_REF] Ando | Effects of polydispersity in bubbly flows[END_REF][START_REF] Fuster | Modeling bubble clusters in compressible liquids[END_REF][START_REF] Fuster | Effect of direct bubble-bubble interactions on linear-wave propagation in bubbly liquids[END_REF][START_REF] Doc | Nonlinear acoustic propagation in bubbly liquids: Multiple scattering, softening and hardening phenomena[END_REF]. Extensions of the theory to account for the effect of direct bubblebubble interactions have been numerical reported in [START_REF] Fuster | Effect of direct bubble-bubble interactions on linear-wave propagation in bubbly liquids[END_REF].

Classical theories accounting for mass transfer effects assume that the system's temperature and pressure fluctuations are related through the equilibrium curve [126]. However, those approaches have difficulties reproducing experimental measurements in the dilute limit [START_REF] Coste | Sound propagation in a liquid with vapour bubbles[END_REF], where the liquid temperature is probably not significantly affected by the presence of the vapor bubbles. Recently, general solutions accounting for mass transfer effects across the interface have been found by [START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF]164]. The theoretical analysis shows that there exists a low frequency limit solution that corresponds to Eq. 56 for noncondensable gas bubbles. However, as we approach the saturation line and phase change effects become important the low frequency limit for the phase speed significantly decreases while the effective attenuation increases with respect to the case of pure non-condensable gas bubbles. Another important effect of mass transfer effects is that ω 0 decreases drastically as the vapor content inside the bubble increases (Eq. 32), reducing the range of applicability of the low frequency limit solution [START_REF] Fuster | Mass transfer effects on linear wave propagation in diluted bubbly liquids[END_REF].

Modeling bubble clusters containing large gas/vapor concentrations

While accurate and consistent in the dilute limit, some problems involving cavitating flows require the general solution of a multiphase flow problem where the dynamic of the gas/vapor pockets is not well represented by the solution of the Rayleigh-Plesset equation. In problems where the gas/vapor phase becomes predominant but it is still not possible to solve all the length scales of the flow it is compulsory to develop models that must versatile enough to deal with the process of inception and bubble collapse in regions where the interface is not well resolved as well as to transport a well resolved interface in those regions where gas/vapor bags are much larger than the grid size. The problem becomes specially challenging if one considers that the limited resolution makes compulsory to model the interaction between the interface, mass transfer effects and turbulence and that large scale gas/vapor structures coexist with small bubbles as it is typically the case in applications involving cavitation around hydrofoils, cavitation in nozzles, etc... Eulerian-Eulerian methods have been mainly proposed to solve the equations for the averaged mixture, being the evolution of the gas/vapor fraction a key ingredient in these models. The evolution of the gas/vapor fraction in each cell is obtained from the solution of Eq. 37 in the quasi-conservative form ∂α ∂t

+ ∇ • (uα) = α∇ • u + (u I -u) • nδ(x I ) ( 57
)
where the two terms appearing on the Right Hand Side (RHS) of Eq. 57 represent the change of void fraction due to the divergence of the averaged field and the contribution to the net gas fraction change due to mass transfer. Note that even in absence of a net flux across the interface, gas compression or expansion can be responsible of a divergent velocity field u. The difficulties to correctly solve the RHS of Eq. 57 in problems involving very different spatial scales have lead to various authors to group and model the effects of mass transfer (evaporation/condensation) and the gas/vapor expansion into a single term

∂α g ∂t + ∇ • (uα g ) = ṁ. (58)
A popular model to compute ṁ is proposed by Schnerr & Sauer [178] who obtain the evolution of the void fraction using

∂α g ∂t + ∇ • (uα g ) = 3α g R dR dt . (59)
The form of this equation is similar to the one derived for disperse bubbly flows, but in this model the bubble radius velocity is obtained in the inertial controlled regime given by Eq. 20. The Schnerr & Sauer model is often expressed in the form of Eq. 58 as

ṁ = 3α g (1 -α g ) R b ρ g ρ l ρ m sign(p v -p) 2 3 |p v -p| ρ l (60)
which, unlike the models proposed above for disperse flows, allows computing the evolution of α g even when the void fraction tends to one if one gives an estimation of the bubble radius. One of the problems of the initial Schnerr model is that it neglects all thermal and mass transfer effects on the dynamic response of the bubble, it assumes that the inertial growth model (Eq. 20) is applicable even for large gas/vapor fractions, and it requires to sed the liquid with initial bubble nuclei whose characteristic size and concentration is typically difficult to control and also to measure experimentally. Although some authors have proposed some remedies to account for thermal effects on the dynamic response of the bubbles by using the Plesset solution already

ṁ+ ṁ- [144] C prod ρ l ρv α l max(0,pv -p) 1/2ρ l U 2 ∞ t∞ C dest αg max(0,p-pv ) 1/2ρ l U 2 ∞ t∞ [124] C prod ρv ρ l α 2 l (1-α l ) t∞ C dest ρv ρ l α l max(0,p-pv ) 1/2ρ l U 2 ∞ t∞ [172] Ccα 2 g (1 -αg) 2 ρ l ρg max(0,pv -p) √ 2πRg T Ceα 2 g (1 -αg) 2 max(0,p-pv ) √ 2πRg T
Table 1: Resume of popular models for the computation of the gas/vapor production and consumption terms.

mentioned above [210], the main limitations of the method remain.

Other models including those proposed by Merkle [144], Singhal et al [183], Kunz et al [124,125] and Saito [172] propose to express ṁ as a competition between a production and a destruction term, merging expansion, compression and phase change processes into a single semi-empirical function that depends on the pressure difference between the actual pressure and the vapor pressure (see Table 1 for an explicit definition of the production and consumption terms). Effectively these models act as the transient mass transfer flux used in single bubble models and tend to keep the interface conditions close to equilibrium conditions. The sensitivity of the results on the details of the model remain a major challenge given the difficulties to calibrate the model in different regions where different mechanisms may control the bubble response. Only for those models including thermal and mass diffusion effects [172] one may expect results to be less sensitive to the expressions given for ṁ given that, as it occurs for the case of single bubbles described above, in a general problem the overall mass flux will be determined by the slowest process between the thermal transport in the surrounding boundary layer, vapor diffusion inside the bubble, and the instantaneous phase change flux modeled through the ad-hoc expression for ṁ.

There exists also models that do not directly solve for the transport equation of α g . Instead, a barotropic EOS is used to relate the averaged density and pressure to implicitly define the vapor fraction as α g = ρ-ρ l ρg-ρ l , where ρ g and ρ l are assumed to be known and constant values (see for example [START_REF] Coutier-Delgosha | Numerical simulation of the unsteady behaviour of cavitating flows[END_REF]). Particularly simple is the approach proposed by Scherr et al [179] who propose to solve the modified Riemann solver and to evaluate the properties whether the averaged density is above or below the saturation vapor density at a given pressure and temperature. These very simple approaches have a number of practical advantages, including the fact that they handle naturally the processes of cavitation inception and mass transfer which are difficult to model. The drawback is that they assume that the system's temperature is uniform and lack of fundamental background for the definition of a rigorous barotropic EOS, on which the results obtained are going to be extremely sensitive.

In addition to the advection equation we need to solve for the mixture conservation equations (Eqs. [START_REF] Coralic | Finite-volume weno scheme for viscous compressible multicomponent flows[END_REF][START_REF] Coste | Sound propagation in a liquid with vapour bubbles[END_REF][START_REF] Coutier-Delgosha | Numerical simulation of the unsteady behaviour of cavitating flows[END_REF]. Various authors solve for the continuity and momentum equation for the mixture assuming both fluids as incompressible substances and neglecting any influence of temperature on the process [124,125,180,197] where a turbulence model is usually added to model the influence of subgrid turbulence.

Models including thermal and compressibility effects have been mainly developed to simulate shocked flows in presence of cavitation and cavitation in cryogenic liquids where thermal effects are known to be relevant.

Regarding models for cryogenic cavitation Tseng et al [194] and Zhu et al [220] solve for the enthalpy equation assuming both gas and liquid as incompressible substances while Gnanaskandan [START_REF] Gnanaskandan | A numerical method to simulate turbulent cavitating flows[END_REF] solves for the internal energy equation to show results on the cavitating flow around a cylinder [START_REF] Gnanaskandan | Large eddy simulation of the transition from sheet to cloud cavitation over a wedge[END_REF][START_REF] Gnanaskandan | Numerical investigation of near-wake characteristics of cavitating flow over a circular cylinder[END_REF]. Saito et al [172] and more recently Goncalves and coworkers [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF][START_REF] Goncalves | Modeling for non isothermal using 4-equation models[END_REF][START_REF] Colombet | On numerical simulation of cavitating flows under thermal regime[END_REF] use formulations that include the equation for the total energy of the mixture (derived from Eq. 4). Most of these models [194,220,[START_REF] Gnanaskandan | A numerical method to simulate turbulent cavitating flows[END_REF]172] still resort to the calculation of the mass transfer flux through the semi-empirical models described above. Only the recent work of Colombet et al [START_REF] Colombet | On numerical simulation of cavitating flows under thermal regime[END_REF] makes a first attempt to impose the condition obtained from the energy balance across the interface (Eq. 12) to obtain the total mass transfer flux from the heat flux at the interface and the enthalpy of vaporization neglecting any temperature gradient inside the gas/vapor phase. Another characteristic of these models is that they require the definition of an EOS for the mixture that defines properties such as the effective sound speed. One possibility is to define an EOS for the mixture using the the averaged pressure and the EOS of each single substance. Goncalves & Patella [86] compares the use of this methodology with an ad-hoc barotropic EOS where the averaged density and pressure are directly related through the barotropic relation proposed in [START_REF] Dellanoy | Two phase flow approach in unsteady cavitation modeling[END_REF]. Remarkably, the authors find that the barotropic EOS performs better than the mixture of EOS for cavitation in standard liquids [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF][START_REF] Goncalves | Modeling for non isothermal using 4-equation models[END_REF] while the barotropic EOS model fails reproducing experimental results when thermal effects are relevant [START_REF] Colombet | On numerical simulation of cavitating flows under thermal regime[END_REF]. This is consistent with the conclusions extracted from the linear oscillation of vapor bubbles presented above where we have concluded that in those regions where transient heat transfer effects control the dynamic response of bubbles a barotropic EOS will have troubles to correctly reproduce the response of the gas/vapor mixture. The interest of these models may extend beyond cryogenic liquids, given that heat diffusion effects are certainly important in all problems where the heat flux controls the overall mass transfer rate in standard liquids in those regions where the gas/vapor pressure remains approximately constant and equal to the vapor pressure.

For high speed cavitating flows we face similar problems to those already described for the DNS of single bubbles (Section 4.3.2). Thus the quality of the solution obtained is strongly determined by the numerical methods used to solve the physical equations and models including the explicit solution of the energy equation, which are critical to correctly predict the appearance of shocks and to capture the correct dynamic response of the interface. The correct choice numerical method becomes especially important in situations where the Mach number is large, which can occur at relatively low velocities when mass transfer effects are important and the effective sound speed of the mixture reaches very low values. The formulation proposed by Saurel and Agbrall [173] already discussed for the case of DNS simulations has been widely adopted in the literature [114,174,176,154] but the uncertainties about the accuracy of these models remain [175].

Turbulence modeling

Modelling turbulent cavitating flows remains as a major challenge due to the difficulties to extend the classical turbulence models extensively used and developped in the context of single phase flows simulations. The two main differences with respect to single phase flows are (i) the generation of turbulence due to the volume changes induced by gas/vapor compressibility and phase change and (ii) the generation of turbulence induced by the presence of an interface itself.

The source of volume generated in a cavitating flow usually leads to large expansion/compression velocities that induce significant turbulence at small scales. Note that the vapor density at vapor pressure is extremely low and therefore the volume of vapor generated is typically of the order of one million times the volume of liquid vaporized.

In the diluted limit bubbles can be seen as microscale agitators that induce velocity perturbations in the flow due to the flow pressure variations. In the linear regime the non-dimensional velocity fluctuations due to volume variations scale as u ∝ ωR0 U∞ R , where U ∞ is a characteristic flow velocity of the process and R is the non-dimensional bubble radius perturbation that is related to the external pressure fluctuations through Eq. 31 as R = F (ω)p . Thus, compared to a single phase flow, bubbles introduce additional velocity fluctuations that in the linear regime can be obtained as u ∝ ωR0 U∞ F (ω)p being especially important when the forcing frequency is close to the bubble resonance frequency, when F (ω) takes very large values. It is also interesting to note that F (ω) at low frequencies is F (ω → 0) ≈ 1 3(1-Y0) (Eq. 33). This implies that as we approach the saturation curve and bubbles fill with vapor (Y 0 → 1) the perturbations induced by the bubble oscillation in the system become more important for a fixed amplitude of the external pressure variation. In highly-nonlinear and out-of-equilibrium systems the analysis becomes more complicated but the overall picture is not expected to change significantly. In addition, the increase on the turbulence level is known to lead to an increase on the transport coefficients and therefore on the effective bubble growth. The coupled feedback between the turbulence and mass transfer has been experimentally observed to lead an auto-catalytic process that induces an exponential growth on the gas fraction in liquids supersaturated with dissolved gas [168]. Diluted bubbly flows are also encountered during the first instants after nucleation and in the collapse stage, where bubbles are known to play a significant role on the effective increase of the turbulence level and the effective transport coefficients in the flow.

The analysis becomes more intrincate as vapor fraction increases given that different mechanisms counter-act to the natural increase of the turbulence level as the void fraction increases. For example in some applications related to hydrodynamic cavitation the increase on the effective drag and pressure losses induced by the presence of bubble clusters usually leads to a decrease on the effective velocity inside the bubble cluster, tending to limit the increase of the effective Reynolds number. In ultrasonic applications the significant damping of pressure waves induced by highly concentrated bubble clusters significantly attenuate the amplitude of pressure perturbations and therefore the level of turbulence from the oscillation of bubbles is expected to saturate or even decrease. One of the major challenges in cavitating flow models is then to be able to capture the balance between mechanims inducing the gas/vapor fraction growth and those limiting it, which involves the resolution of both the large scales associated to the overall structure of the bubble cluster and the small structures developing inside the multiphase flow and probably controlling the dissipation mechanisms.

In addition to the volume source induced by the presence of a highly compressible substance undergoing phase change, the presence of an interface itself influences the type of structures present in the flow as well as the dissipation mechanisms. The development of physical models capturing these effects at the subgrid scale are scarce due to the intrinsic challenges when modelling multiphase flows and the difficulties to validate them with reliable DNS and experimental data. To get further insight about the strong coupling between the dynamic response of an interface and turbulence let's consider the generation of vorticity in an academic two-dimensional incompressible multiphase flow. In this case the evolution of the total vorticity in a given system can be expressed as [136,[START_REF] Brøns | Vorticity generation and conservation for two-dimensional interfaces and boundaries[END_REF]

d dt ωdA = ν∇ω • ndS - b a S ω ds, (61)
where S ω is the density of a vorticity source located at the interface defined in absence of tangencial velocity jump as

S ω = ∂ ∂s p ρ + κ (u • n) (u • t). (62)
The first term contains the generation of vorticity due to variations on the surface tension forces (e.g. curvature) along the bubble interface, normal viscous stresses and mass transfer effects (these two last effects usually play a minor role). The second term depends on the interface curvature κ, the tangential velocity to the interface and the normal velocity jump. Because the normal velocity jump is directly related to the mass transfer flux across the interface, this expression reveals that mass transfer effects have an impact on the generation and destruction of vorticity in a turbulent multiphase flow. Thus, the presence of an interface significantly modifies the generation of turbulence in a given system even in the incompressible limit, being surface tension and mass transfer effects two key mechanisms that significantly influence the characteristics of the overall flow.

How to model the interaction between turbulence and an interface and the impact of these effects on the overall averaged quantities is not easy given the difficulties to propose theoretical models based on physical principles but also due to the significant and usually unknown dissipation introduced by a given numerical method when simulating a turbulent multiphase flow problem [START_REF] Fuster | An energy preserving formulation for the simulation of multiphase turbulent flows[END_REF]. We summarize next some of the approaches that have been reported until now to model turbulence in cavitating flows.

One solution is to use the natural diffusion of numerical schemes as an implicit turbulence model which ideally captures the influence of the subgrid scales on the resolved structures. Using this approach, typically known as Monocally Integrated LES (MILES) [START_REF] Fureby | Monotonically integrated large eddy simulation of free shear flows[END_REF][START_REF] Fureby | Large eddy simulation of highreynolds-number free and wall-bounded flows[END_REF], there is no need to include any physical model capturing the influence of unresolved structures on the numerical solution. In the context of cavitating flows some examples include those in [178,153,179,[START_REF] Bensow | Implicit les predictions of the cavitating flow on a propeller[END_REF]102,[START_REF] Colombet | On numerical simulation of cavitating flows under thermal regime[END_REF].

Another classical alternative is the use of RANS models proposed as simple extension of already existing models for single phase flow. Various authors [124,125,184,144,[START_REF] Ahuja | Simulations of cavitating flows using hybrid unstructured meshes[END_REF]200,[START_REF] Andriotis | Vortex flow and cavitation in diesel injector nozzles[END_REF]194,220]

∂ρ m k ∂t + ∇ • (ρ m ku) = ∇ • µ t P r t,k ∇k + P -ρ ∂ρ m ∂t + ∇ • (ρ m u) = ∇ • µ t P r t, ∇ + [C 1 P -C 2 ρ ] k
This model has been shown to provide reasonable approximations of the wall pressure distribution in steady-state and transient sheet-and super-cavity flows [124,125,200] but the density and stress distributions obtained typically exhibit a high sensitivity to the model parameters. Saito [172] computes the turbulent viscosity using the modification of the Baldwin-Lomax model [START_REF] Baldwin | Thin-layer approximation and algebraic model for separated turbulentflows[END_REF] proposed by . This model requires the definition of a turbulent heat conductivity κ t defined as κ t = µ t /P r t . The model is validated for cavitating flow around a body of hemi-sphere/cylinder geometry and by numerical simulations of 2D and 3D unsteady cavitating flow around a NACA0015 hydrofoil.

The k-ω SST model proposed by Menter [142], where the turbulence kinetic energy and specific dissipation rate are computed from the solution of two PDEs, has become also popular to simulate hydrodynamic cavitating flows [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF]169,155]. In particular the recent work of Pendar et al [155] shows interesting comparisons between the predictions of this model and the experimental results of [START_REF] Brandner | An experimental investigation of cloud cavitation about a sphere[END_REF] and LES simulations, showing that the model is able to provide time-averaged values of the unsteady vapor shedding of the cavity behind the sphere, but it is unable to capture the instabilities of the cavity accurately.

Despite the improvements on the quality of numerical results it is generally admitted that it is difficult to reproduce experimental observations using RANS models given that they neglect any influence of the multiphase nature of the flow in the turbulence model. Coutier-Delgosha et al [START_REF] Coutier-Delgosha | Numerical simulation of the unsteady behaviour of cavitating flows[END_REF] propose a modified k-model based on the correction of the turbulent viscosity, such that µ t = f (ρ)C µ k 2 / , with C µ = 0.085 and f (ρ) = ρ g + α n g (ρ l -ρ g ) n > 1. This empirical correction is based on the observation that the standard kmodel over-estimates the viscous dissipation introduced in cavitating flows supressing physical cyclic behaviors observed in a venturi section. Similar findings have been reported by other authors (see for example [START_REF] Gnanaskandan | A numerical method to simulate turbulent cavitating flows[END_REF][START_REF] Gnanaskandan | Numerical investigation of near-wake characteristics of cavitating flow over a circular cylinder[END_REF]). Although the definition of a modified viscosity circumvent the problem for a specific application, this fact reveals the problems of single phase turbulence models to correctly capture the physical mechanisms induced by the interaction between an interface and a turbulent flow already discussed.

Due to the difficulties of RANS models to reproduce unsteady cavitating flows mechanisms, in recent years many researchers have been testing LES models obtaining a significant improvement in the predictions compared to RANS models. Some examples include simulations around a sphere [START_REF] Gnanaskandan | A numerical method to simulate turbulent cavitating flows[END_REF][START_REF] Gnanaskandan | Numerical investigation of near-wake characteristics of cavitating flow over a circular cylinder[END_REF]155] and a disk [169], the flow around an hydrofoil [199,170,106,118,102] and simulations of the flow generated by a propeller [START_REF] Bensow | Implicit les predictions of the cavitating flow on a propeller[END_REF]. In these models the subgrid scale turbulent viscosity, µ k , is obtained by a "Local Eddy-Viscosity" model in which the influence of the presence of multiple phases on the problem is still neglected. Most of these works are only straightforward extensions of the methodologies followed for the simulation of single phase flows including the classical Dynamic Smagonisky model [START_REF] Gnanaskandan | A numerical method to simulate turbulent cavitating flows[END_REF][START_REF] Gnanaskandan | Numerical investigation of near-wake characteristics of cavitating flow over a circular cylinder[END_REF], implicit LES methods [START_REF] Bensow | Implicit les predictions of the cavitating flow on a propeller[END_REF]102], and the "one equation eddy viscosity model" [170,118,155]. Among the various results it is worth mentioning the work of Pendar et al [155] who report LES simulations of supercavitation over a sphere at a constant Reynolds number of 1.5 × 10 6 and a broad range of cavitation numbers (0.36 < σ < 1) obtaining structures that are relatively close to those experimentally reported in [START_REF] Brandner | An experimental investigation of cloud cavitation about a sphere[END_REF] compared to those previously reported using RANS simulations.

Nucleation

Another important process influencing the characteristics of cavitating flows is the process of bubble nucleation. The amount of liquid vaporized or the rate of bubble nuclei activated control the characteristics of the entire cavitating flow but how to model these effects is still object of an intense research. This work is not intended to provide a review of nucleation theories but it is still worth mentioning some challenges when modeling these processes. Different theories and mechanisms have been proposed by various authors to explain the process of bubble nucleation [112,[START_REF] Caupin | Cavitation in water: a review[END_REF] but the experimental validation of these theories is difficult and in most of the applications it is not easy to discriminate among the various proposed mechanisms the one that eventually controls the process of nucleation in a given application.

Homogeneous nucleation models are designed to predict phase transition from pure substances that enter in a metastable state. These models are characterized by the prediction of the spinodal curve which establishes the limits of metastability beyond which the presence of single substance is unstable and phase separation occurs naturally. The problem is that the lack of consensus about how to theoretically obtain the spinodal curve and the difficulties to experimentally obtain the EOS in a metastable state makes difficult to validate various different approaches proposed in the literature [START_REF] Debenedetti | Metastable liquids: concepts and principles[END_REF]. To measure the cavitation threshold in water has been particularly challenging given that different experimental techniques provide different negative pressures beyond which bubble nucleation occur. Only using mineral inclusions scientists have been able to obtain values consistent with the Classical Nucleation Theory (-140 MPa) [219]. Other techniques including high intensity ultrasounds in ultrapurified and degassed water have revealed that water naturally nucleates at lower negative values being the interpretation of these results is still topic of an intense research [START_REF] Azouzi | A coherent picture of water at extreme negative pressure[END_REF].

Heterogeneous nucleation models justify the process of bubble inception by the growth of pre-existing bubble nuclei that are naturally present in the system [START_REF] Dean | The formation of bubbles[END_REF][START_REF] Harvey | Bubble formation in animals. i. physical factors[END_REF][START_REF] Harvey | Bubble formation in animals. ii. gas nuclei and their distribution in blood and tissues[END_REF]. In this type of models, continuum equations are still valid although the typical length scales required to capture the dynamic response of pre-existing germs are very small (below 1 µm). Blake's radius relates the size of the bubble nuclei in a liquid and the maximum tensions the liquid can stand before inception occurs for a spherical nucleus,

p c = p vap - 4σ 3R c , (63)
R c R 0 = √ 3 1 + p 0 -p vap 2σ R 0 1/2 , ( 64
)
where σ is the surface tension and R c is the critical radius for which the maximum negative pressure is reached. This result reveals the possibility to find bubbly liquids in mechanical equilibrium where the liquid is in a metastable state (e.g. under negative pressure). The expression above also allows obtaining the minimum nuclei size R 0,min above which pre-existing bubbles will nucleate for a known p min below the vapor pressure. Extensions of this reasoning have been proposed in the literature to capture the influence of the walls which is important to prevent the nuclei from dissolution [START_REF] Apfel | The role of impurities in cavitation-threshold determination[END_REF][START_REF] Atchley | The crevice model of bubble nucleation[END_REF]147] or polydispersity [START_REF] Fuster | Stability of bubbly liquids and its connection to the process of cavitation inception[END_REF].

All heterogeneous nucleation theories assume that the bubble nuclei are in mechanical equilibrium with the surrounding medium at every instant before inception occurs. In reality dynamic effects can be relevant in applications involving ultrasounds and shock waves where the excitation frequencies can be of the order or above the bubble resonance frequency of pre-existing germs, being the excitation frequency an important parameter in the process of nucleation [START_REF] Esche | Untersuchung der schwingungskavitation in flüssigkeiten[END_REF][START_REF] Crum | Acoustic cavitation thresholds in water[END_REF]. A rough estimation for the upper bound of the range of applicability of these approaches can be obtained by imposing that a bubble will only equilibrate its pressure with a varying external pressure fluctuation of the surrounding medium if the excitation frequency ω is smaller than the bubble resonance frequency, ω 0 ∼ 1 R0 p0 ρ l . This implies that among the overall distribution of bubble nuclei pre-existing in a liquid, only those within the range R 0,min (p min ) < R 0 < R 0,max (ω), where R 0,max (ω) ∼ 1 ω p0 ρ l , will fulfill the conditions for the Blake's estimation to hold valid as an activation criterion. This result must be taken cautiously as nucleation usually occurs close to the walls and the hypothesis of sphericity is questionable but it sheds light into the relevance of dynamic processes on the process of bubble nucleation.

Finally, another problem of heterogeneous cavitation models is that even when they have been widely used to qualitatively describe the process of bubble nucleation, to use them as a quantitative tool require to know parameters such as the size and geometry of pre-existing nuclei, which are expected to dynamically evolve in time and to depend on parameters such as the surface roughness that can be only controlled in laboratory designed experiments [START_REF] Bremond | Cavitation on surfaces[END_REF][START_REF] Bremond | Controlled multibubble surface cavitation[END_REF]. To predict in real applications the density of nucleation sites and the rate of depletion of these sites depends, among other factors, on the influx of bubble nuclei in the system, the influence of the walls roughness and hydrophobicity properties, the presence of impurities and mass diffusion effects that induce memory effects on the process of bubble nucleation [212,[START_REF] Caupin | Cavitation in water: a review[END_REF]. Some authors propose the use of Molecular Dynamic simulations to predict the natural formation of sub-micron bubbles in the flow [143] and recent works have investigated the mechanisms of nanobubble formation and stabilization in flat surfaces (see [132] for a review) but in practice there is not yet a clear strategy about how to develop and validate models able to reproduce the overall process of heterogeneous nucleation.

Conclusions and perspectives

The development of numerical methods and models for the simulation of bubble clusters described above is making possible to explore problems of critical importance in fundamental science and engineering applications. Current numerical techniques are a useful tool to improve our understanding about the intrincate mechanisms involving pressure variations and bubbles but certainly further developments are needed in order to capture the rich variety of physical phenomena usually taking place in applications involving fast bubble dynamic responses and cavitating flows.

Many of the models used for DNS of collapsing bubbles assume that the response of the bubble is adiabatic (high Peclet number regime) neglecting mass and heat transfer across the interface. These codes have shown to be a useful tool for the investigation of the interaction of pressure waves and bubbles, the peak pressures and shear stresses generated by collapsing bubbles near walls and to characterize the turbulent structures generated during the collapse of a bubble. However, they are still limited to regimes where heat and mass transfer effects are negligible. This is problematic as it has been shown in numerical studies accounting for mass transfer assuming spherical symmetry. Part of the reason why mass transfer models have not been included in three-dimensional Eulerian-Eulerian solvers is that mass transfer significantly increase the complexity of the numerical methods required to solve the equations. How to correctly impose the jump conditions established by the balances across an infinitely thin interface remain a challenge in these methods. In addition Riemann solvers and pressure based solvers require expressions for the sound speed in cells containing liquid and vapor mixtures, for which there are not well-established expressions when phase change effects are relevant. The development of numerical codes able to handle phase change process across the interface will require the validation of codes in problems related to cavitating flows where mass transfer effects become relevant either by using well-controlled experiments or theoretical solutions.

Another fundamental difficulty on the DNS of real cavitating bubbles is how to model interfacial phenomena. The bubble drag and mass transfer processes across the interface are only two examples of known processes to be influenced by the presence of solid impurities and surfactants accumulating at the interface but how to model these effects in continuum models is an important modelling issue that has not been solved yet. The dynamics of the triple contact point is also a major modelling issue for the DNS of the process of bubble dynamics close to walls which has implications understanding the process of heterogeneous nucleation and also reproducing the dynamics of bubbles in contact with solid walls. Models for large bubble clusters in the dilute limit have been shown to relatively well capture various observations such as the propagation of acoustic and shock waves in bubbly liquids. However, modeling cavitating flows at large void fractions still face important modeling issues that have not been correctly addressed yet. For instance a large portion of models assume that the mixture response can be modeled using a barotropic equation of state, while numerical studies for single bubbles and also for bubble clusters in the dilute limit show that mechanisms such as the dynamic response of the bubbles and the heat and mass transfer effects through the interface usually invalidate the use of a barotropic EOS for the mixture directly relating the averaged pressure and density. The use of correct EOS for the mixture determining the effective sound speed of the mixture is challenging given that, as shown in studies in the diluted limit, in most situations we cannot assume that the bubble interface is equilibrium with the surrounding and therefore the effective sound speed is a complex function. To revise the use of EOS imposes additional restrictions on the numerical methods that can be used: In those regions where mass transfer effects become important the sound speed of the mixture is significantly reduced and therefore it is likely to encounter supersonic flows [START_REF] Ganesh | Bubbly shock propagation as a mechanism for sheet-to-cloud transition of partial cavities[END_REF]. The use of numerical schemes able to correctly capture shock-waves seems then compulsory to correctly solve flows where phase change processes are important.

Other modeling difficulties discussed in this review include the development of reliable models for bubble nucleation, mass transfer processes and turbulence, the last two processes closely interacting with each other. The difficulties of averaged models to capture the experimental observations have been reported since long time to measure the turbulence intensity in cavitating flows [START_REF] Gopalan | Flow structure and modeling issues in the closure region of attached cavitation[END_REF] and also the density and stress distributions obtained in hydrodynamic applications [180]. Current lines of research for large scale bubble clusters are mainly focused on the use of large computational power able to carry on LES simulations of cavitating flows. Although these simulation have started to reproduce effects experimentally observed [155], the consistency and the range of applicability of the models and numerical methods used require revision in order to correctly account for the presence of multiple phases and their effect on large scale quantities (e.g. averaged void fraction, pressure distributions in surrounding objects, shear stresses...). To clarify the role of phenomena at small (unresolved) scales on the large scales in multiphase flows seems critical in order to develop advanced models able to accurately predict the response of bubble clusters in real applications.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Real part (continuous line) and imaginary part (dashed line) of the transfer function Φ as function of the vapor content Y 0 for an air/water-vapor bubble

Fig. 6 :

 6 Fig. 6: Numerical simulation of the development of large scale vortices generated by a bubble cloud raising in a supersaturated liquid.

  use a turbulence closure based on the kmodel with standard wall functions[113]

  Senocak Shyy [180] use a nonequilibrium version of the k-model reaching similar conclusions: while the predictions of pressure are not significantly influenced by the turbulence model other quantities, such as wall shear stress and velocity profiles, can be sensitive to it. Information about the application of this model for different test cases is given inShyy et al [182] and other variants of the kmodel include those fromJohansen et al [107], who propose a filtered based k-model.

Acknowledgements The author would like to acknowledge the exchanges and support of Professor Cesar Dopazo and the useful discussions with Maurice Rossi.

The authors declare that they have no conflict of interest.