# Replacing Detonation by Compressed Balloon Approaches in Finite Element Models

Authors: Pierre Legrand; S. Kerampran; M. Arrigoni
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

The evaluation of blast effects from malicious or accidental detonation of an explosive device is really challenging especially on large buildings. Indeed, the time and space scales of the explosion together with the chemical reactions and fluid mechanics make the numerical model really difficult to achieve acceptable structural design. Nevertheless, finite element methods and especially Arbitrary Lagrangian Eulerian (ALE) have been extensively used in the past few decades with some simplifications. Among them, the replacement of the explosive event by a compressed balloon of detonation products has been proven useful in numerous different situations. Unfortunately, the ALE algorithm does not achieve a proper energy balance through the numerical integration of the discrete scheme; this important drawback is not compensated by the use of the classical compressed balloon approach. The paper focuses on increasing the radius of the equivalent ideal gas balloon in order to achieve better energy balance and thus better results at later stages of the blast wave propagation.


---

## Body

## 1. Introduction

The evaluation of blast effects from the malicious or accidental detonation of an explosive device on structures is a multiphysics issue involving widely different time and space scales. Indeed, the detonation process of high explosives implies chemical reactions propagating within the solid explosive for a few microseconds. This reaction will create detonation products, in other words high pressure and temperature gases which will first react with themselves to create secondary reactions called afterburning. These gases then will expand and also react with the surrounding air and create a blast wave propagating outwards in the air for a few milliseconds, a possible fireball, and high pressure and temperature gradients. An extensive description of the detonation process can be found in [1]. Finally structural components will be impacted by the blast wave and may collapse over a few seconds. All those phenomena come with different space scales from the cm of the explosive sources to the dozens of meters of the structure.In practical use, when there is a need for shorter calculation times, cost reduction, and structural optimization, dealing with all phenomena is challenging, especially on a commonly available high performance computer (HPC). In order to speed up the calculation, a lot of effort has been made in the last 50 years to find suitable methods to evaluate the behavior of structures loaded by blast waves in a reasonable computational time. Among them, one of the most efficient methods consists in calculating the blast wave without taking the actual detonation process into consideration. In fact, for structural design, the pressure-time curve and specifically both overpressure and impulse transmitted to the structure are the most important parameters [2, 3]. In spite of neglecting the chemical detonation and so close-in effects (afterburning, fireball) the proposed method leads to proper time-history pressure curves and thus adequate structural design.Several approaches can be used to evaluate the properties of a blast wave from a given explosive charge. Taylor [4], and several authors after him (see, for example, [5–7]), managed to obtain adequate analytical blast wave models for strong explosions by making two main assumptions: (1) all explosive energy is concentrated in one central source point. (2) The blast overpressure should be much greater than the ambient pressure. They managed to establish a relation between the maximum pressure Pmax and the radial position of the shock R given by (1)Pmax=0.155meiR−3=0.155eiRm−1/3−3,where ei is the internal energy released per unit mass at one source point and m the total mass of the explosive. The factor R·m−1/3 is called the scaled distance and denoted Z hereafter. Unfortunately, neither of the assumptions of this model is fulfilled in the case of explosions implying conventional explosives. In fact, for the intensively studied TNT, both conditions are fulfilled in a fairly narrow scaled distance range around Z = 1 m·kg−1/3. At larger scaled distances, the overpressure is no more than ten times the ambient pressure and, close to the charge, the scaled distance is less than ten times the explosive diameter and so the point source assumption is no longer appropriate. Alternatively, empirical or semiempirical approaches have been used to obtain reliable loading values [8–10]. The technical manual from the US Army, UFC-3-340-2 [8], is mainly used as the state of the art in structural design offices and constitutes the most widely used experimental database for both civilian and military infrastructure design since it offers the most complete open database of blast wave parameters. Those empirical solutions have the benefit of relying on experimental data but have several drawbacks. The first one is the lack of reliable data for scale distance smaller than 0.5 m·kg−1/3. Indeed, it is more than questionable to use these databases for close-in structures [11, 12]. This is mainly due to the chemical reactions between the detonation products and the surrounding air. Secondly, data are available only for the incident blast wave and possible reflections cannot be easily taken into account. This prevents any of these solutions from being applied in confined or urban space. Finally, only simple geometries for both the charge and the structure can be considered (spherical charge and plan wall).Another widespread solution is the use of explicit numerical software. Nevertheless, in spite of all efforts made to enhance finite elements and finite volumes solvers, modeling different aspects of physics at different time scales remains a considerable challenge. In particular, a time step lower than the characteristic time of the detonation from a few orders of magnitude is required to obtain a proper model of the detonation process. Furthermore, an explicit numerical method makes calculation for each and every time step; thus, a small time step directly increases the overall computational time. Some hybrid solution using the UFC-3-340-2 data is added to numerical simulation in order to obtain a rapid loading for a given structure [13–15]. These methods are based on empirical data and suffer from the same problem as them, namely, they do not account for waves reflection and recombination in a confined or urban area and full simulation of the detonation process are needed.According to the ideal detonation process, the reaction front propagates within the high explosive at a speed of about 7000 m·s−1 for TNT [16, 17]. Just behind this reaction front, the Chapman-Jouguet [18] model predicts a pressure of PCJ = 21 GPa. The Chapman-Jouguet wave velocity is denoted as DCJ. In fact, the time step will be about less than a microsecond which leads to calculation times above the acceptable range for a commonly available HPC. Instead of losing accuracy and trying to poorly represent the detonation process, it is possible to only model the relaxation of the detonation products. In numerical software, the detonation products can be initially modeled as a homogeneous volume of high pressure and temperature gas, referred to as a balloon.For simplicity’s sake, the chemical composition of the detonation products is most of the time considered as air and they are not allowed to react with either themselves or the surrounding air. The multiphysics detonation problem is reduced to a fluid mechanics problem, namely, the relaxation of high pressure and temperature gases. An equation of state needs to be chosen to describe the behavior of the relaxation of the gas mixture within the balloon. Several choices are available. The most used equation of state in commercial software is the John Wilkins Lee (JWL) equation of state (2):(2)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV.The parameters of this equation are explosive dependent and have been extensively studied. The values given, for TNT at a density of 1630 kg·m−3, by Dobratz [16] and Souers [17] yield results in relatively good agreement with the UFC-3-340-2. Nevertheless, the JWL equation of state induces a large amount of computational effort by requiring high initial pressure state within the detonation products. These high pressure gradients lead towards small element sizes (less than 1 mm).In fact, for current HPC, the number of finite elements should not be higher than 5 million, which represents a volume of 10 cm3 in the case of a mesh element size of 1 mm, and is clearly not enough to model a full-scale infrastructure. In spite of the non-representation of the detonation, the balloon filled with JWL gas is not suitable for a fast and accurate numerical model. As soon as the 1950s, Brode [19] proposed to simplify the model and consider the ideal gas law to model the evolution of the detonation products. The ideal gas equation of state has been used over the past couple of decades to represent the relaxation of the detonation products in numerical simulations of realistic events [20, 21] or laboratory scale experiments [22]. All these authors have used ideal gas with different values of pressure, density, and energy for representing the physical state of the detonation products.In 2010s, Blanc [23] showed that the key parameters for monitoring the balloon-induced blast wave are density and specific energy of the explosive source. For Dobratz’s TNT parameters, this yields ρb = ρT·NT = 1630 kg·m−3 and e = eT·NT = 4.3 MJ·kg−1. These parameters account for the best accuracy for all relevant blast wave parameters (overpressure, impulse, and arrival time) with respect to the ideal detonation model. Other sets of parameters can be used to describe the relaxation of the high pressure gas inside the balloon thanks to the ideal gas equation of state (Chapman-Jouguet parameters: [24], isochoric compression [19], isothermic compression, etc.). The main approaches are presented in Table 1.Table 1
Compressed balloon possible sets of parameters.
Chapman-Jouguet parameters
Isochoric compression
Density-energy approach
Pressure (GPa)
21
21
2.8
Internal energy (MJ·kg−1)
32
4.3
4.3
Density (kg·m−3)
1630
12209
1630In this work, the approach of [23] is selected since it gives the better agreement with respect to the ideal detonation process, and relatively acceptable physical values for the density and pressure of the balloon gas in its initial state. As a consequence, the requirements for the element size are less constraining, but the resulting blast wave properties could be underestimated, especially in the near field. Another benefit of reducing the high pressure and energy gradient at the early stage of the relaxation is the mitigation of one unbalanced total energy, which is one of the remaining issues of the finite element model [25]. In fact, to deal with discontinuities in pressure, density, and energy, finite element software needs to add artificial viscosity [26] to smooth discontinuities in the pressure, energy, mass, and material speed in order to perform the derivation of the pressure function needed for finite element analysis. Also referred to as numerical damping, this will artificially reduce the total energy, especially at early stages of the blast wave propagation because of the large pressure gradients. Then it leads to higher errors made by commercial software.This work addresses this issue from an end-user standpoint, by increasing the radius of the compressed balloon up to ten times the initial radius of the explosive. This main novelty is leading towards the following:(i)
A better energy balance for finite element algorithm.(ii)
A better accuracy with experimental data from the UFC-3-340-2, [8] at later stage of the blast wave propagation.(iii)
A reduced computational time for engineering applications.First, the numerical model is presented and the ideal gas equation of state for the detonation products is compared with the classical JWL equation of state. In a second part, the influence of the balloon radius is investigated and all models are evaluated with respect to experimental data found in the UFC-3-340-2. All numerical results are compared to the chosen reference (the UFC-3-340-2), an experimental database accepted by authorities for structural design. A mesh sensitivity analysis is performed to validate the results and finally the impact on the energy balance is evaluated.


