Theory of combustion in disordered media
Mauro Schiulaz Christopher R. Laumann Alexander V. Balatsky Boris Z. Spivak


INT-PUB-09-024



Department of Physics, Boston University, Boston, MA, 02215, USA

Institute for Materials Science, Los Alamos, NM 87545, USA
Nordita, Center for Quantum Materials, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

Department of Physics, University of Washington, Seattle, WA 98195, USA


The conventional theory of combustion describes systems where all of the parameters are spatially homogeneous. 
On the other hand, combustion  in disordered explosives has long been known to occur after local regions of the material, called “hot spots”, are ignited. 
In this article we show that a system of randomly distributed hot spots exhibits a dynamic phase transition, which, depending on parameters of the system  can be either first or second order. 
These two regimes are separated by a tri-critical point. 
We also show that on the qualitative level the phase diagram of the system is universal. 
It is valid in two and three dimensions, in the cases when the hot spots interact either by heat or sound waves and in a broad range of microscopic disorder models.


64.60.ah, 64.60.Kw, 64.60.Ht, 65.60.+a

Theory of combustion in disordered media
    Boris Z. Spivak
    December 30, 2023
========================================








§ INTRODUCTION
 


The microscopic mechanism of ignition of explosives is thermal in origin: it is related to exothermic reactions whose rates are quickly increasing functions of  temperature.
The theory of burning in uniform media was developed in Refs.  <cit.>.
Usually two regimes are considered:
in the detonation regime heat propagates supersonically via shock waves,
while in the deflagration regime  (“combustion”) it propagates subsonically as determined by heat diffusion, and can be described by a nonlinear heat conduction equation.

The conventional theory of combustion describes systems where the parameters of the material are spatially homogeneous.
Generally however, in solid explosives these parameters are random, sample specific functions of coordinates.
Moreover, the explosion in disordered explosives has long been known to occur after local regions of the material, called “hot spots”, are ignited by various processes
(see, e.g., Refs. <cit.>).
In this case most of the predictions of the theory of burning are quite different from the conventional case <cit.>.
If the temperature near hot spot i exceeds the critical value T_ci, a local explosion begins and the hot spot burns and releases energy Q_i.
This heat propagates through the neutral medium to other hot spots which it may in turn ignite.
If the hot spots are distributed sufficiently far from each other, the burning time of a hot spot is negligible compared to the time of intersite heat propagation and the dynamics are governed by the heat equation with pointlike sources:

    (∂_t-κ/C∂_ r^2)T( r,t)=∑_iQ_i/Cδ( r- r_i)δ(t-t_i)-T-T_0/τ

where t_i is the time of explosion of hot spot i. 

Generally, both T_ci and Q_i are random quantities.
We have shown in Ref. <cit.> that in the deflagration regime a two-dimensional system of randomly distributed hot spots exhibits a 
dynamic phase transition. 
Once started, an explosion either is able to propagate through the entire sample or stops after burning only a finite fraction of the system.
Depending on the parameters of the system, the phase transition can be either a first-order transition or a second-order transition, which are separated by a tricritical point.
In the second-order regime, the transition lies in the percolative universality class.
The microscopic structure of hot spots is currently under debate (see, for example, a corresponding discussion in Ref. <cit.>.) 
Different models yield different correlations between T_ci and Q_i. The analysis in Ref. <cit.> was performed for the case where the critical temperatures of the hot spots and the energies they release are proportional to each other Q_i∼ T_ci. 

In this article we argue that the results obtained in Ref. <cit.> are in a sense universal.
They remain qualitatively the same for a much broader class of disorder models (as characterized by the correlations between T_ci and Q_i).
They are also valid both in the three-dimensional (3D) case and in the case where a 2D layer of an explosive is embedded into a 3D neutral environment.
Perhaps most strikingly, they are valid even if the hot spots interact via  pressure waves rather than diffusive heat waves as we discuss next.

Under certain circumstances explosives can be ignited by the gentlest of mechanical blows.
This suggests that in these cases the ignition of the hot spots is controlled by  
local pressure rather than the local temperature (see, for example, Ref. <cit.>).  
That is, hot spot i is ignited when the local pressure reaches its critical pressure, P_ci, which is randomly distributed.  
Currently, the microscopic mechanism of ignition of hot spots by pressure is under debate. 
One possible mechanism is that the explosive contains gas bubbles which can be  compressed.   
This compression leads to rapid adiabatic heating which in turn ignites the bubble.
In the regime of pressure-controlled ignition, the hot spots interact via sound or weak shock waves.
This regime may be viewed as intermediate between deflagration and detonation.
We show below that in this case, on the  qualitative level, the picture of explosion obtained in Ref. <cit.> for the deflagration regime  remains the same: depending on the values of the parameters the system exhibits either a second- or first-order dynamical phase transition, separated by a tricritical point. 

