Acta Mech DOI 10.1007/s00707-015-1334-2

O R I G I NA L PA P E R

M. M. Rahman · J. H. Merkin · I. Pop

Mixed convection boundary-layer flow past a vertical flat plate with a convective boundary condition

Received: 14 November 2014 / Revised: 13 January 2015 © Springer-Verlag Wien 2015

Abstract The mixed convection boundary-layer flow on a vertical surface with an applied convective boundary condition is considered. Specific forms for the outer flow and surface heat transfer parameter are taken to reduce the problem a similarity system, which is seen to involve three parameters: m, the exponent of the outer flow; λ, the mixed convection parameter and B, the Biot number, as well as the Prandtl number. m = 15 is found to be a transitional case with different behaviour depending on whether m > 15 or m < 15 . For m > 15 , there is a critical value λc with solutions only for λ ≥ λc , a range of values of λ where there are dual solutions and the upper solution branch continuing into aiding flow. For 0 < m < 15 , there is a value of B where there is a change in behaviour from that seen for m > 15 to solutions terminating at a finite value of λ > 0 and continuing to large values of |λ|. For m = 0 (uniform outer flow) only this latter behaviour is seen to arise. When m < 0, though still in the range when there is a solution to the Falkner–Skan, λ = 0, problem, the behaviour is similar to that seen for m > 15 . 1 Introduction Mixed convection flows, or combined forced and free convection boundary-layer flows, arise in many transport processes both naturally and in engineering applications. They play an important role, for example, in atmospheric boundary-layer flows, heat exchangers, solar collectors, nuclear reactors and electronic equipment. Such processes occur when the effects of buoyancy forces in forced convection or the effects of forced flow in free convection become significant. According to Chen and Armaly [1], the domain of the mixed convection regime is generally defined as the region, a ≤ Gr/Ren ≤ b, where Gr is the Grashof number, Re is the Reynolds number, a and b are, respectively, the lower and the upper bounds of the domain, and n is a constant, which depends on the flow configuration and the surface heating conditions. The buoyancy parameter Gr/Ren provides a measure of the influence of free convection in comparison with that of forced convection on the flow. Outside the mixed convection region a ≤ Gr/Ren ≤ b, either a pure forced convection or a pure M. M. Rahman Department of Mathematics and Statistics, College of Science, Sultan Qaboos University, P.O. Box 36, P.C. 123 Al-Khod, Muscat, Sultanate of Oman E-mail: [email protected] J. H. Merkin (B) Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK E-mail: [email protected] I. Pop Department of Mathematics, Babe¸s-Bolyai University, 400084 Cluj-Napoca, Romania E-mail: [email protected]

M. M. Rahman et al.

free convection analysis can be used to describe accurately the flow and the temperature fields. The interaction between forced and free convection is especially pronounced in situations where the forced flow velocity is low and/or the temperature differences are large. Buoyancy forces also play a significant role in instigating flow instabilities. They can be responsible for either delaying or speeding up the transition from laminar to turbulent flow. Many analyses of mixed convection flow of a viscous and incompressible fluid over a vertical surface have already been performed. Analytical and numerical solutions for the temperature and the velocity fields have been obtained both for prescribed wall temperatures and for prescribed wall heat fluxes. A comprehensive description of the theoretical work prior to 1987 for both laminar and turbulent mixed convection boundarylayer flows for selected flow geometries has been given in a review paper by Chen and Armaly [1] and in the books by Gebhart et al. [2], Schlichting and Gersten [3], Pop and Ingham [4] and Bejan [5], for example. A wide range of geometrical configurations and heating conditions has been treated for both natural convection and cases where there is an interaction between the buoyancy-induced convection and some external imposed flow (mixed convection). The majority of these previous studies have assumed either some prescribed wall temperature or some prescribed wall heat flux distribution. However, there is a class of boundary-layer flows and heat transfer problems in which the surface heat flux depends on the local surface temperature. Perhaps the simplest case of this is when there is a linear relation between the surface heat transfer and surface temperature. This situation arises in conjugate heat transfer problems, see, for example, Merkin and Pop [6], and also when there is Newtonian heating of the convecting fluid from the surface, this latter case being discussed in a detail by Merkin [7]. More recently, heat transfer problems involving convective boundary conditions have been investigated by Aziz [8], Ishak [9] and Magyari [10] for Blasius flow. A similar analysis was applied to Blasius flow and Sakiadis (stretching) flow with radiation effects by Bataller [11]. Yao et al. [12] have investigated the heat transfer of a viscous fluid flow over a permeable stretching/shrinking sheet with a convective boundary condition. It was found in [12] that the convective boundary condition results in a temperature slip at the wall and this temperature slip is greatly affected by the mass transfer parameter, the Prandtl number and the parameters associated with the stretching/shrinking of the bounding surface. It was also shown in [12] that the temperature profiles within the fluid are considerably different to those seen when there is a prescribed wall temperature. In a large number of problems in fluid mechanics and heat transfer, the equations expressing the basic physical laws are partial differential equations, with the aim being to find exact (or numerical) solutions to these equations. Solutions of certain sets of partial differential equations occurring in applied fields can be found quite readily in spite of the failure of the more common classical methods to yield results. Notable among such solutions are those that have been obtained by employing transformations of the variables that reduce the system of partial differential equations to a system of ordinary differential equations. These solutions are generally designated as similarity solutions, see Hansen [13]. Such similarity solutions for the problem under consideration here have recently been presented in a series of papers by Merkin and Mahmood [14], Mahmood and Merkin [15,16], Rhida [18] and Merkin and Pop [17] for the steady mixed convection boundary-layer flow over an impermeable vertical flat plate and for impermeable vertical cylinders. The forced convection heat transfer resulting from the flow of a uniform stream over a flat surface on which there is a convective boundary condition has been recently analysed by Merkin and Pop [19]. In the previous papers [8–10], it was assumed that the convective heat transfer parameter h f associated with the hot surface depended on x, where x measures distance along the surface, so that problem could be reduced to similarity form. In [19], it was assumed that this heat transfer parameter h f is a constant, with the result that the temperature profiles and overall heat transfer characteristics evolve as the solution develops from the leading edge. We mention here also the paper by Makinde and Olanrewaju [20] on buoyancy effects on thermal boundary layer over a vertical plate with a convective surface boundary condition. Pantokratoras [21], in a recently published paper, has studied the steady laminar boundary-layer flow along a vertical stationary wall with convective surface boundary condition. The heat transfer coefficient is considered either constant or variable along the plate, and the problem is either non-similar or similar. The results were obtained with the direct numerical solution of the governing equations. It was found that there are differences in the results between the non-similar and the similar case at low values of the convective parameter, but as this parameter increases the differences decrease and the flow tends to the classical natural convection along a vertical isothermal plate. The present analysis is an extension of this study by Pantokratoras [21] aiming to investigate the role of the governing parameters in the boundary-layer flow and heat transfer from a heated or cooled vertical semi-infinite flat surface in the mixed convection regime. In order to reduce the basic partial differential equations to a similarity form, we assume that the external (inviscid) velocity Ue (x) and convective

