Numerical Investigation into CO_𝟐 Emission, O_𝟐 Depletion, and Thermal Decomposition in a Reacting Slab

Abstract

The emission of carbon dioxide (CO₂) is closely associated with oxygen (O₂) depletion, and thermal decomposition in a reacting stockpile of combustible materials like fossil fuels (e.g., coal, oil, and natural gas). Moreover, it is understood that proper assessment of the emission levels provides a crucial reference point for other assessment tools like climate change indicators and mitigation strategies. In this paper, a nonlinear mathematical model for estimating the CO₂ emission, O₂ depletion, and thermal stability of a reacting slab is presented and tackled numerically using a semi-implicit finite-difference scheme. It is assumed that the slab surface is subjected to a symmetrical convective heat and mass exchange with the ambient. Both numerical and graphical results are presented and discussed quantitatively with respect to various parameters embedded in the problem.

1. Introduction

Studies relating to transient heating of combustible materials due to exothermic oxidation chemical reactions are extremely important and have a wide range of applications in industry, engineering, and environmental science [1]. For instance, fossil fuels (coal, oil, and natural gas) account for 85% of world’s primary energy supply, 70% of worlds electricity and heat generation and over 94% of energy for transportation [2]. The production and use of these combustible materials contribute up to 80% of CO₂ emission. Given expected increases in global population, economic growth, and energy demand, a continuous rise in emissions is expected unless fundamental technology changes occur in global energy systems which are currently dominated by fossil fuels. The CO₂ pollution is the principal human cause of global warming and climate change [3]. Meanwhile, for proper assessment of the CO₂ emission and O₂ depletion levels together with their impact on both the environment and life on Earth, knowledge of the mathematical models of these complex chemical systems is essential. These provide crucial reference points for other assessment tools like thermal stability of the materials, climate change indicators, and mitigation strategies. It may also help in developing medium to long-term action plans for climate change research and reliable design of the systems [4].

An extensive review of detailed chemical kinetic models for the heating-up of combustible materials is given by Simmie [5]. His review considered post-1994 work and focuses on the modeling of hydrocarbon fuel oxidation in the gas phase by detailed chemical kinetics and those experiments which validate them. Moreover, thermal combustion analysis has received much attention in the literature [6–8]. Several studies have been directed towards obtaining critical conditions for thermal ignition to occur in the form of a critical value for the Frank-Kamenetskii parameter [9]. Usually, chemical processes include many, up to a several hundred, intermediate elementary reactions [10]. For example, in combustion science, it is very common to use complex multistep reaction mechanisms to predict the oxidation of hydrocarbons [11]. However, the use of one-step decomposition kinetics clearly simplifies the complicated chemistry involved in the problem but is both practical and necessary without additional information about the individual decomposition reaction steps [12, 13]. Meanwhile, analytical solutions of the partial differential equations governing transient heating of the combustible material undergoing oxidation reactions are usually impossible or extremely difficult to obtain. The exothermic nature of such reactions leads to complex nonlinear transient interaction of heat conduction, mass diffusion, and chemical reactions, resulting in steep concentration and temperature gradients [14]. In such circumstances, a better understanding of the system behavior can only be accomplished by conducting numerical simulations to capture the frontal behavior of the processes.

The basic objective of this study is to provide a numerical estimate for the thermal stability together with the rate of CO₂ emission and O₂ depletion in transient heating of a slab of combustible material in the presence of convective heat and mass exchange with the ambient at the slab surface. The mathematical formulation of the problem is established in section two. In section three, the semi-implicit finite difference technique is implemented to tackle the problem. Both numerical and graphical results are presented and discussed quantitatively with respect to various parameters embedded in the system in Section 4.

2. Mathematical Model

