Onset of Linear and Nonlinear Thermosolutal Convection with Soret and Dufour Effects in a Porous Collector under a Uniform Magnetic Field

Abstract

The present paper reports on an analytical and numerical study of combined Soret and Dufour effects on thermosolutal convection in a horizontal porous cavity saturated with an electrically conducting binary fluid under a magnetic field. The horizontal walls of the system are subject to vertical uniform fluxes of heat and mass, whereas the vertical walls are assumed to be adiabatic and impermeable. The main governing parameters of the problem are the Rayleigh, the Hartmann, the Soret, the Dufour and the Lewis numbers, the buoyancy ratio, the enclosure aspect ratio, and the normalized porosity of the porous medium. An asymptotic parallel flow approximation is applied to determine the onset of subcritical nonlinear convection. In addition, a linear stability analysis is performed to predict explicitly the thresholds for the onset of stationary, overstable and oscillatory convection, and the Hopf bifurcation as functions of the governing parameters. The combined effect of a magnetic field, Soret and Dufour parameters have a noticeable influence on the intensity of the convective flow, the heat and mass transfer rates, and the thresholds of linear convection. It is found that the imposition of a magnetic field delays the onset of convection and its intensification can lead to the total suppression of the convective currents. The heat transfer rate increases with the Dufour number and decreases with the Soret number and vice versa for the mass transfer rate.

1. Introduction

The problem of thermosolutal natural convection in enclosures filled with saturated porous media with or without the presence of a magnetic field had been the subject of numerous recent and past studies. The interest rose from the occurrence of the phenomenon in many engineering applications such as geothermal energy, diffusion of moisture in fibrous insulations, food processing, drying processes, spread of pollutants in soil, solar ponds, crystal growth in fluids, and metal casting [1,2,3,4,5].

Thermosolutal natural convection was studied widely in square or rectangular cavities at different thermal and solutal boundary conditions as reported in the literature [6,7,8,9,10,11,12,13,14,15], where the authors examined numerically and performed scale analyses of the effect of the permeability ratio and other governing parameters on the flow structure, isotherms, isoconcentrations, and on heat and mass transfer rates in two- and three-dimensional double diffusive convection generated by combined thermal and solutal gradients in a bi-layered porous enclosures. The authors showed that the permeability of the two porous layers and the variation of the ratios of thermal conductivities and mass diffusivities had a significant influence on the double diffusive convection, the flow structure, and on the heat and mass transfer rates. The onset of double-diffusive convection in a shallow enclosure subject to vertical gradients of heat and solute was investigated analytically and numerically by Mamou and Vasseur [16]. On the basis of the parallel flow approximation, the onset of convection from the rest state was subcritical for opposing flow when the Lewis number was greater than unity. Multiple solutions and travelling waves were found to exist. Kalla et al. [17] presented analytical and numerical studies of double-diffusive natural convection in a rectangular enclosure filled with fluid or porous medium for both aiding and opposing flow cases. The authors demonstrated the existence of multiple convective solutions for a wide range of the governing parameters. A stability analysis of double-diffusion convection in a horizontal porous enclosure subject to different thermal and solutal boundary conditions was performed by Mamou [18]. It was found that an increase in the porosity of the porous media delayed the appearance of oscillatory flows. Rebhi et al. [19,20] recently carried out analytical and numerical studies of convective flows induced in vertical and horizontal porous cavities filled with a binary mixture using the Dupuit–Darcy model. Studying the linear stability of the parallel flow approximation, the authors predicted explicitly and implicitly the onset of convective motion and the triggering point of Hopf bifurcation characterizing the transition from convective steady to oscillatory states. The form drag was found to have a significant influence on all the critical conditions, and on the heat and mass transfer rates. In addition, the authors demonstrated numerically that multiple solutions were possible for a given set of the governing parameters.

Changhao and Payne [21] presented a mathematical study on the thermosolutal convection in a porous medium where the Darcy model was employed. The authors established a continuous dependence of the flow solution on the Soret effect. Theoretical and numerical analysis of Soret-driven convection in a horizontal porous layer saturated by an n-component mixture was investigated by Mutshler and Mojtabi [22]. In the first part, an analytical and numerical study of the onset of Soret driven convection was presented. The study was based on the classical Darcy–Boussinesq equations, which admitted a mechanical solution associated with the pure double-diffusive regime. In the second part, the analytical solution for the unicellular flow was obtained, and the separation was expressed in terms of the Lewis number, the separation ratio, the cross-diffusion coefficient and the Rayleigh number. Benano-Molly et al. [23] investigated the effect of Soret coefficient within a rectangular porous medium saturated by a binary fluid mixture when the thermal and solutal buoyancy forces were opposing each other. It was shown that, when the solutal buoyancy force ratio was negligible, the theory represented well the solute behavior. Mansour et al. [24] studied the Soret effect on double diffusive convection and on heat and mass transfer rates in a square cavity. The heat transfer rate was found to be significantly affected by the Soret effect.

Furthermore, Joly et al. [25] presented an analytical and numerical study of the influence of the Soret effect on the onset of convection in a vertical porous cavity saturated with a binary mixture. The vertical walls were subjected to uniform heat fluxes. The Brinkman-extended Darcy model was used to solve the governing equations. The results indicated that the critical Rayleigh number depended strongly upon of the control parameters such as the aspect ratio of the cavity, the Darcy and the Lewis numbers. Gaikwad et al. [26] made an analysis of thermosolutal convection in a horizontal anisotropic saturated porous layer with Soret effect. The heat and mass transfer rates increased with the anisotropy parameters and the Lewis number; in addition, the heat transfer increased with the negative Soret parameter while it decreased with the positive one. A reverse trend was found for the mass transfer rate. Malashetty et al. [27] presented a numerical investigation of thermosolutal convection in a porous layer saturated by a couple–stress fluid with Soret effect. Linear and weak nonlinear stability analyses were performed. The heat and mass transfer rates decreased with increasing the Taylor number and the couple–stress parameter, while both increased with increasing the solute Rayleigh number. The heat transfer rate decreased with increasing the Lewis number while the mass transfer rate increased significantly.

In addition, Mojtabi et al. [28] carried out an analytical and numerical analysis of the species separation in a parallelepipedic cell filled by a binary mixture. Constant velocity was imposed on the top and/or the bottom plate of the cavity. The fixed velocity was obtained from the superposition of the flow generated by the velocity of the wall under weightless condition and thermoconvective flow under gravity only. 2D and 3D direct numerical simulations were performed using a finite difference method in order to corroborate the analytical results. The authors showed that the effective species separation admitted a partial optimum as a function of the velocity ratio of the moving walls and the wall velocity. The aspect ratio in the y-direction had an effect on the species separation. The same problem was considered by Mojtabi [29] taking into account the possibility of greatly improving the species separation of a binary mixture in weightlessness by using a rectangular cavity with opposite tall walls moving at equal but opposite optimal velocities. The authors observed that, for a fixed temperature difference, the species separation was optimal for an optimum thickness. The species separation decreased sharply when the thickness decreased. For mixtures with a negative thermo-diffusion coefficient, the heaviest component migrated towards the upper part of the column and the lightest one toward the lower part. The loss of stability of the configuration led to a brutal homogeneity of the binary solution.

Tai and Char [30] analyzed the Soret and Dufour effects on free convective flow of non-Newtonian fluids across a porous medium with thermal radiation. The results indicated that, for aiding flows, the local Nusselt number increased with increasing the power–law index and the Soret number or decreased with the radiation parameter and the Dufour number. Er-Raki et al. [31] analyzed the effect of the Soret effect in a porous cavity. The results showed the significant effect of the Soret parameter on the vertical boundary layer thickness for aiding and opposing flows (N < 0 and N > 0), the boundary layer thickness increased when the Soret parameters increased. Lakshmi et al. [32] investigated the effect of Soret and Dufour diffusion on natural convection in a saturated porous medium. The results indicated that the Nusselt number increased linearly with the increase of Dufour parameter for aiding buoyancy. However, it decreased nonlinearly with the Lewis number. In addition, it increased when the Soret parameter increased, and the Sherwood number increased nonlinearly with increasing Lewis number and decreased linearly with increasing Soret and Dufour parameters. Tsai and Huang [33] discussed numerically the Soret and Dufour effects on natural convection flow over a vertical plate with a power-low heat flux embedded within a porous media. Partha et al. [34] analyzed the effects of Soret and Dufour on thermosolutal convection in a non-Darcy electrically conducting fluid saturated porous medium for both aiding and opposing buoyancy forces. It was observed that the Soret and Dufour effects influenced strongly the heat and mass transfer rates. The magnetic field parameter reduced the heat and mass transfer coefficients.

Teamah [35] performed a numerical study of double-diffusive convective flow in a rectangular enclosure. The vertical walls of the enclosure were subject to constant temperatures and concentration and a uniform horizontal magnetic field, whereas the upper and lower surfaces were insulated and impermeable. The results showed that the heat and mass transfer rates and the flow characteristics within the enclosure depended strongly on the strength of the magnetic field and on the heat generation or absorption effects. In addition, the presence of a heat source or sink slightly reduced the average Sherwood number. Maatki et al. [36] performed a numerical study of double diffusive convection in a cubic enclosure filled with a binary mixture subject to a magnetic field. The results were presented in terms of flow structures, temperature and concentration distributions, and the average Nusselt and Sherwood numbers. The results of this investigation indicated that increasing the intensity of the magnetic field caused a monotonic reduction of intensities of the three-dimensional main transverse flows when the flow was thermally dominated, but there was a significant intensification of three-dimensional flow with a multi-cell structure of secondary flow when the flow was solutably dominated. Costa et al. [37] investigated numerically a natural convection in square porous enclosures under a magnetic field. The Darcy flow model was used. The results were presented and analyzed in terms of streamlines and isotherms. It was shown that the effect of the induced magnetic field resulted in a reduction of the convective flow and the heat transfer rates inside the cavity. Recently, a more complete asymptotic and numerical investigation of double-diffusive convection in the presence of a magnetic field was performed by Rebhi et al. [38,39]. The study was generalized using Neumann and Dirichlet boundary condition types for temperature and solute concentrations, and various convective modes were identified.

