D Fuster

F Montel

Mass transfer effects on linear wave propagation in diluted bubbly liquids

This article investigates the importance of mass transfer effects in the effective acoustic properties of diluted bubbly liquids. The classical theory for wave propagation in bubbly liquids for pure gas bubbles is extended to capture the influence of mass transfer on the effective phase speed and attenuation of the system. The vaporization flux is shown to be important for systems close to saturation conditions and low frequencies. We derive a general expression for the transfer function that relates bubble radius and pressure changes solving the linear version of the conservation equations inside, outside and at the bubble interface. Simplified expressions for various limiting situations are derived in order to get further insight about the validity of common assumptions typically applied in bubble dynamic models. The relevance of the vapor content, the mass transfer flux across the interface and the effect of the bubble interface temperature variations is discussed in terms of characteristic nondimensional numbers. Finally we derive the various conditions that must be satisfied in order to reach the low frequency limit solutions where the phase speed does no longer depend on the forcing frequency.

Introduction

The effect of bubbles on the process of wave propagation in liquids is an open scientific problem with many applications fields such as ultrasonic fluid flow monitoring and geophysics [START_REF] Kuster | Velocity and attenuation of seismic waves in twophase media: Part i. theoretical formulations[END_REF][START_REF] Lynnworth | Ultrasonic measurements for process control: theory, techniques, applications[END_REF]. One family of models typically used for wave propagation in bubbly flows applies for dilute systems. Most of these models [START_REF] Van Wijngaarden | On the equations of motion for mixtures of liquid and gas bubbles[END_REF][START_REF] Chapman | Thermal effects in the free oscillations of gas bubbles[END_REF][START_REF] Prosperetti | Thermal effects and damping mechanisms in the forced radial oscillations of gas bubbles in liquids[END_REF][START_REF] Sangani | A pairwise interaction theory for determining the linear acoustic properties of dilute bubbly liquids[END_REF][START_REF] Zhang | Momentum and energy equations for disperse two-phase flows and their closure for dilute suspensions[END_REF][START_REF] Ando | Improvement of acoustic theory of ultrasonic waves in dilute bubbly liquids[END_REF]) succeed reproducing experimental results for frequencies below resonance in situations where the amount of vapor is negligible [START_REF] Silberman | Sound velocity and attenuation in bubbly mixtures measured in standing wave tubes[END_REF][START_REF] Cheyne | Phase velocity measurements in bubbly liquids using a fiber optic laser interferometer[END_REF][START_REF] Wilson | Phase speed and attenuation in bubbly liquids inferred from impedance measurements near the individual bubble resonance frequency[END_REF][START_REF] Leroy | Sound velocity and attenuation in bubbly gels measured by transmission experiments[END_REF]). These theories have been recently improved at frequencies above the bubble natural frequency. For instance, high frequency corrections have been proposed by [START_REF] Ando | Improvement of acoustic theory of ultrasonic waves in dilute bubbly liquids[END_REF], while [START_REF] Fuster | Effect of direct bubble-bubble interactions on linear-wave propagation in bubbly liquids[END_REF] propose correction terms in order to capture direct bubble-bubble interactions based on the results obtained from a full non-linear model [START_REF] Fuster | Modeling bubble clusters in compressible liquids[END_REF].

Theoretical models for linear wave propagation in bubbly liquids require to model the response of a single bubble to an external pressure excitation. This response is obtained through the linearization of the Rayleigh-Plesset equation and the solution of the mass, momentum and energy conservation equations inside the bubble. Expressing the amplitude of the bubble radius oscillation as a function of the external pressure amplitude allows us to define a bubble resonance frequency and damping factor, see [START_REF] Ainslie | Review of scattering and extinction cross-sections, damping factors, and resonance frequencies of a spherical gas bubble[END_REF] for a review. These parameters depend on the so-called transfer function, which relates bubble pressure and volume changes in the linear regime. [START_REF] Prosperetti | Nonlinear bubble dynamics[END_REF] propose an expression capturing the heat transfer exchange between the bubble and its environment assuming that the bubble interface temperature is constant and neglecting mass transfer effects across the interface. This function can be used to define an effective state equation for the bubble interior of the type p b V γ eff b = C, where p b is the bubble pressure, V b is the bubble volume and γ eff is a frequency dependent effective coefficient that tends to recover the isothermal limit for low frequencies (γ eff = 1) and the adiabatic response for very large frequencies, γ eff = γ, where γ is the gas polytropic coefficient.

Wave propagation in systems where mass transfer effects become relevant have been mainly investigated for large void fractions far from the diluted limit [START_REF] Mecredy | The effects of nonequilibrium heat, mass and momentum transfer on two-phase sound speed[END_REF][START_REF] Kieffer | Sound speed in liquid-gas mixtures: Water-air and water-steam[END_REF][START_REF] Ardron | Acoustic wave propagation in a flowing liquid-vapour mixture[END_REF][START_REF] Landau | Fluid Mechanics[END_REF]. For low enough frequencies, [START_REF] Kieffer | Sound speed in liquid-gas mixtures: Water-air and water-steam[END_REF] and [START_REF] Landau | Fluid Mechanics[END_REF] derive limiting expressions assuming that the far field pressure and temperature are related through the saturation conditions. Models accounting for transient heat, mass and momentum transport have been proposed by [START_REF] Mecredy | The effects of nonequilibrium heat, mass and momentum transfer on two-phase sound speed[END_REF], [START_REF] Ardron | Acoustic wave propagation in a flowing liquid-vapour mixture[END_REF] and more recently by [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF].

In diluted systems, mass transfer is expected to influence the local bubble response without significantly influencing the far field temperature as the liquid plays the role of an infinite energy reservoir for the bubble oscillation. In this case, one needs to derive corrected expressions to relate bubble radius and external pressure changes. [START_REF] Hao | The dynamics of vapor bubbles in acoustic pressure fields[END_REF] consider the problem of mass transfer in pure vapor bubbles by assuming saturation inside the bubble at every instant. [START_REF] Preston | A reduced order model of diffusive effects on the dynamics of bubbles[END_REF] solve the linearized problem of the dynamics of air/vapor bubbles for situations where transient vaporization effects and the heat transport in the liquid boundary layer is not a controlling mechanism. [START_REF] Prosperetti | Vapor bubbles in flow and acoustic fields[END_REF] also discuss interesting phenomena induced by mass transfer on the dynamics of bubbles and [START_REF] Prosperetti | The speed of sound in a gas-vapor bubbly liquid[END_REF] presents a simplified model for the influence of mass transfer on the sound-speed of a gas-vapor bubbly liquid that reveals that mass transfer effects play an important role of the propagation of waves at low frequencies.

The process of wave propagation near saturation conditions for diluted systems is less investigated experimentally. [START_REF] Coste | Sound propagation in a liquid with vapour bubbles[END_REF] report evidence of a strong decrease of the sound velocity at low void fractions when approaching the saturation curve in diethylether. Unfortunatelly, the study mainly reports data for large void fractions, where the diluted limit conditions are not met.

In this work we extend the classical linear theory for disperse bubbly liquids to include mass transfer effects across the interface and the diffusion of heat in the liquid. The goals of this study are twofold. First, we want to evaluate the influence of the mass transfer process in the bubble dynamic response and to discuss typical assumptions used in the literature. Second, we want to quantify the importance of mass transfer effects on the effective acoustic properties of the medium. To that end, the full model that accounts for transient mass transfer effects is presented in Section 2. In Section 3 we present the linearization procedure followed to derive the expressions for the bubble resonance frequency, damping factor, mass transfer flux, interface temperature variations and the effective sound phase speed and attenuation of the gas-vapor-liquid mixture (in the following we will simply name them phase speed and attenuation). In Section 4 we present a summary of the model. Section 5 presents some numerical results to gain further insight about mass transfer effects in the process of wave propagation, discussing the relevance of various mechanisms on the bubble response as a function of the vapor content and frequency. Section 6 derives low frequency limiting solutions and the conclusions are presented in Section 7.

Full model

The linearized equations for wave propagation in bubbly media can be derived from basic principles using the conservation equations applied to the averaged mixture. Typically the separation of scales assumption is used to write the averaged equations for the large scale wave propagation problem. Thus, neglecting the influence of liquid viscosity at large scales, we write the mass and momentum conservation equations as

1 ρc 2 Dp Dt + ∇ • v = ∂β ∂t , (2.1) ρ Dv Dt = -∇p, (2.2)
where t is the time coordinate, v is the fluid velocity, p is the pressure and ρ is the average density defined in terms of the void fraction β, the liquid's density, ρ l , and bubble's density ρ b as ρ = (1 -β)ρ l + βρ b . The void fraction β is defined using a probabilistic function for the bubble equilibrium radius, f (a), and the number of bubbles per unit volume, n,

β = 4 3 πn ∞ 0 a 3 f (a)da.
(2.3)

The system of equations above requires to derive an equation to relate external pressure changes with the local bubble radius variation. This equation is found by solving the conservation equations inside and outside the bubble at the local scale. If we impose spherical symmetry for the bubble oscillation, the basic equations for a system of N species can be written as

Dρ Dt + ρ r 2 ∂(v r r 2 ) ∂r = 0, (2.4) ρ DY i Dt = - 1 r 2 ∂ ∂r r 2 j diff i , (2.5) ρ Dv r Dt = - ∂p ∂r + 1 r 2 ∂r 2 τ rr ∂r - τ θθ + τ φφ r , (2.6) ρ De Dt = - p r 2 ∂(v r r 2 ) ∂r - 1 r 2 ∂(r 2 q r ) ∂r + φ v , (2.7)
where r is the radial coordinate, Y i is the mass fraction of the ith component and j diff i is the diffusive mass flux, p is the pressure, T is temperature, τ is the viscous stress tensor, e is the specific internal energy, q r is the radial heat flux and φ v is the viscous dissipation function. These equations apply for both the gas/vapor mixture inside the bubble and the liquid surrounding it.

The ideal gas equation is a good approximation of the real equation of state for both gas and vapor components inside the bubble in standard conditions. In the liquid, the definition of the sound speed in a pure liquid is typically used to relate pressure and density variations. Using the Fourier's and Fick's law to express the diffusive heat and mass flux and neglecting the enthalpy difference between the different components of the mixture at a given temperature we obtain the following set of equations that apply both, inside and outside the bubble,

Dρ m Dt + ρ m r 2 ∂(v m r 2 ) ∂r = 0, (2.8) ρ m DY i,m Dt = 1 r 2 ∂ ∂r r 2 D M i/m ∂Y i,m ∂r ,
(2.9)

ρ m Dv m Dt = - ∂p m ∂r ,
(2.10)

