Stability of the horizontal throughflow of a power-law fluid in a double-diffusive porous layer under convective boundary conditions

International Journal of Thermal Sciences 146 (2019) 106098

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: http://www.elsevier.com/locate/ijts

Stability of the horizontal throughflow of a power-law fluid in a double-diffusive porous layer under convective boundary conditions Seema Kumari *, P.V.S.N. Murthy Department of Mathematics, Indian Institute of Technology Kharagpur, WB 721302, India

A R T I C L E I N F O

A B S T R A C T

Keywords: Porous medium Power-law fluid Biot numbers Viscous dissipation Double-diffusion Horizontal throughflow Stationary and oscillatory instability

Double-diffusive convective instability of a horizontal througflow in a power-law fluid-saturated porous layer lying between infinitely extended two plane horizontal boundaries is investigated by considering internal heating associated with the viscous dissipation. These boundaries are assumed to be impermeable, isosolutal and are subjected to the third kind thermal boundary conditions which are characterized by different Biot numbers B0 ; B1 . The linear stability analysis of basic flow for pseudoplastic and dilatant fluids is analyzed for the longitudinal and transverse rolls in both aiding and opposing buoyancy cases. The eigenvalue problem is solved numerically using the bvp4c in Matlab. It is shown that the longitudinal rolls are the most unstable mode for dilatant fluids, whereas the transverse rolls are the most unstable mode of disturbance for pseudoplastic fluids. In both these cases, an increase in either of these Biot numbers is enhancing the stability of power-law fluids. Increasing P�eclet number destabilized the transverse rolls for the pseudoplastic fluid while the longitudinal rolls are stabilized in case of the dilatant fluid. The effect of viscous dissipation on the onset of instability is also presented. A dis­ cussion on the onset of oscillatory convective instability (ω 6¼ 0) in this double-diffusive process, under third kind thermal boundary conditions is also presented. Further, the instability of power-law fluid is discussed for two limiting cases of thermal boundary conditions; the first one is B0 ¼ 0; B1 →∞ which corresponds to the constant heat flux at the bottom boundary while the second one is B0 →∞;B1 →∞, that corresponds to isothermal walls. As Pe ¼ 0 is a singularity for the eigenvalue problem, a detailed asymptotic analysis of Pe→0 is presented for all these cases.

1. Introduction Convection in non-Newtonian fluid-saturated porous medium has gained more attention due to its growing applications in food process­ ing, enhanced oil recovery and polymer industries, etc. However, the study of non-Newtonian fluid behavior is a more complex area of research. These fluids are frequently treated using the power-law model, viscoelastic fluid model, Maxwell fluid model etc. [1–3]. Considering a non-zero base flow of power-law fluid, Barletta and Nield [4] analyzed the convective instability of the Horten-Rogers-Lapwood problem. Also, using the more general Robin thermal boundary conditions (third kind conditions), this problem was extended by Alves and Barletta [5]. In this extension, all possible limiting cases of these Robin conditions and the convective instability of the power-law fluid was discussed in detail. Considering a vertical throughflow, the instability of power-law fluid was examined by Barletta and Storesletten [6], where the boundary planes are assumed to be permeable and isothermal. More recently,

Alves and Barletta [7] analyzed the transition of convective instability into absolute instability, with three-dimensional disturbances in both pseudoplastic and dilatant fluids. The results indicated that the pseu­ doplastic and Newtonian fluids have the same behavior whereas dilatant fluids show a different behavior. Nield [8] investigated the onset of double-diffusive convective instability in a horizontal porous channel. It was found that the oscil­ latory instability occurs when a strongly destabilizing thermal gradient opposed the stabilizing solutal gradient. Considering inclined thermal and solutal gradients along the horizontal boundaries, the same problem was extended by Nield et al. [9]. The double-diffusive fingering in a fluid-saturated porous medium was explained by Chen and Chen [10], where they used the Galerkin-finite difference method to solve the boundary value problem. Considering the Maxwell fluid model, Wang and Tan [11] studied the effect of relaxation time on the double-diffusive convective instability in a fluid-saturated porous layer. In case of non-Newtonian fluids, the viscous dissipation might become

* Corresponding author. E-mail addresses: [email protected] (S. Kumari), [email protected] (P.V.S.N. Murthy). https://doi.org/10.1016/j.ijthermalsci.2019.106098 Received 28 February 2019; Received in revised form 3 September 2019; Accepted 5 September 2019 Available online 18 September 2019 1290-0729/© 2019 Elsevier Masson SAS. All rights reserved.

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

obtained with n < 1, while n > 1 describes dilatant fluids. The gener­ alized Darcy’s law for power-law fluid-saturated porous medium may be written as

μ� �� ��n K�

u

1

u¼

(2)

rP:

In the above, P is the dynamic pressure and K� is a generalized permeability of power-law fluid defined as K� ¼

� �ðnþ1Þ=2 1 � nε �n 50K ; 2Ct 3n þ 1 3ε

(3)

where ε is porosity, K is the intrinsic permeability [18] and Ct is tortu­ osity factor equal to constant 25/12 [19]. For Newtonian fluid (n ¼ 1), is equal to Kμ , where μ is dynamic viscosity of Newtonian fluid.

μ�

K�

Fig. 1. Schematic diagram of the power-law fluid-saturated porous channel with horizontal throughflow.

2.1. Governing equations The power-law fluid-saturated porous medium is assumed to be homogeneous, isotropic, and also it is in local thermal equilibrium. The linear Oberbeck-Boussinesq approximation is employed in the general­ ized Darcy’s equation which governs the fluid flow in the porous channel. Under these assumptions, the governing equations for doublediffusive convection with viscous dissipation may be written as

significant because of heat generation within the medium. For Newto­ nian fluids, the convective instability of base flow with the effect of viscous dissipation in a porous medium has been explained by many researchers [12–14]. Considering the effect of viscous dissipation, Celli and Barletta [15] investigated the onset of thermal instability of power-law fluid in a horizontal porous layer confined between imper­ meable isothermal upper boundary and adiabatic lower boundary. The combined effect of viscous dissipation and double-diffusion on the instability of Newtonian fluid flow through porous media was analyzed by Barletta and Nield [16], where boundary planes are assumed to be impermeable, with adiabatic lower boundary and isothermal upper wall. In this study, they considered the disturbances in the form of oblique rolls. The purpose of the present study is to investigate the significance of third kind thermal boundary conditions on the onset of double-diffusive convective instability of a horizontal througflow of power-law fluids in a horizontal porous layer with viscous dissipation in the medium. This eigenvalue problem is solved numerically using the bvp4c routine in Matlab. The onset of instability for both longitudinal and transverse rolls is investigated in both aiding and opposing buoyancy cases. The doublediffusive convective instability is also discussed by considering two limiting cases of third kind thermal boundary conditions. Also, the sin­ gular solution of the eigenvalue problem at Pe ¼ 0 is examined.

μ� �� ��n K�

u

1

u¼

ρ0 g ½βT ðT

rP

T0 Þ þ βC ðC

(5)

C0 Þ�;

σ

∂T μ� �� ��nþ1 2 u ; þ u:rT ¼ χ r T þ � ∂t K ρ0 c p

(6)

φ

∂C 2 þ u:rC ¼ D r C; ∂t

(7)

where βT and βC are thermal and solutal expansion coefficients of the fluid, σ is the ratio between the volumetric heat capacity of the porous medium and that of the fluid and ρ0 is fluid density at the reference temperature T0 and reference concentration C0 respectively. In the above, ΔT is the temperature difference and ΔC is the concentration difference across the channel, χ is the effective thermal diffusivity of the porous medium, cp is the specific heat of the fluid, φ and D are porosity and mass diffusivity in the medium, respectively. The last term in energy balance equation is the modified viscous dissipation term for nonNewtonian power-law fluid [15,20]. The corresponding boundary con­ ditions are written as

2. Mathematical model We consider the Ostwald-de Waele ([17,18]) power-law fluid-­ saturated plane horizontal porous channel of thickness H bounded be­ tween two impermeable and isosolutal boundary walls situated at z ¼ 0 and z ¼ H (see Fig. 1). Convective thermal boundary conditions are considered at both horizontal boundary walls which induces a vertical thermal gradient in the channel. Both these lower and upper walls ex­ change heat between the internal and external fluid environment, hav­ ing reference temperature T0 þ ΔT and T0 at z ¼ 0 and z ¼ H, such that T0 þ ΔT > T0 (refer to Fig. 1). Here, q0 is the wall heat flux. The grav­ itational acceleration g ¼ gbe z is oriented in the opposite direction of the vertical z-axis. Also ðu; v; wÞ are the components of seepage velocity u in x; y and z directions respectively. For the non-Newtonian power-law fluid, the shear stress τ is a function of shear strain rate γ_ which is given by

