Research Article  Open Access
Analysis of Blast Wave Parameters Depending on Air Mesh Size
Abstract
Results of numerical simulations of explosion events greatly depend on the mesh size. Since these simulations demand large amounts of processing time, it is necessary to identify an optimal mesh size that will speed up the calculation and give adequate results. To obtain optimal mesh sizes for further largescale numerical simulations of blast wave interactions with overpasses, mesh size convergence tests were conducted for incident and reflected blast waves for close range bursts (up to 5 m). Ansys Autodyn hydrocode software was used for blast modelling in axisymmetric environment for incident pressures and in a 3D environment for reflected pressures. In the axisymmetric environment only the blast wave propagation through the air was considered, and in 3D environment blast wave interaction and reflection of a rigid surface were considered. Analysis showed that numerical results greatly depend on the mesh size and Richardson extrapolation was used for extrapolating optimal mesh size for considered blast scenarios.
1. Introduction
Today's global society faces many challenges that have not been present before. Not only are military facilities exposed to terrorist threats, but also many civilian buildings and structures are in danger. Most of these potential targets cannot be protected since they are in public use and have to be accessible. This paper is part of a wider research of the impact of explosions on overpasses that cross over the public roads [1, 2]. The imagined scenario assumes that the vehicle, in which explosive is placed, is parked underneath the overpass and then detonated. Since it is impossible to carry out the tests in the 1:1 scale and handling explosive is very dangerous and requires special conditions, the most practical tool is computer simulation. In both cases, carrying out field tests or computer simulations, it is very demanding task. Due to the very nature of the explosion process and many parameters that affect measured effects, even the field tests cannot provide accurate and repeatable results. On the other hand, in the computer simulation some assumptions and idealizations must be adopted. Very important issues in these simulations are hourglass effect, material properties, fragmentation of material, and fireball effect [3] but the most important task is to determine adequate mesh size. Coarse mesh may result in useless results, but too fine mesh may occupy computer resources for hours and often days in case of complex problems. The purpose of this paper is to offer guidelines for selection of appropriate air mesh size. Two scenarios were considered for the purpose of further clarifying the influence of air mesh size on simulation results and enabling further development of blast loading analysis. The first scenario is ideal freeair explosion where incident pressures were calculated and the second is specific blast scenario where reflected pressures from the overpass superstructure were calculated. For this latter case, the convergence tests were carried out.
Several researchers have studied the mesh size effect, individually or as a part of broader research on using numerical tools, to study blast wave propagation and interaction with structures. Chapman et al. [4] conducted a comprehensive parametric study of model grid size to obtain optimal mesh size as a part of their research on the use of Autodyn 2D for simulation of blast wave interaction with structures. Their results were compared to those from existing simple experimental tests, and they concluded that Autodyn 2D is a suitable tool for blast wave investigation. Luccioni et al. [5] studied influence of mesh size on blast wave propagation as a part of their research on urban blast. They concluded that a finer air mesh (10 cm) gives better results, but this can be too expensive in the terms of computational time. Alternatively, a coarser air mesh (50 cm) may be used to obtain qualitative results for comparison. Shi et al. [6] studied the relationship between the scaled distance and blast load parameters in a series of numerical simulations with different mesh sizes and compared the results with those from TM51300 [7]. They found that the pressure and impulse are less sensitive to mesh size with an increase of the scaled distance; they also proposed a method for the adjustment of peak pressures for coarse meshes. Nam et al. [3] determined the maximum element size to ensure that analytical results are independent of mesh size; this eliminates load gradient differences and reduces the gap between explosive and internal energy. Deng and Jin [8] conducted a parametric study of freeair blast wave propagation as a part of their research on bridge damage induced by blast loading. Simulation results indicated that a 5 cm air mesh size guarantees effective blast loading and reliable computational results. Wilkinson et al. [9] in their study on modelling closein detonations in the near field compared numerical results of overpressure obtained by Autodyn and Air3D. Mesh sensitivity analysis was conducted and results clearly indicated that mesh dependency of calculated results is existent. Shin et al. [10–12] performed numerical 1D, 2D, and 3D simulation of closein detonation of high explosives in order to validate CFD code for determining incident and reflected overpressure and impulses. It was concluded that mesh refinement studies are needed to establish cell sizes for air blast calculations and that a cell size of R/500 should be sufficient to predict overpressure histories in the near field, where R is the distance from the centre of the charge to the monitoring location. In addition, it was determined that the afterburning has a little effect on peak overpressure but can increase impulse if sufficient oxygen is present. Shin et al. [12–14] as a part of their research on blast parameters of detonation of spherical high explosives in freeair performed mesh size analysis on a wide range of scaled distances. Analysis produced recommendation of new polynomials and design charts for detonation parameters of spherical TNT charges. It was concluded that cell size of R/500 is sufficient for calculating pressures in the scaled distance range of 8 m/kg^{1/3}, similar to previous research [10].
Referred numerical simulations show that the accuracy of numerical results is strongly dependent on mesh size and blast scenario—a mesh size acceptable for one blast scenario might not be suitable for another. This is a clear indicator that a numerical mesh size convergence test needs to be carried out prior to largescale numerical simulations.
With respect to these results and recommendations, for the purpose of this research, it was necessary to determine the amount of explosive which is likely to be placed under the overpass. Table 1 gives crude estimations of possible TNT quantities that can be placed in common vehicles and used for explosive deployment under a bridge [21, 22]. The first three quantities are assumed to be the most probable and may offer some hope of bridge survival; detonation of larger TNT quantities would for sure cause imminent bridge collapse. Physical characteristics of vehicles, which dictate the shape, and confinement of the charge, detonation point, and consequent flying debris (fragmentation) are not taken into consideration, and only transport capacity was used in order to estimate the amount of explosive.
 
