# Physics-constrained machine learning for thermal turbulence modelling at low
Prandtl numbers

Matilde Fiore  Lilla Koloszar  Clyde Fare  Miguel Alfonso Mendez  Matthieu
Duponcheel  Yann Bartosiewicz  Von Karman Institute for Fluid Dynamics, Sint-
Genesius-Rode 1640, Belgium  IBM Research Europe, The Hartree Centre, United
Kingdom  Institute of Mechanics, Materials and Civil Engineering (IMMC),
Université catholique de Louvain (UCLouvain), Place du Levant 2, 1348 Louvain-
la-Neuve, Belgium

###### Abstract

Liquid metals play a central role in new generation liquid metal cooled
nuclear reactors, for which numerical investigations require the use of
appropriate thermal turbulence models for low Prandtl number fluids. Given the
limitations of traditional modelling approaches and the increasing
availability of high-fidelity data for this class of fluids, we propose a
Machine Learning strategy for the modelling of the turbulent heat flux. A
comprehensive algebraic mathematical structure is derived and physical
constraints are imposed to ensure attractive properties promoting
applicability, robustness and stability. The closure coefficients of the model
are predicted by an Artificial Neural Network (ANN) which is trained with DNS
data at different Prandtl numbers. The validity of the approach was verified
through a priori and a posteriori validation for two and three-dimensional
liquid metal flows. The model provides a complete vectorial representation of
the turbulent heat flux and the predictions fit the DNS data in a wide range
of Prandtl numbers (Pr=0.01-0.71). The comparison with other existing thermal
models shows that the methodology is very promising.

###### keywords:

Low-Prandtl , Turbulence , Heat Flux , Artificial Neural Networks

††journal: Journal Name

##  1 Introduction

Heat transfer in liquid metals is fundamental in many engineering
applications, including high performance cooling systems in new generation
nuclear reactors. For such systems, the thermal hydraulics conditions
establishing during normal operation and transients must be carefully analyzed
for economic and safety reasons. Specifically, the safe design of liquid metal
cooled nuclear reactors requires the accurate description of the temperature
distributions in both solid and fluid regions to evaluate the efficiency of
the cooling system and prevent phenomena of solidification or thermal fatigue
[1].

Experimental investigations of liquid metal flows are particularly challenging
due to their operating conditions, the metal opaqueness and its chemical
aggressiveness [2, 3]. Therefore, numerical investigations play a dominant
role during design and testing to provide accurate descriptions of the
thermal-hydraulics conditions with reasonable computational costs. These
requirements often collide with the limited availability of RANS thermal
models in commercial codes: the Reynolds analogy [4] is usually provided to
model the turbulent heat flux. This simple model is based on a simple gradient
hypothesis, and on the definition of a constant turbulent Prandtl
number111Recommended values for the turbulent Prandtl numbers are 0.85 for
near unity Prandtl numbers and around 2.0 for low Prandtl number fluids [5]..
The approach is widely accepted in presence of forced convection of common
fluids, as water or air, having near unity or higher than unity Prandtl
numbers. However, at the low Prandtl numbers characterizing liquid metals, the
similarity hypothesis between momentum and thermal turbulence is less
justified. Low Prandtl number fluids then require more sophisticated modelling
that should go beyond the assumptions of similarity and isotropy.

Researchers developed a variety of thermal turbulence models for low Prandtl
numbers [6, 7, 8, 9], which gave satisfactory predictions in several flows of
academic interest. However, the range of Prandtl numbers employed for their
validation is too narrow to expect that these models generalise well in
industrial scenarios. In addition, none of the formulations accurately matches
all the components of the turbulent heat flux, especially those which are
normal to the mean temperature gradient, providing an incomplete
representation of thermal turbulence [10].

Clearly, the increasing amount of DNS data for many Prandtl numbers and flow
configurations represents a valuable source of information that should be
efficiently used to construct more accurate models. The information should be
extracted from data, summarised and translated in a mathematical form. Recent
developments in Machine Learning provide structured regression methods whose
use is emerging in the context of turbulence modelling. Data-driven methods
are employed to construct advanced representations of the turbulence
statistics [11, 12, 13, 14, 15, 16] or to introduce corrections for the
production/destruction terms appearing in the transport equations [17, 18,
19]. The input-output mappings derived with these techniques are essentially
local and algebraic, in order to be mesh independent.

Among the available techniques, artificial neural networks (ANNs) have become
very popular to solve high dimensional and highly non-linear regression
problems which are common in turbulence modelling. Compared to traditional
approaches, the data-driven approach relies on data and, as such, it opens
enormous challenges related to the generalization and the practical
applicability of the developed formulation. If not properly trained and
regularized, data-driven models can lead to unexpected or nonphysical results
when applied far from the conditions involved in the training. The model
robustness needs to be increased by embedding essential physical and
mathematical properties indirectly or by construction. Examples of essential
properties are the invariance under rotation of the coordinate system, the
realizability of the Reynolds stresses, or the consistency with the second law
of thermodynamics. Moreover, the predicted fields must have smooth
derivatives, especially in the vicinity of the walls, to ensure the stable
integration of RANS transport equations. In the special case of a thermal
turbulence model, the coupling with its momentum counterpart must be
considered as well, to develop a robust formulation that is not too sensitive
to inaccuracies of the momentum modelling.

Low Prandtl number fluid flows pose additional challenges to the data-driven
approach for the heat flux modelling. First of all, at low Prandtl the thermal
diffusion effects are significant and the locality hypothesis, for which the
heat flux at a given point is a function of local flow statistics, is less
justified. This means that an accurate algebraic closure valid in the low
Prandtl number range is difficult to develop, as non-local effects should be
reproduced. Moreover, the near wall behaviour of the thermal turbulence
statistics is highly sensitive to the thermal boundary conditions at low
Prandtl numbers [20], which have instead minimal influence on the mean
temperature distribution. Hence, there is the need of including some thermal
turbulence statistics in the parametrization to achieve a good fit. This
choice would tend to make the model non-linear with respect to the thermal
quantities, which is against the linearity principle expressed by Pope [21]
for passive scalar turbulence.

In this work, we present a deep-learning approach for the modelling of the
turbulent heat flux that embeds essential physical and mathematical properties
and aims at being valid in a wide range of Prandtl numbers and flow
configurations. Section 2 presents the mathematical structure of the model,
the details of the artificial neural network and its training, while section 3
briefly introduces existing thermal low Prandtl based turbulence models used
for comparison. The database employed to train the model is described in
section 4, while section 5 presents the auxiliary models applied for the
validation of the data-driven closure. Section 6 outlines the results of the
ANN training, the a priori and a posteriori validation and the limitations of
the current formulation. Such limitations are then explained in detail in
section 7, where a sensitivity analysis is proposed to attribute the
discrepancy to the most critical parameters. Finally, conclusions and outlooks
are provided in the last section.

##  2 Mathematical formulation

A model for the turbulent heat flux $\overline{\mathbf{u}\theta}$ is required
to close the averaged energy equation for incompressible turbulent flows:

$\begin{split}\frac{\partial T}{\partial{t}}+U_{j}\frac{\partial
T}{\partial{x_{j}}}=\alpha\frac{\partial^{2}T}{\partial{x_{j}}\partial
x_{j}}-\frac{\partial\overline{u_{j}\theta}}{\partial{x_{j}}}.\end{split}$ (1)

having used Einsten’s notation of implied summation, denoting as $U_{j}$ the
$j$-th component of the mean velocity and referring the reader to the list of
symbols for the remaining symbols. The data-driven approach for
$\overline{\mathbf{u}\theta}$ is based on the development of a general
mathematical structure that is presented in section 2.1. The closure
coefficients of the algebraic expression are then modelled with an artificial
neural network, whose structure is defined in section 2.5. The regression
problem consists then in finding the network parameters that lead the model
fitting the available DNS data presented in section 4. The optimization of the
network parameters is guided by an appropriate objective function that is
given in section 2.4.

###  2.1 General assumptions

The algebraic model for the turbulent heat flux $\overline{\mathbf{u}\theta}$
was derived starting from a general functional relationship involving tensors,
vectors and scalars of interest:

$\overline{\mathbf{u}\theta}=f\left(\mathbf{b},\mathbf{S},\mathbf{\Omega},\nabla
T,\mathbf{g},k,\epsilon,k_{\theta},\epsilon_{\theta},\alpha,\nu\right),$ (2)

where tensors $\mathbf{b}$, $\mathbf{S}$ and $\mathbf{\Omega}$ are defined as:

$b_{ij}=\frac{\overline{u_{i}u_{j}}}{k}-\frac{2}{3}\delta_{ij},$ (3)

$S_{ij}=\frac{1}{2}\left(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial
U_{j}}{\partial x_{i}}\right),$ (4)

$\Omega_{ij}=\frac{1}{2}\left(\frac{\partial U_{i}}{\partial
x_{j}}-\frac{\partial U_{j}}{\partial x_{i}}\right).$ (5)

The ansatz (2) implicitly assumes locality and instantaneity as a direct
consequence of the algebraic structure. From a numerical point of view, this
assumption is convenient, as it avoids the expenses of a second order closure
model for the heat flux. On the other hand, if transport equations are solved
for $k_{\theta}$ and $\epsilon_{\theta}$, historical and non-local effects can
be still represented, though at an isotropic level. Equilibrium turbulence is
also assumed, since the gradients of the turbulent statistics ($k$,
$k_{\theta}$ and $b_{ij}$) are not involved in the parametrization. This
requires that thermal turbulence production, destruction and redistribution
mechanisms are more significant than transport terms, as it is assumed in many
other works [22, 23, 24, 25] on this topic222Clearly, this assumption is not
satisfied in the near wall region, and asks the algebraic model to
artificially mimic transport effects through the other parameters involved in
the parametrization.. In addition, if the analysis is limited to forced
convection flows, the dependence on the gravity vector can be eliminated.

Since the model is algebraic and explicit, it could be also written as a
product of a dispersion tensor $\mathbf{D}$ with the mean temperature
gradient:

$\overline{\mathbf{u}\theta}=-\mathbf{D}\nabla T,$ (6)

