A simplified approach to modelling blasts in computational fluid dynamics (CFD)

D.Mohotti ,K.Wijesooriy ,S.Weckert

a School of Engineering and Information Technology, The University of New South Wales, Canberra, ACT, 2600, Australia

b Sensors and Effectors Division, Defence Science and Technology Group, Edinburgh, SA, 5111, Australia

Keywords: Blast loads Computational fluid dynamics Explosions Numerical simulations

ABSTRACT This paper presents a time-efficient numerical approach to modelling high explosive (HE) blastwave propagation using Computational Fluid Dynamics (CFD).One of the main issues of using conventional CFD modelling in high explosive simulations is the ability to accurately define the initial blastwave properties that arise from the ignition and consequent explosion.Specialised codes often employ Jones-Wilkins-Lee (JWL) or similar equation of state (EOS) to simulate blasts.However,most available CFD codes are limited in terms of EOS modelling.They are restrictive to the Ideal Gas Law (IGL) for compressible flows,which is generally unsuitable for blast simulations.To this end,this paper presents a numerical approach to simulate blastwave propagation for any generic CFD code using the IGL EOS.A new method known as the Input Cavity Method (ICM) is defined where input conditions of the high explosives are given in the form of pressure,velocity and temperature time-history curves.These time history curves are input at a certain distance from the centre of the charge.It is shown that the ICM numerical method can accurately predict over-pressure and impulse time history at measured locations for the incident,reflective and complex multiple reflection scenarios with high numerical accuracy compared to experimental measurements.The ICM is compared to the Pressure Bubble Method(PBM),a common approach to replicating initial conditions for a high explosive in Finite Volume modelling.It is shown that the ICM outperforms the PBM on multiple fronts,such as peak values and overall overpressure curve shape.Finally,the paper also presents the importance of choosing an appropriate solver between the Pressure Based Solver (PBS) and Density-Based Solver (DBS) and provides the advantages and disadvantages of either choice.In general,it is shown that the PBS can resolve and capture the interactions of blastwaves to a higher degree of resolution than the DBS.This is achieved at a much higher computational cost,showing that the DBS is much preferred for quick turnarounds.

1.Introduction

The effects of blastwaves and their consequent damage to structures have been an increasingly popular research topic in the past decade [1-3].The blastwaves can arise from various sources ranging from deliberate explosions to accidental gas or vapour cloud explosives.Experimental methods have been at the forefront of evaluating the effects of blast blastwaves where Kingery and Bulmash [4]were the first to produce comprehensive analyses for air blast parameters undertaken for TNT explosives ranging from 1 to 4×106kg.The outcomes of Kingery and Bulmash are the basis of the UFC 3-340-02 [5]code of practice for the effects of blasts on structures.The code contains fundamental parameters recorded from the experiments relevant to structural design,such as incident and reflected pressure,incident and reflected impulse and time duration of the positive and negative phases of the blastwave.

With the rapid improvement of computational hardware over the last decade,numerical simulations in the form of Finite Element(FE)analysis for blastwaves have been widely investigated in recent literature.The main advantage of a numerical approach is the ability to perform multiple parametric studies,which would otherwise be expensive and cumbersome in an experimental setting.At the forefront of numerical simulations for blasts,Eulerian and Arbitrary Lagrangian-Eulerian (ALE) methods are the most prominent type of modelling undertaken via software such as LS-DYNA [6].ALE modelling is capable of incorporating both fluid and structural dynamics where the shockwave propagation and structural dynamics are coupled,resulting in complex fluidstructure interaction (FSI) capabilities.Thus,Eulerian methods have been proven to reliably predict blastwaves in numerical modelling in recent literature over the last few decades [7].

Computational Fluid Dynamics (CFD) is a numerical method that has recently gained traction in various fields of application for fluid flow problems.Most CFD codes employ the Finite Volume(FV)method based on the Navier-Stokes equation of flow to predict the fluid flow.With recent advancements in computational hardware over the past few decades,CFD is increasingly being used in many fields of study outside of its traditional scope.Recently,CFD has been actively used in Wind Engineering applications [8,9],previously dominated by experimental methods.With the availability of new computational hardware coupled with efficient modelling techniques,CFD has transitioned from a research-focused technique to a practical design approach.

Recently,CFD has also been employed in shockwave modelling of explosives in a few computational codes,where the advantages of the Finite Volume(FV)based analysis have been identified over other techniques used.Shockwave modelling is not a novelty in CFD,given that supersonic and hypersonic flows are commonly studied for external aircraft flows[10].However,modelling of blast and its propagation within the Finite Volume(FV)technique is not common.Commercial and research tools such as BlastFoam [11]and Viper:Blast [12]are available and have shown capabilities of modelling complex blast scenarios.One of the drawbacks of CFD modelling of blastwaves is that a limited amount of literature is only available compared to the more popular ALE or Eulerian methods.However,literature pertaining to blast modelling in CFD has recently been on an upsurge.One of the earliest and successful implementations of CFD in predicting blastwaves was conducted by Hansen et al.[13]where the software FLACS [14]was used.The study aimed to showcase the validity of the CFD approach,and to that end,several blast scenarios relating to gas explosions were modelled[13].Li et al.[15]presented a study where two CFD codes were used to predict internal and external pressures from a vented gas explosion.Their study used FLACS [14]for the reaction and pressure predictions for an internal area,whilst far-field external pressure predictions were made in FLUENT [16].It was noted that FLACS provided reasonable estimates of pressure time histories in the close vicinity of the explosion whilst providing relatively poor results in the far-field region,which was why FLUENT was implemented [15].The study showcased that the adopted method of using two CFD solvers provided reasonable estimates of overpressure at the near and far-field regions.Li et al.[17]also studied Boiling Liquid Expanding Vapour Explosion (BLEVE),where nearfield and far-field validations showed that CFD was capable of predicting peak pressures.This study was further enhanced where more BLEVE validations were conducted for complex interaction scenarios,and results were used to develop an analytical model[18].Molkov et al.[19]studied a complex scenario of a blast wave generated from a hydrogen tank rupture in a fire using CFD.Their study aimed to validate both pressure and fireball size and shape against experimental results,and it was shown that the combustion model could achieve this.In comparison,only a limited number of studies have been performed on high explosives in the application of CFD simulations [20-22].Ofengeim and Drikakis [20]applied CFD to simulate blast wave propagation over a cylinder where the effects of blast wave duration on pressure profile,viscous-inviscid interaction and loads on the cylinder were discussed.Chun [21]used CFD to simulate blast loads on a flexible tent for TNT explosives.To consider fluid-structural interaction (FSI) effects,Chun combined the CFD application with the deflection of the flexible tent to calculate blast wave propagation over the tent and blast loads on the tent.Finally,Ngo et al.[22]included one case study that used CFD to predict the blast effects of TNT on a building's facade.An 8-story building was modelled,and pressure distribution on the facade was analysed.The peak pressure in excess of the prescribed probability of failure was calculated to predict the failure of the building.From the literature above it can be seen that CFD can be used to model high explosive blasts.However,most of the recent studies conducted thus far are focused on gas type explosions or induced liquid expanding vapour explosions.Conversely,the studies conducted for high explosives are outdated and do not highlight the method by which the initial explosive is achieved.In the case of Chun[21]a moving wall that replicates the propagation of a shockwave front is modelled.The moving of the wall (which is stopped at some distance),creates the shock propagation through the domain.Ofengeim and Drikakis [20]uses a pressurised zone which is ruptured and in Ngo et al.[22]the details are obscure.

