Boundary Condition Effects on Laminar Natural Convection of Power

Dec 4, 2013 - Department of Mechanical Engineering, Karadeniz Technical University, Trabzon 61080 Turkey. ‡. School of Mechanical and Systems ...
1 downloads 0 Views 4MB Size
Article pubs.acs.org/IECR

Boundary Condition Effects on Laminar Natural Convection of Power-Law Fluids in a Square Enclosure Heated from below with Differentially Heated Horizontal Walls Osman Turan,† Frank Fotso-Choupe,‡ Jiawei Lai,‡ Robert J. Poole,§ and Nilanjan Chakraborty*,‡ †

Department of Mechanical Engineering, Karadeniz Technical University, Trabzon 61080 Turkey School of Mechanical and Systems Engineering, Newcastle University, Claremont Road, Newcastle-Upon-Tyne NE1 7RU, United Kingdom § School of Engineering, University of Liverpool, Brownlow Hill, Liverpool L69 3GH, United Kingdom ‡

ABSTRACT: Two-dimensional steady-state laminar natural convection of power-law fluids in square enclosures with differentially heated horizontal walls (heated from below) subjected to constant wall temperature (CWT) and constant wall heat flux (CHWF) boundary conditions has been analyzed in detail based on computational simulations and a detailed scaling analysis. The effects of power-law exponent n ranging from 0.6 to 1.8 on the thermal transport have been investigated for nominal values of Rayleigh number in the range 103−105 and a Prandtl number range of 10−105. It is found that the mean Nusselt number Nu increases with increasing (decreasing) values of Rayleigh number (power-law exponent) for both CWT and CWHF configurations because of the strengthening of convective transport. In contrast, the mean Nusselt number Nu remains insensitive to changes in Prandtl number Pr in the range 10−105. It has been found that the Nu values for the CWHF configuration remain smaller than the corresponding values in the case of CWT boundary condition (bc) for shear-thinning fluids, whereas Nu in the CWHF configuration for large values of n remains greater than the corresponding values in the case of CWT bc (for identical values of the power-law exponent and nominal Rayleigh and Prandtl numbers). We find that the steady two-dimensional convection in this configuration is realized in a narrower parameter range in the CWT bc than in the CWHF bc. Underpinned by a scaling analysis, new correlations have been proposed for the mean Nusselt number Nu for both CWT and CWHF boundary conditions and these correlations are shown to capture the computational results satisfactorily for the entire range of power-law exponents and nominal Rayleigh and Prandtl numbers considered here.

1. INTRODUCTION Natural convection within rectangular enclosures has been extensively analyzed in the literature for Newtonian fluids,1−3 and interested readers are referred to refs 2 and 3 for a detailed review. In comparison, relatively limited effort has been directed to the analysis of natural convection on non-Newtonian fluids obeying the power-law model of viscosity.4−21 Kim et al.5 demonstrated that the effects of convective transport strengthen with decreasing power-law exponent n based on transient computational simulations of natural convection in square enclosures with differentially heated vertical sidewalls subjected to constant wall temperature (CWT). Similar qualitative behavior has been demonstrated by Lamsaadi et al.6 for natural convection of power-law fluids in shallow rectangular enclosures with differentially heated horizontal walls subjected to constant wall heat flux (CWHF). The same problem has recently been revisited by Alloui et al.,7 who indicated that any nonzero value of Rayleigh number is sufficient for the onset of convection in a shallow enclosure for shear-thickening fluids (i.e., n > 1), whereas convection sets in only once a critical Rayleigh number is surpassed in the case of shear-thinning (n < 1) and Newtonian fluids (n = 1). The onset of convection in strongly temperaturedependent power-law fluids in rectangular enclosures was analyzed in detail by Solomatov and Barr8 in the context of geological applications (i.e., onset of convection in planetary mantles). Albaalbaki and Khayat9 used the Carreau−Bird model © 2013 American Chemical Society

of viscosity (similar to the power-law model but with limiting zero and infinite shear rate viscosities) to analyze threedimensional pattern formation and the onset of convection in generalized Newtonian fluids (GNF), which also demonstrated the effects of rolls, squares, and hexagons on the rate of heat transfer and the influences of Prandtl number on this behavior. Ohta et al.10 analyzed transient natural convection of shearthinning fluids in rectangular enclosures with differentially heated horizontal walls and adiabatic side walls, which demonstrated an increasing trend of mean Nusselt number Nu with increasing extent of shear-thinning. This behavior has been found to be consistent with both experimental and numerical analyses in the same configuration by Inaba et al.11 Similar observations have been made by Khezzar et al.,12 who analyzed the effects of inclination of rectangular enclosures on natural convection of power-law fluids. Natural convection of power-law fluids in rectangular enclosures with differentially heated vertical side walls subjected to CWHF was analyzed by Lamsaadi et al.13,14 in the high Prandtl number limit for both tall13 and shallow enclosures14 where the mean Nusselt number becomes dependent only on the nominal Received: Revised: Accepted: Published: 456

July 22, 2013 November 26, 2013 December 4, 2013 December 4, 2013 dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Rayleigh number Ra and the power-law index n. The same qualitative behavior was reported by Turan et al.,15,16 who also proposed correlations for natural convection of power-law fluids in square enclosures with differentially heated vertical side walls for both CWT15,16 and CWHF16 boundary conditions (bcs). Recently, Turan et al.17 also demonstrated the effects of aspect ratio (i.e., ratio of enclosure height to length) on natural convection of power-law fluids in rectangular enclosures with differentially heated vertical sidewalls subjected to both CWT and CWHF bcs and proposed correlations for the mean Nusselt number Nu. The analysis by Kaddiri et al.18 concentrated on the effects of temperature dependence of viscous properties, whereas computational simulation results were used by Turan et al.19 to propose a correlation for Nu, which was shown to predict the heat transfer behavior accurately for a range of different values of nominal Rayleigh and Prandtl numbers for power-law fluids. Barth and Carey20 used more sophisticated inelastic models, incorporating limiting viscosities at low and high shear rates, to simulate three-dimensional natural convection in cubic enclosures with differentially heated horizontal walls but a linear variation in temperature is used for the side walls to match the experimental conditions of Leung et al.21 Most existing analyses of natural convection of power-law fluids have been carried out for configurations involving two opposite walls subjected to either CWT 5,9−12,15−17 or CWHF6,13,14,16−19 bcs, whereas the other two walls are considered to be adiabatic. Turan et al.16,17,22 demonstrated that the mean Nusselt number Nu for natural convection of power-law and Newtonian fluids in rectangular enclosures with differentially heated vertical sidewalls can be significantly different for CWT and CWHF bcs. In addition, the recent analysis of Turan et al.16 demonstrated that the mean Nusselt number for natural convection of power-law fluids in square enclosures with differentially heated horizontal walls subjected to the CWHF boundary condition (bc) can be significantly different in comparison to the corresponding values obtained for the differentially heated vertical sidewall configuration for same nominal values of Rayleigh and Prandtl numbers. The effects of boundary conditions on natural convection of power-law fluids in rectangular enclosures subjected to differentially heated horizontal walls are yet to be analyzed in detail in the existing literature. Recently, Turan et al.23 analyzed the effects of CWT and CWHF bcs for natural convection of yield-stress fluids obeying the Bingham model in square enclosures with differentially heated horizontal walls with the bottom wall subjected to a higher temperature. However, no such analysis has been carried out for power-law fluids. In the present analysis, the effects of CWT and

CWHF bcs, for differentially heated horizontal walls with the bottom wall subjected to a higher temperature, on natural convection of power-law fluids in square enclosures have been analyzed in detail by carrying out two-dimensional simulations for power-law exponents n ranging from 0.6 to 1.8 for a range of different values of nominal Rayleigh, Ra, and Prandtl, Pr, numbers; Ra = 103 − 105 and Pr = 10 − 105. A value of the power-law exponent n smaller than 0.6 produces a fluid that is extremely shearthinning, similar in some ways to some yield stress fluids. On the other hand, convection becomes extremely weak and, therefore, cannot impart any influence on heat transfer within the enclosure for n ≥ 1.8. Moreover, the range of n considered here remains either comparable to or greater than the range of power-law exponent analyzed in previous studies by the other authors.5,6,13,14,18 The current analysis is primarily concerned with the flow conditions that are sufficiently beyond the critical condition for the onset of fluid motion in the enclosure. However, some limited simulations are undertaken close to this boundary to test the semianalytical results of Alloui et al.7 In this respect, the main objectives of the present analysis are as follows: 1. To demonstrate and explain the differences in the heat transfer behavior of power-law fluids in square enclosures due to CWT and CWHF bcs for the differentially heated horizontal walls configuration. 2. To propose correlations for the mean Nusselt number Nu for natural convection of power-law fluids in the differentially heated horizontal walls configuration. The rest of the paper will be organized as follows. The necessary mathematical background and information related to the numerical implementation will be provided in the next two sections. This will be followed by an overview of existing scaling analysis to elucidate the effects of the power-law exponent and Rayleigh and Prandtl numbers on the mean Nusselt number. Following this, results will be presented and subsequently discussed. Finally, our main findings will be summarized and conclusions will be drawn.

