1 Introduction

The thermal convective instability in a layer of Newtonian fluid-saturated porous media is a classical problem in convective heat transfer, and now it is a well-understood phenomenon. Different effects have been taken into account on this study, and the developments are well documented in the literature^[1-4]. The non-Newtonian fluid flows and heat transfer in porous media have been a topic of interest in recent years due to their importance in geophysics, ceramic processing, bioengineering, filtration, liquid composite molding, polymer engineering, and oil reservoir engineering^[5-6]. To a large extent, the performance of an oil reservoir depends on the physical nature of the crude oil presented in the reservoir. The light crude oil is essentially Newtonian, but the heavy crude oil is found to exhibit a non-Newtonian fluid behavior. In particular, some oil sand contains waxy crude at shallow depth in reservoirs, which is considered to be a viscoelastic fluid^[7-8]. The Darcy-Bénard problem for viscoelastic fluids has received considerable attention in the literature because the study can give necessary information about the mobility control in the oil displacement mechanism, which can improve the efficiency of oil recovery^[9-17]. It has been revealed that the convection onset in a viscoelastic fluid-saturated porous medium is oscillatory instead of stationary, which depends on the fluid elasticity. Recently, various types of flow problems for the Oldroyd-B and Maxwell viscoelastic fluids have been analyzed^[18-20].

The double diffusive convection in porous media is another interesting topic of research, and has attracted many researchers from various fields of science and engineering due to its wide applications, e.g., the solidification of binary mixtures, the migration of the solutes in water-saturated soils, and the geophysical system, electro-chemistry, and the migration of moisture through the air contained in fibrous insulation. The study reveals a variety of interesting convective phenomena, which are not observed in a single component fluid-saturated porous medium. For example, the exchange principle of stabilities is not always valid, and the conductive base state may become unstable through a growing oscillatory mode. Newtonian fluids have been widely studied on this topic^[21-22]. However, the double diffusive convection in a non-Newtonian fluid-saturated porous medium has received only limited attention in the literature^[23-25].

Nonetheless, many fluid dynamical systems occurring in nature and engineering applications involve more than two stratifying agencies with different molecular diffusivities, in which a multicomponent convection-diffusion is bound to occur. As a first step towards the understanding of the complex multicomponent convection, Griffiths^[26], Pearlstein et al.^[27], Straughan and Walker^[28], and Straughan and Tracey^[29] studied the triple diffusive convection in a fluid layer both experimentally and theoretically. Rudraiah and Vortmeyer^[30], Poulikakos^[31], and Griffiths^[26] studied the counterpart in a layer of porous media. Since the results obtained by these authors were incomplete in the vein of Pearlstein et al.^[27], Tracey^[32] reconsidered the problem, presented systematically the results similar to those of Pearlstein et al.^[27], and performed the nonlinear stability analysis with the energy method. Rionero^[33] considered a triply convective diffusive fluid mixture, which saturated a horizontal porous layer. Rionero^[34] considered the triple diffusive convection in porous media, obtained the conditions for inhibiting the convective onset, and performed a global nonlinear stability analysis. Ghalambaz et al.^[35] theoretically studied the triple diffusive convection in a square porous cavity.

The studies on the triple diffusive convection in porous media are mainly focused on Newtonian fluids, and a corresponding problem for non-Newtonian fluids is in the much-to-be desired state. Although Zhao et al.^[36] investigated the triply diffusive convection in a Maxwell fluid-saturated porous layer and obtained the criterion for the onset of stationary and oscillatory convection, their study is silent in revealing some of the unusual behaviors that the system is capable of supporting, which have important implications on the instability of the system. The interest to the present work is of two folds. One is to study the linear instability theory in detail and unveil some remarkable departures from those of single and double diffusive cases. The other is to perform a weakly nonlinear stability analysis to understand the stability of bifurcating the equilibrium solution by deriving a cubic Landau equation. The heat and mass transfer is also discussed in terms of the Nusselt numbers.

2 Mathematical formulation