ρ m c p,m DT m Dt = α V T m Dp m Dt + 1 r 2 ∂ ∂r r 2 κ m ∂T m ∂r + 4 3 µ m ∂v m ∂r - v m r 2 , (2.11)
where c p,m is the average specific heat, α V is the thermal dilatation coefficient which is α V = 1/T m for an ideal gas and approximately zero for liquids, κ m is the average conductivity and D M i/m is the diffusion coefficient of the ith component in the mixture. The subindex m is used to denote the average mixture properties. When the equations are applied inside the bubble the subindex m will be replaced by b, when applied into the liquid, we will use the subindex l.

The boundary conditions used to solve the system above can be found for example in [START_REF] Hauke | Dynamics of a single cavitating and reacting bubble[END_REF]. At the bubble center, spherical symmetry imposes the radial gradients to be zero, (i.e. ∂T b ∂r = 0). The bubble velocity is also set to zero, v b (r = 0, t) = 0.

Far from the bubble, pressure, temperature and all the species concentrations in the liquid are assumed to be known (p l,∞ , T 0 , Y i,l,0 ).

The local balances at the infinitely thin interface relate the liquid and gas properties at both sides of the interface. The mass balance across the bubble is

J = ρ g ( Ṙ -v I b ) = ρ l ( Ṙ -v I l ) at r = R, (2.12)
where J is the total evaporation mass flux across the interface, R is the bubble radius, Ṙ is the interface velocity, and v I b and v I l denote respectively the fluid interface velocity in the bubble and in the liquid. The momentum balance accounting for mass transfer effects reads

J v I b -v I l = p I l -p I b + τ rr,l -τ rr,b + 2σ R at r = R, (2.13)
where σ is the surface tension coefficient.

The total flux of the ith component across the interface is calculated as the sum of the advective flux and the diffusive flux. As both fluxes must be equal at r = R, we obtain an equation for the conservation of the ith component across the interface as

-JY i,l -ρ l D M i/l ∂Y i,l ∂r = -JY i,b -ρ g D M i/b ∂Y i,b ∂r at r = R.
(2.14)

The energy balance at the interface establishes the relation among the energy fluxes in the liquid/gas boundaries. This equation can be approximated by using the total latent heat of evaporation, ∆H vap , as

κ l ∂T l ∂r = κ g ∂T g ∂r + N i=1 J i ∆H vap,i at r = R, (2.15)
where J i is the flux across the interface of the ith component and ∆H vap,i is the enthalpy of vaporization related to the phase change. Finally, we assume a continuous temperature profile across the interface so that the interface temperature is equal for both phases

(T g (r = R) = T l (r = R) = T int ).
To evaluate the flux of every component across the interface J i that finally determines the total flux J = N i=1 J i we can either impose equilibrium conditions at the interface at every instant, or assume that the flux is proportional to the difference between the equilibrium state and the current state. In fact, using the kinetic theory of gases, it is possible to obtain the total flux of vapor across the interface using the Hertz-Knudsen-Langmuir expression [START_REF] Knudsen | Maximum rate of vaporization of mercury[END_REF][START_REF] Hertz | Über die Verdunstug der Flüssigkeiten, Inbesondere des Quecksilbers im lufteren Räume (on the evaporation of fluids, especially of mercury, in vacuum spaces)[END_REF],

J vap = α evap (p I eq -p I b,vap ) 2πr vap T int (2.16)
where α evap is the accommodation coefficient, which is a measure of the ratio of the molecules hitting the interface that change the phase, p I b,vap is the partial pressure of vapor at the interface, r vap is the vapor's perfect gas constant and p I eq is the equilibrium pressure at the interface conditions. For the vapor, the equilibrium pressure can be obtained from the Clausius-Clapeyron relation, dp eq dT = ∆H vap p eq r vap T 2 .

(2.17)

It is interesting to remark that the system of equations above simplifies in cases where the interface is assumed to be at equilibrium at every instant. As we will prove along the theoretical development, for low enough frequencies, transient mass transfer effects (Eq. 2.16) are not relevant. In this situation the total flux is given by diffusion assuming equilibrium conditions at the interface at every instant. Analogously, the model could be easily extended for soluble gases by using the equilibrium pressure (or concentration) given by the Henry's law. It must be noted that in this work we will not consider the presence of soluble gases in the liquid given that the mass and heat flux related to phase change of soluble components is usually negligible compared with that of vapor.

Linearized solutions

The set of equations described above can be simplified assuming that all variables oscillate around an equilibrium state with a given frequency ω. Thus, for a given variable y, we look for solutions of the form y = y 0 (1 + ∆ye ıωt ), where ∆y is a complex quantity to be determined. For the linear solution to be valid, ∆y ≪ 1.

At the local scale, the momentum and continuity equations in the liquid can be rewritten as the Rayleigh-Plesset equation, which relates pressure (or potential) variations induced by the far field with the local response of the bubble. Neglecting compressibility effects the modified Rayleigh-Plesset equation accounting for mass transfer effects reads

(Prosperetti 1982) R R - R J ρ l + 3 2 Ṙ - J ρ l 2 -2 J ρ l Ṙ - J ρ l = p I l -p l,∞ ρ l , (3.1)
where p I l can be expressed in terms of the properties inside the bubble using Eq. 2.13. In the following, we will consider the system of a bubble containing an inmiscible gas and a vapor with an equilibrium concentration Y 0 in a monocomponent liquid. This is representative of most of the gas/liquid systems as the most important flux across the interface is expected to be given by liquid vaporization. In the linear limit the flux across the interface is J = J 0 ∆Re ıωt , where J 0 is a complex quantity to be determined. Neglecting non-linear terms we simplify the modified Rayleigh-Plesset equation as

-R 2 0 ω 2 ∆R - R 0 J 0 ıω ρ l ∆R = ∆R ρ l -p b,0 Φ + 2σ R 0 -4µ l ıω - p l,0 ∆p ∞ ρ l , (3.2)
where Φ is a complex function to be obtained that is typically named transfer function and it serves to relate bubble radius oscillations and internal bubble pressure variations, such that ∆p b = -Φ∆R.

The equation above can be rearranged to express the bubble radius variations in terms of the far field pressure variations as

∆R = - 1 ω 2 0 -ω 2 + 2iδω ∆p ∞ p l,0 ρ l R 2 0 , (3.3)
where

δ = 2µ l ρ l R 2 0 + p g,0 2ωρ l R 2 0 ℑ(Φ) + Re(J 0 ) 2ρ l R 0 , (3.4) ω 2 0 = p g0 ρ l R 2 0 Re(Φ) - 2σ p g,0 R 0 + ωℑ(J 0 ) ρ l R 0 . (3.5)
Eqs. 3.4 and 3.5 require the evaluation of the transfer function Φ and the flux across the interface J 0 . Both variables can be obtained from the solution of the conservation equations inside the bubble. Following a development similar to that of [START_REF] Commander | Linear pressure waves in bubbly liquids: Comparison between theory and experiments[END_REF], we consider a perfect gas and uniform pressure inside the bubble. The continuity equation inside the bubble can be rewritten using the energy equation as

ṗb γp b + ∇ • v b - γ -1 γp b κ b ∇T b = 0. (3.6)
Integrating the equation from the bubble center to a given distance r, the radial velocity inside the bubble is given by

v b (r) = 1 γp (γ -1)κ b ∂T ∂r - 1 3 ṗb r . (3.7)
This equation is evaluated at the interface to find an expression for the bubble's internal pressure,

ṗb = 3 R (γ -1)κ b ∂T b ∂r r=R - 3 R γp b v I b . (3.8)
Using Eq. 2.12 to express the gas/vapor velocity at the interface (v I b ) as a function of the total evaporation flux and the interface velocity and linearizing we obtain

∆p b = -3(γ -1)ıPe -1 b γ γ -1 ∂∆T b ∂ζ ζ=1 -3γ∆R -3γıJ * 0 ∆R, (3.9) where ζ = r/R is the nondimensional radial distance, Pe b = ωR 2 0 D T b
is the bubble's Peclet number defined using the bubble's thermal diffusitivity D T b , and the nondimensional mass transfer flux is

J * 0 = J 0 ρ b,0 R 0 ω . (3.10)
Assuming that only one component (denoted with the subindex w) is vaporized and that the vapor pressure solely depends on temperature, we can evaluate the mass transfer flux from the linearization of Eq. 2.16 to obtain

J * 0 ∆R = J * max ∂p eq ∂T T 0 p eq ∆T b (ζ = 1) -∆p b + ∆Y (ζ = 1) , (3.11)
where J * max is a nondimensional flux defined as

J * max ≡ α evap p eq,0 ρ b,0 R 0 ω √ 2πr w T 0 .
(3.12) Equation 3.11 requires to know the vapor mass fraction variations and the temperature variations at the interface. To obtain ∆Y (ζ = 1) we solve the linear transport equation for the vapor inside the bubble. The transport equation (Eq. 2.9) is written in a nondimensional form using the Sherwood number, Sh D = ωR 2 0 D w/b , where D w/b is the diffusion coefficient of the vapor inside the bubble mixture,

Sh D ∆Y + i 1 ζ 2 ∂ ∂ζ ζ 2 ∂∆Y ∂ζ = 0. (3.13) Its general solution, ∆Y (ζ) = C D sinh √ Sh D ıζ ζ , (3.14)
is particularized using the boundary condition for component w at the interface (Eq. 2.14)

-J * 0 ∆R = -J * 0 Y 0 ∆R -Sh -1 D Y 0 ∂∆Y ∂ζ ζ=1 , (3.15)
where we have imposed Y w,l = 1. Thus, we can determine the constant C D to write the evolution of the mass fraction at the interface as

∆Y (ζ = 1) = J * 0 Sh D √ ıSh D coth( √ ıSh D ) -1 1 -Y 0 Y 0 ∆R. (3.16)
The temperature variations at the interface are obtained using the energy equation inside the bubble (Eq. 2.11). The nondimensional linear version of this equation,

∆T b -∆p b γ -1 γ + ıPe -1 b 1 ζ 2 ∂ ∂ζ ζ 2 ∂∆T b ∂ζ = 0, (3.17)
has the following general solution,

∆T b (ζ) = ∆p b γ -1 γ + C 1 exp - √ Pe b ıζ ζ + C 2 exp √ Pe b ıζ ζ . (3.18)
Because the temperature is finite at ζ = 0, we find that C 1 = -C 2 = -C 0 /2. Thus, the interface's temperature variation is given by

∆T b (ζ = 1) = ∆p b γ -1 γ + C 0 sinh Pe b ı . (3.19)
Plugging Eqs. 3.16 and 3.19 into Eq. 3.11 we obtain the following explicit relation for the mass flux across the interface

J * 0 ∆R = J * c ∆H * vap -1 ∆p b + ∆H * vap γ γ -1 C 0 sinh Pe b ı , (3.20)
where we have defined a nondimensional enthalpy of vaporization as ∆H * vap = ∆Hvap c p,b T0 and a characteristic nondimensional flux

