Continuum Mech. Thermodyn. (2005) 17: 17–27 Digital Object Identifier (DOI) 10.1007/s00161-004-0184-2

Nonlinear convection of a viscoelastic fluid with inclined temperature gradient P.N. Kaloni, J.X. Lou Department of Mathematics and Statistics, University of Windsor, Windsor, Ontario, Canada N9B 3P4 Received March 12, 2004 / Accepted April 19, 2004 Published online September 17, 2004 – © Springer-Verlag 2004 Communicated by B. Straughan

Abstract. The energy method is used to analyze the viscoelastic fluid convection problem in a thin horizontal layer, subjected to an applied inclined temperature gradient. The boundaries are considered to be rigid and perfectly conducting. Both linear and nonlinear stability analyses are carried out. The eigenvalue problem is solved by the Chebyshev Tau-QZ method and comparisons are reported between the results of the linear theory and energy stability theory. Key words: viscoelastic fluid, energy method, inclined temperature gradient, Chebyshev Tau-QZ method PACS: 47.20 Ky, 47-27 Te, 83.60 Wc 1 Introduction To study the hydrodynamic stability problems, there are two prevailing theories: the linearized stability theory and the energy stability theory. The linearized theory provides sufficient condition for the disturbance of the basic flow to be unstable, while energy stability theory provides sufficient condition for the disturbance to be asymptotically stable (in the Lyapunov sense). In cases, however, when the operator representing the linear system of governing equations happens to be self adjoint, the stability predictions, based upon the two theories, coincide (cf. Davis [1], Galdi and Straughan [2]). In more general situations, both theories have certain short comings. Thus, there is always the possibility of instability below the critical values predicted by linearized theory (often called subcritical instability or metastability) which has to be determined by other methods. On the other hand, the standard energy theory, since it determines the lower bound below which the disturbances are always stable, usually predicts very conservative results, sometimes quite low from the actual experimental values. In order to improve upon these results the generalized energy and Lyapunov methods are being proven attractive in th recent years. The books by Drazin and Reid [3] and Straughan [4] give very good expositions of both kind of theories. In the present paper we discuss the linear and nonlinear stability of convection induced by an applied temperature gradient inclined to the vertical in a shallow horizontal layer of a viscoelastic fluid. We represent the viscoelastic fluid by Jeffrey’s model [5]. In the case of a viscous fluid the above problem has been discussed by Weber [6], Nield [7] and Kaloni and Qiao [8]. While Weber and Nield gave results for the linear stability method, Kaloni and Qiao presented the nonlinear stability discussions. In the case of viscoelastic fluids the thermal convection problem of a viscoelastic fluid heated from below has been studied extensively by the linear stability method (see [9] and other references therein). The situation when the temperature gradient is inclined to the vertical has, however, only been studied more recently [10]. Correspondence to: P.N. Kaloni (e-mail: [email protected])

18

P.N. Kaloni, J.X. Lou

Our purpose here is to discuss the stability of the inclined temperature gradient problem in a viscoelastic fluid. We present results for both linearized and energy methods and compare them with each other. The eigenvalue calculations, in both cases, are carried out by the Chebyshev Tau-QZ method, following Dongarra et al. [11]. We discuss the effect of elasticity and compare results with the viscous and Maxwell fluid model. We also discuss the effect of the Prandtl number on the stability character in both cases. One of the conclusions drawn from our work is the possible occurence of subcritical instability (bifurcation), suggesting that the solution may become unstable at a value lower than that determined by the linear theory.

2 Governing equations We assume that the viscoelastic fluid is confined between two horizontal planes at a distance d apart. A Cartesian coordinate system (˜ x, y˜, z˜) is chosen such that origin is in the middle of the layer and z˜-axis is vertically upward. The fluid layer is taken to be confined in a cavity to be large in its horizontal extent, the x ˜y˜-plane, and small of the thickness d in the z˜ direction, It is also assumed that the ratio of the height to the length of the layer is considerably small so that the lateral end effects do not influence the motion in the horizontal central part. A constant horizontal temperature gradient β is imposed in the x ˜ direction, and a constant temperature difference ∆T is maintained between the two planes. On using the Boussinesq approximation, the basic equations take the form:


 ˜ ∂u ˆ ˜ u ˜ P˜ + ∇ ˜ · τ˜ − ρ0 1 − α T˜ − T0 g k, ˜ ·∇ ˜ = −∇ ρ0 + ρ0 u (1) ∂ t˜  ∂ T˜
