Weakly Nonlinear Stability Analysis of a Nanofluid in a Horizontal Porous Layer Using a Multidomain Spectral Collocation Method

Abstract

In this chapter, we present a weakly nonlinear stability analysis of the flow of a nanofluid in a porous medium with stress-free boundary conditions. Some previous studies have investigated cross-diffusion in a nanofluid layer although in most cases these studies mostly deal with linear stability analysis. It is important to study the nonlinear stability in flows subject to cross-diffusion due to the wide range of applications where such flows arise such as in hydrothermal growth, compact heat exchanges, the solidification of binary mixtures, geophysical systems, solar pond, etc. Here we consider flow between parallel plates with an applied magnetic field and zero nanoparticle flux at the boundaries. A truncated Fourier series is introduced reducing the flow equations to a Lorenz-type system of nonlinear evolution equations. The multidomain spectral method is used to solve the equations that describe the growth of the convection amplitudes. The solutions are obtained as sets of trajectories in the phase space. Some interesting spiral trajectories and their sensitivity to the Rayleigh number are given.

1. Introduction

The enhancement of thermal conductivity of a fluid is a matter of supreme interest to engineers due to the important applications of fluids in heat transfer processes. Natural and forced convection plays an important role in heat transfer processes due to continuous molecular movements in fluid. Recent studies show that the suspension of solid nanoparticles in a fluid can substantially improve the fluid’s thermophysical properties, including thermal conductivity.

The term nanofluid describes a liquid containing a suspension of nanometer sized 1–100 nm solid particles [1]. Examples of commonly used nanoparticles include metallic particles such as Al, Cu and Ag, and oxides such as Al₂O₃ and CuO. The base fluid is often a common liquid such as water, ethylene, glycol, or oil. The enhancement of thermochemical properties of a fluid due to the addition of nanoparticles has been observed in experimental studies such as in [2, 3]. Researchers have investigated the influence of seven slip mechanisms, namely, inertia, Brownian diffusion, thermophoresis, diffusiophoresis, magnus effect, fluid drainage, and gravity in nanofluids. It has been shown that, in the absence of turbulence, the most significant among these mechanisms are the Brownian diffusion and thermophoresis.

The classical Rayleigh-Benard convection problem in a heated horizontal layer has been extensively studied in the literature. Among recent studies on nanofluids, Tzou [4] studied the thermal instability and natural convection in nanofluid flow using an eigenfunction expansion method. Narayana et al. [5, 6] studied convection and the stability of a Maxwell fluid in a porous medium. Yadav et al. [7] investigated thermal instability of a rotating nanofluid layer. The studies by Kuznetsov and Nield [8, 9, 10, 11] focused on thermal instability in a porous layer saturated with a nanofluid. They investigated the onset of instability in a horizontal porous layer using a model for the nanofluid that incorporated particle Brownian motion and thermophoresis. Related studies with various assumptions on the geometry and flow structure have been made by [12, 13, 14, 15]. In the last few decades, researchers have also investigated thermal instability in a horizontal nanofluid layer subject to an applied magnetic field [16, 17]. The effects of a magnetic field on convection and the onset of instability have important applications in problems such as in cooling systems, pumps, magnetohydrodynamics and generators. The experimental study by Heris et al. [18] showed that thermal efficiency could be achieved by subjecting the flow to a magnetic field. The studies by Ghasemi et al. [19] and Hamad et al. [20] focused on the flow behavior and heat transfer in an electrically conducting nanofluid under the influence of a magnetic field and subject to Brownian diffusion and thermophoresis. They used a water-based nanofluid containing different types of nanoparticles such as copper, alumina and silver in their numerical simulations. Related studies of interest include [21, 22, 23, 24]. Rana et al. [25] studied thermal convection in a Walters (Model B) fluid in a porous medium. They showed that a magnetic field may introduce oscillatory instability modes and acts to stabilize the system.

In this chapter, we give a weakly nonlinear stability analysis of a nanofluid layer with an applied magnetic field, stress free boundary conditions and under the assumption of zero nanoparticle flux at the boundary. The studies by Kuznetsov and Nield [9] and Nield and Kuznetsov [10, 11] investigated cross-diffusion in a nanofluid layer. However, these studies mostly presented a linear stability analysis. It is important to study the nonlinear regime for a nanofluid flow subject to cross-diffusion due to the wide range of applications where such flows may arise. Typical examples may be found in hydrothermal growth, compact heat exchanges, solidification of binary mixtures, geophysical systems, and so on. Hence, with this in mind, we studied the finite amplitude convection in a nanofluid flows subject to cross-diffusion. By introducing a truncated Fourier series, a Lorenz-type system of seven nonlinear differential equations is obtained. The recent multidomain spectral method is used to solve the nonlinear equations. This method is accurate and very easy to implement compared to older methods such as finite difference methods. An analysis of heat and mass transfer for different parameters such as the Prandtl number, the Dufour and thermophoresis is presented.

2. Mathematical formulation

Consider viscous incompressible MHD nanofluid flow in an infinitely extended horizontal porous layer, confined between two boundaries at $z=0$ and $z=h$ . The layer is heated from below and cooled from above, see Figure 1. A Cartesian frame of reference is chosen in which the z-axis is vertically upward. The boundaries are perfectly conducting. The temperature at the lower and upper walls is $T_c$ and $T_h$ , respectively with $T_h>T_c$ . The Oberbeck-Boussinesq approximation and the Darcy law are assumed to be applicable. The continuity equation, momentum equation, energy equation, concentration equation and volumetric fraction nanoparticle equation, which describe the above configuration in dimensionless form, are given as

$$\nabla \cdot V=0,$$ E1

$$\frac{\mathit{Da}}{\mathit{Pr}}\frac{\partial V}{\partial t}=-\nabla P+\mathit{Da}{\nabla }^{2}V-V+\mathit{QV}{\widehat{e}}_z-\mathit{Rm}{\widehat{e}}_z+\mathit{RaT}{\widehat{e}}_z+\mathit{RsC}{\widehat{e}}_z-\mathit{Rn\varphi }{\widehat{e}}_z,$$ E2

