Abstract
The linear hydrodynamics of Rayleigh convection in a horizontal nanofluid layer heated from below was studied. The hydrodynamic stability of the fluid layer bounded by two horizontal perfect thermal conducting walls was extended to analyze steady and oscillatory convection, and the role played by thermophoresis. Experimental data of TiO2 particle–based nanofluid was used to discuss the stability of the fluid layer. Results on the relationship between thermal and volume faction Rayleigh numbers are used to discuss experiments in nanofluid Rayleigh convection, while the absence of thermophoresis in the model equations was also considered. For this nanofluid, steady convection sets in at critical wavenumber ac = 3.12, but thermal RT and nanoparticle volume fraction RV Rayleigh numbers are given by an implicit relationship. For the onset of oscillatory convection, the wavenumber is also obtained from an implicit equation involving RT and RV. Results are discussed in terms of physical dimensionless parameters of the system like the Lewis and Prandtl numbers. This work complements the earlier efforts of Tzou and more recently by Nield and Kuznetsov.
Introduction
Heat transfer enhancement by nanofluids has called the attention of researchers since several years ago due to industrial and lab applications.1–4 Most of these applications are related to heat transfer enhancement with a number of theoretical5,6 and experimental results reported. For example, Garoosi et al. 7 used the nanofluid Buongiorno’s 8 model to investigate the heat transfer process in heat exchangers, Mahian et al. 2 made a wide review on the application of nanofluids to solar collectors and solar water heaters, Rashidi et al. 9 made a theoretical study about heat transfer enhancement using nanofluids in applications involving rotating flow, and in 2012 Sivashanmugam 3 provided an overview on the usage of nanofluids for heat transfer through the compilation of some theoretical and experimental works. Although nanofluids find most of its applications in heat transfer enhancement, Wong and De Leon 4 envisaged some other interesting applications in biomedicine and food.
Rayleigh convection phenomenon in nanofluids has been considered since several years ago. Buongiorno, 8 for example, made a review on convective heat transport by nanofluids. Later, Tzou10,11 formulated the problem of Rayleigh convection for fixed temperature and particle volume fraction for different mechanical boundary conditions. However, Tzou10,11 only considered steady convective motions in his theoretical studies. Next, Nield and Kuznetsov 12 extended the results of Tzou10,11 by adding a diffusion equation for the nanoparticles concentration implying an extra destabilizing agency. More recently, Nield and Kuznetsov 13 approached the thermoconvective instability of a nanofluid layer for a different set of boundary conditions for the particle volumetric fraction. In that case, 13 the new boundary conditions were based on zero mass flux at the walls. Difficulties to measure and keep constant volumetric fractions at the walls justify the change in the boundary conditions proposed by Nield and Kuznetsov. 13 Despite the reasonable modification in the volumetric fraction boundary conditions, a number of terms, accounting for cross thermal-nanoparticle volumetric fraction effects and the basic state, were omitted in the working model so that their results on the hydrodynamic stability of the system were not complete.
Experimental reports on the convective heat transport enhancement justify, in part, the present investigation to include nanofluids in the Rayleigh convection. A number of investigations have been conducted to study the thermal conductivity of nanofluids under diverse conditions (e.g. see the review of Das et al. 14 ). On the other hand, experiments were performed on natural convection in nanofluids,15–20 where the interactions of physical properties and mechanisms are analyzed, which motivated this theoretical investigation too. Some experimental observations about the nanofluids’ heat transfer enhancement effects on the Rayleigh convection are contradictory: in some cases, fluid instability is favored by nanofluids15,17 while in others the opposite effect was found,17–20 among others. More recently, in a numerical investigation, Haddad et al. 21 studied the role played by Brownian motions and thermophoresis in the heat transfer enhancement in Rayleigh convection for CuO-water nanofluids and found that enhancement and deterioration of heat transfer are possible for different nanoparticle volume fractions. Further understanding of the hydrodynamics of nanofluids layers is then needed.
This study focuses on the linear thermal hydrodynamic stability of a nanofluid layer heated from below. In this approach, the dimensionless groups, for the problem in hand, have been reorganized to concentrate the information in more tractable parameters. Thermophoretic influence is also investigated, in view of previous experimental and theoretical investigations, aiming to explain heat transfer enhancement and/or deterioration. Besides, missing terms of previous works have been considered for the hydrodynamics of the fluid layer, so that new results and an extended discussion could be made. On the other hand, some interesting features of the basic state are briefly discussed in terms of its nonlinearity too. The Galerkin method was used to approach the eigenvalue problem for the thermal Rayleigh number so that analytical expressions for the onset of steady and oscillatory convection could be retrieved. A numerical analysis was also conducted in order to confirm the analytical calculations, and a very good agreement was found. The boundary value problem was then studied by the shooting method.
The article is organized as follows. In the section of governing equations, the mathematical formulation of the problem is given along with proper nondimensionalization of the equations of the problem while the subsequent section presents the boundary conditions for the problem in hand. The unperturbed temperature and nanoparticle volumetric fraction profiles are exposed in the basic state section. Sections “The perturbation equations” and “Linear stability analysis” are devoted to provide perturbation of the linear hydrodynamic stability of the nanofluid layer. Results and discussion, and conclusions are presented in the last two sections, respectively. The analytical calculations in this article were performed using MapleTM.
The governing equations
Because the aim of this research is to understand the role of nanofluids in the Rayleigh convection problem, for the study of the linear thermoconvective hydrodynamic stability in a fluid layer with suspended nanoparticles, the arrangement shown in Figure 1 is then considered. Since in some Rayleigh convection experiments working nanofluids behave as Newtonian fluids 20 for low particle volume fractions, this feature is assumed here too. Agglomeration, coagulation, and particle cohesion and adhesion to the walls are omitted in this theoretical study since lab nanofluids preparation techniques may reduce these effects.15,20

