1. Introduction
Turbulence is a multiscale phenomenon, and resolving all turbulent scales directly remains prohibitively expensive at high, practically relevant Reynolds numbers (Choi & Moin Reference Choi and Moin2012; Yang & Griffin Reference Yang and Griffin2021; Li, Yang & Kunz Reference Li, Yang and Kunz2022). Consequently, turbulence modelling is indispensable for predictive simulations. The overarching goal of turbulence modelling is to construct a mapping from the resolved scales to the unresolved ones. Conventional approaches achieve this by calibrating a model against a limited set of canonical flows. The underlying philosophy – or perhaps more accurately, the hope – is that once calibrated, the model will retain predictive capability for flows outside the calibration set, provided that these flows can in some sense be composed from those in the calibration set. This ability to extend predictive performance beyond the training domain to flows that share similar physics with the training data is what we refer to here as cross-scenario generalisation.
A familiar example is Reynolds-averaged Navier–Stokes (RANS) modelling, whose objective is to establish a mapping between the resolved mean flow and the unresolved turbulent stresses. Canonical models such as the one-equation Spalart–Allmaras (SA) model and the two-equation
$k$
–
$\omega$
shear stress transport model are calibrated against two-dimensional wakes and jets, as well as flat-plate boundary layers or the logarithmic law of the wall (LoW) (Spalart & Allmaras Reference Spalart and Allmaras1992; Menter Reference Menter1994). In practice, real-world flows often involve complex interactions among wakes, jets and boundary layers, and calibrations derived from canonical flows do not necessarily transfer directly. Nevertheless, the SA and
$k$
–
$\omega$
shear stress transport models are known to yield reasonably accurate predictions across a wide range of flows, thereby demonstrating cross-scenario generalisation.
The same philosophy underlies near-wall turbulence modelling in large-eddy simulations (LES). The aim of wall modelling is to relate the resolved turbulence at some distance from the wall to the unresolved wall stress. The widely used equilibrium wall model, for instance, is calibrated against the logarithmic LoW. While real flows may deviate from the ideal log-law behaviour, equilibrium wall models nonetheless provide robust predictions across a variety of configurations (Hanjalic Reference Hanjalic2005; Durbin Reference Durbin2018; Mani & Dorgan Reference Mani and Dorgan2023), again illustrating cross-scenario generalisation.
We note, however, that calibrations for conventional turbulence models – particularly those most widely deployed in production computational fluid dynamics (CFD) codes – are almost exclusively restricted to equilibrium flows. By conventional wisdom, one might expect that incorporating non-equilibrium or separated flows into the calibration dataset would improve model performance for such flows. This has long been recognised as a major challenge in turbulence modelling. Classical closures such as the
$k$
–
$\epsilon$
model were developed at a time when high-fidelity reference data for non-equilibrium flows were scarce. With advances in high-performance computing, however, direct numerical simulation (DNS) and wall-resolved LES (WRLES) now provide a wealth of high-quality data (Moin & Mahesh Reference Moin and Mahesh1998; Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008; Graham et al. Reference Graham2016). The bottleneck is no longer data availability but rather the limited descriptive power of traditional closures: with only a handful of adjustable parameters, these models struggle to represent the spatially evolving dynamics characteristic of non-equilibrium flows.
The advent of machine learning (ML) offers a potential solution. Modern tools such as artificial neural networks contain orders of magnitude more adjustable parameters than conventional turbulence models. This expressive capacity makes it possible, in principle, to fit complex data, including data from non-equilibrium and separated flows, at a level of fidelity unattainable with traditional closures.
The earliest efforts in data-driven turbulence modelling focused on a priori training, i.e. fitting high-fidelity data in an offline manner. These studies emphasised the incorporation of invariances – such as Galilean invariance and rotational invariance – to reduce the amount of data required for training (Ling, Kurzawski & Templeton Reference Ling, Kurzawski and Templeton2016; Wang, Wu & Xiao Reference Wang, Wu and Xiao2017; Wu, Xiao & Paterson Reference Wu, Xiao and Paterson2018). They also pursued model interpretability by employing white-box techniques such as symbolic regression (Weatheritt & Sandberg Reference Weatheritt and Sandberg2016). These considerations remain important, as underscored by more recent studies (Shan & Zhang Reference Shan and Zhang2025; Zhang, Xia & Kou Reference Zhang, Xia and Kou2025b ). Despite these advances, early offline-trained models often lacked robustness: in some cases, they produced unphysical results and even degraded accuracy for flows within the training dataset (Guo, Xia & Chen Reference Guo, Xia and Chen2021; Li et al. Reference Li, Bin, Huang and Yang2024). Recognising this limitation, Duraisamy (Reference Duraisamy2021) drew a distinction between a priori training and model-consistent training. He argued that ML-based corrections must remain compatible with the underlying turbulence framework, rather than being designed solely to fit DNS data offline. A number of frameworks have since embraced this philosophy (Zhao et al. Reference Zhao, Akolekar, Weatheritt, Michelassi and Sandberg2020; MacArt, Sirignano & Freund Reference MacArt, Sirignano and Freund2021; Bae & Koumoutsakos Reference Bae and Koumoutsakos2022; Fang et al. Reference Fang, Zhao, Waschkowski, Ooi and Sandberg2023; Vadrot et al. Reference Vadrot, Yang, Bae and Abkar2023b ; Zhou et al. Reference Zhou, Zhang, He and Yang2025). A prominent example is the field inversion and ML (FIML) framework, which also has the advantage of requiring only sparse reference data (Singh, Medida & Duraisamy Reference Singh, Medida and Duraisamy2017; Holland et al. Reference Holland, Baeder and Duraisamy2019b ). The procedure involves two steps. In the first step – field inversion – an optimisation problem is solved to spatially modify a component of the baseline model at each grid point so that the adjusted model matches the reference data. In the second step – machine learning – an ML algorithm is employed to relate user-selected flow features to the inverted model coefficient. Subsequent refinements sought to improve stability and generalisability, e.g. by bounding the input features (Srivastava & Duraisamy Reference Srivastava and Duraisamy2021) and by constraining the spatial region of influence of the learned corrections (Srivastava et al. Reference Srivastava, Rumsey, Coleman and Wang2024; Wu, Zhang & Zhang Reference Wu, Zhang and Zhang2025). These refinements proved fruitful: the resulting models demonstrated substantial improvements not only for the training flows but also for flows closely resembling the training dataset, thereby achieving what may be termed scenario-specific generalisation.
Nevertheless, the goal of turbulence modelling is cross-scenario generalisation. Criticism of ML-augmented turbulence models intensified when practitioners attempted to apply them outside the training domain and found that the augmentations could be destructive (Rumsey, Coleman & Wang Reference Rumsey, Coleman and Wang2022). In particular, some models failed to preserve fundamental calibrations, most notably the LoW (Chen et al. Reference Chen, Bin, Yang, Shi, Abkar and Park2023; Vadrot et al. Reference Vadrot, Yang and Abkar2023a ). These issues fuelled broader scepticism: Spalart (Reference Spalart2023) argued that ML has yet to succeed in turbulence modelling, and that stronger guidance is necessary.
Recent studies have attempted to explain this lack of cross-scenario generalisation. Bin et al. (Reference Bin, Chen, Huang and Yang2022) and Yang et al. (Reference Yang, Shan, Yang and Zhang2025) showed that the problem arises from forgetting, i.e. the failure to preserve calibrations already embedded in the baseline model. Follow-up work by Bin et al. (Reference Bin, Hu, Li, Grauer and Yang2024a , Reference Bin, Huang, Kunz and Yangb ) proposed remedies in the form of a rubber-band approach or progressive ML, in which constraints are imposed to ensure that data-driven corrections do not disrupt the LoW – identified as the most fundamental calibration for turbulence models (Spalart Reference Spalart2015). This perspective has since been adopted in several studies (Huang, Yang & Kunz Reference Huang, Yang and Kunz2019; Amarloo et al. Reference Amarloo, Rincón, Reclari and Abkar2023, Reference Amarloo, Zehtabiyan-Rezaie and Abkar2024; Rincón et al. Reference Rincón, Amarloo, Reclari, Yang and Abkar2023, Reference Rincón, Reclari, Yang and Abkar2025; Zhang et al. Reference Zhang, Zhou, Yang and He2025a ), achieving varying levels of cross-scenario generalisation. However, a notable limitation of the previous rubber-band approaches in Bin et al. (Reference Bin, Hu, Li, Grauer and Yang2024a , Reference Bin, Huang, Kunz and Yangb ), among others, is the treatment of the model constant. Many data-driven turbulence models learn local augmentations as functions of local flow features at each grid point. In contrast, existing progressive-ML approaches have largely been limited to updating model constants globally. As a result, it remains unclear whether progressive ML alone can resolve the cross-scenario generalisation issues documented in e.g. Rumsey et al. (Reference Rumsey, Coleman and Wang2022) and Wu et al. (Reference Wu, Zhang and Zhang2025).
Here, we address this gap by extending progressive ML to the FIML setting. Specifically, we combine the LoW-preserving ‘rubber-band’ SA formulation of Bin et al. (Reference Bin, Huang, Kunz and Yang2024b
) with the FIML paradigm of Parish & Duraisamy (Reference Parish and Duraisamy2016), enabling a spatially varying SA coefficient field to be learned and deployed. To isolate the effect of progressive ML, i.e. LoW preservation, we compare two strategies: a conventional, unconstrained FIML approach, and a constrained FIML approach that explicitly preserves the LoW. Training is restricted to periodic hills, while extrapolation tests span a diverse set of configurations, including channel flows at high Reynolds numbers, periodic hills with varying slopes, the backward-facing step, and the three-dimensional BeVERLI hill at a
$30^\circ$
angle of rotation. As conventional FIML has previously demonstrated only scenario-specific generalisation, this study highlights the important role of the LoW in cross-scenario generalisation. Here, cross-scenario generalisation is the transfer of learned flow physics from a training flow to a target flow that shares similar flow physics. Accordingly, training on periodic hills is not expected to make the model universally accurate for all non-equilibrium flows. Rather, the expectation is improved performance on flows that share physics with the periodic-hill training data (e.g. other periodic hills, the backward-facing step and the BeVERLI hill), while preserving baseline calibrations.
The remainder of this paper is organised as follows. Section 2 presents the methodology. Section 3 discusses the results. Finally, § 4 summarises the main findings and offers concluding remarks.
2. Methodology
This section outlines the methodology. The overall workflow is summarised in figure 1. The framework comprises four main components: field inversion, ML, deployment of the augmented model within a CFD solver, and assessment of model performance. Each component is described in detail in the following subsections. We begin by introducing the baseline one-equation SA model and the associated LoW constraint in § 2.1. The field inversion procedure based on physics-informed neural networks (PINNs) is presented in § 2.2, followed by the feedforward neural network (FNN) employed for ML in § 2.3. Details of the model deployment and validation strategy are given in § 2.4.
The CFD information used in this study falls into two categories. First, high-fidelity DNS/WRLES reference data are used for field inversion and for testing. These include the periodic-hill DNS database of Xiao et al. (Reference Xiao, Wu, Laizet and Duan2020), the channel-flow DNS of Lee & Moser (Reference Lee and Moser2015) at
$ \textit{Re}_{\tau }=5200$
, the canonical LoW scaling reported on the NASA TMR website (Rumsey, Vogel & Pederson Reference Rumsey, Vogel and Pederson2026) at
$ \textit{Re}_{\tau }\approx 10^6$
, and the BeVERLI WRLES data of Nair et al. (Reference Nair, Wadhai, Kunz and Yang2026). Second, the low-fidelity CFD data used for model deployment are RANS calculations generated in the present work. Their numerical fidelity is discussed in Appendix C, where we summarise the reported quality of the reference datasets, and document additional convergence evidence for representative low-fidelity RANS calculations.
Flow chart of the FIML framework. Training data consist of three periodic-hill cases, while testing is carried out on plane-channel flows, two-dimensional periodic hills, the two-dimensional backward-facing step, and the three-dimensional BeVERLI hill.
Figure 1 Long description
The figure summarises the field inversion and machine learning workflow. Training uses three periodic-hill cases. Testing covers plane channel flow, two-dimensional periodic hills, a two-dimensional backward-facing step, and the three-dimensional BeVERLI hill.
2.1. The SA model and the LoW constraint
The baseline model employed in this work is the one-equation SA turbulence model (Spalart & Allmaras Reference Spalart and Allmaras1992; Bin, Huang & Yang Reference Bin, Huang and Yang2023). The governing transport equation for the modified eddy viscosity
$\tilde {\nu }$
is
\begin{align} \frac {\mathrm{D}\tilde {\nu }}{\mathrm{D}t} = c_{b1}\tilde {S}\tilde {\nu } - c_{w1} f_w \left (\frac {\tilde {\nu }}{d}\right )^2 + \frac {1}{\sigma }\left [ \frac {\partial }{\partial x_{\kern-1pt j}}\left ((\nu +\tilde {\nu })\frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}}\right ) + c_{b2}\frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}}\frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}} \right ]\!, \end{align}
where
$\mathrm{D}/\mathrm{D}t$
denotes the material derivative. The turbulent eddy viscosity is defined as
$\nu _t=\tilde {\nu }f_{\nu 1}$
, and the Reynolds-stress tensor is defined here as
$R_{\textit{ij}}\equiv -\overline {u_i'u_{\kern-1pt j}'}$
, so that the anisotropic part under the Boussinesq hypothesis is
$R_{\textit{ij}}^d=2\nu _t S_{\textit{ij}}$
. The modified strain-rate magnitude is
$\tilde {S}=\varOmega +\tilde {\nu }f_{\nu 2}/(\kappa ^2 \text{d}^2)$
, where
$\varOmega$
is the vorticity magnitude,
$d$
is the distance to the nearest wall, and
$\kappa =0.41$
is the von Kármán constant. The auxiliary functions
$f_{\nu 1}$
and
$f_{\nu 2}$
take the forms
\begin{align} f_{\nu 1}=\left [1-\exp \left (-\frac {\chi }{\kappa A}\right )\right ]^2\quad (A=17),\quad f_{\nu 2}=1-\frac {\chi }{1+\chi f_{\nu 1}}, \end{align}
where
$\chi =\tilde {\nu }/\nu$
, and
$f_{\nu 1}$
is the damping function in Kawai & Larsson (Reference Kawai and Larsson2012) and Yang & Lv (Reference Yang and Lv2018).
The present form differs from the original
$f_{\nu 1}$
in Spalart & Allmaras (Reference Spalart and Allmaras1992), which does not yield the correct near-wall asymptotic scaling as
$y\to 0$
. The modified form in (2.2) ensures the proper scaling
$f_{\nu 1}\sim y^2$
while remaining otherwise consistent with the original formulation.
In this work, the model constant
$c_{b1}$
is treated as a local variable to be determined in a data-driven manner, following Singh et al. (Reference Singh, Medida and Duraisamy2017), Rumsey et al. (Reference Rumsey, Coleman and Wang2022) and Srivastava et al. (Reference Srivastava, Rumsey, Coleman and Wang2024).
The treatment of
$c_{w1}$
,
$\sigma$
and
$f_w$
when
$c_{b1}$
is varied distinguishes the constrained approach from the conventional unconstrained one.
In the unconstrained formulation, the remaining ingredients retain their baseline SA forms. Specifically,
$c_{w1}=3.28$
,
$c_{b2}=0.622$
,
$\sigma =2/3$
and
\begin{align} f_w=g{\left [\frac {1+2^6}{g^6+2^6}\right ]}^{1/6},\qquad g=r+0.3(r^6-r),\qquad r=\min {\left [\frac {\tilde {\nu }}{\tilde {S}\kappa ^2 \text{d}^2},\,10\right ]}. \end{align}
In contrast, the constrained approach varies
$c_{b1}$
,
$c_{w1}$
,
$\sigma$
, and
$f_w$
in a coupled manner so that calibration against the LoW is preserved. To derive the constraint that preserves the LoW in the logarithmic region, one invokes the two scalings
$\tilde {\nu }=\kappa d$
and
$\tilde {S}=1/{\kappa d}$
. These two scalings lead to
$r=1$
and
$f_w=1$
. Consequently, the
$\tilde {\nu }$
transport equation gives
\begin{align} c_{w1}=\frac {c_{b1}}{\kappa ^2}+\frac {1+c_{b2}}{\sigma }, \end{align}
with
$\sigma =2/3$
and
$c_{b2}=0.622$
. More details regarding the LoW constraint can be found in Bin et al. (Reference Bin, Huang, Kunz and Yang2024b
). In addition to (2.4), Bin et al. (Reference Bin, Huang, Kunz and Yang2024b
) expressed
$f_w$
as a function of
$c_{b1}$
and
$r$
, following equation (8) of Bin et al. (Reference Bin, Huang, Kunz and Yang2024b
):
\begin{align} \begin{split} &f_w(r)=\frac {c_{b1}}{\kappa ^2 c_{w1}}\frac {1}{r} +\frac {1+c_{b2}}{\sigma c_{w1}}\frac {F_1^2}{F_0^2} +\frac {1}{\sigma c_{w1}}\frac {F_2}{F_0^2}, \quad 0.4 \leqslant r \leqslant 1,\\[5pt]&f_w(r)=B\,\textrm {tanh}\left (\frac {r-1}{C}\right )+1,\quad r\gt 1,\\[5pt]&B=10^{2c_{s2}-1}-1,\quad C=10^{4c_{s1}-1}/5, \end{split} \end{align}
where
$c_{s1}=0.25$
and
$c_{s2}=0.65$
are model constants. For
$0\leqslant r\lt 0.4$
, the function is linearly continued from
$f_w(0)=0$
to the value at
$r=0.4$
, following Bin et al. (Reference Bin, Huang, Kunz and Yang2024b
). The quantities
$F_0$
,
$F_1$
and
$F_2$
are obtained from the channel-flow DNS data (Graham et al. Reference Graham2016):
\begin{align} F_0=\frac {\nu _t^{+, \textit{DNS}}}{d^+},\quad F_1=\frac {\text{d}\nu _t^{+, \textit{DNS}}}{\text{d}y^+},\quad F_2=\left (1+\nu _t^{+, \textit{DNS}}\right )\frac {\text{d}^2\nu _t^{+, \textit{DNS}}}{\text{d}y^{+2}}. \end{align}
Equations (2.4)–(2.6) are enforced during both the PINN-based field inversion and model deployment. Unlike approaches that rely on shielding, the present method applies the learned correction throughout the domain without region-specific treatment.
We summarise the two settings. In this work, ‘unconstrained’ refers to the conventional one-field FIML baseline, rather than to a fully unconstrained closure in a mathematical sense. Specifically, we invert only
$c_{b1}(\boldsymbol{x})$
while retaining the baseline SA functional forms for the remaining ingredients. This is consistent with common practice in prior FIML studies, where one augments a selected source term – e.g. via a learned
$\beta$
field multiplying production (Parish & Duraisamy Reference Parish and Duraisamy2016; Holland et al. Reference Holland, Baeder and Duraisamy2019a
,
Reference Holland, Baeder and Duraisamyb
; Rumsey et al. Reference Rumsey, Coleman and Wang2022; Srivastava et al. Reference Srivastava, Rumsey, Coleman and Wang2024) – while leaving the rest of the closure in its baseline form. This choice imposes an implicit closure assumption because once
$c_{b1}$
is modified, the unchanged remaining ingredients define a particular closure manifold. In the constrained framework, changes in
$c_{b1}$
induce changes in other model ingredients through the LoW constraint, but these induced changes are deterministically slaved to
$c_{b1}$
so that the number of independent degrees of freedom remains 1 in both frameworks. The comparison is therefore between one freely varying field and one LoW-constrained one-degree-of-freedom manifold, not between models with different numbers of independently tunable closure functions. This controlled comparison isolates the effect of LoW preservation on generalisation relative to conventional one-field FIML more cleanly than a comparison in which
$c_{w1}$
,
$\sigma$
and
$f_w$
are adjusted independently. Finally, the LoW constraint is introduced as an internal consistency requirement, not as a universal calibration for complex flows: channel flow is used as a canonical reference to construct a mapping that preserves the near-wall log-layer asymptote when
$c_{b1}$
is modified.
2.2. Field inversion using a PINN
The goal of field inversion is to identify a spatially varying distribution of
$c_{b1}$
such that the augmented SA model reproduces reference data with high fidelity. In this study, the reference data are the velocity fields over periodic hills with slopes
$\alpha =0.5$
,
$0.8$
and
$1.2$
, taken from the DNS database of Xiao et al. (Reference Xiao, Wu, Laizet and Duan2020). The periodic-hill configuration is a widely used benchmark in turbulence modelling, and details of the geometry can be found in Xiao et al. (Reference Xiao, Wu, Laizet and Duan2020) and are not repeated here for brevity.
Schematic of the FIML framework. A PINN is employed for field inversion, and an FNN is used for ML. The SA model serves as the baseline, with ML augmentation applied to
$c_{b1}$
.
Figure 2 Long description
The schematic shows how reference flow data are used in a physics-informed neural network for field inversion. The inferred spatial coefficient field is then learned by a feedforward neural network and embedded into the Spalart-Allmaras model as an augmentation to c_b1.
We employ a PINN (Raissi et al. Reference Raissi, Perdikaris and Karniadakis2019, Reference Raissi, Yazdani and Karniadakis2020; Cai et al. Reference Cai, Mao, Wang, Yin and Karniadakis2021; Molnar & Grauer Reference Molnar and Grauer2022; Molnar et al. Reference Molnar, Venkatakrishnan, Schmidt, Sipkens and Grauer2023; Song et al. Reference Song, Cao, Liao and Zhang2025) for field inversion. Figure 2(a) illustrates the PINN set-up. The network architecture consists of eight hidden layers, each with 60 neurons, following Patel et al. (Reference Patel, Mons, Marquet and Rigas2024). The inputs are the spatial coordinates
$(x,y)$
, while the outputs include the velocity components
$u$
and
$v$
, the pressure
$p$
, the logarithm of the modified eddy viscosity field
$\log (\tilde {\nu })$
, and the spatially varying
$c_{b1}$
. The logarithm of
$\tilde {\nu }$
is used as an output rather than
$\tilde {\nu }$
because using
$\tilde {\nu }$
directly was found to cause convergence difficulties. All variables are non-dimensionalised using the hill height (
$H$
) and bulk flow rate (
$U_b$
). The loss function comprises residuals of the governing equations, residuals from the boundary conditions, discrepancies with reference data, and regularisation. Here, regularisation is introduced to mitigate non-uniqueness and to promote a stable solution. For the present problem, the physical losses correspond to the RANS equations in dimensionless form: the continuity equation, the momentum equations and the transport equation for
$\tilde {\nu }$
:
\begin{align} L_{\textit{mass}} = \frac {\partial U_{\kern-1pt j}}{\partial x_{\kern-1pt j}}, \\[-28pt] \nonumber \end{align}
\begin{align} L_{i,\textit{momentum}} = U_{\kern-1pt j}\frac {\partial U_i}{\partial x_{\kern-1pt j}} + \frac {\partial p}{\partial x_i} - \frac {\partial }{\partial x_{\kern-1pt j}}\left [\left (\nu +\nu _t\right ) \left (\frac {\partial U_i}{\partial x_{\kern-1pt j} }+\frac {\partial U_{\kern-1pt j}}{\partial x_i }\right )\right ] \!, \\[-28pt] \nonumber \end{align}
\begin{align} L_{\tilde {\nu }} = \underbrace {U_{\kern-1pt j} \frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}}}_{\textit{transport}} - \underbrace {c_{b1} \tilde {S} \tilde {\nu } \vphantom {\left (\frac {\tilde {\nu }}{d}\right )^2}}_{\textit{production} } + \underbrace {c_{w1} f_w \left (\frac {\tilde {\nu }}{d}\right )^2}_{\textit{destruction}} - \underbrace {\frac {1}{\sigma } \left [ \frac {\partial }{\partial x_{\kern-1pt j}} \left ((\nu + \tilde {\nu }) \frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}}\right ) + c_{b2} \frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}} \frac {\partial \tilde {\nu }}{\partial x_{\kern-1pt j}} \right ]}_{\textit{diffusion}}. \\[-12pt] \nonumber \end{align}
Here, the subscript
$i$
equals 1 or
$2$
, and summation over repeated indices is implied. The boundary conditions at the upper and lower walls impose no-slip conditions for velocity, a fixed value of zero for
$\tilde {\nu }$
, and zero gradient for pressure. Periodic boundary conditions are applied to the inlet and outlet, leading to
\begin{align} \begin{aligned} &L_{U_i,\textit{no slip}} = U_{i,w\textit{all}}, \quad L_{\tilde {\nu },w\textit{all}} = \tilde {\nu }_{w\textit{all}}, \quad L_{p,w\textit{all}} = \left .\frac {\partial p}{\partial x_n}\right |_{w\textit{all}}, \\[5pt]&L_{U_i,\textit{periodic}} = U_{i,\textit{inlet}} - U_{i,\textit{outlet}}, \quad L_{\tilde {\nu },\textit{periodic}} = \tilde {\nu }_{\textit{inlet}} - \tilde {\nu }_{\textit{outlet}}, \\[5pt]& L_{p,\textit{periodic}} = p_{\textit{inlet}} - p_{\textit{outlet}}. \end{aligned} \end{align}
The data loss penalises deviations between the PINN-predicted velocity field and the DNS reference data:
\begin{align} L_{i,\textit{data}} = U_i - U_{i,\textit{DNS}}. \end{align}
A penalty term is also included to limit deviations of
$c_{b1}$
from its baseline value:
\begin{align} L_{c_{b1}} = c_{b1} - 0.1355. \end{align}
The total loss is defined as the mean squared error (MSE) of all contributions:
\begin{align} L & = \frac {1}{\mathcal{V}} \int _V \left ( \sum _{i=1}^2L_{i,\textit{data}}^2 + L_{\textit{mass}}^2 + \sum _{i=1}^2L_{i,\textit{momentum}}^2 + L_{\tilde {\nu }}^2 +\lambda L_{{cb_1}}^2\right ) \text{d}V' \nonumber\\[4pt]& \quad + \frac {1}{\mathcal{O}} \int _O \left (\sum _{i=1}^2L_{u_i,\textit{no slip}}^2 +L^2_{\tilde {\nu },w\textit{all}}+ L_{p,w\textit{all}}^2+\sum _{i=1}^2L^2_{u_i,\textit{periodic}} \right. \nonumber\\[4pt]& \quad \left. +L^2_{\tilde {\nu },\textit{periodic}}+L_{p,\textit{periodic}}^2 \,\right ) \text{d}l'. \end{align}
Here, the regularisation weight is
$\lambda =0.0001$
following Holland et al. (Reference Holland, Baeder and Duraisamy2019b
), the volume integral
$\int _V$
covers the entire flow domain, and the surface integral
$\int _O$
covers all boundaries;
$V'$
and
$l'$
are dummy variables, and
$\mathcal{V}$
and
$\mathcal{O}$
denote the areas of the fluid and boundary surfaces, respectively.
2.3. Machine learning via FNN
The objective of the ML step in the FIML framework is to learn a functional relationship between selected local flow features and the spatial distribution of
$c_{b1}$
obtained from field inversion. Following Singh et al. (Reference Singh, Medida and Duraisamy2017), Rumsey et al. (Reference Rumsey, Coleman and Wang2022) and Srivastava et al. (Reference Srivastava, Rumsey, Coleman and Wang2024), we employ an FNN for this purpose. A schematic of the architecture is shown in figure 2(b). The FNN consists of two hidden layers, each with 40 neurons, and uses the hyperbolic tangent (tanh) activation function. The target output is the inverted
$c_{b1}$
field. The inputs are chosen to represent key aspects of the local flow state: the ratio of the production term
$P$
to the destruction term
$D$
with
$P=c_{b1}\tilde {S}\tilde {\nu }$
and
$D=c_{w1}f_w(\tilde {\nu }/d)^2$
, the viscosity ratio
$\chi =\tilde {\nu }/\nu$
, the near-wall damping function
$f_w$
, and the rotational Reynolds number
$\varOmega d_{w}^2/\nu$
, where
$d_w$
denotes the wall distance. These features are physically motivated:
$|P|/|D|$
characterises the local balance between modelled production and destruction;
$\chi$
measures the turbulence level relative to molecular diffusion;
$f_w$
reflects departures from equilibrium in the near-wall model behaviour; and
$\varOmega d_w^2/\nu$
has been used as an indicator of flow separation. Further discussion regarding the discriminative power, the non-dimensional consistency, and invariance considerations of these input features can be found in Srivastava & Duraisamy (Reference Srivastava and Duraisamy2021) and Srivastava et al. (Reference Srivastava, Rumsey, Coleman and Wang2024). Input feature engineering remains an active research area. While alternative feature sets may be more effective for broader classes of flows, the focus of the present work is to assess whether enforcing the LoW constraint enables improved generalisation. In the constrained formulation,
$c_{b1}$
depends on
$f_w$
, while
$f_w$
depends on
$c_{b1}$
through the LoW constraint. Such coupling is common in RANS modelling and is resolved by the iterative solution procedure. For example, in the baseline SA model, the
$f_w$
function appears in the
$\tilde {\nu }$
transport equation, yet it depends on
$\chi =\tilde {\nu }/\nu$
. The use of model-internal quantities such as production–destruction ratios and damping functions as ML inputs also follows prior FIML practice (Parish & Duraisamy Reference Parish and Duraisamy2016; Singh et al. Reference Singh, Medida and Duraisamy2017; Holland et al. Reference Holland, Baeder and Duraisamy2019b
; Rumsey et al. Reference Rumsey, Coleman and Wang2022; Srivastava et al. Reference Srivastava, Rumsey, Coleman and Wang2024). In the present tests, we have not observed numerical inconsistency or performance degradation attributable to this coupled input–output structure when the model is deployed self-consistently. Nevertheless, because
$f_w$
and
$c_{b1}$
are algebraically coupled in the constrained formulation, the possibility that the network partly learns this imposed relation should be examined further in future work. To reduce extrapolation in feature space, we adopt the bounded re-parametrisations proposed in Srivastava & Duraisamy (Reference Srivastava and Duraisamy2021) and Srivastava et al. (Reference Srivastava, Rumsey, Coleman and Wang2024):
$|P|/(|D|+|P|)$
,
$100/(\chi +100)$
,
$1/(f_w+1)$
and
$1000\nu /(d_w^2\varOmega +1000\nu )$
. The network is implemented in the PyTorch library (Paszke et al. Reference Paszke2019) and trained using the Adam optimiser. The resulting model achieves a coefficient of determination
$R^2 \approx 0.98$
.
The standard two-step FIML workflow may exhibit FI–ML mismatch: regression errors can be amplified by the nonlinear RANS solver, and the feature distribution may shift because features are recomputed from the iterated solution during deployment. Here, the learned model is deployed self-consistently, and we mitigate feature-space extrapolation using bounded feature re-parametrisations (Srivastava & Duraisamy Reference Srivastava and Duraisamy2021; Srivastava et al. Reference Srivastava, Rumsey, Coleman and Wang2024). An alternative is integrated (end-to-end) training through the coupled solver (Holland et al. Reference Holland, Baeder and Duraisamy2019a , Reference Holland, Baeder and Duraisamyb ; Srivastava & Duraisamy Reference Srivastava and Duraisamy2021; Gupta & Duraisamy Reference Gupta and Duraisamy2026) at substantially higher computational cost; this alternative is not pursued here, since our main focus concerns LoW preservation rather than the choice of two-step versus integrated training.
2.4. Model deployment and testing
Once training is complete, the neural network is embedded into a RANS solver for offline testing. The solver used in this work is OpenFOAM (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998). The implementation is model-consistent: the FNN is integrated directly into the solver, with its inputs computed at each iteration, passed through the network to generate an updated
$c_{b1}$
field, then fed back into the RANS equations following prior work (Singh et al. Reference Singh, Medida and Duraisamy2017; Holland et al. Reference Holland, Baeder and Duraisamy2019b
; Rumsey et al. Reference Rumsey, Coleman and Wang2022; Srivastava et al. Reference Srivastava, Rumsey, Coleman and Wang2024). Additional information on the reference datasets and on the numerical fidelity of the low-fidelity RANS calculations is collected in Appendix C.
Flow configurations. (a) The three-dimensional BeVERLI hill flow configuration. The size of the domain is
$L_x\times L_y\times L_z=36.13H\times 9.9H\times 4H$
, with distance
$6.73H$
from the precursor inlet to the leading edge of the hill. The domain has a slip top boundary, a no-slip bottom wall, and a periodic boundary condition in the transverse direction. (b) An illustration of the BeVERLI hill geometry. The width of the hill is
$W=5H$
. The hill has a flat top and super-elliptic corners.
Figure 3 Long description
The figure has two schematic panels. One shows the computational domain around the BeVERLI hill, including inlet distance, domain dimensions, wall and top boundary conditions, and periodic spanwise direction. The other shows the BeVERLI hill geometry, including its width, flat top, and rounded super-elliptic corners.
We use periodic-hill cases with
$\alpha =0.5$
,
$0.8$
and
$1.2$
for training. The ML-augmented model is then assessed on a series of canonical flows. We first present the results for periodic-hill cases at
$\alpha =1.0$
,
$1.2$
and
$1.5$
. The
$\alpha =1.2$
case is part of the training dataset, whereas the
$\alpha =1.0$
and
$\alpha =1.5$
cases correspond to interpolation and extrapolation relative to the training set, respectively. Next, we consider two plane-channel validations. These configurations lie entirely outside the training dataset, and feature geometries that are distinct from those in training. Nevertheless, they belong to the calibration set of the baseline model. The
$ \textit{Re}_{\tau }=5200$
case is benchmarked against the DNS of Lee & Moser (Reference Lee and Moser2015), whereas the
$ \textit{Re}_{\tau }\approx 10^6$
case is benchmarked against the canonical reference scaling reported on the NASA TMR website (Rumsey et al. Reference Rumsey, Vogel and Pederson2026). Finally, we assess performance on more challenging separated flows: the backward-facing step at
$ \textit{Re}_h=5600$
, and the three-dimensional BeVERLI hill at hill-height Reynolds number
$ \textit{Re}_H=15{\,}000$
. The geometry of the backward-facing step is standard and is not elaborated here for brevity, with reference data taken from Bin et al. (Reference Bin, Park, Lv and Yang2024c
). The BeVERLI hill geometry is shown in figure 3, with further details available in Lowe et al. (Reference Lowe, Roy, Devenport, Borgoltz, Grzyb, Shanmugam, Borole and Gargiulo2024) and Roy et al. (Reference Roy, Lowe, Devenport, Borgoltz, Grzyb, Borole, Shanmugam, HogePatil and Gargiulo2024); the high-fidelity reference data are provided by the WRLES of Nair et al. (Reference Nair, Wadhai, Kunz and Yang2026). Both cases are absent from the training dataset and differ markedly in geometry from the training configurations. Moreover, the baseline SA model performs poorly for these flows, particularly in predicting separation-bubble size and downstream skin-friction distribution over the three-dimensional hill. These two test cases therefore constitute the most stringent assessments of the generalisability of the ML-augmented model.
The overarching objective is to evaluate whether training on one separated flow – absent from the baseline model’s calibration dataset – can enhance predictive capability across a broader class of separated-flow configurations.
3. Results
In this section, we evaluate the performance of the ML-augmented model across the validation cases introduced in § 2.4. The results are presented in increasing order of complexity, beginning with periodic-hill flows and followed by channel flow. We then turn to more challenging separated flows, namely the backward-facing step and the three-dimensional BeVERLI hill. For each case, predictions of the present LoW-constrained FIML model and the conventional unconstrained FIML model are compared against one another, as well as against DNS, WRLES, experimental data or canonical reference scaling, as appropriate. The comparisons focus on mean velocity profiles, skin-friction distributions, and separation characteristics. Particular attention is given to assessing the extent to which training on a single separated flow enhances predictive capability across a broader class of separated-flow configurations.
3.1. Periodic hill with
$\alpha = 1.2$
Streamwise velocity profiles for the
$\alpha =1.2$
periodic-hill case. Here, ‘Uncons.’ and ‘Cons.’ refer to unconstrained and constrained, respectively.
Figure 4 Long description
The figure compares streamwise velocity profiles at several streamwise stations for DNS, the baseline Spalart-Allmaras model, unconstrained FIML, and constrained FIML. The baseline model differs noticeably from DNS in the near-wall and separated-flow regions. Both FIML variants improve the velocity prediction, with the unconstrained FIML showing the largest reduction in error for this training-set case.
We begin with the periodic-hill case at
$\alpha = 1.2$
. This case is included in the training dataset. Figure 4 presents the streamwise velocity profiles at multiple streamwise locations. The baseline SA model shows noticeable deviations from the DNS results, particularly in the near-wall region. The constrained and unconstrained FIML models yield slightly and substantially improved velocity predictions, respectively. To quantify the prediction error relative to the DNS reference data, we define the MSE as
\begin{align} \textit{MSE} =\frac {1}{V} \int _V \left |\boldsymbol{U}_{\textit{mod}} - \boldsymbol{U}_{\textit{ref}}\right |^2/U_b^2 \, \text{d}V', \end{align}
where
$\boldsymbol{U}$
is the velocity vector, and the subscripts ‘mod’ and ‘ref’ denote the model prediction and reference data, respectively. The MSE values for the baseline SA model, the unconstrained FIML model and the constrained FIML model are
$0.0119$
,
$0.0031$
and
$0.0071$
, respectively. The constrained FIML model shows less improvement than the unconstrained FIML model. This is expected for flows that are well represented in the training dataset, since the unconstrained FIML model has greater descriptive power than its constrained counterpart. Because this case is part of the training set, any remaining discrepancies can be attributed to the limited discriminative power of the input features and the limited expressiveness of the neural network (Srivastava Reference Srivastava2022; Gupta & Duraisamy Reference Gupta and Duraisamy2026).
Comparison of
$c_{b1}$
and
$\nu _t$
for the periodic-hill case at
$\alpha = 1.2$
obtained using (a,d) the baseline SA model, (b,e) the unconstrained FIML model, and (c,f) the constrained FIML model. The dashed lines in (b,c) highlight approximately the baseline value of
$c_{b1}$
.
Figure 5 Long description
The multi-panel contour figure compares baseline, unconstrained FIML, and constrained FIML predictions for the spatially varying production coefficient c_b1 and eddy viscosity nu_t. The unconstrained model enhances c_b1 mainly in the shear layer from the hill crest, whereas the constrained model modifies c_b1 more around and outside that shear-layer region. The eddy-viscosity fields remain qualitatively similar.
Before proceeding, we examine the
$c_{b1}$
and
$\nu _t$
fields in figure 5. Both the constrained and unconstrained FIML models produce enhanced eddy viscosity upstream of the hill, with the constrained FIML model predicting a larger increase. Although the resulting eddy-viscosity fields are similar, the underlying modifications that produce them differ between the constrained and unconstrained FIML models. Specifically, the unconstrained FIML model enhances
$c_{b1}$
within the shear layer emanating from the hill crest, whereas the constrained FIML model modifies
$c_{b1}$
primarily outside the shear-layer region. It is also worth noting that
$c_{b1}$
is not necessarily positive. While a globally negative
$c_{b1}$
would imply negative turbulence production and lead to numerical issues, locally negative values of
$c_{b1}$
are not problematic. Indeed, such values commonly arise when clipping is not enforced (Yan et al. Reference Yan, Li, Zhang and Chen2022). From a physical standpoint, negative turbulence production has also been reported in DNS studies of separated flows (Abe et al. Reference Abe, Mizobuchi, Matsuo and Spalart2012; Balakumar & Park Reference Balakumar and Park2015; Cimarelli et al. Reference Cimarelli, Leonforte, De Angelis, Crivellini and Angeli2019). Despite this difference, both models appear to identify similar flow structures, as evidenced by the dashed contour lines in figures 5(b,c), which indicate approximately the baseline
$c_{b1}$
values.
3.2. Periodic hills with
$\alpha =1.0$
and 1.5
We next evaluate the performance of the ML-augmented models for periodic hills with
$\alpha =1.0$
and
$1.5$
. These cases represent interpolation and extrapolation relative to the training dataset, respectively, while retaining geometric similarity to the training flows.
Figure 6 presents the streamwise velocity results. The trends are similar to those observed for the
$\alpha = 1.2$
case. Both the constrained and unconstrained FIML models improve upon the baseline SA predictions, with the unconstrained FIML model providing the greater improvement. Quantitatively, at
$\alpha = 1.0$
, the baseline SA, unconstrained FIML and constrained FIML models yield MSE values
$0.0080$
,
$0.0026$
and
$0.0045$
, respectively. At
$\alpha = 1.5$
, the corresponding MSE values are
$0.0179$
,
$0.0047$
and
$0.0112$
, respectively.
Streamwise velocity profiles of the periodic hill for slopes (a)
$\alpha =1.0$
and (b)
$\alpha =1.5$
.
Figure 6 Long description
The figure compares velocity profiles for two periodic-hill geometries that are not identical to the training cases. For both slopes, the baseline model shows larger deviations from the DNS reference than the FIML models. Both constrained and unconstrained FIML improve the profiles, with the unconstrained model giving the lower velocity mean-squared error for these geometrically similar cases.
3.3. Channel flows
We next consider plane-channel flows. These cases are outside the training dataset but are part of the calibration set of the baseline model. Figure 7(a) shows velocity profiles at
$ \textit{Re}_\tau =5200$
, benchmarked against the DNS of Lee & Moser (Reference Lee and Moser2015). All three models are fairly accurate, although the unconstrained FIML exhibits slight deterioration compared with the baseline SA and the constrained FIML in the neighbourhood of
$y^+\approx O(100)$
. This discrepancy becomes much more evident only in the very-high-Reynolds-number case, shown in figure 7(b) for
$ \textit{Re}_\tau \approx 10^6$
and compared against the canonical reference scaling. Here, the constrained FIML reproduces the log LoW by construction, whereas the unconstrained FIML produces an incorrect log-law slope. This outcome is not unexpected. The ML models, if not otherwise constrained, are known to suffer from catastrophic forgetting, whereby previously learned behaviour is lost when the model is retrained. In this case, the unconstrained FIML fails to retain the LoW. Appendix D extends this comparison by including channel-flow results at several Reynolds numbers between
$ \textit{Re}_{\tau }=5200$
and
$ \textit{Re}_\tau=10^6$
, and showing that the centreline-velocity deviation from the canonical LoW grows gradually with Reynolds number for the unconstrained model.
Velocity profiles of plane-channel flow at (a)
$ \textit{Re}_\tau =5200$
and (b)
$ \textit{Re}_\tau \approx 10^6$
. These validations are outside the training dataset but are part of the calibration set of the original SA model. The
$ \textit{Re}_\tau =5200$
case is compared against DNS, whereas the
$ \textit{Re}_\tau \approx 10^6$
case is compared against canonical reference scaling.
Figure 7 Long description
The two panels compare channel-flow velocity profiles predicted by the baseline model, unconstrained FIML, constrained FIML, and reference data or reference scaling. At Re_tau 5200 all models are close to the reference. At Re_tau about 10^6, the constrained FIML and baseline model preserve the logarithmic law-of-the-wall slope, while the unconstrained FIML produces an incorrect slope.
3.4. Backward-facing step
We next examine the backward-facing step. Figure 8 shows the skin-friction coefficient
$C_{\kern-1pt f}$
, comparing the baseline SA model, the constrained and unconstrained FIML models, and DNS reference data. The baseline SA model overpredicts the separation-bubble size. The unconstrained FIML degrades performance, producing excessively strong reverse flow immediately downstream of the step, and predicting reattachment too early. The constrained FIML, by contrast, offers moderate improvements: while it does not eliminate the baseline errors in the near-step region, it avoids further degradation and predicts the reattachment location with high accuracy. For a more quantitative assessment, we define the percentage error in the predicted separation-bubble size as
\begin{align} \textit{Err}_{\textit{BFS}} = \frac {|x_{\textit{reattachment,mod}} - x_{\textit{reattachment,ref}}|}{x_{\textit{reattachment,ref}}}. \end{align}
According to this metric, the baseline SA, unconstrained FIML and constrained FIML models yield errors
$15.67\,\%$
,
$21.41\,\%$
and
$1.57\,\%$
, respectively. The degraded performance of the unconstrained FIML model is consistent with a well-known limitation of ML-based turbulence models, namely their poor extrapolation outside the training range. The constrained FIML model alleviates this generalisation issue. In the present case, it more accurately captures the boundary layer recovery following separation, which leads to a highly accurate prediction of the separation-bubble size. We emphasise, however, that this level of accuracy is likely fortuitous and should not be expected for other types of separated flows. The main point here is that the constrained FIML model does not degrade the performance of the baseline model, and instead provides a slight improvement even when applied outside the training dataset.
Skin-friction coefficient
$C_{\kern-1pt f}$
for the backward-facing step.
Figure 8 Long description
The plot compares skin-friction coefficient along the wall for DNS, the baseline Spalart-Allmaras model, unconstrained FIML, and constrained FIML. The baseline model overpredicts the separation-bubble size. The unconstrained model degrades the prediction and reattaches too early, while the constrained model avoids that degradation and predicts a reattachment location close to DNS.
Comparison of the predicted
$c_{b1}$
fields for the backward-facing step: (a) unconstrained FIML and (b) constrained FIML. The dashed lines indicate approximately the baseline values of
$c_{b1}$
.
Figure 9 Long description
The two contour panels show spatial distributions of the learned c_b1 correction for backward-facing-step flow. The unconstrained and constrained models identify broadly similar separated-flow structures, marked by dashed contours near the baseline c_b1 value, but the unconstrained model modifies the shear layer more directly while the constrained model distributes changes around the shear-layer region.
Before proceeding, we examine the predicted
$c_{b1}$
fields. Figure 9 shows the spatial distributions of the learned corrections from the unconstrained and constrained FIML models. The results are similar to those in figure 5: the unconstrained FIML model modifies the eddy viscosity within the shear layer, whereas the constrained FIML model produces changes around the shear layer. At the same time, both models identify the same flow structures, as indicated by the dashed contour lines in the figure. In addition, both the constrained and unconstrained FIML models predict enhanced
$c_{b1}$
slightly above the top wall.
3.5. The BeVERLI hill
Finally, we present the results for the BeVERLI hill. Figure 10 presents the skin-friction distributions, shown along
$z=0$
(the centreline) and
$z=-H/4$
. Due to the three-dimensional nature of the flow, we define the skin-friction coefficient using the magnitude of the wall shear stress following Sundarraj et al. (Reference Sundarraj2023), Lowe et al. (Reference Lowe, Roy, Devenport, Borgoltz, Grzyb, Shanmugam, Borole and Gargiulo2024) and Roy et al. (Reference Roy, Lowe, Devenport, Borgoltz, Grzyb, Borole, Shanmugam, HogePatil and Gargiulo2024). All models capture the skin-friction drop near the hill crest (in the neighbourhood of
$x=0$
), although the magnitude of the drop is underpredicted relative to the high-fidelity WRLES reference. The primary region of interest is
$x/H\gtrsim 1$
, where the flow separates. The baseline SA model reproduces the separation on the lee side of the hill and the associated skin-friction drop at
$x/H \approx 1.4$
in both planes, but it predicts an unphysical rapid recovery, and fails to capture the secondary drop of the skin friction at
$x/H \approx 3$
. The unconstrained FIML predictions are qualitatively similar to those of the baseline model, but instead of predicting an overly fast recovery, they predict an overly slow recovery. In contrast, the constrained FIML significantly improves agreement with WRLES, accurately predicting both the location and magnitude of the secondary drop in the skin friction. For a more quantitative measure of the error, we define
\begin{align} {\textit{ERR}_{\textit{bhill}} =\frac { \int _{-H}^{4H} \left | C_{f,\textit{model}} - C_{f,\textit{WRLES}} \right | \, \text{d}x } {\int _{-H}^{4H} \left | C_{f,\textit{WRLES}} \right | \, \text{d}x}.} \end{align}
The integration range
$-1\leqslant x/H \leqslant 4$
covers a region that is slightly wider than the flow-separation region. At
$z=0$
, the baseline SA, unconstrained FIML and constrained FIML yield errors
$34.52\,\%$
,
$40.74\,\%$
and
$27.10\,\%$
, respectively. At
$z=-H/4$
, the baseline SA, unconstrained FIML and constrained FIML yield errors
$45.68\,\%$
,
$58.26\,\%$
and
$30.75\,\%$
, respectively.
Skin-friction coefficient for the BeVERLI hill at rotation angle
$30^\circ$
along (a)
$z=0$
(centreline) and (b)
$z=-H/4$
.
Figure 10 Long description
The two panels compare skin-friction coefficient along the BeVERLI hill at the centreline z = 0 and at z = -H/4. WRLES, baseline Spalart-Allmaras, unconstrained FIML, and constrained FIML are shown. The baseline and unconstrained models miss part of the downstream recovery behaviour, whereas the constrained model better captures the secondary skin-friction drop and reduces the integrated error.
Figure 11 shows contours of the skin-friction coefficient and representative streamlines. The baseline SA model and the unconstrained FIML model both predict an extended low-shear-stress region on the leeward surface of the hill, which is consistent with the recovery errors seen in figure 10. By contrast, the constrained FIML model predicts a wall shear distribution that is much more consistent with the WRLES reference, including a more faithful lee-side skin-friction drop and downstream recovery. This result is expected because the periodic-hill training cases and the BeVERLI test case share certain qualitative flow physics. In the present paper, ‘shared physics’ is meant qualitatively rather than in the sense of exact geometric similarity. The periodic hills, backward-facing step and BeVERLI hill all exhibit adverse-pressure-gradient separation, lee-side shear-layer development, and downstream recovery, even though their geometries differ substantially. This interpretation is consistent with the BeVERLI results, and together with those in § 3.4, supports the view that preserving the LoW helps learned behaviour to transfer more robustly between separated flows with related separation and recovery physics.
Contours of skin-friction coefficient and representative streamlines for the BeVERLI hill: (a) baseline SA, (b) unconstrained FIML, (c) constrained FIML, and (d) reference WRLES.
Figure 11 Long description
The four panels compare baseline Spalart-Allmaras, unconstrained FIML, constrained FIML, and WRLES reference results on the BeVERLI hill. Each panel shows skin-friction contours with representative streamlines. The constrained FIML panel most closely resembles the WRLES reference, including the lee-side low-shear region and downstream recovery.
3.6. Summary
We summarise the results in table 1. Three categories of test cases are considered: first, flows that are closely represented in the training dataset, namely the periodic-hill cases; second, a flow that is part of the calibration set of the baseline model, namely channel flow; and third, flows outside the training dataset with different geometries. The unconstrained FIML model provides substantial improvements for flows similar to those in the training dataset, but degrades the baseline model outside this range, even violating the LoW, which is part of the calibration of the baseline model. The constrained FIML model addresses this lack of generalisation. Although it achieves smaller improvements for flows within the training dataset, it provides more generalisable improvements for the backward-facing step and the BeVERLI hill, while preserving the calibration of the baseline model.
Summary of the test results. Here, Uncons. and Cons. denote unconstrained and constrained, respectively.
Table 1 Long description
The table lists eight test cases, their relationship to the training set, the error metric, and results for the baseline Spalart-Allmaras model, unconstrained FIML, and constrained FIML. For periodic hills, unconstrained FIML gives the lowest velocity MSE. For channel flow at Re_tau about 10^6, the unconstrained model has an incorrect log-law slope of 0.34, while the baseline and constrained models retain 0.41. For the backward-facing step and BeVERLI hill, constrained FIML gives the lowest reported errors.
4. Conclusions
This study examined the role of preserving the law of the wall (LoW) in cross-scenario generalisation of data-driven turbulence closures. Within the field inversion and machine learning (FIML) framework applied to the Spalart–Allmaras model, we compared a conventional, unconstrained FIML approach with a constrained formulation that enforces LoW preservation. Training was performed on periodic-hill flows, and validation included plane-channel flow, a backward-facing step, and the three-dimensional BeVERLI hill.
The results show that LoW preservation is largely immaterial when extrapolation remains close to the training distribution, as in the periodic-hill cases. In contrast, when applied outside the training distribution, the benefit of enforcing the constraint becomes clear. In the very-high-
$ \textit{Re}$
channel case, unconstrained FIML violates logarithmic scaling in the overlap region, and yields incorrect wall shear stress, whereas the constrained formulation retains the correct near-wall behaviour by construction. In separated flows such as the backward-facing step and BeVERLI hill, unconstrained FIML fails to improve upon the baseline model and can even degrade it, whereas the constrained formulation yields consistent gains, including more accurate reattachment predictions. Importantly, these gains are not due to recalibrating the LoW, since the baseline model already contains such a calibration. Rather, LoW preservation acts as a guardrail that prevents learned corrections from corrupting baseline near-wall equilibrium behaviour, enabling physics learned from the training data to transfer as a net improvement outside the training distribution.
Overall, these comparisons provide evidence that LoW preservation is an important ingredient for achieving cross-scenario generalisation in ML-augmented turbulence modelling. The constraint helps to prevent learned modifications from erasing fundamental baseline calibrations, and thereby improves extrapolation relative to unconstrained FIML while maintaining baseline performance on canonical attached flows. This is a step towards data-driven closures that remain reliable when deployed beyond the training distribution.
We conclude by noting several limitations and directions for future work. The evidence presented here is based on a specific FIML setting and a finite set of training and validation flows. Preserving the LoW is not claimed to universally improve all non-equilibrium or three-dimensional separated flows. Its role is to prevent degradation of baseline near-wall equilibrium behaviour and thereby improve extrapolation relative to unconstrained FIML. An important open question is when an equilibrium-flow constraint acts as a beneficial guardrail, and when it becomes restrictive for non-equilibrium or separated-flow adaptation. In the present comparison, the constrained and unconstrained FIML models have the same number of independently learned fields, so the LoW constraint does not reduce the nominal number of tunable functions. A more material loss of adaptability may arise in higher-dimensional augmented closures, where preserving coefficient relations necessarily restricts the admissible combinations of independently tunable model ingredients. Input-feature design remains an active research area, and can constrain performance even on training cases. Finally, while the PINN-based inversion is validated here against an adjoint baseline for representative periodic-hill cases, broader validation and comparisons across inversion strategies are worthwhile. Future work will explore improved feature sets, diagnostics for coupled input–output relations, additional physically motivated constraints beyond the LoW, and extensions to other turbulence closures and higher-Reynolds-number configurations.
Funding
X.I.A.Y. acknowledges AFOSR grant no. FA9550-23-1-0272, with Dr G. Abate as the technical monitor. R.F.K. acknowledges ONR grant N00014-24-1-2170, with P. Chang and J. Young as technical monitors. Computational resources were provided through Penn State ROAR and a United States Department of Defense Frontier project of the High Performance Computing Modernization Program.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request. The trained neural network for the constrained FIML model used in this study is available at https://github.com/Constrained-FIML/constrained_NN.
Appendix A. The PINN verification
This appendix provides supplementary evidence of the convergence of the PINN-based field inversion used in this work, and compares the inferred fields with those obtained from an adjoint-based inversion. The objective is to document the convergence behaviour of the PINN training, and to assess the plausibility of the inferred spatially varying coefficient field through comparison with a conventional partial differential equation (PDE) constrained inversion approach.
We show results for one periodic-hill case, and focus on the unconstrained inversion. The quantities shown are the reconstructed
$u$
,
$v$
and
$c_{b1}$
fields, where
$c_{b1}$
is the SA production coefficient inferred through field inversion. The results for the other two periodic-hill cases exhibit similar behaviour. Training the PINN for a single geometry took less than 3 hours on an RTX 4090 workstation. The adjoint solver used for comparison is DAFoam (He et al. Reference He, Mader, Martins and Maki2018, Reference He, Mader, Martins and Maki2020). The corresponding field-inversion process took approximately 36 core-hours to complete. These numbers are reported only to give a practical sense of the resources required by the two implementations as they presently exist. We do not re-implement the two approaches on a common hardware platform in this work, so this is not intended as a hardware-normalised cost comparison, particularly since one measure is GPU wall-clock time and the other is CPU core-hours.
Figure 12 reports the evolution of the training loss as a function of training iteration. In the present case, the total loss decreases and plateaus at approximately
$4\times 10^{-3}$
after 10 000 epochs, and remains at essentially the same level from 10 000 to 30 000 epochs.
Training history of the PINN: evolution of the loss as a function of training iteration.
Figure 12 Long description
The plot shows the PINN loss as a function of training iteration for a representative periodic-hill inversion. The total loss decreases during training and reaches a plateau of about 4 times 10 to the minus 3 after roughly 10,000 epochs, remaining nearly unchanged through 30,000 epochs.
Unconstrained field inversion: comparison between PINN-based inversion and adjoint-based inversion for a representative periodic-hill case. Panels show
$u$
,
$v$
and
$c_{b1}$
: (a–c) adjoint-based field inversion, (d–f) PINN-based field inversion.
Figure 13 Long description
The multi-panel contour figure compares adjoint-based and PINN-based field inversion. The top row shows adjoint-based fields and the bottom row shows PINN-based fields. The quantities are streamwise velocity u, wall-normal velocity v, and the inferred c_b1 field. The two inversion strategies produce broadly similar fields and comparable agreement with DNS.
To compare against a conventional inversion approach, figure 13 summarises the inferred fields from the PINN-based inversion and the adjoint-based inversion. The figure shows contour plots of
$u$
,
$v$
and
$c_{b1}$
obtained using the PINN inversion and the adjoint inversion. These comparisons provide a direct assessment of whether the PINN inversion yields fields consistent with those obtained from a conventional PDE-constrained inversion approach for this configuration. We observe that the inferred fields are broadly similar. We also compare the inverted velocity fields against the DNS velocity field through the MSE. For the representative case shown in figure 13, the MSE relative to DNS is
$0.003$
for the PINN-based inversion, and
$0.006$
for the adjoint-based inversion, indicating that both inversion strategies recover velocity fields close to DNS for this case.
Despite the favourable comparison obtained here, we do not claim that the PINN-based inversion is generally equivalent to adjoint-based inversion for arbitrary flows. In the present case, the mean-flow fields for the periodic hill are smooth and do not contain sharp gradients, which likely makes them easier to represent with a moderately sized neural network.
Appendix B. Summary of data used in the present study
This appendix summarises the reference datasets and flow parameters used in the present study in table 2. Several canonical turbulent-flow datasets are employed to train and assess the proposed framework. The periodic-hill cases with
$\alpha =0.5$
,
$0.8$
and
$1.2$
serve as the primary training cases for the inverse problem. Additional periodic-hill configurations, plane-channel flow, the backward-facing step, and the BeVERLI hill are used to evaluate cross-scenario generalisation. For the channel-flow validations, the
$ \textit{Re}_{\tau }\approx 5200$
case uses DNS reference data, whereas the
$ \textit{Re}_{\tau }\approx 10^6$
case uses canonical reference scaling rather than DNS.
Summary of the datasets used in the present study.
Table 2 Long description
The table lists the reference flow cases used in the study, their Reynolds-number definitions, their role in training or validation, and the reference data source. The training set consists of periodic hills with alpha values 0.5, 0.8, and 1.2 at Re_b = 5600. Validation cases include periodic hills with alpha values 1.0 and 1.5, plane channel flows at Re_tau about 5200 and 10^6, a backward-facing step at Re_b about 5600, and a BeVERLI hill at Re_H about 15000.
Appendix C. Reference data quality and numerical verification
The CFD information used in this paper consists of high-fidelity DNS/WRLES reference datasets and low-fidelity RANS calculations generated for model deployment. The high-fidelity reference datasets used here are widely used comparison cases in turbulence-model verification and validation, but the detailed verification of those simulations is documented in the original sources rather than in the present study. Our goal here is therefore to summarise the relevant reported evidence, and to document the numerical fidelity of the low-fidelity RANS calculations generated in the present work. For the periodic hills, the DNS database of Xiao et al. (Reference Xiao, Wu, Laizet and Duan2020) is itself accompanied by mesh-sensitivity and validation studies against benchmark computations for the same configuration. The reported meshes contain
$704$
–
$832$
points in the streamwise direction,
$385$
points in the wall-normal direction, and
$128$
points in the spanwise direction, with first-cell height
$\Delta y_{w\textit{all}}/H \approx 6.22\times 10^{-3}$
. Xiao et al. (Reference Xiao, Wu, Laizet and Duan2020) noted that this resolution is comparable to established benchmark simulations at the same Reynolds number, and documented that the mean-flow statistics are insensitive to further mesh refinement for the cases considered. For the BeVERLI hill, the high-fidelity reference is the WRLES dataset of Nair et al. (Reference Nair, Wadhai, Kunz and Yang2026). The BeVERLI configuration has also been used extensively in recent experimental and RANS-validation studies (Lowe et al. Reference Lowe, Roy, Devenport, Borgoltz, Grzyb, Shanmugam, Borole and Gargiulo2024; Roy et al. Reference Roy, Lowe, Devenport, Borgoltz, Grzyb, Borole, Shanmugam, HogePatil and Gargiulo2024; Nair et al. Reference Nair, Wadhai, Kunz and Yang2026). The WRLES use a structured surface-orthogonal grid with
$940\times 95\times 662$
points (58.3 million total). At the inlet, the reported spacings are
$\Delta x^+=18$
,
$\Delta z^+=12$
and
$\Delta y^+\lt 1$
at the first off-wall point.
The low-fidelity CFD calculations used for model deployment are the RANS simulations generated in the present work. For these cases, we use standard benchmark meshes drawn from the literature or NASA TMR resources. In particular, the periodic-hill deployment calculations use the benchmark meshes of Xiao et al. (Reference Xiao, Wu, Laizet and Duan2020) together with the same standardised steady-RANS solver settings documented by McConkey, Yee & Lien (Reference McConkey, Yee and Lien2021), and for the BeVERLI hill, the RANS grid matches the surface-orthogonal WRLES grid described above (Nair et al. Reference Nair, Wadhai, Kunz and Yang2026). In cases where a DNS-quality or benchmark mesh is already used for a RANS solve, the grid is substantially finer than what is typically required for RANS. Even so, we include representative convergence evidence for the cases in which the present work relies most directly on its own deployment meshes.
Representative grid-convergence study for two periodic-hill RANS calculations. Streamwise-velocity profiles are shown on coarse, baseline and fine meshes obtained by coarsening each direction by factor
$1.5$
, and refining by factor
$2$
. (a) Periodic hill with
$\alpha =1.0$
. (b) Periodic hill with
$\alpha =1.5$
.
Figure 14 Long description
The two panels compare streamwise velocity profiles computed on coarse, baseline, and fine meshes for periodic-hill cases with alpha equal to 1.0 and 1.5. The profiles collapse closely across mesh resolutions, indicating that the baseline mesh is adequate for the deployment calculations.
Representative iterative convergence history for the periodic-hill case at
$\alpha =1.0$
on the fine mesh. Residuals of
$U_1$
,
$U_2$
,
$p$
and
$\tilde {\nu }$
all fall below
$10^{-10}$
.
Figure 15 Long description
The plot shows residual histories for the velocity components U_1 and U_2, pressure p, and modified eddy viscosity nu_tilde on the fine mesh. All residuals fall below 10 to the minus 10, documenting tight iterative convergence for the representative deployment calculation.
For grid convergence, figure 14 reports representative streamwise-velocity profiles for two periodic-hill RANS calculations using coarse, baseline and fine meshes obtained by coarsening each direction by factor
$1.5$
, and refining by factor
$2$
. The profiles collapse closely across the three grids, confirming that the baseline mesh is adequate for the present deployment calculations. For iterative convergence, figure 15 reports a representative residual history for the periodic-hill case at
$\alpha =1.0$
on the fine mesh. With residuals reduced below
$10^{-10}$
, this is tighter than the
$10^{-6}$
level reported by McConkey et al. (Reference McConkey, Yee and Lien2021) as generally sufficient for the corresponding benchmark cases, and is in line with their most converged cases, many of which fall below
$10^{-8}$
. The other deployment cases exhibit similar convergence behaviour.
Appendix D. Additional channel-flow diagnostic
In this appendix, we include channel-flow results at
$ \textit{Re}_{\tau }=5200$
,
$10^4$
,
$10^5$
and
$10^6$
. Figure 16 shows the velocity profiles for the baseline SA model and the unconstrained FIML model. We quantify the departure from the canonical LoW through the difference between the predicted centreline velocity and the value of the log law evaluated at the same wall-normal location. This quantity increases gradually from
$1.06$
to
$2.59$
,
$5.50$
and
$8.41$
as the Reynolds number increases from 5200 to
$10^4$
,
$10^5$
and
$10^6$
. Therefore, the unconstrained model does not fail abruptly at a single Reynolds number; rather, it progressively drifts away from the baseline channel-flow calibration as the Reynolds number increases.
Additional channel-flow diagnostic for the unconstrained model at
$ \textit{Re}_{\tau }=5200$
,
$10^4$
,
$10^5$
and
$10^6$
. The solid black line denotes the canonical LoW. The dashed lines show unconstrained FIML results, and the solid lines show constrained FIML results.
Figure 16 Long description
The figure compares velocity profiles for the baseline model and unconstrained FIML at Re_tau values of 5200, 10^4, 10^5, and 10^6, with the canonical law of the wall shown for reference. The departure of the unconstrained model from the canonical law grows gradually with Reynolds number, while the constrained model preserves the log-law behaviour.
Appendix E. Sensitivity analysis of the training data
In this appendix, we examine the sensitivity of the constrained FIML model to the training dataset. For this purpose, the neural network is trained using only the periodic-hill case at
$\alpha = 1.2$
. Since the periodic-hill predictions exhibit only minor differences, and since the constrained FIML model preserves the LoW by construction, we restrict attention to the backward-facing step and the BeVERLI hill.
Skin-friction coefficient
$C_{\kern-1pt f}$
for the backward-facing step.
Figure 17 Long description
The appendix figure shows skin-friction coefficient for the backward-facing step when the constrained FIML model is trained only on the alpha equal to 1.2 periodic-hill case. The result remains similar to the model trained on all three periodic-hill cases, indicating limited sensitivity to the size of the periodic-hill training set for this test.
Skin-friction coefficient of the BeVERLI hill at
$30^\circ$
rotation along (a)
$z=0$
(centreline) and (b)
$z=-H/4$
.
Figure 18 Long description
The appendix figure shows BeVERLI-hill skin-friction coefficient at z = 0 and z = -H/4 when the constrained FIML model is trained only on the alpha equal to 1.2 periodic-hill case. The predictions are similar to those obtained from training on all three periodic-hill cases, indicating limited sensitivity to the training subset for this test.
Figure 17 shows the skin-friction results for the backward-facing step, and figure 18 shows the corresponding results for the BeVERLI hill. The constrained FIML model trained solely on the
$\alpha = 1.2$
periodic-hill case performs very similarly to the model trained on the three periodic-hill cases with
$\alpha = 0.5$
,
$0.8$
and
$1.2$
.
cb1
Lx×Ly×Lz=36.13H×9.9H×4H
6.73H
W=5H
α=1.2
cb1
νt
α=1.2
cb1
α=1.0
α=1.5
Reτ=5200
Reτ≈106
Reτ=5200
Reτ≈106
Cf
cb1
cb1
30∘
z=0
z=−H/4
u
v
cb1
1.5
2
α=1.0
α=1.5
α=1.0
U1
U2
p
ν~
10−10
Reτ=5200
104
105
106
Cf
30∘
z=0
z=−H/4
Open access
