1.1. Injecting Reynolds stress tensor decompositions used in data-driven closure modelling frameworks

The steady-state, RANS momentum equations for an incompressible, Newtonian fluid are given by (Reynolds Reference Reynolds1895)

(1.1) \begin{equation} U_i \frac {\partial U_{\!j}}{\partial x_i} = -\frac {1}{\rho } \frac {\partial P}{\partial x_j} + \nu \frac {\partial ^2 U_{\!j}}{\partial x_i \partial x_i} -\frac {\partial \tau _{\textit{ij}}}{\partial x_i} , \end{equation}

where $U_{\!j}$ is the mean velocity, $P$ is the mean pressure, $\rho$ is the fluid density, $\nu$ is the kinematic viscosity and $\tau _{\textit{ij}}$ is the Reynolds stress tensor.

The continuity equation also applies to this flow

(1.2) \begin{equation} \frac {\partial U_i}{\partial x_i} = 0 . \end{equation}

Together, (1.1) and (1.2) are unclosed. In the most general case, there are four equations (continuity + 3 momentum), and 10 unknowns: $P$ , 3 components of $U_i$ and 6 components of $\tau _{\textit{ij}}$ . The goal of turbulence closure modelling is to express $\tau _{\textit{ij}}$ in terms of $P$ and $U_i$ .

The Reynolds stress tensor can be decomposed into isotropic (hydrostatic) and anisotropic (deviatoric) components

(1.3) \begin{equation} \tau _{\textit{ij}} = \frac {2}{3}k\delta _{\textit{ij}} + a_{\textit{ij}} , \end{equation}

where $k = ({1}/{2})\tau _{kk}$ is half the trace of the Reynolds stress tensor, and $a_{\textit{ij}}$ is the anisotropy tensor. The isotropic component ( $({2}/{3})k\delta _{\textit{ij}}$ ) in (1.3) can be absorbed into the pressure gradient term in (1.1) to form a modified pressure

(1.4) \begin{equation} U_i \frac {\partial U_{\!j}}{\partial x_i} = -\frac {1}{\rho } \frac {\partial }{\partial x_j} \left (P+\frac {2}{3}\rho k\right ) + \nu \frac {\partial ^2 U_{\!j}}{\partial x_i \partial x_i} -\frac {\partial a_{\textit{ij}}}{\partial x_i} . \end{equation}

Equation (1.4) is the form of the RANS momentum equation commonly used in RANS turbulence closure modelling. Several data-driven closure modelling frameworks are based on modelling the Reynolds stress tensor itself, or it’s divergence (as in Brener et al. (Reference Brener, Cruz, Macedo and Thompson2022)), which imply the use of (1.1). In our work, we model the Reynolds stress anisotropy tensor $a_{\textit{ij}}$ , implying the use of (1.4). In either case, the term ‘injection’ refers to using a model prediction as the closure term in the momentum equation and numerically solving the momentum equation with this predicted closure term.

In an eddy viscosity hypothesis, the anisotropy tensor $a_{\textit{ij}}$ is postulated to be a function of the mean strain rate $S_{\textit{ij}}$ and rotation rate $R_{\textit{ij}}$ tensors

(1.5) \begin{align} a_{\textit{ij}} &= a_{\textit{ij}}(S_{\textit{ij}}, R_{\textit{ij}}) , \end{align}

(1.6) \begin{align} S_{\textit{ij}} &= \frac {1}{2}\left (\frac {\partial U_i}{\partial x_j} + \frac {\partial U_{\!j}}{\partial x_i}\right )\! , \end{align}

(1.7) \begin{align} R_{\textit{ij}} &= \frac {1}{2}\left (\frac {\partial U_i}{\partial x_j} - \frac {\partial U_{\!j}}{\partial x_i}\right)\!. \end{align}

The eddy viscosity hypothesis implies the following important constraints:

1. (i) The anisotropy tensor can be predicted locally and instantaneously. In (1.5), there is no temporal dependence.

2. (ii) The anisotropic turbulent stresses are caused entirely by mean velocity gradients.

There are many examples of turbulent flows where these assumptions are violated (Pope Reference Pope2000). For example, history effects occur in boundary layers (Bobke et al. Reference Bobke, Vinuesa, Örlü and Schlatter2017). Nevertheless these approximations remain widely used in eddy viscosity modelling. The most common eddy viscosity hypothesis is the linear eddy viscosity hypothesis

(1.8) \begin{equation} a_{\textit{ij}} = - 2 \nu _t S_{\textit{ij}} , \end{equation}

where $\nu _t$ is the eddy viscosity, a scalar. This linear hypothesis draws direct analogy from the stress–strain rate relation for a Newtonian fluid. Along with the important constraints implied by invoking an eddy viscosity hypothesis, the linear eddy viscosity hypothesis implies the following:

1. (i) The anisotropy tensor is aligned with the mean strain rate tensor.

2. (ii) The mapping between mean strain rate and anisotropic stress is isotropic, in that it can be represented using a single scalar ( $\nu _t$ ).

More general nonlinear eddy viscosity hypotheses have been used in several models. The primary advantage of these models is that they permit a misalignment of the principal axes of $a_{\textit{ij}}$ and $S_{\textit{ij}}$ , which occurs for even simple flows. Applying Cayley–Hamilton theorem (Hamilton Reference Hamilton1853; Cayley Reference Cayley1858) to (1.5), Pope (Reference Pope1975) derived the most general expression for a nonlinear eddy viscosity model

(1.9) \begin{equation} a_{\textit{ij}} = 2k \sum _{n=1}^{10} g_n \hat {T}^{(n)}_{\textit{ij}} , \end{equation}

where $n=1,2,\ldots 10$ indexes the scalar coefficients $g_n$ , and the following basis tensors:

(1.10) \begin{align} & \hat {T}^{(1)}_{\textit{ij}} =\hat {S}_{\textit{ij}}, && \!\! \hat {T}^{(6)}_{\textit{ij}}= \hat {R}_{ik} \hat {R}_{kl} \hat {S}_{lj}+ \hat {S}_{ik} \hat {R}_{kl} \hat {R}_{lj} - \tfrac {2}{3} \hat {S}_{kl}\hat {R}_{lm}\hat {R}_{mk} \delta _{\textit{ij}} ,\nonumber \\ & \hat {T}^{(2)}_{\textit{ij}} =\hat {S}_{ik} \hat {R}_{kj} - \hat {R}_{ik}\hat {S}_{kj}, && \!\! \hat {T}^{(7)}_{\textit{ij}} = \hat {R}_{ik} \hat {S}_{kl} \hat {R}_{lm}\hat {R}_{mj} - \hat {R}_{ik}\hat {R}_{kl} \hat {S}_{lm} \hat {R}_{mj},\nonumber \\ & \hat {T}^{(3)}_{\textit{ij}} =\hat {S}_{ik}\hat {S}_{kj} - \tfrac {1}{3}\hat {S}_{kl}\hat {S}_{lk} \delta _{\textit{ij}}, && \!\! \hat {T}^{(8)}_{\textit{ij}} = \hat {S}_{ik}\hat {R}_{kl}\hat {S}_{lm}\hat {S}_{mj} - \hat {S}_{ik}\hat {S}_{kl} \hat {R}_{lm} \hat {S}_{mj}, \nonumber \\ & \hat {T}^{(4)}_{\textit{ij}} \!=\! \hat {R}_{ik}\hat {R}_{kj} \!-\! \tfrac {1}{3}\hat {R}_{kl}\hat {R}_{lk}\delta _{\textit{ij}}, && \!\! \hat {T}^{(9)}_{\textit{ij}} \!=\! \hat {R}_{ik} \hat {R}_{kl}\hat {S}_{lm}\hat {S}_{mj} \!+\! \hat {S}_{ik} \hat {S}_{kl}\hat {R}_{lm}\hat {R}_{mj} \!-\! \tfrac {2}{3} \hat {S}_{kl}\hat {S}_{lm} \hat {R}_{mo}\hat {R}_{ok} \delta _{\textit{ij}}, \nonumber \\ & \hat {T}^{(5)}_{\textit{ij}} = \hat {R}_{ik} \hat {S}_{kl} \hat {S}_{lj} - \hat {S}_{ik} \hat {S}_{kl}\hat {R}_{lj}, && \!\! \hat {T}^{(10)}_{\textit{ij}} = \hat {R}_{ik}\hat {S}_{kl}\hat {S}_{lm} \hat {R}_{mo} \hat {R}_{oj} - \hat {R}_{ik}\hat {R}_{kl}\hat {S}_{lm} \hat {S}_{mo} \hat {R}_{oj}\ . \end{align}

The tensors $\hat {S}_{\textit{ij}}$ and $\hat {R}_{\textit{ij}}$ above are the non-dimensionalised strain rate and rotation rate tensors. In Pope’s original work, these tensors are given by

(1.11) \begin{align} \hat {S}_{\textit{ij}} &= \frac {k}{\varepsilon } S_{\textit{ij}} , \end{align}

(1.12) \begin{align} \hat {R}_{\textit{ij}} &= \frac {k}{\varepsilon } R_{\textit{ij}}. \end{align}

However, other normalisation constants are possible, and in several nonlinear eddy viscosity models, different basis tensors use different normalisation constants. To our knowledge, varying the basis tensor normalisation has not been explored in data-driven turbulence closure modelling.

With Pope’s general expression for the anisotropy tensor, (1.4) becomes

(1.13) \begin{equation} U_i \frac {\partial U_{\!j}}{\partial x_i} = -\frac {1}{\rho } \frac {\partial }{\partial x_j} \left (P+\frac {2}{3}\rho k\right ) + \nu \frac {\partial ^2 U_{\!j}}{\partial x_i \partial x_i} -\frac {\partial }{\partial x_i}\left (2k \sum _{n=1}^{10} g_n \hat {T}^{(n)}_{\textit{ij}}\right )\! . \end{equation}

This is the form of the momentum equation used in several studies which aim to augment the closure relationship via machine learning. For example, in Ling et al.’s TBNN investigation (Ling et al. Reference Ling, Kurzawski and Templeton2016), (1.13) was used. Kaandorp (Reference Kaandorp2018) and Kaandorp & Dwight (Reference Kaandorp and Dwight2020) also used this form of the momentum equation. The TBNN-based approaches are advantageous for several reasons. The TBNN architecture is motivated by a tensor algebra-based argument that if the anisotropy tensor is to be expanded in terms of the mean strain and rotation rate tensors (a dominant approach in eddy viscosity modelling), the TBNN is the most general form of this expansion (Pope Reference Pope1975). Additionally, other techniques to incorporate equivariance are cumbersome, requiring learning and predicting in an invariant eigenframe of some tensorial quantity (Wu et al. Reference Wu, Xiao and Paterson2018; Brener et al. Reference Brener, Cruz, Macedo and Thompson2022). The TBNN-based architectures do not require computing eigenvalues of the strain rate tensor at every evaluation data point. Lastly, this closure term is highly expressive, in that ten different combinations of $S_{\textit{ij}}$ and $R_{\textit{ij}}$ can be used to represent the anisotropy tensor. However, a major disadvantage with numerically solving (1.13) is that the closure term is entirely explicit, greatly reducing numerical stability and conditioning (Brener et al. Reference Brener, Cruz, Thompson and Anjos2021). For this reason, Kaandorp & Dwight (Reference Kaandorp and Dwight2020) needed to implement a blending function, which blends the fully explicit closure term in (1.13) with the more stable implicit closure term treatment made possible with a linear eddy viscosity hypothesis. Here, implicit treatment refers to the way the velocity discretisation matrix is constructed in the numerical solver. When a term is treated implicitly, it contributes to the $U_i$ coefficient matrix. In contrast, explicit treatment refers to treating a term as a fixed source term in the discretised equation. Assuming the closure term takes the form $a_{\textit{ij}} = -2 \nu _t S_{\textit{ij}}$ , (1.4) can be written as

(1.14) \begin{equation} U_i \frac {\partial U_{\!j}}{\partial x_i} = -\frac {1}{\rho } \frac {\partial }{\partial x_j} \left (P+\frac {2}{3}\rho k\right ) + \frac {\partial }{\partial x_i} \left [(\nu + \nu _t) \frac {\partial U_{\!j}}{\partial x_i}\right ]\! . \end{equation}

