{"id":174646,"date":"2025-09-22T19:59:45","date_gmt":"2025-09-22T19:59:45","guid":{"rendered":"https:\/\/www.newsbeep.com\/us\/174646\/"},"modified":"2025-09-22T19:59:45","modified_gmt":"2025-09-22T19:59:45","slug":"realisability-informed-machine-learning-for-turbulence-anisotropy-mappings-journal-of-fluid-mechanics","status":"publish","type":"post","link":"https:\/\/www.newsbeep.com\/us\/174646\/","title":{"rendered":"Realisability-informed machine learning for turbulence anisotropy mappings | Journal of Fluid Mechanics"},"content":{"rendered":"

It is prohibitively expensive to resolve all relevant scales of turbulence for industrially relevant flows. Even with increasing computational capacity, Kolmogorov microscale-resolving techniques such as direct numerical simulation (DNS) will be out of reach for decades (Slotnick et\u00a0al. Reference Slotnick, Khodadoust, Alonso, Darmofal, Gropp, Lurie and Marvriplis2014<\/a>). In order to enable the practical simulation of turbulent flows, a variety of techniques are currently in use, such as Reynolds-averaged Navier Stokes (RANS), detached eddy simulation (DES) and large eddy simulation (LES). Each technique comes with its own advantages. The user\u2019s available computational resources are the primary consideration. The RANS technique, which models turbulence as a single-scale phenomenon, remains the most popular industrial technique due to the computational costs of scale-resolving simulations (Witherden\u00a0& Jameson Reference Witherden and Jameson2017<\/a>).<\/p>\n

In recent years, a new menu item for industrial turbulence modelling is emerging. Machine learning has been used to augment RANS, DES, LES and DNS, with various objectives, such as accelerating simulation times (Kochkov et\u00a0al. Reference Kochkov, Smith, Alieva, Wang, Brenner and Hoyer2021<\/a>), inferring optimal coefficient fields (Singh, Duraisamy\u00a0& Zhang Reference Singh, Duraisamy and Zhang2017<\/a>), model calibration (Matai\u00a0& Durbin Reference Matai and Durbin2019<\/a>) and turbulence model augmentation (Duraisamy, Iaccarino\u00a0& Xiao Reference Duraisamy, Iaccarino and Xiao2019<\/a>; Brunton, Noack\u00a0& Koumoutsakos Reference Brunton, Noack and Koumoutsakos2020<\/a>). For LES, emerging research on training in situ and within a differentiable solver offers the prospect of a generalisable subgrid model (List, Chen\u00a0& Thuerey Reference List, Chen and Thuerey2022<\/a>; Sirignano\u00a0& MacArt Reference Sirignano and MacArt2023<\/a>). It is clear that for LES, allowing the model to interact with a solver for an unsteady flow vastly improves the generalisability of the learned closure relationship (List et\u00a0al. Reference List, Chen and Thuerey2022<\/a>). Within the context of steady-state RANS simulations, a promising new menu item is the ability to train a specialised turbulence model with a given industrial dataset. Significant attention has been given to the development of machine learning-augmented closure models. Specifically, flow-specific sensitisation of the stable but often inaccurate linear eddy viscosity relationship via machine learning is an area of major interest. In their seminal work, Ling, Kurzawski\u00a0& Templeton (Reference Ling, Kurzawski and Templeton2016<\/a>) proposed a neural network architecture based on a tensor basis expansion of the anisotropy tensor, referred to as a tensor basis neural network (TBNN). This architecture was extended to random forests by Kaandorp (Reference Kaandorp2018<\/a>) and Kaandorp\u00a0& Dwight (Reference Kaandorp and Dwight2020<\/a>). The TBNN architecture has been used in several studies, such as by Song et\u00a0al. (Reference Song, Zhang, Wang, Ye and Huang2019<\/a>) and Zhang et\u00a0al. (Reference Zhang, Song, Ye, Wang, Huang, An and Chen2019<\/a>). Further modifications to the TBNN framework within the context of simple channel flows have been proposed by Cai et\u00a0al. (Reference Cai, Angeli, Martinez, Damblin and Lucor2022<\/a>, Reference Cai, Angeli, Martinez, Damblin and Lucor2024<\/a>). Specifically, Cai et\u00a0al. studied issues related to non-unique mappings between the closure term and input features for plane channel flow, an issue reported by Liu et\u00a0al. (Reference Liu, Fang, Rolfo, Moulinec and Emerson2021<\/a>). This issue primarily occurs when a small input feature set is used, such as the 5 invariants (several of which are zero for two-dimensional flows) used in Ling et\u00a0al.\u2019s original TBNN (Ling et\u00a0al. Reference Ling, Kurzawski and Templeton2016<\/a>). Several strategies have been proposed to address the issue of a non-unique mapping, including incorporation of additional input features by Wu, Xiao\u00a0& Paterson (Reference Wu, Xiao and Paterson2018<\/a>), and ensembling through a divide and conquer approach by Man et\u00a0al. (Reference Man, Jadidi, Keshmiri, Yin and Mahmoudi2023<\/a>). In the present investigation, we address the issue related to non-uniqueness of the mapping via use of a rich input feature set. In an effort to make these mappings more transparent, Mandler\u00a0& Weigand (Reference Mandler and Weigand2023<\/a>) analysed predictions made via an anisotropy mapping using input feature importance metrics such as Shapley additive explanation values. Interpretability analysis sheds light on which RANS features contain more information about different flow physics, which helps guide future input feature selection (Mandler\u00a0& Weigand Reference Mandler and Weigand2023<\/a>). The TBNN-type models are not the only architecture in use\u00a0\u2013 others include the eigenvalue reconstruction technique proposed by Wu et\u00a0al. (Reference Wu, Xiao and Paterson2018<\/a>, Reference Wu, Xiao, Sun and Wang2019a
\n<\/a>), and the Reynolds force vector approach proposed by Cruz et\u00a0al. (Reference Cruz, Thompson, Sampaio and Bacchi2019<\/a>), and further investigated by Brener et\u00a0al. (Reference Brener, Cruz, Macedo and Thompson2022<\/a>) and Amarloo, Forooghi\u00a0& Abkar (Reference Amarloo, Forooghi and Abkar2022<\/a>).<\/p>\n

Here, we consider the specific problem of training data-driven anisotropy mappings between RANS and a higher fidelity closure term from LES or DNS. Several open questions remain in the area of machine learning anisotropy mappings for RANS, with perhaps the most popular question being whether a \u2018universal turbulence model\u2019 could be produced via machine learning (Duraisamy Reference Duraisamy2021<\/a>; Spalart Reference Spalart2023<\/a>). However, when it comes to generalisability of these mappings, the no-free-lunch theorem is at play as demonstrated in McConkey et\u00a0al. (Reference McConkey, Yee and Lien2022b
\n<\/a>). Accurate sensitisation of the anisotropy mapping for a given flow comes at the cost of generalisation to completely new flows. Nevertheless, numerous studies show that the sensitised mapping generalises well within the same flow type (Wu et\u00a0al. Reference Wu, Xiao and Paterson2018<\/a>; Kaandorp\u00a0& Dwight Reference Kaandorp and Dwight2020<\/a>; McConkey et\u00a0al. Reference McConkey, Yee and Lien2022b
\n<\/a>; Man et\u00a0al. Reference Man, Jadidi, Keshmiri, Yin and Mahmoudi2023<\/a>). While this lack of new flow generalisability means that a machine learning-augmented turbulence model is off the menu for some applications, it remains on the menu for many industrial applications. For example, an industrial user with LES data can leverage these techniques to develop an augmented RANS turbulence model sensitised to a flow of interest. Our opinion is that this lack of generalisability is due to a lack of sufficient data. Despite numerous datasets for this purpose being available, we believe the entire space of possible mappings between RANS input features and higher fidelity closure terms is still sparsely covered by available data. For the foreseeable future, these mappings will therefore be limited to flow-specific sensitisation.<\/p>\n

While the generalisability question is perhaps the most popular one, several other key issues remain in training anisotropy mappings for RANS. How can predicted quantities be injected in a stable and well-conditioned manner? How can eddy viscosity-based approaches be united with TBNN-type approaches? How can we incorporate alternatives to the traditional fully connected layers in a TBNN? In this investigation, we address several of these remaining questions. We formulate an injection framework which unites TBNN-type model architectures with the more stable optimal eddy viscosity-based techniques. This unification enables stable use of the equivariance-enforcing TBNN within the momentum equation, without the use of stabilising blending factors (Kaandorp\u00a0& Dwight Reference Kaandorp and Dwight2020<\/a>). The proposed framework also has the advantage of producing a well-conditioned solution without the use of an optimal eddy viscosity, an often-unstable quantity that is difficult to predict via a machine learning model. We also target the problem of producing realisable predictions via a TBNN-type architecture via a physics-based loss function penalty (viz., incorporation of a learning bias in the framework). To further improve the flexibility and representational capacity of TBNN-type models, we investigate the inclusion of a Kolmogorov\u2013Arnold network (KAN) (Liu et\u00a0al. Reference Liu, Wang, Vaidya, Ruehle, Halverson, Solja\u010di\u0107 and Hou2024<\/a>) into the framework to replace the multi-layer perception in the TBNN. We also further investigate the invariant input feature sets commonly used for TBNN architectures, and provide new insights as to which input features are appropriate for use in flows with certain zero gradient directions. We demonstrate good generalisation performance of the \u2018realisability-informed\u2019 TBNN. We also demonstrate that the realisability-informed loss function greatly reduces non-realisable predictions when generalising.<\/p>\n

The present work unites TBNN-type frameworks (for example Ling et\u00a0al. Reference Ling, Kurzawski and Templeton2016<\/a>; Kaandorp Reference Kaandorp2018<\/a>; Kaandorp\u00a0& Dwight Reference Kaandorp and Dwight2020<\/a>; and Man et\u00a0al. Reference Man, Jadidi, Keshmiri, Yin and Mahmoudi2023<\/a>) with optimal eddy viscosity frameworks (e.g. Wu et\u00a0al. Reference Wu, Xiao and Paterson2018<\/a>; Brener et\u00a0al. Reference Brener, Cruz, Thompson and Anjos2021<\/a>; and McConkey et\u00a0al. Reference McConkey, Yee and Lien2022a
\n<\/a>). The advantage here is maintaining a stable injection environment (Wu et\u00a0al. Reference Wu, Sun, Laizet and Xiao2019b
\n<\/a>; Brener et\u00a0al. Reference Brener, Cruz, Thompson and Anjos2021<\/a>), while also retaining the simplicity, elegance and implicit equivariance of the TBNN architecture. Additionally, the present work is the first to implement a way to inform the TBNN of physical realisability during the training process. Whereas most existing techniques to enforce realisability of the neural network involve an ad-hoc post-processing step, our technique leverages physical realisability as an additional training target, thereby embedding an additional physics-based (learning) bias into the model. This idea extends the learning bias proposed by Riccius, Agrawal\u00a0& Koutsourelakis (Reference Riccius, Agrawal and Koutsourelakis2023<\/a>). Lastly, despite widespread use of minimal integrity basis input features for two-dimensional flows, and flows through a square duct, the present investigation is the first to systematically examine these input features for flow through a square duct. This examination leads to several unsuitable input features being identified, and a codebase for investigating other flows of interest.<\/p>\n

This manuscript is organised as follows. Section 1<\/a> describes the novel techniques in detail, including the injection framework\u00a0((1.1<\/a>)\u2013(1.2<\/a>)), realisability-informed loss function (1.3<\/a>) and improved TBNN architecture\u00a0(1.4<\/a>). Details on the datasets, input features, hyperparameters, implementation and code availability are given in \u00a7 1.5<\/a>. Two primary results are presented in \u00a7 2<\/a>: generalisation tests of the realisability-informed TBNN\u00a0(2.1<\/a>), and an examination of how realisability-informed training produces more realisable predictions when generalising\u00a0(2.2<\/a>). Conclusions and future work are discussed in \u00a7 3<\/a>.<\/p>\n

1. Methodology<\/p>\n

1.1. Injecting Reynolds stress tensor decompositions used in data-driven closure modelling frameworks<\/p>\n

The steady-state, RANS momentum equations for an incompressible, Newtonian fluid are given by (Reynolds Reference Reynolds1895<\/a>)<\/p>\n

(1.1)<\/p>\n

<\/p>\n

\\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}<\/p>\n

where <\/p>\n

<\/p>\n

$U_{\\!j}$<\/p>\n

is the mean velocity, <\/p>\n

<\/p>\n

$P$<\/p>\n

is the mean pressure, <\/p>\n

<\/p>\n

$\\rho$<\/p>\n

is the fluid density, <\/p>\n

<\/p>\n

$\\nu$<\/p>\n

is the kinematic viscosity and <\/p>\n

<\/p>\n

$\\tau _{\\textit{ij}}$<\/p>\n

is the Reynolds stress tensor.<\/p>\n

The continuity equation also applies to this flow<\/p>\n

(1.2)<\/p>\n

<\/p>\n

\\begin{equation} \\frac {\\partial U_i}{\\partial x_i} = 0 . \\end{equation}<\/p>\n

Together,\u00a0(1.1<\/a>) and\u00a0(1.2<\/a>) are unclosed. In the most general case, there are four equations (continuity + 3 momentum), and 10 unknowns: <\/p>\n

<\/p>\n

$P$<\/p>\n

, 3 components of <\/p>\n

<\/p>\n

$U_i$<\/p>\n

and 6 components of <\/p>\n

<\/p>\n

$\\tau _{\\textit{ij}}$<\/p>\n

. The goal of turbulence closure modelling is to express <\/p>\n

<\/p>\n

$\\tau _{\\textit{ij}}$<\/p>\n

in terms of <\/p>\n

<\/p>\n

$P$<\/p>\n

and <\/p>\n

<\/p>\n

$U_i$<\/p>\n

.<\/p>\n

The Reynolds stress tensor can be decomposed into isotropic (hydrostatic) and anisotropic (deviatoric) components<\/p>\n

(1.3)<\/p>\n

<\/p>\n

\\begin{equation} \\tau _{\\textit{ij}} = \\frac {2}{3}k\\delta _{\\textit{ij}} + a_{\\textit{ij}} , \\end{equation}<\/p>\n

where <\/p>\n

<\/p>\n

$k = ({1}\/{2})\\tau _{kk}$<\/p>\n

is half the trace of the Reynolds stress tensor, and <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

is the anisotropy tensor. The isotropic component (<\/p>\n

<\/p>\n

$({2}\/{3})k\\delta _{\\textit{ij}}$<\/p>\n

) in\u00a0(1.3<\/a>) can be absorbed into the pressure gradient term in\u00a0(1.1<\/a>) to form a modified pressure<\/p>\n

(1.4)<\/p>\n

<\/p>\n

\\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}<\/p>\n