Mixed convection boundary-layer flow

temperature Tf (x) take specific functional forms. The resulting system of ordinary differential equations subject to the corresponding boundary conditions is solved numerically using the function bvp4c from MATLAB for representative values of the dimensionless parameters. Since there is almost no experimental data, the choice of the values of the parameters was dictated by the values chosen by previous investigators.

2 Equations We consider the steady, mixed convection boundary-layer flow over a vertical flat surface in a viscous and incompressible fluid. The surface is located at y = 0, where y is the coordinate measured in the direction normal to surface. We assume that the velocity of the external (inviscid) flow is Ue (x) and is at a constant temperature T∞ , where x is the coordinate measured along the surface. We assume that the flow takes place on the right-hand side of the surface with the left-hand side of the surface being heated or cooled by convection from a hot or cold fluid bath at temperature Tf (x), providing a heat transfer coefficient h f (x). Here, Tf (x) > T∞ corresponds to assisting flow, and Tf (x) < T∞ corresponds to opposing flow. Under these assumptions and the standard Boussinesq approximation, the governing equations are described in the Cartesian coordinates x and y as ∂u ∂v + = 0, ∂x ∂y u

(1)

dUe ∂u ∂u ∂ 2u +v = gβ(T − T∞ ) + Ue +ν 2, ∂x ∂y dx ∂y

(2)

∂T ∂T ∂2T +v =α 2. ∂x ∂y ∂y

(3)

u

Equations (1–3) are subject to the boundary conditions v = 0, u = 0, −k

∂T = h f (Tf − T ) ∂y

on y = 0, u → Ue (x), T → T∞

as y → ∞.

(4)

Here, u and v are the velocity components along the x and y axes, T is the fluid temperature within the boundary layer, g is the acceleration due to gravity, α is the thermal diffusivity of the fluid, β is the coefficient of thermal expansion, k is the thermal conductivity of the surface and ν is the kinematic viscosity of the fluid. Our aim here is to reduce Eqs. (1–4) to similarity form. In order to do so, we assume that Ue and Tf take particular functional forms, namely Ue (x) = ax m , Tf (x) − T∞ = T0 x 2m−1 for some exponent m, where a is a positive constant and T0 is a characteristic temperature scale, with T0 > 0 for assisting flow (heated surface) and T0 < 0 for opposing flow (cooled surface). We then write, for m  = −1,  ψ=

2νa m+1

 1/2  (m + 1)a 1/2 (m−1)/2 T − T∞ x (m+1)/2 f (η), = θ (η), η = y x , Tf − T∞ 2ν

(5)

where ψ is the stream function defined in the usual way. Applying transformation (5) in Eqs. (1–3) leads to 2m (1 − f 2 ) = 0, m+1 1  2(2m − 1)  θ + f θ − f θ = 0, σ m+1

f  + λθ + f f  +

(6) (7)

subject to the boundary conditions that f (0) = 0, f  (0) = 0, θ  (0) = −B [1 − θ (0)],

f  → 1, θ → 0 as η → ∞,

(8)

where primes denote differentiation with respect to η. In order that boundary conditions (8) are independent of x, and hence, the problem is reduced fully to similarity form, we need a specific form for h f (x), namely h f (x) = h 0 x (m−1)/2

where h 0 > 0 is a constant.

(9)

M. M. Rahman et al.

Here, σ = ν/α is the Prandtl number, λ is the mixed convection parameter and B is the Biot number (or surface convection parameter), defined as  1/2 gβT0 2 h0 2ν λ= . (10) , B = a2 m + 1 k a(m + 1) We note that λ > 0 corresponds to assisting flow, λ < 0 to opposing flow and λ = 0 gives forced convection. The quantities of particular interest are the skin friction coefficient Cf and the surface temperature Tw , defined by   μ ∂u Cf (x) = , Tw (x) − T∞ = (Tf − T∞ ) θ (x, 0) = T0 x 2m−1 θw . (11) 2 ρUe ∂ y y=0 From (5), we have

 Cf =

ν(m + 1) 2a

1/2

x −(m+1)/2 f  (0), θw = 1 +

θ  (0) , B

(12)

