Thermal convection in a higher-gradient Navier–Stokes fluid

We discuss models for flow in a class of generalized Navier–Stokes equations. The work concentrates on producing models for thermal convection, analysing these in detail, and deriving critical Rayleigh and wave numbers for the onset of convective fluid motion. In addition to linear instability theory we present a careful analysis of fully nonlinear stability theory. The theories analysed all possess a bi-Laplacian term in addition to the normal spatial derivative term. The theories discussed are Stokes couple stress theory, dipolar fluid theory, Green–Naghdi theory, Fried–Gurtin–Musesti theory, and a second theory of Fried and Gurtin. We show that the Stokes couple stress theory and the Fried–Gurtin–Musesti theory involve the same partial differential equations while those of Green–Naghdi and dipolar theory are similar. However, we concentrate on boundary conditions which are crucial to understand all five theories and their differences.


Introduction
There has been much recent interest in models in fluid mechanics where the equations are generalizations of the Navier-Stokes equations. Of particular interest to the present article are a class where instead of just having the viscous term which involves the Laplacian of the velocity field, there is also present a higher spatial derivative term, namely one involving the bi-Laplacian operator of the velocity field.
In classical viscous fluid dynamics the Cauchy stress tensor, T i j , depends linearly on the symmetric part of the velocity gradient, where v(x, t) denotes the velocity field at a point x at time t, and v i, j ∂v i /∂ x j . For situations where the molecular structure of the fluid may involve long molecules or a suspension of such, or antisymmetric effects may be important, theories have been developed which generalize the classical approach and dictate that the Cauchy stress tensor depends not just upon the velocity gradient, but also upon the spatial derivative of the velocity gradient, i.e. upon v i, jk ≡ ∂ 2 v i /∂ x j ∂ x k . This article concentrates on analysing five separate theories where the total stress tensor depends also upon the second derivatives of velocity, but the final equations are either the same, or very similar. We generalize these models to include temperature effects, thereby encompassing the area of thermal convection. The models we investigate are the dipolar fluid theory of Bleustein and Green [1] with the dipolar inertia of Green and Naghdi [2], the couple stress theory of Stokes [3], the incompressible fluid model of Green and Naghdi [4], the viscous fluid theory of Fried and Gurtin [5] with the important higher-order stress contribution of Musesti [6], and the model of Fried and Gurtin [5] involving velocity gradients in the kinetic energy. To incorporate the temperature field we employ the energy balance law in each case together with a Boussinesq approximation, see, e.g. [7], to obtain a tractable set of equations in each case.
It is worth drawing attention to the fact that thermal convection studies involving a Navier-Stokes fluid are topics of much recent interest and importance in real life, see, e.g. [8][9][10][11][12][13][14]. It is interesting to wonder if these effects could be much richer with hyperstress effects introduced by the presence of higher spatial gradients in the momentum equation. In addition, thermal convection in a fluid is yielding very novel and important results in the field of renewable energy, see [15].
One of the key developments of the present work is to pay particular attention to the correct form of boundary conditions which must be used when higher-order derivatives are present. A set of partial differential equations is only useful if one can completely and correctly prescribe suitable boundary conditions which are applicable to real-life situations. Since for all five models we investigate the momentum equation involves the bi-Laplacian term 2 v i , being the three-dimensional Laplace operator, finding the correct boundary conditions to employ is far from a trivial matter. In this regard we rely on the fundamental work of Fried and Gurtin [5] regarding boundary conditions, although the extension of Musesti [6] to the higher-order stress tensor is crucial. a e-mail: brian.straughan@durham.ac.uk (corresponding author) We now present the five models for thermal convection in a generalized Navier-Stokes fluid. We then specialize these equations to analyse the problem of Bénard convection where a layer of fluid is heated from below and subsequently convective motion may ensue. The stability problem for this is analysed in detail from both a linear instability perspective, but also from a detailed nonlinear viewpoint using energy stability techniques. The boundary conditions are broken down into three main physically relevant cases, and detailed numerical results are presented and interpreted.

Generalized Navier-Stokes models for thermal convection
We now present five generalizations of the Navier-Stokes equations and then extend them in a form including the balance of energy equation for the temperature, where a Boussinesq approximation is utilized, cf. [7].

