Numerical study of spike characteristics due to the motions of a non-spherical rebounding bubble

52

2016,28(1):52-65

DOI: 10.1016/S1001-6058(16)60607-8

Numerical study of spike characteristics due to the motions of a non-spherical rebounding bubble* Jia-xia WANG (王加夏)1, Zhi ZONG (宗智)1, Lei SUN (孙雷)1, Zhang-rui LI (李章锐)2, Ming-zuo JIANG (姜明佐)3 1. Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116023, China, E-mail: [email protected] 2. Marine Design and Research Institute of China, Shanghai 200011, China 3. Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian 116023, China (Received February 18, 2014, Revised September 10, 2015) Abstract: The boundary integral method (BIM) is used to simulate the 3-D gas bubble, generated within the two bubble pulsation periods in proximity to a free surface in an inviscid, incompressible and irrotational flow. The present method is well validated by comparing the calculated shapes of the bubble and the free surface with both the experimental results and the numerical ones obtained by the Axisymmetric BIM code. The expansion, the collapse of the gas bubble and the further evolution of the rebounding non-spherical bubble are simulated. The various variation patterns of the free surface spike and the bubble centroid for different standoff distances, the buoyancy parameters and the strength parameters are obtained to reveal the nonlinear interaction between the bubble and the free surface. The amplitude of the second maximum bubble volume and the four typical patterns of the bubble jet and the free surface spike are examined in the context of the standoff distance. The large buoyancy is used to elevate the spray dome rather than the free surface spike. Key words: boundary integral method, toroidal bubble, free surface spike

Introduction The motion of the bubble near the free surface is an important research topic in science, engineering and military fields. Its studies include experiments, theoretical investigations and numerical simulations. Owing to the large nonlinear deformation of the free surface, the interaction of a bubble with a free surface is a difficult problem. During the collapse, a jet pene0F

* Project supported by the National Natural Science Foundation of China (Grant Nos. 51221961, 51279030) the National Key Basic Research Development Program of China (973 Program, Grant Nos. 2013CB036101, 2010CB832704) and the Fundamental Research Funds for the Central Universities (Grant No. L2012016). Biography: Jia-xia WANG (1988-), Male, Ph. D. Candidate Corresponding author: Zhi Zong, E-mail:[email protected]

trates the bubble, which also increases the difficulty of the research especially when the spherical bubble transforms to a toroidal one. In the last three decades, the development of numerical methods helps to solve the difficulties. Some important researches in this area were made by Blake, Wang, Klaseboer et al.. Blake et al.[1] used an approximate integral equation and the boundary integral equation approach to simulate the expansion and the collapse of a vapor bubble with a free surface. Wang et al.[2] investigated the nonlinear interaction of a buoyant vapor bubble close to a free surface numerically, to reveal three collapse patterns in accordance with the buoyancy force. Klaseboer et al.[3] used the BIM to describe the behavior of an oscillating bubble near a resilient/rigid and floating structure. By comparison with experimental results, the BIM proves to be a powerful and accurate numerical tool and is extensively used to simulate the interaction between a bubble and a free surface. Wilkerson[4] applied it on the simulation of the 3-D bubble dynami-

53

cs in 1990. Now the 3-D BIM is widely used in the simulation of bubble dynamics. Zhang et al.[5], Pearson et al.[6], and Li et al.[7] simulated two bubbles beneath the free surface. Klaseboer et al.[8] took into account the viscous effects on the bubble evolution. Wang et al.[9,10] simulated the non-spherical bubble dynamics in a compressible liquid. The subsequent studies include Pearson et al.[11], Zhang et al.[12], Wang et al.[13], Klaseboer et al.[14], Li et al.[15]. In these studies, the BIM is used to simulate the expansion and the contraction of the bubble successfully. The spark discharge generated bubble, observed in experiments[1,16], has one pulsation before disappearing. However, for large-scale underwater explosion bubbles, the existence of non-condensable gases ensures that the bubble will rebound[2,9,12]. At the beginning of the rebounding phase, the high pressure in the bubble was believed to cause the hogging phenomenon[17,18]. The rebounding of the toroidal bubble was identified by using the Axisymmetric BIM in Ref.[3], which is worth further studying. Nevertheless, the 3-D BIM was rarely employed to study the non-spherical rebounding bubble and the continuous free surface spike, because the code tends to break down in the midst of calculation[1]. Thus it is a challenging task to simulate the interaction between a bubble and a free surface within the two pulsation periods using the 3-D BIM. In this work, we mainly focus on the interaction between a gas bubble and a free surface. Due to the strong nonlinear interaction between a bubble and a free surface, the non-uniform free surface mesh, the smoothing scheme[12] and the elastic mesh technique (EMT)[13] method are used to tackle the problems of large free surface deformation, numerical instability and mesh distortions, respectively. Moreover, the numerical surgical cut is used to realize the transformation from a singly connected bubble to a doubly connected one[12]. With the help of these techniques, the BIM codes could be used without breaking down and could run well during the calculation of the highly nonlinear problem. In the present paper, the characteristics of the free surface spike and the movement of the bubble centeriod are revealed, and the effects of the standoff distance, the buoyancy parameter and the strength parameter are identified. 1. Fundamental formulations The fundamental assumption is that the fluid in the studied domain is incompressible and inviscid. Thus, a potential theory is adopted here[1-3]. The surface tension effect is neglected according to Blake[1], where it was pointed out that the surface tension is insignificant over much of the bubble’s lifetime for large bubbles (of size of millimeter or larger). Let Φ ( x, y, z ) be the potential function in the fluid domain Ω with the boundaries of the bubble surface Sb

