Olivier Louisnard

Analytical expressions for primary Bjerknes force on inertial cavitation bubbles

Keywords: numbers: 47, 55, dd, 43, 35, Ei

The primary Bjerknes force is responsible for the quick translational motion of radially oscillating bubbles in a sound field. The problem is classical in the case of small-amplitude oscillations, for which an analytical expression of the force can be easily obtained, and predicts attraction of subresonant bubbles by pressure antinodes. But for high-amplitude sound fields, the bubbles undergo large amplitude nonlinear oscillations, so that no analytical expression of the force is available in this case. The bubble dynamics is approximated on physical grounds, following the method of Hilgenfeldt et al. [J. Fluid Mech., 365, 171 (1998)], but carefully accounting for surface tension. The analytical expression of the maximum radius of the bubble is recovered, the time of maximum expansion is noticeably refined, and an estimation of the collapse-time is found. An analytical expression for the time-varying bubble volume is deduced, and the Bjerknes force is obtained in closed form. The result is valid for any shape of the sound field, including purely standing or purely traveling waves, and is ready to use in a theoretical model of bubble clouds evolution. Besides, the well-known sign inversion of the Bjerknes force for large standing waves is recovered and the inversion threshold in the parameter space is obtained analytically. The results are in good agreement with numerical simulation and allow a quantitative assessment of the physical parameters effect. It is found that either reducing surface tension, or increasing the static pressure, should produce a widening of the bubble-free region near high-amplitude pressure antinodes.

I. INTRODUCTION

When excited by a sinusoidal sound field, gas bubbles undergo radial oscillations. Most of the practical applications of this phenomenon, known as acoustic cavitation, use high-amplitude sound fields, of typical amplitude greater than the static pressure, so that the liquid is under tension for some part of the cycle. In such conditions, whatever the frequency, two distinct dynamic bubble behaviors can be clearly divided by the so-called Blake threshold [1][2][3][4]: in the tension phase, very small bubbles are retained to grow by surface tension. Conversely, larger ones suffer an explosive expansion followed by a violent collapse, responsible for chemical [START_REF] Suslick | proceedings of the NATO Advanced Study Institute on Sonoluminescence and Sonoluminescence[END_REF], mechanical effects [START_REF] Lauterborn | [END_REF]7] and sonoluminescence [8][9][10]. The latter oscillation regime is known as "inertial cavitation".

Bubbles in liquids experience various hydrodynamic forces. The buoyancy force is the most familiar one, and is the pressure force that an sphere of liquid replacing the bubble would experience. This remains true in a accelerating liquid [START_REF] Magnaudet | ASME Fluids Engineering Division Summer Meeting[END_REF], and the generalized buoyancy force experienced by the bubble is -V ∇P where P(r, t) is the pressure that would exist at the center of the bubble if it were absent, and V (t) the bubble volume. For a bubble oscillating radially in a sound field, both P(r, t) and V (t) are oscillatory quantities so that the time-average of the product over one cycle is not zero. The bubbles experiences therefore a net force known as "primary Bjerknes * Electronic address: louisnar@enstimac.fr force" [START_REF] Bjerknes | Fields of Force[END_REF][START_REF] Blake | [END_REF]:

F B = -V ∇P (1)
The Bjerknes force can be easily calculated from the knowledge of both the shape of the sound field and the bubble dynamics. A classical result is that for lowamplitude standing waves, sub-resonant bubbles are attracted by pressure antinodes, while bubbles larger than resonant size are repelled [14][15][16]. For the case of strong driving pressures, sub-resonant inertial bubbles can also be attracted by pressure antinodes, which constitutes the basic principle of SBSL levitation cells [8,17]. However it has been shown by numerical calculations that above a given threshold, the primary Bjerknes force on sub-resonant inertial bubbles undergoes a sign change [18]. This behavior is due to the resonance-like response curve (termed as "giant resonance" by Lauterborn and co-workers [START_REF] Lauterborn | proceedings of the NATO Advanced Study Institute on Sonoluminescence and Sonoluminescence[END_REF]) of the bubble just above the Blake threshold, which is a physical consequence of the effect of surface tension. Experiments indeed demonstrate that above a certain driving level, no bubbles are visible in the neighborhood of large pressure antinodes [START_REF] Mettin | [END_REF].

