An Energetic Model for Detonation of Granulated Solid Propellants

Abstract

Unexpected detonation of granular solid energetic materials is a key safety issue in the propellants manufacturing industry. In this work, a model developed for the characterization of the early stages of the detonation process of granular solid energetic materials is presented. The model relies on a two-phase approach which considers the conservation equations of mass, momentum, and energy and constitutive relations for mass generation, gas-solid particle interaction, interphase heat transfer, and particle-particle stress. The work considers an extension of approximated Riemann solvers and Total Variation Diminishing (TVD) schemes to the solid phase for the numerical integration of the problem. The results obtained with this model show a good agreement with data available in the literature and confirm the potential of the numerical schemes applied to this type of model. The results also permit to assess the effectiveness of different numerical schemes to predict the early stages of this transient combustion process.

1. Introduction

Unexpected detonation of granular solid propellants is a key safety issue. Gun tubes experience pressure build-up during unstable burning in the chamber leading to high mechanical stresses which act on the grains leading to an increase of the chamber pressure, causing the explosion of the gun tube [1]. Moreover, ignition and subsequent combustion of granular solid propellants are also complex mechanisms of high importance in the manufacturing process of many energetic materials, not only because of already mentioned safety issues, but also due to the design of propulsion systems. The penetration of the hot gases into the voids promotes the convective heating of the granules as well as the granule compaction and the rapid pressurization of the system. In case the ignition is not correctly controlled, the process can generate compressive waves which compact the un-reacting region of the grain bed. Under strong pressure waves, the propellant bed might rapidly collapse, leading to a compressive reaction and an explosion. Therefore, it is of major importance to properly understand and model this problem.

The theoretical modeling of this problem has focused on the interests of many studies in the last decades and different approaches can be found in the literature. Some of the first researchers facing the theoretical modeling of this kind of problem were Kuo et al. [2] who developed in their work a fixed-bed model, usually referred to as the KVS (Kuo, Vichenevetsky, Summerfield) model, to describe the transient problem of flame propagation in a packed fixed-bed of granular energetic materials. They used a combination of an explicit finite difference scheme and the method of characteristics to numerically solve the problems. Koo and Kuo [3] compared the model predictions with experimental firing tests with satisfactory results.

On the other hand, Markatos [4] proposed a different approach based on the previous work of Spalding [5]. In that case, the description of the problem considered the governing equations on the basis that mass, momentum, and energy fluxes are balanced over control volumes occupied by space-sharing interspersed continua. According to that concept, distinct phases are present within the same space, their shares of space being measured by their volume fractions. The approach was used to simulate the accelerating flame front for beds of granulated solid energetic materials fixed in a rigid enclosure with satisfactory results [4].

Another type of approach is that developed by Gough [6] and adopted by Oton-Martinez et al. [7] among others, which is based on the formal averaging of the solid propellant and the surrounding gas in the spatial domain considered. This approach relies on the conservation equations of the solid and gas phases where the phases are considered to be a homogeneous continuum. This means that both the gas and the propellant grains are defined to have a uniform value of density and porosity within each fundamental volume cell. Thus, the heterogeneous flow composed of two interacting continua can be described by appropriately defined averages of flow properties.

Regarding the experimental data to feed numerical models, Sandusky [8] provided constitutive relations of compaction obtained from quasi-static compaction experiments with bed diameters approaching that of gun breeches performed on five cannon propellants. The granular materials tested were single-base, double-base, triple-base propellants as well as a nitramine-composite propellant. They found that the dominant compaction mechanism was plastic deformation with the virtual absence of fracture. Kooker et al. [9] experimentally studied the convective ignition process of granular solid propellant beds with the help of a previously built apparatus that transformed a confined bed to a planar hot gas ignition wave. The investigation examined convective ignition behavior of five different granular solid propellants with the same flow conditions: single-base, double-base, triple-base, and nitramine composites.

As for the numerical methods used in the modeling of detonation of granular solid energetic materials, most of the authors used finite-difference schemes [10] although with distinct approaches. Hoffman and Krier [11] and Krier and Kezerle [12] used the MacCormack method [13] with a predictor-corrector scheme. Sheu et al. [14,15,16] and Krier et al. [17] also used a finite-different scheme with the Lax–Wendroff approach [18] whereas Baer and Nunziato [19] and Butler [20] applied the method of lines [21,22]. Markatos [4] used an implicit finite difference algorithm called IPSA (interphase slip algorithm) which anticipates the effects of a change in the local property of one phase on the properties of the other phase at the same location.

To the author’s knowledge, finite-volume approximate Riemann solvers [23] or Total Variation Diminishing (TVD) schemes [24] have not been previously applied to the resolution of problems describing the detonation process of the granular solid energetic materials or explosives. Notwithstanding, these numerical methods are widely used in other types of fluid mechanic problems with good results. Among others, approximate Riemann solvers are heavily used in transient transonic or supersonic problems as well as in the case of transient contact discontinuity problems [25,26], whereas TVD schemes have shown good results in incompressible fluid flow problems [27].