Note: 1 pound = 0.45 kg. 
2. Blast Wave Parameters
Explosion is an event of sudden and rapid energy release accompanied by a large pressure and temperature gradient on the wave front. It is very important to notice that blast parameters for freeair explosion are much more different from parameters which dictate the influence of explosion on engineering structure, in this case the overpass.
The blast wave is described by several parameters, but the main one is scaled distance, Z, which depends on charge weight, W [kg], and charge distance to target, R [m]. Scaled distance is calculated by the following expression:
Use of TNT as a referent explosive for determining Z is widely accepted. The first step in quantifying the blast waves from a source other than TNT is to transform the charge mass into the equivalent mass of TNT. In general, the TNT equivalent represents the mass of TNT that would result in an explosion with the same effect as the unit weight of explosive under consideration. It is defined as the ratio of the mass of TNT to the mass of the explosive resulting in the same magnitude of blast wave (or impulse, pressure, and energy) at the same radial distance for each charge, which assumes the scaling laws of Sachs and Hopkins. Common way to make conversion is to multiply the explosive mass by the conversion factor determined using one of several existing methods (the specific energy based concept, pressure based concept, impulse based concept, ChapmanJouguet state based concept, etc.).
The primary reason for choosing TNT as the referent explosive is that most of the available experimental data regarding the characteristics of blast waves were collected using TNT. Several methods exist for determining the explosive characteristics of different explosives, but they do not yield the same values for the TNT equivalent. These values depend on the characteristic parameter of the blast wave, the geometry of the load, and the distance from the explosive charge [23]. The TNT equivalence of terroristmanufactured explosive material is difficult to define precisely because of the variability of its formulation and the quality control of its manufacture. Table 2 gives the TNT equivalents for the most commonly used explosives [24].
 
Note: 1 pound = 0.45 kg. 
Pressuretime history for a blast wave is commonly described by the Friedlander equation [24]:where represents the blast wave overpressure, b is the waveform parameter, t_{0} is the positive phase duration, and t is considered time. Peak overpressure is the maximum blast pressure above normal atmospheric pressure (p_{0} = 101.325 kPa), which is dependent on the distance from the device. It can be calculated by numerous analytical formulations developed by researchers [25–27] who were able to carry out experiments on blast wave propagation. Differences in results exist in peak overpressures for near and farfield explosions, and they are incorporated in these formulations.
As the blast wave travels away from the detonation point, it undergoes reflection when the forward moving air molecules are brought to rest and further compressed upon meeting an obstacle. Combining peak overpressure and dynamic pressure, Rankine and Hugoniot derived the following equation for reflected overpressure:where represents the atmospheric pressure.
Further important parameters include specific impulse, time of arrival, and duration of positive and negative phase. Specific impulse represents the area beneath the pressuretime curve from arrival time to the end of the positive phase. Specific impulse is given by the following equation:where represents blast wave time of arrival, and is the duration of the positive phase. These specific blast times as well as other blast parameters can be obtained from US Army manual UFC 334002 [7].
3. Numerical Simulations
Numerical models developed for this investigation were created in Ansys Autodyn hydrocode software [28], which is intended for fluid dynamic analysis. Autodyn is a fully threedimensional explicit finite difference program which can utilize several different numerical techniques (Eulerian, Lagrangian, ALE, SPH, and Shell) to optimize the analysis of nonlinear dynamic problems. The structural analysis requires the ability to simulate both fluid/gas behaviour and structural response.
3.1. Numerical Techniques
The Lagrange processor is used for modelling solid continua and structures and operates on a structured (IJK) numerical mesh of quadrilateral (2D) or bricktype elements (3D) [29]. The numerical mesh moves and distorts with the material motion; no transport of material occurs from cell to cell. The advantage of such a scheme is that the motion of material is tracked very accurately, and the material interface and free surfaces are clearly defined. The primary disadvantage of the Lagrange formulation is that for severe material deformations or flow the numerical mesh will also become highly distorted, with attendant loss of accuracy and efficiency or outright failure of calculation [30].
The Euler processor is used for modelling fluids, gases, and large distortions. These processors include firstorder and secondorder accuracy schemes. A control volume method is used to solve the equations that govern conservation of mass, momentum, and energy. Numerical mesh is fixed in space, and material flows through it. The advantage of such a scheme is that large material flows and distortions can be easily treated. The disadvantage is that material interfaces and free surfaces are not naturally calculated, and sophisticated techniques must be utilized to track material interfaces. Additionally, historydependent material behaviour is more difficult to track.
ALE (Arbitrary Lagrange Euler) is a hybrid processor wherein the numerical mesh moves and distorts according to user specification. An additional computational step is employed to move the grid and remap the solution onto a new grid. ALE is an extension of the Lagrange method and combines the best features of both Lagrange and Euler methods. Due to advantages of both Lagrange and Eulerian solver, ALE solver was used for model calculation.
Modelling with Ansys using explicit time integration is limited by the CourantFriedrichsLevy condition [31]. This condition implies that the time step is limited so that a disturbance (e.g., stress wave) cannot travel further than the smallest characteristic element dimension in the mesh in a single time step. Thus the time step criterion for solution stability iswhere dt is the time increment, dx and dy are characteristic dimensions of an element, and c is the local material sound speed in an element. Based on this condition one can conclude that smaller mesh size implies smaller time steps and consequently longer calculation time.
3.2. Material Models
The numerical models used in this study include three different materials: environment, explosive charge, and overpass superstructure. Environment was modelled as the air, which transmits pressure waves, presented as an ideal gas, and the pressure was related to energy using following expression:where ρ is air density, e is specific internal energy, and γ is a constant (i.e., ratio of specific heats). This expression is one of the simplest forms of equation of state for gases. Gaseous materials have no ability to support any kind of stresses, and consequently no strength or failure relations are associated with this material.
Explosive charge was modelled as a spherical TNT charge. Similar to the air model, no strength or failure relations are associated with explosive material. Detonation is a rapid process that converts explosive material into gaseous products, and it is usually completed at the start of a given simulation. TNT explosive was modelled with the JonesWilkinsLee equation of state [32] using following expression:where p is hydrostatic pressure, χ is specific volume (1/ρ), and e is specific internal energy, and A, R_{1}, B, R_{2}, and ω are constants experimentally determined for several common explosives [33]. The shape of TNT charge in this study is spherical in order to ensure uniform propagation of the pressure waves and thus the clarity of results. Material properties of TNT and air are given in Table 3.
(a)  
 