Equation\u00a0(1.4<\/a>) 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\u2019s divergence (as in Brener et\u00a0al. (Reference Brener, Cruz, Macedo and Thompson2022<\/a>)), which imply the use of\u00a0(1.1<\/a>). In our work, we model the Reynolds stress anisotropy tensor <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

, implying the use of\u00a0(1.4<\/a>). In either case, the term \u2018injection\u2019 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.<\/p>\n

In an eddy viscosity hypothesis, the anisotropy tensor <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

is postulated to be a function of the mean strain rate <\/p>\n

<\/p>\n

$S_{\\textit{ij}}$<\/p>\n

and rotation rate <\/p>\n

<\/p>\n

$R_{\\textit{ij}}$<\/p>\n

tensors<\/p>\n

(1.5)<\/p>\n

<\/p>\n

\\begin{align} a_{\\textit{ij}} &= a_{\\textit{ij}}(S_{\\textit{ij}}, R_{\\textit{ij}}) , \\end{align}<\/p>\n

(1.6)<\/p>\n

<\/p>\n

\\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}<\/p>\n

(1.7)<\/p>\n

<\/p>\n

\\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}<\/p>\n

The eddy viscosity hypothesis implies the following important constraints:<\/p>\n

(i) The anisotropy tensor can be predicted locally and instantaneously. In\u00a0(1.5<\/a>), there is no temporal dependence.<\/p>\n

(ii) The anisotropic turbulent stresses are caused entirely by mean velocity gradients.<\/p>\n

There are many examples of turbulent flows where these assumptions are violated (Pope Reference Pope2000<\/a>). For example, history effects occur in boundary layers (Bobke et\u00a0al. Reference Bobke, Vinuesa, \u00d6rl\u00fc and Schlatter2017<\/a>). Nevertheless these approximations remain widely used in eddy viscosity modelling. The most common eddy viscosity hypothesis is the linear eddy viscosity hypothesis<\/p>\n

(1.8)<\/p>\n

<\/p>\n

\\begin{equation} a_{\\textit{ij}} = – 2 \\nu _t S_{\\textit{ij}} , \\end{equation}<\/p>\n

where <\/p>\n

<\/p>\n

$\\nu _t$<\/p>\n

is the eddy viscosity, a scalar. This linear hypothesis draws direct analogy from the stress\u2013strain 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:<\/p>\n

(i) The anisotropy tensor is aligned with the mean strain rate tensor.<\/p>\n

(ii) The mapping between mean strain rate and anisotropic stress is isotropic, in that it can be represented using a single scalar (<\/p>\n

<\/p>\n

$\\nu _t$<\/p>\n

).<\/p>\n

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 <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

and <\/p>\n

<\/p>\n

$S_{\\textit{ij}}$<\/p>\n

, which occurs for even simple flows. Applying Cayley\u2013Hamilton theorem (Hamilton Reference Hamilton1853<\/a>; Cayley Reference Cayley1858<\/a>) to\u00a0(1.5<\/a>), Pope (Reference Pope1975<\/a>) derived the most general expression for a nonlinear eddy viscosity model<\/p>\n

(1.9)<\/p>\n

<\/p>\n

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

where <\/p>\n

<\/p>\n

$n=1,2,\\ldots 10$<\/p>\n

indexes the scalar coefficients <\/p>\n

<\/p>\n

$g_n$<\/p>\n

, and the following basis tensors:<\/p>\n

(1.10)<\/p>\n

<\/p>\n

\\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}<\/p>\n

The tensors <\/p>\n

<\/p>\n

$\\hat {S}_{\\textit{ij}}$<\/p>\n

and <\/p>\n

<\/p>\n

$\\hat {R}_{\\textit{ij}}$<\/p>\n

above are the non-dimensionalised strain rate and rotation rate tensors. In Pope\u2019s original work, these tensors are given by<\/p>\n

(1.11)<\/p>\n

<\/p>\n

\\begin{align} \\hat {S}_{\\textit{ij}} &= \\frac {k}{\\varepsilon } S_{\\textit{ij}} , \\end{align}<\/p>\n

(1.12)<\/p>\n

<\/p>\n

\\begin{align} \\hat {R}_{\\textit{ij}} &= \\frac {k}{\\varepsilon } R_{\\textit{ij}}. \\end{align}<\/p>\n

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.<\/p>\n

With Pope\u2019s general expression for the anisotropy tensor,\u00a0(1.4<\/a>) becomes<\/p>\n

(1.13)<\/p>\n

<\/p>\n

\\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}<\/p>\n

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\u00a0al.\u2019s TBNN investigation (Ling et\u00a0al. Reference Ling, Kurzawski and Templeton2016<\/a>),\u00a0(1.13<\/a>) was used. Kaandorp (Reference Kaandorp2018<\/a>) and Kaandorp\u00a0& Dwight (Reference Kaandorp and Dwight2020<\/a>) 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<\/a>). Additionally, other techniques to incorporate equivariance are cumbersome, requiring learning and predicting in an invariant eigenframe of some tensorial quantity (Wu et\u00a0al. Reference Wu, Xiao and Paterson2018<\/a>; Brener et\u00a0al. Reference Brener, Cruz, Macedo and Thompson2022<\/a>). 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 <\/p>\n

<\/p>\n

$S_{\\textit{ij}}$<\/p>\n

and <\/p>\n

<\/p>\n

$R_{\\textit{ij}}$<\/p>\n

can be used to represent the anisotropy tensor. However, a major disadvantage with numerically solving\u00a0(1.13<\/a>) is that the closure term is entirely explicit, greatly reducing numerical stability and conditioning (Brener et\u00a0al. Reference Brener, Cruz, Thompson and Anjos2021<\/a>). For this reason, Kaandorp\u00a0& Dwight (Reference Kaandorp and Dwight2020<\/a>) needed to implement a blending function, which blends the fully explicit closure term in\u00a0(1.13<\/a>) 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 <\/p>\n

<\/p>\n

$U_i$<\/p>\n

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 <\/p>\n

<\/p>\n

$a_{\\textit{ij}} = -2 \\nu _t S_{\\textit{ij}}$<\/p>\n

,\u00a0(1.4<\/a>) can be written as<\/p>\n

(1.14)<\/p>\n

<\/p>\n

\\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}<\/p>\n

Equation\u00a0(1.14<\/a>) has the major advantage of increasing diagonal dominance of the <\/p>\n

<\/p>\n

$U_{\\!j}$<\/p>\n

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.<\/p>\n

In the present work, we propose the following hybrid treatment of the closure term:<\/p>\n

(1.15)<\/p>\n

<\/p>\n

\\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}<\/p>\n

In\u00a0(1.15<\/a>), the <\/p>\n

<\/p>\n

$n=1$<\/p>\n

(linear) term has been separated and receives implicit treatment, while the remaining <\/p>\n

<\/p>\n

$n = 2,3,\\ldots 10$<\/p>\n

terms grant an opportunity for a machine learning model to provide rich representation of the nonlinear part of <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

. As we will discuss, the separation of the linear term requires special treatment during training and closure term injection.<\/p>\n

1.2. Conditioning analysis<\/p>\n

Various decompositions of the Reynolds stress tensor were investigated by Brener et\u00a0al.\u2019s conditioning analysis (Brener et\u00a0al. Reference Brener, Cruz, Thompson and Anjos2021<\/a>). 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\u00a0al. (Reference Brener, Cruz, Thompson and Anjos2021<\/a>) 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:<\/p>\n

(1.16)<\/p>\n

<\/p>\n

\\begin{equation} a_{\\textit{ij}} = -2\\nu ^{{R}}_{t} S^{\\theta}_{\\textit{ij}} + a^{\\perp}_{\\textit{ij}} \\end{equation}<\/p>\n

also achieves a well-conditioned solution. We follow the nomenclature of Duraisamy et\u00a0al. (Reference Duraisamy, Iaccarino and Xiao2019<\/a>) in that the <\/p>\n

<\/p>\n

$\\theta$<\/p>\n

superscript indicates a quantity that comes from a high-fidelity source such as DNS. The <\/p>\n

<\/p>\n

${{R}}$<\/p>\n

superscript indicates a quantity taken from the corresponding baseline RANS simulation.<\/p>\n

The decomposition in\u00a0(1.16<\/a>) 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 \u2018turned off\u2019 in unstable situations. It also is more interpretable \u2013 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\u00a0al. Reference McConkey, Yee and Lien2022a
\n<\/a>), this quantity is often unstable and difficult to predict via a machine learning model.<\/p>\n

<\/p>\n

<\/p>\n

