Arch Appl Mech (2011) 81:1697–1715 DOI 10.1007/s00419-011-0512-5

O R I G I NA L

I. S. Shivakumara · Jinho Lee · S. Suresh Kumar · N. Devaraju

Linear and nonlinear stability of double diffusive convection in a couple stress fluid–saturated porous layer

Received: 18 September 2010 / Accepted: 13 January 2011 / Published online: 4 February 2011 © Springer-Verlag 2011

Abstract The linear and nonlinear stability of double diffusive convection in a layer of couple stress fluid– saturated porous medium is theoretically investigated in this work. Applying the linear stability theory, the criterion for the onset of steady and oscillatory convection is obtained. Emphasizing the presence of couple stresses, it is shown that their effect is to delay the onset of convection and oscillatory convection always occurs at a lower value of the Rayleigh number at which steady convection sets in. The nonlinear stability analysis is carried out by constructing a system of nonlinear autonomous ordinary differential equations using a truncated representation of Fourier series method and also employing modified perturbation theory with the help of self-adjoint operator technique. The results obtained from these two methods are found to complement each other. Besides, heat and mass transport are calculated in terms of Nusselt numbers. In addition, the transient behavior of Nusselt numbers is analyzed by solving the nonlinear system of ordinary differential equations numerically using the Runge–Kutta–Gill method. Streamlines, isotherms, and isohalines are also displayed. Keywords Double diffusive convection · Couple stress fluid · Porous medium · Nusselt number

1 Introduction Double diffusive convection is a type of instability that occurs in a fluid that possesses two opposing densityaltering components with differing molecular diffusivities, such as heat and salt or any two solute concentrations. The marked difference between single and double diffusive systems lies in the verity that in double I. S. Shivakumara · J. Lee (B) School of Mechanical Engineering, Yonsei University, Seoul 120-749, South Korea, E-mail: [email protected] Present Address: I. S. Shivakumara · N. Devaraju UGC-Centre for Advanced studies in Fluid Mechanics, Department of Mathematics, Bangalore University, Bangalore 560 001, India E-mail: [email protected] N. Devaraju E-mail: [email protected] S. Suresh Kumar Department of Mathematics, Siddaganga Institute of Technology, Tumkur 572103, India E-mail: [email protected]

1698

I. S. Shivakumara et al.

diffusive systems convection can occur even when the system is hydrostatically stable if the diffusivities of the two diffusing fields are widely different. The study of double diffusive convection has received much attention over the years due to its numerous fundamental and industrial applications. Some examples of interest can be found in the fields of astrophysics, oceanography, chemical engineering, metallurgy and materials science, geophysics, and geology. Excellent reviews of these studies have been reported by Turner [1–3], Huppert and Turner [4], and Platten and Legros [5]. Double diffusive convection in a porous medium has also received due attention in the literature because of its importance in various fields such as high-quality crystal production, liquid gas storage, migration of moisture in fibrous insulation, transport of contaminants in saturated soil, solidification of molten alloys, and geothermally heated lakes and magmas. A review of this field is given by Trevisan and Bejan [6] and also documented in the books by Vafai [7,8] and Nield and Bejan [9]. It is well known that non-Newtonian fluids occur in a wide variety of industrial and natural flows, ranging from polymer processing and construction of oil wells, processing of food stuffs to lava and mud flows. There exist different types of non-Newtonian fluids. In particular, the theory of polar fluids has received wider attention in recent years because the traditional Newtonian fluids cannot precisely describe the characteristics of the fluid flow involved therein. These fluids deform and produce a spin field due to the microrotation of suspended particles forming micropolar fluid developed by Eringen [10]. The micropolar fluids take care of local effects arising from microstructure and as well as the intrinsic motions of microfluidics. The spin field due to microrotation of freely suspended particles set up an anti-symmetric stress, known as couple stress, and thus forming couple-stress fluid. Thus, couple-stress fluid, according to Eringen [10], is a particular case of micropolar fluid when microrotation balances with the natural vorticity of the fluid. The couple-stress fluid has distinct features, such as polar effects and whose microstructure is mechanically significant. Besides, couple stresses are found to appear in noticeable magnitude in fluids with large molecules. For such a special kind of non-Newtonian fluids, the constitutive equations are given by Stokes [11] and the theory proposed by him is the simplest one for microfluids, which allows polar effects such as the presence of couple stress, body couples, and nonsymmetric tensors. Thermal convective instability in either a layer of couple-stress fluid or a layer of couple stress fluid– saturated porous medium has drawn some attention of researchers in the recent past. The effect of rotation on a couple-stress fluid heated from below in a porous medium has been considered by Sharma et al. [12]. Sharma and Sharma [13] have investigated the effect of suspended particles on electrically conducting couple-stress fluid heated uniformly from below under the influence of uniform rotation and magnetic field. The effect of thermal/gravity modulation on the onset of Rayleigh–Benard convection in a layer of couple-stress fluid has been discussed by Malashetty and Basavaraja [14]. A few studies have also addressed the effect of an additional diffusing component on thermal convective instability in a layer of couple-stress fluid. An analytical study of linear and nonlinear double diffusive convection with Soret effect in couple-stress liquids has been considered by Malashetty et al. [15]. Gaikwad et al. [16] have investigated linear and nonlinear double diffusive convection by considering both Soret and Dufour effects in a layer of couple-stress liquid. Recently, Shivakumara [17] has discussed the effect of various nonuniform basic temperature gradients on the onset of convection in a layer of couple stress fluid–saturated porous layer, while Shivakumara et al. [18] have analyzed the effect of Coriolis force on both linear and weakly nonlinear thermal convection in a couple stress fluid–saturated rotating rigid porous layer. Double diffusive convection encountered in many practical problems in porous media often involves different types of dissolved substances of chemicals that are freely suspended in fluid–saturated porous media and they will be executing microrotation forming micropolar fluid. The presence of these suspended particles plays a major role in mixing processes. Although double diffusive convection in Newtonian fluid–saturated porous media has been studied extensively, the problem considering the above facts has not received due attention in the literature. The intent of the present paper is therefore to study double diffusive convection in non-Newtonian fluid–saturated porous media. In the present paper, Stokes [11] couple-stress model is being used in analyzing the problem. The linear stability theory is studied using normal mode technique, while the nonlinear stability analysis is performed using the minimal representation of double Fourier series as well as perturbation theory with the help of self-adjoint operator technique. The stability of finite amplitude steady solution is discussed and the condition for sub-critical instability is obtained. Besides, heat and mass transport are calculated in terms of Nusselt numbers. Also, the transient behavior of Nusselt numbers is analyzed by solving the nonlinear system of ordinary differential equations numerically using the Runge–Kutta–Gill method. Streamlines, isotherms, and isohalines are presented, and the effect of couple stress on the same is brought out.

