Modelling of thermal wall boundary conditions with temperature-dependent material properties for use in RANS

International Journal of Heat and Fluid Flow 80 (2019) 108495

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Modelling of thermal wall boundary conditions with temperature-dependent material properties for use in RANS

T

H. Steiner , C. Irrenfried ⁎

Institute of Fluid Mechanics and Heat Transfer, Graz University of Technology, Inffeldgasse 25F, Graz 8010, Austria

ARTICLE INFO

ABSTRACT

Keywords: DNS Thermal boundary conditions Turbulent heat transfer Turbulent Prandtl number P-Function

The present work extends a recently proposed P-function based model for describing the near-wall variation of temperature in forced convective turbulent flow to the case with temperature-dependent material properties. The extension essentially modifies the model formulations for describing the local variation of the turbulent mixing length and the turbulent Prandtl number. Direct Numerical Simulations (DNS) and experimental measurements are carried to provide comprehensive validation data for a wide range of Reynolds numbers, considering molecular Prandtl numbers well beyond unity. The observed good agreement of the predictions with the DNS data and experiments proves the present extended model as a well-suited approach for prescribing reliable thermal boundary conditions in Reynolds Averaged Navier-Stokes (RANS) simulations, assuming temperaturedependent material properties.

1. Introduction The computational modelling of heat transfer is of high importance for many industrial engineering applications. The predicted heat transfer rates in RANS-type simulations strongly depend on an accurate description of the near wall temperature variation. A big challenge arises particularly for molecular Prandtl numbers well above unity, as typically met with organic cooling liquids or oils, where the frequently assumed Reynolds analogy between the momentum and heat transfer breaks down. The P-function based approach has become a very popular concept for modelling the thermal boundary conditions in RANS. Although not explicitely termed as such, this concept was already introduced in the pioneering studies of Prandtl (1910) and Kármán (1939), who related the heat to the momentum transfer for liquids with non-unity Prandtl numbers, using universal laws for the near-wall velocity. A widely applied formulation was proposed by Jayatilleke (1969), who empirically modified and calibrated the analytical P-function of Spalding (1967) using experimental data, especially for molecular Prandtl numbers far beyond unity Pr ≫ 1. Considering the same high Prandtl number regime, Malin (1987) derived a formulation consistent with a logarithmic temperature law originally proposed by Kader and Yaglom (1972). The applicable y+-range and accuracy of the P-function based concept was recently markedly advanced by Irrenfried and Steiner (2017) including appropriate submodels for the near-wall behavior of the

⁎

eddy-viscosity and the turbulent Prandtl number into Spalding’s original analytical formulation. The present work further extends this approach to temperature-dependent material properties, which may lead to considerable variations in the local Reynolds- and Prandtl numbers. The extended model is validated against DNS data and further assessed against experimental data for high Reynolds numbers not amenable to DNS. 2. Mathematical formulation The P-function is introduced relating the enthalpy to the axial velocity in wall coordinates as +

(1)

= PrT, (w+ + P ),

= w¯ /w , non-dimensionalized with the wall friction velocity with and enthalpy, w = ( w/ w )1/2 and h = qw / w w , respectively. The Pfunction is defined as w+

P=

w+ 0

(1 + ) µT

+ qtot

µ

PrT , Pr

+

+ tot

PrT ,

µT

PrT

µ

1 dw˜ +, (2)

involving the local molecular Prandtl number Pr = µcp/ , and its turbulent counterpart PrT = µT / aT , defined as the ratio of the turbulent transfer coefficients for momentum and heat. Considering turbulent pipe flow with constant wall heat flux, the analytically obtained total

Corresponding author. E-mail address: [email protected] (H. Steiner).

https://doi.org/10.1016/j.ijheatfluidflow.2019.108495 Received 19 April 2019; Received in revised form 15 October 2019; Accepted 19 October 2019 0142-727X/ © 2019 Elsevier Inc. All rights reserved.

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

Nomenclature

Greek symbols

A a Am, A0 b C C˜

γ1, γ2 δ Δ ϵ ε η ηθ κ λ μ μe ρ τ φ χ