$$\frac{\partial T}{\partial t}+V\cdot \nabla T={\nabla }^{2}T+\frac{N_B}{\mathit{Les}}\nabla \varphi \cdot \nabla T+\frac{N_A{N}_B}{\mathit{Les}}\nabla T\cdot \nabla T+\mathit{Du}{\nabla }^{2}C,$$ E3

$$\frac{\partial C}{\partial t}+V\cdot \nabla C=\frac{1}{\mathit{Le}}{\nabla }^{2}C+\mathit{Sr}{\nabla }^{2}T,$$ E4

$$\frac{1}{\sigma }\frac{\partial \varphi }{\partial t}+\frac{1}{\epsilon }V\cdot \nabla \varphi =\frac{1}{\mathit{Les}}{\nabla }^2\varphi +\frac{N_A}{\mathit{Les}}{\nabla }^{2}T,$$ E5

subject to the boundary conditions

$$V=0,T=1,C=1\frac{\partial \varphi }{\partial z}+N_A\frac{\partial T}{\partial z}=0\mathrm{at}z=0,$$ E6

$$V=0,T=0,C=0\frac{\partial \varphi }{\partial z}+N_A\frac{\partial T}{\partial z}=0\mathrm{at}z=1,$$ E7

where $V$ is the fluid velocity, $T$ is the temperature, $C$ is the solute concentration and $\varphi$ is the volumetric fraction of nanoparticles. The dimensionless parameters are the Darcy number (modified by the viscosity ratio) $\mathit{Da}$ , Prandtl number $\mathit{Pr}$ , Hartmann-Darcy number $Q$ , thermal Rayleigh-Darcy number $\mathit{Ra}$ , nanoparticle Rayleigh number $\mathit{Rn}$ and the basic density Rayleigh number $\mathit{Rm}$ . The parameter $N_A$ is a modified diffusivity ratio, $\mathit{Le}$ is the Lewis number, $\mathit{Rs}$ is solutal Rayleigh number, $N_B$ is a modified nanoparticle density increment and $\mathit{Du}$ is a modified Dufour parameter. The parameter $\mathit{Les}$ is the thermo-nanofluid Lewis number, $\nu$ is the kinematic viscosity and $\mathit{Sr}$ is a modified Soret parameter. These parameters have the form

$$\mathit{Da}=\frac{\tilde{\mu }K}{\mu h^2},\mathit{Pr}=\frac{\mu }{{\rho }_f{\alpha }_m},Q=\frac{\delta {\mathbf{B}}_0^{2}K}{\mu },\mathit{Ra}=\frac{{\rho }_f\mathit{\text{βKhg}}\left(T_h^{\ast }-T_c^{\ast }\right)}{{\mathit{\mu \alpha }}_m},\mathit{Les}=\frac{{\alpha }_m}{D_B},$$ E8

$$\mathit{Rn}=\frac{\left({\rho }_p-{\rho }_f\right)\left({\varphi }_1^{\ast }-{\varphi }_0^{\ast }\right)\mathit{gKh}}{{\mathit{\mu \alpha }}_m},\mathit{Rm}=\frac{{\rho }_p{\varphi }_0^{\ast }+\left(1-{\varphi }_0^{\ast }\right){\rho }_f\mathit{gKh}}{{\mathit{\mu \alpha }}_m},\mathit{Le}=\frac{{\alpha }_m}{D_S},$$ E9

$$N_B=\frac{\epsilon {\left(\mathit{\rho c}\right)}_p}{{\left(\mathit{\rho c}\right)}_f}\left({\varphi }_1^{\ast }-{\varphi }_0^{\ast }\right),N_A=\frac{D_T\left(T_h^{\ast }-T_c^{\ast }\right)}{D_B{T}_c^{\ast }\left({\varphi }_1^{\ast }-{\varphi }_0^{\ast }\right)},\mathit{Rs}=\frac{{\rho }_f\mathit{\text{βKhg}}\left(C_h^{\ast }-C_c^{\ast }\right)}{{\mathit{\mu \alpha }}_m},$$ E10

$$\mathit{Du}=\frac{\sigma D_{\mathit{TC}}\left(C_h^{\ast }-C_c^{\ast }\right)}{{\alpha }_m\left(T_h^{\ast }-T_c^{\ast }\right)},\mathit{Sr}=\frac{\sigma D_{\mathit{CT}}\left(T_h^{\ast }-T_c^{\ast }\right)}{{\alpha }_m\left(C_h^{\ast }-C_c^{\ast }\right)},$$ E11

where ${\rho }_f,{\rho }_p,\tilde{\mu },{\beta }_1,{\beta }_2,{\kappa }_m,\delta$ , $\epsilon$ and $K$ are the fluid density, nanoparticle density, effective viscosity of porous medium, thermal volumetric expansion coefficient of the fluid, solutal volumetric expansion coefficient, the thermal conductivity of porous medium, the electrical conductivity, the porosity, and permeability of porous medium, respectively. The gravitational acceleration is denoted by $g$ and $D_B$ is the Brownian diffusion coefficient, $D_T$ is the thermophoresis diffusion coefficient, $D_S$ is the solutal diffusion coefficient, $D_{\mathit{TC}}$ is the Dufour parameter and $D_{\mathit{CT}}$ is the Soret parameter. The heat capacity of the fluid is ${\left(\mathit{\rho c}\right)}_f$ , ${\left(\mathit{\rho c}\right)}_p$ is the effective heat capacity of the nanoparticle, ${\left(\mathit{\rho c}\right)}_m$ is the effective heat capacity of the porous medium and ${\mathbf{B}}_0$ is the uniform magnetic field strength.

Figure 1.

A schematic diagram of the problem.

