# On New High Order Quasilinearization Approaches to the Nonlinear Model of Catalytic Reaction in a Flat Particle

Authors: S. S. Motsa; O. D. Makinde; S. Shateyi
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

A novel computational approach known as pseudospectral quasilinearization (SQLM) is employed to tackle the two-point boundary value problem describing the reactivity behaviour
of porous catalyst particles subject to both internal mass
concentration gradients and temperature gradients, in
endothermic or exothermic catalytic reactions. A comparison with
the numerical results generated using the inbuilt MATLAB boundary
value solver,bvp4c, for different values of the governing
physical parameters is performed and an excellent agreement is
achieved. A systematic way of improving the convergence of the
SQLM is also presented.


---

## Body

## 1. Introduction

In many engineering and industrial applications, catalytic processes in chemical reactors are often considered to be very useful. This induces particular attention to the study of catalytic reactions at the single-particle level [1]. Moreover, the problem of how the intraparticle diffusion of molecules would modify the overall reaction behaviour of porous catalyst particles had been studied over nearly a quarter of a century [2–4]. Majority of chemical reactions are accompanied by heat transfer effects; they either release or absorb heat. This can lead to appreciable increase (or decrease) of temperature toward the particle centre [5–7]. Since chemical reaction rates vary rapidly (exponentially) with temperature, this effect could radically change the behaviour of the catalyst particles from that which we would otherwise expect. Analysis of chemical kinetics with diffusion effects usually leads to solving highly nonlinear differential equations. Detailed reviews of mathematical models describing reactions in a porous catalyst particle can be found in [8]. Assuming a flat geometry for the particle and that conductive heat transfer is negligible compared to convective heat transfer, Hlavácek et al. [9] derived a dimensionless nonlinear two-point boundary value problem for catalytic reaction in a flat particle as
(1)d2cdy2-λcexp[γβ(1-c)1+β(1-c)]=0,
with boundary condition
(2)dcdy(0)=0,c(1)=1,
where c is the reactant concentration, y is a coordinate measured along the particle, λ is the Thiele modulus or the reaction rate parameter, δ is the activation energy parameter that expresses the sensitivity of the reaction rate to temperature, and β is the heat evolution parameter that shows the maximum temperature variation which could exist within the particle relative to the boundary temperature. The main objective of the current research is to solve the nonlinear problem described in (1)-(2) using the pseudospectral quasilinearization method (SQLM). This method is formed by blending the pseudospectral methods [10–12] with the quadratically convergent Bellman and Kallaba quasilinearization method, [13]. We also present, for the first time, hybrid schemes which converge much faster than the quasilinearization method. We demonstrate the capability of the SQLM in tackling highly nonlinear two-point boundary value problems. Numerical and graphical results obtained using the SQLM approach are validated through comparison with numerical results generated using the inbuilt MATLAB boundary value solver, bvp4c, for different values of the governing physical parameters.


## 2. Solution Methods