and the free surface S f . A Cartesian coordinate is used. According to the Green’s second theorem, the governing equation can be expressed as 

∫ φ (Q) s

∂G ( P, Q) ∂φ (Q)  G ( P, Q)  d s + c ( P)φ ( P) = 0 − ∂n ∂n  (1)

where P = ( x, y, z ) is the field point, Q = ( x′, y ′, z ′) is the source point, c( P) is the solid angle at the field point P , G ( P, Q) = 1/ Q − P = 1/ rPQ is the free space green function, n is the unit normal vector pointing outward, s is the sum of the bubble surface and the free surface, rPQ is the distance between point P and point Q , and γ is the dimensionless vertical distance of the initial location of the bubble centroid from the free surface. The coordinate system and the geometry used to simulate the interaction between the bubble and the free surface are shown in Fig.1.

Fig.1 The coordinate system and geometry used to simulate the interaction between the bubble and the free surface

The kinetic boundary condition is Dr = ∇φ Dt

(2)

where r = ( x, y, z ) is the position vector, and D / Dt is the material time derivative. The dynamic boundary conditions of the bubble and the free surface are expressed by the nonlinear Bernoulli equation Dφ P∞ 1 P 2 = + ∇φ − − g z Dt ρ 2 ρ

(3)

where P∞ = Patm + ρ gH is the pressure at the initial location of the bubble centroid, Patm is the atmospheric pressure, ρ is the density of the fluid, H is the depth of the TNT explosives, and g is the gravitational acceleration. For the free surface, the pressure

54

P = Patm , as for the bubble, the pressure P contains

two parts: the constant vapor pressure Pc and that of the non-condensing gas that conforms to the adiabatic law. Owing to the continuity condition of the interface pressure, P = Pc + Pg 0 (V0 / V )2 , where Pg 0 and V0 are the initial pressure and volume of the bubble, respectively, and V is the transient volume of the bubble. The parameter λ is the ratio of the specific heat and is equal to 1.4 for the spark-generated bubble[15], and 1.25 for that generated by the underwater explosion with explosive TNT[3]. For underwater explosion bubbles, the initial conditions can be given as: W  Pg 0 = 7.8 × 108−3λ    V0 

0.543 m and 72.22 MPa, respectively, and the corresponding dimensionless parameters are γ = 1.73 , δ = 0.533, ε = 362.4 . These parameters will be used in the subsequent cases. 2. Numerical scheme The boundaries of the calculated regions are discretized into triangular elements. The triangular mapping of elements in the physical space onto the parametric space is shown in Fig.2.

λ

(4)

1/ 3

 W  Rm = 3.38 ×    H +10 

(5)

where W is the weight of the charge, V0 = 4 π R03 / 3 , and R0 is the initial bubble radius and can be given as[3] 7.8 × 108−3λ 1  3W    λ − 1  4 π Rm3  ∆P

λ

  R 3 − 3 λ   R 3 1 −  0   =  0  − 1   Rm    Rm  (6)

All variables are non-dimensionalized with respect to the length scale Rm (the maximum bubble radius), the pressure scale ∆P = P∞ − Pc , the time scale Rm ρ / ∆P , and the potential scale Rm ∆P / ρ . The non-dimensionalized forms of the unsteady Bernoulli equations are: Dφ 1 2 = ∇φ − δ 2 ( z − γ ) Dt 2 Dφ 1 2 V  = 1 + ∇φ − δ 2 z − ε  0  Dt 2 V 

(7)

Fig.2 The element mapping for the bubble surface

The linear interpolation functions associated with each node for the position and the potential are given with the coordinates (ξ ,η ) as: 3

rm = ∑ N i rm, i 3

φm = ∑ N iφm, i

(8)

Equations (7) and (8) are the Bernoulli equations on the free surface and the bubble, respectively. In the two equations, δ = ρ gRm / ∆P is the buoyancy parameter, ε = Pg 0 / ∆P is the strength parameter, and

