Modeling the wall shear stress in large-eddy simulation using graph neural networks

Abstract

As the Reynolds number increases, the large-eddy simulation (LES) of complex flows becomes increasingly intractable because near-wall turbulent structures become increasingly small. Wall modeling reduces the computational requirements of LES by enabling the use of coarser cells at the walls. This paper presents a machine-learning methodology to develop data-driven wall-shear-stress models that can directly operate, a posteriori, on the unstructured grid of the simulation. The model architecture is based on graph neural networks. The model is trained on a database which includes fully developed boundary layers, adverse pressure gradients, separated boundary layers, and laminar–turbulent transition. The relevance of the trained model is verified a posteriori for the simulation of a channel flow, a backward-facing step and a linear blade cascade.

Impact Statement

Numerical simulation is an important tool for design optimization in various industries, notably in the transport and aeronautics market. The development of graph-based machine-learning wall-shear-stress models for large-eddy simulation could improve the accuracy of practical industrial simulations, paving the way for a reduction in the energy consumption of aircraft, cars, and ships, resulting in lower emissions and noise.

1. Introduction

This paper investigates the modeling of the wall shear stress (WSS) in large-eddy simulations (LESs) using a graph neural network (GNN) that directly operates on the mesh of the simulation. In LESs of turbulent flows, wall modeling enables the use of coarser meshes, that are not able to resolve the velocity gradients in the inner viscous and buffer regions of the near-wall boundary layers. This ability relieves the grid-point and time-step requirements of LES as the Reynolds number increases (Yang and Griffin, 2021), and is key to the LES of complex industrial flows. Indeed, the use of a no-slip boundary condition requires an accurate resolution of the entire boundary layer, including the viscous sublayer. While LESs do not require the direct computation of the smaller turbulence length scales, even the most energetic integral length scales decrease in size near the walls hence preventing further coarsening in the wall-normal direction and constraining the time step of the simulation. This is often intractable in complex flows and in flows at high Reynolds number, as the number of grid points required to resolve the near-wall region becomes prohibitive. Wall models have been reviewed in Larsson et al. (2016) and Bose and Park (2018). The most common and classical approach is the algebraic equilibrium wall model, which imposes the WSS as a Neumann-type boundary condition by enforcing the law of the wall locally (Deardorff, 1970; Schumann, 1975; Piomelli et al., 1989). This assumption is relevant in fully developed boundary layers, for which the law of the wall holds in the mean, but cannot generally be assumed to hold in complex flows, notably in the presence of strong adverse pressure gradients, curvature, separated, or transitional boundary layers.

Several alternative wall models have been proposed in the literature. Bose and Moin (2014) and Bae et al. (2019) suggested a slip wall model that provides a boundary condition for the wall-parallel velocity. Chung and Pullin (2009), Inoue and Pullin (2011), Saito et al. (2012), and Gao et al. (2019) have coupled a stretched-vortex subgrid-scale model with a virtual wall model that terminates the LES domain at a finite distance from the wall. A large class of methods models the WSS using a strategy that combines LES and Reynolds-average Navier–Stokes (RANS) simulation. This includes zonal methods that solve the RANS or thin-boundary-layer equations (Cabot, 1995; Balaras et al., 1996; Baurle et al., 2003; Davidson and Peng, 2003; Temmerman et al., 2005; Piomelli, 2008). This is computationally expensive since the grid resolution requirements are not necessarily lower than would be in a wall-resolved LES, although the boundary-layer equations are less expensive to resolve. The integral wall model introduced by Yang et al. (2015) and extended by Catchirayer et al. (2018) use a parameterized velocity profile to avoid the numerical integration of the boundary-layer equations. Wall models based on control theory have also been suggested, in which the WSS is imposed so as to lead to the target mean velocity profile, either given by the law of the law or a coupled RANS model (Nicoud et al., 2001; Templeton et al., 2006; 2008). Recently, the use machine-learning methodologies to model the WSS has finally been proposed by several authors. Yang et al. (2019) developed a data-driven wall-stress model in a turbulent channel flow. Huang et al. (2019) compared data-driven and physics-based approaches in a spanwise rotating channel. Zhou et al. (2021) assessed a data-driven wall model for the turbulent flow over periodic hills. Lozano-Durán and Bae (2020) applied a self-critical machine-learning wall model to a wing-fuselage juncture flow, Zangeneh (2021) developed a data-driven wall model for supersonic separated flows. Bhaskaran et al. (2021) has trained a wall model for a boundary layer in the presence of shock-boundary layer interaction. Bae and Koumoutsakos (2022) assessed a multi-agent reinforcement learning approach to wall modeling in a turbulent channel flow. Radhakrishnan et al. (2021) developed a wall model based on gradient boosted decision trees. Dupuy et al. (2022) investigated the data-driven wall modeling of turbulent separated flows. These efforts may be seen as part of a broader movement to introduce data-driven methods to the field of turbulence simulation and modeling (Duraisamy et al., 2019; Brunton et al., 2020). This has been motivated by the recent emergence of data-driven approaches with capability beyond those of handcrafted models in various complex tasks, such as image classification, speech recognition, natural language processing, short-term weather prediction, or protein structure prediction (LeCun et al., 2015; Jumper et al., 2021; Ravuri et al., 2021). Supervised machine-learning wall-shear-stress models have proven capable of discriminating and operating in non-equilibrium flow regions, provided that similar flow phenomena were encountered in the training dataset (Dupuy et al., 2022).

The studies on machine-learning wall-shear-stress modeling have demonstrated the advantages of going beyond algebraic formulations and capturing topological flow features to improve the wall friction prediction. However, a challenge with nonlocal models is defining this topological extraction in a general manner that applies to any complex geometry. As an example, multipoint stencils might exit the flow domain in small nooks or sharp concave corners. The present study proposes a general formulation that is adaptable to any mesh, without any concern for how to implement it geometrically, by exploiting GNNs that rely on the data mesh itself. This work directly allows the use of data on unstructured meshes similarly to Lozano-Durán and Bae (2020), which is the first study to apply a machine-learning wall-shear-stress model in an unstructured-grid simulation. In their study, the input data from the flow simulation are interpolated on a local Cartesian grid to feed the inputs of the wall-shear-stress model, arranged in a predefined stellated arrangement with uniform spacing. By contrast, the present machine-learning model is a topologically aware wall-shear-stress model that directly operates on the unstructured mesh of the simulation. This is achieved by using a GNN wherein the mesh of the simulation is represented as a graph. Neural networks operating on graphs were first used by Gori et al. (2005) and have demonstrated their effectiveness in various tasks, for instance, visual scene understanding, multi-agent systems, relational reasoning, or graph classification (Scarselli et al., 2008; Bruna et al., 2013; Sukhbaatar et al., 2016; Bronstein et al., 2017; Gilmer et al., 2017; Raposo et al., 2017; Santoro et al., 2017; Kipf et al., 2018). GNNs were also used by various authors to learn the dynamics of physical systems, including fluid flows (Battaglia et al., 2016; Chang et al., 2016; Watters et al., 2017; Sanchez-Gonzalez et al., 2018; Van Steenkiste et al., 2018; Pfaff et al., 2020; Chen et al., 2021; Gao et al., 2022). A review may be found in Battaglia et al. (2018) or Zhou et al. (2020). The relevance of GNNs to wall-shear-stress modeling can be hypothesized from its potential ability to detect geometric features, such as angles and curvature. In addition, the architecture of GNNs encodes some physical assumptions that are relevant for a wall-shear-stress model and therefore may ease the training process and increase the generalizability of the model, such as translational invariance and an inherent bias toward spatial locality.

