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

This work presents a Direct Ghost Fluid Method (DGFM) as part of a two-ﬂuid numerical framework suitable to model explosive gas and water ﬂows 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 diﬀusion may not be eﬀective. Numerical diﬀusion occurs due to nonphysical diﬀused density at material interfaces, which creates spurious pressure oscillations and signiﬁcantly degrades the quality of the numerical results. To eliminate or minimize numerical diﬀusion, sharp interface methods having no mixed elements may be used in multiﬂuid ﬂow 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 ﬂuid. The spurious pressure oscillations near the material interface are therefore minimized. One-, two-, and three-dimensional computational ﬂuid 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.


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. is is particularly complex in proximity to a structure and/or free surface where cavitation can be generated by reflected waves. is poses a difficult modeling problem. e 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 highpressure gaseous materials is completed so that there is no external heat source throughout the simulation [1,2]. e 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][6][7][8]. Compared to earlier methods [4,[9][10][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].
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]. e 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. is method uses marker points to track the interface. e location of marker points is decided via the solution of fluid equations [15][16][17]. As a cutcell 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. e error associated with distortion of the mesh may cause the degradation of numerical results and the sudden breakdown of the computation [19][20][21]. is 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. e Direct Ghost Fluid Method using direct extrapolation of density and tangential velocity across the material interface is developed for this purpose.
is 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. e 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. e 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. e 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.

Governing Equations.
e Euler equations in the conservative form are used as the fluid governing equations. e Arbitrary Lagrangian-Eulerian (ALE) form of governing equations are listed in the following equations.
where the terms , and E are density, velocity, convection velocity, identity tensor, and total energy per unit mass, respectively. ere 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. ere are also source terms if body forces exist in the model configuration. e difference between velocities of fluid flow, v → , and velocities of the moving grid, v → grid , is the convection velocity, m → . e convection velocity equals the velocity of fluid flow in a pure Eulerian description because the grid is fixed in space. e 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.
is system is closed by the equation of state. e 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]. e JWL EOS is where A, B, R 1 , R 2 , and ω are experimentally determined coefficients and ρ, e int , 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. e stiffened gas EOS is used to model the thermodynamic properties of water: where c B and B are experimentally determined coefficients. We have also used the Tait EOS: where N, B, and A are experimentally determined coefficients. e 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].
e Discontinuous Galerkin (DG) method [26][27][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. is 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]. e 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. is relaxes the element-by-element continuity requirement. e remaining terms are treated using the standard Galerkin method. e 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]. e 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. e high accuracy of the RKDG method helps to produce lessdiffusive 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. e 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.

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]. e GFM as implemented by Osher and Fedkiw used the Level Set Method (LSM) to perform this task. e 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]. e material interface is embedded in a spatially fixed grid. e transport of the LS equations is performed in an Eulerian manner [29,30]. e transport and reinitialization equations are zϕ zτ 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. e 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]. e 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]. e Hamilton-Jacobi WENO scheme and threestep 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 → . e 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 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. e 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]. e 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. e spurious pressure oscillation at the interface is therefore significantly reduced. Extrapolation of density in the 1D implementation of DGFM is illustrated in Figure 2. e density is extended from the real fluid into the ghost fluid by numerically solving the following advection equation: where n → � ∇ϕ/|∇ϕ| is obtained from the LS function [29]. A few iterations, say 25, are required to be performed. e 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 → � − n y i + n x j. e real velocity field is thus decomposed into normal and tangential velocity fields Tangential velocity of the real area is then extended to the ghost area as v t ′ by iterating equation (9) so that total velocity of the ghost fluid is a combination of real normal velocity and ghost tangential velocity 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, From this, we can obtain t → 2 based on n → and t → 1 so that these three vectors constitute a special orthonormal basis. e real velocity field is thus decomposed into one normal and two tangential velocity Both tangential velocities of the real area are then extended to the ghost area as v t1 ′ and v t2 ′ 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 Extrapolation of density and tangential velocity in 2D and 3D implementation of DGFM is illustrated in Figure 3.
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 → � u i → + vj + wk, the velocity of the real area is extended to the ghost area as (9). e real velocity field v → is used to obtain real normal velocity v n n → � , respectively. Total ghost fluid velocity is thus the sum of real normal velocity and ghost tangential velocity as v In the extrapolation from real gas to ghost gas, the density at element i − 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]. e 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.