We consider a stockpile of combustible material in a rectangular slab. It is assumed that the slab is undergoing an 𝑛th-order oxidation chemical reaction, and its surface is subjected to a symmetrical convective heat and mass exchange with the ambient (see Figure 1). The complicated chemistry involved in this problem can be simplified by assuming a one-step finite-rate irreversible reaction between the combustible material (hydrocarbon) and the oxygen of the air; that is,C_𝑖H_𝑗+𝑗𝑖+4O₂⟶𝑖CO₂+𝑗2H₂O+Heat. (2.1) Following [1, 5, 6, 12, 14], the nonlinear partial differential equations describing temperature, oxygen, and carbon dioxide concentration in the combustible material can be written as𝜌𝑐_𝑝𝜕𝑇𝜕𝑡𝜕=𝑘²𝑇𝜕𝑦²+𝑄𝐴𝐾𝑇𝜈𝑙^𝑚𝐶^𝑛−𝐸exp,𝑅𝑇𝜕𝐶𝜕𝑡𝜕=𝐷²𝐶𝜕𝑦²−𝐴𝐾𝑇𝜈𝑙^𝑚𝐶^𝑛−𝐸exp,𝑅𝑇𝜕𝑃𝜕𝑡𝜕=𝛾²𝑃𝜕𝑦²+𝐴𝐾𝑇𝜈𝑙^𝑚𝐶^𝑛−𝐸exp,𝑅𝑇 (2.2) with initial and boundary conditions as follows:𝑇𝑦,0=𝑇₀,𝐶=1𝑦,02𝐶_𝑎,𝑃𝑦,0=0,𝜕𝑇𝜕𝑦0,𝑡=𝜕𝐶𝜕𝑦0,𝑡=𝜕𝑃𝜕𝑦0,𝑡=0,−𝑘𝜕𝑇𝜕𝑦𝑎,𝑡=ℎ₁𝑇𝑎,𝑡−𝑇_𝑎,−𝐷𝜕𝐶𝜕𝑦𝑎,𝑡=ℎ₂𝐶𝑎,𝑡−𝐶_𝑎,−𝛾𝜕𝑃𝜕𝑦𝑎,𝑡=ℎ₃𝑃𝑎,𝑡−𝑃_𝑎, (2.3) where 𝑇 is the absolute temperature, 𝐶 is the (depletion) oxygen concentration, 𝑃 is the carbon dioxide emission concentration, 𝑇_𝑎 is the ambient temperature, 𝐶_𝑎 is the oxygen concentration in the surrounding air, 𝑃_𝑎 is the carbon dioxide concentration in the surrounding air, 𝑡 is the time, 𝑇₀ is the slab initial temperature, the initial depletion of oxygen in the slab is zero, 𝜌 is the density, 𝑐_𝑝 specific heat at constant pressure, 𝑘 is the thermal conductivity of the reacting slab, 𝐷 is the diffusivity of oxygen in the slab, 𝛾 is the diffusivity of carbon dioxide in the slab, 𝑄 is the exothermicity, 𝐴 is the rate constant, 𝐸 is the activation energy, 𝑅 is the universal gas constant, 𝑙 is the Planck number, 𝐾 is the Boltzmann constant, 𝜈 is the vibration frequency, 𝑎 is the slab half width, 𝑦 is the distance measured transverse direction, ℎ₁ is the coefficient of heat transfer between the slab and its surroundings, ℎ₂ is the coefficient of oxygen transfer between the slab and its surroundings, ℎ₃ is the coefficient of carbon dioxide transfer between the slab and its surroundings, 𝑛 is the order of exothermic chemical reaction, and 𝑚 is the numerical exponent such that 𝑚∈{−2,0,0.5}. The three values taken by the parameter 𝑚 represent the numerical exponent for sensitized, Arrhenius and Bimolecular kinetics, respectively, (see [1, 6]). We introduce the following dimensionless variables into (2.2)-(2.3):𝐸𝜃=𝑇−𝑇₀𝑅𝑇²₀,𝜃_𝑎=𝐸𝑇_𝑎−𝑇₀𝑅𝑇²₀𝐶,Φ=𝐶_𝑎𝑃,Ψ=𝑃_𝑎,Bi₁=𝑎ℎ₁𝑘,Bi₂=𝑎ℎ₂𝐷,Bi₃=𝑎ℎ₂𝛾,𝛽₁=𝜌𝑐_𝑝𝑅𝑇²₀𝑄𝐸𝐶_𝑎,𝛽₂=𝜌𝑐_𝑝𝑅𝑇²₀𝑄𝐸𝑃_𝑎,𝜆=𝐾𝑇₀𝜈𝑙^𝑚𝑄𝐴𝐸𝑎²𝐶^𝑛_𝑎𝑘𝑅𝑇²₀−𝐸exp,𝑅𝑇𝑦=𝑦𝑎𝑘,𝑡=𝑡𝜌𝑐_𝑝𝑎²,𝜀=𝑅𝑇₀𝐸,𝛼=𝐷𝜌𝑐_𝑝𝑘,𝜎=𝛾𝜌𝑐_𝑝𝑘, (2.4) and we obtain the dimensionless governing equations as𝜕𝜃=𝜕𝜕𝑡²𝜃𝜕𝑦²+𝜆(1+𝜀𝜃)^𝑚Φ^𝑛𝜃exp,1+𝜀𝜃𝜕Φ𝜕𝜕𝑡=𝛼²Φ𝜕𝑦²−𝜆𝛽₁(1+𝜀𝜃)^𝑚Φ^𝑛𝜃exp,1+𝜀𝜃𝜕Ψ𝜕𝜕𝑡=𝜎²Ψ𝜕𝑦²+𝜆𝛽₂(1+𝜀𝜃)^𝑚Φ^𝑛𝜃exp.1+𝜀𝜃 (2.5) The corresponding initial and boundary conditions then become𝜃𝑦,0=0,Φ𝑦,0=0.5,Ψ𝑦,0=0,𝜕𝜃()=𝜕𝑦0,𝑡𝜕Φ()=𝜕𝑦0,𝑡𝜕Ψ()𝜕𝑦0,𝑡=0,𝜕𝜃()𝜕𝑦1,𝑡=−Bi₁𝜃(1,𝑡)−𝜃_𝑎,𝜕Φ()𝜕𝑦1,𝑡=−Bi₂[Φ(1,𝑡)−1],𝜕Ψ()𝜕𝑦1,𝑡=−Bi₃[Ψ(1,𝑡)−1], (2.6) where 𝜆,𝜀,𝛽₁,𝛽₂,𝛼,𝜎,Bi₁,Bi₂,andBi₃ represent the Frank-Kamenetskii parameter, activation energy parameter, oxygen consumption rate parameter, carbon dioxide emission rate parameter, oxygen diffusivity parameter, carbon dioxide diffusivity parameter, the thermal Biot number, oxygen Biot number, and carbon dioxide Biot number, respectively. A body of material releasing heat to its surroundings may achieve a safe steady state where the temperature of the body reaches some moderate value and stabilizes. However, when the rate of heat generation in the material exceeds the rate of heat loss to the surroundings, then ignition can occur. In the following sections, (2.5)-(2.6) are solved numerically using a semi-implicit finite difference scheme.

