Investigation of free-field underwater explosion with Eulerian finite element method

Ocean Engineering 166 (2018) 182–190

Contents lists available at ScienceDirect

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

Investigation of free-field underwater explosion with Eulerian finite element method

T

Yunlong Liu, A-Man Zhang∗, Zhaoli Tian, Shiping Wang College of Shipbuilding Engineering, Harbin Engineering University, Harbin, 150001, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Underwater explosion Eulerian finite element method Shockwave Bubble pulsation

Underwater explosion in free field is investigated in this paper numerically. With the multiphase interface captured by the Volume of Fluid (VOF) method, the Eulerian finite element method (EFEM) for underwater explosion is established. The MUSCL algorithm to advect the element centered variables is modified to maintain its conservativeness in the axisymmetric model, and the second order Early Time Approximation (ETA2) is adopted to deal with the non-reflecting boundary condition of the computational domain. By comparing with the experiment results, the numerical model is proved in both shockwave and bubble oscillation simulation. With the numerical model, 2 underwater explosion cases with different initial depths are simulated and the near-field pressure characteristics are analyzed. The near-field bubble pulsation load exhibits noteworthy asymmetry around the bubble, which is attributed to the 2 high pressure regions above and beneath the bubble, respectively. By comparing the evolution of bubbles with different buoyancy parameters, it is found that the increase of the buoyancy parameter can delay the penetration and extend the collapse of the bubble.

1. Introduction The underwater explosion is an important research topic in both military and civilian fields. The pressure load caused by the underwater explosion can lead to much more serious damage on a nearby structure than the same explosive detonated in air because the water is more capable in transmitting the shockwave load. Furthermore, the bubble consisting of the gaseous productions will move violently in the water after the shockwave radiation. If the explosive is placed at a proper location, the bubble may impact the nearby structure with its highspeed liquid jet and the secondary bubble pulsation load. In the past few decades, many researches have devoted themselves into the studies on the underwater explosion. Cole (Cole, 1948) summarized the proceeding works on the underwater explosion and established the main framework of this subject, which had been a significant guidance to the following studies. Because of the large ratios in both time and space between the shockwave phase and the bubble motion phase, there're great difficulties in simulating the whole process of the underwater explosion continuously, especially in times that the computational resources are extremely limited. The shockwave phase of the underwater explosion is used to being predicted by the empirical formulas fitted by the experiments data (Cole, 1948; Geers and Hunter, 2002; Price, 1979). For the far-field pressure profile, many formulas have enough

∗

accuracy for the engineering application because the disturbance of the fluid is so small that the problem can be easily linearized. However, when the near-field pressure is concerned, most empirical formulas lose accuracy because the nonlinearity arises (Jin et al., 2017; Ming et al., 2016). Bubble motion is one of the most important aspects in which the underwater explosion differs from the air blast. The bubble oscillates in the water and radiates bubble pulsation load outward every time it reaches its minimum volume. Meanwhile, the bubble will interact with the nearby boundaries and develop high speed jet. Extensive studies have been carried out on the bubble dynamics. Geers and Hunter established a spherical bubble pulsation model based on the Doubly Asymptotic Approximation method considering the linear wave effects of the internal gas and the outside water (Geers, 1978; Geers and Park, 2005). Their works have been widely used in predicting the far-field underwater explosion load including the bubble oscillation phase. For the non-spherical bubble dynamics, boundary element method has been employed as one of the most popular methods (Blake and Gibson, 1987; Klaseboer et al., 2005a; Liu et al., 2016a; Pearson et al., 2004; Wang and Khoo, 2004; Wang et al., 2003; Wang, 2013; Zhang and Liu, 2015). However, the incompressible or weak compressible assumption of the fluid limits it in the simulation of the shockwave phase, and the energy lose mechanism of the bubble oscillation is difficult to simulate

Corresponding author. E-mail address: [email protected] (A.-M. Zhang).

https://doi.org/10.1016/j.oceaneng.2018.08.001 Received 10 October 2017; Received in revised form 26 May 2018; Accepted 3 August 2018 Available online 13 August 2018 0029-8018/ © 2018 Elsevier Ltd. All rights reserved.

Ocean Engineering 166 (2018) 182–190

Y. Liu et al.

accurately. Liu et al. adopted the finite volume method with the interface captured by front tracking method to study the collapse of a bubble near a solid wall (Liu et al., 2017). Wang et al. used the Discontinuous Galerkin Method to consider the compressibility of the fluid and preserved the sharp pressure profile generated by the jet impact (Long-kan et al., 2016). Liu et al. presented a simplification model to predict the far-field bubble oscillation considering the asymmetry caused by the migration and the jet impact of the bubble (Liu et al., 2014). Numerical investigation was carried out by Barras et al. with the ALE method on the pulsating bubble phenomenon of the underwater explosion (Barras et al., 2012). Zhang et al. developed a fully coupled fluid-structure interaction numerical model to investigate the damage of an underwater explosion bubble on the submerged structures (Zhang et al., 2017). However, these works seldom simulated the shockwave and the bubble oscillation phase continuously. The Eulerian finite element method (EFEM) is an effective approach in computational fluid dynamics which has been widely used in the hydrocodes (Benson, 1992a; Benson and Okazawa, 2004; Tian et al., 2018). The operator split technique makes it flexible to be combined with other numerical techniques. In this paper, the underwater explosion in the free-field is simulated with the EFEM. The paper is organized as follows. Firstly, the basic theoretical and numerical models are presented, and the special treatments of the model is elaborated in details. Then 2 cases of underwater explosion are investigated with the numerical model which is validated by the experiment. Then some conclusions are drawn to offer reference to the relevant studies.