In most if not all of the past studies on thermosolutal convection in porous enclosures, based on the above literature survey, the Soret and Dufour effects under the influence of a magnetic field on the unsteady double diffusive natural convection and the determination of the thresholds for the onset of subcritical, oscillatory, and stationary convections have not been addressed yet. Therefore, the present investigation focuses on the examination of the combined effects of Soret and Dufour with the presence of a magnetic field on the unsteady double diffusive natural convection, and on the structure and intensity of the convective flow, and on the heat and mass transfer characteristics inside a horizontal rectangular enclosure saturated with an electrically conducting binary mixture. The numerical confirmation of the stable analytical results is also presented. The critical Rayleigh number characterizing the onset of subcritical convection $R_{TC}^{sub}$ is determined analytically using the parallel flow approximation. The linear stability analysis is also conducted to predict numerically the thresholds for the onset of supercritical, $R_{TC}^{sup}$, overstable, $R_{TC}^{over}$, and oscillatory, $R_{TC}^{osc}$, convection for a wide range of the governing parameters. In addition, the linear stability of the predicted convective state is performed to predict the threshold of Hopf’s bifurcation, $R_{TC}^{Hopf}$, which marks the transition from steady convective flow to oscillatory behavior. Within the present analysis, a whole picture is drawn about various convection modes from rest state to unsteady convection modes within their delineated regions of existence.

2. Geometry and Governing Equation

The physical problem considered in the present work consists of a horizontal saturated porous layer as sketched in Figure 1. The layer has an aspect ratio $A=L\prime /H\prime$, where $H\prime$ is the height and $L\prime$ is the width. Thermo-diffusion and Dufour effects are considered. The porous layer is considered as isotropic and homogeneous, filled with an electrically conducting binary fluid mixture. A magnetic field with uniform strength ${{\overrightarrow{B}}^{\prime }}_0$ is applied in the vertical direction. All the boundaries are assumed to be hydrodynamically impermeable. No external electric field is assumed to exist, and the Hall effect of magneto-hydrodynamics is supposed to be negligible. The magnetic Reynolds number is supposed to be very small so that the induced magnetic field is negligible in front of the applied magnetic field. Both the viscous dissipation and inertial terms are assumed to be negligible. The convective flow in the enclosure is assumed to be laminar, two-dimensional, incompressible and Newtonian. The mixture physical properties are assumed independent of temperature and solute except the density in the buoyancy term, which is supposed to vary according to the Boussinesq approximation. Thus, the density variation with temperature and concentration is described by the linearized state equation, $\rho ={\rho }_0\left[1-{\beta }_T(T^{\prime }-T_0)+{\beta }_S(S^{\prime }-S_{0.})\right]$, where ${\rho }_0$ is the fluid mixture density at temperature $T^{\prime }=T_0$ and mass fraction $S^{\prime }=S_0$, and ${\beta }_T$ and ${\beta }_S$ are the thermal and concentration expansion coefficients, respectively.

According to Gray and Giorgini [40], the Boussinesq approximation was linearly developed with the assumption of small variation in the fluid thermophysical properties over a valid range of temperature, concentration and pressure. The application of the approximation for a given fluid with a reasonable accuracy is valid only within a predetermined range of temperature, concentration and pressure variation.

The equations relating the fluxes of heat, $q^{\prime }$, and matter, ${j^{\prime }}_{ST}$, to the thermal and solute gradients within the binary fluid mixture are given respectively by $q^{\prime }=-k\nabla T^{\prime }-D_{TS}\nabla S^{\prime }$ and ${j^{\prime }}_{ST}^{}=-D\nabla S\prime -D_{ST}\nabla T\prime$, where $k$ and $D$ are the thermal conductivity and the mass diffusivity of the saturated porous medium, respectively. The properties $D_{TS}$ and $D_{ST}$ are respectively the Dufour and the Soret diffusion coefficients.

Assuming a well packed porous medium, so that the Brinkman and form drag effects are negligible, the flow is governed by the Darcy law and the averaged energy and species convection-diffusion equations, which are stated as follows:

$$\frac{\partial u^{\prime }}{\partial x^{\prime }}+\frac{\partial v^{\prime }}{\partial y^{\prime }}=0$$

(1)

$$\frac{\mu }{K}{\overrightarrow{V}}^{\prime }=\nabla p^{\prime }+\rho \overrightarrow{g}+{\overrightarrow{J}}^{\prime }\times {{\overrightarrow{B}}^{\prime }}_0$$

(2)

$$\sigma \frac{\partial T^{\prime }}{\partial t^{\prime }}+{\overrightarrow{V}}^{\prime }\cdot \nabla T^{\prime }=\alpha {\nabla }^2{T}^{\prime }+D_{TS}{\nabla }^2{S}^{\prime }$$

(3)

$$\varphi \frac{\partial S^{\prime }}{\partial t^{\prime }}+{\overrightarrow{V}}^{\prime }\cdot \nabla S^{\prime }=D{\nabla }^2{S}^{\prime }+D_{ST}{\nabla }^2{T}^{\prime }$$

(4)

$${\overrightarrow{J}}^{\prime }=\kappa ({\overrightarrow{E}}^{\prime }+{\overrightarrow{V}}^{\prime }\times {{\overrightarrow{B}}^{\prime }}_0)\begin{array}{cc},& \nabla \cdot {\overrightarrow{J}}^{\prime }=0\end{array}$$

(5)

where ${\overrightarrow{V}}^{\prime }$ is the Darcy velocity vector, $\overrightarrow{g}$ is the gravitational acceleration vector, $\mu$ is the fluid dynamic viscosity, $p^{\prime }$ is the hydrodynamics pressure, $K$ is the porous medium permeability, $\varphi$ is the porosity of the porous medium, ${\overrightarrow{J}}^{\prime }$ is the electric current density, $\kappa$ is the electrical conductivity, ${\overrightarrow{E}}^{\prime }$ is the electrical field magnitude, and $\sigma ={(\rho C)}_p/{(\rho C)}_f$ is the saturated porous medium to fluid heat capacity ratio, with ${(\rho C)}_f$ and ${(\rho C)}_p$ being the heat capacities of the fluid and the saturated porous medium, respectively, and $\alpha$ is the thermal diffusivity of the fluid mixture.

The dimensionless form of the above equations is worked out using the following dimensionless variables:

$\left(x,y\right)=\left(x^{\prime },y^{\prime }\right)/H^{\prime }$, $\left(u,v\right)=\left(u^{\prime },v^{\prime }\right)/U^{*}$, $t=t^{\prime }\alpha /{H^{\prime }}^2$, $P=P^{\prime }/P^{*}$, $\mathsf{\Psi }={\mathsf{\Psi }}^{\prime }/\alpha$, $T=\left(T^{\prime }-T_0\right)/∆T^{*}$ and $S=\left(S^{\prime }-S_0\right)/∆S^{*}$,

where $U^{*}$, $P^{*}$, $∆T^{*}$, and $∆S^{*}$ are the characteristic velocity, pressure, temperature, and solutal scales, defined as:

$U^{*}=\frac{\alpha }{H^{\prime }}$, $P^{*}=\frac{\mathrm{Pr}}{\mathrm{Da}}{\rho }_0{U}^{*2}$, $∆T^{*}=\frac{q^{\prime }{H}^{\prime }}{k}$ and $∆S^{*}=\frac{{j^{\prime }}_{ST}{H}^{\prime }}{D}$,

where $\mathrm{Da}$ is the Darcy number and $\mathrm{Pr}$ is the Prandtl number.

To satisfy the continuity equation, a stream function $\mathsf{\Psi }$ is defined as function of the fluid velocity components such that $u=\partial \mathsf{\Psi }/\partial y$ and $v=-\partial \mathsf{\Psi }/\partial x$.

The dimensionless governing equations for the problem under consideration are obtained as:

$${\nabla }^2\mathsf{\Psi }+{\mathrm{Ha}}^2\frac{{\partial }^2\mathsf{\Psi }}{\partial y^2}=-R_T\frac{\partial }{\partial x}\left(T+\phi S\right)$$

(6)

$$\frac{\partial T}{\partial t}+\frac{\partial \mathsf{\Psi }}{\partial y}\frac{\partial T}{\partial x}-\frac{\partial \mathsf{\Psi }}{\partial x}\frac{\partial T}{\partial y}={\nabla }^2\left(T+D_{u}S\right)$$

(7)

$$\epsilon \frac{\partial S}{\partial t}+\frac{\partial \mathsf{\Psi }}{\partial y}\frac{\partial S}{\partial x}-\frac{\partial \mathsf{\Psi }}{\partial x}\frac{\partial S}{\partial y}={\mathrm{Le}}^{-1}{\nabla }^2\left(S+S_{r}T\right)$$

(8)

The hydrodynamic boundary conditions for the Darcy model are given by:

$$x=\pm \frac{A}{2},u=\mathsf{\Psi }=0\mathrm{and}y=\pm \frac{1}{2},v=\mathsf{\Psi }=0$$

(9)

The thermal and solutal boundary conditions are given by:

$$y=\pm \frac{1}{2},\frac{\partial T}{\partial y}+D_u\frac{\partial S}{\partial y}=-1,\frac{\partial S}{\partial y}+S_r\frac{\partial T}{\partial y}=-1$$

(10)

$$x=\pm \frac{A}{2},\frac{\partial T}{\partial x}=\frac{\partial S}{\partial x}=0$$

(11)

for the impermeable and adiabatic end walls.

From Equations (6)–(8), eight controlling parameters emerged from the dimensionless formation; these are:

$\mathrm{Ha}=B_0\sqrt{\frac{\kappa K}{\mu }}$, $R_T=\frac{\rho g{\beta }_T∆T^{*}{H^{\prime }}^3}{\alpha \mu }$, $\phi =\frac{{\beta }_S}{{\beta }_T}\frac{∆S^{*}}{∆T^{*}}$, $\epsilon =\frac{\varphi }{\sigma }$, $\mathrm{Le}=\frac{\alpha }{D}$, $A=\frac{L^{\prime }}{H^{\prime }}$, $D_u=\frac{D_{ST}∆S^{*}}{\alpha ∆T^{*}}$ and $S_r=\frac{D_{TS}∆T^{*}}{D∆S^{*}}$

which are, respectively, the Hartmann number, Rayleigh number, buoyancy ratio, normalized porosity, Lewis number, aspect ratio, and Dufour and Soret parameters.

The heat and mass transfer rates at any station, $x$, are expressed in terms of the Nusselt, $\mathrm{Nu}$, and Sherwood, $\mathrm{Sh}$, numbers defined as:

$$\mathrm{Nu}{(x)}^{-1}=∆T+D_u∆S=\left(T_{(x,-1/2)}-T_{(x,+1/2)}\right)+D_u\left(S_{(x,-1/2)}-S_{(x,+1/2)}\right)$$

(12)

$$\mathrm{Sh}{(x)}^{-1}=∆S+S_r∆T=\left(S_{(x,-1/2)}-S_{(x,+1/2)}\right)+S_r\left(T_{(x,-1/2)}-T_{(x,+1/2)}\right)$$

(13)

The average Nusselt and Sherwood numbers are computed using the following integral:

$${\mathrm{Nu}}_{\mathrm{m}}=A^{-1}{\displaystyle {\int }_{-A/2}^{A/2}\mathrm{Nu}}(x)\mathrm{d}x$$

(14)

$${\mathrm{Sh}}_{\mathrm{m}}=A^{-1}{\displaystyle {\int }_{-A/2}^{A/2}\mathrm{Sh}}(x)\mathrm{d}x$$

(15)

3. Numerical Solutions

The full governing equations associated with their respective boundary conditions are solved numerically using the finite difference method. The equations are discretized using a second-order central finite difference scheme in time and space on a uniform grid. The discretized governing equations are converted into a system of algebraic equations. The energy and species equations, Equations (7) and (8), are solved iteratively in a time-accurate mode using the alternating-direction implicit method (A.D.I). The system of algebraic equations for each variable is solved by means of the Thomas algorithm. The momentum equation, Equation (6), is solved using the well-known S.O.R. (successive-over-relaxation method). At each new time step, the S.O.R iterative procedure is repeated until the convergence criterion $\left|{\displaystyle {\sum }_{i,j}^m{\mathsf{\Psi }}_{i,j}^{k+1}-{\displaystyle {\sum }_{i,j}^m{\mathsf{\Psi }}_{i,j}^k}}\right|/\left|{\displaystyle {\sum }_{i,j}^m{\mathsf{\Psi }}_{i,j}^{k+1}}\right|\le {10}^{-6}$ is satisfied, where ${\mathsf{\Psi }}_{i,j}^k$ is the stream function value at the node $(i,j)$ at the iteration, $k^{th}$. The results presented in this study are obtained with a grid of 300 × 200, which is judged to be adequate after a grid sensitivity study.

The numerical algorithm is tested for the classical natural convection problem in a rectangular cavity. The grid sensitivity test is performed using a grid size ranging from 50 × 100 to 250 × 300, and the results are listed in

Table 1. Grid sensitivity study for $A=10$, $R_T=100$, $\mathrm{Le}=2$, $\phi =-0.5$, $S_r=D_u=0.1$ and $\mathrm{Ha}=0$.

${\mathit{N}}_{\mathit{x}}\times {\mathit{N}}_{\mathit{y}}$	Numerical Solution					Analytical Solution
	50 × 100	100 × 150	150 × 200	200 × 250	250 × 300
${\mathsf{\Psi }}_0$	3.0327	3.0285	3.0269	3.0262	3.0259	3.0255
Error (%)	1.66	0.10	0.04	0.03	0.02	Reference
$\mathrm{Nu}$	3.1639	3.1496	3.1431	3.1403	3.1389	3.1371
Error (%)	0.85	0.62	0.19	0.10	0.05	Reference
${\mathrm{Nu}}_{\mathrm{m}}$	3.0230	3.0100	3.0042	3.0018	3.0007	…
$\mathrm{Sh}$	4.4160	4.4106	4.4010	4.3954	4.3823	4.3834
Error (%)	0.64	0.52	0.30	0.22	0.11	Reference
${\mathrm{Sh}}_{\mathrm{m}}$	4.2405	4.2346	4.2257	4.2206	4.2179	…

for a rectangular enclosure heated and salted from below having an aspect ratio of $A=10$. The numerical results are computed for $R_T=100$, $\mathrm{Le}=2$, $\phi =-0.5$, $S_r=D_u=0.1$ and $\mathrm{Ha}=0$. Compared to the accurate analytical solution reported in the next section, the results indicated that, for the refined grid, there is a relative difference of 0.02% for the stream function, which is less than 0.05% for the heat transfer rate and nearly 0.11% for the mass transfer rate. Thus, a grid size of 250 × 300 is found to be adequate to simulate accurately the convective flow. The results obtained from the present method are also compared to those reported in the literature as shown in

Table 2. Comparison of ${\mathsf{\Psi }}_0$, $\mathrm{Nu}$, ${\mathrm{Nu}}_{\mathrm{m}}$, $\mathrm{Sh}$ and ${\mathrm{Sh}}_{\mathrm{m}}$ with some previous studies numerical results for $A=10$, $R_T=100$ and Ha = 0.

	$\mathbf{L}\mathbf{e}=10$$,\mathit{\phi }=-0.24$$\mathbf{and}{\mathit{S}}_{\mathit{r}}={\mathit{D}}_{\mathit{u}}=0$			$\mathbf{L}\mathbf{e}=2$$,\mathit{\phi }=1$$,{\mathit{S}}_{\mathit{r}}=-0.1$$\mathbf{and}{\mathit{D}}_{\mathit{u}}=0.1$
	Mamou [18]	Present Study	Present study versus Mamou [18]	Attia et al. [41]	Present Study	Present Study versus Attia et al. [41]
${\mathsf{\Psi }}_0$	3.685	3.652	0.89%	4.6986	4.6987	0.002%
$\mathrm{Nu}$	3.734	3.716	0.48	4.2454	4.2418	0.084
${\mathrm{Nu}}_{\mathrm{m}}$	3.597	3.571	0.72	…	4.1142	…
${\mathrm{Nu}}_{\mathrm{m}}$	6.029	6.031	0.03	5.7574	5.7652	0.1353
${\mathrm{Sh}}_{\mathrm{m}}$	6.726	6.728	0.05	…	5.7314	…

for $A=10$, $R_T=100$ and Ha = 0.

Numerical results are presented in Figure 2a–f for different values of the governing parameters. Streamlines, isotherms, and isoconcentrations are presented in the figure from top to bottom. The results indicate, independently of the governing parameters and for a shallow cavity, $A>>1$, the flow in the central part of the enclosure is parallel, while the temperature and concentration are linearly stratified in the horizontal direction. The analytical solution, developed in Section 4, relies on these numerical observations.

4. Analytical Solution

In this section, an analytical solution is obtained for a thin porous layer, which is equivalent to a cavity with a large aspect ratio $(A>>1)$. In this case, even though the Rayleigh–Bénard multi-cellular flow could be encountered, unicellular cell flow filling up the entire cavity is possible and deemed to be stable even for a high Rayleigh number. In this situation, under steady conditions and a bit away from the end walls, the streamlines remain parallel to the horizontal walls. The flow pattern is assumed to be parallel in the x-direction such that $\mathsf{\Psi }(x,y)\approx \mathsf{\Psi }(y)$. Therefore, the governing Equations (6)–(8) can be considerably simplified under this assumption $\mathsf{\Psi }(x,y)\approx \mathsf{\Psi }(y)$, $T(x,y)\approx C_{T}x+{\theta }_T(y)$, and $S(x,y)\approx C_{S}x+{\theta }_S(y)$, where $C_T$ and $C_S$ are respectively unknown constants of temperature and concentration gradients in the x-direction.

Using these approximations, the governing Equations (6)–(8) are reduced to the following set of ordinary differential equations:

$$\left(1+{\mathrm{Ha}}^2\right)\frac{{\mathrm{d}}^2\mathsf{\Psi }}{\mathrm{d}{y}^2}=-R_T\left(C_T+\phi C_S\right)$$

(16)

$$\frac{{\mathrm{d}}^2{\theta }_T}{\mathrm{d}{y}^2}+D_u\frac{{\mathrm{d}}^2{\theta }_S}{\mathrm{d}{y}^2}=C_T\frac{\mathrm{d}\mathsf{\Psi }}{\mathrm{d}y}$$

(17)

$$\frac{{\mathrm{d}}^2{\theta }_S}{\mathrm{d}{y}^2}+S_r\frac{{\mathrm{d}}^2{\theta }_T}{\mathrm{d}{y}^2}=C_S\mathrm{Le}\frac{\mathrm{d}\mathsf{\Psi }}{\mathrm{d}y}$$

(18)

The solution of Equation (16) satisfying the boundary conditions, Equations (9)–(11), is given by:

$$\mathsf{\Psi }(y)={\mathsf{\Psi }}_0(1-4y^2)$$

(19)

From Equation (19), the velocity field, $u$, derived from the stream function is given by the following expression:

$$u(y)=-8{\mathsf{\Psi }}_{0}y$$

(20)

where ${\mathsf{\Psi }}_0=R_T(C_T+\phi C_S)/[8(1+{\mathrm{Ha}}^2)]$ is the stream function value at the center of the enclosure.

From Equations (17) and (18), the temperature and concentration profiles are obtained as follows:

$$T\left(x,y\right)=C_{T}x+\frac{C_T-C_T\mathrm{Le}{D}_u}{3\left(1-D_u{S}_r\right)}{\mathsf{\Psi }}_0\left(3y-4y^3\right)-a_{{}_T}y$$

(21)

$$S\left(x,y\right)=C_{S}x+\frac{C_S\mathrm{Le}-C_T{S}_r}{3\left(1-D_u{S}_r\right)}{\mathsf{\Psi }}_0\left(3y-4y^3\right)-a_{{}_S}y$$

(22)

where the constants $a_T$ and $a_S$ are defined as:

$$a_T=\frac{1-D_u}{1-D_u{S}_r}\mathrm{and}{a}_S=\frac{1-S_r}{1-D_u{S}_r}$$

(23)

where $C_T$ and $C_S$ are obtained by considering an arbitrary control volume as displayed in Figure 1. The energy and solute balances, in this control volume, are performed to find the constants $C_T$ and $C_S$. The balances’ analysis yields a zero neat heat and solute exchange across any vertical cross section of the porous layer; this is expressed by the following integral:

$${\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\left(\frac{\partial T}{\partial x}+D_u\frac{\partial S}{\partial x}\right)\mathrm{d}y+{\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\frac{\partial \mathsf{\Psi }}{\partial y}T\mathrm{d}y=0}}$$

(24)

$${\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\left(\frac{\partial S}{\partial x}+S_r\frac{\partial T}{\partial x}\right)\mathrm{d}y+\mathrm{Le}{\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\frac{\partial \mathsf{\Psi }}{\partial y}S\mathrm{d}y}}=0$$

