## Abstract

The existing practice of designing air-blast-resistant structures relies on the ideal gas model. But this model predicts the maximum value of the reflection coefficient (ratio of the reflected to the incident pressure) to be 8, whereas it can go up to 20 or more as reported in the literature. To address this discrepancy, air medium is modelled as a real gas instead of an ideal gas, where the effect of intermolecular forces, vibration, dissociation, electronic excitation and ionization are included. Ranges of peak over-pressure are identified where the ideal gas assumption cannot be used. Differences in impulse transmitted to the free-standing plates of different mass owing to relaxing of the ideal gas assumption and consideration of the real gas model are evaluated. Impulse transmitted to the structures for constant and variable back pressure (VBP) is also compared considering the real gas model. The result shows that for high-intensity shock, the ideal gas model under-predicts impulse transmitted to heavy plates but over-predicts the same for light-weight plates. Impulse transmitted to light-weight plates is also overestimated if VBP is neglected. The implications of this research are substantial for designing high-intensity air-blast mitigating structures, which if not considered properly, may lead to compromise in structural performance.

## 1. Introduction

The beneficial effect of fluid–structure interaction (FSI) in mitigating explosion-induced shock wave loads was identified by Taylor [1]. Taylor's theory [1] considers medium to be acoustic and thus limits its applicability to low pressure regimes (over-pressure less than 1 atm) or far-field analysis where the effects of nonlinear compressibility of medium can be ignored [2–4]. In the study of shock waves in air [2,5–8] or in water [9,10], it is observed that consideration of nonlinear compressibility over acoustic approximation for the medium reduces the impulse transmitted to a lightweight structure, but increases in the case of a heavy structure.

In all these theories, the guidelines for designing a blast-resistant structure are primarily governed by a parameter called the reflection coefficient (*C*_{R}), which is the ratio of reflected pressure to the incident pressure. Kambouchev *et al.* [2], considering the ideal gas equation of state (EOS) for air, showed that the reflection coefficient has a minimum value of 2 (pertaining to acoustic limit) and a maximum value of 8. It should be noted that other theories [6–8] for air-explosion-induced shock wave loading were also derived using the ideal gas assumption. However, in the literature [11,12], it is reported that the value of the reflection coefficient may even go up to 20 or more. Brode [12] reported that for a peak over-pressure higher than 6.9 atm, the reflection coefficient measured by the ideal gas EOS is increasingly in error and should not be used. Guzas & Earls [13] also restricted the application of ideal gas EOS to a peak over-pressure below 6.9 atm in their proposed methodology for calculating load on structures subjected to air blast. Baker [11] reported that non-ideal behaviour, such as dissociation and ionization of gases, is responsible for such high values of for the reflection coefficient. Details of these experiments designed to measure the reflection coefficient can be found in Baker [11], pp. 8,9 and 130–132.

With the increase in terrorist activities around the world, it is imperative that better mitigation measures are found for design of infrastructure. Following the terrorist bomb attack on the Alfred P. Murrah building in downtown Oklahoma City, USA, the Federal Emergency Management Agency (FEMA) [14] prescribed a minimum stand-off distance of 10 ft (3.048 m) for an explosive charge of 15 kg and 70 ft (21 m) for an explosive charge of 4535 kg. In this regard, it should be mentioned that a large truck can typically contain 4535 kg (10 000 lb) or more equivalent of TNT, whereas vehicle bombs employing sedans to vans can contain 227–1814 kg (500–4000 lb) of TNT equivalent explosive material [14]. Vehicle bombs nearer than the threshold stand-off distance prescribed by FEMA can easily result in a situation where the incident shock strength can exceed even 500 atm, as estimated by the Kingery–Bulmash blast indicator (based on data from explosive tests using charge weights from less than 1 to over 400 000 kg and having a range of applicability of 0.05–40 m). Protection against such high-intensity blasts often requires sacrificial foam claddings, composite panels, fibre-reinforced polymer retrofit system and/or mixed concrete, etc., to lessen the effects. In the analysis of these structures, the loading curves are estimated from the FSI theories [2,5–8] to avoid complexities associated with interaction between fluids and structures. For example, Vaziri & Hutchinson [15], Main & Gazonas [16], Aleyaasin *et al.* [17] used empirical formula of Kambouchev *et al.* [2] to calculate load on structures for FSI analysis. However, for the shock strengths considered in these works (200 atm by Vaziri & Hutchinson [15] or 600 atm by Aleyaasin *et al.* [17]), the ideal gas assumption is no longer valid and may lead to possible discrepancies in designing blast-mitigating structures.

Typically, in an event of near-field, high-intensity explosion in air, the gas is strongly compressed and heated up by the shock wave, resulting in extremely high temperatures and pressures [18]. It is well known from the literature [18–20] that at high pressure and temperature, the ideal gas assumption (applicable below 1000 K) loses its basis because it only considers translational and rotational degrees of freedom of molecules. At high pressure and temperature intermolecular forces, vibration, dissociation, electronic excitation and ionization become important physical characteristics of a medium. So, the physics of high-intensity shock wave propagation in air demands development of a model which is able to capture these effects at high pressure and temperature. Thus, a real gas EOS is considered in this manuscript for modelling air medium.