To develop the iteration schemes, we begin by rewriting (1) as
(3)c′′+f(c)=0,
where f(c) is the nonlinear component of (1). Given that cδ is the initial approximation of the solution of (3), we introduce the following coupled system:
(4)c′′+f(cδ)+(c-cδ)∂f∂c(cδ)+g(c,cδ)=0,(5)g(c,cδ)=f(c)-f(cδ)-(c-cδ)∂f∂c(cδ).
We write (4) as
(6)c′′+c∂f∂c(cδ)+g(c,cδ)=Φ(cδ),
where
(7)Φ(cδ)=cδ∂f∂c(cδ)-f(cδ).We use the quasilinearization method (QLM) of Bellman and Kalaba [13] to solve (6). The QLM determines the (i+1)th iterative approximation ci+1 as the solution of the differential equation
(8)ci+1′′+ci+1∂f∂c(cδ)+g(ci,cδ)+(ci+1-ci)∂g∂c(ci,cδ)=Φ(cδ),
which can be written as
(9)ci+1′′+[∂f∂c(cδ)+∂g∂c(ci,cδ)]ci+1=ci∂g∂c(ci,cδ)-g(ci)+Φ(cδ),
subject to
(10)ci+1′=0,ci+1(1)=1.
We assume that c0 is obtained as a solution of the linear part of (6) given by(11)c0′′+c0∂f∂c(cδ)=Φ(cδ),
which yields the iteration scheme
(12)cr+1′′+cr+1∂f∂c(cr)=Φ(cr),
which is the standard QLM iteration scheme for solving (1).Wheni=0 in (9), we can approximate c as
(13)c≈c1.
Thus, setting i=0 in (9), we obtain
(14)c1′′+[∂f∂c(cδ)+∂g∂c(c0,cδ)]c1=c0∂g∂c(c0,cδ)-g(c0,cδ)+Φ(cδ),
which yields the iteration scheme
(15)cr+1′′+[∂f∂c(cr)+∂g∂c(cr+1(0),cr)]cr+1=cr+1(0)∂g∂c(cr+1(0),cr)-g(cr+1(0),cr)+Φ(cr),
where cr+1(0) is the solution of
(16)cr+1′′(0)+∂f∂c(cr)cr+1(0)=Φ(cr).The general iteration scheme obtained by settingi=m (m≥2) in (9), hereinafter referred to as scheme-m, is
(17)cr+1′′+[∂f∂c(cr)+∂g∂c(cr+1(m-1),cr)]cr+1=cr+1(m-1)∂g∂c(cr+1(m-1),cr)-g(cr+1(m-1),cr)+Φ(cr),
where cr+1(m-1) is obtained as a solution of
(18)cr+1′′(m-1)+[∂f∂c(cr)+∂g∂c(cr+1(m-2),cr)]cr+1(m-1)=cr+1(m-2)∂g∂c(cr+1(m-2),cr)-g(cr+1(m-2),cr)+Φ(cr).The iterative schemes (12) and (17) can easily be solved using numerical methods such as finite differences, finite elements, Runge-Kutta based shooting methods, or collocation methods. In this work, we employ the Chebyshev spectral collocation method. This method is based on approximating the unknown functions by the Chebyshev interpolating polynomials in such a way that they are collocated at the Gauss-Lobatto points defined as
(19)zj=cosπjN,j=0,1,…,N,
where N is the number of collocation points used (see, e.g., [10, 12]). In order to implement the method, the physical region [0,1] is transformed into the region [-1,1] using the mapping
(20)y=z+12,-1≤z≤1.
The second derivative of ci at the collocation points is represented as
(21)d2cidy2=4∑k=0N‍Dkj2ci(zk),j=0,1,…,N,
where D is the Chebyshev spectral differentiation matrix (see, e.g., [10, 12]). Substituting (19)–(21) in (12), for example, results in the matrix equation
(22)ArCr+1=Φr,
in which Ar is (N+1)×(N+1) square matrix and Cr+1 and Φr are (N+1)×1 column vectors defined by
(23)C=[c(z0),c(z1),…,c(zN-1),c(zN)]T,Φr=[Φ(z0),Φ(z1),…,Φ(zN-1),Φ(zN)]T,Ar=D2+ar.
In the above definitions, ar is a diagonal matrix of size (N+1)×(N+1) which is defined as
(24)ar=∂f∂c(cr)
and is evaluated at the collocation points. After modifying the matrix system (21) to incorporate boundary conditions, the solution is obtained by solving (22) iteratively, starting from a suitable initial approximation.


## 3. Results and Discussion