Quantitative agreement between theory and experiment has been found in the case of linear or quasi-linear oscillations [15]. Particle simulations [START_REF] Mettin | [END_REF]21] were also found in excellent agreement with recent experiments involving inertial bubbles [START_REF] Koch | Proceedings of the IEEE International Ultrasonics Symposium[END_REF]. While the Bjerknes force can be calculated analytically for linear bubble oscillations, only numerical results can yet be found for inertial bubbles [18,[START_REF] Mettin | Oscillations, Waves and Interactions[END_REF]. An analytical expression for the latter would first be helpful in particle or continuum models, describing the self-organization of bubbles, in order to get more efficient calculations. Furthermore, analytical results allow a direct assessment of the sensitivity of the force to the physical parameters, and the establishment of scaling laws. These two objectives motivated this study.

Owing to the strong nonlinearity of the bubble dynamics equations, inertial cavitation has long been thought intractable analytically, up to the seminal papers of Löfstedt et al. [START_REF] Löfstedt | [END_REF] and Hilgenfeldt et al. [1], who demonstrated that several terms of the Rayleigh-Plesset equation (RP) could be neglected during the explosive expansion of the bubble. This theoretical breakthrough allowed to obtain scaling laws for the maximum radius of the bubble and the time of maximum expansion. In this paper, we closely follow the approach of Hilgenfeldt et al. [1] and refine their analytical solutions in order to account more precisely for the effect of surface tension. The approximate dynamics found are then used to obtain an analytical expression of the bubble volume. The latter are then conveniently recast in order to obtain the Bjerknes force (1) in closed form, in any acoustic field, including the two extreme cases of traveling and standing waves. Finally, in the latter case, we seek an approximate expression of the Bjerknes force inversion threshold, evidencing the role of surface tension.

II. PRIMARY BJERKNES FORCE

A. Acoustic field

We assume that the acoustic field in the liquid is monoharmonic at angular frequency ω, and defined in any point r by P(r, t) = P (r) cos [ωt + φ(r)] .

(

) 2
This expression may represent a traveling wave, a standing wave, or any combination of both. We also define the pressure gradient in general form as

∂P ∂x i (r, t) = G i (r) cos [ωt + ψ i (r)] , (3)
where the fields G i and ψ i can be expressed as functions of P and φ once the acoustic field is known. The following two extreme cases deserve special consideration:

• for a standing wave, φ(r) = φ 0 , so that G i (r) = ∂P/∂x i and ψ i (r) = φ 0 ,

• for a traveling wave, P (r) = P 0 and φ(r) = -k.r so that G i (r) = k i P 0 and ψ i (r) = φ(r) -π/2.

B. Bubble model

The radial oscillations of a gas bubble in a liquid under the action of the sound field can be described by the Rayleigh-Plesset (RP) equation [1,[25][26][27]:

R R + 3 2 Ṙ2 = 1 ρ p g + R c l dp g dt -4µ Ṙ R - 2σ R -(p 0 + P(t)) , (4)
where p 0 is the hydrostatic pressure, p g (t) is the gas pressure, ρ, µ and c l are the density, viscosity and sound speed of the liquid, respectively, and σ is the surface tension. The ambient radius of the bubble R 0 is the radius that would have the gas the in absence of the sound field. Time is non-dimensionalized by the angular frequency ω, and in order to obtain a formulation consistent with Ref. 1, we set

p 0 + P(r, t) = p 0 (1 -p cos x), (5)
so that

p = P (r)/p 0 (6) x = ωt + φ(r) -π. (7)
Using x as the time-variable, and non-dimensionalizing pressure with p 0 , equation ( 4) can be written as:

RR ′′ + 3 2 R ′2 = R 2 res 3 p * g + Rω c l dp * g dx - 4µω p 0 R ′ R -α S R 0 R + p cos x -1 , (8)
where primed variables denotes d/dx,

R res = ω -1 (3p 0 /ρ) 1/2 (9)
is the resonance radius, and

α S = 2σ/p 0 R 0 (10)
is the dimensionless Laplace tension. Several models can be used for the bubble internal pressure p g , [10,[28][29][START_REF] Prosperetti | proceedings of the NATO Advanced Study Institute on Sonoluminescence and Sonoluminescence[END_REF][START_REF] Storey | [END_REF][32]. As will be seen below, we are mainly interested here in the expansion phase of the bubble, during which the density of the gas in the bubble remains weak, so that the precise choice of the thermal bubble interior's model is unimportant. However, in order to assess the validity of the approximate expressions developed hereafter, simulations will be performed by using the Keller equation [33,34]. The bubble interior is modeled by using a thermal diffusion layer following Ref. 32, neglecting water evaporation and condensation through the bubble wall. In the remaining part of the paper, we will consider air bubbles in water (σ = 0.072 N.m -1 , p 0 = 101300 Pa, ρ = 1000 kg.m -3 , c l = 1498 m.s -1 , µ = 10 -3 Pa.s).

C. The Bjerknes force

The primary Bjerknes force acting on a bubble is defined as

F B = -V (t)∇P , (11)
where V (t) is the instantaneous bubble volume. The average is taken over one acoustic period, so that, using Eq. ( 3):

F B i = -G i (r) 1 T T 0 V (t) cos [ωt + ψ i (r)] dt. (12)
Using the dimensionless time x defined by (7) and the periodicity of V , the latter expression becomes

F B i = G i (r) 1 2π 2π 0 V (x) cos [x -φ(r) + ψ i (r)] dx.
(13) The generic problem is therefore to obtain an approximate analytical expression for the integral

I = 1 2π 2π 0 V (x) cos (x -x 0 ) dx, (14)
valid for any bubble dynamics, and for any value of x 0 . The problem can be easily solved for small-amplitude linear oscillations [16]. Here, we focus on the case of inertial oscillations, that is for any combination of parameters (p, R 0 ) above the Blake threshold. The special cases of standing waves and traveling waves can be simply recovered by setting, respectively, x 0 = 0 and x 0 = π/2.

III. APPROXIMATE EXPRESSIONS A. Bubble radius

The method used to obtain analytical formula for the bubble radius are mainly inspired from the approach of Hilgenfeldt et al. [1]. For self-consistency, we will recall in this section the main lines of the method, and, where convenient, specify the refinements obtained by our approach.

Figure 1 displays the dimensionless bubble radius (1a), bubble volume (1b), and driving pressure (1c) in a typical case of inertial cavitation of (f = 20 kHz, R 0 = 3 µm, and p = 1.4). With the choice of the dimensionless timevariable Eq. ( 5), x = 0 represents the time of maximum tension of the liquid. We set

x + = acos 1 p , (15)
and we denote by x m the time of maximum expansion of the bubble, and by x c the time of its maximum compression (see Fig. 1). It is shown in Ref. 1 that, during the expansion phase and most of the collapse phase, the dominant terms in the right-hand-side of Rayleigh equation are the driving term p cos x-1 and also the surface tension term α S R 0 /R for ambient radii just above the Blake threshold. Following Ref. 1, we neglect the dependence of the surface tension term in R , and replace α S R 0 /R by α S /K(p), where K(p) will be determined later. The approximate Rayleigh equation becomes

-x + 1 -p cos x x R R 0      R R 0      3 x + x m x c
RR ′′ + 3 2 R ′2 = R 2 res 3 p cos x -1 + α S K(p) . (16)
We set, for further use

A = 1 + α S K(p) . (17)
Besides, noting that

RR ′′ + R ′2 = 1/2 d 2 (R 2 ) dx 2 , (18)
the right-hand side of Rayleigh equation can be written in two different forms:

RR ′′ + 3 2 R ′2 = 1 2 d 2 (R 2 ) dx 2 + 1 2 Ṙ2 = 3 4 d 2 (R 2 ) dx 2 - 1 2 R R. Numerical simulations show that Ṙ2 ≫ R R on the interval [-x + , x + ], while Ṙ2 ≪ R R holds on the inter- val [x + , x m ] [1]
. Additionally, we found that the latter property still holds in fact during almost all the collapse, except in its ultimate phase, where the gas and acoustic terms become significant again. This could be expected since the main part of the collapse is inertially driven and that Ṙ becomes significant only when the liquid has acquired enough kinetic energy. We therefore obtain the following equations for the bubble radius, over the interval [-x + , x c ]:

d 2 (R 2 ) dx 2 = 4 9 R 2 res (p cos x -A) on [-x + , x + ], (19) d 2 (R 2 ) dx 2 = 2 3 R 2 res (p cos x -A) on [x + , x c ]. ( 20
)
These equations are the same as the ones of Hilgenfeldt et al. [1], except that the validity of the second is extended up to x c . The first equation can be solved with the initial condition R(-x + ) = ζR 0 ,where ζ ≃ 1.6 and Ṙ(-x + ) ≃ R(-x + ) [1]. The second equation is solved by requiring continuity of R(x) and R ′ (x) at x = x + .

Integrating both equations twice, we obtain:

R 2 -(x) = 4 9 R 2 res 1 -p cos x + p(x + x + ) sin x + - A 2 (x + x + ) 2 + ζ 2 R 2 0 [1 + 2 (x + x + )] , (21)
and

R 2 + (x) = 2 3 R 2 res 1 -p cos x + p x 3 + x + sin x + - A 2 x 2 + x 2 + + 2 3 x + x + ζ 2 R 2 0 [1 + 2 (x + x + )] . (22)
The point (x m , R max ) of maximum expansion is obtained by setting d(R 2 + )/dx = 0, so that x m is given in implicit form by

p sin x m -x m + 1 3 (p sin x + -x + ) - α S K(p) x m + 1 3 x + + 3ζ 2 R 0 R res 2 = 0, (23)
and R max reads

R 2 max = R 2 0 f (p, x m ) + R 2 res g(p, x m ) - 2 3 α S K(p) h(p, x m ) , (24)
where

f (p, x m ) = ζ 2 [1 + 2 (x m + x + )] , (25)
g(p, x m ) = 2 3 1 -p cos x m + p x m 3 + x + sin x + - 1 2 x 2 m + x 2 + + 2 3 x + x m , (26)
h(p, x m ) = 1 2 x 2 m + x 2 + + 2 3 x + x m . (27)
In order to obtain x m , the implicit equation ( 23) should be solved. To avoid this, Hilgenfeldt and co-workers [1] developed this equation near π/2 at first order, neglecting on the one hand α S /K(p), and also (R 0 /R res ) 2 , which is appropriate for driving the bubble at low frequencies. They obtain

x m0 = p + 1 3 (p sin x + -x + ) , (28)
which can be further simplified as x m = p, if p is small enough. Plugging the latter into Eqs. ( 24)-( 27), they obtain an expression of R max which depends on R 0 only through the α S term in [START_REF] Löfstedt | [END_REF]. The expression of K(p) is then determined by using the fact, confirmed numerically, that the maximum of the response curve (R max /R 0 )(R 0 ) is obtained for an ambient radius R C 0 very close to the Blake threshold

∂ ∂R 0 R max (p, R 0 ) R 0 = 0 for R 0 = R c 0 = 4 √ 3 9 σ p 0 1 p -1 . ( 29
)
This scheme yields a good approximation for R max , which was the main objective of Hilgenfeldt and co-workers [1], but the approximation (28) of x m yields a rather large error (see dotted line in Fig. 2 and[41]). Since the value of the integral ( 14) is found to be very sensitive to the precise location of x m , we seek a better approximation. We therefore revert to the original equations ( 23)- (27). The main difficulty lies in the presence of the α S term in [START_REF] Mettin | Oscillations, Waves and Interactions[END_REF], which makes rigorously x m a function of both p and R 0 . Thus R max depends not only on R 0 through α S but also through x m in the expressions of f , g and h.

condition (29) therefore becomes more complex, and should be solved simultaneously with [START_REF] Mettin | Oscillations, Waves and Interactions[END_REF]. We initially followed this complex process, but finally found that a better approximation of x m could be obtained by using a simple trick. First, as was done in Ref. 35, we neglect the α S term in [START_REF] Mettin | Oscillations, Waves and Interactions[END_REF] and develop the latter near π/2, but up to second order :

x m1 = π 2 - 1 p + 1 p 1 + 2p × x m0 - π 2 + 3ζ 2 R 0 R res 2 1/2 . ( 30
)
For low frequency driving, R 0 ≪ R res , and x m1 depends only slightly on R 0 . We then plug [START_REF] Prosperetti | proceedings of the NATO Advanced Study Institute on Sonoluminescence and Sonoluminescence[END_REF] in Eqs ( 24)-( 27) and express the condition (29), neglecting ∂x m /∂R 0 , to obtain

K 1 (p) = x 2 m1 + x 2 + + 2 3 x + x m1 g(p, x m1 ) 9 4 √ 3 (p -1). ( 31
)
We now expand again (23) near π/2, but keeping the α S term, in which we set K = K 1 (p), to obtain

x m2 = π 2 - A 1 p + 1 p A 2 1 + 2p × x m0 -A 1 π 2 + (1 -A 1 ) x + 3 + 3ζ 2 R 0 R res 2 1/2 , (32)
where

A 1 = 1 + α S K 1 (p) . ( 33
)
Setting A 1 = 1 in x m2 , that is, neglecting the effect of surface tension, the result of Ref. 35, Eq. ( 30), is recovered.

Figure 2 represents the variations of x m for a bubble of ambient radius R 0 = 1 µm (fig. 2.a) and R 0 = 3 µm (fig. 2.b) in water. The thick solid lines are the exact value obtained numerically, and the thin solid lines represent x m2 . The agreement is seen to be excellent, although a noticeable difference can be seen for R 0 = 1 µm, which originates from the over-simplification done when accounting for surface tension by the simple term α S /K(p) in Eq. ( 16). Also shown is the approximation x m1 (dash-dotted line), which does not take surface tension into account. This clearly introduces an noticeable error on x m , reasonably corrected by Eq. (32). Finally, the approximation x m = p proposed in Ref. [1] is displayed (dotted line).

Finally, the approximation of R max can then easily be obtained by plugging an approximation of x m into Eq. ( 24). This was done in Ref. 35 using x m1 , and an excellent agreement was found. The gain brought by using x m2 instead of x m1 in (24) remains unimportant, and for brevity, we do not present the comparison between the analytical and numerical expressions of R max here.

B. Bubble volume

Approximations of the bubble volume could readily be obtained from the approximations (21), [START_REF] Koch | Proceedings of the IEEE International Ultrasonics Symposium[END_REF] of the bubble radius. However, such expressions do not yield analytical expressions of the integral (14) in closed form, and further approximations are therefore required. First, we consider frequencies low enough to have R res ≫ R 0 , so that the ζ term can be safely neglected in equations ( 21)- [START_REF] Koch | Proceedings of the IEEE International Ultrasonics Symposium[END_REF]. 32); Dash-dotted line: xm 1 from Eq. (30) (ref. [35]); dotted line : xm = p (ref. [1]). Thick dashed line: xc calculated from numerical solutions of equation ( 4); Thin dashed line: xc from Eq. ( 43). The results are calculated for a bubble of ambient radius R0 = 1 µm (a), and R0 = 3 µm (b).

Approximate expression on [-x+, x+]

Numerical simulations demonstrate that R -is almost linear between 0 and x + (Fig. 1.a), which suggests that Eq. ( 21) is almost a perfect square in this interval. We then develop the cos term in (21) near x = 0 up to second order, and write the result as

R 2 -(x) = 4 9 R 2 res p -A 2 x + p sin x + -Ax + p -A 2 + px + sin x + - A 2 x 2 + + 1 -p - (p sin x + -Ax + ) 2 2(p -A) (34)
For R -to be linear in x, the constant term in the bracket must be negligible, so that R -can be simplified as

R -(x) = 2 3 R res p -A 2 x - Ax + -p sin x + p -A . (35)
The expression of the bubble volume on [-x + , x + ] therefore reads

V -(x) = 2 3 R res p -A 2 3 (x -x 1 ) 3 , (36)
where

x 1 = Ax + -p sin x + p -A , (37)
which allows to calculate integral ( 14) in closed form.

Approximate expression on [x+, xc]

Using equations ( 23) and ( 24), it can be easily checked that, setting y = x -x m , the expression [START_REF] Koch | Proceedings of the IEEE International Ultrasonics Symposium[END_REF] of R + can be recast as:

R 2 + = R 2 max + 2 3 R 2 res L(y), (38)
where

L(y) = 2p cos x m sin 2 y 2 -A y 2 2 +p sin x m (sin y -y) . (39)
The bubble volume on [x + , x c ] becomes therefore

V + = R 3 max 1 + 2 3 R res R max 2 L(y) 3/2 , (40)
which unfortunately does not yield an explicit integration of ( 14). Further progress can be done by noting that, from Eq. ( 24), R res and R max are of the same order of magnitude, and that from (39), L(y) = O(y 2 ) near y = 0. Equation ( 40) can therefore be approximated by

V + = R 3 max 1 + R res R max 2 L(y) + O(y 4 ) . (41)
Thus, to the same order of approximation, L(y) can be replaced by any equivalent expression up to order 4 in y, and the choice must be directed by the ability of V + cos(x -x 0 ) to be integrable in closed form. We therefore choose to set y 2 /2 = sin 2 (y/2) + O(y 4 ) and sin y -y = -1/6 sin 3 y + O(y 5 ) in Eq. ( 39) to finally obtain:

V + (x) = R 3 max + R max R 2 res 2(p cos x m -A) sin 2 y 2 - 1 6 p sin x m sin 3 y + O(y 4 ), (42)
which can now yield an explicit expression for integral (14). It can further be noted that neglecting the sin 3 y term in the square bracket, and setting sin 2 y 2 ≃ y 2 /4, V + is found to be zero for

y c = x c -x m = R max R res 2 A -p cos x m 1/2 , ( 43
)
which constitutes a simple approximation of the collapsetime. The comparison between this expression and the exact instant of minimum radius is visible in Fig. 2 (dashed lines). Here again, an excellent agreement is found, but deteriorates toward small bubble radii.

IV. BJERKNES FORCE A. Analytical expression

With the expressions of the bubble volume ( 36) and (42) at hand, the integral ( 14) can be calculated in analytical form, keeping the contribution of the integrand only in the intervals [0, x + ] and [x + , x c ], since V can be neglected in the other regions (see Fig. 1.b). The integral is thus the sum of the two contributions:

I = I -+ I + , (44)
where

I -= x+ 0 V -(x) cos (x -x 0 ) dx , I + = xc x+ V + (x) cos (x -x 0 ) dx.
Using the approximate expressions ( 36) and ( 42) of the bubble volume, integration yields

I -= 8 27 R 3 res p -A 2 3/2 × ∆x(∆x 2 -6) sin(x + -x 0 ) +3(∆x 2 -2) cos(x + -x 0 ) +x 1 (6 -x 2 1 ) sin x 0 + 3(2 -x 2 1 ) cos x 0 , ( 45
) with ∆x = x + -x 1 .
The contribution I + reads

I + = R 3 max [sin(x c -x 0 ) -sin(x + -x 0 )] +R max R 2 res 1 4 (p cos x m -A) × f 2 (y c ) -f 2 (y + ) - 1 192 p sin x m [f 3 (y c ) -f 3 (y + )] , (46)
where f 2 (y) = 4 sin(y -y 0 ) -sin(2y -y 0 ) -2y cos y 0 , f 3 (y) = 2 cos(2y + y 0 ) + cos(4y -y 0 ) +12y sin y 0 -6 cos(2y -y 0 ), and

y 0 = x 0 -x m , y + = x + -x m .
The value of I from ( 44)-( 46) is displayed in figure 3 (thick lines) for R 0 = 1 µm (Fig. 3a), 3 µm (Fig. 3b), and 6 µm (Fig. 3c), in the case of a standing wave (x 0 = 0), for drivings ranging from the Blake threshold to p = 2.5. In order to get a clear picture, I is drawn in logarithmic scale, the solid part of the curves representing a positive sign and the dashed part a negative sign. The thin lines are the results obtained by solving (4) and calculating (14) numerically, for f = 20 kHz. It is seen that an excellent agreement is obtained, except for R 0 = 1 µm (Fig 3a). Particularly, the point of inversion of the Bjerknes force is shifted toward large drivings. This feature originates from the errors induced on the values of x m , x c (see Fig. 2) and R max for small ambient radii, by replacing the surface tension in the RP equation by α S /K(p) in (16). It should be noticed that even the small errors visible on the curves of Fig. 2.a yields large differences on the estimation of I. This could be expected since the phase between V and cos(x -x 0 ) crucially influences the value of integral I.

B. Bjerknes force inversion threshold in standing waves

We consider the case of a standing wave x 0 = 0, and look for an approximate locus in the parameter space where the Bjerknes force changes sign. Summing equa-tions (45) and (46), it is seen that integral ( 14) is zero for

a 3 X 3 + a 1 X + a 0 = 0 (47)
where

X = R max R res ,
and the coefficients a i depend on

• x + , which is just acos(1/p),

• x 1 , which from (37) depends on p and α S ,

• x c , which from (43) depends on p, x m , X and α S ,

• x m , which from (32), only depends on p and α S , for R 0 /R res ≪ 1.

Furthermore, looking at Eq. ( 24), for R 0 ≪ R res , X = R max /R res can be written as

X = g(p, x m ) - 2 3 α S K(p) h(p, x m ) 1/2 , (48)
and from Eqs. ( 25)-( 27) and ( 30)-( 33), still under the assumption R 0 ≪ R res , the terms g, h in the above equation depend on R 0 only through α S . We conclude that, provided that R 0 ≪ R res , equation (47) becomes frequency independent, and can in fact be written in implicit form as

I(α S , p) = 0. ( 49
)
This equation can easily be solved for α S (p), in order to find the approximate, frequency-independent threshold for inversion of the Bjerknes force. The solution is presented in the inset of Fig. (4). Below the curve, I > 0, so that the Bjerknes force attracts the bubble toward pressure antinodes, while it becomes repulsive above.

From α S = 2σ/(p 0 R 0 ), the inversion threshold can also be plotted in the (R 0 , p) plane in the case of water at ambient pressure (σ = 0.072 N.m -1 , p 0 = 101300 Pa). The result is displayed in Fig. 4 (thick solid line) and compared to the exact inversion thresholds calculated from numerical simulation for three driving frequencies 20 kHz (dash-dotted line), 40 kHz (dashed line), and 80 kHz (thin solid line). The labels on the two latter curves represent the value of R 0 /R res . It is seen that the above procedure yields a good estimation of the inversion threshold, up to R 0 /R res = 0.1, above which it starts to diverge from the exact value. The reasons for this disagreement comes from the neglected R 0 /R res term in all expressions, and also from the fact that for increasing frequency, the bubble rebounds become more important, so that the bubble dynamics for x > x c also contributes to expression (14). Besides, a cascade of period-doubling bifurcations and chaos [START_REF] Lauterborn | proceedings of the NATO Advanced Study Institute on Sonoluminescence and Sonoluminescence[END_REF]36,37] appear in some cases (and are responsible for the noisy oscillations on the 80 kHz curve), so that the correct averaging of the Bjerknes force in such cases should be carried out over more than a single acoustic period. We did not pursue further this issue, since analytical predictions for these bifurcations are out of the scope of the present paper. Marginally, it can be seen that the inversion threshold in the (α S , p) plane is almost linear, so that the following linear fit (represented by a dashed line in the inset of Fig. 4) can be proposed for practical applications:

p = 0.269 α S + 1.62. ( 50
)
These results suggest that the inversion threshold is independent of frequency, and of the properties of the gas and liquid other than surface tension, as long as R 0 /R res ≪ 1. This astonishing result originates from the fact that the Bjerknes force mainly depends on the expansion phase of the bubble, which, within the approximations leading to Eq. ( 16), is merely governed by the driving pressure amplitude and surface tension. The reasonably good agreement found in Figs. 234partially supports this analysis.

In order to further investigate this issue, we first recalculated the three inversion thresholds of Fig. 4 (f = 20, 40, 80 kHz), replacing the thermal model of Ref. 32 by an isothermal behavior for the bubble interior. Figure 5 displays the results obtained (thick solid lines) and recalls the thresholds calculated in Fig. 4 (thin solid lines). It can be seen that the thresholds slightly diverge for increasing R 0 but remain almost indistinguishable for 0 /R res < 0.15 We also repeated the calculations with the thermal model of Ref. 32, but for argon bubbles (not shown), and found a negligible deviation from the air curves. We therefore conclude that the detailed bubble interior has a very weak influence on the expansion phase, at least for low enough values of R 0 /R res , so that Eq. ( 50) indeed constitutes a gas-independent law, within its range of validity. Another issue is the sensitivity of the results to the liquid viscosity. The latter has been neglected in the analytical approach, when approximating the RP equation (8) by Eq. ( 16). The good agreement found in Fig. 4 between analytical and numerical results, calculated for water at ambient temperature (µ = 10 -3 Pa.s), suggests that for such low values, viscosity indeed plays a minor role during the bubble expansion. One should however check whether it is still the case for larger viscosities. We therefore repeated the calculation of the inversion threshold for viscosities 10 and 20 times larger than the one of water (Fig. 6, thick dashed line and thick dashdotted line). It is clearly seen that the threshold increases noticeably with viscosity. Conversely, we also checked that the result was unaffected by decreasing the viscosity below the water's one, by computing the threshold for µ = 0.1 µ water (thick solid line). This indicates that viscous friction plays a non-negligible role in the bubble expansion for viscosities above some critical value. As already mentioned in Ref. 1, increasing viscosity decreases R max , and we also checked that it decreases x m too, so that, strictly speaking, the Bjerknes force and its inversion threshold are viscosity-dependent. Following our results, this influence is negligible for viscosities near or lower than the water's one, but for slightly larger values, the viscous term should be kept in the RP equation.

V. DISCUSSION

Important conclusions can be drawn from these results. Figure 4 shows that the inversion thresholds for all frequencies (thin lines) asymptotically merge with the Blake threshold (thick dashed line) for small bubble radii, and in reasonable agreement with the analytical approximation (thick solid line). Thus, as the driving pressure reaches, say 1.8 bar, the range of ambient radii of inertial bubbles attracted toward the antinode is suddenly reduced, with an upper limit lower than 2 µm. This explains why a well-defined bubble-free region can be observed around the pressure antinode for high amplitude standing waves [21]. The range of attracted bubbles is however not void, which suggests that the zone around the antinode could still be filled with inertial bubbles, of ambient radii very close to the Blake threshold, but too small to be visible. As noticed in Ref. [18], in a high amplitude standing wave, the Bjerknes force acts as a sorter of inertial bubbles, leaving the smallest ones approaching or even reaching the pressure antinodes. The advantage of the present analysis is that it yields, through Eq. ( 49), or its simpler form (50), an explicit classification of the bubble sizes as a function of the local acoustic pressure, parametrized by the ratio σ/p 0 .

As the increasingly small bubbles approach the pressure antinode, they may coalesce or quickly grow by rectified diffusion [35]. Increasing their size, they may again enter the repulsion zone in the (R 0 , p) plane and move back again. This picture is still complicated by the potential appearance of surface instabilities. Thus, the apparently void region observed around large pressure antinodes may be in fact the locus of the complex evolution of very small bubbles, of sizes close to the Blake threshold.

Finally, it is seen from the inset of Fig. 4 that decreasing α S lowers the driving at which the pressure antinode becomes repulsive. The dimensionless parameter α S can be varied experimentally by modifying the surface tension σ (for example adding ionic salts or surfactants), or by changing the static pressure p 0 . The present results suggest that, for identical bubble ambient radii, the Bjerknes force would become repulsive for lower drivings, when either decreasing σ or increasing p 0 . This should have an observable effect on the size of the bubble-free region around the pressure antinode. However, it should be noted that surface tension also plays a crucial role for bubbles surface instabilities [10,38,39], and also for rectified diffusion [35], through the same dimensionless parameter α S . Thus, changing α S may also directly influence these two processes, with probable consequences on the bubble cloud behavior. The present result just demonstrates that surface tension can influence the shape of the bubble cloud through its direct effect on the bubble dynamics, and on the primary Bjerknes force.

Figure 6 also indicates that the size of the bubblefree region around the pressure antinode would decrease noticeably when increasing viscosity slightly above the one of water. As mentioned in Ref. 1, this may be easily achieved experimentally by adding glycerin in water.

Here again, such a macroscopic effect is mediated by the sensitivity of the bubble dynamics to the physical properties. To account analytically for this dependence on viscosity, the viscous term should be kept in the Rayleigh equation, which renders the approximation scheme more involved. A generalization of our analytical results to this case may be addressed in a future study.

Finally, it is highly probable that the same effect of surface tension could be observed on the secondary Bjerknes force, as suggested by numerical simulations [40]. The extension of the present analytical method to the latter effect is difficult, first because the expression of the secondary Bjerknes force also involves the bubbles velocities, which are much more sensitive to approximations than the bubble radius itself, and secondly because the dynamics equation of the two bubbles must be coupled by a radiation term.

FIG. 1 :

 1 FIG. 1: (a) Dimensionless bubble radius R/R0; (b) Dimensionless bubble volume (R/R0) 3 ; (c) Dimensionless driving pressure 1 -p cos x. The case considered is a 3 µm air bubble in water and p = 1.4. The times -x+ and x+ are the two instants of zero-crossing of the driving pressure, xm is the time of maximum expansion of the bubble, and xc the time of maximum compression. The dashed curve in (a) represents the approximate dynamics given by Eqs. (21), (22). The dashed line in (b) is the final approximation of the bubble volume (36)-(42).

FIG. 2 :

 2 FIG.2: Thick solid line: xm calculated from numerical solutions of equation (4); Thin solid line: xm 2 from Eq. (32); Dash-dotted line: xm 1 from Eq. (30) (ref.[35]); dotted line : xm = p (ref.[1]). Thick dashed line: xc calculated from numerical solutions of equation (4); Thin dashed line: xc from Eq. (43). The results are calculated for a bubble of ambient radius R0 = 1 µm (a), and R0 = 3 µm (b).

FIG. 3 :

 3 FIG. 3: Thick lines: value of I predicted by approximation (44)-(46), for x0 = 0 (standing wave), for an air bubble in water, of radius R0 = 1 µm (a), 3 µm (b) and 6 µm (c). The solid parts of the curves correspond to I > 0 and the dashed parts to I < 0.

FIG. 4 :

 4 FIG. 4: Threshold of Bjerknes force inversion in the (R0, p) plane, for a bubble in water in ambient conditions (σ = 0.072 N.m -1 , p0 = 101300 Pa, µ = 10 -3 Pa.s). The region I > 0 corresponds to attraction by the pressure antinode, and I < 0 to repulsion. The thin lines are calculated from numerical simulations of the RP equation. Thin solid line: f = 80 kHz; dashed line: f = 40 kHz; dash-dotted line:f = 20 kHz. The labels on the curves indicate the ratio R0/Rres (triangles: f = 80 kHz; filled circles: f = 40 kHz). Thick solid line: universal threshold calculated from approximate dynamics by solving (49). Thick dashed line: Blake threshold. The inset represent the solution of (49) in the (αS, p) plane.

FIG. 5 :

 5 FIG.5: Same as Fig.4. The thin solid lines are the numerical curves of Fig.4(f = 20, 40, 80 kHz). The thick solid lines are calculated in the same conditions, except that the gas behavior is considered isothermal. The thin dashed line is the analytical threshold calculated from Eq. (49).

FIG. 6 :

 6 FIG.6: Same as Fig.4, only for f = 20 kHz, and for different liquid viscosities. Thin solid line: water (same as dash-dotted line of Fig.4); Thick solid line: µ = 0.1 µwater; Thick dashed line: µ = 10 µwater; Thick dash-dotted line µ = 20 µwater. The thin dashed line is the analytical threshold calculated from Eq. (49).


O Louisnard

J Gonzalez-Garcia

I Tudela

J Klima

V Saez

Y Vargas-Hernandez

FEM simulation of a sono-reactor accounting for vibrations of the boundaries

Keywords: finite elements, ultrasound, standing waves, fluid-structure interaction, acoustic cavitation, PACS: 43.20.Ks PACS: 43.20.Tb PACS: 43.40.At PACS: 43.35.Ei

The chemical effects of acoustic cavitation are obtained in sono-reactors built-up from a vessel and an ultrasonic source. In this paper, simulations of an existing sono-reactor are carried out, using a linear acoustics model, accounting for the vibrations of the solid walls. The available frequency range of the generator (19 kHz-21 kHz) is systematically scanned. Global quantities are plotted as a function of frequency in order to obtain response curves, exhibiting several resonance peaks. In absence of the precise knowledge of the bubbles size distribution and spatial location, the attenuation coefficient of the wave is taken as a variable, but spatially uniform parameter, and its influence is studied. The concepts of acoustic energy, intensity, active power, and source impedance are recalled, along with the general balance equation for acoustic energy, which is used as a convergence check of the simulations. It is shown that the interface between the liquid and the solid walls cannot be correctly represented by the simple approximations of either infinitely soft, or infinitely hard boundaries. Moreover, the liquid-solid coupling allows the cooling jacket to receive a noticeable part of the input power, although it is not in direct contact with the sonotrode. It may therefore undergo cavitation and this feature opens the perspective to design sono-reactors which avoid direct contact between the working liquid and the sonotrode. Besides, the possibility to shift the main pressure antinode far from the sonotrode area by exciting a resonance of the system is examined.

Introduction

When a liquid is irradiated by a high-power ultrasonic wave, numerous radially oscillating micronsized bubbles appear. The phenomenon is known as Email address: louisnar@enstimac.fr (O. Louisnard). 1 Corresponding author acoustic cavitation [START_REF] Leighton | The acoustic bubble[END_REF][START_REF]Sonochemistry and Sonoluminescence[END_REF]. The strong collapse, following the explosive expansion of these bubbles, induces extreme conditions inside or near the bubbles, which are responsible for a specific chemistry, known as sonochemistry [START_REF] Suslick | Ultrasound. Its chemical, Physical and Biological effects[END_REF].

Various experimental devices can be used to produce this phenomenon. The most common one is a horn transducer, diving in the liquid and creating large acoustic pressure in its vicinity. Other systems involve a transducer with larger emitting area, from the bottom (sometimes referred as "cup-horn reactors"), or various transducers sticked to plane walls, the latter system being commonly known as "ultrasonic bath".

There is a common thought that thin horn transducers emit attenuated spherical waves, decaying rapidly as the distance from the tip increases, while large-area transducers rather yield standing waves, by reflexion of the emitted waves on the vessel wall [START_REF] Yasui | [END_REF]. In the latter configuration bubbles gather near the points where the acoustic amplitude reaches a critical value at which the Bjerknes force undergoes a sign change [5], and form Jellyfish-like structures [START_REF] Mettin | Bubble and Particle Dynamics in Acoustic Fields: Modern Trends and Applications, Research Signpost[END_REF][START_REF] Lauterborn | [END_REF], while, in the former, a dense bubble cloud forms just below the sonotrode tip [START_REF] Yasui | [END_REF]8].

