# Incorporating basic RANS calibrations in existing machine-learned turbulence
modeling

Jiaqi J.L. Li111PhD Candidate, Mechanical Engineering Pennsylvania State
University, State College, Pennsylvania, USA, 16802  Yuanwei Bin222PhD
Candidate, Mechanical Engineering Pennsylvania State University, State
College, Pennsylvania, USA, 16802  George P. Huang333Professor, Mechanical
Engineering Wright State University, Dayton, Ohio, USA, 45435  Xiang I.A.
Yang444Assistant Professor, Mechanical Engineering Pennsylvania State
University, State College, Pennsylvania, USA, 16802

###### Abstract

This work aims to incorporate basic calibrations of Reynolds-averaged Navier-
Stokes (RANS) models as part of machine learning (ML) frameworks. The ML
frameworks considered are tensor-basis neural network (TBNN), physics-informed
machine learning (PIML), and field inversion & machine learning (FIML) in J.
Fluid Mech., 2016, 807, 155-166, Phys. Rev. Fluids, 2017, 2(3), 034603 and J.
Comp. Phys., 2016, 305, 758-774, and the baseline RANS models are the one-
equation Spalart-Allmaras model, the two-equation $k$-$\omega$ model, and the
seven-equation Reynolds stress transport models. ML frameworks are trained
against plane channel flow and shear-layer flow data. We compare the ML
frameworks and study whether the machine-learned augmentations are detrimental
outside the training set. The findings are summarized as follows. The
augmentations due to TBNN are detrimental. PIML leads to augmentations that
are beneficial inside the training dataset but detrimental outside it. These
results are not affected by the baseline RANS model. FIML’s augmentations to
the two eddy viscosity models, where an inner-layer treatment already exists,
are largely neutral. Its augmentation to the seven-equation model, where an
inner-layer treatment does not exist, improves the mean flow prediction in a
channel. Furthermore, these FIML augmentations are mostly non-detrimental
outside the training dataset. In addition to reporting these results, the
paper offers physical explanations of the results. Last, we note that the
conclusions drawn here are confined to the ML frameworks and the flows
considered in this study. More detailed comparative studies and validation &
verification studies are needed to account for developments in recent years.

##  1 Nomenclature

$b_{ij}$ = | anisotropy tensor part of the normalized Reynolds stress tensor
---|---
$C_{b1}$, $C_{b2}$, $C_{w1}$, $C_{w2}$, $C_{w3}$ = | Spalart-Allmaras model constants
$d$ = | distance to the wall
$D_{ij}$ = | diffusion tensor in the Reynolds stress transport equations
${\rm D}/{\rm D}t$ = | material derivative
$f_{\nu 1}$, $f_{\nu 2}$, $f_{w}$, $g$ = | Spalart-Allmaras model functions
$I_{2}$, $I_{3}$ = | identity matrices
$k$ = | turbulent kinetic energy
$\mathcal{P}_{ij}$ = | production tensor in the Reynolds stress transport equations
$r$ = | Spalart-Allmaras model variable
$R_{ij}$ = | Reynolds stress tensor
$Re_{\tau}$ = | friction Reynolds number
$S_{ij}$, $\boldsymbol{S}$ = | strain rate tensor
$u$, $U$ = | velocity
$x_{i}$ = |  $i$-th Cartesian direction
$y$ = | vertical direction, wall-normal direction
$\beta$ = | model correction
$\theta$, $\theta^{*}$ = | Wilcox $k$-$\omega$ model constants
$\gamma$ = | model constant in the Wilcox $k$-$\omega$ model
$\delta$ = | channel half height
$\delta_{ij}$ = | Kronecker delta tensor
$\epsilon_{ij}$ = | dissipation tensor in the Reynolds stress transport equations
$\epsilon$ = | dissipation rate
$\kappa$ = | von Kármán constant
$\mu$ = | molecular dynamic viscosity
$\mu_{t}$ = | turbulent eddy viscosity
$\nu$ = | molecular kinematic viscosity
$\nu_{t}$, $\tilde{\nu}$ = | turbulent eddy viscosity
$\Phi_{ij}$ = | pressure-strain correlation tensor in the Reynolds stress transport equations
$\sigma$, $\sigma^{*}$ = | Spalart-Allmaras model constants
$\tau_{ij}$ = | Reynolds stress tensor
$\chi$ = | Spalart-Allmaraos model variable
$\omega$ = | specific dissipation rate
$\Omega_{ij}$, $\boldsymbol{R}$ = | rotation tensor

##  2 Introduction

Due to the high computational cost associated with scale-resolving tools [1,
2, 3], as well as the limited computing power, cost-effective, non-scale-
resolving tools like Reynolds averaged Navier Stokes (RANS) will remain the
workhorse for fluids engineering for the foreseeable future. RANS solves for
the mean flow while modeling the entirety of turbulence. Turbulence modeling
in the context of RANS is a topic with a long history [4, 5, 6]. Early models
are algebraic [7, 8] and have limited applicable ranges. The work in the 80s
and 90s led to a series of transport RANS models that are widely used today.
Notable examples include the one-equation Spalart-Allmaras model [9], the
various two-equation $k$-$\omega$ model [10, 11], the seven-equation full
Reynolds stress model [12, 13], among others [14, 15, 16]. These RANS models,
or variations thereof, are readily available in many commercial codes and are
the industrial standards. Although off-the-shelf RANS models are known to have
difficulty handling separation, boundary-layer three-dimensionality, and
unsteadiness, among others [17, 18, 19], they generally produce results that
are considered “reasonable” in the sense that they are usually not apparently
unphysical.

The story is different regarding machine-learning models. The use of machine-
learning tools for turbulence modeling dates back to as early as 2002 [20],
but their extensive use was not until the work in Refs. [21, 22] in 2015 and
2016. Since then, numerous researchers have delved into data-based approaches
for turbulence modeling, resulting in promising models that have demonstrated
exceptional accuracy for flows within their training sets [23, 24, 25].
However, akin to other machine-learning applications, these models face
limitations in their generalization capabilities. Worse still, when applied to
flows outside the training set, machine-learned augmentations are often
detrimental [26, 27], leading to unphysical results. This lack of robustness
prevents the use of machine-learning models in an engineering setting, which,
in turn, undermines the potential impact of endeavors aimed at enhancing
learning efficiency [24, 28] and capturing the elusive flow physics absent in
conventional models [25, 29, 30].

Understanding why machine-learning models do not generalize is crucial.
Exploring the answer to the question can help identify the underlying issues
and guide the development of machine-learning models. While a definitive
answer to this question is challenging, insights can be gained from the
discussions in Refs. [26, 27, 31, 32, 33, 34]. The basic calibrations of
empirical RANS models, such as decaying isotropic homogeneous turbulence,
zero-pressure-gradient flat-plate boundary layers, and free shear flows, are
the fundamental building blocks that make up a turbulent flow. These basic
calibrations enable empirical RANS models to generalize beyond their
calibration datasets. Based on this observation, one can surmise that machine-
learning models struggle to generalize precisely because they do not preserve
these basic calibrations. Incorporating the basic calibrations of RANS models
into machine-learning approaches could potentially grant them the
generalization abilities observed in empirical RANS models. An important
implication of this explanation is that solely focusing on capturing missing
physics in conventional empirical models, as attempted in Refs. [22, 35], will
not automatically grant machine-learning models the generalizability of
empirical models.

Attempting to incorporate all basic calibrations of empirical RANS models into
machine-learning models is excessively ambitious and unnecessary. Among the
basic calibrations of conventional empirical models, the law of the wall
(LoW), i.e., the scaling of the mean flow as a function of the distance from
the wall-normal coordinate in the viscous and logarithmic layer, is regarded
by many as the most important [31, 36]. Spalart argued that “the law of the
wall is the most useful and trusted part of our knowledge of turbulence in the
field of turbulence modeling” [36]. Menter and company further emphasized
that: while certain aspects such as decaying homogeneous isotropic turbulence
and free-shear flows may be disregarded, preserving the law of the wall is
imperative [31]. Therefore, in order to enable generalization, it is essential
to incorporate the law of the wall as an integral component of machine-
learning models. This work focuses precisely on the incorporation of the law
of the wall into machine-learning models.

