Finite Amplitude Oscillatory Convection of Binary Mixture Kept in a Porous Medium

Abstract

In the present study, the double-diffusive oscillatory convection of binary mixture, ${}^{3}He$–${}^{4}He$, in porous medium heated from below and cooled from above was investigated with stress-free boundary conditions. The Darcy model was employed in the governing system of perturbed equations. An attempt was made, for the first time, to solve these equations by using the nonlinear analysis-based truncated Fourier series. The influence of the Rayleigh number (R), the separation ratio ($\psi$) due to the Soret effect, the Lewis number ($Le$), and the porosity number ($\chi$) on the field variables were investigated using the finite amplitudes. From the linear stability analysis, expressions for the parameters, namely, R and wavenumbers, were obtained, corresponding to the bifurcations such as pitchfork bifurcation, Hopf bifurcation, Takens–Bogdnanov bifurcation and co-dimension two bifurcation. The results reveal that the local Nusselt number ($N_L$) increases with R. The total energy is enhanced for all increasing values of R. The deformation in the basic cylindrical rolls and the flow rate are enhanced with R. The trajectory of heat flow was studied using the heatlines concept. The influence of R on the flow topology is depicted graphically. It is observed that the intensity of heat transfer and the local entropy generation are increased as R increases.

1. Introduction

The convection of a binary mixture, kept in a porous medium, provides a rich analysis for the study of non-equilibrium phenomena. The antagonistic relationship between the flow-driven by buoyancy effects, namely the density difference caused by the temperature and concentration variations, is often the novel aspect of convection studies. Such phenomena are generally referred to as the natural convection of combined heat and mass transfer, thermohaline or double-diffusive convection. In many engineering applications and problems related to seawater and mantle in the crust of Earth, such double-diffusive convective flow often occurs. Some more examples include the distribution of chemical pollutants in water-saturated soil, grain storage facilities and the migration of moisture through the air contained in insulated fibers. The oscillatory instability is a characteristic feature of double-diffusive convection. In different branches of fluid mechanics and condensed matter physics, the convective systems which exhibit oscillatory convection as a first instability show the interesting physical phenomena. In simple fluids, only the oscillatory instability can occur as secondary bifurcation. Since higher bifurcations are more difficult to describe theoretically, it will be useful to have a real physical system that shows oscillatory instability as the first bifurcation.

The oscillatory instability depends on the stabilizing effect of imposed concentration gradient, such a gradient can also be developed in response to the applied temperature gradient. The first instability with the spatial structure of the convection rolls vary from the stationary bifurcation to the oscillatory bifurcation as the control parameters of the system keeps changing. Specifically, when the superfluid is heated from below in horizontal convective system, it becomes unstable with the oscillatory convection [1]. It is interesting to observe that oscillatory instability occurs along with the concentration perturbations, which stabilizes the system and the temperature perturbations destabilize the convective system. At a certain critical value of temperature gradient, the relaxation time of the temperature fluctuations is much larger than that of the concentration fluctuations.

Further research concerning the normal fluid mixture subjected to the vertical temperature gradient was investigated by Steinberg [2] using the linear stability analysis. He pointed out the existence of oscillatory instabilities when the mixture is heated from below. Later, the influence of two-fluid effects near the onset of convection in ${}^{3}He$–${}^{4}He$ mixture was investigated by Steinberg and Brand [3]. They observed that the potential chemical perturbations reduced the stability of system. Similar flow–visualization experiments in the ethanol-water mixture were performed by Kolodner et al. [4]. They have shown that the conducting state became unstable due to the existence of oscillatory convection. In this convective system, apart from the stability, the occurrence of forward and backward bifurcations over the particular range of were shown by Guenter and Ingo [5] using heat transport measurements. In the ethanol-water and ${}^{3}He$–${}^{4}He$ mixtures, Barbara et al. [6] experimentally presented the data for the critical Rayleigh number and neutral frequency as a function of separation ratio. Based on the sign and magnitude of Soret coefficient, Brand and Steinberg [7] showed that when a layer of two miscible fluids embedded in a porous medium was heated from below or above, any one of two bifurcations will occur, namely, stationary or oscillatory instabilities. Later, Steinberg and Brand [8] illustrated a similar result for a layer of two miscible fluids in a porous medium, that is heated from below or above and with a rapid chemical reaction. Further, the amplitude equation was derived by Brand and Steinberg [9] near the threshold for both stationary and oscillatory instabilities. In a simple fluid, they have shown when heating is conducted from below for the Rayleigh-Bénard convection, the first instability is always stationary. Later, depending on the importance of occurrence of parameters in these systems, Steinberg and Brand [10] predicted the existence of forwarding bifurcation for stationary convection, as well as forwarding and inverse bifurcations for oscillatory convection.

It can be observed from the above studies that the natural convection due to the double diffusion alone has been widely studied and well-reported, while very few studies have been devoted to the double-diffusive convection in porous media. Thermohaline convection is an example of double-diffusive convection. Double-diffusive convection in porous media is the focus of an extensive study due to its significance in predicting the movement of groundwater in aquifers and the method of energy extraction from geothermal reservoirs. Brand and Hohenberg [11] derived an amplitude equation using the weakly nonlinear analysis for the double-diffusive convection and discussed the experimental realizability in alcohol-water and binary-fluid (${}^{3}He$–${}^{4}He$) mixtures. For such systems, Rehberg and Ahlers [12] elaborated on the occurrence of steady or oscillatory bifurcations based on the mean temperature. The thresholds for finite-amplitude, oscillatory and monotonic convective instabilities are calculated both analytically and numerically by Bahloul et al. [13] in terms of the governing parameters of these systems. Using a Galerkin method, Bahloul et al. [14] investigated the nonlinear state of ideal straight rolls arising in these systems. They observed that the stability regions of these rolls are restricted by the cross roll, zigzag and oscillatory instabilities.

For a horizontal layer saturated in a porous medium, the oscillatory instability may be possible when a strongly stabilizing solute gradient is opposed by a destabilizing thermal gradient (Nield [15]). An order of magnitude predictions for the total heat and mass transfer rates in such systems and their valid domains are shown to be in agreement with the results produced by discrete numerical experiments (Trevisan and Bejan [16,17]). Near the onset, when the negative separation ratio ($\psi$) decreases, the traveling wave bifurcation branch is shifted towards the lower values of R without qualitative changes as long as the Lewis number ($Le$) is small and the normalized porosity is large (Augustin et al. [18]). The heat and concentration flux values were obtained experimentally by Griffiths [19] through an interface between two-fluid layers and with different temperature values and salt concentrations. He showed that the differences in horizontal properties in each convecting layer lead to an increased heat flow through the two-layer device due to the bending of interface. Hounsou et al. [20] studied the criterion of the appearance of stationary convection in a porous horizontal layer saturated by a binary mixture of ferro-fluids, which is heated from below. The effect of Darcy number, a ratio of viscosities, the magnetic and binary parameters are illustrated by them on the appearance and size of the convection cells. Recently Hu and Zhang [21] discovered a new oscillatory instability with a positive separation ratio. This ratio showed the interactions among several factors and provided the analytical criteria to predict the onset of Rayleigh -Bénard convection under the ideal boundary conditions. Rashmi and Murthy [22] have established the nonlinear stability theory for a uniform flow along a Brinkman porous layer with two impermeable and isothermal horizontal boundaries. The cross-diffusion effect is also considered as a factor contributing to the convective instability in the medium.

The nonlinear theory of Rayleigh-Bénard convection was carried out by Kuo [23] using the Fourier analysis of perturbations which are valid for all the Rayleigh numbers $R\ge R_{cs}$ (critical stationary convection Rayleigh number) and he compared his theoretical results with the experimental observations. For the convection in porous media, the finite amplitude convection was studied by Palm et al. [24] and their theoretical results showed very good agreement with the experiments. Rameshwar et al. [25,26] and Rawoof Syeed et al. [27], have investigated the steady-state finite amplitude cellular convection for the double diffusive physical models. It can be noted that the research carried by Rameshwar et al. [26] is based on stationary convection with the eigenvalue, $p=0$ and the present research work is based on oscillatory convection with eigenvalue, $p=i\omega$. The dynamic behaviors of stationary convection and oscillatory convection are different in nature. For stationary convection (steady) model, the approximate solutions until $O\left({ϵ}^8\right)$ ($ϵ$ is the expansion parameter) are computed, but in the present oscillatory convection (unsteady) model, due to the high complexity of analytical expressions of the eigenfunctions, we are able compute the approximate solutions until $O\left({ϵ}^5\right)$.

Much less attention is paid to the nonlinear studies of the considered complex physical model, while near the onset of oscillatory convection, only weakly nonlinear studies exist in the literature. Thus, the objectives of the present study were as follows:

• Investigate theoretically the linear and nonlinear behaviour of the Rayleigh-Bénard convection system of a superfluid mixture, ${}^{3}He$–${}^{4}He$, kept in a porous medium;

• Obtain the critical values of control parameters such as R, from the linear stability analysis for the onset of stationary and oscillatory convection;

• Find the R values for different bifurcations viz. pitchfork bifurcation, Hopf bifurcation, Takens–Bogdanov bifurcation and co-dimension two bifurcation;

• Solve the nonlinear partial differential equations using the perturbation method proposed by Kuo [23], until the $O\left({ϵ}^5\right)$ and obtain the approximate solutions (eigenfunctions) to analyse the time-dependent nonlinear behavior of the convective system;

• Find the local ($N_L$) and time–dependent averaged ($Nu$) Nusselt numbers on the hot wall to understand the development of heat flow and rate of heat transfer, respectively;

• Obtain the cellular pattern of the fluid flow, hot regions (isotherms), concentration regions (isoconcentrations) from the eigenfunctions related to streamfunction, temperature and concentration, respectively;

• Study the heatline patterns of the flow by using the heatfunction;

• Study the effect of physical parameters on the entropy generation.

The manuscript is arranged, further, into nine sections. Section 2 lists the hydrodynamic governing equations related to the double-diffusive flow in a porous medium under the Boussinesq approximation. For this mathematical model, in Section 3, the linear stability analysis is presented. In Section 4, the method of solution for this model by using the perturbation method is elaborated. The approximate solutions for the field variables are also presented in this section. In Section 5, the approximate analytical expression for $Nu$ is given. In Section 6, the streamlines and isotherm patterns of the flow are presented. Section 7 and Section 8 deal with the heatline patterns of the flow by using the heat function and the entropy generation, respectively. Finally, the conclusions of the present study are given in Section 9.

2. Mathematical Formulation

The binary mixture ${}^{3}He$–${}^{4}He$, placed in a sparsely packed porous medium of thickness d was considered. This layer is parallel to the horizontal $X^{{}^{\prime }}{Y}^{{}^{\prime }}$-plane with a very large horizontal extension that is kept in the gravitational field $\overrightarrow{g}=-g{\widehat{e}}_z$. The layer has its interfaces at the vertical coordinates $Z^{{}^{\prime }}=0$ and $Z^{{}^{\prime }}=d$. A static temperature difference across the layer was assumed to be imposed. The temperature on the bottom plate $(Z^{{}^{\prime }}=0)$ was assumed to be $T^{{}^{\prime }}=T_0^{{}^{\prime }}+▵T^{{}^{\prime }}$ and on the top plate $(Z^{{}^{\prime }}=d)$ the temperature is $T^{{}^{\prime }}=T_0^{{}^{\prime }}$. The basic equations were taken according to the Boussinesq approximation and the Darcy model. Based on the above assumptions, the momentum equation is given by [7]:

(1)$\begin{array}{c}\hfill \frac{1}{\epsilon }\frac{\partial {\overrightarrow{V}}^{\prime }}{\partial t^{{}^{\prime }}}+\mathbf{a}{\overrightarrow{V}}^{\prime }\left|{\overrightarrow{V}}^{\prime }\right|=-\frac{{\nabla }^{{}^{\prime }}{p}^{{}^{\prime }}}{{\rho }_0}+\frac{\overrightarrow{g}\rho }{{\rho }_0}-\frac{\nu }{K}{\overrightarrow{V}}^{\prime }.\end{array}$

In above equation, the quantity ‘a’ measures the microscale inertial effects such as thermal and particle diffusion anisotropies. The inertial term ${\overrightarrow{V}}^{\prime }\left|{\overrightarrow{V}}^{\prime }\right|$ is associated with ‘a’. In the present study, the term ‘a’ is neglected since the contribution of anisotropy is assumed to be very small [17,28]. The other notations are explained in the Nomenclature. The other hydrodynamic basic equations for the conservation of mass, the concentration, and the temperature are given by [7,29]:

(2)$\begin{array}{c}\hfill {\nabla }^{{}^{\prime }}.{\overrightarrow{V}}^{\prime }=0,\end{array}$

