# A Direct Ghost Fluid Method for Modeling Explosive Gas and Water Flows

Authors: Nan Si; Jinwon Park; Alan J. Brown
Journal: Shock and Vibration
 (2022)

---

## Abstract

This work presents a Direct Ghost Fluid Method (DGFM) as part of a two-fluid numerical framework suitable to model explosive gas and water flows resulting from underwater explosion (UNDEX). Due to the presence of explosive gas and water with shock waves in the modeling domain, classic Eulerian methods with inherent diffusion may not be effective. Numerical diffusion occurs due to nonphysical diffused density at material interfaces, which creates spurious pressure oscillations and significantly degrades the quality of the numerical results. To eliminate or minimize numerical diffusion, sharp interface methods having no mixed elements may be used in multifluid flow computations. The Direct Ghost Fluid Method (DGFM) described in this paper uses direct extrapolation of density (vice pressure) and tangential velocity from real to ghost fluid. The spurious pressure oscillations near the material interface are therefore minimized. One-, two-, and three-dimensional computational fluid dynamics (CFD) solvers that have DGFM as an essential part in their framework to model UNDEX interface conditions are developed, explored, and applied to the simulation of a series of benchmark problems. Excellent agreement is obtained among the simulations, the analytical solutions, and the experiments.


---

## Body

## 1. Introduction

An underwater explosion (UNDEX) creates a high-pressure gas bubble and a complex fluid flow that consists of explosive gas and water with shock waves. This is particularly complex in proximity to a structure and/or free surface where cavitation can be generated by reflected waves. This poses a difficult modeling problem. The following assumptions are applied to this problem: (1) the explosive gas and water are individually homogeneous and separated by a well-defined material interface; (2) the resulting fluid flow is compressible as well as inviscid since the time scale of UNDEX is extremely small so that diffusion can be neglected; and (3) the explosion is adiabatic because the simulation starts after the chemical reaction which turns the explosives into high-pressure gaseous materials is completed so that there is no external heat source throughout the simulation [1, 2]. The fluid flow may be modeled using Euler equations, and the thermodynamic properties of the fluids across the material interface may be described by their respective equations of state (EOS), but the direct application of classical Eulerian methods with inherent numerical diffusion may not be effective [3]. Numerical diffusion causes nonphysical density differences at the material interface to diffuse from one fluid to the other [4]. Mixed elements that simultaneously possess properties of both fluids within one element are generated [5]. Since pressure is determined by a thermodynamic relation, diffused density may produce spurious pressure oscillations in the region near the material interface with a discontinuity of the EOS. Numerical diffusion can be a significant source of nonphysical results and often cause run-time errors [4]. To eliminate or minimize numerical diffusion, multiple types of sharp interface methods such as the Ghost Fluid Method (GFM) and the Lagrangian Interface Method (LIM) have proven useful in multifluid flow computations [3].Fedkiw developed the GFM to numerically simulate multifluid flows that include gas-gas, gas-liquid, and liquid-gas interactions [5–8]. Compared to earlier methods [4, 9–11], the GFM is more robust and easier to extend to multidimensions without geometric complexities. However, it may not be useable when the material interface interacts with a strong shock as in UNDEX fluid flows [12, 13]. For explosive gas and water flows where the density may vary rapidly, direct application of the original GFM may be less efficient and generate numerical inaccuracies. Liu et al. stated that interface conditions set by the original GFM may still cause spurious pressure oscillations at the material interface [3, 12]. In early studies by the authors, the original GFM failed to accurately simulate a gas-water shock tube flow because a severe density overshoot occurred early as shown in Figure 1. It is believed that the process of first extrapolating entropy across the material interface and then restoring density as in the original GFM results in run-time errors and thus failure of the computation [14].Figure 1
GFM shock tube density overshoot [14].Liu et al. presented a Modified Ghost Fluid Method (MGFM) using local shock values to correct the affected cell values in the region near the material interface [3, 12]. The basis of the MGFM is similar to the Godunov Prediction-Correction Method [4]. By analytically applying local shock values, the MGFM has wider application for various compressible fluid flows, but with inherent disadvantages: the extension to multidimensions is not trivial because of geometric complexities; and the use of general equations of state may require complex analytic computations comparable with the original GFM.In addition to the ghost fluid types of sharp interface methods, there are other Eulerian types of methods, such as the front-tracking method. This method uses marker points to track the interface. The location of marker points is decided via the solution of fluid equations [15–17]. As a cut-cell type of method, it comes with a problem that the elements cut by interfaces may contain small subcells, resulting in extremely small time step size in an explicit time integration scheme in order to meet the Courant–Friedrichs–Levy (CFL) requirement [18]. Also, the complicated front-tracking method is not always effective especially when a large interface deformation exists [18].Unlike the Eulerian methods, the Lagrangian Interface Method (LIM) explicitly tracks the motion of the interface by defining it as a Lagrangian surface. Since the material interface is treated in a Lagrangian manner, the LIM provides accurate interface motion and does not create mixed elements, which may be a source of error as in Eulerian methods. However, if the LIM is used to simulate cases where the material interface is subjected to large deformation as in UNDEX fluid flows, the LIM may experience severe mesh distortion since the mesh is fixed on the fluid throughout the computation. The error associated with distortion of the mesh may cause the degradation of numerical results and the sudden breakdown of the computation [19–21]. This makes the application of the LIM to the simulation of UNDEX fluid flows difficult. A variant GFM with a problem-dependent interface treatment suitable in the analysis of explosive gas and water flows would be more suitable. The Direct Ghost Fluid Method using direct extrapolation of density and tangential velocity across the material interface is developed for this purpose. This approach overcomes the difficulties experienced in the original GFM, when it treats multiphase flows in UNDEX problems. It avoids numerical inaccuracies and density overshoot when treating UNDEX flows, where density varies rapidly, or the material interface interacts with strong shock waves. The objective of this approach is to extend the applicability of the GFM to the analysis of UNDEX fluid flows in which the original GFM may be less effective. The scheme must also effectively be applied to fluid-structure interaction and consider a free-surface boundary. In this research, mathematical details regarding imposition of the DGFM, especially three-dimensional separation of normal and tangential velocities before extrapolation, is first developed. The effectiveness of multiple methods for this, both from the literature and devised by the authors, is tested. Multiple cases with different configurations are simulated using the solver with DGFM to further prove its capability.


## 2. Numerical Method

### 2.1. Governing Equations

The Euler equations in the conservative form are used as the fluid governing equations. The Arbitrary Lagrangian–Eulerian (ALE) form of governing equations are listed in the following equations.(1)∂ρ∂t+∇⋅ρm⟶=0,(2)∂ρv⟶∂t+∇.ρv⟶⊗m⟶+pI==0,(3)∂ρE∂t+∇⋅ρEm⟶+pv⟶=0,where the terms ρ,v⟶,m⟶=v⟶−v⟶grid,I=, and E are density, velocity, convection velocity, identity tensor, and total energy per unit mass, respectively. There are source terms on the right-hand side of equations (1)–(3) when a 1D cylindrical, 1D spherical, or 2D axisymmetric analysis is being performed. There are also source terms if body forces exist in the model configuration. The difference between velocities of fluid flow, v⟶, and velocities of the moving grid, v⟶grid, is the convection velocity, m⟶. The convection velocity equals the velocity of fluid flow in a pure Eulerian description because the grid is fixed in space. The convection velocity equals zero; however, in a pure Lagrangian description, because the grid moves together with the fluid flow, the ALE description stays between the Eulerian and the Lagrangian descriptions.This system is closed by the equation of state. The Johns–Wilkins–Lee (JWL) EOS is primarily used to model the explosive gas, but the ideal gas law is used for a comparison to some analytic results [2, 3, 6, 22]. The JWL EOS is(4)p=A1−ωηR1e−R1/η+B1−ωηR2e−R2/η+ωρeint,where A,B,R1,R2, and ω are experimentally determined coefficients and ρ, eint, and η are density, internal energy per unit mass, and the ratio of current density to the initial density, respectively. Liu demonstrated that JWL EOS can predict peak pressures approximating that of detonation experiments, while the ideal gas law can predict good results only beyond the region of 8 ∼ 10 times the charge radius [1, 3]. Despite this, the ideal gas law is often used in analytic solutions because of its common use in gas dynamics [23, 24] and its commonality with the stiffened gas EOS and Tait EOS used for water.The stiffened gas EOS is used to model the thermodynamic properties of water:(5)p=γB−1ρeint−γBB,where γB and B are experimentally determined coefficients. We have also used the Tait EOS:(6)p=N−1ρeint−NB−A,where N,B, and A are experimentally determined coefficients. The Grüneisen EOS is a more general EOS, capable of modeling a larger class of fluids. When the density change is small, the Grünesien EOS reduces to the stiffened gas EOS [22, 25].


### 2.2. Discretization

The Runge–Kutta Discontinuous Galerkin (RKDG) method is chosen as the numerical scheme to solve equations (1)–(3). The Discontinuous Galerkin (DG) method [26–28] discretizes the equations in space and the Runge–Kutta (RK) method discretizes in time. P2 Legendre polynomials are chosen as the basis functions in the DG method. This makes the spatial discretization formally third order accurate [14]. A three-step RK method is utilized whose formal order of accuracy is compatible with that of the spatial discretization [14, 28]. The DG method combines advantages of the Finite Element Method (FEM) and the Finite Volume Method (FVM). A local weak form provides weak connectivity between elements through the boundary integration of flux terms similar to the FVM. This relaxes the element-by-element continuity requirement. The remaining terms are treated using the standard Galerkin method. The p2 Legendre polynomial basis functions provide high-order accuracy and diagonalize the mass matrix without lumping. Simple extension of numerical flux functions allows variable discontinuities. Decoupled time derivatives between elements enable straightforward implementation of time marching algorithms, such as the explicit RK method. Details of one- and two-dimensional DG discretizations are provided in [14]. The three-dimensional DG discretization, verification, and validation are documented in a preceding publication (Si and Brown, 2021).In summary, the RKDG method as implemented here provides high accuracy and weak connectivity between elements and is easily extendable to multidimensions. Because of these characteristics, the RKDG method is competent in compressible flow analysis especially with strong shocks. The high accuracy of the RKDG method helps to produce less-diffusive numerical results compared to low-order approximation methods and helps to minimize inherent numerical diffusion, which is a source of error in Eulerian fluid methods. The order-of-accuracy verification of the 1D solver using a smooth configuration with sinusoidal initial conditions is described in [14]. A good match between observed and formal order of accuracy is indicated.


### 2.3. Multiphase and the Direct Ghost Fluid Method (DGFM)

