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

Abstract

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 (0≤𝑚≤10), the flow behaviour index (0.6≤𝑛≤1.4), and the Rayleigh number (0<Ra≤10⁵), 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 [4–11] 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 10³⁰ [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: 𝜇^′_𝑎=𝑘_𝑇2𝜕𝑢^′𝜕𝑥^′2+𝜕𝑣^′𝜕𝑦^′2+𝜕𝑢^′𝜕𝑦^′+𝜕𝑣^′𝜕𝑥^{′2(𝑛𝑇−1)/2}, (1) 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 (𝑛_𝑇≈constant=𝑛) compared to that of 𝑘_𝑇 [14], which can be described by the following Frank-Kamenetskii exponential law [26]: 𝑘_𝑇=𝑘𝑒^{−𝑏(𝑇′−𝑇′𝑟)}, (2) 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𝑛=1 the behavior is Newtonian and the consistency is just the viscosity. For 0<𝑛<1, the effective viscosity decreases with the shear rate and the behaviour is shear thinning (or pseudoplastic). Conversely, for 𝑛>1, 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 𝐻^′, 𝐻^′2/𝛼, 𝛼/𝐻^′, 𝛼/𝐻^′2, 𝑞^′𝐻^′/𝜆 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: 𝜕Ω+𝜕()𝜕𝑡𝑢Ω+𝜕()𝜕𝑥𝑣Ω𝜇𝜕𝑦=Pr_𝑎𝜕²Ω𝜕𝑥²+𝜕²Ω𝜕𝑦²+2𝜕𝜇_𝑎𝜕𝑥𝜕Ω+𝜕𝑥𝜕𝜇_𝑎𝜕𝑦𝜕Ω𝜕𝑦+𝑆_Ω𝜕𝑇+𝜕()𝜕𝑡𝑢𝑇+𝜕()𝜕𝑥𝑣𝑇=𝜕𝜕𝑦²𝑇𝜕𝑥²+𝜕²𝑇𝜕𝑦²,𝜕²𝜓𝜕𝑥²+𝜕²𝜓𝜕𝑦²=−Ω, (3)(4)(5) where 𝑢=𝜕𝜓𝜕𝑦,𝑣=−𝜕𝜓𝜕𝑥;Ω=𝜕𝑣−𝜕𝑥𝜕𝑢,𝜇𝜕𝑦_𝑎=𝑒^−𝑚𝑇2𝜕𝑢𝜕𝑥²+𝜕𝑣𝜕𝑦²+𝜕𝑢+𝜕𝑦𝜕𝑣𝜕𝑥^{2(𝑛−1)/2},𝑆_Ω𝜕=Pr²𝜇_𝑎𝜕𝑥²−𝜕²𝜇_𝑎𝜕𝑦²𝜕𝑢+𝜕𝑦𝜕𝑣𝜕𝜕𝑥−2²𝜇_𝑎𝜕𝑥𝜕𝑦𝜕𝑢−𝜕𝑥𝜕𝑣𝜕𝑦+PrRa𝜕𝑇.𝜕𝑥 (6)(7)(8)

2.3. Boundary Conditions

For the present problem, the appropriate nondimensional boundary conditions are 𝑢=𝑣=𝜓=𝜕𝑇𝜕𝑥=0,for𝑥=0and1,𝑢=𝑣=𝜓=𝜕𝑇𝜕𝑦+1=0,for𝑦=0and1. (9)

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 𝑚=𝑏Δ𝑇^∗1=−𝑘_𝑇𝑑𝑘_𝑇𝑘𝑑𝑇=−𝑑ln_𝑇/𝑘,𝐻𝑑𝑇Pr=𝑘/𝜌^′2−2𝑛𝛼^2−𝑛,Ra=𝑔𝛽𝐻^′2𝑛+2𝑞^′𝑘/𝜌𝛼^𝑛𝜆. (10)

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 (𝑥)=𝑞Nu^′𝐻^′𝜆Δ𝑇^′=1=1Δ𝑇𝑇(𝑥,0)−𝑇(𝑥,1) (11) and its mean value expressed as follows: Nu=¹₀Nu(𝑥)𝑑𝑥. (12)

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: 𝜕𝐻𝜕𝑦=𝑢𝑇−𝜕𝑇𝜕𝑥;−𝜕𝐻𝜕𝑥=𝑣𝑇−𝜕𝑇,𝜕𝑦 (13) whose derivation, with respect to 𝑥 and 𝑦, and combination give rise to 𝜕²𝐻𝜕𝑥²+𝜕²𝐻𝜕𝑦²=−𝜕𝑣𝑇+𝜕𝑥𝜕𝑢𝑇.𝜕𝑦 (14)

To obtain the boundary conditions associated with (14), an integration of (13), along the four cavity walls, is necessary, which gives 𝐻()𝐻()()𝐻()𝐻()()0,𝑦=𝐻0,0,for𝑥=0,𝑥,1=𝐻0,1−𝑥,for𝑦=0,1,𝑦=𝐻1,1,for𝑥=1,𝑥,0=𝐻1,0+1−𝑥,for𝑦=0. (15)

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 𝐻(0,0)=0 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 81×81 has been required for obtaining adequate results since, as shown in Table 1, a refinement to 101×101 leads to maximum differences of 1.32% and 1.52% in terms of |𝜓_max| and Nu, respectively. At each time step, 𝛿𝑡, which has been chosen between 10⁻⁷ and 10⁻⁴ (depending on the values of 𝑛 and Ra), the convergence criterion ∑_𝑖,𝑗|𝜓^𝑘+1_𝑖,𝑗−𝜓^𝑘_𝑖,𝑗|/∑_𝑖,𝑗|𝜓^𝑘+1_𝑖,𝑗|<10⁻⁴ 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 (𝑚=5, Pr→∞).

With the Ostwald power-law model, the dimensionless viscosity (7) tends towards infinity at the cavity corners (where the velocity gradients are nil) for 0<𝑛<1, 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 (𝑦=0) and leaving it through the top one (𝑦=1). 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, Ra, that vary in the ranges 0≤𝑚≤10, 0.6≤𝑛≤1.4, and 0<Ra≤10⁵, respectively. Note that the values of 𝑛 are selected so that they can cover the shear-thinning (0<𝑛<1), Newtonian (𝑛=1), and shear-thickening (𝑛>1) cases. The role of the Prandtl number, Pr, is still ignored since, for the non-Newtonian fluids considered in this study, Pr 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 Ra, and all the corresponding results are valid in the limit of Pr→∞.

4.1. Onset of Buoyancy-Driven Convection

In Figure 2 is displayed the critical Rayleigh number, Ra_𝑐, that is, the threshold of the convection onset, versus 𝑛 for various values of 𝑚. This figure shows that, for any value of 𝑚, Ra_𝑐 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 Ra_𝑐, 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 Ra𝑐 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 𝑚=0 and 10, which correspond to constant and temperature-dependent viscosity, respectively, for various values of 𝑛 and Ra, which have been chosen dependently on those assigned to 𝑛, since Ra_𝑐 changes with this parameter.

(a)

(b)

(c)

Figure 3 

Streamlines (a), isotherms (b), and heatlines (c) for 𝑛=0.6, 𝑚=0 (solid line), and 𝑚=10 (dashdot line).

(a)

(b)

(c)

Figure 4 

Streamlines (a), isotherms (b), and heatlines (c) for 𝑛=1.0, 𝑚=0 (solid line), and 𝑚=10 (dashdot line).

(a)

(b)

(c)

Figure 5 

Streamlines (a), isotherms (b), and heatlines (c) for 𝑛=1.4, 𝑚=0 (solid line), and 𝑚=10 (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 Ra. In fact, in the case of 𝑚=0, the flow spreads to the whole region of 0<𝑦<1, but for 𝑚>0, 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 Ra, 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 Ra, respectively, inhibits such an effect.

At the same time, compared to the isoviscous case (𝑚=0), the isotherms exhibit a strong distortion in the rheological sublayer confining convection for 𝑚=10, 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 Ra 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 (𝑚=0), the heatlines corresponding to the case 𝑚=10 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 Ra encourage such a fact.

4.3. Flow Intensity and Mean Heat Transfer Rate

Figures 6 and 7 display the evolutions of the flow intensity, |𝜓_max|, and the mean heat transfer rate, Nu, with Ra, for the selected values of 𝑚 and 𝑛. First of all, it is obvious to see that |𝜓_max| and Nu are increasing functions of Ra, 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, |𝜓_max| and Nu 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 Nu, expressing the degradation of the heat transfer, is expected from the analysis of the heatlines patterns.

Figure 6 

Evolution of |𝜓max| with Ra for 𝑚=0 (solid line) and 𝑚=10 (dashdot line) and various values of 𝑛.

Figure 7 

Evolution of Nu with Ra for 𝑚=0 (solid line) and 𝑚=10 (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.

Nomenclature

𝑏∶	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
Nu:	Local Nusselt number
Nu:	Mean Nusselt number
Pr:	Generalised Prandtl number
𝑞′:	Constant density of heat flux
Ra:	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.

Superscript

′:	Dimensional variables.

Subscripts

𝑎:	Effective variable
max:	Maximum value
𝑟:	Reference value taken at the cavity centre.