The basic state is the time independent solution of Eqs. (1)–(5). Solving these equations with boundary conditions, we obtain

$$T_b=1-z,C_b=1-z,{\varphi }_b={\varphi }_0+N_{A}z.$$ E12

3. Weakly nonlinear stability analysis

In this section, we restrict the analysis to the case of two-dimensional disturbances. We define the stream function $\Psi$ by the equations

$$u=\frac{\partial \Psi }{\partial z},w=-\frac{\partial \Psi }{\partial x}.$$

Eqs. (1)–(5) may now be simplified by introducing the truncated Fourier series

$${\Psi }^{\prime }=A_{11}\text{sin}\mathit{\alpha x}\text{sin}\mathit{\pi z},T^{\prime }=B_{11}\text{cos}\mathit{\alpha x}\text{sin}\mathit{\pi z}+B_{02}\text{sin}2\mathit{\pi z},$$ E13

$$C^{\prime }=C_{11}\text{cos}\mathit{\alpha x}\text{sin}\mathit{\pi z}+C_{02}\text{sin}2\mathit{\pi z},{\varphi }^{\prime }=-N_A\left(D_{11}\text{cos}\mathit{\alpha x}\text{sin}\mathit{\pi z}+D_{02}\text{sin}2\mathit{\pi z}\right),$$ E14

where $A_{11},B_{11},B_{02},C_{11},C_{02},D_{11}$ and $D_{02}$ are amplitudes that depend on time. This leads to the Lorenz-type system of nonlinear ordinary differential equations

$${\dot{Y}}_1=\frac{\mathit{Pr}}{\mathit{Da}}\left\{-{\mathit{BY}}_1-N\left(Y_2+Y_4\right)+\frac{N_A\mathit{Rn}}{R}{Y}_6\right\}$$ E15

$${\dot{Y}}_2={\mathit{RY}}_1-Y_2-{\mathit{DuY}}_4-Y_1{Y}_3$$ E16

$${\dot{Y}}_3=\frac{1}{2}{Y}_1{Y}_2-G\left(Y_3+{\mathit{DuY}}_5\right)$$ E17

$${\dot{Y}}_4={\mathit{RY}}_1-\frac{1}{\mathit{Les}}{Y}_4-{\mathit{SrY}}_2-Y_1{Y}_5$$ E18

$${\dot{Y}}_5=\frac{1}{2}{Y}_1{Y}_2-G\left(\frac{1}{\mathit{Les}}{Y}_5+{\mathit{SrY}}_3\right)$$ E19

$${\dot{Y}}_6=\frac{\mathit{\sigma R}}{\epsilon }{Y}_1-\frac{\sigma }{\mathit{Le}}\left(Y_6-Y_2\right)-\frac{N_A\sigma }{\epsilon }{Y}_1{Y}_7$$ E20

$${\dot{Y}}_7=\frac{N_A\sigma }{2\epsilon }{Y}_1{Y}_7-\frac{G}{\mathit{Le}}\left(Y_7-Y_3\right)$$ E21

subject to $Y_n\left(0\right)=Y_n^0$ for $n=1,2,\cdots ,7.$ The following variables have been introduced in the equations above:

$$\begin{array}{c}{Y}_1=\frac{\mathit{\alpha \pi }}{\gamma }{A}_{11},Y_2=-\pi {\mathit{RB}}_{11},Y_3=-\pi {\mathit{RB}}_{02},Y_4=-\pi {\mathit{RC}}_{11},Y_5=-\pi {\mathit{RC}}_{20},\\ Y_6=-\pi {\mathit{RD}}_{11},Y_7=-\pi {\mathit{RD}}_{20},t^{\ast }=\mathit{\gamma t},R=\frac{{\alpha }^2}{{\gamma }^3}\mathit{Ra},G=\frac{4{\pi }^2}{\gamma }\text{and}N=\frac{\mathit{Rs}}{\mathit{Ra}},\end{array}$$

$$B=\frac{\mathit{Da}{\gamma }^2+\gamma -{\alpha }^{2}Q}{{\gamma }^2}.$$

Eqs. (15)–(21) give an approximate description of the full dimensional nonlinear system. An analytical solution of the system of nonlinear ordinary differential Eqs. (15)–(21) is not possible for the general time variable $t$ . However, it is possible to discuss the stability of the nonlinear system of equations. The system of equations is uniformly bounded in time and dissipative in the phase space. We can easily show that

$$\sum _{i=1}^7\frac{\partial {\dot{Y}}_i}{\partial Y_i}=-\left[\frac{\mathit{DaB}}{\mathit{Pr}}+1+G+{\mathit{Les}}^{-1}+{\mathit{\text{GLes}}}^{-1}+\sigma {\mathit{Le}}^{-1}+{\mathit{GLe}}^{-1}\right].$$ E22

This is always true if $B\ge 0$ . As has been shown in previous studies, the trajectories may be attracted to a fixed point, limit cycle or other attractor. For a set of initial points in the phase space occupying a region $V\left(0\right)$ at time $t=0$ , after a time $t>0$ , the end point of the corresponding trajectories fills a volume

$$V\left(t\right)=V\left(0\right)\text{exp}\left\{-\left[\frac{\mathit{DaB}}{\mathit{Pr}}+1+G+{\mathit{Les}}^{-1}+{\mathit{\text{GLes}}}^{-1}+\sigma {\mathit{Le}}^{-1}+{\mathit{GLe}}^{-1}\right]t\right\}.$$ E23

Eq. (23) shows that the volume decays exponentially with time. Further, it can be noted that the system of Eqs. (15)–(21) are invariant under the transformation

$$S\left(Y_1,Y_2,Y_3,Y_4,Y_5,Y_6,Y_7\right)\to -S\left(Y_1,Y_2,Y_3,Y_4,Y_5,Y_6,Y_7\right).$$ E24

