Mahamadou Adama

Olivier Coutier-Delgosha

Daniel Buisine

M Adama Maiga

A new cavitation model based on bubble-bubble interactions

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Cavitation consists in the development of vapor bubbles in a liquid submitted to a local pressure drop. This phenomenon is usually considered as nearly isothermal in water at ambient temperature. Cavitation is inevitable in many industrial devices that involve high speed liquid flows, such as naval and submarine propellers, rocket engine pumps, water turbines, and hydraulic systems in general. It has several detrimental effects, such as noise, vibrations, and erosion on solid surfaces, the latter resulting from the local high-speed jets and intense pressure waves generated during the bubble collapses. The mechanisms of bubble implosion and subsequent erosion have been extensively studied in the last decades [see, for example, [START_REF] Philipp | Cavitation erosion by single laser-produced bubbles[END_REF] or [START_REF] Brujan | Collapse of micrometer-sized cavitation bubbles near a rigid boundary[END_REF]], and the damage due to cavitation was measured, for example, by [START_REF] Fortes-Patella | Cavitation damage measurement by 3D laser profilometry[END_REF] with three-dimensional laser profilometry. The study of cavitation is generally a challenging task, especially because of its problematic effects on high technology devices. On the scientific level, the influence of parameters such as interfacial tension, viscosity, quality of the liquid, surface properties [START_REF] Ezzatneshan | Study of surface wettability effect on cavitation inception by implementation of the lattice Boltzmann method[END_REF], and particle-bubble interactions [START_REF] Li | Transient interaction between a particle and an attached bubble with an application to cavitation in siltladen flow[END_REF] on cavitation a) Author to whom correspondence should be addressed: ocoutier@vt.edu.

Tel.: +1(540)2313168.

inception and bubble collapse is still an open question in hydrodynamics.

Investigation of cavitating flows is strongly connected to the development of specific experimental techniques, which enabled improving the understanding of this complex phenomenon in the last two decades [START_REF] Arakeri | Viscous effects on the position of cavitation separation from smooth bodies[END_REF][START_REF] Dowson | Cavitation in bearings[END_REF][START_REF] Lush | Visualisation of the cavitating flow in a venturi type duct using high-speed cine photography[END_REF][START_REF] Franc | Attached cavitation and the boundary layer: Experimental investigation and numerical treatment[END_REF][START_REF] Kubota | A new modelling of cavitating flows: A numerical study of unsteady cavitation on a hydrofoil section[END_REF]Tassin Leger and Ceccio, 1998;Tassin Leger et al., 1998;[START_REF] Astolfi | An experimental investigation of cavitation inception and development on a two-dimensional eppler hydrofoil[END_REF][START_REF] Arndt | Cavitation in vortical flows[END_REF][START_REF] Stutz | X-ray measurements within unsteady cavitation[END_REF][START_REF] Coutier-Delgosha | Analysis of cavitating flow structure by experimental and numerical investigations[END_REF][START_REF] Ceccio | Friction drag reduction of external flows with bubble and gas injection[END_REF]. It is now well-known that the behaviour of cavitating flows is greatly influenced by several factors, among which the gaseous micro-inclusions contained in the liquid, which are called cavitation nuclei. The dynamics of a population of these cavitation nuclei in a water tunnel has been investigated, for example, by [START_REF] Liu | Cavitation nuclei population dynamics in a water tunnel[END_REF].

When cavitation occurs, the microscopic air or vapour bubbles contained in the liquid progressively grow and become macroscopic bubbles mainly composed of vapour. However, at cavitation inception, [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF] shows that a multitude of bubbles appear and coalesce before growing. At this stage of incipient cavitation development, emissions of sound and light, which are the signatures of bubble collapses, are reported by [START_REF] Buogo | Implosion of an underwater spark-generated bubble and acoustic energy evaluation using the Rayleigh model[END_REF] and [START_REF] Takahashi | Experimental study on cavitation starting at and flow characteristics close to the point of separation[END_REF]. The very premature collapse of small bubbles is thus confirmed by these authors. These studies suggest that the development of cavitation does not consist in a linear growth of bubbles, but rather in a complex process of bubble interactions, which involves multiple production and collapse of microscopic bubbles. The process may be similar in visco-elastic materials that mimic human tissues, as reported, for example, recently by [START_REF] Oguri | Cavitation bubble nucleation induced by shockbubble interaction in a gelatin gel[END_REF].

Several types of cavitation models can be found in the literature. They can be roughly classified into two different groups: the two-fluid models and the single fluid models. The first category treats the liquid and vapour phases as two different fluids [START_REF] Grogger | Calculation of the cavitating flow in venturi geometries using two fluid model[END_REF][START_REF] Grogger | Calculation of the cavitating flow in venturi geometries using two fluid model[END_REF]Sauer and[START_REF] Sauer | Unsteady cavitating flow: A new cavitation model based on a modified front capturing method and bubble dynamics[END_REF]. They are seldom used for cavitating flows because the mass and momentum exchange terms between the two phases are usually not available. On the other hand, the single fluid models-also called homogeneous models-treat both phases as a single continuous, homogeneous fluid, with an average density. They are based on the classical balance equations for compressible or pseudo-compressible flows. They can be divided into three groups, according to the closure equation of the problem: (i) the barotropic models use a postulated equation of state that links the density to the pressure [START_REF] Delannoy | Two phase flow approach in unsteady cavitation modeling[END_REF][START_REF] Coutier-Delgosha | Numerical simulations in unsteady cavitating flows[END_REF][START_REF] Koop | Design of twisted cavitating hydrofoil using a barotropic flow method[END_REF], (ii) the models based on a transport equation for the volume fraction of vapour, which includes a source term for vaporization and condensation processes [START_REF] Ishii | Thermo-Fluid Dynamic Theory of Two-Phase Flow[END_REF][START_REF] Saito | Numerical analysis if unsteady vaporous cavitating flow around a hydrofoil[END_REF], (iii) the models that focus on bubble evolution, which are usually based on the Rayleigh-Plesset equation [START_REF] Plesset | Dynamics of cavitation bubbles[END_REF] or the Gilmore equation [START_REF] Gilmore | The growth or collapse of a spherical bubble in a viscous compressible liquid[END_REF]. Since they are based on the equations that govern the evolution of the cavitation nuclei, this third group of models may be the most appropriate to predict the phenomena involved in cavitation inception. Various expressions derived from the initial Rayleigh-Plesset formulation have been proposed in the last decades, to consider additional effects like viscosity, bubble/bubble interactions, or surface tension [see, for example, [START_REF] Mancas | Evolution of spherical cavitation bubbles: Parametric and closed-form Solutions[END_REF]].

Indeed, according to [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF], bubble-bubble interactions, including the collapse of microscopic bubbles at the early stage of cavitation development, should be considered for accurate modelling of cavitation inception. However, although such a phenomenon was clearly demonstrated by experimental studies, most of the current models do not include these effects. Typically, the equation proposed by [START_REF] Gilmore | The growth or collapse of a spherical bubble in a viscous compressible liquid[END_REF] provides the evolution of a unique bubble submitted to a pressure variation. Such a single-bubble model enables to predict correctly the evolution of one or several bubbles that survive during the first stage of cavitation expansion. However, it does not take into account the vapour volume transiently contained in all the bubbles that have disappeared during that time. Most of the studies devoted to the cloud of bubbles [START_REF] Chahine | Pressure field generated by the collective collapse of cavitation bubbles[END_REF][START_REF] D'agostino | Linearized dynamics of spherical bubbly clouds[END_REF][START_REF] Kubota | A new modelling of cavitating flows: A numerical study of unsteady cavitation on a hydrofoil section[END_REF] make the hypothesis of a constant number of bubbles. This assumption was also used by [START_REF] Kubota | A new modelling of cavitating flows: A numerical study of unsteady cavitation on a hydrofoil section[END_REF] who proposed an initial cloud of n bubbles which all undergo the same expansion. In an attempt to generalize this approach, [START_REF] Chen | Modeling hydrodynamics nonequilibrium in cavitating flows[END_REF] considered a variable number of bubbles with the hypothesis of no interfacial tension and without any exchange of volume between the bubbles. [START_REF] Oguz | A generalization of the impulse and virial theorems with an application to bubble oscillations[END_REF] also studied the hydrodynamic interactions of N bodies by impulse and virial theorems. They applied this approach to the mutual interaction of two bubbles during free and forced oscillations. Their results enable to explain the occurrence of stable bubble clusters initiated by acoustic cavitation. More recently, [START_REF] Du | A numerical model for the evolution of internal structure of cavitation cloud[END_REF] have developed a cavitation model that includes phenomena of bubble breakup inside a cloud of bubbles, which provided simulations in agreement with experiments for configurations of developed cloud cavitation.