Gas explosions are usually modelled as real gasses where the chemical combustion processes are readily available in commercial codes[19,23],and only stoichiometry needs to be defined.For high explosives,equation of state(EOS)modelling is generally employed in Eulerian modelling,usually defined by the Jones-Wilkins-Lee(JWL) EOS.The development of these initial conditions is often not implemented in commercial CFD codes.Thus,researchers resort to other means such as using a moving a wall [20]or rupturing a pressurised zone [21]as discussed previously.This means that the inability to define initial input parameters correctly leads to poor predictions of the blastwave propagation throughout the domain.Thus,a reliable approach to input parameters for blasts resulting from explosives is necessary in CFD simulations.

The aim of this study is to investigate the capability of CFD simulations to model blastwaves within a domain and provide a complete framework for modelling.This includes a newly proposed input method for high explosive blast simulations intended to be consistent,repeatable,relatively easy to implement and applicable to any CFD code.In this study,the commercial code ANSYS FLUENT[16]is used where validation was performed based on three experimental case studies;(a)Case study 1-1 kg of Nitromethane where incident pressure was measured at distances 1.74,2.05 and 2.35 m(b)Case study 2-1 kg of Composition-B where incident and reflected pressure was measured at a distance 1.5 m,(C)Case study 3 -experiments performed by Caitlin et al.[24]where complex reflections of overpressure and impulse were measured.For the CFD analysis,two input methods are investigated for the blastwave creation within the domain.The first method is defined as the“Input Cavity”method introduced in this paper as a new modelling technique.The second method presents a refined version of the“Pressurised Bubble” method first introduced by Brode [25].The paper mainly investigates and presents the essential considerations that are required to be taken when modelling in CFD,such as the Equation of State modelling and choice of solver type.

2.Experimental setup and acquisition

The incident and reflected pressure data were obtained from three experimental tests which were utilised to validate the CFD simulations presented in this manuscript.The experiments all utilised single charge explosives where different charge configurations are considered in terms of type,weight,and standoff distance.The three experimental programs are given below,where the first two experiments were unique to this study,whilst the third was adopted from Caitlin et al.[24]as it captures complex multiple reflections of blastwaves.

Case study 1 presents the incident pressure measurement where the blast from a 1 kg Spherical Nitromethane charge is concerned,with a TNT equivalency of 1.0 [26].Non-sensitised Nitromethane was used as the high explosive material and was contained in a spherical casing with a diameter of 120 mm manufactured with 3D printed Nylon.The experimental setup is shown in the schematic in Fig.1,where the location of the explosive source and the pressure gauges are located to measure the incident pressure.Pressure gauges are placed and measured at distances of 1.74,2.05 and 2.35 m away from the blast source.Fig.2 also shows the ignition and initial gaseous fireball development of the explosives captured from a high-speed camera.

Fig.1.Test setup for Case A: 1 kg Nitromethane.

Fig.2.Detonation gases emanating from a free air burst of a spherical Nitromethane charge [2]: (a) Initial detonation;(b) Expansion of fireball gas;(c) Visible shockwave propagation;(d) Further progression of fireball and shockwave (Frame-rate=6,095 fps).

Case study 2 was performed to measure both the incident and the reflected pressures of a 1 kg Composition-B charge encapsulated in a spherical container of diameter 103.6 mm.The experimental configuration and setup are presented in Fig.3,where the charge is suspended from a steel frame above the target plate.The distance from the charge to the plate is 1.5 m,and the pressure sensors are embedded in the plate.The plate is placed on top of the blast rig,as shown.In addition,incident pressure was measured at 1.5,2.0 and 2.5 m from the Composition-B charge in a separate experiment (not displayed in Fig.3).

Fig.3.Test setup for Composition-B 1 kg for reflected pressure.

Case study 3,as mentioned previously,was performed by Caitlin et al.[24]where a single case was used from a series of experimental tests.The experimental setup(Fig.4)consists of three rows of concrete blocks (0.6×0.6 m2square cross-section).In Caitlin et al.[24],multiple configurations of block sizes and gap widths between rows were studied to understand the effects of energy reflections of blastwaves due to these changes in geometries.The aim of this manuscript does not coincide with that of Caitlin et al.[24]and instead would only choose one model case scenario to determine the capability of the proposed CFD method to accurately predict the complex blast propagation and reflections within the domain.Thus,Experimental Test 1 is considered by Caitlin et al.[24],which used a 0.324 kg PE-4 explosive (TNT equivalency of 0.46 kg),and the experimental layout is displayed as a schematic in Fig.4.Finally,as displayed,the location at which pressure is measured is between the first two blocks (Point “A”) and last two blocks (Point “B”).

Fig.4.Schematic view of the experimental setup Caitlin et al.[24].

3.Background for numerical modelling

The state of matter,also known as the Equation of State (EOS),plays a vital role in accurately determining the blastwave properties.In the LS-DYNA Eulerian modelling approach,the employment of the JWL EOS is a common approach for modelling explosives.The libraries of different parameters are provided for the different types of high explosives.However,in most generic commercial codes of CFD,the JWL EOS is not available,and the modification of the EOS can be restrictive.Thus,this section discusses the limitations of existing EOS modelling in typical commercial codes and provides a reliable alternate method of modelling explosives for any generic CFD code.The details of the approach are explained in the following sections,where an understanding of the JWL equation needs to be elaborated and compared with the Ideal Gas Law (IGL) to demonstrate the presented methods’ validity.Also,comparisons will be drawn to Eulerian modelling where applicable through the following sections.