area thermal diffusivity parameters parameter parameter parameter constant friction coefficient specific heat capacity diameter parameter blending functions enthalpy friction enthalpy, qw/(ρwτ) turbulent kinetic energy prefactor length mixing length Nusselt number P-function turbulent Peclet number, μT/μPrw Prandtl number heat flux radius friction Reynolds number, wτD/ν enthalpy thickness Reynolds number, wτΔ2/ν momentum thickness Reynolds number, wτδ2/ν scaling factor temperature radial, azimuthal and axial velocity Kolmogorov velocity, (νϵ)1/4 friction velocity, (τw/ρw)1/2 wall distance axial direction

Cμ cf cp D E f1, f2 h hτ k K L lm Nu P PeT Pr q r Reτ Re 2 Re 2 S T u, v, w w wτ y z

tot

+

=

+ qtot =

tot

=1

2

w

1 qtot = qw

Re

,

()′ fluctuations ()+, ()*, ()• non-dimensionalized wall units value at y → ∞ ()∞ ()b bulk ()i inner ()mod modeled outer ()o ()tot total turbulent ()T ()w value at the wall () Reynolds average Abbreviations DC DNS DV RANS

y

4 wb+Re

0

1

(1

2y˜+ Re

2y + /Re

) w dy˜ +

+

,

Am+ = A0 [1

(4)

PrT =

exp

y+ µb*

Am+

.

(6)

The bulk viscosity ratio µb* = µb / µ w , computed with the cross-sectional average

µb =

A

µdA

bA02 2

,

(8)

1 + CPeT

(CPeT )2 1

2

exp

1 CPeT

1 2

,

(9)

with

(5)

using a van Driest-type dependence of the mixing length on y+ , written as

lm+ = y + 1

exp( y+ / C )]0.5 , C =

with the constants A0 and b, which provides the correct near wall 0. asymptotics, µT / µ y +3 for y+ The variation of the local turbulent Prandtl number is modeled as

respectively, where y+ = w w y /µ w denote the non-dimensional distance from the wall, Re = w w D/ µ w the wall friction Reynolds number, and wb+ = wb/ w the non-dimensional bulk velocity, equivalent to the volumetric flow rate. The local eddy viscosity ratio is computed from

µT / w +2 dw+ = lm µ µ/ µ w dy+

Cases with constant fluid properties Direct Numerical Simulations Cases with temperature-dependent fluid properties Reynolds Averaged Navier-Stokes

has been introduced in the damping term of (6) accounting for a temperature-dependent deviation of the molecular viscosity from the wall value. As suggested by Grifoll and Giralt (2000), the damping parameter Am+ is not set constant, but computed as a function of y+ as

(3) +

dissipation rate, ui / xj ui / xj dynamic scale thermal scale von Kármán constant thermal conductivity dynamic viscosity effective dynamic viscosity density shear stress azimuthal direction ethalpy difference

Subscripts/Superscript

shear stress and heat flux read

y+

parameters momentum thickness enthalpy thickness

1

=

2

=

µb* PrT , + 0.166 Pr0.7 w 2

1 PrT ,

1

(10)

(11)

and the model parameters C = 1.5, PrT, = 0.85 . These parameters differ from the proposition given by a previous study of Irrenfried and Steiner (2017), where they were calibrated solely for constant fluid properties. The present modified parameter setting was proposed so as

(7) 2

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

to be well applicable to both constant and variable fluid property cases, instead. The considered temperature-dependent variation of the fluid properties enters the formulation (9) through the turbulent Peclet number, defined as PeT = µT / µPrw, involving the local eddy viscosity ratio obtained from (5). The asymptotically approached wall value, which is determined by the model parameter γ1, is modified as well including the bulk viscosity ratio into the definition (10). The definite integral of (2) is evaluated upon substitution of all y+-dependent inputs (3), (4), (5), and (9), as defined above, relating the axial velocity w+ to the wall distance y+ through

dw+ = dy+

2 1+

the unheated equalization section the temperature is supposed to recover a uniform distribution, representing the enthalpy flow equivalent mixing temperature Tm,2. Constant wall heat flux boundary conditions are realized based on Ohm’s law, by electrically short circuiting the heated test section of axial length LHT = 2 m with a high current transformer. The data were measured for the same working fluid in dedicated experiments described in more detail in Steiner et al. (2017). The skin-friction coefficient is determined from the axial pressure drop as

cf =

+ tot

µ/µw

1 + 4 lm+2

+ tot

/

.

D 2 b wb

,

(13)

where the difference Δp is measured with pressure probes, axially located at a distance Δzp. The bulk velocity wb is obtained from the mass flow rate m provided by the mass-flowmeter. The Nusselt number is measured at the end of the heated test section using the expression