Bubble surface

z

Gravity

2a

r

o Initial explosive

Computational domain a

2. Theoretical and numerical method 2.1. Eulerian finite element method with operator split In this paper, the fluid field of underwater explosion is concerned. If the explosive is spherical and other fluid boundaries are ignored, the fluid flow motion in the gravity field can be treated as axisymmetric. Therefore, a cylindrical coordinate system is established with its origin locating at the center of the explosive and the revolving axis pointing to the opposite direction of the gravity, as shown in Fig. 1 r and z denote the radial and the axial coordinates, respectively. The computational domain is taken as a rectangle sized by a × 2a. The EFEM is adopted in this paper to simulate the underwater explosion with proper efficiency and accuracy. The detailed main framework of this method can refer to the published papers (Benson, 1992a, 1992b; 2008; Benson and Okazawa, 2004). The Reynolds number of the fluid field around the underwater explosion is usually high enough to ignore the viscosity (Best, 2002; Klaseboer et al., 2005b; Liu et al., 2016b; Wang and Khoo, 2004; Zhang and Liu, 2015). Then the stress tensor can be simplified as a scalar of pressure and satisfies the momentum conservation equations,

∂U + ∇p + (u˙ ⋅∇) U = ρ g, ∂t

Fig. 1. Configuration of the underwater explosion bubble motion in free-field.

term in Eq. (2), it is based on the perspective of Lagrange and easy to be solved with traditional FEM. In the present cylindrical coordinate system, it can be written in the form of weighted residual method with v denoting the weight function

∬ r (ρu¨ v − p∇v) dΩ = − ∫ rpnˆ vdΓ + ∬ r gvdΩ + ∬ pv∇rdΩ, Ω

Ω

where the integral by parts and the Gauss-Green formula are used to eliminate the calculation of the gradient of p (Wu and Gu, 2012). nˆ is the unit normal vector at the boundary Γ of the integral domain Ω pointing outward of Ω. Taking Φ as the shape function of the element and choose it as the weight function, Eq. (4) can be rewritten in the semi-discretized form,

(1)

∑ ∫ (rρΦM ΦN )dΩu¨ N = ∫ (r gΦM + rp∇ΦM )dΩ− ∫ rpnˆ ΦM dΓ N

Ω

Ω

+

Γ

∬ pΦM ∇rdΩ, Ω

(5)

¨ N is the unknown nodal achere M and N are the nodal number and u celeration vector calculated by solving Eq. (5) efficiently with the lumped mass matrix. Then the nodal velocity and displacement are advanced explicitly with the second order central difference method,

(2)

1

1

u˙ n + 2 = u˙ n − 2 + Δt u ¨ n,

and

∂U + (u˙ ⋅∇) U = 0. ∂t

Ω

(4)

where u˙ is the velocity vector and U = ρ u˙ is the momentum. ρ and p are the density and the pressure of the fluid. g is the gravity acceleration. The third term in Eq. (1) is the convection term and requires upwind scheme to be solved which is hard to implement in traditional FEM. Thus, with the operator split technique, Eq. (1) is decomposed into two equations,

∂U + ∇p = ρ g ∂t

Γ

1

un + 1 = un + Δt u˙ n + 2 , (3)

(6) (7)

where the superscript stands for the number of the increment, Δt is the time increment which is limited by Δt < L/ C0 , L is the smallest edge length of the elements, and C0 is the dilatational wave speed. After the

After these two equations are solved, the increments of U are added together to update the real U. Because of the absence of the convection 183

Ocean Engineering 166 (2018) 182–190

Y. Liu et al.

conservativeness requires that the integration over each element with the linear distribution defined by the slopes equals to the integration with the constant distribution. However, the MUSCL algorithm must be modified to maintain conservativeness in the axisymmetric model because of the coefficient r in the volumetric integration in the cylindrical coordinate system. Thus, for an element centered variable ϕ, a correction ϕc must be considered to preserve the conservativeness

nodal position is calculated, the density and the internal energy of each element can be updated as follows,

∂ρ + ρ∇⋅u˙ = 0, ∂t

(8)

∂ρEm + p∇⋅u˙ = 0, ∂t

(9)

where Em is the internal energy per unit mass and the fluid is assumed to be adiabat because the time period of the underwater explosion is short enough to neglect the thermal conduction. The equation of state (EOS) of the fluid outside of the bubble is chosen as the Tait equation (Qiu et al., 2008),