Although the literature in shock wave propagation in real gases is abundant [21–25], these have not been used for determination of reflected pressure as well as determination of impulse transmitted to the structure with an objective of constructing performance charts for blast-resistant design. In all the existing theories [1,2,5–8] on the basis of which the performance charts for air-blast-resistant designs are prepared, modelling of air is done as a polytropic fluid (specific heat at constant volume is constant) which considers *only* translational and rotational motion of molecules to be responsible for internal energy with its concomitant effect restricting the value of *C*_{R} to 8 and not 20 or more as reported in the literature [3,11]. Thus, it can be argued that current performance charts of air-blast-resistant design of structures based on acoustic or ideal gas assumptions are questionable [3] and an upper limit of *C*_{R}≤8 may serve as a bottleneck for design of structures.

In this manuscript, an effort is made to highlight the necessity of preparing design guidelines for mitigation of high-intensity explosion in air considering real gas behaviour. So, a non-polytropic van der Waals EOS is used to model shock wave propagation. This includes the effect of not only the intermolecular forces, but also molecular vibration, dissociation, electronic excitation and ionization in addition to translational and rotational motion [18]. However, reaction between gases, coupling of rotational and vibrational degrees of freedom and anharmonicity of vibration are ignored in this manuscript. A seven-species air system (N_{2}, O_{2}, N, O, N^{+}, O^{+} and e^{−}) is modelled with reference composition 79% N_{2} and 21% O_{2} by volume.

The outline of the manuscript is as follows: existing FSI theories are discussed in §2. Section 3 contains details of the real gas model. A problem statement along with shock wave modelling methodology in a real gas is provided in §4. Section 5 discusses the effects of considering the real gas over ideal gas assumption. Finally, in §6, conclusions are discussed.

## 2. Existing fluid–structure interaction theories

Here, the state-of-the-art models of FSI are summarized. It should be noted that this manuscript deals only with constant shock, and thus FSI theories related to constant shock are discussed.

### (a) Acoustic models

