1 Introduction

Bubble frequently exists in multiphase flow,and results in discontinuousness in the liquid. As one hot topic for fluid mechanics,the motion characteristics of bubble have been widely studied in propeller cavitation^[1],underwater explosion^{[2, 3]},micro bubble anti-drag^[4],etc. With the perturbation induced by growth and contraction of bubble,the mechanical energy of bubble will transfer to the sound energy in compressible fluid. For a ship,there is a large quantity of bubbles on the surface of the propeller^[1],even of the hull^[5],and in the wake flow^[6]. The sound pressure radiation from these bubbles will break the stealth of the ship and interfere the performance of sonar^[7]. Therefore,the research about the motion characteristics and the induced sound pressure is valued for practical engineering.

When the velocity of the bubble wall is comparable to acoustic speed,the compression fluid will have great influence on the motion of bubble^{[8, 9]}. Based on the incompressible model,with neglecting the energy loss since the process of bubble motion^[10],the pulsation amplitude of bubble and the sound pressure at any point of the fluid will not be reduced with time. Since sound propagates to far-filed as the form of acoustic wave,the total energy of fluid including kinetic energy,potential energy,and acoustic energy is not conservational.

For the spherical bubble in compressible fluid,the motion can be solved by the equations in Refs. ^{[8, 9]},where compression is taken into account. For non-spherical bubble,the motion is not able to be obtained with analysis,and thus numerical simulation is adopted. Wang and Blake^{[11, 12]} expanded the spherical symmetric model^{[8, 9]} to the axisymmetric model to simulate the motion of non-spherical bubble induced by acoustic pressure in compressible fluid. The fluid was divided into an incompressible part in near-field and a compressible part in far-field. But the model in Refs. ^{[11, 12]} was too complex,and the compression was not considered for the whole fluid. With the idea from the doubly asymptotic approximation (DAA) method^{[13, 14, 15]}, which matches the plane wave and potential flow in the frequency domain,Zhang^[16] deduced the boundary integral equation in wake compressible fluid with the match of bubble motion in prophase and anaphase. With the method in Ref. ^[16],Wang calculated the motion of explosion bubble in compressible fluid.

For the far-field sound radiation induced by single spherical bubble in free field,as the approximation of acoustical simulation,bubble is simplified to a monopole source,and the sound radiation is just caused by the volume pulsation of the bubble^{[17, 18]}. For non-spherical bubble,the solution procedure for motion and noise radiation is separated in the researches presented. The dynamic characteristics of bubble are solved based on the incompressible nu- merical model,and the boundary integral equation is used to calculate the sound pressure distribution in the fluid. The motion of cavitation in incompressible vortex flow was solved by Choi and Georges^{[19, 20]},and the sound pressure radiated to fluid was obtained by the bound- ary integral equation corresponding to the Laplace equation. But the compression of fluid and the retardance of acoustic propagation were not considered. Lu and Qi^[7] and Qi and Lu^[21] proposed the mixture boundary element method (BEM) which calculated the motion of bubble with the incompressible potential flow theory and the sound radiation pressure with a fixed vir- tual surface. The complex moving boundary problem about the solution of sound pressure was simplified to a fixed boundary problem. For the superposition of the perturbation in Refs. [7, 21],the mesh of the virtual surface depended on the relative locations of the observer and the bubble,and the time step was related to the element size. However,since the re-mesh of virtual surface was taken for different observers which were not convenient for multiple observers,the reduction of the element size would increase the calculation costing. Jamaluddin and Ball^[22] and Turangan and Jamaluddin^[23] utilized the free Lagrange method (FLM) to solve the col- lapse of the bubble induced by shock wave. The far-field radiation pressure was obtained by the Kirchhoff formula^{[24, 25]} and the Ffowcs Williams-Hawkings (FW-H) formula^{[26, 27]} for fixed boundaries. The instability of FLM at high wave numbers could be avoided by adding artificial surface tension on the particle at the interface. For the simulation of free field,the artificial non-reflecting boundary condition was needed at proper locations. Compared with the BEM, the cost of FLM is higher.

In this paper,based on the compressible fluid theory,the three-dimension BEM is used to solve the motion of bubble in free field. The low order evolution of the Kirchhoff retarded potential formula is used to obtain the velocity potential distribution on a fixed radiation surface which arranges at the near-field of the bubble. With the radiation surface treated as the new noise source,the sound radiation pressure in the unbounded-field is calculated by the Kirchhoff and FW-H formulas for fixed boundaries. The match is implemented for prophase and anaphase with the help of Laplace transformation to obtain the boundary integral equation in compressible fluid. On the other hand,with the consideration about the retardance and superimposition of the sound pressure at any point in the fluid,the match between the near- field and far-field at the two sides of the radiation surface is implemented. As above,the radiation noise during the bubble motion process in the fluid is obtained. 2 Basic theory 2.1 Boundary integral equation in compressible fluid