p = ρEm (γ − 1) − γPw ,

∫ r [ϕ + ϕc + (p − pe)⋅∇ϕ]ds = ∫ rϕe ds, Se

(10)

p dp = ⎜⎛Em + ⎟⎞ (γ − 1) ρ⎠ dρ ⎝

ϕc

Se

Se

(16)

For the structured mesh with the element edges parallel to the axis and bounded by r = [r1, r2] and z = [z1, z2], ϕc can be simply calculated with r

The pressure of the internal gaseous products of the explosive are modeled with the Jones-Wilkens-Lee (JWL) equation (Lee et al., 1968),

ϕc =

1 er ⋅∇ϕ 1 1 ⎛ r 3 − rc r 2⎞ . 2 rc (r1 − r2 ) ⎝ 3 ⎠ r2

(17)

In the present model, the momentum is stored at the nodes. Thus, it cannot be transported with the MUSCL algorithm. The half index shift algorithm (Benson, 1992b, 2008) is adopted in this paper to transport the momentum.

(12)

where A, B, R1, R2 and ω are the material constants; ρ0 and ρ are the densities of the explosive and the gaseous products. e1 and e2 are defined as

( ) = exp ( − R ) .

∫ r ds = ∫ rϕe ds − ∫ r [ϕe + (p − pe)⋅∇ϕ]ds. Se

(11)

ωρ ⎞ ωρ ⎞ ⎛ p = A ⎜⎛1 − ⎟ e1 + B ⎜1 − ⎟ e2 + ωρEm R2 ρ0 ⎠ R1 ρ0 ⎠ ⎝ ⎝

(15)

where the subscript ϕe is the value at the element center of element Se, ∇ϕ is the gradient of ϕ at the element center determined by the calculated slopes, p is the coordinate vector of the integration point and pe is the coordinate of the element center. Then ϕc can be expressed as

where γ is the ratio of specific heat, and Pw is the reference pressure. Combining Eq. (9) with Eq. (10), we can easily obtain the dilatational wave speed C0

C02 =

Se

2.2. Non-reflecting boundary condition

ρ

e1 = exp − ρ0 R1 e2

ρ0 ρ

2

For the underwater explosion simulation, the non-reflecting boundary condition is essential for the accuracy. The non-reflecting boundary condition should absorb the outward pressure wave effectively to avoid the spurious reflection. At the same time, the long-term effect on the computational domain must be minimized to preserve the accurate oscillation of the bubble. To archive the non-reflecting effect, a proper surface pressure must be applied to the boundary of the computational domain and taken into account in the second integration at the right side of Eq. (5). According to the Bernoulli equation for irrotational flow, we can decompose the total pressure into three parts:

(13)

Similarly, the dilatation wave speed C for the gaseous products can be written as

ρ0 Ae1 Be2 ⎞ C02 = ω (1 + ω) ⎜⎛Em − − (AR1 e1 + BR2 e2). ⎟ + 2 R ρ ρ R ρ 1 0 0 2⎠ ⎝

(14)

The material constants used in this paper are listed in Table 1. After the mesh and the variables are advanced forward in a time increment, the mesh is moved back to its original position to eliminate the element distortion and to solve Eq. (3). This procedure refers to the Eulerian phase which considers the advection quantities between adjacent elements during the moving back of the mesh. The convected volume is calculated with the Volume of Fluid (VOF) (Benson, 1992a; Hirt and Nichols, 1981) method with the prescribed velocity calculated in the Lagrangian phase. Then the element centered variables, such as the transported mass and internal energy, are updated subsequently. To archive a second order accuracy, the MUSCL algorithm (Benson, 1992a; Van Leer, 1977) is adopted to transport the element centered variables. Conventionally, the MUSCL algorithm takes the element centered variable as the average value of the element. Then, by introducing a second order distribution among the adjacent elements in each direction, the gradient of the variable is calculated as the slope at the element center with specific slope limiters. Theoretically, the

p = p∞ −

TNT productions

Pw MPa

γ

ρ0 Kg/m3

330.9

7.15

1000

A GPa 371.2

B GPa 3.2

R1 4.15

R2 0.95

ω 0.3

E0 MJ/kg 4.3

(18)

where the first term p∞ = patm + ρgh is the undisturbed pressure at the same depth. The second term is nonlinear and decreases along with R−2and R is the distance to the center of the bubble. If the computational domain is taken large enough, this term is negligible for the outer boundary. The third term pd in Eq. (18) is the transient term. Based on the linearized acoustic assumption which is reasonable when the computational domain is sufficient large, the dynamic pressure pd and the material velocity u˙ satisfy the second order Early-Time Approximation (ETA2) equation (Felippa, 1980)

pd + κp c

∫ pd dt = ρcu˙ ⋅e,

(19)

