The propagation of planar decaying shock and blast waves in non-uniform channels is investigated with the use of a two-equation approximation of the g...

20 downloads
2241 Views
5MB Size

No documents

ORIGINAL ARTICLE

On the propagation of decaying planar shock and blast waves through non-uniform channels J. T. Peace1

· F. K. Lu1

Received: 11 October 2017 / Revised: 19 February 2018 / Accepted: 13 March 2018 © Springer-Verlag GmbH Germany, part of Springer Nature 2018

Abstract The propagation of planar decaying shock and blast waves in non-uniform channels is investigated with the use of a twoequation approximation of the generalized CCW theory. The effects of flow non-uniformity for the cases of an arbitrary strength decaying shock and blast wave in the strong shock limit are considered. Unlike the original CCW theory, the twoequation approximation takes into account the effects of initial temporal flow gradients in the flow properties behind the shock as the shock encounters an area change. A generalized order-of-magnitude analysis is carried out to analyze under which conditions the classical area–Mach ( A–M) relation and two-equation approximation are valid given a time constant of decay for the flow properties behind the shock. It is shown that the two-equation approximation extends the applicability of the CCW theory to problems where flow non-uniformity behind the shock is orders of magnitude above that for appropriate use of the A–M relation. The behavior of the two-equation solution is presented for converging and diverging channels and compared against the A–M relation. It is shown that the second-order approximation and A–M relation have good agreement for converging geometries, such that the influence of flow non-uniformity behind the shock is negligible compared to the effects of changing area. Alternatively, the two-equation approximation is shown to be strongly dependent on the initial magnitude of flow non-uniformity in diverging geometries. Further, in diverging geometries, the inclusion of flow non-uniformity yields shock solutions that tend toward an acoustic wave faster than that predicted by the A–M relation. Keywords Shock wave · Blast wave · CCW theory · Flow non-uniformity

List of symbols a Sound speed (m/s) A Channel cross-sectional area (m2 ) Starting channel height (m) h0 M Mach number n Geometric index p Pressure (N/m2 ) ∂ p/∂t + ρa∂u/∂t (kg/m s3 ) Q1 ∂ Q 1 /∂t (kg/m s4 ) Q2 t Time coordinate (s) Communicated by D. Zeitoun and A. Higgins. An abridged version of this paper was presented at the 22nd International Shock Interaction Symposium, July 4–8, 2016, Glasgow, UK.

B 1

J. T. Peace [email protected] Aerodynamics Research Center, Mechanical and Aerospace Engineering Department, University of Texas at Arlington, Arlington, TX 76019, USA

u W x α β γ θ ρ τ Subscripts 0 s

Velocity (m/s) Shock wave velocity (m/s) Axial coordinate (m) Coefficient function Coefficient function Ratio of specific heats Channel half-angle Density (kg/m3 ) Time constant of flow property decay (s)

Initial or undisturbed gas state Property on the shock

1 Introduction The propagation of shock waves through non-uniform channels with varying cross-sections has been studied analytically from the 1950s. Notably, Chester [1] provided a linearized

123

J. T. Peace, F. K. Lu

solution to the problem of a transmitted shock wave in a channel of non-uniform cross-sectional area, specifically, channels that feature monotonically increasing or decreasing sections. The solution allowed for an approximate determination of the change in shock strength over a finite length of channel of arbitrary shape. A significant result was that the shock strength, averaged at any time over the flow crosssectional area, is proportional to the change in area of the channel. A similar theoretical treatment by Chisnell [2] yielded a closed-form approximate expression for the change in shock strength as a function of channel area by using a steadystate analysis. This was achieved by neglecting the reflected disturbances generated by the shock. The final major theoretical treatment was that of Whitham [3]. Whitham’s theory allowed for the computation of the shock motion without directly determining the flow field following the shock. The shock motion was determined by integrating the Rankine–Hugoniot (RH) shock jump conditions on the forward propagating characteristic, which yielded an explicit relationship for the channel area and the local shock Mach number. The shock strength and channel area relation obtained is exact with those of Chester and Chisnell. Whitham’s theory has been widely used in analyzing the motion of shock waves in various geometries with accurate results [4–7]. Collectively, the work of Chester, Chisnell, and Whitham is known as the CCW theory. It is noted that the results of CCW theory provide the basis for formulating Whitham’s theory of geometrical shock dynamics [8]. A comprehensive review and comparison of geometrical shock dynamics with kinematic models for shock wave propagation in various geometries can be found in [9]. In the case of blast waves, or shock motion with a nonuniform flow field following the shock, the application of CCW theory is not entirely appropriate. This is due to CCW theory neglecting gradients in the flow field just behind a shock wave, which strictly prohibits its application to initially uniformly propagating shocks, or commonly referred to as freely propagating shocks. In regard to this aspect, Best [10] provided a theoretical framework for the generalization of CCW theory. Best reconsidered the motion of a shock through a channel of varying cross section and described it by an infinite sequence of ordinary differential equations. Truncation was used to show the varying degree of approximation between the original CCW theory and the generalized CCW theory with the inclusion of higher-order terms for cylindrical and spherical shock motion. In this manner, a non-uniform flow field following the shock can be taken into account. Best applied the generalized theory successfully to the propagation of underwater blast waves. The theory was shown to have excellent agreement with the approximate Kirkwood–Bethe method [11].

123

A0

Ms

A(x)

(a)

A0

Ms

A(x)

(b) Fig. 1 a Planar shock wave traveling with a shock Mach number of Ms entering converging channel and b diverging channel

The focus of the current work is to understand the propagation of decaying planar shock and blast waves incident on non-uniform channels using the generalized theory of Best. This gas dynamic process is shown schematically in Fig. 1 for both converging and diverging planar geometries. It is desired to analyze the case of an arbitrary strength shock that is being overtaken by an unsteady rarefaction wave. For this particular wave process, the overtaking rarefaction creates a non-uniformity in the flow behind the shock that can largely influence the shock wave propagation dynamics as it traverses through a change in channel area. The original CCW theory is unable to account for this non-uniformity that is present in decaying shock and blast wave profiles. For the case of a decaying shock wave, the shock profile has been analytically described by Sharma et al. [12] and will be used in conjunction with the generalized CCW theory. Similarly, the classical Sedov–von Neumann–Taylor (S–vN– T) [13–15] blast wave will be analyzed with the generalized theory. The present analysis will enable the description of blast wave dynamics for arbitrary strength shocks and in the strong shock limit. The dynamics of such planar waves in non-uniform channels have yet to be discussed in the literature when accounting for the non-uniform flow behind the wave. Lastly, the results will be compared with the results of planar shock waves to demonstrate the difference in shock and blast propagation through non-uniform channels.

