About this Journal Submit a Manuscript Table of Contents
ISRN Thermodynamics
Volume 2012 (2012), Article ID 614712, 10 pages
Research Article

Rayleigh-Bénard Convection of Non-Newtonian Power-Law Fluids with Temperature-Dependent Viscosity

1Laboratory of Flows and Transfers Modelling (LAMET), Physics Department, Faculty of Sciences and Technologies, Sultan Moulay Slimane University, BP 523, Beni-Mellal, Morocco
2Laboratory of Fluid Mechanics and Energetics (LMFE), Physics Department, Faculty of Sciences Semlalia, Cadi Ayyad University, BP 2390, Marrakech, Morocco

Received 30 September 2012; Accepted 16 October 2012

Academic Editors: G. L. Aranovich, P. Espeau, and H. Hirao

Copyright © 2012 Mourad Kaddiri et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Two-dimensional steady-state Rayleigh-Bénard convection of thermodependent power-law fluids confined in a square cavity, heated from the bottom and cooled on the top with uniform heat fluxes, has been conducted numerically using a finite difference technique. The effects of the governing parameters, which are the Pearson number , the flow behaviour index , and the Rayleigh number , on the flow onset, flow structure, and heat transfer have been examined. The heatlines concept has been used to explain the heat transfer deterioration due to temperature-dependent viscosity effect that m expresses.

1. Introduction

Rayleigh-Bénard convection is a flow that occurs beyond a critical temperature gradient, called convection threshold, when a fluid layer is bounded below by a hot wall and above by a cold one. This phenomenon, which consists in general of counterrotating cells, whose number depends on the layer aspect ratio, has been the subject of many investigations in the past, due to its importance in various physical situations, natural or industrial, such as geophysics, astrophysics, meteorology, thermal storage, crystal growth, thermal comfort in buildings, and electronic cooling [1]. In addition to that, it is particularly suitable, from a scientific point of view, for studying and understanding dynamic instability, bifurcation, or chaotic behavior in fluids. The research undertaken in this area, which is of experimental, numerical, or theoretical nature, covers both Newtonian [2, 3] and non-Newtonian [411] fluids, but with an abundance of works for the former although the latter are frequently encountered in the engineering practice.

It is, however, to note that, in most of studies dealing with Rayleigh-Bénard convection, the fluid properties are customarily considered as constant with respect to temperature variations, except for density in the buoyancy force term of the momentum equation which obeys the Boussinesq approximation. According to Gray and Giorgini [12], this is valid only for a small temperature difference between the thermally active walls. For many liquids, such as petroleum oils, glycerine, silicone fluid, and some molten salts, the viscosity variation with temperature is often much stronger than that of the other properties; for example, in earth mantle and magma chamber temperature differences bring viscosity variations of up to 1030 [13]. This means that at large temperature differences the above-mentioned assumption could fail and lead to unrealistic predictions, as proved experimentally by Scirocco et al. [14] and numerically by Shin and Cho [15], whose results, of local Nusselt numbers for a polyacrylamide (Separan AP-273) solution, show 70%–300% heat transfer enhancement over those of a constant-viscosity fluid. In fact, according to Scirocco et al. [14], the viscosity decreases notably near the hot wall leading to an increase of the axial velocity in this region, through temperature-viscosity coupling, which results in a reduction of the thermal boundary layer thickness and, consequently, in a heat transfer improvement. Therefore, to be realistic, the dependency of the viscosity on the temperature (thermodependency in other words) is to be included in the flow and heat transfer analysis.

Some authors, namely, Stengel et al. [16], Richter et al. [17], Busse and Frick [18], White [19], Bottaro et al. [20], Davaille and Jaupart [21], Zaranek and Parmentier [22], Kim and Choi [23], and many others, have considered the thermodependency effect in their investigations where the fluid behaviour is Newtonian, but to our best knowledge, only very few studies on Rayleigh-Bénard convection in non-Newtonian fluids have taken into account such an effect. Among them is the relatively recent study conducted by Solomatov and Barr [24, 25] who examined numerically such an effect, for a power-law non-Newtonian fluids, on the onset of the Rayleigh-Bénard convection. Using the Frank-Kamenetskii approximation [26], for the variation of the viscosity with temperature, these authors investigated how the shape of the perturbation affects the onset of convection and suggested preliminary constraints on the dependence of the critical Rayleigh number on the amplitude of the most effective perturbations for a broad range of viscosity parameters. Although they believe that their study has shed some light upon the finite-amplitude onset of convection for power-law viscosity fluids, the problem is far from being well understood and needs to be investigated in the future both numerically and theoretically.