˜ T˜ = κ∇ ˜ 2 T˜, ˜ ·∇ + u (2) ∂ t˜ ˜ ·u ˜ = 0. ∇ (3) ˜ the velocity vector, τ˜ the stress tensor, T˜ the temperature, α the coefficient of thermal Here, ρ0 is the density, u expansion, T0 is the ambient temperature and κ the thermal diffusivity. The boundary conditions are: d ∆T − βx ˜ at z = ± , T˜ = T0 − (±) 2 2 where β is a constant. The constitutive equation for Jeffrey’s viscoelastic fluid model is [5]:    ˙     ∂ τ˜
∂γ ˜
˜ ˜ ˜ ˜ ˙ ˜ · ∇ τ˜ = η γ ˜ ·∇ γ τ˜ + λ1 + u ˜ + λ2 + u ˜˙ , ∂ t˜ ∂ t˜ u ˜ = v˜ = w ˜ = 0,

where

(4)

(5)


T ˜ u + ∇˜ ˜u γ ˜˙ = ∇˜

˜1 = λ ˜ 2 and λ ˜ 2 = 0 respectively. We nonThe Newtonian and Maxwell models are obtained by setting λ dimensionalize the quantities as follows: x=

˜ x , d

u=

d ˜, u κ

t=

κt˜ d2

τ =

d2 τ˜ , ηκ

d2 ˜ κ˜ gαd3 ρ0 ˜ P = ε= (T − T0 ), P, λ = 2λ 1, ηκ ηκ d gαd3 ρ0 ∆T gαd4 ρ0 β η , RV = Pr = , RH = , κρ0 ηκ ηκ

T =

˜2 λ , ˜1 λ

where the coordinate vector x = (x, y, z), velocity vector u = (u, v, w), stress tensor τ = [(τxx , τxy , τxz ), (τxy , τyy , τyz ), (τxz , τyz , τzz )], g is the gravitational constant, P r is the Prandtl number, λ is the Deborah number, ε is the ratio of retardation to relaxation time, RV is the vertical Rayleigh number, and RH is the horizontal Rayleigh number. The non dimensional governing equations then take the form:   ∂u ,P r−1 + (u · ∇) u = −∇P + ∇ · τ + T k (6) ∂t

Nonlinear convection of a viscoelastic fluid with inclined temperature gradient

19

∂T + (u · ∇) T = ∇2 T, ∂t ∇ · u = 0,     ∂τ ∂ γ˙ τ +λ + (u · ∇) τ = γ˙ + ελ + (u · ∇) γ˙ , ∂t ∂t

(7) (8) (9)

where: T γ˙ = ∇u + (∇u) , and k = (k1 , k2 , k3 ) = (0, 0, 1). The non-dimensional forms of boundary conditions become:

u = v = w = 0,

T = − (±) RV /2 − RH x

at

(10) 1 z=± . 2

(11)

We look for a steady state solution of the form: us = u(z), v s = v(z), ws = 0, ps = p(x, y, z), τ s = τ s (z), T s = T (z) − RH x On substituting there expression in equations (6) - (9) and applying (11), we find that the basic steady state solution (us , T s , P s , τ s ) is given by: v s = 0, ws = 0, us = RH f (z), T s = −RH x − Rv z + RH 2 g(z), s τxz = RH s τxx

=

s τyy

df (z) , dz s s s = τzz = τxy = τyz = 0,

(12)

s

∇P = ∇ · τ s + T s k, where

1 1 (z − 4z 3 ); (7z − 40z 3 + 48z 5 ) g(z) = 24 5760 and where we have imposed the requirement that there is no net horizontal flux:  12  12 s u dz = v s dz = 0. f (z) =

− 12

− 12