Let us consider a submerged rigid plate of mass per unit area *m*_{p}=*ρ*_{p}*h*_{p} (*ρ*_{p} and *h*_{p} represent density and thickness of the solid plate, respectively) having two wetted surfaces on the left or front side and right or back side of the plate. A plane pressure wave is incident on the front side with *p*_{s} as the peak over-pressure, i.e. *p*^{i}(*x*,*t*)=*p*_{s}. The incident wave is partially reflected back owing to translation of the plate. Pressure on the front side of the plate, *p*_{f}, can be obtained as the sum of the incident (*p*^{i}) and the reflected wave (*p*^{r}),
*φ*(*t*) is the unknown functional for the reflected wave. Balance of mass and momentum yields the velocity of the plate (*p*_{s}/*ρ*_{0}*c*_{0} and *p*_{s}*φ*(*t*)/*ρ*_{0}*c*_{0} represent the particle velocities owing to incident and reflected waves, respectively. Here, *ρ*_{0} and *c*_{0} are the initial density and sound speed of the medium. Combining equations (2.1) and (2.2) and eliminating *φ*(*t*), we obtain
*p*_{f} on the front side and *p*_{b} on the back side of the plate. So, the equation of motion of the plate can be written as

#### (i) Taylor's model [1]

Taylor's model assumes that the pressure on the back side of the plate remains constant, i.e. *p*_{b}(*t*)=0. So, the equation of motion according to Taylor's theory can be written from equation (2.3) and (2.4) as
*u*(0)=0 and *t* gives front-side pressure as
*t*_{i} is the duration of pressure pulse. FSI parameter for the constant back pressure (CBP) case is given as

#### (ii) Liu–Young model [4]

In the Liu–Young (LY) model [4], it is assumed that the receding motion of the plate will cause pressure to rise at the back side of the plate. Thus, mass and momentum balance on the back side of the plate give
*p*_{f} with equation (2.3) and *p*_{b} with equation (2.7) yields
*C*_{R}) is always 2.

### (b) Fluid–structure interaction theories based on ideal gas model

Typically, these theories consider nonlinear compressibility of the medium. A system of fluid is characterized by two equations of state: caloric and thermodynamic. In an ideal gas model, the caloric EOS is
*ρ* is the current density of the fluid, *T* is the temperature, *p* is the pressure, *e* is the internal energy, *R* is the specific gas constant and *c*_{v} is the heat capacity at constant volume. For a polytropic fluid, the contribution of only translational and rotational degrees of freedom to internal energy is considered and thus *c*_{v} turns out to be constant. This yields a value of 1.4 for the isentropic exponent *γ*, the ratio of specific heat at constant volume (*c*_{v}) to the specific heat at constant pressure (*c*_{p}).

#### (i) Kambouchev–Noels–Radovitzky model [2]

The Kambouchev–Noels–Radovitzky (KNR) model uses ideal gas EOS to account for nonlinear compressibility. This FSI theory asymptotically satisfies Taylor's theory and presents an empirical expression to estimate transmitted impulse as
*β*^{CBP}_{s} is the FSI parameter using CBP assumption, *C*_{R} is the reflection coefficient, *I*_{p} is the transmitted impulse and *I*_{i} is the incident impulse.

The FSI parameter can be expressed as

Shock velocity (*U*_{1}) can be obtained from the formula
*c*_{0} and *p*_{0} are the speed of sound and pressure in ambient conditions.

*ρ*_{1} is the density downstream of incident shock and can be expressed as
*f*_{R} (refer equation (2.13)) can be written as
*C*_{R}) is given as

Once the relative impulse transmitted to the plate is evaluated, the design pressure to be applied to the plate on the front side can be calculated as [15,17]
*et al.* [2] used the same constant back pressure (CBP) assumption as used by Taylor, i.e *p*_{b}=0.

#### (ii) Peng–Zhang–Gogos–Gazonas model [7]

This model removes the CBP assumption of the KNR model. The flow behind the plate is modelled by considering the resistance provided by air on the back side of the plate. The pressure rise on the back side of the plate can be calculated analytically from the receding velocity of the plate (*u*_{b}) as

Here, it should be noted that for the theories based on ideal gas the reflection coefficient has a range 2≤*C*_{R}≤8, where the lower bound corresponds to the linear acoustic range.

## 3. Real gas model

In this manuscript, the van der Waals EOS is used to model air as a real gas, consisting of 79% nitrogen and 21% oxygen, in which both the gaseous constituents are symbolized as *A*_{2} in their undissociated state (*A*_{2}∈*N*_{2},*O*_{2}). In the van der Waals EOS, pressure (*p*) of real gas is expressed as [26,27]
*a*_{A2} is the pressure/attraction correction and *b*_{A2} is the volume/repulsion correction of the constituents. Here, *α*_{A2} is the volume fraction of the constituents. *Z*_{r})_{A2} is the compressibility factor for any molecular species and can be expressed as a function of mass fraction of dissociated molecules (*α*_{d}) and ionized atoms (*α*_{1})
*α*_{d}) as [28]
*ρ*_{d} and *T*_{d} are the characteristic density and temperature of a constituent gas, respectively.

Mass fraction of ionic species can be obtained from equilibrium of ionization (*α*_{1}) and is given by the modified Saha equation [19] as
*T*_{1} is the characteristic temperature of ionization of a constituent and *c*_{1} is a constant given as
*m*_{e}, *k* and *h* denote the electron mass, Boltzmann's constant and Planck's constant, respectively.

Equation (3.1) can also be written as
*Z*_{c})_{A2} is the expression for the compressiblity factor found in the text books on thermal science where the presence of intermolecular forces (attraction and repulsion) are only considered. In this work, compressibility is updated to introduce ionization and dissociation. So the compressibility factor can be redefined as

The caloric EOS can be written as a sum of partition contribution of each degree of freedom to the internal energy [23,24]
*e*_{trans} is the internal energy contribution owing to molecular translation and is expressed as
*e*_{rot}) is given as
*T*_{R} is the characteristic temperature for rotation.

The vibrational energy of a diatomic molecule, considering it as a harmonic oscillator, can be approximated as
*T*_{vib} is the characteristic temperature of vibration.

Internal energy owing to compressibility of van der Waals gas is given as

Contributions of electronic excitation to internal energy (*e*_{elec}) are due to the individual contribution of undissociated molecules (*e*_{elA2}), atoms produced after dissociation (*e*_{elA}) and ionized atoms (*e*_{elA+}). It can be expressed as
*e*_{elA2}, *e*_{elA}, *e*_{elA+}) is given by a general formula as
*g*_{m} refers to statistical weight or degeneracy factor of the energy levels and *T*_{m} is the characteristic temperature of the corresponding degeneracy level. *T*_{m} can be obtained from the relation, *T*_{m}=*e*_{m}/*k*. *e*_{m} is the internal energy of any electronic level and can be obtained from spectroscopic measurements.

Internal energy contributions to cause dissociation of molecules and ionization of atoms are given, respectively, as
*T*_{d} and *T*_{1} are the characteristic temperature owing to dissociation and ionization, respectively.

The standard values of constants in equations (3.1)–(3.16), as obtained from numerous references, are assembled and provided in table 1 (see appendix A).

## 4. Problem statement

The problem of uniform/constant normal shock propagation and reflection from a fixed and movable rigid wall is considered in this section. A semi-analytical model with real gas EOS is presented which is also used to validate a numerical model. Impulses transmitted to the free-standing plates for different plate mass, shock intensity and backing conditions are calculated to highlight the differences among existing theories which use ideal gas EOS. Given the primary objective of the manuscript is to clarify the discrepancy between the maximum value of reflection coefficient reported and the one used in existing FSI theories, consideration of exponential shock loading (typically used to represent blast loading) is outside the scope of this work.

### (a) Semi-analytical solution of uniform shock reflection from fixed rigid wall

The solution of constant shock wave reflection from a fixed rigid wall gives reflected over-pressure (*p*_{r}) for a given incident over-pressure (*p*_{s}) or in other words the reflection coefficient (*C*_{R}=*p*_{r}/*p*_{s}). This requires solving of Rankine–Hugoniot–Jump (RHJ) conditions [18] along with thermal (equation (3.1)) and caloric (equation (3.9)) EOS. Assuming an inviscid medium, RHJ conditions (which are mass, momentum and energy conservation across the discontinuity/shock front) can be written as
*U* represents the material shock velocity and *u* represents the particle velocity. Flow state ahead and behind the shock are represented by index *i* and *j,* respectively. Internal energy, pressure, temperature and fluid density are symbolized by *e*, *p*, *T* and *ρ*, respectively. Three different cases are considered to study the implication of the energy contribution of different degrees of freedom on reflection coefficient. In this regard, the energy contribution of different degrees of freedom to the caloric EOS are varied. However, the thermodynamic EOS is the same in all the three models.

*Model I*: this model considers the contributions from the translational and rotational degrees of freedom to the internal energy. In addition to these, contributions from the compressibility of a van der Waals gas are also considered. Thus, the caloric equation of state in this case turns out to be
*e*_{comp} in the reflection coefficient in comparison with that of an ideal gas which only considers translational and rotational degree of freedom.

*Model II*: in this model, the contribution from vibrational degree of freedom is considered in addition to the degrees of freedom already considered in model I. So, in this model, the caloric EOS is
*e*_{vib} on the reflection coefficient.

*Model III:* contributions from all degrees of freedom are considered in model III and hence the caloric EOS is the same as equation (3.9). This model will be able to highlight the effect of ionization and dissociation on the reflection coefficient, when compared with model I or II.

Here, it should be noted that the caloric EOS contains transcendental functions and unlike the ideal gas case, RHJ conditions cannot be solved directly. Thereby, a Newton–Raphson iterative scheme is used. The solution procedure requires, solving of the RHJ condition along with both caloric and thermodynamic EOS applied between shocked ( *j*=1 in equation (4.1)) and unshocked (*i*=0 in equation (4.1)) zones, to determine incident shock quantities (*ρ*_{1},*U*_{1},*u*_{1},*e*_{1},*T*_{1}) for a given peak over-pressure (*p*_{s}=*p*_{1}−*p*_{0}). Air properties of the unshocked zone are chosen to be that of the ambient condition, with the upstream particle velocity (*u*_{0}) taken as zero. Once the incident shock quantities are obtained, the reflected shock quantities (*ρ*_{2},*U*_{2},*e*_{2},*p*_{2},*T*_{2}) can also be obtained by again solving the RHJ condition along with the EOS (both caloric and thermodynamic), applied between reflected ( *j*=2 in equation (4.1)) and incident (*i*=1 in equation (4.1)) shock. Wall boundary is modelled by considering particle velocity along the downstream of reflected shock to be zero, i.e. *u*_{2}=0. For the first peak over-pressure, ideal gas solutions (both incident and reflected) are used as the initial estimate for the solution, whereas from next step onwards the initial estimate is taken as the converged solution of the previous step. The converged solution satisfies a set of equations within the specified tolerance, i.e. the norm of the residuals should be less than 10^{−6}.

### (b) Numerical model

Governing equations for an inviscid compressible fluid in a one-dimensional infinitesimal closed thermodynamic system are solved considering a Lagrangian framework. An artificial viscosity based on the von Neumann–Richtmyer algorithm is used for shock capturing. Artificial dissipative terms are introduced into the governing equation to smear out the discontinuity/shock over a certain width (generally taken as grid spacing of the spatial discretization which is negligible compared with its domain length) which eventually converts the discontinuous RHJ equations to a set of continuous governing equations valid for the entire domain. Compressibility of the medium is modelled using both ideal gas and real gas EOS to study the difference in impulse transmitted to the plates of different mass subjected to various shock intensities. Different backing conditions of the plate (CBP and VBP) are also considered in this study.

#### (i) Governing equations

*Kinematic relation:* material velocity (*V* _{x}) and acceleration (*A*_{x}) are given as
*x* and *X* represent the position of any fluid particle in the current and reference configuration, respectively, and *t* is time.

*Mass conservation*: mass conservation in one space dimension (no change in area between the undeformed and deformed state), according to infinitesimal deformation theory, can be written as
*ρ*_{0} and *ρ* are the initial and current density of any particle in the fluid, respectively.

*Momentum conservation*: an artificial viscosity is subtracted from pressure for numerical stability, and thus linear momentum conservation reduces to
*p* is the pressure and *A*_{x} is the acceleration measured in a material coordinate system. *Q* is artificial viscosity and can be expressed as
*K*_{1} and *K*_{2} are the artificial viscosity coefficients. Δ is the width of the smeared shock which is approximated as one grid spacing to avoid numerical instabilities. *D* is the rate of deformation and can be expressed as *F*, can be expressed as *F*=∂*x*/∂*X*. *F*, i.e. *c*) in an ideal gas is
*c*) in a real gas can be obtained from its definition as
*H*=*e*+*p*/*ρ*, represents specific enthalpy. The detailed derivation of the expressions *χ*=(∂*e*/∂*ρ*)_{T}, *κ*=(∂*H*/∂*p*)_{T}, *c*_{p}=(∂*H*/∂*T*)_{p} and *c*_{v}=(∂*e*/∂*T*)_{ρ} is provided in appendix B.