Therefore, in order to contribute to a better understanding of the thermodependence effect on buoyancy-driven convection in confined non-Newtonian media, a numerical study is carried out to investigate Rayleigh-Bénard convection flow and heat transfer in a square cavity containing non-Newtonian power-law fluids, whose effective viscosity depends strongly on temperature. Although the present work deals with such kinds of fluids, it is different from that of Solomatov and Barr [24, 25] insofar as it considers boundary conditions of Newman type by imposing uniform heat fluxes at the horizontal sides, situation that hasnever been consideredbefore.

2. Mathematical Formulation

2.1. Problem Statement and Viscosity Model

The geometry under consideration is sketched in Figure 1. It consists of a two-dimensional square enclosure of size . The vertical walls are thermally insulated while the horizontal ones are subjected to uniform density of heat flux , with the bottom heated and the top cooled. The non-Newtonian fluids considered here are those whose rheological behaviour can be represented by the power-law model, due to Ostwald-de Waele, which, in terms of laminar effective viscosity, can be expressed as: where the two empirical parameters, and , stand for the flow behavior and consistency indices, respectively. They are, in general, functions of the temperature, but in most of cases the temperature dependence of is very weak compared to that of [14], which can be described by the following Frank-Kamenetskii exponential law [26]: where is an exponent related to the flow energy activation and the universal gas constant, called the thermodependency coefficient and is the consistency index at the reference temperature . Note that (2) is an approximation of the original Arrhenius exponential law and reflects the viscosity decrease with the temperature.

Figure 1: Schematic view of the geometry and coordinates system.

On the other hand, for the behavior is Newtonian and the consistency is just the viscosity. For , the effective viscosity decreases with the shear rate and the behaviour is shear thinning (or pseudoplastic). Conversely, for , the viscosity increases with the shear rate and the behaviour is shear thickening (or dilatant).

Although the empirical power-law model does not converge towards a Newtonian behaviour in the limit of zero and maximum shear rates, it presents the advantage of being simple and mathematically tractable. Besides, it is a useful fit to the observed data and can often provide good quantitative results over a relatively large range of the shear rate, which justifies its use in most theoretical investigations and engineering purposes.

2.2. Governing Equations and Boundary Conditions

On the basis of the assumptions commonly adopted in natural convection problems and using the characteristic scales , , , , and , which correspond respectively, to length, time, velocity, vorticity, temperature, and stream function, the dimensionless governing equations for Boussinesq-temperature-dependent viscosity fluids, written in terms of vorticity, , temperature, , and stream function, , are as follows: where

2.3. Boundary Conditions

For the present problem, the appropriate nondimensional boundary conditions are

Note that such a formulation (, , and ) reduces the number of equations and unknowns by eliminating the pressure, which is without interest in the considered problem, but it presents the disadvantage of unknown at the boundaries, difficulty that has been overcome by adopting the Woods approximation [27].

2.4. Governing Parameters

In addition to the flow behaviour index, , three other dimensionless parameters appear in the governing equations. These are the Pearson, Prandtl, and Rayleigh numbers defined, respectively, as

The Pearson number [28], which is a new dimensionless quantity taking place in this study, measures the effect of temperature change on the effective viscosity.

2.5. Heat Transfer

The steady solution has been used to calculate the local Nusselt number defined as and its mean value expressed as follows:

2.6. Heatlines Formulation

The visualization of the paths followed by the heat flow through the enclosure requires the use of the heatlines concept, which consists of lines of constant heat function, , that are defined, according to Kimura and Bejan [29], from the following equations: whose derivation, with respect to and , and combination give rise to

To obtain the boundary conditions associated with (14), an integration of (13), along the four cavity walls, is necessary, which gives

