Minimum mesh design cr iter ia for blast wave development and structural response – MMALE method 1

The Multi Material Arbitrary Lagrange Euler (MMALE) method is widely used method for numerical investigation of structural response under blast loading. However the method is very demanding for use at the other hand. In this paper are presented the results of the detailed numerical investigation in order to simplify some decisions contributing to the accuracy and efficiency of this model. The influence of mesh properties (particularly mesh size, its biasing and distance of the boundary (DoB) conditions from the deforming structure) on blast wave loading parameters and structural response is investigated in detail and based on the results minimum mesh design criteria is proposed. The results obtained are presented as a function of the scaled distance and additionally related to the radius of the charge. Validation studies were also done successfully.


INTRODUCTION
Examination of blast loading by explosive devices is of great importance in both, design and protection of different military equipment.Beside this, probability of terroristic attack and exposure to blast loading of public structures like bridges, metros, embassies, and other governmental buildings, is dramatically increased especially in the last two decades.Therefore, it is necessary to improve the response of such structures to this high-intensity, short-duration loads.However, preparation and performance of blast loading experiments using high energetic materials is extremely expensive and time consuming even for small scale tests in laboratory conditions, making the numerical analysis the most valuable tool for examination.On the other hand, precise modeling of blast loading of structures is one of the most complex dynamic analyses.However, considerable effort has been made in the last two decades in the field of both, numerical approximation Latin American Journal of Solids and Structures 11 (2014) 1999-2017 of blast loads and numerical response of structures to this kind of loads.Being two separated branches in the past they are now interactively coupled in the most advanced finite element solvers like Abaqus, AUTODYN, LS-DYNA.
The simplest numerical approximations of blast loads are triangular, rectangular, or bilinear pressure-time curves, more often used in the past but also at present (Biggs (1964), Jama et al. (2009), Krauthammer and Ku (1996), Louca et al. (1996), Tavakoli and Kiakojouri (2014)), and in analytical methods (Biggs (1964), Jones (1997), Smith and Hetherington (1994)).More precise, build-in numerical representation of blast loads appeared when Randers-Pehrson and Bannister (1997) implemented previous work of Kingery and Bulmash (1984) (ConWep blast model Bruce and Jon (1991)) in LS-DYNA (LSTC (2013a, b)).Shortly after that authors started using it for simulation in different scenarios, comparing their experimental results with the numerical.Xu and Lu (2006) applied ConWep to study the spallation of reinforced concrete, Neuberger et al. (2007a), to examine scaling of air blast loaded plates, Børvik et al. (2009) investigated the blast response of ISO container, Hussein et al. (2011) investigated the response of Retrofitted RC columns.Nearly at the same time, a new Arbitrary Lagrange Euler (ALE) method for modeling explosions became widely available, attacking the reality of the physical process itself.This method allows explosive and air to be represented as separated materials with a possibility to shear the same elements in the model (Alia and Souli (2006), Souli et al. (2000)).All the processes, starting with the detonation, blast wave formation and its interaction with the surrounding structures through Fluid Structure Interaction (FSI) algorithm were joined together in this model.This contributed to a more precise approximation of blast loads allowing more complicated real scenarios to be modeled and more optimized design to be achieved.Using ALE method, Zhao et al. (2012) investigated the response of reinforced concrete containment of an nuclear power plant under internal blast loading, Neuberger et al. (2007a, b) applied it to study armor plate response under air blast and buried charges detonation.Zakrisson et al. (2012), Zakrisson et al. (2011) investigated the plate response under blast loading modeling cylindrical explosive geometry in a steel pot according to NATOs standard AEP-55) requirement.Soutis et al. (2011) studied the response of GLARE panels to air blast loading using both, ALE and ConWep methods.In other studies ALE method was also used in variety of different applications and purposes: Chafi et al. (2009), Chung Kim Yuen et al. (2012), Fox et al. (2011), Jayasinghe et al. (2013), Langdon et al. (2010), Liu et al. (2012), Ma et al. (2013), Pi et al. (2012), Spranghers et al. (2013), Tai et al. (2011).In order to efficiently represent FSI on higher scaled distances, other techniques like coupling of the empirical ConWep model with ALE (Todd (2010)) or ALE mapping techniques are used (Luccioni et al. (2006), Zakrisson et al. (2011)).Some simplifications based on compressibility effect in FSI were also proposed by Kambouchev et al. (2006) and applied by Lin et al. (2013).Few authors (Barsotti et al. (2012), Genevieve and Robert (2008), Wang et al. (2005), Xu and Liu (2008)) used Smooth Particle Hydrodynamics (SPH) method to model blast response.The corpuscular approach was also proposed (Olovsson et al. (2010)) and successfully applied (Børvik et al. (2011)) in blast response scenario.However, MMALE method remains one of the most used methods for blast response analysis.Despite the fact that the accuracy of the numerical results is strongly influenced by the geometrical mesh properties, quite large differences can be observed in the above mentioned studies.These can lead to unreliable results or unnecessary large and nonefficient models on the other hand.As reported in Zukas and Scheffler (2000), given the same task of air blast loaded silo door, four users provided quite different results (reaching even 80% difference among each other) for maximum displacement.Therefore, special care should be given on the most influential numerical parameters in order achieve betters uniqueness and deliver proper loading to the structure while keeping the model cost effective.
In this paper, explosion of 6 kg TNT in free air was simulated using empirical ConWep model, axisymmetric MMALE, and 3D MMALE models using different mesh sizes.The results for incident blast wave parameters corresponding to close range scaled distances of Z=(0.11-0.385)m/kg 1/3 were compared with the experimental results available in US Army Defense Design Manual UFC-3-340-02 (UFC ( 2008)) (previously known as TM-5 1300).Numerical asymptotic values of the incident blast wave parameters are presented as a function of number of elements per radius of the charge.After the validation of free air model, the same was upgraded adding a steel plate.This enabled FSI to appear, allowing reflected blast wave parameters to be tracked.The examination was continued further, investigating the influence of DoB conditions from the loaded structure, and the influence of mesh size and its biasing on the reflected blast loading parameters, and consequently on structural displacement.Based on the results, recommendations for more efficient modeling of blast wave phenomena using MMALE method are also given.

