Microstructural heterogeneity drives reaction initiation in granular materials

Granular compaction has many important applications including processing pharmaceutical tablets,¹ metal powder consolidation,² diamond synthesis,³ and energetic materials safety.⁴ The behavior of reactive granular materials under dynamic compaction has been the subject of numerous experimental and theoretical studies.^5–8 The added complexity in understanding reactive powder mixtures lies in the need to model and predict thermal localization in critical areas of the material, known as “hotspots,” where ignition occurs.

Although the physics of hotspots is not yet fully understood,⁹ observations confirm the hypothesis that reaction initiation occurs in the vicinity of macropores.¹⁰ This is consistent with the fact that areas with high plastic deformations cause large thermal dissipation and consequently extreme temperature fluctuations, which trigger chemical reactions. Other dynamic phenomena, such as shear banding¹¹ and elastoplastic deformation concentration due to material heterogeneities,⁹ are also likely mechanisms contributing to hotspot formation. Furthermore, the grain size distribution has been shown to play a major role in the sensitivity of hotspot formation.⁴ 

In light of this knowledge, phenomenological¹² and computational⁹ models of reaction initiation have been developed. The most common modeling approach is to use continuum models with thermodynamic closures and equations of state that are consistent with physics. Such continuum models account only for averaged microscopic dynamics, ignoring localization events that are due to microstructural defects. Since accurate physical measurements of reaction initiation are elusive, one has to rely on numerical modeling and simulation to study hotspot formation.

Hotspot formation in energetic materials has been extensively analyzed using (i) molecular dynamics simulations studying the thermochemical properties of hotspot formation, (ii) mesoscale pore collapse models of either a single pore or a few microstructural defects, characterizing the effect of the pore size and morphology on reaction initiation, and (iii) multiscale models that combine microscale dynamics with continuum models.

Molecular dynamic simulations^13,14 are essential to understand the underlying chemistry of hotspot formation. Yet, their computational cost renders them impractical for making quantitative predictions at the device scale.

Pore-scale numerical simulations have been used to characterize the criticality of hotspot formation.^10,15–17 The early work on isotropic (i.e., one-dimensional/radial) compaction of a viscoplastic spherical shell^18,19 has been extended to two- and three-dimensional simulations of a single macropore collapse.^15,20,21

Multiscale simulations incorporate this mesoscale dynamics into continuum models by adding source terms, such as the phase pressure difference⁸ or a reaction rate model,²² to the latter. Another multiscale approach²³ is to build a surrogate model from multiple mesoscale pore collapse simulations and use it as part of a continuum scale simulation.

While the current trend in modeling hotspot formation is to incorporate as much physics as possible, the predictive uncertainty of the resulting models (on any scale) is largely ignored. Specifically, continuum models are largely phenomenological and often neglect the subscale dynamics that is believed to be essential for hotspot formation. At any scale, models involve unknown parameters that are adjusted to fit experimental data, without necessarily reflecting the underlying physics. Although there is consensus that hotspots cause reaction initiation, they occur over wide time scales (10^–5–10^–3 s), space scales (≈0.1–10 μm), and temperature ranges (≥700 K).²⁴ This makes it impossible to obtain accurate measurements and reliably validate a physical hypothesis. Consequently, many questions about modeling assumptions remain unanswered.

Is a single solid phase viscoplastic Euler model accurate enough to model granular compaction,^22,25 or is the gas phase essential for capturing the dynamics of hotspot formation?²⁶ Do closure terms in the two-phase models honor the underlying physics? Are equilibrium thermodynamic assumptions appropriate in a physical problem that is (clearly) very far from equilibrium? How accurate is the nonconservative compaction law in, e.g., the Baer-Nunziato model? Is the atomistic scale relevant in studying hotspot formation?

The absence of direct observational data does not allow one to discriminate between alternative modeling hypotheses and, hence, to answer conclusively these basic questions. For instance, two models might predict different temperatures and reaction rates, but be equivalent within a range of irreducible uncertainties. Thus, increasing modeling complexity and improving numerical accuracy do not necessarily improve the predictive power when the model parameters are uncertain.²⁷ 

Yet, most material properties, e.g., porosity, pore size distribution, and heat capacity, are heterogeneous, and their point values are knowable only probabilistically. The values of some characteristics, such as viscosity and yield strength, are up to 100% uncertain in energetic materials like HMX.⁹ Solid surface properties are usually assumed the same across the solid-gas interface, often without experimental validation.

These ubiquitous structural (model) and parametric uncertainties argue for the use of Occam's razor: A model that honors experimental observations with the fewest modeling assumptions is preferable. Guided by this principle, we show that treating initial porosity as random is sufficient to capture hotspot formation using a standard continuum model described below. The Monte Carlo solution of this problem reveals the high sensitivity of the probability of reaction initiation to the heterogeneity of initial porosity, which is consistent with experimental observations.

