On the linear convective stability of a TiO 2 particle–based nanofluid layer heated from below

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 TiO₂ 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 a_c = 3.12, but thermal R_T and nanoparticle volume fraction R_V Rayleigh numbers are given by an implicit relationship. For the onset of oscillatory convection, the wavenumber is also obtained from an implicit equation involving R_T and R_V. 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.

Keywords

Nanofluid thermal convection hydrodynamic stability

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 theoretical^5,6 and experimental results reported. For example, Garoosi et al.⁷ used the nanofluid Buongiorno’s⁸ model to investigate the heat transfer process in heat exchangers, Mahian et al.² made a wide review on the application of nanofluids to solar collectors and solar water heaters, Rashidi et al.⁹ made a theoretical study about heat transfer enhancement using nanofluids in applications involving rotating flow, and in 2012 Sivashanmugam³ 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⁴ envisaged some other interesting applications in biomedicine and food.

Rayleigh convection phenomenon in nanofluids has been considered since several years ago. Buongiorno,⁸ for example, made a review on convective heat transport by nanofluids. Later, Tzou^10,11 formulated the problem of Rayleigh convection for fixed temperature and particle volume fraction for different mechanical boundary conditions. However, Tzou^10,11 only considered steady convective motions in his theoretical studies. Next, Nield and Kuznetsov¹² extended the results of Tzou^10,11 by adding a diffusion equation for the nanoparticles concentration implying an extra destabilizing agency. More recently, Nield and Kuznetsov¹³ approached the thermoconvective instability of a nanofluid layer for a different set of boundary conditions for the particle volumetric fraction. In that case,¹³ 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.¹³ 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.¹⁴). 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 nanofluids^15,17 while in others the opposite effect was found,^17–20 among others. More recently, in a numerical investigation, Haddad et al.²¹ 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 Maple^TM.

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²⁰ 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

Figure 1.

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 Tzou^10,11 and more recently by Nield and Kuznetsov¹³ 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¹³ 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

\nabla^{2} u^{*} = 0

(1)

ρ_{0} (\frac{\partial u^{*}}{\partial t^{*}} + u^{*} \cdot \nabla u^{*}) = - \nabla p^{*} + μ \nabla^{2} u^{*} + ρ g k

(2)

ρ_{0} c_{0} (\frac{\partial T^{*}}{\partial t^{*}} + u^{*} \cdot \nabla T^{*}) = k \nabla^{2} T^{*} - c_{p} J^{*} \cdot \nabla T^{*}

(3)

\frac{\partial n^{*}}{\partial t^{*}} + u^{*} \cdot \nabla n^{*} = - \frac{1}{ρ_{p}} \nabla \cdot J^{*}

(4)

where J^* represents the mass flux of particles⁸ and is defined as

J^{*} = - ρ_{p} D_{B} \nabla n^{*} - ρ_{p} \frac{D_{T}}{T_{T}^{*}} \nabla T^{*}

(5)

The above governing equations (3) and (4) account for the Brownian and thermophoretic effect through the coupling terms indicated by the diffusion coefficients D_B and D_T, respectively. In equations (1)–(5), u^* is the fluid velocity, p^* is the pressure, T^* is the temperature, n^* is the nanoparticle volumetric fraction, ρ₀ is the base fluid density, µ is the fluid dynamic viscosity, ρ is the nanofluid density, c₀ is the fluid specific heat, k is the fluid thermal conductivity, c_p is the nanoparticle specific heat, ρ_p is the nanoparticle density, t^* is the time, g is the acceleration due to gravity, and $T_{T}^{*}$ is the temperature at the top wall. The * indicates the dimensional form of variables and parameters. The Boussinesq approximation applies to this problem, as pointed out by Wen and Ding,¹⁹ which implies that density changes are only important in the buoyancy term in the momentum balance equation (2) and in the case of natural convection in nanofluids the following density equation of state proposed by Tzou^10,11 is considered

ρ = ρ_{p} n^{*} + ρ_{0} (1 - n^{*}) [1 - β (T^{*} - T_{T}^{*})]

(6)

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^* – n₀)/n₀ and $T = (T * - T_{T}^{*}) / (T_{B}^{*} - T_{T}^{*})$ , respectively. Since parameters related to particle volumetric fraction may not be easily estimated, the choice for the dimensionless form of n^*, with n₀ as a reference value (like the average volume fraction for example), is reasonable. $T_{B}^{*}$ is the temperature at the bottom wall. Thus, in its dimensionless form equations (1)–(4) are rewritten as

