Mixed Convection Boundary Layer Stagnation-Point Flow of a Jeffery Fluid Past a Permeable Vertical Flat Plate Mohammad M. Rahman

Mixed Convection Boundary Layer Stagnation-Point Flow of a Jeffery Fluid Past a Permeable Vertical Flat Plate

Mohammad M. Rahman^aand Ioan Pop^b

aDepartment of Mathematics and Statistics, College of Science, Sultan Qaboos University, P.O.

Box 36, P.C. 123 Al-Khod, Muscat, Sultanate of Oman

bDepartment of Mathematics, Faculty of Mathematics and Computer Science, Babe¸s-Bolyai University, 400084 Cluj-Napoca, Romania

Reprint requests to I. P.; E-mail:popm.ioan@yahoo.co.uk

Z. Naturforsch.69a, 687 – 696 (2014) / DOI: 10.5560/ZNA.2014-0065

Received June 4, 2014 / revised August 18, 2014 / published online November 5, 2014

This paper analyzes the combined effects of buoyancy force, mass flux, and variable surface tem- perature on the stagnation point flow and heat transfer due to a Jeffery fluid over a vertical surface.

The governing nonlinear partial differential equations are transformed into a system of coupled non- linear ordinary differential equations using similarity transformations and then solved numerically using the function bvp4c from computer algebra software Matlab. Numerical results are obtained for skin friction coefficient, Nusselt number as well as dimensionless velocity and temperature profiles for various values of the controlling parameters namely mixed convection parameterλ, mass flux parameters, elastic parameter (Deborah number)γ, and the ratio of relaxation and retardation time parameterλ1. The results indicate that dual solutions exist in a certain range of the mixed convection and mass flux parameters. In order to establish the physically realizable of these solutions, a sta- bility analysis has also been performed. The results indicate that mixed convection and mass flux significantly affects the nature of the solutions, skin friction, and Nusselt number of a Jeffery fluid.

Key words:Mixed Convection; Boundary Layer; Jeffery Fluid; Stagnation-Point Flow; Dual Solutions.

1. Introduction

In non-Newtonian fluids the flow properties espe- cially the fluid viscosity differs from those of New- tonian fluids. In these fluids the viscosity depends on shear rate or shear rate history (Bird et al. [1]; Cole- man and Noll [2] and Harris [3]). Although there could be some non-Newtonian fluids with shear-independent viscosity, they still exhibit normal stress differences or other non-Newtonian behaviour. In non-Newtonian fluids, the constitutive relationships between the stress tensor and the rate of strain tensor are much com- plicated in comparison to the classical Navier–Stokes equations. Due to the extensive applications in science, industry, and technology (for example, molten poly- mer, starch suspensions, productions of toiletries and paints, thermal oil recovery, slurry transportation, food processing, and characterizing animal blood, etc.) a va- riety of non-Newtonian fluid models can be found in the open literature. A Jeffery fluid [4] is a kind of non- Newtonian fluid that is capable of describing the char-

acteristics of relaxation and retardation times and has attracted attentions by many investigators ([5–11], and the references therein) because of its simplicity. This fluid model is able to describe the characteristics of re- laxation and retardation times which arise in complex polymeric flows. Furthermore the Jeffrey type model utilizes time derivatives rather than convected deriva- tives, which greatly facilitates numerical simulations.

This flow model provides an elegant formulation for simulating retardation and relaxation effects arising in non-Newtonian flows. There were recently published several paper of Jeffery fluid, such as, for example, Hayat et al. [12,13], Shehzad et al. [14–16], etc.

The steady viscous flow in the neighbourhood of a stagnation point on a fixed wall, first studied by Hiemenz [17], is very important for many practical ap- plications in the modern industry. The stagnation flow occurs whenever the flow impinges on any solid ob- ject and the local velocity of the fluid at the stagna- tion point is zero. A large number of calculations on this problem including different effect into the govern-

© 2014 Verlag der Zeitschrift für Naturforschung, Tübingen·http://znaturforsch.com