3 Stability analysis We now perturb the steady-state solution as follows: u = us + u
, T = T s + θ, P = P s + P
, τ = τs + τ
On substituting (13) in (6) - (9) and dropping the primes on the fluctuating variables we obtain:   −1 ∂u Pr + (u · ∇) u = −∇P + ∇ · τ + θk ∂t − P r−1 [(us · ∇) u + (u · ∇) us ] , ∂θ + (u · ∇) θ = ∇2 θ − (us · ∇) θ − (u · ∇) T s , ∂t ∇ · u = 0,   ∂ ˙ − λ (us · ∇) (εγ˙ − τ ) λ (εγ˙ − τ ) + (u · ∇) (εγ˙ − τ ) = (τ − γ) ∂t − λ (u · ∇) (εγ˙s − τ s ),

(13)

(14) (15) (16)

(17)

where:

T γ˙s = ∇us + (∇us ) , and where us , T s and τ s are given by (12). The corresponding boundary conditions become: 1 u = 0, θ = 0, at z = ± . 2

(18)

20

P.N. Kaloni, J.X. Lou

3.1 Linear stability analysis We first discuss the linear stability analysis. In this case we linearize equations (14) -(17) with respect to the variables ( u, θ, P , τ ) and use the solenoidal vector format (φ, ψ) for u as:   2 ∂ φ ∂ψ ∂2φ ∂ψ ∂2φ ∂2φ . (19) (u, v, w) = + , − , − 2 − 2 ∂x∂z ∂y ∂y∂z ∂x ∂x ∂ y This representation helps to reduce the size of the matrix as noted below. We substitute expressions (19) into the linearized form of above equations (14)- (17), and take the vertical component of double curl and single curl of (14), together with (15) -(17). These equations in component form become:  
 ∂ 2
˜ 2  d2 f (z) ∂
˜ 2  ∂ 2
˜ 2  Rh ˜ 2 τzz ˜ 2θ + ∂ ∇ f (z) ∇ φ + ∇ ∇ ∇ φ + ∇ ∇ φ − ∂t Pr ∂x dz 2 ∂x ∂z     2 2 3 ∂ ∂ ∂ τ τ τ ∂ ∂ xz yz xx ˜ 2 τxz − ˜ 2 τyz − + ∇ ∇ + − ∂x ∂z 2 ∂y ∂z 2 ∂x2 ∂z ∂ 3 τxy ∂ 3 τyy −2 − 2 = 0, (20) ∂x∂y∂z ∂y ∂z   ∂
˜ 2  df (z) ∂
˜ 2  ∂ 2 τxy ∂ 2 τxy ∂
˜ 2  Rh ∇ ψ + −f (z) ∇ ψ + ∇ φ + − − ∂t Pr ∂x dz ∂y ∂y 2 ∂x2 ∂ 2 τxz ∂ 2 τxx ∂ 2 τyy ∂ 2 τyz + = 0, (21) + − − ∂x∂y ∂x∂y ∂x∂z ∂y∂z     2 ∂θ 2 dg (z) ˜ 2 φ + RH f (z) ∂θ − ∂ψ − ∂ φ = 0, ∇ (22) − ∇2 θ + Rv − RH ∂t dz ∂x ∂y ∂x∂z   ∂τxx ∂3ψ ∂3φ ∂2ψ ∂4φ ∂τxx λ + − 2 − 2 + τxx − 2λ + λRH f (z) 2 2 ∂t ∂x ∂z∂t ∂x∂y∂t ∂x ∂z ∂x∂y ∂x     4 df (z) ∂ 2 ψ ∂ φ ∂3ψ ∂3φ + 2λ (2 − 1) RH − 2λRH f (z) + + ∂x3 ∂z ∂x2 ∂y dz ∂y∂z ∂x∂z 2 df (z) df (z) ∂ ˜ 2 − 2λ τxz − 2λRH (23) ∇ φ = 0, dz dz ∂x   ∂3ψ ∂3φ ∂2ψ ∂4φ ∂τyy ∂τyy + τyy − 2λ − −2 2 +2 + λRH f (z) λ 2 ∂t ∂y ∂z∂t ∂x∂y∂t ∂y ∂z ∂x∂y ∂x   3 4 ∂ φ ∂ ψ − = 0, (24) + 2λRH f (z) ∂x2 ∂y ∂x∂y 2 ∂z ∂τzz ∂2
˜ 2  ∂
˜2  ∂τzz λ + τzz + 2λ ∇ φ +2 ∇ φ + λRH f (z) ∂t ∂z∂t ∂z ∂x df (z) ∂
˜ 2  ∂2
˜ 2  ∇ φ + 2λ (1 − ) RH ∇ φ = 0, + 2λRH f (z) (25) ∂x∂z dz ∂x ∂3ψ ∂3φ ∂4φ ∂τxy ∂3ψ ∂2ψ λ + τxy − 2λ + λ 2 − λ 2 − 2 + ∂t ∂x∂y∂z∂t ∂x ∂t ∂y ∂t ∂x∂y∂z ∂x2   ∂3ψ ∂2ψ ∂τxy ∂3ψ ∂4φ + 3 − 2 − + λRH f (z) − 2 2 2 ∂y ∂x ∂x ∂y∂z ∂x ∂x ∂y   3  
2 ∂ φ ∂ ˜2  df (z) ∂ ψ + λRH (2 − 1) ∇ −  (26) − φ − τ yz = 0 dz ∂y∂z 2 ∂x∂z 2 ∂y ∂τxz ∂3ψ ∂
˜2  ∂2
˜ 2  ∂4φ ∂3φ λ + τxz + λ ∇ φ − λ − λ + ∇ φ − ∂t ∂x∂t ∂x∂z 2 ∂t ∂y∂z∂t ∂x ∂x∂z 2    2 2
4 3 ∂ ∂ ψ ∂τxz ˜ 2φ −  ∂ φ −  ∂ ψ + λRH f (x) + 2 ∇ − ∂y∂z ∂x ∂x ∂x2 ∂z 2 ∂x∂y∂z   
 3 ∂2ψ df (z) ∂ ˜2 ∂ φ + λRH −τzz − 2 ∇ φ + (1 − ) − z ∂z ∂y 2 ∂z ∂x∂y