Figure 1 

Geometry of the problem.

3. Numerical Solution

Our numerical algorithm is based on the semi-implicit finite difference scheme [15–18]. The implicit terms are taken at the intermediate time level (𝑁+𝜉) where 0≤𝜉≤1. The algorithm employed in [18] uses 𝜉=1/2, we will, however, follow the formulation in [15–17] and take 𝜉=1 in this paper so that we can use larger time steps. In fact being nearly fully implicit, our numerical algorithm presented in this paper is conjectured to work for 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. We approximate both the second and first spatial derivatives with second-order central differences. The equations corresponding to the first and last grid points are modified to incorporate the boundary conditions. The semi-implicit schemes for the temperature, O₂ concentration, and CO₂ concentration, respectively, read𝜃^(𝑁+1)−𝜃^(𝑁)=𝜕Δ𝑡𝜕𝑦²²𝜃^(𝑁+𝜉)+𝜆(1+𝜀𝜃)^𝑚Φ^𝑛𝜃exp1+𝜀𝜃^(𝑁),Φ^(𝑁+1)−Φ^(𝑁)𝜕Δ𝑡=𝛼²𝜕𝑦²Φ^(𝑁+𝜉)−𝜆𝛽₁(1+𝜀𝜃)^𝑚Φ^𝑛𝜃exp1+𝜀𝜃^(𝑁),Ψ^(𝑁+1)−Ψ^(𝑁)𝜕Δ𝑡=𝜎²𝜕𝑦²Ψ^(𝑁+𝜉)+𝜆𝛽₂(1+𝜀𝜃)^𝑚Φ^𝑛𝜃exp1+𝜀𝜃^(𝑁). (3.1) The equations for 𝜃^(𝑁+1), Φ^(𝑁+1), and Ψ^(𝑁+1), thus, become−𝑟𝜉𝜃^(𝑁+1)_𝑗−1+(1+2𝑟𝜉)𝜃_𝑗^(𝑁+1)−𝑟𝜉𝜃^(𝑁+1)_𝑗+1=−𝑟(1−𝜉)𝜃^(𝑁)_𝑗−1+(1−2𝑟(1−𝜉))𝜃_𝑗^(𝑁)−𝑟(1−𝜉)𝜃^(𝑁)_𝑗+1+𝜆Δ𝑡1+𝜀𝜃_𝑗^𝑚Φ^𝑛_𝑗𝜃exp^𝑗1+𝜀𝜃_𝑗^(𝑁),−𝑟𝜉𝛼Φ^(𝑁+1)_𝑗−1+(1+2𝑟𝜉𝛼)Φ_𝑗^(𝑁+1)−𝑟𝜉𝛼Φ^(𝑁+1)_𝑗+1=−𝑟(1−𝜉)𝛼Φ^(𝑁)_𝑗−1+(1−2𝑟(1−𝜉)𝛼)Φ_𝑗^(𝑁)−𝑟(1−𝜉)𝛼Φ^(𝑁)_𝑗+1−𝜆Δ𝑡𝛽₁1+𝜀𝜃_𝑗^𝑚Φ^𝑛_𝑗𝜃exp^𝑗1+𝜀𝜃_𝑗^(𝑁),−𝑟𝜉𝜎Ψ^(𝑁+1)_𝑗−1+(1+2𝑟𝜉𝜎)Ψ_𝑗^(𝑁+1)−𝑟𝜉𝜎Ψ^(𝑁+1)_𝑗+1=−𝑟(1−𝜉)𝜎Ψ^(𝑁)_𝑗−1+(1−2𝑟(1−𝜉)𝜎)Ψ_𝑗^(𝑁)−𝑟(1−𝜉)𝜎Ψ^(𝑁)_𝑗+1+𝜆Δ𝑡𝛽₂1+𝜀𝜃_𝑗^𝑚Φ^𝑛_𝑗𝜃exp^𝑗1+𝜀𝜃_𝑗^(𝑁), (3.2) where 𝑟=Δ𝑡/Δ𝑦². The solution procedures reduce to inversion of tridiagonal matrices. The schemes (3.2) were checked for consistency. For 𝜉=1, these are first-order accurate in time but second order in space. The schemes in [18] have 𝜉=1/2 which improves the accuracy in time to second order. We use 𝜉=1 here so that we are free to choose larger time steps and still converge to the steady solutions. As already conjectured, our algorithm works for any value of the time step! The code was checked for both spatial and temporal convergence. In particular, solutions calculated from, say Δ𝑡=1, using 200 time steps are exactly the same as those after 40 time steps with Δ𝑡=5 or those after 20,000 time steps with Δ𝑡=0.01. Similarly solutions using Δ𝑦=0.02 converge to the same results as those say for Δ𝑡=0.025 or Δ𝑡=0.05, and so forth. Our code, thus, runs extremely fast, and; hence, we can easily obtain and, thus, present all our results at steady state using nearly insignificant computational times.

