# Effects of Thermocapillarity and Thermal Radiation on Flow and Heat Transfer in a Thin Liquid Film on an Unsteady Stretching Sheet

Authors: R. C. Aziz; I. Hashim; S. Abbasbandy
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper examines the effects of thermocapillarity and thermal radiation on the boundary layer flow and heat transfer in a thin film on an unsteady stretching sheet with nonuniform heat source/sink. The governing
partial differential equations are converted into ordinary differential equations by a similarity transformation and then are solved by using the homotopy analysis method (HAM). The effects of the radiation parameter, the
thermocapillarity number, and the temperature-dependent parameter in this
study are discussed and presented graphically via velocity and temperature
profiles.


---

## Body

## 1. Introduction

The analysis of heat transfer of boundary layer flow with thermal radiation is important in electrical power generation, astrophysical flows, solar power technology, space vehicle reentry, and other industrial engineering processes.Wang [1] was the pioneer in investigating the hydrodynamics of a flow in a thin liquid film on an unsteady stretching surface. Later, Andersson et al. [2] studied the heat transfer characteristics of the hydrodynamical problem solved by Wang [1]. Liu and Andersson [3] examined the problem with a more general form of prescribed temperature variation of a stretching sheet. Wang [4] investigated the same problem of Andersson et al. [2], presenting analytic solutions. Several researchers have extended Wang’s [1] classical problem, taking into consideration non-Newtonian [5–9], thermocapillarity [10–13], and magnetic effects [14, 15].Wang [4] was probably the first to analyze the flow and heat transfer in a thin liquid film on an unsteady stretching surface using the homotopy analysis method (HAM). The HAM is a general analytic method for obtaining series solutions for various types of nonlinear equations [16]. Several studies have successfully applied HAM to various nonlinear problems in science and engineering [17–22]. Studies by Dandapat et al. [10, 11] reported that thermocapillarity generates surface-tension gradients along the horizontal interface between the passive gas and the liquid film. Both Abd El-Aziz [23] and Mahmoud [24] have studied the influence of thermal radiation on the flow and heat transfer on an unsteady stretching sheet.Motivated by these studies, in this paper we extend the results of Wang [4] and Mahmoud [24] for thin film flow on an unsteady stretching sheet with combined effects of thermocapillarity, thermal radiation, and internal heating. Furthermore, we employ the similarity transformation introduced by Wang [4] to transform the extent of the independent variable into a finite range 0-1. The solutions reached using HAM are presented and implications discussed.


## 2. Problem Formulation