\nabla^{2} u = 0

(7)

\begin{matrix} P r^{- 1} (\frac{\partial u}{\partial t} + u \cdot \nabla u) \\ = - \nabla p + \nabla^{2} u - (n_{0} - 1) R_{T} T k - R_{V} n k \\ + R_{T} n_{0} Tn k - (\frac{h^{3} g}{κ ν} + R_{V}) k \end{matrix}

(8)

\frac{\partial T}{\partial t} + u \cdot \nabla T = \nabla^{2} T + \frac{N_{B}}{Le} \nabla n \cdot \nabla T + \frac{N_{A} N_{B}}{Le} \nabla T \cdot \nabla T

(9)

\frac{\partial n}{\partial t} + u \cdot \nabla n = \frac{1}{Le} \nabla^{2} n + \frac{N_{A}}{Le} \nabla^{2} T

(10)

where the dimensionless groups are defined as follows

\Pr = \frac{ν}{κ} N_{A} = \frac{D_{T} Δ T^{*}}{D_{B} T_{T}^{*} n_{0}}

R_{T} = \frac{h^{3} β g Δ T^{*}}{κ ν} N_{B} = \frac{ρ_{p} c_{p} n_{0}}{ρ_{0} c_{0}}

Le = \frac{κ}{D_{B}} R_{V} = \frac{h^{3} (ρ_{p} - ρ_{0}) g n_{0}}{κ μ}

where the * has been dropped out to represent dimensionless variables in equations (7)–(10). Pr is Prandtl number, R_T is the thermal Rayleigh number, Le is the Lewis number, N_A is the ratio of diffusion coefficients (thermophoretic/Brownian motions), N_B is the ratio of volumetric heat capacities (nanoparticles/base fluid), and R_V 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, Tzou^10,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,¹³ the absence of a cross effect source in the buoyancy terms of their equation (10) and the presence of a complex dimensionless group R_m are observed. In this work, only thermal and volume fraction Rayleigh numbers are identified.

Terms –n₀R_TT k + R_T n₀Tn k in equation (8) come from the products in the density state equation (6) and are analogous to those in the model equations of Tzou.^10,11 These terms may also be the source of changes in the instability of the fluid which introduce interesting differences with the familiar double diffusion problem.²² Next, as the nanofluid layer stability is considered, a set of small perturbations shall be introduced into the dimensionless governing equations (7)–(10). These perturbations are

\begin{matrix} u = u_{1} (x, y, z, t) T = T_{b} (z) + T_{1} (x, y, z, t) n = n_{b} (z) \\ + n_{1} (x, y, z, t) \\ p = p_{b} (z) + p_{1} (x, y, z, t) \end{matrix}

(11)

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.¹³ In dimensional form this is

\begin{matrix} J \cdot k = 0 at z^{*} = 0, h \\ (- ρ_{p} D_{B} \nabla n^{*} - ρ_{p} \frac{D_{T}}{T_{T}^{*}} \nabla T^{*}) \cdot k = 0 at z^{*} = 0, h \end{matrix}

(12)

According to Nield and Kuznetsov¹³ 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

0 = \frac{d n_{b}}{dz} + N_{A} \frac{d T_{b}}{dz} at z = 0, 1 respectively

(13)

T_{b} = 0 at z = 0, 1 respectively

(14)

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 p_b, T_b, and n_b

0 = - \frac{d p_{b}}{dz} - (n_{0} - 1) R_{T} T_{b} - R_{V} n_{b} - (\frac{h^{3} g}{κ ν} + R_{V}) + R_{T} n_{0} T_{b} n_{b}

(15)

0 = \frac{d^{2} T_{b}}{d z^{2}} + \frac{N_{B}}{Le} \frac{d n_{b}}{dz} \frac{d T_{b}}{dz} + \frac{N_{A} N_{B}}{Le} {(\frac{d T_{b}}{dz})}^{2}

(16)

0 = \frac{1}{Le} \frac{d^{2} n_{b}}{d z^{2}} + \frac{N_{A}}{Le} \frac{d^{2} T_{b}}{d z^{2}}

(17)

Most attention is given to the temperature profile T_b and to the volumetric fraction profile n_b. Because T_b and n_b are independent of p_b, these are first calculated from equations (16) and (17) and later substituted into equation (15) from which the solution for p_b is straight. Since pressure and velocity fields are later split, p_b 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 T_b and n_b may exist. Also, the linear stability of two different basic states, for temperature and volumetric fraction, have been analyzed by Nield and Kuznetsov^12,13 and Tzou,^10,11 but there is still another set of possible basic states for this system. These are

