# The Thermal Explosion Revisited

## 1 Introduction

We revisit in this Note the classical problem of a thermal explosion in a long circular cylindrical vessel containing an exothermically reacting gas at rest. It was formulated and solved under some assumptions by D. A. Frank–Kamenetsky ,  (see also ). In the original formulation it was assumed that the wall of the vessel is ideally conducting, so that the gas temperature at the boundary is equal to the temperature of the ambient medium. Such boundary condition made the problem cylindrically symmetric, and the symmetry essentially simplified the solution. In the present Note the problem is modified in the following way: the wall is partially isolated so that the symmetry is lost. The critical values of the radius of the vessel are determined numerically. Especially instructive results were obtained for the cases when the isolated parts of the boundary are distributed periodically with large angular frequency, and the isolated part of the boundary is large. The scaling laws for the critical values were found. For large angular frequencies it was found that there exists an axisymmetric core of the temperature distribution which occupies a major part of the vessel. The conditions in this core at the critical case were found to be subcritical.

## 2 Mathematical Problem Formulation

Assume that a gas at rest is enclosed in a long cylindrical vessel of radius $`r_0`$. An exothermic reaction is going in the gas with the thermal effect $`Q`$ per unit mass of reacted gas. For the reaction rate the Arrhenius law is assumed with the activation energy $`E`$. If the thermal effect and the activation energy are large, it can be shown (see ) that an ‘intermediate-asymptotic’ steady state regime is achieved. For this regime the gas consumption in the reaction can be neglected, and, from the other side, the temperature distribution in the vessel is steady. Applying the Frank–Kamenetsky large activation energy approximation, a non-linear equation for dimensionless reduced temperature $`𝕦`$ is obtained:

$$\mathrm{\Delta }u+\lambda ^2e^u=0,$$
(2.1)   
where

$$u=\frac{(TT_0)E}{RT_0^2},$$
(2.2)   
the Laplace operator $`\mathrm{\Delta }`$ is related to dimensionless variables $`\rho =r/r_0`$, $`\theta `$; $`r,\theta `$ are the polar coordinates. The constant $`\lambda `$ is

$$\lambda =\frac{r_0}{l},l=\left(e^{\frac{E}{RT_0}}\kappa RT_0^2c/QE\sigma (T_0)\right)^{1/2}.$$
(2.3)   
Here $`T`$ is the absolute temperature, $`T_0`$–the temperature of the ambient medium, $`R`$–the universal gas constant, $`\kappa `$–the molecular thermal diffusivity, $`𝕔`$–heat capacity of gas per unit volume, $`\sigma (T)`$ is the pre-exponential factor in the Arrhenius reaction rate expression: a slow function of temperature.

In the classical problem formulation it was assumed that the whole wall of the vessel is ideally heat conducting, so that the gas temperature at the boundary is equal to the temperature of the ambient medium. This gives a Dirichlet condition for the equation (2.1):

$$u(1,\theta )=0.$$
(2.4)   
The boundary value problem formulation under this condition is axisymmetric, and this was important for obtaining the analytic solution in an explicit form. D. A. Frank–Kamenetsky showed that the solution to the problem (2.1), (2.4) does exist for $`\lambda \lambda _{cr}=\sqrt{2}`$ only. Physically it means that a quiet, non-explosive proceeding of the reaction is possible only if the radius of the vessel is less than a critical one: $`r_0(r_0)_{cr}=\sqrt{2}l`$. This condition is known as the condition of the thermal explosion.

In the present Note the following modification of the problem (2.1), (2.4) is proposed: Only a fraction $`\alpha `$ of the wall is heat conducting, while the fraction $`1\alpha `$ is thermally isolated. The simplest formulation corresponds to a mixed problem (Figure 1,a)

$$\begin{array}{ccc}\hfill u(1,\theta )=0& \text{at}\hfill & 0\theta 2\pi \alpha \hfill \\ \hfill _\rho u(1,\theta )=0& \text{at}\hfill & 2\pi \alpha <\theta 2\pi .\hfill \end{array}$$
(2.5)   
More interesting is the case when the isolated part of the wall is not concentrated on a single arc, but is distributed periodically (Figure 1,b): the boundary $`\rho =1`$ is divided into $`N`$ segments

$$\frac{2\pi }{N}(1+k)\theta \frac{2\pi }{N}k,k=0,1,\mathrm{},N1.$$
(2.6)   
The fraction $`\alpha `$ of each segment is left heat conducting, while the fraction $`1\alpha `$ becomes isolated. In this case the mixed boundary condition at $`\rho =1`$ has the form:

$$\begin{array}{ccc}\hfill u(1,\theta )=0& \text{at}\hfill & \frac{2\pi }{N}k\theta \frac{2\pi \alpha }{N}+\frac{2\pi }{N}k,\hfill \\ \hfill _\rho u(1,\theta )=0& \text{at}\hfill & \frac{2\pi \alpha }{N}+\frac{2\pi }{N}k<\theta \frac{2\pi }{N}(k+1).\hfill \end{array}$$
(2.7)   
The central question addressed in the present Note is: What are the asymptotic laws for the critical radius if $`N\mathrm{}`$ and $`\alpha 0`$, i.e. the period of the boundary condition tends to zero, and the isolation is close to a complete one. Remember: for the complete isolation the critical radius is equal to zero.

## 3 The Numerical Method

In order to answer the questions posed above, we must numerically evaluate the critical values $`\lambda _{cr}`$ for fixed $`N`$ and $`\alpha `$. There are two special items in determining $`\lambda _{cr}`$: the singularity of the linearized equation (2.1):

$$\mathrm{\Delta }\delta u+\lambda ^2e^u\delta u=0$$
(3.1)   
at the critical value of $`\lambda =\lambda _{cr}`$, and the dependence of the values of $`\lambda _{cr}`$ obtained by discretization upon the number of points, $`n`$ used to discretize each dimension. Here $`\delta u`$ is the perturbation of the solution.

For fixed values of $`N`$, $`n`$, and $`\alpha `$, we determined a trajectory of solutions $`u`$ versus $`\lambda `$ by solving equation (2.1) with a Newton–Raphson method. This method requires the solution of the linearized equation (3.1) which becomes singular at $`\lambda =\lambda _{cr}`$. In practice it prevents us from approaching the critical point closely. Therefore to make an accurate determination of $`\lambda _{cr}`$ an extrapolation procedure was used. It was assumed that for $`\lambda `$ approaching $`\lambda _{cr}`$ a parabolic approximation is valid:

$$\lambda _{cr}^2\lambda ^2=C(uu_0)^2.$$
(3.2)   
The parameters $`\lambda _{cr}`$, $`C`$ and $`u_0`$ were determined to fit the last 10 points on the trajectory $`u`$ versus $`\lambda `$ approaching $`\lambda _{cr}`$. The approximation (3.2) happened to be satisfactory. Typically the fit (3.2) is accurate to a few parts in $`10^6`$.

The above procedure yielded an estimate for $`\lambda _{cr}`$ as a function of $`N`$, $`\alpha `$ and the number of discretization points $`n`$. It is natural to remove the dependence of $`\lambda _{cr}`$ of non-physical, computational parameter $`n`$. For this purpose another extrapolation was used. In our numerical approximations the second-order accurate discretizations of the operators was employed. If the boundary conditions were smooth, the solution $`u`$ would approach a limit with an error of the order of $`1/n^2`$. But the boundary conditions are non-smooth, the derivative of the solution is discontinuous at the boundary, and it causes the order of the approximation to decrease. Extensive numerical calculations have shown $`\lambda _{cr}^2`$ to vary linearly with $`1/n`$. Therefore, the following iterative procedure was used: the value $`\lambda _{cr}(N,\alpha ,n)`$ in (3.2) was calculated for three different values of $`n`$. These three values are then fit to a linear function $`a+b/n`$. If the fit was poor, as might happen if the values of $`n`$ were too small, the procedure is repeated with larger values of $`n`$ and so on, until a satisfactory fitting was obtained. The extrapolation $`n\mathrm{}`$ is simply the value of the coefficient $`a`$.

## 4 Results of the Numerical Analysis

The results obtained by numerical solution are represented on Figures 2–4. Three instructive properties are revealed.

(i) On the graph of Figure 2 the values of $`\lambda _{cr}^2`$ are presented for growing values of $`N`$ as the functions of $`1/\alpha `$. It is seen that the critical value $`\lambda _{cr}`$, i.e. the critical radius of the vessel for large $`N`$ is practically insensitive to $`\alpha `$ up to $`\alpha `$ very small. For small $`N`$ the dependence of $`\lambda _{cr}`$ on $`\alpha `$ is strong. Clearly, for any $`N`$, $`\lambda _{cr}=0`$ for $`\alpha =0`$, but it is instructive that for instance, for $`N=256`$ when only $`1/512`$ ($`0.2`$ percent) part of the boundary is heat conducting, the critical value of the radius is only 4 percent less than the critical radius for wholly heat conducting wall.

