2816
Ind. Eng. Chem. Res. 2004, 43, 28162829
CFD Simulation of Heat Transfer in Turbulent Pipe Flow Mahesh T. Dhotre and Jyeshtharaj B. Joshi* Institute of Chemical Technology, University of Mumbai, Matunga, Mumbai400 019, India
In the present work, the flow pattern in pipe flow has been simulated using a lowReynoldsnumber k model. The model has been extended to predict the heattransfer coefficient in regions of both high and low Prandtl number. The computational fluid dynamic (CFD) approach has been shown to be superior to the conventional semiempirical approach for high turbulent Prandtl numbers and in cases where the molecular properties such as diffusivity and viscosity steeply change in the vicinity of the wall. 1. Introduction It is wellknown that the rates of transfer of heat, mass, and momentum between a solid wall and a fluid in turbulent flow can, under suitable conditions, be characterized by an eddy diffusion coefficient. The concept of an eddy diffusion coefficient, however, does not provide a physical basis to illuminate the nature of turbulence; it does provide a means of calculating massor heattransfer rates from momentumtransfer rates. The calculation procedure involves determining the momentum eddy diffusivity distribution from the slope of the velocity profiles, taking the eddy diffusivity of mass and heat to be equal or proportional to the eddy diffusivity of momentum, and then integrating the transport equations to obtain a temperature or concentration profile. The real success of any such analysis depends to a considerable extent on how carefully the eddy diffusivity variation is chosen near the wall. Calculations of this type have been reported in the literature, and the principal difference among most of them is the starting point of the velocity distribution. In view of the importance of the accurate description of the velocity profile and the eddy diffusivity near the wall, an attempt is made in the present work to employ computational fluid dynamics (CFD) for this purpose. In recent years, lowReynoldsnumber k models of turbulence have been widely used in numerical simulations because of their ability to resolve the nearwall region. The lowReynoldsnumber k modeling approach incorporates either a walldamping effect, or a direct effect of molecular viscosity, or both in the turbulence transport equations. Fairly complete reviews of lowReynoldsnumber k models of turbulent shear flows have been presented by Patel et al.,1 Hrenya et al.,2 and Thakre and Joshi.3,4 Nearwall turbulence models or lowReynoldsnumber models that attempt to describe the relative influence of molecular and turbulent viscosities have been developed for singlephase flows. Thakre and Joshi3 have analyzed 12 different lowReynoldsnumber k models for the case of a singlephase pipe flow. For this purpose, they used four criteria, three of which were accurate predictions of the radial variation of axial velocity, the turbulent kinetic energy, and the eddy diffusivity (compared with nearwall experimental data of Durst et al.5). The fourth * To whom correspondence should be addressed. Tel.: 0091222414 5616. Fax: 0091222414 5614. Email:
[email protected] udct.org.
criterion was rather stringent and stipulated the necessary condition of the overall energy balance, i.e., the volume integral of must be equal to the energyinput rate (the pressure drop multiplied by the volumetric flow rate). All four criteria were found to be satisfied by the model of Lai and So6 (LSO), as shown by Thakre and Joshi.3 The authors successfully applied Lai and So6 model to predict the heattransfer coefficient in singlephase pipe flows. Because the heattransfer process depends entirely on the flow pattern, a model that gives good flow predictions can be expected to give good heattransfer predictions as well. However, this statement should be in the form of a question. There are a few possible reasons for this limitation: (1) Few attempts have been made to extend the flow knowledge from lowReynoldsnumber models to heat transfer. (2) The accuracy of modeling the energydissipation equation (which directly affects the flow quantities and thereby heat transfer)is limited. (3) The selection of the turbulent Prandtl number, Prt (one value or radial profile), is difficult. The Prt concept is widely used in an eddy diffusivity model, in which the eddy diffusivity for the momentum is evaluated by either a mixinglength model or the k model and the eddy diffusivity for heat is evaluated using Prtbased relationships. The knowledge of Prt is a central problem of all theoretical considerations concerning turbulent heat transfer in the duct flows. A large number of models that address the estimation of Prt for such situations have been published. Reviews of the existing work can be found in Reynolds,7 Kays,8 and more recently Churchill.9 In view of the above considerations, it was thought desirable to analyze the reported turbulent Prandtl number relationships and understand turbulent heat transfer using a lowReynoldsnumber k model over a wide range of Prandtl and Reynolds numbers. Toward this end, we selected the Prt formulation proposed by Yakhot and Orszag10 and compared the results with the experimental data of Monin and Yaglom,11 Skupinski et al.,12 Fuchs,13 and Kader.14 2. Previous Work After the pioneering beginning by Reynolds,15 a large number of analytical solutions have been published until fairly recently.16 In all of these cases, correlations were given for eddy diffusivity together with boundary conditions. The mass or heat balance equation was solved analytically, and the resulting predictions for the mass and heattransfer coefficients were given in the
10.1021/ie0342311 CCC: $27.50 © 2004 American Chemical Society Published on Web 04/24/2004
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2817
form of correlations. Asymptotic results have also been presented for high values of Sc and Pr. These types of studies spanned over 125 years; a brief summary of them, in chronological order, is presented in Table 1, and further details are provided by Thakre and Joshi.4 Churchill and coworkers9,2830 investigated this problem comprehensively and presented new formulation for the prediction of turbulent flow and convection. Churchill28 proposed interesting simplified models and integral solutions for the case of fully developed convection in terms of the local fraction of transport due to turbulence. A very general and accurate correlating equation for the turbulent shear stress allowed the use of these integral formulations in the derivation of even more accurate numerical solutions for the velocity distribution in fully developed flow. This implementation for numerical calculations necessarily invokes some empiricism, but generally less than that required by prior models and formulations. Heng et al.29 used models proposed by Churchill28 to compute improved values for the Nusselt number for fully developed convection in uniformly heated round tubes over a wide range of Reynolds numbers. The results for two limiting (Pr ) 0, Pr f 0) and one intermediate (Pr ≈ 0.87) values were presented and found to be more accurate than the previous predictions. Yu et al.30 used models proposed by Churchill28 and presented differential models for fully developed turbulent flow and convection. They considered two cases, one for uniform heating and another more complex case of uniform temperature. They solved the differential model for all values of Pr and a wide range of Reynolds numbers. They found that the predictions (for the case of uniform heating) are slightly more accurate than those of Heng et al.29 Even though the latter authors had used the same basic model, Yu et al.30 obtained rapid and complete convergence by incorporating a stepwise procedure for solving the differential equations, whereas Heng et al.29 had evaluated integrals by quadrature. The predicted rates of heat transfer were found to be relatively insensitive to the particular empirical expressions used for the turbulent Prandtl number, which is the only significant source of uncertainty. Churchill9 showed the eddy viscosity at any location in fully developed turbulent pipe flow to be simply a ratio of the shear stress due to the timeaveraged turbulent fluctuations in the velocity to that due to the molecular motion. The eddy conductivity in fully developed turbulent convection was similarly shown to be a ratio of the corresponding contributions to the radial heat flux density. Churchill9 reexpressed the turbulent Prandtl number in terms of stresses and used different Prt expressions proposed by Notter and Sleicher,25 Kays,8 Jischa and Rieke,31 and Yakhot et al.10 He found that the predictions of Kays8 and Yakhot et al.10 were in fair agreement for small values of Pr, but those of Jischa and Rieke31 and Notter and Sleicher25 were much higher. Further, he found the prediction of the Nusselt number to be insensitive to the turbulent Prandtl number and pointed out a need for improved measurements of flow parameters and correlating equations for turbulent Prandtl number. In a totally different approach, Fortuin et al.16 extended the random surface renewal (RSR) model of Fortuin and Klijn32 to describe heat and masstransfer processes in turbulent pipe flow assuming that the tube
wall is covered by a mosaic of fluid elements of laminar flow and random ages. This socalled extended random surface renewal model (ERSR) was used to derive a relationship between the friction factor and the mean value of the randomly distributed ages of the fluid elements at the tube wall, showing that the age distribution and the mean age of the fluid elements at the tube wall agree quantitatively with the experimental results obtained from the velocity signals measured with a laser Doppler anemometer. This information was then used to calculate the timeaveraged radial profiles of the axial velocity, the temperature, and the concentration in the wall region; to derive equations for the quantitative prediction of heat and masstransfer coefficients; and to provide a basis for explaining the analogy between momentum, heat, and mass transfer in turbulent pipe flow. Figure 1 shows a comparison of predictions of various theories. It can be seen that large variations exist in the predictions, both among themselves and as compared with the experimental data. The above discussion (also refer to Thakre and Joshi4) gives an account for these variations. They are mainly due to the particular selection of eddy diffusivity profiles and the other simplifications in the solution procedure. In view of the importance of accurate descriptions of the velocity profile and the eddy diffusivity near the wall, various attempts have been made to employ computational fluid dynamics (CFD) for the prediction of heat transfer. These efforts can be classified as (a) eddy viscosity models, (b) Reynoldsstress model (RSM), (c) large eddy simulation (LES), and (d) direct numerical simulation (DNS). The widely used eddy viscosity models are of three types: k (Launder and Spalding33), kτ (Speziale et al.34), and kω (Wilcox35). Sarkar and So36 carried out a critical evaluation of these models against direct numerical simulation for six flow cases. They found that the k models that are asymptotically consistent perform better than the kτ and kω models. Moreover, the kτ and kω models exhibit no real numerical advantage over the k models. There are two main methods for the nearwall treatment of these models: (i) the wall function approach (Launder and Spalding37) and (ii) lowReynoldsnumber modeling. In most industrial CFD simulations, a wall function approach is used, which bridges the viscous sublayer and makes use of the available knowledge in the logarithmic profile region of the boundary layer. This method allows for the use of much coarser nearwall grids, which also benefits the cell aspect ratio. On the downside, wall functions impose strict limitations on the gridgeneration process, as they impose an upper limit on the grid density so that it remains in the logarithmic part of the boundary layer. Experience has shown that this is a severe limitation that has negatively impacted a large number of industrial CFD simulations. Another more promising approach is lowReynoldsnumber modeling. Jones and Launder38 were the first to propose a lowReynoldsnumber k model for nearwall turbulence, which was then followed by a number of similar k models. Fairly complete reviews of the development of such twoequation nearwall turbulence closures have been provided by Patel et al.,1 Hyrenya et al.,2 and Thakre and Joshi.4 Thakre and Joshi4 presented a systematic analysis of 12 versions of lowReynoldsnumber k models for heat transfer. The model by Lai and So6 was shown
momentum and energy transported at equal mass rates between the bulk of the fluid and the wall of the tube by the oscillatory radial motion of the eddies eddies penetrate only to a finite distance from the wall, and transport of momentum and energy through the remaining distance occurs by a linear process of molecular diffusion transport in the buffer layer assumed to occur by both molecular and eddy diffusivity; universal velocity profile used for the estimation of νT; closed form obtained by neglecting molecular diffusion in the turbulent core as in von Karman analysis, transport in the buffer layer assumed to occur by both molecular and eddy diffusivity; universal velocity profile used for the estimation of νT following Prandtl, assumed that νT ) 0 in the viscous sublayer and turbulent effects predominate in finite distance from the wall; distance can be estimated by the point at Prt ) 1 doubted truth of existence of a sublaminar layer for molecular diffusion and attempted to model by introducing into the laminar sublayer a small appropriate amount of eddy motion that decreases very rapidly to zero at the wall; implemented by introducing a correlation for νT in the laminar sublayer (y+ < 5) in terms of eddy diffusivity varying as the cubic power of the normal wall distance (a) eddy diffusivities for momentum (νT) and heat transfer (νH) assumed equal; variations of the shear stress (τ) and heattransfer flux (q) across the tube assumed to have a negligible effect on the velocity and temperature distributions (b) contribution of molecular shear stress and heat transfer can be neglected in the region sufficiently far from the wall
Reynolds15
Deissler22
Lin et al.21
Levich20
Martinelli19
von Karman18
Prandtl17
postulate
researcher
Table 1. Previous Work: Analytical Approach
) )
1  y+ R+ 1  y+ R+
( (
+ ) 1 + q[1  exp(q)] y + ) 2.78 q ) n2u+y+
( )
+ ) 1 + ) 0.4y+ y+ 1 can be estimated by the point at Prt ) 1 y+ 3 + ) 1 + 14.5 + ) 1 + 0.2y+  0.959 + + ) 0.4y
+ ) 0.4y+
+ ) 0.2y+
0 e y+ e 26 y+ > 26
0 < y+ e 5 5 < y+ < 30 y+ > 30
0 e y+ e y+ 1 y+ 1 > 30
0 e y+ e 5 5 e y+ e30 y+ > 30
0 e y+ e 5 5 e y+ e 30 y+ > 30
+ ) 1 + ) 0.2y+ + ) 0.4y+

