# Turbulence closure for the compressible RANS equations assisted by field-inversion and machine-learning

Bartolomeo Fanizza, Pedro Stefanin Volpiani, Florent Renac, Denis Sipp

► To cite this version:

Bartolomeo Fanizza, Pedro Stefanin Volpiani, Florent Renac, Denis Sipp. Turbulence closure for the compressible RANS equations assisted by field-inversion and machine-learning. ETMM15, Sep 2025, Dubrovnik, Croatia. pp.202-207. hal-05333743

HAL Id: hal-05333743

https://hal.science/hal-05333743v1

Submitted on 27 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# TURBULENCE CLOSURE FOR THE COMPRESSIBLE RANS EQUATIONS ASSISTED BY FIELD-INVERSION & MACHINE-LEARNING

Bartolomeo Fanizza¹, Pedro Stefanin Volpiani¹, Florent Renac¹, Denis Sipp¹

¹ Aerodynamics, Aeroelasticity and Acoustics Department (DAAA), ONERA

bartolomeo.fanizza@onera.fr

## Abstract

This work presents a field-inversion (FI) and machine learning (ML)-augmented turbulence closure model for Reynolds-averaged Navier-Stokes (RANS) equations, targeting shock-wave/turbulent boundary-layer interaction (SWTBLI) flows. A database of optimal corrections is generated via adjoint-based FI across representative SWTBLI configurations. An artificial neural network (ANN) is then trained to regress these corrections for unseen cases.

Special emphasis is placed on input feature design and on the training loss definition: rotational and Galilean invariances are enforced, and a hybrid local/non-local damping function ensures smooth transitions at turbulent/non-turbulent interfaces. The ANN is trained with a custom loss function that promotes robustness and invariance to noisy inputs.

The database is used to perform both interpolation and extrapolation scenarios, enabling a comprehensive assessment of the ANN's capabilities and limitations. Full exploitation of the available data ensures that the resulting FI-ML augmented model improves turbulence closure across weakly to strongly separated SWTBLI flows, yielding enhanced predictions of the interaction zone and accurate boundary-layer resolution in both interpolative and extrapolative regimes.

## 1 Introduction

Reynolds-averaged Navier-Stokes (RANS) simulations are crucial in aerospace due to their low computational cost, making them widely used in aircraft and turbo-machinery design. However, they rely on mathematical closures for turbulence effects, limiting their accuracy in flows with separation, high streamline curvature, strong pressure gradients, or intense turbulent mixing.

In recent years, the integration of ML techniques into turbulence modeling has garnered significant attention. These state-of-the-art approaches have demonstrated remarkable potential in enhancing RANS models, positioning ML-augmented models as a powerful tool in computational fluid dynamics (CFD). Early studies, such as those by Ling and Templeton (2016) and Wang et al. (2017), have applied ML algorithms to directly learn unresolved turbulent

quantities or their discrepancies in the RANS equations from high-fidelity training data. An alternative approach, employed by Parish and Duraisamy (2016) and Volpiani et al. (2021), refines existing transport-equation-based turbulence models—such as $k-\epsilon$, $k-\omega$ or Spalart-Allmaras —by solving field-inversion (FI) problems and training ML models on the resulting correction fields.

In this study, we extend the Field Inversion-Machine Learning (FI-ML) framework introduced by Volpiani et al. (2021), which focused on correcting the Reynolds force vector in incompressible RANS equations, to address compressible flow regimes. The enhanced methodology is applied to a class of SWTBLI configurations, encompassing cases with varying degrees of flow separation, from mild to severe, which are known to pose significant challenges for traditional RANS models. Training cases are selected using efficient distance-based metrics to ensure both representativeness of the flow physics and minimal reliance on large training datasets. The relevance of input features (IFs) within the ML model is quantified using a suite of complementary evaluation metrics. Furthermore, the generalization capabilities of the FI-ML model are examined through both data-dependent and data-agnostic measures, and potential correlations between these metrics and a posteriori predictive performance are explored.

## 2 Baseline model problem and discretization

### RANS equations

In the following, we consider the discretization of the steady-state compressible Reynolds-averaged Navier-Stokes (RANS) equations using a Discontinuous Galerkin (DG) formulation. Let $\Omega \subset \mathbb{R}^d$ denote a bounded domain, where $d$ is the spatial dimension. The governing equations are given by:

$$ \nabla \cdot \mathbf{f}_c(\mathbf{u}) - \nabla \cdot \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \mathbf{s}(\mathbf{u}, \nabla \mathbf{u}) \quad (1) $$

in $\Omega$, subject to appropriate boundary conditions on $\partial\Omega$. The vector $\mathbf{u} = (\rho, \rho\mathbf{v}^\top, \rho E, \rho\tilde{\nu})^\top$ represents the conservative variables where $\rho$ is the density, $\mathbf{v} = (u, v, w)^\top$ is the velocity vector, and $E = e + \frac{1}{2}\|\mathbf{v}\|^2$
is the total specific energy, with e denoting the specific internal energy. The nonlinear convective and diffusive fluxes in Eq.(1) are defined respectively as:

$$ \mathbf{f}_c(\mathbf{u}) = \begin{pmatrix} \rho \mathbf{v}^\top \\ \rho \mathbf{v} \mathbf{v}^\top + p \mathbf{I} \\ \rho H \mathbf{v}^\top \\ \rho \tilde{\nu} \mathbf{v}^\top \end{pmatrix}, $$

and

$$ \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \begin{pmatrix} 0 \\ \tau + \tau_r \\ \mathbf{v}^\top(\tau + \tau_r) - \mathbf{q}^\top - \mathbf{q}_t^\top \\ \frac{1}{\sigma} \zeta \nabla \tilde{\nu}^\top \end{pmatrix}. $$

Here, $H = E + \frac{p}{\rho}$ is the total specific enthalpy, and $p$ is the pressure, defined by the ideal gas law $p = (\gamma - 1)\rho e$, where $\gamma$ is the ratio of specific heats. The viscous stress tensor $\tau$ and the conductive heat flux vector $\mathbf{q}$ account for molecular transport, while the Reynolds stress tensor (RST) $\tau_r$ and turbulent heat flux $\mathbf{q}_t$ are modeled using a Boussinesq-like approximation:

$$ \begin{aligned} \tau_r &= \mu_t \left( -\frac{2}{3}(\nabla \cdot \mathbf{v})\mathbf{I} + \nabla \mathbf{v} + \nabla \mathbf{v}^\top \right), \\ \mathbf{q}_t &= -\frac{\mu_t C_p}{Pr_t} \nabla T, \end{aligned} $$

where $\mu_t$ is the turbulent dynamic viscosity, $C_p$ is the specific heat at constant pressure, and $Pr_t = 0.9$ is the turbulent Prandtl number.

### Spalart-Allmaras one-equation model

In this work, the Spalart-Allmaras (SA) one-equation turbulence model, originally proposed by Spalart and Allmaras (1992), is employed. The turbulent dynamic viscosity is expressed through the transported variable $\tilde{\nu}$ as:

$$ \mu_t = \rho \tilde{\nu} f_{v1}, \quad f_{v1} = \frac{\chi^3}{\chi^3 + c_{v1}^3}, \quad \chi = \frac{\rho \tilde{\nu}}{\mu}. \quad (2) $$

The source terms in Eq.(1) are given by $s = (0, 0, 0, s_{\tilde{\nu}})^\top$, where

$$ s_{\tilde{\nu}}(\mathbf{u}, \nabla \mathbf{u}) = c_{b1} \tilde{S} \rho \tilde{\nu} - c_{w1} f_w \frac{\rho \tilde{\nu}^2}{d_w^2} + \frac{c_{b2} \rho}{\sigma} \nabla \tilde{\nu} \cdot \nabla \tilde{\nu}, $$

with $\tilde{S} = S + \bar{S}$, $S = \|\nabla \times \mathbf{v}\|_2$, and $\bar{S} = \frac{\tilde{\nu} f_{v2}}{\kappa^2 d_w^2}$. Here, $d_w$ denotes the distance to the wall and $\xi = \mu + \rho \tilde{\nu}$ is the effective diffusive coefficient. The auxiliary functions are defined as:

$$ \begin{aligned} f_{v2} &= 1 - \frac{\chi}{1 + \chi f_{v1}}, & f_w &= g \left( \frac{1 + c_{w3}^6}{g^6 + c_{w3}^6} \right)^{\frac{1}{6}}, \\ g &= r + c_{w2} (r^6 - r), \\ r &= \frac{\tilde{\nu}}{\tilde{S} \kappa^2 d_w^2}, & c_{w1} &= \frac{c_{b1}}{\kappa^2} + \frac{1 + c_{b2}}{\sigma}. \end{aligned} $$