## 2. The Compressed Balloon Model

Altair Hyperworks§R suite is used to build all numerical models and the explicit finite element code. Radioss§R is used for the simulations. The solver used is call legacy and it is a classical Arbitrary Lagrangian Eulerian (ALE) finite element/finite volume mixed solver. The conservation equations for mass and energy are solved at the center of each element and the momentum conservation equation is solved at each node. The first model, the most used, takes JWL equation of state as presented in the introduction together with the chosen parameters for TNT from [16]; it is hereafter referred to as the “JWL model.” All other gases are modeled by a hydrodynamic fluid material using a polynomial equation of state (3).(3)P=Co+C1μ+C2μ2+C3μ3+C4+C5μe,where e is the energy per unit volume, µ is the dilation coefficient = ρ/ρ0 − 1, and P is the pressure.By choosingC0 = C1 = C2 = C3 = 0 and C4 = C5 = 0.4, the hydrodynamic polynomial law describes an ideal diatomic gas with a specific heat ratio γ = 1.4. In order to model the ambient state of the air, the density and the energy per unit volume are, respectively, set to ρair = 1.225 kg·m−3 and eair = 0.25 MJ·m−3 with respect to (3). As detailed in the introduction, for the compressed balloon approach, the density and the internal energy per unit volume of the balloon are set to the TNT values: eballoon = 7000 MJ·m−3 and ρballoon = 1630 kg·m−3. Parameter beta is defined by equation (4). It is the ratio of the balloon radius over the radius of the TNT equivalent spherical charge. This ratio will be varied between one and ten in this study:(4)β=RballonRTHN.The total mass and energy of the explosive source should not be modified by the increase of the balloon radius. Thus, the densityρβ and the specific energy Eβ are reduced by the coefficient β3:(5)ρβ=ETNTβ3,(6)Eβ=ETNTβ3.Two types of mesh are studied, with equivalent element size. The first one is a Cartesian mesh, which is of a more convenient use in structural design offices since it is easier to create and better fits with most fluid-structure interfaces. The second one is a radial structured mesh. Even if it is harder to obtain, the radial mesh is more accurate [27] and more adapted to immersed fluid-structure interfaces. Both radial and Cartesian meshes used for the balloon model are presented in Figure 1.Figure 1
Cartesian (a) and radial (b) mesh structure for the numerical simulation.The Cartesian mesh is straightforward; each element is a regular square and the mesh is not impacted by the increase inβ. On the opposite, the radial mesh is in two parts: the first one is Cartesian from the center of the charge to two-thirds of the explosive radius and the rest is a radial mesh. This radial mesh is adapted for each beta in order to reduce the mesh distortion at the transition between the Cartesian and radial parts. This has no consequence on the blast wave propagation because the mesh size in the propagation direction is kept identical in all β models.An axisymmetrical boundary condition is enforced on the vertical axe to ensure the spherical propagation of the blast wave. Likewise, a symmetry condition is applied on the horizontal axis in order to reduce the computational domain to one-fourth of the space.The scaled mesh size is defined as the actual mesh over the mass of equivalent high explosive at the cubicle root:(7)dxscaled=dxm1/3.The results of pressure as a function of time are highly dependent on the scale mesh size [25]. In practical use, 5 mm·kg−1/3 is among the smallest possible scaled mesh sizes for acceptable CPU time for civil engineering applications (around 5 million elements). It is then used in this section for the element size.In order to evaluate the simplified approach of theβ balloon, two parameters are considered, maximum overpressure and positive impulse. Indeed, both play an important part in the structural behavior of the loaded structure, one being usually more significant than the other depending on the natural frequencies of the structure. The overpressure is defined as the maximum of the pressure history and the positive impulse is the integration of the positive part of the pressure history curve defined more precisely by(8)I=∫t1t2PTdt,with t1 and t2 arbitrarily chosen as follows in order to avoid numerical noise: t1 the time where P>0.01Pmax and t2 the time where P<−0.01Pmax.Furthermore, two relative difference ratios are defined with respect to the UFC-3-340-2 in(9)δP=pmaxmodel−pmaxUFCpmaxUFC,(10)δI=Imaxmodel−ImaxUFCImaxUFC.Figures2 and 3 show, respectively, the evolution of δI and δp as functions of the scaled distance. Gauges have been placed in the numerical simulation to record the pressure history every 20 cm from the center of the charge and the relevant parameters have been extracted and compared with values given by the UFC-3-340-2. Furthermore, the finite volume method computed with the python module clawpack developed by Washington University [28] is also compared with the finite element calculation. All δI curves show a high gradient for scaled distances lower than 1 m·kg−1/3, which could be explained by the difficulty to get experimental data so close to the explosive charge [29].Figure 2
Relative difference with UFC-3-340-2 for impulseδI, for all models as a function of scaled distance Z.Figure 3
Relative differences with UFC-3-340-2 for overpressure, for all models as function of scaled distance.For scaled distances larger than 1.5 m·kg−1/3, δI curves tend to stabilize and exhibit a clear dependency on the mesh type (Cartesian or radial). As expected, the radial mesh gave better results and all δI for radial meshes are smaller than those for Cartesian meshes. The same trend is seen on the δp graphs but all differences between models are smaller.On both pressure and impulse, finite volumes give better results. Unfortunately, the blast-structure interaction between Lagrangian finite element representing the structural part and finite volume representing the air domain remains challenging for all commercial software. Structural design is hardly possible with finite volume involved for the air model, especially if nonlinear material and displacement are taken into account.


