The Coriolis Effect on Thermal Convection in a Rotating Sparsely Packed Porous Layer in Presence of Cross-Diffusion

Abstract

The effect of rotation and cross-diffusion on convection in a horizontal sparsely packed porous layer in a thermally conducting fluid is studied using linear stability theory. The normal mode method is employed to formulate the eigenvalue problem for the given model. One-term Galerkin weighted residual method solves the eigenvalue problem for free-free boundaries. The eigenvalue problem is solved for rigid-free and rigid-rigid boundaries using the BVP4c routine in MATLAB R2020b. The critical values of the Rayleigh number and corresponding wave number for different prescribed values of other physical parameters are analyzed. It is observed that the Taylor number and Solutal Rayleigh number significantly influence the stability characteristics of the system. In contrast, the Soret parameter, Darcy number, Dufour parameter, and Lewis number destabilize the system. The critical values of wave number for different prescribed values of other physical parameters are also analyzed. It is found that critical wave number does not depend on the Soret parameter, Lewis number, Dufour parameter, and solutal Rayleigh number; hence critical wave number has no impact on the size of convection cells. Further critical wave number acts as an increasing function of Taylor number, so the size of convection cells decreases, and the size of convection cells increases because of Darcy number.

1. Introduction

Rotating convection in a sparsely packed porous layer, which is heated from below, has important applications in geophysics and geophysical fluid dynamics. Horton and Rogers [1] and Lapwood [2] were the first who did the experimental analysis of convective instabilities in a porous layer in the absence of rotation. The effect of rotation on Rayleigh-Benard convection (RBC) is studied by Tagare [3], and it is observed that while limiting the case of Prandtl number, Hopf bifurcation is not present. Gupta et al. [4] studied the RBC with rotation and magnetic field. Tagare et al. [5] have studied linear analysis and non-linear stability analysis of RBC of rotating fluids. A normal mode approach has been used to get the critical Rayleigh number for the modulated case by Om et al. [6]. They concluded that in the presence of modulation, the Taylor number could reduce the onset of convection. Novi et al. [7] deduced the effect of rotation with the tilted axis on RBC numerically. King et al. [8] studied the RBC for open and closed rotating cavities. Bhadauria et al. [9] studied the weakly non-linear analysis in a rotating porous medium and showed that as the Taylor number increases, the Nusselt number decreases.

The thermo-diffusion effect or Soret contributes to mass fluxes due to temperature gradients. Similarly, the diffusion-thermal effect or Dufour effect contributes to thermal energy flux due to concentration gradients. There may be a small effect of cross-diffusion, but when it is present in double-diffusive convections, they are more important as they have a significant impact on hydrodynamics stability compared to their contributions to the buoyancy of fluids. Dufour and Soret parameters have been found to appreciably affect the flow field in the mixed convection boundaries layer where the vertical surface is embedded in the porous layer. Dufour and Soret parameters are explained in many practical applications, such as geoscience, bioengineering, and chemical engineering. Many people studied convection with cross-diffusion. For example, Venkatesh and Pranesh [10] observed the effect of Dufour and Soret parameters on double-diffusive convection and deduced that the system is unstable if the Soret parameter increases. Kim et al. [11] studied the convection of nanofluid in the presence of Soret and Dufour effects. They observed that if Soret and Dufour parameters are included in the analysis for the nanofluid behavior, the system becomes unstable and heat transfer in nanofluid with the presence of the Soret effect is more significant than normal nanofluid. Hu et al. [12] studied the effect of the Soret parameter on Poissullie–Rayleigh–Benard convection. Stevens et al. [13] observed heat transport in rotating RBC. Convection in rotating fluids with Soret and Dufour parameters is explained by Duba et al. [14], Lewis number and Dufour parameter increase the heat transport strength. In contrast, the Soret parameter decreases the mass transport rate. Khalid et al. [15] studied the effect of Soret and Dufour parameters on magneto-convection. They used the Galerkin weighted residual method for solving the Eigenvalue problem.

The Soret parameter is useful in isotopes separation [16]. Niche et al. [17] investigated the Dufour and Soret effects on unsteady double-diffusive natural convection using the finite volume method. They concluded that with the Dufour coefficient, heat and mass transport increases. Gaikwad and Kamble [18] studied the effects of cross-diffusion on rotating anisotropic porous layer and deduced that the Dufour parameter could stabilize the system. Non-linear convection in couple stress fluids and Soret parameter is discussed by Malashetty et al. [19]. It has been deduced that heat and mass transfer of the system can be suppressed if the value of the Soret parameter increases. So, it is quite clear from the above discussion that we should not neglect the Soret and Dufour effects in double-diffusive convection.

Rotation in the sparsely packed porous layer plays a vital role in convective instability in geophysics, especially in the analysis of the interior part of Earth where molten liquid iron and other metals are electrically conducting. Rotating convection with purely internal heating on the horizontal porous medium is studied by Yadav et al. [20] and it is observed that while increasing the value of rotating parameter inhibit the onset of convection. Ravi et al. [21] studied the effect of cross-diffusion parameters on primary and secondary thermo-convective instabilities using the NWS equation and Lorentz equation. The effect of rotation on the viscoelastic fluid in the porous layer is examined by Rana and Kango [22], and they noticed that viscoelasticity and rotation increase the oscillatory mode. However, compressibility delays the onset of thermal instability. Later many researchers such as Malashetty and Swamy [23], Malashetty et al. [24,25,26], Mahajan and Sharma [27], Chand et al. [28], Yadav [29], Rana et al. [30,31], and Mikhailenko et al. [32], Mikhailenko and Sheremet [33] studied the effect of rotation with a different physical model. The system can be unstable because of differential diffusion and cross-diffusion. It plays a more important role as the system is significantly influenced by hydrodynamic stability in comparison to the buoyancy of fluids. The effect of the cross-diffusion parameter is observed in many real-life applications such as geosciences, bioengineering [34,35,36,37,38], geothermal heating from below for oceans, and chemical engineering. The Soret effect can be seen in isotope separation and a mixture of light molecular weight of gases.