The physical configuration is shown in Fig. 1. We consider a horizontal layer of an incompressible Maxwell fluid-saturated Darcy porous medium with the thickness d, which contains three stratifying agencies, i.e., the temperature T and two solute concentrations S_i (i = 1, 2). The lower and upper impermeable boundaries of the porous layer are kept at different constant temperatures, T_L and T_U (< T_L), and solute concentrations, S_iL and S_iU (< S_iL), respectively, with the sign convention that ΔT (= T_L − T_U) > 0 when the temperature is destabilizing and ΔS_i(= S_iL -S_iU) > 0 when a solute concentration is stabilizing. A Cartesian coordinate system is chosen such that the z-axis is vertically upwards and the x-axis is horizontal. The gravity is acting vertically downwards with the constant acceleration g = −g, where is the unit vector in the vertical direction. The equation of states is given by

(1)

where ρ is the fluid density, and ρ₀ is the density at the reference temperature and solute concentration. By analogy with the constitutive equation of the Maxwell fluid, a modified Darcy-Maxwell model is used to describe the flow in a porous medium. The governing nonlinear stability equations in dimensionless form are

(2)

(3)

(4)

(5)

where q = (u, v, w) is the velocity, p is the pressure, R_T = α_TgΔTKd/(νκ_T) is the thermal Darcy-Rayleigh number, R_Si = α_SigΔS_iKd/(νκ_T) (i = 1, 2) is the solute Darcy-Rayleigh number of the ith component, Λ₁ = λ₁κT/(d²ε) is the stress relaxation time parameter, Pr_D = νd²ε²/(κ_TK) is the Darcy-Prandtl number, τ_i = κS_i/κ_T (i = 1, 2) is the ratio of diffusivities, A = M/ε, in which ε is the porosity, K is the permeability of the porous medium, λ₁ is the relaxation time, M is the ratio of heat capacities, κ_T is the thermal diffusivity, κ_S1 and κ_S2 are the solute analogs of κ_T, α_T is the thermal expansion coefficient, and α_S1 and α_S2 are the solute analogs of α_T. It should be noted that the considered modified Darcy-Maxwell model includes the classical viscous Newtonian fluid, which is taken as a special case for Λ₁ = 0.

The following transformations are used in non-dimensionalizing the governing equations:

(6)

The boundary conditions are

(7)

(8)

(9)

The steady basic state is quiescent and considered as follows:

(10)

The basic state temperature, solute concentration, and pressure distributions are

(11)

where p₀ is the pressure at z=0. Since the temperature and solute concentrations vary linearly with respect to the vertical coordinate z, perturbing the basic state and following the standard stability analysis procedure, we can obtain the following nonlinear stability equations:

(12)

where is the linear differential operator given by

ζ = (ψ, T, S₁, S₂)^T, J(·, ·) stands for the Jacobian with respect to x and z, z, ψ(x, z, t) is the stream function, , and is is the Laplacian operator.

The boundary conditions now become

(13)

3 Weakly nonlinear stability analysis

A weakly nonlinear stability analysis is carried out to analyze the stability of the bifurcating periodic solution with the perturbation method. Accordingly, the dependent variables ψ, T, S₁, and S₂ and the thermal Darcy-Rayleigh number R_T are expanded in power series of a small perturbation parameter χ (≪ 1)^[37-38] as follows:

(14)

while the other parameters Pr_D, A, Λ₁, τ₁, τ₂, R_S1, and R_S2 are taken as given. At each stage in the expansion, a column vector may be defined by

The scaling for the time variable t is allowed such that

In Eq. (14), R_T1 is of no significance because it becomes zero due to the symmetry when the solvability condition is imposed. Substituting Eq. (14) into Eq. (12) and equating like powers of χ yield a series of linear partial differential equations at each order.

3.1 First-order system: linear instability analysis

At the leading order in χ, the equations are linear and homogeneous. Therefore, the first-order problem reduces to the linear instability problem for overstability, and we have

(15)

The solution of Eq. (15) satisfying the boundary conditions is assumed as follows:

(16)

(17)

(18)

(19)

where the overbar denotes the complex conjugate. The amplitudes (A₁ − D₁) and (A₁ − D₁) are allowed to vary over the slow time s. The relation between the amplitudes is obtained by substituting Eqs. (16)-(19) into Eq. (15) as follows:

(20)

where δ² = α² + π². The amplitude B₁ remains to be undetermined at this stage, and it will be determined from the solvability condition of the O(χ³) equation. In Eq. (15), the first equation gives an expression for the thermal Darcy-Rayleigh number for the occurrence of oscillatory convection, and this can be written (after clearing the complex quantities from the denominator) as follows:

(21)

and ω² satisfies

(22)

where

In the above equations,

Equation (22) shows that, for a suitable combination of the dimensionless parameters Pr_D, A, Λ₁, τ₁, τ₂, R_S1, and R_S2, it is possible to have either two or three different real positive values of ω² at the same wave number α. In that case, for each one of these frequency values (ω² > 0), there is a corresponding real value of the thermal Darcy-Rayleigh number on the oscillatory neutral curve.

For a singly diffusive case (R_S1 = R_S2 = 0), it is observed that oscillatory convection occurs, i.e.,

(23)

Thus, the necessary condition for the occurrence of oscillatory convection is

(24)

For a doubly diffusive case (e.g., R_S2=0), we obtain a dispersion relation quadratic in ω² as follows:

(25)

where

Since ω²>0 for the occurrence of oscillatory convection, a careful glance at Eq. (25) provides the necessary conditions as follows:

(26)

The condition ω = 0 corresponds to the stationary onset. Thus,

(27)

is the thermal Darcy-Rayleigh number, above which the layer is unstable. It is observed that Eq. (27) is independent of the viscoelastic parameter, which indicates that the effect of viscoelasticity appears only in the case of time dependent motion. We note that R_T^S attains its critical value at α_c = π and the critical thermal Darcy-Rayleigh number for the stationary onset is

(28)

We can obtain important information about the neutral stability curves in the (α, R_T)-plane by locating the bifurcation points, at which the steady and oscillatory neutral curves meet. These will occur on the steady neutral curve at the wave number α_b, for which ω = 0 is a root of Eq. (22). Thus, a₄(α_b) = 0 or equivalently

(29)

where . For any chosen parametric values, the critical value of R_T⁰ with respect to the wave number, which is denoted by R^S_Tc, is determined as follows. Equation (22) is solved first to determine the positive values of ω². If there are no positive values of ω², no oscillatory convection is possible. If there is only one positive value of ω², the critical value of R_T⁰ with respect to the wave number is calculated numerically from Eq. (21). If there are two or more positive values of ω², the least of R_T⁰ amongst these two ω² is retained to find the critical value of R_T⁰ with respect to the wave number.

3.2 Second-order system

At order χ², the equation is inhomogeneous and given by

(30)

where

(31)

(32)

(33)

The above relations suggest that the stream function, temperature, and solute concentration fields involve the terms of the frequency 2ω, and are independent of the fast time scale t. Thus, the stream function, temperature, and solute concentration fields at second-order can be articulated as follows:

(34)

(35)

(36)

(37)

where (ψ₂₀, ψ₂₂), (T₂₀, T₂₂), and (S₁₂₀, S₁₂₂, S₂₂₀, S₂₂₂) are, respectively, the stream function, temperature, and solute concentration fields, and they are independent of the fast time scale. The coefficients ψ₂₀, ψ₂₂, T₂₀, T₂₂, S₁₂₀, S₁₂₂, S₂₂₀, and S₂₂₂ are related to the amplitude at order χ as follows:

(38)

Therefore, the solution at this order is

(39)

3.3 Third-order system

Now, the equation at order χ³ is

(40)

where

(41)

(42)

(43)

(44)

The third-order equations have the solution as follows:

(45)

(46)

(47)

(48)

The right-hand-side terms of Eq. (40) have been evaluated from the previously known solutions at orders χ and χ². The algebra involved in finding the solutions at this order is straightforward. Therefore, only the results are presented here. We derive the solvability condition for Eq. (40), which is in the form of a first-order nonlinear ordinary differential equation (cubic Landau equation) for the unknown complex amplitude B₁ as follows:

(49)

where

(50)

(51)

In the above equations,

Let B₁ in the phase-amplitude form be as follows:

(52)

Let

Then, substituting the above equations into Eq. (49) yields

(53)

(54)

where ph(·) represents the phase shift. The magnitude and direction of the periodic convective solution and also the frequency shift are determined in Eq. (49). The nature of bifurcation depends on the sign of the quantity as follows:

(55)

If Ω > 0, the bifurcation is supercritical and stable. If Ω < 0, the bifurcation is subcritical and unstable. The temporal evolution of |B₁| can be expressed as a function of the initial amplitude B₀^[39] as follows:

(56)