2. MATHEMATICAL BACKGROUND 2.1. Governing Equations and Boundary Conditions. According to the power-law model, the viscous stress tensor τij is given by ⎛ e e ⎞(n − 1)/2 eij τij = K ⎜ kl kl ⎟ ⎝ 2 ⎠

(1)

where eij = (∂ui/∂xj + ∂uj/∂xi) is the strain rate tensor, uj is the jth component of velocity, K is the consistency, and n is the power-law

Figure 1. Schematic diagram of the simulation domain in the (a) CWT configuration and the (b) CWHF configuration. 457

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

index. In the present study, natural convection of power-law fluids in a square enclosure (of dimension L) with differentially heated horizontal walls has been analyzed for different values of n and the nominal values of Rayleigh and Prandtl numbers for both CWT and CWHF bcs. Schematic diagrams of the simulation domains for both the CWT and CWHF configurations are presented in Figure 1. Here, the thermophysical properties such as consistency K, thermal conductivity k, and specific heat capacity cp are taken to be constant and independent of temperature following several previous analyses.1−7,9−17,19−23 For steady flow of incompressible power-law fluids, the conservation equations for nondimesnional mass, momentum, and energy in steady state can be written using tensor notation (i.e., x1 = x is the horizontal direction and x2 = y is the vertical direction) as Nondimensional mass conservation equation

∂ui+ =0 ∂xi+

where g is the acceleration due to gravity and β is the thermal expansion coefficient. Equations 6i and 6ii will be used for the remainder of this paper; they are the same as those used by Ng and Hartnett24 and Lamsaadi et al.6,13,14 Of course, as for all such dimensionless numbers for inelastic non-Newtonian fluids, these definitions are not unique and other, essentially equally valid, definitions can be used. Equations 2−4 are solved in conjunction with following bcs: Velocity boundary conditions u+1 = 0 and u+2 = 0 at x+2 = 0 and x+2 = 1.0 due to no-slip condition and impenetrability on horizontal walls. u+2 = 0 and u+1 = 0 at x+1 = 0 and x+1 = 0 due to no-slip condition and impenetrability of vertical walls. Temperature boundary conditions ∂θ/∂x+1 = 0 at x+1 = 0 and x+1 = 1.0 due to adiabatic vertical walls. θ = 1 (∂θ/∂x+2 = −1) at x+2 = 0 at the hot horizontal wall for CWT (CWHF) bc. θ = 0 (∂θ/∂x+2 = −1) at x+2 = 1 at the cold horizontal wall for CWT (CWHF) bc. According to Buckingham’s Π theorem, the mean Nusselt number Nu = ∫ L0 (Nu/L) dx1 for CWT (CWHF) bc for natural convection of power-law fluids in square enclosures can be expressed as: Nu = f CWT(RaCWT, Pr, n) (Nu = f CWHF(RaCWHF, Pr, n)), where the local Nusselt number Nu is given by

(2)

Nondimensional momentum conservation equations ∂u + uj+ i+ ∂xj

∂τij+ ∂P+ = − + + δi2RaPrθ + + ∂xj ∂xi

(3)

Nondimensional energy conservation equation uj+

∂θ ∂ 2θ + = ∂xj ∂xj+∂xj+

Nu =

(4)

where the spatial coordinates xi, velocity components ui, pressure P, shear stress τij, and temperature T are nondimensionalised in the following manner: x u P ui+ = i P+ = xi+ = i L Uref ρUref 2 τijL (T − Tref ) θ= τij+ = ραUref ΔTref (5)

gβ ΔTL2n + 1 α n(K /ρ) ⎛K⎞ Pr = ⎜ ⎟α n − 2L2 − 2n ⎝ ρ⎠

g β q L2n + 2 (K /ρ) α nk ⎛K⎞ Pr = ⎜ ⎟α n − 2L2 − 2n ⎝ ρ⎠

RaCWHF =

GrCWT =

h = −k

q ∂T 1 × = ∂x 2 (Tx2 = 0 − Tx2 = L) (Tx2 = 0 − Tx2 = L) (8)

where subscript “wf” refers to the condition of the fluid in contact with the wall. 2.2. Numerical Implementation. A finite-volume formulation has been used to solve the coupled conservation equations of mass, momentum, and energy where the diffusive and convective terms are discretised using second-order central difference and second-order up-wind schemes respectively. A commercial FLUENT code, which was used successfully in a number of previous analyses15−17,22,23 involving non-Newtonian fluids has been used here. The well-known SIMPLE (semiimplicit method for pressure-linked equations) algorithm25 has been used to solve the coupled system of mass, momentum, and energy conservation equations. A convergence criterion of 10−7 was used for all the relative (scaled) residuals. There is a possibility of singularity in the power-law model when shear rate approaches zero. Therefore, it is a standard practice to truncate the maximum and minimum permissible values of viscosity μ for power-law fluid simulations to avoid this similarity. These choices should not affect the results as long as the range of μ values obtained from the simulation remains within the predetermined minimum and maximum limits. Here, the minimum and maximum permissible values of μ are taken to be μmin = (10−4 μn=1) and μmax = (104 μn=1), respectively, where μn=1 is the viscosity of the Newtonian fluid for the same nominal values of Rayleigh and Prandtl numbers. For all the simulations considered here, the viscosity μ remained with (10−4 μn=1) and (104 μn=1), and thus the results remain independent of the choices of μmin and μmax.

gβ ΔTL4n − 1 (K /ρ)2 α 2n − 2 (6i)

GrCWHF =

(7)

where the heat transfer coefficient h is defined as

where ρ is the fluid density, Uref = α/L is the reference velocity scale, ΔTref is a reference temperature difference, α = k/(ρcp) is the thermal diffusivity. For the CWT (CWHF) configuration ΔTref can be taken to be equal to (TH − TC)(qL/k), where q is the wall heat flux. The reference temperature Tref in eq 5 is taken to be the temperature at the geometrical center of the domain (i.e., Tref = Tcen) for the CWHF configuration.6,13,14,16−19 In contrast, Tref is taken to be the cold wall temperature TC for the CWT configuration.5,9−12,15−17 In eqs 2−4 Ra = (Gr Pr) and Pr are nominal Rayleigh and Prandtl numbers, respectively, and Gr is the nominal Grashof number. For CWT and CWHF bcs the Rayleigh, Grashof, and Prandtl numbers can be defined as RaCWT =

hL k

g β q L4n (K /ρ)2 α 2n − 2k (6ii) 458

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Table 1. Scaling Estimates of Wall Heat Flux q, Characteristic Velocity ϑ, Thermal Boundary Layer Thicknesses δth, Nusselt Number Nu, Effective Viscosity μeff, Effective Grashof Number Greff, and Effective Rayleigh Number Raeff. quantity

CWT boundary condition

CWHF boundary condition

wall heat flux q

ΔT q∼k ∼ hΔT δth

characteristic velocity scale ϑ

ϑ∼

thermal boundary layer thickness δth

⎡ ⎛ KL(gβ ΔTL)(n /2) − 1 ⎞1/(n + 1)⎤ 1 ⎢ ⎥ ⎜ ⎟⎟ δth ∼ min⎢L , ⎥ ρ f1 (RaCWT , Pr, n) ⎜⎝ ⎠ ⎣ ⎦

ΔT q∼k ∼ hΔT δth

gβ ΔTL

ϑ∼

⎡ ⎤ 1 L ⎥ δth ∼ min⎢L , (n + 1)/((n /2) + 2) 1 ( /2) /2 1/(( /2) 2) − n − n n + ⎢⎣ f ⎥⎦ (RaCWHF ) Pr 2 ⎡ ⎤ 1 L ⎥ ∼ min⎢L , ⎢⎣ f3 (RaCWHF , Pr , n) (RaCWHF1 − (n /2)Pr −n /2)1/((n /2) + 2) ⎥⎦

⎡ ⎛ ⎞⎤ 1 L ⎜ ⎟⎥ ∼ min⎢L , ⎢⎣ f1 (RaCWT , Pr , n) ⎜⎝ (RaCWT 2 − nPr −n)1/(2(n + 1)) ⎟⎠⎥⎦ Nusselt number Nu

