SimEx: A Tool for the Rapid Evaluation of the Effects of Explosions

Abstract

The dynamic response of structural elements subjected to blast loading is a problem of growing interest in the field of defense and security. In this work, a novel computational tool for the rapid evaluation of the effects of explosions, hereafter referred to as SimEx, is presented and discussed. The classical correlations for the reference chemical (1 kg of TNT) and nuclear (${10}^6$ kg of TNT) explosions, both spherical and hemispherical, are used together with the blast wave scaling laws and the International Standard Atmosphere (ISA) to compute the dynamic response of Single-Degree-of-Freedom (SDOF) systems subject to blast loading. The underlying simplifications in the analysis of the structural response follow the directives established by UFC 3-340-02 and the Protective Design Center Technical Reports of the US Army Corps of Engineers. This offers useful estimates with a low computational cost that enable in particular the computation of damage diagrams in the Charge Weight–Standoff distance (CW–S) space for the rapid screening of component (or building) damage levels. SimEx is a computer application based on Matlab and developed by the Fluid Mechanics Research Group at University Carlos III of Madrid (UC3M). It has been successfully used for both teaching and research purposes in the Degree in Security Engineering, taught to the future Guardia Civil officers at the Spanish University Center of the Civil Guard (CUGC). This dual use has allowed the development of the application well beyond its initial objective, testing on one hand the implemented capacities by undergraduate cadets with the end-user profile, and implementing new functionalities and utilities by Masters and PhD students. With this experience, the application has been continuously growing since its initial inception in 2014 both at a visual and a functional level, including new effects in the propagation of the blast waves, such as clearing and confinement, and incorporating new calculation assistants, such as those for the thermochemical analysis of explosive mixtures; crater formation; fragment mass distributions, ejection speeds and ballistic trajectories; and the statistical evaluation of damage to people due to overpressure, body projection, and fragment injuries.

1. Introduction

Unlike the slow energy release exhibited by deflagrations, the instantaneous energy deposition associated with the detonation of a high explosive produces an extremely rapid increase in temperature and pressure due to the sudden release of heat, light, and gases [1]. The gases produced by the explosion, initially at extremely high temperatures and pressures, expand abruptly against the surrounding atmosphere, vigorously pushing away any other object that may be found in their path. This gives rise to the two most notable effects of explosions: the aerial, or blast, wave [2], and the projection of shell fragments or other items (i.e., secondary fragments) located in the surroundings of the charge [3]. If the explosive device is located at a ground level, a fraction of its energy is effectively coupled to the ground, generating seismic waves and a well distinguished surface crater that results from the ejection of the shattered ground materials in direct contact with the charge [4]. Quantifying these phenomena and assessing their effect on the environment, including structural elements, vehicles, objects, or people located around the blast site, is a highly complex task that requires a thorough knowledge of the physical-chemistry of explosions [1,5,6,7] and their dynamic interactions with nearby structures [8] or the human body [9].

As a result of the growing terrorist threat experienced in the last few decades [10], estimating the effects of explosions has become a critical issue in the design, protection, and restoration of buildings and infrastructures, both civil and military [11]. However, this task is far from trivial, in that it involves transient compressible flows, nonlinear structural response, and highly dynamic fluid–structure interactions. These phenomena can be described with some accuracy using multiphysics computational tools, also known as hydrocodes [12], such as Ansys Autodyn, LS-Dyna, or Abaqus, based on the explicit finite element method [13]. In the simulations, all the critical components are modeled, including the detonation of the explosive charge, the resulting blast wave, the induced dynamic loads, and the nonlinear structural response. However, the enormous computational effort required to complete detailed computational analyses, which includes not only the calculation time itself, but also complex pre- and post-processing stages, remains a critical issue. For instance, simulating the effect of an explosive charge on a full-scale bridge may require more than 10 million finite elements [14]. For this reason, most engineering analyses still make use of simplified models for determining the explosive loads and estimating the resulting dynamic structural response in a timely manner. This enables the fast computation of damage diagrams in the Charge Weight–Standoff distance (CW–S) space, of utility to determine the level of protection provided by an input structural component loaded by blast from an input equivalent TNT charge weight and standoff [15].

In this regard, the American Unified Facilities Criteria UFC 3-340-02 [16], which supersedes the former ARMY TM 5-1300, establishes the requirements imposed by the US Department of Defense in the tasks of planning, design, construction, maintenance, restoration, and modernization of those facilities that must be protected against explosive threats. In the absence of similar regulations in other countries, UFC 3-340-02 [16] is widely used by engineers and contractors outside the US, as it provides a valuable guide for calculating the effects of blast-induced dynamic loads, including step-by-step procedures for the analysis and design of buildings to resist the effects of explosions.

To facilitate the application of the procedures set forth in the UFC 3-340-02 [16], as well as other analyses established in classic references of explosives engineering [3,5,6,9,17,18,19], fast evaluation software tools have been developed that incorporate the vast amount of data available as tables or graphs in the literature [7]. For instance, the United States Army Corps of Engineers (USACE) has developed and provides support for a series of software packages related to the design of explosion-resistant buildings [20]. Those tools were developed with public funding, and therefore there are regulations that restrict distributing those products outside of the United States. In addition, given the critical nature of this knowledge, access to these packages is severely limited to US government agencies and their contractors, with use only authorized to US citizens.

The inability to access these software packages motivated the authors to develop their own computational toolbox for the rapid evaluation of the effects of explosions. The result was the SimEx platform to be presented in this work. Conceived initially for educational purposes, the main goal was to develop a virtual software platform with an easy and intuitive Graphical User Interface (GUI) to be used in the computer lab sessions of the Explosion Dynamics course of the Degree in Security Engineering, taught at the University Center of the Civil Guard (CUGC) in Aranjuez, Spain. The Civil Guard is the oldest and biggest law enforcement agency in Spain. Of a military nature, its competencies include delinquency prevention, crime investigation, counter-terrorism operations, coastline and border security, dignitary and infrastructure protection, as well as traffic, environment or weapons and explosives control using the latest research techniques. The paradigm of the Civil Guard’s capacity is its outstanding role in the defeat of the terrorist group ETA, the longest-running terrorist group in Europe and the best technically prepared. In this context, the main target of the Degree in Security Engineering is the training of Guardia Civil cadets (i.e., the Guardia Civil’s future officer leadership) in the development, integration, and management of last generation civil security systems.

The purpose of SimEx was initially limited to the blast damage assessment on simple structural elements [21], such as beams, columns, pillars, or walls, following the Single-Degree-of-Freedom (SDOF) system analysis established by UFC 3-340-02 [16]. The tool has been successfully used since its initial inception in 2014 in both the computer lab sessions of the Explosion Dynamics course, and as a research tool for the development of a number of Bachelor and Master’s theses on explosion dynamics and blast effects. This double use as end-users and software developers by the Civil Guard cadets and students from other UC3M degrees has enabled the development of the application well beyond the initially planned objectives [22]. As a result, the current version of SimEx incorporates advanced topics in blast wave propagation, such as the prediction of cleared blast pressure loads due to the generation of rarefaction waves, as well as confined blast loading in vented structures [23]. It also includes several other calculation assistants for the thermochemical analysis of explosive mixtures [5,7,24]; crater formation [4,6,25]; fragment mass distributions, ejection speeds and ballistic trajectories [3,26,27,28]; and the statistical evaluation of damage to people due to overpressure, body projection and fragment injuries [9,29,30].

2. SimEx Capabilities

This section presents the current capabilities of SimEx, starting with the main interface used for computing the dynamic response of SDOF systems subjected to blast loading, and following with the description of the remaining calculation assistants.

2.1. Single-Degree-of-Freedom System Analysis

In many situations of practical interest, the response of structural elements to blast loading can be reduced, in first approximation, to that of an equivalent spring-mass SDOF system. As sketched in Figure 1, this system is made up of a concentrated mass subject to external forcing and a nonlinear weightless spring representing the resistance of the structure against deformation [8]. The mass of the equivalent system is based on the component mass, the dynamic load is imposed by the blast wave, and the spring stiffness and yield strain on the component structural stiffness and load capacity. Generally, a small viscous damping is also included to account for all energy dissipated during the dynamic response that is not accounted by the spring-mass system, such as slip and friction at joints and supports, material cracking, or concrete reinforcement bond slip [31].