ing equations have been done. Excellent review of ex- isting theoretical and experimental work on this sub- ject can be found in the books by Gebhart et al. [18], White [19], Shlichting and Gersten [20], Pop and In- gham [21], Bejan [22], etc. There are also several in- teresting published papers on the stagnation point flow as can be seen in Merkin and Mahmood [23], Ishak et al. [24], etc. On the other hand, it should be men- tioned that mixed convection flows, or the combination of both forced and free convection, arise in many trans- port processes in engineering devices and in nature, in- cluding solar receivers exposed to wind currents, elec- tronic devices cooled by fans, nuclear reactors cooled during emergency shutdown, heat exchanges placed in a low-velocity environment, etc. Such processes occur when the effect of the buoyancy force in forced con- vection or the effect of a forced flow in free convection becomes a significant feature. It is already well estab- lished that thermal buoyancy forces play a significant role in forced convection heat transfer when the flow velocity is relatively low and the temperature differ- ence between the surface and the free stream is rel- atively large. Buoyancy forces also play a significant role in the onset of flow instabilities, as they can be responsible for either delaying or speeding up the tran- sition from laminar to turbulent flow (see Chen and Ar- maly [25]). These authors have shown that the domain of the mixed convection regime is generally defined as the region a≤Gr/Reⁿ≤b, where Gr is the Grashof number, Re is the Reynolds number, a andb are, re- spectively, the lower and the upper bounds of the do- main, and nis a constant which depends on the flow configuration and the surface heating conditions. The buoyancy parameter Gr/Reⁿprovides a measure of the influence of free convection in comparison to that of forced convection on the flow. Outside the mixed con- vection region,a≤Gr/Reⁿ≤b, either a pure forced convection or a pure free convection analysis can be used to describe the flow and the temperature field ac- curately.

The objective of this paper is to study the steady mixed convection boundary layer flow and heat trans- fer of a Jeffery fluid near the stagnation point on a per- meable vertical surface. Particular attention is given to deriving numerical results for the critical/turning points which determine the range of existence of mul- tiple (dual) upper and lower branch solutions. In this respect, a stability analysis of the dual solutions has been performed. This new type of flow is essentially

a backward flow as discussed by Goldstein [26]. Fur- ther, it is worth mentioning that suction or injection of a fluid through the bounding surface, as, for example, in mass transfer cooling, can significantly change the flow field and, as a consequence, affect the heat trans- fer rate at the plate. In general, suction tends to increase the skin friction and heat transfer coefficients, whereas injection acts in the opposite manner. Injection or with- drawal of a fluid through a porous bounding heated or cooled wall is of general interest in practical problems involving boundary layer control applications such as film cooling, polymer fiber coating, coating of wires, etc. We mention to this end that the following relevant very recently studies for the mixed convection flow may be also mentioned: Hayat et al. [27–29], Awais et al. [30], Shehzad et al. [31], etc.

2. Basic Equations

Following Bird et al. [1], the Cauchy stress tensorΞΞΞ of a Jeffery non-Newtonian fluid takes the form

ΞΞΞ=−pI+S, (1)

where

S= µ

1+λ1 (γ˙+λ2γ)¨ (2) and

γ˙=∇V+ (∇V)^T,

¨ γ= d

dt(γ) =˙ ∂

∂t(γ˙) + (V·∇)γ˙. (3) HereSis the extra stress tensor,pis the pressure,Iis the density tensor,Vis the velocity vector,µis the dy- namic viscosity,λ1is the ratio of relaxation to retarda- tion time,λ2is the retardation time, ˙γis the shear rate, and a dot above a quantity denotes the material time derivative. Following these relations, we consider the two-dimensional boundary layer flow of a Jeffery fluid near the stagnation point over a permeable vertical sur- face coinciding with the plane y=0, the flow being confined toy>0, where theycoordinate is measured in the normal direction to the vertical surface, as shown in Figure1. It is assumed that the velocity distribu- tion of the external flow (ambient fluid) isue(x) =ax, wherexis the coordinate measured along surface, and a is a positive constant (a>0). It is also assumed that the constant mass flux velocity isv0 withv0<0

x

y 0

T∞

u_e

T_w (x) > T x

∞

v₀ g u

u v₀

(x) < T∞

T_w T_w(x) < T∞

T∞

x

y

v₀ 0 u u

u_e g

v₀ u

T_w(x) < T∞

(a) (b)

Fig. 1. Flow configuration and coordinate system (a) assisting flow, (b) opposing flow.

for suction andv0>0 for injection or withdrawal of the fluid, respectively. Also, the temperature of plate is T_w(x), while the uniform temperature of the ambi- ent fluid is T_∞. Under these conditions, the boundary layer equations, which govern this problem, are (see Qasim [11])

∂u

∂x+∂v

∂y=0, (4)

u∂u

∂x+v∂u

∂y=uedue

dx + υ

1+λ₁

"

∂²u

∂y²+λ2

u ^∂3u

∂x∂y²+v^∂3u

∂y³

−^∂_∂x^u∂_∂y^2u₂ +^∂u

∂y ∂²u

∂x∂y

!#

+gβ(T−T_∞), (5)

u∂T

∂x+v∂T

∂y =α∂²T

∂y², (6)

subject to the bounday conditions (see Garg and Ra- jagopal [32])

v=v0, u=0, T=Tw(x) =T_∞+T0x at y=0, u=ue(x)→ax, ∂u

∂y→0, T→T_∞ as y→∞, (7) whereuandvare the velocity componets alongxand yaxes,T is the fluid temperature,T0is a temperature characteristic withT₀>0 for a heated plate (assisting flow) andT0<0 for a cooled plate (opposing flow) sur- face, respectively,α is the thermal diffusivity,gis the acceleration due to gravity, and β is the thermal ex- pansion coefficient. Hence the initial studies for mixed convection flows should be highlighted. Some relevant studies may be mentioned as follows.