Schematics for the problem of Rayleigh convection in a horizontal fluid layer heated from below and with suspended particles. The bounding walls are perfect thermal conductors. The wall thicknesses are shown for visual purposes only since in the present analysis these are not considered in the model equations. Also, convection cells were sketched in (a) to depict the horizontal velocity profile (left) and the thermal boundary layer (right) in (b).
The system consists of an infinite horizontal fluid layer heated from below and bounded by two perfect thermal conducting parallel walls located at z* = 0,h as sketched in Figure 1. The model used by Tzou10,11 and more recently by Nield and Kuznetsov 13 which includes Brownian motions and thermophoretic effect is used since these effects may help to give more insight on the heat transfer enhancement mechanism. Due to missing terms in the model equations of Nield and Kuznetsov 13 and due to a different, particle volumetric fraction, boundary condition used by Tzou,10,11 new results are obtained from the present revision and extension.
The model for the problem of natural convection in a fluid layer consists of governing equations for the mass conservation, the momentum balance, the heat conduction, and the nanoparticle volume fraction. In dimensional form, these are
where
The above governing equations (3) and (4) account for the Brownian and thermophoretic effect through the coupling terms indicated by the diffusion coefficients DB and DT, respectively. In equations (1)–(5),
where β is the volumetric expansion coefficient. The density equation (6) shows nonlinearities on temperature T* and particle volume fraction n* which can be translated as a cross dependence between these two variables. From equation (6), further modifications to the fluid layer instabilities are envisaged. Next, equations (1)–(6) shall be studied in dimensionless form in order to obtain groups of parameters that may help to explain the fluid hydrodynamics. Then, for the nondimensionalization process, the following scales are used: h for lengths, h2/κ for time, κ/h for fluid velocity, and κµ/h2 for pressure, with κ as the nanofluid thermal diffusivity. Finally, the nanoparticle volume fraction and the temperature shall be written in dimensionless form as n = (n* – n0)/n0 and
where the dimensionless groups are defined as follows
where the * has been dropped out to represent dimensionless variables in equations (7)–(10). Pr is Prandtl number, RT is the thermal Rayleigh number, Le is the Lewis number, NA is the ratio of diffusion coefficients (thermophoretic/Brownian motions), NB is the ratio of volumetric heat capacities (nanoparticles/base fluid), and RV is the particles Rayleigh number. Next, dimensionless groups in equations (7)–(10) have been organized in a different way to previous authors,10–13 despite the present n* dimensionless form. The main idea is the usage of a smaller number of dimensional groups from which physical insight derived from the analytical calculations, and numerical validation, could benefit. For example, Tzou10,11 defined the thermal Rayleigh number, in his equation (16), in terms of the square of the thermal diffusivity and introduced other parameters, probably, of low significance to his investigations. In the work of Nield and Kuznetsov, 13 the absence of a cross effect source in the buoyancy terms of their equation (10) and the presence of a complex dimensionless group Rm are observed. In this work, only thermal and volume fraction Rayleigh numbers are identified.
Terms –n0RTT
After substitution of the perturbed variables given in equation (11) into equations (7)–(10), two systems of equations are obtained. One for the basic state in which the fluid is at rest and temperature and volume fraction depend only on z, and another one corresponding to the perturbed governing equations.
Boundary conditions
The nanofluid layer is bounded by two infinite horizontal rigid and non-permeable walls. Then, thermal and mechanical nature of the bounding walls determines the fluid hydrodynamic instability. Based on the previous statement particle volume fraction, temperature and fluid velocity should satisfy proper boundary conditions at z*= 0, h in dimensional form according to the scheme in Figure 1.
For the nanoparticles volume fraction, the boundary condition is obtained by considering zero mass flux of particles at the walls. 13 In dimensional form this is
According to Nield and Kuznetsov 13 boundary conditions, equation (12) may be more tractable in comparison to other studies.10–12 Zero flux mass at the walls could be, experimentally, implemented in an easier way than keeping constant volume fraction at the boundaries. Therefore, Rayleigh number dependence on the other parameters of the system changes. Here, it is considered that the bounding walls are perfect thermal conductors so that fixed temperature boundary conditions are used. Since the walls are rigid too, the fluid velocity in these places must be zero because of the adherence condition.
In particular, for the basic state indicated with the subscript b, the boundary conditions in dimensionless form are
The basic state
In the basic state indicated with the subscript b, the fluid is at rest and temperature and particle volumetric fraction depend on z only. Thus, the basic state is determined by the following system of ordinary differential equations for pb, Tb, and nb
Most attention is given to the temperature profile Tb and to the volumetric fraction profile nb. Because Tb and nb are independent of pb, these are first calculated from equations (16) and (17) and later substituted into equation (15) from which the solution for pb is straight. Since pressure and velocity fields are later split, pb shall be omitted, but it can be obtained from authors upon request. Notice that system of equations (16) and (17) is nonlinear and different solutions for Tb and nb may exist. Also, the linear stability of two different basic states, for temperature and volumetric fraction, have been analyzed by Nield and Kuznetsov12,13 and Tzou,10,11 but there is still another set of possible basic states for this system. These are
where C1 through C4 are integration constants. The subscript a = 1,2 in equation (20) so that equations (18) and (19) generate two different temperature profiles. Basic states in equations (18)–(20) are remarkable since it means that the nanofluid may have different temperature and volumetric fraction profiles before perturbation. In other words, it is possible that for some basic states the nanofluid layer destabilizes at lower critical conditions. However, in a brief comparison between the basic states shown in equations (18)–(20) and those previously used by Nield and Kuznetsov12,13 and Tzou,10,11 for a TiO2 particle–based nanofluid (see Table 1), the volume fraction profile equations (18)–(19) showed non-physical behavior. Further comparisons for the temperature profiles10,11 have shown almost linear behavior for small volumetric fractions and this temperature basic state becomes nonlinear as the volumetric fraction increases. The same deviation from linear behavior was found for the volumetric fraction basic state.
kB is the Boltzmann constant,
Parameters βB, DT, and DB were calculated through equations in the work by Buongiorno. 8
Water data were taken.
Δn* = 57–19 was taken for demonstration based on the data of Wen and Ding. 19
This is the highest temperature value for 0.57% used in Wen and Ding. 19
Figure 2 shows plots of nb of previous works for the case of TiO2 nanofluid. Plots in Figure 2 are important for understanding of previous approaches and since it represents the initial volumetric fraction distribution different hydrodynamics could be obtained after perturbation. As the nature of suspended nanoparticles, preparation of the nanofluid, and other experimental conditions are unknown, profiles in Figure 2 could still have physical importance. Here, the main concern is with the basic state reported by Nield and Kuznetsov 13 because of its simplicity and need of further comparison in view of the present extension to their results. Then, Tb and nb are expressed as
where