The fluid flow, modeled as an unsteady, two-dimensional, incompressible viscous laminar flow on a horizontal thin elastic sheet, emerges from a narrow slot at the origin of a Cartesian coordinate system. Fluid motion and heat transfer arise in the stretching of the horizontal elastic sheet. It is assumed that the elastic sheet has internal heat generation/absorption and that the sheet temperature varies with the coordinatex and time t. Under these assumptions, the governing conservation equations of mass, momentum, and energy at unsteady state can be expressed as(2.1)∂u∂x+∂v∂y=0,(2.2)∂u∂t+u∂u∂x+v∂u∂y=ν∂2u∂y2,(2.3)ρCp(∂T∂t+u∂T∂x+v∂T∂y)=κ∂2T∂y2+Q-∂qr∂y,
subject to(2.4)u=U,v=0,T=Tsaty=0,μ∂u∂y=∂σ∂x,∂T∂y=0,v=dhdtaty=h,
where u and v are the velocity components of the fluid in the x- and y-directions, t is the time, T is the temperature, ν is the kinematic viscosity, ρ is the density, κ is the thermal diffusivity, Cp is the specific heat at constant pressure, μ is the viscosity, and h(t) is the uniform thickness of the liquid film. Note that the surface tension varies linearly with temperature, defined as σ=σ0[1-δ(T-T0)], where δ is a positive fluid property. The velocity of the stretching surface is defined as U=bx/(1-αt), with α and b as positive constants. The Q term is the heat generated (>0) or absorbed (<0) per unit volume and is defined as(2.5)Q=(κUxν)B*(T-To),
where B* is the temperature-dependent heat generation/absorption, that is, B* is positive in the case of the elastic sheet’s generation of heat and negative in the case of the sheet’s absorption of heat from the fluid flow. The radiative heat flux qr under Rosseland approximation (Brewster, 1992, cited in [23]) is given in the form(2.6)qr=-4σ*3k1∂T4∂y,
where σ* is the Stefan Boltzmann constant and k1 is the mean absorption coefficient. We assume that the temperature difference within the flow is sufficiently small for T4 to be expressed as a linear function of temperature. This is accomplished by expanding T4 in a Taylor series about T0 and neglecting higher-order terms, thus(2.7)T4≅4T03T-3T04.
According to Wang [4], the surface temperature Ts of the stretching sheet varies with the distance x from the slot and time t in the form(2.8)Ts=To-Trefbx22ν(1-αt)-3/2,
where To is the temperature at the slit, Tref is the constant reference temperature for all t<1/α. The surface of the planar liquid film is assumed to be smooth and free of surface waves while viscous shear stress and heat flux are assumed to vanish at the adiabatic free surface. Liu and Andersson [3] developed similarity transformations, which are given as(2.9)ψ=βx[νb1-αt]1/2f(η),T=To-Tref[dx2ν(1-αt)3/2]θ(η),η=1β[bν(1-αt)]1/2y,
where β is the dimensionless film thickness and ψ(x,y) is the stream function defined by(2.10)u=∂ψ∂y=bx1-αtf′(η),v=-∂ψ∂x=-(νb1-αt)1/2βf(η),
where a prime denotes differentiation with respect to η. Apparently, (2.10) have already satisfied (2.1). Consequently, (2.1)–(2.4) are transformed to the following nonlinear boundary value problem:(2.11)f′′′+γ(ff″-12Sηf″-(f′)2-Sf′)=0,(2.12)1Pr(1+NR)θ″+γ(fθ′-2f′θ-12Sηθ′-32Sθ+1PrB*θ)=0,
subject to(2.13)f(0)=0,f′(0)=1,θ(0)=1,(2.14)f(1)=12S,f″(1)=Mθ(1),θ′(1)=0,
where a prime denoting differentiation with respect to η, S=α/2 is the dimensionless measure of unsteadiness, Pr is the Prandtl number, NR is the radiation parameter defined as NR=16σ*T03/κk1, and γ=β2 is an unknown constant which must be determined as a part of the present problem. The thermocapillarity number M is given as(2.15)M=βδσ0Trefμbν.
It is worth noting that the momentum boundary layer problem defined by (2.11) subject to the relevant boundary conditions (2.13) is decoupled from the thermal boundary layer problem, while the temperature field θ(η) is on the other hand coupled to the velocity field. The most important characteristics of flow and heat transfer are shear stress τs and heat flux qs of the stretching sheet, which are defined as(2.16)τs=μ(∂u∂y)y=0,qs=-κ(∂T∂y)y=0,
where μ is the fluid dynamic viscosity. The local skin-friction coefficient Cf and the local Nusselt number Nux can be defined as(2.17)Cf=2τsρu2,Nux=-xqsκTref.
Thus, the skin friction and the rate of heat transfer for fluid flow in a thin film can be expressed as(2.18)12CfRex1/2=1βf″(0),NuxRex-1/2=dx2βν(1-αt)1/2θ′(0),
where Rex=Ux/ν is the local Reynolds number.


## 3. Solution Approach