In order to solve (4) – (6) along with the boundary conditions (7), we introduce the following similarity variables:

ψ=√

aνx f⁰(η),

θ(η) = (T−T_∞)/(T_w−T_∞), η=p

a/νy, (8) whereψ is the stream function which is defined in the usual way asu=∂ ψ/∂yandv=−∂ ψ/∂x. Thus, we have

u=ax f⁰(η), v=−√

aνf(η), (9)

where prime denotes differentiation with respect toη.

Substituting (8) into (5) and (6), the following set of ordinary differential equations results in

f⁰⁰⁰+ (1+λ1)(f f⁰⁰+1−f⁰²)

+γ(f⁰⁰²−f f⁰⁰⁰⁰) + (1+λ1)λ θ=0, (10) 1

Prθ⁰⁰+fθ⁰−f⁰θ=0, (11)

and the boundary conditions (4) become f(0) =s, f⁰(0) =0, θ(0) =1,

f⁰(η)→1, f⁰⁰(η)→0, θ(η)→0 as η→∞. (12) Here Pr=ν/αis the Prandtl number,λis the constant mixed convection parameter,γis the elastic parameter (or Deborah number), andsis the constant suction(s>

0)or injection(s<0)parameter, which are defined as λ=gβ(Tw−T_∞)x³/ν²

(uex)²/ν² =Grx

Re²_x, γ=aλ2, s=− v0

(aν)^1/2,

(13)

with Grx=gβ(T_w−T_∞)x³/ν²being the local Grashof number, and Rex=u_e(x)x/ν is the local Reynolds number. We notice thatλ =0 corresponds to forced convection flow, λ >0 corresponds to assisting flow (heated plate), andλ<0 corresponds to opposing flow (cooled plate), respectively.

Physical quantities of interest are the skin friction coefficientCfand the local Nusselt number Nux, which are defined as

Cf= τw

ρu²_e, Nux= xq_w

k(T_w−T_∞), (14) whereτwandqware the wall skin friction and the heat transfer from the plate, which are given by

τw= µ 1+λ₁

∂u

∂y+λ₂

u ∂²u

∂y∂x+v∂²u

∂y²

y=0

, q_w=−k

∂T

∂y

y=0

.

(15)

Using variables (8), we get Re^1/2x C_f= 1

1+λ1

f⁰⁰(0)−γs f⁰⁰⁰(0) , Re^−1/2x Nux=−θ⁰(0).

(16)

3. Flow Stability

To investigate the stability of the problem let us con- sider the corresponding unsteady problems of (5) and (6) as follows:

∂u

∂t +u∂u

∂x+v∂u

∂y=uedue

dx + υ

1+λ1

"

∂²u

∂y²+λ2

∂³u

∂y²∂t+u_∂x∂^∂3u_y2+v^∂_∂y^3u₃

−^∂u

∂x∂²u

∂y²+^∂u

∂y ∂²u

∂x∂y

!#

+gβ(T−T_∞),

(17)

∂T

∂t +u∂T

∂x +v∂T

∂y =α∂²T

∂y² , (18)

subject to the boundary conditions

t<0 : u=0, v=v₀, T=T_∞ for anyx,y, t≥0 : u=0, v=v₀, T=T_w(x) =T_∞+T₀x

at y=0, (19)

u=u_e(x)→ax, ∂u

∂y→0, T→T_∞ as y→∞.

Let us consider the following transformations:

u=ax f⁰(η,τ), v=−√

aυf(η,τ), θ(η,τ) = T−T_∞

T_w−T_∞, τ=at, η= ra

υy, (20) whereτ=atis a dimensionless time variable. Substi- tuting (20) into (17) and (18), we obtain

∂³f

∂ η³+ (1+λ1)

"

f∂²f

∂ η²+1− ∂f

∂ η 2#

−(1+λ1) ∂²f

∂ η ∂ τ +γ

"

∂²f

∂ η^{2 2}

−f∂⁴f

∂ η⁴+ ∂⁴f

∂ η³∂ τ

#

+ (1+λ1)λ θ=0, (21) 1

Pr

∂²θ

∂ η²+f∂ θ

∂ η −∂f

∂ ηθ−∂ θ

∂ τ =0. (22)

The boundary conditions (19) become f(0,τ) =s, ∂f

∂ η =0, θ(0,τ) =1,

∂f

∂ η(η,τ)→1, ∂²f

∂ η²(η,τ)→0, θ(η,τ)→0 as η→∞.

(23)

To test stability of the steady flow solution f(η) = f0(η)andθ(η) =θ0(η)satisfying the boundary-value problem (10) – (12), we write (see Weidman et al. [33], Postelnicu and Pop [34], and Ro¸sca and Pop [35,36])