While the above observations are indeed common in numerous experiments, it does not rule out that different behaviors could be observed in specific geometric cases. Following an interesting optimization procedure, Klima and co-workers [9] have demonstrated that standing waves, materialized by a cavitation cloud in the center of the vessel, could also be obtained in an experiment with a thin horn transducer, provided that the geometry of the vessel and the liquid height were carefully tuned. This result conforms to theoretical acoustics, which indeed states that any liquid bounded by solid walls or free liquid-gas interfaces presents normal vibration modes at definite frequencies, and resonates when excited at these frequencies [START_REF] Morse | Theoretical acoustics[END_REF]. Some available commercial systems exploit this property, and adjust automatically the working frequency (in a reasonable range), in order to possibly excite a normal mode of the system. One of the motivations of this paper is to assess by simulation the existence of the resonance modes of a given sono-reactor.

The question arises if such resonant behaviors can be reasonably predicted by linear acoustics. At first sight, it might sound doubtful, since high-amplitude waves are subject to nonlinear phenomena, even in homogeneous liquids. Moreover, the presence of cavitation bubbles modifies the acoustic properties of the effective medium, decreasing the effective sound velocity, and introducing dispersion and nonlinear phenomena [START_REF] Van Wijngaarden | [END_REF]12,13]. However, in most cases, the bubbles are concentrated in relatively small regions of the liquid, and one could expect that linear acoustics may at least give a qualitative idea of the acoustic field and the approximate location of the various resonance frequencies, in the first step of the design of a sono-reactor.

