Theory for Electrochemical Heat Sources and Exothermic Explosions: Akbari-Ganji’s Method

1. Introduction

In explosive applications, the safety of combustible materials during transport and storage is essential. Self-ignition is a common occurrence with certain materials. This interior healing occurs when an explosive chemical is heated to a temperature where the decomposition process begins to cause significant exothermic effects. A thermal runaway effect occurs when the temperature rises rapidly, resulting in a fast thermal breakdown. Understanding the components that control this occurrence in many industrial processes is essential.

Many authors have provided models to explain thermal runaway and abusive behaviour. Hatchard et al. [1] make the first attempts to explain the exothermic reaction kinetics. In the works of Peng et al. [3,4] and Kim et al. [2] a simplified reaction-diffusion model is presented. There is also a discussion of an electrochemical model of lithium-ion battery nail absorption [5]. Wang et al. [6] present a catastrophe theoretic approach using a reduced ordinary differential equation model as a foundation. In [7], this model has been expanded. One must identify the primary exothermic chemical reactions to describe the thermal runaway of LIBs. Rising temperature can be used to describe the overall mechanism that causes a thermal runaway [1,3,4,8,9].

The equilibrium voltage and the method by which the equilibrium voltage is derived from the temperature are crucial factors in the electrochemical heat source [10]. The mechanisms of thermal abuse inside the Lithiyam- iron battery are directly related to modelling the thermal runaway and exothermic heat sources. As the temperature rises inside a battery, many exothermic chemical reactions may occur. If there is insufficient heat transmission to the environment [11] or if the rate of heat creation is higher than the heat dissipation rate, this could produce heat that accumulates inside the cell and speeds up the chemical reaction between the components of the cell. Any external factors, including outer heating, high or low voltage/current, nail piercing, exterior shorts, and others, might cause a temperature increase. Cells can be destroyed by a thermal runaway that can be brought on by an air leak, cigarette smoke, gas combustion, fires, etc. In 1930, Semonov, Zeldovith, and Frank Kamenetskii described this behaviour first, and their pioneering contributions are presented in [12]. Frank-Kamenetskii also proposed the steady-state theory of thermal explosion. This idea has been applied to various combustible material geometries in the literature. For the infinite slab, Boddington et al. [13,14,15] reviewed the case of two-step parallel exothermic processes, and they extended their research to include the geometries of the sphere and the circular cylinder.

Graham-Eagle et al. [16] studied a system of exothermic processes co-occurring. The critical values of the Frank-Kamenetskii parameter and the maximum temperature were calculated using a variational technique. The analysis of simultaneous processes [17] was expanded by the same authors to exothermic and endothermic reactions. Gelfand et al. [18] numerically determined the parameter's values.

Ajadi and Gol'dshtein [19] used a three-step kinetics of reactions model including steps for initiation, transmission, and cessation. The computation was performed using an approximation of the effective activation energy. The critical values for several geometries, including infinite square, rod, and cube, have been determined by Balakrishnan et al. [20] employing the finite difference method.

Makinde et al. [21] used the perturbation technique and Hermite-Pade approximants to arrive at the analytical solutions to the governing nonlinear boundary-value problem. They also explored the fundamental characteristics of the thermodynamic field, such as heat and bifurcations. Ananthasamy et al. [22] derived the analytical expression of temperature in exothermic explosions in a slab by solving the nonlinear equation using HAM. It is time-consuming since this contains an infinite number of convergence-control parameters.

Er-Riani and Chetehouna [23] apply the homotopy perturbation method for solving the steady-state nonlinear equation in an exothermic chemical reaction. This method was introduced by He [24,25,26]. In this communication, we derive the simple analytical expression for the temperature field for three kinds of reactions by solving the nonlinear equaftion using Akbari-Ganji’s method.

2. Problem Formulation and Analysis

Considering the steady-state of an exothermic chemical reaction in a combustible slab with the potential of heat loss to the environment. Frank-Kamenetskii [12] first proposed the classical formulation of this problem. The heat balance equation for steady-state conditions is given as follows:

$$k\frac{d^{2}T\left(Y\right)}{d{Y}^2}+Q{C}_{0}A{\left(\frac{KT\left(Y\right)}{vh}\right)}^m{e}^{\frac{-E}{RT\left(Y\right)}}=0$$

(1)

The boundary conditions are

$$\frac{dT}{dY}=0\mathrm{at}Y=0,$$

(2)

$$T=T_0\mathrm{at}Y=a,$$

(3)

The material's thermal conductivity is denoted by the letters k, while T represents the absolute temperature. The other parameters have the same standard meaning as in the reference [12]. Additionally, the numerical exponent for Arrhenius, bimolecular, and sensitised kinetics is m = 2, 0, 1/2. To reduce the complexity, we make the nonlinear Eq. (1) into dimensionless form by defining the following dimensionless parameters.

$$\theta =\frac{E(T-T_0)}{R{T}_0^2},\epsilon =\frac{R{T}_0}{E},y=\frac{Y}{a},\lambda =\frac{QEA{a}^2{C}_0{K}^m{T}_0^{m-2}{e}^{\frac{E}{R{T}_0}}}{v^m{h}^{m}Rk}$$

(4)

where $\theta$is the dimensionless temperature field, $\lambda$ represents the Frank-Kamenetskii constant, $\epsilon$ denotes the activation energy. Using Eq. (4), the Eq. (1) reduces to the following dimensionless form.

$$\frac{d^2\theta \left(y\right)}{d{y}^2}+\lambda {(1+\epsilon \theta \left(y\right))}^m{e}^{\left(\frac{\theta \left(y\right)}{1+\epsilon \theta \left(y\right)}\right)}=0$$

(5)

The boundary conditions (2) and (3) can be reduced as folows:

$$\frac{d\theta }{dy}=0\mathrm{at}y=0,$$

(6)

$$\theta =0\mathrm{at}y=1,$$

(7)

3. Analytical Expression of the Temperature using Akbari-Ganji’s method

Recently, a variety of asymptotic techniques for solving nonlinear differential equations have been established such as the Adomian decomposition [27,28], Variational iteration [29,30], Taylor series [31,32] and Akbari-Ganji method (AGM) [33,34,35,36,37,38,39,40]. Among these techniques, AGM might be recognised as a useful algebraic (semi-analytic) method of resolving such problems. According to the AGM, a solution function with unidentified constant coefficients is supposed to satisfy the differential equation and the initial conditions. Then, the unknown coefficients are calculated using algebraic equations derived from initial conditions and derivatives.

We can assume that the trial solution of the Eq. (5) is

$$\theta \left(y\right)={\sum }_{i=0}^2{\theta }_i{y}^i={\theta }_0+{\theta }_{1}y+{\theta }_2{y}^2$$

(8)

where${\theta }_0,{\theta }_1\mathrm{a}\mathrm{n}\mathrm{d}{\theta }_2$ are constants. We obtain the constant ${\theta }_{0}and{\theta }_1$ using the boundary conditions (6) and (7) as follows:

$${\theta }_1=0,{\theta }_0=-{\theta }_2$$

(9)

Now define the function $G$by

$$G\left(y\right)=\frac{d^2\theta }{d{y}^2}+\lambda {(1+\epsilon \theta )}^m{e}^{\left(\frac{\theta }{1+\epsilon \theta }\right)}=0$$

(10)

Using Eq. (8), the Eq. (10) at $y=0$becomes

$$G\left(y=0\right)={2\theta }_2+\lambda {(1-\epsilon {\theta }_2)}^m{e}^{\left(\frac{{-\theta }_2}{1-\epsilon {\theta }_2}\right)}=0$$

(11)

Using eq. ( 9), the eq. (8) can be rewritten as follows:

$$\theta \left(y\right)={\theta }_2\left(y^2-1\right)$$

(12)

The parameter ${\theta }_2$ is obtained by solving the nonlinear equation

$${2\theta }_2+\lambda {(1-\epsilon {\theta }_2)}^m{e}^{\left(\frac{{-\theta }_2}{1-\epsilon {\theta }_2}\right)}=0$$

(13)

