Predicting specific impulse distributions for spherical explosives in the extreme near-field using a Gaussian function

Accurate quantification of the blast load arising from detonation of a high explosive has applications in transport security, infrastructure assessment and defence. In order to design efficient and safe protective systems in such aggressive environments, it is of critical importance to understand the magnitude and distribution of loading on a structural component located close to an explosive charge. In particular, peak specific impulse is the primary parameter that governs structural deformation under short-duration loading. Within this so-called extreme near-field region, existing semi-empirical methods are known to be inaccurate, and high-fidelity numerical schemes are generally hampered by a lack of available experimental validation data. As such, the blast protection community is not currently equipped with a satisfactory fast-running tool for load prediction in the near-field. In this article, a validated computational model is used to develop a suite of numerical near-field blast load distributions, which are shown to follow a similar normalised shape. This forms the basis of the data-driven predictive model developed herein: a Gaussian function is fit to the normalised loading distributions, and a power law is used to calculate the magnitude of the curve according to established scaling laws. The predictive method is rigorously assessed against the existing numerical dataset, and is validated against new test models and available experimental data. High levels of agreement are demonstrated throughout, with typical variations of <5% between experiment/model and prediction. The new approach presented in this article allows the analyst to rapidly compute the distribution of specific impulse across the loaded face of a wide range of target sizes and near-field scaled distances and provides a benchmark for data-driven modelling approaches to capture blast loading phenomena in more complex scenarios.


Introduction and context
Accurate appraisal of the viability of structures and protective systems following the close-in detonation of a high explosive is only possible if the analyst has an accurate model detailing both the distribution and magnitude of the imparted load. In these so-called extreme near-field scenarios, the loading is highly localised, whilst deformation and damage occur on timescales which are orders of magnitude longer than the duration of loading. It is clear, therefore, that the temporal form of the loading is of lesser importance, yet a description of the spatial distribution of the peak loading, that is, peak specific impulse, is critical.
Quantification of blast loading became a particular topic of interest to the engineering community in the years following the second world war. One such study by Kingery (1966), later revisited by Kingery and Bulmash (1984) -hereafter abbreviated as KB -compiled existing blast pressure measurements and used geometric scaling laws, formulated independently by Hopkinson (1915) and Cranz (1926), to develop a semi-empirical tool for predicting normally reflected and incident blast wave parameters for hemispherical surface bursts and spherical free air bursts. Relationships are provided from effectively contact explosions, out to Z ≈ 40 m/kg 1/3 , where Z is scaled distance, given by S W / 1/3 , where S is the distance from the centre of the explosive to the target (in m, commonly termed stand-off distance), and W is the mass of explosive (in kg). However, as noted by Esparza (1986), direct measurements of blast wave parameters in close proximity to the charge were either 'non-existent or very few', and near-field semi-empirical predictive data were inferred from non direct measurements such as smoke trails (Dewey, 1964) or rudimentary numerical analyses (Brode, 1955).
The KB semi-empirical predictions have been shown to be accurate for medium to far-field scaled distances (Rigby et al., 2014), however in a review of uncertainty and conservatism in blast load predictive models, Bogosian et al. (2002) found that the KB curves 'substantially underpredicted the recorded data for Z m kg < 2.4 1 3 / / '. Furthermore, numerical modelling studies have shown there to be order-of-magnitude differences between the semi-empirical and computational blast parameters in the near-field (Shin et al., 2015). Whilst the KB approach and similar (semi-) empirical methods are attractive in terms of their very low computational demand and ease of use, a demonstrated lack of accuracy in the near-field renders them unsuitable when considering the effects of close-in detonations.
Alternatively, experimental methods can be used to directly measure near-field reflected specific impulse and specific impulse distributions, for example, impulse plugs (Huffington and Ewing 1985;Nansteel et al., 2013), Hopkinson pressure bars (Cloete and Nurick, 2016;Edwards et al., 1992;Piehler et al., 2009, Rigby et al., 2015aTyas et al., 2016), and flush-mounted pressure gauges (Aune et al., 2016). Near-field blast load measurements are onerous for a number of reasons. Firstly, the high magnitude of loading necessitates the use of robust support structures and protective housing for sensitive equipment (which itself is required to record in the MHz frequency range). Secondly, the measurements themselves are highly variable owing to the presence of surface instabilities in the early stages of expansion of the detonation products (Rigby et al., 2020). It is therefore not practical to develop a predictive approach based solely on physical testing, however experimental data remains a fundamental requirement for validation of numerical modelling schemes. Finite element (FE) and computational fluid dynamics (CFD) approaches have been used recently to investigate near-field blast loading, with the results generally demonstrating good agreement with experimental data where it is available (Rigby et al., 2018;Shin et al., 2014a;Whittaker et al., 2019). Despite research into near-field blast loading currently being hampered by a lack of well-controlled experimental validation data (Tyas, 2019), FE/CFD analyses can provide data at considerably higher spatial and temporal resolution than experimental studies and are therefore suitable tools with which to develop a refined predictive approach. However, physics-based models have a relatively high computational demand, and are unsuitable for probabilistic, risk-based analyses.
An outstanding challenge for the blast protection community is the provision of techniques for predicting near-field specific impulse distributions in a way which is both accurate and computationally inexpensive. An appropriate technique, which is adopted in the current study, is to use validated CFD analyses to create a dataset from which a predictive method can be developed. This process, known as data-driven surrogate modelling, is widely used in the field of computational dynamics, for example, Loy et al. (2017), Both et al. (2019), but has not yet been successfully applied in the field of near-field blast load prediction.
In this article, a regression analysis is performed on normalised specific impulse curves and a new predictive equation is derived, given by the product of a power law and Gaussian function. The results and equations are presented as scaled parameters for a range of near-field scaled distances, allowing the method to be used for any sized explosive charge. The accuracy of this new approach is rigorously tested against the existing CFD dataset, unseen CFD analyses and available experimental data from the literature. It is demonstrated that localised peak specific impulse and area-integrated total impulse (for a range of target sizes) are typically predicted to within 10%.