Nu ∼

hL L L ∼ ∼ f1 k δth δ

when Nu > 1,

Nu ∼

Nu ∼ (RaCWT 2 − nPr −n)1/(2(n + 1))f1 (RaCWT , Pr , n) effective viscosity μeff∼K(ϑ/δ)n−1

g β q δth L k

hL L L ∼ ∼ f2 k δth δ

when Nu > 1,

Nu ∼ (RaCWHF1 − (n /2)Pr −n /2)1/((n /2) + 2) × f3 (RaCWHF , Pr , n)

⎛ K ⎞2/(n + 1)⎛ (gβ ΔTL)3(n − 1)/(2(n + 1)) ⎞ ⎜ ⎟ μeff ∼ ρ⎜ ⎟ ⎜ ⎟ ⎝ ρ⎠ L(n − 1)/(n + 1) ⎝ ⎠

⎛ K ⎞5/(n + 4)⎛ gβq ⎞(3n − 3)/(n + 4) (2n − 2)/(n + 4) ⎟ ⎜ μeff ∼ ρ⎜ ⎟ L ⎝ k ⎠ ⎝ ρ⎠ × [f2 (RaCWHF , Pr , n)](2n

effective Grashof number Greff

GrCWT,eff =

effective Rayleigh number Raeff

RaCWT,eff =

ρ2 gβ ΔTL3 μeff

2

∼ RaCWT(4 − 2n)/(n + 1)Pr (−2n)/(n + 1)

GrCWHF,eff =

ρ2 gβ ΔTL3 ⎛ μeff c p ⎞ ⎜⎜ ⎟⎟ ∼ RaCWT(5 − n)/(2n + 2)Pr (1 − n)/(2n + 2) μeff 2 ⎝ k ⎠

RaCWHF,eff =

ρ2 gβqL4 μeff 2 k

2

)( )

+ 3n − 5 / n + 4

∼ Ra CWHF(10 − 5n)/(n + 4)Pr (−5n)/(n + 4)

ρ2 c pgβqL4 ⎛ μeff c p ⎞ ⎟ ∼ Ra CWHF(7 − 2n)/(n + 4)Pr (2 − 2n)/(n + 4) ⎜ μeff k 2 ⎜⎝ k ⎟⎠

