The Critical Condition for Thermal Explosion in an Isoperibolic System

Knowing the conditions for a system to undergo thermal explosion is of utmost importance for many applications. We present a critical condition that accounts for reactant consumption and covers most practical situations, including low activation energy reactions. Our solution applies to cylindrical reactors of any radius to height ratio. In the case of films, it is shown that thermal explosion is virtually impossible. We also introduce a new criterion to define the boundary of thermal runaway based on heat balance. This new definition of criticality allows us to check the accuracy of the nonstationary model to describe the critical condition. The non-stationary model is the base of most approaches.


Introduction
Knowledge of the critical condition to predict the occurrence of a thermal explosion is pertinent for a large number of applications. 1 It is useful when evaluating the chemical risk associated to the thermal stability of compounds or reacting mixtures and in preventing ignition during the storage and transportation of hazardous materials or in chemical reactors, [2][3][4][5][6][7] in determining the conditions for pyrotechnic reactions to occur, [8][9][10] to determine munitions cook-off temperatures, 11 in preventing thermal explosion in batteries 12,13 and, in general, in establishing the ignition condition for chemical engineering processes. 14,15 Over the last two decades, combustion has become a very versatile route for the synthesis of intermetallic or oxide powders and compacts. [16][17][18][19][20][21] Combustion synthesis takes advantage of the heat evolved during the chemical reaction of precursor substances to obtain materials that would otherwise imply high-temperature processing techniques. 16 Therefore, highly exothermic reactions of high adiabatic temperatures above 1800 K are needed. 16 Recently, combustion synthesis of functional oxides from precursor salts (nitrates or metalorganic precursors) has attracted much attention due to the growing interest in these materials. 19,22 The decomposition of oxide precursor salts may entail much lower adiabatic temperatures, i.e., of only several hundred kelvins. In other words, the range of thermodynamic and kinetic parameters involved in combustion synthesis has expanded considerably. One striking application is the combustion of bimetallic foils used for joining temperature-sensitive or dissimilar materials. 23,24 Also, a promising new low-temperature fabrication route to synthetize high-performance metal oxide thin-films based on combustion has recently attracted much attention. 25 In these particular applications, heat dissipation to the substrate is very important because it can hinder combustion by avoiding the overheating needed to set a thermal runaway. Consequently, film thickness is a critical parameter that has been analyzed both experimentally 24,26 and theoretically. 27,28 Knowledge of the critical condition for combustion is also important for intermetallic reactions and thermites (metal-metal oxides mixtures) that are used in some pyrotechnic devices such as igniters. These reactions evolve at a very high rate without forming a gas phase product. 29 We will analyze an isoperibolic system, i.e., the sample is placed in a vessel, the walls of which are kept at a constant temperature T f . The system undergoes thermal explosion, when reactants are heated to a temperature where the reaction becomes locally unstable. 18,19,27,30 Frank-Kamenetskii solution 31 is probably the best known critical condition for thermal explosion and it has been extensively applied to describe combustion in condensed matter. 14,32,33 In its derivation, reactant consumption is neglected and it is assumed that the activation energy is infinite. As we will see, the first assumption is accurate in some instances where thermal explosion arises at the very early stages of the reaction but for relatively low enthalpy reactions, such as the decomposition of oxide precursor salts, this may be a poor approximation. 34,35 A realistic model should take into account reactant consumption and finite values of the activation energy. 36,37 Frank-Kamenetskii and other authors improved the analysis by taking into account reactant consumption 31,[34][35][36][38][39][40][41][42][43] and finite values of the activation energy. 36,42 Analytical solutions are limited to idealized geometries that allow the system to be reduced to a one-dimensional (1D) model. A common approach to reduce a cylindrical vessel to a 1D model is to assume a spherical vessel with the same volume; 6 we will show that this is quite a crude approximation. Several authors have numerically explored two-dimensional models (2D) 44,45 but a general analytical solution for more realistic geometries is still lacking.
In this paper, we explore numerically the accuracy of existing critical conditions for a range of parameters that covers nearly all practical cases. We show that most models provide an accurate prediction for highly exothermic reactions and high activation energies, but they exhibit significant inaccuracies for low enthalpy or low activation energy reactions. We introduce minor modifications to the existing models to overcome these limitations. This result allows us to determine the critical thickness for thermal explosion to occur in a film. This analysis shows that thermal explosion is virtually impossible to achieve for thin-films. We also expand the critical condition to account for a more realistic vessel geometry. As a result of that, we obtain an analytical solution for the sample's critical mass needed for a thermal explosion to occur.
The structure of this paper is as follows. In the second section, we introduce a 2D model that accounts for heat transport by means of diffusion and heat generation due to a thermally activated exothermic reaction. In the third section we compare numerical results against experimental data to assess the ability of the numerical model to describe real situations. In the fourth section we introduce a new definition of criticality and we analyze the accuracy of several runaway critical conditions. Afterwards, we introduce an improved critical condition for 1D and 2D geometries, respectively. The accuracy of these critical conditions is tested against the numerical and experimental data introduced in the third section. The final section is a practical guide to determine whether a thermal runaway will occur or not.

