# Study on Load Characteristics of Underwater Explosion Using RKDG-LS-DGF and BEM

Authors: Longkan Wang; Hailong Chen; Xi Ye; Xiongliang Yao
Journal: Shock and Vibration
 (2015)

---

## Abstract

This study numerically investigates the underwater explosion bubble dynamics during the process from detonation to bubble jet with a hybrid algorithm based on Runge Kutta Discontinuous Galerkin-Level Set-Direct Ghost Fluid method (RKDG-LS-DGF) and boundary element method (BEM). RKDG-LS-DGF method is used to simulate the physical process from the detonation of a column charge at one end to the formation of a nonspherical initial bubble and the process of bubble jet. And BEM is adopted to simulate bubble pulsing characteristics. In addition, the numerical results are compared with the experimental results to verify the feasibility of the numerical method. It is found that, during the detonation process of a column charge, the detonation product experiences a shape change from an initial ellipsoid into a sphere during expansion. After the detonation, the bubble experiences expansion and contraction and develops a jet. The jet threads through the bubble in the opposite direction to gravity and induces a high-pressure region. Subsequently, the pressure of this region decreases when the bubble reexpands after being penetrated by the jet. The numerical results agree well with the experimental data, which proves that axisymmetric RKDG-LS-DGF method and BEM are successfully combined to simulate the whole process of underwater explosion.


---

## Body

## 1. Introduction

Various loads can be generated by underwater explosion, including strong discontinuous shock wave, high-pressure gas bubble impulse, and high-speed jet during the bubble collapse, which may cause severe damage to the structures of naval ships [1]. Theory analysis and experiments are usually combined to study the damage characteristics of underwater explosion bubble. However, the theory analysis is limited to studying some simple and regular structures [2–4], and the disadvantages of experiments are nonrepeatability and a lack of flexibility [5–7]. With the rapid development of bubble dynamics and computational technology, the numerical simulation of underwater explosion bubble has been investigated recently. Rayleigh [8] established the motion equation of a spherical bubble in the incompressible flow field. The growth and collapse of a single vapor bubble near free surface were investigated numerically using BEM by Blake and Gibson [9], and the validity of the numerical results was verified by the comparison with experimental results. In Wang et al.’s research, a Vortex-Ring model [10] was successfully established to solve the numerical problems when the bubble transformed from a simply connected domain to doubly connected domain. Qi et al. [11] applied BEM in the simulation of dynamic characteristic of two neighboring cavity bubbles. Li et al. [12] proposed that the motions of a gas bubble in proximity to a free surface with and without buoyancy force, as well as in shallow water, can be simulated based on a numerical time integration coupled with three-dimensional boundary integral spatial solution. Zhang et al. [13–16] investigated the bubble dynamics in free field and both vicinity to free surface and rigid boundary, combining the numerical simulation and the experimental results.However, some numerical technologies, such as Vortex-Ring model, Vortex sheet model, and numerical fairing, needed to be introduced as modification when the traditional boundary element method (BEM) is used to simulate the bubble jet, which will increase the difficulty and complexity of calculation. Besides, such numerical treatment as smoothing may have negative influence on the accuracy of numerical results. The Level Set method of the Runge Kutta Discontinuous Galerkin-Level Set-Direct Ghost Fluid method (RKDG-LS-DGF), which is an implicit interface tracking method, has an advantage to solve the problem of complex interfaces and the cases when their topology changes. In addition, the simulation of bubble jet can proceed continuously without stop and the introduction of other models. Ever since Reed and Hill [17] carried out the Discontinuous Galerkin (DG) method to solve a neutron transport problem, the RKDG-LS-DGF method has been improved and adopted to solve the problems of compressible fluids [18–20]. Qiu et al. [21–23] used the RKDG finite element method in the simulation of 2D compressible multiphase flow and cavitation problems. With this method, Liu captured the phenomenon of air-water interface motion and cavitation. A coupled solution approach was presented for numerically simulating a near-field early-time underwater explosion (UNDEX) by Park [24] and several test cases were examined and assessed by comparing the RKDG-DGF results with analytic solutions and experimental data. Recently, however, most researchers have focused on the improvement of algorithm for the RKDG-LS-DGF method and the early-time simulation of underwater explosion. There are barely any researches about the simulations of bubble pulsing and bubble jet of underwater explosion using RKDG-LS-DGF method.In this paper, RKDG-LS-DGF and BEM are combined to fully utilize their advantages to simulate the whole process of underwater explosion, including shock wave propagation, bubble pulsing, and jet formation. Firstly, in consideration of the advantages of Level Set method to solve complex interaction, RKDG-LS-DGF numerical model of column charge is established to simulate the physical process from the detonation of a column charge at one end to the formation of a nonspherical initial bubble. Subsequently, in order to improve the computational efficiency and save computation time, BEM is adopted to simulate bubble pulsing, taking the previous results calculated with RKDG-LS-DGF as initial condition. Finally, to avoid manual interference caused by using BEM, such as mesh cutting and numerical smoothing, the bubble jet is simulated by an axisymmetric RKDG-LS-DGF method and the previous results calculated with BEM are treated as initial condition. The results from numerical simulation are compared with the experimental results to verify the feasibility of the hybrid algorithm and the numerical model.


## 2. Theoretical Background and Numerical Method

### 2.1. Basic Equations