τ ¼ μ* ð_γÞn ; � �n 1 ¼ μ* �γ_ � γ_

(4)

r:u ¼ 0;

z¼0:

w ¼ 0;

z¼H:

w ¼ 0;

k k

∂T þ h0 ðT ∂z

∂T þ h1 ðT ∂z

ðT0 þ ΔTÞÞ ¼ q0 ; T0 Þ ¼

q0 ;

C ¼ C0 þ △C;

C ¼ C0 ; (8)

where h0 and h1 are convection heat transfer coefficients and k is ther­ mal conductivity. The reference temperature difference (ΔT) is defined as ΔT ¼ q0kH. 2.2. Dimensionless formulation Consider the following non-dimensionalization for the physical variables as:

(1)

� �n 1 where μ� �γ_ � is the apparent viscosity and μ� represents consistency factor of the fluid and n is the power-law index. In case of a Newtonian fluid, power-law index n ¼ 1 and pseudoplastic fluid behavior is 2

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 2. Variation of critical Rayleigh number Rac with power-law index n in aiding buoyancy ðN ¼ 1Þ for Pe ¼ 0:5; Ge ¼ 0:01 with different values of Le. Table 1 Comparison of the present results with the data available in Alves and Barletta [5] for Ge ¼ 0; N ¼ 0; n ¼ 1. Alves and Barletta [5]

Λ

Present results (bvp4c)

Rac

∂C 1 þ u:rC ¼ r2 C; Le ∂t

z¼0:

w ¼ 0;

z¼1:

w ¼ 0;

(13)

∂T þ B0 T ¼ 1 þ B0 ; ∂z ∂T þ B1 T ¼ 1; ∂z

C ¼ 1;

(14)

Longitudinal rolls

Transverse rolls

Rac

ac

Rac

ac

3.4 � 10 4

12.000044

2.6 � 10 4

non-Newtonian fluid is defined as

appear in Eq. (11) and Eq. (13) are due to solute transport in the system. The Lewis number Le is the ratio of thermal and solutal diffusivity (Le ¼

12.000073

C ¼ 0;

where be z is unit vector along z direction. The Rayleigh number Ra for

B0 ¼ 10 12 ; B1 ¼ 10 12

12

B0 ¼ 0:1; B1 ¼ 0:1

15.3912

15.391169

1.201175

15.391170

1.201174

B0 ¼ 1; B1 ¼ 1

22.352

22.352284

2.056626

22.352284

2.056656

B0 ¼ 10; B1 ¼ 10

33.858

33.860135

2.879147

33.860135

2.879148

B1 ¼ h1kH respectively.

B0 ¼ 0:1; B1 ¼ 10

25.281

25.281196

2.240794

25.281196

2.240794

2.3. Basic solution

B0 ¼ 1; B1 ¼ 10

27.960

27.960389

2.479446

27.960389

2.479446

Consider the uniform horizontal througflow along x and y axes in the porous channel, which is the steady state solution for the flow filed. This is represented as

number Ge is defined as

D)

χ

gβH cp .

ρ0 gβT △TK� Hn , μ� χ n

Λ ¼ φσ and the Gebhart

The non-dimensional numbers N and Le

and buoyancy ratio N ¼ ββc ΔC ΔT is the ratio of thermal and solutal T

buoyancy forces. B0 and B1 are the Biot numbers defined as B0 ¼ h0kH and

(15)

uB ¼ ðPe cos α; Pe sin α; 0Þ;

where the subscript B represents basic flow field, Pe ¼ Uχo H is the P�eclet

1 H H χ 2 ðx; y; zÞ→ðx; y; zÞ; u ¼ ðu; v; wÞ→u; t→t; H r→r; H 2 r →r2 H χ χ σH2 T

T0 C C0 →T; →C: ΔT ΔC

number with Uo as the dimensional basic flow velocity and α is an angle of inclination of this basic velocity made with the x-axis. Corresponding to this basic velocity profile, we obtain the temperature profile which depends on the pair of Biot numbers and the concentration profile which is a linear function of z only, as:

(9)

The pressure eliminated governing Eqs. (4)–(7) and the corre­ sponding boundary conditions (8) using the above non-dimensional variable may be written as r:u ¼ 0;

(10)

� �n 1 r � �u� u ¼ Rar � ðT þ N CÞb ez;

(11)

∂T Ge� �nþ1 þ u:rT ¼ r2 T þ �u� ; Ra ∂t

(12)

TB ¼ 1

z

Ge Penþ1 ð ¼1

z:

2 þ B0 ð 2 þ zÞz þ B1 ð 1 þ zÞð1 þ z þ B0 zÞÞ ; 2ðB0 þ B1 þ B0 B1 ÞRa

CB (16)

The basic temperature profile for all limiting cases of the Biot numbers, with no viscous dissipation in the medium coincides with the one obtained by Alves and Barletta [5]. When we consider non-zero Gebhart number (Ge 6¼ 0), the basic temperature profile has a 3

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 3. Variation of critical Rayleigh number Rac with power-law index n for transverse rolls in aiding buoyancy ðN ¼ 1Þ at Pe ¼ 10; Ge ¼ 0:01.

Fig. 4. Variation of critical Rayleigh number Rac with power-law index n for longitudinal rolls in aiding buoyancy ðN ¼ 1Þ at Pe ¼ 10; Ge ¼ 0:01.

singularity at B0 ¼ 0; B1 ¼ 0 (constant heat flux boundary at z ¼ 0 and z ¼ H). Also, this case leads to an ill-posed problem (Eqs. (10) to (14)).

� ∂V n sin2 α þ cos2 α ∂z � � Ra ∂θ ∂Γ ¼� n 1 þN ; �Pej ∂y ∂y

∂W ∂y

3. Linear stability analysis The onset of instability of this system is investigated by superposing small amplitude disturbances on the basic velocity, temperature and concentration profiles as u ¼ uB þ ε U; T ¼ TB þ ε θ; C ¼ CB þ ε Γ

� ∂U ∂W n cos2 α þ sin2 α þ ðn ∂z ∂x � � Ra ∂θ ∂Γ ¼ � n 1 þN ; �Pej ∂x ∂x

(17)

where ε≪1 is a small perturbation parameter and U ¼ ðU; V; WÞ. Substituting expressions given in Eq. (17) into Eqs.(10)–(14) and neglecting the non-linear terms of order ε2 , we obtain the system of equations in terms of the perturbation quantities as

∂U ∂V ∂W þ þ ¼ 0; ∂x ∂y ∂z

ðn

� ∂V n sin2 α þ cos2 α ∂x � ∂V ¼ 0; ∂y

1Þsin α cos α

∂U ∂z (19)

1Þ

∂V sin α cos α ∂z

� ∂U n cos2 α þ sin2 α þ ðn ∂y

(20) � 1Þsin α cos α

∂U ∂x

(21)

(18)

4

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 5. Variation of critical Rayleigh number Rac with power-law index n in opposing buoyancy ðN ¼

0:1Þ for Pe ¼ 0:5; Ge ¼ 0:01 with different values of Le.

Fig. 6. Variation of critical Rayleigh number Rac with power-law index n for transverse rolls in opposing buoyancy ðN ¼

39.4784175 and corresponding wave number is ac ¼ 3:141593. For two dimensional solution of Eqs.(18)–(24), we consider arbitrary distur­ bances in the velocity, temperature and concentration profiles as invariant along the y-direction (are functions of ðx; z; tÞ only). The angle of inclination α is varied from 0 to 2π. Eqs. (19) and (21) enable us to have an explicit relation between U and V as

∂θ ∂θ ∂θ dTB þ Pe cos α þ Pe sin α þ W ∂t ∂x ∂y dz ¼ r2 θ þ

ðn þ 1ÞGe Pen ðU cos α þ V sin αÞ; Ra

∂Γ ∂Γ ∂Γ dCB 1 þ Pe cos α þ Pe sin α þ W ¼ r2 Γ; Le ∂t ∂x ∂y dz

(22) (23)

V¼

and the boundary conditions are z¼0:

W ¼ 0;

z¼1:

W ¼ 0;

∂θ þ B0 θ ¼ 0; ∂z ∂θ þ B1 θ ¼ 0; ∂z

Γ ¼ 0;