ΔT is the characteristic temperature difference between hot and cold walls and f1 and f 2 are the functions of Rayleigh number, Prandtl number, and power-law index, which relate the thermal boundary-layer thickness δth to the hydrodynamic boundary-layer thickness δ as δ/δth ∼ f1(RaCWT, Pr, n) (δ/δth ∼ f 2(RaCWHF, Pr, n)) + 2) . Table 1 for CWT (CWHF) bc with f 3 being f 3 = f (n+1)/(n/2 2 shows that the mean Nusselt number Nu attains a value equal to unity when δth ∼ L. These scaling predictions provide useful insight into the anticipated behavior of Nu in response to variations of Ra, Pr, and n. Moreover, the Nusselt number scaling presented in Table 1 suggests that Nu is expected to decrease with increasing n for a given value of Ra, whereas Nu increases with increasing Ra for a given value of n for n < 2. Table 1 indicates that the effective values of Ra become increasingly larger than their nominal values for decreasing values of n for shear-thinning fluids. This suggests that the magnitudes of effective Rayleigh number may attain large values for small values of n so that a steady two-dimensional laminar solution may not exist, whereas a steady laminar solution can be obtained for the same set of nominal values of Rayleigh and Prandtl numbers for a higher value of n. Thus, a threshold value Rathr can be expected for the effective Rayleigh number such that a steady two-dimensional solution does not exist when RaCWT,eff > RaCWT,thr and RaCWHF,eff > RaCWHF,thr. A number of simulations have been carried out for different values ofRaCWT (RaCWHF), Pr, and n, and it has been found that a converged steady solution cannot be obtained when RaCWT,eff > 1 × 106 and RaCWHF,eff ∼ Ra(7−2n)/(n+4)Pr(2−2n)/(n+4) > 2 × 106 (unsteady numerical solutions can be obtained however), which essentially suggests that steady two-dimensional solutions do not exist for the following condition:

2.3. Computational Grid and Bench-Marking. The simulations have been carried out using three different nonuniform meshes, M1 (100 × 100), M2 (160 × 160), and M3 (200 × 200). Detailed information regarding these meshes has been provided in ref 19 and is thus is not repeated here. The numerical uncertainty levels for the mean Nusselt number Nu and maximum nondimensional horizontal velocity component Umax on the vertical midplane x1/L = 0.5 remain smaller than 1% for the three different meshes for n = 0.6−1.8 (see Table 2 in ref 19). Based on this finding, the simulations have been carried out using mesh M2 for the sake of high accuracy and computational efficiency. It has been demonstrated earlier that the present numerical scheme has been validated with respect to the results obtained for natural convection of Newtonian fluids in square enclosures with differentially heated vertical constant temperature walls,1 and an excellent agreement was achieved (see Table 3 of ref 26). Moreover, the variation of Nu/Nun = 1 with n in square enclosures with differentially heated vertical walls subjected to CWT bc has also been compared to the corresponding results presented by Kim et al.,5 and an excellent agreement was obtained. Interested readers are referred to Figure 2 of ref 15 for further information.

3. SCALING ANALYSIS AND REGIMES OF THERMAL TRANSPORT Detailed scaling analyses are performed in refs 15−17 and 19 to elucidate the anticipated effects of Rayleigh number, Prandtl number, and power-law index on the Nusselt number for natural convection of power-law fluids in square enclosures and, thus, are not repeated here for the sake of conciseness. The scaling estimates of wall heat flux q, characteristic velocity scale ϑ, thermal boundary layer thickness δth, Nusselt number Nu, effective viscosity μeff, effective Grashof number Greff, and effective Rayleigh number Raeff, according to refs 15−17 and 19, are summarized in Table 1, where

RaCWT > [1 × 106Pr (n − 1)/(2n + 2)](2n + 2)/(5 − n) 459

(9i)

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Table 2. Value of Nominal Rayleigh Number Ra, above Which Nu Deviates from Unity in the Third Decimal Place (i.e., Nu ≥ 1.001) for Different Values of n and Pr CWHF

CWT

n

Pr = 10

Pr = 102

Pr = 103

Pr = 104

Pr = 10

Pr = 102

Pr = 103

Pr =104

0.6 0.7 0.8 0.9 1.0 1.2 1.4 1.6 1.8

600 800 1100 1400 1700 1300 1000 800 700

600 800 1100 1400 1700 1300 1000 800 700

600 800 1100 1400 1700 1300 1000 800 700

600 800 1100 1400 1700 1300 1000 800 700

700 1000 1500 2000 2500 3300 4300 5600 7300

700 1000 1500 2000 2500 3300 4300 5600 7300

700 1000 1500 2000 2500 3300 4300 5600 7300

700 1000 1500 2000 2500 3300 4300 5600 7300

Table 3. Value of Critical Rayleigh Number Racrit, above Which Nu Deviates from Unity in the Third Decimal Place (i.e., Nu ≥ 1.001) for Different Values of n and Pr CWHF

CWT

n

Pr = 10

Pr = 102

Pr = 103

Pr = 104

Pr = 10

Pr = 102

Pr = 103

Pr = 104

0.6 0.7 0.8 0.9 1.0 1.2 1.4 1.6 1.8

4 751 3 861 3 198 2 396 1 700 476 153 56 24

7 091 5 180 3 875 2 633 1 700 398 108 34 13

10 583 6 950 4 694 2 892 1 700 334 77 21 7

15 795 9 325 5 687 3 177 1 700 279 55 12 3

10 890 7 627 5 768 3 872 2 500 985 438 216 116

14 521 9 345 6 555 4 114 2 500 887 362 166 84

19 365 11 450 7 449 4 371 2 500 799 299 127 60

25 823 14 030 8 466 4 644 2 500 719 246 98 43

RaCWHF > [2 × 106Pr (2n − 2)/(n + 4)] (n + 4)/(7 − 2n)

fluid motion in the enclosure. However, a criterion based on velocity or stream function will affect the mean Nusselt number differently for different values of n, so it is not clear whether the degree of arbitrariness can be eliminated if a condition is devised based on the velocity or stream function (e.g., the maximum value of the stream function may not be associated with the near-wall region, so the maximum value of stream function may not offer much information regarding wall heat flux). Moreover, the findings based on a criterion based on the velocity or stream function will be qualitatively similar to those indicated in Tables 2 and 3. In the present analysis, we are interested in identifying the regime where the thermal transport takes place primarily due to conduction by using the criterion of Nu < 1.001, and this serves the main purpose of the current analysis. Moreover, the critical Rayleigh number RaCWT,crit (RaCWHF,crit) for the effective Rayleigh numbers RaCWT,eff (RaCWHF,eff) decreases with increasing n for both CWT and CWHF bcs, which is qualitatively consistent with the semianalytical results of Alloui et al.7 The criterion above which convection starts to play a key role in the heat transfer rate (which is reflected in Nu > 1) can be expressed as follows: for n ≤ 1

(9ii)

Moreover, it has been found that the mean Nusselt number Nu remains unity for RaCWT,eff ∼ RaCWT(5−n)/(2n+2)Pr(1−n)/(2n+2) < Racrit and RaCWHF,eff ∼ RaCWHF(7−2n)/(n+4)Pr(2−2n)/(n+4) < Racrit, where Racrit is the critical Rayleigh number below which the thermal transport is principally conduction-driven. A series of simulations revealed that Nu values deviate from unity in the third decimal place (i.e., Nu ≥ 1.001) for RaCWT > 2500n2.5 (RaCWT > 2500e1.35(n−1)) and RaCWHF > 1700n2 (RaCWHF > 1700n−3/2) for n ≤ 1 (n > 1) for the CWT and CWHF bcs irrespective of the value of Prandtl number as shown in Table 2. The onset of convection for RaCWT > 2500 for n = 1.0 is entirely consistent with Ternik et al.27 for the same configuration. The precise numerical value of this onset Rayleigh number is likely to be dependent on the exact numerical scheme used. For the current analysis we are interested in an accurate order-ofmagnitude estimation of Racrit (however, the critical Rayleigh number reported by Ternik et al.27 is only 3.5% greater than 2500 for the CWT bc) that is sufficient to identify the regime boundary between conduction and convection driven transports for the current analysis. The critical Rayleigh number obtained from simulation data is summarized in Table 3, which shows that the nominal Rayleigh number above which Nu ≥ 1.001 increases (decreases) with n in the case of CWT (CWHF) bc, and further analyses are needed to explain this important difference in behavior. A detailed analysis of the precise onset of fluid motion in this configuration is beyond the scope of the current analysis and the values of Ra (Racrit) in Table 2 (Table 3) are only indicative in nature because Nu ≥ 1.001 is not a sufficiently accurate proof of the onset of fluid motion. Perhaps a criterion based on the velocity or stream function may be more appropriate for identifying the onset of

RaCWT,eff ∼ Ra(5 − n)/(2n + 2)Pr (1 − n)/(2n + 2) > RaCWT,crit = (2500n2.5)(5 − n)/(2n + 2) Pr(1 − n)/(2n + 2) (10i)

for n ≤ 1 RaCWT,eff ∼ Ra(5 − n)/(2n + 2)Pr (1 − n)/(2n + 2) > RaCWT,crit = (2500e1.35(n − 1))(5 − n)/(2n + 2)Pr(1 − n)/(2n + 2) (10ii) 460

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 2. (a) Different regimes of convection for n = 0.6, (b) temporal evolution of Nu with dimensionless time τ = αt/L2 at Pr = 50, n = 0.6 for cases A (RaCWT (RaCWHF) = 500), B (RaCWT (RaCWHF) = 1.5 × 104), C (RaCWT (RaCWHF) = 5 × 104), and D (RaCWT (RaCWHF) = 1 × 105). The boundaries of the regimes for the CWT and CWHF configuration are shown by red and blue lines, respectively (the dashed line highlights the side wall case).

for n ≤ 1

for n ≤ 1 RaCWHF,eff ∼ Ra(7 − 2n)/(n + 4)Pr (2 − 2n)/(n + 4) 2 (7 − 2n)/(n + 4)

> RaCWHF,crit = (1700n )

[1 × 106 Pr (n − 1)/(2n + 2)](2n + 2)/(5 − n) > RaCWT > 2500n2.5 (2 − 2n)/(n + 4)

Pr

(11i)

(10iii)

for n > 1

for n > 1 RaCWHF,eff ∼ Ra(7 − 2n)/(n + 4)Pr (2 − 2n)/(n + 4)

[1 × 106 Pr (n − 1)/(2n + 2)](2n + 2)/(5 − n) > RaCWT > 2500e1.35(n − 1) (11ii)

> RaCWHF,crit = (1700n−3/2)(7 − 2n)/(n + 4)Pr(2 − 2n)/(n + 4)

for n ≤ 1

(10iv)

The above discussion suggests that the steady two-dimensional laminar convection can be realized in the current configuration when the following condition is met:

[2 × 106 Pr (2n − 2)/(n + 4)] (n + 4)/(7 − 2n) > RaCWHF > 1700n2 (11iii) 461

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 3. Variations of nondimensional temperature (a) θCWT = (T − Tcen)/(TH − TC) and (b) θCWHF = ((T − Tcen)k)/qL along the vertical midplane (i.e., x1/L = 0.5) for RaCWT(RaCWHF) = 1 × 104, 5 × 104, and 1 × 105 at Pr = 104.

given by [1 × 106 Pr(n−1)/(2n+2)](2n+2)/(5−n) > RaCWT and RaCWT,eff > RaCWT,crit ([2 × 106 Pr(2n−2)/(n+4)](n+4)/(7−2n) > RaCWHF and RaCWHF,eff > RaCWHF,crit) in Figure 2a is termed as the “steady two-dimensional laminar convection regime”. As steady twodimensional laminar solutions do not exist for RaCWT > [1 × 106 Pr(n−1)/(2n+2)](2n+2)/(5−n) (RaCWHF > [2 × 106 Pr(2n−2)/(n+4)](n+4)/(7−2n)), the corresponding regime is referred to as the “unsteady convection regime”. The validity of the above regime diagram can be substantiated from cases A, B, C, and D, shown in Figure 2. For cases A and B, Nu is a constant (i.e., case A: Nu = 1.0) as predicted by the regime

for n > 1 [2 × 106 Pr(2n − 2)/(n + 4)](n + 4)/(7 − 2n) > RaCWHF > 1700n−3/2 (11iv)

Equations 10i−11iv are partially based on scaling arguments and, thus, are valid strictly only in an order-of-magnitude of sense. The conditions given by eqs 10i−11iv are highlighted on a regime diagram shown in Figure 2a. For RaCWT,eff < RaCWT,crit (RaCWHF,eff < RaCWHF,crit), the heat transfer takes place principally due to thermal conduction, and this regime has been termed as the “conduction-dominated regime” in Figure 2a. The region 462

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 4. Variations of nondimensional horizontal velocity component U for (a) CWT and (b) CWHF boundary conditions along the vertical midplane (i.e., x1/L = 0.5) for RaCWT = RaCWHF = 1 × 104, 5 × 104, and 1 × 105 at Pr = 104.

diagram. For case C, a transient simulation has been carried out for the CWT configuration and the temporal evolution of Nu shows a complex temporal oscillation, but Nu for the CWHF configuration exhibits a constant value (see Figure 2b). The transient simulation for case D yielded a complex temporal oscillation of Nu for both CWT and CWHF bcs. The regime diagram boundaries are based on scaling arguments, so these boundaries should not be treated strictly but need to be considered only in an order-of-magnitude sense. The regimes of convection for power-law fluids in the case of differentially heated sidewalls configuration are also shown in

Figure 2a (taken from refs 15 and 16). Interested readers are referred to refs 15 and 16, for the criteria demarcating the regimes in the differentially heated sidewall (SW) configuration. Figure 2a indicates that the conduction-dominated regime can be obtained for a wider parameter range for the differentially heated horizontal wall (HW) configuration than in the SW configuration. The fluid motion initiates as soon as a temperature difference is maintained between the walls in the SW configuration. In contrast, convection starts only once a critical Rayleigh number is surpassed in the HW configuration for n ≤ 1 .6−9 Thus, convection within the enclosure starts to yield Nu > 1 for a higher 463

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 5. Contours of (a) nondimensional temperature θCWT, (b) nondimensional streamlines (Ψ = ψ/α), and (c) normalized viscosity (μ/μn=1) for n = 0.8, and 1.8 at Pr = 104.

4. RESULTS AND DISCUSSION Distributions of Nondimensional Temperature and Velocity Fields. The variations of nondimensional temperature along the vertical midplane for RaCWT = (RaCWHF =) 1 × 104, 5 × 104, and 1 × 105 at Pr = 104 are shown in Figure 3a and b, respectively, for n ranging from 0.6 to 1.8. Only the results where steady state two-dimensional solutions exist are shown in Figure 3. The temperature and velocity variations are identical for the range of Prandtl number Pr analyzed here (10−105) and are, therefore, not explicitly shown. Figure 3 shows that the variation of nondimensional temperature in the vertical direction x2 becomes increasingly nonlinear with increasing (decreasing) values of Rayleigh number (power-law exponent n) for a given set of values of power-law exponent (RaCWT or RaCWHF) and Pr. The linear temperature variation in the vertical direction is indicative of a pure-conduction solution and the extent of nonlinearity in the temperature distribution in the vertical direction increases with increasing convection strength within the enclosure. Figure 3a and b show that the maximum magnitude of θCWT remains 0.5 because the wall temperatures are kept constant, whereas the maximum magnitude of θCWHF decreases with decreasing (increasing) n (RaCWHF) in the CWHF configuration.

value of Rayleigh number in the HW configuration than in the SW configuration. However, steady two-dimensional convection in the HW configuration is inherently more unstable than in the SW configuration because of the instability induced by the presence of heavier cold fluid directly on top of lighter hot fluid. As a result, steady two-dimensional convection can be obtained over a significantly wider range of parameters in the SW configuration than in the current HW configuration. The present analysis focuses on steady state two-dimensional convection of power-law fluids; thus, the simulations within the steady two-dimensional laminar convection regime have been carried out in the steady state mode. The steady state solutions for a Newtonian fluid (i.e., n = 1) for a given set of values of Rayleigh and Prandtl numbers have been used as the initial conditions for power-law fluids. Thus, the simulations involving power-law fluids were initiated from an already bifurcated solution for the power-law fluids, avoiding the issue of tracking the bifurcated solution branch from below. Moreover, the inherent round-off and truncation errors in the numerical scheme provided sufficient perturbation for initiating convection from the pure conduction solution in the current configuration for Newtonian fluids. 464

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 6. Variation of the mean Nusselt number Nu with nominal Rayleigh number RaCWT (RaCWHF) for different values of power law index n for (a) Pr = 100, (b) Pr = 103, (c) Pr = 104, and (d) Pr = 105 (red line highlights CWT case).

The nondimensional temperatures θCWT and θCWHF can be estimated in the following manner using scaling estimates provided in Table 1 θCWT ∼ O(1)

of U = u1L/α along the vertical direction x2 for both CWT and CWHF bcs, which are presented in Figure 4a and b, respectively. It is evident from Figure 4a and b that the magnitude of U increases with increasing (decreasing) Rayleigh number (powerlaw exponent n) for a given set of values of power-law exponent (RaCWT or RaCWHF) and Pr. For a square enclosure, the orders of magnitude of horizontal and vertical velocity components remain of the same order as governed by continuity (i.e., u1/L ∼ u2/L); thus, the characteristic velocity scale for this configuration can be estimated as

⎛ qδ k ⎞ ⎛δ ⎞ θCWHF ∼ O⎜ th ⎟ ∼ O⎜ th ⎟ ⎝ L⎠ ⎝ qLk ⎠

∼ f2−(n + 1)/((n /2) + 2) (RaCWHF1 − (n /2)Pr −n /2)−1/((n /2) + 2) (12)

Equation 12 suggests that θCWHF is expected to decrease with increasing (decreasing) RaCWHF (power-law exponent) for a given set of values of n (RaCWHF) as observed in Figure 3a and b. This indicates that for large (small) values of n (RaCWHF), the distribution of θCWHF for the CWHF configuration will approach the corresponding temperature distribution for the CWT bc for identical numerical values of nominal Ra and Pr. This convergence can be confirmed from Figure 3a and b, which show that the distributions of θCWHF and θCWHF approach each other for Rayleigh number equal to 1 × 104 and for n = 1.8. The smaller temperature difference between the horizontal walls in the CWHF configuration than that in the corresponding CWT configuration gives rise to weaker buoyancy-driven flow in the CWHF configuration than in the corresponding CWT configuration. This can be further substantiated by the variation

U ∼ u1L /α ∼

RaCWTPr

(for CWT)

(13i)

U ∼ u1L /α n

1

∼ RaCWHF( 2 + 2 ) (for CWHF)

( 2n + 2)Pr( 2n + 1) ( 2n + 2)f −( 2n + 12 ) ( 2n + 1) 2

(13ii)

According to eqs 13i and 13ii, U is expected to increase with increasing Rayleigh number in both CWT and CWHF bcs as observed from Figure 4a and b. The effective viscosity μeff increases with increasing n for a given set of nominal values of Ra and Pr according to Table 1. Moreover, (n/2 + 1/2)/(2 + n/2) remains smaller than 0.5 (i.e., (n/2 + 1/2)/(2 + n/2) < 0.5) 465

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

for 0.6 ≤ n ≤ 1.8; thus, U is expected to assume greater values in the CWT configuration than in the CWHF configuration according to eqs 13i−13ii. This behavior can indeed be verified from Figure 4a and b. Moreover, Table 1 indicates that the effective Grashof number increases with decreasing n for a given set of values of nominal Rayleigh and Prandtl numbers (i.e., RaCWT (RaCWHF) and Pr). As the effective Grashof number indicates the ratio of the strengths of buoyancy to viscous forces, an increasing trend of effective Grashof number with decreasing n indicates strengthening of buoyancy forces in comparison to viscous forces. The strengthening (weakening) of convection (viscous resistance) with decreasing n and increasing Rayleigh number is further supported by the contours of θCWT, nondimensional stream function Ψ = ψ/α (i.e., ψ is the dimensional stream function), normalized viscosity μ/μn=1 in the CWT configuration shown in Figure 5 for RaCWT = 103 and 5 × 104 and n = 0.8 and 1.8 at Pr = 104. Similar qualitative distributions of θCWHF, Ψ, and μ/μn=1 have been observed in the CWHF configuration and, thus, are not shown here. The distribution of μ/μn=1 in Figure 5 is qualitatively consistent with previous findings by Albaalbaki and Khayat.9 Figure 5 reveals that the viscosity in shear-thinning (shear-thickening) fluids assumes smaller (greater) values than in Newtonian fluids for the major part of the flow domain; thus, the effects of viscous resistance (advective transport) strengthen (weaken) with increasing n. Thus, the isotherms become increasingly curved and the magnitude of Ψ increases with decreasing (increasing) n (nominal Rayleigh number) for a given value of nominal Rayleigh number (n). For example, at RaCWT = 103 for n = 0.8 and 1.8, the isotherms remain parallel to the horizontal boundaries, indicating a stable thermal stratification characteristic of pure conduction. Moreover, Figure 5 demonstrates that the isotherms remain almost parallel to horizontal boundaries in n = 1.8 for RaCWT ≤ 104, where Ψ is negligible. Thus, the thermal transport becomes increasingly conduction-driven for decreasing (increasing) values of nominal Rayleigh number (power-law exponent). Behavior of the mean Nusselt number and differences between CWT and CWHF boundary conditions. Scaling relations shown in Table 1 indicate that a decrease in δth leads to an increase in Nu in both CWT and CWHF configurations, which is consistent with the expected increase in the rate of heat transfer with the strengthening of thermal advection. The variation of Nu with RaCWT (RaCWHF) for different values of n for Pr = 102, 103, 104, and 105 are shown in Figure 6, which indicates that Nu increases with decreasing (increasing) n (nominal Rayleigh number) for both CWT and CWHF bcs. This behavior is consistent with the scaling estimates of Nu for 0.6 ≤ n ≤ 1.8. Figure 6 further demonstrates that Nu for Ra > 5 × 103 assumes higher values in the CWT configuration than in the CWHF configuration for n ≤ 1 for all values Prandtl numbers considered here. However, for shear-thickening fluids Nu for the CWT bc are found to be smaller than in the CWHF configuration for some combination of RaCWHF and Pr, whereas for other combinations of Rayleigh and Prandtl numbers, Nu in the CWT bc is found to be greater than the corresponding value in the CWHF configuration. This behavior can be explained as follows. The exponent of Ra CWHF in Ra CWHF,eff ∼ RaCWHF(7−2n)/(n+4)Pr(2−2n)/(n+4) assumes greater positive values than the exponent of RaCWT in RaCWT,eff ∼ RaCWT(5−n)/(2n+2)Pr(1−n)/(2n+2), whereas the exponent of Pr is more negative for the CWHF case than in the CWT configuration for shear-thickening fluids.

The difference between exponents of Pr increases with increasing n for shear-thickening fluids. The smaller positive exponent of Rayleigh number in RaCWT,eff than in RaCWHF,eff overcomes the effects of smaller negative exponent of Pr in the effective Rayleigh number expression for the CWT configuration than in the CWHF configuration, to give rise to a smaller effective Rayleigh number in the CWT configuration than in the CWHF bc for small values of nominal Rayleigh number in shear-thickening fluids with 1 ≤ n ≤ 1.4 for a given set of values of Pr and n. This difference suggests that the convective transport is stronger in the CWHF bc than in the CWT bc for small values of nominal

Figure 7. Variation of the mean Nusselt number Nu with Prandtl number Pr for (a) n = 0.6, (b) n = 1.0, and (c) n = 1.8 for both CWT (red) and CWHF boundary conditions. 466

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 8. Variations of Nu with the power-law index, n ≤ 1 for CWT (▽) and CWHF (○) for different values of Pr and RaCWT (RaCWHF) along with the predictions of eqs 15i−15ix for CWT () and CWHF (- - -) boundary conditions.

Rayleigh number in shear-thickening fluids with n ≤ 1.4, which is reflected in the smaller value of Nu in the CWT bc than in the CWHF configuration. However, the smaller negative exponent of Pr in RaCWT,eff than in RaCWHF,eff overcomes the effects of smaller positive exponent of Rayleigh number in RaCWT,eff than in RaCWHF,eff for high values of nominal Rayleigh number for 1.0 ≤ n ≤ 1.4. This leads to a greater value of effective Rayleigh number in the CWT configuration than in the CWHF bc for relatively large values of nominal Rayleigh number for a given value of Pr and n in the case of 1 ≤ n ≤ 1.4. In turn giving rise to higher values of Nu in the CWT bc than in the CWHF configuration for relatively large values of Ra for a given set of values of Pr and n. For n > 1.4 the effects of the smaller positive exponent of Rayleigh number in RaCWT,eff than in RaCWT,eff overcomes the effects of smaller negative exponent of Pr in RaCWT,eff than in RaCWHF,eff to give rise to a smaller effective Rayleigh number in the CWT configuration than in the CWHF bc for a given set of values of nominal Ra and Pr. Thus, the effects of convection are found to be weaker in the CWT configuration

than in the CWHF bc. Furthermore, Nu in the CWHF configuration assumes greater values than in the case of the CWT bc for a given set of nominal values of Rayleigh and Prandtl numbers for shear-thickening fluids with n > 1.4. According to Table 1, RaCWT,eff > RaCWHF,eff for a given set of Rayleigh and Prandtl numbers for the shear-thinning fluids, which leads to higher values of Nu in the CWT configuration than in the case of CWHF bc due to stronger convective transport in the case of CWT bc. The exponent of Rayleigh number in the mean Nusselt number scaling (see Table 1) is smaller in the CWHF bc than in the CWT configuration (0.2 in CWHF and 0.25 in CWT) for Newtonian fluids, which gives rise to a greater value of Nu for the CWT bc than in the CWHF configuration for a given set of values of Rayleigh and Prandtl numbers within the parameter range considered here. The differences between the mean Nusselt number Nu for the CWT and CWHF configurations for shearthinning and Newtonian fluids in this configuration is found to be in qualitative agreement with previous findings in the context of differentially heated vertical side walls.15,16 467

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 9. Variations of Nu with the power-law index, n > 1 for CWT (▽) and CWHF (○) for different values of Pr and RaCWT (RaCWHF) along with the predictions of eqs 15i−15ix for CWT () and CWHF (- - -) boundary conditions.

Effects of Nominal Prandtl Number Pr. The variations of Nu with Pr for different values of RaCWT (RaCWHF) and n are shown in Figure 7 for the conditions under which steady twodimensional laminar solutions exist. Figure 7 shows that Pr has a marginal influence on Nu for a given set of values of nominal Rayleigh number and power-law exponent. That Nu is independent of Pr for Pr ≫ 1 is broadly consistent with earlier findings for differentially heated vertical9,15,16,28,29 and horizontal19,27 walls. Recent findings of Ternik and Rudolph28,29 and Ternik et al.27 also confirms that Nu remains independent of Pr for both differentially heated vertical and horizontal wall configurations for nanofluids even for a smaller range of Prandtl number (i.e., 1 < Pr < 10) than the values considered here. For Pr ≫ 1 the hydrodynamic boundary-layer thickness is thicker than the thermal boundary layer thickness (i.e., δ ≫ δth); thus, a considerable amount of unheated fluid within the viscous boundary layer remains unaffected by the presence of the thermal boundary layer.30 The order-of-magnitude analysis

by Bejan30 demonstrates that the balance between viscous and buoyancy forces governs the transport within δth for Pr ≫ 1, whereas the relative balances of inertial and viscous forces govern the transport behavior outside δth within the hydrodynamic boundary layer. Thus, a change in Prandtl number acts to modify the relative balance between viscous and inertial forces in natural convection, whereas the thermal transport within δth is marginally affected for Pr ≫ 1 fluids. As a result, Nu remains insensitive to the variation of Prandtl number for both Newtonian and power-law non-Newtonian fluids for Pr ≫ 1. Correlations for Mean Nusselt Number for Both CWT and CWHF Boundary Conditions. As shown in Table 1, the mean Nusselt number can be taken to scale with Nu ∼ (Ra C W T 2 − n Pr − n ) 1 / ( 2 ( n+ 1 ) ) f 1 (Ra C W T , Pr, n) and Nu ∼ (RaCWHF1−(n/2)Pr−n/2) 1/((n/2)+2))f 3(RaCWHF, Pr, n) and recently Turan et al.23 proposed correlations for Nu for Newtonian fluids in the following manner for 105 ≥ RaCWT ≥ 5 × 10 3 (10 5 ≥ RaCWHF ≥ 5 × 10 3) 468

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 10. Variations of Nu with the power-law index n for both horizontal wall and side wall case for different values of RaCWT (RaCWHF) at Pr = 100.