the latter expression relating the surface heat flux −θ  (0) to the wall temperature θw . From (12), we see that, when B → ∞, we recover the classical solution for a prescribed surface temperature. 3 Numerical technique Following Ro¸sca and Pop [23] and Rahman et al. [24–28], the system of ordinary differential equations (6– 7) subject to boundary conditions (8) was solved numerically using the function bvp4c from MATLAB for representative values of the dimensionless parameters. Since there is almost no experimental data, the choice of the values of the parameters 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 descriptions of the respective figures. The code bvp4c is developed using a finite-difference method that implements the three-stage Lobatto IIIa formula. This is a collocation method with fourth-order accuracy. In this approach, the ordinary differential equations (6–7) are first reduced to a system of first order by introducing new variables. The mesh selection and error control are based on the residual of the continuous solution. The relative error tolerance was set at 10−7 . Since the present problem may have more than one solution (dual, upper and lower branch solutions), a ‘good’ initial guess is necessary. The ‘infinity’ condition in boundary conditions (8) is implemented at a finite value η = η∞ . We started the computations at small value of η∞ , for example, η∞ = 20.0, and then subsequently increased the value of η∞ until the outer boundary conditions are satisfied to sufficient accuracy. In our approach, we chose a suitable finite value of η∞ = 20.0 for the upper branch (first) solutions and η∞ in the range 40–50 for the lower branch (second) solutions. Examples of solving boundary-value problems using bvp4c can be found in the book by Shampine et al. [29] or through online tutorial by Shampine et al. [30]. 4 Results For forced convection λ = 0, Eq. (6) with the boundary conditions for f in (8) is the classical Falkner–Skan wedge flow [3,22,31–33]. This problem is usually characterised by the parameter β = 2m/(m + 1) rather than m. In this case, see [32] in particular, there is a critical value βc of β, with βc  −0.19884 (or m c  −0.09043), with no solutions for β < βc , dual solutions for βc < β < 0 and only a single solution for β ≥ 0. In Fig. 1a, we plot the corresponding values of θ  (0) against m, and the existence of a critical value at m = m c and dual solutions for m > m c are not particularly clear in the figure. Perhaps the more interesting feature is the existence of singularity in the solution at m = m 0  0.1058, shown in the figure by a broken line, where θ  (0) becomes infinite. The Falkner–Skan solutions pass smoothly through m 0 with f  (0)  0.510463 at this value. The nature of the singularity as m → m 0 is discussed in [35], where we put m = m 0 + , θ = −1 h and look for a solution valid for small . This results in an eigenvalue problem for m 0 at leading order, namely 2m 0 (1 − f 2 ) = 0, m0 + 1 f (0) = f  (0) = 0, h  (0) = B h(0),

f  + f f  +

1  2(2m 0 − 1)  h + f h − f h = 0, σ m0 + 1 f  → 1 h → 0 as η → ∞.

(13) (14)

Mixed convection boundary-layer flow

(a) 6 4

θ (0)

2

0

-2

-4

-6

(b)

0

0.2

0.4

0.6

m

0.8

1

8

10

0.2 0.18

m0

0.16 0.14 0.12 0.1 0.08 0.06

0

2

4

6

B Fig. 1 Forced convection, λ = 0: a a plot of θ  (0) against m obtained from the numerical solution of (6–8) for B = 1.0. The singularity at m = m 0  0.1058 is shown by the broken line. b A plot of m 0 against B obtained from the numerical solution of Eqs. (13, 14), the asymptotic limit for B large given in [35] is shown by a broken line

We also impose the additional condition h(0) = 1 to force a nontrivial solution for h. A graph of m 0 against B obtained from the numerical solution of (13, 14) is plotted in Fig. 1b. For B large, the boundary condition reduces to θ (0) = 1, the case treated in [35], giving m 0  0.07072, shown in Fig. 1b by a broken line. When B = 0, h  (0) = 0 giving m 0 = 15 . The existence of a critical value for this case suggests that there could well be critical values in the more general problem (6–8) and this feature motivates our discussion. Numerical solutions to Eqs. (6, 7) subject to boundary conditions (8) were obtained for various values of the physical parameters λ, m and the Biot number B taking the Prandtl number σ = 1 throughout. We concentrate mainly on obtaining the conditions under which dual (upper and lower branch) solutions may exist, i.e. in the opposing flow regime, when a convective boundary condition applies on the surface. We start by noting that, when m = 15 , Eq. (7) can be integrated to get, on satisfying the outer boundary conditions, θ  + σ f θ = 0.

(15)

From (15), θ  (0) = 0, and hence, from boundary conditions (8), θ (0) = 1. Thus, for this value of m, the solution is independent of the Biot number B. In Fig. 2, we plot f  (0) against λ, noting that θ (0) ≡ 1. This

M. M. Rahman et al.

0.9 0.8 0.7

f (0)

0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1

-0.5

-0.4

-0.3

-0.2

-0.1

0

λ Fig. 2 A plot of the skin friction f  (0) against λ for the transitional case m = 15 with Prandtl number σ = 1.0, obtained from the numerical solution of (6–8). For this value of m, θ(0) ≡ 1 and the solution is independent of B

figure indicates a critical value λc  −0.5189 of λ, a range of negative λ > λc where there are dual solutions. The upper solution branch continues to large positive values of λ. The critical value λc is, following [34,36], in this case determined by nontrivial solution to the homogeneous system φ  + λc h + f c φ  + f c φ −

2   f φ = 0, 3 c

1  h + f c h + θc φ = 0, σ

(16)

obtained by making a linear perturbation (φ, h) to Eqs. (6, 7) subject to φ(0) = φ  (0) = 0, h(0) = 0, φ  → 0, h → 0

as η → ∞, φ  (0) = 1,

(17)

where ( f c , θc ) is the solution to (6–8) at the critical value λc and the final condition in (17) is applied to force a nontrivial solution. Equations (16) have the solution φ = A f  , h = A θ  for some constant A  = 0. This solution satisfies Eq. (16) and boundary conditions (17) provided f  (0) = 0. Thus, for this case, we can find the critical value by solving Eqs. (6, 7) with m = 15 subject to (8) and the additional constraint that f  (0) = 0. In our numerical integration to obtain Fig. 2, we noted that f  (0) passed through zero as the solution passed through the critical point. We find that this value of m = 15 is a transitional value with different behaviour of the solution depending on whether m > 15 or m < 15 . We treat these cases separately starting by taking m > 15 . 4.1 m>1/5 In Figs. 3 and 4, we plot the reduced skin friction coefficient f  (0), as defined in (12), and the surface temperature θw against λ for m = 1 (planar stagnation-point flow) and m = 0.5 (Tf being constant) taking B = 4.0 (Fig. 3) and B = 0.5 (Fig. 4). These figures show that, for opposing flow, λ < 0, there is a critical value λc of λ limiting the range of existence of a similarity solution to λc ≤ λ and a range of λ where there are dual solutions for λc < λ. The value of λc depends on both m and B, increasing as m is increased (for a given B) and decreasing as B is increased (for a given m). The upper solution branch, as characterised in Fig. 3a or Fig. 4a, continues to large positive values of λ (not really seen in these figures), whereas the lower branch solutions appear to terminate at a finite negative value of λ with θw becoming large. The lower branch solutions terminate at a finite, nonzero value λ0 of λ with f  (0) approaching zero from below and θ (0) approaching finite, nonzero values. This is, perhaps, best seen in Fig. 2 and the plots in Figs. 3 and 4 for m = 0.5. These values of λ0 depend on both B and m.