(3)$\begin{array}{c}\hfill \frac{\partial C^{{}^{\prime }}}{\partial t^{{}^{\prime }}}+({\overrightarrow{V}}^{\prime }.{\nabla }^{{}^{\prime }})C^{{}^{\prime }}=D{\nabla }^{{}^{\prime }2}{C}^{{}^{\prime }}-\frac{{\kappa }_T}{T_0^{{}^{\prime }}}D{\nabla }^{{}^{\prime }2}{T}^{{}^{\prime }},\end{array}$

(4)$\begin{array}{c}\hfill \frac{\rho C_p^{{}^{\prime }}}{{\left(\rho C_p^{{}^{\prime }}\right)}_{liq}}\frac{\partial T^{{}^{\prime }}}{\partial t^{{}^{\prime }}}+({\overrightarrow{V}}^{\prime }.{\nabla }^{{}^{\prime }})T^{{}^{\prime }}=\kappa {\nabla }^{{}^{\prime }2}{T}^{{}^{\prime }},\end{array}$

where ${\kappa }_T$ is the Soret coefficient, $\rho C_p^{\prime }/{\left(\rho C_p^{\prime }\right)}_{liq}\simeq 1$ for liquids; additionally, $\rho ={\rho }_0[1-{\beta }_1(T-T_0)-{\beta }_2(C-C_0)]$. Near the onset of convection, Equations (1)–(4) give rise to the static solution as

(5)$\begin{array}{c}\hfill \overrightarrow{V_s^{\prime }}=\overrightarrow{0},T_s^{{}^{\prime }}=T_0-\Lambda Z^{\prime },C_s^{{}^{\prime }}=C_0-\left(\frac{{\kappa }_T}{T_0}\right)\Lambda Z^{\prime },\end{array}$

where $\Lambda$ is the applied temperature gradient. The Equations (1)–(4) are non-dimensionalized using the scales listed below,

$\begin{array}{c}\hfill X=\frac{X^{{}^{\prime }}}{d},Y=\frac{Y^{{}^{\prime }}}{d},Z=\frac{Z^{{}^{\prime }}}{d},t=\frac{\kappa }{d^2}{t}^{{}^{\prime }},u=\frac{d}{\kappa }{u}^{{}^{\prime }},\end{array}$

(6)$\begin{array}{c}\hfill v=\frac{d}{\kappa }{v}^{{}^{\prime }},w=\frac{d}{\kappa }{w}^{{}^{\prime }},T=\frac{{\beta }_{1}gKd}{\nu \kappa }{T}^{{}^{\prime }},C=\frac{T_0{\beta }_{1}gKd}{\nu \kappa {\kappa }_T}{C}^{{}^{\prime }}.\end{array}$

The perturbed velocity, temperature and concentration quantities are given by

(7)$\begin{array}{c}\hfill \overrightarrow{V^{\prime }}=\overrightarrow{V_s^{\prime }}+\overrightarrow{V^{\prime }}*,T^{\prime }=T_s^{\prime }+{\theta }^{\prime }*,C^{\prime }=C_s^{\prime }+C^{\prime }*.\end{array}$

Using Equations (5)–(7), the basic Equations (1)–(4) are transformed into the perturbed dimensionless equations. For convenience, the asterisk is removed from $\overrightarrow{V^{*}}$, ${\theta }^{*}$ and $C^{*}$. The perturbed dimensionless equations are given by

(8)$\begin{array}{c}\hfill (1+\chi \frac{\partial }{\partial t}){▿}^{2}w-{▿}_h^2\theta -\psi {▿}_h^{2}C=0,\end{array}$

(9)$\begin{array}{c}\hfill (\frac{\partial }{\partial t}+\overrightarrow{V}.▿-{▿}^2)\theta -Rw=0,\end{array}$

(10)$\begin{array}{c}\hfill (\frac{\partial }{\partial t}+\overrightarrow{V}.▿-Le{▿}^2)C-Rw+Le{▿}^2\theta =0,\end{array}$

where

${▿}_h^2=\frac{{\partial }^2}{\partial X^2}+\frac{{\partial }^2}{\partial Y^2},{▿}^2=\frac{{\partial }^2}{\partial X^2}+\frac{{\partial }^2}{\partial Y^2}+\frac{{\partial }^2}{\partial Z^2}.$

The control physical parameters in Equations (8)–(10) are given by R = ${\beta }_{1}gdK\Lambda /\kappa \nu$ which accounts for buoyancy effects, $\chi =K\kappa /\nu \epsilon$d relates to the porous medium, $\psi =-\left({\beta }_2{\kappa }_T\right)/\left({\beta }_1{T}_0\right)$, is the separation ratio due to the Soret effect and $Le=D/\kappa$ is the ratio of mass diffusivity to thermal diffusivity. As $Le$ decreases, the thermal diffusivity increases and accordingly the mass diffusivity decreases. In the linear part of Equations (8)–(10), $\theta$ and C are eliminated to obtain the following form:

(11)$\begin{array}{c}\hfill \mathcal{L}w=\mathcal{N},\end{array}$

$\begin{array}{ccc}\hfill \mathcal{L}& =& {\nabla }^2(\frac{\partial }{\partial t}-{\nabla }^2)(\frac{\partial }{\partial t}-Le{\nabla }^2)(1+\chi \frac{\partial }{\partial t})\hfill \\ & & +R{{\nabla }_h}^2[\psi Le{\nabla }^2-\psi (\frac{\partial }{\partial t}-{\nabla }^2)-(\frac{\partial }{\partial t}-Le{\nabla }^2)],\hfill \end{array}$

and

$\begin{array}{ccc}\hfill \mathcal{N}& =& [\psi Le{\nabla }^2-(\frac{\partial }{\partial t}-Le{\nabla }^2)]{\nabla }_h^2(\overrightarrow{V}.\nabla )\theta \hfill \\ & & -\psi (\frac{\partial }{\partial t}-{\nabla }^2){\nabla }_h^2(\overrightarrow{V}.\nabla )C.\hfill \end{array}$

The horizontal planes $Z=0$ and $Z=1$ are assumed to be kept at constant boundary conditions and hence

(12)$\begin{array}{c}\hfill \theta =C=0\mathrm{at}Z=0\mathrm{and}Z=1,\mathrm{for}\mathrm{all}X,Y,\end{array}$

and also the normal component of the velocity should vanish on the horizontal planes, i.e.,

(13)$\begin{array}{c}\hfill w=0\mathrm{at}Z=0\mathrm{and}Z=1,\mathrm{for}\mathrm{all}X,Y.\end{array}$

The conditions (12) and (13) are applicable for any combination of top–bottom boundaries such as free–free or rigid–rigid or rigid–free or free–rigid. The stress-free boundary conditions [30] are considered and are given by

(14)$\begin{array}{c}\hfill \frac{{\partial }^{2}w}{\partial Z^2}=0\mathrm{at}Z=0\mathrm{and}Z=1,\mathrm{for}\mathrm{all}X,Y.\end{array}$

Since the selected physical system is a double-diffusive system, it is unstable to either stationary convection or oscillatory convection at the onset.

3. Linear Stability Analysis

The initial perturbations of the system are quite small during the onset of convection. Consequently, the nonlinear terms are much smaller than the linear terms. Thus, by neglecting the nonlinear contributions, Equation (11) leads to a linear differential equation as $\mathcal{L}w=0.$ This procedure is known as linearization and is described by

(15)$\begin{array}{ccc}\hfill & & \{{\nabla }^2(\frac{\partial }{\partial t}-{\nabla }^2)(\frac{\partial }{\partial t}-Le{\nabla }^2)(1+\chi \frac{\partial }{\partial t}).{\nabla }^2\hfill \\ & & +R{{\nabla }_h}^2[\psi Le{\nabla }^2-\psi (\frac{\partial }{\partial t}-{\nabla }^2)-(\frac{\partial }{\partial t}-Le{\nabla }^2)]\}w=0.\hfill \end{array}$

The linearization around the basic state of momentum, heat and concentration Equations (8)–(10) led to Equation (15). The study of linear stability analysis is carried out by the normal mode procedure. Substituting $w(X,Z,t)=W\left(Z\right)e^{(iaX+pt)}$ in the Equation (15), we get

(16)$\begin{array}{c}\hfill R=\frac{{\mathit{d}}_2\left({\mathit{d}}_2+p\right)\left(Le{\mathit{d}}_2+p\right)\left(1+\chi p\right)}{a^2\left(\psi Le{d}_2+\psi \left(d_2+p\right)+Le{d}_2+p\right)},\end{array}$

where $d_2=a^2+{\pi }^2$, the quantity a is the horizontal wavenumber of the basic cell and p stands for the growth rate of perturbation. If the temperature gradient is very small, conduction will arise and when the gradient crosses the critical value, the convection rolls will start to appear. The basic state of the flow is defined by the absence of motion.

3.1. Stationary Convection $(p=0)$

The considered flow is unstable to stationary convection or oscillatory convection. The classification of the instabilities, such as stationary instability and oscillatory instability, can be studied using the eigenvalue value p. If $p=0$ stationary convection is obtained and oscillatory convection is obtained if p = i$\omega$. The Rayleigh number for stationary convection $(R=R_s)$ is given by the condition $(p=0)$, i.e.,

$\begin{array}{c}\hfill R_s=\frac{Le{d_2}^2}{\left(\psi Le+\psi +Le\right)a^2}.\end{array}$

The critical value of $R_s$ is obtained from $\partial R_s/\partial a=0.$ The critical wavenumber is given by $a^2$ = ${a^2}_{cs}={\pi }^2$ and the critical Rayleigh number for stationary convection ($R_{cs}$) is (Brand and Steinberg [7]):

(17)$\begin{array}{c}\hfill R_{cs}=\frac{4{\pi }^2}{(\psi +\frac{\psi }{Le}+1)}.\end{array}$

The minimum value of $R_{cs}$ for $\psi =0$ and $Le=1$ is $R_{cs}=4{\pi }^2$ with $a=\pi$. The $R_s$ is obtained for a vanishing eigenvalue p, which gives the pitchfork bifurcation.

3.2. Oscillatory Convection $(p=i\omega )$

In Equation (16), it can be noted that R is expressed in a complex form. Since the values of R are always real, the imaginary part of R must be zero. Thus by equating the imaginary part of R to zero we get

(18)$\begin{array}{c}\hfill {\omega }^2=\frac{-{d_2}^2\left(\left({Le}^2\chi \psi +{Le}^2\chi +Le\chi \psi \right)d_2+\left({Le}^2+Le+1\right)\psi +{Le}^2\right)}{{(d_2\chi +\psi +1)}^2}.\end{array}$

From Equations (16) and (18), the R value for oscillatory convection, $R=R_o$, is obtained as (Brand and Steinberg, [7])

(19)$\begin{array}{c}\hfill R_o=\frac{{d_2}^2\left(d_2\chi +1\right)\left(1+Le\right)\left(Le\chi d_2+1\right)}{\left(d_2\chi +\psi +1\right)a^2}.\end{array}$

The $R_0$ can be obtained from the relationship $\partial R_o/\partial a=0$ and the purely imaginary eigenvalue, $p=i\omega$, gives the Hopf bifurcation. The Takens–Bogdanov bifurcation point is obtained when $R_s\left(a_s\right)=R_0\left(a_0\right),a_s=a_0$. The co-dimension two bifurcation is obtained when $R_{cs}\left(a_{cs}\right)=R_{co}\left(a_{co}\right),a_{cs}\ne a_{os}$.

4. Method of Solution

The occurrence of oscillatory convection depends on the physical parameters of the system. The double-diffusive fluid is confined between two horizontal boundaries, in which the bottom and top boundaries are uniformly heated and cooled respectively, the cellular convection begins at critical value $R_{co}$. Initially, the convection has a definite form and scales [7] which are given by solutions of the linearized equations with boundary conditions [31]. As the value of R crosses $R_{co}$, the flow becomes unstable and changes from laminar to turbulent. For a wide range of R the flow remains stable and laminar, followed by an unstable turbulent flow at a much higher temperature difference between the lower and upper plates. The $ϵ$ expansion parameter is chosen as [23]:

(20)$\begin{array}{c}\hfill {ϵ}^2=\frac{R-R_{co}}{R}.\end{array}$

Note that $ϵ$ is less than one for all values of $\mathit{R}$. As in the case of regular perturbation studies, the solution of Equations (8)–(10) are expanded in powers of $ϵ$ as

(21)$\begin{array}{c}\hfill f=ϵf_1+{ϵ}^2{f}_2+{ϵ}^3{f}_3+{ϵ}^4{f}_4+{ϵ}^5{f}_5+{ϵ}^6{f}_6+\cdots ,\end{array}$

where $f=f(u,w,C,\theta ).$ The R, based on Equation (20) is defined by

