## Abstract

This paper reports the results of a numerical study of natural convective heat transfer from a vertical isothermal surface to pseudoplastic and dilatant fluids. The analysis calculates the average Nusselt number in the laminar regime when the surface is exposed to an otherwise quiescent fluid. Because the solution utilizes constitutive equations that are valid over the entire shear rate range, the results map the behavior regardless of the shear rates that exist in the flow field. The Nusselt number is shown to approach Newtonian values when the shear rates throughout the flow field are predominantly either in the zero or in the high shear rate Newtonian regions of the flow curve. When they are principally in the power law regime, and if the fluid is also strongly non-Newtonian, then the Nusselt number approaches power law values. For all other cases, it is seen to attain intermediate values. Furthermore, a shear rate parameter is identified that determines the shear rate regime where the system is operating. The average Nusselt number is presented in both graphical and tabular forms over a broad range of system parameters.

## Introduction

Natural convection from plane, vertical surfaces with Newtonian fluids is a fundamental problem that has been investigated both numerically and experimentally for various surface boundary conditions. For example, Ostrach [1] calculated the Nusselt number in the laminar regime for the isothermal boundary condition. He solved the governing equations by introducing a similarity transformation and then evaluated the resulting system of coupled ordinary differential equations numerically. Other such solutions from various sources have been compiled in books such as Gebhart et al. [2] who provide detailed analyses of many buoyancy-induced transport phenomena with both Newtonian and non-Newtonian fluids.

Natural convection with non-Newtonian fluids has also received considerable attention because of its importance in industrial applications where processes with fluids such as molten plastics, slurries, paints, blood, etc., are common. A review of published data for natural convection with non-Newtonian fluids may be found in Shenoy and Mashelkar [3] and in Irvine and Capobianchi [4], among other sources. In such cases, the fluid's rheology significantly complicates the problem because the apparent viscosity at any location in the flow field becomes, as a minimum, a function of the local shear rate, thus presenting a nonlinear problem. Consequently, a constitutive equation that is appropriate for the specific fluid under investigation must be used and models have been developed to describe various fluid behaviors. Examples include the power law (i.e., Ostwald–de Waele), modified power law (MPL), and extended modified power law (EMPL) models for pseudoplastic and dilatant fluids, the Carreau–Yasuda, Cross, and Ellis models for pseudoplastic fluids, the Casson and Herschel–Bulkley models for viscoplastic fluids, etc. The theory leading to such models is documented in a number of sources such as Macosko [5]. Regardless of the model selected, the additional parameters required to define the rheology typically make general solutions difficult to attain. Because the current work deals specifically with natural convection with pseudoplastic and dilatant fluids, it is important to understand their rheology and how they are modeled in order to provide a proper context for the cited works that follow.

Pseudoplastic (i.e., shear thinning) and dilatant (i.e., shear thickening) fluids are purely viscous, time-independent, non-Newtonian fluids with no initial yield stress. The former have an apparent viscosity that decreases with increasing shear rate whereas it increases in the latter. Their rheology is characterized by a flow curve (i.e., a curve of the apparent viscosity, ηa, versus shear rate, $γ·$) that generally has five shear rate regimes as shown for a pseudoplastic fluid in Curve C of Fig. 1. Near zero shear rates (region I, the zero shear rate Newtonian region), the apparent viscosity approaches a constant value, η0. As the shear rate is increased, the viscosity transitions smoothly (region II, the low shear rate transition region) until it reaches power law behavior (region III, the power law region). Increasing the shear rate beyond region III produces another smooth transition (region IV, the high shear rate transition region) from power law behavior to a different constant viscosity, η, when sufficiently high shear rates are attained (region V, the high shear rate Newtonian region). Thus, characterization of the rheology contributes at least four parameters to the problem: the limiting constant viscosities, η0 and η, and the two power law parameters, the fluid consistency, K, and the flow index, n. A larger and possibly different parameter set may be required depending on the particular choice of model selected for the fluid's rheology.

Fig. 1
Fig. 1
Close modal

Many of the published solutions for natural convection with pseudoplastic and dilatant fluids utilize the power law equation, ηa = K$γ·n-1$, to model the apparent viscosity. Examples include Acrivos [6] whose analytical solution considered natural convection from various isothermal surfaces. He calculated the Nusselt number by utilizing a similarity transformation for the limiting case when the Prandtl number approaches infinity. Tien [7] used an integral technique to obtain an approximate solution for the Nusselt number at large Prandtl numbers. He considered the laminar regime of a vertical flat surface with three different boundary conditions: isothermal, constant heat flux, and arbitrary temperature distribution. Kawase and Ulbrecht [8] used an integral technique to obtain an approximate solution at large Prandtl numbers for laminar and turbulent free convection from an isothermal vertical surface. Huang and Chen [9] utilized a local similarity solution to obtain local and average Nusselt numbers over a broad Prandtl number range for vertical surfaces with both isothermal and constant heat flux boundary conditions. It should be emphasized that such solutions are strictly valid only if the shear rates throughout the flow field are predominantly in region III of the flow curve (see Fig. 1). Furthermore, as will be demonstrated in this study, there is a second, additional criterion that must be satisfied for power law solutions to apply; the fluid must be strongly non-Newtonian (i.e., η must differ sufficiently from η0). This will be explained in detail and quantified in the subsequent sections.

A notable study that utilizes a more complete, albeit approximate, model for the apparent viscosity is that of Moulic and Yao [10]. They numerically calculated both the Nusselt number and the skin friction coefficient for the laminar regime of an isothermal vertical surface. Their simplified model for the flow curve included only regions I, III, and V and featured abrupt transitions between these regimes by extending them to where they intersect. Such a model, while computationally convenient, has two important disadvantages. First, for cases where the shear rates are predominantly in region II or IV, the absence of the transition regions means that the local apparent viscosity will either be under- or overestimated depending on the particular transition region within which the shear rates reside. This will inevitably lead to inaccuracies in the results in such instances. Second, as mentioned above, power law results are only approached when the shear rates are chiefly in region III and if the fluid is also strongly non-Newtonian. However, when the fluid is weakly non-Newtonian, it may be shown that the transition regions cause the rheology in region III to behave with an apparent K and n that differ from the physical values of K and n. Power law results are then never achieved and solutions always reside between the power law and Newtonian asymptotes even when all the shear rates are in region III. In contrast, the rheology in region III predicted by a model with abrupt transitions always behaves with the physical values of n and K even for weakly non-Newtonian fluids.