(25)

Substituting the temperature, concentration and velocity profiles into Equations (24) and (25) and after performing the integration, it is readily found that the constant gradients of temperature and concentration along the x-direction, $C_T$ and $C_S$, are respectively expressed by:

$$C_T=\frac{8b^2(a_T-a_S{D}_u\mathrm{Le}){\mathsf{\Psi }}_0+4b{\mathrm{Le}}^2{\mathsf{\Psi }}_0^3}{3\left[(2b+{\mathsf{\Psi }}_0^2)(2b+{\mathrm{Le}}^2{\mathsf{\Psi }}_0^2)-D_u{S}_r{(2b-\mathrm{Le}{\mathsf{\Psi }}_0^2)}^2\right]}$$

(26)

$$C_S=\frac{8b^2(a_S\mathrm{Le}-a_T{S}_r){\mathsf{\Psi }}_0+4b\mathrm{Le}{\mathsf{\Psi }}_0^3}{3\left[(2b+{\mathsf{\Psi }}_0^2)(2b+{\mathrm{Le}}^2{\mathsf{\Psi }}_0^2)-D_u{S}_r{(2b-\mathrm{Le}{\mathsf{\Psi }}_0^2)}^2\right]}$$

(27)

Substituting the expressions of $C_T$ and $C_S$, Equations (26) and (27), into the expression of ${\mathsf{\Psi }}_0$, the following fourth-order polynomial equation is obtained:

$${\mathrm{Le}}^2{\mathsf{\Psi }}_0^4+2b_{0}a{\mathsf{\Psi }}_0^2+4b^2-\frac{4b_0^2{R}_T^{0}c}{1+{\mathrm{Ha}}^2}=0$$

(28)

where:

$$\begin{array}{l}a=\mathrm{Le}(\mathrm{Le}+2D_u{S}_r)-\mathrm{Le}(\mathrm{Le}+\phi )R_T/12+1\\ b=b_0(1-D_u{S}_r)\begin{array}{cc},& b_0=15/16\begin{array}{cc},& R_T^0=R_T/12\end{array}\end{array}\\ c=D_u{S}_r(\mathrm{Le}+\phi )-(\mathrm{Le}+1)(D_u\phi +S_r)+\phi \mathrm{Le}+1\end{array}\}$$

(29)

The present analytical solution is assumed to be steady and nonlinear. Therefore, it is possible to examine its existence within given governing parameters’ intervals. The only critical conditions that can be identified by the present analytical solution are the thresholds of subcritical and supercritical convection, when they do exist. Supercritical convection onset occurs at zero flow amplitude, ${\mathsf{\Psi }}_0=0$, thus the corresponding critical Rayleigh number is obtained as:

$$R_{TC}^{sup}=\frac{(1+{\mathrm{Ha}}^2)(1-D_u{S}_r)R^{sup}}{a_T(1-\phi S_r)-a_S\mathrm{Le}(D_u-\phi )}$$

(30)

For an infinite horizontal layer, the constant $R^{sup}$ is computed accurately and is given by $R^{sup}=12$. This result is independent of the type of thermal and solutal boundary conditions imposed on the horizontal walls of the system.

Subcritical convection occurs usually at finite amplitude convection. The critical Rayleigh number, $R_{TC}^{sub}$, marking the bifurcation point, is obtained by deriving ${\mathsf{\Psi }}_0$ in Equation (29) with respect to $R_T$ and then performing the limit $\mathrm{d}{\mathsf{\Psi }}_0/\mathrm{d}{R}_T\to \infty$ at the saddle–node point. After some algebra, the critical Rayleigh number is obtained implicitly as follows:

$$R_{TC}^{sub}=\frac{(1+{\mathrm{Ha}}^2)\left[{\mathrm{Le}}^2{\mathsf{\Psi }}_{0C}^2+b_0\mathrm{Le}(\mathrm{Le}+2D_u{S}_r)+b_0\right]}{b_0\mathrm{Le}(\mathrm{Le}+\phi )}$$

(31)

where ${\mathsf{\Psi }}_{0C}$ is the critical stream function value located at the saddle-node point, which can be computed from:

$${\mathsf{\Psi }}_{0C}=\pm \frac{1}{8{\mathrm{Le}}^2}{\left[-d_1+\sqrt{d_1^2-4{\mathrm{Le}}^2{d}_2}\right]}^{1/2}$$

(32)

with:

$$\begin{array}{c}\hfill {{}_{}^{}}_{}^{}{d}_1=\frac{4b_0\mathrm{Le}}{(\mathrm{Le}+\phi )}\left[1+\phi \mathrm{Le}-(D_u+\phi S_r)(\mathrm{Le}+1)+D_u{S}_r(\mathrm{Le}+\phi )\right]\\ \hfill d_2=\frac{4b_0^2}{\mathrm{Le}(\mathrm{Le}+\phi )}(1+\mathrm{Le}(\mathrm{Le}+2D_u{S}_r)[1+\phi \mathrm{Le}-(D_u+\phi S_r)(\mathrm{Le}+1)+\\ \hfill \begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}{\begin{array}{cc}& \end{array}}_{}^{}{}_{}^{}& & & \end{array}& & & \end{array}& & & \end{array}& & & \end{array}{D}_u{S}_r(\mathrm{Le}+\phi )]-4b^2\end{array}\}$$

(33)

The final two useful engineering parameters are the Nusselt and Sherwood numbers, which express the heat and mass transfer rates through the system, and they are obtained by substituting Equations (21) and (22) into (12) and (13) as follows:

$$\mathrm{Nu}=\frac{3}{3-2{\mathsf{\Psi }}_0{C}_T}\mathrm{and}\mathrm{Sh}=\frac{3}{3-2\mathrm{Le}{\mathsf{\Psi }}_0{C}_S}$$

(34)

The critical Rayleigh number expressions derived in Equations (30) and (31) are similar to those reported by Attia et al. [41] expression for Ha = 0.

5. Stability Analysis

The stability of both motionless and convective states is now considered. The total convective unsteady solution consists of a basic steady-state solution $({\mathsf{\Psi }}_b{,}^{}{T}_b{,}^{}{S}_b)$ and an infinitesimal dynamic perturbation $({\mathsf{\Psi }}_p{,}^{}{T}_p{,}^{}{S}_p)$. The total solution is then expressed by the following equations:

$$\begin{array}{l}\mathsf{\Psi }(x,y,t)={\mathsf{\Psi }}_b(x,y)+{\mathsf{\Psi }}_p(x,y,t)\\ T(x,y,t)=T_b(x,y)+T_p(x,y,t)\\ S(x,y,t)=S_b(x,y)+S_p(x,y,t)\end{array}\}$$

(35)

In general, the basic solution, $({\mathsf{\Psi }}_b{,}^{}{T}_b{,}^{}{S}_b)$, can represent the pure diffusive state solution (${\mathsf{\Psi }}_b$= 0, $T_b=-a_{T}y$, $S_b=-a_{S}y$) or the steady-state convective solution (${\mathsf{\Psi }}_b(x,y)=\mathsf{\Psi }(y)$, $T_b(x,y)=C_{T}x+{\theta }_T(y)$ and $S_b(x,y)=C_{S}x+{\theta }_S(y)$) as predicted by the parallel flow approximation. Thus, the stability of both basic solutions is studied to get the whole picture of the convective system stability.

In an infinite porous layer, the perturbation profiles (${\mathsf{\Psi }}_p$, $T_p$, $S_p$) in space and time can be expressed as follows:

$$\begin{array}{l}{\mathsf{\Psi }}_p(x,y,t)={\psi }_0{e}^{pt+ikx}F(y)\\ T_p(x,y,t)={\theta }_0{e}^{pt+ikx}G(y)\\ S_p(x,y,t)={\varphi }_0{e}^{pt+ikx}H(y)\end{array}\}$$

(36)

where $p=p_r+i{p}_i$ is a complex number, where its real part, $p_r$, expresses the growth rate of the perturbation and the imaginary part, $p_i$, expresses the circular frequency, $k$ is the wave number, and $F(y)$, $G(y)$ and $H(y)$ are one-dimensional space eigenfunctions describing the perturbation profiles ${\mathsf{\Psi }}_p$,

$T_p$, and $S_p$, where ${\psi }_0$, ${\theta }_0$ and ${\varphi }_0$ are unknown infinitesimal amplitudes.

Substituting Equations (35) and (36) into Equations (6)–(8) and neglecting the second-order nonlinear terms, the linearized stability equations are obtained as follows:

$$(1+{\mathrm{Ha}}^2)D^{2}f-k^{2}f=-ik{R}_T(g+h)$$

(37)

$$pg+ikD{\mathsf{\Psi }}_{b}g-ikD{T}_{b}f+C_{T}Df=(D^2-k^2)(g+D_{u}h)$$

(38)

$$\epsilon ph+ikD{\mathsf{\Psi }}_{b}h-ikD{S}_{b}f+C_{S}Df={\mathrm{Le}}^{-1}(D^2-k^2)(h+S_{r}g)$$

(39)

where $D=\mathrm{d}/\mathrm{d}y$, $f={\psi }_{0}F$, $g={\theta }_{0}G$ and $h={\varphi }_{0}H$.

The boundary conditions for the perturbations are given by:

$$y=\pm \frac{1}{2},f=0,Dg=0,Dh=0$$

(40)

The above linear equations system, Equations (37)–(39), subject to boundary conditions in Equation (40), is solved numerically using a finite element method based on the cubic Hermite elements. The discretized linear equations are assembled into a global matrix system, which is obtained as follows:

$$\left[\begin{array}{ccc}\left[{\mathcal{K}}_{\psi }\right]\hfill & \hfill -R_T\left[{\mathcal{B}}_{\psi }\right]& \hfill {}_{}^{}-R_T\phi \left[{\mathcal{B}}_{\psi }\right]\\ \left[{\mathcal{B}}_{\theta }\right]\hfill & \hfill {}_{}^{}\begin{array}{cc}\hfill & \hfill \end{array}\left[{\mathcal{K}}_{\theta }\right]& \hfill {{}_{}^{}}_{}^{}{}^{}{D}_u\left[{\mathcal{L}}_{\theta }\right]\\ \left[{\mathcal{B}}_{\varphi }\right]\hfill & \hfill {{}_{}^{}}_{}^{}\frac{S_r}{\mathrm{Le}}\left[{\mathcal{L}}_{\varphi }\right]& \hfill {{}_{}^{}}_{}^{}{}_{}^{}\frac{1}{\mathrm{Le}}\left[{\mathcal{K}}_{\varphi }\right]\end{array}\right]\left\{\begin{array}{c}f\\ g\\ h\end{array}\right\}=p\left[\begin{array}{ccc}0\hfill & 0& 0\\ 0\hfill & -\left[{\mathcal{M}}_{\theta }\right]& 0\\ 0\hfill & 0& -\epsilon \left[{\mathcal{M}}_{\varphi }\right]\end{array}\right]\left\{\begin{array}{c}f\\ g\\ h\end{array}\right\}$$