BLAST WAVE PARAMETERS
When explosion occurs in free air, the gaseous products are rapidly expanding out of the initially occupied volume, creating a blast wave.This incident wave has initial velocity close to the detonation velocity of the explosive (6-10 km/s) (Zukas and Walters (2002)).As it is shown on Figure 1, the incident pressure (also called side on pressure) instantly arrives at a certain distance in space at arrival time t A , with its peak value P SO , after which exponentially decreases.The area under pressure-time curve represents the specific impulse from which only the positive one is considered responsible for structural deformation.The negative specific impulse is often neglected due to its small size.and Structures 11 (2014) 1999-2017 The simplest form of blast wave in free air is described with Friedlander wave equation (1); in which P 0 is the atmospheric pressure, P SO is the peak positive pressure, t is current time and t d is time duration of positive pressure.
If the blast wave strikes an object on its way, it reflects and it delivers reflected pressure P R , which could be twice to eight times stronger than the incident one.This appears because the particles at the front of the blast wave are stopped by the structure but they are still forced to move forward by the particles coming from behind.

NUMERICAL MODELS
Explosion of 6 kg spherical in shape trinitrotoluene (TNT) in free air was modeled using the explicit code LS-DYNA, with three different numerical models, pure Lagrangian model (ConWep), axisymmetric MMALE model, and 3D MMALE model.Blast wave parameters (pressure and specific impulse) in all models were tracked at each 0.1 m distance starting from 0.2 to 0.7 m.This corresponds to close range scaled distance 0.11 m/kg 1/3 < Z < 0.385 m/kg 1/3 , where Z=R/m 1/3 , "R" is the distance from the center of the charge to the target (stand-off distance) and "m" is the equivalent mass of TNT charge.

Pure Lagrangian model (ConWep)
The simplest numerical way of modeling and investigating design of components under blast loading is pure Lagrangian model.With this model only the component under investigation is modeled with Lagrangian element formulation and the ConWep function is called to apply the blast pressure over the elements.The ConWep function ( 2) is based on blast loading equations proposed by Kingery and Bulmash (1984), which are nothing else but empirical approximation of experimental results.
In this equation ! is the angle of incidence of the blast wave relative to the reflecting surface.This is numerically the cheapest method for evaluation of blast response of structures.Main drawback is that the method is unable to represent the interaction between the blast wave and the structure when complicated geometries are investigated.This means that the blast wave cannot be focused due to geometry and also the main structure cannot be "hidden" by the sacrificial one, if it is not physically connected to it.Another limitation of this method is that the empirical equations underlying the spherical air blast are valid for the range of scaled distance 0.147 m/kg 1/3 < Z < 40 m/kg 1/3 (LSTC (2013a)).