The unidentified parameter ${\theta }_2$ in the equation (13) can be obtained with the Ying Buzu technique [41]. The value of this parameter was also determined using the regular false method and the secant algorithm. The Ying Buzu algorithm approach provides a convergent solution very quickly.

4. Previous analytical results

4.1. Homotopy analysis method

Anathaswamy et al. [11] used the homotopy analysis method (HAM) to solve Eq. (5) with boundary conditions (6) and (7). They obtained that the analytical expression for the temperature field

$$\theta \left(y,\lambda ,\epsilon ,m\right)=\frac{\lambda \left(1-y^2\right)}{2}-h\left[\frac{-{\lambda }^2\left(1+\epsilon m\right)}{2}\left(\frac{y^2}{2}-\frac{y^4}{12}\right)-\frac{{\lambda }^3\epsilon \left(1-m\right)}{4}\left(\frac{y^2}{2}-\frac{y^4}{6}+\frac{y^6}{30}\right)+\frac{{5\lambda }^3\left(1+\epsilon m\right)}{24}+\frac{{11\lambda }^3\epsilon \left(1-m\right)}{120}\right]$$

(14)

It takes a very long time since the convergence-control parameter h are unbounded.

4.2. Perturbation method

Makinde et al [10] solved Eq. (5) with boundary conditions (6) and (7) using the perturbation approach. According to their computations, the temperature was:

$$\theta \left(y,\lambda ,\epsilon ,m=-2\right)=\frac{\lambda (1-y^2)}{2}-\frac{{\lambda }^2}{24}\left(y^2-1\right)\left(y^5-5\right)\left(2\epsilon -1\right)+\frac{{\lambda }^3}{360}\left(y^2-1\right)(-11y^4\epsilon +2y^4+11y^4{\epsilon }^2+64y^2\epsilon -13y^2-64y^2{\epsilon }^2+221{\epsilon }^2-221\epsilon +47)$$

(15)

$$\theta \left(y,\lambda ,\epsilon ,m=0\right)=\frac{\lambda (1-y^2)}{2}+\frac{{\lambda }^2}{24}\left(y^2-1\right)\left(y^5-5\right)+\frac{{\lambda }^3}{360}\left(y^2-1\right)(-2y^4+3y^4\epsilon -12y^2\epsilon +13y^2+33\epsilon -47)$$

(16)

$$\theta \left(y,\lambda ,\epsilon ,m=0.5\right)=\frac{\lambda (1-y^2)}{2}-\frac{{\lambda }^2}{48}\left(y^2-1\right)\left(y^5-5\right)(2+\epsilon )+\frac{{\lambda }^3}{1440}\left(y^2-1\right)(4y^4\epsilon -8y^4+y^4{\epsilon }^2+4y^2\epsilon +52y^2+y^2{\epsilon }^2-14{\epsilon }^2-56\epsilon -188)$$

(17)

These approximate solutions are obtained by carrying out expansions in terms of a small parameter $\left(\epsilon \right)$. In addition to this, perturbation method is not always to generate a continuous family of solutions in terms of the small parameter.

5. Discussion

Equation (12) represents the simplest analytical expression of the temperature profile. The thermal decomposition of the reacting combustible material depends on the parameters $\lambda \epsilon andm$, which are of great importance for applications in explosives handling and industrial safety.

5.1. Numerical Simulation

This section validates the above theoretical results using numerical simulation for physically realistic values of various embedded parameters. The function bvp4c in Matlab/Scilab software, which solves non-linear boundary value problems for ordinary differential equations, is used to solve these equations numerically. The present (AGM) and previous (HAM, PM) analytical results are compared to this numerical solution in Table 1, Table 2, Table 3 and Table 4. The maximum average relative error between our result and the simulation result is 1.44. But the maximum error between numerical results and HAM and PM is 1.86% and 4.92%.

6. Limiting case

When the activation energy $\left(\epsilon \right)$ or Arrhenius kinetics (m) is very small, the equation (5) becomes

$$\frac{d^2\theta }{d{y}^2}+\lambda e^{\theta }=0$$

(18)

In this case, the temperature becomes,