*Energy conservation*: in order to keep the formulation simple, flow is assumed to be non-conductive and non-radiative. So, the energy conservation equation in one space dimension becomes
*e*=*e*(*ρ*,*T*). Thereby, differentiation of the term with respect to time (*t*) yields

*Equation of state*: both ideal gas and real gas EOS are used in this manuscript for comparative study.

Ideal gas: combining the thermodynamic EOS (equation (2.12)) and the caloric EOS (equation (2.11)) for an ideal gas yields

Real gas: in order to model air as a real gas, the thermodynamic and caloric EOS presented in §3 are used. In numerical calculations, the contribution from all the degrees of freedom is taken into account. So, the numerical model and mode lIII of semi-analytical calculation (§4*a*) are comparable.

#### (ii) Discretization of governing equations

Discretization of governing equations is done by the finite difference method proposed in a seminal work of von Neumann & Richtmyer [29]. The domain is discretized into equally spaced *M*+1 grid points. Coordinates of the *m*th grid point can be written as *X*_{m}=(*m*−*M*)Δ*X*, where Δ*X* is the grid spacing. Temporal discretization (*t*_{j}=*t*_{j−1}+Δ*t*_{j}, *j*=0,1,2,…) is done with variable time-steps (Δ*t*_{j}), calculated based on CFL condition as
*α* is a regulation parameter and typically varies within the limits 0<*α*≤0.5, depending on shock strength and grid size. An explicit Newmark time integration scheme is used for time marching. Detail of numerical method and discretization of the equations (equation (4.4)–(4.7), (4.8), (4.10) and (4.13)) to model air as an ideal gas, can be found in reference [2] (see appendix of reference [2]) and reference [29] (see section VI of reference [29]). However, for a real gas model, implementation of the thermodynamic EOS (equation (3.1)), energy equation (equation (4.12)) and sound speed (equation (4.9)) is different and is as follows.