w

(µ / µ w ) 2

1 p 2 zp

(12)

3. Direct numerical simulations

Nu =

D

=

w

Direct Numerical Simulations (DNS) of dynamically and thermally fully developed turbulent pipe flow with constant wall heat flux were carried out for model validation. The simulations considered three cases with temperature-dependent fluid properties (indicated by “DV”), associated with molecular Prandtl number Prw = µ w cp, w / w = 20 for Reynolds numbers Re = 360/500, and with Prw = 50 for Re = 360, always based on wall conditions. For the actually assumed heat transfer oil, the temperature variation significantly influences the viscosity, as indicated by Fig. 1, showing the mean radial variation of the viscosity relative to the wall value, as obtained from the DNS for the considered cases. Corresponding cases with constant material properties (indicated by ″DC″) were simulated as well for comparison. The present DNS solves the conservation equation for mass, momentum and energy directly in cylindrical coordinates using a 4th order accurate Finite-Volume discretization in space and 2nd order accurate Adams-Bashforth discretization in time. The axial length of the computational domain is five diameters, imposing periodic boundary conditions at the in- and outlet. The computational mesh consists of 256 × 512 × 1024 cells in radial, azimuthal and axial direction, respectively. Fig. 2 shows the local resolution, measured in terms of the smallest relevant dynamic and thermal scales, the Kolmogorov scale η and Batchelor scale = (Pr) 1/2 , respectively, for the most demanding case, associated with Prw = 50 . The shown maximum values + + rmax / = 3.1, (r + )max / = 8.933 and z max / = 7.110 for the radial, azimuthal and axial direction, respectively, remain acceptably close to unity, and are well comparable to DNS, which were presented in literature by other authors, see, e.g. Zonta et al. (2012), Lee et al. (2013).

qw D w (Tw

Tm,2)

.

(14)

The wall heat flux is determined here from the heat-up of the fluid along the electrically heated test section, written as

qw =

mcp, b (Tm,2

Tm,1)

D LHT

.

(15)

The required wall temperature Tw, as well as the mean mixing temperatures, Tm,1 and Tm,2, are measured using PT-100 temperature sensors. The accuracy of these sensors, being T = ±0.2 K, mainly determines the uncertainty of the measured Nusselt number. Accordingly, the error in the experimental values for Nu ranges from 4% to 20%, depending on the measured temperature difference, (Tw Tm,2), input into the denominator of (14). The axial section, which is considered by the present DNS is located at the end of the heated test section, as indicated in Fig. 3. Two different experimental series are used for model validation, keeping the molecular Prandtl number constant, with Prw = 20 and Prw = 50, respectively, based on the wall conditions at the end of the heated test section, while varying the bulk Reynolds numbers Reb = 4m / µb D from the laminar-turbulent transition regime up to fully turbulent conditions. 5. Results and discussion The significant viscosity variations indicated in Fig. 1 have a strong influence on the turbulence statistics, as illustrated in Figs. 4 and 5 for the momentum and heat flux budgets, rewritten as + tot

4. Experimental setup for validation For an assessment of the present model also in high Reynolds number cases, which are beyond the scope of DNS, an experimental data base was acquired for validating in particular the predicted skinfriction coefficient and the Nusselt numbers. The specially designed test facility used for evaluation data acquisition basically follows an experimental heated pipe flow setup presented in literature by Ghajar and Tam (1994). The complete experimental loop is shown in Fig. 3. It essentially consists of a reservoir, a rotation speed controlled circulation pump, a mass-flowmeter, the test section, and a heat exchanger. The operating fluid is thermally conditioned by the water cooled heat exchanger and an electrical heating rod inside the reservoir. The test section consists of a long straight cylindrical pipe made of stainless-steel with an inner diameter di = D = 12 mm and an outer diameter da = 15 mm . It entails three sub-sections. Along the unheated entrance section and the heated test section, the flow is supposed to reach dynamically and thermally fully developed conditions, respectively. Along

=

µ w+ µ w y+ lam

µ w+ µw y+ µ

u+ w+ w turb

Fig. 1. Mean relative molecular viscosity vs. y+ . 3

(16)

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

Fig. 4. Comparison of the laminar/turbulent shear stress for the DC and DV cases. Fig. 2. Spatial resolution normalized by the Kolmogorov scale η and the smallest turbulent thermal scale ηθ.