Finally, the solution of (14) yields the values of , in the computational domain, whose contour plots provide the heatline patterns. Note that only the differences between the values of are required instead of its intrinsic ones, which offers the possibility to choose as an arbitrary reference value for .

3. Numerics

The two-dimensional governing equations have been discretised using the second order central finite difference methodology with a regular mesh size. The integration of (3) and (4) has been performed with the alternating direction implicit method (ADI). To satisfy the mass conservation, (5) has been solved by a point successive over relaxation method (PSOR) with an optimum relaxation factor calculated by the Frankel formula [27]. A grid of has been required for obtaining adequate results since, as shown in Table 1, a refinement to leads to maximum differences of 1.32% and 1.52% in terms of and , respectively. At each time step, , which has been chosen between and (depending on the values of and ), the convergence criterion has been satisfied for , where is the value of the stream function at the node for the kth iteration level.

Table 1: Preliminary tests on the grid size effect (, ).

With the Ostwald power-law model, the dimensionless viscosity (7) tends towards infinity at the cavity corners (where the velocity gradients are nil) for , which has made impossible direct numerical computations. However, this difficulty has been surmounted by using average values for the viscosity at the corners, which has led to possible and stable computations.

As an additional check of the results accuracy, an energy balance has been systematically verified for the system at each numerical code running. Thus, the overall heat transfer, through each horizontal plane has been evaluated and compared with the quantity of heat furnished to the system through the bottom wall () and leaving it through the top one (). For all the results reported here, the energy balance has been satisfied within 0.1% as a maximum difference.

Moreover, the present computational code has been validated against Benchmark solutions obtained, in the case of the classical Rayleigh-Bénard problem in square cavity, by Ouertatani et al. [30]. Comparative results are summarized in Table 2 where it can be seen that the agreement is quite good; the maximum difference does not exceed 0.6%.

Table 2: Validation of the numerical code for Pr = 0.71.

The validation of the present computational code is verified against the existing results for Rayleigh-Bénard convection in a square cavity filled with air [30] and is shown in Table 2.

4. Results and Discussion

In this section are presented and discussed the results corresponding to the effects of the three dimensionless parameters defined previously, namely, the Pearson number, , the flow behaviour index, , and the Rayleigh number, , that vary in the ranges , , and , respectively. Note that the values of are selected so that they can cover the shear-thinning , Newtonian , and shear-thickening cases. The role of the Prandtl number, , is still ignored since, for the non-Newtonian fluids considered in this study, is much larger than unity [31], and an increase of this parameter, in such a situation, makes the contribution of the inertia terms (on the left-hand side of (3)) negligible and the fluid flow and heat transfer unchanged, as was reported by Ozoe and Churchill [4], Lamsaadi et al. [8], and many others. Therefore, the real parameters governing the study are , , and , and all the corresponding results are valid in the limit of .

4.1. Onset of Buoyancy-Driven Convection

In Figure 2 is displayed the critical Rayleigh number, , that is, the threshold of the convection onset, versus for various values of . This figure shows that, for any value of , is an increasing function of , which means that the shear-thinning behaviour anticipates the convection, while the shear-thickening one produces the opposite effect. Moreover, it can be seen that, for a given , an increase of leads to a decrease of , expressing the fact that the thermodependency tends to make precocious the convection. However, this effect gets less and less pronounced with a decreasing . This seems obvious, since the shear-thinning behaviour tends to reduce the effective viscosity and, therefore, to limit its variation with temperature (reduction of the thermodependency effect). In fact, with less viscous fluid, such as water, the viscosity is weakly temperature dependent, whilst with more viscous fluid, such as honey, the viscosity is highly temperature dependent.

Figure 2: Evolution of with , for two values of .
4.2. Analysis of Streamlines, Isotherms, and Heatlines Patterns

Let us, now, examine the explicit effect of a temperature-dependent viscosity on mass flow, temperature, and heat flow fields. With the assumed form of (7), numerical solutions, presented in terms of streamlines (a), isotherms (b), and heatlines (c) (Figures 3, 4, and 5), have been generated using and 10, which correspond to constant and temperature-dependent viscosity, respectively, for various values of and , which have been chosen dependently on those assigned to , since changes with this parameter.

