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 large-scale 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 free-air 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 TM5-1300 [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 free-air 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 close-in 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 close-in 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 free-air 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/kg1/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 large-scale 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.
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 free-air 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, Chapman-Jouguet 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 terrorist-manufactured 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].
Pressure-time history for a blast wave is commonly described by the Friedlander equation [24]:where represents the blast wave overpressure, b is the waveform parameter, t0 is the positive phase duration, and t is considered time. Peak overpressure is the maximum blast pressure above normal atmospheric pressure (p0 = 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 far-field 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 pressure-time 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 3-340-02 [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 three-dimensional 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 (I-J-K) numerical mesh of quadrilateral (2D) or brick-type 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 first-order and second-order 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, history-dependent 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 Courant-Friedrichs-Levy 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 Jones-Wilkins-Lee equation of state [32] using following expression:where p is hydrostatic pressure, χ is specific volume (1/ρ), and e is specific internal energy, and A, R1, B, R2, 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.
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 free-air 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 free-air explosion is created. This model is composed of circular TNT charge, surrounding air, and outflow boundary, each with its properties (see Table 3) (free-air 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 free-air 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.
(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.
Free-air 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.
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 AT-Blast [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 one-quarter 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].
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 free-air axisymmetric numerical simulation were compared with results acquired by AT-Blast, the software for analytical determination of the referent blast parameters. AT-Blast 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 free-field 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 curve-fitting techniques to represent the data with high-order polynomial equations. Another reason why AT-Blast 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 AT-Blast, ConWep [37], another widely used software for estimating blast properties, is not available to non-US 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 AT-Blast 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 AT-Blast 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 AT-Blast value. Pressure histories for different mesh sizes relative to AT-Blast 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.
(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 AT-Blast.
(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 AT-Blast. 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 AT-Blast values was 8.57% for 115 kg and 0.08% for 230 kg. For 680 kg of TNT closest reflected pressure in comparison to AT-Blast was obtained for 12.5 mm mesh size with difference of 4.17% (Table 8). Reflected overpressure histories in relation to AT-Blast 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. Pressure-time 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.
(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 AT-Blast. 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 higher-order 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, p-th 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 second-order 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.
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 [f1(1 − GCI12/100%), f1(1 + GCI12/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 AT-Blast 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 close-in and intermediate scaled distances from 0 m/kg1/3 to 1 m/kg1/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 pressure-structure 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 AT-Blast could provide an adequate approximation of blast wave pressures for both incident and reflected pressures. AT-Blast 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 AT-Blast 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 AT-Blast 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 AT-Blast 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%.
If results from Table 11 are analysed it can be seen that for some equations calculated overpressure is higher than from AT-Blast, and in other cases overpressure is lower than from AT-Blast. The same can be said for FEM results. Exact value of overpressures (ps) for analysed charge quantities is obtained by Swisdak’s equation (see (13)) but that is because it is based on the Kingery-Bulmash equation and parameters, and AT-Blast internal equations are exactly those defined by Kingery-Bulmash:where Z is the scaled distance and A1, B1, C1, D1, E1, F1, and G1 are simplified Kingery air blast coefficients given in Table 12.
Due to large discrepancy between overpressures calculated using different expressions and AT-Blast (and numerical) results it is hard to evaluate overall performance of AT-Blast. 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 free-field 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 AT-Blast in order to validate the software as a tool for quick pressure and impulse assessment.
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 AT-Blast, 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 AT-Blast 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 (free-air and reflected scenario), coarse mesh led to stretching of the pressure-time 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 free-air 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.