(22)$\begin{array}{c}\hfill R=\frac{R_{co}}{1-{ϵ}^2},\end{array}$

Expanding Equation (22) in the power series of $ϵ$ or by applying the finite formula, the following expansion is obtained

(23)$\begin{array}{c}\hfill R=R_{co}+R_{os}({ϵ}^2+{ϵ}^4+{ϵ}^6+....+{ϵ}^{2s}),\end{array}$

where

(24)$\begin{array}{c}\hfill R_{os}=\frac{R_{co}}{1-{ϵ}^{2s}},s=1,2,3\dots \end{array}$

By introducing Equations (21) and (23) in Equation (11), and equating to zero the coefficients of $ϵ,{ϵ}^2,{ϵ}^3,\cdots$, a sequence of non-homogeneous partial differential equations are obtained and are given in the following form:

(25)$\begin{array}{c}\hfill {\mathcal{L}}_1{w}_1+R_{co}{\mathcal{L}}_2{w}_1=0,\end{array}$

(26)$\begin{array}{c}\hfill {\mathcal{L}}_1{w}_2+R_{co}{\mathcal{L}}_2{w}_2={\mathcal{N}}_1.\end{array}$

In general,

(27)$\begin{array}{c}\hfill {\mathcal{L}}_1{w}_i+R_{co}{\mathcal{L}}_2{w}_i+R_{os}{\mathcal{L}}_2(w_{i-2}+w_{i-4}+w_{i-6}......)={\mathcal{N}}_{i-1},\mathrm{for}i\ge 2.\end{array}$

where

(28)$\begin{array}{c}\hfill {\mathcal{L}}_1={\nabla }^2\left(\frac{\partial }{\partial t}-{\nabla }^2\right)\left(\frac{\partial }{\partial t}-Le{\nabla }^2\right)\left(1+\chi \frac{\partial }{\partial t}\right),\end{array}$

and

(29)$\begin{array}{c}\hfill {\mathcal{L}}_2={{\nabla }_h}^2\left[\psi Le{\nabla }^2-\psi \left(\frac{\partial }{\partial t}-{\nabla }^2\right)-\left(\frac{\partial }{\partial t}-Le{\nabla }^2\right)\right].\end{array}$

The values of ${\theta }_i$ and $C_i$ are calculated from $w_i$, ($i=$1, 2, 3 … n) using the spectral Equations (9) and (10). The spectral temperature equations are given by

(30)$\begin{array}{c}\hfill \left(\frac{\partial }{\partial t}-{\nabla }^2\right){\theta }_1=R_{co}{w}_1,\end{array}$

(31)$\begin{array}{c}\hfill \left(\frac{\partial }{\partial t}-{\nabla }^2\right){\theta }_2+({\overrightarrow{V}}_1.\nabla ){\theta }_1=R_{co}{w}_2.\end{array}$

In general,

(32)$\begin{array}{c}\hfill \left(\frac{\partial }{\partial t}-{\nabla }^2\right){\theta }_i+\sum _{l=1}^{i-1}({\overrightarrow{V}}_{\mathit{l}}.\nabla ){\theta }_{i-l}=R_{co}{w}_i+R_{os}(w_{i-2}+w_{i-4}+w_{i-6}...),\mathrm{for}i\ge 2.\end{array}$

Similarly, the concentration is given by the auxiliary equations

(33)$\begin{array}{c}\hfill \left(\frac{\partial }{\partial t}-Le{\nabla }^2\right)C_1=R_{co}{w}_1-Le{\nabla }^2{\theta }_1,\end{array}$

(34)$\begin{array}{c}\hfill \left(\frac{\partial }{\partial t}-Le{\nabla }^2\right)C_2+({\overrightarrow{V}}_1.\nabla )C_1=R_{co}{w}_2-Le{\nabla }^2{\theta }_2.\end{array}$

In general,

(35)$\begin{array}{c}\hfill \left(\frac{\partial }{\partial t}-Le{\nabla }^2\right)C_i+\sum _{l=1}^{i-1}({\overrightarrow{V}}_{\mathit{l}}.\nabla )C_{i-l}=R_{co}{w}_i+R_{os}(w_{i-2}+w_{i-4}+w_{i-6}...)-Le{\nabla }^2{\theta }_i,\mathrm{for}i\ge 2.\end{array}$

4.1. Approximate Solutions

The approximate solutions of $w,\theta$ and C are calculated for the oscillatory convection. When both of the boundaries are stress-free, all the space functions of $w,\theta$ and C modes are considered to be the sine and cosine functions. Thus, from Equations (25), (30) and (33), the following first-order approximate solutions are considered:

(36)$\begin{array}{ccc}\hfill & & w_1=A_{1}sin\omega tcosaXsin\pi Z,\hfill \\ & & {\theta }_1=\frac{R_{co}{A}_1\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z}{{d_2}^2+{\omega }^2},\hfill \\ & & C_1=\frac{R_{co}{A}_1\left({Le}^2{d_2}^2\omega +Le{d_2}^2\omega +{d_2}^2\omega +{\omega }^3\right)sin\omega tcosaXsin\pi Z}{\left({Le}^2{d_2}^2+{\omega }^2\right)\left({d_2}^2+{\omega }^2\right)}\hfill \\ & & +\frac{R_{co}{A}_1{d_2}^3\left({Le}^2+Le\right)cos\omega tcosaXsin\pi Z}{\left({Le}^2{d_2}^2+{\omega }^2\right)\left({d_2}^2+{\omega }^2\right)}.\hfill \end{array}$

Generally, the terms given in Equation (21) are written as

(37)$\begin{array}{ccc}\hfill & & w_j=A_{j}sin\omega tcosaXsin\pi Z\hfill \\ & & +\sum _{p=0,q=0,r=0}^{\infty }{w}_{pqr}^{\left(j\right)}\left(\omega sinr\omega t+d_{2}cosr\omega t\right)cospaXsinq\pi Z,\hfill \end{array}$

(38)$\begin{array}{ccc}\hfill & & {\theta }_j=\frac{R_{co}{A}_j\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z}{{d_2}^2+{\omega }^2}\hfill \\ & & +\sum _{p=0,q=0,r=0}^{\infty }{\theta }_{pqr}^{\left(j\right)}\left(\omega sinr\omega t+d_{2}cosr\omega t\right)cospaXsinq\pi Z,\hfill \end{array}$

(39)$\begin{array}{ccc}\hfill & & C_j=\frac{R_{co}{A}_j\left({Le}^2{d_2}^2\omega +Le{d_2}^2\omega +{d_2}^2\omega +{\omega }^3\right)sin\omega tcosaXsin\pi Z}{\left({Le}^2{d_2}^2+{\omega }^2\right)\left({d_2}^2+{\omega }^2\right)}\hfill \\ & & +\frac{R_{co}{A}_j{d_2}^3\left({Le}^2+Le\right)cos\omega tcosaXsin\pi Z}{\left({Le}^2{d_2}^2+{\omega }^2\right)\left({d_2}^2+{\omega }^2\right)}\hfill \\ & & +\sum _{p=0,q=0,r=0}^{\infty }{H}_{pqr}^{\left(j\right)}\left(\omega sinr\omega t+d_{2}cosr\omega t\right)cospaXsinq\pi Z,\hfill \end{array}$

where $w_{pqr}^{\left(j\right)}$, ${\theta }_{pqr}^{\left(j\right)}$, $H_{pqr}^{\left(j\right)},j=1,2,\cdots 7$ are non-linear functions of the amplitudes $A_1$, $A_2$, $A_3$. The term $A_{j}sin\left(\omega t\right)$ represents the lateral oscillation of the roll pattern with amplitude $A_j$ and frequency $\omega$ which is caused by oscillatory instability. Using the fact that the coefficient of every power of $ϵ$ must vanish, the unknown functions $w_{pqr}^{\left(j\right)}$, ${\theta }_{pqr}^{\left(j\right)}$, $H_{pqr}^{\left(j\right)}$ are found by replacing the Equations (21) and (23) in Equation (11). We have obtained the velocity, temperature, and concentration Fourier modes $w_{111}^{\left(1\right)}$, ${\theta }_{111}^{\left(1\right)}$, $H_{111}^{\left(1\right)}$ and $w_{111}^{\left(2\right)}$, ${\theta }_{111}^{\left(2\right)}$, ${\theta }_{202}^{\left(2\right)}$, $H_{111}^{\left(2\right)}$, $H_{202}^{\left(2\right)}$ and $w_{111}^{\left(3\right)}$, $w_{113}^{\left(3\right)}$, $w_{133}^{\left(3\right)}$, $w_{131}^{\left(3\right)}$, ${\theta }_{111}^{\left(3\right)}$, ${\theta }_{113}^{\left(3\right)}$, ${\theta }_{133}^{\left(3\right)}$, ${\theta }_{131}^{\left(3\right)}$, $H_{111}^{\left(3\right)}$, $H_{113}^{\left(3\right)}$, $H_{133}^{\left(3\right)}$, $H_{131}^{\left(3\right)}$ and $w_{244}^{\left(4\right)}$, $w_{242}^{\left(4\right)}$, $w_{224}^{\left(4\right)}$, $w_{222}^{\left(4\right)}$, $w_{133}^{\left(4\right)}$, $w_{131}^{\left(4\right)}$, ${\theta }_{244}^{\left(4\right)}$, ${\theta }_{242}^{\left(4\right)}$, ${\theta }_{224}^{\left(4\right)}$, ${\theta }_{222}^{\left(4\right)}$, ${\theta }_{133}^{\left(4\right)}$, ${\theta }_{131}^{\left(4\right)}$, $H_{244}^{\left(4\right)}$, $H_{242}^{\left(4\right)}$, $H_{224}^{\left(4\right)}$, $H_{222}^{\left(4\right)}$, $H_{133}^{\left(4\right)}$, $H_{131}^{\left(4\right)}$ respectively, at first, second, third and fourth orders.

Figure 1a,b are plotted for ${\theta }_{113}^{\left(3\right)}$ and ${\theta }_{133}^{\left(3\right)}$ (third order temperature modes) with t changing from 0 to 1000 and for fixed values of $\psi$, $Le$ and $\chi$. Here it is noticed that functions ${\theta }_{113}^{\left(3\right)}$ and ${\theta }_{133}^{\left(3\right)}$ are continuous functions of time and depend on the physical parameters, which depict some important changes in the behavior of the flow. Additionally, Figure 1a,b show the projection of limit cycles in the phase plane ${\theta }_{113}^{\left(3\right)}$–${\theta }_{133}^{\left(3\right)}$ for $R=R_{co}$ and $R=3R_{co}$. The oscillatory patterns corresponding to these limit cycles form the standing waves. The time period and amplitude of these limit cycles increase as R increases. For a given R the trajectory oscillates on either side of the origin in the plane (${\theta }_{113}^{\left(3\right)}$, ${\theta }_{133}^{\left(3\right)}$). To attain the maximum magnitude in the negative plane (${\theta }_{113}^{\left(3\right)}$, ${\theta }_{133}^{\left(3\right)}$), the trajectory takes less time, while the opposite behavior of the trajectory is observed for the positive plane. In Figure 1a, at $R=R_{co}$, it is observed that near the origin, the trajectory circulates more as time increases in the positive and negative planes of (${\theta }_{113}^{\left(3\right)}$, ${\theta }_{133}^{\left(3\right)}$). This result shows the occurrence of the onset of chaos. Similarly, for $R=3R_{co}$, the intensity of the chaotic behavior of the flow increases, since the number of circulations of trajectory increases near ${\theta }_{113}^{\left(3\right)}$ = ${\theta }_{133}^{\left(3\right)}$ = 0 (Figure 1b). Figure 2a–c show the vertical velocity eigenfunctions as a function of horizontal length X at $Z=1/2$ and $3R_{co}$ for different values of physical parameters. From these figures, it can be observed that with the increasing values of $\psi$ or $Le$ or $\chi$ the modulus of a maximum value of $w\left(X\right)$ decreases. Figure 3a–c display the vertical velocity eigenfunctions as a function of depth Z at $X=1/2$ and $3R_{co}$ for different values of other physical parameters. These figures exhibit that with increasing values of $\psi$ or $Le$ or $\chi$ the maximum value of $w\left(Z\right)$ decreases. The vertical velocity of the particle is increased up to the first half of the depth and decreased in the second half of the depth due to physical quantities such as acceleration, force and momentum. Thus, the fully developed flow profile occurs when $R\ge R_{co}$. On the hot ($Z=0$) and cold ($Z=1$) walls the vertical velocity $w\left(Z\right)=0$. Each solid line in Figure 3a–c, shows the parabolic shape. The $w\left(Z\right)$ attains its maximum value at the midpoint of the Z-axis, since viscosity is assumed as constant.

4.2. Evaluation of Amplitude $A_1$

