Analytic computations of nonideal corrections to blast wave overpressure predictions at high altitudes



[1] In this paper, the effects of temperature, pressure, winds, moisture, and molecular content on the propagation of blast waves at high altitudes are investigated. These cause refractions and attenuations which modify the recorded ground overpressures from the ideal predictions. By coupling these effects together, the nonideal corrections to the overpressures are estimated by applying approximations which are dependent on the angle of propagation of the blast wave.

1 Introduction

[2] The phenomenon of nonideal air blast propagation has been known for many years and had been the subject of many investigations over the years [e.g., Gilbert, 1962; Reed, 1972a; Reed, 1972b]. A blast wave is typically characterized by a nearly instantaneous pressure rise along its front as shown in Figure 1.

Figure 1.

Typical blast wave overpressure profile.

[3] Depending on the size of the source of the blast wave, the excess pressure or overpressure can reach values which are multiples of the atmospheric pressure in regions close to the source. Hence, the atmosphere ahead of the shock is heated irreversibly by the passage of the wave and thus travels at a speed greater than the local velocity of sound. In this near-field regime, the blast velocity is dependent on the size of the burst, the proximity to the source and the existing meteorological conditions. However, due to spherical spreading, the overpressure behind the front eventually decreases and reaches a value comparable to that of a local sound wave. Thereafter, in the far-field regime, the blast wave propagates at the sound speed. In this report, we will be concerned with waves in the far-field regime.

[4] Another difference between sound waves and blast waves concerns their corresponding frequencies. Unlike a sound wave which typically has a single frequency, a blast wave is a superposition of waves with different frequencies. Furthermore, in the near-field regime, the blast frequencies change as a function of distance during the expansion of the source. Thereafter, the blast wave resides wholly outside the source products. It can therefore be assumed that in the far-field regime when the blast wave becomes a sound wave, the component frequencies remain constant with distance.

[5] Now the propagation of sound waves in the atmosphere can be mathematically analyzed. This is because the wavefront characteristics allow the use of Huygens' construction for determining future frontal positions. This concept, furnished by Lord Rayleigh, although entirely artificial, allows the location of wavefronts to be constructed using velocity vectors called “sound rays.” This is known as the ray-trace method [Reed and Church, 1963; Gilbert, 1962] the advantages of which include its ease and flexibility of use. It is, however, not possible to use this method to predict the overpressures of diffracted waves. It is also not very accurate for predicting overpressures for low frequency waves which is a characteristic of high-yield explosions.

[6] Typically, blast codes predict blast overpressures using computational fluid dynamics computations in most regions including close to the source [e.g., Rose et al., 2005; Lohner et al., 2004]. These often fail in the far-field or are too dissipative to give meaningful results in this regime. Furthermore, these calculations are often done in an ideal atmosphere characterized by a homogenous temperature, pressure, and with no winds or moisture content. Hence, sound rays propagate in straight lines with overpressures decaying with distance from the source. However, in reality, vertical temperature, wind, and moisture gradients exist which alter the sound speed and cause refractions, reflections, and diffractions of the sound rays. These changes to the sound speeds and directions alter the paths and times of arrival of the waves from that computed in the ideal case and hence alter the predicted overpressures. Also, often omitted in the ideal scenario is the acoustic attenuation of the waves due to the absorption of sound frequencies by the air molecules. This reduces the predicted overpressures. These effects have been observed experimentally in numerous tests [Reed and Church, 1963; Haskell et al., 1952].

[7] In this paper, we will consider nonideal air blast propagation in the far-field regime which is typically characterized by overpressure values of less than about 7 kPa (~1 psi) [Reed and Church, 1963; Haskell et al., 1952]. An analytic method of estimating the effects of the more realistic nonideal conditions on the ideal overpressure computations in the far-field regime is presented. Hence, the fidelity of the blast code computations can be maintained while the nonideal corrections are estimated using the method discussed in this paper. In section 2, the effects of temperature gradients, wind gradients, winds, and moisture content on sound speeds are respectively discussed, and in section 3, the corresponding effects on overpressure predictions are mathematically analyzed. In section 3, the effects of acoustic attenuations are also presented and are coupled to the refractive effects to result in a single equation for which the most accurate solutions are sought and reported. In section 4, an example scenario is considered, and the results are presented and discussed in section 5. A summary of the main conclusions is provided in section 6. In Appendices A and B, more details of the calculations are provided, while in Appendix C, some error analysis is conducted on the mathematical approximations introduced in this paper.

