Blast wave kinematics: theory, experiments, and applications

Measurements of the time of arrival of shock waves from explosions can serve as powerful markers of the evolution of the shock front for determining crucial parameters driving the blast. Using standard theoretical tools and a simple ansatz for solving the hydrodynamics equations, a general expression for the Mach number of the shock front is derived. Dimensionless coordinates are introduced allowing a straightforward visualization and direct comparison of blast waves produced by a variety of explosions, including chemical, nuclear, and laser-induced plasmas. The results are validated by determining the yield of a wide range of explosions, using data from gram-size charges to thermonuclear tests.


I. INTRODUCTION
Recent large-scale industrial accidents such as those in Tianjin (2015; 173 deaths) and Beirut (2020; 218 deaths) provide stark illustrations of the devastating potential of explosions. In addition to the tragic loss of human lives, the latter caused an estimated $15B in property damage: complete destruction of buildings extended to a few hundred metres from the source of the explosion, and broken glass and debris was observed at distances up to 3 km from the explosion center, encompassing an area with more than 750,000 inhabitants [1]. Clearly, in order for engineers to design structures for resilience against explosions, the properties of the blast wave must be known both relatively close to (where the structure should be designed to avoid/limit progressive and disproportionate collapse) and relatively far from the source (where the majority of injuries are caused by either lacerations from airborne glass fragments or by damage to hearing from failed glass panels [2]).
Knowledge of the arrival time of a blast wave at various distances from the source enables a radius-time relationship to be developed, from which other key parameters such as peak pressure can be derived [3,4]. Thus, the ability to determine this relationship a priori, from a known explosive yield, will provide vital information on the properties of the blast wave as it propagates. Further, a well-defined relationship that is valid for any distance permits the yield of an explosive to be determined through inverse analysis [5].
This article presents a description of the propagation of a shock wave produced by an explosion in free air, an extension of the standard strong-shock solution to its later phase transitioning into an acoustic wave, and the applications of the results for estimating the yield of a wide variety of explosions as well as the method is outlined for its future application.