The study of the Soret and Dufour effect on rotating convection in a sparsely packed porous layer for the realistic boundary conditions is of tremendous importance because it may be used as a fundamental mechanism for contaminant transport in groundwater, biochemical engineering, petroleum industry, oceanography, chemical engineering or oceans experience geothermal heating from below. It is observed that different kinds of boundaries play a very significant role in the Rayleigh–Benard problems in the case of the onset of convection. However, attention has not been given to such issues with cross-diffusion and rotation in the sparsely packed porous layer. So, there is sufficient space for further study. There are numerous practical applications of such study: petroleum industry, and geophysics including rotating, biochemical engineering, oceans experience geothermal heating from below, and many others.

This literature analysis shows that no work has been examined to study the cross-diffusion effect on thermohaline rotating convection in a sparsely packed porous layer for the realistic boundary conditions. Therefore, this paper studied the effect of rotation and cross-diffusion on convection in a sparsely packed porous medium for the realistic boundary conditions (rigid–rigid and rigid–free boundaries) and free–free boundaries and made a mathematical model of system with rotating fluid salted and heated from below and cross-diffusion factors. In Section 2, the relevant governing equations are discussed. Section 3 deals with linear stability analysis. The method of solution is described in Section 4. The results and discussion of obtained data are written in Section 5. Furthermore, finally, the conclusions are listed in Section 6.

2. Governing Equations

We consider thermally conducting fluid in a sparsely packed porous medium, which is placed between two infinitely parallel horizontal layers at z = 0 and z = d and kept rotating at a constant angular velocity, as shown in Figure 1. We used Cartesian coordinates (x, y, z) where x and y are horizontal coordinate and z is vertical coordinate. All the fluid physical properties are assumed to be constant, except for the density in the buoyancy term. The governing equations in the presence of cross-diffusion parameters for rotating fluid in a sparsely packed porous medium are considered such that the Oberbeck–Boussinesq approximation is valid.

First let us assume that

• Brinkman’s law holds;
• Viscous dissipation can be neglected;
• Local thermal equilibrium between solid phase and fluid phase holds;
• The Oberbeck–Boussinesq approximation can be applied.

The governing equations with Boussinesq approximation, are:

$$\nabla \cdot \overline{V}=0$$

(1)

$$\frac{{\mathsf{\rho }}_0}{\mathsf{\varphi }}\left(\frac{\partial \overline{V}}{\partial t}+\frac{1}{\mathsf{\varphi }}\left(\overline{V}\cdot \nabla \right)\overline{V}+2\left(\overline{\mathsf{\Omega }}\times \overline{V}\right)\right)=-\nabla P+\mathsf{\rho }\overline{g}-\frac{\mathsf{\mu }}{K}\overline{V}+{\mathsf{\mu }}_e{\nabla }^2\overline{V}$$

(2)

$$M\frac{\partial T}{\partial t}+\left(\overline{V}\cdot \nabla \right)T=k_{11}{\nabla }^{2}T+k_{12}{\nabla }^{2}S$$

(3)

$$\mathsf{\varphi }\frac{\partial S}{\partial t}+\left(\overline{V}\cdot \nabla \right)S=k_{22}{\nabla }^{2}S+k_{21}{\nabla }^{2}T$$

(4)

$$\mathsf{\rho }={\mathsf{\rho }}_0\left(1-\mathsf{\alpha }\left(T-T_0\right)+{\mathsf{\beta }}_s\left(S-S_0\right)\right)$$

(5)

where, $\overline{V}$ is velocity vector, ρ is density, ϕ is porosity, t is time, µ is viscosity of fluid, K is permeability of porous medium, μ_e is effective fluid viscosity, M is dimensionless heat capacity, g is acceleration due to gravity, T is temperature, P is pressure, k₁₁ is thermal diffusivity, k₁₂ and k₂₁ are Dufour and Soret coefficient, k₂₂ is mass diffusivity, S is solutal concentration, β_s is solutal expansion coefficient, α is thermal expansion coefficient. The density of fluid depends linearly on salinity and temperature. The thermal boundary conditions and solutal boundary conditions are given by:

$$\begin{array}{lll}T=T_0,& S=S_0,& \mathrm{on} z=d,\\ T=T_0+\Delta T,& S=S_0+\Delta S,& \mathrm{on} z=0.\end{array}\}$$

(6)

2.1. Basic State

Let us assume that basic state solution is time independent and only dependent on z−direction, hence basic state solutions can be written in form of

$${\overline{V}}_b=0,P=P_b\left(z\right),\mathsf{\rho }={\mathsf{\rho }}_b\left(z\right),S=S_b\left(z\right),T=T_b\left(z\right)$$

(7)

Putting the value of Equation (7) into Equations (1)–(4), we get

$$\nabla \cdot {\overline{V}}_b=0$$

(8)

$$\nabla P_b\left(z\right)-{\mathsf{\rho }}_b\overline{g}=0$$

(9)

$$k_{11}{\nabla }^2{T}_b\left(z\right)+k_{12}{\nabla }^2{S}_b\left(z\right)=0$$

(10)

$$k_{22}{\nabla }^2{S}_b\left(z\right)+k_{21}{\nabla }^2{T}_b\left(z\right)=0$$

(11)

Substituting Equation (6) in Equations (10) and (11), we get the solution of basic temperature state and basic concentration state are

$$T_b=T_0+\Delta T\left(1-\frac{z}{d}\right)$$

(12)

