Mathematical Problems in EngineeringVolume 2011 (2011), Article ID 489217, 19 pagesdoi:10.1155/2011/489217

Research Article

On New Numerical Techniques for the MHD Flow Past a Shrinking Sheet with Heat and Mass Transfer in the Presence of a Chemical Reaction

Z. G. Makukula,¹ P. Sibanda,¹ S. S. Motsa,¹ and S. Shateyi²

¹School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa²Department of Mathematics and Applied Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, South Africa

Received 29 July 2011; Revised 22 August 2011; Accepted 23 August 2011

Academic Editor: Oleg V. Gendelman

Copyright © 2011 Z. G. Makukula 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.

We use recent innovative solution techniques to investigate the problem of MHD viscous flow due to a shrinking sheet with a chemical reaction. A comparison is made of the convergence rates, ease of use, and expensiveness (the number of iterations required to give convergent results) of three seminumerical techniques in solving systems of nonlinear boundary value problems. The results were validated using a multistep, multimethod approach comprising the use of the shooting method, the Matlab bvp4c numerical routine, and with results in the literature.

1. Introduction

Boundary layer flow over a stretching surface occurs in several engineering processes such as hot rolling, wire drawing, and glass-fibre production. Materials that are manufactured by extrusion processes and heat-treated substances proceeding between a feed roll and a wind-up roll can be classified as a continuously stretching surface [1–3]. A shrinking film is useful in the packaging of bulk products since it can be unwrapped easily with adequate heat [4–7]. Shrinking problems can also be applied to the study of capillary effects in small pores and the hydraulic properties of agricultural clay soils [8]. Studies of flow due to a shrinking sheet with heat transfer and/or mass transfer have been considered by, among others, [7, 9].

In recent years, several analytical or semianalytical methods have been proposed and used to find solutions to most nonlinear equations. These methods include the Adomian decomposition method (ADM) [10, 11], differential transform method (DTM) [12], variational iteration method (VIM) [13], homotopy analysis method (HAM) [14–17], and Homotopy perturbation method (HPM) [18–23].

Motsa and Shateyi [24] obtained a numerical solution of magnetohydrodynamic (MHD) and rotating flow over a porous shrinking sheet by the new approach known as spectral homotopy analysis method (SHAM). Muhaimin et al. [5] studied magnetohydrodynamic viscous flow due to a shrinking sheet in the presence of suction. The study found out that the shrinking of the sheet has a substantial effect on the flow field and, thus, on the heat and mass transfer rate from the sheet to the fluid.

In this paper we provide a qualitative assessment of key features of three recent seminumerical techniques, namely, the successive linearisation method (SLM), the spectral-homotopy analysis method (SHAM), and the improved spectral-homotopy analysis method (ISHAM). The two methods were introduced and used by Motsa and his coworkers (see Motsa et al. [25, 26] and Makukula et al. [27–30]) to solve nonlinear boundary value problems. In Motsa et al. [25, 26, 29] the SHAM approach was tested on simple one-dimensional nonlinear boundary value problems. Later, Makukula et al. [28, 30, 31] extended the application of the SHAM to a system of two coupled nonlinear equations that model the von Kármán fluid flow problem. The SLM method was applied on one-dimensional nonlinear differential equations in Makukula et al. [27]. In this study we solve the nonlinear equations that govern the shrinking sheet problem for purposes of evaluating the efficiency of each method with regards to speed of convergence, ease of use, and expensiveness (in terms of the number of iterations required to give convergent results). We introduce the ISHAM as a method that is meant to improve the accuracy of the standard SHAM approach. The governing equations for the problem are a rather formidable system of three nonlinear differential equations in three unknowns. Parametric study of the effect of different parameters is made and the results compared with previous findings in the literature (see Noor et al. [6], Mohd and Hashim [7], and Muhaimin et al. [5]). The solutions are further compared with results obtained using the shooting method and the bvp4c solver, which is based on Runge-Kutta fourth-order schemes.

2. Mathematical Formulation

We investigate the effect of chemical reaction, heat and mass transfer on nonlinear MHD boundary layer past a porous shrinking sheet with suction. The governing boundary layer equations of momentum, energy, and mass diffusion in terms of the velocity components , , and are (see Muhaimin et al. [5])