2 Effects of Nonideal Phenomena on Sound Speeds

[8] As mentioned in section 1, nonideal phenomena affect the propagation of blast waves through the atmosphere and hence cause changes to the overpressures predicted in an ideal atmosphere. These effects are considered in turn in the following subsections.

2.1 Atmospheric Temperature

[9] The temperature of the atmosphere varies with height about sea level. The U.S. ANSI standard atmospheric profile [U. S. Standard Atmosphere, 1976] is shown in Figure 2.

Figure 2.

U.S. standard atmospheric temperature [U. S. Standard Atmosphere, 1976].

[10] It can be seen that within the lower 10 km of the atmosphere or troposphere, the temperature decreases with height with a lapse rate of −6.5 K/km. There is a further steeper decrease until at 15 km, where the temperature remains constant for about 5 km. This region of constant temperature is called the tropopause. Above the tropopause is the stratosphere where the temperature increases with height for about 30 km at which point the ionosphere is reached and the temperature again decreases with height.

[11] In this paper, we will focus on blast propagation within the troposphere although the results obtained are readily applied to all regions of the atmosphere. It can be shown from molecular theory the sound speed V in m s-1 in the atmosphere is related to the ambient temperature T in Kelvin by equation (1).

display math(1)

[12] Sound speed increases with temperature and hence variation of atmospheric temperature with height implies a sound speed profile in the atmosphere. Therefore, if we consider for example, a sound ray propagating downward at an angle in the troposphere, from Snell's law of refraction, it will be bent toward the vertical. Likewise, a sound ray propagating upward would be bent away from the vertical.

[13] In an ideal atmosphere, these rays follow straight lines with constant angles, and therefore, there result differences in the path lengths and also ground ranges (distance from where the ray impacts the ground to the vertical). These differences imply that the nonideal overpressures would differ from the ideal values. In section 3, Snell's law of refraction will be used to estimate the raypaths and ground ranges and a method of estimating the corresponding nonideal corrections to the overpressure will be discussed.

2.2 Wind

[14] As was the case with temperature, the presence of winds in the atmosphere alters the sound speed profile. In this case, the profile is no longer isotropic (i.e., the same for all azimuthal angles) as was the case for temperature. Instead, for a given direction, the sound speed Vwind can be obtained by adding the component of the wind velocity w in that direction to its value Vno-wind in the absence of wind as determined by equation (1).

display math(2)

[15] Hence, in the presence of a positive wind gradient in the troposphere (i.e., wind speed increasing with height), the sound rays propagating downward and downwind will be bent further away from the vertical than in the no-wind scenario and vice versa for rays traveling downward and upwind.

2.3 Moisture

[16] The presence of water vapor in the atmosphere also causes changes to the sound speed profile. Water vapor has a lower density than air hence moist air is lighter than dry air (when it condenses, the reverse is true). Since the sound speed is inversely proportional to the square root of the density, then the sound speed increases with the presence of moisture in the atmosphere. By considering the perfect gas equation of state and Dalton's law of partial pressures, it can be shown that the sound speed in moist air Vmoist is given by equation (3),

display math(3)

where Vdry is the sound speed in dry air, ε is the partial pressure of water vapor, and p is the ambient atmospheric pressure. The maximum saturation pressure of water vapor in moist is a function of the ambient temperature and is given by

display math(4)

where T is the temperature in Kelvin.

2.4 Acoustic Attenuation

[17] As explained in section 1, as sound propagates through a nonideal atmosphere, the air molecules absorb certain frequencies which results in an attenuation of the wave and a reduction in its overpressure. The acoustic attenuation of a sound ray is typically parameterized by an attenuation coefficient α per unit distance traveled such that the change in overpressure Δp with distance s is given by