Arbitrary Language Euler (ALE) models
Lagrangian formulation becomes almost useless when very large deformation processes are simulated (Souli et al. (2000)), leading to a very small time steps and possible solution failure.For Latin American Journal of Solids and Structures 11 (2014) 1999-2017 such cases an alternative ALE method was developed in which the mesh does not necessary follows the material.In this method Lagrangian motion of the nodes is computed at every time step, followed by a possible advection stage in which the mesh is either not advected (pure Lagrangian formulation), advected to the original shape (pure Eulerian formulation), or advected to some more advantageous shape (somewhat between Lagrangian and Eulerian).In the current model the explosive charge and the air were modeled with MMALE formulation in which one element can possible contain both ALE materials.Advection between the steps in all MMALE models was controlled using the modified Van Leer advection algorithm.Structural material in the model was represented with Lagrangian formulation.The coupling between ALE fluids and Lagrangian plate was defined with penalty-based FSI algorithm which conserves energy.

Geometry and mesh
In order to study the influence of mesh size on the incident blast loading parameters the square box model with size of 0.9 m was represented with axisymmetric 2D and 3D models (Figure 2) with different mesh sizes of: 25, 50, 100, 200, 400, 800, and 1600 elements per side (corresponding to element lengths of: 36, 18, 9, 4.5, 2.25, 1.125, 0.5625 mm) for 2D axisymmetric model and 25, 50, and 100 elements per side for 3D model.The charge radius of 96.58 mm (6 kg of TNT) was represented with approximately 3, 6, 11, 22, 43, 86, and 172 elements with the finest regular mesh, leading to a massive models of over 2.5 million elements.Despite the fact that spherical mesh better represents the blast wave in free air (Kakogiannis et al. (2010)), rectangular mesh was used for all the models because it is better in solving the leakage problem when FSI is involved (Olovsson and Souli ), Swee et al. ( 2012)) (at least in the particular case of flat plate, placed in parallel with the mesh of the fluids).Spherical shape of the explosive was represented using the initial volume fraction geometry card available in LS-DYNA.The fluids in all 3D mod-els were represented with eight-node brick elements with one integration point and default viscous form of hourglass control, while the structural plate was modeled with eight-node fully integrated solid elements.In the axisymmetric MMALE model axisymmetric area weighted element formulation was used.Boundary conditions as depicted on Figure 2 were applied in all MMALE models.The reflected blast wave parameters were investigated for stand-off distance of 0.3 m (scaled distance Z=0.165 m/kg 1/3 ) using one-fourth symmetry 3D models.To examine in more detail the influence of DoB conditions on structural deformation, and to relate the minimum required distance at which the boundary conditions can be placed with some physical property of the model, two set of simulations were performed.First set was performed at each 0.1 m stand-off distance from 0.1 m to 0.7 m while keeping the mass of the charge (m = 6 kg) constant.In the second set of simulations the stand-off distance was kept constant (R = 0.3 m) for different mass values of the explosive charge (m = 0.4, 1, 8 and 24 kg).An example representing the series of models for case of R = 0.3 m and m = 8 kg is shown on Figure 3.
After the selection of the proper DoB conditions, simulations were performed with different mesh densities of: 20, 25, 50, and 100 elements per side (≈ 3, 4, 8, and 16 elements per radius of the charge) to study the mesh size influence.In each model, plate sides were represented with two Lagrangian elements per one Eulerian, while five elements through the thickness (t = 10 mm) were kept for each model.Mesh biasing influence was examined using 3D model of 25 elements per side and three different mesh biasing factors of 2, 4, and 6.Finally, the most optimal values for DoB conditions, mesh size and mesh biasing were used in structural response validation study.The explosive charge in all ALE models, was represented with *MAT_HIGH_EXPLOSIVE_BURN in combination with Jones-Wilkins-Lee (JWL) equation of state (EOS); which calculates the blast pressure as a function of relative volume != !!/!, and internal energy E, for an explosive element.In this equation A, B, R 1 , R 2 , and ! are parameters related to the explosive material and can be found in most of the explosive textbooks.They were taken from reference (Zukas and Walters (2002)) for TNT high explosive and are given in Table 1.

