# Computational Modelling of Thermal Stability in a Reactive Slab with Reactant Consumption

Authors: O. D. Makinde
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

This paper investigates both the transient and the steady state of a one-stepnth-order oxidation exothermic reaction in a slab of combustible material with an insulated lower surface and an isothermal upper surface, taking into consideration reactant consumption. The nonlinear partial differential equation governing the transient reaction-diffusion problem is solved numerically using a semidiscretization finite difference technique. The steady-state problem is solved using a perturbation technique together with a special type of the Hermite-Padé approximants. Graphical results are presented and discussed quantitatively with respect to various embedded parameters controlling the systems. The crucial roles played by the boundary conditions in determining the thermal ignition criticality are demonstrated.


---

## Body

## 1.Introduction

Analysis of possible development of runaway at production, storage, and use of a chemical product and subsequent choice of measures that can prevent an accident or mitigate its consequences are the important tasks of reaction hazards assessment [1]. The kinetic model of a reaction that describes heat generation plays a crucial role in the complete explosion model. Although thermal explosion has received much attention in the literature, the vast majority of investigations have been concerned with homogeneous boundary conditions ranging from the infinite Biot number case of a constant surface temperature [2] through a range of Biot numbers to zero [3]. In most of these studies, the determination of critical conditions that separate explosive and nonexplosive domains of a proceeding reaction and evaluation of induction period of an explosion if it appears has been the main focus. Moreover, two main approaches are used for obtaining the necessary data for critical condition. The first approach is based on direct determination of the explosion characteristics by means of explosive experiments. This approach is very expensive, dangerous, and time consuming [4]. The second approach involves the application of mathematical theory of thermal combustion. It involves derivation of appropriate mathematical models that allow major chemical and physical processes within an exothermic reacting system to be taken into account and application of sophisticated analytical and numerical techniques to tackle the problem. Several authors have analysed theoretically the problem of thermal explosion in a reacting slab; they include Zaturska and Banks [5], Bebernes and Eberly [6] and Makinde [7]. In all these earlier studies, the combined effect of asymmetric boundary condition and reactant consumption on the thermal ignition criticality has not been properly reported.In the present study, the theoretical analysis of Makinde [7] is extended to include the effect of asymmetric boundary condition on both transient and steady-state exothermic nth-order oxidation reaction in a slab of combustible material with reactant consumption. This paper is organized as follows; firstly, the governing partial differential equations for oxidation reactions are presented and solved numerically using the semidiscretization finite difference technique known as method of lines. Secondly, the steady-state problem is tackled using perturbation technique coupled with a special type of the Hermite-Padé approximant in order to obtain the thermal criticality conditions in the system. Pertinent results are presented graphically and discussed quantitatively.


## 2.Mathematical Model

We consider the transient problem of combustible material undergoing annth-order oxidation chemical reaction in a slab with insulated lower surface (see Figure 1). The complicated chemistry involved in this problem may be simplified by assuming a one-step finite-rate irreversible reaction given by
(2.1)CiHj+(i+j4)O2⟶iCO2+j2H2O+HeatFigure 1
Geometry of the problem.[Combustible material + Oxygen→ Heat + Carbon dioxide + Water].The dimensionless equations that describe the physical situation are given by [2, 3, 5–7]
(2.2)∂θ∂t=∂2θ∂y2+λ(1+εθ)mϕneθ/(1+εθ),(2.3)α∂ϕ∂t=∂2ϕ∂y2-λβ(1+εθ)mϕneθ/(1+εθ),
with initial and boundary conditions as
(2.4)θ(y,0)=0,ϕ(y,0)=0,(2.5)θ(1,t)=θw,ϕ(1,t)=1,(2.6)∂θ∂y(0,t)=0,ϕ(0,t)=1,
where λ is the Frank-Kamenetskii parameter, ε is the activation energy parameter, β is oxygen consumption rate parameter, α oxygen diffusivity parameter, n is the order of exothermic chemical reaction, and m is the numerical exponent given such that m={-2,0,1/2} represent numerical exponent for sensitised, Arrhenius, and bimolecular kinetics, respectively (see [6]). Equations (2.2)–(2.6) are obtained after introducing the dimensionless variables and quantities into the governing energy balance and concentration equations; that is,
(2.7)ϕ=(C-C0)(Cw-C0),θ=E(T-T0)RT02,θw=E(Tw-T0)RT02,y=y-a,t=kt-ρcpa2,ε=RT0E,α=kDρcp,β=kRT02DQE(Cw-C0),λ=QAEa2(Cw-C0)nkRT02[KT0υl]me-(E/(RT0)),
where T is the absolute temperature, t- is the time,Tw is the slab upper surface temperature, Cw is the slab surface oxygen concentration, T0 is the slab initial temperature, C0 is the initial oxygen concentration in the material, θ is the dimensionless temperature, θw is the slab upper surface dimensionless temperature, ϕ is the dimensionless oxygen concentration, (x,y) represents the Cartesian coordinates, ρ is the density, cp specific heat at constant pressure, k is the thermal conductivity of the material, 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, K is the Boltzmann constant, υ is the vibration frequency, a is the slab width, y- is the distance measured transverse direction, and D is the diffusivity of oxygen in the material.