The remainder of this paper is structured as follows. In Sec. <ref> we introduce a model in which hot spots are activated by heat (deflagration) and discuss its phase diagram. In Sec. <ref> we do the same for a model in which hot spots are activated by pressure only. In Sec. <ref> we describe in detail the numerical procedures we employ, and how the data are analyzed. Finally, in Sec. <ref> we summarize our results and discuss some still open questions.






§ COMBUSTION IN THE DEFLAGRATION REGIME



In the deflagration regime, the energy released by the hot spots propagates via heat diffusion.
The analysis that follows is explained in more detailed in Ref. <cit.>.
If the time it takes individual hot spots to burn is short compared to the time of the inter-hot spot heat propagation, then the combustion is described by the equation

    (∂_t-κ/C∂_ r^2)T( r,t)=∑_iQ_i/Cδ( r- r_i)δ(t-t_i)-T-T_0/τ

Here t_i is the time at which the temperature at the ith hot spot reaches its critical value T( r= r_i)=T_ci, and r_i is the position of the spot. 
C is the heat capacity per unit volume, and κ is the heat conductance. T_0 is the temperature of the environment, while the dissipation time τ sets the strength of the coupling to the environment. 
For simplicity of notation, we set T_0 = 0 for the remainder of the paper.

The heat released from the explosion of a single hot spot at the origin of a uniform medium propagates as,

    T(𝐫, t)    = Q/C/(4π D t)^d/2exp(-r^2/4Dt - t/τ),

where d is the dimensionality of space,   D = κ/C is the diffusion constant, and r is the distance from the origin.
This wave ignites hot spots at position 𝐫_i if the local temperature rises above the critical temperature T_ci.
At a given position r, the passing wave attains its maximum temperature

    T_max(r,Q)   =   2^d Q/(√(16 r^2 D τ+d^2 D^2 τ^2)-d D τ)^d/2
       ×   exp(-√(16 r^2+d^2 Dτ)/2√(Dτ)).

at time

    t^* = d τ/4(√(4/d^2r^2/l^2 + 1) - 1),

where l = √(Dτ) is the dissipation length of the system.

To be concrete we use a distribution function of  T_ci in the following form:

    P(T_ci)∝
    	T_ci^α(T_max-T_ci)^α    0≤ T_ci≤ T_max,
    
    	0    otherwise,

where we normalize the units of temperature such that T_max = 1 and have chosen α = 4 for all simulations, to provide the soft gap necessary for stability of the sample.

Generally, the values of T_ci and Q_i are correlated.  To illustrate the universality of the phase diagram we consider below three models for the correlation, wherein Q_i is either linearly, quadratically, or inversely proportional to T_ci:
 Q_i=BCa^3T_ci,
Q_i=Ba^3CT_ci^2/T_max, or Q_i=BT_max^2a^3C(T_ci+U)^-1. 
Here C is the heat capacity, a is the average distance between  the hot spots, U is a cutoff temperature, and B is a dimensionless constant characterizing the strength of the hot spot's explosions. 
 
 The relation between the dissipation length l  the average distance between the hot spots a, and the parameter B determines the dynamics of the burning.
If l≪ a, then each hot spot can receive heat only from its closest neighbors and the effects of heat accumulation from multiple explosions may be neglected.
It has been shown in Ref.  <cit.> that in this limit, the theory of the heat propagation in the system of hot spots can be reduced to   percolation theory.







Thus  for B>B_c, the explosion percolates through the sample, while for B<B_c  the explosion propagates up to the correlation radius

    ξ∼ |B-B_c|^-ν.

In this regime, the explosion activates only a finite number of hot spots, and

    N̅_exp∼1/|B-B_c|^γ

increases as a power law of B-B_c.  
The value of the exponents ν and γ are given by percolation theory.
 
At B>B_c the explosion propagates along the percolating cluster and burns a finite fraction of the system. We call this the explosive phase.
As a result the system exhibits a second order phase transition.

This picture changes for l≫ a: in this regime heat from multiple sources accumulates, and no mapping to percolation theory is possible: rather, the system displays a first-order transition, in which the explosion travels either to distances of order a only or through the entire sample. No diverging length scale can be identified at the transition point. The two regimes are separated by a tricritical point.







