Thermal decomposition analysis in a sphere of combustible materials

Abstract

In this article, we look at spontaneous combustion due to exothermic chemical reaction taking place within a stockpile of combustible material. The model includes mass and energy balance equations in a spherical domain. The complicated chemical reaction is simplified by considering a one-dimensional process. The differential equations governing the problem are solved using semi-implicit finite difference method. The effects of kinetic parameters embedded within the system are analyzed and the results are expressed graphically and discussed accordingly.

Introduction

Spontaneous combustion taking place in a stockpile of combustible material is due to oxygen trapped within the stockpile reacting with the carbon containing material. The stockpile of combustible material includes solid fuels such as coal, cotton, hay, wood, wool, and other carbon or hydrocarbon containing materials. Figure 1 illustrates how self-ignited fires in stockpiles of hay can cause hazards to fauna and flora.

Figure 1.

Self-ignited cylindrical hays.

Physical factors or parameters that influence spontaneous combustion or self-ignition process, to mention just a few, are combustible material particle size, volume-to-surface ratio, porosity, thermal conductivity, density, and heat capacity including convection in the surrounding.^1,2 The exothermic chemical reaction or low-temperature oxidation causes emission of heat and greenhouse gases such as carbon dioxide and carbon monoxide.^3,4 Initially, the temperature of the system is assumed to be equal to that of the ambient, but as the exothermic chemical reaction continues the temperature of the system increases too. If the heat generated is not appropriately transferred to the surroundings, the temperature of the system increases such that the rate of reaction increases exponentially according to Arrhenius law.^5,6 When the heat released due to exothermic chemical reaction is not dissipated to the surrounding environment, the reaction rate becomes faster and this may result in thermal runaway phenomena to a point where spontaneous ignition occurs.^7,8 The heat produced in a stockpile of combustible material determines temperature gradients across the stockpile. These temperature gradients cause existence of natural convection, which pushes air into the stockpile to supply oxygen for the chemical reaction, but at the same time it cools the system with air at ambient temperature. ⁹ The chemical reaction taking place in this case is complicated and it includes many radicals.^10,11 The chemical kinetics model of the reaction of oxygen with carbon containing materials is reviewed in Williams, ¹² and this process involves nonlinear interactions of reacting species. The process is described by nonlinear differential equations that govern the chemical reaction of the system. The governing equations are solved numerically, and in this case, a semi-implicit finite difference method (FDM) is used. A one-step decomposition kinetics is applied to solve the complicated chemical reactions of the system as indicated in Legodi and Makinde ¹³ and Sadiq and Merkin. ¹⁴ In this article, the investigation of self-ignition in a spherical domain of combustible materials is investigated. A similar study was investigated in a reactive slab in Makinde et al., ¹⁵ and in Chinyoka and Makinde, ¹⁶ the study was done using a long cylindrical pipe. In section “Mathematical formulation,” we consider the mathematical formulation and the relevant numerical methods are described in section “Numerical method.” Detailed discussion of results and graphical solutions are presented in section “Results and discussion.”

Mathematical formulation

A spherical domain of combustible material with constant thermal conductivity k is considered. A one-step finite irreversible Arrhenius kinetics with convective heat loss to the ambient is assumed, and it is described by the formula^8,15,16

$$C_i{H}_j+(i+\frac{j}{4}){\mathrm{O}}_2\to i\mathrm{C}{\mathrm{O}}_2+\left(\frac{j}{2}\right){\mathrm{H}}_2\mathrm{O}+\mathrm{heat}$$

(1)

Convective heat loss at the surface of the sphere follows Newton’s law of cooling and it is expressed by $-(h/k)[T-T_b]$ , where h is heat exchange/transfer coefficient, and it is a measure of how many calories flow across the boundary per unit temperature difference per second per centimeter, k is the thermal conductivity of the sphere and it is also a measure of how well the sphere conducts heat. T is the sphere temperature and $T_b$ is the ambient temperature. Figure 2 illustrates the geometry of the problem.

Figure 2.

Geometry of the problem.

Following^15,16 the nonlinear partial differential equations (PDEs) describing temperature, oxygen concentration and carbon dioxide emission in the combustible material can be written as