\nFigures\u00a01<\/a> and 2<\/a> demonstrate a conditioning test similar to the tests conducted by Brener et\u00a0al. (Reference Brener, Cruz, Thompson and Anjos2021<\/a>). 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\u00a01<\/a> and 2<\/a>, we can see that the proposed decomposition of the Reynolds stress tensor\u00a0(1.16<\/a>) 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\u00a0al. Reference Brener, Cruz, Thompson and Anjos2021<\/a>). We further confirm the findings of Brener et\u00a0al. (Reference Brener, Cruz, Thompson and Anjos2021<\/a>) 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\u00a0al. (Reference Wu, Xiao, Sun and Wang2019a
\n<\/a>) that implicit treatment helps address the ill-conditioning issue, but using the DNS strain rate tensor <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

to calculate the explicitly injected <\/p>\n

<\/p>\n

$a^{\\perp}_{\\textit{ij}}$<\/p>\n

.<\/p>\n

1.3. Realisability-informed training<\/p>\n

The Reynolds stress tensor is symmetric positive semidefinite. A set of constraints on the non-dimensional anisotropy tensor <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

arise from this property, as determined by Banerjee et\u00a0al. (Reference Banerjee, Krahl, Durst and Zenger2007<\/a>). These constraints are<\/p>\n

(1.17)<\/p>\n

<\/p>\n

\\begin{align} -\\frac {1}{3} \\leqslant b_{\\textit{ij}} &\\leqslant \\frac {2}{3}, \\,\\, i= j, \\end{align}<\/p>\n

(1.18)<\/p>\n

<\/p>\n

\\begin{align} -\\frac {1}{2} \\leqslant b_{\\textit{ij}} &\\leqslant \\frac {1}{2}, \\,\\, i\\neq j, \\end{align}<\/p>\n

(1.19)<\/p>\n

<\/p>\n

\\begin{align} \\lambda _1 &\\geqslant \\frac {3|\\lambda _2| – \\lambda _2}{2}, \\end{align}<\/p>\n

(1.20)<\/p>\n

<\/p>\n

\\begin{align} \\lambda _1 &\\leqslant \\frac {1}{3}- \\lambda _2, \\\\[6pt] \\nonumber \\end{align}<\/p>\n

where the non-dimensional anisotropy tensor <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

is calculated by<\/p>\n

(1.21)<\/p>\n

<\/p>\n

\\begin{equation} b_{\\textit{ij}} = \\frac {a_{\\textit{ij}}}{2k} , \\end{equation}<\/p>\n

and the eigenvalues of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

are given by <\/p>\n

<\/p>\n

$\\lambda _1 \\geqslant \\lambda _2 \\geqslant \\lambda _3$<\/p>\n

.<\/p>\n

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 <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\varepsilon$<\/p>\n

model, <\/p>\n

<\/p>\n

$\\varepsilon$<\/p>\n

is the TKE dissipation rate (Launder\u00a0& Spalding Reference Launder and Spalding1974<\/a>), <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

(Wilcox Reference Wilcox1988<\/a>) and <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model (Menter Reference Menter1994<\/a>; Menter, Kuntz\u00a0& Langtry Reference Menter, Kuntz and Langtry2003<\/a>) do not guarantee physical realisability.<\/p>\n

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\u2019s 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\u00a0al. (Reference Ling, Kurzawski and Templeton2016<\/a>). 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\u00a0al. Reference Jiang, Vinuesa, Chen, Mi, Laima and Li2021<\/a>).<\/p>\n

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\u00a0& Karniadakis Reference Raissi, Perdikaris and Karniadakis2019<\/a>), we term the use of this loss function \u2018realisability-informed training\u2019. Whereas PINNs encourage the model to learn physics by penalising violations of conservation laws, realisability-informed TBNNs learn to predict physically realisable anisotropy tensors.<\/p>\n

The realisability penalty <\/p>\n

<\/p>\n

$\\mathcal{R}(b_{\\textit{ij}})$<\/p>\n

is given as follows:<\/p>\n

(1.22)<\/p>\n

<\/p>\n

\\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}<\/p>\n

where the <\/p>\n

<\/p>\n

$\\sim$<\/p>\n

above a symbol denotes a model prediction for the quantity associated with the symbol.<\/p>\n

Equation\u00a0(1.22<\/a>) can be thought of as the mean squared violation in the components of <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

, plus the mean squared violation in the eigenvalues of <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

. To help visualise this penalty function, figure\u00a03<\/a> shows the penalties incurred by violating various realisability constraints on the anisotropy tensor.<\/p>\n

<\/p>\n

Figure\u00a03. Realisability penalty for the (a) diagonal components of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

, and (b) off-diagonal components of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

.<\/p>\n

It should be noted that in the same way that a PINN\u2019s 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 <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

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 <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

value. The merits of this approach are demonstrated in \u00a7 2.2<\/a>.<\/p>\n

Another important metric in the loss function is the error-based penalty <\/p>\n

<\/p>\n

$\\mathcal{E}(\\tilde b_{\\textit{ij}})$<\/p>\n

. 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.<\/p>\n

(i) Since <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

is a symmetric tensor, we propose to sum the squared errors as follows:<\/p>\n

(1.23)<\/p>\n

<\/p>\n

\\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}<\/p>\n

Calculating the squared error in this way avoids double penalising the off-diagonal components, a situation which arises when summing over all components of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

.<\/p>\n

(ii) Although the TBNN model predicts <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

, we propose to use a loss function based on the error in <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

. Near the wall, <\/p>\n

<\/p>\n

$\\hat {T}^{(n)}_{\\textit{ij}}\\rightarrow 0$<\/p>\n

, but <\/p>\n

<\/p>\n

$b_{\\textit{ij}} \\nrightarrow 0$<\/p>\n

. The vanishing of <\/p>\n

<\/p>\n

$\\hat {T}^{(n)}_{\\textit{ij}}$<\/p>\n

with non-vanishing <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

causes instabilities during training, leading to <\/p>\n

<\/p>\n

$g_n \\rightarrow \\infty$<\/p>\n

here. Since <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

is the tensor injected into the momentum equation, its accurate prediction should be the focus of the training process. To dimensionalise <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

, the turbulent kinetic energy <\/p>\n

<\/p>\n

$k$<\/p>\n

must be used. This leads to the following error-based loss:<\/p>\n

(1.24)<\/p>\n

<\/p>\n

\\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}<\/p>\n

or stated more simply
\n(1.25)<\/p>\n

<\/p>\n

\\begin{equation} \\mathcal{E}(\\tilde a_{\\textit{ij}})=(2k^\\theta )^2 \\mathcal{E}(\\tilde b_{\\textit{ij}}) . \\end{equation}<\/p>\n

The reason for the use of <\/p>\n

<\/p>\n

$k^\\theta$<\/p>\n

in\u00a0(1.24<\/a>) and\u00a0(1.25<\/a>) is discussed in \u00a7 1.4<\/a>.<\/p>\n

The final loss function includes an error-based metric and a realisability-violation penalty. This loss function <\/p>\n

<\/p>\n

$\\mathcal{L}$<\/p>\n

is given by<\/p>\n

(1.26)<\/p>\n

<\/p>\n

\\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}<\/p>\n

where <\/p>\n

<\/p>\n

$\\alpha$<\/p>\n

is a factor used to control the relative importance of the realisability penalty. Here, <\/p>\n

<\/p>\n

$N$<\/p>\n

is the total number of points in the dataset. The points are indexed by <\/p>\n

<\/p>\n

$p = 1,2,\\ldots ,N$<\/p>\n

. The cases in the dataset are indexed <\/p>\n

<\/p>\n

$k=1,2,\\ldots ,s$<\/p>\n

. <\/p>\n

<\/p>\n

${\\mathcal I}_{C_k}(p)$<\/p>\n

is an indicator function<\/p>\n

(1.27)<\/p>\n

<\/p>\n

\\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}<\/p>\n

For a given point <\/p>\n

<\/p>\n

$p$<\/p>\n

, <\/p>\n

<\/p>\n

${\\mathcal I}_{C_k}(p)$<\/p>\n

selects all points which come from the same case as <\/p>\n

<\/p>\n

$p$<\/p>\n

. The denominator for all points from the same case <\/p>\n

<\/p>\n

$C_k$<\/p>\n

is the same \u2013 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<\/p>\n

(1.28)<\/p>\n

<\/p>\n

\\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}<\/p>\n

where <\/p>\n

<\/p>\n

$|C_k|$<\/p>\n

is the cardinality of the case <\/p>\n

<\/p>\n

$C_k$<\/p>\n

(viz., the number of points in the <\/p>\n

<\/p>\n

$k$<\/p>\n

th case <\/p>\n

<\/p>\n

$C_k$<\/p>\n

). Here, <\/p>\n

<\/p>\n

$Z_k$<\/p>\n

is simply the average of the Frobenius norm <\/p>\n

<\/p>\n

$\\| \\boldsymbol{\\cdot } \\|_F$<\/p>\n

of the anisotropy tensor <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

for all points <\/p>\n

<\/p>\n

$p$<\/p>\n

(viz., <\/p>\n

<\/p>\n

$a_{\\textit{ij}}(p)$<\/p>\n

) in case <\/p>\n

<\/p>\n

$C_k$<\/p>\n

. The objective of this denominator in\u00a0(1.26<\/a>) is to promote a more balanced regression problem, since data points from various cases or flow types may have <\/p>\n

<\/p>\n

$||a_{\\textit{ij}}||$<\/p>\n

that differ by orders of magnitude. Normalisation on a case-by-case basis is made in an effort to normalise all <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

error magnitudes to a similar scale.<\/p>\n

In this study, we use <\/p>\n

<\/p>\n

$\\alpha =10^2$<\/p>\n

to encode a high preference for realisable results in \u00a7 2<\/a>. Lower values of <\/p>\n

<\/p>\n

$\\alpha$<\/p>\n

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 <\/p>\n

<\/p>\n

$(2k^\\theta )^2$<\/p>\n

is used to formulate the loss function in terms of predicting <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

rather than <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

. However, <\/p>\n

<\/p>\n

$\\mathcal{R}(\\tilde b_{\\textit{ij}})$<\/p>\n

is also multiplied by <\/p>\n

<\/p>\n

$(2k^\\theta )^2$<\/p>\n

to ensure that the realisability penalty and mean-squared error in <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

are of similar scales.<\/p>\n

1.4. Neural network architecture<\/p>\n

Motivated by improving the training and injection stability, as well as generalisability, we propose several modifications to the original TBNN (Ling et\u00a0al. Reference Ling, Kurzawski and Templeton2016<\/a>).<\/p>\n

The original TBNN is shown in figure\u00a04<\/a>. At training time, this network predicts the non-dimensional anisotropy tensor <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

. All basis tensors used in this prediction at training time come from RANS, and during training the prediction <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

is evaluated against a known value of <\/p>\n

<\/p>\n

$b^\\theta _{\\textit{ij}}$<\/p>\n

from a high-fidelity simulation. At injection time, the network is used in this same configuration to predict <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

. Here, <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

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 <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

, to obtain an updated estimate for the turbulent kinetic energy <\/p>\n

<\/p>\n

$k$<\/p>\n

, and therefore an updated estimate for <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}} = 2 k \\tilde b_{\\textit{ij}}$<\/p>\n

.<\/p>\n

<\/p>\n

The modified TBNN is shown in figure\u00a05<\/a>. 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<\/p>\n

(1.29)<\/p>\n

<\/p>\n

\\begin{equation} \\hat {T}^{(1)}_{\\textit{ij}} = \\frac {k^{{R}}}{\\varepsilon ^{{R}}}S^{{R}}_{\\textit{ij}} , \\end{equation}<\/p>\n

our modified TBNN uses<\/p>\n

(1.30)<\/p>\n

<\/p>\n

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

<\/p>\n

Figure\u00a05. Model architecture and training configuration proposed in the present investigation. Note that during training, <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

is used in <\/p>\n