The present study calculates the average Nusselt number for natural convection from an isothermal, vertical flat surface to quiescent pseudoplastic or dilatant fluids in the laminar regime. Unlike prior studies, the current work utilizes constitutive equations that span the entire shear rate range and that include all five shear rate regimes. The results are therefore valid for whatever shear rates may exist in the flow field and demonstrate appropriate asymptotic behavior. Furthermore, a dimensionless shear rate parameter is identified that establishes the shear rate regime where the system is operating. Solutions are reported for values in all shear rate regimes including the transition regions where no prior data appear to exist. In addition, results are provided for weakly non-Newtonian fluids where again there appear to be no published solutions. These also permit the criterion for applying power law solutions to be quantified. The data are provided in both graphical and tabular forms over a broad range of system parameters.

## Analysis

### Problem Description and Analytical Model.

The geometry considered in this study, along with the definitions of the coordinates, x and y, and the corresponding velocity components, u and v, is shown in Fig. 2. A vertical, flat surface with constant and uniform temperature, Tw, is exposed to a semi-infinite expanse of quiescent fluid. Free convective heat transfer is allowed to occur between the surface and the fluid with no other heat transfer modes active. The fluid's temperature far from the wall is T and its density is ρ. All fluid properties are considered constant except for (1) viscosity, which is a function of shear rate as defined below and (2) density, which is assumed to vary only in those changes that drive the flow (i.e., the Boussinesq approximation). The surface is of infinite extent in the z-direction so that the problem is two-dimensional in x and y. Furthermore, the portion of the surface that is selected for study is wholly within the laminar regime of the boundary layer and extends in the x-direction a distance, L, from the surface's leading edge. The average Nusselt number was calculated along this section of the surface under steady-state conditions assuming the typical boundary layer approximations (e.g., see Ref. [2]), negligible viscous dissipation and pressure work, and no ambient pressure gradients.

Fig. 2
Fig. 2
Close modal
The fluid under investigation was either pseudoplastic or dilatant and its rheology was modeled using the EMPL models introduced by Capobianchi [11]. For pseudoplastic fluids, the EMPL expression for the apparent viscosity is given by
$ηa=η0-η∞1+η0Kγ·1-n+η∞$
(1)
The analogous EMPL expression for dilatant fluids is defined as
$ηa=η∞-η01+η∞Kγ·1-n+η0$
(2)

These models were selected because they describe the apparent viscosity throughout all five regions of the flow curve while requiring only the two limiting Newtonian viscosities, η0 and η, and the power law parameters, K and n, in their definitions. They are thus comparatively simple relative to other models that span the entire shear rate regime. More importantly, because they explicitly contain the power law parameters, they allow the investigation of the criteria necessary for power law solutions to apply.

These models were derived using the Dunleavy and Middleman [12] MPL model for pseudoplastic fluids and the analogous MPL model for dilatant fluids proposed by Capobianchi and Irvine [13]. The MPL models are only valid in regions I, II, and III and are therefore applicable to problems with strongly non-Newtonian fluids when all the shear rates are within those regions. It may be shown that for such fluids, the low shear rate transition regions predicted by the EMPL models match those of the MPL models. Furthermore, Eq. (1) has the same functional form as the Cross [14] model for pseudoplastic fluids which is also valid over the entire shear rate regime, although it was derived by a different approach and utilizes a different set of parameters. Thus, the accuracy of the EMPL equations is similar to that of the MPL and the Cross models.

Consider a differential parcel of fluid within the thermal boundary layer at location x where the thermal and hydrodynamic boundary layer thicknesses are δT and δH, respectively, noting that at all locations δTδH (see Fig. 2). With the above assumptions, the continuity, energy, and x-direction momentum equations may be shown to be, respectively
$∂u∂x+∂v∂y=0$
(3)
$u∂T∂x+v∂T∂y=α∂2T∂y2$
(4)
$u∂u∂x+v∂u∂y=1ρ∞∂∂y(ηa∂u∂y)+gβ(T-T∞)$
(5)
where T, α, g, and β are the temperature, thermal diffusivity, gravitational acceleration, and volumetric thermal expansion coefficient, respectively. The associated boundary conditions may be expressed as
$T(x,y=0)=TW T(x,y→∞)=T∞$
(6)
$u(x,y=0)=0 u(x,y→∞)=0$
(7)
$v(x,y=0)=0 ∂v∂y|(x,y→∞)=0$
(8)
ηa in Eq. (5) is given by either Eq. (1) or (2) depending on the specific type of fluid under investigation. However, the shear rate must be defined for the particular problem. For an incompressible, purely viscous, non-Newtonian fluid in simple shear, application of the typical boundary layer approximations reduce the general expression for the shear rate to
$γ·=|∂u∂y|$
(9)
Finally, the local convection coefficient, h(x), and the local Nusselt number, Nu(x), may be obtained by performing a local energy balance at the surface, which yields
$h(x)=-kTW-T∞∂T∂y|x,y=0; Nu(x)=h(x)xk=-x(TW-T∞)∂T∂y|x,y=0$
(10)
where k is the thermal conductivity of the fluid. h(x) may then be utilized to evaluate an average convection coefficient, hL, from x = 0 to L, and a corresponding average Nusselt number, NuL, based on L, giving
$hL=1L∫x=0Lh(x)dx→NuL=hLLk=-1TW-T∞∫x=0L∂T∂y|x,y=0dx$
(11)
Equations (1)–(9) allow the calculation of the velocity and temperature fields, and the latter may then be used to evaluate NuL from Eq. (11). However, to generalize the problem and to reduce the number of parameters (the above equation set contains a total of eleven parameters), it is useful to solve the equation set in dimensionless form. Thus, the following dimensionless variables are defined:
$x+=xL y+=yL$
(12)
$T+=T-T∞TW-T∞ u+=uu* v+=vu*$
(13)
$ηa+=ηaη* γ·+=γ·γ·*$
(14)

where the “+” superscript denotes a dimensionless quantity. Here, η* is a characteristic EMPL viscosity that depends on the fluid type, u* is a characteristic x-direction velocity, and $γ·*$ is a characteristic shear rate. Each of these must be defined before proceeding.