display math(5)

[18] The attenuation coefficient α can be defined as a function of the sound frequency ν and ambient pressure p,

display math(6)

where k is a constant, various values of which abound in the literature. As can be seen from equation (5), the amount of attenuation is dependent on the path traveled by the sound wave. Hence, to determine the reduction in overpressure in nonideal conditions, equation (5) should be integrated along the refracted paths arising from the sound speed profiles originated by the temperature, wind, and moisture gradients already discussed.

2.5 Atmospheric Pressure

[19] As mentioned earlier, in ideal overpressure computations, the atmospheric pressure is held fixed at a constant value (typically the value at sea level p0 = 1.013e5 Pa) when in reality, the atmospheric pressure decreases with height. The nonideal correction to the overpressure is, however, trivial and is given by

display math(7)

where Δpideal is the ideal overpressure computed at the constant atmospheric pressure p0 and Δpnonideal is the nonideal overpressure computed at the nonideal atmospheric pressure p.

[20] Now the atmospheric pressure p in a standard nonideal atmosphere can be approximated by

display math(8)

where b is a constant and h is the height above sea level. Note that this correction to the atmospheric pressure will be applied to the definition of the attenuation coefficient in equation (6).

3 Effects on Overpressure Predictions

[21] Now that some contributing factors to the nonideal propagation of blast waves have been identified, we now turn to ray theory to mathematically quantify the corresponding effects on the calculated overpressures. As was mentioned in section 2, the sound speed profiles arising from temperature, winds, and moisture gradients cause the sound rays to refract in the atmosphere. In section 3.1, Snell's law will be used to determine the refracted paths and an estimate will be made of the effects on the overpressure on the ground. In section 3.2, the frequency-dependent acoustic attenuation along these refracted raypaths would then be calculated and coupled to the refraction correction to give the total nonideal correction to the ideal overpressures.

3.1 Refraction

[22] Consider a blast wave propagating radially in all directions at a height of burst h0 in a sound profile of negative gradient (e.g., the troposphere). In this scenario, rays initially traveling upward would be bent away from the horizontal and rays traveling downward would be bent toward the horizontal. This is illustrated in Figure 3.

Figure 3.

Illustration of Snell's law.

[23] If a ray at an initial inclination angle of θ0 to the horizontal with a speed V0 moves into a region where the sound speed V1, it is bent to an angle θ1 given by

display math(9)

[24] This is Snell's law of refraction. Therefore, in an atmosphere of negative sound speed gradient, the rays propagating upward are refracted further away from the horizontal and never reach the ground. For rays propagating downward, there is a critical initial inclination angle θC below which the ray would never reach the ground but will be refracted back upward. This happens when

display math(10)

[25] As shown in Figure 4, this defines a shadow zone within which the raypath representation we have employed here predicts that the ground overpressures are zero where in reality they would be nonzero due to diffraction processes which have been neglected here. In Section 3.1.2 we will consider how to estimate the overpressures in this zone.

Figure 4.

Shadow zone (not to scale).

[26] Outside the shadow zone, where θ0 > θC, the refracted rays do reach the ground. The horizontal distance r or the ground range from the vertical at which they impact the ground is given by

display math(11)

where β is the sound speed gradient and the subscripts 0 and 1 respectively refer to parameters at the start and end of the raypaths. The total distance s traveled by the ray from the source to the ground is given by

display math(12)

[27] The above analysis can be extended to rays traveling through multiple layers of the atmosphere with varying values of the sound speed gradient β.

3.1.1 Refraction Factor, rf

[28] In order to estimate the change in overpressure at a ground range r due to nonideal refractive effects, the concentrations of the rays on the ground can be determined from both the ideal and nonideal scenario and the ratio used as a measure of the overpressure correction. The concentration of the rays at a ground range r is given by the square root of the derivative of the inclination angle with r (see Appendix A). Therefore, the refraction factor rf can be evaluated as

display math(13)

[29] Expressions for the derivatives are given in Appendix A. Note that similar methods have been used in literature [e.g., Gilbert, 1962; Slater, 1998].

3.1.2 Shadow Zones and Diffractive Overpressures

