# Finite Volume Method for the Launch Safety of Energetic Materials

Authors: Yanchun Yu; Han Yan; Weidong Chen; Jingxin Ma
Journal: Shock and Vibration
 (2021)

---

## Abstract

It is important to examine the ignition of energetic materials for launch safety. Given that there is a paucity of experimental tests, numerical simulations are important for analysing energetic materials. A computer program based on the finite volume method and viscoelastic statistical crack mechanics model is developed to study the ignition of energetic materials. The trends of temperature and stress of energetic materials subjected to projectile base pressure are studied by numerical examples. The results are compared with those in an extant study, which verified the correctness of the proposed method. Additionally, the relationships between the temperature increase and nonimpact ignition of energetic materials were analysed. The results show that when the temperature at the bottom of the explosive rises to a certain value, it will cause the explosive to ignite. This research has significance to the study of the base gap of explosives, and it provides a reference for launch safety evaluation of energetic materials.


---

## Body

## 1. Introduction

Launch safety has important research value [1–6]. Zhang Linjun studied the effect of thermal aging on launch safety of RDX-based aluminized and pressed explosive charge [7]. Xiao Wei studied the launch safety of RDX-based aluminized explosive [8]. Teng Wanting studied launch safety of thermobaric explosive [9]. Domestic and foreign scholars agree that the base gap of the explosive is one of the important factors of explosion in advance. Myers Thomas F studied the effect of base gap on setback-shock sensitivities of cast composition B and TNT with the NSWC setback simulator [10]. Starkenberg et al. studied the problem of compressive heating ignition; several explosives have been tested using the artillery setback simulator [11]. Rui et al. revealed the mechanism of explosion in the chamber of the propellant charge and proposed a method for evaluating the safety of the propellant charge launch [12]. Yan et al. studied the effect of bottom gap of the explosive on launch safety [13]. When the explosive is launched, the bottom of the charge is subjected to a high-pressure load; the residual gas at the bottom of the explosive is in adiabatic compression, which causes the temperature at the bottom of the explosive to rise and it can easily cause an explosion [14]. Therefore, the effect of temperature rise on ignition of the explosive has important research significance for explosion in the chamber and launch safety.The constitutive model is important for accurately describing the performances of the explosive. The damage constitutive model of energetic materials has been researched by many scholars. Dienes et al. of Los Alamos National Laboratory in the United States conducted a large number of dynamic damage mechanical property experiments and XDT experiments for polymer-bonded explosives (PBXs) [15]. Based on the mechanism of “hot spots,” Dienes established and gradually improved the statistical mesodamage mechanical model. The statistical crack mechanics (SCRAM) model is a representative mesodamage constitutive model for examining the damage of energetic materials. Addessio and Johnson established the Iso-SCRAM model, which is an isotropic statistical mesodamage constitutive model, based on the SCRAM model [16]. Bonnett and Butler established the Visco-SCRAM model (the viscoelastic statistical crack constitutive model) based on the SCRAM and Iso-SCRAM model [17]. The novelty of this model is that the generalized viscoelastic body is coupled with the microcrack body to describe the viscoelastic mechanical response of PBX explosives, and the microcrack body is used to describe the evolution of damage. The viscoelastic statistical crack constitutive model is suitable for analysing dynamic damage of PBX explosives subjected to nearly isotropic stress, such as collision and impact [18]. Therefore, the viscoelastic statistical crack model is used to investigate the ignition of energetic materials in this study.Currently, the Visco-SCRAM model is usually embedded into finite element software [19–22]. Therefore, a computer program based on the finite volume method and viscoelastic statistical crack mechanics model is developed for the launch safety of the explosive in this study.The finite volume method (FVM) is widely used in computational fluid dynamics [23–28]. In recent years, the finite volume method has been increasingly applied in the field of solid mechanics [29, 30]. FVM is established on the basis of conservation of physical quantities, which is one of the advantages of this method. Compared with the finite element method (FEM), the finite volume method has specific physical significance; it is easy to solve and program, and it is suitable for solving complex boundary problems. Moreover, explicit integral algorithm of time is used in FVM, which is more suitable for solving highly nonlinear problems such as explosions, impacts, and collisions. Therefore, the finite volume method is used to study the launch safety of explosive in this study.


## 2. Finite Volume Method

### 2.1. Discretization Equation of the Finite Volume Method

