Acta Mech 223, 967–983 (2012) DOI 10.1007/s00707-012-0616-1

M. S. Malashetty · A. A. Hill · Mahantesh Swamy

Double diffusive convection in a viscoelastic fluid-saturated porous layer using a thermal non-equilibrium model

Received: 25 May 2011 / Published online: 7 February 2012 © Springer-Verlag 2012

Abstract The stability of a binary viscoelastic fluid-saturated porous layer which is heated from below is studied, where the fluid and solid phases are not in local thermal equilibrium. The modified Darcy-Oldroyd model is employed as a momentum equation, with the fluid and solid phase temperature fields modelled separately. It is found that the inter-phase heat transfer coefficient has a significant effect on the stability of the system. Competition between the processes of viscous relaxation and thermal diffusion cause the convection to set in through oscillatory rather than stationary instability, with the viscoelastic parameters inhibiting the onset of convection. In the case of weakly non-linear theory, both steady and unsteady cases are considered. In the unsteady case the transient behaviour of the Nusselt and Sherwood numbers is investigated. The effect of thermal non-equilibrium on heat and mass transfer is also discussed.

1 Introduction The problem of double diffusive convection in porous media has attracted considerable interest during the last few decades due to its wide range of applications, from the solidification of binary mixtures to the migration of solutes in water-saturated soils. Other important examples include geophysical systems, electro-chemistry and the migration of moisture through air contained in fibrous insulation. A comprehensive review of the literature concerning double diffusive natural convection in fluid-saturated porous media may be found in the books [1–4]. In most investigations a state of local thermal equilibrium (LTE) between the fluid and solid phases is assumed. This is where the temperature gradient at any location between the two phases is assumed to be negligible. For many practical applications, such as the drying and freezing of foods, microwave heating and rapid heat transfers, the assumption of local thermal equilibrium is inadequate, making the motivation behind the study of local thermal non-equilibrium (LTNE) effects highly relevant [2]. An extensive review of work on the LTNE model is given in [1–3]. More recent contributions include [5–11]. This work is supported by UGC New Delhi, under the Special Assistance Programme DRS Phase II. M. S. Malashetty Department of Mathematics, Gulbarga University, Jnana Ganga Campus, Gulbarga 585106, India E-mail: [email protected] Tel.: +91-08472-245633/250086 Fax: +91-08472-245927 A. A. Hill School of Mathematical Sciences, University of Nottingham, Nottingham NG7 2RD, UK M. Swamy (B) Department of Mathematics, Government College, Sedam Road, Gulbarga 585105, India E-mail: [email protected]

968

M. S. Malashetty et al.

Bories [12] introduced a two-field model for the energy equation and the same is used by Nield and Bejan [2] (and Nield and Kuznetsov [13]) as the simplest way in which LTNE may be modelled. Instead of having a single energy equation, which describes the common temperature of the saturated porous media, two equations are utilised for the fluid and solid phases separately. In this two-field model, the energy equations are coupled by terms which account for the heat transfer between the phases. Viscoelastic fluid flow in porous media is of interest for many engineering fields such as enhanced oil recovery, paper and textile coating, and composite manufacturing processes. A proper description of the viscoelastic fluid flow through fibrous media is a very important input in the modelling of composite manufacturing processes. The general linear viscoelastic model of fluid can describe some of the time-dependent motions of polymeric fluids, although it describes a restricted class of flows with very small displacement gradients. Modelling of polymeric flow in porous media has motivated many studies that essentially focus on the numerical simulation of viscoelastic fluid flow in a specific pore geometry model, such as capillary tubes, undulating tubes, packs of spheres or cylinders. A good review of these studies can be found in Skartsis et al. [14]. Recently, interest in viscoelastic flows through porous media has grown considerably, due largely to the demands of such diverse fields as biorheology, geophysics, chemical and petroleum industries. Convective instability thresholds for viscoelastic fluids [9,15], binary viscoelastic fluids [16–19] and viscoelastic fluids in porous media [20–25] are also widespread in the literature. It has been demonstrated that the arising instability may be stationary or oscillatory, depending on the binary and viscoelastic properties of the fluid. At present, the convective instability problem for a binary viscoelastic fluid in a porous medium has not been given much attention. Recently, Wang and Tan [26] performed a stability analysis of double diffusive convection of Maxwell fluid in a porous medium, and Malashetty et al. [9] performed a linear and non-linear analysis of double diffusive convection in a viscoelastic fluid-saturated porous medium. In this paper, we perform linear and weakly non-linear stability analyses of double diffusive convection in a binary viscoelastic fluid-saturated porous layer, with the assumption that the fluid and solid phases are not in local thermal equilibrium. Our objective in this paper is to study how the onset criterion for steady, oscillatory and finite amplitude modes, and heat and mass transfer, are affected by local thermal non-equilibrium condition. 2 Mathematical formulation Consider a viscoelastic fluid occupying the three-dimensional layer {(x, y) ∈ R2 } × {z ∈ (0, d)}, with a vertically downward gravity force g acting upon it. A uniform adverse temperature gradient T = Tl − Tu , for temperature T , and a stabilising concentration gradient S = Sl − Su , for concentration S, where Tl > Tu and Sl > Su are maintained between the lower and upper surfaces. A Cartesian reference frame is chosen with the origin in the lower boundary and the z-axis vertically upwards. The Cauchy stress tensor T for an Oldroyd-B fluid is given by the constitutive equation [27]   DS DA T = − pI + S, S + λ1 = μ A + λ2 , (1) Dt Dt where p is the hydrostatic pressure, I the identity tensor, μ the viscosity of the fluid, S the extra stress tensor, and λ1 and λ2 are constant relaxation and retardation times, respectively. A = ∇q + ∇qT is the strain-rate tensor, q is the velocity vector, and   DS ∂ (2) = + q · ∇ S − S(∇q) − (∇q)T S, Dt ∂t   DA ∂ (3) = + q · ∇ A − A(∇q) − (∇q)T A. Dt ∂t It should be noted that this model includes the classical viscous Newtonian fluid as a special case for λ1 = λ2 = 0, and a Maxwell fluid when λ2 = 0. It is well known that in the flow of viscous Newtonian fluid at a low speed through a porous medium the pressure drop caused by the frictional drag is directly proportional to velocity, which is Darcy’s law. By analogy with Oldroyd-B constitutive relationships, the following phenomenological model, which relates pressure drop and velocity for a viscoelastic fluid in an anisotropic porous medium, is given by [28]     μ ∂ ∂ 1 + λ1 (4) ∇p = − 1 + λ2 qD , ∂t K ∂t