[30] As we saw earlier, the use of the raypath representation for sound waves leads to shadow zones for inclination angles less than a critical angle θC. At these angles, the rays are refracted back upward before they reach the ground and the predicted overpressures on the ground are therefore zero. However, in reality, overpressures will be recorded in these zones due to diffraction processes which the representation we use here neglect. Experimental data from Reed and Church [1963] show that the recorded diffractive overpressures were correlated with the ground range and vary approximately as

display math(14)

where ΔpC is the overpressure for the ray with inclination angle θC and ground range rC (i.e., at the edge of the shadow zone) and Δpdiff is the overpressure for a ray in the shadow zone at ground range rdiff.

3.2 Acoustic Attenuation

[31] The change in the overpressure due to acoustic attenuation is dependent on the distance traveled in the atmosphere s and is obtained by integrating equation (5). The resulting attenuation factor af is given by

display math(15)

where p and s have already been defined in equations (8) and (12), respectively.

3.2.1 Coupled Refraction and Attenuation

[32] In an ideal atmosphere, the integral equation (15) is analytic since the integration path is linear. However, to couple the refraction and attenuation effects, the definition of s given in equation (12) for the refractive paths can be used. Unfortunately, the resulting indefinite integral does not have a complete analytic solution (see Appendix B). Following standard approaches, the integral can be evaluated numerically. We propose an alternative method which applies approximations in different angular regimes. A full discussion can be found in Appendix B. Here we present the method and the main results.

[33] For initial inclination angles θ0 < 0.7 rad, by making use of the fact that

display math(16)

the integral becomes analytic and can be evaluated to give

display math(17)

where DawF is the Dawson's function and D and E are combinations of variables defined in Appendix B. This solution is the small-angle approximation.

[34] For θ > 1.3 rad, using the approximation

display math(18)

the integral was evaluated to give

display math(19)

where A and B are combinations of variables defined in Appendix B. This will be termed the large-angle approximation.

[35] Note that the 0.7 rad and 1.3 rad were chosen as the limits for the small-angle and large-angle approximations because at these angles, the errors made in the approximations in equations (16) and (18) are less than 1.3%.

[36] For raypaths with 0.7 rad < θ < 1.3 rad, the integral is not analytic so a standard numerical approach will have to be used. This can be achieved by splitting up the paths into linear segments and summing the integral over each segment. This will be termed the straight-line approximation. So for raypath approximated by a straight line between the heights h1 and h0 and constant inclination angle θ0, the contribution to the attenuation factor from this segment is (see Appendix B)

display math(20)

3.2.2 Frequency Decomposition

[37] Figure 1 illustrates the generic form of a blast wave which is typically characterized by a peak value for the overpressure Δppeak and a positive pressure phase duration τ. Scaling relations, obtainable from Glasstone and Dolan [1997] and Reed [1972a], can be used to determine the value of τ for any blast wave given the parameters of a reference blast wave [Reed, 1972a].

[38] Approximating the generic waveform to an isosceles triangle of height Δppeak and base τ, we can express the wave as a Fourier series sum of N waves such that the wave is now represented mathematically as

display math(21)

where the fundamental frequency νf = 1/2τ and hence the Nth wave in the series has frequency N νf.

[39] Now if the computed attenuation factor for a sine wave with a frequency of νf is denoted affund, then it can be shown that after undergoing attenuation, the wave described by equation (21) becomes

display math(22)

[40] Summing this series as done by Reed [1972b] yields the following approximation for the total attenuation factor af in terms of affund:

display math(23)

4 Example

[41] Now that we can compute the refraction and attenuation factors, the total nonideal correction (or coupled factor) cf is simply af • rf. We now examine in detail a specific example for which the aim is to compute cf over a range of inclination angles (i.e., ground ranges). The scenario we consider has the characteristics listed in Table 1. The value of k used was obtained from Reed [1972b] while the values of b and p0 are those attributed to standard atmospheres. Also presented in Table 1 are some derived variables and results. These are shown in bold.