n_{b 1} = - \frac{Le C_{2} z}{\sqrt{Le C_{2} N_{B}}} + \frac{\exp [\frac{(C_{3} + z) \sqrt{Le C_{2} N_{B}}}{Le}]}{C_{2} N_{B}} + C_{4}

(18)

n_{b 1} = \frac{Le C_{2} z}{\sqrt{Le C_{2} N_{B}}} + \frac{\exp [- \frac{(C_{3} + z) \sqrt{Le C_{2} N_{B}}}{Le}]}{C_{2} N_{B}} + C_{4}

(19)

T_{b} = \frac{Le}{N_{A} N_{B}} \ln (\frac{d^{2} n_{ba}}{d z^{2}}) - \frac{n_{ba}}{N_{A}} + C_{1}

(20)

where C₁ through C₄ 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 Kuznetsov^12,13 and Tzou,^10,11 for a TiO₂ particle–based nanofluid (see Table 1), the volume fraction profile equations (18)–(19) showed non-physical behavior. Further comparisons for the temperature profiles^10,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.

Table 1.

Data for TiO₂ particles–based nanofluids.^8,19,23,24

Parameter	Value	Parameter	Value
κ	1.52 × 10^–7 m²/s	Δn^*	0.38%^b
k_B	1.38064852 × 10^–23	$\bar{n_{b}}$	0.57%^c
$\bar{T}$	300 K	ΔT^*	2.3 K^d
d_p	14 nm	$T_{T}^{}$*	297.7 K
µ	8.33 × 10^–4 Pa s	β_B	8.28 × 10^–2
k	0.61497 W/m K	D_T	3.95 × 10^–8 m²/s
k_p	≈0.59–1.1 W/m K	D_B	3.77 × 10^–11 m²/s
ρ_p(c_p)_p	2.94 MJ/m³ K	Le	4022.88
ρ ₀(c_p)₀	4.05 MJ/m³ K	N_A	21.30
ρ_nf	≈9.96 × 10² kg/m^3a	N_B	27.54 × 10^–2

k_B is the Boltzmann constant, $\bar{T}$ is the average nanofluid temperature, d_p is the particle diameter, k_p is the particle thermal conductivity, ρ_nf is the nanofluid density, Δn* is a volumetric fraction difference, and ΔT* is a temperature difference.

Parameters β_B, D_T, and D_B were calculated through equations in the work by Buongiorno.⁸

Water data were taken.

Δn* = 57–19 was taken for demonstration based on the data of Wen and Ding.¹⁹

$\bar{n_{b}}$ was taken for demonstration based on the data of Wen and Ding.¹⁹

This is the highest temperature value for 0.57% used in Wen and Ding.¹⁹

Figure 2 shows plots of n_b of previous works for the case of TiO₂ 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¹³ because of its simplicity and need of further comparison in view of the present extension to their results. Then, T_b and n_b are expressed as

T_{b} = 1 - z

(21)

n_{b} = (\bar{n_{b}} - \frac{N_{A}}{2}) + N_{A} z

(22)

where $\bar{n_{b}}$ is the average volume fraction which was omitted in previous studies. In fact, $\bar{n_{b}}$ had to be introduced since an integration constant in the solution to n_b could not be determined through the boundary conditions.

Figure 2.

Plot of the volumetric fraction basic states used by Nield and Kuznetsov¹³ and Tzou^10,11 for a TiO¹¹ nanofluid according to the data in Table 1. Notice that although both behave linearly, the boundary conditions change the basic state. The solid line corresponds to equation (22).¹³

The perturbation equations

The perturbed governing equations, for u₁, p₁, T₁, and n₁, follow after the set of perturbations given in equation (11) is substituted into equations (7)–(10). In these set of equations, pressure and velocity fields can split by operating twice $\nabla x \nabla x$ on the momentum balance equation which is a very convenient simplification. For the linear stability problem, nonlinear terms are omitted and the following model equations are obtained

\begin{matrix} P r^{- 1} \frac{\partial \nabla^{2} w_{1}}{\partial t} = \nabla^{4} w_{1} + R_{T} \nabla_{⊥}^{2} T_{1} \\ - R_{V} \nabla_{⊥}^{2} n_{1} + R_{T} n_{0} \nabla_{⊥}^{2} (n_{b} T_{1}) + \\ R_{T} n_{0} \nabla_{⊥}^{2} (T_{b} n_{1}) - R_{T} n_{0} \nabla_{⊥}^{2} T_{1} \end{matrix}