An expression for $γ·*$ may be obtained by nondimensionalizing Eq. (9) using the variables defined in Eqs. (12)–(14). Substituting these definitions into Eq. (9) gives
$γ·=|∂u∂y|=u*L|∂u+∂y+|∴ let γ·*=u*L→γ·+=γ·γ·*=|∂u+∂y+|$
(15)
Next, the η* functions are specified by applying the same definitions utilized by Capobianchi [11]. Specifically, the EMPL models, Eqs. (1) and (2), were used, replacing the power law apparent viscosity, K$γ·n-1$, in the denominators by the characteristic power law viscosity, ηP, defined below. For pseudolastic and dilatant fluids this yielded, respectively
$η*=η0-η∞1+η0ηP+η∞$
(16)
$η*=η∞-η01+η∞ηP+η0$
(17)
where
$ηP=Kγ·*n-1$
(18)

By doing so, η* equals the EMPL apparent viscosity corresponding to the given $γ·*$.

What remains, then, is to define u*. Selection of u* is complicated by the fact that there is no imposed velocity in the problem, the flow being driven by buoyancy forces that are balanced by shear forces. A suitable value for u* may be obtained by performing an order-of-magnitude analysis (OMA) on the governing equations to determine the order of magnitude of u. Such a study was performed by Capobianchi and Aziz [15] who showed that the order of magnitude of u is given by
$u~(Bo1+Pr*)1/2αL=u*$
(19)
where Bo is the Boussinesq number and where the Prandtl number from the Capobianchi and Aziz study is here replaced by Pr*, the Prandtl number based on η*. These are defined as
$Pr*=η*/ρ∞α; Bo=gβL3(TW-T∞)α2$
(20)
With the above terms fully specified, the governing equations, Eqs. (3)–(5), and their boundary conditions, Eqs. (6)–(8), may now be nondimensionalized giving, respectively, after simplification
$∂u+∂x++∂v+∂y+=0$
(21)
$u+∂T+∂x++v+∂T+∂y+=(1+Pr*Bo)1/2∂2T+∂y+2$
(22)
$u+∂u+∂x++v+∂u+∂y+=Pr*(1+Pr*Bo)1/2∂∂y+(ηa+∂u+∂y+) +Bo(1+Pr*Bo)T+$
(23)
$T+(x+,y+=0)=1 T+(x+,y+→∞)=0$
(24)
$u+(x+,y+=0)=0 u+(x+,y+→∞)=0$
(25)
$v+(x+,y+=0)=0 ∂v+∂y+|(x+,y+→∞)=0$
(26)
Next, the definitions of the dimensionless variables are introduced into Eq. (11) giving NuL in terms of the gradient of the dimensionless temperature profile at the surface of the wall
$NuL=hLLk=-∫x+=01∂T+∂y+|x+,y+=0dx+$
(27)
The expressions for the apparent viscosity, Eqs. (1) and (2) are now nondimensionalized following the procedure used by Capobianchi [11]. First, the MPL shear rate parameters, β0 and β, the EMPL shear rate parameter, β*, and the limiting Newtonian viscosity ratio, R*, are defined
$β0=η0ηP; β∞=η∞ηP; β*=1+β01+β∞; R*=η∞η0$
(28)
These are substituted into the equations for η*, Eqs. (16) and (17), giving the following expressions for pseudoplastic and dilatant fluids, respectively:
$η*=η0β*$
(29)
$η*=β*R*η0$
(30)

Note that β* establishes the value of the characteristic viscosity and thus the corresponding characteristic shear rate where the system is operating. It was therefore termed the EMPL shear rate parameter.

For pseudoplastic fluids, the valid range of β* is from 1 to 1/R*, with values at these endpoints indicating that the system is operating in the low shear rate and the high shear rate Newtonian regions, respectively. The analogous β*-range for dilatant fluids is 1/R* to 1. Thus, when β* is at one of these endpoints, solutions should be expected to reach the corresponding Newtonian values. Furthermore, for cases where the fluid is strongly non-Newtonian, Capobianchi [11] showed that the characteristic viscosity approaches the power law viscosity at the midpoint of the log10(β*) range. The system is then operating in the power law region, region III, and solutions should thus approach power law results. For all other cases, they should be expected to reside between the power law and Newtonian asymptotes.

Next, Eqs. (28)–(30) are introduced into Eqs. (1) and (2) along with the definitions of the other dimensionless variables yielding, for pseudoplastic and dilatant fluids, respectively
$ηa+=β*{1+β∞γ·+1-n1+β0γ·+1-n}$
(31)
$ηa+=1β*{1+β0γ·+1-n1+β∞γ·+1-n}$
(32)

where the expression for $γ·+$ is given by Eq. (15).

Also, the Prandtl numbers based on η0 and on η are defined as
$Pr0=η0/ρ∞α; Pr∞=η∞/ρ∞α=R*Pr0$
(33)
Then, for each fluid type, the expression for Pr*, Eq. (20), may be written in terms Pr0 by substituting Eq. (33) with Eq. (29) or (30) into Eq. (20), yielding the expression for Pr* for pseudoplastic and dilatant fluids, respectively, as
$Pr*=Pr0β*$
(34)
$Pr*=β*R*Pr0$
(35)

Thus, noting that β0 and β may be calculated from β* and R*, there remain five parameters when the problem is cast in dimensionless form: β*, R*, n, Pr0, and Bo. The equation set was solved numerically using a finite volume method over a broad range of these five parameters as described below.

Prior to proceeding, one final piece of information is necessary to perform the computations efficiently. Because of the definition of x+, Eq. (12), the domain in the x+-direction is of fixed, finite size, 0 ≤ x+ ≤ 1. However, it is semi-infinite in the y+-direction. Thus, the computational upper limit of y+ must be specified and chosen sufficiently large so that T+ and u+ are essentially unchanged from their quiescent values at this boundary. Since δH depends in some fashion on Pr* and Bo, selecting an acceptable maximum y+ for the worst-case condition and then utilizing that domain for all cases means that often numerous unnecessary nodes will be included. To preclude this problem, the y+-domain should be selected as a function of δH. This requires knowledge of the functional relationship between δH and the parameters Pr* and Bo which may be found by performing an OMA. The Capobianchi and Aziz [15] study provides the order of magnitude of δH/L as well as that of δT/L as
$δHL~(1+Pr*Bo)1/4[Pr*1/2+1];δTL~(1+Pr*Bo)1/4;δHδT~[Pr*1/2+1]$
(36)