J * c = J * max 1 + J * max ShD √ ıShD coth( √ ıShD)-1 1-Y0 Y0 . (3.21)
Note that this quantity provides a quantification about the influence of transient mass transfer effects. In fact, we can write Eq. 3.21 using the characteristic nondimensional mass flux, J * c,eq that is obtained by assuming that equilibrium conditions prevail at the interface at every instant (see appendix A for details about the derivation of J * c,eq ),

J * c = J * max 1 + J * max (J * c,eq ) -1 . (3.22)
When the characteristic mass flux related to transient effects is much faster than the one at equilibrium conditions (J * max ≫ J * c,eq ), the mass transfer is controlled by the vapor diffusion inside the bubble: J * c ≈ J * c,eq . In the opposite limit (J * max ≪ J * c,eq ), transient mass transfer effects control the mass flux across the interface:

J * c ≈ J * max .
To obtain the value of the integration constant C 0 in Eq. 3.20 one needs to solve the energy equation in the liquid (Eq. 2.11). In its linearized form this equation reads

∆T l = -ıPe -1 l 1 ζ 2 ∂ ∂ζ ζ 2 ∂∆T l ∂ζ , (3.23)
and its general solution is

∆T l (ζ) = C 3 exp - √ Pe l ıζ ζ + C 4 exp √ Pe l ıζ ζ , (3.24)
where Pe l =

ωR 2 0 D T l
is the liquid's Peclet number. Note that consistent with the derivation of the far field equations, the effect of viscous dissipation as well as thermal dilatation on the temperature variations have been neglected. The particular solution of this equation requires to evaluate the liquid's temperature at some distance ζ ∞ . In the dilute limit, one can consider that the temperature far from the bubble is constant and equal to T 0 . In this case we find C 4 = 0.

The value of the free constants C 3 and C 0 in Eqs. 3.19 and 3.24 are obtained applying the boundary conditions at the interface. First, we assume that the interface's temperature has to be continuous. Thus we write,

C 3 exp -Pe l ı = ∆p b γ -1 γ + C 0 sinh Pe b i . (3.25)
In addition, the energy balance across the interface (Eq. 2.15) establishes that

∂∆T l ∂ζ ζ=1 = κ b κ l ∂∆T b ∂ζ ζ=1 + J * 0 Pe b ∆H * vap ∆R. (3.26)
Replacing the corresponding expressions for the temperature derivatives and using Eq. 3.25 to express C 3 as a function of C 0 we obtain

C 0 = - γ-1 γ (1 -∆T I c ) sinh √ Pe b i ∆p b (3.27)
where the complex quantity ∆T I c is defined as

∆T I c = κ b κ l √ Pe b ı coth( √ Pe b ı) -1 + J * c Pe b ∆H * vap γ γ-1 1 + √ Pe l ı + κ b κ l √ Pe b ı coth( √ Pe b ı) -1 + J * c Pe b ∆H * vap 2 γ γ-1 . (3.28)
Upon substitution of the integration constant into Eq. 3.25 we can represent the temperature variations at the interface as

∆T b (ζ = 1) = γ -1 γ ∆T I c ∆p b . (3.29)
The temperature gradient inside the bubble evaluated at the interface is obtained from Eq. 3.19,

∂∆T b ∂ζ ζ=1 = - γ -1 γ (1 -∆T I c ) Pe b ı coth( Pe b ı) -1 ∆p b , (3.30)
and the mass flux across the interface (Eq. 3.20) is expressed as a function of the bubble pressure variations as

J * 0 ∆R = J * c ∆H * vap ∆T I c -1 ∆p b . (3.31)
Replacing the expressions for the temperature gradient and mass flux (Eqs. 3.30 and 3.31) into Eq. 3.9 we finally obtain the transfer function Φ as

Φ = 3γ 1 -3(γ -1)ıPe -1 b √ Pe b ı coth( √ Pe b ı) -1 (1 -∆T I c ) + 3γıJ * c ∆H * vap ∆T I c -1 (3.32)
This transfer function is a generalization of the classical expression for pure gas bubbles [START_REF] Prosperetti | Nonlinear bubble dynamics[END_REF]) accounting for interface temperature variations and mass transfer effects across the interface. For instance, it is easy to check that we recover the classical solution for pure gas bubbles by imposing no net mass flux across the interface (J * c = 0) and by assuming that the interface temperature variations are zero (∆T I c = 0).

Figure 1 depicts the real and imaginary part of the transfer function obtained for an air/vapor bubble for various water vapor contents at constant atmospheric pressure. At large frequencies the bubbles tend to recover the adiabatic response irrespective of the vapor content. At low frequencies, the transfer function recovers the isothermal limit for pure gas bubbles. As the amount of vapor increases, the real part takes values below 1 at low frequencies. This value corresponds to the limiting case in which equilibrium conditions prevail inside the bubble at every instant. It can be proven that this result can also be found assuming that vapor and gas are both ideal gases at constant temperature (see appendix B). The imaginary part is also influenced by the vapor content. In general The vapor content for the various curves are Y 0 = 0, 0.03, 0.27, 0.50, 0.74, 0.99. the imaginary part increases with the vapor content irrespective of the Peclet number although the influence of the vapor content on the imaginary part is especially notorious at low frequencies. Consistent with experimental observations, the influence of vaporization is completely negligible in gas/vapor bubbles at standard laboratory conditions (25 degree Celsius and 1 atmosphere) where the vapor void fraction is negligible compared with that of a gas.

Once we have obtained expressions for the transfer function (Eq. 3.32) and the mass flux across the interface (Eq. 3.31) we follow the classical procedure to linearize the averaged continuity and momentum equations (Eqs. 2.1 and 2.2) written as a wave equation (3.33) where R = -ω 2 R 0 ∆R is the interface acceleration. Rewriting Eq. 3.33 as a wave equation for p l,∞ , the complex sound speed, c m , is given by

of the form 1 c 2 eff ∂ 2 p l,∞ ∂t 2 -△p l,∞ = 4πρ ∞ 0 a 2 Rf (a)da,
1 c 2 m = 1 c 2 + 4πn ∞ 0 af (a)da ω 2 0 -ω 2 + 2iδω (3.34)
from which we obtain the phase velocity, c ph , and attenuation, Q -1 , as

c ph = Re 1 c m -1 , (3.35) Q -1 = -20 log 10 (e)ℑ ω c m . (3.36)
Using the expressions for the bubble resonant frequency and damping constant (Eqs. 3.4-3.5) we obtain the phase speed (Eq. 3.35) and attenuation (Eq. 3.36).

Summary of equations

As a summary of the full model equations for a bubbly liquid with vapor and an inmiscible gas, we solve for the averaged complex sound speed

1 c 2 m = 1 c 2 + 4πn ∞ 0 af (a)da ω 2 0 -ω 2 + 2iδω , (4.1)
in order to obtain the phase speed and attenuation as

c ph = Re 1 c m -1 , (4.2) Q -1 = -20 log 10 (e)ℑ ω c m . (4.3)
The second term in Eq. 4.1 represents the influence of the bubble oscillation on the large scale wave propagation problem and it can be obtained in the linear regime by solving the conservation equations at the local scale (single bubble and its surrounding liquid).

From the linearization of the Rayleigh-Plesset equation (e.g. continuity and momentum equation in the surrounding liquid) we find the following expressions for the bubble resonant frequency and the damping coefficient

δ = 2µ l ρ l R 2 0 + p g,0 2ωρ l R 2 0 ℑ(Φ) + 1 2 Re(J * 0 ) ρ b ρ l ω, (4.4
)

ω 2 0 = p g0 ρ l R 2 0 Re(Φ) - 2σ p g,0 R 0 + ℑ(J * 0 ) ρ b ρ l ω 2 . (4.5)
These properties require to evaluate the nondimensional mass transfer flux J * 0 and the transfer function Φ that relates bubble's pressure and volume changes. From the solution of the conservation equations inside the bubble and the energy equation in the surrounding liquid we find the following expressions for these two quantities reveals that the total flux across the interface can be controlled either by vapor diffusion inside the bubble, J * c,eq , or by transient vaporization effects across the interface, J * max . The diffusion flux