0:1Þ at Pe ¼ 10; Ge ¼ 0:01.

ð1 nÞsin α cos α U: n sin2 α þ cos2 α

(25)

Using the above Eq. (25), the governing Eqs.(18)–(23) with bound­ ary conditions become � � n ∂U ∂W Ra ∂θ ∂Γ � ¼ þN ; (26) 2 n 1 2 � ∂x ∂x ∂x n sin α þ cos α ∂z Pej

(24)

Γ ¼ 0:

Eqs. (19) and (20) have singularity at Pe ¼ 0. An asymptotic analysis of Pe→0 for longitudinal and transverse rolls with convective boundary conditions and solute transport in a system is described in section 8. In case of pure thermal convection with Pe→0, the present numerical re­ sults for isothermal boundaries coincide with the one of generalized Horton-Rogers-Lapwood Problem for power-law fluid. In this case, the critical Rayleigh number Rac of a Newtonian fluid (n ¼ 1) is equal to

∂θ ∂θ dTB ∂2 θ ∂2 θ ðn þ 1ÞGe Pen þ Pe cos α þ W ¼ 2þ 2þ ðU cos α þ V sin αÞ; ∂t ∂x dz ∂x ∂z Ra

(27)

5

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 7. Variation of critical Rayleigh number Rac with power-law index n for longitudinal rolls in opposing buoyancy ðN ¼

0:1Þ at Pe ¼ 10; Ge ¼ 0:01.

Fig. 8. Variation of critical Rayleigh number Rac with Biot number B0 in aiding buoyancy N ¼ 1 for Le ¼ 0:1; Pe ¼ 5; Ge ¼ 0:01 with different values of B1 .

�

�

∂Γ ∂Γ dCB 1 ∂2 Γ ∂2 Γ þ ; þ Pe cos α þ W ¼ Le ∂x2 ∂z2 ∂t ∂x dz z¼0:

W ¼ 0;

z¼1:

W ¼ 0;

∂θ þ B0 θ ¼ 0; ∂z ∂θ þ B1 θ ¼ 0; ∂z

Γ ¼ 0;

(28)

∂θ ∂θ þ Pe cos α ∂t ∂x

þ (29)

Γ ¼ 0:

∂Γ ∂Γ þ Pe cos α ∂t ∂x

Introducing the stream function ψ which is defined as U¼

∂ψ ; ∂z

W¼

∂ψ ; ∂x

∂ψ dTB ∂2 θ ∂2 θ ¼ 2þ 2 ∂x dz ∂x ∂z

(30)

� � ðn þ 1ÞGe Pen ∂ψ cos α ; ∂z n sin2 α þ cos2 α Ra �

the system of partial differential equations and boundary conditions are written as � � ∂2 ψ n ∂2 ψ Ra ∂θ ∂Γ � þ ¼ þ N ; (31) �Pejn 1 ∂x ∂x2 n sin2 α þ cos2 α ∂z2 ∂x

�

∂ψ dCB 1 ∂2 Γ ∂2 Γ þ ; ¼ ∂x dz Le ∂x2 ∂z2

∂ψ ∂θ ¼ 0; þ B0 θ ¼ 0; Γ ¼ 0; ∂x ∂z ∂ψ ∂θ z¼1: ¼ 0; þ B1 θ ¼ 0; Γ ¼ 0: ∂x ∂z z¼0:

(32) (33)

(34)

The marginal stability of a system is described by using wavelike solution of disturbance Eqs.(31)–(34) as consider in Ref. [16]. � � ψ ¼ R if ðzÞeiðaxþtðω a Pe cosαÞÞ ; θ ¼ R hðzÞeiðaxþtðω a Pe cosαÞÞ ; Γ � (35) ¼ R gðzÞeiðaxþtðω a Pe cosαÞÞ ;

6

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 9. Variation of critical Rayleigh number Rac with Biot number B0 in opposing buoyancy N ¼

0:1 for Le ¼ 0:1; Pe ¼ 5; Ge ¼ 0:01 with different values of B1 .

Fig. 10. Variation of critical Rayleigh number Rac with Biot number B1 in aiding buoyancy N ¼ 1 for Le ¼ 0:1; Pe ¼ 5; Ge ¼ 0:01 with different values of B0 .

where real parameter ω is an angular frequency, a is the wave number and R denotes the real part of the expression. Substituting these ex­ pressions given in Eq. (35) into Eqs.(31)–(34), we obtain the eigenvalue problem as, n f 00 n sin2 α þ cos2 α h00

g

00

a2 h

2

a g

z¼0: z¼1:

iωh

a2 f ¼

dTB af ¼ dz

Ra a � n 1 ðh þ N gÞ; �Pej � � iðn þ 1ÞGe Pen cos α 0 f ; 2 Ra n sin α þ cos2 α

dCB i ω Le g ¼ a Le f ; dz f ¼ 0; f ¼ 0;

0

h þ B0 h ¼ 0; g ¼ 0; 0 h þ B1 h ¼ 0; g ¼ 0:

3.1. Longitudinal and transverse rolls The system of coupled ordinary differential equations involving complex valued functions reduces to a real system of ordinary differ­ ential equation for longitudinal rolls ðα ¼ π=2Þ with ω ¼ 0. These equations are

(36)

a2 f ¼

h00

a2 h ¼ a f

g00

a2 g ¼ a Le f

(37) (38) (39)

Ra a � n 1 ðh þ N gÞ; �Pej

f 00

dTB ; dz

0

dCB ; dz

z ¼ 0 : f ¼ 0; h þ B0 h ¼ 0; g ¼ 0; 0 z ¼ 1 : f ¼ 0; h þ B1 h ¼ 0; g ¼ 0:

In wavelike disturbances defined by Eq. (35), the longitudinal rolls corresponds to α ¼ π=2 and the transverse rolls are obtained with α ¼ 0.

(40) (41) (42) (43)

This real system of equations allows one to determine the value of 7

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 11. Variation of critical Rayleigh number Rac with Biot number B1 in opposing buoyancy N ¼

0:1 for Le ¼ 0:1; Pe ¼ 5; Ge ¼ 0:01 with different values of B0 .

Fig. 12. Variation of critical angular frequency ωc with n in aiding buoyancy N ¼ 1 for Le ¼ 5; Pe ¼ 10; Ge ¼ 0:01.

real parameter Ra corresponding to a. For transverse rolls (α ¼ 0), with ω 6¼ 0, the resulting eigenvalue problem is given by Eqs.(36)–(39) which involves complex functions fðzÞ, hðzÞ and gðzÞ. The real and imaginary parts of these equations are separated, which yields a six second order ordinary differential equations. These equations along with the bound­ ary conditions and the normalization condition are solved to determine Ra and ω.

R2014b. To gain higher-order accuracy, the relative and absolute tolerance are considered as 10 6 and 10 10 respectively. The results obtained for both the longitudinal and transverse rolls are depicted in Figs. 2–19. 5. Result and discussion The double-diffusive convective instability in a power-law fluidsaturated porous medium with horizontal throughflow considering the effect of viscous dissipation and third kind thermal boundary conditions is governed by the diffusivity ratio Le, buoyancy ratio N, Biot numbers B0 and B1 , the P�eclet number Pe and viscous dissipation parameter Ge. The following physically realistic range of these parameters is consid­ ered; 0 � n � 2 ([6]), 1 < N � 1 and 0 < Le � 10 ([21]), 0 < Pe � 10 ([14]), 0 � Ge � 0:1 ([13,14]) and 0 < B0 ; B1 � 100 ([5]). The limiting cases B0 ¼ 0; B1 →∞ and B0 →∞; B1 →∞ are discussed in section 6 and 7 respectively. The numerical results obtained here for the pure thermal convection problem with convective thermal boundary conditions for Ge ¼ 0 and n ¼ 1 are in excellent agreement with the results presented in Ref. [5], this comparison is shown in (Table 1). In

4. Numerical solution The eigenvalue problem defined by Eqs.(36)–(39) is solved using the bvp4c routine in MATLAB R2014b. In this numerical method, we need to convert the governing equations and boundary conditions into a system of first-order ordinary differential equations. For a non-trivial solution of this problem, we consider the normalization condition 0 f ð0Þ ¼ 1. Using this normalization condition, we determine the values of unknown real parameters Ra and ω. After obtaining eigenvalue Ra, for every assigned value of input parameters ðN;Le;Ge;Pe;B0 ;B1 Þ, the critical value of Rac and the corresponding wavenumber ac are obtained by minimizing the function RaðaÞ using indexmin command in MATLAB 8

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 13. Variation of critical angular frequency ωc with n in opposing buoyancy N ¼