## 3. Results

In order to evaluate theβ-methods which could be used for civil engineering application, all data are extracted from the radial structured mesh at a scaled mesh size of 5 mm·kg−1/3. The Central Processing Unit (CPU) time of the computation of these models is reasonable for parametric studies (30 min for 10 ms of blast wave propagation on 16 CPU) and the radial mesh is commonly used for 3D industrial applications.First, Figures2 and 3 will be extended to all beta values from one to ten and the relative discrepancies between the beta models and UFC-3-340-2 (δ see equations (9) and (10)) will be studied, leading to a critical distance d(β) where the beta models become more accurate than the models using JWL equation of state for the description of the relaxation of detonation products, for a given mesh. In a second time, the mean values of δ for all scaled distances are plotted for both mesh types to get a more global idea of the influence of β in the scaled distance range of interest.


### 3.1. Influence of Balloon Radius

Figure4(a) shows that, close to the explosion center, δP quickly increases with increasing values of β. This is to be expected due to the lack of the proper distribution of the pressure and density. In the case of an explosive charge, once the blast wave is created all the mass and energy are concentrated within a few centimeters behind the shock discontinuity. On the opposite, with balloon models, the mass and the energy of the gaseous products of detonation are spread uniformly in the entire balloon volume and it takes a few-balloon radius to concentrate it on the shock front. At larger scaled distances, all δP are smaller than the δP of the model with JWL model. On Figure 4(b) iso-values of relative error to UFC-3-340-2 are extrapolated for intermediate values of mesh size and beta.Figure 4
Pressure relative differenceδP with UFC-3-340-2 for all beta models as function of distance (a) or distance and β (b).
(a)
(b)Fromβ = 3 to β = 10, all δp curves cross the δp(JWL) curve. This implies that there is a distance from which the beta model becomes more accurate for the predicted peak overpressure than the JWL one, for the considered mesh size. This scaled distance is hereafter denoted d and plotted in Figure 5.Figure 5
Scale distancedP as function of β.For scaled distances greater than 1.5 m·kg−1/3, the beta balloon models presented in this study give an overpressure closer to the one given in the UFC-3-340-2. This means that if a concrete wall is far enough from the explosive center, the beta balloon model will give better results than the model using JWL equation of state for the relaxation of the detonation products. There are no points before β = 3 in Figure 5 because for βi 3 beta balloon models are better than JWL model even close to the explosive center.Concerning the impulse on Figure6, high positive δI is observed in the vicinity of the high explosive, which means that the impulse is severely overestimated before about 0.8 m·kg−1/3. Furthermore, the impulse results are more complicated to analyze due to the numerical damping changes the blast wave shape. See Figure 7, where the pressure-time curves and impulse-time (computed as defined in 2) of β = 5 and JWL model are plotted for two scaled distances of 1.4 m·kg−1/3 and 2.6 m·kg−1/3. Indeed, the rising time (defined as the time from ambient pressure to maximum pressure) is increased between both measurements at 1.4 and 2.6 m·kg−1/3. The model considering the JWL equation of state exhibits a larger increase than the β = 5 model. Thus the total area under the pressure-time curve (the impulse) is greater and closer to values given by the reference UFC-3-340-2 for the JWL equation of state. This accuracy is artificial because the shape of the profile is changed and so is the loading on the structure. In order to get a meaningful value of distance dI for the structural behavior, the correlation between the increased rising time and the structural design should be studied. Considerations such as materials or natural frequencies are not taken into account in the present work and thus changes in the blast profile are not addressed by the authors in this study.Figure 6
Impulse relative differenceδI with UFC-3-340-2 for all β models.Figure 7
Temporal evolution of overpressure and impulse for both JWL andβ = 5 models at two scaled distances: 1.4 and 2.6 m·kg−1/3.In some practical applications, especially for confined explosions, there is no particular distance of interest and the blast wave should be accurate in the entire air domain. For example, in a room full of equipment, the peak overpressure and the impulse should be as accurate as possible at every scaled distance of each standoff point. This is why the mean of the relative differencesδi and δp is studied in the next section.