The finite volume method is established on the basis of conservation equations. The continuum should satisfy three conservation equations:(1)
Mass conservation equation:(1)DρDt+ρ∂vi∂xi=0,whereρ denotes the density and vi (i = 1, 2, 3) denotes the velocity components of the particle. In this study, the finite volume method automatically satisfies the mass conservation equation.(2)
Momentum conservation equation:(2)ρDviDt=ρbi+∂σji∂xj,wherebi denotes the force per unit mass and σij (i = 1, 2, 3) denotes the components of the Cauchy stress tensor.(3)
Energy conservation equation:(3)ρDωintDt=Dijσij=ε˙ijσij,whereωint denotes the internal energy per unit mass; Dij denotes the components of the deformation rate tensor; and ε˙ij denotes the strain rate.The analysis of the explosive based on the finite volume method is a coupling analysis process of mass conservation equation, momentum conservation equation, and energy conservation equation.The unstructured grid is used in the finite volume method. The unstructured grid is characterized by irregular grids, such as triangles, polygons, and tetrahedrons, which are difficult to create, but they are highly adaptable to the complex solution area and highly nonlinear problems.A two-dimensional unstructured grid is used in this study, and the grid division is shown in Figure1. The control volume V is surrounded by 3-4-5-6-7-8-9-10-11-12, and k is the centre of the control volume.Figure 1
Grid division.According to Newton's second law, the inertial force∫Vρu¨dV is the sum of the body force ∫VfbdV and surface force Fs. The mathematical expression is as follows:(4)∫Vρu¨dV=∫VfbdV+FS.The body force is ignored in the following equations. The component of forceFs on the axis is Fi (i = 1, 2, 3), and it can be represented by the stress tensor as follows:(5)Fi=∫SσiknkdS,where nk (k = 1, 2, 3) denotes the component of outer normal vector n of the surface S. For two-dimensional problems, equations (4) and (5) are as follows:(6)∫Sρu¨idS=Fi=∫Lσijnjdl=∫Lσi1dy−σi2dxi=1,2,where n = nj (j = 1, 2) denotes the outer normal vector of curve L and curve L is the boundary of the control volume. For the two-dimensional problems, equation (6) is as follows:(7)∫Sρu¨xdxdy=∫Lσxxdy−σxydx,(8)∫Sρu¨ydxdy=∫Lσxydy−σyydx,where u¨x and u¨y denote the components of acceleration in each direction and σ denotes the density of the control volume.To solve the integral on the right side of (7) and (8), the stress in the grid is assumed to be constant. The integrals on the left side of (7) and (8) are solved by centralized quality method. Based on the assumption, the final discrete equations are shown in (9) and (10). The detailed discrete process can be found in [31].(9)Mmu¨xm=∑n=1nm3σxxnamn−∑n=1nm3σxynbmn+∑n=1nm4σxxncmon−∑n=1nm4σxyndmon+∫a⟶m⟶bX¯dl,(10)Mmu¨ym=∑n=1nm3σxynamn−∑n=1nm3σyynbmn+∑n=1nm4σxyncmon−∑n=1nm4σyyndmon+∫a⟶m⟶bY¯dl,where Mm denotes the mass of control volume centred on node m and u¨xm and u¨ym denote the acceleration components of the node m in x and y directions, respectively. Furthermore, nm3 denotes the total number of triangular elements associated with node m, and nm4 denotes the total number of the quadrangle elements associated with node m. Additionally, X¯ and Y¯ denote the linear density of the external force on the surface in x and y directions, respectively; amn, bmn, cmon, and dmon are coefficients, the subscript m or “mo” denotes the control volume centred on node m or “mo,” and the subscript n denotes the nth grid in the control volume centred on node m or “mo.”


### 2.2. Comparison and Verification of the Finite Volume Method and Finite Element Method

To verify the accuracy of the FVM proposed in this study, the results of FVM are compared with the results of FEM by an example.Taking a thin steel plate as an example, the model of the thin steel plate is 1 m × 0.02 m and the steel plate has a fixed bearing at both ends, as shown in Figure2. The upper surface of the steel plate is subjected to a uniform load; it is a simplified triangular impact load, as follows:(11)Pt=−11.55t+10.63MPa/mt≤0.92ms0t>0.92ms.Figure 2
Model of the thin steel plate.The constitutive model of the steel plate is an elastoplastic linear hardening model, as follows:(12)σ=Eεε≤εsσ=σs+Etε−εsε>εs,where σ is the stress, σs is the yield stress, Et is the Tangent modulus, ε is the strain, and εs is the yield strain.The density is 7830 kg/m3, elasticity modulus is 2.07 × 1011 Pa, tangent modulus is 5.0 × 109 Pa, yield stress is 4.0 × 108 Pa, and Poisson's ratio is 0.3. The values of the material parameters are obtained from the material library of AUTODYN. The transient dynamic response of the steel plate is calculated.The finite element method has been widely accepted for solving the problems of transient dynamic response, and the calculation accuracy is high. Therefore, it is feasible to compare the results of the FVM with that of FEM. LS-DYNA is used to calculate the transient dynamic response of the steel plate. The results of the two methods are compared, and the time-history curves of the displacement and stress of the model at points A, B, and C are shown in Figures3 and 4.Figure 3
Comparison of the displacement between FEM and FVM. (a) Time history curve of the displacement at point A. (b) Time history curve of the displacement at point B. (c) Time history curve of the displacement at point C.
(a)(b)(c)Figure 4
Comparison of the stress between FEM and FVM. (a) Time history curve of the stress at point B. (b) Time history curve of the stress at point C. (c) Time history curve of the stress at point D.
(a)(b)(c)It is seen that the results of the FVM are almost consistent with those of the FEM. The correctness of FVM proposed in this study is verified. Near the extreme value of the stress, the results of FVM show a small numerical oscillation. On the one hand, FVM uses explicit integration algorithm and FEM uses implicit integration algorithm. Due to the difference between explicit algorithm and implicit algorithm, the total number of time steps of explicit algorithm is much more than that of implicit algorithm, so the results of FVM show more details than those of FEM. On the other hand, the stress propagates in the structure as waves, and when the structure is subjected to similar impact loads, the physical quantities such as stress will change suddenly, so the physical quantities on both sides of the wave front are discontinuous and the solution of differential equation is unstable.To verify the calculation accuracy of the FVM, the relative error of the displacement and stress values calculated by the two methods is analysed. The results are shown in Tables1 and 2.Table 1
Relative error of displacement at point B.
Time (s)(FEM) displacement in theY direction (m)(FVM) displacement in theY direction (m)Relative errors of displacement (%)0.00300−0.03009−0.030000.2990.00301−0.02995−0.029870.2670.00302−0.02982−0.029730.3020.00303−0.02968−0.029600.2700.00304−0.02955−0.029460.3050.00305−0.02941−0.029330.2720.00306−0.02927−0.029190.2730.00307−0.02914−0.029050.3090.00308−0.02900−0.028920.2760.00309−0.02886−0.028780.277Table 2
Relative error of stress at point B.
Time (s)(FEM) equivalent stress (MPa)(FVM) equivalent stress (MPa)Relative errors of stress (%)0.00097176.87179.841.650.00098162.53164.231.040.00099152.06152.940.580.00101141.77141.960.000.00102127.26127.280.130.00103110.48110.570.080.0010492.8392.920.100.0010577.0977.060.040.0010663.9763.870.160.0010754.5654.600.07The results show that the relative error of the two methods is small, which meets the engineering requirements; thus, the calculation accuracy of the FVM is further verified.