Note: 1 inch = 2.54 cm; 1 pound = 0.45 kg; 1 psi = 6.89 kPa; 1 K = 272.15°C; †JWL = JonesWilkinsLee.  
(b)  
 
Note: 1 inch = 2.54 cm; 1 pound = 453.59 g; 1 psi = 6.89 kPa; 1 K = 272.15°C. 
In order to reduce the calculation time and get clear results, rigid material was chosen for the slab. This is a special feature of Autodyn software which may be used if the stresses and/or deformations of target object are not of interest. Insensitivity of the results on the rigid slab mesh size was confirmed by conducting several convergence tests and therefore is eliminated as an influential parameter in determining the dependency of blast wave parameters on air mesh size. Ground surface on which TNT charge was placed in 3D simulation before detonation was modelled as rigid boundary in order to provide full reflection of the blast wave and also eliminate ground mesh size as an influential parameter. Axisymmetric freeair simulations were performed using Eulerian solver, and 3D simulations were done using a coupled Lagrangian/Eulerian solver as it allows the combination of the best aspects of Lagrangian and Eulerian solvers.
3.3. Parametric Study
Modelling process of the impact of the explosion on overpass structure is divided into two parts. First, 2D FE model of freeair explosion is created. This model is composed of circular TNT charge, surrounding air, and outflow boundary, each with its properties (see Table 3) (freeair axisymmetric model). Radius of the TNT charge is determined by the specific reference density and chosen weight while radius of the air surrounding may be arbitrary defined (in this case 500 cm in order to shorten the calculation time). Using axial symmetry, only one circle sector may be modelled (Figure 2) and calculated which greatly reduces the calculation time. The results from this FE model are used to calculate freeair incident pressures. Second part is generation of the 3D model by revolving a 2D section 360° degrees and remapping the incident pressures to the overpass structure and rigid ground surface in order to compute reflected pressures (Figure 3).
Since the explosive charge is part of the model, it should be investigated whether its own mesh size (different from air mesh size) affects the results of the simulation. It was considered if the refinements of the charge mesh influence the level of energy release during detonation, time of blast wave arrival, or incident pressures, especially in 3D simulations, and should be taken into consideration. Due to this concern, additional simulations were conducted in order to study and determine the influence of charge mesh size. Simulations were conducted by varying air and charge mesh size for each of the three chosen quantities of explosive. Three charge weights, as mentioned earlier, 115 kg, 230 kg, and 680 kg of TNT (different line type), three air mesh sizes (different point type), and eight charge mesh sizes were considered. Numerical model is axisymmetric as shown in Figure 2, consisting of two separate parts, TNT charge and air surrounding. Simulation results are shown in Figure 1. Figure 1(a) shows normalized pressures while Figure 1(b) shows normalized total blast energy. Table 4 shows the values used for normalization, for each of the three air mesh sizes and three TNT quantities. These values were obtained in separate simulation.
 
Note: 1 inch = 25.4 mm; 1 pound = 0.45 kg; 1 psi = 6.89 kPa. 
(a) Normalized total blast energy versus charge mesh size
(b) Normalized pressure versus charge mesh size
Normalization is conducted in regard to the numerical model with the same charge and air mesh size; i.e., referent models were those with 5 mm, 10 mm, and 20 mm as both charge and air mesh sizes. Comparison of normalized total energies indicates that differences are minimal for all three observed parameters (under 5%), with somewhat larger but still acceptable differences in normalized pressures (under 10%). Due to minimal differences, charge meshes to coincide with air mesh in every simulation. This conclusion also has an impact on accelerating the modelling process because there is no need to mesh separately air and charge part of the model.
Freeair axisymmetric numerical simulation of explosion was performed using twelve different air mesh sizes (Table 5). Here, only the TNT in the air environment was modelled, and the wave propagation was observed. On the edge of the air environment, an outflow boundary condition was implemented to avoid the blast wave reflection and amplification; that is, to allow free gas outflow from designated air environment. Observed air environment is 5 m radius, with gauge points placed at equal distances (30.5 cm) in order to compare results for different scaled distances; air environment and gauge position are shown in Figure 2.
 