From the Equation (26), it can be observed ${\mathcal{N}}_1=0$, which leads to $\mathcal{L}$$\mathit{w}$${}_2$ = 0. The unknown functions $w_2,{\theta }_2$ and $C_2$ are calculated from the Equations (26), (31) and (34), respectively, and which are given by

(40)$\begin{array}{ccc}\hfill & & w_2=A_{2}sin\omega tcosaXsin\pi Z,\hfill \\ & & {\theta }_2=\frac{R_{co}{A}_2}{(d_2+{\omega }^2)}\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z\hfill \\ & & -\frac{1}{8}\frac{R_{co}{A_1}^2\left(4{\pi }^4\left(\omega sin2\omega t+d_{2}cos2\omega t\right)-2{\pi }^2\omega \left(\omega cos2\omega t+d_{2}sin2\omega t\right)\right)sin\left(2\pi z\right)}{\pi \left({d2}^2+{\omega }^2\right)\left(4{\pi }^4+{\omega }^2\right)}\hfill \\ & & -\frac{1}{8}\frac{R_{co}{A_1}^2\left(4{\pi }^{4}d2+d2{\omega }^2\right)sin\left(2\pi z\right)}{\pi \left({d2}^2+{\omega }^2\right)\left(4{\pi }^4+{\omega }^2\right)},\hfill \\ & & C_2=\frac{R_{co}{A}_2\left({Le}^2{d_2}^2\left(\omega sin\omega t+d_{2}cos\omega t\right)+Le{d_2}^2\left(\omega sin\omega t+d_{2}cos\omega t\right)\right)}{\left({Le}^2{d_2}^2+{\omega }^2\right)\left({d_2}^2+{\omega }^2\right)}cosaXsin\pi Z\hfill \\ & & +\frac{R_{co}{A}_2\left({d_2}^2\omega sin\omega t+{\omega }^{3}sin\omega t\right)}{\left({Le}^2{d_2}^2+{\omega }^2\right)\left({d_2}^2+{\omega }^2\right)}cosaXsin\pi Z\hfill \\ & & +K_1{{\mathit{A}}_{\mathit{1}}}^2\left(\omega sin2\omega t+d_{2}cos2\omega t\right)sin2\pi Z+K_2{{\mathit{A}}_{\mathit{1}}}^{2}sin2\pi Z,\hfill \end{array}$

where the coefficients $K_1$ and $K_2$ in Equation (40) depend on a, $\psi$ and $Le$. For $i=3$, Equation (27) gives

(41)$\begin{array}{c}\hfill \left({\mathcal{L}}_1+R_{co}{\mathcal{L}}_2\right)w_3-R_{os}{\mathcal{L}}_2{w}_1={\mathcal{N}}_2,\end{array}$

The solvable condition for Equation (41) determines ${\mathit{A}}_{\mathit{1}}$. Since ${\mathit{A}}_{\mathit{1}}$ is a complicated expression, for the special case $\psi =0,\chi =1$ and $Le=1$ this amplitude is given by

(42)$\begin{array}{c}\hfill A_1=\frac{5\sqrt{-R_{os}\left(158.97tan\left(14.97\right)t-9710.07\right)\left(15.25tan\left(14.96\right)t-15.25\right)}}{15.90tan\left(14.96\right)t-971.00}.\end{array}$

Similarly the unknown functions $w_3,{\theta }_3$ and $C_3$ are obtained from the Equations (27), (32) and (35), respectively and they are given by

(43)$\begin{array}{ccc}\hfill & & w_3=A_3{D}_{1}sin\omega tcosaXsin\pi Z\hfill \\ & & +A_1^3{D}_2\left(\omega sin3\omega t+d_{2}cos3\omega t\right)cosaXsin3\pi Z\hfill \\ & & +A_1^3{D}_3\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin3\pi Z\hfill \\ & & +A_1^3{D}_4\left(\omega sin3\omega t+d_{2}cos3\omega t\right)cosaXsin\pi Z,\hfill \\ & & {\theta }_3={\mathit{A}}_{\mathit{3}}{D}_5\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_6\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin3\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_7\left(\omega sin3\omega t+d_{2}cos3\omega t\right)cosaXsin\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_8\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z\hfill \\ & & +{\mathit{A}}_{\mathit{1}}{\mathit{A}}_{\mathit{2}}{D}_9\left(\omega sin2\omega t+d_{2}cos2\omega t\right)sin2\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_{10}\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z\hfill \\ & & +{\mathit{A}}_{\mathit{1}}{\mathit{A}}_{\mathit{2}}{D}_{11}sin2\pi Z,\hfill \end{array}$

(44)$\begin{array}{c}\hfill \\ & & C_3={\mathit{A}}_{\mathit{3}}{D}_{12}\left(\omega sin\omega t+d_{2}cos\omega t\right)cosaXsin\pi Z\hfill \\ & & +{\mathit{A}}_{\mathit{1}}{\mathit{A}}_{\mathit{2}}{D}_{13}\left(\omega sin2\omega t+d_{2}cos2\omega t\right)sin2\pi Z\hfill \\ & & +{\mathit{A}}_{\mathit{1}}{\mathit{A}}_{\mathit{2}}{D}_{14}sin2\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_{15}\left(\omega sin3\omega t+d_{2}sin3\omega t\right)cosaXsin3\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_{16}\left(\omega sin\omega t+d_{2}sin\omega t\right)cosaXsin3\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_{17}\left(\omega sin3\omega t+d_{2}sin3\omega t\right)cosaXsin\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_{18}\left(\omega sin\omega t+d_{2}sin\omega t\right)cosaXsin\pi Z\hfill \\ & & +{{\mathit{A}}_{\mathit{1}}}^3{D}_{19}\left(\omega sin\omega t+d_{2}sin\omega t\right)cosaXsin\pi Z.\hfill \end{array}$

The coefficients $D_i(i=1to19)$ in Equations (43) and (44) depend on $a,\psi \mathrm{and}Le$ i.e., $D_i=D_i(a,\psi ,Le)$. Proceeding with the procedure as above and by substituting $i=4$ in the Equation (27) we have obtained $A_2=0$ at the fourth order of $ϵ$. As the order of $ϵ$ increases, calculations become more lengthy. By substituting $i=5$ in the Equation (27) and applying the solvability condition, we have computed the amplitude $A_3.$ Since the $A_3$ expression is a lengthy, it is not shown here to conserve space.

5. Convective Heat Transport

The heat transfer coefficient in terms of the local Nusselt number ($N_L$) is defined by [32,33]:

(45)$\begin{array}{c}\hfill N_L=\frac{\partial T}{\partial n}\end{array}$

where n denotes the normal direction on a plane. The average heat transport is measured by the time-dependent averaged Nusselt number ($Nu$), that depends on t and is given by

(46)$\begin{array}{c}\hfill Nu=\overline{wT}-\frac{\partial \overline{T}}{\partial Z},\end{array}$

here the bar indicates the horizontal mean. The Equation (46) defined on the hot wall ($Z=0$) as [25]

$\begin{array}{c}\hfill Nu=\overline{wT}-\frac{\partial \overline{T}}{\partial Z}=\frac{1}{L}{\int }_0^L\left(wT-\frac{\partial T}{\partial Z}\right){|}_{Z=0}dX,\end{array}$

where L represents the normalized horizontal width of the cell.

5.1. Local Nusselt Number ($N_L$)

Figure 4a–d represent the influence of the physical parameters R, $\psi$, $Le$ and $\chi$ on $N_L$ with respect to X for a fixed value of t. The number of peaks of $N_L$ depend on the variation of R only. The location of minimum and maximum values of $N_L$ depends on $R,\psi ,Le$ and $\chi$. At the onset $R=47.65$, the distribution line of $N_L$ is parallel to X-axis (constant) and for $R>47.65$, the distribution lines of $N_L$ attain the maximum values depending on R, but not on X. The maximum values of $N_L$ increase as R increases. This outcome demonstrates the development of heat flow ($N_L$) on the hot wall. The number of peaks is unaffected by changes in the physical parameters $\psi$ or $Le$ or $\chi$. The location of a maximum value of $N_L$ depends on $\psi$ or $Le$ or $\chi$ and also on X. The individual effect of these parameters is same on $N_L$, that is as $\psi$ or $Le$ or $\chi$ decreases the maximum value of $N_L$ increases.

5.2. Time-Dependent Averaged Nusselt Number ($Nu$)

Let $N{u}^{\left(2\right)}(s=1)$ and $N{u}^{\left(4\right)}(s=2)$ denote the approximations of the Nusselt number in the second and fourth orders, respectively, and they are given as

(47)$\begin{array}{ccc}\hfill & & N{u}^{\left(2\right)}-1=\frac{{A_1}^2{ϵ}^2{R}_{co}\left(4{\pi }^2\omega \left(2{\pi }^2+d_2\right)sin2\omega t+2{\pi }^2\left(8{\pi }^2{d}_2-{\omega }^2\right){cos}^2\omega t\right)}{\left(32{\pi }^4+8{\omega }^2\right)\left(d_2^2+{\omega }^2\right)}\hfill \\ & & +\frac{{A_1}^2{ϵ}^2{R}_{co}\left(-4{\pi }^4{d}_2+{\omega }^2\left(4{\pi }^2+d_2\right)\right)}{\left(32{\pi }^4+8{\omega }^2\right)\left(d_2^2+{\omega }^2\right)},\hfill \end{array}$

and

(48)$\begin{array}{ccc}\hfill & & {N{u}^{\left(4\right)}}_{at\omega =0}=-3.142\times {10}^{-12}{{\mathit{A}}_{\mathit{1}}}^4+6.283\times {10}^{-12}\frac{{A_1}^4{R}_{co}}{R}-3.142\times {10}^{-12}\frac{{A_1}^4{R_{co}}^2}{R^2}\hfill \\ & & +0.258R_{os}{A_1}^2-0.516\frac{{A_1}^2{\mathit{R}}_{\mathit{os}}{R}_{co}}{R}+0.258\frac{{A_1}^2{R}_{os}{R_{co}}^2}{R^2}+0.375\frac{{\mathit{R}}_{\mathit{co}}{A_1}^2}{{\pi }^2+a^2}\hfill \\ & & -0.375\frac{{R_{co}}^2{A_1}^2}{\left({\pi }^2+a^2\right)R}+45.199A_1{A}_3-90.397\frac{A_1{A}_3{R}_{co}}{R}+45.19832981\frac{A_1{A}_3{R_{co}}^2}{R^2}\hfill \\ & & +2.227-2.454\frac{R_{co}}{R}+1.227\frac{{R_{co}}^2}{R^2}.\hfill \end{array}$

The relation between $Nu$ and R can be expressed by a power law $Nu\propto R^r$ as $Nu=a{(R/R_{co})}^r$. From Table 1 it can be seen the exponent r lies between 0.21–0.384, for $-0.03<\psi <-0.01$. The results show that the power law exponent r increases as $\psi$ decreases. In Table 2, the power law relation $NuvsR/R_{co}$ is shown, for $0.04<Le<0.06$ and r lies between 0.38–0.43. This result shows that r increases with decreasing $Le$. Table 3 shows the power law relationship between $Nu$ and $R/R_{co}$ for $0.2<\chi <0.4$, with r ranging between 0.2–0.38. This result shows that r increases with decreasing $\chi$.

Figure 5a–c show the dependence of $Nu$ on R for different values of $\psi ,Le$ and $\chi$. These figures show as the physical parameters $\psi ,Le$ and $\chi$ increase, the value of $Nu$ on the hot wall decays. Thus the presence of a physical parameters suppresses the heat transfer rate, and accordingly, the boundary layer thickness increases on the hot wall. Whereas the augmenting values of R favor the appearance of convection and turbulence. The same physical interpretation can be explained from the power law equations. Figure 5a–c show the variation of $Nu$ with respect to R for different values other parameters. The corresponding power law equations are presented in the third and fourth columns of Table 1, Table 2 and Table 3. From the power law equations, it can be predicted that as the power law exponent increases the intensity of the fluid motion as well as the heat transfer rate are increased and the nature of the convection cell becomes more complex. Figure 5d is plotted for $Nu$ against R with different values of t. This figure shows that, as time increases for a constant value of R (constant temperature gradient), $Nu$ increases, i.e., as time increases, the flow increases on the hot wall.