Linear and nonlinear stability

1699

2 Mathematical formulation We consider a binary couple stress fluid–saturated horizontal porous layer of thickness d in which the density ρ depends on two different stratifying agencies, namely temperature and solute concentration. The lower surface is held at constant temperature T0 + T and constant solute concentration S0 + S, while the upper surface is held at T0 and S0 with T > 0 and S > 0. A Cartesian coordinate system (x, y, z) is chosen such that the origin is at the bottom of the porous layer and the z-axis is vertically upward in a gravitational field. Under the Boussinesq approximation, the equation of continuity and the equation of state are, respectively, ∇ · q = 0

(1)

ρ = ρ0 {1 − αt (T − T0 ) + αs (S − S0 )} .

(2)

and

Here, q = (u, v, w) is the velocity vector, T the temperature, S the solute concentration, αt the thermal expansion coefficient, αs the solute analog of αt , and ρ0 the reference density. The equation of motion of an incompressible couple-stress fluid in the absence of body couple and a porous medium is given by Stokes [11]     ∂ q f ρ0 (3) + q f · ∇ q f = −∇ p + ρ g + μ∇ 2 q f − μc ∇ 4 q f ∂t where q f is the velocity of the couple-stress fluid in the absence of porous medium, p the pressure, μ the dynamic viscosity, μc the material constant responsible for the couple-stress property known as the couplestress viscosity, and g the acceleration due to gravity. In the case of polar fluids, the action of one part of the body on its neighborhood cannot be represented by a force alone but rather by a force and couple. The last term on the right-hand side of Eq. (3) represents the effect of couple stresses in an incompressible fluid. When μc = 0, Eq. (3) reduces to the Navier–Stokes equation. The equation of motion describing the flow of couple-stress fluid through porous media can be viewed as follows. It is a known fact that the Darcy equation can be derived from the Navier–Stokes equation by statistical averages and simplifications of the complicated microscopic flow picture (Hubbert [19] and De Wiest [20]). Following the same procedure, the equation of motion for couple-stress fluid through porous media can be obtained from Eq. (3) and with the inclusion of inertial effect the equation reads as [17,18]   1 ∂ q 1 1 ρ0 q · ∇) q = −∇ P + ρ g − (μ − μc ∇ 2 ) q (4) + 2 ( ε ∂t ε k where q = εq f , ε the porosity, P the effective pressure, and k the permeability of the porous medium. Equation (4) is indeed the averaged equation applicable for couple-stress fluid flow through porous media. We note that the above equation reduces to the Darcy–Lapwood equation, when μc = 0. Under the circumstances, Eq. (4) can be regarded as modified Darcy–Lapwood equation for the flow of couple-stress fluid through porous media. The temperature and solute concentration equations are, respectively, given by ∂T + ( q · ∇) T = κ∇ 2 T ∂t

(5)

∂S + ( q · ∇) S = κs ∇ 2 S ∂t

(6)

Ah and ε

where Ah = (ρ0 c)m /(ρ0 c p ) f = [(1−ε)(ρ0 c)s +ε(ρ0 c p ) f ]/(ρ0 c p ) f the ratio of heat capacities, c the specific heat, c p the specific heat at constant pressure, κ the thermal diffusivity, and κs the solute analog of κ. The subscripts m, s, and f refer to the porous medium, solid, and fluid, respectively. The basic state is quiescent and is given by qb = (0, 0, 0) ,

P = Pb (z) , T = Tb (z) , S = Sb (z) , ρ = ρb (z) .

(7)

1700

I. S. Shivakumara et al.

Substituting Eq. (7) in Eqs. (4)–(6) and solving the equations using the appropriate boundary conditions, we obtain  z Tb (z) = T0 + T 1 − (8a) d  z (8b) Sb (z) = S0 + S 1 − d    2 z 1 Pb (z) = P0 − ρ0 g z+ (αt T − αs S) − 2z (8c) 2 d where P0 is the pressure at z = 0. It may be noted that the basic heat and solute concentration distributions are linear in z and the pressure distribution is of no consequence here as we are eliminating the same. On the basic state, we superimpose perturbations in the form q = q  , T = Tb + θ  , S = Tb + S  , P = Pb + P  , ρ = ρb + ρ 

(9)

where primes indicate perturbed quantities. Substituting Eq. (9) into the governing equations and using the basic state equations, we obtain the following stability equations: 

    1 ∂ q 1   1 ρ0 + 2 q · ∇ q = −∇ P  + αt gθ  kˆ − αs S  g kˆ − μ − μc ∇ 2 q ε ∂t ε k      ∂θ    + q · ∇ θ + q · ∇ Tb = κ∇ 2 θ  Ah ∂t    ∂ S   ε + q · ∇ S  + q · ∇ Sb = κs ∇ 2 S  . ∂t

(10) (11) (12)

We consider only two-dimensional motion and introduce stream function ψ through q = (u  , 0, w  ) =



∂ψ ∂ψ , 0, − ∂z ∂x

.

(13)

Eliminating the pressure term from the momentum equation by operating curl, and nondimensionalizing ψ by κ, x and z by d, time by d 2 Ah /κ, temperature by T , and solute concentration by S, we arrive at the following equations (after ignoring the primes for simplicity): 

    ∂ ∂θ ∂S 1 2 J ψ, ∇ 2 ψ + 1 − c ∇ + RsD + ∇ 2 ψ = −R D Pr D M ∂t ∂x ∂x Pr D  ∂ ∂ψ − ∇2 θ = − + J (ψ, θ ) ∂t ∂x  1 ∂ ∂ψ − τ ∇2 S = − + J (ψ, S) M ∂t ∂x 1

(14) (15) (16)

∂ f ∂g where J ( f, g) = ∂∂ xf ∂g ∂z − ∂z ∂ x is the Jacobian, R D = αt gT dk/νκ is the Darcy–Rayleigh number, RsD = αs gSdk/νκ is the solute Darcy–Rayleigh number, τ = κs /κ is the ratio of diffusivities, Pr D = νε2 d 2 /κk is the modified Darcy–Prandtl number, M = Ah /ε is a nondimensional group, c = μc /μd 2 is the couple-stress parameter, and ν = μ/ρ0 is the kinematic viscosity. The boundaries are assumed to be impermeable with vanishing couple stress and perfect conductors of both heat and solute concentration. Accordingly, the boundary conditions are:

∂ 2ψ = 0 at z = 0, 1. ∂z 2 θ = 0 = S at z = 0, 1.

ψ=

(17a) (17b)

Linear and nonlinear stability

1701

3 Linear stability theory In this section, we discuss the linear stability analysis that is very useful in the local nonlinear stability analysis discussed in the next section. In the linear stability theory, we neglect the Jacobians in Eqs. (14)–(16) and eliminate θ and S from Eq. (14) to obtain an equation for ψ in the form      ∂ ∂ 1 ∂ 1 2 2 2 + 1 − c ∇ −∇ − τ ∇ ∇ 2ψ Pr D M ∂t ∂t M ∂t  2  2 1 ∂ ∂ 2 ∂ ψ 2 ∂ ψ − τ∇ −∇ − RsD . (18) = RD M ∂t ∂x2 ∂t ∂x2 Assuming the solution for ψ in the form ψ = eσ t sin (αx) sin (π z)

(19)

where σ is the growth rate which is a complex quantity in general and α is the horizontal wave number. Substituting Eq. (19) in Eq. (18), we get an expression for the Darcy–Rayleigh number in the form      σ + δ2 M σ δ2 2 2  RsD + 2 RD =  + δ η σ + δ (20) α M Pr D σ + Mτ δ 2 where δ 2 = π 2 + α 2 is the total wave number and η = 1/δ 2 + c is a representative of the viscosity of the fluid. The onset of instability can be examined by setting the real part of σ equal to zero so that σ = iω. Substituting σ = iω in Eq. (20) and clearing the complex quantities from the denominator, it can be written in the form    M ω2 + Mτ δ 4 δ2 ω2  RsD + 2 δ 4 η − RD =  2 (21) + iωδ 2 N α M Pr D ω + M 2 τ 2 δ4 where δ2 M (Mτ − 1)  RsD + 2 α ω2 + M 2 τ 2 δ 4



N=

1 + δ2η . M Pr D

(22)

Since R D is a physical quantity, Eq. (21) implies either ω = 0 or N = 0, and accordingly, we can obtain the condition for the occurrence of stationary and oscillatory convection. 3.1 Stationary convection (ω = 0) (S)

Substituting ω = 0 in Eq. (21) gives the Darcy–Rayleigh number R D for the onset of stationary convection in the form R sD =

RsD δ6η + 2. τ α

We note that R sD attains its minimum value at α 2 = αc2 , where      − π 2 c + 1 + π 2 c + 1 9π 2 c + 1 . αc2 = 4 c

(23)

(24)

From Eq. (24), it is observed that the critical wave number αc is independent of an additional diffusing component (i.e., RsD ) but depends on the couple-stress parameter c . The critical wave number calculated from Eq. (24) for different values of c is tabulated in Table 1. We note that an increase in c is to decrease the critical wave number and hence its effect is to increase the size of convection cells. When c = 0 (i.e., in the absence of couple stresses), Eq. (23) reduces to R sD =

RsD δ4 + 2 τ α

(25)

1702

I. S. Shivakumara et al.

Table 1 Values of αc2 for different values of c

c

0

0.5

1.0

2.0

2.5

3.0

αc2

9.8722

5.241

5.094

5.016

5.000

4.989

which coincides with the result of Nield and Bejan [9]. Equation (25) attains its minimum value at αc = π and the minimum value is s RDmin =

RsD + 4π 2 . τ

(26)

s If RsD = 0 in Eq. (26), then we get back the classical value RDmin = 4π 2 for the case of Newtonian fluid saturating a porous medium.

3.2 Oscillatory convection (ω  = 0, N = 0) It is a well established fact that the presence of additional constraints like rotation and magnetic field as well as solute concentration gradient supports oscillatory convection. The oscillatory onset corresponds to N = 0 (ω  = 0), and this condition gives the dispersion relation ω2 =

M (1 − Mτ ) α 2 RsD − M 2 δ 4 τ 2 . δ 2 (η + 1/M Pr D )

(27)

Eliminating ω2 from Eq. (21) using Eq. (27) and noting N = 0, we obtain an expression for oscillatory Rayleigh number in the form   M 2 RsD δ 6 (Mτ + 1) R oD = (28) + (Pr D η + τ ) . Pr D α 2 (M Pr D η + 1) We note that R oD attains its minimum value at α 2 = αc2 , where  5  4  3  2   a5 δc2 + a4 δc2 + a3 δc2 + a2 δc2 + a1 δc2 + a0 = 0.

(29)

Here,   2 τ 1 a 5 = 2 c +

c + Pr D M Pr D       1 τ 1 1 a4 = 4 − 3π 2 c +

c + + c +

c + M Pr D Pr D M Pr D M Pr D       τ 1 1 1 2 − 6π 2 c + + 2 − 2π 2 c +

c + a 3 = c + Pr D M Pr D M Pr D M Pr D     RsD (Mτ − 1) τ 1 − 3π 2 c + + 1 − 4π 2 c + a2 = M Pr D (Mτ + 1) Pr D M Pr D   RsD (Mτ − 1) +1 a1 = −2π 2 Pr D (Mτ + 1) RsD (Mτ − 1) π 4 and δc2 = αc2 + π 2 a0 = Pr D (Mτ + 1) From Eq. (29), it is seen that αc2 depends not only on c but also on RsD , Pr D , τ , and M. As Pr D → ∞, which is usually the case, Eqs. (28) and (29), respectively, reduce to R oD = M RsD +

δ 6 η (Mτ + 1) α2

(30)

Linear and nonlinear stability

and

1703

 2  2  2   

c δc2 + 1 2 c δc2 + 1 − 3π 2 c δc2 − 2π 2 = 0.

(31)

2

  2 Since c δc2 + 1  = 0, it follows that  2    2 c δc2 + 1 − 3π 2 c δc2 − 2π 2 = 0

(32)

and on simplification we find that this equation is the same as Eq. (24). That is, when Pr D → ∞, the critical wave number remains the same for the onset of both steady and oscillatory convection. However, to find the critical wave number for arbitrary values of Pr D , we have to solve Eq. (29). Since ω2 > 0, which is the requirement for the occurrence of oscillatory convection, the necessary conditions for the occurrence of oscillatory convection are obtained from Eq. (27) and are given by  1 δ 6 τ 2 ηM 1 τ< 1+ and RsD > 2 . (33) M α (1 − Mτ ) M Pr D η From the above equation, it is evident that the oscillatory convection is possible only if the solute Darcy– Rayleigh number RsD exceeds a threshold value that in turn depends on the couple-stress parameter as well. Using Eqs. (23) and (28), it is interesting to note that Eq. (27) can be expressed as ω2 =

 s  M Pr D τ α 2 R D − R oD . 2 δ (M Pr D η + 1 + Mτ )

(34)

Hence, if oscillatory convection is possible, it always occurs at a lower value of R sD at which stationary convection occurs. 4 Nonlinear stability analysis Although the linear stability analysis gives the criterion for the onset of convection, it fails to provide information on the possibility of occurring sub-critical instability and also the rate of heat and mass transfer. Hence, in this section, we consider the weakly nonlinear stability analysis to obtain the additional information. In the immediate vicinity of the linear stability boundary, we could develop a weakly nonlinear analysis. We perform the weakly nonlinear stability analysis using the following methods: (i) Truncated representation of double Fourier series and (ii) Modified perturbation theory. 4.1 Truncated representation of double Fourier series In this method, the dependent variables are expanded in double Fourier series and construct a set of coupled nonlinear ordinary differential equations with solutions that mimic the essential features of solution to the full partial differential Eqs. (14)–(16). Moreover, the representation is just that which reproduces results obtained by modified perturbation theory to second order in amplitude, for the full problem, and qualitatively accurate for larger amplitudes. The weak nonlinear analysis has been carried out by considering only one term for stream function (ψ) and two terms for temperature (θ ) and solute concentration (S) distributions. Following Rudraiah et al. [21] and Shivakumara and Sumithra [22] studies, ψ, θ , and S are represented in the following form: √ 2 2δ ψ= (35) A(t ∗ )sin (αx) sin (π z) α √ C(t ∗ ) 2 2 B(t ∗ )cos (αx) sin (π z) − sin (2π z) (36) θ= δ π √ 2 2 E(t ∗ ) S= D(t ∗ )cos (αx) sin (π z) − sin (2π z) (37) δ π

1704

I. S. Shivakumara et al.

where t ∗ = δ 2 t and the amplitudes A(t ∗ ) − E(t ∗ ) are to be determined from the dynamics of the system. Substituting Eq. (35a, b, c) into Eqs. (14)–(16) and equating the coefficients of like terms lead to the following nonlinear system of autonomous nonlinear ordinary differential equations:   α 2 RsD α2 R D ˙ A = M Pr D −η A − B+ D (38a) δ6 δ6 B˙ = A (C − 1) − B (38b) ˙ C = β (AB + C) (38c) ˙ D = M [−τ D + A (E − 1)] (38d) E˙ = −Mβ [AD + τ E] (38e) where β = 4π 2 /δ 2 and the dot above the quantities denotes the differentiation with respect to time t ∗ . The above nonlinear system of ordinary differential equations is not amenable to analytical treatment and we have to solve it numerically. We note that the above system of equations is uniformly bounded in time and possesses many properties of the full problem. First, the divergence of the flow in a five-dimensional phase space, ∂ A˙ ∂ B˙ ∂ C˙ ∂ D˙ ∂ E˙ + + + + = − (M Pr D η+1 + β + Mτ + β Mτ ) , ∂A ∂B ∂C ∂D ∂E

(39)

is always negative, and therefore, the system is bounded and dissipative. As a result, the solutions are attracted to a set of measure zero in the phase space and this may be a fixed point, a limit cycle, or a strange attractor. Besides, from Eq. (38a–e), it can also be concluded 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     V t ∗ = V (0) ex p − (M Pr D η+1 + β + Mτ + β Mτ ) t ∗ . (40) Equation (40) indicates that the volume decreases with time. Moreover, increase in the value of couple-stress parameter and the Prandtl number is to enhance the dissipation. Further, we observe that the equations have an important symmetry, for they are invariant under the transformation (A, B, C, D, E) → (−A, − B, C, − D, E) .

(41)

Taking a cue from the above observations, we now look into the possibility of an analytical solution. For the steady case, Eq. (38a–e) can be solved in a closed form and such a study is useful because it predicts the possibility of the occurrence of sub-critical instability. Setting the left-hand sides of Eq. (38a–e) equal to zero, we arrive at ηδ 6 A + R D B − RsD D = 0 α2 A (1 − C) + B = 0 AB + C = 0 A (1 − E) + τ E = 0 AD + τ E = 0

(42a) (42b) (42c) (42d) (42e)

On solving for the amplitudes in terms of A, we get B=−

A2 A , C = , A2 + 1 A2 + 1

D=−

τA , 2 A + τ2

E=

A2 . A2 + τ 2

(43)

Substituting for B and D in the first equation of (42a), we get   2      RsD η˜ A2 + η˜ 1 + τ 2 − R D + RsD τ A2 + τ 2 η˜ − R D + =0 τ where η˜ = ηδ 6 /α 2 .

(44)

Linear and nonlinear stability

1705

The condition that should be satisfied is that the discriminant of Eq. (44) is non-negative. Hence, the limiting condition for real A2 is    2  RsD 2 2 η˜ 1 + τ − R D + RsD τ − 4ητ ˜ = 0. η˜ − R D + τ

(45)

f

The solution of this equation for R D is denoted by R D , called the finite amplitude Rayleigh number, which characterizes the onset of finite amplitude steady motions. The finite amplitude Rayleigh number is given by f RD

=





  2 2 RsD τ + η˜ 1 − τ .

(46)

f

The second solution for R D leads to negative values for R D and is therefore omitted. Equation (46) gives f f meaningful results only when RsD > 0 and τ < 1. We note that R D attains its minimum value, RDmin at α 2 = αc2 , where

αc2 =

     − π 2 c + 1 + π 2 c + 1 9π 2 c + 1 4 c

(47)

and this is the same as Eq. (27). For a single component fluid (i.e., RsD =0), from Eq. (46) we find that f

R D = ηδ 6 /α 2 = R sD

(48)

and note that sub-critical motions are not possible. Since we are dealing with weakly nonlinear stability analysis, the amplitudes are assumed to be small. Accordingly, we can expand R D in powers of A2 (A2 << 1) in the form R D = R sD + R sD2 A2 + · · · .

(49)

Substituting Eq. (49) into Eq. (45), and collecting the coefficients of different powers of A2 , we obtain the following results. At zeroth order in A2 : R sD =

RsD δ6η + 2. τ α

(50)

This is just the linear stability analysis result obtained in Sect. 3.1. At the first order in (A2 ): R sD2

  RsD τ 2 − 1 ηδ 6 = 2 + . α τ3

(51)

This is the first nontrivial finite amplitude Rayleigh number and note that R sD2 may be either positive or negative. The finite amplitude solution is said to be stable (i.e., supercritical) if R sD2 > 0 and unstable (i.e., sub-critical) if R sD2 < 0 when ω2 < 0. From Eq. (51), it is clear that the finite amplitude solution is stable if   τ > 1 or RsD < ηδ 6 τ 3 / 1 − τ 2 α 2 and τ < 1.

(52)

In the absence of an additional diffusing component (i.e., RsD = 0), we find that R sD2 = R sD , and hence, sub-critical instability is not possible.

1706

I. S. Shivakumara et al.

4.2 Modified perturbation theory In this section, modified perturbation theory is employed to the full partial differential Eqs. (14)–(16) to investigate the stability of steady finite amplitude solutions in the neighborhood of steady onset at R sD . Such a study also helps in knowing the accuracy of the truncated model used in the previous section to study the nonlinear stability of the problem. For steady case, Eqs. (14)–(16) become     ∂θ ∂S 1 J ψ, ∇ 2 ψ + RsD + 1 − c ∇ 2 ∇ 2 ψ = −R D ∂x ∂x Pr D ∂ψ ∇ 2θ = − J (ψ, θ ) ∂x ∂ψ − J (ψ, S) . τ ∇2 S = ∂x

(53) (54) (55)

We expand the dependent variables ψ, θ , and S and one of the Rayleigh numbers R D in terms of a small parameter ε p , such that ψ = ε p ψ1 + ε2p ψ2 + · · ·

(56a)

θ = ε p θ1 + ε2p θ2 + · · ·

(56b)

ε p S1 + ε2p S2

(56c)

S=

+ ···

and R D = R sD + ε p R sD1 + ε2p R sD2 + . . . . At each stage in the expansion, we define a column vector ⎡ ⎤ ψn n = ⎣ θn ⎦ . Sn

(57)

(58)

Substituting Eqs. (56a–c) and (57) into (58)–(55) then at leading order in ε p , the equations are linear, homogeneous and can be written in the operator form 1 = 0 L

(59)

where L is a self-adjoint differential operator and is given by ⎡ ∂ s ∂ L1 rD ∂ x rsD ∂ x ⎢ s ∂ s ∇2 0 L = ⎣ −r D ∂ x r D rsD ∂∂x

and

0

−rsD

⎤ ⎥ ⎦

(60)

τ ∇2

⎡

⎤ ψ1  1 = ⎣ θ1 ⎦  S1

(61)

with L1 =

2  α2  α2 2 s s α ∇ = R and r = R . 1 −

, r c sD sD D D δ6 δ6 δ6

This is just a two-dimensional form of the eigenvalue problem and has been already discussed in Sect. 3.1 and for which the eigenvalue is s = rsD /τ + η rD

(62a)

Linear and nonlinear stability

1707

and the eigenfunction is √ ⎡ sin (αx) sin (π z) ⎤ 2 2δ ⎣ α 1 = cos (αx) sin (π z) ⎦ .  δ2 α α cos (αx) sin (π z) δ2

(62b)

The normalization chosen above is for subsequent convenience. The operator L is self-adjoint, and hence, we have the identity      1T Ln =   nT L1 = 0 

(63)

for all values of n. Here, the angular brackets denote averages over the domain 0 ≤ x ≤ 2π/α and 0 ≤ z ≤ 1. To the second order in ε p , the equations are linear, inhomogeneous and are given by ⎡ ⎢ ⎢  L 2 = ⎢ ⎢ ⎣

 ⎤ J ψ1 , ∇ 2 ψ 1 ⎥ ⎥ ⎥. s −r D J (ψ1 , θ1 ) ⎥ ⎦ rsD J (ψ1 , S 1 )

s ∂θ1 + −r D1 ∂x

α2 Pr D δ 6

(64)

  From Eq. (62), it follows that J ψ1 , ∇ 2 ψ1 = 0, and the solvability condition on this inhomogeneous equation s = 0. If we impose the orthogonality is obtained by applying Eq. (63) to Eq. (64). Then, it follows that r D1 condition   1T n = 0, (n  = 1) , (65) then the solution of Eq. (64) is found to be ⎡

⎤ 1 1 ⎥ 2 = ⎢  ⎣ − π sin (2π z) ⎦ . − τ13 sin (2π z)

(66)

Now, at the third order in ε p , we have √ ⎤ 2 2α s r sin sin z) (αx) (π ⎥ ⎢ √ sδ D2 ⎢ 4 2δr cos (αx) sin (π z) cos (2π z) ⎥ . ⎦ ⎣ √ D − 4 τ 32δ rsD cos (αx) sin (π z) cos (2π z)

⎡ L ψ 3 =

(67)

As before the solvability condition yields, from Eq. (63) and Eq. (67), s s r D2 = rD −

rsD . τ3

(68)

That is, R sD2

  RsD τ 2 − 1 ηδ 6 = 2 + . α τ3

(69)

We note that Eq. (69) is the same as Eq. (51). Thus, the finite amplitude solution obtained from the truncated system in the previous section is identical to the one arrived at from the full two-dimensional problem.

1708

I. S. Shivakumara et al.

5 Convective heat and mass transport The vigor of double diffusive convection can be measured in terms of Nusselt numbers N u (thermal) and N u s (solute). In the steady state, the vertical heat (or solute concentration) flux is independent of the vertical coordinate, and it can be evaluated as HT = −κ <

∂T >z=0 ∂z

(70)

where the angular brackets correspond to a horizontal average. Since ∂ T /∂z is composed of the constant gradient plus the change in the mean field, the Nusselt number is calculated as Nu =

HT d = 1 + 2C. κT

(71)

Substituting for the amplitude C in terms of A, we obtain Nu = 1 +

2 A2 . A2 + 1

(72)

where A2 can be evaluated from Eq. (44). Similarly, the mass transport is calculated in terms of solute Nusselt number (Sherwood number) and it is given by N u s = 1 + 2E.

(73)

Substituting for the amplitude E in terms of A, we obtain N us = 1 + At R D = R sD ,

A2 = 0 or A2 =

2 A2 . A2 + τ 2

  RsD 1 − τ 2 ητ ˜

(74)

− τ 2.

(75)

We note that A2 = 0 corresponds to the conduction state and in this case both N u = 1 = N us .

(76)

f

At R D = R D ,

 A = 2

RsD (1 − τ 2 ) − τ 2. η˜

(77)

As RsD → ∞, A → ∞ and we note that both N u and N u s will approximate to 3. 6 Results and discussion Both linear and nonlinear stability of double diffusive convection in a layer of couple stress fluid–saturated porous medium is investigated. The linear theory is used to obtain the criterion for the onset of stationary and oscillatory convection. The weakly nonlinear stability analysis is carried out by constructing a fifth-order system of nonlinear ordinary differential equations which possesses both periodic and steady solutions and allows sub-critical instability. In addition, modified perturbation theory is used to solve full two-dimensional nonlinear equations to obtain the finite amplitude solution. Besides, heat and mass transfer are quantified in terms of Nusselt numbers. The behavior of the system is analyzed through Darcy–Rayleigh number as a function of solute Darcy–Rayleigh number RsD , the diffusivity ratio τ , the couple-stress parameter c , the nondimensional group M, the Darcy–Prandtl number Pr D , and the wave number α.

Linear and nonlinear stability

1709

_ _ _ _ Direct mode

Oscillatory mode

4

10

RD

PrD = 10,100

3

10

0.2

1

10

a

Fig. 1 Neutral stability curves for different values of Pr D with RsD = 10, τ = 0.01, M = 1.25, and c = 2

Based on the expression for Rayleigh numbers at which the steady and oscillatory convection occur, the neutral stability curves are plotted in Fig. 1 for the fixed values of RsD = 10, τ = 0.01, M = 1.25 and c = 2 with Pr D = 10 and 100. From this figure, we note that the neutral curves are connected in a topological sense and the oscillatory mode is the preferred mode of instability. Also, the oscillatory neutral curves coalesce for the values of Prandtl number considered, indicating insignificance of the fluid conductivity on the onset of oscillatory convection. This characteristic is in conformity with the classic results observed by Baines and Gill [23] in the study of double diffusive instabilities involving Newtonian fluids. It may also be noted that the critical wave number is the same for both steady and oscillatory convection as noticed earlier for large values of Pr D . Figure 2 summarizes the stability regions in the (R D , RsD )-plane, for different values of c with Pr D = 10, τ = 0.32, and M = 1.25. In the figure, the portion of each stability boundary lying to the left of the discontinuity in slope corresponds to steady onset (ω = 0) while to the right, the onset is of oscillatory type (ω2 > 0). From the figure, it is evident that an increase in the value of c is to increase the region of stability and also the point of discontinuity in slope. That is, the presence of couple stresses is to delay the onset of 1.5 1

2000

0.5

RD

Λc = 0

1000

-1000

500

-500 0

-1000

-500

00

1000

RsD 500

-1000

-2000

Fig. 2 Stability boundaries for different values of c with Pr D = 10, τ = 0.32, M = 1.25

1000

1710

I. S. Shivakumara et al.

400

300

=200

2

200

100 100

50 20

0 0.1

Fig. 3 Variation of

ωc2

1

10

as a function of c for different values of RsD with Pr D = 10, τ = 0.01, and M = 1.25 400

RDo min

m RDmin

R Df m in

300

RDmin

Λc = 0.3 R Df m in

200

0.2

R Df m in

0.1 100 -1 10

0

10

R sD

1

10

2

10

f

m o Fig. 4 Variation of RDmin , RDmin and RDmin as a function of RsD for different values c with Pr D = 10, τ = 0.01, M = 1.25

steady and oscillatory convection. At ∗ = RsD = RsD

Mτ 2 α 2 δ 6 η (1 − Mτ )

(78)

the slope of the (R D , RsD )-plot changes as does the preferred mode of instability. From Eq. (78), it can be seen ∗ increases with an increase in the value of couple-stress parameter. The above condition is that the value of RsD obtained for the case of Pr D → ∞ since the variation in Pr D has negligible effect on the onset of oscillatory convection. The critical wave numbers obtained from Eq. (29) for different values of c and RsD with Pr D = 10, τ = 0.01, and M = 1.25 are substituted back in Eq. (27) to calculate the critical frequency of oscillations. Figure 3 depicts the variation of ωc2 as a function of c for different values of RsD . From the figure, we note that increase in c is to decrease monotonically the frequency of oscillations and hence its effect is to suppress the oscillations. Whereas increase in RsD is to increase the value of ωc2 . To understand the possibility of occurring finite amplitude sub-critical instability, the values of f f s o RDmin , RDmin , and RDmin (i.e., minimum values of R sD , R oD and R D obtained with respect to the wave number)

Linear and nonlinear stability

1711

2.5

3.0

(a) .........

_ _ _ _ Nu s Nu

Nus

Λc=0.1

Nu

(b) 0.7

2.5

2.0

0.2

Nu,

τ = 0.81 2.0

Nus

0.3

Nu,

Nus

1.5

1.5

1.0 100

200

RD

300

1.0 100

400

200

RD

300

400

Nu 2.5

-----

Nus (c)

2.0

Nu,

Nus R

sD

= 10

1.5

60 30 1.0 100

200

RD

300

400

Fig. 5 Variation of Nusselt numbers for Pr D = 10, M = 1.25 different values of (a) c with RsD = 10, τ = 0.81 (b) τ with RsD = 10, c = 0.1 (c) RsD with c = 0.1, τ = 0.81

are plotted as functions of RsD for different values of c with Pr D =10, τ = 0.01, M = 1.25 and are shown in Fig. 4. From this figure, it is seen that an increase in c is to increase the value of R D and to reinforce the stability on the system. It is also evident that, although the onset of convection is of oscillatory type according to the linear theory, the finite amplitude steady motions (i.e., sub-critical motions) exist for values o of the Rayleigh number even smaller than RDmin . Thus, once the system becomes unstable, whether it is to o o infinitesimal perturbations at R D = RDmin , or to finite amplitude perturbations at R D < RDmin , the developed state of convection will be steady. f As RsD → ∞, R sD , R oD , and R D asymptotically become 1 RsD τ (Pr D η + τ ) R oD → RsD (M Pr D η + 1) R sD →

f

R D → τ RsD .

(79a) (79b) (79c)

1712

I. S. Shivakumara et al.

4

3.5

3

Nu 2.5

2

0

2

t*

4

6

8

10

12

14

16

5 4.5 4

3.5

Nus

3 2.5

2 1.5

0

2.5

5

t*

7.5

10

12.5

15

17.5

20

Fig. 6 Variation of a N u and b N u s with time for M = 1.25, Pr D = 10, R D = 4, 400, RsD = 10, τ = 0.32 c = 0.1 (dotted line), and c = 0.5 (dashed lines)

In the study of double diffusive convection, the determination of heat and mass transport across the porous layer plays a very important role. The onset of convection, as the thermal Darcy–Rayleigh number is increased, is more rapidly detected by its effect on the heat and mass transfer. The quantity of heat and mass transfer across the layer are given by N u and N u s , respectively. Figure 5a–c indicate the effect of R D ψ on the thermal (N u) and solute (N u s ) Nusselt numbers ψ for different values of couple-stress parameter c , the diffusivity ratio τ and the solute Darcy–Rayleigh number RsD , respectively. In each of these cases, we observe that as R D increases, the heat and mass transfer increase sharply with the solute Nusselt number is above the thermal Nusselt number and finally reach asymptotic values with further increase in R D . We also note that the effect of increasing c is to decrease the values of N u and N u s (see Fig. 5a) because increase in c is to delay the onset of convection. It is interesting to note that both N u and N u s decrease with an increase in RsD and decrease in τ only up to a certain value of R D and exceeding which the trend reverses (see Figure. 5b, c). This may be attributed to the fact that increase in the thermal Rayleigh number is to mix the solute and redistribute it so that the interior layers of the fluid are more neutrally stratified. Due to this, the inhibiting effect of solute gradient is greatly reduced and hence the fluid convects more heat and mass when RsD is increased.

Linear and nonlinear stability

1713

4.5

4

3.5

Nu

3

2.5

2

1.5

2

0

4

t*

6

8

10

12

14

5 4.5 4

3.5

Nus 3 2.5 2

1.5

0

2

4

t*

6

8

10

12

14

Fig. 7 Variation of a N u and b N u s with time for M = 1.25, Pr D = 10, R D = 4, 400, τ = 0.32, c = 0.1, RsD = 10 (dottedlines), and RsD = 100 (dashed lines)

The transient behavior of thermal and solute Nusselt numbers has been considered by solving numerically an autonomous nonlinear system of ordinary differential Eq. (42a–e) using the Runge-Kutta-Gill method with appropriate initial conditions. The results obtained are depicted in Figs. 6a, b and 7a, b for different values of

c ψ and RsD , respectively, for M = 1.25, Pr D = 10, R D = 4, 400, and τ = 0.32. Figure 6a, b, respectively, show the variation of N u and N u s with time for two values of c = 0.1 and 0.5, when RsD = 10. From the Figures, it is noted that both N u and N u s oscillate initially with time t ∗ and finally reaches a steady-state value with further increase in the value of t ∗ . In addition, increase in the value of couple-stress parameter is to decrease the amplitude of the oscillations of heat flux, but to increase the amplitude of the oscillations of mass flux. The variation of N u and N u s with time for two values of RsD = 10 and 100, when c = 0.1, is presented in Fig. 7a, b, respectively, and we observe a similar effect as noticed in the previous case by increase in the value of c .

1714

I. S. Shivakumara et al.

1

(a) 0.9 0.8 0.7

Z

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

1.2

X 1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

(c)

(b)

Z

Z

1

0

0.2

0.4

0.6

X

0.8

1

1.2

0

0

0.2

0.4

0.6

0.8

1

1.2

X

Fig. 8 Contour plots of a streamlines, b isotherms and c isohalines for RsD = 10, τ = 0.81, R D = 300, c = 0.1 (dotted lines) and 0.3 (dashed lines)

Figure 8a–c show the streamlines, isotherms, and isohalines, respectively, for two values of c = 0.1 and 0.3 when RsD = 10, R D = 300, and τ = 0.81. We have confined our attention to the rectangular region 0 ≤ x ≤ π/α and 0 ≤ z ≤ 1 since the cell-scale α is fixed for a given choice of physical parameters. The streamlines showed in Fig. 8a reveal that an increase in the value of c is to increase the size of convection cells. The rectangular streamlines pattern gets deformed into circular ones as the value of ψ is increased. The contours of isotherms (Fig. 8a) and isohalines (Fig. 8b) are more or less vertical in the centre of the cell, while they are horizontal near the lower and upper boundaries and flat at x = 0 and π/α. Increase in the value of c is to diminish the anvil-shaped plumes. 7 Conclusions The important results of linear and nonlinear double diffusive convection in couple stress fluid–saturated porous medium may be summarized as follows. (1) The critical wave number is independent of solute Darcy–Rayleigh number and remains same for the onset of steady, oscillatory when Pr D → ∞ and finite amplitude steady convection, but alters with couple-stress parameter c . Increase in c is to decrease the critical wave number.

Linear and nonlinear stability

1715

(2) The oscillatory neutral curves coalesce for different values of Prandtl number, indicating insignificance of the fluid conductivity on the onset of oscillatory convection. (3) Increase in the value of c is to delay the onset of convection and also to suppress the frequency of oscillations. Whereas increase in RsD is to increase the frequency of oscillations. (4) The finite amplitude critical Rayleigh number is less than the critical Rayleigh number for stationary and oscillatory convection indicating the existence of sub-critical instability. (5) Heat and mass transfer decreases with an increase in the value of couple-stress parameter c , but the diffusivity ratio τ, and the solute Darcy–Rayleigh number RsD shows a dual effect on the same. (6) Increase in c and RsD is to decrease the amplitude of oscillations of heat flux and to increase the amplitude of the oscillations of mass flux. Acknowledgments One of the authors (ISS) wishes to thank the Brain Korea 21 (BK21) Program of the School of Mechanical Engineering, Yonsei University, Seoul, Korea, for inviting him as a visiting Professor and also the Bangalore University for sanctioning sabbatical leave. One of the authors (SSK) would like to thank the Director, Principal, and Management of Siddaganga Institute of Technology, Tumkur, for their encouragement. We thank the reviewer for useful suggestions.

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23.

Turner, J.S.: Buoyancy Effects in Fluids. Cambridge University Press, London (1973) Turner, J.S.: Double diffusive phenomena. Annu. Rev. Fluid Mech. 6, 37 (1974) Turner, J.S.: Multicomponent convection. Annu. Rev. Fluid Mech. 17, 11 (1985) Huppert, H.E., Turner, J.S.: Double diffusive convection. J. Fluid Mech. 106, 299 (1981) Platten, J., Legros, J.C.: Convection in Liquids. Springer, Berlin (1984) Trevisan, O.V., Bejan, A.: Combined heat and mass transfer by natural convection in a porous medium 20, 315 (1990) Vafai, K.: Handbook of Porous Media. Marcel Dekker, New York (2000) Vafai, K.: Handbook of Porous Media. 2nd edn. Taylor and Francis(CRC), Boca Raton (2005) Nield, D.A., Bejan, A.: Convection in Porous Media. 3rd edn. Springer-Verlag, New York (2006) Eringen, A.C.: Theory of micropolar fluids. J. Math. Mech. 16, 1 (1966) Stokes, V.K.: Couple stresses in fluids. Phys. Fluids 91, 1709 (1966) Sharma, R.C., Sunil, Pal, M.: On couple-stress fluid heated from below in porous medium in the presence of rotation. Appl. Mech. Engg. 5(4), 883 (2000) Sharma, R.C., Sharma, M.: Effect of suspended particles on couple-stress fluid heated from below in the presence of rotation and magnetic field. Indian J. Pure Appl. Math. 35(8), 973 (2004) Malashetty, M.S., Basavaraja, D.: Effect of thermal/gravity modulation on the onset of Rayleigh–Benard convection in a couple stress fluid. Int. J. Trans. Phenomena 7, 31 (2005) Malashetty, M.S., Gaikwad, S.N., Swamy, M.: An analytical study of linear and non-linear double diffusive convection with Soret effect in couple stress liquids. Int. J. Therm. Sci. 45(9), 897 (2006) Gaikwad, S.N., Malashetty, M.S., Ramaprasad, K.: An analytical study of linear and non-linear double diffusive convection with Soret and Dofour effects in couple stress liquids. Int. J. Nonlinear Mech. 42(7), 903 (2007) Shivakumara, I.S.: Onset of convection in a couple-stress fluid–saturated porous medium: effects of non-uniform temperature gradients. Arch. Appl. Mech. 80, 949 (2010) Shivakumara, I.S., Suresh Kumar, S., Devaraju, N.: Coriolis effect on thermal convection in a couple-stress fluid–saturated rigid porous layer. Arch. Appl. Mech. (2010) doi:10.1007/s00419-010-0425-8 Hubbert, M.K.: Darcy’s law and the field equations of the flow of underground fluids. Trans. Am. Inst. Min. Met. Engrs. 207, 222 (1956) De Wiest, R.J.M.: Geohydrology. Wiley, New York (1965) Rudraiah, N., Srimani, P.K., Friedrich, R.: Finite amplitude convection in a two- component fluid saturated porous layer. Int. J. Heat Mass Transf. 25, 715 (1982) Shivakumara, I.S., Sumithra, R.: Non-Darcian effects on double diffusive convection in a sparsely packed porous medium. Acta Mech. 132, 113 (1999) Baines, P.G., Gill, A.E.: On thermohaline convection with linear gradients. J. Fluid Mech. 37, 289 (1969)