f(η,τ) =f₀(η) +e^{−σ τ}F(η,τ),

θ(η,τ) =θ0(η) +e^{−σ τ}G(η,τ), (24) whereσis an unknown eigenvalue parameter,F(η,τ) andG(η,τ)are small relative tof0(η)andθ0(η). Sub- stituting (24) into (21) and (22), we obtain the follow- ing linearized problem:

∂³F

∂ η³+ (1+λ1)

f₀∂²F

∂ η²+F∂²f₀

∂ η² −2∂f₀

∂ η

∂F

∂ η

−(1+λ1) ∂²F

∂ η ∂ τ −σ∂F

∂ η

+γ

2∂²f0

∂ η²

∂²F

∂ η²−f0∂⁴F

∂ η⁴−F∂⁴f0

∂ η⁴

+γ ∂⁴F

∂ η³∂ τ−σ∂³F

∂ η³

+ (1+λ₁)λG=0,

(25)

1 Pr

∂²G

∂ η²+f0∂G

∂ η +F∂ θ0

∂ η −∂F

∂ ηθ0

− ∂f₀

∂ η −σ

G−∂G

∂ τ =0,

(26)

along with the boundary conditions F(0,τ) =0, ∂F

∂ η(0,τ) =0, G(0,τ) =0,

∂F

∂ η(η,τ)→0, ∂²F

∂ η²(η,τ)→0, G(η,τ)→0 as η→∞.

(27)

As suggested by Weidman et al. [33], we investigate the stability of the steady flow and heat transfer so- lution f0(η) andθ0(η)by setting τ=0, and hence F =F₀(η)andG=G₀(η)in (10) and (11) to iden- tify initial growth or decay of the solution (24). To test our numerical procedure, we have to solve the linear eigenvalue problem

(1−γ σ)F₀⁰⁰⁰+ (1+λ1)

f0F₀⁰⁰+F0f₀⁰⁰− 2f₀⁰−σ F₀⁰ +γ2f₀⁰⁰F₀⁰⁰−f₀F₀⁰⁰⁰⁰−F₀f₀⁰⁰⁰⁰

+ (1+λ₁)λG₀=0, (28) 1

PrG⁰⁰₀+f₀G⁰₀+F₀θ₀⁰−F₀⁰θ0− f₀⁰−σG₀=0, (29) along with the boundary conditions

F₀(0) =0, F₀⁰(0) =0, G₀(0) =0, F₀⁰(η)→0, F₀⁰⁰(η)→0,

G(η)→0 as η→∞.

(30)

For particular values of Pr,λ,s,γ, andλ1, the sta- bility of the corresponding steady flow solutions f0(η) andθ0(η)is determined by the smallest eigenvalueσ. As it has been done by Harris et al. [37], the range of possible eigenvalues can be determined by relax- ing a boundary condition onF0(η)orG0(η). For the present problem, we relax the condition thatF₀⁰⁰(η)→ 0 asη→∞. Thus, for a fixed value ofσwe solve the problem (28) – (30) along with the new boundary con- ditionF₀⁰⁰(0) =1.

4. Numerical Technique

Following Ro¸sca and Pop [35,36] and Rahman et al. [38,39] the system of ordinary differential equa- tions (10) and (11) subject to the boundary conditions

(12) is solved numerically using the function bvp4c from Matlab for different values of the controlling pa- rameters. The numerical simulations are carried out for various values of the physical parameters mixed con- vection parameterλ, mass flux parameters, elastic pa- rameter (or Deborah number)γ, and the ratio of re- laxation and retardation time parameterλ1for a fixed value of the Prandtl number Pr. Because of the lack of experimental data, the choice of the values of the pa- rameters was dictated by the values chosen by previous investigators. The value of the Prandtl number is set equal to 1 throughout the paper. The values of the other parameters are mentioned in the description of the re- spective figures. The code bvp4c is developed using fi- nite difference method that implements the three-stage Lobatto IIIa formula, which is a collocation method with forth-order accuracy. In this approach, the ordi- nary differential (10) and (11) are reduced to a first- order system by introducing new variables. The mesh selection and error control are based on the residual of the continuous solution. The relative error tolerance has been set to 10⁻⁷. Because the present problem may have more than one solution (dual, upper and lower branch solutions), a ‘good’ initial guess is necessary.

The ‘infinity’η→∞in the boundary condition (9) is replaced by a finite valueη=η∞. We started the com- putation at small value, for example,η=5, then sub- sequently increased the value ofη until the boundary conditions are verified. In this method, we have chosen a suitable finite value ofη→∞, namelyη=η_∞=20 for the upper branch (first) solution and η =η∞ in the range 40 – 50 for the lower branch (second) solu- tion. Examples of solving boundary value problems by bvp4c can be found in the book by Shampine et al. [40]