Figure 6 shows the time variation of $Nu$ for distinct values of $\psi$ and for fixed values of other remaining parameters. From Figure 6a it can be observed that for $\psi =-0.03$ and $R\approx R_{co}$, when $t=0$, $Nu$ takes the large values and as t increases $Nu$ oscillates with decreasing amplitude. Figure 6b is plotted for $\psi =-0.02$ and $R\approx R_{co}$. This figure shows that $Nu$ oscillates with two different amplitudes periodically as t increases. Indeed, during the oscillation, $Nu$ decreases to a value close to $1.096$ when $\psi =-0.01$ and $R=41.5137494$ with calculated $R_{co}=41.47227710$ (Figure 6c). From Figure 6a–c, it is noted that with a negative reduction of $\psi$, the oscillatory frequency and amplitude decreased, and the time of the oscillatory period increased. In Figure 6d, the $Nu$ approximations, with respect to time t, for the amplitudes $A_1(s=1)$ and $A_3(s=2)$ are plotted. The amplitudes $A_1$ and $A_3$ are obtained at the orders $O\left({ϵ}^3\right)$ and $O\left({ϵ}^5\right)$, respectively. The time-dependent averaged Nusselt number ($Nu$) showed a similar trend to that obtained by Brand and Steinberg [9]. The present results (solid lines) are compared with those results obtained by Brand and Steinberg [9]. At the $O\left({ϵ}^3\right)$, the approximate average $Nu$ shows a similar trend as those obtained by Brand and Steinberg [9]. Figure 6e shows the dependence of vertical velocity on time, t and R with fixed values of $Le=0.04$, $\psi =-0.03$ and $\chi =0.02$. This figure exhibits, as t and R increase the amplitude and the frequency of the oscillations increase.

Changes in kinetic energy (KE) and potential energy (PE) related to $R/R_{co}$ is shown in Figure 7a–d for distinct values of $\psi ,Le,\chi$ and t, respectively. From these figures, it can be observed that KE is equal to PE when $R=R_{co}$. For $R/R_{co}>1$, initially, KE dominates PE, but this trend is reversed after a certain value of $R/R_{co}$. The total energy increases for all increasing values of $R/R_{co}$. The KE and PE decrease as $\psi$ or $Le$ or $\chi$ increases. From Figure 7a–c, the total energy of the system can be observed to increase as $\psi$ or $Le$ or $\chi$ decreases. The time history of total energy is explained in Figure 7d. The results show that the total energy is enhanced as time increases and KE always dominates PE. The influence of all physical parameters on the total energy and the rate of heat transfer is observed to be similar.

6. Distortion of Streamlines and Isotherms

The fluid flow displayed using the stream function $\Psi (X,Z)$ is obtained from the u and w velocity components. For two-dimensional flow, the relationship between $\Psi (X,Z)$ and velocity components is given by [34]

$u=-\frac{\partial \Psi }{\partial Z}\mathrm{and}w=\frac{\partial \Psi }{\partial X},$

which yields the following single equation and by utilizing the accompanying limit conditions for $\Psi$ and removing the arbitrary integration constants (Kuo [23]):

$\Psi ={\nabla }^2\Psi =0\mathrm{at}X=0,\pi /a\mathrm{and}Z=0,1.$

The family of streamlines gives a clear view of the entire flow domain and its important characteristics. Whereas the points with uniform temperature connected with lines that are called as the isotherms. The snapshots of the flow field and the hot regions close to the beginning of oscillatory convection are shown in terms of streamlines and isotherms, respectively, for distinct values of R, $\psi$, $Le$, $\chi$ and t.

The combined heat and mass transfer phenomena are extended to the convective heat transfer principle of heatlines along with the graphical presentation of convective mass transfer. A representative set of streamlines, isotherms and isoconcentration lines is presented in Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 for the variation in R, $\psi$, $Le$, t and $\chi$. The circulation strength is assumed to be positive in the anticlockwise direction and negative in the clockwise direction. The snapshot of streamlines, isotherms and isoconcentration lines for increasing values of R with $Le=0.06$, $\psi =-0.01$, $\chi$= $0.4$ and $t=0.1$ are displayed in Figure 8a–l. The fluid flow is represented at the beginning by cellular patterns (Figure 8a). The cell lying in $0\le X\le 1$ (Figure 8a) has the absolute maximum and minimum circulation strength as $0.008$ and $-0.008$, respectively. From Figure 8b, near the start of oscillatory convection, isotherms depict the effect on the temperature field with absolute maximum and minimum values $-0.00986$ and $-1.01015$, respectively. Isotherms are observed to be nearly horizontal parallel lines, suggesting that conduction carries much of the heat transfer. Figure 8c is also a plot of the isoconcentration lines that are calculated when $R=R_{co}$, and the absolute maximum and minimum values of these isoconcentration lines are $0.3063$ and $-1.3066$, respectively. Figure 8d–f illustrate streamlines, isotherms and isoconcentration lines, respectively, for $R=3R_{co}$ for the same values of other physical parameters as those defined in Figure 8a–c. The temperature, concentration and gravitational buoyancy forces act together and change the fluid flow pattern, as shown by the streamlines (Figure 8d). The bicellular patterns become slightly deformed (Figure 8d) and, for a cell lying between $0\le X\le 1$, the absolute maximum and absolute minimum of circulation strength are $0.2168$ and $-0.21683$, respectively. The isotherms in Figure 8e unambiguously indicate the presence of a convective mode of heat flow inside the fluid with absolute maximum and minimum values as $-0.00537$ and $-1.0256$, respectively. Figure 8f plots the isoconcentration lines, indicating the absolute maximum and minimum values as $159.461$ and $-160.497$, respectively. Figure 8g, is plotted for $R=5R_{co}$. As the thermal buoyancy increases, two small vortices, namely B and $B^{\prime }$, near the top right and left bottom plates are generated with the equal circulation strength of $0.026$. The original cell becomes more deformed as R changes from $R=3R_{co}$ to $R=5R_{co}$ (Figure 8g). For a cell lying between $0\le X\le 1$, the absolute maximum circulation strength is $0.20509$ and the minimum is $-0.20512$. The isotherms in Figure 8h exhibit, inside the fluid layer, the convective mode of heat transfer creates the oval-shaped circular patterns in the middle region with absolute maximum and minimum values of $0.35803$ and $-1.35824$, respectively. Figure 8i is a plot of the isoconcentration lines and their the absolute maximum and minimum values are $264.976$ and $-266.014$, respectively. Figure 8j, is plotted for $R=7R_{co}$. When R is increased from $5R_{co}$ to $7R_{co}$, due to the increased thermal buoyancy, the basic cell encounters more deformation, as shown in Figure 8j. Additionally, the two small generated vortices, B and $B^{\prime }$, are grown in size and shape with the increased circulation strength of $0.110$. Thus, the basic cell divided into two more new vortices, namely A and $A^{\prime }$ near the top and bottom of the basic cell with the same circulation strength of $-0.241$. The isotherms that are shown in Figure 8k demonstrate the strong convective heat flow in the fluid layer, with the absolute maximum and minimum values of $0.86877$ and $-1.86911$, respectively, which produces the circular oval patterns in the middle of the domain. Figure 8l plots the isoconcentration lines, with the absolute maximum and minimum values of $364.572$ and $-365.63$, respectively. From Figure 8a–l, it is noted that when R is increased from $R_{co}$ to $7R_{co}$, oscillatory flow moves very quickly and the transfer of mass increases due to thermal buoyancy.

Figure 9 illustrates the streamlines, isotherms and isoconcentrations for distinct values of $\psi$ with a fixed set of other parameters, say, $R=3R_{co}$, $Le=0.06$, $\chi =0.4$ and $t=0.1$. Figure 9a–c are plotted for $\psi$ = $-0.02$. In Figure 9a the primary concentric roll patterns of streamlines are tilted towards the right side by retaining the flow orientation in the considered range of $0\le X\le 1$. The absolute maximum and minimum values of circulation strength are $0.03585$ and $-0.03585$, respectively. Figure 9d–f are plotted for $\psi$ = $-0.03$, with the same fixed set of parameters and contour lines defined in Figure 9a–c. In Figure 9d, the flow pattern in the range $0\le X\le 1$ becomes more deformed by creating the two additional vortices B and $B^{\prime }$ at the top and bottom limits, respectively, and on either side of the primary cell. The circulation strength of the two vortices, B and $B^{\prime }$ is $-0.012$. In the area, $0\le X\le 1$, circulation strength has absolute maximum and minimum values of $0.03556$ and $-0.03555$, respectively. From Figure 9a,d, we can observe that the streamline deformation is increased when $\psi$ values decrease. This implies that the effect of an increase in $\psi$ values stabilizes the convective system. The flow of heat transfer is depicted using the isotherms and is shown in Figure 9b,e for $\psi$ = $-0.02$ and $-0.03$, respectively. The isotherms appeared as plumes in the range of $0\le X\le 5$, as shown in Figure 9b. The boundary layer thickness increased when $\psi$ decreased from $-0.02$ to $-0.03$. The isoconcentrations (concentration field) are plotted in Figure 9c,f for $\psi =-0.0$2 and $-0.03$, respectively. The isoconcentrations occurred in the form of vortices in the range of $0\le X\le 5$. The size of vortex structures increased as $\psi$ decreased. The boundary layer of the concentration region increases as $\psi$ decreases from $-0.02$ to $-0.03$. For a particular $\psi$ value, the distribution of the concentration field is higher than the temperature field due to Soret effect.

The distributions of velocity, temperature and concentration fields are plotted in Figure 10a–f for distinct values of $Le$ with $R=3R_{co}$, $\psi =-0.01$, $\chi =0.4$ and $t=0.1$. In Figure 10a, the streamlines are plotted for $Le=0.05$. In the range of $0\le X\le 1$ the basic primary cell has been deformed and two new vortices, namely, B and $B^{\prime }$, are generated at the top left and right bottoms, respectively, of the basic primary cell with a circulation strength of $-0.027$. The absolute maximum circulation strength of the basic cell, namely A, is 0.036. In Figure 10d, the streamlines are plotted for the decreased value of $Le$ i.e., $Le$ = $0.04$, with the same parameters defined to plot the contour lines as shown in Figure 10a. The generated two vortices, namely, B and $B^{\prime }$ in Figure 10a become enlarged with the increased circulation strength of $-0.006$, as shown in Figure 10d. The primary cell is split into two vortices, namely, A and $A^{\prime }$, with the absolute maximum circulation strength of the basic cell $0.04$. As $Le$ decreases, the thermal diffusivity increases and accordingly, the mass diffusivity decreases. As $Le$ decreases from $0.05$ to $0.04$, the streamlines become more deformed. This implies that the decreasing $Le$ destabilizes the convective system.

The flow of heat transfer is depicted by using the isotherms and is shown in Figure 10b,e for $Le$ = $0.05$ and $0.04$, respectively. It is observed that the isotherms in the hot regions get enlarged as $Le$ decreases. Hence, we can conclude that, with the decreasing values $Le$, the system becomes unstable. Figure 10c is a plot of the isoconcentration lines for $Le$ = $0.05$ and the absolute maximum and minimum values of isoconcentrations are $118.373$ and $-119.436$, respectively. Figure 10f is a plot of the isoconcentration for $Le$ = $0.04$, and the absolute maximum and minimum values of isoconcentrations are $81.9809$ and $-82.9856$, respectively. Figure 10c,f show $Le=0.05$ and $0.04$, respectively. The size of the vortex patterns decreases as $Le$ decreases. Finally, as $Le$ decreases, the thickness of thermal boundary layer increases, and the boundary layer thickness of the concentration field decreases. Hence, the increasing values of $Le$ stabilize the system.

Figure 11 illustrates the pattern of streamlines, isotherms and isoconcentration lines for distinct values of t with $R=3R_{co}$, $\psi =-0.01$, $\chi =0.4$ and $Le=0.06$. Figure 11a represents the pattern of streamlines for $t=1.1$ and have the absolute maximum and absolute minimum values of circulation strength $1.07184$ and $-1.07222$, respectively, in the range of $0\le X\le 1$. The streamlines are of parallel concentric roll patterns, thus indicating the stable fluid flow at $t=1.1$. Figure 11d shows the pattern of streamlines for t = $1.6$ for the same values of other physical parameters that are considered in Figure 11a. The absolute maximum and absolute minimum values of circulation strength are $0.2809$ and $-0.28082$, respectively. The stable primary cells have been deformed and two new vortices are generated inside the primary cell, namely, B and $B^{\prime }$, with an absolute minimum value of circulation strength as $-0.257$. The isotherms plotted in Figure 11b,e are for the same physical parameters that are considered in Figure 11a,d, respectively. Figure 11b shows the isotherms that appeared in the form of plumes in the middle region of the temperature field. As t increases from $1.1$ to $1.6$, the size of the plume structure increases vertically. This result suggests that the heat transfer is enhanced as time increases for the fixed values of the physical parameters. Figure 11c is a plot of the isoconcentration lines for t = $1.1$ and the absolute maximum and minimum values of isoconcentrations are $203.197$ and $-204.238$, respectively. Figure 11f plots the isoconcentration lines for the same values of physical parameters as those defined in Figure 11c and for t = $1.6$. The absolute maximum and minimum values of isoconcentrations are $290.001$ and $-291.049$, respectively. As t increases from $0.1$ to $1.6$, the convective mass transport of the fluid increases. Finally, the boundary layer thickness of both the temperature and concentration fields decrease as time increases.