Nonlinear convection of a viscoelastic fluid with inclined temperature gradient

21

d2 f (z)
˜ 2  ∇ φ =0 dz ∂3ψ ∂
˜2  ∂2
˜ 2  ∂4φ ∂3φ ∂τyz + τyz + λ ∇ φ − λ − λ + ∇ φ − λ ∂t ∂y∂t ∂y∂z 2 ∂t ∂x∂z∂t ∂y ∂y∂z 2    2 2
4 3 ∂ ψ ∂τyz ∂ ˜ 2φ −  ∂ φ +  ∂ ψ ∇ + + λRH f (x) + ∂x∂z ∂x ∂x∂y ∂x∂y∂z 2 ∂x2 ∂z  2  df (z) ∂ ψ ∂3φ + (1 − ) λRH =0 − z ∂x2 ∂x∂y∂z 2 2 ∂2 ∂2 ∂2 ˜2 = ∂ + ∂ , where ∇2 = + + , and ∇ ∂x2 ∂y 2 ∂z 2 ∂x2 ∂y 2 We look for the solution of the above equations in the form − λ (1 − ) RH

T

[φ, ψ, θ, τxx , τyy , τzz , τxy , τxz , τyz ] = X exp [(kx + ly) i + σt] ,

(27)

(28)

(29)

T

where X = [φ(z), ψ(z), θ(z), τxx (z), τyy (z), τzz (z), τxy (z), τxz (z), τyz (z)] . To match the domain of the Chebyshev Tau-QZ method, we reset the present domain from [− 12 , 12 ] to [−1, 1] with coordinate transformation of z to 2z. Substituting the expression from (29) into equations (20) - (28), and applying the appropriate boundary conditions we obtain a set of eigenvalue equations. We next expand N +2 the solutions for different variables in a finite series of Chebyshev polynomials as (φ, ψ, θ, τxx , · · · ) = k=0 (φ, ψ, θ, τxx , · · · )Tk (z), where Tk (z) in the Chebyshev polynomial of degree k defined as: Tk (z) = cos(k arccosz), −1 ≤ z ≤ 1, k = 0, 1, 2, · · · .

(30)

