1 Introduction

The buoyancy-driven flow instability in a layer of fluid-saturated porous media has received increasing attention over the past fifty years because of its numerous applications in geophysics and energy-related systems^[1-5]. The porous media involved in many engineering and industrial applications are usually anisotropic in their mechanical and thermal properties, and several investigators have addressed these effects on the stability characteristics of the system. Castinel and Combarnous^[6] were the first to investigate the effect of the anisotropy in the permeability of a porous medium on the natural convection in the porous medium. Epherre^[7] extended this work to the case of anisotropy in both thermal diffusivity and permeability. Subsequently, several studies were undertaken, covering various effects on this problem, and the results have been reported in the open literature^[8-11].

Most of the investigations on the thermal convection in anisotropic porous media are for Newtonian fluids, and seldom for non-Newtonian fluids. However, the thermal convection encountered in many engineering applications, e.g., geophysics, material processing, petroleum, chemical and nuclear industries, reservoir engineering, and bioengineering, exhibits non-Newtonian characteristics. This has encouraged researchers to consider non-Newtonian fluids in their investigations. Among different kinds of non-Newtonian fluids, viscoelastic fluids are found to be of considerable importance in various engineering applications^[12]. Alishaev and Mirzadjanzade^[13] investigated the viscoelastic flows in porous media for the calculations of delay phenomenon in the filtration theory. Rudraiah and Kaloni^[14] provided a review on some of the constitutive equations of non-Newtonian fluids flowing through porous media. By applying the linear stability theory, Rudraiah et al.^[15] examined the onset of convection in an Oldroyd-B fluid-saturated horizontal porous layer heated from below. Shenoy^[16] gave a comprehensive review on the non-Newtonian fluids and heat transfer in porous media. Kim et al.^[17] studied the thermal instability in a porous layer saturated with a viscoelastic fluid. Malashetty et al.^[18] and Shivakumara et al.^[19] analyzed the effects of local thermal non-equilibrium on the convection onset in a viscoelastic fluid-saturated porous layer. Sheu et al.^[20] investigated the buoyancy-induced convection in a viscoelastic fluid-saturated porous medium. Wang and Tan^[21] used both linear and nonlinear stability theories to investigate double diffusive convection with the modified Darcy-Maxwell model with the Soret effect. Recently, various types of flow problems for the Oldroyd-B and Maxwell viscoelastic fluids have been analyzed^[22-24]. Considering another non-Newtonian fluid known as the Casson fluid, Makinde and Eegunjobi^[25] discussed the thermally radiating magnetohydrodynamics slip flow in a micro channel filled with Casson fluid-saturated porous media, and Makinde and Rundora^[26] analyzed the unsteady mixed convection flow in a permeable wall channel with a reactive Casson fluid-saturated porous medium.

The investigations on the thermal convection in a viscoelastic fluid-saturated porous medium are mainly dispensed with isotropic porous media except the study of Malashetty and Swamy^[12], where the convection onset in a layer of viscoelastic liquid-saturated anisotropic Darcy porous media was investigated. Since the stability of the viscoelastic fluid-saturated anisotropic porous layer takes the form of overstable motion only, it is of interest to perform the nonlinear stability analysis and quantify the role of anisotropy and viscoelastic parameters on the same layer. The goal of the present paper is to investigate the nonlinear stability of thermal convection in an Oldroyd-B fluid-saturated anisotropic porous layer with the perturbation method. The stability of bifurcating equilibrium solutions is discussed by deriving cubic Landau equations. The results of the linear instability analysis are delineated. Besides, the consequence of the viscoelastic and anisotropy parameters on the variation of the Nusselt number with respect to the Darcy-Rayleigh number is examined.

2 Mathematical formulation

The physical configuration is as shown in Fig. 1. We consider an infinite horizontal anisotropic porous layer heated from below and saturated with a viscoelastic fluid of an Oldroyd-B type confined between the impermeable planes z=0 and z=d in the presence of gravity. The anisotropy in both thermal diffusivity and permeability is considered. A Cartesian coordinate system (x, y, z) is selected such that the origin is located at the lower boundary and the z-axis is measured vertically upward. The lower and upper boundaries are maintained at T₀ +ΔT (Δ T>0) and T₀, respectively.

The basic governing equations under the Boussinesq approximation are^[13-17]

(1)

(2)

(3)

(4)

where q=(u, v, w) is the velocity vector, p is the pressure, µ is the fluid viscosity, λ₁ is the stress relaxation time, λ₂ is the strain retardation time, g is the gravitational acceleration, ρ is the fluid density, T is the temperature, β is the thermal expansion coefficient, ρ₀ is the reference density at T=T₀, and

(5)