The three-dimension model for the bubble simulation in the free field is shown in Fig. 1. The buoyance direction is set along the z-axis. The origin of the Cartesian coordinate system is located at the center of the initial spherical bubble. The z-axis points the opposite direction to gravity,that is to say,the jet direction is along the z-axis. The fluid around the bubble is assumed as the irrotational and compressible ideal fluid. The motion of the bubble is supposed as an adiabatic procedure,and the surface tension is ignored.

In the whole fluid field,the velocity potential ϕ satisfies the linear wave-motion equation

the boundary condition at the surface of the bubble is given by and the initial condition is given by with the assumption about small amplitude perturbation,where the sound speed c approximates the sound speed c_∞ at the infinity fluid,u is the normal velocity of the bubble surface,and r is the location vector. The Green function for the Eq. (1) is defined as^[28] where G is the Green function,δ is the Dirac Delta function,r_p is the location vector of the field point,r_p is the location vector of the source point,t is the time of the field point,and τ′ = t − r_pq/c is the source time. Solving Eqs. (1) and (4) with Green’s second identity,the Kirchhoff retarded potential formula is given by In Eq. (5),α is the solid angle,nq is the unit outward normal vector at the source point,and v is the velocity vector of the bubble surface. To solve the motion of the bubble,the field point and the source point are all arranged on the surface of the bubble. When time is small, the physical quantities on the bubble surface are just related with the minute spherical area around^[13]. Then,the prophase approximation can be given by where κ is the curvature of the bubble surface. With the movement of the bubble,the physical quantities on the bubble surface are related with the whole surface^[15]. Expanding Eq. (5) with Taylor’s expansion after the Laplace transform,and reserving the first-order of rpq/c,the anaphase approximation is given by where a is the acceleration of the bubble surface. For convenience,Eq. (7) is written as where ℑ,Θ,and ℜ are the coefficient matrices corresponding to Eq. (7). Matching Eqs. (6) and (8) with the Laplace transform,the boundary integral equation in compressible fluid is given by^[16] The detail derivation process can be found in Ref. ^[16]. In Eq. (9),κ is the local curvature matrix of the bubble surface,Π = (Θ⁻¹ℑ − κ)(E −Θ⁻¹ℜΘ⁻¹ℑ)⁻¹,and E is the identity matrix. When the sound speed is set to infinite,Eq. (9) can be simplified to the boundary integral equation in incompressible fluid. The velocity potential on the bubble surface is updated with the Bernoulli equation in compressible fluid as follows: where P_∞ and ρ_∞ are the pressure and the density of the fluid in infinite,respectively. The relationship between the inner pressure and the volume of the bubble is given by where P₀ and V₀ are the initial inner pressure and the volume of the bubble,respectively,P_c is the saturated vapor pressure,and γ is the specific heat of gas whose value is set to be 1.25 in this paper.

To make the equations dimensionless,denote the maximum radius of the bubble in the incompressible fluid R_m,the pressure difference ΔP = P_∞−P_c,the density ρ_∞,and the sound speed c∞. The buoyance parameter describing jet and rising phenomena is given by

and the intensive parameter in the bubble is With Eqs. (9)-(11),the motion of the bubble in the compressible fluid is obtained.

The valid of the boundary integral equation (9) to solve the motion of the bubble is verified as follows. The motion of single spherical bubble in free field is solved with the method presented in this paper as follow:

The result is compared with that of the Rayleigh-Plesset equation^[29]. The modified motion equation in compressible fluid presented by Prospertti and Lezzi^{[8, 9]} is as follows: The dimensionless parameters are R′ ₀ = 0.163,ε = 77.4,and ω² = 0 (R'⁰ is the non-dimensional initial radius).

In Fig. 2,the results calculated by the boundary integral equation in this paper is matched with the results solved by the equation presented by Prospertti and Lezzi in Refs. ^{[8, 9]}. With the consideration for compression,the maximum radius and the inner pressure of the bubble decrease with the cycle count.

For the non-spherical bubble,the numerical simulation in this paper for the bubble motion is stopped at the moment when the jet impacts the other side of the bubble. After reaching the maximum volume,the bubble begins to contract. In this paper,the process from the moment that the maximum volume appears to the moment that the jet touches the other side of the bubble is defined as the contraction stage. Figure 3(a) is the result for the experiment about the explosion bubble with 1.5 g PENT and 1m depth^[30]. Figure 3(b) is the numerical simulation result corresponding to the experiment. The quantities in Fig. 3 are non-dimensional. Thecolored contour in Fig. 3(b) represents the velocity potential distribution on the surface of the bubble. The result of the numerical method agrees well with the experimental data.

2.2 Solution for radiation noise induced by bubble motion

The sound radiation pressure induced by the growth and contraction of the bubble is solved by the Kirchhoff formulation and FW-H formulaes. A fixed and closed radiation surface is arranged to enclose the bubble,and the radiation nodes are set on the radiation surface,as shown in Fig. 4.

Substituting the velocity potential,the normal velocity,and the normal acceleration of the bubble surface solved by Eqs. (9)-(10) into the low-order expanded formula of the retarded potential equation (7) yields the velocity potential distribution on the radiation surface. In Eq. (7),the field point p is set on the radiation nodes,the source q is set on the bubble surface, and the solid angle α = 4π. The radiation surface is treated as a new source,and the sound pressure in the fluid is calculated by the Kirchhoff and FW-H formulas.

The Kirchhoff formula on the fixed radiation surface is given by^{[24, 25]}

In Eq. (16),ϕ is the velocity potential in the unbounded acoustic field out of the radiation surface,and satisfies the linear wave-motion equation (1). n_c is the unit normal vector pointing outward on the radiation surface. r = x − y,where x is the location of the observer and y is the location of the radiation node. r = |x − y| is the distance between the observer and the radiation node. The radiation time τ′ = t − τ,where t is the time at the observer,and τ = r/c_∞ is the propagation time,i.e.,the retarded time. The velocity potential of the observer at time t is superimposed about the perturbation induced by the radiation nodes at τ′ after the propagation with τ.

The velocity potential at any location in the fluid out of the radiation surface can be solved by the Kirchhoff formula (16). Then,the sound pressure is acquired. Equation (16) indicates that the solution of the Kirchhoff formula needs the history of the velocity potential and its derivative respect to bubble’s time. All the parameters mentioned above can be obtained by Eq. (7). Since the Kirchhoff formula requires perfect satisfaction for the linear wave equation, the radius of the radiation surface would be large enough to enclose all nonlinear factors in the fluid. Moreover,to gain the precise velocity potential on the radiation surface,Eq. (7) requires the radiation surface to be set near the bubble surface since with the incensement of τ =r_pq/c_∞, neglecting high orders of τ will cause an obvious error. Meanwhile,when the retarded time of propagation increases,the velocity potential on the radiation surface at some time is difficultly calculated with Eq. (7). The time interpolation and superimposition for the velocity potential on the radiation surface is required,i.e.,the solving procedure becomes complex. As mentioned above,when using the Kirchhoff formula to solve the noise field induced by the bubble,the radiation surface is expected to arrange at a proper location.

The arrangement of radiation surface with the FW-H formula is flexible compared with the Kirchhoff formula. Since the radiation surface does not require to enclose the nonlinear factor in the fluid,the radiation surface can just be arranged at the near-field of the bubble^[22]. The FW-H formula is given by^[26]

where r,r,τ′,x,and y are the same as those in the Kirchhoff formula. ϑ = ρvc ·nc,ρc and vc are the fluid density and the velocity at radiation surface,respectively. ζ = (P_i · n_c + ρv_cu_c)·r, where Pi and uc are the pressure tensor and the normal velocity at the radiation surface.

The FW-H formula can acquire the sound pressure directly at any location in the fluid out of the radiation surface. The solution for Eq. (17) requires the history of the pressure tensor, density,velocity,and normal velocity at the radiation surface. The physical quantities at the radiation surface mentioned above can be obtained with Eq. (7) and the difference of the velocity potential in time and space.

The retarded time of the observer and radiation nodes is defined as the propagation time of acoustic wave simply. To obtain the history of the sound pressure at observer,the velocity potential solved by Eq. (7) at the radiation nodes is interpolated for time. After that,the velocity potential or sound pressure at observer is calculated with Eq. (16) or Eq. (17),and the contribution to perturbation at observer induced by each radiation node is superimposed at observer’s time t. Since the requirement for the known quantities of the Kirchhoff formula on the radiation surface is lower than that of the FW-H formula,the Kirchhoff formula is easier and simpler to be implemented. But as mentioned above,it should be paid attention to that the Kirchhoff formula is sensitive to the arrangement of the radiation surface by which the size of the radiation surface should be appropriate while the FW-H formula is not sensitive to the arrangement of the radiation surface. 3 Numerical results 3.1 Verification for acoustical pressure calculation

With the BEM in compressible fluid and Kirchhoff or FW-H formulas,the motion of the bubble and its radiation noise is obtained. To verify the method in this paper,the far-field radiation induced by single bubble is solved analytically at first. Then,the result is compared with that calculated by the method in this paper. The non-dimensional initial parameters for spherical bubble are R'⁰ = 0.35 and ε = 10.24. The non-dimensional observer location is x′ = 0, y′ = 0,and z′ = 15. The analytical formula for sound radiation pressure is shown as follows^[17]:

where R is the bubble radius obtained from Eq. (15),rm is the distance between observer and the center of the bubble,and (r_m − R)/c_∞ represents the retarded time.

The non-dimensional results for comparison are shown in Fig. 5. The sound pressure is non- dimensional with the reference pressure ΔP = P_∞−P_c mentioned above. The far-field acoustic radiation calculated by the Kirchhoff and FW-H formulas agrees well with the analytical result. With different radii of the radiation surface,the two formulas can all obtain the correct solution. Since the fluid around the bubble satisfies the linear wave equation (1) at any point with the method presented in this paper,the nonlinear effect is neglected in the fluid. When the radius of the radiation surface is 1.1 times of bubble’s maximum radius,i.e.,the radiation surface is very close to the bubble surface,the result calculated by the Kirchhoff formula agrees well with the analytical result. However,the FW-H formula needs more known quantities than the Kirchhoff formula,which needs more time and space in the computation. Therefore,the numerical error of the FW-H formula is higher than that of the Kirchhoff formula.

3.2 Noise characteristics of non-spherical bubble

For the convenience about implementation and the numerical precision,the sound radiation induced by non-spherical bubble is solved by the Kirchhoff formula. The non-dimensional initial parameters of the bubble are R'⁰ = 0.13,ε = 142.63,and ω² = 0.1. The calculation is stopped when the jet touches the other side of the bubble.

The sound radiation pressure induced by the spherical bubble and the non-spherical bubble with the same initial radii and intensive parameters are shown in Fig. 6. The non-dimensional location of observer is x′ = 0,y′ = 0,and z′ = −3. At the initial phase,since the jet has not been generated and the shape of the bubble is near the sphere,the sound pressure induced by the non-spherical bubble is similar to that of the spherical bubble. In the contraction phase, after the bubble reaches the maximum radius,the bubble rises obviously and the jet at the bottom is formed.

In Fig. 7,the volume variation speeds of the spherical bubble and the non-spherical bubble are given. At the early stage of contraction,the acceleration of volume’s variation for the non- spherical bubble is similar to that for the spherical bubble,while at the anaphase of contraction, after the forming of jet,the acceleration of volume’s variation for the non-spherical bubble is lower than that of the spherical bubble. For far-field radiation,the contribution to the sound pressure from acceleration of volume’s variation is more obvious^[17]. Otherwise,the time history of the equivalent center on the z-axis of the non-spherical bubble is shown in Fig. 8,considering the buoyance,the bubble rises with time,and the speed of rising is quicker after the forming of jet. That is to say compared with the spherical bubble without buoyance,the non-spherical one is farer from the bottom observer. As mentioned above,the far-field amplitude of the sound pressure for the non-spherical bubble is less than that for the spherical bubble when the observer is at the bottom of the bubble after the forming of jet.

In Fig. 9,the comparison curves for the sound radiation pressure at different observers are given. The non-dimensional locations of the observers are shown in Table 1. Figure 9(a) presents the sound radiation pressure at the observers which have different distances from the bubble’s initial center (0,0,0). The retardance from the propagation of acoustic wave is obvious. When the distance increases,with the straggling of the sound energy,the amplitude of the sound pressure at the observer is reduced apparently for t > 0.2. Figure 9(b) presents the sound radiation pressure at the observers around the bubble with the same distances from the initial center.

The contours for various speeds of the velocity potential at the radiation surface are shown in Fig. 10. The red area in dotted box represents higher variation speeds for the velocity potential,while the blue area represents lower ones. When 0.2 < t < 1.0,with the growth and increase in the bubble,the red areas in the contours increase as well. The bubble reaches the maximum volume near t = 1,and begins to contract. The maximum variation speed of the velocity potential gets the top of the radiation surface. At the anaphase of contraction,the bubble begins to generate the jet at the bottom,and the red area also moves to the bottom of the bubble. With the motion of the jet and the increase in the bubble,the red area rises as well.

To analyze Fig. 9(b) associated with Fig. 10,before contraction,the bubble still approximates sphere,and its rising displacement is small. Therefore,the sound radiation pressured at different observers are similar. In pace with the contraction of the bubble,when t > 1,since the jet forms at the anaphase,the rising speed is higher,and the difference of the sound pressure between the observers becomes obvious. At the early stage of contraction,the observer D near the topof the bubble has the highest amplitude of sound pressure,and the amplitude of observer E at the side of bubble is just less than that of observer D,while the amplitude of observer A near the bottom of the bubble is the lowest. At the anaphase of contraction,with the generation and moving up of the high speed jet,which also cause a higher rising speed of the bubble,the red area in Fig. 10 will lift,while the observer which has the highest amplitude of sound pressure is firstly the bottom observer A,then moves to the side observer E,and finally reaches the top observer D. That is to say,the existence of jet will affect the sound pressure around the bubble,and the point with the highest amplitude of sound pressure will move with the motion of the jet tip. 3.3 Influence of initial parameters on radiation characteristics of bubble

The curves in Fig. 11 represent the sound radiation pressure under different initial parame- ters. The detail setting parameters are shown in Table 2.

For convenience,we discuss the relation between R'⁰ and ε in the incompressible fluid. Substituting Eq. (11) into Eq. (14),we can get the expression as follows:

When the intensive parameter or the initial radius is changed,the initial inner pressure of the bubble will be different as well,which will affect the motion characteristics of the bubble,even the radiation sound pressure induced by the bubble.

The comparison curves for sound pressure with different buoyance parameters at observer A are shown in Fig. 11(a). At the stage of growth,since the bubble still approximates sphere,the influence of the buoyance parameter on the acoustic radiation is low. As shown in Fig. 12(a), at the early stage of contraction,since the jet tip velocity of the non-spherical bubble increases with the buoyance parameter,but not significant,the perturbation increases a little. Therefore, the induced sound pressure is a little greater with higher buoyance parameters at this stage, and will be greater obviously with the process of contraction. For the anaphase of contraction with the apparent jet,the tip velocity reduces with the increment of the buoyance parameter, and the perturbation to the fluid around also decreases. Thus,higher buoyance parameters will result in lower amplitude of sound pressure.

The comparison curves in Fig. 11(b) represent the sound radiation pressure at observer A under different intensive parameters (initial radii). Since larger intensive parameters lead to larger inner pressure and faster expansion speed for the growth of the bubble in the initial stage, more severely volume pulsation and higher amplitude of sound pressure will be generated. As shown in Fig. 12(b),at the stage of contraction,the jet tip speed is higher with larger intensive parameters,so is the perturbation to fluid around. Therefore,larger intensive parameters will induce higher amplitude of sound pressure at the contraction stage.

4 Conclusions

Based on the compressible fluid theory,the three-dimension BEM is utilized to solve the motion of the bubble. Associating with the arrangement of the radiation surface,the Kirchhoff formula and the FW-H formula are adopted to calculate the sound radiation pressure in the growth and contraction phase of the bubble. The comparison for the noise radiation of non- spherical bubble and spherical bubble at different location in compressible fluid is presented. The influences of different initial radii,intensive parameters,and buoyance parameters on radiation field for non-spherical bubble are analyzed. From the results in this paper,we can present some useful conclusions.

Since the acceleration of non-spherical bubble’s volume variation is less than that of spherical bubble at the contraction stage,the far-field sound pressure induced by the former is lower than that by the later. With the increment of the distance between observer and the initial center of the bubble,the sound pressure in the fluid is reduced,and the retardance effect for the propagation of sound wave is obvious.

In the phase of growth,the sound pressure at observers which have the same distance from the initial center of the bubble in the fluid is similar. Since the rising of the bubble,the maximum amplitude of the sound pressure in the fluid appears near the top of the bubble at the early stage of contraction. At the anaphase of contraction,with the obvious jet,the maximum amplitude of sound pressure moves from the location below the bottom of the bubble to the location above the top of the bubble.

For the buoyance parameter,in the stage of growth,the influence of the buoyance parameter on noise radiation is small; in the early stage of contraction,with the enlargement of the buoyance parameter,the amplitude of sound pressure increases a little; in the anaphase of contraction,with obvious jet,smaller buoyance parameters will generate higher amplitude. Higher intensive parameter (less than the initial radius) will induce greater amplitude of sound pressure in the whole procedure of the bubble motion because of the greater initial inner pressure and more severe oscillation of the bubble’s volume.