On following the procedure outlined in [11,12], we then obtain the following eigenvalue system: σBX = AX, where X = [φ0 , · · · , φN +2 , · · · , (τxx )0 , · · · , (τxx )N +2 , · · · ], and matrices B, A are given:   −4α2 D 2 + α4 I 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0  0 α2 I    0 0 I 0 0 0 0 0 0     4λk 2 D 2λklI 0 λI 0 0 0 0 0    −2λklI 0 0 λI 0 0 0 0  4λl2 D B=    −4λα2 D 0 0 0 0 λI 0 0 0     4λklD λ(l2 − k 2 )I 0 0 0 0 λI 0 0     −iλk(4D 2 + α2 I) −i2λlD 0 0 0 0 0 λI 0  −iλl(4D 2 + α2 I) i2λkD 0 0 0 0 0 0 λI   A11 0 A13 A14 A15 A16 A17 A18 A19  A21 A22 0 A24 A25 0 A27 A28 A29     A31 A32 A33 0 0 0 0 0 0     A41 A42 0 A44 0 0 0 0 0     A=  A51 A52 0 0 A55 0 0 0 0  .  A61 A62 0 0 0 A66 0 0 0     A71 A72 0 0 0 0 A77 0 A79     A81 A82 0 0 0 A86 0 A88 0  A91 A92 0 0 0 0 0 0 A99 Here α2 = k 2 + l2 and variables Ai,j are given as   1 2 1 1 A11 = iRH k α (−Z 3 + Z)D 2 + α4 Z 3 − α2 (α2 − 24)Z , 12 48 48 A13 = P rα2 I,

A14 = −2P rk 2 D, A15 = −2P rl2 D, A16 = 2P rα2 D,  2    A1,7 = −4P rklD, A18 = iP rk 4D + α2 I , A19 = iP rl 4D 2 + α2 I ,     1 1 3 1 1 A22 = iRH kα2 Z − Z , A21 = iRH lα2 − Z 2 + I , 8 24 48 48

(31)

(32)

(33)

22

P.N. Kaloni, J.X. Lou

A24 = P rklI,

A25 = −P rklI,

A27 = P r(l2 − k 2 )I,

A28 = −2iP rlD, 1 A32 = −lRH I, A33 = 4D 2 − α2 I + i kRH (Z 3 − Z) A29 = 2iP rkD, 48   1 4 1 2 7 2 2 2 A31 = RH α − Z + Z + I + RV α I + i2kRH D, 384 192 5760 1 1 A41 = −4k 2 D + i λRH k 3 (Z 3 − Z)D, A42 = −2klI + i λRH k 2 l(Z 3 − Z), 12 24 1 1 3 2 2 A51 = −4l D + i λRH kl (Z − Z)D, A52 = 2klI + i λRH k 2 l(−Z 3 + Z), 24 24   1 1 3 2 2 A61 = 4α D + i λRH kα −Z + Z , A71 = −4klD + i λRH k 2 l(Z 3 − Z)D, 24 12 1 1 3 2 2 2 2 A72 = (k − l )I + i λRH (l − k )(Z − Z), A82 = RH kl(Z 3 − Z)D + i2lD, 48 24   1 1 1 (Z 3 − Z)D 2 + α2 (Z 3 − Z) + (1 − )λα2 RH Z + ik(4D 2 + α2 I), A81 = λRH k 2 12 48 2   1 1 1 A91 = λRH kl (Z 3 − Z)D 2 + α2 Z 3 − α2 Z + il(4D 2 + α2 I), 12 48 48   1 1 1 λRH k 2 (−Z 3 + Z)D − i2kD, A79 = A86 = λRH − Z 2 + I A92 = 24 8 24 1 A44 = A55 = A66 = A77 = A88 = A99 = −I + i λRH k(Z 3 − Z), 48 The boundary conditions become: dφ = θ = 0, at z = ±1 φ= dz In the above I is N × N the identity matrix, and N × N block matrices D and Z are defined as [11]: 2
ci π [< Tj , Ti >], 2 D 2 = [Di,j ] = ci2π [< Tj

, Ti >], m Z m = [Zi,j ] = ci2π [< Z m Tj , Ti >],

D = [Di,j ] =

(34)

i, j = 0, 1, 2, · · · , N + 2, i, j = 0, 1, 2, · · · , N + 2,

(35)

i, j = 0, 1, 2, · · · , N + 2, m = 1, 2, 3,