are the effective thermal diffusivity and the inverse permeability tensors, respectively, whose principal axes are associated with the coordinate system. A horizontal isotropy in the permeability and thermal diffusivity is assumed and considered, i.e.,

where k_h and α_h are the permeability and the thermal diffusivity in the horizontal and directions, respectively, while k_v and α_v are the corresponding values in the vertical direction, respectively.

The steady basic state is quiescent, and is given by

(6)

To study the stability of the basic state, the perturbations on the basic state are superimposed as follows:

(7)

where primes designate the perturbed quantities. Substituting Eq.(7) into Eqs.(1)-(4), eliminating the pressure term from Eq.(2) by operating curl, introducing the stream function ψ(x, z, t) through

(8)

and rendering the resulting equations to dimensionless form by using d, d² /α_v, α_v, and ΔT as the units of the length, time, stream function, and temperature, respectively, we have

(9)

where

(10)

and J(., .) stands for the Jacobian with respect to x and z. Here, R_D is the Darcy-Rayleigh number defined by

Λ₁ is the relaxation parameter defined by Λ₁ =λ₁α_v/d², Λ₂ is the retardation parameter defined by Λ₂ =λ₂α_v /d², η is the thermal anisotropy parameter defined by η =α_h /α_v, and ξ is the mechanical anisotropy parameter defined by ξ =k_h /k_v.

Since the boundaries are impermeable and isothermal, the appropriate boundary conditions are

(11)

3 Nonlinear stability analysis

The nonlinear stability analysis near the convection threshold is performed with the perturbation method. The cubic Landau equations are derived for stationary and oscillatory convection modes. Such a study helps in analyzing the stability of the bifurcating equilibrium solutions (subcritical/supercritical) and also in estimating the convective rate of heat transfer. Accordingly, the dependent variables ψ, T, and R_D are expanded in the power series of a small perturbation parameter χ(« 1) as follows^[27-28]:

(12)

The other parameters Λ₁, Λ₂, ξ, and η are taken as given, and R_Dc is the critical value at the threshold as the case may be. A slow time scale s is introduced, i.e., s=χ²t. The operator is replaced, depending on the nature of the bifurcating solutions.

3.1 Bifurcation of the stationary solution

In this case, R_D =R_Dc^s. Substituting Eq.(12) and into Eq.(9) and equating the coefficients of different powers of χ lead to a sequence of equations.

For the first-order power of χ, the resulting stability equations are homogeneous, and are

(13)

where

The eigenvalue and eigenfunctions of this problem are

(14)

(15)

The undetermined amplitudes A₁ and B₁ are related by

(16)

where c=ηa²+π². The eigenvalue R_D^s attains the critical value at , and the corresponding critical value is

(17)

which is free from the viscoelastic parameters and coincides with the Newtonian case^[9]. For the isotropic case (ξ =η), a_c =π, and R_Dc^s =4π², which are the known exact values^[2].

For the second-order power of χ, the stability equations are inhomogeneous, and are

(18)

where . The solution of the above system of equations is

(19)

For the third-order power of χ, the stability equations become

(20)

where

The solution of the above equations is

(21)

The solvability condition has been derived for Eq.(20), which is in the form of the first-order nonlinear ordinary differential equation (the cubic Landau equation) for the unknown amplitude B₁ as follows:

(22)

where

(23)

(24)

For the steady state, the amplitude is found to be

(25)

Equation (25) is independent of viscoelastic parameters while depends on anisotropy parameters. Note that R_D2 >0, which indicates that the stationary bifurcation is always supercritical (stable).

The convective heat transfer is determined in terms of the area-averaged thermal Nusselt number. The Nusselt number is defined by

(26)

By substituting the critical values into the above equation, we have

(27)

If η =ξ, Eq.(27) coincides with the Newtonian case^[29].

3.2 Bifurcation of the oscillatory solution

In this case, R_D =R_Dc^o. A minor alteration of the method applied in the earlier section is used to find out the bifurcation of the oscillatory convection. The time derivative is not zero in the present case, and is replaced by .

For the first-order power of χ, the equations reduce to the linear instability problem for overstability as follows:

(28)

The eigenfunctions are

(29)

where the overline denotes the complex conjugate, ω and a are taken to be the critical conditions associated with the oscillatory onset. The amplitudes A₁ and B₁ are functions of the slow time scale, and are related by

(30)

The eigenvalue is found to be

(31)

where

(32)

Equation (32) shows that the oscillatory convection is not possible if Λ₁ < Λ₂. It is seen that R_D^o attains its critical value R_Dc^o at a²=a_c², where

(33)

which is independent of Λ₁, and the corresponding expression for R_Dc^o is

(34)

For the second-order power of χ, the equations are inhomogeneous and found to be

(35)

where

(36)

Equation (36) suggests that the stream function and temperature should contain the terms involving the frequency 2ω. Based on this fact, the second-order stream function and temperature can be expressed as follows:

(37)

The solution of the second-order problem is now found to be

(38)

For the third-order power of χ, the stability equations are

(39)

where

The third-order problem has the solution as follows:

(40)

Equation (39) gives the following cubic Landau equation that explains the temporal variation of B₁ of the convection cell:

(41)

where

(42)

(43)

From Eq.(42), the following relation can be obtained:

(44)

(45)

where

and ph(.) represents the phase shift. The temporal evolution of |B₁| can be expressed as a function of the initial amplitude B₀ as follows:

(46)

From the above equation, it is seen that |B₁|~ B₀ exp (p_rs) as s→ -∞ and |B₁|→ 0, just as the linear theory, but as s→∞, which is independent of the value of B₀. For the post-transient state, Eq.(45) yields an expression for the amplitude as follows:

(47)

If Ω >0, the bifurcation is supercritical. If Ω < 0, the bifurcation is subcritical. This can be achieved by evaluating the expression for Ω for various values of the physical parameters since there is no simple way to analyze this expression. For this case, the area and time-averaged thermal Nusselt number can be represented by using Eq.(11) as follows:

(48)

With Eqs.(37) and (38), we can rewrite Eq.(48) as follows:

(49)

4 Results and discussion

The effects of the anisotropy in permeability and thermal diffusivity on the nonlinear stability of thermal convection in a horizontal porous layer saturated by an Oldroyd-B fluid are investigated. Since the considered nonlinear stability analysis is based on the linear instability analysis, the results of the linear instability theory are also discussed. Although the stationary convection boundary depends on anisotropy parameters while is independent of viscoelastic parameters, it concurs with the Newtonian fluid-saturated anisotropic porous layer when the base flow is quiescent. The oscillatory convection boundary, however, depends on viscoelastic parameters, e.g., mechanical and thermal anisotropy parameters.

The neutral stability curves on the (a, R_D)-plane for different values of the stress relaxation parameter Λ₁, strain retardation parameter Λ₂, mechanical anisotropy parameter ξ, and thermal anisotropy parameter η are presented in Figs. 2 and 3. The region underneath the neutral curve corresponds to the stability region, above which it is unstable. It is observed that the effects of increasing Λ₁ (see Fig. 2(a)) and ξ (see Fig. 3(a)) as well as decreasing Λ₂ (see Fig. 2(b)) are to decrease the stability region, while the effect of increasing η (see Fig. 3(b)) is to increase the stability region. Besides, the oscillatory neutral stability curves shift towards lower values of the wavenumber when Λ₁ and ξ increase, which indicates that the cell width at the critical state increases while Λ₂ and η decrease.

The critical Darcy-Rayleigh number R_Dc and the corresponding critical oscillation frequency ω_c are obtained for various physical parameters. The results are presented in Figs. 4-6.

Figures 4(a) and 4(b), respectively, show the variations of R_Dc and ω_c as functions of Λ₂ for different values of Λ₁, where η/ξ =0.5. It is noted that increasing Λ₁ is to advance the onset of oscillatory convection for any fixed value of Λ₂. This may be attributed to the fact that increasing the relaxation parameter ceases the stickiness of the viscoelastic fluid and hence the effect of friction will be reduced so that the convection sets in at lower values of R_Dc. On the contrary, increasing Λ₂ delays the onset of oscillatory convection for a fixed value of Λ₁ because increasing Λ₂ amounts to increasing the time taken by the fluid element to respond to the applied stress. Further inspection of Fig. 4(a) reveals that the range of the values of Λ₂, within which the oscillatory convection possibly increases with increasing Λ₁. In other words, for a fixed value of Λ₁, there exists a threshold value Λ₂^* which divides the boundary of regimes between the oscillatory and stationary convection. Initially, convection begins in the form of an oscillatory mode. As the value of Λ₂ reaches Λ₂^*, convection ceases to be oscillatory and stationary convection becomes the preferred mode of instability. The value of Λ₂^* depends on other physical parameters as well. The critical frequency ω_c shown in Fig. 4(b) exhibits that it decreases with increasing Λ₂ and increases with increasing Λ₁ due to the increase in the elasticity of the fluid.

