12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

HEAT FLUX MODELING IN WALL-BOUNDED TURBULENT FLOWS BY
REYNOLDS STRESS MODELS WITHIN RANS AND SENSITIZED-RANS
FRAMEWORKS

Maximilian Bopp, Sebastian Wegt, Louis Krüger, Suad JakirlićInstitute of Fluid Mechanics and Aerodynamics, Technical University of Darmstadt
Alarich-Weiss-Str. 10, 64287 Darmstadt, Germany
bopp@sla.tu-darmstadt.de, jakirlic@sla.tu-darmstadt.de

Francesco Secchi, Bettina Frohnapfel
Institute of Fluid Mechanics, Karlsruhe Institute of Technology (KIT)
Kaiserstrasse 10, 76131 Karlsruhe, Germany
francesco.secchi@kit.edu, bettina.frohnapfel@kit.de

ABSTRACT

mon underlying problem seems to be the variation of the turbulent Prandtl number (Prt ) in the wall proximity, which is not
at all accounted for, e.g., in the most-widely used model formulation for u0j θ 0 , the so-called gradient diffusion approach- representing the straightforward proportionality to the mean
temperature gradient - in conjunction with the Boussinesq approximation:

A comparative assessment of different modeling approaches for the heat flux u0j θ 0 , with the gradient of which
representing the turbulence-related heat source in the temperature transport equation, is in focus of the present computational
study. The convective heat transfer under consideration takes
place within two benchmark configurations: a fully-developed
channel flow and a pipe-discharging jet impinging perpendicularly onto a heated wall. The heat flux models’ performanceis investigated in conjunction with a conventional Reynoldsstress model and its scale-resolving variant employed withina time-accurate RANS (Reynolds-Averaged-Navier-Stokes)
simulation. In the Sensitized-RANS framework employing
the scale-resolving Improved Instability-Sensitive Reynolds
Stress Model (IISRSM) (Jakirlić & Maduta (2015)), the turbulent correlations are employed within the set of relevant
momentum equations as ’sub-scale’ quantities. The scalar
flux models represent the so-called gradient approach with the
diffusion-like coefficient formulated in terms of both eddyviscosity and the full Reynolds stress tensor. The IISRSMrelated ’sub-scale’ turbulent length and time scales constituting
these relationships are determined by utilizing the unresolved
residual part of the entire anisotropic sub-scale Reynoldsstress tensor, as well as the corresponding turbulence kinetic
energy and its viscous dissipation rate.

−u0j θ = αt

∂Θ
∂xj

with

αt =

νt
Prt

(1)

Here, the eddy diffusivity of the temperature field αt is expressed commonly in terms of the eddy viscosity νt and theconstant turbulent Prandtl number Prt . This kind of formulation leads actually to another major weakness: a simple formulation like this is not capable of capturing turbulence anisotropy residing in modeled diffusive fluxes due tothe scalar nature of the eddy viscosity. Accordingly, not the
Reynolds stress component being mostly affected by the wall
proximity is going to be applied within this proportionality coefficient, but the entire turbulence kinetic energy (k = 0.5ui ui ).
The present work focuses therefore on studying various algebraic expressions for the heat flux correlation −u0j θ 0 aiming
at a comparative assessment of their predictive performanceswith respect to the resulting near-wall temperature field and associated Nusselt number distributions. This background modelof turbulence is a near-wall Reynolds stress model appliedin the context of both computational frameworks, the conventional RANS and the so-called sensitized, eddy-resolving
RANS approach (Jakirlić & Maduta (2015)).

INTRODUCTIONA crucial disadvantage of heat transfer predictions whenusing Reynolds-Averaged Navier-Stokes (RANS) models istheir lack of accuracy in even weakly complex flow configurations. Even in conjunction with the models accounting forthe near-wall turbulence, the majority of the wall heat transfer
models fails to capture correctly the heat transport in the diffusive sub-layer and the thermal buffer layer, Šarić & Basara(2018). The influence of modeled heat transport in RANS simulations is therefore of great importance. Even the heat transfer
studies in simpler flow configurations with exclusively gaseous
media by using near-wall turbulence models are often reported
to deliver wrong Nusselt number predictions compared to corresponding experimental data, Chang & Shyu (2000). A com-

COMPUTATIONAL MODEL
The equations of motion governing the underlying incompressible velocity field within an unsteady conventional/sensitized RANS framework read:
 


∂
∂ (Ui )∂∂Ui
(Ui ) +U j=ν− u0i u0j
∂t
∂xj
∂xj
∂xj

1

