Dynamics of cavitation bubbles in acoustic field near the rigid wall

Ocean Engineering 109 (2015) 507–516

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Dynamics of cavitation bubbles in acoustic field near the rigid wall Ye Xia,n, Yao Xiongliangb, Han Ruib a b

Marine Design & Research Institute of China, Shanghai, China Harbin Engineering University, Harbin, China

art ic l e i nf o

a b s t r a c t

Article history: Received 10 April 2015 Accepted 21 September 2015

Dynamics of multiple cavitation bubbles near the rigid wall subjected to the incident acoustic wave with different arrangements is researched in this paper. Wave equation is used to approximate the prophase and anaphase motion of cavitation bubble to consider the compressibility of fluid. The modified Bernoulli equation is deduced as the boundary condition for cavitation bubble in acoustic field. Influences of compressibility, acoustic wave frequency and amplitude, incident angle, and the bubble arrangements on the dynamics of cavitation bubbles are presented. Numerical results show that the compressibility of fluid will impact the intensity of motion. The dynamics of bubbles is affected by the incident amplitude and frequency, and limited by the rigid wall. The different incident angles and arrangements of bubbles will change the acting pressure phase and interaction between bubbles and rigid wall, which will lead to the motion advancement or delay. The radius difference will increase with the large phase difference and the strong rejection. The jet will be oblique to the incident direction, as well as the rigid wall and bubbles around. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Multiple cavitation bubbles Acoustic wave Rigid wall Compressibility

1. Introduction Cavitation bubble in acoustic field will oscillate periodically due to the pressure gradient in space. During the oscillation, the high pressure environment is appeared inside the bubble. At the same time, noise is radiated induced by the bubble motion. And for some cases, there will be light emitting from bubble, which called sonoluminescence (Crum and Reynolds, 1985; Gaitan et al., 1992, Gaitan and Crum 1990) as well known. Since the severe contraction, the compressibility of fluid is significant to the dynamics of bubble. Prosperetti and Lezzi (1986a, 1986b) made the correction to consider compressibility of fluid for bubble dynamics. Otherwise, there are many improved equations for bubble motion with compressible correction (Keller and Miksis, 1980; Flynn, 1975a, 1975b; Lastman and Wentzell, 1981, 1982). Bubble dynamics near the wall of structure underwater has been researched in many fields ( Zhang et al., 2013a, 2013b, 2015a, 2015b). Due to the existence of wall will lead to the non-spherical deformation for cavitation bubble such as jet or splitting, the dynamics of cavitation bubble is not able to be solved analytically. Therefore, the numerical method is adopted to simulate the complex motion of bubble (Park, 2008; Qiu et al., 2008; Blake and Gibson, 1981; Wang, 2005; Liu, 2014; Zhang and Liu, 2015). As well n

Corresponding author. Tel.:+86 13654508241. E-mail address: [email protected] (X. Ye).

http://dx.doi.org/10.1016/j.oceaneng.2015.09.045 0029-8018/& 2015 Elsevier Ltd. All rights reserved.

known, BEM (Blake and Gibson, 1981; Wang, 2005; Liu, 2014; Zhang et al., 2015a) is able to reduce the problem dimension using the boundary integral equation with the advantages in precision and efficiency. There are many presented results about bubble dynamics using BEM. Michael (2007) analyzed the motion stability of acoustic cavitation bubble with classical BEM. And the heat exchange, diffusive equilibrium and the chemical reaction during the motion are considered. However, the classical BEM is based on the potential theory to solve the Laplace equation without take the compressibility of fluid into account. To overcome this disadvantage of classical BEM, Wang (Wang and Blake, 2010, 2011) utilized perturbation expression to the velocity potential which satisfies the Laplace equation and wave equation respectively. Finally, the corrected Bernoulli equation considering compressibility was obtained by the matching of velocity potential in space. In this paper, in order to consider the compressibility of fluid, the bubble motion is divided into prophase and anaphase in timedomain, and approximated by wave equation respectively. Then the approximations in the two time phase are matched with the help of Laplace transformation. And the corrected boundary integral equation considering compressibility is obtained. Although using the BEM, the method is quite different from that presented by Wang and Blake (2010, 2011). With BEM using the corrected boundary integral equation, Ye and Zhang et al. have presented the dynamics of single cavitation bubble in free field (Yao et al., 2013) or near rigid wall (Ye et al., 2014a, 2014b) subject to incident

508

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

Fig. 1. Numerical model.

acoustic wave with axisymmetric model or three-dimensional model. However, the axisymmetric model is not able to solve the multiple cavitation bubbles with arbitrary arrangement and arbitrary incident angle. Thus, the three-dimensional model is utilized in this paper. The dynamics of multiple cavitation bubbles with different arrangements near rigid wall in sound field are analyzed. And the influences of compressibility, incident angle, incident frequency and incident amplitude on the cavitation bubbles are presented.

1 ,  ∬SðtÞ φn ðr q ; tÞdS: c

ð5Þ

,

where ξ is solid angle; n q is unit outward normal vector at the source point; r pq is the distance between field point and source point. For convenience, Eq. (5) is also rewritten as follow: ,

,

,

cΛφ þ cℑφ n  Θφ nt ¼ 0:

ð6Þ

2. Basic theory