0:1 for Le ¼ 5; Pe ¼ 10; Ge ¼ 0:01.

Fig. 14. Variation of critical Rayleigh number Rac with n in aiding buoyancy N ¼ 1 for Le ¼ 1; Pe ¼ 10 with different values of Ge.

this pure thermal convection case with Ge ¼ 0, Newtonian fluids are equally unstable (same value of Rac ) for both longitudinal and trans­ verse rolls. In the following subsections, we present the significance of convective boundary conditions on the onset of double-diffusive convective instability for the (a) aiding buoyancy (considering N ¼ 1) (b) opposing buoyancy (N ¼ 0:1).

N, Le and Ge. When Ge ¼ 0 and no solute concentration present in the system, the stability of Newtonian fluid is independent of the angle of inclination, α. In that case, transverse and longitudinal rolls are equiv­ alent, represented by an equal value of Rac . A close comparison between these figures indicates the impact of the diffusivity ratio on the insta­ bility of flow, at a fixed value of Pe ¼ 0:5. The destabilizing character of fluid is increased when thermal diffusivity dominates over the solutal diffusivity. These figures also explain the impact of Biot numbers B0 and B1 on the convective instability of power-law fluid in a porous medium. It reveals how the imperfect heat transfer at the boundaries (i.e., finite heat transfer coefficient), affect the instability of basis state solution. The pseudoplastic fluids (n < 1) exhibit non-uniform trend for the transverse rolls while the longitudinal rolls show a continuous decrease in the critical Rayleigh number for all sets of values of B0 and B1 . In case of B0 ¼ 0:1; B1 ¼ 10 and B0 ¼ 10; B1 ¼ 0:1, the power-law fluids are equally unstable at the low shear rate, which can be noticed from the equal values of Rac . The stability of both the longitudinal and transverse rolls decreased with the dominance of thermal diffusivity over the sol­ utal diffusivity, this is seen in all four sets of values considered for the

5.1. Aiding buoyancy The effect of buoyancy ratio under convective boundary conditions on the onset of instability of power-law fluid through the porous channel with horizontal throughflow is shown for low diffusivity ratio (Le ¼ 0:1) in Fig. 2(a) and for high diffusivity ratio (Le ¼ 5) in Fig. 2(b). The letter “L” stands for the longitudinal rolls and “T” for the transverse rolls. It is noticed that the pseudoplastic and dilatant fluids exhibit different behaviors in aiding buoyancy situation. Transverse rolls (α ¼ 0) are the most unstable mode of disturbance for the pseudoplastic fluids (n < 1), and the longitudinal rolls (α ¼ π=2) for dilatant fluids (n > 1). In case of a Newtonian fluid (n ¼ 1), instability of the basic state solution for longitudinal and transverse rolls depends on the angle of inclination α, 9

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 15. Variation of critical Rayleigh number Rac with n in opposing buoyancy N ¼

0:1 for Le ¼ 1; Pe ¼ 10 with different values of Ge.

Table 2 Critical value of Rayleigh number and corresponding wave number for transverse roll at Pe ¼ 10; Le ¼ 1; B0 ¼ B1 ¼ 100 with effect of viscous dissipation. n

0.5 0.7 1 1.5 2

N ¼ 1; Ge ¼ 10

N ¼ 1; Ge ¼ 0:01

6

4.496039 8.255166 19.54853 76.56882 285.3300

N¼

0:1; Ge ¼ 10

N¼

6

0:1; Ge ¼ 0:01

ac

Rac

ac

Rac

ac

Rac

ac

2.628599 2.859276 3.125988 3.459496 3.717513

4.4959776 8.2550804 19.548385 76.568379 285.32868

2.628544 2.859280 3.126001 3.459509 3.717514

9.8535035 18.109892 42.930559 168.35623 627.90138

2.613052 2.842345 3.107418 3.438910 3.695333

9.8533568 18.109696 42.930222 168.35530 627.89838

2.612894 2.842256 3.107423 3.438906 3.695339

Table 3 Critical value of Rayleigh number and corresponding wave number for longitudinal Roll at Pe ¼ 10; Le ¼ 1; B0 ¼ B1 ¼ 100 with effect of viscous dissipation. n

N ¼ 1; Ge ¼ 10

0.5 0.7 1 1.5 2

6.1817905 9.7974789 19.548539 61.817909 195.48539

N ¼ 1; Ge ¼ 0:01

6

N¼

0:1; Ge ¼ 10

N¼

6

0:1; Ge ¼ 0:01

ac

Rac

ac

Rac

ac

Rac

ac

3.125972 3.125977 3.125980 3.125976 3.125976

6.18176514 9.79743729 19.5484570 61.8176493 195.484570

3.125990 3.125988 3.125985 3.125983 3.125985

13.5758345 21.5162483 42.9305591 135.758347 429.305591

3.107414 3.107412 3.107414 3.107414 3.107414

13.57577736 21.51615803 42.93037953 135.7577795 429.3037939

3.107422 3.107421 3.107423 3.107424 3.107423

Biot numbers. Variation of the critical Rayleigh number Rac with power-law index n for both the transverse and longitudinal rolls are presented for high shear rate (Pe ¼ 10) in Figs. 3 and 4 respectively. The value of Rac tends to zero uniformly as n tends to 0. At the high shear rate, the apparent viscosity of the pseudoplastic fluids tends to zero and it is infinite for dilatant fluid. From these figures, it can also be noticed that the longi­ tudinal rolls are the most unstable mode of disturbance for dilatant fluids. The effect of Lewis number on the critical value of thermal Rayleigh number Rac and power-law index n at high shear rate are also shown in Figs. 3 and 4. It is observed that increasing the Lewis number decreased Rac and it leads to increased destabilization of the basic state solution in the horizontal porous channel.

3 in opposing buoyancy as compared from aiding buoyancy. The reason is that the fluid flow is getting stabilize by the solutal buoyancy component while thermal buoyancy, which is oppositely directed is to destabilize this action. It is explained by the comparison of Figs. 2 and 5 for different pairs of Biot number in aiding and opposing buoyancy cases respectively. To describe the combined effect of thermal and solutal diffusivity on the instability of base flow at the low shear rate (Pe ¼ 0:5) is plotted in Fig. 5(a) and (b). From these two figures, it is inferred that in opposing buoyancy situation, the dominance of thermal diffusivity over solutal diffusivity will stabilize the basic state solution in both longitu­ dinal and transverse rolls. However in this scenario, at Le ¼ 5 and B0 ¼ 10; B1 ¼ 10, longitudinal and transverse rolls are not equally unstable for a Newtonian fluid, having different values of Rac . The variation of Rac for all n at high shear rate (Pe ¼ 10) is illustrated in Figs. 6 and 7 for transverse and longitudinal rolls respectively for N ¼ 0:1. In this case also, the value of Rac tends to zero from a very high value, as the power-law index approaches the pseudoplastic limit (n→0) for both longitudinal and transverse rolls. But the slope of these curves is higher compared to those in the case of aiding buoyancy, at the value of n ¼ 1þ . It means that at the high shear rate, the basic solution is more

5.2. Opposing buoyancy In opposing buoyancy, which is contrary to the case of aiding buoyancy (N < 0), Fig. 5(a) and (b) demonstrates the convective insta­ bility of power-law fluid for both longitudinal and transverse rolls. For Pe ¼ 0:5 and Ge ¼ 0:01, the flow is getting stable as n is varies from 0 to 10

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 16. Variation of critical Rayleigh number Rac with n for constant heat flux boundary for aiding buoyancy (N ¼ 1) at Ge ¼ 0:01.

Fig. 17. Variation of critical Rayleigh number Rac with n for constant heat flux boundary for opposing buoyancy (N ¼

0:1) at Ge ¼ 0:01.

Table 4 Critical value of Rayleigh number and wave number for transverse rolls at Pe ¼ 10; Le ¼ 1 for B0 ¼ 0; B1 →∞. n

0.5 0.7 1 1.5 2

N ¼ 1; Ge ¼ 10

3.5665560 6.7593464 16.509731 66.734280 253.65363

N ¼ 1; Ge ¼ 0:01

6

N¼

0:1; Ge ¼ 10

N¼

6

0:1; Ge ¼ 0:01.

ac

Rac

ac

Rac

ac

Rac

ac

2.184062 2.404155 2.663136 2.992197 3.249532

3.4809879 6.6255157 16.246396 65.913775 251.08332