APOLLO Blastsimulator
Numerical analyses were performed using APOLLO Blastsimulator (hereafter shortened to Apollo), an explicit CFD software specialising in the simulation of dynamic flow problems; focussing on detonations, blast waves and gas dynamic (Fraunhofer, 2018). Apollo solves the conservation equations for transient flows of compressible, inviscid and non-heat conducting, inert or chemically reacting fluid mixtures. The HLL-type Reimann solver is particularly suited for problems involving rapid changes in pressure and density, and high speed flows, and is second order accurate owing to the tri-linear reconstruction of cell-centred conservative variables.
Apollo is equipped with an inbuilt dynamic mesh adaptation (DMA) strategy, as well as a 1D solver, 1D-to-3D mapping and 3D-to-3D staged mapping (see Figure 1). The user specifies a 'zone length', L, which represents the element size at the coarsest resolution level, and a maximum resolution level, N max , such that the element size ('cell length'), h, at any resolution level, N, is given as h L N = 2 − , that is, elements at each successive resolution level are half the size of those in the previous resolution level, in all dimensions. The 'ultimate cell length', that is, the element size at the highest resolution level, is therefore given as h L N = 2 − max .  DMA is used to locally refine specific areas of the domain which require a finer spatial resolution, for example, those elements which comprise, and are close to, the shock front, or elements at the interface between ambient air and embedded objects. Levels are adapted automatically based on the time-dependent, zone-specific gradients of pressure, temperature, density, Mach number and concentrations (Klomfass, 2018). This has the benefit of directing computational resource appropriately whilst reducing demand in areas of the mesh where it is not needed. Additionally, staged mapping can be used to define a series of computational sub-domains (or 'stages') with increased mesh resolution. The analysis will transition from one stage to the next once a predefined condition has been met, typically relating to the position of the shock front. Finally, Apollo has the option to progressively expand a computational domain according to a pre-defined velocity-time function (e.g. the rate of expansion of a freely expanding spherical blast wave), however this is not usually required for close-in detonations.
PE4 is a plastic explosive comprising 87% RDX and 13% mineral oil binder (Tyas, 2019). The explosives were modelled as 100 g PE4 spheres throughout the mesh sensitivity study, validation, and data harvesting aspect of this study, using Apollo's in-built model for C4 explosive, as the two explosives are of nominally identical composition (Bogosian et al., 2016). The pressure-densitytemperature relationship of the post-detonation explosive products is given by the Jones Wilkins Lee equation of state (Lee et al., 1968): where ρ and ρ 0 are density and initial density, R is a gas constant ( R = 365 J/(kg·K) for PE4), and T is absolute temperature. The constants C 1 , C 2 , R 1 and R 2 are assigned the values 734.60 GPa, 8.86 GPa, 4.79 (−) and 1.06 (−) respectively for PE4. The air is modelled as a thermally perfect gas, p RT = ρ , such that under ambient conditions (288 °K), ambient pressure is p = 101.3 kPa. Ambient pressure has been subtracted from all pressure histories in this article to present over-pressures and impulses only. Afterburn has been shown to influence the development of loading in near-field blast scenarios (Tyas et al., 2016), and so was modelled using the Klomfass Afterburning (KAB) Model. The KAB model defines an empirical relationship for global combustion rate based on volume and rate of change of volume of the explosive fireball, whereby the total mass of combustion products generated in a time-step is obtained (Fraunhofer, 2018). In all numerical simulations the explosives were centrally detonated and the mass of the detonator was not included. The Chapman-Jouguet detonation model was used, with the progression rate a function of: local burning velocity, defined by the Chapman-Jouguet condition; local detonation product density; and local un-reacted explosive density. This model allows the detonation front to propagate in a physically meaningful manner (Fraunhofer, 2018).