We obtain the possible stationary points of the nonlinear system of equations by setting ${\dot{Y}}_i=0$ for $i=1,2,\cdots ,7$ . One of these stationary points is $Y_i=0$ and by linearizing about this point, we obtain the Jacobian matrix

$$A=\left(\begin{array}{lllllll}-\frac{\mathit{PrB}}{\mathit{Da}}& -\frac{\mathit{Pr}}{\mathit{Da}}& 0& \frac{\mathit{PrN}}{\mathit{Da}}& 0& \frac{{\mathit{PrN}}_A\mathit{Rn}}{\mathit{DaR}}& 0\\ R& -1& 0& -\mathit{Du}& 0& 0& 0\\ 0& 0& -G& 0& -\mathit{GDu}& 0& 0\\ R& -\mathit{Sr}& 0& -{\mathit{Les}}^{-1}& 0& 0& 0\\ 0& 0& -\mathit{GSr}& 0& -{\mathit{\text{GLes}}}^{-1}& 0& 0\\ \frac{\mathit{\sigma R}}{\epsilon }& \sigma {\mathit{Le}}^{-1}& 0& 0& 0& -\sigma {\mathit{Le}}^{-1}& 0\\ 0& 0& {\mathit{GLe}}^{-1}& 0& 0& 0& -{\mathit{GLe}}^{-1}\end{array}\right).$$ E25

The eigenvalues of the above matrix depend on the various parameters. For the specific parameters $R={10}^3,\mathit{Da}=20,\mathit{Pr}=10,N=25,\mathit{Du}=0.2,\mathit{Sr}=3,\mathit{Les}=10,\mathit{Le}=5,\sigma =0.05,G=3 \text{and}\epsilon =0.04$ , the characteristic polynomial is

$$P\left(\lambda \right)={\lambda }^7+21.41{\lambda }^6+651.5{\lambda }^5+5391.7{\lambda }^4+12772.232{\lambda }^3-370.962{\lambda }^2-2996.712\lambda +545.8$$

with eigenvalues

$$\begin{array}{l}{\lambda }_1=0.2955056985, {\lambda }_2=0.2402382976, {\lambda }_3=-0.6139990637, {\lambda }_4=-4.886000936,\\ {\lambda }_5=-5, {\lambda }_6=-5.7228719981-21.9033954659i, {\lambda }_7=-5.7228719981+21.9033954659i.\end{array}$$

This stationary point is a saddle point. Nonetheless, because the eigenvalues depend on various parameters, we cannot make general conclusions as to the stability of the system. We note, however, that if we denote the trace of the matrix $A$ by $T$ and the determinant $d$ , then

$$T=-\left(\frac{\mathit{DaB}}{\mathit{Pr}}+1+G+{\mathit{Les}}^{-1}+{\mathit{\text{GLes}}}^{-1}+\sigma {\mathit{Le}}^{-1}+{\mathit{GLe}}^{-1}\right),$$ E26

and

$$d=\frac{\sigma \text{Pr}{G}^3}{\mathit{\text{DaLeLes}}}\left(-\frac{\mathit{\text{DuNRSr}}}{\mathit{Le}}-\frac{\mathit{\text{DuNRSr}}}{\epsilon }+\frac{\mathit{\text{BDuSr}}}{\mathit{Le}}+\frac{\mathit{\text{DuRSr}}}{\mathit{Le}}\right.\left.+\frac{\mathit{NR}}{\mathit{\text{LeLes}}}+\frac{\mathit{NR}}{\mathit{\text{Lesε}}}-\frac{B}{\mathit{\text{LeLes}}}-\frac{R}{\mathit{\text{LeLes}}}\right).$$ E27

The trace is always negative, but the sign of determinant depends on the parameter values. If $d<0$ then

$$\left(1-N\right)\mathit{\text{εDuSrRLes}}+\left(\mathit{\epsilon B}-\mathit{\text{NRLe}}\right)\mathit{\text{DuSrLes}}+\left(\epsilon +\mathit{Le}\right)\mathit{NLe}<\left(B+R\right)\epsilon ,$$ E28

suggesting a saddle point.

4. Method of solution

To study the influence of various physical parameters on the average Nusselt and Sherwood numbers, we solved the nonlinear system of Eqs. (15)–(21) numerically using the multidomain spectral collocation method. This is a novel technique for solving nonlinear initial value problems and parabolic equations with large time domains. It has been suggested in the literature that the method gives better accuracy compared to other methods such as finite difference and Runge-Kutta methods [26]. To apply the multidomain spectral collocation to the nonlinear system of equations, we first divide the interval $\left[0,T\right]$ into subintervals ${\mathrm{\Omega }}_i=\left[t_{i-1},t_i\right]$ for $i=1,2,\cdots ,p$ . The transformation

$$t=\frac{t_i-t_{i-1}}{2}\tau +\frac{t_i+t_{i-1}}{2}$$ E29

is used to transform each subinterval ${\mathrm{\Omega }}_i$ into the interval $\left[-1,1\right]$ . The system of Eqs. (15)–(21) can be written in the form

$$\frac{{\mathit{dY}}_1^i}{\mathit{dt}}=\frac{\mathit{Pr}}{\mathit{Da}}\left\{-{\mathit{BY}}_1^i-N\left(Y_2^i+Y_4\right)+\frac{N_A\mathit{Rn}}{\mathit{Ra}}{Y}_6^i\right\},$$ E30

$$\frac{{\mathit{dY}}_2^i}{\mathit{dt}}={\mathit{RY}}_1^i-Y_2^i-{\mathit{DuY}}_4^i-Y_1^i{Y}_3^i,$$ E31

$$\frac{{\mathit{dY}}_3^i}{\mathit{dt}}=\frac{1}{2}{Y}_1^i{Y}_2^i-G\left(Y_3^i+{\mathit{DuY}}_5^i\right),$$ E32

$$\frac{{\mathit{dY}}_4^i}{\mathit{dt}}={\mathit{RY}}_1^i-\frac{1}{\mathit{Les}}{Y}_4^i-{\mathit{SrY}}_2^i-Y_1^i{Y}_5^i,$$ E33