$$\begin{array}{l}\rho c_p\frac{\partial T}{\partial \overline{t}}=k\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \overline{r}}\left({\overline{r}}^2\frac{\partial T}{\partial \overline{r}}\right)+QA{\left(\frac{KT}{vl}\right)}^m{\left(C-C_0\right)}^{n}exp\left(\frac{-E}{RT}\right)\\ \frac{\partial C}{\partial \overline{t}}=D\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \overline{r}}\left({\overline{r}}^2\frac{\partial C}{\partial \overline{r}}\right)-A{\left(\frac{KT}{vl}\right)}^m{\left(C-C_0\right)}^{n}exp\left(\frac{-E}{RT}\right)\\ \frac{\partial P}{\partial \overline{t}}=\gamma \frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \overline{r}}\left({\overline{r}}^2\frac{\partial P}{\partial \overline{r}}\right)+A{\left(\frac{KT}{vl}\right)}^m{\left(C-C_0\right)}^{n}exp\left(\frac{-E}{RT}\right)\end{array}\}$$

(2)

with initial and boundary conditions as follows

$$\begin{array}{c}T(\overline{r},0)=T_0,C(\overline{r},0)=0.5C_b,P(\overline{r},0)=0\hfill \\ \frac{\partial T}{\partial \overline{r}}(0,\overline{t})=\frac{\partial C}{\partial \overline{y}}(0,\overline{t})=\frac{\partial P}{\partial \overline{y}}(0,\overline{t})=0\hfill \\ \frac{\partial T}{\partial \overline{r}}(a,\overline{t})=-\frac{h_1}{k}[T(a,\overline{t})-T_b]\hfill \\ D\frac{\partial C}{\partial \overline{r}}(a,\overline{t})=-h_2[C(a,\overline{t})-C_b]\hfill \\ \gamma \frac{\partial P}{\partial \overline{r}}(a,\overline{t})=-h_3[P(a,\overline{t})-P_b]\hfill \end{array}$$

(3)

where T is the absolute temperature, C is the oxygen concentration, P is the carbon dioxide emission concentration, $T_b$ is the ambient temperature, $C_b$ is the oxygen concentration in the surrounding air, $P_b$ is the carbon dioxide concentration in the surrounding air, $\overline{t}$ is the time, $T_0$ is the slab initial temperature, $C_0$ is the initial concentration of oxygen in the sphere, $\rho$ is the density, $c_p$ is the specific heat at constant pressure, D is the diffusivity of oxygen in the sphere, $\gamma$ is the diffusivity of carbon dioxide in the sphere, Q is the exothermicity, A is the rate constant, E is the activation energy, R is the universal gas constant, l is the Planck number, v is the vibration frequency, K is the Boltzmann constant, $a=\overline{r}$ is the radial distance, $h_1$ is the coefficient of heat transfer between the sphere and its surroundings, $h_2$ is the coefficient of oxygen transfer between the sphere and its surroundings, $h_3$ is the coefficient of carbon dioxide transfer between the sphere and its surroundings, n is the order of exothermic chemical reaction, and $m\in \{-2,0,0.5\}$ is the numerical exponent such that $m=-2$ represents the sensitized kinetics (laser/light induced), Arrhenius kinetics is represented by $m=0$ and $m=0.5$ is for bimolecular kinetics, as indicated in previous studies.^3,4,8,15,16 The following dimensionless parameters are introduced into the set of equations (2) and (3)

$$\begin{array}{l}\theta =\frac{E(T-T_0)}{R{T}_0^2},{\theta }_b=\frac{E(T_b-T_0)}{R{T}_0^2},\mathrm{\Phi }=\frac{C-C_0}{C_b-C_0},\mathrm{\Psi }=\frac{P}{P_b}\\ B{i}_1=\frac{a{h}_1}{k},B{i}_2=\frac{a{h}_2}{D},B{i}_3=\frac{a{h}_3}{\gamma }\\ {\beta }_1=\frac{\rho c_{p}R{T}_0^2}{QE{C}_b},\lambda ={\left(\frac{KT}{vl}\right)}^m\frac{QAE{a}^2{C}_b^n}{kR{T}_0^2}exp\left(\frac{-E}{RT}\right)\\ r=\frac{\overline{r}}{a},t=\frac{k\overline{t}}{c_p\rho a^2},\epsilon =\frac{R{T}_0^2}{E},\alpha =\frac{D\rho c_p}{k},\sigma =\frac{\gamma \rho c_p}{k}\end{array}\}$$