Equation (36) was used to determine the order of magnitude of δH/L and this quantity utilized to establish the y+-domain as described in the next section. Equation (36) demonstrates appropriate asymptotic behavior. When Pr*≪1, the hydrodynamic boundary layer should approach the thermal boundary layer so that δH/δT → 1, as properly suggested by Eq. (36). And when Pr *≫ 1, δH/δT → Pr*1/2 as computed by Bejan [16] and as Eq. (36) correctly predicts.

### Numerical Model and Algorithm.

A finite volume scheme was utilized to solve the above set of dimensionless equations numerically. Prior to discretization, the size of the domain in the y+-direction was established as ten times the quantity δH/L, with δH/L evaluated from Eq. (36) for the particular Pr* and Bo numbers being calculated. The values and y+-direction gradients of u+, v+, and T+ at this border were then checked to ensure that the domain was adequately sized. The domain was discretized using a rectangular grid with zero-thickness control volumes at the boundaries. While uniform grid spacing was utilized in the x+-direction, an expanding grid was used in the y+-direction to improve computational efficiency. This is because u+ and T+ remain essentially unchanged from their quiescent values for y+δH/L. Therefore, to ensure accuracy in calculating the T+-gradient at the surface (required to find NuL) while precluding an unnecessarily fine mesh far from the surface, uniform grid spacing was selected from y+ = 0 to δH/L and expanded linearly thereafter. The spacing in the uniform grid portion of the mesh was determined by first calculating the order of magnitude of δT at the x+-location of the first row of nodes (i.e., by replacing L in Eq. (36) by the x-coordinate of the first row of nodes), and then dividing this quantity by a user-supplied factor. With this scheme, a grid refinement in the x+-direction automatically refined the grid in the y+-direction ensuring sufficiently fine grid spacing near the wall for all rows of nodes.

The above set of dimensionless equations was discretized in finite volume form following the work of Patankar [17]. Because the boundary layer assumptions make x+ a one-way coordinate, the domain was computed using a marching scheme starting at the first row of nodes and progressing one row at a time in the positive x+-direction. When discretized for a single row of nodes, the dimensionless energy and momentum equations with boundary conditions, Eqs. (22)–(26), each yields a set of algebraic equations that is tridiagonal when cast in matrix form. They may thus be conveniently solved using the tridiagonal matrix algorithm (see Ref. [17] for details), bearing in mind that iteration is nevertheless required because of the nonlinearity of the problem. In contrast, the dimensionless continuity equation, Eq. (21), when discretized for a single row of nodes produces an explicit algebraic equation for each node.

Finally, care must be exercised in the numerical evaluation of NuL from Eq. (27) because the discontinuity in T+ at x+ = 0 yields an infinite y+-direction gradient at that location. The numerical integration was therefore performed using a trapezoidal quadrature formula that excluded the low-limit endpoint of the domain (Press et al. [18]). The numerical consequences of this problem were thus precluded.

The discretized set of dimensionless equations was then solved with the following algorithm:

1. (1)

Read the input parameters: β*-range and step size, R*, n, Pr0, Bo, and the mesh and convergence control parameters.

2. (2)

Loop for each value of β* in the β*-range:

• (a)

Calculate the run parameters (i.e., Pr*, δH/L, etc.) and the mesh parameters (i.e., node spacing, etc.)

• (b)

Loop for each row of nodes:

• (i)

Guess the u+ distribution.

• (ii)

Calculate the v+ distribution from the continuity equation using the latest u+ distribution.

• (iii)

Calculate the T+ distribution from the energy equation using the latest u+ and v+ distributions.

• (iv)

Calculate the u+ distribution iteratively from the momentum equation using the latest v+ and T+ distributions.

• (v)

Have the u+, v+, and T+ distributions converged?

1. (1)

2. (2)

Yes: calculate Nu(x) and write it to file along with the u+, v+, and T+ distributions.

• (c)

Next row of nodes…

• (d)

Calculate NuL and write it to file.

3. (3)

Next β*…

What follows are the results of the computations along with a discussion of the salient findings.

## Results and Discussion

### Program Validation.

Prior to performing the computations, the program was validated. First, a mesh refinement study was performed to ensure that results were essentially independent of further mesh refinements. A 500 × 10 × 10.0 mesh was found to be adequate, with successive mesh refinements producing changes no greater than 0.05%. Here, the three values that specify the mesh indicate the following: first quantity—the number of control volumes in the x+-direction; second quantity—the number of control volumes in the y+-direction within the thermal boundary layer at the first row of nodes, and; third quantity—the y+-direction gradient of the control volume expansion in the region y+ > δH/L.

Next, with the mesh defined, the accuracy of the code was verified by calculating cases for which analytical solutions exist. First, Newtonian asymptotes were calculated at the endpoints of the β*-range for various values of Pr* and Bo with n = 0.50 and 1.50. These were compared to the Newtonian values from Ostrach [1] and those reported in Gebhart et al. [2]. Computed values were lower in all cases and were within 0.12% when compared to both studies except at log10(Pr*) = −2.00 where they were within 0.77% when compared to Ostrach, and within 0.22% when compared to Gebhart et al.

Second, power law asymptotes were computed at the midpoint of the log10(β*)-range with strongly non-Newtonian fluids (i.e., |log10(R*)| = 10.00). These were compared to Acrivos′ [6] solution which is valid for $Pr∧$ and $Gr∧$1/(n+1)$Pr∧$2n/(3n+1) ≫ 1, where $Pr∧$ and $Gr∧$ are power law Prandtl and Grashof numbers, respectively, defined by Acrivos (Pr* and Bo may each be shown to be functions of both $Pr∧$ and $Gr∧$). All computed values were found to be within 1.2% of Acrivos′ data, generally improving with increasing $Pr∧$ and $Gr∧$. Results are likely closer than 1.2% because evaluating Acrivos′ values requires the use of data that are reported in graphical form thus limiting the accuracy of the computation.

### Numerical Solution.