To numerically treat the material interface using the GFM, the motion of the interface must be tracked at each time step during the computation [8, 29]. The GFM as implemented by Osher and Fedkiw used the Level Set Method (LSM) to perform this task. The LSM provides the geometry of the interface such as its curvature and the normal vector at different locations [29]. Because of its Eulerian nature, the implementation of the LSM does not increase geometric difficulty in multidimensions [29, 30]. The material interface is embedded in a spatially fixed grid. The transport of the LS equations is performed in an Eulerian manner [29, 30]. The transport and reinitialization equations are(7)∂ϕ∂t+v⟶⋅∇ϕ=0,(8)∂ϕ∂τ+Sϕ0∇ϕ−1=0,where ϕ is the level set (LS) function and Sϕ is a Heaviside function. Equation (7) implies that the interface is transported by the flow velocity. The function ϕ is the signed distance function ϕ=±d. It is set in the initial condition as the shortest normal distance (d) from any point in the domain to the material interface. In the LS function, a positive distance ϕ=+d indicates that a point resides in the water, a negative distance ϕ=−d indicates that a point resides in the explosive gas, and a zero distance ϕ=0 indicates that a point resides on the interface [29, 30]. The time-varying shape of the material interface is therefore evolved through the numerical solution of equation (7). Since the LS function ϕ must be transported by a real nonuniform flow field, the LS function may become distorted if not frequently reinitialized to keep the LS function real [29]. A nonuniform flow field generates noise that is not appropriate because the numerical results with distorted LS function may become unreasonable at subsequent time steps [29, 31]. The Hamilton–Jacobi WENO scheme and three-step RK method are used to numerically solve this transport equation and reinitialize the distorted LS function during the computation [32].In the GFM, elements possess unique conservative properties (i.e., mass, momentum and energy in both the real fluid and the ghost fluid [29, 30]). Once the conservative properties in the ghost fluid are determined, the two-fluid flows are solved separately using the RKDG method. In the solution of the fluid, there is no spatial concern for the interface. Solution of the fluid equations using the RKDG method provides the flow velocity v⟶. The LS function ϕ is transported in the subsequent time step using the computed flow velocity with the output of the fluid calculation in the last time step. Since pressure and normal velocity are continuous at the interface except for the initial step when pressure is discontinuous across the interface, real pressure and normal velocity are copied directly into ghost values in a simple node-to-node manner [21, 30] as described in Figure 2. In the original GFM, based on the isentropic assumption, entropy is copied by a one-sided constant extrapolation across the interface, which then allows the density to be calculated indirectly.Figure 2
Extrapolation of density in 1D implementation of DGFM for explosive gas and water flow simulations [6].Figure2 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].In the DGFM implementation, density is used instead of the entropy as in the original GFM [29] to avoid overshoot as shown in Figure 1. Density and tangential velocity is directly extrapolated from the real water to the ghost water (i.e., the area occupied by real gas) rather than using an indirect variable extrapolation from entropy. The density and tangential velocity are extrapolated using the isobaric fix technique from the real gas to the ghost gas (i.e., the area occupied by real water) minimizing the overheating error [29]. The isobaric fix technique is often included in the implementation of the GFM to minimize nonphysical temperature and density overshoot when a compressible flow interacts with a stiff medium [7, 8]. To distinguish this procedure from the indirect variable extrapolation in the original GFM, this direct extrapolation method is called the Direct Ghost Fluid Method. Diffusion of density is decreased by the DGFM across the material interface between the explosive gas and the water. The spurious pressure oscillation at the interface is therefore significantly reduced. Extrapolation of density in the 1D implementation of DGFM is illustrated in Figure 2. The density is extended from the real fluid into the ghost fluid by numerically solving the following advection equation:(9)∂ρ∂τ±n⟶⋅∇ρ=0,where n⟶=∇ϕ/∇ϕ is obtained from the LS function [29]. A few iterations, say 25, are required to be performed. The positive sign is used when density and tangential velocity are extended from the real gas (ϕ<0) to the ghost gas (ϕ>0), while the negative sign is used when they are extended from real water (ϕ>0) to ghost water (ϕ<0) [6].It is worth noting that there is an important detail in the implementation of the variable extrapolation by solving equation (9). In a 2D simulation, it is straightforward that the tangential vector is obtained by rotating the normal vector so that t⟶=−nyi^+nxj^. The real velocity field is thus decomposed into normal and tangential velocity fields vn=v⟶⋅n⟶ and vt=v⟶⋅t⟶. Tangential velocity of the real area is then extended to the ghost area as vt′ by iterating equation (9) so that total velocity of the ghost fluid is a combination of real normal velocity and ghost tangential velocity as v⟶ghost=vnn⟶+vt′t⟶. In a 3D simulation, however, there exists a pair of orthonormal bases to span the tangential plane that is perpendicular to the normal vector n⟶; for example,(10)t⟶1=nznx2+nz2i^+0j^+−nxnx2+nz2k^.From this, we can obtaint⟶2 based on n⟶ and t⟶1 so that these three vectors constitute a special orthonormal basis. The real velocity field is thus decomposed into one normal and two tangential velocity fields vn=v⟶⋅n⟶, vt1=v⟶⋅t⟶1 and vt2=v⟶⋅t⟶2. Both tangential velocities of the real area are then extended to the ghost area as vt1′ and vt2′ in the same fashion. Total velocity of the ghost fluid is thus a combination of real normal velocity and two ghost tangential velocities as v⟶ghost=vnn⟶+vt1′t⟶1+vt2′t⟶2. Extrapolation of density and tangential velocity in 2D and 3D implementation of DGFM is illustrated in Figure 3.Figure 3
Extrapolation of density and tangential velocity in 2D and 3D implementation of DGFM.An alternative way of implementing the DGFM in a 3D simulation is achieved by applying a basis-free projection method [6]. From the real velocity field v⟶=ui⟶+vj^+wk^, the velocity of the real area is extended to the ghost area as v⟶′=u′i⟶+v′j^+w′k^ by iterating equation (9). The real velocity field v⟶ is used to obtain real normal velocity vnn⟶=v⟶⋅n⟶n⟶ while the extrapolated velocity field v⟶′ is used to obtain ghost tangential velocity v⟶′−vn′n⟶=v⟶′−v⟶′⋅n⟶n⟶, respectively. Total ghost fluid velocity is thus the sum of real normal velocity and ghost tangential velocity as v⟶ghost=vnn⟶+v⟶′−vn′n⟶[6].In the extrapolation from real gas to ghost gas, the density at elementi−1 in the real gas is extrapolated to both element i in the real gas and element i+1,i+2,i+3, …, in the ghost gas if a 1D simulation is considered. Unlike the original GFM, the DGFM does not require an entropy calculation so the use of JWL EOS for the simulation of explosive gas does not cause any difficulty. In the extension of DGFM to multidimensions, complicated shock analysis in the MGFM and other gas dynamic-based methods is avoided [9, 11, 34]. The combination of the RKDG method and the DGFM is called the RKDG-DGFM. Assessment of the RKDG-DGFM is made by performing numerical simulations of multiple cases and comparing the results with analytic solutions, experiments, and numerical simulations using different schemes.


## 2.1. Governing Equations

The Euler equations in the conservative form are used as the fluid governing equations. The Arbitrary Lagrangian–Eulerian (ALE) form of governing equations are listed in the following equations.(1)∂ρ∂t+∇⋅ρm⟶=0,(2)∂ρv⟶∂t+∇.ρv⟶⊗m⟶+pI==0,(3)∂ρE∂t+∇⋅ρEm⟶+pv⟶=0,where the terms ρ,v⟶,m⟶=v⟶−v⟶grid,I=, and E are density, velocity, convection velocity, identity tensor, and total energy per unit mass, respectively. There are source terms on the right-hand side of equations (1)–(3) when a 1D cylindrical, 1D spherical, or 2D axisymmetric analysis is being performed. There are also source terms if body forces exist in the model configuration. The difference between velocities of fluid flow, v⟶, and velocities of the moving grid, v⟶grid, is the convection velocity, m⟶. The convection velocity equals the velocity of fluid flow in a pure Eulerian description because the grid is fixed in space. The convection velocity equals zero; however, in a pure Lagrangian description, because the grid moves together with the fluid flow, the ALE description stays between the Eulerian and the Lagrangian descriptions.This system is closed by the equation of state. The Johns–Wilkins–Lee (JWL) EOS is primarily used to model the explosive gas, but the ideal gas law is used for a comparison to some analytic results [2, 3, 6, 22]. The JWL EOS is(4)p=A1−ωηR1e−R1/η+B1−ωηR2e−R2/η+ωρeint,where A,B,R1,R2, and ω are experimentally determined coefficients and ρ, eint, and η are density, internal energy per unit mass, and the ratio of current density to the initial density, respectively. Liu demonstrated that JWL EOS can predict peak pressures approximating that of detonation experiments, while the ideal gas law can predict good results only beyond the region of 8 ∼ 10 times the charge radius [1, 3]. Despite this, the ideal gas law is often used in analytic solutions because of its common use in gas dynamics [23, 24] and its commonality with the stiffened gas EOS and Tait EOS used for water.The stiffened gas EOS is used to model the thermodynamic properties of water:(5)p=γB−1ρeint−γBB,where γB and B are experimentally determined coefficients. We have also used the Tait EOS:(6)p=N−1ρeint−NB−A,where N,B, and A are experimentally determined coefficients. The Grüneisen EOS is a more general EOS, capable of modeling a larger class of fluids. When the density change is small, the Grünesien EOS reduces to the stiffened gas EOS [22, 25].


## 2.2. Discretization

The Runge–Kutta Discontinuous Galerkin (RKDG) method is chosen as the numerical scheme to solve equations (1)–(3). The Discontinuous Galerkin (DG) method [26–28] discretizes the equations in space and the Runge–Kutta (RK) method discretizes in time. P2 Legendre polynomials are chosen as the basis functions in the DG method. This makes the spatial discretization formally third order accurate [14]. A three-step RK method is utilized whose formal order of accuracy is compatible with that of the spatial discretization [14, 28]. The DG method combines advantages of the Finite Element Method (FEM) and the Finite Volume Method (FVM). A local weak form provides weak connectivity between elements through the boundary integration of flux terms similar to the FVM. This relaxes the element-by-element continuity requirement. The remaining terms are treated using the standard Galerkin method. The p2 Legendre polynomial basis functions provide high-order accuracy and diagonalize the mass matrix without lumping. Simple extension of numerical flux functions allows variable discontinuities. Decoupled time derivatives between elements enable straightforward implementation of time marching algorithms, such as the explicit RK method. Details of one- and two-dimensional DG discretizations are provided in [14]. The three-dimensional DG discretization, verification, and validation are documented in a preceding publication (Si and Brown, 2021).In summary, the RKDG method as implemented here provides high accuracy and weak connectivity between elements and is easily extendable to multidimensions. Because of these characteristics, the RKDG method is competent in compressible flow analysis especially with strong shocks. The high accuracy of the RKDG method helps to produce less-diffusive numerical results compared to low-order approximation methods and helps to minimize inherent numerical diffusion, which is a source of error in Eulerian fluid methods. The order-of-accuracy verification of the 1D solver using a smooth configuration with sinusoidal initial conditions is described in [14]. A good match between observed and formal order of accuracy is indicated.