(41)

where $\left[{\mathcal{B}}_{\psi }\right]$, $\left[{\mathcal{B}}_{\theta }\right]$, $\left[{\mathcal{B}}_{\varphi }\right]$, $\left[{\mathcal{M}}_{\theta }\right]$, $\left[{\mathcal{M}}_{\varphi }\right]$, $\left[{\mathcal{K}}_{\theta }\right]$, $\left[{\mathcal{K}}_{\varphi }\right]$, $\left[{\mathcal{K}}_{\psi }\right]$, $\left[{\mathcal{L}}_{\theta }\right]$, and $\left[{\mathcal{L}}_{\varphi }\right]$ are square matrixes of dimension $m\times m$, where $m=2N_{ey}+1$ ($N_{ey}$ is the number of elements in the y-direction) is the number of degrees of freedom within the domain $-1/2\le y\le 1/2$. The vectors $\left\{f\right\}$, $\left\{g\right\}$ and $\left\{h\right\}$ are the stream function, temperature and solute eigenvectors of dimension $m$. The associated elementary matrices are obtained as follows:

$$\begin{array}{l}{\left[{\mathcal{K}}_{\psi }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}\left((1+{\mathrm{Ha}}^2)\frac{2}{∆y^{{}^{}}}\frac{\partial {\mathcal{N}}_j}{\partial {\eta }^{{}^{}}}\frac{\partial {\mathcal{N}}_i}{\partial {\eta }^{{}^{}}}+k^2\frac{∆y}{2^{{}^{}}}{\mathcal{N}}_j{\mathcal{N}}_i\right)}\mathrm{d}\eta \begin{array}{cc},& {\left[{\mathcal{M}}_{\theta }\right]}^e={\left[{\mathcal{M}}_{\varphi }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}{\mathcal{N}}_j{\mathcal{N}}_i\frac{∆y}{2}\mathrm{d}\eta ,}\end{array}\\ \begin{array}{cc}& \end{array}{\left[{\mathcal{B}}_{\psi }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}k{\mathcal{N}}_j{\mathcal{N}}_i\frac{∆y}{2}\mathrm{d}\eta }\begin{array}{cc},& {\left[{\mathcal{B}}_{\theta }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}\left(C_T\frac{\partial {\mathcal{N}}_j}{\partial \eta }-ik\left(-a_T+\frac{\partial T_b}{\partial \eta }\right){\mathcal{N}}_j\right){\mathcal{N}}_i\mathrm{d}\eta ,}\end{array}\\ \begin{array}{cc}& \end{array}{\left[{\mathcal{B}}_{\varphi }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}\left(C_S\frac{\partial {\mathcal{N}}_j}{\partial \eta }-ik\left(-a_S+\frac{\partial S_b}{\partial \eta }\right){\mathcal{N}}_j\right){\mathcal{N}}_i\mathrm{d}\eta ,}\\ \begin{array}{cc}& \end{array}{\left[{\mathcal{K}}_{\theta }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}\left(\frac{\partial {\mathcal{N}}_j}{\partial {\eta }^{{}^{}}}\frac{\partial {\mathcal{N}}_i}{\partial {\eta }^{{}^{}}}\frac{2}{∆y}+k^2\frac{∆y}{2^{{}^{}}}{\mathcal{N}}_j{\mathcal{N}}_i+ik\frac{\partial {\mathsf{\Psi }}_b}{\partial {\eta }^{{}^{}}}{\mathcal{N}}_j{\mathcal{N}}_i\right)}\mathrm{d}\eta ,\\ \begin{array}{cc}& \end{array}{\left[{\mathcal{K}}_{\varphi }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}\left(\frac{\partial {\mathcal{N}}_j}{\partial {\eta }^{{}^{}}}\frac{\partial {\mathcal{N}}_i}{\partial {\eta }^{{}^{}}}\frac{2}{∆y^{{}^{}}}+k^2\frac{∆y}{2^{{}^{}}}{\mathcal{N}}_j{\mathcal{N}}_i+ik\mathrm{Le}\frac{\partial {\mathsf{\Psi }}_b}{\partial {\eta }^{{}^{}}}{\mathcal{N}}_j{\mathcal{N}}_i\right)}\mathrm{d}\eta ,\\ \begin{array}{cc}& \end{array}{\left[{\mathcal{L}}_{\theta }\right]}^e={\left[{\mathcal{L}}_{\varphi }\right]}^e={\displaystyle \underset{-1}{\overset{1}{\int }}\left(\frac{\partial {\mathcal{N}}_j}{\partial {\eta }^{{}^{}}}\frac{\partial {\mathcal{N}}_i}{\partial {\eta }^{{}^{}}}\frac{2}{∆y^{{}^{}}}+k^2\frac{∆y}{2^{{}^{}}}{\mathcal{N}}_j{\mathcal{N}}_i\right)}\mathrm{d}\eta \end{array}\}$$

(42)

where ${\mathcal{N}}_i$ is a function of $\eta$ representing the cubic Hermite element test functions. The EigPack double precision routines package for solving generalized eigenvalue problems of the form $[A]\left\{x\right\}=\lambda [B]\left\{x\right\}$ are used to solve the above eigenvalue problem. A similar numerical procedure, using the IMSL library, was used in the past by Mamou and Vasseur [16].

5.1. Stability of the Rest State

The stability response to small perturbations imposed on the quiescent state, ${\mathsf{\Psi }}_b=0$, $T_b=-a_{T}y$ and $S_b=-a_{S}y$ (Equations (6)–(8)) are now considered. The methodology for obtaining the thresholds of various types of convective modes is described hereafter. The eigenvalue problem, Equation (41), is valid for all governing parameter values. To explicitly determine the thresholds of stationary and oscillatory convection, the Galerkin method is the most suitable technique to use provided that the eigenvectors of a given perturbation are determined accurately through the numerical analysis. The eigenvectors $f$, $g$, and $h$ obtained from Equation (41) are used as the weighing functions. For the rest state solution and using F, G and H as weight functions, the Galerkin integration of Equations (37)–(39) leads to the following scalar linear equations:

$${\psi }_0(1+{\mathrm{Ha}}^2){\mathcal{K}}_{\psi }=R_T({\theta }_0+{\varphi }_0\phi )\mathcal{B}$$

(43)

$$p{\theta }_0{\mathcal{M}}_{\theta }+{\psi }_0{a}_T{\mathcal{L}}_{\theta }=-({\theta }_0+{\varphi }_0{D}_u){\mathcal{K}}_{\theta }$$

(44)

$$p\epsilon {\varphi }_0{\mathcal{M}}_{\varphi }+{\psi }_0{a}_S{\mathcal{L}}_{\varphi }=-{\mathrm{Le}}^{-1}({\varphi }_0+{\theta }_0{S}_r){\mathcal{K}}_{\varphi }$$

(45)

where $\mathcal{B}$, ${\mathcal{M}}_{\theta }$, ${\mathcal{M}}_{\varphi }$, ${\mathcal{K}}_{\theta }$, ${\mathcal{K}}_{\varphi }$, ${\mathcal{K}}_{\psi }$, ${\mathcal{L}}_{\theta }$ and ${\mathcal{L}}_{\varphi }$ are constants, which can be computed from the following Galerkin integrals:

$\mathcal{B}={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\frac{\partial G}{\partial y}F\mathrm{d}y}$, ${\mathcal{K}}_{\psi }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{\left(\frac{\partial F}{\partial y}\right)}^2}\mathrm{d}y$, ${\mathcal{K}}_{\theta }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{\left(\frac{\partial G}{\partial y}\right)}^2}\mathrm{d}y$, ${\mathcal{K}}_{\varphi }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{\left(\frac{\partial H}{\partial y}\right)}^2}\mathrm{d}y$, ${\mathcal{L}}_{\theta }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\frac{\partial F}{\partial x}G\mathrm{d}y}$, ${\mathcal{L}}_{\varphi }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}\frac{\partial F}{\partial x}H}\mathrm{d}y$, ${\mathcal{M}}_{\theta }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{G}^2\mathrm{d}y}$, ${\mathcal{M}}_{\varphi }={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{H}^2\mathrm{d}y}$

with: ${\mathcal{K}}_{\theta }={\mathcal{K}}_{\varphi }=\mathcal{K}$, ${\mathcal{L}}_{\theta }={\mathcal{L}}_{\varphi }=\mathcal{L}$ and ${\mathcal{M}}_{\theta }={\mathcal{M}}_{\varphi }=\mathcal{M}$.

Substituting Equations (44) and (45) into Equation (43), we readily arrive at the dispersion relationship, which is stated as follows:

$$(1+{\mathrm{Ha}}^2)\epsilon \mathrm{Le}{p}^2-\gamma p_{1}p-{\gamma }^2{p}_0=0$$

(46)

where

$$\begin{array}{l}{\begin{array}{cc}& \end{array}}_{}^{}{p}_1=R_T^0\mathrm{Le}(a_T\epsilon +a_S\phi )-(1+{\mathrm{Ha}}^2)(\epsilon \mathrm{Le}+1)\\ p_0=R_T^0\left[a_T(1-\phi S_r)-a_S(D_u-\phi )\right]-(1+{\mathrm{Ha}}^2)(1-D_u{S}_r)\\ {\begin{array}{cc}& \end{array}}_{}^{}{}_{}^{}{R}_T^0=\frac{R_T}{R^{sup}},\begin{array}{cc}& R^{sup}\end{array}=\frac{{\mathcal{K}}_{\psi }\mathcal{K}}{\mathcal{B}\mathcal{L}},\begin{array}{cc}& \gamma =\frac{\mathcal{K}}{\mathcal{M}}\end{array}\end{array}\}$$

(47)