where

$\mathbf{D}=\mathcal{F}\left(\mathbf{b},\mathbf{S},\mathbf{\Omega},\nabla
T,k,\epsilon,k_{\theta},\epsilon_{\theta},\alpha,\nu\right).$ (7)

For the sake of simplicity, we can assume that the anisotropic part of the
dispersion tensor $\mathbf{D}$ depends on the anisotropy of the momentum
field, not on the temperature distribution. Hence, it is assumed that the
temperature gradient does not affect the anisotropic part of $\mathbf{D}$.
Then, the functional relationship simplifies to:

$\mathbf{D}=\mathcal{F}\left(\mathbf{b},\mathbf{S},\mathbf{\Omega},||\nabla
T||,k,\epsilon,k_{\theta},\epsilon_{\theta},\alpha,\nu\right).$ (8)

having limited the attention to the invariants of $\nabla T$, i.e. its norm
$||\nabla T||$.

###  2.2 Invariance Properties

A model for the turbulent heat flux must be invariant under rotation of the
coordinate system. This property can be analytically expressed as:

$\mathbf{Q}\mathbf{D}\mathbf{Q}^{T}=\mathcal{F}\Bigl{(}\mathbf{Q}\mathbf{b}\mathbf{Q}^{T},\mathbf{Q}\mathbf{S}\mathbf{Q}^{T},\mathbf{Q}\mathbf{\Omega}\mathbf{Q}^{T},||\nabla
T||,k,\epsilon,k_{\theta},\epsilon_{\theta},\alpha,\nu\Bigr{)},$ (9)

where $\mathbf{Q}$ is any rotation matrix. Hence, the compliance with the
invariance property requires constraining the functional form in eq. (8).

Let us consider a general set of tensors $\mathbf{T}^{i}$, formed through
products of $\mathbf{b}$, $\mathbf{S}$ and $\mathbf{\Omega}$. A linear mapping
from $\mathbf{D}$ to $\mathbf{T}^{i}$ would satisfy the rotational invariance:

$\mathbf{Q}\mathbf{D}\mathbf{Q}^{T}=\sum_{i=1}^{n}c_{i}\mathbf{Q}\mathbf{T}^{i}\mathbf{Q}^{T},$
(10)

if the coefficients $c_{1},...,c_{n}$ are isotropic functions, i.e. functions
of the invariants of the tensorial set described by eq.(8). Clearly,
infinitely many tensor products can be formed from tensors $\mathbf{b}$,
$\mathbf{S}$ and $\mathbf{\Omega}$, thus leading to an infinite series
expansion and, accordingly, an infinite number of tensorial invariants.

However, tensor representation theory shows that, for a finite tensorial set,
there exists a finite integrity basis of independent invariants under the
proper orthogonal group [26]. Given the minimum tensor basis of a certain
degree, higher order tensor products can be expressed as polynomials of the
basis tensors, whose coefficients depend on the invariant basis (see also
[27]). This property bounds the parameter space and make the regression
problem more tractable.

The complete derivation of the tensor basis and the corresponding invariants
under the assumption of 2-D flow is given in Appendix Appendix. Briefly, it
can be shown that the minimal basis of tensors consists of the following set:

$\begin{split}\mathbf{I},\mathbf{b},\mathbf{S},\mathbf{\Omega},\mathbf{bS},\mathbf{S\Omega},\mathbf{b\Omega},\mathbf{b\Omega
S},\end{split}$ (11)

and the corresponding minimal basis of invariants is given by

$\\{\mathbf{b_{2}}\\},\\{\mathbf{b}^{2}\\},\\{\mathbf{S}^{2}\\},\\{\mathbf{\Omega}^{2}\\},\\{\mathbf{bS}\\},\\{\mathbf{bS\Omega}\\},$

where the notation $\\{\dot{\\}}$ indicates the tensor trace and
$\mathbf{b}_{2}$ denotes the projection of $\mathbf{b}$ on the $x-y$ plane of
the flow, i.e.

$\mathbf{b}_{2}=\begin{pmatrix}b_{11}&b_{12}&0\\\ b_{12}&b_{22}&0\\\
0&0&0\end{pmatrix}\,.$ (12)

The coefficients $c_{i}$ in eq.(10) will then be functions of the following
isotropic quantities:

$\begin{split}c_{i}=f_{i}\Bigl{(}\\{\mathbf{b_{2}}\\},\\{\mathbf{b}^{2}\\},\\{\mathbf{S}^{2}\\},\\{\mathbf{\Omega}^{2}\\},\\{\mathbf{bS}\\},\\{\mathbf{bS\Omega}\\},||\nabla
T||,k,\epsilon,k_{\theta},\epsilon_{\theta},\alpha,\nu\Bigr{)}.\end{split}$
(13)

The number of parameters involved in eq.(13) was further reduced by appylying
the Buckingham Theorem to the group of physical scalar quantities. The
approach leads to the definition of 10 independent dimensionless groups, which
are indicated in the first column of Table 1. Then, the coefficients $c_{i}$
are expressed as

$c_{i}=f_{i}(\pi_{i},Re_{t},Pr).$ (14)

The reader should note that the dimensionless groups contain temperature
dependent quantities such as $||\nabla T||$, $k_{\theta}$ and
$\epsilon_{\theta}$. In particular, the parameter $\pi_{7}$ is view as
critical, as the dependency of the coefficients $c_{i}$ on $\pi_{7}$ will
alter the linear relationship between $\overline{\mathbf{u}\theta}$ and
$\nabla T$, at the basis of the most of the simple gradient and generalized
gradient assumptions.

###  2.3 Consistency with the second law of thermodynamics

A thermal turbulence model is realizable if it satisfies the second law of
thermodynamics, which requires that the real part of the eigenvalues of
$\mathbf{D}$ are positive [28]. This property prevents counter gradient heat
fluxes in forced convection regimes. To satisfy this requirement, $\mathbf{D}$
is expressed as:

$\mathbf{D}=\left[(\mathbf{A}+\mathbf{A}^{T})(\mathbf{A}^{T}+\mathbf{A})+\frac{k}{\epsilon^{0.5}}(\mathbf{W}-\mathbf{W}^{T})\right],$
(15)

where $\mathbf{A}$ and $\mathbf{W}$ can be written as combination of the basis
tensors $\mathbf{T}^{i}$:

$\displaystyle\mathbf{A}$ $\displaystyle=\sum_{i=1}^{n}a_{i}\mathbf{T}^{i},$
$\displaystyle\mathbf{W}$ $\displaystyle=\sum_{i=1}^{n}w_{i}\mathbf{T}^{i},$
(16)

and the coefficients $a_{i}$ and $w_{i}$ have the functional form in eq. (14).
The formulation is then rotational invariant according to the discussion in
the previous section. Eq.(15) enforces a Cholesky decomposition for the
symmetric part of $\mathbf{D}$, to enforce the real part of its eigenvalues to
be positive (see Roger et al. [29]). On the other hand, the skew-symmetric
part has imaginary eigenvalues, which do not affect the net energy exchange.

It is worth noting that the invariants $\pi_{i}$ indicated in the first column
of Table 1 are expressed in dimensionless form to make the data-driven
inference valid for arbitrary flow conditions. Accordingly, the modelled
coefficients $a_{i}$ and $w_{i}$ will be dimensionless. Hence, the tensors
constituting the basis in eq.(11) are multiplied by powers of $k$ and
$\epsilon$ to get a dimensional heat flux according to333Note that the choice
of the powers of $k$ and $\epsilon$ is not restrictive for the
parametrization, since each tensor can be multiplied by any function of the
invariants indicated in Table 1. eq.(6), (15) and (16). The resulting
definition of the tensors $\mathbf{T}^{i}$ is given in the second column of
Table 1.

Table 1: Basis tensors and invariants. Invariant Basis |  Tensor Basis
---|---
$\pi_{1}=\frac{k^{2}}{\epsilon^{2}}\\{\mathbf{S}^{2}\\}$, $\pi_{2}=\frac{k^{2}}{\epsilon^{2}}\\{\mathbf{\Omega}^{2}\\}$, $\pi_{3}=\frac{1}{\\{\mathbf{b}\\}}$, $\pi_{4}=\frac{k}{\epsilon}\\{\mathbf{bS}\\}$, $\pi_{5}=\\{\mathbf{b}_{2}\\}$, $\pi_{6}=\\{\mathbf{bS\Omega}\\}\frac{k^{2}}{\epsilon^{2}}$, $\pi_{7}=\frac{||\nabla T||}{\sqrt{k_{\theta}}}\frac{k^{3/2}}{\epsilon}$, $\pi_{8}=\frac{k_{\theta}\epsilon}{k\epsilon_{\theta}}$, $Re_{t}=\frac{k^{2}}{\epsilon\nu}$, $Pr=\frac{\nu}{\alpha}$ |  $\mathbf{T_{1}}=\frac{k}{\epsilon^{1/2}}\mathbf{I}$, $\mathbf{T_{2}}=\frac{k}{\epsilon^{1/2}}\mathbf{b}$, $\mathbf{T_{3}}=\frac{k^{2}}{\epsilon^{3/2}}\mathbf{S}$, $\mathbf{T_{4}}=\frac{k^{2}}{\epsilon^{3/2}}\mathbf{\Omega}$, $\mathbf{T_{5}}=\frac{k^{2}}{\epsilon^{3/2}}\mathbf{bS}$, $\mathbf{T_{6}}=\frac{k^{2}}{\epsilon^{3/2}}\mathbf{b\Omega}$ , $\mathbf{T_{7}}=\frac{k^{3}}{\epsilon^{5/2}}\mathbf{S\Omega}$, $\mathbf{T_{8}}=\frac{k^{3}}{\epsilon^{5/2}}\mathbf{bS\Omega}$

###  2.4 Smoothness requirement

The smoothness of the predicted heat flux field is an important requirement,
because it enters the averaged energy equation through its divergence.
Specifically, eq. (1) shows that the smoothness of the thermal field (and the
stability of the numerical integration) relies on the smoothness of the
predicted turbulent heat flux. Based on such considerations, the following
loss function was employed for the training:

$\begin{split}\mathcal{L}=\frac{1}{N}\left(\sum_{i=1}^{N}\sum_{j=1}^{3}\Bigl{(}\hat{q}_{i,j}-q_{i,j}\Bigr{)}^{2}\right)+\frac{\lambda}{N}\left(\sum_{i=1}^{N}\sum_{j,k=1}^{3}\left|\frac{\partial\hat{q}_{i,j}}{\partial
x_{k}}-\frac{\partial q_{i,j}}{\partial x_{k}}\right|\Delta
x_{k}\right),\end{split}$ (17)

where $i\in[1,\dots N]$ is the index spanning across the $N$ data points
contained in each mini-batch, $q_{i,j}=\overline{u_{j}\theta}$ is the
prediction of the ANN and $\hat{q}_{i,j}$ is the corresponding flux provided
by DNS data. Both values $q_{i,j}$ and $\hat{q}_{i,j}$ are normalized with
respect to the maximum DNS value achieved in each flow configuration. The
partial derivatives appearing on the right-hand side are numerically
approximated with finite differences to approximate a mean integral of the
absolute error made in the turbulent heat flux derivatives. The regularizing
parameter $\lambda$ is thus a mesh independent hyperparameter that was
empirically set to 10.0.

###  2.5 Structure of the artificial neural network and training

The coefficients $a_{i}$ and $w_{i}$ in eq.(16) are modelled by an ANN which
takes the input vector $\mathbf{X}$ constituted by the invariants indicated in
Table 1:

$\mathbf{X}=\left[\pi_{1},...,\pi_{8},Re_{t},Pr\right],$ (18)

and gives in output the coefficients $a_{i}$ and $w_{i}$:

$\mathbf{Y}=\left[a_{1},...,a_{8},w_{1},...,w_{8}\right],$ (19)

through a general function $\mathcal{M}$ defined by the network architecture
and its weights $\mathbf{\mathcal{W}}$:

$\mathbf{Y}=\mathcal{M}(\mathbf{X};\mathbf{\mathcal{W}}).$ (20)

The model weights $\mathbf{\mathcal{W}}$ are determined with a stochastic
gradient based optimization of the loss function employing the backpropagation
[30] to compute the gradients.

The structure of the proposed ANN developed in Pytorch444See
https://pytorch.org/ is depicted in Figure 2: the network consists of two
input layers for the invariant basis and the molecular Prandtl number,
respectively. Six hidden layers are set for the first branch of the network,
which has rectified linear unit activation functions. In order to enable the
ANN to leverage on products and polynomials of the invariants $\pi_{i}$, the
inputs are logarithmically transformed and the outputs are passed through an
hyperbolic tangent activation in the last hidden layer:

$\tanh(x)=\frac{\left(\exp{x}-\exp(-x)\right)}{\left(\exp{x}+\exp(-x)\right)}.$
(21)

The second branch of the network consists of two hidden layers, each with a
rectified linear unit and hyperbolic tangent activation functions. The two
branches of the network merge through a multiplication of the outputs
performed in the final layer. Then, the turbulent heat flux is calculated
following eq. (15), (16).

The network is trained using the Adaptive Moment Estimation (ADAM) optimizer,
which is a classic minibatch gradient descend with momentum estimates and a
constant learning rate that here is set to 0.001. In this work, the batches
are constructed by sampling subpartitions of the domains belonging to
different test cases. This allows to compute the output’s spatial derivatives,
as required by the cost function (17).

The network was trained for 12000 epochs with several choices for the number
of units in the hidden layers. Figure 1 shows the evolution of the loss
function during the training evaluated over the training a test datasets. It
can be seen that, with hidden layers of 25 and 50 units, the neural network
tends to underfit both the sets of data. By increasing the number of units
from 100 to 200, the training loss shows a slight decreases, though the
validation loss increases significantly. Hence, the structure involving 100
neurons for each hidden layer was selected as the best compromise to prevent
both underfitting and overfitting.

(a) Training loss

(b) Validation loss

Figure 1: Evolution of the loss function evaluated on training and test
datasets for networks with different number of units in the hidden layers.

###  2.6 Uncertainty quantification

The Gaussian Stochastic Weight Averaging (SWAG) algorithm [31] was implemented
and used during the training to quantify the epistemic uncertainty of the data
driven model. The method relies on the assumption that the model uncertainty
is highly correlated to the trajectory followed by the optimizer in the weight
space. The simplest probabilistic model for the weights is given by a
multivariate Gaussian distribution of mean $\mu$ and covariance $\Gamma$:

$\mathcal{W}\sim\mathcal{N}(\mu,\Gamma),$ (22)

During the training, the approximate Gaussian posterior distribution for the
network weights $\mathcal{W}$ is constructed by collecting first and second
moments of the updated parameters. Herein, the covariance matrix $\Gamma$ is
assumed to be diagonal to simplify the implementation of the SWAG555This
strong assumption implies that the weights are independent on each other.. The
Gaussian approximation is constructed in the last 2000 epochs of the training.
After the training, model averaging is performed by sampling the parameters
from the resulting Gaussian distribution:

$\mathcal{W}=\mu+\frac{1}{\sqrt{2}}\Gamma^{1/2}z_{1}\text{ where
}z_{1}\sim\mathcal{N}(0,I_{d}).$ (23)

The averaging provides the mean prediction and the standard deviation of the
network output, which can be used to construct confidence intervals for the
predicted components of the turbulent heat flux. This information could give
important indications about the effectiveness of the learning process and its
convergence, as well as the nature of the prediction error observed during
testing. It could also highlight possible under-fitting and/or over-fitting
problems, thus suggesting improvements for the model structure and the
training methodology. Further details about the SWAG approach for uncertainty
quantification can be found in [31].

Figure 2: Structure of the artificial neural network used to predict the
turbulent heat flux as a function of the molecular Prandtl number and the
basis of invariants indicated in Table 1.

###  2.7 Noise addition and corrections

The algebraic structure of the model given by eq.(6) is gradient-driven by
definition, and thus it looses its theoretical justification where the heat
flux does not correlate with the temperature gradient. This means that the
regression problem is ill-posed in some vanishing gradient regions where DNS
simulations give finite values of the turbulent heat flux. One example is the
region close to the centerline in case of non-isothermal turbulent channel
flow: as shown in Figure 7, the streamwise component of the heat flux reaches
a non-zero value at the centerline, although the temperature gradient is zero.
The behaviour of the heat flux in this region is a result of the local
equilibrium between redistribution, dissipation and transport effects, which
require a differential model to be properly represented. The ANN would tend to
compensate the flaws of this mathematical structure by increasing the
magnitude of $\mathbf{D}$ in the vanishing gradient regions, to obtain a
better fitting of the heat flux components. This issue has been mentioned by
Sotgiu et al. [32] when developing their data-driven algebraic model for the
Reynolds stresses: the authors fixed it by eliminating the training points for
which the trained scalar functions were higher than a given threshold. In the
framework of this work, two strategies were jointly applied (during and after
the training) to mitigate such an unphysical behaviour:

  * 1.

Noise addition: during the training, $||\nabla T||$ is perturbed by adding
noise from a Gaussian distribution of variance equal to 0.05% of the maximum
temperature gradient observed in each flow configuration:

$\widetilde{||\nabla T||}=||\nabla T||+\mathcal{N}(0,5\cdot
10^{-4})\cdot\max_{x,y,z}||\nabla T||$ (24)

This approach is aimed at reducing the sensitivity of the network to the
parameter $\pi_{7}$, which is critical, since it could alter the linear
dependency of $\overline{\mathbf{u}\theta}$ on its driving force $\nabla T$.

  * 2.

Correction: the symmetric part ($\mathbf{A}$) of the predicted dispersion
tensor ($\mathbf{D}$) is a posteriori corrected to prevent extremely large
values of $\mathbf{D}$:

$\widetilde{A}_{i,j}=\begin{cases}\max\left(A_{ij},20\cdot\alpha_{t}\right),&\text{if
}i=j,\\\
\min\left(\max\left(A_{ij},-\sqrt{\prod_{k}A_{kk}}\right),\sqrt{\prod_{k}A_{kk}}\right),&\text{otherwise}\end{cases}$
(25)

where $\alpha_{t}$ is the turbulent thermal diffusivity calculated as
specified by Manservisi [6] and $\prod_{k}A_{kk}$ denotes the product of the
diagonal entries in $\mathbf{A}$. Note that the correction applied on the main
diagonal is equivalent to limit the minimum turbulent Prandtl number. Under
the assumption of two-dimensional flow, this correction does not alter the
conditions on the eigenvalues required to satisfy the second law of
thermodynamics666In the more general, three-dimensional case, the correction
for the off-diagonal entries of $\mathbf{A}$ should be implemented to satisfy
the conditions $\\{\mathbf{A}\\}^{2}-\\{\mathbf{A}^{2}\\}>0$ and
$\det\mathbf{A}>0$, to ensure that the real part of the eigenvalues of
$\mathbf{D}$ are non-negative..

##  3 Some thermal turbulence models for low Prandtl numbers

Current literature provides thermal turbulence models for low Prandtl number
which go beyond the Reynolds Analogy either by adding terms to the
constitutive relationship of the heat flux $\overline{\mathbf{u}\theta}$ or
modifying the expression of the eddy diffusivity, e.g. introducing thermal or
mixed time scales rather than dynamical ones. An overview of these low-Prandl
based closures is given by Shams et al.[10]. Some of these formulations have
been used in the present work for the validation of the data-driven thermal
model:

  * 1.

The Manservisi and Meneghini model [6] and the Kays correlation [33] are based
on a standard gradient diffusion hypothesis:

$\overline{u_{i}\theta}=-\alpha_{t}\frac{\partial T}{\partial x_{i}},$ (26)

where $\alpha_{t}$ is modelled by Manservisi and Meneghini [6] as a function
of $k_{\theta}$ and $\epsilon_{\theta}$, which are computed by solving two
additional transport equations. This model was developed to be applied in
combination with a $k$-$\epsilon$ type model that is detailed in reference
[6]. The Kays correlation [33] calculates $\alpha_{t}$ based on a variable
turbulent Prandtl number, which is a function of $Re_{t}$ and $Pr$:

$\displaystyle\alpha_{t}$ $\displaystyle=\frac{\nu_{t}}{Pr_{t}},$
$\displaystyle Pr_{t}$ $\displaystyle=0.85+\frac{0.7}{Re_{t}Pr}.$ (27)

  * 2.

Shams et al. [7] developed an algebraic model for $\overline{u_{i}\theta}$:

$\overline{u_{i}\theta}=-C_{t0}\frac{k}{\epsilon}\left(C_{t1}\overline{u_{i}u_{j}}\frac{\partial
T}{\partial x_{j}}+C_{t2}\overline{u_{j}\theta}\frac{\partial U_{i}}{\partial
x_{j}}+C_{t3}\beta\overline{\theta^{2}}g_{i}\right)+C_{t4}b_{ij}\overline{u_{j}\theta},$
(28)

where the values of the closure constants $C_{t0}$-$C_{t4}$ as well ad the
transport equation for the thermal variance $\overline{\theta^{2}}$ can be
found in [7]. Among the different AHFM formulations, this work considers the
AHFM-EB that was calibrated to be applied in combination with the Elliptic
Blending model (EBM) implemented in the code Saturne 777Seehttps://www.code-
saturne.org/cms/..

##  4 Description of the database

The database used to train the artificial neural network consists of DNS
datasets related to different forced convection flows at various $Re$ and $Pr$
numbers. The momentum and thermal statistics were provided by the authors
after time averaging the simulation results in a time interval related to the
characteristic time of the flows [34, 35, 36]. The details of the datasets are
summarized in Table 2. Most of these databases concern standard, two-
dimensional flow configurations simulated at various Reynolds and Prandtl
numbers and with different thermal boundary conditions. In particular, the
databases of Kawamura et al. [36], Tiselji et al. [34] and Bricteux et al.
[37, 38] consist of DNS simulations of non-isothermal turbulent channel flow
at Reynolds number ranging from $Re_{\tau}=180$ to 2000 and Prandtl number
from 0.01 to 0.71. Li et al. [39] provided the DNS data of a non-isothermal
turbulent boundary layer flow simulation at $Re_{\theta}$=800, and Prandtl
numbers of 0.71 and 0.2.
Beside these fundamental flows, two additional three-dimensional flow
configurations were used for training and test to promote the model
adaptability to real world applications. Specifically, the database of Oder et
al. [35] refers to a three-dimensional confined backward facing step flow at
$Re_{b}=3200$ and with an expansion ratio of 2.25. The forced convective heat
transfer for this flow configuration was simulated at Prandtl numbers of 0.1
and 0.005, which is typical of sodium. Moreover, we include the dataset by
Angeli et al. [40] which considers the convective heat transfer around
cylindrical rods arranged in bundles. This kind of non-isothermal flow is of
relevant interest in the nuclear reactor community, as it reproduces the
cooling of reactor core fuel assemblies. This flow reaches $Re_{\tau}=550$ in
the lattice subchannels and the $Pr=0.031$ was considered as representative of
liquid Lead-Bismuth Eutectic (LBE).
Table 2 also specifies the databases used for the training and the validation
of the data-driven model. The allocation of training and test data is based on
the completeness of the turbulent statistics provided (which must allow to
calculate all the invariants $\pi_{i}$ and basis tensors $\mathbf{T}_{i}$) and
the spatial resolution of the different datasets (which should allow the
accurate computation of the first order derivatives of the temperature and
velocity). The overall number of training points taken from the different
databases amounts to 425400 out of 500000.

Table 2: Available DNS databases for forced convection at different $Re$ and
$Pr$ numbers

Author | Flow | Reynolds | Prandtl | Usage
---|---|---|---|---
and Ref. | Configuration |  number* | number |
Kawamura et al.[36] | Channel flow |  $Re_{\tau}$=180-640 | 0.025-0.71 | training/test
Li et al. [39] | Turbulent BL |  $Re_{\theta}$=830 | 0.2, 0.71 | training/test
Tiselj et al. [34] | Channel flow |  $Re_{\tau}$=180-590 | 0.01 | training/test
Oder et al. [35] | BFS |  $Re_{b}$=3200 | 0.1,0.005 | training/test
Angeli et al. [40] | Rode Bundle Flow |  $Re_{\tau}=$550 | 0.031 | test
Bricteux et al. [37, 38] | Channel flow |  $Re_{\tau}=$180-2000 | 0.01-0.025 | test

  * *

The complete definition of the Reynolds number for each flow configuration can
be found in the related references.

##  5 Auxiliary turbulence models employed for the validation

The full RANS simulation of a forced convection flow requires combining a
thermal turbulence model with a momentum turbulence model computing the
Reynolds stresses. Due to the use of high-fidelity data for the training and
the intrinsic dependency of the thermal statistics on the momentum ones, a
significant sensitivity of the thermal model to the combined momentum
turbulence model is expected. Hence, the choice of the momentum turbulence
model is critical for the validation of the thermal closure and it could
largely affect the accuracy of the results achieved.

Two different wall-resolved momentum treatments were combined to the thermal
model for its validation in OpenFoam888See https://www.openfoam.com/:

  * 1.

The Elliptic Blending Reynolds Stress Model (EBRSM) developed by Manceau [41].
This is a second order model solving an additional elliptic equation for a
blending parameter accounting for blockage effects caused by the walls. This
model ensures an advanced representation of the Reynolds stress tensor and its
anisotropic part $\mathbf{b}$, which is involved in both the invariant basis
and the tensor basis. In particular, Figure 3 shows the agreement of the in-
house implementation of the EBRSM with the DNS data for a turbulent channel
flow at $Re_{\tau}=395$.

  * 2.

The Launder-Sharma $k-\epsilon$ model, which applies the eddy viscosity
concept and the Boussinesq approximation to compute the Reynolds stress
tensor. As a consequence, the anisotropic part $\mathbf{b}$ is roughly
approximated.

The isotropic thermal quantities $k_{\theta}$ and $\epsilon_{\theta}$
indicated in Table 1 are computed by solving the differential transport
equations proposed by Manservisi et al. [6] for liquid metal flow simulations.

Figure 3: Distribution of the Reynolds stresses obtained by applying the
Elliptic Blending model [41] in case of a turbulent channel flow at
$Re_{\tau}=395$.

##  6 Results

This section presents the results of the training and the validation of the
data-driven model. The a priori validation consists in calculating the heat
flux predictions over high-fidelity (DNS) input data excluded from the
training dataset. The a posteriori validation involves the integration of the
model into the CFD code (OpenFoam) and the subsequent resolution of the
transport equations using the data-driven thermal turbulence closure.

###  6.1 Dispersion tensor

Figures 4 and 5 illustrate the components of the dispersion tensor
$\mathbf{D}$ obtained after the training. The tensor is computed over the
corresponding DNS input data for turbulent channel flow and rode bundle flow.
Specifically, Figure 4 shows the behaviour of the components obtained with and
without noise addition. The noise addition leads to more reasonable
distributions of the tensor components and mitigates the steep increase of
$\mathbf{D}$ in the vanishing gradient regions. The skew-symmetric part of
$\mathbf{D}$ is preserved by the data-driven formulation, as it can be seen by
comparing the components $D_{xy}$ and $D_{yx}$ in Figure 4. The magnitude of
the tensor decreases with the decrease of the Prandtl number, following the
damping of thermal turbulence by conduction shown in the high fidelity
simulations [42]. Figure 5 highlights that the predicted components of
$\mathbf{D}$ naturally satisfy all the symmetries characterising the rode
bundle flow configuration, thanks to the invariance properties enforced into
the model structure. The anisotropy of $\mathbf{D}$ can be visualised by
representing its eigenvalues on a barycentric map [43], to show the distance
of the given anisotropic state to the limiting states of turbulence. Given
$\psi_{i}$ the eigenvalues of the anisotropic part of $\mathbf{D}$, these
states are:

  * 1.

one-component ($x_{1c}$) for which
$\psi_{i}=\frac{2}{3},-\frac{1}{3},-\frac{1}{3}$ and the turbulent
fluctuations only exist along one direction;

  * 2.

two-components ($x_{2c}$), for which
$\psi_{i}=\frac{1}{6},\frac{1}{6},-\frac{1}{3}$ and the turbulent fluctuations
exist along two directions with equal magnitude;

  * 3.

isotropic ($x_{3c}$) for which $\psi_{i}$ are all zero.

Figure 6 shows the anisotropic part of $\mathbf{D}$ mapped in barycentric
coordinates in case of turbulent channel flow at Prandtl numbers of 0.71,
0.025 and 0.01. The behaviour of the anisotropic part along the channel
resembles the one of the Reynolds stresses for the same flow configuration. It
can be seen that passing from $Pr=0.71$ to $Pr=0.025$ and 0.01, the
anisotropic part tends to depart from the one-component state. This is in
agreement with the expected behaviour of thermal turbulence at low Prandtl
numbers, which shows a lower degree of anisotropy and redistribution due to
the significant molecular diffusion effects.

(a) $Pr=0.71$

(b) $Pr=0.025$

Figure 4: Components of the dispersion tensor predicted by the ANN after the
training for a non-isothermal channel flow at $Re_{\tau}=395$. The solid lines
indicate the results obtained by perturbing the parameter $\pi_{7}$. Dashed
lines represent the results obtained without noise addition. Figure 5:
Contours of the components of the dispersion tensor $\mathbf{D}$ by the ANN
using the DNS input data provided by Angeli et al. [40]

(a) $Re_{\tau}=640$

(b) $Re_{\tau}=395$

(c) $Re_{\tau}=180$

Figure 6: Visualization of the anisotropy of $\mathbf{D}$ on the barycentric
map as predicted by the ANN using the available DNS data of Kawamura et al.
[36] (non-isothermal turbulent channel flow) as input data.

###  6.2 A priori validation