## 2.1. Discretization Equation of the Finite Volume Method

The finite volume method is established on the basis of conservation equations. The continuum should satisfy three conservation equations:(1)
Mass conservation equation:(1)DρDt+ρ∂vi∂xi=0,whereρ denotes the density and vi (i = 1, 2, 3) denotes the velocity components of the particle. In this study, the finite volume method automatically satisfies the mass conservation equation.(2)
Momentum conservation equation:(2)ρDviDt=ρbi+∂σji∂xj,wherebi denotes the force per unit mass and σij (i = 1, 2, 3) denotes the components of the Cauchy stress tensor.(3)
Energy conservation equation:(3)ρDωintDt=Dijσij=ε˙ijσij,whereωint denotes the internal energy per unit mass; Dij denotes the components of the deformation rate tensor; and ε˙ij denotes the strain rate.The analysis of the explosive based on the finite volume method is a coupling analysis process of mass conservation equation, momentum conservation equation, and energy conservation equation.The unstructured grid is used in the finite volume method. The unstructured grid is characterized by irregular grids, such as triangles, polygons, and tetrahedrons, which are difficult to create, but they are highly adaptable to the complex solution area and highly nonlinear problems.A two-dimensional unstructured grid is used in this study, and the grid division is shown in Figure1. The control volume V is surrounded by 3-4-5-6-7-8-9-10-11-12, and k is the centre of the control volume.Figure 1
Grid division.According to Newton's second law, the inertial force∫Vρu¨dV is the sum of the body force ∫VfbdV and surface force Fs. The mathematical expression is as follows:(4)∫Vρu¨dV=∫VfbdV+FS.The body force is ignored in the following equations. The component of forceFs on the axis is Fi (i = 1, 2, 3), and it can be represented by the stress tensor as follows:(5)Fi=∫SσiknkdS,where nk (k = 1, 2, 3) denotes the component of outer normal vector n of the surface S. For two-dimensional problems, equations (4) and (5) are as follows:(6)∫Sρu¨idS=Fi=∫Lσijnjdl=∫Lσi1dy−σi2dxi=1,2,where n = nj (j = 1, 2) denotes the outer normal vector of curve L and curve L is the boundary of the control volume. For the two-dimensional problems, equation (6) is as follows:(7)∫Sρu¨xdxdy=∫Lσxxdy−σxydx,(8)∫Sρu¨ydxdy=∫Lσxydy−σyydx,where u¨x and u¨y denote the components of acceleration in each direction and σ denotes the density of the control volume.To solve the integral on the right side of (7) and (8), the stress in the grid is assumed to be constant. The integrals on the left side of (7) and (8) are solved by centralized quality method. Based on the assumption, the final discrete equations are shown in (9) and (10). The detailed discrete process can be found in [31].(9)Mmu¨xm=∑n=1nm3σxxnamn−∑n=1nm3σxynbmn+∑n=1nm4σxxncmon−∑n=1nm4σxyndmon+∫a⟶m⟶bX¯dl,(10)Mmu¨ym=∑n=1nm3σxynamn−∑n=1nm3σyynbmn+∑n=1nm4σxyncmon−∑n=1nm4σyyndmon+∫a⟶m⟶bY¯dl,where Mm denotes the mass of control volume centred on node m and u¨xm and u¨ym denote the acceleration components of the node m in x and y directions, respectively. Furthermore, nm3 denotes the total number of triangular elements associated with node m, and nm4 denotes the total number of the quadrangle elements associated with node m. Additionally, X¯ and Y¯ denote the linear density of the external force on the surface in x and y directions, respectively; amn, bmn, cmon, and dmon are coefficients, the subscript m or “mo” denotes the control volume centred on node m or “mo,” and the subscript n denotes the nth grid in the control volume centred on node m or “mo.”


## 2.2. Comparison and Verification of the Finite Volume Method and Finite Element Method

