## Abstract

A. E. Green, FRS and P. M. Naghdi developed a new theory of continuum mechanics based on an entropy identity rather than an entropy inequality. In particular, within the framework of this theory, they developed a new set of equations to describe viscous flow. The new theory additionally involves vorticity and spin of vorticity. We here develop the theory of Green and Naghdi to be applicable to thermal convection in a fluid in which is suspended a collection of minute metallic-like particles. Thus, we develop a non-Newtonian theory we believe capable of describing a nanofluid. Numerical results are presented for copper oxide or aluminium oxide particles in water or in ethylene glycol. Such combinations are used in real nanofluid suspensions.

## 1. Introduction

Green & Naghdi (1991, 1995*a*–*c*, 1996) develop an analysis in a rational way to produce fully consistent theories of continuous media that incorporate thermal pulse transmission in a very logical manner. In particular, Green & Naghdi (1996) developed such a theory to describe nonlinear fluid behaviour.

In the papers cited above, Green & Naghdi developed a theory for describing the behaviour of a continuous body that relies on an entropy balance law rather than an entropy inequality. Their thermodynamics introduces a quantity *T*, which is the ‘empirical’ temperature, and a ‘thermal displacement variable’,
1.1
The paper of Green & Naghdi (1991) developed the thermodynamics of a rigid heat conductor employing their thermal displacement variable, *α*, and their balance of entropy equality. This very novel approach led to theories in which heat may propagate as a wave and, in particular, they developed the ideas of type II and type III rigid heat conductors, the classical parabolic theory being classified as type I. Type II essentially is heat transport without energy dissipation, whereas type III still allows finite speed of propagation, but incorporates dissipation. Extensions of type II and type III theories to continuous bodies that move have been highly successful, particularly in the context of thermoelasticity (e.g. Chandrasekharaiah 1998; Hetnarsky & Ignaczak 1999). The theory of Green & Naghdi (1996) develops a novel theory for fluids that involves vorticity and spin of vorticity. This introduces higher spatial gradients into the equations than those of the Navier–Stokes theory and so is likely to be relevant where non-Newtonian fluid behaviour is expected. They allow two temperatures and are motivated by turbulence and the work of Marshall & Naghdi (1989*a*,*b*). However, Green & Naghdi (1996) point out that their theory is of relevance in its own right. It is the goal of this paper to modify the theory of Green & Naghdi (1996) to be applicable to describe a nanofluid suspension. Such suspensions are very important in the current heat transfer industry (e.g. Maïga *et al.* 2005; Vadasz *et al.* 2005; Vadasz 2006; Kulkarni *et al.* 2008). To incorporate the suspension and solvent thermal characteristics, we adapt a non-thermal equilibrium idea originally proposed for flow in porous media where the solid and fluid phases may have different temperatures. Specific details of this idea in thermal convection in a porous medium may be found in Nield & Bejan (2006, p. 31), Nield (2002), Nield & Kuznetsov (1999), Banu & Rees (2002) and Malashetty *et al.* (2005), and in the nonlinear case in Straughan (2006).

Current theories for nanofluid behaviour (see Savino & Paterna 2008; Tzou 2008; Kuznetsov & Nield 2010, in press) all assume the fluid to be Newtonian. However, there is clear evidence in the literature (e.g. Kwak & Kim 2005) that a nanofluid suspension does not behave like a Newtonian (linearly viscous) fluid. Kwak & Kim (2005) demonstrate for a CuO suspension in ethylene glycol that the zero shear viscosity varies highly nonlinearly as a function of particle volume fraction and the same suspension exhibits very strong shear thinning. Thus, we believe there is a need for a theory where the fluid is non-Newtonian. By non-Newtonian we do not mean viscoelastic. A CuO nanofluid suspension contains particles of the shape of a prolate spheroid of aspect ratio 3 (see Kwak & Kim 2005), and such molecular fluids are well known to display non-Newtonian behaviour (e.g. Travis *et al.* 1997). For example, in Poiseuille flow, a flattened velocity profile is observed rather than the parabolic one for a Newtonian fluid. By including vorticity and spin of vorticity in the theory, Green & Naghdi (1996) achieve such a flattened profile. Non-Newtonian behaviour in fluid suspensions has been known for many years (see Frohlich & Sack 1946).