(4)

to obtain the following dimensionless governing equations

$$\begin{array}{l}\frac{\partial \theta }{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial T}{\partial r}\right)+\lambda (1+\epsilon \theta {)}^m{\mathrm{\Phi }}^{n}exp\left(\frac{\mathrm{\theta }}{1+\epsilon \theta }\right)\\ \frac{\partial \mathrm{\Phi }}{\partial t}=\alpha \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial \mathrm{\Phi }}{\partial r}\right)-\lambda {\beta }_1(1+\epsilon \theta {)}^m{\mathrm{\Phi }}^{n}exp\left(\frac{\mathrm{\theta }}{1+\epsilon \theta }\right)\\ \frac{\partial \mathrm{\Psi }}{\partial t}=\sigma \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial \mathrm{\Psi }}{\partial r}\right)+\lambda {\beta }_2(1+\epsilon \theta {)}^m{\mathrm{\Phi }}^{n}exp\left(\frac{\mathrm{\theta }}{1+\epsilon \theta }\right)\end{array}\}$$

(5)

The corresponding initial and boundary conditions then become

$$\begin{array}{c}\multicolumn{1}{c}{\theta (r,0)=0,\Phi (r,0)=0.5,\Psi (r,0)=0}\\ \multicolumn{1}{c}{\frac{\partial \theta }{\partial r}(0,t)=\frac{\partial \Phi }{\partial r}(0,t)=\frac{\partial \Psi }{\partial r}(0,t)=0}\\ \multicolumn{1}{c}{\frac{\partial \theta }{\partial r}(1,t)=-B{\mathrm{i}}_1[\theta (1,t)-{\theta }_b]}\\ \multicolumn{1}{c}{\frac{\partial \Phi }{\partial r}(1,t)=-B{\mathrm{i}}_2[\Phi (1,t)-1]}\\ \multicolumn{1}{c}{\frac{\partial \Psi }{\partial r}(1,t)=-B{\mathrm{i}}_3[\Psi (1,t)-1]}\end{array}\}$$

(6)

where $\lambda$ is the Frank-Kamenetskii parameter, r is the sphere dimensionless radial distance, $\epsilon$ is the activation energy parameter, ${\beta }_1$ is the oxygen consumption rate parameter, ${\beta }_2$ is the carbon dioxide emission rate parameter, $\alpha$ is the oxygen diffusivity parameter, $\sigma$ is the carbon dioxide diffusivity parameter, and $B{i}_1$ , $B{i}_2$ , and $B{i}_3$ represent the thermal Biot number, oxygen Biot number, and carbon dioxide Biot number, respectively. The dimensionless heat and mass transfer rates at the sphere surface are expressed in terms of Nusselt and Sherwood numbers, respectively, as^3,4,10

$$\mathrm{Nu}=-\frac{d\theta }{\mathrm{dy}},S{h}_1=\frac{d\Phi }{\mathrm{dy}},S{h}_2=-\frac{d\Psi }{\mathrm{dy}},\mathrm{at}y=1$$

(7)

The Nusselt number values were worked out from the process of numerical computation. The values of the Nusselt number are the same as those for the Sherwood numbers.

Numerical method

The numerical algorithm is based on the semi-implicit finite difference scheme, following Makinde et al. ¹⁵ and Chinyoka and Makinde. ¹⁶ The implicit terms are taken at the intermediate time level (N + ξ) where $0\le \xi \le 1$ . The algorithm employed in this article uses $\xi =1$ to allow the use of larger time steps, and it is therefore possible to work with any value of the time step. The discretization of the governing equations is based on a linear Cartesian mesh and uniform grid on which finite differences are taken. The PDEs governing the problem are approximated with second-order central differences for both the second and first spatial derivatives. The equations corresponding to the first and last grid points are modified to incorporate the boundary conditions.

The transformation for each component of the PDE where j represents position, N time, and $\Delta r$ mesh spacing, is done as follows

$$\frac{\partial \theta }{\partial t}=\frac{{\theta }_j^{N+1}-{\theta }_j^N}{\Delta t}$$

