Long-wave instability of flow with temperature dependent fluid properties down a heated incline

International Journal of Engineering Science 70 (2013) 73–90

Contents lists available at SciVerse ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Long-wave instability of flow with temperature dependent fluid properties down a heated incline J.P. Pascal a,⇑, N. Gonputh a, S.J.D. D’Alessio b a b

Department of Mathematics, Ryerson University, Toronto, Ontario, Canada M5B 2K3 Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1

a r t i c l e

i n f o

Article history: Received in revised form 6 April 2013 Available online 5 June 2013 Keywords: Variable fluid properties Marangoni effect Inclined flow Interfacial instability

a b s t r a c t We present an analytical study which investigates the effect of temperature dependence in fluid properties on the interfacial instability of flow down a heated incline. Along with temperature variation in surface tension we consider variable mass density, viscosity, thermal conductivity and specific heat. A linear stability analysis is carried out which yields the critical conditions for the onset of instability in long-wave perturbations. Results are obtained for the particular case when there is variation only in surface tension, density and specific heat, and in the case with negligible and high rate of heat transfer across the free surface. For the general case, asymptotic expansions are implemented based on the assumed smallness of the variation with temperature in viscosity and thermal conductivity, or on weak heat transfer across the free surface. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The flow of a thin fluid layer over an inclined surface continues to be an active area of research. Despite the numerous investigations devoted to this problem, voids in our understanding still remain. This study addresses one such void. Inclined free-surface flow is subject to inertial instability which leads to the formation of surface waves as first observed experimentally by Kapitza and Kapitza (1949). The critical conditions for the onset of instability of the uniform flow can be obtained from a linearization of the perturbation equations followed by asymptotic expansions with respect to small wavenumbers. This analysis was originally carried out by Benjamin (1957) and Yih (1963). These studies were followed by many undertakings focussing on different aspects of the flow which are summarized in the recent reviews by Kalliadasis, Ruyer-Quil, Scheid, and Velarde (2012) and Craster and Matar (2009). Over the years several models have been advanced to simulate these flows. Key models include the integral-boundary-layer model originally proposed by Shkadov (1967) and the shallow-water model of Balmforth and Mandre (2004). Although these models are successful in many respects, they are plagued by the failure to reproduce the critical conditions for the onset of instability as predicted by Orr–Sommerfeld calculations (Benjamin, 1957; Yih, 1963), more recent experiments (Liu, Paul, & Gollub, 1993; Liu, Schneider, & Gollub, 1995) and direct numerical simulations (Ramaswamy, Chippada, & Joo, 1996). This flaw was rectified by Ruyer-Quil and Manneville (2002) using a weighted-residual technique. For non-isothermal flows the Marangoni effect couples with the inertial force to destabilize the flow. A long-wave asymptotic analysis of the linearized perturbation equations provides a Marangoni number dependent critical Reynolds number (D’Alessio, Pascal, Jasmine, & Ogden, 2010; Trevelyan, Scheid, Ruyer-Quil, & Kalliadasis, 2007). Studies of non-isothermal flows usually focus on the variation of surface tension with temperature, but neglect the variation of other fluid properties

⇑ Corresponding author. E-mail addresses: [email protected] (J.P. Pascal), [email protected] (N. Gonputh), [email protected] (S.J.D. D’Alessio). 0020-7225/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijengsci.2013.05.003

74

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

under the assumption that these variations, being small, will not significantly impact the onset of instability. Indeed, Ogden, D’Alessio, and Pascal (2011) have recently demonstrated that if the variations in mass density and viscosity are small, and in particular of the same order of magnitude as the wavenumber of unstable perturbations, then the critical Reynolds number for the onset of instability is unaffected by these variations. Larger variations have been considered by Goussis and Kelly (1985) and later by Hwang and Weng (1988). However, these investigations only consider the temperature variation in viscosity, and furthermore assume a constant temperature at the free surface which eliminates the Marangoni effect. Kabova and Kuznetsov (2002) include variable viscosity and variable surface tension; however, their investigation is limited to only the steady-state problem. In this study we take all fluid properties including surface tension, mass density, viscosity, thermal conductivity and specific heat to be temperature dependent. We perform a linear stability analysis which yields the critical conditions for the onset of instability and investigate the effect of variations in fluid properties as a result of temperature changes. The paper is organized as follows. In the next section we present the mathematical formulation, governing equations and boundary conditions. Then in Section 3 we perform a linear stability analysis and highlight the results along with some discussion in Section 4. We summarize the key points in the concluding section. 2. Governing equations The flow configuration that we consider consists of a two-dimensional laminar flow down an impermeable even surface inclined at an angle b with the horizontal. An (x,z) coordinate system is oriented so that the x-axis points along the incline in the downhill direction with the z-axis pointing in the upward normal direction. The x and z components of the fluid velocity are denoted by u and w, respectively, and the position of the free surface of the fluid layer is given by z = h(x,t). The solid substrate is maintained at a constant temperature Tb which is taken to be greater than that of the ambient gas above the fluid layer, denoted by Ta. Due to the fact that we are investigating a flow with variable fluid properties, we begin by considering the general conservation equations for fluid flow (Spurk & Aksel, 2008). For the conservation of mass and momentum we have respectively

Dq þ qr  v ¼ 0; Dt Dv q ¼ b þ r  s; Dt D where Dt denotes the two-dimensional material derivative, q is the mass density, s is the stress tensor, v = (u,w) and b = (qg sinb,  qg cosb) with g being the acceleration due to gravity. The conservation of energy equation is given by

q

  De @ @T @v i þ sji ¼ K ; Dt @xi @xi @xj

where e is the internal energy per unit mass, T is the temperature, K is the thermal conductivity, (v1,v2) = v, (x1,x2) = (x,z) and the repeated indices signify summation. The Newtonian constitutive relation states that

s ¼ pI þ l1 ðtrEÞI þ 2lE; where p is the pressure, l1 and l are the viscosity coefficients, E is the rate of deformation tensor and I is the identity tensor. Rajagopal (2013) has recently shown that it is in fact not appropriate to make the ‘‘Stokes assumption’’ and relate the viscosity coefficients as l1 ¼  23 l. However, in the present investigation we are dealing with a gravity-driven flow and consequently employ the Boussinesq approximation. As such, we assume the density to be

q ¼ q0  a^ ðT  T a Þ; ^ is a positive parameter and q0 is the reference density at temperature Ta. The density variation is assumed to be where a significant only in the gravitational terms and as a result mass conservation gives that r  v = trE is negligible. The constitutive relation thus reduces to the one for an incompressible fluid given by

s ¼ pI þ 2lE: The conservation equations can then be expressed as

@u @w þ ¼ 0; @x @z      Du @p @ @u @ @u @w þ ; q ¼  þ g q sin b þ 2l l þ Dt @x @x @x @z @z @x