Air
The air is best approximated as an ideal gas for which *MAT_NULL was used in combination with linear polynomial equation of state; in which ! = 1.4 is the ratio of the specific heats, ! is current density, ρ !=1.29 kg/m 3 is initial density and E=250 kJ represents initial internal energy of the air at atmospheric pressure of 1 bar.

Structural material
Structural plate (diameter D = 0.5 m, plate thickness t = 10 mm) in all the models is represented with bilinear material model and parameters representing RHA steel were taken from Neuberger et al. (2007a).

Incident (Side on) blast wave parameters
The change of incident blast wave parameters (pressure and specific impulse) with change of mesh density, are presented on  Pressure contours for three different mesh sizes at the same computational time are shown on Figure 5.It can be observed that as the mesh is refined the blast wave is getting its spherical shape despite the rectangular shape of the mesh.It is also visible on the figures with coarser mesh (Figure 5a  mum values along the axes.As the mesh is refined, the pressure reaches nearly the same value in all directions (Figure 5c).This is more precisely shown on Figure 6.Comparison of the axisymmetric MM ALE model and 3D MM ALE model with the experimental results is shown on Figure 7.The results are here presented as a function of number of elements per radius of the charge.It is clear from Figure 7 that we need considerably refined mesh in order to reach experimentally obtained peak pressures while the values of specific impulses can be approached using coarser meshes (See Figure 7).However, in impulsive loading regime the specific impulse is responsible for structural deformation and not the peak pressure (Westine et al. (1975)), which allows us more efficient representation of even complicated 3D models on personal computers.Number of elements per radius of the charge needed to achieve a value of 90% of the exact numerical solution for both, incident pressure and impulse is presented on Figure 8 as a function of scaled distance.The values of the exact numerical solutions were estimated using the finest three solutions according to Grid Convergence Index (GCI) method given by Roache (1998).It can be observed from the same figure that depending on the scaled distance, various mesh densities are needed to achieve the 90% of the exact numerical solution.For the specific impulse 65 elements per charge radii were needed near the surface of the charge, while the number of elements decreased to 4 at the scale distance of 0.165 m/kg 1/3 and then stabilized at nearly 17 elements per charge radii for higher scaled distances.

Reflected blast wave parameters
Beside the material parameters and advection algorithm scheme, mainly the mesh size has the greatest influence on incident blast wave parameters.However, when some structure is investigated under blast loading, additional parameters like DoB conditions of ALE domain, FSI algorithm, and also the ratio between Lagrangian and ALE mesh sizes can have considerable influence on the reflected loading parameters, and consequently on the structural deformation.
-Influence of the DoB conditions On Figure 9 are shown the results for reflected pressure and reflected specific impulse profiles for the most loaded element located at the corner of the one-fourth symmetry models.It is clearly visible that the DoB conditions have great influence on blast loading parameters especially when they are too close to the object of interest (cases 0.325 and 0.4 m).As we move the boundary conditions away from the plate they are losing the influence and after 0.6 m box case we are getting the same results for pressure and specific impulse (See Figure 9).This appears because ambient pressure applied at the boundaries, influences the FSI and makes it unstable, when boundaries are close to the FSI location.-Mesh size influence With defined DoB (0.6 m box model) simulations were performed with different mesh sizes of: 20, 25, 50, and 100 elements per side (≈ 3, 4, 8, and 16 elements per radius of the charge) to study the mesh size influence on reflected blast wave parameters.The simulation results are presented on Figure 10.It is clear that reflected specific impulse approaches the experimental value as the mesh is refined.However, it should be noted here, that the curves are not wholly comparable, because the area of the corner element considerably differs among the models.

Structural deformation
Numerical factors that have influence on blast loading parameters will consequently influence the deformation of the loaded object.On Figure 12a, time-displacement curves are presented for the models with different DoB conditions from the plate.Similar to blast loading parameters, timedisplacement curves are changing until 0.6 m box case while no difference is visible with their further moving away (See Figure 12a).On Figure 12b, the influence of mesh density on structural displacement of the plate is presented.It can be observed here, that structural displacement curves are not in consistence with the blast loading parameters in this case.This is due to area differences of the centrally loaded elements.Maximum displacement is growing as the mesh is refined from 3 to 8 el./radius of the charge while no significant difference is visible with further mesh refinement.On Figure 12c) mesh biasing influence on structural deformation is presented.The results for minimum required distance from both sets of simulations (described in part 3.3.)are related to the radius of the charge and plotted on Figure 13 as a function of the scaled distance.It can be observed that the results from both sets are quite different, which confirms the complexity of numerous parameters that affect the minimum value of the DoB conditions.However, it can be concluded that DoB conditions higher than four times the radius of the charge, can assure a model in which boundary conditions have no influence on the structural response.