### 3.2. Average Discrepancies

For a given mesh size of 5 mm·kg−1/3, Figures 8 and 9 show the average performance of all balloon models for, respectively, impulse and overpressure. The performance is defined as the average value of all relative differences δP or δI for all scaled distances from 0.5 to 5m·kg1/3. The value β = 0 corresponds to the JWL equation of state model.Figure 8
Impulse average error with UFC-3-340-2 as a function of balloon sizeβ.Figure 9
Pressure average error with UFC-3-340-2 as a function of balloon sizeβ.The behavior for both meshes is really different regarding the blast wave impulse; see Figure8. Indeed, the radial mesh shows increase relative differences with respect to UFC-3-340-2 for increasing beta while the other, the Cartesian, shows decreasing relative differences with UFC-3-340-2.On the opposite, both curves look similar for the overpressure; see Figure9. Both show an extremum for values of β depending on the kind of mesh considered: β = 2 for the Cartesian mesh and β = 5 for the radial one. This results from the trade-off between poor accuracy at early stage of the simulation and accuracy of the big balloon at larger scaled distances. In conclusion for this section, there is an optimal β value for a given mesh (element size and mesh type). In the discussion, these results will be extended to other element sizes.


### 3.3. Computational Time

The main parameter ruling the calculation time is the time step. In RadiossR, as in most finite element software, a needed condition must be fulfilled in order to ensure numerical stability. This condition was established by Courant, Friedrich, and Levy (CFL) [30]. Unfortunately, the CFL condition is required to ensure the stability of the numerical scheme but it is not sufficient and it is possible that the time step should be decreased even further in order to allow the calculation to run until the end. According to Radioss R documentation, the stability condition is calculated for each element following the CFL following equation (11). Then the smaller time step ∆t is used for the whole calculation.(11)Δt=kΔlc+v−w,with ∆t being the time step, ∆l the element size, v the material velocity, w the ALE grid velocity, and k the CFL coefficient 1.c is the speed of sound, calculated in (12).(12)c=∂P∂p,with ρ being the density and P the pressure.At early stages of the blast wave propagation, results obtained with smallβ models, as well as with the model considering the JWL equation of state (β = 0), exhibit larger overpressures than large β models. This leads to smaller time steps. Figure 10 shows the specific example of two radial balloon models for β = 1 and β = 10. Before few hundredths of milliseconds, the time step is smaller for β = 1. However, at later stages of the blast wave propagation, the trend is reversed and the β = 10 model time step decreases due to higher maximum overpressures than for small β balloon at the same time. The doted curves show the mean value of the time step that account for both parts of the calculation.Figure 10
Time step temporal evolution for bothβ = 1 and β = 10 balloon models. Radial model.The later stage (after 4 ms) supersedes the early one and the mean time step (dotted line) is smaller for beta = 10; thus the calculation lasts up to twice as long as the calculation with beta = 1. The k-parameter presented in equation (11) should be reduced to 0.1 to ensure numerical stability due to mesh distortion of the radial mesh. Nevertheless, large beta balloon models decrease pressure gradients and thus allow the use of greater k than for the model with JWL equation of state or small beta balloon models. In Figure 11, the increase in k is taken into account for the calculation of the mean time step in the radial mesh. On the contrary, the Cartesian mesh exhibits no numerical stability issues preventing the calculation from running until the end.Figure 11
Evolution of the mean time step as a function of the balloon size for Cartesian and radial meshes.The evolution of the mean time step and thus the CPU time is therefore different for both meshes. Indeed, broadly speaking, the increase in the CFL condition’s parameterk presented previously benefits only the radial mesh. Then, the larger the balloon, the higher the time step and thus the lower the calculation time. On the opposite, this parameter is kept as high as possible (k = 0.9) for the Cartesian mesh and thus the trend is reversed and the larger the balloon, the lower the time step. To conclude this part, the amplitude of variation of δP and δI is really large in the first few dozens of cm·kg−1/3. Nevertheless, β balloon models give better overpressure with respect to the UFC-3-340-2 after a scaled distance d between 0.8 and 1.5 m·kg−1/3. The average values of δP and δI give results within 10% for all beta values and cannot be used to evaluate the β model. The first goal of using these methods was to simplify, and thus fasten, the calculation. This is obtained by increasing the numerical stability for the radial mesh. Optimal values of β can be found based on the maximum overpressure. In this respect, δP is the smallest between beta = 5 and beta = 6 for the Cartesian mesh and around beta = 2 for the radial mesh. This result is important from a practical point of view, when trying to fasten or increase the accuracy of calculation. Unfortunately, the scaled mesh size is the key parameter for calculation which deals with blast wave effects on large structures and the beta balloon solution has to be validated for a wider range of reduced mesh size.