Equation (1.14) has the major advantage of increasing diagonal dominance of the $U_{\!j}$ coefficient matrix obtained from discretisation of this equation, via the eddy viscosity. However, this closure framework only permits the inaccurate linear eddy viscosity closure approximation.

In the present work, we propose the following hybrid treatment of the closure term:

(1.15) \begin{equation} U_i \frac {\partial U_{\!j}}{\partial x_i} = -\frac {1}{\rho } \frac {\partial }{\partial x_j} \left (P+\frac {2}{3}\rho k\right ) + \frac {\partial }{\partial x_i} \left [(\nu + \nu _t) \frac {\partial U_{\!j}}{\partial x_i}\right ] -\frac {\partial }{\partial x_i}\left (2k \sum _{n=2}^{10} g_n \hat {T}^{(n)}_{\textit{ij}}\right )\! . \end{equation}

In (1.15), the $n=1$ (linear) term has been separated and receives implicit treatment, while the remaining $n = 2,3,\ldots 10$ terms grant an opportunity for a machine learning model to provide rich representation of the nonlinear part of $a_{\textit{ij}}$ . As we will discuss, the separation of the linear term requires special treatment during training and closure term injection.

1.2. Conditioning analysis

Various decompositions of the Reynolds stress tensor were investigated by Brener et al.’s conditioning analysis (Brener et al. Reference Brener, Cruz, Thompson and Anjos2021). Conditioning analysis for data-driven turbulence closure frameworks is important, since ill-conditioned momentum equations have the potential to amplify errors in the predicted closure term. Brener et al. (Reference Brener, Cruz, Thompson and Anjos2021) concluded that an optimal eddy viscosity approach is necessary to achieve a well-conditioned solution, since it incorporates information about the DNS mean velocity field. In the present work, we demonstrate that the following closure decomposition:

(1.16) \begin{equation} a_{\textit{ij}} = -2\nu ^{{R}}_{t} S^{\theta}_{\textit{ij}} + a^{\perp}_{\textit{ij}} \end{equation}

also achieves a well-conditioned solution. We follow the nomenclature of Duraisamy et al. (Reference Duraisamy, Iaccarino and Xiao2019) in that the $\theta$ superscript indicates a quantity that comes from a high-fidelity source such as DNS. The ${{R}}$ superscript indicates a quantity taken from the corresponding baseline RANS simulation.

The decomposition in (1.16) permits an augmented turbulence closure framework that treats the machine learning correction only in an explicit term in the momentum equation. Separating the machine learning model prediction has several advantages. It allows the model correction to be easily ‘turned off’ in unstable situations. It also is more interpretable – rather than correcting both the eddy viscosity and also injecting an explicit correction term, the correction is contained entirely within an explicit term in the momentum equation. Lastly, it avoids the necessity for an optimal eddy viscosity to be computed from high-fidelity data. Although there are methods to increase the practicality of computing the optimal eddy viscosity (McConkey et al. Reference McConkey, Yee and Lien2022a ), this quantity is often unstable and difficult to predict via a machine learning model.

Figure 1. Conditioning test comparing the errors in $U_1$ after injecting: (a) the DNS anisotropy tensor fully explicitly, (b) the optimal eddy viscosity (implicitly) and the remaining nonlinear part of the anisotropy tensor calculated using $S^\theta _{\textit{ij}}$ (explicitly), (c) the eddy viscosity estimated by the $k$ - $\omega$ shear stress transport (SST) model (implicitly) and the remaining nonlinear part of the anisotropy tensor calculated using $S^\theta _{\textit{ij}}$ (explicitly) and (d) the eddy viscosity estimated by the $k$ - $\omega$ SST model (implicitly), and the remaining nonlinear part of the anisotropy tensor calculated using $S^{{R}}_{\textit{ij}}$ (explicitly), $k$ is the turbulent kinetic energy (TKE), $\omega$ is the TKE specific dissipation rate and RSME is the Root mean squared error.

Figure 2. Conditioning test comparing the errors in $U_2$ after injecting: (a) the DNS anisotropy tensor fully explicitly, (b) the optimal eddy viscosity (implicitly) and the remaining nonlinear part of the anisotropy tensor calculated using $S^\theta _{\textit{ij}}$ (explicitly), (c) the eddy viscosity estimated by the $k$ - $\omega$ SST model (implicitly) and the remaining nonlinear part of the anisotropy tensor calculated using $S^\theta _{\textit{ij}}$ (explicitly) and (d) the eddy viscosity estimated by the $k$ - $\omega$ SST model (implicitly), and the remaining nonlinear part of the anisotropy tensor calculated using $S^{{R}}_{\textit{ij}}$ (explicitly).

Figures 1 and 2 demonstrate a conditioning test similar to the tests conducted by Brener et al. (Reference Brener, Cruz, Thompson and Anjos2021). Several different decompositions of the DNS Reynolds stress tensor are injected into a RANS simulation, to identify which decompositions permit a well-conditioned solution. Comparing panels (b) and (c) in both figures 1 and 2, we can see that the proposed decomposition of the Reynolds stress tensor (1.16) achieves an equally well-conditioned solution as the optimal eddy viscosity framework. To our knowledge, this is the first result in the literature showing that a well-conditioned solution can be achieved without the use of an optimal eddy viscosity to incorporate information about the DNS velocity field (Brener et al. Reference Brener, Cruz, Thompson and Anjos2021). We further confirm the findings of Brener et al. (Reference Brener, Cruz, Thompson and Anjos2021) with respect to the requirement that the closure decomposition includes information about the DNS mean velocity field in order to achieve a well-conditioned solution. We confirm the notion from Wu et al. (Reference Wu, Xiao, Sun and Wang2019a ) that implicit treatment helps address the ill-conditioning issue, but using the DNS strain rate tensor $S^\theta _{\textit{ij}}$ to calculate the explicitly injected $a^{\perp}_{\textit{ij}}$ .

1.3. Realisability-informed training

The Reynolds stress tensor is symmetric positive semidefinite. A set of constraints on the non-dimensional anisotropy tensor $b_{\textit{ij}}$ arise from this property, as determined by Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007). These constraints are

