2.1. Data preparation

The fields in the numerical database are processed to prepare the data for the machine-learning wall models. The goal of this preprocessing step is to produce fields that are similar in some ways to the fields of a wall-modeled LES. Namely, the fields are filtered to attenuate large frequencies that cannot be resolved in a wall-modeled LES and resampled onto coarse tetrahedral meshes that could be used for WMLES computations. The present machine-learning methodology does not make any assumption regarding the nature of the function that should be learned by the model and can only operate in the range of mesh resolution that was seen during training. This operating range effectively determines the range of Reynolds number that can be addressed a posteriori at moderate computational cost. Since the exact mesh resolution required for unseen test cases is not known, it is important to consider a wide range of mesh refinements to train the model. For each simulation, 12 tetrahedral WMLES meshes with varying mesh resolution are generated using the mesh adaptation library MMG3D (Dobrzynski and Frey, Reference Dobrzynski and Frey2008; Dapogny et al., Reference Dapogny, Dobrzynski and Frey2014; Balarac et al., Reference Balarac, Basile, Bénard, Bordeu, Chapelier, Cirrottola, Caumon, Dapogny, Frey, Froehly, Ghigliotti, Laraufie, Lartigue, Legentil, Mercier, Moureau, Nardoni, Pertant and Zakari2021), with nominal edge length in the range $ 25\le {e}^{+}\le 80 $ , in wall units, where the classical wall-unit scaling is based on the friction velocity and the kinematic viscosity at the wall, namely $ {e}^{+}={eu}_{\tau }/{\nu}_{\omega } $ , with $ e $ the nominal edge length, $ {u}_{\tau }={\left(\tau /{\rho}_{\omega}\right)}^{0.5} $ the friction velocity, $ \tau $ the wall shear stress, $ {\rho}_{\omega } $ the wall density, and $ {\nu}_{\omega } $ the wall kinematic viscosity. This range of $ {e}^{+} $ implies that the first point off the wall of the WMLES mesh is typically within the logarithmic layer or the upper part of the buffer layer, in the case of a fully developed turbulent boundary layer. It is suitable for instance in turbomachinery-flow simulations (Leonard et al., Reference Leonard, Sanjose, Moreau and Duchaine2016; Dombard et al., Reference Dombard, Duchaine, Gicquel, Odier, Leroy, Buffaz, Le-Guyader, Démolis, Richard and Grosnickel2020; Odier et al., Reference Odier, Thacker, Harnieh, Staffelbach, Gicquel, Duchaine, Rosa and Müller2021). Indeed, while real turbomachines are typically not instrumented for boundary-layer measurements, experimental measurements, and high-fidelity LESs show that the friction Reynolds number based on boundary-layer thickness is in the range of 600–1000 in academic linear blade cascades with realistic operating conditions in terms of Reynolds and Mach numbers (Arts et al., Reference Arts, Rouvroit and Rutherford1990; Ma et al., Reference Ma, Ottavy, Lu, Francis and Gao2011; Gao et al., Reference Gao, Zambonini, Boudet, Ottavy, Lu and Shao2015; Zambonini et al., Reference Zambonini, Ottavy and Kriegseis2017; Dupuy et al., Reference Dupuy, Gicquel, Odier, Duchaine and Arts2020). However, this implies that the trained models are not applicable for cell sizes larger than $ {\Delta}^{+}=80 $ , which would make their use in flows with a very large friction Reynolds number ( $ {10}^4 $ or above) impractical. It would in principle be possible to train the model on a wider range of mesh resolutions. However, it should be noted that this range would remain limited by the Reynolds number of the high-fidelity simulations in the training database. For instance, using a direct numerical simulations of a channel with a friction Reynolds number of 10,000 (Hoyas et al., Reference Hoyas, Oberlack, Alcántara-Ávila, Kraheberger and Laux2022), the present method could be used to train a model that can operate up to a cell size in the order of a thousand wall unit, but not an order of magnitude more.

To filter the data, the instantaneous fields of the numerical database are first linearly interpolated on fine tetrahedral meshes. The resulting fields are used to compute the filtered variables on the coarse tetrahedral meshes. A surface filter is used for wall quantities while a volume filter is used outside the walls:

1. 1. The surface filter ( $ {\underset{\dot{\mkern6mu}}{-}}^S $ ) is defined as

   (1) $$ {\overline{\phi}}^S\left({p}_0\right)=\sum \limits_{p\in {\mathfrak{w}}_{p_0}^S}{C}_{p_0}^S{S}_p{e}^{-\frac{{\boldsymbol{r}}_p-{\boldsymbol{r}}_{p_0}}{2{\left({\sigma}_{p_0}^S\right)}^2}}\phi (p), $$
   with $ {\overline{\phi}}^S\left({p}_0\right) $ a filtered variable on a point $ {p}_0 $ of the target coarse mesh, $ \phi (p) $ the corresponding unfiltered variable on a point $ p $ of the source fine mesh, $ {C}_{p_0}^S $ a normalization constant, $ {S}_p $ the nodal area associated with node $ p $ on the source fine mesh, $ {\boldsymbol{r}}_p $ and $ {\boldsymbol{r}}_{p_0} $ the position vector associated with point $ p $ and $ {p}_0 $ , respectively, $ {\sigma}_{p_0}^S={\tilde{\Delta}}_{p_0}^S/\sqrt{12} $ the standard deviations of the surface Gaussian kernel, $ {\tilde{\Delta}}_{p_0}^S={\tilde{S}}_{p_0}^{1/2} $ the local filter width, $ {\tilde{S}}_{p_0} $ the nodal area of node $ {p}_0 $ on the target coarse mesh and $ {\mathfrak{w}}_{p_0}^S=\left\{p\in {\mathfrak{P}}_{\mathfrak{f}}:\left(\parallel {\boldsymbol{r}}_p-{\boldsymbol{r}}_{p_0}\parallel \le {\tilde{\Delta}}_{p_0}^S\right)\wedge \left({D}_{\mathfrak{W},p}=0\right)\right\} $ the isotropic window of the surface filter, where $ {\mathfrak{P}}_{\mathfrak{f}} $ is the set of points of the source mesh and $ {D}_{\mathfrak{W},p} $ the shortest distance between $ p $ and the target walls $ \mathfrak{W} $ .

   1. 2. The volume filter ( $ \underset{\dot{\mkern6mu}}{-} $ ) is defined as

   (2) $$ \overline{\phi}\left({p}_0\right)=\sum \limits_{p\in {\mathfrak{w}}_{p_0}}{C}_{p_0}{V}_p{e}^{-\frac{{\boldsymbol{r}}_p-{\boldsymbol{r}}_{p_0}}{2{\left({\sigma}_{p_0}\right)}^2}}\phi (p), $$
   with $ {C}_{p_0} $ a normalization constant, $ {V}_p $ the nodal volume associated with node $ p $ on the source fine mesh, $ {\sigma}_{p_0}={\tilde{\Delta}}_{p_0}/\sqrt{12} $ the standard deviations of the volume Gaussian kernel, $ {\tilde{\Delta}}_{p_0}={\tilde{V}}_{p_0}^{1/3} $ the volume filter width, $ {\tilde{V}}_{p_0} $ the nodal volume of node $ {p}_0 $ on the target coarse mesh and $ {\mathfrak{w}}_{p_0}=\left\{p\in {\mathfrak{P}}_{\mathfrak{f}}:\left(\parallel {\boldsymbol{r}}_p-{\boldsymbol{r}}_{p_0}\parallel \le {\tilde{\Delta}}_{p_0}\right)\wedge \left(|{D}_{\mathfrak{W},p}-{D}_{\mathfrak{W},{p}_0}|\le \frac{1}{2}{\tilde{\Delta}}_{p_0}\right)\right\} $ the window of the volume filter, restricted in the wall-normal direction. The volume Gaussian filter is corrected to compensate any bias on the mean profiles induced by the filter. The corrected filtered field may be expressed as $ {\overline{\phi}}^{\ast }=\overline{\phi}+\left\langle \phi \right\rangle -\left\langle \overline{\phi}\right\rangle $ , where $ \left\langle \cdot \right\rangle $ denotes a time average.