γ = H / Rm is the depth parameter. For instance, a weight of 100 kg TNT explosive is initiated at a depth of 10 m, the maximum bubble radius Rm = 5.78 m , the initial bubble radius R0 and pressure P0 are

(10)

i =1

where N1 = ξ , N 2 = η , N1 = 1 − ξ − η , the subscript “ m ” means the mth element and the subscript “ i ” denotes the ith node in Fig.2. The global coordinates of the element should be transformed into the local coordinates to calculate the surface of the triangular element 1 1−η

λ

(9)

i =1

∫ dS = ∫ ∫

0 0

S

J dξ dη

(11)

where J is the Jacobian of the transformation from the global Cartesian coordinates to the local coordinates (ξ ,η ) on the triangle. So Eq.(1) is transformed to M

1 1−η

∑∫ ∫ m =1

0 0

M

∂  1  ∂n  rPQ

1 1−η

∑∫ ∫ m =1

0 0

1 rPQ

 3  ∑ N iφm , i (Q) J m dξ dη + c ( P)φ ( P) =  i =1 3

∑N i =1

i

∂ [φm, i (Q)] J m dξ dη ∂n

(12)

55

where M is the total number of the discrete free surfaces, and rPQ = [( X Q − X P ) 2 + (YQ − YP ) 2 + ( Z Q − Z P )]1/ 2 is the distance between point P and point Q . Equation (12) can be simply written as M

3

M

3

c ( P)φ ( P ) + ∑∑ H m, iφm, i (Q) = ∑∑ Gm, i m =1 i =1

m =1 i =1

∂ φm, i (Q) ∂n (13)

Here, we define: H m, i = ∫

1 1−η

∫

0 0

Gm, i = ∫

1 1−η

∫

0 0

∂  1  ∂n  rPQ

  N i J m dξ dη 

1 J m N i dξ dη rPQ

(14)

Fig.3 The non-uniform triangulation mesh of the free surface

(15)

Equation (13) should be rearranged with the connected nodes by elements. And a matrix equation relating the velocity potentials and the normal velocities of all nodes is formed through two influence matrices H and G , thus  ∂φ  [ H ]N × N [φ ]N ×1 = [G ]N × N    ∂n  N ×1

(16)

where N is the total number of the nodes on the bubble and the free surface. Singularity appears on the diagonal terms of H and G . Polar coordinates transformed from Cartesian coordinates can be used to solve the problem of singularity of the diagonal elements of matrix G [4]. The diagonal elements of matrix H are the solid angles, and the 4π rule[14] is adopted to calculate the diagonal terms for a closed surface. It means that the diagonal terms of the influence matrix H can be obtained by subtracting the other off-diagonal elements of a row from 4π as follows H i, i = 4π − ∑ H i, j

(17)

j ≠i

Another problem is the calculation of the solid angle of a free surface. The free surface model used here is an open surface, thus the 4π rule method cannot be used to calculate the influence matrix H . So a direct approach is adopted based on the definition of the solid angle, which means the area subtended by the linear elements on a unit sphere. More details can be found in Li et al.[7,15]. In order to improve the computational speed and maintain the accuracy, a non-uniform mesh on the free

Fig.4 Comparison between experimental results (in solid line), axisymmetric BIM simulations (in dashed line) and numerical ones (in the deformed mesh) obtained using the present 3-D BIM

surface, as shown in Fig.3, is used with a densely distributed mesh in the region above the bubble and a coarse mesh in the far field. A uniform mesh is first formed and then the following two steps are taken for the transformation to the non-uniform mesh. Firstly, a zoom is used for all the meshes according to α

α

d d x = x  , y = y  , 0 ≤ d ≤ D D D

(18)

Secondly, the central meshes are compacted according to:

56

Fig.5 Comparison between axisymmetric BIM simulations (in thick line) and present 3-D BIM simulation results (in the deformed mesh) β

  d β   d   d β x =  x −  1  x   +  1  x ,   d    d 2   d 

∆ t2 =

β

  d   d1  β y   +   y ,   d 2   d 

0 < d1 ≤ d ≤ d 2 < D

∆φ 1 2 V  max 1 + ∇φ − δ 2 z − ε  0  2 V 

(19)

0 < d1 ≤ d ≤ d 2 < D β   d1  y = y −  d  

∆ t1 =

(20)

where α , β are the transform parameters with the 2 2 values of 0.3 and 0.2, respectively, D = xmax + ymax

