## Abstract

Compressed-gas-driven shock tubes have become popular as a laboratory-scale replacement for field blast tests. The well-known initial structure of the Riemann problem eventually evolves into a shock structure thought to resemble a Friedlander wave, although this remains to be demonstrated theoretically. In this paper, we develop a semi-analytical model to predict the key characteristics of pseudo blast waves forming in a shock tube: location where the wave first forms, peak over-pressure, decay time and impulse. The approach is based on combining the solutions of the two different types of wave interactions that arise in the shock tube after the family of rarefaction waves in the Riemann solution interacts with the closed end of the tube. The results of the analytical model are verified against numerical simulations obtained with a finite volume method. The model furnishes a rational approach to relate shock tube parameters to desired blast wave characteristics, and thus constitutes a useful tool for the design of shock tubes for blast testing.

## 1. Introduction

In addition to their many uses in the science and engineering of compressible flows, shock tubes have recently become popular as devices to generate shock waves whose structure resembles a real blast wave such as produced by the detonation of explosive charges (e.g. [1–10]). Lab-scale tests using shock tubes are convenient surrogates for field blast tests as they provide better precision and variable control, increased safety and reduced cost [8].

The basic functioning of a shock tube in its traditional configuration is well understood (e.g. [11]): pressurized gas in a driver section is suddenly released by the rupture of a membrane, leading to a well-known shock structure comprising a shock wave, a contact discontinuity and a family of rarefaction waves [12]. What is observed experimentally is that, for sufficiently long shock tubes, this wave structure eventually evolves into a shape resembling an air blast wave. What has received less attention is the process by which this happens and the relation between shock tube parameters and the characteristics of the resulting ‘blast wave’. A remaining open question [13] is whether the pseudo blast wave that forms in the shock tube actually corresponds to the commonly accepted quasi-exponential Friedlander wave form [14,15]:
*p*(*t*) as a function of time after arrival of the wave at a given location; Δ*p* is the peak overpressure, *τ* the decay time, *α* the wave form parameter and *p*_{0} the ambient pressure.

In fact, it has been well established that, despite their popularity as blast test surrogates, shock tubes sometimes fail to generate true blast waves [6]. Instead, if the shock tube parameters are not designed ‘correctly’, the shock waves obtained have a tendency to adopt trapezoidal shapes where the peak pressure plateaus before it decays as blast waves do [1].

Knowing the initial location where the blast wave forms is important for properly placing the target in animal studies of blast-induced Traumatic Brain Injury (TBI) [5]. However, estimates of the location of the onset of the blast wave are mostly empirical. In Refs [4,8], Bass and co-workers recognize the uncertainty on the location where the Friedlander wave forms and adopt the ‘rule of thumb’ for shock design that the driven length to diameter ratio should be greater than 10 to ensure that the wave is planar. In Ref. [6], Courtney *et al.* also describe an approach in the design of shock tubes for blast testing (in their case for combustion-driven shock tubes) based on empirical rules rather than prescriptive mathematical expressions. Specifically, they mention previous work suggesting ‘…that choosing the length of the driven section to be about 60 diameters resulted in a desirable blast loading profile’ in combustion-driven shock tubes (rifles), and also mention that ‘the length of most of the driven sections tested were about ten times the length of the driver sections, which is consistent with early published shock tube designs’ [6], p. 0451112. The inconsistency and uncertainty in the empirical determination of proper shock tube parameters therefore makes it clear that there is a need for developing a rational approach for designing shock tubes for laboratory-scale blast testing. It is important to note that our contribution in this paper does not address the case of combustion-driven shock tubes.

Previous efforts to analyse the formation of Friedlander waves in shock tubes have been based on computational fluid dynamics (CFD) analysis. CFD simulations have shown that the Riemann shock wave structure eventually leads to the formation of a wave pressure profile which looks like a Friedlander blast wave. The accuracy of CFD numerical tools has been validated by comparing the time evolution of the pressure at specific location in the tube with experimental pressure sensor data [4]. CFD has also been instrumental in describing the flow conditions as the wave approaches and goes past the open end of the tube [7,16].

A recent review by Needham *et al*. [17] extensively identifies all of the limitations in shock-tube blast testing for TBI and provides some general guidance on how computational approaches can be used to better design experimental tests.

In this paper, we develop an analytical model that provides a functional relationship between shock tube parameters (driver section length, pressure and driver gas) and the sought blast wave characteristics in closed mathematical form. As in many other areas of engineering and science, an analytical description, if available, is highly desirable in the design and configuration of the testing device, as it eliminates the need for costly parametric studies using numerical solution methods. Moreover, analytical models enable the immediate solution of the design (inverse) problem, i.e. how to configure the problem inputs to obtain a desired set of outcomes. In the specific case of laboratory-scale blast testing, this translates to: what is the required set of shock tube configuration parameters (initial pressure, driver section length and gas type) to produce a surrogate blast wave with given characteristics (free-field overpressure and decay time, or alternatively explosive energy and offset), and what is the minimum required tube length?