range
0 e y+ e 11.5 y+ > 11.5
)
νt ν
+ ) 1 + ) 0.4y+

(
+ ) 1 +
(
) )
0.4
Nu )
x2f Re Pr
)
)
0.5
[ ( ) ]}
ln Re + ln Pr + ln x2 / f + (1 + ln 0.2)
(
0.248Re Pr xf(Pr)0.75 π
f Re Pr x 2 Nu ) x2f + F(Pr)
Nu )
{
x1f (T
(1/Prt)
(
1 + 5Pr (2/f)0.5 + 5(Pr  1) + 5 ln 6
Re Prxf/2
1 1 δ+ 1 + 1Pr Re(f/2) Re Pr(f/2)
TW  T C Re Pr W  TM Nu ) Pr Re f Pr 5 + ln 1 + 5 + 2.5F ln Prt Prt 60 2
Nu )
Nu )
Nu ) Re Pr(f/2)
resulting equations for mass and heattransfer coefficients
used powerlaw expression for νT/ν near the wall
equation found to agree well with data over a wide range of Pr and Sc values (0.53200)
analysis found to be in agreement with experimental data for heat transfer in liquid metals
applicable for Pr ) 0.530
improvement over Prandtl treatment for Pr > 1 due to allowance made for the existence of turbulence down to y+ ) 5, rather than y+ ) 8.7, as by Prandtl
0.73 < Pr < 40 improvement over Reynolds treatment for Pr > 1 due to allowance made for the existence of turbulence down to y+ ) 8.7
applicable only for Pr ) 1
remarks
2818 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004
assumed shear stress and heat flux to be linear functions of the radial position within the tube; assumed eddy diffusivities for heat transfer and momentum transfer to be equal; maintained the form of the Prandtl analogy; values of constants varied to account for the turbulent transport in the laminar sublayer and the variation in thickness of the laminar sublayer with respect to the Reynolds number eddy viscosity distribution obtained in terms of increasing powers of y+ from the wall; away from the wall, von Karman’s velocity distribution function retained; after velocity and νT distributions defined, diffusionconvection equation solved numerically to obtain expressions for local masstransfer rates without reference to any hypothetical mechanism of turbulence, chose a simple, empirical function for the eddy diffusivity that is constrained to go to zero as y f 0 and to agree with measurements in the region where measurements are reliable, y+ > 10 logarithmic velocity distribution law does not hold for the region close to the wall (i.e., y+ < 30); proposed a model assuming the eddy viscosity to be proportional to the distance from the wall, shear stress to be constant, and eddy viscosity and velocity to be continuous for the whole range of pipe section presented a method to calculate heattransfer coefficients from a more fundamental approach using eddy viscosity data in the integration of the energy equation
Friend and Metzner23
Edwards and Smith27
Muzushina and Ogino26
Notter and Sleicher25
Wasan and Wilke24
postulate
researcher
Table 1. (Continued)
+ 3
+ 4
a(y+)3  b(y+)4
)
νt ν
)
St ) 0.827b2/3(f/2)0.5Pr(2/3)
( 
valid for Pr ) 0.6100
+ ) 5.23 × 104(y+)3
valid for Pr > 100
valid for Re ) 3000(3 × 106) Pr ) 0.5600
remarks
St ) 0.0149Re0.12 Pr(2/3) Nu ) 5 + 0.016ReaPrb a ) 0.88  0.24/(4 + Pr) b ) 0.33 + 0.5 exp(0.6Pr)

0 < y+ < 45
Re Pr(f/2)
()
Re Pr(f/2) f 1/2(Pr  1) 1.2 + 11.8 2 Pr1/2
, 1 + xf/2[13.0(Pr)0.8  13.0] 0.2 e Pr e 2 Re Pr(f/2) Nu ) , 1 + xf/2[13.8(Pr)0.71  13.0] 2.0 e Pr e 100 Nu ) 0.058xf/2Re Pr0.34, 100 < Pr < 10000
Nu )
Nu )
resulting equations for mass and heattransfer coefficients
0 e y+ e y+ 1 + + y+ 1 e y e y2 + y2 e y+ e R+
0.00090y+3 (1 + 0.0067y+2)0.5
0 e y+ e 20

range
+ ) by+3 y+ y+ + ) 1  + 1 2.5 R + ) 0.07R+ R ) tube radius
+ )
1  a(y )  b(y ) a ) 4.16 × 104 b ) 15.15 × 106
+ )

