Development of a new empirical formula for prediction of triple point path

This paper develops a new empirical formula for the prediction of the triple point path in irregular shock reflection cases. Numerical simulations using a two-dimensional axisymmetric multi-material arbitrary Lagrangian–Eulerian formulation are employed to obtain the data of fluid density. Using the data of fluid density and nodal coordinates, the gradients of fluid density are determined and then used to generate numerical schlieren images. Based on these images, the triple point paths are derived and compared with the models of the Unified Facilities Criteria (UFC) and Natural Resources Defense Council (NRDC) as well as two models from the open literature. It is found that the numerically derived triple point paths are in good agreement with those predicted by a recently published model in the open literature for the typical ground range of shock wave propagation of up to 6 m. Considering the whole distance range, it is found that the agreement of different models of the triple point path with the numerical ones depends on the considered blast scenario, i.e., the scaled charge height. For small-scaled charge heights, the model of the UFC and the recently published model in the open literature are in better agreement with the numerical results than the other two models, whereas the NRDC model has the best agreement with the numerical results for large-scaled charge heights. Based on the numerical results, a new empirical formula is proposed for the prediction of the triple point path, which is valid for a wide range of the scaled charge heights from 0.5 to 3.5 m/kg1/3 and scaled ground distances up to 15 m/kg1/3.


Introduction
Blast parameters such as peak overpressure and maximum impulse can be calculated from the overpressure-time history. They play a crucial role in the blast-resistant design of structures interacting with shock waves. If a height of burst (HOB) explosion is concerned, the boundary conditions around the shock wave significantly affect the emerging overpressure-time history (incident blast, ground reflected blast, or Mach stem) and hence influence the protection strategy. Figure 1a illustrates the blast environment in an HOB explosion. The charge has a mass of W and is detonated Communicated  at a height of H c above the ground. The angle of incidence α is defined in Fig. 1a, where point A depicts a point of observation on the ground at a distance R from ground zero (GZ). Depending on whether the incident and reflected waves have merged to form a Mach stem or not, the ground range is subdivided into a Mach or a regular reflection region. The triple point (T) is the intersection point of the incident wave (I), reflected wave (R), and Mach stem (M). The origin of the triple point path is denoted as OTP. The quantities α min and R 0 are the minimum angle of incidence and the distance for the Mach stem formation, respectively. As the Mach stem travels along the ground surface (x-direction), the triple point rises in the y-direction and the height of the triple point H T increases. More details on the Mach reflection phenomena can be found, for example, in [1][2][3].
Regarding an HOB explosion, the gauge below the triple point (e.g., point B in Fig. 1a) will experience a single peak (Mach stem) in the overpressure-time history (Fig. 1b), whereas two successive peaks will appear in the overpressure-time history (Fig. 1c) measured at the gauge above the triple point (e.g., point C in Fig. 1a). These two peaks indicate the arrival of the incident and reflected waves, respectively. Thus, significantly different blast loads could be exerted on structures located in different regions. Therefore, an accurate prediction of the triple point path for different blast scenarios is of great importance.  [5], Natural Resources Defense Council (NRDC) [6], and Boutillier et al. [7]. UFC published ten curves of the triple point path, for example, Figure 2-13 of UFC 3-340-02 [4]. The scaled charge height varies from 0.4 to 2.8 m/kg 1/3 . The major drawback of the UFC model is the lack of a formula allowing extrapolation of the triple point prediction for scenarios other than the ten published curves. In addition, the origin of the triple point path is not provided. In order to overcome these two issues, Slavik [8] implemented an engineering model via the function *LOAD_BLAST_ENHANCE in LS-DYNA [9]. Note that all functions used in this study and mentioned in this paper are standard routines of the used software package LS-DYNA.
for R > R 0 . This formula is based on the observations from a variety of experiments using spherical charges. As depicted in Fig. 1a, the distance R 0 from ground zero to the origin of the triple point path can be calculated as: Kinney and Graham [5] proposed an empirical formula to determine the minimum angle of incidence as: where M i is the Mach number of the incident wave. In order to calculate the distance R 0 without information of α min , Needham [2] proposed an empirical formula as  [2,7]. It should be noted that a high explosive is about twice as effective at forming an air blast compared to a nuclear detonation. NRDC proposed also a model of the triple point path. It is expressed as: The quantity R 0 is the scaled quantity of R 0 , i.e.: It should be noted that the NRDC model is valid for the scaled charge height from 0 to 8 m/kg 1/3 and for charge masses of TNT ranging from 0.1 kt to 25,000 kt. The used units for length and mass are metre and kiloton, respectively. Boutillier et al. [7] conducted new experiments in order to develop a new empirical formula to predict the triple point path for smooth and hard reflecting surfaces. It is given as: The quantity Eq TNT is the TNT equivalence of C4. It is noted that this formula is valid for the scaled charge height from 0.246 to 1.729 m/kg 1/3 . For this equation to be valid, length and mass are to be given in centimetre and kilogram, respectively.