$$\frac{{\mathit{dY}}_5^i}{\mathit{dt}}=\frac{1}{2}{Y}_1^i{Y}_2^i-G\left(\frac{1}{\mathit{Les}}{Y}_5^i+{\mathit{SrY}}_3^i\right),$$ E34

$$\frac{{\mathit{dY}}_6^i}{\mathit{dt}}=\frac{\mathit{\sigma R}}{\epsilon }{Y}_1^i-\frac{\sigma }{\mathit{Le}}\left(Y_6^i-Y_2^i\right)-\frac{N_A\sigma }{\epsilon }{Y}_1^i{Y}_7^i,$$ E35

$$\frac{{\mathit{dY}}_7^i}{\mathit{dt}}=\frac{N_A\sigma }{2\epsilon }{Y}_1^i{Y}_7^i-\frac{G}{\mathit{Le}}\left(Y_7^i-Y_3^i\right),$$ E36

subject to

$$Y_n^i\left(t_{i-1}\right)=Y_n^{i-1}\left(t_{i-1}\right)\text{for}n=1,2,\cdots ,7.$$ E37

The first step in using the multidomain spectral collocation method (MDSCM) concerns the quasilinearization of Eqs. (30)–(36) leading to a system of equations in the form

$$\sum _{n=1}^7{a}_{\left(j,n\right)r}^i{Y}_{n,r+1}^i-\frac{{\mathit{dY}}_{j,r+1}^i}{\mathit{dt}}=R_{\mathit{jr}}^i,$$ E38

subject to

$$Y_{n,r+1}^i\left(t_{i-1}\right)=Y_{n,r+1}^{i-1}\left(t_{i-1}\right)\text{for}n=1,2,\cdots ,7.$$ E39

where $a_{\left(j,n\right)r}^i$ and $R_{\mathit{jr}}^i$ for $j=1,2,\dots ,7$ are given in the Appendix. Having linearized the equations, the second step is to integrate Eqs. (30)–(36). To this end, we use the Gauss-Lobatto nodes

$${\tau }_j^i=\text{cos}\frac{\mathit{\pi j}}{\mathit{Nc}},\text{for}j=0,1,\cdots ,\mathit{Nc}.$$ E40

We approximate the derivatives of the unknown functions $Y_{n,r+1}^i\left(t\right)$ at the collocation points by

$$\frac{{\mathit{dY}}_{n,r+1}^i}{\mathit{dt}}\left({\tau }_j^i\right)={\sum }_{k=0}^{\mathit{Nc}}{\mathbf{D}}_{\mathit{jk}}{Y}_{n,r+1}^i\left({\tau }_j^j\right)={\left[{\mathbf{DU}}_{n,r+1}^i\right]}_j,$$ E41

where $\mathbf{D}=2D/\left(t_i-t_{i-1}\right)$ , $D$ is the Chebyshev differentiation matrix and

$${\mathbf{U}}_{n,r+1}^i={\left(Y_{n,r+1}^i\left({\tau }_0^i\right),\cdots Y_{n,r+1}^i\left({\tau }_{\mathit{Nc}}^i\right)\right)}^T,$$

is a vector of the unknown functions at the collocation points. Substituting Eq. (41) into Eqs. (38) and reducing the result into matrix form, we obtain

$$\begin{array}{lll}{\mathbf{AU}}_{n,r+1}^i& =& {\mathbf{R}}_n^i,\\ {\mathbf{U}}_{n,r+1}^i\left({\tau }_{N_c}^{i-1}\right)& =& {\mathbf{U}}_n^i\left({\tau }_{N_c}^{i-1}\right),n=1,2,\cdots ,7.\end{array}$$ E42

where the matrices $\mathbf{A}=\left[A_{\mathit{ij}}\right]$ and ${\mathbf{R}}_n^i$ are given in the Appendix.

5. Heat and mass transfer

The study of heat and mass transfer in a horizontal nanofluid layer heated from below and cooled from above has important engineering applications. We define the rate of heat transfer by the average Nusselt number $\mathit{Nu}\left(t\right)$ where

$$\mathit{Nu}\left(t\right)=1+{\left[\frac{\frac{\alpha }{2\pi }{\int }_0^{\frac{2\pi }{\alpha }}\frac{\partial T}{\partial z}\mathit{dx}}{\frac{\alpha }{2\pi }{\int }_0^{\frac{2\pi }{\alpha }}\frac{\partial T_b}{\partial z}\mathit{dx}}\right]}_{z=0}+\mathit{Du}\left\{1+{\left[\frac{\frac{\alpha }{2\pi }{\int }_0^{\frac{2\pi }{\alpha }}\frac{\partial C}{\partial z}\mathit{dx}}{\frac{\alpha }{2\pi }{\int }_0^{\frac{2\pi }{\alpha }}\frac{\partial C_b}{\partial z}\mathit{dx}}\right]}_{z=0}\right\}.$$ E43

Substituting Eqs. (12) and (13) into Eq. (43), we obtain

$$\mathit{Nu}\left(t\right)=1+\frac{2}{R}{Y}_3+\mathit{Du}\left(1+\frac{2}{R}{Y}_5\right).$$ E44

Similarly, the rate of mass transfer stated in terms of the average Sherwood number is

$$\mathit{Sh}\left(t\right)=1+\frac{2}{R}{Y}_5+\mathit{Sr}\left(1+\frac{2}{R}{Y}_3\right)$$ E45

6. Results and discussion

We have studied the weakly nonlinear instability of nanofluid flow in a horizontal layer with stress free boundary conditions. For numerical simulations, the parameter values were chosen from the literature on nanofluid flow such as [4, 7]. In the literature, the critical Rayleigh number is found when the Darcy number is very large. In this study, we investigated the critical Rayleigh number for low Darcy numbers.