With the code validated, runs were performed to evaluate the behavior of NuL as a function of n, R*, β*, Pr0, and Bo. The selected ranges for these parameters were n = 0.50–1.50 in increments of 0.10, noting that n = 1.00 implies a Newtonian fluid; |log10(R*)| = 1, 2, 3, and 4 with negative values indicating pseudoplastic fluids; log10(β*) = 0 to log10(1/R*) for pesudoplastic fluids with the range reversed for dilatant fluids, the ranges being divided into 100 equal increments; log10(Pr0) = −2.00 to 4.00 in increments of 1.00, and; log10(Bo) selected for each combination of the other parameters such that the corresponding Rayleigh number, Ra* = Bo/Pr*, remained in the range 4.00 ≤ log10(Ra*) ≤ 9.00. This requirement ensures that the flow is sufficiently energetic for the boundary layer assumptions to be valid while remaining within the laminar regime. Considering all of the combinations of the above parameters, more than 59,000 values of NuL were computed in order to map its behavior.

The results of these calculations are shown graphically in Figs. 3 and 4 for pseudoplastic fluids, and in Figs. 5 and 6 for dilatant fluids, which are plots of NuL versus Pr*, recalling that Pr* is a function of β* (see Eqs. (34) and (35)). Corresponding points on Figs. 3 and 5 have been annotated to describe some of the important features of the solution. Also, representative numerical values of NuL are provided in Tables 1 and 2 for convenience.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal
Table 1

Representative numerical values of NuL for pseudoplastic fluids

NuL

log10(R*) = −4.00

log10(R*) = −2.00

log10(β*)
log10(Bo)log10(Pr0)n0.8001.6002.0002.4003.2000.4000.8001.0001.2001.600
4.000.000.507.39107.82757.90107.93607.95806.73817.27387.40217.48257.5646
0.707.06937.67957.81147.88667.94826.40897.02587.21957.35877.5231
0.906.76177.48687.67747.79897.92076.14826.75276.98347.17177.4403
−1.000.507.86647.96477.97967.98647.99047.66947.80617.83307.84897.8644
0.707.73997.91847.95357.97297.98797.49427.71487.77227.81057.8526
0.907.58227.84277.90427.94217.97917.31517.57687.66437.73157.8214
6.002.000.5013.49919.39621.46522.78723.82710.02912.80813.91914.84016.150
0.7011.60917.07019.53521.45523.4878.992611.43212.61613.73715.630
0.9010.40915.23217.70219.88222.8018.337110.38011.48412.62514.905
1.000.5020.27923.61224.29424.63424.85216.57519.41420.23020.79121.421
0.7018.10022.31323.46324.16224.75614.87917.83618.97919.89221.097
0.9016.38720.82022.32123.36824.49313.72016.34217.56218.67720.492
8.004.000.5015.55326.61933.30240.01549.68310.83514.41315.91317.21519.243
0.7012.70521.03826.51232.83345.7869.478512.45413.99215.52418.367
0.9011.09917.71622.18727.57240.7818.677511.06212.42613.90917.188
3.000.5029.07646.33354.46560.78067.04120.31926.32928.61030.47333.155
0.7023.07936.82944.81552.79264.47917.23922.49425.05727.50531.749
0.9019.66330.66937.57045.17859.83915.42119.58221.90524.38229.653

NuL

log10(R*) = −4.00

log10(R*) = −2.00

log10(β*)
log10(Bo)log10(Pr0)n0.8001.6002.0002.4003.2000.4000.8001.0001.2001.600
4.000.000.507.39107.82757.90107.93607.95806.73817.27387.40217.48257.5646
0.707.06937.67957.81147.88667.94826.40897.02587.21957.35877.5231
0.906.76177.48687.67747.79897.92076.14826.75276.98347.17177.4403
−1.000.507.86647.96477.97967.98647.99047.66947.80617.83307.84897.8644
0.707.73997.91847.95357.97297.98797.49427.71487.77227.81057.8526
0.907.58227.84277.90427.94217.97917.31517.57687.66437.73157.8214
6.002.000.5013.49919.39621.46522.78723.82710.02912.80813.91914.84016.150
0.7011.60917.07019.53521.45523.4878.992611.43212.61613.73715.630
0.9010.40915.23217.70219.88222.8018.337110.38011.48412.62514.905
1.000.5020.27923.61224.29424.63424.85216.57519.41420.23020.79121.421
0.7018.10022.31323.46324.16224.75614.87917.83618.97919.89221.097
0.9016.38720.82022.32123.36824.49313.72016.34217.56218.67720.492
8.004.000.5015.55326.61933.30240.01549.68310.83514.41315.91317.21519.243
0.7012.70521.03826.51232.83345.7869.478512.45413.99215.52418.367
0.9011.09917.71622.18727.57240.7818.677511.06212.42613.90917.188
3.000.5029.07646.33354.46560.78067.04120.31926.32928.61030.47333.155
0.7023.07936.82944.81552.79264.47917.23922.49425.05727.50531.749
0.9019.66330.66937.57045.17859.83915.42119.58221.90524.38229.653
Table 2

Representative numerical values of NuL for dilatant fluids

NuL

log10(R*) = 4.00

log10(R*) = 2.00

log10(β*)
log10(Bo)log10(Pr0)n−3.200−2.400−2.000−1.600−0.800−1.600−1.200−1.000−0.800−0.400
8.000.001.5029.51619.67916.06813.1318.926136.90430.22927.69025.54922.392
1.3032.39321.48817.45214.1629.414340.15332.81729.88027.32123.272
1.1036.12423.98719.37715.60110.09443.90936.25732.90829.85624.644
−1.001.5044.57531.83426.45921.87715.17852.81145.91343.10740.64936.905
1.3048.98735.03229.02423.86216.15557.09049.68646.39343.35338.249
1.1054.23439.33432.55326.63117.55261.51154.41650.79947.20940.455
10.002.001.5030.87219.98616.21813.2058.947140.01631.70128.75826.35622.936
1.3034.26321.93617.67414.2739.444544.22834.77731.29228.36423.906
1.1038.90124.70119.73715.78310.14249.46339.07234.93531.32625.446
1.001.5051.79033.66927.36522.32815.31067.03653.56948.86345.08439.908
1.3058.28337.53030.29024.50416.33574.90259.34453.61148.81641.641
1.1067.19043.06834.49827.63417.82084.68567.52760.59954.52344.592
12.004.001.5030.95920.00516.22713.2108.948440.24131.80028.83026.41222.977
1.3034.41121.96917.69014.2819.446844.58834.93631.40528.44823.960
1.1039.17824.76719.77015.79910.14650.06439.35535.13431.46925.524
3.001.5052.31333.78327.42022.35615.31868.35954.17649.31245.43540.174
1.3059.14137.72530.38624.55216.34876.93660.27254.27849.31541.965
1.1068.73343.44434.68627.72817.84587.96569.10261.71755.32845.034