$$S_b=S_0+\Delta S\left(1-\frac{z}{d}\right)$$

(13)

2.2. Perturbed State

We now superpose small perturbations in the form

$$\begin{array}{l}\overline{V}={\overline{V}}_b+{\overline{V}}^{\prime },P=P_b\left(z\right)+P^{\prime },T=T_b\left(z\right)+T^{\prime },\\ S=S_b\left(z\right)+S^{\prime },\mathsf{\rho }={\mathsf{\rho }}_b\left(z\right)+{\mathsf{\rho }}^{\prime }\end{array}$$

(14)

Here primes indicate the perturbations. Now substituting Equation (14) into Equations (1)–(5), we obtain

$$\nabla \cdot {\overline{V}}^{\prime }=0$$

(15)

$${\mathsf{\rho }}_0\left(\frac{1}{\mathsf{\varphi }}\frac{\partial {\overline{V}}^{\prime }}{\partial t}+\frac{1}{{\mathsf{\varphi }}^2}\left({\overline{V}}^{\prime }\cdot \nabla \right){\overline{V}}^{\prime }+\frac{2}{\mathsf{\varphi }}\left(\overline{\mathsf{\Omega }}\times {\overline{V}}^{\prime }\right)\right)=-\nabla P^{\prime }+\left(\mathsf{\alpha }{T}^{\prime }-{\mathsf{\beta }}_S{S}^{\prime }\right){\mathsf{\rho }}_0\overline{g}-\frac{\mathsf{\mu }}{K}{\overline{V}}^{\prime }+{\mathsf{\mu }}_e{\nabla }^2{\overline{V}}^{\prime }$$

(16)

$$M\frac{\partial T^{\prime }}{\partial t}+\left({\overline{V}}^{\prime }\cdot \nabla \right)T^{\prime }+w^{\prime }\frac{\partial T_b}{\partial z}=k_{11}{\nabla }^2{T}^{\prime }+k_{12}{\nabla }^2{S}^{\prime }$$

(17)

$$\mathsf{\varphi }\frac{\partial S^{\prime }}{\partial t}+\left({\overline{V}}^{\prime }\cdot \nabla \right)S^{\prime }+w^{\prime }\frac{\partial S_b}{\partial z}=k_{22}{\nabla }^2{S}^{\prime }+k_{21}{\nabla }^2{T}^{\prime }$$

(18)

We are introducing dimensionless variables as follows

$$x^{*}=\frac{x^{\prime }}{d}, y^{*}=\frac{y^{\prime }}{d}, z^{*}=\frac{z^{\prime }}{d}, t^{*}=\frac{t}{\frac{M{d}^2}{ k_{11}}}, V^{*}=\frac{V^{\prime }}{\frac{k_{11}}{Md}}, T^{*}=\frac{T^{\prime }}{\Delta T}, S^{*}=\frac{S^{\prime }}{\Delta S}, P^{*}=\frac{P^{\prime }}{{\mathsf{\rho }}_0{k}_{11}^2{M}^{-2}{d}^{-2}}$$

The non-dimensional (after omitting the asterisk*) system of governing equations are

$$\left(\frac{1}{M^2\mathsf{\varphi }\mathit{Pr}}\frac{\partial }{\partial t}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}{\nabla }^2\right)\overline{V}-\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\left(\overline{V}\times {\overline{e}}_z\right)-\left(RaT-RsS\right){\overline{e}}_z=-\frac{\nabla P}{M^2\mathit{Pr}}-\frac{1}{M^2{\mathsf{\varphi }}^2\mathit{Pr}}\left(\overline{V}\cdot \nabla \right)\overline{V}$$

(19)

$$\left(\frac{\partial }{\partial t}-{\nabla }^2\right)T-Du\frac{Rs}{Ra}{\nabla }^{2}S-\frac{w}{M}=-\frac{1}{M}\left(\overline{V}\cdot \nabla \right)T$$

(20)

$$\left(\frac{\mathsf{\varphi }}{M}\frac{\partial }{\partial t}-\frac{1}{Le}{\nabla }^2\right)S-Sr\frac{Ra}{Rs}{\nabla }^{2}T-\frac{w}{M}=-\frac{1}{M}\left(\overline{V}\cdot \nabla \right)S$$

(21)

Here $\mathit{Pr}=\frac{\mathsf{\mu }}{{\mathsf{\rho }}_0{k}_{11}}$ is Prandtl number, $Ta={\left(\frac{2{\mathsf{\rho }}_0\mathsf{\Omega }{d}^2}{\mathsf{\mu }}\right)}^2$ is Taylor number, $Ra=\frac{{\mathsf{\rho }}_0\mathsf{\alpha }g\Delta T{d}^3}{\mathsf{\mu }{k}_{11}}$ is Rayleigh number, $Le=\frac{k_{11}}{k_{22}}$ is Lewis number, $Rs=\frac{{\mathsf{\rho }}_0{\mathsf{\beta }}_s\Delta S{d}^3}{\mathsf{\mu }{k}_{11}}$ is solutal Rayleigh number, $Da=\frac{K}{d^2}$ is Darcy number, $\mathsf{\Lambda }=\frac{{\mathsf{\mu }}_e}{\mathsf{\mu }},$ $Du=\frac{k_{12}\mathsf{\alpha }}{k_{11}{\mathsf{\beta }}_s}$ is Dufour parameter, $Sr=\frac{k_{21}{\mathsf{\beta }}_s}{k_{11}\mathsf{\alpha }}$ is Soret parameter, w is z-component of velocity vector, and ${\nabla }^2$ is Laplacian operator.

Operating curl on Equation (19), we get