(1.17) \begin{align} -\frac {1}{3} \leqslant b_{\textit{ij}} &\leqslant \frac {2}{3}, \,\, i= j, \end{align}

(1.18) \begin{align} -\frac {1}{2} \leqslant b_{\textit{ij}} &\leqslant \frac {1}{2}, \,\, i\neq j, \end{align}

(1.19) \begin{align} \lambda _1 &\geqslant \frac {3|\lambda _2| - \lambda _2}{2}, \end{align}

(1.20) \begin{align} \lambda _1 &\leqslant \frac {1}{3}- \lambda _2, \\[6pt] \nonumber \end{align}

where the non-dimensional anisotropy tensor $b_{\textit{ij}}$ is calculated by

(1.21) \begin{equation} b_{\textit{ij}} = \frac {a_{\textit{ij}}}{2k} , \end{equation}

and the eigenvalues of $b_{\textit{ij}}$ are given by $\lambda _1 \geqslant \lambda _2 \geqslant \lambda _3$ .

A given Reynolds stress tensor is physically realisable if it satisfies these constraints. While the physical realisability of the closure term may seem an important constraint for turbulence models, many commonly used turbulence models such as the $k$ - $\varepsilon$ model, $\varepsilon$ is the TKE dissipation rate (Launder & Spalding Reference Launder and Spalding1974), $k$ - $\omega$ (Wilcox Reference Wilcox1988) and $k$ - $\omega$ SST model (Menter Reference Menter1994; Menter, Kuntz & Langtry Reference Menter, Kuntz and Langtry2003) do not guarantee physical realisability.

Based on the widespread acceptance and popularity of non-realisable turbulence models, it is fair to say that realisability is not a hard constraint on a new turbulence model. Nevertheless, the true Reynolds stress tensor is realisable, and if a machine learning model is able to learn to predict realisable closure terms, it may be more physically accurate. Unfortunately, Pope’s tensor basis expansion for the anisotropy tensor (and machine learning architectures based on this expansion) do not provide a means to achieve realisability. To enforce realisability, a variety of ad hoc strategies have been used, including in the original TBNN paper by Ling et al. (Reference Ling, Kurzawski and Templeton2016). Most of these strategies involve postprocessing predictions by the TBNN, such as shrinking the predicted anisotropy tensor in certain directions until it is physically realisable (Jiang et al. Reference Jiang, Vinuesa, Chen, Mi, Laima and Li2021).

In the present work, we propose including a penalty for violating realisability constraints in the loss function. In a similar spirit of physics-informed neural networks (PINNs) (Raissi, Perdikaris & Karniadakis Reference Raissi, Perdikaris and Karniadakis2019), we term the use of this loss function ‘realisability-informed training’. Whereas PINNs encourage the model to learn physics by penalising violations of conservation laws, realisability-informed TBNNs learn to predict physically realisable anisotropy tensors.

The realisability penalty $\mathcal{R}(b_{\textit{ij}})$ is given as follows:

(1.22) \begin{align} \mathcal{R}(\tilde b_{\textit{ij}}) &= \frac {1}{6}\left \lbrace \sum _{i=j} \left ( {\textrm{max}}\left [\tilde b_{\textit{ij}}-\frac {2}{3}, - \left (\tilde b_{\textit{ij}}+\frac {1}{3}\right )\!, 0\right ] \right )^2 \right . \nonumber \\ &\quad +\, \left . \sum _{i\neq j} \left ({\textrm{max}}\left [\tilde b_{\textit{ij}}-\frac {1}{2}, - \left (\tilde b_{\textit{ij}}+\frac {1}{2}\right )\!, 0\right ]\right )^2 \right \rbrace \nonumber \\ &\quad +\, \frac {1}{2}\left \lbrace \left ({\textrm{max}}\left [\frac {3|\lambda _2| - |\lambda _2|}{2} - \lambda _1, 0 \right ]\right )^2 \right . \nonumber \\ &\quad +\, \left . \left ({\textrm{max}}\left [\lambda _1 - \left ( \frac {1}{3}-\lambda _2\right )\!,0 \right ] \right )^2 \right \rbrace \!, \end{align}

where the $\sim$ above a symbol denotes a model prediction for the quantity associated with the symbol.

Equation (1.22) can be thought of as the mean squared violation in the components of $\tilde b_{\textit{ij}}$ , plus the mean squared violation in the eigenvalues of $\tilde b_{\textit{ij}}$ . To help visualise this penalty function, figure 3 shows the penalties incurred by violating various realisability constraints on the anisotropy tensor.

Figure 3. Realisability penalty for the (a) diagonal components of $b_{\textit{ij}}$ , and (b) off-diagonal components of $b_{\textit{ij}}$ .

It should be noted that in the same way that a PINN’s prediction cannot be guaranteed to satisfy a conservation law, a realisability-informed TBNN cannot be guaranteed to predict a physically realisable anisotropy tensor. The goal is that the incorporation of a physics-based loss at training time will encode into the model a tendency to predict physically realisable anisotropy tensors. At training time, when the model predicts a $b_{\textit{ij}}$ component outside the realisability zone, then it is penalised in two ways: the error in this prediction will be non-zero (since all label data are realisable), and the realisability penalty will be non-zero. Therefore, the realisability penalty serves as an additional driving factor towards the realisable, true $b_{\textit{ij}}$ value. The merits of this approach are demonstrated in § 2.2.

Another important metric in the loss function is the error-based penalty $\mathcal{E}(\tilde b_{\textit{ij}})$ . When the model predicts a certain anisotropy tensor, it is evaluated against a known high-fidelity anisotropy tensor available in the training dataset. Typically, mean-squared error loss functions are used to train machine learning-augmented closure models. However, we propose the following modifications to the loss function of a TBNN.

1. (i) Since $b_{\textit{ij}}$ is a symmetric tensor, we propose to sum the squared errors as follows:

   (1.23) \begin{equation} \mathcal{E}(\tilde b_{\textit{ij}}) = \frac {1}{6}\left \lbrace \sum _{\substack {\textit{ij} \in \lbrace 11, 12, 13,\\ 22,23,33\rbrace }} \big(\tilde b_{\textit{ij}} - b^\theta _{\textit{ij}}\big)^2 \right \rbrace \!. \end{equation}
   Calculating the squared error in this way avoids double penalising the off-diagonal components, a situation which arises when summing over all components of $b_{\textit{ij}}$ .

