# Explicit Solutions of a Gravity-Induced Film Flow along a Convectively Heated Vertical Wall

Authors: Ammarah Raees; Hang Xu
Journal: The Scientific World Journal
 (2013)

---

## Abstract

The gravity-driven film flow has been analyzed along a vertical
wall subjected to a convective boundary condition. The Boussinesq approximation is applied
to simplify the buoyancy term, and similarity transformations are used on the mathematical
model of the problem under consideration, to obtain a set of coupled ordinary differential
equations. Then the reduced equations are solved explicitly by using homotopy analysis
method (HAM). The resulting solutions are investigated for heat transfer effects on velocity
and temperature profiles.


---

## Body

## 1. Introduction

The flow domain, described by thin film flow, in one of the dimensions is much smaller than the flow in the other one or two dimensions. By utilizing this fact, a set of simplified modeling equations can possibly be deduced from the Navier Stokes equations. Gravity-driven laminar flow problems including the thin film flow have significant practical applications in many fields like industrial and chemical engineering, coating flows, biofluids, microfluidic engineering, and medicine. The vast majority of the investigations on different falling film flow phenomena have been studied over the years. Fulford [1] has described in detail a variety of concepts to analyze the thin film flow procedure. However, Andersson and Ytrehus [2, 3] studied the diffusion from a vertical wall into an accelerating falling liquid film and gave the Falkner-Skan solutions for gravity-driven film flows. On the other hand, Sparrow et al. [4] described the combined forced and free convection in boundary layer flow about the nonisothermal body subjected to a uniform free stream velocity, and also he gave the criteria for cataloging flows as purely forced, purely free, and mixed. A different approach was adopted by Andersson and Irgens [5, 6], namely, to divide the accelerating film flow into a developing viscous boundary layer and an external inviscid free stream. They furthermore demonstrated that a similarity transformation exists, such that the boundary layer momentum equation for power-law fluids is exactly transformed into a Falkner-Skan type ordinary differential equation. Andersson et al. [7] also investigated the effects of high Prandtl number and temperature differences in the laminar film flow with combined and natural convection. He also concluded that for the vast majority of the parameter combinations the resulting velocity profiles f′(η) increased monotonically from zero at the surface to unity in the free stream.Our motivation to do the present work is to investigate the heat transfer due to the gravity-driven laminar film flow over the convectively heated surface. Further, both aiding and opposing buoyancy are considered to see its effect on the film flow and heat transfer. The homotopy analysis method (HAM; see [8–13]) is implied to solve the considered problem, and explicit solutions with high precision are also obtained. To our knowledge, this is the first time to explore the explicit solutions for this particular gravity-driven film flow problem with convective boundary condition. Moreover, the squared residual has been calculated which shows the correctness of our obtained explicit solutions.


## 2. Mathematical Formulation and Analysis