Table 1. Example Scenario Used in Section 3
Yield of burst (kg HE)0.5
Height of burst (km)10
Moisture contentMaximum vapor pressure (see equation (4))
Wind speed (direction θ0 = 0 rad) (m s-1)35
p0 (Pa) [U. S. Standard Atmosphere, 1976]1.013e5
b (km-1) [U. S. Standard Atmosphere, 1976]0.119
k (kg m-2) [Reed, 1972b]2e-6
νf (Hz)61.9
Critical angle θC (rad)0.48
Range at edge of shadow zone rC/km40.9

[42] Using these parameters and the equations discussed in preceding sections, the refraction, attenuation, and total coupled factors were computed and are shown as a function of initial inclination angle in Figures 5, 6, and 8, respectively. Note that all these results are for the downwind direction although as highlighted later there is not much difference to the upwind results.

Figure 5.

Refraction corrections, rf.

Figure 6.

Attenuation corrections, af.

[43] In Figure 6, the results for the three attenuation approximations are displayed and these results are combined with the corresponding refraction factors in Figure 5 to give the total coupled corrections in Figure 8. Note that for the straight-line approximation, the path was split into 40 segments of equal length, i.e., 250 m. Figure 9 shows, as a function of ground range, the ideal peak overpressures as estimated using Brode's formula [Brode, 1987], with the nonideal peak overpressures predicted by our method.

5 General Discussion

5.1 Trends

[44] From Figure 5, it can be seen that as the inclination angle decreases from π/2, the refraction factor decreases slowly from a value close to 1 to about 0.75 at the edge of the shadow zone beyond which it decreases rapidly as would be expected of the diffracted waves which contribute in this region. It should be noted that most of the refractive correction is due to the temperature profile, the effects of the winds and moisture being smaller in comparison.

[45] Figure 6 shows the three approximations used for the attenuation computations for angles θ0 > θC (for the shadow zone, the paths of the diffracted waves are assumed to be the ideal paths, so the straight-line approximation is used here). The results from the straight-line approximation are shown over all the considered angles for comparison with the analytic approximations. For the straight-line approximation, the integration path was split into 40 segments. Figure 7 displays the results from straight-line integrations of 10, 20, and 40 segments in comparison to the small-angle approximation. As can be seen, the straight-line approximation converges as the number of segments is increased.

Figure 7.

Convergence of straight-line approximation.