In the present study, attention is particularly focused on the early stage of cavitation inception. A new model is proposed to take explicitly into account the interactions between two bubbles and their subsequent growth or collapse. Indeed, in case a single bubble in a control volume is submitted to a decreased pressure, which results in liquid vaporization at the interface, the bubble obviously receives the entire created vapour volume. Conversely, with at least two bubbles of different size in a control volume, the interaction between the bubbles, due to their different surface tensions, influences their evolution: it mitigates the expansion of the small one, while it increases the growth of the big one. This effect is called "volume exchange" hereafter in this paper.

The present model is an extension of the preliminary work from Adama Maiga et al. (2014). It enables to take into account the phenomenon of very premature disappearance of small bubbles and more generally includes possible interactions between nearby bubbles through the previously defined exchanges of volume. These ones induce a possible instability between neighbouring bubbles what leads the small bubbles to disappear for the benefit of the big bubbles. In the present approach, the biggest bubble inside a control volume V c is followed in its movement governed by the local flow velocity, taking into account the contributions of the small bubbles also located inside V c . This is a major originality of this model, compared with the approach reported by [START_REF] Plesset | Dynamics of cavitation bubbles[END_REF] and [START_REF] Gilmore | The growth or collapse of a spherical bubble in a viscous compressible liquid[END_REF], which focuses on a single bubble behaviour to predict its volume variation according to the pressure evolution. Using the divergence of the local velocity, the present model enables to predict the bubble evolution as well as the local pressure.

This paper is divided into three parts. In Sec. II, the model is presented in detail. Section III is dedicated to the model validation, using experimental measurements reported by [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]. In Sec. IV, the influence of the collapses of the smallest bubbles and volume exchanges on the bubble evolutions is analyzed, and the pressure predicted by the model is discussed.

II. MODEL DESCRIPTION

A. Theoretical approach

The possible scenario of the typical evolution of a cloud of vapour bubbles submitted to a pressure wave (tension followed by pressure re-increase) is represented schematically in Fig. 1. At the beginning, the cloud of bubbles expands rapidly. Some of the small bubbles disappear between times t 0 , which corresponds to cavitation inception, and t 1 . Indeed, surface tension generates pressure differences at the bubble interfaces, which induces a significant contraction of the small bubbles. At the same time, the volume of the remaining bubbles has increased. Between times t 1 and t 2 , the bubbles continue to grow without any new collapse.

Between t 2 and t 4 , the cloud is subjected to an opposite process of contraction. At the beginning, between times t 2 and t 3 , the compression of the small bubbles is accelerated, while some of the big bubbles still continue to grow. Between times t 3 and t 4 , all the bubbles almost disappear. From this schematic evolution, it can be postulated that for a given expansion there is a critical radius above which all bubbles will increase (when the volume is submitted to expansion) and below which any bubble will decrease (when the volume is submitted to contraction). This critical radius depends on various parameters, such as the internal pressure inside the bubbles, the surface tension, the number of bubbles, and the distance between them. It is thus a function of local conditions, rather than a unique value.

Ideally, a large number of bubbles should be considered. The purpose of this study is to build a simplified model that includes the mechanisms of bubble/bubble interactions by the exchange of volume. It has two objectives: (i) to determine the equation that governs the radius of the biggest local bubble, taking into account the interactions with the nearby bubbles, and (ii) to determine the local pressure. Hereafter the assumptions, the definition of the control volume, the liquid dynamics, the interface dynamics, and the equation for the volume conservation are detailed successively.

B. Assumptions

The following assumptions are made: (i) For the purpose of bubble modelling, the vorticity is assumed to be zero, so the liquid flow is supposed to be potential, and bubbles can be considered as sources. (ii) The bubble deformation due to the influence of the other bubbles is neglected. So, all bubbles are supposed to remain spherical. The flow around bubbles is supposed to be purely radial, and the evaporation or condensation is assumed to occur at the bubble interfaces, like in the equations of Rayleigh-Plesset or Gilmore. The slipping condition is thus not completely respected; the potential of the system is simply the sum of the separate potentials due to individual bubbles.

(iii) In the proposed model, the movement of bubbles with regard to the fluid is not considered, so the drag is zero.

In addition, it is assumed that inertia forces and gravity effects can be neglected. (iv) It is also supposed that thermal effects and phenomena of degassing play only a minor role in the prediction of big bubbles' evolution, which is the main objective of this work.

C. Control volume and sequences

As shown in Fig. 2, the study is performed inside a control volume V c , which contains two bubbles. Subscript 2 will be hereafter referred to the smallest bubble, while subscript 1 indicates the biggest one. The size of V c is determined by the bubble concentration of cavitation nuclei denoted η. The control volume V c is temporarily fixed, until one of the two initial bubbles collapses. Volume V c is submitted to a uniform average expansion resulting from the divergence of the average velocity field divU. The presence of two bubbles results in the additional possibility of volume exchange between them. Consequently, the small bubble may shrink and implode, or not, according to the conditions of expansion. If the small bubble disappears, a new control volume is defined and enlarged to contain a new nearby bubble. Then, a new sequence of calculation is started.

The value of the new control volume V n+1 c at sequence n + 1 is given by Eq. (2.1), which expresses the conservation of the volume fraction of vapour between sequences n and n + 1. In Eq. (2.1), V n c is the control volume, V n v is the vapour volume at the end of the previous sequence n, and V n+1 v is the new volume of vapour, i.e., the volume of the big bubble increased by one of the new small bubbles

V n v V n c = V n+1 v V n+1 c .
(2.1)

D. Dynamics of the liquid

Figure 3 schematises the respective positions of bubbles 1 and 2, with Ox being the axis that joins their centre and 2d being the distance between the centres. The potential of the system can thus be expressed as

F(x, y, z) = - 1 4π q 1 X F -X 1 + q 2 X F -X 2 , ( 2.2)
withX F = (x, y, z), X 1 = (-d, 0, 0), and X 2 = (d, 0, 0). It can be remarked that if the velocity of the bubbles due to their mutual interaction is neglected, the potential becomes identical to the one given by [START_REF] Oguz | A generalization of the impulse and virial theorems with an application to bubble oscillations[END_REF].

The velocity of the fluid on the line between the centres is

v(x) = 1 4π q 1 (d + x) 2 - q 2 (d -x) 2 , (2.3) with -d + R 1 < x < d -R 2 ,
where R 1 and R 2 are the radii of bubbles 1 and 2, respectively, and q 1 , q 2 are their expansion flow rate, respectively. The Navier-Stokes equation in the x direction can be written as

∂v ∂t + v ∂v ∂x = - 1 ρ L ∂p ∂x , (2.4)
where p is the local static pressure. Note that the viscosity term in the Navier-Stokes equations vanishes when the mass and momentum conservation is applied in the liquid around the bubble, as the conservation mass requires that the velocity field has an expression of the general form v(x) = F(t)/x 2 [for example, see [START_REF] Brennen | Cavitation and Bubble Dynamics[END_REF]]. By integrating Eq. (2.4) and using the velocity at the interface of bubble i:

v i = Ṙi = q i 4πR 2 i , with Ṙi = ∂R i
∂t , the following expression is obtained:

      p ρ L + 1 32π 2 q 1 (d + x) 2 - q 2 (d -x) 2 2      x 2 x 1 = 1 4π q1 (d + x) + q2 (d -x) x 2 x 1 .
(2.5)

Note that v i = Ṙi is an exact condition at the interface of bubble i in the case of zero mass transport across the interface. In our model, phase changes (vaporization and condensation) are supposed to occur at the bubble interfaces, but this expression of the velocity is still a good approximation [START_REF] Brennen | Cavitation and Bubble Dynamics[END_REF].

Equation (2.5) gives eventually the difference of pressure in the liquid at the bubble interfaces

p 1 -p 2 ρ L = q1 4πR 1 - q1 4π(D 12 -R 2 ) - q2 4πR 2 + q2 4π(D 12 -R 1 ) - q 2 1 32π 2 R 4 1 + q 2 2 32π 2 R 4 2 + q 2 1 32π 2 (D 12 -R 2 ) 4 - q 2 2 32π 2 (D 12 -R 1 ) 4 - q 1 q 2 16π 2 R 2 1 R 2 2 × R 2 1 (D 12 -R 2 ) 2 - R 2 2 (D 12 -R 1 ) 2 , (2.6)
with ρ L being the density of liquid,

q i = 4π • R 2 i Ṙi , and qi = ∂q i ∂t = 4π 2R i Ṙ2 i + R 2 i Ri .
(2.7)

E. Dynamics of the interface

The equation of equilibrium at the interface of a bubble can be expressed as

p i = p v + p g i -4µ L Ṙi R i -2 σ R i , (2.8a)
with i = 1, 2, Ṙi being the interface velocity of bubble i, σ being the surface tension, µ L being the dynamic viscosity, and p i being the pressure in the liquid at the interface of the bubbles.

The pressure inside the bubbles equals the sum of the vapour pressure p v and the gas pressure p g i . The latter is modelled by the barotropic law of Laplace

