# Analytical Solution of Nonlinear Dynamics of

Authors: Felicia Shirly Peace; Narmatha Sathiyaseelan; Lakshmanan Rajendran
Journal: International Journal of Chemical Engineering
 (2014)

---

## Abstract

A mathematical model of the dynamics of the self-ignition of a reaction-diffusion system is studied in this paper. An approximate analytical method (modified Adomian decomposition method) is used to solve nonlinear differential equations under steady-state condition. Analytical expressions for concentrations of the gas reactant and the temperature have been derived for Lewis number (Le) and parametersβ, γ, and ϕ2. Furthermore, in this work, the numerical simulation of the problem is also reported using MATLAB program. An agreement between analytical and numerical results is noted.


---

## Body

## 1. Introduction

Nonlinear dynamical phenomena in combustion process are an active area of experimental and theoretical research. Mathematical models that describe this phenomenon can be considered as nonlinear dynamical systems. The dynamic characterization of such models was developed by Continillo et al. [1]. The detailed numerical simulation of autoignition of coal stockpiles leads to the observation of steady regimes. To better investigate this phenomenon, two simplified distributed-parameter models were discussed which incorporate heat conduction, mass diffusion, and one-step Arrhenius exothermic chemical reaction. Both model equations were solved with straight forward finite-difference schemes [2]. The problem of spontaneous ignition of coal stockpiles is challenging for safety implications and for its theoretical complexity: a spontaneous combustion reaction takes place in a bed of solid fuel, while flow, driven by natural convection generated by the onset of temperature gradients within the pile, occurs. Coal stockpiles self-ignite when reaction of coal with oxygen present in the atmosphere generates heat, that is, not efficiently removed toward the external ambient [3]. Continillo et al. [4, 5] have analyzed the self-combustion of coal piles in the absence of natural convection. The three main phenomena in the self-ignition of coal stockpile are convection, reaction, and diffusion. On the other hand, Continillo et al. [6] studied the dynamic behavior of a two-dimensional coal pile also by accounting for natural convection.As part of a comprehensive study of self-heating of coal stockpiles, a simple mathematical model has been developed. To our knowledge, no rigorous analytical expressions of gas reactant (Y) and temperature (T) have been derived for all possible values of parameters under steady-state conditions. The purpose of this paper is to derive approximate analytical expressions for gas reactant concentration and temperature using the modified Adomian decomposition method.


## 2. Mathematical Formulation of the Boundary Value Problem

The nonlinear differential equations are those of a distributed-parameter dynamic model of heterogeneous reaction in a one-dimensional layer. The gaseous reactant diffuses through the reacting medium and a first order one-step exothermic chemical reaction takes place. The reaction rate depends on the temperature through the Arrhenius exponential. The Arrhenius rate equation is a mathematical expression which relates the rate constant of a chemical reaction to the exponential value of the temperature. The model nonlinear equations in dimensionless form are [1]
(1)∂Y∂t=Le∂2Y∂x2-ϕ2Yexp⁡(-γT),∂T∂t=∂2T∂x2+βϕ2Yexp⁡(-γT),
where Y is the concentration of the gas reactant, T is the temperature, Le is the Lewis number (the ratio between mass and heat diffusivities),  β  is the dimensionless heat of reaction, ϕ is the thermal Thiele modulus (the ratio of the time scale of the limiting transport mechanism to the time scale of intrinsic reaction kinetics [7]), and γ is the dimensionless activation energy (minimum amount of energy between reactant molecules for effective collisions between them). The boundary conditions are
(2)T(0,t)=T(1,t)=1,Y(0,t)=1,∂Y∂x|x=1=0fort>0.
Under steady-state condition, the equations become
(3)d2Ydx2-ϕ2LeYexp⁡(-γT)=0,d2Tdx2+βϕ2Yexp⁡(-γT)=0
with boundary conditions
(4)T=1,Y=1,atx=0,T=1,dYdx=0atx=1.


