Spherical solutions of an underwater explosion bubble

The evolution of the 1D explosion bubble flow field out to the first bubble minimum is examined in detail using four different models. The most detailed is based on the Euler equations and accounts for the internal bubble fluid motion, while the simplest links a potential water solution to a stationary, Isentropic bubble model. Comparison of the different models with experimental data provides insight into the influence of compressibility and internal bubble dynamics on the behavior of the explosion bubble.


Introduction
Simulation of underwater explosions is an essential component of platform vulnerability and weapon lethality assessments. Success at this task requires accurate prediction of target loading, which primarily occurs as a consequence of the initial shock and subsequent bubble collapse. These damage mechanisms occur over a wide range of time scales, from submillisecond for shock response to hundreds of milliseconds for response to bubble pulses. Unfortunately, three-dimensional Eular simulations of such events, particularly bubble collapse, tax the ability of supercomputers. This has resulted in the use of less complete models for water and explosion products (e.g., Van Tuyl and Collins [11], Szymczak et al. [10] and Chahine [1]). The objective of this paper is to examine the consequence of these simplifications by studying the one-dimensional (1D) explosion bubble.
The 1D explosion exhibits two important damage mechanisms: the initial shock, and subsequent pressure pulses from bubble collapse and rebound. The limited size of this problem makes it numerically tractable, allowing mesh converged solutions to be achieved to provide insight into mesh resolution requirements in 2D and 3D. Shock and interface fitting is also easily implemented in this single-dimensional setting. This problem does not exhibit bubble jetting, a form of collapse that induces a directed water jet. However, it is reasonable to assume that accurate treatment of the 1D problem is a necessary prerequisite for the prediction of jetting.
This paper examines four different numerical descriptions of the underwater explosion. The simplest, termed the Incompressible Model, assumes water to be incompressible and the bubble to be a stationary, homogenous gas, whose pressure changes isentropically with bubble volume. The second approach, termed the Isentropic Euler Model, includes a bubble description similar to that of the Incompressible Model, but introduces the compressibility of the water. The third approach, termed the Uniform Euler Model, adds an Euler description of the bubble, and starts the calculation with uniform bubble properties. The fourth approach, termed Non-Uniform Euler Model, is based on nonuniform initial bubble conditions, derived from a Taylor wave solution to the gaseous explosive products.
The flow field evolution predicted by these four models is examined and compared. Also, solution sensitivity to changes in water and explosion product equations of state is considered. The paper concludes by comparing calculations with experimental data. Comparisons are made only through the first bubble minimum; multidimensional dissipative mechanisms influencing bubble evolution beyond the first minimum are outside to scope of this paper, but are addressed elsewhere (e.g., Geers and Hunter [5]).

Solution methodology
Two types of material descriptions are considered in this paper: compressible (Euler) and incompressible. These methods have been assembled into the four different models mentioned in the Introduction.This section details the equations and numerical methods of the compressible and incompressible descriptions, while the following section defines the specifics of the individual models.

Compressible (Euler) calculations
The water and the gas bubble are described using the Euler equations in spherical coordinates: where To ensure a high quality solution, the numerical method avoids mixed cells that contain gas and water. This is achieved using an Arbitrary Lagrangian/Eulerian (ALE) approach which moves the gas/water interface node at the computed interface velocity. A consequence of this methodology is that the number of grid points is constant throughout the calculation, in both the bubble and the surrounding water. Equations (1) are solved using a second order Godunov method based on the technique described by Collela [2]. However, this method has been modified to contain both a Lagrangian step and a re-map step, as is described in Appendix A. Fitting the shock to the outer edge of the computational mesh during the shock phase allows the shock jump to be explicitly calculated rather than captured by the numerical scheme.

Incompressible calculations
In the Incompressible Model, the water is described using the potential solution for an expanding sphere (e.g., Karamcheti [7]). The pressure on the surface of the bubble is given by: where r(t) is the bubble radius and p ∞ is the ambient pressure. Equation (4) is coupled to the equation of state for the bubble, which in this study is taken to be of the JWL form. By virtue of the isentropic assumption, the equation of state provides pressure as a function of density, which in turn is a function of initial bubble mass (M ) and the bubble volume. This reduces the bubble pressure to: Combining Eqs (3) and (4) yields a differential equation which can be integrated numerically, starting with a known radius r 0 and dr/dt = 0.