(2)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
The Reynolds-stress tensor u0i u0j , mimicking the subscalestress tensor in this RANS-based eddy-resolving method, isdetermined by a second-moment closure formulation by solving the following transport equation:Du0i u0j
Dt= Pi j + Φi j − εihj +0
∂
0.5Dνi j + Dui j∂ xkoutcome in even simple flow configurations. The absence ofan explicit dependence on the Reynolds stresses u0i u0j makesthese models incapable of dealing with turbulence anisotropy,
which is especially pronounced in the wall vicinity. Accordingly, the isotropic nature of the turbulent diffusivity implies
obvious alignment of the turbulent scalar flux and the meanscalar gradient. This is clearly not the case in any complex
flow configuration. A more extensive closure formulation is in
focus of the present modeling activities, representing a generalization towards the gradient–transport hypothesis by using aturbulent diffusivity tensor Di j (see e.g. Abe & Suga (2001)):(3)

The source terms on the RHS of the Reynolds-stress transport equation describe the exactly treated production rate Pi j ,the redistribution process Φi j , the viscous dissipation correla(ν+u0 )tion εi j and various diffusion processes Di j
. The molecular diffusion Dνi j , similar as the production rate Pi j , is to be
exactly treated; all other processes have to be adequately modeled. Unlike a conventional RANS Reynolds-stress model, being incapable of describing any spectral dynamics of the turbulence field, the presently utilized IIS-RSM, representing its
scale-resolving version, is appropriately sensitized towards adequate capturing of turbulent fluctuations. The latter feature is
introduced along with the SAS (Scale-Adaptive Simulation)
methodology, proposed by Menter & Egorov (2010). The IISRSM model accounts for an additional production term whichis introduced into the equation governing the specific homogeneous dissipation rate ω h = ε h /k, with ε h = ε − 0.5Dνk andk = 0.5ui ui .


Dω h
Dt




=

IIS−RSMDω hDt−u0i θ 0 = Di j

Di j = kτ Cθ ,1∂Θ∂xj

u0i u0jk+Cθ ,2u0i u0k u0k u0j
k2

!
(8)

with τ = k/ε h representing the turbulent time scale. In thecase of Cθ ,2 = 0 the model reduces to the well-known GeneralGradient Diffusion Hypothesis (GGDH). The above model isan anisotropy-reflecting algebraic model capable of dealing
with complex flow straining.

+ PIIS−RSM

(4)

RSM

CASE DESCRIPTION AND COMPUTATIONALDETAILSν ∂Θ
− u0j θ 0
Pr ∂ x j

The presently considered flow configurations are concerned with the convective heat transfer in a fully-developed
channel flow at Reτ = 180 (reference DNS by Horiuti (1992)
and Kasagi et al. (1992)) and a round impinging jet onto asmooth, heated wall at Re = 10000 (reference DNS by Secchi
et al. (2022)).
The Reynolds stress model equations are implemented into the finite-volume-based open source toolbox
OpenFOAM R (Open Source Field Operation and Manipulation) with which all simulations are performed. Both the conventional Reynolds-stress model and its instability-sensitizedversion (Jakirlić & Maduta (2015), see also Jakirlić & Maduta(2016)) are applied in the time-accurate RANS framework.
The background model formulation is based on the so-called
homogeneous dissipation rate concept, proposed by Jakirlić &Hanjalić (2002). A controlled adaptive time step ensures a
Courant number smaller than one in the entire solution domain. The discretization of the convective terms in the momentum and temperature equations is performed by applying
the 2nd order accurate Central Differencing Scheme (CDS).
The flow domain is meshed by using the OpenFOAM’s meshutility denoted as ’blockMesh’. The mesh is generated with
an appropriate grading towards the walls, providing the wallnext computational node situated well within the viscous sublayer, exhibiting a dimensionless wall distance y+ substantially smaller than one. This is especially important for thecorrect calculation of the sub-scale heat flux in the immediate
wall proximity.
Fig. 1 illustrates schematically the convective heat transfer within a fully-developed channel flow. The uniform temperature profile over the channel cross-section is introduced
into the flow and subjected to continuous heat removal throughthe walls by a constant heat flux, leading to a constant temperature value at both walls. The solution domain adopted for the(5)

The turbulent heat flux −u0j θ 0 is commonly modeled by
the so-called SGDH formulation (Simple Gradient Diffusion
Hypothesis). In models of this type, the scalar fluxes are related to the mean scalar-field gradients via a scalar turbulent
diffusivity (Γt ).

−u0j θ 0 = Γt

(7)

with

