Laminar Natural Convection From an Isothermal Vertical Surface to Pseudoplastic and Dilatant Fluids

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, η₀. 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, η₀ 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.

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 η₀). 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.

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, T_w, 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.

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, η₀ 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.

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.

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

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 log₁₀(β*) 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.

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

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

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.

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 Nu_L) 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 Nu_L 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, Pr₀, 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)

          No: return to Step 2.b.ii.

       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 Nu_L and write it to file.

3. (3)

   Next β*…

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

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 log₁₀(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 log₁₀(β*)-range with strongly non-Newtonian fluids (i.e., |log₁₀(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.

With the code validated, runs were performed to evaluate the behavior of Nu_L as a function of n, R*, β*, Pr₀, 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; |log₁₀(R*)| = 1, 2, 3, and 4 with negative values indicating pseudoplastic fluids; log₁₀(β*) = 0 to log₁₀(1/R*) for pesudoplastic fluids with the range reversed for dilatant fluids, the ranges being divided into 100 equal increments; log₁₀(Pr₀) = −2.00 to 4.00 in increments of 1.00, and; log₁₀(Bo) selected for each combination of the other parameters such that the corresponding Rayleigh number, Ra* = Bo/Pr*, remained in the range 4.00 ≤ log₁₀(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 Nu_L 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 Nu_L 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 Nu_L are provided in Tables 1 and 2 for convenience.

Each Bo value produces a set of plots containing four families of curves each starting at Pr* = Pr₀, labeled “A”, and ending at Pr* = R* Pr₀ = Pr_∞. The heavy lines labeled “B” in each set are the Newtonian solutions (n = 1.00). Thus, the value of Pr₀ 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* = Pr₀, 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, Pr₀. For example, the families labeled “C” in Figs. 3 and 5 are for fluids with |log₁₀(R*)| = 3.00 (i.e., log₁₀(Pr_∞) = 1.00 and log₁₀(Pr₀) = 4.00 in Fig. (3), and log₁₀(Pr_∞) = 7.00 and log₁₀(Pr₀) = 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*, Nu_L 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 log₁₀(β*) 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 |log₁₀(R*)| = 4.00 is sufficiently large. To investigate whether this is true, Nu_L for these two points was also calculated at |log₁₀(R*)| = 2.00, 6.00, 8.00, and 10.00. In each case, the corresponding value of Pr₀ was shifted to ensure that the desired Pr* was in the middle of the log₁₀(β*) range. The fractional changes were then computed using the values at |log₁₀(R*)| = 10.00 as the basis for comparison, since this value of R* was found to produce power law asymptotes as noted earlier. The Nu_L 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 |log₁₀(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 log₁₀(β*) range may be expected to be within a few percent of power law behavior if the fluid complies with |log₁₀(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., log₁₀(R*) = 0.00, n = 1.00) with a viscosity of 1000 Ns/m². Curve B is the power law model, η_a = K$γ·n-1$⁠, for a fluid with n = 0.50 and K = 5000Nsⁿ/m². It should be noted that in log–log coordinates, the slope of this model is n−1 and its intercept is log₁₀(K). Curve C is the EMPL model, Eq. (1), for a strongly non-Newtonian, pseudoplastic fluid with log₁₀(R*) = −4.00. It otherwise has the same K and n as those in Curve B and the same η₀ 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 log₁₀(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 |log₁₀(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 log₁₀(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) |log₁₀(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 Nu_L value, 33.302, and the corresponding power law asymptote, 34.352 (taken as the Nu_L value when |log₁₀(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.

A larger |log₁₀(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. Nu_L should therefore be expected to be closer to its power law asymptote. The results show this to be the case: Nu_L = 34.235 when |log₁₀(R*)| = 6.00 and increases to 34.341 for |log₁₀(R*)| = 8.00. In contrast, a lower |log₁₀(R*)| would produce the opposite effects: a narrower region III, and a slope shallower than n−1. Thus, Nu_L should be expected to depart more from the associated power law value as it in fact does: Nu_L = 28.610 for |log₁₀(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 log₁₀(β*) range.

The current study maps the behavior of Nu_L 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, η₀ and η_∞, and the power law parameters, K and n, in their definitions. The calculated Nu_L 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.

Nu_L 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 η₀). Calculated values suggested that the second condition is approached when |log₁₀(R*)|≳4.00. Nu_L values were then found to be within a few percent of power law values and approached closer as |log₁₀(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 Nu_L values for weakly non-Newtonian fluids where no published data appear to exist. These data therefore provide a broad mapping of the behavior of Nu_L for the given problem with these types of fluids.

Nomenclature

Bo =

Boussinesq number

g =

gravitational acceleration (m/s²)

$Gr∧$ =

power law Grashof number defined by Acrivos [6]

h(x) =

local convection coefficient (W/m² K)

h_L =

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

k =

fluid's thermal conductivity (W/m K)

K =

power law fluid consistency (N · sⁿ/m²)

L =

surface length (m)

n =

power law flow index

Nu(x) =

local Nusselt number

Nu_L =

average Nusselt number based on h_L and L

P_∞ =

ambient pressure (Pa)

$Pr∧$ =

power law Prandtl number defined by Acrivos [6]

Pr₀ =

Prandtl number based on η₀

Pr_∞ =

Prandtl number based on η_∞

Pr* =

Prandtl number based on η*

R* =

limiting Newtonian viscosity ratio, η_∞/η₀

Ra* =

Rayleigh number based on η*

T =

temperature (K)

T_W =

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

α =

thermal diffusivity (m²/s)

β =

volumetric thermal expansion coefficient (1/K)

β₀ =

MPL shear rate parameter based on η₀

β_∞ =

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/m²)

η_P =

characteristic power law viscosity (N · s/m²)

η₀ =

zero shear rate Newtonian viscosity (N · s/m²)

η_∞ =

high shear rate Newtonian viscosity (N · s/m²)

η* =

characteristic EMPL viscosity (N · s/m²)

ρ_∞ =

fluid density far from the surface (kg/m³)

Superscript

+ =

denotes a dimensionless quantity