Double diffusive viscoelastic porous convection

969

where K is the permeability of the porous medium and q D is the Darcian velocity which is related to the usual (i.e. volume averaged over a volume element consisting of fluid only in the pores) velocity vector q by q D = εq, with ε being the porosity. We note that for λ1 = λ2 = 0, (4) simplifies to Darcy’s law for the flow of a viscous Newtonian fluid through an isotropic porous medium. Thus, (4) can be regarded as an approximate form of an empirical momentum equation for the flow of Oldroyd-B fluid through an isotropic porous medium. Under consideration of the balance of forces acting on a volume element of fluid, the local volume average balance of linear momentum is given by dq ρ = −∇ p + ρgk + ∇ · S + r, (5) dt where d/dt is the material time derivative, ρ is the density, k = (0, 0, 1), and r is the Darcy resistance for an Oldroyd-B fluid in the porous medium. Since the pressure gradient in (4) can also be interpreted as a measure of the resistance to flow in the bulk of the porous medium, and r is a measure of the flow resistance offered by the solid matrix, r can be inferred from (4) to satisfy     με ∂ ∂ 1 + λ1 r=− 1 + λ2 q. (6) ∂t K ∂t Substituting (6) into (5), we obtain      ∂ ∂ dq με 1 + λ1 ρ + ∇ p − ρgk − ∇ · S = − 1 + λ2 q. (7) ∂t dt K ∂t For the Darcy model, neglecting the inertia term q · ∇q and the viscous term ∇ · S, (7) can be simplified (after dropping the suffix D on q for simplicity)      ρ0 ∂q μ ∂ ∂ + ∇ p − ρgk = − 1 + λ2 q (8) 1 + λ1 ∂t ε ∂t K ∂t with reference density ρ0 . In modelling the energy equation for a fluid-saturated porous system, there are two different theoretical approaches which can be used. In the first model, the fluid and solid structures are assumed to be in local thermal equilibrium. This assumption is satisfactory for small-pore media such as geothermal reservoirs and fibrous insulations, where the temperature differences between the phases are small. In the second approach, the fluid and solid structures are assumed to be in thermal non-equilibrium, where there is large temperature difference between the fluid and solid phases. If the temperatures difference between phases is a very important safety parameter (e.g. fixed bed nuclear propulsion systems and nuclear reactor modelling), the thermal non-equilibrium model is indispensable. Adopting a thermal non-equilibrium model, the fluid and solid phase temperatures are modelled separately [2], such that ∂Tf ε(ρ0 c) f (9) + (ρ0 c) f (q · ∇)T f = εk f ∇ 2 T f + h(Ts − T f ), ∂t ∂ Ts (1 − ε)(ρ0 c)s (10) = (1 − ε)ks ∇ 2 Ts − h(Ts − T f ), ∂t where c is the specific heat and k is the thermal conductivity, with the subscripts f and s denoting fluid and solid phase, respectively. In this two-field model, the energy equations are coupled by the inter-phase heat transfer coefficient h, which accounts for the heat lost to, or gained from, the other phase. This coefficient depends on the nature of the porous matrix and the saturating fluid and has been the subject of intense experimental interest. Large values of h correspond to a rapid transfer of heat between the phases (LTE), whereas small values of h give rise to relatively strong thermal non-equilibrium effects. In Eqs. (9)–(10), T f and Ts are intrinsic averages of the temperature fields which allows one to set T f = Ts = Tb whenever the boundary of the porous medium is maintained at the temperature Tb . The equations of continuity, state and solute concentration are ∇ · q = 0, (11)   (12) ρ = ρ0 1 − βT (T f − Tl ) + β S (S − Sl ) , ∂S (13) ε + (q · ∇) S = χ ∇ 2 S, ∂t where χ , βT and β S denote solute diffusivity and thermal and solutal coefficients of expansion, respectively.

970

M. S. Malashetty et al.

Assuming the basic state is quiescent, we superimpose infinitesimal perturbations on it (denoted by  s) and eliminate the pressure term by taking the double curl of (8). By adopting the following transformations: x = dx∗ , t =

(ρ0 c) f d 2 ∗ εk f t , q = q∗ , kf (ρ0 c) f d

T f = (T )T f∗ , Ts = (T )Ts∗ , S  = (S)S ∗ , the non-dimensional perturbation equations (dropping the asterisks for simplicity) are given by      1 ∂ ∂ ∂ (∇ 2 w) − RaT ∇12 T f + Ra S ∇12 S + 1 + λ2 ∇ 2 w = 0, 1 + λ1 ∂t Pr D ∂t ∂t   ∂ − ∇ 2 T f + (q · ∇) T f = w + H (Ts − T f ), ∂t   ∂ α − ∇ 2 Ts = γ H (T f − Ts ), ∂t   ∂ 1 2 − ∇ S + (q · ∇) S − w = 0, ∂t Le

(14) (15) (16) (17)