Mesh sensitivity and validation
Mesh sensitivity analysis. Prior to validating Apollo results against available experimental data, a mesh sensitivity study was conducted with the aims of determining the required element size to achieve convergence and identifying suitable combinations of zone length and resolution level.
A series of simulations were run for a 100 g PE4 sphere (24.3 mm radius) located at a stand-off distance, S , of 0.05 m from the centre of the charge to a rigid reflecting wall. This corresponds to a scaled distance of Z = 0.108 1 3 m / kg / , which is the shortest scaled distance to be used in the subsequent modelling study. The domain was 0.35 0.35 0.35 × × m, and quarter symmetry was utilised, with symmetry planes located in the directions orthogonal to the reflecting wall, originating at the centre of the charge. Numerical pressure gauges were placed along the rigid surface out to a distance of 0.28 m from the point normal to the explosive. Two hundred gauges in total were placed at evenly-spaced angles of incidence, θ , where the angle of incidence is defined as the angle between the outward normal of the reflecting surface and the direct vector from the explosive charge to that point (Rigby et al., 2015b). The gauge which is normal to the explosive therefore has an angle of incidence of θ = 0  , and the most remote gauge has an angle of incidence of θ = 80  . A total of 13 models were run in the initial sensitivity analysis, with details of the meshes provided in Table 1. Apollo's auto-staging was used throughout.
Results from the mesh sensitivity analysis are shown in Figure 2, where each sub-figure shows: peak specific impulse for the normal gauge; area-integrated impulse out to a radius of 0.28 m from Table 1. Ultimate cell length (element size at highest resolution level) and number of elements (between charge centre and normal gauge location) for initial mesh sensitivity study, Z = 0.108 m/kg 1/3 . (a) Ultimate cell length (mm).