Numerical derivation of triple point path
This section aims to develop a new empirical formula to predict the triple point path. Firstly, numerical simulations are conducted with the software LS-DYNA to calculate the fluid density field for different HOB scenarios. Secondly, a corresponding sequence of the numerical schlieren images is generated by using the data of fluid density and nodal coordinates obtained from numerical simulations. These images are employed to obtain the triple point path. Thirdly, a curvefitting method is used to derive the formula to predict the triple point path.

LS-DYNA simulations
Two-dimensional (2D) axisymmetric multi-material arbitrary Lagrangian-Eulerian (MM-ALE) simulations are run The spherical form of the charge is modelled via the function *INITIAL_VOLUME_FRACTION_GEOMETRY. At first, a one-dimensional (1D) numerical simulation is conducted with an element size of 0.4 mm. This element size provides 131 elements across the charge radius of 1 kg TNT. It is worth pointing out that a sufficiently small element size is vital to accurately model the explosive energy release. As a rule of thumb, at least ten elements across the charge radius should be used for the meshing of the explosive charge [10]. However, at least 50 elements across the charge radius are recommended for modern high-performance computing. In addition, a mesh sensitivity study is performed and three element sizes of 1.6 mm, 0.8 mm, and 0.4 mm are used for the study. Figure 3 depicts the time histories of the fluid density at the gauges, which are located at a distance R of 0.75 m and 0.8 m from the charge, respectively. It is observed that the numerical results of the fluid density using all three element sizes are nearly identical. The element size of 0.4 mm is selected for the 1D numerical simulations.
The simulation terminates before the shock wave impinges on the ground surface. Then, the 1D numerical results are mapped to the 2D axisymmetric model (Fig. 2). The use of the mapping techniques is to enable highly fine meshes in the 1D geometry and comparatively coarse meshes in the 2D geometry. They are employed based on two considerations. On the one hand, the accuracy of the numerical results can be ensured. On the other hand, the computational demand is still manageable. In order to evaluate the effect of the advection method on the numerical results, predictive simulations are conducted with all four advection methods, i.e., the control parameter METH is set to 1, 2, − 2, and 3 in the function *CONTROL_ALE. Based on the observations on the density histories at several gauges, METH − 2 is adopted for the 1D numerical simulations. This second-order accurate Van Leer with half index shift (HIS) advection method is recommended by the software LS-DYNA since the function *MAT_HIGH_EXPLOSIVE_BURN is used in the 1D model. The default method (METH 2) is selected for the 2D numerical simulations because all advection methods achieve nearly identical results. Furthermore, the function *MAT_NULL is used instead for the 2D numerical simulations owing to the mapping technique from 1D to 2D. Compared to METH 2, the monotonicity condition in METH − 2 is relaxed during the advection process, which aims to better preserve the material interfaces for materials defined by *MAT_HIGH_EXPLOSIVE_BURN [9].
The high explosive charge is modelled via *MAT_HIGH_EXPLOSIVE_BURN with the Jones-Wilkins-Lee (JWL) equation of state (EoS). The pertinent material parameters of TNT are listed in Table 1 for the strength model and EoS, where ρ e is the density of the explosive material, D is the velocity of detonation, p CJ is the Chapman-Jouguet pressure, A, B, R 1 , R 2 , and ω are material constants, and e 0 is the specific internal energy. The air is modelled via *MAT_NULL and regarded as an ideal gas with a linear polynomial EoS.
In order to verify the 2D axisymmetric numerical model, a mesh sensitivity study is conducted, in which three element sizes of 4 mm, 8 mm, and 16 mm are used. Figure 4