$$\frac{\partial \theta }{\partial r}=\frac{1}{2\Delta r}[({\theta }_{j+1}^{N+1}-{\theta }_{j-1}^{N+1})+({\theta }_{j+1}^N-{\theta }_{j-1}^N)]$$

and

$$\frac{{\partial }^2\theta }{\partial r^2}=\frac{1}{\Delta r^2}\left[(\xi )\left({\theta }_{j+1}^{N+1}-2{\theta }_j^{N+1}+{\theta }_{j-1}^{N+1}\right)+(1-\xi )\left({\theta }_{j+1}^N-2{\theta }_j^N+{\theta }_{j-1}^N\right)\right]$$

Similarly, we obtain expressions for

$$\frac{\partial \Phi }{\partial t};\frac{\partial \Psi }{\partial t};\frac{\partial \Phi }{\partial r};\frac{\partial \Psi }{\partial r};\frac{{\partial }^2\Phi }{\partial r^2};\frac{{\partial }^2\Psi }{\partial r^2}$$

It follows that the semi-implicit FDM for temperature, oxygen concentration, and carbon dioxide concentration are, respectively, expressed as follows

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\theta }^{N+1}-{\theta }^N}{\Delta t}=\frac{{\partial }^2}{\partial r^2}{\theta }^{N+\xi }+\frac{2}{r}\frac{\partial }{\partial r}{\theta }^N+\lambda {\left[{(1+\epsilon \theta )}^m{\Phi }^n\exp \left(\frac{\theta }{1+\epsilon \theta }\right)\right]}^N}\\ \multicolumn{1}{c}{\frac{{\Phi }^{N+1}-{\Phi }^N}{\Delta t}=\alpha [\frac{{\partial }^2}{\partial r^2}{\Phi }^{N+\xi }+\frac{2}{r}\frac{\partial }{\partial r}{\Phi }^N]-\lambda {\beta }_1{\left[{(1+\epsilon \theta )}^m{\Phi }^n\exp \left(\frac{\theta }{1+\epsilon \theta }\right)\right]}^N}\\ \multicolumn{1}{c}{\frac{{\Psi }^{N+1}-{\Psi }^N}{\Delta t}=\sigma [\frac{{\partial }^2}{\partial r^2}{\Psi }^{N+\xi }+\frac{2}{r}\frac{\partial }{\partial r}{\Psi }^N]+\lambda {\beta }_2{\left[{(1+\epsilon \theta )}^m{\Phi }^n\exp \left(\frac{\theta }{1+\epsilon \theta }\right)\right]}^N}\end{array}\}$$

(8)

The terms at the new time, $(N+1)$ , are placed on the left-hand side. These are the unknowns. The known values at the present time, $(N)$ , are moved to the right-hand side resulting in the following

$$\begin{array}{l}-\xi g{\theta }_{j+1}^{N+1}+\left(1+2\xi g\right){\theta }_j^{N+1}-\xi g{\theta }_{j-1}^{N+1}=-g(1-\xi ){\theta }_{j+1}^N+\left[1-2g(1-\xi )\right]{\theta }_j^N-g(1-\xi ){\theta }_{j-1}^N+\frac{1}{2r_j^2}g\left({\theta }_{j+1}^N-{\theta }_{j-1}^N\right)+\lambda \Delta t{\left[{(1+\epsilon \theta )}^m{\Phi }^{n}exp\left(\frac{\theta }{1+\epsilon \theta }\right)\right]}^N\\ -\xi g\alpha {\Phi }_{j+1}^{N+1}+\left(1+2\xi g\right)\alpha {\Phi }_j^{N+1}-\xi g\alpha {\Phi }_{j-1}^{N+1}=-g(1-\xi )\alpha {\Phi }_{j+1}^N+[1-2g(1-\xi )]\alpha {\Phi }_j^N-g(1-\xi )\alpha {\Phi }_{j-1}^N+\frac{\alpha }{2r_j^2}g\left({\Phi }_{j+1}^N-{\Phi }_{j-1}^N\right)-\lambda {\beta }_1\Delta t{\left[{(1+\epsilon \theta )}^m{\Phi }^{n}exp\left(\frac{\theta }{1+\epsilon \theta }\right)\right]}^N\\ -\xi g\sigma {\Psi }_{j+1}^{N+1}+(1+2\xi g)\sigma {\Psi }_j^{N+1}-\xi g\sigma {\Psi }_{j-1}^{N+1}=-g(1-\xi )\sigma {\Psi }_{j+1}^N+[1-2g(1-\xi )]\sigma {\Psi }_j^N-g(1-\xi )\sigma {\Psi }_{j-1}^N+\frac{\sigma }{2r_j^2}g\left({\Psi }_{j+1}^N-{\Psi }_{j-1}^N\right)+\lambda {\beta }_2\Delta t{\left[{(1+\epsilon \theta )}^m{\Phi }^{n}exp\left(\frac{\theta }{1+\epsilon \theta }\right)\right]}^N\end{array}\}$$