Even in the restricted frame of linear acoustics, there remains the problem of the suitable boundary conditions. Precedent studies generally treat the vessel boundaries either as infinitely rigid walls [14], or infinitely soft [9]. Solid boundaries may in fact be thin enough to vibrate and deform, and these effects must be considered, if one seeks precisely the resonant modes of the system. Simulations of a sono-reactor accounting for vibrations of the vessel wall have been recently published [15], and show the effect of the vessel wall thickness in a rectangular cell for two frequencies, 100 kHz and 140 kHz.

Here we present the response of an experimental sono-reactor, including a cooling jacket, excited by a cup-horn type transducer. The simulations presented hereafter account for the vibration of the boundaries, and possible transmission of sound into the cooling jacket. The frequency is systematically varied in the range available for the generator used, the resonances of the whole structure are sought and characterized by calculating several global quantities as functions of frequency.

Besides, we recall the concepts of acoustic energy density, acoustic intensity, active power, and source impedance, which are sometimes confusingly used. It is first recalled, on rigorous arguments, that if no wave escapes from the sono-reactor, the active power is nothing else than the power dissipated in the liquid. The latter can be estimated by recording the initial temperature increase, which forms the basis of the calorimetric method [16].

Equations of the problem

Linear acoustics in the liquid

The linear propagation of an acoustic wave is described by the linearization of the Euler equations [START_REF] Morse | Theoretical acoustics[END_REF]:

1 ρ l c 2 l ∂p ∂t + ∇.v = 0, (1)
ρ l ∂v ∂t = -∇p, (2)
where v(x, t) is the liquid velocity field associated to the acoustic wave, p(x, t) the local acoustic pressure, and ρ l , c l are the density and the sound speed of the liquid, respectively.

We assume mono-harmonic waves at angular frequency ω. In what follows, we will use the complex notation p(x, y, z, t) = 1 2 P (x, y, z)e iωt + P * (x, y, z)e -iωt ,

v(x, y, z, t) = 1 2 V(x, y, z)e iωt + V * (x, y, z)e -iωt . (3)
Equations ( 1), (2) transform to

iωP + ρ l c 2 l ∇.V = 0, ( 5
) iωρ l V = -∇P, (6)
and eliminating V between the two equations yields the well-known Helmholtz equation

∇ 2 P + k 2 l P = 0, (7)
where k l = ω/c l is the wavenumber.

Boundary conditions for the liquid

The free surface in the reactor, if any, is generally represented by an infinitely soft boundary, that is p = 0, or P = 0 in the case of mono-harmonic waves.

On the radiating surface of the transducer, one generally assumes that the normal displacement is known. This is an over-simplification, and a correct representation would involve the modelization of all parts of the transducer, including the piezoceramics, and an appropriate coupling with the fluid. Here, in a first approach, we will however use this simplification. Noting U 0 the complex amplitude of the transducer displacement, and using (6), the boundary condition for the fluid reads

∇P.n = ρ l ω 2 U 0 , (8)
where n is the normal pointing outward the liquid.

In past studies, the walls of the vessel enclosing the liquid have been either also considered as infinitely soft (P = 0) [9], or infinitely rigid [14,17,18]. The latter condition imposes a zero normal velocity, which, for mono-harmonic waves, using [START_REF] Mettin | Bubble and Particle Dynamics in Acoustic Fields: Modern Trends and Applications, Research Signpost[END_REF], translates to ∇P.n = 0. To relax any of these two simple approximations, the vibration of the solid parts must be properly accounted for.

Vibrations of the solid

Neglecting the volumic forces, the vibrations of an elastic solid are given by

ρ s ∂ 2 u s ∂t 2 = ∇.σ, (9)
where ρ s is the solid density. u s (x, y, z, t) is the displacement field and σ(x, y, z, t) the elastic stress tensor given by

σ = Eν (1 -2ν)(1 + ν) (Tr ) I + E 1 + ν , ( 10
)
where E is the Young modulus, ν the Poisson ratio, I the identity tensor, Tr the trace operator, and (x, y, z, t) the strain tensor:

= 1 2 ∇u s + T ∇u s . (11)
For mono-harmonic vibrations, the displacement field is set as

u s (x, y, z, t) = 1 2 U s (x, y, z)e iωt + U s * (x, y, z)e -iωt , (12)
so that Eq. ( 9) becomes

-ρ s ω 2 U s = ∇.Σ, (13)
where

Σ = Eν (1 -2ν)(1 + ν) (TrE) I + E 1 + ν E, (14)
E = 1 2 ∇U s + T ∇U s . (15)

Liquid-solid interface conditions

The coupling between the liquid and solid vibrations are obtained by writing two interface conditions on such surfaces. The first one is a cinematic condition, and just states that the displacement is continuous on each side of the interface, so that the condition

v.n = ∂u s ∂t .n (16)
must hold, n being the normal unit vector pointing outward the liquid. For mono-harmonic waves, derivating this relation and using [START_REF] Mettin | Bubble and Particle Dynamics in Acoustic Fields: Modern Trends and Applications, Research Signpost[END_REF] yields

∇P.n = ρ l ω 2 U s .n, (17)
which is used as a boundary condition for linear acoustics in the liquid.

The second condition is a dynamic one, and states that the normal force per unit area exerted on the solid boundary in contact with the liquid is just the pressure of the liquid, which reads

Σ.n = -P n. ( 18
)
This continuity equation is used as a boundary conditions for equation (9) governing the vibration of the solid.

Other boundary conditions for the solid

On the boundaries of the solid in contact with air (typically for the outer boundaries of the vessel), the three components of the stress are ascribed to zero.

Energy consideration

Conservation of acoustic energy

If attenuation due to energy absorption is ignored, the equation of acoustic energy conservation can be readily deduced from Eqs. (1), (2), and reads [START_REF] Morse | Theoretical acoustics[END_REF]

d dt V e a dV = S -I.n dS, ( 19
)
where V is an arbitrary volume of fluid, S its boundary. The quantity

e a = 1 2 ρ l v 2 + 1 2 p 2 ρ l c 2 l ( 20
)
is the acoustic energy density (in W/m 3 ), which is the sum of the kinetic energy and the potential compressional energy of the liquid, and

I = pv (21)
is the acoustic intensity (in W/m 2 ). Both e a and I are local quantities, and it should be recalled that I is a vector field, which streamlines represent the path followed by acoustic energy. Equation (19) just expresses that the variations of acoustic energy in a volume of fluid, results from the difference between the energy fluxes entering and leaving this volume.

From a mechanical point of view, this equation is just the consequence of the kinetic energy theorem: the surface integral is the power of external pressure forces on S, the integral of 1 2 ρ l v 2 is the total kinetic energy in V , and the volume integral of 1 2 p 2 /(ρ l c 2 l ) is the opposite of the power of internal pressure forces (see also appendix A).

Active power, attenuation and energy dissipation

We consider the case of a sono-reactor containing a volume V of liquid, excited by the vibrating surface S sonotrode of a sonotrode, and closed by boundaries S boundaries , which may be a liquid-air surface or the internal boundary of the reactor's walls. It is shown on formal grounds in appendix A that, for monoharmonic waves in steady state, the period-averaged energy conservation for the liquid reads:

P active = P diss + P bound , (22)
where