or through online tutorial by Shampine et al. [41].

It is good to mention that forλ=λ1=γ=s=0, (7) f⁰⁰⁰+f f⁰⁰+1−f⁰²=0 along with the boundary condi- tions, (12) f(0) =0, f⁰(0) =0, f⁰(∞)→1 coincides with the classical Falkner–Skan wedge flows (4 – 71) with boundary conditions (4 – 72) of White [19]. For m=1, White [19] reportedf⁰⁰(0) =1.23259 which ex- actly matches with our calculated value.

5. Results and Discussion

The numerical simulations of (10) and (11) subject to the boundary conditions (12) are carried out for var- ious values of the physical parametersλ,s,γ, andλ1

when Pr=1 for obtaining the condition under which the

dual (upper and lower branch) solutions for the steady mixed convection flow of a Jeffery fluid over a vertical flat plate with variable surface temperature may exist.

In Figures 2 to 9, we have investigated the vari- ation of f⁰⁰(0) and−θ⁰(0) with λ for different val- ues of s, γ, and λ1 keeping Pr=1 fixed. These fig- ures clearly demonstrate the parameter space for the existance of the dual solutions. Figures 2and3 show that for an opposing flow (λ <0)dual solutions (up- per and lower branch solutions) can be found for the studied parameters when the cirical mass flux param- eter sc≥1.0477. Fors<sc there exits only one so- lution that is the upper branch (solid line) solution.

From Figures4to5we found that for a fixed value of the mass flux parametersdual solutions can be found when λc<λ <0, one solution whenλ =λc, and no solutions whenλ>λc, whereλcis the critical value of

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

−12

−10

−8

−6

−4

−2 0 2 4

s f ′′(0) s_c = 1.0477

upper branch solution lower branch solution

Fig. 2 (colour online). Values of f⁰⁰(0)versussfor Pr=1, λ=−1,γ=0.2,λ₁=0.2.

−40 −35 −30 −25 −20 −15 −10 −5 0

−6

−5

−4

−3

−2

−1 0 1 2 3

λ f ′′(0)

s = 1.5, 2, 2.5

λ_c = −12.5906 λ_c = −22.2357

λ_c = −38.4160 lower branch solution upper branch solution

Fig. 4 (colour online). Values of f⁰⁰(0)versusλfor variouss when Pr=1,γ=0.2,λ₁=0.2.

λ. Thus, forλ <λcthe full Navier–Stokes equations and energy equation need to be solved. The critical val- ues ofλ areλc=−12.5906,−22.2357,−38.4160 for s=1.5, 2, 2.5, respectively. The|λ_c|increases with the increase ofs. Therefore, suction at the surface broad- ens the parameter space inλfor the existence of the so- lution. Figure4reveals that at anyλ station within the domainλc<λ<0 the values of f⁰⁰(0)for the upper branch solution increase with the increase ofswhereas as opposite result is found for the lower branch (dotted line) solution. For a fixed value ofsthe upper branch solution increases with the increase of the mixed con- vection parameter λ for an opposing flow. The cor- responding lower branch solution first decreases with the increase ofλ whenλc≤λ <λ1(say) then it be- gins to increase with the further increase ofλ. On the other hand, for the upper branch solution the values

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

−60

−50

−40

−30

−20

−10 0 10

s

−θ′(0)

s_c = 1.0477

upper branch solution

lower branch solution

Fig. 3 (colour online). Values of−θ⁰(0)versussfor Pr=1, λ=−1,γ=0.2,λ1=0.2.

−40 −35 −30 −25 −20 −15 −10 −5 0

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2 2.5 3 3.5

λ

−θ′(0)

λ_c = −38.4160

upper branch solution

lower branch solution

λ_c = −12.5906

λ_c = −22.2357 s = 1.5, 2, 2.5

Fig. 5 (colour online). Values of−θ⁰(0)versusλfor various swhen Pr=1,γ=0.2,λ₁=0.2.

of−θ⁰(0)(Fig.5) show a similar behaviour to that of the corresponding f⁰⁰(0). But for the lower branch so- lution these values decrease very rapidly with the in- crease of λ. It is good to mention that the values of

−θ⁰(0)for the lower branch solution are found very large compared to the corresponding values of the up- per branch solution in the vicinity of λ =0. To vi- sulalize both the solution branches we have truncated the lower branch solution as a consiquence it shows asymptotic behaviour as λ →0. Otherwise with the same scale the upper branch solutions will not be dis- tinguishable. From the stability analysis it is found that the upper solution branch is stable and physically re- alizable whereas the lower solution branch is unstable and physically not realizable. In Table1we have calcu- lated the smallest eigenvalues (σ)for the upper branch solutions for an opposing flow for various values of the

−45 −40 −35 −30 −25 −20 −15 −10 −5 0

−6

−5

−4

−3

−2

−1 0 1 2 3

λ f ′′(0)

