Modeling mesoscale energy localization in shocked HMX, Part II: training machine-learned surrogate models for void shape and void–void interaction effects

Abstract

Surrogate models for hotspot ignition and growth rates were presented in Part I (Nassar et al., Shock Waves 29(4):537–558, 2018), where the hotspots were formed by the collapse of single cylindrical voids. Such isolated cylindrical voids are idealizations of the void morphology in real meso-structures. This paper therefore investigates the effect of non-cylindrical void shapes and void–void interactions on hotspot ignition and growth. Surrogate models capturing these effects are constructed using a Bayesian Kriging approach. The training data for machine learning the surrogates are derived from reactive void collapse simulations spanning the parameter space of void aspect ratio, void orientation $(\theta )$, and void fraction $(\phi )$. The resulting surrogate models portray strong dependence of the ignition and growth rates on void aspect ratio and orientation, particularly when they are oriented at acute angles with respect to the imposed shock. The surrogate models for void interaction effects show significant changes in hotspot ignition and growth rates as the void fraction increases. The paper elucidates the physics of hotspot evolution in void fields due to the creation and interaction of multiple hotspots. The results from this work will be useful not only for constructing meso-informed macroscale models of HMX, but also for understanding the physics of void–void interactions and sensitivity due to void shape and orientation.

Appendices

Appendix 1: Governing equations

The hyperbolic conservation laws for mass, momentum, and energy are solved:

$$ \frac{\partial \rho }{\partial t} + \frac{{\partial \left( {\rho u_{i} } \right)}}{{\partial x_{i} }} = 0, $$

(28)

$$ \frac{{\partial \left( {\rho u_{i} } \right)}}{\partial t} + \frac{{\partial \left( {\rho u_{i} u_{j} - \sigma_{ij} } \right)}}{{\partial x_{j} }} = 0, $$

(29)

and

$$ \frac{{\partial \left( {\rho E} \right)}}{\partial t} + \frac{{\partial \left( {\rho Eu_{j} - \sigma_{ij} u_{i} } \right)}}{{\partial x_{j} }} = {\dot{\mathcal{E}}}. $$

(30)

The source term ${\dot{\mathcal{E}}}$ in (30) is the rise in specific internal energy of the system due to heat released in the decomposition of solid HMX into gaseous reaction products, and $\sigma_{ij}$ is the Cauchy stress tensor and is decomposed into volumetric and deviatoric component, i.e.,

$$ \sigma_{ij} = S_{ij} - p\delta_{ij} . $$

(31)

The deviatoric stress tensor, $S_{ij}$, is evolved using the following evolution equation.

$$ \frac{{\partial \left( {\rho S_{ij} } \right)}}{\partial t} + \frac{{\partial \left( {\rho S_{ij} u_{k} } \right)}}{{\partial x_{k} }} + \frac{2}{3}\rho GD_{kk} \delta_{ij} - 2\rho GD_{ij} + \rho S_{ik} \Omega_{kj} - \rho \Omega_{ik} S_{kj}= 0, $$

(32)

where $D_{ij}$ is the strain rate tensor $\Omega_{ij}$ is the spin tensor, and G is the shear modulus of material. First, the deviatoric stresses are evolved using an elastic response and then mapped back to the yield surface using a radial return algorithm [41]. The yield surface is given by the function $f = S_{\text{e}} - \sigma_{y}$, where $S_{\text{e}} = \sqrt {\frac{3}{2}\left( {S_{ij} S_{ij} } \right)}$. The yield strength, $\sigma_{y}$, is taken to be a constant and set to $260\;{\text{MPa}}$ [42] for HMX, i.e., hardening, and visco-plastic effects are neglected in the mesoscale computational models. A detailed description of the governing equations and radial return algorithm is provided in previous work [18, 19, 23,24,25].

The chemical species are evolved in time by solving the conservation equation:

$$ \frac{{\partial \rho Y_{i} }}{\partial t} + {\text{div}}\left( {\rho \vec{V} \;Y_{i} } \right) = \dot{Y}_{i} , $$

(33)

where $Y_{i}$ is the mass fraction of the $i{\text{th}}$ species and $\dot{Y}_{i}$ is the production rate source term for the $i{\text{th}}$ species. The numerical stiffness in solving the reactive set of equations is circumvented by using a Strang operator splitting approach [43], where first the advection of species is performed using the flow time step to obtain predicted species values:

$$ \frac{{\partial \rho Y_{i}^{*} }}{\partial t} + {\text{div}}\left( {\rho \vec{V}^{n} \; Y_{i}^{*} } \right) = 0. $$

(34)

In the second step, the evolution of the species mass fraction due to chemical reactions is calculated:

$$ \frac{{{\text{d}}Y_{i}^{n + 1} }}{{{\text{d}}t}} = \dot{Y}_{i}^{*n} . $$

(35)

The species evolution (35) is advanced in time using a fifth-order Runge–Kutta–Fehlberg method [44], which uses an internal adaptive time-stepping scheme to deal with the stiffness of the chemical kinetic equations.

Appendix 2: Constitutive and reaction models

The pressure in the mesoscale in (31) is obtained from a Birch–Murnaghan equation of state [45], which can be written in the general Mie–Gruneisen form as:

$$ p\left( {\rho ,e} \right) = p_{k} \left( \rho \right) + \rho \varGamma_{\text{s}} \left[ {e - e_{k} \left( \rho \right)} \right], $$

(36)

where

$$ p_{k} \left( \rho \right) = \frac{3}{2}K_{T0} \left[ {\left( {\frac{\rho }{{\rho_{0} }}} \right)^{ - 7/3} - \left( {\frac{\rho }{{\rho_{0} }}} \right)^{ - 5/3} } \right]\left[ {1 + \frac{3}{4}\left( {K_{T0}^{{\prime }} - 4} \right)\left[ {\left( {\frac{\rho }{{\rho_{0} }}} \right)^{ - 2/3} } \right] - 1} \right]. $$

(37)

Void collapse under shock loading can lead to the melting of HMX; therefore, thermal softening of HMX is modeled using the Kraut–Kennedy relation, $T_{\text{m}} = T_{{{\text{m}}0}} \left( {1 + a\frac{\Delta V}{{V_{0} }}} \right)$, with model parameters provided in the work of Menikoff et al. [42]. Once the temperature exceeds the melting point of HMX, the deviatoric strength terms are set to zero. Furthermore, the specific heat of HMX is known to change significantly with temperature. The variation of specific heat is modeled as a function of temperature as suggested in [42].

The chemical decomposition of HMX is modeled using a three-step mechanism proposed by Tarver et al. [46]. A detailed description of the implementation of the three-step model in the current numerical framework is presented in a previous work [7]. Here, a brief overview of the reaction model and its implementation is provided.

Chemical decomposition of HMX takes place in three steps involving four different species:

$$ {\text{Reaction}}\;1{:}\quad {\text{HMX}}\,({\text{C}}_{4} {\text{H}}_{8} {\text{N}}_{8} {\text{O}}_{8} ) \to {\text{fragments}}\, ({\text{CH}}_{2} {\text{NNO}}_{2} ) $$

(38)

$$ {\text{Reaction}}\;2{:}\quad {\text{fragments}}\, ({\text{CH}}_{2} {\text{NNO}}_{2} ) \to {\text{intermediate}}\;{\text{gases}}\, ({\text{CH}}_{2} {\text{O}},{\text{N}}_{2} {\text{O}},{\text{HCN}},{\text{HNO}}_{2} ) $$

(39)

and

$$ {\text{Reaction}}\;3{:}\quad 2 \times {\text{intermediate}}\;{\text{gases}}\; ({\text{CH}}_{2} {\text{O}},{\text{N}}_{2} {\text{O}},{\text{HCN}},{\text{HNO}}_{2} ) \to {\text{final}}\;{\text{gases}}\; ( {\text{N}}_{2} ,{\text{H}}_{2} {\text{O}},{\text{CO}}_{2} ,{\text{CO)}} . $$

(40)

The solid HMX (species 1, mass fraction $Y_{1}$) under high temperature decomposes into fragments (species 2, $Y_{2}$). The fragments are further decomposed to intermediate gases (species 3, $Y_{3}$) which are later converted to the final gases (species 4, $Y_{4}$) through exothermic reactions leading to high temperatures in the hotspot. In the absence of information about the equations of state for the intermediate and the final products, in this work, it is assumed that these intermediates and the products follow the same cold curves as the bulk material.

The change in temperature because of the chemical decomposition of HMX is calculated by solving the evolution equation:

$$ \rho C_{p} \dot{T} = \dot{Q}_{\text{R}} + k\Delta T , $$

(41)

where $\rho$ is the density of HMX, $C_{p}$ is the specific heat of HMX, $T$ is the temperature, $k$ is the thermal conductivity of HMX, and $\dot{Q}_{\text{R}}$ is the total heat release rate because of the chemical reaction. The values of $C_{p}$, λ, and $\dot{Q}_{\text{R}}$ are obtained from the work of Tarver et al. [46]. The source term in (30) is computed by setting ${\dot{\mathcal{E}}} = C_{v} \dot{T}$, where $C_{v}$ is the specific heat of HMX at constant volume.