Continuous simulation of the whole process of underwater explosion based on Eulerian finite element approach

Applied Ocean Research 80 (2018) 125–135

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Continuous simulation of the whole process of underwater explosion based on Eulerian finite element approach

T

⁎

W.T. Liua, F.R. Minga, , A.M. Zhanga, X.H. Miaoa,b, Y.L. Liua a b

College of Shipbuilding Engineering, Harbin Engineering University, Harbin, Heilongjiang, China Navy Academy, Beijing, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Underwater explosion Shock wave Bubble motion Eulerian finite element formulation Continuous simulation

Underwater explosion is an extremely complex hydrodynamic problem and it is always divided into two stages and studied individually, i.e., shock wave stage and bubble motion stage. Actually underwater explosion is a continuous process and there are intensive interactions between the two stages in some cases, e.g. underwater explosion in the vicinity of the wall. In this paper, the Eulerian finite element formulation considering the fluid compressibility is applied to continuously simulate the process of underwater explosion based on the subroutine of ABAQUS. Firstly, the governing equation and its solution by the operator splitting approach are discussed. Then the free-field underwater explosion is continuously simulated, including shock wave propagation and the bubble expansion, contraction, collapse, jet and rebound. The applicability and reliability of this method are discussed by comparing with experimental results from Klaseboer [5]. Furthermore, the traditional treatments of the simulation initialized at the bubble motion stage are discussed. Finally, the underwater explosion near wall is investigated regarding the interactions of shock wave and bubble motion and the importance of the continuous simulation of underwater explosion on the bubble motion characteristics is summarized.

1. Introduction The main phenomena of underwater explosion includes shock wave formation and propagation, bubble pulsation and migration [1]. Generally, considering the differences of time sequence and time scale, the process of underwater explosion is usually divided into two stages, i.e., the shock wave stage and bubble pulsation stage, and studied individually [2–6]. At the former stage, the duration of shock wave is in milliseconds and the peak pressure can be up to the level of GPa [7]. This stage is usually featured by strong nonlinearity and the compressibility of the fluid should be considered. As the predecessors, Cole and Zamyshlyaev had systematically investigated and summarized the characteristics of the shock wave generation, propagation and the pressure history by theoretical and experimental methods [1,8]. Following their great contributions, some scholars have also carried out a series of researches on shock waves and their damage to structures by different numerical methods, e.g. Smoothed Particle Hydrodynamics [7,9–11], the modified Ghost Fluid Method [12] and so on. The understanding of this stage has been relatively in-depth. Following the shock wave propagation, the high-pressure bubble formed by the detonation products begins to expand and lasts several hundred milliseconds. That is, it comes into the bubble pulsation stage.

⁎

At this stage, the pressure of the bubble is of the magnitude of MPa and the interface will change complicatedly, especially when the jet forms under the ambient pressure difference [5,13]. There are two main popular methods for studying bubble dynamics. One is the boundary element method (BEM) [14] based on potential flow theory, which is a kind of numerical method for solving integral equations efficiently. The dimension is reduced and thus it is widely used in the field of bubble dynamics. Another is the classical method in computational fluid dynamics, i.e., the discrete solution method. Generally, nonlinear Eulerian equations solved by spatial discretization method are used to describe the flow field. The mentioned methods include finite difference method, finite volume method and finite element method. Based on the above two kinds of methods, the characteristics of underwater explosion bubbles under different boundary conditions are studied, such as free field [5,15–17], near free surface [18–22], near rigid wall [23–26] and under more complex boundaries [27–30]. Furthermore, whipping motions and other response characteristics of structures subjected to underwater explosion bubbles are also analyzed [31–34]. Undoubtedly, the above two stages of underwater explosion are a continuous process. In fact, the bubble pulsation stage has some memory effects [35] on the initial conditions, namely the explosive charge shape and the detonating style would affect the following bubble

Corresponding author. E-mail address: [email protected] (F.R. Ming).

https://doi.org/10.1016/j.apor.2018.08.016 Received 13 April 2018; Received in revised form 18 July 2018; Accepted 22 August 2018 0141-1187/ © 2018 Elsevier Ltd. All rights reserved.

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

milliseconds), so it is assumed that the fluid is adiabatic and the thermal diffusion effect is ignored in the present work. Eqs. (1)–(3) are the governing equations for describing the flow field of underwater explosion. Generally, there are two ways to solve Eqs. (1)–(3): one is to solve the non-symmetric system of equations directly and the other one is referred to as an operator splitting in [41,42], where the Lagrangian motion is decoupled from the additional mesh motion. The latter one is used in this paper because of its high efficiency, where the calculation is divided into two phase: Lagrange phase and transport phase in each time step. In the first phase, the mesh moves with the material flow, hence this phase is called as Lagrange phase. The changes in density, velocity and energy due to the actions of external and internal forces in Lagrange phase are calculated by the following equilibrium equations:

dynamic behavior. For the investigations of bubble dynamics beginning from the second stage, the memory effects are difficult to take into account. Some common tricks for the initial condition of the highpressure bubble, such as the employ of empirical formulas of spherical charge or the bubble motion equations [5,27,36], are not rigorous for the accurate simulations. On the one hand, before bubble expansion and collapse, the shock wave has propagated in water and the pressure and velocity fields different from that in still water are caused. The initial velocity of the bubble and its surrounding field is ignored compared to the bubble from the shock wave stage. On the other hand, the interactions of the bubble with the shock waves are not covered, especially for the case of the explosion near different boundaries [37]. In order to accurately understand the consecutive phenomena and pressure field of underwater explosion, it is necessary to simulate the whole process of underwater explosion continuously. Therefore, this paper aims to conduct the continuous simulation of a complete process of underwater explosion which is initiated from a charge detonation. This is different from the extensive studies of bubble dynamics that are directly initiated from high-pressure bubble. Due to the violent and large deformation, the high compressibility and the moving interface, the numerical simulation of underwater explosion is challenging for traditional Lagrange formulations. In this paper, the Eulerian finite element approach [38,39] based on ABAQUS user subroutines is applied to study the consecutive process of underwater explosion, including the shock wave and the bubble motion. The Eulerian finite element method is attractive for simulating underwater explosion process because it can overcome the mesh excessive deformation and distortion. The governing equations are solved by operator splitting for computational efficiency. Also, volume of fraction (VOF) method is combined to track the material interface. This paper is presented as follows, in Section 2, the governing equations are described. The operator splitting method used to solve mass, momentum and energy conservation in the Eulerian finite element formulation is discussed. In Section 3, the whole process of free-field underwater explosion is simulated, including shock propagation and bubble pulsation. Then the numerical results are compared with empirical formulas and experimental results [5] to validate the reliability. In Section 4, the results from the continuous simulation are compared with the simulation of underwater explosion bubble initialized from a high-pressure bubble. In Section 5, the underwater explosion near wall are focused and finally some conclusions from the continuous simulation are enclosed.