The ability of a GNN to model the WSS is investigated in this paper using both a priori tests, based on filtered numerical data, and a posteriori tests, based on LESs implementing such models. A database of seven direct numerical simulations or high-fidelity wall-resolved LESs is used to train the model and validate the model: two channel-flow simulations, the simulation of the flow in a diffuser, two simulations of flows over backward-facing steps (BFSs) and two simulations of flows in a linear blade cascade. This provides a variety of flow phenomena, including fully developed boundary layers, adverse pressure gradients, separated boundary layers, and laminar–turbulent transition. Instantaneous snapshots from each simulation are filtered onto coarse meshes representative of a LES with a wall-shear-stress model, of which only near-wall nodes are used by the model for its prediction. The analysis is restricted to coarse tetrahedral meshes, for their ubiquity in simulations of complex industrial flows, although the method could in principle be applied to any type of mesh. An Encode-Process-Decode architecture (Battaglia et al., 2018) is used for the GNN, with no global features to ensure the spatial locality of the prediction. The inputs of the model are the displacement vector of the mesh edges and the three components of velocity at each node, scaled to ensure the Mach number invariance of the model at the incompressible limit, and augmented by arbitrary rotations and reflections of the coordinate system. The model is implemented a posteriori by coupling the compressible massively parallel flow solver AVBP (Schönfeld and Rudgyard, 1999 to a code dedicated to inference, using a message passing interface (MPI; Serhani et al., 2022). A channel flow, a BFS and a linear blade cascade configuration are simulated for these a posteriori tests. In each case, the simulations are compared to the corresponding simulation using an algebraic model based on the law of the wall to assess the relevance of the model.

The paper is organized as follows. The dataset and the preparation of the data for the machine-learning model is presented in Section 2. The architecture of the GNN is described in Section 3. The results are analyzed and discussed in Section 4.

2. Dataset and Data Preparation

2.1. Dataset

The machine-learning wall-shear-stress model developed in this paper relies on a database of direct numerical simulations or high-fidelity wall-resolved LESs which have been produced by various research groups (Figure 1). Namely, the database includes the simulations of: a fully developed channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=180 $ (CF1), performed at Imperial College London (Agostini and Vincent, 2020) using the high-order flux reconstruction method of Huynh (2007) in the CFD code PyFR (Witherden et al., 2014); a fully developed channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ (CF2), performed at the Polytechnic University of Madrid (Del Álamo and Jiménez, 2003; Lozano-Durán and Jiménez, 2014; 2015) using a hybrid Fourier–Chebyshev spectral method; a three-dimensional diffuser (3DD) corresponding to the geometry “Diffuser 1” of Cherry et al. (2008), performed at Barcelona Supercomputing Center (Ercoftac, 2022) using the low-dissipation finite element scheme of Lehmkuhl et al. (2019); a BFS, performed at CERFACS (Pouech et al., 2019; Pouech et al., 2021) using a cell-vertex finite-element method (Schönfeld and Rudgyard, 1999 with second-order accurate convection and diffusion schemes (Lax and Wendroff, 1960); a curved BFS based on the geometry of Disotell and Rumsey (Disotell and Rumsey, 2017; Alaya et al., 2020), hereafter referred to as adverse-pressure-gradient (APG) simulation and performed at the University of Bergamo (Ercoftac, 2022) using discontinuous Galerkin method and a fifth order linearly implicit Rosenbrock scheme (Bassi et al., 2015; 2016); and a NACA 65-009 blade cascade (N65) such as studied experimentally by Ma et al. (2011) and Zambonini et al. (2017), performed using a cell-vertex finite-element method (Schönfeld and Rudgyard, 1999 and a two-step Taylor–Galerkin scheme (Colin and Rudgyard, 2000). The N65 case is subdivided into two sub-configurations that differ only by the inlet boundary condition: a simulation with an incidence angle of 4° (N65a) and a simulation with an incidence angle of 7° (N65b). The two simulations are described in Appendix C. The database is heterogeneous in terms of flow solvers and numerical methods since the 3DD, BFS, N65a, and N65b simulations were performed using finite-element methods, the CF1 and APG simulations using high-order methods and the CF2 simulation using a spectral method. We believe that data-driven modeling approaches should demonstrate their ability to learn from such heterogeneous databases because heterogeneous databases can more easily scale and cover a large diversity of flow configurations. The numerical parameters of each simulation are summarized in Table 1.

Figure 1. Geometry of the six simulations in the database. The blue regions are separated in the mean.

Table 1. Numerical parameters of the numerical simulation included in the a priori database.

Note. The acronym FE denotes a finite element method, DG a discontinuous Galerkin method, and $ {\mathrm{\mathbb{P}}}_n $ indicates the use of $ n $ th order polynomials in the case of high-order methods. Note that due to the use of high-order polynomial, the effective number of solution points is ( $ \times $ ) 15 M in the CF1 simulation and ( $ + $ ) 960 M in the APG simulation.

A wide variety of flow behavior is represented by the dataset. The CF1 and CF2 simulations are focused on providing data on fully developed boundary layers. The 3DD and BFS simulations include fully separated flow regions. The APG simulation includes regions of incipient separation, that is with a backflow probability between 1 and 20% (Simpson, 1989). The N65a and N65b simulations include transitional boundary layers, through small-scale flow separation, on the blade surface, a light adverse pressure gradient as well as a large region of unsteady three-dimensional corner separation on the junction between the endwall and the suction side of the blade. For a given model to produce accurate predictions in all simulations, an ability to discriminate between these various flow regimes will thus be required.

2.2. Data preparation

The instantaneous fields of the a priori database need to be preprocessed for use by the machine-learning procedure. Namely, the fields are filtered onto coarse meshes representative of a LES with a wall-shear-stress model, as is customary in a priori tests for LES (Leonard, 1974; Sagaut, 2006). We focus on tetrahedral meshes, for their simplicity, their ability to represent arbitrarily complex domains and widespread use in simulations of complex industrial flows. The preparation of the data for the machine-learning algorithm is represented in Figure 1 and is as follows:

1. 1. For each simulation, a series of coarse tetrahedral meshes are constructed using the mesh adaptation library MMG3D (Dobrzynski and Frey, 2008; Dapogny et al., 2014). All coarse tetrahedral meshes have a uniform target refinement, although there are local variations of cell sizes resulting from the meshing algorithm, since regular tetrahedra do not fill space (as may be seen in Figure 2). To select the refinement level, a particular flow region is used to define the nominal WSS of each simulation. Those are: the channel walls in the simulations CF1 and CF2, the bottom wall upstream of the expansion in the 3DD simulation, the boundary layer upstream of the step in the APG and BFS simulations, and the boundary layer upstream of the blade in the N65a and N65b simulations. Coarse tetrahedral meshes are produced for 12 different refinement levels for each simulation, with a nominal edge length in wall units, in the regions listed above, in the range $ 25\hskip0.35em \le \hskip0.35em {e}^{+}\hskip0.35em \le \hskip0.70em 80 $ . This implies that the trained GNN can only operate on meshes that are within this range of mesh refinement in terms of scaled edge length.

2. 2. For each coarse tetrahedral mesh, the neighborhood of the target walls is selected, by preserving only nodes separated from the target walls $ \mathfrak{W} $ by $ {N}_H $ edges or less. The resulting tetrahedral meshes are conformal and referred to as coarse near-wall meshes hereafter. Unless otherwise specified, $ {N}_H $  = 3 is used.

3. 3. A fine tetrahedral mesh is produced, associated with a nominal edge length $ {e}^{+}=4 $ . The instantaneous fields of the a priori database are interpolated linearly onto this fine mesh. The resulting fields are used as source data for the filtering operation.

4. 4. The instantaneous data are filtered onto and at the resolution of the coarse near-wall meshes. This low-pass filtering operation attenuates large frequencies, associated with turbulent scales that cannot be resolved using the target computational grid. A surface filter is used to filter the WSS and a volume filter is used to filter the velocity components in the bulk of the domain. We denote ( $ \overline{\cdot} $ ) the volume filter associated with the local filter width $ {\tilde{\Delta}}_{p_0}={\tilde{V}}_{p_0}^{1/3} $ and ( $ {\overline{\cdot}}^S $ ) the surface filter associated with the local filter width $ {\tilde{\Delta}}_{p_0}^S={\tilde{S}}_{p_0}^{1/2} $ , with $ {\tilde{V}}_{p_0} $ and $ {\tilde{S}}_{p_0} $ , respectively, the nodal volume of node $ {p}_0 $ and the nodal area of $ {p}_0 $ on the surface $ \mathfrak{W} $ . The volume and surface filters are truncated Gaussian filters weighted by the local cell volume or face area,

(1) $$ \overline{\phi}\left({p}_0\right)=\sum \limits_{p\in {\mathfrak{w}}_{p_0}}{C}_{p_0}{V}_p{e}^{-\frac{{\mathbf{r}}_p-{\mathbf{r}}_{p_0}}{2{\left({\sigma}_{p_0}\right)}^2}}\phi (p), $$

(2) $$ {\overline{\phi}}^S\left({p}_0\right)=\sum \limits_{p\in {\mathfrak{w}}_{p_0}^S}{C}_{p_0}^S{S}_p{e}^{-\frac{{\mathbf{r}}_p-{\mathbf{r}}_{p_0}}{2{\left({\sigma}_{p_0}^S\right)}^2}}\phi (p), $$

with $ {C}_{p_0} $ and $ {C}_{p_0}^S $ normalization constants, selected such that $ \overline{a}={\overline{a}}^S=a $ for any constant $ a $ , $ {V}_p $ , and $ {S}_p $ , respectively, the nodal volume and face area associated with the source fine mesh and $ {\sigma}_{p_0}={\tilde{\Delta}}_{p_0}/\sqrt{12} $ and $ {\sigma}_{p_0}^S={\tilde{\Delta}}_{p_0}^S/\sqrt{12} $ , respectively, the standard deviations of the volume and surface Gaussian kernels, following Leonard (1974). The volume filter window $ {\mathfrak{w}}_{p_0} $ associated with a point $ {p}_0 $ of the target coarse near-wall mesh is restricted in the wall-normal direction to avoid biasing the mean velocity profile:

(3) $$ {\mathfrak{w}}_{p_0}=\left\{p\in {\mathfrak{P}}_{\mathfrak{f}}:\parallel {\boldsymbol{r}}_p-{\boldsymbol{r}}_{p_0}\parallel \hskip0.35em \le \hskip0.35em {\tilde{\Delta}}_{p_0}\wedge |{D}_{\mathfrak{W},p}-{D}_{\mathfrak{W},{p}_0}|\hskip0.35em \le \hskip0.35em \frac{1}{2}{\tilde{\Delta}}_{1,{p}_0}\right\}, $$

with $ {\mathfrak{P}}_{\mathfrak{f}} $ the set of points of the source mesh, $ {\boldsymbol{r}}_p $ the position vector associated with point $ p $ , and $ {D}_{\mathfrak{W},p} $ the shortest distance between $ p $ and the target walls $ \mathfrak{W} $ . The window $ {\mathfrak{w}}_{p_0}^S=\left\{p\in {\mathfrak{P}}_{\mathfrak{f}}:\parallel {\boldsymbol{r}}_p-{\boldsymbol{r}}_{p_0}\parallel \le {\tilde{\Delta}}_{p_0}^S\wedge {D}_{\mathfrak{W},p}=0\right\} $ of the surface filter is isotropic.

1. 5. The Gaussian filter defined above leads to a small bias on the velocity profile. In particular, the mean Gaussian-filtered velocity is slightly lower than mean unfiltered velocity in the case of a fully developed boundary layer. This bias is compensated for by correcting the filtered velocity as follows:

(4) $$ {\overline{\boldsymbol{u}}}^{\ast }=\overline{\boldsymbol{u}}+\left\langle \boldsymbol{u}\right\rangle -\left\langle \overline{\boldsymbol{u}}\right\rangle, $$

with $ {\overline{\boldsymbol{u}}}^{\ast } $ the corrected filtered velocity and where $ \left\langle \boldsymbol{u}\right\rangle $ and $ \left\langle \overline{\boldsymbol{u}}\right\rangle $ are assessed using the time-average fields from each simulation.

1. 6. The resulting data are partitioned into contiguous chunks of similar size using the graph partitioning library METIS (Karypis and Kumar, 1998). This helps producing mini-batches of consistent size to fit in the dynamic memory of the graphics processing unit (GPU) during training. In addition, subdividing the data into smaller chunks allows for more diversity in each mini-batch, by aggregating small regions of different simulations, which reduces noise in the optimization process.

Figure 2. Preprocessing steps required to prepare the data for the training procedure.

3. Graph Neural Network

The wall-shear-stress model developed in this paper is a GNN, and more specifically uses an Encode-Process-Decode architecture (Battaglia et al., 2018). The GNNs used in the present study are feedforward deep learning models that aim to exploit the structure of the data in high dimension to predict the shear stress at the wall. Indeed, GNNs can involve a relatively large number of points in the prediction, which immerses the learning problem in a high-dimensional space. The high-dimensional nature of the learning problem may reduce the inherent ambiguities of the WSS prediction, by allowing the model to discriminate between different types of local flow configurations. For instance, a minimal capacity would be needed to allow the model to identify that a boundary layer is separated (Dupuy et al., 2022). In addition, the use of information in a region for the prediction might allow the model to take into account the space displacement associated with regions of high correlations with the WSS (Boxho et al., 2022). More specifically, GNNs are used for their ability to work directly with unstructured mesh data, without requiring the prior interpolation of the data. This ability is important for LESs in complex geometries, with angles, corners, or curvature. Finally, the architecture of the GNNs is physically well suited to the modeling of the WSS. First, the model is biased toward spatial locality, such that nearby grid points will more easily influence the prediction than grid points that are far away. Second, the model is invariant to translation as the same learned function is applied everywhere in the domain. In addition to the physical assumptions related to the model architecture, the present learning process ensures the equivariance of the model to Galilean transformations, orthogonal transformations such as rotations or reflections, and “changes of scales,” hereafter referred to as Mach number equivariance. We furthermore assume that the flow is isothermal incompressible and obeys the incompressible Navier–Stokes equations without external body force. This is consistent with the selection of simulations in the database of Section 2.1, which involves neither significant temperature variations nor compressibility effects.

3.1. Model architecture

Each sample mesh produced in Section 2.2 is represented as a symmetric directed graph $ \mathfrak{G}=\left(\mathfrak{V},\mathfrak{E}\right) $ , where $ \mathfrak{V} $ is the set of nodes and $ \mathfrak{E} $ the set of directed edges. The number of nodes is denoted $ {n}_v $ and the number of edges $ {n}_e $ . Let $ {v}_i\hskip0.35em \in \hskip0.35em \mathfrak{V} $ be the $ i $ th node of the graph and $ {e}_k=\left({v}_{S(k)},{v}_{R(k)}\right)\in \mathfrak{E} $ an edge pointing from $ {v}_{S(k)} $ to $ {v}_{R(k)} $ . The input of the model is the ordered pair $ \boldsymbol{X}=\left(\boldsymbol{V},\boldsymbol{E}\right) $ , where $ \boldsymbol{V}\in {\mathrm{\mathbb{R}}}^{n_v\times {d}_v} $ is the node input features, whose $ i $ th element $ {V}_i\hskip0.35em \in \hskip0.35em {\mathrm{\mathbb{R}}}^{d_v} $ represents the input features at node $ {v}_i $ , and $ \boldsymbol{E}\in {\mathrm{\mathbb{R}}}^{n_e\times {d}_e} $ the edge input features, whose $ k $ th element $ {E}_k\hskip0.35em \in \hskip0.35em {\mathrm{\mathbb{R}}}^{d_e} $ represents the input features at edge $ {e}_k $ . The output of the model is the node vector $ \boldsymbol{Y}\in {\mathrm{\mathbb{R}}}^{n_v\times 1} $ . With an Encode-Process-Decode architecture (Battaglia et al., 2018), the output of the GNN may be expressed as

(5) $$ \boldsymbol{Y}={f}_{\gamma}\left(\boldsymbol{X}\right)={f}_{\delta}\hskip0.35em \circ \hskip0.35em {f}_{\pi}^N\circ \hskip0.35em \dots \hskip0.35em \circ \hskip0.35em {f}_{\pi}^2\hskip0.35em \circ \hskip0.35em {f}_{\pi}^1\hskip0.35em \circ \hskip0.35em {f}_{\varepsilon}\left(\boldsymbol{X}\right), $$

where $ {f}_{\varepsilon },{f}_{\pi}^1,{f}_{\pi}^2,\dots, {f}_{\pi}^N $ and $ {f}_{\delta } $ are the encoder, processor, and decoder, respectively (Figure 3). The architecture of the model is described below:

Figure 3. Graphical representation of the Encode-Process-Decode architecture.

3.1.1. Encoder

The input features are encoded into a latent space of size $ {d}_L=128 $ for both nodes and edges, $ {\boldsymbol{L}}^{\mathbf{0}}=\left({\hat{\boldsymbol{V}}}^{\mathbf{0}},{\hat{\boldsymbol{E}}}^{\mathbf{0}}\right)={f}_{\varepsilon}\left(\boldsymbol{X}\right) $ , with $ {\hat{\boldsymbol{V}}}^{\mathbf{0}}\in {\mathrm{\mathbb{R}}}^{n_v\times {d}_L} $ and $ {\hat{\boldsymbol{E}}}^{\mathbf{0}}\in {\mathrm{\mathbb{R}}}^{n_e\times {d}_L} $ . The encoding is performed independently at each node or edge by $ {\hat{V}}_i^0={f}_{\varepsilon}^V\left({V}_i\right) $ and $ {\hat{E}}_k^0={f}_{\varepsilon}^E\left({E}_k\right) $ , where $ {f}_{\varepsilon}^V $ and $ {f}_{\varepsilon}^E $ are two multilayer perceptrons.

3.1.2. Processor

The latent space is transformed by the successive application of $ N=4 $ message-passing steps, with residual connections. The latent space after $ n $ steps, $ {\boldsymbol{L}}^{\boldsymbol{n}}=\left({\hat{\boldsymbol{V}}}^{\boldsymbol{n}},{\hat{\boldsymbol{E}}}^{\boldsymbol{n}}\right)={\boldsymbol{L}}^{\boldsymbol{n}-\mathbf{1}}+{f}_{\pi}\left({\boldsymbol{L}}^{\boldsymbol{n}-\mathbf{1}}\right) $ , is given by $ {\hat{V}}_i^n={\hat{V}}_i^{n-1}+{f}_{\pi}^V\left({\hat{V}}_i^{n-1},{\sum}_{k\mid R(k)=i}{\hat{E}}_k^n\right) $ and $ {\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 $ {f}_{\pi}^V $ and $ {f}_{\pi}^E $ are two multilayer perceptrons.

3.1.3. Decoder

The latent space after the last message-passing step is decoded back to physical space at nodes only, $ {Y}_i={f}_{\delta}^V\left({\hat{V}}_i^N\right) $ , with $ {f}_{\delta}^V $ a multilayer perceptron.

The same encoding, processing, and decoding functions are learned and applied everywhere in the graph. Hence, the message-passing steps in the processor will propagate the information through the edge and nodes to predict the WSS on every wall nodes.

The multilayer perceptrons $ {f}_{\varepsilon}^V $ , $ {f}_{\varepsilon}^E $ , $ {f}_{\pi}^V $ _, and $ {f}_{\pi}^E $ each have $ {n}_{\mathrm{\ell}}=2 $ layers composed of $ {d}_L $ neural units and followed by layer normalization (Ba et al., 2016). The output of the $ i $ th neuron of hidden layer $ \mathrm{\ell} $ is given by

(6) $$ {\overline{z}}_{\mathrm{\ell}}^{(i)}={\overline{w}}_{\mathrm{\ell}}^{(i)}\frac{z_{\mathrm{\ell}}^{(i)}-{\left\langle {z}_{\mathrm{\ell}}^{(i)}\right\rangle}_{\mathfrak{F}}}{\sqrt{{\left\langle {\left({z}_{\mathrm{\ell}}^{(i)}-{\left\langle {z}_{\mathrm{\ell}}^{(i)}\right\rangle}_{\mathfrak{F}}\right)}^2\right\rangle}_{\mathfrak{F}}}}+{\overline{b}}_{\mathrm{\ell}}^{(i)}, $$

where $ {\left\langle \cdot \right\rangle}_{\mathfrak{F}} $ denotes an average across features and with

(7) $$ {z}_{\mathrm{\ell}}^{(i)}=\sigma \left(\sum \limits_j{w}_{\mathrm{\ell}}^{\left(i,j\right)}{\overline{z}}_{\mathrm{\ell}-1}^{(j)}+{b}_{\mathrm{\ell}}^{(i)}\right), $$

where the weights and biases $ {\overline{w}}_{\mathrm{\ell}}^{(i)} $ , $ {w}_{\mathrm{\ell}}^{\left(i,j\right)} $ , $ {\overline{b}}_{\mathrm{\ell}}^{(i)} $ , and $ {b}_{\mathrm{\ell}}^{(i)} $ are learnable parameters and $ \sigma $ the activation function, namely a rectified linear unit $ \sigma (x)=\max \left(0,x\right) $ in the present work. The multilayer perceptron $ {f}_{\delta}^V $ has $ {n}_{\mathrm{\ell}}+1 $ hidden layers, the first $ {n}_{\mathrm{\ell}} $ layers being as described above and the last layer a linear layer with one neural unit,

(8) $$ {z}_{n_{\mathrm{\ell}}}=\sum \limits_j{w}_{n_{\mathrm{\ell}}}^{(j)}{z}_{\left({n}_{\mathrm{\ell}}-1\right)}^{(j)}+{b}_{n_{\mathrm{\ell}}}. $$

The weights of the network are initialized using a truncated normal distribution, discarding and redrawing values beyond two standard deviations, and trained using the Adam optimizer (Kingma and Ba, 2015) with a base learning rate of 0.001. The loss function is the mean squared error (MSE) between the predicted WSS and the reference WSS,

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

where $ {\left\langle \cdot \right\rangle}_{\mathfrak{A}} $ denotes an average across graph nodes and samples and $ \chi $ a target boundary mask, with $ \chi =1 $ on the target wall nodes $ \mathfrak{W} $ and $ \chi =0 $ elsewhere. The training time is 20 hr using a single NVIDIA A30 GPU unit.

3.2. Feature selection

The input of the GNN WSS model is given by $ {d}_v=4 $ node input features and $ {d}_e=4 $ edge input features. The node input features are the target boundary mask $ \chi $ and the three components of the scaled velocity vector $ {\overset{\check{}}{\overline{\boldsymbol{u}}}}^{\ast } $ , defined as

(10) $$ {\overset{\check{}}{\overline{\boldsymbol{u}}}}^{\ast }={\overline{\boldsymbol{u}}}^{\ast }{\left\langle \parallel \mathbf{e}\parallel \right\rangle}_{\mathfrak{G}}/{\left\langle \nu \right\rangle}_{\mathfrak{G}}, $$

where $ {\left\langle \parallel \boldsymbol{e}\parallel \right\rangle}_{\mathfrak{G}} $ is the average edge length of a given graph $ \mathfrak{G} $ and $ {\left\langle \nu \right\rangle}_{\mathfrak{G}} $ the average kinematic viscosity. Note that the average viscosity is included in the scaling since, although all the flows considered are incompressible, some of the numerical methods used are based on the compressible Navier–Stokes equations, which results in negligible but nonzero variations of viscosity. Following Pfaff et al. (2020), the $ {d}_e=4 $ edge input features are the three components of the scaled displacement vector of the edge $ \overset{}{\check{\boldsymbol{e}}} $ , with $ \overset{}{\check{\boldsymbol{e}}}=\boldsymbol{e}/{\left\langle \parallel \boldsymbol{e}\parallel \right\rangle}_{\mathfrak{G}} $ , and its magnitude $ \parallel \overset{}{\check{\boldsymbol{e}}}\parallel $ . The edge and node input features are taken on the entire graph of the neural network. However, only a limited region will influence the prediction of the WSS at a given node. This region is the receptive field of the model, governed by the number of message-passing steps $ N $ as represented schematically in Figure 4. In accordance with the scaling of the input velocity (10), the target WSS is scaled as

(11) $$ {\overline{\tau}}_{\mathrm{ref}}^S={\overline{\tau}}_{\mathrm{ref}}^S{\left({\left\langle \parallel \boldsymbol{e}\parallel \right\rangle}_{\mathfrak{G}}/{\left\langle \nu \right\rangle}_{\mathfrak{G}}\right)}^2/{\left\langle \rho \right\rangle}_{\mathfrak{G}}, $$

where $ {\left\langle \rho \right\rangle}_{\mathfrak{G}} $ is the average density and with $ \tau $ the norm of the shear stress vector $ \boldsymbol{\tau} =\Sigma \cdot {\boldsymbol{e}}_{\boldsymbol{n}}-({\boldsymbol{e}}_{\boldsymbol{n}}\cdot \Sigma \cdot {\boldsymbol{e}}_{\boldsymbol{n}}){\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.

Figure 4. Effect of the number of message-passing steps $ N $ on the receptive field of the graph neural network (red shaded area) for the prediction of the wall shear stress on a given target node (yellow node). (a) $ N=1 $ . (b) $ N=2 $ . (c) $ N=3 $ . (d) $ N=4 $ .

The equivariance of the model under several transformations is ensured by the selection of input features, as reported in Table 2. The orthogonal equivariance of the model, that is its equivariance under rotations and reflections, is ensured by data augmentation: for each sample graph, a random three-dimensional rotation and reflection is applied onto the input edge displacement vectors and velocity vectors. To ensure that the predictions of the model do not depend on the coordinate system, the components of the scaled edge displacement vector $ \overset{}{\check{\boldsymbol{e}}} $ and scaled velocity vector $ {\overset{\check{}}{\overline{\boldsymbol{u}}}}^{\ast } $ are expressed in a basis arbitrarily rotated and reflected for each graph. The Galilean invariance of the model is ensured by using the wall as an intrinsic velocity reference, that is, by using as input feature, the relative velocity of the flow with respect to the wall velocity. The Mach number equivariance of the model, or namely the equivariance of the model to a change of scale of velocity, density, viscosity, or length, is ensured by the scaling ( $ \overset{\smile }{\cdot } $ ) of the input and target features. Indeed, the scaled variables tend to a constant as the Mach number of the flow tends to zero, at constant Reynolds number and mesh resolution (Paolucci, 1982). In other words, for an isothermal incompressible flow without external body force, the scaled variables depend on the Reynolds number of the flow, which governs the flow physics, but are unaffected by a change of Mach number, which is not a relevant physical parameter. This scaling is a natural extension of the scaling used by Dupuy et al. (2022) for a Cartesian grid. It is purely local and only involves quantities that can be computed in a LES. Enforcing the equivariance of the model under these transformations increases the sample efficiency of each graph in the numerical dataset and improves the generalization of the model to different configurations.

Table 2. Strategies used to enforce the equivariance of the machine-learning model to various transformations.

4. Results

4.1. A priori tests

This section presents a priori tests of the machine-learning wall-shear-stress model, based solely on the numerical database described in Section 2.1. The model is trained using the simulations CF1, CF2, 3DD, and N65a. In each training simulation, a subset of the data is isolated for validation purpose. The split between training and validation datasets is performed independently for each simulation. In the CF1, CF2, and 3DD simulations, the snapshots associated with the later time steps are reserved for validation purpose. The BFS and APG simulations are only used as test simulations, to test the performance of the model on configurations not seen during training. For the N65 configuration, training is exclusively performed on the sub-configuration N65a, that is with an incidence angle $ I=4{}^{\circ} $ , while testing is exclusively performed on the sub-configuration N65b, that is with an incidence angle $ I=7{}^{\circ} $ . This is summarized in Table 3.

Table 3. Simulations used for the training and validation of the graph neural network WSS model, for a priori tests and for a posteriori tests.

In order to assess a priori the relevance of the GNN WSS model, the behavior of the model in each simulation is first analyzed using $ {u}^{+}\left({y}^{+}\right) $ graph wherein the predicted local WSS is used to compute the wall-unit scaling ( $ {}^{+} $ ). The results are reported in Figure 5 along with corresponding data based on the reference WSS. Agreement or departure of both reference data and model prediction from fully developed equilibrium wall turbulence is emphasized by comparing the instantaneous data points to the analytical profile of Reichardt (1951),

(12) $$ {u}^{+}\left({y}^{+}\right)=\frac{1}{\kappa}\log \left(1+\kappa {y}^{+}\right)+7.8\left[1-\exp \left(\frac{-{y}^{+}}{11}\right)-\frac{y^{+}}{11}\exp \left(\frac{-{y}^{+}}{3}\right)\right], $$

with $ \kappa =0.41 $ the von Kármán constant. It should be emphasized that, except in the case of the channel-flow simulations, Reichardt’s law is not a reference but rather as a baseline, indicative of what a model with insufficient capacity would predict. In the two channel-flow simulations (CF1 and CF2), the reference data points are scattered around the law of the wall, with a large variance associated with velocity and WSS fluctuations. This behavior is reproduced by the model, but the variance of the data points around Reichardt’s law is clearly underestimated. In the 3DD simulation, most data points are scattered below Reichardt’s law in the $ {u}^{+}\left({y}^{+}\right) $ graph, which is consistent with the presence of separated regions. The GNN WSS model reproduces this behavior, although the predictions are more biased toward the law of the wall and extreme points well above Reichardt’s law are missing at low $ {y}^{+} $ . The presence of data points that deviates from the law of the wall is also correctly predicted by the model in the test BFS and APG simulations, suggesting reliable generalization to other flow configurations with separation. Finally, the NACA 65-009 blade simulation (N65b) is most saliently associated with a transitional boundary layer on the blade surface, and most data points are found above Reichardt’s law on the $ {u}^{+}\left({y}^{+}\right) $ graph. This is qualitatively reproduced by the model but with a clearly lower variance. Overall, the model thus leads to markedly different distributions of $ {u}^{+}\left({y}^{+}\right) $ in each simulation, indicating an ability to locally discriminate, to some extent, between the various type of boundary layers, fully developed, separated, or transitional, present in the numerical database.

Figure 5. A priori validation: Norm of the scaled tangential velocity $ {u}^{+} $ as a function of the scaled distance to the wall $ {y}^{+} $ in the CF1, CF2, 3DD, BFS, APG, and N65b datasets using the local reference wall shear stress (left) or the prediction of the graph neural network WSS model (right) to compute the wall unit scaling ( $ {}^{+} $ ). The red line is Reichardt’s law, given by equation (12).

Several quantitative measures of the agreement between the predicted WSS and corresponding local instantaneous reference WSS are given in Table 4. First, we use the concordance coefficient of Lin (1989) between the predicted WSS $ Y $ and the reference WSS $ {\tau}_{\mathrm{ref}} $ , defined as

(13) $$ {C}_{\tau }=\frac{{\left\langle {\tau}_{\mathrm{ref}}Y\right\rangle}_{\mathfrak{A}}-{\left\langle {\tau}_{\mathrm{ref}}\right\rangle}_{\mathfrak{A}}{\left\langle Y\right\rangle}_{\mathfrak{A}}}{{\left\langle {\tau}_{\mathrm{ref}}^2\right\rangle}_{\mathfrak{A}}-{\left\langle {\tau}_{\mathrm{ref}}\right\rangle}_{\mathfrak{A}}^2+{\left\langle {Y}^2\right\rangle}_{\mathfrak{A}}-{\left\langle Y\right\rangle}_{\mathfrak{A}}^2+{\left({\left\langle {\tau}_{\mathrm{ref}}\right\rangle}_{\mathfrak{A}}-{\left\langle Y\right\rangle}_{\mathfrak{A}}\right)}^2}. $$

Table 4. A priori validation: Integral measures of the disagreement between the reference wall shear stress and the prediction of the graph neural network WSS model in the CF1, CF2, 3DD, BFS, APG, and N65b datasets.

In addition, the explained variance $ {R}_2 $ is used to compare the variance of the residual of the model to the underlying WSS variance:

(14) $$ {R}_2=1-\frac{{\left\langle {\left({\tau}_{\mathrm{ref}}-Y\right)}^2\right\rangle}_{\mathfrak{A}}}{{\left\langle {\left({\tau}_{\mathrm{ref}}-\left\langle {\tau}_{\mathrm{ref}}\right\rangle \right)}^2\right\rangle}_{\mathfrak{A}}}. $$

The Reichardt-residual explained variance similarly compares the model residual to the residual that would be obtained using Reichardt’s law to compute the WSS,

(15) $$ {R}_{2,\mathrm{Reichardt}}=1-\frac{{\left\langle {\left({\tau}_{\mathrm{ref}}-Y\right)}^2\right\rangle}_{\mathfrak{A}}}{{\left\langle {\left({\tau}_{\mathrm{ref}}-{\tau}_{\mathrm{Reichardt}}\right)}^2\right\rangle}_{\mathfrak{A}}}. $$

The concordance coefficient $ {C}_{\tau } $ is equal to zero if there is no linear correlation between the model and the reference and equal to $ 1 $ if and only if the predicted WSS is equal to the reference WSS for each time step. The explained variance $ {R}_2 $ is equal to zero if the variance of the residual is equal to the reference shear stress variance and positive if it is lower. In other words, the coefficient $ {R}_2 $ compares the model to a constant model that always returns the mean value of the data. The Reichardt-residual explained variance $ {R}_{2,\mathrm{Reichardt}} $ instead compares the model to Reichardt’s law (12). Namely, $ {R}_{2,\mathrm{Reichardt}} $ is equal to zero is the variance of the residual is equal to the variance of the residual of a model based on Reichardt’s law, positive if it is lower, and negative if Reichardt’s law is more accurate than the GNN WSS model. This is particularly relevant in flows where the analytical profile of Reichardt (1951) is suitable, such as fully developed turbulent channel flows. As for $ {C}_{\tau } $ , both $ {R}_2 $ and $ {R}_{2,\mathrm{Reichardt}} $ are equal to $ 1 $ for an ideal model. Note that the values obtained for different datasets can only be compared with great care, since these three indicators depend not only on the intrinsic performance of the model but also on the underlying distribution of the data. For instance, the value of $ {C}_{\tau } $ associated with Reichardt’s law is larger in the 3DD and BFS simulations than in the TCF1 and TCF2 simulations, although the model performs worse in these two simulations. In all simulations, the GNN WSS model agrees well with the reference WSS and improves by all measures the results compared to Reichardt’s law baseline. Even in the case of a channel flow (case CF1 or CF2), for which the law of the wall is most well suited, the GNN WSS model reduces the residual variance of a model based on Reichardt’s law by around a third. Figure 6 shows that this variance reduction is consistently obtained for a wide range of distances to the wall. In the 3DD, BFS, and N65b simulations on the other hand, variations of model performance with $ {y}^{+} $ may be observed. In particular, very low value of $ {y}^{+} $ are associated with a very large value of $ {R}_{2,\mathrm{Reichardt}} $ because those corresponds to separated flow regions that are not accurately described by Reichardt’s law.

Figure 6. A priori validation: Reichardt-residual explained variance ( $ {R}_{2,\mathrm{Reichardt}} $ ) in the CF1, CF2, 3DD, BFS, APG, and N65b datasets conditioned on the scaled distance to the wall $ {y}^{+} $ .

Finally, the instantaneous agreement between the predicted and reference WSS is reflected in a improvement of the time- and space-average predicted WSS, as may be seen Figure 7. Note that the profiles of WSS in Figure 7 have not been scaled in order to explicitly show that the model can operate on flows that have widely different scales, thanks to the scaling used for the model inputs. In accordance with the instantaneous results, a law-of-the-wall model based on Reichardt’s law greatly underestimates the mean WSS in the BFS simulation after $ x/{h}_s=2 $ and overestimates the mean WSS in the N65b simulation upstream of the blade and near the leading edge. In the APG simulation, the mean WSS is overestimated in the region of positive curvature ( $ x/{R}_{\mathrm{max}}<1 $ ) and slightly underestimated in the region of negative curvature ( $ x/{R}_{\mathrm{max}}>1 $ ). In each case, the GNN WSS model improves the average WSS prediction and leads to predictions that are close to the reference value. This demonstrates that the model is able to leverage the instantaneous spatial information to improve the mean WSS prediction, even when generalizing to cases that are not seen in the training dataset. A parametric study of the influence of the number of message-passing steps $ N $ on this ability has been performed. This parameter is related to the capacity of the GNN since it controls the size of the receptive field of the GNN WSS model (4). The results, provided in Appendix A, demonstrate that a number of message-passing steps $ N\hskip0.35em \ge \hskip0.35em 3 $ is required to properly predict the WSS in the various flows included in our dataset. This is intuitively understood by the fact that the GNN WSS model needs to use information in a relatively large region to immerse the learning problem in the high dimensional space required to reduce the ambiguities associated with the WSS prediction and discriminate between various flow configurations. The value $ N=4 $ used in the present study is thus suitable.

Figure 7. A priori validation: Average prediction of a model based on Reichardt’s law and of the graph neural network WSS model. The average is performed in both time and spanwise direction in the BFS simulation (a), the APG simulation (b), and the N65b simulation (c) on the blade surface. ( $ \dagger $ ) On the blade surface of the N65b simulation, the mean wall shear stress is for clarity reported with positive values for the suction side of the blade and negative values for the pressure side of the blade.

4.2. A posteriori tests

This section investigates the machine-learning wall-shear-stress model a posteriori, that is, by performing simulations implementing the GNN. All simulations presented in this section are performed using a compressible, unstructured and massively parallel flow solver (AVBP; Schönfeld and Rudgyard, 1999), coupled to a code dedicated to the inference of the GNN. The coupling approach is described in Serhani et al. (2022), but the computational cost of the model is also briefly discussed in Appendix B. At each time step, data are extracted from the fields of the simulation in the vicinity of walls to create near-wall graphs, sent to the GNN using the MPI. The inference is distributed to a number of GPU according to a specific data partitioning, independent of the CPU partitioning. The GPU partitions have an overlap corresponding to the number of message-passing steps $ N $ of the model, to ensure that each partition contains the receptive field required for the WSS computation (Figure 7), and accordingly that the model prediction is independent of the number of partitions. The resulting predicted WSS is then sent back to main flow solver to impose the boundary condition. The numerical method used in all simulations is similar and relies on the resolution of compressible LES equations with the ideal gas equation of state and an explicit time stepping scheme.

4.2.1. Governing equations

The governing equations of the LESs are:

(16) $$ {\partial}_t\rho +{\partial}_j\left(\rho {U}_j\right)=0, $$

(17) $$ {\partial}_t\left(\rho {U}_i\right)+{\partial}_j\left(\rho {U}_j{U}_i\right)=-{\partial}_iP+{\partial}_j{\Sigma}_{ij}, $$

(18) $$ {\partial}_t\left(\rho E\right)+{\partial}_j\left(\rho {U}_jH\right)=-{\partial}_j{Q}_j-{\partial}_j\left({U}_jP\right)+{\partial}_j\left({U}_j{\Sigma}_{ij}\right), $$

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

with $ \rho $ the filtered density, $ t $ the time, $ \boldsymbol{U} $ the Favre-filtered velocity, $ \boldsymbol{x} $ the Cartesian coordinate, $ P $ the filtered pressure, $ \boldsymbol{Q} $ the conductive heat flux, $ E $ the total energy per unit mass, $ H=E+P/\rho $ the total enthalpy per unit mass, and $ r $ the specific ideal gas constant. Density and temperature, which do not experience significant variations in our test cases, are hence resolved nonetheless by the numerical method. The fluid is air. No external body force or heat source are taken into account. Temperature is computed using tabulated data from Stull and Prophet (1971) based on internal energy $ e=E-\frac{1}{2}{U}_i{U}_i $ . 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, thus, given by

(20) $$ {\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 $ {\delta}_{ij} $ denotes the Kronecker delta.

4.2.2. Channel flow

Several test cases have been selected to verify a posteriori the relevance of the GNN WSS model. First, a simulation of a fully developed turbulent channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ is performed. The goal is to ensure that the model is able to match in a channel flow the performance of classical algebraic wall stress models. The domain is a rectangular channel of size $ 12.6{h}_c\times 2{h}_c\times 6.3{h}_c $ , where $ {h}_c $ is the half-height of the channel. The mesh is fully tetrahedral and contains 2 million grid points. It is referred to as the original mesh (O) hereafter. The mean edge length is $ {\Delta}^{+}=50 $ in wall units, and the mean first off-the-wall grid point height is $ {y}^{+}=37 $ , which is within the logarithmic layer. A cell-vertex finite-element method is used to discretize equations (16)–(19). The convective scheme is a two-step Taylor–Galerkin scheme with third-order spatial and temporal accuracy (Colin and Rudgyard, 2000), while the diffusive scheme is a centered second-order scheme. The subgrid-scale viscosity is computed using the Smagorinsky’s model (Smagorinsky, 1963). The profiles of mean streamwise velocity and standard deviation of streamwise velocity, given in Figure 8, clearly show that a simulation with a GNN WSS model leads to similar results to an algebraic model based on the law of the wall. For both models, there is, however, a non-negligible discrepancy between the LESs and the reference profiles of Hoyas and Jiménez (2008). In particular, there is a mismatch between the logarithmic layers in the velocity profile. This mismatch is well documented in the literature for shear stress models and is due to numerical and subgrid-modeling errors in the first few grid points (Larsson et al., 2016; Bose and Park, 2018). In the context of algebraic wall modeling, Kawai and Larsson (2012) proposed a method to address the issue of logarithmic-layer mismatch by using as input the velocity farther from the wall, at least three cells within the large-eddy-simulation grid, rather than the velocity at the first off-the-wall point.

Figure 8. A posteriori validation: (a) Mean streamwise velocity and (b) standard deviation of streamwise velocity in large-eddy simulations of a channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ with an algebraic wall stress model and a graph neural network WSS model. The unfiltered direct numerical simulation profile of Hoyas and Jiménez (2008) is given for comparison.

The approach of Kawai and Larsson (2012) can be adapted to the case of a GNN WSS model, by arguing that since the velocity close to the wall is more affected by numerical and subgrid-scale modeling errors, only the velocity far from the wall should be considered for the prediction of the GNN. This hypothesis is tested by training a dedicated GNN WSS model that only uses the grid points three neighbors away from the wall for its prediction, that is the grid points that are separated from the target walls $ \mathfrak{W} $ by exactly three edges. The a posteriori results obtained with this model in the case of a channel flow are given in Figure 9. The logarithmic-layer mismatch is greatly reduced using the GNN WSS model that ignores the grid points too close to the wall. This explicitly demonstrates that the approach of Kawai and Larsson (2012), well acknowledged for algebraic wall models, can be readily adapted to the case of a machine-learning wall-shear-stress model. However, this approach has the significant drawback of withholding valuable near-wall data from the model, that although not necessary in a channel flow may be more critical in more complex flows and geometries.

Figure 9. A posteriori validation: (a) Mean streamwise velocity and (b) standard deviation of streamwise velocity in large-eddy simulations of a channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ with a graph neural network WSS model without restrictions (GNN model) and a graph neural network WSS model that only uses the grid points at least three neighbors away from the wall for its prediction (GNN model, far inputs). The unfiltered direct numerical simulation profile of Hoyas and Jiménez (2008) is given for comparison.

The grid sensitivity of the results is assessed by performing a posteriori LESs of the channel flow for two additional levels of refinement, namely a finer mesh (F) associated with a mean first off-the-wall grid point height $ {y}^{+}=27 $ and a coarser mesh (C) associated with a mean first off-the-wall grid point height $ {y}^{+}=54 $ . With the F mesh (finer mesh), the GNN WSS model lead to a more accurate velocity profile than the law-of-the-wall model (Figure 10a). This is likely due to the fact that the majority of first off-the-wall grid points are located within the buffer layer, thus outside the domain of validity of the logarithmic law of the wall. With the C mesh (coarser mesh), the GNN WSS model and the law-of-the-wall model lead to almost identical results (Figure 10b). The logarithmic-layer mismatch observed in the velocity profile can, as with the original mesh O, be removed by applying the approach proposed by Kawai and Larsson (2012) to the GNN WSS model with the mesh C. However, this solution does not lead to a significant improvement with the mesh F, which exhibits a less pronounced logarithmic-layer mismatch.

Figure 10. A posteriori validation: Mean streamwise velocity and standard deviation of streamwise velocity in large-eddy simulations of a channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ with the (a) meshes F (finer mesh) and (b) C (coarser mesh). The unfiltered direct numerical simulation profile of Hoyas and Jiménez (2008) is given for comparison.

Finally, the simulation of a fully developed turbulent channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=2000 $ is performed to verify that the GNN WSS model can operate for a different Reynolds number, outside of the range of friction Reynolds numbers seen during training ( $ {\mathit{\operatorname{Re}}}_{\tau }=180 $ and 950). The geometry and the numerical method are the same as for the LES at $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ . The mesh refinement is similar to the mesh refinement of the original mesh (O) used at $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ . Namely, the mean first off-the-wall grid point height is $ {y}^{+}=37 $ . The profiles of mean streamwise velocity and standard deviation of streamwise velocity are given in Figure 11. The results demonstrate that the GNN WSS model performs as well as an algebraic model based on the equilibrium law of the wall at $ {\mathit{\operatorname{Re}}}_{\tau }=2000 $ . As in the simulation at $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ , the velocity profile features a logarithmic-layer mismatch in both simulations. The logarithmic-layer mismatch can be removed by applying the approach proposed by Kawai and Larsson (2012) to the GNN WSS model (Figure 11), as described above for the case $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ . For clarity, we stress that although we did not perform simulations to demonstrate it, the approach of Kawai and Larsson (2012) could be applied to remove the logarithmic-layer mismatch with an algebraic law-of-the-wall model, as was proposed and performed in the literature Kawai and Larsson (2012). Overall, the performance of the present GNN WSS model is thus satisfactory for the a posteriori simulation of a channel flow at both $ {\mathit{\operatorname{Re}}}_{\tau }=950 $ and $ {\mathit{\operatorname{Re}}}_{\tau }=2000 $ .

Figure 11. A posteriori validation: (a) Mean streamwise velocity and (b) standard deviation of streamwise velocity in large-eddy simulations of a channel flow at friction Reynolds number $ {\mathit{\operatorname{Re}}}_{\tau }=2000 $ with an algebraic wall stress model, a graph neural network WSS model without restrictions (GNN model) and a graph neural network WSS model that only uses the grid points at least three neighbors away from the wall for its prediction (GNN model, far inputs). The unfiltered direct numerical simulation profile of Hoyas and Jiménez (2008) is given for comparison.

4.2.3. Backward-facing step

The LES of the flow over a BFS using a GNN WSS model is addressed in this section. The configuration is similar to the simulation BFS described in Section 2.1, which was not included in the training dataset. The expansion ratio $ ER=\left({L}_y+{h}_s\right)/{L}_y $ of the step, where $ {L}_y $ is the inlet channel height and $ {h}_s $ the step height, is equal to 1.2. The Reynolds number based on $ {h}_s $ and the inflow bulk velocity is $ {\mathit{\operatorname{Re}}}_h=5100 $ . The inlet is located 10 $ {h}_s $ before the step. At the inlet, the mean streamwise velocity profile of the wall-resolved simulation of Dupuy et al. (2022) is imposed with the characteristic nonreflecting boundary condition of Poinsot and Lele (1992). In addition, isotropic fluctuations are imposed following the turbulence kinetic energy profile of the wall-resolved simulation using the approach of Kraichnan (1970). At the outlet, located 20 $ {h}_s $ after the step, a constant pressure of 1 bar is imposed. The spanwise length of the domain is 12 $ {h}_s $ . Periodic boundary conditions are used on the side walls. A symmetric boundary condition, that is a zero shear stress and a zero wall-normal velocity, is imposed at the upper boundary of the computational domain. The convective and diffusive schemes are second-order accurate (Lax and Wendroff, 1960). The subgrid-scale viscosity is computed using the Sigma model (Nicoud et al., 2011). The mesh contains 0.9 million grid points. The mean first off-the-wall grid point height is $ {y}^{+}=41 $ in wall units, based on the WSS of the boundary layer upstream of the step to compute the wall-unit scaling, and in the range $ {y}^{+}=20 $ – $ 25 $ in the separated region downstream of the step. However, a specific refinement is performed farther from the wall, near the interface between the inlet boundary layer and the separated region. A close-up view of the mesh around the refined region is given in Figure 12.

Figure 12. Cross-section of the mesh around the refined region in the $ x $ – $ y $ plane.

Wall-modeled LESs performed using an algebraic shear stress model based on the law of the wall and GNN WSS models are compared in Figure 13. Simulations are performed with up to $ N=5 $ message-passing steps in order to study a posteriori the influence of the number of steps on the ability of the model to detect separation. The GNN WSS model is only applied on the bottom wall after the step to guarantee that all simulations have the same velocity profile upstream of the step. Overall, the mean velocity downstream of the step is mostly governed by the influence of the shear layer induced by the step. The influence of the GNN WSS model may nevertheless be observed near the onset of the mean backflow region ( $ x/{h}_s $ = 0–2). In particular, the presence of a region of positive streamwise velocity at the bottom corner of the step is predicted with $ N=2 $ or above, in consistency with the reference high-fidelity simulation, whereas it is absent in the simulation with an algebraic wall stress model. As may be seen in Figure 14, the velocity profiles from $ x/{h}_s=3 $ to further downstream are almost identical regardless of the wall-shear-stress model used.

Figure 13. A posteriori validation: Mean streamwise velocity in the flow over a backward-facing step according to (a) the reference wall-resolved simulation, (b) a large-eddy simulation with an algebraic wall stress model, and (c–g) graph neural network WSS models. The contour lines denote the level sets −4, −2, 0, 4, 8, 12, and 16 m/s. There are no values within the first large-eddy-simulation cell because the large-eddy simulations do not provide a physical velocity at the wall.

Figure 14. A posteriori validation: Profile of mean streamwise velocity in the flow over a backward-facing step at the location $ x/{h}_s=1 $ , 3, 7, 9, scaled by the free-stream velocity $ {u}_0 $ . The horizontal solid line gives the height of the first point off the wall.

The profile of mean WSS along the wall is reported in Figure 15. With the algebraic wall stress model, the predicted WSS is accurate near the step ( $ x/{h}_s<2 $ ) and far from the step ( $ x/{h}_s>12 $ ) but underestimated in the backflow and reattachment regions. The GNN WSS models improve the prediction in these non-equilibrium regions, by increasing the predicted WSS. The increase in WSS is negligible with $ N=1 $ but substantial already with $ N=2 $ . The improvements from $ N=3 $ to $ N=5 $ are rather small. Even with $ N=5 $ , the peak WSS remain underestimated compared to the reference high-fidelity simulation. It is difficult to disentangle the many possible sources of errors in LESs and determine the contribution of the wall-shear-stress model, the subgrid-scale model and the numerical scheme to the remaining discrepancy with the reference simulation. Note in addition that since the present shear-stress modeling paradigm only considers the magnitude of the WSS, the shear stress vector is assumed to be aligned with the tangential velocity at the first point above the wall. This assumption is unrealistic near the reattachment point, and leads to a significant underestimation of the reattachment point location in all LESs (Figure 16). The MSE between the shear stress of the LESs and the reference wall-resolved simulation is given in Table 5. The table summarizes the fact that, for this particular simulation, the model associated with $ N=4 $ message-passing steps leads to more accurate results than the algebraic model based on the law of the wall. Hence, these a posteriori results are consistent with the a priori tests of Section 4.1 and confirms that a GNN WSS model can lead to large improvements over a model based on the law of the wall for turbulent separated flows.

Figure 15. A posteriori validation: Mean wall shear stress profile in the flow over a backward-facing step.

Figure 16. A posteriori validation: Mean streamwise component of the wall shear stress vector in the flow over a backward-facing step. The numerical results of Le et al. (1997) are given for reference.

Table 5. A posteriori validation: Integral measures of the disagreement between the reference wall shear stress in the flow over a backward-facing step and the wall shear stress predicted by large-eddy simulations with an algebraic wall stress model and graph neural network WSS models.

4.2.4. Linear blade cascade

Finally, the simulation of a linear blade cascade has been performed with an algebraic model based on the law of the wall and a GNN WSS model. The configuration is similar to the simulation N65b described in Section 2.1 and Appendix A, which was not included in the training dataset. The blade geometry is a linearly extruded NACA 65-009 airfoil profile (Abbott et al., 1945), of chord $ c=0.15 $ m, attached to flat endwalls with turbulent incoming boundary layers. The geometric parameters of the blade and operating point conditions are reported in Table 6. The Reynolds number of this test case is representative of typical compressor blades. The flow is subsonic, with a low freestream Mach number. The incidence angle is $ I=7{}^{\circ} $ . The computational domain of the simulation is periodic in the pitchwise direction, implying that each blade is subject to the same flow field. In the spanwise direction, only half a blade has been simulated with a symmetry condition imposed at mid-span. These assumptions have been checked in the experimental campaigns and were shown to be sufficiently valid (Ma et al., 2011). We denote $ x $ the axial direction, $ y $ the pitchwise direction, and $ z $ the wall-normal direction. The origin of the coordinate system $ \left(x,y,z\right) $ is placed at the intersection of the leading edge of the blade and the endwall. A curvilinear coordinate $ {s}_1 $ is also defined on the blade surface, where $ {s}_1 $ is zero at the leading edge of the blade, positive on the suction side of the blade and negative on the pressure side of the blade. The inlet plane is located at $ {x}_{\mathrm{inlet}}/c=-3.65 $ and the outlet plane at $ {x}_{\mathrm{outlet}}/c=2.98 $ . The pitchwise length of the domain is $ s=0.89c $ . At the inlet, the experimental velocity profile of Zambonini et al. (2017) is imposed using the characteristic nonreflecting boundary condition of Poinsot and Lele (1992). The nonreflecting boundary condition of Granet et al. (2010) is used to impose a constant pressure of 1 standard atmosphere at the outlet. At the top boundary of the domain, located $ h=1.23c $ above the wall, there is a symmetric boundary condition which imposes the normal velocity and the normal gradients of velocity, temperature, and pressure to zero. A rectangular step of width 0.017 $ c $ and height 0.015 $ c $ , located 0.085 $ c $ after the inlet, triggers the transition of the inlet endwall boundary layer. A rectangular step is also placed on the suction side and pressure side of the blade surface, to reproduce the effect of the sandpapers used in the experiments and trigger turbulence transition. The blade steps are sized similarly to the sandpapers present in the experiments. However, the inlet step does not correspond to a feature present in the wind tunnel; the experimental incoming endwall turbulent boundary layer is originating from further upstream and hence cannot be reproduced numerically without simulating the entire wind tunnel.

Table 6. Geometric parameters and operating point conditions of the linear cascade (left) and schematic representation of the notations used (right).

Two LESs are performed: a simulation using an algebraic shear stress model based on the law of the wall and a simulation using the GNN WSS model. In both cases, a second-order accurate scheme is used for convection (Lax and Wendroff, 1960) and diffusion. The subgrid-scale viscosity is computed using the Sigma model (Nicoud et al., 2011). The mesh of the LESs is fully tetrahedral and contains 15 million grid points in total. The mesh is more refined around the leading edge in order to resolve the effect of the tripping steps (Figure 17). Outside this region, the mesh refinement near the blade is uniform. In wall units, the edge length varies in the range $ {\Delta}^{+}=40 $ – $ {\Delta}^{+}=90 $ on the blade surface, and accordingly the first off-the-wall grid point height on the blade surface varies from $ {y}^{+}=30 $ to $ {y}^{+}=67 $ due to the spatial variation of the WSS.

Figure 17. Cross-section of the mesh around the refined region in the $ x $ – $ z $ plane.

A comparison of the mean static pressure coefficient on the blade surface obtained in the reference wall-resolved simulation and in the LESs is given in Figure 18. The static pressure coefficient $ {C}_p $ is defined as

(21) $$ {C}_p=\frac{p-{p}_{\infty }}{p_{s,\infty }-{p}_{\infty }}, $$

where $ p $ is the local static pressure and $ {p}_{s,\infty } $ and $ {p}_{\infty } $ the reference inlet stagnation and static pressures, respectively. The profile of the pressure coefficient on the pressure side of the blade is well reproduced by the LESs both at $ z/h=5.4 $ %, near the endwall, and at $ z/h=46 $ %, near mid-height. At $ z/h=5.4 $ %, the plateau of static pressure coefficient in the aft part of the suction side, associated with the presence of a corner separation region, is shorter in the two LESs, in particular with the GNN WSS model. The contours of mean axial velocity (Figure 19) confirms that the corner separation region on the suction side occurs closer to the trailing edge in the two LESs than in the reference wall-resolved simulation. The coefficient of stagnation pressure loss $ {\omega}_s=\left({p}_{s,\infty }-{p}_s\right)/\left({p}_{s,\infty }-{p}_{\infty}\right) $ downstream of the trailing edge of the two LESs is compared to that of the reference wall-resolved simulation in Figure 20. In the two LESs, the wake has a wider footprint and the influence of the corner separation region does not extend as far away from the endwall.

Figure 18. A posteriori validation: Mean static pressure coefficient on the blade close to the endwall (a) and at mid-height (b). The upper branch corresponds to the pressure side of the blade and the bottom branch corresponds to the suction side of the blade.

Figure 19. A posteriori validation: Mean axial velocity on the plane $ z/h=3 $ % in the flow around the NACA 65-009 blade according to (a) the reference wall-resolved simulation, (b) a large-eddy simulation with an algebraic wall stress model, and (c) the graph neural network WSS model.

Figure 20. A posteriori validation: Stagnation pressure loss coefficient $ {\omega}_s $ at 36.3% axial chord downstream from trailing edge in (a) the reference wall-resolved simulation, (b) a large-eddy simulation with an algebraic wall stress model, and (c) a large-eddy simulation with the graph neural network WSS model.

The profile of mean WSS along the blade is reported in Figure 21. The mean WSS is overestimated by the LES with the algebraic wall stress model on the pressure side and on the suction side in the accelerating region near the leading edge. The GNN WSS model leads to a lower WSS amplitude in these regions, which leads to a clear improvement compared to the baseline algebraic model. However, the WSS is overestimated in the region of corner separation by the LES with the GNN WSS model, while the simulation with a law-of-the-wall algebraic model leads to a more accurate prediction in that region. This clearly demonstrates that a priori improvements do not always materialize into a posteriori improvements in all locations. This can be attributed to a discrepancy between the fields resulting from the filtering of direct numerical simulation data and the fields of LESs. It is difficult to assess whether this discrepancy is the responsibility of the wall-shear-stress model or other types of numerical or modeling errors. Overall, we thus consider that the generalizability of the model in the case of compressor blade cascades is still uncertain.

Figure 21. A posteriori validation: Mean wall shear stress profile on the blade surface, averaged in the spanwise direction.

5. Conclusion

A machine-learning wall-shear-stress model has been developed based on reference simulations of the flow in a smooth channel, in a 3DD and in a linear blade cascade. The model directly operates on the unstructured mesh of the simulation and is thus directly applicable in LESs of complex industrial flows. The relevance of the model has been verified a priori using simulations of the flow over a BFS not used during training and the simulation of a linear cascade at a larger incidence angle than used during training. The model showed a good agreement with the reference WSS and an ability to discriminate between the fully developed turbulent boundary layers, separated flow regions, and transitional boundary layers. The ability of the model to detect separated regions and improve the predictions of a model based on the law of the wall has been demonstrated in an a posteriori setting for the BFS simulation. Furthermore, the model was also able to match the performance of a model based on the law of the wall in a channel flow, for which it is well suited, although it is subject to the same classical issue of logarithmic-layer mismatch. The GNN WSS model also improves the WSS prediction on the pressure side of a blade in a linear cascade, and on the suction side of the blade in the region directly after a tripping step. However, the GNN WSS model is not strictly superior to an algebraic equilibrium wall model in that case, since the WSS is overestimated in the region of corner separation on the suction side of the blade. Overall, this paper demonstrates that graph-based models are a promising tool to address the problem of wall-shear-stress modeling. Improvements to the machine-learning methodology would require a better understanding of the transition from a priori tests to a posteriori tests, and of the relationship between the wall-shear-stress model, the subgrid-scale model and the numerical scheme. This study shows the value of data-driven efforts toward wall-shear-stress modeling in realistic solvers. This justifies further efforts to assemble even more comprehensive and varied databases, which are critical to produce more accurate and robust models.