⎛ Pr ⎞n ⎟ Nu = aRaCWT m⎜ ⎝ 1 + Pr ⎠

and

⎛ Pr ⎞n ⎟ Nu = aRaCWHFm⎜ ⎝ 1 + Pr ⎠ (14i)

Nu = Nu 2

0.014 Nu 2 = 0.289RaCWHF

The values of the coefficients a, m, and n were determined by an iterative minimization function of a commercial software package which provides the following values: a = 0.178 a = 0.289

m = 0.269 m = 0.214

n = 0.02 n = 0.017

(for CWT) (for CWHF)

when

Nu1 = 0.178RaCWT 0.019

(14ii)

Nu = 1

(14iii)

when

Nu1 ≤ 1

when

Nu 2 ≤ 1

(15ii)

where bCWT and bCWHF are the correlation parameters which can be expressed as

bCWT = c1RaCWTc2Pr c3

bCWHF = c1RaCWHFc2Pr c3

(15iii)

where c1, c2, and c3 are given by for the CWT case

c1 = 0.994, c 2 = 0.091, and c3 = 0.051 3

5 × 10 ≤ RaCWT ≤ 1 × 10

Pr 0.21 (1 + Pr )0.02

5

and

5 × 103 ≤ RaCWT ≤ 1 × 104