NuL

log10(R*) = 4.00

log10(R*) = 2.00

log10(β*)
log10(Bo)log10(Pr0)n−3.200−2.400−2.000−1.600−0.800−1.600−1.200−1.000−0.800−0.400
8.000.001.5029.51619.67916.06813.1318.926136.90430.22927.69025.54922.392
1.3032.39321.48817.45214.1629.414340.15332.81729.88027.32123.272
1.1036.12423.98719.37715.60110.09443.90936.25732.90829.85624.644
−1.001.5044.57531.83426.45921.87715.17852.81145.91343.10740.64936.905
1.3048.98735.03229.02423.86216.15557.09049.68646.39343.35338.249
1.1054.23439.33432.55326.63117.55261.51154.41650.79947.20940.455
10.002.001.5030.87219.98616.21813.2058.947140.01631.70128.75826.35622.936
1.3034.26321.93617.67414.2739.444544.22834.77731.29228.36423.906
1.1038.90124.70119.73715.78310.14249.46339.07234.93531.32625.446
1.001.5051.79033.66927.36522.32815.31067.03653.56948.86345.08439.908
1.3058.28337.53030.29024.50416.33574.90259.34453.61148.81641.641
1.1067.19043.06834.49827.63417.82084.68567.52760.59954.52344.592
12.004.001.5030.95920.00516.22713.2108.948440.24131.80028.83026.41222.977
1.3034.41121.96917.69014.2819.446844.58834.93631.40528.44823.960
1.1039.17824.76719.77015.79910.14650.06439.35535.13431.46925.524
3.001.5052.31333.78327.42022.35615.31868.35954.17649.31245.43540.174
1.3059.14137.72530.38624.55216.34876.93660.27254.27849.31541.965
1.1068.73343.44434.68627.72817.84587.96569.10261.71755.32845.034

Each Bo value produces a set of plots containing four families of curves each starting at Pr* = Pr0, labeled “A”, and ending at Pr* = R* Pr0 = Pr. The heavy lines labeled “B” in each set are the Newtonian solutions (n = 1.00). Thus, the value of Pr0 locates the families on the Newtonian fluid curves and all families begin and end at Newtonian asymptotes. This should be expected since the endpoints are at the limits of the shear rate parameter, β*, indicating that the system is operating in either region I or V depending on the endpoint. Thus, at the endpoint labeled “A” where Pr* = Pr0, the characteristic shear rate is in the zero shear rate Newtonian region, whereas it is in the high shear rate Newtonian region at the opposing endpoint where Pr* = Pr. Also, the corresponding value of R* is obtained as the ratio of the end value of Pr* = Pr, to its start value, Pr0. For example, the families labeled “C” in Figs. 3 and 5 are for fluids with |log10(R*)| = 3.00 (i.e., log10(Pr) = 1.00 and log10(Pr0) = 4.00 in Fig. (3), and log10(Pr) = 7.00 and log10(Pr0) = 4.00 in Fig. (5)).

There are five individual curves in each family that correspond to five different values of n. For pseudoplastic fluids, the upper curve is for n = 0.50 and the subsequent curves moving towards the Newtonian curve are for n = 0.60, 0.70, 0.80, and 0.90. For dilatant fluids, the lowest curve is for n = 1.50 and the subsequent curves moving towards the Newtonian curve are for n = 1.40, 1.30, 1.20, and 1.10. At any Pr*, NuL for pseudoplastic fluids is always greater than the corresponding Newtonian value except, of course, at the endpoints where it is equal. The opposite is true for dilatant fluids.

As stated earlier, when the shear rate parameter, β*, is at the midpoint of the log10(β*) range, the system may be shown to be operating in region III. Thus, systems with β* between the midpoint and an endpoint indicate that their characteristic shear rates are in the corresponding transition region. Here solutions are seen to be between the Newtonian and power law asymptotes, smoothly approaching the asymptotic values. At the midpoint, solutions will approach power law asymptotes only if the fluid is then also strongly non-Newtonian. The physical reasons for this will be explained in detail below. For example, the points on the n = 0.50 and 1.50 curves at the locations labeled “D” would approach the power law asymptotes at their corresponding Pr* if |log10(R*)| = 4.00 is sufficiently large. To investigate whether this is true, NuL for these two points was also calculated at |log10(R*)| = 2.00, 6.00, 8.00, and 10.00. In each case, the corresponding value of Pr0 was shifted to ensure that the desired Pr* was in the middle of the log10(β*) range. The fractional changes were then computed using the values at |log10(R*)| = 10.00 as the basis for comparison, since this value of R* was found to produce power law asymptotes as noted earlier. The NuL values were found to be within −16.8%, −3.1%, −0.4%, and −0.03% for pseudoplastic fluids, and within 4.6%, 0.5%, 0.05%, and 0.004% for dilatant fluids, for |log10(R*)| = 2.00, 4.00, 6.00, and 8.00, respectively. Hence, although a more complete investigation is required to quantify more concretely the limiting value of R* where power law behavior occurs, it appears that systems operating at the midpoint of the log10(β*) range may be expected to be within a few percent of power law behavior if the fluid complies with |log10(R*)|≳4.00. Otherwise, power law values are simply never reached.

To understand this behavior, it is useful to examine how R* affects the rheology of the fluid when all other properties are kept constant. Figure 1 shows the flow curves of four related fluids. Curve A is a Newtonian fluid (i.e., log10(R*) = 0.00, n = 1.00) with a viscosity of 1000 Ns/m2. Curve B is the power law model, ηa = K$γ·n-1$, for a fluid with n = 0.50 and K = 5000Nsn/m2. It should be noted that in log–log coordinates, the slope of this model is n−1 and its intercept is log10(K). Curve C is the EMPL model, Eq. (1), for a strongly non-Newtonian, pseudoplastic fluid with log10(R*) = −4.00. It otherwise has the same K and n as those in Curve B and the same η0 as the Newtonian viscosity of Curve A. It may be observed that region III of Curve C is coincident with Curve B and therefore has the same slope and intercept (i.e., values of n and K) as Curve B. In other words, shear rates that are in region III of Curve C will produce the same apparent viscosity as that predicted by the power law model.

