Cross Diffusion Effect on Linear and Nonlinear Double Diffusive Convection in a Viscoelastic Fluid Saturated Porous Layer with Internal Heat Source

Abstract

Double diffusive convection in a binary viscoelastic fluid saturated porous layer in the presence of a cross diffusion effect and an internal heat source is studied analytically using linear and nonlinear stability analysis. The linear stability theory is based on the normal mode technique, while the nonlinear theory is based on a minimal representation of truncated double Fourier series. The modified Darcy law for the viscoelastic fluid of the Oldroyd type is considered to model the momentum equation. The onset criterion for stationary and oscillatory convection and steady heat and mass transfer have been obtained analytically using linear and nonlinear theory, respectively. The combined effect of an internal heat source and cross diffusion is investigated. The effects of Dufour, Soret, internal heat, relaxation and retardation time, Lewis number and concentration Rayleigh number on stationary, oscillatory, and heat and mass transport are depicted graphically. Heat and mass transfer are presented graphically in terms of Nusselt and Sherwood numbers, respectively. It is reported that the stationary and oscillatory convection are significantly influenced with variation of Soret and Defour parameters. An increment of the internal heat parameter has a destabilizing effect as well as enhancing the heat transfer process. On the other hand, an increment of internal heat parameter has a variable effect on mass transfer. It is found that there is a critical value for the thermal Rayleigh number, below which increasing internal heat decreases the Sherwood number, while above it increasing the internal heat increases the Sherwood number.

1. Introduction

Mathematical modeling of double-diffusive free convection in a fluid-drenched porous module has increased in popularity over the last few decades due to its numerous applications in diverse contexts, including nuclear waste repository, impurities within water resources, atmospheric effluents, material processing, oil reservoir modeling, geothermal fields, and biochemical processes. The thermal instability of a Newtonian fluid-drenched porous module was initially investigated by [1] and by Lapwood [2] using linear stability analysis. Their results were further authenticated experimentally by Elder [3] and Katto [4]. On the other hand, Ingham and Pop [5], Vafai [6] and Nield and Bejan [7] conducted comprehensive reviews of pertinent literature and the fundamental applications of thermal performance in a porous medium. External regulation of convection streams is important for enhancing or diminishing the flow performance in physical systems. At present, several techniques are used for controlling the convective instability of such systems, including those based on gravity, rotation, thermal modulation and imposed magnetic field. In their studies, Gershuni et al. [8] and Gresho and Sani [9] investigated the effect of gravity modulation in a differentially heated fluid module. Their findings indicate that the undesirable buoyancy-driven flows in a differentially heated fluid layer can be effectively stabilized by appropriate vertical vibrations. Similarly, Malashetty and Padmavathi [10] investigated the stability solution of a fluid-drenched porous module depending on the amplitude of vertical vibrations. These authors found that low-frequency vertical vibrations significantly affected the system stability limits. The magnetic field effect on double-diffusive convection in a micropolar nanofluid was established by Abidi et al. [11], who found that the heat and mass transfer are influenced by the increasing Hartman number. More recently, Kumar et al. [12] examined a weak nonlinear heat transfer in a viscoelastic fluid-saturated porous medium under the effect of rotational modulation and found that a modified Prandtl number suspends the heat transfer process. Flow and heat transfer in non-Newtonian nanofluids over porous surfaces was studied by Maleki et al. [13], who found these nanofluids to be superior to their Newtonian counterparts for injection and impermeable plates. The authors also noted that, in injection, the heat transfer can be decreased by utilizing non-Newtonian nanofluids. In an earlier study, Govender [14] discussed the effect of g-jitter on convection in a differentially heated homogeneous porous module and predicted the exact changeover point from the synchronous to the subharmonic solution by plotting Mathieu’s stability charts. These results were subsequently confirmed by Saravanan and Purusothaman [15], who also studied the effects of anisotropy parameters of a porous module subjected to vertical vibrations using the Brinkman model. Thermal convection instability in a mushy layer subject to vertical vibration has been examined by Srivastava and Bhadauria [16] using linear stability analysis. These authors discussed the influence of pertinent parameters on mushy layer stability. A large number of fundamental studies on the stability of a rotating fluid-saturated porous module has also been conducted due to its diverse industrial requirements and applications. Govender [17] and Malashetty and Swamy [18] provided a detailed overview of rotating fluid-saturated porous module characteristics and applications. On the other hand, Maleki et al. [19] focused on the influences of radiation and slip boundary conditions on heat transfer over a porous plate. They also examined the effect of nanoparticles on nanofluid flow and heat transfer and noted significant changes in heat generation, viscous dissipation and radiation, which led to a decrease in heat transfer.

Some fluids found in foodstuffs, medicines, polymers and paints are of intricate rheological nature, and can thus be considered viscoelastic fluids for the purpose of modeling, in which case a viscoelastic model rather than the Newtonian model should be adopted. Thermal instability in a viscoelastic fluid was initially considered by Herbert [20] and Green [21]. Several studies in this domain focused on buoyancy-driven thermal instability in a viscoelastic fluid-drenched porous module, including those published by Kim et al. [22], Malashetty et al. [23], Naryana et al. [24], Bhadauria and Kiran [25], and Altawallbeh et al. [26]. The stability of double-diffusive convection of a Maxwell fluid in a porous module has also been widely investigated over the last few decades. For instance, Awad et al. [27], Wang and Tan [28], and Zhao et al. [29] found that the Maxwell fluid properties exert a significant influence on system stability. There are numerous practical situations in which the convective stability of the porous module can be compromised when subjected to an internal heat source, as is the case of electrometallurgical development for treating spent nuclear fuel, radioactive decay of unstable isotopes and exothermic reaction in porous materials. The onset of convection in a porous module with an interior heat source was examined by Parthiban and Parthiban [30], Bhadauria [31], Altawallbeh et al. [32], and Badruddin [33], among other authors. Recently, Altawallbeh et al. [34] examined the influence of a magnetic field on the instability of a fluid layer saturated with a viscoelastic fluid with an internal heat source. These authors found that the magnetic field has a stabilizing effect in such cases. Similarly, Rana and Khurana [35] studied the thermoconvective instability of a water-based nanofluid layer in the presence of a magnetic field when the fluid and nanoparticles are in a local thermal non-equilibrium. They found that the isotherms for both fluid and particle phases are flattened near the wall, suggesting an increase in the convection within the nanofluid layer. The effects of internal heat on the instability of an inclined anisotropic porous layer were investigated by Storesletten and Rees [36]. They found that the onset criteria in the porous layer were influenced by both the anisotropy parameter and the inclination angle. Maleki et al. [37] studied the influence of heat generation and viscous dissipation on heat transfer and the fluid flow of pseudo-plastic nanofluid over a moving permeable plate. These authors observed that, for both Newtonian and non-Newtonian nanofluids, the local Nusselt number decreased as the injection parameter, heat generation or volume fraction of nanoparticles increased.