The model
The classical theory for ignition 1,14,31,40,46,47 involves heat generated locally by a chemical reaction and heat propagation through the sample. We ignore the effect of reactive gas depletion or the evolution of the system parameters during the reaction.
These are reasonable approximations to predict the onset of a thermal runaway but may result in a poor prediction of the evolution of the reaction for those systems whose parameters evolve during the reaction. For instance, in the occurrence of gas depletion, the reaction front may extinguish. 48 Therefore, our analysis would provide necessary, but not sufficient condition, for a thermal runaway to take place. The model describes heat propagation in a homogeneous medium. It is a reasonable approximation to describe the macroscopic behavior of front propagation provided that the reaction front is wider than the medium inhomogeneities 1 or if the time of heat exchange in the heterogeneous medium is much shorter than the time of chemical reaction. 49 We also overlook heat losses from convection and radiation. Under these assumptions, and for a cylindrical vessel, the system is reduced to a 2D model where the rate of temperature change has two contributions, a conduction term and a reaction term where ρ is the density, q is the specific heat of reaction and λ is thermal conductivity, z and r are the axial and radial coordinates, respectively (see Fig. 1) and α(r,z,t) is the degree of conversion (α = 0 untransformed, α = 1 totally transformed). We assume that initially the sample is in thermal equilibrium with the furnace walls and that the degree of conversion is zero throughout the sample. The boundary conditions are: where H is the sample height, R is the vessel's inner radius and T f is the vessel temperature. The last condition states that the heat flux is null at the top of the sample (no convention, no radiation losses). One may think that the latter assumption is a very rough approach to deal with decomposition reactions where a significant amount of heat is lost through the gases evolved. However, when gases evolve, the heat capacity of the sample diminishes. If the gas is thermalized with the sample, the heat evolved is almost completely compensated with the loss of heat capacity, so the final heat balance is negligible with respect to heat losses through the crucible walls. 50 The set of boundary conditions, Eq. (2), also describes a closed vessel full of reactant provided that H is the vessel's half height (due to symmetry, for a cylindrical vessel the boundary condition As for the reaction, we assume a n-order reaction kinetics to account for the reactant consumption. Moreover, the reaction is thermally activated and we suppose an Arrhenius temperature dependence, i.e.: ( where A and E A are the pre-exponential constant and the activation energy of the reaction rate constant, respectively, and R G is the universal gas constant.
This model has been successfully applied to describe combustion in solids, 32 thermal explosion in munitions 11,51 and front propagation speed in pyrotechnic systems. 8 It has also been shown that it provides a reasonable description for those systems whose parameters evolve with temperature and can be used in solid-gas reactions provided that gas exchange is fast enough and that heat transfer through the gases is negligible. 52

Experimental and numerical results
To check the ability of the numerical model to describe real processes, we have analyzed the decomposition of a metalorganic precursor in the form of powders: yttrium trifluoroacetate, Y(CF 3 COO) 3 (Y(TFA) 3 ). Preparation details are given in. 53 The reaction course is monitored by thermogravimetry (TG). TG records the mass of a sample, m(t), when it is submitted to a controlled temperature program. The averaged conversion fraction at any time can be easily calculated as: 54 where m in is the initial sample mass and m fin is the final mass. The curves of the transformed fraction versus time for different initial sample masses are shown in Fig. 2. One can observe an evolution of the kinetics with the sample mass that is a characteristic feature of overheating due to the heat released by the chemical reaction. 28 Also, when the mass changes from 35 to 38 mg the curves develop an abrupt step that reveals the occurrence of a thermal runaway. 26,28,55,56 Conversely, all curves measured at 280ºC are smooth. Thus, there is a critical temperature and a critical size above which thermal runaway occurs.
The parameters used to simulate the decomposition of Y(TFA) 3 are summarized in Table 1. These parameters have been determined experimentally. 28 The activation energy and pre-exponential term have been determined from several TG measurements performed at different heating rates by means of isoconversional kinetic methods. 54,[57][58][59] In Fig. 3 we illustrate the calculated ) (t  for different temperatures and sample masses.
The evolution of the kinetics with temperature and sample size is reproduced by the numerical simulations. At 290ºC we observe that an abrupt step develops in the 41-42 mg mass range; thus the predicted critical mass is in fair agreement with that determined from the experiments. Also, in agreement with experiments, no thermal runaway is observed at 280ºC.
In contrast, the numerical simulation fails to provide an accurate description of the reaction course. The reason is that we have assumed a single-step first-order reaction model, which is a very rough approximation to the actual decomposition kinetics of Y(TFA) 3 . So, there is a fair agreement between the critical masses and temperatures determined from experiments and simulations despite the fact that the reaction course depends on the reaction model. This apparent contradiction is solved if we consider that the thermal runaway occurs at the early stages of the transformation; consequently, a first-order reaction model would correctly predict the threshold provided that it describes accurately enough the first stages of the reaction.
Also, it is important to note that isothermal experiments are far from ideal; the furnace does not reach the isotherm instantaneously but it takes more than 10 minutes to fully stabilize the temperature. Since, during this time lapse, the reaction progresses, a compromise has been taken between furnace thermalization and null initial sample transformation: the experiment onset is taken when the furnace temperature is 4ºC below the programmed one. Though, this strategy is far from compensating the artifacts due to deviations from the programmed temperature. Indeed, near the thermal runaway threshold, the reaction course is very sensitive to temperature fluctuations around the programmed temperature.
The occurrence of a thermal runaway is controlled by two competing phenomena: the exothermic reaction that tends to increase the local temperature, and heat dissipation that lowers the temperature through heat conduction. 1 The heat power released by the reaction is given by where V is the sample's volume. The heat power dissipated is, where the closed surface S is the boundary between the sample and the vessel.
Semenov 60,61 showed that, in the subcritical parameter region, a stable low temperature state exists where heat generation is compensated by heat dissipation, I Gen = I Dis . In homogeneously heated systems, thermal explosion occurs when this low temperature state is unstable, i.e., heat generated by the reaction exceeds the heat removed through the walls of the vessel, I Gen > I Dis . Hence, the occurrence of a thermal runaway can be revealed from the value of the ratio between the heat generation and the heat removal rates, ϕ: In Fig. 4 we have plotted the evolution of  for two different masses. For a sample below the critical mass (22 mg), after a transient period,  evolves to 1, i.e., it reaches the subcritical regime where heat generation is compensated by heat removal.
The transient period is due to the build-up of the temperature gradient related to heat removal. Conversely, for a sample above the critical mass (56 mg), after the transient period, the sample experiences an abrupt increase of heat generation that is responsible for the thermal runaway. Eventually, due to the reactant consumption, the heat generated diminishes. Thus, when thermal runaway occurs, the evolution of  goes through a local maximum.
An example of a simulation where thermal runaway takes place is given in Movie 1 and Fig. 5. Movie 1 shows the evolution of α, vertical and horizontal axes correspond to the z and r coordinates, respectively (see Fig 1); so, the top left corner corresponds to the top center point in the sample. Given that when thermal runaway occurs the decomposition rate increases considerably, frames are taken at constant degree of transformation intervals; otherwise, combustion would be like a flash in the movie.
In the supercritical region the system approaches the adiabatic temperature T AD , i.e., the temperature the sample would reach under adiabatic conditions, 31 One can observe that, locally, the temperature approaches the adiabatic one (T AD = 602ºC). Still, the maximum temperature is 30ºC below T AD . This difference is due to the low reaction enthalpy. Also, due to the local overheating, a temperature front is formed that propagates through the sample, this feature is clearer in Movie 1.
From Movie 1, we can observe that the combustion front is initiated at the top center and not at the bottom of the sample where it is in thermal contact with the furnace. The reason is that the vessel walls are a thermal sink that allows rapid heat dissipation thus preventing local overheating and, consequently, the formation of a combustion front there. For a sample mass near the critical value, the front is always initiated at the farthest location from the vessel walls, i.e., the top center of the sample, and propagates towards the vessel walls (for a closed vessel full of reactant the combustion front would initiate right in the center of the system).

Dimensionless systems
The model introduced in the previous section leads to a partial differential equation (PDE) coupled to a non-linear ordinary differential equation (ODE). As a result, thermal explosion is a complex system to treat both analytically and numerically. To simplify the analysis, simple geometries are assumed and the model is rewritten in a dimensionless form: 31,32 where N=0, 1 and 2 corresponds to a thin-film, an infinite cylinder and a sphere, respectively. The boundary conditions are: The dimensionless temperature θ, time τ and space coordinate ξ are defined as: ) where, t i is the thermal runaway induction time scale, i.e., a time scale for the time elapsed before the thermal runaway, 39 and d i is a scale of the width of the zone where the reaction rate is significant: 31,32 where d is a characteristic length. For a spherical vessel or an infinite cylinder it is the radius (d=R) and for a thin-film is the film thickness (d=H). Note that Todes parameter is the dimensionless adiabatic temperature.

Critical condition in the limit of no reactant consumption (θ T →∞)
Historically, two different approaches have been followed to simplify the model: the stationary and the non-stationary approximations. 31,32 In the stationary approach, the condition under which a stationary state is impossible establishes the critical condition. For a given symmetry the boundary problem can be solved to obtain the critical value. 31,[62][63][64] In particular, if we neglect the reactant consumption and we assume an infinite activation energy, the critical condition depends only on the Frank-Kamenetskii parameter (see appendix A): where subscript cr stands for the parameter value at the threshold of ignition, and C (1D) is a constant that depends on the system's geometry: for a spherical vessel, for an infinite cylinder and 878 . 0 for a thin film. The superscript (1D) stands for one-dimensional heat propagation. Above this critical value, C   , ignition occurs.
The non-stationary approximation consist in describing the evolution of a spaceaveraged temperature  and transformation degree  . It has been shown that this approach allows to reduce PDE systems that exhibit a low-dimensional dynamics to an equivalent ODE system. 65 Under this approximation, the dimensionless system, Eqs. (9) and (10), becomes: 2,31,40,66 where B is a constant that depends on the geometry of the system and the initial values of  and  are zero. Parameter B can be deduced by comparing the critical conditions delivered by the stationary and the non-stationary approximations in the limit case of infinite activation energy and no reactant consumption; Eq. (16) has no analytical solution and its dynamics is controlled by two time scales; a time scale related to heat generation, t i , and a time scale related to heat dissipation, t D (see appendix B): So the critical condition can be expressed in terms of the ratio Eqs. (15) and (18)

Criterion to state the occurrence of a thermal runaway
For θ T = ∞, the transition between the subcritical and supercritical states can be established rigorously (see Appendix C); in the subcritical regime the system evolves to a stable steady state, while in the supercritical regime, the temperature is unbounded and it increases steadily with time. Conversely, for finite values of θ T due to reactant consumption the system always converges to a stable steady state. Thus, when reactant consumption is taken into account, Semenov's criteria based on the boundedness of the supercritical state can no longer be used. 36,43,67 For highly exothermic reactions and high activation energies (θ T → ∞ and ε → 0) the separation between subcritical and supercritical regimes is sharp. The ignition temperature is well below the temperature of the supercritical state. In the subcritical state, the reaction heat is balanced by heat dissipation while in the supercritical state the reaction behaves nearly adiabatically: the heat evolved is some orders of magnitude larger than the heat dissipated. The time scale is also very different. It is some orders of magnitude smaller in the supercritical regime than in the subcritical regime. So, the particular criterion used to determine whether runaway occurs or not has little relevance.
However, things are different when the heat of reaction or the activation energy is relatively low (typically below 10 6 J/kg). Frank-Kamenetskii 31 showed that the difference between the ignition and furnace temperatures is approximately, Thus, for low activation energies, i T  will be high. On the other hand, for small values of the reaction enthalpy, the adiabatic temperature rise, AD T  , will also be small. Thus, in both cases the difference between the ignition and adiabatic temperatures will be relatively small and, therefore, the transition between the subcritical and supercritical regimes would be smooth. So, the critical condition for a thermal runaway depends on the particular criterion used to establish the boundary between both regimes. 42 Note, that  is the ratio between the two terms that determine the evolution of the spaceaveraged temperature in Eq. (16). This is as expected, because these two terms are directly related to the contribution of heat generation and heat dissipation.
Our criterion is independent of the particular choice of the scaling used to derive the dimensionless system. Moreover, our criterion is equivalent to the often used criterion of Adler and Enig (see Appendix D).
Now, we will analyze how the particular choice of the runaway criterion affects the critical condition. First, we will determine the parameter range to be explored to account for most real cases. Most bonding energies are between 50 and 500 kJ/mol (bond energies bellow 50 kJ/mol are characteristic of weak interactions such as Van der Waals bonding). Assuming T f = 600 K, a variation of E A between 50 and 500 kJ/mol results in a variation of ε between 0.01 and 0.1. Indeed there is no thermal runaway for ε values above 1/4 (Appendix C). 35 As for the Todes parameter, Morbidelli and Varma 73 showed that for first order reactions, no thermal runaway occurs for θ T values below 10 when ε > 0.15 (this limit goes down to θ T =4 for ε = 0). So we will limit our numerical analysis to ε values between 0 and 0.1 and θ T >10.
In Fig. 6  Morbidelli et al. 43 confirmed this result and extended its validity to parametric sensitivity based on any system control parameter. This conclusion is also apparent in

Accuracy of the non-stationary model
The non-stationary model is the starting point for most analytical approaches. So to check its accuracy we will compare it to the exact PDE system, Eqs.
(1)-(3). In the following, to determine δ cr we will use the energy balance criteria because it can be applied to both systems. The results are plotted in Figs. 7, 8 and 9. In Fig. 7 we have plotted δ cr as a function of θ T for ε = 0 while in Fig. 8 we have plotted δ cr as a function of θ T for θ T =100. In A noteworthy property of the ODE system, see Eq. (16), is that geometry only affects the numerical value of parameter δ cr by a constant factor 1/B. Therefore, the functional dependence on ε and θ T is independent of the particular geometry. To verify this property we will consider two (1D) geometries; a thin film and an infinite cylinder.
From Eq. (15) it can be stated that, in the limit θ T = ∞ and ε = 0, ) 1  Fig. 9 confirms this agreement for the whole range of of ε and θ T values studied. Thus, the contribution of the geometry can be separated from the rest of parameters. We will take advantage of this property to derive the runaway critical condition for a finite cylinder.

