A coupled finite difference and material point method for simulation of gas–solid flow of thermobaric explosions

A coupled finite difference and material point method (CFDMP) was presented to simulate the dynamic of the gas–solid flow of thermobaric explosions. The coupled method takes advantages of both the two methods where the finite difference method captures the wave structures of gas phase and the material point method tracks the movement of solid phase. Particle collision is easily realized in MPM framework. By using a uniform grid system, the grid nodes have both the two phases’ information and their interactions are directly calculated. The simulation of thermobaric explosions in different containers is carried out and the results validate the accuracy of CFDMP program. The particle collision reduce the compression degree of the solid phase in the initial expansion stage but slightly influence the pressure history at the gauge. The process of particle movement in the container is demonstrate and the wall-reflected wave enhances the burning rate of particles.


Introduction
Thermobaric Explosives (TBX), commonly consisting of high explosives (HE) and metal powder (MP), have been widely applied in attacking targets hidden in caves and tunnels due to its high destructive ability compared with conventional explosives. The detonation of the HE creates a blast wave dispersing the MP into the ambient air and forms a gas-solid reactive flow. Numerical simulation is a powerful tool to demonstrate the evolution process of the thermobaric explosion field and some significant works have been done by Khasainov, Kuhl and Balakrishnan et al. [1][2][3]. They presented a two-phase Eulerian continuum model which combines the gas-dynamic conservation laws for the gas phase with the dilute continuum conservation laws for the solid phase. Inter-phase mass, momentum and energy interactions are prescribed by the empirically-based phenomenological relations [4]. Good agreements between numerical and experimental results have been demonstrated [5].
The solution of the two phase's equations is carried out in Eulerian-Eulerian framework. For the stability of the calculation, a phenomenological constitutive model [3,6] of the dispersed solid phase needs to be introduced to keep the equations hyperbolic. However, the drastically changing sound speed of the solid still presents a challenge to numerical methods. Moreover, The trajectory of particles is difficult to be described in Eulerian framework, which is important for describing the evolution of gassolid flow field. Therefore, the Eulerian-Lagrangian framework is another available way, such as particle-fluid method [7,8] and CFD-DEM [9]. Nevertheless, the continuity of the solid phase's displacement field is difficult to be guaranteed and the collisions between particles will cause a lot of EMCEME 2020 IOP Conf. Series: Earth and Environmental Science 692 (2021) 022091 IOP Publishing doi: 10.1088/1755-1315/692/2/022091 2 computation. The material point method (MPM) shows good computational properties in dealing with the problems about granular materials [10], but the coupling scheme of CFD and MPM applied to twophase flow still needs to be studied.
A coupled finite difference and material point method (CFDMP) for simulating the gas-solid flow of thermobaric explosions is presented in this paper. Firstly, we make a brief review of the two-phase model. After that, the solution schemes of two phases' conservation equations solved by FDM and MPM are described respectively. The collision of solid phase's particles is considered in MPM framework. The calculation of inter-phase interactions is realized by a uniform grid system of CFDMP. Then, the confined explosion experiments carried out by Kuhl et al. [2] are numerically investigated to validate the CFDMP. Finally, the effects of particle collision and container's size on the burning rate are discussed.

Conservation laws
In the two-phase model [2], the gas phase contains the detonation products, the ambient air, and their combustion products. The limit of large Reynolds and Peclet numbers in explosion situation make the effects of molecular diffusion and heat conduction negligible. The dynamics conservation equations for the gas phase are summarized as follows: where ρ, p, and u represent the gas density, pressure, and velocity vector respectively, / 2 E e = + ⋅ u u denotes the total energy of the gas phase and e denotes the specific energy. Q  represents the energy release rate due to reaction. i F and i b  denote the mass fractions and the density change rate of the component i. s σ , s f  and s q  represent burning rate, drag force and heat exchange respectively.
The solid phase refers to the metal powder cloud and is taken as an Eulerian continuum [11]. For the thermobaric explosion field, the volume fraction of the solid phase is very small in most situation and the collisions between the particles could be ignored. Neglecting pressure, the conservation laws for solid phase are summarized as follows: where σ and v represent the solid phase's density and velocity vectors. / 2 s s E e = + ⋅ v v denotes the total energy and s s s e C T = denotes the heat energy. Cs and Ts represent the material's specific heat and the solid phase's temperature respectively. Because the particles are incompressible, the specific energy associated with pressure disappears.