Figure 12 shows the streamline patterns, isotherms and isoconcentration lines for distinct fixed values of $\chi$ with $R=3R_{co}$, $\psi =-0.01$, $t=0.1$ and $Le=0.06$. Figure 12a, illustrates the streamlines for $\chi$ = $0.3$. The primary basic cell has been deformed and is pushed into the middle of the region by the two generated vortices, namely, B and $B^{\prime }$ near the top and bottom boundaries, respectively, in the considered range of $0\le X\le 1$ with the circulation strength of $-0.015$. In Figure 12d, the streamlines are drawn for the same values of physical parameters as defined in Figure 12a and $\chi$ = $0.2$. In the considered range of $0\le X\le 1$, the two vortices B and $B^{\prime }$, shown in Figure 12a, increased in size with the circulation strength of $-0.007$ (Figure 12d). The basic cell is split into two new vortices, namely, A and $A^{\prime }$, with a maximum circulation strength of $0.004$. This implies that the decreasing $\chi$ destabilizes the convective system. The flow of heat transfer is depicted through isotherms and is shown in Figure 12b,e for $\chi$ = $0.3$ and $0.2$, respectively. It can be observed that, for the decreased values of $\chi$, the hot regions increased. Hence, as $\chi$ decreases, the system becomes unstable. Figure 12c is a plot of the isoconcentration lines for $\chi$ = $0.3$, the absolute maximum and minimum values of isoconcentrations are $143.211$ and $-144.223$, respectively. For $\chi$ = $0.2$, in Figure 12f, the plot of the isoconcentrations are shown with the absolute maximum and minimum values are $125.82$ and $-126.828$, respectively. This figure reveals that, as $\chi$ decreases from $0.4$ to $0.2$, mass transfer decreases due to less permeability of the porous medium.

Topology of Flow

Fluid flow topology is explored based on the Euler number, ${\zeta }^{{}^{\prime }}$. As described by Jana et al. [35], on the surface, ${\zeta }^{{}^{\prime }}$ is the aggregate of the Poincare indices of the critical points. The invariance of topology relation is defined by

(49)$\begin{array}{c}\hfill NE-(NH+\frac{1}{2}NP)={\zeta }^{{}^{\prime }},\end{array}$

where the number of elliptic, hyperbolic and parabolic points are denoted by $NE$, $NH$ and $NP$, respectively [36,37]. The vorticity contours for $R\simeq R_{co}$, $\psi$ = −0.01 and $Le$ = 0.06 are shown in Figure 13a. For a fixed interval of X, topological rule given in Equation (49) is satisfied with $NE=2,NH=2$ and $NP=0$. A similar study has been carried out for vorticity contours with $R=3R_{co}$ depicted in Figure 13b. For $R=3R_{co}$, Equation (49) is also satisfied by the vorticity contours as exhibited in Figure 13b by $NE$ = 8, $NH$ = 8 and $NP$ = 0. The vorticity contours appeared in the form of vortices (elliptical shape) for $R=R_{co}$, and they are known as primary contours. Each pair of consecutive vortices has the counter-rotating form. The value of vorticity on a line is constant. The value of vorticity decreases as it moves away from the center of a primary cell, i.e., the circulation of fluid is less near the center. When R increased from $R_{co}$ to $3R_{co}$, two new vortices showed their presence at the top and bottom walls and on either side of the primary vortex. Due to the presence of these new vortices, the primary vortex became deformed and, inside the primary vorticity contour, two new vortices emerged. Thus, the strength of vorticity and the velocity of the fluid increases as R increases.

7. Heatfunction

The convective heat transfer phenomenon is represented by heatlines in the present study. Morega and Bejan [38] successfully used the concept of heatlines to visualize the fluid flow. Various researchers [39,40,41,42,43,44,45,46] enlarged this concept to a wide range applications of natural convective systems. The nondimensional heatfunction, $H^{*}$ is defined as

(50)$\begin{array}{c}\hfill \frac{\partial H^{*}}{\partial X}=wT-\frac{\partial T}{\partial Z},\end{array}$

(51)$\begin{array}{c}\hfill -\frac{\partial H^{*}}{\partial Z}=uT-\frac{\partial T}{\partial X},\end{array}$

where $T=T_s+\theta$ and $T_s=T_0-Z$ denotes the static temperature. Here $T_0$ denotes the arbitrary reference temperature. It should be noticed that different values of $T_0$ lead to different slopes of $H^{*}$[40,41]. The $H^{*}$ defined in Equations (50) and (51) is similar to the one defined by Kimura and Bejan [42] but with opposite signs. Eliminating the temperature gradients from the above Equations (50) and (51) by cross differentiation, the following Poisson-type equation is obtained:

(52)$\begin{array}{c}\hfill \frac{{\partial }^2{H}^{*}}{\partial X^2}+\frac{{\partial }^2{H}^{*}}{\partial Z^2}=\frac{\partial \left(wT\right)}{\partial X}-\frac{\partial \left(uT\right)}{\partial Z}.\end{array}$

The following boundary conditions for $H^{*}$ are obtained by integrating the Equations (50) and (51):

(53)$\begin{array}{c}\hfill atZ=0,0\le X\le \pi /a:H^{*}(X,0)=H^{*}(0,0)+{\int }_0^X\left(wT-\frac{\partial T}{\partial Z}\right)dX,\end{array}$

(54)$\begin{array}{c}\hfill atZ=1,0\le X\le \pi /a:H^{*}(X,1)=H^{*}(0,1)+{\int }_0^X\left(wT-\frac{\partial T}{\partial Z}\right)dX,\end{array}$

(55)$\begin{array}{c}\hfill atX=0,0\le Z\le 1:H^{*}(0,Z)=H^{*}(0,0)-{\int }_0^Z\left(uT-\frac{\partial T}{\partial Z}\right)dZ,\end{array}$

(56)$\begin{array}{c}\hfill atX=\pi /a,0\le Z\le 1:H^{*}(\pi /a,Z)=H^{*}(\pi /a,0)-{\int }_0^Z\left(uT-\frac{\partial T}{\partial Z}\right)dZ.\end{array}$

The integral expressions that exist in the boundary conditions (53)–(56) represent the average Nusselt numbers on the corresponding boundaries. The total heat transfer across the system in terms of $H^{*}$ is plotted graphically in Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18.

Figure 14a–d show the pattern of heatlines for $R\simeq R_{co}$, $R=$$3R_{co}$, $5R_{co}$ and $7R_{co}$ respectively, for $\psi =-0.03$, $Le=0.04$, $t=0.8$ and $\chi =0.2$. From Figure 14a, the heatline contours within the domain are found to be normal to $Z=0$ and $Z=1$ lines due to conduction dominant heat transfer. The absolute maximum and absolute minimum values of heatlines are $4.88944$ and $-0.01944$, respectively, in the considered range $0\le X\le 5$. In the neighborhood of the $X=0$ line, the heatlines are vertically parallel. The curvature at the central part of each heatline increases with X and the heatlines are similar to that of parabolic structure. It represents the nonlinear propagation of heat transfer that occurs at $R\simeq R_{co}$. Hence the transition takes place from the conduction to the convection state at $R\simeq R_{co}$. Figure 14b shows the heatlines for $R=3R_{co}$ and the absolute maximum and minimum values of heatlines are $4.79219$ and $-0.36944$, respectively. The heatlines shown in Figure 14a are converted to curves and closed paths as depicted in Figure 14b when $R\simeq R_{co}$ increased to $R=3R_{co}$. Figure 14c illustrates the heatlines for $R=5R_{co}$ and have the absolute maximum and minimum $4.99725$ and $-0.63857$, respectively. Figure 14d shows the heatlines for $R=7R_{co}$, which have the absolute maximum and minimum $5.21091$ and $-0.8888$, respectively. This figure shows that some of the heatlines shown in Figure 14c are divided into closed paths (vortices) and curves. It can also be noted from these figures that the size of closed paths increases as R increases. Thus, from Figure 14a–d, we can conclude that the heat transfer intensity increases with R.

Looking at Figure 15a,b, the influence of $\psi$ is analyzed on the heat transfer for fixed values of other physical parameters, $R=3R_{co}$, $Le=0.04$, $t=0.8$ and $\chi =0.2$. For $\psi =-0.02$, the heatlines drawn in Figure 15a, have the absolute maximum and the absolute minimum values $4.79076$ and $-0.19418$, respectively, in the assumed range $0\le X\le 5$. Figure 15b shows the heatlines for $\psi =-0.01$ and the absolute maximum and absolute minimum values are $4.89789$ and $-0.01382$, respectively. Figure 15a shows the pattern of heatlines, namely, the combination of curves and vortices with vortices located in the central part of the considered domain. Figure 15b shows that the deformation of the curved lines and the size of the vortices are reduced compared to those lines shown in Figure 15a. Thus, Figure 15a,b show the decreasing intensity of heat transfer as $\psi$ increases from −0.02 to −0.01.

Figure 16a,b illustrate the visualization of heat flow via heatlines for different $Le$ and for the fixed values of other physical parameters, $R=3R_{co}$, $\psi =-0.03$, $t=0.8$ and $\chi =0.2$. Figure 16a illustrates heatlines for $Le=0.05$ and has the absolute maximum and the absolute minimum $4.77309$ and $-0.32969$, respectively. Figure 16b shows heatlines for $Le=0.06$ and have the absolute maximum and the absolute minimum $4.76157$ and $-0.29866$, respectively, in the considered range of X. The heat flow patterns of Figure 16a,b are similar to those in Figure 15a,b, respectively. The deformation of the curved lines and the size of vortices decreases as $Le$ increases from 0.05 to 0.06. Thus, as $Le$ increases the thermal diffusivity decreases.

The snapshot of heatlines for fixed $R=3R_{co}$, $Le=0.04$, $\psi =-0.03$ and $\chi =0.2$ for different t are shown in Figure 17a,b. Figure 17a shows the pattern of heatlines for $t=0.1$, which have a structure similar to that of heatlines for $R\simeq R_{co}$ with the absolute maximum and absolute minimum being $4.9932$ and $-0.00544$, respectively. Figure 17b illustrates the heatlines for $t=0.5$. The absolute maximum and the absolute minimum are $4.78354$ and $-0.23072$, respectively, in the considered range of X. Some heatlines shown in Figure 17a are converted into curves and closed paths as t increases. Thus, as t increases the intensity of heat flow also increases.

Figure 18a,b show the visualization of heat flow via heatlines for different values of $\chi$ and for fixed values of other parameters $R=3R_{co}$, $\psi =-0.03$, $Le=0.04$ and $t=0.8$. Figure 18a illustrates heatlines for $\chi =0.3$ and have the absolute maximum and absolute minimum $4.71885$ and $-0.19442$, respectively. Figure 18b illustrates heatlines for $\chi =0.4$ with the absolute maximum and absolute minimum $4.71914$ and $-0.1101$, respectively, in the considered range $0\le X\le 5$. Figure 18a,b show that, as $\chi$ decreases, the intensity of the heat flow increases.

8. Entropy Generation

It is important to understand the efficient use of energy resources, together with the minimum degradation of energy, i.e., the minimum generation of entropy should be efficient during heat transfer. In the present work, entropy generation is related to the irreversible existence of heat transfer and viscous effects within the fluid and the fluid-solid interfaces. The non-dimensional form of the local entropy generation rate ($S_{gen}$) is calculated by the addition of the heat transfer irreversibility ($S_{HTI}$) and the fluid friction irreversibility ($S_{FFI}$) as discussed in [47,48,49] and is expressed as:

(57)$\begin{array}{c}\hfill S_{gen}=S_{HTI}+S_{FFI},\end{array}$

(58)$\begin{array}{c}\hfill S_{HTI}={\left(\frac{\partial T}{\partial X}\right)}^2+{\left(\frac{\partial T}{\partial Z}\right)}^2,\end{array}$

(59)$\begin{array}{c}\hfill S_{FFI}=\varphi \left[u^2+w^2+Da\left\{2{\left(\frac{\partial u}{\partial X}\right)}^2+2{\left(\frac{\partial w}{\partial Z}\right)}^2+{\left(\frac{\partial u}{\partial Z}+\frac{\partial w}{\partial X}\right)}^2\right\}\right],\end{array}$