and

1 × 104 < RaCWT ≤ 1 × 105 469

and

(15iv)

for

n>1

c1 = 0.606, c 2 = 0.104, and c3 = 0.039 (15i)

for

n≤1

c1 = 0.021, c 2 = 0.463, and c3 = 0.049

2 − n ⎞1/(2(n + 1)) ⎛ Ra CWT ⎟ ebCWT(n − 1) > 1 ×⎜ n ⎝ Pr ⎠

Nu = 1

Pr 0.217 (1 + Pr )0.017

⎛ Ra 1 − n /2 ⎞1/(n /2 + 1) ⎟ × ⎜⎜ CWHF ebCWHF(n − 1) > 1 n /2 ⎟ Pr ⎝ ⎠

Equations 14i−14iii satisfactorily captures the variation of Nu in the range given by 5 × 103 ≤ RaCWT ≤ 105, 5 × 103 ≤ RaCWHF ≤ 105, and 10 ≤ Pr ≤ 105 for Newtonian fluids. Thus, the correlation for Nu for power-law fluids should be proposed in such a manner that limn→1 Nu = 0.178Ra CWT0.269[Pr/(1+Pr)]0.02 and limn→1 Nu = 0.289RaCWHF0.214 [ Pr/(1+Pr)]0.017 are satisfied in the CWT and CWHF configurations, respectively. On the basis of the aforementioned limiting conditions, a correlation for Nu is proposed here for one roll flow pattern in a square enclosure in the following manner: Nu = Nu1