In contrast, Curve D is for a weakly non-Newtonian fluid with log10(R*) = −1.00. All other properties are the same as those in Curve C. The slope in region III of Curve D fails to match that of the power law model. It is shallower such that it corresponds to a larger value of n and a lower value of K than the corresponding physical values. The opposite would be true for weakly non-Newtonian dilatant fluids. This occurs because region III is insufficiently wide to allow the low shear rate transition region, region II, to establish power law behavior before the high shear rate transition region, region IV, begins transitioning back to Newtonian behavior. Regions II, III, and IV then blend together into a continuous, smooth transition between regions I and V.

It should be emphasized that this behavior is physical rather than being an artifact of the chosen model. Consider, for example, a strongly non-Newtonian, pseudoplastic fluid. As |log10(R*)| is decreased and approaches 0 while all else remains constant, the fluid's rheology approaches that of a Newtonian fluid for which the slope in region III is 0 and, consequently, n = 1.00. Thus, while the apparent n for a strongly non-Newtonian fluid approaches the physical n, it must increase towards 1.00 as log10(R*) → 0.00 in order for the fluid to approach Newtonian behavior.

As a specific illustration of these points, consider the solution that is on the curve for n = 0.50 located at Point D in Fig. 3. The corresponding dimensionless shear rate profiles at the midpoint and at the end of the plate (i.e., at locations x+ = 0.500 and 1.000) are plotted in Fig. 7 along with the associated dimensionless apparent viscosities from the EMPL model, Eq. (31), and from the power law model. It may be observed that both requirements for reaching power law behavior are approached: (1) |log10(R*)| is sufficiently large to produce an apparent viscosity in region III that approaches that of the power law model (i.e., the slope in region III approaches n−1) and (2) the shear rates within the flow field are predominantly within region III. The latter is true throughout most of the profile except in the local region of the cusp (i.e., corresponding to the y+-location where u+ reaches its maximum), which stretches into the lower shear rate regimes, and immediately adjacent to the wall where the shear rates extend slightly into the leading edge of region IV. The difference between the resulting NuL value, 33.302, and the corresponding power law asymptote, 34.352 (taken as the NuL value when |log10(R*)| = 10.00), may be attributed to these departures of the shear rate from region III and to any difference that may exist between the actual slope in region III and n−1.

Fig. 7
Fig. 7
Close modal

A larger |log10(R*)| would make region III wider thus encompassing more of the shear rates within the flow field, and would drive the slope closer to n−1. NuL should therefore be expected to be closer to its power law asymptote. The results show this to be the case: NuL = 34.235 when |log10(R*)| = 6.00 and increases to 34.341 for |log10(R*)| = 8.00. In contrast, a lower |log10(R*)| would produce the opposite effects: a narrower region III, and a slope shallower than n−1. Thus, NuL should be expected to depart more from the associated power law value as it in fact does: NuL = 28.610 for |log10(R*)| = 2.00 (i.e., Point E in Fig. 4).

This then is the physical behavior that explains why weakly non-Newtonian fluids never reach power law behavior even when the system is operating in the middle of the log10(β*) range.

## Conclusions

The current study maps the behavior of NuL over a broad range of parameters for laminar natural convection between an isothermal vertical surface and quiescent pseudoplastic and dilatant fluids. The apparent viscosity was modeled using the EMPL models, Eqs. (1) and (2), which span the entire shear rate regime while requiring only the limiting Newtonian viscosities, η0 and η, and the power law parameters, K and n, in their definitions. The calculated NuL values are therefore valid for whatever shear rates exist in the flow field. Furthermore, a shear rate parameter was defined, β*, whose value indicates the characteristic shear rate for the particular flow situation and, consequently, the flow regime where the system is operating. Results were provided in all shear rate regimes including the two transition regions, regions II and IV, where no data seem to exist in the literature.

NuL values were found to range between the Newtonian and the power law asymptotes, the former occurring when the characteristic shear rate was either in region I or V. However, one of the key findings was that two criteria must be satisfied in order to reach power law results: the shear rates must be predominantly in region III, and the fluid must be strongly non-Newtonian (i.e., η must differ sufficiently from η0). Calculated values suggested that the second condition is approached when |log10(R*)|≳4.00. NuL values were then found to be within a few percent of power law values and approached closer as |log10(R*)| was increased. This behavior is physical and a consequence of the blending of regions II, III, and IV into a continuous transition region when the fluid is weakly non-Newtonian. Thus, this study also provided NuL values for weakly non-Newtonian fluids where no published data appear to exist. These data therefore provide a broad mapping of the behavior of NuL for the given problem with these types of fluids.

### Nomenclature

Nomenclature

• Bo =

Boussinesq number

•
• g =

gravitational acceleration (m/s2)

•
• $Gr∧$ =

power law Grashof number defined by Acrivos [6]

•
• h(x) =

local convection coefficient (W/m2 K)

•
• hL =

average convection coefficient from x = 0 to L (W/m2 K)

•
• k =

fluid's thermal conductivity (W/m K)

•
• K =

power law fluid consistency (N · sn/m2)

•
• L =

surface length (m)

•
• n =

power law flow index

•
• Nu(x) =

local Nusselt number

•
• NuL =

average Nusselt number based on hL and L

•
• P =

ambient pressure (Pa)

•
• $Pr∧$ =

power law Prandtl number defined by Acrivos [6]

•
• Pr0 =

Prandtl number based on η0

•
• Pr =

Prandtl number based on η

•
• Pr* =

Prandtl number based on η*

•
• R* =

limiting Newtonian viscosity ratio, η/η0

•
• Ra* =

Rayleigh number based on η*

•
• T =

temperature (K)

•
• TW =

surface temperature (K)

•
• T =

fluid temperature far from the surface (K)

•
• u =

velocity component in the x-direction (m/s)

•
• u* =

characteristic x-direction velocity (m/s)

•
• v =

velocity component in the y-direction (m/s)