## 3.Numerical Procedure

Here we employed the method of lines as our solution technique [8]. The governing equations (2.2)-(2.3) with the initial and boundary conditions (2.4)–(2.6) are transformed into a system of ODEs using finite differences for the spatial derivatives. Let Δy=1/N, yi=(i-1)Δy, 1≤i≤N+1; and θi(t), ϕi(t) represents θ(yi,t), ϕ(yi,t) respectively, then the semidiscrete system for the problem reads
(3.1)dθidt=1(Δy)2(θi+1-2θi+θi-1)+λ(1+εθi)mϕin[eθi/(1+εθi)],(3.2)αdϕidt=1(Δy)2(ϕi+1-2ϕi+ϕi-1)-λβ(1+εθi)mϕin[eθi/(1+εθi)],
with initial conditions
(3.3)θi(0)=0,ϕi(0)=0,1≤i≤N+1.
The first and last grid points are modified to incorporate the boundary conditions; that is,
(3.4)θ1=θ2,θN+1=θw,(3.5)ϕ1=1,ϕN+1=1.
The MAPLE program is employed to solve (3.1)–(3.5) using a fourth-order Runge-Kutta method.


## 4.Steady-State Analysis

A body of chemically reacting material releasing heat to its surroundings may achieve a safe steady state where the temperature of the body reaches some moderate value and stabilizes. Once a steady state is attained, (2.2)–(2.6) then become
(4.1)d2θdy2+λ(1+εθ)mϕneθ/(1+εθ)=0,(4.2)d2ϕdy2-λβ(1+εθ)mϕneθ/(1+εθ)=0,
with
(4.3)θ(1)=θw,ϕ(1)=1,(4.4)dθdy(0)=0,ϕ(0)=1.
The nonlinear nature of (4.1)–(4.4) precludes its exact solution. However, it is convenient to form a power series expansion in the Frank-Kamenetskii parameter; λ, that is,
(4.5)θ=∑i=0∞θiλi,ϕ=∑i=0∞ϕiλi.
Substituting the solution series (4.5) into (4.1)–(4.4) and collecting the coefficients of like powers of λ, we obtained and solved the equations for the coefficients of solution series iteratively. The solutions for the temperature and the oxygen concentration in the slab are given as
(4.6)θ(y)=θw-λ2(1+εθw)meθw/(1+εθw)(y2-1)+O(λ2),θ(y)=1+λβ2(1+εθw)meθw/(1+εθw)y(y-1)+O(λ2).
Using MAPLE, we obtained the first few terms of the above solution series. It is well known that this power series solution is valid for very small parameter values. However, by using the Hermite-Padé approximation technique [9], the usability of the solution series is extended beyond small parameter values as illustrated in the following section.


## 5.Thermal Criticality Determination

