- Email: [email protected]
Application of homotopy perturbation method in nonlinear heat conduction and convection equations
Physics Letters A 360 (2007) 570–573 www.elsevier.com/locate/pla
Application of homotopy perturbation method in nonlinear heat conduction and convection equations A. Rajabi ∗ , D.D. Ganji, H. Taherian Department of Mechanical Engineering, Mazandaran University, Babol, PO Box 484, Iran Received 20 June 2006; received in revised form 13 August 2006; accepted 23 August 2006 Available online 12 September 2006 Communicated by A.R. Bishop
Abstract One of the newest analytical methods to solve nonlinear equations is the application of both homotopy and perturbation techniques. In this Letter, homotopy perturbation method (HPM), which does not need small parameters, is compared with the perturbation method in the field of heat transfer. The justification for using this method is the difficulties and limitations of perturbation and homotopy when used individually. In this research, homotopy perturbation method is used to solve for the temperature distribution in lumped system of combined convection–radiation and also a nonlinear equation of the steady conduction in a slab with variable thermal conductivity. © 2006 Elsevier B.V. All rights reserved. Keywords: Homotopy perturbation method; Conduction and convection heat transfer
1. Introduction Most scientific problems such as heat transfer are inherently of nonlinearity. We know that except a limited number of these problems, most of them do not have analytical solution. Therefore, these nonlinear equations should be solved using other methods. Some of them are solved using numerical techniques and some are solved using the analytical method of perturbation. In the numerical method, stability and convergence should be considered so as to avoid divergence or inappropriate results. In the analytical perturbation method, we should exert the small parameter in the equation. Therefore, finding the small parameter and exerting it into the equation are difficulties of this method. Since there are some limitations with the common perturbation method, and also because the basis of the common perturbation method is upon the existence of a small parameter, developing the method for different applications is very difficult. * Corresponding author.
E-mail address: [email protected] (A. Rajabi). 0375-9601/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2006.08.079
Therefore, many different methods have recently introduced some ways to eliminate the small parameter, such as artificial parameter method introduced by He [1,2], the homotopy analysis method by Liao [3,4], the variational iteration method by He [5–7]. One of the semi-exact methods is the homotopy perturbation method [8–15]. The applications of this method in different fields of nonlinear equations, integro-differential equations, Laplace transform, fluid mechanics and heat transfer have been studied by Cai [16], Cveticanin [17], El-Shahed [18], Abbasbandy [19], Siddiqui [20,21] and Ganji [22–24]. Thermal properties of most metallic materials used in engineering are not constant but, are dependent on temperature. In some cases such as cooling of gas turbine blades (Fig. 1), nuclear reactors, etc., a combination of different kinds of heat transfer (conduction, convection and radiation) is present; Therefore, in this Letter, two types of heat transfer equations are studied via HPM. In the first type, the nonlinear equation of the cooling of a lumped system by combined convection and radiation is solved through the two methods: homotopy perturbation method and the common perturbation method. In the second type, the steady conduction in a slab with the variable thermal conductivity is solved using HPM.
A. Rajabi et al. / Physics Letters A 360 (2007) 570–573
571
Nomenclature HPM k k2 L p PM t T
homotopy perturbation method thermal conductivity thermal conductivity at temperature T2 later heat length parameter of homotopy perturbation method time temperature
length dimensionless length
x X
Greek symbols constant small parameter dimensionless time dimensionless temperature
β ε τ θ
changing from u0 to u(r). We consider ν, as the following: ν = ν0 + pν1 + p 2 ν2 ,
(7)
and the best approximation for the solution is: u = lim ν = ν0 + ν1 + ν2 + · · · . p→1
(8)
The above convergence is discussed in [16,17]. 3. The application of HPM in heat transfer equations Fig. 1. Kinds of heat transfer in gas turbine blades.
Example 1. Consider the cooling of a lumped system by combined convection and radiation [25]:
2. Basic idea of homotopy perturbation method In this Letter, we apply the homotopy perturbation method to the discussed problems. To illustrate the basic ideas of the new method, we consider the following nonlinear differential equation, A(u) − f (r) = 0,
r ∈ Ω.
With boundary conditions of: ∂u = 0, r ∈ Γ. B u, ∂n
(1)
(2)
(3)
Where L stands for the linear and N for the nonlinear part. Homotopy perturbation structure is shown as the following equation: H (ν, p) = L(ν) − L(u0 ) + pL(u0 ) + p N (ν) − f (r) = 0. (4) Where, ν(r, p) : Ω × [0, 1] → R.
(5)
Obviously, using Eq. (4) we have: H (ν, 0) = L(ν) − L(u0 ) = 0, H (ν, 1) = A(ν) − f (r) = 0,
(9)
3.1. Homotopy perturbation method In Eq. (9), assuming θa = θs = 0 and separating the linear and nonlinear parts of the equation: dθ + θ + εθ 4 = 0. (10) dτ Now we apply homotopy perturbation method to Eq. (10). L(ν) − L(θ0 ) + pL(θ0 ) + p εν 4 = 0. (11)
Where A(u) is defined as follows: A(u) = L(u) + N(u).
dθ + (θ − θa ) + ε θ 4 − θs4 = 0. dτ
(6)
where p ∈ [0, 1] is an embedding parameter and u0 is the first approximation that satisfies the boundary condition. The process of changes in p from zero to unity is that of ν(r, p)
Where, L(ν) =
dν + ν, dτ
L(θ0 ) =
dθ0 + θ0 . dτ
(12)
0 Assuming dθ dτ = θ0 = 0 and substituting ν from Eq. (7) into Eq. (11) and some simplification and rearranging based on powers of p-terms, we have:
dν0 + ν0 = 0, dτ τ = 0, ν0 = 1, dν 1 + ν1 + εν04 = 0, p1 : dτ τ = 0, ν1 = 0, dν2 + ν2 + 4εν03 ν1 = 0, p2 : dτ τ = 0, θ2 = 0. p0 :
(13) (14) (15) (16) (17) (18)
572
A. Rajabi et al. / Physics Letters A 360 (2007) 570–573
Solving Eqs. (13), (15) and (17) considering appropriate initial conditions, we have: ν0 = e−τ , ε ν1 = e−4τ − e−τ , 3 2 ν2 = ε 2 e−τ − 2e−4τ + e−7τ . 9 So θ will be generally as follows: θ = Lim ν0 + pν1 + p 2 ν2
(19) (20) (21)
2 1 = e−τ + ε e−4τ − e−τ + ε 2 e−τ − 2e−4τ + e−7τ . (22) 3 9 3.2. Perturbation method (PM) For weak radiation regime (ε 1), let us assume a regular perturbation expansion and calculate the first three terms. Thus we assume: θ = θ0 + εθ1 + ε θ2 + · · · .
(23)
Substituting Eq. (23) into Eq. (10) and collecting the terms based on powers of ε as 0, 1, 2, . . . gives: dθ0 dτ τ = 0, dθ1 ε1 : dτ τ = 0, dθ2 ε2 : dτ τ = 0, ε0 :
+ θ0 = 0,
(24)
θ0 = 1,
(25)
+ θ1 + θ04 = 0, θ1 = 0, + θ2 + 4θ03 θ1 = 0, θ2 = 0.
(26) (27) (28)
1 −4τ − e−τ , e 3 −4τ − 2e + e−7τ .
X = 0,
θ = 1,
X = 1,
θ = 0.
2 1 θ = e−τ + ε e−4τ − e−τ + ε 2 e−τ − 2e−4τ + e−7τ . 3 9
x = 0,
T = T1 ,
x = L,
T = T2 .
(33)
(36)
(37)
(38) where L(ν) =
d 2ν , dX 2
L(θ0 ) =
d 2 θ0 . dX 2
(39)
2
Assuming ddXθ20 = 0 and substituting ν from Eq. (7) into Eq. (38) and some simplification and rearranging based on powers of pterms, we have: d 2 ν0 = 0, dX 2
(40)
ν0 = 1,
X = 1,
(41)
p1 :
(42)
ν0 = 0, d 2 ν0 d 2 ν1 dν0 2 + ε ν + = 0, 0 dX dX 2 dX 2
X = 0,
ν1 = 0,
X = 1,
ν1 = 0,
(31) Example 2. Consider one-dimensional conduction in a slab of thickness L made of a material with temperature dependent thermal conductivity k. Assuming that the two faces are maintained at uniform temperatures T1 and T2 , T1 > T2 , the governing equation and the boundary conditions are d dT (32) k = 0, dx dx
(35)
Now we apply homotopy perturbation method to Eq. (36). d 2θ dθ 2 L(ν) − L(θ0 ) + pL(θ0 ) + p εθ +ε = 0, dX dX 2
X = 0, (30)
X=
into Eqs. (32) and (33), we obtain: d dθ (1 + εθ ) = 0, dX dX
p0 :
θ1 =
2 θ2 = e−τ 9 The three term expansion in Eq. (23) now becomes:
T − T2 , T 1 − T2
x , L k 1 − k2 ε = β(T1 − T2 ) = k2
(29)
The solutions of Eqs. (24), (26) and (28) are: θ0 = e−τ ,
where k2 is the thermal conductivity at temperature T2 and β is a constant. Using Eq. (34) and introducing the dimensionless quantities of: θ=
p→1
2
Let the thermal conductivity vary linearly with temperature such that k = k2 1 + β(T − T2 ) , (34)
p2 :
d 2 ν1 d 2 ν0 d 2 ν2 dν0 dν1 + ε ν + ν + 2 = 0, 0 1 dX dX dX 2 dX 2 dX 2
X = 0,
ν2 = 0,
X = 1,
ν2 = 0.
(43) (44)
(45)
The solutions of Eqs. (40), (42) and (44) are: −ε 2 ε ν0 = −X + 1, ν1 = X + X, 2 2 2 1 ε ν2 = ε 2 − X 3 + X 2 − X. 2 2
(46)
A. Rajabi et al. / Physics Letters A 360 (2007) 570–573
So θ will be generally as follows:
θ = −X + 1 +
1 ε2 −ε 2 ε X + X + ε2 − X 3 + X 2 − X. 2 2 2 2
(47) The perturbation method leads to the same results as the HPM. 4. Conclusion Comparing the two methods of HPM and PM shows that they have got nearly the same answers. As shown in Eqs. (7) and (8), the homotopy perturbation method does not need a small parameter. Finally, it has been attempted to show the capabilities and wide-range applications of the homotopy perturbation method in comparison with the previous ones in solving heat transfer problems.
[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
References [1] J.H. He, Non-perturbative methods for strongly nonlinear problems, Berlin, dissertation.de-Verlag im Internet GmbH, 2006. [2] J.H. He, Int. J. Mod. Phys. B 20 (10) (2006) 1141. [3] S.J. Liao, Int. J. Non-Linear Mech. 303 (1995) 371.
[22] [23] [24] [25]
573
S.J. Liao, Eng. Anal. Bound. Elm. 202 (1997) 91. J.H. He, J. Comput. Math. Appl. Mech. Eng. 167 (1998) 57. J.H. He, Comput. Math. Appl. Mech. Eng. 167 (1998) 69. J.H. He, Int. J. Non-Linear Mech. 344 (1999) 699. J.H. He, J. Comput. Math. Appl. Mech. Eng. 17 (8) (1999) 257. J.H. He, Int. J. Non-Linear Mech. 351 (2000) 37. J.H. He, Phys. Lett. A 350 (1–2) (2006) 87. J.H. He, X.H. Wu, Chaos Solitons Fractals 29 (1) (2006) 108. J.H. He, Phys. Lett. A 347 (4–6) (2005) 228. J.H. He, Chaos Solitons Fractals 26 (3) (2005) 827. J.H. He, Chaos Solitons Fractals 26 (3) (2005) 695. J.H. He, Int. J. Non-Linear Sci. Numer. Simul. 6 (2) (2005) 207. X.C. Cai, W.Y. Wu, M.S. Li, Int. J. Non-Linear. Sci. Numer. Simul. 7 (1) (2006) 109. L. Cveticanin, Chaos Solitons Fractals 30 (5) (2006) 1221. M. El-Shahed, Int. J. Non-Linear Sci. Numer. Simul. 6 (2) (2005) 163. S. Abbasbandy, Chaos Solitons Fractals 30 (5) (2006) 1206. A.M. Siddiqui, R. Mahmood, Q.K. Ghori, Int. J. Non-Linear Sci. Numer. Simul. 7 (1) (2006) 7. A.M. Siddiqui, M. Ahmed, Q.K. Ghori, Int. J. Non-Linear Sci. Numer. Simul. 7 (1) (2006) 15. D.D. Ganji, M. Rafei, Phys. Lett. A 356 (2) (2006) 131. D.D. Ganji, A. Rajabi, Int. Commun. Heat Mass Transfer 33 (3) (2006) 391. D.D. Ganji, Phys. Lett. A 355 (4–5) (2006) 337. A. Aziz, G. Hamad, Int. J. Mech. Eng. Ed. 5 (1997) 167.
Application of homotopy perturbation method in nonlinear heat conduction and convection equations A. Rajabi ∗ , D.D. Ganji, H. Taherian Department of Mechanical Engineering, Mazandaran University, Babol, PO Box 484, Iran Received 20 June 2006; received in revised form 13 August 2006; accepted 23 August 2006 Available online 12 September 2006 Communicated by A.R. Bishop
Abstract One of the newest analytical methods to solve nonlinear equations is the application of both homotopy and perturbation techniques. In this Letter, homotopy perturbation method (HPM), which does not need small parameters, is compared with the perturbation method in the field of heat transfer. The justification for using this method is the difficulties and limitations of perturbation and homotopy when used individually. In this research, homotopy perturbation method is used to solve for the temperature distribution in lumped system of combined convection–radiation and also a nonlinear equation of the steady conduction in a slab with variable thermal conductivity. © 2006 Elsevier B.V. All rights reserved. Keywords: Homotopy perturbation method; Conduction and convection heat transfer
1. Introduction Most scientific problems such as heat transfer are inherently of nonlinearity. We know that except a limited number of these problems, most of them do not have analytical solution. Therefore, these nonlinear equations should be solved using other methods. Some of them are solved using numerical techniques and some are solved using the analytical method of perturbation. In the numerical method, stability and convergence should be considered so as to avoid divergence or inappropriate results. In the analytical perturbation method, we should exert the small parameter in the equation. Therefore, finding the small parameter and exerting it into the equation are difficulties of this method. Since there are some limitations with the common perturbation method, and also because the basis of the common perturbation method is upon the existence of a small parameter, developing the method for different applications is very difficult. * Corresponding author.
E-mail address: [email protected] (A. Rajabi). 0375-9601/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2006.08.079
Therefore, many different methods have recently introduced some ways to eliminate the small parameter, such as artificial parameter method introduced by He [1,2], the homotopy analysis method by Liao [3,4], the variational iteration method by He [5–7]. One of the semi-exact methods is the homotopy perturbation method [8–15]. The applications of this method in different fields of nonlinear equations, integro-differential equations, Laplace transform, fluid mechanics and heat transfer have been studied by Cai [16], Cveticanin [17], El-Shahed [18], Abbasbandy [19], Siddiqui [20,21] and Ganji [22–24]. Thermal properties of most metallic materials used in engineering are not constant but, are dependent on temperature. In some cases such as cooling of gas turbine blades (Fig. 1), nuclear reactors, etc., a combination of different kinds of heat transfer (conduction, convection and radiation) is present; Therefore, in this Letter, two types of heat transfer equations are studied via HPM. In the first type, the nonlinear equation of the cooling of a lumped system by combined convection and radiation is solved through the two methods: homotopy perturbation method and the common perturbation method. In the second type, the steady conduction in a slab with the variable thermal conductivity is solved using HPM.
A. Rajabi et al. / Physics Letters A 360 (2007) 570–573
571
Nomenclature HPM k k2 L p PM t T
homotopy perturbation method thermal conductivity thermal conductivity at temperature T2 later heat length parameter of homotopy perturbation method time temperature
length dimensionless length
x X
Greek symbols constant small parameter dimensionless time dimensionless temperature
β ε τ θ
changing from u0 to u(r). We consider ν, as the following: ν = ν0 + pν1 + p 2 ν2 ,
(7)
and the best approximation for the solution is: u = lim ν = ν0 + ν1 + ν2 + · · · . p→1
(8)
The above convergence is discussed in [16,17]. 3. The application of HPM in heat transfer equations Fig. 1. Kinds of heat transfer in gas turbine blades.
Example 1. Consider the cooling of a lumped system by combined convection and radiation [25]:
2. Basic idea of homotopy perturbation method In this Letter, we apply the homotopy perturbation method to the discussed problems. To illustrate the basic ideas of the new method, we consider the following nonlinear differential equation, A(u) − f (r) = 0,
r ∈ Ω.
With boundary conditions of: ∂u = 0, r ∈ Γ. B u, ∂n
(1)
(2)
(3)
Where L stands for the linear and N for the nonlinear part. Homotopy perturbation structure is shown as the following equation: H (ν, p) = L(ν) − L(u0 ) + pL(u0 ) + p N (ν) − f (r) = 0. (4) Where, ν(r, p) : Ω × [0, 1] → R.
(5)
Obviously, using Eq. (4) we have: H (ν, 0) = L(ν) − L(u0 ) = 0, H (ν, 1) = A(ν) − f (r) = 0,
(9)
3.1. Homotopy perturbation method In Eq. (9), assuming θa = θs = 0 and separating the linear and nonlinear parts of the equation: dθ + θ + εθ 4 = 0. (10) dτ Now we apply homotopy perturbation method to Eq. (10). L(ν) − L(θ0 ) + pL(θ0 ) + p εν 4 = 0. (11)
Where A(u) is defined as follows: A(u) = L(u) + N(u).
dθ + (θ − θa ) + ε θ 4 − θs4 = 0. dτ
(6)
where p ∈ [0, 1] is an embedding parameter and u0 is the first approximation that satisfies the boundary condition. The process of changes in p from zero to unity is that of ν(r, p)
Where, L(ν) =
dν + ν, dτ
L(θ0 ) =
dθ0 + θ0 . dτ
(12)
0 Assuming dθ dτ = θ0 = 0 and substituting ν from Eq. (7) into Eq. (11) and some simplification and rearranging based on powers of p-terms, we have:
dν0 + ν0 = 0, dτ τ = 0, ν0 = 1, dν 1 + ν1 + εν04 = 0, p1 : dτ τ = 0, ν1 = 0, dν2 + ν2 + 4εν03 ν1 = 0, p2 : dτ τ = 0, θ2 = 0. p0 :
(13) (14) (15) (16) (17) (18)
572
A. Rajabi et al. / Physics Letters A 360 (2007) 570–573
Solving Eqs. (13), (15) and (17) considering appropriate initial conditions, we have: ν0 = e−τ , ε ν1 = e−4τ − e−τ , 3 2 ν2 = ε 2 e−τ − 2e−4τ + e−7τ . 9 So θ will be generally as follows: θ = Lim ν0 + pν1 + p 2 ν2
(19) (20) (21)
2 1 = e−τ + ε e−4τ − e−τ + ε 2 e−τ − 2e−4τ + e−7τ . (22) 3 9 3.2. Perturbation method (PM) For weak radiation regime (ε 1), let us assume a regular perturbation expansion and calculate the first three terms. Thus we assume: θ = θ0 + εθ1 + ε θ2 + · · · .
(23)
Substituting Eq. (23) into Eq. (10) and collecting the terms based on powers of ε as 0, 1, 2, . . . gives: dθ0 dτ τ = 0, dθ1 ε1 : dτ τ = 0, dθ2 ε2 : dτ τ = 0, ε0 :
+ θ0 = 0,
(24)
θ0 = 1,
(25)
+ θ1 + θ04 = 0, θ1 = 0, + θ2 + 4θ03 θ1 = 0, θ2 = 0.
(26) (27) (28)
1 −4τ − e−τ , e 3 −4τ − 2e + e−7τ .
X = 0,
θ = 1,
X = 1,
θ = 0.
2 1 θ = e−τ + ε e−4τ − e−τ + ε 2 e−τ − 2e−4τ + e−7τ . 3 9
x = 0,
T = T1 ,
x = L,
T = T2 .
(33)
(36)
(37)
(38) where L(ν) =
d 2ν , dX 2
L(θ0 ) =
d 2 θ0 . dX 2
(39)
2
Assuming ddXθ20 = 0 and substituting ν from Eq. (7) into Eq. (38) and some simplification and rearranging based on powers of pterms, we have: d 2 ν0 = 0, dX 2
(40)
ν0 = 1,
X = 1,
(41)
p1 :
(42)
ν0 = 0, d 2 ν0 d 2 ν1 dν0 2 + ε ν + = 0, 0 dX dX 2 dX 2
X = 0,
ν1 = 0,
X = 1,
ν1 = 0,
(31) Example 2. Consider one-dimensional conduction in a slab of thickness L made of a material with temperature dependent thermal conductivity k. Assuming that the two faces are maintained at uniform temperatures T1 and T2 , T1 > T2 , the governing equation and the boundary conditions are d dT (32) k = 0, dx dx
(35)
Now we apply homotopy perturbation method to Eq. (36). d 2θ dθ 2 L(ν) − L(θ0 ) + pL(θ0 ) + p εθ +ε = 0, dX dX 2
X = 0, (30)
X=
into Eqs. (32) and (33), we obtain: d dθ (1 + εθ ) = 0, dX dX
p0 :
θ1 =
2 θ2 = e−τ 9 The three term expansion in Eq. (23) now becomes:
T − T2 , T 1 − T2
x , L k 1 − k2 ε = β(T1 − T2 ) = k2
(29)
The solutions of Eqs. (24), (26) and (28) are: θ0 = e−τ ,
where k2 is the thermal conductivity at temperature T2 and β is a constant. Using Eq. (34) and introducing the dimensionless quantities of: θ=
p→1
2
Let the thermal conductivity vary linearly with temperature such that k = k2 1 + β(T − T2 ) , (34)
p2 :
d 2 ν1 d 2 ν0 d 2 ν2 dν0 dν1 + ε ν + ν + 2 = 0, 0 1 dX dX dX 2 dX 2 dX 2
X = 0,
ν2 = 0,
X = 1,
ν2 = 0.
(43) (44)
(45)
The solutions of Eqs. (40), (42) and (44) are: −ε 2 ε ν0 = −X + 1, ν1 = X + X, 2 2 2 1 ε ν2 = ε 2 − X 3 + X 2 − X. 2 2
(46)
A. Rajabi et al. / Physics Letters A 360 (2007) 570–573
So θ will be generally as follows:
θ = −X + 1 +
1 ε2 −ε 2 ε X + X + ε2 − X 3 + X 2 − X. 2 2 2 2
(47) The perturbation method leads to the same results as the HPM. 4. Conclusion Comparing the two methods of HPM and PM shows that they have got nearly the same answers. As shown in Eqs. (7) and (8), the homotopy perturbation method does not need a small parameter. Finally, it has been attempted to show the capabilities and wide-range applications of the homotopy perturbation method in comparison with the previous ones in solving heat transfer problems.
[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
References [1] J.H. He, Non-perturbative methods for strongly nonlinear problems, Berlin, dissertation.de-Verlag im Internet GmbH, 2006. [2] J.H. He, Int. J. Mod. Phys. B 20 (10) (2006) 1141. [3] S.J. Liao, Int. J. Non-Linear Mech. 303 (1995) 371.
[22] [23] [24] [25]
573
S.J. Liao, Eng. Anal. Bound. Elm. 202 (1997) 91. J.H. He, J. Comput. Math. Appl. Mech. Eng. 167 (1998) 57. J.H. He, Comput. Math. Appl. Mech. Eng. 167 (1998) 69. J.H. He, Int. J. Non-Linear Mech. 344 (1999) 699. J.H. He, J. Comput. Math. Appl. Mech. Eng. 17 (8) (1999) 257. J.H. He, Int. J. Non-Linear Mech. 351 (2000) 37. J.H. He, Phys. Lett. A 350 (1–2) (2006) 87. J.H. He, X.H. Wu, Chaos Solitons Fractals 29 (1) (2006) 108. J.H. He, Phys. Lett. A 347 (4–6) (2005) 228. J.H. He, Chaos Solitons Fractals 26 (3) (2005) 827. J.H. He, Chaos Solitons Fractals 26 (3) (2005) 695. J.H. He, Int. J. Non-Linear Sci. Numer. Simul. 6 (2) (2005) 207. X.C. Cai, W.Y. Wu, M.S. Li, Int. J. Non-Linear. Sci. Numer. Simul. 7 (1) (2006) 109. L. Cveticanin, Chaos Solitons Fractals 30 (5) (2006) 1221. M. El-Shahed, Int. J. Non-Linear Sci. Numer. Simul. 6 (2) (2005) 163. S. Abbasbandy, Chaos Solitons Fractals 30 (5) (2006) 1206. A.M. Siddiqui, R. Mahmood, Q.K. Ghori, Int. J. Non-Linear Sci. Numer. Simul. 7 (1) (2006) 7. A.M. Siddiqui, M. Ahmed, Q.K. Ghori, Int. J. Non-Linear Sci. Numer. Simul. 7 (1) (2006) 15. D.D. Ganji, M. Rafei, Phys. Lett. A 356 (2) (2006) 131. D.D. Ganji, A. Rajabi, Int. Commun. Heat Mass Transfer 33 (3) (2006) 391. D.D. Ganji, Phys. Lett. A 355 (4–5) (2006) 337. A. Aziz, G. Hamad, Int. J. Mech. Eng. Ed. 5 (1997) 167.