Fig. 3. Schematic of the experimental facility.

+ qtot =

1 Prw

cp, w w

cp

qlam

+

y+

1 Prw

w

q

cp, w

+

cp

y+

u+

+

w

qturb

Fig. 5. Comparison of the laminar/turbulent heat flux for the DC and DV cases.

(17)

turbulent heat flux is still not enhanced like in the DC cases, where the transition from the diffusive sublayer to the fully turbulent inner layer is shifted towards the wall for increasing Prandtl number, as indicated by comparison of DC360/20 with DC360/50 in Fig. 5. In the corresponding variable fluid property cases DV360/20 and DV360/50 the transition is consistently shifted towards the inner core flow region instead. The already noted reduced turbulent mixing motion, in particular of the radial component, evidently overcompensates the increased enthalpy fluctuations. This effectively increases the thickness of the diffusive sublayer as well. The adequate representation of the significantly modified flux budgets, as observed for the DV cases, in the modelling of the turbulent

respectively. The non-linear contributions arising from the fluctuating molecular transport coefficients turned out to be insignificantly small, and are therefore not considered. The increase of the molecular viscosity towards the colder core flow region reduces the turbulent mixing motion in the upper buffer layer, as indicated by the considerably decreased contribution from the turbulent shear stress around y+ 10 . The consistently increased laminar stress component makes the thickening of the viscous sublayer as compared to the constant property case (DC) evident. Furthermore, the increased viscosity implies an increased local molecular Prandtl number, which basically leads to steeper ethalpy gradients, due to the locally increased thermal resistance, and, hence, to stronger turbulent ethalpy fluctuations χ′. The resulting 4

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

transport coefficients for momentum and heat is of key importance for the accuracy of model predictions. Due to the fairly low Reynolds number range considered in this study, the model parameters required for the mixing length computed from Eq. (6) were chosen somewhat different from the commonly used standard setting, being, e.g., for the Kármán constant around = 0.4 . This aspect has also been extensively discussed in previous experimental studies (Durst et al. (1995)) and DNS studies (Wu and Moin (2008)). The presently prescribed model parameters are summarized in Table 1. Furthermore, it is recommended to employ the commonly used standard setting for the mixing length modeling, when a fully turbulent flow is considered with a Reynolds number above Reτ ≥ 750. Fig. 6 compares the predicted turbulent viscosity ratio against DNS data for the cases with Re = 360/Prw = 50, and Re = 500/Prw = 20, respectively. For DV360/50, where the viscosity most significantly increases with the distance from the heated wall, the eddy-viscosity is notably reduced, as shown by comparison with the corresponding DNS results for DC360/50, where the wall values of the fluid properties are assumed throughout the entire domain. The model based on the presently proposed mixing length formulation (6)-(8) obviously captures this reduction very well, including the correct description of the near wall asymptotics, where the turbulent viscosity ratio scales with y+3 . For the higher Reynolds number case DV500/20 the reduction of the eddy viscosity is less pronounced, which is due to the smaller wall superheat Tw Tm, leading to less variation of the material properties. The present model produced good agreement also for this case. The reliable modelling of the mixing length consistently translates into accurate predictions of the axial velocity, computed based on Eq. (12). This is highlighted in Fig. 7, comparing the results for the cases with Re = 360/Prw = 50 and Re = 500/Prw = 20 . Very next to the wall, where the local temperature is close to the wall reference value, equally prescribed for both the DC and the DV cases, the variations are on top of each other. Inside the buffer layer the effectively thickened viscous sublayer is clearly indicated by the shift of the DNS profiles to higher y+. Including the bulk viscosity into Eq. (6) the present model is capable of reproducing the constant and the variable DNS results in the buffer region very accurately. Some deviations from the DNS results appear in the central core region, where the assumption of a linearly increasing mixing length does not hold. The observed reduction of the eddy-viscosity near the wall due to the increase of the molecular viscosity also notably affects the near wall behavior of the turbulent Prandtl number. According to its definition, it can be computed from the statistically averaged DNS resluts as

PrT =

u +w + d + / dy + . u + + dw +/dy+