$$\theta \left(y\right)=log\left[\frac{n-n{Tanh}^2\left(0.5y\sqrt{n}\right)}{2\lambda }\right]$$

(19)

where n is obtained from the equation

$$n-n{Tanh}^2\left(0.5\sqrt{n}\right)=2\lambda$$

(20)

We can notice that the average error percentage between AGM and an exact limiting case result (eqn.(19)) does not exceed 1.8% for the slab.

Table 5. Comparison of our approximate analytical result (12) with exact result (19) for the limiting case.

$y$	$\lambda =0.1$			$\lambda =0.3$			$\lambda =0.5$
	Exact solution	AGM Eq. (12)	Error%	Exact solution	AGM Eq. (12)	Error%	Exact solution	AGM Eq. (12)	Error%
0	0.0522	0.0527	0.9578	0.1733	0.1795	3.5776	0.3290	0.3474	5.5927
0.2	0.0501	0.0505	0.7984	0.1661	0.1689	1.6857	0.3148	0.3246	3.1131
0.4	0.0436	0.0439	0.6881	0.1443	0.1440	0.2079	0.2728	0.2768	1.4663
0.6	0.0329	0.0330	0.3039	0.1085	0.1087	0.1843	0.2040	0.2060	0.9804
0.8	0.0180	0.0180	0.0000	0.0590	0.0590	0.0000	0.1102	0.1102	0.0000
1	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
	Average error (%)		0.4580	Average error (%)		0.9426	Average error (%)		1.8587

Table 5. Comparison of our approximate analytical result (12) with exact result (19) for the limiting case.

$y$	$\lambda =0.1$			$\lambda =0.3$			$\lambda =0.5$
	Exact solution	AGM Eq. (12)	Error%	Exact solution	AGM Eq. (12)	Error%	Exact solution	AGM Eq. (12)	Error%
0	0.0522	0.0527	0.9578	0.1733	0.1795	3.5776	0.3290	0.3474	5.5927
0.2	0.0501	0.0505	0.7984	0.1661	0.1689	1.6857	0.3148	0.3246	3.1131
0.4	0.0436	0.0439	0.6881	0.1443	0.1440	0.2079	0.2728	0.2768	1.4663
0.6	0.0329	0.0330	0.3039	0.1085	0.1087	0.1843	0.2040	0.2060	0.9804
0.8	0.0180	0.0180	0.0000	0.0590	0.0590	0.0000	0.1102	0.1102	0.0000
1	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
	Average error (%)		0.4580	Average error (%)		0.9426	Average error (%)		1.8587

7. Influence of the parameters on temperature

7.1. Effect of the Frank-Kamenetskii parameter $\lambda$ on temperature

Figure 1 illustrates the effects of the Frank-Kamenetskii parameter on a temperature profile. Figure 2 shows that an increase in Frank-Kamenetskii leads to an increase in the rate of an exothermic reaction. The maximum temperature along the slab is at the centre line, and the minimum is at the slab surface.. The slab temperature will invariably rise as a result of this.

7.2. Effect of the numerical exponent m on temperature

The effects of numerical exponent m on temperature profiles are shown in Figure 2 (a). The figure shows that the temperature increases as the numerical exponent number m increases. Moreover, the tables and figures show that in a bimolecular (m = 0.5) reaction, a thermal explosion occurs faster than in the Arrhenius (m = 0), sensitised (m = -2) reactions.

7.3. Effect of the activation energy parameter (ε) on temperature

A similar effect of temperature is observed with increasing values of activation energy parameter (ε) (Figure 2-b). With rising values of the activation energy, a similar effect of temperature enhancement is observed. Increasing values of ε imply that the reacting slab's activation energy is insufficient, and thus the reacting slab's volatility characteristic is significantly reduced.

8. Conclusions

We studied the exothermic explosion of a viscous combustible in the slab under various laws of reactions.The exothermic chemical reaction in a slab of flammable material is considered. The expected outcomes reveal that this method provides an outstanding and highly accurate approximation of the solution of this nonlinear system. The effects of the parameters Frank-Kamenetskii, numerical exponent of temperature and activation energy on temperature profiles are discussed. This method can be applied to infinite cylindrical and spherical geometries.