Interactions
where ρs and ds denote the particle density and diameter respectively, σB is the Stefan-Boltzmann constant, ε is the emissivity, and λg is the thermal conductivity of the gas. The ignition temperature Tign is specified as the melting point of aluminum (932 K) [2]. Res, CD, and Nu represent the Reynolds number, drag force coefficient, and Nusselt number of the particles, given by: where cp and μ denote the specific heat at a constant pressure and the dynamic viscosity of the gas phase respectively. The particle burning time ts is described as: 2 2 , 150sec/ cm

Combustion
In Kuhl et al. [12], the field consists of five components: air, two fuels (detonation products of the HE and active particles), and their combustion products. The specific internal energy of each component is fitted with a quadratic function of temperature, and the afterburning combustion process is illustrated by a Le Chatelier diagram. However, there may be some difficulties with this model when the components of the detonation products (such as CO, CO2, and H2O) react with the metal particles. Therefore, it is more common to consider combustion between the more basic components.
Taking the aluminum (Al) as a typical metal, we consider the reaction occurs between follow basic components: C, CO, CO2, H2, H2O, Al, Al2O3, O2. In this paper, the infinite rate model is used. The model means that when the components involved in the combustion process enter a computational cell they are completely consumed in stoichiometric proportions in a time step. We set n1, n2, n3…n8 represent the mole densities of the components before the combustion. Generally, the components are not in the state of chemical equilibrium due to the convection function and the reaction occurs. The following reactions are considered: The mole densities after the combustion are respectively denoted by m1, m2, m3…m8. So the energy release rate and density change rate are calculated as ( ) where M i and ∆H i represent the relative molecular weight and formation enthalpy of the component i, the values can refer to the JANAF table.

Equation of state
The equation of state (EOS) of the detonation products can be described by Jones-Wilkins-Lee (JWL) equation [13], which is where V represents relative specific volume, ρ 0 is the initial density of explosive and e is the specific internal energy. A, B, C, R 1 , R 2 and ω are the parameters. As the volume goes up, Eq.(16) degenerate into 0 e p V ωρ = (17) which is ideal gas equation and can describes the behavior of the ambient gas. It should be noticed that the combustion reaction proceeds in the condition that the detonation products expand and mix with the ambient air fully. So the combustion products are mostly at low pressure and the equation of state could be also described by Eq (17). However, the parameter ω should be decided by its composition and temperature. Considering the parameter ω has a relationship with the specific heat ratio γ 1 ω γ = − , the EOS of combustion products can be expressed as Due to the JWL parameters' different functions, we can combine Eq. (16) and Eq. (18) and get Where A, B, R1 and R2 are relevant to the detonation products and constant. ω is calculated according to Eq. (18). In this way, the behaviors of the detonation products, the air, and the combustion products are described uniformly.

Wave propagation method for gas phase
The conservation equations of the gas phase are solved in FDM framework. Splitting the source terms, Eqs. (1)-(4) are summarized as follows: For simplicity, only one-dimensional case is discussed here. The discretization of the variables is based on a Godunov's scheme, in which the cell average value n The wave propagation method is applied and illustrated in Figure 1. The cell averaged value 1 n j + Q are calculated by the resulting waves S k from the cell interfaces x j-1/2 and x j+1/2 , expressed as where λ k is the sound speed of the resulting wave and W k is the discontinuity value of the adjacent solution areas which include The limitations of the wave speeds λ -=min(λ, 0) at x j-1/2 and λ + =min(λ, 0) at x j+1/2 mean that the information propagates along characteristic paths at the correct wave speeds. Obviously, the discretization belongs to a class of upwind schemes. The wave speeds and the jumps can be calculated by the Roe approximate Riemann solver in the original version [14][15][16], which can be easily replaced by the Harten-Lax-van Leer-Contact (HLLC) scheme to ensure the flux is nonnegative [17]. The multidimensional case and the high order version of Eq. (23) can refer to the paper [14].
where m p and x sp denote the mass and the spatial coordinates of the material point p respectively, and N p is the total number of material points.
Due to the convection term of velocity, Eq. (24) cannot be discretized directly in the MPM framework. Thus, we introduce the variable η which represents the relative mass at current moment and ranges from 0 to 1. Then, the mass equation (24) is transformed into another equivalent form, expressed as: So the MPM framework of the solid phase is achieved. It can be found that the gradient of the shape function is unnecessary. Therefore, the standard MPM is more appropriate than the generalized interpolation material point method (GIMP) due to its computational efficiency and the unit decomposition condition on the boundary.