Note: 1 inch = 25.4 mm. 
The 3D model binds the air environment with a rigid slab, which simulates bridge superstructure and ground level for blast wave reflection and amplification. The Eulerian multimaterial solver was used for accurate simulation of hot gas expansion and interaction with structure primarily because there are no grid distortions during fluid flow through cells, and it allows large deformations, mixing of initially separate materials, and larger time steps. Its use requires more computations per cycle, finer zoning for similar accuracy, and extra cells for potential flow regions.
Parametric 3D numerical simulation of a ground explosion was performed using nine different air mesh sizes (Table 6). The mesh size of 800 mm was discarded because it was too coarse and gave unrealistically small pressures in comparison to ATBlast [34]. Mesh sizes smaller than 3.13 mm were impractical because of high requirements for computer hardware components. The ground explosion model and gauge scheme are shown in Figure 3. The ideally rigid concrete slab was placed 5 m above the detonation point on the ground, which is typical clearing underneath the overpass. In addition, using the symmetry, only onequarter of the numerical model presented in Figure 3 was modelled and simulated. Additional shortening of calculation time was achieved by remapping the final solution of axisymmetric analysis for selected charge quantities (115 kg, 230 kg, and 680 kg) as an input (load) to the 3D environment [35].
 
Note: 1 inch = 25.4 mm. 
4. Comparison of Results and Discussion
Numerical simulation results using coarse and fine air mesh sizes varied considerably. This indicates that mesh size convergence tests are needed to obtain the optimal mesh size for an accurate simulation of blast wave propagation. Detonation products influence blast wave pressures and consequently impulses; that is, they are within the blast fireball. The pressures depend on the type of explosive as each explosive has different detonation products. The exact influence of the fireball on the pressures could not be determined separately because additional analysis would be needed in which data are collected outside the fireball radius and compared to those obtained within the fireball.
4.1. Incident Pressures
Results obtained by freeair axisymmetric numerical simulation were compared with results acquired by ATBlast, the software for analytical determination of the referent blast parameters. ATBlast uses equations for blast estimation of overpressures at different ranges which were developed by Kingery and Bulmash [36]. These equations are widely accepted as engineering predictions for determining freefield pressures and loads on structures. Kingery and Bulmash have compiled results of explosive tests in which they detonated charges weighting from less than 1 kg to over 400,000 kg and used curvefitting techniques to represent the data with highorder polynomial equations. Another reason why ATBlast was chosen as the reference for the blast parameter comparison was that the software is publicly available, subjected to ITAR (International Traffic in Arms Regulations) controls. Unlike ATBlast, ConWep [37], another widely used software for estimating blast properties, is not available to nonUS citizens or organizations; therefore, it was not used as an additional control for the numerical results.
It can be seen that decreasing the mesh size allowed maximum incident pressure to converge to the referent value. With smaller mesh sizes (from 0.39 to 12.5 mm) maximum incident pressures exceeded values obtained by ATBlast and can be accounted for by inaccuracies in the experimental pressure readings (the measuring instruments were influenced by very high pressures, and the blast wave duration and time of arrival were also very short, i.e., only a few milliseconds). The relative difference in maximum incident pressures measured at the last gauge point, 5 m from detonation point, between smaller mesh sizes and values generated by ATBlast can clearly be seen from Table 7. The closest incident peak pressure was obtained for the 25 mm mesh size (marked by bold in Table 7), with about 2% maximum difference in relation to the ATBlast value. Pressure histories for different mesh sizes relative to ATBlast pressure values are plotted in Figure 4. It can be seen that the rate of pressure amplification from ambient to peak pressure becomes slower and the shape of the curve becomes flatter with the increase of the mesh size. The calculated incident peak pressure was also smaller with the larger mesh sizes.
 
Note: 1 inch = 25.4 mm; 1 pound = 0.45 kg; 1 psi = 6.89 kPa. 
(a) 115 kg TNT (compact car trunk)
(b) 230 kg TNT (trunk of a large car)
(c) 680 kg TNT (closed van)
If the positive peak incident pressure in relation to scaled distance is considered, its value decreases with the increase of scaled distance. This is in accordance with the initial assumption that with the increase of distance from detonation point to structure surface pressures exerted on structure surfaces are decreasing. Figure 5 shows that this is consistent for all air mesh sizes. With smaller air mesh size, pressures became closer to the referent values obtained by ATBlast.
(a) 115 kg TNT (compact car trunk)
(b) 230 kg TNT (trunk of a large car)
(c) 680 kg TNT (closed van)
4.2. Reflected Pressures
Results obtained by 3D numerical simulation were also compared with the results acquired by ATBlast. In terms of absolute values, the comparison gave the closest reflected pressures for 25 mm mesh size in an analysis for 115 kg and 230 kg of TNT; here, the minimum difference in relation to ATBlast values was 8.57% for 115 kg and 0.08% for 230 kg. For 680 kg of TNT closest reflected pressure in comparison to ATBlast was obtained for 12.5 mm mesh size with difference of 4.17% (Table 8). Reflected overpressure histories in relation to ATBlast reflected overpressure values are plotted in Figure 6. Pressure amplification from ambient to peak reflected overpressure has a smaller rate of increase, and the shape of the curve becomes flatter with the increase of the mesh size as for the incident pressures. Pressuretime curve has additional pressure peak after the initial peak, which is caused by blast wave reflection and it can be clearly seen in Figure 6. Calculated reflected overpressure also becomes smaller with larger mesh sizes.
 