p g i = p 0 + 2σ R i0 -p v R i0 R i 3γ , (2.8b)
where p 0 is the reference pressure. So, the difference of pressure in the liquid at the interfaces of the bubbles is given by Eq. (2.9). In the calculations hereafter, it will be supposed that the partial pressure of the gas in the bubbles is initially zero

p 1 -p 2 = 2σ 1 R 2 - 1 R 1 + 4µ L Ṙ2 R 2 - Ṙ1 R 1 + p g 1 -p g 2 .
(2.9)

F. Volume conservation

In addition to the volume variations due to the exchanges between the two bubbles, each of them is submitted to the expansion of the control volume V c , which is related to the average velocity divergence divU of the fluid contained in V c . Hereafter, q m denotes the global flow rate due to this expansion, q i is the expansion/contraction flow rate of bubble i, and q e is the exchange flow rate between the two bubbles. The following relations can be written as:

q 1 + q 2 = q m ,
(2.10) (2.11)

q 1 -q 2 = q e ,
q m = S c U • n dS = V c divU dV = V c • divU, (2.12)
where U is the velocity vector. So, the volume fraction of vapour α(t)created during the sequence equals

α(t) = q m dt V c = divUdt. (2.13)
The conservation of the total volume of bubbles gives

R 3 1 + R 3 2 = R 3 10 + R 3 20 + 3 4π t t 0 q m dt, (2.14)
where R 10 and R 20 are the initial radii of the bubbles and t 0 is the time at the beginning of the sequence.

G. Equation for the radius R 1

From Eqs. (2.6), (2.7), (2.9), and (2.10), the equation governing the radius R 1 is obtained

q1 4πR 1 1 + R 1 R 2 - R 1 D 12 -R 1 - R 1 D 12 -R 2 - q 2 1 32π 2 R 4 1 1 - R 4 1 R 4 2 + R 4 1 (D 12 -R 1 ) 4 - R 4 1 (D 12 -R 2 ) 4 + 2R 2 1 (D 12 -R 1 ) 2 - 2R 4 1 R 2 2 (D 12 -R 2 ) 2 + µ L q 1 π ρ L R 3 1 1 + R 3 1 R 3 2 - q 1 q m 16π 2 R 4 2 1 + R 2 2 (D 12 -R 2 ) 2 - R 4 2 (D 12 -R 1 ) 4 - R 4 2 R 2 1 (D 12 -R 1 ) 2 = qm 4πR 2 1 - R 2 D 12 -R 1 - q 2 m 32π 2 R 4 2 1 - R 4 2 (D 12 -R 1 ) 4 + µ L q m π ρ L R 3 2 - 2σ ρ L R 1 1 - R 1 R 2 + 1 ρ L p g 1 -p g 2 , (2.15a)
while the equation for the radius R 2 is as follows:

q2 4πR 2 1 + R 2 R 1 - R 2 D 12 -R 2 - R 2 D 12 -R 1 - q 2 2 32π 2 R 4 2 1 - R 4 2 R 4 1 + R 4 2 (D 12 -R 2 ) 4 - R 4 2 (D 12 -R 1 ) 4 + 2R 2 2 (D 12 -R 2 ) 2 - 2R 4 2 R 2 1 (D 12 -R 1 ) 2 + µ L q 2 π ρ L R 3 2 1 + R 3 2 R 3 1 - q 2 q m 16π 2 R 4 1 1 + R 2 1 (D 12 -R 1 ) 2 - R 4 1 (D 12 -R 2 ) 4 - R 4 1 R 2 2 (D 12 -R 2 ) 2 = qm 4πR 1 1 - R 1 D 12 -R 2 - q 2 m 32π 2 R 4 1 1 - R 4 1 (D 12 -R 2 ) 4 + µ L q m π ρ L R 3 1 - 2σ ρ L R 2 1 - R 2 R 1 + 1 ρ L p g 2 -p g 1 . (2.15b)
Only one of the two equations (2.15a) or (2.15b) is eventually used and coupled with Eq. (2.14) for the numerical resolution. It will be denoted Eq. (2.15) hereafter.

If divU is known, Eqs. (2.14) and (2.15) lead to the determination of the evolution of radii R 1 and R 2 . If all terms that depend on D 12 are neglected in the left-hand term of Eq. (2.15), the remaining term is similar to the one reported by [START_REF] Plesset | Dynamics of cavitation bubbles[END_REF]. However, the right-hand side term depends only on q m , i.e., on the expansion rate, whereas it depends on the pressure given in [START_REF] Plesset | Dynamics of cavitation bubbles[END_REF].

Equations (2.14) and (2.15) show that the evolutions of bubbles 1 and 2 are strongly connected, whereas in the previous model proposed by [START_REF] Oguz | A generalization of the impulse and virial theorems with an application to bubble oscillations[END_REF], each bubble evolution is governed by a linearized Rayleigh-Plesset equation where the pressure field is modified to take into account the influence of the other bubble. In addition, the present expression of the bubble evolution as a function of the velocity divergence ensures a better conservation of the total volume of vapor created in the control volume V c , which is mandatory to model the volume exchange between the two bubbles.

In the particular case divU = 0, the present model enables to reproduce the condensation of the small bubbles for the benefit of the big ones, whose number progressively decreases until only one remains. This well-known phenomenon is called Ostwald ripening. Conversely, the cavitation models based on the evolution of a single bubble lead in such a case to stable isolated bubbles.

H. Equation for the pressure

The local pressure p at the point M depends on the distance between M and the bubbles. So, the mean pressure should be obtained by a technique of homogenization. In the present study, it is supposed that this averaged pressure is the pressure p c at equal distance between the two bubbles. With this hypothesis and since the distance 2d between bubbles is much larger than the bubble radius, by taking into account the dynamics in the liquid and at one bubble interface, two equations are available to determine the pressure p c , one for each bubble. The one associated with the big bubble is used here because the objective is to characterize the flow according to the visible two-phase structure, i.e., the big bubbles. The following equation is thus obtained:

p c ρ L = q1 4πR 1 R 1 D 12 -R 1 -1 + q 2 1 32π 2 R 4 1 × 1 - 64R 4 1 D 4 12 + R 4 1 (D 12 -R 1 ) 4 + 2R 2 1 (D 12 -R 1 ) 2 + q 1 q m 16π 2 R 4 1 32R 4 1 D 4 12 - R 4 1 (D 12 -R 1 ) 4 - R 2 1 (D 12 -R 1 ) 2 - q 1 π ρ L R 3 1 + qm 4πD 12 2 - D 12 D 12 -R 1 - 2σ ρ L R 1 + p g 1 ρ L .
(2.16)

Apart from the terms that depend on D 12 , this equation is the same as the Rayleigh-Plesset one, rewritten as a function of q 1 expansion flow rate of the big bubble. In the present model, Eq. (2.16) is used to calculate the pressure.

I. Initial conditions

Equation (2.14) is a second order equation that requires two initial conditions. In the initial state, the radius R 10 is given, and R 20 = kR 10 where the parameter k is given.

The two initial velocities Ṙi0 are calculated using Eq. (2.17). Ṙi0 is the sum of two contributions: one is derived from linearization of Eq. (2.15) without the terms q m and the second one is due to the expansion associated with divU 0 the velocity divergence at time t = 0

Ṙi0 = ω 1 ω 2 ω 2 -ω 1 R g -R i0 exp(ω 2 t) + V ci • divU 0 4πR 2 10
, (2.17)

with

ω 1 = -2 Re - 4 Re 2 + 2 We , ω 2 = -2 Re + 4 Re 2 + 2 We , We = ρ L R 3 g σ •T 2 c the Weber number, Re = ρ L R 2 g µ L T c
the Reynolds number, T c a characteristic time, R g the average initial radius of the bubbles, divU 0 the initial velocity divergence, R 10 the initial radius of the big bubble, and V ci the control volume containing the bubble i. The control volume containing the big bubble is

V c1 = 1 η , (2.18)
where η is the concentration of cavitation nuclei. In the calculations performed in Sec. III, T c = 1 µs and R g = 1 µm. The control volume containing the small bubble is

V c2 = R 3 20 R 3 10 • V c1 = k 3 0 η . (2.19)
So, the initial control volume is

V 0 c = 1 η 1 + k 3 0 .
(2.20)

Assuming here for the sake of simplicity a cubic control volume where each bubble is located at the centre of each half control volume, the distance between the bubble centres is

D 12 = 1 2 3 V 0 c . (2.21)
Depending on the conditions of expansion, the small bubble may disappear or not. A new sequence is started only after the small bubble has disappeared. It is considered that it disappeared if its radius reaches one tenth of its initial value, which corresponds to a volume 1000 times smaller than the initial one. To start a new step of the calculation, a new control volume is defined and enlarged to contain a new nearby bubble, which requires transition conditions.

J. Transition between sequences n and n + 1

At the transition between two sequences, the radius R 1 is constant. The radius R 2 is initialized as a function of the previous radius R 1 , according to the following condition:

R n+1 20 = kR n 1 , (2.22)
where the parameter k defines the radius of the smallest bubble that interacts with the big one: effects of smaller bubbles are neglected. For example, if k = 0.5 in the calculations, only the bubbles whose volume is smaller than 12.5% of the big bubble volume are neglected.

The initial velocity Ṙn+1 1 of the big bubble must also be updated at the beginning of the new sequence, since the control volume has changed. The best condition, which preserves the respective roles played by the small and the big bubbles, has been found to be the following equation:

Ṙn+1 1 V n+1 c 2γ exp /3 = Ṙn 1 V n c 2γ exp /3 , (2.23)
with γ exp = -1 in the expansion phase and γ exp = 1 during the contraction phase. So, from one sequence to the next, the radius R 1 is unchanged, the radius R n+1 20 is calculated using Eq. (2.22), V n+1 c is calculated using Eq. (2.1), and the velocity Ṙn+1 1 is calculated using Eq. (2.23).

III. MODEL VALIDATION

To validate the model, comparisons with the experimental results are performed hereafter. For this purpose, the time evolution of the velocity divergence divU must be introduced into the model, which subsequently predicts (i) the size evolution of the bubbles and (ii) the pressure variation. So, it requires the experimental determination of the velocity field or volume fraction of vapour (to obtain divU), the bubble size (to be compared to the evolution of the big bubble in the model), and the bubble collapses (to be compared to the behavior of the small bubble in the model) in a cavitating flow. Complexity of such flows makes very hazardous the measurements of these different unknowns: indeed, opacity of cavitation areas induces significant difficulties in flow visualization, and light reflection on bubble interfaces is a major drawback in the application of PIV (Particle Image Velocimetry) for velocity measurements. In addition, estimations of void fractions have been based until recently on observations, which did not enable to obtain such data with reasonable precision.

The validation of the model is performed hereafter based on the experimental study reported by [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]. To the best of our knowledge, this experiment is the only one in the literature where the radius of the cavitation bubbles and the pressure are both measured as a function of time at cavitation inception. The experiment is realized in a water basin made of stainless steel with glass windows on its four sides (Fig. 4).

The shock wave which is at the origin of the pressure drop is generated with an electromagnetic shock-wave generator at the bottom of the container. Two configurations are considered: (i) free passage of the shock wave in the negative x direction and (ii) reflection of the shock wave from a rigid reflector located at the top (x = 0), inducing a combination of wave propagations in the negative and positive directions. The pressure is recorded with a fiber optic hydrophone, and the resulting pressure evolution is modeled by the following multiparameter curve fitted in a least-square sense to the sampled data reported in [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]:

p ± = 1 2       p 10 exp(-α 1 (t ± x/c)) + p 20 exp(-α 2 (t ± x/c)) * cos       3 i=0 ν i ((t ± x/c)) i 2π(t ± x/c) + φ             * [1 + tanh(n 1 (t ± x/c))] * [1 + tanh(n 2 (∆ -(t ± x/c)))]. (3.1)
With c being the speed of sound, other constants are given in Table I.

The liquid used for the experiment is degassed water characterized by an oxygen content of 3 mg/l. The initial concentration of nuclei is not reported. With some simplifying assumptions, the liquid velocity in the x direction is modelled with a plane wave assumption

u l = - 1 ρ l c p -(t, x) + 1 ρ l c p + (t, x), (3.2)
where p -and p + denote the pressure waves travelling from negative x and from positive x, respectively. The velocity divergence divU can be derived from this liquid velocity

divU = ∂u l ∂x = - 1 ρ l c ∂p -(t, x) ∂x + 1 ρ l c ∂p + (t, x) ∂x . (3.3)
The final expressions of divU and its time derivative are given in the Appendix.

A. Case of a single pressure wave

In the first part of Ohl's experiment, the study is conducted without the reflector. Thus, in the expressions of the velocity and velocity divergence, only the pressure wave moving towards the negative x will be considered. Figure 5 shows the time evolutions of (i) the non-dimensional pressure measured in the experiment and (ii) the non-dimensional velocity divergence divU based on Eq. (3.3). The reference values are the maximum pressure and divergence values.

Figure 5 exhibits two successive steps: (i) the initial pressure wave results in a large pressure peak, which is correlated with a negative velocity divergence. This overpressure results in the collapse or fragmentation of the big nuclei. This first period where divU < 0 is thus characterized by a contraction of the vapour volume (condensation). (ii) The second step where divU > 0 (after t = 4 µs) corresponds to an expansion of the vapour volume (vaporization).

Based on a theoretical analysis, Ohl shows that only the germs smaller than 1.5 µm can survive to the passage of the initial pressure wave. For comparison of the model with the experimental measurements, a bubble of initial radius R 0 = 1.25 µm is thus considered hereafter. The influence of the first contraction phase (divU < 0 for t < 4 µs) on these germs can be neglected, so divU can be considered as equal to zero in this first phase.

qm and q m are obtained using Eq. (2.12), using the experimental velocity divergence and its partial time derivative. These two quantities allow determining the evolutions of the bubble radii [Eq. (2.15)] and comparing them with the experimental observations and measurements. Figure 6 displays the time evolution of the bubble radii R 1 and R 2 predicted by the model for a nuclei concentration of 15 000 germs/cm 3 (this value will be discussed below) and two different values for k: 0.5 and 0.8, respectively.

The results exhibit successive implosions of all small bubbles in the vicinity of the big one (42 and 15 collapses in the first 4 µs for k = 0.5 and 0.8, respectively). It reveals that even bubbles slightly smaller than the big one cannot develop under the influence of the big bubble. This result is consistent with Ohl's experimental observations, which is a major improvement compared with previous models such as Gilmore or modified Rayleigh-Plesset, which were found unable to reproduce the experimental observations [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]and Adama Maiga et al., 2014). Indeed, in the first version of the present model (Adama Maiga et al., 2014), some interactions were neglected, and although some bubble collapses were predicted at the first moments of the experiment, a small bubble eventually managed to develop close to the big one. So, the present result confirms that all interactions between bubbles should be considered in order to model correctly the first stage of cavitation development. The different expansions of the big bubble obtained for k = 0.5 and k = 0.8 are due to the number of collapses of small bubbles in the two calculations, before the expansion. Indeed, Adama Maiga et al. (2014) have shown in a previous study that the implosion of the small bubbles in the first stage of cavitation tends to favor the development of the big local bubbles.

Figure 7 displays the time evolutions of the big bubble radius R 1 predicted by the model for different values of the initial concentration of nuclei. The predicted radii are compared to the values measured in the experiments. The evolutions given by the Gilmore and Rayleigh-Plesset equations are also drawn on the same graph.

The results provided by the new model show that the initial concentration of nuclei significantly influences the bubble evolution. The lower the initial concentration is, i.e., the fewer the germs are present in the liquid, the larger the big bubble becomes. Indeed, previous studies (Adama Maiga et al., 2014;and Adama Maiga et al., 2015) have already shown that germs enhance the interactions between bubbles, which eventually slow down the bubble development. The results show that the best agreement between the experimental and numerical radii is obtained for an initial concentration η 0 = 15 000 nuclei/cm 3 . This value is lower than the one considered previously by [START_REF] Yuan | Modeling and computation of unsteady cavitation flows in injection nozzles[END_REF], while it is also slightly higher than those measured in [START_REF] Maeda | Laser holography measurement of bubble population in cavitation cloud on a foil section[END_REF] or [START_REF] Liu | Cavitation nuclei population dynamics in a water tunnel[END_REF]. However, it should be noted that in these studies, the germs smaller than about ten microns could not be detected. In addition, in the present configuration, the initial pressure peak at the beginning of the experiments tends to fragment the biggest germs and thus to increase the number of nuclei. The value η 0 = 15 000 nuclei/cm 3 will be applied in all the calculations presented hereafter in the paper.

In this case, the new model predicts very well the explosive growth phase of the bubble: the start of the expansion, the radius evolution, and the maximum value exhibit a nice agreement with the experimental observations. Conversely, the radius derived from the Gilmore and Rayleigh-Plesset equations are in poor agreement with the measured one. Both models predict a very early growth of the bubble, and also a very early decrease in the radius for t > 6 µs. It confirms that the models based on the evolution of a single bubble are not appropriate to reproduce the early stages of cavitation inception.

In the experimental study, Ohl observed that the cavitation bubbles become visible at t ≈ 4 µs, but without really explaining why they appear at that time. In the present model, the increase in the vapor volume requires a positive divU, which happens precisely in the experiments at t ≈ 4 µs (see Fig. 5). The good agreement between the experimental and simulated radius evolutions suggests that determining the time evolution of the bubble as a function of divU may be the best option for a correct prediction of cavitation inception. It can also explain why the previous models, which relate the bubble evolution to the local pressure, are not able to predict successfully the growth of the bubbles at the right time.