## 2.3. Multiphase and the Direct Ghost Fluid Method (DGFM)

To numerically treat the material interface using the GFM, the motion of the interface must be tracked at each time step during the computation [8, 29]. The GFM as implemented by Osher and Fedkiw used the Level Set Method (LSM) to perform this task. The LSM provides the geometry of the interface such as its curvature and the normal vector at different locations [29]. Because of its Eulerian nature, the implementation of the LSM does not increase geometric difficulty in multidimensions [29, 30]. The material interface is embedded in a spatially fixed grid. The transport of the LS equations is performed in an Eulerian manner [29, 30]. The transport and reinitialization equations are(7)∂ϕ∂t+v⟶⋅∇ϕ=0,(8)∂ϕ∂τ+Sϕ0∇ϕ−1=0,where ϕ is the level set (LS) function and Sϕ is a Heaviside function. Equation (7) implies that the interface is transported by the flow velocity. The function ϕ is the signed distance function ϕ=±d. It is set in the initial condition as the shortest normal distance (d) from any point in the domain to the material interface. In the LS function, a positive distance ϕ=+d indicates that a point resides in the water, a negative distance ϕ=−d indicates that a point resides in the explosive gas, and a zero distance ϕ=0 indicates that a point resides on the interface [29, 30]. The time-varying shape of the material interface is therefore evolved through the numerical solution of equation (7). Since the LS function ϕ must be transported by a real nonuniform flow field, the LS function may become distorted if not frequently reinitialized to keep the LS function real [29]. A nonuniform flow field generates noise that is not appropriate because the numerical results with distorted LS function may become unreasonable at subsequent time steps [29, 31]. The Hamilton–Jacobi WENO scheme and three-step RK method are used to numerically solve this transport equation and reinitialize the distorted LS function during the computation [32].In the GFM, elements possess unique conservative properties (i.e., mass, momentum and energy in both the real fluid and the ghost fluid [29, 30]). Once the conservative properties in the ghost fluid are determined, the two-fluid flows are solved separately using the RKDG method. In the solution of the fluid, there is no spatial concern for the interface. Solution of the fluid equations using the RKDG method provides the flow velocity v⟶. The LS function ϕ is transported in the subsequent time step using the computed flow velocity with the output of the fluid calculation in the last time step. Since pressure and normal velocity are continuous at the interface except for the initial step when pressure is discontinuous across the interface, real pressure and normal velocity are copied directly into ghost values in a simple node-to-node manner [21, 30] as described in Figure 2. In the original GFM, based on the isentropic assumption, entropy is copied by a one-sided constant extrapolation across the interface, which then allows the density to be calculated indirectly.Figure 2
Extrapolation of density in 1D implementation of DGFM for explosive gas and water flow simulations [6].Figure2 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].In the DGFM implementation, density is used instead of the entropy as in the original GFM [29] to avoid overshoot as shown in Figure 1. Density and tangential velocity is directly extrapolated from the real water to the ghost water (i.e., the area occupied by real gas) rather than using an indirect variable extrapolation from entropy. The density and tangential velocity are extrapolated using the isobaric fix technique from the real gas to the ghost gas (i.e., the area occupied by real water) minimizing the overheating error [29]. The isobaric fix technique is often included in the implementation of the GFM to minimize nonphysical temperature and density overshoot when a compressible flow interacts with a stiff medium [7, 8]. To distinguish this procedure from the indirect variable extrapolation in the original GFM, this direct extrapolation method is called the Direct Ghost Fluid Method. Diffusion of density is decreased by the DGFM across the material interface between the explosive gas and the water. The spurious pressure oscillation at the interface is therefore significantly reduced. Extrapolation of density in the 1D implementation of DGFM is illustrated in Figure 2. The density is extended from the real fluid into the ghost fluid by numerically solving the following advection equation:(9)∂ρ∂τ±n⟶⋅∇ρ=0,where n⟶=∇ϕ/∇ϕ is obtained from the LS function [29]. A few iterations, say 25, are required to be performed. The positive sign is used when density and tangential velocity are extended from the real gas (ϕ<0) to the ghost gas (ϕ>0), while the negative sign is used when they are extended from real water (ϕ>0) to ghost water (ϕ<0) [6].It is worth noting that there is an important detail in the implementation of the variable extrapolation by solving equation (9). In a 2D simulation, it is straightforward that the tangential vector is obtained by rotating the normal vector so that t⟶=−nyi^+nxj^. The real velocity field is thus decomposed into normal and tangential velocity fields vn=v⟶⋅n⟶ and vt=v⟶⋅t⟶. Tangential velocity of the real area is then extended to the ghost area as vt′ by iterating equation (9) so that total velocity of the ghost fluid is a combination of real normal velocity and ghost tangential velocity as v⟶ghost=vnn⟶+vt′t⟶. In a 3D simulation, however, there exists a pair of orthonormal bases to span the tangential plane that is perpendicular to the normal vector n⟶; for example,(10)t⟶1=nznx2+nz2i^+0j^+−nxnx2+nz2k^.From this, we can obtaint⟶2 based on n⟶ and t⟶1 so that these three vectors constitute a special orthonormal basis. The real velocity field is thus decomposed into one normal and two tangential velocity fields vn=v⟶⋅n⟶, vt1=v⟶⋅t⟶1 and vt2=v⟶⋅t⟶2. Both tangential velocities of the real area are then extended to the ghost area as vt1′ and vt2′ in the same fashion. Total velocity of the ghost fluid is thus a combination of real normal velocity and two ghost tangential velocities as v⟶ghost=vnn⟶+vt1′t⟶1+vt2′t⟶2. Extrapolation of density and tangential velocity in 2D and 3D implementation of DGFM is illustrated in Figure 3.Figure 3
Extrapolation of density and tangential velocity in 2D and 3D implementation of DGFM.An alternative way of implementing the DGFM in a 3D simulation is achieved by applying a basis-free projection method [6]. From the real velocity field v⟶=ui⟶+vj^+wk^, the velocity of the real area is extended to the ghost area as v⟶′=u′i⟶+v′j^+w′k^ by iterating equation (9). The real velocity field v⟶ is used to obtain real normal velocity vnn⟶=v⟶⋅n⟶n⟶ while the extrapolated velocity field v⟶′ is used to obtain ghost tangential velocity v⟶′−vn′n⟶=v⟶′−v⟶′⋅n⟶n⟶, respectively. Total ghost fluid velocity is thus the sum of real normal velocity and ghost tangential velocity as v⟶ghost=vnn⟶+v⟶′−vn′n⟶[6].In the extrapolation from real gas to ghost gas, the density at elementi−1 in the real gas is extrapolated to both element i in the real gas and element i+1,i+2,i+3, …, in the ghost gas if a 1D simulation is considered. Unlike the original GFM, the DGFM does not require an entropy calculation so the use of JWL EOS for the simulation of explosive gas does not cause any difficulty. In the extension of DGFM to multidimensions, complicated shock analysis in the MGFM and other gas dynamic-based methods is avoided [9, 11, 34]. The combination of the RKDG method and the DGFM is called the RKDG-DGFM. Assessment of the RKDG-DGFM is made by performing numerical simulations of multiple cases and comparing the results with analytic solutions, experiments, and numerical simulations using different schemes.


## 3. Assessment

This section describes the application and assessment of the RKDG-DGFM solvers. RKDG-DGFM results are compared with analytic solutions and other reference results [14]. Cases that are simulated using the developed CFD solver are summarized in Table 1.Table 1
Cases simulated by RKDG-DGFM.
Case no.Description11D shock tube problem21D spherical problem32D explosion problem43D explosion problem5–6Dynamics of explosive gas bubble7Explosion inside tube made of rigid wall8Explosion inside tube made of deformable wall


### 3.1. One-Dimensional Cases

Several 1D Cartesian cases are simulated. Analytic solutions are obtained using Riemann solutions [35] if they exist. RKDG calculations use CFL = 0.2 and 200 uniform/nonuniform elements in one direction. The material interface is identified by the LSM and treated by the DGFM. In our comparisons to the analytic solutions, the Tait EOS with N=7.15, B=3.31×109dyne/cm2 and ρ0=1.0g/cm3 is used to model the water. The ideal gas law with γ=1.4 is used to model the explosive gas. Two of these cases are chosen to be presented in this paper with their results illustrated.Case 1.
This case considers a gas-water shock tube problem previously solved by Qiu et al. [36]. A shock tube problem can be generally understood as a 1D tube initially filled by fluid of different states, that is, high- and low-pressure regions, as illustrated in Figure 4, and then the boundary is lifted instantaneously so that the fluids of different states interact. More details regarding the physics of Case 1 are illustrated in Figure 5 and results are shown in Figure 6.Figure 4
A typical shock tube problem.Figure 5
Initial conditions of Case1.Figure 6
RKDG-DGFM results of Case1 at 0.16×10−3 seconds.Case 2.
results are representative of the other 1D planar cases. This example is explored using the 1D two-fluid RKDG-DGFM, which treats the moving interfaces conservatively. The initial conditions are shown in Figure5. The 0,100 cm domain is discretized with 200 uniform mesh. The 1D RKDG-DGFM simulation is run at 0.16×10−3 seconds. A CFL equal to 2 is used in the simulation. The material interface is initially located at the center of the domain. It is expected that the rarefaction fan and shock wave do not reach their respective boundaries within this time, so the type of boundary conditions does not affect the solution in a hyperbolic system.
Qiu et al.’s method corrects the cell values near the interface by analytically computing local multimedium Riemann problems at the interface [36]. This is similar to the Godunov-type Prediction Correction Method [4], which requires local Riemann solutions to correct the affected cell values in the region near the interface. A good summary on the Godunov-type Prediction Correction Method is found in [4]. The RKDG-DGFM results, the analytic solutions, and Qiu et al.’s results are compared in Figure 6. At 0.16×10−3 second, the shock front and material interface are resolved at 86cm and 55cm, respectively. The RKDG-DGFM resolves the shock front with only 3 elements. Small density oscillations are observed in the region near the interface, but no pressure oscillations occur. Compared to Qiu et al.’s results taken from [36], the RKDG-DGFM provides less-diffusive and oscillatory density profiles.Case 3.
This case considers a spherical UNDEX flow that has been extensively used in method developments [4, 9, 11, 37, 38]. It was first studied by Flore and Holt using the Random Choice Method (RCM) [38]. Cooke and Chen employed the subcell resolution method for this example [11]. Wardlaw and Mair employed an ALE-based interface method [37] to solve this problem. The JWL EOS is used to model the explosive gas. The initial conditions are given in Figure 7. This case is initialized by a TNT-produced gas bubble with radius r0=16cm surrounded by water at rest.
The domain0,250cm is discretized with 200 nonuniform elements and the RKDG-DGFM simulation is run at 0.5×10−3 seconds and 1.0×10−3 seconds, respectively. A CFL equal to 2 is used in the simulation. A symmetry condition is imposed on the left while an outflow condition is imposed on the right. Figure 8 compares the RKDG-DGFM results with Wardlaw’s results. The pressure and density profiles are provided.
Figure8 shows that the two methods predict very similar pressure profiles at both final times. At 0.5×10−3 seconds, the material interface and shock front are located at 41cm and 120cm, respectively. The shock peak pressure is 1.9×109dyne/cm2. At 1.0×103 seconds, the material interface and shock front are located at 52cm and 205cm, respectively. The peak pressure is 0.75×109dyne/cm2. The locations of the material interface and shock fronts and the peak pressures compare well to those in [37].
Additional cases were evaluated to assess the applicability of the RKDG-DGFM for simple high-pressure gas and water flows, which are similar to TNT-produced gas and water flows in UNDEX events. These were presented in [14]. The high-pressure gas produced strong shocks traveling into the water and simultaneously rarefaction fans moving toward the origin. The RKDG-DGFM captured the shock fronts and material interfaces at the correct locations and magnitudes. Some density diffusion was observed in the region near the material interface but was less than the comparison methods. No spurious pressure oscillations causing failures occurred within the domains. Compared to the original GFM and the MGFM, the RKDG-DGFM did not increase difficulty or complexity in the implementation and computation of the JWL EOS for modeling the explosive gas. Since the RKDG-DGFM is based on the Eulerian description, no special concern for mesh adjustment is required as in an ALE-based interface method [19].Figure 7
Initial conditions of Case2.Figure 8
RKDG-DGFM results of Case2 at 0.5×10−3 and 1.0×10−3 seconds.