In this work, a predictive model is developed as a tool for the characterization of the detonation process of granular solid energetic materials under shock tube conditions. The tool relies on a two-phase model, which considers a set of equations of continuity, momentum, and energy for the solid and the gas phase, to describe the transient behavior of this combustion process. Regarding the numerical methods, the tool considers the use of approximate Riemann solvers and TVD schemes to face the resolution of the system of partial differential equations of the model. The article is structured as follows: firstly, the physical model and the considered constitutive relations are presented. In the next section, the numerical schemes used for the integration of the model are described. Afterwards, the results provided by the model are compared with data available in the literature and the validation of the model is discussed. The effectiveness of different numerical schemes when facing the description of the early stages of the detonation of granular solid energetic materials is also discussed. Finally, some conclusions are drawn, and future research lines are outlined.

2. Physical Model Description

The approach proposed in this work relies on the conservation equations of mass, energy, and momentum of both gas phase and solid phase to model the physical problem of combustion of granulated solid energetic materials. In this case, the solid phase is considered as incompressible. In the present work, the model is generalized in order to be applied to different solid granulated energetic materials. The basic guidelines proposed by Hoffman and Krier [11] together with the considerations described by Gidaspow [28] are included in this model. Gidaspow [28] suggests including the time derivative of the porosity into the energy equation in order to satisfy the inequality of Clausius-Duhem. However, all space and time derivatives of the porosity, which are present in the momentum and energy equations, respectively, are non-conservative terms that cannot be solved with an explicit numerical method for hyperbolic systems. Thus, in the present work, an alternative model is presented where some constitutive laws are modified, and a re-organization of the conservation equations is performed to obtain a conservative form of the system of equations. That allows the use of numerical schemes for conservative fluxes with a finite-volume approach. This way, updated numerical methods for compressible flow can be applied for the numerical resolution [7]. Moreover, an extension of approximated Riemann solvers to the solid phase was considered for the numerical integration of the problem. Another key aspect in these models is the inclusion of a pressure into the solid, which is defined as the addition of the gas pressure and the inter-granular stress, which introduces an important stiffness from the numerical point of view. All in all, the differential equation system considered in this work can be expressed as:

$$\frac{\partial \left(\alpha {\rho }_g\right)}{\partial t}+\overrightarrow{\nabla }\left(\alpha {\rho }_g{\overrightarrow{u}}_g\right)=\dot{\Gamma },$$

(1)

$$\frac{\partial \left(\alpha {\rho }_g{\overrightarrow{u}}_g\right)}{\partial t}+\overrightarrow{\nabla }\left(\alpha {\rho }_g{\overrightarrow{u}}_g\otimes {\overrightarrow{u}}_g+\alpha p_g\overrightarrow{I}\right)=p_g\overrightarrow{\nabla }\alpha +\dot{\Gamma }{\overrightarrow{u}}_p-{\overrightarrow{F}}_D,$$

(2)

$$\begin{array}{c}\frac{\partial \left(\alpha {\rho }_g{E}_g\right)}{\partial t}+\overrightarrow{\nabla }\left(\alpha {\rho }_g{\overrightarrow{u}}_g\left(E_g+\frac{p_g}{{\rho }_g}\right)\right)=-p_g\frac{\partial \alpha }{\partial t}+\dot{\Gamma }\left(E_{chem}^g+\frac{{\left|{\overrightarrow{u}}_p\right|}^2}{2}\right)-{\overrightarrow{u}}_p{\overrightarrow{F}}_D-\dot{Q}\end{array},$$

(3)

$$\frac{\partial \left(\left(1-\alpha \right){\rho }_p\right)}{\partial t}+\overrightarrow{\nabla }\left(\left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\right)=-\dot{\Gamma },$$

(4)

$$\begin{array}{c}\frac{\partial \left(\left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\right)}{\partial t}+\overrightarrow{\nabla }\left(\left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\otimes {\overrightarrow{u}}_p+\left(1-\alpha \right)p_p\overrightarrow{I}\right)=-p_p\overrightarrow{\nabla }\alpha -\dot{\Gamma }{\overrightarrow{u}}_p+{\overrightarrow{F}}_D\end{array},$$

(5)

$$\begin{array}{rr}\hfill \frac{\partial \left(\left(1-\alpha \right){\rho }_p{E}_p\right)}{\partial t}+\overrightarrow{\nabla }\left(\left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\left(E_p+\frac{p_p}{{\rho }_p}\right)\right)& \hfill =p_g\frac{\partial \alpha }{\partial t}+\dot{\Gamma }\left(E_{chem}^p+\frac{{\left|{\overrightarrow{u}}_p\right|}^2}{2}\right)+{\overrightarrow{u}}_p{\overrightarrow{F}}_D+\dot{Q}\end{array},$$

(6)

$$p_p=p_g+{\tau }_p,$$

(7)