when

(15v)

for

n>1

(15vi)

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 11. The variation of nondimensional temperature θ along the vertical midplane x1/L = 0.5 (horizontal midplane x2/L = 0.5) for different n values (n ≤ 1) in the differentially heated horizontal wall (vertical sidewall) configuration for both CWHF (first column) and CWT (second column) cases at Pr = 100. (a) RaCWT (RaCWHF) = 5 × 103 and (b) RaCWT (RaCWHF) = 5 × 104.

shown in Figures 8 and 9 for n < 1 and n > 1 fluids, respectively. The predictions for RaCWT = 103 and RaCWHF = 103 are not shown here because Nu remains equal to unity for all values of n and Pr considered here. Figures 8 and 9 show that eqs 15i−15ix satisfactorily predict the qualitative and quantitative behavior of Nu in both n < 1 and n > 1 fluids. Although eqs 15i−15ix satisfactorily capture the variation of Nu with the power-law exponent and Rayleigh and Prandtl numbers, it overpredicts the value of Nu for n = 0.6 in the CWT configuration. Despite this overprediction, the observations from Figures 8 and 9 indicate that eqs 15i−15ix can be used for the accurate estimation of Nu in twodimensional laminar steady-state natural convection of power-law fluids in the following parameter range: 5 × 103 ≤ RaCWT ≤ 105 (5 × 103 ≤ RaCWHF ≤ 105), 10 ≤ Pr ≤ 105, and 0.6 ≤ n ≤ 1.8. Differences between Differentially Heated Horizontal Wall and Vertical Sidewall Configurations. The variations

for the CWHF case c1 = 0.345, c 2 = 0.129, and c3 = 0.103 5 × 103 ≤ RaCWHF ≤ 1 × 105

and

n≤1

c1 = 0.014, c 2 = 0.458, and c3 = 0.097 5 × 103 ≤ RaCWHF ≤ 1 × 104

and

1 × 10 < RaCWHF ≤ 1 × 10

5

and

(15vii)

for

n>1

c1 = 0.259, c 2 = 0.153, and c3 = 0.073 4

for

(15viii)

for

n>1

(15ix)

As intended, the expression of Nu becomes exactly equal to an existing correlation (eqs 14i−14iii) for Newtonian fluids when n is taken to be unity. The performance of this correlation for the aforementioned range of Rayleigh and Prandtl numbers is 470

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

Figure 12. The variation of nondimensional temperature θ along the vertical midplane x1/L = 0.5 (horizontal midplane x2/L = 0.5) for different n values (n ≥ 1) in the differentially heated horizontal wall (vertical sidewall) configuration for both CWHF (first column) and CWT (second column) cases at Pr = 100. (a) RaCWT (RaCWHF) = 5 × 103 and (b) RaCWT (RaCWHF) = 5 × 104.

values of Rayleigh number for the CWHF bc. On the other hand, the values of Nu for both SW and HW configurations remain comparable for Newtonian and shear-thickening fluids when the Rayleigh number is high (see 5 × 104 and 1 × 105 cases in Figure 10) for both CWT and CWHF bcs. The mean Nusselt number Nu assumes greater values in the SW configuration for shear thinning fluids (i.e., n < 1) than in the HW configuration for a given set of n and Pr for high values of RaCWT in the CWT bc. It is also worth reiterating that the steady twodimensional solution has been obtained for a greater range of values of n in the SW configuration than in the HW configuration for a given set of values of Ra and Pr. The Nu correlations for the SW configuration for both CWT and CWHF bcs were proposed elsewhere by the present authors, and interested readers are referred to refs 15 and 16 for further information. The variation of θCWT (θCWHF) along the vertical midplane x1/L = 0.5 (horizontal midplane x2/L = 0.5) for different values

of Nu with n for a given set of values of Ra and Pr for CWT and CWHF bcs for both HW and SW configurations (data for SW taken from refs 15 and 16 for CWT and CWHF bcs, respectively) are shown in Figure 10, which demonstrates that Nu decreases with increasing values of n for both configurations. It is evident from Figure 10 that in the CWT bc Nu attains greater values in the SW configuration than in the HW configuration for both shear-thinning and shear-thickening fluids for small values of Rayleigh number (e.g., RaCWHF = 5 × 103). By contrast, Nu in the CWHF bc assumes greater values in the HW configuration for Newtonian and shear-thinning fluids for small values of Rayleigh number (e.g., RaCWHF = 5 × 103), whereas Nu assumes higher values in the SW configuration than in the HW configuration for shear-thickening fluids. However, Nu assumes higher values in the SW configuration than in the HW configuration for shear-thinning fluids, and the difference between the Nu values becomes more noticeable for small values of n and high 471

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

of RaCWT (RaCWHF) for n ≤ 1 and n ≥ 1 are shown in Figures 11 and 12, respectively, in the HW (SW) configuration for both CWT and CWHF bcs in order to understand and explain the differences in the Nu values between the SW and HW configurations. Figures 11 and 12 show that in the CWT bc the temperature gradient next to the active walls (hot and cold walls) for the HW configuration is smaller than in the SW configuration for shear-thinning fluids for the range of Rayleigh number considered here and also in Newtonian and shear-thickening fluids for small values of nominal Rayleigh number (see RaCWT = 5 × 103 and 104). However, the distributions of temperature between the active walls for both the HW and SW configurations remain similar for Newtonian and shear-thickening fluids for high values of nominal Rayleigh number (see RaCWT = 5 × 104 and 105). High values of thermal gradient next to the active walls, indicating a thin thermal boundary layer adjacent to the walls, lead to higher values of Nu ∼ L/δth for the SW configuration than in the HW configuration for shear-thinning fluids for the range of RaCWT considered here and also for shear-thickening fluids for small values of nominal Rayleigh number (see RaCWT = 5 × 103 and 104). The comparable thermal boundary layer thicknesses for the HW and SW configurations lead to similar values of Nu ∼ L/δth for high values of nominal Rayleigh number (see RaCWT = 5 × 104 and 105) in the case of the CWT bc. The temperature differences between the active walls (hot and cold wall) ΔT in the CWHF boundary configuration can be scaled as ΔT ∼ qδth/k; thus, a large (small) ΔT indicates a large (small) thermal boundary-layer thickness δth for the CWHF bc. Figures 11 and 12 show that the temperature difference between the active walls for the HW configuration is smaller (greater) than in the SW configuration for Newtonian and shear-thinning fluids (for shear-thickening fluids) for small values of Rayleigh number (e.g., RaCWHF = 5 × 103) leading to NuHW > NuSW (NuHW < NuSW) for n ≤ 1 (n > 1) (see Figure 10). However, the temperature difference between the active walls for both HW and SW configurations remains comparable for high values of Rayleigh number (see RaCWHF = 5 × 104 and 105) in the case of shear-thickening fluids (n > 1, especially the n = 1.8 case). Moreover, the temperature difference for the HW configuration is greater than in the SW configuration for shear-thinning fluids at RaCWHF = 5 × 104, which gives rise to NuHW < NuSW in Figure 10.

configuration for other combinations of Rayleigh and Prandtl numbers. The mean Nusselt number Nu has been found to be insensitive to the variation of nominal Prandtl number Pr range given by Pr = 10 −105 for both CWT and CWHF configurations. The differences between the mean Nusselt number in response to the bcs have been explained based on scaling and detailed physical arguments. Finally, new correlations for Nu for powerlaw fluids with n ranging from 0.6 to 1.8 have been proposed for both CWT and CWHF bcs guided by a detailed scaling analysis. The proposed correlations are shown to satisfactorily capture the variation of Nu with RaCWHF (RaCWT), Pr, and n for all cases considered in this analysis. It has also been demonstrated that Nu assumes greater values in the SW configuration than in the HW configuration for both shear-thinning and shear-thickening fluids for small values of Rayleigh number in the CWT bc. In contrast, Nu for shearthinning fluids (shear thickening fluids (i.e., n > 1)) in the differentially heated horizontal lower wall configuration remains greater (smaller) than the corresponding values obtained for the differentially heated vertical sidewall configuration for small values of RaCWHF in the case of CWHF bc. However, Nu attains higher values in the differentially heated vertical sidewall configuration than in the differentially heated horizontal wall configuration for shear-thinning fluids for high values of Rayleigh number in the CWHF bc. The values of Nu for both differentially heated vertical sidewall and horizontal lower wall configurations remain comparable for Newtonian and shear-thickening fluids for high values of Rayleigh number for both CWT and CWHF bcs. In contrast, greater values of Nu have been obtained in the SW configuration for shear thinning fluids than in the HW configuration for a given set of n and Pr for high values of RaCWT in the CWT bc. Detailed physical explanations have been provided for the observed behavior. Although the temperature dependence of thermophysical properties such as consistency K, power-law exponent n, and thermal conductivity k are not expected to modify the qualitative nature of the physics reported here,18,31 the temperature dependence of thermophysical properties might be necessary for accurate quantitative predictions of real non-Newtonian fluids. The effects of temperature dependence of thermophysical properties on the natural convection of power-law fluids will form the basis of a future study.



