Optimization of the G & H bubble model

Abstract. A spherical model for the bubble created by an underwater explosion is enhanced to account, in approximate fashion, for the effect of bubble distortion on translation and dilation. The enhancement consists of introducing artificial drag in the form C|ν|P , whereν is translation velocity, and C, P produce an optimum fit to empirical formulas for the second dilation maximum and first two translation jumps. The recommended values are C = 0.4 andP = 1 for charge weights between 100 lb–1000 lb and and depths exceeding 200 ft.


Introduction
Recently, Geers and Hunter developed a bubble bubble model for a very deep to moderately deep underwater explosion that integrates the shockwave and oscillation phases of the motion [2].The model was then specialized to the case of no bubble distortion, thereby limiting the motion to dilation plus translation.The specialized model produced bubble-motion histories exhibiting predictive capability superior to that of earlier bubble models.
The principal limitation of the specialized model is that its neglect of bubble distortion results in the over-prediction of bubble translation, which in turn affects bubble dilation.Hicks reviously attempted to remedy such deficiency by introducing the hydrodynamic drag force F D (t) = πa 2 (t) • 1  2 ρ D(t), in which a(t) is the bubble radius, ρ is the density of the liquid surrounding the bubble, and D(t) = C D u2 (t), where u(t) is the translation velocity of the bubble [3].Hicks recommended C D = 2.25, which is an order of magnitude larger than the typical value for a rigid sphere at high Reynolds number.Hunter and Geers found that neither this nor any other C D -value simultaneously yields accurate dilation and accurate translation at moderate depths [4].
In principle, the above limitation of the specialized model can be relieved by using the boundary element method to implement the more general model.However, this increases the number of degrees of freedom by one or two orders of magnitude.Furthermore, the resulting non-spherical collapse process is difficult to treat [1].Hence, it is worthwhile to seek an approximate technique that accurately accounts for the effects of bubble distortion on translation and dilation within the confines of the specialized model.This is the motivation for the work reported here.

Equations of motion for the specialized bubble model
With a(t) and u(t) as the spherical bubble's radius and vertical displacement, with φ D (t) and φ T (t) as the dilation and translation velocity potentials for the external liquid at the bubble surface, and with ϕ T (t) as the translation velocity potential for the internal gas at the bubble surface, the oscillation-phase equations of motion for the specialized bubble model, including the drag function D(t), are [ In these equations, c and c g (t) are the speeds of sound in the liquid and the gas, ρ g (t) is the density of the gas, , is the uniform component of pressure in the gas defined by its equation of state, P atm is atmospheric pressure, g is the acceleration due to gravity, and d is the initial depth of the explosive charge.In the explicit numerical integration of Eq. ( 1), the right sides of (1c) and (1d) are used for φD and φT on the right sides of (1a) and (1b).
If c (t) and ρ g (t) are both set to zero, and D(t) is taken as C D u2 (t), Eq. ( 1) reduce to the equations of Hicks [3]: These equations clearly exhibit the coupling between dilation and translation.Of particular interest are the terms gu (in the first equation) and 2ga (in the second).The first accounts for the change in ambient pressure as the bubble rises, and the second accounts for the buoyancy force that causes the bubble to rise.Although insights provided by Eq. ( 2) are valuable, the equations are not accurate predictors of bubble motion [2].
Shown in Fig. 1 are response histories for bubble dilation and translation produced by Eq. ( 1) for a 1000-lb charge of TNT detonated at a depth of 250 ft; one pair of histories pertains to C D = 0 and the other pair to C D = 2.25.The horizontal and vertical lines in Fig. 1(a) mark empirical values for maximum bubble radius and time of minimum bubble radius, respectively.The asterisks and circles in Fig. 1(b) mark empirical values for bubble depth at the times of bubble maximum and bubble minimum, respectively.These values are based on extensive experimental data [3,5,6].In [5] Snay characterizes the data beyond the third bubble minimum as "sparse and uncertain".
The histories in Fig. 1 demonstrate that the use of C D = 0 over-predicts both dilation and translation, while the use of C D = 2.25 under-predicts these responses.This implies that values between these extremes might work better.In fact, a drag formula that admits velocity to a power other than two might be better still.These considerations are only pertinent to motion during and after the first bubble minimum, as the bubble undergoes negligible translation before that time.