### 3.2. Multidimensional Cases

An alternative method to assess the effectiveness of the RKDG-DGFM in multidimensions is necessary because analytic solutions are generally unavailable [35]. In this paper, results on a fine mesh of the same level produced by 2D and 1D cylindrical solvers are used as a baseline for comparison. Results on a fine mesh of same level produced by 3D, 2D axisymmetric and 1D spherical solvers are also used as a baseline for comparison.Case 4.
considers a cylindrical explosive gas and water flow modeled by JWL and Tait EOS, respectively. The initial conditions are shown in Figure9. The 2D domain −100,100×−100,100cm is discretized with 200×200 uniform elements.
The center of the gas cylinder with radiusr0=50cm is located at the origin. The 2D and 1D cylindrical RKDG-DGFM simulations are both run at 0.1×10−3 seconds. It is expected that the shock front does not reach the boundary, so the boundary shape and type do not affect the disturbed region by the end of this time range in this hyperbolic problem. Numerical results are compared in Figure 10, which indicates a good match between 2D and 1D cylindrical solutions.
Contours of density and pressure are plotted in Figure11. Results are consistent with the xy-plots of the 2D results in Figure 10. No contour distortion occurs and symmetry is good.Figure 9
Initial conditions of Case3.Figure 10
RKDG-DGFM results of Case3 at 0.1×10−3 seconds.Figure 11
Density and pressure contours of Case3 at 0.1×10−3 seconds.Case 5.
explores an UNDEX scenario with its configuration entirely symmetric about the vertical axis [23, 35, 39]. The JWL EOS and Tait EOS model the explosive gas and surrounding flow of water, respectively. The initial conditions are illustrated in Figure 12. The 3D domain of −200,200×−200,200×−200,200cm is discretized with 200×200×200 uniform elements.
The center of the gas sphere with radiusr0=100cm is located at the origin. The 3D, 2D axisymmetric and 1D spherical RKDG-DGFM simulations are all run at 0.3×10−3 seconds. Again, it is expected that the shock front does not reach the boundary by the end of this time, so the boundary shape and type do not affect the disturbed region in this hyperbolic problem. Numerical results are compared in Figure 13, which indicates a good match among 3D, 2D axisymmetric and 1D spherical solutions.
Contours of density and pressure produced by 2D axisymmetric and 3D solvers are plotted in Figures14 and 15, respectively. The multidimensional results are consistent with their respective xy-plots in Figure 13. No contour distortion occurs and symmetry is good.Figure 12
Initial conditions of Case4.Figure 13
RKDG-DGFM results of Case4 at 0.3×10−3 seconds.Figure 14
2D axisymmetric density and pressure contours of Case4 at 0.3×10−3 seconds.Figure 15
3D density and pressure contours of Case4 at 0.3×10−3 seconds.


### 3.3. Practical UNDEX Applications

The performance of RKDG-DGFM solvers assessed in Cases1–4 demonstrates the potential applicability of this numerical framework for UNDEX simulation. In this section, we apply the RKDG-DGFM to more practical UNDEX applications, which include the gas bubble, shock-bubble interaction, and cavitation. Experimental data and reference results from previous works are used to assess the RKDG-DGFM results.


#### 3.3.1. UNDEX Explosion Bubble

Cases5 and 6 study the behavior of a spherical gas bubble, which is an important phenomenon in UNDEX events. As illustrated in Figure 16, following the detonation of an underwater explosion (UNDEX), the detonation wave proceeds quickly to the boundary of the explosive charge converting it (almost instantaneously) into a sphere or cylinder of high-pressure gas. A shock wave travels outward from this boundary into the water ahead of the expanding gas bubble much as we have seen in the 1D and 2D shock cases. As the bubble expands and the gas pressure inside the bubble drops to the hydrostatic level, the outward expansion slows until the expansion ceases after some inertial overshoot past pressure equilibrium. The now higher pressure in the surrounding water collapses the bubble back towards its original diameter, again with overshoot. As the bubble collapses to a minimum radius, it produces another outward pressure pulse. Several cycles follow at reduced strength due to the dissipation of energy [1, 40, 41].Figure 16
UNDEX gas dynamics [40].Case 6.
uses a spherical bubble to model the open-water UNDEX event. Since the explosion occurs deeply in the water and the duration of the simulation is relatively not very long, the body force term (gravity in this problem) is neglected. We therefore assume spherical symmetry and model the fluid flow using 1D spherical RKDG-DGFM. JWL and Stiffened EOSs are used to simulate the TNT-produced gas and the water, respectively. Wardlaw modeled this problem using an ALE-based method with a highly resolved mesh [34, 37] and his results are used for comparison.
Wardlaw’s case considers an UNDEX with28kg of TNT explosive surrounded by water at a depth of 178m. The TNT charge is replaced by a gas bubble with equivalent initial volume and initial internal energy of the original High Explosive (HE). The initial bubble radius r0 is determined by r0=3V0/4π1/3, where V0 is the initial volume of the explosive charge. The 1D spherical domain 0,10000cm is discretized with 2000 nonuniform elements with a refinement towards the charge center. The center of the gas bubble with r0=16cm is located at the origin. The initial conditions are specified in Figure 17. The 1D spherical RKDG-DGFM simulation is run at 0.15 seconds. The same boundary conditions as in Case 2 are imposed on the two ends, respectively.
Figure18 shows the calculated bubble radius-time history and the pressure at 121cm from the origin. The bubble radius is estimated using the sign change of the LS function ϕ. The maximum bubble radius is approximately 220cm at 0.07 seconds. The minimum bubble radius is approximately 30cm at 0.13 second. These results compare well with Wardlaw’s results. Table 2 lists Wardlaw’s results, the RKDG-DGFM results, and their difference. The small difference may be because Wardlaw identifies the material interface in a Lagrangian manner while the RKDG-DGFM uses the sign change of the LS function ϕ.Figure 17
Initial conditions of Case5.Figure 18
RKDG-DGFM results of Case5 at 0.15 seconds.Table 2
Comparison between simulations performed by RKDG-DGFM solver and Wardlaw.
WardlawRKDG-DGFMDifference (%)Max. bubble radius (cm)220221+0.45Bubble period (sec.)0.1320.130−1.50Peak pressure (dyne/cm2)1.9E91.9E90.0Case 7.
considers experiments conducted by Swift and Decius in 1950 [42]. They conducted measurements on the oscillations and maximum radii of gas bubbles from underwater explosions [42]. These experiments provide a useful bubble radius-time history for a 300g charge detonated at a depth of 91.46m. The experiments continue up to the second maximum radius. The initial conditions are shown in Figure 19. Brett also investigated this case using a nonlinear explicit FE-DYNA2D code [43] and Miller simulated this case using an OpenFOAM solver developed based on a cell-centered finite volume method where isothermal equations of state are added and thus solution of conservation of energy is skipped [44].
The domain0,9146cm is discretized with 2000 nonuniform elements. Compared to Case 5, a smaller initial radius r0=3.529cm corresponding to a 300g TNT charge is used. The 1D spherical RKDG-DGFM simulation is run at 0.045 seconds. The RKDG-DGFM results are compared with the experimental data and DYNA2D results in Figure 20. Both the RKDG-DGFM and DYNA2D predict the experimental data well for the first bubble pulse but deviate during the second pulse. Both simulations fail to predict enough energy loss that actually happens prior to and in the second bubble pulse. Unmodeled mechanisms of energy loss include viscous effects and thus turbulence and heat loss.Figure 19
Initial conditions of Case6.Figure 20
RKDG-DGFM results of Case6 at 0.045 seconds.


#### 3.3.2. Wall Interaction and Cavitation

