# Solution of Boundary Layer Problems with Heat Transfer by Optimal Homotopy Asymptotic Method

Authors: H. Ullah; S. Islam; M. Idrees; M. Arif
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Application of Optimal Homotopy Asymptotic Method (OHAM), a new analytic approximate technique for treatment of Falkner-Skan equations with heat transfer, has been applied in this work. To see the efficiency of the method, we consider Falkner-Skan equations with heat transfer. It provides us with a convenient way to control the convergence of approximate solutions when it is compared with other methods of solution found in the literature as finite difference (N. S. Asaithambi, 1997) and shooting method (Cebeci and Keller, 1971). The obtained solutions show that OHAM is effective, simpler, easier, and explicit.


---

## Body

## 1. Introduction

Most of the problems in engineering sciences are nonlinear, particularly some of heat transfer problems. Limited analytic methods are presented for the solution of such problems in the literature. Therefore, the researcher’s profound attention is to hunt some new analytic methods for the solution of the problems. Most of the methods like Adomian Decomposition Method (ADM) [1], Variational Iteration Method (VIM) [2], Differential Transform Method (DTM) [3], Radial basis function [4] and Homotopy Perturbation Method (HPM) [5], are used for the solution of weakly nonlinear problems, and limited for strongly nonlinear problems. For the solution of the strongly nonlinear problems the perturbation methods were studied [6–8]. These methods comprise a small parameter which cannot be found easily. To overcome this issue, some new analytic methods such as Artificial Parameters Method [9], Homotopy Analysis Method (HAM) [10], and Homotopy Perturbation Method (HPM) [5] were introduced. These methods pooled the homotopy with the perturbation techniques. Recently, Marinca et al. introduced Optimal Homotopy Asymptotic Method (OHAM) [11–15] for the solution of nonlinear problems which made the perturbation methods independent of the assumption of small parameters.The Falkner-Skan equation has been considered in the last forty years due to its importance in the boundary layer theory. The boundary layer theory plays a vital role in the diverse area of engineering and scientific applications. The solution of the Falkner-Skan equation has been studied numerically first by Hartree [16]. Smith and Cebeci [17, 18] solved this equation by shooting method. Maksyn [19] solved the Falkner-Skan equation by analytic approximation. Asithambi [20–22] found its solution by finite differences, Liao [23] applied homotopy analysis to solve Falkner-Skan equation, and recently Vera [24] found its solution by Fourier series. An important case is the Blasius equation. This problem was solved by Rosales and Valencia [25] using Fourier series. Boyd [26] found the solution of Falkner-Skan equation by numerical method. An enormous amount of research work has been invested in the study of nonlinear boundary value problems [27–38]. In this paper, we will deal with the Falkner-Skan equations with heat transfer, a nonlinear boundary value problem [24] in different forms.The motivation of this paper is to enhance OHAM for the solution of Falkner-Skan equation with heat transfer. In [11–15], OHAM has been proved to be useful for obtaining an approximate solution of nonlinear boundary value problems. In this work, we have proved that OHAM is also useful and reliable for the solution of the Falkner-Skan equation with heat transfer, hence showing its validity and great potential for the solution of transient physical phenomenon in science and engineering.In the succeeding section, the basic idea of OHAM [11–15] is formulated for the solution of boundary value problems arising in heat transfer. In Section 3, the effectiveness of the enhanced formulation of OHAM for Falkner-Skan equation with heat transfer has been studied. Two special cases [24] of Falkner-Skan equation with heat transfer problems have been analyzed.


## 2. Basic Mathematical Theory of OHAM