Assessment
is 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. 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. e 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 × 10 9 dyne/cm 2 and ρ 0 � 1.0g/cm 3 is used to model the water. e ideal gas law with c � 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.
is 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. is example is explored using the 1D two-fluid RKDG-DGFM, which treats the moving interfaces conservatively. e initial conditions are shown in Figure 5. e [0, 100] cm domain is discretized with 200 uniform mesh. e 1D RKDG-DGFM simulation is run at 0.16 × 10 − 3 seconds. A CFL equal to 2 is used in the simulation. e 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.     [36]. is 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]. e 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. e 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.
is 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. e JWL EOS is used to model the explosive gas. e initial conditions are given in Figure 7. is case is initialized by a TNT-produced gas bubble with radius r 0 � 16cm surrounded by water at rest. e domain [0, 250]cm 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. e pressure and density profiles are provided. Figure 8 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. e shock peak pressure is 1.9 × 10 9 dyne/cm 2 . At 1.0 × 10 3 seconds, the material interface and shock front are located at 52cm and 205cm, respectively. e peak pressure is 0.75 × 10 9 dyne/cm 2 . e 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. ese were presented in [14]. e high-pressure gas produced strong shocks traveling into the water and simultaneously rarefaction fans moving toward the origin. e 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].

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. e center of the gas cylinder with radius r 0 � 50cm is located at the origin. e 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 Figure 11. Results are consistent with the xy-plots of the 2D results in Figure 10. No contour distortion occurs and symmetry is good. 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 Figures 14 and  15, respectively. e multidimensional results are consistent with their respective xy-plots in Figure 13. No contour distortion occurs and symmetry is good.

Practical UNDEX Applications.
e performance of RKDG-DGFM solvers assessed in Cases 1-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.