(ii) For large $`N`$, starting, say, from $`N=32`$, there exists an internal core $`0\rho \rho _{}`$ where the solution is close to axisymmetric one (see Figure 3). The value $`\rho _{}`$ was selected so that $`|u_{\mathrm{max}}(\rho _{},\theta )u_{\mathrm{min}}(\rho _{},\theta )|<10^4`$. Introducing the mean value $`u_{}=(u_{\mathrm{max}}+u_{\mathrm{min}})/2`$ we notice, that for $`0\rho \rho _{}`$ the solution is close to axisymmetric, so that the equation (2.1) and the boundary condition at $`\rho =\rho _{}`$ assume the form:

$$\frac{1}{\rho }\frac{d}{d\rho }\rho \frac{du}{d\rho }+\lambda ^2e^u=0,u=u_{}\text{ at }\rho =\rho _{}.$$
(4.1)   
Transforming the variables

$$u=u_{}+v,R=\frac{\rho }{\rho _{}}$$
(4.2)   
we reduce the problem to a classic one:

$$\frac{1}{R}\frac{d}{dR}R\frac{dv}{dR}+\mathrm{\Lambda }^2e^v=0,v(1,\theta )=0$$
(4.3)   
where $`0R1`$, $`\mathrm{\Lambda }^2=\lambda ^2\rho _{}^2e^u_{}`$. We calculate now the values of $`\mathrm{\Lambda }_{cr}^2=\lambda _{cr}^2\rho _{}^2e^u_{}`$ where $`\lambda _{cr}`$ is the critical value obtained in previous calculations. The graphs $`\mathrm{\Lambda }_{cr}^2`$ as the functions of $`1/\alpha `$ for different $`N`$ are presented on Figure 4. It can be seen that the values of $`\mathrm{\Lambda }_{cr}^2`$ are always less than $`2`$ (within the limits of our numerical accuracy). This means that the ‘rugged’ boundary layer near the wall $`\rho =1`$ controls the approach to criticality. The thickness of this layer is of the order of the length of the segment $`2\pi /N`$. The angular derivative $`_\theta u`$ in the boundary layer is large. The value $`u_{}`$ decreases as $`N`$ increases.

(iii) The intermediate power laws are observed for $`\lambda _{cr}^2`$ at large $`N`$ and small $`\alpha `$:

$$\lambda _{cr}^2=S(N)\alpha ^{t(N)}.$$
(4.4)   
The values of $`S(N)`$ and $`t(N)`$ for various $`N`$ are given in the Table.

## 5 Conclusion

Non-axisymmetric modification of the problem of thermal explosion in a cylindrical vessel is formulated. The boundary is partly isolated, and only partly ideally conducting. Special attention is paid to the case of periodic distribution of the isolated and conducting parts. The critical values of radius and other relevant properties are obtained numerically.

It is shown that for the period small in comparison with the vessel radius the critical value of radius of the vessel is practically insensitive to the relative size of the open area of the wall up to its very small values. The temperature distribution in the central core is axisymmetric and subcritical even at globally critical conditions: the criticality is due to a thin boundary layer near the wall where the temperature distribution is highly non-axisymmetric. Intermediate power laws are obtained for the critical radius as the function of the relative open area of the wall.

This work was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy under Contract DE-AC-03-76-SF00098, and in part by the National Science Foundation under Grants DMS94-14631 and DMS-2732710.

Table

| $`N`$ | $`S(N)`$ | $`t(N)`$ |
| --- | --- | --- |
| 32 | 2.03 | 0.10 |
| 64 | 2.03 | 0.055 |
| 128 | 2.02 | 0.028 |
| 256 | 2.01 | 0.015 |

Figure Captions

Figure 1. A fraction of the wall is isolated. (a) Isolated (5/6) and conducting (1/6) parts are connected. (b) Isolated and conducting parts are distributed periodically.

Figure 2. Dimensionless critical radius as the function of conducting fraction $`\alpha `$ for different angular frequencies. It is seen that at large frequencies the critical radius is practically $`\alpha `$-independent up to very small values of the conducting fraction $`\alpha `$.

Figure 3. The solution reveals an axisymmetric internal central core and ‘rugged’ boundary layer ($`N=32`$, $`\alpha =1/32`$).

Figure 4. The temperature distribution in the internal central core is a subcritical one: $`\mathrm{\Lambda }^2<2`$.