Note: 1 inch = 25.4 mm; 1 pound = 0.45 kg; 1 psi = 6.89 kPa. 
(a) 115 kg TNT (compact car trunk)
(b) 230 kg TNT (trunk of a large car)
(c) 680 kg TNT (closed van)
4.3. Time of Arrival
Mesh size also influenced the blast wave time of arrival. Some authors [6] have concluded that the arrival time is independent of the mesh size, and if we consider that the time in blast problems is measured in milliseconds, differences arising from analyses with different mesh sizes exist. Tendency is that the smaller mesh sizes shorten the time required for the blast wave front to arrive to the designated measuring point. If the mesh size is sufficiently small, it no longer influences the arrival time, as can be seen in Figure 7. For larger mesh sizes for which pressure history is flatter, the time of maximum incident pressure is taken as the time of blast wave arrival.
(a) Reflected wave
(b) Incident wave
If a difference of ±15% was considered, all arrival times were within this range. Results presented in Figure 7 show that the numerical analysis provided earlier arrival times in comparison to ATBlast. Similar to incident pressures, there was a tendency for smaller mesh sizes to shorten the arrival time; if the mesh size was sufficiently small, it no longer had an influence on the arrival time.
5. Result Convergence
Individual application of numerical methods differs in the level of accuracy with which a numerical model can reflect the physical phenomena. An analysis with the help of complex numerical models is particularly applicable if there is no closed analytical solution and if the experiments are impossible or insufficient. In all such cases, the question arises about the ability of such analysis to correctly predict the results of the physical process in question. The first recommended procedure within the verification process is the mesh refinement study. Its main objective is to evaluate the error of discretization and to check if the applied mesh is sufficiently refined.
5.1. Richardson Extrapolation
The mesh refinement study is conducted based on a comparison of the results for at least two but usually three meshes [29]. The Richardson extrapolation serves as a higherorder estimate of the evaluated quantity. A quantity f calculated for a mesh characterized by parameter (mesh size) h can be expressed using Taylor’s theorem. In practice, the Richardson extrapolation is generalized for any, also noninteger, pth order approximations and the mesh refinement ratio r and is considered as p + 1 order approximation. The exact solution of extrapolation can be described as an asymptotic solution for a mesh with element dimension h approaching zero. Estimate of the asymptotic solution is given by where denotes asymptotic solution for h approaching 0 (extrapolated value); and secondorder approximations of ; r is refinement ratio and p is order of convergence. Refinement ratio, r, can be set as arbitrary when using unstructured meshes or as a constant if regular meshes are considered. There is no recommendation about refinement ratio selection. In this study refinement ratio is set as a constant value of 2 which means that the finer meshes considered are generated by dividing node spacing in all directions into halves. Order of convergence, p, can be calculated no matter what refinement ratio is selected, constant or arbitrary. For structured meshes order of convergence can be calculated using (9) and for unstructured meshes using (10):where and . If constant refinement ratio is implemented in equation for calculating order of convergence for unstructured meshes (implemented in (10)), equal value as for equation for structured mesh should be obtained. All parameters used in extrapolation and obtained results are listed in Tables 9 and 10.
 
Note: 1 pound = 0.45 kg; 1 psi = 6.89 kPa. 
 
Note: 1 pound = 0.45 kg; 1 psi = 6.89 kPa. 
Given in a percentage manner, the grid convergence index (GCI) can be considered as a relative error bound showing how the solution calculated for the finest mesh is far from the asymptotic value. It predicts how much the solution would change with further refinement of the mesh. The smaller the value of the GCI, the better. This indicates that the computed solution is within the asymptotic range. The safety factors are arbitrarily set based on the accumulated experience on computation fluid dynamics (CFD) calculations [38–40]. The safety factor represents 95% confidence for the estimated error bound. This assumption can be expressed as the following statement: There is 95% confidence that the converged solution is within the range [f_{1}(1 − GCI_{12}/100%), f_{1}(1 + GCI_{12}/100%)]. The GCI is defined aswhere is a safety factor. The recommended CFD values of the safety factor are
= 3.0 when two meshes are considered,
= 1.25 for three and more meshes.
Quantity ε defines relative difference between subsequent solutions:
Due to lack of experimental investigation and reliance on ATBlast values, Richardson extrapolation was introduced to determine convergence of the simulations. Pressure values () were obtained through numerical simulations, using different mesh sizes with a constant refinement ratio (r) of 2. With known solutions, order of convergence (p) was determined using (9), based on three subsequent pressure values. Extrapolated pressures ( and ) were calculated using pressures from two subsequent meshes ( and ), order of convergence (p), and refinement ratio (r). Table 9 lists the pressure extrapolated values for axisymmetric simulations and Table 10 for 3D simulations together with estimated error bounds, GCI, refrainment ratio, and order of convergence. Extrapolated values were obtained by using three smallest air mesh sizes, 1.56 mm, 0.78 mm, and 0.39 mm for axisymmetric and 12.5 mm, 6.25 mm, and 3.13 mm for 3D simulations, respectively, following the previously stated calculation procedure given by [41].
With the reduction of the mesh size, pressures converge to a value dependent on the charge mass, as can be seen in Figure 8. This fact should be taken with caution because smaller mesh sizes substantially prolong calculation time and demand large quantities of computational memory. Mesh sizes which produce pressure difference in comparison to extrapolated pressure less than 10% are considered sufficiently accurate for engineering use. The analysis was based on closein and intermediate scaled distances from 0 m/kg^{1/3} to 1 m/kg^{1/3} that coincide with the previously mentioned bridge superstructure distance from the road surface. If lower bound is considered, appropriate incident pressure values are obtained for 12.5 mm air mesh size for 115 kg, 230 kg, and 680 kg of TNT, respectively.
Gauge point of interest was located in the lower plane of bridge superstructure, at a distance of 5 m from ground plane, which is equivalent to scaled distances for each TNT quantity: 1.03 m/ for 115 kg, 0.82 m/ for 230 kg, and 0.57 m/ for 680 kg. Using Richardson extrapolation and accepting 10% upper and lower bound it was possible to obtain acceptable air mesh sizes for reflected pressures (Figure 9). Analysis resulted in air mesh sizes of 12.5 mm for 115 kg TNT and 25 mm for 230 kg and 680 kg TNT being adequate for pressurestructure interaction analysis. If it is considered that it is practically impossible to determine absolutely accurate blast pressures, using experimental measurements or numerical simulation, according to this analysis ATBlast could provide an adequate approximation of blast wave pressures for both incident and reflected pressures. ATBlast calculated pressures are very close or are within the value of lower bound (10%) of extrapolated pressures, clearly depicted in Figures 8 and 9, which is yet another argument for using ATBlast for quick assessment of blast pressures.
5.2. Empirical and Analytical Expressions
Empirical and analytical expressions are given in literature for estimating blast overpressure. Some of the expressions are based on the analysis of experimentally acquired pressure data while other are product of purely mathematical algorithms. Expressions are given for calculating overpressures generated from spherical or hemispherical airbursts. Based on the parameters used for numerical simulations (stand of distance of 5 m and charge weight of 115 kg, 230 kg, and 640 kg) calculation of analytical and empirical expressions for incident overpressure estimation by different authors was conducted. Calculated results were used for numerical and ATBlast overpressure evaluation. The list of analytical expressions was given by Ullah et al. (2017) [42]. Table 11 gives some perspective on calculated values and differences of expressions with ATBlast and FEM simulations. Based on the results given in Table 11 calculations provide a wide variety of overpressure values with differences ranging from 0.1% to 166%.
 