and DV cases generally indicates lower asymptotic wall values of PrT for the latter. The decrease in the turbulent shear stress is apparently not accompanied by an equivalent decrease in the turbulent heat flux, so that the constant K determining the limit of PrT in Eq. (19)is effectively reduced. Due to the modification of the model parameter γ1 with the bulk viscosity ratio µb* , the proposed model formulation (9) is capable to reflect this trend in the near wall asymptotics of PrT, shown in Fig. 8 for the DV cases. Fig. 9 compares the relative enthalpy variations obtained from the DNS for constant and variable material property against the corresponding model predictions from Eq. (1). The already noted increase of the thermal resistance near the wall due to the local increase of the molecular Prandtl number in the DV cases is clearly indicated by the steeper gradients around y+ 2 . Together with the reduced turbulent convective mixing, being also caused by the local increase in molecular viscosity, this finally produces consideraby higher wall superheats, as seen from the higher levels of relative enthalpy, which are reached in the inner bulk flow region for the DV cases. For the higher Reynolds the bulk levels of relative enthalpy are somewhat reduced, as expected due to the more intense turbulent motion. The predictions of the present model agree evidently very well with the DNS data also for the cases with temperature-dependent material properties. Differently from many other P-function based approaches the present model is not restricted to the inertial subrange, associated with the log-law. It apparently provides a very accurate description of the whole near-wall region, including the viscose/diffusive and the buffer layers. 5.1. Comparison against experimental data Due to the difficulties in measuring accurately the enthalpy variation in the near wall region, only very little reliable measurement data are available for validation. Hollingsworth et al. (1989) published near wall temperature measurements of a fully turbulent boundary layer developping at the bottom of an open channel flow. Despite the differences to the presently considered fully developped channel flow configuration, these data could still be adopted for model validation. The comparability of the model predictions with the data was achieved, assuming the wall friction velocity and enthalpy, wτ and hτ, respectively, so as to provide the momentum and enthalpy thickness based Reynolds numbers, which actually characterize the heated turbulent boundary layer flow conditions in the experiments. They are defined as

Re

u +w + Prw = K (Prw) Prw . 0u+ +

y

0

y

bw

=

µb

2

, Re

2

bw

=

2

(20)

µb

with the momentum and enthalpy thickness written as

(18)

The so obtained radial variation of the turbulent Prandtl number is qualitatively very similar between the DC and DV cases, as seen in Fig. 8. Inside the fully turbulent inner core flow region the turbulent Prandtl number remains around unity, where the DV cases tend to attain a slightly lower level. In the viscous/diffusive subrange all cases show a considerable increase towards the wall, which is in good accordance with various former studies, e.g., Kays and Crawford (1980), Schwertfirm and Manhart (2007) and Irrenfried and Steiner (2017). The shown asymptotic behavior is also consistent with analytical con0, siderations based on Taylor series expansions in the limit y+ which yield

lim PrT = lim + +

2

2

=

2

=

0

w w

0

w h w hw

1

w w

dy ,

h dy, h

respectively. Both length scales were equivalently computed with the present pipe flow model, assuming the mean centerline conditions as free-stream conditions, varying the parameters wτ and hτ until the targeted experimental setting associated with the data of Hollingsworth et al. (1989), Reδ2 ≈ 1400 and ReΔ2 ≈ 300, was obtained. The considered fluid was water, and the temperature variation

(19)

Table 1 Model parameters for turbulent mixing length .

This expression suggests an increase of the asymptotic near-wall limit for increasing molecular Prandtl number at the wall Prw, albeit the analysis does not reveal the actual dependence of the prefactor K on Prw. The present DNS results as well as the availabe DNS data from literature support this tendency. A cross comparison between the DC

Case Re = 360 Re = 500 Reτ ≥ 750

5

κ 0.34 0.36 0.40

A0 33 33 33

b 0.0006 0.0008 0.0010

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

Fig. 6. Local eddy viscosity ratio vs. y+ .

Fig. 8. Turbulent Prandtl number vs. y+ .

Figs. 11 and 12, the model predicts both parameters in fairly good agreement with the experimental data. The maximum relative error remains beneath 7% for the whole range of bulk Reynolds numbers considered in the present experiments. The model performs evidently very well also in the case of higher Reynolds numbers, which are markedly beyond the scope of the present DNS. A minor trend of discrepancies can be seen for the highest Prandtl number Prw = 50 for increasing the Reynolds number, still the deviation stay within the error margins of the employed measurement approach. The predictions for both parameters cf and Nu, as obtained from the present DNS with temperature-dependent material properties, are also indicated in the comparison at the lower end of the shown range of Reynolds numbers. The DNS data are evidently in very good agreement with the experiments. Fig. 7. Non-dimensional axial velocity w+ vs. y+ for Re = 360/500 (upper/ lower subfigure).