If the system properties are properly defined, the deflection of the spring-mass system, $x\left(t\right)$, will reproduce the deflection of a characteristic point on the actual system (e.g., the maximum deflection). The system properties required for the determination of the maximum deflection are the effective mass of the equivalent SDOF system, $M_e$, the effective viscous damping, $C_e$, the effective resistance function, $R_e\left(x\right)$, and the effective load history acting on the system, $F_e\left(t\right)$. To systematize the calculations, the effective properties are obtained using dimensionless transformation factors that multiply the actual properties of the blast-loaded component, respectively, M, C, $R\left(x\right)$, and $F\left(t\right)$ [32]. These factors are obtained from energy conservation arguments in order to guarantee that the equivalent SDOF system has the same work, kinetic, and strain energies as the real component for the same deflection when it responds in a given, assumed mode shape, typically the fundamental vibrational mode of the system [31].

In the analysis of blast-loaded SDOF systems, it is therefore of prime importance to identify the fundamental vibrational mode of the structural element. This procedure is not trivial, since obtaining the fundamental mode can entail certain difficulties, in which case its shape must be approximated in some way [32]. To determine the equivalent properties of the SDOF system, it is also necessary to determine the type of structure (beam, pillar, frame, etc.) and how the load is applied (typically, a uniform load is assumed). The elastic behavior of the material is often modeled as perfect elasto-plastic, probably the simplest of all nonlinear material models. This assumes that the initial response follows a linear elastic behavior described by an apparent elastic constant K, but once the yield strain is reached, $x\ge x_{\mathrm{u}}$, the material behaves as plastic, flowing at a constant stress with an ultimate resistance $R_{\mathrm{u}}=K{x}_{\mathrm{u}}$, i.e.,