J * 0 = J * c 1 -∆H * vap ∆T I c Φ, (4.6) Φ = 3γ 1 -3(γ -1)ıPe -1 b √ Pe b ı coth( √ Pe b ı) -1 (1 -∆T I c ) + 3γıJ * c ∆H * vap ∆T I c -1 , ( 4
J * c,eq = √ ıSh D coth( √ ıSh D ) -1 Sh D Y 0 1 -Y 0 , ( 4
J * max = α evap p eq,0 ρ b,0 R 0 ω √ 2πr w T 0 . (4.10)
Finally, the nondimensional interface temperature variation is a complex function depending on the heat and mass transfer processes taking place at the local scale and its general solution is

∆T I c = κ b κ l √ Pe b ı coth( √ Pe b ı) -1 + J * c Pe b ∆H * vap γ γ-1 1 + √ Pe l ı + κ b κ l √ Pe b ı coth( √ Pe b ı) -1 + J * c Pe b ∆H * vap 2 γ γ-1
, (4.11)

where Pe l =

ωR 2 0 D T l .
In the following sections we discuss the importance of mass transfer in the bubble dynamic response and the wave propagation properties of the effective medium as well as simplified solutions for limiting situations.

The influence of mass transfer effects on the acoustic properties of the effective medium

In this section we evaluate the importance of mass transfer effects in a system by comparing the results obtained with the full model with those provided by the classical linear theory that neglects mass transfer effects. We consider a mixture of monodisperse bubbles with average radius 100 µm and bubble concentration of β 0 = 10 -4 (vol/vol) measured at standard temperature (25 degrees Celsius) and pressure (1 atmosphere). We take the physical properties of an air/water system. Because the exact value of the accommodation coefficient is not clear yet [START_REF] Gumerov | Determination of the Accomodation Coefficient Using Vapor/Gas Bubble Dynamics in an Acoustic Field[END_REF] in this work we have taken α evap = 0.35 [START_REF] Yasui | Alternative model of single sonoluminiscence[END_REF][START_REF] Hauke | Dynamics of a single cavitating and reacting bubble[END_REF][START_REF] Fuster | Influence of accommodation coefficient on nonlinear bubble oscillations[END_REF], where we have neglected any influence of the temperature on this parameter.

Figure 2 contains the effective phase velocity and attenuation as a function of the vapor content and frequency obtained with the full model. The vapor content is controlled by varying the ambient temperature from 25 degrees Celsius to the normal boiling point of water (100 degrees Celsius) at constant ambient pressure. Under these conditions, the classical theory that neglects mass transfer effects does not predict any significant influence of the vapor content (or ambient temperature) on the results obtained. However, the results obtained with the full model reveal that mass transfer effects become relevant at low frequencies for bubbles with a large vapor content. At low frequencies the phase velocity is significantly reduced with respect to the value predicted when neglecting mass transfer effects whereas attenuation values increase various order of magnitude. The analysis of low frequency asymptotic solutions reached for low frequencies is postponed to the next section.

The differences in the formulation of the full model including mass transfer effects with respect to the classical formulation for pure gas bubbles can be condensed in two nondimensional quantities: J * 0 as defined in Eq. 3.10, which is interpreted as the ratio between the interface velocity induced by mass transfer J0 ρ b ∆R and the actual interface velocity R 0 ω∆R, and ∆T I c as defined in Eq. 3.28, which is a measure of the nondimensional temperature variations at the interface (Eq. 3.29). When J * 0 is close to one mass transfer controls of the evolution of the bubble radius, whereas when it is smaller than one the gas/vapor expansion controls the bubble radius evolution.

Figure 3 shows the modulus and the phase lag with respect to the bubble radius oscillations of the nondimensional mass transfer flux as a function of the Sherwood number and the vapor content (varied changing the ambient temperature from 25 to 100 degrees Celsius at 1 atmosphere). For a given vapor content, the evaporation flux reaches an asymptotic limit for low Sherwood numbers (low frequencies) which corresponds to the limit obtained assuming thermodynamic equilibrium inside the bubble at every instant (appendix B). This low frequency limit is also the maximum mass transfer flux across the interface for a given vapor content. Only for systems approaching saturation conditions the nondimensional mass transfer flux becomes of order one, meaning that evaporation and condensation processes take control of the dynamic response of the interface. On the contrary, for Sh D >> 1 the bubble has almost no time to respond to pressure waves and the evaporation flux becomes negligible. In this high frequency regime neither phase speed nor attenuation are significantly influenced by mass transfer effects and we can conclude that the bubble dynamic response is governed by the gas/vapor expansion/compression. The relevance of the interface temperature variations are captured by the nondimensional quantity ∆T I c . Figure 4 shows ∆T I c as a function of the bubble's Peclet number, the nondimensional enthalpy of vaporization and the amount of vapor, which ultimately controls the mass flux across the interface. As expected, the interface temperature variations are negligible when the vapor content is low, irrespective of the forcing frequency and enthalpy of vaporization. In this case, the energy exchange due to evaporation/condensation is not significant. As we approach saturation conditions, the mass flux across the interface and the associated energy exchange make the interface temperature variations important. The sensitivity is maximal for a given enthalpy of vaporization that slightly depends on the forcing frequency. For reference, in the particular case of water the nondimensional enthalpy of vaporization, ∆H * vap , takes values of order one, meaning that the interface temperature variations become important near saturation conditions. Regarding the frequency, the interface temperature variations are not important at low frequencies for both gas and vapor bubbles as the bubble has enough time to equilibrate its temper- Influence of the vapor content on the phase velocity (top) and attenuation (bottom) curves obtained from the set of equations summarized in Section 4 for gas/vapor bubbles at 1 atm and various ambient temperatures. The vapor content for each condition is Y 0 = 0.99, 0.96, 0.89, 0.55, 0.03. For reference, the red line represents the solution of pure air bubbles (Y 0 = 0) at 1 atm and 25 degrees Celsius. Mass transfer effects significantly modify both phase speed and attenuation curves for frequencies below resonance. In the low frequency limit, mass transfer decreases the phase velocity and increases attenuation.

ature with the surrounding liquid at every instant. As we increase the frequency (and therefore the Peclet number) the interface temperature variations are more important for bubbles with a low vapor content although in all cases the nondimensional values are close to zero. Thus, we can safely assume that the interface temperature is constant for gas bubbles in water irrespective of the forcing frequency. For bubbles with a large vapor content the situation is different because the interface temperature is influenced by the mass transfer flux. As this flux decreases when increasing the forcing frequency, there exists a range of intermediate frequencies for which the interface temperature variations become important in order to correctly predict the bubble radial oscillations. In order to gain further insight into the applicability of the constant interface temperature assumption Figure 5 compares the solution obtained with the full model with the solution obtained by impossing that the interface temperature remains constant (∆T I c = 0) still considering the mass transfer flux across the interface. The results are shown for a gas/vapor bubble at 25 degree Celsius (approximately 0.03% of vapor content) and 99.9 degree celsius (Y 0 = 0.99). While at low ambient temperatures (low vapor content) the results obtained with ∆T I c = 0 accurately represent the solution obtained with the full model, we do observe significant differences when approaching saturation conditions. In this case the vaporization flux is significant and it is important to include interface temperature variations for an accurate estimation of the phase speed and attenuation at frequencies below the resonance frequency. Only if we decrease frequencies further down, the bubble response can be assumed to be isothermal and the phase velocities obtained assuming constant temperature accurately represent the solution of the full model. Note that this is not the case of the attenuation, where the interface temperature variations seem to have a significant influence on the predictions and must then be always consid- ered in the predictions.

Low frequency limits

We have seen in the previous section that mass transfer effects are especially important for low frequency excitations. In this case both the bubble response and the process of wave propagation are significantly influenced by the mass and energy exchange across the interface. In order to gain further insight into the influence of mass transfer at low frequencies we derive limiting expressions for particular regimes. For frequencies well below the natural frequency and assuming monodisperse mixtures, it is readily shown that [START_REF] Commander | Linear pressure waves in bubbly liquids: Comparison between theory and experiments[END_REF])

1 c 2 m = 1 c 2 + 4π nR 0 ω 2 0 1 -2ıδω ω 2 0 . ( 6
Figure 5: Color online. Phase speed (top) and attenuation curves (bottom) for an air/water system with 100 µm bubbles with concentration β 0 = 10 -4 at 25 C (red lines, Y 0 = 0.03) and 99.9 C (blue lines, Y 0 = 0.99). The continuous lines represent the solution of the full model. The dashed lines represent the solution assuming a constant interface temperature. Both solutions overlap for systems containing gas bubbles whereas we observe significant differences between them for low frequency waves propagating in systems with bubbles with a large vapor content.

The limiting expressions for ω 0 and δ depend on the evaporation flux and the transfer function. The nondimensional flux across the interface, given by Eq. 3.31, is J * 0 = -J * c Φ. In the general case J * c has to be obtained from Eq. 3.22, which in the low frequency limit tends to

J * c ≈ J * max 1 -3J * max ı 1-Y0 Y0 when Sh D < 1. (6.2)
This expression can be further simplified when J * max > 1 3 Y0 1-Y0 , which can be written as a condition for the bubble's Peclet number,

Pe b < Pe I c = 3α evap √ 2π (1 -Y 0 ) R 0 √ r w T D T b .
(6.3) The right hand side (Pe I c ) is identified with a critical Peclet number below which transient mass transfer effects are negligible. In this case equation 6.2 reduces to

J * c ≈ 1 3 ı Y 0 1 -Y 0 when Sh D < 1 ∧ Pe b < Pe I c , (6.4)
which is equivalent to the low frequency solution obtained when assuming equilibrium conditions at the interface (see appendix A). Thus, we write the mass transfer flux across the interface and the transfer function as

J * 0 ≈ - 1 3 ı Y 0 1 -Y 0 Φ, (6.5) Φ ≈ 3 (1 -Y 0 ) 1 -(1 -Y 0 ) ∆T I c γ-1 γ -Y 0 ∆H * vap ∆T I c (6.6)
which are valid when Sh D < 1 ∧ Pe b < min(1, Pe I c ). Note that to derive the expression for the transfer function we have approximated the function M (Pe) = √ Peı coth( √ Peı) -1 by M (Pe) ≈ Pe 3 ı, which is a reasonable assumption for Pe b < 1 (Figure 6).

Both the mass transfer flux and the transfer function are a function of the interface temperature variations given by Eq. 3.28. For low frequencies, Pe b < 1, Eq. 3.28 simplifies to

∆T I c ≈ κ b κ l 1 3 ı 1 + Y0 1-Y0 ∆H * vap γ γ-1 Pe b 1 + √ Pe l ı + κ b κ l 1 3 ı 1 + Y0 1-Y0 ∆H * vap 2 γ γ-1 Pe b (6.7) when Sh D < 1 ∧ Pe b < min(1, Pe I c )
When the interface temperature is isothermal, ∆T I c = 0, which is true for Pe b → 0, the nondimensional mass flux across the interface and the transfer function for bubbles containing gas and vapor are

J * 0 = Y 0 ı (6.8) Φ ≈ 3(1 -Y 0 ). (6.9)
Using these limiting expressions the final expressions for the phase speed and attenuation in the low frequency limit are given by

c 2 ph = c 2 1 + c 2 β0ρ l (1-Y0)p0
, (6.10)

Q -1 = 20 log 10 (e) cβ 0 δω 2 ρ 2 l R 2 0 3p 2 0 (1 -Y 0 ) 2 1 + c 2 β0ρ l (1-Y0)p0
. (6.11) For low void fractions (β 0 < p0(1-Y0)

ρ l c 2
), the phase speed tends to recover the liquid's sound speed. On the contrary, for large void fractions the phase speed is shown to be proportional to the gas content inside the bubble.

These limiting expressions are only expected to represent reasonably well the solution for dilute systems when Sh D < 1 and Pe b < min(1, Pe I c ) and also when the interface temperature remains constant. This last condition introduces an additional constrain for Pe b in order to reach the low frequency limits obtained above. From Eq. 6.7, we conclude that we can only assume constant interface temperature when

Pe b ≪ Pe II b,c = 3 1 -Y 0 Y 0 κ l κ b γ -1 γ 1 (∆H * vap ) 2 .
(6.12)

Thus, the frequency required to reach the low frequency limit is expected to be lower in systems with a low liquid thermal conductivity, a large content of vapor in the bubbles and a large enthalpy of vaporization. This effect can be clearly seen in Figure 7 where we compare the solution of systems with different non-dimensional values of the enthalpy of vaporization with the solution of a model where the interface temperature is assumed to be isothermal. As expected, all systems tend to the same low frequency limit for the phase speed, although we need to low the forcing frequency as the enthalpy of vaporization increases. In particular, in the case of the highest enthalpy of vaporization tested here (twice the enthalpy of vaporization of water) the curves tend to follow the low frequency limit for gas bubbles at frequencies between 10 3 and 10 4 Hz. At these frequencies, the vaporization flux is limited by heat diffusion from the bulk towards the interface and the bubble dynamic response is mainly controlled by the expansion/compression rather than mass transfer effects. For lower frequencies, heat has enough time to diffuse and vaporization controls the effective compressibility of the medium. This remark is important if one wants to guarantee the validity of the low frequency limiting solutions in systems near saturation conditions, where it is not sufficient to satisfy that the excitation frequency is below the bubble's resonance frequency.

It is also interesting to compare the low frequency limit solutions obtained here with the classical solution for gas bubbles proposed by [START_REF] Wood | A Textbook of Sound[END_REF]. Figure 8 compares the results obtained with the full model at 1 Hz as a function of the void fraction and the vapor content with the results obtained from Wood's theory. As it can be seen the phase speed is significantly lower for bubbles with a large vapor content than for pure gas bubbles. The effect of the vapor content can be also seen in Figure 9 where we see how both phase speed and attenuation converge to the solution for pure gas bubbles with vapor content tending to zero. For the phase velocity Wood's formula provides a good approximation of the exact solution only when the void fraction is sufficiently low. In the case of attenuation (Figure 9), δ strongly depends on the vapor content increasing . For reference we include the solution of pure gas bubbles and the approximate solution provided by Eqs. 6.10-6.11.

several orders of magnitude compared to the case of pure gas bubbles.

It is also worth mentioning that the limiting solutions of Eqs. 6.10-6.11 differ from the solutions obtained by [START_REF] Landau | Fluid Mechanics[END_REF] for pure liquid/vapor bubbly liquids, who derive an asymptotic solution considering that pressure and temperature at the large scale are coupled through the Clausius-Clapeyron equation. In this study we restrict ourselves to situations where the temperature far from the bubble remains constant, which is a reasonable assumption for dilute systems containing small amount of gas but it is not applicable for pure vapor bubbles and low frequencies where the heating/cooling of the bulk liquid may be relevant. Although the current model for pure vapor bubbles may be extended accounting for large scale thermal effects for the sake of simplicity we postpone the study of this situation for future investigation. In any case, it is important to note that the limiting solutions reported by [START_REF] Landau | Fluid Mechanics[END_REF] may be very difficult to obtain experimentally in linear regimes [START_REF] Coste | Sound propagation in a liquid with vapour bubbles[END_REF]). On one side, very small amounts of gas make the vapor diffusion inside the bubble to be the controlling mechanism determining the vaporization flux. As we have seen above, the frequency threshold below which equilibrium conditions prevail inside the bubble at every instant tend to zero for vapor bubbles, which would make difficult to reach the low frequency limit in these situations. On the other side, for the particular case of pure vapor bubbles, the amplitude required to keep the linear regime valid tends to zero for low frequencies. This effect can be clearly seen in Figure 10, where we represent the effect of the vapor content on the sensitivity factor S defined from Eq. 6.13 as .13) This factor represents the nondimensional amplification factor of the bubble radius oscillation with respect to the external pressure excitation. Given that a pure vapor bubble is unable to reach equilibrium conditions when we perturb the pressure, the bubble tends to infinitely grow (or shrink) as we decrease the frequency. Small amounts of gas serve to kill this singular behavior, although the sensitivity of the bubble radius oscillation against the pressure pulse is still significantly influenced by the amount of vapor. Sensitivity factor (as defined in Eq. 6.13) as a function of the vapor content and frequency for a 100 µm gas/vapor bubble. The vapor content is changed by modifying the ambient temperature. As the amount of vapor increases, the bubble radial oscillation is more important for a constant nondimensional pressure wave amplitude.