Mixed convection boundary-layer flow

(a) 1.4 1.2 1 0.8

f (0)

0.6

1.0

0.4

0.5

0.2 0 -0.2 -0.4 -0.6 -0.8

-2.5

-2

-1.5

-1

-0.5

0

-0.5

0

λ

(b)

1.1

1.05

θw

1

0.95

0.5

0.9

1.0 0.85

0.8

-2.5

-2

-1.5

-1

λ Fig. 3 Plots of a the skin friction f  (0) and b the wall temperature θw against λ for B = 4.0 and m = 1, 0.5 and Prandtl number σ = 1.0, obtained from the numerical solution of (6–8)

In Fig. 5, we give velocity f  and temperature θ plots representative of the opposing flow regime when there are dual solutions, here for λ = −1.0, B = 4.0, m = 1.0. The lower branch solution is shown by a broken line. The upper branch solution is monotone, with f  increasing from zero at the wall, where f  (0) = 0.8159, to its outer value and θ decreasing from its wall value of θw = 0.8428 to zero at large distances. However, the lower branch solutions start with f  decreasing, here f  (0) = −0.3446, before attaining a minimum value of approximately −0.2850 and then increasing to its outer value. The reverse is the case for θ , which starts by increasing θw = 1.0704, θ  (0) = 0.2816, reaching a maximum of approximately 1.3094 before decreasing to zero. The boundary-layer thickness is greater for the lower branch solutions. These plots are typical of the velocity and temperature profiles for opposing flow when m > 15 . We now consider the critical values λc of λ, seen in Figs. 3 and 4, in more detail. 4.1.1 Critical values An important feature of our results for opposing flow is the existence of a critical value λc of the mixed convection parameter λ limiting the range of existence of a similarity solution, with the values of λc being dependent on both m and B (as well as the Prandtl number). We can calculate λc numerically using the approach described in [34,36], see also Eqs. (16, 17). In Fig. 6, we plot λc against B for m = 0.5, 1.0, and the values of

M. M. Rahman et al.

(a) 1.4 1.2 1

1.0

0.8

0.5

f (0)

0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -5

-4

-3

-2

-1

0

1

0

1

λ

(b) 1.1 1 0.9

θw

0.8 0.7 0.6 0.5

0.5 1.0

0.4 0.3 -5

-4

-3

-2

-1

λ Fig. 4 Plots of a the skin friction f  (0) and b the wall temperature θw against λ for B = 0.5 and m = 1, 0.5 and Prandtl number σ = 1.0, obtained from the numerical solution of (6–8)

λc for B = 0.5 and B = 4.0 agree well with those seen in Figs. 3 and 4. The figure shows that λc < 0 for all B and that λc increases as B increases, approaching a finite, nonzero limit for B large. We have already noted that, as B → ∞, we recover the prescribed surface temperature, θ (0) = 1, case. For this case, λc  −1.1936 for m = 1.0 and λc  −0.5357 for m = 0.5. These values are shown in Fig. 6 by broken lines and are the values being approached by λc as B becomes large. Figure 6 suggests that λc becomes large and negative as B decreases to zero. To describe this behaviour, we put θ = B θ. This leaves Eqs. (6, 7) unaltered except that now, the buoyancy term in (6) becomes (B λ) θ . Boundary conditions (8) give θ  (0) = −1+ B θ (0), becoming θ  (0) = −1 as B → 0, leading to the prescribed wall heat flux case at leading order. By calculating the critical values for this case, we have λc ∼ −2.3618 B −1 + · · · for m = 1.0,

λc ∼ −1.4847 B −1 + · · · for m = 0.5.

(18)

In Fig. 7, we plot λc against m for B = 0.5 and B = 4.0. We see again that λc < 0 for all m and that the values of λc approach a constant value, dependent on B, as m increases, more clearly seen for B = 4.0. Also, both curves start at m = 0.2 with the same value of λc  −0.5189 quoted above. For m large, Eqs. (6, 7) reduce to f  + λ θ + f f  + 2(1 − f 2 ) = 0,

1  θ + f θ  − 4 f  θ = 0, σ

(19)

Mixed convection boundary-layer flow

(a)

1 0.8 0.6

f

0.4 0.2 0 -0.2 -0.4

0

1

2

3

0

1

2

3

η

4

5

6

4

5

6

(b) 1.4 1.2 1

θ

0.8 0.6 0.4 0.2 0

η

Fig. 5 a Velocity and b temperature profile plots in the opposing flow regime for λ = −1.0, m = 1.0 and B = 4.0, the lower branch solution is shown by a broken line, obtained from the numerical solution of (6–8)

still subject to boundary conditions (8). The solution to Eqs. (19, 8) has a structure similar to that seen in Figs. 3 and 4; in particular, there is again a critical value λc , dependent on B. We can calculate these critical values as for the general case, finding that λc → −5.7090 for B = 4.0,

λc → −12.1196 for B = 0.5 as m → ∞

(20)