Using the numerical analysis, the threshold of stationary convection is obtained when the marginal stability occurs, and this corresponds to a zero growth rate $(p=0)$. The linear global system, Equations (43)–(45), can be rearranged to yield the following eigenvalue problem:

$$\left(E-\lambda I\right)F=0$$

(48)

where $\lambda ={\left[{\mathcal{K}}_{\psi }\right]}^{-1}\left[{\mathcal{B}}_{\psi }\right]{\left[\mathcal{K}\right]}^{-1}\left[\mathcal{B}\right]$ and $\lambda =\frac{1-D_u{S}_r}{\left[a_T(1-\phi S_r)-a_S\mathrm{Le}(D_u-\phi )\right]}\frac{1}{R_{{}_T}}$.

In this way, Equation (48) yields $m$ eigenvalues that are classified in ascending order as ${\lambda }_1\le {\lambda }_2\dots \le {\lambda }_m$, where ${\lambda }_m$ is the maximum eigenvalue, which represents the minimum critical thermal Rayleigh number $R_{TC}^{sup}=1/{\lambda }_m$ for the onset of stationary convection, when $R_{TC}^{sup}$ is positive. The critical value $R_{TC}^{sup}$ is referred to hereafter as the supercritical Rayleigh number. The lowest eigenvalue, ${\lambda }_1$, represents the highest critical thermal Rayleigh number, $R_{TC}^{sup}=1/{\lambda }_1$, when $R_{TC}^{sup}$ is negative. Both numerical and Galerkin analyses yield a threshold for the stationary convection, which can be expressed as:

$$R_{TC}^{sup}=\frac{(1+{\mathrm{Ha}}^2)(1-D_u{S}_r)R^{sup}}{a_T(1-\phi S_r)-a_S\mathrm{Le}(D_u-\phi )}$$

(49)

Furthermore, the threshold for the onset of oscillatory convection is obtained when the real part of the eigenvalue, $p=p_r+i{p}_i$, becomes zero ($p_r=0$) but with a finite imaginary part ($p_i\ne 0$), i.e., $p_1\ne 0$. Numerically, the overstable convection threshold can be obtained for given values of the governing parameters. However, the Galerkin method can lead to an explicit expression of the overstable critical Rayleigh number, $R_{TC}^{over}$, as a function of the governing parameters as:

$$R_{TC}^{over}=\frac{(1+{\mathrm{Ha}}^2)(\epsilon \mathrm{Le}+1)R^{sup}}{\mathrm{Le}(a_T\epsilon +a_S\phi )}$$

(50)

The overstable regime is known to exist up to an upper limit of the Rayleigh number, $R_{TC}^{osc}$, where the oscillatory convective regime vanishes. The critical point where the frequency vanishes is obtained from $p_1^2+4p_0=0$. The expression of $R_{TC}^{osc}$ is then developed as:

$$R_{TC}^{osc}\frac{\left[\epsilon \mathrm{Le}\left(a_S\mathrm{Le}(D_u-\phi )-a_T(1-\phi S_r)\right)+2\mathrm{Le}(a_T\epsilon +a_S\phi )(\epsilon \mathrm{Le}+1)\pm \sqrt{\mathsf{\Lambda }}\right]R^{sup}}{2\mathrm{Le}{(a_T\epsilon +a_S\phi )}^2/(1+{\mathrm{Ha}}^2)}$$

(51)

where $\mathsf{\Lambda }={\epsilon }^2{\mathrm{Le}}^2{\left[a_T(1-\phi S_r)-a_S\mathrm{Le}(D_u-\phi )\right]}^2-4\epsilon {\mathrm{Le}}^2(a_T\epsilon +a_S\phi ){\mathsf{\Lambda }}_0$ and ${\mathsf{\Lambda }}_0=(\epsilon \mathrm{Le}+1)\left[a_T(1-\phi S_r)-a_S\mathrm{Le}(D_u-\phi )\right]+\mathrm{Le}(a_T\epsilon +a_S\phi )(1-D_u{S}_r)$.

The above three thresholds are the essential parameters to determine the nature of the linear stability convection, whether it is monotonic, oscillatory, or stationary.

5.2. Stability of the Convective State: Hopf Bifurcation

A stability analysis of the basic convective steady state predicted by the parallel flow approximation is now investigated. At a higher Rayleigh number, the flow intensity becomes significant, and the transition to turbulent flows is a common occurrence. Before reaching a fully turbulent flow, the transitional flows undergo a sequence of oscillatory flow behavior starting from a well-organized oscillatory flow mode and then evolving toward multi-frequencies’ oscillations, then to quasi-periodic flow and then to chaotic flow. The first appearance of periodic oscillatory flows, known as Hopf bifurcation, occurs at a threshold, $R_{TC}^{Hopf}$. The thresholds can be determined by running a stability analysis of the steady state convective flow. Following the numerical stability analysis described above, the Hopf bifurcation threshold, for a given wave number, is obtained when the real part of the growth rate $p_r$ is nil, or there is a transition from negative to positive values. The optimal value of the thresholds is obtained by optimization with respect to the perturbation wave number, $k$.

Some results of the linear stability analysis for an infinite horizontal porous layer are presented in

Table 3. Effects of the Soret, S_r, Dufour, $D_u$, and the Hartmann, $\mathrm{Ha}$, parameters on the threshold of the Hopf bifurcation in an infinite horizontal layer for $\mathrm{Le}=2$, $\phi =-0.5$ and $\epsilon =0.2$.

$\mathbf{H}\mathbf{a}=0$$,{\mathit{D}}_{\mathit{u}}=0$				$\mathbf{H}\mathbf{a}=0$$,{\mathit{D}}_{\mathit{u}}=0$				${\mathit{D}}_{\mathit{u}}=0$$,{\mathit{S}}_{\mathit{r}}=0$
${\mathit{S}}_{\mathit{r}}$	${\mathit{R}}_{\mathit{T}\mathit{C}}^{\mathit{H}\mathit{o}\mathit{p}\mathit{f}}$	${\mathit{A}}_{\mathit{C}}$	${\mathit{f}}_{\mathit{r}}$	${\mathit{D}}_{\mathit{u}}$	${\mathit{R}}_{\mathit{T}\mathit{C}}^{\mathit{H}\mathit{o}\mathit{p}\mathit{f}}$	${\mathit{A}}_{\mathit{C}}$	${\mathit{f}}_{\mathit{r}}$	$\mathbf{H}\mathbf{a}$	${\mathit{R}}_{\mathit{T}\mathit{C}}^{\mathit{H}\mathit{o}\mathit{p}\mathit{f}}$	${\mathit{A}}_{\mathit{C}}$	${\mathit{f}}_{\mathit{r}}$
1.0	477.18	1.23	20.07	1.0	…	…	…	0.0	301.97	2.57	7.77
0.6	452.64	1.49	15.88	0.6	…	…	…	1.0	414.82	2.25	7.20
0.3	397.58	1.99	11.34	0.3	315.74	3.00	6.51	2.0	712.90	1.90	6.64
0.0	301.97	2.57	7.77	0.0	301.97	2.57	7.77	5.0	2615.07	1.53	5.97
−0.3	217.05	2.88	5.78	–0.3	282.38	2.16	9.12	10.0	9249.76	1.40	5.70
−0.6	153.32	3.08	4.32	–0.6	252.55	1.73	10.80	50.0	220,863.58	1.34	5.55
−1.0	…	…	…	–1.0	208.85	1.40	12.48	100.0	882,124.38	1.35	5.55

. The table illustrates the influence of Soret and Dufour ($S_r$ and $D_u$) effect and the Hartmann parameter, $\mathrm{Ha}$, on the perturbation wavelength, $A_C=2\pi /k$, and oscillation frequency, $f_r=p_i/2\pi$, at the onset of Hopf bifurcation, $R_{TC}^{Hopf}$. The results are obtained for $\epsilon =0.2$. As shown in the table, the decrease of the Soret, $S_r$, and Dufour, $D_u$, parameters below, $S_r=D_u=1$, causes the onset of Hopf bifurcation, $R_{TC}^{Hopf}$, to decrease significantly. It follows that the steady convective flow is destabilized earlier with the decreasing of $S_r$ and $D_u$. The increase of the $\mathrm{Ha}$ parameter has a strong stabilizing effect and induces a reduction in the wavelength and the oscillatory frequency at $D_u=S_r=0$.

6. Results and Discussion

In this paper, analytical and numerical studies have been investigated to examine the combined Soret and Dufour effects, with the presence of a magnetic field, on the double diffusive natural convection, the flow structure, and on the heat and mass transfer rates. The numerical solutions are obtained for a buoyancy ratio $\phi =-0.5$, an aspect ratio $A=10$, and Lewis number $\mathrm{Le}=10$. The other parameter ranges are: $8\le R_T\le {10}^5$, $-1\le S_r\le 1$, $-1\le D_u\le 1$ and $0\le \mathrm{Ha}\le {10}^2$.

Figure 3 displays a comparison between the numerical and the parallel flow solutions for the stream function, $\mathsf{\Psi }$, the horizontal velocity component, $u$, the temperature, $T$, and concentration, $S$, profiles, respectively at the mid-width of the layer for $R_T=500$ and different values of Hartmann number, Ha. The profiles show an anti-symmetric trend with respect to the mid-horizontal plane of the enclosure. The parallel flow prediction presented by solid lines, Equations (19)–(22), is seen to be in an excellent agreement with the numerical solution of the full governing equations, which is depicted by empty circles, thus demonstrating the validity of the parallel flow approximation. Figure 3a shows that, when the Hartman number decreases, the flow intensity increases and therefore strengthens the convective flow inside the cavity. A similar behavior is confirmed by the vertical profile of the horizontal velocity component, $u$, as illustrated in Figure 3b. The effect of the Hartmann number, Ha, on the temperature profile across the porous layer is illustrated in Figure 3c. As anticipated, the temperature difference between the horizontal walls decreases with decreasing Ha and causes a significant increase in the convective heat transfer rate. The magnetic field effect on the concentration profile is depicted in Figure 3d, and it looks similar to the effect of the magnetic field observed on the temperature profile.