To verify the accuracy of the FVM proposed in this study, the results of FVM are compared with the results of FEM by an example.Taking a thin steel plate as an example, the model of the thin steel plate is 1 m × 0.02 m and the steel plate has a fixed bearing at both ends, as shown in Figure2. The upper surface of the steel plate is subjected to a uniform load; it is a simplified triangular impact load, as follows:(11)Pt=−11.55t+10.63MPa/mt≤0.92ms0t>0.92ms.Figure 2
Model of the thin steel plate.The constitutive model of the steel plate is an elastoplastic linear hardening model, as follows:(12)σ=Eεε≤εsσ=σs+Etε−εsε>εs,where σ is the stress, σs is the yield stress, Et is the Tangent modulus, ε is the strain, and εs is the yield strain.The density is 7830 kg/m3, elasticity modulus is 2.07 × 1011 Pa, tangent modulus is 5.0 × 109 Pa, yield stress is 4.0 × 108 Pa, and Poisson's ratio is 0.3. The values of the material parameters are obtained from the material library of AUTODYN. The transient dynamic response of the steel plate is calculated.The finite element method has been widely accepted for solving the problems of transient dynamic response, and the calculation accuracy is high. Therefore, it is feasible to compare the results of the FVM with that of FEM. LS-DYNA is used to calculate the transient dynamic response of the steel plate. The results of the two methods are compared, and the time-history curves of the displacement and stress of the model at points A, B, and C are shown in Figures3 and 4.Figure 3
Comparison of the displacement between FEM and FVM. (a) Time history curve of the displacement at point A. (b) Time history curve of the displacement at point B. (c) Time history curve of the displacement at point C.
(a)(b)(c)Figure 4
Comparison of the stress between FEM and FVM. (a) Time history curve of the stress at point B. (b) Time history curve of the stress at point C. (c) Time history curve of the stress at point D.
(a)(b)(c)It is seen that the results of the FVM are almost consistent with those of the FEM. The correctness of FVM proposed in this study is verified. Near the extreme value of the stress, the results of FVM show a small numerical oscillation. On the one hand, FVM uses explicit integration algorithm and FEM uses implicit integration algorithm. Due to the difference between explicit algorithm and implicit algorithm, the total number of time steps of explicit algorithm is much more than that of implicit algorithm, so the results of FVM show more details than those of FEM. On the other hand, the stress propagates in the structure as waves, and when the structure is subjected to similar impact loads, the physical quantities such as stress will change suddenly, so the physical quantities on both sides of the wave front are discontinuous and the solution of differential equation is unstable.To verify the calculation accuracy of the FVM, the relative error of the displacement and stress values calculated by the two methods is analysed. The results are shown in Tables1 and 2.Table 1
Relative error of displacement at point B.
Time (s)(FEM) displacement in theY direction (m)(FVM) displacement in theY direction (m)Relative errors of displacement (%)0.00300−0.03009−0.030000.2990.00301−0.02995−0.029870.2670.00302−0.02982−0.029730.3020.00303−0.02968−0.029600.2700.00304−0.02955−0.029460.3050.00305−0.02941−0.029330.2720.00306−0.02927−0.029190.2730.00307−0.02914−0.029050.3090.00308−0.02900−0.028920.2760.00309−0.02886−0.028780.277Table 2
Relative error of stress at point B.
Time (s)(FEM) equivalent stress (MPa)(FVM) equivalent stress (MPa)Relative errors of stress (%)0.00097176.87179.841.650.00098162.53164.231.040.00099152.06152.940.580.00101141.77141.960.000.00102127.26127.280.130.00103110.48110.570.080.0010492.8392.920.100.0010577.0977.060.040.0010663.9763.870.160.0010754.5654.600.07The results show that the relative error of the two methods is small, which meets the engineering requirements; thus, the calculation accuracy of the FVM is further verified.


## 3. Viscoelastic Statistical Crack Model

The viscoelastic statistical crack model is suitable for the analysis of the dynamic damage of PBX explosives subjected to nearly isotropic stress such as collision and impact. This model is used to study the ignition of energetic materials in this study. The specific theory of the model is shown in extant studies [32, 33]. Only a brief equation of the model is provided here.The relation between the deviator stress rateS˙ij and deviator strain rate e˙ij for the viscoelastic statistical crack model is as follows:(13)S˙ij=2Ge˙ij−∑n=1NSijn/τn−3c/a2c˙/aSij1+c/a3,where Sij denotes the deviator stress, G denotes the shear modulus of the elastic element, Sijn denotes the deviator stress of the nth viscoelastic body, τn denotes the relaxation time of the nth viscoelastic body, c denotes the mean crack radius, c˙ denotes crack growth rate, and a denotes an initial parameter. The deviator stress rate of each viscoelastic element is as follows:(14)S˙ijn=2Gne˙ij−Sijnτn−GnG3ca2c˙aSij+ca3S˙ij,where Gn denotes the shear modulus of the nth elastic body.To describe the evolution and growth of cracks based on a study by Dienes, it is assumed that the crack growth rate is related to the stress intensity factor. The crack growth can be expressed as follows:(15)c˙=νmaxKK1mK<K′c˙=νmax1−K0K2K≥K′,where(16)K′=K01+2m,K1=K01+2m1+m21/m,where c˙ denotes crack growth rate, m denotes the factor of crack speed growth, vmax denotes the terminal crack speed, and K0 denotes the fracture toughness. Furthermore, K denotes the stress intensity factor, and K1 is a constant.