2 Governing equations 2.1 Area–Mach relation The area–Mach (A–M) relation is a classical result obtained by Chester [1], Chisnell [2], and Whitham [3] from linearization of the governing equations of inviscid and quasi one-dimensional flow with no heat transfer or mass addition.

On the propagation of decaying planar shock and blast waves through non-uniform channels

0.4

Substituting the RH shock jump conditions on the forward propagating C+ characteristic yields the ordinary differential equation 1 A (x) dM =− dx g(M) A(x)

0.3

(1) 0.2

where

g(M) =

M 2 1+ −1 γ +1

M2

1 − μ2 μ

1 + 2μ +

1 M2

0.1

0.0

(2) and μ2 =

1

2

3

4

5

6

7

8

9

10

Fig. 2 Measure of coincidence between C+ characteristic and shock in x–t space

(γ − 1) M 2 + 2 2γ M 2 − (γ − 1)

(3)

It is noted that through linearization, the accuracy of the above result is dependent upon the condition that [10] A (x) 1 A(x) ρa 2 u

1 − a + u ∂ p + ρa ∂u a0 M ∂t ∂t

p(x)

p(x)

|Q1| = 0 |Q1| > 0 Ms

(4)

where p, u, ρ, and a represent the post-shock state and are given by the RH shock jump conditions. The above inequality necessitates that the area change dominates the motion of the shock wave compared to the unsteadiness in pressure and velocity following the shock. On the right-hand side of (4), the term |1 − (a + u)/(a0 M)| is a measure of the coincidence between the C+ characteristic and the shock wave in x–t space. This term is exactly equal to zero for M = 1, and approaches 0.309, 0.274, and 0.261 as M → ∞ for γ = 5/3, 7/5, and 4/3, respectively. The behavior of this term is shown in Fig. 2. Similarly, the term |∂ p/∂t + ρa∂u/∂t| is a measure of the non-uniformity of the flow behind the shock and is exactly zero for a freely propagating shock, namely, uniform planar shock motion. This type of shock motion is shown schematically in Fig. 3a. However, as the shock moves through an area change, the changing shock strength perturbs the immediate flow behind the shock such that |∂ p/∂t + ρa∂u/∂t| = 0. Han and Yin [16] used Chester’s small perturbation theory to demonstrate that the term is usually very small along the C+ characteristics originating from the uniform region behind the shock. Therefore, it is these conditions that justify the assumptions of (1), and make the A–M relation a good approximation for describing the dynamics and strength of an initially uniformly propagating shock wave as it traverses through a non-uniform area change.

(a)

Ms x

(b)

x

Fig. 3 a Pressure profile of uniformly propagating shock wave and b decaying shock wave with flow non-uniformity

2.2 Higher-order equations For the cases of decaying shock and blast wave propagation, such that the wave is continuously being overtaken by an unsteady rarefaction wave as shown in Fig. 3b, the term |∂ p/∂t + ρa∂u/∂t| can be very large, thereby violating the inequality relation in (4) that permits use of the A–M relation. This is especially true if the temporal and spatial decay rates of such waves are rapid. The theory of shock dynamics must therefore be extended to take into account the non-uniformity of the flow following the shock wave. This generalization was carried out by Best [10]. Best reconsidered the propagation of a shock wave down a channel of varying cross section and demonstrated that the motion of the shock wave is governed by an infinite sequence of ordinary differential equations given by 1 A (x) dM (5) =− + f (M) Q 1 dx g (M) A(x) n 2 dQ n ∂ ρa u A (x) =− dx ∂t a + u A(x) n i n ∂ 1 Q n+1−i + i ∂t a + u i=1

123

J. T. Peace, F. K. Lu

∂(ρa) ∂u ∂ n−1 ∂(ρa) ∂u − ∂t ∂t ∂ x ∂ x ∂t 1 1 − Q n+1 + a+u a0 M

+

(6)

for n = 1, 2, 3, . . ., where Q n = ∂tn−1 (∂ p/∂t + ρa∂u/∂t)

f (M) = (γ + 1) (γ − 1) M 2 + 2 ·

(γ − 1) M 2 − (γ + 1) ν M + 2 ·

2ρ0 a03 M 2M 2 + (γ + 1) ν M − 2 ·

(γ + 1) ν + 2 (γ − 1) M 3 −1 + (γ + 1) ν M 2 + 4M

(7)

(8)

and

ν2 =

(γ − 1) M 2 + 2 2γ M 2 − γ + 1 (γ + 1) M 2

(9)

The behavior of both g(M) and f (M) is shown in Figs. 17 and 18 of the Appendix. These equations have been recast from that reported by Best for convenience to include ordinary derivatives of the spatial coordinate x, as opposed to time t. By using the general Leibniz rule combined with the characteristic form of the governing equations, Best obtained explicit formulas for the nth-order time and space partial derivatives of p, u, ρ, and a. These formulas are not repeated here; however, it is noted that the expressions were reduced in terms of the RH shock jump conditions and channel area profile. In order to compute a solution, truncation is required to obtain a mathematically closed set of ordinary differential equations in terms of the Mach number M and variable Q n . As a first-order approximation for n = 0, (6) becomes null and therefore truncating Q 1 in (5) yields the classical CCW A–M relation. Moreover, for higher-order approximations, the term Q n+1 is truncated to yield a closed set of n + 1 coupled ordinary differential equations in the variables M and Q n , such that each additional equation is coupled to its successor. This is the same closure scheme of [10], where details regarding the convergence of such a closure scheme were discussed. For this work, only a second-order approximation is considered and the resulting equation set is used to describe the motion of decaying shock waves and blast waves in non-uniform channels. Therefore, for n = 1, the resulting equations consists of (5) and

123

2 ∂ ρa u dQ 1 =− dx ∂t a + u ∂ ∂u + − (ρa) ∂x ∂t

A (x) ∂ ∂u + (ρa) A(x) ∂t ∂x ∂ 1 Q1 ∂t a + u

(10)