When the rate of heat generation in the reacting slab exceeds the rate of heat loss to the surroundings, then ignition can occur. Hence, the evaluation of critical regimes that separate the regions of explosive and nonexplosive ways of chemical reactions is extremely important from the application point of view. In order to achieve this goal, we employ a simple technique of series summation and improvement based on the generalization of Padé approximation technique (Baker and Graves-Morris [10]) and may be described as follows. Let
(5.1)UN(λ)=∑i=0Naiλi+O(λN+1),asλ⟶0,
be a given partial sum. It is important to note here that (5.1) can be used to approximate any output of the solution of the problem under investigation (e.g., the series for the wall heat flux parameter in terms of the Nusselt number Nu=-dθ/dy at y=1), since everything can be Taylor expanded in the given small parameter. Assume the U(λ) is a local representation of an algebraic function of λ in the context of nonlinear problems; we construct a multivariate series expression of the form
(5.2)Fd(λ,U)=∑m=1d∑j=0mfm-j,jλm-jUj,
of degree d≥2, such that
(5.3)∂Fd∂U(0,0)=1,Fd(λ,UN)=O(λN+1),asλ⟶0.
The requirement (5.3) yields f0,1=1 ensuring that the polynomial Fd has only one root which vanishes at λ=0 and reduces the problem to a system of N linear equations for the unknown coefficients of Fd. The entries of the underlying matrix depend only on the N given coefficients ai  in (5.1); consequently, we take N=(d2+3d-2)/2, so that the number of equations equals the number of unknowns. The polynomial Fd is a special type of the Hermite-Padé approximant [7, 9, 11] and is then investigated for bifurcation and criticality conditions using the Newton diagram [12]. The chief merit of this method is its ability to reveal the solution branches, criticality values as well as extending the usability of the power series solution beyond small parameter values.


## 6.Results and Discussion

Computational results in Table1 illustrate the rapid convergence of Hermite-Padé approximation procedure highlighted in the above section with gradual increase in the number of series coefficients utilized for the approximants. In Table 2, we observed that the magnitude of thermal ignition criticality (λc) increases with an increase in the parameter values of, β,n,ε and a decrease in the slab upper surface temperature parameter θw. Consequently, a delay in the development of thermal runaway in the reacting slab will be experienced, hence, enhances thermal stability of the system. Moreover, it is noteworthy from Table 2 that thermal ignition occur faster in a bimolecular (m=0.5) type of exothermic oxidation reaction as compared to the Arrhenius (m=0) and sensitised (m=-2) type of reaction.Table 1
Computations showing the criticality procedure rapid convergence (β=n=1,θw=ε=0).
d
N
N
u
=
-
θ
′
(
1
)
λ
c
N
2
4
2.293830707
1.07039427340
3
8
2.289550271
1.07040213355
4
13
2.289549038
1.07040198565
5
19
2.289549037
1.07040198791
6
26
2.289549037
1.07040198791Table 2
Computations showing thermal ignition criticality for different parameter values.
β
n
θ
w
m
ε
N
u
=
-
θ
′
(
1
)
λ
c
0.1
1
0.0
0.0
0.00
2.025223122
0.8941780502
0.5
1
0.0
0.0
0.00
2.134510046
0.9640287456
1.0
1
0.0
0.0
0.00
2.289549037
1.0704019879
0.1
3
0.0
0.0
0.00
2.074810972
0.9265325131
0.1
5
0.0
0.0
0.00
2.122974429
0.9600824416
0.1
1
0.5
0.0
0.00
2.025223122
0.5423464026
0.1
1
1.0
0.0
0.00
2.025223122
0.3289497214
0.1
1
0.0
0.5
0.01
2.057311239
0.8993850447
0.1
1
0.0
−2.0
0.01
2.112233348
0.9231510103
0.1
1
0.0
0.0
0.01
2.068063491
0.9040392715
0.1
1
0.0
0.0
0.10
2.600639145
1.0110740529A slice of the bifurcation diagram for0≤ε≪1 in the (λ,Nu) plane is shown in Figure 2. It represents the qualitative change in the thermal system as parameter λ increases. In particular, for 0≤ε≪1,β>0, and n>0, there is a critical value λc (a turning point) such that, for 0<λ<λc, there are two solutions (labeled I and II). The upper and lower solution branches occur due to nonlinearity in model equations for energy and concentration balance. When λ>λc the system has no real solution and displays a classical form indicating thermal runaway. As exothermic reaction due to oxidation chemical kinetics increases, the slab temperature increases uncontrollably until it ignites.Figure 2
A slice of approximate bifurcation diagram in the (λ, Nu(β=n=1, θw=ε=0)) plane.


### 6.1.Effect of Various Parameters on Temperature and Oxygen Concentration Profiles