(23)

\frac{\partial T_{1}}{\partial t} + w_{1} \frac{d T_{b}}{dz} = \nabla^{2} T_{1} + \frac{N_{B}}{Le} (\frac{d n_{b}}{dz} \frac{\partial T_{1}}{\partial z} + \frac{d T_{b}}{dz} \frac{\partial n_{1}}{\partial z}) + \frac{2 N_{A} N_{B}}{Le} \frac{d T_{b}}{dz} \frac{\partial T_{1}}{\partial z}

(24)

\frac{\partial n_{1}}{\partial t} + w_{1} \frac{d n_{b}}{dz} = \frac{1}{Le} \nabla^{2} n_{1} + \frac{N_{A}}{Le} \nabla^{2} T_{1}

(25)

where w₁ 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, N_A, and N_B 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, w₁, T₁, and n₁ are rewritten as

(w_{1}, T_{1}, n_{1}) = [W (z), Θ (z), Φ (z)] \exp [σ t + i (a_{x} x + a_{y} y)]

(26)

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

\frac{i ω}{\Pr} (D^{2} - a^{2}) W = {(D^{2} - a^{2})}^{2} W - (1 - n_{0} + n_{0} n_{b}) R_{T} a^{2} Θ + (R_{V} - R_{T} n_{0} T_{b}) a^{2} Φ

(27)

i ω Θ + W \frac{d T_{b}}{dz} = (D^{2} - a^{2}) Θ + \frac{N_{B}}{Le} (\frac{d n_{b}}{dz} D Θ + \frac{d T_{b}}{dz} D Φ) + \frac{2 N_{A} N_{B}}{Le} \frac{d T_{b}}{dz} D Θ

(28)

i ω + W \frac{d n_{b}}{dz} = \frac{1}{Le} (D^{2} - a^{2}) Φ + \frac{N_{A}}{Le} (D^{2} - a^{2}) Θ

(29)

where D = d/dz and $a^{2} = a_{x}^{2} + a_{y}^{2}$ is the perturbations wavenumber. Further calculations of the stability of the fluid layer are concerned with the critical values of the Rayleigh number so that σ_R = 0 from now on. For the case of convection in a nanofluid layer confined between perfect thermal conducting walls, the boundary conditions are

W = DW = D Φ + N_{A} D Θ = Θ = 0 at z = 0, 1

(30)

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.¹³ 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.²⁵ Some Fortran routines from the book of Press et al.²⁵ 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

W = \sum_{p = 1}^{N} z^{1 + p} {(1 - z)}^{1 + p}, = \sum_{p = 1}^{N} z^{p} {(1 - z)}^{p}, = \sum_{p = 1}^{N} - N_{A} z^{p} {(1 - z)}^{p}

(31)

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

| \begin{matrix} - \frac{i ω}{\Pr} 〈 W_{i} L W_{j} 〉 + 〈 W_{i} L^{2} W_{j} 〉 - R_{T} a^{2} 〈 W_{i} (1 - n_{0} + n_{0} n_{b}) Θ_{j} 〉 a^{2} 〈 W_{i} (R_{V} - R_{T} n_{0} T_{b}) Φ_{j} 〉 \\ - \frac{d T_{b}}{dz} 〈 Θ_{i} W_{j} 〉 - i ω 〈 Θ_{i} Θ_{j} 〉 + 〈 Θ_{i} L Θ_{j} 〉 + A_{0} 〈 Θ_{i} D Θ_{j} 〉 \frac{N_{B}}{Le} \frac{d T_{b}}{dz} 〈 Θ_{i} D Φ_{j} 〉 \\ - \frac{d n_{b}}{dz} 〈 Φ_{i} W_{j} 〉 \frac{N_{A}}{Le} 〈 Φ_{i} L Θ_{j} 〉 - i ω 〈 Φ_{i} Φ_{j} 〉 + \frac{1}{Le} 〈 Φ_{i} L Φ_{j} 〉 \end{matrix} | = 0

(32)

where

L = D^{2} - a^{2}

A_{0} = \frac{N_{B}}{Le} \frac{d n_{b}}{dz} + 2 \frac{N_{A} N_{B}}{Le} \frac{d T_{b}}{dz}

〈 \dots 〉 = \int_{0}^{1} \dots dz

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