( xmax and ymax are the maximum abscissa and ordinate values of the free surface boundary), d1 and d 2 are the upper and lower limit values of the dense mesh region, respectively. A singly connected bubble is transformed to a doubly connected one when the liquid jet impacts on the opposite face of the bubble at the end of the collapse phase. An idealized impact model is adopted where the initial contact between the fluid jet with the bubble wall occurs at just one point to form a circulation in the fluid domain. Zhang[12] extended the axisymmetric model of the jet impact to a 3-D vortex ring model, where the ring is not necessarily planar. Thus, the total potential is decomposed into two parts: one is a potential associated with the circulation generated by the impact, and the other is a remainder. In order to tackle the numerical instability and the severe mesh distortion, a smoothing scheme[12] and the EMT technique[13] should be employed. The non-uniform time step ∆t is chosen for each of the iteration procedure:

∆φ max

1 2 ∇φ − δ 2 ( z − γ ) 2

λ

(21)

(22)

(23)

∆ t = min{∆ t1 , ∆ t2 }

where ∆Φ is a constant, and in this work it is equal to 0.02. 3. Verification and validation To validate our numerical model, we compare the 3-D BIM results with the experiment results in the existing literature and the results obtained by the axisymmetric BIM code. The first comparison is made against the commonly used experimental results obtained at γ = 0.98 by Blake et al.[1] and the results obtained by the axisymmetric code given by Wang[2] during the expansion and contraction phases. For vapor bubbles, the initial value of the potential on the bubble can be obtained using the formulas given by Blake[1]:

φ0 = − R0

t0 = 3

3  2 ∆ P  Rm    − 1 3 ρ  R0   

3ρ 5 Ba  , 2∆P  6

3  2

(24)

(25)

where a = ( R0 / Rm )3 and Ba is the incomplete beta

57

function. The comparison between the experiment results and the 3-D BIM results are shown in Fig.4. The thick line in Fig.4 represents the experiment results, while the mesh gives us the profile obtained by our code. As the non-uniform time step is used here, it is impossible to match the time with the experiment exactly. The solid lines above the free surface in every figure are used to compare the height of the free surface clearly. Generally, the 3-D results are in good agreement with the axisymmetric BIM results and the experiment results, except that the axisymmetric results have a lower free surface elevation which we mainly focus on. It means that during the expansion and contraction phases, the present model is more in accordance with the experiment. However, the 3-D code gives a little faster contraction than the experiments in Fig.4(c), which may be due to the potential theory and the match of the initial conditions. The second comparison is made against the results obtained by the axisymmetric code given by Klaseboer[3] with parameters γ = 1.1 , ε = 500 and

buoyancy parameter of δ = 0.533 , and a dimensionless initial radius of R0 = 0.0939 . The shapes of the bubble and the free surface at some representative times are shown in Fig.6 and some corresponding pressure contours are plotted in Fig.7, and Z means the vertical position of the bubble and the free surface.

δ 2 = 0.25 . This comparison focuses on the rebounding phase of the toroidal bubble. The thick red lines in Fig.5 represent the results obtained by the axisymmetric BIM code, while the mesh is the profile obtained by the present BIM code. A very good agreement is observed between the results obtained by the axisymmetrical BIM code and the 3D BIM results during the expansion and contraction phases, as shown in Figs.5(a) and 5(b). Figs.5(c)-5(f) show the rebounding of the toroidal bubble. For the evolution of the shapes of the bubble and the free surface, on the whole, these two sets of results match well in the general shape, except that our calculation gives a slower and smaller growth than the axisymmetrical BIM results. This discrepancy may result from the smoothing scheme and the EMT technique used in the 3-D code. 4. Numerical results Five cases are shown in this part. The first two cases represent two typical situations for different standoff distances. The evolution of the shapes of the bubble and the free surface and the spike characteristics are given and discussed in details. In the last three cases, the effects of the standoff distance, the buoyancy parameter and the strength parameter are discussed, respectively. All variables in the subsequent sections are in dimensionless form. 4.1 The representative tendency 1: For a large standoff distance In this section, the bubble is initiated at a depth of γ = 1.73 , with a strength parameter of ε = 362.4 , a

Fig.6 For γ = 1.73 , the evolution of the shapes of the bubble and free surface at representative times

The bubble remains remarkably spherical throughout the expansion phase, as shown in Figs.6(a)-6(b) and there is only a small hump on the free surface. At the time t = 0.066 in Fig.7(a), the high pressure promotes the bubble expansion until the bubble reaches the maximum volume, and the pressure around the bubble is uniformly distributed with the pressure below the bubble being higher than that above it. At the time t = 1.700 , a broad upward directed jet is seen at the bottom of the bubble, which is induced by the buoyancy force, larger than the Bjerknes force. The intermediate pressure region is evident below the bubble as

58

minimum volume at the time t = 1.767 with a high pressure region above the bubble, as shown in Fig.7(d). The high pressure will elevate the free surface in the rebounding phase. At the time t = 4.160 , the highest and broadly distributed spray dome is formed on the free surface. The free surface looks like a “plateau”, because the symmetry axis region of the free surface rises more slowly than the region close by. This kind of “plateau” shape of the free surface can also be seen in the experiment of Zhang[16].