(
+ ) 1 +
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2819
2820 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004
Figure 1. Comparison of the predictions of Nusselt number for Re ) 10 000 using different approaches: (1) Reynolds,15 (2) Prandtl,17 (3) von Karman,18 (4) Colburn,89 (5) Deissler,22 (6) Wasan and Wilke,24 (7) Notter and Sleicher,25 (8) Friend and Metzner,23 (9) Edward and Smith,27 (10) Kays.8
to have the best overall predictive ability for heat transfer in turbulent pipe flow. The success of the Lai and So model was attributed to its success in predicting flow parameters such as mean axial velocity, turbulent kinetic energy, and energy dissipation rate. Further, they recommended using Prt variation near the wall to improve the predictive ability of this model. “Reynoldsstress” or “secondmoment” turbulence models do not rely on the turbulent viscosity concept, but instead incorporate transport equations for secondorder velocity and temperature correlations. Launder39 provides a comprehensive discussion of these models. The complexity of Reynoldsstress models led to the development of truncated forms known as “algebraic stress models” (ASMs). Despite some successes, it is still believed that the most reliable prediction methods are those based on a secondmoment closure. The reason is that the turbulent interactions that generate the Reynolds stresses and heat fluxes can be treated with less empiricism. Moreover, for those processes that cannot be so handled, a more rational and systematic set of approximations can be derived. Although much progress has been achieved in recent years in the modeling of the Reynoldsstress transport equations, the modeling of the scalar field, on the other hand, is still rather primitive. This is because the development of the model for the heat flux depends largely on the availability and correctness of turbulent stresses predicted by the Reynoldsstress model. Turbulent heattransfer process models using highReynoldsnumber secondmoment closures have been developed by meteorological fluid dynamicists.4042 On the other hand, applications of similar secondmoment turbulence closures to engineering heattransfer problems have been attempted by Baughn et al.43 and Launder and Samaraweera,44 among others. Prud’homme and Elghobashi45 proposed a secondorder closure for momentum and an algebraic stress model for heat transfer. In all such secondmoment closure models, the principal hurdle has been the shortage of reliable and relatively accurate nearwall heat flux measurements, which has contributed to the slow development of the subject of nearwall turbulence modeling for heat fluxes. Recently, the model was extended by Yoo and So46 to calculate isothermal, variabledensity flows in a suddenexpansion pipe. In their approach, the flow and turbu
lence fields were resolved by a lowReynoldsnumber secondmoment closure. The scalar flux equation was closed by highReynoldsnumber models, and the nearwall scalar fluxes were evaluated assuming a constant turbulent Schmidt number. This is one way to handle the scalar fluxtransport equations for the nearwall flow, even though the approach is known to be quite inappropriate for most turbulent heat and masstransfer problems of engineering interest (Launder,47 Yoo and So46). The reason for this appears to be that, so far, no suitable nearwall secondmoment closure for scalar flux transport has been developed. This is due, in part, to a lack of detailed nearwall scalar flux measurements and, also in part, to the unavailability of an asymptotically correct nearwall Reynoldsstress model. Lai and So48 developed a nearwall Reynoldsstress turbulence model that can correctly predict the anisotropy of the turbulent normal stresses. The success of that model motivated Lai and So6 to extend the momentum approach to model turbulent heat transport near a wall. The modeling approach was similar to that outlined in Lai and So48 and was based on the limiting wall behavior of the heat fluxtransport equations. In this way, the modeled equation was valid all the way to the wall and the assumptions of a temperature wall function and a constant turbulent Prandtl number were not required. It should be noted, however, that detailed and accurate experimental documentation of wall turbulent flow is presently not available and the modeling of the dissipation rate of temperature variance is quite immature, as even the highReynoldsnumber version of the modeled equation is not well developed. Attempts to evaluate the various models have been published in the literature. In one such attempt, Martinuzzi and Pollard49 compared six different turbulence models (two lowReynoldsnumber k models and four algebraic stress models) and found that, for fully developed pipe flow, the lowReynoldsnumber k model provided the best predictions. In a similar attempt Pollard and Martinuzzi50 compared seven different turbulence models (twolowReynoldsnumber models and five Reynoldsstress models) and found that the lowReynoldsnumber k model provided predictions that were generally in closer accord with various experimental data. Thakre and Joshi3,51 found that the Lai and So k model and the Reynoldsstress model gave overall good agreement with the experimental data for the cases of three Prandtl numbers. A comparative analysis between the k and Reynoldsstress models for the prediction of the mean temperature and the Nusselt number showed the applicability of k model for heat transfer to be relatively superior. Even though the Reynoldsstress models neglect the assumptions of isotropy and the turbulent Prandtl number, their applicability for heat transfer needs reconsideration. This is mainly because of the speculation that the heat flux transport is more complicated than momentum transport and is more likely to be influenced by two or more time scales instead of one. A second reason might inappropriate nearwall modeling of the dissipation and pressure scrambling terms. The overall discrepancy observed in both the k and Reynoldsstress models for heat transfer can be attributed to our limited understanding of the dissipation term and also the lack of detailed nearwall temperature and scalar flux measurements at higher Prandtl numbers.
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2821
The first application of a large eddy simulation (LES) was made by Deardorff,52 who simulated turbulent channel flow at an indefinitely large Reynolds number. In this pioneering work, Deardorff showed that threedimensional computation of turbulence (at least for simple flows) is feasible. He did not fully resolve the turbulence but used a subgridscale (SGS) model to describe the behavior of the smallest scales. This use of subgridscale models made it possible to simulate turbulent channel flow at very high Reynolds number, despite the low resolution. Using only 6720 grid points, Deardorff was able to predict several features of turbulent channel flow with a fair amount of success. Of particular significance was the demonstration of the potential of LES for use in basic studies of turbulence. Following Deardorff’s work, Schumann53,54 also calculated turbulent channel flow and extended the method to cylindrical geometries (annuli). He developed a finitedifference method for simulations of turbulent flows in plane channels and annular pipe flows. His results agreed rather well with the experimental values. It should be noted that Schumann used up to 10 times more grid points than Deardorff, as well as an improved subgridscale model. In addition to dividing SGS stresses into a locally isotropic part and an inhomogeneous part, he employed a separate partial differential equation for the SGS turbulent kinetic energy. However, the added differential equation did not improve the results over the calculations in which only an eddy viscosity model was used (Schumann54). Grotzbach and Schumann55 extended their channel flow calculations to account for temperature fluctuations and heat transfer. Later extensions by Grotzbach included calculations of secondary flows in partly roughened channels, inclusion of buoyancy effects, and liquid metal flows in plane channels and annuli. Moin and Kim56 used LES for the simulation of channel flow at a Reynolds number of 13 800, but unlike Schumann, they did not model the walllayer dynamics. Instead, they extended the calculations down to the wall. They resolved the wall layer in their LES of channel flow. Their results were found to be in agreement with experiments and have been used as input data for numerous studies of organized structure in wallbounded turbulent flows. The arrival of direct numerical simulation (DNS) has shed some light on detailed flow structures in the nearwall region. Further, the detailed energy budgets derived from DNS data provide a route toward the modeling of wall turbulence. Orszag and Kells57 performed the first direct numerical simulation of channel flow. They performed a study of the nonlinear stages of the transition from laminar to turbulent channel flow. Orszag and Patera58 simulated turbulent flow in a channel by using the Orszag and Kells57 algorithm with the fully explicit AdamsBashforth scheme to time advance the nonlinear terms. The Orszag and Kells57 algorithm was found to be satisfactory as long as it was used for highReynoldsnumber flows, for which the behavior at the wall is not critical. Marcus59 developed a modification of the Orszag and Kells57 algorithm that incorporates viscous effects into the pressure model and satisfies continuity at the walls. Marcus59 used his algorithm to simulate TaylorCouette flow (where, for this problem, the behavior at the wall is critical to the entire flow field) and obtained excellent agreement with experimental measurements of wavy Taylor vortex flow.
Kim, Moin, and Moser60 applied DNS to investigate fully developed turbulent flow between two parallel plates. At a Reynolds number of 3300 using nearly 4 × 106 grid points, they computed various statistical parameters for the flow field. The general characteristics of the turbulence statistics showed good agreement with the experimental results of Ecjelmann61 and Kreplin and Eckelmann,62 except for some flow quantities in the wall layer. Recently, using same code, a similar DNS was performed at a larger Reynolds number of 7900. Together with the results of the previous DNS of Kim, Moin, and Moser,60 Antonia et al.63 focused on lowReynoldsnumber effects in the inner region of the turbulent flow. The DNS and experimental results reported in their paper agreed well and indeed showed significant lowReynoldsnumber effects. In addition to these two computations, several other direct simulations of wallbounded turbulent flows have been reported in the literature. For instance, some of these include studies of plane channel flow by Lyons et al.,64 developing turbulent boundary layer on a flat plate by Spalart,65 turbulent flow in a rotating channel by Kristoffersen and Andersson,66 and turbulent flow through a square duct by Gavrilakis.67 All of these computations have one commonality in that they considered geometries with either a rectangular cross section or a flat plate. Unger et al.68 used a secondorderaccurate finitedifference method for the DNS of fully developed turbulent pipe flow, with Reynolds numbers up to 5300. They obtained good agreement with experiments by Westerweel et al.69 for the secondorder statistics and fair agreement for higherorder statistics. Zhang et al.70 used a spectral method in their DNS of pipe flow with Reynolds numbers up to 4000. Their code uses Fourier transforms in the axial and azimuthal directions and spectral elements with Jacobi polynomials in the radial direction. Eggels et al.71 carried out a DNS to study fully developed turbulent pipe flow at a Reynolds number Re ) 7000. They found that the mean velocity profile in the pipe fails to conform to the accepted law of the wall, in contrast to the channel flow. This result confirms earlier observations reported in the literature. Many attempts have been made to apply DNS for fully developed turbulent channel flow because of its simple geometry and fundamental nature to understand the convective heat transfer between fluid and a solid wall. The first attempt in this area was made by Kim and Moin72 for Pr ) 0.1, 0.71, and 2.0 with Ret ) 180. Later, Lyons et al.64 carried out a similar DNS for Pr ) 1.0 with a lower Reynolds number of Re* ) 150. The two walls of their channel were kept at different temperatures. Kasagi et al.73 and Kasagi and Ohtsubo74 carried out DNS for Pr ) 0.71 and 0.025 with a slightly lower Reynolds number of Re* ) 150. They obtained the aggregate of the transport equations for the temperature variance, turbulent heat fluxes, and their dissipations. Papavassiliou and Hanratty75 carried out a DNS of channel flow for the velocity field and used a Lagrangian method to describe turbulent transport of the scalar field. They used tracking of heat markers that are released instantaneously from a point source on a wall boundary in a numerically simulated turbulent channel flow. This methodology was used to describe the transfer of heat. The calculated fully developed temperature profile was compared with a direct numerical simulation. Further, they pointed out that DNS solutions of
2822 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 Table 2. Previous Work: CFD Approach researcher(s)
geometry
Reynolds number(s)
grid(s)
method
Jones and Launder38
pipe