with the constants c0 = 2, cn = 1, for n = 1, 2, · · · . Moreover the dashes over T denote the derivatives and <, > represents the inner product defined in the weighted space L2 (−1, 1) as  1 fg √ < f, g >= dz (36) 1 − z2 −1 We follow the procedure and algorithms, as described in [10–12], and employ the LAPACK library for solving the generalized eigenvalue equations. We then follow the algorithum for determining neutral stability curves. The critical Rayleigh number RV C with critical wave number αc = kc2 + lc2 in this case, can be defined as: RV C = min RV (P r, RH , λ, ε, k, l) k,l

(37)

Numerical results are discussed in the next section. 3.2 Nonlinear stability analysis To discuss the nonlinear stability analysis, we define an energy function as follows: 1 β1 β2 ˙ 2  u 2 +  θ 2 +λ  (τ − εγ) E(t) = (38) 2P r 2 2 where β1 is the positive coupling parameter, and β2 is constant to be determined later on. On multiplying (14) ˙ and after using the boundary conditions and divergence theorem, we find by u, (15) by θ and (17) by (τ − εγ), 1 d  θ 2 = − < θ,i θ,i > − < ui T,si θ > (39) 2 dt

Nonlinear convection of a viscoelastic fluid with inclined temperature gradient

23

1 d 1  u 2 =< τij,j ui > + < θki ui > − < ui uj usi,j > 2P r dt Pr 1 d ˙ 2 = − < (γ˙ ij − τij )(εγ˙ ij − τij ) > λ  (τ − εγ) 2 dt s s − λ < uk (εγ˙ ij − τij ),k (εγ˙ ij − τij ) >

(40)

(41) 2

Here V denots a typical periodicity cell, < · > denotes the integration over V , and  ·  denots the L (V ) norm. The system of equations (39), (40) and (41), along with (38), can be put in the form [4]: 1 dE =I −D 2 dt

(42)

where 1 < ui uj usi,j > −β1 < ui T,si θ > Pr s s + β2 (1 + ε) < γ˙ ij τij > −β2 λ < uk (εγ˙ ij − τij ),k (εγ˙ ij − τij ) >,

I =< τij,j ui > + < θki ui > − 2

2

D = β1  ∇θ  +β2 ε  γ˙ ij  +β2  τij 

2

(43) (44)

Following [4], we define: I H D where H is the space of admissible solutions. On combining (42) - (45) we can infer m = max

dE ≤ −D(1 − m). dt which in turn, for 0 ≤ m ≤ 1, and on using the Poincar´ e inequality, becomes

(45)

(46)

dE (47) ≤ −2π 2 (1 − m) min(1, P r)E. dt Inequality (47) clearly indicates that for 0 ≤ m ≤ 1, E(t) → 0 at least exponentially as t → ∞. We now use the calculus of variation to find the maximum problem at the critical argument m = 1. The associated Euler–Lagrange equations (45) become: 1 s uk γ˙ ik 4β2 εui,jj + (1 − 2(1 + ε)β2 )τij,j + ki θ − − β1 θT,is Pr s s s s − λβ2 (εγ˙ mn − τmn )(εγ˙ mn − τmn ),i + 2λεβ2 (uk (εγ˙ im − τim ),k ),m = π,i (48) s ki ui − β1 ui T,i + 2β1 θ,ii = 0 (49) ui,i = 0 s s − ui,j + λβ2 uk (εγ˙ ij − τij ),k + (1 + ε)β2 (ui,j + uj,i ) − 2β2 τij = 0

(50) (51)

where π is the Lagrange multiplier introduced since u is solenoidal. We solve τij from (51) and substitute them into (48). As before we match the domain of Chebyshev Tau-QZ method, by resetting the present domain from [− 12 , 12 ] to [−1, 1] with coordinate transformation of z to 2z. Performing the standard normal mode analysis, by assuming the solution in the form: [u, v, w, θ, π] = [u(z), v(z), w(z), θ(z), π(z)] exp[i(kx + ly)]

(52)

and expanding the variables in finite series of Chebyshev polynomials, and proceeding exactly in the mammer as in the previous section, we obtain following eigenvalue system: R V B e X e = Ae X e where X e = [u, v, w, θ, π]T , and B e , Ae are given as:  00 0 0 0 0  Be =  0 0 0  0 0 β1 I 00 0