II. THEORETICAL DESCRIPTION OF THE BLAST WAVE
Let us model the shock wave produced by an explosion in free air as a sphere of time-dependent radius R. A reflection factor can be used to extend the results in this section to explosions in the vicinity of surfaces and those produced by hemispherical charges. The energy E 0 of the explosion is assumed to be released instantaneously and in a minuscule volume in air of undisturbed ambient conditions of atmospheric pressure P 0 and density ρ 0 . Conservation of energy and the equation of state of an ideal gas can be used to write the energy released by the explosion in terms of the kinetic and thermal energy of the gas contained within a radius R in the form where r represents a radial coordinate measuring the distance from the center of the explosion to the shock front R. The factor γ is the heat capacity ratio, assumed to be unaffected by the passing of the shock; its value for air in normal conditions described as a diatomic gas is γ = 1.4. The radial velocity u, pressure P , and density ρ of the air behind the shock front satisfy well-known hydrodynamics equations, which must be solved to determine their radial dependence before performing the integration in (1).
The PDE system describing the motion, continuity, and equation of state of the fluid are respectively given by subject to the boundary conditions at the shock front given by the Rankine-Hugoniot relations. In terms of the Mach number of shock front M S = a −1 0 (dR/dt), these relations are where a 0 = γP 0 /ρ 0 is the speed of sound at ambient conditions. Let us characterize the motion of the shock front by introducing the dimensionless variables where η specifies the distance from the explosion center (η = 0) to the shock front (η = 1); whereas λ characterizes the speed of the shock front from high Mach number (λ → 0) to the ambient speed of sound (λ = 1). Let us now write the ratios of the three quantities of interest in terms of the new variables as Note that for the strong-shock regime (M S 1) Taylor's definitions [6] are recovered. Using the functions (9), the energy equation (1) can be rewritten as where we have introduced the dimensionless scaled distance z = R/R 0 , which measures distance in units of the explosion characteristic length R 0 = (E 0 /P 0 ) 1/3 . Notice that z differs from the standard scaled distance Z = R/W 1/3 used in blast engineering; the latter normalizes the distance by the cubic root of the mass W of the explosive charge, whereas z removes cumbersome units and, more importantly, eliminates sometimes problematic TNT equivalence of different explosive materials. The function K(λ) is defined as In the limit R → ∞ the blast wave decays to an acoustic wave (λ → 1), hence the constant K 1 ≡ K(1) in (10) corresponds to the boundary value of K in the far field. The decay of the blast wave can be parametrized by the auxiliary function that describes how the speed of the shock front decreases as it moves away from the explosion center. From the energy equation (10), it follows that this auxiliary function and K(λ) are related by the ordinary differential equation This relation implies that the auxiliary function must satisfy the boundary conditions ζ(0) = 1, ζ(1) = 0. The simplest description of the blast decay that allows for an analytical description of the Mach number of the shock front and satisfies the boundary conditions is the linear decay ζ(λ) = 1−λ. Numerical analysis and experimental observations suggest that the decay is nonlinear, with the Mach number decaying more rapidly at early times. In this work we intend to provide an approximate description of the phenomena; therefore, the linear choice will suffice. Since the boundary conditions are satisfied, our approximate description will match the exact solutions in the early and late regimes, whereas some small deviation can appear in the mid-range where the strong shock transitions to the acoustic wave. In Sec. IV we will see that the linear ansatz provides an accurate description of the blast wave for all ranges. The linear form of the auxiliary function ζ(λ) leads to a simple solution of (13) given by where the integration constant has been chosen so that K 0 denotes K(λ) evaluated at λ = 0. This solution allows inverting the energy equation (10) to write the the Mach number in terms of the scaled distance as where we have introduced the dimensionless scaled time τ = a 0 t/R 0 . Another reason for using these dimensionless variables (τ, z) is that they allow direct comparison of a wide range of experiments independent of the yield of the explosion under consideration. This enables us to visualize the results from gram-sized charges to megaton yields from thermonuclear explosions in the same plot, as is done in this article. Notice that λ in (8), the auxiliary function (12), and its linear form can also be used to write a nonlinear differential equation for M S (z). This equation is of the Bernoulli type so that it can be analytically solved; its solution is again given by (15), which confirms the mathematical self-consistency of the system. The solution for the Mach number (15) shows that only the numerical value of the function K(λ) (11) at λ = 0 is needed for fully describing the propagation of the shock front. This observation in turn implies that the solutions of the hydrodynamics functions φ(η, λ), ψ(η, λ), and f (η, λ) are necessary only at λ = 0, which significantly simplifies the ODE system (2)(3)(4). Using the definitions (9), the solution to the system (2-4) with boundary conditions given by the Rankine-Hugoniot relations (5-7) at λ = 0 is where the exponents κ i , i = 1, . . . , 4 are only functions of the heat capacity ratio γ: The three functions in terms of the dimensionless radial coordinate are shown in Fig. 1. We can now use these solutions in the definition of K(λ) to determine K 0 in the form where the heat capacity ratio for air has been used since we have assumed the explosion to take place in free air. Once this value is determined, the Mach-number equation (15) can be used to describe the growth of the spherical shock front as a function of the distance from the explosion center. The general solution of (15) is shown in Fig. 2 together with the strong-shock solution discussed in Sec. III and the acoustic wave that the general solution must asymptotically approach. Given the analytical form of the Mach number (15), the Rankine-Hugoniot relations can be used to write a simple expression for the peak hydrostatic overpressure behind the shock front as where the last form is relevant for chemical explosions. The energy of the explosion has been related to the mass of a charge by E 0 = E m W , where E m is the specific energy per unit mass that characterizes the chemical energy converted into kinetic and thermal energy after the explosion. For example, considering a TNT explosion (E m ≈ 4.3 MJ/kg) the oversimplified expression (21) leads to an overpressure barely distinguishable from the Brode formula for spherical blasts [9].