For completeness, we give the values of the remaining involved constants: $\sigma = \frac{2}{3}$, $c_{b1} = 0.1355$, $c_{b2} = 0.622$, $\kappa = 0.41$, $c_{w2} = 0.3$, $c_{w3} = 2$, and $c_{v1} = 7.1$. Negative values of the eddy viscosity in Eq. (2) are limited using a negative version of the SA model as advocated in Allmaras et al. (2012).

We also include the correction proposed by Spalart and Garbaruk (2020) for shear flows at low Reynolds numbers on the constant $c_{w2}$, changed to the following

$$ c_{w2,LRe} = c_{w4} + \frac{c_{w5}}{\left[ \left( \frac{\chi}{40} \right) + 1 \right]^2} $$

with $c_{w4} = 0.21$ and $c_{w5} = 1.5$. The effect of $c_{w2,LRe}$ vanishes for large Reynolds-number values and therefore the calibration of the original model is still preserved for high Reynolds number shear flows.

### Discontinuous Galerkin formulation

The DG method consists in defining a discrete weak formulation of the problem (1) with initial and boundary conditions. The domain $\Omega$ is partitioned into a shape-regular mesh $\Omega_h$ consisting of nonoverlapping and nonempty elements $\kappa$. We look for approximate solutions in the function space of discontinuous polynomials of total degree at most $p$ in order to obtain a $(p+1)$-order accuracy in the spacial discretization. The numerical solution of the problem (1) is sought under the form

$$ \mathbf{u}_h(\mathbf{x}) = \sum_{l=1}^{N_p} \phi_{\kappa}^l(\mathbf{x}) \mathbf{U}_{\kappa}^l, \quad \forall \mathbf{x} \in \kappa, \kappa \in \Omega_h, $$

where $(\mathbf{U}_{\kappa}^l)_{1 \le l \le N_p}$ are the degrees of freedom (DOFs) in the element $\kappa$. The subset $(\phi_{\kappa}^l)_{1 \le l \le N_p}$ constitutes a hierarchical and orthogonal modal basis of the function space of discontinuous polynomials, restricted onto the element $\kappa$ and $N_p = \frac{(p+d)!}{p!d!}$ is its dimension. A polynomial representation of degree $p = 1$ is employed to achieve second-order spatial accuracy throughout the numerical frameworks considered in this study.

## 3 Field-Inversion for RANS simulation

In the FI framework, we introduce a corrective parameter into the data-augmented (DA) compressible RANS equations:

$$ \nabla \cdot \mathbf{f}_c(\mathbf{u}) - \nabla \cdot \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \mathbf{s}(\mathbf{u}, \nabla \mathbf{u}) + \rho \mathbf{f}_{DA}, $$

where the term, $\mathbf{f}_{DA} = (0, \mathbf{f}_{xyz}, \mathbf{v}^\top \mathbf{f}_{xyz})$, is introduced to account for modeling errors associated with the divergence of the RST. The DOFs of the corrective term $\mathbf{f}_{xyz}$ are optimized using a discrete-adjoint approach, as detailed in the work of Fanizza et al. (2025). The objective function of the FI problem is a $L^2$-norm of the mean flow velocity fields of the DA-RANS solution and a reference solution, weighted by the wall distance in order to focus the correction in the shear region. A $L^2$-norm penalty term is also included in the objective function order to penalize intrusive correction and regularize the optima space.
# 4 Machine-Learning

## Architecture & Training Process

In the present framework, fully connected ANNs are employed to approximate a nonlinear mapping $\mathcal{N}$ from the so-called input features (IFs) space of size $N_I$, denoted by $\eta = \{\eta_i\}_{1\le i\le N_I}$, to the optimal forcing terms $f_{xyz}^{\text{opt}}$ identified by the FI framework. The resulting prediction is written as $f_{xyz}^{\text{NN}} = \mathcal{N}(\eta)$.