Fig.8 Time history of the height of the spray dome, bubble volume and centroid. The time corresponding to capital letters A to D is t = 0.976, 1.767, 2.96 and 4.16, respectively

Fig.7 For γ = 1.73 , pressure contours at representative times

shown in Fig.7(c). Once the jet penetrates the bubble wall, a doubly connected toroidal bubble is formed. Due to the lower pressure inside the bubble, this toroidal bubble goes on collapsing and reaches the second

Figure 8(a) gives the time history of the height of the spray dome, and Fig.8(b) shows the curves of the bubble volume and the bubble centroid. Hs is the height of the spike, and Z 0 is the bubble centriod. The height grows with the expansion of the spherical bubble until the time A, when the spray dome reaches its first peak, and the spherical bubble reaches the maximum volume. Then the volume decreases and the free surface spike falls until the time B, when the spray dome reaches its valley, the bubble reaches the second minimum volume and its centroid moves upward sharply in view of the momentum theorem. Then the bubble begins rebounding. The toroidal bubble reaches the secondary maximum volume at the time C and falls afterward. At the time D, the height of the spray dome reaches the maximum. During the rebounding

59

phase the height of the spray dome continues to grow. This is the first motion tendency of the free surface spike and the bubble centroid.

pressure contours are plotted in Fig.11.

Fig.9 Time history of the displacement of the upper and the lower bubble points and the Kelvin impulse

Figure 9(a) gives the time history of the upper and lower points. The intersection point of the upper and lower poles means the forming of the toroidal bubble. When the toroidal bubble emerges, there are no upper and lower points anymore. So this just appears in the first expansion and contraction phases. During the expansion phase, the upper and lower points are symmetrical. In the contraction phase, the lower point moves upward rapidly under the buoyancy force, and the upper point moves downward slowly under the combined action of the buoyancy and Bjerknes forces. The vertical coordinate of the two line’s intersection is 0.78, and this value will be used in Part 4.2. We note, from Fig.9(b), that the sign of the Kelvin impulse is positive, which means that the bubble migration and the jet direction are directed to the free surface. I is the Kelvin impulse of bubble. 4.2 The representative tendency 2: For a small standoff distance In this case, the evolution of the bubble is initialized at γ = 0.8 below the free surface. The other parameters are the same as in Part 4.1. The profiles of the bubble below the free surface at some representative times are shown in Fig.10 and some corresponding

Fig.10 For γ = 0.8 , the evolution of the shapes of the bubble and the free surface at representative times

By comparing Fig.10(a) with Fig.6(a), it is found that during the early expansion phase there is little difference between these two figures. At the time t = 0.707, the bubble is elongated along the axis of symmetry during the later stage of the expansion phase. When the bubble reaches the maximum volume, the top bubble is entrained into the base of the raised free surface, pushed by a high pressure region, as shown in Fig.11(a), beneath the bubble. During the collapse phase until the time t = 1.337 , a downward directed jet is formed at the top bubble, the region of the free surface close to the axis of symmetry rises continually, and the far region falls. Thus a free surface spike forms. The buoyancy force is equal to that given in Part 4.1, while the jet direction is different, which means that the shorter standoff distance, the larger Bjerknes force is generated by the free surface. From the pressure contours in Fig.11(b), an intermediate

60

sponding motion of the bubble centroid. During the whole process, the height of the spike grows continually. The bubble centroid moves to the free surface until the time E and goes away until the second minimum volume at the time G. The bubble centroid is repelled by the free surface sharply in view of the momentum theorem at the time H, then the bubble moves upward again under the effect of the buoyancy force. This is the second motion tendency of the free surface spike and the bubble centroid. Comparing Fig.8(b) with Fig.12(b), it can be found that there is no significant change of the maximum bubble volume when γ decreases from 1.73 to 0.8. Only the second maximum volume at the time C is 2.588 for γ = 1.73 , while it is 3.51 for γ = 0.8 at the time I. This difference is due to the standoff distance for the second bubble pulsation: for γ = 1.73 , this distance of the toroidal bubble is 0.95, while for γ = 0.8 it is 1.01. It means that for γ = 0.8 the standoff point at the beginning of the rebounding phase is far away from the free surface.

Fig.11 For γ = 0.8 , pressure contours at representative times

pressure is seen above the bubble, which drives the spike to grow continually. After the jet penetrates the bubble wall, the toroidal bubble continues collapsing to the second minimum volume. At the time t = 1.464 , the high pressure around the toroidal bubble promotes the rebounding of the bubble and the elevation of the free surface, as shown in Fig.11(c). At the time t = 1.680, the spray dome looks like a crown and keeps on growing stably. In the later stage of the rebounding phase at the time t = 2.285 , the spray dome becomes unstable, followed by the venting phenomenon. This kind of “crown” shape of the free surface can also be seen in the experiment of Zhang[16]. Figure 12(a) gives the time history of the height of the spray dome, and Fig.12(b) shows the variation of the non-dimensional bubble volume and the corre-