where $u,$w and T are the known non-dimensional velocity and temperature functions given in Section 2 and $\varphi =\left(\mu T_0^2{\nu }^2\right)/(K^4{\beta }_1^2{g}^2{d}^2$) is the irreversibility distribution ratio, which is initially minimized numerically (fixed as ${10}^{-2}$) and the physical parameter $Da$ was fixed as $1.2\times {10}^{-5}$. Al-Hadhrami et al. [50] proposed three different models, namely, Darcy, Darcy–Freshman and Brinkman–Brinkman, for viscous dissipation flow through porous media. Hooman and Gurgenci [51] compared these three models. They concluded that these models were effectively showed the same results only for small values of $Da$. The viscous dissipation model proposed by Baytas [52] is employed to obtain Equation (57). Figure 19a–f show the results of entropy generation for different R. At $R=R_{co}$, the pattern of $S_{HTI}$ is the combination of vortices and curves. The vortices are converted into the wavy pattern as R increases. At $R=R_{co}$, the pattern of $S_{FFI}$ appeared in the form of vortices for $R=R_{co}$ (Figure 19b). For $R=5R_{co}$, the above pattern is converted into the combination of elliptical and deformed vortices. Each deformed vortex encloses five other vortices (Figure 19d). When R is increased to $9R_{co}$, the patterns become more deformed, as shown in Figure 19f. It is observed that, for $R\simeq R_{co}$ (Figure 19a,b), the maximum entropy generations due to $S_{HTI}$ and $S_{FFI}$ are 0.03098 and 0.00032, respectively. The maximum entropy generation due to $S_{HTI}$ is 3.69005 when $R=5R_{co}$ and 18.6827 when $R=9R_{co}$ (Figure 19c,e). Similarly, the maximum entropy generation due to $S_{FFI}$ is 0.5214 when $R=5R_{co}$ and 4.86893 when $R=9R_{co}$, i.e., for small values of R, the $S_{HTI}$ and $S_{FFI}$ are small and they increase with R. This is because, at low R, the viscous force dominates and the heat transfer is mainly due to conduction. As already mentioned, the $S_{gen}$ is the sum of the two contributions, namely, $S_{HTI}$ and $S_{FFI}$. Thus, it can be concluded that $S_{gen}$ increases as R increases.

9. Conclusions

Theoretically, near and above the onset of oscillatory convection, the nonlinear state of the Rayleigh–Bénard system of a superfluid mixture, ${}^{3}He$–${}^{4}He$ kept in a porous medium was investigated in the present study. As functions of the control parameters relative to ${}^{3}He$–${}^{4}He$, the oscillation frequency, the linear growth rate of instability and the spatial flow pattern were studied. The results are visualized for the study of the essence of fluid flow and heat transfer rate by employing streamlines, isotherms, heatlines, kinetic energy, potential energy and entropy production due to heat transfer and fluid friction irreversibility. The results revealed that the maximum value of the local Nusselt number ($N_L$) increases as the Rayleigh number (R) and time (t) increase. As the separation ratio ($\psi$), Lewis number ($Le$) and porosity number ($\chi$) decrease, the maximum value of $N_L$ increases. From time-dependent averaged Nusselt number ($Nu$) calculations, it is noted that the heat flux occurs more as $\psi$, $Le$ and $\chi$ decrease. The total energy of the system increases as $\psi$, $Le$ and $\chi$ decrease. By observing the streamline patterns, isotherms, isoconcentration lines and heatlines patterns, it can be said that the flow field encountered more deformation and the intensity of heat flow increased as R and t increased in comparison with other parameters. It is also observed that, as $\psi$, $Le$ and $\chi$ decrease, the flow field encounters more deformation and the intensity of heat flow increases. The streamlines are of parallel concentric rolls near the onset but, as $\psi$ increases, the absolute maximum and minimum values of the circulation strength of rolls decrease and increase, respectively. At the onset of oscillatory convection, the isotherms are of almost horizontal and parallel lines. Entropy generation is observed to increase with R.

Author Contributions

Conceptualization, Y.R. and G.S.; methodology, Y.R.; software, D.L.; validation, Y.R., G.S. and D.L.; writing—original draft preparation, G.S.; writing—review and editing, Y.R. and D.L.; funding acquisition, D.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No applicate.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

$A_i$	Amplitude
a	Wavenumber
$C^{{}^{\prime }}$	Concentration field
$C_p$	Specific heat at constant pressure
$C_s^{{}^{\prime }}$	Static concentration
D	Diffusion coefficient for particle diffusion
$Da$	Darcy number
d	Depth of horizontal layer
$\overrightarrow{g}$	Gravitational field
H*	Heatfunction
K	Permeability of porous medium
$Le$	Lewis number
N	Rate of transport of solute
$N_L$	Local Nusselt number
$Nu$	Time-dependent averaged Nusselt number
$N{u}^{\left(2\right)},N{u}^{\left(4\right)},N{u}^{\left(6\right)}$	Approximate Nusselt numbers of second,
	fourth and sixth orders, respectively
$q_{cs}$	Critical wavenumber
p	Pressure
R	Rayleigh number
$R_{cs}$	Critical Rayleigh number of stationary convection
$R_o$	Rayleigh number of oscillatory convection
$R_s$	Rayleigh number of stationary convection
$T^{{}^{\prime }}$	Temperature field
$T_0^{{}^{\prime }}$	Reference temperature
$T_s^{{}^{\prime }}$	Static temperature
${\overrightarrow{V}}^{\prime }=(u^{{}^{\prime }},v^{{}^{\prime }},w^{{}^{\prime }})$	Velocity field
${\overrightarrow{V}}_s^{\prime }$	Static velocity field
$t^{\prime }$	Time
$x^{\prime },y^{\prime },z^{\prime }$	Cartesian coordinates
Greek symbols
$\Delta T^{{}^{\prime }}$	Temperature gradient between top and bottom layers
$\beta$	Adverse temperature gradient
$\kappa$	Thermal diffusivity
${\kappa }_T$	Thermal diffusion ratio
$\omega$	Frequency of oscillations
${\theta }^{{}^{\prime }}$	Perturbed temperature
$\rho$	Density
${\rho }_0$	Reference mass density
$\nu$	Kinematic viscosity
${\beta }_1$	Thermal expansion coefficient
${\beta }_2$	Mass expansion coefficient
$\psi$	Separation ratio
$\epsilon$	Porosity
$ϵ$	Expansion parameter
$\Lambda$	Applied temperature gradient
$\chi$	Porosity number
Superscript
${}^{\prime }$	Variables in dimensional form

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Variation of [Forumla omitted. See PDF.] with respect to [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.].

Figure 2. Solutions of Z-direction velocity component (w) profile against X-direction at [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.].

Figure 3. Solutions of Z-direction velocity component (w) profile against Z-direction at [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.].

Figure 4. Variation of [Forumla omitted. See PDF.] with respect to X with R = 3[Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different R (b) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for different [Forumla omitted. See PDF.].

Figure 5. Dependence of [Forumla omitted. See PDF.] on R for (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 6. Dependence of [Forumla omitted. See PDF.] on t with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and (d) Comparison of present work with numerical results: Dotted lines represent numerical values obtained by Brand and Steiberg [9] and solid curves represent analytical results of present work. (e) Measurement of the vertical component of the velocity (W) for different R with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 7. Energy (Kinetic energy (KE) and potential energy (PE)) with respect to [Forumla omitted. See PDF.]. (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 8. Streamlines (left), isotherms (middle) and isoconcentration lines (right) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.] (a–c), R = [Forumla omitted. See PDF.] (d–f), [Forumla omitted. See PDF.] (g–i), [Forumla omitted. See PDF.] (j–l).

Figure 9. Streamlines (left), isotherms (middle) and isoconcentration lines (right) for [Forumla omitted. See PDF.] (a–c) and [Forumla omitted. See PDF.] (d–f) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], R = [Forumla omitted. See PDF.].

Figure 10. Streamlines (left), isotherms (middle) and isoconcentration lines (right) for [Forumla omitted. See PDF.] (a–c) and [Forumla omitted. See PDF.] (d–f) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and R = [Forumla omitted. See PDF.].

Figure 11. Streamlines (left), isotherms (middle) and isoconcentration lines (right) for [Forumla omitted. See PDF.] (a–c) and 1.6 (d–f) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], R = [Forumla omitted. See PDF.].

Figure 12. Streamlines (left), isotherms (middle) and isoconcentration lines (right) for [Forumla omitted. See PDF.] (a–c) and [Forumla omitted. See PDF.] (d–f) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], R = [Forumla omitted. See PDF.].

Figure 13. Vorticity lines for (a) [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], (b) R = [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 14. Heatlines for [Forumla omitted. See PDF.] (a), R = [Forumla omitted. See PDF.] (b), [Forumla omitted. See PDF.] (c), and [Forumla omitted. See PDF.] (d) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 15. Heatlines for [Forumla omitted. See PDF.] (a) and [Forumla omitted. See PDF.] (b) at [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], R = [Forumla omitted. See PDF.].

Figure 16. Heatlines for [Forumla omitted. See PDF.] (a) and [Forumla omitted. See PDF.] (b) at [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and R = [Forumla omitted. See PDF.].

Figure 17. Heatlines at [Forumla omitted. See PDF.] (a) and [Forumla omitted. See PDF.] (b) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], R = [Forumla omitted. See PDF.].

Figure 18. Heatlines of [Forumla omitted. See PDF.] (a) and [Forumla omitted. See PDF.] (b) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], R = [Forumla omitted. See PDF.].

Figure 19. Local entropy generation due to heat transfer (SHTI) and fluid friction (SFFI) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (a,b) [Forumla omitted. See PDF.] (c,d) [Forumla omitted. See PDF.] and (e,f) [Forumla omitted. See PDF.].

Figure 19. Local entropy generation due to heat transfer (SHTI) and fluid friction (SFFI) for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], (a,b) [Forumla omitted. See PDF.] (c,d) [Forumla omitted. See PDF.] and (e,f) [Forumla omitted. See PDF.].

Table 1

Power law equations for different $\psi$.

$\mathit{\psi }$	${\mathit{R}}_{\mathbf{co}}$	$\mathit{Nu}=\mathit{a}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathit{co}}}\right)}^{\mathit{r}}$	$\mathit{Nu}=\mathit{a}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathit{co}}}\right)}^{\mathit{r}}+\mathit{c}$
$-0.01$	47.44498	$1.0025{\left(\frac{R}{R_{co}}\right)}^{0.2197}$	$-0.157{\left(\frac{R}{R_{co}}\right)}^{-1.845}+1.158$
$-0.02$	47.54581	$1.0073{\left(\frac{R}{R_{co}}\right)}^{0.2994}$	$-0.134{\left(\frac{R}{R_{co}}\right)}^{-3.772}+1.136$
$-0.03$	47.64702	$1.0121{\left(\frac{R}{R_{co}}\right)}^{0.3835}$	$-0.154{\left(\frac{R}{R_{co}}\right)}^{-4.88}+1.156$

Table 2

Power law equation for different $Le$.

$\mathit{Le}$	${\mathit{R}}_{\mathit{co}}$	$\mathit{Nu}=\mathit{a}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathit{co}}}\right)}^{\mathit{r}}$	$\mathit{Nu}=\mathit{a}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathit{co}}}\right)}^{\mathit{r}}+\mathit{c}$
$0.04$	47.64702	$1.01296{\left(\frac{R}{R_{co}}\right)}^{0.4277}$	$-0.154{\left(\frac{R}{R_{co}}\right)}^{-4.88}+1.156$
$0.05$	49.72228	$1.0126{\left(\frac{R}{R_{co}}\right)}^{0.3991}$	$-0.1969{\left(\frac{R}{R_{co}}\right)}^{-3.52}+1.2$
$0.06$	51.7553	$1.0121{\left(\frac{R}{R_{co}}\right)}^{0.3835}$	$-0.279{\left(\frac{R}{R_{co}}\right)}^{-2.421}+1.283$

Table 3

Power law equations for different $\chi$.

$\mathit{\chi }$	${\mathit{R}}_{\mathit{co}}$	$\mathit{Nu}=\mathit{a}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathit{co}}}\right)}^{\mathit{r}}$	$\mathit{Nu}=\mathit{a}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathit{co}}}\right)}^{\mathit{r}}+\mathit{c}$
0.2	47.64702	$1.0121{\left(\frac{R}{R_{co}}\right)}^{0.3835}$	$-0.1547{\left(\frac{R}{R_{co}}\right)}^{-4.88}+1.156$
0.3	50.61914	$1.0081{\left(\frac{R}{R_{co}}\right)}^{0.2582}$	$-0.1278{\left(\frac{R}{R_{co}}\right)}^{-3.622}+1.129$
0.4	53.57611	$1.0062{\left(\frac{R}{R_{co}}\right)}^{0.2071}$	$-0.11807{\left(\frac{R}{R_{co}}\right)}^{-2.905}+1.119$

References