We focus on three well-established modeling frameworks: tensor basis neural
networks (TBNN) introduced in Ref. [22], physics-informed machine learning
(PIML) proposed in Ref. [37], and field inversion and machine learning first
formulated in Ref. [23]. Below, we provide a summary of these three
frameworks, with a more comprehensive discussion deferred to Section 3.2. The
objective of TBNN is to model the anisotropic part of the Reynolds stress,
denoted as $b_{ij}$. Ling et al. [22] invoked the closed-form expression for
$b_{ij}$ in Ref. [38], which expresses $b_{ij}$ as a function of the strain
rate tensor $S_{ij}$ and the rotational rate tensor $\Omega_{ij}$. A deep
neural network is employed to learn the coefficients in the expression as a
function of the invariants of $S_{ij}$ and $\Omega_{ij}$. Consequently, TBNN
exhibits Galilean invariance. PIML, on the other hand, focuses on modeling the
discrepancy between the Reynolds stress predicted by RANS and that observed in
DNS. This is achieved through a random forest regressor, which takes as inputs
all locally available flow information. Notably, PIML also maintains Galilean
invariance when evolving invariant inputs only. The objective of FIML is to
augment the auxiliary equation in RANS models such that the resulting model
yields more accurate predictions of some given quantity of interest. The
method consists of field inversion followed by machine learning. During the
field inversion step, an augmentation is learned by minimizing a cost function
through model-consistent training. During the machine learning step, a neural
network is trained to relate the augmentation obtained from the field
inversion step to the flow variables.

TBNN, PIML, and FIML have had further developments since the early
publications [22, 37, 23]. Here, we review some of these developments. Parish
et al. [39] argued that instead of modeling the anisotropic part of the
Reynolds stress, modeling the error in the Reynolds stress anisotropy is a
more reasonable choice for TBNN. Xu et al. and Xie et al. [40, 41] extended
the framework of TBNN to LES modeling. Wu et al. [35] explored velocity
propagation for PIML. Zhang et al. [42] utilized the ensemble Kalman method to
handle sparse data for PIML. Singh et al. [24] incorporated adjoint-based
optimization into FIML, which significantly improved the training efficiency.

In addition to the aforementioned frameworks, other machine-learning
frameworks for RANS modeling have also emerged. A notable example is the
symbolic regression method in Refs. [43, 44, 45]. There, an evolutionary
algorithm drives algebraic forms of the Reynolds stress anisotropic tensor.
The method is appealing because it gives mathematically interpretable forms.
Yet another framework is progressive machine learning [32, 46]. The method
aims to preserve the basic calibrations of a baseline model. Applying this
framework, Bin et al. obtained a data-enabled recalibration of the SA model,
maintaining the basic calibrations (e.g., flat plate boundary layer and free
shear flows) while improving the model’s predictions of flow separation.
Considering the multitude of existing methods, implementing all of them would
be extremely time-consuming and may not provide substantial additional
insights at this stage. Therefore, a comparative study of methods other than
TBNN, PIML, and FIML is left for future investigation.

Besides the ML frameworks, we must also pick baseline models. This study
focuses on the one-equation Spalart-Allmaras (SA) model [9], the two-equation
Wilcox $k$-$\omega$ model [10], and the SSG seven-equation full Reynolds
stress model (FRSM) [12]. Both the SA model and the Wilcox $k$-$\omega$ model
are extensively used in the industry and are available in many software codes
[47, 48]. The Speziale-Sarkar-Gatski (SSG)-FRSM, while less commonly used, is
implemented in STARCCM+, OpenFOAM, and ANSYS Fluent. We postpone further
details of these models to section 3.1.

The remainder of this paper is organized as follows. Details of the baseline
RANS models and the machine-learning frameworks are presented in Sec. 3. The
results are shown in Sec. 4. Last, we conclude in section 5

##  3 Methodology

We present the details of the baseline RANS models in section 3.1 and the
details of the machine-learning frameworks in section 3.2. The details of the
training and testing data are presented in section 3.3.

###  3.1 Baseline RANS models

####  3.1.1 Spalart-Allmaras model

The Spalart-Allmaras (SA) model in Ref. [9] is a one-equation model. It solves
a transport equation for $\tilde{\nu}=\nu_{t}/f_{\nu 1}$, where $\nu_{t}$ is
the eddy viscosity, and $f_{\nu 1}$ is a function of the dimensionless
parameter $\chi\equiv\tilde{\nu}/\nu$. The transport equation for
$\tilde{\nu}$ is given by:

$\frac{\rm D\tilde{\nu}}{{\rm
D}t}=c_{b1}\tilde{S}\tilde{\nu}-c_{w1}f_{w}\left(\frac{\tilde{\nu}}{d}\right)^{2}+\frac{1}{\sigma}\left[\frac{\partial}{\partial
x_{j}}\left((\nu+\tilde{\nu})\frac{\partial\tilde{\nu}}{\partial
x_{j}}\right)+c_{b2}\frac{\partial\tilde{\nu}}{\partial
x_{j}}\frac{\partial\tilde{\nu}}{\partial x_{j}}\right],$ (1)

where ${\rm D}/{\rm D}t$ is the material derivative,
$\tilde{S}=\Omega+\tilde{\nu}/(\kappa^{2}d^{2})f_{\nu 2}$, $\kappa=0.41$ is
the von Kármán constant, $\Omega$ is the vorticity magnitude,
$c_{w1}=c_{b1}/\kappa^{2}+(1+c_{b2})/\sigma$, and $d$ is the distance to the
closest wall. The functions $f_{\nu 1}$, $f_{\nu 2}$ are

$f_{\nu 1}=\frac{\chi^{3}}{\chi^{3}+c_{\nu 1}^{3}},~{}~{}f_{\nu
2}=1-\frac{\chi}{1+\chi f_{\nu 1}};$ (2)

The function $f_{w}$, which depends on
$r=\tilde{\nu}/(\tilde{S}\kappa^{2}d^{2})$, is given by:

$f_{w}(r)=g\left(\frac{1+c_{w3}^{6}}{g^{6}+c_{w3}^{6}}\right)^{1/6},~{}~{}g=r+c_{w2}(r^{6}-r);$
(3)

The six model constants are: $\sigma=2/3$, $c_{b1}=0.1355$, $c_{b2}=0,622$,
$c_{w2}=0.3$, $c_{w3}=2$, and $c_{\nu 1}=7.1$.

####  3.1.2 Wilcox $k$-$\omega$ model

The Wilcox $k$-$\omega$ model in [10] is a two-equation model. It solves two
equations for $k$ and $\omega$:

$\begin{split}\frac{{\rm D}k}{{\rm D}t}&=\tau_{ij}\frac{\partial
u_{i}}{\partial x_{j}}-\theta^{*}\omega k+\frac{\partial}{\partial
x_{j}}\left[\left(\nu+\sigma^{*}\nu_{t}\right)\frac{\partial k}{\partial
x_{j}}\right],\\\ \frac{\rm D\omega}{{\rm
D}t}&=\frac{\gamma\omega}{k}\tau_{ij}\frac{\partial u_{i}}{\partial
x_{j}}-\theta\omega^{2}+\frac{\partial}{\partial
x_{j}}\left[\left(\nu+\sigma\nu_{t}\right)\frac{\partial\omega}{\partial
x_{j}}\right],\end{split}$ (4)

with

$\nu_{t}=k/\omega,$ (5)

and