4. Results and Discussion

Unless otherwise stated, we employ the parameter values:𝑚=0.5,𝑛=1,𝜃_𝑎𝛽=0.1,𝛼=1,𝜎=1,₁=1,𝛽₂=1,Bi₁=1,Bi₂=1,Bi₃=1,𝜀=0.1,Δ𝑦=0.01,Δ𝑡=1,𝑡=200. (4.1) These will be the default values in this work, and; hence, in any graph where any of these parameters is not explicitly mentioned, it will be understood that such parameters take on the default values.

4.1. Transient and Steady Flow Profiles

We display the transient solutions in Figures 2, 3, and 4. At the given parameter values, Figure 2 shows a transient increase in temperature until a steady-state is reached. A similar scenario obtains in Figure 4 in which a transient increase of carbon dioxide emission is observed until a steady state is reached. An opposite situation is noticed in Figure 3 where a decrease in oxygen concentration is observed with increasing time until a steady-state concentration is attained. These results are consistent with intuition regarding exothermic oxidation reactions.

Figure 2 

Transient and steady-state temperature profiles.

Figure 3 

Transient and steady state oxygen profiles.

Figure 4 

Transient and steady state carbon dioxide profiles.

4.1.1. Parameter Dependance of Solutions in Steady State

It is understood from Figures 2–4 that, at the default parameter values, solutions have reached steady state at, say, times 𝑡≥40. All the solutions at 𝑡=200 given below will, thus, be understood to be steady solutions. The dependence of solutions on parameter 𝑛 is illustrated in Figures 5, 6, and 7. As 𝑛 increases, the results show a decrease in both the temperature and the CO₂ emission. The reduced oxidation reactions mean lower oxygen consumption, hence, lead to a corresponding increase in O₂ concentration.