where e is the unit vector pointing to the direction that the acoustic wave propagating, and κp is the local mean curvature of the surface S which is taken as the wave front in this paper. As we can see that, the second term at the left side of eq. (19) is the effect of the curvature. Because the boundary of the computational domain is not always coinciding with the wave front surface, the calculation of the curvature κp is important to the accuracy. If the bubble moves in an infinite field, no physical reflection wave should be observed at the non-reflecting boundary. Thus, all the incident waves come directly from the bubble and κp can be taken as 1/R for simplicity, where R is the distance to the bubble center. Considering the size of the computational domain should

Table 1 Material constants for Tait equation and JWL equation (Dobratz and Crawford, 1981; Qiu et al., 2008). Water

1 ρ u˙ 2 + pd , 2

ρ0 Kg/m3 1630

184

Ocean Engineering 166 (2018) 182–190

Y. Liu et al.

0.02

be much larger than the size of the bubble, this simplification has a second order accuracy.

0.019 0.018

3. Results and discussion

0.017

In free field, the jet evolution due to the buoyancy is an important feature except for the pulsation. This feature can be represented by the non-dimensional buoyancy parameter (Klaseboer et al., 2005a)

0.015 0.014

(20)

0.013

When δ is 0, there is no buoyancy effect and the bubble tends to pulsate spherically. With δ increasing, an upward jet would develop during the collapsing phase of a bubble with a great δ. In this paper, 2 groups of cases are investigated with different buoyancy parameters.

0.012 0.011 0.2

θ=

where W is the explosive weight with the unit of kilogram. Specifically, Rm = 0.988 m for this case. Thus, the buoyancy parameter δ is small and equals 0.22. Because Rm is much smaller than the tank size and the initial depth, effects of the free surface and the tank wall are ignored in the simulation. The computational domain is taken as a rectangle with a = 5 Rm as shown in Fig. 1, and is discretized into square elements with the length Le = 0.01 Rm. The simulation is performed on a PC with Core i7-4790 (3.6 GHz) and 16 GB memery. OpenMP is used in the code to accelerate the simulation. The simulation takes about 5 h. The pressure contour of the numerical results during the shockwave phase is shown in Fig. 2. The results show that the shockwave radiated from the detonated explosive propagates outward spherically. A sharp front of the shockwave is preserved in the simulation. The shockwave reaches the boundary of the computational domain at t = 3.04 ms. Little spurious reflection wave emerges, which is attributed to the non-reflecting boundary condition based on ETA2 as shown in Fig. 2(e). 2 pressure sensors are placed 3 m away from the explosive at the same depth in the

0.5

(22)

p / Pa 2.72ms

z/m

1.94ms

z/m r/m

pr − p∞ , pi − p∞

p / Pa

1.48ms

z/m

z/m

r/m

0.45

where pi and pr are the amplitudes of the incident pressure and the reflected pressure at the boundary. Fig. 3 shows the change of θ with Rm/a. We can see that the reflection parameter θ is linearly related to Rm/a. Considering the incident wave pi is nearly proportional to 1/a, the computational results of the reflected wave converge to θ with the same rate as 1/a2. The pressure profiles along time are compared in Fig. 4 for shockwave which are from the present model with different element size Le, the Geers-Hunter model (Geers and Hunter, 2002; Geers and Park, 2005) and the experiment. We can see that with the decrease of Le, the shock wave profile from the present model gets more sharp and converges to the experimental one. In spite of that the numerical model smears the shockwave a little, the profile of the shockwave pressure calculated by the numerical model with Le = 0.01Rm agrees with the experimental results well as shown in Fig. 4. Obviously, the GeersHunter model based on the linear assumption fails in predicting the arrival time of shockwave because the wave speed is assumed constant

p / Pa

1.02ms

0.4

simulation and the experiment, respectively. To demonstrate the convergence of the numerical model with respect to the domain size a, 7 simulation results are compared together in which a is taken from 2Rm to 5Rm. Here, we define a non-dimensional parameter θ to represent the reflection of the present non-reflecting boundary:

(21)

p / Pa

0.35

Fig. 3. Change of the reflection parameter θ with Rm/a.

Firstly, cases with small buoyancy parameters are considered. In this case, a TNT explosive weighted 0.5 Kg is detonated at the depth of 10 m. An experiment is carried out in a cylindrical water tank with a radius of 24 m and a depth of 20 m for comparison. The maximum radius of the bubble can be calculated with the empirical formula (Klaseboer et al., 2005a)

p / Pa

0.3

Rm /a

3.1. Model validation and asymmetry of bubble pulsation load

W ⎞1/3 Rm = 3.38⎛ , ⎝ h + 10 ⎠

0.25

r/m

3.67ms

z/m

δ=

ρgRm . p∞

0.016

r/m

r/m

Fig. 2. Pressure contour of the numerical results for the shockwave propagation with a TNT explosive weighted 0.5 Kg detonated at the depth of 10 m. The units of length and pressure are m and Pa, respectively. The blue line represents the bubble wall. 185

Ocean Engineering 166 (2018) 182–190

Y. Liu et al.

12

106 Experiment EFEM, Le = 0.01 R m