3.1. Equation of state (EOS) modelling in CFD

Ideal Gas Law (IGL) is commonly used for modelling compressible flows in CFD.It covers a broad scope of flow scenarios,including the modelling of shockwaves[10]for supersonic[27]and hypersonic flows [28]of external aerodynamics.The IGL can be fundamentally written as given in Eq.(1),whereRis the ideal gas is constant andnis the amount of substance.The Pressure (P),Volume (v) and Temperature (T) of the gas and its behaviour are interrelated as that described in Eq.(1).

The assumption of IGL(that is,perfect gas conditions)does not hold true for scenarios of high pressure,temperature and density such as that encountered in a high explosive blast scenario.The detonation of chemical reactants occurs within microseconds and consequently gives rise to high pressures,temperature and density.As presented in Fig.5(left),there are two visible zones when a high explosive detonates;(1) the near field -the fireball region in the closest vicinity to the high explosive and (2) the far field -the region further away,which is air in atmospheric conditions.In this fireball region (as shown in Fig.5),the IGL relationship does not hold true between the variables of pressure,temperature and density (Eq.(1)).Therefore,to simulate the detonation and the resulting expanding detonation products in the area surrounding of the source,the Jones-Wilkins-Lee (JWL) EOS is often used as the EOS model [29].In conjunction,the parameters at the Chapman-Jouguet state are also required to model the detonation of the explosive[29].The explosive at the centre(Fig.5)and its behaviour can be defined by the JWL,including the surrounding expanded gases.The pressure and volume(density and mass)of the reactants and the behaviour of the gas are described as given in Eq.(2),whereA，B，ω，R1andR2are material specific constants,whileEis the detonation energy per unit volume andVis the relative volume or otherwiseV=ρsol/ρa(bǔ)irwhere ρsolis the density of the detonating material (example: 1630 kg/m3for TNT) and ρa(bǔ)iris the density of air.Parameters for TNT and Comp-B which were used in this study are presented in Table 4.

Table 4 JWL parameters used for 1D-Eulerian model.

Fig.5.Gaseous mediums of an explosion,on the left schematic on the right actual capturing the gas fireball for a nitromethane charge.

In most commercial CFD codes such as ANSYS FLUENT,the JWL EOS is unavailable.The available IGL EOS cannot model the initial explosion and the detonation products(fireball gas).Modelling this complex chemical reaction is outside the scope of this study,as mentioned in Section 1 and also disadvantageous as it can add to the complexity and computation cost of the overall simulation.The identification of two zones and the inability of the IGL to model the detonation products was identified.However,the IGL can still be used in the region further away from the fireball.The manuscript presents a method by which the IGL can be used to define the blastwave generated by the explosive.For this,the transition from JWL to IGL needs to be identified.This can be recognised by comparing the nuances of the two EOS explained in the following section.

3.2. Comparison of the IGL and JWL EOS

