The reflection of a blast wave by a very intense explosion

We demonstrate that the geometric similarity of Taylor’s blast wave persists beyond reflection from an ideal surface. Upon impacting the surface, the spherical symmetry of the blast wave is lost but its cylindrical symmetry endures. As the flow acquires dependence on a second spatial dimension, an analytic solution of the Euler equations becomes elusive. However, the preservation of axisymmetry, geometric similarity and planar symmetry in the presence of a mirror-like surface causes all flow solutions to collapse when scaled by the height of burst (HOB) and the shock arrival time at the surface. The scaled blast volume for any yield, HOB and ambient air density follows a single universal trajectory for all scaled time, both before and after reflection.


Introduction
The Taylor [1,2], von Neumann & Richtmyer [3] and Sedov [4] solutions for a self-similar blast wave in the strong-shock limit have been used for eight decades to estimate the yields of nuclear tests and to explain the behaviour of supernovae [5], stellar wind bubbles [6] and other high-energy phenomena that produce strong shock waves [7,8]. During the atmospheric-testing era, weapon scientists employed Taylor's equation to estimate the yields of 210 atmospheric tests by analysing high-speed camera films [9,10]. In most of these tests, the devices were detonated sufficiently high above the ground to ignore the shock interaction with the surface. In at least 14 of these events, however, the blast wave reached the ground too soon for Taylor's equation to correctly apply to the data.       distances in these cases appear to exceed Taylor's prediction due to reflected energy reinforcing the upward-moving shock. Consequently, their yields appear to be up to 40% higher than the yields determined by radiochemical techniques, light-curve analyses and other methods. The Mohawk event, for example, exhibited a large discrepancy in the yield inferred by Edgerton, Germeshausen and Grier, Inc. (EG&G), when comparing the blast volume with the light output; i.e. the blast volume suggested a yield of approximately 360 kilotons whereas the light output indicated approximately 248 kilotons [11]. EG&G reasoned that spots on the fireball surface and the appearance of a bulbous jet, seen in figure 1, corrupted the light measurements; hence they rejected the light-curve data and reported the higher yield. EG&G, however, did not account for the reflection of the blast wave from the ground surface, which significantly increased the blast volume. The purpose of this paper is to lay out the additional reflection equations, which must be applied to film data in order to properly account for reflected blast energy. Throughout this paper, we use the term 'blast' to refer to the shock front, rather than the fireball surface, although the latter is typically easier to identify on the films.

Dimensional analysis
Taylor's blast wave [1,2] for air with zero ambient energy and pressure is governed by only two parameters possessing physical dimensions; the energy of the point source  and the density of the air into which the energy is deposited Combining Y a and ρ a to eliminate M results in a single fixed relationship between L and T . In Taylor's infinite atmosphere, where no intrinsic length or time scales exist, L and T remain locked together (L 5 ∝ T 2 ), such that neither can abide independent of the other. The blast radius must therefore scale with its own time of propagation, where the similarity constant S γ is solely a function of γ . The fixed relation between L and T means that the blast wave is completely specified either by its radius or time of propagation. The presence of a boundary introduces a third dimensional parameter into the system, H[=]L, along with a time scale, t H [=]T associated with the new parameter. For a detonation at distance H above a perfectly reflecting (mirror-like) surface, the blast wave will reach the surface at time For t > t H , the fixed relation between L and T is broken, the spherical symmetry is lost and R s ceases to parametrize the flow. However, the flow remains axisymmetric and the blast can still be fully characterized by its net volume, regardless of the complexity of the Mach stem and triple point emerging from the reflection. Since the ideal surface acts as a plane of symmetry, the flow behaves as if there were two identical blasts, separated a distance 2H apart; i.e. a blast and its mirror image [12]. Prior to reflection, the blast volume can be written in terms of (2.3) For t/t H > 1, the blast wave gradually transforms from a sphere into a hemisphere, as energy is redirected back up from the surface. In the long-time limit, all of the blast energy is contained in a hemispherical volume. This long-term hemispherical volume, containing energy Y a could also where Here we have introduced Ω to parametrize the geometric transformation of the blast from a sphere into a hemisphere, implicitly embedding all the intermediate complexities of the Mach stem and triple point associated with the shock-boundary interaction. The Θ parameter accounts for the energy reflected from the surface, including the time it takes for this energy to propagate back through the blast region and reach the perimeter. Since Ω and Θ obey the same limits, it is convenient to combine them into a single parameter At t = t H , the blast volume is The goal of this paper is to determine Ψ by performing computer simulations to directly calculate V.