## 3.1. Influence of Balloon Radius

Figure4(a) shows that, close to the explosion center, δP quickly increases with increasing values of β. This is to be expected due to the lack of the proper distribution of the pressure and density. In the case of an explosive charge, once the blast wave is created all the mass and energy are concentrated within a few centimeters behind the shock discontinuity. On the opposite, with balloon models, the mass and the energy of the gaseous products of detonation are spread uniformly in the entire balloon volume and it takes a few-balloon radius to concentrate it on the shock front. At larger scaled distances, all δP are smaller than the δP of the model with JWL model. On Figure 4(b) iso-values of relative error to UFC-3-340-2 are extrapolated for intermediate values of mesh size and beta.Figure 4
Pressure relative differenceδP with UFC-3-340-2 for all beta models as function of distance (a) or distance and β (b).
(a)
(b)Fromβ = 3 to β = 10, all δp curves cross the δp(JWL) curve. This implies that there is a distance from which the beta model becomes more accurate for the predicted peak overpressure than the JWL one, for the considered mesh size. This scaled distance is hereafter denoted d and plotted in Figure 5.Figure 5
Scale distancedP as function of β.For scaled distances greater than 1.5 m·kg−1/3, the beta balloon models presented in this study give an overpressure closer to the one given in the UFC-3-340-2. This means that if a concrete wall is far enough from the explosive center, the beta balloon model will give better results than the model using JWL equation of state for the relaxation of the detonation products. There are no points before β = 3 in Figure 5 because for βi 3 beta balloon models are better than JWL model even close to the explosive center.Concerning the impulse on Figure6, high positive δI is observed in the vicinity of the high explosive, which means that the impulse is severely overestimated before about 0.8 m·kg−1/3. Furthermore, the impulse results are more complicated to analyze due to the numerical damping changes the blast wave shape. See Figure 7, where the pressure-time curves and impulse-time (computed as defined in 2) of β = 5 and JWL model are plotted for two scaled distances of 1.4 m·kg−1/3 and 2.6 m·kg−1/3. Indeed, the rising time (defined as the time from ambient pressure to maximum pressure) is increased between both measurements at 1.4 and 2.6 m·kg−1/3. The model considering the JWL equation of state exhibits a larger increase than the β = 5 model. Thus the total area under the pressure-time curve (the impulse) is greater and closer to values given by the reference UFC-3-340-2 for the JWL equation of state. This accuracy is artificial because the shape of the profile is changed and so is the loading on the structure. In order to get a meaningful value of distance dI for the structural behavior, the correlation between the increased rising time and the structural design should be studied. Considerations such as materials or natural frequencies are not taken into account in the present work and thus changes in the blast profile are not addressed by the authors in this study.Figure 6
Impulse relative differenceδI with UFC-3-340-2 for all β models.Figure 7
Temporal evolution of overpressure and impulse for both JWL andβ = 5 models at two scaled distances: 1.4 and 2.6 m·kg−1/3.In some practical applications, especially for confined explosions, there is no particular distance of interest and the blast wave should be accurate in the entire air domain. For example, in a room full of equipment, the peak overpressure and the impulse should be as accurate as possible at every scaled distance of each standoff point. This is why the mean of the relative differencesδi and δp is studied in the next section.


## 3.2. Average Discrepancies

For a given mesh size of 5 mm·kg−1/3, Figures 8 and 9 show the average performance of all balloon models for, respectively, impulse and overpressure. The performance is defined as the average value of all relative differences δP or δI for all scaled distances from 0.5 to 5m·kg1/3. The value β = 0 corresponds to the JWL equation of state model.Figure 8
Impulse average error with UFC-3-340-2 as a function of balloon sizeβ.Figure 9
Pressure average error with UFC-3-340-2 as a function of balloon sizeβ.The behavior for both meshes is really different regarding the blast wave impulse; see Figure8. Indeed, the radial mesh shows increase relative differences with respect to UFC-3-340-2 for increasing beta while the other, the Cartesian, shows decreasing relative differences with UFC-3-340-2.On the opposite, both curves look similar for the overpressure; see Figure9. Both show an extremum for values of β depending on the kind of mesh considered: β = 2 for the Cartesian mesh and β = 5 for the radial one. This results from the trade-off between poor accuracy at early stage of the simulation and accuracy of the big balloon at larger scaled distances. In conclusion for this section, there is an optimal β value for a given mesh (element size and mesh type). In the discussion, these results will be extended to other element sizes.


## 3.3. Computational Time