Consider the two-dimensional laminar film flow of Newtonian fluid along a vertical surface. That vertical surface is heated or cooled from below by convection from a fluid of temperatureTf, while T∞ is the temperature outside the boundary layers. To describe this laminar film flow we must add a gravity term (body force term) to the momentum equation, and then using the boundary layer, assumptions will give rise to the following partial differential equation:
(1)∂u∂x+∂υ∂y=0,(2)ρ0[u∂u∂x+υ∂u∂y]=ρg+μ∂2u∂y2,(3)ρ0cp[u∂T∂x+υ∂T∂y]=k∂2T∂y2,
subject to boundary conditions
(4)u=0,υ=0,-k∂T∂y=hf(Tf-T),aty=0,u⟶U(x),T⟶T∞,asy⟶∞.
Here u is vertical velocity component along vertical x-axis and υ is horizontal velocity component along horizontal y-axis. Whereas cp is specific heat at constant pressure, ρ is the fluid density, μ is dynamic viscosity, k is thermal conductivity, g is the gravitational constant, hf(x) is the heat transfer coefficient due to Tf, and U(x) is free stream velocity. According to the classical boundary-layer approximations, it is assumed that stream-wise diffusion of momentum and longitudinal heat conduction is negligible in (2) and (3). Also all the properties of fluid are kept constant and the internal dissipation of energy is neglected except the density, (ρ(T)). For this purpose the Boussinesq approximation is used on the buoyancy term in the momentum equation to represent the density variation as
(5)ρ=ρ0[1-α(T-T0)],
where ρ0 is the density at the arbitrary reference temperature T0 and α is the coefficient of thermal expansion. Since the frictionless flow between the viscous boundary layer and the free streamline bordering the constant-pressure atmosphere is considered to be irrotational with downward velocity, U(x) is quasi-one-dimensional. Assuming the infinite film thickness at the entrance x=0, the simple free stream solution
(6)U(x)=2gx
is readily derived by integrating the one-dimensional version of (2), given as
(7)UdUdx=g.
We now introduce the following similarity transformations:
(8)ψ(η)=f(η)·(4Uνx3)1/2,η=y·(3U4νx)1/2,θ(η)=T-T∞Tf-T∞,
where ν is kinematic viscosity and ψ is the stream function defined as
(9)u=∂ψ∂y,υ=-∂ψ∂x.
The following set of equations is formed by substituting the similarity transformations from (8) to (1)–(3) and thus can be rewritten as
(10)f′′′(η)+f(η)f′′(η)+23(1-f′(η)2)-43λθ=0,θ′′(η)+Prf(η)θ′(η)=0
with the associated boundary conditions given as
(11)f(0)=0,f′(0)=0,f′(∞)=1,θ(∞)=0,θ′(0)=-γ(1-θ(0)),
where Pr is the Prandtl number, γ is reduced heat transfer parameter, and λ is the dimensionless temperature difference defined as λ≡Grx/Rex2=α(Tf-T∞)/2. Since the dimensionless parameters related to the measure of convection velocity are defined as the local Grashof number, Grx=αg(Tf-T0)x3/ν2, and local Reynolds number, Rex=Ux/ν.


## 3. Explicit Solutions by Homotopy Analysis Method

