Advances in Mathematical Physics

Volume 2013 (2013), Article ID 350810, 6 pages

http://dx.doi.org/10.1155/2013/350810

Research Article

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

S. S. Motsa,¹ O. D. Makinde,² and S. Shateyi³

¹School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa

²Faculty of Military Science, Stellenbosch University, Private Bag X2, Saldanha 7395, South Africa

³Department of Mathematics & Applied Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, South Africa

Received 28 September 2013; Accepted 8 November 2013

Academic Editor: Raseelo Joel Moitsheki

Copyright © 2013 S. S. Motsa et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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

with boundary condition

where is the reactant concentration, 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

where is the nonlinear component of (1). Given that is the initial approximation of the solution of (3), we introduce the following coupled system:

We write (4) as

where

We use the quasilinearization method (QLM) of Bellman and Kalaba [13] to solve (6). The QLM determines the th iterative approximation as the solution of the differential equation

which can be written as

subject to

We assume that is obtained as a solution of the linear part of (6) given by

which yields the iteration scheme

which is the standard QLM iteration scheme for solving (1).

When in (9), we can approximate as

Thus, setting in (9), we obtain

which yields the iteration scheme

where is the solution of

The general iteration scheme obtained by setting () in (9), hereinafter referred to as scheme-, is

where is obtained as a solution of

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

where is the number of collocation points used (see, e.g., [10, 12]). In order to implement the method, the physical region is transformed into the region using the mapping

The second derivative of at the collocation points is represented as

where is the Chebyshev spectral differentiation matrix (see, e.g., [10, 12]). Substituting (19)–(21) in (12), for example, results in the matrix equation

in which is square matrix and and are column vectors defined by

In the above definitions, is a diagonal matrix of size which is defined as

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 concentration 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 . For illustration purposes, results are presented for the first three iteration schemes obtained by setting .

We begin by presenting results for scheme-0 which corresponds to the pseudospectral quasilinearization method (SQLM).

Table 1 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 and the heat evolution parameter is . 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 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 for the selected parameters. This proves the accuracy and efficiency of the proposed SQLM approach.

Table 1: Comparison of the present SQLM results for for , 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.06446821

Table 2 gives a comparison of the present SQLM results for for , against the bvp4c numerical results for various values of . Table 2 indicates that the concentration at 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 for for , 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.47282850

In Table 3, we show the comparison of the present SQLM results for for , against the bvp4c numerical results for various values of . Table 3 indicates that the concentration at 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 for for , 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.57812877

Table 4: Comparison of the present SQLM results for for , 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.97164690

Figure 1 depicts a slice of bifurcation diagram in the -plane for a given value of . 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 at against when at various values of .

Figure 2: Comparison of the 8th order SQLM approximate solution for the concentration profile against the bvp4c numerical results when , with different values of .

Figure 3: Comparison of the 8th order SQLM approximate solution for the concentration profile against the bvp4c numerical results when , with different values of .

Figure 4: Comparison of the 8th order SQLM approximate solution for the concentration profile against the bvp4c numerical results when , with different values of .

Table 5 gives the values of 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 of .

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.

Conflict of Interests

The authors declare that there is no conflict of interests.

References

J. Kierzenka and L. F. Shampine, “A BVP solver based on residual control and the MATLAB PSE,” ACM Transactions on Mathematical Software, vol. 27, no. 3, pp. 299–316, 2001.

G. Damkohler, “Ubertemperatur in kontaktkornern. (Excess temperature in catalyst grains.),” Zeitschrift für Physikalische Chemie, vol. 193, pp. 16–28, 1943.

D. N. Kim and Y. G. Kim, “An experimental study of multiple steady states in a porous catalyst due to phase transition,” Journal of Chemical Engineering of Japan, vol. 14, pp. 311–317, 1981.

D. H. Kim and Y. G. Kim, “Simulation of multiple steady states in a porous catalyst due to phase transition,” Journal of Chemical Engineering of Japan, vol. 14, no. 4, pp. 318–322, 1981.

D. N. Jaguste and S. K. Bhatia, “Partial internal wetting of catalyst particles: hysteresis effects,” AIChE Journal, vol. 37, pp. 650–660, 1991.

O. D. Makinde, “Exothermic explosions in a Slab: a case study of series summation technique,” International Communications in Heat and Mass Transfer, vol. 31, no. 8, pp. 1227–1231, 2004.

O. D. Makinde, “Strong exothermic explosions in a cylindrical pipe: a case study of series summation technique,” Mechanics Research Communications, vol. 32, pp. 191–195, 2005.

L. D. Landau and Y. M. Lifshits, Statistical Physics, Nauka, Moscow, Russia, 1964, (Russian).

V. Hlavácek, M. Marek, and M. Kubícek, “Modelling of chemical reactors, X: multiple solutions of enthalpy and mass balances for a catalytic reaction within a porous catalyst particle,” Chemical Engineering Science, vol. 23, pp. 1083–1097, 1968.

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, Germany, 1988.

D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, Pa, USA, 1977.

L. N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, Pa, USA, 2000.

R. E. Bellman and R. E. Kalaba, Quasilinearization and Nonlinear Boundary-Value Problems, Elsevier, New York, NY, USA, 1965.

J. Kierzenka and L. Shampine, “A BVP solver based on residual control and the Matlab PSE,” ACM Transactions on Mathematical Software, vol. 27, pp. 299–316, 2001.

L. F. Shampine, I. Gladwell, and S. Thompson, Solving ODEs with MATLAB, Cambridge University Press, Cambridge, UK, 2003.

Y. Lin, J. A. Enszer, and M. A. Stadtherr, “Enclosing all solutions of two-point boundary value problems for ODEs,” Computers & Chemical Engineering, vol. 32, pp. 1714–1725, 2008.