2. (ii) Although the TBNN model predicts $b_{\textit{ij}}$ , we propose to use a loss function based on the error in $a_{\textit{ij}}$ . Near the wall, $\hat {T}^{(n)}_{\textit{ij}}\rightarrow 0$ , but $b_{\textit{ij}} \nrightarrow 0$ . The vanishing of $\hat {T}^{(n)}_{\textit{ij}}$ with non-vanishing $b_{\textit{ij}}$ causes instabilities during training, leading to $g_n \rightarrow \infty$ here. Since $a_{\textit{ij}}$ is the tensor injected into the momentum equation, its accurate prediction should be the focus of the training process. To dimensionalise $b_{\textit{ij}}$ , the turbulent kinetic energy $k$ must be used. This leads to the following error-based loss:

   (1.24) \begin{equation} \mathcal{E}(\tilde a_{\textit{ij}}) = \frac {1}{6}\left \lbrace \sum _{\substack { \textit{ij} \in \lbrace 11, 12, 13,\\ 22,23,33\rbrace }} \big( 2k^\theta \tilde b_{\textit{ij}} - 2k^\theta b^\theta _{\textit{ij}}\big)^2 \right \rbrace = \frac {2(k^\theta )^2}{3}\left \lbrace \sum _{\substack { \textit{ij} \in \lbrace 11, 12, 13,\\ 22,23,33\rbrace }} \big(\tilde b_{\textit{ij}} - b^\theta _{\textit{ij}}\big)^2 \right \rbrace\! , \end{equation}
   or stated more simply
   (1.25) \begin{equation} \mathcal{E}(\tilde a_{\textit{ij}})=(2k^\theta )^2 \mathcal{E}(\tilde b_{\textit{ij}}) . \end{equation}
   The reason for the use of $k^\theta$ in (1.24) and (1.25) is discussed in § 1.4.

The final loss function includes an error-based metric and a realisability-violation penalty. This loss function $\mathcal{L}$ is given by

(1.26) \begin{equation} \mathcal{L} = \frac {1}{N}\sum _{p} \frac {\big(2k^\theta _p\big)^2}{\displaystyle \sum _{k=1}^s Z_k^2 {\mathcal I}_{C_k}(p)}\big( \mathcal{E}\big(\tilde b_{ij,p}\big)+\alpha \mathcal R\big(\tilde b_{ij,p}\big)\big), \end{equation}

where $\alpha$ is a factor used to control the relative importance of the realisability penalty. Here, $N$ is the total number of points in the dataset. The points are indexed by $p = 1,2,\ldots ,N$ . The cases in the dataset are indexed $k=1,2,\ldots ,s$ . ${\mathcal I}_{C_k}(p)$ is an indicator function

(1.27) \begin{equation} {\mathcal I}_{C_k}(p) = \begin{cases} 1,\quad {\textrm{if}} p \in C_k,\\ 0 \quad {\textrm{otherwise}}. \end{cases} \end{equation}

For a given point $p$ , ${\mathcal I}_{C_k}(p)$ selects all points which come from the same case as $p$ . The denominator for all points from the same case $C_k$ is the same – normalisation is applied on a case-by-case basis. A given point is normalised by the mean-squared Frobenius norm of the anisotropy tensor over all the points from the case it comes from. The mean Frobenius norm over a case is given by

(1.28) \begin{equation} Z_k = \frac {1}{|C_k|}\sum _{p \in C_k} \| a_{\textit{ij}}(p)\|_F\ ,\qquad k=1,2,\ldots ,s\ , \end{equation}

where $|C_k|$ is the cardinality of the case $C_k$ (viz., the number of points in the $k$ th case $C_k$ ). Here, $Z_k$ is simply the average of the Frobenius norm $\| \boldsymbol{\cdot } \|_F$ of the anisotropy tensor $a_{\textit{ij}}$ for all points $p$ (viz., $a_{\textit{ij}}(p)$ ) in case $C_k$ . The objective of this denominator in (1.26) is to promote a more balanced regression problem, since data points from various cases or flow types may have $||a_{\textit{ij}}||$ that differ by orders of magnitude. Normalisation on a case-by-case basis is made in an effort to normalise all $a_{\textit{ij}}$ error magnitudes to a similar scale.

In this study, we use $\alpha =10^2$ to encode a high preference for realisable results in § 2. Lower values of $\alpha$ will reduce the penalty applied to realisability violations, which may be necessary for flows in which the anisotropy tensor is difficult to predict via a TBNN. As discussed, the multiplicative term $(2k^\theta )^2$ is used to formulate the loss function in terms of predicting $a_{\textit{ij}}$ rather than $b_{\textit{ij}}$ . However, $\mathcal{R}(\tilde b_{\textit{ij}})$ is also multiplied by $(2k^\theta )^2$ to ensure that the realisability penalty and mean-squared error in $a_{\textit{ij}}$ are of similar scales.

1.4. Neural network architecture

Motivated by improving the training and injection stability, as well as generalisability, we propose several modifications to the original TBNN (Ling et al. Reference Ling, Kurzawski and Templeton2016).

The original TBNN is shown in figure 4. At training time, this network predicts the non-dimensional anisotropy tensor $b_{\textit{ij}}$ . All basis tensors used in this prediction at training time come from RANS, and during training the prediction $\tilde b_{\textit{ij}}$ is evaluated against a known value of $b^\theta _{\textit{ij}}$ from a high-fidelity simulation. At injection time, the network is used in this same configuration to predict $\tilde b_{\textit{ij}}$ . Here, $\tilde b_{\textit{ij}}$ is injected into a coupled system of equations consisting of the continuity/momentum equations (explicit injection), as well as the turbulence transport equations. This system of equations is iterated around a fixed $\tilde b_{\textit{ij}}$ , to obtain an updated estimate for the turbulent kinetic energy $k$ , and therefore an updated estimate for $\tilde a_{\textit{ij}} = 2 k \tilde b_{\textit{ij}}$ .

Figure 4. Model architecture and training configuration for the original TBNN proposed by Ling et al. (Reference Ling, Kurzawski and Templeton2016).