The homotopy analysis method has been employed here to give the explicit solutions of the nonlinear coupled differential equation. As it is shown by(A.11) (see the Appendix), the solution series for f(η) and θ(η) can be expressed as
(12)f(η)=f0(η)+∑m=1+∞‍fm(η),θ(η)=θ0(η)+∑m=1+∞θm(η),
where fm(η) and θm(η) are higher order deformation derivatives. Now we can solve the first several higher order deformation equations (A.12)–(A.14) to obtain fm(η) and θm(η), which is defined as
(13)fm(η)=A0,0m+∑j=12m+2∑i=02m+2‍σi,jmAi,jmηiexp-jβη,θm(η)=∑j=12m+2‍∑i=02m+2‍σi,jmBi,jmηiexp-jβη,
where β is the given positive constant and σi,jm is the coefficient defined as
(14)σi,jm={1,0≤i≤2m,1≤j≤2m+2-χi-2⌊i+12⌋,0,otherwise.
Here χi is given by (A.17) and ⌊·⌋ represents the integer floor function. Substituting the above expression into (A.12)–(A.14), we have the recursive coefficients Ai,jm and Bi,jm for m≥1 determined as
(15)Ai,jm=χmσi,jm-1Ai,jm-1+∑q=i2m+2‍Γq,jmμj,iq,Bi,jm=χmδi,jm-1Bi,jm-1+∑q=i2m+2‍Ψq,jmδj,iq,
for 2≤j≤2m+2 and 0≤i≤2m,
(16)Ai,1m=χmσi,1m-1Ai,1m-1+∑q=max[0,i-1]2m‍Δq,1mμ1,iq,Bi,1m=χmδi,1m-1Bi,1m-1+∑q=max[0,i-1]2m‍Ωq,1mδ1,iq,A0,0m=χmA0,0m-1+∑j=22m+2‍∑q=12m+2‍Γq,jm[(j-1)μj,0q-1βμj,1q]-∑q=02m‍1βΔq,1mμ1,1q+∑j=22m+2‍(j-1)Γ0,jmμj,00,A0,1m=χmσ0,1m-1A0,1m-1+∑j=22m+2‍∑q=12m+2‍1βΓq,jmμj,1q+∑q=02m‍1βΔq,1mμ1,1q-∑j=22m+2‍jΓq,jmμj,0q,B0,1m=χmσ0,1m-1B0,1m-1-∑j=22m+2‍∑q=02m‍Ψq,jmδj,0q,A0,jm=0,
where
(17)δn,jq={q!j!1(2β)q-j+2,0≤j≤q,n=1,12β(q+1),j=q+1,n=1,q!j!12β(1[(n+1)β]q-j+1-1[(n-1)β]q-j+1),0≤j≤q,n>1,0,j=q+1,n>1,μn,jq={∑i=jq+1‍δn,jqβi-j+1,0≤j≤q+1,n=1,∑i=jq‍δn,jq(nβ)i-j+1,0≤j≤q,n>1,0,j=q+1,n>1.
Other coefficients involved in the above recursive formulae are given as
(18)Γi,jm=ℏf(A~i,jm+B~i,jm+E~i,jm+Ji,jm-1-43λσi,jm-1Bi,jm-1),Δi,1m=ℏf(Ji,1m-1-43λσi,1m-1Bi,1m-1+B~i,1m),Ψi,jm=ℏθ(H~i,jm-1+PrF~i,jm+PrK~i,jm),Ωi,1m=ℏθ(Hi,1m-1+PrK~i,1m),
where ℏf and ℏθ are the convergence control parameters, and
(19)Di,jm=σi+1,jmAi+1,jm(i+1)-(jβ)σi,jmAi,jm,Ei,jm=σi+1,jmBi+1,jm(i+1)-(jβ)σi,jmBi,jm,Gi,jm=(i+2)(i+1)σi+2,jmAi+2,jm-(2jβ)(i+1)σi+1,jmAi+1,jm+(jβ)2σi,jmAi,jm,Hi,jm=(i+2)(i+1)σi+2,jmBi+2,jm-(2jβ)(i+1)σi+1,jmBi+1,jm+(jβ)2σi,jmBi,jm,Ji,jm=(i+3)(i+2)(i+1)σi+3,jmAi+3,jm-(3jβ)(i+2)(i+1)σi+2,jmAi+2,jm+3(i+1)(jβ)2σi+1,jmAi+1,jm-(jβ)3σi,jmAi,jm,Ki,jm=(i+3)(i+2)(i+1)σi+3,jmBi+3,jm-(3jβ)(i+2)(i+1)σi+2,jmBi+2,jm+3(i+1)(jβ)2σi+1,jmBi+1,jm-(jβ)3σi,jmBi,jm.
Also,
(20)A~i,jm=∑n=0m-1‍∑s=max[1,j-2n-2]min[2m-2n,j-1]‍∑r=max[1,i-2n-2]min[2m-2n,i]‍σr,sm-1-nAr,sm-1-nGi-r,j-sn,B~i,jm=∑n=max[j-1,i-1]m-1‍A0,0m-1-nGi,jn,E~i,jm=∑n=0m-1‍∑s=max[1,j-2n-2]min[2m-2n,j-1]‍∑r=max[1,i-2n-2]min[2m-2n,i]‍Dr,sm-1-nDi-r,j-sn,F~i,jm=∑n=0m-1‍∑s=max[1,j-2n-2]min[2m-2n,j-1]‍∑r=max[1,i-2n-2]min[2m-2n,i]‍σr,sm-1-nAr,sm-1-nEi-r,j-sn,K~i,jm=∑n=max[j-1,i-1]m-1‍A0,0m-1-nEi,jn.
Using all the above recursive formulae and setting A0,00=1, A0,10=-1/β, A0,20=1/β, B0,10=1, and B0,20=-β/(γ+2β), we can calculate all the coefficients, and thus purely explicit analytical solution for f(η) and θ(η) can be obtained for m=1,2,3,….


## 4. Results and Discussions

It is well known that convergence of the HAM series, given by(A.1), not only depends upon the proper choice of initial guess or linear operator but it also relies on the proper value of convergence-control parameter. Basically we can select the appropriate values of auxiliary parameters by two methods: first by plotting the ℏ-curves [10] and then choosing its value in the corresponding valid regions of ℏ and second by determining the minimum of the squared residual [9] of the governing equation which gives the appropriate value of convergence-control parameter. However, one can also accelerate the convergence of homotopy series by using the homotopy-Pade technique or homotopy-iterative approach [9] for the problems with strong nonlinearity.In order to check the validity of our gained explicit solutions, we apply the second approach by defining the discrete squared residual forf(η) and θ(η) as
(21)Errfm=1k+1∑j=0k‍[Nf(∑n=0m‍fn)]2dηErrθm=1k+1∑j=0k‍[Nf(∑n=0m‍θn)]2dη,
where k is an integer. Nf and Nθ denote the nonlinear operators given by (A.7). Here we use k=40 to gain the computational error. It is pertinent to mention here that using discrete squared residual reduces the CPU time, suggested in [9, 14] rather than finding exact square residual of the governing equations. Especially, it is convenient to use it when there is more than one unknown auxiliary parameter involved.Further, one can improve the convergence of obtained results by introducing more convergence-control parameters in the frame of HAM. In particular, we have used another convergence-control parameter induced in the linear operator and initial guess denoted byβ. The other two convergence-control parameters are ℏf and ℏθ. The optimal values for the β, ℏf, and ℏθ are determined by evaluating the minimum error of (21). In our case the value of convergence-control parameters is taken to be ℏf=-1/2 and ℏθ=-1/2. Also, here the selected value for third parameter is β=5, which has greatly accelerated the convergence of our series solution by decreasing the average square residual. The reduction in the error becomes quite slow as compared to our present case if we change the value of β, and even in some cases the results are not convergent. So, we can always reduce the error and can control the convergence of HAM analytical approximations by choosing the proper values of convergence-control parameters.Obviously, for given order of approximationm, if the value of Errfm and Errθm is smaller, the better is the approximation. To confirm this, we calculated the error for f(η) given in Tables 1 and 2 for different values of γ and λ. Similarly the error for θ(η) is shown in Tables 3 and 4. It is also noted that Errfm and Errθm decrease more quickly in case of favorable buoyancy, which approached 1.6×10-5 and 2.1×10-5 at 40th-order approximation for several values of γ. The reduction in the error for each case with the increasing number of iterations shows the reliability of our analytic results. Also the accuracy of these results can be improved far more by increasing the order of approximation.Table 1
Average squared residual errors forErrfm with ℏf=-1/2, ℏθ=-1/2, and β=5 in the case of Pr=1 and λ=1.
kth order
γ
=
1
γ
=
5
γ
=
10
γ
=
15
γ
=
20
5
2.2555
×
1
0
-
1
2.48030
×
1
0
-
1
2.56497
×
1
0
-
1
2.60104
×
1
0
-
1
2.62083
×
1
0
-
1
10
6.91281
×
1
0
-
2
8.82655
×
1
0
-
2
9.43782
×
1
0
-
2
9.66925
×
1
0
-
2
9.78824
×
1
0
-
2
20
5.043001
×
1
0
-
3
1.00880
×
1
0
-
2
1.12832
×
1
0
-
2
1.17021
×
1
0
-
2
1.19141
×
1
0
-
2
30
8.19099
×
1
0
-
4
2.58524
×
1
0
-
3
3.04192
×
1
0
-
3
3.21354
×
1
0
-
3
3.30340
×
1
0
-
3
40
1.69486
×
1
0
-
4
9.09989
×
1
0
-
4
1.14358
×
1
0
-
3
1.23538
×
1
0
-
3
1.28426
×
1
0
-
3Table 2
Average squared residual errors forErrfm with ℏf=-1/2, ℏθ=-1/2, and β=5 in the case of Pr=1 and λ=-1.
kth order
γ
=
1
γ
=
5
γ
=
10
γ
=
15
γ
=
20
5
1.11846
×
1
0
-
1
9.74131
×
1
0
-
2
9.26109
×
1
0
-
2
9.06737
×
1
0
-
2
8.96424
×
1
0
-
2
10
1.09109
×
1
0
-
2
6.67928
×
1
0
-
3
5.51535
×
1
0
-
3
5.09535
×
1
0
-
3
4.88484
×
1
0
-
3
20
3.02345
×
1
0
-
4
6.20106
×
1
0
-
4
6.00372
×
1
0
-
4
5.82802
×
1
0
-
4
5.71778
×
1
0
-
4
30
8.14679
×
1
0
-
5
1.45809
×
1
0
-
4
1.45809
×
1
0
-
4
1.35611
×
1
0
-
4
1.32648
×
1
0
-
4
40
9.69675
×
1
0
-
6
1.70644
×
1
0
-
5
1.70644
×
1
0
-
5
1.58480
×
1
0
-
5
1.54868
×
1
0
-
5Table 3
Average squared residual errors forErrθm with ℏf=-1/2, ℏθ=-1/2, and β=5 in the case of Pr=1 and λ=1.
k
th order
γ
=
1
γ
=
5
γ
=
10
γ
=
15
γ
=
20
5
5.31587
×
1
0
-
2
9.32074
×
1
0
-
2
1.19486
×
1
0
-
1
1.33022
×
1
0
-
1
1.41151
×
1
0
-
1
10
1.81487
×
1
0
-
2
4.38441
×
1
0
-
2
5.34957
×
1
0
-
2
5.68101
×
1
0
-
2
5.83709
×
1
0
-
2
20
2.81659
×
1
0
-
3
7.36339
×
1
0
-
3
8.33899
×
1
0
-
3
8.64952
×
1
0
-
3
8.80023
×
1
0
-
3
30
4.96569
×
1
0
-
4
1.62656
×
1
0
-
3
1.91199
×
1
0
-
3
2.01813
×
1
0
-
3
2.07356
×
1
0
-
3
40
9.44952
×
1
0
-
5
4.78105
×
1
0
-
4
6.01700
×
1
0
-
4
6.50358
×
1
0
-
4
6.76279
×
1
0
-
4Table 4
Average squared residual errors forErrθm with ℏf=-1/2, ℏθ=-1/2, and β=5 in the case of Pr=1 and λ=-1.
k
th order
γ
=
1
γ
=
5
γ
=
10
γ
=
15
γ
=
20
5
5.01248
×
1
0
-
2
8.58652
×
1
0
-
2
1.10077
×
1
0
-
1
1.22631
×
1
0
-
1
1.30186
×
1
0
-
1
10
1.37880
×
1
0
-
2
3.10149
×
1
0
-
2
3.73286
×
1
0
-
2
3.92840
×
1
0
-
2
4.01227
×
1
0
-
3
20
1.10969
×
1
0
-
3
2.08900
×
1
0
-
3
2.24156
×
1
0
-
3
2.25752
×
1
0
-
3
2.25533
×
1
0
-
3
30
1.00475
×
1
0
-
4
1.92774
×
1
0
-
4
1.86875
×
1
0
-
4
1.80412
×
1
0
-
4
1.76338
×
1
0
-
4
40
8.43581
×
1
0
-
6
2.22668
×
1
0
-
5
2.18639
×
1
0
-
5
2.11068
×
1
0
-
5
2.06069
×
1
0
-
5Moreover, the velocity and temperature profiles are plotted for favorable (λ=-1) and unfavorable (λ=+1) buoyancy with the variation of γ. Figure 1 shows that in the case of favorable buoyancy, as the value of γ increases the velocity component increases with the increase in η but reduces by unfavorable buoyancy as presented in Figure 2. Further the temperature profile observed from Figures 3 and 4 describes the effect of γ with the temperature differences (|λ|=1). From these figures we can examine that as η increases the temperature decreases, while the increase in the reduced heat transfer coefficient γ increases θ(η) for both aiding and opposing buoyancy.Figure 1
Velocity profilef′(η) for different values of γ where λ=-1 and Pr=1.Figure 2
Velocity profilef′(η) for different values of γ where λ=1 and Pr=1.Figure 3
Temperature profileθ(η) for different values of γ where λ=1 and Pr=1.Figure 4
Temperature profileθ(η) for different values of γ where λ=-1 and Pr=1.


## 5. Conclusions

The explicit solutions are obtained in this paper using HAM for the gravity-driven film flow over a vertical impermeable sheet with convective boundary condition. By choosing the appropriate value for the convergence-control parameters we obtain the error forf(η) and θ(η). The decrease in the error with the increase in the number of approximations depicts the validity of our analytic solutions. According to the authors’ view about such kind of analytical solutions, they have never been presented in the literature before. Finally, this is due to homotopy analysis method that we are able to give quite accurate explicit solutions by choosing the proper values of convergence-control parameters, base functions, initial approximations, and linear operators.



---

*Source: 475939-2013-12-25.xml*