where Q 1 = ∂ p/∂t + ρa∂u/∂t, is a measure of flow nonuniformity behind the shock. Thus, with this second-order approximation, the role of flow non-uniformity on the motion of the shock, or pressure and velocity gradients behind the shock, is taken into consideration. In essence, this secondorder approximation allows for an analytical description of waves in non-uniform channels where the initial value of Q 1 is non-zero. In the event that the initial condition Q 1,0 = 0, namely, the gradients in the flow field following the shock are exactly zero, the shock motion is simply that of a uniformly propagating shock as shown in Fig. 3a. However, if Q 1,0 < 0, the shock motion is representative of one that is being overtaken by a rarefaction wave resulting in the propagation of a decaying shock as shown in Fig. 3b. It is also worth noting that for decaying shock and blast waves, the rarefaction following the shock is continuously overtaking the wave causing a decrease in the strength and velocity of the wave as it propagates. Similarly, when the shock is subjected to a nonuniform area change, the effects of the changing area alters the wave dynamics. This process is a coupling between the effects of the following rarefaction and the change in channel cross-sectional area. The current work does not exactly treat this coupling in closed form. However, the effects of the overtaking rarefaction enters the analysis through the initial condition imposed on Q 1 . This aspect will be further discussed in the following sections. As previously mentioned, it is desirable to use the secondorder approximation of (10) in the context of the strong shock limit to analyze the behavior of the S–vN–T blast wave in non-uniform channels. In the development of the S–vN–T blast wave theory, the pressure ahead of the wave is assumed to be negligible compared to the pressure behind the wave. As a result, the Mach number becomes a meaningless quantity and focus must be redirected to the shock wave velocity W = a0 M. Applying this consideration and the strong shock limit to (10), the equations can be recast into the form 1 A (x) dW (11) =− + s (W ) Q 1 dx r (W ) A(x) 2 dQ 1 ∂ ρa u A (x) ∂ ∂u =− + (ρa) dx ∂t a + u A(x) ∂t ∂x 1 ∂ ∂u ∂ (12) Q1 − + (ρa) ∂x ∂t ∂t a + u where r (W ) =

2+λ

λ + 2 (γ − 1) λ2 W

(13)

On the propagation of decaying planar shock and blast waves through non-uniform channels

and s(W ) =

(γ − 1) (γ + 1) (λ − γ + 1) 2ρ0 W 2 (λ + 2) λ + 2 (γ − 1)

validity of the above equations, it is convenient to define the quantity (14)

with λ2 = 2γ (γ − 1)

A (x) dM dM 2 Q 1 + α2 (M) Q 1 + α3 (M) A(x) dx dx 2 dM A (x) + α5 (M) + α4 (M) (18) dx A(x)

Λ = α1 (M)

(15)

Hence, in the strong shock limit for blast waves, two coupled ordinary differential equations are obtained in the shock velocity W and the non-uniformity variable Q 1 . Before proceeding, it is desired to investigate the nature of (10) and (12) and provide a formal criterion regarding their validity. Substitution of the first-order partial time and space derivatives is required to reveal the complete structure of these equations. By using the formulas obtained by Best [10] together with the RH shock jump conditions for (10) and the strong shock limit jump conditions for (12), the following equations are obtained: dM 2 dQ 1 A (x) dM = α1 (M) Q 1 + α2 (M) Q 1 + α3 (M) dx A(x) dx dx 2 dM A (x) + α5 (M) + α4 (M) (16) dx A(x) dQ 1 A (x) dW dW 2 = β1 (W ) Q 1 + β2 (W ) Q 1 + β3 (W ) dx A(x) dx dx 2 dW A (x) + β5 (W ) + β4 (W ) (17) dx A(x)

In this form, αi has the functional dependence αi = αi ( p, u, ρ, a, d p/dM, du/dM, dρ/dM) on the shock, while βi in the strong shock limit has the functional dependence

Thus, by truncation, the accuracy of the two-equation approximation is dependent upon the condition that |Λ|

a + u ∂ ∂ p ∂u 1 1 − + ρa a+u a0 M ∂t ∂t ∂t

(19)

Note that an analogous inequality relation exists for the strong shock limit relation (17). It is now clear that the accuracy of (16) and (17) necessitates that the combined effects of area change and first-order temporal derivatives of flow properties behind the shock dominate the dynamics of the shock compared to the effects of second-order temporal derivatives of flow properties. Further, the righthand side of (19) is scaled by the measure of coincidence between the C+ characteristic and shock wave in x–t space. The behavior of this term was previously discussed in regard to the A–M inequality relation (4). The second term, |∂/∂t(∂ p/∂t + ρa∂u/∂t)|, represents the first-order temporal derivative of flow non-uniformity behind the shock, which can be rewritten as |∂ Q 1 /∂t|. The physical interpretation of this quantity is rendered difficult given the dependence on second-order temporal derivatives of pressure and gas velocity behind the shock. Despite this complexity, the criterion remains true for the validity of a two-equation approximation for shock dynamics governed by an area change and flow non-uniformity following the shock. Further investigation regarding the behavior and estimates for the orders of magnitude of both the remaining terms and the truncated terms is provided in the following sections.

βi = βi ( p, u, a, d p/dW )

3 Shock and blast wave propagation

on the shock. Explicit formulas were obtained within a Mathematica environment for αi and βi . For brevity, these relations are not reported here. However, the behavior of the αi coefficient functions is shown in Figs. 19, 20, 21, 22 and 23 of the Appendix. It should be noted that these figures do not include real gas effects such as vibrational excitation, dissociation, or ionization; therefore, the results at higher Mach numbers in the forgoing analysis should be taken with caution. At this point, it is evident that the shock wave motion for a second-order approximation is governed by two coupled first-order, nonlinear ordinary differential equations. For the purposes of establishing a formal criterion to investigate the

3.1 Decaying planar shock wave For application of (5) and (16) to planar decaying shock waves, it is required to know the behavior of such waves in constant cross-sectional area channels. Sharma et al. [12] provided an analytical solution for the case of a shock wave that is being overtaken by a rarefaction wave. The solution exactly satisfies the governing equations and associated boundary conditions, but only approximately satisfies the entropy and particle velocity jump boundary conditions. Despite this aspect, the resulting error was determined to be very small even for strong shock waves. The nature of the solution is

123

J. T. Peace, F. K. Lu

1011

3.0

1010

2.5

109

2.0 108

1.5

1.0

107 106

0

20

40

60

80

100

Fig. 4 Planar shock wave Mach number decay versus x/x0