The ML training process involves minimizing a custom loss function composed of three terms: (i) a normalized mean-square error (MSE) between the predicted output $f_{xyz}^{\text{NN}}$ and the reference $f_{xyz}^{\text{opt}}$; (ii) an $L^2$ weight decay (WD) regularization term to constrain the complexity of the model parameters; and (iii) an adaptive noise stabilization (adaNS) term, introduced to reduce the model's sensitivity to input noise. The adaNS term is formulated as an MSE between the nominal ML prediction $f_{xyz}(\eta)$ and its counterpart evaluated with perturbed inputs, $f_{xyz}(\eta + \delta\eta)$. Its contribution is modulated by a hyperparameter $\alpha_{\text{NS}}(\chi_t)$, which is adapted based on the turbulent intensity $\chi_t = \frac{\mu_t}{\mu}$. This dynamic weighting ensures strong penalization in laminar or weakly turbulent regions, while maintaining a mild regularization in highly turbulent zones. Such an adaptive formulation improves the robustness of the model to small perturbations in the IFs throughout the entire flow field. Notably, using a fixed penalty strength—such as in the WD term—either results in under-regularization in sensitive regions or imposes excessive constraints, ultimately degrading predictive performance.

Finally, since the IFs are non-dimensional by construction, a non-trainable layer, at the end of the architecture, rescales the prediction using the standard deviation of the target values computed from the training dataset.

## Input Features formulation

A critical assessment of the IFs design process is conducted, incorporating considerations from both ML and classical turbulence modeling methodologies, as reviewed by Spalart (2015). To enhance the generalization capability of the model, both local and invariant IFs are formulated. Furthermore, a damping function is integrated into the feature normalization scheme, defined as $\eta(\alpha, \beta) = \frac{\alpha}{\alpha+\beta+D_{\infty}\epsilon_t}$, to ensure the safe deactivation of features in non-turbulent regions. This formulation also prevents unintended activations due to ill-posed expressions, particularly when both numerator and denominator quantities approach zero. The damping function

$$ \epsilon_t = \left( \frac{40}{\chi+1} \right)^2 \qquad (3) $$

is inspired by the low Reynolds number correction in the SA turbulence model, when both quantities $\alpha$ and $\beta$ approach zero, $\epsilon_t$ attains a finite positive value, thereby regularizing the expression. Conversely, in

turbulent regions, $\epsilon_t$ remains negligibly small. Dimensional consistency is restored via the inclusion of the free-stream scaling parameter $D_{\infty}$. Although this approach is not strictly optimal, it ensures that the feature exhibits non-local dependence on the mean flow variables only when inactive, while recovering an asymptotically local behavior in regions where the feature becomes active, thereby endowing it with a mixed local/non-local formulation.