<\/p>\n

$\\hat {T}^{(1)}_{\\textit{ij}}$<\/p>\n

, whereas all other quantities come from the original RANS simulation.<\/p>\n

Further, while the original TBNN calculates the linear (denoted by superscript <\/p>\n

<\/p>\n

${\\textrm{L}}$<\/p>\n

) component of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

as<\/p>\n

(1.31)<\/p>\n

<\/p>\n

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

our modified TBNN uses<\/p>\n

(1.32)<\/p>\n

<\/p>\n

\\begin{equation} b^{\\textit{L}}_{\\textit{ij}} = – \\hat {T}^{(1)}_{\\textit{ij}} = – \\frac {\\nu _t}{k}S_{\\textit{ij}} , \\end{equation}<\/p>\n

where the superscript <\/p>\n

<\/p>\n

${{R}}$<\/p>\n

denotes a quantity that comes from the original RANS simulation.<\/p>\n

These changes are motivated by the following.<\/p>\n

(i) At injection time, we use implicit treatment of the linear term <\/p>\n

<\/p>\n

$\\hat {T}^{(1)}_{\\textit{ij}}$<\/p>\n

to formulate\u00a0(1.15<\/a>) in a stable manner. After injection, <\/p>\n

<\/p>\n

$S_{\\textit{ij}}$<\/p>\n

will continue to evolve. In a similar spirit as optimal eddy viscosity frameworks, we therefore use <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

to compute the linear component at training time. In optimal eddy viscosity based frameworks, using <\/p>\n

<\/p>\n

$S_{\\textit{ij}} = S^\\theta _{\\textit{ij}}$<\/p>\n

at training time helps drive <\/p>\n

<\/p>\n

$U_{\\!j}\\rightarrow U^\\theta _{\\!j}$<\/p>\n

at injection time (the cause of this behaviour is currently unknown). In Ling et\u00a0al. (Reference Ling, Kurzawski and Templeton2016<\/a>), all basis tensors (and therefore <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

) remain fixed after injection. We also fixed <\/p>\n

<\/p>\n

$\\nu _t$<\/p>\n

= <\/p>\n

<\/p>\n

$\\nu ^{{R}}_t$<\/p>\n

at injection time (viz., no further evolution of the eddy viscosity is permitted).<\/p>\n

(ii) At training time, we dimensionalise <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

using <\/p>\n

<\/p>\n

$k^\\theta$<\/p>\n

: <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}} = 2k^\\theta \\tilde b_{\\textit{ij}}$<\/p>\n

. At evaluation time, we do not have <\/p>\n

<\/p>\n

$k^\\theta$<\/p>\n

. However, the need for <\/p>\n

<\/p>\n

$k^\\theta$<\/p>\n

is avoided, since we only use the nonlinear part of <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

at test time. The reason we only need <\/p>\n

<\/p>\n

$\\tilde b^{\\perp}_{\\textit{ij}}$<\/p>\n

at test time is that the training process has been designed to use the linear part of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

estimated by the RANS turbulence model, and augment this by the TBNN\u2019s equivariant prediction for <\/p>\n

<\/p>\n

$\\tilde b^{\\perp}_{\\textit{ij}}$<\/p>\n

.<\/p>\n

(iii) Using <\/p>\n

<\/p>\n

$\\nu _t\/k$<\/p>\n

to normalise the basis tensors and fixing <\/p>\n

<\/p>\n

$g_1=-1$<\/p>\n

in\u00a0(1.32<\/a>) results in the RANS prediction for <\/p>\n

<\/p>\n

$b^{\\textit{L}}_{\\textit{ij}}$<\/p>\n

being implicitly used in the TBNN. Therefore, the TBNN learns to correct <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

using <\/p>\n

<\/p>\n

$b^{\\perp}_{\\textit{ij}}$<\/p>\n

in a way that allows a realisability-informed training process, and fully implicit treatment of the linear term at injection time.<\/p>\n

Lastly, the use of <\/p>\n

<\/p>\n

$k^\\theta$<\/p>\n

to dimensionalise <\/p>\n

<\/p>\n

$\\tilde b_{\\textit{ij}}$<\/p>\n

is also enabled by our use of a separate neural network to correct <\/p>\n

<\/p>\n

$k^{{R}}$<\/p>\n

at injection time. This neural network is called the <\/p>\n

<\/p>\n

$k$<\/p>\n

-correcting neural network (KCNN), and is a simple fully connected feed-forward neural network that predicts a single output scalar <\/p>\n

<\/p>\n

$\\varDelta$<\/p>\n

(Figure 6<\/a>)<\/p>\n

(1.33)<\/p>\n

<\/p>\n

\\begin{equation} \\varDelta = \\log \\left (\\frac {k^\\theta }{k^{{R}}}\\right )\\! , \\end{equation}<\/p>\n

such that at injection time, an updated estimate for <\/p>\n

<\/p>\n

$k$<\/p>\n

can be obtained <\/p>\n

<\/p>\n

$\\tilde k = {e}^\\varDelta k^{{R}}$<\/p>\n

, without the need to re-couple the turbulence transport equations. The KCNN shares the same input features as the TBNN.<\/p>\n

<\/p>\n

Figure\u00a06. Test-time injection configuration proposed in the present investigation. The hatching over <\/p>\n

<\/p>\n

$g_1$<\/p>\n

and <\/p>\n

<\/p>\n

$\\hat {T}^{(1)}_{\\textit{ij}}$<\/p>\n

indicates that they are not used at injection time.<\/p>\n

Together, the KCNN and TBNN predict the anisotropy tensor in the following manner:<\/p>\n

(1.34)<\/p>\n

<\/p>\n

\\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}<\/p>\n

With the linear component of <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}}$<\/p>\n

being treated implicitly during injection, the entire closure framework is summarised as explicit injection of the following term into the momentum equation:<\/p>\n

(1.35)<\/p>\n

<\/p>\n

\\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}<\/p>\n

1.4.1. Tensor basis Kolmogorov\u2013Arnold network (TBKAN)<\/p>\n

The tensor basis KAN shown in figure\u00a07<\/a> in the training configuration replaces the multi-layer perceptron in the modified TBNN with a KAN introduced by Liu et\u00a0al. (Reference Liu, Wang, Vaidya, Ruehle, Halverson, Solja\u010di\u0107 and Hou2024<\/a>). The Kolmogorov\u2013Arnold 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\u00a0al. Reference Liu, Wang, Vaidya, Ruehle, Halverson, Solja\u010di\u0107 and Hou2024<\/a>). 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.<\/p>\n

In the TBKAN architecture, KAN replaces the hidden layers of the standard TBNN, while the anisotropy mapping portion remains unchanged. The TBKAN\u2019s output layer is designed to predict the coefficients of the tensor basis.<\/p>\n

<\/p>\n

Figure\u00a07. Model architecture and training configuration for TBKAN which replaces the multi-layer perceptron in the modified TBNN (see figure\u00a05<\/a>) with a KAN.<\/p>\n

1.5. Machine learning procedure<\/p>\n

The KCNN, TBNN and TBKAN architectures proposed in \u00a7 1.4<\/a> were implemented in PyTorch, along with the proposed realisability-informed loss function\u00a0(1.26<\/a>). These models were trained using an open-source dataset for data-driven turbulence modelling (McConkey, Yee\u00a0& Lien Reference McConkey, Yee and Lien2021<\/a>). 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\u00a0& Kalia Reference McConkey and Kalia2024<\/a>).<\/p>\n

1.5.1. Datasets<\/p>\n

The training flows consist of flow over periodic hills (Xiao et\u00a0al. Reference Xiao, Wu, Laizet and Duan2020<\/a>), flow through a square duct (Pinelli et\u00a0al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010<\/a>) and flow over a flat plate with zero pressure gradient (Rumsey Reference Rumsey2021<\/a>). These flows are selected because they contain several challenging physical phenomena for RANS, including separation, reattachment and Prandtl\u2019s 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\u00a0\u2013 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.<\/p>\n

\nTable\u00a01<\/a> 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<\/a>; McConkey et\u00a0al. Reference McConkey, Yee and Lien2022b
\n<\/a>; Man et\u00a0al. Reference Man, Jadidi, Keshmiri, Yin and Mahmoudi2023<\/a>). 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.<\/p>\n

Table\u00a01. Datasets used for training, validation and testing. Varied parameters are defined in \u00a7 2.1<\/a>.<\/p>\n

<\/p>\n

1.5.2. Input features<\/p>\n

The input features here are all derived from the baseline RANS <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST simulation. The input features form the vector <\/p>\n

<\/p>\n

$x^{{R}}_m$<\/p>\n

. The superscript <\/p>\n

<\/p>\n

${{R}}$<\/p>\n

has been dropped in this section to avoid crowded notation, but it applies to all quantities discussed in \u00a7 1.5.2<\/a>.<\/p>\n

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\u00a0al. Reference Wu, Xiao and Paterson2018<\/a>). We emphasise that all scalars must be Galilean invariant \u2013 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.<\/p>\n

We use a mixture of heuristic scalars and scalars systematically generated from a minimal integrity basis (Wu et\u00a0al. Reference Wu, Xiao and Paterson2018<\/a>). We use the following heuristic scalars: <\/p>\n

(1.36)<\/p>\n

<\/p>\n

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

(1.37)<\/p>\n

<\/p>\n

\\begin{align} q_2 &= \\frac {k}{\\varepsilon } \\sqrt {\\sum _{i} \\sum _{\\!j} |S_{\\textit{ij}}|^2}, \\end{align}<\/p>\n

(1.38)<\/p>\n

<\/p>\n

\\begin{align} q_3 &= \\frac {\\sqrt {\\displaystyle \\sum _{i} \\displaystyle \\sum _{\\!j} |\\tau _{\\textit{ij}}|^2}}{k} , \\end{align}<\/p>\n

(1.39)<\/p>\n

<\/p>\n

\\begin{align} q_4 &= \\frac {\\sqrt {k}}{0.09\\omega y_w} , \\end{align}<\/p>\n

(1.40)<\/p>\n

<\/p>\n

\\begin{align} q_5 &= \\frac {500\\nu }{y_w^2 \\omega } , \\end{align}<\/p>\n

(1.41)<\/p>\n

<\/p>\n

\\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}<\/p>\n

with <\/p>\n

<\/p>\n

$\\varepsilon = 0.09 \\omega$<\/p>\n

, and <\/p>\n

<\/p>\n

$y_w$<\/p>\n

is the distance to the nearest wall.<\/p>\n

These input features correspond to: the wall-distance-based Reynolds number (<\/p>\n

<\/p>\n

$q_1$<\/p>\n

), the ratio of turbulent time scale to mean strain time scale (<\/p>\n

<\/p>\n

$q_2$<\/p>\n

), the ratio of total Reynolds stress to <\/p>\n

<\/p>\n

$k$<\/p>\n

(<\/p>\n

<\/p>\n

$q_3$<\/p>\n

) and different blending scalars used within the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model (<\/p>\n

<\/p>\n

$q_4$<\/p>\n

, <\/p>\n

<\/p>\n

$q_5$<\/p>\n

, <\/p>\n

<\/p>\n

$q_6$<\/p>\n

). The full list of integrity basis tensors is given in Appendix\u00a0A<\/a>. We use the first invariant (the trace, <\/p>\n

<\/p>\n

$A_{\\textit{ii}}$<\/p>\n

) of the following tensors: <\/p>\n

<\/p>\n

$B^{(1)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(3)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(4)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(6)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(7)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(16)}_{\\textit{ij}}$<\/p>\n

and <\/p>\n

<\/p>\n

$B^{(35)}_{\\textit{ij}}$<\/p>\n

. We use the second invariant, <\/p>\n

<\/p>\n