\begin{matrix} 0 = - [\frac{A_{2} N_{A}^{2}}{567, 000} + \frac{28}{27} \frac{(a^{2} - 336) (1 + Le) A_{1}}{a^{2} \Pr}] ω^{2} - A_{1} R_{T} - (1 + Le) A_{1} N_{A} R_{V} \\ + (\frac{Le N_{A}}{2} + 1 - \bar{n_{b}}) A_{1} R_{T} n_{0} + \frac{28}{27} \frac{A_{1} A_{2} (a^{2} + 10)}{a^{2}} - \frac{i ω^{3}}{567, 000} \frac{(a^{2} + 12) N_{A}^{2}}{\Pr} - i ω \\ [\frac{[R_{T} - (1 - \bar{n_{b}}) R_{T} n_{0} + R_{V} N_{A}] a^{2} N_{A}^{2}}{588, 000} - \frac{28}{27} \frac{A_{1} A_{2} (1 + Le)}{a^{2}} - \frac{28}{27} \frac{A_{1} (a^{2} + 12) (a^{2} + 10)}{a^{2} \Pr}] \end{matrix}

(33)

where

A_{1} = \frac{(a^{2} + 10) a^{2} N_{A}^{2}}{588, 000 Le}

(34)

A_{2} = a^{4} + 24 a^{2} + 504

(35)

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 $\bar{n_{b}}$ . However, there have been experimental reports indicating that particle volume fraction gradients may or may not be an important destabilizing mechanism just at the beginning of convective motions.^19,20 As expected and according to equation (33), both R_T and R_V compete to destabilize the fluid layer but, also, a number of physical parameters may influence this competition as well.

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

[1 - \frac{Le n_{0} N_{A}}{2} + (\bar{n_{b}} - 1) n_{0}] R_{T}^{s} + (Le + 1) R_{V}^{s} = \frac{28 (a^{4} + 24 a^{2} + 504) (a^{2} + 10)}{27 a^{2}}

(36)

where Brownian motions and thermophoretic effect can be identified by Le and N_A. Notice that the ratio of volumetric heat capacities N_B 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 a_c = 3.12 and

[1 - \frac{Le n_{0} N_{A}}{2} + (\bar{n_{b}} - 1) n_{0}] R_{T}^{s} + (Le + 1) R_{V}^{s} = 1750

(37)

In other words, since thermal Rayleigh number $R_{T}^{s}$ is a function of the wavenumber a and convective motions will set in at the smallest value of $R_{T}^{s}$ , then equation (36) is derived against a to find the minimum a which is the critical value. When critical wavenumber a_c is achieved, Rayleigh numbers are critical too as in equation (37). Notice that the coefficient of $R_{T}^{s}$ is negative, due to typical magnitudes of Le and N_A (see Table 1, for example), and in equation (37) and as a consequence $R_{V}^{s}$ is, for this case, larger than $R_{T}^{s}$ .

Oscillatory convection

In order to study the onset of oscillatory convection, R_T is isolated from equation (33). This can be written symbolically as

R_{T}^{o} = f_{1} + i ω f_{2}

(38)

where f₁ and f₂ depend on ω, a, R_V, Pr, Le, N_A, n₀, and $\bar{n_{b}}$ . Since ω ≠ 0 for oscillatory convection, it follows that f₂ = 0. Thus, the following equation for ω is obtained

\begin{matrix} 0 = [112 A_{2} A_{3} ⪻ + (a^{2} + 12) (a^{2} + 10) (56 N_{A} n_{0} + 112 A_{3})] \\ L e^{2} ω^{2} + (54 N_{A} n_{0} + 108 A_{3}) (a^{2} + 10) a^{2} N_{A} L e^{2} \Pr R_{V} \\ - [56 (1 + Le) n_{0} N_{A} Le - 112 A_{3}] A_{2} {(a^{2} + 10)}^{2} \Pr \\ - (a^{2} + 12) {(a^{2} + 10)}^{3} (56 N_{A} n_{0} Le - 112 A_{3}) \end{matrix}

(39)

where

A_{3} = (\bar{n_{b}} - 1) n_{0} + 1

(40)

Notice that A₁, A₂, A₃ > 0 provided that the wavenumber is nonzero and $\bar{n_{b}} < 1$ too. In this way, equation (39) implies that there is a possible value for the frequency of oscillation ω if the following condition holds