(9)

where $g=(\Delta t/\Delta r^2)$ . The system of equations resulting from equation (8) forms a tri-diagonal system, which can be solved using any numerical software such as MATLAB and MAPLE. Graphical solutions for ${\theta }^{N+\xi }$ , ${\Phi }^{N+\xi }$ , and ${\Psi }^{N+\xi }$ are given in the following section.

Results and discussion

Here, we present effects of various kinetic parameters on temperature, oxygen concentration, and carbon dioxide emission. Unless otherwise stated, the following parameters are used

$$m=0.5,n=1,\sigma =1,B{i}_1=1,B{i}_2=1,B{i}_3=1,{\beta }_1=1,{\beta }_2=1,\lambda =1,\alpha =1,\epsilon =0.1,t=10$$

Transient and steady flow profiles

In this case, temperature, carbon dioxide emission, and oxygen depletion profiles at various time intervals are provided. Figure 3 displays the transient increase in temperature with varying times until a steady state is attained. The same scenario is observed in Figure 4, where there is transient increase in carbon dioxide emission until steady state is reached. The general trend is that temperature and carbon dioxide emission profiles are highest at the center of the sphere and lowest at the surface thereof. This is due to continued exothermic chemical reaction that keeps on producing heat and carbon dioxide to increase high levels of these products inside the sphere where combustion takes place. At the surface of the sphere, heat and mass transfer occur, where heat is lost to the ambient and carbon dioxide diffuses out of the sphere to the ambient. A different observation is shown in Figure 5. In this case, we see a decrease in oxygen concentration as time increases, until a steady-state situation is attained. The reason for this behavior is that oxygen is consumed during combustion and its concentration is reduced but more concentration thereof is toward the surface of the sphere and diffuses into the sphere to enhance exothermic chemical reaction.

Figure 3.

Effect of increasing t on temperature.

Figure 4.

Effect of increasing t on carbon dioxide emission.

Figure 5.

Effect of increasing t on oxygen concentration.

Effects of thermo-physical parameters variation on sphere temperature profiles

Figures 6–13 show the effects of various parameters on temperature of the sphere. We observe that the temperature increases with corresponding increases in the parameters, n, $B{i}_1$ , and ${\beta }_1$ . The parameter n is the order of reaction and from the set of equations (5), we see that the index n relates to the oxygen (reactant) concentration. So increasing the n will reduce the rate of reaction according to the rate law expression. $B{i}_1$ (thermal Biot number) is related to the heat transfer resistances and hence is related to the material composition of the sphere. Its increase will increase heat transfer resistance and thus, the exothermic chemical reaction is reduced. The decrease in temperature of the sphere due to thermal Biot number is also attributed to the action of convective cooling at the surface of the sphere. The oxygen consumption rate parameter ${\beta }_1$ can be physically increased by reducing the oxygen supply to surrounding air, which will also discourage exothermic chemical reaction. The results are shown in Figures 6–8, respectively. The decrease in the profiles indicates a slowdown of the exothermic chemical reaction to reduce temperature increase and hence retention of the thermal stability of the system. A different scenario is observed in Figures 9–13, where an increase in $\alpha$ (carbon dioxide diffusivity), $\lambda$ (reaction rate), $\epsilon$ (activation energy), $B{i}_2$ (oxygen Biot number), and m (numerical index for kinetics type) gives a corresponding increase in temperature profiles. An increase in $B{i}_2$ implies an increase in the supply of oxygen from the surrounding environment to support the exothermic chemical reaction process, leading to more generation of internal heat in the system that corresponds to temperature raise of the sphere. Note that the temperature profile for biomolecular kinetics (m = 0.5) is the highest and the sensitized kinetics (m = −2) gives the lowest profile. This shows that thermal ignition occurs faster in bimolecular type of exothermic chemical reaction as compared to Arrhenius (m = 0) and sensitized kinetics. These parameters encourage exothermic chemical reaction to occur and if the raise in temperature is not well investigated, self-ignited fires may result.