Let us consider the following differential equation:(1)ℒ(u(t))+h(t)+𝒩(u(t))=0,
along with boundary conditions of the form
(2)ℬ(u,dudt)=0,
where ℒ is the linear operator, u(t) is an unknown function, h(t) is a known function, 𝒩(u(t)) is a nonlinear differential operator, and ℬ is a boundary operator.According to OHAM, one can construct an optimal homotopyϕ(t,q): Ω×[0,1]→ℜ which satisfies
(3)(1-q)[ℒ(ϕ(t,q))+h(t)]=H(q)[ℒ(ϕ(t,q))+h(t)+𝒩(ϕ(t,q))],(4)ℬ(ϕ(t,q),∂ϕ(t,q)∂t)=0,
where q∈[0,1] is an embedding parameter, ϕ(t,q) is an unknown function, and H(q) is a nonzero auxiliary function. The auxiliary function H(q) is nonzero for q≠0 and H(0)=0. Equation (3) is the structure of OHAM homotopy.It is defined that(5)q=0⟹ϕ(t,0)=u0(t),q=1⟹ϕ(t,1)=u(t),
respectively. Thus, as q varies from 0 to1, the solution ϕ(t,q) varies from u0(t) to u(t), where u0(t) is obtained from (1) and (2) for p=0 as follows:
(6)ℒ(u0(t))+h(t)=0,ℬ(u0,du0dt)=0.Next, we choose auxiliary functionH(q) in the form
(7)H(q)=qC1+q2C2+q3C3+⋯,
where C1,C2,C3,… are constants and can be found latter.To obtain an approximate solution, we expandϕ(t,q,Ci) by Taylor’s series about q in the following form:
(8)ϕ(t,q,C1,C2,…,Ci)=u0(t)+∑k=1∞uk(t,C1,C2,…,Ci)qk,i=1,2,….Now substituting (8) into (1) and (2) and equating the coefficient of like powers of q, we obtain the zeroth-order problem given by (6), the first- and second-order problems given by (9)–(11), respectively, and the general governing equations for uk(t) given by (11):
(9)ℒ(u1(t))=C1𝒩0(u0(t)),ℬ(u1,du1dt)=0,(10)ℒ(u2(t))-ℒ(u1(t))=C1𝒩0(u0(t))+C1[ℒ(u1(t))+𝒩1(u0(t),u1(t))],ℬ(u2,du2dt)=0,(11)ℒ(uk(t))-ℒ(uk-1(t))=Ck𝒩0(u0(t))+∑i=1k-1Ci[ℒ(uk-i(t))+𝒩k-i(u0(t),u1(t),…,uk-i(t))],ℬ(uk,dukdt)=0,k=2,3,…,
where 𝒩k-i(u0(t),u1(t),…,uk-i(t)) is the coefficient of qk-i in the expansion series of 𝒩(ϕ(t,q)) about the embedding parameter p as follows:
(12)𝒩(ϕ(t,q,Ci))=𝒩0(u0(t))+∑k≥1𝒩k(u0,u1,u2,…,uk)qk,i=1,2,3,….It should be underscored that theuk for k≥0 are governed by the linear equations with linear boundary conditions that come from the original problem, which can be easily solved.It has been observed that the convergence of the series (8) depends upon the auxiliary constants C1,C2,…. If it is convergent at q=1, one has
(13)u~(t,C1,C2,…,Ci)=u0(t)+∑k≥1uk(t,C1,C2,…,Ci).Substituting (13) into (1), it results in the following expression for residual:
(14)R(t,C1,C2,…,Ci)=ℒ(u~(t,C1,C2,…,Ci))+h(t)+𝒩(u~(t,C1,C2,…,Ci)).IfR(t,C1,C2,…,Ci)=0, then u~(t,C1,C2,…,Ci) is the exact solution of the problem. Generally it does not happen, especially in nonlinear problems.For the determinations of auxiliary constants,Ci,i=1,2,…,m, there are different methods like Galerkin’s method, Ritz method, least squares method, and collocation method. One can apply the method of  least squares as follows:
(15)J(C1,C2,…,Cm)=∫abR2(t,C1,C2,C3…,Cm)dt,
where a and b are two values, depending on the nature of the given problem.The auxiliary constantsCi,i=1,2,…,m can be optimally found from
(16)∂J∂C1=∂J∂C2=⋯=∂J∂Cm=0.Themth order approximate solution can be obtained bythese constants. The constants Ci can also be determined by another method as follows:
(17)R(k1,C1,C2,…,Cm)=R(k2,C1,C2,…,Cm)=⋯=R(km,C1,C2,…,Cm)=0,i=1,2,…,m.The convergence of OHAM is directly proportional to the number of optimal constantsC1,C2,… which is determined by (16).It is easy to observe [13] that the Homotopy Perturbation Method (HPM) proposed by He [4] is a special case of (3) when H(q)=-q, and on the other hand, the Homotopy Analysis Method (HAM) proposed by Liao [11] is another special case of (3) when H(q)=qℏ where ℏ is chosen from “ℏ-curves” [12].


