Open access peer-reviewed chapter

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

By Osman A.I. Noreldin, Precious Sibanda and Sabyasachi Mondal

Submitted: May 29th 2017Reviewed: September 19th 2017Published: December 20th 2017

DOI: 10.5772/intechopen.71066

Downloaded: 548


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.


  • nonlinear instability
  • nanofluid flow
  • porous medium
  • multidomain spectral collocation method

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 Al2O3 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=0and 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 Tcand Th, respectively with Th>Tc. 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


subject to the boundary conditions


where Vis the fluid velocity, Tis the temperature, Cis the solute concentration and ϕis the volumetric fraction of nanoparticles. The dimensionless parameters are the Darcy number (modified by the viscosity ratio) Da, Prandtl number Pr, Hartmann-Darcy number Q, thermal Rayleigh-Darcy number Ra, nanoparticle Rayleigh number Rnand the basic density Rayleigh number Rm. The parameter NAis a modified diffusivity ratio, Leis the Lewis number, Rsis solutal Rayleigh number, NBis a modified nanoparticle density increment and Duis a modified Dufour parameter. The parameter Lesis the thermo-nanofluid Lewis number, νis the kinematic viscosity and Sris a modified Soret parameter. These parameters have the form


where ρf,ρp,μ˜,β1,β2,κm,δ, εand Kare 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 gand DBis the Brownian diffusion coefficient, DTis the thermophoresis diffusion coefficient, DSis the solutal diffusion coefficient, DTCis the Dufour parameter and DCTis the Soret parameter. The heat capacity of the fluid is ρcf, ρcpis the effective heat capacity of the nanoparticle, ρcmis the effective heat capacity of the porous medium and B0is 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


3. Weakly nonlinear stability analysis

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


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


where A11,B11,B02,C11,C02,D11and D02are amplitudes that depend on time. This leads to the Lorenz-type system of nonlinear ordinary differential equations


subject to Yn0=Yn0for n=1,2,,7.The following variables have been introduced in the equations above:


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


This is always true if B0. 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 V0at time t=0, after a time t>0, the end point of the corresponding trajectories fills a volume


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


We obtain the possible stationary points of the nonlinear system of equations by setting Ẏi=0for i=1,2,,7. One of these stationary points is Yi=0and by linearizing about this point, we obtain the Jacobian matrix


The eigenvalues of the above matrix depend on the various parameters. For the specific parameters R=103,Da=20,Pr=10,N=25,Du=0.2,Sr=3,Les=10,Le=5,σ=0.05,G=3 andε=0.04, the characteristic polynomial is


with eigenvalues

λ1=0.2955056985, λ2=0.2402382976, λ3=0.6139990637, λ4=4.886000936,λ5=5, λ6=5.722871998121.9033954659i, λ7=5.7228719981+21.9033954659i.

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 Aby Tand the determinant d, then




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


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 0Tinto subintervals Ωi=ti1tifor i=1,2,,p. The transformation


is used to transform each subinterval Ωiinto the interval 11. The system of Eqs. (15)(21) can be written in the form


subject to


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


subject to


where ajnriand Rjrifor j=1,2,,7are 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


We approximate the derivatives of the unknown functions Yn,r+1itat the collocation points by


where D=2D/titi1, Dis the Chebyshev differentiation matrix and


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


where the matrices A=Aijand Rniare 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 Nutwhere


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


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


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 tmax=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 24. 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 numberNuand (b) the Sherwood numberShforDa=0.05,Le=2,Du=0.2,ε=0.04,σ=0.05,Les=100,Rn=5,Ra=1000and various values of the Dufour and Soret parameters.

Figure 3.

The effect of Prandtl numberPron (a) the Nusselt numberNuand (b) the Sherwood numberShwhenDa=0.05,Le=2,Du=0.2,ε=0.04,σ=0.05,Les=100,Rn=5andRa=1000.

Figure 4.

The effect of Lewis number on (a) the Nusselt numberNuand (b) the Sherwood numberShwhenDa=0.05,Du=0.2,ε=0.04,σ=0.05,Les=100,Rn=5andRa=1000.