ð1Þ ð2Þ

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

     Dw @p @ @w @ @u @w þ ; ¼   g q cos b þ 2l l þ Dt @z @z @z @x @z @x     Dðc TÞ @ @T @ @T þ ; K K q p ¼ Dt @x @x @z @z

q

75

ð3Þ ð4Þ

where cp is the specific heat and the viscous dissipation terms in the energy equation have been discarded in accordance with the Boussinesq approximation. The other fluid properties are also assumed to vary linearly with temperature as follows

l ¼ l0  ^kðT  T a Þ; cp ¼ cp0 þ ^SðT  T a Þ; ^ ðT  T a Þ; K ¼ K0 þ K

r ¼ r0  cðT  T a Þ; ^ and ^ k; K where r is the surface tension, c; ^ S are positive parameters measuring the rate of change with respect to temperature, and l0, cp0, K0 and r0 are reference values at T = Ta of the viscosity, specific heat, thermal conductivity and surface tension, respectively. We next scale the equations relative to the uniform and steady isothermal flow corresponding to the reference fluid properties. For the length scale we choose the Nusselt thickness, H, resulting from a constant isothermal volume flux, Q, which is given by

 H¼

1=3 3l0 Q : g q0 sin b

The pressure is scaled using q0U2 where U = Q/H is the velocity scale, while the corresponding time scale is taken to be H/U. Lastly, the scaled temperature difference is T⁄ = (T  Ta)/DT, with DT = Tb  Ta. The dimensionless conservation equations, with the asterisk removed from the scaled temperature for notational convenience, then become

@u @w þ ¼ 0; @x @z     Du @p @ @u @ @u @T @u @T @w þ k ¼ Re þ 3ð1  aTÞ þ ð1  kTÞ ð1  kTÞ k ; Re Dt @x @x @x @z @z @x @x @z @x Re

    Dw @p @ @w @ @w @T @u @T @w þ k ¼ Re  3 cot bð1  aTÞ þ ð1  kTÞ ð1  kTÞ k ; Dt @z @x @x @z @z @x @z @z @z

PrRe

      D @ @T @ @T ð1 þ S=DT r ÞT þ ST 2 ¼ þ ; ð1 þ KTÞ ð1 þ KTÞ Dt @x @x @z @z

ð5Þ ð6Þ

ð7Þ

ð8Þ

where Re = q0UH/l0, the Prandtl number is Pr = l0cp0/K0, the relative temperature difference is given by DTr = (Tb  Ta)/Ta ^ DT=K 0 and S ¼ ^ ^ DT=q0 , k ¼ ^ and a ¼ a kDT=l0 , K ¼ K SDT=cp0 are the scaled gradients with respect to temperature of the density, viscosity, thermal conductivity and specific heat, respectively. The parameter Re is defined in terms of the reference fluid properties to which the temperature variations are added, and is the Reynolds number for the corresponding isothermal flow. Our aim is to investigate how the gradients of the fluid properties affect the critical value of this Reynolds number, and thus determine the effect of temperature variation on the stability of the flow. At the surface of the fluid layer we assume the ambient gas to be dynamically passive. Consequently, the force exerted by the flow is balanced by surface tension, and the dynamic conditions at the surface (z = h(x,t)) are

p¼

2ð1  kTÞ ReF



@h @x

2

! @u @w @h @u @h @w ðWe  MaTÞ @ 2 h þ    @x @z @x @z @x @x @x2 F 3=2

ð9Þ

and

    pffiffiffi@T @h @T  @u @w @h @u ¼ ð1  kTÞ G 4 ; MaRe F þ þ @x @x @z @z @x @x @x

ð10Þ

where

F ¼1þ



2 @h ; @x

G¼1



2 @h ; @x

We = r0/(q0U2H) is the Weber number and Ma = c DT/(q0U2H) is the Marangoni number which represents the scaled rate of change of surface tension with respect to temperature.

76

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

The transfer of heat across the free surface z = h(x,t) is assumed to be governed by Newton’s Law of Cooling which in nondimensional form is expressed as

  pffiffiffi @T @h @T ; Bi F T ¼ ð1 þ KTÞ  @z @x @x

ð11Þ

where Bi = agH/K0 is the Biot number with ag denoting the heat transfer coefficient across the liquid–air interface. Evaporation of the fluid is assumed to be negligible, and thus conservation of mass at the free surface leads to the following kinematic condition

w¼

@h @h þu : @t @x

At the interface between the fluid layer and the impermeable bottom the tangential and normal fluid velocity components are zero. This results in the no-slip and impermeability conditions

u ¼ w ¼ 0 at z ¼ 0:

ð12Þ

Finally, we have the bottom temperature condition given by

T ¼ 1 at z ¼ 0:

ð13Þ

3. Stability analysis We investigate the stability of a steady uniform flow in the streamwise direction given by h  1, w  0, u = us(z), p = ps(z) and T = Ts(z). The equilibrium temperature distribution, Ts(z), is the solution to the boundary-value problem given by

  d dT s ¼ 0; ð1 þ KT s Þ dz dz dT s þ BiT s ¼ 0 at z ¼ 1; ð1 þ KT s Þ dz T s ð0Þ ¼ 1:

ð14Þ ð15Þ ð16Þ

The exact solution to (14)–(16) is

T s ðzÞ ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 a  bz  ; K

where

a¼

ð1 þ KÞ2

K2

;

b¼

2Bið1 þ BiÞ

K2

"sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # ðð1 þ KÞ2  1Þ 1þ 1 : ð1 þ BiÞ2

For special cases in the limit as K ? 0 this solution simplifies to: Ts(z)  1 when Bi = 0 and Ts(z) = 1  z as Bi ? 1. The equilibrium velocity, us(z), on the other hand, satisfies the following boundary-value problem

  d dus þ 3ð1  aT s Þ ¼ 0; ð1  kT s Þ dz dz dus ¼ 0 at z ¼ 1; dz us ð0Þ ¼ 0:

ð17Þ ð18Þ ð19Þ

Using the known temperature, Ts(z), the exact solution to (17)–(19) has been found to be

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi! pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi A  k a  bz a pffiffiffi us ðzÞ ¼ a0 ln þ a1 z  z2 þ a2 ð a  bz  aÞ þ a3 ½ða  bzÞ3=2  a3=2 ; k Ak a

where

a0 ¼

2AB 2

þ

6ACð1 þ Ka Þ 2 2

bk b k 2ð1 þ Ka Þ 4aA a3 ¼  2 2; 2 b k 3b k



4aA4 2 5

b k

;

A¼1þ

a1 ¼ k ; K

3Að1 þ Ka Þ 6Cð1 þ Ka Þ 2B 4aA3 2a ðC þ 2aÞ  ; a2 ¼ þ  ; 2 2 bk b2 k4 bk bk b k  a  2a A2 B¼3 1þ ða  bÞ3=2 ; C ¼ 2  a: þ K b k