In this section, we present the results for the solution of the reactant concentrationc(y) which is the solution of the governing nonlinear boundary value problem (1). In order to assess the accuracy of the proposed iteration methods, the present numerical results were compared against results generated using the MATLAB routine bvp4c. The MATLAB routine bvp4c is based on an adaptive Lobatto quadrature scheme [14, 15]. The present results were also compared against other results from the literature that were generated using other methods. Unless otherwise specified, the number of collocation points used in computing the results presented here is N=50. For illustration purposes, results are presented for the first three iteration schemes obtained by setting m=0,1,2.We begin by presenting results for scheme-0 which corresponds to the pseudospectral quasilinearization method (SQLM).Table1 gives a comparison between the SQLM results and the results of Lin et al. [16], who used interval analysis and the validated solver for parametric ordinary differential equations (VSPODE) to produce guaranteed bounds of the solution to (1) for different values of the reactive rate parameter λ, when the activation energy parameter is γ=20 and the heat evolution parameter is β=0.4. Depending on the value of λ, there exist different solutions to (1) as indicated in Table 1. The results are also compared with the bvp4c numerical results. It can be seen from Table 1 that the SQLM results lie within the interval where the true solution for c(y) is expected to be found. We also observe that the SQLM results converge very rapidly to the bvp4c results. It takes only four iterations to achieve an exact match that is accurate to order 10-8 for the selected parameters. This proves the accuracy and efficiency of the proposed SQLM approach.Table 1
Comparison of the present SQLM results forc(0) for β=0.4, γ=20 against the bvp4c numerical results and the VSPODE enclosures of the true solution reported in [16] for various values of λ.
λ
VSPODE enclosure [16]
2nd order QLM
4th order QLM
bvp4c
0.05
[0.970345, 0.970346]
0.97034564
0.97034556
0.97034556
0.1
[0.922680, 0.922681]
0.92275692
0.92268041
0.92268041
[0.505872, 0.505873]
0.50586098
0.50587258
0.50587258
[0.064468, 0.064469]
0.06379374
0.06446821
0.06446821Table2 gives a comparison of the present SQLM results for c(0) for β=0.4, γ=12 against the bvp4c numerical results for various values of λ. Table 2 indicates that the concentration at y=0 gradually decreases with an increase in λ. It can also be seen from the table that there is good agreement between the SQLM and the bvp4c numerical results. The SQLM results converge fully to the numerical results only at 3rd order of approximation. Again, this shows the efficiency of the proposed SQLM approach.Table 2
Comparison of the present SQLM results forc(0) for β=0.4, γ=12 against the bvp4c numerical results for various values of λ.
λ
2nd order
3rd order
4th order
bvp4c
0.04
0.97861570
0.97861566
0.97861566
0.97861566
0.08
0.95387927
0.95387919
0.95387919
0.95387919
0.12
0.92454876
0.92454710
0.92454710
0.92454710
0.16
0.88854336
0.88852609
0.88852609
0.88852609
0.20
0.84197392
0.84188249
0.84188249
0.84188249
0.24
0.77617737
0.77590844
0.77590839
0.77590839
0.28
0.66647431
0.66638660
0.66638660
0.66638660
0.32
0.47282464
0.47282850
0.47282850
0.47282850In Table3, we show the comparison of the present SQLM results for c(0) for γ=12, λ=0.3 against the bvp4c numerical results for various values of β. Table 3 indicates that the concentration at y=0 gradually decreases with an increase in β. It can also be seen from the table that there is good agreement between the SQLM and the bvp4c numerical results. The SQLM results converge fully to the numerical results only at 3rd order of approximation. Again, this shows the efficiency of the proposed SQLM approach. Table 4 gives a comparison between the current SQLM results and the bvp4c results for fixed values of β and λ and various values of γ. It must be pointed out here that a possible advantage of the proposed approach over the bvp4c approach is that it yields accurate results using fewer grid points, iterations, and function evaluations. The higher order scheme further improves the efficiency of the QLM by requiring fewer iterations to achieve better accuracy.Table 3
Comparison of the present SQLM results forc(0) for λ=0.3, γ=12 against the bvp4c numerical results for various values of β.
β
2nd order
3rd order
4th order
bvp4c
0.10
0.84786938
0.84788701
0.84788701
0.84788701
0.15
0.83525321
0.83529531
0.83529531
0.83529531
0.20
0.81920551
0.81926631
0.81926631
0.81926631
0.25
0.79771484
0.79776696
0.79776696
0.79776696
0.30
0.76640940
0.76641380
0.76641380
0.76641380
0.35
0.71295057
0.71291150
0.71291150
0.71291150
0.40
0.57807777
0.57812876
0.57812877
0.57812877Table 4
Comparison of the present SQLM results forc(0) for β=0.4, λ=0.05 against the bvp4c numerical results for various values of γ.
γ
2nd order
3rd order
4th order
bvp4c
2
0.97511266
0.97511280
0.97511280
0.97511280
4
0.97469451
0.97469524
0.97469524
0.97469524
6
0.97425504
0.97425591
0.97425591
0.97425591
8
0.97379354
0.97379273
0.97379273
0.97379273
10
0.97330579
0.97330326
0.97330326
0.97330326
12
0.97278478
0.97278473
0.97278473
0.97278473
14
0.97230716
0.97223388
0.97223388
0.97223388
16
0.97281606
0.97164723
0.97164690
0.97164690Figure1 depicts a slice of bifurcation diagram in the (λ,c(0))-plane for a given value of β>0. It demonstrates the variation in the reacting species concentration along the flat particle centreline with increasing intensity of destructive reaction rate (λ). This result is in perfect agreement with the one reported by Hlavácek et al. [9]. Meanwhile, it is observed that the species concentration along the particle centreline generally decreases with an increase in the reaction rate parameter, since the reactant is consumed. It is worth mentioning that a further decrease in the species concentration is observed with an increase in the reacting species activation energy (γ). Figures 2, 3, and 4 illustrate the transverse variation in the reacting species concentration across the flat particle. It is seen that the concentration is highest at the flat particle surface and decreases transversely with minimum value along the flat particle centreline. Moreover, it is noteworthy that the reaction species concentration decreases with a combined increase in the heat evolution parameter (β), reaction rate parameter (λ), and the activation energy parameter (γ).Figure 1
Bifurcation diagrams for concentration aty=0 against λ when β=0.3 at various values of γ.Figure 2
Comparison of the 8th order SQLM approximate solution for the concentration profile against thebvp4c numerical results when λ=0.4, γ=12 with different values of β.Figure 3
Comparison of the 8th order SQLM approximate solution for the concentration profile against thebvp4c numerical results when β=0.1, λ=0.4 with different values of γ.Figure 4
Comparison of the 8th order SQLM approximate solution for the concentration profile against thebvp4c numerical results when β=0.1, γ=10 with different values of λ.Table5 gives the values of c(0) computed using the proposed higher order hybrid schemes: scheme-0, scheme-1, and scheme-2. It can be seen from Table 5 that there is an improvement in the convergence speed in moving from scheme-0 to scheme-1 then to scheme-2. It takes 4 iterations to obtain a result that converges fully to 18 decimal digits for scheme-0 compared to just 1 iteration when using scheme-2. This demonstrates the benefit in using the proposed higher order schemes.Table 5
Comparison between the iteration schemes for the computed values ofc(0).
Iter.
Scheme-0
Scheme-1
Scheme-2
1
0.862901516195264896
0.973668118653808797
0.973303269931342285
2
0.973668118653808797
0.973303264463730095
0.973303264463730094
3
0.973303269931342285
0.973303264463730094
0.973303264463730094
4
0.973303264463730095
0.973303264463730094
0.973303264463730094
5
0.973303264463730094
0.973303264463730094
0.973303264463730094
6
0.973303264463730094
0.973303264463730094
0.973303264463730094


## 4. Conclusion

In this study, a highly accurate hybrid method that blends the pseudospectral methods of differentiation and quasilinearization techniques was presented. The method of solution, named as the pseudospectral quasilinearization method (SQLM), was employed to tackle the two-point boundary value problem describing the reactivity behaviour of porous catalyst particles subject to both internal mass concentration gradients and temperature gradients, in endothermic or exothermic catalytic reactions. The validity of the SQLM was established by a comparison with the numerical results generated using the inbuilt MATLAB boundary value solver,bvp4c, for different values of the governing physical parameters. A systematic way of improving the convergence of the SQLM based was also presented. From this preliminary investigation, more numerical experimentation will be conducted on other examples of nonlinear BVPs from other applications in science and engineering.



---

*Source: 350810-2013-11-24.xml*