10

EFEM, Le = 0.02 R m EFEM, L = 0.03 R e

8

p/ Pa

in the experiment. Compared with the Geers-Hunter model, the present model preserves more accuracy in predicting the bubble period. Furthermore, the bubble pulsation load is treated as spherical symmetric in the Geers-Hunter model, while it's not true for the near-field. It's easy to find the center of the bubble pulsation pressure wave is at about [0, 0.5] because the upward migration of the bubble due to the buoyancy. The pressure profiles in the 3 directions, which are upward, downward, and horizontal directions, at t = 136 ms and 137 ms are drawn in Fig. 7 to investigate the distribution of the bubble pulsation load in different directions. Apparently, the bubble pulsation load at the downward direction is greater than the others, which can be attributed to the high-pressure region beneath the bubble before the penetration which breaks the symmetry of the pressure field around the bubble. Furthermore, the bubble plays a role in obstructing the pressure wave propagating upward. Previous study (Liu et al., 2014)has pointed that the asymmetric difference between the pressure peaks in different directions is a second order infinitesimal with respect to R−1 and decays with propagation time or distance rapidly. Thus, in the far-field problems, the difference can be ignored and the Geers-Hunter model can predict accurately. We also should notice that there is a smaller pressure peak at the left side of the bubble pulsation pressure peak, which is only found in pressure profiles in the upward direction. This pressure peak is caused by the high-pressure region above the bubble after the penetration. Because the bubble obstructs the propagation downward of this pressure wave, the pressure profile in this direction drops fastest after bubble pulsation load passed. As the initial depth of the explosive is increased to 15 m, the simulation results are shown in Fig. 8. Because the explosive depth has little influence on the shockwave phase and the bubble expanding phase except the bubble pulsation period, Fig. 8 only presents the bubble collapsing from its maximum volume at t = 55.3 ms which is 11.4 ms earlier than the previous case. Similar bubble dynamics is shown despite the preceding in time. The jet penetrates the bubble at t = 112.1 ms which is 22.6 ms earlier than the first case. As shown in Eq. (21), the maximum radius of the bubble decreases with the increasing of the depth. Due to the bigger ambient pressure in the deeper water, the maximum radius of the bubble is smaller, too. The changes of the equivalent radius of the bubble and the pressure 3 m away from the explosive are drawn in Fig. 9 for the comparison of these two cases. It's easy to find that the bubble radiates a pressure peak each time it reaches its minimum volume. The first one is the shockwave generated right after the detonation. The second one is the bubble pulsation load generated when the bubble collapses to its minimum volume. Each pressure peak delays a little compared to the change of the radius due to the propagation from the bubble to the sensor. In this section, the buoyancy parameters of the 2 cases are 0.220 and 0.189, respectively. With these small buoyancy parameters, we can see that the jet develops in the final stage of the collapsing phase and the penetration occurs near the moment when the bubble reaches its minimum volume. Next, we will increase the weight of the explosive to investigate the bubble dynamics with large buoyancy parameter.

m

EFEM, Le = 0.04 R m Geers-Hunter model

6

4

2

0 1.5

2

2.5

t/ s

3 10-3

Fig. 4. Comparison of the shock wave profiles of the case of Fig. 2 from experiment, the present model with different Le and the Geers-Hunter model.

and equals to the sound speed. Actually, the shockwave speed is faster than the sound speed in the near field of underwater explosion. At t = 2.05 ms, a small fluctuation is observed only in the experiment results. This is very likely caused by the wave structure of the internal gas. After the explosive is detonated at the center, a detonation wave propagates outward. When the detonation is completed and the detonation wave reaches the explosive-water interface, a rarefaction wave and a shock wave travel inward and outward from the interface. The outward wave makes up the main part of the shock wave. The inward rarefaction wave will reflect at the bubble center and radiate a pressure disturbance into the surrounding flow field which refers to the pressure fluctuation after the initial shockwave. By contrast, the detonation process is not included in the numerical model so that the pressure profile is decreasing smoothly without the fluctuation. The bubble evolution and the surrounding pressure field are shown in Fig. 5. The explosive is detonated and the gas bubble starts expanding at t = 0 ms. High pressure shockwave and bubble load radiate from the gas bubble outward. The bubble keeps expanding and reaches its maximum volume at t = 66.7 ms and starts collapsing due to the pressure difference between the internal gas and the ambient fluid. Attributing to the gradient of the ambient pressure caused by the gravity, the lower side of the bubble subjects to greater ambient pressure than the upper side of the bubble. Thus, the lower side of the bubble collapses faster. Consequently, the water converges together beneath the bubble and produces a high-pressure region, which accelerates the collapse of the bottom of the bubble further. Therefore, a high speed upward jet emerges during the collapsing phase and penetrates the bubble into toroidal at t = 134.8 ms. At the same time, a highpressure region emerges at the toroidal eye due to the jet impact, which exists during the whole rebound phase of the bubble. The bubble reaches its minimum volume during the collapsing phase at t = 135.0 ms, and the internal gas pressure reaches its peak meanwhile. A bubble pulsation pressure wave radiates away from the pressed bubble. Then the bubble starts rebounding, and the second pulsation period begins. The bubble pulsation loads of the numerical and experimental results are compared in Fig. 6. The results show good agreement in general except that the pressure peak is a little wider in the simulation. The arrival time of the pulsation load is determined by both the wave speed and the bubble period. As the distance of the pressure sensor from the bubble is relatively near, the error of the bubble period would overide that of the wave speed, and the effect of the wave speed can be neglected. The arrival time of the bubble pulsation load is 136.7 ms in the numerical results with an error of 0.2% compared with the 137.0 ms