## 4. Macroscopic Volume Heating Model of a Viscoelastic Statistical Cracked Body

The macroscopic volume heating model includes chemical thermal decomposition term and mechanical term which describe viscosity, crack damage, and the changes of adiabatic volume. Chemical thermal decomposition is the model of Arrhenius first-order reaction kinetics. For the overall heat change in the nonimpact ignition process, the overall heat conduction term in the formula can be ignored in this study. Therefore, time-dependent temperature change rate equation can be simplified as follows:(17)T˙=−γTε˙ii+IρCVW˙ve+W˙cr+ΔHCVZexp−EART,where T denotes the temperature, γ denotes the Gruneisen coefficient, ε˙ii denotes the strain rate, Idenotes the conversion constant of inelastic power into thermal power, CV denotes the constant volume heat capacity, W˙ve denotes the viscous power, W˙cr denotes the crack damage power, ΔH,EA,andZ, respectively, denote the decomposition heat per unit mass, activation energy, and preexponential factor for the thermal decomposition reaction of energetic materials, and R denotes the universal gas constant. The W˙ve and W˙cr are as follows:(18)W˙ve=∑n=1NSijne˙ijve=∑n=1NSijnSijn2Gntn,W˙cr=Sije˙ijcr=12G3ca2c˙aSijSij+ca3S˙ijSij.The temperature change of viscoelastic statistical cracked materials in an element is described in equation (17). All the variables in the equation, with the exception of temperature, can be solved by the viscoelastic statistical crack model. Therefore, this equation can be regarded as an ordinary differential equation with temperature as a variable, and the temperature of an element can be obtained by integrating time.


## 5. Equation of State

Shock equation of state or JWL equation of state is often used in the study of explosives. The solid JWL equation of state is used in this study as follows:(19)p=A1−wR1Ve−RL1V+B1−wR2Ve−RL2V+wEV,where p denotes the pressure, E denotes the internal energy per unit initial volume, V denotes the relative volume, and A, B, R1, R2, and w denote the material parameters.


## 6. Dynamic Response Analysis of Explosives Based on the Finite Volume Method

The numerical model of explosive is established by FVM. It is based on the conservation equation of mass, momentum, and energy. The momentum conservation equation is discretized by the finite volume method, and then, the state variables are related by three conservation equations, constitutive equation, and equation of state. The detailed process is described as follows:(1)
Mesh and calculate the initial variables of the grids.(2)
According to the initial conditions and boundary conditions, calculate the initial velocity and displacement of the grid node based on the discrete equation of FVM.(3)
Update velocity and displacement of the grid node, and calculate the strain rate of the grid element at the next time step.(4)
Calculate the deviator stress of the grid element at the next time step based on the viscoelastic statistical crack model. Calculate the internal energy and pressure of the grid element at the next time step according to the coupling energy conservation equation and equation of state. Calculate the stress of the grid element at the next time by Cauchy stress. Calculate the temperature of the grid element at the next time step by the macroscopic volume heating model.(5)
Bring the updated stress into the discrete equation to calculate the force and acceleration of the grid node at the next time step, and calculate the velocity and displacement of the grid node at the next time step. Then, go to the next cycle.The calculation flow chart is shown in Figure5.Figure 5
Calculation flow chart.


## 7. Numerical Calculation 1

Based on FVM, the viscoelastic statistical crack constitutive model is used to calculate the stress and temperature of energetic materials subjected to projectile base pressure. The relationship between temperature rise and ignition of explosive is analysed.The calculation model of energetic materials is an axisymmetric model, where theX-axis is the symmetry axis, as shown in Figure 6. The total length of the axis of symmetry is 56.7 cm, and the maximum radius of the model is 6.3 cm. The load is applied at the bottom of the explosive (the far-right side of the model, x = 567 mm), the front of the explosive (the far-left side of the model) is free boundary, and the X-axis is a symmetric boundary condition (UY = 0). The explosive model is divided into 7849 triangular elements and 4133 nodes, as shown in Figure 7.Figure 6
Calculation model.Figure 7
Division of elements.PBX9501 explosive is used in this study. The relevant parameters of PBX9501 explosive are provided in extant studies [32, 33]. The viscoelastic parameters, statistical crack parameters, relaxation parameters, and thermodynamic parameters of PBX9501 explosive are provided in Tables 3–7 [32, 33].Table 3
Viscoelasticity parameters.
G1/PaG2/PaG3/PaG4/PaG5/Pa9.440×1081.738×1085.212×1089.085×1086.875×108G1 denotes the shear modulus of the 1st elastic element, G2 denotes the shear modulus of the 2nd elastic element, and so on.Table 4
Crack parameters of the explosive.
vmC0/ma/mvmax/m⋅s−1k0/Pa⋅m1/20.3100.000030.0013005 × 105v denotes Poisson's ratio, m denotes the factor of crack speed growth, c0 denotes the mean initial crack radius, a denotes the initial model parameter, vmax denotes the terminal crack speed, and K0 denotes the fracture toughness.Table 5
Thermodynamic parameters of the explosive.
ki/W/m⋅Kcv/J⋅Kg−1⋅K−1ρ/kg/m3ΔH/J/kgZ/m⋅s−1EA/R(K)0.5120018105.5 × 1065 × 1092.652 × 104Cv denotes the constant volume heat capacity; ρ denotes the density; ki denotes the thermal conductivity; ΔH denotes the decomposition heat per unit mass for chemicals; Z denotes the preexponential factor; EA denotes the activation energy; and R denotes the universal gas constant.Table 6
Relaxation parameters.
1/r1/s1/r2/s1/r3/s1/r4/s1/r5/s07.32×1037.32×1047.32×1057.32×106τ denotes the relaxation time. The parameters of the equation of state are shown in Table 7 [34].Table 7
Material parameters.
A/M barB/M barR1R2w7.781 × 102−5.031 × 10−211.31.130.8938A, B, R1, R2, and w denote the material parameters.The projectile base pressure subjected to the bottom of the model is shown in Figure8.Figure 8
Curve of the projectile base pressure.