The concurrent process of heat and mass transfer in a fluid can be studied by applying the theory of temperature gradient, concentration gradient and material composition, which causes density variances and thus implicitly induces the fluid to the flow. Heat flux generated due to the variance in concentration is denoted as thermal diffusion or the Soret effect, and mass flux resulting from the variance in temperature is referred to as the Dufour effect, both of which are highly relevant to chemical engineering, petrology, hydrology, disposal of nuclear waste, geosciences, and so on. The cross-diffusion effects in a porous module have been extensively studied. For example, Taslim and Narusawa [38] established a connection between the cross-diffusion problem and double-diffusion problem in a porous module. On the other hand, Rudaraiah and Siddheshwar [39] evaluated nonlinear double-diffusive convection in a porous module subject to Soret and Dufour effects. Many researchers have also investigated the effect of cross-diffusion coefficients on the double-diffusion convection in a porous module. Specifically, Gaikwad et al. [40] focused on the influence of cross-diffusion on linear and nonlinear double-diffusive convection in a fluid-saturated anisotropic porous layer. The Soret effect on double-diffusive convection in a couple-stress fluid-saturated porous medium was reported by Malashetty et al. [41], who found that it may lead to system instability. Recently, Altawallbeh et al. [42] investigated the effect of the Soret parameter on the linear and nonlinear stability of a porous layer when the fluid and solid phase are in a thermal non-equilibrium state.

The literature review presented above reveals a paucity of studies on the double-diffusion driven convective instability in a non-Newtonian fluid-drenched porous module with an interior heat generation subjected to cross-diffusion effects. This gap in extant literature has motivated the present study, the aim of which is to investigate the combined effect of cross-diffusion and an interior heat source on a viscoelastic fluid-drenched isotropic porous module. In this work, the onset criteria for the double-diffusion driven convective instability are assessed by performing both linear and nonlinear analyses. The results yielded by this study will contribute to the better understanding of the onset criterion for the double-diffusion convection in a viscoelastic fluid-drenched isotropic porous module with an interior heat generation subjected to the cross-diffusion effects, which is a frequently encountered phenomenon in different systems and industries.

2. Mathematical Formulation

The physical model under consideration is a viscoelastic fluid saturated horizontal isotropic porous layer with an internal heat source, confined between two parallel horizontal planes at $z=0$ and $z=d$ with a distance d apart. A Cartesian frame of reference is chosen in such a way that the origin lies on the lower plane and the z-axis as vertically upward, where the gravity force $\mathbf{g}$ is acting vertically downward. Adverse temperature and concentration gradients are applied across the porous layer, and the lower and upper planes are kept at temperature $T0+\Delta T$, concentration $C_0+\Delta C$, and $T_0,C_0$, respectively, where $\Delta T$ and $\Delta C$ are temperature difference and concentration difference between walls, respectively, as shown in Figure 1. In the momentum equation, the modified Darcy law has been applied for an Oldroyd type viscoelastic fluid. Both Soret and Dufour terms take to account in the temperature and concentration equations. In addition, the internal heat source parameter is involved in the temperature equation. Finally, the Oberbeck–Boussinesq approximation has been applied to account for the effect of density variations. Under these assumptions, the governing equations are given as [23]:

$$\begin{array}{ccc}& & \nabla .\mathbf{q}=0\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & \left({\overline{\lambda }}_1\frac{\partial }{\partial t}+1\right)\left(\frac{{\rho }_0}{\delta }\frac{\partial \mathbf{q}}{\partial t}+\nabla P-\rho \mathbf{g}\right)+\frac{\mu }{K}\left({\overline{\lambda }}_2\frac{\partial }{\partial t}+1\right)\mathbf{q}=0\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & \gamma \frac{\partial T}{\partial t}+\mathbf{q}.\nabla T={\kappa }_T{\nabla }^{2}T+K_1{\nabla }^{2}C+Q\left(T-T_0\right)\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & ϵ\frac{\partial C}{\partial t}+\mathbf{q}.\nabla C={\kappa }_s{\nabla }^{2}C+K_2{\nabla }^{2}T\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & \rho ={\rho }_0[1-{\alpha }_T(T-T_0)+{\alpha }_s(C-C_0)]\hfill \end{array}$$

(5)

subject to the following boundary conditions:

$$\begin{array}{ccc}\hfill T& =& T_0+▵T,C=C_0+▵C\mathrm{at}z=0\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill T& =& T_0,C=C_0\mathrm{at}z=d\hfill \end{array}$$

(7)

where ${\overline{\lambda }}_1$ is the relaxation time parameter, ${\overline{\lambda }}_2$ is the retardation time parameter, $\mathbf{q}$ is velocity $(u,v,w)$, $\mu$ is the dynamic viscosity, Q is the internal heat source, K is the permeability, ${\kappa }_T$ is the thermal diffusivity, T is temperature, ${\alpha }_T$ is the thermal expansion coefficient, ${\alpha }_s$ is the concentration expansion coefficient, ${\kappa }_s$ is the concentration diffusivity of the fluid, $\gamma$ is the ratio of heat capacities, $K_1$ and $K_2$ are the cross diffusion due to C and T components, respectively, $\rho$ is the density, $\mathbf{g}=(0,0,-g)$ is gravitational acceleration, and ${\rho }_0$ is the reference density, for the viscoelastic fluid of Oldroyd type ${\lambda }_1\ge {\lambda }_2$. When ${\lambda }_2=0$, the Maxwell viscoelastic model is recovered. When ${\lambda }_1={\lambda }_2=0$, it reduces to a Newtonian fluid.

2.1. Basic State

Now, for the basic state, the fluid layer is stable and assumed to be quiescent. Therefore it has no motion and the velocity vector is set equal to zero, where the heat and mass are transferred by conduction only. Thus, the basic state is given by:

$${\mathbf{q}}_b=(0,0,0),P=P_b\left(z\right),T=T_b\left(z\right),C=C_b\left(z\right),\rho ={\rho }_b\left(z\right)$$

(8)

Substituting Equation (8) into Equations (1)–(4), we obtain:

$$\begin{array}{ccc}& & \frac{d{P}_b}{dz}=-\rho g\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & {\kappa }_T\frac{d^2{T}_b}{d{z}^2}+K_1\frac{d^2{C}_b}{d{z}^2}+Q(T_b-T_0)=0\hfill \end{array}$$

(10)

$$\begin{array}{ccc}& & {\kappa }_s\frac{d^2{C}_b}{d{z}^2}+K_2\frac{d^2{T}_b}{d{z}^2}=0\hfill \end{array}$$

(11)

where b refers to the basic state. The basic state temperature and concentration satisfy Equations (10) and (11) with boundary conditions $T=T_0+\Delta T,C=C_0+\Delta C$ at $z=0$, and $T=T_0,C=C_0$ at $z=d$.

By solving the resulted equations, the conduction state solution is obtained, and is given by:

$$\begin{array}{ccc}\hfill T_b\left(z\right)& =& \frac{\Delta Tsin\left(A_1(d-z)\right)}{sin\left(A_{1}d\right)}+T_0\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill C_b\left(z\right)& =& \frac{K_2\Delta T}{{\kappa }_s}\frac{sin\left(A_1(d-z)\right)}{sin\left(A_{1}Q\right)}+\left(1-\frac{z}{d}\right)\left(\Delta C+\frac{K_2\Delta T}{{\kappa }_s}\right)+C_0\hfill \end{array}$$

(13)

where $A_1=\sqrt{\frac{Q{\kappa }_s}{{\kappa }_T{\kappa }_s-K_1{K}_2}}$. Equations (12) and (13) represent the heat and mass transfer in conduction state only, where the layer is stable. This means no motion or convection exist in the layer and the transferred heat and mass occur by conduction only.

2.2. Perturbed State

Now, to study the stability of the system, we superimpose the infinitesimal perturbations on the basic state. Thus, the basic state is perturbed with thermal perturbation as follows:

$$\mathbf{q}={\mathbf{q}}_b+\acute{\mathbf{q}},T=T_b+\acute{T},C=C_b+\acute{C},P=P_b+\acute{P},\rho ={\rho }_b+\acute{\rho }$$