Figure 3: Streamlines (a), isotherms (b), and heatlines (c) for , (solid line), and (dashdot line).
Figure 4: Streamlines (a), isotherms (b), and heatlines (c) for , (solid line), and (dashdot line).
Figure 5: Streamlines (a), isotherms (b), and heatlines (c) for , (solid line), and (dashdot line).

As can be seen, the flow remains, in general, unicellular and clockwise, but loses its symmetry with an increasing , in the limit of the explored values of and . In fact, in the case of , the flow spreads to the whole region of , but for , the streamlines are crowded in the region neighbouring the left lower corner, which means that the flow is intensified as a result of the viscosity decrease in such a region. This gives rise to a rheological sublayer where the convection is concentrated (strongly developed). In contrast, in the zone near the cold upper and vertical right walls the flow is suppressed due to the increase of the viscosity, which leads to a stagnant region: this is the conductive lid regime. This phenomenon (effect) becomes pronounced with an increasing and a decreasing , which implies that temperature-dependent viscosity effect manifests itself only when convection is weak, which is quite obvious since strong convection, that occurs for small and high values of and , respectively, inhibits such an effect.

At the same time, compared to the isoviscous case (), the isotherms exhibit a strong distortion in the rheological sublayer confining convection for , which affects seriously their centrosymmetry, but tend to be less distorted (i.e., the temperature gradient tends to weaken) in the stagnant region especially near the cold wall where the viscosity is strong. This tendency takes importance with an increasing and a decreasing as explained previously.

On the other hand, in order to have a microscopic description of the heat transfer process, which is different from the conventional Nusselt number that describes macroscopically the heat transfer rate, a heatlines analysis is required. Hence, with comparison with the isoviscous case (), the heatlines corresponding to the case present more distortion, which indicates that the path followed by the heat flow to reach the cold wall is more complicated in the rheological sublayer, and, consequently, the heat transfer is expected to be deteriorated in such a situation. Note, here also, that an increase of and a decrease of encourage such a fact.

4.3. Flow Intensity and Mean Heat Transfer Rate

Figures 6 and 7 display the evolutions of the flow intensity, , and the mean heat transfer rate, , with , for the selected values of and . First of all, it is obvious to see that and are increasing functions of , as a result of the positive contribution of the buoyancy effects in promoting convection. In contrast, an increase of leads to the opposite effect, since it is well known that the shear-thinning behaviour enhances natural convection heat transfer, while the shear-thickening one reduces it. In addition to that, and are seen to increase and decrease with , respectively. These tendencies are related to the fact that with a temperature-dependent viscosity, the convective activities are facilitated near the hot wall (i.e., in the rheological sublayer) and suppressed out off this region, especially near the cold wall. Besides, the decrease of , expressing the degradation of the heat transfer, is expected from the analysis of the heatlines patterns.

Figure 6: Evolution of with for (solid line) and (dashdot line) and various values of .
Figure 7: Evolution of with for (solid line) and (dashdot line) and various values of .

5. Conclusion

A numerical investigation of steady buoyancy convection in a square cavity filled with temperature-dependent viscosity non-Newtonian power-law fluids has been performed. Both horizontal sides of the enclosure are submitted to uniform heat flux while its vertical ones are insulated. The exponential model, due to Frank-Kamenetski, for the viscosity variations with the temperature, has been adopted. The study has been focused particularly on the effects of the thermodependency of the viscosity on the flow and thermal fields and heat transfer.

It emerges from such a study that the viscosity variations with temperature act to reduce the convective zone thickness, since the circulation remains confined in the region of low viscosity, while in the layer of high viscosity the flow slow downs giving rise to a conductive lid regime that reduces notably the heat transfer, depending on the articulated effects of the rheological behaviour and the buoyancy forces, as confirmed by the heatlines analysis.

The new model for Rayleigh-Benard convection, developed and analysed in this study, can find application, for instance, in geophysics (convection in robe). In addition it may provide an insight into the physics of the three-dimensional problem.