In Fig. <ref> we present phase diagrams of the system obtained by numerical simulations of the three different correlation models for the T_ci and Q_i.  
Thus we argue that on the qualitative level, this diagram remains the same for a broad range of correlations and in different dimensionalities.
 
One may also wonder whether the results would change if one considers the more realistic situation of a two-dimensional array of hot spots embedded in a three-dimensional environment. 
In this case, the surrounding environment absorbs heat from the system, thereby providing a source of dissipation even when τ=∞. 
This is not enough to make the transition become second-order, though: as we show in Sec. <ref>, in the absence of dissipation the transition is first-order even in this case (albeit the discontinuity is weaker than in the purely two dimensional model). Therefore, one needs to actively drain heat from the system, in order to observe the second-order transition.



§ COMBUSTION IN THE PRESSURE-MEDIATED REGIME
 


In this section we consider the case where the explosion of a hot spot creates a spherical sound wave which propagates through a neutral media igniting other hot spots. 
In the absence of the dissipation, the system of hot spots can be described by a sound wave equation for the pressure P( r,t):

    (∂^2_t -c^2∂^2_ r)P( r,t)= ∑_iQ_ie^-( r- r_i)^2/r_0^2δ(t-t_i).

As in the deflagration regime, t_i is the time at which the pressure P( r_i,t) is larger than the critical pressure P_ci, while c is the speed of sound in the medium and Q_i is the magnitude of the pressure impulse (rather than a heat released). 

Even at late times, sound waves are more sensitive to the initial profile of an impulse than heat waves are. 
At distances larger than the spot radius r_0 in 3D, a wave created by an explosion of a hot spot has the form of a spherical shell with width of order r_0.
In Eq. (<ref>), we take into account this spatial structure by using a Gaussian pressure impulse of width r_0  rather than a delta function.
This has the added advantage that it reduces simulation errors due to the discretization of space (which modifies the high-momentum dispersion of the waves). 

Equation (<ref>) describes spherical sound waves propagating from some localized sources. 
Such waves are well known to have a different character in two and three dimensions <cit.>. 
In d=3, assuming spherical symmetry around a hot spot located at r=0, we can rewrite Eq. (<ref>) as

    ∂^2 P/∂ t^2(r,t)=c^21/r^2∂/∂ r(r^2∂ P/∂ r).

This equation can be solved by defining ϕ(r,t)≡ rP(r,t), 





with

    P(r,t)=1/rϕ(r-ct),

where ϕ is a function, determined by  the initial conditions. For the Gaussian source term in Eq. (<ref>), we have ϕ(r-ct)=Qr_0^2/2cexp(-(r-ct)^2/r_0^2). 
Such a solution describes a way propagating outwards from the center, without changing its qualitative shape: in particular, only a Gaussian small residual pressure is left behind the wave front.
No such procedure is possible in two dimensions. Rather, it can be shown <cit.> that, in this case, an outgoing solution of the wave equation with radial symmetry can be written as

    P(r,t)=∫_-∞^ct-rdξϕ(ξ)/√((ct-ξ)^2-r^2),

where again the function ϕ is determined by the initial conditions. In this case, the wave has no well defined backward front: it can be checked that the pressure decays in time as

    P∼1/t

for r≪ ct, i.e., far behind the wave front.

Despite this difference, we find numerically that the phase diagram for the explosion model with sound waves is qualitatively similar to the one with deflagration. We use a distribution function of  P_ci analogous to the one used in the deflagration regime,

    P(P_ci)∝
    	P_ci^α(P_max-P_ci)^α    0≤ P_ci≤ P_max,
    
    	0    otherwise,

 where once again P_max = 1 and α = 4.

 In Fig. <ref> we present results of numerical simulations of Eq. (<ref>)  in both d=2 and d=3, with relations Q_i=BcP_ci/a, Q_i= Bc^2 P_ci^2/(a^2 P_max), and Q_i=BcP_max/a (P_ci+U)^-1, with the cutoff pressure U=P_max/2.  
 In all these cases, as a function of the coefficient B the system of the hot spots exhibits a first order phase transition.  

Equation (<ref>) does not take into account sound wave dissipation which is controlled by heat conduction. 
Since the propagating spherical waves have a single spatial scale (the width) of order r_0, the decay length of each wave is of order l∼ r_0^-2. 
Therefore, in this limit Eq. (<ref>) should be multiplied by a factor exp(-r/l).
If l≪ a, then a hot spot can be ignited only by waves generated by neighbors.  In this case, similarly to the case of the deflagration regime,  the problem can be reduced to a percolation problem, and the system exhibits a second-order dynamical phase transition. 
The only difference from the deflagration regime is that we have to introduce the critical pressures igniting the hot spots P_ci instead of critical temperatures. 
Thus, on the qualitative level, the phase diagram of the system is similar to that presented in Fig. <ref> for the deflagration regime.
This means that the structure of the phase diagram does not depend on the details of how energy is transported within the system, but rather is of universal character.