(14)

where the primes denote that the quantities are infinitesimal perturbations. Now, we substitute Equation (14) into Equations (1)–(5), and using basic state Equations (9)–(11) we obtain:

$$\begin{array}{ccc}& & \nabla .{\mathbf{q}}^{\prime }=0\hfill \\ & & \left(1+{\overline{\lambda }}_1\frac{\partial }{\partial t}\right)\left(\nabla P^{\prime }+\frac{{\rho }_0}{\delta }\frac{\partial {\mathbf{q}}^{\prime }}{\partial t}-{\rho }_0{\beta }_{t}g{T}^{\prime }+{\rho }_0{\beta }_{s}g{C}^{\prime }\right)\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & +\frac{\mu }{K}\left(1+{\overline{\lambda }}_2\frac{\partial }{\partial t}\right){\mathbf{q}}^{\prime }=0\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & \gamma \frac{\partial T^{\prime }}{\partial t}+({\mathbf{q}}^{\prime }.\nabla )T^{\prime }+w^{\prime }\frac{d{T}_b}{dz}={\kappa }_T{\nabla }^2{T}^{\prime }+K_1{\nabla }^2{C}^{\prime }+Q{T}^{\prime }\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& & ϵ\frac{\partial C^{\prime }}{\partial t}+({\mathbf{q}}^{\prime }.\nabla )C^{\prime }+w^{\prime }\frac{d{C}_b}{dz}={\kappa }_s{\nabla }^2{C}^{\prime }+K_2{\nabla }^2{T}^{\prime }\hfill \end{array}$$

(18)

Eliminate the pressure from Equation (16) by operating the curl twice to obtain:

$$\begin{array}{ccc}\hfill & & \left(1+{\overline{\lambda }}_1\frac{\partial }{\partial t}\right)\left[\frac{{\rho }_0}{\delta }\frac{\partial }{\partial t}\left(\frac{\partial u^{\prime }}{\partial z}-\frac{\partial w^{\prime }}{\partial x}\right)+{\rho }_{0}g{\beta }_T\frac{\partial T^{\prime }}{\partial x}-{\rho }_0{\beta }_{s}g\frac{\partial C^{\prime }}{\partial x}\right]\hfill \\ & & +\frac{\mu }{K}\left(1+{\overline{\lambda }}_2\frac{\partial }{\partial t}\right)\left[\frac{\partial u^{\prime }}{\partial z}-\frac{\partial w^{\prime }}{\partial x}\right]=0\hfill \end{array}$$

(19)

Now, introduce the stream function $u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x}$, and use the following transformation:

$$\begin{array}{ccc}\hfill & & (x^{*},z^{*})=(\frac{x}{d},\frac{z}{d}),t^{*}=\frac{t{\kappa }_{T_z}}{d^2},{\psi }^{*}=\frac{\psi }{{\kappa }_{T_z}},T^{*}=\frac{T^{\prime }}{▵T},S^{*}=\frac{C^{\prime }}{▵C},\hfill \\ & & {\overline{\lambda }}_1=\frac{{\lambda }_1{d}^2}{{\kappa }_T},{\overline{\lambda }}_2=\frac{{\lambda }_2{d}^2}{{\kappa }_T}\hfill \end{array}$$

(20)

into Equations (15), (17)–(19); we obtain:

$$\begin{array}{ccc}& & \left[\left(1+{\lambda }_1\frac{\partial }{\partial t^{*}}\right)\frac{1}{Va}\frac{\partial }{\partial t^{*}}{\nabla }^2+\left(1+{\lambda }_2\frac{\partial }{\partial t^{*}}\right){\nabla }^2\right]{\psi }^{*}+\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \times \left[R{a}_T\frac{\partial T^{*}}{\partial x^{*}}-R{a}_s\frac{\partial S^{*}}{\partial x^{*}}\right]=0\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & \left[\gamma \frac{\partial }{\partial t^{*}}-{\nabla }^2-R_i\right]T^{*}=\frac{\partial {\psi }^{*}}{\partial x^{*}}f\left(z^{*}\right)+\frac{\partial ({\psi }^{*},T^{*})}{\partial (x^{*},z^{*})}+Du\frac{R{a}_s}{R{a}_T}{\nabla }^2{S}^{*}\hfill \end{array}$$

(22)

$$\begin{array}{ccc}& & \left[ϵ\frac{\partial }{\partial t^{*}}-\frac{1}{Le}{\nabla }^2\right]S^{*}=\frac{\partial {\psi }^{*}}{\partial x^{*}}b\left(z^{*}\right)+\frac{\partial ({\psi }^{*},S^{*})}{\partial (x^{*},z^{*})}+S_r\frac{R{a}_T}{R{a}_s}{\nabla }^2{T}^{*}\hfill \end{array}$$

(23)

where $Va=\frac{Pr}{Da}$ is the Vadasz number, $Da=\frac{K}{d^2}$ is the Darcy number, $Pr=\frac{v}{k_T}$ is the Prandtl number, $R{a}_T=\frac{{\beta }_{T}gK▵Td}{\nu {\kappa }_T}$ the thermal Rayleigh number, $R{a}_s=\frac{{\beta }_{s}gK▵Cd}{\nu {\kappa }_T}$ is the concentration Rayleigh number, $R_i=\frac{Q{d}^2}{{\kappa }_T}$ is the internal Rayleigh number, $\nu =\mu /{\rho }_0$ is the kinematic viscosity, $Le=\frac{{\kappa }_T}{{\kappa }_s}$ is the Lewis number, $S_r=\frac{K_2{\alpha }_s}{{\kappa }_T{\alpha }_T}$ is the Soret parameter, $Du=\frac{K_1{\alpha }_T}{{\kappa }_T{\alpha }_s}$ is the Defour parameter, ${\lambda }_1$ is the relaxation parameter, and ${\lambda }_2$ is the retardation parameter; ${\displaystyle f\left(z^{*}\right)=\frac{d{T}_b^{*}}{d{z}^{*}}}$ and ${\displaystyle b\left(z^{*}\right)=\frac{d{C}_b^{*}}{d{z}^{*}}}$. $T_b^{*}$ and $S_b^{*}$ are given in dimensionless form as:

$$\begin{array}{ccc}\hfill T_b^{*}\left(z\right)& =& \frac{sin\left(A_2(1-z)\right)}{sin\left(A_2\right)}\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill S_b^{*}\left(z\right)& =& \frac{-LeSrR{a}_{T}sin\left(A_2(1-z)\right)}{sin\left(A_2\right)}+\left(1-z\right)\left(\frac{LeSrR{a}_T}{R{a}_s}+1\right)\hfill \end{array}$$

(25)

where ${\displaystyle A_2=\sqrt{\frac{R_i}{1-DuSrLe}}}$.

The viscoelastic character depends on the relaxation time parameter ${\lambda }_1$ and the retardation time parameter ${\lambda }_2$, where these two parameters give the viscoelastic fluid its viscosity and elasticity. The smaller the relaxation parameter, the more fluid the material appears. ${\lambda }_1$ is of order one for most viscoelastic fluids. For some special cases, ${\lambda }_1$ is most likely in the range $[0.1,2]$ and ${\lambda }_2$ is in the range $[0.1,1]$ for a diluted polymeric solution. At present, there are no experimental data for the ranges of the parameters, and hence a wide range of parameters were considered. Some parameters, such as $Le$ and $Da$, depend on the fluid properties [43].