Modification of the drag formula
Seeking to improve model performance without increasing computational effort, we replaced the drag formula with the ad hoc relation Then we searched for optimum values of C and P over practical ranges of charge weight W and detonation depth d.We based our search on accurate computation of the wave field generated by an explosion bubble rather than on accurate prediction of bubble motion per se.Specifically, emphasis was placed on accurately predicting the first and second pressure pulses generated by bubble collapse/rebound.Bubble dilation is the primary contributor to these pulses, translation is a secondary contributor, and distortion is a minor contributor.The pulse characteristics are significantly affected by wave effects in both the external liquid and internal gas, as well as by bubble depth [4].
In accordance with the above criterion, the search employed an error measure computed with the following procedure: 1.For specified values of W , d and C, P , perform a bubble simulation with the G&H equations (Eq.( 1)) and record the values a(t 2 ), u(t 1 ), u(t 2 ), and u(t 3 ), where t n is the time of the n th bubble maximum.2. For the same values of W and d, calculate a(t 2 ), u(t 2 ) − u(t 1 ), and u(t 3 ) − u(t 2 ) from the empirical formulas of Snay [5].

With the results produced in the preceding steps, compute the errors
Here, ā and h are the equilibrium radius and initial head, respectively, given by where a c is the charge radius, K c is the energetic constant for the charge material [2], and γ is the ratio of specific heats for the bubble gas.4. With Eq. ( 4), compute the error magnitude  Fits with polynomials of higher order were also carried out, with little improvement.B. C opt and P opt were considered functions of the parameter η = h/ā, where, from Eq. ( 5), h is a function of d and ā is a function of W and d.With the 20 η-values as points on a horizontal η-axis, Cand P -values were sought for each η-value that minimized ε.Then the curves C opt (η) and P opt (η) were determined by performing least-squares fits of the quadratic representation c 0 + c 1 η + c 2 η 2 .Fits with higher-order polynomials were also carried out, to little effect.C. C opt and P opt were considered constants.The ε-values were grouped according to their C, P pairs, and the average and maximum values of ε for each C, P pair were recorded in Table 1.

Computational results
The most successful of the optimization procedures was the last.Table 1 produced by that procedure is shown below.We see that the errors in the vicinity of C = 0.4, P = 1.0 are the smallest.As the maximum error for C = 0.4, P = 0.8 is slightly smaller than that for C = 0.4, P = 1.0, we performed additional simulations for C = 0.4 and P < 0.8. Figure 2 shows representative translation histories generated by these simulations.We see that a sufficiently small value of P causes large negative translation.For all of the simulations performed, a minimum P -value of unity was safe in this regard.Hence, the recommended drag formula for all cases in the ranges  (7) Shown in Fig. 3 are response histories for bubble radius and translation produced by the G&H equations with Eq. (7) for a 1000-lb charge of TNT detonated at a depth of 250 ft.We see that overall agreement with experimental data is better than that seen in Fig. 1.

Fig. 1 .
Fig. 1.Dilation and translation histories for a 1000-lb charge of TNT detonated at a depth of 250 ft, calculated with two different hydrodynamic drag coefficients.C D = 0.0 (dashed), C D = 2.25 (solid).

Fig. 2 .
Fig. 2. Translation histories for a 1000-lb charge of TNT detonated at a depth of 250 ft, calculated with C = 1.0 and various values of P .Thirty-six hundred G&H simulations and corresponding Snay calculations were performed for W = 100, 400, 700, 1000 lb; d = 200, 250, 300, 400, 500 ft; C = 0.2 (0.2) 3.0; P = 0.8 (0.2) 3.0.Three different optimization procedures were pursued, as follows: A. C opt and P opt were considered functions of W and d.With the 20 W, d pairs as points on a horizontal W -d plane, the Cand P -values were sought for each W, d pair that minimized ε.Then the surfaces C opt (W, d) and P opt (W, d) were determined by performing least-squares fits of the quadratic representation c 0 + c 10 W + c 01 d + c 20 W 2 + c 11 W d + c 02 d 2 .Fits with polynomials of higher order were also carried out, with little improvement.B. C opt and P opt were considered functions of the parameter η = h/ā, where, from Eq. (5), h is a function of d and ā is a function of W and d.With the 20 η-values as points on a horizontal η-axis, Cand P -values were sought for each η-value that minimized ε.Then the curves C opt (η) and P opt (η) were determined by performing least-squares fits of the quadratic representation c 0 + c 1 η + c 2 η 2 .Fits with higher-order polynomials were also carried out, to little effect.C. C opt and P opt were considered constants.The ε-values were grouped according to their C, P pairs, and the average and maximum values of ε for each C, P pair were recorded in Table1.