This filtering operation verifies the properties of (1) the conservation of constants $ {\overline{a}}^{\ast }=a $ , for any constant $ a $ ; (2) linearity $ {\overline{\phi +\psi}}^{\ast }={\overline{\phi}}^{\ast }+{\overline{\psi}}^{\ast } $ , for any $ \phi $ and $ \psi $ ; and (3) DNS-convergence, $ {\lim}_{{\tilde{\Delta}}_{p_0}\to 0}\overline{\phi}=\phi $ , that is, the filter has no effect in the limit of an arbitrarily small filter size. It should, however, be noted that the filter does not commute with spatial derivation, since it is a discrete approximation of a truncated Gaussian filter (Sagaut, Reference Sagaut2006).

The filtered fields are partitioned to produce contiguous chunks of consistent size that can be used to train the model. The chunks only include nodes that are separated from the target walls $ \mathfrak{W} $ by $ {N}_H=3 $ edges or less. The partitioning is performed using the library METIS (Karypis and Kumar, Reference Karypis and Kumar1998).

3.1. Mach-number equivariance

First, we give some theoretical background on the behavior of low Mach number flows, that will be relevant to construct the input features. We will show how, under some assumptions, modifying the scales of velocity, density, viscosity, and length does not change the underlying flow physics and thus should accordingly not alter the behavior of the model.

Consider a flow modeled using the compressible Navier–Stokes equations without body forces or heat sources and the ideal gas equation of state,

(3) $$ \frac{\partial \rho }{\partial t}+\frac{\partial \rho {u}_j}{\partial {x}_j}=0, $$

(4) $$ \frac{\partial \rho {u}_i}{\partial t}+\frac{\partial \rho {u}_j{u}_i}{\partial {x}_j}=-\frac{\partial p}{\partial {x}_i}+\frac{\partial {\sigma}_{ij}}{\partial {x}_j}, $$

(5) $$ \frac{\partial p}{\partial t}+\frac{\partial {u}_jp}{\partial {x}_j}=-\left(\gamma -1\right)\frac{\partial {q}_j}{\partial {x}_j}+p\left(1-\gamma \right)\frac{\partial {u}_j}{\partial {x}_j}, $$

(6) $$ p= r\rho T, $$

where $ \rho $ is the density, $ t $ the time, $ p $ the pressure, $ \gamma $ the adiabatic index of the fluid, $ r $ is the ideal gas-specific constant, $ {u}_i $ the $ i $ th component of velocity, and $ {x}_i $ the Cartesian coordinate in $ i $ th direction. The shear-stress tensor $ \sigma $ and the conductive heat flux $ \boldsymbol{q} $ are assumed to be of the form $ {\sigma}_{ij}=\mu (T)\left[\left({\partial}_j{u}_i+{\partial}_i{u}_j\right)-\left(2/3\right){\partial}_k{u}_k{\delta}_{ij}\right] $ and $ {q}_j=-\lambda (T){\partial}_jT $ respectively, where the dynamic viscosity $ \mu (T) $ and the heat conductivity $ \lambda (T) $ are functions of temperature. Dissipation has been neglected in the pressure evolution equation (5), as this term vanishes in the low Mach number limit (Paolucci, Reference Paolucci1982). The dynamic viscosity $ \mu (T) $ is monotonous increasing in terms of $ T $ and separable, $ \mu (kT)={h}^{-1}(k)\mu (T) $ , which is for instance the case for a power law $ \mu (T)={\mu}_0{\left(T/{T}_0\right)}^k $ . The thermal conductivity is given by $ \lambda (T)=\mu (T){C}_p/\mathit{\Pr} $ , with $ {C}_p $ the isobaric heat capacity of the fluid and $ \mathit{\Pr} $ the Prandtl number of the fluid, both assumed constant. Suppose that this flow is characterized by a density scale $ {\rho}^b $ , a velocity scale $ {u}^b $ , and a pressure scale $ {p}^b $ . The nondimensional numbers associated with the flow are the Reynolds number $ \mathit{\operatorname{Re}}={\rho}^b{u}^b{x}^b/\mu \left({T}^b\right) $ , the Prandtl number $ \mathit{\Pr}=\mu \left({T}^b\right){C}_p/\lambda \left({T}^b\right) $ , and the Mach number $ Ma={u}^b/{c}^b $ , with $ {x}^b $ a length scale characterizing the geometry, and $ {c}^b=\sqrt{\gamma {rT}^b} $ the typical speed of sound.

At the low Mach number limit, modifying the scales of the flows does not modify the flow in a nondimensional sense, provided that the Reynolds and Prandtl number are kept constant. Indeed, consider another similar flow with a different Mach number and different scales of length, velocity, and temperature but the same Reynolds and Prandtl number, as follows:

$$ {\displaystyle \begin{array}{l}{x^{\prime}}^b=\beta {x}^b;\hskip11.6em {\rho^{\prime}}^b=\left(1/\alpha \right){\rho}^b;\\ {}{u^{\prime}}^b=\left(\alpha \varpi /\beta \right){u}^b;\hskip7em {T^{\prime}}^b=h\left(\varpi \right){T}^b;\\ {}{\mu^{\prime}}^b=\mu \left({T^{\prime}}^b\right)={\varpi \mu}^b;\hskip4.8em {\lambda^{\prime}}^b={\mu^{\prime}}^b{C}_p/P{r}^{\prime }={\varpi \lambda}^b;\\ {}R{e}^{\prime }=\mathit{\operatorname{Re}};\hskip11em Ma^{\prime }=\left(\alpha \varpi /\beta \right) Ma,\\ {}P{r}^{\prime }=\mathit{\Pr};\end{array}} $$

where $ \alpha $ , $ \beta $ , and $ \varpi $ are constant scalars characterizing the transformation. At low Mach number, the nondimensionalized density and velocity are independent of the Mach number while the Mach number dependence of pressure cannot be neglected. Namely, $ u/{u}^b\approx {\hat{u}}_0 $ , $ \rho /{\rho}^b\approx {\hat{\rho}}_0 $ and $ p/{p}^b\approx {\hat{p}}_0+{Ma}^2{\hat{p}}_1 $ , where $ {\hat{u}}_0 $ , $ {\hat{\rho}}_0 $ , $ {\hat{p}}_0 $ , and $ {\hat{p}}_1 $ do not depend on the Mach number (Lions and Moulden, Reference Lions and Moulden1996; Meister, Reference Meister1999; Munz et al., Reference Munz, Dumbser and Zucchini2003). The zeroth-order nondimensionalized pressure $ {\hat{p}}_0 $ can be shown to be constant in space by injecting these asymptotic developments into the Navier–Stokes equations. It is therefore useful to decompose pressure in a thermodynamical pressure $ {p}_0={p}^b{\hat{p}}_0 $ and a mechanical pressure $ {p}_1=p-{p}_0={p}^b{Ma}^2{\hat{p}}_1 $ . Hence, the thermodynamical and mechanical pressures are given by $ {p}_0^{\prime }={p^{\prime}}^b{\hat{p}}_0=\left(h\left(\varpi \right)/\alpha \right){p}_0 $ and $ {p}_1^{\prime }={p^{\prime}}^b{\alpha}^2{Ma}^2{\hat{p}}_1=\left({\alpha \varpi}^2/{\beta}^2\right){p}_1 $ .

Based on these insights, let us consider the change of variables $ x=\left(1/\beta \right){x}^{\prime } $ , $ t=\left(\alpha \varpi /{\beta}^2\right){t}^{\prime } $ , $ u={u}^{\prime}\beta /\left(\alpha \varpi \right) $ , $ \rho =\alpha {\rho}^{\prime } $ , $ {p}_0=\left(\alpha /h\left(\varpi \right)\right){p}_0^{\prime } $ , $ {p}_1={p}_1^{\prime }{\beta}^2/\left({\alpha \varpi}^2\right) $ , $ T=\left(1/h\left(\varpi \right)\right){T}^{\prime } $ , $ \mu ={\mu}^{\prime }/\varpi $ , and $ \lambda ={\lambda}^{\prime }/\varpi $ . Hereafter, we refer to this change of variable as a Mach-number transformation, since this change of variable modifies the Mach number of the flow. Introducing the change of variables in equations (3)–(6) leads to:

(7) $$ \frac{\partial {\rho}^{\prime }}{\partial {t}^{\prime }}+\frac{\partial {\rho}^{\prime }{u}_j^{\prime }}{\partial {x}_j^{\prime }}=0, $$

(8) $$ \frac{\partial {\rho}^{\prime }{u}_i^{\prime }}{\partial {t}^{\prime }}+\frac{\partial {\rho}^{\prime }{u}_j^{\prime }{u}_i^{\prime }}{\partial {x}_j^{\prime }}=-\frac{\partial {p}_{1^{\prime }}}{\partial {x}_i^{\prime }}+\frac{\partial {\sigma}_{ij}^{\prime }}{\partial {x}_j^{\prime }}, $$