Numerical schlieren images
Based on the data of fluid density (ρ) and nodal coordinates (x, y), a corresponding sequence of numerical schlieren images is generated. These data are obtained from numerical simulations at predefined tracer points in the ALE mesh. Since the tracer points are defined as fixed in space (via TRACK 1 in the function *DATABASE_TRACER), the tracer point coordinates are identical to the nodal coordinates. The distance between the tracer points is 2 cm. Hence, in total 0.3 million tracer points are defined. They are used to determine the gradient of fluid density, which can be used to define the shock profiles. The magnitude |∇ρ| of the gradient of fluid density is calculated as: Furthermore, a nonlinear shading function Φ is used to accentuate the weak flow features, where |∇ρ| max is the maximum value of the gradient of fluid density and k is a constant used for the flow visualization, in which the greyscale fringe levels in the numerical schlieren images are accurately adjusted. The darker the image, the larger the gradient of fluid density. More description of the numerical schlieren algorithm can be found in [13]. Figure 6 demonstrates how the triple point can be optically detected by means of numerical schlieren images for the scenario of W 1 kg TNT and H c 1 m. In this case, a value of k 20 is used for the flow visualization. At t 0.5 ms (Fig. 6a), the incident wave (I) propagates spherically. Once the shock wave impinges on the ground surface, a reflected wave (R) from the ground surface is generated. At a distance of 1.0 m from ground zero, the Mach stem begins to form, i.e., at t 1.1 ms (Fig. 6b). It is indicated by the intersection of the incident (I) and reflected (R) waves at the point OTP on the ground surface (see zoomed view in Fig. 6b). As the shock wave further propagates, the triple point (T) rises ( Fig. 6c-f). Below the triple point (T), the Mach stem (M) is developed by the coalescence of the incident (I) and reflected (R) waves. In addition, not only the primary incident (I 1 ) and reflected (R 1 ) waves but also the secondary incident (I 2 ) and reflected (R 2 ) waves are evident in the flow field (Fig. 6d-f). Applying this method to consecutive numerical schlieren images, the triple point path is derived. This method implies that a Mach stem only becomes visible once the selected distance between the tracer points is smaller than or equal to 2 cm.
Using the concept of comparison in [7], the numerically derived triple point paths are illustrated in Fig. 7 (labelled as "Schlieren"), compared to the other four models of triple point path (labelled as "Kinney and Graham", UFC, NRDC, and "Boutillier et al.", respectively). Six different blast scenarios are involved, in which the charge mass (W 1 kg TNT) is detonated at a height H c from Considering the whole distance range, the numerical triple point path for H c 0.5 m (Fig. 7a) agrees well with the ones predicted by the models of UFC and Boutillier et al. For H c 1.5-3.0 m (Fig. 7c-f), it agrees reasonably well with the ones predicted by the NRDC model. For H c 1.0 m (Fig. 7b), it lies between the ones predicted by the Boutillier et al. and NRDC models. It is noted that, for this charge height, the triple point paths predicted by the models of Boutillier et al. and Kinney and Graham are identical. In general, the NRDC model predicts a triple point path serving as the lower bound, whereas the UFC (Kinney and Graham) model serves as the upper bound for small (large)-scaled charge heights. The dif- ference between these models of the triple point path is most likely due to the different properties (roughness and hardness) of the ground surface at the test sites of the respective experiments used for the development of the models [7]. It is worth pointing out that the models of Kinney and Graham, UFC, and NRDC had not provided any information on the roughness and hardness of the reflecting surface [7].
However, the location of the triple point path depends not only on the charge height and mass but also on the roughness and hardness of the reflecting surface. It was revealed in extensive experiments (e.g., in [14][15][16]) that high surface roughness delayed the Mach stem formation and lowered the height of the triple point. It is noted that as H c becomes larger (e.g., 3.0 m, Fig. 7f), the numerical triple point path is likely In summary, the observations mentioned above reveal that the agreement of the models of the triple point path with the numerical ones depends on the considered blast scenario, i.e., the scaled charge height. For small-scaled charge heights, the models of UFC and Boutillier et al. are in better agreement with the numerical results than the other two models (Kinney and Graham, NRDC), whereas the NRDC model has the best agreement with the numerical results for large-scaled charge heights. These differences confirm the intention of this study, i.e., that a new tool for the prediction of the triple point path is required, which is applicable to a wide range of ground distances and charge heights.