To compare the two EOSs,a plot of Pressure (P) vs Relative Volume (V) is displayed in Fig.6,where the respective curves are built based on Eqs.(1)and(2).The plots are constructed with TNT as the example explosive and the derivations necessary to construct Fig.6,from Eqs.(1) and (2) is provided in.Appendix A The parameters for the JWL TNT equation are well established and were obtained from Ref.[30].Figs.6(a) and 6(b) are the same plots but with an adjusted Pressure scale on theY-axes.Observing Fig.6(a),it can be said that there is a significant disparity between the JWL and IGL equations atV＜5.Decreasing Relative Volume on theX-axes indicates that ρa(bǔ)ir→ρTNTand this is expected as the density of air closer to the detonation product must be higher due to the fireball.The IGL for air at these low relative volumes do not capture the same Pressure variation as that of the JWL.This indicates that the IGL EOS cannot be appropriately used to capture the same Pressure-Density relationship closer to the explosive source;hence the disparity observed at lower Relative Volume values.Fig.6(b)shows the same graph,but with compressed Y-Axes values and as observed,it can be seen that at higher values ofV,the JWL EOS equation tends to merge with IGL EOS at aroundV＞10 (that is,ρa(bǔ)ir＜168 kg·m-3）.A closer observation of the JWL equation (Eq.(2)) clearly displays that the first and second term derivatives,at higher relative volumes,tend to zero,leaving only the third term.Upon closer inspection,the third term of the JWL equation is,in essence,the IGL equation which indicates that the behaviour of the gas further away from the source of the explosion follows the IGL properties of air.

Fig.6.Comparison of JWL (TNT) and Ideal Gas Law (IGL): (a) Full view and (b) Squished Y-axes.

To summarise,at the initial detonation,high pressure and temperature,which consequently results in high density,are experienced by the air(ρa(bǔ)ir)in the closest vicinity of the explosive.In the case of TNT (as the example used here),ρa(bǔ)ir→ρTNT,which indicates the lower relative volumes depicted in Fig.6 and the resulting high pressures.At low relative volumes(or inversely high ρa(bǔ)ir) it is clear that the IGL cannot make accurate pressure predictions.However,at the relative volume of 10 and higher(Fig.6(b)),IGL for air and TNT start to overlap.This visible transition of the JWL leading to IGL for air is an important feature that is exploited when defining the input method of the blastwave in the CFD code.

4.Numerical model development

This section explains the numerical modelling aspect in CFD and the methods by which the blastwave can be implemented with respect to the two different approaches described in this study.The first method is known as the Input Cavity method(ICM),a concept introduced by the authors in this study.The second method is known as the Pressurised Bubble method(PBM).Both methods will be explained in detail,with the main precedence given to the former,the input cavity method.

4.1. The input cavity method (ICM)

The Input Cavity method (ICM) can be visualised from the depiction in Fig.7,where two scenarios are shown.Fig.7(a)depicts a scenario where the numerical simulation models the detonation of the explosive whereby a blast pressure propagates within the domain.This straightforward approach is normally modelled,especially within a Eulerian domain.Fig.7(b) shows an empty cavity where the explosive and the resulting gas bubble(up to some distance) are located.This volume is taken out,and instead,as shown,an inlet is provided where time histories of Pressure,Velocity and Density/Temperature of the blastwave are defined,depicting the ICM.This procedure is explained in the succeeding paragraphs.

Fig.7.The input cavity method idealised: (a) A normal numerical simulation;(b) The cavity method.

This is a similar concept to what is performed in typical Eulerian modelling(in softwares such as LS-DYNA)where a 1D-model of the explosive is performed and then mapped to the 2D or 3D domain.The reason for such an approach is twofold.Firstly,modelling the explosion in 2D,and more so in 3D domain,is numerically demanding.Secondly,the 1D model is simpler and has a higher mesh resolution.This allows the prediction of the blast pressure(the range that provides high pressure with complex fire ball creation) from the source up to the mapping distance at a finer precision,which may not be computationally efficient in 2D or 3D setting as finer mesh requirements may ramp up the computational time.

The Input Cavity method for the CFD simulation bears a similar resemblance to the 1D mapping technique explained earlier.The procedure,however,is not straightforward in a generic CFD code such as ANSYS FLUENT as user interaction is required,and the placement of the inlet is critical due to the properties of the IGL as discussed in Section 3.The inlet properties for the CFD simulation are determined via a fast-running 1D-Eulerian model of the explosive,where at the pre-determined location of the inlet(for the CFD simulation),the time history curves of Pressure,Velocity and Density are recorded in the 1D simulation.

The location of the inlet for the CFD simulation must be selected with caution since the IGL EOS available in FLUENT cannot capture the complex gaseous behaviour near the explosive source,as explained in Section 3.1.The distance of the inlet,from the centre of the explosive charge,must be selected at a location beyond the radius of the detonation products fireball(Fig.5).That is,a location where the JWL EOS and IGL EOS coincide in the Pressure vs Relative Volume plot (Fig.6),as described in Section 3.This can be determined by plotting either the relative volume vs time or density vs time history curves at measured distances from the source of the charge.An example of this is given in Fig.8 where the density vs time curves measured at various distances from the source of 1 kg of TNT is displayed.Three curves at distances 25,50 and 75 cm from the centre of the charge were selected as probable locations outside the fireball region where the density time histories are presented.A few observations can be made where(1)the location of the peak of the curve shifts further away with recording distance which is selfexplanatory and expected.(2) The density peak is higher for the location recorded closer to the charge,which is once again expected as the air closer to the charge will experience higher temperatures.Observing the density curve at 25 cm it can be seen that the peak reaches a value of around 25 kg/m3which,when converted to relative volume corresponds to a value ofV=65.Observing Fig.6,it can be seen that the IGL overlaps JWL(for TNT)at relative volume values greater than ten and hence,creating a cavity with a 25 cm blast radius(or greater)for 1 kg of TNT where input parameters are extracted would be appropriate

Fig.8.Comparison of density vs time curves at 25,50 and 75 cm away from a 1 kg blast source.

Once the cavity and inlet size is established,it is important to import the correct input parameters to perform the simulation in the CFD domain.Three input components are necessary to describe the blastwave at the inlet of the CFD simulation: Static pressure,Velocity (dynamic pressure) and Density (or temperature).This study used the 1D-Eulerian model in LS-DYNA to generate the input data.In FLUENT,the variation of temperature must be input;thus,the density time history input needs to be converted to temperature.For this,the relationship of pressure,density and temperature(T) must be idealised in the Ideal Gas Law (IGL) equation.The classical IGL equation can be rearranged such that the temperature is a function of the absolute pressure and density,as displayed in Eq.(3).WhereRis the specific gas constant and takes a value of 287 J·kg-1·K-1for air.

Thus,in summary,for theInput Cavity Method,three input parameters are required for the conservation of momentum in ANSYS FLUENT and are ideally described via a Velocity inlet:

I.Static pressure

II.Velocity (or Dynamic pressure)

III.Temperature profile

A sample of the Static Pressure,Velocity and Temperature profile is presented in Fig.9 for Nitromethane 1 kg measured at an input cavity location of 1.74m,which is used in results Section 6.1.To summarise the application of ICM,Table 1 presents a step-bystep procedure to determine the parameters and input profiles required to detonate 1 kg of Nitromethane.

Fig.9.Example of defined inlet profiles required for ICM for a 1 kg Nitromethane charge with an inlet radius of 1.74 m as that presented in results,Section 6.1:(a)Static Pressure;(b)Velocity and (c) Temperature.

Table 1 Steps required to employ ICM.

4.2. The pressurised bubble method (PBM)

The Pressurised Bubble Method (PBM) used for this study employs the Isochoric conditions as described by Blanc et al.[32].The PBM(as shown in Fig.10)is centralised on the amount of energy a pressurised bubble releases when it “ruptures” and can be calculated as given in Eq.(4).WhereEbis the energy in the bubble,rbis the radius of the bubble(assuming it’s a sphere),ρbdensity of the pressurised air in the bubble,ebenergy of the bubble,ρ0density of the air in the surrounding vicinity (for example atmospheric conditions) ande0the energy of air (0.21 MJ/kg at atmospheric conditions).

Fig.10.Idealised Pressure Bubble method.

The energetic equivalence thus requires that

wherePbis the pressure in the bubble andP0is the pressure in the vicinity of the bubble (for example atmospheric conditions).

Observing Eqs.(4)or(5)it is clear that three parameters(eb，ρb，&rborPb)are unknown and it is required to be chosen to solve the problem where the pressurised balloon has TNT equivalency.The only requirement is such thatPb＞P0and it must be stated that the balloon size doesn’t need to be the same size as the TNT charge used [32].The user has the choice to either choose the size of the balloon and then calculate the resulting pressure or define the pressure in the balloon and then calculate the consequent volume/size based on Eq.(4).The Isochoric approach suggests that the energy in the pressure bubble（eb) is equal to the energy of TNT(eTNT) and maintain Isochoric conditions within the bubble.The density of the bubble(ρb)needs to be calculated in accordance with Eq.(6).The adiabatic constant,γ is taken as 1.4 (for air) as the pressurised bubble constitutes of air,as mentioned in Blanc et al.[32]).

In FLUENT,the pressurised bubble method is relatively easy to implement where the domain is split into two regions as idealised in Fig.10.The volume where the gas is pressurised needs to be patched with high pressure and temperature until the density inside the bubble is achieved to that what is acquired in Eq.(6).A usual mistake made by users is not adjusting this density and only applying pressure to an appropriate volume as calculated by Eqs.(4) and (5).This is often referred to as Isothermal conditions [32]which leads to poorer results.