In this section we apply HAM to solve system (2.11)–(2.14). We assume that the solutions of f(η) and θ(η) can be expressed by a set of base functions {ηm∣m=0,1,2,…} with the following forms:(3.1)f(η)=∑m=0+∞amηm,θ(η)=∑m=0+∞cmηm,
where am and cm are constants. Under the rule of solution expression given by (2.11) and (2.12), subject to the boundary conditions (2.13) and (2.14), it is straightforward to choose(3.2)f0(η)=η+3S-6-Mθ0(η)4η2+2-S+Mθ0(η)4η3,θ0(η)=1,
as the initial guesses of f(η) and θ(η). The auxiliary linear operators Lf=∂3/∂η3 and Lθ=∂2/∂η2 are chosen with the following properties:(3.3)Lf[C1+C2η+C3η2]=0,Lθ[C4+C5η]=0,
where C1, C2, C3, C4, and C5 are constants of integration. From (2.11) and (2.12), the nonlinear operators are defined as(3.4)Nf[F(η,q),Γ(q)]=F′′′+Γ[FF′′-12SηF′′-(F′)2-SF′],Nθ[F(η,q),Θ(η,q),Γ(q)]=1Pr(1+NR)Θ″+Γ[FΘ′-2F′Θ-12SηΘ′-32SΘ+1PrB*Θ],
where F(η,q) and Θ(η,q) are both unknown functions of η and q, while Γ is a function dependent on q. The prime denotes differentiation with respect to η. Let ℏf and ℏθ denote the nonzero auxiliary parameters, whereas Hf and Hθ denote nonzero auxiliary functions. The zero-order deformation equation can be constructed as(3.5)(1-q)Lf[F(η,q)-f0(η)]=qℏfHfNf[F(η,q),Γ(q)],(1-q)Lθ[Θ(η,q)-θ0(η)]=qℏθHθNθ[F(η,q),Θ(η,q),Γ(q)],
subject to the boundary conditions(3.6)F(0,q)=0,F′(0,q)=1,Θ(0,q)=1,F(1,q)=12S,F′′(1,q)=Mθ(1),Θ′(1,q)=0,
where q is an embedding parameter. From (3.2), it is straightforward to show that when q=0, the solutions of (3.5) and (3.6) are(3.7)F(η,0)=f0(η),Θ(η,0)=θ0(η).
Since ℏf,ℏθ≠0 and Hf,Hθ≠0 when q=1, (3.5) and (3.6) are equivalent to (2.11)–(2.14), respectively, provided that(3.8)F(η,1)=f(η),Θ(η,1)=θ(η),Γ(1)=γ.
Thus, as q increases from 0 to 1, F(η,q) and Θ(η,q) vary from the initial guesses f0(η) and θ0(η) to the solutions f(η) and θ(η) in (2.11)–(2.14). Γ also varies from the initial guess(3.9)Γ(0)=γ0,
to the time-scale parameter γ. Applying the Taylor series to (3.7) and (3.9), F(η,q), Θ(η,q), and Γ(q) can be expanded as series of q,(3.10)F(η,q)=f0(η)+∑m=1+∞fm(η)qm,Θ(η,q)=Θ0(η)+∑m=1+∞θm(η)qm,Γ(q)=γ0+∑m=1+∞γmqm,
where(3.11)fm(η)=1m![∂mF(η,q)∂qm]q=0,θm(η)=1m![∂mΘ(η,q)∂qm]q=0,γm=1m![∂mΓ(q)∂qm]q=0.
Thus, using (3.8), we have(3.12)f(η)=f0(η)+∑m=1+∞fm(η),θ(η)=θ0(η)+∑m=1+∞θm(η),γ=γ0+∑m=1+∞γm.
By differentiating (3.5) m times with respect to q, then setting q=0, and finally dividing by m!, the mth-order deformation equations are obtained:(3.13)Lf[fm(η)-χmfm-1(η)]=ℏfHf(η)R1,m(η),Lθ[θm(η)-χmθm-1(η)]=ℏθHθ(η)R2,m(η),
subject to the boundary conditions(3.14)fm(0)=0,fm′(0)=0,θm(0)=0,fm(1)=0,fm′′(1)=Mθm(1),θm′(1)=0,
for m≥1, where(3.15)R1,m(η)=fm-1′′′+∑n=0m-1γm-1-n∑i=0n(fifn-i′′-fi′fn-i′)-12Sη∑n=0m-1γnfm-1-n′′-S∑n=0m-1γnfm-1-n′,R2,m(η)=1Prθm-1″+NRPrθm-1″+∑n=0m-1γm-1-n∑i=0n(fiθn-1′-2fn-i′θi)-12Sη∑n=0m-1γnθm-1-n′-32S∑n=0m-1γnθm-1-n+1PrB*∑n=0m-1γnθm-1-n,χm={1m>1,0m=1.
Let fm* and θm* denote the particular solutions of (3.13) as(3.16)fm*=∭0ηℏfHf(s)R1,m(s)dsdηdη+χmfm-1,θm*=∬0ηℏθHθ(s)R2,m(s)dsdη+χmθm-1.
The general solutions for (3.3) are(3.17)fm(η)=fm*+C1+C2η+C3η2,θm(η)=θm*+C4+C5η.
It should be pointed out that fm* and θm* contain the unknown parameter γm-1, which should be determined along with C1,C2, C3, C4, and C5 by the boundary conditions (3.14). We solve (3.13) for m=1,2,3,… successfully, and the mth-order approximations of f(η),θ(η), and γ are given, respectively, by(3.18)f(η)≈∑n=0mfn(η),θ(η)≈∑n=0mθn(η),γ≈∑n=0m-1γn.


## 4. Convergence of the HAM Solution

We note the presence of the auxiliary parametersℏf and ℏθ in (3.13). Liao [16, 25] indicates that the convergence and rate of approximation of such series depend on the values of ℏf and ℏθ. For the purpose of choosing the proper values for ℏf and ℏθ, we set S=1.4, Pr=1, M=1, NR=1, and B*=0.05. Figure 1(a) shows the variation of γ with ℏf using the tenth-order approximation. From the figure, it is seen that convergent result can be obtained when -0.8≤ℏf≤-0.35. Thus, we can choose an appropriate value for ℏf in this range to get convergent solution for γ. In Figure 1(b), we observe that convergent results of f″(0) can be obtained by choosing a value of ℏf in the range of -0.78≤ℏf≤-0.16 under the tenth-order approximation. From (2.13), we can see that ℏθ depends on ℏf. Figure 1(c) shows the appropriate value for ℏθ in the range of -0.75≤ℏθ≤-0.25 when ℏf=-0.6 under the tenth-order approximation. To see the accuracy of the solution, we define Ef(η) and Eθ(η) as the residual errors of the system. Figure 2 shows the residual errors of the HAM solution in the case ℏf=-0.6 and ℏθ=-0.4.(a)ℏf-curve of f″(0), (b) ℏf-curve of γ, and (c) ℏθ-curve of θ′(0) (ℏf=-0.6) using 10th-order HAM approximation for the case of S=1.4, Pr=1, M=1, NR=1, and B*=0.05.
(a)(b)(c)The residual errors for the solution using 15th-order HAM approximation for the case ofS=1.4, Pr=1, M=1, NR=1, and B*=0.05.
(a)(b)


## 5. Results and Discussion

We solved (3.13) and (3.14) with the aid of Maple, a symbolic computation software. The auxiliary functions Hf and Hθ in (3.13) were set to be equal to 1 in all calculations reported in this paper.In order to assess the accuracy of the numerical method, Table1 compares our present results with those of previous investigations. We can claim that our numerical results are in excellent agreement with those of [4, 14, 24] under some limiting cases.Table 1
Values of dimensionless film thicknessβ and skin friction f″(0) when Pr=1, M=0, NR=0, ℏf=-0.6, and ℏθ=-0.4 for different values of unsteadiness parameter S. The value of f″(0)/β in Wang [4] and the present work is the same as that of Abel et al. [14] despite different similarity transformations used.
Wang [4]Abel et al. [14]Mahmoud [24]Present workSβf′′(0)/ββf′′(0)βf′′(0)βf′′(0)/β0.45.122490−1.3077854.981455−1.134098——5.126821−1.0407650.63.131250−1.1951553.131710−1.195128——3.131668−1.1936540.82.151990−1.2457952.151990−1.2458052.1519950−1.2458102.151994−1.2457931.01.543617−1.2777621.543617−1.277769——1.543616−1.2777681.21.127780−1.2791771.127780−1.2791711.1277815−1.2791701.127780−1.2791721.40.821032−1.2335490.821033−1.233545——0.821032−1.2335491.60.567173−1.1149370.576176−1.114941——0.576173−1.1149371.80.356389−0.8674140.356390−0.867416——0.356389−0.867414In Table2, a good agreement was found in comparisons of values of free surface temperature θ(1) between [24] and present work.Table 2
Values of free surface temperatureθ(1) using 10th-order HAM approximation for several values of M, NR with Pr=0.1, B*=0, and S=1.2.
MNRθ(1)Mahmoud [24]Present work00.20.866810.8668440.10.20.851840.8542570.50.20.788150.8101030.100.826930.8296060.10.20.851840.8542570.11.00.905940.907657The effects of the radiation parameterNR on the film thickness β, surface shear stress f″(0), free surface temperature θ(1), and heat flux -θ′(0) at a specific thermocapillarity number M are shown in Table 3. It is concluded that the film thickness and the free surface temperature increase as the radiation parameter increases and consequently the surface shear stress and dimensionless heat flux decrease.Table 3
Variations ofγ(η), β, f′′(0), θ(1), and -θ′(0) using 10th-order HAM approximation when Pr=1, M=1, B*=0.05, S=1.4, and NR is varied.
NRγ(η)βf′′(0)θ(1)-θ′(0)01.0089411.004460−1.2011790.3077141.8418990.51.1036551.050550−1.2524140.3961791.5121001.01.1736161.083335−1.2898240.4617401.2998091.51.2272401.107809−1.3188630.5132591.1471572.01.2702511.127054−1.3415330.5542681.031198Table4 shows the effects of thermocapillarity number on the film thickness β, surface shear stress f″(0), free surface temperature θ(1), and heat flux -θ′(0). We agree with [24] that the thermocapillarity parameter has the effect of enhancing the velocity and the local Nusselt number while decreasing the temperature and the surface shear stress.Table 4
Variations ofγ(η), β, f′′(0), θ(1), and -θ′(0) using 10th-order HAM approximation when Pr=1, NR=1, B*=0.05, S=1.4, and M is varied.
Mγ(η)βf′′(0)θ(1)-θ′(0)00.6740940.821032−1.0127740.6187390.8734920.50.9690910.984424−1.1724200.5197001.1403711.01.1736161.083335−1.2998240.4617401.2998091.51.3500541.161918−1.3821590.4195341.4258352.01.4999331.224717−1.4597570.3880861.524753The heat absorption sink(B*<0) leads to a decrease in the thermal boundary layer  whereas the boundary layer thickness increases with increase in B*. These effects can be observed in Table 5. As the temperature-dependent parameter B* increases for both cases (heat source or heat sink), the free surface temperature increases while the surface shear stress and wall heat flux decrease.Table 5
Variations ofγ(η), β, f′′(0), θ(1), and -θ′(0) using 10th-order HAM approximation when Pr=1, NR=1, M=1, S=1.4, and B* is varied.
B*γ(η)βf′′(0)θ(1)-θ′(0)−1.01.1084641.052836−1.2541290.3993461.485058−0.51.1372351.066412−1.2700270.4269651.40120901.1700671.081696−1.2879060.4583541.3094860.51.2061361.099152−1.3083520.4946301.2078631.01.2530281.119387−1.3321920.5373951.093403The effects of the different values of radiation parameterNR on the temperature profiles are depicted in Figure 3. By fixing M=1, S=1.4, Pr=1, and B*=0.05, the thin film flow swings from slight deceleration to higher velocity while the temperature increases as NR increases.Effects of the radiation parameterNR on (a) f′(η) and (b) θ(η) using 10th-order HAM approximation for the case of ℏf=-0.6, ℏθ=-0.8, S=1.4, Pr=1, M=1, and B*=0.05.
(a)(b)Thermocapillarity produced an outward flow along the free surface. Figure4 shows that by increasing the value of the thermocapillarity number, the velocity decreases until it arrives at a minimum and then increases to its free surface value while the temperature consistently cools down.Effects ofM on (a) f′(η) and (b) θ(η) using 10th-order HAM approximation for the case of ℏf=-0.6, ℏθ=-0.8, S=1.4, NR=1, Pr=1, and B*=0.05.
(a)(b)Figure5 shows that when the temperature-dependent parameter B* increases, the velocity decreases while the temperature increases. This is due to the fact that the internal heat source enhances or damps heat transport.Effects ofB* on (a) f′(η) and (b) θ(η) using 10th-order HAM approximation for the case of ℏf=-0.6, ℏθ=-0.8, S=1.4, M=1.0, NR=1, and Pr=1.
(a)(b)The unsteadiness parameterS has the effect of increasing both velocity profiles and temperature distribution at a specific NR, B*, and M, as shown in Figure 6. The effects of Prandtl number Pr on temperature distribution is presented in Figure 7. As Pr increases, the velocity profile increases and the temperature distribution decreases. The higher the Prandtl number, the cooler the temperature of the fluid flow.Effects ofS on (a) f′(η) and (b) θ(η) using 10th-order HAM approximation for the case of ℏf=-0.6, ℏθ=-0.8, M=1.0, NR=1, Pr=1, and B*=0.05.
(a)(b)Effects ofPr on (a) f′(η) and (b) θ(η) using 10th-order HAM approximation for the case of ℏf=-0.6, ℏθ=-0.8, M=1.0, NR=1, S=1.4, and B*=0.05.
(a)(b)


## 6. Concluding Remarks

The effects of thermal radiation and thermocapillarity in a thin liquid film on an unsteady stretching sheet with nonuniform heat source/sink was analyzed successfully by means of the homotopy analysis method (HAM). With the presence of internal heat generation/absorption, the radiation parameter plays a significant role in controlling the temperature of the fluid flow by enhancing the temperature of fluid flow and decreasing the dimensionless heat flux. The thermocapillarity enhances the velocity and the dimensionless heat flux while reducing the temperature and the surface shear stress.



---

*Source: 127320-2012-02-08.xml*