## 3.   Application of OHAM to Falkner-Skan Equations with Heat Transfer

To demonstrate the effectiveness of OHAM formulation, two models are studied.Model 1 (see [23]).
When an incompressible fluid passes in the vicinity of solid boundaries, the Navier-Stokes equations may be reduced drastically into the boundary layer equations:(18)u∂u∂x+v∂v∂y=UedUedx+v∂2u∂2y,∂u∂x+∂v∂y=0,
where Ue(x) is the free stream velocity, u and v are velocity components in x- and y-directions, and υ is the kinematic viscosity. In case of two-dimensional flow, the incompressible boundary layer flow over a wedge, when the free stream velocity is of the form Ue(x)=Kxm, is the following similarity transformation:
(19)u(x,y)=Ue(x)f′(η),η=y(m+1)K2x(m-1)/2,υ(η)=T-T∞Tw-T∞.
Using (19) into (18), we obtained the Falkner-Skan equation
(20)d3fdη3+fd2fdη2+β(1-(dfdη)2)=0,
along with boundary conditions
(21)f(0)=0,f′(0)=0,f′(5)=1.
According to (1), we have
(22)ℒ(f(η))=f′′′(η)+β,h(η)=0,𝒩(f(η))=f(η)f′′(η)-β(f′(η))2.
The boundary conditions are(23)f(0)=0,f′(0)=0,f′(5)=1.
Applying the method formulation mentioned in Section2 leads to the following.
Zeroth-Order Problem.Consider
(24)f0′′′(η)+β=0,(25)f0(0)=0,f0′(0)=0,f0′(5)=1,
from which we obtain
(26)f0(η)=160(6+75β)η2-1060βη3.
First-Order Problem.Consider
(27)f1′′′(η)=β+C1(β-β(f0′)2+f0f0′′+f0′′′)+f0′′′,(28)f1(0)=0,f1′(0)=0,f1′(5)=0.
Its solution is(29)f1(η)=(-5C148-5βC116+1825β2C1576+625β3C196)η2+(C13000+23βC13000+17β2C1480-5β3C148)η5+(-βC1900-11β2C1900+β3C148)η6+(β2C11260-β3C1840)η7.
Second-Order Problem.Consider
(30)f2′′′(η)=(1+C1)f1′′′+C2(β+f0f0′′-βf0′+f0′′′)+C1(f1f0′′-2βf0′f1′+f0f1′′),
with BC
(31)f2(0)=0,f2′(0)=0,f2′(5)=0,
whose solution is
(32)f2(η)=(-5C148-5βC116+1825β2C1576+625β3C196-95C124032-395βC12504-538225β2C12145152+2093125β3C12145152+3971875β4C1241472+4296875β5C1248384-5C248-5βC216+1825β2C2576+625β3C296)η2+(C13000+23βC13000+17β2C1480-5β3C148-13C1236000-41βC1224000+899β2C1217280+491β3C122304-125β4C121728-625576β5C12+C23000+23βC23000+17β2C2480-5β3C248)η5+(-βC1900-11β2C1900+β3C148+βC1221600-151β2C1214400-203β3C1210368-5256β4C12+125β5C121152-βC2900-11β2C2900+β3C248)η6+(β2C11260-β3C1840+β2C121260-1840β3C12+β2C21260-β3C2840)η7+(11C125040000+761βC1210080000+3181β2C124032000+463β3C12322560-85β4C128064+125β5C1216128)η8+(-βC12100800-253β2C121134000-17089β3C1218144000+2561β4C12725760-55β5C1224192)η9+(β2C1270875+β3C127000-227β4C12567000+β5C124320)η10+(-β3C12155925+19β4C121247400-β5C12118800)η11.
Adding (25), (27), and (29), we obtain
(33)f(η)=A2η2+A3η3+A5η5+A6η6+A7η7+A8η8+A9η9+A10η10+A11η11,
where
(34)A2=(110+5β4-5C124-5βC18+1825β2C1288+625β3C148-95C124032-395βC12504-538225β2C12145152+2093125β3C12145152+3971875β4C1241472+4296875β5C1248384-5C248-5βC216+1825β2C2576+625β3C296),A3=(-β6),A5=(C11500+23βC11500+17β2C1240-5β3C124-13C1236000-41βC1224000+899β2C1217280+491β3C122304-125β4C121728-625576β5C12+C23000+23βC23000+17β2C2480-5β3C248),A6=(-βC1450-11β2C1450+β3C124+βC1221600-151β2C1214400-203β3C1210368-5256β4C12+125β5C121152-βC2900-11β2C2900+β3C248),A7=(β2C1630-β3C1420+β2C121260-1840β3C12+β2C21260-β3C2840),A8=(11C125040000+761βC1210080000+3181β2C124032000+463β3C12322560-85β4C128064+125β5C1216128),A9=(-βC12100800-253β2C121134000-17089β3C1218144000+2561β4C12725760-55β5C1224192),A10=(β2C1270875+β3C127000-227β4C12567000+β5C124320),A11=(-β3C12155925+19β4C121247400-β5C12118800).
For the computation of the constantsC1 and C2 applying the method of least square mentioned in (14)–(16), we get
(35)C1=-0.08144083042145557,C2=-0.05004674822211045,forβ=1.
Putting these values in (30), we obtained the approximate solution of the form
(36)f(η)=6.67768×10-1η2-1.677×10-1η3+7.00357×10-3η5-1.20593×10-4η6+8.18634×10-5η7-3.23784×10-6η8+5.32977×10-7η9-7.89596×10-8η10+2.65857×10-9η11.
Now consider the energy equation of an incompressible fluid which passes through the vicinity of the solid boundaries:(37)u∂T∂x+v∂T∂y=a∂2T∂y2.
Upon using the transformation in (19) into (37), we obtained
(38)d2υdη2+Prfdυdη=0,
with boundary conditions
(39)υ(0)=1,υ(5)=0,
where “a” is the thermal diffusivity and Pr is the prandtl number.
Applying the method formulated in Section2 leads to the following.
Zeroth-Order Problem.Consider
(40)υ0′′(η)=0,υ(0)=1,υ(5)=0.
Its solution is(41)υ0(η)=5-η5.
First-Order Problem.Consider
(42)υ1′′(η)=PrC1f0(η)υ0′(η)+(1+C1)υ0′′(η),υ(0)=0,υ(5)=0,
whose solution is(43)
υ
1
(
η
)
=
(
5
P
r
C
1
24
+
25
16
P
r
β
C
1
)
η
-
(
P
r
C
1
600
-
1
48
P
r
β
C
1
)
η
4
+
1
600
P
r
β
C
1
η
5
.Second-Order Problem.Consider
(44)υ2′′(η)=Pr(C2f0(η)+C1f1(η))υ0(η)+(1+C1)υ1′′(η)+PrC1f0(η)υ1′(η)+C2υ0′′(η),υ2(0)=0,υ2(5)=0.We obtain the following solution:(45)υ2(η)=(5PrC124+2516PrβC1+65PrC124032+125P2rC124032+7564PrβC12+625P2rβC122688+970625Prβ2C12145152+3125P2rβ2C1218144+15625Prβ3C121512+5PrC224+2516PrβC2)η+(-PrC1600-148PrβC1+PrC1214400+1576P2rC12-164PrβC12+5144P2rβC12-365Prβ2C126912+125768P2rβ2C12-125Prβ3C121152-PrC2600-148PrβC2)η4+(1600PrβC1+1600PrβC12-1576P2rβC12-5384P2rβ2C12+1600PrβC2)η5+(-PrC12630000-P2rC1263000-23PrβC12630000-P2rβC122520-17Prβ2C12100800-5P2rβ2C122016+Prβ3C122016)η7+(PrβC12252000+P2rβC1228800+11Prβ2C12252000+P2rβ2C122304-Prβ3C1213440)η8+(-Prβ2C12453600-P2rβ2C1251840+Prβ3C12302400)η9.Adding (41), (43), and (45), we obtain
(46)υ(η)=1+B1η+B4η4+B5η5+B7η7+B8η8+B9η9,
where
(47)B1=(-15+5PrC112+258PrβC1+65PrC124032+125P2rC124032+7564PrβC12+625P2rβC122688+970625Prβ2C12145152+3125P2rβ2C1218144+15625Prβ3C121512+5PrC224+2516PrβC2),B4=(-PrC1300-124PrβC1+PrC1214400+1576P2rC12-164PrβC12+5144P2rβC12-365Prβ2C126912+125768P2rβ2C12-125Prβ3C121152-PrC2600-148PrβC2),B5=(1300PrβC1+1600PrβC12-1576P2rβC12-5384P2rβ2C12+1600PrβC2),B7=(-PrC12630000-P2rC1263000-23PrβC12630000-P2rβC122520-17Prβ2C12100800-5P2rβ2C122016+Prβ3C122016),B8=(PrβC12252000+P2rβC1228800+11Prβ2C12252000+P2rβ2C122304-Prβ3C1213440),B9=(-Prβ2C12453600-P2rβ2C1251840+Prβ3C12302400).Using (46) in (38) and applying the method of least square, we obtain
(48)C1=-0.027238050474420454,C2=-0.016903606167477613.Substituting these values in (46) for β=1 and Pr=10, we obtain
(49)υ(η)=1-1.29659η+2.95285×10-2η4-2.27213×10-3η5-2.12478×10-4η5+3.45784×10-6η8-1.42298×10-6η9.Model 2 (see [23]).
In case ofβ=0 in (52), we obtain the Blasius equation which is the famous equation of fluid dynamics and represent the problem of an incompressible fluid that passes through on a semi-infinite flat plate.
One has(50)d3fdη3+fd2fdη2=0,(51)f(0)=0,f′(0)=0,f′(5)=1.
According to (1), we define the operators
(52)ℒ(f(η))=f′′′(η),g(η)=0,𝒩(f(η))=f(η)f′′(η),
where f′′(η) and f′′′(η) represent the second and third derivatives of f(η) with respect to η.Applying the method formulated in Section2 leads to the following.Zeroth-Order Problem.Consider
(53)f0′′′(η)=0,(54)f0(0)=0,f0′(0)=0,f0′(5)=1.Its solution is(55)f0(η)=η210.First-Order Problem. Consider
(56)f1′′′(η)=-C1f0(η)f0′′(η)-f0′′′(η)-C1f0′′′(η),f1(0)=0,f1′(0)=0,f1′(5)=0,whose solution is(57)f1(η)=-5C148η2+C13000η5.Second-Order Problem.Consider
(58)f2′′′(η)=-C2f0(f0)′′-C1f1(η)f0′′(η)-C1f0(η)f1′′(η)-C2f0′′′(η)-(1+C1)f1′′′(η),f2(0)=0,f2′(0)=0,f2′(5)=0.
We obtain the following solution:
(59)f2(η)=(-5C148-95C124032-5C248)η2+(C13000-13C1236000+C23000)η5+11C125040000η8.Third-Order Problem. Consider
(60)f3′′′(η)=-C3f0(η)f0′′(η)-C2f1(η)f0′′(η)-C1f2(η)f0′′(η)-C2f0(η)f1′′(η)-C1f1(η)f1′′(η)-C1f0(η)f2′′(η)-C3f0′′′(η)-C2f1′′′(η)-(1+C1)f2′′′(η),f3(0)=0,f3′(0)=0,f3′(5)=0.
Its solution is
(61)f3(η)=(-5C148-95C122016+7915C13193536-5C248-95C1C22016-5C348)η2+(C13000-13C1218000-631C134032000+C23000-13C1C218000+C33000)η5+(11C122520000-11C134480000+11C1C22520000)η8+C1366528000η9.
From (53), (57), (59), and (61), we obtain
(62)f(η)=D2η2+D5η5+D9η9+D11η11,
where
(63)D2=110-5C116-95C121344+7915C13193536-5C224-95C1C22016-5C348,D5=C11000-13C1212000-631C134032000+C21500-13C1C218000+C33000,D9=11C121680000-11C134480000+11C1C22520000,D11=C1366528000.
Using (62) in (50) and applying the technique as discussed in (14)–(16), we obtain
(64)C1=-0.7047085843229471,C2=0.4265066662188469,C3=0.0092786353498473.
Substituting these values in (62), we have
(65)f(η)=1.95147×10-1η2-6.83434×10-4η5+2.79896×10-6η8-5.26047×10-9η11.


