Viscoelastic Natural Convection

Heat convection occurs in natural and industrial processes due to the presence of temperature gradients which may appear in any direction with respect to the vertical, which is determined by the direction of gravity. In this case, natural convection is the fluid motion that occurs due to the buoyancy of liquid particles when they have a density difference with respect the surrounding fluid. Here, it is of interest the particular problem of natural convection between two horizontal parallel flat walls. This simple geometry brings about the possibility to understand the fundamental physics of convection. The results obtained from the research of this system may be used as basis to understand others which include, for example, a more complex geometry and a more complex fluid internal structure. Even though it is part of our every day life (it is observed in the atmosphere, in the kitchen, etc.), the theoretical description of natural convection was not done before 1916 when Rayleigh [53] made calculations under the approximation of frictionlesswalls. Jeffreys [27] was the first to calculate the case including friction in the walls. The linear theory can be found in the monograph by Chandrasekhar [7]. It was believed that the patterns (hexagons) observed in the Benard convection (see Fig. 1, in Chapter 2 of [7] and the references at the end of the chapter) were the same as those of natural convection between two horizontal walls. However, it has been shown theoretically and experimentally that the preferred patterns are different. It was shown for the first time theoretically by Pearson [45] that convection may occur in the absence of gravity assuming thermocapillary effects at the free surface of a liquid layer subjected to a perpendicular temperature gradient. The patterns seen in the experiments done by Benard in the year 1900, are in fact only the result of thermocapillarity. The reason why gravity effects were not important is that the thickness of the liquid layer was so small in those experiments that the buoyancy effects can be neglected. As will be shown presently, the Rayleigh number, representative of the buoyancy force in natural convection, depends on the forth power of the thickness of the liquid layer and the Marangoni number, representing thermocapillary effects, depends on the second power of the thickness. This was not realized for more than fifty years, even after the publication of the paper by Pearson (as seen in the monograph by Chandrasekhar). Natural convection may present hexagonal patterns only when non