The main parameter ruling the calculation time is the time step. In RadiossR, as in most finite element software, a needed condition must be fulfilled in order to ensure numerical stability. This condition was established by Courant, Friedrich, and Levy (CFL) [30]. Unfortunately, the CFL condition is required to ensure the stability of the numerical scheme but it is not sufficient and it is possible that the time step should be decreased even further in order to allow the calculation to run until the end. According to Radioss R documentation, the stability condition is calculated for each element following the CFL following equation (11). Then the smaller time step ∆t is used for the whole calculation.(11)Δt=kΔlc+v−w,with ∆t being the time step, ∆l the element size, v the material velocity, w the ALE grid velocity, and k the CFL coefficient 1.c is the speed of sound, calculated in (12).(12)c=∂P∂p,with ρ being the density and P the pressure.At early stages of the blast wave propagation, results obtained with smallβ models, as well as with the model considering the JWL equation of state (β = 0), exhibit larger overpressures than large β models. This leads to smaller time steps. Figure 10 shows the specific example of two radial balloon models for β = 1 and β = 10. Before few hundredths of milliseconds, the time step is smaller for β = 1. However, at later stages of the blast wave propagation, the trend is reversed and the β = 10 model time step decreases due to higher maximum overpressures than for small β balloon at the same time. The doted curves show the mean value of the time step that account for both parts of the calculation.Figure 10
Time step temporal evolution for bothβ = 1 and β = 10 balloon models. Radial model.The later stage (after 4 ms) supersedes the early one and the mean time step (dotted line) is smaller for beta = 10; thus the calculation lasts up to twice as long as the calculation with beta = 1. The k-parameter presented in equation (11) should be reduced to 0.1 to ensure numerical stability due to mesh distortion of the radial mesh. Nevertheless, large beta balloon models decrease pressure gradients and thus allow the use of greater k than for the model with JWL equation of state or small beta balloon models. In Figure 11, the increase in k is taken into account for the calculation of the mean time step in the radial mesh. On the contrary, the Cartesian mesh exhibits no numerical stability issues preventing the calculation from running until the end.Figure 11
Evolution of the mean time step as a function of the balloon size for Cartesian and radial meshes.The evolution of the mean time step and thus the CPU time is therefore different for both meshes. Indeed, broadly speaking, the increase in the CFL condition’s parameterk presented previously benefits only the radial mesh. Then, the larger the balloon, the higher the time step and thus the lower the calculation time. On the opposite, this parameter is kept as high as possible (k = 0.9) for the Cartesian mesh and thus the trend is reversed and the larger the balloon, the lower the time step. To conclude this part, the amplitude of variation of δP and δI is really large in the first few dozens of cm·kg−1/3. Nevertheless, β balloon models give better overpressure with respect to the UFC-3-340-2 after a scaled distance d between 0.8 and 1.5 m·kg−1/3. The average values of δP and δI give results within 10% for all beta values and cannot be used to evaluate the β model. The first goal of using these methods was to simplify, and thus fasten, the calculation. This is obtained by increasing the numerical stability for the radial mesh. Optimal values of β can be found based on the maximum overpressure. In this respect, δP is the smallest between beta = 5 and beta = 6 for the Cartesian mesh and around beta = 2 for the radial mesh. This result is important from a practical point of view, when trying to fasten or increase the accuracy of calculation. Unfortunately, the scaled mesh size is the key parameter for calculation which deals with blast wave effects on large structures and the beta balloon solution has to be validated for a wider range of reduced mesh size.


## 4. Discussion

### 4.1. Mesh Sensitivity Analysis

All blast wave results are highly dependent on the element size. All presented methods try to enable a larger mesh size for a given accuracy and/or a better accuracy for a given mesh size. A mesh sensitivity analysis is therefore required to evaluate the performances of the proposed beta balloon model with reduced element size from 1 mm·kg−1/3 up to 15 mm·k−1/3.Firstly, the mean value ofδP is plotted, see Figure 12, for all the following scaled element sizes: 1 mm·kg−1/3; 2 mm·kg−1/3; 5 mm·kg−1/3; 10 mm·kg−1/3; and 15 mm·kg−1/3. In Figure 12(b), iso-values of relative error to UFC-3-340-2 are extrapolated for intermediate values of mesh size and beta.Figure 12
Relative peak of overpressure difference with UFC-3-340-2 for all beta parameter and all studied scaled mesh sizes.
(a)
(b)On one side, for the smallest scaled mesh size, the numerical damping is limited and smallβ models are more accurate by about 20%. On the other, for larger element sizes, the numerical damping is more important and the gain from large β is noticeable with an ideal size of β = 5. As an example, starting at point (β = 1, dx = 7 mm−1/3), the iso-value of the average relative difference with the UFC-3-340-2 shows an extremum in scaled mesh size of 14 mm−1/3. As a rough estimate, twice the scaled mesh size means a twice longer time step and so a reduced computational time. On this model, the total elapsed time of the computation is reduced to 80% as fast without losing any average accuracy. The differences between 1 m·kg−1/3 and 15 m·kg−1/3 decrease from almost 30% for β = 1 to less than 10% for β = 10 which indicates that models with large values of β are less sensitive to element size. In fact, on one side, at least 5 elements are required within the blast wavelength in order to represent properly the blast wave. On the other side, balloons with small values of β create the blast wave with smaller wavelength than balloons with large values of β due to increasing wave length with scaled distances [8, 9]. Thus, small β models require smaller elements to represent the blast profile than large β models and so they are more sensitive to element size.Figure13 shows the pressure as a function of time for all mesh sizes at a given scaled distance of 3 m·kg−1/3. The artificial viscosity is more visible for larger mesh sizes. As a consequence, in spite of the larger mesh size, the impulse will increase and the relative difference with respect to UFC-3-340-2 will decrease. Yet, this is only due to numerical artifacts that should not lead to hasty conclusions.Figure 13
Overpressure at 3 m·kg−1/3 for beta = 3, obtained with various reduced mesh sizes.Figure14 is an extension of Figure 5 and shows the distance d where the balloon model becomes more accurate than the JWL model for the same scaled mesh size. All iso-values are computed using the linear extrapolation method. “Never better” means the distance is higher than 3.7 m·kg−1/3, and thus the JWL model is more accurate than every β model, or δp(JWL) < δp(β = [1, …, 10]). On the other side, “always better” means δp(β = [1, …, 10]) < δp(JWL).Figure 14
Scaled distance d as function of the beta parameter and the reduced mesh size dx.
(a)
(b)The blue domain in Figure14 shows the parameters leading to a more accurate overpressure estimation obtained by the balloon model with respect to the JWL one. For dx = 1 mm·kg−1/3, only the beta = 1 balloon is more accurate than the JWL model so there is only one point on the curve in Figure 11. For all other mesh sizes, the distance d increases with the mesh size. Nevertheless, for mesh sizes larger than 5 mm·kg−1/3, distance d is less than 1.5 m·kg−1/3, which is a likely standoff scaled distance for the first target wall.The beta balloon model performs better when the mesh is not good enough and so it takes over the large element size or the mesh distortion to ensure better accuracy and shorter calculation times.