consistent with figure 7. In figure 8, we plot λc against B for the large m limit (19). The figure shows that λc approaches a finite, negative value of approximately −4.969 (large B limit) as B increases and becomes increasingly larger as B is decreased, appearing to become infinite as B approaches 0.2. We now consider the case when m < 15 . 4.2 m<1/5 For this case, we start by considering the results for m = 0.1. In Fig. 9, we plot f  (0) and θw against λ for B = 4.0 and for B = 1.0. For B = 4.0, the behaviour of the solution is similar to that seen previously for m > 15 , in that there is a solution for all the aiding flow, λ > 0, regime and critical value λc  −0.0273 limiting solutions in the opposing flow regime to λ ≥ λc (though this value of λc is somewhat smaller than those seen

M. M. Rahman et al.

-1

0.5 -2 -3

1.0

λc

-4 -5 -6 -7 -8 -9 -10

0

1

2

3

4

5

B Fig. 6 Plots of the critical values λc of the mixed convection parameter λ against B for m = 1.0, 0.5. The asymptotic values for large B are shown by the broken lines 0

-2

0.5

λc

-4

-6

-8

4.0

-10

-12

0

2

4

6

m

8

10

12

14

Fig. 7 Plots of the critical values λc of the mixed convection parameter λ against m for B = 4.0, 0.5

previously). However, for B = 1.0, the nature of the solution is very different. Here, there is no limit to the solution for opposing flow (note that we obtained our numerical results for very much larger negative values of λ than are plotted in Fig. 9). The solution now terminates in aiding flow at a finite value of λ  0.28195. A significant difference between the solutions for B = 4.0 and B = 1.0 is that, in the latter case, θw < 0 in general, becoming positive only at λ  0.2306 giving a very small range of λ where θw > 0. For B = 4.0, θw > 0 throughout. We were unable to continue our numerical solutions past this point and did not detect a critical value or dual solutions. We next consider the case when m = 0 (uniform outer flow) giving plots of f  (0) and θw against λ in Fig. 10 for B = 4.0 and B = 0.5. For this value of m, we see behaviour similar to that seen previously for m = 0.1, B = 1.0, in that the solution terminates at a finite, positive value λ0 of λ, dependent on B with λ0  0.1790 for B = 0.5 and λ0  0.0781 for B = 4.0, and continues to large, negative values of λ, i.e. there is a finite limit to the range of λ for solutions for aiding flows and no limit on λ for opposing flows. We now consider a negative value for m, limiting attention to the situation when there is a forced convection, λ = 0 solution, i.e. to m > m c = −0.09043, taking m = −0.05 as a representative case. In Fig. 11, we show plots of f  (0) and θw against λ for B = 0.5 and B = 4.0. We again find critical values and dual solutions in

Mixed convection boundary-layer flow

-4 -6 -8 -10

λc

-12 -14 -16 -18 -20 -22 0

1

2

3

4

5

6

7

8

B Fig. 8 A plot of the critical values λc of the mixed convection parameter λ against B for the large m limit given by Eqs. (19, 8)

the opposing flow regime, limiting the range of existence of possible solutions, with the upper solution branch continuing to large, positive values of λ. This situation is similar to that seen for m > 15 , comparing Fig. 11 with Figs. 3 and 4. A feature of the problem for m < 15 is the possibility of having solutions for λ < 0 and |λ|  1, as in Figs. 9 and 10. We derive a solution in this limit by putting λ = −μ, μ > 0 and looking for a solution for μ large by writing f = (B μ)1/5 F, θ = B 4/5 μ−1/5 H, ζ = (B μ)1/5 η, (B > 0).

(21)

Applying transformation (21) in Eqs. (6–8) gives 2m 2m F 2 + (B μ)−4/5 = 0, m+1 m+1 1  2(2m − 1)  H + FH − F H = 0, σ m+1

F  − H + FF −

(22)

subject to F(0) = F  (0) = 0, H  (0) = −1 + B 4/5 μ−1/5 H (0), F  → (B μ)−2/5 , H → 0

as ζ → ∞.

(23)

where primes now denote differentiation with respect to ζ . The leading-order problem is, from (22, 23), F  − H + FF −

2m 1  2(2m − 1)  F 2 = 0, H + FH − F H = 0, m+1 σ m+1

(24)

subject to F(0) = F  (0) = 0, H  (0) = −1,

F  , H → 0 as ζ → ∞.

(25)

We plot F  (0) and H (0) obtained from the numerical solution of (24, 25) in Fig. 12. We see that F  (0) > 0 and that H (0) changes from positive to negative at m  0.01. Also, both F  (0) and |H (0)| become large as m → 0.2. We could expect a singularity as m → 15 for, with this value of m, Eq. (24) integrates to H  + σ F H = c0 , where c0 is a constant which cannot be chosen so as to satisfy both the outer boundary condition and H  (0) = −1. Our numerical solution terminated at a finite value m  −0.072 where F  (0) → 0 and H (0) approached a positive limit in a singular way (though not obvious from the figure), thus putting a limit on the range of m for this asymptotic solution to apply. From (21) f  (0) ∼ B 3/5 |λ|3/5 F  (0), θw ∼ B 4/5 |λ|−1/5 H (0) as λ → −∞,

(26)

M. M. Rahman et al.

5

(a)

4.5 4

1.0

f (0)

3.5

4.0

3 2.5 2 1.5 1 0.5 0 -2

-1.5

-1

-0.5

0

0.5

1

1.5

2

1.5

2

λ

(b) 6 5 4

θw

3

4.0

2 1 0 -1

1.0

-2 -3 -2

-1.5

-1

-0.5

0

0.5

1

λ Fig. 9 Plots of a the skin friction f  (0) and b the wall temperature θw against λ for m = 0.1 and B = 4.0, 1.0 and Prandtl number σ = 1.0, obtained from the numerical solution of (6–8)

where the values of F  (0) and H (0) depend on m and are given in Fig. 12. Thus, the values f  (0) increase more rapidly with |λ| for larger values of B and that θw decreases to zero slowly as |λ| is increased, consistent with the plots shown in Fig. 10. The results presented above suggest that we again examine the critical values λc . 4.2.1 Critical values We first consider the range 0 < m < 15 . In Fig. 13, we plot the values of λc against m for B = 4.0, and in Fig. 14, the values of λc are plotted against B for m = 0.1. Figure 13 shows that λc approaches the previously quoted value of λc  −0.5189 as m → 0.2. For this value of B, the graph terminates at m = m 0  0.0845 in a singular way, with λc → 0 and with, at λ = λc , the values of θ (0) → ∞ and f  (0) approaching a finite, nonzero limit. Figure 14 indicates that λc approaches a finite (prescribed temperature) limit as B → ∞ and that, for this particular value of m, λc → 0 at B  1.74 with f  (0) approaching a finite, nonzero limit of 0.6658 and θ (0) becoming increasingly larger. These results, together with the plots in Fig. 10 for B = 1.0, suggest that, at least for 0 < m < 15 , there are values of the parameters m, B, λ where there is a transition from having a solution for all aiding flow and only a limited range for opposing flow to having a solution for all opposing flow and only a limited range for aiding flow.

Mixed convection boundary-layer flow