3.2. Influence of buoyancy parameter on bubble dynamics Because gravity is negligable in the detonation and the propagation of the shockwave, the buoyancy parameter has little effect on the shockwave phase. Thus, only the bubble evolution is considered in this section. The weight of the explosive increases from 0.5 kg to 5 kg and 50 kg at the depth of 15 m to be simulated with different buoyancy parameters. The other quantities also need to be non-dimensionalized to be compared together with different buoyancy parameters. According to the published papers on bubble dynamics (Klaseboer et al., 2005a; Zhang and Liu, 2015), the length scale, density scale and pressure scale are taken as Rm, ρ and p∞, respectively. Thus, the time scale can be calculated with 186

Ocean Engineering 166 (2018) 182–190

Y. Liu et al. p / Pa

p / Pa

p / Pa

r/m

z/m r/m

r/m

p / Pa

(d) t = 66.7ms

z/m

Charge

(c) t = 20.0ms

z/m

(b) t = 2.0ms

z/m

(a) t = 0.0ms

p / Pa

p / Pa

p / Pa (e) t = 120.0ms

r/m

(f) t = 130.0ms

p / Pa (g) t = 134.0ms

(h) t = 134.8ms

z/m

z/m

z/m

z/m

Jet impact region

Liquid jet High pressure region

r/m

r/m

r/m

p / Pa

p / Pa

p / Pa

p / Pa (l) t = 150.0ms

z/m

z/m

(k) t = 140.0ms

z/m

(j) t = 137.0ms

z/m

(i) t = 135.7ms

r/m

Bubble pulsation pressure wave

r/m

r/m

r/m

r/m

Fig. 5. Bubble evolution of the explosive weighted 0.5 kg detonated at the depth of 10 m. The contour color represents the pressure with the unit of Pa, and the array of arrows represents the fluid velocity.

3

106

6

5

EFEM

4

Geers-Hunter model

3.5

2

3

p / Pa

p/ Pa

Upward direction at t = 136ms Downward direction at t = 136ms Horizontal direction at t = 136ms Upward direction at t = 137ms Downward direction at t = 137ms Horizontal direction at t = 137ms

4.5

Experiment

2.5

x 10

1.5

2.5 2

1

1.5

0.5

0.5

1

0

0 0.125

0.13

0.135

0.14

1

1.5

2

2.5

3

3.5

4

Distance from the bubble center / m

0.145

t/ s

Fig. 7. Pressure profiles along the distance from the bubble center in different directions at t = 136 ms and 137 ms, respectively.

Fig. 6. Comparison of the bubble pulsation pressure profiles of the case of Fig. 2 from experiment, the present model and the Geers-Hunter model.

187

Ocean Engineering 166 (2018) 182–190

Y. Liu et al. p / Pa

p / Pa

p / Pa (a) t = 55.3ms

(b) t = 107.0ms

p / Pa (c) t = 111.5ms

(d) t = 112.1ms

z/m

z/m

z/m

z/m

Jet impact region

Liquid jet High pressure region

r/m

r/m

r/m p / Pa

p / Pa

(f) t = 118.0ms

z/m

z/m

(e) t = 115.0ms

z/m

(d) t = 112.9ms

p / Pa (g) t = 125.0ms

z/m

p / Pa

r/m

Bubble pulsation pressure wave

r/m

r/m

r/m

r/m

Fig. 8. Bubble evolution of the explosive weighted 0.5 kg detonated at the depth of 15 m. The contour color represents the pressure with the unit of Pa, and the array of arrows represents the fluid velocity. 6

x 10 10 1

drawn in Fig. 11 in which the solid triangles denote the moments of the penetration and the vertical lines indicate the moment imaged in Fig. 10. We can see that the upward jet develops earlier in the case with a greater buoyancy parameter from the first column of Fig. 10. Because the asymmetric velocity component caused by the jet development occupys part of the kinetic energy, the pulsation process is slowed down compared with the case with lower buoyancy parameter. Thus, the moment when the bubble reaches its minimum radius delays along with the increase of the buoyancy parameter as shown in both Figs. 10 and 11. Meanwhile, the minimum radius of the bubble also increases with the buoyancy parameter, which indicates that the bubble pulsation load would be smaller in a case with a greater buoyancy parameter.

9

R/ m

0.6

0.4

Shockwave

Rm with h = 10m

7