### 7.1. Results and Discussion

The programs are developed to calculate the stresses and temperatures of the explosive subjected to projectile base pressure by the finite volume method; the initial temperature of the explosive is set to 300K, and the results are shown in Figures9 and 10.Figure 9
Time-history curve of temperature and stress. (a) A temperature curve at the bottom of the explosive. (b) A stress curve at the bottom of the explosive. (c) A temperature curve at the middle of the explosive. (d) A stress curve at the middle of the explosive. (e) A temperature curve at the top of the explosive. (f) A stress curve at the top of the explosive.
(a)(b)(c)(d)(e)(f)Figure 10
Stress and temperature nephogram at different times. (a) Temperature nephogram att=12 ms. (b) Stress nephogram at t = 12 ms. (c) Temperature nephogram at t = 18 ms. (d) Stress nephogram at t = 18 ms.
(a)(b)(c)(d)Figures9(b), 9(d), and 9(f) show that the change trend of stress is generally consistent with that of the projectile base pressure and the magnitude of the stress is the same as the magnitude of the projectile base pressure. Figure 10 shows that the maximum stress and temperature appear at the bottom of the explosive model and the stress and temperature decrease along the X axis of the explosive.


### 7.2. Verification of Calculation Results

In order to verify the correctness of the results, the results are compared with those from the existing research. Wang Hao calculated the change of temperature and stress of explosive charge subjected to the projectile base pressure [35]. The calculation model of explosive charge is shown in Figure 2. In [35], the model is composed of three parts: fuse, explosive charge, and projectile body. The total length of the projectile is 72.65 cm, the caliber is 155 mm, the total length of explosive charge is 57.7 cm, and the maximum radius is 6.1 cm. It constitutes press-packed aluminum-containing explosive, and the density is 1620 kg/m3. The shape of the explosive in [35] is similar to the calculation model in Figure 7 in this article. The projectile base pressure is applied to the bottom of the explosive (at the far-right side of the model), as shown in Figure 2 in [35], and the curve of the pressure is similar to the pressure in Figure 8 in this study. The constitutive model of explosive is piecewise linear plasticity, and the temperature and stress of explosive charge under the pressure are calculated by DYNA 2D in [35].Wang Hao's conclusion is that the change trend of stress at element A (the bottom of the explosive model) as shown in Figure2 in [35] is consistent with that of the projectile base pressure; the maximum stress appears at element A (the bottom of the explosive model); stress and temperature decrease along the axial direction of the explosive as shown in Figure 2 in this article and Figure 2 in [35]. Therefore, the conclusion in this study is consistent with the conclusion in [35]; it also verifies the correctness of the calculation results in this study.


### 7.3. Analysis of the Relationship between  Temperature Rise and Ignition of the Explosive

The base gap of the explosive charge has always been considered as one of the main factors affecting the bore premature. The trapped air in the base gap is compressed by high overload in chamber during launch, which causes the explosive at the bottom of the projectile to form high compressive stress and form high temperature hot spots. The high-temperature hot spots will heat the adjacent explosive layer, and the explosive layer adjacent to the air will form high-temperature hot spots. At the same time, the explosive will move axially, producing great friction and impact, which will cause the explosives to decompose, burn, or even detonate [36]. Therefore, it is important to study the influence of temperature rise at the bottom of the explosive on explosion in the chamber.To further study the effect of temperature rise on ignition of the explosive, the initial temperature at the bottom of the explosive is set to 450 K and 500 K, respectively, while the temperature at other parts of the explosive is still 300 K. The temperature and stress of the explosive are calculated. The results are shown in Figures11–14.Figure 11
Initial temperature at the bottom of the explosive is 450 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 12
Initial temperature at the bottom of explosive is 500 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 13
Initial temperature at the bottom of explosive is 450 K. (a) Temperature nephogram att = 18 ms. (b) Pressure nephogram at t = 18 ms.
(a)(b)Figure 14
Initial temperature at the bottom of the explosive is 500 K. (a) Temperature nephogram att = 13.8 ms. (b) Pressure nephogram at t = 12 ms.
(a)(b)The results indicate that when the initial temperature at the bottom of the explosive is 450 K, finally, the temperature at the bottom of the explosive increases to 570 K. Furthermore, when the initial temperature at the bottom of the explosive is 500 K, the result shows that when the computation time is 14.7 ms, the temperature suddenly increases sharply and the value of temperature increases infinitely. Based on the chemical kinetic Arrhenius equation used in this study, the criterion of ignition temperature is dT/dt ⟶∞ [37], which indicates that the temperature at the bottom of explosive has reached the ignition temperature and the explosive is ignited. The results show that when the temperature at the bottom of the explosive rises to a certain value, it will cause ignition. In fact, temperature rise can trigger ignition, which has been verified by Dienes and Kershner [33] and Zuo [38] based on the statistical crack model (SCRAM). If there is a base gap in the explosive charge and when the base gap of the explosive charge is compressed, the temperature at the bottom of the explosive charge will rise; when the temperature rises to a certain value, the explosive will ignite.