self-similar with respect to the variable z = t/h, where h is a Lagrangian mass coordinate and has the relationship, ∂ x/∂h = 1/ρ. Moreover, the solution is obtained by considering a piston in a confined channel that begins to move forward with a designated velocity into an undisturbed gas that is initially at rest. The instantaneous piston motion at t = −tp from a distance x = −xp creates a uniform shock wave with Mach number M that begins to traverse the length of the channel. At a later instant, the piston suddenly stops at t = 0 and distance x = 0, propagating an unsteady rarefaction wave into the disturbed uniform gas behind the shock wave. At a time defined by t0 and distance x0 , the leading characteristic of the unsteady rarefaction catches the shock wave and begins to overtake the shock. Sharma et al. [12] detailed the behavior of the decay from the initial shock strength to an eventual acoustic wave using a derived shock propagation law. This solution is shown for an initial shock Mach number of 3 for various values of γ in Fig. 4. Although not shown in Fig. 4, the far-field behavior as x/x0 → ∞ approaches that of an acoustic wave where M = 1 and ps / p0 = 1 as a result of the Riemann invariant. Despite not showing the far-field solution, note that at x/x0 = 10, the quantity |∇ p/∇ p0 | evaluated on the shock has approached 0.004, 0.028, and 0.052 for γ = 5/3, 7/5, and 4/3, respectively. This is an indication that the spatial decay rate is very rapid in the immediate near-field region of the decaying shock solution. Although the spatial decay rate is shown to be rapid only in the near field, the defined quantity Q 1 is the main parameter of interest for the current work. Evaluation of Q 1 = ∂ p/∂t + ρa∂u/∂t on the shock is shown in Fig. 5. It is evident that despite having a more gradual decay beyond x/x0 = 100, the magnitude of |Q 1 | is relatively large even at x/x0 = 500 for all values of γ . This result also demonstrates the potential

123

0

100

200

300

400

500

Fig. 5 |Q 1 | = |∂ p/∂t + ρa∂u/∂t| evaluated on the decaying shock versus x/x0

invalidity of the A–M relation for describing the dynamics of a planar decaying shock wave that is subjected to a given cross-sectional area change. The magnitude of |Q 1 | in this particular case would likely violate the inequality condition (4) that permits the use of the A–M relation. This will be elaborated on in the next section.

3.2 Sedov–von Neumann–Taylor planar blast wave Similar to the previous section, it is required to know the dynamics of the S–vN–T planar blast wave in a constant cross-sectional area channel for application of (11) and (17). The classical planar blast wave solution was achieved by similarity arguments upon which a self-similar shock motion was derived as the result of an intense explosion subject to the strong shock limit. The solution is readily obtained by introducing an initial charge energy E 0 at the origin of the channel and applying the derived similarity relations. In the far-field where the pressure ahead of the shock is no longer negligible, the strong shock limiting assumption begins to break down, making the blast solution invalid. More important, the similarity begins to break down once p0 is introduced into the defining variables of the physical problem. It is for this reason that the solution of [12] was used in this study to analyze both arbitrary strength decaying shocks and shocks in the strong shock limit. For t > 0, the subsequent blast motion was obtained in closed form by Sedov [13] for a constant initial gas density. However, for this work, the solution procedure outlined by Kamm [17] was implemented as a numerical routine as this method provides a simple integration scheme for evaluating the nondimensional energy constant of the similarity solution. For an initial charge energy of 1 MJ and gas density of ρ0 = 1 kg/m3 , the blast motion and strength are shown in Fig. 6 for various

On the propagation of decaying planar shock and blast waves through non-uniform channels

8

1013

7

1012

6

1011

5

1010

4

109

3

108

2

107

1 0

106

0

20

40

60

80

100

Fig. 6 Planar S–vN–T blast wave velocity decay versus x/x0

values of γ , respectively. In this representation, x0 is defined as the distance travelled by the blast following a time of t0 = 1 ms in order to overcome the singularity at the origin. Similar to the decaying shock wave, the S–vN–T blast wave has a very rapid decay behavior in the near field. In the case of planar blast waves at x/x0 = 10, the quantity |∇ ps /∇ ps,0 | evaluated on the shock has approached 0.01 for all values of γ . This is again an indication that the spatial decay rate is very rapid in the near-field region of the blast wave solution. As mentioned previously, the defined quantity Q 1 is the main parameter of interest for the current work. Therefore, Q 1 is once again evaluated on the shock for the S–vN–T blast solution and the results are shown in Fig. 7. It is evident that despite having a more gradual decay beyond x/x0 = 50, the magnitude of |Q 1 | is very high even at x/x0 = 500 for all values of γ . This result also demonstrates the potential invalidity of the A–M relation for describing the dynamics of a planar blast wave when traversing a cross-sectional area change. The next section details an order-of-magnitude analysis to determine under which conditions the A–M relation and two-equation approximation are valid.

4 Eﬀects of truncation It is desired to examine the conditions for which the criteria in (4) and (19) remain true. This requires assessing the relative order of magnitude for each term contained in these criteria, namely the channel area properties A (x)/A(x), and general flow non-uniformity properties Q 1 and Q 2 . For the channel sectional area properties, a generalized relationship exists given a constant half-angle and starting channel height. This geometrical relation is given by

0

100

200

300

400

500

Fig. 7 |Q 1 | = |∂ p/∂t + ρa∂u/∂t| evaluated on the decaying blast versus x/x0

2n tan θ A (x) = A(x) h 0 + 2x tan θ

(20)

where θ and h 0 correspond to the half-angle of the area change and starting channel height, respectively. In this expression, n represents a geometric index equaling unity for a two-dimensional finite width channel, and 2 for a conical area change in a tube. In this form, a converging geometry corresponds to a half-angle on the range − π/2 < θ < 0 radians, which is constrained over the domain 0 ≤ x ≤ −h 0 /2 tan θ , namely, from the origin of the non-uniform area change to the apex of the converging geometry. Further, a diverging geometry corresponds to a half-angle on the range 0 < θ < π/2, with the domain 0 ≤ x < ∞. For a converging geometry, it is evident from (20) that the minimum of |A (x)/A(x)| occurs at the origin of the area change when x = 0, and equal to 2n tan θ/h 0 . Further, in the limit x → −h 0 /2 tan θ , the value of |A (x)/A(x)| → ∞. Conversely, for a diverging geometry, the maximum of |A (x)/A(x)| occurs at the origin of the area change when x = 0, and is also equal to 2n tan θ/h 0 . Moreover, as x → ∞ the value of |A (x)/A(x)| → 0. A few comments can be made about the nature of the converging and diverging geometries discussed above. It is obvious that in the case of converging geometries, the magnitude of |A (x)/A(x)| grows as the shock traverses through a constricting area profile. As a result, the inequality relation in (4) is more likely to be satisfied in shock dynamic problems where such geometries are encountered. This makes the approximate A–M relation especially useful for the widely studied shock implosion problem as |A (x)/A(x)| grows unbounded as the shock moves toward the focal point at the center of the implosion [18–20]. This will remain true for decaying shock and blast motion with an initially finite