The method of solution described in Section 4 was used to solve Eqs. (15)–(21). All computations are carried out up to a value of maximum time $t_{\mathit{max}}=1$ , and solutions are obtained using initial conditions selected in the neighborhood of stationary points. Periodic solution sets were obtained for the system of nonlinear equations. We determined the rate of heat and mass transfer as functions of time for different parameter values. The results are shown in Figures 2–4. Figure 2 shows the effect of the Dufour and Soret parameters on the Nusselt and Sherwood numbers with time $t$ . Figure 2(a) shows how the heat transfer coefficient changes with both the Dufour parameter and time. The heat transfer coefficient increases with the Dufour parameter but eventually settles to a steady value with time. In Figure 2(b), the Soret parameter is similarly shown to enhance the mass transfer coefficient. We investigated the effect of the Prandtl and Lewis numbers (see Figures 3 and 4). An increase in the Lewis number enhances both heat and mass transfer in a nanofluid layer heated from below. However, Figure 3 shows that increasing the Prandtl number reduces the amplitude of oscillatory heat and mass transfer. The Prandtl number can lead to both positive and negative contributions to the Nusselt and Sherwood numbers. It is interesting to note that our investigation shows that the magnetic field parameter has very little effect on the heat and mass transfer for this type of flow.

Figure 2.

The effect of cross-diffusive parameters on (a) the Nusselt number $\mathit{Nu}$ and (b) the Sherwood number $\mathit{Sh}$ for $\mathit{Da}=0.05$ , $\mathit{Le}=2$ , $\mathit{Du}=0.2$ , $\epsilon =0.04$ , $\sigma =0.05$ , $\mathit{Les}=100$ , $\mathit{Rn}=5$ , $\mathit{Ra}=1000$ and various values of the Dufour and Soret parameters.

Figure 3.

The effect of Prandtl number $\mathit{Pr}$ on (a) the Nusselt number $\mathit{Nu}$ and (b) the Sherwood number $\mathit{Sh}$ when $\mathit{Da}=0.05$ , $\mathit{Le}=2$ , $\mathit{Du}=0.2$ , $\epsilon =0.04$ , $\sigma =0.05$ , $\mathit{Les}=100$ , $\mathit{Rn}=5$ and $\mathit{Ra}=1000$ .

Figure 4.

The effect of Lewis number on (a) the Nusselt number $\mathit{Nu}$ and (b) the Sherwood number $\mathit{Sh}$ when $\mathit{Da}=0.05$ , $\mathit{Du}=0.2$ , $\epsilon =0.04$ , $\sigma =0.05$ , $\mathit{Les}=100$ , $\mathit{Rn}=5$ and $\mathit{Ra}=1000$ .

Figures 5–11 show the effect of the Rayleigh number on the trajectories projected onto the $\left(Y_i,Y_j\right)$ phase planes. The solution sets provide a visual representation of the system’s behavior with every phase point on the phase space representing the physical state of the system. The convective solution sets for different values of $R$ have been presented with the trajectories projected onto the $\left(Y_i,Y_j\right)$ phase planes. These trajectories spiral toward the fixed point for Rayleigh numbers from ${10}^2$ to ${10}^4$ . The solution sets give spiral phase portraits as $R$ increases and for the high Rayleigh numbers, the trajectories spiral many times before they reach a fixed point.

Figure 5.