Figure 7 shows the results of the neural network training for cases of non-
isothermal turbulent channel flow at different Prandtl numbers
($Pr=0.01-0.71$). It is clear that the ANN is able to fit the DNS data in this
wide range of Prandtl numbers, also matching the correct near wall behaviour
for both the components of the turbulent heat flux. The accuracy of the
predictions and the uncertainty ranges computed with the SWAG algorithm are
quite similar for both training and test data, showing that the network does
not overfit the training data.

The results of the training in the mid-plane for the non-isothermal backward
facing step flow are depicted in Figure 8. For this database, the portion of
data employed for the validation was randomly selected among the rectangular
partitions in which the entire flow field is divided to construct batches of
65 x 95 points. Also for this flow configuration, the network satisfactorily
represents both the components of the turbulent heat flux along the step. The
predicted field of the streamwise heat flux shows small inaccuracies close to
the lower wall, where the data-driven model slightly overestimates the
extension of the thermal separation region characterised by a change of sign
in the streamwise turbulent heat flux. The lower accuracy detected in this
region is probably due to the limitations of the algebraic mathematical
structure, which is inherently local and makes it hard to follow the abrupt
evolution of the thermal field after the step.

(a) $Pr=0.71$

(b) $Pr=0.025$

(c) $Pr=0.01$

Figure 7: Streamwise and wall-normal components of the turbulent heat flux
predicted by the ANN after the training. Predictions in the training range and
validation range are depicted with dashed and dotted lines, respectively.
Shaded areas represent the uncertainty ranges of the heat flux predictions
estimated with dropout. DNS data [36] are indicated with solid lines.

The data-driven model was then used to predict the turbulent heat flux over
the DNS data of the rode bundle flow [40], which was not used for the model
training. The turbulent heat flux predictions in the $x$, $y$ and $z$
directions are depicted in Figures 10, 10 and 11, where they are compared with
the corresponding DNS data and the values obtained by applying the Kays
correlation introduced in section 3. The data driven model slightly
underestimates the peak values achieved by the cross-flow heat flux in the
gaps, though its distribution agrees very well with the DNS counterpart and it
is much better than the distribution given by the Kays correlation, for which
the highest peaks are located close to the bundle walls. The current model
also well predicts the distribution of the streamwise turbulent heat flux
(Figure 11), which is identically zero for the Kays correlation due to the
standard gradient assumption.

(a) Stream-wise heat flux

(b) Wall-normal heat flux

Figure 8: Stream-wise and wall-normal components of the turbulent heat flux
predicted by the ANN after the training for a non-isothermal, backward-facing
step flow at $Pr=0.1$. Predictions of training and validation data are
depicted with blue and red dashed lines, respectively. DNS data [35] are
represented with solid lines.

Figure 9: Contours of the x-component of the turbulent heat flux for the rode
bundle flow. Left to right: DNS data [40], predictions of the data-driven
model and predictions of the Kays correlation.

Figure 10: Contours of the y-component of the turbulent heat flux for the rode
bundle flow. Left to right: DNS data [40], predictions of the data-driven
model and predictions of the Kays correlation. Figure 11: Contours of the
z-component of the turbulent heat flux for the rode bundle flow. Left to
right: DNS data [40], predictions of the data-driven model and predictions of
the Kays correlation.

###  6.3 A posteriori validation

The a posteriori validation of the data driven model was carried out in
OpenFoam 6.0. The structure of the artificial neural network and its trained
parameters were implemented into an existing library for thermal turbulence
models. A similar implementation procedure is described in detail in the
recent work of Maulik et al. [44]. The test cases selected for the a
posteriori validation consist of non-isothermal channel flows at $Re_{\tau}$
up to 2000 and $Pr$ numbers of 0.025 and 0.01, and the backward facing step
flow indicated in Table 2 at Prandtl number equal to 0.1.

(a) $Pr=0.01$

(b) $Pr=0.025$

Figure 12: Comparison of the wall-normal heat flux computed with the data-
driven model with the corresponding DNS [36] and LES data [38] for simulations
of non-isothermal turbulent channel flows at $Re_{\tau}$ up to 2000.

(a) $T$

(b) $k_{\theta}$

(c) $\overline{u\theta}$

(d) $\overline{v\theta}$

Figure 13: Comparison of the results achieved with the data-driven model (ML),
the Manservisi model (MM), the Kays correlation and the AHFM model at
$Re_{\tau}$=395 and $Pr$=0.01.

(a) $T$

(b) $k_{\theta}$

(c) $\overline{u\theta}$

(d) $\overline{v\theta}$

Figure 14: Comparison of the results achieved with the data-driven model (ML),
the Manservisi model (MM), the Kays correlation and the AHFM model at
$Re_{\tau}$=395 and $Pr$=0.025 as thermal turbulence models, and the EBM
(left) and the Launder-Sharma $k$-$\epsilon$ (right) as momentum turbulence
models.

Figure 12 shows the wall-normal turbulent heat flux predicted by the channel
flow simulations carried out at different Reynolds and Prandtl numbers with
the EBRSM as momentum turbulence model. Even if the ANN was trained up to
$Re_{\tau}=640$ at $Pr=0.025$ and up to $Re_{\tau}=590$ at $Pr=0.01$ (see
Table 2), the model is able to extrapolate well at higher $Re_{\tau}$ and
reproduces the correct heat fluxes at $Re_{\tau}=2000$. These generalisation
properties are very important for the practical applicability of the
formulation, as industrial problems usually involve much higher Reynolds
numbers than those characterizing the current training database. It can be
noted that the prediction is slightly less accurate for $Pr=0.01$ than for
$Pr=0.025$ at $Re_{\tau}=2000$. This shows that the conditions of high
Reynolds and very low Prandtl numbers are still the most critical for the
data-driven model and that the current low-Reynolds database should be
integrated with higher Reynolds number flows to further improve the agreement
of the network predictions in these regimes.

Figures 13 and 14 (left) depicts the results of the channel flow simulations
for $Re_{\tau}=395$ and $Pr=0.025$ and999Note that these conditions were
excluded from the training dataset to ensure a meaningful validation. 0.01 and
compare them with those obtained with the models presented in section 3. Note
that these conditions were excluded from the training dataset to ensure a
meaningful validation.

The data-driven model shows to be very accurate for this kind of flow in terms
of all the thermal statistics, except for the stream-wise component of the
heat flux, which is more accurately predicted by the AHFM at $Pr=0.01$. The
predictions are more accurate than the models based on standard gradient
assumptions [6, 33], as well as the anisotropic formulation developed by Shams
et al. [7]. This confirms the potential of ANNs in the field of turbulence
modelling, and shows that a wide exploration of the parameter space is useful
to develop accurate and advanced thermal turbulence closures.

Nevertheless, as mentioned in section 3, each of these thermal models was
developed to be applied in combination with a certain momentum turbulence
model, which differs from the current EBRSM. Hence, another comparison of the
models was carried out over a less accurate momentum field computed with the
Launder-Sharma $k-\epsilon$ model. The results at $Pr=0.025$ and
$Re_{\tau}=395$ are shown in Figure 14 (right): the data driven model
overestimates the wall-normal heat flux by almost 50%, being less robust than
the other thermal models which show the ability to better adapt to different
momentum treatments. This lack of robustness is due to the anisotropic
character of the data-driven formulation and its complexity, i.e. the wide
range of parameters accounted (Table 1) which are not all accurately estimated
by the auxiliary turbulence models and directly propagate their inaccuracies
to the thermal turbulence predictions.

The second test case was aimed at reproducing the backward facing step flow
simulated by Oder et al. [35] at $Pr=0.1$. The original DNS simulation showed
a secondary flow recirculation in the spanwise direction due to the presence
of the lateral walls. However, the spanwise velocity components in the
midplane are practically zero and the flow can be considered locally two-
dimensional. Based on this consideration, the flow on the midplane was
simulated in OpenFoam with a two-dimensional geometry. The results of the
simulations are depicted in Figure 15. As shown by figures 15a, the secondary
circulation is clearly not captured by the two-dimensional setup, leading to
significantly different vertical velocity distributions, especially close to
the step. Although this velocity component is small compared to the horizontal
one, its distribution affects the strain and rotational tensors involved in
the parametrization (eq. (8)). As a consequence, the turbulent heat flux
obtained at convergence (figures 15c and 15d) differs from its DNS
counterpart, especially near the lower wall. In particular, the region of
thermal separation, as well as the components of the heat flux in the
reattachment zone are overestimated. Given $T$ the DNS temperature field and
$\hat{T}$ the one obtained with RANS simulation, the deviation from the
reference DNS data (Figure 15b) is computed by integrating the pointwise
relative error over the vertical axis:

$e=\frac{1}{L_{y}}\int_{y_{min}}^{y_{min}+L_{y}}\frac{\hat{T}-T}{T}dy,$ (29)

where $L_{y}$ indicates the height of the channel after the step. The
behaviour of the relative error along the streamwise direction is shown in
Figure 16, which compares the level of accuracy of the data-driven model and
the model of Manservisi. Both models lead to deviations of around 20% very
close to the step. However, the Manservisi model leads to more accurate
temperature predictions in the separation zone ($0<x<12$) than the data-driven
model.

Such inaccuracies again underline the sensitivity of the model to the quality
of the momentum field, which is a result of the training performed using high-
fidelity data. This sensitivity places limits to the applicability of the
model and its reliability in cases where the computed momentum field is
affected by high inaccuracies in the second order statistics.

(a) Vertical velocity

(b) Temperature

(c) Streamwise turbulent heat flux

(d) Wall normal heat flux

Figure 15: Results of the simulation of a backward-facing step flow at
$Pr$=0.1 with the elliptic blending model of Manceau et al. [41] and the data-
driven thermal turbulence model (ML). Comparison with the corresponding DNS
data [35]. Figure 16: Top: comparison of the temperature profiles obtained
with the data-driven model (blue) and the Manservisi model (green) with the
reference DNS data (red) [35]. Bottom: Integral deviation of the temperature
profiles obtained with the models from the reference DNS data.

##  7 Model interpretation