The perturbation equations
The perturbed governing equations, for
where w1 is the fluid velocity vertical component. Equations (23)–(24) show the coupling of temperature and particle volumetric fraction which can be easily identified in the right-hand side terms by its product with the basic states. Notice that parameters Le, NA, and NB identify the Brownian and thermophoretic effects, and despite the nonlinear source of the right-hand side terms of equation (23), only linearized contributions appear. Convective motions in the nanofluid layer are assumed to occur as regular patterns so it is reasonable to write the perturbed variables using normal modes. Thus, w1, T1, and n1 are rewritten as
where σ is a complex parameter whose real part σR is the perturbations growth rate and the imaginary part ω is the perturbations frequency of oscillation. The model equations (23)–(25) are then reduced to the following system of ordinary differential equations
where D = d/dz and
Equations (27)–(30) represent an eigenvalue problem for the thermal and particle volume fraction Rayleigh numbers that shall be analytically and numerically studied in the following sections.
Linear stability analysis
The linear hydrodynamic stability analysis is now investigated for a plausible situation in an experimental setup or application. The aim is then to determine the critical conditions at which convective motions in the nanofluid layer set in. In other words, critical values for the Rayleigh and wavenumber shall be determined in this section.
The present stability analysis extends the results of previous studies. 13 Analytical and numerical techniques are used to give certainty to the findings of this study. In the case of numerical calculations, the eigenvalue problem given by equations (27)–(30) was treated as follows. First, the differential equations are integrated by the Runge–Kutta method while the boundary problem was solved with help of the shooting method. 25 Some Fortran routines from the book of Press et al. 25 were used to perform the numerical computations. The shooting method uses trials, starting at the first boundary, to check consistency between values of dependent variables at the second boundary with those obtained after integration of the differential equations. Numerical results are used to validate the analytical formulas for the hydrodynamics while analytical results would also help to guide the numerical analysis.
Then, the Galerkin method is used to approach the eigenvalue problem given by equations (27)–(30) because of its good accuracy at low-order approximations.26,27 Also, the Galerkin method allows to find an analytical expression of the Rayleigh number without the need of solving the system of differential equations. Therefore, the following set of trial functions, for W, Θ, and Φ, is considered
which also satisfy the corresponding boundary conditions. Next, trial functions equation (31) are substituted into the perturbed model equations (27)–(30) to form the residual which is also made orthogonal to the respective trial functions. A matrix of coefficients is then formed and its determinant gives a necessary and sufficient condition to be satisfied. In this way, this condition is given in the form of the following determinant
where
The i, j counters run from 1 to N, the order of approximation. The size of the matrix in equation (32) increases with N and consequently the analytical results become larger along with a higher number of roots for the Rayleigh number. For the first approximation, N = 1, good accuracy can be obtained from the determinant equation (32). This is
where
Oscillatory and steady convection can then be investigated from equation (33). Although thermal and particle volumetric fraction forces may induce convective motions in the fluid layer, the first one may be a stronger destabilizing mechanism than the last one. The previous statement may apply for experimental conditions with small
Steady convection
For steady convection, there is no temporal dependency so that ω = 0 is set in equation (33). In this way, equation (29) reduces to the following relationship
where Brownian motions and thermophoretic effect can be identified by Le and NA. Notice that the ratio of volumetric heat capacities NB do not contribute to steady or oscillatory convection, as seen in equation (33). Next, the onset of convection can be determined by minimization of the thermal Rayleigh number as function of the wavenumber in equation (36) which leads to ac = 3.12 and
In other words, since thermal Rayleigh number
Oscillatory convection
In order to study the onset of oscillatory convection, RT is isolated from equation (33). This can be written symbolically as
where f1 and f2 depend on ω, a, RV, Pr, Le, NA, n0, and
where
Notice that A1, A2, A3 > 0 provided that the wavenumber is nonzero and
Next, the thermal Rayleigh number for oscillatory convection
where ω2 and ω4 should be obtained from equation (39).
In the absence of thermophoresis
The absence of thermophoresis in nanofluid convection is now investigated to discuss the effect of it on the nanofluid instability. Then, the nanoparticle flux is reduced to
in which some further simplifications were made due to the new basic state (nb) B . Since the thermophoretic effect was eliminated, the boundary conditions for (nb) B were modified too. This is
Boundary condition equation (45) leads to a constant particle volume fraction profile. On the other hand, the momentum balance equation (27) remains the same because this is not affected by
At the first order of approximation, with the Galerkin method applied on the eigenvalue problem equations (27), (43), and (44), it is obtained
The absence of RV in equation (48) is noticeable, but the effect of particle volume fraction still remains through
and the critical wavenumber is acB = 3.12, again. On the other hand, in the absence of thermophoresis oscillatory convection vanishes. ω = 0 was expected because of the only iRTω product in equation (48).
Results and discussion
The linear stability of a nanofluid layer heated from below was investigated for the case of perfect thermal conducting walls. The results are given in terms of dimensionless parameters of the system like the thermal and particle volume fraction Rayleigh numbers, the wavenumber, the frequency of oscillation, the Prandtl number, the average particle volume fraction, a reference value for the particle volume fraction, and a set of parameters to represent the double diffusion across the fluid layer. The particular case of a TiO2 nanofluid was addressed in order to give some perspective to the present results.
The results were derived from a simplified nondimensionalization process and a revisited linear analysis of the natural convection problem.10–13 It was found that both steady and oscillatory convection are possible destabilizing mechanisms of the fluid layer and that the average
Steady and oscillatory convection
The present analysis shows that steady and oscillatory convection is possible, for positive RT and RV. It has to be mentioned that considered values for
For the case of steady convection, the linear relationship between the thermal and volumetric fraction Rayleigh numbers is graphically shown in Figure 3. The corresponding critical wavenumber ac = 3.12 could be determined explicitly which helps to visualize the competition between thermal and volume fraction destabilizing forces. The ratio of diffusion coefficients NA in equation (37) has two sources: the trial function and the thermophoretic and Brownian motion contribution terms in equations (3) and (4), so that confusion should be avoided.