The inherence of the model is its self-adjustment to thescales residing in the unresolved residual motion, by interplaying with the underlying grid resolution, Menter & Egorov(2010). Herewith, the development of turbulent fluctuations
is enabled. Unlike the initial SAS proposal employing the
von Karman length scale (LvK ∝ ∇Ui /∇2Ui ) as a triggeringparameter to enable the resolving mode, the presently utilized
term is modeled only in terms of the second velocity derivative(∇2Ui ), as originated from the equation governing the integral
length scale, Rotta (1972). By doing so the sensitivity of themodel in capturing the turbulence fluctuations is enhanced enabling the use of coarser grid resolutions.
The temperature equation used resembles the most common form of the energy equation for an ideal gas with constant heat capacities, thermal conductivity and density, small
temperature variations and without consideration of thermalradiation (Chang & Shyu (2000)):
∂ (Θ)
∂
∂ (Θ)
+U j
=
∂t
∂xj
∂xj

∂Θ
∂xj

with Γt =Cµ k2
νtor Γt =
Prt
Prt ε

(6)

As it is well-known the models of this type, commonly
applied in terms of a constant Prandtl number, result in a poor2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Sensitized RANS simulations with the IIS-RSM corresponds
typically to a regular hexahedron, confined by bottom and upper walls, with the dimensions Lx × Ly × Lz = 2πh × 2h × πh
and h representing the channel half height. The relevant numerical mesh consists of Nx × Ny × Nz = 60 × 64 × 32 cells,
resulting in 122880 cells in total. The conventional RANS
computations are performed on a two-dimensional flow domain (Lx × Ly = 2πh × 2h) consisting of Nx × Ny = 60 × 64
cells. The height of the wall-next grid cell is ∆y+ ≈ 0.9. The
resolution of the uniformly distributed grid cells in streamwise and spanwise directions corresponds to ∆x+ = 19.0 and∆z+ = 18.0 respectively. The fully-developed flow conditions
are provided by applying periodic inlet-outlet boundary conditions.

Figure 2: Schematic of the impinging jet configuration

Figure 1: Schematic of the fully-developed channel flow

Fig. 2 schematically shows the flow domain adopted forthe simulations of the round impinging jet. The jet structure corresponds to a fully-developed pipe flow at a diameter D-based Reynolds number of Re = 10000. The pipe-likeround jet entering the flow domain is realized by a separate
computation of a pipe segment by applying the periodic inlet/outlet conditions. The flow field thus generated is subsequently mapped onto the domain inlet. The jet generated inthis way impinges perpendicularly on a partially heated wall
situated at 2D distance from the pipe exit, with boundary conditions corresponding to a constant wall temperature. The solution domain having the form of a truncated cylinder, confined in the normal direction by the bottom wall and an upper confinement plate, extends to 12D in the radial direction.
The grid detail in the radial x − y plane is shown in Fig. 3.
The number of the cells within the horizontal x − y plane of
Nx−y = 46400 multiplied by the Nz = 138 in the normal direction amounts 6.4 Million cells in total. The height of the walladjacent grid cell is ∆z+ = 0.06 − 0.94 (the lower values relate
to the immediate impinging area and the higher values to theregion close to the configuration outlet). The grid resolution in
the radial direction corresponds to ∆r+ = 4.76 − 139.54. The
same three-dimensional solution domain and grid resolution
have been applied for both conventional and Sensitized RANS
simulations.

Figure 3: Detail of the numerical grid in the radial x − y plane
for the impining jet configurationby the IISRSM. In both cases the temperature difference between the working fluid and the heated/cooled wall is 10oC.
The temperature field is scaled to have zero values (temperature minimum) at both channel walls and the pipe exit and unit
values (temperature maximum) at the channel symmetry plane
and at the impinging wall. The scale-resolving capabilities of
the IISRSM are clearly evident at both flow fields exhibiting
turbulent fluctuations. The dynamics of both velocity and temperature fields is mainly driven by their resolved fractions generated via the convective terms in the corresponding transport
equations. However, in the wall proximity, the sub-scale turbulent momentum and heat fluxes are adequately enhanced andare therefore crucial for capturing correctly the near-wall effects. In the following individual sections some time-averaged
mean flow and turbulence properties will be discussed alongwith the reference DNS results.

Channel flow
Figures 6 illustrate the mean velocity and Reynolds
stress component profiles obtained by both RANS-RSM andIISRSM across the fully-developed channel flow. Apart froma slight underprediction of the peak value of the streamwise
stress component u2 and a certain overprediction of the spanwise stress component w2 , the results obtained show very good

RESULTS AND DISCUSSIONAs an illustration, Figures 4 and 5 display the instantaneous velocity and temperature fields in both considered configurations, plane channel and round impinging jet, obtained3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