The data-driven model was further analysed to have a clearer view of its
performance and limitations. The model was interpreted by applying the
Integrated Gradient method [45] which is able to attribute the variation from
a baseline to a target solution to all the parameters involved in the data-
driven model. The method is based on the integration of the gradients of the
outputs $\mathbf{Y}$ with respect to the inputs $\mathbf{X}$ over a linear
path from the baseline to the target solutions parametrized by the coordinate
$\gamma$:

$\centering\underbrace{\mathbf{I_{G}}}_{\text{Attribution value
matrix}}=(\underbrace{\mathbf{X^{\prime}}}_{\begin{subarray}{c}\text{Target}\\\
\text{input}\\\
\text{vector}\end{subarray}}-\underbrace{\mathbf{X}}_{\begin{subarray}{c}\text{Baseline}\\\
\text{input}\\\
\text{vector}\end{subarray}})\int_{\gamma=0}^{1}\frac{\partial\overbrace{\mathbf{Y}}^{\text{Output
vector}}(\mathbf{X^{\prime}}+\gamma(\mathbf{X}-\mathbf{X^{\prime}}))}{\partial\mathbf{X}}d\gamma\@add@centering$
(30)

where $\mathbf{I_{G}}$ is the attribution value matrix, giving an estimate of
the sensitivity of the data-driven function on the input parameters, while
$\mathbf{X}$ and $\mathbf{X^{\prime}}$ are the baseline and target vectors of
the inputs. In the present case, the baseline and the target solutions are
represented by the channel flow predictions of $\mathbf{D}$ ($Pr=0.025$)
obtained with the EBRSM and the $k-\epsilon$ model as momentum turbulence
models (Figure 14).

For this particular flow configuration, the component $D_{yy}$ is the only
affecting the temperature field, as the flow is fully developed. For this
component, the highest attribution values assigned to the parameters to
explain the variation from the baseline to the target solutions are indicated
in Figure 17, for different wall distances. It is clear that the most
important tensors and coefficients are $\mathbf{T}_{1}$, $\mathbf{T}_{2}$
and101010This evidence indicates that the data-driven formulation is
effectively an extension of the Daly and Harlow [46] and Higher Order
Generalised Gradient Diffusion [42] approaches. $a_{1}$, $a_{2}$. The same
analysis is then repeated for the coefficients $a_{1}$-$a_{8}$ to attribute
their variations to the parameters $\pi_{i}$. The results for $a_{1}$ and
$a_{2}$ are shown in Figure 18. For both the coefficients, the variation is
mostly explained by the parameters $\pi_{3}$ and $\pi_{6}$, whose definition
is reported in Table 1. Indeed, the parameter $\pi_{3}$ changes significantly
with the change of the momentum turbulence model, since the $k-\epsilon$ model
provides a rough approximation of the anisotropic part of the Reynolds
stresses. $\pi_{6}$ is identically zero for the target solution as a direct
consequence of the Boussinesq approximation. Hence, this analysis of the data-
driven model proved its sensitivity to the Reynolds stress anisotropy, thus
explaining why the heat flux predictions significantly deteriorates when the
thermal turbulence model is coupled with a momentum turbulence model relying
on the Boussinesq approximation.

Figure 17: Maximum 10 attribution values assigned with the integrated gradient
method to the coefficients $a_{i}$ and the tensors $\mathbf{T}^{i}$ involved
in the parametrization, at different distances from the wall.

Figure 18: Attributions assigned with the integrated gradient method to the
input parameters of the ANN for the coefficients $a_{1}$ and $a_{2}$.

##  8 Conclusions and Perspectives

This work presents a data-driven approach for the modelling of the turbulent
heat flux with special focus to low Prandtl number flows. The algebraic
structure of the model embeds invariance and realizability properties, and the
training methodology ensured smooth predictions and good generalisation
properties. The training database covered a wide range of Prandtl numbers
($Pr=0.01-0.71$) to achieve a general and unified formulation.

The a priori and a posteriori validation concerned different two-dimensional
and three-dimensional flow configurations in the low Prandtl number regime.
The results showed that the data-driven model performs well both within and
outside the range of conditions explored in the training data. When compared
to other thermal turbulence models for low Prandtl numbers, the model proved
to be very accurate and satisfactorily reproduced the different components of
the turbulent heat flux.

Due to the use of DNS data during training, the data-driven formulation is
quite sensitive to the accuracy of the momentum field and the momentum
turbulence model used in combination with the thermal one. The data data-
driven formulation seems not robust enough against the defects of the momentum
modelling, especially in case of complex flow conditions involving separation
and reattachment, where even the most sophisticated models lack of accuracy.
Hence, further efforts will be moved towards the development of a more robust
formulation for its employment in large scale problems. The robustness could
be practically enhanced by acting on the training methodology, e.g. by
perturbing the quality of DNS input data applying popular RANS approximations
(e.g. Boussinesq approximation) to reconstruct the Reynolds stress tensor.
This training approach would make the model less sensitive to the accuracy of
the momentum turbulence statistics.

This higher robustness would be at the expense of introducing compensation
errors in the model structure, which are counterproductive to gain a better
understanding of thermal turbulence, especially at low Prandtl numbers. Hence,
the future of the research in this field will move forward in two distinct
directions. On one side, machine learning and data-driven techniques will be
used to develop new closures, explore new solutions and gain knowledge from
the formulations retrieved. A promising field of research would be extending
the use of data-driven methods to develop and calibrate differential models,
which are able to introduce non-local and historical effects and do not suffer
for ill-posedness in the vanishing gradient regions. A differential approach
will thus contribute to improve the behaviour of the model in presence of
wall-detached and separated flows, where the turbulence production is
typically negligible and transport and redistribution effects govern the
dynamics of turbulence. On the other hand, we will seek to find a good
compromise between the accuracy and the practical applicability of the thermal
model and ensure its compatibility with the other sub-models involved in
complex heat transfer simulations.

## Appendix

The derivation of the tensorial set given by equation (11) comes from the
repeated use of the Cayley-Hamilton theorem, which for a two dimensional
tensors can be written as:

$\mathbf{M}^{2}=\\{\mathbf{M}\\}\mathbf{M}-\frac{1}{2}(\\{\mathbf{M}\\}^{2}-\\{\mathbf{M}^{2}\\})\mathbf{I}_{2}$
(31)

and in case of three dimensional tensors as :

$\mathbf{M}^{3}=\\{\mathbf{M}\\}\mathbf{M}^{2}-\frac{1}{2}(\\{\mathbf{M}\\}^{2}-\\{\mathbf{M}^{2}\\})\mathbf{M}+\det(\mathbf{M})\mathbf{I}$
(32)

Given the tensor $M_{ij}^{r}$ with $r=1,...,n$, an invariant can be expressed,
without loose of generality, in the form:

$\pi=\beta_{j_{1},k_{1},j_{2},k_{2},...,j_{l},k_{l}}M_{j_{1},k_{1}}^{s_{1}}M_{j_{2},k_{2}}^{s_{2}}...M_{j_{l},k_{l}}^{s_{l}}$
(33)

where $s_{1},s_{2},...,s_{l}$ are integers (not necessarly all different)
chosen from $1,2,...,\mu$ and
$\beta_{j_{1},k_{1},j_{2},k_{2},...,j_{l},k_{l}}$ are a set of numerical
coefficients. Due to the required invariance of $\pi$ under any orthogonal
transformation, $\beta_{j_{1},k_{1},j_{2},k_{2},...,j_{l},k_{l}}$ would be the
components of an isotropic tensor. Hence,
$\beta_{j_{1},k_{1},j_{2},k_{2},...,j_{l},k_{l}}$ can be expressed as a
combination of the type:

$\delta_{i_{\alpha},i_{\beta}}\delta_{i_{\gamma},i_{\delta}}...\delta_{i_{\sigma},i_{\tau}}$
(34)

As a result, $\pi$ would be expressed as polynomial of terms of type:

$\Pi_{ii}=M_{i,k}^{s_{1}}M_{k,l}^{s_{2}}...M_{n,i}^{s_{m}}$ (35)

which is the trace of an arbitrary product of tensors within the defined
tensorial set. However, these terms constitute a not finite set. The minimal
integrity basis for $\Pi_{ii}$ will be derived with the aid of eq. (31). In
particular, it will be shown that the trace of any tensor product can be
expressed as a polynomial in traces of matrix products, each of them
satisfying the following properties:

  1. (i)

It is a product of factors of the form $\mathbf{M}$ and $\mathbf{M}^{2}$.
Indeed, factors higher than $\mathbf{M}^{2}$ are reducible, as it can be shown
by multiplying eq.(31) by $\mathbf{M}$ and taking the trace:

$\\{\mathbf{M}^{3}\\}=\\{\mathbf{M}\\}\\{\mathbf{M}^{2}\\}-\frac{1}{2}\left(\\{\mathbf{M}\\}^{2}-\\{\mathbf{M}^{2}\\}\right)\\{\mathbf{M}\\}$
(36)

  2. (ii)

If it contains the factor $\mathbf{M}^{2}$ it has no other factor. This can be
demonstrated by multiplying eq. (31) for any tensor $\mathbf{A}$ and taking
the trace:

$\\{\mathbf{M}^{2}\mathbf{A}\\}=\\{\mathbf{M}\\}\\{\mathbf{A}\\}-\frac{1}{2}\left(\\{\mathbf{M}\\}^{2}-\\{\mathbf{M}^{2}\\}\right)\\{\mathbf{A}\\}$
(37)

Hence, the trace of each product of the form $\mathbf{M}^{2}\mathbf{A}$ is
reducible.

  3. (iii)

The first and the last factors of the product are not the same tensors. It
follows by substituting $\mathbf{M}=\mathbf{A}+\mathbf{B}$ in eq. (31):

$\mathbf{AB}+\mathbf{BA}=\\{\mathbf{B}\\}\mathbf{A}+\\{\mathbf{A}\\}\mathbf{B}-\left(\\{\mathbf{A}\\}\\{\mathbf{B}\\}-\\{\mathbf{AB}\\}\right)\mathbf{I}$
(38)