(1)$R\left(x\right)=\left\{\begin{array}{c}\hfill Kx \mathrm{for} \left|x\right|<x_{\mathrm{u}}\\ \hfill R_{\mathrm{u}} \mathrm{for} \left|x\right|\ge x_{\mathrm{u}}\end{array}\right.$

Although more complex models could be used, they are not considered here due to the heavy simplifications introduced in the formulation of the problem.

The mass transformation factor, $K_M$, is defined as the ratio between the equivalent mass $M_e$ and the real mass M of the blast-loaded component; the load transformation factor $K_L$ is defined as the ratio between the equivalent load $F_e\left(t\right)$ and the actual load $F\left(t\right)$, and usually coincides with the resistance and damping transformation factors; and finally the load-mass factor $K_{LM}$ is defined as the ratio between the mass factor and the load factor

(2)$K_M=\frac{M_e}{M} ;K_L=\frac{F_e\left(t\right)}{F\left(t\right)}=\frac{R_e\left(x\right)}{R\left(x\right)}=\frac{C_e}{C} ; K_{LM}=\frac{K_M}{K_L}=\frac{M_e}{M}·\frac{F\left(t\right)}{F_e\left(t\right)}$

Although all these factors are easy to obtain, even through analytical expressions in some cases, most of them can be found tabulated in the UFC-3-340-02 [16].

The linear momentum equation for the equivalent SDOF system then takes the form [32]

(3)$K_{LM}M\ddot{x}+C\dot{x}+R\left(x\right)=F\left(t\right)$

where, as previously discussed, C represents the viscous damping constant of the blast-loaded component. This constant is often specified as a small percentage, z, of the critical viscous damping, $C=(z/100)C_{\mathrm{cr}}$, with a damping coefficient $z=2$ being a good value when not otherwise known (for further details see [31]). Note, however, that damping has very little effect on the maximum displacement, which typically occurs during the first cycle of oscillation, so the actual value of z is not of major relevance. The inhomogeneous term, $F\left(t\right)$, appearing on the right-hand side of Equation (3) represents the dynamic load associated with the blast wave, to be discussed in Section 2.1.1 below.

SimEx provides an easy and intuitive GUI environment for the study of the dynamic response to blast loadings of a variety of structural elements that can be modeled as SDOF systems. Figure 2 shows the main SimEx interface, divided into three calculation assistants for the three basic elements that make up the SDOF system: a module for calculating the properties of the blast wave (forcing term, $F\left(t\right)$), a module for calculating the equivalent mechanical properties (resistance term, $R\left(t\right)$), and a module for the numerical integration of the problem, which includes the post-processing of the results and their graphic representation in the form of displacements, forces, and deformation diagrams (see the bottom plots of Figure 2) and of CW–S damage charts, to be discussed in Section 3.3.

As a final remark, it is important to note that, following standard practice, the SDOF analysis carried out by SimEx uses the load defined in terms of pressure, $F\left(t\right)=p^{\prime }\left(t\right)$ ($\mathrm{Pa}$), so that both the mass M ($\mathrm{kg}/{\mathrm{m}}^2$), the damping coefficient C ($\mathrm{kg}/({\mathrm{m}}^2 \mathrm{s})$) and the ultimate resistance $R_{\mathrm{u}}$ ($\mathrm{Pa}$) must all be introduced as distributed values per unit surface (p.u.s.) in the different calculation assistants.

2.1.1. Forcing Term

As previously discussed, the blast wave overpressure defined in Equation (4) below can be used directly in Equation (3) as forcing term, $F\left(t\right)=p^{\prime }\left(t\right)$, as long as the analysis is formulated per unit surface and uses distributed masses and forces. In order to determine the blast parameters (arrival time, peak overpressure, positive phase duration, impulse per unit area, waveform parameter, etc.), classical correlations [1,2,17,18,19,33,34] in terms of scaled distance are used together with the scaling laws for spherical or hemispherical blast waves [1,17,35,36], which allow their evaluation for arbitrary CW–S pairs. It is interesting to note that the standoff distance is defined as the minimum distance from the charge to the structural element under study (e.g., a wall). However, the actual distance to a given point of that element, e.g., the centroid (or geometric center), which may be considered the most representative point of the structure, may be slightly different due to the incidence angle being larger than 0 at that point.

The local atmospheric pressure, $p_a$, and temperature, $T_a$, are determined using the International Civil Aviation Organization (ICAO) Standard Atmosphere (ISA) [37] with a temperature offset (ISA$\pm \Delta T$). The user must specify the geopotential height, in meters, and the non-standard offset temperature $\pm \Delta T$, although arbitrary ambient temperature and pressure can also be introduced directly [38]. TNT is used as reference explosive, although the results can be extrapolated to other compositions using either the equivalence tables included in SimEx for selected explosives [39], or the thermochemical calculation assistant, to be presented in Section 2.2.1, for less conventional formulations or explosive mixtures.

To estimate the dynamic load exerted by the blast wave, the angle of incidence of the incoming shock wave must be considered, the worst-case conditions being usually those of normal incidence. UFC 3-340-02 [16] contains scaled magnitude data for both spherical and hemispherical blast waves. It also provides methods to calculate the properties of the blast wave with different incidence angles, including both ordinary and Mach reflections for oblique shocks. The time evolution of the blast wave overpressure $p^{\prime }\left(t^{\prime }\right)$ at a fixed distance, d, sufficiently far from the charge (at least, larger than the fireball scaled distance) is approximated using the modified Friedlander’s equation, which captures also the negative overpressure phase [1,17,40]

(4)$p^{\prime }\left(t^{\prime }\right)=p\left(t^{\prime }\right)-p_1=p^{\circ }\left(1-\frac{t^{\prime }}{t_d}\right)\exp \left(-\alpha \frac{t^{\prime }}{t_d}\right)$

where $p^{\circ }=p_2-p_1$ represents the peak overpressure measured from the undisturbed atmospheric pressure $p_1=p_a$, with $p_2$ denoting the peak post-shock pressure, $t^{\prime }=t-t_a$ is time measured from the blast arrival time, $t_d$ is the positive phase duration, and $\alpha$ is the waveform parameter, closely related to the impulse per unit area of the positive phase $I/A={\int }_0^{t_d}{p}^{\prime }\left(t^{\prime }\right) \mathrm{d}{t}^{\prime }$ (area under the positive phase of the overpressure-time curve) according to $I/A=p^{\circ } t_d\left[1/\alpha -(1-e^{-\alpha })/{\alpha }^2\right]$. SimEx performs by default the complete integration of the Friedlander waveform, but the equivalent triangular pressure pulse can also be used without significant errors [32]. This simplified waveform has the same maximum peak overpressure, $p^{\circ }$, but a fictitious positive duration computed in terms of the total positive impulse and the peak over pressure, $t_d=2(I/A)/p^{\circ }$.

The “Blast wave” calculation assistant allows the activation of the effects of clearing and confined explosions, which increases the computational capabilities to more realistic situations. The clearing effect takes into account the time required for reflected pressures to clear a solid wall that has received the impact of a blast wave as a result of the propagation of rarefaction waves from the edges of the wall. In the case of confined explosions, SimEx implements the procedure outlined in UFC 3-340-02 [16] to estimate the gas phase peak overpressure and duration of the equivalent triangular pressure pulse in terms of the chamber’s total vent area and free volume. These effects can be activated on the lower part of the “Blast wave” calculation assistant.

2.1.2. Resistance Term

The “Resistance” calculation assistant provides a means to define the equivalent mechanical properties (i.e., structural mass, damping coefficient, and structural strength) of the SDOF system under study modeled as a perfectly elasto-plastic system with elastic stiffness K until the yield strain, as given in Equation (1). The characteristic length, L, of the structural element must also be provided, as it is required to determine the maximum rotation angle at its boundaries, often referred to as support rotation, $\theta$. For the equivalent SDOF system, the assistant computes the fundamental natural period, $T_n=2\pi \sqrt{K_{\mathrm{LM}}M/K}$, the critical damping, $C_{\mathrm{cr}}=2\sqrt{K_{\mathrm{LM}}KM}$, and the deflection at which plastic deformation initiates in the system, $x_{\mathrm{u}}$. Direct access to calculation assistants that compute the equivalent properties (M, K, $K_{\mathrm{LM}}$, $R_{\mathrm{u}}$) required for the calculations is also provided for various types of systems. Currently, standard European wide flange “metal beams” [41] and reinforced “concrete beams” are included (see Section 3.2), although it could be possible to incorporate additional assistants for other elements, such as metal panels/plates, open-web steel joists, reinforced concrete slabs, reinforced/unreinforced masonry, or wood panels/beams. The metal beams assistant also provides the possibility of studying custom (i.e., non-normalized) profiles and materials in order to widen the computation capabilities.

2.1.3. Numerical Integration

Once the characteristics of the equivalent SDOF system have been defined, the resulting ordinary differential equation that models the transient nonlinear response of the equivalent structural system (3) must be integrated numerically. The integration module implements the two numerical methods recommended by UFC-3-340-02 [16], namely the “Acceleration-Impulse-Extrapolation Method” and the “Average Acceleration Method” [16], which can be selected from a drop-down menu. Text boxes are also included to set the initial conditions (displacement and initial speed, which are zero by default) as well as the final integration time. Since both numerical methods use constant time steps, a sufficiently short time increment, typically of the order of a few percentage of the natural period or the positive phase duration (usually, fractions of a millisecond), should be used in order to ensure the numerical convergence of the integration.

2.1.4. Post-Processing

After integration, three plots appear in a pop-up window and a summary table is provided at the bottom left corner of the main window. The left plot shows the instantaneous displacement (solid line) and the permanent displacement, or deformation (dashed line). The central plot shows the temporal variation of the forcing term (i.e., the blast pressure wave, solid line) together with the resistance strength of the SDOF system (dashed line). The right plot shows the displacement–resistance graph, in which it is possible to determine more clearly whether permanent deformations occur or not. Finally, the table of results shows the maximum displacement obtained, $x_{\max }$, along with two damage indicators: the ductility ratio, $\mu =x_{\max }/x_{\mathrm{u}}$, defined as the ratio of the peak deflection to the ultimate elastic deflection, and the maximum support rotation, $\theta$, whose calculation depends on the type of structure under study.

By integrating different combinations of charge weights and standoff distances for the same structural element, damage level diagrams can be rapidly obtained in the CW–S distance space. SimEx has a function for it located in the central part of the integrator module. One can select the range of charge weights and standoff distances, the number of intermediate values and the type of damage in terms of the quantitative indicators $\mu$ and $\theta$ [15]. From the two quantitative indicators, the structural damage level can be classified qualitatively into: superficial, moderate, heavy, hazardous failure, and blowout, with response limit boundaries between these levels denoted respectively by B1 (superficial to moderate), B2 (moderate to heavy), B3 (heavy to hazardous failure), and B4 (hazardous failure to blowout). Convenient limits for the boundaries of component damage levels for common structural components in terms of $\mu$ and $\theta$ are provided in [15]. An example of a damage level diagram for the façade of a conventional building subject to blast loading computed with SimEx will be presented in Section 3.3.

2.2. Other Calculation Assistants

The main SimEx interface gives access to several other calculation assistants. These include: a module for the calculation of the theoretical (i.e., thermochemical) properties of explosives and explosive mixtures; a module for estimating the initial velocity, mass distribution and ballistic trajectories of primary fragments; a crater formation calculator; and a module for estimating damage to people, including both primary and tertiary injuries. The fragment assistant also provides estimations of the secondary injuries due to the impact of primary fragments on people. In this section, we shall briefly present and discuss the above-mentioned assistants.

2.2.1. Assistant for the Calculation of the Thermodynamic Properties of Explosives

For the calculation of the theoretical thermodynamic properties of explosives and explosive mixtures, SimEx includes an extensive database of pure CHNO propellants and explosives extracted from Kinney and Graham [1], updated with data from Meyer [7] and Akhavan [5] for more recent explosives. From the properties of pure explosives, the thermochemical assistant estimates the properties of explosive mixtures formed by two or more components by specifying the mass fractions and the density of the mixture.

First, it computes the apparent chemical formula of the explosive mixture along with its molecular weight and maximum density. For the calculation of the decomposition reaction in nominal products, which provides the heat of explosion and the volume of gases generated, one can choose different calculation hypotheses: Kamlet–Jacobs (KJ), Kistiakowsky–Wilson (KW), Modified Kistiakowsky–Wilson (modified KW), Springall–Roberts (SR), or chemical equilibrium [5]. In the latter case, SimEx determines the composition of the product mixture following the chemical equilibrium approach considering a constant–volume explosion transformation that uses the ideal gas Equation of State (EoS) for the products according to the norm UNE 31-002-94 [42], as illustrated in Figure 3.

More complex computations based on the European Standard EN 13631-15 [43], which use the semi-empirical Becker–Kistiakowsky–Wilson (BKW) EoS [44,45,46] or the Heuzé (H9) EoS [47] for the products, are also supported in the last version of SimEx. As sample results of these computations, Table 1 shows the detonation properties obtained by SimEx for different explosive mixtures (see Table 2 for its composition) compared with the results reported in the European Standard EN 13631-15 [43], and obtained with the W-DETCOM code [48,49], which computes directly the Chapman–Jouguet state.

The equilibrium calculations are carried out using Combustion Toolbox (CT), an in-house thermochemical equilibrium package developed at UC3M [24,50]. CT determines the equilibrium composition of the product mixture through the Gibbs free energy minimization method by using Lagrange multipliers combined with a multidimensional Newton–Raphson method. The thermodynamic properties (specific heat, enthalpy, and entropy) are computed as a function of temperature derived from NASA’s 9-coefficient polynomial fits for combustion of ideal and non-ideal gases and condensed phases.

From the resulting composition of the product mixture at equilibrium, the assistant computes the volume of gases generated, the heat of explosion, the Gurney constant, the detonation pressure, the detonation velocity, and the explosive force (or power index). To estimate the detonation pressure and velocity, the approximate expressions of Kamlet & Jacobs [51,52] are used, whereas the explosive force is estimated using the well-known Berthelot approximation [1]. These data are subsequently used to calculate the TNT equivalent of the explosive composition under study.

2.2.2. Crater

SimEx also has an assistant for the direct and inverse calculation of craters based on the classical correlations for craters reviewed by Cooper [6] (see also Refs. [4,25]), whose interface is shown in Figure 4. With this assistant, one can calculate the radius of the crater generated by the detonation of a certain amount of a given explosive at a certain height above the ground, considering different types of soil. It is also possible to calculate the explosive charge required to produce a crater of a certain size, which may be useful for the forensic analysis of explosions [53]. Buried craters are not yet included in the assistant, but could be incorporated in future versions following the work of Westine [54], as reviewed by Baker et al. [55].

2.2.3. Primary Fragments

SimEx incorporates assistants for calculating the mass distribution, ejection velocity, and ballistic trajectory of primary fragments. The corresponding interfaces are shown in Figure 5, Figure 6 and Figure 7. The fragment size distribution is estimated using Mott’s statistical theory for fragmentation of steel cylindrical shells [3,26,27,28], as suggested by UFC-3-340-03 [16]. As shown in Figure 5, this model determines the average number of fragments and their average weight. It also provides the size of the largest fragment corresponding to a given Confidence Level (CL). SimEx also includes a ballistic trajectory assistant for primary fragments that, in addition to the flight path, provides the flight time, velocity, and maximum distance, as illustrated in Figure 6. The initial velocity of primary fragments is computed using Gurney’s analysis [56] for cylindrical, spherical, and symmetrical/asymmetrical sandwich charges. Although this analysis assumes that all fragments have the same the initial velocity, given the different fragment sizes, both their initial kinetic energy and their subsequent aerodynamic deceleration are different. The assistant thus includes an initial aerodynamic deceleration chart, shown in Figure 7, that provides the fraction of the initial velocity achieved at a certain distance, given the fragment mass and material, and the local air density, specified through the ISA ±$\Delta T$ model. The aerodynamic assistants assume spherical fragments with a variable drag coefficient for all Mach numbers [57], although the model could be extended to account for more realistic (i.e., irregular) fragment shapes in future versions [58]. The results of these models are also used to estimate the lethality risk by the impact of primary fragment in the event of a strike on a person, which is found to depend on the speed and the mass of the fragment, as illustrated by Figure 5.

Non-tabulated explosives or explosive mixtures can also be considered, with the Gurney constant being computed by the thermochemical assistant presented in Section 2.2.1. In this case, the user must select a “custom” explosive, and the thermochemical assistant will open to specify the desired explosive composition. Once the wizard is closed, the Gurney constant is automatically exported to the fragment wizard.

2.2.4. Damage to People

SimEx includes an assistant for estimating damage to people using the widely accepted probit (probability unit) functions [58,59] provided by the TNO’s Green Book [9] and summarized in Table 3. For each type of injury or cause of death (eardrum rupture, lung injury, etc.), a probit function is defined that depends on the blast parameters: side-on, dynamic or reflected peak overpressure (depending on the body position), impulse per unit area, etc. For primary injuries, lethality due to lung damage is evaluated together with the probability of eardrum rupture. For tertiary injuries, lethality is evaluated for shock-induced body displacement and subsequent direct impact, either with the head or the whole body [29].

The appearance of the interface is shown in Figure 8. All necessary parameters can be selected on the left: size, type, and geometry of the explosive charge, as well as the body position relative to the incoming pressure wave, which determines whether side-on, dynamic, or reflected pressure is used to compute the peak overpressure and impulse. The rest of the window presents the results both numerically and graphically, using overpressure–impulse diagrams on the left and CW–S diagrams on the right, with primary injuries shown above and tertiary injuries below. Overpressure–impulse diagrams display the characteristic overpressure–impulse–distance curve for the selected charge weight to facilitate the interpretation of results [60], while CW–S diagrams include a diagonal dashed line indicating the approximated position of the fireball radius, corresponding roughly to an scaled distance $Z=d/W^{1/3}=1$ m. Above this line, the Freidlander waveform is not valid, and the blast wave parameters are increasingly imprecise [1].

3. Example of Application: Façade of a Building under Blast Loading

To illustrate the capabilities of SimEx, this section presents a preliminary study to asses the ability of a conventional three-story steel frame building, such as the one shown in Figure 9, to resist three different combinations of charge weight, W, and standoff distance, d, preserving a similar scaled distance, $Z=d/W^{1/3}$. The three CW–S combinations are summarized in Table 4. For simplicity, we assume mean sea level ISA conditions for all the calculations. For illustrative purposes, the figures quoted below show results corresponding to the first floor of the building (hereafter referred to as Level 1) and Case 2 conditions. That is, we shall consider as reference conditions a ground explosion of 150 kg of TNT at a 20 m standoff distance from the front façade of the building, as depicted in Figure 9a.

3.1. Incident Load

As previously discussed, SimEx allows the user to enter directly the desired CW–S combination to define the incident blast load. Figure 2 shows the results corresponding to the reference conditions (Level 1, Case 2). For a more detailed analysis of the load induced by the blast wave, the “Blast Wave” calculation assistant shown in Figure 10 allows a fast evaluation of all blast parameters as a function of the standoff distance. To this end, the user must provide the following input data: the ground distance from the explosion to the point of calculation, d, the elevation of the explosive charge, $h_c$, and the elevation of the calculation point, $h_0$, both measured from the ground.

For $h_c=0$, a hemispherical surface burst computed from Kingery and Bulmash parameters for TNT [18] is considered, although other correlations for hemispherical explosions [1] could also be selected. For $h_c>0$, hemispherical or spherical blasts are both available, letting the user decide what is the best option based on the height of burst. The code does not include correlations for more complex configurations, such as air bursts producing regions of regular and Mach reflections that eventually modify the incident shock wave. The user must also introduce the angle formed by the normal to the structural element at the point of calculation with the horizontal projection of the line joining the center of the explosion with that point, $\delta$, which is identically zero in our case studies if we assume a symmetric configuration with a pillar in the center of the front façade. These distances and angles are employed for simplicity in obtaining in-field measurements.

With these data, the wizard is able to compute the real distance and incidence angle, thereby providing the peak overpressure, $p^{\circ }$, the impulse per unit area, $I/A$, the duration of the positive phase, $t_d$, the blast arrival time, $t_a$, the average speed of the pressure front, $\sigma =d_{\mathrm{real}}/t_a$, the positive phase length, $L_w$, and the waveform parameter, $\alpha$. The results are presented in a table for several standoff distances, d, which also gives the real distances, $d_{\mathrm{real}}$, and angles of incidence, ${\delta }_{\mathrm{real}}$. The maximum and minimum distances that appear in the table can be easily modified by the user, who can select any intermediate value using a slider bar to compute the blast parameters at a fixed specified distance. A button has also been included to graphically represent the variation of any of the blast parameters as a function of the distance to the center of the explosion. The results are also exportable as a "comma-separated-value" format for further postprocessing.

For more qualitative information, two exportable graphs are presented in the lower part. The graph on the left displays the time evolution of the overpressure at a fixed horizontal distance. The user can change this distance easily with the slider bar. All the characteristics of the blast wave are shown for the particular distance chosen by the user. The graph on the right represents the maximum overpressure and the impulse per unit area as a function of the horizontal distance. As previously indicated, the range of distances is also adjustable by the user. Using the "Export and exit" button, the module is closed and the weight and type of explosive, the distance to the charge, and the real angle of incidence to be used in the integration of the SDOF system are exported to the main SimEx module. Figure 10 shows the calculation of the blast parameters for an explosive charge of 150 kg of TNT on a point at a height of 5.5 m above the horizontal, i.e., the geometric center of the façade of the first floor, corresponding to the reference case (Level 1, Case 2). Other distances are also included in the top table, showing how the angle of incidence tends to become normal as the charge moves away from the target.

3.2. Estimation of the Equivalent SDOF System Response

To study the structural response to an explosive charge, it is necessary to know in detail the type of construction. However, when using a simplified SDOF model, the study can be simplified and generalized for many different cases. In the present example, we will analyze a façade structure like the one in Figure 9b, composed of equally spaced pillars and an outer enclosure wall.

The first element that receives the blast wave is the enclosure of the façade. This, in turn, transmits the load to the rest of the structure. In most constructions, the façade is only an enclosure without structural function (glass façades, brick, etc.). In first approximation, it can be considered that the exterior enclosure transmits the full load received directly to the pillars. The pillars are structural elements whose integrity is considered critical. It will therefore be the first element to be studied since the protection of the supporting structure is pivotal to avoid the potential collapse of the building. The enclosure can be considered as a secondary element in most constructions and therefore a significantly higher level of damage than in primary elements can be allowed.

Figure 9c shows the simplest element in which the façade is to be divided. Each pillar receives loads from a part of the façade corresponding to the distance between pillars and the height between floors. The load generated by the explosion is applied to the pillars crosswise, so they behave in first approximation as bending elements. For the calculation of the equivalent properties, the beam assistants available in SimEx are employed. Either for metal or concrete beams, the length corresponds to the height between floors, while the span is the spacing between pillars. In the case of pillars, the boundary condition between floors is that of embedment on both sides, whereas a free condition is preferred at the roof. As a result, we use fixed-fixed conditions for Levels 0 and 1 and cantilever (or fixed-free) for Level 2. The presence of a roof diaphragm element may require additional considerations regarding the boundary condition at the roof top, but we prefer to use a fixed-free boundary condition for the second floor both for simplicity and for illustrating the effect of considering different boundary conditions on different floors.

In the case of metal beams, it is only necessary to indicate the standard shape of the profile and the size. SimEx uses European cross-section profiles HEB, IPE, and IPN in accordance with Euronorm 53–62 (DIN 1025) [41]. Figure 11 shows the result for a HEB 340 profile with a length of 3 m and a separation between pillars of 5 m. The assistant uses standardized profiles, so if a non-existent measure is introduced, it corrects down to the nearest lower normalized profile. However, it is also possible to select custom profiles and materials. In this case, the area, first moment of area about the bending axis, moment of inertia about the bending axis, density, Young’s modulus, and resistance must be provided by the user. Once the structural properties have been introduced, closing the assistant incorporates the computed data into the main SimEx interface. Figure 2 shows the result for the case under study. It should be noted that the additional enclosure mass supported by the pillar when flexed must also be included in the mass of the equivalent SDOF system in the main interface.

If a rectangular reinforced concrete pillar is considered, SimEx requires that the external measurements b and h (perpendicular and parallel to the direction of application of the load, respectively) be introduced. In addition, the properties of the reinforcement should be indicated in a simplified manner, that is, interior spacing, $d_c$, and reinforcement area, $A_s=n\pi d_{\mathrm{bar}}^2/4$, where n represents the number of steel reinforced bars per side. Figure 12 shows results for a pillar of 45 × 45 cm${}^2$ with 5 A36 steel reinforcement bars of # 7 (approximately 22.5 mm in diameter) per side, for a length of 3 m and a spacing between pillars of 5 m. The distance $d_c$ must be estimated according to the constructive detail. In this particular case, it is assumed that the reinforcement centers are located at 4 cm from the edge, resulting in an interior reinforcement spacing of $d_c=37$ cm.

It is worth noting that neglecting axial load can be considered a conservative approach, particularly in the case of columns or pillars. These elements are initially subjected to a significant compression load due to the weight of the supported structure, which reduces the tensile stresses caused by bending. This simplification constitutes a first approximation in the study of the structural response. For a more detailed analysis, the wall should be the next element to be analyzed in order to assure that it is able to fully transmit the blast load to the load-bearing element. If the wall was made of concrete, this could be done using the concrete beam assistant with $b=S$. In this case, the mass of the element under study would be the total mass of the equivalent SDOF system. However, in the case considered here of load-bearing elements (beams or columns/pillars), the total mass can be significantly larger than the mass of the element.

3.3. SDOF System Integration and CW–S Damage Diagrams

Once the user sets the explosive charge and the properties of the equivalent SDOF system, SimEx is ready to integrate the resulting mathematical problem. Figure 2 shows the results for the case of a HEB 340 profile with a 5 m span between pillars. The main results are the maximum deflection, $x_{\max }$, the ductility ratio, $\mu$, and the maximum rotation angle, $\theta$. The two latter parameters are used as indicators to quantify the component damage levels [15]. Assuming that the Level of Protection (LOP) required is very low, in case of a hot rolled compact steel shape for the columns, according to [15], the allowable component damage is heavy (response between B2–B3).

For fixed values of the structural parameters, a parametric sweep can be carried out in CW–S space to obtain damage diagrams such as the ones shown in Figure 13. To this end, it is enough to indicate in the assistant the charge weight and standoff distance ranges to be analyzed and the number of intervals to be used for each parameter. In addition, the desired damage level criteria must be indicated to separate the zones. Figure 2 shows characteristic values of $\mu$ and $\theta$ for metallic elements, although other values could be selected from [15] for other structural elements and materials. Note that CW–S damage diagrams are presented both in linear and log-log scales.

As can be seen, the CW–S damage diagrams shown in Figure 13 include three points corresponding to the three cases considered in Table 4. As the three scaled distances are almost equal, then the damage levels are also very similar, although differences in real distances and incidence angles make them grow from superficial-moderate (B1) to (almost) moderate-heavy (B2) for increasing charge weights and standoff distances. According to the PDC-TR 06-08 [15], a superficial damage level implies “no visible permanent damage”, whereas a moderate damage level implies “some permanent deflection” that generally can be repaired. By way of contrast, a heavy damage is associated with “significant permanent deflections” that cause the component to be unrepairable.

To summarize the results obtained in the different case studies, Table 5 reports the incident blast load parameters and the corresponding component damage indicators per floor for Cases 1, 2, and 3. The reference case (Level 1, Case 2) and the worst-case scenario (Level 2, Case 3) are both highlighted for clarity. As can be seen, damage levels are significantly higher in the upper floor (Level 2) as a result of the lowest rigidity imposed by the cantilever boundary condition at the roof top, resulting in heavy damage levels for cases 2 and 3.

3.4. Crater, Fragments, and Damage to People

Figure 4 presents an estimation of the crater generated in the reference case on a sandstone soil, with an approximated radius of $1.6$ m. For surface bursts, $\mathrm{HOB}=0$ m, as the one considered here, the equivalent charge radius is irrelevant, as it is only used to determine the dimensionless height of burst, which is identically zero in our example. The figure also shows that, for above-surface bursts, $\mathrm{HOB}>0$ m, the crater radius is significantly smaller for the same amount of explosive due to the air cushion that exists between the load and the ground, which reduces to a great extent the pressure that reaches the ground surface [6].

Figure 5 shows the interface of the fragment assistant using the input data of the reference case. For the application of Mott’s statistical theory for fragmentation of steel cylindrical shells [3,26,27,28], the explosive charge is approximated to a cylinder of approximately 50 cm diameter surrounded by a steel fragmentation shell with a mass of the order of about 13% of the charge and a thickness of 2 mm.

Finally, Figure 8 shows the calculating assistant for estimating damage to people in the reference case. As an illustrative example, the figure presents the results of lethality due to different types of injuries at a distance of 20 m from the origin of the explosion, assuming the worst-case scenario of an average person located close to the façade of the building being attacked. In the pressure-impulse graphs, representative distances are indicated using red dots plotted along the characteristic overpressure–impulse–distance curve [60]. As can be seen, at 20 m standoff distance, lethality due to lung damage or whole-body projection is negligible, but large primary fragments (e.g., CL 99%) may still produce secondary injuries with fatal results, as indicated by Figure 5.

4. Conclusions

SimEx is a computational tool that allows a rapid and easy estimation of the effects of explosions on structural elements and their damage to people. It has been developed in accordance with the specifications of American standard UFC-3-340-02 and other widely accepted directives published in the open literature. It provides assistants for the calculation of the blast-wave load; SDOF dynamic response, including the calculation of the equivalent structural properties of standardized metal and reinforced concrete beams; thermodynamic properties of explosive mixtures; crater formation; projection of primary fragments; and damage to people.

After presenting the main calculating assistants, a preliminary study has been presented to illustrate the full capabilities of SimEx in the assessment of the ability of a building to resist a given explosive charge. The analysis enables the determination of component damage levels for the main structural components, and a further study of the reference case has led to the computation of CW–S damage diagrams for a pillar of the first floor. These diagrams are very useful to provide design guidelines for those facilities that must be protected against explosive threats.

Although still under development, SimEx is being successfully used for research and teaching activities at the Spanish University Center of the Civil Guard. Due to its advanced stage of maturation, it could also be used in other areas within the Army and Law enforcement Agencies involved in the fight against terrorism and the design of blast resistant buildings and structures.

Author Contributions

Conceptualization, J.S.-M. and M.V.; methodology, J.S.-M. and M.V.; software, J.S.-M. and A.C.; validation, all authors; formal analysis, J.S.-M. and M.V.; investigation, all authors; resources, C.H. and M.V.; data curation, J.S.-M. and A.C.; writing—original draft preparation, J.S.-M., A.C. and M.V.; writing—review and editing, C.H. and M.V.; visualization, J.S.-M. and A.C.; supervision, C.H. and M.V.; project administration, M.V.; funding acquisition, C.H. and M.V. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowldgments

The authors would like to thank Henar Miguelez from University Carlos III of Madrid and Col. Fernando Moure from the University Center of the Civil Guard for their continuous support, as well as Com. Miguel Ángel Albeniz from the SEDEX-NRBQ (EOD-CBRN) service of the Civil Guard for many enlightening discussions. Fruitful discussions with Lina López from the School of Mines at Technical University of Madrid are also gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

SimEx License & Distribution

SimEx is a closed-source proprietary software that may be licensed by the copyright holders, UC3M & Guardia Civil, under specific conditions. Please contact the corresponding author for further information.

Abbreviations

The following abbreviations are used in this manuscript:

BKW	Becker–Kistiakowsky–Wilson EoS
CL	Confidence Level
CT	Combustion Toolbox
CUGC	Centro Universitario de la Guardia Civil
CW–S	Charge Weight–Standoff
EoS	Equation of State
GUI	Graphical User Interface
H9	Heuzé EoS
HOB	Height of Burst
IED	Improvised Explosive Device
ISA	International Standard Atmosphere
LOP	Level of Protection
PDC	Protective Design Center
SDOF	Single Degree of Freedom
SEDEX-NRBQ	Explosive Ordnance Disposal (EOD) and CBRN Defense Service
UC3M	University Carlos III of Madrid
UFC	Unified Facilities Criteria
US	United States
USACE	United States Army Corps of Engineers

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Sketch of the equivalent SDOF system showing the different terms involved in its mathematical description. Left: forcing term; right: resistance term; center: equivalent spring-mass SDOF system and its associated differential equation.

Figure 2. Main interface of SimEx showing the “Blast wave”, “Resistance”, and “Integration” assistants for the computation of the structural response of perfect elasto-plastic SDOF systems under blast loading. The access buttons to the other calculation assistants are seen under the top toolbar. The bottom plots show the post-processing pop-up window that displays the results of the numerical integration in terms of displacements, forces, and deformation diagrams (for a detailed discussion of these diagrams, see Section 2.1.4).

Figure 3. Interface of the assistant for the calculation of the theoretical thermodynamic properties of explosives and explosive mixtures.

Figure 4. Interface of the assistant for the calculation of craters. HOB denotes the height of burst.

Figure 5. Interface of the primary fragment mass distribution and lethality assistant.

Figure 6. Interface of the primary fragment calculation assistant showing the ballistic fragment trajectory, flight time, velocity, and maximum distance charts. Fragments are assumed spherical.

Figure 7. Interface of the primary fragment calculation assistant showing the initial deceleration chart, which provides the fraction of the initial velocity, [Forumla omitted. See PDF.], achieved at a certain distance (contour lines), given the fragment mass, [Forumla omitted. See PDF.], and material (e.g., lead), and the atmospheric conditions (e.g., ISA mean sea level). Fragments are assumed spherical.

Figure 8. Interface of the assistant for estimating blast-induced damage to people. The CW–S and atmospheric data, along withe the body position relative to the incoming pressure wave, are introduced in the top-left corner, the blast wave parameters and the statistical damage indicators for the chosen CW–S combination appear in the bottom left corner. The right plots represent graphically the statistical damage indicators in the form of overpressure-impulse and CW–S diagrams. Both show the conditions corresponding to the specified CW–S combination with a solid red dot, while the CW–S diagrams include also a diagonal dashed line indicating the approximated position of the fireball radius. Above this line, the Freidlander waveform is not valid, and the blast wave parameters are increasingly imprecise [1].

Figure 9. Schematic diagram of the three-story building under study, composed of equally spaced pillars and an outer enclosure wall, including: (a) the distances and angles used for the different floor levels ([Forumla omitted. See PDF.], 1, and 2), including the standoff distance, d, the real distance to the midpoint of the different levels, [Forumla omitted. See PDF.], and the corresponding angles of incidence, [Forumla omitted. See PDF.]; (b) schematic of the façade constructive details and dimensions; and (c) diagram of the equivalent façade element used in the SDOF analysis. L[Forumla omitted. See PDF.] denotes the height of Level i, representing the length of the pillars, and S is the spacing between pillars, representing the tributary loaded width.

Figure 10. Interface of the Blast Wave calculation assistant for a charge weight of 150 kg of TNT at the ISA mean sea level, showing the variation of the blast parameters with the standoff distance from the front façade (top table). The lower part of the assistant shows the blast parameters calculated at a point located at [Forumla omitted. See PDF.] m standoff distance and [Forumla omitted. See PDF.] m above the charge.

Figure 11. Metal beam calculation assistant showing results for a HEB 340 pillar with a length of 3 m and a spacing between pillars of 5 m. Note that, even though a HEB 350 is requested, which is not included in the norm, the assistant corrects down to the nearest normalized value, HEB 340.

Figure 12. Reinforced concrete beam calculation assistant showing results for a pillar of 45 × 45 cm[Forumla omitted. See PDF.] with a length of 3 m and a spacing between pillars of 5 m. The pillar is reinforced using 5 A36 steel reinforcement bars of 22.5 mm of diameter per side spaced apart 37 cm.

Figure 13. CW–S linear (a) and log-log (b) damage diagrams for reflected blast load on the façade of the first floor (Level 1): Case 1 (◯), Case 2 (✧), Case 3 (△).

Table 1

Comparison of the calculated temperature at constant volume, T, detonation pressure, $p_{\mathrm{CJ}}$, detonation velocity, $v_{\mathrm{CJ}}$, heat release at constant volume, $Q_v$, and explosive force, $F_e$, with the results provided by the European Standard EN 13631-15 [43] and by the thermochemical code W-DETCOM [49] for different explosive mixtures using the BKW–S EoS.

Explosive	Source	T [K]	${\mathit{p}}_{\mathbf{CJ}}$ [GPa]	${\mathit{v}}_{\mathbf{CJ}}$ [m/s]	${\mathit{Q}}_{\mathit{v}}$ [kJ/kg]	${\mathit{F}}_{\mathit{e}}$ [kJ/kg]
ANFO	CT	2592	7.14	5353	3845	943
	EN 13631-15	2586	-	-	3820	945
	W-DETCOM ${}^1$	2919	6.62	5326	3849	-
ANFO-Al	CT	3026	7.38	5442	4666	1009
	EN 13631-15	3060	-	-	4642	1020
	W-DETCOM ${}^1$	3370	6.55	5215	4655	-
Emulsion	CT	2112	15.3	6549	3263	766
	EN 13631-15	2099	-	-	3236	771
	W-DETCOM ${}^1$	2438	13.9	6758	3214	-
Dinamite I	CT	4173	25.03	7960	6452	1147
	EN 13631-15	4130	-	-	6338	1138
Dinamite II	CT	3165	23.58	7729	5049	987
	EN 13631-15	3151	-	-	4989	984

¹ Calculation performed assuming Chapman-Jouguet detonation.

Table 2

Composition [mass %], density, and oxygen balance of different explosive mixtures tested.

Component	ANFO	ANFO-Al	Emulsion	Dinamite I	Dinamite II
Aluminium	-	5	-	-	-
Ammonium nitrate	94	91	80	-	49
Cellulose	-	-	-	-	3
2,4-Dinitrotoluene	-	-	-	-	4
Nitrocellulose 12%	-	-	10	-	4
Nitroglycerin	-	-	-	45	20
Nitroglycol	-	-	-	45	20
Fuel oil	6	4	7	-	-
Sodium nitrate	-	-	5	-	-
Water	-	-	8	-	-
Density [kg/m${}^3$]	850	850	1300	1500	1500
Oxygen balance [%]	$-1.7$	$0.08$	$-5.57$	$-2.26$	$0.84$

Table 3

Probit functions used to estimate the probability of different types of primary and tertiary injuries. Pr is the probit value, $p^{\circ }$ [Pa] the peak overpressure, $p_{\mathrm{ef}}^{\circ }$ [Pa] the maximum effective overpressure, depending on the relative orientation of the person with respect to the shock wave, $p_1$ [Pa] the atmospheric pressure, $I/A$ [Pa · s] the impulse per unit area and m [kg] the weight of the person [9].

Effect	Probit Function
Primary injuries
Eardrum rupture	$\mathrm{Pr}=-12.6+1.52\ln p^{\circ }$
Death due to lung damage	$\mathrm{Pr}=5-5.74\ln \left(\frac{4.2}{p_{\mathrm{ef}}^{\circ }/p_1}+\frac{1.3}{i/\left(p_1^{1/2}{m}^{1/3}\right)}\right)$
Tertiary injuries
Death due to displacement and whole-body impact	$\mathrm{Pr}=5-2.44\ln \left(\frac{7380}{p^{\circ }}+\frac{1.3 \times {10}^9}{p^{\circ }i}\right)$
Death due to displacement and skull impact	${\displaystyle \mathrm{Pr}=5-8.49\ln \left(\frac{2430}{p^{\circ }}+\frac{4\times {10}^8}{p^{\circ }i}\right)}$

Table 4

Standoff distance, d, explosive charge, W, and scaled distance, Z, of the different case studies. The reference case is shown in blue.

Variables	Case 1	Case 2	Case 3
d (m)	12	20	25
W (kg)	30	150	300
Z (m/kg${}^{1/3}$)	3.86	3.76	3.73

Table 5

Incident load parameters and component damage indicators per floor. According to the PDC-TR-06-08 [15], the response limits for hot rolled structural steel can be defined in terms of the ductility ratio, $\mu$, and support rotation angle, $\theta$, as follows: B1—superficial $\{\mu ,\theta \}=\{1,-\}$; B2—moderate $\{\mu ,\theta \}=\{3,3^{\circ }\}$; B3—heavy $\{\mu ,\theta \}=\{12,{10}^{\circ }\}$; B4—hazardous $\{\mu ,\theta \}=\{25,{20}^{\circ }\}$. The reference case and worst-case scenario are indicated in blue and gray, respectively.

Level	Type	Variables	Case 1	Case 2	Case 3
0	Incident load parameters	$\Delta p$ (kPa)	168.30	182.50	186.80
		$I/A$ ($\mathrm{kPa}·\mathrm{ms}$)	406.70	724.40	922.20
		$d_{\mathrm{real}}$ (m)	12.17	20.10	25.08
		${\delta }_{\mathrm{real}}$ (deg)	9.46	5.71	4.57
	Damage level indicators	$\mu$ (-)	1.60	3.26	4.40
		$\theta$ (deg)	0.19	0.39	0.53
1	Incident load parameters	$\Delta p$ (kPa)	139.40	169.50	178.00
		$I/A$ ($\mathrm{kPa}·\mathrm{ms}$)	349.90	688.10	893.00
		$d_{\mathrm{real}}$ (m)	13.20	20.74	25.60
		${\delta }_{\mathrm{real}}$ (deg)	24.62	15.38	12.41
	Damage level indicators	$\mu$ (-)	0.90	1.64	2.10
		$\theta$ (deg)	0.08	0.15	0.19
2	Incident load parameters	$\Delta p$ (kPa)	110.20	152.00	165.40
		$I/A$ ($\mathrm{kPa}·\mathrm{ms}$)	293.00	630.50	845.50
		$d_{\mathrm{real}}$ (m)	14.71	21.73	26.41
		${\delta }_{\mathrm{real}}$ (deg)	35.31	23.03	18.78
	Damage level indicators	$\mu$ (-)	1.67	5.85	9.26
		$\theta$ (deg)	0.87	3.05	4.83

References

1. Kinney, G.F.; Graham, K.J. Explosive Shocks in Air; 2nd ed. Springer: New York, NY, USA, 1985.

2. Brode, H.L. Numerical solutions of spherical blast waves. J. Appl. Phys.; 1955; 26, pp. 766-775. [DOI: https://dx.doi.org/10.1063/1.1722085]

3. Grady, D. Fragmentation of Ring and Shells, the Legacy of N.F. Mott; Springer: Berlin/Heidelberg, Germany, 2006.

4. Ambrosini, R.D.; Luccioni, B.M. Craters produced by explosions on the soil surface. J. Appl. Mec.; 2006; 73, pp. 890-900. [DOI: https://dx.doi.org/10.1115/1.2173283]

5. Akhavan, J. The Chemistry of Explosives; 3rd ed. Royal Society of Chemistry: Cambridge, UK, 2011.

6. Cooper, P.W. Explosives Enginering; Wiley-VCH, Inc.: New York, NY, USA, 1996.

7. Meyer, R.; Köhler, J.; Homburg, A. Explosives; 6th ed. Wiley-VCH Verlag GmbH: Weinheim, Germany, 2007.

8. Ngo, T.; Mendis, P.; Gupta, A.; Ramsay, J. Blast loading and blast effects on structures—An overview. Electron. J. Struct. Eng.; 2007; 1, pp. 76-91. [DOI: https://dx.doi.org/10.56748/ejse.671]

9. TNO Green Book Methods for the Determination of Possible Damage to People and Objects Resulting from Releases of Hazardous Materials (CPR 16E); CIP-data of the Royal Library: The Hage, The Netherlands, 1992.

10. START. National Consortium for the Study of Terrorism and Responses to Terrorism. Global Terror-Ism Database. 2017; Available online: https://www.start.umd.edu/gtd (accessed on 3 September 2022).

11. US Department of Homeland Security (DHS). Reference Manual to Mitigate Potential Terrorist At-Tacks Against Buildings; 2nd ed. Federal Emergency Management Agency: Washington, DC, USA, 2011.

12. Luccioni, B.; Ambrosini, D.; Danesi, R. Blast load assessment using hydrocodes. Eng. Struct.; 2006; 28, pp. 1736-1744. [DOI: https://dx.doi.org/10.1016/j.engstruct.2006.02.016]

13. Anderson, C.E., Jr. An overview of the theory of hydrocodes. Int. J. Impact Eng.; 1987; 5, pp. 33-59. [DOI: https://dx.doi.org/10.1016/0734-743X(87)90029-7]

14. Han, Y.; Liu, H. Finite element simulation of medium-range blast loading using LS-DYNA. Shock Vib.; 2015; 2015, 631493. [DOI: https://dx.doi.org/10.1155/2015/631493]

15. Protective Design Center Technical Report PDC-TR 06-08; Single Degree of Freedom Structural Response Limits for Antiterrorism Design US Army Corps of Engineers: Washington, DC, USA, 2008.

16. Unified Facilities Criteria (UFC) 3-340-02; Structures to Resist the Effects of Accidental Explosions US Department of Defense: Washington, DC, USA, 2008.

17. Baker, W.E. Explosions in Air; University of Texas Press: Austin, TX, USA, 1973.

18. Kingery, C.N.; Bulmash, G. Airblast Parameters from TNT Spherical Air Burst and Hemispherical Surface Burst; US Technical Report ARBRL-TR-02555 Aberdeen Proving Ground, Ballistic Research Laboratories (BRL): Aberdeen, MD, USA, 1984.

19. Needham, C.E. Blast Waves; Springer: Berlin/Heidelberg, Germany, 2010.

20. Protective Design Center Software, US Army Corps of Engineers. Available online: https://www.nwo.usace.army.mil/About/Centers-of-Expertise/Protective-Design-Center/PDC-Software (accessed on 3 September 2022).

21. Sánchez-Monreal, J.; Vera, M. SimEx: SIMulador de EXplosiones mediante sistemas SDOF. Proceedings of the III Na-tional Congress of R&D in Defense and Security (DESEi+d 2015); Naval Military School of Marín, Pontevedra, Spain, 19–20 November 2015; pp. 993-1000.

22. Sánchez-Monreal, J.; Vera, M. SIMulador de EXplosiones mediante sistemas SDOF (SimEx): Desarrollo de la aplicación y nuevas funcionalidades. Proceedings of the IV National Congress of R&D in Defense and Security (DESEi+d 2016); University Defense Center, General Air Academy, San Javier, Spain, 16–18 November 2016; pp. 949-955. Available online: https://www.tecnologiaeinnovacion.defensa.gob.es/../Actas_DESEi+d2016_1.pdf" (accessed on 3 September 2022).

23. Sánchez-Monreal, J.; Vera, M. SimEx: Una herramienta para la evaluación rápida de los efectos de explosiones. Bol. Obs. Tecnol. Def.; 2017; 55, pp. 23-27. Available online: https://publicaciones.defensa.gob.es/../bot_55.pdf (accessed on 3 September 2022).

24. Cuadra, A.; Huete, C.; Vera, M. Combustion Toolbox: A MATLAB-GUI Based Open-Source Tool for Solving Gaseous Combustion Problems (v0.9.7). Zenodo: Geneve, Switzerland, 2022.

25. Ambrosini, R.D.; Luccioni, B.M.; Danesi, R.F.; Riera, J.D.; Rocha, M.M. Size of craters produced by explo-sive charges on or above the ground surface. Shock Waves; 2002; 12, pp. 69-78. [DOI: https://dx.doi.org/10.1007/s00193-002-0136-3]

26. Mott, N.F. A Theory of Fragmentation of Shells and Bombs; AC4613 Ministry of Supply: Great Britain, UK, 1943.

27. Mott, N.F. Fragmentation of H.E. Shells: A Theoretical Formula for the Distribution of Weights of Fragments; AC3642 Ministry of Supply: Great Britain, UK, 1943.

28. Mott, N.F. Fragmentation of shell cases. Proc. R. Soc.; 1947; A189, pp. 300-308.

29. Sánchez-Monreal, J.; Llamazares-Mendo, F.; Vera, M. Estimación de daños de explosiones sobre personas. Proceedings of the V National Congress of R&D in Defense and Security (DESEi+d 2017); Toledo, Spain, 22–24 November 2017; pp. 1689-1697. Available online: https://publicaciones.defensa.gob.es/../actas_v_congreso_id_2017.pdf (accessed on 3 September 2022).

30. González Duperón, E.; Martínez Gil, M.; Giménez Francés, B.; González Ferradás, E. Analysis of effects of detonation of explosive substances on humans from characteristic curves explosive mass-distance from the origin of the explosion. Process. Saf. Prog.; 2016; 35, pp. 233-240. [DOI: https://dx.doi.org/10.1002/prs.11794]

31. Protective Design Center Technical Report PDC TR 06-01; Rev 1, Methodology manual for the Single-Degree-of-Freedom Blast Effects Design Spreadsheets (SBEDS) US Army Corps of Engineers: Washington, DC, USA, 2008.

32. Biggs, J.M. Introduction to Structural Dynamics; McGraw-Hill College: New York, NY, USA, 1964.

33. Brode, H.L. Blast wave from a spherical charge. Phys. Fluids; 1959; 2, pp. 217-229. [DOI: https://dx.doi.org/10.1063/1.1705911]

34. Needham, C.E.; Crepeau, J.E. The DNA Nuclear Blast Standard (1 kT); Report No. DNA 5648-T Systems, Science and Software, Inc.: Alburquerque, Spain, 1981.

35. Hopkinson, B. British Ordnance Board Minutes 13565; The National Archives: Richmond, UK, 1915.

36. Sachs, R.G. The Dependence of Blast on Ambient Pressure and Temperature; Technical Report 466 Ballistic Research Laboratories: Aberdeen, MD, USA, 1944.

37. ISO 2533:1975 Standard Atmosphere; International Organization for Standardization: Geneva, Switzerland, 1975.

38. Gallo, E. Quasi static atmospheric model for aircraft trajectory prediction and flight simulation. arXiv; 2021; arXiv: 2101.10744

39. International Ammunition Technical Gudeline (IATG) 01.80 Formulae for Ammunition Management; UNODA Office for Disarmament Affairs: New York, NY, USA, 2013.

40. Friedlander, F.G. The diffraction of sound pulses I. Diffraction by a semi-infinite plane. Proc. R. Soc. London Ser. Math. Phys. Sci.; 1946; 186, pp. 322-344.

41. Euronorm 53–62 European Wide Flange Beams; CEN: Brussels, Belgium, 1993.

42. Norma UNE 31-002-94 Cálculo de las Principales Características Teóricas de los Explosivos; Asociación Española de Normalización UNE: Biddeford, Maine, 1994.

43. European Standard EN 13631-15 Explosives for Civil Use-High Explosives—Part 15: Calculation of the Thermodynamic Properties; UNE: Biddeford, Maine, 2005.

44. Arnold, W.; Rottenkolber, E.; Hartmann, T. DRAGON—The German Thermo-Chemical Code Based on the Becker-Kistiakowsky-Wilson Equation of State. Propellants Explos. Pyrotech.; 2022; e202100329. [DOI: https://dx.doi.org/10.1002/prep.202100329]

45. Hobbs, M.L.; Baer, M.R. Nonideal thermoequilibrium calculations using a large product species data base. Shock Waves; 1992; 2, pp. 177-187. [DOI: https://dx.doi.org/10.1007/BF01414640]

46. Mader, C.L. Detonation Properties of Condensed Explosives Computed Using the Becker-Kistiakowsky-Wilson Equation of State; Report LA-2900 Los Alamos Scientific Laboratory: Los Alamos, NM, USA, 1963.

47. Heuzé, O. Equations of state of detonation products: Influence of the repulsive intermolecular potential. Phys. Rev. A; 1986; 34, pp. 428-432. [DOI: https://dx.doi.org/10.1103/PhysRevA.34.428] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9897266]

48. López, L.M. Evaluación de la Energía de los Explosivos Mediante Modelos Termodinámicos de Detonación. Ph.D. Thesis; Universidad Politecnica de Madrid: ETSI Minas, Madrid, Spain, 2003.

49. Sanchidrián, J.A.; López, L.M. Calculation of the energy of explosives with a partial reaction model. Comparison with cylinder test data. Propellants Explos. Pyrotech.; 2006; 31, pp. 24-32. [DOI: https://dx.doi.org/10.1002/prep.200600003]

50. Cuadra, A.; Huete, C.; Vera, M. Desarrollo de un código termoquímico para la evaluación de las propiedades teóricas de explosivos (CT-EXPLO) y la estimación del rendimiento de motores cohete (CT-ROCKET). Proceedings of the IX National Congress of R&D in Defense and Se-curity (DESEi+d 2022); Army’s “General Morillo” Base, Pontevedra, Spain, 15–17 November 2022; accepted

51. Kamlet, M.J.; Jacobs, S.J. Chemistry of detonations. I. A simple method for calculating detonation proper-ties of CHNO explosives. J. Chem. Phys.; 1968; 48, pp. 23-55. [DOI: https://dx.doi.org/10.1063/1.1667908]

52. Politzer, P.; Murray, J.S. Some perspectives on estimating detonation properties of C, H, N, O compounds. Cent. Eur. J. Energetic Mater.; 2011; 8, pp. 209-220.

53. Beveridge, A. Forensic Investigation of Explosions; 2nd ed. CRC Press: Boca Raton, FL, USA, Tyalor & Francis Group: Boca Raton, FL, USA, 2012.

54. Westine, P.S. Explosive cratering. J. Terramechanics; 1970; 7, pp. 9-19. [DOI: https://dx.doi.org/10.1016/0022-4898(70)90130-8]

55. Baker, W.E.; Westine, P.S.; Dodge, F.T. Similarity Methods in Engineering Dynamics: Theory and Practice of Scale Modeling; Fundamental Studies in Engineering No. 12 Elsevier Science Limited: Amsterdam, The Netherlands, 1991.

56. Gurney, R.W. The Initial Velocities of Fragments from Bombs, Shells and Grenades; Report No 405 Aberdeen Proving Ground, Ballistic Research Laboratories (BRL): Aberdeen, MD, USA, 1943.

57. Carter, R.T.; Jandir, P.S.; Kress, M.E. Estimating the drag coefficients of meteorites for all Mach number regimes. Proceedings of the 40th Lunar and Planetary Science Conference, (Lunar and Planetary Science XL); The Wood-Lands, TX, USA, 23–27 March 2009.

58. Häring, I. Damage Analysis I: Probit Functions and Probability Distributions. Risk Analysis and Management: Engineering Resilience; Springer: Singapore, 2015; pp. 187-204.

59. Finney, D.J. Probit Analysis; 3rd ed. Cambridge University Press: Cambridge, UK, 1971.

60. Alonso, F.D.; Ferradás, E.G.; Pérez, J.F.S.; Aznar, A.M.; Gimeno, J.R.; Alonso, J.M. Characteristic overpressure–impulse–distance curves for the detonation of explosives, pyrotechnics or unstable substances. J. Loss Prev. Process. Ind.; 2006; 19, pp. 724-728. [DOI: https://dx.doi.org/10.1016/j.jlp.2006.06.001]