Now, Equations (21)–(23) are solved for stress free, isothermal, and isohaline scenarios. Therefore the boundary conditions are expressed as:

$${\psi }^{*}=\frac{{\partial }^2{\psi }^{*}}{\partial z^{*2}}=T^{*}=S^{*}=0\mathrm{at}{z}^{*}=0\mathrm{and}{z}^{*}=1$$

(26)

For the next steps the asterisks are removed, and $\gamma$ and $ϵ$ are set equal to unity for simplicity.

3. Linear Stability Theory

Stability means stability with respect to all infinitesimal disturbances. Thus, the reaction of the system with all disturbances should be examined. This can be established by considering an arbitrary disturbance as a superposition of certain basic possible modes and testing the stability of the system with all of these modes [44]. In the present problem, an arbitrary disturbance in terms of periodic waves will be analyzed. Therefore we assume the solution to be periodic waves as given in Equation (27).

In this section, the linear stability analysis is discussed. The thresholds of marginal and oscillatory convection are predicted. Equations (21)–(23) form an eigenvalue problem subjected to the boundary conditions given in Equation (26). For the onset of convection, the nonlinear terms in Equations (21)–(23) must be ignored and the resulting linear system is sufficient to study the onset of convection. Thus, to make the linear stability study, we neglect the nonlinear terms in Equations (21)–(23), and assume the solutions to be periodic waves of the form:

$$\begin{array}{ccc}\hfill \left[\begin{array}{c}\psi \\ T\\ S\end{array}\right]& =& e^{\sigma t}\left[\begin{array}{c}\mathsf{\Psi }sin\alpha \pi x\\ \mathsf{\Theta }cos\alpha \pi x\\ \mathsf{\Phi }cos\alpha \pi x\end{array}\right]sin\pi z\hfill \end{array}$$

(27)

where ($\mathsf{\Psi },\mathsf{\Theta },\mathsf{\Phi }$) are the amplitudes of ($\psi ,T,S$), respectively. $\sigma ={\sigma }_r+i{\sigma }_i$ is a growth rate and in general a complex quantity and $\alpha$ is a horizontal wave number. $\alpha$ is the wave number that corresponds to disturbances. The reaction of the system to a disturbance can be studied in terms of system reaction to all disturbances of all associated wave numbers $\alpha$. This means that the stability will depend on its stability as a result of disturbances of all wave numbers $\alpha$ and instability will follow from its instability as a result of disturbances of even one wave number. The infinitesimal perturbations may grow or damp depending on the values of the growth rate $\sigma$. Now, substitute Equation (27) into the linearized form of Equations (21)–(23) to obtain:

$$\begin{array}{cc}& \left[\left(1+{\lambda }_1\sigma \right)\frac{\sigma }{Va}{a}^2+{\delta }^2\left(1+{\lambda }_2\sigma \right)\right]\mathsf{\Psi }\end{array}$$

$$\begin{array}{cc}& +\left(1+{\lambda }_1\sigma \right)\left[\alpha \pi \left(R{a}_T\mathsf{\Theta }-R{a}_s\mathsf{\Phi }\right)\right]=0\end{array}$$

(28)

$$\begin{array}{cc}& \left(\sigma +{\delta }^2-R_i\right)\mathsf{\Theta }-2\alpha \pi F\mathsf{\Psi }+{\delta }^{2}Du\frac{R{a}_s}{R{a}_T}\mathsf{\Phi }=0\end{array}$$

(29)

$$\begin{array}{cc}& \left(\sigma +\frac{{\delta }^2}{Le}\right)\mathsf{\Phi }-2\alpha \pi B\mathsf{\Psi }+a^2{S}_r\frac{R{a}_T}{R{a}_s}\mathsf{\Theta }=0\end{array}$$

(30)

where ${\delta }^2={\pi }^2(1+{\alpha }^2)$. Now, change Equations (28)–(30) to a matrix form:

$$\begin{array}{ccc}\hfill \left[\begin{array}{ccc}A& R{a}_T& -R{a}_s\\ -2\alpha \pi F& \sigma +{\delta }^2-R_i& {\delta }^{2}Du\frac{R{a}_s}{R{a}_T}\\ -2\alpha \pi B& {\delta }^2{S}_r\frac{R{a}_T}{R{a}_s}& \sigma +\frac{{\delta }^2}{Le}\end{array}\right]\left[\begin{array}{c}\mathsf{\Psi }\\ \mathsf{\Theta }\\ \mathsf{\Phi }\end{array}\right]& =& \left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]\hfill \end{array}$$

(31)

where:

$$\begin{array}{ccc}\hfill & & A=\frac{{\delta }^2}{\alpha \pi }\left(\frac{\sigma }{Va}+\frac{1+{\lambda }_2\sigma }{1+{\lambda }_1\sigma }\right),F={\int }_0^{1}f\left(z\right){sin}^2\left(\pi z\right)dz,\hfill \\ & & B={\int }_0^{1}b\left(z\right){sin}^2\left(\pi z\right)dz.\hfill \end{array}$$

For the non-trivial solution of Equation (31), the determinant must equal to zero, which implies an expression of the thermal Rayleigh number, and is given by:

$$\begin{array}{ccc}\hfill R{a}_T& =& \frac{1}{\mathsf{\Omega }}\left[-A(DuLeSr{\delta }^4-Le{\delta }^2\sigma -{\delta }^4+Le{R}_i\sigma -Le{\sigma }^2+R_i{\delta }^2-{\delta }^2\sigma )\right.\hfill \\ & +& \left.\alpha \pi BR{a}_{s}Le(Du{\delta }^2+{\delta }^2-R_i+\sigma )\right]\hfill \end{array}$$

(32)

where $\mathsf{\Omega }=\alpha \pi F(LeSr{\delta }^2+Le\sigma +{\delta }^2)$.

3.1. Marginal Stationary State

The growth rate sigma in general is a complex quantity and is given as $\sigma ={\sigma }_r+{\sigma }_i$. For (${\sigma }_r<0$), the system is stable, and for $({\sigma }_r>0)$, the system is unstable. For the existence of marginal curves (stationary or oscillatory) or neutral stability curves, the real part of sigma must be zero $({\sigma }_r=0)$. Now, for the existence of a marginal stationary state the imaginary part must also equal zero $({\sigma }_i=0)$. As a result, sigma becomes equal to zero [32]. Now, for the stationary Rayleigh number, substitute $\sigma =0$ into Equation (32). Then, the stationary Rayleigh number is given as:

$$\begin{array}{ccc}\hfill R{a}_T^{st}& =& ((1/4)R_i-\phi {\pi }^2)((-\phi ){\delta }^6+R_i{\delta }^4+Le{\pi }^{2}R{a}_s{\alpha }^2(R_i-(Du+1){\delta }^2))\hfill \\ & /& \left[{\alpha }^2{\pi }^2((-(1/4)R_{i}Sr(Du+1)L{e}^2+(LeSr+1)\phi {\pi }^2){\delta }^2\right.\hfill \\ & +& \left.(1/4)L{e}^2{R}_i^{2}Sr)\right]\hfill \end{array}$$

(33)

where $\phi =1-DuSrLe$.

It is noticeable that the expression of the stationary Rayleigh number that appears in Equation (33) is independent of the viscoelastic parameters, the relaxation time parameter ${\lambda }_1$ and the retardation time parameter ${\lambda }_2$, where these two parameters give the viscoelastic fluid its viscosity and elasticity. Thus, when these two parameters do not appear in the expression of the stationary Rayleigh number, it means that the onset of convection via the stationary state will not be affected by these parameters. Therefore, the present fluid layer behaves as a Newtonian binary fluid layer. This is because the fluid in the base state has no flow and the viscoelasticity will not affect the onset of stationary convection when the flow is steady and weak.