B. Bubbles submitted to combined positive and negative pressure waves

In the second part of the Ohl's experiment, a reflector was placed at the top of the setup, so the liquid in the container is submitted to both the incoming wave (towards negative x) and the reflected wave (towards positive x) with a small time delay between the two contributions. Time t = 0 is the time where the first wave reaches the reflector. In his study, Ohl has focused the analysis on an horizontal thin band located at about 1.5 mm below the reflector, where the the bubble size evolutions are recorded at 1 × 10 6 frames/s. The author observes some large differences in the maximum radius of bubbles located as close as 0.3 mm to each other and also some premature collapses of the smallest bubbles. Both phenomena cannot be reproduced with the Gilmore equation [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]. In the present modelling, it is assumed that the reflector is ideal, i.e., the pressure wave is entirely reflected. In addition, the pressure waves resulting from the expansions and collapses of the neighboring bubbles are neglected. Thus, as in the first part of the study, the flow velocity divergence is given by Eq. (3.3) with both terms included (Fig. 8). In his experiments, Ohl reported the evolution of two specific bubbles. Their position is 1.267 mm and 1.573 mm below the reflector, respectively, and they will be denoted bubble a and bubble b hereafter. For the modelling of these two bubbles, two options can be considered as follows:

(1) As the positions of the bubbles are different, it can be assumed that the bubbles a and b will not interact with each other. In this case (denoted case #1 hereafter), they are considered as two separate bubbles, each of them having its own control volume with a small neighboring bubble. Their evolutions are thus disconnected and should be very similar, since only their position is slightly different. Bubble 1 is considered at the position of bubble a and bubble 2 at the position of bubble b. (2) In case #2, both bubbles are located in the same control volume, so the two bubbles see the same pressure wave at the same time. In the model, bubble 1 will be bubble a and bubble 2 will be bubble b, as it was observed in the experiments that bubble b is smaller than bubble a, from the first instant it is detected to the maximum expansion of both bubbles.

Figure 9 shows the evolution of the two bubbles in case #1, for several values of the initial radius R 10 of bubble 1, from 1.25 µm (value applied previously) up to 5 µm. The experimental data reported by Ohl are also indicated.

The start of the expansions was detected at t = 2.8 µs and t = 3.8 µs for bubbles a and b, respectively, but the framerate was 1 image per microsecond, so the uncertainty on these values is significant. Here the model predicts in all simulations a start of the two bubble expansions around t = 3 µs, which is thus consistent with the observations. The slight delay between the two bubbles is due to their difference of location (bubble a being about 0.3 mm closer to the reflector, so it sees the reflected wave earlier than bubble b). Compared with the bubble of Sec. III A (without the reflector), which was expanding at 4 µs after the passage of the shock wave, the growth of the bubbles happens significantly earlier here because of the combined incoming and reflected waves.

For R 10 = 1.25 µm, the maximum radius of bubble a is bigger than those for R 0 = 3 µm and 5 µm and significantly overestimated, compared with the experiments. Similar to the case studied in Sec. III A, this difference in the bubble expansion is due to the number of collapses of small bubbles before the expansion phase: 36, 18, and 12 collapses for R 10 = 1.25, 3, and 5 µm, respectively. The best agreement with the experiments for the radius evolution of bubble 1 is obtained for R 10 = 5 µm. The presence of germs of this size, even though theoretically the initial pressure peak tends to make them collapse or fragment, can be explained by the first passage of the incident low pressure wave, after the pressure peak, which can promote the development of germs larger than 1.5 µm.

Conversely, the model does not reproduce the smaller size of bubble b reported by [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]. The radius of the second bubble is thus significantly overestimated. In addition, the behavior of both bubbles after t = 6 µs is not consistent with the experiments: in the model, both bubbles shrink significantly before stabilizing, while in the experiments, all bubbles at the location of bubble a continue to expand until t = 13 µs, and bubbles located close to bubble b collapse and rebound.

Figure 10 shows the results in case #2, where a single calculation is performed with bubbles a and b in the same control volume. Now a close agreement with the experimental data is obtained for both bubbles with the initial germ size R 10 = 5 µm: the maximum size of bubble 2 is about half the radius of bubble 1, similar to bubbles b and a in the experiments. Bubble 2 collapses at t ≈ 6 µs, which is also consistent with the experimental observations at the location of bubble b. It confirms that the mutual interactions between bubbles located in the thin band studied by Ohl are of primary importance in the dynamics of the bubbles and, more particularly, the disparities in bubble expansion and the collapse of multiple small bubbles in the vicinity of bigger bubbles that continue to grow.

This validation shows that the new model can be coupled with a computational fluid dynamics (CFD) flow calculation in the future: in each cell of the domain, the local velocity divergence would enable to determine the time evolution of the bubbles located inside, including the rate of collapses of the smallest bubbles.

IV. DISCUSSION

In this section, the analysis is focused on the effect of the size of the initial nuclei on the bubble evolutions, and the local pressure predicted by the model at mid-distance between the two bubbles. For this purpose, the model is applied to the configuration presented in Sec. III A, i.e., the single pressure wave leading to bubbles growing at t ≈ 4 µs in Ohl's experiment.

A. Influence of the nuclei size on the bubble growth

In Sec. III, it has been shown that small initial germs usually lead to larger expansions of the big bubble. This phenomenon is studied here in more details by considering three sizes of initial nuclei, namely, R 10 = 1 µm, 10 µm, and 100 µm. Figure 11 represents the time evolutions of the big bubble radius obtained for these values. The results basically confirm that small nuclei lead to larger bubble expansion rates: for R 10 = 1 µm, the maximum radius reached by bubble 1 is about 400 µm, while it is only 90 µm for R 10 = 10 µm. For R 10 = 100 µm, the bubble only slightly expands up to 120 µm.

To discuss these differences, it should be noted first that the initial void fraction varies significantly between the three cases, as the concentration of nuclei is kept constant η 0 = 15 000 nuclei/cm 3 . The resulting void fractions are α 0 = 6 × 10 -6 %, 6 × 10 -3 %, and 6%, respectively. In the two first cases, interactions between bubbles in neighbouring control volumes can certainly be neglected, while in the latter, such interactions may not be negligible, although they are not considered in the present modelling. However, the local interaction between the two bubbles is the main reason for the different results. Figure 12 shows the time evolution of the small bubble, for the three sizes of initial nucleus. The number of collapses is 47 for R 10 = 1 µm, 9 for R 10 = 10 µm, and only 2 for R 10 = 100 µm. It confirms the conclusion of Adama Maiga et al. (2014) with a simplified version of the present model, who showed that the interactions between bubbles and the collapses of small bubbles result in volume exchanges between bubbles, which favor the expansion of the local big bubbles.

To understand the mechanisms of bubble expansion, Fig. 13 represents the time evolutions of the flow rates q m (due to the expansion/contraction of the control volume, according to the sign of divU) and q e (due to the volume exchange between the two bubbles) defined in Sec. II. The highest value of q m is used to obtain non-dimensional values of all flow rates. Note that the q e curves for R 10 = 1 µm and R 10 = 10 µm are almost superimposed in Fig. 13,i.e., the evolutions of q e in the two cases look similar at this scale.

For t < 4 µs, the divergence is set to zero in the calculations: q m = 0 so q 1 = -q 2 , which means that any variation in the volume of bubble 1 is related to the same opposite variation for bubble 2, due to their mutual interaction. During this time, the volume of the big bubble increases, while the volume of the small bubbles decreases until their collapse (see Figs. 11 and 12). It results in a positive q e according to its definition [Eq. (2.11)]. For 4 µs < t < 6.2 µs, divU has become positive (see Fig. 5), which promotes the expansion of both bubbles. Indeed, bubble 2 grows explosively, like bubble 1, for R 10 = 1 µm and R 10 = 10 µm. Conversely, it continues to shrink (but more slowly) in the case of R 10 = 100 µm (see Fig. 12). During that time, q e is negative for the two smallest values of R 10 , which means that bubble 2 expands faster than bubble 1 (q 2 > q 1 ), while it remains positive (but smaller than previously) for R 10 = 100 µm, so the interaction between the two bubbles in this case is weakly affected by the evolution of divU. Finally, for t > 6.2, divU becomes negative after the pressure has re-increased, and a contraction phase is started, where q m is negative. It generally favors the compression of the bubbles, but bubble 1 almost stabilizes in the three cases (i.e., q 1 ≈ 0, see Fig. 11), while bubble 2 does shrink and even collapse (see Fig. 12), so q 2 is negative and q e becomes positive, as can be shown in Fig. 13. During that period, the big bubble volume remains nearly constant, although q m is negative, because of the interaction with the small bubble.

These general trends suggest that during the expansion phases, the variations in divU only slightly affect the growth of the big nuclei, which is mainly governed by their interaction FIG. 13. Evolution of the non-dimensional flow rates q m (due to expansion/contraction) and q e (due to the volume exchange between the two bubbles). with the neighboring bubbles, while it triggers the expansion of the small nuclei. Conversely, at the start of the compression phase, the behavior of the bubbles is not influenced anymore by the nucleus initial size: in all cases, the biggest bubbles will survive and keep almost the same size.

To get some more information about the influence of the nucleus size on the bubble/bubble interactions, the relative exchange flow rates q R e are shown in Fig. 14 q R e = q e V nucleus = dV e dt

V nucleus = dV R e dt , (4.1)
with V e being the "volume exchange" between bubble 1 and bubble 2 and V R e being the relative "volume exchange." Figure 14 enables to analyze the q e variations for the three initial sizes of nuclei, although the order of magnitudes of the flow rates is strongly different, due to the large differences between the nuclei volumes (V nucleus R10 = 100 µm /V nucleus R10 = 1 µm = 10 6 ).

Figure 14 confirms that for t < 4 µs, the q e evolutions depend strongly on the nucleus size: q e generally decreases for R 10 = 1 µm, it is almost stable for R 10 = 10 µm (not speaking about the multiple collapses), and it continuously increases for R 10 = 100 µm. These differences can also be seen in Fig. 11: as q e = 2q 1 (since q m = 0), the q e evolution is correlated with the slope of the big bubble radius evolution. It suggests that when divU = 0, nuclei of all sizes can expand, if they interact with smaller nearby nuclei. In this case, their expansion rate depends on their size: the bigger they are (or they become), the slower they grow.

For 4 µs < t < 6.2 µs, we can see now the differences between the two smallest nuclei: at t = 4 µs, q e decreases very suddenly, and its negative peak is much more intense for R 10 = 1 µm than for R 10 = 10 µm. It means than q 2 becomes much higher in the latter case, as can be seen in Fig. 12, where the expansion of bubble 2 is the most explosive. So, it confirms that during the expansion phases (divU > 0), the bubble/bubble interactions promote the fast expansion of the smallest nuclei, while the biggest ones are only slightly influenced (see the case R 10 = 100 µm, where the q e evolution is not significantly altered after t = 4 µs).

B. Local pressure

In addition to the time evolution of the bubbles' diameter, the model also provides Eq. (2.16) to determine the pressure p c at equal distance between the two bubbles. The time evolution of p c is shown in Fig. 15, together with the pressure measured in the experiments [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF]. Four successive steps can be observed: (1) From t = 0 to 2 µs, the measured pressure is positive.

During this step, the model predicts several intense pressure peaks, which are the signatures of the collapses of the small bubbles. In each new sequence, the initial small bubble is proportional to the big one, which means that each new small bubble considered in the calculation is a little bit bigger. As a result, the intensity of the pressure peaks also increases during that period, so the collapse of the small bubble is more and more violent. These multiple collapses are consistent with the experimental observations, as no bubble growth was observed during that period. (2) For 2 µs < t < 4 µs, the measured pressure becomes negative. As described previously, all existing models based on the pressure evolution would predict an intense vaporization, while the present model still predicts a systematic collapse of all small bubbles and only a slight growth of the big one, since divU is still negative during that period of time. This is also in agreement with the observations reported by Ohl. Indeed, cavitation inception does not necessarily happen when the pressure decreases below the vapor pressure. A delay in cavitation inception is observed, for example, in Venturi type sections or hydrofoil suction side, where cavitation does not develop at the Venturi throat or foil leading edge, even if the local pressure is lower than the vapor pressure, but slightly downstream.

(3) The large pressure drop at t ≈ 4 µs is related to the dynamics of the big bubble. Indeed, the pressure Eq. (2.16) depends on the radius of the two bubbles: it is similar to the classical equations that give the evolution of the bubble as a function of pressure, but here it provides the pressure evolution. As the big bubble begins its explosive growth at t = 4 µs, the local pressure necessarily drops at the same time. Conversely, the ambient pressure measured in the experiments has become negative at t = 2 µs, much earlier. (4) During the last 4 µs of the simulation, p c increases until a new bubble collapses. In this phase, apart from the pressure peak due to this collapse, the pressure evolution is fairly close to the ambient pressure measured in the experiment.

Obviously, the local pressure evolution predicted by the model is significantly different from the ambient pressure measured far from the bubbles. In the present case, numerous pressure waves due to the bubble collapses are obtained, which is of primary importance for the modeling of the damage and erosion due to cavitation. After the start of the big bubble expansion, the pressure in the middle of the two bubbles follows the evolution of the ambient pressure. So, in the very first instants of cavitation inception as well as later during the cavitation development, large amplitude pressure variations resulting from the multiple collapses of small bubbles should be taken into account, in addition to the ambient pressure evolution, to calculate the instantaneous local pressure in the liquid and eventually evaluate the flow aggressiveness and the resulting erosion.

The results also confirm that water can undergo negative pressure without cavitation [START_REF] Briggs | Limiting negative pressure in water[END_REF][START_REF] Zheng | Liquids at large negative pressures: Water at the homogeneous nucleation limit[END_REF]. The pressure drop measured by Ohl is about 80 bars, and the local negative peak predicted by the model at the beginning of the expansion phase is even higher. When the initial concentration of nuclei is reduced in the simulations, the intensity of this peak is reduced. This tendency is in agreement with the well-known observation in the cavitation experiments-the purer the water is, the more intense the pressure drop must be to initiate cavitation.

V. CONCLUSION

In this study, a new model for cavitating flows is proposed. Its originality consists in taking into account the interactions between bubbles located close to each other, which was called "exchange of volume." More specifically, the model predicts the dynamic behaviour of two non-identical bubbles located in a control volume. Both bubbles are submitted to an expansion/contraction law derived from the local flow velocity divergence. Depending on their mutual interaction, the small bubble may shrink and thereby contribute to the expansion of the nearby big bubble, or grow and thus slow down the big bubble development. If the small bubble eventually collapses, the control volume is redefined to include a new one, and the next sequence of simulation is initiated. The evolution of the big bubble during the successive sequences provides the information of the local variation in the volume fraction of vapour, and the pressure evolution at mid-distance between the two bubbles is also predicted.

Validation of the proposed model was performed in the configuration of the experiment conducted by [START_REF] Ohl | Cavitation inception following shock wave passage[END_REF] for the analysis of cavitation inception. It focuses on the first microseconds of the expansion of one or several bubbles. Comparison of the time evolution of the bubble radius predicted by the model with the one measured by Ohl have shown a close agreement, while the Rayleigh-Plesset or Gilmore equations failed to reproduce some major features observed in the experiments, like the time of bubble explosive growth. A more complex test case was also studied, where nearby bubbles grow at different rates, before some of them collapse while the other survive. This behavior was also well reproduced by the new model, which demonstrates the capability of this approach to simulate the very first moments of vapour development.

The present simulations generally show that bubble/bubble interactions are a primary mechanism of cavitation inception. Indeed, the development of the bubbles depends not only on the local velocity divergence resulting from the flow conditions but also on the influence of the nearby bubbles. Even if a positive local divU does promote vaporization, it will mostly trigger the expansion of the smallest nuclei, while the big ones may be weakly affected, due to the preponderant influence of the smaller nearby bubbles. In the compression phase, i.e., when divU becomes negative, the bubble/bubble interactions remain crucial, leading to the collapse of the smallest bubbles, while the biggest one can survive and even sometimes continue to expand. Unlike the expansion phases, the behavior of the bubbles at this stage is not influenced anymore by the nucleus initial size.

The local pressure resulting from these interactions exhibits multiple peaks which are the signature of the collapses of the small bubbles. When no bubble implosion does occur, like at the beginning of the compression phases, the local pressure is consistent with the evolution of the ambient pressure. It suggests that the large amplitude pressure variations resulting form the local collapses of small bubbles should be considered, in addition to the ambient pressure evolution, to calculate the instantaneous local pressure in the liquid and eventually evaluate the flow aggressiveness and the subsequent erosion.

The initial concentration of nuclei in the water is taken into account by the model, which is not the case in many previous cavitation models. This is a noticeable advantage, since it was, for example, written in the final report of the 23th ITTC (International Towing Tank Conference) that the development of new predictive methods for cavitation inception taking into account the water quality is strongly recommended [START_REF] Mehmet | The specialist committee on water quality and cavitation, Final report and recommendations to the 23rd international towing tank conference[END_REF].

APPENDIX: DERIVATION OF THE VELOCITY DIVERGENCE

The pressure equation is

p ± = 1 2      

(A1)

The constants are given in Table I.

The liquid velocity u l , divU, and

∂ ∂t (divU) = ∂ 2 u l ∂t∂x = - 1 ρ l c ∂ 2 p -(t, x) ∂t∂x + 1 ρ l c ∂ 2 p + (t, x) ∂t∂x , ( A4
)
where p -and p + denote the pressure waves traveling from negative x or from positive x, respectively. Using the velocity divergence, its time partial derivative, and the control volume, we obtain qm and q m . To simplify the calculations of the velocity divergence and its time partial derivative, one poses X a 1 = p 10 exp(-α 1 (t + ax/c)), (A5)

X a 2 = p 20 exp(-α 2 (t + ax/c)), (A6)

X a 3 = cos       3 i=0 ν i ((t + ax/c)) i 2π(t + ax/c) + φ       , ( A7
)
X a 4 = 1 + tanh(n 1 (t + ax/c)), (A8)

X a 5 = 1 + tanh(n 2 (∆ -(t + ax/c))), ( A9
)
where a isor +. Then, it is simple to calculate ∂t∂x (where a stands foror + and i = 1,2, . . . 5), which are easy to calculate, p a , u l , divU, and ∂ ∂t (divU) can be determined.

Adama Maiga, M., Coutier-Delgosha, O., and Buisine, D., "Interactions by exchange of volume between two spherical bubbles," Int. J. Appl. Mech. 06(06), 14500 (2014). Adama Maiga, M., Coutier-Delgosha, O., and Buisine, D., "Cavitation in a hydraulic system: The influence of the distributor geometry on cavitation inception and study of the interactions between bubbles," Int. J. Engine Res. 17(5), 543 (2015).

FIG. 1 .

 1 FIG. 1. A possible scenario for the evolution of a cloud of bubbles.

FIG. 2 .

 2 FIG. 2. Definition of a new control volume at the end of one calculation sequence.

FIG. 3 .

 3 FIG. 3. Schematic position of the two bubbles.

FIG. 4 .

 4 FIG. 4. Experimental setup. Reproduced with permission from "Cavitation inception following shock wave passage," Phys. Fluids 14(10), 3512-3521 (2002). Copyright 2002 AIP Publishing LLC.

  FIG. 5. Evolution of the pressure and the velocity divergence in the OHL experiments.

FIG. 6 .

 6 FIG. 6. Evolution of the radii R 1 and R 2 .

FIG. 7 .

 7 FIG. 7. Evolution of the big bubble radius: comparison between the present model, the experiments, the Gilmore model, and the Rayleigh-Plesset model.

FIG. 9 .

 9 FIG. 9. Model prediction in case #1.

FIG. 11 .

 11 FIG. 11. Evolution of the big bubble radius for different initial radii.

FIG. 14 .

 14 FIG. 14. Evolution of the relative q e for 3 values of the initial nuclei.

FIG. 15 .

 15 FIG.15. Pressure evolution obtained (i) with the model for R 10 = 1 µm and the initial concentration η 0 = 15 000 nuclei/cm 3 and (ii) in the experiments.

p

  10 exp(-α 1 (t ± x/c)) + p 20 exp(-α 2 (t ± x/c)) * cos (t ± x/c)) i 2π(t ± x/c) + φ tanh(n 1 (t ± x/c))] * [1 + tanh(n 2 (∆ -(t ± x/c)))].