Curve-fitting method
The triple point path is likely to follow a parabolic-like shape. If the charge height increases or the charge mass decreases, the triple point path flattens. Furthermore, as observed from the UFC curves, the scaled height of the triple point decays exponentially with the scaled charge height. The model of Kinney where Based on the numerical data for blast scenarios, i.e., charge mass W 1 kg TNT and charge height H c 0.5-3.5 m, the coefficients a, b, c, d, and f of (18) are determined using the MATLAB Curve Fitting Toolbox™ [17], i.e., a 0.2104, b 0.6003, c 0.248, d 2.534, and f 0.009393. The coefficient of determination is the square of the correlation between the response values and predicted response values. In this case, it is equal to 0.9998, which indicates that there is a good correlation between the numerically derived values of the heights of the triple point and the ones predicted by (18). It should be noted that only the data of triple points before the shock wave impinges on the boundaries are used for the development of the empirical formula. Figure 8 illustrates the fitted surface versus the numerical data used for the calibration. The scaled charge height H c /W 1/3 varies from 0.5 to 3.5 m/kg 1/3 . The scaled distance R/W 1/3 varies from 0 to 15 m/kg 1/3 . It is worth mentioning that two assumptions are made during the derivation of the new empirical formula. One is the scaled height of the triple point decays exponentially with the scaled charge height, which is based on the observation from the UFC curves. Another assumption is that the scaled height of the triple point has a parabolic relationship with the scaled distance from ground zero, which is proposed by Boutillier et al. [7].
It is worth pointing out that two main differences exist between this new empirical formula and the one of Boutillier et al. Firstly, Boutillier et al. [7] used own experimental data to develop the empirical formula. This study uses data from numerical simulations. This provides more data points for the curve fitting. The formula of Boutillier et al. is valid for the scaled charge height up to 1.729 m/kg 1/3 and for scaled ground distances up to 6 m/kg 1/3 . By virtue of numerical simulations, the application range of (18) is extended to 3.5 m/kg 1/3 for the scaled charge height and 15 m/kg 1/3 for the scaled ground distance. Secondly, for the formula The triple point paths predicted by (18) are also plotted in Fig. 7 (labelled as CFM) for the different blast scenarios. As observed in Fig. 7, the fitted curves are in good agreement with the numerically derived triple point paths. In order to evaluate the reliability of the formula developed in this study, a possibility would be the use of the experimental data, especially for large ground distances. However, such experimental data are scarce in the open literature.

Summary
This paper proposes a new empirical formula to predict the triple point path in HOB explosion scenarios. Two-dimensional axisymmetric multi-material arbitrary Lagrangian-Eulerian numerical simulations are conducted to calculate the fluid density field. Using the data of fluid density and nodal coordinates, the gradients of fluid density are determined and used to generate the numerical schlieren images. Based on these images, the triple point paths are derived, which are in good agreement with those predicted by the model of Boutillier et al. [7] for the typical ground range (up to 6 m) of the shock wave propagation used in this model. Considering the whole distance range, it can be seen that the agreement of different models of the triple point path with the numerical ones depends on the considered blast scenario, i.e., the scaled charge height. For small-scaled charge heights, the models of UFC and Boutillier et al. are in better agreement with the numerical results than the other two models (Kinney and Graham, NRDC), whereas the NRDC model has the best agreement with the numerical results for large-scaled charge heights. Based on the numerical data of triple points, a new empirical formula is proposed to derive the triple point path. It is valid for the range of scaled charge heights from 0.5 to 3.5 m/kg 1/3 and scaled distances up to 15 m/kg 1/3 . For general use, this formula is given in terms of the scaled quantities. Unlike the available models of the triple point path, the dependencies of the scaled height of the triple point both on the scaled charge height and on the scaled distance from ground zero are incorporated into a single equation.