In the isothermal limit it can be shown that the above expression recovers the familiar parabolic profile. Lastly, the equilibrium pressure, ps(z), can be obtained by solving the initial-value problem

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

dps ¼ 3 cot bð1  aT s Þ; dz ps ð1Þ ¼ 0

Re

77

ð20Þ ð21Þ

and has the solution

ps ðzÞ ¼

3 cot b  a 2a cot b ð1  zÞ þ 1þ ½ða  bÞ3=2  ða  bzÞ3=2 : Re bRe K

In order to determine the stability of the equilibrium flow we impose infinitesimal perturbations and obtain the flow parameter range for which the perturbations are amplified. Denoting the perturbations by tildes, the perturbed flow is given by

~ ðx; z; tÞ; u ¼ us ðzÞ þ u

~ z; tÞ; w ¼ wðx;

~ðx; z; tÞ; p ¼ ps ðzÞ þ p

~ z; tÞ; T ¼ T s ðzÞ þ Tðx;

h ¼ 1 þ gðx; tÞ;

where g is the perturbation in the thickness of the fluid layer. We next introduce the perturbed equilibrium flow in the governing equations and linearize in the perturbations. In order to separate the time dependence we consider solutions of the form

~ gÞ ¼ ðu ^ ~ ; w; ~ p ~; T; ^ ðzÞ; wðzÞ; ^ ^ðzÞ; TðzÞ; ðu g^ ÞeikðxctÞ ; p where k is a real positive quantity representing the wavenumber of the perturbation and c is a complex quantity where the real part is the phase speed of the perturbation while the imaginary part multiplied by k is the temporal growth rate. The zdependent amplitudes of the perturbations are governed by the equations

^ þ iku ^ ¼ 0; Dw

ð22Þ

^ 2 us  kDus DT^  ikkwDT ^ ^ þ wDu ^ s  ¼ ikRep ^ þ k2 ðkT s  1Þu ^ þ D½ð1  kT s ÞDu ^   kTD ^ s  3aT; Re½ikðus  cÞu

ð23Þ

2

^ s  kDT s Dw; ^ ¼ ReDp ^ þ 3a cot b T^  k ð1  kT s Þw ^ þ D½ð1  kT s ÞDw ^  ikkTDu ^ ikReðus  cÞw

ð24Þ

^ ^ s  ¼ k2 ð1 þ KT s ÞT^ þ D2 ½ð1 þ KT s ÞT; PrReð1 þ S=DT r þ 2ST s Þ½ikðus  cÞT^ þ wDT

ð25Þ

where D is a differential operator denoting differentiation with respect to z. The boundary conditions at z = 1 are

2 ^ þ k2 ðWe  MaT s Þg ^; ð1  kT s ÞDw Re 2 ^ þ ikwÞ ^ ¼ ikMaReðT^ þ g ^ D us þ Du ^ DT s Þ; ð1  kT s Þðg ^ þg ^ D½ð1 þ KT s ÞDT s þ BiT s  þ BiT^ ¼ 0; D½ð1 þ KT s ÞT

ð28Þ

^ ¼ ikðus  cÞg ^; w

ð29Þ

^ ¼ g ^ Dps þ p

ð26Þ ð27Þ

while at z = 0 the conditions are

^ ^ ð0Þ ¼ wð0Þ ^ ¼ 0: u ¼ Tð0Þ

ð30Þ

Eqs. (22)–(30) are posed as an eigenvalue problem with c being the eigenvalue possessing characteristic values. Based on the expectation that small wavenumber perturbations are the most unstable, we solve the problem asymptotically as k ? 0 by expanding the perturbations in the following series

^ ¼ u0 ðzÞ þ ku1 ðzÞ þ Oðk2 Þ; u ^ ¼ w0 ðzÞ þ kw1 ðzÞ þ Oðk2 Þ; w ^ ¼ p0 ðzÞ þ kp1 ðzÞ þ Oðk2 Þ; p 2 T^ ¼ T 0 ðzÞ þ kT 1 ðzÞ þ Oðk Þ;

g^ ¼ g0 þ kg1 þ Oðk2 Þ; 2

c ¼ c0 þ kc1 þ Oðk Þ: Since the growth rate of the perturbation is kIðcÞ the bifurcating branch of neutral stability is given by IðcÞ ¼ 0. Integrating (22) and imposing (30) leads to

^ ¼ ik w

Z

z

^ ðfÞdf; u 0

then, condition (29) becomes

78

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

ik

Z

1

^ ðfÞdf þ ikðus ð1Þ  cÞg ^ ¼ 0: u

0

Introducing the above expansions in k into this equation, and employing the fact that the eigenvalue problem can be normalized by taking g0 = 1 and g1 = 0, we obtain

c0 ¼

Z

1

u0 ðfÞdf þ us ð1Þ

0

and

c1 ¼

Z

1

u1 ðfÞdf:

0

pffiffiffiffiffiffiffi Examining Eqs. (22)–(30) it can be seen that the terms in even powers of k are not multiplied by ið¼ 1Þ, and consequently the solution for u0(z) must be real. As a result, the condition for neutral stability can be determined once the solution for u1(z) is obtained, and is given by

Z

1

Iðu1 ðfÞÞdf ¼ 0:

ð31Þ

0

Owing to the complicated steady-state solutions for Ts,us,ws and ps, exact solutions for u1(z) have only been obtained for certain special cases. One such case is obtained by setting Bi = 0, while another consists of letting K = k = 0. An exact solution for u1(z) can also be obtained if we set K = 0 and assume that the rate of heat transfer at the free surface is sufficiently high to render the temperature difference between the fluid surface and ambient medium negligible. In implementing the expansions in k all parameters are assumed to be O(1), with the exception of the Weber number. If We is assumed to be sufficiently large so that k2We = O(1) as k ? 0, then our small-wavenumber theory provides a parabolic approximation for the neutral stability curve in the (Re,k) plane. Case I: Bi = 0 We first present the case Bi = 0. Here, the steady-state solutions simplify to

T s ðzÞ  1;

us ðzÞ ¼

  3ð1  aÞ z2 ; z ð1  kÞ 2

ps ðzÞ ¼

3ð1  aÞ cot b ð1  zÞ: Re

Introducing the asymptotic expansions with respect to k into the system (22)–(30) and collecting terms having the same power of k leads to a hierarchy of problems at various orders. The leading-order problem governing T0, w0, u0 and p0 is given by

D2 T 0 ¼ 0;

DT 0 ð1Þ ¼ 0;

Dw0 ¼ 0;

w0 ð0Þ ¼ 0;

D2 u0 ¼ 0; Dp0 ¼ 0;

T 0 ð0Þ ¼ 0;