where is the thermal conductivity of the fluid, B₀ is the magnetic field, is the thermal viscosity, is the permeability of the porous medium, is the rate of chemical reaction, is the kinetic viscosity, is the dynamic viscosity, and is the electrical conductivity.

The applicable boundary conditions are where is the shrinking constant and is the suction velocity. The cases and correspond to shrinking sheets in the - and -directions, respectively.

Using the similarity transformations (see Sajid and Hayat [32]): (2.1) are transformed to the system of nonlinear equations subject to where is the Prandtl number, is the Schmidt number, is the porosity, and is the chemical reaction parameter. We remark that (2.4) can be solved independently of equations of (2.5)-(2.6) for , but the solutions for and directly depend on the solution for . To demonstrate how robust the proposed methods of solution are, the system of (2.4)-(2.5) is solved simultaneously in the next section. Solving the equations simultaneously is also important when conducting the parametric study because some of the governing parameters such as and affect all three unknown variables.

3. Solution Methods

We solve (2.4)–(2.6) using three recent innovative semi-numerical methods. Validation of the results is done by further solving the equations numerically using a shooting method and the Matlab bvp4c solver. For the last two methods we used a tolerance of .

We begin by transforming the domain to , using the domain truncation method, the domain is first approximated by the computational domain , where is a fixed length that is taken to be larger than the thickness of the boundary layer. The domain is then transformed to using the algebraic mapping

3.1. The Successive Linearisation Method (SLM)

The successive linearisation method (see Makukula et al. [27, 28]) is used to solve (2.4)–(2.7). The starting point is to assume that the independent variables , , and may be expanded as where , and , are unknown functions and , , and are approximations that are obtained by recursively solving the linear part of the equation that results from substituting (3.2) in the governing equations (2.4)–(2.7). Substituting (3.2) in the governing equations (2.4)–(2.7) gives where the coefficient parameters , , (), , , and are defined as

Starting from the initial approximations which are chosen to satisfy the boundary conditions (2.7), the subsequent solutions for , , and , , are obtained by successively solving the linearized form of (3.3) which are subject to the boundary conditions Once each solution for , , and () has been found from iteratively solving (3.6)-(3.7) for each , the approximate solutions for , , and are obtained as In coming up with (3.8), we have assumed that

Equations (3.6)-(3.7) are integrated using the Chebyshev spectral collocation method (Canuto et al. [33] and Trefethen [34]). The unknown functions are defined by the Chebyshev interpolating polynomials with the Gauss-Lobatto points defined as where is the number of collocation points used. The unknown functions , , and are approximated at the collocation points by where is the th Chebyshev polynomial defined as The derivatives at the collocation points are represented as where is the order of differentiation and with being the Chebyshev spectral differentiation matrix. Substituting (3.13) in (3.6)-(3.7) leads to the matrix equation where is a square matrix and and are column vectors defined by

In the above definitions, , and , () are diagonal matrices of size . After modifying the matrix system (3.14) to incorporate the boundary conditions, the solution is obtained as

3.2. Spectral-Homotopy Analysis Method (SHAM)

The spectral-homotopy analysis method (SHAM) has been used by Motsa et al. [25, 26]. It is also convenient to first ensure that the boundary conditions are made homogeneous by using the transformations where , and , are chosen to satisfy the boundary conditions (2.7) of the governing equations (2.4)–(2.6). From (3.1) and the chain rule, we have that Substituting (3.1) and (3.17)-(3.18) in the governing equations and boundary conditions gives where prime now denotes derivative with respect to and