0 0 β1 I 0 0

 0 0  0  0 0

(53)

(54)

24

P.N. Kaloni, J.X. Lou



0 Ae11  0 Ae22  e 0 Ae =   A31  −β1 RH I 0 −ikI −ilI

 Ae13 −β1 RH I −ikI 0 0 −ilI   e e A33 A3,4 −2D   Ae4,3 Ae4,4 0  −2D 0 0

(55)

and where α2 = k 2 + l2 and variables Aei,j are given as   3 1 e 2 (4D 2 − α2 I) A11 = β2 (1 − ) − (1 + ) + 2 2β2   1 1 1 e 2 A13 = λRH β2 (1 − ) − (1 − ) (ZD + I) − RH Z 2 + RH I 2 8P r 24P r   1 3 Ae22 = β2 (1 − )2 − (1 + ) + (4D 2 − α2 I) 2 2β2 1 1 1 Ae31 = λ(1 − )RH [β2 (3 − 1) + 1] ZD − RH Z 2 + RH I 8P r 24P r  2 3 1 1 2 (4D 2 − α2 I) − β2 λ2 RH (1 − )2 Z 2 Ae33 = β2 (1 − )2 − (1 + ) + 2 2β2 4 1 + i β2 λRH k(1 − 2 )Z 4 1 e e 2 β1 RH A34 = A34 = (15Z 4 − 30Z 2 + 7)I − I 5760 Ae44 = −2β1 (4D 2 − α2 I) The boundary conditions are given as: u = v = w = θ = 0,

at

z = ±1

(56)

Note again that in the above I is an N × N identity matrix,and N × N block matrices D and Z are defined in the previous section and N is the order of Chebyshev polynomials. We now follow the procedure and algorithms, as described in [10–12], and employ the subroutine ZGGEV of LAPACK library for solving the generalized eigenvalue equations. We obtain a set of eigenvalue equations (53) for RV . They will be real since they come from a symmetric analysis. First, we consider the smallest positive real eigenvalue of the generalized eigenvalue equations (53) as RV . The critical energy vertical Rayleigh number RV is defined by RV = max min RV (RH , P r, λ, , k, l, β1 ) (β1 ) (k,l)

(57)

The maximization and minimization of (57) is carried out by means of the double Golden section search routine. We fix β2 = 12 so that the generalized eigenvalue equations (53) will reduce exactly to equations of the viscous fluid energy method, when  = 1. 4 Numerical results and discussion From [11,12] and our previous paper [10], we note the rapid convergence of Chebyshev polynomials. The calculation with N = 20 in the series of Chebyshev polynomials is thus satisfactory for our purpose, and we take, the order, N = 20 in our following calculations. We first report the results by the energy and linearized methods when RH = 0, and λ = 0.1. Figure 1 shows the plot between RV critical and  (elasticity parameter) for different Prandtl numbers. We note that in the viscous case ( = 1.0) the results are independent of P r and both the energy and linear stability methods give an identical critical value for RV . For the viscoelastic fluid case, we find that while in the linear stability method, results are dependent on P r, they do not depend in P r in the energy method. Also in the linear stability theory, the critical RV decreases as  decreases and also critical RV decreases as P r increases. The energy calculations always give lower critical RV values, and as  decreases, the critical RV decreases at a much faster rate than in the linear stability method. In the Maxwell model case, for example, the critical RV for P r = 100, in the linear stability case, is 226.685, it is zero and independent of P r in the energy stability case.

Nonlinear convection of a viscoelastic fluid with inclined temperature gradient

25





#'





3U OLQHDU