Thermodependency coefficient
:Acceleration due to gravity
:Heat function
:Height or width of the enclosure
:Consistency index for a power-law fluid at the reference temperature
:Pearson number
:Flow behaviour index for a power-law fluid at the reference temperature
:Local Nusselt number
:Mean Nusselt number
Pr:Generalised Prandtl number
:Constant density of heat flux
:Generalised Rayleigh number
:Dimensionless temperature
:Reference temperature at the centre of the enclosure
:Characteristic temperature
:Dimensionless horizontal and vertical velocities
:Dimensionless horizontal and vertical coordinates.
Greek Symbols
:Thermal diffusivity of fluid at the reference temperature
:Thermal expansion coefficient of fluid at the reference temperature
:Thermal conductivity of fluid at the reference temperature
:Dynamic viscosity for a Newtonian fluid at the reference temperature
:Dimensionless effective viscosity of fluid
:Density of fluid at the reference temperature
:Dimensionless vorticity
:Dimensionless stream function.
:Dimensional variables.
:Effective variable
:Maximum value
:Reference value taken at the cavity centre.


  1. B. Gebhart, Y. Jaluria, R. L. Mahajan, and B. Sammakia, Buoyancy-Induced Flows and Transport, Hemisphere, Washington, DC, USA, 1985.
  2. S. W. Churchill, Heat Exchanger Design Handbook, section 2.5.8., Hemisphere Publishing Corporation, New York, NY, USA, 1983.
  3. A. Bejan, Convection Heat Transfer, John Wiley Sons, New York, NY, USA, 1984.
  4. H. Ozoe and S. W. Churchill, “Hydrodynamic stability and natural convection in Ostwald-De Waele and Ellis fluids: the development of a numerical solution,” American Institute of Chemical Engineers Journal, vol. 18, pp. 1196–1207, 1972.
  5. H. M. Park and D. H. Ryu, “Rayleigh-Bérnard convection of viscoelastic fluids in finited domains,” Journal of Non-Newtonian Fluid Mechanics, vol. 98, no. 2-3, pp. 169–184, 2001. View at Publisher · View at Google Scholar · View at Scopus
  6. M. Ohta, M. Ohta, M. Akiyoshi, and E. Obata, “A numerical study on natural convective heat transfer of pseudoplastic fluids in a square cavity,” Numerical Heat Transfer A, vol. 41, no. 4, pp. 357–372, 2002. View at Publisher · View at Google Scholar · View at Scopus
  7. H. Inaba, C. Dai, and A. Horibe, “Natural convection heat transfer of microemulsion phase-change-material slurry in rectangular cavities heated from below and cooled from above,” International Journal of Heat and Mass Transfer, vol. 46, no. 23, pp. 4427–4438, 2003. View at Publisher · View at Google Scholar · View at Scopus
  8. M. Lamsaadi, M. Naïmi, and M. Hasnaoui, “Natural convection of non-Newtonian power law fluids in a shallow horizontal rectangular cavity uniformly heated from below,” Heat and Mass Transfer, vol. 41, no. 3, pp. 239–249, 2005. View at Publisher · View at Google Scholar · View at Scopus
  9. J. Zhang, D. Vola, and I. A. Frigaard, “Yield stress effects on Rayleigh-Bénard convection,” Journal of Fluid Mechanics, vol. 566, pp. 389–419, 2006. View at Publisher · View at Google Scholar · View at Scopus
  10. N. J. Balmforth and A. C. Rust, “Weakly nonlinear viscoplastic convection,” Journal of Non-Newtonian Fluid Mechanics, vol. 158, no. 1–3, pp. 36–45, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. A. Vikhansky, “Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers,” Physics of Fluids, vol. 21, no. 10, Article ID 103103, 2009. View at Publisher · View at Google Scholar · View at Scopus
  12. D. D. Gray and A. Giorgini, “The validity of the boussinesq approximation for liquids and gases,” International Journal of Heat and Mass Transfer, vol. 19, no. 5, pp. 545–551, 1976. View at Scopus
  13. D. A. Yuen, W. R. Peltier, and G. Schubert, “On the existence of a second scale of convection in the upper mantle,” Geophysical Journal, Royal Astronomical Society, vol. 65, no. 1, pp. 171–190, 1981. View at Scopus
  14. V. Scirocco, R. Devienne, and M. Lebouché, “Ecoulement laminaire et transfert de chaleur pour un fluide pseudo-plastique dans la zone d'entrée d'un tube,” International Journal of Heat and Mass Transfer, vol. 28, no. 1, pp. 91–99, 1985. View at Scopus
  15. S. Shin and Y. I. Cho, “Laminar heat transfer in a rectangular duct with a non-Newtonian fluid with temperature-dependent viscosity,” International Journal of Heat and Mass Transfer, vol. 37, no. 1, pp. 19–30, 1994. View at Scopus
  16. K. C. Stengel, D. S. Oliver, and J. R. Booker, “Onset of convection in a variable-viscosity fluid,” Journal of Fluid Mechanics, vol. 120, pp. 411–431, 1982.
  17. F. M. Richter, H. C. Nataf, and S. Daly, “Heat transfer and horizontally averaged temperature of convection with large viscosity variations,” Fluid Mechanics, vol. 129, pp. 173–192, 1983.
  18. F. H. Busse and H. Frick, “Square-pattern convection in fluids with strongly temperature-dependent viscosity,” Journal of Fluid Mechanics, vol. 150, pp. 451–465, 1985. View at Scopus
  19. D. B. White, “The plan forms and onset of convection with a temperature-dependent viscosity,” Journal of Fluid Mechanics, vol. 191, pp. 247–286, 1988.
  20. A. Bottaro, P. Metzener, and M. Matalon, “Onset and two-dimensional patterns of convection with strongly temperature-dependent viscosity,” Physics of Fluids A, vol. 4, no. 4, pp. 655–663, 1992. View at Scopus
  21. A. Davaille and C. Jaupart, “Onset of thermal convection in fluids with temperature-dependent viscosity: application to the oceanic mantle,” Journal of Geophysical Research, vol. 99, no. 10, pp. 19–866, 1994. View at Scopus
  22. S. E. Zaranek and E. M. Parmentier, “The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres,” Earth and Planetary Science Letters, vol. 224, no. 3-4, pp. 371–386, 2004. View at Publisher · View at Google Scholar · View at Scopus
  23. M. C. Kim and C. K. Choi, “The onset of buoyancy-driven convection in fluid layers with temperature-dependent viscosity,” Physics of the Earth and Planetary Interiors, vol. 155, no. 1-2, pp. 42–47, 2006. View at Publisher · View at Google Scholar · View at Scopus
  24. V. S. Solomatov and A. C. Barr, “Onset of convection in fluids with strongly temperature-dependent, power-law viscosity,” Physics of the Earth and Planetary Interiors, vol. 155, no. 1-2, pp. 140–145, 2006. View at Publisher · View at Google Scholar · View at Scopus
  25. V. S. Solomatov and A. C. Barr, “Onset of convection in fluids with strongly temperature-dependent, power-law viscosity. 2. Dependence on the initial perturbation,” Physics of the Earth and Planetary Interiors, vol. 165, no. 1-2, pp. 1–13, 2007. View at Publisher · View at Google Scholar · View at Scopus
  26. N. J. Balmforth and A. Provenzale, Geophysical Aspects of Non-Newtonian Fluid Mechanics, vol. 582 of Liberal National Party, Springer, 2001.
  27. P. J. Roache, Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, NM, USA, 1982.
  28. C. Métivier and C. Nouar, “Linear stability of the Rayleigh-Bénard Poiseuille flow for thermodependent viscoplastic fluids,” Journal of Non-Newtonian Fluid Mechanics, vol. 163, no. 1–3, pp. 1–8, 2009. View at Publisher · View at Google Scholar · View at Scopus
  29. S. Kimura and A. Bejan, “The heatline visualization of convective heat transfer,” ASME Journal of Heat Transfer, vol. 105, pp. 916–919, 1983.
  30. N. Ouertatani, N. Ben Cheikh, B. Ben Beya, and T. Lili, “Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure,” Comptes Rendus, vol. 336, no. 5, pp. 464–470, 2008. View at Publisher · View at Google Scholar · View at Scopus
  31. Y. I. Cho and J. P. Harnett, “Non-newtonian fluids in circular pipe flow,” Advances in Heat Transfer, vol. 15, pp. 59–141, 1982. View at Publisher · View at Google Scholar · View at Scopus