Fig.12 Time history of the height of the spike, bubble volume and centroid. The time corresponding to the capital letter E to I are t = 0.653, 0.72, 1.464, 1.514, 2.19, respectively

The trajectories of the bubble upper and lower points versus the time are shown in Fig.13(a). In the later stage of the expansion phase, the upper pole is greater than 1 as it gets sucked up into the free surface. Under the repelling effect of the Bjerknes force, the

61

upper point moves downward quickly during the contraction phase. The joint point of the upper and lower points is lower than the initial position, and the vertical coordinate value is ‒0.21. This position will affect the second pulsation of the toroidal bubble. As for the time history curve of the Kelvin impulse shown in Fig.13(b), the greatest magnitude is reached when the bubble is close to the free surface. The value of this curve is negative, which means that the bubble migration and the jet formation are directed away from the free surface.

4.3 The effect of standoff distance The dynamic behavior of the bubble pulsation near a free surface is determined by the buoyancy parameter, the strength parameter and the Bjerknes force related to the standoff distance. The motion of a single bubble generated beneath a free surface is considered. In order to study the effects of the standoff distance, the buoyancy and the strength parameters are kept unchanged as 0.533 and 362.4, respectively. As there is no significant change during the expansion phase, we only consider the four typical patterns of the bubble during the collapse phase, as shown in Fig.14. For the interaction between a bubble and a free surface, generally, the downward directed jet emerges due to the Bjerknes force, which is inversely proportional to the standoff distance, the upward directed jet is generated because of the buoyancy force, which is proportional to the buoyancy parameter. The two jets are in a reverse tendency: if one is large enough, the other one will disappear. With the standoff distance reduced, it can be seen from Fig.14 that the Bjerknes force grows, the downward directed jet enhances, the height of the free surface spike rises, the formation of the toroidal bubble is delayed and the width of the spray dome narrows.

Fig.13 Time history of the displacement of the upper and the lower bubble points and the Kelvin impulse Fig.15 Free surface hump profiles for various standoff distances

Fig.14 Jet profiles for various standoff distances

Figure 15 shows the four typical shapes of the bubble and the free surface for various standoff distances during the rebounding phase. Figure 15(a) gives the “plateau” shape, Fig.15(b) gives the “plateau” with a cusp, Fig.15(c) gives the “plateau” with a hump, and Fig.15(d) gives the “crown” shape or the “plateau” with a water column. The time history of the height of the spike for different standoff distances ranging from 0.8 to 2.5 is shown in Fig.16(a). When γ > 1 , the height of the spike increases during the expansion phase, decreases when the bubble collapses, and rises again when the toroidal bubble rebounds, where the motion tendency looks similar to that in Fig.8(a). When γ < 1 , the trend

62

first bubble pulsation decreases, which is mainly determined by the motion of the upper bubble pole. Figure 17 shows the volume of the bubbles considered in Fig.16(b). As the standoff distance γ increases, the expansion bubble becomes in greater extent spherical and the maximum non-dimensional volume tends to be 4π / 3 . The crucial difference happens in the second period of the bubble pulsation, where the second maximum volume increases with the reduction of γ . The reason for this phenomenon is discussed in Part 4.2.

Fig.16 Time history of the height of the spike and the displacement of the upper and lower bubble points for different standoff distances

Fig.18 Time history of the bubble centroid and the Kelvin impulse for different standoff distances

Fig.17 Time history of the bubble volumes for different standoff distances

of the height curve changes, and it goes up straight forward as shown in Fig.12(a). When γ = 1 , the tendency of the height curve is transitional, it does not have the second smallest elevation as shown in Fig.8(a), or it grows continuously as shown in Fig.12(a). As Fig.16(b) illuminates, during the expansion phase, the outward displacements in both upper and lower points increase initially. Once the collapse phase begins, the upper pole will change sharply because of the repulsion effect of the free surface. With the standoff distance γ decreasing, the repulsion on the upper bubble pole becomes severe, while γ has little influence on the lower bubble pole, the period of the

The corresponding motions of the bubble centroid are shown in Fig.18(a). The bubble centroid generally moves to the free surface during the early stage of the expansion phase. During the collapse phase, when γ ≥ 1 , the centroid continues moving upward, while the others move downward at the same time. In the rebounding phase, the bubble centroid keeps moving to the free surface. With regard to the centroid movement in these two phases, it is found that at about t = 1.5 in the early stage of the rebounding phase, the centroid moves more sharply than at the other times in view of the momentum theorem. The difference of the movement directions can be explained by the sign change of the Kelvin impulse shown in Fig.18(b). For γ ≥ 1 , the Kelvin impulse is negative in the early stage of expansion and turns positive in the subsequent phases, which indicates the