Due to the high velocity on the bubble surface, the Mach number in near flow field would be very large during the early stage of expansion and bubble jet. Besides, pressure wave would be generated by the bubble. In this case, detonation products and the surrounding fluid field should be considered as inviscid, irrotational, and compressible fluid. Fluid dynamics equations for 2D compressible axisymmetric flow can be expressed as [24, 25](1)Ut+∇·FU=SsymU,where(2)U=ρρuρvρe,FU=ρuρvρu2+pρvuρuvρv2+pρue+puρve+pv,Ssym=-1xρuρu2ρuvuρe+p,where t, ρ, e, and p denote time, density, energy, and pressure, respectively; x and y indicate the direction along axis; u and v are velocity components in x direction and y direction, respectively. The state equations of different materials are listed as follows.(1) The state equation of ideal gas is given by [24, 25](3)ρe=pγ-1,where γ is the specific heat capacity of gas taken as 1.4 in the simulation of underwater explosion.(2) The Jones-Wilins-Lee (JWL) state equation is used for explosive gas [24, 25]:(4)P=A1-ωρR1ρ0e-R1ρ/ρ0+B1-ωρR2ρ0e-R2ρ/ρ0+ωρe,where A=3.712×1011Pa, B=3.231×109Pa, R1=4.15, R2=0.95, and ω=0.3; ρ0 is the initial density which is equal to 1630kg/m3; E0 is the initial internal energy which is equal to 4.29×106J/m3.(3) The Tait state equation of water is given by [24, 25](5)ρe=p+NB¯N-1,where B¯=B-A, N=7.15, A=1.0×105Pa, B=3.31×108Pa, and ρ0=1000.0kg/m3.


### 2.2. RKDG-LS-DGF Method

Column charge detonation and bubble jet are simulated using axisymmetric RKDG-LS-DGF method and three numerical methods are involved concretely. Discontinuous Galerkin method is applied to solve fluid motion equations [22, 24, 26]. Level Set method has an advantage in capturing moving interfaces of multimaterial flow [27]. These two methods are well combined by Direct Ghost Fluid method [24].The computational domain of fluid is denoted asΩ. Fluid motion equation, that is, (1), is solved by Discontinuous Galerkin method in this paper. After multiplying (1) by the test function w(x,y) on both sides, integrating over a cell Ωj, and integrating by parts, the equations can be derived as(6)ddt∫KUwx,ydΩ+∑e∈∂K∫eFU·n→wx,ydΓ-∫KFU∇wx,ydΩ=0,where n→ is the external normal vector.Line integral and surface integral are obtained by Gaussian quadrature formulas [26](7)∫eFU·n→wx,ydΓ≈∑l=1LωlFUxel,yel,t·n→wxel,yele,∫KFU·∇wx,ydΩ≈∑j=1MωjFUxKj,yKj,t·∇wxKj,yKjK.Owing to the fact that the flux at the boundary of the cell is discontinuous,F(U(xel,yel,t))·n→ is replaced by numerical flux F(U^(xl,yl,t))·n→. Besides, the exact solution U and the test function w are substituted with the approximate solution U^ and w^, respectively. Based on these numerical approximations, substitute (7) into (6) and we can get(8)ddt∫KU^w^x,ydΩ+∑e∈∂K∑l=1LωlFU^xel,yel,t·n→wxel,yele-∑j=1MωjFU^xKj,yKj,t·∇w^xKj,yKjK=0,where Lax-Friedrichs flux is used to get the numerical flux F(U^(xl,yl,t))·n→, given as(9)FU^xel,yel,t·n→=12FU^-xel,yel,t+FU^+xel,yel,t·n→-α2U^+xel,yel,t-U^-xel,yel,t,where α is the maximum eigenvalue of Jacobian matrix ∂F(U)/∂U.Inside the quadrilateral elementxi-1/2,xi+1/2×yj-1/2,yj+1/2, the approximate solution U^(x,y,t) can be expressed as(10)U^x,y,t=U-t+Uxtξix+Uytηjy+Uxytξixηjy+Uxxtξix2-13+Uyytηjy2-13,where ξi(x)=x-xi/Δxi/2, ηj(y)=y-yj/Δyj/2, Δxi=xi+1/2-xi-1/2, and Δyj=yj+1/2-yj-1/2.Semidiscrete formulation is obtained by substituting (9) and (10) into (8), written as(11)Ut=LU,where LU is the discrete operator of spatial derivatives. A third-order TVD-Runge-Kutta time discretization method is used to get the discrete difference scheme as shown below:(12)U1=Un+ΔtLUn,U2=34Un+14U1+ΔtLU1,Un+1=13Un+23U2+ΔtLU2.


### 2.3. BEM Based on Potential Flow Theory

During the early stage of underwater explosion, the pressure inside the bubble is dozens of times to hundreds of times larger than the hydrostatic pressure around the bubble and the velocity of bubble is pretty big. In this case, the viscosity of fluid can be neglected. As a result, the fluid can be assumed to be inviscid, irrotational, and incompressible and satisfies the Laplace Equation, which is given as [10](13)∇2Ψ=0.On the basis of Green Equation, (13) is transformed and the boundary integration equation, which is satisfied by the potential Ψ, can be obtained [10]:(14)λΨp=∬S⁡∂Ψq∂n1rpq-Ψq∂∂n1rpqdS,where λ, S, and n denote the solid angle, the boundary of the fluid filed, and the normal vector pointing out of the fluid field, respectively; p and q denote the field point and the source point on the boundary, respectively. The numerical calculation method of (14) is elaborated in [10].The kinematic boundary condition and the dynamic boundary condition of bubble surface are expressed as [28](15)dΨdt=u22+P∞ρ-Pρ-gz,(16)drdt=∇Ψx,y,z,where p∞ denotes the pressure far away from the initial charge center in the horizontal direction; ρ is the fluid density; p is the pressure inside the bubble; r is the vector of node position on the bubble surface.The pressurep at any point in the flow domain is obtained through the Bernoulli equation, in which the velocity ∇Ψ and ∂Ψ/∂t are evaluated using the indirect boundary element method [28] and the auxiliary function method [29–31], respectively.