TABLE I .

 I Values of constants in Eq. (3.1).

	P 01 = 2.275 × 10 7 Pa	P 02 = 2.77 × 10 6 Pa	α 1 = 1.767 × 10 6 s 1


Miyoung Han

Pierre Senellart

Stéphane Bressan

Huayu Wu

Routing an Autonomous Taxi with Reinforcement Learning

Singapore's vision of a Smart Nation encompasses the development of effective and efficient means of transportation. The government's target is to leverage new technologies to create services for a demand-driven intelligent transportation model including personal vehicles, public transport, and taxis. Singapore's government is strongly encouraging and supporting research and development of technologies for autonomous vehicles in general and autonomous taxis in particular. The design and implementation of intelligent routing algorithms is one of the keys to the deployment of autonomous taxis. In this paper we demonstrate that a reinforcement learning algorithm of the Q-learning family, based on a customized exploration and exploitation strategy, is able to learn optimal actions for the routing autonomous taxis in a real scenario at the scale of the city of Singapore with pick-up and drop-off events for a fleet of one thousand taxis.

INTRODUCTION

Singapore's Prime Minister Office outlines its vision of a Smart Nation as the harnessing of information and communication technologies "to support better living, create more opportunities, and support stronger communities" 1 . The vision encompasses the development of effective and efficient means of transportation and aims at the development of a demand-driven intelligent transportation model including personal vehicles, public transports, and taxis. In particular, Singapore's government is strongly encouraging and supporting research and development of technologies for autonomous vehicles (AV). Its agencies, together with academia and industry, are involved in learning and understanding opportunities and challenges for autonomous vehicle technologies and autonomous taxis. For example, the Land Transport Authority of Singapore [START_REF] Tan Cheon | Autonomous vehicles, next stop: Singapore[END_REF] has been working with A*STAR's Institute for Infocomm Research to evaluate their 1 http://www.pmo.gov.sg/smartnation Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org. AV Proof-of-Concept on-road trials2 since last year. "Unlike autonomous vehicle trials elsewhere, Singapore is focusing on applying the technology to public buses, freight carriers, autonomous taxis and utility operations such as road sweepers." [START_REF] Tan | Driverless vehicles hit the road in trials around singapore[END_REF]. Self-driving riders have already been experienced by the public at Gardens by the Bay last December 3 .