S = ∆R ∆p ∞ = - 1 ω 2 0 -ω 2 + 2iδω p l,0 ρ l R 2 0 . ( 6

Conclusions

This manuscript presents a generalization of the classical linear theory for dilute bubbly liquids that solves for the transient mass transfer effects across the interface as well as the heat transport in the liquid boundary layer. The proposed linear model converges to the solution of the classical linear theory for pure gas bubbles when the vapor content inside the vapor is negligible. Mass transfer effects are shown to play an important role in the phase speed and attenuation curves at frequencies below the bubble resonant frequency. At low frequencies the phase speed is significantly reduced while attenuation can increase various orders of magnitude with respect to the solution obtained when neglecting mass transfer effects.

Simplified solutions are derived for low frequency excitations. The conditions to be satisfied for these solutions to be valid are also obtained. We show that to reach the low frequency limit for the phase speed it is not sufficient to guarantee that the forcing frequency is lower than the bubble resonance frequency. For a general case, the Sherwood number, Sh D = , must be smaller than the minimum among two characteristic Peclet numbers, Pe I c and Pe II c , in order to guarantee that thermodynamic equilibrium prevails at every instant. In systems with low liquid conductivity, large enthalpy of vaporization and large vapor content, the frequencies required to reach the isothermal limit for situations with a large vapor content can be extremely low. In this low frequency regime the bubble response is isothermal and it is possible to derive expressions that capture relatively well the influence of the vapor content on the acoustic properties of the effective medium. Thus, in systems close to saturation conditions and for frequencies below resonance, it is important to solve for the mass transfer flux across the bubble interface, the heat transport equation in the liquid and to account for the interface temperature variations in order to accurately 8. Appendix A: Mass transfer flux assuming equilibrium conditions at the interface It is commonly assumed that equilibrium conditions prevail at the interface at every instant. Under these conditions, the partial water vapor pressure at the interface is given at every instant by the saturation pressure at the interface's temperature, p eq (T 9. Appendix B: Low frequency limit of the mass transfer flux assuming thermodynamic equilibrium inside the bubble

The limiting solution of the mass transfer flux for low frequencies can be also obtained by assuming that: (i) the temperature of the bubble remains constant at every instant and equal to the reference temperature, (ii) the vapor pressure inside the bubble is maintained at every instant. Thus, the water vapor mass can be obtained as m w = p w V r w T (9.1)

Because only the bubble volume changes, the flux is obtained as

J = ṁw 4πR 2 b = p w r w T Ṙb . (9.2)
Linearizing the expression above we find (9.3) which in nondimensional form reads (9.4)

J 0 ∆R = m w V ωıR 0 ∆R,
J * 0 = J 0 R b ρ b,0 ω = Y 0 ı.

  ∂Y i,b ∂r = ∂ρ b ∂r =

Figure 1 :

 1 Figure 1: Color online. Influence of the vapor content on the real (top) and imaginary (bottom) part of the transfer function (Eq. 3.32) for various Peclet numbers. The values of the nondimensional quantities J * c , ∆H * vap , ∆T I c are obtained for a 100 µm air/water vapor bubble at 1 atmosphere and varying temperature (from 25 to 100 degrees Celsius).The vapor content for the various curves are Y 0 = 0, 0.03, 0.27, 0.50, 0.74, 0.99.

  b T0 , and J * c and ∆T I c are the characteristic nondimensional flux and the nondimensional interface temperature variation. The general expression for the nondimensional flux J *

  .9) depends on the Sherwood number Sh D = ωR 2 0 D w/b and the amount of vapor Y 0 . The characteristic transient mass flux is

Figure 2 :

 2 Figure2: Color online. Influence of the vapor content on the phase velocity (top) and attenuation (bottom) curves obtained from the set of equations summarized in Section 4 for gas/vapor bubbles at 1 atm and various ambient temperatures. The vapor content for each condition is Y 0 = 0.99, 0.96, 0.89, 0.55, 0.03. For reference, the red line represents the solution of pure air bubbles (Y 0 = 0) at 1 atm and 25 degrees Celsius. Mass transfer effects significantly modify both phase speed and attenuation curves for frequencies below resonance. In the low frequency limit, mass transfer decreases the phase velocity and increases attenuation.

Figure 3 :

 3 Figure 3: Color online. Top: Nondimensional evaporation flux (Eq. 4.6) as a function of the vapor content (varied by changing the ambient temperature at a constant pressure equal to 1 atm) and the Sherwood number. Bottom: Phase lag between the bubble radius and flux oscillation as a function of the vapor content and the Sherwood number.

Figure 4 :

 4 Figure 4: Color online. Dimensionless interface temperature variations (Eq. 4.11) as a function of the nondimensional enthalpy of vaporization, ∆H * vap , and the bubble Peclet number, Pe b , for a 100 µm bubble with large vapor content Y 0 = 0.99 (top) and low vapor content Y 0 = 0.05 (bottom).

Figure 6 :

 6 Figure 6: Color online. Real and imaginary parts of 3Pe -1 M , where M (Pe) = √ Peı coth( √ Peı) -1, as a function of Pe.

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Color online. Phase speed predictions for systems with 100 µm bubbles with concentration β 0 = 10 -4 at 1 atmosphere. The ambient temperature is fitted to obtain a vapor molar fraction of 0.95. The dashed black line represents the predictions of the model accounting for mass transfer effects that imposes constant temperature at the interface. The solutions obtained from the full model are showed for different values of the enthalpy of vaporization which are given with respect to the enthalpy of vaporization of water, ∆H 0 vap , such that ∆H r = ∆Hvap ∆H 0 vap. For reference, the continuous black line represents the predictions of the model for a pure gas bubble, which neglects mass transfer effects across the interface.