(a) 2.5

4.0

f (0)

2

0.5

1.5

1

0.5 -10

-8

-6

-4

-2

0

-2

0

λ

(b) 0.9 0.8 0.7

4.0

θw

0.6 0.5 0.4 0.3

0.5

0.2 0.1 -10

-8

-6

-4

λ Fig. 10 Plots of a the skin friction f  (0) and b the wall temperature θw against λ for m = 0 and B = 4.0, 0.5 and Prandtl number σ = 1.0, obtained from the numerical solution of (6–8)

The critical values λc for m = −0.05 are given in Fig. 15 plotted against B. We see that λc approaches a nonzero, finite limit (the prescribed wall temperature value of −0.3199, shown in the figure by a broken line) as B → ∞ and becomes large and negative as B → 0. The picture here is similar in many respects to the plots of λc shown in Fig. 6. 4.3 Free convection limit, solution for λ large So far, we have been concerned mainly with solutions for opposing flow, concentrating on critical values and dual solutions. We now briefly examine aiding flow, λ > 0, in particular the solution for λ large, the free convection limit. To get a solution valid for λ large, we follow closely that derived above in Eqs. (21–24) for λ < 0, |λ|  1. This leads us to make the transformation, on assuming that B > 0 and is of O(1), as in (21) f = (λ B)1/5 F, θ = λ−1/5 B 4/5 H, ζ = (λ B)1/5 η.

(27)

When we apply (27) in Eqs. (6–8), we obtain, at leading order, F  + H + FF −

2m F 2 = 0, m+1

1  2(2m − 1)  H + FH − F H = 0, σ m+1

(28)

M. M. Rahman et al.

(a) 0.3

f (0)

0.2

4.0 0.5

0.1

0

-0.1

-0.2 -0.6

-0.5

-0.4

-0.3

λ

-0.2

-0.1

0

0.1

-0.1

0

0.1

(b) 0.95 0.9

4.0

0.85

θw

0.8 0.75 0.7 0.65

0.5

0.6 0.55 0.5 -0.6

-0.5

-0.4

-0.3

-0.2

λ Fig. 11 Plots of a the skin friction f  (0) and b the wall temperature θw against λ for m = −0.05 and B = 4.0, 0.5 and Prandtl number σ = 1.0, obtained from the numerical solution of (6–8)

subject to F(0) = F  (0) = 0, H  (0) = −1,

F  → 0, H → 0

as ζ → ∞,

(29)

where primes now denote differentiation with respect to ζ . In Fig. 16, we plot F  (0) and H (0) against m obtained from a numerical integration of (28, 29). Both F  (0) and H (0) become large as m → 15 , again suggesting that m = 15 is a transitional value of m. Compare this figure with Fig. 12 for the λ < 0, |λ|  1 limit. From (27) f  (0) ∼ λ3/5 B 3/5 F  (0) + · · · , θ (0) ∼ λ−1/5 B 4/5 H (0) + · · · as λ → ∞.

(30)

In Fig. 17, we plot λ−3/5 f  (0) and λ1/5 θ (0) against λ obtained from solving (6–8) numerically with m = 1.0 and B = 4.0, 0.5. The respective large λ limits are 0.4395 and 4.5918 for B = 4.0 and 1.5304 and 0.8700 for B = 0.5. We see that these asymptotic limits are being approached extremely slowly, more clearly for B = 0.5 rather than for B = 4.0. This might be expected as the corrections are relatively large, only of O((λ B)−2/5 ), thus requiring very much larger values of λ before these limits are attained (as was borne out in our numerical integrations of (6–8) for larger λ).

Mixed convection boundary-layer flow

(a) 14 12

F (0)

10 8 6 4 2 0

-0.05

0

0.05

0.1

0.15

0.2

-0.05

0

0.05

0.1

0.15

0.2

m

(b) 0

H(0)

-5

-10

-15

-20

-25

-30

m

Fig. 12 Plots of a F  (0) and b H (0) against m for the large μ, large negative λ, limit given by Eqs. (24, 25) 0 -0.05 -0.1

λc

-0.15 -0.2 -0.25 -0.3 -0.35 -0.4 -0.45 0.08

0.1

0.12

0.14

m

0.16

0.18

0.2

Fig. 13 A plot of the critical values λc of the mixed convection parameter λ for B = 4.0 against m in the 0 < m <

1 5

range

M. M. Rahman et al.

0 -0.005 -0.01 -0.015

λc

-0.02 -0.025 -0.03 -0.035 -0.04 -0.045 -0.05

1

2

3

4

5

6

7

8

9

B Fig. 14 A plot of the critical values λc of the mixed convection parameter λ against B for m = 0.1 -0.2 -0.4 -0.6

λc

-0.8 -1 -1.2 -1.4 -1.6

0

0.5

1

1.5

2

2.5

3

3.5

4

B Fig. 15 A plot of the critical values λc of the mixed convection parameter λ against B for m = −0.05. The asymptotic value for large B is shown by the broken line

4.3.1 The limit as m → 1/5 To determine the nature of the singularity seen in Fig. 16 as m → 15 , we put m = and look for a solution valid for δ small by putting F = δ −1/5 g,

1 5

+ δ, where 0 < δ 1,

H = δ −4/5 h, ξ = δ −1/5 η.

(31)

Equations (28, 29) become  1 25δ 1+ + · · · g 2 = 0, 3 6   1  25δ  h +gh + 1− + · · · g  h = 0, σ 6

g  + h + g g  −

(32)

where primes now denote differentiation with respect to ξ , subject to g(0) = g  (0) = 0, h  (0) = −δ,

g  → 0, h → 0

as ξ → ∞.

(33)

Mixed convection boundary-layer flow

14 12 10 8 6 4

H(0) 2

F (0)

0 0.2

0.3

0.4

0.5

0.6

0.7

m

0.8

0.9

1

Fig. 16 The free convection, λ  1 limit: plots of a F  (0) and b H (0) against m obtained from the numerical solution of (28, 29)

Equations (32, 33) suggest looking for a solution by expanding g(ξ ; δ) = g0 (ξ ) + δ g1 (ξ ) + · · · , h(ξ ; δ) = h 0 (ξ ) + δ h 1 (ξ ) + · · · .