Accuracy of the existing critical conditions
Several authors have derived analytical critical conditions using approximate solutions of Eq. (16). For instance Thomas 34 derived a criterion that accounts for the reactant consumption (finite values of θ T ), where n is the reaction order. A formally identical solution was obtained by Kassoy et al. 35 Notice than in the limit of large values of θ T Eqs. (21) and (22) In Fig. 7 we have plotted the prediction of δ cr calculated from all these condition but Lacey's condition (it diverges in the limit ε → 0). From Two aspects affect the observed differences: the accuracy of the approximations assumed in the derivation of the model and the criterion taken to establish the occurrence of a thermal runaway. Note that despite the fact that all thermal runaway criteria converge to the same value for θ T = ∞ (Fig. 6), deviations of the different models are quite significant for θ T = ∞ (see Fig. 9). Hence, the main source of the deviations is related to the accuracy of the models. For this reason, the finest development is the one that gives the most accurate prediction.

Thermal runaway condition for 1D geometry: analysis of combustion in thin-films
Despite the fact that Babushok et al. critical condition is the most accurate model, deviations are still relevant for low reaction enthalpies or low activation energies. Also, the previous critical conditions exhibit a similar dependence on θ T that qualitatively reproduces the observed dependence on θ T . As for the dependence on parameter ε, it can be shown, (see Appendix C) that in the limit of θ T → ∞ the critical condition is approximately    To check the accuracy of our approximation, we have calculated the percent error in the calculation of δ cr . The result is plotted in Fig 10 where it is apparent that our critical condition delivers an accurate prediction within the whole range of parameters; the relative error in the calculation of the critical thickness is less than 5% and below 3% for ε < 0.05. An inaccuracy below 3% is negligible, especially when the approximate nature of the model itself and the uncertainties of the system parameters are taken into account.
It is well-known that a critical size exists below which no thermal runaway occurs. 1,61 Now, we are in the position to evaluate the critical thickness for thermal explosion to occur in films when they are heated homogenously. In particular, we will analyze the possibility to synthetize metal oxide thin-films via volume combustion synthesis. From Figs. 6 and 7 it is apparent that the value of δ cr is between 1 and 2.
Once we now δ cr , we can determine the critical film thickness from Eq. (14): where a is the thermal diffusivity ( c a    ) and t R is a reaction time scale, Now, we are going to estimate the physically reasonable lower bound of the critical thickness. For solids K J/m 10 3 Combustion criterion for a finite size cylindrical vessel.
As we have seen in section 4.4, the contribution of geometry can be separated from that of the rest of parameters. Thus, Eq. (25) is still valid for a 2D geometry but we need to determine the geometrical factor for a finite cylinder, C (2D) .
The contribution of the geometry is related to heat transport through diffusion.
Diffusion can be approximately addressed by considering that two paths compete in parallel; an axial path and a radial path. Thus the diffusion length is, and the diffusion time for a finite cylinder is given by, Combining Eqs. (17) and (29)  gives an accurate estimation of the geometrical factor for a finite cylinder. Note, that as expected, this parameter depends only on the cylinder aspect ratio, H/R. Also, in Fig. 11 we have plotted the usual approximation of assuming a sphere of the same volume, so that the 2D model is reduced to a 1D model. 6 This approach is equivalent to assume a constant C (2D) = 3.32. From Fig. 11  The steps to determine if a thermal runaway will occur are summarized in Figure   12.
Finally, from Eqs. (25), (28) and (31) Eq. (32) has been determined assuming that the vessel containing the sample is not covered by a lid. The case of a sample that fills a closed vessel is equivalent but with H being half the sample height.
We have applied Eq.  We have extended the applicability of this critical condition to a cylindrical vessel without any restriction on its aspect ratio. This approach leads to a prediction of the minimum sample mass or sample size.
 The sample aspect ratio plays a crucial role since thermal explosion is harder to achieve when one of the dimensions vanishes. In particular, it is concluded that thermal explosion of submicrometric films is virtually impossible.
 Numerical simulations were compared to experiments carried out on metalorganic powders that decompose inside a crucible. Although numerical simulations fail to describe accurately the complete reaction course, they are able to predict the temperature and mass thresholds for thermal explosion.