For an infinite aspect ratio layer, Figure 4 and Figure 5 show the dependence of the stream function at the center, ${\mathsf{\Psi }}_0$, Nusselt number, $\mathrm{Nu}$, and Sherwood number, $\mathrm{Sh}$, on $R_T$, $S_r$ and $D_u$ at $\mathrm{Ha}=1$. A good agreement is observed between the parallel flow approximation presented by solid lines and the numerical solution displayed by solid symbols. The solid lines correspond to stable branches and the dot-dot-dashed lines to unstable ones, which could not be sustained numerically. It is noticed that the transformation $\mathsf{\Psi }\to -\mathsf{\Psi }$ does not alter the governing equations and the boundary conditions. For this reason, the results are presented only in the first and fourth quadrant, while those in the second and third quadrant being symmetrical are omitted. The results show that, for $S_r=1^{}(D_u=-1)$ and $R_T>530$$(R_T>260)$, the strength of convection becomes quite large, giving rise eventually to oscillating flow. From a physical point of view, the transition to oscillatory or chaotic flows is unavoidable for large values of the Rayleigh number; thus, a further numerical investigation is performed to examine the flow behavior at a high value of Rayleigh number, and the thresholds for transition or Hopf bifurcation is determined. The results clearly indicate that, for given values of the Soret, $S_r$, and the Dufour, $D_u$, parameters, there exist a supercritical or a subcritical Rayleigh numbers ($R_{TC}^{sup}$ or $R_{TC}^{sub}$) for the onset of unicellular convection. In the case of $S_r=1$$(D_u=-1)$, Figure 4a and Figure 5a show that the onset of convection occurs at $R_{TC}^{sup}=16$($R_{TC}^{sup}=8$), below the critical value, and the solution is purely conductive (${\mathsf{\Psi }}_0=0$, $\mathrm{Nu}=\mathrm{Sh}=1$), regardless of the perturbation amplitude magnitude. In the case of $S_r<0.6$ $(D_u>-0.3)$, Figure 4a and Figure 5a show that any decrease of the Soret (Dufour) parameter induces a convection below the threshold of stationary convection, $R_{TC}^{sup}$, which demonstrates the existence of subcritical convection, which is triggered at a subcritical Rayleigh number, $R_{TC}^{sub}$, as a function of the parameters $S_r$ and $D_u$. In the case where $S_r=D_u=0$, the bifurcation curve, predicted by the parallel flow theory, indicates that the onset of motion occurs at a subcritical Rayleigh number of $R_{TC}^{sub}=72$ at which: ${\mathsf{\Psi }}_0=0.976$, $\mathrm{Nu}=1.390$, and $\mathrm{Sh}=2.266$. The strength of convection, ${\mathsf{\Psi }}_0$, increases monotonously with the Rayleigh number, $R_T$, and it becomes more significant when the Soret and Dufour numbers take negative values. However, for positives values, it experiences an inverse trend. From Figure 4b,c and Figure 5b,c, it is clear that, for large values of $R_T$, both $\mathrm{Nu}$ and $\mathrm{Sh}$ tend asymptotically toward a constant value, $\mathrm{Nu}=\mathrm{Sh}\to 6.0$, according to Equation (34) and this is independent of the Soret and Dufour parameters. Usually, when $\mathrm{Nu}$ and $\mathrm{Sh}$ hit the asymptotic values, oscillatory flows prevail and early transition to turbulence occurs.

A more complete view of the effects of $\mathrm{Ha}$, $S_r$, and $D_u$ on ${\mathsf{\Psi }}_0$, $\mathrm{Nu}$, and $\mathrm{Sh}$ is presented in Figure 6 and Figure 7 for $R_T=500$. An excellent agreement between the parallel flow analysis and numerical results is observed, within the range of the governing parameters considered here. In the absence of the magnetic field Ha = 0 and for $D_u=S_r=0$, the resulting curves are the same as those obtained by Bourich et al. [1] while investigating natural convection in a shallow porous cavity under a magnetic field modeled according to the Brinkman–Darcy model. The graphs clearly illustrate the effect of the Hartmann number, $\mathrm{Ha}$, which is having a pre-dominant effect on the strength of convection, ${\mathsf{\Psi }}_0$, and on the heat and mass transfer rates ($\mathrm{Nu}$ and $\mathrm{Sh}$). In general, it is seen from Figure 6 and Figure 7 that, for a given Soret (Dufour) number, an increase in the Hartmann number reduces considerably the strength of the convective motion in the layer. Consequently, as can be observed from Figure 6b,c and Figure 7b,c, both heat and mass transfer rates are considerably inhibited with the increase of the Hartmann number. The plots also show the effect of Soret (Dufour) on the existence of supercritical and subcritical convection, which are similar to those discussed in Figure 4 and Figure 5. As can be seen from Figure 6a and Figure 7a, in the absence of the Soret and Dufour effects $(S_r=D_u=0)$, the resulting bifurcation diagram indicates the existence of a subcritical pitchfork bifurcation. The upper stable branch (solid line) is connected to the lower unstable one (dashed line) at a turning point (saddle-node) occurring at $\mathrm{Ha}={\mathrm{Ha}}_c^{sub}=3.59$. The case of supercritical convection is given in Figure 6a and Figure 7a, as exemplified by the results obtained for $D_u=-1$, for which $\mathrm{Ha}={\mathrm{Ha}}_c^{sup}=7.84$ (Figure 6a) and $\mathrm{Ha}={\mathrm{Ha}}_c^{sup}=11.13$ (Figure 7a). Above these values $(\mathrm{Ha}\ge {\mathrm{Ha}}_c^{sup})$, the rest state is eventually reached, for which heat and mass transfer rates are ruled purely by conduction, ${\mathsf{\Psi }}_0\to 0$ and $\mathrm{Nu}(\mathrm{Sh})\to 1$.

Figure 8 shows the variation of the critical Hartmann numbers, ${\mathrm{Ha}}_c^{sub}$ and ${\mathrm{Ha}}_c^{sup}$, as function of $S_r$ and $D_u$. It is found that, upon decreasing (increasing) the Soret (Dufour) parameters ranging between $S_r=-1$ to 1 ($D_u=1$ to −1), the range of the critical Hartmann number decreases sharply. The graph indicates that, within the range $-1\le S_r\le 0.46$$(-0.24\le D_u\le 1)$, the subcritical Hartmann number is the lowest critical Hartmann number $({\mathrm{Ha}}_c^{sub}\le {\mathrm{Ha}}_c^{sup})$, such that the onset of motion is subcritical, as exemplified in Figure 6 for $S_r=-0.3$ and in Figure 7 for $D_u=0.3$. At $S_r=0.46$$(D_u=-0.24)$, it is found that ${\mathrm{Ha}}_c^{sub}={\mathrm{Ha}}_c^{sup}=5.297$$({\mathrm{Ha}}_c^{sub}={\mathrm{Ha}}_c^{sup}=5.385)$. Upon increasing the value of the Soret and Dufour parameters, it is observed that, in the range $0.46\le S_r\le 1$, the Hartmann threshold is the supercritical one, ${\mathrm{Ha}}_c^{sup}$, as exemplified in Figure 6a at $S_r=0.6$ and in Figure 7a at $D_u=-0.3$.

The influence of the Soret and Dufour parameters on the thresholds of subcritical, $R_{TC}^{sub}$, and supercritical, $R_{TC}^{sup}$, convection is depicted in Figure 9a,b for various values of $\mathrm{Ha}$. The subcritical Rayleigh number is evaluated from the analytical solution, Equation (31), by calculating numerically which value of $R_T$ leads to a zero inverse derivative of ${\mathsf{\Psi }}_0$ with respect to $R_T$. The results, depicted by solid lines in the graphs, correspond to the thresholds of a subcritical unicellular finite amplitude convective regime $({\mathsf{\Psi }}_0\ne 0)$. The dashed lines are the prediction of the linear stability theory, Equation (49). It is noted, from Figure 9a,b, that, upon increasing the value of the Hartmann number, Ha, the thresholds $R_{TC}^{sub}$ and $R_{TC}^{sup}$ increase since the magnetic field effect becomes more and more stabilizing. As shown in Figure 9a,b, for a given value of the Hartmann number, upon decreasing the Soret $(0.4\le S_r\le -1)$ and Dufour $(1\le D_u\le -0.2)$ numbers, it is seen that the supercritical Rayleigh number, $R_{TC}^{sub}$, tends asymptotically towards $R_{TC}^{sup}$ as $S_r\to 0.4$$(D_u\to -0.2)$. This limit $(R_{TC}^{sub}=R_{TC}^{sup})$ is independent of the Hartmann number, Ha. The graph also indicates the existence of a supercritical bifurcation laying within the range $0.4\le S_r\le 1$$(-0.2\le D_u\le -1)$. On the other hand, upon decreasing the Soret and Dufour parameters below $S_r=1$ and $D_u=-1$, the threshold, $R_{TC}^{sup}$, increases sharply toward, $R_{TC}^{sup}\to \infty$.

The stability diagram, as predicted by the linear stability analyses in terms of critical Rayleigh numbers $R_{TC}^{Hopf}$, $R_{TC}^{sub}$, $R_{TC}^{over}$, $R_{TC}^{osc}$, and $R_{TC}^{sup}$ (Equations (31), (49)–(51)), is depicted in Figure 10. As already mentioned, for an increase of Soret parameter above $S_r=-1$, it is seen that the subcritical Rayleigh number decreases monotonically towards $R_{TC}^{sub}\to R_{TC}^{sup}=40$ at $S_r=0.4$. Upon increasing the value of the Soret parameter above $S_r=0.4$, it is observed that the onset of steady motion is supercritical, $R_{TC}^{sup}$. Numerically, it is found that, upon decreasing the value of the Rayleigh number, below these critical Rayleigh numbers ($R_{TC}^{sub}$ and $R_{TC}^{sup}$), the convective flow remains at rest, which corresponds to the stable diffusive regime in which all perturbation decays in time, region (I). The graph also indicates that the critical Rayleigh number for the onset of supercritical, $R_{TC}^{sup}$, convection decreases sharply toward a steady finite amplitude convective regime ($R_{TC}^{sup}\to R_{TC}^{sub}=86.02$), as the Soret number is decreased below zero, $(S_r<0)$. On the other hand, upon increasing the value of the Soret number above zero, $(S_r>0)$, the onset of supercritical convection decreases monotonously toward $R_{TC}^{sup}\to 0$. In general, it is observed that $R_{TC}^{sup}\to \infty$ when $S_r\to 0$. This results from the non-existence of a stationary convection solution. The critical Rayleigh number, $R_{TC}^{Hopf}$, at which a Hopf bifurcation occurs, decreases considerably upon decreasing the Soret number toward $S_r=-0.83$, the value at which the condition $R_{TC}^{Hopf}\to R_{TC}^{sub}=128.64$ is reached. In region (II), the linear theory predicts a stable rest state while a finite amplitude convection is possible according to the nonlinear theory. In region (III), the system is unstable, so any arbitrary dynamic perturbation can initiate a convective flow. For a given value of the normalized porosity equal to 0.35, the linear stability theory predicts the possibility of the existence of oscillating flow, within the overstable regime (zone (IV)) which is delineated by the hatched area (i.e., delineated by $R_{TC}^{over}$ and $R_{TC}^{osc}$, Equations (48) and (49)). In the overstable region, when initiated from the rest state, convection is amplified in an oscillatory way. At $S_r=-0.18$, it is found that $R_{TC}^{sub}=R_{TC}^{over}=R_{TC}^{osc}=R_{TC}^{sup}=84.72$, where a co-dimension-two point occurs. In zone (V), where $R_T\ge R_{TC}^{Hopf}$, the convective flow is oscillatory.