Important consideration needs to be made on the pressurised bubble method,which works only for 3D domains as the energy stored (volume of the bubble) is defined in terms of volumetric calculations.To apply the pressure bubble method for the 2D domain,anAxisymmetricboundary condition is used.In FLUENT the Axisymmetric solver needs to be enabled and would solve the 2D Axisymmetric equations of flow instead of the 2D Cartesian form Ref.[33].The effect of axisymmetric boundary conditions can be better visualised using Fig.11.Fig.11(a) shows the 2D domain modelled in FLUENT with the axisymmetric boundary condition.Fig.11(b) shows the idealisation of applying the axisymmetric boundary condition and how FLUENT interprets.If the boundary was defined as symmetry instead(which is a common mistake),the domain is improperly defined,causing results that are not relevant to the idealised conditions.

Fig.11.Pressure Bubble method example of 2D case with axisymmetric boundary: (a) Real model;(b) Idealised model when using axisymmetric boundary.

Table 2 Calculation steps for 1 kg of Comp-B to apply PBM method.

To summarise the application of PBM Table 2 presents a step-bystep procedure to determine the parameters required to detonate 1 kg of Composition-B charge.The bubble pressure(Pb)needs to be pre-defined by the user as shown in Step 2 with a requirement such thatPb＞P0.

5.Numerical modelling settings

This chapter presents the numerical modelling settings and simulation steps employed in FLUENT for both the Input Cavity method and the Pressurised Bubble method.Both methods follow common general settings with minor differences mainly due to the difference in initial conditions used to generate the blastwave as explained in Section 4.

5.1. Global settings

All simulations were solved in transient conditions where transient time stepping was selected in ANSYS FLUENT.Also,the energy equation was turned on as the flow is compressible and energy conservation is essential in simulations where flow speeds exceed Mach 0.3.When modelling turbulence,the inviscid option was chosen.Flow generated following a blast event is in the supersonic region,and its high-pressure interaction with objects renders the viscous effects insignificant to the simulation's overall flow.Air material properties were selected and represented as ideal-gas conditions including the equation of state.For pressure measurements,“Static Pressure” needs to be selected and residual convergence was set to the default of 10-3.Residual convergence of 10-3was deemed sufficient based on the accuracy of incident pressure and impulse curves obtained.Also,employing a tighter residual of convergence would require a smaller time step,thus prolonging the simulation with diminishing returns to accuracy.Finally,the time step size was chosen based onCFL＜0.5 where the solution converged within transient iterations.

5.2. Pressure Based Simulations (PBS) and Density Based Simulations (DBS)

The Pressure Based solver has always been used for incompressible flows which are usuallyM＜0.3 and situations where the pressure is not sensitive to density changes and/or temperature changes.The Density Based solver is mainly reserved for highly compressible flows where shockwaves are concerned and are often employed to solve such problems [33].The theoretical difference between the two solvers is down to the method by which the flow equations are solved.ThePressure Basedsolver (PBS) has a segregated solver and a coupled solver of which,thePISOsegregated solver is recommended when using the Pressure Based solver for a compressible flow.In cases where convergence is harder to achieve,the COUPLED solver is recommended as it simultaneously solves mass and momentum equations in a coupled manner.Whilst it offers better convergence,the nature of coupled algorithm leads to longer solving time.TheDensity Basedsolver(DBS)solves Mass,Momentum,Energy and Species in a coupled manner from the outset.The only difference is the formulation of Implicit and Explicit simulation methods.In general,Implicit is preferred where time-step control is required whilst Explicit is used for high Mach-Flows with shockwaves.

In recent iterations of ANSYS FLUENT,through development,both types of solvers are capable of handling simulations outside the scope of its original implementation.Thus,the choice of solver is not as straight forward as expected and the recommendation given by FLUENT states that the solver has to be chosen specific to the problem at hand[33].Given that both solvers are applicable to be used in a wide range of flow simulation types,the choice of the solver from the user’s perspectives then mainly depends on the solution time required and the detail of flow features resolved.From the onset of this work the Density Based solver was used with explicit time stepping as the default solver.This was chosen because of the nature of blast simulation where the Density Based solver is preferred.Explicit is chosen because the time-step size can be controlled based on the Courant number (CFL) which requires minimal user input and aCFL＜0.5 is used to ensure stability in the simulation.A comparison between the PBS and DBS will be made in the results section.Also,the difference in settings used for the DBS and PBS are highlighted in Table 3 as shown below.

5.3. Boundary conditions

Modelling the domain for open air free blast tests can be idealised as presented in Fig.12 which shows the sequence of steps taken to model for a 2D domain.Fig.12(a)presents the complete 3D domain model and Fig.12(b) shows 1/8th of the domain where three surfaces are marked as symmetrical.This leads to computational time savings and also allows for better mesh refinement.The 2D model (Fig.12(c)) is a plane cut across the 1/8th sphere model and has two symmetries a blast source and outlet.In total,five boundary conditions are used for the FLUENT simulations.

Fig.12.Modelling CFD domain: (a) Idealised full model;(b) 3D model which 1/8th of a sphere and (c) 2D model which is a cut plane of the 1/8th.

I.A velocity inlet is necessary where the Input Cavity method is opted and is used to describe the input variables of momentum and thermal quantities.For this,the static pressure,velocity and temperature profiles are used.

II.The Pressure Outlet defines the location at which pressure exit is defined and usually needs to be situated reasonably far from the location of interest.The location of the pressure outlet must be such that it’s not too far away to cause the domain to be unusually large but not close enough to have reverse flow effects.

III.The wall definition is used for surfaces where reflections are anticipated.Because the flow is described as inviscid,wall effects on flow are not calculated.The wall boundary by default,has infinite rigidity and has to be used in locations where pressure reflections are measured.

IV.The Symmetry boundary is used when the physical geometry of interest,and the expected pattern of the flow/thermal solution,has mirror symmetry.However,caution needs to be taken when using symmetry such that if,flow is not parallel to the boundary it may act as a wall and reflect thus causing simulation errors within the domain.That is,both conditions of geometry and flow symmetry need to be maintained.A symmetric boundary condition cannot be defined if the geometry is only symmetric whilst the flow condition/pattern isn’t.

V.Axis is used to define axisymmetric boundary condition for 2D flows where 3D idealisation is required as explained previously.

VI.A pressure inlet(Optional)can be opted to be used instead of a velocity inlet however this requires a static pressure,temperature and dynamic pressure as opposed to velocity profile as that required for a velocity inlet.The dynamic pressure can be calculated but is one extra calculation step which can be avoided if used the velocity inlet.

Table 3 FLUENT settings for simulations.

6.Results