upper banch solution

λ_c = −40.712

λ_c = −22.2357 λ_c = −14.294

γ = 0.01, 0.2, 1 lower banch solution γ = 0.01, 0.2, 1

Fig. 6 (colour online). Values of f⁰⁰(0)versusλfor variousγ when Pr=1,λ₁=0.2,s=2.

−25 −20 −15 −10 −5 0

−10

−8

−6

−4

−2 0 2 4

λ f ′′(0)

upper banch solution lower banch solution

λ₁ = 0.01, 1, 3 λ_c = −23.2491

λ_c = −19.7666

λ₁ = 0.01, 1, 3 λ_c = −17.381

Fig. 8 (colour online). Values of f⁰⁰(0)versusλ for various λ₁when Pr=1,γ=0.2,s=2.

parametersγandλ1. Due to the smallness of the thick- ness of the hydrodynamic and thermal boundary layers thicknessess it is extreamly difficult to find the small- est eigenvalues for the lower branch solution, hence not calculated here.

Figures 6 to 7 illustrate the effects of the elastic parameter, i.e. Deborah numberγ and λ on the val- ues of f⁰⁰(0)and−θ⁰(0), respectively. Small values of Table 1. Smallest eigenvalues (σ)for the upper branch solu- tion for different values ofγandλ₁when Pr=1,λ =−1, s=2.

γ λ1 σ

0.01 0.2 0.2

0.2 0.2 0.01

0.54464 0.20240 0.09120

−45 −40 −35 −30 −25 −20 −15 −10 −5 0

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2 2.5 3

λ

−θ′(0)

γ = 0.01, 0.2, 1

upper banch solution lower banch solution λ_c = − 40.712

− 22.2357 = λ_c

λ_c = −14.294

Fig. 7 (colour online). Values of−θ⁰(0)versusλfor various γwhen Pr=1,λ1=0.2,s=2.

−25 −20 −15 −10 −5 0

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2 2.5 3

λ

−θ′(0)

upper banch solution λ_c = −23.2491

λ₁ = 3, 1, 0.01

lower banch solution λ_c = −19.7666

λ_c = −17.381

Fig. 9 (colour online). Values of−θ⁰(0)versusλfor various λ₁when Pr=1,γ=0.2,s=2.

γ(≤1)characterize the liquid-like behaviour whereas the solid-like behaviour is associated with the large values of γ(>1). Here we have considered only the small values of γ. These figures reveal that values of

|λc|increase with the increase of the elastic parameter.

Figure6further demonstrates that values of f⁰⁰(0)for the upper branch solution increase with the increase of γwhenλc≤λ<λc₁(not precisely determined) and de- crease whenλc₁≤λ<0 with the increase of the elas- tic parameter. For the lower branch solution these val- ues increase with the increase ofγ. From Figure7we found that the effect ofγ on the values of −θ⁰(0)for the upper branch solution is similar to its effect on the values of f⁰⁰(0). The lower branch solution decreases very rapidly with the increase of the elastic parameter (Deborah number) γ, and the ratio of relaxation and retardation time parameterλ1.

0 1 2 3 4 5 6

−1.5

−1

−0.5 0 0.5 1

η f ′(η)

s = 1.5, 2, 2.5 s = 1.5, 2, 2.5

upper branch solution

lower branch solution

Fig. 10 (colour online). Velocity profiles for variousswhen Pr=1,λ=−1,γ=0.2,λ₁=0.2.

0 1 2 3 4 5 6 7

−1.5

−1

−0.5 0 0.5 1

η f ′(η)

γ = 1, 0.2, 0.01

upper banch solution

lower banch solution

Fig. 12 (colour online). Velocity profiles for variousγwhen Pr=1,λ=−1,λ₁=0.2,s=2.

The variation off⁰⁰(0)and−θ⁰(0)againstλ for dif- ferent values of the ratio of relaxation to retardation time parameterλ1are illustrated in Figures8to9, re- spectively. These figures show that values of criticalλ are−23.2491,−19.7666,−17.381 forλ1=0.01, 1, 3, respectively. Thus, |λ_c|decrease with the incraese of λ1. The values of f⁰⁰(0)for the upper branch solution decrease with the increase ofλ1whenλc≤λ≤ −8.63.

Outside this domain a reverse trend is observed. On the other hand the lower branch solution decreases with the increase ofλ1 for all values of λc≤λ ≤0. The values of the reduced Nusselt number−θ⁰(0)for the upper branch solution (Fig.9) increase with the de- crease ofλ1 whenλc≤λ ≤ −12.31 whereas an op- posite trend is observed when−12.31<λ ≤0. This figure further shows that the lower branch solution de- creases with the decrease ofλ1.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

0 0.5 1 1.5 2 2.5 3 3.5

η

θ(η)

upper branch solution lower branch solution

s = 1.5, 2, 2.5