In an UNDEX event with surrounding wall structures to interact with the fluid flow, the transient shock pressure loading on the structure causes an inward acceleration of the structure, which continues until the structure begins to move faster than the surrounding fluid. The fluid next to the structure is in tension and a local cavitation region develops at the fluid-structure interface. After reaching its maximum inward velocity, the structure begins to slow down and reverse direction or rebound. The local cavitation region closes and reloads the structure [1, 41, 45]. The bubble/shock cavitation mechanism associated with wall interaction is depicted in Figure 21. This section explores the performance of the RKDG-DGFM framework in its simulation of near-field UNDEX where the gas-water-cavitation-structure interact.Figure 21
Wall cavitation mechanism [46].Brett performed small-scale experiments showing that cavitation reloading and bubble effects in near-field UNDEX are as important as the primary shock impact [47]. To assess the cavitation formation, collapse, and reloading, Wardlaw suggested a pressure cutoff model in which cavitation is assumed to occur when pressure of the water falls below the vaporing pressure [46].(11)p=MAXppredict,plimit.In the cutoff cavitation model, the low-pressure limitplimit is usually 0.05 bars [47]. Brett’s research indicated that this limit hardly has any impact on the simulation results if the value is positive and much smaller than 1 bar [46]. For ease of implementation, we use this pressure cutoff model.Case 8.
investigates an internal explosion inside a water-filled tube subjected to an explosion of3g PETN internally. This experiment was performed by Sandusky et al. [48, 49]. The 2.8g PETN plus 0.2g detonator is located at the center of this tube that is made of aluminum. A thin plastic sheet is used to seal the bottom of the tube. The top of the tube is left open. A pressure gauge is installed at the inner center of the tube wall to measure the pressure-time history. Using results from the experiment, Sandusky et al. assessed the applicability of a coupled Eulerian–Lagrangian DYSMAS code developed by the Naval Surface Warfare Center (NSWC) Indian Division [48]. This problem is also simulated by Chambers [49] and Wardlaw [46] using DYSMAS and GEMINI-DYNA_N codes. A pressure cutoff model is used in both codes to predict the formation and collapse of cavitation [46, 49].
The experimental arrangement and its corresponding initial conditions in the mathematical model are shown in Figure22. Both 3D and 2D axisymmetric simulations can be performed since the configuration is axisymmetric about the z-axis and thus the fluid computation needs to only treat the sweeping plane. In the 2D axisymmetric simulation, the domain of 0,4.415×−8.9,8.9cm is discretized using 40×160 uniform elements. In the 3D simulation, the rectangular domain of −6,6×−8.9,8.9×−6,6cm is discretized using 120×178×120 uniform elements. The embedded boundary method is used to treat the wall boundary condition. Details of the embedded boundary method will be addressed in a subsequent paper. The JWL EOS with PETN parameters is used to model the thermodynamics of the explosive gas. The Tait EOS is used to model the thermodynamics of the water [9, 11, 46, 48].
Figure22 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].Figure 22
Experimental arrangement and its modeling. (a) Experimental arrangements. (b) Initial conditions.
(a)(b)case 9.
first assumes that the wall is rigid instead of deformable. The explosion occurs inside this rigid wall. In the 2D axisymmetric simulation, a symmetry condition is applied along thez-axis. Both the top and bottom are treated using a nonreflecting boundary condition (NRBC) to minimize the interference of reflection waves from these boundaries. The initial conditions of 2D axisymmetric and 3D simulations are set up and shown in Figure 23. The RKDG-DGFM simulation is run at 1.5×10−4 seconds. The time history of pressure at the inner center of the rigid wall is shown in Figure 24. Pressure contours at various times are plotted in Figure 25 and Figure 26 corresponding to 2D axisymmetric and 3D simulations, respectively.
Figure25 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].
The peak pressures at the inner center of the wall in the 2D simulation when the shock wave first arrives and later reloads after the cavitation collapse, that is, at1.5×10−5 and 1.2×10−4 seconds, are predicted to be 6.6×108Pa and 3.2×108Pa. These values are close to the reference results taken from [46]. The pressure when the cavitation collapses and reloads the rigid wall predicted by 3D simulation deviates from the 2D and DYSMAS simulations because the geometry of wall boundary is not perfectly conformed by the 3D fluid mesh in the embedded boundary method. The domain in the 3D fluid has flat sides between nodes at the fluid-wall boundary, while in reality it is curved. The embedded structural wall penetrates the fluid mesh and it penetrates at different angles along the perimeter. This leads to an asymmetry along the perimeter about the center axis which can be seen in the pressure contour at 1.25×10−4 seconds in Figure 26. This problem is alleviated as the fluid mesh is refined. The pressure-time histories produced by simulations using DYSMAS and RDKG-DGFM show excellent agreement for most time.
Pressure contours indicate that the shock wave, reflected wave, gas bubble and the interactions between them are dominant before the cavitation cutoff at around4×10−5 seconds. The shock wave loads and reflects off the rigid wall and interacts with the explosive gas bubble at around 3×10−5 seconds. An expansion wave is created through the interaction between the reflected wave and the explosive gas bubble. This expansion wave travels back towards the rigid wall and significantly reduces the pressure on the wall. Water close to the rigid wall is cavitated due to this low-pressure region. This local cavitation is formed at the inner center of the rigid wall at approximately 4×10−5 seconds. Pressure of the surrounding water is still high and eventually rushes into the cavitation zone, collapsing it at around 1.2×10−4 seconds and the tube wall is reloaded. These results indicate that the shock-bubble interaction is the main mechanism of formation and collapse of cavitation in the rigid wall case. A successful implementation of DGFM to treat the gas-water interface is crucial in the framework of the RKDG-DGFM fluid solver.
In Case 8, the rigid wall tube from Case 7 is replaced by a deformable wall tube made of AL5086 aluminum. Assuming that the structure density is constant and the entire process is adiabatic, the wall deformation can be modeled by the Lagrangian momentum equation:(12)ρDu⟶Dt=∇⋅σ=,where the stress tensor σ= is expanded using a stress-strain relationship for linear elastic materials. Typical properties of AL5086 are used in the calculation. The wall deformation is solved using a finite element method, while the ALE version of the fluid governing equations is solved using RKDG-DGFM with a moving wall condition. A partitioned FSI scheme is used to couple the solution of fluid and structure. Reference [14] provides more details for this FSI framework.
The RKDG-DGFM code is run at a final time of1.5×10−4 seconds. The pressure-time history at the inner center of the wall is shown in Figure 27. The pressure contours at various times are plotted in Figure 28. It is observed from the pressure contours that cavitation occurs at about 2.8×10−5 seconds in the water close to the wall, which is earlier than that in the rigid wall scenario. The reason is that this cavitation forms by a different mechanism that does not happen in the rigid wall case. The wall deforms because of the primary shock impact. The deformed wall generates an extra volume to be filled by the surrounding water and reduces the density and pressure of the water. This volume increase also makes the initial shock strength lower than that of the rigid wall. At around 3×10−5 seconds, local cavitation occurs midway between the tube and gas bubble due to the interaction between gas bubble and the reflected wave from the wall. Subsequently, the interaction generates an expansion wave traveling back toward the tube wall. These observations indicate that the first cavitation results from the lower wall pressure and the second cavitation results from the shock-bubble interaction as in the rigid wall tube case.
The pressure-time history from the RKDG-DGFM simulation matches closely with the experiment and Wardlaw’s work. The predicted peak pressure6.05×109dyne/cm2 is somewhat lower than the experimental value 6.25×109dyne/cm2. In the experiment, the first cavitation occurs at 3×10−5 seconds and then pressure reloads at 6×10−5 seconds. The RKDG-DGFM predicts that the first cavitation occurs at around the same moment, but the pressure reloads at an earlier time of 4.7×10−5 seconds. This is likely due to the constitutive model used in this work, which underpredicts the extent of wall deformation. The deformation of the wall in reality enters the plasticity region while the constitutive model used in this simulation is the hyperelastic model. As a result, it takes the surrounding water less time to rush into the cavitated region and collapse the cavitation. Unlike the DYSMAS calculation, a small cavitation region forms and collapses at the wall during between 4.7×10−5 and 5.8×10−5 seconds.
Figure28 are reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].
This 2D axisymmetric simulation uses an ALE method as the coupling approach between the fluid and the structural domain. This method can achieve high accuracy at the fluid-structure interface as the fluid mesh is always body-conforming. Special care, however, needs to be taken with the mesh-smoothing and remapping algorithm when the motion and deformation of the structural object are very large. In problems such as UNDEX, structural cracking and rupture might occur which cause a topological challenge that the ALE method cannot overcome. The embedded boundary method, therefore, will eventually be chosen as the coupling approach between the fluid and the structural domain. A 3D simulation of this case using the embedded boundary method will be performed and discussed in a subsequent paper.Figure 23
Initial conditions of Case 7 (left: 2D simulation; right: 3D simulation).Figure 24
Pressure-time history at the inner center of the rigid wall.Figure 25
Pressure contours inside the rigid wall tube (2D axisymmetric simulation).Figure 26
Pressure contours inside the rigid wall tube (3D simulation).Figure 27
Pressure-time history at the inner center of deformable wall.Figure 28
Pressure contours inside the deformable wall tube of half domain.


## 3.1. One-Dimensional Cases

Several 1D Cartesian cases are simulated. Analytic solutions are obtained using Riemann solutions [35] if they exist. RKDG calculations use CFL = 0.2 and 200 uniform/nonuniform elements in one direction. The material interface is identified by the LSM and treated by the DGFM. In our comparisons to the analytic solutions, the Tait EOS with N=7.15, B=3.31×109dyne/cm2 and ρ0=1.0g/cm3 is used to model the water. The ideal gas law with γ=1.4 is used to model the explosive gas. Two of these cases are chosen to be presented in this paper with their results illustrated.Case 1.
This case considers a gas-water shock tube problem previously solved by Qiu et al. [36]. A shock tube problem can be generally understood as a 1D tube initially filled by fluid of different states, that is, high- and low-pressure regions, as illustrated in Figure 4, and then the boundary is lifted instantaneously so that the fluids of different states interact. More details regarding the physics of Case 1 are illustrated in Figure 5 and results are shown in Figure 6.Figure 4
A typical shock tube problem.Figure 5
Initial conditions of Case1.Figure 6
RKDG-DGFM results of Case1 at 0.16×10−3 seconds.Case 2.
results are representative of the other 1D planar cases. This example is explored using the 1D two-fluid RKDG-DGFM, which treats the moving interfaces conservatively. The initial conditions are shown in Figure5. The 0,100 cm domain is discretized with 200 uniform mesh. The 1D RKDG-DGFM simulation is run at 0.16×10−3 seconds. A CFL equal to 2 is used in the simulation. The material interface is initially located at the center of the domain. It is expected that the rarefaction fan and shock wave do not reach their respective boundaries within this time, so the type of boundary conditions does not affect the solution in a hyperbolic system.
Qiu et al.’s method corrects the cell values near the interface by analytically computing local multimedium Riemann problems at the interface [36]. This is similar to the Godunov-type Prediction Correction Method [4], which requires local Riemann solutions to correct the affected cell values in the region near the interface. A good summary on the Godunov-type Prediction Correction Method is found in [4]. The RKDG-DGFM results, the analytic solutions, and Qiu et al.’s results are compared in Figure 6. At 0.16×10−3 second, the shock front and material interface are resolved at 86cm and 55cm, respectively. The RKDG-DGFM resolves the shock front with only 3 elements. Small density oscillations are observed in the region near the interface, but no pressure oscillations occur. Compared to Qiu et al.’s results taken from [36], the RKDG-DGFM provides less-diffusive and oscillatory density profiles.Case 3.
This case considers a spherical UNDEX flow that has been extensively used in method developments [4, 9, 11, 37, 38]. It was first studied by Flore and Holt using the Random Choice Method (RCM) [38]. Cooke and Chen employed the subcell resolution method for this example [11]. Wardlaw and Mair employed an ALE-based interface method [37] to solve this problem. The JWL EOS is used to model the explosive gas. The initial conditions are given in Figure 7. This case is initialized by a TNT-produced gas bubble with radius r0=16cm surrounded by water at rest.
The domain0,250cm is discretized with 200 nonuniform elements and the RKDG-DGFM simulation is run at 0.5×10−3 seconds and 1.0×10−3 seconds, respectively. A CFL equal to 2 is used in the simulation. A symmetry condition is imposed on the left while an outflow condition is imposed on the right. Figure 8 compares the RKDG-DGFM results with Wardlaw’s results. The pressure and density profiles are provided.
Figure8 shows that the two methods predict very similar pressure profiles at both final times. At 0.5×10−3 seconds, the material interface and shock front are located at 41cm and 120cm, respectively. The shock peak pressure is 1.9×109dyne/cm2. At 1.0×103 seconds, the material interface and shock front are located at 52cm and 205cm, respectively. The peak pressure is 0.75×109dyne/cm2. The locations of the material interface and shock fronts and the peak pressures compare well to those in [37].
Additional cases were evaluated to assess the applicability of the RKDG-DGFM for simple high-pressure gas and water flows, which are similar to TNT-produced gas and water flows in UNDEX events. These were presented in [14]. The high-pressure gas produced strong shocks traveling into the water and simultaneously rarefaction fans moving toward the origin. The RKDG-DGFM captured the shock fronts and material interfaces at the correct locations and magnitudes. Some density diffusion was observed in the region near the material interface but was less than the comparison methods. No spurious pressure oscillations causing failures occurred within the domains. Compared to the original GFM and the MGFM, the RKDG-DGFM did not increase difficulty or complexity in the implementation and computation of the JWL EOS for modeling the explosive gas. Since the RKDG-DGFM is based on the Eulerian description, no special concern for mesh adjustment is required as in an ALE-based interface method [19].Figure 7
Initial conditions of Case2.Figure 8
RKDG-DGFM results of Case2 at 0.5×10−3 and 1.0×10−3 seconds.