123

J. T. Peace, F. K. Lu

magnitude of flow non-uniformity following the shock. This aspect will be discussed further in the following section. In diverging geometries, the magnitude of |A (x)/A(x)| characterizes the opposite behavior, such that the area function tends toward zero with increasing distance from the origin of the area change. Thus, caution must be taken in such geometries to ensure the A–M relation satisfies (4), as relatively minor flow non-uniformity following the shock could be enough to render the approximation invalid. For this reason, it is useful to approximate the relative order of magnitude of Q 1 and Q 2 for the sake of establishing a range of validity for the approximate shock dynamic expressions in (1) and (10). As previously stated, the measures of flow non-uniformity Q 1 and Q 2 are dependent on temporal partial derivatives of pressure, gas velocity, density, and sound speed behind the shock. Therefore, as a first-order approximation, the relaxation in pressure behind the shock at any point in the flow field can be characterized with the simple Friedlander waveform [21] t e−t/τ (21) p(xi , t) ≈ ps 1 − τ

100

109 1010

10-3

1011

10-4

1012

10-5 10-6

1013

1

2

3

4

5

6

7

8

9

10

Fig. 8 Contour plot of |Q 1 | (kg/m s3 ) versus shock Mach number and time constant of flow non-uniformity τ for air at standard temperature and pressure

100

Further, the density profile, although not entirely similar to the pressure profile, can be estimated using an exponential decay

10-5

1026 10

107 108 109

-2

1010

10

1012 1013 1014

-3

10

1015

10-4

10-6

1016 1017 1018 1019

1

2

3

4

5

6

7

8

9

10

(23)

In these approximate expressions, the subscript s represents the post-shock state and given by the RH shock jump conditions. Note, the above expressions merely approximate the behavior of flow properties following a decaying shock or blast wave, and are thus used only for the sake of argument at estimating the relative order of magnitude for both Q 1 and Q 2 . Assuming an ideal gas, a 2 = γ p/ρ, evaluation of the appropriate temporal derivatives of flow properties at t = 0 yields (24) (25)

Thus, the estimates for flow non-uniformity scales with the properties in the post-shock state and the time constant asso-

123

108

10-2

10-1

1 Q 1 ≈ − (2 ps + ρs as u s ) τ 1 ps u s 1 + ρs as u s Q 2 ≈ 2 3 ps + γ τ as 2

107

10-1

In this form, τ represents the duration of the over-pressure which can be regarded as a time constant for flow nonuniformity behind the shock. In a similar manner, the transient velocity profile at a given point can be estimated with a linear decay expressed as t (22) u(xi , t) ≈ u s 1 − τ

ρ(xi , t) ≈ ρs e−t/τ

106

Fig. 9 Contour plot of |Q 2 | (kg/m s4 ) versus shock Mach number and time constant of flow non-uniformity τ for air at standard temperature and pressure

ciated with the relaxation of flow quantities behind the shock. The more important result regarding these approximations is that Q i ∝ (−1)i /τ i , or equivalently, |Q i | ∝ 1/τ i . In most shock decay and blast wave problems, the time constant of relaxation behind the shock is less than unity. Therefore, it is to be expected that |Q 2 | > |Q 1 | under most conditions. This aspect is what motivates the investigation of the criterion in (19) because the term involving Q 2 is truncated for mathematical closure of the remaining coupled equations. Figures 8 and 9 show contour plots of |Q 1 | and |Q 2 | for ambient air (γ = 7/5) at standard temperature and pressure, respectively.

On the propagation of decaying planar shock and blast waves through non-uniform channels

4.1 One-equation approximation It is now possible to assess the condition under which the classical A–M relation is appropriate for shock motion with a varying degree of non-uniform flow Q 1 in the immediate post-shock flow field. Defining the quantity, Ψ =

1 ρa 2 u

a+u 1− Q1 a0 M

(26)

where ρ, u, and a correspond to the post-shock state and are given by the RH shock jump conditions, it then follows that the criterion for validity of the A–M relation can be reduced to A (x) (27) A(x) |Ψ | Figure 10 is a contour plot of |Ψ | with respect to the shock Mach number and time constant of relaxation for air (γ = 7/5) at standard temperature and pressure. Therefore, given the value of |A (x)/A(x)| at any location in the area change, Fig. 10 can be used to determine if |A (x)/A(x)| |Ψ |. A reasonable approximation would be the condition such that |Ψ |/|A (x)/A(x)| ≤ 0.1. Effectively, this would ensure the A–M relation is used under conditions where the shock motion is primarily influenced by an area change as opposed to the interaction with flow non-uniformity behind the shock. As an example, consider the case where |A (x)/A(x)| = 1 and a shock Mach number of 5. Using the reasoning above, the minimum time constant of decay τ for appropriate use of the A–M relation under these conditions would be restricted to approximately 100 ms. Note that the time constant of decay

4.2 Two-equation approximation In a similar manner, the two-equation approximation can be assessed to determine under which conditions (19) remains true. For this purpose, it is convenient to make use of (18) to define the ratio a + u Q 2 1 (28) 1− Ω= a+u a0 M Λ in which case the criterion for validity of the two-equation approximation reduces to

100

10-2

10-1

10

10-2

10

10-3

|Ω| 1

0

102 103

10-5 104

1

2

3

4

5

6

(29)

-1

101

10-4

10-6

is application and problem specific. However, in most laboratory conditions, a time constant on the order of 100 ms would require establishing decaying shock and blast profiles in linear channels or tubes of impractically long length. An equivalent argument for the specified conditions in the context of flow non-uniformity can also be made. Hence, for |A (x)/A(x)| = 1 and a shock Mach number of 5, the maximum magnitude of flow non-uniformity |Q 1 | behind the shock for appropriate use of the A–M relation would be restricted to approximately O(107 ) kg/m s3 . With respect to both the decaying shock solution of Sharma and the S–vN–T blast wave solution, this magnitude of flow non-uniformity only exists in the far-field region of the solution where the gradients of flow properties are significantly more relaxed. Therefore, under most laboratory conditions, use of the A–M relation for decaying shock and blast motion in non-uniform area channels would be inappropriate making higher-order approximations necessary. The exception would be the converging shock or implosion problem, whereby the magnitude of |A (x)/A(x)| is significantly larger than |Ψ | as the shock converges toward the focus of a contoured geometry or center of the implosion.

7