The effect of the thermal Rayleigh number, $R_T$, on the flow intensity and the heat and mass transfer rates for $\mathrm{Ha}=1$, $D_u=S_r=0.1$ and $\epsilon =0.28$, is depicted in Figure 11. The curves presented in the graphs show the present analytical and numerical nonlinear solutions. The numerical solutions of the full governing equations, obtained for $A=20$, are shown by circles (solid circles for stable solution and empty circles with dashed lines for unsteady solution). In the stable regime, a good agreement is observed between the analytical and numerical solutions, and the flow structure consists of a unicellular flow filling up the entire cavity (i.e., the flow in the core region of the layer is parallel) as exemplified by the streamline patterns included in the Figure 11a. In this figure, ${\mathsf{\Psi }}_0$, is the flow intensity at the center of the enclosure for both analytical and numerical stable solutions, with ${\mathsf{\Psi }}_0^{\tau }$ being the numerical averaged flow intensity over a period of time of the oscillation. In the Figure 11b, $\mathrm{Nu}$ and $\mathrm{Sh}$ are the Nusselt and Sherwood numbers at the center of the enclosure and ${\mathrm{Nu}}^{\tau }$ and ${\mathrm{Sh}}^{\tau }$ are the time-and space-averaged values, respectively. For the values of the governing parameters considered here, the thresholds of bifurcations, as predicted by Equations (31) and (49)–(51) and by the linear stability theory discussed in Section 5, are given by $R_{TC}^{sub}=78.47$, $R_{TC}^{over}=93.60$, $R_{TC}^{osc}=142.96$, $R_{TC}^{sup}=174.24$ and $R_{TC}^{Hopf}=475.3$. The numerical results presented in Figure 11a,b indicate that, below the subcritical Rayleigh number ($R_T\le R_{TC}^{sub}$), region (I), the rest state prevails $({\mathsf{\Psi }}_0\to 0{,}_{}^{}\mathrm{Nu}(\mathrm{Sh})\to 1)$. In region (II), delineated by, $R_{TC}^{sub}\le R_T\le R_{TC}^{over}$, the parallel flow theory predicts the existence of two solutions (one stable and one unstable). In region (III), delineated by $R_{TC}^{over}\le R_T\le R_{TC}^{osc}$, the overstable regime where any perturbation grows in an oscillatory manner. The overstable regime extends to an upper limit $R_{TC}^{osc}$, where the oscillation frequency vanishes. Region (IV), where $R_{TC}^{osc}\le R_T\le R_{TC}^{sup}$, represents the stationary convection regime, while region (V), where $R_T\ge R_{TC}^{Hopf}$, corresponds to the oscillatory finite amplitude convection that occurs right above the threshold of Hopf bifurcation. A typical example of the oscillatory periodic flow, near the thresholds of Hopf bifurcation, is presented by comparing the linear stability analysis results to the nonlinear numerical prediction obtained right above the onset of Hopf bifurcation. The onset of oscillatory motion, predicted by the linear stability theory described in Section 5.2, occurs at $R_{TC}^{Hopf}=475.3$ with a corresponding perturbation wavelength of $A_C=2.37$ and frequency $f_r=7.01$. This point is confirmed by the numerical results obtained from the full governing equations developed in Section 3 for a cavity with an aspect ratio of $A=20$. As expected, right below the threshold $R_T<R_{TC}^{Hopf}=475.3$, at $R_T=470$, the numerical solution is steady, and just above this value $(R_T>R_{TC}^{Hopf}=475.3)$ for $R_T=520$, the solution is found to be oscillatory in a periodic manner, indicating an existence of a single oscillation convective mode, as predicted by the linear stability analysis. The results presented in the graph, by empty circle and dashed lines, correspond to the time-averaged values, ${\mathsf{\Psi }}_0^{\tau }$, ${\mathrm{Nu}}^{\tau }$ and ${\mathrm{Sh}}^{\tau }$(see Figure 12d). It is noticed that the increase of the Rayleigh number well above the threshold of the Hopf bifurcation leads to a periodic and then to chaotic oscillatory convective flows. For $R_T=520$, the flow remains unicellular, but the parallel nature of the streamlines is slightly broken, which indicates the existence of small vortices’ layers traveling along the horizontal wall, as depicted in Figure 12d for $R_T=520$. A similar trend has been reported by Bahloul et al. [42] for the case of both double-diffusive and Soret induced convection. The time evolution of ${\mathsf{\Psi }}_0$, $\mathrm{Nu}$, and $\mathrm{Sh}$, obtained for $R_T=520$, is illustrated in Figure 12a–c for τ = 1.0.

Snapshots of the perturbation profiles of the stream function, ${\mathsf{\Psi }}_p$, temperature, $T_p$, and solute, $S_p$, during a period of oscillation, points (1)–(7), are presented in Figure 13a–c. At each given time step, the perturbation field is computed as $F_P=F^{\tau }-1/\tau {\displaystyle \underset{\tau }{\overset{t+\tau }{\int }}Fdt}$, where $F$ stands for $\mathsf{\Psi }$, $T$, and $S$. As illustrated in Figure 13, the convective perturbation patterns are exemplified by two layers consisting of a series of small counter-rotating vortices traveling along the horizontal wall from left to right near the bottom wall and from right to left near the top wall. The two vortex layers are seen to travel in opposite directions, leading to a temporal and partial merging and separation of the vortices. The vortices are seen to become weak as they approach the end walls of the enclosure, and restore progressively their strength later on as they quit the end walls. The patterns of the vortices are quite similar to those predicted by the stability analysis. At the onset of Hopf bifurcation, there exist two symmetrical solutions at the same circular frequency $p_i$ (imaginary part of $p$) but with opposite signs, see Figure 14. The numerical solution indicates that the critical wavelength and the oscillation frequency are given by $A_C=1.25$ and $f_r=7.31$, respectively, which are very close to the values $A_C=2.37$ and $f_r=7.01$ predicted by the linear stability theory. The dynamics of the incipient perturbation flow patterns are illustrated in Figure 14, and once superposed, they appear to be similar to those obtained numerically in Figure 13.

Figure 15 shows the influence of the Hartmann number on the critical Rayleigh numbers for $D_u=S_r=0.1$ and 0.28. In Figure 15a, the different regions are outlined by the thresholds of subcritical, $R_{TC}^{sub}$, overstable, $R_{TC}^{over}$, oscillatory, $R_{TC}^{osc}$, stationary convection, $R_{TC}^{sup}$, and Hopf bifurcation, $R_{TC}^{Hopf}$, in which different flow behaviors may occur. Zones (I) to (V) are equivalent to those defined in Figure 11. As expected, the graph indicates that the onset of Hopf bifurcation, $R_{TC}^{Hopf}$, increases sharply upon decreasing the value of the Hartmann number. It follows that the steady parallel flow is more and more stabilized as Ha increases. A similar trend is observed for the evolution of $R_{TC}^{sub}$, $R_{TC}^{over}$, $R_{TC}^{osc}$, and $R_{TC}^{sup}$ with the Hartmann number variation. On the other hand, for small values of $\mathrm{Ha}$$(\mathrm{Ha}\to 0)$, the thresholds of bifurcations tend toward constant values. The influence of the Hartmann number, Ha, on the perturbation wavelength, $A_C$, and oscillation frequency, $f_r$, at the onset of Hopf bifurcation, $R_{TC}^{Hopf}$, is illustrated in Figure 15b. As can be seen, the magnetic field has a strong effect on the wavelength and oscillatory frequency. It is found that the wavelength and frequency decrease monotonously towards asymptotic values as $\mathrm{Ha}\to \infty$.

7. Conclusions

In this paper, we studied the combined influence of magnetic field, Dufour and Soret effects on the onset of double diffusive convection and on the Hopf bifurcation within an electrically conducting binary mixture, confined inside a horizontal porous cavity. The horizontal walls were subjected to uniform fluxes of heat and mass, whereas the vertical walls were assumed to be adiabatic and impermeable.

An asymptotic analytical solution was developed on the basis of the parallel flow approximation, and a numerical solution was obtained by solving the full governing equations. For shallow enclosure, an excellent agreement was obtained between the analytical and numerical results with the range of the governing parameters considered in the present study.

It was found that the convective flow intensity and the heat and mass transfer rates decreased as the value of the Hartmann number increased, and, for $\mathrm{Ha}\to \infty$, the rest state was eventually reached and the heat and mass transfer rates were essentially driven by conduction, ${\mathsf{\Psi }}_0\to 0$ and $\mathrm{Nu}(\mathrm{Sh})\to 1$. Thus, when the Rayleigh number was very large $(R_T\to \infty )$, both heat and mass transfer rates tended asymptotically toward a constant value $\mathrm{Nu}(\mathrm{Sh})\to 6.0$, independently of the values of Soret, Dufour, and Hartmann parameters.

The Soret, Dufour, and Hartmann numbers had a strong influence on the critical Rayleigh numbers, $R_{TC}^{sup}$, $R_{TC}^{over}$ and $R_{TC}^{osc}$, for the onset of supercritical and oscillatory convection as predicted by the linear stability theory, namely.

The thresholds for subcritical, oscillatory, and stationary convection were obtained explicitly as functions of the governing parameters.

The analysis showed that the Soret, Dufour and Hartmann parameters had a significant influence on heat and mass transfer rates, streamlines, isotherms, isoconcentrations lines and on the critical Rayleigh number, $R_{TC}^{Hopf}$, of Hopf bifurcation.