## 4. Results and Discussions

The formulation presented in Section2 provides highly accurate solutions for the problems demonstrated in Section 3. We have used Mathematica 7 for most of our computational work. In Table 1, we have presented the initial slope f′′(0) for different values of β obtained by OHAM. Table 1 shows that the results obtained by OHAM are in excellent agreement with the results found in the literature. Table 1 shows a benchmark for the initial slope f′′(0) found by different authors. It is found that the method present in this work is very good and provides the same values that optimized numerical methods with which it is compared. Figure 1 shows the variation of the function f′(η) against η for different values of β with OHAM, while Figure 2 shows the variation of the function υ′(η) with respect to η for different values of β. Table 2 shows the solution of the Blasius equation obtained by the present method. In order to verify the accuracy of the present method, we have compared the results obtained by OHAM to the results available in the literature and found an excellent agreement. Figures 2 and 3 show the variation of f′(η) with respect to η for the Blasius equation which is identical to results in the literature [41].Table 1
Initial slopef′′(0) obtained by OHAM for different values of β.
β
Hartree [16]
Asaithambi [22]
Asaithambi [21]
Salama [39]
Zhang [40]
Vera [24]
OHAM
2
1.687
1.687222
1.687218
1.687218
1.687218
1.687218
1.67054
1
1.233
1.23589
1.232588
1.232588
1.232587
1.232587
1.23567
0.5
0.927
0.927682
0.927680
0.927680
0.927680
0.927680
0.92697
−0.1
0.319
0.319270
0.319270
0.319270
0.319270
0.319270
0.31927Table 2
Solution of the Blasius equation by OHAM.
η
f
(
η
)
f
′
(
η
)
f
′
′
(
η
)
0.0
0.000000000
0.000000000
3.90295
×
10
-
1
0.5
4.87655
×
10
-
2
1.94934
×
10
-
1
3.88588
×
10
-
1
1.0
1.94467
×
10
-
1
3.86900
×
10
-
1
3.76782
×
10
-
1
1.5
4.33963
×
10
-
1
5.68522
×
10
-
1
3.45926
×
10
-
1
2.0
7.59425
×
10
-
1
7.28781
×
10
-
1
2.90680
×
10
-
1
2.5
1.15707
×
10
0
8.55368
×
10
-
1
2.12781
×
10
-
1
3.0
1.60768
×
10
0
9.39646
×
10
-
1
1.24115
×
10
-
1
3.5
2.08955
×
10
0
9.81346
×
10
-
1
4.67757
×
10
-
2
4.0
2.58389
×
10
0
9.92571
×
10
-
1
5.82245
×
10
-
3
4.5
3.08065
×
10
0
9.94745
×
10
-
1
8.42702
×
10
-
3
5.0
3.57943
×
10
0
1.0000000000
6.18202
×
10
-
4Figure 1
f
′
(
η
) for different values of β.Figure 2
υ
′
(
η
) for different values of β.Figure 3
Plot off′(η) with respect to η.


## 5. Conclusion

In this work, we have seen the effectiveness of OHAM [11–15] to Falkner-Skan, Energy and Blasius equations. By applying the basic idea of OHAM to Falken-Skan, Energy and Blasius equations, we found that it is simpler in applicability and, more convenient to control convergence and involved less computational overhead. Therefore, OHAM shows its validity and great potential for the solution Falken-Skan, Energy and Blasius equations, with heat transfer problems arising in science and engineering.



---

*Source: 324869-2013-09-17.xml*