Green-Naghdi model
Green and Naghdi [4] derive a theory which involves two temperatures. This theory is adapted by Straughan [16] to be applicable to convection in a nanofluid. In this case the two temperatures are independent and are that associated with the fluid itself and to the nanoparticles, which are in suspension in the fluid. The basic momentum equation is Eq. (57) of Green and Naghdi [4], and this has form For an incompressible fluid the density ρ is constant and this equation is accompanied with the equation of continuity v i,i 0.
In (1) μ is the dynamic viscosity of the fluid, ν μ/ρ is the kinematic viscosity of the fluid, f i is an externally supplied body force, p(x, t) is the pressure, is the Laplacian. We employ standard indicial notation throughout together with the Einstein summation convention. For example, the divergence of the velocity field is for a function T depending upon x, t.
In this work we allow for temperature effects and suppose for each model that the density in the body force is linear in the temperature field, i.e.
where ρ 0 is a constant, T 0 is a reference temperature, and α is the coefficient of expansion of the fluid. The density is assumed constant everywhere else and then with a Boussinesq approximation, see [7], we may obtain the Green-Naghdi equations for thermal convection as In these equations κ is the thermal diffusivity, g is the gravity which is assumed acting in the negative z-direction, and k (0, 0, 1). An analogous model for thermal convection employing the Pavlovskii [17] equations is Hereμ is another constant.

The dipolar fluid
The dipolar fluid was introduced by Bleustein and Green [1] and most subsequent use involves the modified inertia coefficient of Green and Naghdi [2]. Thermal convection in a dipolar fluid is studied by Straughan [18] and we include it briefly for completeness, see also [19,20]. The equations for thermal convection in a dipolar fluid are given in [18] as (2.2), (2.12) and (2.14) and are In these equationsv i is the material derivative of v i ,d 2 ,γ are constants, and theγ term plays no role in a thermal convection analysis. The term (k j)i is a dipolar stress given by for constants h 1 , h 2 , h 3 with h 1 , h 3 > 0, where (kj) denotes the symmetric part in k and j.
To see the relevance of (6) to the work presented here we insert (7) into (6) 1 , and we expand thed 2 terms on the left and rearrange noting that a term of form −d 2 (v k, j v k, j ) , i arises which may be incorporated into the pressure gradient. In this way we may rewrite (6) as There is clearly a similarity between (8) and (4). The differences are that the coefficients of thed 2 and 2 terms are independent, the presence of theγ term, and the term on the left of form −d 2 v j,i v j . The latter term is very important in a nonlinear analysis.

The couple stress theory of V. K. Stokes
Condiff and Dahler [21] develop equations for polar fluids where the structural aspects of the molecules in the fluids are important. They relate this to an antisymmetric stress. Stokes [3] introduces the symmetric part of the velocity gradient, d i j into the Cauchy stress, but also introduces the skew-symmetric part of the velocity gradient, into the theory. He employs a curvature twist-rate tensor, K i j ω j,i , where ω i is the vorticity, and in addition to the Cauchy stress he requires a couple stress tensor. His momentum equation has a term of form 2 v i present. Devi and Mahajan [22] derive a Boussinesq approximation form of the Stokes couple stress theory to incorporate temperature effects. Their equations have form where ν > 0 is a constant. The form of this system of partial differential equations is exactly the same as what one finds from Navier-Stokes theory except that there is now present the higher derivative term −ν 2 v i . Devi and Mahajan [22] and Mahajan and Nandal [23] develop linear instability and nonlinear stability analyses for (9) and for an analogous system with a heat source, respectively. Both sets of writers prescribe boundary conditions on v i and T , but no discussion is included initially on further boundary conditions, a topic we return to in Sect. 4.