•
• x, y, z =

coordinate directions (m)

### Greek Symbols

Greek Symbols

• α =

thermal diffusivity (m2/s)

•
• β =

volumetric thermal expansion coefficient (1/K)

•
• β0 =

MPL shear rate parameter based on η0

•
• β =

MPL shear rate parameter based on η

•
• β* =

EMPL shear rate parameter

•
• $γ·$ =

shear rate (1/s)

•
• $γ·*$ =

characteristic shear rate (1/s)

•
• δH =

hydrodynamic boundary layer thickness (m)

•
• δT =

thermal boundary layer thickness (m)

•
• ηa =

apparent viscosity (N · s/m2)

•
• ηP =

characteristic power law viscosity (N · s/m2)

•
• η0 =

zero shear rate Newtonian viscosity (N · s/m2)

•
• η =

high shear rate Newtonian viscosity (N · s/m2)

•
• η* =

characteristic EMPL viscosity (N · s/m2)

•
• ρ =

fluid density far from the surface (kg/m3)

### Superscript

Superscript

• + =

denotes a dimensionless quantity

## References

1.
Ostrach
,
S.
,
1953
, “
An Analysis of Laminar Free-Convection Flow and Heat Transfer About a Flat Plate Parallel to the Direction of the Generating Body Force
,” National Advisory Committee for Aeronautics, Report No. 1111, pp.
63
79
.
2.
Gebhart
,
B.
,
Jaluria
,
Y.
,
Mahajan
,
R. L.
, and
Sammakia
,
B.
,
1988
,
Buoyancy-Induced Flows and Transport
, Reference Edition,
Hemisphere Publishing Corp.
,
New York
, pp.
1
132
, 857–890, Chaps. 1–3, 16.
3.
Shenoy
,
A. V.
, and
Mashelkar
,
R. A.
,
1982
, “
Thermal Convection in Non-Newtonian Fluids
,”
,
15
, pp.
143
225
.10.1016/S0065-2717(08)70174-6
4.
Irvine
,
T. F.
, Jr.
, and
Capobianchi
,
M.
,
2005
, “
Non-Newtonian Fluids—Heat Transfer
,”
The CRC Handbook of Mechanical Engineering
, 2nd ed.,
F.
Kreith
and
Y.
Goswami
, eds.,
Boca Raton
, pp.
4-269
4-278
, Chap. 4.9.
5.
Macosko
,
C. W.
,
1994
,
Rheology Principles, Measurements, and Applications
,
Wiley-VHC, Inc.
,
New York
, pp.
1
175
, Part I, Chaps. 1–4.
6.
Acrivos
,
A.
,
1960
, “
A Theoretical Analysis of Laminar Natural Convection Heat Transfer to Non-Newtonian Fluids
,”
Am. Inst. Chem. Eng. J.
,
6
, pp.
584
590
.10.1002/aic.690060416
7.
Tien
,
C.
,
1967
, “
Laminar Natural Convection Heat Transfer From Vertical Plate to Power-Law Fluid
”,
Appl. Sci. Res.
,
17
(
3
), pp.
233
248
.10.1007/BF00386093
8.
Kawase
,
Y.
, and
Ulbrecht
,
J. J.
,
1984
, “
Approximate Solution to the Natural Convection Heat Transfer From a Vertical Plate
,”
Int. Commun. Heat Mass Transfer
,
11
(
2
), pp.
143
155
.10.1016/0735-1933(84)90018-6
9.
Huang
,
M.-J.
, and
Chen
,
C.-K.
,
1990
, “
Local Similarity Solutions of Free Convective Heat Transfer From a Vertical Plate to Non-Newtonian Power Law Fluids
,”
Int. J. Heat Mass Transfer
,
33
(
1
), pp.
119
125
.10.1016/0017-9310(90)90146-L
10.
Moulic
,
S. G.
, and
Yao
,
L. S.
,
2009
, “
Non-Newtonian Natural Convection Along a Vertical Flat Plate With Uniform Surface Temperature
,”
ASME J. Heat Transfer
,
131
, p.
062501
.10.1115/1.3090810
11.
Capobianchi
,
M.
,
2008
, “
Pressure Drop Predictions for Laminar Flows of Extended Modified Power Law Fluids in Rectangular Ducts
,”
Int. J. Heat Mass Transfer
,
51
, pp.
1393
1401
.10.1016/j.ijheatmasstransfer.2007.11.019
12.
Dunleavy
,
J. E.
, Jr.
, and
Middleman
,
S.
,
1966
, “
Correlation of Shear Behavior of Solutions of Polyisobutylene
,”
Trans. Soc. Rheol.
,
10
(
1
), pp.
157
168
.10.1122/1.549055
13.
Capobianchi
,
M.
, and
Irvine
,
T. F.
, Jr.
,
1992
, “
Predictions of Pressure Drop and Heat Transfer in Concentric Annular Ducts With Modified Power Law Fluids
,”
Heat Mass Transfer
,
27
(
4
), pp.
209
215
.10.1007/BF01589918
14.
Cross
,
M. M.
,
1965
, “
Rheology of Non-Newtonian Fluids: A New Flow Equation for Pseudoplastic Systems
,”
J. Colloid Sci.
,
20
, pp.
417
437
.10.1016/0095-8522(65)90022-X
15.
Capobianchi
,
M.
, and
Aziz
,
A.
,
2012
, “
A Scale Analysis for Natural Convective Flows Over Vertical Surfaces
,”
Int. J. Therm. Sci.
,
54
, pp.
82
88
.10.1016/j.ijthermalsci.2011.11.009
16.
Bejan
,
A.
,
1984
,
Convection Heat Transfer
,
John Wiley and Sons, Inc.
,
New York
, pp.
114
122
.
17.
Patankar
,
S. V.
,
1980
,
Numerical Heat Transfer and Fluid Flow
,
Hemisphere Publishing Company
,
New York
, pp.
11
109
, Chaps. 2–5.
18.
Press
,
W. H.
,
Teukolsky
,
S. A.
,
Vetterling
,
W. T.
, and
Flannery
,
B. P.
,
1992
,
Numerical Recipes in FORTRAN: The Art of Scientific Computing
, 2nd ed.,
Cambridge University Press
,
Cambridge, UK
, pp.
129
130
, Chap. 4.