Structural response validation
Beside validation of the blast loading parameters, validation study of the structural response was also performed.For that purpose numerical simulations were done and the results compared with the experimental tests results published by Neuberger et al. (2007a).Three experimental cases with different close range stand-off distances (see Table 2) were simulated, and the results for maximum normalized deflection compared.Numerical models were prepared as explained in part 3.3, applying the findings of the study; MMALE materials were represented with a box of size R+4R TNT , meshed with 10 elements per radius of the charge, additionally refined toward the center of the charge with a biasing factor of 4. On Figure 14, numerical normalized deflections for all three cases are plotted versus time and corresponding experimental values are added for better comparison.

APPLICATION TO OTHER SCENARIOS
In reality, more complicated structures like sandwich plates, military vehicles and civil buildings are often analyzed under blast loading.So the decision where to put the boundary conditions and is it necessary to cover all the structure within MMALE can be challenging.To answer that question two different structures representing a sandwich plate (Figure 15a) and simplified vehicle chassis (Figure 16a) were analyzed under blast loading of 8kg TNT charge and constant stand-off distance 0.3m, while changing the DoB conditions.The results are presented on Figure 15b) and on Figure 16b).It can be observed that in both cases the minimum required DoB conditions is 0.4m.The difference is that the sandwich plate has to be all covered with MMALE materials, thus minimum DoB conditions should be measured from the upper plate (Figure 15a) while only a part of the structure can be covered in the case of simplified vehicle (or building structure) and the DoB conditions should be measured from the bottom plate of the vehicle chassis.This is due to the mutual proximity and also connectivity between different parts of particular structure.In the case of sandwich plate, both plates are locally deformed.In the case of blast response of simplified vehicle chassis, the bottom plate is locally deformed while the roof experiences only the global movement of the vehicle.

SUMMARY AND CONCLUSION
In this paper the influence of the numerical parameters on blast wave loading and structural deformation was investigated in detail.Numerical simulations with different mesh densities were performed for spherical TNT charge in free air.As a result the asymptotic values of incident pressure and specific impulse were obtained for different scaled distances and compared with the experimental results.From these results the optimal number of elements per charge radius was also determined for different scaled distances.The investigation was continued examining the influence of DoB conditions from the object of interest, mesh size, and mesh biasing on the reflected blast loading parameters and structural displacement.It was find out that combination of these parameters has great influence on the results as well as on the size and efficiency of the numerical model.The proper combination of these parameters can greatly increase the accuracy and decrease the computational time.According to the results obtained minimum mesh design criteria is proposed.Models with 10 el./radius of the charge, boundary conditions placed away from the deforming structure at a minimum distance of four charge radii (4R TNT ) additionally refined at the place of the charge with mesh biasing factor of 3 to 4 are recommended.This should yield reasonably accurate results in a relatively short computational time.

Figure 4 ,
for four different stand-off distances using axisymmetric MM-ALE model.For better comparison, corresponding experimental values taken from US Army defense design manual (UFC (2008)) and empirical ConWep data from the simulations are also added.It can be noticed that numerical values are approaching the experimental ones as the Latin American Journal of Solids and Structures 11 (2014) 1999-2017 mesh is refined.It can also be noticed that the empirical model in LS-DYNA overestimates blast loading parameters for scaled distances Z < 0.147 m/kg 1/3 .
and 5b) that pressure contours have maximum values along the diagonals and mini-Latin American Journal of Solids and Structures 11 (2014) 1999-2017

Figure 9 :
Figure 9: Influence of DoB conditions on the reflected loading parameters.

Figure 10 :
Figure 10: Mesh size influence on the reflected blast loading parameters.

Figure 11 :
Figure 11: Mesh biasing influence on the reflected loading parameters.

Figure 12 :
Figure 12: Influence of; a) DoB conditions b) mesh size and c) mesh biasing on plate displacement.

Figure 13 :
Figure 13: Minimum DoB vs charge radius as a function of scaled distance.

Figure 15 :
Figure 15: Application to other structures -Sandwich plate.

Figure 16 :
Figure 16: Application to other structures -Simplified (vehicle chassis or building).

Table 1 :
Material properties and JWL parameters of TNT.