8

9

10

Fig. 10 Contour plot of |Ψ | (1/m) versus shock Mach number and time constant of flow non-uniformity τ for air at standard temperature and pressure

Note that the order of magnitude of Λ is dependent on dM/dx, area profile properties A (x)/A(x), coefficient functions αi (M), and flow non-uniformity Q 1 . The order of magnitude for dM/dx is estimated by using (5) and the expressions for g(M), f (M), and the approximation for Q 1 in (24). Likewise, the order of magnitude for Λ is evaluated by using the expressions for αi (M) and Q 1 at a specified A (x)/A(x). Figure 11 is a contour plot of |Ω| with respect to shock Mach number and time constant of relaxation in air (γ = 7/5) at standard temperature and pressure for A (x)/A(x) = 1. A reasonable use of the two-equation approximation is under conditions for which |Ω| ≤ 0.1, namely, when the truncated term involving Q 2 negligibly contributes to dQ 1 /dx in (16). Using this reasoning, and the Mach 5 example for assessing the A–M relation, the time

123

J. T. Peace, F. K. Lu

100

5 Results and discussion 10

10-1

10-7 10-6

10-2

10-5 10-4

10-3

10

10-4 10-5 10-6

-8

0.5

1

2

0.4

3

0.3

4

0.2

0.1

5

-3

10-2 0.05

6

7

8

9

10

Fig. 11 Contour plot of |Ω| versus shock Mach number and time constant of flow non-uniformity τ for air at standard temperature and pressure and A (x)/A(x) = 1

constant of decay can be extended from 100 ms to approximately 0.1 ms with the two-equation approximation. Hence, for |A (x)/A(x)| = 1 and a shock Mach number of 5, the maximum magnitude of flow non-uniformity |Q 1 | behind the shock for appropriate use of the two-equation approximation would be extended to approximately O(1010 ) kg/m s3 . This is an important result that demonstrates the degree to which the two-equation approximation can account for non-uniform flow following the shock and still accurately represent the dynamics of a shock in a non-uniform area change. In the context of the incident Mach 3 decaying shock solution of Sharma shown in Fig. 5, this order of magnitude of flow non-uniformity takes place within the proximity of the near-field solution. Alternatively, for the planar, 1-MJ S–vN–T blast wave solution in Fig. 7, this order of magnitude in flow non-uniformity occurs in the mid-field region. Therefore, it can be concluded that there exist practical problems regarding decaying shock and blast motion for which the two-equation approximation is appropriate. Another important feature of the two-equation approximation shown in Fig. 11 is the behavior in the weak shock limit. From the contour lines of |Ω|, it is clear that in the limit of weak shocks, the equations can be applied to problems with significantly smaller values of τ without exceeding |Ω| ≤ 0.1. This result is to be expected as the measure in coincidence between the C+ characteristic and shock approaches zero, which in turn cancels the effects of flow non-uniformity behind the shock. Therefore, the twoequation approximation can properly be used in the weak shock limit to describe the dynamics of decaying shocks in non-uniform area changes.

123

It is desired to analyze and discuss the nature of the solutions to the coupled equations (5) and (16) and (11) and (17), for the case of converging and diverging channels for both arbitrary strength shocks and shocks in the strong shock limit, respectively. For the remainder of this study, the channel (n = 1) under consideration consists of a channel featuring a 1 cm starting channel height with either a converging or diverging profile characterized by a half-angle of θ = 1◦ . For the initially uniform planar shock motion, it follows that Q 1,0 ≡ 0 and a direct comparison can be made between the A–M relation and the second-order approximation. In essence, making Q 1,0 = 0 allows for isolating the effects of Q 1 on the solution behavior. Figure 12 shows the behavior of planar shock waves with initial Mach numbers of 2 and 4 for both converging and diverging channels. It is clear that for converging channels, the A–M relation and the secondorder approximation have very good agreement. This is to be expected as the converging channel area is the primary source that influences the strength of the shock. Additionally, in the case of Q 1,0 = 0, the contribution of |Q 1 | in the converging geometry becomes negligible when scaled by the difference in coincidence of the C+ characteristic and shock, and the 1/ρau 2 term in (26). However, for the diverging channel, the two solutions begin to diverge. This is caused by the initial increase in |Q 1 | as the shock traverses the diverging channel. Figure 13 better demonstrates this deviation for diverging channels over the range 1 ≤ A(x)/A0 ≤ 103 . It is evident that the second-order approximation approaches an acoustic wave much faster for both initial shock Mach numbers. In the weak shock limit, the two solutions will be identical as the coincidence between the C+ characteristic and shock approaches zero, rendering the effects of flow non-uniformity irrelevant. However, as the shock strength is increased, the deviation is more pronounced as indicated by the solution for M0 = 4 in Fig. 13. As expected, this deviation is attributed to the introduction of Q 1 into the twoequation approximation. Although the flow non-uniformity initial condition Q 1,0 = 0 is being imposed, once the wave begins to propagate into the diverging channel, the value of Q 1 becomes non-zero on the shock as dQ 1 /dx is dependent on both A (x)/A(x) and dM/dx. Note that the effects of re-reflected disturbances are not accounted for in this study. Moreover, the degree to which these disturbances influence the shock trajectory as opposed to the changing area and |Q 1 | may become more significant at large distances from the start of the area change. Further analysis, similar to that carried out by Milton [22], would be required to formally assess this aspect of decaying shock and blast propagation in non-uniform area changes.

On the propagation of decaying planar shock and blast waves through non-uniform channels

25

9 A-M 2-Eq.

8

20

7 6

Diverging

Converging

Diverging

Converging

A-M 2-Eq.

15

5 4

10

3 2

5

1 0 10-2

10-1

100

101

Fig. 12 Comparison of A–M relation and second-order approximation for initial shock Mach numbers of 2 and 4 in converging and diverging channels

10-1

100

101

102

Fig. 14 Comparison of A–M relation and second-order approximation for initial shock velocities of 5 and 10 km/s in the strong shock limit for converging and diverging channels

9

3.0 A-M 2-Eq.

2.5

A-M

8 7 6

2.0

5

1.5

4

Diverging

3

1.0

2

0.5 0.0 100

0 10-2

102

1 101

102

103

Fig. 13 Comparison of A–M relation and second-order approximation for initial shock Mach numbers of 2 and 4 in diverging channels