### 2.4. A Hybrid Algorithm Based on RKDG-LS-DGF Method and BEM

Euler equations and Laplace Equation are used as the fluid motion governing equations in RKDG-LS-DGF method and BEM, respectively. The Laplace Equation is a special form of continuity equation and momentum equation in Euler equations. During the first pulsing phase and bubble jet, the time needed to let the Maher number reach 0.1 only accounts for 0.1% of the total time, and the maximum Maher number that appears on the top point of bubble jet flow merely reaches about 0.12 [32]. It indicates that the fluid compressibility only slightly affected the calculation in the pulsing phase. As a result, the underwater explosion bubble dynamics is numerically investigated in simulation bubble pulsing and bubble jet with a hybrid algorithm based on RKDG-LS-DGF method and BEM, the flowchart of which is shown in Figure 1. Parameters include node coordinate, velocity, density, pressure, and energy in the numerical simulation. The hybrid algorithm can be divided into two steps.Figure 1
Flowchart of the whole process of column charge underwater explosion using the hybrid algorithm based on RKDG-LS-DGF method and BEM.The First Step. During the process of information-delivery from RKDG-LS-DGF method to BEM, the calculation results with RKDG-LS-DGF method are taken as the initial conditions for the calculation using BEM. The initial node coordinates and velocity vectors at the boundaries in the calculation using BEM are obtained from those at the air-water interfaces got in the calculation using RKDG-LS-DGF method. And the initial pressure of bubble in the calculation using BEM is set as the average pressure of gas obtained from the calculation using RKDG-LS-DGF. Based on the initial conditions given above and BEM bubble model presented in Section 2.3, bubble position, velocity vectors, and pressure inside can be further updated. To be specific, the dynamic boundary condition of bubble surface, that is, (16), is used to calculate the bubble position. As for velocity vector of bubble surface, the normal velocity component can be got from boundary integration equation (14) and the tangential velocity component can be got through differential processing of the velocity potential. The state equation of ideal gas, that is, (3), is used to calculate the pressure inside the bubble.The Second Step. During the process of information-delivery from BEM to RKDG-LS-DGF method, the calculation results using BEM are taken as the initial conditions for the calculation using RKDG-LS-DGF method. For fluid field, the velocity and pressure calculated using BEM are regarded as the initial velocity vw and initial pressure pw for the calculation using RKDG-LS-DGF method. Initial density ρw can be got by(17)ρw=ρrpwB0+11/N,where ρr is the reference density taken as 1000kg/m3; pw denotes the initial pressure of water; B0 is used to limit the maximum change of density and the initial maximum pressure is set as 2×107Pa; N is density ratio index taken as 7 in this paper.The initial pressure of gaspa can be obtained by the relation between the pressure inside bubble p and the bubble volume V; pa can be expressed as(18)Pa=Pc+P0V0Vaγ,where Pc, P0, and V0 denote the pressure of saturated condensable vapor, the initial pressure, and the density of the bubble; γ is the specific heat capacity and it is set as 1.25 for the bubble generated by TNT explosion; the initial density ρa is obtained by ρa=W/V and the initial velocity vector va can be got by using the velocity vector linear interpolation within the fluid field around the bubble.


## 2.1. Basic Equations

Due to the high velocity on the bubble surface, the Mach number in near flow field would be very large during the early stage of expansion and bubble jet. Besides, pressure wave would be generated by the bubble. In this case, detonation products and the surrounding fluid field should be considered as inviscid, irrotational, and compressible fluid. Fluid dynamics equations for 2D compressible axisymmetric flow can be expressed as [24, 25](1)Ut+∇·FU=SsymU,where(2)U=ρρuρvρe,FU=ρuρvρu2+pρvuρuvρv2+pρue+puρve+pv,Ssym=-1xρuρu2ρuvuρe+p,where t, ρ, e, and p denote time, density, energy, and pressure, respectively; x and y indicate the direction along axis; u and v are velocity components in x direction and y direction, respectively. The state equations of different materials are listed as follows.(1) The state equation of ideal gas is given by [24, 25](3)ρe=pγ-1,where γ is the specific heat capacity of gas taken as 1.4 in the simulation of underwater explosion.(2) The Jones-Wilins-Lee (JWL) state equation is used for explosive gas [24, 25]:(4)P=A1-ωρR1ρ0e-R1ρ/ρ0+B1-ωρR2ρ0e-R2ρ/ρ0+ωρe,where A=3.712×1011Pa, B=3.231×109Pa, R1=4.15, R2=0.95, and ω=0.3; ρ0 is the initial density which is equal to 1630kg/m3; E0 is the initial internal energy which is equal to 4.29×106J/m3.(3) The Tait state equation of water is given by [24, 25](5)ρe=p+NB¯N-1,where B¯=B-A, N=7.15, A=1.0×105Pa, B=3.31×108Pa, and ρ0=1000.0kg/m3.


## 2.2. RKDG-LS-DGF Method