## 2. Basic equations

The basic equations of Green & Naghdi (1996) are the balance of linear momentum, balance of mass and balances of entropy for the temperatures *θ*_{H} and *θ*_{T}, namely,
2.1
Throughout, standard indicial notation is used with a repeated roman index denoting summation from 1 to 3, a superposed dot denotes the material derivative *d*/*dt*=∂/∂*t*+*v*_{i}∂/∂*x*_{i} and Δ is the Laplacian. In equations (2.1), *ρ*,*v*_{i},*b*_{i} and *p* are the density, velocity, body force and pressure. The density will be taken to be constant everywhere except in the body force term, where it will be allowed to depend on the temperatures *θ*_{H},*θ*_{T}. The coefficients *μ* and *μ*_{1} are the kinematic viscosity of the fluid and another constant reflecting the geometry of the particles and the interaction with the fluid. In fact, in the appendices to Bleustein & Green (1967) and Green & Rivlin (1964), they derive an expression for the kinetic energy of a system of particles as a function of the velocity of the centroid and the derivative of this velocity. While neither the article of Green & Rivlin (1964) nor that of Bleustein & Green (1967) has a direct bearing on the fluid theory of Green & Naghdi (1996), their procedure leads to a kinetic energy that can be equated to the kinetic energy of the current fluid. The relation between the dipolar inertia coefficient *ρμ*_{1}/*μ* and the coefficient 2*μ*_{1} is explored for both the dipolar fluid of Bleustein & Green (1967) and the fluid of Green & Naghdi (1996) by Jordan & Puri (2002), Quintanilla & Straughan (2005) and Puri & Jordan (2006). The quantities *η*,*s*,*ξ* and *p*_{i} (with sub/superscripts *H*, *T*) are entropy, external supply of entropy, intrinsic entropy supply (which is determined by constitutive theory and thermodynamics) and the entropy flux vectors. In our work, we associate the temperatures *θ*_{H} and *θ*_{T} with the temperature of the fluid and of the solid suspension particles, respectively. Thus, we identify *θ*_{H}=*θ*_{F} and *θ*_{T}=*θ*_{S}, where *F* and *S* signify fluid and solid, respectively. The same convention is used for *η*_{H}, etc. In fact, the identification of *θ*_{F} and *θ*_{S} has already been used in nanofluid theory by Kuznetsov & Nield (2010, in press), who developed a theory for a nanofluid suspension in a porous medium. However, their theory is not based on Green–Naghdi thermodynamics; rather, it is based on a linear viscous structure.

Green & Naghdi (1996) assume that the Helmhotlz free energy function *ψ* has the form
2.2
with *c*_{H},*c*_{T} positive constants, and the entropies and entropy fluxes have the form
2.3
and
2.4
for positive constants *κ*_{H},*κ*_{T},*θ*_{0}, with *θ*_{0} being some reference temperature. The intrinsic entropy supply functions have the form
2.5
and
2.6
Here, *D*_{ij}=(*v*_{i,j}+*v*_{j,i})/2, *P*_{ij}=−Δ*v*_{i,j} and *ϕ* in their theory is constant. It is very important to note, however, that Green & Naghdi (1996) observe that for some purposes *ϕ* could depend on temperatures.

By substituting equations (2.2)–(2.6) into equations (2.1)_{3,4}, Green & Naghdi (1996) derive their entropy balance eqns (62) and (63). Here we do not adopt equations (2.4)–(2.6) with *θ*_{0}. Indeed, this does not reduce to Fourier theory for the heat flux. Instead, we argue *θ*_{0} should be replaced by *θ*_{H} or *θ*_{T} as appropriate. If we now replace *H* by *F* and *T* by *S*, then we assume, instead of equations (2.4)–(2.6), the equations
2.7
2.8
2.9
It is important to note that we also assume
2.10
for *h*>0 constant. (We employ *k*_{F} and *k*_{S} in equations (2.7)–(2.9) because these are thermal conductivities and *k* is a more usual notation.)

With equations (2.7)–(2.10), one may then show that equations (2.1)_{3,4} reduce to
2.11
and
2.12
In the momentum equation (2.1)_{1}, we assume that the body force term is gravity with *ρ* depending on *θ*_{F},*θ*_{S} so that
2.13
where *g* is gravity, **k**=(0,0,1), *α*_{F}, *α*_{S} are the thermal expansion coefficients of the fluid and solid, respectively, and are reference (constant) temperatures, *ρ*_{0} being a constant.

The momentum and conservation of mass equations then become
2.14
and
2.15
where *ν*=*μ*/*ρ*_{0}, *ν*_{1}=*μ*_{1}/*ρ*_{0} and *p* absorbs the constant terms from equation (2.13).

Our basic equations governing the non-isothermal fluid flow are then equations (2.14), (2.15), (2.11) and (2.12).

One may then ask whether the solid particles really do contribute to the buoyancy in equation (2.13). Thus, to be consistent with other approaches to nanofluids, we also consider the equation
2.16
where *ϕ* is the volume fraction of particles in the suspension. The modifications that arise for thermal convection by employing equation (2.16) rather than equation (2.13) are described below. (In fact, very little difference is noticed in the numerical results.)

## 3. Thermal convection

Suppose now the fluid occupies the horizontal layer *z*∈(0,*d*) with gravity acting downward. Then, equations (2.14), (2.15), (2.11) and (2.12) hold on the domain . The boundaries are assumed fixed with the temperatures maintained at constants so that
3.1
with *T*_{L}>*T*_{U}.

Equations (2.14), (2.15), (2.11) and (2.12) then admit the steady (conduction) solution
3.2
where *β* is the temperature gradient,
3.3

The pressure is then found from equation (2.14).

To investigate the instability of solution (3.2), we let (*u*_{i},*π*,*θ*_{F},*θ*_{S}) be perturbations to the basic state and put , , in equations (2.14), (2.15), (2.11) and (2.12). In what follows, *θ*_{F} and *θ*_{S} are always the perturbation variables. One finds the following equations for (*u*_{i},*π*,*θ*_{F},*θ*_{S}):
3.4
where *d*_{ij}=(*u*_{i,j}+*u*_{j,i})/2, *p*_{ij}=−Δ*u*_{i,j} and *w*=*u*_{3}.

Equations (3.4) are now non-dimensionalized with the scalings
The Rayleigh number *Ra*=*R*^{2} is defined by
Equations (3.4) in non-dimensional form become, dropping the *’s on *x*_{i} and *t*,
3.5

One may show that, if one adopts the density relation (2.16), then, with *ρ*_{0} replaced by *ρ*_{N}=*ρ*_{S}*ϕ*+*ρ*_{F}(1−*ϕ*), we again arrive at equations (3.5), but Ra must now be replaced by Ra/(1−*ϕ*). The definitions of the nanofluid Rayleigh number given in §5 and the subsequent values in §6 must then be modified by the value (1−*ϕ*)^{−1}.

## 4. Linear instability

Equations (3.5) are now linearized, and then we put
We take curlcurl of what remains from equation (3.5)_{1} in order to eliminate *π*. In this manner, we are led to solve the system of equations
4.1
where Δ*=∂^{2}/∂*x*^{2}+∂^{2}/∂*y*^{2} is the horizontal Laplacian. This is an eigenvalue problem for *σ* on the spatial domain .

The boundary conditions are those of two fixed surfaces at constant temperatures, so we have 4.2

However, owing to the higher order terms in the Green–Naghdi momentum equation, we require an extra boundary condition. For this, we follow Green & Naghdi (1996, eqn (74)). To do this, we introduce their tensor , where *p*_{1} is an arbitrary scalar and *N*_{ij}=*ϵ*_{irs}*u*_{s,rj}, i.e. ∂/∂*x*_{j} of (*curl* **u**)_{i}. Then, the couple at the walls is , where **n**=(0,0,1) at *z*=1 and **n**=(0,0,−1) at *z*=0. The couple in the horizontal directions *x*,*y* is supposed zero. Thus,
**e**^{1},**e**^{2} being standard basis vectors. If we write **u**=(*u*,*v*,*w*), this leads us to the equations
4.3
But in the whole fluid domain *u*_{x}+*v*_{y}+*w*_{z}=0, so
4.4
also on *z*=0,1. We may differentiate equation (4.3) in the *x* and *y* directions, and so
4.5
Equations (4.5) and (4.3) together yield the boundary condition
4.6
We now follow the standard procedure and put *w*=*W*(*z*)*f*(*x*,*y*), *θ*_{F}=*Θ*_{F}(*z*)*f*(*x*,*y*) and *θ*_{S}=*Θ*_{S}(*z*)*f*(*x*,*y*), where *f* is a plane-tiling function such that
*a* being a wavenumber. Explicit forms for *f* are given in, for example, Chandrasekhar (1981, pp. 43–52).

If we denote by *D*=*d*/*dz*, then equations (4.1) reduce to
4.7
where *α*=*α*_{S}/*α*_{F}, *k*=*k*_{F}/*k*_{S} and *κ*=*κ*_{F}/*κ*_{S}. (If and *H*=0, equations (4.7) are formally equivalent to the equations for double-diffusive convection in a fluid layer heated and salted from below.) System (4.7) is to be solved for the eigenvalues *σ* subject to the boundary conditions
4.8

System (4.7) and (4.8) is solved numerically by a *D*^{2} Chebyshev tau method (see Dongarra *et al.* 1996), and so we define *χ*=(*D*^{2}−*a*^{2})*W*, *Φ*=(*D*^{2}−*a*^{2})*χ* and solve the matrix system *A***x**=*σB***x**, where the matrices *A* and *B* are given by
and
with the boundary conditions (4.8) added to the last two rows of each of the blocks of *A*. (The last two rows of each of the blocks of *B* are set equal to zero.) The vector **x** is given by
where *W*_{i} are the Chebyshev coefficients in the expression *T*_{i} being the Chebyshev polynomials. A similar interpretation holds for .

## 5. Numerical results

We now report numerical results for the instability analysis presented in §4. We employed four sets of nanofluid suspensions. Namely, those where the fluid was either water or ethylene glycol and where the particles are CuO or Al_{2}O_{3}. Values of *α*=*α*_{S}/*α*_{F}, *k*=*k*_{F}/*k*_{S}, *κ*=*k*_{F}*c*_{S}/*k*_{S}*c*_{F} and *Pr* were obtained from Chandrasekhar (1981, p. 66), Dow (2009), Accuratus (2009) and Engineering Toolbox (2009). The values used were

— H

_{2}O–CuO,*α*=0.0797,*k*=1.494×10^{−3},*Pr*=7.02,*κ*=1.373×10^{−4};— H

_{2}O–Al_{2}O_{3},*α*=0.04,*k*=1.712×10^{−2},*Pr*=7.02,*κ*=3.596×10^{−3};— ethylene glycol–CuO,

*α*=0.02538,*k*=6.426×10^{−4},*Pr*=159.25,*κ*=1.01889×10^{−4};— ethylene glycol–Al

_{2}O_{3},*α*=0.0129,*k*=7.3629×10^{−3},*Pr*=159.25,*κ*=2.668×10^{−3}.

The numerical routine calculates the minimum value of *a* where instability will commence, i.e. that value of Ra for which *σ*_{r}=0, where *σ*=*σ*_{r}+*iσ*_{1}. For all of the parameter values we investigated, *σ* was found to be real at the instability transition. Thus, instability is by stationary convection and overstability is not witnessed (at least for the parameter values we addressed).

We point out that the Rayleigh number we have defined is based on the fluid properties only. Since it is well known that the thermal conductivity of a nanofluid may be considerably higher than that of the solvent (see Xuan & Roetzel 2000; Putra *et al.* 2003; Xuan *et al.* 2003; Kwak & Kim 2005; Hwang *et al.* 2007; Kim *et al.* 2007; Wong & Kurma 2008; Masoumi *et al.* 2009), and the recent work of Kwak & Kim (2005) and Masoumi *et al.* (2009) shows that the effective viscosity of the nanofluid may also be significantly different from that of the base fluid, we wish to give consideration to a Rayleigh number which might more accurately reflect the effective properties of the nanofluid suspension itself. In this regard, we use the thermal conductivity expression in Kwak & Kim (2005, eqn (6)), where
5.1
where *k*_{eff} is the effective thermal conductivity, is the ratio of thermal conductivities of the solid particles and the base fluid and *ϕ* is the volume fraction of the particles. For the viscosity variation, there are several possibilities. Kwak & Kim (2005) quote that the zero shear viscosity is given by
where *p* is the aspect ratio. The intrinsic viscosity though has the form
For an aspect ratio of 3, we have and [*η*]_{0}=2.18. For the onset of thermal convection instability, we do not expect strong shear flow and so we might expect a viscosity correction like 2.18/0.39=5.5897. The thermal conductivity correction depends on the volume concentration *ϕ*. For all four combinations of fluid/particle we considered, we find a correction factor of 1.06 for a 4 per cent concentration and 1.12 for a 2 per cent concentration. Thus, we define a nanofluid Rayleigh number, (where *J* denotes the percentage by volume of particles), by
or
Alternatively, if we consider the viscosity correction (7) quoted by Hwang *et al.* (2007),
then, at a 4 per cent concentration, *μ*_{eff}=3.41864*μ*_{F}. Thus, for a 4 per cent concentration, we might consider another nanofluid Rayleigh number, *Ra*_{1}, defined by

These values are tabulated along with Ra in table 1.

We present numerical results only for the H_{2}O–CuO case, the other combinations being similar. We do not have accurate values for *H* nor for . Therefore, we estimate these and present results for various values.

## 6. Conclusions

We have presented a modification of the fundamental theory of viscous flow presented by Green & Naghdi (1996). We have then derived a detailed analysis of linearized instability in the problem of thermal convection in a layer of fluid heated from below. An interpretation of the model has been made in the context of a nanofluid suspension where the different temperatures are ascribed to the fluid and solid particles separately. The coefficients appearing in the equations have been represented by using real data for four nanofluid systems that use either water or ethylene glycol as the solvent, and CuO or Al_{2}O_{3} as the particles. The interaction coefficient *H* and the Green–Naghdi extra ‘viscosity’ coefficient have been estimated using various values, in the case of from 10^{−4} to 10^{−1}.

The findings are interesting. If we use the nanofluid Rayleigh numbers Ra, Ra or Ra_{1}, then we see that, for small values, we may certainly obtain a large reduction in the Rayleigh number when compared with that for a classical Newtonian fluid with no particles, i.e. *Ra*≈1707. This is completely in line with the findings of Tzou (2008), although we stress that we are incorporating non-Newtonian effects in the nanofluid owing to the presence of very small metallic oxide particles. We believe this is the first time this non-Newtonian effect has been modelled theoretically. The reduction of the Rayleigh number, very significantly in some cases (see when from 1707 to 277), means convective motion occurs much more easily in the nanofluid suspension. This, in turn, means heat transfer occurs more readily, and this agrees with the perceived use of a nanofluid in a heat transfer device.

We might note from table 1 that for *H* small, as decreases, Ra is approaching its value for a Newtonian fluid of 1707.

If we adopt the modified definition of *ρ* as given in equation (2.16), then the Rayleigh number values are slightly higher, but only by the factor (1−*ϕ*)^{−1}, although the numerical values obtained will actually be smaller since then *α* is effectively zero.

## Acknowledgements

I am indebted to two anonymous referees for their careful reading of the manuscript and for their helpful comments that have led to marked improvements in this article.

## Footnotes

- Received October 5, 2009.
- Accepted January 5, 2010.

- © 2010 The Royal Society