Fried-Gurtin-Musesti theory
A very clear development of a viscous flow theory with application to small length scales is due to Fried and Gurtin [5]. Given the current interest in microfluid dynamics we believe this model is highly relevant. Fried and Gurtin [5] employ a Cauchy stress but since they also include second gradients of the velocity field they introduce a third-order tensor, G i jk , they call a hyperstress.
The beauty of the work of Fried and Gurtin [5] is that they include a lucid and full development of boundary conditions. With higher-order derivatives appearing in the governing evolutionary equations the subject of boundary conditions is vital and the work of Fried and Gurtin [5] tackles this problem concisely. The constitutive theory for Fried and Gurtin's hyperstress is completed in an inspiring article by Musesti [6]. We here develop a theory for thermal convection in the Fried and Gurtin [5] and Musesti [6] model paying particular attention to the correct form of boundary conditions. The governing equations for incompressible, isothermal flow are given by Fried and Gurtin [5], and these are completed with the form for the hyperstress of Musesti [6], Eq. (17). We here apply a Boussinesq approximation to these equations and suppose the density ρ is a linear function of temperature when it appears in the body force term. In this manner we derive the following system of equations appropriate to thermal convection in an incompressible viscous fluid incorporating length scale effects, whereξ > 0 is a constant. The constantξ > 0 is that of Musesti [6], but we have divided by ρ 0 .
Remark We could add a linear term in T , i to the momentum Eq. (10) 1 , as appears in the dipolar fluid equations. However, this term does not affect the stability analysis and we omit it. If a nonlinear constitutive theory was adopted where products of terms involving T , i were present then the third-order tensor may lead to some interesting effects, cf. the explanation by Leslie [24] of the Lehmann effect. However, we are concentrating here on a constitutive theory linear in the dependent variables.
Existence results for a solution to the isothermal equations corresponding to (10) are presented in detail in [25] and in [26].
In this article we give a complete linear instability and nonlinear energy stability analysis for the Bénard problem according to (10). Importantly, we take account fully of the boundary conditions suggested by Fried and Gurtin [5], although we incorporate modifications necessitated by the article of Musesti [6]. We stress this issue of the boundary conditions, and we essentially analyse five classes. These are the strong adherence, the weak adherence, and the general adherence boundary conditions of Fried and Gurtin [5]. In addition we discuss the case of stress-free boundary conditions, and an illustrative case which allows an analytical solution.

Fried-Gurtin theory with generalized kinetic energy
Fried and Gurtin [5] present another theory in which the kinetic energy of the fluid contains in addition a term of form β|∇v| 2 /2, see equation (16) of their work. This, incorporating temperature via a Boussinesq approximation will result in a system of equations of form where we assume β > 0 is constant.

Basic state and perturbation equations for Bénard convection
We observe that (9) and (10) are exactly the same, apart from the notation of the coefficient of 2 v i . In what follows we replace the coefficient of 2 v i byξ in all five systems (4), (8), (9), (10) and (11). For consistency we replace μ 1 /ρ 0 μ byd 2 in (4), and we replaceμ byd 2 in (5). Suppose now each of the five systems holds in the horizontal layer where β is the temperature gradient, namely In each case the steady pressure is found as a function of z from the momentum equation.
wherep is the steady-state pressure. The equations for (u i , θ, φ, π) are derived and we non-dimensionalize with the scales The Rayleigh number Ra is The continuity and temperature perturbation equations are the same for all five systems and in non-dimensional form they are and The perturbation momentum equations are: and Pavlovskii model and Dipolar model and Fried-Gurtin-Musesti or Stokes couple stress model and Fried-Gurtin second model We obtain the appropriate set of non-dimensional perturbation equations for each of the five models by choosing the correct equation from (15)- (19) to accompany (13) and (14).