Then, if we multiply eq. (38) by $\mathbf{A}$ and take the trace, we get:

$\begin{split}\\{\mathbf{ABA}\\}+\\{\mathbf{BA}^{2}\\}=\\{\mathbf{B}\\}\\{\mathbf{A}\\}+\\{\mathbf{A}\\}\\{\mathbf{BA}\\}-\bigl{(}\\{\mathbf{B}\\}\\{\mathbf{A}\\}\\\
-\\{\mathbf{AB}\\}\bigr{)}\\{\mathbf{A}\\}\end{split}$ (39)

Since the trace of a tensor product is invariant under the cyclic permutation
of its factors, then $\\{\mathbf{BA}^{2}\\}=\\{\mathbf{A}^{2}\mathbf{B}\\}$.
Due to the property 2, the product $=\\{\mathbf{A}^{2}\mathbf{B}\\}$ is
reducible. Hence, from eq. (39) $\\{\mathbf{ABA}\\}$ is also reducible.

  4. (iv)

Not two factors of the product are the same. This condition directly arises
from conditions 2 and 3. Indeed, by multiplying eq. (38) by $\mathbf{AC}$ and
taking its trace, we get:

$\begin{split}\\{\mathbf{ABAC}\\}+\\{\mathbf{B}\mathbf{A}^{2}\mathbf{C}\\}=\\{\mathbf{B}\\}\\{\mathbf{A}^{2}\mathbf{C}\\}+\\{\mathbf{A}\\}\\{\mathbf{BAC}\\}-\\\
\left(\\{\mathbf{A}\\}\\{\mathbf{B}\\}-\\{\mathbf{AB}\\}\right)\\{\mathbf{AC}\\}\end{split}$
(40)

Since $\\{\mathbf{BA}^{2}\mathbf{C}\\}=\\{\mathbf{A}^{2}\mathbf{CB}\\}$, which
is reducible due to property 2, then $\\{\mathbf{BA}^{2}\mathbf{C}\\}$ is also
reducible.

  5. (v)

The maximum total degree of the product is four. Indeed, it can be proved that
traces of tensor products of degree higher than four are reducible. Let’s
consider eq. (38), multiplied by $\mathbf{C}$:

$\mathbf{ABC}+\mathbf{BAC}=\\{\mathbf{B}\\}\mathbf{AC}+\\{\mathbf{A}\\}\mathbf{BC}-\left(\\{\mathbf{A}\\}\\{\mathbf{B}\\}-\\{\mathbf{AB}\\}\right)\mathbf{C}$
(41)

Taking the trace of eq. (41), we get:

$\\{\left(\mathbf{AB}+\mathbf{BA}\right)\mathbf{C}\\}\doteq 0$ (42)

where the notation $\doteq 0$ means that the trace of the product on the left
hand side is reducible. Let’s substitute $\mathbf{A}$ with $\mathbf{AM}$ and
$\mathbf{B}$ with $\mathbf{BM}$ in eq. (42):

$\\{\left(\mathbf{AMBN}+\mathbf{BNAM}\right)\mathbf{C}\\}\doteq 0$ (43)

This means that the trace of a tensor product is equivalent to the negative of
the trace of the same tensor product having two adjacent factors interchanged,
i.e. the two traces differ for polynomials of traces tensors of lower degree.
$\\{\mathbf{BNAMC}\\}$ can be obtained from $\mathbf{AMBNC}$ by an even number
of interchanges (4) of adjacent factors. Hence, $\\{\mathbf{BNAMC}\\}$ is
equivalent to $\\{\mathbf{AMBNC}\\}$, which is then reducible:

$\\{\mathbf{AMBNC}\\}\doteq 0$ (44)

In particular, we are dealing with tensors $\mathbf{b}$, $\mathbf{S}$ and
$\mathbf{\Omega}$ defined by eq. (3), (4) and (5). Under the assumption od 2-D
flow, these tensors have the following form:

$\mathbf{b}=\begin{pmatrix}b_{11}&b_{12}&0\\\ b_{12}&b_{22}&0\\\
0&0&b_{33}\end{pmatrix}$ $\mathbf{S}=\begin{pmatrix}S_{11}&S_{12}&0\\\
S_{12}&S_{22}&0\\\ 0&0&0\end{pmatrix}$
$\mathbf{\Omega}=\begin{pmatrix}0&\Omega_{12}&0\\\ -\Omega_{12}&0&0\\\
0&0&0\end{pmatrix}$

where $\\{\mathbf{b}\\}$, $\\{\mathbf{S}\\}$ and $\\{\mathbf{\Omega}\\}$ are
nil. It can be noted that $\mathbf{b}$ is not 2-D, hence it satisfies the
Cayley-Hamilton theorem written for 3-D tensors (eq.(32)). On the other hand,
we can easily verify that products of $\mathbf{b}$ with other 2-D tensors are
equal to the products of these latter with its corresponding 2-D projection,
$\mathbf{b}_{2}$:

$\mathbf{b}_{2}=\begin{pmatrix}b_{11}&b_{12}&0\\\ b_{12}&b_{22}&0\\\
0&0&0\end{pmatrix}$ (45)

where

$\mathbf{b}_{2}=\mathbf{b}+\\{\mathbf{b}_{2}\\}\mathbf{I}_{2}$ (46)

Hence, without loss of generality, we can reformulate the problem by seeking
the minimal basis of invariants of the following 2 tensorial sets:

$\mathbf{b}\text{ (3D) }$ $\mathbf{b}_{2},\mathbf{S},\mathbf{\Omega}\text{
(2D) }$

From eq. (46) it follows that:

$\\{\mathbf{b}^{2}\\}=\\{\mathbf{b}_{2}^{2}\\}+\\{\mathbf{b}_{2}\\}^{2}$ (47)

Hence, $\mathbf{b}_{2}^{2}$ can be expressed as a polynomial in
$\\{\mathbf{b}^{2}\\}$, $\\{\mathbf{b_{2}}\\}$. By multiplying eq. (46) by
$\mathbf{b}^{2}$ and taking its trace, we also get:

$\\{\mathbf{b}^{3}\\}=\\{\mathbf{b}_{2}^{3}\\}-\\{\mathbf{b}_{2}\\}^{3}$ (48)

$\\{\mathbf{b}_{2}^{3}\\}$ is reducible due to eq. (48), hence
$\\{\mathbf{b}^{3}\\}$ can be also removed from the basis.
Based on the above results and the conditions i-v, the following basis of
invariants can be constructed:

$\\{\mathbf{b_{2}}\\},\\{\mathbf{b}^{2}\\},\\{\mathbf{S}^{2}\\},\\{\mathbf{\Omega}^{2}\\},\\{\mathbf{bS}\\},\\{\mathbf{bS\Omega}\\},\\{\mathbf{b\Omega
S}\\},\\{\mathbf{\Omega bS}\\},\\{\mathbf{\Omega
Sb}\\},\\{\mathbf{Sb\Omega}\\},\\{\mathbf{S\Omega b}\\}$

This basis can be further reduced. Indeed, eq. (41) can be used to derive six
similar relationships for the traces of the products given by all the possible
permutations of $\mathbf{b}$, $\mathbf{S}$ and $\mathbf{\Omega}$:

$\\{\mathbf{bS\Omega}\\}+\\{\mathbf{Sb\Omega}\\}\doteq 0$ $\\{\mathbf{\Omega
bS}\\}+\\{\mathbf{b\Omega S}\\}\doteq 0$ $\\{\mathbf{S\Omega
b}\\}+\\{\mathbf{\Omega Sb}\\}\doteq 0$ $\\{\mathbf{b\Omega
S}\\}+\\{\mathbf{bS\Omega}\\}\doteq 0$
$\\{\mathbf{Sb\Omega}\\}+\\{\mathbf{S\Omega b}\\}\doteq 0$ $\\{\mathbf{\Omega
Sb}\\}+\\{\mathbf{\Omega bS}\\}\doteq 0$

The resulting homogeneous system has rank 5, hence there is only one
independent trace among the traces of products of degree 3.
Based on such considerations, the invariant basis becomes:

$\\{\mathbf{b_{2}}\\},\\{\mathbf{b}^{2}\\},\\{\mathbf{S}^{2}\\},\\{\mathbf{\Omega}^{2}\\},\\{\mathbf{bS}\\},\\{\mathbf{bS\Omega}\\}$

Similar considerations can be made to derive the basis of independent tensor
products given the set $\mathbf{b}$, $\mathbf{S}$ and $\mathbf{\Omega}$. In
particular, the basis will be given by:

$\begin{split}\mathbf{I},\mathbf{b},\mathbf{S},\mathbf{\Omega},\mathbf{bS},\mathbf{S\Omega},\mathbf{b\Omega},\mathbf{\Omega
SW}\end{split}$

## Nomenclature

$\mathbf{U}$ | [m/s] | Mean velocity
---|---|---
$T$ | [K] | Mean Temperature
$k$ |  [m2/s2] | Turbulent kinetic energy
$\epsilon$ |  [m2/s3] | Turbulent dissipation rate
$k_{\theta}$ |  [K2] | Thermal variance
$\epsilon_{\theta}$ |  [K2/s] | Thermal dissipation rate
$\overline{\mathbf{u}\mathbf{u}}$ |  [m2/s2] | Reynolds stress
$\overline{\mathbf{u}\theta}$ | [mK/s] | Turbulent heat flux
|  |

$\mathbf{g}$ |  [m/s2] | Gravity vector
---|---|---
$\nu$ |  [m2/s] | Molecular viscosity
$\alpha$ |  [m2/s] | Molecular diffusivity
$\alpha_{T}$ |  [m2/s] | Turbulent diffusivity
$\delta$ | [m] | Half channel width
$\gamma$ | [-] | Blending parameter
$\mathbf{I}$ | [-] | Identity tensor
|  |

## Acknowledgments

M. Fiore is supported by a F.R.S.-FNRS FRIA grant and gratefully acknowledges
Prof. Diego Angeli, Prof. Iztok Tiselj and Dr. Jure Oder for providing their
datasets.