63

change in the direction of the fluid momentum of the bubble. 4.4 The effect of the buoyancy parameter When the bubble volume is large enough, for example, the underwater bubble in a large-scale explosion, the buoyancy force may have a significant effect on the motion of the bubble and the free surface spike. In this part, in order to investigate the effect of the buoyancy force, the buoyancy parameter δ ranges from 0 to 0.53, the standoff distance γ is 0.8 below the free surface and the strength parameter ε is still equal to 362.4.

Fig.19 Time history of the spike height with different buoyancy parameters ranging from 0 to 0.53

As Fig.19 illustrates, the spike height grows with the decrease of the buoyancy parameter δ . Generally speaking, the buoyancy force is upward, and the large buoyancy will promote the bubble to move upward. Thus the height of the spray dome should be high for a large buoyancy. The reason for this phenomenon as shown in Fig.19 is further explained in Figs.20 and 21. For the constant standoff distance γ = 0.8 , the motion tendency of the spike height curve in Fig.19 looks similar to that in Fig.12(a), which means that the variation of the buoyancy parameter δ does not change the motion tendency of the time history curve of the spike height.

Fig.21 The representative profiles of the bubble and the free surface for various buoyancy parameters

According to Fig.20, both the width and the height of the spray dome increase as the buoyancy parameter δ increases. H d is the height of the spray dome, and Wd is the width of the spray dome. The representative shapes of the bubble and the free surface at about time t = 1.9 are shown in Fig.21. The height of the free surface spike at δ = 0 is higher than that at δ = 0.533 , and the height and the width of the spray dome at δ = 0 are obviously smaller than those at δ = 0.533 . It means that a large buoyancy is used to elevate the spray dome rather than the free surface spike.

Fig.20 Time history of the width and the height of the spray dome with different buoyancy parameters ranging from 0 to 0.53

Fig.22 Time history of the bubble centroid with different buoyancy parameters

64

The corresponding motion of the bubble centroid is shown in Fig.22. With the buoyancy parameter δ increasing, the upward movement of the bubble centroid enhances. The motion tendency of the bubble centroid looks like that in Fig.12(b), which means that the variation of the buoyancy parameter δ does not change the motion tendency of the time history curve of the bubble centroid. At about time t = 1.5 , the bubble centroid decreases sharply with the decrease of δ. 4.5 The effect of strength parameter In this part, in order to investigate the effect of the initial pressure, simulations are carried out with the strength parameter ε ranging from 200 to 500, the standoff distance γ of 1.73 below the free surface and the buoyancy parameter δ of 0.533.

Fig.23 Time history of the height of the spike and the bubble centroid with different strength parameters

Figure 23(a) shows the time history of the height of the spike. For the constant standoff distance γ = 1.73, the motion tendency of the height curve looks like that in Fig.8(a). The related bubble centroids in the above cases are depicted in Fig.23(b). The motion tendency of the bubble center looks similar to that in Fig.12(b). With the strength parameter going up, some numerical characteristics including the height of the free surface spike, the valley of the free surface, and the upward movement of the bubble centroid are increased. It means that the variation of the strength pa-

rameter ε does not change the motion tendency of the time history curve of the height of the spike and the bubble centroid. 5. Discussions and conclusions The BIM is a robust and commonly used method, but the 3-D BIM has not been given enough attention in simulating the evolution of the free surface spike characteristics induced by the motions of a non-spherical rebounding bubble. The present 3-D BIM result is well compared with both validated experimental data and axisymmetric BIM results. Many significant physical features are seized by our numerical simulation. During the first bubble pulsation period, the present code gives results more close to the experiment data than the axisymmetric BIM. A difference between the two numerical methods is seen during the rebounding phase, which may result from the different smoothing scheme and the EMT technique. Two different spike and centroid motion tendencies, the variations of the bubble volume amplitude during the rebounding phase and the numerical spike characteristics are caught and discussed. The effects of the standoff distance, the buoyancy force and the strength parameter on the dynamic behavior of the bubble pulsation near a free surface are analyzed. (1) When the standoff distance is large, the first bubble centroid motion tendency is the sharp motions toward the free surface in the early stage of the rebounding. For the free surface, as the first spike motion tendency, it rises with the expansion of the bubble, falls with collapse and rises again during the rebounding phase. The free surface looks like a “plateau” as the symmetry axis region of the free surface rises more slowly than the region close by. (2) When the bubble is close to the free surface, one sees the second bubble centroid motion tendency, drastically far away from the free surface. As the second spike motion tendency, it rises all the time, the free surface looks like a crown, and the boundary of the spray dome is more instable than the water column. As can be seen in our simulation, γ = 1 is the demarcation line of the two free surface tendencies. It means that, from the free surface spike motion tendency, a preliminary judgment can be made whether the standoff distance is larger than 1. (3) The first maximum volume is in direct proportion to the standoff distance γ , while the second maximum volume is inversely proportional to the standoff distance γ due to the difference of the secondary standoff distance for the toroidal bubble. (4) Four typical patterns of the bubble jet during the contraction phase and four typical profiles of the free surface during the rebounding phase for the various standoff distances are shown. It means that, from