Introduction
Heat convection occurs in natural and industrial processes due to the presence of temperature gradients which may appear in any direction with respect to the vertical, which is determined by the direction of gravity. In this case, natural convection is the fluid motion that occurs due to the buoyancy of liquid particles when they have a density difference with respect the surrounding fluid. Here, it is of interest the particular problem of natural convection between two horizontal parallel flat walls. This simple geometry brings about the possibility to understand the fundamental physics of convection. The results obtained from the research of this system may be used as basis to understand others which include, for example, a more complex geometry and a more complex fluid internal structure. Even though it is part of our every day life (it is observed in the atmosphere, in the kitchen, etc.), the theoretical description of natural convection was not done before 1916 when Rayleigh [53] made calculations under the approximation of frictionless walls. Jeffreys [27] was the first to calculate the case including friction in the walls. The linear theory can be found in the monograph by Chandrasekhar [7]. It was believed that the patterns (hexagons) observed in the Bénard convection (see Fig. 1, in Chapter 2 of [7] and the references at the end of the chapter) were the same as those of natural convection between two horizontal walls. However, it has been shown theoretically and experimentally that the preferred patterns are different. It was shown for the first time theoretically by Pearson [45] that convection may occur in the absence of gravity assuming thermocapillary effects at the free surface of a liquid layer subjected to a perpendicular temperature gradient. The patterns seen in the experiments done by Bénard in the year 1900, are in fact only the result of thermocapillarity. The reason why gravity effects were not important is that the thickness of the liquid layer was so small in those experiments that the buoyancy effects can be neglected. As will be shown presently, the Rayleigh number, representative of the buoyancy force in natural convection, depends on the forth power of the thickness of the liquid layer and the Marangoni number, representing thermocapillary effects, depends on the second power of the thickness. This was not realized for more than fifty years, even after the publication of the paper by Pearson (as seen in the monograph by Chandrasekhar). Natural convection may present hexagonal patterns only when non Boussinesq effects [52] occur, like temperature dependent viscosity [57] which is important when temperature gradients are very large. The Boussinesq approximation strictly assumes that all the physical parameters are constant in the balance of mass, momentum and energy equations, except in the buoyancy term in which the density may change with respect to the temperature. Any change from this assumption is called non Boussinesq approximation. When the thickness of the layer increases, gravity and thermocapillary effects can be included at the same time [40]. This will not be the subject of the present review. Here, the thickness of the layer is assumed large enough so that thermocapillary effects can be neglected.
The effects of non linearity in Newtonian fluids convection were taken into account by Malkus and Veronis [33] and Veronis [65] using the so called weakly non linear approximation, that is, the Rayleigh number is above but near to the critical Rayleigh number. The small difference between them, divided by the critical one, is used as an expansion parameter of the variables. The patterns which may appear in non linear convection were investigated by Segel and Stuart [57] and Stuart [61]. The method presented in these papers is still used in the literature. That is, to make an expansion of the variables in powers of the small parameter, including normal modes (separation of horizontal space variables in complex exponential form) of the solutions of the non linear equations. With this method, an ordinary non linear differential equation (or set of equations), the Landau equation, is obtained for the time dependent evolution of the amplitude of the convection cells. Landau used this equation to explain the transition to turbulent flow [31], but never explained how to calculate it. For a scaled complex A(t), the equation is: In some cases, the walls are considered friction free (free-free case, if both walls have no friction). One reason to make this assumption is that the nonlinear problem simplifies considerably. Another one is that the results may have relevance in convection phenomena in planetary and stellar atmospheres. In any way, it is possible that the qualitative results are similar to those of convection between walls with friction, mainly when the interest is on pattern formation. This simplification has also been used in convection of viscoelastic fluids. To describe the nonlinear envelope of the convection cells spatial modulation, it is possible to obtain a non linear partial differential equation by means of the multiple scales approximation [3], as done by Newell and Whitehead [39] and Segel [56]. This equation is called the Newell-Whitehead-Segel (NWS) equation. For a scaled A(X,Y,T), it is: Here, X, Y and T are the scaled horizontal coordinates and time, respectively. In the absence of space modulation it reduces to the Landau Equation 1. It is used to understand the non linear instability of convection flow. However, it has been found that this equation also appears in the description of many different physical phenomena. The non linear stability of convection rolls depends on the magnitude of the coefficients of the equation. If the possibility of the appearance of square or hexagonal patterns is of interest, then the stability of two coupled or three coupled NWS equations have to be investigated. They are obtained from the coupling of modes having different directions (see [22] and [23]).
The shear stress tensor of Newtonian fluids have a linear constitutive relation with respect to the shear rate tensor. The constitutive equation of that relation has as constant of proportionality the dynamic viscosity of the fluid, that is Here, the shear rate tensor is Any fluid whose stress tensor has a different constitutive relation, or equation, with respect to the shear rate tensor is called non Newtonian. That relation might have an algebraic or differential form. Here, only natural convection of viscoelastic fluids will be discussed [4,9] as non Newtonian flows. These fluids are defined by constitutive equations which include complex differential operators. They also include relaxation and retardation times. The physical reason can be explained by the internal structure of the fluids. They can be made of polymer melts or polymeric solutions in some liquids. In a hydrostatic state, the large polymeric chains take the shape of minimum energy. When shear is applied to the melt or solution, the polymeric chains deform with the flow and then they are extended or deformed according to the energy transferred by the shear stress. This also has influence on the applied shear itself and on the shear stress. When the shear stress disappears, the deformed polymeric chains return to take the form of minimum energy, carrying liquid with them. This will take a time to come to an end, which is represented by the so called retardation time. On the other hand, there are cases when shear stresses also take some time to vanish, which is represented by the so called relaxation time. It is possible to find fluids described by constitutive equations with both relaxation and retardation times. The observation of these viscoelastic effects depend on different factors like the percentage of the polymeric solution and the rigidity of the macromolecules.
A simple viscoelastic model is the incompressible second order fluid [10,16,34]. Assuming τ ij as the shear-stress tensor, the constitutive equation is: and for a tensor P ij and where D Dt is the Lagrange or material time derivative. The time derivative in Equation 6 is called the lower-convected time derivative, in contrast to the following upper convected time derivative and to the corrotational time derivative where the rotation rate tensor is These time derivatives can be written in one formula as where the time derivatives correspond to the upper convected for a = 1, the corrotational for a = 0 and the lower convected for a = −1, respectively [47]. These time derivatives are invariant under a change of reference frame. In Equations 3 and 5 η 0 is the viscosity and in Equation 5 β and γ are material constants. The second order model Equation 5 has limitations in representing fluid motion. It is an approximation for slow motion with small shear rate [4]. Linear and nonlinear convection of second order fluids has been investigated by Dávalos and Manero [12] for solid walls under the fixed heat flux boundary condition. The same fluid has been investigated looking for the possibility of chaotic motion (aperiodic and sensitive to initial conditions [28]) by [58] for the case of free boundaries and fixed temperature boundary condition.
The Maxwell model [4] is used to describe motion where it is possible to have shear stress relaxation. The constitutive equation of this model is: where λ is the relaxation time. A characteristic of this equation is that for λ small the fluid nearly behaves as Newtonian. For large λ it tends to behave as an elastic solid as can be seen if e ij is considered as the time derivative of the strain. In the limit of very large λ, the approximate equation is integrated in time to get Hook's law, that is, the stress is proportional to the strain. This constitutive equation has three versions, the upper convected, the lower convected and the corrotational Maxwell models, depending on the time derivative selected to describe the fluid behavior. The natural convection of the Maxwell fluid has been investigated by Vest and Arpaci [66] for free-free and solid-solid walls with fixed temperature. Sokolov and Tanner [59] investigated the linear problem of the Maxwell fluid, among other viscoelastic fluids, using an integral form of the stress tensor. The non linear problem has been investigated for free-free boundaries by Van Der Borght et al. [64], using the upper convected time derivative. Brand and Zielinska [5] show that nonlinear traveling waves appear for different Prandtl numbers in a convecting Maxwell fluid with free-free walls. The Prandtl number Pr is the ratio of the kinematic viscosity over the thermal diffusivity. The chaotic behavior of convection of a Maxwell fluid has been investigated by Khayat [29]. The effect of the thickness and thermal conductivity of the walls has been taken into account in the linear convection of a Maxwell fluid by Pérez-Reyes and Dávalos-Orozco [46].
The Oldroyd's fluid model [4,41] includes, apart from a relaxation time, a retardation time.
The linear version of this model is called the Jeffreys model (but the non linear model is sometimes called by this name). The constitutive equation is where λ 1 is the retardation time. Notice that when λ 1 = 0 this contitutive equation reduces to that of a Maxwell fluid. Therefore, a number of papers which investigate the convection in Oldroyd fluids also include results of the Maxwell fluid. When λ = 0, the equation reduces to that of the second-order fluid with a zero coefficient γ. Linear convection of Oldroyd fluids has been investigated by Green [21], Takashima [62], Kolkka and Ierley [30], Martínez-Mardones and Pérez-García [35] and Dávalos-Orozco and Vázquez-Luis [14] for free upper surface deformation. Nonlinear calculations of the Oldroyd fluid where done first by Rosenblat [55] for free-free boundaries. The non linear problem of solid-solid and solid-free boundaries was investigated by Park and Lee [43,44]. Nonlinear problems were investigated by Martínez-Mardones et al. for oscillatory and stationary convection [36], to study the stability of patterns in convection [37] and to investigate the convective and absolute instabilities by means of amplitude equations [38].
The following section presents the balance equations suitable for natural convection. Section 3 is an introduction to Newtonian fluids convection. The Sections 4, 5, and 6 correspond to reviews of convection of second-order, Maxwell and Oldroyd fluids, respectively. Finally, some conclusions are given in the last Section 7.