5. CONCLUSIONS Two-dimensional steady-state laminar natural convection of power-law fluids in square enclosures with horizontal walls heated from below has been numerically analyzed for the powerlaw index n in the range 0.6−1.8 and the Prandtl number range Pr = 10 to 105 for RaCWT = 103−105 and RaCWHF = 103−105 for both CWT and CWHF bcs. The effects of power-law exponent and nominal Rayleigh and Prandtl numbers on heat and momentum transport have been analyzed in detail. It has been found that the mean Nusselt number Nu increases with increasing (decreasing) values of nominal Rayleigh number (power-law exponent) due to strengthening of convective transport in both bcs. The mean Nusselt number Nu for the CWT bc remains greater than in the CWHF bc for a given set of numerical values of power-law exponent, nominal Rayleigh and Prandtl numbers in the case of Newtonian and shear-thinning fluids. However, mean Nusselt number Nu for the CWT bc are found to be smaller than in the CWHF configuration for a particular combination of RaCWHF and Pr for shear-thickening fluids, whereas Nu in the CWT bc is found to be greater than the corresponding value in the CWHF

AUTHOR INFORMATION

Corresponding Author

*N. Chakraborty. Phone: +44(0) 191 208 3570. Fax: +44(0) 191 208 8600. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



REFERENCES

(1) de Vahl Davis, G. Natural convection of air in a square cavity: A bench mark numerical solution. Int. J. Numer. Methods Fluids 1983, 3, 249. (2) Catton, I.; Ayyaswamy, P. S.; Clever, R. M. Natural convection flow in a finite, rectangular slot arbitrarily oriented with respect to the gravity vector. Int. J. Heat Mass Transfer 1974, 17, 173. (3) Ostrach, S. Natural convection in enclosures. J. Heat Transfer 1988, 110, 1175. (4) Ozoe, H.; Churchill, S. W. Hydrodynamic stability and natural convection in Ostwald-De Waele and Ellis fluids: the development of a numerical solution. AIChE J. 1972, 18, 1196.

472

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473

Industrial & Engineering Chemistry Research

Article

(5) Kim, G. B.; Hyun, J. M.; Kwak, H. S. Transient buoyant convection of a power law non-Newtonian fluid in an enclosure. Int. J. Heat Mass Transfer 2003, 46, 360. (6) Lamsaadi, M.; Naïmi, M.; Hasnaoui, M. Natural convection of nonNewtonian power law fluids in a shallow horizontal rectangular cavity uniformly heated from below. Heat Mass Transfer 2005, 41, 239. (7) Alloui, Z.; Ben Khelifa, N.; Beji, H.; Vasseur, P.; Guizani, A. The onset of convection of power-law fluids in a shallow cavity heated from below by a constant heat flux. J. Non-Newtonian Fluid Mech. 2013, 196, 70. (8) Solomatov, V.; Barr, A. C. Onset of convection in fluids with strongly temperature-dependent power-law viscosity. Phys. Earth Planet. Inter. 2006, 155, 140. (9) Albaalbaki, B.; Khayat, R. E. Pattern selection in the thermal convection of non-Newtonian fluids. J. Fluid Mech. 2011, 668, 500. (10) Ohta, M.; Akiyoshi, M.; Obata, E. A numerical study on natural convective heat transfer of pseudo-plastic fluids in a square cavity. Numer. Heat Transfer, Part A 2002, 41, 357. (11) Inaba, H.; Daib, C.; Horibe, A. Natural convection heat transfer of micro-emulsion phase-change materialslurry in rectangular cavities heated from below and cooled from above. Int. J. Heat Mass Transfer 2003, 46, 4427. (12) Khezzar, L.; Signiner, D.; Vinogradov, I. Natural convection of power law fluids in inclined cavities. Int. J. Thermal Sci. 2011, 53, 8−17. (13) Lamsaadi, M.; Naïmi, M.; Hasnaoui, M.; Mamou, M. Natural convection in a vertical rectangular cavity filled with a non-Newtonian power law fluid and subjected to a horizontal temperature gradient. Numer. Heat Transfer, Part A 2006, 49, 969. (14) Lamsaadi, M.; Naïmi, M.; Hasnaoui, M. Natural convection heat transfer in shallow horizontal rectangular enclosures uniformly heated from the side and filled with non-Newtonian power law fluids. Energy Convers. Manage. 2006, 47, 2535. (15) Turan, O.; Sachdeva, A.; Chakraborty, N.; Poole, R. J. Laminar natural convection of power-law fluids in a square enclosure with differentially heated side walls subjected to constant wall temperatures. J. Non-Newtonian Fluid Mech. 2011, 166, 1049. (16) Turan, O.; Sachdeva, A.; Poole, R. J.; Chakraborty, N. Laminar natural convection of power-law fluids in a square enclosure with differentially heated side walls subjected to constant wall heat flux. J. Heat Transfer Transactions of ASME 2012, 134, 122504. (17) Turan, O.; Sachdeva, A.; Poole, R. J.; Chakraborty, N. Effects of aspect ratio on natural convection of power-law fluids in rectangular enclosures with differentially heated vertical sidewalls. Int. J. Heat Mass Transfer 2013, 60, 722. (18) Kaddiri, M.; Naïmi, M.; Raji, A.; Hasnaoui, M. Rayleigh-Benard convection of non-Newtonian power-law fluids temperature-dependent viscosity. ISRN Thermodyn. 2012, 614712. (19) Turan, O.; Lai, J.; Poole, R. J.; Chakraborty, N. Laminar natural convection of power-law fluids in a square enclosure submitted from below to a uniform heat flux density. J. Non-Newtonian Fluid Mech. 2013, 199, 80. (20) Barth, W. L.; Carey, G. F. On a natural-convection benchmark problem in non-Newtonian fluids. Num. Heat Transfer, Part B 2006, 50, 193. (21) Leung, W. H.; Hollands, K. G. T.; Brunger, A. P. On a physicallyrealizable Benchmark problem in internal natural convection. Int. J. Heat Mass Transfer 1998, 41, 3817. (22) Turan, O.; Poole, R. J.; Chakraborty, N. Influences of boundary conditions on the aspect ratio effects in laminar natural convection in rectangular enclosures with differentially heated side walls. Int. J. Heat Fluid Flow 2012, 33, 131. (23) Turan, O.; Poole, R. J.; Chakraborty, N. Boundary condition effects on natural convection of Bingham fluids in a square enclosure with differentially heated horizontal walls. J. Comput. Thermal Sci. 2012, 4, 77. (24) Ng, M. L.; Hartnett, J. P. Natural convection in power law fluids. Int. Commun. Heat Mass Transfer 1986, 13, 115. (25) Patankar, S. V. Numerical Heat Transfer and Fluid Flow; Hemisphere: Washington, DC, 1980.

(26) Turan, O.; Chakraborty, N.; Poole, R. J. Laminar natural convection of Bingham fluids in a square enclosure with differentially heated side walls. J. Non-Newtonian Fluid Mech. 2010, 165, 901. (27) Ternik, P.; Rudolf, R.; Zunic, Z. Numerical study of RayleighBénard natural convection heat transfer characteristics of water-based Au nanofluids. Mater. Technol. 2013, 2, 211. (28) Ternik, P.; Rudolf, R. Laminar natural convection of nonNewtonian nanofluids in a square enclosure with differentially heated side walls. Int. J. Simul. Model. 2013, 12, 5. (29) Ternik, P.; Rudolf, R. Conduction and coznvection heat transfer characteristics of water-based Au nanofluids in a square cavity with differentially heated side walls subjected to constant temperatures. Therm. Sci. 2013, No. 10.2298/TSCI130604082T. (30) Bejan, A. Convection heat transfer; John Wiley Sons Inc.: New York, 1984. (31) Emery, A. F.; Lee, J. W. The effects of property variations on natural convection in a square cavity. J. Heat Transfer 1990, 121, 57.

473

dx.doi.org/10.1021/ie4023447 | Ind. Eng. Chem. Res. 2014, 53, 456−473