The modified TBNN is shown in figure 5. This TBNN relies on the same tensor basis expansion as the original TBNN. However, the linear term has been modified in this expansion. Whereas the original TBNN uses

(1.29) \begin{equation} \hat {T}^{(1)}_{\textit{ij}} = \frac {k^{{R}}}{\varepsilon ^{{R}}}S^{{R}}_{\textit{ij}} , \end{equation}

our modified TBNN uses

(1.30) \begin{equation} \hat {T}^{(1)}_{\textit{ij}} = \frac {\nu _t^{{R}}}{k^\theta }S^\theta _{\textit{ij}} . \end{equation}

Figure 5. Model architecture and training configuration proposed in the present investigation. Note that during training, $S^\theta _{\textit{ij}}$ is used in $\hat {T}^{(1)}_{\textit{ij}}$ , whereas all other quantities come from the original RANS simulation.

Further, while the original TBNN calculates the linear (denoted by superscript ${\textrm{L}}$ ) component of $b_{\textit{ij}}$ as

(1.31) \begin{equation} b^{\textit{L}}_{\textit{ij}} = g_1 \hat {T}^{(1)}_{\textit{ij}} = g_1\frac {k}{\varepsilon }S_{\textit{ij}} , \end{equation}

our modified TBNN uses

(1.32) \begin{equation} b^{\textit{L}}_{\textit{ij}} = - \hat {T}^{(1)}_{\textit{ij}} = - \frac {\nu _t}{k}S_{\textit{ij}} , \end{equation}

where the superscript ${{R}}$ denotes a quantity that comes from the original RANS simulation.

These changes are motivated by the following.

1. (i) At injection time, we use implicit treatment of the linear term $\hat {T}^{(1)}_{\textit{ij}}$ to formulate (1.15) in a stable manner. After injection, $S_{\textit{ij}}$ will continue to evolve. In a similar spirit as optimal eddy viscosity frameworks, we therefore use $S^\theta _{\textit{ij}}$ to compute the linear component at training time. In optimal eddy viscosity based frameworks, using $S_{\textit{ij}} = S^\theta _{\textit{ij}}$ at training time helps drive $U_{\!j}\rightarrow U^\theta _{\!j}$ at injection time (the cause of this behaviour is currently unknown). In Ling et al. (Reference Ling, Kurzawski and Templeton2016), all basis tensors (and therefore $b_{\textit{ij}}$ ) remain fixed after injection. We also fixed $\nu _t$ = $\nu ^{{R}}_t$ at injection time (viz., no further evolution of the eddy viscosity is permitted).

2. (ii) At training time, we dimensionalise $\tilde b_{\textit{ij}}$ using $k^\theta$ : $\tilde a_{\textit{ij}} = 2k^\theta \tilde b_{\textit{ij}}$ . At evaluation time, we do not have $k^\theta$ . However, the need for $k^\theta$ is avoided, since we only use the nonlinear part of $\tilde b_{\textit{ij}}$ at test time. The reason we only need $\tilde b^{\perp}_{\textit{ij}}$ at test time is that the training process has been designed to use the linear part of $b_{\textit{ij}}$ estimated by the RANS turbulence model, and augment this by the TBNN’s equivariant prediction for $\tilde b^{\perp}_{\textit{ij}}$ .

3. (iii) Using $\nu _t/k$ to normalise the basis tensors and fixing $g_1=-1$ in (1.32) results in the RANS prediction for $b^{\textit{L}}_{\textit{ij}}$ being implicitly used in the TBNN. Therefore, the TBNN learns to correct $b_{\textit{ij}}$ using $b^{\perp}_{\textit{ij}}$ in a way that allows a realisability-informed training process, and fully implicit treatment of the linear term at injection time.

Lastly, the use of $k^\theta$ to dimensionalise $\tilde b_{\textit{ij}}$ is also enabled by our use of a separate neural network to correct $k^{{R}}$ at injection time. This neural network is called the $k$ -correcting neural network (KCNN), and is a simple fully connected feed-forward neural network that predicts a single output scalar $\varDelta$ (Figure 6)

(1.33) \begin{equation} \varDelta = \log \left (\frac {k^\theta }{k^{{R}}}\right )\! , \end{equation}

such that at injection time, an updated estimate for $k$ can be obtained $\tilde k = {e}^\varDelta k^{{R}}$ , without the need to re-couple the turbulence transport equations. The KCNN shares the same input features as the TBNN.

Figure 6. Test-time injection configuration proposed in the present investigation. The hatching over $g_1$ and $\hat {T}^{(1)}_{\textit{ij}}$ indicates that they are not used at injection time.

Together, the KCNN and TBNN predict the anisotropy tensor in the following manner:

(1.34) \begin{equation} \tilde a_{\textit{ij}} = 2 \big(e^\varDelta k^{{R}} \big) \sum _{n=1}^{10}g_{n} \hat {T}^{(n)}_{\textit{ij}} . \end{equation}

With the linear component of $\tilde a_{\textit{ij}}$ being treated implicitly during injection, the entire closure framework is summarised as explicit injection of the following term into the momentum equation:

(1.35) \begin{equation} \tilde a^{\perp}_{\textit{ij}} = 2 \big(e^\varDelta k^{{R}}\big) \sum _{n=2}^{10} g_{n} \hat {T}^{(n){{R}}}_{\textit{ij}} . \end{equation}

1.4.1. Tensor basis Kolmogorov–Arnold network (TBKAN)

The tensor basis KAN shown in figure 7 in the training configuration replaces the multi-layer perceptron in the modified TBNN with a KAN introduced by Liu et al. (Reference Liu, Wang, Vaidya, Ruehle, Halverson, Soljačić and Hou2024). The Kolmogorov–Arnold representation theorem states that any continuous multivariate function can be expressed as a composition of continuous univariate functions. KANs are based on this theorem, replacing the typical linear weight matrices in neural networks with learnable one-dimensional (1-D) functions (Liu et al. Reference Liu, Wang, Vaidya, Ruehle, Halverson, Soljačić and Hou2024). These functions are parameterised using splines, offering a flexible and computationally efficient approach to represent continuous functions. The KANs utilise B-splines, which are piecewise polynomial functions, to model local variations in data. Each spline segment corresponds to a polynomial function, and their piecewise nature allows KANs to approximate the local variations in the data during training. This adaptability enables KANs to capture intricate functional relationships more effectively, merging the advantages of B-splines with the traditional neural network framework, thereby enhancing both accuracy and interpretability.