Appendix A. Stationary approach.
As we have seen, in the subcritical parameter region, a stable low temperature state exists where heat generation is compensated by heat dissipation: so, the time derivative of the temperature in Eq.
where k is an integration constant. The value of k can be determined from the boundary condition, Thermal explosion occurs when there is no solution for the stationary state. From shown that this particular equation has been solved in astrophysical problems and that the critical condition for a sphere is δ cr = 3.32.

Appendix B. Non-stationary approach.
Under the non-stationary approach and assuming 1  Appendix C. Approximate critical condition in the limit θ T → ∞.
In the limit case θ T → ∞, Eq. (16) reduces to: 31,40 Eq. (C.1) is equivalent to neglecting reactant consumption. Eq. (C.1) has two different asymptotic solutions: a stationary one that corresponds to the subcritical solution and a divergent one that is related to the supercritical regime. Therefore, the combustion criterion corresponds to the boundary between these two regimes. This rigorous criterion cannot be applied in the case of finite values of θ T because reactant consumption ensures that the system will always evolve to a stationary solution.
The stationary solution of Eq. (C.1) is given by Likewise in Appendix B, the function ) ( f has a local minimum (see Fig. B (C.5) If we substitute Eq. (C.4) into Eq. (C.5) and we perform a Taylor series expansion for small values of ε we obtain: And finally, in a first order approximation the critical condition is: (C.7)

Appendix D. Equivalence between energy balance and Adler and Enig criteria.
According to Adler and Enig criterion, 2 ϕ Ratio between the heat generation and the heat dissipation rates, Eq. (7).           Thermal diffusivity, m 2 /s 1.3×10 -5