where $\rho$ is the density, $\alpha$ is the gas porosity, $\overrightarrow{u}$ is the velocity vector, $\overrightarrow{I}$ is the identity vector, $\dot{\Gamma }$ is the mass generation, $p$ is the pressure, ${\tau }_p$ is the intergranular stress, ${\overrightarrow{F}}_D$ is the drag force, and $\dot{Q}$ is the interfacial heat transfer. The chemical energy of the gas phase is defined as,

$$E_{chem}^g=E_{chem}-E_{chem}^p,$$

(8)

where $g$ and $p$ stand for the gas and the particle, respectively.

According to Hoffman and Krier [11], intergranular stress formulation does not prevent particle compaction below minimum compaction. To consider that limitation, they proposed a modification of the momentum equation of particles, so they cannot be accelerated once they have reached that porosity value. If that limitation is imposed, the Equation (5) can be rewritten to:

$$\begin{array}{c}\frac{\partial \left(\left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\right)}{\partial t}+\overrightarrow{\nabla }\left(\left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\otimes {\overrightarrow{u}}_p+\left(1-\alpha \right)p_p\overrightarrow{I}\left(1-f\left(\alpha \right)\right)\right)=\\ -p_p\frac{\partial \left(\alpha \left(1-f\left(\alpha \right)\right)\right)}{\partial x}-\dot{\Gamma }{\overrightarrow{u}}_p+\left(1-f\left(\alpha \right)\right){\overrightarrow{F}}_D\end{array}.$$

(9)

The term $f\left(\alpha \right)$ is defined as,