This section displays the results for test cases outlined in Section 2 where (a) free air blasts for Nitromethane and Composition-B is measured,(b) reflected pressure is measured for Composition-B and (c) multiple reflections of pressure and impulse are measured for PE-4.All simulations are 2D and were performed on a standard desktop computer with 4 cores and wherever the time for a simulation is shown,it is based on clock minutes or hours.These validation studies are followed by parametric studies where comparisons of the(a)Input Cavity Method(ICM)and Pressure Bubble Method(PBM)are made and the(b)choice of solver(PBS or DBS)is critically evaluated.Table 4 summarises the JWL parameters for TNT and Composition-B which were used for the 1D-Eulerian models in LS-DYNA,to generate the inlet profiles for the ICM in CFD.For the multiple reflections test where PE-4 is used,the TNT equivalency provided by Caitlin et al.[24]is used.All models presented are evaluated on a fine mesh setting where further refinement has diminishing returns on solution accuracy.

6.1. Free air blast tests for nitromethane (1 kg) and Composition-B(1 kg)

The overpressures for the free air blast experimental tests for Nitromethane were measured at distances of 1.74,2.05 and 2.35 m as mentioned in Section 2.For the CFD analysis,the input cavity method was used and the location of the inlet was selected to be at 1.74 m from the centre of the charge.This approach was opted because the experimental results were available at this point in order to validate the 1D-Eulerian model from which the input data is obtained for the CFD model.The 1D-Eulerian overpressure model results at 1.74 m is displayed in Fig.13(a)and can be seen to provide good correlation with the experimental outcome.The overpressure,temperature and velocity curves were obtained from the 1D-Eulerian model for the input into the CFD domain at 1.74 m.The size of the 2D domain for this particular analysis was 6 m in radius with a cut-out for the velocity inlet at 1.74 m as described previously.A total of 150×103quad cells were employed for the simulation with a maximum skewness of 0.1 which is a high quality mesh,as displayed in Fig.13(b).Also as observed,the mesh stretches from fine to coarse radially in the domain thus concentrating the finer cells closer to the inlet.TheDensity Basedsolver(DBS)was used and the simulation only took ten clock minutes to run.

Fig.13.(a) Input curve from 1D-Eulerian simulation and (b) typical 2D mesh for air blast CFD simulation.

The overpressure time history and impulse time history plots are presented in Figs.14(a) and 14(b) at 2.05 m and 2.35 m respectively for 1 kg of Nitromethane.As observed at 2.05 m(Fig.14(a)) the peak pressure and the decay for the curve look identical for the numerical and experimental (EXP) results where the difference in peak pressures is less than 5%.Also,the decay of the curve shows very good correlation between the experimental and numerical results.Some discrepancies in the negative phase of pressure is captured in the numerical simulation.This perhaps could be attributed to the inability of modelling the combustion of gases in the 1D-Eulerian method,where the output of this simulation is used as the input for the CFD ICM.Similarly observing the impulse time history,it can be seen that at 2.05 m the numerical curve can predict the peak value.Past the point of roughly 3.2 ms,the numerical curve follows a downwards trend in comparison to the experimental.This difference is mainly attributed to the two positive peaks recorded later in the experimental record (secondary and tertiary shocks).These are not captured in the numerical simulation due to the 1D-ALE mapping as the input condition(these are 3D effects from over-expansion and subsequent contraction of the detonation products).At 2.35 m,the numerical model presents a larger impulse reading in comparison to the experimental,especially peak values.This can be attributed towards the differences observed in the slope of the overpressure curve once it’s achieved its peak.Over time,the numerical curve falls below the experimental curve and this is because of the positive phase pressure registered in the experimental results between 5 and 7 ms.In general,it can be said that CFD simulation technique in the 2D domain was in good agreement with the experimental results and presents capability in predicting both overpressure and impulse time histories to a good accuracy.The capability of the CFD simulation is its ability to propagate the shock/blast wave within the domain.The accuracy of the CFD model is highly dependent on the input curves used and observing Fig.13(a)a difference in curve smoothness is observed.Furthermore,the effects of signal noise must also be taken into account.

Fig.14.Overpressure/Impulse vs Time graphs for 1 kg Nitromethane at (a) 2.05 m and (b) 2.35 m.

Fig.15.Overpressure/Impulse vs time graphs for 1 kg Composition-B at (a) 1.5 m and(b) 2.0 m and (c) 2.5 m.

Modelling of the Composition-B 1 kg explosive was conducted similar to that for the Nitromethane case where the only difference was the inlet conditions obtained from a 1D-Eulerian model.The location of the inlet was 1 m from the blast source where input curves are derived from the 1D-Eulerian model.The graphs for the overpressure inlet curve and the mesh aren’t shown here for brevity as they are similar to those presented in Fig.13 for the Composition-B explosive.The overpressure and impulse time history plots are presented in Figs.15(a)-15(c)at a distance of 1.5,2.0 and 2.5 m from the source of the blast.As observed at all distances,the overpressure predictions closely correlate with the experimental results in terms of peak and decay of the curve at all measured points.The minor deviations observed in the decay curve at 2.0 m and 2.5 m can consequently affect the impulse time history and this is observed in Figs.15(b) and 15(c).However,in general,the impulse time history curves show good correlation with the experimental results thus showcasing the ability of the CFD numerical simulation to be capable of simulating the effects of blastwaves to a high degree of accuracy.

6.2. Single charge reflection for Composition-B and comparison of ICM and PBM

The reflected pressure for Composition-B (1 kg) explosive was measured at a distance of 1.5 m as mentioned in Section 2,where the setup is shown in Fig.3.To replicate this in the numerical domain in a 2D manner the domain was first idealised as a 3D case and then progressively simplified as shown from Figs.16(a)to 16(c).This model consisted of high-quality hexahedral mesh (maximum skewness of 0.5)and is as displayed in Fig.16(d)where the reflected pressure is measured at the location marked next to the intersection of the wall and symmetry boundary (Fig.16(d)).For the ICM method,the inlet was provided at 1mdistance similar to that described in Section 6.1 for the free air blast tests.For the PBM method,a pressurised bubble with radius 0.37 m is patched with a pressure of 1×107Pa and temperature of 7400 K was used to simulate a 1 kg Composition-B charge as described in Table 2.The numerical domain and its configuration are as depicted in Fig.11 where the patched zone is clearly visible,and the reflection wall is in the direction of the propagating wave.The axisymmetric boundary is defined to simulate the 3D idealisation as that described in Section 4.2.