Now, in the case with no internal heat source, Equation (33) becomes:

$$\begin{array}{c}\hfill R{a}_T^{st}=\frac{Le(Du+1)}{(LeSr+1)}R{a}_s-\frac{{\pi }^2{({\alpha }^2+1)}^2(DuLeSr-1)}{\left({\alpha }^2(LeSr+1)\right)}\end{array}$$

(34)

which recovers the results obtained by Malashetty et al. [45].

When the cross diffusion effects and internal heat source are absent, the result by Malashetty et al. [23] is recovered, when an isotropic porous layer is considered. Thus, Equation (33) becomes:

$$\begin{array}{c}\hfill R{a}_T^{st}=\frac{{\pi }^2{({\alpha }^2+1)}^2}{{\alpha }^2}+LeR{a}_s\end{array}$$

(35)

For $R{a}_s=0$, we have $R{a}_T^{st}=\frac{{\pi }^2}{{\alpha }^2}{(1+{\alpha }^2)}^2$, and have the critical value $R{a}_T^{st}=4{\pi }^2$ for ${\alpha }_c^2=1$, which recovers the results obtained by Horton and Rogers [1] and Lapwood [2] for a single component convection in a porous medium.

The critical Rayleigh number for steady state convection in the presence of an internal heat source, the Soret and Dufour effects, is calculated from Equation (33) with variations of the parameters.

3.2. Marginal Oscillatory State

For the existence of a marginal oscillatory state, the real part must equal zero $({\sigma }_r=0)$, but the imaginary part does not equal zero $({\sigma }_i\ne 0)$. Thus, in Equation (32) we set $\sigma =i\sigma$ and cleared up the complex quantities from the denominator to obtain:

$$\begin{array}{ccc}\hfill R{a}_T^{Osc}& =& \omega /2Va({\lambda }_1^2{\sigma }^2+1){\alpha }^2{\pi }^2\left[B^2({\sigma }^2+{(Du{\delta }^2+{\delta }^2-R_i)}^2)R{a}_s^2\right.\hfill \\ & -& \left(2({\sigma }^2+(Sr{\delta }^2+\frac{1}{Le}{\delta }^2)(Du{\delta }^2+{\delta }^2-R_i))\right)FBR{a}_s\hfill \\ & +& \left.({\sigma }^2+{(Sr{\delta }^2+\frac{1}{Le}{\delta }^2)}^2)F^2\right]\hfill \end{array}$$

(36)

where $\omega$ is not given here for brevity.

From the previous equation, it follows that ${\sigma }_i=0$ gives a stationary convection and $({\sigma }_i\ne 0)$ gives an oscillatory convection, which gives the dispersion relation of the form:

$$b_1{\left({\sigma }^2\right)}^2+b_2\left({\sigma }^2\right)+b_3=0$$

(37)

where:

$$\begin{array}{c}\hfill b_1=((Sr-1){\delta }^2+R_i){\lambda }_1-Va{\lambda }_2\end{array}$$

$$\begin{array}{ccc}\hfill b_2& =& (-BLe{\pi }^{2}Va{\alpha }^2{\lambda }_1^2(DuLe-LeSr+Le-1){\delta }^2+BL{e}^2{\pi }^2{R}_{i}Va{\alpha }^2{\lambda }_1^2)R{a}_s\hfill \\ & -& {\lambda }_1^2(LeSr+1)(DuLeSr-1){\delta }^8+{\lambda }_1(DuL{e}^{2}SrVa{\lambda }_2+L{e}^{2}SrVa{\lambda }_2\hfill \\ & -& Le{R}_{i}Sr{\lambda }_1+LeSrVa{\lambda }_2-R_i{\lambda }_1+Va{\lambda }_2){\delta }^6-L{e}^2(R_{i}SrVa{\lambda }_1{\lambda }_2\hfill \\ & -& SrVa{\lambda }_1+SrVa{\lambda }_2+Va{\lambda }_1-Va{\lambda }_2+Sr-1){\delta }^4\hfill \\ & +& L{e}^2(R_{i}Va{\lambda }_1-R_{i}Va{\lambda }_2-R_i+Va){\delta }^2\hfill \end{array}$$

$$\begin{array}{ccc}\hfill b_3& =& (-BLe{\pi }^{2}Va{\alpha }^2(DuLe-LeSr+Le-1){\delta }^2+BL{e}^2{\pi }^{2}Va{\alpha }^2{R}_i)R{a}_s\hfill \\ & +& (Va{\lambda }_1-Va{\lambda }_2-1)(LeSr+1)(DuLeSr-1){\delta }^8+(DuL{e}^{2}SrVa\hfill \\ & +& Le{R}_{i}SrVa{\lambda }_1-Le{R}_{i}SrVa{\lambda }_2+L{e}^{2}SrVa-Le{R}_{i}Sr+LeSrVa\hfill \\ & +& R_{i}Va{\lambda }_1-R_{i}Va{\lambda }_2-R_i+Va){\delta }^6-{\delta }^{4}L{e}^2{R}_{i}SrVa\hfill \end{array}$$

Now, the convection sets in via oscillatory mode if $\sigma >0$. To find an oscillatory marginal solution, we use Equation (36). For the value of ${\sigma }^2$, Equation (37) is solved and substituted in Equation (36) for a marginal oscillatory convection. However, Equation (37) is quadratic in ${\sigma }^2$, which can give more than one positive real root for some fixed values of parameters. Thus, first we determine the number of solutions of Equation (37). If there are no positive solutions, then no oscillatory instability is possible. If there are two positive solutions, then the minimum (over ${\alpha }^2$) of Equation (36) gives the oscillatory Rayleigh number with ${\sigma }^2$ given by Equation (37).

4. A Weak Nonlinear Theory

The aim of the linear theory discussed in Section 3 is to determine the critical values of the thermal Rayleigh number below which the system is stable, and above which it is unstable. Thus, when the flow with the thermal Rayleigh number is grater than critical values $Ra>R{a}_T^c$, the linear stability analysis will not give any idea about the flow field and become invalid, which is reasonable, because the nonlinear terms were neglected from the governing equations in the linear study. Linear stability theory does not give any information about the convection amplitudes, which is necessary to study the heat and mass transfer process.

In this section, a nonlinear stability analysis will be used to study the nonlinear effect on the flow field. The nonlinear study is performed involving a truncated Fourier series representations for stream function, temperature, and concentration. Nonlinear analysis with two terms of double Fourier series is enough to understand the physical mechanism for the rate of heat and mass transfer and to give more representation for nonlinear terms in the mathematical model with less mathematical analysis. In addition, it will be a step forward toward the study and generalization of the full nonlinear problem. Additionally, it is to be noted that the effect of nonlinearity is to distort the temperature and concentration fields through the interaction of $\psi$ and T, and $\psi$ and S, respectively. As a result, a component of the form $sin\left(2\pi z\right)$ will be generated [42].

For simplicity, the case of two-dimensional rolls is considered, and thus all physical quantities are made as independent of y. By eliminating pressure from Equation (16), introducing the stream function $u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x}$, using the transformations given in Equation (20), and dropping the asterisks, Equations (15)–(18) become:

$$\begin{array}{ccc}& & \left[\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\frac{1}{Va}\frac{\partial }{\partial t}{\nabla }^2+\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){\nabla }^2\right]\psi +\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \times \left[R{a}_T\frac{\partial T}{\partial x}-R{a}_s\frac{\partial S}{\partial x}\right]=0\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& & \left[\gamma \frac{\partial }{\partial t}-{\nabla }^2-R_i\right]T=\frac{\partial \psi }{\partial x}f\left(z\right)+\frac{\partial (\psi ,T)}{\partial (x,z)}+Du\frac{R{a}_s}{R{a}_T}{\nabla }^{2}S\hfill \end{array}$$

(39)

$$\begin{array}{ccc}& & \left[ϵ\frac{\partial }{\partial t}-\frac{1}{Le}{\nabla }^2\right]S=\frac{\partial \psi }{\partial x}b\left(z\right)+\frac{\partial (\psi ,S)}{\partial (x,z)}+S_r\frac{R{a}_T}{R{a}_s}{\nabla }^{2}T\hfill \end{array}$$

(40)

Now consider a minimal Fourier series with one term in the stream function and two terms in the temperature and concentration fields in order to obtain some effect of non-linearity as given below [32]:

$$\begin{array}{ccc}\hfill \psi & =& A_1\left(t\right)sin\left(\alpha x\right)sin\left(\pi z\right)\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill T& =& B_1\left(t\right)cos\left(\alpha x\right)sin\left(\pi z\right)+B_2\left(t\right)sin\left(2\pi z\right)\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill S& =& C_1\left(t\right)cos\left(\alpha x\right)sin\left(\pi z\right)+C_2\left(t\right)sin\left(2\pi z\right)\hfill \end{array}$$

(43)

where the amplitudes $A_1\left(t\right),B_1\left(t\right),B_2\left(t\right),C_1\left(t\right)$ and $C_2\left(t\right)$ are functions of time and are to be determined. Substituting Equations (41)–(43) into Equations (38)–(40) and equating the coefficients of like terms of the resulting equations, a system of nonlinear differential equations is obtained and is given as follows:

$$\begin{array}{ccc}& & \frac{d{A}_1\left(t\right)}{dt}=\frac{-Va}{\left({\delta }^2\right)(Va{\lambda }_2+1)}\left[\frac{d}{dt}\left(B_1\left(t\right)\right)R{a}_T\alpha {\lambda }_1-\frac{d}{dt}\left(C_1\left(t\right)\right)R{a}_s\alpha {\lambda }_1\right.\hfill \end{array}$$

$$\begin{array}{ccc}& & \left.+{\pi }^2{A}_1\left(t\right)+R{a}_T{B}_1\left(t\right)\alpha -R{a}_s{C}_1\left(t\right)\alpha +A_1\left(t\right){\alpha }^2\right]\hfill \\ & & \frac{d{B}_1\left(t\right)}{dt}=({\delta }^2-R_i)B_1\left(t\right)-A_1\left(t\right)\alpha (-B_2\left(t\right)\pi +2F)\hfill \end{array}$$

(44)

$$\begin{array}{ccc}& & -Du{C}_1\left(t\right)R{a}_s({\pi }^2-{\alpha }^2)/R{a}_T\hfill \end{array}$$

(45)

$$\begin{array}{ccc}& & \frac{d{B}_2\left(t\right)}{dt}=(4{\pi }^2-R_i)B_2\left(t\right)-\pi B_1\left(t\right)A_1\left(t\right)\alpha -\frac{4DuR{a}_s{C}_2\left(t\right){\pi }^2}{R{a}_T}\hfill \end{array}$$

(46)

$$\begin{array}{ccc}& & \frac{d{C}_1\left(t\right)}{dt}=\frac{-{\delta }^2{C}_1\left(t\right)}{Le}-\frac{SrR{a}_T{\delta }^2{B}_1\left(t\right)}{R{a}_s}+A_1\left(t\right)\alpha (-C_2\left(t\right)\pi +2B)\hfill \end{array}$$

(47)

$$\begin{array}{ccc}& & \frac{d{C}_2\left(t\right)}{dt}=A_1\left(t\right)\pi C_1\left(t\right)\alpha -\frac{4{\pi }^2{C}_2\left(t\right)}{Le}-\frac{4B_2\left(t\right){\pi }^{2}R{a}_{T}Sr}{R{a}_s}\hfill \end{array}$$

(48)

An analytical solution to the nonlinear system of autonomous differential equations is not suitable for the general time-dependent variables. A numerical solution for the nonlinear system of autonomous differential equations can be established. Then the amplitude functions $A_1\left(t\right)$, $B_1\left(t\right)$, $B_2\left(t\right)$, $C_1\left(t\right)$ and $C_2\left(t\right)$ can be determined. Therefore, the Nusselt number and Sherwood number can be obtained as functions of time.

The steady analysis is performed by setting the left-hand side of Equations (44)–(48) as equal to zero,

$$\begin{array}{ccc}& & -A_1{\delta }^2-B_{1}R{a}_T\alpha +C_{1}R{a}_s\alpha =0\hfill \end{array}$$

(49)

$$\begin{array}{ccc}& & (\alpha (-B_2\pi +2F)A_1-B_1({\delta }^2-R_i))R{a}_T-Du{C}_1{\delta }^{2}R{a}_s=0\hfill \end{array}$$

(50)

$$\begin{array}{ccc}& & ((-4{\pi }^2+R_i)B_2+\pi B_1{A}_1\alpha )R{a}_T-4DuR{a}_s{C}_2{\pi }^2=0\hfill \end{array}$$

(51)

$$\begin{array}{ccc}& & (Le\alpha (-C_2\pi +2B)A_1-C_1{\delta }^2)R{a}_s-Le{B}_{1}Sr{\delta }^{2}R{a}_T=0\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & -A_1{C}_{1}LeR{a}_s\alpha +4B_{2}Le\pi R{a}_{T}Sr+4C_2\pi R{a}_s=0\hfill \end{array}$$

(53)

with the consideration that all of the amplitudes are constants, since they are independent of time for the steady-state solution. We write $B_1,B_2,C_1$ and $C_2$ in terms of $A_1$ using Equations (49)–(53), and substitute these into Equation (49), with $x=A_1^2$ to obtain:

$$x_1{x}^2+x_{2}x+x_3=0$$

(54)

where:

$$\begin{array}{ccc}\hfill x_1& =& -L{e}^2{\pi }^2{\alpha }^4{\delta }^2\hfill \\ \hfill x_2& =& (-8{\alpha }^{4}FL{e}^2(Sr-1){\pi }^2-2{\alpha }^{4}FL{e}^2{R}_i)R{a}_T+8{\alpha }^{4}BLe{\pi }^2(DuLe-1)R{a}_s\hfill \\ & +& 4{\alpha }^2{\delta }^2(2DuL{e}^{2}Sr{\delta }^2+L{e}^2{\delta }^2-L{e}^2{R}_i+{\delta }^2){\pi }^2+{\alpha }^{2}L{e}^2{R}_i{\delta }^2(-{\delta }^2+R_i)\hfill \\ \hfill x_3& =& -\left(8(4DuLe{\pi }^{2}Sr-4{\pi }^2+R_i)\right)BLe{\alpha }^2(Du{\delta }^2+{\delta }^2-R_i)R{a}_s\hfill \\ & -& (DuLeSr{\delta }^4-2FLeR{a}_{T}Sr{\alpha }^2-2FR{a}_T{\alpha }^2-{\delta }^4+R_i{\delta }^2)\hfill \\ & \times & \left(4(4DuLe{\pi }^{2}Sr-4{\pi }^2+R_i)\right){\delta }^2\hfill \end{array}$$