s = 1.5, 2, 2.5

Fig. 11 (colour online). Temperature profiles for various s when Pr=1,λ=−1,γ=0.2,λ₁=0.2.

0 1 2 3 4 5 6

0 0.5 1 1.5 2 2.5 3 3.5 4

η

θ(η)

upper banch solution

γ = 0.01, 0.2, 1

γ = 0.01, 0.2, 1

lower banch solution

Fig. 13 (colour online). Temperature profiles for variousγ when Pr=1,λ=−1,s=2,λ₁=0.2.

0 1 2 3 4 5 6

−2.5

−2

−1.5

−1

−0.5 0 0.5 1

η f ′(η)

upper banch solution λ₁ = 0.01, 1, 4

lower banch solution

Fig. 14 (colour online). Velocity profiles for variousλ₁when Pr=1,λ=−1,γ=0.2,s=2.

The variations of the different controlling parame- ters on the fluid velocity and temperature distributions within the boundary layer are shown graphically in Figures10to15. The effects of the mass flux parameter on these functions are illustrated in Figures 10to11.

We found that the fluid velocity for the upper branch solution increases with the increase of the mass flux parameter while the lower branch solution decreases in the vicinity of the surface with the increase of s.

The thickness of the hydodynamic boundary layer de- creases with the increase ofsas expected for the up- per branch solution. The thickness of the upper branch solution is lower than the corresponding thickness of the lower branch solution. This gives the physical ac- ceptance of the upper branch solution. The tempera- ture profile within the boundary layer decreases, as a consequence the thickness of the thermal boundary layer also decreases with the increase ofsfor the up- per branch solution. Although the thickness of the ther- mal boundary layer for the lower branch solution de- creases with the increase ofsit can not be accepted as it is higher than the corresponding value of the upper branch solution.

The variation of the elastic parameter (Deborah number) on the Jeffery fluid velocity and temperature distributions are investiagted in Figures 12to13, re- spectively. Figure 12shows that the fluid velocity as well as the hydrodynamic boundary layer thickness for the upper branch solution increases with the decrease of γ. This is due to the fact that the greater the Deb- orah number, the more solid the material; the smaller the Deborah number, the more fluid it is. The lower branch solution although it decreases with the decrease

0 1 2 3 4 5 6

0 0.5 1 1.5 2 2.5 3

η

θ(η)

upper banch solution λ₁ = 0.01, 1, 4

λ₁ = 0.01, 1, 4

lower banch solution

Fig. 15 (colour online). Temperature profiles for variousλ1

when Pr=1,λ=−1,s=2,γ=0.2.

ofγwithin 0≤η≤ηcrit, it then further increases with the decrease ofγ. The temperature of the Jeffery fluid as well as the thickness of the thermal boundary layer increases with the increase ofγ as expected.

Figures 14 to 15 demonstrate the variation of fluid velocity and temperature distributions within the boundary layer for various values ofλ1and fixed val- ues of the other parameters. These figures clearly show that momentum as well as thermal boundary layer thicknessess decrease with the increase ofλ1.

6. Conclusions

In this paper we investigated the steady mixed con- vective boundary layer flow and heat transfer charac- teristics of a Jeffery fluid considering variable surface temperature over a vertical flat plate in the presence of mass flux at the surface. The numerical simulation was carried out to investigate the existence of the dual so- lutions. The critical suction and mixed convection pa- rameters for an opposing flow have been identified for the existence of the dual solutions. From this thorough investigation the major points can be summarized as follows:

(i) Dual solutions are found when the cirical mass flux parametersc≥1.0477. Fors<scthere exits only one solution that is the upper branch solu- tion.

(ii) Forλ >λcthe solutions have both the upper and lower branches, one solution whenλ=λcand no solution whenλ <λc.

(iii) The upper branch solution is found stable hence physically acceptable.

(iv) The elastic parameter (or Deborah number) β significantly influences the flow and heat trans- fer characteristics within the boundary layer flow of a Jeffery fluid.

(v) Intensification of the ratio of the relaxation to retardation times parameter λ1 decreases the mometum and thermal boundary layer thick- nessess.

[1] R. B. Bird, R. C. Armstrong, and O. Hassager, Dynam- ics of Polymeric Liquids, Vol 1, Fluid Mechanics, Wi- ley, New York 1977.

[2] B. D. Coleman and W. Noll, Arch. Ration. Mech. An.

56, 355 (1974).

[3] J. Harris, Rheology and Non-Newtonian Flow, Long- man, London 1977.

[4] G. B. Jeffery, Proc. R. Soc. Lon. Ser.-A 102, 161 (1922).

[5] T. Hayat, N. Ahmad, and N. Ali, Commun. Nonlinear Sci.13, 1581 (2008).

[6] M. Kothandapani and S. Srinivas, Int. J. Nonlinear Mech.43, 915 (2008).