Equations of balance of momentum, mass and energy
Here, the basic equations of balance of momentum, mass and energy for an incompressible fluid are presented. In vector form, they are The dimensional variables are defined as follows. ρ is the density, u * = (u * , v * , w * ) is the velocity vector, p * is the pressure, τ τ τ * is the stress tensor which satisfies one of the constitutive equations presented above. T * is the temperature, C V is the specific heat at constant volume and X F is the heat conductivity of the fluid. Use is made of ∇ * = (∂/∂x * , ∂/∂y * , ∂/∂z * ). g = (0, 0, −g) = −gk is the vector of the acceleration of gravity with g its magnitude andk a unit vector in the direction opposite to gravity. Equation 15 means that the fluid is incompressible and that any geometric change of a fluid element volume in one direction is reflected in the other the directions in such a way that the volume is preserved according to this equation.
If the thickness and conductivity of the walls are taken into account, the temperature in each wall satisfies the equation where T W is the temperature of one of the walls (T L for the lower wall and T U for the upper wall). ρ W , C VW and X W are the density, specific heat at constant volume and heat conductivity of one of the walls (ρ L , C VL , X LW for the lower wall and ρ U , C VU , X UW for the upper wall). The variables are subjected to boundary conditions. The velocity has two types of conditions: for friction free walls and for solid walls with friction. They are, respectively: n · u * = 0 and n · ∇ * u * = n · ∇ * v * = 0 at z * = z 1 and z * = z 1 + d f ree boundary (18) u * = 0 at z * = z * 1 and z * = z * 1 + d solid boundary where n is the unit normal vector to one of the walls, z 1 is a particular position of the lower wall in the z-axis and d is the thickness of the fluid layer. The conditions in the first line of Equation 18 mean that the fluid can not penetrate the wall and that the wall does not present any shear due to the absence of friction. The condition of the second line means that the fluid sticks to the wall due to friction.
The temperature satisfies the boundary conditions of fixed temperature and fixed heat at the walls, respectively, where q 0 is a constant heat flux normal and through one of the walls.
If the thickness and heat conductivity of the walls are taken into account, the temperature has to satisfy the conditions The equations and boundary conditions can be made non dimensional by means of representative magnitudes for each of the dependent and independent variables. For example, the distance is scaled by the thickness of the fluid layer d or a fraction of it, the time is scaled with d 2 /κ, where the thermal diffusivity is κ = X F /ρ 0 C V , the velocity with κ/d, the pressure and the stress tensor with ρ 0 κ 2 /d 2 . ρ 0 is a representative density of the fluid. The temperature is made non dimensional with a characteristic temperature difference or with a quantity proportional to a temperature difference. The time can also be scaled with d 2 /ν, where the kinematic viscosity is ν = η 0 /ρ 0 and the velocity with ν/d. Then, the pressure and the stress tensor can be scaled in two ways, by means of ρ 0 κν/d 2 or by ρ 0 ν 2 /d 2 . The difference stems on the importance given to the mass diffusion time d 2 /ν or to the heat diffusion time d 2 /κ. It is assumed that before a perturbation is applied to the Equations 14 to 16 the system is in a hydrostatic state and that the variables satisfy The solution of these two equations will be the main pressure P 0 and the main temperature profile T * 0 of the system before perturbation. Here, ρ 0 is a reference density at the reference temperature T * R which depends on the boundary conditions. α T is the coefficient of volumetric thermal expansion of the fluid. These two solutions of Equations 21 and 22 are subtracted from Equations 14 to 16 after introducing a perturbation on the system. In non dimensional form, the equations of the perturbation are 1 Pr The non dimensionalization was based on the heat diffusion time and the scaling of the pressure and shear stress with ρ 0 κν/d 2 . u, p, τ, θ and θ W are the perturbations of velocity, pressure, shear stress, fluid temperature and walls temperature (θ L and θ U for the lower and upper walls), respectively. R = gα T ΔTd 3 /νκ is the Rayleigh number and Pr = ν/κ is the Prandtl number. ΔT is a representative temperature difference. κ W is the thermal diffusivity of one of the walls (κ L for the lower wall and κ U for the upper wall).
The last term in the left hand side of Equation 25 appears due to the use of the linear temperature solution of Equation 22. If the temperature only depends on z * in the form Here, a 1 is a constant which is proportional to a temperature difference or an equivalent if the heat flux is used. In these equations, the Boussinesq approximation has been taken into account, that is, in Equations 14 to 16 the density was assumed constant and equal to ρ 0 everywhere except in the term ρg where it changes with temperature. The other parameters of the fluid and wall are also assumed as constant. These conditions are satisfied when the temperature gradients are small enough.
The boundary conditions of the perturbations in non dimensional form are n · u = 0 and n · ∇u = n · ∇v = 0 at z = z 1 and z = z 1 + 1 f ree boundary (31) u = 0 at z = z 1 and z = z 1 + 1 solid boundary Here, D L and D U are the ratios of the thickness of the lower and upper walls over the fluid layer thickness, respectively. The meaning of the conditions Equation 32 is that the original temperature and heat flux at the boundary remain the same when θ = 0 and n · ∇θ = 0. The same can be said from the first two Equations 33, that is, the temperature below the lower wall and the temperature above the upper wall stay the same after applying the perturbation.