${1}\/{2} ( (A_{\\textit{ii}} )^2 – A_{\\textit{ij}}A_{\\textit{ji}} )$<\/p>\n

, of the following tensors: <\/p>\n

<\/p>\n

$B^{(3)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(6)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(7)}_{\\textit{ij}}$<\/p>\n

, <\/p>\n

<\/p>\n

$B^{(8)}_{\\textit{ij}}$<\/p>\n

. These input features were hand picked from the full set of 94 invariants listed in Appendix\u00a0A<\/a> (Wu et\u00a0al. Reference Wu, Xiao and Paterson2018<\/a>; McConkey et\u00a0al. Reference McConkey, Yee and Lien2022a
\n<\/a>). As discussed in McConkey et\u00a0al. (Reference McConkey, Yee and Lien2022a
\n<\/a>), 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 <\/p>\n

<\/p>\n

$\\partial () \/ \\partial x_j$<\/p>\n

are to be zero. We have performed a systematic investigation using a symbolic math toolbox (sympy (Meurer et\u00a0al. Reference Meurer2017<\/a>)) to determine which invariants are zero for the duct case, and general 2-D flows. We make the results and code available in Appendix\u00a0A<\/a> and on Github (McConkey Reference McConkey2023<\/a>), respectively. Input features used in this investigation were selected based on the results in Appendix\u00a0A<\/a>. Therefore, the input features are not uniformly zero on any of the considered flows.<\/p>\n

To ensure all input features are of the same magnitude during training, they are scaled according to the following formula:<\/p>\n

(1.42)<\/p>\n

<\/p>\n

\\begin{equation} x^{{R}}_m = \\frac {X^{{R}}_m – \\mu _{m}}{\\sigma _m} , \\end{equation}<\/p>\n

where <\/p>\n

<\/p>\n

$x^{{R}}_m$<\/p>\n

is the input feature vector for the neural network, <\/p>\n

<\/p>\n

$X^{{R}}_m$<\/p>\n

is the raw input feature vector from the RANS simulation, <\/p>\n

<\/p>\n

$\\mu _m$<\/p>\n

is a vector containing the mean of each input feature over the entire training dataset and <\/p>\n

<\/p>\n

$\\sigma _{m}$<\/p>\n

is a vector containing the standard deviation of each input feature over the entire training dataset.<\/p>\n

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.<\/p>\n

1.5.3. Hyperparameters and training procedure<\/p>\n

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\u00a0& Le Reference Ramachandran, Zoph and Le2017<\/a>). 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\u00a02<\/a> 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.<\/p>\n

Table\u00a02. Hyperparameters selected for each model.<\/p>\n

<\/p>\n

<\/p>\n

Figure\u00a08. 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.<\/p>\n

The AMSGrad Adam optimiser (Reddi, Kale\u00a0& Kumar Reference Reddi, Kale and Kumar2018<\/a>) was found to achieve better performance than the standard Adam optimiser (Kingma\u00a0& Ba Reference Kingma and Ba2015<\/a>) for training TBNNs. The learning rate and number of epochs for each optimiser is given in table\u00a02<\/a>. 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\u00a08<\/a>.<\/p>\n

1.5.4. Hyperparameter tuning for TBKAN<\/p>\n

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 <\/p>\n

<\/p>\n

$[6, 9, 10]$<\/p>\n

, 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.<\/p>\n

Hyperparameter tuning focused on refining the grid size (<\/p>\n

<\/p>\n

$ g$<\/p>\n

), network width (<\/p>\n

<\/p>\n

$ w$<\/p>\n

), 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 <\/p>\n

<\/p>\n

$k$<\/p>\n

of the splines fixed at a value of three (viz., <\/p>\n

<\/p>\n

$k=3$<\/p>\n

corresponding to cubic splines). This choice was found to provide the best balance between flexibility and computational efficiency.<\/p>\n

The final hyperparameters for the best-performing model were as follows: architecture <\/p>\n

<\/p>\n

$[6, 9, 10]$<\/p>\n

, a depth of 2 layers, a grid size of 8 control points and a learning rate of <\/p>\n

<\/p>\n

$ 4.9 \\times 10^{-3}$<\/p>\n

. 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\u2019s adaptability in predicting finer details of the anisotropy tensor. The AMSGrad Adam optimiser (Reddi et\u00a0al. Reference Reddi, Kale and Kumar2018<\/a>) 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 <\/p>\n

<\/p>\n

$w=9$<\/p>\n

, <\/p>\n

<\/p>\n

$g=8$<\/p>\n

and <\/p>\n

<\/p>\n

$k=3$<\/p>\n

(cubic splines), are shown in figure\u00a09<\/a>.<\/p>\n

<\/p>\n

Figure\u00a09. 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 <\/p>\n

<\/p>\n

$w=9$<\/p>\n

, a grid size of <\/p>\n

<\/p>\n

$g=8$<\/p>\n

and a polynomial order of <\/p>\n

<\/p>\n

$k=3$<\/p>\n

(cubic B-spline).<\/p>\n

1.6. Computational costs<\/p>\n

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\u20135 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.<\/p>\n

2. Results for proposed TBNN and TBKAN architectures<\/p>\n

2.1. Generalisation tests<\/p>\n

It was of interest to determine how well the trained models generalise to unseen variations of their training flows. This section demonstrates generalisation results for flow over a flat plate with zero pressure gradient, flow through a square duct and flow over periodic hills.<\/p>\n

For all cases, the original RANS solution was generated using OpenFOAM v2212, assuming an isothermal, incompressible and Newtonian fluid. Simulation parameters such as solver, schemes and solution methodology for the zero pressure gradient flat plate case were identical to those discussed in McConkey et\u00a0al. (Reference McConkey, Yee and Lien2021<\/a>).<\/p>\n

It should be noted that the <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

predictions shown for the TBNN\/KCNN use <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

for predicting the linear part of the anisotropy tensor. The nature of the proposed TBNN training process is to utilise <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

during training, so that the remaining nonlinear part can be extracted during injection.<\/p>\n

2.1.1. Flat plate (TBNN)<\/p>\n

<\/p>\n

Figure\u00a010. Computational domain for the zero pressure gradient flat plate boundary layer case.<\/p>\n

This case features a developing turbulent boundary layer on a flat plate with zero pressure gradient, based on the NASA \u20182DZP\u2019 validation case (Rumsey Reference Rumsey2021<\/a>). Figure\u00a010<\/a> shows the domain for the flat plate case. The NASA-provided meshes are sufficient to resolve the viscous sublayer region, with a total number of cells <\/p>\n

<\/p>\n

$ N \\approx 2\\,00\\,000$<\/p>\n

. However, this mesh was further refined to increase the number of solution data points available for training and testing. The goal of this case is to learn how the anisotropy tensor evolves in a turbulent boundary layer, therefore substantial mesh refinement was required to generate data points in this region. The total number of cells in the mesh is <\/p>\n

<\/p>\n

$N=4\\,673\\,130$<\/p>\n

. The plate-length Reynolds number is <\/p>\n

<\/p>\n

${\\textit{Re}}_L= 5(10)^6$<\/p>\n

to match the NASA reference data. The reference data for this case consist of a series of wall-normal profiles, for various <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta$<\/p>\n

, defined as<\/p>\n

(2.1)<\/p>\n

<\/p>\n

\\begin{equation} Re_\\theta = \\frac {U_\\infty \\theta }{\\nu } , \\end{equation}<\/p>\n

where the momentum thickness <\/p>\n

<\/p>\n

$\\theta$<\/p>\n

is given as<\/p>\n

(2.2)<\/p>\n

<\/p>\n

\\begin{equation} \\theta = \\int _{0}^\\infty \\frac {U_{1}}{U_\\infty } \\left (1-\\frac {U_1}{U_\\infty }\\right ) {\\textrm{d}} x_2 , \\end{equation}<\/p>\n

and <\/p>\n

<\/p>\n

$U_\\infty$<\/p>\n

is the free-stream velocity. The following boundary conditions and fluid properties are used for the domain in figure\u00a010<\/a>. At the inlet boundary, <\/p>\n

<\/p>\n

$U_{\\!j}$<\/p>\n

, <\/p>\n

<\/p>\n

$k$<\/p>\n

and <\/p>\n

<\/p>\n

$\\omega$<\/p>\n

are uniform: <\/p>\n

<\/p>\n

$U_{\\!j} = (69.4, 0, 0) {\\textrm{m s}}^{-1}$<\/p>\n

, <\/p>\n

<\/p>\n

$k=1.08(10)^{-3} \\ {\\textrm{m}}^2\\, {\\textrm{s}}^{-2}$<\/p>\n

, <\/p>\n

<\/p>\n

$\\omega = 8675 \\ {\\textrm{s}}^{-1}$<\/p>\n

and <\/p>\n

<\/p>\n

$P$<\/p>\n

is zero normal gradient. At the outlet, <\/p>\n

<\/p>\n

$P$<\/p>\n

is zero, and all other variables are zero normal gradient. At the symmetry plane, all variables are zero normal gradient, and normal velocity is zero. At the top plane, all flow variables are zero normal gradient. At the no-slip wall, <\/p>\n

<\/p>\n

$U_{\\!j}=0$<\/p>\n

, <\/p>\n

<\/p>\n

$k=0$<\/p>\n

, <\/p>\n

<\/p>\n

$\\omega =6\\nu \/\\beta _1 y^2$<\/p>\n

(here, <\/p>\n

<\/p>\n

$\\beta _1=0.075$<\/p>\n

) and <\/p>\n

<\/p>\n

$P$<\/p>\n

is zero normal gradient. A kinematic viscosity of <\/p>\n

<\/p>\n

$\\nu = 1.388(10)^{-5} \\ {\\textrm{m}}\\, {\\textrm{s}}^{-2}$<\/p>\n

was used.<\/p>\n

The DNS reference data for a developing turbulent boundary layer come from Schlatter\u00a0& \u00d6rl\u00fc (Reference Schlatter and \u00d6rl\u00fc2010<\/a>). This dataset contains a variety of turbulent boundary layer profiles, at various <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta$<\/p>\n

shown in table\u00a01<\/a>. As discussed in \u00a7 1.5.1<\/a>, the TBNN model was trained on various <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta$<\/p>\n

, with <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta =3630$<\/p>\n

serving as a hold-out test case. The results shown in this section are for this hold-out test case.<\/p>\n

The baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model performs well for the test case flow in terms of predicting the mean velocity profiles. Figure\u00a011<\/a> shows a sample mean velocity profile predicted by the baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, and figure\u00a012<\/a> shows the predicted evolution of <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta$<\/p>\n

along the plate. The excellent performance of the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model demonstrated by figures\u00a011<\/a> and 12<\/a> is expected. This case features a fully attached boundary layer with zero pressure gradient, which is one of the fundamental calibration scenarios for RANS models. The zero pressure gradient turbulent boundary layer is considered a \u2018solved problem\u2019 for modern RANS models (Spalart Reference Spalart2023<\/a>).<\/p>\n

<\/p>\n

<\/p>\n

Figure\u00a012. Momentum thickness Reynolds number growth along the flat plate as predicted by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, and the reference DNS data from Schlatter\u00a0& \u00d6rl\u00fc (Reference Schlatter and \u00d6rl\u00fc2010<\/a>).<\/p>\n

While the mean velocity profile is predicted well, figure\u00a013<\/a> shows that the evolution of the near-wall anisotropy tensor is not predicted well. For this reason, the model architecture discussed in \u00a7 1.4<\/a> was used to correct the anisotropy tensor in the near-wall region. This test also aims to determine whether the input feature set is sufficiently expressive to enable predicting the evolution of the anisotropy tensor within a boundary layer. Figure\u00a013<\/a> shows the wall-normal profiles of various anisotropy tensor components predicted by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, the ML-augmented <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model and the reference DNS simulation for the hold-out test case. As discussed in \u00a7 1.5.1<\/a>, these models were trained on flat plate data at various values of <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta$<\/p>\n