(34)

The leading-order problem is, on integrating the equation for h 0 and satisfying the boundary conditions, 1 2 g , h 0 + σ g0 h 0 = 0, 3 0 g0 (0) = g  (0) = h 0 (0) = 0, g0 → 0, h 0 → 0 as ξ → ∞.

g0 + h 0 + g0 g0 −

(35)

Equation (35) is a homogeneous problem. Now suppose (g 0 , h 0 ) is a solution to this problem with g 0 (0) = 1. Numerical integration for this case gives h 0 (0) = 0.87606 for σ = 1. In general, g0 (0) = a0 for some constant a0  = 0 and the general solution (g0 , h 0 ) can be obtained from (g 0 , h 0 ) by the transformation 1/3

g0 = a0

4/3

g 0 , h 0 = a0

1/3

h 0 , ξ = a0

ξ.

(36)

At O(δ), we need only to consider the equation for h 1 , which is 1  25  h 1 + g0 h 1 + g1 h 0 + g1 h 0 + g0 h 1 = g h0, σ 6 0

(37)

subject to g1 (0) = g1 (0), h 1 (0) = −1,

g1 → 0, h 1 → 0 as ξ → ∞.

(38)

Equation (37) can be integrated to give, on satisfying the outer boundary condition, 25 1  h 1 + g0 h 1 + g1 h 0 = − σ 6

 ξ

∞

g0 h 0 dξ

5/3

25 a0 =− 6

 ξ

∞

g 0 h 0 dξ ,

(39)

on using transformation (36). Applying the boundary condition on ξ = 0 then gives 5/3  ∞ 25 a0 1 g 0 h 0 dξ . (40) = σ 6 0 ∞ Our numerical integration gives, for σ = 1, 0 g 0 h 0 dξ = 0.73072 and hence a0 = 0.51272. The above analysis can be applied equally well in the λ < 0, |λ|  1 limit. The results shown in Fig. 12 also indicate

M. M. Rahman et al.

(a)

3

2.5

λ−3/5 f (0)

2

1.5

0.5 1

0.5

0

4.0

0

5

10

15

20

25

30

35

40

30

35

40

λ

(b) 1.6 1.4

4.0

λ1/5 θ(0)

1.2 1 0.8 0.6

0.5

0.4 0.2 0

0

5

10

15

20

25

λ Fig. 17 Plots of a λ−3/5 f  (0), b λ1/5 θ(0) for m = 1.0 and B = 0.5, 4.0 and σ = 1.0 obtained the numerical solution of (6–8)

a singularity as m → 15 . We now put m = 15 − δ, and at least to the order considered, we then need only to change the sign of h 0 with Eqs. (37, 39) staying the same. From (31, 36), we then have F  (0) ∼ 0.5127 δ −3/5 + · · · ,

H (0) ∼ ±0.3595 δ −4/5 + · · · as δ → 0,

(41)

where δ = |m + 15 |. Expressions (41) clearly show the nature of the singularity arising in the solution to (28, 29), or to Eq. (24), as m → 15 and are consistent with the values plotted in Figs. 12 and 16. 5 Conclusions We have considered the mixed convection boundary-layer flow on a vertical surface when there is a convective boundary condition relating the surface temperature to the surface heat flux. By taking specific forms for the outer flow and surface heat transfer parameter, we were able to reduce the problem to a system of similarity Eqs. (6–8). Apart from the Prandtl number (which we took as unity throughout), the problem involved three parameters: m, the exponent of the outer flow; λ, the mixed convection parameter and B, the Biot number, defined in (10).

Mixed convection boundary-layer flow

We identified m = 15 as a transitional value for the exponent m. For m > 15 , the behaviour of the solution is similar to that seen previously in mixed convection similarity systems, in that there is a critical value λc < 0 with solutions only for λ ≥ λc , a range of values of λ < 0 where there are dual solutions and the upper solution branch continuing into aiding flow to large values of λ, see Figs. 3 and 4. We calculated these critical values in terms of B, figure 6, showing that λc progressed smoothly from the prescribed surface heat flux case (small B) to the prescribed surface temperature case (large B), and in terms of m, Fig. 7, with the values of λc starting at the same value for m = 15 and decreasing to the corresponding large m limit, Fig. 8. Similar behaviour was seen when m = 15 , Fig. 2, but now the solution is independent of B and has zero wall heat flux. For 0 < m < 15 , there is a value B0 of B, dependent on m, where there is a change in behaviour from that seen for m > 15 to the case where the solution terminates at a finite value of λ in the aiding flow regime but continues to large values of |λ| in opposing flow, Figs. 9 and 13. The asymptotic solution for the large |λ| limit was derived, see particularly expression (26) and Fig. 12, appearing to terminate at a finite value m  −0.072 of m. For m = 0 (uniform outer flow), our numerical results indicate that only this latter behaviour arises, Fig. 10, with the solution continuing to large, negative values of λ. When m < 0, though still in the range when there is a solution to the Falkner–Skan, λ = 0, problem, the behaviour is again similar to that seen for m > 15 , Figs. 11 and 15, though the temperature distribution can be very different, Fig. 11b. One property of our similarity solutions, particularly of the lower branch solutions, is the existence of a region within the boundary layer where the fluid temperature is above that on the wall, see Fig. 5b for example. This appears to be an inherent feature of our present model (as well as other related models). Such a situation has been discussed in detail by Steinrück [37] where a justification for having above-wall temperatures in similarity solutions is presented as well as a discussion of the stability of such solutions.

