The Critical Conditions for Thermal Explosion in a System Heated at a Constant Rate

We have analyzed the condition needed for thermal explosion to occur in a solid sample when the temperature of the vessel walls is raised at a constant rate. We have developed a dimensionless model that allows its direct comparison with an isoperibolic system (constant vessel wall temperature). We have obtained an analytical expression for the critical condition as a function of the system parameters. Our solution takes into account reactant consumption and covers different geometries: thin film, finite and infinite cylinder. The critical condition has been validated with numerical simulations and experiments. We show that, compared to the isoperibolic system, thermal explosion is a little bit more difficult to achieve under constant heating conditions. Besides, we show that thermal explosion on submicrometric films is nearly impossible.


Introduction
Knowledge of the critical conditions for the occurrence of thermal explosion (sample mass, container geometry, activation energy, thermal conductivity, temperature) allows to determine: a) the chemical risk associated with the storage and transportation of hazardous materials or in running chemical reactors [1][2][3][4][5][6], b) the conditions for pyrotechnic reactions to occur [7][8][9], c) munitions cookoff temperatures [10] and d) in general the ignition condition in chemical engineering processes [11,12].
Furthermore, the so-called "combustion synthesis" takes advantage of local heating related to the heat evolved during the chemical reaction to obtain materials at low processing temperature [13,14]. For instance, synthesis of functional metal oxide thin-films via combustion synthesis has attracted great attention because it would allow the use of low-temperature substrates and it would be a promising route towards the development of large-area and low-cost printed electronics [15]. Unfortunately, under isothermal conditions or during slow heating ramps, fast heat dissipation to the substrate hinders thermal explosion by avoiding the local overheating needed for a thermal runaway to occur [16][17][18][19][20][21][22]. Several authors have suggested that, during rapid heating ramps, combustion may be achieved in thin films [14,16]. So, the conditions leading to thermal runaway during heating at a constant rate is of interest for many applications.
In [23][24][25][26][27][28] thermal explosion under dynamic conditions is analysed numerically and theoretically. The parameter range explored numerically is quite limited. Besides, to reduce the original system to an ordinary differential equation (ODE) system, analytical approaches neglect the temperature and transformation degree distribution. As we will see, this approximation is too crude to provide an accurate description of the critical condition. As Boddington et al. [29] noted for the constant heating case: "The development of a realistic model for the ignition experiment is a formidable problem." The aim of the present work is to provide a thermal runaway condition for a system that is heated uniformly at a constant rate. To be more precise, the reactant is placed in a vessel whose walls are kept at a temperature, T f , that is raised at a constant rate. We consider a homogeneous solid sample where heat dissipation occurs mainly through heat diffusion towards the vessel walls. We will show that this model correctly describes the observed behavior of one metalorganic precursor. We will use a new dimensionless system of equations that is equivalent to that introduced by Frank-Kamenetskii for an isoperibolic system (a system where T f remains constant with time) [30,31]. Thus, our approach allows to determine whether thermal runaway during constant heating is easier or not than during an isothermal temperature program. We will obtain a combustion condition that covers different geometries and a wide parameter range that accounts for most practical cases. Finally, we will analyze the possibility of combustion in thin films.