P active = S sonotrode - 1 2 (P V * ).n dS, (23)
P bound = S boundaries 1 2 (P V * ).n dS, (24)
P diss = V α |P | 2 ρ l c l dV. ( 25
)
-P active is the power entering the liquid through the sonotrode area, and is termed as active power, -P bound is the power lost through the other (passive) boundaries of the liquid. This quantity vanishes for either infinitely soft boundaries, or nondissipative unconstrained elastic solid walls (see appendix A). This will be the case in this study so that P bound = 0. -P diss is the power dissipated in the medium, and, as shown in appendix A, is related to the attenuation coefficient α described below. Eq. ( 22) is a generalization of Eq. ( 19), applied to mono-harmonic waves, when dissipative phenomena come into play within the liquid. In absence of dissipation, P diss is zero, and ( 22) predicts a zero active power. In the opposite case, integral (25) is positive, and equation (22) states that the active power entering through the sonotrode is dissipated in the medium, if it cannot flow out through the other boundaries. This relation, along with expressions ( 23) and [START_REF] Mettin | Sixth International Symposium on Cavitation -CAV2006 (paper 75[END_REF], will be useful to check the validity of the simulation results.

Dissipation in acoustic waves occurs for various reasons: viscosity of the medium, diffusion of compressional heat energy at finite rate, and relaxation delay in the case of gases. For linear mono-harmonic wave, this effects can be represented by a complex wavenumber (or sound velocity)

k = k r -iα, (26)
where α > 0 is the attenuation coefficient, and k r = ω/c l . In the case of cavitation, dissipation occurs at the level of each bubble, and originates from viscous friction, heat diffusion in the bubble, and acoustic radiation [19,20,21,22]. The corresponding sound velocity and attenuation factor of linear acoustic waves in a bubbly liquid can be expressed easily, once the bubbles size-distribution function is known [12]. The latter theory has been used by several authors, but setting an arbitrary Gaussian-shaped distribution of bubble sizes [17,18,23]. Moreover, it is well known that, owing to Bjerknes forces, the cavitation bubbles are not homogeneously distributed in the whole liquid, but rather arrange in complex localized structures [START_REF] Mettin | Bubble and Particle Dynamics in Acoustic Fields: Modern Trends and Applications, Research Signpost[END_REF]. Owing to the lack of theory predicting the bubble spatial and size distribution that would allow the correct calculation of the attenuation factor, we chose to treat the latter quantity as a variable, but spatially uniform, parameter in the simulations presented hereafter.

Source impedance, plane wave and baffled piston hypothesis

The active power in the energy conservation equation ( 22) can be interestingly connected to the concept of source impedance. The general term acoustic impedance depicts the ratio P/V . Assuming that the sonotrode area vibrates uniformly perpendicular to its plane, the source impedance can be defined as:

Z = 1 S sonotrode V 0 S sonotrode P dS, (27)
where V 0 = iωU 0 is the complex uniform source velocity. From this relation and ( 23), it is readily seen that the active power can be written as

P active = 1 2 (Z)|V 0 | 2 S sonotrode . (28)
It is therefore seen from ( 22) that the real part of the source impedance originates either from the power dissipated in the medium (P diss = 0), or from the waves escaping the medium, for example in open space (P bound = 0). When a vibrating area emits plane traveling waves (PTW) in open space, it can be easily shown that the source impedance is just the specific impedance of the medium:

Z = ρ l c l (29)
This hypothesis becomes valid when the source radius a is much larger than the wavelength (ka 1), but we emphasize that its validity is restricted for waves propagating in open space, without any reflecting boundary. It is therefore expected, and will be checked numerically, that this approximation definitely yields unrealistic results for a sono-reactor, which is by essence a closed geometry.

In our case, the dimensionless number ka is approximately 1.3 so that the plane wave hypothesis is anyway not expected to hold. A more general expression of the source impedance for waves emitted in open space can be obtained under the hypothesis of a baffled piston, which depicts a circular vibrating source, embedded in a rigid plane, emitting in half-space. The source impedance reads in this case:

Z = ρ l c l 1 - J 1 (2ka) 2ka + iM 1 (2ka) , (30)
where J 1 and M 1 are, respectively, the Bessel function of the first kind, and the Struve function [START_REF] Morse | Theoretical acoustics[END_REF][START_REF] Morse | Methods of theoretical physics[END_REF]. Plane traveling waves and spherical waves are recovered for ka → ∞ and ka → 0, respectively. Here again, it must be recalled that this expression of the source impedance is only valid for waves emitted in an infinite domain.

The problem solved

Experimental configuration studied

The experimental setup simulated is represented in Figure 1. Since it presents axial symmetry, only one half of a cut plane is represented. The sonotrode (diameter 30 mm) emits from below in a liquid in contact with air on its upper surface, and limited laterally by thin glass walls. The sono-reactor also presents a cooling jacket to maintain the temperature constant all along an experiment. A Teflon ring prevents the fluid to leak downward and maintains the whole setup. The liquid in the reactor and in the cooling jacket is water with properties ρ l = 1000 kg.m -3 and c l = 1500 m.s -1 . The elastic properties are E = 73 GPa, ν = 0.17 for glass and E = 0.5 GPa, ν = 0.46 for Teflon. The liquid height is fixed to h = 72 mm, counted from the vibrating surface of the sonotrode. The latter is assumed to vibrate vertically with an amplitude of 10 -6 m.

The sonotrode is constituted internally from a sandwich transducer excited by a variable frequency generator (20 kHz -100 W, Undatim) [14]. The system is able to scan the emission frequency between 19 kHz and 21 kHz, and records the frequency for which the electrical impedance on the transducer terminals is minimal. This allows experimentally to work at a frequency at which the system presents a resonance, if any.

Simulation strategy

The problem defined in section 2 is solved with the COMSOL software, varying the frequency by small steps in the range 19 kHz -21 kHz. Scalar quantities representative of the global vibrations of the system are computed and plotted as functions of frequency, in order to get a clear picture of the system's resonances.

The computed quantities can be chosen among the following: -the average acoustic energy stored in the liquid:

E a,l = V 1 4 ρ l |V| 2 + 1 4 |P | 2 ρ l c 2 l dV, (31)
-the average acoustic energy stored in the solid walls, either Teflon and glass. The latter quantity can be written as

E a,s = V walls 1 4 ρ s ω 2 |U s | 2 + 1 2 (Σ : E * ) dV, (32)
where the first term represents the kinetic energy density, and the second is the deformation energy density of the solid.

-the average acoustic pressure in the liquid:

P av = 1 V V |P | dV, (33)
-the power dissipated in the liquid

P diss = V α |P | 2 ρ l c l dV, (34)
where V may be the volume of the working liquid, the volume of the cooling liquid in the jacket, or the sum of both, -the source impedance, given by Eq. ( 27) -the active power entering the liquid through the sonotrode Eq. ( 23) As seen above, from Eq. ( 22) with P bound = 0, the two quantities P diss and P active should be equal, which is checked for each simulation as a convergence test. Besides, a mesh convergence study has been performed.

Results and discussion.

Response curves

Figure 2 displays the mean acoustic pressure in the reactor as a function of frequency, with zero attenuation, in the following cases: cooling jacket filled with water (thick solid line), cooling jacket empty (thin solid line), internal walls considered as hard boundaries (dashed line), or as soft boundaries (dash-dotted line). In the two latter cases, the vibration of the glass and Teflon walls are not accounted for, and it is clearly seen that the response curves differ significantly from the ones obtained by properly taking into account the walls vibration.

The peaks observed in the response curve correspond to global resonances of the whole mechanical system formed by the liquid coupled with the solid walls. In absence of dissipation, the amplitude of these peaks should be infinite, which would traduce an infinite accumulation of acoustic energy in the system. Thus, the finite maxima of the peaks appearing on Fig. 2 are only due to the discretization of the frequency axis. It is interesting to note that the response curves, with the cooling jacket respectively filled and empty, present approximately the same set of peaks, shifted in frequency, except for peak II which is specific to the filled cooling jacket, which already demonstrates the strong influence of the operating conditions on the response of the system.

We now focus on the case where the cooling jacket is filled with water (Fig. 2, thick solid line). 20610 Hz). White zones correspond to large positive acoustic pressures, while black ones correspond to large negative acoustic pressures. Both black or white regions therefore indicate the loci of pressure antinodes, while gray ones represent pressure nodes. The deformed boundary (magnified) is also displayed, in order to illustrate the coupling between the vibrations of the fluid and the solid. The lat-ter effect is remarkable in the breathing behavior of the cooling jacket: for ωt = 0, the negative pressure sucks the glass walls, while the positive pressure repels them for ωt = π. Similar behaviors have been evidenced in Ref. [15]. The uniform gray color observed for ωt = π/2 corresponds to the phase at which the acoustic pressure crosses zero uniformly in the reactor. At this time, the velocity field reaches its peak value both in the liquid and the solid. This behavior is well-known for standing waves, and Eqs. ( 22), (23) show that in absence of dissipation, pressure and velocity are indeed in phase quadrature, so that, from Eq. (A.5), the acoustic intensity is zero everywhere. It is thus retrieved that in absence of dissipation the acoustic power is zero.

And interesting feature observable in Fig. 3 is the presence of a pressure antinode inside the cooling jacket, which suggests that cavitation can appear in a cavity insulated from the sonotrode by a glass wall. This is indeed often observed experimentally, and we emphasize that this effect demonstrates clearly the role of the wall vibrations, and the need of their proper representation in reliable simulations.

Another important point is that, at this resonant frequency, the pressure antinode appears far from the sonotrode. The emitting area of sonotrodes are known to suffer erosion, because cavitation near solid surfaces produces liquid jets and shock waves, responsible for the destructive effects [START_REF] Lauterborn | [END_REF]. The present simulation suggests that, exciting the sono-reactor near a resonant frequency could help to avoid this problem, as was evidenced in ref. [9].

To pursue this issue, it is instructive to examine the pressure field obtained far from a resonance, for example at f = 19500 Hz (see the thick solid line in Fig. 2). Figure 4 shows that in this case, the pressure antinode is located near the sonotrode, and that the pressure amplitude rapidly decreases with the distance from the latter. This behavior is commonly observed experimentally, and the decay of acoustic pressure far from the sonotrode is generally attributed to strong attenuation of the wave by cavitation bubbles. Since no attenuation was introduced in the simulation of Fig. 4, the present result demonstrates that this is not necessary the case. The pressure field observed rather indicates that the sonotrode mainly emits a diverging wave, which has strongly decayed in amplitude when it reaches the boundaries (see also section 5.4). Nevertheless, as seen in Fig. 4, it is still strong enough to excite the breathing mode of the cooling jacket.