[46] Figure 8 shows the combination of Figures 5 and 6 to give the total coupled correction factor, while in Figure 9, a comparison is made of the ideal overpressures (as calculated using Brode's formula) with the nonideal values computed using the computed corrections. This shows the large discrepancy expected in the shadow zone (beyond a ground range of about 40 km) as a result of the nonideal diffractive processes which are approximately accounted for in this region. Note that the linear relationship between the logarithms of the overpressures and the ranges in the shadow zone as shown in Figure 9 is as expected from the crude inverse square law approximation used in this region (see equation (14)). Further work in this field will attempt to determine a better correlation for the nonideal overpressures in this zone.

Figure 8.

Coupled corrections, cf.

Figure 9.

Ideal and non ideal overpressures as a function of ground range.

5.2 Error Analysis

[47] In practice, it is possible to make the straight-line approximation more accurate by increasing the number of segments over which the summed integrations are carried out. Though this may be computationally intensive (especially if a multitude of paths over long ranges are to be computed), it is worthwhile investigating the number of segments required to make the straight-line approximation approach the values predicted by the two analytic approximations.

[48] In Appendix C, this is done by estimating the errors in the integrals resulting from the three approximations. From this, it is then possible to determine the number of line segments required for the straight-line approximation errors to be commensurate with those of the small-angle and large-angle approximations.

[49] For the scenario considered in section 4, it was found that this occurs for both the small-angle approximation and large-angle approximation when the segment lengths are less than 5 m (i.e., greater than 20,000 segments over a 10 km path).

6 Conclusions

[50] The corrections to ideal ground overpressure predictions have been estimated by the consideration of the nonideal conditions arising as a result of variations in atmospheric temperature, pressure, moisture content, and winds. These computations were applied in the far-field regime where blast waves approximate sound waves. Both atmospheric refractions and acoustic attenuations were modeled and their effects coupled together for a computation of a composite total nonideal overpressure correction. In two different angular regimes, analytic approximations were applied to obtain the most accurate solutions with a standard numerical approach applied elsewhere. Finally, the errors in the approximations were estimated and compared.

Appendix A: Estimating the Refraction Factor

[51] In Figure A1 are shown the raypaths for downward propagating rays in both ideal and nonideal conditions, the former being straight lines and the latter being curved refracted paths.

Figure A1.

Raypaths in ideal and nonideal atmospheres.

[52] Now consider a beam of rays of angular width around a ray with an initial inclination angle θ. On the ground, this beam would describe the concentric area shown in Figure A2, the area dA of which is given by

display math(A1)

where for the ideal scenario θ and s are, respectively, the inclination angle and length of the straight line connecting the source to where the beam impacts the ground.

Figure A2.

Differential ground area, dA.

[53] To compare the areas dA at a given ground range r, the values for θ and s necessarily coincide for the ideal and nonideal scenarios, and hence the ratio of the areas is equivalent to the ratios of the differential ground ranges dr. Now since the overpressure is the square root of the intensity and intensity is inversely proportional to the beam area, therefore an estimate of the refraction correction is given by

display math(A2)

Appendix B: Estimating the Attenuation Correction

[54] The attenuation factor af was defined in equation (15) as a function of an integral over the path length s. This is presented again in equation (B1).

display math(B1)

[55] Using the expression for s in equation (12) to perform a change of variables to the inclination angle θ

display math(B2)

and using the definition given for p in equation (8),

display math(B3)

as well as the expression for the altitude h,

display math(B4)

we obtain the following expression:

display math(B5)

[56] Now this integral is of the form

display math(B6)

where A and B are constants and are given by

display math(B7)

[57] For values of θ < 0.7 rad

display math(B8)

[58] This reduces equation (B6) to the form

display math(B9)

where D and E are given by

display math(B10)

[59] Equation (B9) has an analytic solution of the form

display math(B11)

where DawF is the Dawson's function or integral, mathematical tables, and packages for which exist. We shall term this solution the small-angle approximation.

[60] For θ > 1.3 rad, if we use the approximation,

display math(B12)

the integral becomes

display math(B13)

which can be solved to give

display math(B14)

[61] This will be termed the large-angle approximation.

[62] For raypaths with 0.7 rad < θ0 < 1.3 rad, the integral is not analytic so a standard numerical approach will have to be used. This can be achieved by splitting up the paths into linear segments and summing the integral over each segment. This will be termed the straight-line approximation.

[63] For a raypath approximated by a straight line between the heights h1 and h0 and constant inclination angle θ0, the path length s in equation (12) is approximated by

display math(B15)

and the integral (now over the height h) becomes

display math(B16)

which when solved gives

display math(B17)

Appendix C: Error Analysis

[64] In order to estimate the errors in the computations resulting from the three approximations, the form of the integral in the exponent of equation (15) was investigated. After a change of variables to the inclination angle θ, this becomes of the form

display math(C1)

[65] In the small-angle approximation, the error in the integrand is

display math(C2)

[66] Likewise, the error in the integrand for the large-angle approximation is

display math(C3)

[67] And for the straight-line approximation where the inclination angle is fixed at θ1, the corresponding error is

display math(C4)

[68] Since in the scenario considered in section 4, θ1 < θ0, it follows that, for all three approximations, this error is a maximum for θ = θ0. Denoting this as errmax and the approximated integral solution as Iapprox, the maximum error in the integrals ERRmax may be computed as shown in equation (C5).

display math(C5)

[69] This can be calculated in turn for each approximation by inserting the appropriate value for errmax. In this way, the estimated error in the straight-line approximation can be compared against the analytic approximations, and the number of line segments required for similar accuracy can be determined.


[70] The author would like to acknowledge the valuable help received from Tim Rose, Mike Kerry, Jig Atwal-Patel, Rhys Edwards, Paul Mattingley, and John Adams during the course of this study. This work was funded by AWE Plc, Aldermaston, Reading, Berkshire RG7 4PR, UK.