(9) $$ \frac{\partial {p}_1^{\prime }}{\partial {t}^{\prime }}+\frac{\partial {u}_j^{\prime }{p}_1^{\prime }}{\partial {x}_j^{\prime }}=-\left(\gamma -1\right)\frac{\alpha^2{\varpi}^2}{h\left(\varpi \right){\beta}^2}\frac{\partial {q}_j}{\partial {x}_j^{\prime }}+\left({p}_1^{\prime}\left(1-\gamma \right)-\gamma {p}_0^{\prime}\frac{\alpha^2{\varpi}^2}{h\left(\varpi \right){\beta}^2}\right)\frac{\partial {u}_j^{\prime }}{\partial {x}_j^{\prime }}-\frac{\alpha^2{\varpi}^2}{h\left(\varpi \right){\beta}^2}\frac{\partial {p}_0^{\prime }}{\partial {t}^{\prime }}, $$

(10) $$ {p}_0^{\prime }+\frac{h\left(\varpi \right){\beta}^2}{\alpha^2{\varpi}^2}{p}_1^{\prime }=r{\rho}^{\prime }{T}^{\prime }. $$

At the limit of a low Mach number, the terms involving $ {p}_1^{\prime } $ in equations (9) and (10) vanish and the system of equations (7)–(10) tends to the low Mach number equations of Paolucci (Reference Paolucci1982):

(11) $$ \frac{\partial {\rho}^{\prime }}{\partial {t}^{\prime }}+\frac{\partial {\rho}^{\prime }{u}_j^{\prime }}{\partial {x}_j^{\prime }}=0, $$

(12) $$ \frac{\partial {\rho}^{\prime }{u}_i^{\prime }}{\partial {t}^{\prime }}+\frac{\partial {\rho}^{\prime }{u}_j^{\prime }{u}_i^{\prime }}{\partial {x}_j^{\prime }}=-\frac{\partial {p}_1^{\prime }}{\partial {x}_i^{\prime }}+\frac{\partial {\sigma}_{ij}^{\prime }}{\partial {x}_j^{\prime }}, $$

(13) $$ 0=-\left(\gamma -1\right)\frac{\partial {q}_j}{\partial {x}_j^{\prime }}-\gamma {p}_{0^{\prime }}\frac{\partial {u}_j^{\prime }}{\partial {x}_j^{\prime }}-\frac{\partial {p}_0^{\prime }}{\partial {t}^{\prime }}, $$

(14) $$ {p}_0^{\prime }=r{\rho}^{\prime }{T}^{\prime }. $$

This system of equations does not depend on any of the scalar parameters $ \alpha $ , $ \beta $ , or $ \varpi $ . Hence, the Mach-number transformation is parameterized by $ \alpha $ , $ \beta $ , and $ \varpi $ does not modify the system of equations. It follows that the solution of the system of equations (11)–(14) with different scales of velocity, density, viscosity, and length can be deduced from such transformation.

To summarize, a relevant corollary with respect to the machine learning model is that under a Mach-number transformation of its the input features,

(15) $$ {\displaystyle \begin{array}{ll}\overset{\check{}}{\boldsymbol{e}}=\boldsymbol{e}/\beta, & \overset{\check{}}{\rho }=\alpha \rho, \\ {}\overset{\check{}}{\boldsymbol{u}}=\boldsymbol{u}\beta /\left(\alpha \varpi \right),& \overset{\check{}}{T}=T/h\left(\varpi \right),\\ {}\overset{\check{}}{\mu }=\mu /\varpi, & \overset{\check{}}{\lambda }=\lambda /\varpi, \end{array}} $$

the model’s prediction should undergo the same transformation. Namely, the norm of the shear stress vector $ \tau =\Sigma \cdot {\boldsymbol{e}}_{\boldsymbol{n}}-\left({\boldsymbol{e}}_{\boldsymbol{n}}\cdot \Sigma \cdot {\boldsymbol{e}}_{\boldsymbol{n}}\right){\boldsymbol{e}}_{\boldsymbol{n}} $ , where $ \Sigma =\mu \left(\nabla \boldsymbol{u}+{\left(\nabla \boldsymbol{u}\right)}^T-\left(2/3\right)\left(\nabla \cdot \boldsymbol{u}\right){\mathbf{I}}_{\mathbf{d}}\right) $ is the viscous stress tensor and $ {\boldsymbol{e}}_{\boldsymbol{n}} $ a unit wall-normal vector, and the norm of the heat flux at the wall should undergo the transformation

(16) $$ \hskip0.5em \overset{\check{}}{\tau }=\alpha \tau {\left(\beta /\left(\alpha \varpi \right)\right)}^2,\hskip1.00em \overset{\check{}}{q}= q\beta /\left(\varpi h\left(\varpi \right)\right). $$

This property will thereupon will be referred to as the Mach-number equivariance of the model. Particular cases of the Mach-number transformation are the $ \alpha $ -transformation, $ \beta $ -transformation, and $ \varpi $ -transformation, which respectively correspond to $ \beta =\varpi =1 $ , $ \alpha =\varpi =1 $ , and $ \alpha =\beta =1 $ . The equivariance of the model under a Mach-number transformation is not exact for real flows, including the flows in the numerical database, as it relies on several approximations, most notably that of a low Mach number and constant adiabatic index. It is nevertheless useful to introduce this approximate equivariance in the learning process to build with limited data a model that can generalize to flows with different scales of velocity, density, viscosity, and length. The a priori results should demonstrate that these approximations are not critical for the performance of the model.

3.2. Strategies for equivariance enforcement

There are at least to two ways to leverage the equivariance of the model to a transformation $ \mathcal{P} $ within the training process: the invariant feature strategy and the data-augmentation strategy. In the invariant-feature strategy, the model is expressed in terms of invariant trainable functions by algorithmically constructing invariant features from the input data $ \boldsymbol{X} $ (Villar et al., Reference Villar, Storey-Fisher, Yao and Blum-Smith2021, Reference Villar, Yao, Hogg, Blum-Smith and Dumitrascu2022). One particular way to construct these invariant features is to use a reference to restrict the input space of the model. For instance, a privileged direction and an intrinsic velocity reference can be used to enforce rotational invariance and Galilean equivariance respectively. For machine-learning wall modeling, the wall may provide such references for the $ \alpha $ - and $ \varpi $ -transformations. In the data-augmentation strategy, the input space of the model is increased by using the transformation $ \mathcal{P} $ to create modified copies of the training samples. For instance, arbitrary rotations and translations can be used to enforce rotational invariance and Galilean equivariance respectively. The two strategies are represented schematically in Figure 1. The invariant-feature strategy is generally more costly than the data augmentation in terms of training time. However, it adds a pre-processing and post-processing step to the model at inference time, which may hinder its implementation in a CFD code.

Figure 1. Schematic representation of the scaling and data augmentation strategies, where the input space is formally decomposed into the dimension of a transformation $ \mathcal{P} $ , with respect to which the model is assumed equivariant, and an intrinsic dimension, which represents of the physical content of the input data, invariant under $ \mathcal{P} $ . The red crosses represent different simulations seen during training or at inference time. The blue areas are the portion of input space seen during training, which for illustrative purpose is assumed to be a sphere in the feature space.

In the present study, the equivariance of the model under an orthogonal transformation, a Galilean transformation, and the $ \alpha $ -, $ \beta $ -, and $ \varpi $ -transformations defined in Section 3.1 are considered. Since the walls are assumed to be isothermal, the dynamic viscosity is constant at the walls assuming that dynamic viscosity is only function of temperature. Accordingly, it is straightforward to define a $ \varpi $ -transformation that scales dynamic viscosity. In contrast, the wall density may vary due to variations of pressure and accordingly, it is difficult to equivocally select a particular value of the wall density to define an $ \alpha $ -transformation that scales density. Indeed, the density at each grid point is used for predictions at several wall locations. In light of the above, a strategy based on data augmentation has been selected to impose the equivariance of the model under a $ \varpi $ -transformation while a strategy based on the scaling of input features has been selected to impose the equivariance of the model under an $ \alpha $ -transformation. Data augmentation is also used to ensure that the prediction of the model does not depend on the coordinate system. Table 2 summarizes the strategy used to enforce the equivariance of the model under each transformation considered:

• • Orthogonal equivariance is ensured by augmenting the dataset with arbitrary three-dimensional rotations and reflections of the axes.

• • Galilean equivariance is ensured by expressing the input velocity relatively to the wall motion.

• • $ \alpha $ -equivariance is ensured by augmenting the dataset with arbitrary $ \alpha $ -transformations, with $ \alpha $ sampled from the log-uniform distribution with range $ \left[{\rho}_{\mathrm{min}}/{\left\langle \rho \right\rangle}_{\mathfrak{W}},{\rho}_{\mathrm{max}}/{\left\langle \rho \right\rangle}_{\mathfrak{W}}\right] $ , where $ {\left\langle \parallel \rho \parallel \right\rangle}_{\mathfrak{W}} $ is the mean wall density and where $ \left[{\rho}_{\mathrm{min}},{\rho}_{\mathrm{max}}\right] $ is the desired density range of the model.

• • $ \beta $ -equivariance is ensured by scaling the input data using the average edge length of the graph, denoted $ {\left\langle \parallel \boldsymbol{e}\parallel \right\rangle}_{\mathfrak{G}} $ , that is by applying a $ \beta $ -transformation with $ \beta ={\left\langle \parallel \boldsymbol{e}\parallel \right\rangle}_{\mathfrak{G}} $ .

• • $ \varpi $ -equivariance is ensured by scaling the input data using the wall dynamic viscosity $ {\mu}_{\mathrm{wall}} $ , that is by applying a $ \varpi $ -transformation with $ \varpi ={\mu}_{\mathrm{wall}} $ .

Table 2. Strategies used to enforce the equivariance of the machine-learning model to various transformations, and references used to compute the scalings

For ease of interpretation, notice that the scaled input velocity may be expressed as $ \overset{\check{}}{\boldsymbol{u}}={\overset{\check{}}{\boldsymbol{u}}}^{+}{\left\langle \parallel {\boldsymbol{e}}^{+}\parallel \right\rangle}_{\mathfrak{G}}/\left({\alpha \rho}_{\mathrm{wall}}\right) $ , the scaled target wall shear stress as $ \overset{\check{}}{\tau }={\left\langle \parallel {\boldsymbol{e}}^{+}\parallel \right\rangle}_{\mathfrak{G}}^2/\left({\alpha \rho}_{\mathrm{wall}}\right) $ and the scaled wall heat flux $ \overset{\check{}}{q} $ is proportional to $ Cp\left(\left(\overset{\check{}}{T}-{\overset{\check{}}{T}}_{\mathrm{wall}}\right)/{\overset{\check{}}{T}}_{\mathrm{wall}}\right){\left\langle \parallel {\boldsymbol{e}}^{+}\parallel \right\rangle}_{\mathfrak{G}}/{T}^{+} $ , where $ {}^{+} $ denotes the classical wall-unit scaling, namely $ {e}^{+}={eu}_{\tau }/{\nu}_{\omega } $ , $ {\boldsymbol{u}}^{+}=\boldsymbol{u}/{u}_{\tau } $ , and $ {T}^{+}={\rho}_{\omega }{C}_p{u}_{\tau}\left(T-{T}_{\mathrm{wall}}\right)/q $ . This scaling uniformizes the length scale of the various simulations included in the training database, but it does not negate differences in terms of mesh resolution, which implies that ideally a wide range of mesh resolution need to be included in the training dataset. Similarly, the scaling uniformizes the viscosity scale of the various simulations included in the training database, but it does not negate differences in terms of viscosity ratio between the wall and bulk flow, which implies that ideally a wide range of such viscosity ratio should be found in the training dataset.

4.1. Baseline wall models

Two algebraic wall models are used as baselines: an uncoupled velocity-temperature algebraic wall model and the coupled velocity-temperature algebraic wall model of Cabrit and Nicoud (Reference Cabrit and Nicoud2009).

The uncoupled algebraic wall model provides a boundary condition for the wall shear stress and the wall conductive heat flux without taking into account the coupling between the velocity and temperature profiles. Namely, the wall shear stress is obtained by solving the incompressible law of the wall,