Also in this case, we have explored the case of a two-dimensional explosive embedded in a three-dimensional environment, finding that also in this case the transition is first order. This further confirms the universal character of our results.




§ NUMERICAL METHODS
 


In order to simulate disordered combustion, we follow and generalize the procedure described in Ref. <cit.>. We generate samples composed of a two- or three-dimensional cubic lattice of pointlike hot spots embedded in a passive medium with uniform diffusivity. 
The lattice spacing is a and we take periodic boundary conditions for all simulations, except for that of the 2D active system embedded in a 3D material, where we take periodic boundary conditions in the two lattice directions and absorbing boundary conditions in the transverse direction.













The heat released on exploding hot spot i, Q_i, is determined by T_c^(i) according to which of the three microscopic disorder models we are simulating:
Q_i = B T_ci,  Q_i = B T_ci^2, or Q_i = B (T_ci + U)^-1. 
The units are such that C = a = 1. 
We have chosen the value U=T_max/2, to cut off the unphysical divergence of Q for T_c→ 0 in model 3. 

For the pressure-mediated model, the samples are generated according to the above description simply replacing the critical temperatures, T_ci, by critical pressures, P_ci. 



 §.§ Numerical simulation of the deflagration model




For the diffusive model, we simulate the heat equation with relaxation, Eq. (<ref>), discretizing time in steps Δ t.  We take a regular lattice of hot spots, with lattice spacing a. 
For each run, we initialize the combustion by making a single randomly chosen site explode, releasing its energy Q_i. 
During each time step, the heat propagates according to Eq. (<ref>);
additionally, when the temperature at an unexploded site i exceeds its local critical temperature, i.e., if T( r_i,t)≥ T_ci, site i explodes and releases energy Q_i.
After exploding, the hot spot at site i is exhausted and will not explode again.
Explicitly, this is implemented by iterating the following procedure at each time step Δ t, for each site i.




  * If site i is still active and , the site explodes, 

    T(𝐫_i,t)→ T(𝐫_i,t)+Q_i,

and becomes exhausted. 


  * The temperature then relaxes diffusively,

    T(𝐫_i,t+Δ t)=T(𝐫_i,t)+∑_j=1^ZT(𝐫_j,t)/Z+1,

where Z is the coordination number of the lattice, and j labels the nearest neighbors of site i. This update rule fixes the value of Δ t:

    Δ t=(Z+1)^2/da^2/4π D.

The procedure is asymptotically correct although at short times it distorts the dynamics. 
However, it allows us to greatly increase the efficiency of the simulations, and we do not expect it to have any significant qualitative impact on the physics of the system. 
Its primary effect is to introduce a systematic deviation in the position of the transition point, but not to change its nature.


  * Finally, the heat dissipation acts:

    T(𝐫_i,t+Δ t)→ T_0+(T(𝐫_i,t+Δ t)-T_0)e^-Δ t/τ.




We end the simulation when the temperature field has decayed to the point that no further sites may explode.
Our primary measurement is then to count the number N_exp of sites which have exploded.


The numerical data are shown in Fig. <ref> for the three-dimensional deflagration model, in Fig. <ref> for the two-dimensional model, and in Fig. <ref> for the two dimensional layer imbedded in a three-dimensional medium.  Since N_exp remains finite in the inert phase, and becomes extensive in the exploded phase, we take the fraction of exploded sites N_exp/L^d as an order parameter for the transition.














  §.§.§ Scaling analysis
 


As our initialization procedure involves randomly sampling an initial site, the final number of exploded sites, N_exp, corresponds to the mean cluster size in the mapping to percolation. 
In order to confirm that the second-order transition lies in the percolative universality class, we perform a finite size scaling of our data, using the following scaling ansatz:

    N_exp(B,L)∝ L^γ/νf(L^1/νB-B_c/B_c),

where L is the linear size of the sample, and γ and ν are scaling exponents. 
In particular, ν governs the divergence of the correlation length ξ at the critical point, while γ describes the divergence of the mean cluster size. 