$R_{ij}=2\nu_{t}\left(S_{ij}-\frac{1}{3}\frac{\partial u_{k}}{\partial
x_{k}}\right)-\frac{2}{3}k\delta_{ij},~{}~{}S_{ij}=\frac{1}{2}\left(\frac{\partial
u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),$ (6)

where $\delta_{ij}$ is the Kronecker delta. The five model constants are:
$\gamma=5/9$, $\theta^{*}=0.09$, $\theta=0.075$, $\sigma^{*}=0.5$, and
$\sigma=0.5$.

####  3.1.3 Full Reynolds Stress Model

FRSM solves the transport equations of the Reynolds stresses

$\centering\frac{{\rm D}R_{ij}}{{\rm
D}t}=\mathcal{P}_{ij}-\epsilon_{ij}+\mathcal{D}_{ij}+\Pi_{ij},\@add@centering$
(7)

where $\mathcal{P}_{ij}$ is the production tensor defined as
$\mathcal{P}_{ij}=-R_{ik}\partial{\overline{u}_{j}}/\partial{x_{k}}-R_{jk}\partial{\overline{u}_{i}}/\partial{x_{k}}$,
$\epsilon_{ij}$ is the dissipation tensor, $\mathcal{D}_{ij}$ is the diffusion
tensor, and $\Pi_{ij}$ is the pressure-strain correlation. The last three
terms on the right-hand side are unclosed. The dissipation tensor is modeled
as [13]

$\epsilon_{ij}=\frac{2}{3}\epsilon\delta_{ij}.$ (8)

The dissipation rate, $\epsilon$, is obtained by solving the following
transport equation

$\frac{{\rm D}\epsilon}{{\rm D}t}=C_{\epsilon
1}\frac{\mathcal{P}\epsilon}{k}-C_{\epsilon
2}\frac{\epsilon^{2}}{k}+\frac{\partial}{\partial
x_{i}}\left[\left(\nu\delta_{ij}+C_{\epsilon}\frac{k}{\epsilon}\langle
u_{i}u_{j}\rangle\right)\frac{\partial\epsilon}{\partial x_{j}}\right],$ (9)

where $\mathcal{P}$ is the production of TKE defined as $\langle
u_{i}u_{j}\rangle\partial u_{i}/\partial x_{j}$; $C_{\epsilon 1}$,
$C_{\epsilon 2}$ and $C_{\epsilon}$ are model constants. The diffusion term
and pressure-strain correlation term are modeled via the Daly-Harlow (DH)
gradient diffusion model [49] and the SSG model [12]:

$\mathcal{D}_{ij}=\frac{\partial}{\partial
x_{k}}\left(\nu\delta_{kl}+C_{DH}\frac{kR_{kl}}{\epsilon}\right)\frac{\partial
R_{ij}}{\partial x_{l}}$ (10)

and

$\displaystyle\Pi_{ij}$
$\displaystyle=-C_{1}\epsilon+C_{1}^{*}\mathcal{P}_{kk}b_{ij}+C_{2}\epsilon\left(b_{ik}b_{kj}-\frac{1}{3}b_{mn}b_{mn}\delta_{ij}\right)+\left(C_{3}-C_{3}^{*}\sqrt{b_{kl}b_{kl}}\right)kS_{ij}$
(11)
$\displaystyle+C_{4}k\left(b_{ik}S_{jk}+b_{jk}S_{ik}-\frac{2}{3}b_{mn}S_{mn}\delta_{ij}\right)+C_{5}k\left(b_{ik}W_{jk}+b_{jk}W_{ik}\right)$

The model constants are:

$\centering\begin{split}C_{\epsilon 1}=1.44,C_{\epsilon
2}=1.83,&C_{\epsilon}=0.15,C_{DH}=0.22\\\
C_{1}=3.4,C_{1}^{*}=1.8,C_{2}=4.2,C_{3}=\frac{4}{5},&C_{3}^{*}=1.3,C_{4}=1.25,C_{5}=0.4.\end{split}\@add@centering$
(12)

Table 2 summarizes the relevant information of the baseline models. Whether
the model is an eddy viscosity model and the number of auxiliary transport
equations are defining features of an empirical RANS model. Whether the model
contains the turbulent kinetic energy and dissipation rate information
determines whether TBNN and PIML are applicable to the baseline model and
therefore are relevant to this study. A wall treatment ensures that the model
captures the behavior of the mean flow in the viscous sublayer and the buffer
layer. For the three models considered here, a wall treatment exists in both
the SA model and the Wilcox $k$-$\omega$ model; there is a wall treatment for
FRSM in the literature [13], but this is not implemented in the present code.

Table 2:  Details of the RANS models. Here, Y and N stand for Yes and No. | Eddy viscosity | Equations | Contains k info | Constains $\epsilon$ info | Wall treatment
---|---|---|---|---|---
SA | Y | 1 | N | N | Y
Wilcox $k$-$\Omega$ | Y | 2 | Y | Y | Y
SSG FRSM | N | 7 | Y | Y | N

We implement the SA and Wilcox $k$-$\omega$ models in Python. A second-order
central difference is employed in space, and an implicit Euler method is used
for time integration. We validate these two implementations against the ones
in OpenFOAM regarding plane channel flow (not shown for brevity). An
implementation of the SSG FRSM model exists in OpenFOAM, and we make use of
that implementation.

###  3.2 Machine-learning frameworks

####  3.2.1 TBNN

Figure 1: Schematic of (a) TBNN and (b) PIML.

TBNN, as introduced in Ref. [22], predicts the Reynolds stress anisotropy
tensor, denoted as $b_{ij}=R_{ij}/2k-(1/3)\delta_{ij}$, following the
formulation outlined in Ref. [38]:

$b_{ij}=\sum_{k=0}^{10}c_{k}(\lambda_{1},\lambda_{2},...,\lambda_{5})T^{(k)},$
(13)

where

$\centering\begin{split}T^{(0)}=\frac{1}{3}\boldsymbol{I_{3}}-\frac{1}{2}\boldsymbol{I_{2}},~{}T^{(1)}=\boldsymbol{S},~{}&T^{(2)}=\boldsymbol{SR}-\boldsymbol{RS},~{}T^{(3)}=\boldsymbol{S^{2}}-\frac{1}{3}\boldsymbol{I}\cdot
Tr(\boldsymbol{S^{2}}),~{}T^{(4)}=\boldsymbol{R^{2}}-\frac{1}{3}\boldsymbol{I}\cdot
Tr(\boldsymbol{R^{2}}),\\\
~{}T^{(5)}=\boldsymbol{RS^{2}}-\boldsymbol{S^{2}R},&~{}T^{(6)}=\boldsymbol{R^{2}S}+\boldsymbol{SR^{2}}-\frac{2}{3}\boldsymbol{I}\cdot
Tr(\boldsymbol{SR^{2}}),~{}T^{(7)}=\boldsymbol{RSR^{2}}-\boldsymbol{R^{2}SR},\\\
~{}T^{(8)}=\boldsymbol{SRS^{2}}-\boldsymbol{S^{2}RS},&~{}T^{(9)}=\boldsymbol{R^{2}S^{2}+S^{2}R^{2}}-\frac{2}{3}\boldsymbol{I}\cdot
Tr(\boldsymbol{S^{2}R^{2}}),~{}T^{(10)}=\boldsymbol{RS^{2}R^{2}}-\boldsymbol{R^{2}S^{2}R},\end{split}\@add@centering$
(14)

and

$\centering\lambda_{1}=Tr(\boldsymbol{S^{2}}),~{}\lambda_{2}=Tr(\boldsymbol{R^{2}}),~{}\lambda_{3}=Tr(\boldsymbol{S^{3}}),~{}\lambda_{4}=Tr(\boldsymbol{R^{2}S}),~{}\lambda_{5}=Tr(\boldsymbol{R^{2}S^{2}}).\@add@centering$
(15)

Here, $T^{(0)}$ is reversed for two-dimensional flow only, $\boldsymbol{S}$ is
the strain rate tensor, $\boldsymbol{R}$ is the rotation tensor,
$\boldsymbol{I_{3}}$ is the identity matrix, and $\boldsymbol{I_{2}}$ is also
the identity matrix but $I_{i3}=0$. For two-dimensional flow, only $T_{0}$,
$T_{1}$, $T_{2}$, $\lambda_{1}$, and $\lambda_{2}$ are non-zero. Figure 1 (a)
shows a schematic of TBNN. The network comprises two input layers: the
invariant input layer and the tensor input layer. While the invariant input
layer can be expanded to include parameters like Reynolds number and distance
from the wall [40, 50], such expansion is not explored in this study. We
adhere to the setup proposed in Ref. [22] when configuring our neural network.
The network is a feedforward neural network. It has 4 hidden layers with 30
neurons in each layer. The activation function is leaky-ReLU. The cost
function is the root-mean-square-error (RMSE) in $b_{ij}$.

####  3.2.2 PIML

Physics-informed inputs are fed to a regression tool to predict the error in
$R_{ij}$. The regression tool can take the form of either a neural network or
a random forest. The training process relies on benchmark DNS data, and in the
work of Wang et al. [37], there are approximately O(10) inputs. However, for
one-dimensional flows, such as channels and temporally evolving shear layers,
the number of inputs reduces to four. The four inputs $\hat{q}_{i}$ and their
normalizations $q_{i}^{*}$ are tabulated in Table 3. Following Ref. [37], the
inputs to the neural network are
$q_{\beta}=\hat{q}_{\beta}/(|\hat{q}_{\beta}|+|q_{\beta}^{*}|)$ to prevent
mathematical singularity.

Concerning the output, PIML in Ref. [37] predicts the magnitude $(\Delta k)$,
shape $(\Delta\xi,\Delta\eta)$, and orientation
$(\Delta\phi_{1},\Delta\phi_{2},\Delta\phi_{3})$ of the object tensor. For
plane channel flow, the output reduces to $(\Delta k,\Delta b_{12})$. The
network is a feedforward neural network. It has 4 hidden layers with 30
neurons in each layer. The activation function is the sigmoid function.

Table 3:  Inputs of PIML. Feature$(q_{i})$ | Description | Raw input $\hat{q}_{i}$ | Normalization factor $q_{i}^{*}$
---|---|---|---
$q_{1}$ | Turbulent kinetic energy | $k$ | $\frac{1}{2}U_{i}U_{i}$
$q_{2}$ | Wall-distance based Reynolds number | $\sqrt{k}d/50\nu$ | NA
$q_{3}$ | Ratio of turbulent time scale to mean strain time scale | $k/\epsilon$ | $1/||\boldsymbol{S}||$
$q_{4}$ | Ratio of total to normal Reynolds stresses | $\lVert\langle u_{i}u_{j}\rangle\lVert$ | $k$

####  3.2.3 FIML

Figure 2: Schematic of FIML.

FIML consists of two steps: field inversion followed by machine learning, as
sketched in Fig. 2. During the field inversion step, an augmentation of the
baseline RANS model is learned. In this study, the augmentation manifests as a
multiplier in front of $f_{v1}$ in the SA model, a multiplier in front of the
production term in the $k$ equation in the Wilcox $k$-$\omega$ model, and a
multiplier in front of the production term in the $\epsilon$ equation in the
SSG FRSM, respectively. The field inversion involves minimizing the following
cost function

$\centering
J_{c}(\beta)=\lVert\phi_{d}-\phi_{r}(\beta)\lVert_{2}^{2}-\lambda\lVert\beta-1.0\lVert_{2}^{2},\@add@centering$
(16)

where $\phi_{d}$ is the truth and $\phi_{r}$ is the RANS prediction. $\phi$
can be any variable. For plane channel flow, $\phi=U$ is the mean velocity
profile. The second term in Eq. (16) is a regularization term. Here, we set
$\lambda=0.001$. To minimize $J_{c}$, we resort to gradient descent

$\centering\beta^{i+1}=\beta^{i}-\alpha\frac{\partial
J}{\partial\beta_{i}},\@add@centering$ (17)

where gradient $\partial{J}/\partial{\beta_{i}}$ can be computed via either an
adjoint method or a finite difference method. Here, we employ the adjoint
method for the SA and the Wilcox $k$-$\omega$ models and the finite difference
method for the SSG FRSM. Next, a regression tool is used to relate flow
features to the augmentation $\beta$. We use a neural network as our
regression tool. The flow features we use for channel flow are:

$f_{SA}=\left\\{\frac{\nu}{\nu_{t}+\nu}\right\\},~{}f_{k\rm{-}\omega}=\left\\{\frac{\nu}{\nu_{t}+\nu},-\frac{\nu}{dp/dx}\frac{dU}{dy}\right\\},~{}f_{\rm
FRSM}=\left\\{\frac{k/\epsilon}{k/\epsilon+\lVert\boldsymbol{S}\lVert},\frac{\sqrt{k}d}{\nu}\right\\},$
(18)

for the SA, the Wilcox $k$-$\omega$, and the SSG FRSM, respectively [26].

Table 4 summarizes the relevant information of the ML frameworks. A priori
training is the method where one trains an intermediate variable to match
high-fidelity data. Model consistent training is the method where one trains
directly to match the quantity of interest (QoI). Duraisamy [51] argued that
model-consistent training is preferred over a prior training because the
correct prediction of an intermediate variable does not necessarily guarantee
the correct prediction of the end QoI. “RANS interface” identifies where a
baseline model interfaces with the machine-learning framework. When the
interface is the Reynolds stress, the baseline model is intact, and whether
one conducts velocity propagation becomes relevant. When the interface is the
transport equation, the machine-learning augmentation becomes an integral part
of the model, and velocity propagation is implied.

Table 4: Details of the ML frameworks.  | Target | Velocity propagation | RANS Interface | Training
---|---|---|---|---
TBNN | Reynolds stress | N | Reynolds stress | a priori
PIML | Error in Reynolds stress | N | Reynolds stress | a priori
FIML | Any QoI | N/A | Transport equation | Model consistent

For all ML frameworks, training is conducted in Python. The training involved
an Adam optimizer. All training does not stop until the cost function cannot
be further reduced in another 300 epochs.

###  3.3 Data

We make use of the plane channel flow DNS data in the Johns Hopkins Turbulence
Database and the UP Madrid database [52, 53, 54] and the mixing layer DNS data
in Ref. [55].

The channel flow data are at the friction Reynolds numbers
$\text{Re}_{\tau}=180$, $395$, $590$, $1000$, $2000$, and $5200$. The details
of the DNSs are tabulated in Table 5. The mean flow data should be highly
accurate according to Refs. [2, 56]. For this study, data at $Re_{\tau}=$180,
590, 1000, and 2000 are used for training; for testing, we use data at
$Re_{\tau}=180$, one of the training flow conditions in the training set,
$Re_{\tau}=395$, a condition that can be obtained by interpolating inside the
training set, and $Re_{\tau}=5200$, a condition that must be obtained through
extrapolation. Interpolation and extrapolation here are with respect to the
flow parameter space, i.e., the Reynolds number, rather than the input space.

Table 5: Details of the channel flow DNS data. Here, the asterisk and plus superscripts indicate normalization by the channel half-width $\delta$ and the viscous length scale $\nu/u_{\tau}$, respectively. Further, the subscript $w$ refers to values at the wall. | Reynolds Number:
---
| Domain:
---
$L_{x}^{*},L_{y}^{*},L_{z}^{*}$
| Grid:
---
$N_{x},N_{y},N_{z}$
| Resolution:
---
$\Delta x^{+},\Delta y_{w}^{+},\Delta z^{+}$
$\text{Re}_{\tau}=180$ | $2\pi,2,4\pi/3$ | $128,129,128$ | $17.7,0.05,5.9$
$\text{Re}_{\tau}=395$ | $4\pi,2,\pi$ | $256,193,192$ | $10.0,0.03,6.5$
$\text{Re}_{\tau}=590$ | $4\pi,2,\pi$ | $384,257,384$ | $9.7,0.04,7.2$
$\text{Re}_{\tau}=1000$ | $8\pi,2,3\pi$ | $2048,512,1536$ | $10.3,0.02,4.6$
$\text{Re}_{\tau}=2000$ | $8\pi,2,3\pi$ | $4096,768,3072$ | $8.2,0.002,4.1$
$\text{Re}_{\tau}=5200$ | $8\pi,2,3\pi$ | $10240,1536,7680$ | $12.7,0.07,6.4$

The temporal shear layer data are from Ref. [55]. Figure 3 is a schematic of
the flow. Two miscible fluids with velocities of equal magnitudes but opposite
signs are brought together. The flow is homogeneous in the two horizontal
directions and evolves in time. The initial velocity profile is given by

$U=-\frac{\Delta
U}{2}\tanh\left(\frac{2(y-y_{0})}{\delta_{(\omega,0)}}\right),$ (19)

where $\Delta U$ represents the difference in the velocity between the two
fluids, $y_{0}=0$ is at the center of the domain, and $\delta_{(\omega,0)}$ is
the initial vorticity thickness defined as $\delta_{\omega}\equiv\Delta
U/(d\langle u\rangle/dy)_{max}$. The initial Reynolds number of the flow is
$Re=\Delta U\delta_{(\omega,0)}/\nu=640$. The momentum thickness is defined as
$\delta_{\theta}\equiv\int_{-\infty}^{\infty}(1/4-\langle u\rangle^{2})dy$.
The flow evolves temporally. We reserve the flow at $t\Delta
U/\delta_{(\omega,0)}=$ 5, 13, 15.5, 31, 49 for training and $t\Delta
U/\delta_{(\omega,0)}=$ 46.5, 57 for testing. Further details about the
computational setup of the case can be found in Refs. [57, 55] and are not
repeated here for brevity.

Figure 3: Schematic of a temporally-evolving shear layer.

##  4 Results

We train against channel flow data and present the results in Sec. 4.1. In
Sec. 4.2, we discuss the TBNN and PIML results. In Sec. 4.3, we apply the FIML
model to other flows to test its generalizability. Last, we train against both
channel flow and mixing layer flow data to include another basic calibration.
The results are presented in section 4.4.

Before we proceed, we would like to remark on the presentation of the results.
For FIML, we present only the mean flow results, as the method does not
concern any intermediate variables. For TBNN and PIML, we first present the
target variable of machine learning and only examine the mean flow if ML
yields good predictions for the target variables. It is not advisable to skip
the machine learning variable(s): if there were errors in the mean flow, it
would not be clear whether the culprit is the ML or velocity propagation.

###  4.1 Channel flow

The results are organized into 3 sections for TBNN, PIML, and FIML,
respectively. Again, we train against data at $Re_{\tau}=180$, 590, 1000, and
2000. Testing results will be shown at three Reynolds numbers:
$Re_{\tau}=180$, one of the training conditions, $Re_{\tau}=395$, a condition
that can be interpolated from the training set, and $Re_{\tau}=5200$, a
condition that must be extrapolated. There is an extended logarithmic layer at
$Re_{\tau}=5200$, allowing us to test if the model has learned the LoW.

####  4.1.1 TBNN

The output $b_{ij}$ requires information on turbulent kinetic energy, which is
not available in the SA model. As a result, TBNN, at least in its form as
presented in Refs. [22, 23], cannot be applied to the SA model.

Figure 4 displays the TBNN-predicted $b_{12}$ at $Re_{\tau}=180$, with results
from two baseline models included for comparison. Notably, the baseline SSG
FRSM model and the baseline Wilcox $k$-$\omega$ model tend to over-predict
$\left|b_{12}\right|$ in the wall layer. TBNN’s augmentation to the Wilcox
$k$-$\omega$ model maintains the overall behavior of $\left|b_{12}\right|$ but
brings the profile closer to the DNS. The augmentation to the SSG FRSM alters
the overall behavior of $b_{12}$ and results in a closer agreement between
RANS and DNS at all $y$ locations. However, an unphysical peak emerges in the
viscous sublayer, which will be discussed in Sec. 4.2. Figure 5 presents
results at $Re_{\tau}=395$ and $5200$, showing similarities to those in Fig. 4
at $Re_{\tau}=180$. The observed oscillations in the profiles at
$Re_{\tau}=5200$ are expected, as noted in [58].

Figure 4: Profiles of $b_{12}$ at $Re_{\tau}=180$. The black solid line is the
DNS. The colored solid lines are the results of the baseline RANS. The dot-
dashed lines are TBNN results. Figure 5: Same as Fig. 4 but at (a)
$Re_{\tau}=395$, (b) $Re_{\tau}=5200$. Again, the colored solid lines are the
results of the baseline RANS models; the dot-dashed lines are TBNN results.

Multiplying $b_{12}$ by $k$ and integrating the mean momentum equation yields
the mean velocity profile:

$U=\int_{0}^{y}\frac{dU}{dy}dy=\int_{0}^{y}\frac{1}{\nu}\left[(\tau_{w}/\rho)(1-y/\delta)+R_{12}\right]~{}dy,$
(20)

where

$R_{12}=2b_{12}k_{\text{Baseline RANS}}.$ (21)

Figure 6 illustrates the mean flow results. Although TBNN leads to improved
$b_{12}$ predictions, these improvements do not necessarily translate into
enhancements in the mean flow, primarily due to disparate values of $k$ in
RANS and DNS, as depicted in Fig. 7.

Figure 6: Velocity profiles by integrating TBNN predicted $R_{12}$ through Eqs
20 and 21. Colored solid lines are the results of the baseline RANS model.
Colored dot-dashed lines are the results of TBNN. (a) $Re_{\tau}=180$, (b)
$Re_{\tau}=395$, (c) $Re_{\tau}=5200$. Figure 7: Turbulent kinetic energy in
DNS and the two baseline RANS models. Here, different colors are used to
distinguish different Reynolds numbers. Solid lines are for DNS, dashed lines
are for $k$-$\omega$, and dotted lines are for FRSM.

####  4.1.2 PIML

Figure 8 illustrates the profiles of $R_{12}$ at $Re_{\tau}=180$. The baseline
$k$-$\omega$ result is already very close to the DNS, while the SSG FRSM
result overestimates the Reynolds shear stress. PIML’s augmentation improves
the $R_{12}$ prediction for both SSG FRSM and Wilcox $k$-$\omega$, resulting
in nearly perfect $R_{12}$ profiles for both models.

Figure 9 displays the results at $Re_{\tau}=395$ and 5200. Although it appears
that PIML’s augmentation enhances the predictions of $R_{12}$ at both Reynolds
numbers, closer inspection in Fig. 10 reveals some inaccuracies. Specifically,
PIML’s prediction of $R_{12}$ is slightly positive in the viscous sublayer at
$Re_{\tau}=5200$, accompanied by unphysical undulations around $y\approx 100$.
Furthermore, $dU^{+}/dy^{+}=1-y/\delta+\left<u^{\prime}v^{\prime}\right>^{+}$
is slightly negative in the wake layer. These unphysical behaviors, as we will
discuss shortly, significantly impact the mean flow.

Figure 8: The Reynolds shear stress $\langle uv\rangle$ at $Re_{\tau}=180$.
Colored solid lines are the results of the baseline RANS model. Colored dot-
dashed lines are the results of PIML. Figure 9: Same as Fig. 8 but at (a)
$Re_{\tau}=395$ and (b) $Re_{\tau}=5200$. Figure 10: A zoom-in view of PIML’s
predictions at $Re_{\tau}=5200$. (a) for $0.2\leq y^{+}\leq 5$, (b) for
$10\leq y^{+}\leq 500$, (c)$\left[(\tau_{w}/\rho)(1-y/\delta)+\langle
uv\rangle^{+}\right]$, i.e., $dU^{+}/dy^{+}$ for $800\leq y^{+}\leq 5200$ .

Integrating the mean momentum equation as per Eq. (20) gives the mean flow.
Figure 11 presents the results. In comparison to the TBNN results in Fig. 6,
the PIML results exhibit more favorable outcomes here. PIML’s augmentation
improves the mean flow predictions of the two baseline models at both
$Re_{\tau}=180$ and 395. For SSG FRSM, PIML’s augmentation acts as a wall
treatment. However, PIML’s augmentation proves detrimental when extrapolating,
as depicted in Fig. 11 (c), due to the errors shown in Fig. 10. This is
consistent with the findings in Ref. [59], where it was demonstrated that
velocity propagation incurs significant errors.

Figure 11: Velocity profiles by integrating PIML predicted $R_{12}$ according
to Eq. 20. Colored solid lines are the results of the baseline RANS models.
Colored dot-dashed lines are the results of PIML. (a) $Re_{\tau}=180$, (b)
$Re_{\tau}=395$, (c) $Re_{\tau}=5200$.

####  4.1.3 FIML

FIML is applicable to all baseline models, and the training specifically
addresses the mean flow itself. Figure 12 presents the mean flow results at
$Re_{\tau}=180$. The baseline SA and Wilcox $k$-$\omega$ models already
provide good mean flow predictions in plane channels. For these two models,
FIML’s augmentations are largely neutral, with slight improvements observed in
the buffer layer. The baseline SSG FRSM lacks wall damping. FIML’s
augmentation leads to a different log law slope. While FIML’s augmentation
brings the RANS profile closer to the DNS profile, the solution does not
achieve the correct layered structure. Figure 13 depicts the results at
$Re_{\tau}=395$ and 5200, showing similarities to those in Fig. 12. For SA and
Wilcox $k$-$\omega$, where there is already a wall treatment, FIML’s
augmentation brings the profile closer to the DNS in the buffer layer. In the
case of SSG FRSM, lacking a wall treatment, FIML changes the log layer slope,
bringing the RANS profiles closer to the DNS. We see that FIML’s augmentations
are non-detrimental. This is different from the results in [26], a discussion
of which is postponed to section 4.2.

Figure 12: Velocity profiles at $Re_{\tau}=180$. Here, colored solid lines
represent baseline results and colored dot-dashed lines represent FIML
results. Figure 13: Same as Fig. 12 but at (a) $Re_{\tau}=395$, (b)
$Re_{\tau}=5200$.

Table 6 summarizes the results in Sec. 4.1. The data-based augmentations due
to TBNN prove to be detrimental both inside and outside the training set. The
augmentations due to PIML are beneficial inside the training set but are
detrimental outside it. The augmentations due to FIML are non-detrimental
inside or outside the training set. However, we see that when the baseline
model lacks the right physics, FIML does not necessarily learn that physics.

Table 6: A summary of channel flow results | SA | Wilcox $k$-$\omega$ | SSG FRSM
---|---|---|---
TBNN | N/A | Improved $b_{12}$, detrimental to $U$ inside and outside the training set
PIML | N/A | Improved $R_{12}$, beneficial and detrimental to $U$ inside and outside the training set
FIML | Neutral in log layer, beneficial in the buffer layer | Beneficial, but did not capture the right physics

###  4.2 Further discussion of TBNN and PIML

For TBNN, the network’s inputs consist of invariants derived from the strain
rate and rotation tensor. In the context of plane channels, the strain rate
tensor and the rotation tensor have only one non-zero component, which is
$dU/dy$. Therefore, any non-zero input to the network is a result of a non-
zero $dU/dy$. In Fig. 14, we illustrate $k/\epsilon\partial{U}/\partial{y}$ in
the baseline FRSM and Wilcox $k$-$\omega$ model. It is evident that
$k/\epsilon\partial{U}/\partial{y}$ increases with $y^{+}$ in the inner layer
and decreases with $y^{+}$ in the outer layer. Consequently, a given
$k/\epsilon\partial{U}/\partial{y}$ value appears twice: once in the inner
layer and a second time in the outer layer. A direct consequence is that the
target variable of TBNN, $b_{12}$, is not a single-valued function of any
inputs derived from $k/\epsilon\partial{U}/\partial{y}$. This is clear from
Fig. 15, where we display $b_{12}$ as a function of
$k/\epsilon\partial{U}/\partial{y}$ at $Re_{\tau}=180$. The results at other
Reynolds numbers are similar and therefore are not shown here for brevity.
However, a neural network, particularly in the form presented in Ref. [22], is
inherently single-valued. When trained to predict $b_{12}$ as a function of
variables derived from $k/\epsilon S_{12}$, the neural network tends to
compromise between multiple $b_{12}$ values. The predictions of the trained
networks, as a function of $k/\epsilon dU/dy$, are demonstrated in Fig. 15.
For FRSM, the predicted $b_{12}$ exhibits non-monotonic behavior with respect
to $k/\epsilon dU/dy$. This non-monotonic behavior is responsible for the two
peaks observed in $b_{12}$ profiles in Figs. 4 and 5. The above explains why
TBNN is detrimental.

Figure 14: $k/\epsilon\partial{U}/\partial{y}$ due to (a) the baseline FRSM
and (b) the baseline Wilcox $k$-$\omega$ model. Figure 15: Desired $b_{12}$
and TBNN predicted $b_{12}$ as a function of
$k/\epsilon\partial{U}/\partial{y}$ for (a) FRSM and (b) $k$-$\omega$. We
indicate the direction of increasing $y^{+}$.

Next, we consider PIML. In the case of plane channel flow, four of the O(10)
inputs to the network are non-zero. Figure 16 visually represents these four
non-zero inputs. Notably, inputs $q_{1}$ and $q_{2}$ exhibit monotonic
behavior as functions of $y^{+}$. This characteristic allows a neural network
to accurately approximate the target variable as a function of either $q_{1}$,
$q_{2}$, or a combination of both. Consequently, the network is capable of
making improvements within the training dataset. However, since the Reynolds
stress must be integrated to obtain the mean flow, PIML is susceptible to
errors stemming from the integration process. As noted by Wu et al., depending
on the chosen integration method, there could be a substantial 35% error in
the mean flow [60]. The above explains why PIML is beneficial inside the
training dataset but detrimental outside.

Figure 16: Normalized inputs of PIML calculated from FRSM (a) $k$, (b)
$\frac{\sqrt{k}d}{50\nu}$, (c) $k/\epsilon$, (d)
$||\overline{u^{\prime}_{i}u^{\prime}_{j}}||$. The normalization process
follows that in Ref. [37] and given in Table 3.

###  4.3 Generalization of FIML

We implement FIML’s augmentation of the $k-\omega$ model in OpenFOAM [48] and
assess its performance in scenarios involving a mixing layer and a backward-
facing step. It is important to note that we do not retrain for these two
cases. The objective is to scrutinize whether FIML’s augmentation proves
detrimental outside the training dataset. The test in this subsection is
limited and more comprehensive tests fall outside the scope of this work. Our
methodology follows the steps outlined in Ref. [26]. We first allow the
baseline model to converge, obtaining the initial guess for the $\beta$ field
from the converged flow field. Subsequently, we let the augmented model
converge. The detailed setups for the two cases can be found in Ref. [61] and
are not reiterated here for brevity.

The results for the spatial mixing layer are depicted in Figs. 17 and 18.
Figure 17 displays the velocity contour of both the baseline model and the
FIML-augmented model. Meanwhile, Fig. 18 illustrates the velocity profiles at
various downstream locations, comparing the experimental data from Ref. [62],
the baseline Wilcox $k$-$\omega$ model, and the FIML-augmented model. The
contour plot in Fig. 17 shows no noticeable differences between the baseline
and FIML. Differences between the two are more clear from the line plots in
Fig. 18. The results of the baseline model are already very accurate. FIML’s
augmentation gives rise to enhanced mixing, which negatively impacts the
results. The effect is, however, not as concerning as in Ref. [26]. Figure 19
shows the results for the backward-facing step case. We plot the skin friction
coefficient $C_{f}$ and the pressure coefficient $C_{p}$ for the baseline
model and the FIML-augmented model. We see that FIML’s augmentation does not
affect the result. Most importantly, it is non-detrimental.

Figure 17: Velocity contours of the mixing-layer flow. (a) Baseline Wilcox
$k$-$\omega$ result, (b) FIML-augmented model result. Figure 18: Velocity
profiles in the shear layer at $x=$ (a) 200mm, (b) 650mm, and (c) 950mm.
Figure 19: (a) Skin friction coefficient $C_{f}$ and (b) pressure coefficient
$C_{p}$ of the backward-facing-step.

###  4.4 Incorporating other basic calibrations

In addition to the LoW, other basic calibrations of RANS models include
decaying isotropic turbulence, shear layers, etc. In this section, we build on
Sec. 4.1 and train against both channel flow and temporally-evolving shear
layers.

The shear flow is unsteady, applying FIML requires adjusting the RANS solution
at many time steps. This process can lead to ill-posedness [63], and therefore
FIML is not pursued here. Furthermore, since TBNN is detrimental when trained
against channel flow data, we do not pursue TBNN here either. In this section,
we will focus on PIML and the baseline $k$-$\omega$ model.

The results for the plane channel are similar to those presented in Sec. 4 and
are not shown here for brevity. We plot the Reynolds shear stress $\langle
uv\rangle$ due to the baseline model, the DNS, and PIML at $t=15.5$ and 57 in
Fig. 20. Here, time is normalized by $\Delta U$ and $\delta_{(\omega,0)}$. We
see that the baseline model underestimates and overestimates $\langle
uv\rangle$ at the early and late stages of the flow and that PIML yields no
improvement. In the following, we explain why there is no improvement. First,
we note that the flow is roughly self-similar, i.e., $U/\Delta U$ is a
function of $y/\delta_{\omega}$ only and not a function of time. This is clear
from Fig. 21. A direct consequence of self-similarity is that the locally
normalized PIML inputs, i.e., $q_{1}$, $q_{2}$, $q_{3}$, and $q_{4}$, do not
depend on time. Next, we consider the error in the baseline model. The error
distribution in the baseline model necessitates PIML to predict a positive and
a negative correction to the baseline model. This is not possible with inputs
that do not vary with time: improvements in the early stage will necessarily
result in degradation in the late stage and vice versa. As the error of the
baseline $k$-$\omega$ already gives a rather balanced error distribution, PIML
did not yield further improvements.

Figure 20: Reynolds stress profile of mixing layer at (a) $t=15.5$, (b)
$t=57$. The black solid line represents the DNS in Ref. [55], the colored
solid line is the result of the baseline $k$-$\omega$ model, and the colored
dot-dashed line is the prediction from PIML. Figure 21: Self-similarity of the
temporal mixing layer: (a) velocity profile at various time steps.

##  5 Concluding remarks

We apply TBNN, PIML, and FIML to the one-equation SA model, the two-equation
$k$-$\omega$ model, and the seven-equation FRSM. TBNN and PIML as presented in
Refs. [22, 37] directly modifies the Reynolds stress, and FIML as presented in
Ref. [23, 28] modifies terms in the auxiliary transport equations. The SA and
the $k$-$\omega$ models are eddy viscosity models, whereas FRSM is a Reynolds
stress model. We train against plane channel flows and shear layer flows. The
SA and the $k$-$\omega$ model already yield accurate mean flow estimates for
plane channel. The SSG FRSM, on the other hand, does not capture the mean flow
in the buffer layer.

TBNN and PIML, at least in the forms as presented in Refs. [22, 37], require
TKE and dissipation for normalization purposes. As such information is not
available in the SA model, TBNN and PIML are not applicable to the SA model.
We first train against plane channel flows at various Reynolds numbers.
Applying TBNN to the $k$-$\omega$ model and the SSG FRSM, we see unphysical
but nevertheless improved predictions of $b_{12}$. However, the improvements
in $b_{12}$ do not translate to the mean flow, and TBNN is detrimental both
inside and outside the training dataset. (This conclusion applies only to the
flows considered here.) Further analysis shows that TBNN yields degrading
results here because the target variable is not a single-valued function of
the inputs in this flow. Next, PIML yields improved $R_{12}$ predictions.
However, due to the errors associated with velocity propagation, PIML is
beneficial only at the Reynolds numbers inside the training dataset and is
detrimental outside. These conclusions are independent of the baseline model.

Training against plane channel flow data at various Reynolds numbers, FIML’s
augmentations are largely neutral in the SA and the $k$-$\omega$ model with
slight improvements in the buffer layer. Its augmentation to the SSG FRSM
leads to a closer agreement between the RANS result and the DNS data, but such
improvement is achieved via varying the von Kármán constant rather than a
layer-structure of the mean flow (i.e., sublayer, buffer layer, log layer).
The generalizability of FIML’s augmentation in the $k$-$\omega$ model is
tested in backward-facing step and mixing-layer scenarios. The augmentation
has essentially no impact on these two flows: we see slightly deteriorated
results for the mixing layer and no noticeable change in the results for the
backward-facing step.

The scope of the present study is limited to three ML frameworks, three
baseline RANS models, and two fluid flows. There have been further
developments since the early publications in Refs. [22, 37, 23] of the three
ML frameworks. Furthermore, there are other ML frameworks [43, 44, 45, 32,
64], other baseline RANS models, and other fluid flows. To enable the use of
machine learning models in fluids engineering, many more comparative studies
and validation & verification studies are yet needed.

## Acknowledgments

Li and Yang acknowledge ONR contract number N000142012315 and AFOSR award
number FA9550-23-1-0272. Bin acknowledges NNSFC grant number 91752202.

## References

  * Choi and Moin [2012] Choi, H., and Moin, P., “Grid-point requirements for large eddy simulation: Chapman’s estimates revisited,” _Phys. Fluids_ , Vol. 24, No. 1, 2012, p. 011702.
  * Yang et al. [2021] Yang, X. I. A., Hong, J., Lee, M., and Huang, X. L., “Grid resolution requirement for resolving rare and high intensity wall-shear stress events in direct numerical simulations,” _Phys. Rev. Fluids_ , Vol. 6, No. 5, 2021, p. 054603.
  * Li et al. [2022] Li, J.-Q. J., Yang, X. I. A., and Kunz, R. F., “Grid-point and time-step requirements for large-eddy simulation and Reynolds-averaged Navier–Stokes of stratified wakes,” _Phys. Fluids_ , Vol. 34, No. 11, 2022, p. 115125.
  * Pope [2000] Pope, S. B., _Turbulent flows_ , Cambridge university press, 2000.
  * Durbin [2018] Durbin, P. A., “Some recent developments in turbulence closure modeling,” _Annu. Rev. Fluid Mech._ , Vol. 50, 2018, pp. 77–103.
  * Kalitzin et al. [2005] Kalitzin, G., Medic, G., Iaccarino, G., and Durbin, P., “Near-wall behavior of RANS turbulence models and implications for wall functions,” _J Comput Phys_ , Vol. 204, No. 1, 2005, pp. 265–291.
  * Smith and Cebeci [1967] Smith, A., and Cebeci, T., “Numerical solution of the turbulent-boundary-layer equations,” Tech. rep., Douglas Aircraft Co., Long Beach, CA, Aircraft Div., 1967.
  * Baldwin and Lomax [1978] Baldwin, B., and Lomax, H., “Thin-layer approximation and algebraic model for separated turbulent flows,” _16th aerospace sciences meeting_ , 1978, p. 257.
  * Spalart and Allmaras [1992] Spalart, P., and Allmaras, S., “A one-equation turbulence model for aerodynamic flows,” _30th Aerospace Sciences Meeting and Exhibit_ , 1992, p. 439.
  * Wilcox [1988] Wilcox, D. C., “Reassessment of the scale-determining equation for advanced turbulence models,” _AIAA J._ , Vol. 26, No. 11, 1988, pp. 1299–1310.
  * Menter [1994] Menter, F. R., “Two-equation eddy-viscosity turbulence models for engineering applications,” _AIAA J._ , Vol. 32, No. 8, 1994, pp. 1598–1605.
  * Speziale et al. [1991] Speziale, C. G., Sarkar, S., and Gatski, T. B., “Modelling the pressure–strain correlation of turbulence: an invariant dynamical systems approach,” _J. Fluid Mech._ , Vol. 227, 1991, pp. 245–272.
  * Launder et al. [1975] Launder, B. E., Reece, G. J., and Rodi, W., “Progress in the development of a Reynolds-stress turbulence closure,” _J. Fluid Mech._ , Vol. 68, No. 3, 1975, pp. 537–566.
  * Chien [1982] Chien, K.-Y., “Predictions of channel and boundary-layer flows with a low-Reynolds-number turbulence model,” _AIAA J._ , Vol. 20, No. 1, 1982, pp. 33–38.
  * Hanjalić et al. [2004] Hanjalić, K., Popovac, M., and Hadžiabdić, M., “A robust near-wall elliptic-relaxation eddy-viscosity turbulence model for CFD,” _Int. J. Heat Fluid Flow_ , Vol. 25, No. 6, 2004, pp. 1047–1051.
  * Durbin [1995] Durbin, P. A., “Separated flow computations with the k-epsilon-v-squared model,” _AIAA J._ , Vol. 33, No. 4, 1995, pp. 659–664.
  * Slotnick et al. [2014] Slotnick, J. P., Khodadoust, A., Alonso, J., Darmofal, D., Gropp, W., Lurie, E., and Mavriplis, D. J., “CFD vision 2030 study: a path to revolutionary computational aerosciences,” , 2014.
  * Slotnick et al. [2011] Slotnick, J., Hannon, J., and Chaffin, M., “Overview of the 1st AIAA CFD high lift prediction workshop,” _49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition_ , 2011, p. 862\.
  * Rumsey et al. [2023] Rumsey, C. L., Slotnick, J. P., and Woeber, C. D., “Fourth High-Lift Prediction/Third Geometry and Mesh Generation Workshops: Overview and Summary,” _J. Aircr._ , 2023, pp. 1–18.
  * Milano and Koumoutsakos [2002] Milano, M., and Koumoutsakos, P., “Neural network modeling for near wall turbulent flow,” _J Comput Phys_ , Vol. 182, No. 1, 2002, pp. 1–26.
  * Duraisamy et al. [2015] Duraisamy, K., Zhang, Z. J., and Singh, A. P., “New approaches in turbulence and transition modeling using data-driven techniques,” _53rd AIAA Aerospace sciences meeting_ , 2015, p. 1284.
  * Ling et al. [2016] Ling, J., Kurzawski, A., and Templeton, J., “Reynolds averaged turbulence modelling using deep neural networks with embedded invariance,” _J. Fluid Mech._ , Vol. 807, 2016, pp. 155–166.
  * Parish and Duraisamy [2016] Parish, E. J., and Duraisamy, K., “A paradigm for data-driven predictive modeling using field inversion and machine learning,” _J Comput Phys_ , Vol. 305, 2016, pp. 758–774.
  * Singh et al. [2017] Singh, A. P., Medida, S., and Duraisamy, K., “Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils,” _AIAA J._ , Vol. 55, No. 7, 2017, pp. 2215–2227.
  * Han et al. [2022] Han, J., Zhou, X.-H., and Xiao, H., “VCNN-e: A vector-cloud neural network with equivariance for emulating Reynolds stress transport equations,” _arXiv preprint arXiv:2201.01287_ , 2022.
  * Rumsey et al. [2022] Rumsey, C. L., Coleman, G. N., and Wang, L., “In search of data-driven improvements to RANS models applied to separated flows,” _AIAA Scitech 2022 Forum_ , 2022, p. 0937.
  * Spalart [2022] Spalart, P. R., “An Old-Fashioned Framework for Machine Learning in Turbulence Modeling,” , 2022.
  * Holland et al. [2019] Holland, J. R., Baeder, J. D., and Duraisamy, K., “Towards integrated field inversion and machine learning with embedded neural networks for RANS modeling,” _AIAA Scitech forum_ , 2019, p. 1884.
  * Shirian and Mani [2022] Shirian, Y., and Mani, A., “Eddy diffusivity operator in homogeneous isotropic turbulence,” _Phys. Rev. Fluids_ , Vol. 7, No. 5, 2022, p. L052601.
  * Mani and Park [2021] Mani, A., and Park, D., “Macroscopic forcing method: a tool for turbulence modeling and analysis of closures,” _Phys. Rev. Fluids_ , Vol. 6, No. 5, 2021, p. 054607.
  * Menter et al. [2019] Menter, F., Lechner, R., and Matyushenko, A., “Best practice: generalized k-$\omega$ two-equation turbulence model in ANSYS CFD (GEKO),” _ANSYS Germany GmbH_ , 2019.
  * Bin et al. [2022] Bin, Y., Chen, L., Huang, G., and Yang, X. I. A., “Progressive, extrapolative machine learning for near-wall turbulence modeling,” _Phys. Rev. Fluids_ , Vol. 7, No. 8, 2022, p. 084610.
  * Bin et al. [2023a] Bin, Y., Huang, G., and Yang, X. I. A., “Data-Enabled Recalibration of the Spalart–Allmaras Model,” _AIAA J._ , 2023a, pp. 1–12.
  * Vadrot et al. [2023] Vadrot, A., Yang, X. I. A., and Abkar, M., “Survey of machine-learning wall models for large-eddy simulation,” _Phys. Rev. Fluids_ , Vol. 8, No. 6, 2023, p. 064603.
  * Wu et al. [2018] Wu, J.-L., Xiao, H., and Paterson, E., “Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework,” _Phys. Rev. Fluids_ , Vol. 3, No. 7, 2018, p. 074602.
  * Spalart [2015] Spalart, P. R., “Philosophies and fallacies in turbulence modeling,” _Prog. Aerosp. Sci._ , Vol. 74, 2015, pp. 1–15.
  * Wang et al. [2017] Wang, J.-X., Wu, J.-L., and Xiao, H., “Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data,” _Phys. Rev. Fluids_ , Vol. 2, No. 3, 2017, p. 034603.
  * Pope [1975] Pope, S. B., “A more general effective-viscosity hypothesis,” _J. Fluid Mech._ , Vol. 72, No. 2, 1975, pp. 331–340.
  * Parish et al. [2023] Parish, E., Ching, D. S., Miller, N. E., Beresh, S. J., and Barone, M. F., “Turbulence modeling for compressible flows using discrepancy tensor-basis neural networks and extrapolation detection,” _AIAA SCITECH 2023 Forum_ , 2023, p. 2126.
  * Xu et al. [2023] Xu, D., Wang, J., Yu, C., and Chen, S., “Artificial-neural-network-based nonlinear algebraic models for large-eddy simulation of compressible wall-bounded turbulence,” _J. Fluid Mech._ , Vol. 960, 2023, p. A4.
  * Xie et al. [2020] Xie, C., Yuan, Z., and Wang, J., “Artificial neural network-based nonlinear algebraic models for large eddy simulation of turbulence,” _Phys. Fluids_ , Vol. 32, No. 11, 2020, p. 115101.
  * Zhang et al. [2019] Zhang, X., Wu, J., Coutier-Delgosha, O., and Xiao, H., “Recent progress in augmenting turbulence models with physics-informed machine learning,” _J. Hydrodyn._ , Vol. 31, 2019, pp. 1153–1158.
  * Weatheritt and Sandberg [2016] Weatheritt, J., and Sandberg, R., “A novel evolutionary algorithm applied to algebraic modifications of the RANS stress–strain relationship,” _J Comput Phys_ , Vol. 325, 2016, pp. 22–37.
  * Zhao et al. [2020] Zhao, Y., Akolekar, H. D., Weatheritt, J., Michelassi, V., and Sandberg, R. D., “RANS turbulence model development using CFD-driven machine learning,” _J Comput Phys_ , Vol. 411, 2020, p. 109413.
  * Fang et al. [2023] Fang, Y., Zhao, Y., Waschkowski, F., Ooi, A. S., and Sandberg, R. D., “Toward More General Turbulence Models via Multicase Computational-Fluid-Dynamics-Driven Training,” _AIAA J._ , Vol. 61, No. 5, 2023, pp. 2100–2115.
  * Bin et al. [2023b] Bin, Y., Huang, G., and Yang, X. I. A., “A data-enabled re-calibration of the Spalart-Allmaras model for general purposes,” _AIAA J._ , 2023b.
  * Fluent [2011] Fluent, A., “Fluent 14.0 user’s guide,” _Ansys Fluent Inc_ , 2011\.
  * Weller et al. [1998] Weller, H. G., Tabor, G., Jasak, H., and Fureby, C., “A tensorial approach to computational continuum mechanics using object-oriented techniques,” _Comput. phys._ , Vol. 12, No. 6, 1998, pp. 620–631.
  * Daly and Harlow [1970] Daly, B. J., and Harlow, F. H., “Transport equations in turbulence,” _Phys. Fluids_ , Vol. 13, No. 11, 1970, pp. 2634–2649.
  * Guo et al. [2021a] Guo, X., Xia, Z., and Chen, S., “Practical framework for data-driven RANS modeling with data augmentation,” _Acta Mech. Sin._ , Vol. 37, No. 12, 2021a, pp. 1748–1756.
  * Duraisamy [2021] Duraisamy, K., “Perspectives on machine learning-augmented Reynolds-averaged and large eddy simulation models of turbulence,” _Phys. Rev. Fluids_ , Vol. 6, No. 5, 2021, p. 050504.
  * Moser et al. [1999] Moser, R. D., Kim, J., and Mansour, N. N., “Direct numerical simulation of turbulent channel flow up to Reτ= 590,” _Phys. Fluids_ , Vol. 11, No. 4, 1999, pp. 943–945.
  * Lee and Moser [2015] Lee, M., and Moser, R. D., “Direct numerical simulation of turbulent channel flow up to $Re_{\tau}=5200$,” _J. Fluid Mech._ , Vol. 774, 2015, pp. 395–415.
  * Graham et al. [2016] Graham, J., Kanov, K., Yang, X. I. A., Lee, M., Malaya, N., Lalescu, C., Burns, R., Eyink, G., Szalay, A., Moser, R., et al., “A web services accessible database of turbulent channel flow and its use for testing a new integral wall model for LES,” _J. Turbul._ , Vol. 17, No. 2, 2016, pp. 181–215.
  * Huang et al. [2021] Huang, X. L., Jain, N., Abkar, M., Kunz, R. F., and Yang, X. I. A., “Determining a priori a RANS model’s applicable range via global epistemic uncertainty quantification,” _Comput. fluids_ , Vol. 230, 2021, p. 105113.
  * Oliver and Moser [2011] Oliver, T. A., and Moser, R. D., “Bayesian uncertainty quantification applied to RANS turbulence models,” _J. Phys. Conf. Ser._ , Vol. 318, No. 4, 2011, p. 042032.
  * Jain et al. [2022] Jain, N., Pham, H. T., Huang, X., Sarkar, S., Yang, X. I. A., and Kunz, R., “Second Moment Closure Modeling and Direct Numerical Simulation of Stratified Shear Layers,” _J. Fluid Eng._ , Vol. 144, No. 4, 2022, p. 041102\.
  * Fang et al. [2020] Fang, R., Sondak, D., Protopapas, P., and Succi, S., “Neural network models for the anisotropic Reynolds stress tensor in turbulent channel flow,” _J. Turbul._ , Vol. 21, No. 9-10, 2020, pp. 525–543.
  * Guo et al. [2021b] Guo, X., Xia, Z., and Chen, S., “Computing mean fields with known Reynolds stresses at steady state,” _Theor. App. Mech._ , Vol. 11, No. 3, 2021b, p. 100244.
  * Wu et al. [2019] Wu, J., Xiao, H., Sun, R., and Wang, Q., “Reynolds-averaged Navier–Stokes equations with explicit data-driven Reynolds stress closure can be ill-conditioned,” _J. Fluid Mech._ , Vol. 869, 2019, pp. 553–586.
  * Rumsey et al. [[retrieved 12 Aug. 2015]] Rumsey, C., Smith, B., and Huang, G., “Turbulence Modeling Resource Website [Online Database],” https://turbmodels.larc.nasa.gov/, [retrieved 12 Aug. 2015].
  * Delville et al. [1989] Delville, J., Bellin, S., Garem, J., and Bonnet, J., “Analysis of structures in a turbulent, plane mixing layer by use of a pseudo flow visualization method based on hot-wire anemometry,” _Advances in Turbulence 2: Proceedings of the Second European Turbulence Conference, Berlin, August 30–September 2, 1988_ , Springer, 1989, pp. 251–256.
  * Srivastava [2022] Srivastava, V., “Generalizable Data-driven Model Augmentations Using Learning and Inference assisted by Feature-space Engineering,” Ph.D. thesis, University of Michigan, 2022.
  * Rincón et al. [2023] Rincón, M. J., Amarloo, A., Reclari, M., Yang, X. I. A., and Abkar, M., “Progressive augmentation of Reynolds stress tensor models for secondary flow prediction by computational fluid dynamics driven surrogate optimisation,” _Int. J. Heat Fluid Flow_ , Vol. 104, 2023, p. 109242.