(55)

The required root of Equation (54) is:

$$x^2=\frac{1}{2x_1}\left[-x_2+\sqrt{x_2^2-4x_1{x}_3}\right]$$

(56)

4.1. Steady Heat and Mass Transport

Heat and mass transport appears as a one of the most important issues in convection in fluids. That is because the onset of convection, as the Rayleigh number is increased, is more readily detected by its effect on the heat and mass transfer. In the basic state, heat and mass transfer are given by conduction only.

If H and J are the rate of heat and mass transfer per unit area, respectively, then [32]:

$$\begin{array}{ccc}\hfill H& =& -{\kappa }_T{〈\frac{\partial T_{total}}{\partial z}〉}_{z=0}-K_1{〈\frac{\partial S_{total}}{\partial z}〉}_{z=0}\hfill \end{array}$$

(57)

$$\begin{array}{ccc}\hfill J& =& -{\kappa }_s{〈\frac{\partial S_{total}}{\partial z}〉}_{z=0}-K_2{〈\frac{\partial T_{total}}{\partial z}〉}_{z=0}\hfill \end{array}$$

(58)

where the angular bracket represents the horizontal average. Additionally,

$$\begin{array}{ccc}\hfill T_{total}& =& T_0-\Delta T\frac{z}{d}+T(x,z,t)\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill S_{total}& =& S_0-\Delta S\frac{z}{d}+S(x,z,t)\hfill \end{array}$$

(60)

Substituting Equations (42) and (43) into Equations (59) and (60), respectively, and using the results in Equations (57) and (58), we obtain:

$$\begin{array}{c}\hfill Nu=1-\frac{2\pi B_2}{\sqrt{R_i/\phi }cot\left(\sqrt{R_i/\phi }\right)}\end{array}$$

(61)

$$\begin{array}{c}\hfill Sh=1-\frac{2\pi C_{2}R{a}_s}{-Le\sqrt{R_i/\phi }cot\left(\sqrt{R_i/\mathsf{\Omega }}\right)R{a}_{T}Sr+LeR{a}_{T}Sr+R{a}_s}\end{array}$$

(62)

where ${\displaystyle \phi =1-DuLeSr}$, $B_1$ is given in Appendix A, $B_2$ and $C_2$ are given in terms of $A_1$ as follows:

$$\begin{array}{ccc}\hfill B_2& =& \left[A_1{B}_{1}R{a}_T\alpha -\left(4\right(-A_{1}Le\alpha R{a}_{T}Sr{\delta }^2(4Le{\pi }^2-Le{R}_i+4{\pi }^2)B_1\right.\hfill \\ & +& 2A_1^{2}L{e}^2{\alpha }^{2}BR{a}_s(4{\pi }^2-R_i)\left)\right)Du\pi R{a}_s/(\pi R{a}_{s}L{e}^2{\alpha }^2(4{\pi }^2-R_i)A_1^2\hfill \\ & -& \left.4\pi R{a}_s{\delta }^2(4DuLe{\pi }^{2}Sr-4{\pi }^2+R_i))\right]\frac{\pi }{R{a}_T(4{\pi }^2-R_i)}\hfill \end{array}$$

(63)

$$\begin{array}{ccc}\hfill C_2& =& A_{1}Le\alpha (-4B_{1}Le{\pi }^{2}R{a}_{T}Sr{\delta }^2+8A_{1}BLe{\pi }^{2}R{a}_s\alpha +B_{1}LeR{a}_T{R}_{i}Sr{\delta }^2\hfill \\ & -& 4B_1{\pi }^{2}R{a}_{T}Sr{\delta }^2-2A_{1}BLeR{a}_s{R}_i\alpha )/(\pi R{a}_s(4A_1^{2}L{e}^2{\pi }^2{\alpha }^2\hfill \\ & -& A_1^{2}L{e}^2{R}_i{\alpha }^2-16DuLe{\pi }^{2}Sr{\delta }^2+16{\pi }^2{\delta }^2-4R_i{\delta }^2\left)\right)\hfill \end{array}$$

(64)

where:

$$\begin{array}{ccc}\hfill B_1& =& 8\alpha (\left(4(BDuLeRas-FR{a}_T)\right)(DuLe{\pi }^{2}Sr-{\pi }^2+(1/4)Ri){\delta }^2\hfill \\ & +& {\alpha }^2((BDuRas+FR{a}_T){\pi }^2-(1/4)R{a}_{T}FRi)L{e}^2{A}_1^2)A_1\hfill \\ & /& R{a}_T\left[(\left(16(DuLe{\pi }^{2}Sr-{\pi }^2+(1/4)Ri)\right)(DuLeSr-1){\delta }^4\right.\hfill \\ & +& (8{\alpha }^2((DuSr+1/2){\pi }^2-(1/8)Ri)A_1^{2}L{e}^2+16DuLe{\pi }^{2}RiSr\hfill \\ & +& \left.(4A_1^2{\alpha }^2-16Ri){\pi }^2+4R{i}^2){\delta }^2+{\alpha }^2((A_1^2{\alpha }^2-4Ri){\pi }^2+R{i}^2)L{e}^2{A}_1^2)\right]\hfill \end{array}$$

5. Results and Discussion

In this study, a porous layer saturated with viscoelastic fluid of Oldroyd type that is heated and salted from below was considered. The combined effects of cross-diffusion and internal heat source on the onset of stationary and oscillatory convection, as well as on the heat and mass transfer, were investigated using linear stability theory based on a normal mode technique, and weak nonlinear theory based on minimal representation of a truncated double Fourier series, respectively. Different values of governing parameters were selected to study instabilities and the convection process.

Equations (33) and (36) represent the critical thermal Rayleigh number for stationary and oscillatory modes, respectively. The onset of convection was investigated by varying the values of the following governing parameters: $R{a}_s,Du,Sr,R_i,{\lambda }_1,{\lambda }_2,Le,Va$. A comparison between these two modes was made to check whether the convection was initiated via the stationary or oscillatory mode. The results were presented graphically in terms of marginal stability curves, whereby the system is stable/unstable below/above the curve. In each case, the following parameter values were adopted: $R{a}_s=200,Du=0.05,$$Sr=0.005,R_i=1,{\lambda }_1=1.0,{\lambda }_2=0.3,Le=1,Va=10$, whereby one parameter was varied at a time to study its effect on system stability.

Figure 2 illustrates the effect of the Dufour parameter $\left(Du\right)$ on the thermal Rayleigh number $\left(R{a}_T\right)$ via stationary and oscillatory modes. Marginal stability curves showed that $R{a}_T$ increases with the increase in $Du$. Moreover, as can be seen from the graph that in the oscillatory mode, an increase in $Du$ increases the marginal curves as well as $R{a}_T$. Hence, the Dufour parameter has a stabilizing effect in both modes.

Figure 3 shows the effect of the Soret parameter $\left(Sr\right)$ on the stationary and oscillatory modes, indicating that an increase in $Sr$ decreases both the stationary marginal curves and $R{a}_T$. The Soret parameter also exerts a significant effect on the oscillatory convection. Specifically, as the value of $Sr$ increases, the critical thermal Rayleigh number decreases, indicating that the Soret parameter has a destabilizing effect on the system.

The effect of the internal Rayleigh number $\left(R_i\right)$ is presented in Figure 4. As can be seen from the graph, for smaller and greater values of $R_i$, the convection is initiated via the stationary and oscillatory mode, respectively. In both modes, $R_i$ has a destabilizing effect.