The trajectories of the system of nonlinear equations projected on the $Y_1,Y_2$ -plane when revised Rayleigh number (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ when $\mathit{Da}=0.05$ , $\mathit{Les}=100$ , $\sigma =0.05$ , $\epsilon =0.04$ , $Q=10$ , $\mathit{Du}=0.2$ , $\mathit{Sr}=0.3$ , $\mathit{Le}=2$ , $N=25$ and $\mathit{Rn}=5$ .

Figure 6.

Trajectories of the system of nonlinear equations projected on the $\left(Y_1,Y_3\right)$ plane when the revised Rayleigh number (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ when $\mathit{Da}=0.05$ , $\mathit{Les}=100$ , $\sigma =0.05$ , $\epsilon =0.04$ , $Q=10$ , $\mathit{Du}=0.2$ , $\mathit{Sr}=0.3$ , $\mathit{Le}=2$ , $N=25$ and $\mathit{Rn}=5$ .

Figure 7.

Trajectories of the system of nonlinear equations projected on the $\left(Y_1,Y_5\right)$ -plane showing the sensitive dependence of the trajectories on the revised Rayleigh number for (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ when $\mathit{Da}=0.05$ , $\mathit{Les}=100$ , $\sigma =0.05$ , $\epsilon =0.04$ , $Q=10$ , $\mathit{Du}=0.2$ , $\mathit{Sr}=0.3$ , $\mathit{Le}=2$ , $N=25$ and $\mathit{Rn}=5$ .

Figure 8.

Trajectories of the system of nonlinear equations projected on the $\left(Y_1,Y_6\right)$ -plane showing the sensitive dependence of the trajectories on the revised Rayleigh number for (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ when $\mathit{Da}=0.05$ , $\mathit{Les}=100$ , $\sigma =0.05$ , $\epsilon =0.04$ , $Q=10$ , $\mathit{Du}=0.2$ , $\mathit{Sr}=0.3$ , $\mathit{Le}=2$ , $N=25$ and $\mathit{Rn}=5.$

Figure 9.

The bifurcations in the three-dimension solution space $\left(Y_1,Y_2,Y_3\right)$ for (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ when $\mathit{Da}=0.05$ , $\mathit{Les}=100$ , $\sigma =0.05$ , $\epsilon =0.04$ , $Q=10$ , $\mathit{Du}=0.2$ , $\mathit{Sr}=0.3$ , $\mathit{Le}=2$ , $N=25$ and $\mathit{Rn}=5$ .

Figure 10.

Flow trajectories and bifurcations in the three-dimensional space $\left(Y_1,Y_2,Y_6\right)$ for Rayleigh numbers (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ .

Figure 11.

Flow trajectories and bifurcations in the three-dimensional space $\left(Y_1,Y_6,Y_7\right)$ for Rayleigh numbers (a) $R={10}^2$ , (b) $R={10}^3$ and (c) $R={10}^4$ .

Figure 12.

The pattern of streamlines (top), isotherms (middle) and isoconcentration (bottom) for different values of the buoyancy ratio N.

Figure 13.

The streamlines (top), isotherms (middle) and isoconcentration (bottom) for different values of the Darcy number Da.

Figures 5–8 show the phase portraits projected onto the $\left(Y_i,Y_j\right)$ - plane correspond to a simple spiral for $R=100$ . As $R$ is increased to 10⁴, the complexity of the trajectories increases leading to certain chaotic forms. Figures 8–11 show the trajectories in the three-dimensional phase space. Here, we observe similar solution sets as in the two-dimensional phase portraits.

Figures 12 and 13 show the streamline, isotherm and isoconcentration contours in the nanofluid flow for different values of the Darcy number and buoyancy ratio. Figure 12 displays the streamlines for various values of the buoyancy ratio term. Two different eddies are observed. The clockwise and anticlockwise flows are shown via negative and positive stream function values, respectively. The anticlockwise rotating flow occupies the largest area of the nanofluid layer.

For low buoyancy ratio parameters, the flow structure is significantly influenced by the buoyancy within the whole enclosure. Increasing the buoyancy ratio causes the boundary layer thickness to become thinner. Also, a high buoyancy ratio changes the flow structure, and this impacts significantly on the concentration field, which builds up a vertical stratification in the enclosure. It is interesting to note that for $N=-25$ , the effect of the solutal buoyancy force is in the opposite direction of the thermal buoyancy force. The isothermal and isoconcentration profiles are situated toward the left wall, while for $N=1$ , the thermal and solutal buoyancy forces are equal. For $N=25$ , the effect of solutal buoyancy force is in the same direction as the thermal buoyancy force. In such cases, the isothermal and isoconcentration contours are mostly toward the right wall.

We observe that when $N=-25$ , the stream function values in the central eddies increase because the thickness of the boundary layer increases with the buoyancy ratio. The streamlines and the flow behavior are affected by the change in the buoyancy ratio, but the flow pattern remains unaltered. As N decreases from 1 to −25, the streamlines become very dense to the left side of nanofluid layer while when N increases from 1 to 25, the streamlines are less so. The buoyancy forces that drive the nanofluid motion are mainly due to the temperature gradient.

Three different types of eddies are observed for the isoconcentration contours when $N=25$ . Of these, two have a clockwise rotation and one is anticlockwise. It is seen that the small eddy at the right bottom edge is diminished as N decreases from 1 to $-25$ . Here, the concentration boundary layer decreases due to increasing N values, hence the buoyancy ratio has a significant influence on the concentration gradient. As the buoyancy ratio N increases from 1 to 25 the isoconcentrations become very dense at the bottom of nanofluid layer.

The effect of the Darcy number on the nanofluid flow in the porous medium is shown in detail in Figure 13. The streamline patterns are similar to those depicted in Figure 12. However, as Da increases from 0.05 to 0.07, the rotation of the streamlines changes. Similarly, the isotherm patterns change with increasing Darcy numbers. The value of the center eddies increases with increasing Da. Increasing Da has the effect of increasing the effective fluid viscosity and reducing the thermal and solutal boundary layers.

7. Conclusion

We have investigated the onset of thermal instability in a horizontal porous layer of infinite extent in a cross-diffusive nanofluid flow. The focus of the study has been on stress free boundary conditions with zero nanoparticle flux at the wall. A multidomain spectral collocation method was used to solve the system of nonlinear evolution equations. As the Rayleigh number increases to 10⁴, the trajectories spiral many times before reaching a fixed point. The nanofluid convection regime is complex for Rayleigh numbers higher than $R={10}^4$ , and the flow pattern presents difficulties in interpreting correctly.

Additionally, a change in system parameters, such as an increase in the flow Lewis number, improves the rate of heat and mass transfer in the nanofluid saturated porous media. The Dufour parameter has the effect of increasing heat transfer, while increasing the Soret parameter increases the rate of mass transfer.

Acknowledgments

The authors are grateful to the University of KwaZulu-Natal and the Claude Leon Foundation, South Africa for financial support.

The terms $a_{\left(j,n\right)r}^i$ and $R_{\mathit{jr}}^i$ for $j=1,2,\dots ,7$ in Eq. (38) are given by

$$\begin{array}{lll}{a}_{\left(1,1\right)r}^i& =& -{\gamma }_{3}B,a_{\left(1,2\right)r}^i=-{\gamma }_3,a_{\left(1,4\right)r}^i=-{\gamma }_{3}N,a_{\left(1,6\right)r}^i=\frac{{\gamma }_3{N}_A\mathit{Rn}}{\mathit{Ra}},\\ a_{\left(1,3\right)r}^i& =& a_{\left(1,5\right)r}^i=a_{\left(1,7\right)r}^i=0,a_{\left(2,1\right)r}^i=R-Y_{3,r}^i,a_{\left(2,2\right)r}^i=-1,a_{\left(2,3\right)r}^i=-Y_{1,r}^i,\\ a_{\left(2,4\right)r}^i& =& -\mathit{Du},a_{\left(2,5\right)r}^i=a_{\left(2,6\right)r}^i=a_{\left(2,7\right)r}^i=0a_{\left(3,1\right)r}^i=-\frac{\mathit{\alpha \pi }}{2}{u}_{2,r}^i,a_{\left(3,2\right)r}^i=\frac{1}{2}{Y}_{1,r}^i,\\ a_{\left(3,3\right)r}^i& =& -G,a_{\left(3,5\right)r}^i=-{\mathit{\text{GDua}}}_{\left(3,4\right)r}^i=a_{\left(3,6\right)r}^i=a_{\left(3,7\right)r}^i=0,a_{\left(4,1\right)r}^i=R-Y_{5,r}^i,\\ a_{\left(4,2\right)r}^i& =& -\mathit{Sr},a_{\left(5,3\right)r}^i=-\mathit{GSr}\\ a_{\left(4,4\right)r}^i& =& -\frac{1}{\mathit{Les}},a_{\left(4,5\right)r}^i=-Y_{1,r}^i,a_{\left(4,3\right)r}^i=a_{\left(4,6\right)r}^i=a_{\left(4,7\right)r}^i=0,a_{\left(5,1\right)r}^i=\frac{1}{2}{Y}_{4,r}^i,\\ a_{\left(5,4\right)r}^i& =& \frac{1}{2}{Y}_{1,r}^i,a_{\left(5,5\right)r}^i=-\frac{G}{\mathit{Les}},a_{\left(5,2\right)r}^i=a_{\left(5,6\right)r}^i=a_{\left(5,7\right)r}^i=0,a_{\left(6,1\right)r}^i={\gamma }_1-{\gamma }_2{Y}_{7,r}^i,\\ a_{\left(6,2\right)r}^i& =& \frac{\sigma }{\mathit{Le}},a_{\left(6,6\right)r}^i=-\frac{\sigma }{\mathit{Le}},a_{\left(6,7\right)r}^i=-{\gamma }_2{Y}_{1,r}^i,a_{\left(6,3\right)r}^i=a_{\left(6,4\right)r}^i=a_{\left(6,5\right)r}^i=0,\\ a_{\left(7,1\right)r}^i& =& \frac{{\gamma }_2}{2}{Y}_{6,r}^i,a_{\left(7,3\right)r}^i=\frac{G}{\mathit{Le}},a_{\left(7,6\right)r}^i=\frac{{\gamma }_2}{2}{Y}_{1,r}^i,a_{\left(7,7\right)r}^i=-\frac{G}{\mathit{Le}},\\ a_{\left(7,2\right)r}^i& =& a_{\left(7,4\right)r}^i=a_{\left(7,5\right)r}^i=0,\\ R_{1r}^i& =& 0,R_{2r}^i=-Y_{1,r}^i{Y}_{3,r}^i{R}_{3r}^i=\frac{1}{2}{Y}_{1,r}^i{Y}_{2,r}^i,R_{4r}^i=-Y_{1,r}^i{Y}_{5,r}^i,R_{5r}^i=\frac{1}{2}{Y}_{1,r}^i{Y}_{4,r}^i,\\ R_{6r}^i& =& -\frac{\sigma N_A}{\epsilon }{Y}_{1,r}^i{Y}_{7,r}^i,R_{7r}^i=\frac{\sigma N_A}{2\epsilon }{Y}_{1,r}^i{Y}_{6,r}^i,\end{array}$$

where ${\gamma }_1=\frac{\mathit{\sigma R}}{\epsilon }$ , ${\gamma }_2=\frac{N_A\sigma }{\epsilon }$ and ${\gamma }_3=\frac{\mathit{Pr}}{\mathit{Da}}$ .

The matrices $A_{\mathit{ij}}$ in Eq. (42) are given by

$$\begin{array}{c}{A}_{\mathit{nn}}=\mathbf{diag}\left(a_{n,n)r}^i\right)-\mathbf{D},A_{12}=\mathbf{diag}\left(a_{\left(1,2\right)r}^i\right),A_{16}=\mathbf{diag}\left(a_{\left(1,6\right)r}^i\right),\\ A_{13}=A_{14}=A_{15}=A_{17}=\mathbf{O},\\ A_{21}=\mathbf{diag}\left(a_{\left(2,1\right)r}^i\right),A_{23}=\mathbf{diag}\left(a_{\left(2,3\right)r}^i\right),A_{24}=\mathbf{diag}\left(a_{\left(2,4\right)r}^i\right),A_{25}=A_{26}=A_{27}=\mathbf{O},\\ A_{31}=\mathbf{diag}\left(a_{\left(3,1\right)r}^i\right),A_{32}=\mathbf{diag}\left(a_{\left(3,2\right)r}^i\right),A_{35}=\mathbf{diag}\left(a_{\left(3,5\right)r}^i\right),A_{34}=A_{36}=A_{37}=\mathbf{O},\\ A_{41}=\mathbf{diag}\left(a_{\left(4,1\right)r}^i\right),A_{42}=\mathbf{diag}\left(a_{\left(4,2\right)r}^i\right),A_{45}=\mathbf{diag}\left(a_{\left(4,5\right)r}^i\right),A_{43}=A_{46}=A_{47}=\mathbf{O},\\ A_{51}=\mathbf{diag}\left(a_{\left(5,1\right)r}^i\right),A_{53}=\mathbf{diag}\left(a_{\left(5,3\right)r}^i\right),A_{54}=\mathbf{diag}\left(a_{\left(5,4\right)r}^i\right),A_{52}=A_{56}=A_{57}=\mathbf{O},\\ A_{61}=\mathbf{diag}\left(a_{\left(6,1\right)r}^i\right),A_{62}=\mathbf{diag}\left(a_{\left(6,2\right)r}^i\right),A_{67}=\mathbf{diag}\left(a_{\left(6,7\right)r}^i\right),A_{63}=A_{64}=A_{65}=\mathbf{O},\\ A_{71}=\mathbf{diag}\left(a_{\left(7,1\right)r}^i\right),A_{73}=\mathbf{diag}\left(a_{\left(7,3\right)r}^i\right),A_{76}=\mathbf{diag}\left(a_{\left(7,6\right)r}^i\right),A_{72}=A_{74}=A_{75}=\mathbf{O},\end{array}$$

where $\mathbf{O}$ is an $\left(N+1\right)\times \left(N+1\right)$ matrix of zeros and diag is an $\left(N+1\right)\times \left(N+1\right)$ diagonal matrix.