The model
The model is based on the classical theory for ignition and front propagation in solid samples [11,21,30,[32][33][34][35]. The heat balance is the result of two opposite effects: heat generation by chemical reaction and heat removal through thermal conduction. We neglect the effect of reactive gas depletion or the evolution of the system parameters during the reaction. These are reasonable approximations to predict the development of a thermal runaway but are not sufficient to accurately describe the whole reaction course. We also assume a homogeneous medium. The last assumption is also valid for heterogeneous systems provided that the time of heat exchange is much shorter than that of the chemical reaction [36]. For a cylindrical vessel, the evolution of the temperature is described by a two dimensional (2D) partial differential equation (PDE): where ρ is density, c is heat capacity, q is the reaction specific heat and λ is thermal conductivity; T is temperature, t is time, z and r are the axial and radial coordinates, respectively (see Fig. 1), and α(r,z,t) is the transformation degree (α = 0 untransformed, α = 1 totally transformed).
We assume that the initial sample temperature, T in, is uniform and low enough to ensure that the degree of transformation 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, that is raised at a constant rate b: The last boundary condition states that the heat flux is null at the top of the sample.
Actually, if the vessel temperature is below 600ºC, the heat lost by radiation is negligible when compared to the heat dissipated by conduction through the crucible walls. As for the heat evolved from the released gases, it is nearly compensated by the loss of heat capacity; so, it has a negligible effect on the heat balance [37]. Finally, due to symmetry, this boundary condition also describes a closed vessel full of reactant provided that H is the vessel half height.
Concerning the reaction, we assume first-order reaction kinetics to account for the reactant consumption and an Arrhenian temperature dependence [38,39]: 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.

Dimensionless system
The analysis of the system becomes much simpler by introducing suitable dimensionless parameters. Previous analytical studies have based their dimensionless model on the critical temperature of the isoperibolic system (static conditions) [24][25][26][27].
To derive a model that is formally similar to the isoperibolic system [21], we have looked for a different normalization.
The only difference with respect to the isoperibolic system is the second and third boundary equations: in the isoperibolic system the vessel temperature remains constant while in the present case it increases steadily at rate b. If overheating due to the reaction is negligible, the reaction time, t R , is determined by the vessel temperature; it is inversely proportional to the rate constant. Under constant heating conditions, it has been shown that the reaction time has the same functional dependence but on T Kis [40,41], where T Kis is the temperature at which the reaction rate is at its maximum, and is given by Kissinger's equation [42,43]: Since the role played by the Kissinger temperature, T Kis , for constant heating rate is similar to the constant vessel temperature for the isoperibolic system, it seems natural to use the same dimensionless temperature, θ, time, τ, and space coordinate, ξ, already defined for the isoperibolic system [30,31] but replacing T f by T Kis . For the particular 1D case of a thin film we thus define: where, t i The length and time scales d i and t i have a precise meaning for the isoperibolic system: d i is the width of the zone were the reaction rate is significant [30,31] and t i is the time elapsed before the thermal runaway [44]. With the above dimensionless variables, Eqs.(1) and (4) become for a thin film: where N = 0 with the boundary conditions: The initial conditions are, where we have set t = 0 when T Kis is reached. We have taken T in = 0 K so that the condition that initially the transformation rate is null is fulfilled.
Eqs. (10)- (13) reveal that the 1D system is fully described by three where d for a thin film is the height (d = H). We will call these parameters the Ahrrenius (ε), Todes (θ T ) and Frank-Kamenetskii (δ) parameters to keep the same nomenclature than for the isoperibolic system.

Generalization of the above equations is straightforward for other 1D
geometries: for a sphere (N = 2 and d = R) and for an infinite cylinder (N = 1, d = R).
The derivation of the 2D system is straightforward from Eqs. (1)-(4) and (7) but, for the sake of simplicity, in the analysis of the critical condition we have preferred to write down only the 1D model. Nevertheless, the 2D simulations presented below for a cylinder or radius R and height H have been obtained from the numerical solution of the dimensionless model. In this case:

Experimental and numerical results
To assess the ability of the numerical model to describe real situations, we have analyzed numerically and experimentally the decomposition of barium trifluoroacetate, . Preparation details are given in ref. [45]. The evolution of the reaction with the temperature is monitored by thermogravimetry (TG). TG measures the evolution of the mass of a sample when it is submitted to a controlled temperature program. The conversion fraction at any time can be easily calculated as [46]: where V is the system's volume, m in is the initial sample mass, m t is the mass at time t and m fin is the final mass.
TG analysis was performed with a Mettler Toledo model TGA851eLF thermobalance. Samples were placed inside uncovered alumina crucibles. The internal and external radii of the crucible were 2.4 and 3 mm, respectively. Measurements were corrected with a blank experiment carried out under identical conditions. Inside the furnace a gas flow rate of 50 mL/min was controlled by mass flow meters. High purity synthetic air was used. Water-saturated air was obtained by bubbling the carrier gas in water at standard temperature and pressure (25ºC, 1 atm). Thermal analysis experiments were performed at 20 K/min. The parameters used to simulate the thermal decomposition Ba(TFA) 2 are summarized in Table 1. All have been measured and the details are given in ref. [47].
We can see in Fig Since Ba(TFA) 2 decomposition is governed by the superposition of two different mechanisms [45], the assumption that the reaction is ruled by a first order reaction is quite inaccurate. Thus, we would expect our numerical model to provide a rather poor prediction of the time evolution of  . In Fig. 2 we have plotted simulated and experimental Ba(TFA) 2 transformation curves when the sample is heated at 20 K/min for different sample masses. Although our model fails to predict  accurately, it does provide a rather good description of the first stages of the reaction, which is sufficient for our purpose of establishing the required conditions for a thermal runaway to develop. In the following, we will predict the critical condition for a thermal runaway to occur.

Critical condition.
As already stated in Section 2, the dimensionless set of equations describing the evolution of a 1D system, Eqs. (10)- (12), depends only on three parameters; namely ε, θ T and δ, Eq. (13). They collect all the physical parameters concerning the chemical reaction, heat transport, size (through d) and heating rate (implicitly through T Kis ).
Depending on the particular ε, θ T and δ values, the system will undergo thermal explosion or not. Our aim is to find the boundary of the ε, θ T and δ region where combustion occurs. This boundary is described by a so called "critical condition" relating these three parameters. In the limit of very exothermic transformations and high activation energies (   T  and 0   ), the critical condition states that δ cr (δ in the boundary) equals a constant value that depends on geometry.
Two general analytical approaches are used to simplify the system. In the quasistationary approach, the subcritical regime is related to the existence of a quasistationary solution where the temperature evolution closely follows the heating rate.
The quasi-stationary solution allows deriving the spatial dependence of α and T in the subcritical regime. Conversely, in the non-stationary approach, the PDE is reduced to an ODE system. The resulting ODE system allows calculating the time evolution of the averaged α and T. Both approaches have been applied to the isoperibolic case in our former paper [21] and are applied here to the constant heating case. The technical developments are detailed in Appendices A and B.

Exact solution
We have numerically integrated the exact model, Eqs. (10)- (12), for the particular case of a thin film. We have calculated the evolution of the heat released by the reaction and the heat dissipated by thermal conduction towards the substrate. Their ratio is a good criterion to determine whether thermal runaway develops or not [21]. We have used this criterion because it can be applied to both the exact model and the non-stationary approach. In particular, for a given set of parameters θ T and ε, we look for the value of δ at the threshold of thermal runaway. The result is shown in Fig. 3 where the value of δ cr is plotted as a function of θ T for several values of ε.
We should emphasize that the range of parameters in Fig. 3 (10 < θ T < 10 4 , 0 < ε < 0.1) covers most practical situations. Most bond 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 Kis = 600 K, a variation of E A between 50 and 500 kJ/mol results in a variation of ε between 0.01 and 0.1. As for θ T the quotient q/c corresponds to the adiabatic temperature rise. The adiabatic temperature rise is typically around 1000 K [13,49]. However, in the case of metalorganic precursors, combustion has been observed in powders for q/c values of few hundred Kelvins [17,45,50,51]. If we assume that the adiabatic temperature rise is at least of the same magnitude than T Kis we obtain a lower bound for θ T of 10 (for a maximum value of ε = 0.1).

Non-stationary approach
The base of the existing analytical solutions is the set of ODE obtained after neglecting the temperature distribution in the system's volume [23,24,26]. Using this approach, in Appendix B we derive the ODE system, where C is a constant that depends on the system's geometry: 32  [30,53,54]. The initial conditions are those of Eq. (12).
In Fig. 3 we compare the critical value of Frank-Kamenetskii parameter, δ cr , obtained from the numerical integration of the exact model, Eqs. (10)- (12), with the value obtained from the ODE system, Eq. (15). From Fig. 3 it is apparent that the ODE system correctly describes the functional dependence of δ cr on θ T and ε but, in contrast with the isoperibolic case [21], it fails to deliver an accurate prediction of the δ cr value.
This is because, in the dynamic case, neglecting the temperature and degree transformation distribution is an approximation to rough to obtain an accurate prediction of the combustion threshold. Since the quasi-stationary model (Appendix A) takes into account the temperature distribution, we derive the critical condition combining analytical results of the stationary model with numerical results.

The critical condition
To simplify the fitting procedure, we take advantage of the results delivered by the quasi-stationary approach (See Appendix A) and the fact that, for any value of ε, δ cr tends to a constant value in the limit   T  . These asymptotic values have been fitted to δ = 4.392/(1+2ε), see inset of Fig. 3. The proportionality constant is the value of δ deduced for ε = 0 in Appendix A: Eq. (16) and the critical condition for an isoperibolic system [21,30] are nearly identical, the main difference is a factor 5 that is related to reactant consumption prior to thermal runaway. At the threshold, thermal runaway develops when the transformation rate at the top center of the sample reaches its maximum value, and this occurs for the Kissinger temperature. Due to the thermal gradient between the center and the vessel walls, the thermal runaway is achieved for a furnace temperature below the Kissinger one (see Appendix A): Since in most real system ε<<1, the furnace critical temperature approaches the Kissinger one. Besides, the initial reactant consumption makes thermal explosion a little bit more difficult to achieve under constant heating conditions for a given temperature (T f for the static case or T Kis for the dynamic case).
The shape of the δ(θ T ) curves can be fitted to the functional dependence where f(ε) has been expanded up to ε 2 term. The best fit to δ cr for a thin film has been obtained with f(ε) = 2(2+ε+30ε 2 ), see Fig. 4.a. Finally, the particular critical condition for the thin film can be generalized to the other geometries by replacing C TF by the corresponding C: This generalization is reasonable because it has been shown, from the non-stationary analysis of an isoperibolic system, that the contribution of geometry can be separated from the rest of parameters [21]. In the present case, the non-stationary analysis agrees approximately with the exact dependence of δ cr on θ T and ε; so, we consider that geometry can also be separated.
To check the accuracy of the critical condition, the prediction of δ cr delivered by Eq. (18) has been compared to the exact value obtained by numerical integration of the PDE system for two 1D geometries: thin film, Fig. 4a, and infinite cylinder, Fig. 4b. In all these cases, the prediction (dashed lines) departs less than 2% from the exact value.
For the isoperibolic system, the critical condition is also valid for the 2D geometry of a finite cylinder [21]. In that case, C must be simply replaced by: To test that this assumption works under constant heating conditions, the exact . The agreement is excellent.

Critical mass for a finite size cylindrical vessel.
Combining Eqs. (18), (19) and the definition of δ, Eq. (13), one can determine the sample critical mass: ) and t R is the reaction time, Eq. (5). Eq. (20) has been obtained assuming that the vessel containing the sample is not covered by a lid. The case of a sample that fills a closed vessel is also described by Eq. (20)