Figure 10 :

 10 Figure10: Color online. Sensitivity factor (as defined in Eq. 6.13) as a function of the vapor content and frequency for a 100 µm gas/vapor bubble. The vapor content is changed by modifying the ambient temperature. As the amount of vapor increases, the bubble radial oscillation is more important for a constant nondimensional pressure wave amplitude.

  b , must be also smaller than one and the Peclet, Pe b =

  Mw whose linearized version is∆p I eq = ∆p b + ∆Y I . (8.1)Replacing ∆p I eq by the expression given by the linear Clausious-Clapeyron relation (Eq. 2. 3.16 to express the vapor fraction variations as a function of the mass flux across the interface, we can rewrite Eq. 8.1 as J * 0,eq ∆R = J * c,eq ∆H * vap ∆T I c -1 ∆p b ,

obtain the phase speed and attenuation.

The authors would like to acknowledge Claudine Boehm and Philippe Berthet for their support and discussions.


D Fuster

J M Conoir

T Colonius

On the effect of direct bubble-bubble interactions on linear wave propagation in bubbly liquids

Keywords: numbers: 43.35.-c, 43.30.-k, 43.25.Yw

We study the influence of bubble-bubble interactions on the propagation of linear acoustic waves in bubbly liquids. Using the full model proposed by Fuster and Colonius [J. Fluid Mech, Vol. 688, pp. 253-289, 2011], numerical simulations reveal that direct bubble-bubble interactions have an appreciable effect for frequencies above the natural resonance frequency of the average size bubble. Based on the new results, a modification of the classical wave propagation theory is proposed. The results obtained are in good agreement with previously reported experimental data where the classical linear theory systematically over predicts the effective attenuation and phase velocity.

INTRODUCTION

Acoustic wave propagation in bubbly liquids is relevant to diverse applications including cavitation reloading of surfaces exposed to underwater explosions, shockwave lithotripsy, and high-intensity focused ultrasound. Since Van Wijngaarden [1] proposed the first equations of motion for a dilute dispersion of gas bubbles in a liquid, several authors have proposed extensions to the original equations. Ainslie & Leighton [2] provide an account of this evolution of the theory of linear wave propagation in bubbly liquids. Most approaches are devoted to deriving models that account for the different types of dissipation present in the system, from the initial works of Chapman and Prosperetti [3,4] to most recent studies [5]. One of the strongest assumptions in these models is that mutual interactions among bubbles are negligible except through their effect on the mixture-averaged flow. For instance, Sangani [6] proposes correction terms for large void fractions that show that these terms tend to decrease the attenuation near the bubble resonance frequency. Neglecting bubble-bubble interactions formally restricts the applicability of the classical linear theory to low concentrations, where bubbles are so far from each other that the near field pressure around a bubble is not felt by the surrounding bubbles [7]. The averaged interspace distance between bubbles, d ij , is expected to scale as dij R0 ≈ Cβ -1/3 , where C is a constant of order unit and β is the void fraction. For void fractions as low as 10 -4 -10 -3 vol/vol the averaged interspace distance is only around 10 times the averaged bubble radius R 0 . In the simplest analysis, we would expect the response of each bubble to be influenced by the instantaneous response of the surrounding bubbles for higher void fractions.

It is possible to find in the literature various experimental works measuring the phase velocity and attenuation of bubbly liquids [8][9][10][11][12]. In most of these works, the comparison of the classical linear theory and the experimental results is not satisfactory for frequencies near resonance when comparing both, attenuation and phase velocity. In general, the experimental investigation of direct bubble-bubble interactions is difficult due to the presence of reverberation [13,14] and the dearth of diagnostics capable of directly measuring the the void fraction and bubble radius distribution [15]. As an alternative one can resort to numerical models where the value of all these variables is given. In this case, the main problem is that it is not computationally tractable, for large numbers of bubbles, to directly resolve both phases. Various volume and phase-averaged continuum approaches have been proposed in the literature. In the most sophisticated of these models, two-way coupled (dynamic) effects of cavitation are considered meaning that the time evolution of spherical bubbles are determined by solving a Rayleigh-Plesset-type (RP) equation whose driving pressure is determined from the local pressure of the continuous phase, and whose result effects the continuous phase by altering the void fraction [1,7]. Various extensions of the RP model adequately account for phase change [16], liquid compressibility, heat and mass transfer within the bubble contents [17]. However, existing phase and volume-averaging approaches require two fundamental assumptions that limit their applicability: low void fraction, and scale separation between the typical bubble size and the length scales associated with the (averaged) flow field. In this paper, we consider a new bubbly flow methodology based on volume-averaged equations that relaxes the scale-separation assumption and directly accounts for bubble-bubble interactions [18]. The numerical results are used here to gain insight into the effect of direct bubble-bubble interactions and to propose correction terms in the classical theory for linear wave propagation in bubbly liquids in order to improve the accuracy of the predictions for frequencies near the bubble natural resonance frequency.

MODEL AND NUMERICAL METHOD

The full equations for the model are detailed in Fuster & Colonius [18]. The model utilizes a volume-averaging approach [1]. We solve for the averaged continuity, momentum and energy equations in addition to the advection equation for every component present in the system. These equations can be written as [19]

∂ρ ∂t + ∇ • (ρu) = 0 ( 1
)
∂ρu ∂t + ∇ • (ρuu + pI) = 0 (2) ∂ρe ∂t + ∇ • (u(ρe + p)) = 0 ( 3
)
where ρ is the averaged density, t is time, u is the averaged mixture velocity, p is the pressure and e is the specific energy. This set of equations is closed by specifying an equation of state. In this case we assume that this equation obeys the following form

pΓ + Π ∞ = ρe - 1 2 ρ||u|| 2 (4)
In this work we consider a single bulk component (i.e. water) and a disperse phase (e.g. air bubbles) with void fraction β. This model differs from previous models in how the void fraction, β, is defined in terms of the bubble radius, and how the average pressure in a computational cell relates to the pressure felt by the corresponding bubbles in that cell. Bubbles are treated as point particles that are advected with the flow. If a polydisperse mixture is desired, then we set the initial radii of each bubble according to a representative probability density function. To compute β, we multiply each individual bubble volume by discrete delta function that smears its value onto the computational grid, and, at each point in the domain, sum this over all bubbles. To treat the dynamics of each bubble, the model introduces a new RP-like equation that accounts for a cluster of neighboring bubbles (i.e. within a computational cell). This equation is a generalization of the incompressible bubble-cluster model of Ilinskii et al. [20] to compressible liquids, and likewise assumes spherical bubbles and potential flow in the vicinity of the bubbles.

The key step in the derivation is to represent the pressure at infinity as seen by the bubbles in the potential flow model in terms of the average pressure of the liquid within a computational cell. From the potential flow solution for the bubble cluster, we may expand the pressure at infinity in terms of the cell pressure plus a correction term that can be estimated based on the bubble dynamics. The correction term thus represents a modification to the local value of pressure associated with the bubble dynamics. If the void fraction is sufficiently small, this term also becomes small and the equations relax to the classical volume-averaged equations. By taking a cell based approach rather than a multiple-scales approach, the model we obtain is grid dependent, much in the same way a Large Eddy Simulation of turbulence depends on the grid spacing. In the limit of the grid spacing going to zero, the LES model switches off, resulting in a direct numerical simulation of the turbulent flow. In the present model, as the grid is refined to the size of an individual bubble, the model exactly represents the dynamics of an isolated spherical bubble. This is verified in detail by comparing against theoretical solutions for this case [18]. At the other extreme, the model is demonstrated to converge to the solution to ensemble-averaged models for low spatial resolution. The computational expense is greatly reduced compared to the corresponding ensemble-averaged approach. Therefore, the current model seems suitable to numerically investigate complex problems where the computational effort of both direct interface capturing and current ensemble-averaged models are prohibitive. This cavitation model has been incorporated into a multi-fluid, compressible flow solver that utilizes advanced WENO shock-and interface-capturing techniques [19]. Details about the implementation of the model on the fluid solved can be found in Ref. [START_REF] Colonius | Proc. 8th Int. Symp. Cav[END_REF].

LINEAR MODELS

Classical linear model

The nonlinear model presented above places no restriction on the amplitude of waves propagating through, and interacting with, the bubbles. However for practical proposes it is interesting to derive linearized solutions for wave propagation in bubbly liquids that can be compared with the solutions of the full nonlinear model.

In the classical theory, it is possible to derive a wave equation for the effective pressure far from any bubble [22]

∇ 2 P ∞ + k 2 m P ∞ = 0, (5)
where P ∞ stands for the amplitude of the effective pressure wave acting far from the bubbles. The value of the effective wavenumber in the mixture, k m , can be obtained as

k 2 m = ω 2 c 2 + 4πω 2 ∞ 0 af (a)da ω 2 0 -ω 2 + 2ıδω , (6)
where ı = √ -1, ω is the angular frequency, c is the speed of sound in the liquid, f (a) is the bubble distribution function and ω 0 and δ are usually identified with the bubble resonance frequency and damping constant obtained from the linearization of the Rayleigh-Plesset equation,

R R 1 - Ṙ c + 3 2 Ṙ2 1 - Ṙ 3c = S. (7)
In this equation S stands for the driving term that can be written as

S = 1 + Ṙ c + R c ∂ ∂t H i + ∂Φ ∞ ∂t + I (8)
where H i the liquid enthalpy at the bubble interface

H i = pi p0 dp ρ
, Φ ∞ represents the potential induced by the background flow field, ∂Φ∞ ∂t = -p∞ p0 dp ρ , and I is the potential from the N surrounding bubbles impacting at the location of the ith bubble x i ,

I = N j =i ∂Φ j (x i , t) ∂t . (9)
For weak perturbations we can simplify the driving term as

S = H i + ∂Φ ∞ ∂t + N j =i ∂Φ j (x i , t) ∂t (10)
which reduces to the expression for the classical driving term when direct bubble-bubble interactions are neglected, i.e.

lim β→0 S = p i -p 0 -P ∞ sin(ωt) ρ . (11)
Thus, imposing a bubble radius oscillations of the form

R i = R i,0 (1 + X 0 e ıωt ), (12)
where X 0 ≪ 1 and R i,0 is the equilibrium bubble radius at the reference pressure p 0 , the expressions for the bubble resonant frequency and damping coefficients are [START_REF] Commander | [END_REF]]

ω 2 0 = p 0 ρR 2 i,0 Re(Υ) - 2σ R i,0 p 0 , (13)
δ = 2µ ρR 2 i,0 + p 0 2ωρR 2 i,0 ℑ(Υ) + ω 2 R i,0 2c , ( 14
)
where σ is the surface tension, ρ and µ are the liquid density and viscosity, and Υ is a complex valued function that relates the pressure and radius changes in the linear regime. This function is obtained from the solution of the conservation equations inside the bubble under proper boundary conditions [23,24]. Using the specific heat ratio γ and the Peclet number Pe = ωR 2 i,0 /D g , D g being the gas thermal diffusivity, the function Υ can be written as

Υ = 3γ 1 -3(γ -1)ıPe -1 √ ıPe coth √ ıPe -1 . (15)
The attenuation and phase velocity of the wave can then be obtained as a function of the complex and imaginary part of the effective velocity [START_REF] Commander | [END_REF].

Extended linear model accounting for direct bubble-bubble interactions

The linear solutions obtained above only consider the interaction of bubbles present in the system through the averaged field. The full solution can be obtained applying the procedure proposed by Keller and Miksis [25] for a single bubble to the entire bubble cluster. In this case, we consider that the background potential Φ ∞ is known and we need to solve for the potentials emitted by each bubble in the cloud. When the incident field is a planar wave propagating along x with angular frequency ω and pressure amplitude P * , then,

Φ ∞ = -(P * /ρω) cos(ω(t -x/c)). (16)
In this situation it is readily shown that in the linear regime, except for frequencies much higher than the natural resonance frequency, the driving term can be written as

S = p i -p 0 -P * sin(ωt) ρ + I, (17)
where P * is constant and I accounts for all the bubbles present in the cluster. Given the number of bubbles present in the system N that interact directly with a given bubble, it is possible to write N Bernoulli equations at the interface of each bubble present in a given volume,

N j=1 ∂Φ j (x i , t) ∂t = - 1 2 Ṙ2 i -H i - ∂Φ ∞ ∂t . (18)
Thus, we obtain a system of N equations that can be solved for the potentials Φ i given the background potential variations ∂Φ∞ ∂t . The N potentials need to be evaluated at every bubble location. Assuming that the bubble is a monopolar source emitting a spherical wave decaying with the inverse of the distance between the ith and jth bubble, d ij , the potential can be written as a function of the retarded time τ

= dij c as Φ j (x i , t) = 1 d ij f (t -τ ) = Φ j (t -τ, x j ) R j (t -τ ) d ij .
(19) This concept has been used by Zeravcic [26] to investigate the collective response of bubble clusters. In this work this solution will be obtained using 3D simulations obtained with the nonlinear model presented above in the limit of small perturbations.