Note: 1 ft = 0.3048 m; 1 pound = 0.45 kg; 1 psi = 6.89 kPa. 
If results from Table 11 are analysed it can be seen that for some equations calculated overpressure is higher than from ATBlast, and in other cases overpressure is lower than from ATBlast. The same can be said for FEM results. Exact value of overpressures (p_{s}) for analysed charge quantities is obtained by Swisdak’s equation (see (13)) but that is because it is based on the KingeryBulmash equation and parameters, and ATBlast internal equations are exactly those defined by KingeryBulmash:where Z is the scaled distance and A_{1}, B_{1}, C_{1}, D_{1}, E_{1}, F_{1}, and G_{1} are simplified Kingery air blast coefficients given in Table 12.

Due to large discrepancy between overpressures calculated using different expressions and ATBlast (and numerical) results it is hard to evaluate overall performance of ATBlast. Expressions are mathematical equations obtained through correlation with experimental data and most of these approaches are limited by extent of overall experimental data set. The tendency of empirical or analytical expressions is that the range of the charge (scaled distance) is decreasing so is the accuracy reduced. This is due to uncertainties in measured experimental data in near field and contact explosions.
6. Validation of Numerical Models
Validation of numerical models was conducted based on comparison of numerical simulation results and experimental measurements. Purpose of this comparison is to establish the degree of error of parametrically determined acceptable air mesh size in obtaining field measured blast pressures. Results are compared, not only with measured but also with estimated pressures given in research listed in Tables 13 and 14. Researchers estimated pressures using CONWEP and/or analytical calculation. Incident pressures were measured using freefield pressure gauges positioned on some distance from explosive charge while reflected pressures were measured by pressure gauges positioned in such a way to reduce influence of material deformation and deterioration to pressure values. Sensors were placed on a rigid containment structure that provided planed support and boundary conditions, for test specimens during experiment. In experimental testing performed by [15, 20] sensors were positioned on the surface of tested specimen and they were not able to measure pressures in all tests. From eight sensors, positioned by [20] just three to four were able to successfully measure blast pressures; the rest of the data was corrupted due to sensor malfunctioning. Table 13 provides comparison of measured and simulated incident pressures, while Table 14 gives comparison of reflected pressures. Pressures were additionally estimated by ATBlast in order to validate the software as a tool for quick pressure and impulse assessment.
 
Note: 1 kg = 0.453592 kg; 1 ft = 0.3048 m; 1 psi = 6.89 kPa; M: measured; E: estimated; FEM: finite element model; ATB: ATBlast; N/A: not available. 
 