Figure 14 compares the A–M relation and the secondorder approximation in the strong shock limit for Q 1,0 = 0 with initial shock wave velocities of 5 and 10 km/s. This figure clearly shows that for converging channels, the A– M relation and second-order approximation have very good agreement. Once again, this is an indication that the influence of Q 1 in converging channels is negligible compared to the effects of the converging channel area. The same is not true for the case of diverging channel areas. In the strong shock limit, the second-order approximation yields a non-physical zero shock velocity at A/A0 ≈ 37.5. Similar results were obtained by Best [10] for spherical and cylindrical diverging geometries in the strong shock limit. This is caused by a divergence in Q 1 for diverging channel area geometries.

Converging

0 10-2

10-1

100

101

102

103

Fig. 15 Comparison of A–M relation and second-order approximation for an initial shock Mach number of 4 for converging and diverging channels at various Q 1,0

For weaker shocks, this is avoided as the strong shock limit assumption is removed and the C+ characteristic and shock trajectory become more coincident. It is also desirable to consider the influence of Q 1,0 on the behavior of the two-equation approximation in both converging and diverging channel geometries. For an initial shock wave Mach number of 4, Fig. 15 shows a comparison of the A–M relation and solutions to (5) and (16) for various values of flow non-uniformity Q 1,0 = {0, −106 , −107 , −108 , −109 } kg/m s3 . For the case of a converging channel, it is evident that the planar shock dynamics become independent of the initial non-uniformity of the flow Q 1,0 when |Ω| ≤ 0.1. As

123

J. T. Peace, F. K. Lu

123

3.0

100

100

10-1

100

100

101

2.5 Converging

Ms

a result, all the solutions of the second-order approximation collapse on the A–M relation. A similar result was obtained by Best [10] for the case of converging cylindrical and spherical shocks in the strong shock limit. Nonetheless, this result has yet to be reported for an arbitrary strength shock in a converging planar geometry. It should be noted that a deviation in the converging channel solution was achieved for values of Q 1,0 that were orders of magnitude above those obtained from the solution of Sharma et al. [12] for a decaying shock wave with an initial Mach number of 3. Moreover, these conditions with such large Q 1,0 corresponded to situations where |Ω| is no longer small, thereby making the inequality relation that permits use of the two-equation approximation invalid. Hence, given the physical considerations of the decay process and limitations of the two-equation approximation, the values were deemed non-physical and disregarded from the analysis. For the diverging channel geometry, the initial value of Q 1,0 significantly influences the rate at which the decaying shock wave approaches an acoustic wave. For these cases, it intuitively follows that the larger the initial value of |Q 1,0 |, the shock wave decays faster as the temporal gradient of flow properties behind the shock are larger. However, the solutions reported in Fig. 15 must be taken with caution as the influence of re-reflected disturbances have not been taken into account. It is likely that the solutions are more accurate for smaller channel expansion area ratios; however, this would require further analysis and possibly experimental data to confirm. For strong shock solutions to (11) and (17), the same behavior was obtained for the arbitrary decaying shocks in converging geometries. The second-order approximation solutions are independent of the initial non-uniformity of the flow Q 1,0 and all the subsequent solutions collapsed on that obtained with the A–M relation. As a result, for an initial shock velocity of 10 km/s, the solutions are exact to the converging channel geometry shown in Fig. 14. Deviations in the solutions for converging geometries where obtained for very large initial values of Q 1,0 ; however, these values were orders of magnitude above those obtained from the S–vN–T blast solution and were once again deemed non-physical. Lastly, the influence of Q 1,0 for the case of diverging channel area was not analyzed due to the non-physical behavior of the solution in the strong shock limit. A final shock propagation case is considered to demonstrate the capability of the second-order approximation in application to a decaying shock wave of arbitrary strength in converging and diverging channels. The solution of Sharma

2.0

1.5

1.0

Sharma A-M 2-Eq. 0

20

Diverging

40

60

80

100

Fig. 16 Comparison of A–M relation and second-order approximation for a decaying shock wave in a non-uniform channel

et al. [12] is used where the description of the wave process has previously been provided. In short, a piston is set in motion and moments later is suddenly stopped such that a uniform planar shock wave with a Mach number of 3 is created. At (t0 , x0 ), the leading characteristic of the rarefaction wave catches and begins to overtake the traveling shock. At x/x0 = 50, the planar shock wave has decayed to a Mach number of M0 = 1.96 with non-uniformity value of Q 1,0 = −1.696 × 108 kg/m s3 . These two parameters are the initial conditions for the second-order approximation of (5) and (16). Similarly, at a location of 50 ≤ x/x0 ≤ 100, the channel features either a converging or diverging section with total area ratios of 0.1 and 10, respectively. The solution to this particular geometry is shown in Fig. 16. It is noted that the solid line indicates the solution of Sharma et al. and is assumed to continuously decay in a uniform channel for x/x0 > 50. As expected, the converging channel case in Fig. 16 is independent of the flow non-uniformity behind the decaying shock wave. Alternatively, for the diverging channel case, a more relaxed solution is obtained from that of the A–M relation as a result of the initial flow non-uniformity condition Q 1,0 . The secondorder solution is believed to be more accurate as a result of taking into account the flow non-uniformity Q 1 behind the shock. However, experimental verification must be made to investigate the complete accuracy of the predicted solution.

On the propagation of decaying planar shock and blast waves through non-uniform channels

6 Conclusions The current work has demonstrated two possible shock wave cases where the classical A–M relation is potentially invalid when the flow following the shock is non-uniform. These cases include the propagation of an arbitrary strength decaying planar shock wave and the propagation of planar blast waves in channels of varying cross-sectional area. In regard to this aspect, a second-order approximate equation set was used to take into account the non-uniformity of flow following the shock. A generalized order-of-magnitude analysis is carried out to analyze under which conditions the classical A–M relation and two-equation approximation are valid given a time constant of decay for the flow properties behind the shock. It was shown that the two-equation approximation is applicable in cases where the time constant of decay is orders of magnitude lower than that for appropriate use of the A–M relation. However, in the weak shock limit, both solutions become identical as the difference in coincidence between the C+ characteristic and shock path in x–t space approaches zero, effectively canceling the effects of flow non-uniformity in the post-shock flow field. It was also shown that the A–M relation and second-order approximation have very good agreement for the case of converging channels. More important, it was determined that the secondorder approximate solution becomes independent of flow non-uniformity following the shock, such that the dynamics of the shock are solely governed by the changing area in the converging channel when |Ω| ≤ 0.1. This is explained by the growth in |A (x)/A(x)| for a converging geometry, which dominates the contributing effects of flow non-uniformity on the shock motion. In the case of diverging channels, the second-order approximation was shown to strongly depend on the magnitude of flow non-uniformity behind the shock. Even for the case of a uniformly propagating shock with Q 1,0 = 0, the second-order solution slightly deviates from the A–M relation. It was shown that the deviation becomes more pronounced for higher values of |Q 1,0 |, such that the two-equation approximation tends toward an acoustic wave faster than the A–M relation. The second-order approximation is considered to be more accurate in diverging channels within a limiting channel length. However, a formal consideration of the re-reflected disturbances behind the wave would be required to assess the overall accuracy of the solution.