Figure 4: Instantaneous velocity and temperature fields in aturbulent channel flow with fixed wall temperature
obtained by the IISRSM

Figure 6: Mean velocity and Reynolds stress component profiles in the fully-developed channel flow

Figure 5: Flow visualization by the Q-criterion, colored by the
instantaneous temperature field in a turbulent impinging jet with fixed wall temperature obtained by
the IISRSM
Figure 7: Mean temperature profile in the fully-developed
channel flow

agreement with the DNS database in all characteristic boundary layer zones, viscous sublayer, buffer zone and the logarithmic region. Considering a passive character of the temperature field, there is no relevant influence on the mean velocity
and turbulence field. Accordingly, the use of different heat
flux models is irrelevant here. Similar good agreement with
the DNS reference is reported with respect to the mean temperature profile, Fig. 7. The influence of the SGDH (Eq. 6)
and GGDH (Eqs. 7-8, with Cθ ,1 = 0.3 and Cθ ,2 = 0) heat
flux model formulations on the temperature profile development is not noticeable. Considering that both heat flux formulations were applied in conjunction with a full Reynolds stress
RANS model, for which a correctly predicted Reynolds stressfield representing the constituent of the relevant gradient diffusion coefficients was used, this outcome is not very surprising.
On the other hand, this was anyhow expected for the IISRSM
application, given the appropriate level of the resolved turbulence. The modeled heat flux relates here only to the residual
sub-scale motion. Accordingly, only the results of the GGDH
formulation is shown in conjunction with the eddy-resolving
RSM version.

impinging jet configurations featured by an intensive alternation in the velocity gradient, from abrupt jet deceleration in
the immediate impingement area toward a prompt 90o flow
skewing and transition into the wall-jet region followed by astrong acceleration and subsequent flow relaxation. It is furthermore well known that a streamline curvature-induced turbulence production takes place here, which represents an inherent feature of Reynolds stress modeling group, but on theother hand is beyond the reach of the eddy-viscosity modelingconcept. Due to the upper confinement plane (see Fig. 2), a
flow reversal of weak intensity is generated above the wall jet
region.
Figures 10 illustrate the computationally, by IISRSM and
RSM, obtained semi-log profile development of the mean velocity magnitude, Reynolds stress components and temperature at selected locations within the flow domain along with
the DNS reference results. The vertical coordinate is logarithmically scaled to emphasize the near-wall behavior of the flow
quantities considered. Herewith, a detailed quantitative assessment of the results quality is enabled. The velocity profiles
are characterized by a flow momentum intensification occurring during the impingement event (at r/D = 0.5 − 1.0) andsubsequent transition (at r/D ≥ 1.0) into the wall jet. After-

Impinging jet
Figures 9 depict the iso-contours of the time-averaged
velocity field including mean flow streamlines. The mean
flow topology reveals a flow pattern typically characterizing

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
in the thermal field results with respect to the application of
the SGDH and GGDH models were not revealed this time,
because the background Reynolds stress model of turbulence
used in conjunction with both heat flux model expressions,provided a qualitatively valuable coefficient prediction in both
gradient-diffusion approaches tested.

ward the flow momentum weakens appropriately. Agreementwith the DNS database is remarkable, especially concerning
the IISRSM results. The RSM-related results shows also a
reasonable level of agreement with DNS. It can be said, considering the large departures of the RSM-specific Reynolds
stress components, that the latter behavior is mainly due to therole of the mean pressure gradient which dominates the balance in the equation of motion, overcoming significantly the
Reynolds stress gradients. The radial (parallel to the impinging plate) ur ur and shear ur uz stress components returned by
the IISRSM follow closely the DNS database. At the crosssection r/D = 0.5, which coincides with the pipe outlet, themaxima of both stress components originate from the nearwall region within the pipe. The correctly predicted double
maximum at the ur ur component, visible at r/D = 1.0 − 2.0,
coincides with the wall-jet boundary characterized by a velocity maximum at z/D = 0.1 and the strong mean shear in the immediate vicinity of the wall at z/D = 0.01. The RSM Reynolds
stress profiles reveal significant overprediction in the direct impingement area at r/D = 0.5 and 1.0. This can be partiallyexplained by the ’standard’ expression of the so-called wallreflection redistribution term that features the present RSM,
which was originally formulated for wall-parallel flows (seeCraft et al. (1993) for a related discussion). The difference inthe performance of two tested heat flux models, SGDH and
GGDH, is visible at the temperature profiles, although a clear
superiority of any of the applied models cannot be clearly established. Small differences in the wall temperature gradientsare reflected in the distribution of the Nusselt number displayed in Fig. 8. A characteristic second peak at r/D ≈ 2,apart of that related to the impingement center (at r/D = 0),
is returned by all models despite an obviously flattened Nudistribution at this position. The RSM-related enhanced peakat r/D = 0 is due to increased turbulence intensity in the immediate impingement location. The IISRSM-specific underprediction requires a more detailed analysis also with respect
to the grid resolution requirements, considering that the standard grid resolution constraints are primarily derived for thewall-parallel flows (see e.g., Hadžiabdć & Hanjalić (2008)).