(17) $$ \left\{\begin{array}{ll}{y}^{+}={u}^{+}& \mathrm{if}\;{y}^{+}<{y}_c^{+},\\ {}\left(1/\kappa \right)\log \left({y}^{+}\right)+C={u}^{+}& \mathrm{if}\;{y}^{+}\ge {y}_c^{+},\end{array}\right. $$

and the wall heat flux is obtained by solving the temperature law

(18) $$ \left\{\begin{array}{ll}{Pry}^{+}={T}^{+},& \mathrm{if}\;{y}^{+}<{y}_c^{+},\\ {}\left({\mathit{\Pr}}_t/\kappa \right)\log \left({y}^{+}\right)+{C}_T={T}^{+},& \mathrm{if}\;{y}^{+}\ge {y}_c^{+},\end{array}\right. $$

with $ \kappa =0.41 $ the von Kármán constant, $ C=5.5 $ a scalar constant, $ {\mathit{\Pr}}_t=0.85 $ the turbulent Prandtl number and

(19) $$ {C}_T={\left(3.85{\mathit{\Pr}}^{1/3}-1.3\right)}^2+2.12\log \left(\mathit{\Pr}\right) $$

a function of the Prandtl number (Kader, Reference Kader1981). The threshold $ {y}_c^{+}=11.445 $ is a critical value of the scaled first point height separating the viscous and inertial sublayers. The scaled wall distance is $ {y}^{+}={yu}_{\tau }/{\nu}_{\omega } $ , the scaled temperature $ {T}^{+}=\left(T-{T}_{\omega}\right)/{T}_{\tau } $ and the scaled velocity $ {u}^{+}=u/{u}_{\tau } $ . The wall shear stress and conductive heat flux are solved sequentially from equations (17) then (18) and the friction velocity and friction temperature definitions, namely $ \tau ={\rho}_{\omega }{u}_{\tau}^2 $ and $ {q}_{\omega }={\rho}_{\omega }{C}_p{u}_{\tau }{T}_{\tau } $ . The uncoupled algebraic model does not take into account the effect of the temperature variations to determine the wall shear stress. It is commonly used and accurate in flows with a negligible coupling between the fields of velocity and temperature, for example, when the temperature variations are not large.

The coupled algebraic wall model of Cabrit and Nicoud (Reference Cabrit and Nicoud2009) is an algebraic model that takes into account the coupling between the fields of velocity and temperature. In the case of a compressible non-reacting flow, the coupled wall model of Cabrit and Nicoud (Reference Cabrit and Nicoud2009) may be expressed by the following system of equations.

• • Velocity equation:

(20) (21) $$ \left\{\begin{array}{ll}{y}^{+}={u}^{+}& \mathrm{if}\;{y}^{+}<{y}_c^{+},\\ {}\frac{1}{K}\log \left({y}^{+}\right)+C=\frac{2}{{\mathit{\Pr}}_t{B}_q}\left(\sqrt{1-{KB}_q}-\sqrt{1-{T}^{+}{B}_q}\right)& \mathrm{if}\;{y}^{+}\ge {y}_c^{+};\end{array}\right. $$

• • Temperature equation:

(22) $$ {\mathit{\Pr}}_t{u}^{+}+K={T}^{+}; $$

with $ {B}_q={T}_{\tau }/{T}_{\omega } $ the anisothermicity factor, and $ K $ a function of the Prandtl number, given by

(23) $$ K={C}_T-{\mathit{\Pr}}_tC+\left(\frac{{\mathit{\Pr}}_t}{\kappa }-2.12\right)\left(1-2\log (20)\right), $$

where $ {C}_T $ is defined in equation (19). To solve this system of equations, a candidate friction velocity is first computed by injecting $ {y}^{+}={yu}_{\tau }/{\nu}_{\omega } $ and $ {T}^{+}=\left(T-{T}_{\omega}\right)/{T}_{\tau } $ in equation (2 1). This candidate value is discarded if the corresponding scaled first point height $ {y}^{+} $ is less than $ {y}_c^{+} $ , in which case $ {u}_{\tau } $ is solved using (20). In any case, the friction temperature $ {T}_{\tau } $ is then computed by injecting $ {u}^{+}=u/{u}_{\tau } $ in equation (22). The wall shear stress and conductive heat flux are deduced from the friction velocity and temperature definitions, namely $ \tau ={\rho}_{\omega }{u}_{\tau}^2 $ and $ {q}_{\omega }={\rho}_{\omega }{C}_p{u}_{\tau }{T}_{\tau } $ . The coupled algebraic model of Cabrit and Nicoud (Reference Cabrit and Nicoud2009) takes into account the effect of the temperature variations to determine the wall shear stress. It is relevant for fully developed turbulent boundary layers with large temperature variations.

4.2. Graph neural network wall model

A graph neural network is a class of artificial neural network designed to process graph-structured data, that is a set of nodes and edges along with, optionally, node-valued or edge-valued fields. In the present work, the “graph” corresponds to the numerical fields of an unstructured-mesh simulation. The nodes of the graph are given by the nodes of the mesh. The edges of the graph are given by the edges of the mesh, or in other words the mesh connectivity. The present graph neural network wall model uses the Encode-Process-Decode architecture of Battaglia et al. (Reference Battaglia, Hamrick, Bapst, Sanchez-Gonzalez, Zambaldi, Malinowski, Tacchetti, Raposo, Santoro, Faulkner, Gulcehre, Song, Ballard, Gilmer, Dahl, Vaswani, Allen, Nash, Langston, Dyer, Heess, Wierstra, Kohli, Botvinick, Vinyals, Li and Pascanu2018) and Dupuy et al. (Reference Dupuy, Odier, Lapeyre and Papadogiannis2023b), represented schematically in Figure 2. This architecture has been selected for its ability to directly operate on the unstructured data of the simulation, without requiring the prior interpolation of the inputs. The method is owing to this property easily applicable in LESs of complex geometries, with angles, corners, or curvature. Furthermore, the architecture implicitly makes some physical assumptions that are relevant for the modeling of the wall shear stress or wall heat flux. Namely, the prediction is invariant under a translation of the computational domain and biased toward locality, meaning that the model will more easily use nearby spatial locations for its prediction than far-away locations. The model will be trained to operate on mesh-based data, and produce a field of the wall shear stress or the wall conductive heat flux. To this end, it will leverage the inherent spatial relationships and dependencies found in the training database. Note that while the present approach does not incorporate specific physical laws or equations as in the case of physics-informed neural networks (PINN, (Raissi et al., Reference Raissi, Perdikaris and Karniadakis2019), it is not purely data-based since it is combined with scaling and data augmentation techniques presented in Section 3.2 to improve the generalizability of the model to other scales of length, velocity, temperature or pressure, and encode relevant physical equivariances. Due to this, the model will encode the physical assumptions of Galilean invariance, orthogonal equivariance that is the equivariance under a rotation or reflection of the computational domain, and Mach-number equivariance that is the equivariance under a change of scale that can be made under the assumption of a low Mach number. We accordingly assume that the effect of the Mach number on the boundary layer is not too large, which excludes hypersonic flows from the domain of applicability of the method. Finally, we assume that the Prandtl number of the flow is constant and that the variations of adiabatic index can be neglected for the temperature range of the flow.

Figure 2. Graphical representation of the Encode-Process-Decode architecture. The scaling (pre-processing) and unscaling (post-processing) steps correspond to the transformations as listed in Table 2.

The model architecture may be described as follows. The input of the model $ \boldsymbol{X} $ is first encoded into a latent space, processed iteratively by $ N $ identical message-passing steps, and finally decoded back to physical space to produce the output $ \boldsymbol{Y} $ :

(24) $$ \boldsymbol{Y}={f}_{\gamma}\left(\boldsymbol{X}\right)={f}_{\delta}\hskip0.35em \circ \hskip0.35em \underset{N\;\mathrm{times}}{\underbrace{f_{\pi}\hskip0.35em \circ \hskip0.35em \dots \hskip0.35em \circ \hskip0.35em {f}_{\pi}\hskip0.35em \circ \hskip0.35em {f}_{\pi }}}\hskip0.35em \circ \hskip0.35em {f}_{\varepsilon}\left(\boldsymbol{X}\right), $$

where $ {f}_{\varepsilon } $ , $ {f}_{\pi } $ , and $ {f}_{\delta } $ are learned functions associated respectively with the encoding, processing, and decoding steps. The model input $ \boldsymbol{X} $ is represented as a symmetric directed graph valued at both edges and nodes. Namely, it is an ordered pair $ \boldsymbol{X}=\left(\boldsymbol{V},\boldsymbol{E}\right) $ , where $ \boldsymbol{V} $ is the vector of the input features at the nodes of the graph and $ \boldsymbol{E} $ the vector of the input features at the edges of the graph. The encoder $ {f}_{\varepsilon } $ is a set of two multilayer perceptrons $ {f}_{\varepsilon}^V $ and $ {f}_{\varepsilon}^E $ that upscale the features at each edge or node independently. Specifically, the encoded feature $ {\hat{V}}_i^0 $ at node $ {v}_i $ is given by $ {\hat{V}}_i^0={f}_{\varepsilon}^V\left({V}_i\right) $ and the encoded feature $ {\hat{E}}_k^0 $ at edge $ {e}_k $ is given by $ {\hat{E}}_k^0={f}_{\varepsilon}^E\left({E}_k\right) $ , where $ {V}_i $ the input feature at node $ {v}_i $ and $ {E}_k $ the input feature at edge $ {e}_k $ . The processor is the core part of the prediction. It consists of the successive application of message-passing steps, which propagate the information through the edges and nodes of the graph. Namely, the message-passing function $ {f}_{\pi } $ is a set of two multilayer perceptrons $ {f}_{\pi}^V $ and $ {f}_{\pi}^E $ with residual connections that updates features at an edge from the nodes it contains and the features at a node from edges pointing to it. Specifically, the updated features at iteration $ n $ are given by

(25) $$ {\hat{V}}_i^n={\hat{V}}_i^{n-1}+{f}_{\pi}^V\left({\hat{V}}_i^{n-1},\sum \limits_{k\mid R(k)=i}{\hat{E}}_k^n\right), $$

(26) $$ {\hat{E}}_k^n={\hat{E}}_k^{n-1}+{f}_{\pi}^E\left({\hat{E}}_k^{n-1},{\hat{V}}_{S(k)}^{n-1},{\hat{V}}_{R(k)}^{n-1}\right), $$

where each edge $ {e}_k $ is assumed to point from node $ {v}_{S(k)} $ to node $ {v}_{R(k)} $ . The same message-passing operations are learned and applied everywhere in the graph. The decoder $ {f}_{\delta } $ is a multilayer perceptron that operates on node features only. The decoded output of the model at node $ {v}_i $ is given by $ {Y}_i={f}_{\delta}\left({\hat{V}}_i^N\right) $ . Following Dupuy et al. (Reference Dupuy, Odier, Lapeyre and Papadogiannis2023b), the size of the latent space is $ {d}_L=128 $ at both nodes and edges and the number of processing steps is $ N=4 $ in the present study. This value implies that edge and node features that are within a distance up to $ N=4 $ neighbors away from a given target wall node will influence the prediction. This value was found to be sufficient to discriminate non-equilibrium regions in Dupuy et al. (Reference Dupuy, Odier, Lapeyre and Papadogiannis2023b). The multilayer perceptrons $ {f}_{\varepsilon}^V $ , $ {f}_{\varepsilon}^E $ , $ {f}_{\pi}^V $ , and $ {f}_{\pi}^E $ are composed of $ {n}_{\mathrm{\ell}}=2 $ layers while the multilayer perceptron $ {f}_{\delta}^V $ is composed of $ {n}_{\mathrm{\ell}}+1 $ hidden layers. Each layer is composed of $ {d}_L $ neural units with a rectified linear unit $ \sigma (x)=\max \left(0,x\right) $ as activation function and is followed by layer normalization (Ba et al., Reference Ba, Kiros and Hinton2016), except for the last layer of the multilayer perceptron $ {f}_{\delta}^V $ which is linear. The weights of the multilayer perceptrons are learned using the Adam optimizer (Kingma and Ba, Reference Kingma and Ba2015) with a base learning rate of 0.001.

Two graph neural network wall models are trained in order to model the wall shear stress and the wall conductive heat flux respectively. In both cases, the following features are included in the input $ \mathbf{X} $ of the model:

• • At the nodes: the three components of the velocity vector, the density, the dynamic viscosity, and a target boundary mask $ \chi $ equal to 1 on the target wall nodes $ \mathfrak{W} $ and 0 otherwise;

• • At the edges: the three components of the edge displacement vector, and the length of the edge.

To provide the geometric information of the mesh to the model, we preferred to give the displacement vector of the edges, instead of the node coordinates, because it implies that the input graph will be the same if the geometry is translated, ensuring the translational invariance of the prediction. The angle or distance to the wall is not given to the model as an input. However, it is possible for the model to reconstruct the distance to the wall of an information at a given node, based on the displacement vectors of the various edges connecting the point to the wall. The rationale is to alleviate the need to manually define robust geometrical metrics that can adapt to any complex geometry, and instead opts to let the network learn the optimal way to extract the necessary information itself. The loss function is the mean squared error (MSE) between the output of the model and the reference wall shear stress or wall conductive heat flux, that is either

(27) $$ {\mathrm{\mathcal{L}}}_{\tau }={\left\langle \chi {\left({\overset{\check{}}{\overline{\tau}}}_{\mathrm{ref}}^S-Y\right)}^2\right\rangle}_{\mathfrak{A}},\mathrm{or} $$

(28) $$ {\mathrm{\mathcal{L}}}_q={\left\langle \chi {\left({\overset{\check{}}{\overline{q}}}_{\mathrm{ref}}^S-Y\right)}^2\right\rangle}_{\mathfrak{A}}, $$

with $ {\left\langle \cdot \right\rangle}_{\mathfrak{A}} $ an average across graph nodes and samples. All the input and reference output variables are scaled and/or augmented as described in Section 3.2. To minimize the loss, the learned function has to be resilient to a change of the mesh. Indeed, the training databases prepared in Section 2.1 aggregates various meshes with different resolutions for a given simulation, and the different geometries associated with the different configurations. Furthermore, there are variations present within a given tetrahedral mesh as well. An accurate prediction of the wall shear stress or wall conductive heat flux thus entails being robust to mesh modifications. This will be confirmed in a priori tests.

6.1. Numerical method

In all simulations, the flow is modeled as a continuous medium in local thermodynamic equilibrium using the compressible LES equations and the ideal gas law. These equations may be expressed as follows:

(33) $$ {\partial}_t\rho +\nabla \cdot \left(\rho \boldsymbol{U}\right)=0, $$

(34) $$ {\partial}_t\left(\rho \boldsymbol{U}\right)+\nabla \cdot \left(\rho \boldsymbol{U}\otimes \boldsymbol{U}\right)=-\nabla P+\nabla \cdot \mathtt{S}, $$

(35) $$ {\partial}_t\left(\rho E\right)+\nabla \cdot \left(\rho \boldsymbol{U}H\right)=-\nabla \cdot \boldsymbol{Q}-\nabla \cdot \left(\boldsymbol{U}P\right)+\nabla \cdot \left(\boldsymbol{U}\cdot \mathtt{S}\right), $$

(36) $$ P= r\rho T, $$

where $ \rho $ is the density of the fluid, $ t $ is the time, $ \boldsymbol{U} $ is the velocity, $ \boldsymbol{x} $ is the Cartesian coordinate, $ P $ is the pressure, $ \boldsymbol{Q} $ is the conductive heat flux, $ E $ is the total energy per unit mass, $ H=E+P/\rho $ is the total enthalpy per unit mass, and $ r $ is the specific ideal gas constant. The temperature is computed using tabulated data from Stull and Prophet (Reference Stull and Prophet1971) based on the internal energy $ e=E-\frac{1}{2}{U}_i{U}_i $ . No external body forces or volumetric heat sources are taken into account. The viscous stress tensor and conductive heat flux are computed assuming a Newtonian fluid under Stokes’ hypothesis and Fourier’s law, while subgrid-scale stresses are modeled using an eddy-viscosity model. The stress tensor $ \Sigma $ is given by:

(37) $$ {\Sigma}_{ij}=\left(\mu +{\mu}_{\mathrm{sgs}}\right)\left(\frac{\partial {U}_i}{\partial {x}_j}+\frac{\partial {U}_j}{\partial {x}_i}-\frac{2}{3}\frac{\partial {U}_k}{\partial {x}_k}{\delta}_{ij}\right), $$

where $ \mu $ is the dynamic viscosity of the fluid, $ {\mu}_{\mathrm{sgs}} $ the subgrid-scale viscosity, and $ {\delta}_{ij} $ denotes the Kronecker delta. The conductive heat flux is computed assuming Fourier’s law and an eddy-diffusivity subgrid-scale contribution:

(38) $$ \boldsymbol{Q}=-\left(\lambda +{\lambda}_{\mathrm{sgs}}\right)\nabla T, $$

where $ T $ is the fluid temperature and $ \lambda $ the thermal conductivity of the fluid, computed using a constant Prandtl number assumption. The subgrid-scale conductivity $ {\lambda}_{\mathrm{sgs}}={C}_p{\mu}_{\mathrm{sgs}}/{\mathit{\Pr}}_{\mathrm{sgs}} $ is computed assuming a constant subgrid-scale Prandtl number $ {\mathit{\Pr}}_{\mathrm{sgs}} $ . The isobaric heat capacity $ {C}_p $ is computed from tabulated data (Stull and Prophet, Reference Stull and Prophet1971).

The numerical simulations are performed by coupling the numerical resolution of equations (33)–(36) via the compressible, unstructured, massively parallel flow solver AVBP (Schönfeld and Rudgyard, Reference Schönfeld and Rudgyard1999), to the graph neural network wall model inference. The coupling, implemented via a message passing interface (MPI), is described in Dupuy et al. (Reference Dupuy, Odier, Lapeyre and Papadogiannis2023b). At each time step, data is extracted from the LES fields to create near-wall graphs and sent to the graph neural network for inference. The inference is executed exclusively on graphics processing units (GPUs), concurrently to parts of the LES schemes ran on CPUs. The calculation is distributed across multiple GPUs using a specific GPU partitioning that is separate from the CPU partitioning. To make sure that the model prediction is not influenced by the partitioning, the GPU partitions are overlapped by the number of message-passing steps $ N $ of the model. This overlap ensures that each partition includes the receptive field required for the prediction of the wall shear stress and wall heat flux. The predictions are then sent back to the main flow solver to set the boundary conditions in the momentum and energy transport equations.

6.2. Symmetrically cooled channel flows

This section addresses the wall-modeled LES of symmetrically cooled channel flows. We first consider the simulation SC2, included in the training database. The operating conditions of the wall-modeled simulations are the same as in the wall-resolved simulation (Appendix A), namely, the temperature ratio between the bulk flow and the walls is 3 and the friction Reynolds number is $ {\mathit{\operatorname{Re}}}_{\tau }=1150 $ . In the wall-modeled LESs, the size of the domain is $ 14h\times 2h\times 5.2h $ , with $ h $ the half-height. The WMLES mesh contains 8 million tetrahedral cells. The average height of the first point off the wall, in wall units, is $ {y}^{+}=40 $ . As in the wall-resolved simulation, the convective scheme is a two-step Taylor–Galerkin scheme with third-order spatial and temporal accuracy (Colin and Rudgyard, Reference Colin and Rudgyard2000) and the diffusive scheme is a centered second-order scheme. The subgrid-scale viscosity model is the Smagorinsky model (Smagorinsky, Reference Smagorinsky1963).

Wall-modeled LESs using the baseline uncoupled and coupled algebraic models, the ACSC and Full graph neural network wall models are compared in Figure 8. Namely, the profiles of mean streamwise velocity, mean temperature, standard deviation of streamwise velocity, and standard deviation of temperature are provided using the classical wall scaling ( $ + $ ). The uncoupled algebraic model leads to a mean temperature profile that is farther from the reference mean temperature profile than the other investigated models, since it neglects the effect of the large temperature variations. The uncoupled algebraic model and the two graph neural network wall models lead to profiles that are close to one another. The profiles obtained using the Full model are slightly closer to the reference wall-resolved profiles of temperature and velocity. With the Full model, the mean temperature profile is very close to the reference profile above $ {y}^{+}=200 $ while the slope of the profile is overestimated closer to the wall. The mean velocity profile is shifted upward compared to the reference profile in the logarithmic region. This behavior is alike the logarithmic layer mismatch that is classically obtained in wall-modeled LES, due to numerical and subgrid-modeling errors in the first few grid points (Kawai and Larsson, Reference Kawai and Larsson2012; Larsson et al., Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016; Bose and Park, Reference Bose and Park2018). This is confirmed by the simulation of an incompressible channel in the absence of thermal effects (Appendix B). The magnitude of the peak standard deviation of streamwise velocity is well predicted by all wall-modeled LESs (Figure 8), but its height is overestimated compared to the reference wall-resolved simulation. The peak standard deviation of temperature is also farther from the wall than in the reference, and its magnitude is underestimated by around 10%. Both the standard deviation of streamwise velocity and temperature are overestimated in the bulk of the channel in all wall-modeled LESs. For all turbulence statistics investigated, the profiles obtained with the two graph neural network wall models and with the baseline coupled algebraic model are very similar. To assess the influence of the subgrid-scale model on this conclusion, we performed simulations of the cooled channel SC2 using the Sigma subgrid-scale model (Nicoud et al., Reference Nicoud, Baya Toda, Cabrit, Bose and Lee2011). The simulations with the Sigma model exhibit a more pronounced logarithmic-layer mismatch compared to the simulations with the Smagorinsky model (Figure 9). This is consistent with studies from the literature, which have found that, in wall-modeled simulations, the Smagorinsky model provides a more accurate wall shear stress prediction than the Sigma model in a channel flow (Blanchard et al., Reference Blanchard, Odier, Gicquel, Cuenot and Nicoud2021). The profiles obtained with the Full graph neural network wall model and the baseline coupled algebraic model are similar to the Smagorinsky model. This confirms that, for the simulation of the symmetrically cooled channel flow SC2 included in the training database, the graph neural network wall models can operate in an a posteriori configuration and reach parity performance with a state-of-the-art algebraic model.

Figure 8. Mean streamwise velocity, standard deviation of streamwise velocity, mean temperature, and standard deviation of temperature in large-eddy simulations of the symmetrically cooled channel flow SC2 with algebraic wall models and graph neural network wall models. The wall-resolved simulation presented in Appendix A is given for comparison.

Figure 9. Mean streamwise velocity, standard deviation of streamwise velocity, mean temperature, and standard deviation of temperature in large-eddy simulations of the symmetrically cooled channel flow SC2 with algebraic wall models and graph neural network wall models, using the Smagorinsky and Sigma subgrid-scale models. The wall-resolved simulation presented in Appendix A is given for comparison.

The graph neural network wall models are now assessed for the wall-modeled LES of the symmetrically cooled channel flow SC3, not included in the training database. As described in Appendix A, the flow SC3 is similar to the flow SC2 but uses air as working fluid. In particular, the dynamic viscosity law, the ideal gas specific constant, and the scales of velocity, density, and pressure are modified compared to the case SC2. The mass flow rate is adjusted in order to preserve a similar mean friction Reynolds number to the simulation SC2. The computation domain, the mesh, and the numerical method used in the wall-modeled LESs are identical to those used for the simulation SC2. The numerical results are given in Figure 10. Similar to the SC2 simulation, the uncoupled algebraic model leads to a more inaccurate mean temperature profile than the coupled algebraic model, which is better suited to flows with large temperature variations. The Full model and the coupled algebraic model lead to very similar profiles for all turbulence statistics investigated. The mean velocity profile is as in the SC2 case overestimated in the logarithmic region but, as opposed to the SC2 case, the mean temperature profile is also overestimated in this region. Moreover, the standard deviation of both streamwise velocity and temperature is overestimated in the center of the channel but as in the SC2 simulation its peak magnitude is reasonably well predicted (less than 10% error). The wall-modeled LES with the ACSC model lead to a less accurate mean velocity profile than the simulation with the Full model. Overall, the Full model seems able to generalize to a symmetrically cooled channel flow, even for a different working fluid than seen during training. To provide a more quantitative assessment of the accuracy of the wall models, Table 4 reports the mean wall shear stress and conductive heat flux of each wall-modeled LES. The Full graph neural network wall model is slightly more accurate than the coupled algebraic model for the prediction of the mean wall shear stress in the SC2 and SC3 simulations, but slightly less accurate for the prediction of the wall conductive heat flux in the SC3 simulation.

Figure 10. Mean streamwise velocity, standard deviation of streamwise velocity, mean temperature, and standard deviation of temperature in large-eddy simulations of the symmetrically cooled channel flow SC3 with algebraic wall models and graph neural network wall models. The wall-resolved simulation presented in Appendix A is given for comparison.

Table 4. Mean nondimensionalized wall shear stress $ \tau /\left(\rho {u}_c^2\right)={\left({u}_{\tau }/{u}_c\right)}^2 $ , where $ {u}_c $ is the centerline velocity, in large-eddy simulations of a channel flow at the friction Reynolds numbers $ {\mathit{\operatorname{Re}}}_{\tau }=395 $ , $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ , and $ {\mathit{\operatorname{Re}}}_{\tau }=2000 $ with an algebraic wall stress model and a machine-learning wall model

Overall, these results demonstrate that the Full model is able to reproduce a level of accuracy that is similar to that of the coupled algebraic model for the simulation of a cooled channel with large temperature variations. This confirms the ability of the graph neural network wall modeling approach to generalize from a priori training to an a posteriori testing configuration for this particular flow. In the next section, we will assess to which extent the graph neural network wall modeling approach can provide an improvement for a more complex simulation that is known to be difficult to simulate using an algebraic model devised for anisothermal fully developed turbulent boundary layers.

6.3. High-pressure turbine blade

This section addresses the wall-modeled LES of the high-pressure turbine blade L89, not included in the training dataset. The L89 blade configuration is an experimental configuration instrumented for boundary-layer measurements with realistic operating conditions in terms of Reynolds and Mach numbers. It is thus in many ways representative of real turbomachinery flow, for which similar measurements are typically not available. The simulation reproduces the test case MUR235 of Arts et al. (Reference Arts, Rouvroit and Rutherford1990). This test case features a complex physics compared to the symmetrically cooled channel flow, that involves transitional and fully turbulent boundary layers with large temperature variations. In the literature, the LES of the test case MUR235, characterized by a large inlet turbulence intensity of 6%, has been found to be particularly challenging compared to other test cases of the VKI LS1989-MUR series (Arts et al., Reference Arts, Rouvroit and Rutherford1990). For instance, an accurate prediction is the wall heat flux on the suction side of the blade is more easily obtained for the MUR129 than the MUR235 test case (Collado Morata et al., Reference Collado Morata, Gourdain, Duchaine and Gicquel2012; Gourdain et al., Reference Gourdain, Gicquel and Collado Morata2012; Segui et al., Reference Segui, Gicquel, Duchaine and de Laborderie2017; Dupuy et al., Reference Dupuy, Gicquel, Odier, Duchaine and Arts2020. Accurate predictions have not been obtained with the resolution of a wall-modeled LES in the literature. This suggests the present MUR235 test case is relevant to assess the performance of wall models in a flow with complex boundary layer dynamics, and can be used as a benchmark for the development of wall models.

The operating conditions of the simulation are reported in Table 5. The flow is representative of a typical high-pressure turbine in terms of Mach and Reynolds numbers. Although the inlet temperature is low compared to the typical temperature downstream of a combustion chamber, the temperature ratio between the inlet flow and the blade is realistic. The blade profile is given by the reinterpolated manufacturing coordinates of Wheeler et al. (Reference Wheeler, Sandberg, Sandham, Pichler, Michelassi and Laskowski2016) to ensure a smooth surface curvature. The chord $ c $ is $ 67.647 $ mm. The computational domain is periodic in the spanwise and pitchwise directions to model an infinite linear cascade of two-dimensional linearly extruded blade profiles. The origin of the Cartesian coordinate system ( $ x $ , $ y $ , $ z $ ) is placed at the leading edge of the blade, where $ x $ is the axial coordinate, $ y $ the pitchwise coordinate, and $ z $ the spanwise coordinate. The curvilinear coordinate $ {s}_1 $ follows the blade profile and is zero at the leading edge, positive on the suction side of the blade, and negative on the pressure side of the blade (Figure 12). The inlet is located at $ {x}_{\mathrm{inlet}}=-0.81c $ and the outlet at $ {x}_{\mathrm{outlet}}=1.48c $ . The spanwise length of the domain is set $ 0.148c $ as this value is consistent with the literature (Bhaskaran and Lele, Reference Bhaskaran and Lele2010; Fransen et al., Reference Fransen, Collado Morata, Duchaine, Gourdain, Gicquel, Vial and Bonneau2011; Pichler et al., Reference Pichler, Kopriva, Laskowski, Michelassi and Sandberg2016) and the sensitivity analysis of Collado Morata et al. (Reference Collado Morata, Gourdain, Duchaine and Gicquel2012). A wall-resolved simulation of the case on a mesh containing 587 million cells, performed by Dupuy et al. (Reference Dupuy, Gicquel, Odier, Duchaine and Arts2020), is used as a reference. In addition, the experimental data of Arts et al. (Reference Arts, Rouvroit and Rutherford1990) is used when available. The WMLES mesh contains 47 million tetrahedral cells in total. The mesh is more refined close to the blade surface and near the trailing edge on the suction side (Figure 11, left). The distribution of the height of the first cell off the wall along the blade surface is given in Figure 11 (right). Convective and diffusion are both discretized using second-order accurate schemes (Lax and Wendroff, Reference Lax and Wendroff1960).

Table 5. Operating point conditions of the MUR235 test cases (Arts et al., Reference Arts, Rouvroit and Rutherford1990)

Figure 11. Left: Side view of the mesh, on the slice $ z/c=0.037 $ . Right: Height of the first cell off the wall along the blade surface, in wall units.

Wall-modeled LESs using the baseline uncoupled and coupled algebraic models, the ACSC and Full graph neural network wall models are performed using the Sigma (Nicoud et al., Reference Nicoud, Baya Toda, Cabrit, Bose and Lee2011) subgrid-scale viscosity model. The intensity of the incident turbulence in the wall-modeled LES is comparable to that of the reference high-fidelity simulation. Indeed, the turbulence intensity of the resolved scales, $ Tu $ , defined as

(39) $$ Tu=\frac{\sqrt{\frac{1}{3}\left(\left\langle {u^{\prime}}_x^2\right\rangle +\left\langle {u^{\prime}}_y^2\right\rangle +\left\langle {u^{\prime}}_z^2\right\rangle \right)}}{\sqrt{{\left\langle {u}_x\right\rangle}^2+{\left\langle {u}_y\right\rangle}^2+{\left\langle {u}_z\right\rangle}^2}}, $$

is 5% near the leading edge in the wall-modeled simulations, compared to 5.6% in the reference high-fidelity simulation. The instantaneous visualizations given in Figure 12 (center) show that the wall-modeled LES reproduces the main flow phenomena associated with the case. The trailing-edge flow is associated with both vortex shedding and the generation of a series of acoustic waves impinging on the suction side of the blade. Furthermore, there is a series of shocks quasi-normal to the blade on the suction side. The shocks induce a rapid transition of the suction-side boundary layer from a transitional state to a fully turbulent state (Dupuy et al., Reference Dupuy, Gicquel, Odier, Duchaine and Arts2020. This rapid turbulence transition is associated with a large increase of the wall shear stress in the region downstream of the shocks ( $ {s}_1/c>0.96 $ ). Figure 13a compares the distribution of the wall shear stress along the blade surface in the different simulations. The baseline uncoupled algebraic model, coupled algebraic model, and the ACSC graph neural network wall model lead to very similar wall shear stress profiles. In the three simulations, the wall shear stress is overestimated on the pressure side near the trailing edge and on the suction side until the shocks. In addition, the increase in wall shear stress associated with the shocks occurs further downstream in all wall-modeled LESs than in the reference high-fidelity simulation, even though the shocks location is accurately predicted. The simulation with the Full graph neural network wall model leads to a more accurate wall shear stress distribution. In particular, the wall shear stress is accurate throughout the pressure side of the blade and the peak wall shear stress after the strongly accelerating region near the leading edge of the blade ( $ {s}_1/c=0.3 $ ) is not overestimated. These findings are affected by the selection of the subgrid-scale model. Indeed, the quasi-normal shocks on suction side of the blade do not provoke the expected turbulent transition of the boundary layer using the Smagorinsky model. This halts the generation of acoustic waves at the trailing edge of the blade (Figure 12, bottom). Accordingly, the rapid increase in wall shear stress near the trailing edge of the blade ( $ {s}_1/c>0.96 $ ) is not reproduced with either the baseline algebraic model or the graph neural network wall model with the Smagorinsky model (Figure 14a). On the rest of the blade, the predictions are relatively similar to those obtained using the Sigma model. Namely, the Full graph neural network wall model predicts more accurately the wall shear stress on the pressure side compared to the coupled algebraic model, and the peak wall shear stress after the strongly accelerating region near $ {s}_1/c=0.3 $ on the suction side. Thus, the Full model is able to leverage the non-equilibrium configurations of the training database to improve the wall shear stress prediction in these regions, which are not fully turbulent.

Figure 12. Instantaneous field of axial velocity in the reference high-fidelity simulation of Dupuy et al. (Reference Dupuy, Gicquel, Odier, Duchaine and Arts2020) (top), the wall-modeled large-eddy simulation using the Sigma model (center), and the wall-modeled large-eddy simulation using the Smagorinsky model (bottom) on the plane $ z/c=0.037 $ . The curvilinear abscissa $ {s}_1/c $ along the blade is given in the figures for reference.

Figure 13. Wall shear stress and conductive heat flux along the blade surface using the Sigma model.

Figure 14. Wall shear stress and conductive heat flux along the blade surface, using the Smagorinsky and Sigma subgrid-scale models.

With regard to the wall heat flux prediction, Figures 13b and 14b, the Full model also improves the prediction compared to the baseline uncoupled algebraic model, coupled algebraic model, and the ACSC graph neural network wall model. With both the Sigma and Smagorinsky models, the peak wall heat flux location correctly occurs at the leading edge with the Full model, whereas it occurs further downstream, below $ {s}_1/c=0.2 $ , in the other wall-modeled LESs. In addition, the wall heat flux is significantly less overestimated on the suction side of the blade, in the transitional region before the shocks $ {s}_1/c<0.9 $ . The wall heat flux prediction is also improved on the pressure side of the blade, below $ {s}_1/c<-0.5 $ , especially using the Sigma subgrid-scale model. However, the simulation with the Full model leads to a large underestimation of the wall heat flux downstream of the shocks on the suction side, using the Sigma model, whereas the other wall-modeled LESs provide a more accurate prediction in that region. The deficiency is related to the incompressible isothermal simulation included in the training dataset for the Full model, since the ACSC model lead to a more accurate prediction in that region (Figure 13b). Overall, the results are nevertheless promising for the prospect of graph neural network wall modeling, since outside of this region the wall heat flux distribution on the blade surface is fairly well predicted compared to the reference simulation, providing a level of accuracy that is not trivial to reproduce even in a wall-resolved simulation (Collado Morata et al., Reference Collado Morata, Gourdain, Duchaine and Gicquel2012; Gourdain et al., Reference Gourdain, Gicquel and Collado Morata2012; Segui et al., Reference Segui, Gicquel, Duchaine and de Laborderie2017; Dupuy et al., Reference Dupuy, Gicquel, Odier, Duchaine and Arts2020). Machine-learning wall modeling provides a framework that has the potential to take into account both non-equilibrium effects and variable-density effects in complex geometries, and this ability could be progessively refined and improved over time following the developments of larger numerical training databases.

For completeness, it is useful to discuss the computational cost associated with such machine-learning wall models. The present model has not been developed with computational cost as a primary concern. For instance, the wall shear stress and the wall conductive heat flux are presently predicted by sequentially applying two independent neural networks, although they certainly share common features or representations in the hidden layers. Moreover, the hyperparameters of each network have been selected without considering computational cost and, as a result, we did not attempt to prune the networks to reduce computational overhead. The present graph neural network wall model has been benchmarked, for the wall-modeled simulation of the LS89 blade. The tests are performed on the Jean–Zay super-computing cluster, which as many other recent super-computing clusters has a hybrid architecture. The hybrid nodes contain 40 CPU cores (Intel Cascade Lake 6248 2.5 GHz) and 4 GPUs (Nvidia Tesla V100 SXM2 32GB). The wall-modeled simulations exploit these hybrid resources by executing the graph neural network prediction exclusively on the GPUS, concurrently to parts of the LES schemes, as discussed in Serhani et al. (Reference Serhani, Xing, Dupuy, Lapeyre and Staffelbach2022). In that context, we may defined the deep-learning overhead as the additional cost induced by the use of the graph neural network wall model on a given hybrid node, which may be calculated by considering the cost of the graph neural network prediction and related to communications, from which is subtracted the cost of overlapped LES computations. Figure 15 gives the overhead of the model for 2–16 hybrid nodes. The deep learning overhead ranges from 25 to 30% of the total iteration cost. Moving forward, further developments should consider computational cost in the model development.

Figure 15. Overhead of the graph neural network wall model per iteration for 2, 4, 6, and 16 hybrid nodes.