Figure 5 

Effects of 𝑛 on temperature.

Figure 6 

Effects of 𝑛 on O2.

Figure 7 

Effects of 𝑛 on CO2.

The dependence of solutions on parameter 𝑚 is illustrated in Figures 8, 9, and 10. As 𝑚 increases, the results show an increase in both the temperature and the CO₂ emission. The increased oxygen consumption from higher oxidation reactions correspondingly decreases O₂ concentration.

Figure 8 

Effects of 𝑚 on temperature.

Figure 9 

Effects of 𝑚 on O2.

Figure 10 

Effects of 𝑚 on CO2.

The dependence of solutions on parameter Bi₁ is illustrated in Figures 11, 12, and 13. As Bi₁ increases, the results show a decrease in both the temperature and the CO₂ emission. The reduced oxidation reactions mean lower oxygen consumption, hence, lead to a corresponding increase in O₂ concentration.

Figure 11 

Effects of Bi1 on temperature.

Figure 12 

Effects of Bi1 on O2.

Figure 13 

Effects of Bi1 on CO2.

The dependence of solutions on parameter Bi₂ is illustrated in Figures 14, 15, and 16. An increase in Bi₂ directly corresponds to an increase in O₂ concentration. This correspondingly increases the oxidation reactions and, hence, increases the temperature and, hence, the CO₂ emission.

Figure 14 

Effects of Bi2 on temperature.

Figure 15 

Effects of Bi2 on O2.

Figure 16 

Effects of Bi2 on CO2.

The dependence of CO₂ emission on Bi₃ is shown in Figure 17. An increase in Bi₃ decreases the CO₂ emission. No noticeable effects are observed in the temperature and the O₂ concentration.

