A new cavitation model based on bubble-bubble interactions

In this paper, a new model based on bubble-bubble interactions is proposed for cavitation. Unlike the well-known existing models (Rayleigh-Plesset, Gilmore), which are derived from the local balance equations in the vicinity of a single cavitation bubble, the present approach is based on the mutual interaction between two spherical bubbles of different sizes. The mass and momentum conservation equations, coupled with the local flow divergence, lead to two equations for the evolution of the bubble radii and one equation for the local pressure. The bubble size variations predicted by the model are found in close agreement with the previous experimental data reported by Ohl [“Cavitation inception following shock wave passage,” Phys. Fluids 14(10), 3512–3521 (2002)]. The distinct radii of bubbles located close to each other, as well as the premature collapse of small bubbles during the initial stage of cavitation inception, are correctly reproduced by the model. The results generally show that bubble/bubble interactions play a primary role in the physics of cavitation inception, which is a preponderant phenomenon in cavitation-induced noise and erosion. The influence of the size of the nuclei on these interactions is discussed. During the expansion phases, the variations in the local flow divergence only slightly affect the growth of the big nuclei, which is mainly governed by their interaction with the neighboring bubbles, while it triggers the expansion of the small nuclei. Conversely, in the compression phase, the behavior of the bubbles is not influenced anymore by the initial size of the nuclei. It is also shown that 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.In this paper, a new model based on bubble-bubble interactions is proposed for cavitation. Unlike the well-known existing models (Rayleigh-Plesset, Gilmore), which are derived from the local balance equations in the vicinity of a single cavitation bubble, the present approach is based on the mutual interaction between two spherical bubbles of different sizes. The mass and momentum conservation equations, coupled with the local flow divergence, lead to two equations for the evolution of the bubble radii and one equation for the local pressure. The bubble size variations predicted by the model are found in close agreement with the previous experimental data reported by Ohl [“Cavitation inception following shock wave passage,” Phys. Fluids 14(10), 3512–3521 (2002)]. The distinct radii of bubbles located close to each other, as well as the premature collapse of small bubbles during the initial stage of cavitation inception, are correctly reproduced by the model. The results generally show that bubble/bubble i...


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, Philipp and Lauterborn (1998) or Brujan and Matsumoto (2012)], and the damage due to cavitation was measured, for example, by Fortes-Patella et al. (2000) 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 (Ezzatneshan, 2017), and particle-bubble interactions (Li et al., 2018) on cavitation a) Author to whom correspondence should be addressed: ocoutier@vt.edu.
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 (Arakeri, 1975;Dowson and Taylor, 1979;Lush and Peters, 1982;Franc and Michel, 1985;Kubota et al., 1992;Tassin Leger and Ceccio, 1998;Tassin Leger et al., 1998;Astolfi et al., 2000;Arndt, 2002;Stutz and Legoupil, 2003;Coutier-Delgosha et al., 2007;and Ceccio, 2010).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 Liu et al. (1993).
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, Ohl (2002) 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 Buogo and Cannelli (2002) and Takahashi et al. (2003).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 Oguri and Ando (2018).
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 (Grogger andAlajbegovic, 1998 andSauer andSchnerr, 2000).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 (Delannoy and Kueny, 1990;Coutier-Delgosha et al., 2003;and Koop and Hoeijmakers, 2006), (ii) the models based on a transport equation for the volume fraction of vapour, which includes a source term for vaporization and condensation processes (Ishii, 1975 andSaito et al., 2003), (iii) the models that focus on bubble evolution, which are usually based on the Rayleigh-Plesset equation (Plesset, 1948) or the Gilmore equation (Gilmore, 1952).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, Mancas and Rosu (2016)].
Indeed, according to Ohl (2002), 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 Gilmore (1952) 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 (Chahine, 1982;D'Agostino and Brennen, 1989;and Kubota et al., 1992) make the hypothesis of a constant number of bubbles.This assumption was also used by Kubota et al. (1992) who proposed an initial cloud of n bubbles which all undergo the same expansion.In an attempt to generalize this approach, Chen and Heister (1996) considered a variable number of bubbles with the hypothesis of no interfacial tension and without any exchange of volume between the bubbles.Oguz and Prosperetti (1990) 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, Du et al. (2016) 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 Plesset (1948) and Gilmore (1952), 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 Ohl (2002).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.

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 (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 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 Oguz and Prosperetti (1990).
The velocity of the fluid on the line between the centres is 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 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 Brennen (1995)].By integrating Eq. (2.4) and using the velocity at the interface of bubble i: ∂t , the following expression is obtained: (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 (Brennen, 1995).
Equation (2.5) gives eventually the difference of pressure in the liquid at the bubble interfaces with ρ L being the density of liquid, (2.7)

E. Dynamics of the interface
The equation of equilibrium at the interface of a bubble can be expressed as 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 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 (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: (2.10) (2.11) where U is the velocity vector.So, the volume fraction of vapour α(t)created during the sequence equals The conservation of the total volume of bubbles gives 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 while the equation for the radius R 2 is as follows: 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 Plesset (1948).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 Plesset (1948).
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 Oguz and Prosperetti (1990), 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: (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 , (2.17) with 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 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 So, the initial control volume is (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 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: 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: 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 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 Ohl (2002).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 Ohl (2002): 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 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 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 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 (Ohl, 2002;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 Yuan et al. (2001), while it is also slightly higher than those measured in Maeda et al. (1991) or Liu et al. (1993).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 (Ohl, 2002).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 Ohl (2002).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 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 (Ohl, 2002).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 (Briggs, 1950;and Zheng et al., 1991).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 Ohl (2002) 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 (Mehmet, 2002).

APPENDIX: DERIVATION OF THE VELOCITY DIVERGENCE
The pressure equation is

(A1)
The constants are given in Table I.
The liquid velocity u l , divU, and ∂ ∂t (divU) are 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  ∂t∂x (where a stands for − or + and i = 1,2, . . .5), which are easy to calculate, p a , u l , divU, and ∂ ∂t (divU) can be determined.

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

FIG. 2 .
FIG. 2. Definition of a new control volume at the end of one calculation sequence.
FIG. 5. Evolution of the pressure and the velocity divergence in the OHL experiments.

FIG. 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. 11 .
FIG. 11.Evolution of the big bubble radius for different initial radii.

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

TABLE I .
Values of constants in Eq. (3.1).