## 3.2. Multidimensional Cases

An alternative method to assess the effectiveness of the RKDG-DGFM in multidimensions is necessary because analytic solutions are generally unavailable [35]. In this paper, results on a fine mesh of the same level produced by 2D and 1D cylindrical solvers are used as a baseline for comparison. Results on a fine mesh of same level produced by 3D, 2D axisymmetric and 1D spherical solvers are also used as a baseline for comparison.Case 4.
considers a cylindrical explosive gas and water flow modeled by JWL and Tait EOS, respectively. The initial conditions are shown in Figure9. The 2D domain −100,100×−100,100cm is discretized with 200×200 uniform elements.
The center of the gas cylinder with radiusr0=50cm is located at the origin. The 2D and 1D cylindrical RKDG-DGFM simulations are both run at 0.1×10−3 seconds. It is expected that the shock front does not reach the boundary, so the boundary shape and type do not affect the disturbed region by the end of this time range in this hyperbolic problem. Numerical results are compared in Figure 10, which indicates a good match between 2D and 1D cylindrical solutions.
Contours of density and pressure are plotted in Figure11. Results are consistent with the xy-plots of the 2D results in Figure 10. No contour distortion occurs and symmetry is good.Figure 9
Initial conditions of Case3.Figure 10
RKDG-DGFM results of Case3 at 0.1×10−3 seconds.Figure 11
Density and pressure contours of Case3 at 0.1×10−3 seconds.Case 5.
explores an UNDEX scenario with its configuration entirely symmetric about the vertical axis [23, 35, 39]. The JWL EOS and Tait EOS model the explosive gas and surrounding flow of water, respectively. The initial conditions are illustrated in Figure 12. The 3D domain of −200,200×−200,200×−200,200cm is discretized with 200×200×200 uniform elements.
The center of the gas sphere with radiusr0=100cm is located at the origin. The 3D, 2D axisymmetric and 1D spherical RKDG-DGFM simulations are all run at 0.3×10−3 seconds. Again, it is expected that the shock front does not reach the boundary by the end of this time, so the boundary shape and type do not affect the disturbed region in this hyperbolic problem. Numerical results are compared in Figure 13, which indicates a good match among 3D, 2D axisymmetric and 1D spherical solutions.
Contours of density and pressure produced by 2D axisymmetric and 3D solvers are plotted in Figures14 and 15, respectively. The multidimensional results are consistent with their respective xy-plots in Figure 13. No contour distortion occurs and symmetry is good.Figure 12
Initial conditions of Case4.Figure 13
RKDG-DGFM results of Case4 at 0.3×10−3 seconds.Figure 14
2D axisymmetric density and pressure contours of Case4 at 0.3×10−3 seconds.Figure 15
3D density and pressure contours of Case4 at 0.3×10−3 seconds.


## 3.3. Practical UNDEX Applications

The performance of RKDG-DGFM solvers assessed in Cases1–4 demonstrates the potential applicability of this numerical framework for UNDEX simulation. In this section, we apply the RKDG-DGFM to more practical UNDEX applications, which include the gas bubble, shock-bubble interaction, and cavitation. Experimental data and reference results from previous works are used to assess the RKDG-DGFM results.


### 3.3.1. UNDEX Explosion Bubble

Cases5 and 6 study the behavior of a spherical gas bubble, which is an important phenomenon in UNDEX events. As illustrated in Figure 16, following the detonation of an underwater explosion (UNDEX), the detonation wave proceeds quickly to the boundary of the explosive charge converting it (almost instantaneously) into a sphere or cylinder of high-pressure gas. A shock wave travels outward from this boundary into the water ahead of the expanding gas bubble much as we have seen in the 1D and 2D shock cases. As the bubble expands and the gas pressure inside the bubble drops to the hydrostatic level, the outward expansion slows until the expansion ceases after some inertial overshoot past pressure equilibrium. The now higher pressure in the surrounding water collapses the bubble back towards its original diameter, again with overshoot. As the bubble collapses to a minimum radius, it produces another outward pressure pulse. Several cycles follow at reduced strength due to the dissipation of energy [1, 40, 41].Figure 16
UNDEX gas dynamics [40].Case 6.
uses a spherical bubble to model the open-water UNDEX event. Since the explosion occurs deeply in the water and the duration of the simulation is relatively not very long, the body force term (gravity in this problem) is neglected. We therefore assume spherical symmetry and model the fluid flow using 1D spherical RKDG-DGFM. JWL and Stiffened EOSs are used to simulate the TNT-produced gas and the water, respectively. Wardlaw modeled this problem using an ALE-based method with a highly resolved mesh [34, 37] and his results are used for comparison.
Wardlaw’s case considers an UNDEX with28kg of TNT explosive surrounded by water at a depth of 178m. The TNT charge is replaced by a gas bubble with equivalent initial volume and initial internal energy of the original High Explosive (HE). The initial bubble radius r0 is determined by r0=3V0/4π1/3, where V0 is the initial volume of the explosive charge. The 1D spherical domain 0,10000cm is discretized with 2000 nonuniform elements with a refinement towards the charge center. The center of the gas bubble with r0=16cm is located at the origin. The initial conditions are specified in Figure 17. The 1D spherical RKDG-DGFM simulation is run at 0.15 seconds. The same boundary conditions as in Case 2 are imposed on the two ends, respectively.
Figure18 shows the calculated bubble radius-time history and the pressure at 121cm from the origin. The bubble radius is estimated using the sign change of the LS function ϕ. The maximum bubble radius is approximately 220cm at 0.07 seconds. The minimum bubble radius is approximately 30cm at 0.13 second. These results compare well with Wardlaw’s results. Table 2 lists Wardlaw’s results, the RKDG-DGFM results, and their difference. The small difference may be because Wardlaw identifies the material interface in a Lagrangian manner while the RKDG-DGFM uses the sign change of the LS function ϕ.Figure 17
Initial conditions of Case5.Figure 18
RKDG-DGFM results of Case5 at 0.15 seconds.Table 2
Comparison between simulations performed by RKDG-DGFM solver and Wardlaw.
WardlawRKDG-DGFMDifference (%)Max. bubble radius (cm)220221+0.45Bubble period (sec.)0.1320.130−1.50Peak pressure (dyne/cm2)1.9E91.9E90.0Case 7.
considers experiments conducted by Swift and Decius in 1950 [42]. They conducted measurements on the oscillations and maximum radii of gas bubbles from underwater explosions [42]. These experiments provide a useful bubble radius-time history for a 300g charge detonated at a depth of 91.46m. The experiments continue up to the second maximum radius. The initial conditions are shown in Figure 19. Brett also investigated this case using a nonlinear explicit FE-DYNA2D code [43] and Miller simulated this case using an OpenFOAM solver developed based on a cell-centered finite volume method where isothermal equations of state are added and thus solution of conservation of energy is skipped [44].
The domain0,9146cm is discretized with 2000 nonuniform elements. Compared to Case 5, a smaller initial radius r0=3.529cm corresponding to a 300g TNT charge is used. The 1D spherical RKDG-DGFM simulation is run at 0.045 seconds. The RKDG-DGFM results are compared with the experimental data and DYNA2D results in Figure 20. Both the RKDG-DGFM and DYNA2D predict the experimental data well for the first bubble pulse but deviate during the second pulse. Both simulations fail to predict enough energy loss that actually happens prior to and in the second bubble pulse. Unmodeled mechanisms of energy loss include viscous effects and thus turbulence and heat loss.Figure 19
Initial conditions of Case6.Figure 20
RKDG-DGFM results of Case6 at 0.045 seconds.


### 3.3.2. Wall Interaction and Cavitation