where ∇12 = ∂ 2 /∂ x 2 + ∂ 2 /∂ y 2 , λ1 = κ f λ1 /d 2 and λ2 = κ f λ2 /d 2 are the dimensionless relaxation and retardation parameters, respectively, with κ f = k f /(ρ0 c) f , Pr D = εμd 2 /(ρ0 κ f K ) is the Darcy-Prandtl number, H = hd 2 /(εk f ) is the dimensionless inter-phase heat transfer coefficient, α = (ρ0 c)s κ f /ks is the ratio of diffusivities, γ = εk f /((1 − ε)ks ) is the porosity modified conductivity ratio, and Le = κ f /χ is the Lewis number. The remaining parameters are the thermal and solute Rayleigh numbers RaT =

ρ0 βT gd(T )K ρ0 β S gd(S)K , Ra S = . εμκ f εμκ f

It is worth mentioning that the Rayleigh number RaT defined above is based on the properties of the fluid whilst = RaLTE T

γ ρ0 βT gd(T )K , RaT = 1+γ μ(εk f + (1 − ε)ks

is the Rayleigh number based on the mean properties of the porous medium, and it is this latter definition, which is used in the thermal equilibrium model. Since the fluid and solid phases are not in thermal equilibrium, the use of appropriate thermal boundary conditions may pose a difficulty. However, the assumption that the solid and fluid phases have equal temperatures at the bounding surfaces helps in overcoming this difficulty. Accordingly, Eqs. (14)–(17) are solved for impermeable isothermal and isohaline boundaries. Hence the boundary conditions for the perturbation variables are given by w = T f = Ts = S = 0 at z = 0, 1.

(18)

The objective of this study is to investigate the effect of three important parameters, namely the two viscoelastic parameters λ1 and λ2 and the interface heat transfer coefficient H , which characterises the absence of local thermal equilibrium in a fluid-saturated porous medium. The viscoelastic character of the liquid mixture appears in the relaxation parameter λ1 (which is also known as the Deborah number) and the retardation parameter λ2 . The Deborah number is a dimensionless number used in rheology to characterise how fluid a material is. The smaller the Deborah number, the more fluid the material appears. The parameter λ2 is zero for a Maxwell fluid, whilst λ2 = λ1 for a Newtonian fluid. The parameter λ1 that relates the relaxation time to the thermal diffusion time is of order one for most viscoelastic fluids. The value for the Deborah number for a dilute polymeric solution is most likely in the range [0.1, 2] and the parameter λ2 in the range [0.1, 1]. At present there are no experimental data available for comparison and therefore we have considered a wide range of values for the parameters. In the absence of local thermal equilibrium between fluid and solid phases in a porous medium, the interface heat transfer coefficient H is the governing parameter. When H → 0 there is almost no transfer of heat

Double diffusive viscoelastic porous convection

971

between the fluid and solid phases, and the properties of solid phase have no significant influence on the onset criterion. When H → ∞ the fluid and solid phases have almost equal temperatures and therefore may be treated a single phase. Between these two extremes H gives rise to a strong non-equilibrium effect. Therefore, a wide range of values for the parameter H is considered, including the two limiting values. 3 Linear instability analysis In this section we predict the thresholds of both marginal and oscillatory convections using linear theory. The eigenvalue problem defined by Eqs. (14)–(17), subject to the boundary conditions (18), is solved using time-dependent periodic disturbances in a horizontal plane, upon assuming that amplitudes are small enough and can be expressed as ⎛ ⎞ ⎛ ⎞ w W (z) ⎜ T f ⎟ ⎜ (z) ⎟ (19) ⎝ T ⎠ = ⎝ (z) ⎠ exp[i(lx + my) + ωt], s (z) S where l and m are the horizontal wavenumbers and ω is the growth rate. Infinitesimal perturbations of the rest state may either damp or grow depending on the value of the parameter ω. Substituting the normal modes (19) into the linearised version of Eqs. (14)–(17) we obtain   ω 1 + λ2 ω + (20) (D 2 − a 2 )W + a 2 RaT − a 2 Ra S  = 0, Pr D 1 + λ1 ω ((D 2 − a 2 ) − ω) + W + H ( − ) = 0, ((D 2 − a 2 ) − αω) + γ H ( − ) = 0,   1 (D 2 − a 2 ) − ω  + W = 0, Le

(21) (22) (23)

where D = d/dz and a 2 = l 2 + m 2 . The boundary conditions (18) now become W = θ = =  = 0 at z = 0, 1.

(24)

We assume the solution to W, , and  of the form W = W0 sin π z, = 0 sin π z,

= 0 sin π z,  = 0 sin π z, which satisfy the boundary conditions (24). Substituting this form into Eqs. (20)–(23) we obtain the matrix equation ⎛ ⎞⎛ ⎞ ⎛ ⎞ M11 −a 2 RaT 0 a 2 Ra S W0 0 ⎜ −1 ω + δ 2 + H ⎟ ⎜ 0 ⎟ ⎜ 0 ⎟ −H 0 ⎜ ⎟⎝ ⎝ 0 ⎠ 0 ⎠ = ⎝ 0 ⎠ , 0 −γ H αω + δ 2 + γ H 0 0 −1 0 0 ω + δ 2 Le−1 where

 M11 =

 ω 1 + λ2 ω 2 + δ Pr D 1 + λ1 ω

and δ 2 = π 2 + a 2 is the total wavenumber. For this matrix equation to have a non-trivial solution, we require   H (αω + δ 2 ) RaT = ω + δ 2 + αω + δ 2 + γ H   2  1 + λ2 ω Ra S δ ω + + × . (25) a 2 Pr D 1 + λ1 ω ω + δ 2 Le−1 The growth rate ω is, in general, a complex quantity such that ω = ωr +iωi . For ωr > 0 the system is unstable. For neutral stability ωr = 0.

972

M. S. Malashetty et al.

3.1 Stationary convection At the monotonic instability threshold ωr = ωi = 0, (25) becomes RaTSt =

(δ 2 + H (1 + γ ))(δ 4 + a 2 LeRa S ) , a 2 (δ 2 + γ H )

(26)

which is identical to the result of Malashetty et al. [8], who studied the case of a Newtonian viscous fluid saturating the porous layer. As Ra S → 0, Eq. (26) yields RaTSt

δ4 = 2 a



 δ 2 + H (1 + γ ) . δ2 + γ H

This result was obtained by Banu and Rees [29] for the local thermal non-equilibrium model, in the absence of a salt field.

3.2 Oscillatory convection At the oscillatory instability threshold ωr = 0 and ωi  = 0, (25) becomes RaOsc T = 1 + iωi 2 ,

(27)

where 1 = b1 + ωi2 b2 and 2 = b3 ωi6 + b4 ωi4 + b5 ωi2 + b6 , with b1 to b6 given in the Appendix. Since RaT is a physical quantity, it must be real. Thus, either ωi = 0 (steady onset) or 2 = 0 (ωi  = 0, oscillatory onset). Combining 2 = 0 with (27) yields RaOsc T

=δ

δ +

2

2

H (δ 4 + δ 2 γ H + α 2 ωi2

(δ 2 + γ H )2 + α 2 ωi2



1 + λ1 λ2 ωi2 αγ H 2 LeRa S 2 + ωi 1 + 2 + 2 × 4 a (1 + λ21 ωi2 ) (δ + γ H )2 + α 2 ωi2 δ + Le2 ωi2

δ2 1 λ1 − λ2 Le2 Ra S . × 4 + − a 2 1 + λ21 ωi2 Pr D δ + Le2 ωi2

(28)

The critical oscillatory Rayleigh number is to be found from (28) by minimising over a 2 , using solutions to 2 = 0 to provide ωi2 . Since 2 is cubic in ωi2 , it can give rise to more than one positive value of ωi2 for fixed values of the parameters. This has an important bearing on the linear instability of a doubly diffusive porous layer saturated with an Oldroyd-B fluid.

4 Weakly non-linear analysis In this section we perform a non-linear analysis using a truncated representation of Fourier series with two terms. Although the linear stability analysis is sufficient for obtaining the instability condition of the motionless solution, with the corresponding eigenfunctions describing qualitatively the convective flow, it does not provide information about the convection amplitudes and the rate of heat transfer. To gain some understanding of the physical mechanism (with the minimum amount of mathematical analyses) we perform a non-linear analysis, which is a step forward towards understanding the full non-linear problem. To simplify the analysis, we confine ourselves to the two-dimensional rolls, so that all the physical quantities are independent of y. Introducing a stream function ψ, such that u  = ∂ψ/∂z and w  = −∂ψ/∂ x

Double diffusive viscoelastic porous convection

973

into Eqs. (8)–(13), eliminating pressure from the momentum equation and non-dimensionalising (using the transforms from Sect. 2 and dropping the asterisks) we obtain the equations      ∂Tf 1 ∂ 2 ∂ ∂S ∂ 1 + λ1 ∇ ψ + RaT − Ra S + 1 + λ2 ∇ 2 ψ = 0, ∂t Pr D ∂t ∂x ∂x ∂t   ∂ ∂(ψ, T ) ∂ψ 2 − ∇ Tf − + = H (Ts − T f ), ∂t ∂(x, z) ∂x   ∂ α − ∇ 2 Ts = γ H (T f − Ts ), ∂t   ∂ 1 2 ∂(ψ, S) ∂ψ − ∇ S− + = 0. ∂t Le ∂(x, z) ∂x

(29) (30) (31) (32)

The linear stability analysis is valid only for infinitesimal perturbations. However, for the flows with significantly large Rayleigh number, i.e. RaT > RaTSt , the temperature gradient will be considerably large. In such situation the thermal perturbations will be of finite amplitude. The linear stability theory cannot be employed for the finite amplitude motions. Therefore one has to take into account the non-linear effects. The first effect of non-linearity is to distort the temperature and concentration fields through the interaction of ψ with T f , Ts and S. The distortion of these fields corresponds to a change in the horizontal mean, i.e. a component of the form sin(2π z) will be generated. Thus, a minimal Fourier series which describes the finite amplitude free convection is given by ψ Tf Ts S

= C1 (t) sin(ax) sin(π z), = C2 (t) cos(ax) sin(π z) + C3 (t) sin(2π z), = C4 (t) cos(ax) sin(π z) + C5 (t) sin(2π z), = C6 (t) cos(ax) sin(π z) + C7 (t) sin(2π z),

where the amplitudes C1 , . . . , C7 are determined from the dynamics of the system. Substituting these definitions into Eqs. (29)–(32), we obtain the following non-linear autonomous system of differential equations dC1 dt dM dt dC2 dt dC3 dt dC4 dt dC5 dt dC6 dt dC7 dt

= M, =−

1 λ1

 M+

 Pr D 2 (δ (C + λ M) + aRa (C + λ N ) − aRa (C + λ N )) , 1 2 T 2 1 1 S 6 1 2 δ2

= N1 = −(aC1 + δ 2 C2 + πaC1 C3 + H (C2 − C4 )), = −4π 2 C3 +

πa C1 C3 + H (C5 − C3 ), 2

1 (−δ 2 C4 + γ H (C2 − C4 )), α 1 = (−4π 2 C5 + γ H (C3 − C5 )), α δ2 = N2 = −aC1 − C6 − πaC1 C7 , Le πa 4π 2 =− C7 + C1 C6 . Le 2

(33)

=

Although system (33) must be solved numerically, we can make some qualitative predictions. The system of equations (33) is uniformly bounded in time and possesses many properties of the full problem. Similarly to (14)–(17), system (33) must be dissipative. Thus, volume in the phase space must contract.

974

M. S. Malashetty et al.

In order to prove volume contraction, we must show that the velocity field has a constant negative divergence. It follows that         7 dC j ∂ dM 1 + Pr D λ2 ∂ 1 + 1+ + =− (δ 2 + 4π 2 ) ∂ M dt ∂C j dt λ1 Le j=1  1 + (δ 2 + 4π 2 + 2H (α + γ )) , α which is always negative and therefore the system is bounded and dissipative. As a result, the trajectories are attracted to a set of measure zero in the phase space; in particular they may be attracted to a fixed point, a limit cycle or, perhaps, a strange attractor. From the above expression we conclude that if a set of initial points in phase space occupies a region V (0) at time t = 0, then after some time t, the end points of the corresponding trajectories will fill a volume       1 + Pr D λ2 1 1 2 2 2 2 V (t) = V (0) exp − + 1+ (δ + 4π ) + (δ + 4π + 2H (α + γ )) t . λ1 Le α It is clear from this expression that the volume decreases with time. In the case of steady motions, system (33) can be solved in closed form by setting the time differentials to be 0. Rearranging in terms of C1 yields −1  2 2 δ (δ + γ H + H ) a 2 x(4π 2 + γ H ) a 2 LeRa S 2 2 = 0, + − a RaT C1 δ + 2 δ2 + γ H 4π 2 + γ H + H δ + a 2 Le2 x where x = C12 /8. The conduction state solution is obtained by setting C1 = 0. If C1  = 0 then we have A1 x 2 + A2 x + A3 = 0. By letting the radical in the required root x vanish, we obtain an expression for the finite amplitude Rayleigh number    1 F 2 −B2 + B2 − 4B1 B3 , Ra = 2B1 which characterises the onset of finite amplitude steady motions. The terms A1 , A2 , A3 , B1 , B2 and B3 are given in the Appendix. 4.1 Heat and mass transport As the Rayleigh number is increased, the onset of convection is more readily detected by its effect on the heat and mass transport. In the basic state, heat and mass transport is by conduction alone. If H˜ and J˜ are the rate of heat and mass transport per unit area, respectively, then     ∂(T f )total ∂ Stotal ˜ ˜ H = −κ f , J = −χ , ∂z ∂z z=0 z=0 where the angular bracket corresponds to a horizontal average and T z Sz (T f )total = Tl − + T f (x, z, t), Stotal = Sl − + S(x, z, t). d d Substituting the minimal Fourier series solutions for T f and S used in the previous section into these definitions yields κ f T χ S H˜ = J˜ = (1 − 2πC3 ), (1 − 2πC6 ). d d The Nusselt and Sherwood numbers are defined by d J˜ d H˜ Sh = = 1 − 2πC3 , = 1 − 2πC6 . Nu = κ f T χ S To understand the transient behaviour of the Nusselt and Sherwood numbers, the autonomous system (33) has been solved numerically using the Runge-Kutta method with appropriate initial conditions.

Double diffusive viscoelastic porous convection

975

5 Results and discussion Analytical expressions for the stationary, oscillatory and finite amplitude critical Rayleigh numbers, which characterise the stability of the system, were obtained in Sects. 3 and 4 and are presented graphically in this section for a range of parameter values. Due to the large parameter space, only the most interesting parameter interactions are presented in this section. Unless stated otherwise, the parameters for Figs. 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10 are H = 100, γ = 0.5, α = 0.25, λ1 = 0.7, λ2 = 0.5, Ra S = 100, Le = 2.5 and Pr D = 2.5. Figures 1, 2, 3, 4 and 5 show the neutral curves for different values of the viscoelastic parameters, interface heat transfer coefficient, thermal diffusivity ratio and solute Rayleigh number. From these figures it is clear that the neutral curves are connected in a topological sense. This connectedness allows the linear stability criteria to be expressed in terms of the critical Rayleigh number, below which the system is stable and unstable above. The points where the overstable solutions branch off from the stationary convection can be easily identified from these figures. We also observe that for smaller values of the wavenumber each curve is a margin of the oscillatory instability and, at some fixed wavenumber, the overstability disappears and the curve forms the margin of stationary convection. The effect of relaxation parameter λ1 on the neutral curves is shown in Fig. 1. We observe that the effect of increasing the relaxation parameter decreases the minimum of the Rayleigh number, indicating that the effect of the relaxation parameter is to advance the onset of convection. Figure 2 shows that an increase in the value of the retardation parameter increases the critical oscillatory Rayleigh number. Furthermore, we observe from Figs. 1 and 2 that all the curves lie below of the case where the fluid in the porous medium is Newtonian, indicating that the oscillatory convection sets in earlier in a system with viscoelastic fluid than a Newtonian fluid system. 800

700

600

500

400

300

0

1

2

3

4

5

6

5

6

Fig. 1 Neutral stability curves for different values of stress relaxation parameter λ1 800

700

600

500

400

300

0

1

2

3

4

Fig. 2 Neutral stability curves for different values of stress relaxation parameter λ2

976

M. S. Malashetty et al.

700

600

500

400

300

200

0

1

2

3

4

5

6

7

Fig. 3 Neutral stability curves for different values of inter-phase heat transfer coefficient H 900

800

700

600

500

400

300

0

2

4

6

8

Fig. 4 Neutral stability curves for different values of diffusivity ratio α 800 700 600 500 400 300 200 100

0

1

2

3

4

5

6

Fig. 5 Neutral stability curves for different values of solute Rayleigh number Ra S

Figure 3 explores the effect of the inter-phase heat transfer coefficient H . We find that for the parameter range chosen, the onset of convection is oscillatory. An increase in H is shown to stabilise the system. In Fig. 4 we display the effect of diffusivity ratio α. In this case an effect similar to that of H is observed.

Double diffusive viscoelastic porous convection

977

700 650 600 550 500 450 400 350 300 0

0.2

0.4

0.6

0.8

1

1.2

Fig. 6 Variation of critical Rayleigh number with strain retardation parameter λ2 for different values of stress relaxation parameter λ1

Fig. 7 Variation of critical Rayleigh number with strain retardation parameter λ2 for different values of porosity modified conductivity ratio γ

Figure 5 displays the effect of solute Rayleigh number Ra S on the marginal stability curves for fixed values of other governing parameters. We observe that when the solute Rayleigh number is small the convection sets in first via stationary instability. As the solute Rayleigh number increases beyond a certain value the convection sets in through oscillatory instability. It is also clear from this figure that the minimum of the Rayleigh number increases with the solute Rayleigh number indicating that the solute Rayleigh number delays the onset of convection in both stationary and oscillatory modes. The variation of the critical oscillatory Rayleigh number with strain retardation parameter λ2 for different parameters is shown in Figs. 6 and 7. We find that the convection sets in initially in the form of oscillatory instability, and as soon as the value of λ2 attains a critical value (say λ∗2 ), the convection ceases to be oscillatory, and steady convection occurs as the first bifurcation. The Rayleigh number then takes the value given by Eq. (26), which characterises the stationary convection and consequently becomes independent of λ2 . In Fig. 6, we display the effect of the relaxation parameter λ1 on the critical oscillatory Rayleigh number. We observe that for each value of λ1 the instability sets in initially in the form of overstable motions and the corresponding critical Rayleigh number increases with λ2 . As soon as λ2 > λ∗2 the convection ceases to be oscillatory and becomes steady. A similar behaviour is observed for the other physical parameters. The effect of the porosity modified conductivity ratio is shown in Fig. 7. We find that the effect of increasing γ is to

978

M. S. Malashetty et al.

Fig. 8 Variation of critical Rayleigh number with H for different values of a λ1 , b λ2

destabilise the system. Furthermore, it is interesting to note that when γ = 0.01, for the parameters chosen, oscillatory motion is not possible. The behaviour of the critical values of the Rayleigh number as functions of the inter-phase heat transfer coefficient H is depicted in Figs. 8, 9 and 10. In general, it is observed that for very small and large values of H the stability criterion is found to be independent of H . However, the effect of H on the stability of the system is significant for intermediate values of H . The physical reason for this is that when H → 0 there is almost no transfer of heat between the fluid and solid phases, and the properties of the solid phase have no significant influence on the onset criterion. When the fluid and solid phases have almost equal temperatures they may be treated a single phase. Between these two extremes H gives rise to a strong non-equilibrium effect. The effect of the viscoelastic parameter λ1 on the critical Rayleigh number is shown in Fig. 8a. We find that increasing the relaxation parameter causes the critical oscillatory Rayleigh number to decrease. Furthermore, it is interesting to note that all curves lie between the stationary and finite amplitude curves. The retardation parameter λ2 follows a similar (reciprocal) pattern of influence in Fig. 8b. Figure 9a indicates the effect of the solute Rayleigh number on the critical Rayleigh number. For small values of solute Rayleigh number (Ra S ≤ 10) the convection sets in via stationary instability. As the value of Ra S is increased further, the finite amplitude motions become significant and the convection sets in through them. Fairly large Ra S (≥100) dampens the finite amplitude motions and the overstable mode becomes the most dangerous mode. The critical Rayleigh number for stationary, oscillatory and finite amplitude convection is found to increase with the solute Rayleigh number, indicating that the presence of additional diffusing component stabilises the system towards the stationary, oscillatory and finite amplitude convection.

Double diffusive viscoelastic porous convection

Fig. 9 Variation of critical Rayleigh number with H for different values of a Ra S , b Le

Fig. 10 Variation of critical Rayleigh number with H for different values of γ

979

980

M. S. Malashetty et al.

Fig. 11 Variation of Nusselt number and Sherwood number with time for different values of a Ra S , b Le

Figure 9b displays the variation of the critical Rayleigh number with H for different values of the Lewis number. It is found that for Le = 2 the oscillatory motions occur prior to the stationary ones when H is very small and very large, whilst the trend reverses for intermediate values of H . However, the finite amplitude motions occur prior to both stationary and oscillatory motions for all values of the Lewis numbers. This figure also reveals that with increasing values of Le the critical Rayleigh number for stationary mode increases whilst that for oscillatory and finite amplitude convection decreases. The Lewis number has very small effect on the finite amplitude convection. Therefore, the Lewis number enhances the double diffusive convection in stationary mode whilst suppresses the oscillatory and finite amplitude convection. The variation of the critical Rayleigh number with H for different values of porosity modified conductivity ratio γ is depicted in Fig. 10. We observe from this figure that for small H the critical Rayleigh number is independent of γ and is close to that of the LTE case. This is due to the fact that for very small values of H there is no significant transfer of heat between the phases, and the onset criterion is not affected by the properties of the solid phase. On the other hand, for large values of H , although the stability criterion is also independent of H , the condition for the onset of convection is based on the mean properties of the medium and therefore the critical Rayleigh number is a function of γ . This figure also indicates that for moderate and large values of H the critical Rayleigh number for each of stationary, oscillatory and finite amplitude convection decreases with increasing values of γ . Therefore, the effect of the porosity modified conductivity ratio is to reduce the stabilising effect of the inter-phase heat transfer coefficient. It is important to note that for sufficiently large values of γ the critical Rayleigh number becomes independent of H . The effect of various parameters on the transient behaviour of the Nusselt and Sherwood numbers is shown in Figs. 11 and 12. We observe that the Nusselt and Sherwood numbers oscillate initially and settle to the stationary

Double diffusive viscoelastic porous convection

981

Fig. 12 Variation of Nusselt number and Sherwood number with time for different values of a H , b α

state (steady state value of 3) as time progresses and that the Nusselt number attains the steady state value earlier than the Sherwood number. This is in conformity with the result obtained for the steady state case. One of the shortcomings of the present non-linear models is that due to the severe truncation of the Fourier series the Nusselt number and the Sherwood number are restricted to the maximum value of 3. Better results can be obtained by including a higher number of terms in the Fourier series expansion of the variables. Unless stated otherwise, the parameters for Figs. 11 and 12 are H = 100, γ = 0.5, α = 0.2, λ1 = 0.7, λ2 = 0.5, Ra S = 1,000, Le = 2.5 and Pr D = 10, and Ra = 5 × RacF . Figure 11a shows that an increase of the solute Rayleigh number increases both heat and mass transfer. In Fig. 11b the effect of the Lewis number is found to increase both heat and mass transfer, with its effect being more pronounced on mass transfer, such that an increase in the Lewis number takes longer to attain the steady state value for the Sherwood number. Figure 12a shows the effect of the interface heat transfer coefficient. We find that its effect is more pronounced on heat transfer compared to the mass transfer. An increase in H decreases the heat transfer whilst it increases the mass transfer slightly. Figure 12b shows the effect of conductivity ratio α on heat mass transfer. We find that an increase in conductivity ratio decreases the heat transfer, whereas its effect on mass transfer is insignificant. Appendix A 1 = b1 + ωi2 b2 , 2 = b3 ωi6 + b4 ωi4 + b5 ωi2 + b6 ,

982

M. S. Malashetty et al.



1 + λ1 λ2 ωi2 LeRa S , + (δ 2 + γ H )2 + α 2 ωi2 δ 4 + Le2 ωi2 a 2 (1 + λ21 ωi2 )



αγ H 2 δ2 1 Le2 Ra S λ1 − λ2 b2 = 1 + 2 + − , a 2 1 + λ21 ωi2 Pr D (δ + γ H )2 + α 2 ωi2 δ 4 + Le2 ωi2   b3 = α 2 δ 2 Le2 λ1 λ1 (δ 2 + H ) + λ2 Pr D , b1 = δ 2 δ 2 +

H (δ 4 + δ 2 γ H + α 2 ωi2

b4 = λ21 (α 2 + Le2 )δ 8 + λ1 [λ1 H (α 2 + Le2 (2γ + 1)) + λ2 Pr D (α 2 + Le2 )]δ 6 + Le2 [α 2 + λ21 γ (γ + 1)H 2 + Pr D (α 2 λ2 + λ1 (2γ H λ2 − α 2 ))]δ 4 + Leα 2 [a 2 Pr D Ra S + Le(H + Pr D (1 − a 2 Ra S λ21 +

γ H2 (α + γ )λ1 λ2 + H (λ2 − λ1 )))]δ 2 − a 2 α 2 λ21 Le2 Pr D Ra S H, α2

b5 = λ21 δ 12 + λ1 (λ1 H (2γ + 1) + λ2 Pr D )δ 10 + [α 2 + Le2 + λ21 H 2 γ (γ + 1) + Pr D (α 2 λ2 + Le2 (λ2 − λ1 ) + λ1 (2γ H λ2 − α 2 ))]δ 8 + [Pr D (a 2 λ21 LeRa S + α 2 (1 + λ2 H + λ1 H (λ2 H γ (γ + α) − α 2 ) + Le2 (1 − a 2 λ21 Ra S + H (2γ + 1)(λ2 − λ1 ))) + H (α 2 + Le2 (2γ + 1))]δ 6 + H [H 2 Le2 γ (γ + 1) + Pr D Le(2a 2 λ21 Ra S γ + Le(2γ − a 2 Ra S (2γ + 1)λ21 + H 2 γ (γ + 1)(λ2 − λ1 )))]δ 4 + Pr D Le[a 2 Ra S (α 2 + H 2 γ (γ + α)λ21 ) + Le(H 2 γ (γ + α) − a 2 Ra S (α 2 + H 2 γ (γ + 1)λ21 ))]δ 2 − a 2 α 2 Le2 H Ra S Pr D , b6 = (1 + Pr D (λ2 − λ1 ))δ 12 + [Pr D + H (1 + 2γ )(1 + Pr D (λ2 − λ1 ))]δ 10   + 2γ H Pr D + H 2 γ (γ + 1)(1 + Pr D (λ2 − λ1 )) δ 8   + a 2 Pr D Ra S Le(1 − Le) + H 2 Pr D γ (γ + α) δ 6 + a 2 H Pr D Ra S Le[(2γ (1 − Le) − Le)δ 2 + H γ ((γ + α) − Le(γ + 1))]δ 2 , 4π 2 + γ H , 4π 2 + γ H + H   Le2 (δ 2 + γ H + H ) A2 = a 2 δ 4 A∗ + + a 4 Le (Ra S A∗ − LeRaT ) , δ2 + γ H     H A3 = δ 2 (δ 4 + a 2 LeRa S ) 1 + 2 − a 2 RaT , δ +γH A1 = a 4 δ 2 Le2 A∗ ;

A∗ =

B1 = a 4 Le4 ,

     Le2 (δ 2 + γ H + H ) H 2 B2 = 2a 2 Le2 δ 4 A∗ − LeRa + a − 1 , S δ2 + γ H 4π 2 + γ H + H B3 = a 4 [a 2 LeRa S (4π 2 + γ H )(δ 2 + γ H ) − δ 4 (4π 2 δ 2 (Le2 − 1) + H 2 (Le + γ (Le − 1))(Le + γ (Le + 1)) + H (Le2 + γ (Le2 − 1))(4π 2 + δ 2 ))]2 [(4π 2 + γ H + H )(δ 2 + γ H )2 ]−1 .

References 1. 2. 3. 4.

Ingham, D.B., Pop, I.: Transport Phenomena in Porous Media Volume III. Elsevier, Amsterdam (2005) Nield, D.A., Bejan, A.: Convection in Porous Media. 3rd edn. Springer, Berlin (2006) Vadasz, P.: Emerging Topics in Heat and Mass Transfer in Porous Media. Springer, Berlin (2008) Vafai, K.: Handbook of Porous Media. 2nd edn. Taylor & Francis, Boca Raton (2005)

Double diffusive viscoelastic porous convection

983

5. Malashetty, M.S., Heera, R.: Linear and non-linear double diffusive convection in a rotating porous layer using a thermal non-equilibrium model. Int. J. Non-Linear Mech. 43, 600–621 (2008) 6. Malashetty, M.S., Heera, R.: The onset of double diffusive convection in a sparsely packed porous layer using a thermal non-equilibrium model. Acta Mech. 204, 1–20 (2009) 7. Malashetty, M.S., Pop, I., Heera, R.: Linear and nonlinear double diffusive convection in a rotating sparsely packed porous layer using a thermal non-equilibrium model. Continuum Mech. Thermodyn. 21, 317–339 (2009) 8. Malashetty, M.S., Swamy, M.S., Heera, R.: Double diffusive convection in a porous layer using a thermal non-equilibrium model. Int. J. Thermal Sci. 47, 1131–1147 (2008) 9. Malashetty, M.S., Swamy, M.S., Heera, R.: The onset of convection in a binary viscoelastic fluid saturated porous layer. Z. Angew. Math. Mech. 89, 356–369 (2009) 10. Postelnicu, A.: The effect of a horizontal pressure gradient on the onset of a Darcy–Benard convection in a thermal nonequilibrium condition. Int. J. Heat Mass Transf. 53, 68–75 (2010) 11. Straughan, B.: Global nonlinear stability in porous convection with a thermal non-equilibrium model. Proc. Roy. Soc. Lond. A 462, 409–418 (2006) 12. Bories, S.A.: Natural convection in porous media. In: Bear, J., Corapcioglu, M.Y. (eds.) Advances in Transport Phenomena in Porous Media, pp. 77–141. Martinus Nijhoff, The Netherlands (1987) 13. Nield, D.A., Kuznetsov, A.V.: The interaction of thermal non-equilibrium and heterogeneous conductivity effects in forced convection in layered porous channels. Int. J. Heat Mass Transf. 44, 4369–4373 (2001) 14. Skartsis, L., Khomani, B., Kardos, J.L.: Polymeric flow through fibrous media. J. Rheology 36, 589–620 (1992) 15. Li, Z., Khayat, R.E.: Finite-amplitude Rayleigh–Benard convection and pattern selection for viscoelastic fluids. J. Fluid Mech. 529, 221–251 (2005) 16. Laroze, D., Martinez-Mardones, J., Bragard, J.: Thermal convection in a rotating binary viscoelastic liquid mixture. Eur. Phys. J. Special Top. 146, 291–300 (2007) 17. Laroze, D., Martinez-Mardones, J., Bragard, J., Perez-Garcia, C.: Realistic rotating convection in a DNA suspension. Physica A 385, 433–438 (2007) 18. Martinez-Mardones, J., Tiemann, R., Walgraef, D.: Rayleigh–Benard convection in a binary viscoelastic fluid. Physica A 283, 233–236 (2000) 19. Martinez-Mardones, J., Tiemann, R., Walgraef, D.: Amplitude equation for stationary convection in a binary viscoelastic fluid. Physica A 327, 29–33 (2003) 20. Bertola, V., Cafaro, E.: Thermal instability of viscoelastic fluids in horizontal porous layers as initial problems. Int. J. Heat Mass Transf. 49, 4003–4012 (2006) 21. Kim, M.C., Lee, S.B., Kim, S., Chung, J.B.: Thermal instability of viscoelastic fluids in porous media. Int. J. Heat Mass Transf. 46, 5065–5072 (2003) 22. Sheu, L.J., Chen, J.H., Chen, H.K., Tam, L.M., Chao, Y.C.: A unified system describing dynamics of chaotic convection. Chaos Solitons Fractals 41, 123–130 (2007) 23. Sheu, L.J., Tam, L.M., Chen, J.H., Chen, H.J., Lin, K., Kang, Y.: Chaotic convection of viscoelastic fluids in porous media. Chaos Solitons Fractals 37, 113–124 (2006) 24. Malashetty, M.S., Swamy, M.: The onset of convection in a viscoelastic liquid saturated anisotropic porous layer. Transp. Porous Media 67, 203–218 (2007) 25. Tan, W., Mausoka, T.: Stability analysis of a Maxwell fluid in a porous medium heated from below. Phys. Lett. A 360, 454–460 (2007) 26. Wang, S., Tan, W.: Stability analysis of double-diffusive convection of Maxwell fluid in a porous medium heated from below. Phys. Lett. A 372, 3046–3050 (2008) 27. Rajagopal, K.R., Bhatnagar, R.K.: Exact solutions for some simple flows of an Oldroyd-B fluid. Acta Mech. 113, 233–239 (1995) 28. Khuzhayorov, B., Auriault, J.L., Royer, P.: Derivation of macroscopic filtration law for transient linear viscoelastic fluid flow in porous media. Int. J. Engng. Sci. 38, 487–504 (2000) 29. Banu, N., Rees, D.A.S.: Onset of Darcy–Benard convection using a thermal non-equilibrium model. Int. J. Heat Mass Transf. 45, 2221–2228 (2002)