Combustion in thin films
Despite the claim that metal organic precursor films can be transformed into oxide films at low processing temperatures thanks to thermal explosion [15], we have been unable to reproduce these results so far [17]. In fact, experiments carried out at moderate heating rates (10-20 K/min) or at isothermal conditions for a number of oxide precursors have always shown that their thin films do not experience combustion despite their powders do [18][19][20]. Furthermore, our theoretical analysis of the isoperibolic case delivered a lower bound for the critical thickness around 400 μm [21].
To overcome this limitation, several authors [14,16] have suggested that combustion may be easier to achieve when a substrate with its precursor film is put on a hot plate set at the processing temperature, i.e., the critical thickness is thought to be lower at very high heating rates. Now we are able to solve this question.
Consider that combustion of a precursor powder is observed by TG when it is Substitution of θ T and t R by their definitions, Eqs. (5) and (13) Since cr  will be the same within a factor of the order of 1 (see the evolution of Fig. 5) when passing from a finite cylinder (powder inside a crucible) to a thin film, we can write the following relationship: Besides, T Kis scales as the logarithm of the heating rate [55,56], Combining Eqs. (24) and (25) Suppose that a powder undergoes thermal explosion at, say, b = 10 K/min and H = 500 μm. For this heating rate, typical values of are around 30 (this corresponds to T Kis = 600 K for an activation energy of 150 kJ/mol). Thus, a thin film of H' = 1 μm will undergo a thermal runaway at a substrate temperature around 1222 K when heated at a rate as high as 4 10 7 K/min. Therefore, low temperature synthesis of thin films is virtually impossible, even when films are heated at very high heating rates.