1. Platten, J.K.; Chavepeyer, G. Nonlinear two dimentional Bénard convection with soret effect: Free boundaries. Int. J. Heat Mass Trans.; 1977; 20, pp. 113-122. [DOI: https://dx.doi.org/10.1016/0017-9310(77)90003-5]

2. Steinberg, V. Oscillatory convective instability in a superfluid mixture. Phys. Rev. A; 1981; 24, pp. 2584-2594. [DOI: https://dx.doi.org/10.1103/PhysRevA.24.2584]

3. Steinberg, V.; Brand, H.R. Two fluid effects in the convective instability of superfluid mixture. Phys. Rev. B; 1983; 28, pp. 1618-1620. [DOI: https://dx.doi.org/10.1103/PhysRevB.28.1618]

4. Paul Kolodner, A.; Passner,; Surko, C.M.; Walden, R.W. Onset of oscillatory convection in a binary fluid mixture. Phys. Rev. Lett.; 1986; 56, pp. 2621-2624. [DOI: https://dx.doi.org/10.1103/PhysRevLett.56.2621]

5. Guenter, A.; Rehberg, I. Convection in a binary mixture heated from below. Phys. Rev. Lett.; 1986; 56, pp. 1373-1376.

6. Barbara, J.A.; Zielinska,; Brand, H.R. Exact solution of the linear stability problem for the onset of convection in binary fluid mixtures. Phys. Rev. A; 1987; 35, pp. 4349-4353.

7. Brand, H.R.; Steinberg, V. Convective instabilities in binary mixture in a porous medium. Phys. A; 1983; 119, pp. 327-338. [DOI: https://dx.doi.org/10.1016/0378-4371(83)90164-4]

8. Steinberg, V.; Brand, H.R. Convective instabilities of binary mixtures with fast chemical reaction in a porous medium. J. Chem. Phys.; 1983; 78, pp. 2655-2660. [DOI: https://dx.doi.org/10.1063/1.445024]

9. Brand, H.R.; Steinberg, V. Nonlinear effects in the convective instability of a binary mixture in a porous medium near threshold. Phys. Lett. A.; 1983; 93, pp. 333-336. [DOI: https://dx.doi.org/10.1016/0375-9601(83)90110-X]

10. Steinberg, V.; Brand, H.R. Amplitude equations for the onset of convection in a reactive mixture in a porous medium. J. Chem. Phys.; 1984; 80, pp. 431-435. [DOI: https://dx.doi.org/10.1063/1.446466]

11. Brand, H.R.; Hohenberg, C. Codimension-2 bifurcations for convection in binary fluid mixtures. Phys. Rev. A; 1984; 30, pp. 2548-2561. [DOI: https://dx.doi.org/10.1103/PhysRevA.30.2548]

12. Rehberg, I.; Ahlers, G. Experimental observation of a codimension-two bifurcation in a binary fluid mixture. phys. Rev. Lett.; 1985; 55, pp. 500-503. [DOI: https://dx.doi.org/10.1103/PhysRevLett.55.500] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10032369]

13. Bahloul, A.; Boutana, N.; Vasseur, P. Double-diffusive and Soret-induced convection in a shallow horizontal porous layer. J. Fluid Mech.; 2003; 491, pp. 325-352. [DOI: https://dx.doi.org/10.1017/S0022112003005524]

14. Umla, R.; Augustin, M.; Huke, B.; Lucke, M. Roll convection of binary fluid mixtures in porous media. J. Fluid Mech.; 2010; 649, pp. 165-186. [DOI: https://dx.doi.org/10.1017/S0022112009993302]

15. Nield, A. Onset of thermohaline convection in a porous medium. Water Resour. Res.; 1968; 4, pp. 553-560. [DOI: https://dx.doi.org/10.1029/WR004i003p00553]

16. Trevisan, V.; Bejan, A. Natural convection with combined heat and mass transfer buoyancy effects in a porous medium. Int. J. Heat Mass Trans.; 1985; 28, pp. 1597-1611. [DOI: https://dx.doi.org/10.1016/0017-9310(85)90261-3]

17. Trevisan, V.; Bejan, A. Mass and heat transfer by high Rayleigh number convection in a porous medium heated from below. Int. J. Heat Mass Trans.; 1987; 30, pp. 2341-2356. [DOI: https://dx.doi.org/10.1016/0017-9310(87)90226-2]

18. Augustin, M.; Umla, R.; Huke, B.; Lcke, M. Stationary and oscillatory convection of binary fluids in a porous medium. Phys. Rev. E; 2010; 82, e056303-10. [DOI: https://dx.doi.org/10.1103/PhysRevE.82.056303]

19. Griffiths, R.W. Layered double-diffusive convection in porous media. J. Fluid Mech.; 1981; 102, pp. 221-248. [DOI: https://dx.doi.org/10.1017/S0022112081002619]

20. Hounsou, P.; Monwanou, A.V.; Miwadinou, C.H.; Chabi Orou, J.B. Stationary convection in a binary mixture of ferrofluids in a porous medium. Fluid Dyn. Res.; 2020; 52, 035501. [DOI: https://dx.doi.org/10.1088/1873-7005/ab8809]

21. Hu, Z.C.; Zhang, X.R. A new oscillatory instability in Rayleigh-Bénard convection of a binary mixture with positive separation ratio. Phys. Fluids; 2021; 33, 054113-12. [DOI: https://dx.doi.org/10.1063/5.0049247]

22. Rashmi, D.; Murthy, P.V.S.N. Nonlinear stability of Prats flow in Brinkman porous medium with Soret effect. Int. J. Thermal. Sci.; 2022; 179, pp. 107-486.

23. Kuo, H.L. Solution of the non-linear equations of the cellular convection and heat transport. J. Fluid Mech.; 1961; 10, pp. 611-630. [DOI: https://dx.doi.org/10.1017/S0022112061000408]

24. Palm, E.; Weber, J.E.; Kvernvold, O. On steady convection in a porous medium. J. Fluid Mech.; 1972; 54, pp. 153-161. [DOI: https://dx.doi.org/10.1017/S002211207200059X]

25. Rameshwar, Y.; Rawoof Sayeed, M.A.; Rani, H.P.; Laroze, D. Finite amplitude cellular convection under the influence of a vertical magnetic field. Int. J. Heat Mass Trans.; 2017; 114, pp. 559-577. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.06.076]

26. Rameshwar, Y.; Srinivas, G.; Laroze, D.; Rawoof Sayeed, M.A.; Rani, H.P. Convective instabilities in binary mixture ³He−⁴He in porous media. Chin. J. Phys.; 2022; 77, pp. 773-803. [DOI: https://dx.doi.org/10.1016/j.cjph.2022.03.010]

27. Rawoof Sayeed, M.A.; Rameshwar, Y. Finite amplitude cellular thermohaline convection. J. Heat Trans.; 2022; 144, 112602. [DOI: https://dx.doi.org/10.1115/1.4055346]

28. Homsy and Sherwood, Convective instabilities in porous media with through flow. AlChE; 1976; 22, pp. 168-174. [DOI: https://dx.doi.org/10.1002/aic.690220121]

29. Landau, L.; Lifshitz, E. Fluid Mechanics; Pergamon Press: Oxford, UK, 1959.

30. Chandrashekar, S. Hydrodynamic and Hydromagnetic Stability; Oxford University: New York, NY, USA, 1961.

31. Pellew, A.; Southwell, R.V. On maintained convective motion in a fluid heated from below. R. Soc. Lond. A; 1940; 176, pp. 312-343.

32. Bhardwaj, S.; Dalal, A.; Patil, S. Influence of wavy wall and non-uniform heating on natural convection heat transfer and entropy generation inside porous complex enclosure. Energy; 2015; 79, pp. 467-481. [DOI: https://dx.doi.org/10.1016/j.energy.2014.11.036]

33. Sparrow, E.M.; Carlson, C.K. Local and average natural convection Nusselt numbers for a uniformly heated, shrouded or unshrouded horizontal plate. Int. J. Heat Mass Trans.; 1986; 29, pp. 369-379. [DOI: https://dx.doi.org/10.1016/0017-9310(86)90207-3]

34. Batchelor, G.K. An Introduction to Fluid Dynamics; Cambridge University: Cambridge, UK, 1993.

35. Jana, S.C.; Metcalfe, G.; Ottino, J.M. Experimental andcomputational studies of mixing in complex Stokes flows: The vortex mixing flow and multicellular cavity flows. J. Fluid Mech.; 1994; 296, pp. 199-249. [DOI: https://dx.doi.org/10.1017/S0022112094001539]

36. Sheu, W.H.; Rani, H.P. Exploration of vortex dynamics for transitional flows in a three-dimensional backward facing step channel. J. Fluid Mech.; 2006; 550, pp. 61-83. [DOI: https://dx.doi.org/10.1017/S0022112005007858]

37. Sheu, W.H.; Rani, H.P.; Tan, T.C.; Tsai, S.F. Multilple states, topology and bifurcations of natural convection in a cubical cavity. Comp. Fluids; 2008; 37, pp. 1011-1028. [DOI: https://dx.doi.org/10.1016/j.compfluid.2007.11.003]

38. Morega, A.I.; Bejan, A. Heatline visualization of forced convection laminar boundry layers. Int. J. Heat Mass Trans.; 1993; 36, pp. 3957-3966. [DOI: https://dx.doi.org/10.1016/0017-9310(93)90146-W]

39. Aggarwal, S.K.; Mahapra, A. Use of heatlines for unsteady buoyancy-driven flow in cylinderical enclosure. ASME J. Heat Trans.; 1989; 111, pp. 576-578. [DOI: https://dx.doi.org/10.1115/1.3250719]

40. Nield, D.A.; Bejan, A. Convection in Porous Media; 5th ed. Springer: Cham, Switzerland, 1999.

41. Costa, V.A.F. Bejan’s heatlines and masslines for convection visualization and analysis. Appl. Mech. Rev.; 2006; 59, pp. 126-145. [DOI: https://dx.doi.org/10.1115/1.2177684]

42. Kimura, S.; Bejan, A. The heatline visualization of convective heat transfer. J. Heat Trans.; 1983; 105, pp. 916-919. [DOI: https://dx.doi.org/10.1115/1.3245684]

43. Basak, T.; Roy, S. Role of Bejan’s heatlines in heat flow visualization and optimal thermal mixing for differentially heated square enclosures. Int. J. Heat Mass Trans.; 2008; 51, pp. 3486-3503. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2007.10.033]

44. Basak, T.; Roy, S.; Matta, A.; Pop, I. Analysis of heatlines for natural convection within porous trapezoidal cnclosure: Effect of uniform and non-uniform heating of bottom wall. Int. J. Heat Mass Trans.; 2010; 53, pp. 5947-5961. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2010.07.026]

45. Varol, Y.; Oztop, H.F.; Mobedi, M.; Pop, I. Visualization of heat flow using Bejan’s heatline due to natural convection of water near 4 °C in thick walled porous cavity. Int. J. Heat Mass Trans.; 2010; 53, pp. 1691-1698. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2010.01.020]

46. Hooman, K.; Gurgenci, H.; Dincer, I. Heatline and energy-flux-vector visualization of natural convection in a porous cavity occupied by a fluid with temperature-dependent viscosity. J. Porous Media; 2009; 12, pp. 265-275. [DOI: https://dx.doi.org/10.1615/JPorMedia.v12.i3.60]

47. Bejan, A. Entropy Generation through Heat and Fluid Flow; Wiley: New York, NY, USA, 1982.

48. Basak, T.; Anandalakshmi, R.; Roy, S.; Pop, I. Role of entropy generation on thermal management due to thermal convection in porous trapezoidal enclosures with isothermal and non-isothermal heating of wall. Int. J. Heat Mass Trans.; 2013; 67, pp. 810-828. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.08.019]

49. Basak, T.; Das, D.; Biswal, P. Heatlines: Modeling, visualization, mixing and thermal management. Prog. Energy Combustion Sci.; 2017; 64, pp. 157-218. [DOI: https://dx.doi.org/10.1016/j.pecs.2017.08.003]

50. Al-Hadhrami, A.K.; Elliott, L.; Ingham, D.B. A New Model for Viscous Dissipation in Porous Media Across a Range of Permeability Values. Trans. Porous Med.; 2003; 53, pp. 117-122. [DOI: https://dx.doi.org/10.1023/A:1023557332542]

51. Homan, K.; Gurgenci, H. Effects of viscous dissipation and boundary conditions on forced convection in a channel occupied by a saturated porous medium. Trans. Porous Med.; 2007; 68, pp. 301-319. [DOI: https://dx.doi.org/10.1007/s11242-006-9049-4]

52. Baytas, A.C. Entropy generation for natural convection in an inclined porous cavity. Int. J. Heat Mass Trans.; 2000; 43, pp. 2089-2099. [DOI: https://dx.doi.org/10.1016/S0017-9310(99)00291-4]