where Λ;ℑ;Θ are the coefficient matrixes. With the Laplace transformation, Eq. (2) and Eq. (6) are matched, and the boundary integral formula is deduced in the compressible fluid expressed as follow (Zhang et al., 2013c):

2.1. Boundary integrand equation in the compressible fluid

φ tt þcΑλ t þ Βc2 φ þ cφ nt þc2 Πφ n ¼ 0:

ð7aÞ

We use wave equation to describe the perturbation fluid induced by the cavitation bubble:

A ¼ Πþκ

ð7bÞ

c2 ∇2 φ ¼ φtt ;

B ¼ Πℑ  1 Λ

ð7cÞ

ð1Þ

to consider the compressibility of fluid. Where, φ is the velocity potential induced by cavitation bubble; c is the sound speed in fluid, here we assume the sound speed is a constant. The fluid is assumed to be irrotational and inviscid. Since the small size of cavitation (for example, 80 μm) and short time scale (for example, 4 μs) considered in this paper, the gravitation is neglected. In the early stage of bubble motion, the cavitation bubble is near sphere. Therefore, the prophase approximation for the motion of cavitation bubble is (Zhang, 2013c):

φt ðR; t Þ þ cκφðR; t Þ þcφn ðR; t Þ ¼ 0:

ð2Þ

where κ and R are the surface curvature and radius, respectively. Solving Eq. (1) with its Green function G,
 r 1 δ t  τ  pq ; ð3Þ G¼ 4π r pq c1 we obtain:
 Z , φ r p; t 

1

Sðτ Þ

0

þ ,



Z

,

1 c

,

,

φðr q ; τÞGn Gφn ðr q ; τÞ  2 φðr q ; τÞGτ

 1 , Gφτ ðr q ; τ Þ dSdτ ¼ 0: 2 c

ð4Þ

,

where r p and r q are the location vector of field point and source point, respectively; t and τ are the time for field point and source ,