Fig.16.Domain creation for numerical simulation: (a) Isometric idealisation;(b) Plan view and quarter model idealisation;(c) Side elevation view of quarter model and (d) side elevation view of model and mesh configuration for quarter model.

The reflected pressure and impulse vs time curves recorded at 1.5 m from the centre of the charge are displayed separately in Fig.17(a) and (b) respectively where comparisons between the experimental,PBM and ICM are made.First observations clearly indicate the magnitude of the reflected pressure where at a similar distance,the free airblast tests displays a peak of roughly 400 kPa(Fig.15(a)) whereas,the reflected pressure is roughly four times this value.In comparison to the experimental (EXP),the ICM and PBM methods show good correlations of pressure over time curves(Fig.17(a)) albeit the latter slightly over-predicted peak pressures.In comparison,the ICM shows closer resemblance to the experimental curve in terms of decay and observation of second peak.This is further reflected in the impulse time history curves(Fig.17(b))where it clearly shows that the ICM is more comparable to the experimental result.In general,the Pressure Bubble Method(PBM) can be used as an alternative and can be found useful as it doesn’t require additional input parameters as that required for the Input Cavity Method (ICM).However,the ICM is superior when it comes to accuracy of the numerical results.

Fig.17.Comparison of reflective quantities of ICM and PBM for (a) pressure and (b) impulse.

The process of a blastwave propagation from the inlet towards the reflective boundary,which is the steel plate on a table (in the experiment),and its reflection can be illustrated through Figs.18(a)to 18(d).As observed the blastwave initially has a high incident pressure Figs.18(a) and 18(b) but slowly diminishes as it reaches further outwards as displayed in Fig.18(c).At the wall reflection occurs (Fig.18(d)) as shown by the marked arrow.The coloured contours display the rise in pressure once again thus explaining the high magnitude reflected pressure at 1.5 m(1800 kPa)compared to the incident pressure(396 kPa in Fig.15(a)).This occurs due to the non-linear superposition of the reflected pressure and the approaching blastwave.

Fig.18.Blastwave propagation: (a) Start of simulation;(b) Blastwave propagation;(c) Blastwave propagation further out loosing magnitude and (d) reflected pressure gaining magnitude as marked.

6.3. Multiple reflections and comparison of Pressure Based Solver(PBS) and Density Based Solver (DBS)

The test setup and experimental procedure performed by Caitlin et al.[24]were described in Section 2 where a 0.324 kg PE-4(0.46 kg TNT) charge was detonated and pressures recorded over a sequence of obstructions.The locations at which pressure was measured were in between the obstructions marked point A and point B in Fig.4.The domain setup and mesh is presented in Fig.19 where boundary conditions are labelled.The Input Cavity Method(ICM)is used and a mesh size of 10 mm is employed with a total of 180,000 hexahedral elements.A 1D-Eulerian model was used to generate the inlet profiles for 0.46 kg of TNT at a distance of 0.75 m from the blast source.For purposes of comparison,both Pressure Based Simulation(PBS)and Density Based Simulation(DBS)results are compared.

The overpressure and impulse time histories at locations A and B are presented in Figs.20(a)to 20(d)respectively.At point A,it can be seen that there are four distinct peaks in the overpressure curve before it tends to the negative phase.This trend is seen for both experimental and numerical results albeit some differences between the PBS and DBS methods.Comparing the overpressure numerical curves (Fig.20(a)) it can be seen that the PBS curve shows good correlation with the experimental.The magnitude of the peaks for overpressure and the location of peaks found in the positive phase are in good agreement.Furthermore,the decay of the curve after the fourth peak into the negative phase is well captured along with the final peak which is recorded in the negative phase.In comparison,the DBS overpressure curve also show peaks at the same locations with some minor differences in terms of peak values.The largest deviation is the slope towards the negative phase and the final peak where over-prediction is observed.These minor discrepancies are reflected in the impulse time history curve (Fig.20(b)) where the DBS curve over-predicts significantly.In comparison,the PBS curve shows better correlation.At point B (Fig.20(c)),the experimental and numerical overpressure curves once again show four peaks but not as distinct as that observed in point A.This is self-explanatory as the first incident wave approaching point B should be lower in magnitude as it is further away.There is however a significant spike recorded in the third and fourth peaks,as this is the reflected wave which incidentally have larger magnitudes due to reasons explained in Section 6.2.Comparing the numerical curves,it is clear that both PBS and DBS show similar first peak magnitudes but show larger responses in the second,third and fourth peaks.In general,the PBS curve predicts closer to the experimental where decay of the overpressure curve and negative phase is concerned.This is confirmed in the impulse time history curve(Fig.20(d))where the PBS impulse time history follows the same trajectory as the DBS but then falls closer in line with the experimental once the peak impulse is reached.

Fig.20.Comparison of multiple reflections for 0.324 kg of PE-4: (a) Overpressure at point A;(b) Impulse at point A;(c) Overpressure at point B and (d) Impulse at point B.

The formation of the peaks can be described by observing the contour plots of the blast wave propagation as that shown in Fig.21 where the left column is PBS and the right DBS.Figs.21(a)to 21(d)displays the time stamps of the multiple reflections that occur in this simulation due to the obstructions in the domain and the location of measurement is as marked by a“spot”.Fig.21(a)shows the approaching peak (peak 1),and this is a self-explanatory observation as it is the approaching flow of the initial detonation/blastwave.The second peak in Fig.21(b)shows that it is a result of the immediate ground reflection.Note that the first initial peak has still not reached the second block before this happens.The third peak as shown in Fig.21(c)is due to the reflection from the second block and explains the time duration between the second and third peak.Also,it can be visualised that the approaching multiple pressure peaks in Fig.21(c) which explains the fourth subsequent peak is caused by the reflection between the second block and ground interface.Finally,Fig.21(d)shows the domain reflections at a time where negative pressure is recorded at the point of measurement.

Fig.21.Comparison of PBS (Left) vs DBS (Right) for HSE case study.Blastwave propagation,referencing point A: (a) First approaching peak;(b) Second approaching peak,from ground;(c) Third peak approaching from reflection of second block and (d) negative phase.