Conclusions
The partial differential equations (PDE) system describing heat transport during an exothermic reaction has been solved for 1D (infinite cylinder and thin film) and for 2D (finite cylinder) geometries when the vessel walls or substrate are heated at a constant rate.
It has been shown that, for large and small systems, ignition is produced close to its walls and at the center of the reactant volume, respectively, and that, when the reactant mass is small enough, no combustion occurs.
A systematic search of the critical condition of combustion has been undertaken for a very broad range of the three adimensional parameters (Ahrrenius, ε, Todes, θ T , and Frank-Kamenetskii, δ) that describe the system dynamics. The exact condition, δ cr (ε,θ T ), has been fitted to an analytical functional dependence that is valid for virtually all experimental situations. It has been shown that, in the limit of large activation energy (ε → 0) and high adiabatic temperature rise (θ T → ∞), the combustion critical condition is formally the same as that obtained for an isoperibolic system except for a factor 5 that indicates that, during constant heating conditions, a thermal runaway is more difficult to develop.
This analytical condition has allowed us to derive a formula giving the critical mass for a cylindrical reactor of finite height. This formula agrees with the exact solution of the PDE system and with experimental results carried out with Ba(TFA) 2 .
In the case of thin films, it has been shown that extremely high heating rates are necessary to achieve combustion. Under this high heating rates, combustion occurs when the substrate reaches a high temperature, so low temperature thermal explosion is virtually impossible in thin films.

Appendix A. Quasi-stationary approach.
To derive the quasi-stationary and non-stationary models, we define a new variable that is the difference between the sample and the furnace temperatures, The evolution of ϑ is obtained by substituting Eq. (A.1) into Eqs. (10)-(12). In the absence of a significant local overheating, the sample temperature will follow steadily the time dependence fixed by the furnace walls, i.e. b t T    / , and, consequently: ( where k is an integration constant that is determined from the first boundary condition, where subscript cr stands for the parameter value at the threshold of ignition, and C=0.878 and 2, and k cr =1.81 and 1 for a thin film and an infinite cylinder, respectively [21].
At the threshold, the thermal runaway will develop when the transformation rate at the center of the sample is at its maximum, i.e., when the center of the sample reaches the Kissinger temperature: Notice from Fig A.2 where B is a geometrical factor [21], Finally, if we substitute ϑ by θ we obtain the non-stationary ODE system,