On the other hand, neither the relaxation time parameter $\left({\lambda }_1\right)$ nor the retardation time parameter $\left({\lambda }_2\right)$ affect the stationary thermal Rayleigh number, as indicated by Equation (33), which is independent of these two parameters. Nonetheless, a significant effect of ${\lambda }_1$ and ${\lambda }_2$ is observed on the onset of convection via the oscillatory mode. Figure 5 and Figure 6 illustrate the variation in $R{a}_T$ in terms of ${\lambda }_1$ and ${\lambda }_2$, respectively. As can be seen from these graphs, increasing ${\lambda }_1$ decreases marginal curves, whereas an increase in ${\lambda }_2$ increases both marginal curves and $R{a}_T$. Consequently, ${\lambda }_1$ has a destabilizing and ${\lambda }_2$ a stabilizing effect. It is noticeable that for ${\lambda }_1>{\lambda }_2$, the marginal curves for oscillatory convection decrease in comparison with the Newtonian case $({\lambda }_1={\lambda }_2)$. Therefore the system with ${\lambda }_1>{\lambda }_2$ is more unstable as compared to a Newtonian fluid.

Figure 7 depicts the influence of the Lewis number $\left(Le\right)$ on the onset of stationary and oscillatory convection, indicating that an increase in $Le$ increases the marginal curves for the stationary mode, thus stabilizing the system. On the other hand, as an increasing $Le$ has an opposite effect on the oscillatory convection, it destabilizes the system. This is because when $Le>1$, the heat diffusivity is higher than mass diffusivity. Thus, destabilizing the mass gradient enhances the onset of oscillatory convection.

The effect of the concentration Rayleigh number $\left(Ras\right)$ is depicted in Figure 8, indicating that, in both modes, an increase in $R{a}_s$ increases the marginal curves as well as the thermal Rayleigh number. In addition, for $R{a}_s\le 50$, the convection is initiated via the oscillatory mode, whereas the stationary mode applies to higher $R{a}_s$ values.

Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 depict the effects of all aforementioned governing parameters on the steady heat and mass transfer processes. In each case, one of the parameters varies and the others are set to the following values: $Ras=250,Ri=1,Sr=0.005,Du=0.05,Le=1$. The parameters ${\lambda }_1$, ${\lambda }_2$ and $Va$ are not considered in the steady heat and mass transport, since the expressions for the Nusselt number $\left(Nu\right)$ and Sherwood number $\left(Sh\right)$ given in Equations (61) and (62) are independent of these parameters. In extant literature, heat and mass transfer have been studied in terms of the Nusselt number and Sherwood number, respectively. The findings yielded indicate that variation in $Nu$ and $Sh$ with the thermal Rayleigh number represents the convection when $R{a}_T>R{a}_T^c$, which is the main concept behind the nonlinear stability theory indicating that the system is unstable.

Figure 9, Figure 10, Figure 11 and Figure 12 illustrate the effect of all of the studied parameters on the heat transfer process in terms of the Nusselt number, and Figure 9 illustrates the effect of the Dufour parameter. It is evident that variation in $Du$ does not exert a significant effect on the heat transfer process. However, as $Du$ increases, a small decrease in $Nu$ can be noted.

Figure 10 depicts the influence of $Sr$ on heat transfer, indicating that an increase in $Sr$ increases $Nu$, due to which heat transfer is enhanced. However, it is worth noting that, for $R{a}_T>600$, the effect of $Sr$ becomes insignificant even though $Nu$ continues to decline.

The effect of the internal Rayleigh number $\left(Ri\right)$ is illustrated in Figure 11, indicating that $Ri$ exerts a considerable influence on the heat transfer, as its increase increases $Nu$ and hence enhances the heat transfer process. The effect of $R{a}_s$ is depicted in Figure 12, demonstrating that heat transfer is suppressed as $R{a}_s$ values increase.

The influence of all aforementioned governing parameters on the mass transfer process is depicted in Figure 13, Figure 14, Figure 15 and Figure 16. The results are presented graphically in terms of the Sherwood number $\left(Sh\right)$. Specifically, Figure 13 shows that variation in $Du$ has a minimal effect on $Sh$ even though the critical thermal Rayleigh number increases as $Du$ increases, which in turn decreases $Sh$.

Figure 14 illustrates variations in $Sh$ as a result of changes in $Sr$. It can be seen that, for $R{a}_T<420$, $Sh$ increases as $Sr$ increases, which in turn enhances mass transfer. Conversely, for $RaT>420$, increasing $Sr$ decreases $Sh$ and thus suppresses mass transfer.

The influence of the internal Rayleigh number $\left(R_i\right)$ on the Sherwood number $\left(Sh\right)$ is depicted in Figure 15. Two regimes in the mass transfer process can be observed. Specifically, below a critical $R{a}_T$ value, increasing $R_i$ decreases $Sh$ and hence suppresses mass transfer, while at $R{a}_T$ above this value increasing $R_i$ increases $Sh$ and thus enhances mass transfer.

Figure 16 presents the effect of $R{a}_s$ on $Sh$, revealing that increasing $R{a}_s$ decreases $Sh$, which in turn suppresses the mass transfer process. When $R{a}_T$ is increased beyond $R{a}_T^c$, a sharp increas in the values of the Nusselt number and the Sherwood number is observed. However, when continuing to increase $R{a}_T$, the change in Nu and Sh becomes insignificant.

6. Conclusions

The combined effect of an internal heat source and cross diffusion on double-diffusive convection in a viscoelastic fluid-saturated porous layer was studied analytically using linear and nonlinear stability analysis. The porous layer is heated and salted from below. In each case, the following parameter values were adopted: $R{a}_s=200,Du=0.05,Sr=0.005,$$R_i=1,{\lambda }_1=1.0,{\lambda }_2=0.3,Le=1,Va=10$ whereby one parameter was varied at a time to study its effect on system stability. For the stationary and oscillatory modes, it was found that increasing the values of the Dufour parameter $\left(Du\right)$, retardation time parameter $\left({\lambda }_2\right)$, Lewis number $\left(Le\right)$, and concentration Rayleigh number $\left(R{a}_s\right)$ stabilizes the system. On the other hand, increasing the values of the Soret parameter $\left(Sr\right)$, internal Rayleigh number $\left(R_i\right)$, and relaxation time parameter $\left({\lambda }_1\right)$ destabilize the system. Steady heat and mass transport are independent of both the relaxation time parameter $\left({\lambda }_1\right)$ and retardation time parameter $\left({\lambda }_2\right)$. The effects of the Dufour parameter $\left(Du\right)$ and concentration Rayleigh number $\left(R{a}_s\right)$ are to suppress the heat and mass transfer processes as they increase. For high values of the thermal Rayleigh number $(R{a}_T>600)$, the effect of $Sr$ becomes insignificant. However, a small decrease in $Nu$ is still noticeable as $Sr$ increases. Increasing the internal Rayleigh number increases the Nusselt number as well as the heat transfer. On the other hand, there is a critical value for the thermal Rayleigh number $\left(R{a}_T\right)$, below which increasing the internal Rayleigh number $\left(R_i\right)$ decreases the Sherwood number $\left(Sh\right)$ and above which increasing the internal Rayleigh number $\left(R_i\right)$ increases the Sherwood number $\left(Sh\right)$. Additionally, there is a critical value for the thermal Rayleigh number $(R{a}_T<420)$, where below this value the heat and mass transfer are enhanced as the Soret parameter $\left(Sr\right)$ increased, while above it they are suppressed.