The analytical model is derived by analysing two types of nonlinear wave interactions that arise as the family of rarefaction waves in the initial Riemann problem reflects from the closed end of the tube to eventually catch up with the shock wave. The first problem is concerned with the reflection of the family of centred rarefaction waves (rarefaction fan) off the closed end of the driver section. The solution of the flow in this so-called buffer region [18] is known to be non-simple and has been derived using Riemann invariants, leading to a second-order hyperbolic partial differential equation (PDE) of the Euler–Poisson type [12]. For arbitrary gas-specific heat ratios, the solution to this PDE can be expressed as a hypergeometric function which gives the position of the interacting waves as a function of time in implicit form. We propose an alternative derivation resulting in an Euler–Darboux PDE whose solution in the important cases of air and helium can be written in terms of much simpler rational functions. This, in turn, is convenient for the next steps in the analysis. The simplified solution is used to find the exit time and position where the reflected head wave intersects the incident tail wave, at which point the head and subsequent reflected waves enter a simple region (straight characteristic lines) with known propagation speed.

The second problem involves the interaction of these reflected waves when they reach the contact discontinuity moving in the direction of the shock. The flow in this second buffer region can be analysed using Courant’s solution for a simple wave interacting with a contact discontinuity [19]. A discrete consideration of equally spaced incident waves combined with a linearized analysis of their propagation velocity furnishes a discrete estimate of the curved geometry of the modified contact interface, as well as the time, position and propagation speed of the simple waves transmitted across it. Finally, this information is used to determine the location and time when the head wave encounters the shock front, which is considered the onset of the blast wave structure. This furnishes a closed form expression relating shock tube parameters (driver’s section length and compressed gas state) to the location of the onset of the blast wave in the tube and its peak intensity. The subsequent evolution of the flow at that location (pressure decay) is estimated by considering the arrival time and pressure of a discrete set of trailing waves, which results in a straightforward semi-analytical procedure to determine the wave decay time and impulse. The model shows that shock dynamics imposes a requirement on driver section length for the blast wave to form. This length is independent of the shock tube diameter, as long as the 1D assumptions hold (e.g. the diameter of the tube has to be much larger than the thickness of the boundary layer, which is generally the case in practice). It bears emphasis that other requirements may apply, e.g. the need for a minimum length to achieve a planar blast wave [6,8] owing to three-dimensional local effects at the ruptured membrane which are affected by the tube diameter and are not considered in the one-dimensional model proposed here.

The structure of the paper is as follows: We discuss the analysis of the different wave interactions and develop analytical estimates of the characteristics of the pressure profile at the onset of the Friedlander wave in §2. We compare our analytical model against a numerical approach that we introduce in §3a. We present numerical and analytical results in §3. In §4, we summarize the main conclusion of the paper.

## 2. Analytical estimation of the characteristics of a Friedlander wave forming in a shock tube

The basic functioning of a shock tube driven by compressed gas is well understood and has been extensively studied (e.g. [11]). A shock tube consists of a driver section containing gas—commonly helium or air—at high pressure, and a driven section, initially separated by a membrane (figure 1*a*). Upon rupture of the membrane, a shock wave propagates towards the driven section followed by a contact discontinuity, whereas a family of rarefaction waves propagates in the opposite direction within the driver section (figure 1*b*). While unperturbed by the end sections of the shock tube, the flow corresponds well to the known Riemann solution [20]. Here, we are interested in the various wave interactions taking place upon the reflection of the head rarefaction wave at the left end of the driver section (see figure 1*a* for a schematic of the configuration). In the subsequent evolution of the flow, the reflected head wave reaching the right-propagating shock marks the onset of the Friedlander wave form in the shock tube. Figure 1 shows snapshots of the shock tube at different times emphasizing the location of key flow characteristics (shock, contact discontinuity, head and tail waves) as well as the evolution of the spatial pressure distribution. To help the analysis, we make heavy use of figure 2 which shows the evolution of the characteristic lines in the *x*–*t* plane.

The basic idea of the analysis is to track the evolution (velocity history) of the head and trailing rarefaction waves in the various stages of their interaction with the flow: (i) propagation towards the left closed end of the tube (figure 1*b*), (ii) interaction with trailing rarefaction waves after reflection at the closed end (figure 1*c*), (iii) propagation towards the moving contact discontinuity (figure 1*d*), (iv) interaction with the contact discontinuity (figure 1*e*) and (v) propagation towards the shock (figure 1*f*).

Next, we focus on the determination of the location in the shock tube where this occurs.

### (a) Determination of the location of the onset of the Friedlander wave form

The determination of the location of the onset of the Friedlander wave form only requires tracking the evolution of the head wave. Immediately after the membrane ruptures (figure 1*b*), the constant speed of the head wave is *a*_{4}, *γ*_{4}, *p*_{4} and *ρ*_{4} are, respectively, the sound speed, the specific heat ratio, the pressure and the density in the (undisturbed region of the) driver section. It is worth emphasizing that the parameters on the right-hand side of this equation constitute basic input settings of the shock tube. The time *t*_{a} at which the head rarefaction wave reaches the closed end of the shock tube, figure 1*c* and point *a* in figure 2, is simply given by:

After time *t*_{a}, the flow in the shock tube departs from the classical Riemann problem and the head rarefaction wave enters a non-simple region: the first buffer region.

#### (i) Analysis of the first buffer region

The first buffer region comprises the zone where the centred fan of rarefaction waves interacts with the closed end of the tube (figures 1*c*, 2 and 3*a*). Riemann [12] originally proposed a mathematical formulation of this problem in the form of an Euler–Poisson PDE governing the time *t*(*r*,*s*) at which two characteristics with coordinates (*r*,*s*) meet inside the buffer region. He obtained a general solution for a gas with arbitrary specific heat ratio *γ*. However, the solution is difficult to use in the analysis of the further evolution of the flow in the shock tube, as it is given in implicit form in terms of a hypergeometric function. Landau & Lifschitz [21] derived an explicit expression of the exit time of the head rarefaction wave from the buffer region for monoatomic (*γ*=1.4) gases, based on a reformulation of the PDE using the Legendre transformation. Here, we propose a third alternative derivation which is limited to the cases of practical interest (air and helium), but has the great advantage of providing a full solution for *t*(*r*,*s*) in explicit form in terms of a simple rational function.