Rm with h = 15m

6

p with h = 10m p with h = 15m

5

p / Pa

8 0.8

4 Bubble pulsation load

3 2

0.2

1 0

0

0.05

0 0.15

0.1

4. Conclusions

t/s

In this paper, the EFEM is employed to investigate the underwater explosion in the cylindrical coordinate system. The MUSCL algorithm is modified to maintain the conservativeness in the axisymmetric model, and the ETA2 is adopted to approach the second order non-reflecting boundary condition. With the presented numerical model, 2 underwater explosion cases with different initial depths are firstly simulated for both of the shockwave phase and the bubble oscillation phase continuously. Compared with the experimental results, the numerical model predicts the underwater explosion load with a good accuracy in both amplitude and time and exhibits the capability to consider the nonlinearity in the near-field for the shockwave propagation. Under the gravity effect, the bubble collapses non-spherically and developes an upward jet. The bubble pulsation load has different characteristics in different directions, especially in the near-field. The 2 high pressure regions caused by the gradient of the ambient pressure and the jet impact are the main reasons that the bubble pulsation load loses its symmetry. Because the asymmetry decays with the distance rapidly, the spherical bubble pulsation load model can predict the pressure profile accurately in far-field. By contrast, the asymmetry must be considered in the near-field bubble pulsation load simulation. Then, the effects of the buoyancy parameter on the bubble evolution are investigated by

Fig. 9. Changes of bubble radius and pressure 3 m away from the explosive weighted 0.5 kg detonated at the depth of 10 m and 15 m, respectively. Table 2 Parameters of cases with different buoyancy parameters. Case

W/Kg

h/m

Rm/m

δ

ts /s

1 2 3

0.5 5 50

15 15 15

0.917 1.977 4.259

0.191 0.280 0.411

0.058 0.126 0.271

ts = Rm

ρ . p∞

(23)

And the corresponding maximum radius Rm, the buoyancy parameter δ and the time scale for non-dimensionalization of these 3 cases are listed in Table 2. For each case, the bubble evolution of 4 moments at the same nondimensional time are compared together in Fig. 10, and the changes of the equavalent radii of the bubbles with non-dimensional time are 188

Ocean Engineering 166 (2018) 182–190

Y. Liu et al.

p/p

p/p

p/p

p/p

z / Rm

z / Rm

z / Rm

z / Rm

(a)

r / Rm

r / Rm

p/p

r / Rm p/p

p/p

r / Rm p/p

r / Rm

z / Rm

z / Rm

z / Rm

z / Rm

(b)

r / Rm

p/p

r / Rm

p/p

p/p

r / Rm

p/p

(c) t / ts= 2.019

z / Rm

t / ts= 1.993

z / Rm

z / Rm

t / ts= 1.958

z / Rm

t / ts= 1.859

r / Rm

r / Rm

r / Rm

r / Rm

Fig. 10. Evolution of bubbles generated by (a) 0.5 kg, (b) 5 kg and (c) 50 kg TNT explosives detonated at a depth of 15 m at t/ts = 1.840, 1.938, 1.973 and 2.088, respectively. 0.8

Acknowledgement (t / t s ) 2 (t / t s ) 3

(t / t s ) 1

0.7

(t / t s ) 4

0.5

The authors would like to acknowledge the support of the National Natural Science Foundation of China (Grant No.51609044, 11672082), the China Postdoctoral Science Foundation (Grant No. 2016M600244) and the Heilongjiang Postdoctoral Fund (Grant No. LBH-Z16043).

0.4

References

R/Rm

0.6

0.3 0.2 0.1 0 1.7

Barras, G., Souli, M., Aquelet, N., Couty, N., 2012. Numerical simulation of underwater explosions using an ALE method. The pulsating bubble phenomena. Ocean Eng. 41, 53–66. Benson, D.J., 1992a. Computational methods in Lagrangian and Eulerian hydrocodes. Comput. Meth. Appl. Mech. Eng. 99, 235–394. Benson, D.J., 1992b. Momentum advection on a staggered mesh. J. Comput. Phys. 100 (1), 143–162. Benson, D.J., 2008. Momentum advection on unstructured staggered quadrilateral meshes. Int. J. Numer. Meth. Eng. 75, 1549–1580. Benson, D.J., Okazawa, S., 2004. Contact in a multi-material Eulerian finite element formulation. Comput. Meth. Appl. Mech. Eng. 193, 4277–4298. Best, J.P., 2002. The Effect of Non-spherical Collapse on Determination of Explosion Bubble Parameters. DSTO Systems Sciences Laboratory, Edinburgh, Australia, pp. 1–25. Blake, J.R., Gibson, D.C., 1987. Cavitation bubbles near boundaries. Annu. Rev. Fluid Mech. 19, 99–123. Cole, R.H., 1948. Underwater Explosion. Princeton University Press, Princeton USA. Dobratz, B.M., Crawford, P.C., 1981. LLNL explosives Handbook: Properties of Chemical Explosives and Explosives Simulants. Lawrence Livermore National Lab., CA, USA. Felippa, C.A., 1980. A family of early time approximations for fluid structure interaction. J. Appl. Mech. 47, 703–708. Geers, T.L., 1978. Doubly asympotic approximation for transient motions of submerged structures. J. Acoust. Soc. Am. 64, 1500–1508.