Natural convection in newtonian fluids
The basics of natural convection of a Newtonian fluid are presented in this section in order to understand how other problems can be solved when including oscillatory and non linear flow. The section starts with the linear problem and later discuss results related with the non linear equations. The system is a fluid layer located between two horizontal and parallel plane walls heated from below or cooled from above. Gravity is in the z-direction. As seen from Equations 23 to 25, the linear equations are 1 Pr In Equation 34 use has been made of Equations 3, 4 and 35. The first boundary conditions used will be those of free-free and fixed temperature at the wall [7]. These are the simplest conditions which show the qualitative behavior of convection in more complex situations.
To eliminate the pressure from the equation it is necessary to apply once the curl operator to Equation 34. This is the equation of the vorticity and its vertical component is independent from the other components of the vorticity vector. Applying the curl one more time, it is possible to obtain an equation for the vertical component of the velocity independent of the other components. The last and the first equations are 1 Pr Here, ∇ 2 ⊥ = (∂ 2 /∂ 2 x, ∂ 2 /∂ 2 y) is the horizontal Laplacian which appears due to the unit vector k in Equation 34. The third component of vorticity is defined by ζ Z = [∇ × u] Z . The three components of ζ are related with the three different elements of the rotation tensor Equation 10 multiplied by 2.
These Equations 37, 38 and 36 are partial differential and their variables can be separated in the form of the so called normal modes is the amplitude of the dependent variable. The wavevector is defined by k = (k x , k y ), k x and k y are its x and y-components and its magnitude is k = k. When the flow is time dependent, σ = s + iω where s is the growth rate and ω is the frequency of oscillation. Then, using normal modes and assuming that the system is in a neutral state where the growth rate is zero (s = 0), the equations are iω Pr The last equation is the equation of continuity and U and V are the z-dependent amplitudes of the x and y-components of the velocity. Θ, W and Z are the z-dependent amplitudes of the temperature and the vertical components of velocity and vorticity, respectively. If the heat diffusion in the wall is taken into account with a temperature amplitude Θ W , Equation 26 becomes In normal modes the boundary conditions Equations 31 to 33 change into those for the amplitude of the variables. They are where the two free boundary conditions n · ∇u = n · ∇v = 0 where replaced by d 2 W/dz 2 = 0 using the z-derivative of the continuity Equation 43 and the x and y-components of the solid boundary condition u = 0 are replaced by dw/dz = 0 using Equation 43. The conditions of the z-component of vorticity Z are obtained from its definition using the x and y-components of the boundary condition u = 0 for a solid wall and the derivative for the free wall. Notice that in the linear problem, in the absence of a source of vorticity (like rotation, for example) for all the conditions investigated here, the vorticity Z = 0. Vorticity can be taken into account in the non linear problem (see for example Pismen [49] and Pérez-Reyes and Dávalos-Orozco [48]). Equations 40 and 42 are independent of the vorticity Z. They can be combined to give The first boundary conditions used will be those of Equations 45 and 46 of free-free and fixed temperature at the wall [7]. These are the simplest conditions, important because they allow to obtain an exact solution of the problem and may help to understand the qualitative behaviour of convection in more complex situations. Using Equations 40 and 42 evaluated at the boundaries, these conditions can be translated into: where D = d/dz. These are satisfied by a solution Here, n is an integer number and A is the amplitude which in the linear problem can not be determined. In this way, substitution in Equation 48 leads to an equation which can be written as The Rayleigh number must be real and therefore the imaginary part should be zero. This condition leads to the solution ω = 0. That is, under the present boundary conditions, the flow can not be oscillatory, it can only be stationary. Thus, the marginal Rayleigh number for stationary convection is From this equation it is clear that n represents the modes of instability that the system may show. If the temperature gradient is increased, the first mode to occur will be that with n = 1.
Here, the interest is to calculate the lowest magnitude of R with respect to k because it is expected that the mode n = 1 will appear first for the wavenumber that minimizes R. Thus, taking the derivative of R with respect to k and solving for k, the minimum is obtained for This critical wavenumber is substituted in Equation 52 to obtain the critical Rayleigh number for free-free walls and fixed temperature This may be interpreted as the minimum non dimensional temperature gradient needed for the beginning of linear convection. The space variables were made non dimensional using the thickness of the layer. Therefore, the result of Equation 53 shows that, under the present conditions, the size (wavelength) of the critical convection cells will be Λ = 2 √ 2d.
Now, the solid-solid conditions and the fixed temperature conditions are used. In this case too, it has been shown that convection should be stationary [7]. The problem is more complicated and numerical methods are needed to calculate approximately the proper value problem for R. From Equation 48 and ω = 0 the equation for W is : Due to the symmetry of the boundary conditions it is possible to have two different solutions.
One is even and the other is odd with respect to boundary conditions located symmetrically with respect to the origin of the z-coordinate. That is, when in Equation 55 to obtain This is a sixth degree equation for m (or a third degree equation for m 2 ) which has to be solved numerically. The solutions for Equation 55 can be written as one of this coeficients A i has to be normalized to one. The m i 's contain the proper values R and k. Three conditions for W are needed at each wall. They are W = DW = D 4 W − 2k 2 D 2 W at z = ±1/2. The last one is a result of the use of the condition for Θ in Equation 40 with ω = 0. As an example of the even and odd proper value problem, the evaluation of the conditions of W at z = ±1/2 will be presented. They are 3 2 Addition and subtraction of both conditions give, respectively The same can be done with the other boundary conditions. The important point is that two sets can be separated, each one made of three conditions: one formed by the addition of the conditions and another one made of the subtraction of the conditions, that is, the even and the odd modes of the proper value problem. It has been shown numerically [7] that the even mode gives the smaller magnitude of the marginal proper value of the Rayleigh number. The odd mode gives a far more larger value and therefore it is very stable in the present conditions of the problem. However, there are situations where the odd mode can be the first unstable one (see Ortiz-Pérez and Dávalos-Orozco [42] and references therein). Recently, Prosperetti [51] has given a very accurate and simple formula for the marginal Rayleigh number by means of an improved numerical Galerkin method. That is The marginal curve plotted from this equation gives a minimum R = 1715.08 at k = 3.114.
These critical values are very near to those calculated by means of very accurate but complex numerical methods. The accepted values are R = 1707.76 at k = 3.117 [7]. From the critical wavenumber it is possible to calculate the size (wavelength) of the cell at onset of convection. That is, Λ = 2πd/3.117, which is smaller than that of the free-free case. This is due to the friction at the walls. Walls friction also stabilizes the system increasing the critical Rayleigh number over two and a half times the value of the free-free case.
Linear convection inside walls with fixed heat flux has been investigated by Jakeman [26]. Hurle et al. [25] have shown that the principle of exchange of instabilities is valid for a number of thermal boundary conditions, that is, oscillatory convection is not possible and ω = 0, including the case of fixed heat flux. Jakeman [26] used a method proposed by Reid and Harris [7,54] to obtain an approximate solution of the proper value of R. This is a kind of Fourier-Galerkin method [17,18]. From the expresion obtained, the critical Rayleigh and wavenumber were calculate analytically by means of a small wavenumber approximation. The reason is that it has been shown numerically in the marginal curves, that the wavenumber of the smallest Rayleigh number tends to zero. The critical Rayleigh number for the free-free case is R = 120 = 5! (k = 0) and for the solid-solid case R = 720 = 6! (k = 0). It is surprising that the critical Rayleigh numbers are less than half the magnitude of those of the fixed temperature case. This can be explained by the form of the temperature boundary condition, that is Dθ = 0. From the condition it is clear that the perturbation heat flux can not be dissipated at the wall and therefore the perturbation remains inside the fluid layer.
This makes the flow more unstable and consequently the critical Rayleigh number is smaller. The linear problem was generalized by Dávalos [11] to include rotation and magnetic field, and obtained explicit formulas for the critical Rayleigh number depending on rotation and magnetic field and a combination of both. Notice that the method by Reid and Harris [54] is very effective and it is still in use in different problems of convection. However, as explained above, it has been improved by Prosperetti [51].
The nonlinear problem for the fixed heat flux approximation has been done by Chapman and Proctor [8]. They improved the methods of calculation with respect to previous papers. The method to obtain a nonlinear evolution equation is to make an scaling of the independent and dependent variables taking into account that, if the Rayleigh number is a little far above criticality (weak nonlinearity), the wavenumber of the convection cell is still small. Consequently, the motion will be very slow because the cell is very large. Then, the scaling used in the nonlinear Equations 23 to 25 is as follows They solve a two-dimensional problem using the stream function by means of which the velocity vector field satisfies automatically Equation 24. The velocity with the scaling is defined as u = (∂ψ/∂z, 0, − ∂ψ/∂X). The stream function is also scaled as ψ = φ. The expansion of the functions in terms of is The reason for this expansion is that the substitution of the scaling in Equations 23 to 25 only shows even powers of . The problem is solved in different stages according to the orders of subjected to the corresponding scaled boundary conditions. The critical value of R is obtained from a solvability condition at O( 2 ). Notice that they locate the walls at z = ±1 and obtain R c = 15/2 and R c = 45 for the free-free and the solid-solid cases, respectively. If the definition of the Rayleigh number includes the temperature gradient it depends on the forth power of the thickness of the layer. Therefore, the Rayleigh number defined here is sixteen times that defined by Chapman and Proctor [8]. The evolution equation is obtained as a solvability condition at O( 4 ). That is, with θ 0 = f (X, τ) This equation is valid for free-free and solid-solid boundary conditions and the constants B and C have to be calculated according to them. Chapman and Proctor found that the patterns of convection cells are rolls but that they are unstable to larger rolls. Therefore, the convection will be made of only one convection roll. An extension of this problem was done by Proctor [50] including the Biot number Bi in the thermal boundary conditions. The Biot number is a non dimensional quantity that represents the heat flux across the interface between the fluid and the wall. The fixed heat flux boundary condition is obtained when the Bi is zero. When Bi is small but finite, the critical convection cells are finite and therefore more realistic.
The effect of the thickness and heat conductivity of the wall on natural convection has been investigated by Cericier et al. [6]. The goal is to be able to obtain more realistic critical values of the Rayleigh number and wavenumber. Here it is necessary to use the thermal conditions Equations 22 and 33 for the main temperature profile and the perturbation of temperature, respectively. The main temperature profile is linear with respect to z but is more complex due to the presence of terms which depend on all the new extra parameters coming from the geometry and properties of the walls. From Equation 33 it is possible to calculate a new condition which only has the amplitude of the fluid temperature perturbation. In the present notation it has the form where the coefficients of Θ in both conditions might be considered the Biot numbers of each wall. Here, q = √ k 2 + iω. Now there are four new parameters which influence the convective instability, the heat conductivities ratio X L and X U and the thicknesses ratio D L and D U . Assuming that X = X L = X U and d = D L = D U the problem has some simplification. Figure 1 shows results for the case when the properties of both walls are the same. Notice that when X increases the critical values in both figures change from the fixed temperature case to the fixed heat flux case. The results are similar to those of Cericier et al. [6] and were plotted using a formula calculated from a low order Galerkin approximation. It is important to observe that in the middle range of X the thickness of the walls play a relevant roll producing large differences between the critical values, for fixed X.
The problem of surface deformation in convection requires lower conditions for free or solid walls and an upper condition of a free deformable surface. The stationary linear problem was first investigated by V. Kh. Izakson and V. I. Yudovich in 1968 and their work is reviewed in [19]. The stationary problem with rotation and a variety of thermal boundary conditions was investigated by Dávalos-Orozco and López-Mariscal [13]. The problem of oscillating convection was first investigated by Benguria and Depassier [2]. They found that when the wall is solid, due to the restriction R/PrG < 1 (discussed presently) it is not possible for oscillatory convection to have a smaller critical Rayleigh number than that of stationary convection with surface deformation. Therefore, only the free wall case presents oscillatory convection as the first unstable one. G = gd 3 /ν 2 is the Galileo number, representative of the surface restoration due to the gravitational force. The deformation of the surface is due to a pressure difference which is balanced by the shear stresses at the fluid surface, that is, the dimensional zero stress jump at the surface When the surface is flat the pressure condition is p * − p ∞ = 0 (no pressure jump), where p ∞ is the pressure of the ambient gas whose viscosity is neglected. This problem requires the kinematic boundary condition of the surface deformation which in two-dimensions and in non dimensional form is which means that a fluid particle remains on the fluid interface and that the time variation of the convected surface deflection moves with the same vertical velocity as that of the fluid particle. η is an extra dependent variable representing the amplitude of the free surface deformation. The non dimensional normal and tangent vectors are defined as where N = (η 2 x + 1) 1/2 and the subindexes mean partial derivative. Other conditions, like Equations 18 to 20, defined using the normal and tangential vectors to the free deformable surface have to be modified with the definitions given in Equations 67 and 68. Equation 66 has to be multiplied by n to obtain the normal stress boundary condition and by t to calculate the tangential stress boundary condition. Here, the problem is assumed two-dimensional, but when it is three-dimensional it is necessary to define another tangential vector to the surface, perpendicular to both n and t. The problem simplifies if z 1 = −1 and the boundary conditions at the free surface are set at z = η (x, t). For the linear problem η (x, t) is assumed small. This gives the possibility to make a Taylor expansion of the variables at the free surface, that is, around z = 0, and simplify the boundary conditions. From this expansion, it is found that the fixed heat flux and the fixed temperature conditions remain the same in the linear problem. In two-dimensional flow it is possible to use the stream function. In this case, as a result of the expansion and use of normal modes, the kinematic, tangential stress and normal stress boundary conditions, respectively, are Ψ is the amplitude of the stream function. The approximations done here are only valid when the Galileo number satisfies R/PrG < 1. The reason is that the density and the temperature perturbations, related by ρ = (R/PrG)θ, should be ρ < θ in order to satisfy the Boussinesq approximation, which, among other things, neglects the variation of density with temperature in front of the inertial term of the balance of momentum equation. Consequently, the approximation is valid when the critical Rayleigh numbers satisfy R C < PrG. In the stationary problem the new parameter is in fact PrG due to the condition Equation 71 (see [19], [13]). However, in the oscillatory problem, Pr is an independent parameter, as seen in Equation 48, but it appears again in the condition Equation 71. Notice that in the limit of PrG → ∞ condition Equation 71 reduces to that of a flat wall. Then, the product PrG has two effects when it is large: 1) it works to guarantee the validity of the Boussinesq approximation under free surface deformation and 2) it works to flatten the free surface deformation.
The problem of oscillatory convection was solved analytically by Benguria and Depassier [2] when it occurs before stationary convection, that is, when the flat wall is free with fixed temperature and the upper deformable surface has fixed heat flux. They found that the cells are very large and took the small wavenumber limit. The critical Rayleigh number is R c = 30 and k c = 0.
Nonlinear waves for the same case of the linear problem, have been investigated by Aspe and Depassier [1] and by Depassier [15]. In the first paper, surface solitary waves of the Korteweg-de Vries (KdV) type were found. In the other one, Depassier found a perturbed Boussinesq evolution equation to describe bidirectional surface waves.