ð1  aÞ ; u0 ð0Þ ¼ 0; ð1  kÞ 3ð1  aÞ cot b p0 ð1Þ ¼ þ W; Re Du0 ð1Þ ¼ 3

where W = k2We and is assumed to be O(1). The respective solutions are easily found to be

T 0 ðzÞ  0;

w0 ðzÞ  0;

u0 ðzÞ ¼ 3

ð1  aÞ z; ð1  kÞ

From the O(k) problem it follows that

Dw1 ¼ iu0 ;

w1 ð0Þ ¼ 0;

and thus

w1 ¼ i

3ð1  aÞ 2 z : 2ð1  kÞ

From condition (29) we then obtain that

c0 ¼

3ð1  aÞ : ð1  kÞ

To order k, the temperature, T1, satisfies

D2 T 1 ¼ 0; DT 1 ð1Þ ¼ 0; T 1 ð0Þ ¼ 0;

p0 ðzÞ ¼

3ð1  aÞ cot b þ W: Re

79

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

and has the solution T1(z)  0, while for u1 we obtain

ð1  kÞD2 u1 ¼ iRe½p0  ðc0  us Þu0  þ ReDus w1 ;

Du1 ð1Þ ¼ 0;

u1 ð0Þ ¼ 0;

having the solution

u1 ¼

  i ð1  aÞ2 z4 z3 z : ½3ð1  aÞ cot b þ WReðz2  2zÞ þ 9iRe  þ 2ð1  kÞ ð1  kÞ3 24 6 3

Substituting the expression for u1 into (31) yields the following the neutral stability condition

1 6ð1  aÞ2 2 ¼ 0: ½3ð1  aÞ cot b þ k ReWe  Re 3ð1  kÞ 5ð1  kÞ3 Consequently, the critical value of Re for the onset of instability becomes

Recrit ¼

5 ð1  kÞ2 cot b : 6 ð1  aÞ

ð32Þ

Case II: K = k = 0 Next, we consider the case K = k = 0. Here, the steady-state solutions simplify to

Ts ¼ 1 

Bi z; 1 þ Bi 2

us ¼

ð3z þ 3Biz þ aBiz  3az  3azBi  6  6Bi þ 3aBi þ 6aÞz ; 2ð1 þ BiÞ

ps ¼

3 cot bðz  1ÞðazBi þ 2 þ 2Bi  2a  aBiÞ : 2Reð1 þ BiÞ

From the leading-order problem we obtain 2

T0 ¼

Bi

ð1 þ BiÞ2

z;

w0 ðzÞ  0;

2

u0 ¼

2

2

ðaBi z2 þ 6 þ 12Bi þ 6Bi  6a  6aBi  3aBi Þz 2

2ð1 þ BiÞ 2

p0 ¼

2

2

2

3 cot baBi z2 þ 3 cot bð2 þ 4Bi þ 2Bi  2a  2aBi  aBi Þ þ 2WReð1 þ 2Bi þ Bi Þ 2Reð1 þ BiÞ 2

c0 ¼

;

2

;

2

24 þ 48Bi þ 24Bi  24a  9aBi  28aBi 8ð1 þ BiÞ2

:

From the O(k) terms we obtain 2

w1 ¼

2

2

2

iz ða Bi z2 þ 12 þ 24 Bi þ 12 Bi  12 a  12 a Bi  6 a Bi Þ 2

8ð1 þ BiÞ

;

while the problems for T1 and u1 are given respectively by

D2 T 1 ¼ Pr Re½DT s w1  iðc0  us ÞT 0 ð1 þ S=DT r þ 2ST s Þ;

DT 1 ð1Þ þ BiT 1 ð1Þ ¼ 0;

T 1 ð0Þ ¼ 0

and

D2 u1 ¼ iRep0  iReðc0  us Þu0 þ ReDus w1 þ 3aT 1 ;

Du1 ð1Þ þ iMaReðT 0 ð1Þ þ DT s ð1ÞÞ ¼ 0;

u1 ð0Þ ¼ 0:

Although we have obtained the solution for u1 in closed form, the expression is too long and it serves little purpose to list it here. However, using this expression in the calculation of the neutral state yields the following formula for Recrit

80

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90 2

2

Recrit ¼ 1008 cot bð11 Bi a  40 Bi þ 40 Bi a  80 Bi þ 40 a  40Þð1 þ BiÞ4 =ð48384 þ 290304 Bi  96768 a 4

6

5

3

2

 718506 Bi a  37863 Bi a  250092 Bi a þ 7548 Bi Pr Sa þ 40968 Bi Pr Sa þ 22320 Bi Pr Sa 2

3

3

4

þ 11160 Bi Pr S1 a þ 29205 Bi Pr S1 a þ 19080 Bi Pr S1 a þ 1740 Bi Pr Sa2  20124 Pr Bi Sa 5

4

5

6

2

 10620 Pr Bi Sa þ 9195 a2 Pr Bi S þ 3738 a2 Pr Bi S þ 486 a2 Pr Bi S  18648 Bi Pr Sa2 6

3

5

4

 1596 Pr Bi Sa  22320 Bi Pr Sa2  5715 Bi Pr S1 a2 þ 3075 a2 Pr Bi S1 þ 3780a2 Pr Bi S1 6

6

5

4

2

þ 465 a2 Pr Bi S1  1575 Pr Bi a S1  9000 Pr Bi a S1  6390 Pr Bi a S1  18045 Bi Pr S1 a2 3

3

2

 11160 Bi Pr S1 a2 þ 120960 Bi Ma þ 20160 Bi Ma  1130364 Bi a  485712 Bi a  1013031 Bi a 2

3

2

6

4

5

þ 80640 Bi Ma þ 967680 Bi þ 725760 Bi þ 48384 Bi þ 725760 Bi þ 290304 Bi þ 48384 a2 4

5

4

5

6

3

þ 80640 Ma Bi þ 20160 Ma Bi þ 172890 a2 Bi þ 53313 a2 Bi þ 7285 a2 Bi þ 313621 Bi a2 2

þ 195408 Bi a2 þ 333783 Bi a2 Þ;

ð33Þ

where S1 = 1 + S/DTr. Setting a = 0 this formula reduces to

Recrit ¼

10ð1 þ BiÞ2 cot b 5Ma Bi þ 12ð1 þ BiÞ2

which corresponds to the non-isothermal result obtained by D’Alessio et al. (2010) and coincides with that obtained by Trevelyan et al. (2007) if the difference in scaling is taken into account. Case III: Bi ? 1 As a final case which admits an exact solution for u1(z) we consider that of an infinitely high rate of heat transfer at the free surface. This is equivalent to prescribing the temperature of the liquid surface to be that of the ambient medium, i.e. T(x,h,t) = 0, and corresponds to taking Bi ? 1 in the energy condition (11). Consequently, the perturbation condition (28) simplifies to T^ ¼ DT s at z = 1. If we also set K = 0 then the steady state solutions reduce to