Dρ = −ρ∇⋅v, Dt ρ

(2)

∂ρe ∂t

(3)

+ ∇⋅(ρev ) = −p∇⋅v,

(6)

Eqs. (4)-(6) can be solved by the explicit Finite Element Method [43], in which the governing equations are discretized based on the Galerkin method, e.g. the semi-discretized form of the momentum equation:

∫Ω ρΦβ Φα dΩ DDvtα = ∫Ω ρΦα f dΩ + ∫Ω p∇Φα dΩ − ∫Γ pΦα dΓ,

(7)

where Φ is the shape function of element, α and β is the number of nodes, Ω is the computational domain, Γ is the boundary of the computational domain. For the discretization field, eight-node hexahedral element with one-point integration and hourglass control [44] are used to solve Eqs. (4)–(6) so as to enhance calculation efficiency and avoid volumetric locking. The explicit formulation advances the solution in time with the central difference method which provides the second-order accuracy. Through the discretized equation Eq. (7), the node velocity will be obtained and thus the node position can be updated. Similarly, the density and energy of each element are updated by Eq. (4) and Eq. (6). In the second phase, which is called the transport or advection phase, the mesh moves from its Lagrangian position calculated in the last phase to its initial spatial position. Meanwhile, state variables of each element are transported from Lagrangian mesh to initial Eulerian mesh (the transport method discussed in the following), namely, solving the advection equation, i.e. Eq. (8) where ϕ is the solution variable [42]. By integrating Eq. (8) in the flow field, due to the time independent of Ω, so Eq. (9) holds, which reveals that the element variables are conserved during the transport phase, that is, the integration of the field variables in flow field are invariant in the transport calculation.

During underwater explosion, the pressure of shock wave is of the order of GPa and the bubble jet velocity may exceed 300 m/s [7,40]. For underwater explosion with high Reynolds number, the viscosity of the fluid is negligible. In this paper, the flow field of underwater explosion is described by inviscid compressible Eulerian equations which are composed of mass, momentum and energy conversation equations, they read:

∂ρv + ∇⋅(ρv ⊗ v ) = −∇p + ρf , ∂t

(5)

Dt

2.1. Eulerian finite element formulation

(1)

Dv = −∇p + ρf , Dt

ρ De = −p∇⋅v.

2. Theoretical background

∂ρ + ∇⋅(ρv ) = 0, ∂t

(4)

∂ϕ + v⋅∇ϕ = 0 ∂t

∫Ω ⎛ ∂∂ϕt ⎝

D + v⋅∇ϕ ⎞ dΩ = Dt ⎠

(8)

∫Ω ϕdΩ = 0

(9)

Before the transport of field variables, the material transport between adjacent elements should be calculated. For the problems of underwater explosion, there may be more than one material in one element. Hence the Eulerian formulation for multi-material is adopted based on VOF method in this paper. The volume of material transported between adjacent elements is computed by solving the VOF equation [22,42].