The initial guesses used are Solving the linear part of the equation system (3.19), that is, subject to will yield the initial SHAM approximate solution. Applying the Chebyshev pseudospectral method on equations (3.22)-(3.23) yields the matrix form whereThe superscript denotes the transpose, is a diagonal matrix, and is an identity matrix of size . The boundary conditions (3.23) are implemented in matrix and vector of equation (3.24). The values of are then determined from the following equation: which provides us with the initial approximation for the solution of the governing equations (3.19). With the initial approximate solution, we then find approximate solutions for the nonlinear equations (3.19). We start by defining the following linear operators: where is the embedding parameter and , , and are unknown functions. The zeroth-order deformation equations are given by where is the nonzero convergence controlling auxiliary parameter and , , and are nonlinear operators given by The th order deformation equations are given by subject to the boundary conditions where Applying the Chebyshev pseudospectral transformation to equations (3.31)–(3.33) gives rise to the matrix equation subject to the boundary conditions where and are as defined in (3.25) and Applying the boundary conditions (3.32) on the right-hand side of (3.35) yields the following recursive formula for higher-order approximations for :

3.3. Improved Spectral-Homotopy Analysis Method (ISHAM)

Details of the improved spectral-homotopy analysis method (ISHAM) can be found in Makukula et al. [30]. The main objective is to improve the convergence rate of the spectral-homotopy analysis method by using an optimal initial approximation. Hence, instead of a random solution choice a systematic approach is employed to find the optimal initial approximation. This is achieved by first assuming that the solutions , , and can be expanded into where , , and are unknown functions whose solutions are obtained using the SHAM approach at the th iteration and , , and () are known from previous iterations. We use the same initial guesses as with the SHAM solution in Sections 3.1 and 3.2. Substituting (3.39) into the governing equations gives subject to the boundary conditions The coefficient parameters , , (, , , and are as defined in equation (3.4). Starting from the initial guesses (3.5), the subsequent solutions , , and () are obtained by recursively solving (3.40) using the SHAM approach. To find the SHAM solutions of (3.40), we start by defining the following linear operators: The zeroth-order deformation equations are given by, , and are nonlinear operators given by The th order deformation equations are subject to the boundary conditions where The initial approximations , , and that are used in the higher-order equations (3.45)–(3.47) are obtained by solving the linear part of (3.40) given by with the boundary conditions In a similar manner, we apply the spectral methods to solve for the initial approximate solutions , , and , and the higher-order deformation equations (3.45)–(3.47) for higher order approximate solutions , , and for . The solutions for , , and are then generated using the solutions for , , and as follows: The approximate solutions for , , and are then obtained by substituting , , and from (3.50) into (3.39), where is the th iteration of the higher-order deformation equation and is the th iteration of the initial approximation.

4. Results and Discussion

Equations (2.4)-(2.6) subject to boundary conditions (2.7) have been solved using three recent semi-numerical techniques as described above. To validate our results, we have compared the skin friction coefficient, the Nusselt number, and the Sherwood number with the theoretical results of Muhaimin et al. [9]. We have further compared our results with the full numerical solutions obtained using the shooting method and the Matlab bvp4c routine. The comparison is given in Tables 1–3.

Table 1: Comparison of the approximate solutions of at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of when s = 3, M = 1, m=1, Sc = 0.62, γ =3, Pr = 1, λ = 0, ℏ=−1, L = 30, and N = 150.

SLM

SHAM

ISHAM

Shooting
bvp4c
Ref. [9]

order

order

order

1
1
3.32068
1
3.30709

3.33858
3.30278
3.30278
3.302776

2
3.30283
2
3.30338

3.30278

3
3.30278
4
3.30279

3.30278

4
3.30278
6
3.30278

3.30278

2
1
3.57292
1
3.56643

3.59883
2.30278
2.30278
3.561553

2
3.56157
2
3.56225

2.30278

3
3.56155
4
3.56157

2.30278

4
3.56155
6
3.56155

2.30278

4
1
4.00770
1
4.00549

4.03789
4.00000
4.00000
4.000000

2
4.00000
2
4.00084

4.00000

3
4.00000
4
4.00002

4.00000

4
4.00000
6
4.00000

4.00000

Table 2: Comparison of the approximate solutions of at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of when s = 3, M = 1, m = 1, Sc = 0.62, γ = 3, Pr = 1, λ = 0, ℏ = −1, L = 30, and N = 150.

SLM

SHAM

ISHAM

Shooting
bvp4c
Ref. [9]

order

order

order

1
1
2.56783
1
2.67511

2.80912
2.66554
2.66554
2.665537

2
2.66485
2
2.66656

2.66554

3
2.66554
4
2.66556

2.66554

4
2.66554
6
2.66554

2.66554

2
1
2.59198
1
2.69041

2.84009
2.68032
2.68032
2.680315

2
2.67987
2
2.68132

2.68032

3
2.68032
4
2.68034

2.68032

4
2.68032
6
2.68032

2.68032

4
1
2.62914
1
2.71316

2.88452
2.70240
2.70240
2.702455

2
2.70215
2
2.70336

2.70240

3
2.70240
4
2.70241

2.70240

4
2.70240
6
2.70240

2.70240

Table 3: Comparison of the approximate solutions of at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of when s = 3, M = 1, m = 1, Sc = 0.62, γ = 3, Pr = 1, λ = 0, ℏ = −1, L = 30, and N = 150.

SLM

SHAM

ISHAM

Shooting
bvp4c
Ref. [9]

order

order

order

1
1
2.39294
1
2.41413

2.43976
2.41029
2.41029
2.410283

2
2.41026
2
2.41085

2.41029

3
2.41029
4
2.41030

2.41029

4
2.41029
6
2.41029

2.41029

2
1
2.40181
1
2.42137

2.44981
2.41700
2.41700
2.417000

2
2.41698
2
2.41764

2.41700

3
2.41700
4
2.41702

2.41700

4
2.41700
6
2.41700

2.41700

4
1
2.41559
1
2.43236

2.46460
2.42722
2.42722
2.427225

2
2.42721
2
2.42797

2.42722

3
2.42722
4
2.42724

2.42722

4
2.42722
6
2.42722

2.42722

Tables 1–3 give values of the skin friction, heat transfer rate, and the mass transfer rate, respectively, for different porosity values. The convergence to the two numerical results of the SLM is achieved at the third order of approximation, at the sixth order for the SHAM, and at second order for the ISHAM. Comparison with results reported in Muhaimin et al. [9] shows an excellent agreement.

Table 1 shows an increase in the surface shear stress with an increase in the porosity parameter . The increase in the skin friction with the porosity may be accounted for by the fact that the velocity gradient increases with porosity (Takhar et al. [35]). Tables 2 and 3 show an increase in the surface heat transfer rate and the mass transfer rate with the porosity parameter for large suction values (), suggesting an increase in temperature and concentration gradients with increasing porosity.

Figure 1 serves two purposes: (a) to give sense of the accuracy of the improved spectral homotopy analysis (ISHAM) by means of a comparison between the numerical results and the second-order improved spectral-homotopy analysis results and (b) to demonstrate the effects of the suction parameter and the Hartmann number on the velocity profiles .

(a)

(b)

Figure 1: On the comparison between the 2nd-order ISHAM solution (figures) and the bvp4c numerical solution (solid line) for f(η) and θ(η) at different values of when , , , , , , and .

Firstly we observe an excellent agreement between the second-order ISHAM and the numerical bvp4c results for all parameter values used. Secondly we note that these results are qualitatively similar to those reported in Noor et al. [6] for the case of one-direction shrinking and show that increasing the suction parameter and the Hartmann number leads to an increase in the velocity. This in turn leads to a decrease in the boundary layer thickness as fluid is sucked out of the flow region.

5. Conclusions

We have successfully solved the nonlinear system of equations governing MHD boundary layer past a porous shrinking sheet with a chemical reaction and suction. We demonstrated three recent innovative methods, namely, the successive linearisation method (SLM), the spectral-homotopy analysis method (SHAM), and the improved spectral-homotopy analysis method (ISHAM), and compared the performance of the three methods with regard to the speed of convergence of the solution (the number of iterations required), computational efficiency, and the ease of application of the method. The results were compared with those obtained using the well-known shooting method and the Matlab bvp4c solver. We found that the ISHAM converged at second order. The magnitude of the parameter values used did not affect its performance under the same conditions with the SLM and SHAM. Nevertheless, the ISHAM does not come cheap in terms of the size of the code and computer time, taking about three times as long as the SLM to compute the same result and about double the time taken with the SHAM. The SLM converged at third order, is easy to implement, and has shown a good level of stability when solving highly nonlinear problems. The SHAM gives good convergence under the same conditions but poor convergence with highly nonlinear problems. It is easy to implement but not as easy as with the SLM.

Results from simulations revealed an excellent agreement between results from the shooting method and the bvp4c. Our findings indicate that the ISHAM is the best approach of the three methods in terms of the accuracy of the results and speed of convergence. Parametric studies for effects of different parameter values in the problems agreed with results present in the literature.

Acknowledgments

The authors wish to acknowledge financial support from the University of KwaZulu-Natal, University of Venda, and the National Research Foundation (NRF).

References

N. Bachok, A. Ishak, and I. Pop, “Unsteady three-dimensional boundary layer flow due to a permeable shrinking sheet,” Applied Mathematics and Mechanics, vol. 31, no. 11, pp. 1421–1428, 2010.

T.-G. Fang, J. Zhang, and S. S. Yao, “Slip magnetohydrodynamic viscous flow over a permeable shrinking sheet,” Chinese Physics Letters, vol. 27, no. 12, article 124702, 2010.

T. Fang, S. Yao, J. Zhang, and A. Aziz, “Viscous flow over a shrinking sheet with a second order slip flow model,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 7, pp. 1831–1842, 2010.

T. Fang and J. Zhang, “Closed-form exact solutions of MHD viscous flow over a shrinking sheet,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 7, pp. 2853–2857, 2009.

Muhaimin, R. Kandasamy, and A. B. Khamis, “Effects of heat and mass transfer on nonlinear MHD boundary layer flow over a shrinking sheet in the presence of suction,” Applied Mathematics and Mechanics, vol. 29, no. 10, pp. 1309–1317, 2008.

N. F. M. Noor, S. Awang Kechil, and I. Hashim, “Simple non-perturbative solution for MHD viscous flow due to a shrinking sheet,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 2, pp. 144–148, 2010.

N. F. Mohd and I. Hashim, “MHD flow and heat transfer adjacent to a permeable shrinking sheet embedded in a porous medium,” Sains Malaysiana, vol. 38, no. 4, pp. 559–565, 2009.

S. Nadeem and A. Hussain, “MHD flow of a viscous fluid on a nonlinear porous shrinking sheet with homotopy analysis method,” Applied Mathematics and Mechanics, vol. 30, no. 12, pp. 1569–1578, 2009.

Muhaimin, R. Kandasamy, I. Hashim, and A. B. Khamis, “On the effect of chemical reaction, heat and mass transfer on nonlinear MHD boundary layer past a porous shrinking sheet with suction,” Theoretical and Applied Mechanics, vol. 36, no. 2, pp. 101–117, 2009.

S. Abbasbandy, “A numerical solution of Blasius equation by Adomian's decomposition method and comparison with homotopy perturbation method,” Chaos, Solitons and Fractals, vol. 31, no. 1, pp. 257–260, 2007.

N. S. Elgazery, “Numerical solution for the Falkner-Skan equation,” Chaos, Solitons and Fractals, vol. 35, no. 4, pp. 738–746, 2008.

B. L. Kuo, “Heat transfer analysis for the Falkner-Skan wedge flow by the differential transformation method,” International Journal of Heat and Mass Transfer, vol. 48, no. 23-24, pp. 5036–5046, 2005.

A. Wazwaz, “The variational iteration method for solving two forms of Blasius equation on a half-infinite domain,” Applied Mathematics and Computation, vol. 188, no. 1, pp. 485–491, 2007.

S. J. Liao, Beyond Perturbation, Introduction to Homotopy Analysis Method, Chapman & Hall, Boca Raton, Fla, USA, 2003.

B. Yao and J. Chen, “A new analytical solution branch for the Blasius equation with a shrinking sheet,” Applied Mathematics and Computation, vol. 215, no. 3, pp. 1146–1153, 2009.

B. Yao and J. Chen, “Series solution to the Falkner-Skan equation with stretching boundary,” Applied Mathematics and Computation, vol. 208, no. 1, pp. 156–164, 2009.

A. A. Joneidi, G. Domairry, M. Babaelahi, and M. Mozaffari, “Analytical treatment on magnetohydrodynamic (MHD) flow and heat transfer due to a stretching hollow cylinder,” International Journal for Numerical Methods in Fluids, vol. 63, no. 5, pp. 548–563, 2010.

A. R. Noiey, N. Haghparast, M. Miansari, and D. D. Ganji, “Application of homotopy perturbation method to the MHD pipe flow of a fourth grade fluid,” Journal of Physics, vol. 96, no. 1, 2008.

M. Jalaal and D. D. Ganji, “An analytical study on motion of a sphere rolling down an inclined plane submerged in a Newtonian fluid,” Powder Technology, vol. 198, no. 1, pp. 82–92, 2010.

M. Jalaal, D. D. Ganji, and G. Ahmadi, “An analytical study on settling of non-spherical particles,” Asia-Pacific Journal of Chemical Engineering, 2010.

M. Jalaal, D. D. Ganji, and G. Ahmadi, “Analytical investigation on acceleration motion of a vertically falling spherical particle in incompressible Newtonian media,” Advanced Powder Technology, vol. 21, no. 3, pp. 298–304, 2010.

M. Jalaal and D. D. Ganji, “On unsteady rolling motion of spheres in inclined tubes filled with incompressible Newtonian fluids,” Advanced Powder Technology, vol. 22, pp. 58–67, 2011.

S. M. Moghimi, D. D. Ganji, H. Bararnia, M. Hosseini, and M. Jalaal, “Homotopy perturbation method for nonlinear MHD Jeffery-Hamel problem,” Computers & Mathematics with Applications, vol. 61, no. 8, pp. 2213–2216, 2011.

S. S. Motsa and S. Shateyi, “A new approach for the solution of three-dimensional magnetohydrodynamic rotating flow over a shrinking sheet,” Mathematical Problems in Engineering, vol. 2010, Article ID 586340, 15 pages, 2010.

S. S. Motsa, P. Sibanda, and S. Shateyi, “A new spectral-homotopy analysis method for solving a nonlinear second order BVP,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 9, pp. 2293–2302, 2010.

S. S. Motsa, P. Sibanda, F. G. Awad, and S. Shateyi, “A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem,” Computers & Fluids, vol. 39, no. 7, pp. 1219–1225, 2010.

Z. G. Makukula, P. Sibanda, and S. Motsa, “On a new solution for the viscoelastic squeezing flow between two parallel plates,” Journal of Advanced Research in Applied Mathematics, vol. 2, no. 4, pp. 31–38, 2010.

Z. G. Makukula, P. Sibanda, and S. S. Motsa, “A note on the solution of the von Kármán equations using series and Chebyshev spectral methods,” Boundary Value Problems, vol. 2010, Article ID 471793, 17 pages, 2010.

Zodwa G. Makukula, Precious Sibanda, and Sandile S. Motsa, “A novel numerical technique for two-dimensional laminar flow between two moving porous walls,” Mathematical Problems in Engineering, vol. 2010, Article ID 528956, 15 pages, 2010.

Z. G. Makukula, P. Sibanda, and S. Motsa, “On a new quasi-linearisation method for heat transfer in a visco-elastic fluid between parallel plates,” in New Aspects of Fluid Mechanics, Heat Transfer and Environment, N. Mastorakis, V. Mladenov, and Z. Bojkovic, Eds., 8th iasme/wseas international conference on fluid mechanics & Aerodynamics (FMA '10), pp. 178–186, Tapei, Taiwan, August 2010.

Z. G. Makukula, P. Sibanda, and S. S. Motsa, “On new solutions for heat transfer in a visco-elastic fluid between parallel plates,” International Journal of Mathematical Models and Methods in Applied Sciences, vol. 4, no. 4, pp. 221–230, 2010.

M. Sajid and T. Hayat, “The application of homotopy analysis method for MHD viscous flow due to a shrinking sheet,” Chaos, Solitons and Fractals, vol. 39, no. 3, pp. 1317–1323, 2009.

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

L. N. Trefethen, Spectral Methods in MATLAB, vol. 10, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2000.

H. S. Takhar, A. J. Chamkha, and G. Nath, “Natural convection on a thin vertical cylinder moving in a high-porosity ambient medium,” International Journal of Engineering Science, vol. 41, no. 16, pp. 1935–1950, 2003.