*Thermodynamic equation of state:*
*α*_{d})^{m+1/2}_{j}]_{A2} and [(*α*_{1})^{m+1/2}_{j}]_{A2}, at any spatio-temporal location, for a given thermodynamic state (

*Energy conservation*: temporal update of temperature at any spatial grid location can be obtained from the energy equation (equation (4.12)) as

*Sound speed*:
*c*_{p}, *c*_{v}, *χ* and *κ*, which are described in detail in appendix B.

#### (iii) Dynamics of the plate

Assuming the free-standing plate to be rigid, it is modelled as a point mass. Acceleration (*A*_{p}) of the plate is determined from the pressure difference between the front and back of the plate as
*m*_{p} is the mass of the plate. *p*_{f} is the front-side pressure and taken as the pressure computed in the mid-position of the last fluid element, i.e. *p*_{f}=*p*^{M−1/2}_{j}. Calculation of pressure on the rear-side of the plate (*p*_{b}), for different backing conditions, is discussed in the following paragraph. Motion of the plate is computed by time integrating equation (4.19) with forward Euler scheme.

#### (iv) Boundary condition

*Fixed rigid plate or wall boundary:* this condition is implemented by maintaining velocity and/or acceleration of the node corresponding to the plate as zero throughout the simulation (i.e. *V*_{x}^{(M+1)}_{j+1/2}=0 and/or *A*_{x}^{(M+1)}_{j+1/2}=0 for *j*=0,1,2,…).

*Free-standing plate:* in this manuscript, two choices of backing conditions are solved: CBP and VBP.

*CBP:* In CBP, pressure at the back side of the plate is kept at constant atmospheric pressure, i.e. *p*_{b}=*p*_{0}. It is worth mentioning that for an ideal gas and a real gas the CBP condition is same.

*VBP:* in VBP, it is assumed that the pressure rises at the back side owing to receding motion of the plate. Air at the back of the plate is modelled as a real gas and pressure rise at that side is estimated by the shock jump relation as
*U*_{b} and *u*_{b} are the shock and particle velocity of fluid at the back side, respectively. Here, it is assumed that the receding velocity of the plate and the adjacent fluid are same. In a real gas, a linear fit to the shock-particle (*U*_{b}=*c*_{0}+*s*_{1}*u*_{b}) velocity relation is required for calculating shock velocity and hence back-side pressure. In order to carry out the calculation for the VBP case, fitting coefficients of the shock-particle velocity relation are required. In this regard the RHJ condition (equation (4.1)) along with the thermal (equation (3.1)) and caloric (equation (3.9)), EOS are solved using the Newton–Raphson method as mentioned in the previous subsection (considering upstream particle velocity to be zero) to obtain shock velocity for a given downstream particle velocity. A linear fit to shock-particle velocity data yields the fitting coefficients to be *c*_{0}=346 and *s*_{1}=1.045.

The numerical algorithm for both these backing conditions can be found in [9] (see appendix in [9]).

#### (v) Initial condition

In this numerical model, a one-dimensional domain is chosen, whose left endpoint is the source point of explosion, right endpoint is the plate and the points in between represent the rest of the domain which is air. Because, in this manuscript, a constant shock is modelled, a constant over-pressure of intended shock strength is provided at the left endpoint throughout the simulation, i.e. *p*_{L}(*X*=0,*t*≥0)=*p*_{s}, where *p*_{s} is the desired over-pressure downstream of the shock. Initial conditions for the plate are displacement *x*_{p}(*X*=*X*_{M},*t*=0)=0 and velocity *V* _{p}(*X*=*X*_{M},*t*=0)=0. For unshocked air particles, the atmospheric condition is used as initial. This yields *x*(*X*,*t*=0)=0, *p*(*X*,*t*=0)=*p*_{0}, *T*(*X*,*t*=0)=*T*_{0}, *ρ*(*X*,*t*=0)=*ρ*_{0} and *V* _{x}(*X*,*t*=0)=0.

## 5. Results and discussion

Here, results based on various models (both the semi-analytical and the numerical model, demonstrated in §4*a* and *b,* respectively) are presented and compared. Validation of the numerical model is done against model III of the semi-analytical calculations. Transmitted impulse, using both the ideal and real gas EOS, is calculated and compared for constant shocks of different strength impinging on plates with various mass and backing conditions. In the numerical calculation, a spatial domain of 7.5 m is modelled. It is discretized with equally spaced *M*=5000 grid points. Here, it should be noted that the domain size is chosen in such a way as to avoid any unwanted reflection. The initial density of air, *ρ*_{0}=1.225 kg m^{−3}, and temperature, *T*_{0}=298 *K*, are chosen as ambient conditions in this manuscript for all calculations. This yields a corresponding ambient pressure (*p*_{0}) of 1.048 bar and initial sound speed (*c*_{0}) of 346 m s^{−1}. Depending upon shock intensity artificial viscosity parameters *K*_{1} and *K*_{2} are varied between 0.05–1.0 and 0.2–2.3, respectively.

### (a) Comparison of different models

The reflection coefficient calculated for planar shock reflection from fixed rigid boundary is compared for various models presented in §4*a* along with the ideal gas model. It can be clearly seen from figure 1 that model I, in which only the compressibility correction term is added in caloric EOS in comparison with the ideal gas model (which only considers translation and rotational degrees of freedom), gives more or less the same reflection coefficient for various peak over-pressure ratios with a maximum value equal to 8. In model II, the contribution of vibrational energy is considered in addition to the compressibility correction term in model I. Figure 1 clearly demonstrates that for peak over-pressure ratios greater than ∼7, the contribution from vibrational degrees of freedom to internal energy is prominent and cannot be neglected. This also suggests that for a shock strength exceeding 7 bar, the ideal gas EOS cannot be used. Here, it should be noted that this limit of applicability of the ideal gas EOS is in line with the suggestion of Brode [12] and Guzas & Earls [13]. However, model II can predict the maximum value of reflection coefficient to be approximately 10. In model III, the contribution from ionization and dissociation are also included. This gives the value of the reflection coefficient as 16.1 for a peak over-pressure of 1000 bar. For a shock strength below 100 bar both models (II and III) give identical results which indicate that in this region the internal energy contributions from dissociation and ionization are negligible. Hence, for a shock strength above 100 bar, neglecting the internal energy contribution from dissociation and ionization may lead to serious errors in structural design.

Validation of the numerical model (shown by crosses in figure 1) is done by comparing numerically calculated reflection coefficients with model III of the semi-analytical models. Numerical simulations are carried out for uniform shocks of various strength reflected from a fixed rigid wall to predict the reflection coefficient. Apart from demonstrating good agreement between the numerical and the approximate calculation, figure 1 essentially supports the value of *C*_{R} reported in literature [3,11]. Significant changes in the reflection coefficient (*C*_{R}) owing to the consideration of a real gas assumption for shock wave propagation in air during high-intensity explosion prompt a reassessment of the existing FSI theories presented in §2.

### (b) Constant shock reflection from a free-standing plate

Transmitted impulse (*j* varying from 0 to 6) is computed. The term *t*_{i} is the duration of the pressure pulse experienced by the plate since impingement. It is chosen to be 0.1 *ms* in this manuscript. Simulations are carried out for constant shock with peak over-pressure ratios (*p*_{s}/*p*_{0}) 5, 50, 250 and 500. In figure 2*a,* percentage changes in transmitted impulse, ((*I*_{preal}−*I*_{pideal})/*I*_{pideal}×100) owing to consideration of real gas EOS instead of ideal behaviour assumption [2] are plotted as a function of acoustic FSI parameter (

It can be observed from the plot (figure 2*a*) corresponding to *p*_{s}/*p*_{0} of 5, the percentage change in transmitted impulse is negligible and an ideal gas assumption is sufficient; but for higher *p*_{s}/*p*_{0} ratios (i.e. 50, 250 and 500), assumption of the real gas EOS is necessary, because the percentage change in between the two modelling assumptions is significant. A decrease in percentage change of normalized transmitted impulse can be observed for lightweight plates (around 22% for *m*_{p}=0.004 kg m^{−2} subjected to *p*_{s}/*p*_{0} of 500), whereas an increase could be observed for heavy plates (around 62% for *m*_{p}=400 kg m^{−2} subjected to *p*_{s}/*p*_{0} of 500). The absolute value of the percentage change for lightweight plates can be observed to be comparatively less than heavy plates because of the reduction in FSI timescale for lightweight plates (because for lightweight plates the reflected over-pressure drops to zero instantaneously after shock impingement). It should be noted that simulation results provided in figure 2*a* are with the CBP assumption (*p*_{b}=*p*_{0}) so as to demonstrate differences with the theory of Kambouchev *et al.* [2] which also considers CBP. However, it is reported [4,7,9] that receding motion of the plate induces significant pressure rise (resulting in shockwave) at the back of the plate (specially for lightweight plates) which influences the mitigation effect of FSI. In order to investigate this phenomenon, real gas behaviour on both the front and the back of the plate is considered and simulations are carried out to demonstrate differences in results between CBP and VBP cases for different *p*_{s}/*p*_{0} and plate masses.

It can be observed from figure 2*b* that significant reduction in relative transmitted impulse, i.e. ratio of transmitted impulse (*I*_{p}) to incident impulse (*I*_{i}=*p*_{s}*t*_{i}), can be obtained for lightweight plates in the VBP case (in comparison with CBP case) for higher *p*_{s}/*p*_{0}. Lightweight plates, upon shock impingement, gain significant momentum which causes a rise in back side pressure, thereby reducing the transmitted impulse. This phenomena reduces as the plate mass increases; in the heavy plate limit (*V* _{p}→0 along with pressure (*p*_{b}→*p*_{0}) at the back side of the plate, which explains why the relative impulse transmission for both CBP and VBP are the same in the heavy plate limit.

In figure 3*a,* normalized over-pressure (*p*_{s}/*p*_{0}=100 with a duration of 0.1 ms, impinging on a plate of mass 0.4 kg m^{−2}. Figure 1*a* clearly demonstrates that VBP condition suffers from higher pressure loading on both the front and the back of plate than CBP, i.e. [*p*_{f}(*t*)]_{VBP}>[*p*_{f}(*t*)]_{CBP} and [*p*_{b}(*t*)]_{VBP}>[*p*_{b}(*t*)]_{CBP}. However, because resistance provided by receding motion of plate is more than CBP, VBP suffers from less net total pressure, and, thereby, impulse transmission is less in the VBP condition in comparison with that of CBP (as shown in figure 2*b*). Reduction of net total pressure in the VBP case also causes reduction in velocity (*V* _{p}(*t*)) and displacement of the plate (*x*_{p}(*t*)) in comparison with that of CBP (as observed from figure 3*b*).

### (c) Performance chart and approximate formula for practical use

Approximate formulae to estimate transmitted impulse for an arbitrary plate mass and shock strength are prepared in such way that it should satisfy the acoustic solution for low-intensity shock and the asymptotic limits [2]. The formula of the KNR model (equation (2.13)) requires coefficients *C*_{R} (equation (2.18)), *β*^{CBP}_{s} (equation (2.14)) and f_{R} (equation (2.17)), which are based on the ideal gas assumption. However, the ideal gas assumption is applicable up to a shock strength of approximately 7 bar, as suggested by figure 1 and also by Baker [11]. Furthermore, the KNR model also neglects the pressure rise owing to receding motion of the plate and uses the CBP assumption. Thus, to address these deficiencies, the approximate formula is modified as
*C*^{real}_{R}, *β*^{real}_{s} and f^{real}_{R} are obtained using the real gas assumption. For different incident shock strengths, the *C*^{real}_{R} value can be obtained from model III of figure 1, whereas *β*^{real}_{s} and f^{real}_{R} can be obtained from the following expressions as

Incident shock velocity (*U*_{1}) can be expressed in a simplified form of linear shock-particle velocity relation as *U*_{1}=*c*_{0}+*s*_{1}*u*_{1}. For a given incident shock strength, particle velocity (*u*_{1}) and density of the shocked zone (*ρ*_{1}) can be obtained as

In the case of the VBP plate, it is assumed that the fluid particle adjacent to the back side of the plate has the same velocity as the plate, i.e. *u*_{b}=*u*_{p}, where *u*_{p} is the velocity of the plate and can be obtained from mass-momentum conservation across the plate (detailed derivation can be found in Ghoshal & Mitra [9]). This yields a cubic polynomial of plate velocity, given as
*u*_{p} is a real positive root of the polynomial. Shock wave velocity behind the plate (*U*_{b}) can be approximated by a linear shock particle velocity relationship as *U*_{b}=*c*_{0}+*s*_{1}*u*_{p} and the density of the shocked air at the back side of the plate can be obtained as
*β*^{real}_{s} is analogous to the FSI parameter of the LY model [4] in the acoustic range (*m*_{p}→0, *ρ*_{b}→*ρ*_{1} and *U*_{b}→*U*_{1} as if no plate is present. At very low shock intensity, the expression (equation (5.1)) reduces to the acoustic value, i.e. when *p*_{s}→0, then *C*_{R}→2, *U*_{1}→*c*_{0}, *I*_{p}/*I*_{i} in equation (5.1)

In figure 4, the empirical expression of transmitted impulse (equation (5.1)) for incident shock over-pressure ratio 5, 50, 250 and 500 is plotted along with their corresponding transmitted impulse obtained from numerical simulation using the VBP case. In can be clearly seen that this expression agrees well with numerically computed values not only in the heavy plate or light plate limit, but also in the range of intermediate plate mass. Thus, this expression reduces the complexities associated with numerical modelling and can be used to predict design transmitted impulse for given shock strength and plate mass.

## 6. Conclusion

This manuscript demonstrates that consideration of real gas modelling of medium could explain the long-standing discrepancy in design guidelines for high-intensity air-explosion-induced shock loading on structures. It is successfully shown that the internal energy contribution from vibration, ionization and dissociation accounts for the difference in reflection coefficient values reported in the literature and existing FSI theories.

It is identified that for a shock strength below approximately 7 bar the ideal gas assumption holds, whereas for shock strength in excess of 7 bar, the internal energy contribution from the degrees of freedom mentioned before should not be entirely neglected with consideration of an ideal gas. It is also identified that for a shock strength less than 100 bar, the contribution from ionization and dissociation to the internal energy is less compared with vibrational energy. However, for a peak over pressure beyond 100 bar the effects of ionization and dissociation play a significant role.

Analysis shows that the percentage change in transmitted-impulse owing to real gas modelling of air instead of ideal gas increases up to 62% in the case of heavy plates and decrease up to 22% for lightweight plates (for a maximum value of shock strength 500 bar considered during simulation). This change in transmitted-impulse for lightweight plates decreases further if the effects of shock wave formation at the back side of the plate is considered, which is relevant for high-intensity explosion-induced shock loads. Analysis also demonstrates advantages of lightweight structures in reducing relative impulse transmission. The results indicate that inaccuracy in obtaining *C*_{R} (using the ideal gas assumption for a peak over-pressure >7 bar) may lead to faulty design which may eventually lead to serious structural casualties. Empirical formulae to predict transmitted impulses are also developed to eliminate complexities of the computational model.

However, in this context, it should also be mentioned that this work can be extended to obtain better and more accurate results by following more accurate physics of shock wave propagation through air (considering the conductive and radiative heat flux and/or viscous flow [18]) at the cost of increased computational overhead and complexity.

## Data accessibility

The manuscript has no data.

## Author contributions

R.G. and N.M. both are responsible for establishing the theory within the paper.

## Funding statement

This work was carried out under the auspices of Naval Research Board, India under award no. NRB-226/HYD/10-11.

## Conflict of interests

We have no competing interests.

## Acknowledgements

Authors are thankful to the anonymous referees for their valuable suggestions.

## Appendix A. Parameters used in seven species real gas model

The values of several parameters used in seven species of real gas models of air, obtained from several articles [19,28,30], are provided in table 1 for ease of reading and reproducibility of the manuscript.

## Appendix B. Thermodynamic relations

Calculation of thermodynamic properties, i.e. *c*_{v}, *c*_{p}, *κ* and *χ*, are provided here.

Specific heat at constant volume (*c*_{v}), from its definition can written as
*c*_{v} can be calculated from the weighted contribution from each constituent expression being
*T*) at constant volume gives

Specific heat at constant pressure (*c*_{p}), by definition is
*c*_{p} can be written as

The expression of *κ* can obtained from the thermodynamic relation, given as

In order to obtain expression for *χ* the following thermodynamic relation is used

Differentiation of each component of the internal energy in the caloric EOS (equation (3.9)) with respect to temperature (*T*) at constant pressure and constant volume is as follows:

Derivatives of the mass fractions, used in the above equations, are as follows:

- Received October 24, 2014.
- Accepted January 27, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.