Figure 6.

Effect of increasing n.

Figure 7.

Effect of increasing $B{i}_1$ .

Figure 8.

Effect of increasing ${\beta }_1$ .

Figure 9.

Effect of increasing $\alpha$ .

Figure 10.

Effect of increasing $\lambda$ .

Figure 11.

Effect of increasing $\epsilon$ .

Figure 12.

Effect of increasing $B{i}_2$ .

Figure 13.

Effect of increasing m.

Effects of thermo-physical parameters variation on sphere oxygen concentration

Effects of parameters variation on oxygen concentration are illustrated in Figures 14–21. Figures 14–17 show the effects of $n,B{i}_2,\alpha ,\mathrm{and}B{i}_1$ on oxygen concentration of the system. We see that an increase in each of the parameters results with an increase in oxygen concentration. These parameters reduce the exothermic chemical reaction, hence less oxygen depletion. In other words, these parameters reduce thermal ignition process of the system. Higher levels of oxygen concentration are attained during highest orders of exothermic chemical reactions compared to lowest orders, as indicated by the effect of n. The diffusivity of oxygen $(\alpha )$ into the system also facilitates more concentration of oxygen in the system. The increment of oxygen concentration with increasing values of thermal Biot number and oxygen Biot number is attributed to the combined effect of convective cooling and continuous oxygen supply from the surrounding environment at the surface of the sphere. We have a different observation given by Figures 18–21, where an increase in parameters ${\beta }_1,\lambda ,\epsilon ,\mathrm{and}m$ shows a decrease in oxygen concentration. The effect of oxygen consumption rate parameter ${\beta }_1$ reduces oxygen concentration by facilitation thermal ignition process. The rate of reaction $\lambda$ , the activation energy $\epsilon$ , and the exothermic chemical reaction type m also reduce the oxygen concentration by accelerating thermal ignition of the system. These parameters encourage exothermic chemical reaction of the system to occur. The more oxygen is consumed, the more the chemical reaction takes place with the production of heat and carbon dioxide.

Figure 14.

Effect of increasing n.

Figure 15.

Effect of increasing $B{i}_2$ .

Figure 16.

Effect of increasing $\alpha$ .

Figure 17.

Effect of increasing $B{i}_1$ .

Figure 18.

Effect of increasing ${\beta }_1$ .

Figure 19.

Effect of increasing $\lambda$ .

Figure 20.

Effect of increasing $\epsilon$ .

Figure 21.

Effect of increasing m.

Effects of thermo-physical parameters variation on sphere carbon dioxide emission

Here, we investigate the variation of parameters effects on carbon dioxide emitted by the sphere of combustible material as illustrated by Figures 22–32. We see from Figures 22–26 the effects of $n,B{i}_1,\sigma ,B{i}_3,\mathrm{and}{\beta }_1$ on carbon dioxide. An increase in each of the parameters stated result with a decrease in carbon dioxide emission. As discussed before, these parameters lesson thermal ignition which results with less carbon dioxide emission coupled with heat production. It is of environmental importance that less carbon dioxide which contributes so much to greenhouse effect is emitted. Figures 27–32 illustrate how the parameters ${\beta }_2,\epsilon ,\alpha ,B{i}_2,\lambda ,\mathrm{and}m$ influence the emission of carbon dioxide. An increase in the parameters gives a corresponding increase in the carbon dioxide emission. These parameters enhance the exothermic chemical reaction that results with using up more of oxygen to react with carbon containing material to produce heat and carbon dioxide that are not good to a healthy environment.

Figure 22.

Effect of increasing n.

Figure 23.

Effect of increasing $B{i}_1$ .

Figure 24.

Effect of increasing $\sigma$ .