III. SEDOV-TAYLOR-VON NEUMANN BLAST WAVE
The famous Sedov-Taylor-von Neumann solution [6-8] assumes a strong shock (P P 0 ), which corresponds to setting λ = 0 and neglecting the thermal energy of the air before the explosion. This is equivalent to solving the blast-wave equation (15) for the early stages of the explosion when z 3 K −1 0 , simplifying the Mach number equation to the reduced form whose solution is shown in Fig. 2 as a straight line of slope 2/5 in the loglog plane. In standard coordinates, we recover the more familiar form whose solution is the well-known STvN blast wave The constant factor for air is which is moderately closer to the exact value S(1.4) = 1.033 than the approximate result S(1.4) = 1.014 found by Chernyi [10]. It should be emphasized that this description of a blast wave is only valid in the early stages of expansion and where the explosion can be assumed to originate as point-source energy release, such as a nuclear explosion or in the mid-range for a chemical explosion. In a later stage, a blast wave will decay and the strongshock approximation will no longer be valid (and in the early stages of a chemical explosion the energy release will not be from a point-source). For a full description of the blast wave, and more crucially, including the transition from a string shock to an acoustic wave we must solve the equation for the general Mach number (15). As shown in (22), the STvN solution is obtained when neglecting the thermal energy of the undisturbed air before the explosion via the strong shock condition (P P 0 ). Similarly, by comparing the general differential equation (15) describing the blast wave and the STvN limit (22), we can write an upper value for the validity of the STvN solution from the general expression (15) in the form For scaled distances higher than z upp deviations from the STvN solution are expected due to the decay of the shock wave. This behavior is independent from the type of explosion; chemical or nuclear.
In the other direction, there is also a lower value z low for the range of validity of the STvN solution for chemical explosions. The solution neglects the mass of the explosive charge W compared to the mass of the surrounding air over which energy has to be transferred. For this reason, there is a minimum distance from the center of the explosion where the mass of the charge can no longer be neglected. Imposing the condition m air 2W , we find where E m is the specific energy per unit mass introduced in the previous section. In standard dimensions, the range of validity of the STvN solution can be written in the form Blast-wave data of the Trinity test. As described by Taylor [28], the fireball data follows the STvN solution. The measurements reported by witnesses of the test from different locations follow the curve in the acoustic regime.
Using scaled distance Z = R/W 1/3 , the range of validity of the STvN solution in air becomes where the specific energy per unit mass E m must be in MJ/kg. Explicit values for TNT, PE4, and ammonium nitrate are presented in Table I.

IV. EXPERIMENTS
As mentioned in the previous section, a very nice property of the dimensionless scaled coordinates (τ, z) is that we can visualize explosions from multiple different yields in a single plot. In this section we consider measurements of the arrival time of the shock front at different distances for a variety of explosions and show how these measurements agree with the results from the previous sections. For nuclear explosions, the units kt and Mt refer to 10 3 and 10 6 tons of TNTe, respectively.

A. Gram-sized explosive charges
The explosion of gram-sized charges offer the possibility of studying the very early stages of a blast as well as the influence of different charge geometries. High-speed cameras allow for recording of the early shock wave and sensitive devices can measure the overpressure without being destroyed by the blast. In recent years, researchers at the University of Sheffield (UoS) Blast and Impact Laboratory have conducted approximately 80 far-field arena tests using hemispheres of PE4 explosive [13][14][15][16][17], and a smaller number of near-field tests using spheres of PE4 [18,19]. The results are shown in the top-left panel of Fig. 4. For comparison, the figure also includes the curve of the ConWep data for 1 kg of TNT [20].