Natural convection in second-order fluids
The methods used in natural convection of Newtonian fluids can also be used in non Newtonian flows. Linear and non linear natural convection of second order fluids was investigated by Siddheshwar and Sri Krishna [58]. They assumed the flow is two-dimensional and used the free-free and fixed temperature boundary conditions. Here, the constitutive Equation 5 is used in the balance of momentum equation.
For the linear problem use is made of normal modes. The amplitudes of the stream function and the temperature are assumed of the form sin nπz which satisfies the boundary conditions. As before, the substitution in the governing equations leads to the formula for the marginal Rayleigh number of a second order fluid Here, Q = γ/ρ 0 d 2 is the elastic parameter. The Rayleigh number is real and the imaginary part should be zero. The only way this is possible is that ω = 0. Therefore, the flow can not be oscillatory and this reduces to Equation 52 for the marginal Rayleigh number and to the critical one of a Newtonian fluid for free-free convection, that is, R C = 27π 4 /4 at k C = π/ √ 2. By means of the energy method for the linear problem, Stastna [60] has shown that, in the solid-solid case and fixed temperature at the walls, the convection can not be oscillatory and ω = 0. Therefore, again the linear critical Rayleigh number and wavenumber are the same as those of the Newtonian fluid.
In their paper, Siddheshwar and Sri Krishna [58] also investigated the possibility of chaotic behaviour to understand the role played by the elastic parameter Q. They use the doble Fourier series method of Veronis [65] to calculate, at third order, a nonlinear system of Lorenz equations [32] used to investigate possible chaotic behavior in convection. In particular, the form selected by Lorenz for the time dependent amplitudes of the stream function and the temperature are Ψ(x, z, t) = X(t) sin kx sin πz (73) Θ(x, z, t) = Y(t) cos kx sin πz + Z(t) sin 2πz which satisfy the boundary conditions. These are used in the equations to obtain the nonlinear coupled Lorenz system of equations for the amplitudes X(t), Y(t) and Z(t) where the q i (1 = 1, · · · , 7) are constants including parameters of the problem. According to [58] the elastic parameter Q appears in the denominator of the constants q 1 and q 2 . In the system of Equations 74 the variables can be scaled to reduce the number of parameters to three. The results of [58] show that random oscillations occur when the parameters Q and Pr are reduced in magnitude. Besides, they found the possibility that the convection becomes chaotic for the magnitudes of the parameters investigated.
When the walls are solid and the heat flux is fixed, results of the nonlinear convective behaviour of a second-order fluid have been obtained by Dávalos and Manero [12]. They used the method of Chapman and Proctor [8] to calculate the evolution equation that describes the instability. A small wavenumber expansion is done like that of Equation 61 and62. However, in contrast to [8], their interest was in a three-dimensional problem and instead of using the stream function, use was made of a function defined by u = ∇ × ∇ × φk. The boundary condition of this function φ at the walls is φ = 0. It is found, by means of the solvability condition at first order, that the critical Rayleigh number is the same as that of the Newtonian case, that is, R C = 720 at k C = 0 and that convection can not be oscillatory. At the next approximation level, the solvability condition gives the evolution equation for the nonlinear convection. The result was surprising, because it was found that the evolution equation is exactly the same as that of the Newtonian fluid, that is, Equation 63 but in three-dimensions. The only difference with Newtonian fluids will be the friction on the walls due to the second order fluid. As explained above, the flow under the fixed heat flux boundary condition is very slow and the convection cell is very large. Therefore, this result may be related with the theorems of Giesekus, Tanner and Huilgol [20,24,63] which say that the creeping flow of a second-order incompressible fluid, under well defined boundary conditions, is kinematically identical to the creeping flow of a Newtonian fluid. The results presented here are a generalization of those theorems for three-dimensional natural convection evolving in time.