The reasoning above allows to obtain numerical solutions for the full problem but it is computationally intensive when N is large. Thus, we derive here correction terms for the classical linear theory that capture the influence of the scattered waves emitted by each individual bubble. To this end, it is interesting to compare the equation used by the classical linear theory (Eqs. 5 and 11) with that obtained from the solution of the individual potentials (Eqs. 16 and 17). As one can clearly see, part of the potential that is emitted by the bubbles is already accounted by considering a coherent planar wave propagating across the bubbly media that is dissipated with a phase velocity influenced by the effective medium. Thus, if one desires to introduce a correction to account for direct bubble-bubble interactions in the classical linear theory we have to look for a correction term of the type

S = p i -p 0 -P ∞ sin(ωt) ρ + I * , (20)
where now I * represents the part of the potentials emitted by the surrounding bubbles not captured by the averaged field equations. The procedure to obtain the exact value of I * is as computational demanding as to solve the full problem (Eq. 17) and therefore it is not interesting for practical purposes. As an alternative, we can obtain an estimation of I * modelling the potential emitted from the jth bubble using the total potential Φ at a given location,

Φ ≡ Φ ∞ + Φ j + N k =j Φ k . ( 21
)
The difference between the total potential and the background potential at the interface of the jth bubble has to satisfy the Bernoulli equation at the interface of the jth bubble

∂ ∂t (Φ(R j ) -Φ ∞ (R j )) = - 1 2 Ṙj 2 - p j -p ∞ ρ . (22)
As previously stated, only part of this potential can be directly attributed to the jth bubble,

Φ(R j ) -Φ ∞ (R j ) = Φ j (R j ) + N k =j Φ k (R j ). (23)
In fact, using the equation above it is possible to rewrite Eq. 22 as an equation for the potential of the jth bubble as

∂Φ j (R j ) ∂t = g - 1 2 Ṙj 2 - p j -p ∞ ρ ( 24
)
where g is defined as

g ≡ ∂ ∂t Φ j (R j ) ∂ ∂t Φ j (R j ) + N k =j Φ k (R j ) . ( 25
)
This function, bounded between zero and one, is unknown for a general case. In the limiting case where the potentials emitted by the surrounding bubbles are more important than the potential generated by the bubble itself, the function becomes zero indicating that direct bubble-bubble interactions are not relevant. On the contrary, when the sum of the potentials emitted by the surrounding bubbles is negligible, g becomes unity, e.g. there is no need to solve for the system of N equations to obtain the potential emitted by each bubble because the coupling effect is negligible.

The relative importance of the potentials emitted by the surrounding bubbles can be characterized in terms of the attenuation induced by the effective medium surrounding the bubble. To compute the total interaction potential impacting the ith bubble we integrate the contribution of all potentials emitted at a distance smaller than L = CR 0 /β 1/3 0 , where C is a free parameter to be determined later. At distances larger than L, the bubbles are assumed to interact through the averaged field. In absence of attenuation we can write,

N j=1 ∂Φ j (x i , t) ∂t ≈ L 0 R 0 r ∂Φ j (x j , t) ∂t 4πr 2 ndr ( 26
)
where n is the number of bubbles per unit volume that can be obtained as,

n = β 0 4/3πR 3 0 . ( 27
)
Eq. 26 is also expected to be applicable in a general case for small distances, where the attenuation of the effective medium is negligible. However, if L takes large values, the effective medium is going to induce an attenuation rate A for the effective contribution of bubbles far away from the ith bubble. If A is known, the effective contribution of the potentials is corrected

N j=1 ∂Φ j (x i , t) ∂t ≈ L 0 R 0 r ∂Φ j (x j , t) ∂t 4πr 2 ne -Ar dr, (28)
and the direct bubble-bubble interaction term correction can then be obtained as

I * = N j =i ∂Φ j (x i , t) ∂t = N j=1 ∂Φ j (x i , t) ∂t - ∂Φ i (x i , t) ∂t . (29)
To illustrate the influence of bubble-bubble interactions, let's now consider an ideal situation where all the bubbles emit the same potential intensity at a given instant (e.g. all the bubbles have the same pressure and expansion velocity at a given instant). Thus, we write Eq. 24 as

∂Φ j (R j ) ∂t = gP 0 . (30)
In this case, for distances smaller than Lω/c < 1, Eq. 28 can be exactly evaluated as

N j=1 ∂Φ j (x i , t) ∂t = 3(e AL -AL -1)β 0 (AR 0 ) 2 e AL + 3(e AL -AL -1)β 0 P 0 . ( 31
)
When L is large (AL ≫ 1), we can approximate the sum of all the potentials emitted by the bubbles at a given location,

N j=1 ∂Φ j (x i , t) ∂t ≈ 3β 0 (R 0 A) 2 + 3β 0 P 0 , (32)
and the interaction correction term can be obtained using Eq. 29 as

I * = -(R 0 A) 2 (R 0 A) 2 + 3β 0 P 0 . ( 33
)
This equation reveals that the interaction term neglected in the classical theory can be expressed in terms of the effective attenuation of the medium. If the effective attenuation at a given frequency is below a critical attenuation A c value,

A c = √ 3β 0 R 0 , (34)
then I * tends to zero when decreasing the attenuation.

For dilute systems, where the attenuation is proportional to the void fraction for frequencies well below the bubble resonance frequency [START_REF] Commander | [END_REF], the result above implies that the term neglected in the classical approach is of order O(β 0 ). When the attenuation is large, the magnitude of the total interaction potential is maximized (|I * | → P 0 ) tending to oppose a significant resistance to the expansion/compression of the bubbles present in the system. In this case it is important to note that, as the effective medium induces a large attenuation on the potentials emitted by the surrounding bubbles, the contribution of bubbles further than few times the averaged interbubble distance is going to be negligible. Thus, Eq. 26 is more suitable to obtain the interaction potentials at small distances. In this case the potential emitted by the jth bubble impacting on the ith bubble at distance r can be obtained introducing the time-lag, τ = dij c , R j = R 0 (1 + X 0 e ıωt e -ıωτ ) as,

∂Φ j (x i , t) ∂t = g ρ p ge Υ - 2σ R 0 + 4µıω X 0 e -ıωτ R 0 r .
(35) For monodisperse mixtures we replace Eq. 35 into Eq. 26 to integrate it from zero to a given distance L. Beyond this distance bubbles are assumed to interact through the averaged field. Thus we obtain

I * = I 0 X 0 e ıωt ( 36
)
where I 0 is zero in the regions of low attenuation (A < A c ) and it becomes

I 0 =- 1 ρ 3β 0 (R 0 ω/c) 2 1 -e -ıLω/c (1 + ıLω/c) p ge Υ - 2σ R 0 + 4µıω . ( 37
)
above the critical attenuation (Eq. 34). Following the same procedure used in the classical linear theory, the modification introduced in the driving term of the Rayleigh-Plesset equation due to direct bubble-bubble interactions gives an additional correction on the bubble's resonance frequency and damping coefficient as

ω 2 0 = p 0 ρR 2 i,0 Re(Υ) - 2σ R i,0 p 0 + Re(I 0 ) R 2 i,0 , (38)
δ = 2µ ρR 2 i,0 + p 0 2ωρR 2 i,0 ℑ(Υ) + ω 2 R i,0 2c + ℑ(I 0 ) 2ωR 2 i,0 . (39)
As stated above, this expression is not closed. It contains an unknown free parameter C and a frequency dependent function g that is bounded between zero (for regimes with small attenuation) and one (for regimes with large attenuation). We note that C can be identified with a non-dimensional interaction distance and therefore we expect it to be approximately constant irrespective of the bubble radius and void fraction C = Lβ 1/3 R0 . Numerical tests, which solve for the full problem, are used in the next section to obtain a better understanding about reasonable values of C and g for a general case.

In the next section we provide numerical and experimental evidence that supports the veracity of the arguments exposed in this section.

NUMERICAL TESTS

We now describe the general simulation setup that is used to solve the full (nonlinear) bubble model. We consider air bubbles in water with a known concentration β 0 and an average radius R 0 defined at a reference pressure, p 0 = 1 atm, and temperature, T 0 = 25 C. The bubbles are seeded according to a known probably distribution function, f (a), that represents the distribution of equilibrium bubble radii in the system. For instance, for the polydisperse case included in this manuscript, we initialize the bubble population with a random distribution of bubbles obeying the following Gaussian distribution

f (a) = 1 b √ 2π exp -(a-R0) 2 2b 2 0.1R 0 ≤ a ≤ 5R 0 0 otherwise (40)
where R 0 = 110µm and b/R 0 = 0.25 (Figure 1). Bubbles are distributed randomly in a rectangular domain of dimensions [-L x /R 0 :

L x /R 0 ] × [-L y /R 0 : L y /R 0 ] × [-L z /R 0 : L z /R 0 ].
A 1D Gaussian pressure pulse is initialized in a bubbly liquid at rest (u = 0) following the equation

p(x) = p 0 1 + ǫ exp -x 2 /s 2 p . (41)
The initial bubble radius is then obtained from Laplace's equation so that bubbles are in equilibrium with the initial pressure. This unphysical initialization procedure is only an artifact that only influences the initial state of bubbles in a small region with a characteristic length of the order of s p , and has no significant impact on the results. For the simulation contained in this work, we have set ǫ = 10 -4 and sp R0 = 1.

To obtain a representative behavior of the averaged bubble propagation, the pressure disturbance is integrated in a plane perpendicular to the wave propagation. Thus, we define the pressure perturbation p ′ at a giving sampling point x i as

p ′ (x i ) = (p(x i , y, z) -p 0 )dydz dydz . ( 42
)
Then, following Ref. [27], the attenuation and speed of sound are obtained using the Fourier Transform of the temporal evolution of p ′ (x i ) at multiple sampling points.

The values for the complex and real part of the complex wavenumber are fitted by least squares using the data from five equally spaced sampling points. In particular, two different sets are used; the first ones use a spacing of 10R 0 chosen to correctly capture the speed of sound in the high frequency range. The second set of points is chosen with spacing 30R 0 to minimize the influence of round-off errors in the low attenuation regions.

Validation test case: Wave propagation results neglecting direct bubble-bubble interactions

One of the advantages of the numerical approach is that it allows us to reproduce ideal situations where the code can be rigorously validated against traditional linear theories. In particular, in this section we will switch off the direct bubble-bubble interaction term and use a coarse mesh in the z direction so that the concentration in every computational cell is given by the averaged concentration. In these conditions, the full model presentented in Ref. [18] is expected to recover the classical linear theory predictions for the attenuation and phase speed of linear waves propagating across a bubbly liquid neglecting direct bubble-bubble interactions. As expected Figure 2 shows good agreement between the theoretical predictions and the simulation results, validating the capability of the code to reproduce the theoretical results predicted by the classical linear theory, in those conditions where the grid does not resolve the scattered waves emitted by the individual bubbles. Both attenuation and phase velocity are very well reproduced for frequencies up to the resonance frequency. The peak on the phase velocity is especially challenging for the numerical method due to the large velocities obtained there, but the match for the attenuation for values above 0.2 dB/cm is remarkably good. The results obtained become noisy as we approach to the maximum frequency that can be numerically solved, which is defined as f max = c/∆x. For the simulations included in this work we set ∆x = 10R 0 . When comparing the maximum frequency with the bubble natural frequency (f N = 1 2π 3γp0 ρR 2 0 ) we conclude FIG. 2. Speed of sound and attenuation of a polydisperse bubble cloud of randomly distributed bubbles obtained from simulation results and comparison with the traditional linear theory neglecting bubble-bubble interactions [START_REF] Commander | [END_REF].