Simulations
We have computed Ψ numerically using the Miranda code [13] to solve the governing Euler equations [14] for the transport of mass, momentum and energy. The Miranda code discretizes spatial derivatives with a tenth-order compact-finite-difference scheme [15] and temporally integrates the equations with a fourth-order Runge-Kutta method [16]. Eighth-order hyperviscosity [17] and hyperconductivity [18] are employed, along with eighth-order spectrallike dealiasing for shock capturing [19]. The blast volume is computed by tagging grid cells as the shock passes through them and then adding up the volumes of all the tagged cells.
The ideal-gas law assumed by Taylor and employed in our simulations is p = (γ − 1)ρe. and air bursts; i.e. we found this value to produce the best overall agreement between yields estimated from blast-wave measurements and radio-chemical methods. Before collecting results, we first verified that Miranda agrees with Taylor's solution, given sufficient grid resolution. Here we present our grid-converged results obtained at a mesh spacing of H/800. Simulation versus theory for our benchmark case is displayed in figure 2.

Results
A useful means for comparing different cases is obtained by integrating (3.1) over all space (4.1) Since pressure becomes uniformly distributed at a sufficient distance behind the blast, the non-dimensional pressure, pV/Y a , provides a natural measure of correlation for our various simulations. Figure 3 displays the non-dimensional pressure from two different simulations, one at an HOB of 10 m and the other at an HOB of 1 km. In non-dimensional coordinates, the two blasts are seen to be identical, except for tiny differences arising from grid-seeded perturbations of the vorticity field. As the shock strikes the surface at an angle, it generates baroclinic vorticity, which in turn becomes Kelvin-Helmholtz unstable, leading to the formation of vortices with lowpressure cores, seen in the image near ground zero. The perturbations to the vorticity strip are generated by the numerical discretization but the subsequent evolution of the vortices is physical. The volume of the blast, however, is unaffected by these details.
In figure 4, we plot scaled volume versus scaled time for a 100 kiloton detonation at three different HOBs spanning two orders of magnitude. The non-dimensionalized blast volumes are seen to collapse to a single curve. This is not too surprising, since we chose our reference time and volume to guarantee that all curves would pass through V/V H = 1 at t/t H = 1; i.e. changing the HOB can make no difference if time and volume are measured in HOB units.  In figure 5, we compare different yields at the same HOB. Once again the blasts appear identical, even though the yields differ by two orders of magnitude. Changes in yield correspond to changes in blast velocity, but all velocities become equal in the chosen units, thus removing differences due to yield. Tiny variations are observed once more in the vortex cores near ground zero but they do not influence the blast volume. Scaled volume versus scaled time for three different yields is shown in figure 6. This single collapsed curve is identical to the one in figure 4.  Although the yields and HOBs were varied in the simulations, they nevertheless all correspond to the same non-dimensional case. Finally, we varied the background density as an additional test of the theory. Figure 7 compares two simulations with atmospheric densities differing by two orders of magnitude, roughly onetenth of an atmosphere to 10 atmospheres.  The blasts appear identical, as they did with changes in yield and HOB. The volume versus time curves for three background densities are displayed in figure 8. The collapsed curve matches the previous cases. We have thus far demonstrated that for any given yield, HOB and air density, the blast volume is a predictable function of time. In our simulations, the yield is a known input and the blast volume is calculated, however in applications of interest this process is reversed; i.e. the blast volume is measured and used to infer the yield. This becomes a straightforward procedure once the fixed relation between scaled volume and scaled time is established. Solving (2.11) for Ψ , we obtain the following result: where τ ≡ ln t t H .
We have plotted the reflection function (4.2) for all cases in figure 9. The undulations in the curve are due to the Mach stem and triple point working their way up the side of the blast wave.
The curve fit has a correlation coefficient of 0.999345, an RMS relative error of 0.00481881 and a maximum relative error of 0.00756.

Conclusion
A universal reflection function exists for ideal blasts rebounding from perfect surfaces. We have computed this function from numerical simulations of Taylor blast waves with different yields, HOBs and ambient air densities. The reflection function relates non-dimensional blast volume to non-dimensional time, causing all cases to collapse onto a single curve. This curve serves as the basis for estimating the yield of tower shots and other nuclear events, wherein the cameras recorded the expanding blast wave after the shock reached the ground. Using the theory presented in this paper, it is possible to reconcile the EG&G data for Mohawk; i.e. if the reflection function is included in the blast wave analysis then the yield comes into agreement with the light-curve data and both methods give approximately 248 kilotons. This theory can be extended to include real-gas effects, radiation, humidity, gravity and energy-absorbing surfaces.
Data accessibility. All relevant data are plotted in the figures and should be reproducible by other codes. Authors' contributions. A.C. is the principal architect of the Miranda code. He helped formulate the theory and assisted with the numerical simulations. J.B. helped refine the theory, performed code verification and carried out most of the numerical simulations. G.S. originated the idea to extend Taylor's solution to include reflections. He performed extensive film analysis and supplied the historical background material that motivated this paper. All authors contributed towards drafting/revising the paper, approve the final version and agree to be accountable for all aspects of the work.