. The results in figure\u00a013<\/a> are designed to test these models on input features from an unseen boundary layer profile, to determine whether the learned anisotropy mapping was generalisable.<\/p>\n

<\/p>\n

<\/p>\n

\nFigure\u00a013<\/a>(a) shows the predicted evolution of <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

in the turbulent boundary layer. The baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predicts <\/p>\n

<\/p>\n

$a_{11}=0$<\/p>\n

, since <\/p>\n

<\/p>\n

$\\partial U_1 \/\\partial x_1 \\approx 0$<\/p>\n

in the boundary layer. However, the DNS data clearly show that <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

is non-zero in the boundary layer. The TBNN\/KCNN model combination is able to correct the <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

term to a high degree of accuracy in the boundary layer on this test case, indicating that the anisotropy mapping for the <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

component generalises well. Similar evolutions of <\/p>\n

<\/p>\n

$a_{22}$<\/p>\n

(figure\u00a013<\/a>
\nc) and <\/p>\n

<\/p>\n

$a_{33}$<\/p>\n

(figure\u00a013<\/a>
\nd) are observed in the DNS data. Again, the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predicts <\/p>\n

<\/p>\n

$a_{11}=a_{22}=a_{33}=0$<\/p>\n

, which is not physically correct. The TBNN\/KCNN models are able to correct the baseline prediction to a high degree of accuracy on this unseen boundary layer profile.<\/p>\n

\nFigure\u00a013<\/a>(b) shows the predicted evolution of <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

. The baseline RANS model predicts the evolution of <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

well, and this is likely the reason that the mean velocity profile of <\/p>\n

<\/p>\n

$U_{1}$<\/p>\n

is predicted well (see figure\u00a011<\/a>). While the TBNN\/KCNN is not needed to correct this off-diagonal component, it is able to correct minor inaccuracies in the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predictions in the buffer region (<\/p>\n

<\/p>\n

$5 \\leqslant y^+ \\leqslant 30$<\/p>\n

). Nevertheless, the baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model achieves a satisfactory accuracy level for this flow. As discussed, this is expected, given that low Reynolds number RANS models are able to predict a zero pressure gradient boundary layer with a high degree of accuracy. Figure\u00a013<\/a> demonstrates that this performance is the result of an accurate prediction of <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model.<\/p>\n

2.1.2. Flat plate (TBKAN)<\/p>\n

The predicted evolution of the various components of the anisotropy tensor for the flat plate case obtained using the TBKAN\/KCNN model combination is displayed in figure\u00a014<\/a>(a)\u201314<\/a>(d) for the hold-out test case <\/p>\n

<\/p>\n

${\\textit{Re}}_\\theta = 3630$<\/p>\n

. These predictions are compared with both the baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model and DNS data from Schlatter\u00a0& \u00d6rl\u00fc (2010). For this test case, the predictive accuracy of the TBKAN\/KCNN model combination for the anisotropy tensor components is compared with that of the baseline model.<\/p>\n

The predicted evolution of <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

is shown in figure\u00a014<\/a>(a). Unlike the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, which predicts <\/p>\n

<\/p>\n

$a_{11}=0$<\/p>\n

due to its linear stress assumptions, the TBKAN captures the non-zero nature of <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

in the turbulent boundary layer. The TBKAN aligns closely with the DNS data, indicating its ability to generalise and represent anisotropy accurately, particularly for components where the baseline model is limited.<\/p>\n

\nFigure\u00a014<\/a>(b) presents the predicted evolution of <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

. The baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predicts this off-diagonal component with reasonable accuracy. However, the TBKAN\/KCNN exhibits better conformance with the reference DNS data in the buffer region (<\/p>\n

<\/p>\n

$5 \\leqslant y^+ \\leqslant 30$<\/p>\n

), correcting minor deviations in the predictions of this quantity provided by the baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model.<\/p>\n

\nFigures\u00a014<\/a>(c) and 14<\/a>(d) display the predictions of <\/p>\n

<\/p>\n

$a_{22}$<\/p>\n

and <\/p>\n

<\/p>\n

$a_{33}$<\/p>\n

, respectively. While the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predicts a value of zero from these two anisotropy stress components, the TBKAN\/KCNN captures their evolution with good accuracy compared with the DNS data. Overall, TBKAN\/KCNN enhances the predictive accuracy of the anisotropic stress components across the boundary layer. While <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

predictions are marginally improved compared with the baseline model, TBKAN\/KCNN shows substantial improvements for the normal stress components. Furthermore, a visual perusal of figures\u00a013<\/a> and 14<\/a> shows that the conformance of the predictions of the anisotropy tensor components with the reference DNS data obtained with TBKAN\/KCNN is marginally worse than that obtained with TBNN\/KCNN. However, it is noted that the TBKAN is simpler than TBNN in this case in the sense that the former network used only one hidden layer with 9 nodes (where the information from the edges encoded in the B-splines are simply accumulated), whereas the latter network used four hidden layers consisting of 20 nodes each (where the information embodied by the linear weights in the edges is transformed by the nonlinear activation function).<\/p>\n

2.1.3. Square duct<\/p>\n

Turbulent flow through a square duct is a challenging case for RANS models, since linear eddy viscosity models cannot predict the secondary flows that occur in the cross-sectional plane. The goal of the square duct test case is to determine whether the proposed closure framework could enable the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model to predict these Prandtl secondary flows (Nikitin, Popelenskaya\u00a0& Stroh Reference Nikitin, Popelenskaya and Stroh2021<\/a>). The square duct DNS dataset generated by Pinelli et\u00a0al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010<\/a>) was used as reference data, and the RANS data from McConkey et\u00a0al. (Reference McConkey, Yee and Lien2021<\/a>) was used.<\/p>\n

<\/p>\n

Figure\u00a015. Computational (a) domain and (b) RANS mesh for the square duct flow.<\/p>\n

\nFigure\u00a015<\/a> shows the computational set-up and mesh for the square duct case. The mesh is designed to achieve <\/p>\n

<\/p>\n

$y^+\\leqslant 1$<\/p>\n

for all square duct cases. As discussed in \u00a7 1.5.1<\/a>, the duct half-height <\/p>\n

<\/p>\n

$H$<\/p>\n

Reynolds number varies between cases, calculated by<\/p>\n

(2.3)<\/p>\n

<\/p>\n

\\begin{equation} Re_H = \\frac {U_b H}{\\nu } , \\end{equation}<\/p>\n

where <\/p>\n

<\/p>\n

$U_b$<\/p>\n

is the mean (bulk) cross-sectional velocity. A kinematic viscosity of <\/p>\n

<\/p>\n

$\\nu =5(10)^{-6}$<\/p>\n

m<\/p>\n

<\/p>\n

$^2$<\/p>\n

s\u22121 was used for all square duct cases. With the geometry fixed as shown in figure\u00a015<\/a>, the bulk velocity was adjusted to vary the Reynolds number. More details on the computational set-up for the square duct case are provided by McConkey et\u00a0al. (Reference McConkey, Yee and Lien2021<\/a>). The boundary conditions are periodic at the inlet\/outlet, and no-slip walls were applied along the sides of the duct. As discussed in \u00a7 1.5.1<\/a>, the modified TBNN was trained on several values of <\/p>\n

<\/p>\n

${\\textit{Re}}_{\\!H}$<\/p>\n

, with <\/p>\n

<\/p>\n

${\\textit{Re}}_{\\!H} = 2000$<\/p>\n

serving as a hold-out test case.<\/p>\n

<\/p>\n

\nFigure\u00a016<\/a> shows the components of the anisotropy tensor <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

predicted by RANS, DNS and the TBNN\/KCNN models for the square duct test case. The <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model is a linear eddy viscosity model, and therefore predicts zero <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

where <\/p>\n

<\/p>\n

$S_{\\textit{ij}}$<\/p>\n

is zero. Figure\u00a016<\/a> shows that <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

, <\/p>\n

<\/p>\n

$a_{22}$<\/p>\n

, <\/p>\n

<\/p>\n

$a_{23}$<\/p>\n

and <\/p>\n

<\/p>\n

$a_{33}$<\/p>\n

are all non-zero in the duct, and that the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model is unable to capture this behaviour. The TBNN\/KCNN models predict an accurate evolution of almost all anisotropy tensor components across the duct cross-section (viz. <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

, <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

, <\/p>\n

<\/p>\n

$a_{13}$<\/p>\n

, <\/p>\n

<\/p>\n

$a_{22}$<\/p>\n

and <\/p>\n

<\/p>\n

$a_{33}$<\/p>\n

are all predicted well on this test case). The anisotropy tensor component <\/p>\n

<\/p>\n

$a_{23}$<\/p>\n

is not predicted well, likely because it is at least an order of magnitude smaller than the other components, and therefore errors in <\/p>\n

<\/p>\n

$a_{23}$<\/p>\n

are not penalised as heavily in the loss function.<\/p>\n

\nFigure\u00a017<\/a> shows the turbulent kinetic energy <\/p>\n

<\/p>\n

$k$<\/p>\n

after being corrected by the KCNN model for the square duct test case. Accurate prediction of <\/p>\n

<\/p>\n

$k$<\/p>\n

is critical to an accurate estimate of <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

, since <\/p>\n

<\/p>\n

$a_{\\textit{ij}}=2kb_{\\textit{ij}}$<\/p>\n

. The <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model generally under-predicts <\/p>\n

<\/p>\n

$k$<\/p>\n

. After correction via the KCNN, the <\/p>\n

<\/p>\n

$k$<\/p>\n

field is predicted well compared with the DNS data. The primary feature in the <\/p>\n

<\/p>\n

$k$<\/p>\n

field that is absent from the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST prediction is the high-<\/p>\n

<\/p>\n

$k$<\/p>\n

region along the sidewalls of the duct. The KCNN introduces a correction to the baseline RANS field, and is able to predict this high-<\/p>\n

<\/p>\n

$k$<\/p>\n

region.<\/p>\n

Ultimately, it is the goal of the proposed framework to improve the estimated mean fields in the RANS simulation. To determine whether the corrected closure term would produce corrected mean velocity fields, the predicted <\/p>\n

<\/p>\n

$\\tilde a^{\\perp}_{\\textit{ij}}$<\/p>\n

was injected into the RANS momentum equation as shown in\u00a0(1.15<\/a>). The momentum and continuity equations converged around the fixed <\/p>\n

<\/p>\n

$\\tilde a^{\\perp}_{\\textit{ij}}$<\/p>\n

until numerical convergence was achieved. In OpenFOAM v2212, a modified version of the Pressure implicit splitting of operators – semi implicit method for pressure linked equations (PIMPLE) solver was implemented for the purpose of this injection. The PIMPLE solver was used to incorporate an unsteady term into the system of equations during iteration, to promote stability. Although this unsteady term affects the solution during convergence, the simulation ultimately achieved a steady-state condition, thereby reducing this unsteady term to zero.<\/p>\n

\nFigure\u00a018<\/a> shows that the TBNN\/KCNN model is able to produce secondary flows after injecting <\/p>\n

<\/p>\n

$\\tilde a^{\\perp}_{\\textit{ij}}$<\/p>\n

into the momentum equation. This a posteriori prediction of the mean field is ultimately the main prediction of interest for a ML-augmented RANS closure framework. Whereas the original <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model does not predict formation of any secondary flows in the duct, figure\u00a018<\/a> shows that the ML-augmented <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predicts corner vortices. However, the in-plane kinetic energy is generally underpredicted by the ML-augmented RANS model.<\/p>\n

<\/p>\n