2.178844 2.399463 2.659138 2.988978 3.246872

5.8996091 11.476022 28.838284 120.43108 468.32183

1.905818 2.075257 2.274399 2.529230 2.731505

5.7606412 11.252992 28.387110 118.98079 463.67972

1.906887 2.076627 2.276034 2.531101 2.733501

stable in the case of the opposing buoyancy. This is seen for both transverse and longitudinal rolls. Also, the longitudinal rolls are seen to be the most unstable mode of disturbance for dilatant fluid and the transverse rolls are the most unstable ones for the pseudoplastic fluid. For higher values of the Lewis number (Le ¼ 5), the stabilizing character of a basic solution increases at the high shear rate (for example Pe ¼ 10).

5.3. Effect of convective boundaries A detailed description of the effect of convective boundary condi­ tions on the instability of power-law fluid is shown in Figs. 8–11 for both aiding and opposing buoyancy cases. At the upper boundary z ¼ H, a smaller value of Biot number B1 indicates that the heat conduction in­ side the fluid-saturated porous medium is faster than the heat 11

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Table 5 Critical value of Rayleigh number and wave number for longitudinal rolls at Pe ¼ 10; Le ¼ 1 for B0 ¼ 0; B1 →∞. n

0.5 0.7 1 1.5 2

N ¼ 1; Ge ¼ 10

5.2208358 8.2744662 16.509731 52.208354 165.09731

N ¼ 1; Ge ¼ 0:01

6

N¼

0:1; Ge ¼ 10

N¼

6

0:1; Ge ¼ 0:01.

ac

Rac

ac

Rac

ac

Rac

ac

2.663149 2.663145 2.663149 2.663150 2.663151

5.1376075 8.1425602 16.246542 51.376079 162.46542

2.659136 2.659150 2.659147 2.659146 2.659147

9.1194659 14.453379 28.838284 91.194662 288.38284

2.274405 2.274397 2.274402 2.274400 2.274400

8.9768656 14.227373 28.387342 89.768659 283.87342

2.276049 2.276044 2.276044 2.276046 2.276046

Fig. 18. Variation of critical Rayleigh number Rac with n for isothermal boundaries for aiding buoyancy (N ¼ 1) at Ge ¼ 0:01.

Fig. 19. Variation of critical Rayleigh number Rac with n for isothermal boundaries for opposing buoyancy (N ¼

convection at the surface of the porous layer, the instability of basic state solution is decreased with the increasing value of B0 , this is shown in Fig. 8(a) for the aiding buoyancy situation. When the heat convection at the surface dominates the effect of heat conduction in the interior of the porous channel (indicated with large value of B1 ), the stabilizing char­ acter of power-law fluid increased, this can be seen from Fig. 8(b) for aiding buoyancy (N ¼ 1). Similar behavior is observed for the opposing buoyancy situation as well, which is evident from Fig. 9(a) and (b). The influence of different values of B0 on convective instability is

0:1) at Ge ¼ 0:01.

shown in Figs. 10 and 11 for aiding and opposing buoyancy respectively, considering Le ¼ 0:1, Pe ¼ 5 and Ge ¼ 0:01. From these figures it can be inferred that, at the lower boundary ðz ¼ 0Þ when heat conduction in­ side the porous channel dominates the effect of convection on the sur­ face of porous layer (B0 ≪1), the curves nature is similar to that has been given in Figs. 8 and 9 in both aiding and opposing buoyancy. It means that interchange the dominance of conduction or convection at both the boundaries is not affecting the instability of the system at both longi­ tudinal and transverse rolls (see Fig. 10 for aiding buoyancy and Fig. 11 12

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Table 6 Critical value of Rayleigh number and wave number for transverse rolls at Pe ¼ 10; Le ¼ 1; N ¼ 1 for B0 →∞; B1 →∞. n

N ¼ 1; Ge ¼ 10

0.5 0.7 1 1.5 2

N ¼ 1; Ge ¼ 0:01

6

4.5476522 8.3430407 19.739192 77.237823 287.62135

N¼

0:1; Ge ¼ 10

N¼

6

0:1; Ge ¼ 0:01.

ac

Rac

ac

Rac

ac

Rac

ac

2.641736 2.873288 3.141333 3.476843 3.735997

4.5476341 8.3430138 19.739058 77.237391 287.62002

2.641752 2.873610 3.141579 3.476845 3.735600

10.1058918 18.5402189 43.8648930 171.639663 639.158557

2.641671 2.873557 3.141607 3.476753 3.735998

10.105758 18.540034 43.864552 171.63870 639.15561

2.641677 2.873563 3.141612 3.476757 3.736000

Table 7 Critical value of Rayleigh number and wave number for longitudinal rolls at Pe ¼ 10; Le ¼ 1; N ¼ 1 for B0 →∞; B1 →∞. n

N ¼ 1; Ge ¼ 10

0.5 0.7 1 1.5 2

N ¼ 1; Ge ¼ 0:01

6

6.2420855 9.8930393 19.739208 62.420858 197.39208

N¼

0:1; Ge ¼ 10

N¼

6

0:1; Ge ¼ 0:01

ac

Rac

ac

Rac

ac

Rac

ac

3.141604 3.141593 3.141588 3.141590 3.141591

6.2420610 9.8929995 19.739130 62.420611 197.39130

3.141595 3.141601 3.141597 3.141601 3.141600

13.8713019 21.9845318 43.8649086 138.713019 438.649084

3.141593 3.141590 3.141592 3.141591 3.141591

13.8712469 21.9844451 43.8647341 138.712469 438.647342

3.141602 3.141599 3.141601 3.141600 3.141601

Table 8 Value of λ0 and λ2 for longitudinal rolls at N ¼ 1, Ge ¼ 0:001 and Le ¼ 1. a ¼ 0:5 λ0 B0 ¼ 0:5;B1 ¼ 0:5

a ¼2 λ2

λ0

a ¼ 3:5 λ2

a ¼5

λ0

λ2

λ0

λ2

49.101688

0.004010

14.027510

0.002923

16.096064

0.012512

21.949467

0.000278

B0 ¼ 1; B1 ¼ 1

74.335781

0.163319

15.300422

0.004320

16.474470

0.011079

22.136461

0.000363

B0 ¼ 10; B1 ¼ 10

171.651171

0.002803

21.392871

18.705223

0.005611

23.417997

0.001485

B0 ¼ 10; B1 ¼ 1

122.242090

0.005318

18.208589

17.533642

0.011863

22.750694

0.003383

17.533643

0.004042

22.750700

B0 ¼ 1; B1 ¼ 10

122.241989

0.032426

Table 9 Value of λ0 and λ2 for longitudinal rolls at N ¼

18.208586

0.032684 0.002880 0.067747

0.001668

0:1, Ge ¼ 0:001 and Le ¼ 1.

a ¼ 0:5

a ¼2

a ¼ 3:5

a ¼5

λ0

λ2

λ0

λ2

λ0

λ2

λ0

λ2

B0 ¼ 0:5;B1 ¼ 0:5

56.471279

0.000643

20.565189

0.009241

28.727128

0.000006

43.293371

0.002516

0.024994

B0 ¼ 1; B1 ¼ 1

92.693737

B0 ¼ 10; B1 ¼ 10

318.221653

B0 ¼ 1; B1 ¼ 10

181.092366

B0 ¼ 10; B1 ¼ 1

181.092266

29.988836

0.000012

44.070317

0.002075

0.230831

41.855678

23.445119

0.000459

38.538207

0.000053

49.728145

0.000913

0.002590

31.016576

0.000754

33.698373

0.000505

46.617045

0.000701

46.617065

0.001577

0.002721

0.027803

31.016584

0.001917

33.698378

0.000381

Table 10 Value of λ0 and λ2 for transverse rolls at N ¼ 1, Ge ¼ 0:001; Le ¼ 1 and a ¼ 3:5. n ¼ 0:5

B0 ¼ 0:5;B1 ¼ 0:5

n ¼1

n ¼2

λ0

λ2

λ0

λ2

λ0

λ2

12.480699

0.000789

16.096051

0.000000

23.303616

12.778703

0.000822

16.474467

0.000000

23.846567

B0 ¼ 10; B1 ¼ 10

14.529074

0.000668

18.705211

0.000000

27.054588

0.004582

B0 ¼ 10; B1 ¼ 1

55.187766

0.051331

17.533639

0.120121

25.376656

2.722981

B0 ¼ 1; B1 ¼ 10

55.187449

17.533645

4.822761

25.376660

39.784024

B0 ¼ 1; B1 ¼ 1

1.507289

for opposing buoyancy).