$$\left(\frac{1}{M^2\mathsf{\varphi }\mathit{Pr}}\frac{\partial }{\partial t}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}{\nabla }^2\right)\nabla \times \overline{V}-\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\nabla \times \left(\overline{V}\times {\overline{e}}_z\right)-\nabla \times \left(RaT-RsS\right){\overline{e}}_z=-\frac{1}{M^2{\mathsf{\varphi }}^2\mathit{Pr}}\left(\nabla \times \left(\overline{V}\cdot \nabla \right)\overline{V}\right)$$

(22)

where, $\nabla \times \overline{V}=\mathsf{\omega }$ (vorticity).

Again applying curl to above resultant Equation (22), we obtain

$$\left(\frac{1}{M^2\mathsf{\varphi }\mathit{Pr}}\frac{\partial }{\partial t}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}{\nabla }^2\right) {\nabla }^2\overline{V}+\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\frac{\partial \mathsf{\omega }}{\partial z}+\nabla \times \nabla \times \left(RaT-RsS\right){\overline{e}}_z=\frac{1}{M^2{\mathsf{\varphi }}^2\mathit{Pr}}\left(\nabla \times \nabla \times \left(\overline{V}\cdot \nabla \right)\overline{V}\right)$$

(23)

Now collecting z-components of Equations (22) and (23) we get the equations as

$$\left(\frac{1}{M^2\mathsf{\varphi }\mathit{Pr}}\frac{\partial }{\partial t}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}{\nabla }^2\right){\mathsf{\omega }}_z-\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\frac{\partial w}{\partial z}=-\frac{1}{M^2{\mathsf{\varphi }}^2\mathit{Pr}}{\overline{e}}_z\cdot \left(\nabla \times \left(\overline{V}\cdot \nabla \right)\overline{V}\right)$$

(24)

$$\left(\frac{1}{M^2\mathsf{\varphi }\mathit{Pr}}\frac{\partial }{\partial t}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}{\nabla }^2\right){\nabla }^{2}w+\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\frac{\partial {\mathsf{\omega }}_z}{\partial z}-Ra{\nabla }_h^{2}T+Rs{\nabla }_h^{2}S=\frac{1}{M^2{\mathsf{\varphi }}^2\mathit{Pr}}{\overline{e}}_z\cdot \left(\nabla \times \nabla \times \left(\overline{V}\cdot \nabla \right)\overline{V}\right)$$

(25)

where ω_z and w are z-components of vorticity and velocity, respectively, and ${\nabla }_h^2$ is horizontal Laplacian operator.

3. Linear Stability Analysis

Let us introduce the normal modes by writing the perturbations in the form of

$$\left(w,{\mathsf{\omega }}_z,T,S\right)=\left(w\left(z\right),\mathsf{\Phi }\left(z\right),T\left(z\right),S\left(z\right)\right)e{ }^{\left(i\left(l_{1}x+l_{2}y\right)+\mathsf{\sigma }t\right)}$$

(26)

here l₁ and l₂ are wave numbers in directions x and y, respectively, and σ is the growth rate. We substitute Equation (26) into Equations (20), (21), (24) and (25) and neglect the nonlinear term, which gives us the following equations as

$$\left(\frac{\mathsf{\sigma }}{M^2\mathsf{\varphi }\cdot Pr}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}\left(D^2-a^2\right)\right)\left(D^2-a^2\right)w\left(z\right)+\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\frac{\partial \mathsf{\Phi }\left(z\right)}{\partial z}+a^{2}RaT\left(z\right)-a^{2}RsS\left(z\right)=0$$

(27)

$$\left(\frac{\mathsf{\sigma }}{M^2\mathsf{\varphi }\cdot Pr}+\frac{1}{M\cdot Da}-\frac{\mathsf{\Lambda }}{M}\left(D^2-a^2\right)\right)\mathsf{\Phi }\left(z\right)-\frac{\sqrt{Ta}}{M\mathsf{\varphi }}\frac{\partial w\left(z\right)}{\partial z}=0$$

(28)

$$\left(\mathsf{\sigma }-\left(D^2-a^2\right)\right)T\left(z\right)-\frac{w}{M}-Du\frac{Rs}{Ra}\left(D^2-a^2\right)S\left(z\right)=0$$

(29)

$$\left(\frac{\mathsf{\sigma }\mathsf{\varphi }}{M}-\frac{1}{Le}\left(D^2-a^2\right)\right)S\left(z\right)-\frac{w}{M}-Sr\frac{Ra}{Rs}\left(D^2-a^2\right)T\left(z\right)=0$$

(30)

where $D=\frac{d}{dz}$ and $a^2=l_1^2+l_2^2$.

In general σ is a complex number and in form of $\mathsf{\sigma }={\mathsf{\sigma }}_r+i{\mathsf{\sigma }}_i$. For σ_r < 0, system is always stable, and for σ_r > 0, system becomes unstable and for neutral stability of system σ_r = 0.

The above eigenvalue problem can be solved for the following boundary conditions:

Case (i) Free–Free boundary conditions:

$$w=D^{2}w=D\mathsf{\Phi }=T=S=0 \mathrm{at} z=0,1$$

(31)

Case (ii) Rigid–Free boundary conditions:

$$\begin{array}{ll}w=D^{2}w=\mathsf{\Phi }=T=S=0& \mathrm{at} z=0,\\ w=Dw=D\mathsf{\Phi }=T=S=0& \mathrm{at} z=1\end{array}\}$$

(32)

Case (iii) Rigid–Rigid boundary conditions:

$$w=Dw=\mathsf{\Phi }=T=S=0 \mathrm{at} z=0,1$$

(33)

4. Method of Solution

4.1. Exact Analytical Solution for Free–Free Boundary

Let us assume that solution is in the form of

$$w=w_0\sin \left(\pi z\right), \mathsf{\Phi }={\mathsf{\Phi }}_0\cos \left(\pi z\right), T=T_0\sin \left(\pi z\right), S=S_0\sin \left(\pi z\right)$$