Figures 511 show the effect of the Rayleigh number on the trajectories projected onto the YiYjphase 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 Rhave been presented with the trajectories projected onto the YiYjphase planes. These trajectories spiral toward the fixed point for Rayleigh numbers from 102to 104. The solution sets give spiral phase portraits as Rincreases 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 theY1,Y2-plane when revised Rayleigh number (a)R=102, (b)R=103and (c)R=104whenDa=0.05,Les=100,σ=0.05,ε=0.04,Q=10,Du=0.2,Sr=0.3,Le=2,N=25andRn=5.

Figure 6.

Trajectories of the system of nonlinear equations projected on theY1Y3plane when the revised Rayleigh number (a)R=102, (b)R=103and (c)R=104whenDa=0.05,Les=100,σ=0.05,ε=0.04,Q=10,Du=0.2,Sr=0.3,Le=2,N=25andRn=5.

Figure 7.

Trajectories of the system of nonlinear equations projected on theY1Y5-plane showing the sensitive dependence of the trajectories on the revised Rayleigh number for (a)R=102, (b)R=103and (c)R=104whenDa=0.05,Les=100,σ=0.05,ε=0.04,Q=10,Du=0.2,Sr=0.3,Le=2,N=25andRn=5.

Figure 8.

Trajectories of the system of nonlinear equations projected on theY1Y6-plane showing the sensitive dependence of the trajectories on the revised Rayleigh number for (a)R=102, (b)R=103and (c)R=104whenDa=0.05,Les=100,σ=0.05,ε=0.04,Q=10,Du=0.2,Sr=0.3,Le=2,N=25andRn=5.

Figure 9.

The bifurcations in the three-dimension solution spaceY1Y2Y3for (a)R=102, (b)R=103and (c)R=104whenDa=0.05,Les=100,σ=0.05,ε=0.04,Q=10,Du=0.2,Sr=0.3,Le=2,N=25andRn=5.

Figure 10.

Flow trajectories and bifurcations in the three-dimensional spaceY1Y2Y6for Rayleigh numbers (a)R=102, (b)R=103and (c)R=104.

Figure 11.

Flow trajectories and bifurcations in the three-dimensional spaceY1Y6Y7for Rayleigh numbers (a)R=102, (b)R=103and (c)R=104.

Figure 12.

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

Figure 13.

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

Figures 58 show the phase portraits projected onto the YiYj- plane correspond to a simple spiral for R=100. As Ris increased to 104, the complexity of the trajectories increases leading to certain chaotic forms. Figures 811 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 Ndecreases from 1 to −25, the streamlines become very dense to the left side of nanofluid layer while when Nincreases 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 Ndecreases from 1 to 25. Here, the concentration boundary layer decreases due to increasing Nvalues, hence the buoyancy ratio has a significant influence on the concentration gradient. As the buoyancy ratio Nincreases 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 Daincreases 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 Dahas 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 104, the trajectories spiral many times before reaching a fixed point. The nanofluid convection regime is complex for Rayleigh numbers higher than R=104, 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.


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

The terms ajnriand Rjrifor j=1,2,,7in Eq. (38) are given by


where γ1=σRε, γ2=NAσεand γ3=PrDa.

The matrices Aijin Eq. (42) are given by


where Ois an N+1×N+1matrix of zeros and diagis an N+1×N+1diagonal matrix.

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Osman A.I. Noreldin, Precious Sibanda and Sabyasachi Mondal (December 20th 2017). Weakly Nonlinear Stability Analysis of a Nanofluid in a Horizontal Porous Layer Using a Multidomain Spectral Collocation Method, Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals, Ricardo López-Ruiz, IntechOpen, DOI: 10.5772/intechopen.71066. Available from:

chapter statistics

548total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Small-Angle Scattering from Mass and Surface Fractals

By Eugen Mircea Anitas

Related Book

First chapter

Some Commonly Used Speech Feature Extraction Algorithms

By Sabur Ajibola Alim and Nahrul Khair Alang Rashid

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us