[7] S. Nadeem and N. S. Akbar, Z. Naturforsch.64a, 713 (2009).

[8] T. Hayat, S. A. Shehzad, M. Qasim, and S. Obaidat, Nucl. Eng. Des.243, 15 (2012).

[9] S. Nadeem, S. Zaheer, and T. Fang, Numer. Algorithms 57, 187 (2011).

[10] T. Hayat and N. Ali, Commun. Nonlinear Sci.13, 1343 (2008).

[11] M. Qasim, Alexandria Eng. J.52, 571 (2013).

[12] T. Hayat, M. Nawaz, and A. A. Hendi, Z. Naturforsch.

66a, 53 (2011).

[13] T. Hayat, S. A. Shehzad, M. Quasim, and A. Alsaedi, Braz. J. Chem. Eng.31, 109 (2014).

[14] S. A. Shehzad, A. Alsaedi, and T. Hayat, Int. J. Heat Mass Trans.55, 3971 (2012).

[15] S. A. Shehzad, M. Qasim, A. Alsaedi, T. Hayat, and M. S. Alhouthuali, Eur. Phys. J. Plus128, 56 (2013).

[16] S. A. Shehzad, A. Alsaedi, T. Hayat, and M. S. Al- houthuali, J. Taiwan Inst. Chem. E.45, 787 (2014).

[17] K. Hiemenz, Dinglers Polytech. J.326, 321 (1911).

[18] B. Gebhart, Y. Jaluria, R. L. Mahajan, and B. Sam- makia, Buoyancy-Induced Flows and Transport, Hemi- sphere, New York 1988.

[19] F. M. White, Viscous Fluid Flow, 2nd edn., McGraw- Hill, New York 1991.

[20] H. Schlichting and K. Gersten, Boundary Layer The- ory, 8th edn., Springer, New York 2000.

[21] I. Pop and D. B. Ingham, Convective Heat Transfer:

Mathematical and Computational Modelling of Vis- cous Fluids and Porous Media, Pergamon, Oxford 2001.

[22] A. Bejan, Convection Heat Transfer, 4th edn., Wiley, New York 2014.

[23] J. H. Merkin and T. Mahmood, J. Appl. Math. Phys.

(ZAMP)40, 51 (1989).

[24] A. Ishak, J. H. Merkin, R. Nazar, and I. Pop, J. Appl.

Math. Phys. (ZAMP)59, 100 (2008).

[25] T. S. Chen and B. F. Armaly, Mixed convection in ex- ternal flow. In: Handbook of Ph Single-Phase Con- vective Heat Transfer, S. Kakaç, R. K. Shah, W. Aung, (eds.), pp. 14.1 – 14.35, Wiley, NewYork 1987.

[26] S. Goldstein, J. Fluid Mech.21, 33 (1965).

[27] T. Hayat, S. A. Shehzad, Ansa Rafique, and M. Y. Ma- lik, Int. J. Numer. Meth. Fl.68, 483 (2012).

[28] T. Hayat, N. Saleem, and A. A. Hendi, Chinese Phys.

Lett.29, 114704 (2012).

[29] T. Hayat, M. Waqas, S. A. Shehzad, and A. Alsaedi, J. Mech.29, 403 (2013).

[30] M. Awais, T. Hayat, and A. Qayyum, Eur. Phys. J. Plus 128, 114 (2013).

[31] S. A. Shehzad, A. Alsaedi, and T. Hayat, Braz. J. Chem.

Eng.30, 897 (2013).

[32] V. K. Garg and K. R. Rajagopal, Mech. Res. Comm.17, 415 (1990).

[33] P. D. Weidman, D. G. Kubitschek, and A. M. J. Davis, Int. J. Eng. Sci.44, 730 (2006).

[34] A. Postelnicu and I. Pop, Appl. Math. Comput. 217, 4359 (2011).

[35] A. V. Ro¸sca and I. Pop, Int. J. Heat Mass Trans.60, 355 (2013).

[36] N. C. Ro¸sca and I. Pop, Int. J. Heat Mass Trans.65, 102 (2013).

[37] S. D. Harris, D. B. Ingham, and I. Pop, Transport Porous Med.77, 267 (2009).

[38] M. M. Rahman, M. A. Al-Lawatia, I. A. Eltayeb, and N. Al-Salti, Int. J. Thermal Sci.57, 172 (2012).

[39] M. M. Rahman, A. V. Ro¸sca, and I. Pop, Int. J. Heat Mass Trans.77, 1133 (2014).

[40] L. F. Shampine, I. Gladwell, and S. Thompson, Solving ODEs with Matlab, Cambridge University Press, Cam- bridge 2003.

[41] L. F. Shampine, M. W. Reichelt, and J. Kierzenka, Solving boundary value problems for ordinary differ- ential equations in MATLAB with bvp4chttp://www.

mathworks.com/bvp_tutorial, 2003.