Figure 3*a* shows a plot in the (*x*,*t*) plane of the characteristics of the head wave, the tail wave and two sample intermediate rarefaction waves as they propagate towards the closed end of the driver section, reflect off the wall and interact with each other in the first buffer region. Prior to entering the first buffer region, rarefaction waves move at constant speed along a straight characteristic. Any arbitrary point *A* inside the buffer region lies at the intersection of two rarefaction waves, a right-going wave along a *C*_{+} characteristic and a left-going wave along a *C*_{−} characteristic. The wave speed for these two waves is given respectively by:
*u*=*u*(*x*,*t*) is the local gas speed and *a*=*a*(*x*,*t*) is the local sound speed. On any given characteristic, *C*_{+} or *C*_{−}, the evolution of these quantities is governed by the Riemann invariants [18] as follows:
*u* and the local sound speed *a* can be conveniently expressed in terms of the Riemann invariants *r* and *s* as follows:
*C*_{+} characteristic, the position *x* at a time *t* of a right-going wave only depends on the value of the Riemann invariant *s* of intersecting *C*_{−} characteristics. Conversely, on a *C*_{−} characteristic, the position *x* at a time *t* of a left-going wave depends only on the Riemann invariants *r* of the intersecting *C*_{+} characteristics. These features of the Riemann invariants motivate the change of variables from space–time (*x*,*t*)-coordinates to characteristic space (*r*,*s*)-coordinates.

We aim to determine the time *t*(*r*,*s*) at which an incoming left-going rarefaction wave meets a right-going rarefaction wave reflected from the rigid wall. From (2.2) and using the invariance of *r* and *s* along *C*_{+} and *C*_{−}, respectively, we obtain:
*x* variable can be eliminated by invoking the equality of mixed particles to arrive at:
*u*) and sound (*a*) speeds in terms of the Riemann invariants using (2.4), this expression can be further simplified into a PDE for the time variable *t* only:
*γ*_{4}=1.4) and λ=2 for helium (*f* and *g*, in the following form [22]:

The functions *f* and *g* are determined from the initial conditions of the problem, which correspond to the time at which any rarefaction wave enters the first buffer region. Left-going rarefaction waves, emitted from the centre *P* and travelling on *C*_{−} characteristics, enter the buffer region when they intersect the (reflected) right-going head wave (in red in figure 3*a*). On this characteristic of Riemann invariant *r*_{0}, the derivative of the position *x* and the time *t* with respect to the invariant *s* are related through (2.5). Using (2.4), this expression reads:
*x*=(*u*−*a*)*t*+*x*_{P}, where *x*_{P} is the location of the centre of the rarefaction fan. Differentiating this expression with respect to *s* and combining with (2.4) yields another relation between the derivatives of the position *x* and time *t* with respect to the invariant *s*:
*x* is eliminated by taking the difference between (2.10) and (2.11) and simplifying, resulting in the following ordinary differential equation (ODE):

Right-going reflected rarefaction waves travelling on *C*_{+} characteristics enter the first buffer region at the time they reflect off the wall. However, due to the hyperbolic nature of the problem, these times depend on the solution of the PDE inside the buffer region itself. This issue is easily circumvented by recourse to the method of images (e.g. [23]), by which the problem can be replaced with the interaction of two centred fans as shown in figure 3*b*. In this case, the *C*_{−} characteristic *s*_{0} can be determined by an ODE analogous to (2.12):

Solving both ODEs, we arrive at an explicit expression of the entry time of the left-going rarefaction waves, *t*(*r*_{0},*s*), and the right-going rarefaction waves, *t*(*r*,*s*_{0}), on the red and blue curves that bounds the buffer region in figure 3*b*, respectively:

The initial conditions (2.14) and (2.15) can now be used to determine the functions *f*(*r*) and *g*(*s*) in (2.9). In the case of air (λ=3), the solution (2.9) can be written as:
*f* and *g* with respect to their single argument. Evaluating (2.16) at the two initial conditions *t*(*r*,*s*_{0}), (2.14), and *t*(*r*_{0},*s*), (2.15), yields two independent ODEs for *f* and *g*, respectively:

Here, we are interested in any solutions *f*(*r*) and *g*(*s*) that satisfy the condition *t*(*r*_{0},*s*_{0})=*t*_{a}, i.e. the first buffer region starts at the time that the head rarefaction wave arrives at the fixed wall (2.1). We may choose any initial conditions for the ODEs (2.17) and (2.18) as long as this condition is satisfied. We choose *f*(*r*_{0})=1, *f*′(*r*_{0})=1, *f*′′(*r*_{0})=(*t*_{a}/2)(*r*_{0}−*s*_{0})^{3}, *g*(*s*_{0})=1, *g*′(*s*_{0})=−1 and *g*′′(*s*_{0})=0. This yields an explicit form for the two ODEs:

Assuming the solutions to (2.19) and (2.20) are, respectively, of the form (*r*−*s*_{0})^{m} and (*r*_{0}−*s*)^{m}, the solutions of the corresponding homogeneous problems are given by {(*r*−*s*_{0})^{3},(*r*−*s*_{0})^{4}} and {(*r*_{0}−*s*)^{3},(*r*_{0}−*s*)^{4}}. The particular solutions can be obtained by inspection as: *a*_{1}(*r*−*s*_{0})^{2}+*a*_{2}(*r*−*s*_{0})+*a*_{3} and *b*_{1}(*r*_{0}−*s*)+*b*_{2}, respectively. The general solution can then be obtained by the method of undetermined coefficients inserting these expressions in the ODEs and solving for the coefficients, leading to:

Using a similar approach, the solution of (2.9) for the case of helium (λ=2) can be obtained in the form:

The head rarefaction wave of Riemann invariant *r*_{0} exits the buffer region when it meets the tail rarefaction wave of Riemann invariant *s*_{t} at time *t*_{b}=*t*(*r*_{0},*s*_{t}) (figure 3*b*). Conversely, *t*_{b} can be viewed as the time the tail rarefaction wave enters the first buffer region. Its value is then easily computed from the initial condition of the Euler–Darboux PDE for the left-going rarefaction waves (2.14) in terms of the Riemann invariants *s*_{t}, *r*_{0} and *s*_{0}. In the Riemann solution, the head wave moves to the left into the undisturbed region of the driver section, denoted region 4 in figure 1*a*. In this region, the local gas velocity is zero and the local sound speed is *a*_{4}. The tail wave moves left from the region behind the contact interface, denoted region 3 in figure 1*a*. In this region, the local gas velocity is *u*_{p} and the local sound speed is *a*_{3}. The Riemann invariants *r*_{0}, *s*_{0} and *s*_{t} are obtained using (2.3):
*a*_{3} and *a*_{4} can be obtained by noticing that the Riemann invariant *r* (equation (2.3)) is constant in the expansion wave. This constant must match the value at the quiescent region, where *u*=0,*a*=*a*_{4}, and evaluates to 2*a*_{4}/(*γ*_{4}−1). The result is *a*_{3}=*a*_{4}−((*γ*_{4}−1)/2)*u*_{p}. Using this relation in equation (2.25), we obtain:
*u* and local sound speed in an expansion wave *a*/*a*_{4}=1−((*γ*_{4}−1)/2)(*u*/*a*_{4}). The denominator can be simplified using the isentropic equation that relates the constant pressure behind the contact interface (*p*_{3}) to the initial pressure in the driver section(*p*_{4}):
*p*_{3}=*p*_{2}, (2.26) reduces to:
*t*_{b} for the head wave (equation (2.27)). However, the full solution is required in our subsequent analysis of trailing waves for the determination of the characteristics of the blast wave (e.g. decay time).