Natural convection in Maxwell fluids
In order to investigate the convection of a Maxwell fluid Equations 12 and 29 have to be used in the balance of momentum equation. The linear problem was investigated by Vest and Arpaci [66] for both free and solid walls and by Sokolov and Tanner [59] who present an integral model for the shear stress tensor which represent a number of non Newtonian fluids, including that of Maxwell. The linear equations in two-dimension use the stream function. In normal modes they are expressed as: where N = (1 + iωL) is a complex constant which depends on the viscoelastic relaxation time and the frequency of oscillation.
The combination of these two equations gives: The free-free linear problem of [66] is illustrative. Assuming that the amplitudes are proportional to sin nπz, the Equation 77 is transformed into a complex algebraic equation where In Equation 75 the real and imaginary parts have to be zero. Then There are two possibilities. 1) From the first ω = 0 and from the second A 3 = 0. 2) From the first ω = 0 with ω 2 = A 2 and from the second, after substitution of ω, From 1), A 3 gives the marginal stationary Rayleigh number for different modes n. The critical value for mode n = 1 has already been given above. From 2), the marginal Rayleigh number for oscillatory convection can be calculated for different modes n. Vest and Arpaci show that, when the relaxation time parameter has a magnitude large enough, it is possible to have oscillatory convection as the first unstable one, with R C smaller than that of stationary convection. Also, they show that an increase of Pr decreases considerably R C making convection very unstable. A similar behavior at criticality can be found for the solid-solid case. However, the solution is far more complex because it has to be solved numerically ensuring that the proper value of the Rayleigh number is real. The frequency is obtained form the numerical solution of the imaginary part and it is substituted into the real part, which still contains the frequency. The marginal Rayleigh number is obtained from the real part.
Variation of the wavenumber leads to the minimum of the Rayleigh number, the critical one, with its corresponding wavenumber and frequency. The conclusions obtained by Vest and Arpaci [66] are that the solid-solid case is more stable than the free-free case but qualitatively the response to the parameters L and Pr is similar.
The problem of a viscoelastic fluid layer with free and deformable surface will be discussed in detail in the section for Oldroyd fluid convection. The Maxwell fluid case is included in that problem.
The effect of the thickness and thermal conductivity of the walls on linear convection of a Maxwell fluid layer was investigated by Pérez-Reyes and Dávalos-Orozco [46]. They found that, by making some algebra, those effects can be included in a kind of Biot number which appears in the thermal boundary conditions of the upper and lower walls. The difference with respect to the Newtonian problem is that here it is necessary to calculated numerically the frequency of oscillation in the same way as explained in the last paragraph for the solid-solid case. The number of parameters in the equations increased and therefore the ratio κ/κ W which appears in the heat diffusion Equation 26 of the walls is assumed equal to one. Besides, it is supposed that the ratios of conductivities and thicknesses of the upper and lower walls are the same. With this in mind, some results are presented in Fig. 2.
Notice that in the figures F is used instead of L, for the relaxation time, and that D is used instead of d. Note that here the curves of R C increase with X instead of decreasing as in the Newtonian case for both magnitudes of D. Then, in this case it is easier to reach a codimension-two point where stationary convection competes with oscillatory convection to be the first unstable one. In the figure, the dashed lines correspond to stationary convection. This codimension-two point depends on the Prandtl number. When Pr increases there is a magnitude after which oscillatory convection is always the first unstable one (see [46]). In contrast, for other magnitudes of the relaxation time the behavior is similar to that of the Newtonian fluid (see Fig. 1) but with far more smaller magnitudes of R C , as seen in Fig. 3.
Thus, the oscillatory flow is very unstable for the new magnitude of F = 100. It is of interest to observe the different reaction the flow instability has with respect to X and D for both magnitudes of F. Nonlinear convection of a Maxwell fluid was investigated by Van Der Borght [64] using the ideal free-free boundary conditions. They calculated the heat dissipation of nonlinear stationary hexagonal convection cells by means of the Nusselt number. It was found that, for a given Rayleigh number, viscoelasticity effects only produce a slightly higher Nusselt number than Newtonian convection. Nonlinear traveling waves in a Maxwell fluid were investigated by Brand and Zielinska [5] using free-free boundary conditions. They obtain one Landau equation for stationary convection and other one for oscillatory convection (see Equation1). They found that standing waves are preferred over traveling waves for Pr < 2.82 at a codimension-two point. They also investigated the wave modulation in space by means of an equation similar to Equation 2 but of higher nonlinear order. Irregular and sensitive to initial conditions behavior of a convecting Maxwell fluid was investigated by Khayat [29]  using the free-free boundary conditions. The variables are written in the form of Equations 73, in addition to those of the shear and the two normal stresses. In this way, instead of three coupled Lorenz Equations 74, he obtains a system of four coupled equations which include as new time dependent variable, a linear combination of the amplitude of the normal stresses difference and the shear stress. In the limit L → 0, the new system reduces to that of Lorenz. He found that above a critical magnitude of the relaxation time L C the flow is oscillatory. For an L below the critical one, the route to chaotic motion is different from that of a Newtonian fluid, even in the case where L is very small. Viscoelasticity produces chaos when the Newtonian fluid still is non chaotic.