6. Limits of the present model and possible further advancements The proposed model is basically devised for attached boundarylayer type flows, where all constitutive model equations are formulated in standard wall units, scaled by the wall friction velocity and temperature. As such, the developed universal description clearly hits its limits in the case of non-equilibrium detached flow, where the wall friction velocity becomes zero at the re-attachment points. Therefore, the logarithmic wall function based posture of boundary condition in RANS generally follows the methodology of Launder and Spalding (1974), which avoids this singularity by assuming the wall friction velocity in terms of the turbulent kinetic energy k, written as

inside the boundary layer was Tw T = 4 K, with a boundary layer averaged Prandtl number Prb = 5.9. Fig. 10 compares the measurements of Hollingsworth et al. (1989) with the present model, represented in wall units. The agreement of the predictions for the axial velocity as well as the temperature is obviously very good. Some deviations occur near the upper limit of the buffer layer, beyound y+ 10, which, however, fall still inside the range of experimental uncertainty of the measured local temperature. Unlike the experiments of Hollingsworth et al. (1989) the measurements carried out with the present facility described in Section 4 only provide global transfer parameters for momentum and heat, such as the skin-friction coefficient and the Nusselt number. As shown in

w = Cµ1/4 k1/2,

(21)

with the constant Cµ = 0.09 . Using this expression as velocity scale 6

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

Fig. 9. Non-dimensional enthalpy

+

vs. y+ .

The wall shear stress is further related to the local streamwise velocity and wall distance through an “effective” viscosity μe, such that w

= µe

w¯ , y

with

µe =

w

ln(E

Cµ1/4 k1/2y 1/4 1/2 w Cµ k y /µ w )

.

(23)

The validity of the relations (21) and (23) basically restricts this concept for prescribing boundary conditions to the log-law region. The limits of the y+ range, where the inverse 1/ k + = w 2/ k is close to Cµ1/2 = 0.3, and hence, the relation (21) may be assumed, is clearly shown by the DNS results in Figs. 13 and 14. Using consistently the rescaled wall units defined in (22) instead of their ( )+-counterparts in the present model formulation, would certainly circumvent the singularity for w = 0, providing wall boundary conditions in terms of w* and χ* dependent on y* inside the log-law region, as typically imposed in RANS with standard wall function boundary conditions. The scope of the model could be further extended towards the viscous near-wall region associated with small y+ < 10 adopting an approach of Bredberg et al. (2000), who blended the high Reynolds number (HRN) fully turbulent contribution with a low Reynold number (LRN) viscous near-wall contribution written as

w* = f1 y * + (1

* = f2 Prw y * + (1

Fig. 10. Predicted velocity and temperature variation compared against measurements of Hollingsworth et al. (1989).

1 1 4k2

w yCµ

µw

, w* =

w 1 1 C4 k2

µ

, * =

hw h qw

1 1 4k2

w Cµ

* . f2 ) HRN

(24) (25)

The blending functions, which reduce either of the terms, are defined as simple exponential function

˜ / µ), f = exp( C˜ Prw µ / µ ) f1 = exp( Cµ T T 2

instead of w = ( w/ w )1/2 , the non-dimensional wall units are rewritten as

y* =

* f1 ) wHRN

(26)

with the model constant C˜ = 3. Figs. 15 and 16 exemplarily show velocity and enthalpy variations predicted by this blending function based approach for the case DV500/ * and 20, respectively. The high Reynolds number contributions, wHRN * , were taken here from the present P-function based model along HRN

(22)

Fig. 11. Skin-friction coefficient vs. bulk Reynolds number Reb for the measurement series with Prw = 20 and Prw = 50 . 7

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

Fig. 12. Nusselt number vs. bulk Reynolds number Reb for the measurement series with Prw = 20 and Prw = 50 .

Fig. 13. Radial variations of inverse of the non-dimensional turbulent kinetic energy from DNS results, DC cases.

Fig. 16. Enthalpy difference predicted from (25) compared against DNS results for the case DV500/20 .

Fig. 14. Radial variations of inverse of the non-dimensional turbulent kinetic energy from DNS results, DV cases.

Fig. 17. Rescaled DNS profiles for axial velocity w• vs rescaled wall distance y• for Reynolds numbers Re = 360 and 500, Prw = 20 .