Column charge detonation and bubble jet are simulated using axisymmetric RKDG-LS-DGF method and three numerical methods are involved concretely. Discontinuous Galerkin method is applied to solve fluid motion equations [22, 24, 26]. Level Set method has an advantage in capturing moving interfaces of multimaterial flow [27]. These two methods are well combined by Direct Ghost Fluid method [24].The computational domain of fluid is denoted asΩ. Fluid motion equation, that is, (1), is solved by Discontinuous Galerkin method in this paper. After multiplying (1) by the test function w(x,y) on both sides, integrating over a cell Ωj, and integrating by parts, the equations can be derived as(6)ddt∫KUwx,ydΩ+∑e∈∂K∫eFU·n→wx,ydΓ-∫KFU∇wx,ydΩ=0,where n→ is the external normal vector.Line integral and surface integral are obtained by Gaussian quadrature formulas [26](7)∫eFU·n→wx,ydΓ≈∑l=1LωlFUxel,yel,t·n→wxel,yele,∫KFU·∇wx,ydΩ≈∑j=1MωjFUxKj,yKj,t·∇wxKj,yKjK.Owing to the fact that the flux at the boundary of the cell is discontinuous,F(U(xel,yel,t))·n→ is replaced by numerical flux F(U^(xl,yl,t))·n→. Besides, the exact solution U and the test function w are substituted with the approximate solution U^ and w^, respectively. Based on these numerical approximations, substitute (7) into (6) and we can get(8)ddt∫KU^w^x,ydΩ+∑e∈∂K∑l=1LωlFU^xel,yel,t·n→wxel,yele-∑j=1MωjFU^xKj,yKj,t·∇w^xKj,yKjK=0,where Lax-Friedrichs flux is used to get the numerical flux F(U^(xl,yl,t))·n→, given as(9)FU^xel,yel,t·n→=12FU^-xel,yel,t+FU^+xel,yel,t·n→-α2U^+xel,yel,t-U^-xel,yel,t,where α is the maximum eigenvalue of Jacobian matrix ∂F(U)/∂U.Inside the quadrilateral elementxi-1/2,xi+1/2×yj-1/2,yj+1/2, the approximate solution U^(x,y,t) can be expressed as(10)U^x,y,t=U-t+Uxtξix+Uytηjy+Uxytξixηjy+Uxxtξix2-13+Uyytηjy2-13,where ξi(x)=x-xi/Δxi/2, ηj(y)=y-yj/Δyj/2, Δxi=xi+1/2-xi-1/2, and Δyj=yj+1/2-yj-1/2.Semidiscrete formulation is obtained by substituting (9) and (10) into (8), written as(11)Ut=LU,where LU is the discrete operator of spatial derivatives. A third-order TVD-Runge-Kutta time discretization method is used to get the discrete difference scheme as shown below:(12)U1=Un+ΔtLUn,U2=34Un+14U1+ΔtLU1,Un+1=13Un+23U2+ΔtLU2.


## 2.3. BEM Based on Potential Flow Theory

During the early stage of underwater explosion, the pressure inside the bubble is dozens of times to hundreds of times larger than the hydrostatic pressure around the bubble and the velocity of bubble is pretty big. In this case, the viscosity of fluid can be neglected. As a result, the fluid can be assumed to be inviscid, irrotational, and incompressible and satisfies the Laplace Equation, which is given as [10](13)∇2Ψ=0.On the basis of Green Equation, (13) is transformed and the boundary integration equation, which is satisfied by the potential Ψ, can be obtained [10]:(14)λΨp=∬S⁡∂Ψq∂n1rpq-Ψq∂∂n1rpqdS,where λ, S, and n denote the solid angle, the boundary of the fluid filed, and the normal vector pointing out of the fluid field, respectively; p and q denote the field point and the source point on the boundary, respectively. The numerical calculation method of (14) is elaborated in [10].The kinematic boundary condition and the dynamic boundary condition of bubble surface are expressed as [28](15)dΨdt=u22+P∞ρ-Pρ-gz,(16)drdt=∇Ψx,y,z,where p∞ denotes the pressure far away from the initial charge center in the horizontal direction; ρ is the fluid density; p is the pressure inside the bubble; r is the vector of node position on the bubble surface.The pressurep at any point in the flow domain is obtained through the Bernoulli equation, in which the velocity ∇Ψ and ∂Ψ/∂t are evaluated using the indirect boundary element method [28] and the auxiliary function method [29–31], respectively.


## 2.4. A Hybrid Algorithm Based on RKDG-LS-DGF Method and BEM

Euler equations and Laplace Equation are used as the fluid motion governing equations in RKDG-LS-DGF method and BEM, respectively. The Laplace Equation is a special form of continuity equation and momentum equation in Euler equations. During the first pulsing phase and bubble jet, the time needed to let the Maher number reach 0.1 only accounts for 0.1% of the total time, and the maximum Maher number that appears on the top point of bubble jet flow merely reaches about 0.12 [32]. It indicates that the fluid compressibility only slightly affected the calculation in the pulsing phase. As a result, the underwater explosion bubble dynamics is numerically investigated in simulation bubble pulsing and bubble jet with a hybrid algorithm based on RKDG-LS-DGF method and BEM, the flowchart of which is shown in Figure 1. Parameters include node coordinate, velocity, density, pressure, and energy in the numerical simulation. The hybrid algorithm can be divided into two steps.Figure 1
Flowchart of the whole process of column charge underwater explosion using the hybrid algorithm based on RKDG-LS-DGF method and BEM.The First Step. During the process of information-delivery from RKDG-LS-DGF method to BEM, the calculation results with RKDG-LS-DGF method are taken as the initial conditions for the calculation using BEM. The initial node coordinates and velocity vectors at the boundaries in the calculation using BEM are obtained from those at the air-water interfaces got in the calculation using RKDG-LS-DGF method. And the initial pressure of bubble in the calculation using BEM is set as the average pressure of gas obtained from the calculation using RKDG-LS-DGF. Based on the initial conditions given above and BEM bubble model presented in Section 2.3, bubble position, velocity vectors, and pressure inside can be further updated. To be specific, the dynamic boundary condition of bubble surface, that is, (16), is used to calculate the bubble position. As for velocity vector of bubble surface, the normal velocity component can be got from boundary integration equation (14) and the tangential velocity component can be got through differential processing of the velocity potential. The state equation of ideal gas, that is, (3), is used to calculate the pressure inside the bubble.The Second Step. During the process of information-delivery from BEM to RKDG-LS-DGF method, the calculation results using BEM are taken as the initial conditions for the calculation using RKDG-LS-DGF method. For fluid field, the velocity and pressure calculated using BEM are regarded as the initial velocity vw and initial pressure pw for the calculation using RKDG-LS-DGF method. Initial density ρw can be got by(17)ρw=ρrpwB0+11/N,where ρr is the reference density taken as 1000kg/m3; pw denotes the initial pressure of water; B0 is used to limit the maximum change of density and the initial maximum pressure is set as 2×107Pa; N is density ratio index taken as 7 in this paper.The initial pressure of gaspa can be obtained by the relation between the pressure inside bubble p and the bubble volume V; pa can be expressed as(18)Pa=Pc+P0V0Vaγ,where Pc, P0, and V0 denote the pressure of saturated condensable vapor, the initial pressure, and the density of the bubble; γ is the specific heat capacity and it is set as 1.25 for the bubble generated by TNT explosion; the initial density ρa is obtained by ρa=W/V and the initial velocity vector va can be got by using the velocity vector linear interpolation within the fluid field around the bubble.