Boundary conditions
This is a key section. We concentrate on the model of Fried-Gurtin-Musesti in Sect. 2.4 and we follow the development of boundary conditions as in [5], although we incorporate the improvements to the constitutive theory given by Musesti [6]. The subject of boundary conditions in thermal convection is one of much recent interest. Great care has to be taken to ensure the correct form of boundary conditions is chosen to properly describe the physical problem in hand, cf. [27][28][29][30][31][32][33][34][35][36][37][38][39].
The governing equations for the perturbation variables are (18), (14) and (13). The solution is supposed periodic in x, y and satisfies a planform shape which tiles the plane, typically having a hexagonal shape, cf. the detailed discussion in [40, pp. 43-52]. The temperature boundary conditions are the usual ones and are θ 0 on z 0, 1.
Care has to be taken with the boundary conditions for u i . The discussion below applies to boundary conditions for the perturbation velocity. To handle the extra complexity of the second gradients of u i , Fried and Gurtin [5] introduce a symmetric Cauchy stress, T i j , which has the classical form whereπ is the fluid pressure. However, it is also necessary to introduce a third-order tensor which Fried and Gurtin [5] refer to as a hyperstress, G i jk , which has form whereπ k is a pressure vector, and G 0 i jk is the extra hyperstress involving only u i, jk . The total stress which arises in Eq. (18) is To describe boundary conditions on u i , Fried and Gurtin [5] introduce the total stress vector t i on the surface of a body. In the Bénard problem case the convection cell V is as shown in Fig. 1. The stress vector is We require the perturbation solution to satisfy a plane tiling pattern in (x, y), and so the boundary conditions on u i really only apply on the parts of ∂ V labelled as 1 and 2 in Fig. 1. However, we sometimes refer to the boundary conditions of Fried and Gurtin [5] on a boundary , of a general body ⊂ R 3 . Fried and Gurtin [5] introduce another vector on the surface of ∂ V as m i G i jk n j n k .
They show that the boundary conditions may, for a general domain with boundary , be divided into four classes. The first three of these are strong adherence, for which where ∂/∂n denotes the outward derivative in the direction of the unit outward normal, then weak adherence, for which and general adherence, for which We shall refer to the boundary conditions (25), (26) and (27) as cases I, II and III. In (27) we interpret as a non-dimensional version of the original parameters, μl, with now being at one's disposal. The fourth class of boundary conditions are those free of stress for which The boundary conditions of main interest for thermal convection are those of strong, weak, and general adherence. We briefly mention the stress-free case, but for a convection cell as in Fig. 1 having stress-free conditions on both 1 and 2 is purely an illustrative situation. We point out that the conditions of strong, weak and general adherence all apply to the situation of no-slip boundary conditions. The stress-free case will be important if we consider a layer fixed at the base but free at the upper surface. However, we do not consider this in detail here. We also do not consider the effects of surface tension, nor boundary slip, although those could be accounted for using the theory of Fried and Gurtin [5].
To see the relevance of the boundary conditions we multiply equation (23) by u i and integrate by parts for a general domain ⊂ R 3 with boundary , to find Under boundary conditions (25), (26) or (27) the boundary terms disappear and we then perform a further integration by parts on the G i jk,k term to obtain where ∂/∂n denotes the normal derivative on while s j ∇ s u i denotes the tangential surface derivative components. In fact, where u α are the surface coordinates, ; α denotes surface covariant differentiation, and a αβ is the fundamental form for the surface . Since u i 0 on , (29) reduces to whereas with general adherence conditions These results are useful in our nonlinear stability analysis. One might proceed directly from (23) by multiplying by u i and integrating over to find This is acceptable if we impose strong adherence boundary conditions where u i 0, ∂u i /∂n 0. It might be tempting to require as an alternative approach the boundary conditions It is convenient for the analysis which follows to employ the equivalent notation (u 1 , u 2 , u 3 ) ≡ (u, v, w). However, if we now specialize to the convection cell of Fig. 1, this means requiring w 0 and w 0,  In Sect. 6 we find the instability problem involves a sixth order differential equation for w and a second order one for θ and so we require eight boundary conditions. We have (36) and two boundary conditions on θ and so there are ten boundary conditions. Thus, the problem is overprescribed. Hence, (34) is not allowed and one must follow a procedure like that of Fried and Gurtin [5]. The tensor G i jk in (22) is given by Musesti [6] as where in this work η 1 , η 2 , η 3 are non-dimensional parameters which satisfy ξ η 1 − η 2 − 4η 3 > 0. This point is revisited in Sect. 5. The surface stress vector t i is given by Fried and Gurtin [5] for a general surface as t i T i j n j − 2G i jk,k n j + G i jk, n j n k n + G i jk K jk − 2K n j n k G i jk , where K jk is the curvature tensor, and K is the mean curvature.