Fig. 15. Axial velocity predicted from (24) compared against DNS results for the case DV500/20 .

Fig. 18. Rescaled DNS profiles for enthalpy difference χ• vs rescaled wall distance y• for Reynolds numbers Re = 360 and 500, Prw = 20 . 8

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried

Abe et al. (1994) introduced a scaling based on the Kolmogorov velocity scale defined as (30)

w = ( )1/4 ,

The non-dimensional wall units based on this scaling would accordingly read

y• =

with its constitutive submodels, all evaluated as functions of y* instead of y+. For the sake of comparability to the DNS results, the wall units predicted by (24) and (25) were transformed back into ( )+-units using the wall friction velocity based on the wall shear stress from assumption (23) as reference scale. Accordingly, the rescaling to ( )+-units written as

w* , S

+ mod

= * S,

(27)

w , mod/w Cµ1/4 k +1/2

,

, mod

w

=

Cµ1/4 w+k +1/2 ln(E Cµ1/4 k +1/2y+ )

(28)

1/2

,

=

h w h¯ qw

ww

(31)

The presently proposed extension of the P-function based model, originally devised and validated for constant fluid properties, is proven to yield acceptably accurate predictions also for cases with significantly variations of the fluid properties, especially of the viscosity inside the near-wall region. As such, the present extended model carries over its major benefits to cases with non-constant fluid properties, providing a reliable continuous description of the complete thermal boundary layer, which allows for an accurate and flexible setting of thermal wall boundary conditions in RANS. Applying a wall-function based concept developed for boundary-layer type flow, the present approach can not be expected to predict accurately non-equilibrium flow conditions, as occuring in detached flow. The comparison of the predicted Nusselt numbers against experiments further indicates the applicability of the model in the range of higher Reynolds number beyond the limits of DNS.

involving a ratio of the wall friction velocities written as

w

•

7. Conclusions

applies a scaling factor

S=

y /µ w , w• = w¯ / w ,

Figs. 17 and 18 exemplarily show rescaled DNS data for velocity and enthalpy difference using the Kolmogorov velocity as reference scale, as computed from the corresponding variations of the local dissipation rate presented in Fig. 19. The scaling evidently provides an universal description, that is limited to the region very next to the wall, as indicated by the collapse of the curves for the two different Reynolds numbers below y• ≈ 10 (corresponds here to y+ 3 ÷ 5 in ( )+-units). This limitation to the viscous near-wall region together with the required knowledge of the local dissipation rate for the Kolmogorov velocity scale essentially prohibit the use of this scaling in the present Pfunction based model, as shown possible with the scaling concept of Launder and Spalding (1974).

Fig. 19. Radial variation of dissipation rate from DNS.

+ + ymod = y *S , wmod =

ww

(29)

where the local DNS-data are substituted for w+ and k +, as dependent on the wall distance y+. The rescaled predictions evidently exhibit notably more discrepancy than their counterparts shown in the former comparisons in Figs. 7 and 9. The differences appearing in the log-law region essentially indicate that the assumptions (21) and (23) applied for the rescaling to ( )+-units do not perfectly match the DNS data, which were substituted into (28) and (29) in an a priori sense for computing the scaling factor S. In RANS (28) and (29) would be rather evaluated with the solution from a transport equation for the turbulent kinetic energy. Very close to the wall, the blending functions quite effectively dampen the errorenous contribution from fully turbulent high Reynolds number component (HRN), so that the laminar solution is predicted fairly well. In summary, it can be stated that applying straight-forwardly the methodology of Launder and Spalding (1974) for defining universal wall units certainly avoids the singularity associated with zero wall shear stress. However, the gain in robustness is achieved at the cost of accuracy, which is due to the limitations of the assumptions to be made for computing the wall shear stress velocity and the underlying wall shear stress at the first near-wall grid point. It is further noted, applying the blending approach given by (24) and (25) basically makes the model applicable also very close to the wall. The weighted superpostion of a laminar and turbulent attached flow solution is certainly not expected to describe accurately nonequilibrium detached flow conditions. Highly non-equilibrium conditions like the flow near reattachment are therefore beyond the scope of the present model and would rather require near-wall low Reynolds number models as, e.g., proposed by Abe et al. (1994, 1995). These low Reynolds number approaches generally solve the governing set of RANS-equations deeply into the viscous/diffusive sublayers down to y+ near unity and below, together with transport equations for turbulent kinetic energy and its dissipation rate, as well as for their thermal counterparts, the temperature variance and its dissipation rate. For avoiding the singularity associated with zero wall shear stress,