0.003488 0.007736

describes the onset of oscillatory convective instability of power-law fluid. In case of longitudinal rolls (α ¼ π=2), the modified angular fre­ quency ω is equal to zero. It means that longitudinal rolls always lead to the onset of stationary convective instability. However in case of transverse rolls, the non-zero value of ω describes oscillatory convective

5.4. Description of angular frequency The non-zero value of modified angular frequency ωð ¼ ω aPecosαÞ 13

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Table 11 Value of λ0 and λ2 for transverse rolls at N ¼

0:1, Ge ¼ 0:001; Le ¼ 1 and a ¼ 3:5.

n ¼ 0:5

n ¼1

λ0

λ2

B0 ¼ 0:5;B1 ¼ 0:5

22.178285

0.000038

23.174265

0.026424

B0 ¼ 1; B1 ¼ 1

λ0

n ¼2 λ2

λ0

28.727208

0.000000

41.692989

0.024825

29.988897

0.000000

43.500755

0.029537

0.000000

55.762179

0.035828

B0 ¼ 10; B1 ¼ 10

29.913392

0.006428

38.538218

B0 ¼ 10; B1 ¼ 1

26.045250

0.993981

33.698399

B0 ¼ 1; B1 ¼ 10

26.045267

0.973416

33.698407

0.441416 0.649832

48.879523 48.879449

λ2

5.354532 1.618745

Fig. 20. Longitudinal rolls: streamlines (solid), isotherms (dashed) and isosolutes (dotted) of dilatant fluid (n ¼ 2) for aiding buoyancy (N ¼ 1) taking Le ¼ 10;Pe ¼ 0:1; Ge ¼ 0:01.

instability, which depends on the value of buoyancy ratio N, Lewis number Le, P�eclet number Pe, power-law index n and Biot numbers B0 ; B1 . The variations of this critical value ωc is plotted against the power-

law index n are shown in Figs. 12 and 13 for aiding and opposing buoyancy, respectively. From these two figures, it can be observed that the critical value of ω for both pseudoplastic and dilatant fluids is 14

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 21. Longitudinal rolls: streamlines (solid), isotherms (dashed) and isosolutes (dotted) of dilatant fluid (n ¼ 2) for opposing buoyancy (N ¼ 10; Pe ¼ 0:1; Ge ¼ 0:01.

0:1) taking Le ¼

insignificant in aiding buoyancy (N) for varying values of Biot numbers. At higher value of B1 , the critical value of angular frequency ω decreased with increasing values of B0 (see Fig. 12(a)). The reverse behavior is seen in the case of higher values of B0 at the lower boundary wall, where the increasing value of B1 increased the value of ωc (see Fig. 12(b)). How­ ever, in both these cases, the critical value of ω increased with increasing value of the power-law index n in opposing buoyancy. In limiting case B0 →∞, B1 →∞ corresponding to isothermal bound­ ary walls, the variation of the critical value of ω with the power-law index n is presented in Kumari and Murthy [22]. The comparison with this study revealed that the value of ωc corresponding to isothermal boundary walls is higher than that with the convective boundary con­ ditions in opposing buoyancy. However, in the case of aiding buoyancy, at the low shear rate the critical value of ω is equal to zero. From Fig. 12, we also observe that at the large value of B0 and B1 , the value of ωc is moving towards the x-axis.

transverse rolls. Here, Le ¼ 1, Pe ¼ 10 are taken, to investigate the effect of high shear rate on the instability of base flow when viscous dissipation effect is introduced on the system. It can be concluded that at high shear rate has a destabilizing impact on the base flow for both longitudinal and transverse rolls, whereas at low shear rate, it has no significant effect on the instability of basic state solution. The reason is that, at the high shear rate, more heat is generated within the porous medium. Variation of the critical Rayleigh number Rac and the corresponding wave number ac for power-law fluids with Pe ¼ 10 is presented for the transverse rolls in Table 2 and for the longitudinal rolls in Table 3 for both aiding and opposing buoyancy. From Table 2 it is observed that, when Le ¼ 1, at high shear rate, Rac decreased with increasing values of Ge. This behavior is observed for transverse (Table 2) and longitudinal rolls (Table 3) in both pseudoplastic and dilatant fluids.

5.5. Effect of viscous dissipation

This limiting case corresponds to the isothermal top wall while the bottom wall is with constant wall flux conditions. In this case, the boundary conditions given in Eq. (14) becomes

10

6. Limiting case B0 ¼ 0; B1 →∞

Figs. 14 and 15 describe the impact of viscous dissipation (Ge ¼ 6 ;0:1) on the instability of power-law fluid for both longitudinal and 15

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

Fig. 22. Transverse rolls: streamlines (solid), isotherms (dashed) and isosolutes (dotted) drawn at t ¼ 0 of pseudoplastic fluid (n ¼ 0:5) for aiding buoyancy (N ¼ 1) taking Le ¼ 10; Pe ¼ 0:1; Ge ¼ 0:01.

z¼0:

∂T ¼ 1; ∂z

w ¼ 0;

z¼1:

w ¼ 0;

T ¼ 0;

C ¼ 1;

by convective boundary conditions B0 6¼ 0 and B1 6¼ 0. Variation of Rac against the power-law index n is shown in Figs. 16 and 17 for aiding and opposing buoyancy respectively. At the low shear rate, Rac is mono­ tonically decreasing with increasing values of power-law index n for longitudinal rolls. But in the case of transverse rolls, it increased first and then decreased. One can also note that at the high shear rate, Rac is monotonically increasing with n for longitudinal and transverse rolls in both aiding and opposing buoyancy situations. This feature shows that both at high and low shear rates, transverse rolls are the most unstable mode of instability for pseudoplastic fluids while the longitudinal rolls are so for the dilatant fluids. From these two figures, it can also be observed that, with constant heat flux bottom boundary condition, the oppositely directed thermal and solutal buoyancy forces (N < 0) stabilized the flow field. The opposite behavior is seen in the aiding buoyancy (N > 0). When thermal diffusivity dominates the solutal diffusivity (Le > 1), the value of Rac decreased with the increasing values of the power-law index n for aiding buoyancy case (N > 0) as seen from Fig. 16. For a fixed value of buoy­ ancy ratio N ¼ 0:1, the value of Rac increased with the increasing values of the diffusivity ratio Le. This phenomenon is seen at the low and high shear rate, is shown in Fig. 17. In this case the effect of viscous

(44)

C ¼ 0:

The basic temperature distribution is TB ¼ 1

z

Ge Penþ1 2 z 2 Ra

� 1 :

(45)

Nothing changes in the formulation of the system of coupled ordi­ nary differential equations, which is still given by Eqs.(36)–(38). On the other hand, as corresponding to Eq. (39), the boundary conditions are written as, z¼0: z¼1:

0

f ¼ 0; h ¼ 0; g ¼ 0; f ¼ 0; h ¼ 0 g ¼ 0:

(46)

This eigenvalue problem is solved using the bvp4c routine in MAT­ LAB R2014b. Results indicate that when n ¼ 1, Ge ¼ 0, and N ¼ 0, the present computed value of Rac ¼ 27:09762502 and ac ¼ 2:326235 co­ incides with that results given in Ref. [5]. It is also observed that the instability of base flow for constant heat flux boundary at the lower wall z ¼ 0 is an intermediate scenario as we compare the instability induced 16

International Journal of Thermal Sciences 146 (2019) 106098

S. Kumari and P.V.S.N. Murthy

Fig. 23. Transverse rolls: streamlines (solid), isotherms (dashed) and isosolutes (dotted) drawn at t ¼ 0 of pseudoplastic fluid (n ¼ 0:5) for opposing buoyancy (N ¼ 0:1) taking Le ¼ 10; Pe ¼ 0:1; Ge ¼ 0:01.

dissipation has been analyzed for N < 0 and N > 0 at Pe ¼ 10. The values of critical Rayleigh number Rac and the corresponding wave­ number ac have been tabulated for aiding and opposing buoyancy for transverse rolls in Table 4 and for longitudinal rolls in Table 5. Corre­ sponding to the high shear rate, the value of Rac decreased with increasing Ge for all values of n. In this case more heat generation due to viscous dissipation destabilized the basic state solution, which is seen with smaller values of Rac for both longitudinal and transverse rolls.

TB ¼ 1

z¼0: z¼1:

T ¼ 1; T ¼ 0;

C¼0 C ¼ 0:

� z2 :

(48)