W = 0.5kg W = 5kg W = 50kg Jet penetration

1.75

1.8

1.85

1.9

t/t

1.95

2

2.05

2.1

s

Fig. 11. Changes of bubble radii with time generated by 0.5 kg, 5 kg and 50 kg TNT explosives detonated at a depth of 15 m.

simulating 2 more cases with different explosive weights. The results show that the increase of the buoyancy parameter could delay the penetration and extend the collapse of the bubble. Because the energy is occupied by the motion of the surrounding flow when the bubble reaches its minimum radius, the jet development would weaken the pulsation load for a larger minimum radius. 189

Ocean Engineering 166 (2018) 182–190

Y. Liu et al.

Ming, F.R., Zhang, A.M., Xue, Y.Z., Wang, S.P., 2016. Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions. Ocean Eng. 117, 359–382. Pearson, A., Cox, E., Blake, J.R., Otto, S.R., 2004. Bubble interactions near a free surface. Eng. Anal. Bound. Elem. 28, 295–313. Price, R.S., 1979. Similitude Equations for Explosives Fired Underwater. NSWC. Qiu, J.X., Liu, T.G., Khoo, B.C., 2008. Simulation of compressible two-medium flow by Runge-Kutta discontinuous galerkin methods with teh ghost fluid method. Commun. Comput. Phys. 3, 479–504. Tian, Z.L., Liu, Y.L., Zhang, A.M., Wang, S.P., 2018. Analysis of breaking and re-closure of a bubble near a free surface based on the Eulerian finite element method. Comput. Fluids 170, 41–52. Van Leer, B., 1977. Towards the ultimate conservative difference scheme IV: a new approach to numerical convection. J. Comput. Phys. 23, 276–299. Wang, C., Khoo, B.C., 2004. An indirect boundary element method for three-dimensional explosion bubbles. J. Comput. Phys. 194, 451–480. 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 (9), 1195–1212. Wang, Q.X., 2013. Non-spherical bubble dynamics of underwater explosions in a compressible fluid. Phys. Fluids 25 (7), 131–144. Wu, S.R., Gu, L., 2012. Introduction to the Explicit Finite Element Method for Nonlinear Transient Dynamics. John Wiley & Sons, Hoboken, New Jersey. Zhang, A.M., Liu, Y.L., 2015. Improved three-dimensional bubble dynamics model based on boundary element method. J. Comput. Phys. 294, 208–223. Zhang, A.M., Wu, W.B., Liu, Y.L., Wang, Q.X., 2017. Nonlinear interaction between underwater explosion bubble and structure based on fully coupled model. Phys. Fluids 29 (8), 082111.

Geers, T.L., Hunter, K.S., 2002. An integrated wave-effects model for an underwater explosion bubble. J. Acoust. Soc. Am. 111 (4), 1584–1601. Geers, T.L., Park, C.K., 2005. Optimization of the G&H bubble model. Shock Vib. 12 (1), 3–8. Hirt, C.W., Nichols, B.D., 1981. Volume of Fluid method for the dynamics of free boundaries. J. Comput. Phys. 39, 201–225. Jin, Z., Yin, C., Chen, Y., Hua, H., 2017. Coupling Runge-Kutta discontinuous Galerkin method to finite element method for compressible multi-phase flow interacting with a deformable sandwich structure. Ocean Eng. 130, 597–610. Klaseboer, E., Hung, K.C., Wang, C.W., Khoo, B.C., 2005a. Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/ rigid structure. J. Fluid Mech. 53 (7), 387–413. Klaseboer, E., Khoo, B.C., Hung, K.C., 2005b. Dynamics of an oscillating bubble near a floating structure. J. Fluid Struct. 21 (4), 395–412. Lee, E.L., Hornig, H.C., Kury, J.W., 1968. Adiabatic Expansion of High Explosive Detonation Products. University of California, Livermore. Liu, L.T., Yao, X.L., Zhang, A.M., Chen, Y.Y., 2017. Numerical analysis of the jet stage of bubble near a solid wall using a front tracking method. Phys. Fluids 29, 012105. Liu, Y.L., Wang, Q.X., Wang, S.P., Zhang, A.M., 2016a. The motion of a 3D toroidal bubble and its interaction with a free surface near an inclined boundary. Phys. Fluids 28 (12), 122101. Liu, Y.L., Wang, S.P., Zhang, A.M., 2016b. Interaction between bubble and air-backed plate with circular hole. Phys. Fluids 28 (6), 1195–1212. 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. Long-kan, W., Zhi-fan, Z., Shi-ping, W., 2016. Pressure characteristics of bubble collapse near a rigid wall in compressible fluid. Appl. Ocean Res. 59, 183–192.

190