In an UNDEX event with surrounding wall structures to interact with the fluid flow, the transient shock pressure loading on the structure causes an inward acceleration of the structure, which continues until the structure begins to move faster than the surrounding fluid. The fluid next to the structure is in tension and a local cavitation region develops at the fluid-structure interface. After reaching its maximum inward velocity, the structure begins to slow down and reverse direction or rebound. The local cavitation region closes and reloads the structure [1, 41, 45]. The bubble/shock cavitation mechanism associated with wall interaction is depicted in Figure 21. This section explores the performance of the RKDG-DGFM framework in its simulation of near-field UNDEX where the gas-water-cavitation-structure interact.Figure 21
Wall cavitation mechanism [46].Brett performed small-scale experiments showing that cavitation reloading and bubble effects in near-field UNDEX are as important as the primary shock impact [47]. To assess the cavitation formation, collapse, and reloading, Wardlaw suggested a pressure cutoff model in which cavitation is assumed to occur when pressure of the water falls below the vaporing pressure [46].(11)p=MAXppredict,plimit.In the cutoff cavitation model, the low-pressure limitplimit is usually 0.05 bars [47]. Brett’s research indicated that this limit hardly has any impact on the simulation results if the value is positive and much smaller than 1 bar [46]. For ease of implementation, we use this pressure cutoff model.Case 8.
investigates an internal explosion inside a water-filled tube subjected to an explosion of3g PETN internally. This experiment was performed by Sandusky et al. [48, 49]. The 2.8g PETN plus 0.2g detonator is located at the center of this tube that is made of aluminum. A thin plastic sheet is used to seal the bottom of the tube. The top of the tube is left open. A pressure gauge is installed at the inner center of the tube wall to measure the pressure-time history. Using results from the experiment, Sandusky et al. assessed the applicability of a coupled Eulerian–Lagrangian DYSMAS code developed by the Naval Surface Warfare Center (NSWC) Indian Division [48]. This problem is also simulated by Chambers [49] and Wardlaw [46] using DYSMAS and GEMINI-DYNA_N codes. A pressure cutoff model is used in both codes to predict the formation and collapse of cavitation [46, 49].
The experimental arrangement and its corresponding initial conditions in the mathematical model are shown in Figure22. Both 3D and 2D axisymmetric simulations can be performed since the configuration is axisymmetric about the z-axis and thus the fluid computation needs to only treat the sweeping plane. In the 2D axisymmetric simulation, the domain of 0,4.415×−8.9,8.9cm is discretized using 40×160 uniform elements. In the 3D simulation, the rectangular domain of −6,6×−8.9,8.9×−6,6cm is discretized using 120×178×120 uniform elements. The embedded boundary method is used to treat the wall boundary condition. Details of the embedded boundary method will be addressed in a subsequent paper. The JWL EOS with PETN parameters is used to model the thermodynamics of the explosive gas. The Tait EOS is used to model the thermodynamics of the water [9, 11, 46, 48].
Figure22 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].Figure 22
Experimental arrangement and its modeling. (a) Experimental arrangements. (b) Initial conditions.
(a)(b)case 9.
first assumes that the wall is rigid instead of deformable. The explosion occurs inside this rigid wall. In the 2D axisymmetric simulation, a symmetry condition is applied along thez-axis. Both the top and bottom are treated using a nonreflecting boundary condition (NRBC) to minimize the interference of reflection waves from these boundaries. The initial conditions of 2D axisymmetric and 3D simulations are set up and shown in Figure 23. The RKDG-DGFM simulation is run at 1.5×10−4 seconds. The time history of pressure at the inner center of the rigid wall is shown in Figure 24. Pressure contours at various times are plotted in Figure 25 and Figure 26 corresponding to 2D axisymmetric and 3D simulations, respectively.
Figure25 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].
The peak pressures at the inner center of the wall in the 2D simulation when the shock wave first arrives and later reloads after the cavitation collapse, that is, at1.5×10−5 and 1.2×10−4 seconds, are predicted to be 6.6×108Pa and 3.2×108Pa. These values are close to the reference results taken from [46]. The pressure when the cavitation collapses and reloads the rigid wall predicted by 3D simulation deviates from the 2D and DYSMAS simulations because the geometry of wall boundary is not perfectly conformed by the 3D fluid mesh in the embedded boundary method. The domain in the 3D fluid has flat sides between nodes at the fluid-wall boundary, while in reality it is curved. The embedded structural wall penetrates the fluid mesh and it penetrates at different angles along the perimeter. This leads to an asymmetry along the perimeter about the center axis which can be seen in the pressure contour at 1.25×10−4 seconds in Figure 26. This problem is alleviated as the fluid mesh is refined. The pressure-time histories produced by simulations using DYSMAS and RDKG-DGFM show excellent agreement for most time.
Pressure contours indicate that the shock wave, reflected wave, gas bubble and the interactions between them are dominant before the cavitation cutoff at around4×10−5 seconds. The shock wave loads and reflects off the rigid wall and interacts with the explosive gas bubble at around 3×10−5 seconds. An expansion wave is created through the interaction between the reflected wave and the explosive gas bubble. This expansion wave travels back towards the rigid wall and significantly reduces the pressure on the wall. Water close to the rigid wall is cavitated due to this low-pressure region. This local cavitation is formed at the inner center of the rigid wall at approximately 4×10−5 seconds. Pressure of the surrounding water is still high and eventually rushes into the cavitation zone, collapsing it at around 1.2×10−4 seconds and the tube wall is reloaded. These results indicate that the shock-bubble interaction is the main mechanism of formation and collapse of cavitation in the rigid wall case. A successful implementation of DGFM to treat the gas-water interface is crucial in the framework of the RKDG-DGFM fluid solver.
In Case 8, the rigid wall tube from Case 7 is replaced by a deformable wall tube made of AL5086 aluminum. Assuming that the structure density is constant and the entire process is adiabatic, the wall deformation can be modeled by the Lagrangian momentum equation:(12)ρDu⟶Dt=∇⋅σ=,where the stress tensor σ= is expanded using a stress-strain relationship for linear elastic materials. Typical properties of AL5086 are used in the calculation. The wall deformation is solved using a finite element method, while the ALE version of the fluid governing equations is solved using RKDG-DGFM with a moving wall condition. A partitioned FSI scheme is used to couple the solution of fluid and structure. Reference [14] provides more details for this FSI framework.
The RKDG-DGFM code is run at a final time of1.5×10−4 seconds. The pressure-time history at the inner center of the wall is shown in Figure 27. The pressure contours at various times are plotted in Figure 28. It is observed from the pressure contours that cavitation occurs at about 2.8×10−5 seconds in the water close to the wall, which is earlier than that in the rigid wall scenario. The reason is that this cavitation forms by a different mechanism that does not happen in the rigid wall case. The wall deforms because of the primary shock impact. The deformed wall generates an extra volume to be filled by the surrounding water and reduces the density and pressure of the water. This volume increase also makes the initial shock strength lower than that of the rigid wall. At around 3×10−5 seconds, local cavitation occurs midway between the tube and gas bubble due to the interaction between gas bubble and the reflected wave from the wall. Subsequently, the interaction generates an expansion wave traveling back toward the tube wall. These observations indicate that the first cavitation results from the lower wall pressure and the second cavitation results from the shock-bubble interaction as in the rigid wall tube case.
The pressure-time history from the RKDG-DGFM simulation matches closely with the experiment and Wardlaw’s work. The predicted peak pressure6.05×109dyne/cm2 is somewhat lower than the experimental value 6.25×109dyne/cm2. In the experiment, the first cavitation occurs at 3×10−5 seconds and then pressure reloads at 6×10−5 seconds. The RKDG-DGFM predicts that the first cavitation occurs at around the same moment, but the pressure reloads at an earlier time of 4.7×10−5 seconds. This is likely due to the constitutive model used in this work, which underpredicts the extent of wall deformation. The deformation of the wall in reality enters the plasticity region while the constitutive model used in this simulation is the hyperelastic model. As a result, it takes the surrounding water less time to rush into the cavitated region and collapse the cavitation. Unlike the DYSMAS calculation, a small cavitation region forms and collapses at the wall during between 4.7×10−5 and 5.8×10−5 seconds.
Figure28 are reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].
This 2D axisymmetric simulation uses an ALE method as the coupling approach between the fluid and the structural domain. This method can achieve high accuracy at the fluid-structure interface as the fluid mesh is always body-conforming. Special care, however, needs to be taken with the mesh-smoothing and remapping algorithm when the motion and deformation of the structural object are very large. In problems such as UNDEX, structural cracking and rupture might occur which cause a topological challenge that the ALE method cannot overcome. The embedded boundary method, therefore, will eventually be chosen as the coupling approach between the fluid and the structural domain. A 3D simulation of this case using the embedded boundary method will be performed and discussed in a subsequent paper.Figure 23
Initial conditions of Case 7 (left: 2D simulation; right: 3D simulation).Figure 24
Pressure-time history at the inner center of the rigid wall.Figure 25
Pressure contours inside the rigid wall tube (2D axisymmetric simulation).Figure 26
Pressure contours inside the rigid wall tube (3D simulation).Figure 27
Pressure-time history at the inner center of deformable wall.Figure 28
Pressure contours inside the deformable wall tube of half domain.


## 3.3.1. UNDEX Explosion Bubble

Cases5 and 6 study the behavior of a spherical gas bubble, which is an important phenomenon in UNDEX events. As illustrated in Figure 16, following the detonation of an underwater explosion (UNDEX), the detonation wave proceeds quickly to the boundary of the explosive charge converting it (almost instantaneously) into a sphere or cylinder of high-pressure gas. A shock wave travels outward from this boundary into the water ahead of the expanding gas bubble much as we have seen in the 1D and 2D shock cases. As the bubble expands and the gas pressure inside the bubble drops to the hydrostatic level, the outward expansion slows until the expansion ceases after some inertial overshoot past pressure equilibrium. The now higher pressure in the surrounding water collapses the bubble back towards its original diameter, again with overshoot. As the bubble collapses to a minimum radius, it produces another outward pressure pulse. Several cycles follow at reduced strength due to the dissipation of energy [1, 40, 41].Figure 16
UNDEX gas dynamics [40].Case 6.
uses a spherical bubble to model the open-water UNDEX event. Since the explosion occurs deeply in the water and the duration of the simulation is relatively not very long, the body force term (gravity in this problem) is neglected. We therefore assume spherical symmetry and model the fluid flow using 1D spherical RKDG-DGFM. JWL and Stiffened EOSs are used to simulate the TNT-produced gas and the water, respectively. Wardlaw modeled this problem using an ALE-based method with a highly resolved mesh [34, 37] and his results are used for comparison.
Wardlaw’s case considers an UNDEX with28kg of TNT explosive surrounded by water at a depth of 178m. The TNT charge is replaced by a gas bubble with equivalent initial volume and initial internal energy of the original High Explosive (HE). The initial bubble radius r0 is determined by r0=3V0/4π1/3, where V0 is the initial volume of the explosive charge. The 1D spherical domain 0,10000cm is discretized with 2000 nonuniform elements with a refinement towards the charge center. The center of the gas bubble with r0=16cm is located at the origin. The initial conditions are specified in Figure 17. The 1D spherical RKDG-DGFM simulation is run at 0.15 seconds. The same boundary conditions as in Case 2 are imposed on the two ends, respectively.
Figure18 shows the calculated bubble radius-time history and the pressure at 121cm from the origin. The bubble radius is estimated using the sign change of the LS function ϕ. The maximum bubble radius is approximately 220cm at 0.07 seconds. The minimum bubble radius is approximately 30cm at 0.13 second. These results compare well with Wardlaw’s results. Table 2 lists Wardlaw’s results, the RKDG-DGFM results, and their difference. The small difference may be because Wardlaw identifies the material interface in a Lagrangian manner while the RKDG-DGFM uses the sign change of the LS function ϕ.Figure 17
Initial conditions of Case5.Figure 18
RKDG-DGFM results of Case5 at 0.15 seconds.Table 2
Comparison between simulations performed by RKDG-DGFM solver and Wardlaw.
WardlawRKDG-DGFMDifference (%)Max. bubble radius (cm)220221+0.45Bubble period (sec.)0.1320.130−1.50Peak pressure (dyne/cm2)1.9E91.9E90.0Case 7.
considers experiments conducted by Swift and Decius in 1950 [42]. They conducted measurements on the oscillations and maximum radii of gas bubbles from underwater explosions [42]. These experiments provide a useful bubble radius-time history for a 300g charge detonated at a depth of 91.46m. The experiments continue up to the second maximum radius. The initial conditions are shown in Figure 19. Brett also investigated this case using a nonlinear explicit FE-DYNA2D code [43] and Miller simulated this case using an OpenFOAM solver developed based on a cell-centered finite volume method where isothermal equations of state are added and thus solution of conservation of energy is skipped [44].
The domain0,9146cm is discretized with 2000 nonuniform elements. Compared to Case 5, a smaller initial radius r0=3.529cm corresponding to a 300g TNT charge is used. The 1D spherical RKDG-DGFM simulation is run at 0.045 seconds. The RKDG-DGFM results are compared with the experimental data and DYNA2D results in Figure 20. Both the RKDG-DGFM and DYNA2D predict the experimental data well for the first bubble pulse but deviate during the second pulse. Both simulations fail to predict enough energy loss that actually happens prior to and in the second bubble pulse. Unmodeled mechanisms of energy loss include viscous effects and thus turbulence and heat loss.Figure 19
Initial conditions of Case6.Figure 20
RKDG-DGFM results of Case6 at 0.045 seconds.