where ρ is the material density, v is the material velocity, p is the pressure, f is the body force, e is the internal energy per unit mass. As indicated in the literature of Cole [1], the temperature at the bubble center after charge detonation can be up to 3000 K. However, the process of underwater explosion is instantaneous (about several

∂Vf ∂t 126

+ v⋅∇Vf = 0

(10)

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

where Vf is the volume fraction of a certain material. In the formulations of solving VOF equation, the material interface in mixed elements, i.e. containing more than one material, is captured firstly by interface tracking algorithms. In the paper, Simple Line Interface Calculation technique is employed [45] for the tracking of material interface. The distribution of materials in all the elements is described according to the redistribution of volume fraction of materials, and the material interface in mixed element is approximated by a plane [46,47]. Then based on the certain material interface, the transports of each material between adjacent elements are calculated. Subsequently, Eq. (8) is solved for the mass, momentum and energy respectively with the initial condition computed in the last phase. That is, the field variables computed in Phase 1 are transported from the deformed mesh onto the initial mesh. A good transport algorithm is accurate, stable, conservative and monotonic. In the transport presented in this paper, the transport algorithm of Monotone Upwind Schemes Conservation Laws (MUSCL) [48–50] is used. For one dimension case, this algorithm integrates Eq. (8) in one space-time mesh for solving the transport of the field variables between the step n and n + 1.

∫t ∫x ⎛ ∂∂ϕt ⎝

+v

∂ϕ ⎞ dx dt = ∂x ⎠

∫x (ϕ |tt

n + 1 )dx n

+

∫t v (ϕ | xx

n + 1 )dt n

= 0.

Table 1 Material parameters and coefficients in JWL EOS for TNT [54].

⎜

R1

R2

ω

e (J/kg)

1630

3.71 × 1011

3.21 × 109

4.15

0.95

0.3

4.29 × 106

(

ρ0 c02 μ ⎡1 + 1 − ⎣ ⎡1 − (s1 − 1) μ − ⎣

γ0 2

) μ−

s 2 μ2 μ+1

−

a 2 μ⎤ 2 ⎦ s3 μ3

2

+ ρ0 (γ0 + aμ) Δe

⎤

(μ + 1)2 ⎦

(13)

In the expansion state, the pressure of water is

p = ρ0 c02 μ + ρ0 (γ0 + aμ) Δe

(14)

where ρ0 is the initial density. μ is defined as μ = ρ / ρ0 − 1. When μ > 0 , the water is in the compressed state, and when μ < 0 , water is in the expanded state. c0 is the sound speed at ambient pressure. Δe is the internal energy increment per unit mass. γ0 is the Gruneisen coefficient. a is the volume correction coefficient. s1, s2 and s3 are the fitting coefficients obtained from Hugoniot impact experiment. The values of the corresponding coefficients are listed in Table 2. For the cavitation in water, the Mie-Gruneisen equation of state should be combined with a cavitation model to describe the cavitation behaviors. Because the cavitation dynamics is not the focus in the present study, only a simple cut-off model [57] is applied to describe the cavitation of water, i.e. Eq. (15).

(11)

p , p > Pv p=⎧ P, ⎨ ⎩ v p < Pv

(15)

where Pv is saturated vapor pressure which is 2300 Pa in this paper. The air is modeled as an ideal gas which satisfies the gamma law equation of state, that is,

p = (γ − 1) ρe

(16)

where γ is the ratio of the specific, and γ = 1.4 in this paper. The initial density of air is 1.205 kg/m3 , and therefore the initial energy of air is 2.10 × 105J/kg , which corresponds to the atmospheric pressure of 101,325 Pa according to the equation of state.

In the governing equations mentioned above, the number of variables to be solved is one more than the number of equations. So in order to enclose the set of equations, the equation of state (EOS) is supplemented. The equations of state for the materials applied in underwater explosion are as follows. During the underwater explosion, shock wave propagation and bubble pulsation are simulated continuously in the present work. Hence, the charge detonation process should be considered, and the standard Jones-Wilkins-Lee(JWL) [53] equation of state is employed for detonation products in this paper. For JWL equation of state, the chemical reaction of explosive is ignored, and it is assumed that underwater explosion begins directly from detonation products. ⎟

B (Pa)

p=

2.2. Equation of state

⎜

A (Pa)

high-speed impact, and the compressibility of the fluid could be considered by the Mie-Gruneisen equation of state. Many studies have studied the underwater explosion problems with Mie-Gruneisen equation of state [55,56]. In this paper, the Mie-Gruneisen equation of state in cubic Hugoniot form [54] is used for water. The pressure of water in the compression state is given by

For a certain element, the field variables are approximated by a linear function, so in Eq. (11), they can be integrated with the secondorder accuracy. The MUSCL transport algorithm is just available to element-centered variables, e.g. the mass and energy. In this paper, the velocity is not centered at the element, but at the node. Therefore, the momentum transport is performed with the Half Index Shift algorithm [51], in which an auxiliary set of element-centered variables is constructed from the momentum and transported subsequently, then the new node-centered velocities are reconstructed from the auxiliary variables. During a single time step of transport, the mesh deformation should be less than the dimensions of elements, that is, the Courant number is less than one in computational fluid dynamics [39]. In the transport phase, the time step Δt is fictitious, namely the time step is not updated. But the condition of Courant number limits the time step in Lagrangian phase. For the three dimensional problems in this paper, the transport calculation can be divided into a sequence of one-dimensional transport steps by operator splitting method [52].

ωη ⎞ − R2 ωη ⎞ − R1 e η + ωηρ0 e e η + B ⎛1 − p = A ⎛1 − R2 ⎠ R 1⎠ ⎝ ⎝

ρ0 (kg/m3)

3. Results and discussions 3.1. Continuous simulation of free-field underwater explosion In this section, based on the Eulerian finite element formulation, the whole process of free-field underwater explosion is simulated continuously, including shock wave propagation and bubble oscillation. The simulation is implemented using 3-dimensional model since the problem of bubble collapse is completely 3-dimensional. The free-field underwater explosion experiment by Klaseboer [5] is representative and adopted by many researchers, so following the Klaseboer’s experiment [5], the initial configurations for free-field underwater explosion in this paper is established for comparison and verification,

⎟

(12)

where η is the ratio of the density of detonation products and the initial density of explosive charge, η = ρ / ρ0 . e is the internal energy per unit mass. A , B , R1, R2 and ω are the coefficients obtained by fitting experimental data. Some material parameters and coefficients in the equation of state for TNT are shown in Table 1. Mie-Gruneisen equation of state usually applies to the problems of

Table 2 Material parameters and coefficients in Mie-Gruneisen EOS for water [54].

127

ρ0 (kg/m3)

c0 (m/s)

γ0

a

s1

s2

s3

998.16

1483

0.4934

0

2.56

−1.986

0.2268

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

Fig. 1. The initial configuration for free-field underwater explosion. It is consistent with the Klaseboer’s experiment [5]. The view is the vertical profile of spherical computational domain. Left: the arrangement of the computational case. Right: the layout of measuring points, where r =0.113, 0.121, 0.145, 0.171, 0.205, 0.237, 0.271, 0.289 m.

Fig. 2. The initial field of hydrostatic pressure.

Fig. 3. The gird generation strategy of a quarter of comuptational domain. (The solid line is the domain boundary line, the thick solid line is the block line, and the dotted line is the grid line.).

which is shown in Fig. 1. The charge Hexocirex is replaced by the standard TNT in the simulation with an equivalent coefficient 1.23. The charge mass in the experiment [5] covers three cases: 10 g, 35 g and 55 g, so it is transformed to the TNT mass: 12.3 g, 43.05 g and 67.65 g in the simulation respectively. Besides, some measuring points are set for analyzing requirement, including 8 measuring points for shock wave pressure and 3 measuring points (PⅠ, PⅡ and PⅢ) for the pressure of bubble oscillation, as displayed in Fig. 1. Note that the depth of the explosive charge below free surface is 3.5 m in this paper which is same as the experimental layout [5]. In order to maintain consistency, the depth of the explosive charge is 3.5 m all the time in the subsequent cases, namely, the effect of the free surface is consistent. Besides, the maximum bubble radius is about 0.5 m in the present work. The depth is about 7 times of the maximum bubble radius, so the effect of the free surface can be ignored. Also, in the experiment, the sidewalls of the pool is treated to effectively diffract the shock waves induced by the explosion [5], hence in the simulation, the non-reflecting boundary condition based on impedance match theory is used for the boundary of the computational domain. This boundary condition should be placed far away from the

Fig. 4. The time histories of bubble radius for different gird sizes.

128

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

Fig. 5. The shock wave propagation during the underwater explosion of the charge mass of 35 g. This is the vertical mid-profile view. Left: pressure field. Right: velocity field. The red line denotes the bubble interface.

computational domain is such large that the dynamic pressure pd and the material velocity v at the boundary can be assumed to be linear, namely pd =ρc(v·e), where e is the unit vector of shock wave propagation direction. Hence, the total pressure exerting on the external boundary is equal to the sum of pd and the hydrostatic pressure. The effect of gravity is also taken account and the gravity acceleration is 9.8 m/s2. As shown in Fig. 2, the initial field of hydrostatic pressure is displayed. Based on OPENMPI commands, the computation is implemented with 60 threads on DELL workstation with Intel Xeon E74890 clocked at 2.8 GHz and RAM of 256 Gb. Firstly, the spatial convergence is discussed with different grid refinements. The gird generation strategy used in this paper is shown in Fig. 3. The size of the gird is uniform in the central spherical block and varies radially in the outer blocks. The number of girds along the radius of the central sphere is noted as N . With a magnification of 1.2, there are four cases of N, namely, N = 34, 42, 50 and 60. Taking the case of the charge mass of 35 g as an example, the time histories of the bubble radius are plotted in Fig. 4. Actually, the bubble shape is not completely symmetric because of the gravity effect, but it is assumed that the bubble is spherical during the first bubble oscillation. The bubble radius is calculated from the bubble volume based on this spherical assumption, as shown in Eq. (17), where the bubble volume is obtained by integrating the volume fraction of the material in the VOF method. Obviously, the histories of the bubble radius of N =50 and N =60 are almost coincide in Fig. 4. Considering gird convergence rate [58], the rate of the first three curves and the latter three curves are 1.3 and 4.2 respectively. The numerical results are convergent obviously and the subsequent calculations are based on the gird size of N =60. In this case with about 4 million elements, the computational cost is about 24 h.

Fig. 6. The time histories of the shock wave pressure.

1/3

3V (t ) R (t ) = ⎛ bubble ⎞ 4π ⎝ ⎠

(17)

Fig. 5 presents the process of shock wave propagation for the charge mass of 35 g. Within microseconds after the charge detonated, the explosive charge becomes the high-pressure and high-temperature detonation products. Subsequently, the detonation products begin to expand and compress the surrounding water, resulting in the propagation of spherical shock wave. In Fig. 5(a), it can be seen that the speed of shock wave propagation is larger than that of bubble expansion. In Fig. 5(b) when shock wave propagates to the free surface, local cavitation occurs due to the superposition of reflection waves and incident waves. A set of measuring points is arranged in the depth of the explosive charge, as shown in Fig. 1, and the shock wave pressure is measured.

Fig. 7. The peak pressure of the shock wave at different distances from the explosive charge.

domain of interest. In the cases studied in this paper, the maximum radius of the bubble at the center of the computational domain is about 0.5 m and the radius of computational domain is 7.5 m. Apparently, the 129

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

Fig. 8. Underwater explosion bubble motion of the 35 g case. This is the vertical mid-profile view. Left: the pressure field. Right: the velocity field. The red line denotes the bubble interface.

The empirical formula for the pressure of underwater explosion shock wave proposed by Zamyshlyayev [8] is as follows.

P (r , t ) =

Pm e−t / θ , 0≤t≤θ ⎧ , θ 1.5 ⎨ 0.368Pm [1 − (t /tp) ], θ < t ≤ tp t ⎩

Pm is the peak pressure, θ is the coefficient of attenuation, tp is the elapsed time between the charge detonation and the maximum volume of the bubble. Pm and θ are calculated respectively from the following formulas [8].

(18)

where is the distance of the measuring point from the explosive charge, 130

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

Fig. 11. The comparison of the bubble motion and flow field pressure at PⅠ between UNDEX bubble and HP bubble.

Fig. 9. The comparisons of bubble radius between the calculations and the experiments [5] for the 10 g, 35 g and 55 g cases. Table 3 The maximum bubble radius Rmax and the first period of oscillation Tp for the 10 g, 35 g and 55 g cases [5].

Rmax (m)

Tp (ms)

Case

Experimental result [5]

Numerical result

Error

10 g 35 g 55 g 10 g 35 g 55 g

0.34 0.498 0.574 53 82 94

0.344 0.518 0.596 53.9 81.8 93.9

1.2% 4.0% 3.8% 1.7% −0.2% −0.1%

Fig. 12. The bubble jet of HP bubble. This is the vertical mid-profile view. Left: the pressure field. Right: the velocity field. The red line denotes the bubble interface.

points between the numerical results and the empirical formula are shown in Fig. 6. The numerical results agree well with that from the empirical formula. Besides, Fig. 7 presents the peak pressure of shock wave at different distances from the explosive charge. It can be seen that the peak pressure declines with the increasing of the distance away from the explosive charge, and the tendencies of the datum from numerical simulation and empirical formula agree well with each other. In Fig. 8, the bubble motion of the underwater explosion with a charge mass of 35 g is shown. Accompanied with the formation of shock wave, the detonation products begin to expand, thus the internal pressure of the bubble decreases. In Fig. 8(a), the bubble reaches the maximum volume at 40 ms, and the internal pressure of the bubble is lower than the ambient pressure. So the bubble begins to shrink continuously, as shown in Fig. 8(b) and (c). In Fig. 8(d), due to the difference pressure at the upper and lower surface of the bubble, an upward jet is formed and the jet velocity exceeds 200 m/s. For the three cases of different charge masses, namely 10 g, 35 g and 55 g, the comparisons of bubble radius during the first oscillation

Fig. 10. The comparisons of the oscillation pressure between the calculations and the experiments [5] for the 35 g case. 1.5

( ) ( )

⎧ 4.41 × 107 ⎪ Pm = ⎨ 7 ⎪5.24 × 10 ⎩

θ=

⎧ 0.45r0

W1/3 r

r 0.45 r0

( )

,

1.13 W1/3 , r

6 ≤ r r0 < 12 , 12 ≤ r r0 < 240

× 10−3, r r0 < 30

⎨ 3.5 r0 lg(r r ) − 0.9 , r r ≥ 30 0 0 c ⎩

(19)

, (20)

where W is the mass of explosive charge, r0 is the radius of the explosive charge. The comparisons of the shock wave pressure at two measuring 131

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

3.2. Simulation of underwater explosion initialized from a high-pressure bubble Traditionally, underwater explosion is always divided into two stages and studied separately. In the study of bubble dynamics, the static high-pressure bubble is usually used as the initial condition of underwater explosion bubble. The initial conditions are widely applied in BEM. In this section, the traditional initial condition is discussed and compared with the results from Eulerian finite element method. Based on many experiments of underwater explosion, researchers summarized the empirical formulas of maximum bubble radius Rmax and internal initial pressure P0 for a spherical TNT charge [1].

Rmax = 3.38 P0 = 1.39 ×

1/3 W H + 10 W γ 105 V , 0

(

)

( )

(21)

where W is the mass of explosive charge, H is the water depth of explosive charge, V0 is the initial volume of the bubble, and γ =1.25 for TNT. In the incompressible, inviscid and irrotational flow field, the equation of motion for spherically symmetric bubble is Eq. (22), i.e. the so-called Rayleigh-Plesset equation [59,60].

3 ρRR¨ + ρR˙ 2 = PB − PH , 2

(22)

where ρ is the density of water, PB is the internal pressure of the bubble, PH is the hydrostatic pressure at the depth of H , overdot means differentiation with respect to time. From Eqs.(21)-(22), the relationship between the bubble initial radius R 0 , Rmax , W and H is obtained [27]. 3

γ

3(1 − γ )

5 ⎛ R 0 ⎞ − 1 = 1.39 × 10 ⎛⎜ 3W ⎞⎟ ⎡1 − ⎛ R 0 ⎞ 3 ⎢ ( 1) R P γ − max H ⎝ ⎠ ⎝ Rmax ⎠ ⎝ 4πR max ⎠ ⎣ ⎜

⎟

⎜

⎟

⎤ ⎥. ⎦

(23)

R 0 is carried out by a root solving procedure, and P0 can be calculated with Eq. (21) subsequently. Taking the free-field underwater explosion case of the charge mass of 35 g as an example, the initial bubble radius R 0 is 41.59 mm and the internal initial pressure P0 is 68.63 MPa. Additionally, the initial density of the bubble is also calculated so that the mass of the bubble is equal to the solid explosive charge, thus the initial density is 141.8 kg/m3 , and the ideal gas equation of state (as mentioned in Section 2.2) is adopted. The static high-pressure bubble (referred to as HP bubble hereafter) with such initial conditions is compared to the underwater explosion bubble (hereafter referred to as UNDEX bubble). The bubble motion and the flow field pressure at PⅠ between UNDEX bubble and HP bubble are displayed in Fig. 11. The maximum radius and the oscillation period of HP bubble are 0.464 m and 76.8 ms respectively, which deviate from the results of UNDEX bubble. In addition to the differences in bubble motion, the oscillation pressure of HP bubble at the minimum volume is also smaller. Fig. 8(d) and Fig. 12 present the moment when the two bubbles collapse to form a jet. As shown the jet shapes are different, and the jet velocities are 209 m/s and 105 m/s, respectively. The HP bubble is relatively smaller and weaker than UNDEX bubble. That is to say, the initial conditions of the high-pressure bubble in underwater explosion are not compatible with the Eulerian finite element method applied in this paper. In the traditional two-stage study, the shock wave stage is ignored in the initial condition of bubble in underwater explosion. However, there is still a peak pressure similar to the shock wave in the pressure curve of the flow field (see Fig. 11, the first peak of the magenta solid line), which is much smaller than the real pressure peak of the shock wave. On the other hand, the traditional initial condition of bubble also ignores the initial velocities of the bubble and the surrounding flow field, and the fluid is incompressible in Rayleigh-Plesset equation. It is obvious that there are some differences between the traditional two-stage research and the real underwater explosion. In addition, it is difficult to

Fig. 13. The initial condition for underwater explosion near wall. The configuration is a continuation of the free-field explosion which is referred to [5].

period between the simulations and the experiments [5] are presented in Fig. 9. As revealed, the numerical and experimental results [5] are in good agreement for all cases. Table 3 lists the maximum bubble radius and the first period of oscillation for the three cases. For the maximum radius of the bubble, the error between calculations and experiments [5] is about 1.2% for the 10 g case, 4.0% for the 35 g case, and 3.8% for the 55 g case, but for the first period of oscillation, the error is 1.7% for the 10 g case, -0.2% for the 35 g case, and -0.1% for the 55 g case. Overall, the numerical results have well represented the bubble oscillation of underwater explosion. When the bubble collapses and reaches the minimum volume, the bubble begins to rebound. Meanwhile, the pressure wave is generated and propagates outward, as shown in Fig. 8(e). Three pressure sensors (PⅠ, PⅡ and PⅢ) are placed as those in the experiment of Klaseboer [5] for measuring the pressure wave resulting from bubble rebound. The time histories of the pressure are plotted in Fig. 10. Overall, they are in good agreements for all the cases except a slight underestimation of the pressure peak for PⅠ and PⅡ. In summary, the continuous simulation of underwater explosion including the shock wave phase and the bubble phase is carried out by the Eulerian finite element approach, which has overcome the shortcomings of the traditional two-stage approach of underwater explosion problem. Considering shock wave propagation and bubble motion characteristics, the entire process of underwater explosion is well simulated, which verified the feasibility and reliability of the Eulerian finite element method in dealing with underwater explosion problem. 132

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

Fig. 14. The initial stage of the underwater explosion near-wall. This is the horizontal mid-profile view. The red line denotes the bubble interface. Left: the pressure field. Right: the velocity component field. Table 4 The maximum bubble radius Rmax and the first period of oscillation Tp of UNDEX bubble and HP bubble for different cases of r .

Rmax (m)

Tp (ms)

r Rmax _ FF

UNDEX bubble

HP bubble

Difference

0.965 0.772 0.579 0.965 0.772 0.579

0.505 0.503 0.502 96.5 98.7 99.9

0.456 0.454 0.454 89.1 91.9 94.0

−9.7% −9.7% −9.6% −7.8% −6.9% −5.9%

Fig. 16. The relative deviation of the maximum bubble radius and the first period of oscillation for different r .

determine a reasonable and correct initial condition of bubble in underwater explosion. Therefore, it is necessary to simulate the whole process of underwater explosion continuously. The Eulerian finite element method presented in this paper, which can deal with the whole process of underwater explosion, is a good choice for simulating underwater explosion. 3.3. Continuous simulation of near-wall underwater explosion Fig. 15. The maximum bubble radius and the first period of oscillation for different r .

When a charge is detonated near wall during underwater explosion, the shock wave will be reflected and the reflected wave will interact with the bubble. The stages of shock wave propagation and the bubble 133

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

bubble has changed due to the interaction of the shock wave and the bubble. For the maximum bubble radius, the differences between UNDEX bubble and HP bubble in free field and near the wall are about 10.3% and 9.7%, respectively. The changes of the difference are slight for the three cases. But for the first period of oscillation, the change of differences between the two bubbles is obvious. As r increases, the difference change gradually becomes larger. When r = 0.5 m, the difference change becomes from 6.0% to 7.8%. In other words, compared with the traditional method, the influence of the continuous simulation of the whole process on the bubble motion is about 2% for the near-wall underwater explosion. In the near-wall underwater explosion, shock wave interacts intensely with bubble motions because of the reflection of the wall. In order to further explore the influence of the wall on bubble motions in underwater explosion, the bubble motion with different distances r is investigated. The underwater explosion of the charge of 35 g near wall is also taken as an example. The distance r increases gradually from 0.5 m by the interval of 0.5 m. Fig. 15 displays the maximum bubble radius and the first period of oscillation of UNDEX bubble for different cases of r . As r increases, Rmax and Tp approach the maximum radius and the first period of the bubble oscillation in the free field (hereafter denoted as Rmax _ FF and Tp _ FF , respectively). The relative deviations of Rmax and Tp for all the cases are calculated respectively, as shown in Fig. 16. The maximum deviation of on Tp can reach more than 20%, but the deviation on Rmax is relatively small, within 3%. Because the maximum bubble volume is mainly determined by the bubble energy, which is determined by the initial internal energy of the explosive, so the influence on Rmax is slight. The slight influence possibly results from the calculation of bubble radius based on the spherical bubble assumption. Moreover, the influence of the wall on the bubble motion decreases rapidly as r increases. When r = 4Rmax _ FF , the influence of the wall on Rmax and Tp is reduced to 2%; when r = 9Rmax _ FF , the influence of the wall decreases to 0.7%. Except the maximum bubble radius and the first period of oscillation, the jet characteristics of the bubble are also observed and compared in Fig. 17, where the jet pattern in different cases at the same moment are plotted. In the near-wall underwater explosion, the bubble is attracted by the wall which results in the deflection of the jet direction. As r increases, the influence of the wall on the bubble motion becomes smaller and smaller. When r = 5m , the jet direction of the bubble is almost vertical upward, which is consistent with that in the free-field underwater explosion. To sum up, the influence of the wall on the first period of bubble oscillation is greater than that on the bubble maximum radius in nearwall underwater explosion, and these influences decrease rapidly as the bubble is away from the wall. With the parameters discussed in this paper, when the distance between the explosive charge and the wall is 9 to 10 times of the maximum bubble radius in the free-field case, the bubble motion is almost consistent with that in the free-field case.

Fig. 17. The jet pattern in different cases at the same moment.

motion are difficult to clearly divide. In addition, in the practical applications, the underwater explosion in fully free-field are very rare, most of them are affected by complex boundaries. Therefore, the underwater explosion near wall is investigated in this section. According to the free-field underwater explosion cases in Section 3.1, the boundary condition is changed and the initial condition for underwater explosion near wall is presented in Fig. 13. The horizontal distance between the explosive charge and the vertical rigid wall is 0.5 m and the explosive charge mass is 43.05 g TNT, while other conditions are same as the free-field case. Fig. 14 shows the initial stage of the near-wall underwater explosion. In Fig. 14(a), at 0.4 ms the reflection wave has not reached the bubble, and the expanding velocity of the bubble surface is about 106 m/s. When the reflected wave hits the bubble, the bubble motion becomes asymmetric. The biggest difference between the expanding velocity of bubble’s front and rear surfaces is about 22 m/s, see in Fig. 14(b), then the difference decreases to 4 m/s at 0.9 ms, as plotted in Fig. 14(d). During this period, the bubble is compressed due to the reflected wave and the subsequent motion of the bubble is affected. Conversely, the bubble motion also affects the propagation of reflected wave. When the reflected wave hits the bubble, the rarefaction wave is formed because of the interface reflection. Because of the superposition of shock waves and reflected waves, the pressure field of the flow field becomes extremely complex. Generally, a strong interaction occurs between the shock wave and the bubble motion in the near-wall underwater explosion. However, in the two-stage approach, the shock wave is not considered exactly with the initial condition of high-pressure bubble, so it is difficult to study this interaction by traditional methods. Therefore, it is important to carry out the continuous simulation of the whole process of underwater explosion. In order to further analyze the importance of the continuous simulation of near-wall underwater explosion, the UNDEX bubble and HP bubble (initial conditions are presented in the previous section) are respectively placed in the same position near the wall with the varying horizontal distances r = 0.3 m, 0.4 m and 0.5 m. Focusing on the comparison of the maximum radius and the period of oscillation between UNDEX bubble and HP bubble, the results are listed in Table 4. Note that actually the bubble motion is non-spherical for r = 0.3 m and r = 0.4 m. However, the bubble radius is still calculated on the basis of spherical assumption for the simplicity of calculation, and the bubble radius can be regarded as a representative variable of bubble volume. For the cases of the underwater explosion in free field and near wall, the difference of motion characteristics between UNDEX bubble and HP

4. Conclusions In the present work, the whole process of underwater explosion is continuously simulated by Eulerian finite element method, that is, the shock wave stage and the bubble motion stage are consecutive. This approach has overcome the shortcomings of the traditional two-stage researches and shows good agreements with the experimental results [5]. Many comparative studies demonstrate that the continuous simulation is more reasonable than the traditional simulation based on the initial conditions of underwater explosion bubble widely used in previous literature. In Eulerian finite element method, high-pressure bubble under traditional treatments is relatively smaller and weaker than underwater explosion bubble, and the errors of bubble motion characteristics even exceed 10%. Therefore, more reasonable and accurate initial conditions for underwater explosion bubble should be studied further. 134

Applied Ocean Research 80 (2018) 125–135

W.T. Liu et al.

Based on the Eulerian finite element method, the underwater explosion process near wall is also continuously simulated. It is found that the reflection wave interacts intensively with the expanding bubble in the vicinity of the wall. For the near-wall underwater explosion, the influence of the continuous simulation on bubble motion is about 2% compared to the simulation of high-pressure bubble with an initial condition (the characteristics of bubble motion in free field is consistent) and the influences on the period of oscillation are relatively larger and more obvious than those on the maximum bubble radius. When the distance between the bubble and the wall is more than 4 times of the maximum bubble radius of free-field case, the influence of the wall on the bubble motion decreases to less than 2%. When it is up to 9 to 10 times, the jet direction is almost vertical upward and it is almost consistent with that in the free-field case.

[22] Z.L. Tian, Y.L. Liu, A.M. Zhang, S.P. Wang, Analysis of breaking and re-closure of a bubble near a free surface based on the Eulerian finite element method, Comput. Fluids 170 (2018) 41–52. [23] E.A. Brujan, G.S. Keen, A. Vogel, J.R. Blake, The final stage of the collapse of a cavitation bubble close to a rigid boundary, Phys. Fluids 14 (2002) 85–92. [24] M. Koch, C. Lechner, F. Reuter, K. Köhler, R. Mettin, W. Lauterborn, Numerical modeling of laser generated cavitation bubbles with the finite volume and volume of fluid method, using OpenFOAM, Comput. Fluids 126 (2016) 71–90. [25] S. Li, Y.B. Li, A.M. Zhang, Numerical analysis of the bubble jet impact on a rigid wall, Appl. Ocean. Res. 50 (2015) 227–236. [26] C. Wang, B.C. Khoo, An indirect boundary element method for three-dimensional explosion bubbles, J. Comput. Phys. 194 (2004) 451–480. [27] E. Klaseboer, B.C. Khoo, K.C. Hung, Dynamics of an oscillating bubble near a floating structure, J. Fluids Struct. 21 (2005) 395–412. [28] Y.L. Liu, Q.X. Wang, S.P. Wang, A.M. Zhang, The motion of a 3D toroidal bubble and its interaction with a free surface near an inclined boundary, Phys. Fluids 28 (2016). [29] K.-K. Kan, J.H. Stuhmiller, P.C. Chan, Simulation of the collapse of an underwater explosion bubble under a circular plate, Shock. Vib. 12 (2005). [30] S. Li, A.M. Zhang, R. Han, Y.Q. Liu, Experimental and numerical study on bubblesphere interaction near a rigid wall, Phys. Fluids 29 (2017) 092102. [31] Z.R. Li, L. Sun, Z. Zong, Numerical analysis of gas bubbles in close proximity to a movable or deformable body, Arch. Appl. Mech. 83 (2013) 1715–1737. [32] A.M. Zhang, W.B. Wu, Y.L. Liu, Q.X. Wang, Nonlinear interaction between underwater explosion bubble and structure based on fully coupled model, Phys. Fluids 29 (2017) 082111. [33] A.M. Zhang, X.L. Yao, J. Li, The interaction of an underwater explosion bubble and an elastic–plastic structure, Appl. Ocean. Res. 30 (2008) 159–171. [34] S.W. Gong, B.C. Khoo, Transient response of stiffened composite submersible hull to underwater explosion bubble, Compos. Struct. 122 (2015) 229–238. [35] A.M. Zhang, S.-p. Wang, C. Huang, B. Wang, Influences of initial and boundary conditions on underwater explosion bubble dynamics, Eur. J. Mech. B/Fluids 42 (2013) 69–91. [36] S. Li, R. Han, A.M. Zhang, Q.X. Wang, Analysis of pressure field generated by a collapsing bubble, Ocean. Eng. 117 (2016) 22–38. [37] B.C. Nie, J.C. Li, H.Q. Zhang, Interaction between reflected shock and bubble in near-wall underwater explosion, Proc. Eng. 126 (2015) 344–348. [38] D.J. Benson, S. Okazawa, Contact in a multi-material Eulerian finite element formulation, Comput. Methods Appl. Mech. Eng. 193 (2004) 4277–4298. [39] D.S. Simulia, A.G. Fallis, ABAQUS documentation. (2013). [40] P. Cui, A.M. Zhang, S. Wang, B.C. Khoo, Ice breaking by a collapsing bubble, J. Fluid Mech. 841 (2018) 287–309. [41] A.J. Chorin, T.J.R. Hughes, M.F. McCracken, J.E. Marsden, Product formulas and numerical algorithms, Commun. Pure Appl. Math. 31 (1978) 205–256. [42] D.J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, Comput. Methods Appl. Mech. Eng. 99 (1992) 235–394. [43] S.R. Wu, L. Gu, Introduction to the Explicit Finite Element Method for Nonlinear Transient Dynamics, John Wiley & Sons, Hoboken, New Jersey, USA, 2012. [44] T. Belytschko, L.P. Bindeman, Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems, Comput. Methods Appl. Mech. Eng. 88 (1991) 311–340. [45] W.F. Noh, P. Woodward, SLIC simple line interface calculation, in: A.I. van de Vooren, P.J. Zandbergen (Eds.), Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics June 28 – July 2, 1976 Twente University, Enschede, Berlin, Heidelberg, Springer Berlin Heidelberg, 1976, pp. 330–340. [46] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152. [47] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [48] B. Van Leer, Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys. 23 (1977) 276–299. [49] T. Buffard, S. Clain, Monoslope and multislope MUSCL methods for unstructured meshes, J. Comput. Phys. 229 (2010) 3745–3776. [50] C. Le Touze, A. Murrone, H. Guillard, Multislope MUSCL method for general unstructured meshes, J. Comput. Phys. 284 (2015) 389–418. [51] D.J. Benson, Momentum advection on a staggered mesh, J. Comput. Phys. 100 (1992) 143–162. [52] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of gas/liquid Multiphase Flows, Cambridge University Press, New York, USA, 2011. [53] B.M. Dobratz, LLNL Explosives Handbook: Properties of Chemical Explosives and Explosives and Explosive Simulants, Lawrence Livermore National Laboratory, CA, United States, 1981. [54] Y.S. Shin, M. Lee, K.Y. Lam, K.S. Yeo, Modeling Mitigation Effects of Watershield on Shock Waves, Shock. Vib. 5 (1998). [55] M. Cheng, K.C. Hung, O.Y. Chong, Numerical study of water mitigation effects on blast wave, Shock. Waves 14 (2005) 217–223. [56] D. Steinberg, Spherical Explosions and the Equation of State of Water, Lawrence Livermore National Lab., CA (USA), 1987. [57] J.E. van Aanhold, G.J. Meijer, P.P.M. Lemmen, Underwater shock response analysis of a floating vessel, Shock. Vib. 5 (1998). [58] S. Marrone, Enhanced SPH Modeling of Free-surface flows With Large Deformations, Sapienza University of Rome, Roma, Italy, 2012. [59] M.S. Plesset, The dynamics of cavitaion bubbles, J. Appl. Mech. 16 (1949) 277–282. [60] L. Rayleigh VIII, On the pressure developed in a liquid during the collapse of a spherical cavity. The London, Edinburgh, and Dublin, Philos. Magaz. J. Sci. 34 (1917) 94–98.

Acknowledgements This paper is supported by the National Natural Science Foundation of China (51609049, 51879053, U1430236), and the Natural Science Foundation of Heilongjiang Province (QC2016061, E2017021). References [1] R.H. Cole, Underwater Explosions, Princeton University Press, Princeton, New Jersey, USA, 1948. [2] A. Mouritz, The effect of underwater explosion shock loading on the flexural properties of GRP laminates, Int. J. Impact Eng. 18 (1996) 129–139. [3] F.R. Ming, P.N. Sun, A.M. Zhang, Investigation on charge parameters of underwater contact explosion based on axisymmetric SPH method, Appl. Math. Mech. 35 (2014) 453–468. [4] G. Chahine, R. Duraiswami, K. Kalumuck, Boundary element method for calculating 2-D and 3-D underwater explosion bubble loading on nearby structures including fluid structure interaction effects, Naval Surface Weapons Center, Dahlgren Division, Weapons Research and Technology Department, Technical Report NSWCDD/TR-93/46, (1996). [5] E. Klaseboer, K.C. Hung, C. Wang, C.W. Wang, B.C. Khoo, P. Boyce, et al., Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure, J. Fluid Mech. 537 (2005) 387–413. [6] G. Barras, M. Souli, N. Aquelet, N. Couty, Numerical simulation of underwater explosions using an ALE method. The pulsating bubble phenomena, Ocean. Eng. 41 (2012) 53–66. [7] F.R. Ming, A.M. Zhang, Y.Z. Xue, S.P. Wang, Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions, Ocean. Eng. 117 (2016) 359–382. [8] B.V. Zamyshlyaev, Y.S. Yakovlev, Dynamic Loads in Underwater Explosion, National Technical Information Service, Washington, USA, 1973. [9] M.B. Liu, G.R. Liu, K.Y. Lam, Z. Zong, Smoothed particle hydrodynamics for numerical simulation of underwater explosion, Comput. Mech. 30 (2003) 106–118. [10] A.M. Zhang, F.R. Ming, S.P. Wang, Coupled SPHS–BEM method for transient fluid–structure interaction and applications in underwater impacts, Appl. Ocean. Res. 43 (2013) 223–233. [11] F.R. Ming, A.M. Zhang, H. Cheng, P.N. Sun, Numerical simulation of a damaged ship cabin flooding in transversal waves with smoothed Particle Hydrodynamics method, Ocean. Eng. 165 (2018) 336–352. [12] W.F. Xie, T.G. Liu, B.C. Khoo, The simulation of cavitating flows induced by underwater shock and free surface interaction, Appl. Numer. Math. 57 (2007) 734–745. [13] C.F. Hung, J.J. Hwangfu, Experimental study of the behaviour of mini-charge underwater explosion bubbles near different boundaries, J. Fluid Mech. 651 (2010) 55–80. [14] W.S. Hall, Boundary Element Method. The Boundary Element Method, Springer Netherlands, Dordrecht, 1994, pp. 61–83. [15] T. Li, S. Wang, S. Li, A.M. Zhang, Numerical investigation of an underwater explosion bubble based on FVM and VOF, Appl. Ocean. Res. 74 (2018) 49–58. [16] S. Ge, C. Zu-yu, L. Yuan, Z. Ming-shou, W. Jian-yu, Experimental and numerical investigation of the centrifugal model for underwater explosion shock wave and bubble pulsation, Ocean. Eng. 142 (2017) 523–531. [17] W. Xiao, A.M. Zhang, S.P. Wang, Investigation of bubble dynamics of underwater explosion based on improved compressible numerical model, Appl. Ocean. Res. 59 (2016) 472–482. [18] A. Pearson, E. Cox, J.R. Blake, S.R. Otto, Bubble interactions near a free surface, Eng. Anal. Bound. Elem. 28 (2004) 295–313. [19] P. Koukouvinis, M. Gavaises, O. Supponen, M. Farhat, Simulation of bubble expansion and collapse in the vicinity of a free surface, Phys. Fluids 28 (2016) 052103. [20] S.T. Miller, H. Jasak, D.A. Boger, E.G. Paterson, A. Nedungadi, A pressure-based, compressible, two-phase flow finite volume method for underwater explosions, Comput. Fluids 87 (2013) 132–143. [21] A.M. Zhang, Y.L. Liu, Improved three-dimensional bubble dynamics model based on boundary element method, J. Comput. Phys. 294 (2015) 208–223.

135