B. Large chemical explosions
Many tests of significant amounts of explosives have been carried out using TNT and ANFO to mimic the effects of kiloton-range nuclear explosions [21,22]. Accidental explosions, such as the Beirut blast [23][24][25], also allow for studies in this range. The data for a selection of explosions in this range is shown in the top-right panel of Fig. 4.

C. Early nuclear explosions
From the first nuclear test (Trinity), nuclear explosions with yields in the dozens of kilotons were abundant during the late 1940s through to the 1950s. Many unclassified technical reports of these tests include information of the pressure measurements at different distances from ground zero [26,27]. In particular, Trinity is the only test for which early data is available and this is in fact what G.I. Taylor used in his second paper [28]; however, the far-field data is missing. General Leslie Groves requested many firsthand accounts describing the reactions of peo-ple who witnessed the Trinity test [29]. The reports by the scientists include information of their location and arrival time of the blast wave, which we have used to map the evolution of the Trinity blast in the far-field region shown in Fig. 3. For all later nuclear tests only midto far-field data is available, whereas early-time measurements at millisecond scales remain unpublished. A team of scientists, historians, and filmmakers at Los Alamos and Livermore National Laboratories are currently working on the restoration and digitization of old nuclear-tests films and it is expected that fireball data will be published in the near future [30].

D. Laser-induced shock waves
Shock waves can be generated by the fast deposition of energy in different materials by laser pulses. A second laser can be used for diagnostics of the produced plasma and some of its properties can be inferred by studying the time evolution of the shock as well as the plasma plume in a variety of geometries. These laser-induced shocks are usually characterized by the STvN solution [31]. Data of a spherical shock produced by a joule-range laser is included in the bottom-left panel of Fig. 4.

E. Thermonuclear explosions
During the Cold War the development of advanced nuclear weapons pushed the yield from kilotons to megaton thermonuclear tests [32][33][34]. The formidable amount of energy released by these explosions allow for reliable measurements only very far from ground zero; however, the high yields lead to short scaled distances and times into the mid-field region. Results from a selection of thermonuclear test are shown in Fig. 4.