## 3. Numerical Model, Results, and Discussion

### 3.1. Numerical Model

The numerical results in this paper are compared with the experimental data in [33]. The numerical model is established according to the experimental model, shown as Figure 2. The column TNT charge, 5.8 g, is located at a depth of 1 m in the square water. The TNT column has a side length h=0.016m and a diameter r=0.008m. The length a and width b of the water area are 0.8 m and 1.0 m, respectively. The coordinate system is marked in Figure 1, and the position of test point A is (0.5,0.5). The quadrilateral meshes are even distributed. There are 200000 grid cells in total and the size of grid cells is 0.002 m.Figure 2
Model of column charge underwater explosion.


### 3.2. Charge Detonation Simulation with RKDG-LS-DGF Method

The RKDG-LS-DGF method is used to simulate the column charge detonation. The pressure nephograms are shown in Figure3. According to the figure, a shock wave with extremely high pressure is generated and it rapidly propagates into the water. The maximum pressure is about 480MPa at t≈0.01ms, shown in Figure 3(a). Besides, the column charge becomes an ellipsoidal detonation product rapidly. With the shock wave propagating into the water, the shock wave peak value decays exponentially with the increase of distance. The pressure decreases by one order of magnitude and rapidly to 36MPa at t≈0.13ms, shown in Figure 3(b), and the shape of the shock wave develops from an ellipsoid into a sphericity. Subsequently, as Figure 3(c) has shown the pressure decreases to 16MPa at t≈0.26ms. The high-pressure gas generated by the explosion forms a spherical bubble whose diameter is about 0.1167 m in the water at this time. In addition, as shown in Figure 3, nonphysical condition, such as oscillation, has not occurred close to the symmetry axis. This phenomenon sufficiently indicates that the validity and feasibility of axisymmetric RKDG-LS-DGF model can be verified to simulate the detonation of column charge.Figure 3
Pressure nephograms after column charge detonation (from left to right the times are (a) 0.01 ms, (b) 0.13 ms, and (c) 0.26 ms).
(a)
(b)
(c)The pressure-time curve at test point A is shown in Figure4, in which the solid line and the dashed line denote the experimental results and the numerical results, respectively. At t≈0.30ms, the shock wave reaches the test point A and the pressure increases to about 13MPa. Subsequently, the pressure in the water decreases exponentially as time goes by. In view of the tendency of the pressure-time curve, it can be seen that the pressure calculated in this paper agrees well with the propagation law of shock wave. Furthermore, by comparing the numerical results with the experimental results, it is obvious that the tendency of the numerical results is in good agreement with that of the experimental results, and the error of the pressure peak value between these two results is less than 10%. In conclusion, the pressure-time curve obtained in the simulation agrees well with that obtained from the experimental data, and it further verifies the validity of the axisymmetric RKDG-LS-DGF numerical model in the numerical simulation of charge detonation.Figure 4
Pressure-time curve at test point A (numerical results are compared with experimental data).


### 3.3. Bubble Pulsing Simulation with BEM