Natural convection in Oldroyd fluids
The Maxwell model of viscoelastic fluids presented above, shows an extreme (violent) mechanical behavior in comparison to other models. Mainly, this occurs when the relaxation time is large and the flow shows a more elastic behavior. The most elementary correction to this model is the Oldroyd fluid model which includes the property of shear rate retardation. That means that the fluid motion relaxes for a time interval even after the shear stresses have been removed. The representative magnitude of that time interval is called retardation time. This characteristic moderates the mechanical behavior of the Oldroyd fluid. The equations of motion require the constitutive Equation 13 or 30 (in non dimensional form). The ratio of the retardation time over the relaxation time E appears as an extra parameter which satisfies E < 1 (see [4], page 352). Notice that when E = 0 the Maxwell constitutive equation is recovered. Therefore, for very small E the behavior of the Oldroyd fluid will be similar to that of the Maxwell one.
The convection with free-free boundary conditions was investigated by Green [21] who obtained an equation similar to Equation 78. In the same way, from the solution of the real and imaginary parts, it is possible to calculate the marginal Rayleigh number and frequency of oscillation. In this case it is also possible to find a magnitude of L and E where the stationary and oscillatory convection have the same Rayleigh numbers, the codimension-two point. The Prandtl number plays an important role in this competition to be the first unstable one. The solid-solid problem was solved numerically by Takashima [62]. He shows that the critical Rayleigh number is decreased by an increase of L and increased by an increase of E. The numerical results show that an increase of Pr decreases drastically the magnitude of R C for the Maxwell fluid, but it is not very important when E > 0. Oscillatory convection is the first unstable one after a critical magnitude L C is reached, which depends on the values of Pr and E. For small Pr, L C is almost the same for any E. However, for large Pr the L C for the Maxwell fluid is notably smaller than that of the Oldroyd fluid. This fluid was also investigated by Sokolov and Tanner [59] using an integral model. In contrast to the papers presented above, Kolkka and Ierley [30] present results including the fixed heat flux boundary condition. They also give some corrections to the results of Vest and Arpaci [66] and Sokolov and Tanner [59]. The qualitative behavior of convection with fixed heat flux is the same, for both free-free and solid-solid boundaries, but with important differences in the magnitude of R C . The codimension-two point also occurs for different parameters. Interesting results have been obtained by Martínez-Mardones and Pérez-García [35] who calculated the codimension-two points with respect to L and E for both the free-free and the solid-solid boundary conditions. Besides, they calculated the dependance these points have on the Prandtl number. They show that for fixed E, the L C of codimension-two point decreases with Pr.
When natural convection occurs with an upper free surface it is every day experience to see that the free surface is deforming due to the impulse of the motion of the liquid coming in the upward direction. The assumption that the free surface is deformable in the linear convection of an Oldroyd viscoelastic fluid was first investigated by Dávalos-Orozco and Vázquez-Luis [14]. Under this new condition, the description of linear convection needs the same Equations 75 to 77 and N = (1 + iωL)/(1 + iωLE). However, the free boundary conditions have to change because the surface deformation produces new viscous stresses due to viscoelasticity. The problem is assumed two-dimensional and the stream function appears in the boundary conditions of the upper free deformable surface. The mechanical boundary conditions Equations 69 and 70 are the same. However, the normal stress boundary condition Equation 71 changes into which includes the viscoelastic factor N. Note that here the reference frame locates the free surface at z = 0, that is z 1 = −1. The advantages of doing this were explained above. The thermal boundary conditions remain the same. Numerical calculations were done for free and solid lower walls. In both cases, the fixed temperature boundary condition was used in the lower wall and the fixed heat flux boundary condition was used in the upper surface. In all the calculations the Prandtl number was set equal to Pr = 1. The goal was to compare with the paper by Benguria and Depassier [2]. Under these conditions, the results were first compared with those of the Newtonian flat surface solid-free convection (R CS = 669, k CS = 2.09 see [2]), the Newtonian deformable surface oscillatory solid-free convection (R ON = 390.8, k CS = 1.76 see [2]) and with the viscoelastic (Oldroyd and Maxwell) flat surface solid-free convection (presented in the figures with dashed lines). The results were calculated for different Galileo numbers G. However, here only some sample results are presented (see [14] for more details). Then, inside this rage, viscoelastic convection with deformable surface is more stable than that of the Newtonian convection with deformable surface. The result was verified with different numerical methods (see [14]). This phenomenon was explained by means of the double role played by the Galileo number as an external body force on convection (like rotation, see [13]) and as restoring force of the surface deformation.
In Figure 5 shown are the results of the free-free case for G = 1000, which represents a larger restoration force of the surface deformation. This results were compared with those of the Newtonian flat surface free-free convection (R CS = 384.7, k CS = 1.76 see [2]), the Newtonian deformable surface oscillatory free-free convection (R ON = 30.0, k CS = 0.0 see [2]) and with the viscoelastic (Oldroyd and Maxwell) flat surface free-free convection (presented in the figures with dashed lines). The behavior of the curves is similar but, except for very small E (nearly Maxwell fluid) and large L, it is found that G has no influence on the instability of the free-free convection under the present conditions. The curve of the Maxwell fluid is the more unstable. The jumps found in the curves of k C are also explained due to the dual role played by G on the instability. The results presented for the solid-free and free-free boundary conditions show the importance that free surface deformation has on the natural convection instability of viscoelastic fluids. The nonlinear problem for an Oldroyd fluid was investigated by Rosenblat [55] for free-free boundary conditions. He found results for the three time derivatives Equations 6, 8 and 9. The weakly nonlinear approximation is used where the Rayleigh number is very near to the critical one. He found conditions for subcriticality and supercriticality calculating an algebraic quantity K which includes non dimensional parameters of all the fluids investigated. The conclusion for stationary convection is that the corotational Oldroyd model has subcritical bifurcation (and therefore is unstable) and that the upper and lower convected Oldroyd models can not have this bifurcation and are stable, as the Newtonian model. For oscillatory convection the problem is more complex and it is resolved numerically with plots of L vs E. However, the results are reviewed as follows. The corotational Oldroyd model has supercritical bifurcation and is stable. The upper and lower convected Oldroyd models have subcritical bifurcation and are unstable. A system of four coupled differential equations is proposed to investigate chaotic behavior which generalize the Lorenz system of Equations 74. He found the possibility of chaotic behavior. The solid-solid and solid-free boundary conditions were used by Park and Lee [43,44]. Important results are that the amplitude of convection and heat dissipation increase with L and for E small. The rigid walls cause more easily the subcritical bifurcation than free walls for the same viscoelastic parameters.
They conclude that Oldroyd fluid viscoelastic convection is characterized both by a Hopf bifurcation (for very large value of L) and a subcritical bifurcation.
Analytical and numerical methods were used by Martínez-Mardones et al. [36] to calculate the nonlinear critical parameters which lead to stationary convection as well as traveling and standing waves. By means of coupled Landau amplitude equations Martínez-Mardones et al. [37] investigated the pattern selection in terms of the viscoelastic parameters. They fix Pr and E and show that increasing L stationary convection changes into standing waves by means of a subcritical bifurcation. The convective and absolute instabilities for the three model time derivatives of the Oldroyd fluid were investigated by Martínez-Mardones et al. [38]. If the group velocity is zero at say k = k 0 and the real part of σ, s, in Equation 39 is positive, it is said that the instability is absolute. In this case, the perturbations grow with time at a fixed point in space. If the perturbations are carried away from the initial point and at that point the perturbation decays with time, the instability is called convective. By means of coupled complex Ginzburg-Landau equations (Landau equations with second derivatives in space and complex coefficients) they investigated problems for which oscillatory convection appears first. Besides, they investigated the effect the group velocity has on oscillatory convection. It is found that the conductive state of the fluid layer is absolutely unstable if L > 0 or E > E C and therefore, when 0 < E < E C , the state is convectively unstable. They also show that there is no traveling wave phenomena when passing from stationary convection to standing waves.

Conclusions
In this chapter many phenomena have been discussed in order to show the variety of problems which can be found in natural convection of Newtonian and viscoelastic fluids. One of the goals was to show that the different boundary conditions may give results which differ considerably from each other. Sometimes, the results are qualitatively the same and this is taken as an advantage to solve "simpler" problems as those corresponding to the linear and nonlinear equations with free-free boundary conditions. A change in the setting of the problem may produce large complications, as in the case of the free-free boundary conditions, but with one of them being deformable. In this case a new parameter appears, the Galileo number G, which complicates not only the number of numerical calculations, but also the physical interpretation of the results, as explained above. As have been shown, the introduction of viscoelasticity complicates even more the physics of convection. Depending on the boundary conditions, there can be stationary and oscillatory cells in linear convection. Nonlinear convection can be stationary but for other magnitudes of the parameters, traveling and standing waves may appear as the stable fluid motion. The problem is to find the conditions and magnitudes of the viscoelastic parameters when a particular convection phenomenon occurs. This is the thrilling part of viscoelastic convection. It is the hope of the present author that this review may motivate a number of readers to work in this rich area of research.