The variations of (R_Dc and ω_c) as functions of η/ξ are shown in Figs. 5 and 6 for different values of Λ₁ (with Λ₂ =0.1) and Λ₂ (with Λ₁ =0.5), respectively. These figures clearly indicate that the effects of increasing η/ξ is to delay the onset of convection and to increase the frequency of oscillations. The increase in η/ξ amounts to either decreasing ξ or increasing η. We note that the decrease in ξ amounts to decreasing the horizontal permeability, which impedes the motion of the fluid in the parallel direction. As a result, the transfer process in the porous medium gets suppressed, and hence higher values of R_Dc are needed for the onset of instability. Moreover, the increase in η amounts to increasing the horizontal thermal diffusivity. Thus, heat can be transported with ease in the porous layer, the horizontal temperature differences in the fluid, which are necessary to maintain convection, are more competently dissipated with increasing η, and higher values of R_Dc are required for the onset of convection. It is intriguing to note that by altering the anisotropy in the permeability and thermal diffusivity, it is possible to control the (augment/suppress) convective instability. From the figures, it is further evident that the increase in Λ₁ and decrease in Λ₂ are to hasten the onset of oscillatory convection. Moreover, the critical frequency increases with increasing Λ₁ while decreases with increasing Λ₂.

The parameters for the boundary separating stationary and oscillatory solutions are estimated. Figures 7(a) and 7(b) show the bifurcation of the stationary and oscillatory solutions in the viscoelastic parameter plane for different values of η and ξ, respectively. The region above each curve corresponds to the system which is unstable under oscillatory convection, and the region below the curve corresponds to the system which is unstable under stationary convection. From Figs. 7(a) and 7(b), for a fixed value of Λ₂/Λ₁ , it is seen that the value of Λ₁, at which codimension-two bifurcation occurs, decreases when η increases (see Fig. 7(a)), while an opposite trend is observed when ξ increases (see Fig. 7(b)). As the value of Λ₂/Λ₁ advances towards 0.9, there is a steep rise in the value of Λ₁.

The stability of the stationary and oscillatory bifurcating solutions is analyzed by deriving the cubic Landau equations for these cases. It is an observable fact that the stationary solution always bifurcates supercritically (see Eq.(25)), while the stability of the oscillatory solution can be understood from the sign of Ω. When Ω >0, the bifurcation is supercritical. When Ω < 0, the bifurcation is subcritical (see Eq.(47)). Hence, the expression Ω is evaluated for a wide range of parametric values at the critical values of oscillatory convection, and is denoted by Ω_c. Figures 8(a) and 8(b) represent the computed values of Ω_c as a function of η/ξ for different values of Λ₁ and Λ₂, respectively. These figures indicate that the oscillatory solution always bifurcates supercritically. Thus, the linear instability analysis provides the necessary and sufficient conditions for instability.

The heat transfer is estimated in terms of the Nusselt number for both stationary and oscillatory cases. For the stationary case, the area-averaged Nusselt number is calculated as a function of R_D for different values of η/ξ. It is seen that the Nusselt number originates from higher values of R_D with increasing η/ξ. The value of the Nusselt number increases with increasing R_D for any fixed value of η/ξ, and the heat transfer decreases with increasing η/ξ (see Fig. 9).

For the oscillatory case, the area and time-averaged Nusselt number Nu is calculated as a function of R_D for various values of the physical parameters. The variations of Nu as a function of R_D for different values of η/ξ as well as Λ₂ (with Λ₁ =1) and Λ₁ (with Λ₂ =0.3) are illustrated in Figs. 10(a) and 10(b), respectively. It is noted that the value of Nu increases when R_D increases for a fixed value of η/ξ. Moreover, the heat transfer increases with increasing Λ₁ while decreases with increasing η/ξ and Λ₂.

5 Conclusions

The nonlinear stability of thermal convection in an Oldroyd-B fluid-saturated anisotropic porous layer is investigated with the perturbation method. The onset of stationary and oscillatory convection is delineated since the nonlinear stability analysis is based on the results of the linear instability analysis. The instability sets in via the oscillatory mode under certain conditions, and the effect of increasing the mechanical and thermal anisotropy parameters is to advance and delay the onset of the oscillatory convection, respectively. A codimension-two bifurcation occurs at well-defined parametric conditions, and the value of the relaxation parameter, at which it occurs, decreases with increasing the thermal anisotropy parameter and decreasing the mechanical anisotropy parameter in the viscoelastic parameter plane. The stability of the stationary and oscillatory cases is discussed by deriving the cubic Landau equations. It is observed that these solutions always bifurcate supercritically. The increases in the value of the relaxation and retardation parameters are to enhance and suppress the time and area-averaged heat transfer, respectively. Besides, the increase in the ratio of the thermal anisotropy parameter to the mechanical anisotropy parameter is to decrease the heat transfer. By tuning the anisotropy of the porous medium, it is possible to control the convective instability of the system.

Acknowledgements The authors wish to thank the reviewers for their useful suggestions which are helpful in improving the paper significantly. One of the authors (K. R. Raghunatha) (SRF) wishes to thank the Department of Science and Technology, New Delhi for granting him a fellowship under the Innovation in Science Pursuit for the Inspired Research (INSPIRE) Program (No. DST/INSPIRE Fellowship/[IF 150253]).