Incompressible Model
The Incompressible Model describes the water using the incompressible formulation of Eq. (4). The bubble is assumed to be a stationary, isentropic gas governed by the JWL equation of state, Eq. (3), and associated parameter set. Starting an incompressible calculation from the detonation conditions (i.e., ρ = 1.63 g/cm 3 , e 0 = 4.29 × 10 10 d cm/g for TNT) results in a period that is too long; it is therefore common practice to apply empirical conditions which effectively remove the energy associated with the shock from the calculation. Typically, a set of empirical formulas derived from experiment and tailored to a gamma law gas bubble model are used to set the initial bubble conditions (e.g., Szymczak et al. [10]). In the present study, the initial conditions are chosen by adjusting the initial bubble radius and setting the initial pressure by isentropic expansion from the explosion conditions. Numerical experiments indicated that an initial bubble radius 1.8 times the actual charge radius produces a solution that closely matches the Non-Uniform bubble period. The motivation for selecting this initial condition is to connect the initial empirical bubble pressure to that of the initial explosive state used in the other models while also assuring a match between the predicted incompressible and compressible bubble period. This initial condition prescription is not being offered as a replacement for tradition empirical starting conditions, which have been verified over a broad range of cases.

Isentropic Euler Model
Here the bubble is similar to the isentropic, stationary bubble used in the Incompressible model while the water is treated as compressible via the Euler equations. Ideally, the bubble gas would be described using Eq. (5); however, to simplify model construction, a single Euler cell defines the bubble. The inner edge of the cell is located at the center of the bubble while the outer edge is located at the gas/water interface and moves at the interface velocity. This model produces a homogeneous bubble; however, the velocity of the bubble gas is not zero. Hence, the internal dynamics of the bubble are not completely eliminated, and the Isentropic Euler and Incompressible models differ by more than just the compressibility of the water. However, an examination of the solution indicates that the gas bubble velocity is small and the bubble kinetic energy is negligible. The initial bubble conditions are uniform and the starting state is that for TNT, the same conditions employed in the Incompressible model as the base state for the isentropic expansion.

Uniform Euler Model
The bubble and water are treated as compressible fluids which are described by the Euler equations. This model differs from the Isentropic Euler Model only with regard to the bubble, that is now described by a large number of cells. The same uniform, TNT conditions are used as the initial bubble state in the Uniform Euler and Isentropic Euler models.

Non-Uniform Euler Model
The Euler solution described in the Uniform Euler Model is applied here and the Uniform and Non-Uniform Euler models differ only with regard to the initial bubble conditions. The Non-Uniform and Uniform initial bubble conditions contain the same total energy, however, for Non-Uniform Euler case the distribution of properties within the bubble is determined using a spherical analog to the Taylor plane wave solution based on the work of Guirguis et al. [6]. The resulting bubble property profile features the highest pressure at the bubble/water interface with a constant, minimum pressure region at the bubble center (see Fig. 1, t = 0 curve).

Model comparisons
The four bubble models are compared for the following conditions: detonation of a 28 kg TNT sphere (radius of 16 cm) at a depth of 91.7 m (ambient pressure of 1.0 × 10 7 d/cm 2 ). The Euler solution mesh size was selected through a convergence study and features 512 points in the bubble and 565 point in the water. The same mesh size is used for both the shock and collapse phase; however, the shock is only tracked for the shock   phase. The required run time on a Pentium II 266 MHz PC is one hour for a shock problem, and two hours for a bubble collapse calculation.

Figures 1-3 trace the evolution of the Non-Uniform
Euler solution through early times, where the shock is the dominant feature in the flow field. This "shock phase" can be divided into three sub-phases. The first is illustrated in Fig. 1 and depicts the initial interaction between the high pressure bubble and the low pressure water. This results in a shock that propagates into the water and an expansion that travels back into the bubble. The sub-phase terminates when the inward moving expansion reaches the origin. Figure 2 illustrates the reflection of the expansion from the origin, which marks the start of the second sub-phase. As expected, an expansion reflects as an expansion, dropping the minimum pressure in the bubble several orders of magnitude. However, the pressure gradient inside the bubble and the accompanying inward flow impedes the outward progress of the expansion. A sharp pressure gradient develops in the bubble that evolves into a shock moving toward the center of the bubble. The second sub-phase ends when that shock reaches the origin.
The third sub-phase begins when the shock formed in the second sub-phase reflects from the origin. The reflected shock propagates to the gas/water interface and interacts with it as is shown in Fig. 3. The interaction generates a reflected shock which moves back into the bubble toward the origin, and a transmitted shock that travels outward from the bubble into the water. The third sub-phase ends when the reflected shock reaches the origin.
The process of shock reflection from the origin and interaction with the interface repeats as the bubble expands. Flow field profiles at later times bear the legacy of these interactions; several pressure pulses in the water can be seen moving outward from the bubble. As the bubble expands, pressure levels inside the bubble decrease, the shock inside the bubble weakens, and the intensities of waves transmitted from the bubble to the water become negligible.
Comparison of the four solutions during the shock phase is shown in Figs 4 and 5. The Uniform Euler Model solution is qualitatively similar to the Non-Uniform Euler Model solution. The initial sub-phase 1 expansion reflects from the origin and subsequently forms an inward moving shock. This shock repeatedly reflects from the origin and interacts with the interface, generating transmitted pressure pulses in the water, several of which are shown in Fig. 5. The principle differences between these solutions are the times at which the events occur. The Uniform Euler solution unfolds at a quicker pace, as can be seen in Fig. 5 from the more advanced position of the transmitted pulse.