The design and implementation of intelligent routing algorithms is one of the keys to the deployment of autonomous taxis. While prior studies have explored reinforcement learning approaches to taxi routing, they have mostly done so with synthetic models and data in small and static state spaces. In this paper we demonstrate that a reinforcement learning algorithm is able to progressively, adaptively, efficiently, and effectively learn optimal actions for the routing to passenger pick-up points of an autonomous taxi in a real scenario at the scale of the city of Singapore.

We present a solution based on a customized exploration and exploitation strategy for an algorithm of the Q-learning family [START_REF] Watkins | Technical note: Q-learning[END_REF]. Q-learning is an approach to solve reinforcement learning [START_REF] Sutton | Introduction to Reinforcement Learning[END_REF] problems. Reinforcement learning is used for learning an optimal policy in a dynamic environment: an agent takes an action in a state, receives a reward, moves to some next state, and repeats this procedure. Initially, the agent has no knowledge of which action has to be taken in a given state. The algorithm has the choice between exploiting its knowledge by choosing the action with highest estimated value and exploring its environment by taking any other action. The trade-off between exploration and exploitation is crucial. We devise an exploration strategy that follows an almost-greedy policy. The strategy first selects a set of candidate actions with long-term reward above a given threshold and then selects, among candidates, the action with the highest probability of finding a passenger.

We evaluate how Q-learning progressively learns and how to balance exploration and exploitation. We quantify the influence of the parameters. We propose an original selection strategy. The evaluation uses a dataset of taxi pickups and drop-offs for a fleet of 1000 taxis for one month in Singapore.

BACKGROUND AND RELATED WORK

Most studies addressing the taxi routing problem focus on providing the fastest route and a sequence of pick-up points [START_REF] Qu | A cost-effective recommender system for taxi drivers[END_REF] by mining historical data [START_REF] Qu | A cost-effective recommender system for taxi drivers[END_REF][START_REF] Yuan | T-drive: Enhancing driving directions with taxi drivers' intelligence[END_REF][START_REF] Yuan | T-finder: A recommender system for finding passengers and vacant taxis[END_REF][START_REF] Zheng | Trajectory data mining: An overview[END_REF]. Yuan et al. [START_REF] Yuan | T-drive: Enhancing driving directions with taxi drivers' intelligence[END_REF] cluster road segments and travel time to build a landmark graph of traffic patterns and time-dependent fastest routes. They then present in [START_REF] Yuan | T-finder: A recommender system for finding passengers and vacant taxis[END_REF] a recommendation system for taxi drivers and passengers by clustering road segments extracted from GPS trajectories. The system recommends a parking place and road segments to taxi drivers and passengers. Qu et al. [START_REF] Qu | A cost-effective recommender system for taxi drivers[END_REF] propose a method to recommend a route. They develop a graph representation of a road network by mining the historical taxi GPS traces and generate an optimal route for finding passengers. Those models rely on the availability of accurate historical data and trajectories. They might not be relevant in dynamic environments such as an autonomous taxi looking for optimal passenger pick-up points.

Reinforcement learning [START_REF] Sutton | Introduction to Reinforcement Learning[END_REF] has the potential to continuously and adaptively learn from interaction with the environment. The algorithm discovers which actions produce the greatest reward by experience and estimates how good it is to take a certain action in a given state. Yet reinforcement learning aims to maximize cumulated reward. The notion of Markov Decision Process (MDP ) underlies much of the work on reinforcement learning, and is at the basis of our work.

Q-learning [START_REF] Watkins | Technical note: Q-learning[END_REF] is widely used because of its computational simplicity. In Q-learning, one does not require a model of transition functions and reward functions but learns directly from observed experience. This is an advantage as learning a model often requires exhaustive exploration that is not suitable for a large state space. Thus, Q-learning seems suitable for taxi routing in a city of the scale of Singapore.