In the TBKAN architecture, KAN replaces the hidden layers of the standard TBNN, while the anisotropy mapping portion remains unchanged. The TBKAN’s output layer is designed to predict the coefficients of the tensor basis.

Figure 7. Model architecture and training configuration for TBKAN which replaces the multi-layer perceptron in the modified TBNN (see figure 5) with a KAN.

1.5. Machine learning procedure

The KCNN, TBNN and TBKAN architectures proposed in § 1.4 were implemented in PyTorch, along with the proposed realisability-informed loss function (1.26). These models were trained using an open-source dataset for data-driven turbulence modelling (McConkey, Yee & Lien Reference McConkey, Yee and Lien2021). The objective of this study is to train and evaluate models trained on various flows, to determine whether the proposed realisability-informed loss function and architecture modifications are significantly beneficial. All code is available on Github (McConkey & Kalia Reference McConkey and Kalia2024).

1.5.1. Datasets

The training flows consist of flow over periodic hills (Xiao et al. Reference Xiao, Wu, Laizet and Duan2020), flow through a square duct (Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010) and flow over a flat plate with zero pressure gradient (Rumsey Reference Rumsey2021). These flows are selected because they contain several challenging physical phenomena for RANS, including separation, reattachment and Prandtl’s secondary flows. The flat plate case is also included, to demonstrate how machine learning can improve the anisotropy estimates within the boundary layer. For each flow type, both a hold-out validation set and a hold-out test set are selected. The validation set is used during training to help guide when to stop training in order to prevent overfitting, but the validation set loss is not back propagated through the network to update weights and biases. While we hold out an entire case for the test set (the usual procedure in data-driven turbulence modelling), we also generate the validation sets by holding out entire cases at a time. We recommend this method for generating validation sets in data-driven turbulence modelling – it is analogous to grouped cross-validation, a practice used in machine learning where several data points come from a single observation. Here, we consider each separate flow case as a single observation, each containing many data points. It is therefore prudent to ensure that two data points from the same observation are not used in both the training and test set.

Table 1 outlines the three training/validation/test splits considered. The objective in splitting the dataset this way is to determine whether a realisability-informed model can generalise to a new case for a given flow. Machine learning-based anisotropy mappings do not generalise well to entirely new flows (Duraisamy Reference Duraisamy2021; McConkey et al. Reference McConkey, Yee and Lien2022b ; Man et al. Reference Man, Jadidi, Keshmiri, Yin and Mahmoudi2023). However, they can be used to dramatically enhance the performance of a RANS simulation for a given flow type. In this same spirit, we aim to test how our modifications improve the generalisability of the learned anisotropy mapping to an unseen flow, albeit within the same class of flow.

Table 1. Datasets used for training, validation and testing. Varied parameters are defined in § 2.1.

1.5.2. Input features

The input features here are all derived from the baseline RANS $k$ - $\omega$ SST simulation. The input features form the vector $x^{{R}}_m$ . The superscript ${{R}}$ has been dropped in this section to avoid crowded notation, but it applies to all quantities discussed in § 1.5.2.

The input features must be Galilean invariant in order to generate an appropriately constrained anisotropy mapping. Most data-driven anisotropy mapping investigations use a mixture of heuristic scalars and scalars systematically generated from a minimal integrity basis for a set of gradient tensors (Wu et al. Reference Wu, Xiao and Paterson2018). We emphasise that all scalars must be Galilean invariant – without this criterion, the RANS equations will lose Galilean invariance. Despite the importance of this constraint, several data-driven anisotropy mappings include scalars like the turbulence intensity, which breaks Galilean invariance.

We use a mixture of heuristic scalars and scalars systematically generated from a minimal integrity basis (Wu et al. Reference Wu, Xiao and Paterson2018). We use the following heuristic scalars:

(1.36) \begin{align} q_1 &= {\textrm{min}}\left (\frac {\sqrt {k} y_w}{50 \nu },2\right )\!, \end{align}

(1.37) \begin{align} q_2 &= \frac {k}{\varepsilon } \sqrt {\sum _{i} \sum _{\!j} |S_{\textit{ij}}|^2}, \end{align}

(1.38) \begin{align} q_3 &= \frac {\sqrt {\displaystyle \sum _{i} \displaystyle \sum _{\!j} |\tau _{\textit{ij}}|^2}}{k} , \end{align}

(1.39) \begin{align} q_4 &= \frac {\sqrt {k}}{0.09\omega y_w} , \end{align}

(1.40) \begin{align} q_5 &= \frac {500\nu }{y_w^2 \omega } , \end{align}

(1.41) \begin{align} q_6 &= {\textrm{min}}\left ({\textrm{max}}\left (q_4,q_5\right )\!,\frac {2.0k}{{\textrm{max}}\left (\dfrac {y_w^2}{\omega }\dfrac {\partial k}{\partial x_i}\dfrac {\partial \omega }{\partial x_i}\right )}\right )\! , \\[6pt] \nonumber \end{align}

with $\varepsilon = 0.09 \omega$ , and $y_w$ is the distance to the nearest wall.