References 1. Chen, T.S., Armaly, B.F. : Mixed convection in external flow. In: Kakaç, S., Shah, R.K., Aung, W. (eds.) Handbook of Single-Phase Convective Heat Transfer, pp. 14.1–14.35. Wiley, New York (1987) 2. Gebhart, B., Jaluria, Y., Mahajan, R.L., Sammakia, B.: Buoyancy-Induced Flows and Transport. Hemisphere, New York (1988) 3. Schlichting, H., Gersten, K. : Boundary Layer Theory. Springer, New York (2000) 4. Pop, I., Ingham, D.B.: Convective Heat Transfer: Mathematical and Computational Viscous Fluids and Porous Media. Pergamon, Oxford (2001) 5. Bejan, A.: Convective Heat Transfer. 3rd edn. Wiley, New York (2013) 6. Merkin, J.H., Pop, I.: Conjugate free convection on a vertical surface. Int. J. Heat Mass Transf. 39, 1527–1534 (1996) 7. Merkin, J.H.: Natural convection boundary-layer flow on a vertical surface with Newtonian heating. Int. J. Heat Fluid Flow 15, 392–398 (1994) 8. Aziz, A.: A similarity solution for laminar thermal boundary layer over flat plate with convective surface boundary condition. Commun. Nonlinear Sci. Numer. Simul. 14, 1064–1068 (2009) 9. Ishak, A.: Similarity solutions for flow and heat transfer over permeable surface with convective boundary conditions. Appl. Math. Comput. 217, 837–842 (2010) 10. Magyari, E.: Comment on “A similarity solution for laminar thermal boundary layer over a flat plate with a convective surface boundary condition” by A. Aziz, Comm. Nonlinear Sci. Numer. Simul. 2009; 14: 1064-8. Commun. Nonlinear Sci. Numer. Simul. 16, 599–601 (2011) 11. Bataller, R.C.: Similarity solutions for flow and heat transfer of a quiescent fluid over a nonlinearly stretching surface. J. Mater. Process. Technol. 203, 176–183 (2008) 12. Yao, S., Fang, T., Zhong, Y.: Heat transfer of a generalized stretching/shrinking wall problem with convective boundary conditions. Commun. Nonlinear Sci. Numer. Simul. 16, 752–760 (2011) 13. Hansen, A.G.: Similarity Analyses of Boundary Value Problems in Engineering. Prentice-Hall, New Jersey (1964) 14. Merkin, J.H., Mahmood, T.: Mixed convection boundary layer similarity solutions: prescribed wall heat flux. J. Appl. Math. Phys. (ZAMP) 20, 51–68 (1989) 15. Mahmood, T., Merkin, J.H.: Mixed convection on a vertical circular cylinder. J. Appl. Math. Phys. (ZAMP) 39, 186– 203 (1988) 16. Mahmood, T., Merkin, J.H.: Similarity solutions in axisymmetric mixed-convection boundary-layer flow. J. Eng. Math. 22, 73–92 (1988) 17. Merkin, J.H., Pop, I.: Mixed convection along a vertical surface: similarity solutions for uniform flow. Fluid Dyn. Res. 30, 233–250 (2002) 18. Rhida, A.: Aiding flows non-unique similarity solutions of mixed-convection boundary-layer equations. J. Appl. Math. Phys. (ZAMP) 47, 341–352 (1996) 19. Merkin, J.H., Pop, I.: The forced convection flow of a uniform stream over a flat surface with a convective surface boundary condition. Commun. Nonlinear Sci. Numer. Simul. 16, 3602–3609 (2011) 20. Makinde, O.D., Olanrewaju, P.O.: Buoyancy effects on thermal boundary layer over a vertical plate with a convective surface boundary condition. J. Fluids Eng. 132, 044502 (2010)

M. M. Rahman et al.

21. Pantokratoras, A.: Buoyancy effects on thermal boundary layer over a vertical plate with a convective surface boundary condition. New results. Int. J. Therm. Sci. 76, 221–224 (2014) 22. White, F.: Viscous Fluid Flow. 3rd edn. McGraw-Hill, New York (2006) 23. Ro¸sca, A.V., Pop, I.: Flow and heat transfer over a vertical permeable stretching/shrinking sheet with a second order slip. Int. J. Heat Mass Transf. 60, 355–364 (2013) 24. Rahman, M.M., Al-Lawatia, M.A., Eltayeb, I.A., Al-Salti, N.: Hydromagnetic slip flow of water based nanofluids past a wedge with convective surface in the presence of heat generation (or) absorption. Int. J. Therm. Sci. 57, 172–182 (2012) 25. Rahman, M.M., Eltayeb, I.A.: Radiative heat transfer in a hydromagnetic nanofluid past a nonlinear stretching surface with convective boundary condition. Meccanica 48, 601–615 (2013) 26. Rahman, M.M., Ro¸sca, A.V., Pop, I.: Boundary layer flow of a nanofluid past a permeable exponentially shrinking surface with convective boundary condition using Buongiorno’s model. Int. J. Numer. Methods Heat Fluid Flow (2014) 27. Rahman, M.M., Ro¸sca, A.V., Pop, I.: Boundary layer flow of a nanofluid past a permeable exponentially shrinking surface with second order slip using Buongiorno’s model. Int. J. Heat Mass Transf. (2014) 28. Rahman, M.M., Grosan, T., Pop, I.: Oblique stagnation-point flow of a nanofluid past a shrinking sheet (under review) (2014b) 29. Shampine, L.F., Reichelt, M.W., Kierzenka, J.: Solving boundary value problems for ordinary differential equations in Matlab with bvp4c. http://www.mathworks.com/bvp_tutorial (2010) 30. Shampine, L.F., Gladwell, I., Thompson, S.: Solving ODEs with MATLAB. Cambridge University Press, Cambridge (2003) 31. Falkner, V.M., Skan, S.W.: Solutions of the boundary layer equations. Philos. Mag. 12, 132–135 (1931) 32. Stewartson, K.: Further solutions of the Falkner-Skan equation. Proc. Camb. Philos. Soc. 50, 454–465 (1954) 33. Rosenhead, L.: Laminar Boundary Layers. Clarendon Press, Oxford (1963) 34. Merkin, J.H., Mahmood, T.: Mixed convection boundary layer similarity solutions: prescribed wall heat flux. J. Appl. Math. Phys. (ZAMP) 20, 51–68 (1989) 35. Merkin, J.H.: Mixed convection in a Falkner–Skan system (in preparation) 36. Merkin, J.H., Pop, I.: Natural convection boundary-layer flow in a porous medium with temperature dependent boundary conditions. Transp. Porous Med. 85, 397–414 (2010) 37. Steinrück, H.: About the physical relevance of similarity solutions of the boundary-layer flow equations describing mixed convection flow along a vertical plate. Fluid Dyn. Res. 32, 1–13 (2003)