### 4.2. Energy Balance

The ALE algorithm used in all presented simulations is of interest because it allows a direct coupling between the fluid and the structure. The momentum equation is solved at mesh nodes so the material speed is also calculated at nodes, as it is with the Lagrangian algorithm. It is then easy to interface the fluid and the structure by just equating the material speed of the two parts in the direction normal to the interface. Unfortunately, the counterpart of this is the non-conversation of the total energy in the fluid part that is related to the discrete integration of the conservation equations made by finite element method [31]. Figure 15 shows the total energy as a function of time for all β models together with the JWL model.Figure 15
Total energy variation as a function of the time. Cartesian mesh.The initial energy is computed on the Cartesian mesh because the mesh is constant from aβ value to another. Nevertheless, no relationship giving the specific (nor per unit volume) energy set in the absorbing condition (noted ABS) has been found, so the value is extracted from the β balloon models. In fact a single element model has been made with only one absorbing cell of 1 m × 1 m; thus the total energy is easily extracted directly per unit volume. eABS = −0.25 MJ·m−3. In equations (13) and (14), all capital E stand for total energy and all small e stand for per unit volume.(13)Etot=∑nEn=Etnt+eairVair+Eabs.with Etnt = 4.3 MJ the energy of 1 kg of TNT, eair = 0.25 MJ·kg−3 the energy of air per unit volume, and EABS = 0.589 MJ the energy of the absorbing elements.Furthermore, energy values for the axisymmetric model are given per radian and only one half of the spatial domain is modeled. So, forβ = 1,(14)Eaxi=Etot2π=15.91MJ⋅rad−1.The difference in the initial energy balance should be really small. Indeed, one of the strong hypotheses of theβ balloon is that the total energy of the high explosive is constant and spread on a larger volume. Furthermore, the total volume of numerical domain is also constant in all simulations. As a consequence, the total amount of air is reduced for larger beta model and thus also the initial total energy. This small variation in the amount of air within the simulation explains the small differences at t = 0 ms in Figure 15.The energy loss due to the ALE algorithm in the first millisecond is more significant for small balloons and for the model considering the JWL equation of state. The loss is directly correlated to the actual energy, pressure, and density inside a specific element. As a consequence, the largeβ balloon spreads more energy and pressure gradients, the loss per element is less important, and thus those models present a better energy balance. The overpressure decay is also less important from 2 ms, which corresponds to a maximum overpressure of about 0.14 MPa. This explains the numerical damping reduction and the stabilization of the energy balance. The presented β balloon models manage to mitigate the numerical damping by spreading the initially high pressure and energy over a larger volume. By reducing the energy per unit volume within a single element, the gradient between two elements is also reduced and thus the energy loss proportional to this gradient is also reduced. The increase in energy balance leads to smaller δ as presented in Section 3, especially at larger scaled distances.


## 4.1. Mesh Sensitivity Analysis

All blast wave results are highly dependent on the element size. All presented methods try to enable a larger mesh size for a given accuracy and/or a better accuracy for a given mesh size. A mesh sensitivity analysis is therefore required to evaluate the performances of the proposed beta balloon model with reduced element size from 1 mm·kg−1/3 up to 15 mm·k−1/3.Firstly, the mean value ofδP is plotted, see Figure 12, for all the following scaled element sizes: 1 mm·kg−1/3; 2 mm·kg−1/3; 5 mm·kg−1/3; 10 mm·kg−1/3; and 15 mm·kg−1/3. In Figure 12(b), iso-values of relative error to UFC-3-340-2 are extrapolated for intermediate values of mesh size and beta.Figure 12
Relative peak of overpressure difference with UFC-3-340-2 for all beta parameter and all studied scaled mesh sizes.
(a)
(b)On one side, for the smallest scaled mesh size, the numerical damping is limited and smallβ models are more accurate by about 20%. On the other, for larger element sizes, the numerical damping is more important and the gain from large β is noticeable with an ideal size of β = 5. As an example, starting at point (β = 1, dx = 7 mm−1/3), the iso-value of the average relative difference with the UFC-3-340-2 shows an extremum in scaled mesh size of 14 mm−1/3. As a rough estimate, twice the scaled mesh size means a twice longer time step and so a reduced computational time. On this model, the total elapsed time of the computation is reduced to 80% as fast without losing any average accuracy. The differences between 1 m·kg−1/3 and 15 m·kg−1/3 decrease from almost 30% for β = 1 to less than 10% for β = 10 which indicates that models with large values of β are less sensitive to element size. In fact, on one side, at least 5 elements are required within the blast wavelength in order to represent properly the blast wave. On the other side, balloons with small values of β create the blast wave with smaller wavelength than balloons with large values of β due to increasing wave length with scaled distances [8, 9]. Thus, small β models require smaller elements to represent the blast profile than large β models and so they are more sensitive to element size.Figure13 shows the pressure as a function of time for all mesh sizes at a given scaled distance of 3 m·kg−1/3. The artificial viscosity is more visible for larger mesh sizes. As a consequence, in spite of the larger mesh size, the impulse will increase and the relative difference with respect to UFC-3-340-2 will decrease. Yet, this is only due to numerical artifacts that should not lead to hasty conclusions.Figure 13
Overpressure at 3 m·kg−1/3 for beta = 3, obtained with various reduced mesh sizes.Figure14 is an extension of Figure 5 and shows the distance d where the balloon model becomes more accurate than the JWL model for the same scaled mesh size. All iso-values are computed using the linear extrapolation method. “Never better” means the distance is higher than 3.7 m·kg−1/3, and thus the JWL model is more accurate than every β model, or δp(JWL) < δp(β = [1, …, 10]). On the other side, “always better” means δp(β = [1, …, 10]) < δp(JWL).Figure 14
Scaled distance d as function of the beta parameter and the reduced mesh size dx.
(a)
(b)The blue domain in Figure14 shows the parameters leading to a more accurate overpressure estimation obtained by the balloon model with respect to the JWL one. For dx = 1 mm·kg−1/3, only the beta = 1 balloon is more accurate than the JWL model so there is only one point on the curve in Figure 11. For all other mesh sizes, the distance d increases with the mesh size. Nevertheless, for mesh sizes larger than 5 mm·kg−1/3, distance d is less than 1.5 m·kg−1/3, which is a likely standoff scaled distance for the first target wall.The beta balloon model performs better when the mesh is not good enough and so it takes over the large element size or the mesh distortion to ensure better accuracy and shorter calculation times.