On account of the advantages of BEM in the simulation of bubble pulsing, such as improving the efficiency and saving computation time, a bubble pulsing BEM model is established to simulate the early stage nonspherical bubble pulsing. The previous results obtained from the calculation using RKDG-LS-DGF are set as the initial condition to simulate bubble pulsing. And the match of initial conditions can be divided into three aspects, which are listed as follows: (1) the node coordinates on the gas-liquid interface att≈0.26ms calculated using the RKDG-LS-DGF method are set as the node coordinates on bubble boundary; (2) the velocity vectors on the gas-liquid interface at t≈0.26ms calculated using the RKDG-LS-DGF method are used as velocity vectors on bubble boundary; (3) the average pressure of gas is set as the initial pressure of bubble at the same time.The process of bubble pulsing is shown in Figure5, in which the left one is the velocity vector distributions and the right one is the pressure nephograms. Additionally, the blank area in the center of the figure is the bubble. The initial condition of bubble pulsing is shown in Figure 5(a). It can be known that the initial pressure of the bubble is larger than that of the surrounding fluid field and the direction of bubble velocity is the same with the normal vector on the surface pointing outward. The process of bubble expansion is shown in Figures 5(b) and 5(c). The bubble expansion is driven by the high pressure inside. And the pressure inside the bubble will gradually decrease as the bubble expands. Finally, the size of the bubble will reach a maximum value at a specific time which is shown in Figure 5(c). Subsequently as shown in Figure 5(d), the bubble is in contraction phase, the reason of which is that the pressure inside the bubble is lower than the static pressure in the surrounding fluid field. Besides, as a result of gravity, the pressure on the upper surface of bubble is evidently smaller than that on the lower surface. The process of bubble collapse is shown in Figure 5(e). The shape of the bubble would no longer maintain a spherical one under the effect of gravity and the lower surface of bubble will collapse into the bubble. Figure 5(f) shows the process of bubble jet. According to the figure, there is a jet generated in the opposite direction of gravity.Figure 5
Bubble pulsing process (the left one is the velocity vector distributions and the right one is the pressure nephograms; from left to right the times are (a)–(c) 0.26 ms, 0.86 ms, and 22.48 ms; (d)–(f) 42.76 ms, 44.43 ms, and 44.69 ms).
(a)
(b)
(c)
(d)
(e)
(f)The curve of bubble radius during bubble pulsing obtained from the numerical calculation is compared with that obtained from the experimental data, shown in Figure6. The comparison indicates that the tendency of the former one is in good agreement with that of the latter one. Besides, the numerical bubble radius at any specific time is nearly the same with the experimental measure. In addition, the maximum bubble radiuses obtained from calculation and that obtained from experiments are 0.278m and 0.270m, respectively. The error between them is less than 5%. The displacement curves of upper surface, lower surface, and the center of bubble are plotted in Figure 7. It is obvious that the displacement curves between calculation and experiments share a similar tendency. The displacement curve of the lower surface becomes steeper during the later stage of collapse with increasing collapse speed, and the bubble ascends quickly under the effect of buoyancy.Figure 6
Comparison curves of bubble radius along with time during bubble pulsing.Figure 7
Comparison curves of bubble upper surface, lower surface, and center position along with time.Through the comparison between calculation and experiment, the following conclusions can be obtained. After the detonation, the bubble experiences expansion and contraction and develops a jet, and the whole process obeys the basic theories of bubble motion. Numerical results, including bubble position, size, shape, and radius, are consistent with the experimental results, which verifies the validity of using BEM in the simulation of bubble pulsing.


### 3.4. Bubble Jet Simulation with RKDG-LS-DGF Method

The bubble distorts drastically during jet phase. Additionally, the fluid field has transformed from a simply connected domain to doubly connected domain since the bubble is penetrated by the jet with a high speed. The temporality velocity potential of the fluid field is no longer a uniform function and has been converted into a multivalue function. Some special numerical techniques, such as Vortex-Ring/sheet model, regulating grids, or boundary nodes, should be imported when traditional BEM is applied to simulate the process of bubble jet. In this case, the compressibility of fluid field cannot be taken into consideration because BEM with Vortex-Ring/sheet model is based on potential flow assumption. However, shock waves would be generated by the bubble when the jet is developed and the bubble volume reaches the maximum in collapse phase. In that case the compressibility of fluid cannot be neglected. Therefore, the RKDG-LS-DGF method is applied to simulate the later stage of bubble jet.The previous results obtained from calculation using BEM are set as the initial conditions in the simulation of bubble jet. The conditions are divided into two aspects. On one hand, the parameters of fluid field should be obtained. The initial velocity and the pressure of the fluid field calculated using BEM are regarded as the initial velocityvw and the initial pressure pw. The initial density ρw can be got by (17). On the other hand, the parameters of gas should be obtained. Equation (18) and equation ρa=W/V are used to get the initial pressure pa and the initial density ρa. The initial velocity vector va can be got by using the velocity vector linear interpolation within the fluid field around the bubble.As shown in Figure8, the process of bubble jet is simulated with RKDG-LS-DGF method. The early stage of bubble jet is shown in Figure 8(a). It can be seen that pressure at the upper side of bubble is higher than that at the lower side, which leads to the formation of pressure gradient. Besides, the lower side of bubble is concave and the bubble is simply connected. Subsequently, the jet threads through the bubble in the opposite direction of gravity. The singly connected domain changes into a multiply connected domain and the shape of the bubble becomes a toroidal one, as shown in Figure 8(b). In addition, the high-pressure region is induced in the jet area and the pressure peak value is about 70 MPa. Then Figures 8(c) and 8(d) show that the bubble experiences expansion and contraction, respectively. It can be known that bubble volume deceases continuously in fluid field. When the bubble radius reaches the minimum value, the pressure inside the bubble is higher than that in the fluid field. As a result, the bubble as shown in Figure 8(d) expands for the second time. It is obvious that the pressure of fluid field surrounding the jet is higher than that of the exterior flow field, and the pressure peak value increases to about 12 MPa at about t≈45.08ms. Consequently, the high speed and the high pressure of the bubble jet are so high that secondary damage to the structure may occur. Additionally, it is evident that the numerical calculation results agree well with the experimental results in terms of motion trend and deformation law. The whole process obeys the basic theory of bubble motion, which verifies the validity and feasibility in using RKDG-LS-DGF to conduct the numerical simulation of bubble jet.Figure 8
Pressure nephograms during bubble jet phase (from left to right the times are (a) 44.69 ms, (b) 44.82 ms, (c) 44.95 ms, and (d) 45.08 ms).
(a)
(b)
(c)
(d)


## 3.1. Numerical Model

The numerical results in this paper are compared with the experimental data in [33]. The numerical model is established according to the experimental model, shown as Figure 2. The column TNT charge, 5.8 g, is located at a depth of 1 m in the square water. The TNT column has a side length h=0.016m and a diameter r=0.008m. The length a and width b of the water area are 0.8 m and 1.0 m, respectively. The coordinate system is marked in Figure 1, and the position of test point A is (0.5,0.5). The quadrilateral meshes are even distributed. There are 200000 grid cells in total and the size of grid cells is 0.002 m.Figure 2
Model of column charge underwater explosion.