## 7.1. Results and Discussion

The programs are developed to calculate the stresses and temperatures of the explosive subjected to projectile base pressure by the finite volume method; the initial temperature of the explosive is set to 300K, and the results are shown in Figures9 and 10.Figure 9
Time-history curve of temperature and stress. (a) A temperature curve at the bottom of the explosive. (b) A stress curve at the bottom of the explosive. (c) A temperature curve at the middle of the explosive. (d) A stress curve at the middle of the explosive. (e) A temperature curve at the top of the explosive. (f) A stress curve at the top of the explosive.
(a)(b)(c)(d)(e)(f)Figure 10
Stress and temperature nephogram at different times. (a) Temperature nephogram att=12 ms. (b) Stress nephogram at t = 12 ms. (c) Temperature nephogram at t = 18 ms. (d) Stress nephogram at t = 18 ms.
(a)(b)(c)(d)Figures9(b), 9(d), and 9(f) show that the change trend of stress is generally consistent with that of the projectile base pressure and the magnitude of the stress is the same as the magnitude of the projectile base pressure. Figure 10 shows that the maximum stress and temperature appear at the bottom of the explosive model and the stress and temperature decrease along the X axis of the explosive.


## 7.2. Verification of Calculation Results

In order to verify the correctness of the results, the results are compared with those from the existing research. Wang Hao calculated the change of temperature and stress of explosive charge subjected to the projectile base pressure [35]. The calculation model of explosive charge is shown in Figure 2. In [35], the model is composed of three parts: fuse, explosive charge, and projectile body. The total length of the projectile is 72.65 cm, the caliber is 155 mm, the total length of explosive charge is 57.7 cm, and the maximum radius is 6.1 cm. It constitutes press-packed aluminum-containing explosive, and the density is 1620 kg/m3. The shape of the explosive in [35] is similar to the calculation model in Figure 7 in this article. The projectile base pressure is applied to the bottom of the explosive (at the far-right side of the model), as shown in Figure 2 in [35], and the curve of the pressure is similar to the pressure in Figure 8 in this study. The constitutive model of explosive is piecewise linear plasticity, and the temperature and stress of explosive charge under the pressure are calculated by DYNA 2D in [35].Wang Hao's conclusion is that the change trend of stress at element A (the bottom of the explosive model) as shown in Figure2 in [35] is consistent with that of the projectile base pressure; the maximum stress appears at element A (the bottom of the explosive model); stress and temperature decrease along the axial direction of the explosive as shown in Figure 2 in this article and Figure 2 in [35]. Therefore, the conclusion in this study is consistent with the conclusion in [35]; it also verifies the correctness of the calculation results in this study.


## 7.3. Analysis of the Relationship between  Temperature Rise and Ignition of the Explosive

The base gap of the explosive charge has always been considered as one of the main factors affecting the bore premature. The trapped air in the base gap is compressed by high overload in chamber during launch, which causes the explosive at the bottom of the projectile to form high compressive stress and form high temperature hot spots. The high-temperature hot spots will heat the adjacent explosive layer, and the explosive layer adjacent to the air will form high-temperature hot spots. At the same time, the explosive will move axially, producing great friction and impact, which will cause the explosives to decompose, burn, or even detonate [36]. Therefore, it is important to study the influence of temperature rise at the bottom of the explosive on explosion in the chamber.To further study the effect of temperature rise on ignition of the explosive, the initial temperature at the bottom of the explosive is set to 450 K and 500 K, respectively, while the temperature at other parts of the explosive is still 300 K. The temperature and stress of the explosive are calculated. The results are shown in Figures11–14.Figure 11
Initial temperature at the bottom of the explosive is 450 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 12
Initial temperature at the bottom of explosive is 500 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 13
Initial temperature at the bottom of explosive is 450 K. (a) Temperature nephogram att = 18 ms. (b) Pressure nephogram at t = 18 ms.
(a)(b)Figure 14
Initial temperature at the bottom of the explosive is 500 K. (a) Temperature nephogram att = 13.8 ms. (b) Pressure nephogram at t = 12 ms.
(a)(b)The results indicate that when the initial temperature at the bottom of the explosive is 450 K, finally, the temperature at the bottom of the explosive increases to 570 K. Furthermore, when the initial temperature at the bottom of the explosive is 500 K, the result shows that when the computation time is 14.7 ms, the temperature suddenly increases sharply and the value of temperature increases infinitely. Based on the chemical kinetic Arrhenius equation used in this study, the criterion of ignition temperature is dT/dt ⟶∞ [37], which indicates that the temperature at the bottom of explosive has reached the ignition temperature and the explosive is ignited. The results show that when the temperature at the bottom of the explosive rises to a certain value, it will cause ignition. In fact, temperature rise can trigger ignition, which has been verified by Dienes and Kershner [33] and Zuo [38] based on the statistical crack model (SCRAM). If there is a base gap in the explosive charge and when the base gap of the explosive charge is compressed, the temperature at the bottom of the explosive charge will rise; when the temperature rises to a certain value, the explosive will ignite.