Figure\u00a018. Velocity vectors and in-plane kinetic energy predicted by (a) the DNS data from Pinelli et\u00a0al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010<\/a>), (b) the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model and (c) injecting the TBNN\/KCNN predictions into the RANS momentum equations. Shown here is the lower left <\/p>\n

<\/p>\n

$x_2\\leqslant 0, \\ x_3\\leqslant 0$<\/p>\n

quadrant.<\/p>\n

To further examine the ability of the ML-augmented <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model to predict secondary flows in the duct test case, profiles of <\/p>\n

<\/p>\n

$U_2$<\/p>\n

and <\/p>\n

<\/p>\n

$U_3$<\/p>\n

are plotted in figure\u00a019<\/a>. While the ML-augmented model is able to produce this nonlinear feature, the corner vortex strengths are reduced compared with the reference DNS data. Both the <\/p>\n

<\/p>\n

$U_2$<\/p>\n

and <\/p>\n

<\/p>\n

$U_3$<\/p>\n

components are under-predicted. Nevertheless, the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model (which predicts <\/p>\n

<\/p>\n

$U_2=U_3=0$<\/p>\n

) has clearly been improved via a ML-augmented correction to the closure term in the momentum equation. From figure\u00a016<\/a>, it would appear that the good prediction of the normal stress anisotropy (the primary mechanism responsible for the streamwise vorticity determined by <\/p>\n

<\/p>\n

$U_2$<\/p>\n

and <\/p>\n

<\/p>\n

$U_3$<\/p>\n

) should provide good predictions of the streamwise vorticity. However, once the secondary flow is set in motion by this normal stress anisotropy, it is the secondary (rather than primary) shear stress component <\/p>\n

<\/p>\n

$a_{23}$<\/p>\n

(generated by the presence of the secondary flow itself) that is required to maintain this flow and from figure\u00a016<\/a>, this secondary component of the shear stress is not well predicted. Therefore, the underprediction of <\/p>\n

<\/p>\n

$U_2$<\/p>\n

and <\/p>\n

<\/p>\n

$U_3$<\/p>\n

is likely due to inaccurate prediction of <\/p>\n

<\/p>\n

$a_{23}$<\/p>\n

.<\/p>\n

<\/p>\n

2.1.4. Rectangular duct test case<\/p>\n

Duct aspect ratio has an influence on the in-plane kinetic energy <\/p>\n

<\/p>\n

${1}\/{2} (U_2^2+U_3^2 )$<\/p>\n

and behaviour of the corner vortices (Vinuesa et\u00a0al. Reference Vinuesa, Noorani, Lozano-Dur\u00e1n, Khoury, Schlatter, Fischer and Nagib2014<\/a>, Reference Vinuesa, Schlatter and Nagib2015<\/a>, Reference Vinuesa, Prus, Schlatter and Nagib2016<\/a>). It was of interest to examine the ability of a model trained on square ducts to generalise to a rectangular duct at similar Reynolds number. The TBNN\/KCNN models were applied to a duct with aspect ratio 3 and <\/p>\n

<\/p>\n

${\\textit{Re}}_{\\!H}\\approx 2600$<\/p>\n

, matching the DNS simulation by Vinuesa, Schlatter\u00a0& Nagib (Reference Vinuesa, Schlatter and Nagib2018<\/a>). This generalisation test was carried out in the same manner as the <\/p>\n

<\/p>\n

${\\textit{Re}}_{\\!H}=2000$<\/p>\n

hold-out test case (i.e. by making a predictive correction to the closure terms in the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST turbulence model, and injecting them back into the momentum equation).<\/p>\n

<\/p>\n

Figure\u00a020. Velocity vectors and in-plane kinetic energy predicted by (a) the DNS data from Pinelli et\u00a0al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010<\/a>), (b) the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model and (c) injecting the TBNN\/KCNN predictions into the RANS momentum equations. Shown here is the lower left <\/p>\n

<\/p>\n

$x_2\\leqslant 0, \\ x_3\\leqslant 0$<\/p>\n

quadrant.<\/p>\n

\nFigure\u00a020<\/a> compares the in-plane kinetic energy and velocity vector field predicted by the baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, and the TBNN\/KCNN augmented <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model. Examining figure\u00a020<\/a> shows that ML-based correction enables the augmented SST model to predict secondary flows. However, in the same manner as the square duct test case, the in-plane kinetic energy is underpredicted. Additionally, the strength of the dominant corner vortex is over-predicted by the ML-augmented model, while the strength of the smaller corner vortex is significantly underpredicted. The magnitude of the in-plane kinetic energy is similar, however, the peak locations of this field are also shifted due to a mismatch in predicted corner vortex shape. Nevertheless, the ML-augmented model shows improvement compared with the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, which fails to predict any corner vortices. In order to promote better generalisation to rectangular geometries, a more extensive training dataset consisting of rectangular ducts would need to be used to train the models.<\/p>\n

2.1.5. Periodic hills<\/p>\n

Flow over periodic hills is used as a popular benchmark case for turbulence modelling given the challenging physics of boundary layer separation in an adverse pressure gradient, reattachment along the bottom wall, and acceleration of the flow before re-entering the domain. For the purpose of data-driven turbulence modelling, a variety of periodic hill data have been made available. In this study, we use the <\/p>\n

<\/p>\n

${\\textit{Re}}_{\\!H}=5600$<\/p>\n

configuration, simulated using DNS by Xiao et\u00a0al. (Reference Xiao, Wu, Laizet and Duan2020<\/a>). Xiao et\u00a0al.\u2019s data were included in McConkey et\u00a0al. (Reference McConkey, Yee and Lien2021<\/a>), which is the primary data source for this study.<\/p>\n

<\/p>\n

Figure\u00a021. Computational domain and mesh for the periodic hill case, with <\/p>\n

<\/p>\n

$\\alpha =1.2$<\/p>\n

.<\/p>\n

The geometry and mesh for the periodic hill case are shown in figure\u00a021<\/a>. For all periodic hill cases, the hill height-based Reynolds number is 5600, calculated by<\/p>\n

(2.4)<\/p>\n

<\/p>\n

\\begin{equation} Re_H = \\frac {U_b H}{\\nu } , \\end{equation}<\/p>\n

where <\/p>\n

<\/p>\n

$U_b$<\/p>\n

is the bulk (mean) velocity at the domain inlet. The hill geometry is varied between cases, based on the hill steepness <\/p>\n

<\/p>\n

$\\alpha$<\/p>\n

. Further details on the computational set-up for the baseline RANS periodic hill simulations are provided in McConkey et\u00a0al. (Reference McConkey, Yee and Lien2021<\/a>). The boundary conditions are periodic at the inlet\/outlet, and no-slip walls at the top and bottom of the domain were imposed. As discussed in \u00a7 1.5.1<\/a>, the TBNN and KCNN models were trained on several hill steepness values, with <\/p>\n

<\/p>\n

$\\alpha =1.2$<\/p>\n

being used as a hold-out test set.<\/p>\n

<\/p>\n

Figure\u00a022. Contours of non-zero <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

components predicted by the DNS data from Xiao et\u00a0al. (Reference Xiao, Wu, Laizet and Duan2020<\/a>), the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model and the TBNN\/KCNN model predictions from the present investigation.<\/p>\n

\nFigure\u00a022<\/a> shows the components of the anisotropy tensor <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

predicted by RANS (<\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST), DNS and the ML-augmented RANS simulation. The improvement in all components of <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

is clear. The baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model under-predicts all components, with severe under-prediction of <\/p>\n

<\/p>\n

$a_{11}$<\/p>\n

, <\/p>\n

<\/p>\n

$a_{22}$<\/p>\n

and <\/p>\n

<\/p>\n

$a_{33}$<\/p>\n

. The <\/p>\n

<\/p>\n

$a_{12}$<\/p>\n

prediction by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model showcases similar trends to the DNS data, but the overall magnitude is lower. However, after correction, key features of all <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

fields are captured when the TBNN\/KCNN augment the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model. In particular, the higher magnitudes of the diagonal <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

(normal stress) components are captured by the augmented model.<\/p>\n

As was done for the square duct test case (\u00a7 2.1.3<\/a>), a modified PIMPLE solver was used to inject the predicted <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}}$<\/p>\n

into the RANS momentum equation for the periodic hill test case. The numerical set-up for the periodic hill injection was identical to the square duct case. Figure\u00a023<\/a> compares the mean velocity fields before and after the corrected closure term is used within the RANS simulation. Figure\u00a024<\/a> compares the errors in the velocity components <\/p>\n

<\/p>\n

$U_1$<\/p>\n

and <\/p>\n

<\/p>\n

$U_2$<\/p>\n

estimated by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model, and the a posteriori (post-injection) TBNN\/KCNN-augmented SST model.<\/p>\n

As seen in figure\u00a023<\/a>, the primary feature of this flow is a recirculation zone which appears immediately after the left hill. The recirculation zone is most clearly visualised by examining the <\/p>\n

<\/p>\n

$U_1$<\/p>\n

fields. The <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model over-predicts the size of this recirculation zone. After correction via injecting the TBNN\/KCNN predictions, the recirculation zone size closely matches the DNS data. In the <\/p>\n

<\/p>\n

$U_2$<\/p>\n

field, a region with <\/p>\n

<\/p>\n

$U_2\\lt 0$<\/p>\n

is seen immediately above this recirculation region. The baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model under-predicts the downward velocity here, leading to delayed reattachment, and a longer recirculation zone. After correction, the magnitude of <\/p>\n

<\/p>\n

$U_2$<\/p>\n

in this shear layer more closely matches the DNS data. On the right hill, the upward acceleration of the flow under the favourable pressure gradient is also under-predicted by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST. Here, the injected <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}}$<\/p>\n

is able to better capture the strength of this upward acceleration.<\/p>\n

\nFigure\u00a024<\/a> more closely examines the improvements offered by augmenting the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model via the TBNN\/KCNN. It can be seen that the overall magnitudes of the errors in <\/p>\n

<\/p>\n

$U_1$<\/p>\n

and <\/p>\n

<\/p>\n

$U_2$<\/p>\n

are significantly reduced after injecting the TBNN\/KCNN model predictions. In particular, error in <\/p>\n

<\/p>\n

$U_1$<\/p>\n

is reduced in the reattachment region along the bottom wall, and the bulk flow above this region. Error in <\/p>\n

<\/p>\n

$U_2$<\/p>\n

is reduced in the previously identified shear layer above the recirculation region, and the accelerating region before the outlet.<\/p>\n

2.2. Impact of realisability-informed training<\/p>\n

To determine the impact that including a realisability-informed penalty has on the training process, the closure term predictions for the three test cases were examined in greater detail. Two loss functions were used: one with <\/p>\n

<\/p>\n

$\\alpha =0$<\/p>\n

(representing no realisability penalty), and one with <\/p>\n

<\/p>\n

$\\alpha =10^2$<\/p>\n

(representing an exaggerated realisability penalty term in the loss function). The objective of this test was to determine whether including the realisability penalty during training promotes better generalisation of the closure mapping to unseen flow variations.<\/p>\n

All TBNN hyperparameters were fixed to those given in table\u00a02<\/a>. Since the realisability-informed training procedure only applies to the TBNN, a perfect prediction of <\/p>\n

<\/p>\n

$k$<\/p>\n

via the KCNN was assumed for calculating error in <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}}$<\/p>\n

.<\/p>\n

\nTable\u00a03<\/a> compares the MSE in <\/p>\n

<\/p>\n

$\\tilde a_{\\textit{ij}}$<\/p>\n

on the hold-out test set, with and without realisability penalties being used in training. It should be noted that similar to the a priori tests in \u00a7 2<\/a>, the linear component used when visualising <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

comes from <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

, as is the configuration when training the TBNN. This linear component is the one used when training the TBNN, and therefore its use here provides the most fair assessment of how the proposed loss function promotes more physically realisable results.<\/p>\n