The computational domain size is set to

that, for air/water systems at atmospheric pressure, we expect the results to be highly influenced by spatial descretization errors for nondimensional frequencies of the order of ten times the bubble resonance frequency ( f fN ≈ 10).

Numerical and theoretical results of wave propagation accounting for direct bubble-bubble interactions

In order to investigate the influence of direct bubblebubble interactions on wave propagation we consider various monodisperse bubble clouds with a different concentration β 0 and bubble radius R 0 . Unlike the previous case, the grid now is refined in order to capture not only the coherent wave, but also the waves scattered by the surrounding bubbles that are a result the bubble radius remains small compared to the grid size and the interspacing between bubbles is larger than the grid size. The grid spacing quoted above and used in the calculations that follow was selected by examining results on a series of progressively finer meshes and verifying an approximate independence of the phase speed and attenuation to grid resolution. As was reported in [18], this establishes that we are adequately capturing the local pressure and velocity fluctuations induced by the bubbles. As the leading-order truncation error in our WENO scheme is dissipative, we also computed attenuation rates for waves propagating on the finest grid in the absence of bubbles, i.e a pure inviscid liquid. As numerical dissipation is the only source of attenuation in this case, this serves to establish a lower bound the attenuation rates that we can reliably assess. The value is indicated on the attenuation plots that follow.

Figure 3 contains the phase speed and attenuation simulation data obtained from simulations for the conditions of Case A (Table I). Theoretical predictions for the classical theory, (LT g = 0) and the extended theory with g = 1 (LT g = 1) for different values of C are included for reference. Recall that we expect the real solution to be bounded between these two limiting cases, where results tend to the classical linear theory predictions for the regions with low attenuation (A ≪ A c = 8.11 dB/cm) whereas the predictions for g = 1 should better represent the regions of large attenuation (A > A c ). The solution in the limit of g = 0 are independent of the non-dimensional interaction distance C. On the contrary, when g > 0, results are sensitive to C. Based on this fact, we use the region with the largest attenuation to calibrate the non-dimensional interaction distance C. We remark that in this region we expect g to be approximately constant and equal to one. The peak values of both, phase speed and attenuation are over For reference, the predictions of the classical linear theory (g = 0) and the results from simulations are included.

predicted by the linear theory up to factor 2 for Case A. The modified linear theory significantly improves the fitting between theory and numerical results fitting the "non-dimensional interaction distance" to C = 2.3. This value falls within the range of the averaged interspacing, which is consistent with the fact that only bubbles nearby directly interact with a given bubble, whereas the rest interact through the averaged field.

Figure 4 depicts the results for the various conditions included in Table I, all of them with constant bubble number density. We can see how the differences between the two limiting cases become clear as the concentration increases for frequencies of the order of the bubble natural resonance frequency. As expected, direct bubble-bubble interactions become less important as void fraction decreases. The maximum errors are significantly reduced down to 20 % for concentrations of the order of 10 -5 and the differences are almost negligible for lower concentrations. The critical concentration threshold for which direct bubble bubble interactions are important therefore coincides which the expected value (β = 10 -4 ) although one must keep in mind that we only observe significant differences for frequencies of the order of the resonance frequency. Finally, the distance required to obtain the best fitting value increases as the void fraction decreases (Table I). This is consistent with the fact that as concentration decreases, the bubble inter-spacing increases. A constant value of C = 2.3 in the three cases significantly improves the fitting between the extended linear theory and numerical simulations in the range of frequencies where the attenuation is above the critical value A c .

To gain more insight into the applicability of the fitting constant value to a wider range of conditions we represent in Figure 5 the nondimensional fitting parameter C = Lβ 1/3 R0 for various tests with slightly higher and lower values of the bubble number density. Remarkably, the value of C remains approximately constant ranging from 2 to 2.5 for all the different conditions tested, which include a very large range of bubble concentrations.

It is interesting to note that some of the observations extracted from Figure 4 have been systematically observed in previous experimental and theoretical studies reported in the literature. For instance, Silberman [9] has observed a systematic over prediction of attenuation predicted by the traditional linear theory on the regions of large attenuation above the resonance frequency, slightly under predicting the attenuation coefficient for frequencies below the attenuation peak. Similar mismatching has been found by Commander and Prosperetti in Ref. [START_REF] Commander | [END_REF] when fitting the data from both, Silberman and Fox et al [8]. This later case is interesting because measurements also include phase speeds experimentally measured in the frequency range where there is a clear mismatch in the attenuation. The theoretical results derived by Sangani also pointed that the attenuation around the resonance frequency is lower than that predicted by the classical theory when introducing correction terms to account for the effect of non-negligible void fractions [6]. Wilson et al [11] and Leroy et al [12] report lower values of the phase velocities than those predicted by the traditional linear theory in the range of frequencies where we find large attenuation. These observations are consistent with the numerical results obtained in this work and with the fact that direct bubble-bubble interactions are probably negligible in the low frequency limit (e.g. g ≈ 0), whereas they have an important influence for frequencies above the resonance frequency.

In an attempt to reproduce this effect observed experimentally, we numerically reproduce the lowest concentration tested by Silberman [9] (β 0 = 3.77 • 10 -4 ), where the bubble size ranges between 0.994 mm and 1.07 mm. These conditions are in the limit of the highest values of void fraction for which it is possible to obtain a grid independent solution with the current numerical model. In Figure 6 we include the theoretical predictions for the two limiting cases corresponding to the classical linear theory (g = 0) and the limit of g = 1, for which we use the value of L previously fitted from numerical results (Lβ 1/3 /R 0 ≈ 2.3)). The agreement between the modified linear theory, numerical simulations and experiments is relatively good in the regions of large attenuation, where the classical linear theory over predicts the attenuation. At the lowest frequencies, where the attenuation is small, the experimental values tend to the g = 0 solution, consistent with the previous discussion. We note that the simulation results are not reliable in this regime because the attenuation values are below our estimate for the numerical dissipation, even on the finest mesh we considered. While finer meshes (or less dissipative schemes) could be used to alleviate this problem, we refrain from such calculations since they would not affect the conclusions drawn.

CONCLUSIONS

This work uses the model proposed by Fuster and Colonius [18] to study the influence of direct bubblebubble interactions on the acoustic properties of the effective medium. Taking advantage of the results obtained with the full nonlinear model that accounts for bubble-bubble interactions, we are able to isolate the influence in the linear regime of waves scattered by individual bubbles on the effective linear acoustic response of the mixture. Consistent with previously reported experimental [9] and theoretical data [6], the numerical results provide evidence that the classical linear theory overpredicts the attenuation and phase velocity above the natural resonance frequency.

A modification of the linear theory is suggested to try to capture these effects. The new theory introduces

LT g=0 LT g=1

FIG. 6. (Color online) Phase velocity and attenuation of a polydisperse bubble cloud of randomly distributed bubbles obtained from simulation results and comparison with the experimental results of Silberman [9], the linear theory predictions of ref [START_REF] Commander | [END_REF] (LT g = 0) and the modified linear theory accounting for direct bubble bubble interactions in the limiting case of g = 1. The dashed horizontal line on the attenuation plots indicates the upper bound for which the numerical attenuation values are reliable, due to numerical dissipation.

a characteristic distance L that we identify with an "interaction distance". In those regions where the attenuation is larger than a critical value defined as A c = √ 3β 0 /R 0 , the individual bubble response is assumed to be influenced by the waves scattered by bubbles at distances smaller than L, whereas the interaction with bubbles placed further takes place through the averaged pressure field. The comparison between the numerical results obtained with the full non-linear model and the extended linear theory reveals that the proposed correction term significantly improves the theoretical predictions when the nondimensional value for the interaction distance, Lβ 1/3 0 /R 0 , is approximately constant (Lβ 1/3 0 /R 0 ≈ 2.3) for void fractions below 10 -3 .

FIG. 1 .

 1 FIG.1. Histogram of the bubble size distribution used in the simulation for validation against the classical linear theory results (results presented in Figure2). R0 = 110µm.

  : 3500] × [-250 : 250] × [-250 : 250], the grid size is ∆x/R 0 × ∆y/R 0 × ∆z/R 0 = 10 × 10 × 500, and the initial void fraction, β 0 = 4 3 πn a 3 f (a)da is set to β 0 = 2.5 10 -4 .

FIG. 3 .

 3 FIG. 3. (Color online) Influence of the fitting parameter L/R0on the phase velocity and attenuation of a monodisperse bubble cloud of randomly distributed air bubbles in water obtained from simulation results for case A (g = 1 in all cases). For reference, the predictions of the classical linear theory (g = 0) and the results from simulations are included.

FIG. 4 .

 4 FIG.4. Phase velocity and attenuation of a monodisperse bubble cloud of randomly distributed air bubbles in water obtained from simulation results for case A (top), B (middle) and C (bottom). The dashed horizontal line on the attenuation plots indicates the upper bound for which the numerical attenuation values are reliable, due to numerical dissipation. For reference, the predictions of the classical linear theory and the the proposed extension of the theoretical results with constant C = 2.3 are included.

FIG. 5 .

 5 FIG. 5. Best fitting value of the nondimensional interaction distance, Lβ 1/3 /R0, for various bubble radius and concentrations for monodisperse bubble clouds. The number of bubbles range from n/n ref = 0.5 to 2. The nondimensional value C = Lβ 1/3 R 0 remains approximately constant irrespective of the concentration and bubble radius.

TABLE I .

 I Bubble radius, concentration and best fitting interaction distance L for the simulation results shown in Figure 4. [-250 : 250], is discretized with a regular Cartesian grid of size ∆x/R ref = ∆y/R ref = ∆z/R ref = 10, where R ref = 30µm. The bubble number density is initially set to a constant value given by n ref = β ref 4/3πR 3 ref where β ref = 2 10 -4 . Monodisperse bubble clouds with R 0 < R ref at constant n ref fulfill two different conditions:

		R0 (µm) fres (kHz)	β0	L/R0
	Case A	30	109	2 10 -4	40
	Case B	15	219	2.5 10 -5	63
	Case C	7.5	438	3.125 10 -6 133
	of direct bubble-bubble interaction. The simulation domain, Lx R ref × Ly R ref × Lz R ref ∈ [-3500 : 3500] × [-250 : 250] ×

This work was supported by the French National Agency for Research (ANR-113S09-007-03) project DiA-MAN.