<table><caption>Table 1: Full IF set employed in the ANN model.</caption><thead><tr><th>Formula</th><th>Additional information</th></tr></thead><tbody><tr><td>η<sub>1</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||S||</div><div style="padding: 0 5px;">√||Ω||<sup>2</sup> + ||S||<sup>2</sup> + S<sub>∞</sub>ε<sub>t</sub></div></div></td><td>S<sub>∞</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||S||<sub>∞</sub></div><div style="padding: 0 5px;">δ<sub>x,inlet</sub></div></div></td></tr><tr><td>η<sub>2</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||Ω||</div><div style="padding: 0 5px;">√||Ω||<sup>2</sup> + ||S||<sup>2</sup> + S<sub>∞</sub>ε<sub>t</sub></div></div></td><td>S<sup>*</sup> = √∂<sub>x<sub>i</sub></sub>u<sub>j</sub>∂<sub>x<sub>i</sub></sub>u<sub>j</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">ν<sub>t</sub>f<sub>x<sub>i</sub>y<sub>j</sub></sub></div><div style="padding: 0 5px;">π<sup>2</sup>d<sub>w</sub></div></div></td></tr><tr><td>η<sub>3</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||S-Ω||</div><div style="padding: 0 5px;">√||Ω-Ω||<sup>2</sup> + ||S||<sup>2</sup>||Ω||<sup>2</sup> + S<sub>∞</sub>ε<sub>t</sub></div></div></td><td>r = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">ν<sub>t</sub></div><div style="padding: 0 5px;">(S<sup>*</sup> + S<sub>∞</sub>ε<sub>t</sub>)π<sup>2</sup>d<sub>w</sub><sup>2</sup></div></div></td></tr><tr><td>η<sub>4</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">ν<sub>t</sub><sup>2</sup></div><div style="padding: 0 5px;">(ν<sub>t</sub> + c<sub>v<sub>1</sub>ν</sub>)<sup>3</sup></div></div></td><td>r<sub>d</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">ν<sub>t</sub>+<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">π</div><div style="padding: 0 5px;">d<sub>w</sub></div></div></div></div></td></tr><tr><td>η<sub>5</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">r</div><div style="padding: 0 5px;">1+r</div></div></td><td>P<sub>t</sub> = c<sub>b1</sub>S̃ρ̃ν̃</td></tr><tr><td>η<sub>6</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">r<sub>d</sub></div><div style="padding: 0 5px;">1+r<sub>d</sub></div></div></td><td>D<sub>t</sub> = c<sub>w1</sub>f<sub>w</sub>ρ̃<sup>2</sup>d<sub>w</sub><sup>2</sup></td></tr><tr><td>η<sub>7</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">P<sub>t</sub></div><div style="padding: 0 5px;">|P<sub>t</sub>| + |D<sub>t</sub>ρ̃ν̃| + S<sub>∞</sub>ρ<sub>∞</sub>ρ̃ν̃<sub>∞</sub>ε<sub>t</sub></div></div></td><td>Π<sub>t</sub> = ρ<sup>c<sub>v<sub>2</sub>σ</sub></sup>∂<sub>x<sub>k</sub></sub>ν̃∂<sub>x<sub>k</sub></sub>ν̃</td></tr><tr><td>η<sub>8</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">D<sub>t</sub></div><div style="padding: 0 5px;">|D<sub>t</sub>| + |D<sub>t</sub>ρ̃ν̃| + S<sub>∞</sub>ρ<sub>∞</sub>ρ̃ν̃<sub>∞</sub>ε<sub>t</sub></div></div></td><td>e = E - <sup>1</sup>⁄<sub>2</sub>||v||<sup>2</sup></td></tr><tr><td>η<sub>9</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">Π<sub>t</sub></div><div style="padding: 0 5px;">|Π<sub>t</sub>| + |D<sub>t</sub>ρ̃ν̃| + S<sub>∞</sub>ρ<sub>∞</sub>ρ̃ν̃<sub>∞</sub>ε<sub>t</sub></div></div></td><td>P<sub>S,t</sub> = c<sub>b1</sub>||∇ × v||<sup>2</sup>ρ̃ν̃</td></tr><tr><td>η<sub>10</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">P<sub>S,t</sub></div><div style="padding: 0 5px;">P<sub>S,t</sub> + S<sub>∞</sub>ρ<sub>∞</sub>ρ̃ν̃<sub>∞</sub>ε<sub>t</sub></div></div></td><td>P<sub>S,t</sub> = c<sub>b1</sub><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">ν<sub>t</sub>f<sub>x<sub>k</sub>y<sub>j</sub></sub></div><div style="padding: 0 5px;">π<sup>2</sup>d<sub>w</sub><sup>2</sup></div></div>ρ̃ν̃</td></tr><tr><td>η<sub>11</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">k<sub>QCR-2013V</sub></div><div style="padding: 0 5px;">k<sub>QCR-2013V</sub> + 0.01ε</div></div></td><td>k<sub>QCR-2013V</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">3</div><div style="padding: 0 5px;">c<sub>2</sub>ρ<sub>∞</sub>ν<sub>t</sub>||Ω||</div></div></td></tr><tr><td>η<sub>12</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">2</div><div style="padding: 0 5px;">k<sub>QCR-2013V</sub></div></div></td><td>c<sub>s</sub> = √<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">P</div><div style="padding: 0 5px;">ρ</div></div></td></tr><tr><td>η<sub>13</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">√k<sub>QCR-2013V</sub></div><div style="padding: 0 5px;">√k<sub>QCR-2013V</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">ν</div><div style="padding: 0 5px;">d<sub>w</sub></div></div> + U<sub>∞</sub>ε<sub>t</sub></div></div></td><td>U<sub>∞</sub> = ||v<sub>∞</sub>||</td></tr><tr><td>η<sub>14</sub> = 2√l<sub>cc</sub></td><td>l<sub>cc</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">d<sub>w</sub></div><div style="padding: 0 5px;">√k<sub>QCR-2013V</sub> + U<sub>∞</sub><sup>2</sup>ε<sub>t</sub></div></div></td></tr></tbody></table>