## 3. Analytical Solution of Nonlinear Dynamics of a Self-Igniting Reaction-Diffusion System under Steady-State Condition Using Modified Adomian Decomposition Method

In the recent years, much attention is devoted to the application of the Adomian decomposition method to the solution of various scientific models [8]. An efficient modification of the standard Adomian decomposition method for solving initial value problem in the second order partial differential equation yields the MADM. The MADM without linearization, perturbation, transformation, or discretization gives an analytical solution in terms of a rapidly convergent infinite power series with easily computable terms. The results show that the rate of convergence of modified Adomian decomposition method is higher than standard Adomian decomposition method [9–13]. Using this method (see Appendix A), we obtain an approximate analytical expression of concentration of gas reactant (Y) and temperature (T) (see Appendix B) as follows:
(5)Y(x)=1+ϕ2Leexp⁡(-γ)(x22-x)+ϕ4Leexp⁡(-2γ)×[1Le(x424-x36+x3)-γβ2(x412-x36+x6)],(6)T(x)=1-βϕ2exp⁡(-γ)(x22-x2)-βϕ4exp⁡(-2γ)×[1Le(x424-x36+x8)-γβ2(x412-x36+x12)].


## 4. Numerical Simulation

Nonlinear diffusion equation (3) for the boundary condition (4) is also solved numerically. We have used the function pdex1 in MATLAB software to solve the initial-boundary value problems for the nonlinear differential equations numerically. This numerical solution is compared with our analytical results in Figures 1–4. Upon comparison, it gives a satisfactory agreement for all values of the dimensionless parameters  β,  γ, and ϕ2. The MATLAB program is also given in Algorithm 1.Algorithm 1:MATLAB program to find the numerical solution of nonlinear differential equation (1) when the time (t) is large (or) (3).
function pdex4
m = 0;
x = linspace(0,1);
t = linspace(0,100000);
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
figure
plot(x,u1(end,:))
title(‘u1(x,t)')
xlabel(‘Distance x')
ylabel(‘u1(x,2)')
%––––––––––––––––––––––––––––––––––
figure
plot(x,u2(end,:))
title(‘u2(x,t)')
xlabel(‘Distance x')
ylabel(‘u2(x,2)')
%––––––––––––––––––––––––––––––––––
function[c,f,s] = pdex4pde(x,t,u,DuDx)
c =[1; 1];
f =[1; 1].* DuDx;
l = 0.233;
p = 35000;
b = 4.287;
g = 13.6;
F = −p*u(1)*exp(−g/u(2))/l;
F1 = b*p*u(1)*exp(−g/u(2));
s =[F; F1];
%––––––––––––––––––––––––––––––––––
function u0 = pdex4ic(x);
u0 =[1; 1];
%––––––––––––––––––––––––––––––––––
function[pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl =[ul(1)−1; ul(2)−1];
ql =[0; 0];
pr =[0; ur(2)−1];
qr =[1; 0];Figure 1
Dimensionless concentrationY versus dimensionless spatial coordinatexusing (5) for Le = 0.233, β = 4.287,  γ = 13.6, and various values of ϕ2. Solid lines represent the analytical solution, whereas the dotted lines are for the numerical solution.Figure 2
Dimensionless temperatureT versus dimensionless spatial coordinate using (6) for Le = 0.233, β = 4.287,  γ =13.6, and various values of ϕ2. Solid lines represent the analytical solution, whereas the dotted lines are for the numerical solution.Figure 3
Dimensionless concentrationY versus dimensionless spatial coordinatexusing (5) for Le = 0.233, β = 4.287,ϕ2 = 1000, and various values of  γ. Solid lines represent the analytical solution, whereas the dotted lines are for the numerical solution.Figure 4
Dimensionless temperatureT versus dimensionless spatial coordinatexusing (6) for Le = 0.233, β = 4.287 ϕ2 = 1000, and various values of  γ. Solid lines represent the analytical solution, whereas the dotted lines are for the numerical solution.


## 5. Discussion

Equations (5) and (6) represent the simplest form of approximate analytical expressions for the concentration of gas reactant and temperature for all values of parameters α,  β,  γ, and ϕ2. The Thiele number (thermal) is the ratio of layer thickness (L) and thermal diffusivity (α). Equation (5) represents the new approximate analytical expression of concentration of gas reactant. The numerical solution is compared with the analytical results in Figures 1–4. These figures represent the analytical and numerical concentration profiles of gas reactant and temperature for different values of parameters  β,γ,and ϕ2. Figure 1 represents the dimensionless concentration Y versus dimensionless spatial coordinate x for ϕ2≤35000. From the figure it is inferred that the value of Y decreases when the value of ϕ2 or layer thickness increases. Figure 2 illustrates the dimensionless concentration Y versus dimensionless spatial coordinate x for ϕ2≤35000 and we infer that the dimensionless temperature increases with the increase in values of layer thickness. Figure 3 represents the dimensionless concentration Y versus dimensionless spatial coordinate x for values of γ≤13.6. From the figure it is inferred that the value of Y decreases when the value of dimensionless activation energy (γ) increases. Figure 4 illustrates the dimensionless temperature T versus dimensionless spatial coordinate x for the values of γ≤13.6; the dimensionless temperature increases with increase in values of dimensionless activation energy (γ). From Figures 2 and 4 it is evident that the maximum value for dimensionless concentration is 1 and that temperature attains its maximum value when the spatial coordinate  x=0.5. Figure 5 confirms the results given by Figures 1 to 4.(a) The normalized three-dimensional concentration of gas reactantY versus dimensionless spatial coordinate x and Thiele number ϕ2, (b) the normalized three-dimensional concentration of gas reactant Y versus dimensionless spatial coordinate x and dimensionless activation energy γ, (c) the normalized three-dimensional temperature T versus dimensionless spatial coordinate x and Thiele number ϕ2, and (d) the normalized three-dimensional temperature T versus dimensionless spatial coordinate x and dimensionless activation energy γ.
(a)
(b)
(c)
(d)Our analytical results are compared with the numerical results for the dimensionless concentrationY in Table 1. The maximum relative error between our analytical results and simulation results for the concentration Y is 1.3%. Also in Table 2, our analytical results are compared with the numerical results for the dimensionless temperature T. Satisfactory agreement is noted. The maximum relative error in this case is 0.4%.Table 1
Comparison between analytical and numerical values of dimensionless concentrationY (Le=0.233, β=4.287, γ=13.6, and ϕ2=30000).
x
Numerical value ofY
Analytical value ofY
% deviation
0
1
1
0
0.2
0.9661
0.9713
0.53
0.4
0.9398
0.9489
0.97
0.6
0.9216
0.9329
1.23
0.8
0.9113
0.9233
1.32
1
0.9082
0.9201
1.31Table 2
Comparison between analytical and numerical values of dimensionless temperatureT (Le=0.233, β=4.287, γ=13.6, and ϕ2=30000).
x
Analytical value ofT
Numerical value ofT
% deviation
0
1
1
0
0.2
1.0155
1.1013
0.27
0.4
1.0235
1.0191
0.42
0.6
1.0233
1.0191
0.41
0.8
1.0152
1.0128
0.24
1
1
1
0


## 6. Conclusion

The system of steady-state nonlinear differential equations in different dynamic modes for the concentration of gas reactant and temperature has been solved analytically. The model investigated the influence of parameters over the temperature and concentration of gas reactant in the dynamic mode. The approximate analytical expression for the steady-state concentration of gas reactant and temperature is obtained using the modified Adomian decomposition method. A satisfactory agreement with the numerical result is noted. The analytical results will be useful to characterize the model predictions for the various values of parameters Thiele number (ϕ2), Lewis number (Le), layer thickness (L), and thermal diffusivity (α).



---

*Source: 825797-2014-04-24.xml*