f ¼ 0; h ¼ 0; g ¼ 0; f ¼ 0; h ¼ 0 g ¼ 0:

(49)

This eigenvalue problem is solved using the bvp4c routine in MAT­ LAB R2014b. Results indicate that when n ¼ 1, Ge ¼ 0, and N ¼ 0, the present computed value of Rac ¼ 39:478414 and ac ¼ 3:14160 coincides with the value given in Eq. (40) in Barletta and Nield. For the case of double-diffusive convective stability, the basic state solution is seen to be more stable with these boundary conditions when compared to that of third kind thermal boundary conditions. Variation of Rac against the power-law index n for selected values of Ge; Le; Pe and N is shown in Figs. 18 and 19 for aiding and opposing buoyancy respectively. In both aiding and opposing buoyancy, the curves show similar behavior, as explained in section 6. The effect of high shear rate and the viscous dissipation (Ge) on the critical value of Ra for both aiding and opposing buoyancy is displayed in Table 6 for transverse rolls

The other limiting case B0 →∞; B1 →∞ is considered to understand the behavior of disturbances on the basic steady state solution. In this case, the third kind thermal boundary conditions given in Eq. (14) become the isothermal conditions, given by w ¼ 0; w ¼ 0;

Ge Penþ1 z 2 Ra

The eigenvalue problem remains the same which is still given by Eqs. (36)–(38) with the boundary conditions,

7. Limiting case B0 →∞; B1 →∞

z ¼ 0; z ¼ H;

zþ

(47)

The basic temperature distribution is given as

17

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

and in Table 7 for longitudinal rolls. For Pe ¼ 10 and Le ¼ 1, it can be observed that the values of Rac decreased with increasing Ge for both longitudinal and transverse rolls when N > 0 and N < 0. In other words, heat generation due to viscous dissipation has a destabilizing effect when the shear rate becomes high.

SðzÞ ¼

From Eq. (16) we also obtain the basic temperature profile as Ge Pe2 ½

2 þ B0 ð 2 þ zÞz B1 ð 1 þ zÞð1 þ z þ B0 zÞ� : 2ðB0 þ B1 þ B0 B1 Þλ

(50)

Substituting the expansions of stream function, temperature and concentration profiles in governing Eqs.(36)–(39), we obtain sets of eigenvalue problems for constant, first and second order-approximation with λ0 , λ1 and λ2 as the components of eigenvalue. The constant order O (1) problem is given by n f 00 0 n sin2 α þ cos2 α

(51)

a2 f0 þ aλ0 ðh0 þ Ng0 Þ ¼ 0;

h00 0

a2 h0

iωh0 þ af0 ¼ 0;

(52)

g00 0

a2 g0

iωLeg0 þ aLef0 ¼ 0;

(53)

z¼0: z¼1:

f0 ¼ 0; f0 ¼ 0;

0

h 0 þ B0 h0 ¼ 0; g0 ¼ 0; 0 h 0 þ B1 h0 ¼ 0; g0 ¼ 0:

From the above analysis, it is evident that the critical values of the thermal Rayleigh number (eigenvalue) and the wave number are ob­ tained using the bvp4c in MatlabR2014b. But the in-built functions NDSolve and FindRoot in Mathematica are more amicable to generate the data for presenting the pattern for the streamlines, isotherms and isosolutes at these critical values. These routines make use of the sixth order explicit Runge-Kutta method and shooting method, respectively. The detailed description of this numerical procedure for longitudinal and transverse rolls is discussed in Barletta and Nield [16]. This nu­ merical method is used in this investigation for graphical representation of streamlines, isotherms and isosolutes. Longitudinal rolls are the favorable mode of instability for Newto­ nian fluids. However, in non-Newtonian fluids, both longitudinal and transverse rolls are equally important. The streamlines, isotherms and isosolutes for longitudinal and transverse rolls with different pairs of Biot numbers are shown in Figs. 20–23 for different power-law index n. The solid, dashed and dotted lines represent streamlines, isotherms and isosolutes respectively. Figs. 20 and 21 show the streamlines, isotherms and isosolutes corresponding to longitudinal rolls of a dilatant fluid with n ¼ 2 in the aiding (N ¼ 1) and opposing buoyancies (N ¼ 0:1). These rolls are chosen because they are the most preferred mode of instability for dilatant fluids. For pseudoplastic fluids (n ¼ 0:5) these streamline, isotherm and isosolute pattern is displayed for transverse rolls in Figs. 22 and 23 for aiding and opposing buoyancies respectively. When heat conduction inside the medium is faster than the heat convection at both the boundary walls (small value of B0 and B1 ), the streamlines are seen to be more concentrated in the center in the aiding buoyancy case. With higher values of B0 and B1 , the streamlines patter is seen as scattered in the whole region and these results are shown for Pe ¼ 0:1. With higher values of Pe, not much significant changes are seen in the pattern of the streamlines, isotherms and isosolutes for both

(54)

λ0 h00 1 þ g00 1

a2 h1

a2 f1 þ aλ1 ðh0 þ N g0 Þ þ aλ0 ðh1 þ N g1 Þ ¼ 0;

� iωh1 þ af1 þ λ1 h00 0

a2 h0

iωh0 þ af0

a2 g1

(55)

�

iGeðn þ 1Þcos2 α 0 f 0 ¼ 0; sin2 α þ cos2 α

(56) (57)

iωLeg1 þ aLef1 ¼ 0; 0

z ¼ 0 : f1 ¼ 0; h 1 þ B0 h1 ¼ 0; g1 ¼ 0; 0 z ¼ 1 : f1 ¼ 0; h 1 þ B1 h1 ¼ 0; g1 ¼ 0;

(58)

Similarly, the second order system for OðPe2 Þ, is n f 00 2 n sin2 α þ cos2 α þ N g2 Þ ¼ 0;

a2 f2 þ aλ2 ðh0 þ N g0 Þ þ aλ1 ðh1 þ N g1 Þ þ aλ0 ðh2 (59) �

�

a2 h2 iωh2 þ af2 þ λ1 h00 1 a2 h1 iωh1 þ af1 þ λ2 h00 0 � iGeðn þ 1Þcos2 α 0 iωh0 þ af0 þ aGef0 SðzÞ þ f 1 ¼ 0; n sin2 α þ cos2 α λ0 h00 2

g00 2

a2 g2

iωLeg2 þ aLef2 ¼ 0;

ðB0ð 2 þ zÞ þ ð1 þ B0ÞB1ð 1 þ zÞ þ B0 z þ B1ð1 þ z þ B0 zÞÞ : 2ðB0 þ B1 þ B0 B1Þ

9. Streamlines, isotherms and isosolutes

The system corresponding to the first order approximation OðPeÞ is obtained as n f 00 1 n sin2 α þ cos2 α

(62)

The values of λ0 , λ1 and λ2 along with the corresponding eigenfunc­ tions are obtained by solving these system of equations simultaneously 0 0 along with the normalization conditions f 0 ð0Þ ¼ 1, f 1 ð0Þ ¼ 0 and 0 f 2 ð0Þ ¼ 0 using the bvp4c in MATLAB R2014b. The relative and abso­ lute tolerance are set as 10 3 and 10 6 , respectively. For the case of longitudinal rolls, the value of λ1 becomes zero in both aiding and opposing buoyancies, which is consistent with the results presented in Barletta and Storesletten [6]. The computed values of λ0 and λ2 are shown in Table 8 and Table 9 for aiding and opposing buoyancy respectively. These values are seen to be independent of power-law index n. The non-zero value of λ2 has no significant effect on the value of λ at Pe→0. In case of transverse rolls, values of λ0 and λ2 are shown in Table 10 for aiding buoyancy and in Table 11 for the opposing buoyancy. In this case also, the values of λ1 and ω are equal to zero. And the eigenvalue λ depends on the viscous dissipation parameter Ge and the power-law index n also. The second-order approximate value of λ2 is not equal to zero for transverse rolls due to the presence of the viscous dissipation in the medium. The numerical computations reveal that when there is no viscous dissipation in the medium, the values of λ1 and λ2 become zero for both longitudinal and transverse rolls. For the Newtonian fluid (n ¼ 1) case, the value of λ0 is equal for both longitu­ dinal and transverse rolls. We also noticed that, as Pe→0, the value of Ra→∞ for the pseudoplastic fluids and Ra→0 for the dilatant fluids. These results are consistent with the results presented in the literature.