(34)

which satisfies the free–free boundary condition Equation (31). On substituting Equation (34) in Equations (27)–(30) we get system of equation in matrix form as

$$\left[\begin{array}{cccc}\left(\frac{-\mathsf{\sigma }{\mathsf{\delta }}^2}{M^2\mathsf{\varphi }\cdot Pr}-\frac{{\mathsf{\delta }}^2}{M\cdot Da}-\frac{\mathsf{\Lambda }{\mathsf{\delta }}^4}{M}\right)& \frac{-\sqrt{Ta}}{M\mathsf{\varphi }}\pi & a^{2}Ra& -a^{2}Rs\\ \frac{-\sqrt{Ta}}{M\mathsf{\varphi }}\pi & \left(\frac{\mathsf{\sigma }}{M^2\mathsf{\varphi }\cdot Pr}+\frac{1}{M\cdot Da}+\frac{\mathsf{\Lambda }{\mathsf{\delta }}^2}{M}\right)& 0& 0\\ -\frac{1}{M}& 0& \mathsf{\sigma }+{\mathsf{\delta }}^2& Du\frac{Rs}{Ra}{\mathsf{\delta }}^2\\ -\frac{1}{M}& 0& Sr\frac{Ra}{Rs}{\mathsf{\delta }}^2& \left(\mathsf{\sigma }\frac{\mathsf{\varphi }}{M}+\frac{{\mathsf{\delta }}^2}{Le}\right)\end{array}\right]\times \left[\begin{array}{c}{w}_0\\ {\mathsf{\Phi }}_0\\ T_0\\ S_0\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

(35)

where ${\mathsf{\delta }}^2={\pi }^2+a^2$ and solving the above matrix for non-trivial solution, we get value of stationary Rayleigh number, oscillatory Rayleigh number and frequency of oscillation ω, as

$$R{a}_{sc}=\frac{{\chi }_{1}Le}{{\chi }_2}Rs+\frac{{\mathsf{\delta }}_{sc}^2{\chi }_3}{a_{sc}^2{\chi }_2}\left({\mathsf{\delta }}_{sc}^4\mathsf{\Lambda }+\frac{{\mathsf{\delta }}_{sc}^2}{Da}+\frac{M{\pi }^{2}Ta}{{\mathsf{\varphi }}^2{\chi }_4}\right)$$

(36)

$$R{a}_{oc}=\frac{I_1}{K},{\mathsf{\omega }}^2=\frac{-I_3+\sqrt{I_3^2-4I_2{I}_4}}{I_2}$$

(37)

where, expressions of Equations (36) and (37) are given in the Appendix A.

When rotation is not present (i.e., Ta = 0) and limiting case of sparsely packed porous medium in Equation (36), we obtain

$$R{a}_{sc}=\frac{{\chi }_{1}Le}{{\chi }_2}Rs+\frac{{\mathsf{\delta }}_{sc}^6{\chi }_3}{a_{sc}^2{\chi }_2}$$

(38)

which is the same as obtained by Ravi et al. [39].

In the absence of cross-diffusion i.e., Soret parameter Sr = 0, Dufour parameter Du = 0 and Solutal Rayleigh number Rs = 0 in Equation (36), we obtain

$$R{a}_{sc}=\frac{M}{a_{sc}^2}\left(\frac{{\mathsf{\delta }}_{sc}^4{\left(\frac{1}{MDa}+\frac{\mathsf{\Lambda }}{M}{\mathsf{\delta }}_{sc}^2\right)}^2+{\mathsf{\delta }}_{sc}^2{\pi }^2\frac{Ta}{{\mathsf{\varphi }}^2}}{\left(\frac{1}{MDa}+{\mathsf{\delta }}_{sc}^2\frac{\mathsf{\Lambda }}{M}\right)}\right)$$

(39)

which is the same as obtained by Babu et al. [40].

Further, when rotation is not present (i.e., Ta = 0) and limiting case of sparsely packed porous layer in Equation (39) we get the classical result of Rayleigh number given in Chandrasekhar [41].

$$R{a}_{sc}=\frac{{\mathsf{\delta }}_{sc}^6}{a_{sc}^2}$$

(40)

4.2. Numerical Solution for Rigid–Free and Rigid–Rigid Boundaries

The analytical solution is not possible for the eigenvalue problem for the realistic boundary conditions. So, we had solved eigenvalue problem numerically. Therefore, we use the bvp4c routine in MATLAB R2020b to solve the eigenvalue problem for rigid–free and rigid–rigid boundaries. We reduced the system of higher order ordinary differential Equations (27)–(30) into a system of first-order ordinary differential equations and put σ = 0. For non-trivial solutions and determining the eigenvalue Ra, we used the normalization condition $w^{\prime }\left(0\right)=1$. To calculate the value of critical Rayleigh number Ra_c and corresponding wave number a_c, we have used the indexmin command in MATLAB R2020b. To gain higher-order accuracy, the absolute and relative tolerance has been taken as 10⁻⁹ and 10⁻⁶ respectively.

To validate our solution method, we have compared the obtained results with those existing in the literature. The current problem can be changed to Chandrasekhar [41] in the absence of a cross-diffusion effect and sparsely packed porous medium.

Table 1. Comparison between the present article’s results with existing results for Rigid–Free boundary conditions for limiting case of sparsely packed porous medium.

Ta	Chandrasekhar [41]		Present Study
	Rigid–Free		Rigid–Free
	Rac	ac	Rac	ac
0	1100.6	2.68	1100.659	2.687
6.25	1108.5	2.68	1107.734	2.697
31.25	1135.9	2.70	1135.688	2.737
62.5	1169.5	2.79	1168.727	2.798
187.5	1291.7	2.97	1290.808	2.975
625	1637.6	3.40	1637.408	3.394
1875	2360.3	4.00	2358.337	4.006
6250	4047.7	4.92	4044.944	4.930

and

Table 2. Comparison between the present article’s results with existing results for Rigid–rigid boundary conditions for limiting case of sparsely packed porous medium.

Ta	Chandrasekhar [41]		Present Study
	Rigid–Rigid		Rigid–Rigid
	Rac	ac	Rac	ac
0	1707.7	3.11	1707.767	3.115
10	1713.0	3.10	1712.679	3.121
100	1756.6	3.15	1756.352	3.162
500	1940.3	3.30	1940.204	3.318
1000	2151.7	3.50	2151.345	3.484
2000	2530.5	3.75	2530.129	3.746
5000	3469.2	4.25	3468.500	4.263
10,000	4713.1	4.80	4712.047	4.788
30,000	8326.4	5.80	8324.614	5.797

show an excellent agreement of our numerical results with the critical Rayleigh number, Ra_c, and corresponding critical wave number, a_c, given in Chandrasekhar [41].

5. Results and Discussion

The numerical results and discussion are presented in this section. In the present analysis, the linear stability analysis has modeled rotating RBC of a sparsely packed porous medium in the presence of the cross-diffusion effect. Most of the previous studies are based on the Darcy model. Therefore, they are relatively well packed with low permeability of porous layer. However, in the medium of sparsely packed, we cannot apply Darcy’s law in its usual form. In a sparsely packed medium, they involve big void spaces giving rise to viscous shear and the Darcy resistance. Physically, one can consider the application of convection in the case of disposal of nuclear waste material in the underground soil below the surface of the sedimentary layer, where the porosity is expected to be 40%−50%, so Darcy’s law in its present form is not correct form to explain the flow field as medium having big void spaces. One-term Galerkin weighted residual method solves the eigenvalue problem with free–free boundaries. Eigenvalue problems with rigid–rigid and free–rigid boundaries are solved with the help of bvp4c in MATLAB R2020b. The linear instability threshold parameters consisting of Rayleigh number, Ra. and corresponding wave number, a, depend on Soret parameter, Sr., Dufour parameter, Du., Solutal Rayleigh number, Rs., Taylor number, Ta., Lewis number, Le. and Darcy number, Da. are shown in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11, and

Table 3. Critical values of Ra and a for Da = 0.01, Du = 0.02, Rs = 200, Le = 5, Λ = 8, M = 0.9, ϕ = 0.9.

Ta	Sr	Free–Free		Rigid–Free		Rigid–Rigid
		Rac	ac	Rac	ac	Rac	ac
10	0.2	5174.744	2.512	6968.268	2.824	9372.595	3.156
	0.4	3386.363	2.512	4557.644	2.824	6127.817	3.156
	0.6	2492.173	2.512	3352.333	2.824	4505.428	3.156
	0.8	1955.659	2.512	2629.146	2.824	3531.995	3.156
100	0.2	5180.744	2.514	6972.604	2.825	9375.175	3.156
	0.4	3390.282	2.514	4560.476	2.825	6129.502	3.156
	0.6	2495.051	2.514	3354.412	2.825	4506.665	3.156
	0.8	1957.912	2.514	2630.774	2.825	3532.964	3.156
1000	0.2	5240.303	2.535	7015.794	2.836	9400.934	3.161
	0.4	3429.178	2.535	4588.682	2.836	6146.324	3.161
	0.6	2523.615	2.535	3375.126	2.836	4519.019	3.161
	0.8	1980.277	2.535	2646.992	2.836	3542.636	3.161

,

Table 4. Critical values of Ra and a for Da = 0.01, Sr = 0.2, Rs = 200, Le = 5, Λ = 8, M = 0.9, ϕ = 0.9.

Ta	Du	Free–Free		Rigid–Free		Rigid–Rigid
		Rac	ac	Rac	ac	Rac	ac
10	0.02	5174.744	2.512	6968.268	2.824	9372.595	3.156
	0.04	5089.545	2.512	6846.467	2.824	9201.726	3.156
	0.06	5004.346	2.512	6724.665	2.824	9030.857	3.156
	0.08	4919.147	2.512	6602.864	2.824	8859.987	3.156
100	0.02	5180.744	2.514	6972.604	2.825	9375.175	3.156
	0.04	5095.423	2.514	6850.714	2.825	9204.253	3.156
	0.06	5010.101	2.514	6728.824	2.825	9033.331	3.156
	0.08	4924.780	2.514	6606.934	2.825	8862.409	3.156
1000	0.02	5240.303	2.535	7015.794	2.836	9400.934	3.161
	0.04	5153.767	2.535	6893.023	2.836	9229.486	3.161
	0.06	5067.230	2.535	6770.252	2.836	9058.038	3.161
	0.08	4980.693	2.535	6647.480	2.836	8886.591	3.161

,

Table 5. Critical values of Ra and a for Da = 0.01, Du = 0.02, Sr = 0.2, Le = 5, Λ = 8, M = 0.9, ϕ = 0.9.

Ta	Rs	Free–Free		Rigid–Free		Rigid–Rigid
		Rac	ac	Rac	ac	Rac	ac
10	100	4919.744	2.512	6713.268	2.824	9117.595	3.156
	150	5047.244	2.512	6840.768	2.824	9245.095	3.156
	200	5174.744	2.512	6968.268	2.824	9372.595	3.156
	250	5302.244	2.512	7095.768	2.824	9500.095	3.156
100	100	4925.744	2.514	6717.604	2.825	9120.175	3.156
	150	5053.244	2.514	6845.104	2.825	9247.675	3.156
	200	5180.744	2.514	6972.604	2.825	9375.175	3.156
	250	5308.244	2.514	7100.104	2.825	9502.675	3.156
1000	100	4985.303	2.535	6760.794	2.836	9145.934	3.161
	150	5112.803	2.535	6888.294	2.836	9273.434	3.161
	200	5240.303	2.535	7015.794	2.836	9400.934	3.161
	250	5367.803	2.535	7143.294	2.836	9528.434	3.161

,

Table 6. Critical values of Ra and a for Da = 0.01, Du = 0.02, Rs = 200, Sr = 0.2, Λ = 8, M = 0.9, ϕ = 0.9.

Ta	Le	Free–Free		Rigid–Free		Rigid–Rigid
		Rac	ac	Rac	ac	Rac	ac
10	3	6261.029	2.512	8521.236	2.824	11,551.179	3.156
	4	5657.537	2.512	7658.476	2.824	10,340.854	3.156
	5	5174.744	2.512	6968.268	2.824	9372.595	3.156
	6	4779.731	2.512	6403.553	2.824	8580.383	3.156
100	3	6268.591	2.514	8526.700	2.825	11,554.429	3.156
	4	5664.231	2.514	7663.313	2.825	10,343.732	3.156
	5	5180.744	2.514	6972.604	2.825	9375.175	3.156
	6	4785.163	2.514	6407.478	2.825	8582.718	3.156
1000	3	6343.648	2.535	8581.129	2.836	11,586.891	3.161
	4	5730.679	2.535	7711.498	2.836	10,372.470	3.161
	5	5240.303	2.535	7015.794	2.836	9400.934	3.161
	6	4839.087	2.535	6446.582	2.836	8606.040	3.161

and

Table 7. Critical values of Ra and a for Sr = 0.2, Du = 0.02, Rs = 200, Le = 5, Λ = 8, M = 0.9, ϕ = 0.9.

Ta	Da	Free–Free		Rigid–Free		Rigid–Rigid
		Rac	ac	Rac	ac	Rac	ac
10	0.001	22,816.666	2.955	25,371.434	3.094	28,350.473	3.237
	0.01	5174.744	2.512	6968.268	2.824	9372.595	3.156
	0.1	3304.991	2.264	5042.815	2.705	7422.587	3.126
	1	3110.685	2.227	4847.188	2.684	7226.586	3.116
100	0.001	22,817.662	2.955	25,372.359	3.094	28,351.296	3.237
	0.01	5180.744	2.514	6972.604	2.825	9375.175	3.156
	0.1	3317.022	2.273	5049.936	2.705	7425.830	3.126
	1	3124.063	2.236	4854.861	2.691	7229.924	3.116
1000	0.001	6343.648	2.535	8581.129	3.094	28,359.530	3.237
	0.01	5730.679	2.535	7711.498	2.836	9400.934	3.161
	0.1	5240.303	2.535	7015.794	2.733	7458.252	3.126
	1	4839.087	2.535	6446.582	2.719	7263.209	3.126

.

Let us fix the values Λ = 8 (see Nield and Bejan [41]), M = 0.9 and ϕ = 0.9. Table 3 shows the critical values of Rayleigh number and wave number for the different values of Ta and Sr and for the fixed values of Du = 0.02, Rs = 200, Le = 5 and Da = 0.01.

Thermo-diffusion effect or Soret effect contributes to mass fluxes, and diffusion-thermal effect or Dufour effect contribute to thermal energy flux, which has a direct impact on the stability of the system of rotating fluid, so we draw the graph between critical Rayleigh number and Taylor number for different values of Soret and Dufour numbers. Similarly, the Lewis number is also an important parameter because it relates to thermal diffusivity and mass diffusivity. For example, salt diffuses 100 times less in the ocean than does heat diffuses. Therefore, mass diffusivity plays a significant role in governing the system, so we plotted the graph for critical Rayleigh number for different values of Lewis number. Solutal Rayleigh number depends on the difference between concentrations of two plates. Therefore, the Solutal Rayleigh number also plays a vital role on the onset of convection. Since the Darcy number depends on the geometry of the parallel plate, the Darcy number affects the onset of convection because the geometry of the similar plate plays its role in governing the equations. Therefore, a graph is plotted for the critical Rayleigh number and the different values of Darcy number. Graphical representation of these values has been given in Figure 2 and Figure 3. Figure 2 provides a visual representation of Rac versus Ta for the fixed values of all other given parameters. From this figure, it is clear that the critical Rayleigh number increases as Ta increases, and hence the Taylor number has a stabilizing effect on the system. The same results were obtained for rotating parameters from Yadav et al. [42]. This effect can be attributed to the fact that Coriolis force, which arises due to rotation, enhances the horizontal motions, simultaneously limiting the vertical activity in the system. Overall, this leads to suppressed convective motion hence stabilizing the system. Furthermore, the behavior of the Soret parameter, i.e., $Sr=\frac{k_{21}{\mathsf{\beta }}_s}{k_{11}\mathsf{\alpha }}$ on the onset of convection, is made clear in Figure 2. From Figure 2, it is clear that the critical Rayleigh number decreases as Sr increases, which means that Sr has a destabilizing effect. This happens because the Soret coefficient k₂₁ or thermo-diffusion coefficient increases for the fixed value of thermal diffusivity. So, the Soret number makes a stronger disturbance in the system, and hence it becomes unstable. So, the onset of convection is in advance because of the Soret parameter. Figure 3 depicts a relation between critical wave number, a_c, and Ta. One can observe from this figure that a_c is an increasing function of Ta. Thus we conclude that the size of convection cells is reduced. However, critical wave number a_c is not changed with an increase in the Soret number, so critical wave number is not dependent on Soret number. These results are confirmed by Yadav et al. [20], where the rotating parameter, Ta, has a stabilizing effect. Therefore, the onset of convection is delayed.

In Table 4, the critical values of Ra and a are represented for the different values of Taylor number and Du and the fixed values of Sr = 0.2, Rs = 200, Le = 5, and Da = 0.01. Visual representation of these values is given in Figure 4 and Figure 5. The critical value of Ra increases as Ta increases and decreases as Du increases (as shown in Figure 4). Hence, an increase in Ta causes stabilization of the system, whereas the behavior of the Dufour parameter, i.e., $Du=\frac{k_{12}\mathsf{\alpha }}{k_{11}{\mathsf{\beta }}_s}$ on the onset of convection, is made clear in Figure 4. From Figure 4, it is clear that the critical Rayleigh number decreases as Du increases, which means that Du has a destabilizing effect. This happens because the Dufour coefficient k₁₂ or diffusion-thermal coefficient increases for the fixed value of thermal diffusivity. So, the Dufour number makes a stronger disturbance in the system, and hence it becomes unstable. This indicates that the onset of convection enhances due to Du. Parameter a_c is an increasing function of Ta (see Figure 5). Thus we conclude that the size of convection cells reduces. When the Taylor number increases, Yadav et al. [43] obtain the same type of observation for critical wavenumbers were the size of convection cells contracts. However, a_c is not changed with an increase in Du, which means that a_c does not depend on Du, so the size of convection cells is independent on the Dufour parameter.

The graph of Ra_c and a_c for different values of Ta and Rs is shown in Table 5 or Figure 6 and Figure 7. Figure 6 shows that the critical value of Ra increases as Rs increases; therefore, the onset of convection is delayed with the solutal Rayleigh number, hence system is stable. This may be interpreted from the definition of solutal Rayleigh number that the value of solutal Rayleigh number directly depends upon the difference between concentration at upper and lower plates. Hence, the difference between concentrations plays a role in determining the system’s stability; therefore, the solutal Rayleigh number helps for delaying the onset of convection. As Ta increases, the critical value of Ra increases (see Figure 6). Hence, Ta also has a stabilizing effect on the system. As we observed in Figure 3 and Figure 5, a_c is an increasing function of Ta (see Figure 7), whereas critical wave number a_c is not changed as Rs increases. So, the size of convection cells does not depend upon the solutal Rayleigh number.

The variation of critical values of Ra and wave number as a function of Taylor number for different values of Le is given in Table 6 or Figure 8 and Figure 9. Figure 8 indicates that Ra_c decreases as Le increases. This shows that Le causes a strong destabilization in the system. The interplay may be explained as the ratio of thermal diffusivity to mass diffusivity is defined as Lewis number. It comes into a fixture when we characterize the fluid flow where heat and mass transfer simultaneously happen. For example, Le >> 1 for fluid flow and gases around 1. For a given thermal diffusivity, a higher Lewis number corresponds to a lower molecular diffusivity, hence the critical Rayleigh number decreases, and therefore Le shows destabilization in the system. Alternately, the fixed value of molecular diffusivity Lewis number directly depends on thermal diffusivity. So thermal diffusivity causes a destabilization effect. Taylor number can stabilize the system as expected from Figure 2, Figure 4 and Figure 6. The critical value of a is not changed as the value of Le increases (see Figure 9). So the size of convection cells is independent on the Lewis number.

Variation of Ra_c and a_c is defined as a function of Ta for the different values of Da as given in Table 7 or Figure 10 and Figure 11. The Darcy number is defined as the permeability ratio of the porous layer to the square of the distance between two parallel plates. The Darcy number increases when permeability is more than the square of the distance. From Figure 10, we observed an increase in the value of Da tends to decrease in the value of Ra_c, so the system is unstable, and the onset of convection enhances for the Darcy number. Moreover, as we observed in Figure 2, Figure 4, Figure 6 and Figure 8, Taylor number Ta can stabilize the system (see Figure 10). Figure 11 shows that the critical value of wave number decreases as Da increases, so the size of convection cells increases when the value of Da increases. These results are the same as the result of Yadav et al. [20], where the size of convection cells increases with an increase in the value of Darcy number Da.

Figure 12 is plotted between critical Rayleigh number Ra_c and solutal Rayleigh number Rs with fixed value of Du = 0.02, Sr = 0.2, Le = 5, Da = 0.01, Λ = 8, ϕ = 0.9, M = 0.9 and Ta = 10. The impact of the solutal Rayleigh number on the onset of convection is plotted in Figure 12. It shows that critical Rayleigh number Ra_c increases when the solutal Rayleigh number Rs increases, so the system is stabilized with an increase in solutal Rayleigh number Rs. This may be interpreted from the definition of solutal Rayleigh number that the value of solutal Rayleigh number directly depends upon the difference between concentration at upper and lower plates. Hence, the difference between concentrations plays a role in determining the system’s stability. The system becomes more stable if the solutal Rayleigh number increases.

6. Conclusions

This paper shows the effect of cross-diffusion on rotating convection in a sparsely packed porous medium. One-term Galerkin weighted residual method solves the eigenvalue problem for free–free boundaries. The eigenvalue problem is solved for rigid–free and rigid–rigid boundaries using the BVP4c routine in MATLAB R2020b. Using the results discussed in the previous section, we now draw general conclusions of the problem:

• The solutal Rayleigh number and Taylor number have stabilizing effect;
• The Soret number, Lewis number, Dufour number and Darcy number have destabilizing influence on the system;
• The critical wave number is an increasing function of Taylor number, so the size of convection cells decreases, and critical wave number is a decreasing function of the Darcy number; hence the size of convection cells increases;
• The critical wave number does not depend on the Soret parameter, Lewis number, Dufour parameter and solutal Rayleigh number;
• From the above obtained results, the system of rigid–rigid boundary is found to be most stable whereas the free–free boundary is found to be least stable.