## References

  * Roelofs [2018] F. Roelofs, Thermal hydraulics aspects of liquid metal cooled nuclear reactors, Woodhead Publishing, 2018.
  * Eckert et al. [2007] S. Eckert, A. Cramer, G. Gerbeth,  Velocity measurement techniques for liquid metal flows,  in: Magnetohydrodynamics, Springer, 2007, pp. 275–294.
  * Schulenberg and Stieglitz [2010] T. Schulenberg, R. Stieglitz,  Flow measurement techniques in heavy liquid metals,  Nuclear Engineering and Design 240 (2010) 2077–2087.
  * Geankoplis [2003] C. J. Geankoplis, Transport processes and separation process principles:(includes unit operations), Prentice Hall Professional Technical Reference, 2003\.
  * Grötzbach [2013] G. Grötzbach,  Challenges in low-Prandtl number heat transfer simulation and modelling,  Nuclear engineering and design 264 (2013) 41–55.
  * Manservisi and Menghini [2014] S. Manservisi, F. Menghini,  A CFD four parameter heat transfer turbulence model for engineering applications in heavy liquid metals,  International Journal of Heat and Mass Transfer 69 (2014) 312–326.
  * Shams et al. [2019] A. Shams, A. De Santis, F. Roelofs,  An overview of the AHFM-NRG formulations for the accurate prediction of turbulent flow and heat transfer in low-Prandtl number flows,  Nuclear Engineering and Design 355 (2019) 110342.
  * Suga and Abe [2000] K. Suga, K. Abe,  Nonlinear eddy viscosity modelling for turbulence and heat transfer near wall and shear-free boundaries,  International journal of heat and fluid flow 21 (2000) 37–48.
  * Carteciano et al. [1997] L. Carteciano, D. Weinberg, U. Müller,  Development and analysis of a turbulence model for buoyant flows,  in: 4th World Conference of Exp. Heat Transfer, Fluid Mechanics and Thermodynamics, Bruxelles, volume 3, 1997, pp. 1339–1347.
  * Shams [2018] A. Shams,  The importance of turbulent heat transfer modelling in low-prandtl fluids,  in: Advances in Thermal Hydraulics (ATH 2018), ANS Winter Meeting Embedded Topical, November, Orlando, FL, 2018\.
  * Ling et al. [2016] J. Ling, A. Kurzawski, J. Templeton,  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance,  Journal of Fluid Mechanics 807 (2016) 155–166.
  * Schmelzer et al. [2019] M. Schmelzer, R. P. Dwight, P. Cinnella,  Machine learning of algebraic stress models using deterministic symbolic regression,  arXiv preprint arXiv:1905.07510 (2019).
  * Zhao et al. [2020] Y. Zhao, H. D. Akolekar, J. Weatheritt, V. Michelassi, R. D. Sandberg,  RANS turbulence model development using CFD-driven machine learning,  Journal of Computational Physics 411 (2020) 109413.
  * Weatheritt and Sandberg [2017] J. Weatheritt, R. Sandberg,  The development of algebraic stress models using a novel evolutionary algorithm,  International Journal of Heat and Fluid Flow 68 (2017) 298–318.
  * Jiang et al. [2020] C. Jiang, J. Mi, S. Laima, H. Li,  A novel algebraic stress model with machine-learning-assisted parameterization,  Energies 13 (2020) 258.
  * Wang et al. [2017] J.-X. Wang, J. Wu, J. Ling, G. Iaccarino, H. Xiao,  A comprehensive physics-informed machine learning framework for predictive turbulence modeling,  arXiv preprint arXiv:1701.07102 (2017).
  * Zhang and Duraisamy [2015] Z. J. Zhang, K. Duraisamy,  Machine learning methods for data-driven turbulence modeling,  in: 22nd AIAA Computational Fluid Dynamics Conference, 2015, p. 2460.
  * Duraisamy and Durbin [2014] K. Duraisamy, P. Durbin,  Transition modeling using data driven approaches,  in: Proceedings of the Summer Program, 2014, p. 427.
  * Diez Sanhueza [2018] R. Diez Sanhueza,  Machine Learning for RANS Turbulence Modelling of Variable Property Flows (2018).
  * Mangeon [2020] G. Mangeon, Modélisation au second ordre des transferts thermiques turbulents pour tous types de conditions aux limites thermiques à la paroi., Ph.D. thesis, Pau, 2020.
  * Pope [1983] S. Pope,  Consistent modeling of scalars in turbulent flows,  The Physics of Fluids 26 (1983) 404–408.
  * So et al. [2004] R. So, L. Jin, T. Gatski,  An explicit algebraic Reynolds stress and heat flux model for incompressible turbulence: Part I Non-isothermal flow,  Theoretical and Computational Fluid Dynamics 17 (2004) 351–376.
  * So and Sommer [1996] R. So, T. Sommer,  An explicit algebraic heat-flux model for the temperature field,  International journal of heat and mass transfer 39 (1996) 455–465.
  * Wikström et al. [2000] P. Wikström, S. Wallin, A. V. Johansson,  Derivation and investigation of a new explicit algebraic model for the passive scalar flux,  Physics of fluids 12 (2000) 688–702.
  * Lazeroms et al. [2013] W. Lazeroms, G. Brethouwer, S. Wallin, A. Johansson,  An explicit algebraic Reynolds-stress and scalar-flux model for stably stratified flows,  Journal of Fluid Mechanics 723 (2013) 91–125.
  * Spencer and Rivlin [1997] A. Spencer, R. Rivlin,  The Theory of Matrix Polynomials and its Application to the Mechanics of Isotropic Continua,  in: Collected Papers of RS Rivlin, Springer, 1997, pp. 1071–1098.
  * Pope [2001] S. B. Pope, Turbulent flows, 2001\.
  * Hadjesfandiari [2014] A. R. Hadjesfandiari,  On the symmetric character of the thermal conductivity tensor,  International Journal of Materials and Structural Integrity 8 (2014) 209–220.
  * Horn et al. [1994] R. A. Horn, R. A. Horn, C. R. Johnson, Topics in matrix analysis, Cambridge university press, 1994.
  * Sathyanarayana [2014] S. Sathyanarayana,  A gentle introduction to backpropagation,  Numeric Insight 7 (2014) 1–15.
  * Maddox et al. [2019] W. Maddox, T. Garipov, P. Izmailov, D. Vetrov, A. G. Wilson,  A simple baseline for bayesian uncertainty in deep learning,  arXiv preprint arXiv:1902.02476 (2019).
  * Sotgiu et al. [2019] C. Sotgiu, B. Weigand, K. Semmler, P. Wellinger,  Towards a general data-driven explicit algebraic Reynolds stress prediction framework,  International Journal of Heat and Fluid Flow 79 (2019) 108454.
  * Kays [1994] W. M. Kays,  Turbulent Prandtl number. Where are we?,  ASME Transactions Journal of Heat Transfer 116 (1994) 284–295.
  * Tiselj et al. [2001] I. Tiselj, R. Bergant, B. Mavko, I. Bajsic, G. Hetsroni,  DNS of turbulent heat transfer in channel flow with heat conduction in the solid wall,  J. Heat Transfer 123 (2001) 849–857.
  * Oder et al. [2019] J. Oder, A. Shams, L. Cizelj, I. Tiselj,  Direct numerical simulation of low-Prandtl fluid flow over a confined backward facing step,  International Journal of Heat and Mass Transfer 142 (2019) 118436.
  * Kawamura et al. [2000] H. Kawamura, H. Abe, K. Shingai,  DNS of turbulence and heat transport in a channel flow with different Reynolds and Prandtl numbers and boundary conditions,  Turbulence, Heat and Mass Transfer 3 (2000) 15–32.
  * Bricteux et al. [2012] L. Bricteux, M. Duponcheel, G. Winckelmans, I. Tiselj, Y. Bartosiewicz,  Direct and large eddy simulation of turbulent heat transfer at very low Prandtl number: Application to lead–bismuth flows,  Nuclear engineering and design 246 (2012) 91–97.
  * Duponcheel et al. [2014] M. Duponcheel, L. Bricteux, M. Manconi, G. Winckelmans, Y. Bartosiewicz,  Assessment of RANS and improved near-wall modeling for forced convection at low Prandtl numbers based on LES up to $Re_{\tau}$= 2000,  International Journal of Heat and Mass Transfer 75 (2014) 470–482.
  * Li et al. [2009] Q. Li, P. Schlatter, L. Brandt, D. S. Henningson,  DNS of a spatially developing turbulent boundary layer with passive scalar transport,  International Journal of Heat and Fluid Flow 30 (2009) 916–929.
  * Angeli et al. [2019] D. Angeli, A. Fregni, E. Stalio,  Direct numerical simulation of turbulent forced and mixed convection of LBE in a bundle of heated rods with P/D= 1.4,  Nuclear Engineering and Design 355 (2019) 110320.
  * Manceau [2005] R. Manceau,  An improved version of the elliptic blending model application to non-rotating and rotating channel flows,  in: Fourth International Symposium on Turbulence and Shear Flow Phenomena, Begel House Inc., 2005\.
  * Abe and Suga [2001] K. Abe, K. Suga,  Towards the development of a Reynolds-averaged algebraic turbulent scalar-flux model,  International Journal of Heat and Fluid Flow 22 (2001) 19–29.
  * Emory and Iaccarino [2014] M. Emory, G. Iaccarino,  Visualizing turbulence anisotropy in the spatial domain with componentality contours,  Center for Turbulence Research Annual Research Briefs (2014) 123–138.
  * Maulik et al. [2021] R. Maulik, H. Sharma, S. Patel, B. Lusch, E. Jennings,  Deploying deep learning in OpenFOAM with TensorFlow,  in: AIAA Scitech 2021 Forum, 2021, p. 1485.
  * Sundararajan et al. [2017] M. Sundararajan, A. Taly, Q. Yan,  Axiomatic attribution for deep networks,  arXiv preprint arXiv:1703.01365 (2017).
  * Daly and Harlow [1970] B. J. Daly, F. H. Harlow,  Transport equations in turbulence,  The Physics of Fluids 13 (1970) 2634–2649.