From the above equation, it follows that |B₁| ~ B₀ exp(R_T2β_rs) when s → −∞ and |B₁| → 0, which is just as the linear theory. However, when s → ∞, which is independent of the value of B₀.

4 Heat and mass transfer

The time and area-averaged thermal Nusselt number and the solute Nusselt numbers (, ) are determined^[23] as follows:

(57)

(58)

(59)

5 Results and discussion

The triple diffusive convection in a layer of Maxwell fluid-saturated porous media has been investigated by performing both linear instability and weakly nonlinear stability analyses. The linear instability analysis has revealed some interesting results under certain conditions, which has not been observed hitherto in the literature. This has been achieved through a systematic study on the topology of neutral stability curves. Based on the weakly nonlinear stability analysis, a cubic Landau equation is derived, and the stability of the oscillatory bifurcating solution is analyzed. The heat and mass transfer is quantified in terms of the Nusselt numbers.

The neutral stability curves in the α, R_T plane are illustrated in Figs. 2 and 3 for different values of the stress relaxation parameter Λ₁, Darcy-Prandtl number Pr_D, and solute Darcy-Rayleigh number R_S1 for the chosen parametric values. These figures show that the stationary and oscillatory neutral curves are connected in a topological sense, which implies the requirement of single critical thermal Darcy-Rayleigh number to recognize the instability characteristics of the system. Moreover, the increases in Λ₁ (see Figs. 2(a) and 3(a)) and Pr_D (Fig. 2(b)) will decrease in the stability region. The stability region gets enlarged with an increase in the positive (i.e., stabilizing) values of R_S1, while an opposite trend can be seen with an increase in the negative values of R_S1 (see Fig. 3(b)).

The evolution of the neutral stability curves is displayed in Figs. 4(a)-4(f) for Pr_D=10, A=2, Λ₁=0.1, τ₁=0.159, τ₂=0.270, and R_S1=-15.92 for positive values of R_S2 varying from 63 to 59.3. These figures exhibit an altogether different behavior, which indicates significant ramifications on the linear instability of the system. It is seen that the oscillatory neutral curve is connected to the stationary neutral curve at two bifurcation points when R_S2=63 (see Fig. 4(a)). The oscillatory neutral curve is skewed slightly towards the lower wave number region, and also pinched on both sides, which gives rise to two values of R_T⁰ for some wave numbers. The upper maximum values of the oscillatory neutral curve lie on top of the minimum value of the stationary neutral curve. The bifurcation points move closer together when R_S2 is decreased to 61 (see Fig. 4(b)) and starts detaching from the stationary neutral curve when R_S2 is decreased further to 60 (see Fig. 4(c)). Till this stage, the linear instability of the system can be easily determined by a single value of R_T. With a further decrease in the value of R_S2, the closed convex oscillatory neutral curve totally detaches from the stationary curve (see Fig. 4(d)), and moves well below the minimum of the stationary neutral curve (see Fig. 4(e)) when R_S2=59.5. The importance of this type of detached oscillatory neutral curves is that it is necessary to have three critical values of R_T to specify the linear instability criteria of the Maxwell fluid-saturated porous layer instead of the usual single value. From Fig. 4(e), it is seen that the linear instability criteria involve three values of R_T, and may be viewed as follows. When R_T < R_T1 and R_T2 < R_T < R_T3, the layer is linearly stable. When R_T1 < R_T < R_T2 and R_T>R_T3, the layer is unstable. A further decrease in the value of R_S2 shows that the closed oscillatory neutral curve collapses to a point and ultimately disappears, which leaves only the stationary neutral curve when R_S2=59.3. Nonetheless, one salient feature that does not carry over from the case of Newtonian fluids is the absence of the heart-shaped oscillatory neutral curve within the maxima at the same thermal Darcy-Rayleigh number and different wave numbers.

Figures 5 and 6 exemplify the stability boundaries for the same parametric values considered in Fig. 4. From Fig. 5, we can see that the stability boundary can be viewed individually in three regions. To the right of the point of (A) (R_S2>60.133 7), oscillatory instability first occurs at a lower value of R_T, then stationary instability occurs, and there is a critical thermal Darcy-Rayleigh number R_Tc. To the left of the cusp (B) (R_S2 < 59.390 8), instability occurs as stationary convection, oscillatory instability does not occur, and again there is one value of R_Tc. The enlarged version of the figure clearly reveals three distinct regions of R_S2, and it is seen that, in some range of R_S2, three values of R_Tc are needed to specify the linear instability criteria.