Effect of attenuation

The orders of magnitude of pressure indicated in Fig. 2, 3, 4 may sound unrealistic, especially near resonance. This is expected since no attenuation was introduced in the simulations. We therefore repeated the computations for different values of the attenuation coefficient: α = 0.05, 0.5, 1 and 5 m -1 . Figure 5 displays the response curves obtained. Far from resonances, the effect is weak, except for the largest value α = 5 m -1 (thick solid line). Near resonances, the mean acoustic pressure decreases strongly with attenuation, and some resonance peaks even disappear when increasing attenuation, especially peaks III and V. We refer the reader interested in a fine analysis of the resonance peaks to appendix B.

Attenuation not only influences the amplitude of the acoustic field, but also the structure of the pressure field. This is evidenced in Figure 6, which represents the spatial distribution of the peak pressure amplitude |P | in the reactor, near resonance peak IV, and for increasing α (from left to right, with different color-scales for each figure). A first structural change occurs between the first two fields, as attested by the contour lines and the shape of the glass walls. However, in the first three cases, the pressure antinode remains in the middle of the reactor. The shapes of the field and boundaries deformation are identical for α = 0.5 and 1. For the highest attenuation coefficient α = 5, the pressure antinode comes back near the transducer, as would be observed far from resonances (like in fig. 4). This last feature therefore chastens the above remark about the utilization of resonances to keep the pressure antinode away from the tip of the sonotrode. It is difficult to make further progress concerning this issue, in absence of a more realistic model for attenuation, which is mainly due to the bubbles. If the resonance field of the left-hand graph of Fig. 6 could be obtained, the bubbles would gather at the center near the pressure antinode, and introduce some attenuation there (and also shift the local sound speed). The key issue is to know whether the resonant field structure predicted by linear acoustics would be just shifted in frequency and attenuated by the presence of the bubbles, or completely destroyed. The experimental results of ref. [9] are encouraging in this way, since a resonant field is still visible in the presence of (also visible) cavitation bubbles. Besides, the simulations and experiments of Mettin and co-workers [START_REF] Mettin | Sixth International Symposium on Cavitation -CAV2006 (paper 75[END_REF] suggest that a steady state is not necessary reached and that a coupled periodic destabilization of the acoustic and bubble fields can occur.

Active power and real source impedance

As seen above, the active power is also the dissipated power and it is therefore expected to vary with attenuation. From Eq. (28), the active power is proportional to the real part of the source impedance, so that we chose to draw (Z)/ρ l c l (Fig. 7). Near resonances, the active power drastically decreases as the attenuation coefficient is increased. Conversely, far from resonances, the opposite occurs: the higher the attenuation, the higher the active power. The latter behavior may sound non-intuitive, but it can be seen from equation ( 25) that the dissipated power indeed increases both with the attenuation coefficient α and the square of the acoustic pressure amplitude |P |.

Far from resonance the amplitude of the sound field |P | varies slightly with the attenuation coefficient, so that the variation of P diss is mainly governed by the variation of α. Conversely, near resonance, the pressure amplitude at the antinode increases drastically, and a slight increase in the attenuation coefficient yields a huge increase of energy dissipated near the antinode, thus increasing the active power.

In mode physical words, near resonance, the liquid stores more acoustic energy, and can therefore dissipate more. The real impedance of a plane traveling wave is materialized in Fig. 7 by the horizontal thick dotted line (Z)/ρ l c l = 1, and the thick solid line represents the case of a baffled piston Eq. (30). It can be seen that the both values obtained are far from being realistic for the whole range of frequencies and attenuation coefficients considered, and mainly overestimates the active power. This could be expected since the two configurations correspond to waves traveling in open-space, so that in these cases, the real part of the source impedance is merely a measure of the mechanical energy transported irreversibly far from the source, and not of the power dissipated in the medium.

Finally, figure 8 displays the ratio of the dissipated power in the cooling liquid to the total active power sent by the sonotrode for α = 0.05, 1 and 5. It is seen that in the whole range of frequencies, a noticeable part of the energy enters the cooling jacket and is dissipated inside. Besides, it can be seen that the curves in Fig. 8 cross some resonances (I and IV) without varying noticeably. This is an indication that such resonances are the result of a collaborative coupling between the working liquid and the liquid in the cooling jacket, mediated by the glass wall.

Imaginary source impedance and field structure

The above results show that neither the real source impedance of the plane traveling wave, nor the one of the baffled piston are adequate to calculate the active power dissipated in the sono-reactor. Nevertheless, it is instructive to have a look at the imaginary part of the source impedance.

Figure 9 displays (Z)/ρ l c l , for various attenuation coefficients (thin lines). Near resonance, this quantity presents a stiff jumps, which corresponds to a π phase change (in the non-dissipative case) of the pressure field as resonance is crossed. This behavior is well-known for harmonic oscillators and can be easily demonstrated in the case of plane standing waves [START_REF] Morse | Theoretical acoustics[END_REF].

Also represented in Fig. 9 is the imaginary impedance for a baffled piston (Eq. 30, thick solid line). It is seen that this line cuts more or less the imaginary impedance curves in two symmetric parts, and even merges with these curves for frequencies far from resonances (see the interval . This means that, in the latter case, the acoustic field is not so far from the one of a baffled piston emitting in open-space (which can actually be seen in Fig. 3), and demonstrates the weak effect of the boundaries for such frequencies. This is a confirmation that far from resonance, the location of the maximal pressure amplitude near the source is not only due to attenuation, but is mainly a consequence of the spreading of the wave in space, as if there were no boundaries. The simulation presented here therefore demonstrate that such a field struc- ture, commonly encountered with ultrasonic horns, is a consequence of the detuning of the system, and could be avoided by adjusting either the frequency or the geometry, as demonstrated experimentally in Ref. [9].

Finally , for a plane traveling wave, it is readily seen from (29) that the imaginary impedance would be zero (thick dotted line in Fig. 9), which is far from being realistic in this case. Also represented is the imaginary impedance for a spherical wave (thick dashed line), which appears as an overestimation of (Z) far from resonances. We take this opportunity to note that, rigorously speaking, owing to the value of ka near unity for classical ultrasonic horns at 20 kHz (here ka = 1.3), the wave structure far from resonance is rather the one of a baffled-piston than a spherical wave field (for which ka → 0), as sometimes mentioned in the literature [START_REF] Yasui | [END_REF].

Vibrations of the solid

It is instructive to compare the acoustic energy stored in the glass walls to the one stored in the liquid. Figure 10 displays the ratio of both quantities, calculated respectively from Eqs. (31) and (32). It can be seen they are of the same order of magnitude for any frequency and attenuation coefficient, which demonstrates again the role of the walls deformation.

An interesting feature appears on Fig. 10: some resonance peaks are still visible on the curves: peak IV completely disappears, which means that this resonance is mainly governed by the liquid. Conversely, peak V is considerably magnified, which attests that it is rather governed by the solid. This is confirmed by a calculation of the normal modes of vibration of the solid walls alone, which exhibits a resonance frequency near 20800 Hz. Peaks I and II are smoothed but still visible, and are found to result from a complex coupling between the working fluid, the cooling fluid and the solid vibrations. Besides, the dependence of the curves on attenuation can be explained as follows. Far from resonance, the boundaries receive a minor part of the emitted wave, mainly because the wave is diverging (see Fig. 4 and section 5.4) and has considerably decreased as it reaches the solid boundaries. The energy stored in the latter is therefore weak in this case, and not very sensitive to attenuation. Conversely, near resonance, some synergistic effect between the solid and the liquid vibrations can come into play. The higher the role played by the walls in the resonance, the higher the energy they store depends on the liquid attenuation factor, because the energy they receive from the source is mediated by the liquid.

Sensitivity to liquid properties

In order to emphasize the influence of the operating conditions on the behavior of the system, Figure 11 compares the response curves obtained when the density of the cooling liquid is decreased from 1000 to 950 kg.m -3 . The five resonance peaks are still visible, but peaks I, IV and V are slightly shifted toward higher frequencies. This suggests that even a slight change in density in the cooling jacket, for example produced by a temperature variation of the cooling liquid, can produce a noticeable detuning of the system, if the latter does not automatically selfadjusts its frequency.

Conclusion

Simulations of a real sono-reactor have been carried out, taking into account the vibrations of the walls. The frequency was systematically scanned in the range available for the generator used, in order to detect the resonances of the global system.

It has been shown that the simple hypothesis classically used at the liquid boundaries, either soft or hard wall, definitely yields unrealistic results, and could not anyway predict the sound field in the cooling jacket. Indeed, the vibrations of the solid walls induce a significant coupling between the reacting liquid and the cooling liquid. Simulation predicts that a noticeable energy fraction can be transmitted to the liquid in the cooling jacket, which may therefore undergo cavitation. This opens the possibility to design sono-reactors in which the working liquid is insulated from the sonotrode by a solid wall, which constitutes an interesting feature for electrochemical reactions, where electrical currents leaking through the metallic sonotrode should be avoided [9].

Besides, it has been shown that, exciting the sonoreactor near a resonance frequency may present two advantages. The first one lies simply in the large amplitude of the acoustic pressures obtained. The second is the possible localization of pressure antinode far from the emitter, which could avoid the well-known issue of cavitation-induced erosion of the sonotrode. This latter feature may however disappears for too strong attenuation.

The concepts of acoustic energy, intensity, source impedance, active and dissipated powers have been clarified on a rigorous development, and illustrated by simulation. Disregarding the acoustic energy transmitted to the surrounding air by boundaries, the active power is nothing else than the dissipated power in the liquid, averaged over one period. This was checked by simulation and indeed represents a good test to check the validity of the numerical results. Besides it was shown that the approximation of plane traveling wave, commonly used in order to link the amplitude of the source to the active power, does not make sense to calculate the latter. Neither does the more refined model of baffled piston, but it yields a reasonably good approximation of the field structure for some frequencies far from resonances.

Finally, a slight variation of the liquid properties shifts noticeably the resonance peaks, which demonstrates the interest of using a self-adjusting generator.

Linear acoustics therefore yields qualitatively interesting results in studying and designing sonoreactors, provided that the vibrations of the solid are accounted for. The main unknown remains the attenuation factor, strongly linked to the bubble repartition in the liquid, which in turn depends on the shape of the acoustic field [START_REF] Mettin | Sixth International Symposium on Cavitation -CAV2006 (paper 75[END_REF]. Coupling both phenomena in a FEM model including the main physics of cavitation remains to be done, and constitutes a challenge for further studies. Among other issues, such a model would allow to assess the stability of the resonances when the bubbles come into play, and the possibility to tune sono-reactors in order to get the cavitation field far from the emitter.

(i) the maximum amplitude is obtained for ω = (ω n -α 2 c 2 l ) 1/2 , so that attenuation yields a slight shift of the resonance peak, compared to the non-damped case for which ω n = ω, (ii) the maximum value of the peak scales as α -1 , (iii) the full-width at half maximum (FWHM) of |P |, that is the width between the frequencies for which |P | = 1 2 |P | max , is ∆f = αc l /π √ 3. By repeating the simulation with idealized hard boundaries for different attenuation coefficients (corresponding to the dashed line of Fig. 2 which was obtained in the non-dissipative case), the above three results were checked successfully. This is not the case however when coupling with the solid is taken into account, and the main discrepancies are described hereafter.

Concerning result (i), the shift in the peak due to attenuation is not visible on the curves of Fig. 5 because they are to small to be noticed at the scale of the figure. Refined simulations show that it indeed exists, but is smaller than the theoretical value calculated from (B.1), which amounts to about 35 Hz in the most dissipative case (α = 5 m -1 near f 20000 Hz). The discrepancy is attributed to the coupling with the solid walls, for which the simple result (B.1) obtained for idealized boundaries may become invalid.

It could be checked that the result (ii) above is approximately true for the main peaks of Fig. 5, and the small discrepancy can also be attributed to the coupling with the solid walls.

Finally, concerning (iii), the width of the resonance peaks was found lower than their theoretical value, here again probably because of the coupling with the solid.

Fig. 1 .

 1 Fig. 1. Geometry studied. Only one-half of the reactor is drawn.

Fig. 2 .

 2 Figure2displays the mean acoustic pressure in the reactor as a function of frequency, with zero attenuation, in the following cases: cooling jacket filled with water (thick solid line), cooling jacket empty (thin solid line), internal walls considered as hard boundaries (dashed line), or as soft boundaries (dash-dotted line). In the two latter cases, the vibration of the glass and Teflon walls are not accounted for, and it is clearly seen that the response curves differ significantly from the ones obtained by properly taking into account the walls vibration.The peaks observed in the response curve correspond to global resonances of the whole mechanical system formed by the liquid coupled with the solid walls. In absence of dissipation, the amplitude of these peaks should be infinite, which would traduce an infinite accumulation of acoustic energy in the system. Thus, the finite maxima of the peaks appearing on Fig.2are only due to the discretization of the frequency axis. It is interesting to note that the response curves, with the cooling jacket respectively filled and empty, present approximately the same set of peaks, shifted in frequency, except for peak II which is specific to the filled cooling jacket, which already demonstrates the strong influence of the operating conditions on the response of the system.We now focus on the case where the cooling jacket is filled with water (Fig.2, thick solid line). Figure 3 displays the instantaneous acoustic pressure p(r, z, t) at ωt = 0, π/2 and π near peak IV (f =

Fig. 3 .

 3 Fig. 3. Pressure field p(r, z, t) and wall deformation near peak IV (f = 20610 Hz) at times ωt = 0, ωt = π/2, ωt = π. The wall displacement is magnified 100 times.

Fig. 4 .

 4 Fig. 4. Pressure field p(r, z, t) and wall deformation for f = 19500 Hz at times ωt = 0, ωt = π/2, ωt = π. The wall displacement is magnified 20000 times.

Fig. 5 .

 5 Fig. 5. Mean acoustic pressure peak amplitude Pav vs. frequency, for various attenuation coefficients. Thin solid line: α = 0. Dashed line: α = 0.05. Dash-dotted line: α = 0.5. Dotted line: α = 1. Thick solid line: α = 5.

α 1 Fig. 6 .

 16 Fig. 6. Repartition of the peak acoustic pressure |P | in the reactor for f = 20610 Hz and increasing attenuation coefficients. The color scales are different for each image.

Fig. 7 .

 7 Fig. 7. Real part of the dimensionless source impedance (Z)/ρ l c l , where Z is calculated from (27), for various attenuation coefficients. The active power sent through the emitter can be deduced from Eq. (28). Solid line: α = 0.05. Dashed line: α = 0.5. Dotted line: α = 1. Dash-dotted line: α = 5. The thick line represents the real impedance for a baffled piston emitting in open space Eq. (30). Thick dotted line: same for plane traveling waves ( (Z)/ρ l c l = 1).

Fig. 8 .

 8 Fig. 8. Ratio of the power dissipated in the cooling jacket to the active power sent by the sonotrode. Solid line: α = 0.05. Dashed line: α = 1. Dash-dotted line: α = 5.

Fig. 9 .

 9 Fig. 9. Imaginary part of the dimensionless source impedance (Z)/ρ l c l , where Z is calculated from (27). Solid line: α = 0.05. Dashed line: α = 0.5. Dotted line: α = 1. Dashdotted line: α = 5. Thick solid line (Z)/ρ l c l calculated for a baffled piston emitting in open space Eq. (30). Thick dashed line: same for a spherical wave. Thick dotted line: same for plane traveling wave ( (Z) = 0).

Fig. 10 .

 10 Fig. 10. Ratio of the acoustic energy stored in the glass walls to the acoustic energy stored in the liquid (including both the working liquid and the liquid in the cooling jacket ). Solid line: α = 0.05. Dashed line: α = 1. Dash-dotted line: α = 5.

Fig. 11 .

 11 Fig. 11. Response curves obtained for a density of the liquid in the cooling jacket 1000 kg.m -3 (solid line) and 950 kg.m -3 (dashed line).

Acknowledgements

The authors thank the Generalidad Valenciana (Spain) for its financial support under Projects AINV07/044 and AE/07/079 and also to COST D32 for the framework. JK gratefully thanks the Ministry of Education, Youth and Sports (MSMT) of the Czech Republic for financial support (grant number 1P05OC074).

Appendix A. Equation of energy conservation

For any mechanical system, the theorem of kinetic energy states that dK dt = P ext + P int , (A.1)

where P ext and P int are the powers of external and internal forces respectively. In the case of an arbitrary volume V of fluid in motion, this equation reads where Π diss > 0 is the dissipation function, and n is the normal unit vector pointing outward the liquid.

In the case of acoustic waves, using the mass conservation equation Eq. ( 1) to express ∇.v in the third integral, we obtain

which reverts to equation (19) in the case of nondissipative waves (Π diss = 0). If mono-harmonic dissipative waves are assumed, the acoustic energy e a reads, with expressions (3), ( 4):

+ oscillating terms, (A.4) so that the left-side of Eq. (A.3) is zero on average over one acoustic period, and the period-average of the acoustic intensity pv is found to read

where denotes the real part of a complex number, so that (A.3) finally becomes

This equation can be applied to the volume V of liquid contained in a sono-reactor, bounded by the emitter surface S sonotrode on one hand, and by passive surfaces S boundaries (either internal surface of a solid wall or free air) on the other hand. In this case, Eq. (A.6) can be written as:

If the boundaries do not transmit any wave to the external air, P bound is zero. This is the case for example for infinitely soft (P = 0) or infinitely rigid boundaries (V = 0). In the more complex case where S boundaries is the internal wall of the solid vessel, its external wall being unconstrained, it can be shown that P bound is still zero for a perfectly elastic solid material.

Besides, the dissipated mechanical energy Π diss can be deduced from the knowledge of the attenuation coefficient α (see Eq. ( 26)). To see that, we start by rewriting the Helmholtz equation ( 7) with a complex wavenumber, as

Extracting ∇P from Eq. ( 6), the set of Eqs. ( 5), ( 6) becomes:

where denotes the imaginary part of a complex number. Multiplying the first by P * and the complex conjugate of the second by V, summing, and taking the real part, we obtain:

Integrating on a volume V and using the divergence theorem, we get:

(A.12) Identifying this equation with equation (A.6), the average dissipated power density is therefore seen to read

Further using the definition of the attenuation coefficient (26), the latter expression can be written as

and is clearly positive. Replacing k r by ω/c l , the dissipated power density finally becomes

The quantity |P | 2 /(ρ l c l ) is sometimes erroneously termed as "intensity". We emphasize that this is true only for plane traveling wave, but is definitely not in other cases. Expression (A.15) remains however valid, and allows to estimate |P | by measuring Π diss , generally equating the latter to ρ l C p (dT /dt) 0 .

Appendix B. Acoustic resonance of closed domains

Theoretical acoustics [START_REF] Morse | Theoretical acoustics[END_REF] states that a domain of fluid closed by rigid (∂P/∂n = 0), soft (P = 0) or reactive boundaries (∂P/∂n + βP = 0) exhibits a set of acoustic resonances at definite frequencies ω n , where ω n /c l are the eigenvalues of the Helmholtz equation [START_REF] Lauterborn | [END_REF]. Near such a resonance (ω ω n ), the complex pressure amplitude is proportional to

where k = ω/c l -iα is the complex wavenumber, accounting for attenuation. Thus, near a resonance, the problem reverts to the classical problem of damped resonant systems, which classical properties can be readily checked from Eq. (B.1), and are recalled hereafter:

List of Figures