T s ðzÞ ¼ 1  z us ðzÞ ¼

3 4k3

½a z2 k2  4 zk2 þ 2

a zk2 þ 2 a zk þ 4 lnðzk þ 1  kÞk  2 lnðzk þ 1  kÞa  4 lnð1  kÞk þ 2 lnð1  kÞa

and

ps ðzÞ ¼

3 cot b ð1  zÞðaz þ 2  aÞ : 2Re

The solutions to the leading-order problem are given by

T 0 ¼ z; u0 ¼

w0 ðzÞ  0; 3 3

2k ð1  k þ k zÞ

½4 lnð1  k þ k zÞk2 þ 4 lnð1  k þ k zÞk  2 lnð1  k þ k zÞa  2 lnð1  k þ k zÞa zk

þ 2 lnð1  kÞa k z þ 2 lnð1  kÞa  4 lnð1  kÞk þ 4 lnð1  kÞk2 þ az2 k2 þ 2azk  4zk2  azk3 þ a z2 k3  2 lnð1  kÞak þ 4 lnð1  k þ kzÞzk2 þ 2 lnð1  k þ k zÞak  4 lnð1  kÞk2 z þ 2 zk3  2z2 k3 ;

p0 ¼

3 cot baz2 þ 2WRe þ 6 cot b  3 cot ba ; 2Re

c0 ¼ 3

6 lnð1  kÞa  12 lnð1  kÞk þ 3ak2 þ 6ak  12 k2 þ 2ak3  6k3

The O(k) problem is governed by

4k4

:

81

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

w1 ¼

3 i

½12 lnð1  k þ kzÞk2 þ 12 lnð1  k þ k zÞk  6 lnð1  k þ k zÞa  4 lnð1  k þ k zÞa zk þ 4 lnð1  kÞa k z 4k4 þ 6 lnð1  kÞa  12 lnð1  kÞk þ 12 lnð1  kÞk2 þ az2 k2 þ 6azk  12 zk2 þ a z2 k3  6 lnð1  kÞa k þ 8 lnð1  k þ k zÞzk2 þ 6 lnð1  k þ k zÞa k  8 lnð1  kÞk2 z  2 z2 k3 ;

D2 T 1 ¼ PrReð1 þ S=DT r þ 2ST s Þ½DT s w1  iðc0  us ÞT 0 ; T 1 ð1Þ ¼ 0; T 1 ð0Þ ¼ 0; D½ð1  kT s ÞDu1  ¼ kT 1 D2 us þ kDus DT 1 þ 3aT 1 þ iRe½p0  ðc0  us Þu0  þ Re w1 Dus ; ð1  kT s ÞDu1 ¼ iMaReðT 0 þ DT s Þ at z ¼ 1;

u1 ð0Þ ¼ 0:

We solved these problems exactly and obtained a formula for the critical Re value. Since the expression is too long to display, we present reduced formulations obtained by setting certain parameters to zero. As such, if we let a = 0 then the expression reduces to



Recrit ¼ 1200k6 cot b 2 lnð1  kÞ þ 2 k þ k2 =½9000k3  45000 lnð1  kÞk2  26700Pr S1 k3  29040 Pr Sk2 þ 12480Pr Sk3 þ 22650Pr S1 k4  38700 k4  12300k5  2700 k6  28800Pr S1 ðlnð1  kÞÞ2 k þ 43200Pr Sðlnð1  kÞÞ2 k þ 50400Pr S1 ðlnð1  kÞÞ2 k2  61440 Pr S lnð1  kÞk þ 72000Pr S1 lnð1  kÞk3 2

 55500Pr S1 lnð1  kÞk2 þ 54000Pr S lnð1  kÞk2 þ 3600 lnð1  kÞk4 þ 48600 ðlnð1  kÞÞ k3  140400 lnð1  kÞk3  124200ðlnð1  kÞÞ2 k2  54000 ðlnð1  kÞÞ2 k þ 25200ðlnð1  kÞÞ3 k2  54000ðlnð1  kÞÞ3 k þ 100Pr S1 k5  180Pr S1 k7  180Pr Sk7 þ 11020Pr Sk4  1458Pr Sk6  660Pr Sk5  1575 Pr S1 k6  21600Pr S1 ðlnð1  kÞÞ2 k3  10800Pr Sðlnð1  kÞÞ2 k3 þ 18000Pr S lnð1  kÞk3  2880Pr S lnð1  kÞk5  3600Pr S1 lnð1  kÞk5  15000Pr S1 lnð1  kÞk4  10500Pr S lnð1  kÞk4  32400Pr Sðlnð1  kÞÞ2 : On the other hand, with k = 0 the expression for Recrit simplifies to

1008ð40  11aÞ cot b : 48384  37863a þ 7285a2 þ 486Pr Sa2 þ 465Pr S1 a2  1575Pr S1 a  1596PrSa

Recrit ¼

We note that since we have k = 0 and K = 0 this result corresponds to that given in (33) with the additional assumption of a high rate of heat transfer at the free surface. Indeed, the formula above is identical to the expression for Recrit given by (33) in the limit as Bi ? 1. General case For the general case we can obtain an analytical expression for the neutral stability state by implementing asymptotic expansions as Bi ? 0 or as K, k ? 0. To accomplish this we must relate the magnitude of the expansion parameters to that of k, as k ? 0. We first consider the case where Bi ? 0 and assume that Bi approaches 0 like k1/n, where n is a positive integer. ^  and Then, to determine which terms in the expansion are negligible given that O(k2) quantities are negligible, we let Bi ¼ Bi ^n where Bi ^ are of order unity. We must thus include all powers of Bi and k that introduce powers of  up to 2n. So, ^ and k k¼k if f(z) denotes an arbitrary function, the expansion for f(z) has the form

f ðzÞ ¼

2n1 X

^ ^ j j þ k f 0j ðzÞBi

j¼0

n1 X

^ j jþn þ Oð2n Þ; f 1j ðzÞBi

j¼0

Clearly, we expect the accuracy of the approximation to improve with n. However, the calculations become very tedious for large n due to the increase in the hierarchy of problems that must be solved. We were able to perform the calculations up to n = 4. In this case, the expansions in terms of Bi and k for the steady-state solutions are

T s ðzÞ ¼

7 X j 8 T sj ðzÞBi þ OðBi Þ; j¼0

7 X j 8 usj ðzÞBi þ OðBi Þ; us ðzÞ ¼ j¼0 7 X j 8 psj ðzÞBi þ OðBi Þ; ps ðzÞ ¼ j¼0

82

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

while for the perturbation amplitudes we have

^ TðzÞ ¼

7 3 X X 2 j j 8 T 0j ðzÞBi þ k T 1j ðzÞBi þ Oðk ; Bi Þ; j¼0