5 × 105
lowRe k
Martinuzzi and Pollard49 Pollard and Martinuzzi50 Lai and So6
pipe
∼120 × 50
10 000380 000
RSM, k
pipe
∼120 × 50
10 000380 000
ASM, k
pipe
nonuniform
21 00071 000
RSM, lowRe k
Wilcox35
pipe
∼100
26005300
k, kω
Sarkar and So36
channel and Couette
nonuniform, 110
1300 and 3000
k, kτ, kω
Thakre and Joshi3,4,51
pipe
nonuniform, ∼120
700050 000
k, RSM
Deardorff52
channel
6760

LES
Schumann54
channel
65 536

LES
Grotzbach and Schumann55 Moin and Kim56
channel
65 576
250 000
LES
channel
516 096
13 800
LES
Kim et al.60
channel
3 962 880
3300
DNS
Rai and Moin
channel
1 392 640

DNS
Lyons et al.64
channel
1 064 960
2262
DNS
Eggels71
pipe
3 145 728
7000
LES, DNS
Zhang et al.70
pipe
4000
DNS
Papavassiliou and Hanratty75
channel
1 064 960
2660
Orlandi and Fatica90
rotating pipe
nonuniform, 786 432
7000
DNS, Lagrangian method for heat transfer DNS
Eulerian heat balance equations are not possible at high Pr given the present status of computers and that the use of a Lagrangian method can overcome this difficulty. Kawamura et al.76 performed DNS for a wider range of Prandtl numbers from Pr ) 0.025 to Pr ) 5.0 with Re* ) 180. Recently, Wikstrom and Johansson77 performed a DNS with the higher Reynolds number value of Re*
remarks first to apply lowReynoldsnumber approach; suggested damping function and constants compared five turbulence models; k model found to give good predictions compared seven turbulence models; k model found to give good predictions asymptotically correct nearwall Reynoldsstress closure and k closure proposed; also verified the notion that Prt is not constant near a wall; if Prt is assumed to be constant, the results obtained are at variance with measurements compared eight lowRe k and kω models for highReynoldsnumber, incompressible, turbulent boundary layers with favorable, zero, and adverse pressure gradients; results obtained underscore the k model’s unsuitability for such flows compared 10 nearwall models; found k models that are asymptotically consistent perform better than kτ and kω models compared 12 k models and 2 RSM, applied energy balance criterion, k model found to give good predictions first application of LES to turbulent bounded flows with coarse grid used and improved the subgridscale model; employed a separate partial differential equation for SGS turbulent kinetic energy, but this did not improve results over the calculations used SGS model and extended for SGS temperature transport resolved nearwall region instead of using earlier artificial boundary conditions to account for effect of nearwall dynamics resolved all essential turbulence scales on the computational grid and used no subgrid model compared influence of a finitedifference vs spectral approach on the statistical results used pseudospectral technique to solve the equation without use of the subgridscale modeling and extended calculations to heat transfer found that the mean velocity profile in the pipe fails to conform to the accepted law of the wall, in contrast to channel flow found that the mean velocity profile in the pipe fails to conform to the accepted law of the wall, in contrast to channel flow; confirms earlier observations reported in the literature found that the use of Lagrangian methods can overcome the difficulty of solving the heat balance equation for high Pr by only use of DNS simulations performed by a finitedifference scheme, secondorder accurate in space and time; when pipe rotates, a degree of drag reduction is achieved in the simulation just as in the experiments
) 265 with the different thermal boundary conditions of a uniform temperature difference. The details of all of these methods are provided in Table 2. The results of such direct numerical simulations, which appear to be convergent and in good agreement with experimental measurements, are yet limited to two fully developed flows, namely, flow between parallel
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2823
plates and flow along a flat plate, and, in both instances, to a fairly low range of fully turbulent motion. Also, the results of DNS are in the form of discrete numerical values and therefore provide no direct functional insight. Because of these limitations and because of the truly great computational requirements, this methodology is not yet of direct utility for design and control. 3. Mathematical Model For steady, isothermal, incompressible, fully developed turbulent pipe flow, the set of governing equations for the flow and heat are given as follows:
Momentum equation 0)
1 d du 1 dp r(ν + νT) r dr dr F dz c
[
] ( )
(1)
Transport equation for the turbulent kinetic energy (k) 0)
[ ( ) ] ()
νT dk du 1 d r ν+  νT r dr σk dr dr
2
D
(2)
Transport equation for the turbulent energy dissipation rate () 0)
{ [( ) ]}
νT d 1 d r ν+ r dr σ dr
+
( )
C1f1νT du k dr
2