Figure 17 

Effects of Bi3 on CO2.

The dependence of solutions on the reaction parameter 𝜆 is illustrated in Figures 18, 19, and 20. An increase in 𝜆 directly corresponds to an increase in the (exothermic) oxidation reactions and, hence, increases the temperature and, hence, the CO₂ emission. The increased oxidation reactions increase oxygen consumption and thus decrease O₂ concentration.

Figure 18 

Effects of 𝜆 on temperature.

Figure 19 

Effects of 𝜆 on O2.

Figure 20 

Effects of 𝜆 on CO2.

The dependence of solutions on parameter 𝛽₁ is shown in Figures 21, 22, and 23. An increase in 𝛽₁ directly corresponds to a decrease in O₂ concentration. This correspondingly decreases the oxidation reactions and, hence, decreases the temperature and also CO₂ emission.

Figure 21 

Effects of 𝛽1 on temperature.

Figure 22 

Effects of 𝛽1 on O2.

Figure 23 

Effects of 𝛽1 on CO2.

The dependence of CO₂ emission on 𝛽₂ is shown in Figure 24. An increase in 𝛽₂ increases the CO₂ emission. No noticeable effects are observed in the temperature and the O₂ concentration.

Figure 24 

Effects of 𝛽2 on CO2.

The dependence of solutions on the parameter 𝜀 is illustrated in Figures 25, 26, and 27. An increase in 𝜀 decreases the temperature and the CO₂ emission. The reduced oxygen consumption leads to an increase in O₂ concentration.

Figure 25 

Effects of 𝜀 on temperature.

Figure 26 

Effects of 𝜀 on O2.

Figure 27 

Effects of 𝜀 on CO2.

The dependence of solutions on the parameter 𝛼 is shown in Figures 28, 29, and 30. An increase in 𝛼 directly corresponds to an increase in O₂ concentration. This correspondingly increases the oxidation reactions and hence increases the temperature and also CO₂ emission.

Figure 28 

Effects of 𝛼 on temperature.

Figure 29 

Effects of 𝛼 on O2.

Figure 30 

Effects of 𝛼 on CO2.

The dependence of CO₂ emission on 𝜎 is shown in Figure 31. An increase in 𝜎 decreases the CO₂ emission. No noticeable effects are observed in the temperature and O₂ concentration.

Figure 31 

Effects of 𝜎 on CO2.

4.2. Thermal Stability and Blowup

In Figures 32, 33, 34, and 35, we plot 𝜃_𝑦(1,𝑡) against 𝜆 for varying values of 𝑛, 𝑚, Bi₁, and Bi₂, respectively. The solutions are given up to the values of 𝜆 at which the onset of blowup in temperature is observed. We notice that parameters that increase the temperature (𝜃) correspondingly increase 𝜃_𝑦(1,𝑡). The results also show that early onset of blowup (and, hence, thermal instability) can be delayed by using higher values of 𝑛, lower values of 𝑚, higher values of Bi₁, lower values of Bi₂, and so forth.

Figure 32 

Effects of variation of 𝑛 on the thermal criticality or blowup in the system.

Figure 33 

Effects of variation of 𝑚 on the thermal criticality or blowup in the system.

Figure 34 

Effects of Bi1 variation on the thermal criticality or blowup in the system.

Figure 35 

Effects of Bi2 variation on the thermal criticality or blowup in the system.

5. Conclusion

We develop an unconditionally stable (works for any time step size) and convergent semi-implicit finite-difference scheme and utilize it to computationally investigate the transient dynamics of CO₂ emission, O₂ depletion, and thermal decomposition in a reacting slab. The solutions show that those processes that increase the oxidation reactions lead to enhanced oxygen depletion as well as increased carbon dioxide emission. The results also show that enhanced thermal stability and reduced carbon dioxide emission are best achieved by cutting down on those processes that would otherwise increase the oxidation reactions.