j¼0

7 3 X X 2 j j 8 ^ ðzÞ ¼ u u0j ðzÞBi þ k u1j ðzÞBi þ Oðk ; Bi Þ; j¼0

j¼0

7 3 X X 2 j j 8 ^ wðzÞ ¼ w0j ðzÞBi þ k w1j ðzÞBi þ Oðk ; Bi Þ; j¼0 7 X

^ðzÞ ¼ p

j¼0

3 X 2 j j 8 p0j ðzÞBi þ k p1j ðzÞBi þ Oðk ; Bi Þ:

j¼0

j¼0

Similarly, the eigenvalue is expanded as

c¼

7 3 X X 2 j j 8 c0j Bi þ k c1j Bi þ Oðk ; Bi Þ: j¼0

j¼0

We normalize the eigenvalue problem by setting g = 1. Lastly, we point out that we already have the expressions for the leading-order terms in Bi for these expansions since they correspond to the solutions for the case corresponding to Bi = 0. ^ ðzÞ and is given by As explained earlier, the neutral stability condition only involves the k terms in the expansion of u

Z

1

h

i 2 3 Iðu10 ðzÞÞ þ Iðu11 ðzÞÞBi þ Iðu12 ðzÞÞBi þ Iðu13 ðzÞÞBi dz ¼ 0:

0

Introducing the asymptotic expansions into system (22)–(30) and collecting like terms in k and Bi again leads to a hierarchy of problems at various orders. We emphasize that since the formula for the neutral stability only contains powers of Bi up to power three, we do not need to consider any of the problems arising from higher powers of Bi, including the leading-order problem in k. We obtained closed-form solutions to all the required problems with the aid of the Maple Computer Algebra System; however, the formulation of the problems and the solutions, including the expression for Recrit are extremely long and it is pointless to present them here. Nevertheless, we were able to utilize the formula for Recrit to generate plots. Linearized with respect to the parameters a, k and K this formula becomes

Recrit ¼

1 cot b 3  ð80640  161280 k þ 80640 a  384100 Bi Pr Sa 672 ð12  10 Bi2 Ma þ 5 Bi Ma þ 15 Bi3 MaÞ2 3

2

2

þ 166360 Bi Pr Sk þ 16320 Bi k SPr þ 154920 Bi Pr Sa  68160 Bi Pr Sk  37200 Bi Pr Sa 2

2

3

3

 18600 Bi Pr S1 a þ 62925 Bi Pr S1 a  27870 Bi Pr S1 k þ 56820 Bi Pr S1 k  130350 Bi Pr S1 a 3

3

2

þ 8160 Bi k S1 Pr þ 100800 Bi Ma þ 33600 Bi Ma  153138 Bi a  77520 Bi a þ 60720 Bi k þ 115329 Bi a 2

2

3

6

2

2

 67200 Bi Ma  84690 Bi k þ 108660 Bi k  25200 Bi k Ma þ 100800 Bi a Ma þ 2800 Bi Ma k 3

4

2

 8400 Bi Ma k þ 238560 Bi a Ma þ 33600 Bi Ma K  33600 a Bi Ma  168000 Bi Ma K 6

5

3

5

3

þ 156240 Bi a Ma  232680 Bi a Ma þ 453600 Bi Ma K þ 42000 Bi k Ma  210840 Bi a Ma 4

 50400 Bi k MaÞ þ Oða2 ; k2 ; K2 Þ: In order to verify our full expression for Recrit we first let Bi = 0 which, as expected, reproduces the formula given by (32) pertaining to the case Bi = 0. Next, setting k = K = 0 we reproduce the expression given by (33), apart from an expected small correction which is O(Bi4) as Bi ? 0. The accuracy of our small Bi approximation can be demonstrated by comparing its variation with Bi against those for which the variation is known. This is illustrated in Fig. 1 for the case K = k = 0. ^n where ^ ^ are O(1), then ^  and k ¼ k ^ and k Next, we take Bi to be O(1), and let k and K approach 0. If we set k ¼ ^ k, K ¼ K k; K an arbitrary function f(z) is expanded as

f ðzÞ ¼

2X n12n1j X j¼0

n1 n1j X X ^ ^ m ^kj jþm þ k ^ m ^kj jþmþn þ Oð2n Þ: f0jm ðzÞK f1jm ðzÞK

m¼0

j¼0 m¼0

With n = 4 the expansion for f(z) in terms of k, K and k then becomes

f ðzÞ ¼

7j 3j 7 X 3 X X X 2 f0jm ðzÞKm kj þ k f1jm ðzÞKm kj þ Oðk ; k8 ; K8 Þ: j¼0 m¼0

j¼0 m¼0

83

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

2.2

2

Recrit/cot

1.8

1.6

1.4

Approx. Bi

0

1.2 formula (33) 1

0

0.05

0.1

0.15

0.2

0.25 Bi

0.3

0.35

0.4

0.45

0.5

Fig. 1. Accuracy of the expansion as Bi ? 0 for the case with k = K = 0 for a = 0.6, S = 0, D Tr = 1, Ma = 1 and Pr = 7.

We implement such an expansion for the steady-state solutions and the eigenvalue problem involving the amplitudes of the perturbations. A hierarchy of problems can be constructed by collecting like powers of the expansion parameters. Neutral stability is found by imposing the condition

Z 0

1

3j 3 X X

Iðu1jm ðzÞÞKm kj dz ¼ 0:

j¼0 m¼0

In order to determine the essential terms in the expansion for u(z) to complete the above equation we must solve the problems associated with powers of k and K up to degree three. An analytical expression for Recrit was obtained from this analysis, but again the formulation is too long to present in its entirety. The linearized version of the formula with respect to a, k and K can be expressed as

Recrit ¼

1  672

cot b

3

2

2  ð80640 þ 33600 K Bi Ma þ 322560 Bi þ 6120 Pr Bi a S 12 Bi þ 5 Ma Bi þ 24 Bi þ 12 2

2

4

3

2

4

 809490 k Bi  124530 kBi þ 2800 k Ma Bi þ 2800kMa Bi þ 16320k Bi Pr S  2600 k Pr Bi S 3

2

3

4

4

 8360kPr Bi S  2880kPr Bi S þ 12380Pr Bi Sa þ 2660Pr Bi Sa  37200Pr Bi Sa  2310kPr Bi S1 2

3

3

2

4

þ 4770k Pr Bi S1 þ 9750Pr Bi aS1  5700kPr Bi S1  11475Pr Bi aS1 þ 2625Pr Bi aS1  18600 Pr BiS1 a 3

2

2

þ 8160 k BiPr S1 þ 33600 MaBi  16800K Bi Ma þ 483840 Bi þ 67200 Ma Bi þ 322560 Bi þ 165618 a Bi 2

2

3

2