To compare our data with percolation theory, we collapse the data at different system sizes for a given disorder model using B_c as a fitting parameter.
For the exponents ν and γ we use the values found in the percolation theory literature: ν=4/3 and γ=43/18 for d=2, and ν=0.87±0.02 and γ=1.7±0.1 for d=3 <cit.>.
This procedure works for τ<τ_c: in the dissipative regime, the transition belongs to the same universality class as percolation theory. 
For τ>τ_c, on the other hand, this scaling procedure fails entirely: data at different system sizes simply fall on top of each other for B<B_c, and are proportional to L^d for B>B_c, with a discontinuous jump for B=B_c. The transition is first order.
The tricritical endpoint τ_c is shown as a red dot in Fig. <ref>.



In order to further confirm the robustness of these results, we also collapse data from across all microscopic disorder models. 
This procedure is performed separately for d=2 and d=3, and the results are shown in Fig. <ref>. 
In order to achieve inter-model collapse, two extra fitting parameters, p and q, are needed. 
The intermodel scaling ansatz is

    N_exp(B,L)∝ p L^γ/νf(qL^1/νB-B_c/B_c).

p and q are model dependent, but system size independent, and are of O(1) for all models (ranging between 0.5 and 4 in the data shown). 
This is remarkable, since the values of B_c can vary by orders of magnitude as a function of model, τ, and spatial dimension. 
For d=2 [Fig. <ref>, (a)], the collapse is good, while for d=3 [Fig. <ref>, (b)] it is still visible, but less remarkable. 
This is most likely because we are restricted to much smaller system sizes in d=3, and therefore finite size corrections to scaling are more important.



 §.§ Numerical simulations of the pressure-mediated regime


For the pressure mediated model, we simulate the wave equation, Eq. (<ref>), by discretization in complete analogy to the simulations of the diffusive model.

We choose the time step Δ t such that

    r≡aΔ t/c=1/√(d),

satisfies the well-known Courant-Friedrichs-Lewy condition,  which is necessary for the simulation of wave equations to be stable <cit.>.

Explicitly, this is implemented by repeating the following procedure at each time step Δ t, for each site i.



  * If site i is still active and , the site explodes, 

    P(𝐫,t)→ P(𝐫,t)+Q_i e^-(𝐫-𝐫_i)^2/r_0^2,

and becomes exhausted.


  * The pressure propagates according to the discretized wave equation,

    P(𝐫_i,t+Δ t)   =r^2∑_j=1^ZP(𝐫_j,t)+(1-Z r^2)P(𝐫_i,t)
       -P(𝐫_i,t-Δ t),

where Z is the coordination number of the lattice, and j labels the nearest neighbors of site i. 


The simulation is terminated when the pressure waves are too weak to activate any unexploded site.









The number of exploded sites N_exp is plotted in Fig. <ref> for the two- and three-dimensional sound waves model, and in Fig. <ref> for the two-dimensional layer embedded in a three-dimensional medium. Since no dissipation is present in the system, the transition is first order in all cases.



§ CONCLUSION




We have considered  disordered explosives containing a system of hot spots, which  can interact either via heat propagation or via sound waves.  We have shown that as a function of  values of parameters  the system exhibits either a first or second order dynamical phase transition, with a tri-critical point dividing the two regimes.  On the qualitative level the phase diagram of the system is universal.  
It is independent of the dimensionality of space, details of the correlation between T_ci and Q_i, and the mechanism of interaction between the hot spots. 
To further corroborate this point we also simulated the detonation of 2D explosive embedded into a 3D host in both the deflagration and the pressure-mediated  regimes. 
The result is again qualitatively consistent with the universal phase diagram.

We believe that the results presented in  Sec. <ref>  for the case of hot spots interacting via sound waves should extend as well to the case where the energy is mediated by weak shock waves, as these also propagate ballistically.
In a more general context our results can be applied to systems which exhibit thermal-instability driven avalanches (see, for example,  Refs. <cit.>).

One of the limitations of the presented results is the assumption that during the burning process the different parts of the sample do not significantly move apart.  We are planning to consider this aspect of the problem in future work.

We also note that such first-order transition has been reported numerically in the detonation regime as well <cit.>. However, the authors of Ref. <cit.> assume the detonation wave to be sufficiently strong to average over all local fluctuations of the system parameters: this can be seen as a mean-field version of the problem addressed in this paper.

Note added in proof - Additionally, a model similar to ours was used to study the dynamics of front propagation in inhomogeneous systems (see, e.g., Refs. <cit.>). These studies are carried in the region of the first-order transition, and do not investigate the nature of the transition itself.


A.B. and B.S. acknowledge financial support from Los Alamos National Laboratory, Contract No. XWA200. 
C.R.L. acknowledges support from the Sloan Foundation through a Sloan Research Fellowship and the NSF through Grant No. PHY-1656234. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.