One of the striking differences between the PBS and DBS are the formation of the peaks in Fig.20.It can be observed that the PBS overpressure curves shows much refined peaks whilst the DBS curves are “blunt” (Figs.20(a) and 20(c)).These differences are better reflected in the contour plots as presented in Fig.21 where it can be seen,at every given time step presented,the PBS simulation has better resolved blastwaves in comparison to the DBS method.It must be noted that both methods used an identical mesh setup.It is clear that as multiple reflections occur and blastwaves interact,the PBS is capable of resolving and capturing the interactions of blastwaves to a higher degree of resolution in comparison to the DBS.This explains the “bluntness” of the curve observed in Figs.20(a) and 20(c) for the DBS.In general,the Pressure Based solver (PBS) has always been used for incompressible flows which are usuallyMach＜0.3 and situations where the pressure is not sensitive to density changes and or temperature changes.The Density Based solver is mainly reserved for highly compressible flows where blastwaves are concerned and are often employed to solve such problems [33].In recent iterations of ANSYS FLUENT,through development,both types of solvers are capable of handling simulations outside the scope of its original implementation.Thus,the choice of solver is not as straight forward as expected and the recommendation given by FLUENT states that the solver has to be chosen specific to the problem at hand[33].Given that both solvers are applicable to be used in a wide range of flow simulation types,the choice of the solver from the user’s perspectives then mainly depends on the solution time required and the detail of flow features resolved.Given that explosions and blastwaves occur at a higher Mach number such that the flow can be considered supersonic,the DBS solver,in theory,is the more suitable of the two types.However,as presented here the PBS solver has been capable of providing higher resolution of the blastwave interactions.Another important point to consider is the solving time required for the two different solvers for the same mesh density.The PBS solver took 30 clock minutes whilst the DBS solver only took a total of 5 clock minutes in comparison.The reason for such a difference is the way the PBS solver runs where in each time step,the flow needs to be resolved and converge through iterations.This drastically increases the solving time required in comparison to the DBS.Also,it was noted that convergence in the PBS to be much harder and required a stringent time step size in comparison to the DBS.This is expected for the PBS as at higher velocity flows,where rapid change of density is encountered,the PBS scheme takes longer to converge which is why ANSYS FLUENT [33,34]recommends the use of DBS for supersonic flows.Thus,the recommended Solver is based on the user’s choice of accuracy vs computational time.Where results are required within a reasonable time frame,the DBS is more suitable and vice versa for the PBS.

7.Conclusions

This manuscript investigated the capability of simulating the effects of blastwave propagation using a Finite Volume (FV)approach for any Computational Fluid Dynamics(CFD)code.In this study,ANSYS FLUENT was used to explore the capability of CFD for this scenario.The study in particular presents a comprehensive numerical modelling technique for CFD and also includes important theoretical background which is required to understand the intricacies of the CFD approach.The modelling techniques presented for the CFD code were all validated using experimental data where cases of both free airblast and groundblast involving reflections are presented.Some of the key findings are summarised below:

· The presented Input Cavity Method(ICM)provides a simplified approach of defining input parameters for the blast.Instead of modelling the reactants and explosion,data is input as a time history function of static pressure,velocity and temperature where all three components are required for conservation of energy.The implementation of Input Cavity Method provided comparable results in both predicting over-pressure and impulse time histories of the blastwave propagation for free air blasts,single reflection and complex multiple reflection scenarios as demonstrated.

· The manuscript also discusses the theoretical background of the ideal gas law (IGL) and its applicability to be used in blastwave modelling.It was discussed that the IGL can be successfully employed given that the relative volume criteria is maintained and for TNT,V?10 needs to be considered to successfully employ the ICM.

· A comparison between the Density Based Solver (DBS) and Pressure Based Solver(PBS)was presented.It was shown that in general the PBS was better at resolving the flow field and in general showed better predictions in terms of over-pressure and impulse time histories.The DBS solver on the other hand loses resolution but is much faster in computational time in comparison to the PBS.Therefore,it was deduced that the choice of the solver used depends on the users’ requirements.

· A comparison between the Pressure Bubble method(PBM) and ICM was made and it was shown that using the correct modelling steps the PBM approach can also be used and gave relatively good outcomes in terms of over-pressure and impulse time histories.Some discrepancies were observed in the decay rate but it was noted that the PBM will be a useful method as it doesn’t require another input mechanism to generate the blastwave.

Finally,the proposed Input Cavity Method (ICM) presented in this study can be used as a consistent and reliable approach for blastwave propagation of high explosives when using CFD to conduct numerical modelling.This study fills a void in existing literature where initial conditions of high explosives now can be replicated in a CFD domain with consistency.It is shown that with the correct input criteria,high numerical accuracy can be obtained and for scenarios which include complex multiple reflections,the CFD approach is beneficial as outcomes can be obtained at a rapid pace.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A

Appendix A provides the derivations necessary for the IGL EOS and JWL EOS to plot Fig.6 as that presented in Section 3.2.The original form of the IGL is rearranged and expressed in terms of internal energy (Uint) which can be expressed as follows:

given that volume can be re-arranged in terms of density and mass Eq.(A2) can be re-written as follows:

where ρ andmare the density and mass respectively.Considering that air,which mainly constitutes of nitrogen and oxygen,which are diatomic gases,the value ofcan be assumed to befor monoatomic gasses[35]).Sinceis defined as a constant per unit mass the final equation can be re-written to formulate as that given in Eq.(A4).

Finally with the definition ofV=ρtnt/ρa(bǔ)ir,Eq.(A4) can be rearranged to Eq.(A5) where pressure is defined as a function of relative volume.The IGL EOS curve was constructed with the assumption that internal energy of air is 0.204 MJ/kg.A plot ofPvsVcan then be obtained as presented in Fig.6.

For the JWL EOS obtaining the curve is relatively straight forward.The JWL equation is reproduced below(Eq.(A6))where a plot of pressure vs relative volume can be conducted.It must be noted that at higher relative volumes the first and second terms tend to zero and is left with the third term only as presented in Eq.(A7)

where all the terms have been previously defined in Section 3.2 andEis the internal energy(can also be defined asUint).As mentioned by Shin [36],in software such as Ansys AUTODYN or LS-DYNA the compression ratio (μ) is calculated which is defined as Eq.(A8).

When the compression ratio reaches -0.99 or falls below this number,the density of TNT (ρa(bǔ)ir) is automatically set to 1 kg/ m3thus Eq.(A7) can be rewritten as Eq.(A9) where

where ω=γ-1,γ is the adiabatic constant and for air is generally taken as 1.3 -1.4,Eis the internal energy (as previously mentioned same asUint).Eq.(A9)is essentially the same equation as that derived in Eq.(A4) and Eq.(A5).