Stability
We concentrate on the Fried-Gurtin-Musesti model and collect the relevant perturbation equations here. The other models will be discussed later in this section. The equations are where π contains both p andπ k contributions. Equations (39) are defined on R 2 × {z ∈ (0, 1)} × {t > 0}. We let V be a period cell for the solution, see figure 1, and we suppose the boundary conditions are of type I, II, III as defined in cases I, II, III in section 4. Let · and (·, ·) be the norm and inner product on L 2 (V ). System (39) together with any of the boundary conditions I, II or III is now written in the form where U (u, v, w, θ, π) T and where A is a linear symmetric operator, L is a linear operator, and N(U) represents the nonlinear terms, acting on a dense domain of the Hilbert space H [L 2 (V )] 4 , cf. [41, pp. 80-87]. Here A diag{1, 1, 1, Pr}, N(U) consists of u j u i, j and u i θ , i , while L is represented by the right hand side of (39).
We multiply (39) 1 by u i and integrate over V and then multiply (39) 3 by θ and integrate over V . After using the boundary conditions I or II we may obtain d dt where ϒ (G 0 i jk , u i, jk ). If boundary conditions III are employed then we must add to the right hand side of (41). The tensor G 0 i jk is given in (37) and we introduce the symmetric and skew-symmetric parts of u i, j , namely This allows one to rewrite ϒ as If we denote by ϒ 1 the pointwise form of (42), i.e. not integrated over V , then we may introduce the deviatoric and spherical parts of the tensors d i j,k and ω i j,k , cf. [5], in i and k as where δ ik d r j,r /3 is the spherical part while d i j,k d i j,k − δ ik d r j,r /3 is the deviatoric part, keeping j fixed, with a similar definition for ω i j,k . Then ϒ 1 may be written as Since the dissipation ϒ 1 is non-negative, Musesti [6, p. 85], one may deduce from (43) that From these we find Musesti [6] employed a clever choice of u to show η 1 − η 2 − 4η 3 ξ ≥ 0. The right hand side of (41) allows us to introduce a bilinear form on H. Let U (u 1 , v 1 , w 1 , θ, π 1 ) T and V (u 2 , v 2 , w 2 , φ, π 2 ) T be solutions to (39) subject to boundary conditions of case I, II or III. Let d i j , ω i j be the symmetric and skew-symmetric parts of u 1 i, j , and let e i j , ζ i j be the symmetric and skew-symmetric parts of u 2 i, j . Then define a bilinear form on H by where the boundary term is present in case III, but is not present for cases I and II. It may be shown that and so L is a symmetric operator on H, and thus the result of Galdi and Straughan [41] holds whereby the linear instability threshold is the same as the nonlinear stability one. We know in this case that (U, N (U )) 0 and so the nonlinear stability is global. Effectively, this may also be seen by examining the linear theory and showing that the growth rate is zero since leads to I mσ 0 where σ is the growth rate in a representation of time like e σ t , and U * is the complex conjugate of U. The resulting energy equation is then manipulated and the critical Rayleigh number for nonlinear energy stability is found using the maximum of the production and dissipation terms and the Euler-Lagrange equations which arise from this are found to be the same as the linear instability equations. However, the key result is that the physics of the nonlinear stability problem is completely captured by calculating the instability threshold from linear instability theory. While we have concentrated on the Fried-Gurtin-Musesti model the equivalence of the linear instability and nonlinear stability thresholds also applies to the Dipolar model, and to the Fried-Gurtin second model. To see this we need to show (U, N (U )) 0 for each of these. For the dipolar model this follows since the terms additional to the Fried-Gurtin-Musesti theory yield upon integration by parts, Note that we only use the fact that u i 0 on the 1 , 2 parts of ∂ V . For the Fried-Gurtin second model the extra terms yield For the Green-Naghdi model by multiplying (15)by u i and integrating over V , then multiplying (14) by θ and integrating over V , we find after integration by parts and use of the boundary conditions d dt The (U, N(U)) term is not zero, but we may derive a conditional result since where a value for the constant c is given in [42]. Provided ϒ ≥ k u 2 for some k > 0, then an energy inequality of form 1 then a continuity argument shows E → 0 exponentially and nonlinear stability follows, cf. [43, pp. 14-16]. It is obvious the restriction on ϒ holds in case I, and for cases II and III it follows provided η 1 > 0.
Remark The above proof does not work for the Pavlovskii model, Eq. (16).