Figure 25.

Effect of increasing $B{i}_3$ .

Figure 26.

Effect of increasing ${\beta }_1$ .

Figure 27.

Effect of increasing ${\beta }_2$ .

Figure 28.

Effect of increasing $\epsilon$ .

Figure 29.

Effect of increasing $\alpha$ .

Figure 30.

Effect of increasing $B{i}_2$ .

Figure 31.

Effect of increasing $\lambda$ .

Figure 32.

Effect of increasing m.

Effects of parameter variation on thermal criticality values or blowups

In this section, we present plots for thermal criticality values, Nusselt number Nu, versus the rate of reaction, Frank-Kamenetskii parameter $\lambda$ . The parameters $\alpha ,\epsilon ,B{i}_2,m,\mathrm{and}B{i}_1$ were varied to establish values that help to control exothermic chemical reactions or spontaneous ignition systems against any explosion. Figures 33–36 show that a combined decrease in the parameters $\alpha ,\epsilon ,B{i}_2,\mathrm{and}m$ enhances thermal stability. This is because the thermal criticality values ( $\lambda$ ) decrease with increasing values of these parameters as indicated in Table 1. Keeping $\alpha$ and $B{i}_2$ very low means that less oxygen from the ambient should be allowed to diffuse into the system to reduce thermal ignition in order to attain thermal stability. The same applies the activation energy $\epsilon$ and reaction kinetics index m. A different scenario is given by Figure 37 where an increase in $B{i}_1$ gives a corresponding increase in thermal criticality values, and this enhances thermal stability of the system. This is due to convective cooling at the surface of the sphere. Numerical value results are presented in Table 1.

Figure 33.

Effect of increasing $\alpha$ .

Figure 34.

Effect of increasing $\epsilon$ .

Figure 35.

Effect of increasing $B{i}_2$ .

Figure 36.

Effect of increasing m.

Table 1.

Computations showing the effects of various thermo-physical parameters on thermal criticality values.

$\alpha$	$\epsilon$	$B{i}_2$	m	$B{i}_1$	Nu	$\lambda$
1	0.1	1	0.5	1	0.88472	28.400
1.3	0.1	1	0.5	1	0.98499	7.3369
1.8	0.1	1	0.5	1	1.21023	4.3671
1	0.8	1	0.5	1	0.76113	7.8003
1	1	1	0.5	1	0.73653	6.4788
1	0.1	1.3	0.5	1	0.98156	8.4520
1	0.1	1.6	0.5	1	1.08804	5.8730
1	0.1	1	0	1	0.93270	85.0959
1	0.1	1	−2	1	0.99477	16817.8513
1	0.1	1	0.5	1.2	0.97089	487.2366
1	0.1	1	0.5	1.3	0.97150	536.8221

$n=1,\sigma =1,B{i}_3=1,{\beta }_1=1,{\beta }_2=1$ in all calculations.

Figure 37.

Effect of increasing $B{i}_1$ .

Conclusion

In this article, we investigated the effects of various parameters embedded within the system of a stockpile of combustible material modeled in a spherical domain. The study is based mainly on theory rather than experimentation. The significance of theoretical investigation of spontaneous combustion of combustible material in a stockpile is that it is cost-effective, easier, and faster than experimental approach to the investigation. The important thing is that much knowledge of kinetic parameters of reacting species is needed to do the theoretical study. In this investigation, we were able to identify kinetic parameters that influence the slowdown of exothermic reaction and so reducing the consumption of oxygen in the system. The slowing down of exothermic chemical reaction was indicated by the decrease in both the temperature and the carbon dioxide profiles, but with the increase in the oxygen profiles. In general, the following parameters, $n,B{i}_2,\alpha ,\mathrm{and}B{i}_1$ , showed the tendency to favor the retention of oxygen. On the contrary, parameters which enhance exothermic chemical reaction of the system were generally identified as ${\beta }_1,\lambda ,\epsilon ,\mathrm{and}m$ . It was also shown that in order to control exothermic chemical reaction in a spherical reactive stockpile from leading to explosions, the values of the following parameters, $\alpha ,\epsilon ,B{i}_2,\mathrm{and}m$ , should be kept low, but the value of $B{i}_1$ should be kept high.