V. APPLICATIONS
One useful application of the results of the previous sections is the determination of the yield, E 0 , of an explosion from a set of (t, R) pairs. It is tempting to simply fit the solution of Eq. (15) to data; nonetheless, there are a few considerations to keep in mind to avoid falling into conceptual traps: 1. One aspect to take into account is the behavior of the curve at different ranges. As shown in Fig. 2, the solution coincides with the STvN line in the short range, meaning that for very early times and short distances the solution might fail to properly describe a chemical explosion; this is not an issue for nuclear explosions, as mentioned in Sec. III.
2. Additionally, the solution in the long range asymptotically approaches the acoustic wave (M S → 1) independent of the energy E 0 . This feature translates into a highly degenerate solution, making the use of long-range-only data unreliable for determining E 0 . This degeneracy is broken in the short range, and for this reason short-range data is crucial for a reliable determination of E 0 .
3. Fitting the solution of Eq. (15) to data in the (t, R) space makes the analysis highly sensitive to the values of the long-range data, where large uncertainties can render the analysis useless. Furthermore, the breaking of degeneracy described in the previous point is negligible on a linear scale. Instead, the fit ought to be carried out in the (log τ, log z) space, where the log-log scale eliminates the problems from the linear scale. robust method is obtained by using emcee, a Python implementation of the affine-invariant ensemble sampler for Markov Chain Monte Carlo (MCMC) [35,36]. Using the combined data sets from Refs. [24] and [25], the resulting posterior probability distribution of the model parameter E 0 is shown in Figure 5. The value E 0 = 514 +41 −43 ton TNTe represents the median of the distribution and the uncertainties are based on the 16th and 84th percentiles of the sample. This value accounts for the fact that the Beirut explosion took place at ground level rather than in free air (assumed in previous sections). The correction is obtained by dividing E 0 by the reflection factor 1.8 (for soil), to account for the enhancement of the shock wave due to the ground-reflected hemisphere and the energy loss due to cratering and ground shock [37]. However, due to the built-up nature of the Port of Beirut and its surroundings, a factor of 1.8 is deemed appropriate. For explosions near sea level, during nuclear tests on the Pacific Proving Grounds it was found that the reflection factor is closer to 1.6 due to extra energy dissipation in the form of large water displacements [38]. Figure 6 shows the data and the corresponding scaling using the value E 0 = 514 ton TNTe. Note that the first plot shows the (log τ, log z) space so the curves are scale independent, whereas the data is scaled. On the contrary, the second plot shows the (log t, log R) space, in which the curves rather than the data are scaled.
The value of E 0 found above is in excellent agreement with other yield determinations of the Beirut explosion, including those performed independently by each of the present authors [24,25] and others. As a more general validation of the method, we have applied it for estimating the yield of a selection of high-explosives and nuclear tests over a wide range of energies, from a few tons of TNTe to the high yields of thermonuclear tests [21,26,27,[32][33][34][39][40][41]. The results of the fit of E 0 for thirteen historical explosions are shown in Figure 7, where the fits are compared to the actual yield in tons of TNTe. As indicated earlier, the accuracy of the parameter fit relies on the availability of data in the early stages. Similarly, the precision of the parameter fit depends on the noise level of the data set. These features are noticeable in the figure for the noisiest data sets corresponding to the tests Bee (Operation Teapot) and Harry (Operation Upshot-Knothole).
The excellent agreement between the fitted and actual yields over several orders of magnitude confirms that (15) provides an acceptable description of the shock front, despite the unrefined approximation of a linear decay of the Mach number. We remark in passing that an evident deviation from the exact description of the Mach-number appears as the decay into the acoustic regime according to (15) does not include the logarithmic dependency found both theoretically [42] and semi-empirically [43].

VI. SUMMARY AND CONCLUSIONS
This article illustrates the results of a general characterization of a blast wave in free air, extendable to other configurations by using a reflection factor. A linear ansatz for the decay of the Mach number of the shock front as it expands allows for analytical solutions of the hydrodynamics equations that lead to a concise expression for the Mach number of shock front in terms of the distance from the explosion center. Despite the unsophisticated approximation for the deceleration of the shock front, the subsequently obtained expressions show an excellent agreement with experimental data.
A simple formula for the Mach number was derived in the form of an ordinary differential equation, whose solution describes the position versus time development of the shock front. Here is where time of arrival measurements can be used for estimating the energy released E 0 , a crucial parameter that determines the shock evolution and the loading developed on obstacles with which it interacts. The general solution found contains the wellknown strong-shock solution as a limit in the early stage of the shock development, beyond this regime the solution describes the transition to an acoustic wave in the far-field. Experimental data from gram-sized charges was used verify the validity of the results and later archival data from large-scale explosions was also employed using dimensionless coordinates for time and distance so that explosions from grams of PE4 to thermonuclear blasts can be visualized in a single diagram. The solution found serves as a generalization of other descriptions of the decaying blast wave, in this case, valid from the early (strong) stage to the asymptotically acoustic behavior at the far field.
A discussion about the validity of the strong-shock solution was presented that can serve a valuable resource for blast engineers. The yield of over a dozen explosions was estimated as way to validate the results found in this work and the main aspects of a fit to time-of-arrival data are discussed. Our results show that one of the key features when fitting the yield to time-of-arrival data is that this must be performed in a log-log space; otherwise, slight errors in far-field data will dramatically affect the estimate of E 0 and can possibly render the analysis useless. This property is due to the highly degenerate nature of the blast-wave solution in far field, where all solutions asymptotically approach to an acoustic wave independent of the yield E 0 . Additionally, when applied to laser-induced shocks, the method outlined in this work becomes a direct diagnostic of the laser energy deposited in the material.