Calculation of the critical Rayleigh number for instability
To find the critical Rayleigh numbers for instability the usual procedure is to linearize (39) and look for solutions like We do this here but since we know σ ∈ R then the linear stability threshold is when σ 0. Next, linearize, put σ 0 and then take curlcurl of (39) 1 to remove π. We retain the third component of the result and thus have to solve the system where * ∂ 2 /∂ x 2 + ∂/∂y 2 . We solve this system numerically using a Chebyshev tau method and to avoid problems with spurious eigenvalues it is better to solve instead the equivalent system where R 2 Ra. The system is to be solved with appropriate boundary conditions corresponding to cases I-III.

Case I boundary conditions
Here while on the lateral walls of the cell V the periodic boundary conditions apply. Thus, w 0, w z 0, u z 0, v z 0, on z 0, 1.
Since u x + v y + w z 0 in V , then on z 0, 1, and so w zz 0 on z 0, 1. Thus, the boundary conditions for case I are w w z w zz 0, z 0, 1.

Case II boundary conditions
In this case where i 1, 2, 3 and m i G i jk n j n k with n (0, 0, 1). For m α 0, α 1, 2, we find on z 0, 1, Differentiate the equation for α 1 by x, differentiate the equation for α 2 by y, and then use the continuity equation to see that on z 0, 1, This is rearranged as Since u i 0 on z 0, 1, for all i 1, 2, 3, the continuity equation shows that w z 0 there. Hence the boundary conditions for case II are These boundary conditions were also employed for convection in a dipolar fluid by Straughan [18], following the suggested boundary conditions by Bleustein and Green [1]. They were also employed by Straughan [16] who prescribed the couple at the wall, for convection in a nanofluid employing the Green and Naghdi [4] model.

Remark
The boundary condition m 3 0 is used to determine the pressure vector component π 3 .

Case III boundary conditions
In this case where i 1, 2, 3 and m i G i jk n j n k −π 3 δ i3 + η 1 u i, 33 By working with m 1 and m 2 and the continuity equation we may use arguments similar to those above to show the general adherence boundary conditions are w w z 0, and (3η 2 + 4η 3 ) * w z − ξw zzz w zz , on z 0, 1. Thus, since w z 0 on z 0, 1, the general adherence boundary conditions are where δ /ξ > 0.

Free surface boundary conditions
For two surfaces free of stress which do not deform the boundary conditions on z 0, 1 are w 0, t i 0, and m i 0.

Illustrative boundary conditions
There are other boundary conditions which have been employed in the literature, see, e.g. [22,23,44], where one poses w 0, w zz 0, w zzzz 0, θ 0, on z 0, 1. While I cannot place these in the Fried and Gurtin [5] context they may prove useful for some purposes in that they lead to an analytical solution.
To see this we write w W (z)h(x, y), θ (z)h(x, y), where h is a function compatible with tiling the plane, and h satisfies * h −a 2 h, a being a wavenumber, cf. [40, pp. 43-52]. Then, equations (45) become where D d/dz. Equations (50) applied on the boundary show D 6 W 0, D 2 0, on z 0, 1, and then this process may be repeated to deduce D 2n W 0, D 2n 0, on z 0, 1, where n 0, 1, 2, . . . Hence W, may be represented as sin series with typical element sin(nπ z). This then reduces (50) to an equation for a determinant which yields where n n 2 π 2 + a 2 . We require the smallest value of Ra and this is seen to be when n 1, so π 2 + a 2 . Note the similarity with the thermal convection problem for a Brinkman porous material, where 2 and 3 are present rather than 3 and 4 as here, cf. [45]. The minimum value of Ra in a 2 is found when The critical value of Ra(ξ ) is then given by employing this value of a 2 crit in (51).

Numerical solution
To solve equations (50) in the case of boundary conditions I-III we employ a Chebyshev tau method, cf. [46]. We solve (50) for the eigenvalues Ra by writing these equations as a system (D 2 − a 2 )W χ, The boundary conditions are written as,
The solution to (52) is written as a sum of Chebyshev polynomials of form This gives rise to a generalized matrix eigenvalue problem of form and the boundary conditions are incorporated into the matrix A by writing them into the appropriate rows of A. The matrix eigenvalue problem is solved by the QZ algorithm of Moler and Stewart [47].