While taxi routing has often been used as the example application for reinforcement learning algorithms, it often remained relegated to toy or small scale examples, as it is the case of the seminal 5×5 grid introduced by the authors of [START_REF] Dietterich | The MAXQ method for hierarchical reinforcement learning[END_REF] and used for experimental purposes by the authors of [START_REF] Diuk | A hierarchical approach to efficient reinforcement learning in deterministic domains[END_REF][START_REF] Ghavamzadeh | A multiagent reinforcement learning algorithm by dynamically merging Markov decision processes[END_REF][START_REF] Ghavamzadeh | Learning to communicate and act using hierarchical reinforcement learning[END_REF][START_REF] Hester | Generalized model learning for reinforcement learning in factored domains[END_REF]. Learning pick-up points is a somewhat new application of reinforcement learning.

Q-LEARNING FOR TAXI ROUTING

We consider an autonomous taxi completely relying on reinforcement learning. The algorithm decides where the autonomous taxi should go in order to pick up passengers by learning the existence probability of passengers from its gathered experience. Reinforcement learning models an agent taking an action a in state s, receiving a reward r, and moving to the next state s . With Q-learning, the estimated value of taking action a in state s, denoted Q(s, a), is updated as:

Q(s, a) := Q(s, a) + α [r + γ max a Q(s , a ) -Q(s, a)] .
Here, α is a positive fraction such that 0 < α 1, the stepsize parameter that influences the rate of learning. When α = 1, the agent considers only the most recent information for learning. If α is properly reduced over time, the function converges [START_REF] Sutton | Introduction to Reinforcement Learning[END_REF]. The discount rate γ (0 γ < 1) determines the present value of future rewards. If γ = 0, the agent is only concerned with the immediate reward. The agent's action influences only the current reward. If γ approaches 1, the agent considers future rewards more strongly.

To maximize total reward, the agent must select the action with highest value (exploitation), but to discover such action it has to try actions not selected before (exploration). The ε-greedy method works as follows: most of the time, select an action with the highest estimated value, but with small probability ε select an action uniformly at random. This exploration enables to experience other actions not taken Algorithm 1 Taxi Routing for Learning Pick-up Points 1: Initialize Q(s, a), existence probability of passengers p 2: repeat 3: repeat 4:

if greedy then 5:

V := { a ∈ A | Q(s, a) max a Q(s, a ) -η } 6: if |V | > 1 then 7:
Select action a with highest probability p 8:

else /* not greedy */ 9:

Select action a uniformly at random 10:

Take action a, obtain reward r, observe next state s 11:

Q(s, a) := Q(s, a) + α [r + γ max a Q(s , a ) -Q(s, a)] 12:
Increment visit count on s 13:

Update existence passenger probability p(s ) 14:

if passenger found in s then 15:

Increment found count on s 16:

s becomes the end of the passenger route from s 17: else 18: s := s 19: until a passenger is found 20: until algorithm converges before and it may increase the greater total reward in the long run because we would discover better actions.

We call an episode the movement of a taxi to an actual pickup point. For the first episode, the taxi located at a random position moves according to the initial oblivious learning policy. The episode ends when the taxi finds a passenger. Then, it moves to the passenger's destination and starts a new episode. As the taxi moves it receives rewards and updates its action-value and the existence probability. The road network is discretized and the movements correspond to steps in the discretized network. At each step, the taxi learns where passengers are likely to be located.

The Taxi Routing algorithm for learning pick-up points is outlined in Algorithm 1. According to the ε-greedy policy, an action a is selected in a given state s.

The action selection rule selects the action with the maximum estimated action value (greedy action). However, with this rule, the algorithm ignores other actions that, although they have slightly lesser value, may lead to a state having a higher chance to pick a passenger up. Hence, instead of selecting one greedy action, we loosen the selection condition by setting a lower bound below the maximum value in order to choose from more potentially valuable candidate actions (Line 5). The candidate actions are compared with existence probabilities of passengers in their corresponding states (Line 7). We later refer to the algorithm with this selection strategy as Q-learning using LB/Proba.

After taking an action, we update the Q-value in the current state s with reward r and next state s . As we visit a new state s , the visit and found counts are incremented and the existence probability of passengers is also recalculated. We repeat this procedure until we find a passenger.

PERFORMANCE EVALUATION

For the sake of simplicity, in this paper, we present the results for a map discretized into cells of 0.01 degree longitude × 0.01 degree latitude (about 1.1km × 1.1km) forming a 38 × 20 grid. At each cell of the grid, eight actions are possible: up, down, right, left, and diagonally. A step is the movement from one cell to an adjacent one. Although such a representation does not capture several natural constraints on the traffic, it is sufficient, with limited loss of generality, to evaluate the effectiveness of the algorithm.

Since popular pick-up points generally depend on the time of the day, we run the experiments for selected time intervals. Here, we present the results for two off-peak hours (12h to 13h and 14h to 15h), but we obtain comparable results for other time slots. At each episode we select 300 passengers according to actual geographical distribution in the given time interval in a dataset of taxi pickups and drop-offs for a fleet of one thousand taxis for one month in Singapore.

We first look into the impact of the step-size parameter α, the discount rate γ, and the probability of exploration ε. We evaluate how these parameters influence the learning performance with ordinary Q-learning. We compare the average number of steps. The average steps are calculated at every 100 episode by dividing the total steps from the first episode to the last by the total number of episodes.

Figures 1a-1b show the average number of steps with different step-size parameter α values for different time intervals. We compare four different α with a fixed γ (= 0.5) and ε (= 0.3). For all the time intervals, as the α is smaller, the average number of steps also decreases. Lower step-size values perform better. This indicates that accumulated experience affects value estimation more significant than recent experience, i.e., that the problem is indeed stochastic.

For the discount rate γ experiment, we fixed α (= 0.5) and ε (= 0.3) and changed the γ. The average number of steps with different γ values for different time intervals are shown in Figures 1c-1d. In Figure 1d, the lowest γ (= 0.25) perform better. This means that immediate rewards are more important than future rewards. In Figure 1c, as episodes continue, a higher γ (= 0.75) is slightly better. In 1c, relatively longer step counts than those of the other time intervals are needed to achieve a goal. In this case, future rewards are more significant than current rewards.

Figures 1e-1f show the average number of steps with different ε values for different time intervals, given α = 0.5 and γ = 0.5. The average number of steps first decreases dramatically and then converges gradually. For all time intervals, when ε is 0.1, the average number of steps is higher than the other cases in early episodes but it dominates after about 30,000 episodes. At the beginning, exploration is more effective and relatively inexpensive. Eventually, sufficient knowledge is accumulated and exploitation is worthy.

We now compare Q-learning using LB/Proba (our algorithm) with ordinary Q-learning. For experiments, we select three parameter values shown in the previous section. The step-size parameter α is set to 0.25 because low learning rate is appropriate to our problem. The probability of exploration ε is set to 0.1. Since loosened selection for maximum action has the effect of exploration, high ε is not needed. We take η = 0.01 to set the lower bound on the maximum action value per state. Two algorithms are compared by varying the discount rate γ value.

In Figures 1g-1j, when γ = 0.75 or 0.5, Q-learning using LB/Proba converges faster than Q-learning. On the other hand, when γ is 0.25 (Figures 1k-1l), Q-learning performs similary well or slightly better than Q-learning using LB/Proba. These experiments show that when the learning rate α is low and the discount rate γ approaches 1, Q-learning using LB/Proba outperforms Q-learning. In other words, it has to accumulate much experience for value prediction and it considers future rewards more strongly. The reason is that Q-learning using LB/Proba depends on existence probability of passengers that requires enough experience and that is more related to long-term high rewards.

CONCLUSIONS AND FUTURE WORK

A*STAR's Institute for Infocomm Research, with its academic partners, is contributing to Singapore's Land Transport Authority's efforts to evaluate the viability of autonomous vehicles. We explored the use of reinforcement learning for autonomous taxi routing. The evaluation results provide empirical support of the effectiveness of reinforcement learning for the task at hand. We are planning more extensive studies as soon as extensive on-road trials are carried out and more data is available to us. We are also currently extending this work to consider a more dynamic environment. We are considering a multi-agent approach to scale to larger and more complex state spaces. We also investigate hierarchies, factored representations, and other extended methods that can improve adaptiveness, effectiveness and efficiency of learning.

Figure 1 :

 1 Figure 1: Average number of steps as the number of episodes increases.

http://www.lta.gov.sg/apps/news/page.aspx?c=2&id= 7dac9bd2-e6d1-448f-ba0e-f31fd13a8c3e

http://www.opengovasia.com/articles/6792-asias-firstfully-operational-autonomous-vehicle-running-at-gardensby-the-bay-singapore

Acknowledgment. This research was supported in part by the Télécom ParisTech Chair on Big Data & Market Insights, and in part by the National Research Foundation, Prime Minister's Office, Singapore under the CREATE program, NUS ref. R-702-005-101-281.