Bubble collapse phase
The evolution of the Non-Uniform Euler Model flow field near the bubble minimum volume is depicted in Fig. 6. The pressure within the bubble is highly Non-Uniform, and a careful examination of flow field data indicates that the moment of maximum bubble pressure does not coincide with the time of minimum volume. Although the maximum bubble pressure occurs near the point of minimum volume, the peak pressure will coincide with a time at which a shock wave inside the bubble converges at the bubble center. Figure 6 also exhibits a water pressure profile with numerous local peaks, which are a consequence of the interaction of the shock wave inside the bubble with the gas/water interface.
The solutions for all four models through bubble collapse are compared in Figs 7, 8 and 9. The bubble radius history is shown in Fig. 7, while Fig. 8 illustrates the bubble-water interface pressure. The bubble period is marked by the minimum volume in Fig. 7 and the pressure peak in Fig. 8. The three Euler models feature significantly different bubble periods, varying from 0.127 to 0.136 seconds. The incompressible solution period has been tailored to match that of the Non-Uniform Euler Model. Note that the maximum bubble radius differs between these two cases. Use of the explosive initial conditions for the incompressible solu-     tion rather than the empirical ones yields a bubble period of 0.176 seconds.
The Non-Uniform and Uniform Euler Model interface pressures oscillate due to the interactions of the shock waves inside the bubble with the gas/water interface. This type of oscillatory interface behavior has been previously observed by Mader [8]. The Incompressible and Isentropic Euler solutions, on the other hand, do not have an internal bubble structure and consequently have smoothly varying interface pressures.
The pressure field at bubble collapse is illustrated in Fig. 9 for all four solutions. The three Euler solutions display similar water pressure levels. However, the Non-Uniform Euler Model and Uniform Euler Model solutions exhibit local pressure peaks generated by the gas/water interaction, while the Isentropic Euler Model and Incompressible Model solutions do not. In addition, the Incompressible Model solution exhibits a pressure far from the bubble which exceeds that of the compressible models.
One method of examining the differences among the four models is to consider the manner in which the energy transmitted to the water is partitioned. In the Euler solutions this energy can be converted to kinetic energy or internal energy via compression. Incompress-ible media do not have internal energy, but there is a potential energy associated with the work required to move water from the space occupied by the expanding bubble. For the Euler solutions, the energy partitioning must be determined by numerical integration, but for the Incompressible Model an analytic expression can be derived: Here the surface pressure, p s is evaluated from Eq. (4) and u = dr/dt. The first term is the potential energy while the second term is the kinetic. The results are shown in Fig. 10, which graphs the fraction of the transmitted energy that is converted into kinetic energy. Comparison of the three Euler models indicates that the bubble period is directly related to the water kinetic energy level. Examination of these solutions also reveals a relationship between initial shock strength and kinetic energy, with the lowest kinetic energy associated with the strongest initial shocks. Presumably, Fig. 10. Fraction of water energy converted to kinetic energy. a higher shock pressure increases the water compression, leaving a smaller amount of energy for the kinetic component that propels the water away from the bubble. As can be seen in Fig. 10, the incompressible, explosive initial condition model has a large kinetic energy, which is consistent with its lengthened period. Use of the empirical initial conditions reduces this energy greatly.
As the bubble collapses, water rushes toward it to fill the space left by its diminishing volume. At the point of minimum volume, this inflow comes to a halt. In the compressible models, some of the water kinetic energy is converted to internal energy by compression and heating; however, in the incompressible case, the only mechanism available to remove energy is to transfer it back to the bubble, which regains its initial energy and pressure. Deceleration of the water in the Incompressible model is accomplished only through a pressure gradient. This accounts for the different functional form of the incompressible pressure visible in Fig. 9, that exhibits higher pressures farther from the bubble.