Two-phase (solid/gas) flow models^5,28 are widely used to describe the dynamics of a granular material under dynamic compaction (Fig. 1). In their original formulation, these seven-equation models combine conservation laws for mass, momentum, and energy in each phase with a nonconservative compaction equation for the volume fraction.⁵ Subsequent studies provide simplified and regularized versions of the Baer-Nunziato equations,^6,29 found the relevant equations of state,^30–32 and added terms (e.g., “configurational stress”) to account for the granular microstructure.^4,28 Despite their usefulness, these equations are known to be numerically unstable.

A single velocity/single pressure approximation⁶ regularizes the equations and reduces the numerical instability.³³ The pressure-equilibrium assumption implies that the interface between the two phases is large. This five-equation model assumes the solid and gas phases to have the same velocity u and pressure p, giving rise to four conservation equations

where the subscripts s and g stand for the solid and gas phases whose densities are ρ_i (i = s and g) and volume fractions are α_i (i = s and g). The solid-gas mixture has density $ρ=αsρs+αgρg$⁠, internal energy $e=αsρses+αgρgeg$ (with e_s and e_g denoting the internal energy of the solid and gas phases, respectively), total energy $E=e+Bs(αs)+u2/2$ (with B_s denoting the configuration energy, a prescribed function of α_s, and the granular or configuration pressure defined by $βs=αsρsdBs/dαs$⁠), and pressure $p=αsps+αspg+βs$ supplemented with the stiffened gas equation of state for the i-th phase $pi=(γi−1)ρiei−γipi,∞$⁠, where γ_i and $pi,∞$ are known material dependent constants. The fifth equation comprising the model is a nonconservative compaction law

Hotspots are believed to occur in large pores where heat dissipation due to large plastic deformations and friction is the highest.³⁴ Therefore, the main quantity of interest (QoI) is the temperature at the solid/gas interface, rather than the average temperature $T¯(x,t)$ that can be deduced from the internal energy of the continuum model (1). The mesoscale temperature, T_μ, is governed by a reaction-diffusion equation

The Laplacian $κ∇2Tμ$ represents heat conduction into the solid bulk material. The rate of dissipation due to friction and permanent plastic deformation, $Q̇D$⁠, is, e.g., given by the deformation power.¹⁹ The exothermal reaction heat generation term, $Q̇R$⁠, follows an Arrhenius law.^12,22 In our analysis of hotspot formation, we focus on the dynamics before the onset of chemical reactions and neglect the reaction dissipation term $Q̇R$⁠.

Our QoI is the pore-surface temperature $Tps(x,t)=Tμ(xsg,t)$⁠, where $xsg$ is a point on the solid-gas interface. We expect that $Tps≫T¯$⁠, which can significantly change predictions of hotspot formation based on the continuum scale model. To account for the discrepancy between the continuum- and pore-scale temperatures, we adopt a reduced-order mesoscale heat model.⁸ This model redistributes the macroscale temperature $T¯(x,t)$ across the granular matrix by assuming a spherical shell mesoscale representation of the pores (Fig. 2).

The heat generated at the pore surface, q_c, depends on plastic deformation and, thus, on the rate of decrease in porosity $∂αg/∂t$⁠. For spherical pores, $αg=43πRg3Ng$⁠, where R_g is the pore radius and N_g is the number of pores per unit volume. The heat generated at the surface dissipates into the bulk of the solid phase by heat conduction, q_s. The resulting temperature profile is given by the solution of the diffusion equation over a spherical shell. This gives the pore surface temperature⁸ 

where λ_s is the thermal conductivity and the parameter n determines the shape of the temperature profile. The heat dissipation boundary layer δ is determined from the average temperature constraint $T¯=1/V∫RgδTdV$ (Fig. 2). The core temperature T_core is determined from the solid-phase Hugoniot relation

where $νs=1/ρs$ and the superscript 0 denotes the preshocked state.

An alternate measure of the probability of reaction initiation is the reaction rate $Ṙ$ related to the pore surface temperature T_ps via the Arrhenius rate law⁴ 

where the exponential prefactor Z and the activation temperature T^* are known for a given material; they are set to Z = 5 × 10^–19 s⁻¹ and T^* = 2.65 × 10⁴ K in our simulations.⁴ In this study, this quantity measures the “potential” average number of reactions triggered per unit time, without actually having reactions to take place. This assumes that the expected number of reactions per unit time is directly related to the probability of initiation at each point (x, t).

Microstructural heterogeneity in granular materials is the source of (i) fast temperature and stress fluctuations in the vicinity of large pores and (ii) an epistemic uncertainty in the key pore size and material property distribution. Accordingly, we treat the initial porosity, $αg(x,t=0)$⁠, as a spatially uncorrelated Gaussian random field, $αg(x,0)∼N(α¯g,σg2)$⁠, such that

where $⟨·⟩$ denotes the ensemble average. This approximation is valid for randomly packed grains with no spatial correlation at the grid size (∼10 μm). Spatially correlated microstructures inferred from real images can also be used.

We use (1)–(3), subject to this random initial condition, to estimate the probability of pore-surface temperature $Tps(x,t)$ exceeding a reaction initiation threshold T_ig

where $FTps$ is the cumulative distribution function (CDF) of the random surface temperature $Tps(x,t)$⁠. We start by assuming that reactions are instantaneously triggered when T_ps exceeds the known threshold T_ig, which represents a zeroth-order approximation of the probability of reaction initiation.