The full IFs set is summarized in table 1, along with the necessary supplementary information required for their implementation. For completeness, $k_{QCR-2013V}$ denotes an approximation of the turbulent kinetic energy employed in the SA-QCR-2017V model, detailed in Mani et al. (2013), with the model constant set to $c_{r2} = 2.5$. The speed of sound is denoted by $c_s$, and $\delta_{x,inlet}$ indicates the boundary layer thickness at the domain inlet. The velocity scale $v^* = v_\infty - v_b$, defined as the difference between the free-stream and wall-adjacent velocities, is used in the normalization factor $D_\infty$ to ensure Galilean invariance. The IFs, adapted from the literature and modified with a damping function, reduce instabilities and spurious activations seen in their original form, as exemplified for $\eta_2$ in fig. 1.

# 5 Numerical results

## FI results

The FI methodology described above is applied to SWTBLI configurations with varying shock deflection angles, $\phi$, while keeping the inlet Mach number fixed at $M = 2.28$. Two momentum-thickness Reynolds numbers are considered: a moderate case at $Re_\theta = 1047$ and a high-Reynolds-number (hRe) case at $Re_\theta = 1990$. The high-fidelity reference data are obtained from the direct numerical simulations (DNS) of Volpiani et al. (2018).
<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison of feature $\eta_2$ in the SWTBLI with strong separation: original formulation from the literature (left), modified version using a constant tolerance with units of inverse time [$T^{-1}$] (middle), and the proposed formulation incorporating the damping function (right). A schematic of key flow features under strong separation is also included.</figcaption></figure>

<table><caption>Table 2: Testing and training cases in the database for the three scenarios. Cases with $\varphi = 0.0°$ represent fully developed turbulent boundary layers without shock-wave interactions.</caption><thead><tr><th>Scenario</th><th></th><th>φ - 0.0°</th><th>φ - 0.0° - hRe</th><th>φ - 5.0°</th><th>φ - 6.5°</th><th>φ - 8.0°</th><th>φ - 9.5°</th><th>φ - 9.5° - hRe</th><th>φ - 11.0°</th></tr></thead><tbody><tr><td rowspan="2">NN1</td><td>Objective:</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td><strong>Interpolate</strong></td><td>Train</td><td>Train</td><td>Train</td></tr><tr><td>R<sup>2</sup>-score</td><td>99.6</td><td>97.9</td><td>99.1</td><td>98.5</td><td>91.9</td><td>96.8</td><td>97.0</td><td>96.3</td></tr><tr><td rowspan="2">NN2</td><td>Objective:</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td><strong>Extrapolation</strong></td></tr><tr><td>R<sup>2</sup>-score</td><td>99.5</td><td>97.8</td><td>98.9</td><td>98.6</td><td>98.5</td><td>95.8</td><td>96.8</td><td>78.9</td></tr><tr><td rowspan="2">NN3</td><td>Objective:</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td><td>Train</td></tr><tr><td>R<sup>2</sup>-score</td><td>99.6</td><td>98.3</td><td>99.0</td><td>98.8</td><td>98.8</td><td>96.5</td><td>98.3</td><td>95.9</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Skin friction coefficient ($C_f$) distributions from DNS, SA-RANS, and DA-RANS simulations for the SWTBLI with $\varphi = 9.5°$.</figcaption></figure>

In fig. 2, the skin friction coefficient $C_f$ obtained from SA-RANS and DA-RANS simulations is compared against DNS reference data for a SWTBLI at $\varphi = 9.5°$, representative of moderate/strong flow separation. The DA-RANS solutions demonstrate notable improvements across the whole SWTBLI database.

## ML a priori results

To enhance stability and facilitate predictions in canonical flow regions, two turbulent boundary layers without shock waves ($\varphi = 0.0°$) at Reynolds numbers matching those of the SWTBLI cases are incorporated into the training data. Three ANNs are subsequently trained to investigate distinct scenarios: NN1 focuses on interpolation at a moderate deflection angle $\varphi = 8.0°$; NN2 addresses extrapolation at the highest deflection angle $\varphi = 11.0°$; and NN3 is trained on the full dataset available. The datasets and corresponding training objectives are summarized in table 2, together with the $R^2$ scores on each test case. All ANNs were

trained with the same hyperparameters: 8 hidden layers of 80 neurons each, identical WD and adaNS tuning coefficients, full-batch L-BFGS optimization over 1500 epochs, and early stopping using 20% validation data. The training set was under-sampled to balance samples between key flow regions (boundary layers, shear, recovery etc.) and free-stream zones.