UNDEX Explosion Bubble.
Cases 5 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. e 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].    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 TNTproduced 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 with 28kg of TNT explosive surrounded by water at a depth of 178m. e TNT charge is replaced by a gas bubble with equivalent initial volume and initial internal energy of the original High Explosive (HE). e initial bubble radius r 0 is determined by r 0 � (3V 0 /4π) 1/3 , where V 0 is the initial volume of the explosive charge. e 1D spherical domain [0, 10000]cm is discretized with 2000 nonuniform elements with a refinement towards the charge center. e center of the gas bubble with r 0 � 16cm is located at the origin. e initial conditions are specified in Figure 17. e 1D spherical RKDG-DGFM simulation is run at 0.15 seconds. e same boundary conditions as in Case 2 are imposed on the two ends, respectively. Figure 18 shows the calculated bubble radius-time history and the pressure at 121cm from the origin. e bubble radius is estimated using the sign change of the LS function ϕ. e maximum bubble radius is approximately 220cm at 0.07 seconds. e minimum bubble radius is approximately 30cm at 0.13 second.
ese results compare well with   Wardlaw's results.  [42]. ey conducted measurements on the oscillations and maximum radii of gas bubbles from underwater explosions [42]. ese experiments provide a useful bubble radius-time history for a 300g charge detonated at a depth of 91.46m. e experiments continue up to the second maximum radius. e 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]. e domain [0, 9146]cm is discretized with 2000 nonuniform elements. Compared to Case 5, a smaller initial radius r 0 � 3.529cm corresponding to a 300g TNT charge is used. e 1D spherical RKDG-DGFM simulation is run at 0.045 seconds. e 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.

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. e fluid next to the structure is in tension and a local cavitation region develops at the fluidstructure interface. After reaching its maximum inward velocity, the structure begins to slow down and reverse direction or rebound. e local cavitation region closes and reloads the structure [1,41,45]. e bubble/shock cavitation mechanism associated with wall interaction is depicted in Figure 21.
is section explores the performance of the RKDG-DGFM framework in its simulation of near-field UNDEX where the gas-water-cavitation-structure interact.
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].
In the cutoff cavitation model, the low-pressure limit p limit 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 of 3g PETN internally. is experiment was performed by Sandusky et al. [48,49]. e 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. e 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]. is 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]. e experimental arrangement and its corresponding initial conditions in the mathematical model are shown in Figure 22. 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.9]cm is discretized using 40 × 160 uniform elements. In the 3D simulation, the rectangular domain of [− 6, 6] × [− 8.9, 8.9] × [− 6, 6]cm is discretized using 120 × 178 × 120 uniform elements. e embedded boundary method is used to treat the wall boundary condition. Details of the embedded boundary method will be addressed in a subsequent paper. e JWL EOS with PETN parameters is used to model the thermodynamics of the explosive gas. e Tait EOS is used to model the thermodynamics of the water [9,11,46,48]. Figure 22 is reproduced from "Application of the Runge-Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019" [33]. case 9. first assumes that the wall is rigid instead of deformable. e explosion occurs inside this rigid wall. In the 2D axisymmetric simulation, a symmetry condition is applied along the z-axis. Both the top and bottom are treated  using a nonreflecting boundary condition (NRBC) to minimize the interference of reflection waves from these boundaries. e initial conditions of 2D axisymmetric and 3D simulations are set up and shown in Figure 23. e RKDG-DGFM simulation is run at 1.5 × 10 − 4 seconds. e 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. Figure 25 is reproduced from "Application of the Runge-Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019" [33].
e 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, at 1.5 × 10 − 5 and 1.2 × 10 − 4 seconds, are predicted to be 6.6 × 10 8 Pa and 3.2 × 10 8 Pa. ese values are close to the reference results taken from [46]. e 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. e domain in the 3D fluid has flat sides between nodes at the fluid-wall boundary, while in reality it is curved. e embedded structural wall penetrates the fluid mesh and it penetrates at different angles along the perimeter. is 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. is problem is alleviated as the fluid mesh is refined. e 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 around 4 × 10 − 5 seconds. e 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. is 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. is 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. ese 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: where the stress tensor σ � is expanded using a stress-strain relationship for linear elastic materials. Typical properties of AL5086 are used in the calculation. e 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.
e RKDG-DGFM code is run at a final time of 1.5 × 10 − 4 seconds. e pressure-time history at the inner center of the wall is shown in Figure 27. e 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. e reason is that this cavitation forms by a different mechanism that does not happen in the rigid wall case. e wall deforms because of the primary shock impact. e deformed wall generates an extra volume to be filled by the surrounding water and reduces the density and pressure of the water. is 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. ese observations indicate that the first cavitation results from the lower wall pressure and the second cavitation results from the shockbubble interaction as in the rigid wall tube case. e pressure-time history from the RKDG-DGFM simulation matches closely with the experiment and Wardlaw's work. e predicted peak pressure 6.05 × 10 9 dyne/cm 2 is somewhat lower than the experimental value 6.25 × 10 9 dyne/cm 2 . In the experiment, the first cavitation occurs at 3 × 10 − 5 seconds and then pressure reloads at 6 × 10 − 5 seconds. e 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.
is is likely due to the constitutive model used in this work, which underpredicts the extent of wall deformation. e 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.   Figure 28 are reproduced from "Application of the Runge-Kutta Discontinuous Galerkin-Direct Ghost Fluid Method to internal explosion inside a water-filled tube, 2019" [33].
is 2D axisymmetric simulation uses an ALE method as the coupling approach between the fluid and the structural domain. is 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. e 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.

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. e 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. e 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. e 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.

Data Availability
e data used to support the findings of this study may be released upon request to the Kevin T. Croft Department of Aerospace and Ocean Engineering, Virginia Polytechnic Institute and State University, which can be contacted at albrown5@vt.edu.

Additional Points
e Direct Ghost Fluid Method is presented and utilized to model UNDEX interface conditions. Simulations are compared with analytic and experimental results in a series of benchmark problems. e DGFM decreases the density diffusion across the interface between the explosive gaseous products and the surrounding water. Spurious pressure oscillations at the material interface are therefore minimized. e DGFM is also successfully applied to the simulation of an explosion inside a tube filled with distilled water.