#### (ii) Subsequent wave interactions until the onset of the Friedlander wave

The head wave exits the first buffer region at point (*x*_{b},*t*_{b}), where *x*_{b} is the location where the head wave meets the tail wave entering the first buffer region (figure 3*a*). Thus, *x*_{b} is the distance travelled by the tail wave at constant speed *u*_{p}−*a*_{3} from 0 to *t*_{b}: *x*_{b}=(*u*_{p}−*a*_{3})*t*_{b}. Subsequently, the head wave enters a uniform region characterized by the flow behind the contact interface, region 3 in figure 1*b*, and travels at constant speed *u*_{p}+*a*_{3} along the characteristic: *x*=(*u*_{p}+*a*_{3})(*t*−*t*_{b})+*x*_{b}. The head rarefaction wave meets the contact interface moving at constant speed *u*_{p} along the line: *x*=*u*_{p}*t*. The intersection point is found to be:

The part of the head rarefaction wave transmitted through the contact interface enters a uniform region fully characterized by the flow behind the shock wave, region 2 in figure 1*b* and travels at constant speed *u*_{p}+*a*_{2} along the characteristic: *x*=(*u*_{p}+*a*_{2})(*t*−*t*_{c})+*x*_{c}. It meets the shock wave travelling at constant speed *W* on a straight line of equation: *x*=*Wt*. The intersection point (*x*_{d},*t*_{d}) gives the position and time at which the head rarefaction wave reaches the shock front, defining the onset of the Friedlander wave forming in the shock tube (figure 1*f*):
*t*_{b} from equation (2.27), this simplifies to:
*W*, *u*_{p} and *a*_{2} are, respectively, the shock speed, the speed of the contact discontinuity and the sound speed behind the shock, which are all determined from the Riemann problem as:
*a*_{1}, *γ*_{1}, *p*_{1} and *ρ*_{1} are, respectively, the sound speed, the specific heat ratio, the pressure and the density in the (undisturbed region of the) driven section.

It should also be kept in mind that equation (2.30) only applies to helium and air. Finally, the sought minimal total length required for a Friedlander wave to form in a shock tube is given by:
*L*,*p*_{2}) cannot be expressed explicitly in terms of driver section parameters (*L*_{1},*p*_{1},*T*_{1},*γ*_{1},*p*_{4},*T*_{4},*γ*_{4}) as, in their most simplified forms, equations (2.30) and (2.31), still depend on the value of the shock pressure *p*_{2} which can only be computed implicitly in terms of problem parameters by the expression [18]:

In the next subsection, we provide semi-analytical estimates of the additional characteristics of the Friedlander wave at the moment it forms in the shock tube.

### (b) Pressure history at the onset location