The inclusion of the adaNS term improves learning performance across various $\varphi$ angles, as shown for NN2 in fig. 3. The figure compares the MSE for different $\varphi$ values using original inputs (solid lines) versus perturbed inputs (dashed lines). When adaNS is included in the training loss (blue lines and circles), the MSE remains closer after perturbation, indicating greater robustness—particularly for the boundary layer at $\varphi = 0.0°$ and the extrapolation case at $\varphi = 11.0°$. Similar trends are observed for higher Reynolds number data and NN1.

<figure><img src="image_4.png" /><figcaption>Figure 3: NN2: Effect of the adaNS penalty on the normalized MSE with and without a perturbation $\delta\eta$.</figcaption></figure>
ML a posteriori results: interpolation & extrapolation scenarios

Both interpolation and extrapolation scenarios are
considered in the a posteriori evaluation of the ML-
RANS solutions. The ANN model is integrated into
the Aghora CFD software, a high-order DG solver for
compressible turbulent flows developed at ONERA by
Renac et al. (2015). The ANN is periodically invoked
to update the prediction of f<sub>xyz</sub><sup>NN</sup> at a fixed frequency
relative to the pseudo time-marching scheme. No no-
table changes occurred for ANN calls every N<sub>T</sub> ≤ 5
iterations.

<figure><img src="image_5.png" /><figcaption>Figure 4: **NN1 (interpolation):** Skin friction $C_f$ and $\frac{p_w}{p_\infty}$ from DNS, SA-RANS, DA-RANS, and ML-RANS simulations at $\varphi = 8.0^\circ$.</figcaption></figure>

In figs. 4 and 5, the skin friction coefficient and wall pressure predicted by the NN1 and NN2-augmented models are shown for the interpolation and extrapolation scenarios, respectively. In both cases, the boundary layer remains well-resolved without degeneration. For the interpolation scenario at $\varphi = 8.0^\circ$ (fig. 4), the prediction of the interaction zone is notably improved, and the recovery region closely matches the DA and DNS results. In the extrapolation scenario at $\varphi = 11.0^\circ$ (fig. 5), the ML-augmented model agrees well with the DA solution in both the recovery region and the initial part of the recirculation bubble. The mismatch near the reattachment point likely arises from differences in the underlying flow physics: the training data mostly include weak separation cases, whereas the extrapolation case features strong separation, which is underrepresented in the training set.

**Global quantities of interest in SWTBLI**

To address the limitations of the training sets for
NN1 and NN2, the NN3 scenario is investigated to an-

<figure><img src="image_6.png" /><figcaption>Figure 5: **NN2 (extrapolation):** Skin friction $C_f$ and wall pressure $\frac{p_w}{p_\infty}$ from DNS, SA-RANS, DA-RANS, and ML-RANS simulations at $\varphi = 11.0^\circ$.</figcaption></figure>

alyze the behavior of flow separation across varying
shock-wave intensities. The trends of the reattachment
($x_R$) and separation ($x_S$) points as functions of deflec-
tion angle $\varphi$ are presented in fig. 6.

<figure><img src="image_7.png" /><figcaption>Figure 6: **NN3:** Separation ($x_S$) and reattachment ($x_R$) locations from DNS (black circles), SA-RANS (red gradients), DA-RANS (blue diamonds), and ML-RANS (green squares) across deflection angles. Filled and open green squares denote training and testing cases, respectively.</figcaption></figure>

Using the NN3 model yields significantly im-
proved predictions of $x_R$ and $x_S$ compared to SA-
RANS, showing closer agreement with DNS reference
data. Beyond the training set, NN3 is evaluated on
interpolation test cases $\varphi = \{6.0^\circ, 7.0^\circ, 9.0^\circ, 10.0^\circ\}$
and extrapolation cases representing mild separation
at $\varphi = \{4.0^\circ, 4.5^\circ\}$ and strong separation at $\varphi =
\{11.5^\circ, 12.0^\circ\}$. The predicted trends for separation
and reattachment points closely follow those observed
in DNS and DA results, outperforming SA-RANS in accuracy. Notably, the model exhibits consistent and coherent predictions even for weak separation scenarios at small $\phi$ values, demonstrating robust performance in cases with nearly absent recirculation bubbles. Additionally, incorporating strong separation cases around $\phi = 11.0^{\circ}$ enriches the model's understanding of strongly separated flows, enabling reliable extrapolated predictions of $x_R$ and $x_S$ at higher deflection angles. In general, ML-based RANS simulations accurately capture the location of the recirculation bubble, closely matching the DA-RANS results even for deflection angles $\phi$ outside the training range, confirming the generalizability of the model.