100

10

1

0.1

2

3

4

5

6

7

8

9

10

7

8

9

10

8

9

10

Fig. 17 g(Ms ) versus shock Mach number

10-1

10-2

10-3

10-4

1

2

3

4

5

6

Fig. 18 a0 p0 | f (Ms )| versus shock Mach number

1.0

0.9

0.8

0.7

0.6

Appendix

1

1

2

3

4

5

6

7

Fig. 19 |α1 (Ms )| versus shock Mach number

See Figs. 17, 18, 19, 20, 21, 22 and 23.

123

J. T. Peace, F. K. Lu

104

102

103 1

10

102 100

101

10-1

1

2

3

4

5

6

7

8

9

10

Fig. 20 |α2 (Ms )| versus shock Mach number

100

1

2

3

4

5

6

7

8

9

10

Fig. 23 |α5 (Ms )|/a0 p0 versus shock Mach number

103

References

102 101 100 10-1 10-2

1

2

3

4

5

6

7

8

9

10

8

9

10

Fig. 21 |α3 (Ms )|/a0 p0 versus shock Mach number

104

103

102

101

1

2

3

4

5

6

7

Fig. 22 |α4 (Ms )|/a0 p0 versus shock Mach number

123

1. Chester, W.: The quasi-cylindrical shock tube. Philos. Mag. 45(371), 1293–1301 (1954). https://doi.org/10.1080/ 14786441208561138 2. Chisnell, R.F.: The motion of a shock wave in a channel, with applications to cylindrical and spherical shock waves. J. Fluid Mech. 2(3), 286–298 (1957). https://doi.org/10.1017/ S0022112057000130 3. Whitham, G.B.: On the propagation of shock waves through regions of non-uniform area or flow. J. Fluid Mech. 4(4), 337–360 (1958). https://doi.org/10.1017/S0022112058000495 4. Russell, D.A.: Shock-wave strengthening by area convergence. J. Fluid Mech. 27(2), 305–314 (1967). https://doi.org/10.1017/ S0022112067000333 5. Setchell, R.E., Storm, E., Sturtevant, B.: An investigation of shock strengthening in a conical convergent channel. J. Fluid Mech. 56(3), 505–522 (1972). https://doi.org/10.1017/S0022112072002484 6. Nettleton, M.A.: Shock attenuation in a ‘gradual’ area expansion. J. Fluid Mech. 60(2), 209–5223 (1973). https://doi.org/10.1017/ S0022112073000121 7. Igra, O., Elperin, T., Falcovitz, J., Zmiri, B.: Shock wave interaction with area changes in ducts. Shock Waves 3(3), 233–238 (1994). https://doi.org/10.1007/BF01414717 8. Whitham, G.B.: A new approach to problems of shock dynamics. Part I. Two dimensional problems. J. Fluid Mech. 2(2), 145–171 (1957). https://doi.org/10.1017/S002211205700004X 9. Ridoux, J., Lardjane, N., Monasse, L., Coulouvrat, F.: Comparison of geometrical shock dynamics and kinematic models for shockwave propagation. Shock Waves 28(2), 401–416 (2017). https:// doi.org/10.1007/s00193-017-0748-2 10. Best, J.P.: A generalization of the theory of geometrical shock dynamics. Shock Waves 1(4), 251–273 (1991). https://doi.org/10. 1007/BF01418882 11. Kirkwood J.G., Bethe H.A.: The pressure wave produced by an underwater explosion. Office of Scientific Research and Development Report No. 588 (1942) 12. Sharma, V.D., Ram, R., Sachdev, P.L.: Uniformly valid analytical solution to the problem of a decaying shock wave. J. Fluid Mech. 185, 153–170 (1987). https://doi.org/10.1017/ S0022112087003124

On the propagation of decaying planar shock and blast waves through non-uniform channels 13. Sedov, L.I.: Similarity and Dimensional Methods in Mechanics, Chapter IV, “One-Dimensional Unsteady Motion of a Gas”. Academic Press, New York (1959). https://doi.org/10.1016/B978-14832-0088-0.50011-6 14. Bethe, H.A., Fuchs, K., Hirschfelder, J.O., Magee, J.L., Peierls, R.E., von Neumann, J.: Blast waves. Los Alamos Scientific Laboratory Report LA-2000 (1947) 15. Taylor, G.I.: The formation of a blast wave by a very intense explosion. Proc. R. Soc. Lond. A 201, 159–174 (1950). https://doi.org/ 10.1098/rspa.1950.0049 16. Han, Z., Yin, X.: Shock Dynamics, Chapter 1, “Relation Between M and A for a Uniform Quiescent Gas Ahead of a Shock Wave”. Kluwer, Beijing (1993). https://doi.org/10.1007/978-94017-2995-6_2 17. Kamm J.R.: Evaluation of the Sedov–von Newman–Taylor blast wave solution. Los Alamos Scientific Laboratory Report LA-UR00-6055 (2000)

18. Takayama, K., Kleine, H., Grönig, H.: An experimental investigation of the stability of converging cylindrical shock waves in air. Exp. Fluids 5(5), 315–322 (1987). https://doi.org/10.1007/ BF00277710 19. Zha, Z., Si, T., Luo, X., Yang, J., Liu, C., Tan, D., Zou, L.: Parametric study of cylindrical converging shock waves generated based on shock dynamics theory. Phys. Fluids 24, 026101 (2012). https:// doi.org/10.1063/1.3682376 20. Liverts, M., Apazidis, N.: Limiting temperatures of spherical shock wave implosion. Phys. Rev. 116, 014501 (2016). https://doi.org/10. 1103/PhysRevLett.116.014501 21. Friedlander, F.G.: The diffraction of sound pulses I. Diffraction by a semi-infinite plane. Proc. R. Soc. A 186(1006), 322–343 (1946). https://doi.org/10.1098/rspa.1946.0046 22. Milton, B.E.: Mach reflection using ray-shock theory. AIAA J. 13(11), 1531–1533 (1975). https://doi.org/10.2514/3.60566

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

123