Numerical results
We now report on results for Eq. (18) of the Fried-Gurtin-Musesti model. However, due to the fact that we have shown that with boundary conditions of type I-III the growth rate is real, the numerical values apply to all four models described in Sect. 3; by this we mean the dipolar fluid, Green and Naghdi [4] fluid, Fried-Gurtin-Musesti fluid, and the second fluid of Fried and Gurtin [5].
Numerical results for boundary conditions of type I and II are given in Tables 1 and 2, and in Figs. 2 and 3. It is seen that for the ξ strong adherence boundary conditions, case I, the critical Rayleigh number increases strongly for ξ small and then approaches a growth (approximately) linear in ξ . The same behaviour is true for weak adherence boundary conditions, case II, although the effect for small ξ is not so pronounced. The behaviour of the critical wavenumber for cases I and II is interesting and very different. For case I the critical value of a 2 increases rapidly for ξ small and then approaches a critical value, whereas for case II the opposite effect is seen. This means that while the ξ term, the hyperviscosity in the terminology of Fried-Gurtin-Musesti, is stabilizing for both cases the convection cell behaviour at the onset of convection is different. Since the wavenumber is inversely proportional to the aspect ratio of the convection cell (width to depth), increasing ξ for case I means the cells become more narrow for fixed depth. For case II boundary conditions the convection cells increase in width as ξ increases. The strong adherence boundary conditions are more stabilizing then the weak ones and appear to concentrate the convection cells more together to produce a smaller pattern when viewed from above. For the general adherence boundary conditions with δ /ξ, Figs. 4 and 5 show the critical Rayleigh number varies from the weak adherence value at δ 0, ( 0), to the strong adherence value as δ increases, ( → ∞). Note the difference in scale of Ra where ξ 0.01 in Fig. 4 whereas ξ 0.1 in Fig. 5. Figure 6 shows the variation in a 2 crit as δ varies, for ξ 0.01 and 0.1. We present one value of Ra crit and a 2 crit for two stress-free surfaces. This is when η 1 5.1, η 2 1, η 3 1, ξ 0.1, and then Ra crit 6889.914, a 2 crit 5.5. We only give one set of values, partly because two stress-free surfaces is not of great interest physically, but also we need values for the Musesti parameters η 1 , η 2 and η 3 (Table 3).

Conclusions
We have reviewed four models for fluid flow in what might be termed generalized Navier-Stokes theory. Generalized in the sense that the stress depends not only on the velocity gradient but also on the second spatial derivatives of velocity. These hydrodynamic models are physically very relevant since, for example, Green and Rivlin [48,49] and Green et al. [50] show that they are important when the molecular structure of the fluid involves long molecules, and they are increasingly important in the expanding microfluidic industry where length scales are very small, see [5].
We have incorporated temperature into all four models and we have investigated in detail the Bénard thermal convection problem for each. The results of linear instability are the same for each model provided we employ the physically realistic boundary conditions suggested by Fried and Gurtin [5]. The nonlinear stability boundaries are the same for three of the models although we cannot rule  The first two columns with Ra, a 2 are for ξ 0.01, whereas the second two columns with Ra, a 2 are for ξ 0.1 out sub-critical instabilities for the model of Green and Naghdi [4]. When higher spatial gradients appear in the time derivative terms then the measure of decay in the stable regime is stronger, being in H 1 as opposed to L 2 in the case where only the material time derivative of the velocity appears. While the instability thresholds are the same for thermal convection in all four classes of generalized Navier-Stokes equations considered here we point out that this will not be the case for more exotic convection problems. For example, if rotation is present, magnetic fields are present, or thermosolutal convection is considered then oscillatory convection will be encountered and those theories which involve terms like − u i,t on the left hand side will witness different critical Rayleigh number and wavenumber behaviours.
We observe that Moon et al. [51,52] have identified interesting attractors and behaviour for ordinary differential equation systems derived from double-diffusive convection using Navier-Stokes theory. It would be interesting to know how the inclusion of higher spatial gradients would influence such behaviour.
Finally, in connection with the currently vital topic of energy creation, ElFatnani et al. [15] describe a novel technique where a ceramic plate is heated and cooled by being placed above a bath of oil which is undergoing thermal convection. The variation of temperature in the ceramic plate creates electricity via the pyroelectric effect. It is interesting to wonder whether a fluid with long molecules, or in a situation where micro length scales dominate, could be involved to improve this technique of generating electricity. This could be a practical use for the Fried-Gurtin-Musesti theory of generalized Navier-Stokes fluids.