## 3.2. Charge Detonation Simulation with RKDG-LS-DGF Method

The RKDG-LS-DGF method is used to simulate the column charge detonation. The pressure nephograms are shown in Figure3. According to the figure, a shock wave with extremely high pressure is generated and it rapidly propagates into the water. The maximum pressure is about 480MPa at t≈0.01ms, shown in Figure 3(a). Besides, the column charge becomes an ellipsoidal detonation product rapidly. With the shock wave propagating into the water, the shock wave peak value decays exponentially with the increase of distance. The pressure decreases by one order of magnitude and rapidly to 36MPa at t≈0.13ms, shown in Figure 3(b), and the shape of the shock wave develops from an ellipsoid into a sphericity. Subsequently, as Figure 3(c) has shown the pressure decreases to 16MPa at t≈0.26ms. The high-pressure gas generated by the explosion forms a spherical bubble whose diameter is about 0.1167 m in the water at this time. In addition, as shown in Figure 3, nonphysical condition, such as oscillation, has not occurred close to the symmetry axis. This phenomenon sufficiently indicates that the validity and feasibility of axisymmetric RKDG-LS-DGF model can be verified to simulate the detonation of column charge.Figure 3
Pressure nephograms after column charge detonation (from left to right the times are (a) 0.01 ms, (b) 0.13 ms, and (c) 0.26 ms).
(a)
(b)
(c)The pressure-time curve at test point A is shown in Figure4, in which the solid line and the dashed line denote the experimental results and the numerical results, respectively. At t≈0.30ms, the shock wave reaches the test point A and the pressure increases to about 13MPa. Subsequently, the pressure in the water decreases exponentially as time goes by. In view of the tendency of the pressure-time curve, it can be seen that the pressure calculated in this paper agrees well with the propagation law of shock wave. Furthermore, by comparing the numerical results with the experimental results, it is obvious that the tendency of the numerical results is in good agreement with that of the experimental results, and the error of the pressure peak value between these two results is less than 10%. In conclusion, the pressure-time curve obtained in the simulation agrees well with that obtained from the experimental data, and it further verifies the validity of the axisymmetric RKDG-LS-DGF numerical model in the numerical simulation of charge detonation.Figure 4
Pressure-time curve at test point A (numerical results are compared with experimental data).


## 3.3. Bubble Pulsing Simulation with BEM

On account of the advantages of BEM in the simulation of bubble pulsing, such as improving the efficiency and saving computation time, a bubble pulsing BEM model is established to simulate the early stage nonspherical bubble pulsing. The previous results obtained from the calculation using RKDG-LS-DGF are set as the initial condition to simulate bubble pulsing. And the match of initial conditions can be divided into three aspects, which are listed as follows: (1) the node coordinates on the gas-liquid interface att≈0.26ms calculated using the RKDG-LS-DGF method are set as the node coordinates on bubble boundary; (2) the velocity vectors on the gas-liquid interface at t≈0.26ms calculated using the RKDG-LS-DGF method are used as velocity vectors on bubble boundary; (3) the average pressure of gas is set as the initial pressure of bubble at the same time.The process of bubble pulsing is shown in Figure5, in which the left one is the velocity vector distributions and the right one is the pressure nephograms. Additionally, the blank area in the center of the figure is the bubble. The initial condition of bubble pulsing is shown in Figure 5(a). It can be known that the initial pressure of the bubble is larger than that of the surrounding fluid field and the direction of bubble velocity is the same with the normal vector on the surface pointing outward. The process of bubble expansion is shown in Figures 5(b) and 5(c). The bubble expansion is driven by the high pressure inside. And the pressure inside the bubble will gradually decrease as the bubble expands. Finally, the size of the bubble will reach a maximum value at a specific time which is shown in Figure 5(c). Subsequently as shown in Figure 5(d), the bubble is in contraction phase, the reason of which is that the pressure inside the bubble is lower than the static pressure in the surrounding fluid field. Besides, as a result of gravity, the pressure on the upper surface of bubble is evidently smaller than that on the lower surface. The process of bubble collapse is shown in Figure 5(e). The shape of the bubble would no longer maintain a spherical one under the effect of gravity and the lower surface of bubble will collapse into the bubble. Figure 5(f) shows the process of bubble jet. According to the figure, there is a jet generated in the opposite direction of gravity.Figure 5
Bubble pulsing process (the left one is the velocity vector distributions and the right one is the pressure nephograms; from left to right the times are (a)–(c) 0.26 ms, 0.86 ms, and 22.48 ms; (d)–(f) 42.76 ms, 44.43 ms, and 44.69 ms).
(a)
(b)
(c)
(d)
(e)
(f)The curve of bubble radius during bubble pulsing obtained from the numerical calculation is compared with that obtained from the experimental data, shown in Figure6. The comparison indicates that the tendency of the former one is in good agreement with that of the latter one. Besides, the numerical bubble radius at any specific time is nearly the same with the experimental measure. In addition, the maximum bubble radiuses obtained from calculation and that obtained from experiments are 0.278m and 0.270m, respectively. The error between them is less than 5%. The displacement curves of upper surface, lower surface, and the center of bubble are plotted in Figure 7. It is obvious that the displacement curves between calculation and experiments share a similar tendency. The displacement curve of the lower surface becomes steeper during the later stage of collapse with increasing collapse speed, and the bubble ascends quickly under the effect of buoyancy.Figure 6
Comparison curves of bubble radius along with time during bubble pulsing.Figure 7
Comparison curves of bubble upper surface, lower surface, and center position along with time.Through the comparison between calculation and experiment, the following conclusions can be obtained. After the detonation, the bubble experiences expansion and contraction and develops a jet, and the whole process obeys the basic theories of bubble motion. Numerical results, including bubble position, size, shape, and radius, are consistent with the experimental results, which verifies the validity of using BEM in the simulation of bubble pulsing.