We use Monte Carlo simulations to estimate these probabilities. In every realization, (1)–(3) are solved with the initial condition $αg(x,0)$ drawn from the distribution $N(α¯g,σg2)$ at every grid-point with no spatial correlation. The simulation domain of length 1.6 cm is discretized with 200 grid points, and N_r = 600 realizations were performed. Every realization i yields a pore-surface temperature $Tpsi(x,t)$⁠. The ensemble average is approximated by the sample average $⟨Tps⟩≈⟨Tps⟩s=∑iTpsi/Nr$⁠, and the probability density function is constructed using a kernel density estimator $fTps(T̂;x,t)≈∑iKh(T̂−Tpsi(x,t))/Nr$⁠, where $T̂$ is a sample space variable and $Kh(·)$ is a Gaussian kernel with bandwidth h.

A Godunov-type HLLC numerical scheme³³ is used, with a velocity boundary condition, $vp≈60$ m/s, on the left boundary (Fig. 1) and a reflecting boundary condition on the right (although the simulation is stopped right before the shock reaches that boundary). The nonconservative and highly nonlinear nature of the compaction law (1e) requires special treatment to preserve the volume fraction positivity.

Figure 3 shows a snapshot of one realization of the two main variables at time t = 12 μs: the solid volume fraction $αs(x,·)$ and the pore surface temperature $Tps(x,·)$⁠. The temperature is higher, particularly at the shock front, for a higher microstructural heterogeneity corresponding to a larger σ_g. Comparison of the gas volume fraction before and after the shock reveals that once the grains are compressed, the porosity becomes relatively homogeneous. Hence, larger pores deform more than smaller ones, thus dissipating more energy. This is consistent with experimental observations.⁹ 

The average pore-surface temperature $⟨Tps⟩s$ is almost independent of the initial microstructural heterogeneity σ_g (Fig. 4, left). However, the average reaction rate $⟨Ṙ⟩s$ increases with σ_g (Fig. 4, right). The increased sensitivity of $⟨Ṙ⟩s$ to σ_g is due to the nonlinearity of the Arrhenius law, which gives exponentially higher weights to higher temperatures.

Both the average pore-surface temperature $⟨Tps(x,t)⟩s$ and the average reaction rate $⟨Ṙ⟩s$ are the highest along the line $x/t=vsh$⁠, where v_sh is the speed of the shock (Fig. 5). This is because the highest rate of thermal generation occurs where the change in porosity is the highest. Behind the shock front, the temperature decreases due to thermal diffusion in the bulk.

Figure 5 shows the average reaction rate and pore surface temperature. The line with the peak reaction rate and pore surface temperature coincides with the shock front where plastic deformations are the highest. The average temperature behind the shock front is higher than the initial temperature due to the increase in pressure after compaction (also see Fig. 3). Temperature fluctuations cause the reaction rate to grow exponentially due to the Arrhenius law. This explains the fluctuations in the sample mean of the reaction rate, which amplify rare extreme temperature fluctuations that are averaged out in $⟨Tps⟩s$⁠. The impact of the increase in temperature behind the shock front is insignificant compared to that of extreme fluctuations.

We posit that the reaction initiation occurs when the temperature T_ps exceeds a certain threshold T_ig and computing the probability of this event, $P[Tps(x,t)>Tig]$ using (6). This probability is shown in Fig. 6 as a function of T_ig (left) and $σ¯g$ (right), for an arbitrary point (x_h, t_h) along the shock front (chosen for illustration). This observation highlights the fact that hotspots are not a function of the average pore surface temperature, but of the tail of the distribution. Furthermore, the dependence on T_ig confirms the intuition that the probability of initiation increases as the ignition threshold is lowered, thus becoming easier to exceed.

The exceedance probability in Fig. 6 provides a partial description of the probability of reaction initiation. Another relevant statistics is the probability of T_ps remaining above the ignition threshold T_ig for a certain time interval $τ=t2−t1$ that is sufficient for the reactions to take hold. In other words, we are concerned with the probability of $Tps>Tig$ during a given time interval τ: $P[Tps(x,t1)>Tig,Tps(x,t2)>Tig]=Ft[τ|Tps>Tig]$⁠, for all x and all permutations of t₁ and t₂. The corresponding CDF in Fig. 7 shows that the time interval during which temperature T_ps remains above the threshold T_ig increases with the microstructural fluctuation coefficient of variation $σ¯g$⁠.

In summary, we demonstrated the dependence of reaction initiation in compacted granular materials on the initial pore size distribution using a fluctuating initial microstructure. Our results show that the amplitude of the fluctuations in the initial porosity affected the distribution of the resulting pore-surface temperature, motivating a probabilistic approach to studying hotspot formation. The corresponding reaction rate and probability of initiation are sensitive to the fluctuation amplitude. This study demonstrates that a probabilistic approach is the simplest and most intuitive way to study reaction initiation in granular materials, given a fundamentally uncertain and spatially heterogeneous material microstructure.