These input features correspond to: the wall-distance-based Reynolds number ( $q_1$ ), the ratio of turbulent time scale to mean strain time scale ( $q_2$ ), the ratio of total Reynolds stress to $k$ ( $q_3$ ) and different blending scalars used within the $k$ - $\omega$ SST model ( $q_4$ , $q_5$ , $q_6$ ). The full list of integrity basis tensors is given in Appendix A. We use the first invariant (the trace, $A_{\textit{ii}}$ ) of the following tensors: $B^{(1)}_{\textit{ij}}$ , $B^{(3)}_{\textit{ij}}$ , $B^{(4)}_{\textit{ij}}$ , $B^{(6)}_{\textit{ij}}$ , $B^{(7)}_{\textit{ij}}$ , $B^{(16)}_{\textit{ij}}$ and $B^{(35)}_{\textit{ij}}$ . We use the second invariant, ${1}/{2} ( (A_{\textit{ii}} )^2 - A_{\textit{ij}}A_{\textit{ji}} )$ , of the following tensors: $B^{(3)}_{\textit{ij}}$ , $B^{(6)}_{\textit{ij}}$ , $B^{(7)}_{\textit{ij}}$ , $B^{(8)}_{\textit{ij}}$ . These input features were hand picked from the full set of 94 invariants listed in Appendix A (Wu et al. Reference Wu, Xiao and Paterson2018; McConkey et al. Reference McConkey, Yee and Lien2022a ). As discussed in McConkey et al. (Reference McConkey, Yee and Lien2022a ), many of the invariants are zero for 2-D flows. However, different conditions cause different invariants to be zero. For example, in the present study, there are a different set of zero invariants for flow through over periodic hills, and flow through a square duct. This difference occurs because different components of $\partial () / \partial x_j$ are to be zero. We have performed a systematic investigation using a symbolic math toolbox (sympy (Meurer et al. Reference Meurer2017)) to determine which invariants are zero for the duct case, and general 2-D flows. We make the results and code available in Appendix A and on Github (McConkey Reference McConkey2023), respectively. Input features used in this investigation were selected based on the results in Appendix A. Therefore, the input features are not uniformly zero on any of the considered flows.

To ensure all input features are of the same magnitude during training, they are scaled according to the following formula:

(1.42) \begin{equation} x^{{R}}_m = \frac {X^{{R}}_m - \mu _{m}}{\sigma _m} , \end{equation}

where $x^{{R}}_m$ is the input feature vector for the neural network, $X^{{R}}_m$ is the raw input feature vector from the RANS simulation, $\mu _m$ is a vector containing the mean of each input feature over the entire training dataset and $\sigma _{m}$ is a vector containing the standard deviation of each input feature over the entire training dataset.

When making predictions on the hold-out validation and test sets, the mean and standard deviation values from the training data are used to avoid data leakage.

1.5.3. Hyperparameters and training procedure

The hyperparameters for each neural network were hand tuned based on validation set performance. The hidden layers for all TBNNs and KCNNs are fully connected, feed-forward layers with Swish activation functions (Ramachandran, Zoph & Le Reference Ramachandran, Zoph and Le2017). The appropriate hyperparameters vary between flows, since each dataset contains a different number of data points, and the anisotropy mapping being learned is distinct. Table 2 shows the final hyperparameters used. All training runs used a mini-batch size of 32, except the periodic hill TBNN run, which used a mini-batch size of 128.

Table 2. Hyperparameters selected for each model.

Figure 8. Loss function vs epoch for the (a) TBNN and (b) KCNN models on the flat plate dataset; the (c) TBNN and (d) KCNN models on the square duct dataset; and the (e) TBNN and (f) KCNN models on the periodic hill dataset.

The AMSGrad Adam optimiser (Reddi, Kale & Kumar Reference Reddi, Kale and Kumar2018) was found to achieve better performance than the standard Adam optimiser (Kingma & Ba Reference Kingma and Ba2015) for training TBNNs. The learning rate and number of epochs for each optimiser is given in table 2. Satisfactory performance was achieved with a constant learning rate; for training on more complex flows we recommend the use of learning rate scheduling to achieve better performance. The training/validation loss curves for each model are shown in figure 8.

1.5.4. Hyperparameter tuning for TBKAN

The performance of the TBKAN was extensively optimised through a combination of systematic hyperparameter tuning and manual adjustments. The architecture of the TBKAN model was configured as $[6, 9, 10]$ , where 6 corresponds to the number of input features, 9 represents the network width (number of neurons per hidden layer) and 10 denotes the output size of the network. The depth of the network, representing the number of hidden layers, was fixed to 2 for the flat plate case.

Hyperparameter tuning focused on refining the grid size ( $ g$ ), network width ( $ w$ ), spline order and input feature combinations. An initial set of 27 runs, conducted with randomly selected configurations, broadly explored the hyperparameter space. These runs were used as a warm start for Bayesian optimisation, which further utilised 143 trials to systematically refine the hyperparameters. The configuration employed a grid size of 8 control points for the B-spline basis representation, with the polynomial order $k$ of the splines fixed at a value of three (viz., $k=3$ corresponding to cubic splines). This choice was found to provide the best balance between flexibility and computational efficiency.

The final hyperparameters for the best-performing model were as follows: architecture $[6, 9, 10]$ , a depth of 2 layers, a grid size of 8 control points and a learning rate of $ 4.9 \times 10^{-3}$ . This configuration achieved a mean-squared error (MSE) of 0.22 on the flat plate case. Input feature selection was refined by fixing three core features while sampling from a broader set to enhance the model’s adaptability in predicting finer details of the anisotropy tensor. The AMSGrad Adam optimiser (Reddi et al. Reference Reddi, Kale and Kumar2018) was used for all training runs, with a mini-batch size of 32, which ensured stable convergence. Training and validation loss curves for the best-performing model for the flat plate, corresponding to a KAN configuration with $w=9$ , $g=8$ and $k=3$ (cubic splines), are shown in figure 9.

Figure 9. Loss function vs epoch for the TBKAN model trained on the flat plate dataset, with the best-performing KAN configuration given by a network width of $w=9$ , a grid size of $g=8$ and a polynomial order of $k=3$ (cubic B-spline).

1.6. Computational costs

The inference-time cost for the proposed methodology varies significantly from case to case. A prediction requires (i) running a baseline RANS simulation, (ii) evaluating the Machine learning (ML) model predictions and (iii) running a corrected RANS simulation. Compared with steps (i) and (iii), the cost of step (ii) (model inference) is negligible. Additionally, since the converged fields from step (i) are used to initialise the simulation in step (iii), the cost of the corrected RANS simulation is also reduced. Generally, we found that combined inference-time costs for steps (i), (ii) and (iii) were between 1.5 and 3 times the cost of the baseline simulation (i.e. steps (i), (ii) and (iii) cost approximately 1.5–5 times as much as step (i)), depending on how much the injected closure fields differ from the original linear eddy viscosity-based field. The training cost also varies depending on the dataset. For the periodic hills and square duct datasets here, the training time was approximately 16 GPU hours on a single NVIDIA RTX 3090 GPU.