\begin{matrix} - (54 N_{A} n_{0} + 108 A_{3}) (a^{2} + 10) a^{2} N_{A} L e^{2} \Pr R_{V} \\ + [56 (1 + Le) n_{0} N_{A} Le - 112 A_{3}] A_{3} \\ {(a^{2} + 10)}^{2} \Pr + (a^{2} + 12) {(a^{2} + 10)}^{3} \\ (56 N_{A} n_{0} Le - 112 A_{3}) > 0 \end{matrix}

(41)

Next, the thermal Rayleigh number for oscillatory convection $R_{T}^{o}$ given in equation (38) reduces to

\begin{matrix} R_{T}^{o} = 2 [- 56 (a^{2} + 12) A_{3} L e^{2} ω^{4} - 54 a^{2} A_{3} N_{A} L e^{2} ⪻ R_{V} ω^{2} \\ + 28 (a^{2} + 10) (n_{0} N_{A} + 2 A_{3}) A_{2} L e^{2} \Pr ω^{2} \\ + 28 (a^{2} + 12) {(a^{2} + 10)}^{2} ((1 + Le) n_{0} N_{A} Le - 2 A_{3}) ω^{2} \\ + (27 n_{0} N_{A} Le - 54 A_{3}) (1 + Le) {(a^{2} + 10)}^{2} a^{2} N_{A} \Pr R_{V} \\ - 28 (n_{0} N_{A} Le - 2 A_{3}) {(a^{2} + 10)}^{3} A_{2} \Pr] \\ {[27 a^{2} \Pr (4 A_{3}^{2} L e^{2} ω^{2} + {(2 A_{3} - n_{0} N_{A} Le)}^{2} {(a^{2} + 10)}^{2})]}^{- 1} \end{matrix}

(42)

where ω² and ω⁴ 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 $J_{B}^{*} = - ρ_{p} D_{B} \nabla n^{*}$ in order to eliminate thermophoretic influence from the model equations. As a consequence, and after following the same above nondimensionalization and perturbation process, the governing equations for temperature and particle volume fraction become

i ω Θ + W \frac{d T_{b}}{dz} = (D^{2} - a^{2}) Θ + \frac{N_{B}}{Le} \frac{d T_{b}}{dz} D Φ_{B}

(43)

i ω Φ_{B} = \frac{1}{Le} (D^{2} - a^{2}) Φ_{B}

(44)

in which some further simplifications were made due to the new basic state (n_b)_B. Since the thermophoretic effect was eliminated, the boundary conditions for (n_b)_B were modified too. This is

\frac{d {(n_{b})}_{B}}{dz} = 0 at z = 0, 1

(45)

{(n_{b})}_{B} = \bar{n_{b}}

(46)

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 $J_{B}^{*}$ . Then, following the Galerkin method, a new trial function for Φ_B is chosen,²⁷ in view of the new boundary conditions

Φ_{B} = \sum_{p = 1}^{N} z^{p - 1}

(47)

At the first order of approximation, with the Galerkin method applied on the eigenvalue problem equations (27), (43), and (44), it is obtained

\begin{matrix} 0 = \frac{ω_{B}^{2}}{18, 900} (A_{2} + \frac{120 Le + (1 + Le) a^{4} + 2 (11 Le + 6) a^{2}}{LePr}) + \frac{A_{3} a^{4} R_{TB}}{19, 600 Le} \\ - \frac{(a^{2} + 10) a^{2} A_{2}}{18, 900 Le} + \frac{(a^{2} + 12) {i ω}_{B}^{3}}{18, 900 \Pr} \\ - [- \frac{A_{3} a^{2} R_{TB}}{19, 600} + \frac{A_{2} (a^{2} Le + a^{2} + 10 Le)}{18, 900 Le} + \frac{(a^{2} + 10) (a^{2} + 12) a^{2}}{18, 900 LePr}] i ω_{B} \end{matrix}

(48)

The absence of R_V in equation (48) is noticeable, but the effect of particle volume fraction still remains through $\bar{n_{b}}$ and n₀. For the steady convection case, R_TB can be easily determined to be

R_{TB}^{s} = \frac{28 (a^{2} + 10) A_{2}}{27 A_{3} a^{2}}

(49)