þ 289089 a Bi  161280 k  33600 Bi Ma a þ 80640 a  33600 K Bi Ma þ 245040 a Bi  33600Bi Ma a 3

4

4

3

3

 9240 Bi Ma a þ 80640 Bi þ 40929Bi a  584400 k Bi  510900 k Bi þ 33600 Ma Bi Þ þ Oða2 ; k2 ; K2 Þ: If we set k = K = 0, the formula for Recrit obtained from the asymptotic analysis as k, K ? 0 reduces to the expression given by (33). When Bi = 0 our approximation for Recrit reduces to

Recrit ¼



80 k þ k2 þ 1 þ k3 cot b 3ð96 k þ 192 k2 þ 263 k3 þ 32  96ka  192 k2 a  263 k3 a  32aÞ

;

84

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

2

1.8

1.6

Recrit/cot

1.4

1.2 Approx. as

0

1

0.8 formula (32)

0.6

0.4

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Fig. 2. Accuracy of the expansion as k ? 0 for the case with Bi = 0 for a = 0.5.

1 approximate solution, Bi=0.5 exact solution, Bi=0.5 approximate solution, Bi=0.75 exact solution, Bi=0.75

0.95

steady state temperature

0.9

0.85

0.8

0.75

0.7

0.65

0

0.1

0.2

0.3

0.4

0.5 z

0.6

0.7

0.8

Fig. 3. Steady-state temperature as a function of z with K = 0.5.

0.9

1

85

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

1.8 approximate solution, = =0.5 exact solution, = =0.5 approximate solution, = =0.3 exact solution, = =0.3

1.6

1.4

steady state velocity

1.2

1

0.8

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5 z

0.6

0.7

0.8

0.9

1

Fig. 4. Steady-state velocity as a function of z with a = 0.5 and Bi = 1.

1.4 =0.5 =1 =1.5

1.2

1

Recrit/cot

0.8

0.6

0.4

0.2

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 5. Scaled critical Reynolds number as a function of k with Bi = 0.25, a = 0.5, S = 1, Ma = 1, DTr = 1 and Pr = 7.

86

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

1.2 Tr=0.1 Tr=0.5 Tr=1

1

Tr=2

Recrit/cot

0.8

0.6

0.4

0.2

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 6. Scaled critical Reynolds number as a function of k with Bi = 0.25, a = 0.5, S = 1, Ma = 1, K = 1 and Pr = 7.

which agrees with the expression given by (32) to within O(k3) as k ? 0. In Fig. 2 we illustrate the accuracy of the asymptotic approximation as k ? 0 for the case Bi = 0 by comparing the variation with k with that of formula (32). As final illustrations of the accuracy of the approximate solutions, contrasted in Figs. 3 and 4 are the asymptotic steadystate temperature and velocity profiles with their exact solutions for various values of the parameters, respectively. 4. Results and discussion As previously specified, Re corresponds to a Reynolds number based on the reference values of the fluid properties. Thus, in order to determine the impact of the variation of the fluid properties with temperature on the stability of the flow, we investigate the dependence of the critical value of this Reynolds number for the onset of instability on the parameters measuring the rate of change with temperature for the various fluid properties. Our assumption is that the mass density and the viscosity decrease with temperature in accordance with the scaled formulae 1  aT and 1  kT, respectively. Since the scaled temperature difference, T, attains a maximum value of 1, in order for these quantities to be positive the parameters a and k must be restricted to positive values less than 1. The results displayed in Figs. 5–9 illustrate the dependence of Recrit on the parameters a, k, K, S and DTr as determined by the asymptotic analysis as Bi ? 0 which allows these parameters to be O(1). Fig. 10 is concerned with identifying the effect of varying Bi over a wide range and was determined by the asymptotic analysis as k, K ? 0 with Bi assumed to be O(1). We begin by noting an interesting feature associated with the expression for Recrit given by (32) for the case Bi = 0. If the values of a and k are such that

1  a ¼ ð1  kÞ2 ; then it follows from (32) that the effects of variable fluid properties cancel and the threshold of instability is the same as that for isothermal flow. In Figs. 5 and 6 we focus on the variation of Recrit with k. Increasing k results in a decrease in viscosity which in turn increases the flow inertia and hence destabilizes the flow. Accordingly, the critical Reynolds number is a decreasing function of k in all cases. Fig. 5 contains the Recrit distribution with k for different values of K, the scaled gradient of thermal conductivity with temperature. For values of k up to approximately 0.5 it is apparent that Recrit increases with K. Therefore, the indication is that in the presence of strong thermal conductivity the heating induced by the substrate stabilizes the flow since it reduces temperature differences and thus weakens thermocapillary effects. For the larger values of k the viscosity is greatly reduced

87

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

1.3 S=0 S=2.5 S=5

1.2

1.1

Recrit / cot

1

0.9

0.8

0.7

0.6

0.5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 7. Scaled critical Reynolds number as a function of a with Bi = 0.25, k = 0.2, K = 1, Ma = 1, D Tr = 1 and Pr = 7.

and the accompanying increase in inertia dominates thermocapillarity to the point where flow instability is unaffected by changes in thermal conductivity. By examining the curves displayed in Fig. 6 we can draw conclusions regarding the effect of the parameter DTr which measures the ratio of the temperature difference between the substrate and the ambient medium to the temperature of the ambient medium. It can be seen that for the smaller values of DTr there is less variation in the critical Reynolds number with k. Small values of D Tr corresponds to small temperature differences which results in a reduced dependence. It is also apparent that the critical Reynolds number quickly approaches a limiting value with increasing D Tr. In Figs. 7 and 8 we illustrate the dependence of flow instability on the density variation parameter a. Generally, a decrease in density reduces flow inertia and has a stabilizing effect. However, because of bottom heating a top-heavy density stratification complicates the situation. The depth fluctuations associated with surface undulations can be affected by the density stratification and can result in the amplification of the undulations. The results in Fig. 7 indicate that for sufficiently large values of the specific heat variation parameter S, increasing a destabilizes the flow. This phenomenon is due to the fact that an increase in specific heat reduces thermal diffusivity thereby steepening temperature gradients which accentuates the density stratification. In Fig. 8 we see that increasing Ma has the same effect. More specifically, sufficiently strong thermocapillarity, as measured by the Marangoni number, can reverse the stabilizing role played by a temperature dependent density. Fig. 9 illustrates the graph of Recrit as a function of Ma for various values of k. Since thermocapillarity is a destabilizing agent, all the Recrit distributions in Fig. 9 decrease with Ma. It is also evident that increasing k results in a diminished rate of change in Recrit with Ma, and ultimately with k = 0.75 the critical Reynolds number is almost independent of Ma. Based on these results we can conclude that if the decrease in viscosity is sufficiently high, then flow instability is dominated by inertia and the thermocapillary effect is negligible. The dependence of the critical Reynolds number on the Biot number, Bi, is illustrated in Fig. 10. Since Bi is a measure of the heat transfer coefficient at the free surface, a value of zero corresponds to the surface being insulated and the steady temperature distribution being uniform. In this case surface tension is constant and the Marangoni phenomenon does not manifest itself. As Bi is increased thermocapillarity becomes more significant and Recrit decreases. In the limit as Bi approaches infinity, however, we expect the surface temperature to approach that of the ambient medium because of Newton’s Law of Cooling and we would again have a uniform surface temperature which neutralizes the Marangoni effect. So, for sufficiently large Bi the critical Reynolds number must increase as the thermocapillary contribution to instability diminishes. The curves plotted in Fig. 10 clearly demonstrate the described behavior. The minimum value of Recrit and the critical value