The onset of the double-diffusive convective instability of power-law fluid at the singular point Pe ¼ 0 has been discussed with the convective thermal boundary conditions (B0 , and B1 6¼ 0) in this section. For the asymptotic analysis, we consider a second order regular perturbation expansion for the functions f, h, g and the parameter λ ¼ PejRan 1 in Pe as j � f ðzÞ ¼ f0 ðzÞ þ f1 ðzÞPe þ f2 ðzÞPe2 þ O Pe3 ;� hðzÞ ¼ h0 ðzÞ þ h1 ðzÞPe þ h2 ðzÞPe2 þ O Pe3 �; 2 3 gðzÞ ¼ g0 ðzÞ þ g1 ðzÞPe þ g2 ðzÞPe � þ O Pe ; 2 3 λ ¼ λ0 þ λ1 Pe þ λ2 Pe þ O Pe :

z

0

h 2 þ B0 h2 ¼ 0; g2 ¼ 0; 0 h 2 þ B1 h2 ¼ 0; g2 ¼ 0:

where,

8. An asymptotic solution for Pe→0

TB ¼ 1

f2 ¼ 0; f2 ¼ 0;

z¼0: z¼1:

a2 h0

(60) (61)

18

S. Kumari and P.V.S.N. Murthy

International Journal of Thermal Sciences 146 (2019) 106098

pseudoplastic and dilatant fluids. In case of opposing buoyancy (N ¼ 0:1), with B0 ¼ 0:1; B1 ¼ 100, streamlines are seen to be pushed downwards. On interchanging the value of B0 ¼ 100 and B1 ¼ 0:1, the streamlines are pushed upward because the boundary effects at the lower and upper walls change the streamlines patterns inside the medium. The combined effect of conduction and convection on the isotherms pattern is shown in Fig. 20 for aiding buoyancy and in Fig. 21 for opposing buoyancy for n ¼ 2. When conduction dominates the effect of convection (B0 ¼ 0:1) at the lower wall, and convection dominates over the conduction (B1 ¼ 100) at the upper wall z ¼ H, the isotherm pat­ terns are mutually opposite in nature. It is so because the basic tem­ perature profile depends on both these Biot numbers B0 and B1 . The isotherm pattern for transverse rolls in pseudoplastic fluid (n ¼ 0:5), is similar to that of the longitudinal rolls in dilatant fluid in both aiding and opposing buoyancy cases. But, we did not observe significant vari­ ations in the isosolute pattern in both aiding and opposing buoyancy cases. However, at B0 ¼ 0:1 and B1 ¼ 100, these isosoutes are concen­ trated in the center and, at B0 ¼ 100 and B1 ¼ 0:1, these are seen as scattered in the entire region.

longitudinal rolls are stabilized by increasing the value of P�eclet number, when both B0 and B1 are non-zero. � The asymptotic analysis for Pe→0, the instability of dilatant fluids are described by vanishing values of apparent viscosity and stability of pseudoplastic fluids are defined by infinite apparent viscosity. References [1] R.V. Dharmadhikari, D.D. Kale, Flow of non-Newtonian fluids through porous media, Chem. Eng. Sci. 40 (3) (1985) 527–528. [2] L.S.de B. Alves, A. Barletta, S. Hirata, M.N. Ouarzazi, Effects of viscous dissipation on the convective instability of viscoelastic mixed convection flows in porous media, Int. J. Heat Mass Transf. 70 (2014) 586–598. [3] S. Wang, W. Tan, Stability analysis of Soret-driven double-diffusive convection of Maxwell fluid in a porous medium, Int. J. Heat Fluid Flow 32 (1) (2011) 88–94. [4] A. Barletta, D.A. Nield, Linear instability of the horizontal throughflow in a plane porous layer saturated by a power-law fluid, Phys. Fluids 23 (1) (2011), 013102. [5] L.S.de B. Alves, A. Barletta, Convective instability of the Darcy–B� enard problem with through flow in a porous layer saturated by a power-law fluid, Int. J. Heat Mass Transf. 62 (2013) 495–506. [6] A. Barletta, L. Storesletten, Linear instability of the vertical throughflow in a horizontal porous layer saturated by a power-law fluid, Int. J. Heat Mass Transf. 99 (2016) 293–302. [7] L.S.de B. Alves, A. Barletta, Convective to absolute instability transition in the Prats flow of a power-law fluid, Int. J. Therm. Sci. 94 (2015) 270–282. [8] D.A. Nield, Onset of thermohaline convection in a porous medium, Water Resour. Res. 4 (3) (1968) 553–560. [9] D.A. Nield, D.M. Manole, J.L. Lage, Convection induced by inclined thermal and solutal gradients in a shallow horizontal layer of a porous medium, J. Fluid Mech. 257 (1993) 559–574. [10] F. Chen, C.F. Chen, Double-diffusive fingering convection in a porous medium, Int. J. Heat Mass Transf. 36 (3) (1993) 793–807. [11] S. Wang, W. Tan, Stability analysis of double-diffusive convection of Maxwell fluid in a porous medium heated from below, Phys. Lett. A 372 (17) (2008) 3046–3050. [12] A. Barletta, E.R. di Schio, L. Storesletten, Convective roll instabilities of vertical throughflow with viscous dissipation in a horizontal porous layer, Transp. Porous Media 81 (3) (2010) 461–477. [13] A. Barletta, D.A. Nield, Instability of Hadley–Prats flow with viscous heating in a horizontal porous layer, Transp. Porous Media 84 (2) (2010) 241–256. [14] A. Barletta, M. Celli, D.A.S. Rees, The onset of convection in a porous layer induced by viscous dissipation: a linear stability analysis, Int. J. Heat Mass Transf. 52 (1–2) (2009) 337–344. [15] M. Celli, A. Barletta, Onset of convection in a non-Newtonian viscous flow through a horizontal porous channel, Int. J. Heat Mass Transf. 117 (2018) 1322–1330. [16] A. Barletta, D.A. Nield, Thermosolutal convective instability and viscous dissipation effect in a fluid-saturated porous medium, Int. J. Heat Mass Transf. 54 (7–8) (2011) 1641–1648. [17] R.R. Kairi, P.V.S.N. Murthy, Effect of melting on mixed convection heat and mass transfer in a non-Newtonian fluid saturated non-Darcy porous medium, J. Heat Transf. 134 (4) (2012), 042601. [18] A.V. Shenoy, Non-Newtonian fluid heat transfer in porous media, Adv. Heat Tran. 24 (1994) 102–191. [19] R.H. Christopher, S. Middleman, Power-law flow through a packed tube, Ind. Eng. Chem. Fundam. 4 (4) (1965) 422–426. [20] M.F. El-Amin, M.A. El-Hakiem, M.A. Mansour, Effects of viscous dissipation on a power-law fluid over plate embedded in a porous medium, Heat Mass Transf. 39 (10) (2003) 807–813. [21] M.-C. Charrier-Mojtabi, B. Elhajjar, A. Mojtabi, Analytical and numerical stability analysis of Soret-driven convection in a horizontal porous layer, Phys. Fluids 19 (12) (2007) 124104. [22] S. Kumari, P.V.S.N. Murthy, Stability of the Horizontal Throughflow in a PowerLaw Fluid Saturated Porous Layer, Transport in Porous Media, 2019, pp. 1–20.

10. Conclusions The double-diffusive convective instability of horizontal through­ flow of the power-law fluid in a horizontal porous channel has been investigated by considering viscous dissipation in the medium and with third kind thermal boundary conditions. Significant observations are listed below: � In both aiding N > 0 and opposing N < 0 buoyancy cases, longitu­ dinal rolls are the most unstable mode of disturbance for dilatant fluid (n > 1), whereas for pseudoplastic fluid (n < 1) the transverse rolls are the most unstable mode of disturbance. � When convection dominates the heat conduction inside the porous channel, it enhanced the stability of power-law fluids. As the in­ tensity of diffusivity ratio is increased, the fluid destabilizing char­ acter increased in the aiding buoyancy. The opposite nature is seen for the opposing buoyancy cases. � The basic state solution for the present physical configuration is more stable with isothermal boundary conditions (B0 →∞; B1 → ∞) compared to that with the general convective boundary conditions. The limiting case corresponding to the constant heat flux at the bottom wall (B0 ¼ 0; B1 →∞), is an intermediate scenario to the above mentioned two cases with respect to the stability of the basic state solution. � When more heat generates inside the system due to viscous dissi­ pation effect, the destabilizing nature of power-law fluid with the convective thermal boundary case is higher compared to the isothermal boundary conditions. � For pseudoplastic fluids, the transverse rolls are destabilized by increasing the P�eclet number. And for the dilatant fluids,

19