and the critical wavenumber is a_cB = 3.12, again. On the other hand, in the absence of thermophoresis oscillatory convection vanishes. ω = 0 was expected because of the only iR_Tω 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 TiO₂ 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 $\bar{n_{b}}$ and reference value n₀ for the particle volume fraction may become important for the fluid layer stability. Notice that $\bar{n_{b}}$ and n₀ are a direct contribution of the volumetric fraction basic state through the buoyancy terms in the momentum balance equation (23). These contributions were early envisaged by the way in which the density equation state equation (6) is formulated, and as a result three buoyancy terms in equation (23) make an interesting difference in comparison to the classical double diffusive problem.

Steady and oscillatory convection

The present analysis shows that steady and oscillatory convection is possible, for positive R_T and R_V. It has to be mentioned that considered values for $\bar{n_{b}}$ and n₀ were small so that stability may be changed for other nanofluids with larger values.

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 a_c = 3.12 could be determined explicitly which helps to visualize the competition between thermal and volume fraction destabilizing forces. The ratio of diffusion coefficients N_A 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.

Figure 3.

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 $R_{T}^{s}$ there is a corresponding critical $R_{V}^{s}$ on the solid line. On this steady curve of criticality, R_V is always smaller than R_T meaning that volumetric fraction convection is more likely to occur than thermal convection. For example, in the experimental work of Wen and Ding,¹⁹ it is reported that particle gradients may trigger secondary fluid motions. Near the locus of criticality, the linear relationship between R_V and R_T is changed as mentioned in Figure 3 since other parallel lines could be plotted for different wavenumbers. According to the present investigation, an interesting fact of steady convection is that the critical wavenumber will never be changed despite the nanofluid properties since a_c = 3.12. However, according to equation (37) plotted in Figure 3, different nanofluid properties would only change the slope of the shown line of criticality. The slope S_s can be easily obtained from equation (37) as

S_{s} = \frac{(Le + 1) N_{A}}{1 - \frac{Le n_{0} N_{A}}{2} + (\bar{n_{b}} - 1) n_{0}}

(50)

The lower values of the Rayleigh particle volume fraction R_V in comparison with R_T may be due to the large density difference between nanoparticles and fluid as pointed out by Wen and Ding.¹⁹ Despite the possible sign change in the term $(\bar{n_{b}} - 1) n_{0}$ , for small particle volume fractions, this may not represent a big contribution due to the commonly larger magnitudes of Le and N_A.

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 $R_{T}^{o}$ expression given in equation (42) and are graphically shown in Figures 4 –6.

Figure 4.

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

Figure 5.

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

Figure 6.

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 $R_{V}^{o}$ . This result would imply that oscillatory convection would always depend on the volume fraction Rayleigh number. The found relationship between a and $R_{V}^{o}$ is plotted in Figure 5. Next, with a and $R_{V}^{o}$ , the thermal Rayleigh number $R_{T}^{o}$ and the frequency of oscillation ω can be easily calculated from equations (42) and (39), respectively.

For the TiO₂ nanoparticle fluid, it was found that in oscillatory convection the curve $R_{T}^{o} - R_{V}^{o}$ is nonlinear. Figure 4 also shows that thermal gradient effect is restricted since, approximately, below a = 3.04, R_T becomes negative and in Figure 6 frequencies of oscillation are restricted due to nonphysical roots of equation (39), as well.

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 $J^{*} = J_{B}^{*} = - ρ_{p} D_{B} \nabla n^{*}$ in the set of governing equations.

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 $\bar{n_{b}}$ and n₀. Besides, the absence of thermophoresis works against oscillatory convection allowing steady convective motions in the fluid. There is also a restriction for $R_{TB}^{s}$ since $\bar{n_{b}} < - 1 / n_{0} + 1$ should be satisfied. Otherwise, the thermal gradient would be negative. This last inequality could be a limitation of the model because smaller particle volume fraction values reduce the problem to that of classical Rayleigh convection which can be easily demonstrated from equation (49) after taking $\bar{n_{b}}, n_{0} ≪ 1$ . That is, a departure from the familiar R_T = 1707 for solely thermal convection for small particle volume fraction is expected.

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 TiO₂ nanoparticle fluid.

For TiO₂ nanoparticle-based fluids steady convection sets in at critical wavenumber a_c = 3.12, but thermal R_T and nanoparticle volume fraction R_V 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 a_c = 3.12, the thermal Rayleigh number R_T is always larger than R_V. 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 TiO₂ nanoparticle volume fractions would change, in opposite direction, the competence with thermal gradients to destabilize the fluid layer. The R_V 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 R_V is smaller than that of R_T 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, R_V decreases and R_T 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 R_T 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.

ORCID iD

Ildebrando Pérez-Reyes

References