## 3.3.2. Wall Interaction and Cavitation

In an UNDEX event with surrounding wall structures to interact with the fluid flow, the transient shock pressure loading on the structure causes an inward acceleration of the structure, which continues until the structure begins to move faster than the surrounding fluid. The fluid next to the structure is in tension and a local cavitation region develops at the fluid-structure interface. After reaching its maximum inward velocity, the structure begins to slow down and reverse direction or rebound. The local cavitation region closes and reloads the structure [1, 41, 45]. The bubble/shock cavitation mechanism associated with wall interaction is depicted in Figure 21. This section explores the performance of the RKDG-DGFM framework in its simulation of near-field UNDEX where the gas-water-cavitation-structure interact.Figure 21
Wall cavitation mechanism [46].Brett performed small-scale experiments showing that cavitation reloading and bubble effects in near-field UNDEX are as important as the primary shock impact [47]. To assess the cavitation formation, collapse, and reloading, Wardlaw suggested a pressure cutoff model in which cavitation is assumed to occur when pressure of the water falls below the vaporing pressure [46].(11)p=MAXppredict,plimit.In the cutoff cavitation model, the low-pressure limitplimit is usually 0.05 bars [47]. Brett’s research indicated that this limit hardly has any impact on the simulation results if the value is positive and much smaller than 1 bar [46]. For ease of implementation, we use this pressure cutoff model.Case 8.
investigates an internal explosion inside a water-filled tube subjected to an explosion of3g PETN internally. This experiment was performed by Sandusky et al. [48, 49]. The 2.8g PETN plus 0.2g detonator is located at the center of this tube that is made of aluminum. A thin plastic sheet is used to seal the bottom of the tube. The top of the tube is left open. A pressure gauge is installed at the inner center of the tube wall to measure the pressure-time history. Using results from the experiment, Sandusky et al. assessed the applicability of a coupled Eulerian–Lagrangian DYSMAS code developed by the Naval Surface Warfare Center (NSWC) Indian Division [48]. This problem is also simulated by Chambers [49] and Wardlaw [46] using DYSMAS and GEMINI-DYNA_N codes. A pressure cutoff model is used in both codes to predict the formation and collapse of cavitation [46, 49].
The experimental arrangement and its corresponding initial conditions in the mathematical model are shown in Figure22. Both 3D and 2D axisymmetric simulations can be performed since the configuration is axisymmetric about the z-axis and thus the fluid computation needs to only treat the sweeping plane. In the 2D axisymmetric simulation, the domain of 0,4.415×−8.9,8.9cm is discretized using 40×160 uniform elements. In the 3D simulation, the rectangular domain of −6,6×−8.9,8.9×−6,6cm is discretized using 120×178×120 uniform elements. The embedded boundary method is used to treat the wall boundary condition. Details of the embedded boundary method will be addressed in a subsequent paper. The JWL EOS with PETN parameters is used to model the thermodynamics of the explosive gas. The Tait EOS is used to model the thermodynamics of the water [9, 11, 46, 48].
Figure22 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].Figure 22
Experimental arrangement and its modeling. (a) Experimental arrangements. (b) Initial conditions.
(a)(b)case 9.
first assumes that the wall is rigid instead of deformable. The explosion occurs inside this rigid wall. In the 2D axisymmetric simulation, a symmetry condition is applied along thez-axis. Both the top and bottom are treated using a nonreflecting boundary condition (NRBC) to minimize the interference of reflection waves from these boundaries. The initial conditions of 2D axisymmetric and 3D simulations are set up and shown in Figure 23. The RKDG-DGFM simulation is run at 1.5×10−4 seconds. The time history of pressure at the inner center of the rigid wall is shown in Figure 24. Pressure contours at various times are plotted in Figure 25 and Figure 26 corresponding to 2D axisymmetric and 3D simulations, respectively.
Figure25 is reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].
The peak pressures at the inner center of the wall in the 2D simulation when the shock wave first arrives and later reloads after the cavitation collapse, that is, at1.5×10−5 and 1.2×10−4 seconds, are predicted to be 6.6×108Pa and 3.2×108Pa. These values are close to the reference results taken from [46]. The pressure when the cavitation collapses and reloads the rigid wall predicted by 3D simulation deviates from the 2D and DYSMAS simulations because the geometry of wall boundary is not perfectly conformed by the 3D fluid mesh in the embedded boundary method. The domain in the 3D fluid has flat sides between nodes at the fluid-wall boundary, while in reality it is curved. The embedded structural wall penetrates the fluid mesh and it penetrates at different angles along the perimeter. This leads to an asymmetry along the perimeter about the center axis which can be seen in the pressure contour at 1.25×10−4 seconds in Figure 26. This problem is alleviated as the fluid mesh is refined. The pressure-time histories produced by simulations using DYSMAS and RDKG-DGFM show excellent agreement for most time.
Pressure contours indicate that the shock wave, reflected wave, gas bubble and the interactions between them are dominant before the cavitation cutoff at around4×10−5 seconds. The shock wave loads and reflects off the rigid wall and interacts with the explosive gas bubble at around 3×10−5 seconds. An expansion wave is created through the interaction between the reflected wave and the explosive gas bubble. This expansion wave travels back towards the rigid wall and significantly reduces the pressure on the wall. Water close to the rigid wall is cavitated due to this low-pressure region. This local cavitation is formed at the inner center of the rigid wall at approximately 4×10−5 seconds. Pressure of the surrounding water is still high and eventually rushes into the cavitation zone, collapsing it at around 1.2×10−4 seconds and the tube wall is reloaded. These results indicate that the shock-bubble interaction is the main mechanism of formation and collapse of cavitation in the rigid wall case. A successful implementation of DGFM to treat the gas-water interface is crucial in the framework of the RKDG-DGFM fluid solver.
In Case 8, the rigid wall tube from Case 7 is replaced by a deformable wall tube made of AL5086 aluminum. Assuming that the structure density is constant and the entire process is adiabatic, the wall deformation can be modeled by the Lagrangian momentum equation:(12)ρDu⟶Dt=∇⋅σ=,where the stress tensor σ= is expanded using a stress-strain relationship for linear elastic materials. Typical properties of AL5086 are used in the calculation. The wall deformation is solved using a finite element method, while the ALE version of the fluid governing equations is solved using RKDG-DGFM with a moving wall condition. A partitioned FSI scheme is used to couple the solution of fluid and structure. Reference [14] provides more details for this FSI framework.
The RKDG-DGFM code is run at a final time of1.5×10−4 seconds. The pressure-time history at the inner center of the wall is shown in Figure 27. The pressure contours at various times are plotted in Figure 28. It is observed from the pressure contours that cavitation occurs at about 2.8×10−5 seconds in the water close to the wall, which is earlier than that in the rigid wall scenario. The reason is that this cavitation forms by a different mechanism that does not happen in the rigid wall case. The wall deforms because of the primary shock impact. The deformed wall generates an extra volume to be filled by the surrounding water and reduces the density and pressure of the water. This volume increase also makes the initial shock strength lower than that of the rigid wall. At around 3×10−5 seconds, local cavitation occurs midway between the tube and gas bubble due to the interaction between gas bubble and the reflected wave from the wall. Subsequently, the interaction generates an expansion wave traveling back toward the tube wall. These observations indicate that the first cavitation results from the lower wall pressure and the second cavitation results from the shock-bubble interaction as in the rigid wall tube case.
The pressure-time history from the RKDG-DGFM simulation matches closely with the experiment and Wardlaw’s work. The predicted peak pressure6.05×109dyne/cm2 is somewhat lower than the experimental value 6.25×109dyne/cm2. In the experiment, the first cavitation occurs at 3×10−5 seconds and then pressure reloads at 6×10−5 seconds. The RKDG-DGFM predicts that the first cavitation occurs at around the same moment, but the pressure reloads at an earlier time of 4.7×10−5 seconds. This is likely due to the constitutive model used in this work, which underpredicts the extent of wall deformation. The deformation of the wall in reality enters the plasticity region while the constitutive model used in this simulation is the hyperelastic model. As a result, it takes the surrounding water less time to rush into the cavitated region and collapse the cavitation. Unlike the DYSMAS calculation, a small cavitation region forms and collapses at the wall during between 4.7×10−5 and 5.8×10−5 seconds.
Figure28 are reproduced from “Application of the Runge–Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019” [33].
This 2D axisymmetric simulation uses an ALE method as the coupling approach between the fluid and the structural domain. This method can achieve high accuracy at the fluid-structure interface as the fluid mesh is always body-conforming. Special care, however, needs to be taken with the mesh-smoothing and remapping algorithm when the motion and deformation of the structural object are very large. In problems such as UNDEX, structural cracking and rupture might occur which cause a topological challenge that the ALE method cannot overcome. The embedded boundary method, therefore, will eventually be chosen as the coupling approach between the fluid and the structural domain. A 3D simulation of this case using the embedded boundary method will be performed and discussed in a subsequent paper.Figure 23
Initial conditions of Case 7 (left: 2D simulation; right: 3D simulation).Figure 24
Pressure-time history at the inner center of the rigid wall.Figure 25
Pressure contours inside the rigid wall tube (2D axisymmetric simulation).Figure 26
Pressure contours inside the rigid wall tube (3D simulation).Figure 27
Pressure-time history at the inner center of deformable wall.Figure 28
Pressure contours inside the deformable wall tube of half domain.


## 4. Conclusions

In this work, single- and multidimensional cases are examined to assess the applicability of the RKDG-DGFM to UNDEX modeling and simulation. Good comparisons of RKDG-DGFM results to analytic solutions and other reference results are achieved. The shock front and material interface are sharply captured at the correct locations and correct magnitudes are predicted. No spurious pressure oscillations, contour distortions, or excessive diffusion occur. Compared to the original GFM, the MGFM, and the LIM, the RKDG-DGFM has wider application and is not difficult to implement in multidimensions.The RKDG-DGFM is also applied to practical UNDEX applications including spherical bubble pulses, cavitating flow, and internal explosions inside rigid as well as deformable tubes filled with distilled water.Details of the fluid mathematical model, the discretization using modal discontinuous Galerkin method, and more simulation results using this hybrid framework of RKDG-DGFM are documented in a preceding paper [50]. Further development of the RKDG-DGFM fluid solver with the embedded boundary method (EBM) and coupling the fluid solver with a structural solver will be documented in a subsequent paper. The fluid solver based on this framework of algorithms is capable of more accurately predicting responses of both fluid and structure in real UNDEX events and eventually serving the purpose of early-stage ship design.



---

*Source: 1627382-2022-04-16.xml*