## 8. Numerical Calculation 2

To further analyze the effect of temperature rise on ignition of the explosive, another explosive model is calculated. The model is cylindrical, and it is an axisymmetric model. TheX-axis is the symmetry axis, as shown in Figure 15. The total length of the axis of symmetry is 40 cm, and the maximum radius of the model is 8 cm. The explosive model is divided into 4574 triangular elements and 2408 nodes. The load is subjected to the far-right side of the model (bottom of the explosive, x = 400 mm), the far-left side of the model is a free boundary, and the X-axis is a symmetric boundary condition (UY = 0). PBX9501 explosive is used in this model. The relevant parameters of PBX9501 are provided in Tables 1 to 5. The load is as follows:(20)Ft=400×106×sin2.0×π/20.0×10−3×tPam0≤t≤10ms.Figure 15
Calculation model of the explosive.


### 8.1. Analysis of the Relationship between  Temperature Rise and Ignition of the Explosive

To further study the effect of temperature rise on ignition of this explosive model, the initial temperature at the bottom of the explosive is set to 350 K and 400 K, while the temperature at other parts of the explosive is still 300 K. Furthermore, the temperature and stress of the explosive are calculated. The results are shown in Figures16–19.Figure 16
Initial temperature at the bottom of explosive is 350 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 17
Initial temperature at the bottom of explosive is 400 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 18
Initial temperature at the bottom of the explosive is 350 K. (a) Stress nephogram att = 5 ms. (b) Temperature nephogram at t = 5 ms.
(a)(b)Figure 19
Initial temperature at the bottom of the explosive is 400 K. (a) Stress nephogram att = 10 ms. (b) Temperature nephogram at t = 10 ms.
(a)(b)The results indicated that when the initial temperature at the bottom of the explosive is 350 K, finally, the temperature at the bottom of the explosive increases to 550 K; when the initial temperature at the bottom of the explosive is 400 K, the results show that when the computation time is close to 10 ms, the temperature suddenly increases sharply and the value of the temperature increases infinitely. Based on the chemical kinetic Arrhenius equation used in this study, the criterion of ignition temperature is dT/dt ⟶∞ [37], which indicates that the temperature at the bottom of the explosive has reached the ignition temperature and the explosive has ignited. The results show that when the temperature at the bottom of the explosive rises to a certain value, it will cause ignition.


## 8.1. Analysis of the Relationship between  Temperature Rise and Ignition of the Explosive

To further study the effect of temperature rise on ignition of this explosive model, the initial temperature at the bottom of the explosive is set to 350 K and 400 K, while the temperature at other parts of the explosive is still 300 K. Furthermore, the temperature and stress of the explosive are calculated. The results are shown in Figures16–19.Figure 16
Initial temperature at the bottom of explosive is 350 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 17
Initial temperature at the bottom of explosive is 400 K. (a) Time-history curve of stress at the bottom of the explosive. (b) Time-history curve of temperature at the bottom of the explosive.
(a)(b)Figure 18
Initial temperature at the bottom of the explosive is 350 K. (a) Stress nephogram att = 5 ms. (b) Temperature nephogram at t = 5 ms.
(a)(b)Figure 19
Initial temperature at the bottom of the explosive is 400 K. (a) Stress nephogram att = 10 ms. (b) Temperature nephogram at t = 10 ms.
(a)(b)The results indicated that when the initial temperature at the bottom of the explosive is 350 K, finally, the temperature at the bottom of the explosive increases to 550 K; when the initial temperature at the bottom of the explosive is 400 K, the results show that when the computation time is close to 10 ms, the temperature suddenly increases sharply and the value of the temperature increases infinitely. Based on the chemical kinetic Arrhenius equation used in this study, the criterion of ignition temperature is dT/dt ⟶∞ [37], which indicates that the temperature at the bottom of the explosive has reached the ignition temperature and the explosive has ignited. The results show that when the temperature at the bottom of the explosive rises to a certain value, it will cause ignition.


## 9. Conclusions

In this study, a computer program based on the finite volume method and viscoelastic statistical crack mechanics model is developed. It was used to calculate the stresses and temperatures of the explosive subjected to projectile base pressure, and the results of the numerical examples indicated that the time-history curve of stress is generally consistent with the curve of the load; the maximum stress and temperature appear at the bottom of the explosive model, and the stress and temperature of the explosive along theX-axis decrease from the bottom to the front of the explosive. The results are consistent with those from the existing research; this verifies the correctness of the results in this paper.Additionally, the effects of temperature on the ignition of the explosive are discussed, and the results of the two examples prove that when the temperature at the bottom of the explosive rises to a certain value, it will cause ignition. If there is a base gap in the explosive charge and when the base gap of the explosive charge is compressed, the temperature at the bottom of the explosive charge will rise; when the temperature rises to a certain value, the explosive will ignite. This study has important theoretical and application value for safety evaluation of energetic materials. Besides, the finite volume method is successfully applied to the field of mechanics of explosion in this study.



---

*Source: 9609557-2021-12-15.xml*