Particle collision algorithm of the material point method.
According to Eq. (40), the movement of the solid phase is only associated with the drag force and does not consider the influence of collisions between the particles. However, the collisions affect the configuration of the solid phase during the initial dispersal stage [3]. The contact algorithm proposed by Bardenhagen et al. [10] for compressible granular materials has higher computing efficiency than discrete element method but will cause a loss of kinetic energy for the solid phase discussed here. Our goal is to develop a new collision algorithm applied in incompressible particles and the kinetic energy of the particles is conserved after the collision.
In this algorithm, each material point is an independent object and represents a parcel of particles. Hence, the collision between particles is treated as the collision between material points. Each node has a point-mapping velocity v Ip , which is expressed as: It can be found that when the node has other associated points v Ip is different from the single-value velocity v I which is calculated as Eq. (31). The algorithm's conception is illustrated in Figure 2. The collision will happen in the node when v Ip and v I satisfy Where n Ip is the unit outward normal vector of the point and calculated by the gradient of the shape function N Ip , Viscous: After the elastic collision, the total kinetic energy c I K is calculated by:

Numerical results and discussion
Kuhl et al. [2,3,5] have carried out a series of thermobaric explosion experiments as well as numerical studies using a 0.5g spherical pentaerythritol tetranitrate (PETN) booster surrounded by an 1.0g aluminum (Al) powder in a cylindrical container, as shown in Figure 5. Two containers with different sizes are considered here. 'Container A' has a 20cm diameter and a 21cm length and 'Container B' has a 30cm diameter and a 30cm length. A gauge point is set on top of the container 5cm away from the center and the pressure histories are used to validate the accuracy of CFDMP. Ignoring the effect of gravity, the shape of the container and the location of the charge shown in Figure 5 indicate that this problem has both axisymmetric and plane-symmetric properties. So a quarter numerical model was applied, as shown in Figure 6. The computational grid is filled with the detonation products of PETN and ambient air. Meanwhile, the Al particles are represented by material points. The grid is uniform with a refined size of 0.625mm and each cell includes four material points. The initial temperature of Al powder is 273K and the mean particle diameter is set as 4.3μm [18]. The parameters of the JWL EOS of the detonation products are given in Table 1. The specific heat capacities of the gas components are fitted with a linear function of temperature Cv(T) = a + bT and the fitting coefficients are given in Table 2. Other physical parameters can be found in the literature [2]. The comparisons of the pressure and impulse histories in different containers are shown in Figure 7. For the containers, the arrival time of the first shock wave and the peak overpressure from the calculation are close to those from the experiment. The multiple subsequent peaks are also captured. There are some Although the particle collision algorithm is applied in the calculation of the solid phase, the original MPM automatically guarantees no penetration between particles. Thus, the influence of the collision algorithm is investigated. The comparison is shown in Figure 8. It is found that the distributions of the particles are different and the material points appear to aggregate without the collision algorithm, which may have affected the particle combustion process.  The influence of the particle collision on the pressure and the combustion is shown in Figure 9. While without collision, the aggregation in high-temperature area means that more particles satisfy the critical ignition condition, thus the peak pressure and the total combustion mass increase. The movement process of the particles in container A is shown in Figure 10. The clouds of CO 2 mass fraction are displayed to distinguish the detonation products and the ambient air. The black spots represented the solid particles. As the figure shows, the particles are immersed and dispersed by the detonation products in the initial stage. Then, the particles exceed the border of the detonation products and enter into the air due to the difference in the physical properties of the two phases. After that, under the function of the wall-reflected waves the particles converge into the center and begin to mix with the detonation products again. This process could generally last for several milliseconds.  Figure 11. It can be found that in initial state (< 90μs) both the combustion processes of particles are the same. After the moment, the burning rate of particles in 'Container A' exceeds that in 'Container B'. The situation is related to the container size. The gas phase's pressure distribution and the solid phase's configuration in Container A at 90μs are described in Figure 12. It can be seen that the shock wave reflects on the container's boundary and interacts with the particle. This process promotes the mixing of the particles with the gas mixture and enhances the combustion rate.

Conclusion
A coupled finite difference and material point method is presented to simulate the dynamic of the gassolid flow of thermobaric explosions. The conservation equations of the two phases are respectively solved by FDM and MPM, which takes full advantages of the two methods. The collision algorithm described in MPM framework can easily deal with the particle-particle interaction. After interpolating the material points' variables to the uniform grid, the interactions between the two phases are directly calculated in the nodes. Numerical simulations of PETN/Al explosion experiments are carried out and the similar pressure curves are obtained by CFDMP program. The particle collision reduce the compression degree of the solid phase in the initial expansion stage but slightly influence the pressure history at the gauge. The evolution curves of Al mass in different containers with the pressure distribution demonstrate that the wall-reflected wave enhances the burning rate of particles.