The sensitivity of the viscoelastic parameter Λ₁ on disconnected oscillatory neutral curves is elucidated in Figs. 7(a) and 7(b). Although Fig. 7(a) makes obviously the requirement of three values of R_Tc to specify the linear instability of the system, a slight deviation in the value of Λ₁ (see Fig. 7(b)) entirely changes the instability mode. That is, the instability ceases to be oscillatory and the convective onset turns out to be stationary, which represents the sufficiency of single value of R_Tc to specify the linear instability criteria. Thus, the viscoelastic property of the fluid significantly affects the instability characteristics of the system.

It is possible to estimate the parameter set for the boundary separating stationary and oscillatory solutions. This behavior is illustrated in Fig. 8 on a (Pr_D, Λ₁) plane for different values of R_S2 obtained from Eqs. (21) and (28). It is seen that, for fixed values of Pr_D and R_S2, there exists a value of Λ₁=Λ₁^*, where R_Tc⁰=R_Tc^S. The value of Λ₁^* decreases with the increases in Pr_D and R_S2. In other words, the critical thermal Darcy-Rayleigh number for the onset of both stationary and oscillatory convection coincides at well-defined parametric values. Therefore, a codimension-two bifurcation occurs. The figure also indicates that the value of the stress relaxation parameter, which represents the crossover between stationary and oscillatory bifurcations, decreases with the increase in the Darcy-Prandtl number. Moreover, the oscillatory region increases with the increase in R_S2, which indicates that the presence of the stronger concentration gradient is to support the oscillatory convection as the preferred instability mode.

By use of a weakly nonlinear stability analysis, a cubic Landau equation is derived, and the stability of the oscillatory bifurcating solution is analyzed, which powerfully depends on the viscoelastic parameter. It is observed that subcritical/supercritical bifurcation is possible, which depends on the choice of the physical parameters. Nonetheless, the bifurcating solution is found to be supercritical (stable) in the absence of additional solute concentration fields. The vigor of triple diffusive convection in the Maxwell fluid-saturated porous medium can be understood by heat and mass transfer. This has been achieved by estimating the thermal and solute Nusselt numbers for oscillatory convection. The time and area-averaged thermal Nusselt number () and solute Nusselt numbers ( and ) are displayed as a function of R_T for different values of Λ₁ and Pr_D in Figs. 9(a) and 9(b), respectively. When R_T increases, the heat and mass transfer increases. When Λ₁ and Pr_D increase, the heat and mass transfer characteristics of oscillatory convection increase. This may be credited to the fact that the increase in Λ₁ is to increase the overstability vibration. Therefore, its effect is to increase the heat and mass transfer.

6 Conclusions

The triple diffusive convection in a Maxwell fluid-saturated porous medium is investigated by use of both linear and weakly nonlinear stability analyses. Moreover, disconnected convex oscillatory neutral curves are found to occur for some choices of the parametric values, which indicates the necessity of three critical thermal Darcy-Rayleigh numbers to specify the linear instability criteria instead of the usual single value. However, it is observed that the instability onset of the motionless basic state cannot occur via the simultaneous passage of two sets of complex conjugate temporal eigenvalues from the left half-plane to the right half-plane at the same critical thermal Darcy-Rayleigh number as observed in the case of Newtonian fluids. The onset of oscillatory convection increases with the increases in the relaxation parameter and the Darcy-Prandtl number. The threshold value of the relaxation parameter, at which codimension-two bifurcation occurs, decreases with the increase in the Darcy-Prandtl number. The stability of the bifurcating oscillatory solution is analyzed by deriving a cubic Landau equation. There is a possibility of bifurcation to be either subcritical (unstable) or supercritical (stable), which depends on the choice of the physical parameters. The variations of the time and area-averaged thermal and solute Nusselt numbers with respect to the thermal Darcy-Rayleigh number are presented. The heat and mass transfer rate increases with the increases in the relaxation parameter and the Darcy-Prandtl number.

Acknowledgements One of the authors (K. R. RAGHUNATHA) wishes to thank the 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]). The authors wish to thank the reviewers for their helpful suggestions.