## 3.4. Bubble Jet Simulation with RKDG-LS-DGF Method

The bubble distorts drastically during jet phase. Additionally, the fluid field has transformed from a simply connected domain to doubly connected domain since the bubble is penetrated by the jet with a high speed. The temporality velocity potential of the fluid field is no longer a uniform function and has been converted into a multivalue function. Some special numerical techniques, such as Vortex-Ring/sheet model, regulating grids, or boundary nodes, should be imported when traditional BEM is applied to simulate the process of bubble jet. In this case, the compressibility of fluid field cannot be taken into consideration because BEM with Vortex-Ring/sheet model is based on potential flow assumption. However, shock waves would be generated by the bubble when the jet is developed and the bubble volume reaches the maximum in collapse phase. In that case the compressibility of fluid cannot be neglected. Therefore, the RKDG-LS-DGF method is applied to simulate the later stage of bubble jet.The previous results obtained from calculation using BEM are set as the initial conditions in the simulation of bubble jet. The conditions are divided into two aspects. On one hand, the parameters of fluid field should be obtained. The initial velocity and the pressure of the fluid field calculated using BEM are regarded as the initial velocityvw and the initial pressure pw. The initial density ρw can be got by (17). On the other hand, the parameters of gas should be obtained. Equation (18) and equation ρa=W/V are used to get the initial pressure pa and the initial density ρa. The initial velocity vector va can be got by using the velocity vector linear interpolation within the fluid field around the bubble.As shown in Figure8, the process of bubble jet is simulated with RKDG-LS-DGF method. The early stage of bubble jet is shown in Figure 8(a). It can be seen that pressure at the upper side of bubble is higher than that at the lower side, which leads to the formation of pressure gradient. Besides, the lower side of bubble is concave and the bubble is simply connected. Subsequently, the jet threads through the bubble in the opposite direction of gravity. The singly connected domain changes into a multiply connected domain and the shape of the bubble becomes a toroidal one, as shown in Figure 8(b). In addition, the high-pressure region is induced in the jet area and the pressure peak value is about 70 MPa. Then Figures 8(c) and 8(d) show that the bubble experiences expansion and contraction, respectively. It can be known that bubble volume deceases continuously in fluid field. When the bubble radius reaches the minimum value, the pressure inside the bubble is higher than that in the fluid field. As a result, the bubble as shown in Figure 8(d) expands for the second time. It is obvious that the pressure of fluid field surrounding the jet is higher than that of the exterior flow field, and the pressure peak value increases to about 12 MPa at about t≈45.08ms. Consequently, the high speed and the high pressure of the bubble jet are so high that secondary damage to the structure may occur. Additionally, it is evident that the numerical calculation results agree well with the experimental results in terms of motion trend and deformation law. The whole process obeys the basic theory of bubble motion, which verifies the validity and feasibility in using RKDG-LS-DGF to conduct the numerical simulation of bubble jet.Figure 8
Pressure nephograms during bubble jet phase (from left to right the times are (a) 44.69 ms, (b) 44.82 ms, (c) 44.95 ms, and (d) 45.08 ms).
(a)
(b)
(c)
(d)


## 4. Conclusions

During the collapse phase, the Mach number at the top of bubble jet would be very large, in which case the compressibility of the fluid has a significant influence on the bubble motion. Therefore, the compressibility should be taken into consideration to improve the accuracy. In this paper a hybrid algorithm based on RKDG-LS-DGF method and BEM is applied to simulate the whole process of such complex physical phenomenon as the column charge underwater explosion. To be specific, the underwater explosion bubble dynamics from detonation to bubble jet is investigated in this paper, including the high-pressure shock wave propagation, the bubble pulsing, and the high-speed jet. Finally, the numerical results are compared with the experimental results to verify the validity and feasibility of the hybrid algorithm. The conclusions are as follows:(1)
A RKDG-LS-DGF model of column charge detonation is established. The results indicate that shock wave peak value decays exponentially with the increase of distance, which is in good accordance with the physics law of shock wave. Moreover, the pressure-time curves obtained from calculation and experiments show good agreements with each other and the error of pressure peak between them is within 10%, which all prove the effectiveness and accuracy of RKDG-LS-DGF method to solve the problem of underwater explosion.(2)
The BEM is applied to simulate the process of bubble expansion, collapse, jet, and so forth, due to its high calculating efficiency in the simulation of bubble pulsing process. The numerical results of bubble pulsing are in good agreement with the experimental results, which obeys basic law of bubble motion.(3)
The process of the jet penetrating the bubble is simulated by an axisymmetric RKDG-LS-DGF method, which avoids the manual interference, such as mesh cutting and smoothing, in the BEM. The jet threads through the bubble in the opposite direction of gravity and induces a high-pressure region. The numerical results agree well with the experimental results and meet the basic laws of bubble jet.(4)
A hybrid algorithm that fully utilizes the advantages of RKDG-LS-DGF method and BEM is applied to simulate the whole process of column charge underwater explosion. The numerical results show good agreement with the experimental results, which verifies the feasibility and effectiveness of the hybrid algorithm. Besides, the combined algorithm is suitable to solve other similar hydrodynamic problems as well.



---

*Source: 165252-2015-07-15.xml*