Declaration of Competing Interest None. Acknowledgments The financial support of the Austrian Research Promotion Agency (FFG), the Virtual Vehicle competence center K+viF, and AVL List GmbH are gratefully acknowledged. The simulations were conducted using the Vienna Scientific Cluster (VSC) and the D-Cluster in Graz. References Abe, K., Kondoh, T., Nagano, Y., 1994. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flowsi. flow field calculations. Int. J. Heat Mass Transf. 37 (1), 139–151. Abe, K., Kondoh, T., Nagano, Y., 1995. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flowsi. Thermal field calculations. Int. J. Heat Mass Transf. 38 (8), 1467–1481. Bredberg, J., Peng, S.-H., Davidson, L., 2000. On the wall boundary condition for computing turbulent heat transfer with k-ω models. Proc. ASME Heat Transfer Division, HTD-Vol. 366-5, Nov. 5–10,2000, Orlando, Florida, USA. American Society of Mechanical Engineers, pp. 243–250. Durst, F., Jovanovic, J., Sender, J., 1995. LDA measurements in the near-wall region of a turbulent pipe flow. J. Fluid Mech. 295, 305–335. Ghajar, A.J., Tam, L.-M., 1994. Heat transfer measurements and correlations in the transition region for a circular tube with three different inlet configurations. Exp. Therm. Fluid Sci. 8 (1), 79–90. Grifoll, J., Giralt, F., 2000. The near wall mixing length formulation revisited. Int. J. Heat

9

International Journal of Heat and Fluid Flow 80 (2019) 108495

H. Steiner and C. Irrenfried Mass Transf. 43 (19), 3743–3746. Hollingsworth, D.K., Kays, W.M., Moffat, R.J., 1989. The measurement and prediction of heat transfer in a turbulent boundary layer in water. Symposium on Turbulent Shear Flows, Stanford, CA. Vol. 2. Irrenfried, C., Steiner, H., 2017. DNS based analytical P-function model for RANS with heat transfer at high Prandtl numbers. Int. J. Heat Fluid Flow 66 (Supplement C), 217–225. Jayatilleke, C.L.V., 1969. The influence of prandtl number and surface roughness on the resistance of the laminar sublayer to momentum and heat transfer. Prog. Heat Mass Transf. 1, 193–321. Kader, B., Yaglom, A., 1972. Heat and mass transfer laws for fully turbulent wall flows. Int. J. Heat Mass Transf. 15 (12), 2329–2351. Kármán, T., 1939. The analogy between fluid friction and heat transfer. Trans. ASME 61, 705–710. Kays, W.M., Crawford, M.E., 1980. Convective Heat and Mass Transfer. Tata McGraw-Hill Education. Launder, B., Spalding, D., 1974. The numerical simulation of turbulent flows. Comput. Methods Appl. Mech.Eng. 3, 269–289.

Lee, J., Jung, S.Y., Sung, H.J., Zaki, T.A., 2013. Effect of wall heating on turbulent boundary layers with temperature-dependent viscosity. J. Fluid Mech. 726, 196–225. Malin, M., 1987. On the calculation of heat transfer rates in fully turbulent wall flows. Appl. Math. Modell. 11 (4), 281–284. Prandtl, L., 1910. Eine Beziehung zwischen Wärmeaustausch and Strömungswiderstand der Flüssigkeiten. Phys. Z. 11, 1072–1078. Schwertfirm, F., Manhart, M., 2007. DNS of passive scalar transport in turbulent channel flow at high schmidt numbers. Int. J. Heat Fluid Flow 28 (6), 1204–1214. Spalding, D., 1967. Monograph on Turbulent Boundary Layers. Imperical College: Mechanical Engineering Department Report TWF/TN 33. Steiner, H., Irrenfried, C., Brenn, G., 2017. Near-wall determination of the turbulent Prandtl number based on experiments, numerical simulation and analytical models. Proc. HEFAT 2017 (Portoroz, Slovenia, July 17–19). Wu, X., Moin, P., 2008. A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow. J. Fluid Mech. 608, 81–112. Zonta, F., Marchioli, C., Soldati, A., 2012. Modulation of turbulence in forced convection by temperature-dependent viscosity. J. Fluid Mech. 697, 150–174.

10