Sheikholeslami

Ganji

. Applications of nanofluid for heat transfer enhancement (Micro and nano technologies). New York: Elsevier, 2017.

Mahian

Kianifar

Kalogirou

et al . A review of the applications of nanofluids in solar energy. Int J Heat Mass Transfer 2013; 57: 582–594.

Sivashanmugam

. Chapter 14: application of nanofluids in heat transfer. In: Kazi

(ed.) An overview of heat transfer phenomena. Rijeka: InTech, 2012, pp. 411–440.

Wong

De Leon

. Applications of nanofluids: current and future. Adv Mech Engineering 2010; 2: 519659.

Sheikholeslami

Rashidi

Ganji

. Effect of nonuniform magnetic field on forced convection heat transfer of Fe₃O₄–water nanofluid. Comp Methods App Mech Eng 2015; 294: 299–312.

Freidoonimehr

Rashidi

Mahmud

. Unsteady MHD free convective flow past a permeable stretching vertical surface in a nano-fluid. Int J Thermal Sci 2015; 87: 136–145.

Garoosi

Jahanshaloo

Rashidi

et al . Numerical simulation of natural convection of the nanofluid in heat exchangers using a Buongiorno model. Appl Mathemat Comput 2015; 254: 183–203.

Buongiorno

. Convective transport in nanofluids. J Heat Transfer 2006; 128: 240.

Rashidi

Abelman

Mehr

. Entropy generation in steady MHD flow due to a rotating porous disk in a nanofluid. Int J Heat Mass Transfer 2013; 62: 515–525.

10.

Tzou

. Thermal instability of nanofluids in natural convection. Int J Heat Mass Transfer 2008; 51: 2967–2979.

11.

Tzou

. Instability of nanofluids in natural convection. J Heat Transfer 2008; 130: 072401.

12.

Nield

Kuznetsov

. The onset of convection in a horizontal nanofluid layer of finite depth. Eur J Mech B 2010; 29: 217–223.

13.

Nield

Kuznetsov

. The onset of convection in a horizontal nanofluid layer of finite depth: a revised model. Int J Heat Mass Transfer 2014; 77: 915–918.

14.

Das

Choi

SUS

Patel

. Heat transfer in nanofluids—a review. Heat Transfer Eng 2006; 27: 3–19.

15.

Rao

Srivastava

. Interferometric study of natural convection in a differentially-heated cavity with Al₂O₃–water based dilute nanofluids. Int J Heat Mass Transfer 2016; 92: 1128–1142.

16.

Yadav

Lee

. Onset of convection in a nanofluid layer confined within a Hele-Shaw cell. J App Fluid Mech 2016; 9: 519–527.

17.

Liu

Chang

et al . Natural convection heat transfer of alumina-water nanofluid in vertical square enclosures: an experimental study. Int J Thermal Sci 2010; 49: 1345–1353.

18.

Peterson

. Experimental studies of natural convection heat transfer of Al₂O₃/DI water nanoparticle suspensions (nanofluids). Adv Mech Eng 2010; 2: 742739.

19.

Wen

Ding

. Formulation of nanofluids for natural convective heat transfer applications. Int J Heat Fluid Flow 2005; 26: 855–864.

20.

Putra

Roetzel

Das

. Natural convection of nanofluids. Heat Mass Transfer 2003; 39: 775–784.

21.

Haddad

Abu-Nada

Oztop

et al . Natural convection in nanofluids: are the thermophoresis and Brownian motion effects significant in nanofluid heat transfer enhancement? Int J Thermal Sci 2012; 57: 152–162.

22.

Turner

. Buoyancy effects in fluids. Cambridge: Cambridge University Press, 1973.

23.

Abdulwahab

Thahab

Dhiaa

. Experimental study of thermophysical properties of TiO₂ nanofluid. Iraqi J Chem Petrol Eng 2016; 17: 1–6.

24.

Saeedian

Mahjour-Shafiei

Shojaee

et al . Specific heat capacity of TiO₂ nanoparticles. J Comp Theor Nanosci 2012; 9: 616–620.

25.

Press

Teukolsky

Vetteling

et al . Numerical recipes: the art of scientific computing. 3rd ed. Cambridge: Cambridge University Press, 2007.

26.

Chandrasekhar

. Hydrodynamic and hydromagnetic stability. New York: Dover Publications, 1981.

27.

Finlayson

. The method of weighted residuals and variational principles (Mathematics in Science and Engineering). New York: Academic Press, 1972.