65

the shape of the spray dome, the standoff distance can be estimated. (5) The buoyancy force and the strength parameter have significant effects on the bubble motion. But neither of them changes the spike motion tendency and the bubble centroid motion tendency. When the buoyancy force is large, the height of the spike decreases. This phenomenon is caused by the fact that the buoyancy force mainly elevates the width and the height of the spray dome rather than the free surface spike. References [1] BLAKE J. R., TAIB B. B. and DOHERTY G. Transient cavities near boundaries. Part 2. Free surface[J]. Journal of Fluid Mechanics, 1987, 181: 197-212. [2] WANG Q., YEO K. S. and KHOO B. C. et al. Strong interaction between a buoyancy bubble and a free surface[J]. Theoretical and Computational Fluid Dynamics, 1996, 8(1): 73-88. [3] KLASEBOER E., KHOO B. C. and HUNG K. C. Dynamics of an oscillating bubble near a floating structure[J]. Journal of Fluids and Structures, 2005, 21(4): 395-412. [4] WILKERSON S. A. A boundary integral approach to three dimensional underwater explosion bubble dynamics[D]. Doctoral Thesis, Baltimore, USA: Johns Hopkins University, 1990. [5] ZHANG Y. L., YEO K. S. and KHOO B. C. et al. Threedimensional computation of bubbles near a free surface[J]. Journal of Computational Physics, 1998, 146(1): 105123. [6] PEARSON A., COX E. and BLAKE J. R., et al. Bubble interactions near a free surface[J]. Engineering Analysis with Boundary Elements, 2004, 28(4): 295-313. [7] LI Zhang-rui, SUN Lei and ZONG Zhi et al. A boundary element method for the simulation of non-spherical bubbles and their interactions near a free surface[J]. Acta Mechanica Sinica, 2012, 28(1): 51-65.

[8] KLASEBOER E., MANICA R. and CHAN D. Y. C. et al. BEM simulations of potential flow with viscous effects as applied to a rising bubble[J]. Engineering Analysis with Boundary Elements, 2011, 35(3): 489-494. [9] WANG Q. Non-spherical bubble dynamics of underwater explosions in a compressible fluid[J]. Physics of Fluids, 2013, 25(7): 072104. [10] WANG Qian-xi, YANG Yuan-xiang and TAN Danielle Sweimann et al. Non-spherical multi-oscillations of a bubble in a compressible liquid[J]. Journal of Hydrodynamics, 2014, 26(6): 848-855. [11] PEARSON A., BLAKE J. R. and OTTO S. R. Jets in bubbles[J]. Journal of Engineering Mathematics, 2004, 48(3-4): 391-412. [12] ZHANG Y. L., YEO K. S. and KHOO B. C. 3D jet impact and toroidal bubbles[J]. Journal of Computational Physics, 2001, 166(2): 336-360. [13] WANG C., KHOO B. C. and YEO K. S. Elastic mesh technique for 3D BIM simulation with an application to underwater explosion bubble dynamics[J]. Computers and fluids, 2003, 32(9): 1195-1212. [14] KLASEBOER E., FERNANDEZ C. R. and KHOO B. C. A note on true desingularisation of boundary integral methods for three-dimensional potential problems[J]. Engineering Analysis With Boundary Elements, 2009, 33(6): 796-801. [15] LI Zhang-rui, SUN Lei and ZONG Zhi et al. Some dynamical characteristics of a non-spherical bubble in proximity to a free surface[J]. Acta Mechanica, 2012, 223(11): 2331-2355. [16] ZHANG A-man, WANG Chao and WANG Shi-ping et al. Experimental study of interaction between bubble and free surface[J]. Acta Physica Sinica, 2012, 61(8): 084701(in Chinese). [17] ZONG Zhi, HE Liang and SUN Long-quan. Numerical study of loading on the surface ship near an underwater explosion bubble[J]. Journal of ship Mechanics, 2008, 12(5): 733-739(in Chinese). [18] ZHANG A., YAO X. and YU X. The dynamics of threedimensional underwater explosion bubble[J]. Journal of Sound and Vibration, 2008, 311(3): 1196-1212.