The effects of various thermophysical parameters on the slab temperature and oxygen concentration profiles are displayed in Figures3, 4, 5, 6, 7, 8, and 9. Generally, the temperature is maximum at the slab lower insulated surface and decreases transversely with minimum value at its upper surface. Meanwhile, the oxygen concentration is lowest along the slab centerline (y=0.5) and maximum at the slab surfaces. This can be attributed to the fact that oxygen is utilized within the slab during exothermic chemical kinetics and fresh supply of oxygen from the surrounding is obtained at the slab surfaces. The evolution of the temperature and oxygen concentration in the slab is illustrated in Figures 3 and 4. It is noteworthy that the slab temperature increases while the oxygen concentration decreases gradually with time until it attains its steady-state value. Once the steady-state value is attained, the slab temperature and oxygen concentration remain the same for a given set of parameter values with respect to a further increase in time. In Figure 5, we observed that the slab temperature is highest during bimolecular reaction (m=0.5) and lowest for sensitized reaction (m=-2), hence confirming the earlier results in Table 2. Consequently, oxygen concentration in the slab is lowest during bimolecular reaction and highest for sensitized reaction. In Figures 6 and 7, we observed that the slab temperature decreases while oxygen concentration increases with an increase in the reaction order index (n) and activation energy parameter (ε). This clearly implies that a higher-order exothermic oxidation chemical reaction will be more thermally stable than a lower one. Figures 8 and 9 illustrated the effect of the Frank-Kamenetskii parameter (λ) and the slab upper surface temperature parameter (θw) on the slab temperature and oxygen concentration. The slab temperature increases while the oxygen concentration decreases with an increase in the parameter values of λ and θw. As λ and θw increase, the oxygen consumption within the slab increases and the slab internal heat generation due to exothermic oxidation reaction increases, this invariably leads to an elevation in the slab temperature.Temperature and oxygen concentration profiles form=0.5, n=2, α=β=λ=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5, n=2, α=β=λ=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles forn=β=λ=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,λ=β=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,n=λ=β=1,θw=0.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,n=β=1,θw=0,ε= 0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,n=λ=β=1,ε=0.1.
(a)
(b)


## 6.1.Effect of Various Parameters on Temperature and Oxygen Concentration Profiles

The effects of various thermophysical parameters on the slab temperature and oxygen concentration profiles are displayed in Figures3, 4, 5, 6, 7, 8, and 9. Generally, the temperature is maximum at the slab lower insulated surface and decreases transversely with minimum value at its upper surface. Meanwhile, the oxygen concentration is lowest along the slab centerline (y=0.5) and maximum at the slab surfaces. This can be attributed to the fact that oxygen is utilized within the slab during exothermic chemical kinetics and fresh supply of oxygen from the surrounding is obtained at the slab surfaces. The evolution of the temperature and oxygen concentration in the slab is illustrated in Figures 3 and 4. It is noteworthy that the slab temperature increases while the oxygen concentration decreases gradually with time until it attains its steady-state value. Once the steady-state value is attained, the slab temperature and oxygen concentration remain the same for a given set of parameter values with respect to a further increase in time. In Figure 5, we observed that the slab temperature is highest during bimolecular reaction (m=0.5) and lowest for sensitized reaction (m=-2), hence confirming the earlier results in Table 2. Consequently, oxygen concentration in the slab is lowest during bimolecular reaction and highest for sensitized reaction. In Figures 6 and 7, we observed that the slab temperature decreases while oxygen concentration increases with an increase in the reaction order index (n) and activation energy parameter (ε). This clearly implies that a higher-order exothermic oxidation chemical reaction will be more thermally stable than a lower one. Figures 8 and 9 illustrated the effect of the Frank-Kamenetskii parameter (λ) and the slab upper surface temperature parameter (θw) on the slab temperature and oxygen concentration. The slab temperature increases while the oxygen concentration decreases with an increase in the parameter values of λ and θw. As λ and θw increase, the oxygen consumption within the slab increases and the slab internal heat generation due to exothermic oxidation reaction increases, this invariably leads to an elevation in the slab temperature.Temperature and oxygen concentration profiles form=0.5, n=2, α=β=λ=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5, n=2, α=β=λ=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles forn=β=λ=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,λ=β=1,θw=0,ε=0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,n=λ=β=1,θw=0.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,n=β=1,θw=0,ε= 0.1.
(a)
(b)Temperature and oxygen concentration profiles form=0.5,n=λ=β=1,ε=0.1.
(a)
(b)


## 7.Conclusions

We have computationally investigated the one-stepnth -order oxidation exothermic reaction in a slab with an insulated lower surface and an isothermal upper surface. The model, which consists of a system of coupled heat and mass transfer differential equations, has been solved numerically using a semi-discretization technique and analytical using a perturbation technique coupled with a special type of the Hermite-Padé approximants. Our results revealed, among others, the thermal ignition criticality conditions and with the right combination of thermophysical parameters controlling the system, the thermal runaway can be prevented.



---

*Source: 239052-2012-07-15.xml*