3U OLQHDU 3U OLQHDU (QHUJ\  











ε

Fig. 1. Variation of critical RV against
for different P r (RH = 0, λ = 0.1) /LQHDU3U  /LQHDU3U  /LQHDU3U  /LQHDU3U  (QHUJ\3U  (QHUJ\3U  (QHUJ\3U  (QHUJ\3U 



59





 



















5+

Fig. 2. Variation of critical RV with RH for different P r (λ = 0.1,
= 1.0 and θ = 90)



26

P.N. Kaloni, J.X. Lou ε=1.0



ε=0.8 ε=0.6 ε=0.4 ε=0.2 ε=0.1



59







 





















5+

Fig. 3. Critical RV with RH for different
(linear method, P r = 10, λ = 0.1 and θ = 90) ε=1.0 

ε=0.8 ε=0.6 ε=0.4 ε=0.2 ε=0.1

59







 



















5+

Fig. 4. Critical RV with RH for different
(energy method, P r = 10, λ = 0.1 and θ = 90)



Nonlinear convection of a viscoelastic fluid with inclined temperature gradient

27

Figure 2 shows the variation of critical RV against RH , the horizontal Rayleigh number, for different values of P r. In this case we note that in the case of the linear stability method, whereas the flow is stabilized (RV increases) with the increase of RH , in the case of the energy method the critical RV increases or remains constant up to RH = 500, and then it decreases as RH increases. In the energy method, for higher values of RH , the convection is thus carried out by the horizontal temperature gradient only. In both cases, we also note that the flow is destabilized by increasing P r number. Figures 3 and 4 give results, for different values of viscoelastic parameter  for P r = 10 and λ = 0.1. It is easily observed from Fig. 3, that the stabilizing behaviour of the viscous fluid case as RH increases diminishes as the elastic effect increases. In fact for  = 0.6 the critical RV remains almost constant as RH increases (up to RH = 2000). However, as the elastic effect is further increased (toward the Maxwell model) RV starts decreasing with RH beyond the RH = 1000 value. This clearly indicates the destabilizing effect of elasticity. Figure 4 shows the corresponding results for the energy method. Here we find that an increase in RH with decreases in , always has a destabilizing effect. As in the case of linearized method, the effect of the increase of elasticity is, of course, destabilizing.

5 Concluding remarks We remark that the present work has, in view of different RV critical values in the two methods, demonstrated the possibility of subcritical instability. The values obtained for RV critical, by two methods, are in some instances far apart, and thus it is likely that an intermediate value may be the right one. It is also possible that by suitably choosing a different energy function E(t), i.e., by employing the generalized energy method or Lyapunov methods, we may be able to obtain the two values closer to each other. A 3-D numerical study may, however, be needed to clearly resolve the issue. Acknowledgements. The work reported in the paper is supported by Grant No A7728 of NSERC of Canada. The authors are grateful for the support thus received. They also thank Prof. Straughan for helpful comments.

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.

Davis, S.H.: On the priciple of exchange of stabilities. Proc. Roy Soc. Lond. A 310, 341 (1969) Galdi, G.P., Straughan, B.: Exchange of stabilities, symmetry and non-linear stability. Arch. Rat. Mech Anal. 89, 211 (1985) Drazin, P.G., Reid, W.H.: Hydrodynamic stability. Cambridge University Press (1982) Straughan, B.: The Energy Method, Stability, and Nonlinear Convection, 2nd edn, Applied Mathematical Sciences 91, (2004) Springer-Verlag Bird, R.B., Hassager, O., Armstrong, R.C., Curtis, C.F.: Dynamics of Polymeric Liquids, Volume 1, Fluid Mechanics. Wiley, New York (1987) Weber, J.E.: On the stability of thermally driven shear flow heated from below J. Fluid Mech. 87, 65 (1978) Nield, D.A.: Convection induced by an inclined temperature gradient in a shallow horizontal layer. Int. J. Heat and Fluid Flow 15, 157 (1994) Kaloni, P.N., Qiao, Z.: On the nonlinear stability of thermally driven shear flow heated from below. Physics of Fluids 8, 639 (1996) Kolkka, R.W., Ierley, G.R.: On the convected linear stability of a viscoelastic Oldroyd B fluid heated from below. J. NonNewtonian Fluid Mechanics 25, 209 (1987) Kaloni, P.N., Lou, J.X.: Stability of Hadley circulations in a Maxwell fluid. J. Non-Newtonian Fluid Mechanics 103, 169 (2002) Dongarra, J.J., Straughan, B., Walker, D.W.: Chebyshev tau-QZ algorithm methods for calculating spectra of hydrodynamic stability problems. Applied Numerical Mathematics 22, 399 (1996) Straughan, B., Walker, D.W.: Two very accurate and efficient methods for computing eigenvalues and eigenfunctions in porous convection problems. J. Computat. Physic. 127, 128 (1996)