Ind. Eng. Chem. Res. 2002, 41, 6393-6401
6393
A Reinterpretation of the Turbulent Prandtl Number† Stuart W. Churchill‡ Department of Chemical Engineering, University of Pennsylvania, 311A Towne Building, 220 South 33rd Street, Philadelphia, Pennsylvania 19104
The eddy viscosity at any location in fully developed turbulent flow in a round pipe is shown by an exact analysis to be simply the ratio of the shear stress due to the time-averaged turbulent fluctuations in the velocity to that due to the molecular motion and thus independent of its heuristic diffusional origin. The eddy conductivity in fully developed turbulent convection is similarly shown to be simply the ratio of the corresponding contributions to the radial heat flux density. However, while providing a quantitative physical rationale for the eddy viscosity and the eddy conductivity, the analysis that leads to this insight at the same time eliminates their conceptual value except in a historical sense. On the other hand, the turbulent Prandtl number, as redefined directly in terms of the time-averaged fluctuations, remains an essential parametric variable. The local fraction of the radial heat flux density due to turbulence (the replacement for the eddy conductivity as a variable) appears, on the basis of both experimental evidence and asymptotic analyses, to be independent of geometry and the thermal boundary condition(s) and thereby a universal function only of the molecular Prandtl number and the local fraction of the shear stress due to the turbulence. It follows that the turbulent Prandtl number shares this restricted functionality. Although the functional and numerical behavior of the turbulent Prandtl number remains to this day incompletely defined either experimentally or theoretically, the rate of heat transfer, as represented by the Nusselt number, is fortuitously insensitive to these uncertainties. Introduction The time-averaged once-integrated differential representation for the conservation of momentum in the radial direction for fully developed turbulent flow of a fluid with invariant physical properties in a round tube may be expressed as
du τ)µ - Fu′v′ dy
τ ) (µ + µt)
du dy
or
(
)
j ) -k
∂T + FcT′v′ ∂y
or
-k 1 +
(1)
(3)
and
The new terms such as -Fu′v′ that arise from time averaging are sometimes called the Reynolds stresses because Reynolds1 was the first to space average the partial differential equations of conservation in order to obtain a more tractable representation. In onedimensional flows, the contribution of the turbulent fluctuations to the local shear stress (the rate of momentum transfer in the negative radial direction), namely, Fu′v′, has generally been modeled in terms of an eddy viscosity µt, or a total viscosity µT, thereby converting eq 1 to
µt du µ 1+ µ dy
have a lesser variance and could be more readily predicted or more simply correlated than -Fu′v′ itself. The corresponding expressions for the transport of energy in the negative radial direction, with the additional postulate of negligible viscous dissipation, are
or
µT
du dy (2)
j ) -(k + kt)
∂T ∂y
(
)
kt ∂T k ∂y
or -kT
∂T (4) ∂y
Equation 4 has generally been reexpressed as
[ ( )( ) ]
j ) -k 1 +
cµ kt µt ∂T k cµt µ ∂y
or
[ ( )( ) ]
-k 1 +
Pr µt ∂T (5) Prt µ ∂y
or as
( )( )( )
kT cµ µT ∂T cµT k µ ∂y
( )( )
Pr µT ∂T PrT µ ∂y
This model was proposed more than a century ago by Boussinesq2 with the hope and expectation that µt would
j ) -k
† This paper is dedicated to William R. Schowalter not only in recognition of his technical and educational contributions but also in appreciation of his leadership in the National Academy of Engineering and the American Institute of Chemical Engineers. His presence in those councils was reassuring at critical junctures. ‡ Tel.: 215-898-5579. Fax: 215-573-2093. E-mail: Churchil@ seas.upenn.edu.
with the hope and expectation that the eddy conductivity kt or the total conductivity kT will be better behaved than FcT′v′ and, in turn, that the turbulent Prandtl number, cµt/kt, and/or the total Prandtl number, cµT/kT, will be even better behaved in terms of variance, prediction, and correlation. The expectations and hopes relative to the behavior of µt, µT, kt, kT, Prt, and PrT are generally
or
10.1021/ie011021k CCC: $22.00 © 2002 American Chemical Society Published on Web 07/30/2002
-k
(6)
6394
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002
accepted as being fulfilled, as attested by their continued usage to this day, although some purists have cast scorn on this entire heuristic diffusional concept because of its lack of a mechanistic rationale. Until recently, the modeling of fully developed turbulent flow has generally been based on eq 2 with an empirical correlating equation for µt. The experimental values of µt upon which such a correlating equation is based are subject to considerable uncertainty because of their determination from (-dp/dx)(r/2)/(du/dy) - µ and in particular the determination of du/dy from measured values of u{y}. The modeling of fully developed turbulent convection has correspondingly been based on eq 5 with an additional correlating equation or constant value for Prt and often the additional postulate of a constant value or a linear variation of the heat flux density with radius. The experimental values of kt upon which a correlating equation for Prt is based are subject to even greater uncertainty than those of µt because of 2 their determination from j/(-∂T/∂y) - k ) [(Fc/r)∫r0 u 2 (∂T/∂x) dr ]/(-∂T/∂y) - k and in particular the determination of ∂T/∂x from T{x,y}. The turbulent Prandtl number, as determined from the ratio of experimental values of µt and kt, is obviously subject to even greater uncertainty, particularly for large values of Pr, for which the temperature profile develops very near the wall, and for small values of Pr, for which kt is very small relative to k. Theoretical results, experimental results, and correlating equations for Prt were recently reviewed objectively and in considerable detail by Kays,3 and to avoid duplication, that review will be accepted with only one minor reservation as a point of reference with respect to the state of the art. Furthermore, the primary objective of the analysis presented here is to reinterpret the results reviewed by Kays rather than to reevaluate them. A New Representation for Turbulent Flow and Convection. Prompted by the appearance of essentially exact computed values of the turbulent shear stress by direct numerical simulation (DNS), Churchill and Chan4 readdressed the question of the merits in terms of modeling, prediction, and correlation of the introduction of the eddy viscosity in eq 1 and were surprised to discover, on the favorable side, that this quantity could be expressed directly in terms of the time average of the fluctuating velocities without recourse to a diffusive mechanism involving a gradient of the time-averaged velocity, while, on the unfavorable side, that this quantity served no unique useful function for flow in a round tube or between parallel plates and was singular at one radial position and negative over an adjacent finite range of the radius in circular concentric annuli. Accordingly, they suggested abandoning the eddy viscosity concept and simply expressing eq 1 in “wall” units, namely, as
1-
y+ du+ ) + (u′v′)+ a+ dy+
(7)
where (u′v′)+ ) -Fu′v′/τw. They further devised a generalized correlating equation for (u′v′)+ as a function of y+ and a+. Churchill5 subsequently concluded that (u′v′)++ ≡ -Fu′v′/τ, namely, the local fraction of the shear stress due to fluctuations in velocity, was an even better dimensionless variable for the modeling of tur-
bulent flow. Equation 1 then takes the form
(
1-
)
y+ du+ [1 - (u′v′)++] ) + + a dy
(8)
If τ were a dependent variable, the expression (u′v′)++ ) -Fu′v′/τ might appear to be the equivalent of a similarity transformation. However, for fully developed flow in a round tube, τ bears a known linear relationship to the independent variable τw. The term -pu′v′/τ is simply a convenient representation for -Fu′v′/τw(1 y/a). Equation 8 may be integrated formally in terms of R ) 1 - y+/a+ to obtain the following expression for the velocity distribution:
u+ )
a+ 2
∫R1 [1 - (u′v′)++] dR2 2
(9)
This integral expression may, in turn, be integrated over the cross section and then reduced by means of integration by parts to obtain the following single integral for the mixed-mean velocity and the Fanning friction factor, namely, f ≡ 2τw/Fum2:
2 1/2 ) um+ ≡ f
()
+
∫01u+ dR2 ) a4 ∫01[1 - (u′v′)++] dR4 (10)
The “unity” term of the integrand of eqs 9 and 10 can, of course, be integrated analytically, but that does not constitute a significant simplification. These integral expressions are only marginally simpler than their counterparts in terms of µt, but the reduction to a single integral for um+ does not follow so directly or obviously in terms of µt. More importantly, eliminating the derivative between eqs 1 and 2 reveals that
µt (u′v′)++ -Fu′v′ ) ) µ τ + Fu′v′ 1 - (u′v′)++
(11)
This is a very significant result because of the revelation that µt/µ is simply the ratio of the contributions of the turbulent fluctuations and the molecular motions to the transport of momentum and thereby completely independent of its heuristic diffusional origin involving the time-mean velocity gradient. As a consequence of eqs 9-11, there does not appear to be any future role for the eddy viscosity except in a historical sense. It should not be inferred from eq 11 that µt is proportional to µ because (u′v′)++ will, in general, be an implicit function of µ. The analogue of (u′v′)++ for the transport of energy is (T′v′)++ ) FcT′v′/j, namely, the local fraction of the radial heat flux density due to the turbulent fluctuations. Introduction of this dimensionless quantity as well as the less familiar counterpart of u+, namely, T+ ≡ k(τwF)1/2(Tw - T)/µjw, into eq 3 results in
dT+ j [1 - (T′v′)++] ) jw dy
(12)
Then from eqs 4 and 12 and the definition of T+, it follows that
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002 6395
kt (T′v′)++ FcT′v′ ) ) k j - FcT′v′ 1 - (T′v′)++ and, in turn, from eqs 5 and 6 that
[
]
(13)
(
)
Prt (u′v′)++ 1 - (T′v′)++ -u′v′ j - FcT′v′ ) ) Pr (T′v′)++ 1 - (u′v′)++ cT′v′ τ + Fu′v′
(14)
and
[
]
PrT 1 - (T′v′)++ τ j - FcT′v′ ) ) Pr 1 - (u′v′)++ j τ + Fu′v′
(15)
Eliminating (T′v′)++ between eqs 14 and 15 indicates that ++
(u′v′) 1 ) PrT Prt
+
1 - (u′v′)++ Pr
(16)
and provides for interchangeability between Prt and PrT. Equations 13-15 are very significant results in the same sense as eq 11 in that kt/k, Prt/Pr, and PrT/Pr are revealed to be simple functions of (u′v′)++ and (T′v′)++ and thereby independent of their heuristic diffusive origin. Thus, the eddy conductivity, as well as the eddy viscosity, no longer appears to have any essential role except possibly in a historical sense. On the other hand, the turbulent Prandtl number and the total Prandtl number, as redefined by eqs 14 and 15, respectively, have continued essential utility even if their diffusionally based names are misleading. The total Prandtl number divided by the Prandtl number is seen from eq 15 to be equal to the fractional transport of energy by molecular motion divided by its counterpart for the transport of momentum, while the turbulent Prandtl number divided by the Prandtl number is seen from eq 14 to be equal to the ratio of the radial transport of momentum due to the turbulent fluctuations to that due to the molecular motion divided by the corresponding ratio for the radial transport of energy. Again, it should not be inferred from eq 13 that kt is proportional to k or from eqs 14 and 15 that Prt and PrT are proportional to Pr because the right-hand side of these expressions will generally be a function of Pr as well (u′v′)+. An Alternative but Slightly Inferior Interpretation. Rather than avoiding the eddy viscosity on the basis of eq 8 or eliminating it from existing models for flow on the basis of eq 11, the latter expression could equally well be interpreted as a confirmation of the intuition of Boussinesq in recognizing or postulating that the transport of momentum by the turbulent fluctuations is a diffusive process in terms of the radial gradient of the time-averaged velocity even though he had no idea that the effective diffusivity might be expressed exactly and completely in terms of a welldefined physical property, namely, -Fu′v′/τ. The same argument applies to the eddy conductivity by virtue of eqs 12 and 13. It then follows that the turbulent Prandtl number is indeed analogous to the true Prandtl number with the diffusional coefficients for the turbulent fluctuations substituted for the molecular ones and thereby properly named. Furthermore, eqs 9 and 10 as well as all of the ensuing analysis herein could be expressed with equal validity in terms of µt and kt rather than of
(u′v′)++ and (T′v′)++. It is remarkable that the critics who were scornful of the eddy viscosity, eddy conductivity, and turbulent Prandtl number concepts because of their lack of a mechanistic rational did not appear to recognize this simple justification. On the other hand, the many analysts who applied these concepts without apology may be given the benefit of the doubt. Having said the above, the eddy viscosity and eddy conductivity, although shown retroactively to be valid concepts, are inferior to the fractions of the shear stress and the radial heat flux density due to turbulence as operative variables in three respects. First, the momentum and energy balances in terms of the fluctuating quantities are more transparent and understandable because of the absence of heuristic quantities. Second, the eddy viscosity is singular and negative in all channels with unequal shear stresses on two bounding surfaces, such as in a concentric circular annulus or in a parallel-plate channel with surfaces of unequal roughness. The quantity (u′v′)++ is also inapplicable in this case, but (u′v′)+ is well behaved. Third, the formal integral solutions for the velocity distribution and the friction factor, as illustrated by eqs 9 and 10, and for the temperature distribution and Nusselt number, as illustrated subsequently, are much simpler in terms of the fractions of the shear stress and the radial heat flux density than in terms of the eddy viscosity and eddy conductivity. Implementation of the New Model for Convection. Substituting for (T′v′)++ in eq 12 from eq 14 results in
[ (
++ j Pr (u′v′) ) 1+ jw Prt 1 - (u′v′)++
)]
dT+ dy+
(17)
Although eq 17 incorporates two parametric variables, namely, Prt/Pr and (u′v′)++ as compared to one in eq 12, namely, (T′v′)++, it is more tractable with respect to numerical and functional evaluations because of the well-documented functionality of (u′v′)++ and the very constrained variance of Prt/Pr except in regimes and regions in which it has a very limited influence. For reasons that are not related to the turbulent Prandtl number, it is convenient, following Reichardt,6 to introduce into eq 17 the perturbation-like quantity γ defined by
1+γ)
() (
)
j y+ j τw ) / 1- + jw τ jw a
thereby obtaining
( )[ (
(1 + γ) 1 -
++ y+ Pr (u′v′) ) 1 + Prt 1 - (u′v′)++ a+
(18)
)]
dT+ (19) dy+
Equation 19 may be integrated formally in terms of R to obtain the following analogue of eq 9:
T+ )
a+ 2
∫R
(1 + γ) dR2
1 2
(
++ Pr (u′v′) 1+ Prt 1 - (u′v′)++
)
(20)
and in turn T+, weighted by u+/um+, can be integrated over the cross section to obtain
6396
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002
2a+ ) Tm+ ≡ Nu
( ) +
∫01T+ uu +
[
dR2 )
m
∫R
a+ 2
∫01
(1 + γ) dR2
[
1 2
++
]] (
Pr (u′v′) Prt 1 - (u′v′)++
1+
)
u+ dR2 (21) u m+
Nu f 0.076343
Equation 20 is more complex than its analogue, eq 9, by virtue of 1 + γ and Prt/Pr, and eq 21 relative to eq 10 by virtue of u+/um+ as well, a mathematical reflection of the greater physical complexity of convection as compared to flow. It may be noted that the temperature distribution, as given by eq 20, is not directly analogous to the velocity distribution as given by eq 9 for Pr ) 1, as has often been postulated, nor even for Pr ) Prt, owing to the nonlinear heat flux density distribution as represented by γ. Implementation of eq 21 to obtain numerical values of Nu, as described by Bo et al.7 and others, is beyond the scope and objective of this analysis, but the examination of eq 21 itself is within that objective because it helps to identity the role of the function labeled Prt in the prediction of Nu. For the special case of uniform heating or cooling at the wall of a round tube,
1+γ)
1 R2
∫0R
2
( )
u+ dR2 um+
8 (1 + γ)2 dR4
∫01 1+
[
]
1 R2
+
( ) +
∫0R uu + TT + 2
m
dR2
Pr Prt*
1/3
f 2
Re
1/2
(25)
T+ )
a+ 2
)
( )
PrT y+ dT+ [1 - (u′v′)++] ) + + Pr a dy
( )
∫R1 (1 + γ)[1 - (u′v′)++] 2
2a+ a+ ) Tm+ ) Nu 2
[
PrT dR2 Pr
( )
∫01 ∫R1 (1 + γ)[1 - (u′v′)]++ 2
(26)
(27)
]
PrT dR2 Pr
( ) +
u dR2 (28) + um
and
The dependence of Nu on both Pr and Prt then stems only from the term Pr/Prt in the denominator of the integral. On the other hand, for the other common thermal boundary conditions, namely, that of a uniform wall temperature,
1+γ)
(
(1 + γ) 1 -
(23)
++ Pr (u′v′) Prt 1 - (u′v′)++
( ) ()
The existence of these three special solutions for Nu is very fortuitous in view of the considerable uncertainty concerning Prt, both functionally and numerically in a general sense, as documented by Kays.3 Surprisingly, a significant reduction in that uncertainty has not occurred in the intervening 7 years since his review, except possibly in the followup by Weigand et al.8 It should be noted, as implied above, that these several conclusions regarding the dependence of Nu on Pr and Prt could also have been drawn if µt/µ were substituted for (u′v′)++/[1 - (u′v′)++] in eqs 17, 19-21, and 23. Equations 19-21 and 23 may be reexpressed in terms of PrT rather than Prt as follows:
(22)
and eq 21 may be reduced by integration by parts to obtain the following single integral:
Nu )
the wall where (u′v′)++/[1 - (u′v′)++] can be approximated by 0.7 (y+/10)3 and γ can be set equal to zero. Then, insofar as Prt approaches a finite value Prt*as y+ f 0, eq 19 can be integrated in closed form to obtain an algebraic expression for T+, from which it follows that
(24)
m
Because of the presence of T+/Tm+ in the integrand of eq 24, γ now depends on Prt, Nu depends on Prt by virtue of γ as well as by its explicit presence, reduction of eq 21 to a single integral by integration by parts is no longer possible, and an iterative scheme is required for the solution of eqs 20, 21, and 24 or their differential counterparts. The dependence of Nu on Prt may be seen to vanish from eqs 21 and 23 for Pr ) 0, for which condition Nu is thereby a function only of the rate of flow and the thermal boundary condition. Prt also vanishes from these two expressions for Pr ≡ Prt insofar as the this value of Prt is invariant with radius, a condition which, remarkably enough, appears to occur to at least a good degree of approximation. However, despite the consequent great simplification of eqs 21 and 23 for this particular value of Pr, Nu remains a function of Pr in the region encompassing that value. For Pr f ∞, the temperature field develops almost completely very near
Nu ) 8/
∫01(1 + γ)[1 - (u′v′)++]
( )
PrT dR4 Pr
(29)
Equations 26-29 appear to have a simpler structure than eqs 19-21 and 23, and they do have some utility for that reason in a qualitative sense. However, this impression of simplicity is somewhat misleading. For example, Nu does not approach zero as Pr f ∞, as might be inferred from eq 29. Also, Prt is a more sensitive variable than PrT, and, perhaps for that reason, all correlating equations appear to have been constructed for Prt rather than PrT. Alternative Models and Variables for the Prediction of Turbulent Flow and Convection. The mixing length l of Prandtl,9 as defined by
Fl 2
(du dy )
2
) -Fu′v′
(30)
has been widely used in the past as an alternative to the eddy viscosity, and a thermal mixing length has sometimes been used as an alternative to the eddy conductivity. Although the use of the mixing length has been fading in recent years, it was frequently cited in the past as being conceptually superior to the eddy viscosity because of its mechanistic origin as an analogue, however strained, of the mean free path in the kinetic theory of gases. In a round tube the mixing
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002 6397
length is readily shown to be related to (u′v′)++ as follows:
(l +)2 )
(
(u′v′)++
)
y+ 1 - + [1 - (u′v′)++]2 a
(31)
Because (u′v′)++ is finite at the center line, l +, and thereby l, may be seen from eq 31 to be unbounded at that location, thereby negating their very basis. This serious failing was apparently overlooked for over half a century because of (1) insufficiently precise data for u′v′ and u in that region, (2) uncritical examination of such data, (3) a fortuitous insensitivity of calculations for f and Nu to values of l in that region, and (4) the valid prediction, based on the experimentally observed linearity of l to y in “the turbulent core near the wall,” of a semilogarithmic dependence of u on y in that region. However, the semilogarithmic relationship has since been derived by speculative and asymptotic dimensional analysis without recourse to a mechanistic model such as the mixing length. Furthermore, the corresponding relationship for the velocity distribution in “the turbulent core near the center line”, which was derived on the basis of a finite value of the mixing length at the center line has since been shown to be in serious error. (See, for example, work by Churchill.5) The mixing length is also singular at the central plane between parallel plates and shares the failure of the eddy viscosity in channels with differing shear stresses on the walls. The κ- model of Launder and Spalding10 was designed to predict µt and thereby turbulent flow, but after 20 years its implementation remains quite arbitrary and approximate overall and essentially correlative rather than predictive near the wall. The κ--u′v′ model was utilized successfully by Hanjalic´ and Launder11 to predict the behavior in an annulus, which an eddy viscosity model cannot do. The κ--u′v′-T′v′ model utilized by Suzuki et al.12 for an annulus might, in principle, predict Prt. However, both of these latter models are quite empirical overall and essentially correlative near the walls. The large eddy simulation (LES) model pioneered by Schumann13 is predictive, in principle, for u′v′ and T′v′ and thereby Prt in the turbulent core but incorporates the eddy viscosity in the subgrid formulation for the region very near the wall. Kawamura et al.14 used the LES model for convection as well as for flow, but they postulated a value of Prt ) 0.85 for the subgrid region rather than predicted it. One methodology for the prediction of u′v′ and T ′v′ and thereby the prediction or avoidance of Prt for the entire cross section of a channel without any empiricism is that of DNS, which was apparently first conceived by Orszag and Patterson15 in the context of homogeneous and isotropic turbulence. However, the great early promise of this concept has not yet been completely fulfilled in that numerical results for confined shear flows are limited primarily to parallel-plate channels and to rates of flow barely above the minimum value of 145 for b+ for fully developed turbulence. The predicted values of several investigators, including Kim et al.,16 Lyons et al.,17 and Rutledge and Sleicher,18 for (u′v′)++ near the wall are in good agreement with experimental
measurements for all rates of flow, but the disappearance for b+ < 300 of the semilogarithmic regime of the time-mean velocity distribution, which characterizes the flow in the “turbulent core near the wall” (30 < y+ < 0.1b+), renders questionable the extrapolation of the results of the DNS calculations for the turbulent core to higher values of b+. The effective restriction of DNS calculations to parallel-plate channels is less serious in that the analogy of MacLeod19 implies that the values of (u′v′)++ and u+ are the same function of y+ and a+ in a round tube as of y+ and b+ in a parallel-plate channel. On the other hand, because of their potential precision, the execution of otherwise identical DNS calculations for these two geometries for a wide range of values of b+ and a+ would provide a more severe test of the accuracy and limitations of that analogy than do the best current experimental measurements. Calculations using DNS for heat transfer have been carried out for a wide range of values of Pr, but the various calculated values of T′v′, as contrasted with those of u′v′, are not in complete agreement with one another or with experimental results. In particular, the values of Sct (the analogue for mass transfer of Prt for heat transfer), as predicted in the pioneering work of Papavassiliou and Hanratty20 using a Lagrangian form of DNS for very large values of Sc (analogous to Pr), contradict the postulate of a limiting value of Sct and thereby of Prt for y+ f 0 that leads to eq 23. Because the most reliable experimental data for heat transfer generally agree with eq 24, including the 1/3 power dependence on Pr for large values of Pr, this discrepancy may be due to (1) the upper limit of Pr ) 100 for all fluids, whereas values of Sc are virtually unlimited in magnitude, (2) a fundamental difference between the turbulent transport of energy and species very near the wall, or (3) the limitation of the pertinent DNS calculations to b+ ) 150. Renormalization group theory (RGT) has been used to predict PrT directly rather than by detailed solution of the equations of conservation of momentum and energy for u′v′ and T′v′ as with DNS. Results of great purported generality are achieved at the price of idealizations of unknown consequence. Yahkot et al.21 and others have devised the following structure for this purpose:
( ) ( ) 1 -R PrT 1 -R Pr
(1+R)/(1+2R)
R/(1+2R) 1 +1+R PrT µ ) 1 µ + µt +1+R Pr or 1 - (u′v′)++ (32)
where
([
(
R) 1+81+
2 d
1/2
)]
)
- 1 /2
(33)
An iterative method of solution is obviously required for PrT. Prt can, in turn, be calculated by means of eq 16. Yahkot et al. proposed a value of 21/3 ) 7 for d in eq 33, but in a “Note Added in Proof” stated that 9/3 ) 3 might be better. Elperin et al.22 and Itazu and Nagano23 both proposed 14/3. On the other hand, an arbitrary value of d ) 24/3 ) 8 produces better agreement with experimental data and the subsequently mentioned empirical correlating equations for Prt. Yahkot et al. argue that eqs 32 and 33 are free of any “experimentally
6398
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002
Table 1. Comparison of Predictions of Prt for a+ ) 5000 and y+ ) 500 Pr Pr(u′v′)++/ [1 - (u′v′)++]a eq 35 (N&S) eq 34b (J&R) eq 36 (Kays) eqs 32 and 33 (Yahkot et al.)c a
0.0001 0.01857
0.001 0.1857
0.01 1.857
0.1 18.57
0.7 130.0
0.8633 160.3
1.0 185.7
10 1857
100 18570)
1201 115.9 38.54 42.98
63.80 12.35 4.619 4.855
3.366 2.000 1.227 1.216
1.003 0.9650 0.8877 0.8957
0.9457 0.8664 0.8553 0.8643
0.9458 0.8633 0.8544 0.8633
0.9463 0.8615 0.8538 0.8627
0.9525 0.8512 0.8504 0.8593
0.9553 0.8501 0.8500 0.8590
(u′v′)++ is based on eq 5 of Bo et al.7
b
With a coefficient of 0.0115. c Based on d ) 8.
adjusted parameters”, the quantity d apparently being exempted as a “physical dimension”. In any event, a number of the idealizations made in the derivations of eqs 32 and 33, and in particular that of homogeneous turbulence, would appear to confine their applicability to the turbulent core (y+ > 30) and perhaps to “the turbulent core near the wall” (30 < y+ < ∼0.1a+). Functional Dependence and Correlating Equations for Prt Abbrecht and Churchill24 determined experimentally the eddy conductivity as well as the eddy viscosity in the developing temperature field generated by a step in wall temperature in the fully developed flow of air in a round tube. Their derived values of kt are somewhat less certain than those of µt because of the necessity of evaluating the longitudinal as well as the radial temperature gradients, but the invariance on the mean of kt and Prt with axial distance is unambiguous evidence of their independence from the temperature field and thereby of the thermal boundary condition. The close correspondence of their values of Prt with those measured by Page et al.25 in fully developed turbulent flow between parallel plates at different uniform temperatures provides not only a confirmation of the independence of Prt from the thermal boundary condition but also an indication of independence from geometry for a+ ) b+. Abbrecht and Churchill thereby concluded that Prt is a function only of µt/µ and Pr. In retrospect, this conclusion may or may not be valid for the regime very near the wall because of the increased uncertainty of the values they determined for µt and kt as y+ f 0. This conclusion, that is, the dependence of Prt on µt/µ [or (u′v′)++] and Pr only, has since been dismissed as both obvious and disputed and apparently has never been proven or disproven conclusively, although it is implied by eq 32 and most empirical correlating equations. Kays,3 with some justification, questioned the applicability of this restricted dependence to the region near the wall and near the center line on the grounds of insufficient evidence and, having overlooked the work of Abbrecht and Churchill, speculated on the basis of a comparison of some separate experimental results for uniform heating and uniform wall temperature that Prt may depend on the thermal boundary condition. The extreme difference in the thermal boundary conditions utilized in the work of Abbrecht and Churchill and in that utilized by Page et al. would appear to be a more critical and sufficient test. Turning to correlating equations for Prt , Jischa and Rieke26 concluded that the dependence of Prt on y+ and Re, while real, is of second-order and that the simple expression
Prt ) 0.85 +
0.015 Pr
(34)
provides a reasonable representation for all values of Pr. Hill27 (also see Notter and Sleicher28) concluded from theoretical considerations involving isotropic turbulence and a linear temperature gradient that Prt approaches inverse proportionality to the equivalent of Pr (u′v′)++/ [1 - (u′v′)++] as Pr f 0. Notter and Sleicher developed the following correlating equation based on that result and the apparently false premise that Prt f 1 as (u′v′)++/[1 - (u′v′)++] f ∞:
Prt ) {1 + φ}/{0.025Pr(u′v′)++/[1 - (u′v′)++] + φ} 10 1+ (35) ++ 35 + (u′v′) /[1 - (u′v′)++]
{
}
where φ ) 90Pr3/2[(u′v′)++/[1 - (u′v′)++]]1/4. Kays found that most experimental and computed values of Prt in the turbulent core (y+ > 30) could be represented satisfactorily with an empirical expression that may be rewritten in terms of (u′v′)++ rather than ut/µ as follows:
Prt ) 0.85 +
[
]
++ 0.7 1 - (u′v′) Pr (u′v′)++
(36)
However, he suggested that 2.0 was a better coefficient than 0.7 for liquid metals (Pr , 1). It is informative to compare the predictions of eqs 32 and 34-36 functionally and numerically. Equations 32 and 35 predict a separate dependence of Prt on (u′v′)++ and Pr, whereas eq 36 predicts a dependence only on the combined variable Pr(u′v′)++/[1 - (u′v′)++]. This result is perhaps explicable in that Hill27 derived this combined dependence as an asymptote for Pr f 0, whereas eq 36 stretches the applicability of that asymptote for all values of Pr. The neglect of a dependence on (u′v′)++/[1 - (u′v′)++] by Jischa and Rieke can perhaps be interpreted as an approximation of this term by a mean value that results in only a slight error insofar as the term 0.85 is dominant. Equation 32 implies the existence of a value of Pr ) Prt ) PrT that is independent of (u′v′)++, thus justifying the reduction of eqs 19-21, 23, and 26-29 to expressions free of Prt and Pr for this condition. On the other hand, both eqs 35 and 36 predict a slightly varying dependence of Prt and PrT on (u′v′)++ for all values of Pr. Representative predictions of Prt by eqs 32 and 3436 are listed in Table 1. Values of d ) 8 in eq 33 and a coefficient of 0.0115 (rather than 0.015) in eq 34 are postulated for this comparison in the interest of consistency. The value of a+ ) 5000 was chosen to represent a moderate rate of flow and the value of y+ ) 500 to represent the outer range of the semilogarthmic regime
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002 6399 Table 2. Comparison of Predictions of PrT for a+ ) 5000 and y+ ) 500 Pr eq 34a (J&R) eqs 32 and 33 (Yahkot et al.)b a
0.0001 0.01867 0.01866
0.001 0.1839 0.1798
0.01 0.9681 0.7388
0.1 0.9223 0.8591
0.7 0.8653 0.8633
0.8633 0.8633 0.8633
1.0 0.8621 0.8634
10 0.8553 0.8635
100 0.8547 0.8635
With a coefficient at 0.0115. b Based on d ) 8.
Table 3. Comparative Dependence of Nu for a+ ) 5000 on Extreme Expressions for Prt and PrT Pr Nu1/Nu2 (Prt)1/(Prt)2 (PrT)2/(PrT)1 a
0.0001 1.000 7.962 1.0001
0.001 1.010 4.025 1.009
0.01 1.120 1.432 1.154
0.10 0.9913 1.003 1.003
0.70 1.0076 1.085 1.085
1.0 1.006 1.094 1.094
10 0.9844 1.119 1.119
100 0.9743 1.124 1.124
Prt based on eq 35 (N&S). b Prt based on eq 34 (J&R).
of the velocity distribution. The agreement for large and moderate values of Pr is very good except for the predictions of Notter and Sleicher, which are significantly higher. The predictions of Kays and Yahkot et al. are in fair agreement for small values of Pr, but those of Jischa and Rieke and of Notter and Sleicher are much higher. The appearance of an almost exact agreement between the predictions of Yahkot et al. and Kays in Figure 2 of the latter is to some extent an artifact of the compressed logarithmic coordinates. The predictions of PrT by eq 34 of Jischa and Rieke are compared in Table 2 with those of Yahkot et al. using d ) 8 in eq 33, for the same conditions as those for Table 1. It is noteworthy that the agreement is very close for Pr g 0.7 and for Pr e 0.001 despite the great structural difference of eqs 32 and 34. This is a somewhat fortuitous consequence, on the one hand, of the dominance of the limiting value of PrT ) 0.85 for even moderate values of Pr and, on the other hand, the approach of PrT to Pr/[1 - (u′v′)++] for Pr f 0 per eq 16 and thereby its independence from the expression used to predict Prt. These same two generalizations apply to the predictions (not shown) of PrT by eqs 35 and 36. Dependence of Nu on Prt The principal applicability of an expression for the prediction or correlation of values of Prt is in the prediction of values of Nu by means of eq 21 or 23 or their equivalent. It may be inferred from the reduced forms of eqs 21 and 23 for Pr ) 0 and Pr ) Prt and the limiting form for Pr f ∞ that the power dependence of Nu on Pr/Prt increases from zero to a maximum value approaching unity at a value of Pr/Prt ) 1 and then decreases to a value of 1/3 as Pr f ∞. The less than linear power dependence of Nu on Pr/Prt for all values of Pr thus fortuitously reduces the uncertainty in Nu arising from the uncertainty in Prt. This reduction is illustrated in Table 3, in which the ratio of the values of Nu calculated by Bo et al.7 for a+ ) 5000 and uniform heating using eqs 34 and 35 for Prt is compared with the corresponding ratios of values of Prt and PrT. The difference in the predicted values of Nu is very small except at Pr ) 0.01, for which it is 12% as compared to corresponding differences of 43% in Prt and 15% in PrT. The 2.6% difference in Nu for Pr ) 100 is a consequence of the 12.5% overprediction of the limiting value of Prt by eq 35, thereby a proof rather than an exception with respect to the insensitivity of Nu to Prt. Almost identical results are obtained for uniform wall temperature and for other values of a+. It should be emphasized that the comparisons in Table 3 were deliberately made using a simplistic expression (eq 34) and a more complex but questionable one (eq 35) for Prt in order to exaggerate the dependence of Nu on the expression used for Prt.
The actual uncertainty in the value of Nu resulting from the use of eq 32 or 36 for y+ > 30 and eq 34 for y+ < 30 would be expected to be far less, probably less than 1 or 2% for all values of Pr and a+. It may be concluded that improved numerical values, predictive equations, and/or correlative equations for Prt, while eagerly awaited in an intrinsic sense, will not greatly impact existing calculated values of Nu such as those of Bo et al.7 for round tubes, Danov et al.29 for parallel-plate channels, and Bo et al.30 for concentric circular annuli. Summary and Conclusions The eddy viscosity ratio µt/µ has been shown to be equal to the ratio of the local fraction of the shear stress rate of momentum transfer due to the turbulent fluctuations, namely, (u′v′)++ ) -Fu′v′/τ, to that due to the molecular motions, namely, 1 - (u′v′)++, and thereby independent of its heuristic diffusive origin. The eddy conductivity ratio kt/k has similarly been shown to be equal to the ratio of the local fraction of the radial heat flux density due to the turbulent fluctuations, namely, (T′v′)++ ) FcT′v′/j, to that due to the molecular motions, namely, 1 - (T′v′)++, and thereby also independent of the heuristic diffusive origin of kt. Despite their justification as physically meaningful quantities in round tubes and parallel-plate channels, the eddy viscosity and eddy conductivity are inferior in behavior and in implementation to (u′v′)++ or (u′v′)+ and T′v′)++ even in those two geometries. Turbulent flow and convection can be described completely and with clarity in terms of (u′v′)++ and (T′v′)++ while wholly avoiding the eddy viscosity, the eddy conductivity, and/or the corresponding mixing lengths. This approach might be feasible with undergraduate students, who have nothing to unlearn, were it not for the pervasive presence of the latter quantities in the literature of the past and present. On the other hand, the turbulent Prandtl number ratio, Prt/Pr, which was originally defined as the ratio of µt/µ to kt/k but which now may be reexpressed in terms of (u′v′)++ and (T′v′)++, as in eq 16, and interpreted as the ratio of the fractions of the local rate of momentum transport due to turbulent fluctuations and molecular motions divided by the corresponding ratio for energy transport, has a continued essential role in the prediction of convection, although its symbol and name may be vestiges of the past. Abbrecht and Churchill24 found experimentally that kt/k is a function only of µt/µ and Pr and thereby independent of the thermal boundary condition. The same dependence and independence obviously extends to the turbulent Prandtl number ratio Prt/Pr. The close
6400
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002
agreement of their results for kt/k and Prt/Pr with those of Page et al.25 for convection between parallel plates at different uniform temperatures confirms this finding and implies that Prt/Pr is independent of geometry as well when expressed in terms of µt/µ ) (u′v′)++/[1 (u′v′)++] and Pr. The only significant advance in the measurement and prediction of Prt since the recent review of Kays3 appears to be the Lagrangian DNS calculations of Papavassiliou and Hanratty20 for very large Schmidt numbers. The latter results contradict the postulate of a finite value of Prt at the wall that is the basis for the experimentally confirmed 1/3 power dependence of Nu on Pr. The explanation of this discrepancy is not yet resolved, but the practical consequences for heat transfer are unimportant because all normal fluids have values of Pr < 100. Predictions of Nu, which are the principal application of predictions for Prt, are fortuitously very insensitive to Prt. While improved measurements, predictions, and correlating equations for Prt are eagerly awaited, their impact on the currently predicted results for Nu can be expected to be only of second-order.
γ ) jτw/jwτ - 1 ) rate of dissipation of turbulent energy κ ) kinetic energy of turbulence µ ) dynamic viscosity µt ) eddy viscosity µT ) µ + µt ) total viscosity F ) specific density τ ) shear stress φ ) 90Pr3/2{(u′v′)++/[1 - (u′v′)++]}1/4 Subscripts m ) mixed-mean value w ) at wall
Acknowledgments The very discerning comments and constructive suggestions of the anonymous reviewers resulted in significant improvement of the manuscript and are gratefully acknowledged. The investigation that led to the unexpected findings reported herein might not have been undertaken without the prompting of Professor Tetsu Fujii following a seminar by the author at Kyushu University.
Nomenclature a ) radius of the round tube a+ ) a(τwF)1/2/µ b ) half-spacing of parallel plates b+ ) b(τwF)1/2/µ c ) heat capacity at constant pressure d ) “arbitrary” coefficient in eq 33 f ) 2τw/Fum2 ) Fanning friction factor j ) heat flux density in the y direction k ) thermal conductivity kt ) eddy conductivity kT ) k + kt ) total conductivity l ) mixing length l + ) l (τwF)1/2/µ Nu ) 2ajw/k(Tw - Tm) ) Nusselt number p ) pressure Pr ) cµ/k ) Prandtl number Prt ) cµt/kt ) Pr(u′v′)++[1 - (T′v′)++]/(T′v′)++[1 - (u′v′)++] ) turbulent Prandtl number PrT ) cµT/kT ) Pr[1 - (T′v′)++]/[1 - (u′v′)++] ) total Prandtl number R ) 1 - y+/a+ Re ) 2a+um+/µ ) Reynolds number Sc ) Schmidt number Sct ) turbulent Schmidt number T ) time-averaged temperature T+ ) k(τwF)1/2/(Tw - T)/µjw T′ ) fluctuating component of temperature T′v′ ) time-averaged value of T′v′ (T′v′)++ ) FcT′v′/j u ) time-averaged x component of velocity u+ ) u(F/τw)1/2 u′ ) fluctuating x component of velocity u′v′ ) time-averaged value of u′v′ (u′v′)+ ) -Fu′v′/τw (u′v′)++ ) -Fu′v′/τ v′ ) fluctuating y component of velocity x ) axial coordinate y ) distance from the wall y+ ) y(τwF)1/2/µ Greek Symbols R ) coefficient defined by eq 33
Literature Cited (1) Reynolds, O. On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion. Proc. R. Soc. London 1985, 186, 123. (2) Boussinesq, J. Essai sur la theorie des eaux courantes. (A Study of the Theory of Water Flows). Mem. Acad. Sci. Fr. 1877, 23, 1. (3) Kays, W. M. Turbulent Prandtl NumbersWhere are We? J. Heat Transfer, Trans. ASME 1994, 116, 284. (4) Churchill, S. W.; Chan, C. Turbulent Flow in Channels in Terms of the Turbulent Shear and Normal Stresses. AIChE J. 1995, 41, 2513. (5) Churchill, S. W. New Simplified Models and Formulations for Turbulent Flow and Convection. AIChE J. 1997, 43, 1125. (6) Reichardt, H. Die Grundlagen des Turbulent Wa¨rmeu¨bertraganges. Arch. Gesamte Wa¨ rmetech. (The Principles of Turbulent Heat Transfer). 1951, 2, 129. (7) Bo, Y.; Ozoe, H.; Churchill, S. W. The Characteristics of Fully Developed Turbulent Convection in a Round Tube. Chem. Eng. Sci. 2001, 56, 1781. (8) Weigand, B.; Ferguson, J. R.; Crawford, M. E. An Extended Kays and Crawford Turbulent Prandtl Number Model. Int. J. Heat Mass Transfer 1997, 40, 4192. (9) Prandtl, L. Eine Beziehung zwischen Wa¨rmeaustausch and Stro¨mungswiderstand der Flu¨ssigkeit. (An Analogy between Heat Transfer and Friction in Fluid Flow). Phys. Z. 1910, 11, 1072. (10) Launder, B. E.; Spalding, D. B. Mathematical Models of Turbulence; Academic Press: York, London, 1972. (11) Hanjalic, K.; Launder, B. E. A Reynolds Stress Model of Turbulence and Its Application to Thin Shear Flows. J. Fluid Mech. 1972, 52, 609. (12) Suzuki, K.; Tohkaku, A.; Sato, T. Liquid Metal Turbulent Heat Transfer in Concentric Annuli. Proceedings of the 8th International Heat Transfer Conference; Hemisphere Publ. Corp.: Washington, DC, 1980; Vol. III, p 969. (13) Schumann, U. Subgrid Scale Model for Finite Difference Simulation of Turbulent Channel Flows in Plane Channels and Annuli. J. Comput Phys. 1975, 18, 376. (14) Kawamura, H.; Nakamura, S.; Satake, S.; Kunagi, T. Large Eddy Simulation of Turbulent Heat Transfer in a Concentric Annulus. Thermal Science & Eng. 1994, 2, No. 2, 16. (15) Orszag, S. A.; Patterson, G. S. Numerical Simulation of Three-Dimensional Homogenous Isotropic Turbulence. Phys. Rev. Lett. 1972, 28, 75. (16) Kim, J.; Moin, P.; Moser, R. Turbulence Statistics in Fully Developed Turbulent Channel Flow at Low Reynolds Number. J. Fluid Mech. 1987, 177, 133.
Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002 6401 (17) Lyons, S. L.; Hanratty, T. J.; McLaughlin, J. B. LargeScale Computer Simulation of Fully Developed Turbulent Channel Flow with Heat Transfer. Int. J. Numer. Methods Fluids 1991, 13, 999. (18) Rutledge, J.; Sleicher, C. A. Direct Simulation of Turbulent Flow and Heat Transfer in a Channel. Part I. Smooth Walls. Int. J. Numer. Methods Fluids 1993, 15, 1051. (19) MacLeod, A. C. Liquid Turbulence in a Gas-Liquid Absorption System. Ph.D. Thesis, Carnegie Institute of Technology, Pittsburgh, PA, 1951. (20) Papavassiliou, D. V.; Hanratty, T. J. Transport of a Passive Scalar in a Turbulent Channel Flow. Int. J. Heat Mass Transfer 1977, 40, 1303. (21) Yahkot, V.; Orszag, S. A.; Yahkot, A. Heat Transfer in Turbulent Fields. 1. Pipe Flow. Int. J. Heat Mass Transfer 1987, 30, 15. (22) Elperin, T.; Kleeorin, N.; Rogachevskii, I. Isotropic and Anisotropic Spectra of Passive Scalar Fluctuations in Turbulent Fluid Flow. Phys. Rev. E 1996, 53, 3431. (23) Itazu, Y.; Nagano, Y. Modeling of Thermal Eddy Diffusivity Based on Renormalization Theory (In Japanese). Trans. Jpn. Soc. Mech. Eng. 1997, B36, 3072. (24) Abbrecht, P. H.; Churchill, S. W. The Thermal Entrance Region in Fully Developed Turbulent Flow. AIChE J. 1960, 6, 268.
(25) Page, F., Jr.; Schlinger, W. G.; Beaux, D. K.; Sage, B. W. Point Values of Eddy Conductivity and Viscosity in Uniform Flow Between Parallel Plates. Ind. Eng. Chem. 1982, 44, 424. (26) Jischa, M.; Rieke, H. B. About the Prediction of Turbulent Prandtl and Schmidt Numbers from Modeled Transport Equation. Int. J. Heat Mass Transfer 1979, 22, 1547. (27) Hill, J. W. Estimation of Eddy Diffusivities in Isotropic Turbulence. Preprints of the 6th Annual AIChE Meeting, San Francisco, 1972. (28) Notter, R. D.; Sleicher, C. A. A Solution to the Turbulent Graetz ProblemsIII. Fully Developed and Entry Region Heat Transfer Rates. Int. J. Heat Mass Transfer 1972, 27, 2073. (29) Danov, S. N.; Arai, N.; Churchill, S. W. Exact Formulations and Nearly Exact Numerical Solutions for Convection in Turbulent Flow Between Parallel Plates. Int. J. Heat Mass Transfer 2000, 43, 2767. (30) Bo, Y.; Ozoe, H.; Churchill, S. W. Nearly Exact Numerical Solutions for Turbulent Convection in an Internally Heated Circular Concentric Annulus. In review.
Received for review December 17, 2001 Revised manuscript received May 29, 2002 Accepted June 6, 2002 IE011021K