Euler Model uncertainties
Euler bubble models are dependent on the equations of state for water and the explosion products. These re-lations contain a number of parameters that must be assigned. Additionally, in the case of water, there are several forms of the equation of state that might be used. Before comparing with experiment, it is useful to examine the impact of these factors on the solution.
To gage the sensitivity of the solution to the explosive JWL model, a different set of coefficients, de- veloped for underwater studies by Fiessler [4] and labeled TNT2 is applied. Here A = 5.484 × 10 12 d/cm 2 , B = 9.375 × 10 10 d/cm 2 , R 1 = 4.94, R 2 = 1.121, ω = 1.2801, e 0 = 5.183 × 10 10 d cm/g. The influence of these changes in equation of state on the Uniform Euler bubble period is shown in Table 1. Analogous changes occur in the other Euler solutions. It is evident that the largest uncertainties are associated with the JWL parameters and initial energy. Table 2 shows a comparison among the four models and the experimental data of Swift and Decius [9] for a case featuring 300 g of TNT at 91.4 m. For all cases the Tait equation of state was used for water, in conjunction with the TNT2 set of JWL coefficients for the detonation products; this set was selected because it produced the best agreement between experiment and the Non-Uniform Euler Model. The results shown in this table are consistent with Figs 7 and 8; the Non-Uniform Euler Model exhibits the shortest period, while the Isentropic Euler Model features the longest. Two incompressible results are shown, labeled "incompressible", in which the initial explosive conditions are applied to the bubble, and "incompressibleadjusted", in which the empirical initial conditions are used. Comparison of the Incompressible and Isentropic Euler cases demonstrates the influence of water compressibility on the solution. Table 3 compares the Non-Uniform Euler Model results to three additional test cases taken from Swift and Decius [9] and Yennie and Arons [12]. Close agreement is obtained in all relevant cases (Case 1 is included only for reference, since the charge is located close enough to the free surface to be influenced by it).

Comparison with experiment
In Fig. 11, computed collapse pressures are compared with the experimental data of Yennie and Arons [12]. The test was conducted using 227 g of TNT at a depth of 152.4 m and the pressure was measured at a location 69.5 cm from the center of the charge. To aid in a comparison of peak pressures, the computed curves have been individually translated by the amount shown in the figure legend. As can be seen from this figure, the incompressible peak pressure is the highest, followed by the Isentropic Euler Model. The Uniform and Non-Uniform Euler solutions are oscillatory in character, reflecting the interaction between the bubble and the water. The close agreement between the incompressible and Euler solutions is somewhat fortuitous. The peak incompressible pressure is sensitive to changes in the bubble period; a 5% increase in the bubble period will double the peak pressure while a similar decrease will halve it. This sensitivity is a consequence of requiring the empirical initial state to be on the initial TNT isentrope and the mass of the bubble to be that of the initial explosive. To change the bubble period requires a change in initial bubble radius, that alters the initial pressure via Eq. (5) and consequently the collapse pressure. The standard empirical methods for generating initial conditions do not levy these requirements on the initial state and thus are free to best adjust the bubble initial radius and pressure independently (e.g., Szymczak et al. [10]).

Conclusions
The spherical explosion bubble contains an internal shock that repeatedly reflects off of the bubble center and interacts with the gas/water interface. This interface interaction generates a transmitted shock that travels outward into the water, and a reflected shock that moves inward towards the bubble center. A consequence of this internal bubble structure is a fluctuating gas/water interface pressure, as well as pressure variation at fixed points throughout the flow field.
All of these features are contained in the Non-Uniform Euler Bubble Model. As the description of the bubble is simplified, the following changes in solution result: 1) Substituting the uniform initial bubble conditions for the non-uniform ones lengthens the bubble period on the order of 5%. 2) Eliminating the internal bubble structure (Isentropic Euler Model) removes the bubble shockinterface interaction phenomenon which leads to smooth bubble-water interface pressures and lengthens the bubble period by about 5%. 3) Eliminating the compressibility of water increases the bubble period by about 25%. For the Euler solutions, the strength of the initial explosion shock plays a key role in partitioning the energy transmitted to the water between kinetic and internal. The stronger the shock, the larger the increase in water internal energy and the smaller the kinetic energy. This diminished kinetic energy decreases the outward velocity of the water, reduces the maximum bubble radius, and shortens the bubble period.
Differences among three water equations of state were tested and all produced bubble periods within 1 percent. Changes in the coefficients and initial energy for the JWL equation of state yielded variations of 5 percent.
The explosion bubble solutions have been compared to data from four experiments, with depths varying from 1.5 m to 179 m. Calculated bubble period and maximum radius were in reasonable agreement with experiment. The computed pressure at a point in the flow field agreed reasonably well with experiment at collaps time. (ii) Convect cell and recompute cell properties, denoted by superscript L, via the following: where p = (p * i−1/2 + p * i+1/2 )/2 and V = cell volume.
2. Remap step.  Here D i+1/2 is the distance cell edge i + 1/2 moves during step ∆t, and quantities are evaluated at the location x = x i+1/2 ± ∆d i+1/2 /2, via the slope of step 2(b). The sign is selected to insure that the upwind properties are used.
(d) Compute re-mapped cell properties using: Here V n+1 is the new cell volume. In the case in which the outer shock is tracked, the outermost cell is set to the ambient state, while its neighbor contains the state behind the shock. An exact Riemann problem is solved across the inner edge of the outer cell; the cell edge velocity is assigned that of the shock. The resulting solution exhibits the entire shock jump across this cell edge.