$$f\left(\alpha \right)=\left\{\begin{array}{ll}\begin{array}{c}0\\ {\left(\frac{{\alpha }_c-\alpha }{{\alpha }_c-{\alpha }_{min}}\right)}^l\\ 1\end{array}\hfill & \begin{array}{c}\alpha >{\alpha }_c\\ {\alpha }_{min}\le \alpha \le {\alpha }_c\\ \alpha <{\alpha }_{min}\end{array}\hfill \end{array}\right.,$$

(10)

where $l$ is a constant value that considers the properties of the particle’s framework, ${\alpha }_{min}$ and ${\alpha }_c$ are the minimum and critical porosity respectively.

Several correlations were chosen to close the proposed model. Among them, drag force, used in momentum and energy conservation equations of gas and solid phases, is defined as:

$$F_D=\frac{{\mu }_g\left(u_g-u_p\right)}{d_p^2}{f}_{pg},$$

(11)

where $d_p$ is the particle diameter, $f_{pg}$ is the drag coefficient presented by Butler et al. [29] in their study,

$$\begin{array}{c}{f}_{pg}=\frac{{(1-\alpha )}^2}{{\alpha }^2}150+3.89{\left(\frac{R{e}_p}{\left(1-\alpha \right)}\right)}^{0.88}\\ {10}^3<R{e}_p<2\cdot {10}^5\end{array}.$$

(12)

This correlation (Equation (11)) is a function of the Reynolds number $(R{e}_p)$ and gas viscosity $\left({\mu }_g\right)$:

$${\mu }_g={\mu }_{g,0}{(T_g/T_{g,0})}^{0.65},$$

(13)

where $T$ is the temperature, ${\mu }_{g,0}$ is the value of viscosity at 2000 K, and $T_{g,0}$ is 2000 K [12].

As for the interfacial heat transfer coefficient, the heat expression form, in the case of having spherical particles $\left(S_p/V_p=6/d_p\right)$, is:

$$\dot{Q}=\frac{6\left(1-\alpha \right)}{d_p}{h}_c\left(T_g-T_p\right).$$

(14)

The convective heat transfer coefficient $\left(h_c\right)$ depends on the Nusselt number ($Nu$) with the following expression:

$$h_c=\frac{k_g}{d_p}Nu.$$

(15)

The convective heat transfer coefficient is a function of the Prandtl Number, $Pr=\frac{c_{p,g}{\mu }_g}{k_g}$, and the gas conductivity $\left(k_g\right)$ defined as:

$$k_g=\frac{15}{4}\frac{R_u}{M_g}{\mu }_g\left(\frac{4}{15}\frac{c_{v,g}{M}_g}{R_u}+\frac{3}{5}\right),$$

(16)

where $c_{v,g}$ is the specific heat of the gas at constant volume, $R_u$ is the universal constant of gases, and $M_g$ is the molecular weight of the gas phase.

The correlation of the Nusselt number proposed by Hoffman and Krier [11], Butler et al. [29], and Butler and Krier [20] has the following form:

$$Nu=0.65R{e}_p^{0.7}P{r}^{0.33}.$$

(17)

The gas mass generation equation assuming that particles are spheres is related to the combustion velocity and has the following form:

$$\dot{\Gamma }=6\left(1-\alpha \right)\frac{S_p}{d_p}\dot{r}{\rho }_p.$$

(18)

The combustion law is given by:

$$\dot{r}=a{p}_g^n,$$

(19)

where $a$ and $n$ are vales that depend on the propellant.

As for the gas equation of state, the equation of Noble-Abel is considered,

$$p_g=\frac{{\rho }_g\left({\gamma }_g-1\right)e_g}{1-{\rho }_g{\eta }_g}.$$

(20)

As previously said, in this model the solid is considered incompressible. However, as presented in Equation (7), there is a solid pressure defined as the addition of gas pressure and the intergranular stress (${\tau }_p$), understood as the resistance of the particles to be compacted under the critical porosity. In this study, the authors consider the expressions proposed in several works, such as the ones from Krier and Kezerle [12], Hoffman and Krier [11], and Gokhale and Krier [30]:

$${\tau }_p=\left\{\begin{array}{ll}\begin{array}{c}\frac{K}{\left(1-\alpha \right)}\left(\frac{1}{\left(1-{\alpha }_c\right)}-\frac{1}{\left(1-\alpha \right)}\right)\\ 0\end{array}\hfill & \begin{array}{c}\alpha \le {\alpha }_c\\ \alpha >{\alpha }_c& ,\end{array}\hfill \end{array}\right.$$

(21)

where K is the bulk modulus.

3. Numerical Method

The physical model includes non-conservative terms. In order to handle this issue, a reconversion of the system was performed in order to have a conservative form of the system of partial differential equations. In particular, the momentum and energy equations for both gas and solid phase are re-organized so that the terms containing partial derivatives of $\alpha$ are resolved separately by including them in the source terms. The final system of equations in a homogenous form (without source terms) could be presented in the conservative form and remains hyperbolic. This way, a finite volume approach with an approximate Riemann solver can be used to numerically solve the problem. In this case, an extension of approximated Riemann solvers to the solid phase was considered for the numerical integration of the problem. Therefore, the multi-dimensional Euler equations used in this work for modeling and describing the physical problem of combustion of energetic materials can be written with the following general conservative form,

$$\frac{\partial U}{\partial t}+\nabla ·\mathcal{H}=S\left(U\right),$$

(22)

where $U$ is the conserved variables vector, $\mathcal{H}=\left(F,G,H\right)$ is the tensor of fluxes and $S$ is the vector that includes the non-conservative terms and the source terms, which may or not include non-conservative terms and can be defined as:

$$U=\left[\begin{array}{c}{q}_1\\ q_2\\ \dots \\ q_m\end{array}\right],\hspace{0.17em}F=\left[\begin{array}{c}{f}_1\\ f_2\\ \dots \\ f_m\end{array}\right],\hspace{0.17em}G=\left[\begin{array}{c}g\\ g_2\\ \dots \\ g_m\end{array}\right],\hspace{0.17em}H=\left[\begin{array}{c}{h}_1\\ h_2\\ \dots \\ h_m\end{array}\right],\hspace{0.17em}S=\left[\begin{array}{c}{s}_1\\ s_2\\ \dots \\ s_m\end{array}\right].$$

(23)

In our case, the conserved variables vector U, the flux tensor $\mathcal{H}$ and the source terms vector are

$$\begin{gathered} U=\left(\begin{array}{c}\alpha {\rho }_g\\ \alpha {\rho }_g{\overrightarrow{u}}_g\\ \alpha {\rho }_g{E}_g\\ \left(1-\alpha \right){\rho }_p\\ \left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\\ \left(1-\alpha \right){\rho }_p{E}_p\end{array}\right),\mathcal{H}=\left(\begin{array}{c}\alpha {\rho }_g{\overrightarrow{u}}_g\\ \alpha {\rho }_g{\overrightarrow{u}}_g\otimes {\overrightarrow{u}}_g+\alpha p_g\overrightarrow{I}\\ \alpha {\rho }_g{\overrightarrow{u}}_g\left(E_g+\frac{p_g}{{\rho }_g}\right)\\ \left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\\ \left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\otimes {\overrightarrow{u}}_p+\left(1-\alpha \right)p_p\overrightarrow{I}\left(1-f\left(\alpha \right)\right)\\ \left(1-\alpha \right){\rho }_p{\overrightarrow{u}}_p\left(E_p+\frac{p_p}{{\rho }_p}\right)\end{array}\right), \\ S\left(U\right)=\left(\begin{array}{c}\dot{\Gamma }\\ p_g\frac{\partial \alpha }{\partial x}+\dot{\Gamma }{\overrightarrow{u}}_p-{\overrightarrow{F}}_D\\ -p_g\frac{\partial \alpha }{\partial t}+\dot{\Gamma }\left(E_{chem}^g+\frac{{\left|{\overrightarrow{u}}_p\right|}^2}{2}\right)-{\overrightarrow{u}}_p{\overrightarrow{F}}_D-\dot{Q}\\ -\dot{\Gamma }\\ -\left(1-f\left(\alpha \right)\right)p_p\frac{\partial \alpha }{\partial x}-\dot{\Gamma }{\overrightarrow{u}}_p+\left(1-f\left(\alpha \right)\right){\overrightarrow{F}}_D\\ p_g\frac{\partial \alpha }{\partial t}+\dot{\Gamma }\left(E_{chem}^p+\frac{{\left|{\overrightarrow{u}}_p\right|}^2}{2}\right)+{\overrightarrow{u}}_p{\overrightarrow{F}}_D+\dot{Q}\end{array}\right). \end{gathered}$$

(24)

The integration of the homogeneous part of the equation system (Equation (22)) in a control volume $\Omega$ yields

$$\frac{\partial }{\partial t}{{\displaystyle \iiint }}_{\Omega }^{}\mathit{U}d\Omega +{{\displaystyle \iint }}_A^{}\mathcal{H}\widehat{n}dA=0,$$

(25)

where A is the boundary of Ω and $\widehat{n}$ is the normal vector to surface A. Considering the first integral as a time-rate of change of the averaging of the conserved variables $\mathit{U}$ and the boundary A formed by N surfaces so that $A={\displaystyle \sum }_{s=1}^N{A}_s$, Equation (25) can be written as,

$$\frac{\partial \mathit{U}}{\partial t}+\frac{1}{\left|\Omega \right|}{\displaystyle \sum }_{s=1}^N{{\displaystyle \iint }}_{A_s}^{}\mathcal{H}\widehat{n}dA=0.$$

(26)

The time derivative of the conserved variables is defined as $\frac{U_j^{n+1}-U_j^n}{\Delta t}$. Therefore, making use of the Rotational property $\mathcal{H}\widehat{n}={\overline{T_s}}^{-1}\mathit{F}\left(\overline{T_s}U\right)$ and considering that the surface integral of the flux is approached by ${{\displaystyle \iint }}_{A_s}^{}\mathcal{H}\widehat{n}dA\approx {\overline{T_s}}^{-1}\mathit{F}\left(\overline{T_s}U\right)A_s$, a finite volume scheme for multiple dimensions in unstructured grids is obtained. Considering the source terms of the nonhomogeneous system,

$${\mathit{U}}_j^{n+1}={\mathit{U}}_j^n-\frac{\Delta t}{\left|\Omega \right|}{\displaystyle \sum }_{s=1}^N{\overline{T_s}}^{-1}\mathit{F}\left(\overline{T_s}\mathit{U}\right)A_s+\Delta t\mathit{S}\left({\mathit{U}}_j^n\right),$$

(27)

where $\overline{T_s}$ is the rotation matrix, ${\overline{T_s}}^{-1}$ its inverse, and $A_s$ is the area of the $s^{th}$ surface bounding volume Ω.

The time step to solve the equation system is calculated as,

$$\Delta t=CFL·min\left[min\left(\frac{\Delta x_j}{\left|{\overrightarrow{u}}_g\right|+c_g}\right),min\left(\frac{\Delta x_j}{\left|{\overrightarrow{u}}_p\right|+c_p}\right)\right],$$

(28)

where $CFL$ stands for a Courant–Friedichs–Lewy-like constant, $c_g$ is the speed of sound of the gas phase,

$$c_g=\sqrt{\frac{{\gamma }_g{p}_g}{{\rho }_g\left(1-\eta {\rho }_g\right)}},$$

(29)

and $c_p$ is the speed of sound of the solid phase,

$$c_p=\sqrt{\frac{K_0}{{\rho }_p}},$$

(30)

where $K_0$ is the initial bulk modulus, whose value was taken from Hoffman and Krier [11], and ${\rho }_p$ is the density of the propellant.

In order to evaluate the fluxes, several approximations such as the first-order upwind method of Godunov, Lax–Friedrichs scheme, and Lax–Wendroff scheme can be found in the existing literature.

These methods require the solution of the Riemann problem which can be highly demanding. Therefore, approximate Riemann solvers are used in this work. Among the numerical schemes found in literature, Rusanov and MacCormack numerical schemes were chosen to solve the problem of energetic material combustion.

3.1. Rusanov Scheme

Rusanov scheme [31] is a particular case of HLL (Harten Lax and van Leer) Riemann solver [23]. In this approach, an approximation for the intercell numerical flux is obtained directly. The approximation consists of finding, for each interface between two cells, the approximate numerical flow that will be an average value of the numerical flows computed for each cell in the previous step time, corrected with a factor. This factor is proportional to the difference of conservative variables of each cell in the previous step time and a positive speed value $S^{+}$, so the numerical flux at each interface is computed as:

$${\mathit{F}}_{j+1/2}^n=\frac{1}{2}\left({\mathit{F}}_L\left(\overline{T_s}{\mathit{U}}_L\right)+{\mathit{F}}_R\left(\overline{T_s}{\mathit{U}}_R\right)-S^{+}\left({\mathit{U}}_R-{\mathit{U}}_L\right)\right),$$

(31)

being ${\mathit{F}}_{j+1/2}^n$ the flux at the interface, ${\mathit{F}}_L$ the convective flow in the left cell in the previous time step, $F\left(\overline{T_s}{U}_L^n\right)$, $F_R$ the convective flow in the right cell, $F\left(\overline{T_s}{U}_R^n\right)$ and $S^{+}$ a positive value of speed calculated by,

$$S^{+}=max\left\{\left|u_{nL}-c_L\right|,\left|u_R-c_R\right|,\left|u_{nL}+c_L\right|,\left|u_R+c_R\right|\right\},$$

(32)

where R and L stand for right and left respectively.

3.2. MacCormack-TVD Scheme

The MacCormack numerical method of second-order in both time and space is presented to solve Navier–Stokes equations developed by Robert MacCormack [13]. According to Lee and Kim [32], the scheme can be written as a variation of the two steps Lax–Wendroff scheme [18], one predictor and one corrector. According to Liang et al. [27], the use of a total variation diminishing (TVD) of the MacCormack scheme will accurately reproduce all flow regimes. Thus, the solution of the system of Equation (22) is obtained by solving three 1-D problems in sequence as proposed in [21],

$$\frac{\partial U}{\partial t}+\frac{\partial \mathit{F}}{\partial x}=S$$

(33)

$$\frac{\partial U}{\partial t}+\frac{\partial \mathit{G}}{\partial y}=R,$$

(34)

$$\frac{\partial U}{\partial t}+\frac{\partial \mathit{H}}{\partial z}=T$$

(35)

The TVD-MacCormack scheme is used to solve consecutively the three 1D hyperbolic equations in each time step. Following [21] and taking Equation (33) as an example the discretization scheme is given by

$${\overline{\mathit{U}}}_j^p={\mathit{U}}_j^n-\frac{\Delta t}{\Delta x}\left({\mathit{F}}_j^n-{\mathit{F}}_{j-1}^n\right)+\Delta t\mathit{S}\left({\mathit{U}}_j^n\right),$$

(36)

$${\overline{\mathit{U}}}_j^c={\mathit{U}}_j^n-\frac{\Delta t}{\Delta x}\left({\mathit{F}}_{j+1}^p-{\mathit{F}}_j^p\right)+\Delta t\mathit{S}\left({\mathit{U}}_j^p\right),$$

(37)

$${\overline{U}}_j^{n+1}=\frac{1}{2}\left({\overline{\mathit{U}}}_j^p+{\overline{\mathit{U}}}_j^c\right)+\left[G\left(r_j^{+}\right)+G\left(r_{j+1}^{-}\right)\right]\Delta {\mathit{U}}_{j+1/2}^n-\left[G\left(r_{j-1}^{+}\right)+G\left(r_j^{-}\right)\right]\Delta {\mathit{U}}_{j-1/2}^n,$$

(38)

where the superscripts p and c denote the predictor and corrector steps respectively, ∆x and ∆t are the spatial and time steps, and $G\left(x\right)$ is a function defined as,

$$G\left(x\right)=0.5C\left[1-\phi \left(x\right)\right],$$

(39)

with $\phi \left(x\right)$ a flow-limiting function such as,

$$\phi \left(x\right)=max\left(0,min\left(2x,1\right)\right),$$

(40)

and $C$ a variable defined as,

$$C=\left\{\begin{array}{cc}CFL·\left(1-CFL\right)& CFL\le 0.5\\ 0.25& CFL>0.5\end{array}\right..$$

(41)

Vectors $\Delta {\mathit{U}}_{j+1/2}^n$ and $\Delta {\mathit{U}}_{j-1/2}^n$ are defined as,

$$\Delta {\mathit{U}}_{j+1/2}^n={\mathit{U}}_{j+1}^n-{\mathit{U}}_j^n,$$

(42)

$$\Delta {\mathit{U}}_{j-1/2}^n={\mathit{U}}_j^n-{\mathit{U}}_{j-1}^n,$$

(43)

and the values of $r_j^{+}$ and $r_j^{-}$,

$$\begin{array}{cc}{r}_j^{+}=\frac{\Delta {\mathit{U}}_{j-1/2}^n·\Delta {\mathit{U}}_{j+1/2}^n}{\Delta {\mathit{U}}_{j+1/2}^n·\Delta {\mathit{U}}_{j+1/2}^n},& r_j^{-}=\frac{\Delta {\mathit{U}}_{j-1/2}^n·\Delta {\mathit{U}}_{j+1/2}^n}{\Delta {\mathit{U}}_{j-1/2}^n·\Delta {\mathit{U}}_{j-1/2}^n}\end{array}.$$

(44)

In this work, this numerical scheme was also extended to the solid phase.

4. Results and Discussion

This section is structured in two parts. Firstly, the validation of the presented model is faced and discussed. Then, the capabilities of different numerical schemes to deal with the physical problem under study are analyzed.

4.1. Model Validation

In the case of detonation of granulated solid propellants, the experimental data available in the open literature are scarce and some key parameters are not often provided. In this study, the validation of the model was carried out in two stages. Firstly, a theoretical two-phase shock tube problem was faced in order to ensure that the model considered was able to characterize the discontinuities linked to shock compressions when friction, interphase heat transfer, and combustion terms were considered. Secondly, an additional validation test was performed by comparing the prediction of the present approach with the data provided by Hoffman and Krier [11] for cyclotetramethylenetetranitramine (HMX).

In the first case, the problem considered a straight tube of 2 m length and 0.1 m height in which there were two distinct zones (left and right) separated by a membrane. At the left side, the initial pressure of the gas was ${10}^6$ Pa whereas at the right side the initial gas pressure was ${10}^5$ Pa. Regarding the rest of the variables, the initial conditions on both sides of the tube were the same: ${\alpha }_0=0.99$, $u_{g0}=u_{p0}=0 m/s, T_{g0}=T_{p0}=298 K$. The ignition temperature was considered to be 300 K, whereas the chemical energy provided for the ignition was $E_{chem}=5.67 \mathrm{x} {10}^6$ J/kg. Particle density and particle diameter considered in this problem were, respectively, 1910 kg/m³ and 0.2 x ${10}^{-3}$ m. Finally, $c_{vp}=1500 \frac{J}{kgK}, c_{vg}=1770 \frac{J}{kgK}$ and ${\gamma }_g=1.333$. The results obtained are gathered in Figure 1 in terms of the spatial distribution of gas density, porosity, velocities and temperatures of the gas and particle phases as well as gas pressure. As shown, the numerical scheme used ensures an efficient resolution of the shock wave and the subsequent rarefaction waves. Besides, the results suggest that the numerical approach proposed works well with the constitutive laws implemented.

The second validation problem faced considered the data provided in [11] and permitted to discuss different physical approaches used in the modeling of the initial stages of the detonation of energetic materials. To perform this analysis, a closed and insulated tube of 0.08 m length was chosen. Therefore, there is no transfer of mass and heat through the walls of the tube with the boundary conditions being zero gradient for all variables despite velocity, which is assumed to be zero at the right and left walls.

The necessary input data for the calculation was obtained from the existing bibliography and it is written in

Table 1. Input data.

Parameter	Value	Units	Parameter	Value	Units
${\alpha }_c$	0.4		$K_0$	$1.38 \mathrm{x} {10}^9$	$\mathrm{Pa}$
${\alpha }_{min}$	0.2595		$E_{chem}$	$5.49 \mathrm{x} {10}^6$	$\mathrm{J}/\mathrm{kg}$
${\rho }_p$	1580.52	$\mathrm{kg}/{\mathrm{m}}^3$	$E_{chem}^p$	$0.22 \mathrm{x} {10}^6$	$\mathrm{J}/\mathrm{kg}$
$c_{vg}$	1776.96	$\mathrm{J}/\mathrm{kgK}$	$\eta$	$1.08$
${\gamma }_g$	1.25		${\mu }_{g,0}$	$4.45$	$\mathrm{Ns}/{\mathrm{m}}^2$
$c_{vp}$	1269.27	$\mathrm{J}/\mathrm{kgK}$	${\alpha }_0$	0.4
$K$	$48.8 \mathrm{x} {10}^6$	$\mathrm{Pa}$	$u_{g0}=u_{p0}$	0	$\left(\mathrm{m}/\mathrm{s}\right)$

. Although all authors consider a small zone of the control volume already initiated, each one applies different values of pressure and temperature as initial conditions of the test. High differences in the results can be obtained depending on the initial values assumed for the variables.

Therefore, and in order to validate the model with literature values, the distributions of pressure, porosity, and gas and particle temperatures provided by Hoffman and Krier [11] at 20 μs, were used as initial conditions for the calculations performed. Simulations were carried out using Rusanov and MacCormack-TVD numerical schemes and their results are compared in the following pages.

Figure 2 shows a comparison of the present model with flame front solution provided by Hoffman and Krier [11] for two cases: the case of the model that includes a modification of the momentum equation (i.e., model formed by the Equations (1)–(4), (6), (7), and (9), Figure 2a) and the case of the model that includes a porosity limiter (Figure 2b). In addition, Figure 3 shows a comparison of the time evolution of porosity distributions with the spatial domain for momentum equation modification. All in all, the results obtained are qualitatively similar to those obtained by Hoffman and Krier [11] with a good agreement between Rusanov and MacCormack numerical schemes. In the case of flame front propagation, the present model with a porosity limiter slightly overestimates this value, as a result, the results of porosity and temperature are displaced in the x-direction with respect the ones provided by Hoffman and Krier [11]. Notwithstanding, this slight overestimation is on the “safe side” as in the case of an accident simulation the code would provide a more adverse scenario than the actual one.

A grid convergence study was performed in order to assess the sensitivity of the mesh size in the results of the numerical model. Figure 4 shows the results of the numerical model of the case of the flame front propagation (similar to that in Figure 2) as a function of number of cells used in the discretization with Rusanov and TVD schemes. As shown, as the number of cells increases, the numerical solution tends to converge towards the reference data of Hoffman and Krier [11].

4.2. Evaluation of the Numerical Strategies Proposed for the Physical Model

This section is structured in two parts depending on how the limitation of the porosity variable is performed. Firstly, Hoffman and Krier [11] system of equations together with the modification of the momentum conservation equation are applied. Afterwards, the calculation is performed limiting the porosity directly in the code in order to prevent this variable reaching values under the minimum accepted.

4.2.1. Modification of Momentum Conservation Equation

The system to be considered is formed by the Equations (1)–(4), (6), (7), and (9) and the solutions are plotted in Figure 5. The results obtained with Rusanov and MacCormack agree well at 40 μs and present slight but not relevant differences at 60 μs for porosity, pressure, and gas temperature. The frame front distributions obtained with both numerical schemes are also very similar to those provided by Hoffman and Krier [11] in their work. However, the modification of the momentum equation provides a solution in which the porosity value is below the minimum physically possible one. Besides, the particle temperature has an unrealistic value. Therefore, a different restriction of some physical variables should be considered.

4.2.2. Porosity Direct Limiter

The minimum value for porosity should be that for which particles are completely packed and thus, they cannot be further compacted. As seen in Section 4.2.1., the modification of the momentum equation provides a solution in which not only the porosity reaches values under the minimum physical value, but also the particle temperature has an unrealistic value. In order to avoid these drawbacks, a different restriction was tested. In this case, the porosity was limited during the calculation directly in the code by preventing the porosity reaching values under a minimum threshold, which was fixed based on the minimum physical value that the solid particles could occupy in the volume of the domain.

The results obtained for porosity, pressure, and gas and particle temperatures are presented in Figure 6. As already discussed in the previous section, the results obtained are qualitatively similar to those obtained by Hoffman and Krier [11] with a good agreement between Rusanov and MacCormack numerical schemes. Besides, the distribution of particle temperature takes physical values. It is worth noting that the application of the MacCormack scheme to the solid phase provided good results. These show that MacCormack-TVD is an appropriate numerical scheme to be used when integrating the solid phase of this kind of problem.

However, the x-locations of the peak of pressure and temperatures at 60 μs show a displacement compared with those presented by Hoffman and Krier [11]. This phenomenon can be seen in the flame front propagation results represented in Figure 2b. The figure shows a good agreement of both numerical schemes with the result of the bibliography [11] up to 50 μs, increasing the difference between the calculation and the results from literature from this time on. These differences could be due, on the one hand, to the initial values of the parameters chosen as initial conditions, which highly determine the x-location of the peak values for all variables. On the other hand, these could rely on the inconsistency of some input parameters obtained from different sources in the literature.

5. Conclusions

A model was developed to study the early stages of the detonation process of granular energetic materials. To solve the system of differential equations for both gas and solid phases, Rusanov and MacCormack-TVD numerical schemes were used. Besides, the first-order Euler method was used to solve the source terms. Two different models were used to calculate pressure, temperature, and porosity distributions. The first model studied considered the modification of the particle momentum governing equation done by Hoffman and Krier [11]. According to the authors, this modification prevents the porosity from reaching values below the minimum value of compaction that packed beds of spherical particles can reach. However, when this model is applied, the temperature of the particle phase increases and reaches values up to 8000 K and the porosity takes values below the minimum. In order to control the particle temperature, the limitation of the porosity is performed directly in the code by preventing the porosity reaching values below the possible minimum. The magnitude of the values obtained for the main variables of interest when applying this model is similar to those found in the literature. Moreover, the results obtained using Rusanov scheme agree well with those resulting from applying MacCormack-TVD numerical scheme. However, in some cases, the distributions of the variables are displaced in the x-direction with respect to those from the literature. This behavior can be seen by observing the x-coordinate position of the peak values obtained for pressure, temperature, and porosity variables. These differences could be due to the initial values of the parameters chosen as initial conditions, which highly determine the x-location of the peak values for all variables or, on the other hand, to the inconsistency of input data spread in many different sources from the literature.

All in all, the results show that MacCormack-TVD is an appropriate numerical scheme to be used when integrating the solid phase of this kind of physical problem. The extension of the numerical scheme will permit to apply this strategy to different solid propellant problems in future work. Finally, it can be concluded that the last model considered represents, with reasonable accuracy, the physical behavior of the propellant combustion for all variables of interest becoming a predictive tool for the characterization of the early stages of the detonation process of granular solid propellants.

Finally, it is worth commenting on the limitations found in the physical model considered in this work. These limitations can be improved in future work. On one hand, the model adopts the Equation of State (EOS) of Nobel-Abel for the gas phase. It would be advisable to use more complex real gas EOS to improve model accuracy. Moreover, the use of the Jones–Wilkins–Lee (JWL) EOS for both the solid and products could be also considered. However, there is an important limitation to obtain the coefficients for this kind of EOS for the reaction products of HMX and other explosives and/or propellants. On the other hand, turbulence is not considered in the present work. This is an improvement the authors plan to undertake in future work as turbulence can enhance the transport processes between phases.