The arrival of rarefaction waves trailing the head wave are responsible for the pressure decay subsequent to the formation of the Friedlander wave. Determining the pressure behind the shock wave in terms of the initial parameters in the driven section (speed of sound: *a*_{1}, pressure: *p*_{1} and density: *ρ*_{1}) is in general difficult because entropy is discontinuous across the shock front. However for weak to medium shock, characterized by a shock strength *z*=(*p*_{2}−*p*_{1})/*p*_{1} below 5, the jump in entropy, asymptotically equivalent to a third power monomial in *z*, may be neglected [24]. The jump in the Riemann invariant *s* is also *O*(*z*^{3}). In practice, this assumption amounts to neglecting the reflection of any rarefaction wave off the shock front. Further, the pressure across a shock wave agrees up to second order in the shock strength with the pressure across a simple wave given by:
*u*—continuously modified by the incoming rarefaction waves—up to second order [19], leading to:

#### (i) Evolution of the local gas speed at the onset location

The computation of the local gas velocity *u* at *L* involves tracking the path of the incoming rarefaction waves. From the time it exits the first buffer region to the time it arrives at the Friedlander onset location, any rarefaction wave travels through the contact discontinuity (figure 2). Because the gas on both sides of this interface are in different thermodynamical states, part of the incoming rarefaction wave is transmitted through, and part of it is reflected off, the contact interface. The reflected part then interacts with the trailing incoming rarefaction waves in a second simple wave interaction region: the second buffer region (figure 4). Unlike the first buffer region, this region consists of non-centred rarefaction waves interacting with a moving boundary, which precludes a simple mathematical treatment. Therefore, we consider only a finite number, *N*+1, of equally spaced left-going rarefaction waves in the initial fan. The velocities of two successive rarefaction waves differ by a fixed increment Δ*U*=(*a*_{4}−(*a*_{3}−*u*_{p}))/*N*; the velocity of the *i*th rarefaction wave is then: *i*th rarefaction wave, the local gas speed and the local sound speed can be expressed in term of *s*_{i}:
*i*th rarefaction wave *r*_{i}=−*s*_{i} reflected off the left end of the tube is:

Due to the discontinuity in sound speed, the Riemann invariants along left-going characteristics are not constant across the contact interface. Ahead of this interface, all these characteristics originate from the uniform region between the contact interface and the shock front (region 2 in figure 1*b*) and share a common Riemann invariant *s*^{R} (see also figure 4). On the right side of any point *c*_{i0} on the contact interface shown in figure 4, *s*^{R} can also be evaluated from the sound speed ahead of the interface *a*^{R}_{ci0} and the local gas speed *u*_{ci0} preserved across the interface:
*c*_{i0} to *c*_{00} where the head rarefaction wave meets the contact interface. To that end, we express the sound speed in terms of temperature: *a*_{3} and the sound speed behind the shock front, *a*_{2}, respectively. The sound speeds on both sides of the contact interface at *c*_{i0} are related by:
*γ*_{1}=*γ*_{4}, (2.39) gives us an explicit constraint on the relationship of the local gas speed and sound speed at all the points on the contact interface.

Expressing the Riemann invariant *r*_{i} (2.36) of the incoming right-going characteristic intersecting the contact discontinuity in *c*_{i0} provides an additional relation between *u*_{ci0}:
*u*_{ci0} which, in turn, is the local gas speed when the *i*th rarefaction wave arrives at the Friedlander onset location, since this rarefaction wave transmitted through the contact discontinuity travels at constant speed until it meets the shock wave. Therefore, the pressure at the onset location *L* (figure 2) due to the arrival of the *i*th rarefaction wave can be estimated using (2.35) as:

#### (ii) Successive arrival time of the rarefaction wave at the onset location

In order to complete the description of the pressure history, the arrival times of the successive rarefaction waves at the blast onset location are needed, which, in turn, requires computing the time and position of all the wave intersections from the moment the *i*th rarefaction wave exits the first buffer region.

As shown in figure 3*a*, rarefaction waves reflecting off the closed end of the shock tube exit the first buffer region when they intersect the tail wave. The explicit analytical solution (2.16) determines the exact time at which this happens. The exit speed of the *i*th rarefaction wave is given by *U*^{+}_{bi}=*u*+*a*, where *u* and *a* are determined from the Riemann invariants *r*_{i} and *s*_{t} using (2.4). To compute the position at which each rarefaction wave exits the first buffer region, we assume that any rarefaction wave travels on a straight path and at the average velocity between any two successive wave intersections. For instance in figure 2, paths: *i*th rarefaction wave exits the first buffer region at time *t*_{bi} with wave speed *U*^{+}_{bi} at position *x*_{bi}. *x*_{bi} can be estimated from the values corresponding to the immediately preceding rarefaction wave: exit position *x*_{bi−1}, exit time *t*_{bi−1} and speed of the tail rarefaction wave *U*^{−}_{bi−1} at *x*_{bi−1} as follows:
*t*_{b} (2.27) and *x*_{b0}=(*u*_{p}−*a*_{3})*t*_{b}, thus completing the recurrence relation (2.42).

From the moment it enters the second buffer region, a right-going rarefaction wave trailing the head wave meets at least one of the discrete left-going rarefaction waves reflecting off the contact discontinuity. Let us consider a right-going rarefaction wave, which at point (*x*_{1},*t*_{1}) is travelling at speed *U*^{+}_{1} and a left-going rarefaction wave, which at point (*x*_{2},*t*_{2}) is travelling at speed *U*^{−}_{2} (figure 5*a*). These two waves meet at point (*x*_{3},*t*_{3}) provided that *x*_{1}≤*x*_{3}≤*x*_{2}, *t*_{1}≤*t*_{3} and *t*_{2}≤*t*_{3}. At this location, the right-going rarefaction is travelling at speed *U*^{+}_{3} and the left-going rarefaction wave travelling at speed *U*^{−}_{3}. We make the assumption that, between (*x*_{1},*t*_{1}) and (*x*_{3},*t*_{3}), the right-going rarefaction wave travels at average speed *U*^{+}_{13}=(*U*^{+}_{1}+*U*^{+}_{3})/2 along the line: *x*=*U*^{+}_{13}(*t*−*t*_{1})+*x*_{1} and that, between (*x*_{2},*t*_{2}) and (*x*_{3},*t*_{3}), the left-going rarefaction wave travels at average speed *U*^{−}_{23}=(*U*^{−}_{2}+*U*^{−}_{3})/2 along the line: *x*=*U*^{−}_{23}(*t*−*t*_{2})+*x*_{2}. It is worth noting that *U*^{+L}_{3} is the speed of the right-going rarefaction wave on the left-hand side of the contact interface computed as the sum of the local gas speed *u*_{3} and the sound speed on the left of the interface

Eventually, any right-going rarefaction wave meets the contact interface and modify its speed. Let us consider a right-going rarefaction wave, which at point (*x*_{1},*t*_{1}) travels at speed *U*^{+}_{1} and the contact discontinuity, which at point (*x*_{2},*t*_{2}) is travelling at speed *u*_{2} (figure 5*b*). The rarefaction wave and the contact interface meet at point (*x*_{3},*t*_{3}) under the conditions: *x*_{1}≤*x*_{2}≤*x*_{3}, *t*_{1}≤*t*_{3} and *t*_{2}≤*t*_{3}. A decrease of speed *u*_{3}≤*u*_{2} of the contact interface results from this interaction. Here also we assume that, between (*x*_{1},*t*_{1}) and (*x*_{3},*t*_{3}), the right-going rarefaction wave travels at average speed *x*=*U*^{+}_{13}(*t*−*t*_{1})+*x*_{1} and that, between (*x*_{2},*t*_{2}) and (*x*_{3},*t*_{3}), the contact interface travels at average speed *u*_{23}=(*u*_{2}+*u*_{3})/2 along the line: *x*=*u*_{23}(*t*−*t*_{2})+*x*_{2}. The intersection point is found to be:

The *i*th right-going rarefaction wave exits the second buffer region at (*x*_{c0i},*t*_{c0i}) and travels at constant speed *t*_{di}:

The determination of the position and time needs to happen in incremental order because it requires the time and position of the previous intersection for both the right-going rarefaction wave and either the left-going rarefaction wave or the contact interface, respectively. The simplified discrete approach we propose for estimating each intersection point (*x*_{cij},*t*_{cij}) inside the second buffer region consists in computing the intersection of each discrete right-going rarefaction wave stemming from the first buffer region, first, with the contact interface and, then after reflection off the contact discontinuity, with its trailing rarefaction waves. The intersection between the left-going reflected rarefaction wave with the immediatelly trailing discretized rarefaction wave in the fan is computed first. We proceed sequentially to next immediately trailing discretized right-going rarefaction wave until the tail wave. The steps described above are repeated for all the discretized rarefaction waves starting from the head wave and progresing sequentially to the next immediatelly trailing discretized right-going rarefaction waves until the tail wave. This procedure is summarized in algorithm 2.

The procedure described above furnishes the position of the onset of the Friedlander wave and the history of the pressure decay at that position. In the next section, we present an approach to obtain the practical characteristics of the blast wave from the pressure history at the onset location.

#### (iii) Estimation of the blast wave characteristics

At the onset location, three characteristics of the blast wave are desired: the peak overpressure, the decay time and the impulse. The peak overpressure, Δ*p*, is the difference between the maximum pressure *p*_{2} and the atmospheric pressure *p*_{0}. The time interval during which the blast pressure is greater than the atmospheric pressure is the decay time, denoted as *τ*. The impulse, *I*, is the positive area under the pressure time curve. The impulse is calculated by integrating the pressure–time curve up to the decay time. We estimate the decay time and impulse from the *N*+1 data points (*t*_{di},*p*_{di}) of arrival times and pressures resulting from the arrival of the *N*+1 rarefaction waves , computed from equations (2.45) and (2.41) (see also figure 2). We noted that by the time the tail wave arrives at the onset location, the pressure has dropped below the atmospheric pressure, *p*_{0}. It is then possible to identify the last time *t*_{dk} for which the pressure *p*_{dk} is still above atmospheric and estimate the time at which the blast pressure equals the atmospheric pressure *t*_{atm} by simple linear interpolation. The decay time is then obtained by subtracting *t*_{atm} from *t*_{d}, leading to:
*I* is estimated by numerical integration of the data points from *t*_{d} to *t*_{atm}. In our particular case, we used the trapezoidal rule and compared with the Friedlander impulse obtained by integrating (1.1) up to the decay time, resulting in:
*α* can be obtained numerically as the root of a nonlinear equation derived from the (2.47).

## 3. Results

In this section, we exercise the model in a number of configurations for a wide range of shock tube input parameters. For each configuration, the model furnishes the peak pressure and the exact onset time and location of the Friedlander wave (algorithm 1). In addition, we compute the pressure history at the onset location for each case with the semi-analytical approach outlined in algorithm 2. In the predictions shown, we employ a total of four discrete rarefaction waves (*N*=3 in algorithm 2). These are the head wave, the tail wave and two intermediate rarefaction waves. The time evolution of the pressure profile at the onset location is then described by the four points: (*p*_{d},*t*_{d}), (*p*_{d1},*t*_{d1}), (*p*_{d2},*t*_{d2}) and (*p*_{d3},*t*_{d3}) as illustrated in figure 2. From this discrete description of the pressure history, we use the method described in §2b(iii) to determine the temporal characteristics of the blast: decay time, impulse and wave form parameter.

We evaluate the quality of these predictions by comparing our model results with detailed CFD simulations for all configurations considered in the case of a compressed-air shock tube. Further, for these particular configurations, we assess to what extent the Friedlander wave form provides a good description of the type of blast waves that form in a shock tube. The versatility of the semi-analytical model is demonstrated in the analysis of helium-driven shock tube configurations. Specifically, we use the model to explore the trade-offs in the shock tube design space including the influence of the driver gas. In an attempt to validate the model, we compare the results of the semi-analytical blast estimates with experimental data of an existing shock tube.

### (a) Verification of analytical estimates using numerical simulation in a compressed air shock tube

The analytical estimates of the various parameters of interest in the formation of a Friedlander wave in a shock tube derived above can be verified using numerical CFD analysis. To this end, we use the Virtual Test Facility [25] to conduct simulations of the evolution of the flow in a shock tube under different conditions and compare the analytical estimates of the quantities of interest with the values furnished by the numerical computations. Specifically, we consider 16 different shock tube configurations with four different driver initial pressures and four different driver lengths. The initial driver pressure values used in simulations are selected such that the blast wave peak incident pressures reach approximately 1.5, 2.0, 3.0 and 4.0 atmospheres. For each of these four cases, the following driver lengths are considered: 0.25, 0.5, 0.75 and 1.0 m. For convenience, the problem is modelled using a two-dimensional grid with unit height. The grid length is chosen so that it fully contains the blast wave at the end of the simulation, thus avoiding confounding effects at the mouth of the shock tube. Simulations are run until the pressure at the onset location decays below atmospheric values (onset of the negative phase of the blast). The total time of the simulation should be greater than the sum of the blast formation time *t*_{d} and the decay time *τ*, which can be obtained from equations (2.29) and (2.46), respectively. The initial conditions in the tube consist of the value of the pressure in the driver section, atmospheric pressure values in the rest of the domain, the initial temperature set to 300 K, the gas density computed from the ideal gas equation of state. The verification is done only for the case of an air-driven shock.

Figure 6 presents a comparison of the analytical and numerical values of the pressure history obtained at the location *L* where the analytical model predicts the onset of the blast wave for all 16 configurations in table 1. In all cases, we observe that the numerical and analytical predictions for the peak overpressure are within 2% of each other. In all cases considered, the three analytical predicted values, (*p*_{di},*t*_{di}), describing the time decay of pressure past the initial peak value fall almost on the top of the numerical curve. For low peak incident pressures: 1.5×10^{5} Pa and 2.0×10^{5} Pa, the pressure evolution with respect to time is accurately described by a piecewise linear interpolation between the four analytical points. For higher peak incident pressures: 3.0×10^{5} Pa and 4.0×10^{5} Pa, corresponding to a stronger shock regime, the pressure profile has higher curvature and consequently four analytical points are insufficient to fully capture the shape of the actual (numerical) pressure profile. However, the individual analytical points do capture the numerical values very accurately.

### (b) Comparison of the numerical blast wave pressure history with the Friedlander mathematical form

In this subsection, we want to quantify how closely the numerical pressure profiles are described by the Friedlander mathematical form. To this end, we fit the Friedlander wave form parameters to the positive phase of the pressure profiles predicted with the simulations as follows. The peak overpressure and the decay time in (1.1) are extracted directly from the numerical pressure profiles at the onset location. The wave form parameter *α* is obtained from the simulations by least-square regression of the numerical pressure history. This furnishes the optimal value of the wave form parameter *α*. The Friedlander wave form parameters for all of the 16 cases are listed in table 2. We do not plot the analytical wave forms as they almost exactly match the numerical profiles. In order to further quantify the match, we compare the impulse from the closed form integral of the analytical wave forms (2.47) with the numerical integral of the simulated cases using the trapezoidal rule. Table 3 shows the numerical values obtained as well as the error. As shown in table 3, the error between the numerically integrated impulse and the impulse obtained by a fit of the blast parameters of the Friedlander impulse ranges from 0.1 to 1.8%. The error increases with the shock strength and does not seem to vary significantly with the tested driver length. For all the tested cases, the error in impulse between the numerical model and the Friedlander mathematical form is less than 2%. This demonstrates that the numerical shock tube model produces Friedlander waves.

### (c) Quantitative comparison of the semi-analytical and numerical blast wave parameters

In §3a, we showed the semi-analytic solution gives an excellent estimate of the pressure history for all the points considered (*N*=4 in our case, figure 2). Here, we estimate the blast wave parameters following the procedure in §3. We compute (2.46) to find the decay time *τ* and make use of (2.47) to find the wave form parameter *α*. The computed blast parameters for the 16 cases are reported in table 2. It can be seen that the error in the quantities of interest increases with the strength of the shock due to the few points considered. For stronger shock intensities, the accuracy could be improved by adding more rarefaction waves to the semi-analytic model, of course as long as the isentropic shock assumption is valid.

### (d) Influence of the driver gas on the peak pressure and location of the onset of a blast

In this section, we take advantage of the versatility of the semi-analytic model and use it to investigate how the peak pressure and the blast onset location are affected by choosing helium as driver gas. Table 1 compares the results for air and helium-driven shock tubes for the 16 different shock tube configurations. We observe that helium-driven shock tubes have higher peak pressures which become relatively higher as the driver pressure increases. In the limit of low driver pressure, the peak overpressures match. In the highest driver pressure considered, the peak pressure is almost twice of the corresponding value for air-driven shock tube. For all the cases considered, the blast location for helium is about three times shorter than air.

### (e) Comparison with experimental data

In order to validate the proposed semi-analytic model, we would require pressure-versus-time data at different locations along the length of the driven section. Ideally, an array of closely spaced sensors would be located around the region where the model predicts the onset of the blast profile. This would allow to experimentally capture the transition from the plateau to the fully formed blast wave and posterior decay, and, therefore, an accurate experimental determination of the onset of the blast. This would then be quantitatively compared to the model predictions under different initial conditions and help validate the model. Unfortunately, we have been unable to find any previous work published in the literature where this particular type of experiment has been conducted. The value of the peak incident overpressure is easier to compare (and easier to compute). For this, it suffices to extract the shock pressure while the first rarefaction wave has not caught up with the shock front. The pressure value corresponds to the Riemann solution.

With these considerations in mind, we applied our model to the experimental conditions in Refs [3,26]. It is important to emphasize that this should not be construed as an attempt to validate the model, but merely to make at least some connection with published blast shock tube data.

The shock tube employed in Refs [3,26] measures 7.12 m and is driven by helium. The driver section has a length of 0.76 m and an initial pressure of 241 kPa. Three pressure sensors are placed along the driven section at distances of 2.80 m (the ‘R-wall’ sensor), 5.57 m (the ‘Trigger’ sensor) and 5.98 m (the ‘Pencil’ sensor) from the closed end of the tube. The time evolution of the pressure at these different locations is reported in [26]. At the ‘R-wall’ location, the pressure profile exhibits the short plateau indicating that the head rarefaction wave has not yet caught up with the shock wave, whereas the pressure profile at the two other locations exhibit fully formed blast wave pressure profiles. We estimate the experimental incident pressure by looking at the value of the plateau pressure at the ‘R-wall’ sensor. The reported experimental profile shows a mean plateau pressure of approximately 83 kPa with high-frequency oscillations of about 11 kPa. The analytical model predicts an incident overpressure of 84 kPa. The model also predicts that the blast wave should form 4.02 m away from the closed end of the tube. Lacking information on the experimental location of the onset of the blast, one can at most check for consistency between the model and the test. Specifically, we observe that the two pressure sensors which show fully formed blast pressure profiles are located downstream (respectively, 5.57 and 5.98 m from the closed end of the tube) of the onset location predicted by the model for this case (4.02 m). For convenience, these results are summarized in table 4.

## 4. Conclusion

In this paper, we presented a semi-analytical approach to predict the characteristics of pseudo blast waves that form in a shock tube from test configuration parameters. The expressions for the blast intensity and onset location are given in closed form as a result of finding the intersection point between the shock and the head rarefaction wave. To further compute the pressure history at that location, we use a semi-analytical approach in which discrete pressure–time values are obtained by tracking the evolution of subsequent rarefaction waves. We show that the blast decay time and impulse can be estimated with a few discrete pressure–time points.

The semi-analytical model was verified against numerical CFD calculations in 16 different configurations (driver length and pressure). It was also shown that both the semi-analytical model and the CFD results can be fitted to a Friedlander wave form, thus confirming the long-standing belief that shock tubes do indeed produce Friedlander-like blast waves.

Among some important practical conclusions of interest to shock tube operations derived from the theoretical analysis, we find that: (i) the maximum incident pressure *p*_{2} is not affected by the driver length *L*_{1}, as expected; (ii) the location of the onset of the blast *L* increases with the driver length *L*_{1}; (iii) for a fixed driver pressure *p*_{4}, the decay time *τ* and specific impulse *I* increase proportionally to the driver length, whereas for a fixed driver length, the decay time and impulse increase with increasing driver pressure; (iv) the wave form parameter is found to be independent of the driver length and to increase with driver pressure; and finally, (v) the formation of the blast wave occurs closer to the diaphragm in the case of helium-driven shock tubes.

It is to be noted that our analytical model that employs the theory of one-dimensional unsteady gas flow(that disregards the effects of viscosity and thermal conductivity [27]) is an idealized representation of the shock tube flow. For instance, the shattering of the diaphragm is a complicated process with the implication that the shock wave does not form instantly [18].

## Authors' contributions

J.F. and R.A.R. provided the motivation for this study. A.F.T. and M.H. derived the semi-analytical model, implemented and ran the numerical CFD model. A.F.T., M.H. and R.A.R. drafted the manuscript. All the authors interpreted the results and gave final permission for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the US Army Research Laboratory and the US Army Research Office through the Institute for Soldier Nanotechnologies, under contract no. W911NF-13-D-0001.

- Received August 28, 2015.
- Accepted December 24, 2015.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.