Resolution level
Zone length ( the normal gauge; and total analysis time (termed 'wall time' in Apollo), all plotted against the ratio of stand-off distance to ultimate cell length. The results suggest that an ultimate cell length of S /16 is sufficient to achieve convergence for both specific and total impulse, whilst maintaining a practical analysis time: a decrease in cell length from S /16 to S / 40 increases the total impulse by <3%, whilst increasing the wall time by approximately a factor of 25. For the situation considered in this sensitivity analysis, it appears as though specific impulse and total impulse are predominantly influenced by ultimate cell length, and zone length has a negligible effect. This suggests that elements on the loaded surface are refined to their highest resolution level for all models considered. The benefits of DMA, however, can be seen when considering wall time: at S /10, a clear reduction in analysis time is observed when going from a zone length of 0.01 m to 0.02 m. Therefore a meshing strategy of 50 mm zone length and resolution level of 4 will be implemented, extending to a 100 mm zone length and resolution level of 5 (maintaining an ultimate cell length of 3.125 mm, i.e. S /16 throughout) to accommodate a number of larger domain sizes required to populate the CFD dataset in subsequent sections of this paper.
Experimental validation. The Characterisation of Blast Loading apparatus (Clarke et al., 2015) was designed to measure spatial and temporal distribution of loading from buried  and free air explosions (Rigby et al., 2015a). The apparatus utilises two perpendicular arrays of Hopkinson (1914) pressure bars (HPBs), mounted such that their faces lie flush with the underside of a large (100 mm thick, 1400 mm diameter) steel plate, which acts as a nominally rigid reflecting surface under which the explosive charges are detonated. 17 HPBs are used in total: one bar is located directly above the charge (which itself is aligned with the centre of the target plate); and four bars are located at each distance of 25 mm, 50 mm, 75 mm and 100 mm from the centre. Numerical integration of the recorded pressure signals enables peak specific impulse distributions over the central 200 mm diameter region of the target to be generated. Two different configurations were used to validate Apollo: 100 g PE4 spheres were used throughout, with three tests at 80 mm stand-off (Rigby et al., 2019a) and three tests at 380 mm stand-off (Rigby et al., 2020), corresponding to scaled distances of 0.172 m/kg 1/3 and 0.819 m/kg 1/3 respectively. Apollo models were run with 50 mm zone length and resolution level of 4 as informed by findings from the mesh sensitivity analysis. The reflecting surface was modelled as a rigid wall and quarter symmetry was used as before. The domain was extended to 1.5 1.5 1.5 × × m to prevent expansion waves from the edge of the domain reaching the gauge locations, and to establish a consistent domain size for use in subsequent models. Numerical overpressure and specific impulse histories at 0, 25 and 50 mm from the target centre are compared against those from one of the 80 mm stand-off experiments in Figure 3. It can be seen that the magnitudes and general form of the numerical pressure and impulse histories match the experiments closely, which gives confidence that Apollo is correctly modelling the mechanisms of normal and oblique reflection in near-field blast scenarios.
Additionally, Figure 4 shows compiled peak specific impulse distributions for the 80 mm standoff and 380 mm stand-off tests with the CFD data evaluated at 0.4° increments of angle of incidence. Also shown are the mean experimental recordings at each measurement location for ease of comparison (three data points at the target centre and 12 data points at each other distance from the centre). Again, the numerical results closely match the experimental data, in both form and magnitude, although the results appear to be in better agreement with the upper bound of the experimental results. Deriving a model based on Apollo CFD results will therefore have a slight inbuilt conservatism, however this is still representative of the typical variations seen in blast loading from close-in detonations, as exhibited here.

Generation of dataset
Apollo was used to generate specific impulse distributions along a rigid reflecting surface, at angles of incidence between 0-80°, with 200 numerical gauges placed at regular increments of angle of incidence as per the validation exercise. The analyses were formed into two groupings. Firstly, an initial dataset was generated for extreme near-field conditions, 0.11-0.21 m/kg 1/3 , at increments of ∼ 0.025 m/kg 1/3 . Here, the domain size was fixed at 1.5 1.5 1.5 × × m, using quarter symmetry, with 100 mm zone length and resolution level of 5. This results in the same ultimate cell length as the validation model and recommendations from the mesh sensitivity study whilst maximising computational efficiency through the use of the DMA feature in Apollo.  / / as 0, 25 and 50 mm distance from the target centre: (a) 0 mm from centre, θ = 0 º , (b) 25 mm from centre, θ =17 º and (c) 50 mm from centre, θ = 32 º .
The dataset was then expanded to incorporate intermediate near-field conditions, with analyses run between 0.24-0.55 m/kg 1/3 , again at increments of ∼ 0.025 m/kg 1/3 . The full dataset comprised 18 analyses in total, with 3600 data points between 0.11-0.55 m/kg 1/3 and 0-80°. For the Z > 0.21 1 3 m / kg / analyses, a separate domain was specified for each model, with equal side lengths of an integer number of zone lengths, given by: , where S and L are stand-off and zone length as previously, such that S tan (80 )  is the distance from the centre of the target to the most remote gauge. The domain length was set a factor of 1.3 greater than this distance (rounded up to the nearest multiple of zone length), again to prevent edge expansion waves from reaching the most remote gauge location during the analysis. As with previous models, quarter symmetry was used, with 100 mm zone length and resolution level 5.
CFD results are shown in detail in Figure 5 for the reduced and full datasets. Here, a Savitzky-Golay filter (Savitzky and Golay, 1964) has been applied to the specific impulse distributions to remove small spurious oscillations induced by DMA, which were more apparent in the analyses at larger scaled distances. The Savitzky-Golay filter is a low-pass data-smoothing method based on local least-squares polynomial approximation, and is particularly suited for removing high frequency noise with the advantage of preserving the original shape and features of the signal (Acharya et al., 2016).

Development of a data-driven predictive approach
Overview This section details the development of a data-driven predictive approach for computing near-field specific impulse distributions from spherical PE4 charges. Firstly, the relationship between peak specific impulse and angle of incidence is examined, and three different curve types are trialled to establish the most suitable candidate for representing the spatial variation of loading. Secondly, the relationship between scaled distance and peak specific impulse is considered and a power law fit is determined. The proposed predictive method is given as the product of these two relationships: spatial variation, and variation with scaled distance. It was noted by Wilson et al. (2018) that 'further research into analysis scenarios at varying angles of incidence are required . . . at a scaled distance less than 0.16 m/kg 1/3 ', hence initial efforts are focussed on the reduced dataset, and then extended to the full dataset.

Spatial variation of normalised impulse
From examination of Figure 5(b) and (d), it appears as though that the relationship between peak specific impulse and angle of incidence is similar for all scaled distances considered in this article. Dividing each through by the respective peak specific impulse at 0  , that is, the normally reflected peak specific impulse, i r , allows the normalised peak specific impulse distribution to be plotted, as shown in Figure 6. Here, it can be seen that the curves form a consistent grouping, particularly for the reduced dataset, which suggests that it is possible to define a single curve to describe the shape of the distributed specific impulse with respect to angle of incidence. Three candidate curve-types have been identified: the semi-analytical trigonometric functions of Henrych (1979), and Randers-Pehrson and Bannister (1997); and the Gaussian function proposed by Pannell et al. (2019). The semi-analytical work of Henrych (1979), later developed by Remennikov et al. (2017), is based on the model of instantaneous detonation of a spherical explosive and subsequent expansion of the gaseous detonation products into a vacuum. Here the loading is assumed to develop entirely from momentum transferred by the detonation products as they impact the loaded face. With some cancelling of constants, the Henrych (1979) equations can be expressed in normalised form: where i is specific impulse divided by normally reflected specific impulse, and θ is angle of incidence as previously. Randers-Pehrson and Bannister (1997) present a trigonometric expression to account for angle of incidence effects, which after some manipulation can be presented as a function of normalised specific impulse and angle of incidence:  coefficients must be calculated separately for each and every angle of incidence under consideration. Furthermore, this relationship is also scaled-distance dependent. Finally, Pannell et al. (2019) presented work on using machine learning to fit Gaussian probability density functions to near-field specific impulse distributions. Here, a modified version of the equations as a function of angle of incidence is presented: where A G is the Gaussian width parameter, which controls the rate of decay of the Gaussian curve with respect to angle of incidence. Two other parameters form a traditional Gaussian function, that is, those which determine the height and position of the centre of the peak respectively. Due to the normalised and rotationally-symmetric nature of the loading considered here, the height coefficient is unit throughout, and the position coefficient is zero as each Gaussian curve is centred around θ = 0  . Whilst the loading is considered rotationally-symmetric, the method can still be used to model off-centre charge placements relative to the target surface, as detailed in section 4.2 Figure 7 shows an assessment of the quality of fit for each of the normalised predictions compared against the reduced CFD dataset. For the Randers-Pehrson and Bannister (1997) model, A R and A S were determined at every angle of incidence (200 values between 0 and 80°) for each scaled distance using the 5th-order and 11th-order polynomial equations for i Z r ( ) and i Z s ( ) provided in Table 4-2 of Shin et al. (2014b). Negligible difference was found between the A R and A S versus θ relationships averaged across all scaled distances, and the A R and A S versus θ relationship for the median scaled distance ( Z = 0.16 1 3 m / kg / ) of the reduced dataset. This low dependency of A R and A S with Z is to be expected for this particular grouping of scaled distances owing to the similarity of the normalised impulse curves as shown in Figure 6(a). Even given this simplification, however, the process for determining A R and A S is relatively involved and required a minimum of 200 computations.
For the Pannell et al. (2019) Gaussian model, A G was determined through a non-linear regression analysis in which the root mean squared error (RMSE) was minimised using a trust region reflective method (Branch et al., 1999). For the reduced dataset, the regression-estimated best fit yielded an A G value of 1829.
Also shown in Figure 7 are the model residuals (CFD value minus model prediction), with mean absolute error (MAE) and RMSE provided for each fit against the mean CFD distribution. The Henrych (1979) and Randers-Pehrson and Bannister (1997) trigonometric identities generally under-predict specific impulse for intermediate angles of incidence (which would result in a degree of unconservatism when considering total impulse) whereas the Gaussian function (Figure 7(c)) can be seen to closely match the CFD data for all angles of incidence. The Gaussian model also exhibits the best quantitative agreement when considering both MAE and RMSE and therefore has been selected to represent the spatial distribution of specific impulse in the predictive methodology proposed in this article.
The Gaussian function was subsequently assessed against the full dataset, Figure 8, again with A G determined through non-linear regression ( A G = 2007 for the full dataset). Although the spread of normalised CFD curves is larger for the full dataset, the Gaussian function can again be seen to represent the mean CFD data to a high degree of accuracy and is therefore suitable for use in the entire range of near-field scaled distances, 0.11 ≤ Z ≤ 0.55 m/kg 1/3 .

Variation with scaled distance and final models
With a method for predicting distribution of normalised peak specific established, it still remains to quantify the variation of peak normally reflected specific impulse with scaled distance. Rather than making use of the existing higher order polynomial expressions of Kingery and Bulmash (1984), Swisdak (1994), or Shin et al. (2014b), the aim here is to develop a simple equation with a small number of coefficients over a more focussed range of scaled distances. Accordingly, the relationship between scaled distance and peak normally reflected specific impulse (i.e. the peak specific impulse at θ = 0  for each scaled distance) is plotted in Figure 9. Here, it can be seen that a power law represents the data well, with statistically robust fits of R 2 = 0.999 and R 2 = 0.995 for the reduced and full datasets respectively. Accordingly, the derived relationship for specific impulse, scaled using Hopkinson (1915) and Cranz (1926)

Evaluation against numerical dataset
This section presents a rigorous assessment of the quality of the proposed predictive method by comparing the results against the full CFD dataset. Equations (5) and (6) were solved for 0 80   ≤ ≤ θ between the limits of 0.11 ≤ Z ≤ 0.21 m/kg 1/3 and 0.11 ≤ Z ≤ 0.55 m/kg 1/3 respectively, in order to evaluate the accuracy of the predicted specific impulses against the entire CFD dataset. Filled contours of scaled peak specific impulse are shown for the reduced dataset in Figure 10(a), and the full dataset in Figure 10(c), with associated residuals shown in Figure 10(b) and (d).
It can be seen that the residuals are generally positive, particularly for the reduced dataset. This is acceptable given the slight conservatism shown in the CFD data itself (see Figure 4). The predictions for the full dataset show slight under-predictions for small scaled distances and θ < 30  , which is deemed acceptable given the availability of the reduced dataset predictions which encompass a smaller range of scaled distances but at higher accuracy. The maximum absolute residuals are 0.922 MPa·ms/kg 1/3 for the reduced dataset, and 1.944 MPa·ms/kg 1/3 for the full dataset, and the mean absolute residuals are 0.198 MPa·ms/kg 1/3 for the reduced dataset, and 0.184 MPa·ms/ kg 1/3 for the full dataset. The residual errors are considerably smaller than the maximum specific impulse (25 MPa·ms/kg 1/3 ), and indicate that the predictive equations are representative of the entire CFD dataset.
Whilst knowledge of the distributed specific impulse is a key requirement for determining structural response (Rigby et al., 2019b), total impulse (integral of specific impulse with respect to area, often termed 'area-integrated' impulse) is a useful metric particularly when considering  (5) and (6) were numerically integrated with respect to area for scaled target radii of 0.10, 0.25, 0.50 and 1.00 m/kg 1/3 , and for scaled distances of 0.11, 0.20, 0.30, 0.40 and 0.50 m/kg 1/3 . These values are provided in Table 2, alongside area-integrated CFD data and percentage differences between the two methods. For the two cases where Z = 0.11 1 3 m / kg / and Z = 0.20 1 3 m / kg / , predictions were generated for both equations (5) and (6) in order to comment on the relative accuracy of the two approaches.
Generally the predictions are in excellent agreement with the CFD results, particularly for Z ≤ 0.20 m/kg 1/3 when using equation (5), where the maximum deviation is 4% from the CFD data and the typical deviation is ± 1%. When the relationship in equation (6) is used to calculate total impulse for Z ≤ 0.20 m/kg 1/3 the results are less accurate, with a maximum deviation of 12% from the CFD data. However, when equation (6) is used to make predictions for Z > 0.20 1 3 m / kg / a similar level of accuracy to equation (5) is obtained, with a maximum deviation of 6% from the CFD data, and a typical deviation of ± 4%. Thus, the predictive methods are considered to be highly representative of the CFD data from which they are derived. (c) (d) Figure 10. Predicted specific impulse distributions: (a) filled contours of scaled peak specific impulse for reduced dataset, 0.11-0.21 m/kg 1/3 , (b) residual errors for reduced dataset, (c) filled contours of scaled peak specific impulse for full dataset, 0.11-0.55 m/kg 1/3 and (d) residual errors for full dataset.

Comparison against unseen CFD data
The predictive model was compared against two sets of unseen CFD data (i.e. data not used in the development of the predictive model). The first set comprises two models: a 5 kg spherical PE4 charge detonated at 0.291 m normal stand-off distance (to charge centre) from a rigid surface ( Z = 0.17 1 3 m / kg / ); and a 250 kg spherical PE4 charge detonated at 2.52 m normal stand-off distance (to charge centre) from a rigid surface ( Z = 0.40 1 3 m / kg / ). Note that these models were not used to form the CFD datasets outlined previously, and hence serve as unseen data for comparative purposes. In the first model, a 2.1 2.1 2.1 × × m domain was discretised into 0.1 m zone lengths with resolution level 4, giving an ultimate cell length of S / 46 . In the second model, a 19.0 19.0 19.0 × × m domain was discretised into 1.0 m zone lengths with resolution level 5, giving an ultimate cell length of S / 80 . Both meshes satisfy the limit of ultimate cell length not exceeding S /16 , as determined in the mesh sensitivity study. Figure 11 shows a comparison of the CFD and predicted specific impulse distributions out to 0.25 m and 2.50 m from the target centre, with select values presented in Table 3 alongside percentage differences between CFD and predicted values. Here, predictions for the Z = 0.17 1 3 m / kg / case were determined using equation (5), and predictions for the Z = 0.40 1 3 m / kg / case were determined using equation (6). It can be seen that the predicted distributions closely match those from the CFD analyses, albeit with the predictive model being slightly conservative by design. The quantitative agreement is excellent throughout: the maximum deviation between CFD and  m / kg / case (θ = 36.7  ). This exercise also validates the use of Hopkinson (1915) and Cranz (1926) scaling.
Detonation of explosive charges aligned with the centre of a target plate presents the worst-case scenario in terms of imparted loading and resulting plate deformation/failure (Chung Kim Yuen & Nurick 2005;Chung Kim Yuen et al., 2016). However, there may be situations where an off-centre detonation is of interest to the analyst, and therefore the purpose of the following exercise is to demonstrate generality of the predictive model in terms of charge placement.
According to the terminology in Figure 12(a), angle of incidence at any point on the target plate can be calculated knowing the lateral distance from the point of interest to the normal impingement point, r , and the normal stand-off distance, S : Where, in the previous examples described in this article, the normal impingement point was coincident with the plate centre, thus r = = 0 θ at that point. The second set of unseen CFD comprises one model: a 5 kg spherical PE4 charge detonated at 291 mm normal stand-off distance (normal scaled distance, Z = 0.17 1 3 m / kg / ) from a 200 mm square plate, with the charge offset from the plate centre by 50 mm in both in-plane directions, as in Figure 12(b). Here, the mesh specifications were the same as those in the previous 5 kg model. Noting that the plate centre is labelled as (0,0) , and the normal impingement point is labelled as (50,50) , r = 70.7 mm and θ = 13.7  at the plate centre. Figure 13 shows a comparison between the CFD results and Gaussian model predictions using the Gaussian model for the off-centre charge example, evaluated on a grid of 1 1 × mm elements. It can be seen that a high level of qualitative agreement is attained, which demonstrates the suitability of the method detailed herein for predicting the loading from non-centrally aligned explosive charges.

Comparison against unseen experimental data
Finally, the experimental data presented in Table 2 of Geretto et al. (2015) is used to assess the accuracy of the method in calculating total impulse. Here spheres of PE4 explosive with masses ranging from 10 to 70 g were detonated at a constant stand-off distance of 100 mm (to charge centre, 0.243 ≤ Z ≤ 0.464 m/kg 1/3 ) from a target plate with exposed area of 200 200 × mm (70.5° angle of incidence to the corner of the plate). Thirty-four experiments were performed altogether in eight different configurations. Total impulse was measured using a horizontal ballistic pendulum. Figure  14 shows the experimental results (individual test data and the mean at each scaled distance) and predictions from numerical integration of equation (6)   of 0.1 0.1 × mm elements. The predictions can be seen to closely match the overall trend of the experimental results. In particular, the predictions seem to better match the upper bound of the experimental data, which again is to be expected on account of the inherent conservatism of the CFD analyses and predictive method.
A direct comparison is provided in Table 4. It can be seen that, with the exception of the furthest scaled distance (16% difference), the predictions are all within 9% of the maximum recorded experimental area-integrated impulse. In all instances, the predictions are greater than the experimental values, confirming the slight conservatism of the method presented in this article. It is tentatively suggested that, in cases where a conservative model is not required, a factor of 0.9 may be applied to predicted specific impulses. Applying this factor here would shift the predictions from closely matching the maximum experimental values to closely matching the experimental mean.

Summary and conclusions
Peak specific impulse is the primary parameter that governs structural deformation under shortduration loading (Rigby et al., 2019b). Predicting the distribution of specific impulse acting on a reflecting surface following detonation of a high explosive charge remains a key challenge to engineers involved in transport security, infrastructure assessment and defence. This is a particular requirement in regions extremely close to the charge, where a historic lack of experimental data has precluded both the development of accurate semi-empirical predictive approaches and validation of high fidelity, physics-based numerical schemes. In this article, a new data-driven predictive model is developed to address a key issue facing the blast protection engineering community: the lack of an accurate and fast-running engineering tool for blast load prediction in the extreme near-field.
Newly-available experimental near-field specific impulse data is used to validate APOLLO Blastsimulator, an explicit CFD software specialising in problems involving blast wave propagation. A CFD dataset is populated with Apollo analyses, and the results are used to derive the novel data-driven surrogate model presented in this article. The approach is formed of two aspects: firstly, a Gaussian function is used to describe the relationship between normalised specific impulse and angle of incidence (shown in this article to form a consistent grouping for the range of scaled distances considered). Secondly, a power law is fit to values of peak normally reflected specific Table 4. Comparison of experimental area-integrated impulses (Geretto et al., 2015) and predictions. impulse, and Hopkinson-Cranz scaling (Cranz, 1926;Hopkinson, 1915) is used to express the predictions in scaled form. Separate equations are derived for the full dataset (0.11 ≤ Z ≤ 0.55 m/kg 1/3 ) and a reduced dataset ( 0.11 ≤ Z ≤ 0.21 m/kg 1/3 ) offering the analyst improved accuracy for extreme close-in detonations. Finally, accuracy of the approach is rigorously assessed against: the full CFD dataset; two new CFD analyses using 5 kg and 250 kg explosives (which also demonstrates the applicability of the implemented scaling laws and the generality of the approach for off-centre charge placement); and experimentally-recorded impulse acting on a 200 200 × mm plate from a range of explosives (10-70 g) placed at 100 mm stand-off (Geretto et al., 2015). The spatial distribution of impulsive loading and the total impulsive load are shown to be in good agreement throughout, with the predictive method demonstrating slight conservatism and variations ranging between 1% and 16% of the experimental/test examples, with a typical variation of <5%.
The simple equations presented herein (three coefficients and three input variables: scaled distance, angle of incidence and charge mass) enable the fast and accurate prediction of near-field blast loads in situations where previous methods are known to be unsuitable or inaccurate. This paper has demonstrated the applicability of data-driven predictive approaches for blast load prediction. The methodology outlined in this article has the potential to be supplemented with more advanced, machine learning-assisted approaches to model a more complex suite of scenarios and input parameters, such as explosive type (TNT equivalence), charge shape effects and blast wave clearing.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Jordan J Pannell gratefully acknowledges the financial support from the Engineering and Physical Sciences Research Council (EPSRC) Doctoral Training Partnership.