point, respectively; v is the velocity of bubble surface. In order to obtain the quantities on the surface, the filed point p and source point q are all arranged on the bubble surface. When the motion time t is sufficiently large to ensure r pq =ðtcÞ{1, the term of order 1/c is hold only, i.e. the terms in Eq. (4) with 1/c2 will be ignored. Taylor expansion is used for Green function G, and also the terms upon order of 1/c2 will be neglected. Then, we get the anaphase approximation for the bubble motion (Zhang et al., 2013c): , ðr p ; tÞ ¼

ξφ

, , r pq U n q  ∬Sðt Þ 3 r pq

, ðr q ; tÞdS  ∬Sðt Þ

φ

,

φn ðr q ; tÞ r pq

dS

,

,



,



,

,

1

Π ¼ ℑ  1 Λ  κ I  ℑ  1 Θℑ  1 Λ

ð7dÞ

In Eq. (7), κ is the local curvature matrix; I is identity matrix. When the sound speed is set to be infinite, Eq.(7a) is simplified to the boundary integral formula for incompressible fluid as the classical BEM to solve Laplace equation. 2.2. Boundary condition on the surface of cavitation bubbles Multiple cavitation bubbles in equilibrium state are released near the rigid wall. A spatial traveling plane wave PI incomes to the rigid wall with various angle θi generating the reflecting wave PR (Fig. 1). For convenience, and without affecting the discussion, the incident wave is set to be paralleled to x–z plane. Although the frequency of incident wave (incident frequency) is high (for example, 250 KHz), the wave length is still much greater than the maximum size of bubble during its motion. Therefore, we assume that the propagation of plane wave will not be impacted by the existence of bubbles, but the motion of bubbles will be affected by the plane wave, i.e. the bubbles are one-way coupled to the incident plane wave. To solve the motion of cavitation bubbles near the rigid wall, the mirror method (Morse and Ingard, 1968) is adopted in this paper. The stable acoustic field formed at one side of rigid wall with incident wave PI and reflecting wave PR can be expressed as:   P I ¼ P a cos ωt þ θ0  kz cos θi  kx sin θi : ð8aÞ   P R ¼ P a cos ωt þ θ0 þ kz cos θi  kx sin θi :

ð8bÞ

    P st ¼ P I þ P R ¼ 2P a cos kz cos θi cos ωt þ θ0  kx sin θi :

ð8cÞ

where Pa is the amplitude of the incident wave; k is wavenumber; ω ¼2πf is the round-frequency, f is the incident frequency, and

,

θ0 is initial phase. We set the total velocity vector u

X. Ye et al. / Ocean Engineering 109 (2015) 507–516 ,

in the fluid field as: ,

,s

u ¼ ∇φ þ u :

ð9Þ

s ,

where u is the velocity induced by the acoustic wave, which ,s

satisfies ρ1 u t þ ∇P st ¼ 0 (Morse and Ingard, 1968). Navier–Stokes equation in compressible inviscid fluid could be simplified to: Z dP , , , u t ¼  ðu U ∇Þu  ∇ : ð10Þ

ρ

Where, P including three parts: pressure induced by bubbles, acoustic field pressure Pst and the pressure P1 without perturbation. Noticing that in the irrotational fluid: 1 , 1 , , , , , ðu U ∇Þu ¼ ∇j u j 2 u  ð∇  u Þ ¼ ∇j u j 2 ; 2 2

ð11Þ

and assuming the density ρ in the fluid is only the function of pressure, which can be expressed as (Morse and Ingard, 1968): P  P 1 ¼ c2 ðρ  ρ1 Þ;

ð12Þ

to consider the compressibility of fluid. Where, ρ1 is the density in fluid without perturbation. Substituting Eqs. (9), (11) and (12) into Eq. (10), with the Taylor expansion we can get: s s ∂ 1 P  P1 , , ð∇φ þ u Þ ¼  ∇j ∇φ þ u j 2  ∇ð Þ: ∂t 2 ρ1

ð13Þ

And Eq. (13) is able to be rewritten as follow with further substitution: 

2 ,s P  P st  P 1 1 1,s ∇ þ φt þ ∇φ þu U ∇φ þ u j 2 ¼ 0: ð14Þ 2 2 ρ1 Since P  Pst  P1 and velocity potential φ induced by bubble will decrease to zero at infinite, the Bernoulli equation in Lagrangian frame for the fluid with cavitation bubble and traveling wave is:  ρ P ¼ P st þ P 1  ρ1 ℘φ þ 1 j∇φ2 : ð15Þ 2 where ℘ is the material derivative to time t. The pressure P at the bubble surface satisfies:
γ V0 P ¼ Pv þ P0 M κ : ð16Þ V where P0, V0 is the initial inner pressure and volume of cavitation bubble; Pv is saturated vapor pressure; γ is the gas specific heat, and is equal to 1.4; M is the surface tension coefficient, and is equal to 0.0728 N m  1. The kinematic material boundary condition on the bubble surface is: ,

509

,

℘l ¼ ð∇φ þ u s ÞjSðt Þ :

ð17Þ

where l is the location vector on the bubble surface. Solving Eqs. (7), (15), (16) and (17), we can get the dynamics of cavitation bubble near the rigid wall driven by stable acoustic field in the compressible fluid. The detail flow is shown as Fig. 2. After setting the initial conditions including the bubble location, inner pressure, velocity potential on the bubble surface, etc., the normal velocity of bubble surface is obtained by solving boundary integral Eq. (7), and the material velocity of bubble surface can be calculated using normal velocity and velocity potential as described by Zhang et al. (2001). To avoid the mesh gather on the discrete bubble surface, we utilize the EMT (Wang, 2003) to correct the material velocity. And the corrected material velocity is used to update the location of bubble surface. The bubble surface need to be smoothed using the method presented by Sun (2007). Then, we get the new velocity potential using Bernoulli Eq. (15) to be the initial condition for next time step. And the velocity potential need to be smoothed as well. The time step can be set as Wang and Blake (2010, 2011). In additional, due to the arbitrary incident angle, the three dimensional numerical model is utilized in this paper.

3. Numerical results and discussion 3.1. Numerical verification We calculate a test to prove the correction of the method presented in this paper. A single cavitation bubble is released in the free-field. Since the existence of sound field with constant sound speed, the results obtained by the boundary integral formula (7) are compared to that obtained by Keller–Miksis equation (Keller and Miksis, 1980) without consider the viscosity. The radius calculated by boundary integral equation is defined as the equivalent radius of sphere who has the same volume with bubble. The results are shown in Fig. 3, and we can find that the method presented in this paper is able to solve the dynamics of cavitation bubble driven by acoustic wave precisely. Otherwise, the comparison between compressible and incompressible results indicates that the influence of compressibility on the bubble motion subjected to acoustic wave is obvious besides the first cycle. With the consideration of compressibility, the time for once expansion and contraction is shortened, and the oscillation amplitude is weakened. In order to prove the method presented in this paper is able to solve the motion of multiple bubbles near the rigid wall, we make the comparison between experiment results (Ni, 2012) and numerical results as shown in Fig. 4, and the lines in figure are the outline of bubble from numerical calculation. In this experiment, the maximum radius of bubble is 6.93 mm, and the initial center of Boundary integral equation

Velocity potential on bubble surface

Initial condition

Normal velocity on bubble surface

Bernoulli equation Updating the velocity potential on the surface

Updating the bubble surface location

Material velocity on the bubble surface EMT

Smooth technology Fig. 2. Numerical simulation flow chart.

510

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

bubble is 7.48 mm apart from the two plates respectively. In the numerical model the initial radius of bubble is 1.03 mm, and the initial inner pressure of bubble is 9.8 MPa. The bubble will expand, split, contract and form jet between two parallel rigid plates. We find that the bubble calculated by numerical method is matched well with the bubble in experiment.

3.2. Double cavitation bubbles Cavitation bubbles A and B in equilibrium state are released near the rigid wall, and the bubble A is arranged at (0, 0, zA), while the bubble B is arranged at (xB, yB, zB). The coordinates of cavitation bubbles are non-dimensional with initial radius R0A ¼R0B ¼ 15 μm. Firstly the dynamics are analyzed with setting zA ¼zB ¼1.5, xB ¼3, yB ¼0, i.e. the bubbles are parallel with the rigid wall. The incident frequency, incident amplitude, initial phase are set as: f ¼224.48 KHz, Pa ¼1.8 atm, θ0 ¼ π, the incident angle is π/4. The velocity and pressure distribution in near-field of the double cavitation bubbles are calculated by method presented in by Yao et al. (2013), Ye et al. (2014a) as shown in Fig. 5. Bubble B offsets to positive x-axis compared to bubble A. As Eq.(8), therefore the phase of sound pressure acting on bubble B delays to that on bubble A, which leads to the motion lag of bubble B.

Fig. 3. Verification for the result from BEM. ρ1 ¼ 998 kg/m, c1 ¼1500 m/s, P1 ¼1 atm, Pv ¼2338 Pa, R0 ¼25 μm, Ps ¼ 1.2 P1, f¼ 80 KHz, the acoustic wave is set as Eq. (8) with x ¼z ¼ 0 and θ0 ¼π. (the analytical results are obtained from Keller– Miksis equation; the compressible and the incompressible results are calculated with Eq. (7)).

At the prophase of motion, since the sound pressure acting on the bubbles is negative, they all expand, as shown in Fig. 5(a). The inner pressure is much greater than that in the fluid. After reaching the maximum radius, they begin to contract. As shown in Fig. 5(b), the two bubbles are all at the prophase of contraction at this moment. The inner pressure is lower than that in the fluid, and the velocity in fluid points inside the two bubbles respectively. Due to the motion lag of bubble B, the velocity around bubble A is higher since the advanced contraction. As shown in Fig. 5(c), the two bubbles are going to generate the jets. Since the motion lag of bubble B, the volume of bubble A is less, which leading to the higher inner pressure. Also, due to the bubble A forming the jet earlier than bubble B, the pressure and velocity around is higher at this moment. As shown in Fig. 5(d), bubble A is at the anaphase of jet stage. With the influence of rigid wall and the vicinity cavitation bubble, the jets are attracted each other while pointing to the rigid wall. At this moment, it indicates clearly that since the phase lag of the sound field, the jet of bubble B delays to that of cavitation A obviously. But the velocity and pressure near the jet of bubble B are all higher than that around the bubble A now, i.e. the jet shooting of bubble B is stronger. Fig. 6 shows the radius of two bubbles and their difference (ΔR ¼RB  RA). And the radius is the equivalent of sphere with the same volume with the bubble. Fig. 7 shows the jet tip velocity of two bubbles. At the prophase of contraction, since the contraction state of bubble A is advanced than bubble B, the jet tip velocity of bubble A is higher. And the radius difference is increased until the jet stage. Then, at the anaphase of jet stage of bubble A, since motion lag and the attraction effect between the two bubbles, the jet tip velocity of bubble A will not increase any more, while the jet tip velocity of bubble B is still increased. And finally the jet tip velocity of bubble B is higher than that of bubble A, which will accelerate the contraction of bubble B. Therefore, the radius difference is decreased. In additional, the time for negative pressure acting on the bubble B is longer leading to the sufficient expansion. Thus the volume of bubble B is greater, and ΔR 40. 3.2.1. Compressibility Now, we discuss the influence of compressibility of fluid on the motion of cavitation bubble. The incident frequency, incident amplitude, initial phase are set as: f ¼224.48 KHz, Pa ¼1.8 atm, θ0 ¼ π, the incident angle is 0. Since the incident angle is 0, there is no motion delay between two bubbles. Thus, the motion characteristics of two bubbles are similar besides the direction of jet. Fig. 8 shows the comparison for

Fig. 4. Comparison for experiment (Ni, 2012) and numerical simulation. Time of the results is t¼ 0.10 ms, 0.75 ms, 1.30 ms, 1.55 ms, 1.60 ms, 1.65 ms, 1.70 ms, 1.75 ms. White lines mean the numerical simulation results.

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

0.036 0.034 0.031 0.029 0.026 0.024 0.021 0.018 0.016 0.014 0.011 0.009

z x

511

t =1.15μs

t =2.08μs

t =3.26μs

t =2.97μs

Fig. 5. Pressure and velocity distribution (f ¼224.48 KHz, Pa ¼ 1.8 atm, θ0 ¼ π, θi ¼π/4, non-dimensional coordinate of bubble A is (0,0,1.5), and non-dimensional coordinate of bubble B is (3,0,1.5), (a)–(c) is at the x–z plane, unit of pressure is MPa, the bubble at the left in the figures is bubble A, and at right is bubble B), (a) t¼ 1.15 μs, (b) t ¼2.08 μs, (c) t¼ 2.97 μs, (d) t¼ 3.26 μs.

Cavitation A Cavitation B ΔR

24 22

120

1

100

18 0.8 16 0.6

14

0.4

12

CavitationA CavitationB

80 60 40

0.2

10 8

Velocity/ms

-1

140

1.2

ΔR

R/μm

20

1.4

20 0

0.5

1

1.5 t/μs

2

2.5

3

0

Fig. 6. History of radius and its difference (f¼ 224.48 KHz, Pa ¼ 1.8 atm, θ0 ¼ π, θi ¼ π/ 4, non-dimensional coordinate of bubble A is (0,0,1.5), and non-dimensional coordinate of bubble B is (3,0,1.5)).

radius and jet tip velocity in compressible fluid and incompressible fluid for one bubble. With the consideration of compressibility, the maximum radius at first expansion is less than that in incompressible fluid, and the time reaching maximum velocity of jet tip is advanced. However, the influence of compressibility is not very obvious since the just once oscillation and the perturbation is not propagated to very far field with finite velocity. 3.2.2. Incident amplitude and frequency Here, the influence of incident frequency on the motion of cavitation bubble is discussed. The incident amplitude, initial phase are set as: Pa ¼1.8 atm, θ0 ¼ π, the incident angle is 0. Fig. 9 shows the radius and jet tip velocity of cavitation bubble driving by different incident frequencies. Since the low incident frequency will lead to slow motion state transformation, the bubble is able to expand sufficiently. Thus, the maximum radius at expansion is increased. And the time for once expansion and contraction is longer. However, since the existence of rigid wall, the motion of bubble will be limited, the greater maximum radius

0

0

0.5

1

1.5

2

2.5

3

3.5

t/μs Fig. 7. Jet tip velocity (f ¼224.48 KHz, Pa ¼1.8 atm, θ0 ¼ π, θi ¼ π/4, non-dimensional coordinate of bubble A is (0,0,1.5),and non-dimensional coordinate of bubble B is (3,0,1.5)).

will cause great radius at jet stage. Due to the large amplitude of expansion for low frequency, the contraction will be severer, which leads to the higher peak value of jet tip velocity. The influence of incident amplitude on the bubble motion is shown in Fig. 10. The incident frequency, initial phase are set as: f¼ 224.48 KHz, θ0 ¼ π, the incident angle is 0. With the increment of incident amplitude, the maximum radius at expansion stage will increase as well and the jet tip velocity is enlarged due to the severe contraction. Since the limitation of rigid wall, the larger maximum radius will cause large radius at jet stage. 3.2.3. Incident angle Fig. 11 and Fig. 12 show the radius difference and the shape comparison at jet stage with different incident angles, respectively. Since the no-vertical incidence, there is a motion lag between the two bubbles. With the enlargement of incident angle, as Eq. (8), the phase lag is also increased, which leads to the more obvious motion lag. Therefore the radius difference is increased. Also,

512

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

25 Compressible Incompressible

250 Jetvelocity/ms -1

Radius/μm

20

15

10

200 150 100

Compressible Incompressible

0

0.5

1

50 0

1.5 t/μs

2

2.5

3

2

2.5

3

3.5

t/μs

Fig. 8. Comparison for compressible and incompressible model, (a)Radius, (b)Jet tip velocity.

30

0.8 1 1.2

140 Jet tip velocity/ms-1

25 Radius/μm

160

0.8 1 1.2

20 15 10

120 100 80 60 40 20

5

0

1

2 t/μs

3

4

0

0

0.5

1

1.5

2 t/μs

2.5

3

3.5

4

Fig. 9. Characteristics of cavitation bubble with different incident frequencies, the legend in figure means the radio to the base frequency f¼ 224.48 KHz, (a) Radius, (b) Jet tip velocity.

as shown in Fig. 12, bubble A is at the end of jet stage, but the jet of bubble B is still developing for the no-vertical incidence cases. And with the larger incident angle, the jet stage of bubble B is earlier. 3.2.4. Different arrangement We set the initial central distance Dc as 2.5, 3 and 3.5, which is non-dimensionalized by initial radius. The radius difference and shape comparison at jet stage with difference initial central distances are shown in Fig. 13 and Fig. 14, respectively. Here, the traveling wave is still no-vertical incidence. And with the increment of the initial central distance, the phase lag of incident wave is obvious, which leads to the apparent motion lag of bubble B. Therefore, as shown in Fig. 13, the radius difference is increased. As shown in Fig. 14, with the reduction of initial central distance, the motion lag is decreased, which leads to more sufficient development of the jet for bubble B when the bubble A nearly reaching the end of jet stage. Therefore, when the initial central distance is less, the reduction degree of radius difference at the jet stage is larger as shown in Fig. 13. Keeping the location of bubble A unchanged, the bubble B is arranged at different vertical distance now. Also, the incident frequency, incident amplitude, initial phase and incident angle are retained as well. Deformation histories for the double bubbles with two typical arrangements are shown in Fig. 15.

Since the bubble A is nearer the rigid wall, it will be rejected by the bubble B and the rigid wall at the expanding stage together. Therefore, the sides near the bubble B and rigid wall become flat. On the other hand, the bubble B is less affected by the rigid wall, which leads to the flat only appear at the side near the bubble A. At the contraction stage, the bubble A is attracted by the bubble B and the rigid wall together, and is elongated. Although the phase of acoustic wave acting on bubble A is advanced, since the influence of bubble B and rigid wall and when bubble B reaches the end of jet stage, bubble A is still at the prophase of jet stage in Fig. 15 (a), i.e. the motion of bubble B is advanced. Even in Fig. 15(b), bubble A is attracted by rigid wall and bubble B along z-axis without forming the jet. And bubble B forms the jet since the only attraction from one side. Also, the jet is deflected to the positive x-axis since the incident wave. Fig. 16 shows the radius difference with different locations of bubble B along z-axis. The other parameters are the same with above cases. With the increment of the vertical distance between bubble B and rigid wall, the limitation from the rigid wall is weakened. Thus, the oscillation amplitude of bubble B will be amplified, which leads to the increment of the radius difference at the expanding stage. At the contraction stage, since the offset of bubble B along positive z-axis, the time forming the jet is ahead. Therefore the time for radius difference reaching the maximum value is forward as well.

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

25

513

160 1.8 1.5 1.2

140 Jet tip velocity/ms-1

Radius/μm

20

15

10

1.8 1.5 1.2

0

120 100 80 60 40 20

1

2

3

0

0

0.5

1

1.5

t/μs

2 t/μs

2.5

3

3.5

4

Fig. 10. Characteristics of cavitation bubble with different incident amplitudes, the legend in figure means the radio to the base amplitude Pa ¼ 1 atm. (a)Radius, (b)Jet tip velocity.

2 2.5 3 3.5

ΔR/μm

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

t/μs Fig. 11. History of the radius difference with different incident angles (f¼ 224.48 KHz, Pa ¼ 1.8 atm, θ0 ¼ π, non-dimensional coordinate of bubble A is (0,0,1.5), and non-dimensional coordinate of bubble B is (3,0,1.5)).

Fig. 13. History of the radius difference with different initial central distances (f ¼224.48 KHz, Pa ¼1.8 atm, θ0 ¼π, θi ¼ π/4, the non-dimensional initial central distances are 2.5,3,3.5, respectively).

2.5

3

3.5

90 45 0 Fig. 12. Shape comparison at the jet stage with different incident angles (f¼ 224.48 KHz, Pa ¼ 1.8 atm, θ0 ¼ π, non-dimensional coordinate of bubble A is (0,0,1.5),and non-dimensional coordinate of bubble B is (3,0,1.5)).

As the large zB, the bubble A is attracted by bubble B and the rigid wall, the reduction of volume is not obvious compared to bubble B. Therefore, the volume of bubble B is less than that of bubble A gradually in time, which leads the radius difference to be negative. With the increment of zB, the jet of bubble B will be developed more sufficiently since the less limitation from the wall, and the radius at the end of jet stage is decreased. Therefore, the  ΔRj is increased as well.

Fig. 14. Shape comparison at the jet stage with different initial central distances (f ¼224.48 KHz, Pa ¼1.8 atm, θ0 ¼π, θi ¼ π/4, the non-dimensional initial central distances are 2.5,3,3.5, respectively).

3.2.5. Jet direction We fix the location of bubble A, and move the bubble B from near to far to get the different initial central distance Dc. The incident frequency, incident amplitude, initial phase and incident angle keep unchanged. We define the jet direction as the direction of jet tip velocity. And the angle between jet direction and positive x-axis is defined as α. The variation of jet direction with initial central distance is shown in Fig. 17.

514

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

Fig. 15. Deformation history of the double bubbles with different arrangement (f ¼224.48 KHz, Pa ¼1.8 atm, θ0 ¼ π, θi ¼ π/4, non-dimensional coordinate of bubble A is (0,0,1.5)), (a)non-dimensional coordinate of bubble B is (3,0,2.5), t¼ 0, 1.94, 2.89, 3.09, 3.22 μs, (b)non-dimensional coordinate of bubble B is (0,0,4.5), t¼ 0, 1.65, 2.79, 3.12 μs.

1 0

zB

ΔR

-1

Bubble 2

Bubble 9

Bubble 8

Bubble 3

Bubble 1

Bubble 7

Bubble 4

Bubble 5

Bubble 6

-2 -3 -4

Cavitation B at (0,0,4.5)

-5 -6

0

0.5

1

1.5

2

2.5

3

3.5

t/μs Fig. 16. Radius difference with different location of bubble B(the arrow in figure means the increment direction of zB, and zB ¼ 1.5,1.6,1.7,1.8,1.9,2,2.1,2.2,2.3,2.4,2.5, the dash line represents the radius difference corresponding to Fig. 5, double dot dash line represents the radius difference corresponding to Fig. 15(a), and dash dot line represents the radius difference corresponding to Fig. 15(b)).

x

Rigid wall

Fig. 18. Schema for the cavitation bubble array (9 cavitation bubbles).

the main factors affect the jet direction will become the wall and incident acoustic wave. When the bubble B is move to infinity, i.e. there is only one bubble in the acoustic wave filed near rigid wall, the jet direction is just affected by the wall and the acoustic wave, and the incline level of jet is lowest now,as shown in Fig. 17 when Dc ¼ 1. As the cases shown in Fig. 12 about the bubbles motion with different incident angle, we get tanα as 0.98, 0.85 and 0.83 when incident angle is 0°, 45° and 90° respectively. Thus, the incident angle has a little influence on the jet direction, and the greater incident angle, the higher incline level of jet.

5 4 tanα

y

3 2 1 3

4

5

6

7 Dc

8

9

10

∞

Fig. 17. Jet direction with different initial central distance.

We find that with the increment of Dc the angle α is increased, which means the incline of jet is decreased. The reason is the interaction between bubble B and jet of bubble A is weakened, and

3.3. Cavitation bubble array Based on the analysis of two cavitation bubbles system, the research about the dynamics of cavitation bubbles array near the rigid wall subjected to acoustic wave is implemented. Here, 9 bubbles in equilibrium state are released near the rigid wall, and the vertical distances between the bubbles and rigid wall are set as 1.5. The initial radius, incident frequency, incident amplitude and incident angle are set to be the same with that in Section 3.2. The

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

515

Fig. 19. Shape history of the cavitation bubble array(f ¼224.48 KHz, Pa ¼ 1.8 atm,θ0 ¼ π, θi ¼π/4), (a)t ¼0 μs, (b)t¼ 2.14 μs, (c)t¼ 2.80 μs, (d)t ¼3.03 μs, (e)t ¼3.39 μs, (f)t¼ 3.54 μs.

20

R/μm

18

16

14

Bubble 1 Bubble 2 Bubble 3 Bubble 5 Bubble 6 Bubble 7

12

10

0

0.5

1

1.5

2

2.5

3

3.5

t/μs Fig. 20. Radius for each bubble in the array (f¼ 224.48 KHz, Pa ¼1.8 atm, θ0 ¼ π, θi ¼ π/4).

arrangement of 9 bubbles is shown in Fig. 18, and the central distances for the vicinity bubbles are all equaling to 3. Fig. 19 shows the deformation history for the bubbles array, and considering the symmetry, the motion characteristics of bubble 2 and 4, 5 and 9, 6 and 8 are same respectively. Fig. 20 shows the radius of each bubble in the array. Bubbles in array expand at first, but with the rejection between bubbles, the amplitude of expanding is just small. As shown in Fig. 20, we can obtain the time for reaching the maximum radius of each bubble. Since the interaction, the motion of the four corner bubbles is advanced to other bubbles, and the bubble 1 is delayed to the others. As the phase of the drive acoustic pressure acting on the bubbles, the motion of bubble 2(4) is advanced than 6(8), and bubble 3 is advanced than 5(9) and 7. Then, bubbles in the array begin to contract in successively. And the jet will be deflected to the rigid wall. As the analysis in Section 3.2, the motion of bubble 2(4) having small x coordinate is advanced to others, and contracts first, as shown in Fig. 19(c)– (e) and Fig. 20. Fig. 19(f) shows the shape for each bubble when bubble 2(4) reaches the end of jet stage. It indicates that the bubble 6(8) with the largest phase lag is also forming the jet. Since the attraction between bubbles; bubble 3,5,7,9, which are attracted

by the bubbles arranged beside them respectively, are still at the prophase of jet stage, although the phase of acoustic pressure acting on these bubbles is equal to or advanced than bubble 6(8). And the central bubble 1 is not able to form the jet. In additional, since the attraction from bubble 1, all the jets of bubbles will be deflected to the central bubble as well. As discussed above, for the no corner bubbles, the motion is delayed since the interaction. Therefore, as shown in Fig. 20, the time for reaching maximum radius is later and the maximum radius is limited. The sequence for non-corner bubbles as the maximum radius from large to small is 7, 5(9) and 3, which can be suggested according to the phase of sound pressure acting on these bubbles. The bubble 1 is at the center of the array, which is affected obviously by the interaction with the latest time for reaching the maximum radius. However, since the longer time in the expanding stage, the maximum radius of bubble 1 is larger than that of bubble 5(9) with the same acting sound pressure phase. For the corner bubble 6(8), since the less effect from the interaction between the bubbles, and the phase lag of acting sound pressure, which leads to the long expanding time, the maximum radius is larger than other bubbles. And the maximum radius of corner bubble 2(4) is less than that of bubble 6(8) even bubble 7, since the advancement of the acting sound pressure phase leading to the short expanding time.

4. Conclusions Wave equation is utilized to approximate the dynamics of cavitation bubble in acoustic field with dividing the motion into prophase and anaphase. The new form of boundary integral equation considering compressibility is presented. And modified Bernoulli equation as the boundary condition for cavitation bubble in acoustic field is also deduced. Multiple bubbles with different arrangements are released near the rigid wall. An incident traveling wave comes to the rigid wall with arbitrary angle, and form a stable acoustic field at one side of rigid wall. Bubbles will expand or contraction since the driving of acoustic wave, and generate the high speed jet pointing to both rigid wall and around bubbles. Also if the acoustic wave is incident oblique to the rigid wall, the jet will be deflected as well. Otherwise, the interaction between rigid wall,

516

X. Ye et al. / Ocean Engineering 109 (2015) 507–516

bubbles and acoustic wave is discussed with different parameters, such as incident frequency, amplitude, angle and arrangements. Some conclusions are summarized as follow: (1) The interaction between the multiple cavitation bubbles near the rigid wall is mainly rejection in the expanding stage, but attraction in the contraction stage. With the different arrangements, interaction between bubbles and rigid wall is varied since the distance, as well as the phase difference of acting acoustic wave on the bubbles surface. (2) Since the consideration of compressibility, the maximum radius at expansion stage will be reduced and the time for once expanding and contraction is shorted. (3) Low incident frequency will slow the motion state transformation, and lead to sufficient expansion with larger maximum radius, which will cause great limitation from rigid wall and make the radius at jet stage large, as well as the high incident amplitude. Also, there is a severe contraction after the great maximum radius, which leads to large jet tip velocity. (4) With the increment of incident angle, the motion delay between bubbles will also be enlarged, as well as the oblique degree of jet of advanced bubble.

Acknowledgement This work is supported by National Natural Science Foundation of China (51479041, U1430236, 51279038).

References Blake, J.R., Gibson, D.C., 1981. Growth and collapse of a vapour cavity near a free surface. J. Fluid. Mech. 111, 123–140. Crum, L.A., Reynolds, G.T., 1985. Sonoluminescence produced by ‘stable’ cavitation. J. Acoust. Soc. Am. 78, 137–139. Flynn, H.G., 1975a. Cavitation dynamics: I. A mathematical formulation (in sound field). J. Acoust. Soc. Am. 57, 1379–1396. Flynn, H.G., 1975b. Cavitation dynamics: II. Free pulsations and models for cavitation bubbles. J. Acoust. Soc. Am. 58, 1160–1170. Gaitan, D.F., Crum, L.A., 1990. Frontiers of Nonlinear Acoustics In: Hamilton, M.F., Blackstock, D.T. (Eds.), Proceedings of the 12th ISNA. Elsevier Applied Science, London. Gaitan, D.F., Crum, L.A., Church, C.C., Roy, R.A., 1992. Sonoluminescence and bubble dynamics for a single,stable,cavitation bubble. J. Acoust. Soc. Am. 91, 3166–3183. Keller, J.B., Miksis, M.J., 1980. Bubble oscillations of large amplitude. J. Acoust. Soc. Am. 68, 628–633. Lastman, G.J., Wentzell, R.A., 1981. Comparison of five models of spherical bubble response in an inviscid compressible liquid. J. Acoust. Soc. Am. 69, 638.

Lastman, G.J., Wentzell, R.A., 1982. On two equations of radial motion of a spherical gas-filled bubble in a compressible liquid. J. Acoust. Soc. Am. 71, 835. Liu, Y.L., Zhang, A.M., Tian, Z.L., 2014. Approximation of underwater explosion bubble by singularities based on BEM. Ocean Eng. 75, 46–52. Michael, L.C., Lindau, O., Blake, J.R., Szeri, A.J., 2007. Shape stability and violent collapse of microbubbles in acoustic traveling waves. Phys. Fluids 19, 047101. Morse, P.M., Ingard, K.U., 1968. Theoretical Acoustics. McGray-Hill Book Company, United States. Ni, B.Y., 2012. Study on the Motion and Loads Characteristics of Underwater Viscous Bubbles (Cavitation) (Ph.D. thesis). Harbin Engineering University, Harbin. Park, J., 2008. A coupled Runge-Kutta Discontinuous Galerkin-Direct Ghost Fluid (RKDG-DGF) Method to Near-field Early Time Underwater Explosion (UNDEX) Simulations (Ph.D. thesis). Virginia Polytechnic Institute and State University, Virginia. Prosperetti, A., Lezzi, A., 1986a. Bubble dynamics in a compressible liquid. Part I. First-order theory. J. Fluid Mech. 168, 457–478. Prosperetti, A., Lezzi, A., 1986b. Bubble dynamics in a compressible liquid. Part II. Second-order theory. J. Fluid Mech. 185, 289–321. Qiu, J.X., Liu, T.G., Khoo, B.C., 2008. Simulations of compressible two-medium flow by Runge-Kutta discontinuous Galerkin methods with the ghost fluid method. Commun. Comput. Phys. 3, 479–504. Sun, H., 2007. A Boundary Element Method Applied to Strongly Nonlinear WaveBody Interaction Problems (Ph.D. thesis). Norwegian University of Science and Technology, Trondheim. Wang, C., Khoo, B.C., Yeo, K.S., 2003. Elastic mesh technique for 3D BIM simulation with an application to underwater explosion bubbles. Comput. Fluids 32, 1195–1212. Wang, Q.X., 2005. Unstructured MEL modeling of unsteady nonlinear ship waves. J. Comput. Phys. 210, 183–224. Wang, Q.X., Blake, J.R., 2010. Non-spherical bubble dynamics in a compressible liquid Part 1: travelling acoustic wave. J. Fluid Mech. 659, 191–224. Wang, Q.X., Blake, J.R., 2011. Non-spherical bubble dynamics in a compressible liquid Part 2: acoustic standing wave. J. Fluid Mech. 679, 559–581. Yao, X.L., Ye, X., Zhang, A.M., 2013. Cavitation bubble in compressible free fluid subjected to the travelling wave. Acta Phys. Sin. 62, 244701. Ye, X., Zhang, A.M., Yao, X.L., 2014a. The motion characteristics of cavitation bubble near the rigid wall with the driving of vertical-incidence acoustic wave. China Ocean Eng. 29, 17–32. Ye, X., Yao, X.L., Sun, L.Q., Wang, B., 2014b. Cavitation bubble in compressible fluid near the rigid wall subjected to the acoustic wave with arbitrary incidence angle in three-dimensional. J. Mech. 307–318, http://dx.doi.org/10.1017/jmech. 2014.77. Zhang, Y.L., Yeo, K.S., Khoo, B.C., Wang, C., 2001. 3D jet impact and toroidal bubbles. J. Comput. Phys. 166, 336–360. Zhang, A.M., Xiao, W., Wang, S.P., 2013a. Experimental investigation of the interaction between a pulsating bubble and a rigid cylinder. Acta Mech. Sin. 29, 503–512. Zhang, A.M., Wang, S.P., Huang, C., Wang, B., 2013b. Influences of initial and boundary conditions on underwater explosion bubble dynamics. Eur. J. Mech. B/Fluids 42, 69–91. Zhang, A.M., Wang, S.P., Wu, G.X., 2013c. Simulation of bubble motion in a compressible liquid based on the three dimensional wave equation. Eng. Anal. Bound. Elem. 37, 1179–1188. Zhang, A.M., Cui, P., Cui, J., Wang, Q.X., 2015a. Experimental study on bubble dynamics subject to buoyancy. J. Fluid. Mech. 776, 137–160. Zhang, A.M., Li, S., Cui, J., 2015b. Study on splitting of a toroidal bubble near a rigid boundary. Phys. Fluids 27, 062102. Zhang, A.M., Liu, Y.L., 2015. Improved three dimensional bubble dynamics model based on boundary element method. J. Comput. Phys. 294, 208–223.