88

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

1.2 Ma=1 Ma=2 Ma=5 Ma=10

1.1

1

Recrit / cot

0.9

0.8

0.7

0.6

0.5

0.4

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 8. Scaled critical Reynolds number as a function of a with Bi = 0.25, k = 0, S = 1, K = 0.5, D Tr = 1 and Pr = 7.

1.1 =0 =0.25 =0.5 =0.75

1 0.9 0.8

Recrit / cot

0.7 0.6 0.5 0.4 0.3 0.2 0.1

0

1

2

3

4

5 Ma

6

7

8

9

10

Fig. 9. Scaled critical Reynolds number as a function of Ma with Bi = 0.25, a = 0.5, S = 1, K = 0.5, DTr = 1 and Pr = 7.

89

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

1.1 =0 =0.25

1

Recrit / cot

0.9

0.8

0.7

0.6

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Bi Fig. 10. Scaled critical Reynolds number as a function of Bi with a = 0.5 with k = 0.2, S = 1, DTr = 1, Ma = 1 and Pr = 7.

of Bi at which it is attained has a complicated dependence on the various parameters in this problem. The cases presented in Fig. 10, for example, suggest that the critical value of Bi increases with K.

5. Conclusions Presented here was an analytical investigation into the impact caused by having temperature-dependent fluid properties on the stability of gravity-driven flow down a heated incline. Although flow down an incline has received much attention, this appears to be the first study which considers variations in all the fluid properties. Through a combination of special cases and asymptotic expansions focussing on long-wave perturbations, our linear stability analysis uncovered the following key results: (i) If Bi = 0, then the critical value for the onset of instability of the Reynolds number based on the reference fluid properties is given by

Recrit ¼

5 ð1  kÞ2 cot b : 6 ð1  aÞ

(ii) If K = k = 0 an exact lengthy expression for Recrit was obtained which reproduces known results when appropriate limits are taken. (iii) If Bi ? 1 another lengthy expression for Recrit was obtained which can confirm the result in (ii) in the limit of large Bi. (iv) Increasing k decreases viscosity which destabilizes the flow since it enhances flow inertia. (v) Increasing a decreases density which, in general, stabilizes the flow since it reduces flow inertia. Exceptions to this occur in extreme cases of large variations in the specific heat and for large Marangoni numbers. (vi) The behavior of Recrit with the heat transfer coefficient is more complicated. As Bi increases from zero Recrit initially decreases and reaches a minimum value; as Bi is increases further Recrit increases. The limiting values of Recrit as Bi ? 0 and Bi ? 1 are given by (i) and (iii), respectively.

90

J.P. Pascal et al. / International Journal of Engineering Science 70 (2013) 73–90

References Kapitza, P. L., & Kapitza, S. P. (1949). Wave flow of thin layers of a viscous fluid. Part III: Experimental study of undulatory flow conditions. Journal of Experimental and Theoretical Physics, 19, 105–120. Benjamin, T. B. (1957). Wave formation in laminar flow down an inclined plane. Journal of Fluid Mechanics, 2, 554–574. Yih, C.-S. (1963). Stability of liquid flow down an inclined plane. Physics of Fluids, 6, 321–334. Kalliadasis, S., Ruyer-Quil, C., Scheid, B., & Velarde, M. G. (2012). Falling liquid films. London, UK: Springer-Verlag. Craster, R. V., & Matar, O. K. (2009). Dynamics and stability of thin liquid films. Reviews of Modern Physics, 81, 1131–1198. Shkadov, V. Y. (1967). Wave conditions in flow of thin layer of a viscous liquid under the action of gravity. Izvestiya Akademii Nauk SSSR, Mekhanika Zhidkosti i Gaza, 1, 43–50. Balmforth, N., & Mandre, S. (2004). Dynamics of roll waves. Journal of Fluid Mechanics, 514, 1–33. Liu, J., Paul, J. D., & Gollub, J. P. (1993). Measurements of the primary instabilities of film flows. Journal of Fluid Mechanics, 250, 69–101. Liu, J., Schneider, B., & Gollub, J. P. (1995). Three-dimensional instabilities of film flows. Physics of Fluids, 7, 55–67. Ramaswamy, B., Chippada, S., & Joo, S. W. (1996). A full-scale numerical study of interfacial instabilities in thin-film flows. Journal of Fluid Mechanics, 325, 163–194. Ruyer-Quil, C., & Manneville, P. (2002). Further accuracy and convergence results on the modeling of flows down inclined planes by weighted-residual approximations. Physics of Fluids, 14, 170–183. Trevelyan, P. M. J., Scheid, B., Ruyer-Quil, C., & Kalliadasis, S. (2007). Heated falling films. Journal of Fluid Mechanics, 592, 295–334. D’Alessio, S. J. D., Pascal, J. P., Jasmine, H. A., & Ogden, K. A. (2010). Film flow over heated wavy inclined surfaces. Journal of Fluid Mechanics, 665, 418–456. Ogden, K. A., D’Alessio, S. J. D., & Pascal, J. P. (2011). Gravity-driven flow over heated, porous, wavy surfaces. Physics of Fluids, 23, 122102. Goussis, D. A., & Kelly, R. E. (1985). Effects of viscosity on the stability of film flow down heated or cooled inclined surfaces: long-wavelength analysis. Physics of Fluids, 28, 3207–3214. Hwang, C.-C., & Weng, C.-I. (1988). Non-linear stability analysis of film flow down a heated or cooled inclined plane with viscosity variation. International Journal of Heat Mass Transfer, 31, 1775–1784. Kabova, Y. O., & Kuznetsov, V. V. (2002). Downward flow of a nonisothermal thin liquid film with variable viscosity. Journal of Applied Mechanics and Technical Physics, 43, 895–901. Spurk, J. H., & Aksel, N. (2008). Fluid mechanics (second ed.). Berlin, Germany: Springer. Rajagopal, K. R. (2013). A new development and interpretation of the Navier–Stokes fluid which reveals why the ‘‘Stokes assumption’’ is inapt. International Journal of Nonlinear Mechanics, 50, 141–151.