## 4.2. Energy Balance

The ALE algorithm used in all presented simulations is of interest because it allows a direct coupling between the fluid and the structure. The momentum equation is solved at mesh nodes so the material speed is also calculated at nodes, as it is with the Lagrangian algorithm. It is then easy to interface the fluid and the structure by just equating the material speed of the two parts in the direction normal to the interface. Unfortunately, the counterpart of this is the non-conversation of the total energy in the fluid part that is related to the discrete integration of the conservation equations made by finite element method [31]. Figure 15 shows the total energy as a function of time for all β models together with the JWL model.Figure 15
Total energy variation as a function of the time. Cartesian mesh.The initial energy is computed on the Cartesian mesh because the mesh is constant from aβ value to another. Nevertheless, no relationship giving the specific (nor per unit volume) energy set in the absorbing condition (noted ABS) has been found, so the value is extracted from the β balloon models. In fact a single element model has been made with only one absorbing cell of 1 m × 1 m; thus the total energy is easily extracted directly per unit volume. eABS = −0.25 MJ·m−3. In equations (13) and (14), all capital E stand for total energy and all small e stand for per unit volume.(13)Etot=∑nEn=Etnt+eairVair+Eabs.with Etnt = 4.3 MJ the energy of 1 kg of TNT, eair = 0.25 MJ·kg−3 the energy of air per unit volume, and EABS = 0.589 MJ the energy of the absorbing elements.Furthermore, energy values for the axisymmetric model are given per radian and only one half of the spatial domain is modeled. So, forβ = 1,(14)Eaxi=Etot2π=15.91MJ⋅rad−1.The difference in the initial energy balance should be really small. Indeed, one of the strong hypotheses of theβ balloon is that the total energy of the high explosive is constant and spread on a larger volume. Furthermore, the total volume of numerical domain is also constant in all simulations. As a consequence, the total amount of air is reduced for larger beta model and thus also the initial total energy. This small variation in the amount of air within the simulation explains the small differences at t = 0 ms in Figure 15.The energy loss due to the ALE algorithm in the first millisecond is more significant for small balloons and for the model considering the JWL equation of state. The loss is directly correlated to the actual energy, pressure, and density inside a specific element. As a consequence, the largeβ balloon spreads more energy and pressure gradients, the loss per element is less important, and thus those models present a better energy balance. The overpressure decay is also less important from 2 ms, which corresponds to a maximum overpressure of about 0.14 MPa. This explains the numerical damping reduction and the stabilization of the energy balance. The presented β balloon models manage to mitigate the numerical damping by spreading the initially high pressure and energy over a larger volume. By reducing the energy per unit volume within a single element, the gradient between two elements is also reduced and thus the energy loss proportional to this gradient is also reduced. The increase in energy balance leads to smaller δ as presented in Section 3, especially at larger scaled distances.


## 5. Conclusion

Based on the initial work of Brode and the latest papers, a compressed balloon has been studied to replace the detonation process of high explosive in numerical simulations of blast waves impacting structures. The previous models did not account for the numerical damping acting in every finite element explicit scheme. The main goal of this paper was to address this issue. To do so, the radius of the ideal gas balloon has been increased up to ten times the initial explosive radius. Two structured meshes were considered: Cartesian and radial. Increasing the size of the compressed balloon reduces all pressure and energy gradients at early stages of the blast wave propagation. The gradient reduction leads to a reduction of artificial viscosity and numerical damping and thus more accurate blast wave at large reduced distances. Indeed, there is a specific scaled distance where the lack of proper equation of state and volume of the detonation product is outweighed by the better conservation of the total energy. After about 1m.kg−1/3, large compressed balloons start to be in better accordance with respect to the UFC-3-340-2. All blast wave numerical simulations are highly dependent on the mesh size, respecting the trend the smaller the size, the more accurate the result. The presented work enables a larger mesh size for the same accuracy and/or a better accuracy for the same mesh size if the β parameter is increased. To be more specific, it is possible to double the mesh size without losing any accuracy on the maximum overpressure. In practical use, an optimum β has been found to be around β = 5 which could allow a reduction of computer time around 80% for equivalent accuracy with respect to the UFC-3-340-2. Furthermore, for an equivalent mesh size, the large beta models reduce the calculation time especially for radial meshes that face a challenge for achieving numerical stability. For radius greater than or equal to six times the initial (Brode) radius, the calculation time is six or seven times faster than the initial calculation with the JWL equation of state used for the relaxation of the detonation gases involved in the modeling of the high explosive charge. In addition, to provide faster computational times, the balloon model does not require any specific knowledge, set of parameters, on the explosive charge except its TNT equivalent.This compressed balloon approach could be applied on finite volume algorithm which is, theoretically, conservative for the total energy and that does not need artificial viscosity. Some of those algorithms are under investigation but the finite volume-finite element coupling between, respectively, Eulerian and Lagrangian parts is challenging, specifically for industrial and commercial software. Another approach could be to increase the initial energy in order to achieve the wanted blast energy at the time the shock front is impacting the structure but it would lead to inaccurate results at other distances.



---

*Source: 1497632-2020-05-12.xml*