C2f22 + E (3) k
where
k2 νT ) Cµfµ
(4)
fµ ) 1  exp(0.0115y+)
(5)
[
( )] [ ( ) ] [ ( )] [ ( ) ]{( ) ( )] ( )}
f1 ) 1 + 1  0.6 exp f2 ) 1 
( ) [
dxk E ) 2νC2f2 k dr
2
Re 104
exp 
Ret 2 exp 9 64
Ret 64
2
(6)
2
(7)
Ret 2 7 C 2 × 64 9 2 k 2 1 2νν 2 dxk  2ν (8) dr 2k y + exp 
All of the lowRe models adopt a damping function, fµ, to account for the direct effect of molecular viscosity on the shear stress near the wall (viscous sublayer and buffer zone). It should be noted that the wall functions used in connection with the standard k model (fµ ) 1) are usually applied in a region y+ > 30. Further, in the lowReynoldsnumber k model, the function f2 is introduced primarily to incorporate the lowReynoldsnumber effect in the destruction term of . The important criterion for the function f2 is that it should force the dissipation term in the equation to vanish at the wall. For highReynoldsnumber flows far from the wall, the function f1 asymptotically approaches the value of 1 in accordance with the highReynoldsnumber form of the model. The k model parameters are Cµ ) 0.09, C1 ) 1.35, C2 ) 1.8, σk ) 1, and σ ) 1.3.
The following thermalenergy equation defines the heattransfer problem
u
[ (
) ]
νT/ν dT dT 1 d 1 ) rν + dx r dr Pr Prt dr
(9)
We consider heat transfer in a pipe with constant heat flux through the wall. Equation 9 was normalized, resulting in the following form: r+ ) r/R, u+ ) u/uτ, T + ) T/T*, where T* ) q/CpFuτ and ν+ ) ν/uτR. Before an attempt was made to analyze the available data at high Prandtl number, it was thought desirable to examine the behavior of the momentum eddy diffusivity at points very close to the wall, i.e., at values of y+ < 5. At high Prandtl number, this is precisely the region in which much of the temperature profile resides. This can be understood by examining the term parentheses in eq 9, namely, 1/Pr + (νT/ν)/Prt. Even though the eddy diffusivity ratio, νT/ν, is very small in this region, the molecular conduction term, 1/Pr, can become of the same order of magnitude or even smaller. Thus, turbulent heat transport (and hence Prt) can be important even though turbulent momentum transport is negligible. The reason for the concern for the variation of νT/ν near the wall is that it is virtually impossible to make experimental measurements in this region. Therefore, the values of Prt in this region must be inferred by indirect methods, e.g., by making complete boundary layer calculations using an eddy diffusivity model, and this cannot be done without accurate data on νT/ν. The Lai and So model provides a convenient way to evaluate νT/ν and permits fairly accurate solution of the momentum equation, i.e., velocity profile. However, the solution to the momentum equation is totally insensitive to values of νT/ν for y+ < 5, quite unlike the situation for the energy equation at high Prandtl number. Consequently, it is necessary to examine νT/ν more carefully in this region before trying to use the lowReynoldsnumber model for the energy equation at high Prandtl number. 3.1. Turbulent Prandtl Number. In the past, a large number of formulations for Prt appeared in the literature. Most of these formulations were inspired by the observation in the 1950s that the experimentally determined heattransfer coefficients for the flow of liquid metal in pipes seemed to fall well below what would be predicted using models such as the Reynolds analogy. It was first suggested by Jenkins78 that this was a result of the relatively high thermal conductivity at very lowPrandtlnumber fluids, and Jenkins78 proposed a model based on the idea that a turbulent “eddy” could lose heat by simple conduction during its motion normal to the mean flow direction, thereby reducing the heat transferred by the turbulent exchange process. Most of the other analyses are based on variations of the same idea, but virtually all contain free constants that need to be evaluated through comparisons with the available experimental data. Three of them can be cited as representative. Experimental values of Prt in the turbulent core for Pr g 0.7 were correlated by Jischa and Rieke31 with expressions such as
Prt ) 0.85 +
0.015 Pr
(10)
Over the purported range of validity of the above equation, Prt varies only from 0.87 to 0.85. An expres
2824 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004
tions for liquid metals. Kays8 found that most experimental values of Prt in the turbulent core could be represented satisfactorily with the empirical equation
Prt )
Figure 2. Performance comparison of the Yakhot equation (eq 12) with the Kays8 equation and the DNS data of Bell,80 Kasagi et al.,81 and Kim and Moin56 as a plot of turbulent Prandtl number versus turbulent Peclet number and the available experimental data on turbulent Prandtl number (s, eq 12;   , Prt ) 0.7/Pet + 0.85)
sion with a broader range is that of Notter and Sleicher25
Prt )
[
1+
1 + 90Pr3/2(µt/µ)1/4
]
10 [0.025Pr(µt/µ) + 90Pr3/2(µt/µ)1/4] 35 + (µt/µ) (11)
The equation is in reasonable accord with eq 10 and with the computed values of Lyons et al.64 for Pr ) 1. However, the methodology was not designed to predict a sharp increase in Prt for large Pr and small µt/µ, in accord with the recently computed values of Papassiliou and Hanratty.79 A more comprehensive model was that of Yakhot et al.10 They presented an analytical solution based on the “renormalization group method” that is apparently devoid of empirical input. The relatively simple equation given by Yakhot et al.10 is worth examining
[
(
][
) (
1  1.1793 Preff 1  1.1793 Pr
(
)
0.65
]
)
1 + 2.1793 Preff 1 + 2.1793 Pr
(
)
0.35
)
1 (1 + νT/ν)
(1 + νT/ν) νT/ν 1 + Prt Pr
Nu )
(
)
Equation 12 covers a very wide range of Prandtl numbers, and at high Pr and high values of νT/ν, it converges to Prt ) 0.85. In this equation, Prt appears as a function of Pr and νT/ν. In Figure 2, eq 12 is shown graphically as Prt versus the product (νT/ν)Pr (i.e., Pet). At high values of (νT/ν)Pr, the curve approaches 0.85. Further, for Pr near 1.00 and higher, there is only slight variation with (νT/ν)Pr, which is, of course, consistent with the observations from the experimental data. At very low values of (νT/ν)Pr, the turbulent Prandtl number becomes very high, consistent with the observa
2Re* Pr T+
(14)
This equation was used for numerical predictions of the Nusselt number using the k model. It was thought desirable to establish the overall energy balance. The total (viscous and turbulent) energy dissipation rate can be calculated from the following equation
exp ) (12)
(13)
This equation is similar in form to an equation suggested by Reynolds,7 but in that case, Prt is an effective average for the entire boundary layer. This result then suggests that Prt is close to a unique function of (νt/ν)Pr. In the present work, eq 12 is used for the expression of Prt. Various experimental data have been analyzed and are plotted in Figure 2, which covers the entire Prandtl number range available from experiments. Some of the experiments were performed with fully developed flow in pipes, and others are from experiments on external flatplate boundary layers with no pressure gradients. The data in the central portion of the graph are primarily the results of experiments with air; the data at low values of (νT/ν)Pr are largely for liquid metals; and the data for (νT/ν)Pr > 100 are for water and higherPrandtlnumber liquids. Note that there appears to be a continuum from the liquid metals to the highPrandtlnumber fluids and that all of the data tend toward about Prt ) 0.85 at high values of (νT/ ν)Pr. To solve eq 9, a model for Prt has to be specified. In the present work, we employ the Yakhot correlation for turbulent Prandtl number. The theory leading to relation 12 determines the turbulent diffusivity in terms of the molecular viscosity and the turbulent viscosity. In particular, it describes the interaction between molecular and turbulent transport, an effect of much significance at low Re and Pr values. Thus, the determination of turbulent heat transfer from eq 12 requires reliable data on turbulent viscosity. The expression for the normalized heattransfer coefficient (Nusselt number) was obtained as follows
where
Preff )
0.7 + 0.85 Pet
∆P × volumetric flow rate mass of liquid in pipe
(15)
Normalizing the above equation ( by uτ3/R and u by uτ), one obtains
exp+ ) fu+3
(16)
where f, the friction factor, is expressed as f ) 0.046(Re)0.2. Re is the Reynolds number based on the diameter and mean velocity of flow in the pipe. The predicted energy dissipation rate consists of viscous and turbulent energy dissipation rates + + ) turbulent + + predicted mean
where
(17)
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2825
( ) ( ) ∂u+ ∂r+
+ + mean ) ν
+ turbulent )ν
2
∂ui ∂ui ∂xj ∂xj
(18)
A comparison between the experimental and the predicted dissipation rates is presented in Table 3. 3.2. Boundary Conditions. Because the flow is axisymmetric, only the boundary conditions at the wall and the symmetry line need to be specified
r ) 0: r ) R:
du dk dT ) ) )0 dr dr dr
u ) k ) 0, q ) κ
dT , dr w ) ν(d2k/dy2)w (19)
Equations 18 along with the boundary conditions (eqs 19) were solved. 3.3. Method of Solution. The solution procedure consisted of two steps: The first step was to solve the momentum equations to obtain the mean axial velocity, turbulent kinetic energy (k), and turbulent energy dissipation rate (). Because the momentum and heat equations present a set of coupled equations, the flow information obtained from the first step was used in the second step, which was to obtain the mean temperature profile. Thus, the solution of the coupled momentumheat model necessitates the solution of four ordinary differential equations for the k model. The governing equations for the k model are ordinary differential equations and therefore can be solved by any iteration scheme for split boundaryvalue problems. A finitevolume technique proposed by Patankar88 was used for the solution of the four coupled nonlinear equations for the k models. The integration was carried out from center to wall. Grid generation is one of the important aspects of the numerical simulation. The robustness of any numerical code depends on the effectiveness and stability of the gridgeneration scheme employed for the investigations. In the present work, a nonuniform grid (in the range of 70175 grid points) was employed, with more than half of the points in the range r/R > 0.9. The predictions obtained from this technique were found to be insensitive to the grids, as doubling the number of grid points changed the solution profiles by less than 1%. Although the influence of grid density on flow is expected, the degree of sensitivity is more for heat transfer. This is especially true for the case of high Prandtl numbers, for which the thickness of the thermal boundary becomes smaller than that of Table 3. Numerical Simulations and the Energy Balance Calculations no. of no. of grid points iterations Reynolds required convergence required energy total + number(s) (range) criteria (range) input (eq 18) 7 400 9 800 22 000 22 500 103 000 149 000 202 000
7585 7585 7585 7585 85175 85175 85175
104105 104105 104105 104105 104105 104105 104105
150250 150250 150250 150250 150250 150250 150250
31.96 32.88 35.61 36.18 41.60 43.16 44.49
28.83 30.08 34.39 35.03 41.53 42.97 43.54
Figure 3. Comparison of mean velocity profiles with the experimental data of Durst et al.5 (1) Universal velocity profile; (2) ERSR model, Fortuin et al.;16 (3) Yu et al.;30 (4) CFD for Reynolds number ) 7442.
Figure 4. Comparison of the shear stress profiles with the experimental data of Eggles,71 the DNS, Churchill,28 and the present CFD model.
the viscous sublayer. Hence, the heat transfer becomes confined to a very small region near the wall. To effectively capture the steep temperature gradients in this region, we selected suitably fine grids. The numerical details are given in Table 3. 4. Results and Discussion 4.1. Comparison of the Numerical Predictions with the Experimental Data. To test the accuracy of predictions as a result of flow simulation, an extensive comparison was made with the experimental data of Durst et al.5 A comparison of the predicted mean axial velocity (Re ) 7442) using the models of Yu et al.30 and Fortuin et al.16 and the present CFD model is shown in Figure 3. It can be seen that there is excellent agreement between the predictions of Yu et al.,30 the CFD model, and the experimental data of Durst et al.5 for the region 0.5 < y+ < 100. The model of Fortuin et al.16 can be seen to show some differences in both the nearwall and core regions. Churchill28 devised a very accurate and general correlating equation for the Reynolds stresses based on asymptotic solutions, experimental data, and results obtained by the DNS solution. Figure 4 shows a comparison between the Churchill28 correlation, the present CFD predictions, the DNS results, and the experimental
2826 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004
Figure 5. Comparison of the predicted Nusselt number with the experimental data of Skupinski et al.12 for Pr ) 0.0153 [4, experimental data; s, (1) Lai and So model with eq 12, (2) Lai and So model with Prt ) 0.85, (3) Yu et al.30].
data of Eggles.71 In the bulk region (r/R > 0.8), the predictions of all three models are very close to each other and agree with the experimental data. In the region 0.08 < r/R < 0.23, the DNS predictions can be seen to be more accurate. However, very close to the wall, the CFD, Churchill,28 and DNS predictions can be seen to be essentially the same. To validate the lowReynoldsnumber k model for heat transfer, fully developed pipe flow temperature measurements with constantwallheatflux boundary conditions were chosen. Unlike the measurements of flow quantities, such as the mean velocity, accurate experimental data on the temperature profiles are scarce. Moreover, there are very few cases in which such measurements are reported for high Prandtl numbers. In the present work, the lowPrandtlnumber experimental data of Skupinski et al.12 for liquid NaK (Pr ) 0.0153) and of Fuchs13 for liquid Na (Pr ) 0.007) and the highPrandtlnumber experimental data of Monin and Yaglom11 (1 < Pr < 106) were used so as to cover a wide range. 4.2. Heat Transfer at Low Prandtl Number. The prediction of turbulent heat transfer in the lowPrandtlnumber range is a crucial test for the model. This is because of the strong influence of the molecular Prandtl number on the value of Prt for fluids with very low Prandtl numbers (liquid metals). Figure 5 shows the comparison between a CFD predictions and the experimental data of Skupinski et al.12 for NaK, (Pr ) 0.0153). The simulations were also carried out using a constant Prt value of 0.85. Figure 6 shows a comparison between experimental data of Fuchs13 for Pr ) 0.007 and the CFD model using eq 12 at constant Prt ()0.85). As can be seen from Figure 6, the agreement between the predicted Nusselt number and the experimental data is good. The inclusion of the Prt correlation (eq 12) can be seen to improve the prediction of the Nusselt number. The use of constant Prt (0.85) decreases the predicted Nusselt number. The same observations were reported by Kays.8 Further, from Figures 5 and 6, close agreement can also be seen between the predictions of the present CFD model and those of Churchill.28 In Figure 7, the CFDpredicted temperature profiles for liquid mercury (Pr ) 0.02), a NaK eutectic mixture (Pr ) 0.029), and air (Pr ) 0.7) are plotted for pipe flow at Re ) 1 49 000. The numerical details are given in Table 1. The agreement can be seen to be excellent.
Figure 6. Comparison of the predicted Nusselt number with the experimental data of Fuchs13 for Pr ) 0.007 (2, experimental data; s, (1) Lai and So model with eq 12, (2) Lai and So model with Prt ) 0.85, (3) Yu et al.30].
Figure 7. Comparison of the mean axial temperature predictions of the k model with the experimental data of Buhr et al.84
4.3. Heat Transfer for High Prandtl Number. At very high Prandtl numbers, the temperature profiles move closer to the wall and almost entirely inside the range of y+ ) 5. However, the flow is still very much governed by a turbulent exchange process as it has already been seen that the eddy diffusivity for heat, (νT/ ν)/Prt, though small, is still substantially greater than the molecular diffusivity for heat, 1/Pr. Additionally, the distance from the wall influences the value of the turbulent Prandtl number, which tends to increase Prt close to the wall. This increase in Prt near the wall is especially important for highPrandtlnumber fluids because of the very thin thermal boundary layer. Outside this thin layer near the wall, the turbulent Prandtl number seems to be constant for Pr > 1. Different model results of the temperature profile for Pr ) 95 are illustrated in Figure 8, together with the experimental data of Kader.14 The present model and the Yu et al.30 predictions are in agreement with the experimental data in the nearwall region. The ERSR model (Fortuin et al.16) overpredicts the experimental data in both the nearwall and core regions. The predictions from all of the models compare well with the experimental data up to y+ < 1, the difference widening with increasing wall distance. Figure 9 depicts the comparison between the CFD predictions and the experimental data of Monin and Yaglom11 (1 < Pr <
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2827
Figure 8. Comparison of the mean axial temperature predictions of the k model with the experimental data of Kader14 for Pr ) 95.
CP ) specific heat at constant pressure, kJ‚kg1‚°C1 D ) term contained in the k equation (dp/dz)c ) constant axial pressure gradient E ) term contained in the equation f ) friction factor fµ, f1, f2 ) damping functions used in the lowReynoldsnumber k and Reynoldsstress models k ) turbulent kinetic energy, m2/s2 k+ ) normalized turbulent kinetic energy, k/uτ2 q ) constant wall heat flux, W/m2 r ) radial component R ) radius of pipe, m T ) mean axial temperature, °C T* ) temperature based on friction velocity, q/CPDuτ, °C u ) axial velocity, m/s ub ) bulk mean axial velocity of the fluid, m/s uτ ) friction velocity, xR(dp/dz)c/(2F), m/s p ) pressure drop, N/m2 y ) normal distance from the wall, R  r, m z ) axial coordinate, m Greek Symbols ) turbulent energy dissipation rate, m2/s3 w ) turbulent energy dissipation rate at the wall, m2/s3 µ ) molecular viscosity of the fluid, kg‚m1‚s1 µt ) turbulent viscosity of the fluid, kg‚m1‚s1 F ) density of the liquid, kg/m3 ν ) molecular kinematic viscosity of the liquid, m2‚s1 νT ) eddy viscosity of the liquid, m2/s νH ) eddy diffusivity for heat, m2/s σk ) turbulent Prandtl number for k σ ) turbulent Prandtl number for Dimensionless Numbers
Figure 9. Comparison for the plot of the Nusselt number Nu as a function of Prandtl number Pr with the experimental data of Monin and Yaglom.11
106). It can be seen that the predictions are in good agreement. 5. Conclusion (1) The lowReynoldsnumber k model has been used to predict the liquid velocities and turbulent viscosity. The agreement between the predictions and the experimental profiles of axial velocity, kinetic energy, and eddy diffusivity was found to be excellent. The model also establishes a good energy balance. (2) The lowReynoldsnumber model can be used effectively for the prediction of the turbulent heat transfer throughout the entire range of experimentally accessible Prandtl numbers using the Yakhot analysis for the turbulent Prandtl number. The model compares well with the experimental data of Monin and Yaglom11 over Prandtl number range of 1 < Pr < 105. (3) The agreement between the predicted Nusselt number and the experimental data of Skupinski et al.12 for Pr ) 0.0153 and Fuchs13 for Pr ) 0.007 was also found to be excellent. In addition, the predicted temperature profiles also showed good agreement with the experimental data and the predictions of Yu et al.30 near the wall. Notation CFD ) computational fluid dynamics, Cµ, C1, C2 ) turbulence parameters in the k model
Pr ) Prandtl number; Pr ) Cpµ/k Prt ) turbulent Prandtl number (eq 12) Pet ) Peclet number, Pet ) (νt/ν)Pr Ret ) turbulent Reynolds number, Ret ) k2/ν Nu ) Nusselt number, Nu ) 2Rh/k Rey ) turbulent Reynolds number based on y, Rey ) xky/ν Re ) Reynolds number based on the mean velocity; Re ) 2Ru/ν Re* ) Reynolds number based on the friction velocity, Re* ) Ruτ/ν T + ) normalized mean axial temperature; T + ) T/T* y+ ) dimensionless wall distance, y+ ) yuτ/ν u+ ) normalized mean axial fluid velocity, u+ ) u/uτ
Literature Cited (1) Patel, V. C.; Rodi, W.; Scheuerer, G. Turbulence Models for Near Wall and Low Reynolds Number Flows. A Review. AIAA J. 1984, 23, 13081319. (2) Hrenya, C. M.; Boilo, E. J.; Chakrabarti, D.; Sinclair, J. L. Comparison of Low Reynolds Number k Turbulence Model in Predicting Fully Developed Pipe Flow. Chem. Eng. Sci. 1995, 12, 19231941. (3) Thakre, S. S.; Joshi, J. B. CFD Modeling of Heat Transfer in Turbulent Pipe Flows. AIChE J. 2000, 46, 17981812. (4) Thakre, S. S.; Joshi, J. B. Momentum, Mass and Heat Transfer in Single Phase Turbulent Flow. Rev. Chem. Eng. 2002, 18, 83293. (5) Durst, F.; Jovonavic, J. Sender J. LDA Measurement in the Near Wall Region of Turbulent Pipe Flow. J. Fluid Mech. 1995, 295, 305335. (6) Lai, Y. G.; So, R. M. C. On Near Wall Turbulent Flow Modeling. J. Fluid Mech. 1990, 221, 641673. (7) Reynolds, A. J. The Prediction of Turbulent Prandtl number and Schmidt Number. Int. J. Heat Mass Transfer 1975, 18, 10551069. (8) Kays, W. M. Turbulent Prandtl NumbersWhere are We? ASME J. Heat Transfer 1994, 116, 195284.
2828 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 (9) Churchill S. W. A Reinterpretation of the Turbulent Prandtl Number. Ind. Eng. Chem. Res. 2002, 41, 63696401. (10) Yakhot, V. S. A.; Orszag, S.; Yakhot, A. Heat Transfer in Turbulent Fluids1.Pipe Flows. Int. J. Heat Mass Transfer 1987, 30, 15. (11) Monin A. S.; Yaglom, A. M. Statistical Fluid Dynamics; MIT Press: Cambridge, MA. 1982; Vol. 1. (12) Skupinski, E.; Tortel, J.; Vautrey, L. Determination de Coefficients de Convection d’un alliage SodiumPottasium dans une Tube Circulaire. Int. J. Heat Mass Transfer 1965, 8, 937951. (13) Fuchs, H. Heat Transfer to Flowing Sodium. Thesis, Institute fur Reaktorforschung, Wurenlingen, Schweiz, EIR, 1973; Bericht Nr. 241. (14) Kader, B. A. Temperature and Concentration Profiles in Fully Turbulent Boundary Layers. Int. J. Heat Mass Transfer 1981, 24 (9), 154. (15) Reynolds, O. On the Extent and Action of the Heating Surface for Steam Boilers. Proc. Manchester Lit. Philos. Soc. 1874, 14, 7. (16) Fortuin J. M. H.; Edward, E. M.; Hamersma, P. J. Transfer Processes in Turbulent Pipe Flow Described by the ERSR model. AIChE J. 1992, 3, 343. (17) Prandtl, L. Bemerkung uber den Warmeubergang im Rohr; 1928; p 29. (18) von Karman, T. The Analogy Between Fluid Friction and Heat Transfer. Trans. ASME 1939, 61 (8), 705. (19) Martinelli, R. C. Heat Transfer to Molten Metals. Trans. ASME 1947, 69 (8), 947. (20) Levich, W. G. Physciochemical Hydrodynamics; Prentice Hall: Upper Saddle River, NJ, 1962. (21) Lin, C. S.; Moulton, R. W.; Putnam, G. M. Mass Transfer between Solid Wall and Fluid Streams Mechanism and Eddy Distribution Relationships in Turbulent Flow. Ind. Eng. Chem. 1953, 45, 636. (22) Deissler, R. G. Analysis of Turbulent Heat Transfer, Mass Transfer and Friction in Smooth Tubes at High Prandtl and Schmidt Numbers; United States National Advisory Committee on Aeronautics (NACA) Report 1210; NACA: Washington, DC, 1955. (23) Friend, W. L.; Metzner, A. B. Turbulent Heat Transfer inside Tubes and the Analogy among Heat, Mass and Momentum Transfer. AIChE J. 1958, 4, 393. (24) Wasan, D. T.; Wilke, C. R. Role of Concentration Level of the Nondiffusing Species in Turbulent GasPhase Mass Transfer at Ordinary Mass Transfer Rates. AIChE J. 1968, 14 (4), 577. (25) Notter, R. D.; Sleicher, C. A. A Solution to the Turbulent Graetz ProblemIII. Fully Developed and Entry Region Heat Transfer Rates. Int. J. Heat Mass Transfer 1972, 27, 2073. (26) Mizushina, T.; Ogino, F. Eddy Viscosity and Universal Velocity Profile in Turbulent Flow in Straight Pipe. J. Chem. Eng. Jpn. 1970, 3 (2), 166. (27) Edwards, M. F.; Smith, R. The Integration of the Energy Equation for Fully Developed Turbulent Pipe Flow. Trans. Inst. Chem. Eng. 1980, 260. (28) Churchill, S. W. New Simplified Models and Formulations for Turbulent Flow and Convection. AIChE J. 1997, 43, 1125. (29) Heng, Ly; Christan, C.; Churchill, S. W. Essentially Exact Characteristics of Turbulent Convection in a Round Tube. Chem. Eng. J. 1998, 71, 163. (30) Yu, B.; Ozoe, H.; Churchill, S. W. The Characteristics of Fully Developed Turbulent Convection in a Round Tube. Chem. Eng. Sci. 2001, 56, 1781. (31) Jischa, M.; Rieke, H. B. In Proceedings of the Advanced Study Institute; Istanbul, Jul 20, 1978. (32) Fortuin, J. M. H.; Klijn, P.J. Drag Reduction and Random Surface Renewal in Turbulent Pipe Flow. Chem. Eng. Sci. 1982, 37, 611. (33) Launder, B. E.; Spalding, D. B. Lectures in Mathematical Models of Turbulence; Academic Press: London, 1972. (34) Speziale, C. G.; Abid, R.; Anderson, E. C. A Critical Evaluation of TwoEquation Models for Near Wall Turbulence. AIAA J. 1992, 30 (2), 324. (35) Wilcox, D. C. Simulations of Transition with a TwoEquation Turbulence Model. AIAA J. 1994, 32, 247255. (36) Sarkar, A.; So, R. M. C. A Critical Evaluation of NearWall Two Equation Models against Direct Numerical Simulation data. Int. J. Heat Fluid Flow 1997, 18, 197208.
(37) Launder, B. E.; Spalding, D. B. Numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 1974, 3, 269289. (38) Jones, W. P.; Launder, B. E. The Prediction of Laminarization with TwoEquation Model of Turbulence. Int. J. Heat Mass Transfer 1972, 15, 301314. (39) Launder, B. E. SecondMoment Closure: Present...and Future? Int. J. Heat Fluid Flow 1989, 10, 282300. (40) Monin, A. S. On the symmetry of turbulence in the surface layer of air. Atm. Oceanic Phys. 1965, 1, 4554. (41) Donaldson, C. D. Calculation of Turbulent Shear Flows for Atmospheric and Vortex Motions. AIAA J. 1971, 10, 412. (42) Zeaman, P.; Lumley, J. L. Buoyancy Turbulent Boundary Layers. In Turbulent Shear FlowI; Springer: Heidelberg, Germany; pp 295306. (43) Baughn, J. W.; Hoffman, M. A.; Launder B. E.; Samaraweera, D. S. A. ThreeDimensional Turbulent Heat Transport in Pipe Flow: Experiment and Model Validation. ASME 1978, 78, WAHT15. (44) Launder, B. E.; Samaraweera, D. S. A. Application of a SecondMoment Turbulence Closure to Heat and Mass Transport in Thin Shear FlowsI. TwoDimensional Transport. Int. J. Heat Mass Transfer 1979, 22, 16311643. (45) Prud’homme, M.; Elghobashi, S. Turbulent Heat Transfer Near the Reattachment of Flow Downstream of a Sudden Pipe Expansion. Numer. Heat Transfer 1989, 10, 349368. (46) Yoo C. J.; So, R. M. C. Variable Density Effects on Axisymmetric SuddenExpansion Flows. Int. J. Heat Mass Transfer 1989, 32, 105120. (47) Launder, B. E. On the Computation of Convective Heat Transfer in Complex Turbulent Flows. J. Heat Transfer 1988, 110, 11121128. (48) Lai, Y. G.; So, R. M. C. NearWall Modelling of Turbulent Heat Fluxes. Int. J. Heat Mass Transfer 1990, 33, 14291440. (49) Martinuzzi, R.; Pollard, A. A Comparative Study of Turbulence Models in Predicting Turbulent PipeFlow. Part I: Algebraic Stress and k Models. AIAA J. 1989, 27, 2936. (50) Pollard, A.; Martinuzzi, R. A Comparative Study of Turbulence Models in Predicting Turbulent PipeFlow. Part II: Reynolds Stress and k Models. AIAA J. 1989, 27, 17141721. (51) Thakre S. S.; Joshi, J. B. A Low Reynolds Stress Reynolds Stress Modeling of Turbulent Pipe Flow: Flow Pattern and Energy Balance. Chem. Eng. Commun. 2002, 189, 759785. (52) Deardorff, J. W. A numerical Study of ThreeDimensional Turbulent Channel Flow at Large Reynolds number. J. Fluid Mech. 1970, 41, 453. (53) Schumann, U. Ein Verfahren zur direkten numerischen Simulation Turbulence Stromungen in Plattenund Ringspaltkanalen and unber seine Anwendung zur Untersung von Turbulenzmodellen. Dissertation, University of Karlsruhe, Karlsruhe, Germany, 1973. (54) Schumann, U. Linear Stability of Finite Difference Equations for ThreeDimensional Flow Problems. J. Comput. Phys. 1975, 18, 465470. (55) Grotzbach G.; Schumann, U. Direct Numerical Simulation of Turbulent Velocity, Pressure, and Temperature Fields in Channel Flows. In Turbulent Shear Flow; Springer; Heidelberg, Germany, 1979; pp 295306. (56) Moin, P.; Kim, J. Numerical Investigation of Turbulent Channel Flow. J. Fluid Mech. 1982, 118, 341377. (57) Orszag, S. A.; Kells, L. C. Transition to Turbulence in Plane Poiseuille and Plane Couette Flow. J. Fluid Mech. 1980, 96, 159205. (58) Orszag, S. A.; Patera, A. T. Secondary Instability of WallBounded Shear Flows. J. Fluid Mech. 1983, 128, 347385. (59) Marcus, P. S. Simulation of TaylorCouette flow. J. Fluid Mech. 1984, 146, 4564. (60) Kim, J.; Moin, P.; Moser, R. Turbulence Statistics in Fully Developed Channel Flow at Low Reynolds Number. J. Fluid Mech. 1987, 177, 133186. (61) Ecjelmann, H. The Structure of the Viscous Sublayer and the Adjacent Wall Region in Turbulent Channel Flow. J. Fluid Mech. 1974, 64, 439459. (62) Kreplin, H.; Eckelmann, H. Behavior of Three Fluctuating Velocity Components in the Wall Region of a Turbulent Channel Flow. Phys. Fluids 1979, 22, 12331239. (63) Antonia, R. A.; Teitel, M.; Kim, J.; Browne, L. W. B. LowReynoldsNumber Effects in a Fully Developed Turbulent Channel flow. J. Fluid Mech. 1992, 236, 579605.
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2829 (64) 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, 9991020. (65) Spalart, P. R. Direct Simulation of a Turbulent Boundary Layer up to Reθ ) 1410. J. Fluid Mech. 1988, 187, 6198. (66) Kristoffersen, R.; Bech, K. H.; Andersson, H. I. Numerical Study of Turbulent Plane Couette Flow at Low Reynolds Number. Appl. Sci. Res. 1993, 51, 337343. (67) Gavrilakis, S. Numerical Simulation of Low Reynolds Number Turbulent Flow through a Straight Square Duct. J. Fluid Mech. 1992, 244, 101129. (68) Unger, F.; Friedrich, R. Large Eddy Simulation of Fully Developed Turbulent Pipe Flow. In Flow Simulation of High Performance Computers I; Hirschel, E. H., Ed.; 1993; Vol. 38, pp 201216. (69) Westerweel, J.; Adrian, R. J.; Eggels, J. G. M.; Nieuwstadt, F. T. M. Measurements with partical image velocimetry of fully developed turbulent pipe at low Reynolds number. In Laser Applications in Fluid Mechanics; 1993. (70) Zhang, Y.; Gandhi, A.; Tomboulides, A. G.; Orszag, S. A. Simulation of pipe flow. In Application of Direct and Large Eddy Simulation to Transition and Turbulence; 1994; pp 17.117.9. (71) Eggels, J. G. M.; Unger, F.; Weiss, M. H.; Westerweel, J.; Adrian, R. J.; Friedrich, R.; Nieuwstadt, F. T. M. J. Fluid Mech. 1994, 268, 175209. (72) Kim, J.; Moin, P. Transport of passive scalars in a turbulent channel flow. I. In Turbulent Shear Flows; Andre, J. C., et al., Eds.; Springer: Berlin, 1989; 8596. (73) Kasagi, N.; Tomita, Y.; Kuroda, A. Direct numerical simulation of passive scalar field in a turbulent channel flow. ASME J. Heat Transfer 1992, 114, 598606. (74) Kasagi, N.; Ohtsubo, Y. Direct Numerical Simulation of Low Prandtl Number Thermal Field in a Turbulent Channel Flow. I. In Turbulent Shear Flows; Durst, F., et al., Eds.; Springer: Berlin, 1993; Vol. 8, pp 97119. (75) Papassiliou, D. V.; Hanratty, T. J. The Use of Lagrangian Methods To Describe Turbulent Transport of Heat From a Wall. Ind. Eng. Chem. Res. 1995, 34, 33593367. (76) Kawamura, H.; Ohsaka, K.; Yamamoto, K. DNS of Turbulent Heat Transfer in Channel Flow with Low to Mediumhigh Prandtl Number Fluid. In Proceedings of the 11th Symposium Turbulent Shear Flows; 1997; Vol. 1, pp 8.78.12. (77) Wikstrom, P. M.; Johansson, A. V. DNS and scalarflux transport modeling in a turbulent channel flow. In Turbulent Heat Transfer; 1998; Vol. 1, pp 6.466.51. (78) Jenkins, R. Proceedings of the Heat Transfer and Fluid Mechanics Institute; Stanford University Press: Stanford, CA, 1951.
(79) Papassiliou, D. V.; Hanratty, T. J. Transport of a Passive Scalar in Turbulent Channel Flow. Int. J. Heat Mass Transfer 1997, 40, 1303. (80) Bell, D. M.; Ferziger, J. H. In Turbulent Boundary Layer DNS with Passive Scalars, in NearWall Turbulent Flows; So, R. M. C., Speziale, C. G., Launder, B. E., Eds.; Elsevier: Amsterdam, 1993; p 327. (81) Kasagi, N.; Tomita, Y.; Kuroda, A. Direct Numerical Simulations of the Passive Scalar Field in TwoDimensional Channel Flow. J. Heat Transfer, Trans. ASME 1992, 114, 598. (82) Roganov, P. S.; Zaboltsky, V. P.; Shishov, E. V.; Leontiev, A. I. Some Aspect of Turbulent Heat Transfer in Accelerated Flows on Permeable Surfaces. Int. J. Heat Mass Transfer 1984, 8, 1251. (83) Fulachier, L. Contribution a Fetude des analogies des champs dynamique et thernique dans une couche limit turbulente. These Doctuer es Sciences, Universite de Provence, Provence, France, 1972. (84) Buhr, H. O.; Carr, A. D.; Balshiser, R. E. Temperature Profiles in Liquid Metals and the Effect of Superimposed Free Convection in Turbulent Flow. Int. J. Heat Mass Transfer 1968, 11, 641. (85) Gowan, R. A.; Smith, J. W. The Effect of the Prandtl Number on Temperature Profiles for Heat Transfer in Turbulent Pipe Flow. Chem. Eng. Sci. 1961, 22, 1701. (86) Hishida, M.; Nagano, Y.; Tagawa M. Transport Processes of Heat and Momentum in Wall Region of Turbulent Pipe. In Flow. Proceedings of the International Heat Transfer Conference; Tien, C. L., et al., Eds.; Hemisphere: Washington, DC, 1986; Vol. 3, p 925. (87) Sheriff, N.; Kane, D. T. Sodium Eddy Diffusivity of Heat Measurements in a Circular Duct. Int. J. Heat Mass Transfer 1981, 24, 205. (88) Patankar, S. V. Numerical Heat Transfer and Fluid Flow; McGrawHill: New York, 1980. (89) Colburn, A. P. A Method of Correlating Forced Convection Heat Transfer Data and a Comparison with Fluid Friction. Trans. AIChE 1933, 29, 174. (90) Sleicher, C. A. Experimental Velocity and Temperature Profiles for Air in Turbulent Pipe Flow. Trans. ASME 1958, 80, 693.
Received for review November 3, 2003 Revised manuscript received March 11, 2004 Accepted March 17, 2004 IE0342311