Note: 1 kg = 0.453592 kg; 1 ft = 0.3048 m; 1 psi = 6.89 kPa; M: measured; E: estimated; FEM: finite element model; ATB: ATBlast. 
In Table 13 it can be seen that the numerical results overestimate measured incident pressures from 1% to 25%. From Table 14 it is impossible to draw general conclusion about reflected pressure variations; nevertheless, pressures are within acceptable error margin from which it can be concluded that numerical models provide sufficiently accurate results. If compared to estimated values by authors and ATBlast, numerical results provide smaller errors, in previously stated range. Estimated values are generally larger than those measured which provides pressure overestimation and accounts for all possible measurement errors. This pressure overestimation can lead to conservative structural design but from the aspect of human life protection in structures of higher importance this is acceptable.
7. Conclusion
No recommendations exist for mesh size selection in blast load numerical simulations since the size of the finite element strongly depends on the blast scenario. For this reason the mesh size convergence tests were carried out for ground explosion scenarios. Tests comprised axisymmetric and 3D blast numerical simulations in which different air mesh sizes were used to compare the results with ATBlast and extrapolated values. A 10% difference of blast parameters in relation to extrapolated values was considered sufficiently accurate for engineering use, and based on this assumption air mesh size recommendations were made for two different blast scenarios.
In both cases (freeair and reflected scenario), coarse mesh led to stretching of the pressuretime profile without a distinct arrival time of the blast front. If the rate of initial pressure rise is not very steep (close to 90°) but rather gradual, it is clear sign that the mesh size is too coarse. Another indication of inadequacy of mesh size is lack of clear secondary pressure peak (see Figure 6). Finer mesh shortened the time required for blast wave front to arrive to the structure, and if the mesh size was sufficiently small, it no longer influenced the arrival time. In the general case of freeair explosion, the largest mesh size, which is adequate to capture all the blast parameters, is 12.5 mm. For the specific blast scenario (detonation of the charge beneath the bridge structure), the largest adequate mesh sizes are 12.5 mm for 115 kg of TNT and 25 mm for 230 and 680 kg of TNT.
Results supported the initial presumption that, with the increase of distance between the detonation point and structure, pressures are decreasing but the sensitivity on the mesh size and scattering of the results are greater for smaller scaled distances. The pressure values for larger scaled distances were less sensitive to air mesh size.
It should be emphasized that this research was carried out for ideal conditions (spherical shape of charge, rigid bridge deck, and rigid ground). Further research should aim to investigate pressure variations in less ideal cases: boxlike charges, nonrigid concrete slab, compressible ground, and explosive confined in car trunk (with offset from ground level).
Data Availability
Authors state that the majority of data acquired from research and input parameters needed for numerical simulations (repeatability of simulations) is listed in tables implemented in the submitted article. If additional data is requested, authors will be more than glad to share it with interested researchers. Authors think that publicly available data contributes to research better quality; it is an opportunity for discussion and exchange of new ideas among researchers.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research presented in this paper is a part of the research project “Behaviour of Medium and Short Length Bridges Subjected to Blast Load” supported by the Faculty of Civil Engineering Osijek, and support for this research is gratefully acknowledged.
References
 H. Draganić, Overpasses Subjected to Blast Load, Faculty of Civil Engineering Osijek, Josip Juraj Strossmayer University of Osijek, Osijek, 2014. View at: Publisher Site
 H. Draganić and D. Varevac, “Numerical simulation of effect of explosive action on overpasses,” Gradjevinar, vol. 69, no. 6, pp. 437–451, 2017. View at: Google Scholar
 J. W. Nam, J. J. Kim, S. B. Kim, N. H. Yi, and K. J. Byun, “A study on mesh size dependency of finite element blast structural analysis induced by nonuniform pressure distribution from high explosive blast wave,” KSCE Journal of Civil Engineering, vol. 12, no. 4, pp. 259–265, 2008. View at: Publisher Site  Google Scholar
 T. C. Chapman, T. A. Rose, and P. D. Smith, “Blast wave simulation using AUTODYN2D: a parametric study,” International Journal of Impact Engineering, vol. 16, no. 56, pp. 777–787, 1995. View at: Publisher Site  Google Scholar
 B. Luccioni, D. Ambrosini, and R. Danesi, “Blast load assessment using hydrocodes,” Engineering Structures, vol. 28, no. 12, pp. 1736–1744, 2006. View at: Publisher Site  Google Scholar
 Y. Shi, Z. Li, and H. Hao, “Mesh size effect in numerical simulation of blast wave propagation and interaction with structures,” Transactions of Tianjin University, vol. 14, no. 6, pp. 396–402, 2008. View at: Publisher Site  Google Scholar
 UFC, Structures to Resist The Effects of Accidental Explosions, 2008.
 R.B. Deng and X.L. Jin, “Numerical Simulation of Bridge Damage under Blast Loads,” WSES Transactions on Computers, vol. 8, p. 10, 2009. View at: Google Scholar
 W. Wilkinson, D. Cormie, J. Shin, and A. Whittaker, “Modelling closein detonations in the nearfield,” in Proceedings of the 15th International Symposium on Interaction of the Effects of Munitions with Structures, Potsdam, Germany, 2013. View at: Google Scholar
 J. Shin, A. S. Whittaker, D. Cormie, and W. Wilkinson, “Numerical modeling of closein detonations of high explosives,” Engineering Structures, vol. 81, pp. 88–97, 2014. View at: Publisher Site  Google Scholar
 J. Shin, A. Whittaker, D. Cormie, and A. Aref, Verification and validation of a CFD code for modeling detonations of high explosives, 2015.
 J. Shin, Airblast effects on civil structures, State University of New York, Buffalo, 2014.
 J. Shin, A. S. Whittaker, and D. Cormie, “Incident and Normally Reflected Overpressure and Impulse for Detonations of Spherical High Explosives in Free Air,” Journal of Structural Engineering (United States), vol. 141, no. 12, 2015. View at: Google Scholar
 J. Shin, A. Whittaker, and D. Cormie, “Estimating incident and reflected airblast parameters: updated design charts,” in Proceedings of the 16th international symposium for the interaction of the effects of munitions with structures (ISIEMS 16), 2015. View at: Google Scholar
 N. Yi, J. J. Kim, T. Han, Y. Cho, and J. H. Lee, “Blastresistant characteristics of ultrahigh strength concrete and reactive powder concrete,” Construction and Building Materials, vol. 28, no. 1, pp. 694–707, 2012. View at: Publisher Site  Google Scholar
 C. Q. Wu, L. Huang, and D. J. Oehlers, “Blast testing of aluminum foamprotected reinforced concrete slabs,” Journal of Performance of Constructed Facilities, vol. 25, no. 5, pp. 464–474, 2011. View at: Publisher Site  Google Scholar
 A. Ghani Razaqpur, A. Tolba, and E. Contestabile, “Blast loading response of reinforced concrete panels reinforced with externally bonded GFRP laminates,” Composites Part B: Engineering, vol. 38, no. 56, pp. 535–546, 2007. View at: Publisher Site  Google Scholar
 T. S. Lok and J. R. Xiao, “Steelfibrereinforced concrete panels exposed to air blast loading,” Proceedings of the Institution of Civil Engineers: Structures and Buildings, vol. 134, no. 4, pp. 319–331, 1999. View at: Publisher Site  Google Scholar
 J. Li, C. Wu, and H. Hao, “Residual loading capacity of ultrahigh performance concrete columns after blast loads,” International Journal of Protective Structures, vol. 6, no. 4, pp. 649–669, 2015. View at: Publisher Site  Google Scholar
 L. Chen, Q. Fang, J. Fan, Y. Zhang, H. Hao, and J. Liu, “Responses of masonry infill walls retrofitted with CFRP, steel wire mesh and laminated bars to blast loadings,” Advances in Structural Engineering, vol. 17, no. 6, pp. 817–836, 2014. View at: Publisher Site  Google Scholar
 A. H. Hameed, “Dynamic Behaviour Of Reinforced Concrete Structures Subjected To External Explosion, A Thesis,” Jumada Elawal, vol. 1428, 2007. View at: Google Scholar
 J. S. Humphreys, “Plastic deformation of impulsively loaded straight clamped beams,” Journal of Applied Mechanics, vol. 32, no. 1, pp. 7–10, 1964. View at: Google Scholar
 A. Freidenberg, A. Aviram, L. K. Stewart, D. Whisler, H. Kim, and G. A. Hegemier, “Demonstration of tailored impact to achieve blastlike loading,” International Journal of Impact Engineering, vol. 71, pp. 97–105, 2014. View at: Publisher Site  Google Scholar
 G. C. Mays and P. D. Smith, “Blast effects on buildings. Design of buildings to optimize resistance to blast loading,” T. Telford, 1995. View at: Google Scholar
 N. Newmark and R. Hansen, Design of blast resistant structures, Shock and vibration Handbook, vol. 3, 1961, Design of blast resistant structures, Shock and vibration Handbook.
 J. Drake, L. Twisdale, R. Frank et al., “Protective Construction Design Manual: Resistance of Structural Elements (Section IX),” Report ESLTR8757, Air Force Engineering and Services Center, Tyndall Air Force, Base, 1989. View at: Google Scholar
 P. D. Smith, J. G. Hetherington, and P. Smith, Blast and ballistic loading of structures, ButterworthHeinemann Oxford, 1994.
 ANSYS, “ANSYS AUTODYN User manual,” in ANSYS Release 14.0, pp. 1–464, Canonsburg, PA, USA, 2010. View at: Google Scholar
 L. Kwasniewski, “Application of grid convergence index in FE computation,” Bulletin of the Polish Academy of Sciences—Technical Sciences, vol. 61, no. 1, pp. 123–128, 2013. View at: Google Scholar
 N. K. Bimbaum, J. Tancreto, and K. Hager, “Calculation of blast loading in the high performance magazine with AUTODYN3D,” DTIC Document, 1994. View at: Google Scholar
 R. Courant, K. Friedrichs, and H. Lewy, “Über die partiellen differenzengleichungen der mathematischen physik,” Mathematische Annalen, vol. 100, no. 1, pp. 32–74, 1928. View at: Publisher Site  Google Scholar  MathSciNet
 G. Baudin and R. Serradeill, “Review of JonesWilkinsLee equation of state,” in Proceedings of the EPJ Web of Conferences, vol. 10, 00021, 2010. View at: Google Scholar
 E. L. Lee, H. C. Hornig, and J. W. Kury, “Adiabatic Expansion Of High Explosive Detonation Products,” Tech. Rep. UCRL50422, California Univ., Livermore. Lawrence Radiation Lab., 1968. View at: Publisher Site  Google Scholar
 ARA, A.T.Blast, in, Applied Research Associates, Inc., Vicksburg, Mississippi, USA, 2013.
 ANSYS, “Remapping Tutorial Revision 4.3,” in Autodyn; Explicite Software for Nonlinear Dynamics, pp. 1–32, Canonsburg, PA, USA, 2005. View at: Google Scholar
 G. K. Schleyer, M. J. Lowak, M. A. Polcyn, and G. S. Langdon, “Experimental investigation of blast wall panels under shock pressure loading,” International Journal of Impact Engineering, vol. 34, no. 6, pp. 1095–1118, 2007. View at: Publisher Site  Google Scholar
 J.A. Calgaro, M. Tschumi, and H. Gulvanessian, Designers' Guide to Eurocode 1: Actions on Bridges, 2010.
 I. Celik and O. Karatekin, “Numerical experiments on application of richardson extrapolation with nonuniform grids,” Journal of Fluids Engineering, vol. 119, no. 3, pp. 584–590, 1997. View at: Publisher Site  Google Scholar
 J. Slater, Examining spatial (grid) convergence, Public tutorial on CFD verification and validation, vol. MS, 86, NASA Glenn Research Centre, 2006.
 L. E. Schwer, “Is your mesh refined enough? Estimating discretization error using GCI,” in LSDYNA Forum, Germany, Bamberg, 2008. View at: Google Scholar
 I. B. Celik, “Procedure for estimation and reporting of uncertainty due to discretization in CFD applications,” Journal of Fluids Engineering, vol. 130, no. 7, pp. 0780011–0780014, 2008. View at: Publisher Site  Google Scholar
 A. Ullah, F. Ahmad, H.W. Jang, S.W. Kim, and J.W. Hong, “Review of analytical and empirical estimations for incident blast pressure,” KSCE Journal of Civil Engineering, pp. 2211–2225, 2017. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Hrvoje Draganić and Damir Varevac. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.