## 6 Conclusions

This work presents a FI-ML augmented turbulence closure for RANS simulations of SWTBLI flows. A database of optimal corrections is generated via adjoint-based field inversion and used to train ANNs for predictive deployment on unseen conditions. Key methodological innovations include the careful design of the IFs set that enforce rotational and Galilean invariance, enhanced by a hybrid local/non-local damping function to ensure numerical stability and prevent spurious activations. Furthermore, the introduction of an adaptive noise stabilization (adaNS) term in the training loss substantially improves model robustness against noisy inputs.

The FI-ML augmented closures accurately predict skin friction and wall pressure in both interpolation and extrapolation scenarios, preserving boundary layer resolution and improving interaction zone predictions. Although some mismatch appears near reattachment in strongly separated flows that are underrepresented in training, overall performance remains robust.

Significantly, the NN3 model, trained on an extended dataset, delivers improved predictions of separation and reattachment locations, closely matching DNS and DA reference data even outside of the training parameter space. This demonstrates the model's ability to generalize across a wide range of deflection angles, capturing both weak and strong separation phenomena and reliably predicting recirculation bubble dynamics.

In summary, the integration of adjoint-based FI with robust ML training yields a turbulence closure model that enhances RANS predictive capabilities for complex SWTBLI flows, offering improved accuracy, robustness, and generalization across varying flow regimes.

## Acknowledgments

This work was supported by ONERA under the internal project TURBINA (Modélisation de la turbulence assistée par intelligence artificielle).

## References

Fanizza B., Volpiani P.S., Renac F., Martin E., & Sipp D. (2025), Adjoint-based optimization for non-linear inverse problems with high-order discretization of the compressible RANS equations, *Applied Mathematical Modelling*, Vol. 142, pp. 115984.

Ling J., Kurzawski A., & Templeton J. (2016), Reynolds averaged turbulence modelling using deep neural networks with embedded invariance, *Journal Of Fluid Mechanics*, Vol. 807, pp. 155-166.

Mani, M., Babcock, D., Winkler, C. and Spalart, P., (2013), Predictions of a supersonic turbulent flow in a square duct. *In 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition* p. 860.

Parish, E., Ching, D.S., Jordan, C., Nicholson, G., Miller, N.E., Beresh, S., Barone, M., Gupta, N. & Duraisamy, K., (2025), Data-driven turbulent Prandtl number modeling for hypersonic shock-boundary-layer interactions. *AIAA Journal*, 63(5), pp.1671-1692.

Parish E., & Duraisamy K. (2016), A paradigm for data-driven predictive modeling using field inversion and machine learning, *Journal Of Computational Physics*, Vol. 305, pp. 758-774.

Renac, F., de la Llave Plata, M., Martin, E., Chapelier, J. B., & Couaillier, V. (2015), Aghora: A High-Order DG Solver for Turbulent Flow Simulations. Springer International Publishing, Vol. 128, pp. 315-335.

Spalart P. (2015), Philosophies and fallacies in turbulence modeling. *Progress In Aerospace Sciences*, 74, pp. 1-15.

Spalart, P., & Allmaras, S. (1992). A one-equation turbulence model for aerodynamic flows. 30th Aerospace Sciences Meeting and Exhibit, 439.

Spalart, P. R., & Garbaruk, A. V. (2020). Correction to the Spalart-Allmaras Turbulence Model, Providing More Accurate Skin Friction. *AIAA Journal*, 58(5), 1903-1905.

Volpiani P., Meyer M., Franceschini L., Dandois J., Renac F., Martin E., Marquet O., & Sipp D. (2021), Machine learning-augmented turbulence modeling for RANS simulations of massively separated flows, *Phys. Rev. Fluids*, Vol. 6, 064607.

Volpiani P., Bernardini M., & Larsson J. (2018), Effects of a nonadiabatic wall on supersonic shock/boundary-layer interactions, *Physical Review Fluids*, Vol. 3, 083401.

Wang J., Wu J., & Xiao H. (2017), Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data, *Phys. Rev. Fluids*, Vol. 2, 034603.