Plot of the linear relationship between critical thermal and particle volumetric fraction Rayleigh numbers given by equation (37). The data provided in Table 1 was used too.
From the curve of criticality shown in Figure 3, it is observed that for a given critical
The lower values of the Rayleigh particle volume fraction RV in comparison with RT may be due to the large density difference between nanoparticles and fluid as pointed out by Wen and Ding.
19
Despite the possible sign change in the term
The case of oscillatory convection is more complex than that for steady convection. Temporal dependence introduces the Prandtl number Pr and the frequency of oscillation ω as new parameters in the fluid instability. Main results on oscillatory stability were obtained from the corresponding Rayleigh number

Plot of the thermal Rayleigh number against the wavenumber given by equation (42). The data provided in Table 1 was used too.

Plot of the particle volume fraction Rayleigh number against the wavenumber given by equation (42). The data provided in Table 1 was used too.

Plot of the frequency of oscillation w against particle volumetric fraction Rayleigh numbers given by equation (39). The data provided in Table 1 was used too.
The critical wavenumber could not be explicitly calculated since minimization of equation (42) leads to an implicit relationship between a and
For the TiO2 nanoparticle fluid, it was found that in oscillatory convection the curve
The role of thermophoresis
An interesting point of discussion, on the light of the present results, is the role played by the thermophoretics. For short, temperature and particle volume fraction variables would partially decouple. This is concluded after setting
However, despite the partial decoupling of the governing equations for temperature and particle volumetric fraction, convection may still occur. Only the thermal Rayleigh number remains in the solvability condition reducing the particle volume fraction influence to the parameters
Conclusion
The problem of convection in a nanofluid layer heated from below was studied both analytically and numerically. The used boundary condition for the particle volume fraction was based in the volume fraction flux involving thermophoretic and Brownian motions. Steady and oscillatory convection results were discussed for the particular case of TiO2 nanoparticle fluid.
For TiO2 nanoparticle-based fluids steady convection sets in at critical wavenumber ac = 3.12, but thermal RT and nanoparticle volume fraction RV Rayleigh numbers cannot be easily determined since these are given by an implicit relationship. For steady convection, it is concluded that the whole fluid instability due to the presence of nanoparticles is strongly modified by it. For the critical steady conditions at ac = 3.12, the thermal Rayleigh number RT is always larger than RV. It could be interpreted as if convective motions were triggered by density differences between the particles and the fluid, with thermal gradients playing a complementary role. It is not easy to determine if larger TiO2 nanoparticle volume fractions would change, in opposite direction, the competence with thermal gradients to destabilize the fluid layer. The RV stronger destabilizing force may ease heat transport because of the generated secondary flow.
For the onset of oscillatory convection critical wavenumber, critical thermal Rayleigh number and critical nanoparticle volume fraction Rayleigh number are also given by an equation so that two of them should be fixed in order to determine the third one. When convection sets in as oscillatory motions, the magnitude of RV is smaller than that of RT but with a nonlinear relationship between them, in comparison with steady convection. Here, fluid oscillations play an interesting role in the fluid instability. The frequency of oscillation ω increases with the wavenumber a, while at the same time, RV decreases and RT increases. The increasing of the frequency of oscillation may be due to the growing number of convective cells as indicated by the wavenumber. On the other hand, the oscillations in the fluid may dissipate the heat so that RT increases as shown in Figure 4. The absence of thermophoresis leads to instability conditions in which particle volume fraction influence is reduced. Further discussion in terms of experimental results may be needed to explore the present linear analysis.
Footnotes
Handling Editor: Rui Lima
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.