Table\u00a03. Comparison of mean squared error and number of non-realisable predictions when training with and without a realisability-informed loss function.<\/p>\n

<\/p>\n

As seen in table\u00a03<\/a>, the realisability-informed loss function significantly reduces realisability violations on unseen flow variations. Even on hold-out test cases, the model is able to predict <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

without any realisability violations. In some cases, a small tradeoff in <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

occurs \u2013 this tradeoff is expected, as for some difficult points the realisability-informed loss function involves a tradeoff between error and realisability. However, in all cases, realisability-informed training also reduces the error in <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

. This error reduction in <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

is expected, since all <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

reference data are realisable. In the case of predicting <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

accurately, the gradients of the realisability penalty <\/p>\n

<\/p>\n

$\\mathcal{R}(b_{\\textit{ij}})$<\/p>\n

further push the predictions towards an accurate prediction of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

, compared with a purely MSE gradient.<\/p>\n

To further visualise the realisability of the TBNN predictions, the barycentric map of Banerjee et\u00a0al. (Reference Banerjee, Krahl, Durst and Zenger2007<\/a>) was used. In the barycentric map, the eigenvalues of a given <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

are mapped into a triangle, the bounds of which represent the limiting realisable behaviours of turbulent fluctuations. This triangle is useful to spatially visualise realisability violations, and types of turbulent flow physics predicted by the baseline turbulence model, the ML-augmented turbulence model and DNS. Further details on the construction of this mapping are given in Banerjee et\u00a0al. (Reference Banerjee, Krahl, Durst and Zenger2007<\/a>) and Emory\u00a0& Iaccarino (Reference Emory and Iaccarino2014<\/a>).<\/p>\n

\nFigure\u00a025<\/a> compares the realisability of the predictions made by a model trained on only a MSE loss function with a model trained on the realisability-informed loss function for the three flows in the present study. All predictions are for the hold-out test set for each flow, representing a generalisation test. The goal of the realisability-informed loss function is to encode a preference for realisable predictions into the model when generalising outside of the training dataset. Figure\u00a025<\/a> shows a clear improvement in the realisability of the TBNN predictions. This visualisation supports the results in table 3 in that a realisability-informed model has significantly lower realisability violations when predicting the anisotropy tensor on a new flow. Nearly all of the predictions from the realisability-informed model fall within the realisable boundaries, while the model trained only on MSE predicts several realisability-violating results when generalising to new cases of complex flows such as the duct and periodic hill cases. Also, the violations of physical realisability for the realisability-informed model (when they do occur) are not as severe (as measured from the magnitude of deviation outside the barycentric map) as those obtained from only a MSE loss function.<\/p>\n

\nFigure\u00a025<\/a> also shows that for all flows in the present study, the original <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model predicts plane-strain turbulence. All flows in the present study have a strain rate tensor which results in at least one zero eigenvalue of <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

, for the linear eddy viscosity approximation (<\/p>\n

<\/p>\n

$b_{\\textit{ij}} = -\\nu _t\/k S_{\\textit{ij}}$<\/p>\n

), and therefore plane-strain turbulence is predicted for all flows by the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model.<\/p>\n

<\/p>\n

Figure\u00a025. Projection of the predicted <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

from the DNS reference data (a,
\ne,i), the <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST model (b, f, j) and the TBNN\/KCNN with (d,h,l) and without (c,g,k) realisability-informed learning onto the barycentric triangle. Panels (a)\u2013(d) show the flat plate data, (e)\u2013(h) show the square duct data and (i)\u2013(l) show the periodic hill data. Points outside the triangle are not realisable.<\/p>\n

As discussed in \u00a7 1.3<\/a>, realisability-informed learning function does not guarantee a fully realisable prediction. Rather, realisability-informed learning encodes a preference for realisable predictions into the model predictions. This avoids the need for ad hoc post-processing of the predicted anisotropy tensor, while also encoding a physics-based (learning) bias into the TBNN. Strict realisability can be further enforced by post-processing any non-realisable predictions by the TBNN, as the anisotropy tensor can still be accessed and assessed for realisability in the proposed framework (albeit, with an evolving linear component).<\/p>\n

An important distinction between the dimensional and non-dimensional anisotropy tensor can also be drawn from the results in this section. In turbulence modelling, it is often the non-dimensional anisotropy tensor <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

that is thought to be of interest \u2013 indeed, it is possible to non-dimensionalise the closure problem (e.g. the generalised eddy viscosity model proposed by Pope (Reference Pope1975<\/a>)). However, in the present study, we show that even with physically realisable DNS data, dimensionalising the anisotropy tensor via <\/p>\n

<\/p>\n

$a_{\\textit{ij}}=2kb_{\\textit{ij}}$<\/p>\n

provides a distinct learning target and loss landscape. The realisability-informed loss function does not, in principle, compete against a simple error-based loss function in terms of predicting the non-dimensional <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

\u2013 non-realisable predictions will also have high error. However, setting the dimensional anisotropy tensor as the target creates a tradeoff between these two objectives, as demonstrated by the results in table\u00a03<\/a>. Since predicting the dimensional anisotropy tensor is favourable for the reasons outlined in \u00a7 1.3<\/a>, further investigation is required to determine the exact source of this tradeoff. This future research area is particularly relevant for learning anisotropy mappings from data, a major area of focus for improving RANS via ML.<\/p>\n

3. Conclusion<\/p>\n

The objectives of this study were to propose a physics-informed loss function for training TBNN anisotropy mappings, a new TBNN architecture and a new injection framework that accommodates implicit treatment of the linear anisotropy component within a TBNN-type architecture. This framework addresses an issue related to the stability of injected TBNN predictions, an issue that has been reported by Kaandorp\u00a0& Dwight (Reference Kaandorp and Dwight2020<\/a>), and consistent with our own experience. This framework also addresses the issue of realisability within the context of predicting the anisotropy tensor from RANS input features.<\/p>\n

The results here indicate that the proposed model architecture generalises well to new flow configurations, and the predicted anisotropy tensor can be injected in a highly stable manner. During the injection procedure in the present investigation, the Computational Fluid Dynamics solver remained stable, even when testing erratic model predictions. While this finding corroborates findings from others that frameworks which leverage implicit treatment of the linear anisotropy component via an eddy viscosity are numerically stable (Wu et\u00a0al. Reference Wu, Sun, Laizet and Xiao2019b
\n<\/a>; Brener et\u00a0al. Reference Brener, Cruz, Thompson and Anjos2021<\/a>; Liu et\u00a0al. Reference Liu, Fang, Rolfo, Moulinec and Emerson2021<\/a>), a novelty in the present investigation is the avoidance of using an optimal eddy viscosity. With the proposed decomposition of <\/p>\n

<\/p>\n

$a_{\\textit{ij}}$<\/p>\n

(\u00a7 1.1<\/a>), a simpler (and more straightforward) baseline <\/p>\n

<\/p>\n

$k$<\/p>\n

–<\/p>\n

<\/p>\n

$\\omega$<\/p>\n

SST eddy viscosity is used. The core modification that allows a well-conditioned solution in this case is the use of <\/p>\n

<\/p>\n

$S^\\theta _{\\textit{ij}}$<\/p>\n

within <\/p>\n

<\/p>\n

$\\hat {T}^{(1)}_{\\textit{ij}}$<\/p>\n

in the TBNN training process. As a result, all corrections induced by the TBNN are contained in an explicit term in the momentum equation. Future work includes investigating how new flow generalisation can be improved by a blending factor that multiplies this explicit term, thereby allowing the correction to be turned off. For example, a statistics-based scalar or a separate machine learning model could be used to predict a blending factor. This blending factor could reduce or eliminate corrections when a test data point departs significantly (is out of distribution) from the training dataset, and erroneous predictions are likely.<\/p>\n

While the realisability-informed loss function does not strictly guarantee physical realisability of the predictions, the results in \u00a7 2.2<\/a> indicate that the model retains a realisability bias when generalising. The realisability-informed loss function is not only applicable to the framework and architecture proposed in this study \u2013 it could be used anytime an anisotropy mapping is generated via machine learning. The use of the realisability-informed loss function in the present investigation promoted better realisability of predictions by a TBNN, but we fully expect the bias induced by this loss function to also be beneficial for tensor basis random forests (e.g. Kaandorp\u00a0& Dwight (Reference Kaandorp and Dwight2020<\/a>)), or non-tensor basis frameworks (e.g. Wu et\u00a0al. (Reference Wu, Xiao and Paterson2018<\/a>)). Since high-quality DNS anisotropy tensor data are realisable, the realisability penalty can be viewed as an additional boost to the loss function gradient towards the true value, when a non-realisable prediction is made. Future work will investigate how the realisability-informed loss function performs with training data generated by LES, which are not guaranteed to be realisable.<\/p>\n

The KCNN used in the proposed framework to correct <\/p>\n

<\/p>\n

$k$<\/p>\n

could also be replaced by coupling the <\/p>\n

<\/p>\n

$k$<\/p>\n

equation to the momentum equation, as exists in the original turbulence model, and the TBNN approach originally proposed by Ling et\u00a0al. (Reference Ling, Kurzawski and Templeton2016<\/a>). At training time, <\/p>\n

<\/p>\n

$k^\\theta$<\/p>\n

is used to dimensionalise <\/p>\n

<\/p>\n

$b_{\\textit{ij}}$<\/p>\n

in the present study, since <\/p>\n

<\/p>\n

$k$<\/p>\n

is corrected via the KCNN. This direct correction produced satisfactory results in the present work, but it is possible that generalisation could be further enhanced by the re-coupling the <\/p>\n

<\/p>\n

$k$<\/p>\n

equation with an updated closure term. This enhanced generalisation would result from the fact that a physics-based coupled equation system is used to correct <\/p>\n

<\/p>\n

$k$<\/p>\n

, rather than a simple multiplicative corrector (the KCNN in the present investigation). However, this coupling of an additional partial differential equation introduces the possibility for instability, an issue which is common in machine learning anisotropy modelling. Future work will investigate the merits of this route.<\/p>\n

Ultimately, this investigation demonstrates that with sufficient modifications, TBNN-type anisotropy mappings can be injected in a stable and well-conditioned manner. Further, with appropriate physics-based loss function penalties, the mapping can be sensitised to more physically informative targets than MSE. Moreover, we provided a preliminary investigation of the utility of KANs for turbulence closure modelling. While TBKAN did not outperform TBNN in this context, they nevertheless demonstrated the potential for capturing the complex relationships in the anisotropy tensor (at least in the flat plate case), suggesting that future research work should be conducted on the use of KANs for turbulent closure modelling applications. While industrial use of machine learning-based anisotropy mappings is currently not widespread, the continual development of techniques which increase the practicality of training and injecting model predictions will help accelerate more widespread use. Machine learning-augmented turbulence closure modelling is an aid that the turbulence modelling community can use to help bridge the current computational gap between RANS and widespread use of LES (Witherden\u00a0& Jameson Reference Witherden and Jameson2017<\/a>).<\/p>\n","protected":false},"excerpt":{"rendered":"It is prohibitively expensive to resolve all relevant scales of turbulence for industrially relevant flows. Even with increasing…\n","protected":false},"author":2,"featured_media":174647,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[45],"tags":[182,181,507,74],"class_list":{"0":"post-174646","1":"post","2":"type-post","3":"status-publish","4":"format-standard","5":"has-post-thumbnail","7":"category-artificial-intelligence","8":"tag-ai","9":"tag-artificial-intelligence","10":"tag-artificialintelligence","11":"tag-technology"},"_links":{"self":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/posts\/174646","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/comments?post=174646"}],"version-history":[{"count":0,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/posts\/174646\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/media\/174647"}],"wp:attachment":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/media?parent=174646"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/categories?post=174646"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/tags?post=174646"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}