Acknowledgments
The financial support of the German Research Foundation (DFG) in the framework of the Collaborative Research Center/Transregio 150 (SFB-TRR 150, Project Number
237267381) is gratefully acknowledged. The authors furthermore would like to thank for the computing time granted on
the Lichtenberg HPC at TU Darmstadt.

REFERENCES
Abe, K. & Suga, K. 2001 Towards the development of a
Reynolds-averaged algebraic turbulent scalar-flux model.
Int. J. Heat and Fluid Flow 22, 19–29.
Chang, K.-C. & Shyu, M.-J. 2000 Revisiting the Reynoldsaveraged energy equation in near-wall turbulence models.
Int. J. Heat and Mass Transfer 43, 665–676.
Craft, T.J., Graham, L.J.W. & Launder, B.E. 1993 Impinging
jet studies for trbulence model assessment - ii. an examination of the performance of four turbulence models. Int. J.
Heat and Mass Transfer 36(10), 2685–2697.
Hadžiabdć, M. & Hanjalić, K. 2008 Vortical structures andheat transfer in a round impinging jet. Journal of Fluid Mechanics 596, 221–260.
Horiuti, K. 1992 Assessment of two-equation models of turbulent passive-scalar diffusion in channel flow. Journal of
Fluid Mechanics 238, 405–433.
Jakirlić, S. & Hanjalić, K. 2002 A new approach to modellingnear-wall turbulence energy and stress dissipation. Journalof Fluid Mechanics 459, 139–166.
Jakirlić, S. & Maduta, R. 2015 Extending the bounds of‘steady’ RANS closures: Toward an instability-sensitive
Reynolds stress model. Int. J. Heat and Fluid Flow 51, 175–
194.
Jakirlić, S. & Maduta, R. 2016 Sensitized-RANS modelling of
turbulence: Resolving turbulence unsteadiness by a (nearwall) Reynolds stress model. In Progress in Wall Turbulence 2 (ed. M. Stanislas et al.), pp. 17–35. Springer.
Kasagi, N., Horiuti, K., Miyake, Y., Miyauchi, T. &Nagano, Y. 1992 Establishment of the direct numerical simulation data bases of turbulent transport phenomena. Co-operative Research No. 02302043. Supportedby the Ministry of Education, Science and Culture,https://thtlab.jp/DNS/CH122 PG.WL3.
Menter, F.R. & Egorov, Y. 2010 The scale-adaptive simulation method for unsteady turbulent flow predictions. part 1:
theory and model description. Flow, Turbulence and Combustion 85 (1), 113–138.
Rotta, J.C. 1972 Turbulente Strömungen. Eine Einführung in
die Theorie und ihre Anwendung, 1st edn. Teubner Verlag
Stuttgart.
Secchi, F., Häber, T., Gatti, D., Schulz, S., Trimis, D., Suntz,
R. & Frohnapfel, B. 2022 Turbulent impinging jets on roughsurfaces. GAMM, Special Issue on Direct Numerical Simulations of Turbulent Flows - Part I 45.
Šarić, S. & Basara, B. 2018 A near-wall model for heat transfer at high Prandtl numbers. In International Conference
on Computational Fluid Dynamics (ICCFD10), Barcelona,
Spain, July 9-13.

Figure 8: Nusselt number distribution at the impinging wall

CONCLUSIONSThe predictive performance of heat flux (u0j θ 0 ) modelingin conjunction with passive scalar transport in a plane channel
and an impinging jet flow configuration was studied presently.
The results obtained by applying the RANS Reynolds stressmodel and its eddy-resolving variant denoted by IISRSM showa high degree of consistency. The turbulence level overprediction pertinent to the central impingement region obtained by
the RSM could be significantly improved by utilizing the corresponding model version with the enhanced capability of accounting for turbulence unsteadiness. Pronounced differences

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

Figure 9: Iso-contours of the time-averaged velocity field with associated mean flow streamlines in the impinging jet configuration
obtained by the IISRSM

Figure 10: Profile development of the mean velocity magnitude, radial and shear stress components and mean temperature in theimpinging jet configuration

6