1. Introduction
Simulation of turbulent flows remains a challenge because of the disparate range of spatial and temporal scales that need to be resolved (Pope Reference Pope2000). An alternate to directly solving the Navier–Stokes equations is to solve its reduced complexity versions. Reynolds-averaged Navier–Stokes (RANS) models solve for the ensemble average or time-average of the true solution. Large eddy simulations (LES) (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Meneveau, Lund & Cabot Reference Meneveau, Lund and Cabot1996; Nicoud & Ducros Reference Nicoud and Ducros1999; Codina Reference Codina2002; Vreman Reference Vreman2004; Bazilevs et al. Reference Bazilevs, Calo, Cottrell, Hughes, Reali and Scovazzi2007; Codina et al. Reference Codina, Principe, Guasch and Badia2007; You & Moin Reference You and Moin2007; Gravemeier et al. Reference Gravemeier, Gee, Kronbichler and Wall2010; Wang & Oberai Reference Wang and Oberai2010; Masud & Calderer Reference Masud and Calderer2011; Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011; Parish & Duraisamy Reference Parish and Duraisamy2017) resolve the spatiotemporal dynamics of the large scales. The cost of LES is, however, still prohibitive near the wall. To alleviate the need of mesh refinement near the wall, boundary conditions are imposed weakly in a wall-modelled (WM)LES (Piomelli & Balaras Reference Piomelli and Balaras2002; Bose & Moin Reference Bose and Moin2014; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019). Alternate approaches to WMLES are hybrid RANS–LES (HRLES) techniques (Strelets Reference Strelets2001; Menter & Egorov Reference Menter and Egorov2005; Fröhlich & Von Terzi Reference Fröhlich and Von Terzi2008; Sagaut, Terracol & Deck Reference Sagaut, Terracol and Deck2013; Shur et al. Reference Shur, Spalart, Strelets and Travin2015; Menter Reference Menter2016) such as detached eddy simulation (DES) (Spalart Reference Spalart2009) and improved delayed detached eddy simulation (IDDES) (Shur et al. Reference Shur, Spalart, Strelets and Travin1999), where the inner layer is solved using RANS and the rest using LES. As a consequence, the cost associated with resolving the near-wall structures in the stream-wise and span-wise direction is no longer present. The cost associated with resolving the wall-normal gradient is still present in the HRLES approaches.
Over the past few decades, various contributions have been made in the development and application of these methods to highly complex problems (e.g. Park & Moin Reference Park and Moin2016; Goc, Bose & Moin Reference Goc, Bose and Moin2020; Iyer & Malik Reference Iyer and Malik2020; Lozano-Duran, Bose & Moin Reference Lozano-Duran, Bose and Moin2020; Goc et al. Reference Goc, Lehmkuhl, Park, Bose and Moin2021; Kiris et al. Reference Kiris, Ghate, Duensing, Browne, Housman, Stich, Kenway, Dos Santos Fernandes and Machado2022). Our view is that because all of the scale-resolving methods are coarse-grained from the Navier–Stokes equations, there must exist a unified view. The partially-averaged Navier Stokes (PANS) approach brings together several turbulence closures of various modelled-to-resolved scale ratios ranging from RANS to Navier–Stokes (direct numerical simulations (DNS)) into one formulation (Girimaji & Abdol-Hamid Reference Girimaji and Abdol-Hamid2005). The behaviour of the PANS equations can be varied smoothly from the RANS equations to the Navier–Stokes (DNS) equations by changing the filter-width control parameters. The unified RANS–LES approach (Heinz Reference Heinz2007; Gopalan, Heinz & Stöllinger Reference Gopalan, Heinz and Stöllinger2013) is an optimal HRLES framework that uses different time scales to switch between the RANS and LES approaches. In pursuit of similar unified models and in an effort to augment existing frameworks, we propose a filtering technique using optimal finite-element projections which: (i) offers a unifying perspective through a common coarse-graining strategy; (ii) provides optimal solutions for the existing coarse-grained methods to improve upon.
The use of filtered DNS data to perform a priori analysis of closure models for RANS and LES is indeed not new. LES models such as the scale similarity or Smagorinsky models have also been frequently evaluated against sub-grid stresses obtained from filtered DNS data (Vreman, Geurts & Kuerten Reference Vreman, Geurts and Kuerten1995; Meneveau & Katz Reference Meneveau and Katz2000; Girimaji & Abdol-Hamid Reference Girimaji and Abdol-Hamid2005; Bou-Zeid et al. Reference Bou-Zeid, Vercauteren, Parlange and Meneveau2008). In most prior studies, filtering is either performed in the Fourier space using the sharp spectral cutoff or Gaussian filters when the problem has periodic directions or the box filter in more complex problems. In case of filters that are applied in spectral space, the filter width remains the same along the periodic directions in which it is applied. However, as observed in most coarse-grained simulations, the filter width can vary considerably. In fact, filter sizes define these methods. For example, in case of a LES of channel flow that is performed on a structured grid, the filter size in the span-wise and stream-wise directions scale with the wall units and can be a constant. However, the filter width in the wall-normal direction can vary from a few wall units near the wall to 
$0.1\delta$ at the centre of the channel or the edge of the boundary layer. In the traditional WMLES, the filter width is approximately of the order of 
$0.1\delta$ throughout in all directions. For HRLES approaches (such as DES and IDDES), the filter width is of the order of 
$0.1\delta$ in the span-wise and stream-wise directions, and similar to LES in the wall-normal direction. The non-uniform filtering requirement in LES and HRLES methods stems from the fact that in both the cases, the wall stress is resolved which requires a near-wall grid that scales with wall units. In addition to the filter size, the type of filter can also change the nature of the solution. For example, both the box and spectrally (sharp) filtered DNS both qualify as synthetic LES solutions. In the case of finite-element projections, the quality of the filter is linked to the order of polynomial used to filter the solution. In this work, we aim to address some of these issues by using finite-element projections which allow for variation in the filter width in the domain and also provide the required flexibility to change the quality of the filter by changing the order of the polynomial.
The idea of projection is at the core of the variational multiscale method (VMS) (Hughes et al. Reference Hughes, Feijóo, Mazzei and Quincy1998; Codina Reference Codina2002; Bazilevs et al. Reference Bazilevs, Calo, Cottrell, Hughes, Reali and Scovazzi2007; Codina et al. Reference Codina, Principe, Guasch and Badia2007; Gravemeier et al. Reference Gravemeier, Gee, Kronbichler and Wall2010; Wang & Oberai Reference Wang and Oberai2010; Masud & Calderer Reference Masud and Calderer2011; Parish & Duraisamy Reference Parish and Duraisamy2017). In VMS, projections are used to formally distinguish the coarse scales from the fine scales. The coarse-scale (filtered) solution that is obtained after the projection operation represents the ‘best’ coarse-grained solution 
$u$ on the coarse space based on some optimality condition, for example, the 
$L_2$-optimality condition. The sharp spectral filter obtained by truncation in Fourier space is also based on the idea of 
$L_2$-projection onto the Fourier basis functions. In this work, we perform 
$L_2$-projections on finite-element basis functions. It is also pertinent to mention that the current idea of optimal projections should not be confused with optimal LES (Langford & Moser Reference Langford and Moser1999). Our work optimally represents the DNS solution 
$u$ on a finite-dimensional coarse space, whereas optimal LES is an ideal LES model that targets accurate single-time multi-point statistics of the coarse solution.
Projected DNS data have been previously used to improve both existing finite-element methods (Pradhan & Duraisamy Reference Pradhan and Duraisamy2021) and turbulence models (Vreman et al. Reference Vreman, Geurts and Kuerten1995; Meneveau & Katz Reference Meneveau and Katz2000; Girimaji & Abdol-Hamid Reference Girimaji and Abdol-Hamid2005; Bou-Zeid et al. Reference Bou-Zeid, Vercauteren, Parlange and Meneveau2008). An important question is whether the applicability of the present filtering method is only restricted to the finite-element method because the functions that are used for projection are the finite-element basis functions. In this paper, however, we employ them as an alternative to the traditional filters for assessing all kinds of methods and not just finite-element methods. As discussed previously, the present approach has advantages in cases where the filter length is anisotropic, varies rapidly or when non-homogeneous directions are present, as in wall-bounded flows. The filtering strategy that has some similarities to the present approach is the differential filter (Germano Reference Germano1986; Najafi-Yazdi, Najafi-Yazdi & Mongeau Reference Najafi-Yazdi, Najafi-Yazdi and Mongeau2015), which consists of a filtering length scale 
$l_p$. This length scale 
$l_p$ can be varied along the domain to have a similar effect.
In the past few years, several efforts have been made to train sub-grid models using machine learning approaches both in an offline and model-consistent setting (Sarghini, De Felice & Santini Reference Sarghini, De Felice and Santini2003; Gamahara & Hattori Reference Gamahara and Hattori2017; Maulik & San Reference Maulik and San2017; Maulik et al. Reference Maulik, San, Rasheed and Vedula2018; Wang et al. Reference Wang, Luo, Li, Tan and Fan2018; Beck, Flad & Munz Reference Beck, Flad and Munz2019; Maulik et al. Reference Maulik, San, Jacob and Crick2019; Xie et al. Reference Xie, Li, Ma and Wang2019a,Reference Xie, Wang, Li, Wan and Chenb; Xie, Wang & Weinan Reference Xie, Wang and Weinan2020). Such data-driven LES models require a filtered form of the DNS, which, in turn, will depend on the filter size. Similarly, in WMLES, the model cannot be trained using the mean solution in the first few grid points where the influence of the slip condition will be observed. By applying projections and obtaining statistics from the optimal solution, one can obtain more reasonable targets to training the model (Beck et al. Reference Beck, Flad and Munz2019; Duraisamy, Iaccarino & Xiao Reference Duraisamy, Iaccarino and Xiao2019; Chung & Freund Reference Chung and Freund2022) on more complex problems. Note that this is a first step towards addressing model consistency (Duraisamy Reference Duraisamy2021). Finite-element projections are not restricted to simple geometries and can be applied to more complex flows, and provide the additional flexibility of choosing polynomial orders for the geometry and the solution independently.
The main objective of the present work is to develop a unified framework that can be used as a lens to quantitatively assess, augment and calibrate a wide range of coarse-grained models. Particular attention is paid to the behaviour of various models in the proximity of the wall, and to ascertain whether scaling relationships exist.
In § 2 of this paper, we describe the procedure of performing 
$L_2$-projection and provide a discussion on the choice of coarse basis functions that will be used. In § 3, we compute filtered solutions for the channel flow problem in the LES, WMLES and HRLES limit, and compute the coarse-scale statistics in each of these cases. In § 4, we show that the slip velocity in case of WMLES is a natural consequence of under-resolution in the wall-normal directions and guiding principles for improved slip wall models are proposed. In § 5, we propose new slip-wall-based wall model forms and evaluate its performance in comparison with traditional WMLES. Perspectives on improved slip-wall models are provided in § 6. Finally, we conclude the work in § 7.
2. Finite-element projection
The goal of this section is to construct a generalised filtering approach that can be used to assess various coarse-grained simulations in problems using localised bases applicable to non-periodic boundary conditions. Further, the meshes can contain anisotropic elements along with the possibility of grid stretching. As a first step, however, we obtain high-resolution data for filtering. For our purpose, we use the channel flow DNS data at friction Reynolds numbers of 
$Re_{\tau }\approx 1000$ and 
$Re_{\tau }\approx 5200$ from the Johns Hopkins Turbulence Database (JHTDB) (Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008; Lee & Moser Reference Lee and Moser2015) and a smaller channel 
$Re_{\tau }\approx 950$ case form the Texas turbulence file server (Hoyas & Jiménez Reference Hoyas and Jiménez2006). To this end, consider the decomposition of the full-order (DNS) solution 
$u$ into coarse and fine scales as
\begin{equation} u = {u_h} + u', \end{equation}
where 
${u_h} \in {\mathcal {V}_h}$ and 
$u' \in \mathcal {V}'$ as shown in figure 1. The vector space of functions 
$\mathcal {V}\equiv {L}^2(\varOmega )$ is the space of square-integrable functions. This space is decomposed as
\begin{equation} \mathcal{V} = {\mathcal{V}_h} \oplus \mathcal{V}', \end{equation}
where 
$\oplus$ represents a direct sum of 
${\mathcal {V}_h}$ and 
$\mathcal {V}'$. Let us also define 
$\mathcal {T}_{h}$ to be a tessellation of domain 
$\varOmega$ into a set of non-overlapping elements, 
$K$, each having a sub-domain 
$\varOmega _K$ and boundary 
$\varGamma _K$. The functional space 
$\mathcal {V}$ is infinite dimensional and must be approximated by a finite-dimensional approximation 
${\mathcal {V}_h}$. The domain and boundary of an element marked by 
$\varOmega _e$ and 
$\varGamma _e$, respectively. In the case of the continuous Galerkin (CG) method, the coarse space basis functions 
${\mathcal {V}_h} \subset C^{0} \cap {L}^2(\varOmega )$ have 
$C^{0}$ continuity everywhere including element boundaries. In the case of the discontinuous Galerkin (DG) methods, the coarse space 
${\mathcal {V}_h}$ is defined as
\begin{equation} {\mathcal{V}_h} \triangleq \{u\in{L_2}(\varOmega):u|_{K}\in P^k(K),K\in\mathcal{T}_{h} \}, \end{equation}
where the space of polynomials up to degree 
$k$ is denoted as 
$P^k$. Defining 
${\mathcal {V}_h}$ in this manner allows for discontinuities in the solution across element boundaries. The DG space is a more richer space in comparision to the space if both the number of elements and polynomial order are kept fixed. Irrespective of the choice of basis functions used (CG or DG), given 
$u$ from the high-fidelity simulation, our goal is to find the optimal representation of 
$u$ in the coarse sub-space 
${\mathcal {V}_h}$. In our case, we use the 
${L}^2$-projection to obtain 
$u_h$ which minimises the value of 
$\|u-u_h\|_2^2$. This problem is equivalent to the problem of finding 
$u_h\in {\mathcal {V}_h}$ such that
\begin{equation} (u,w_h)=(u_h,w_h) \quad \forall \, {w_h} \in {\mathcal{V}_h}. \end{equation}
where 
$(\cdot,\cdot )$ denotes the 
$L_2$ inner product and 
$w_h$ denotes a coarse-space weighting function. In the case of CG basis functions, the mass matrix is global and a large matrix needs to be inverted to obtain the final filtered solution. The DG mass matrix on the other hand is local to the element and lends itself to easy parallelisation.
Schematic of the projection of the DNS solution 
$u$ on coarse finite-element spaces 
$\mathcal {V}_h$ to obtain the 
$L_2$-optimal LES, WMLES or HRLES solution 
$u_h$.
As 
$u_h$ and 
$w_h$ are finite dimensional, their inner product can be computed precisely using quadrature rules. Here 
$(u,w_h)$ requires special care because 
$u$ is extremely high dimensional in comparison with 
$w_h$. The high-dimensionality of 
$u$ is restricted by the size of DNS which exists on a very fine mesh capable of resolving the Kolmogorov scales 
$O(\eta )$. To compute this term precisely, we interpolate the coarse-scale basis functions and the DNS solution on a very fine mesh of the size of 
$O(\eta )$ and apply numerical integration to compute the inner products. The size of the numerical integration mesh is adjusted until the final projected solution is independent of the numerical integration mesh size. Additional details on the method to compute the 
$L_2$-procedure are given in Appendix A.
The final comment is on the imposition of the near-wall behaviour of the coarse space. There are two choices: (i) project on a space which strongly satisfies the boundary condition at the nodal points; or (ii) keep the boundary degrees-of-freedom (DOFs) free and make no such assumptions. The second choice appears more reasonable because when the solution is coarse grained in the wall-normal direction, the solution might no longer satisfy the no-slip boundary conditions strongly. This is especially true for WMLES where the coarse-grained solution no longer satisfies the no-slip boundary condition and slip is observed. However, as the grid is refined near the wall, the no-slip boundary condition is naturally satisfied.
3. Application to channel flow
As a first step towards obtaining the projected DNS solution for the channel flow problem, we discuss the effect of the choice of the coarse-space basis functions on the coarse-scale solution obtained after the projection operation. Depending on the coarse-space basis, the projected solution can be either a low-dimensional compressed representation of the original solution or a spatially filtered version of it. The low-dimensional compressed representation is obtained when the coarse basis is tailored using data or existing analytical solutions. To ensure that the projection step leads to a more general spatial filtering approach, non-tailored basis functions commonly used in the finite-element method are used. The projection operation onto these coarse finite-element grids will lead to filtering.
The resulting coarse solution after filtering might be considerably different from the DNS solution due to truncation of the high-frequency components present originally in the DNS solution. In the near-wall region, the effect of projection can vary with the size of the filter in each direction. One manifestation of under-resolution in the wall-normal direction is the occurrence of a slip velocity. This slip-velocity can, in fact, be tracked down to the mean profile itself. To resolve the mean solution, a near-wall grid spacing of 
$\Delta y^+ \approx 1$ is required in the wall-normal direction. However, if a grid size of 
$\Delta y \approx 0.1\delta$ is used, even the mean solution can no longer be resolved and a slip velocity at the wall will be observed. This is true unless the solution is artificially forced to go from a large value to zero over just one grid point. A solution to make the coarse-scale solution satisfy the no-slip boundary condition is to enrich the coarse space with a tailored basis (Krank & Wall Reference Krank and Wall2016). As a consequence, the tailored basis mimics the mean profile between the wall and the first grid point and ensures that the no-slip is satisfied collectively by the coarse non-tailored basis and the enriched tailored basis.
To define the coarse space, we first construct a finite-element mesh and chose the polynomial order of the basis functions. The idea here is that by selecting the grid and the polynomial order of the basis functions, we are enforcing our desired filter size distribution. A variety of coarse spaces have been generated as listed in table 1. The ‘A’-type grids are the DNS grids on which the high-resolution solution 
$u$ exists. Two different DNS solutions at friction Reynolds numbers of 
$Re_\tau \approx 950$ and 
$Re_\tau \approx 1000$ are used and their corresponding grids are marked as 
$A1$ and 
$A2$, respectively. The 
$Re_\tau \approx 950$ solution (Del Alamo et al. Reference Del Alamo, Jiménez, Zandonade and Moser2004; Hoyas & Jiménez Reference Hoyas and Jiménez2006) is obtained from a relatively smaller domain having a stream-wise size of 
$L_x \approx 2{\rm \pi} \delta$ and a span-wise sizes of 
$L_z \approx {\rm \pi}\delta$, whereas the 
$Re_\tau \approx 1000$ solution (Perlman et al. Reference Perlman, Burns, Li and Meneveau2007; Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008; Lee & Moser Reference Lee and Moser2015) is obtained as a cutout from a simulation performed on a larger domain.
Summary of mesh parameters. Here, 
$\varDelta _x^+$, 
$\varDelta _y^+$ and 
$\varDelta _z^+$ are the effective grid sizes in different directions 
$\varDelta _x$, 
$\varDelta _y$ and 
$\varDelta _z$ normalised with wall units, 
$\delta$ is the half-channel height, 
$N_x$, 
$N_x$ and 
$N_z$ represent the number of DOFs in the stream-wise, wall-normal and span-wise directions respectively, 
$p$ is order of polynomial used, 
$SR$ is the stretching ratio used to generate the grid. The effective grid sizes 
$\varDelta _x$, 
$\varDelta _y$ and 
$\varDelta _z$ for the finite-element grid are defined as 
$\varDelta _x = \varDelta ^e_x/p$, 
$\varDelta_y = \varDelta_y^e/p$ and 
$\varDelta_z = \varDelta_z^e/p$, respectively. The quantities 
$\varDelta ^e_x$, 
$\varDelta ^e_y$ and 
$\varDelta ^e_z$ represent the actual element sizes in the finite-element mesh.
The ‘B’-type grids, on the other hand, are tailored for performing wall-resolved LES (figure 2). As a result of the size of the largest energy containing eddies scaling with the distance from the wall (Yang & Griffin Reference Yang and Griffin2021) outside the viscous sub-layer, a mesh resolution of 
$\Delta y \approx 0.1\delta$–
$0.25 \delta$ is used for the ‘B’-type grids at the centre of the channel. Similarly, the ‘C’-type grids are tailored for performing WMLES simulations using the wall-stress- or the slip-wall-based approaches. The ‘D’-type grids are more suitable for assessing the WMLES branch of the HRLES methods (Shur et al. Reference Shur, Spalart, Strelets and Travin2008). For ‘D’-type, the resolution in the stream-wise and the span-wise direction is similar to the ‘C’-type grid. However, in the wall-normal direction, a grid spacing similar to ‘B’-type grid has been assumed, i.e. 
$\Delta y^+ \approx 0.1$–1 in the near-wall region and 
$\Delta y \approx 0.1\delta$–
$0.25 \delta$ at the centre of the channel. The ‘E’-type grid is an extremely coarse grid with resolutions of 
$\Delta x \approx 0.35\delta$, 
$\Delta y \approx 0.334\delta$ and 
$\Delta z \approx 0.35\delta$ in the stream-wise, span-wise and the wall-normal directions, respectively. For all the grid types, the mesh is uniform in the stream-wise and the span-wise directions. However, in the wall-normal direction, the mesh has been stretched geometrically for cases ‘B’ and ‘D’. In case of ‘C’- and ‘E’-type grids, uniform mesh is assumed in the wall-normal directions as well. For each type of grid, two different polynomial orders 
$p=1,2$ are used to construct the projection coarse-space. The stretch rates (SRs) and the polynomial orders for different cases have been summarised in table 1.
Near-wall grids used for cases B, C, D and E.
Figure 3 shows the mean and second-order statistics computed using the projected solution for the cases B1, C1, D1 and E1. The goal is to compare the optimal solutions for different grids that correspond to different coarse-grained approaches, except for E1, which is an extremely coarse mesh and is not suitable for any existing method. From figure 3(a), it can be observed that the mean velocity is well-resolved for cases B1 and D1. For case C1, which represents an optimal WMLES solution, the mean velocity is well-resolved only after the first grid point, i.e. 
$y/\delta > 0.05$. As can be observed in figure 3(a), case E1 is extremely coarse and fails to resolve the mean velocity until the outer limit of the log-layer is reached. For these cases with wall-normal under-resolution, the effect of under-resolution results in slip velocity 
$u_s$ at the wall, which can be calculated by evaluating 
$u_h$ at the wall. The magnitude of the mean stream-wise slip velocity 
$\langle u_s \rangle ^+$ was found to increase with the under-resolution, i.e. 
$\langle u_s \rangle ^+ \approx 6$ for case C1 to 
$\langle u_s \rangle ^+ \approx 12$ for case E1. It can also be observed that all the methods except E1 resolve the second-order statistics outside the inner layer. Inside the near-wall region, only B1 is capable of accurately resolving the turbulence stresses. Among the second-order statistics, the effect of filtering is most strongly felt on the wall-normal fluctuations. It can be observed for cases C1 and D1 that the wall-normal fluctuations far away from the inner layer are under-represented even when stream-wise and the span-wise fluctuations as close to the DNS solution.
Comparison of mean and second-order statistics for cases B1, C1, D1 and E1. The symbols in all the plots correspond to the value at the nodal point. In (a), the solutions for cases C1 (WMLES) and E1 (extremely coarse), are interpolated to the DNS mesh using the coarse finite element basis functions near the wall to show the slip effects.
Figure 4 shows the stream-wise velocity energy spectra of the projected solution in the span-wise and stream-wise directions for cases A2, B4, C2 and D4 at two different wall-normal locations. The choice of cases plotted here is based on the most suitable mesh sizes for the various methods. As can be seen from figures 4(a) and 4(b), the large scales are well-represented at the centre of the channel (
${y}/{\delta } \approx 1.0$) both in the stream-wise and the span-wise directions, for all the methods. However, as can be seen in figure 4(c), the large scales are not represented accurately by the C2 (WMLES) and D4 (HRLES) cases in the near-wall region (
$y^+ \approx 15$). However, in the span-wise direction (figure 4d), the large scales are relatively well represented in D4 in comparison with C2. In the following section, the effect of the projection is discussed individually for each type of mesh.
Comparison of the energy spectra for the stream-wise velocity component at the near-wall region and the centre of channel for cases A2, B4, C2 and D4.
3.1. The wall-resolved LES limit
Coarse-scale statistics for B2, B3 and B4 are provided in the top row of figure 5. It can be observed that all the cases perform well in resolving the mean profile. Similar trends are observed in figures 5(b) and 5(c), where the cases B3 and B4 only slightly outperform the case B2 in resolving the second-order statistics. This suggests that the sensitivity of the coarse-scale statistics to the wall-normal SR is not as high as the sensitivity to the mesh resolution in the stream-wise and the span-wise directions.
Comparison of optimal projections: (a–c) wall-resolved cases; (d–f) WMLES; (g–i) HRLES.
3.2. The WMLES limit
Coarse-scale statistics for C2, C3 and C4 are provided in the middle row of figure 5. In cases C2, C3 and C4, the mesh resolutions in wall units is much larger than unity and cannot resolve the mean velocity profile accurately near the wall. As a consequence, a slip velocity 
$u_s$ is observed for all the cases after the projection step. It was also observed that the magnitude of the mean slip velocity 
$\langle u_s \rangle ^+$ is highest for the case with the maximum wall-normal under-resolution, i.e. C2. This magnitude goes down as the wall-normal mesh is refined from C2 to C3 or C2 to C4. As expected, the second-order coarse-scale statistics are only accurate outside the inner layer towards the centre of the channel. It can also be observed that the stream-wise and span-wise velocity fluctuations computed using the coarse solution do not go to zero near the wall. The wall-normal velocity fluctuations, however, go to zero at the wall. This near-wall behaviour of the coarse scales is consistent with existing WMLES simulations in the literature, e.g. Wang, Hu & Zheng (Reference Wang, Hu and Zheng2020) and Kawai & Larsson (Reference Kawai and Larsson2012).
3.3. The HRLES limit
Cases D2, D3 and D4 represent the grids for HRLES methods such as IDDES (Shur et al. Reference Shur, Spalart, Strelets and Travin2008). Coarse-scale statistics for D2, D3 and D4 are provided in the bottom row of figure 5. The mean profile is resolved accurately. The second-order statistics in the region outside of 
$10$–
$20\,\%$ of the boundary layer is almost identical to that of the ‘C’-type grids in the middle row. However, near the wall, unlike the ‘C’-type WMLES cases, all the velocity fluctuations go to zero due to the no-slip condition being satisfied by the HRLES cases, and they are under-represented when compared with ‘B’-type LES cases. These observations are consistent with previous results from the literature (Friess & Davidson Reference Friess and Davidson2020).
In this section, we have used the labels LES, WMLES and HRLES to distinguish between the various methods. However, the difference between the latter two is subtle. The traditional WMLES method essentially uses RANS knowledge to compute the wall stress at the wall and can also be called a HRLES approach. Similarly, the HRLES approaches (such as IDDES) by virtue of solving the RANS equations near the wall reduce the computational cost associated with resolving the wall and can also be considered a WMLES. However, the context in which the labels WMLES and HRLES have been used in this paper is based on whether these models integrate to the wall or not inside a single domain, in other words if the size of the filter in the wall-normal direction is large or not in wall units.
It is worthwhile to mention that the results presented in this section appear to be more accurate in comparison with those in the literature. For instance, a large stretching ratio of 1.33 has been used for some of the meshes which do not induce a significant error in the filtered solution, however, this SR is more than the suggested limit for many methods. In addition, no log-layer mismatch (LLM) was observed in any of these cases. The resolution considered here for WMLES and HRLES is of the order of 10 points per semi-channel height 
$\delta$ in the stream-wise and the span-wise directions, which is coarse compared with the guidelines for these approaches. A comparison between the optimal solution for the C3 case and a WMLES solution using the traditional wall-stress-based approach computed on the same grid is presented in figure 6. The mean solution for both cases begins to deviate from the DNS at similar locations. The resolved turbulent shear stress is under-represented near the wall, starting almost identically, however, differing in their peaks. Similarly, the wall-normal velocity fluctuations are almost identical. On the other hand, the velocity inside the first element is slightly under-predicted in the traditional WMLES approach in comparison to projected DNS. In addition, both the stream-wise and span-wise velocity fluctuations reveal an overshoot near the second off-wall grid point in the traditional WMLES method.
Comparison of projected DNS and solution from the traditional WMLES method at similar resolutions.
One possible reason for the discrepancies between the results from the true simulation and the filtered DNS is the lack of accurate closures. With a poor model, the LLM may persist until a DNS-like resolution is reached. This is, in principle, similar to attributing inaccuracies in a wall-resolved LES with a standard Smagorinsky model when a dynamic Smagorinsky model might yield near-optimal performance. However, using the optimal projection framework presented here, it is now possible to perform an analysis of the closure terms and evaluate modelling errors. By reducing the modelling errors, the goal is to force the solution to reach a near-optimal state. Ideally, we would have wanted to improve all three approaches using our optimal projection framework. However, to have a compact presentation, we only consider evaluating the modelling errors in the slip-wall-based WMLES models and improve its a posteriori performance.
4. Analysis of slip-based wall models
The slip velocity at the wall in WMLES is related to under-resolution in the wall-normal direction. In this section, we seek to quantify this slip velocity to ensure that the resulting model generalises well to different Reynolds numbers. To understand the Reynolds number dependence, DNS from two different friction Reynolds numbers of 
$Re_{\tau } \approx 1000$ and 
$Re_{\tau } \approx 5200$ are used. As mentioned earlier, the computation of the 3-D projection of the 
$Re_{\tau } \approx 5200$ case by sequential 1-D projections in the wall-normal, stream-wise and span-wise directions is computationally expensive. To ensure computational efficiency and utilising the fact that this is a near-wall phenomenon, we project the DNS solution on uniform elements of size 
$\varDelta _e$ with polynomial basis functions. This is equivalent to performing a full 3-D projection on a DG finite-element solution space. The element shares its bottom face with the wall of the channel to mimic a near-wall grid. By moving the position of this element on the wall surface, different realisations of the slip velocity and coarse solution gradients in the wall-normal direction can be obtained. This is possible due to the statistical homogeneity present in the stream-wise and span-wise directions.
For each realisation, the projection of the DNS solution on the finite-dimensional DG space leads to filtering of the DNS solution as shown in figure 7. Figure 7(a) shows the contour of the DNS solution of the stream-wise velocity component for a sample 3-D element. The projected DNS solution for the same element is shown in figure 7(b). The projected DNS solution does not satisfy the no-slip boundary condition at the wall and does not contain the fine-scale information present in the original DNS solution. The goal is to assess the slip-wall-based wall model proposed by Bose & Moin (Reference Bose and Moin2014) and Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019)
\begin{equation} u_s = C_w \Delta \frac{\partial u_h}{\partial n},\end{equation}
and obtain an estimate of the model coefficient 
$C_w$. The coarse field 
$u_h$ can be obtained by either projecting the stream-wise, the span-wise or the wall-normal velocity fields. As a result, the value of 
$C_w$ obtained is tied to the velocity component that is used for projection. The computation of 
$C_w$ using (4.1) requires the computation of 
$u_s$ and the pre-multiplied gradient 
$\Delta ({\partial u_h}/{\partial n})$. The slip velocity 
$u_s$ is obtained by evaluating the coarse-scale solution at the wall as shown in figure 7(c). The pre-multiplied gradient 
$\Delta ({\partial u_h}/{\partial n})$ is obtained by computing the derivatives of the coarse scale in the wall-normal direction and multiplying with the normalised resolution 
$\varDelta = {{\varDelta _e}/{p}}$ as shown in figure 7(d). However, this results in an over-determined system for 
$C_w$. In general, the value of 
$C_w$ is also expected to change with the filter size 
$\varDelta$ used for the projection operation. The problem of this system being over-determined is solved by performing a least-squares minimisation over many such realisations until convergence in the estimates of 
$C_w$ was obtained. To solve the problem of the model coefficient 
$C_w$ being dependent on the filter size 
$\varDelta$, we perform dimensional analysis. Other parameters that could affect 
$C_w$ are (a) the order of polynomial used for projection 
$p$, the viscosity 
$\nu$ and (b) the wall stress 
$\tau _w$. After non-dimensionalisation, the following model form for 
$C_w$ can be obtained: 
$C_w = g_{p}(\varDelta ^+),$ where 
$g_{p}$ is a function of the grid resolution normalised with wall units 
$\varDelta ^+$ and the subscript 
$p$ denotes the coarse space polynomial order used for projection. The parameter 
$\varDelta ^+$ can be considered to be an indicator of the near-wall grid resolution. Similarly, the order of the numerical method can be encoded in 
$p$. Higher 
$p$ implies that a more accurate numerical method has been used to compute the LES solution. However, this implies that for every polynomial order 
$p$ we have to learn a new function. In addition, the numerical methods used to perform LES might work sub-optimally and the exact order might not be preserved. Hence, it is necessary that the effect of the numerical method be parameterised through a model constant similar to the Smagorinsky model coefficient 
$C_s$.
Filtered solution inside an sample element obtained by 3-D projection of near-wall 
$Re_{\tau } \approx 5200$ channel data.
Before investigating the slip velocity due to the full 3-D projection of the DNS solution, it is important to consider the contribution from the mean stream-wise velocity profile itself. As a first step, we apply 1-D projection to the Reichardt profile (Reichardt Reference Reichardt1951) which describes the mean profile in the inner layer. As the mean solution is invariant in the stream-wise and span-wise direction for a channel, the 3-D projection is reduced to a 1-D projection in the wall-normal direction only. Figure 8(a) shows the estimates for 
$C_w$ obtained from the Reichardt profile for different orders of projection. It can be observed that the 
$C_w$ profiles for different orders are distinct even after normalisation of the element size 
$\varDelta _e$ by 
$p$ to obtain 
$\varDelta = {\varDelta _e}/{p}$.
(a) The model coefficient 
$C_w$ computed using Reichardt profile by projecting on different polynomial basis; (b) 
$\lambda$-normalised version of Reichardt profiles; (c) 
$\lambda$-normalised 
$C_w$, i.e. 
$C_{w,\lambda }$ computed by 3-D projection of DNS on different polynomial spaces compared with 1-D projection of the Reichardt profile.
Inspired by the Smagorinsky model, which consists of a model constant 
$C_s$ that pre-multiplies the grid size in the final model form, we introduce a new model constant 
$\lambda$ which in its inverted form, i.e. 
${1}/{\lambda }$ pre-multiplies the grid-size in our proposed model. Figure 8(b) shows the estimated valued for 
$C_{w,\lambda }$ for different polynomial orders along with the 
$\lambda$ values for which all the curves collapse to the 
$p=1$ curve with 
$\lambda = 1$. As a result, it is possible to learn just one curve and parameterise it with an additional factor 
$\lambda$ to obtain the 
$C_w$ curves for different cases, i.e. 
$C_w = {C_{w,\lambda }}/{\lambda } = g_{1}(\varDelta ^+/\lambda )/{\lambda },$ where 
$C_{w,\lambda }= {\lambda p u_s}/{\varDelta _e ({\partial u_h}/{\partial n})}$. It is also important to check if similar relations also hold true for the 3-D projected solution.
Figure 8(c) compares the 
$C_{w,\lambda }$ curves obtained through the 3-D projection of DNS solutions to that obtained using the Reichardt mean profile for two different projection orders 
$p=1,3$. To obtain these curves, large variations in the element sizes have been considered. We use element sizes with 
$\varDelta _e \approx 0.011\delta$–
$0.28 \delta$ for projecting the 
$Re_{\tau } \approx 1000$ data and element sizes with 
$\varDelta _e \approx 0.033\delta$–
$0.30 \delta$ for projecting the 
$Re_{\tau } \approx 5200$ data. The effective filter sizes corresponding to these grids can be approximated by normalising the element size with 
$p$ to obtain 
$\varDelta = {\varDelta _e}/{p}$ . By re-using the 
$\lambda$ values from figure 8(b), similar collapse in the 
$C_{w,\lambda }$ curves was also obtained for the 3-D projection cases for the stream-wise and the span-wise velocity components. For each projection order, results for two different friction Reynolds numbers of 
$Re_\tau \approx 1000$ and 
$Re_\tau \approx 5200$ are plotted. The results indicate that for different polynomial orders, the 
$C_{w,\lambda }$ estimates for different 
$Re_{\tau }$ at a particular 
$\varDelta ^+$ are same, suggesting that 
$C_{w,\lambda }$ is a universal function of 
$\varDelta ^+$. Further, the 
$C_{w,\lambda }$ values obtained through 1-D projections of the mean profile are already good approximations to that obtained through the 3-D projections of the DNS solution at moderate resolutions. However, at higher 
$\varDelta ^+$, there appears to be a minor discrepancy between the two profiles in the form of a constant shift. The 
$C_{w,\lambda }$ for the span-wise velocity component is found to be negative and has a slight 
$\varDelta ^+$ dependence. Similar to the stream-wise velocity, the 
$C_{w,\lambda }$ curves for different 
$Re_\tau$ suggest a 
$\varDelta ^+$ dependence in the span-wise direction as well. The 
$C_w$ in the wall-normal direction is approximately zero, and does not depend on the mesh resolution. This suggests that the wall-normal slip can be set to zero without the loss of any generalisability. This also suggests that the large scales, which are typically resolved in a WMLES simulation, can only slide along the wall but cannot penetrate it. Finally, to obtain a single model form that works for different projection orders, we re-introduce the 
$\lambda$ factor. In the next section, we use the insights gained in this section to improve the performance of the slip wall model by Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) and Bose & Moin (Reference Bose and Moin2014) on the channel flow problem.
5. Towards accurate slip-wall models
While the state-of-the-art dynamic slip-wall model by Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) is found to be better in comparison with the case with no wall model, it is found to be lacking in accuracy when compared with the traditional WMLES approach. In addition, the slip-wall model has been reported to suffer from instability issues when used with certain high-order methods (Carton de Wiart & Murman Reference Carton de Wiart and Murman2017). Thus, there is a need to improve the stability and performance of slip-wall models on canonical turbulent flow problems before it can be confidently used in more complex flows. Indeed, it is recognised that one disadvantage of the traditional approach is that unlike the dynamic slip-wall model, it requires a priori specification of tunable coefficients. The authors are of the opinion, however, that tunable coefficients should not be used as a reason to replace the traditional WMLES approach which has been shown to perform well across a wider range of problems. To this end, we try to use data from existing WMLES simulations and our optimal projection techniques to improve the performance of existing slip-wall models to the level of traditional WMLES for the channel flow problem.
As we observed in § 3, if the grid resolution is sufficiently coarse, a slip velocity is present at the wall. Hence, it is expected that even the solution from the traditional wall-stress-based WMLES will have a slip velocity at the wall. Given the good performance of the traditional WMLES approach for the channel flow problem, it is also expected that the universal relationship given in figure 8 should also hold true for the traditional WMLES approach.
Traditional WMLES solutions were computed using a DG solver with 
$p=3$ discretisation on different meshes using two different sub-grid models: (i) a constant coefficient Smagorinsky model with 
$C_s = 0.12$; (ii) Vreman (Reference Vreman2004) model. Figure 9 shows the comparison of 
$C_{w,\lambda }$ computed using 1-D projection of the Reichardt profile to that computed using the solutions obtained using the traditional WMLES approach. To compute 
$C_{w}$ for a traditional WMLES solution, the solution and its wall-normal gradients are evaluated at the wall to obtain the slip velocity and the pre-multiplied wall-normal gradient. Finally, a least-squares fit is performed to obtain a single value of 
$C_w$. While computing 
$C_w$, the size of the element 
$\varDelta _e$ is required. However, for all the traditional wall-model cases, the size of the element varies in each direction unlike the grids used for projection of DNS. As a first attempt, 
$\varDelta _e$ is taken to be the size of the element in the wall-normal direction. Finally, an optimal value of 
$\lambda$ is found such that the curves collapse asymptotically. By changing the value of 
$\lambda$ only the slope of the asymptotic part of the 
$C_{w,\lambda }$ curve can be changed. However, when the slope of the 
$C_{w,\lambda }$ curve in the asymptotic part was made parallel to the 
$C_{w,\lambda }$ curve obtained for the Reichardt profile by projecting on the 
$p=1$ basis functions, the intercepts were also found to match. This can be seen in figure 9 where the profiles appear identical at large resolutions. However, small discrepancies exist near the lower resolution limit (i.e. the wall-resolved LES limit), suggesting that either the traditional WMLES approach is inaccurate or the sub-grid model is not accurate. Figure 9 also suggests that a universal slip-wall model form exists irrespective of the sub-grid model or the numerical scheme as long as 
$\lambda$ is known.
The 
$\lambda$-normalised 
$C_w$ versus normalised grid-size 
${\varDelta ^+_e}/{p \lambda }$. The 
$\lambda$-normalised 
$C_w$ are computed by 1-D projection of the Reichardt profile and compared with the same obtained using the traditional WMLES solution. The plots marked by ‘FGP’ use only the explicit sub-grid models inside the first element and gradually change to implicit LES outside the first element.
Even if 
$\lambda$ is known prior to the simulation or is dynamically determined, a model for 
$C_{w,\lambda }$ which takes as input the normalised grid size 
${\varDelta ^+_e}/{p \lambda }$ is not useful. This is because the value of 
$\varDelta ^+$ is not known unless the wall stress is also known. One option is to use the traditional wall model to obtain the friction velocity 
$u_{\tau }$ to compute 
$\varDelta ^+$ Whitmore et al. (Reference Whitmore, Griffin, Bose and Moin2021). A better choice would be to represent the slip-wall model coefficient 
$C_{w,\lambda }$ as a function of the mean slip velocity 
$\langle u_s \rangle$ based Reynolds number, i.e. 
$Re_{slip} = {\langle u_s \rangle \varDelta }/{p \lambda \nu }$, as a consequence of which the wall stress will no longer be required to predict 
$C_{w,\lambda }$. Figure 10 shows the 
$\lambda$-normalised slip-wall model coefficient 
$C_{w,\lambda }$ as a function of the slip-velocity (mean stream-wise)-based Reynolds number. As can be observed in figure 10, a universality in the model form similar to the curves in figure 9 also exists in the case when the slip-wall Reynolds number is used as a feature in place of the normalised grid size. In addition, the curves were found to collapse to the 
$p=1$ Reichardt curve for exactly the same value of 
$\lambda$ used in the case of 
$C_{w,\lambda }$ versus 
$\varDelta ^+$. In addition to the plots for 
$C_{w,\lambda }$ for the various traditional approach obtained using various sub-grid models, a model fit is also provided in the figure 9. This fit can be used as a model to specify 
$C_{w}$ at the wall as a function of the slip-wall Reynolds number once 
$\lambda$ is known.
The 
$\lambda$-normalised 
$C_w$ versus the slip velocity 
$u_s$-based Reynolds number 
${\langle u_s \rangle \varDelta }/{p \lambda \nu }$. The 
$\lambda$-normalised 
$C_w$ are computed by 1-D projection of the Reichardt profile and compared with the same obtained using the traditional WMLES solution. The plots marked by ‘FGP’ use the explicit sub-grid models inside the first element only and gradually change to implicit LES outside the first element.
As a first step, we apply the 
$C_w$ computed using the traditional approach and apply it as a slip-boundary condition to check whether the traditional WMLES results can be recreated with the slip boundary condition. At this stage, we are applying the same value of 
$C_w$ for the stream-wise and the span-wise components. In this implementation, it is assumed that the there is no transpiration, i.e. no flow through the wall. To apply the slip-wall boundary condition, we first use the slip velocity components 
$u_{s,i}$ at the wall to compute the wall-normal derivatives of the velocity components 
$u_{h,i}$,
\begin{equation} \frac{\partial u_{h,i}}{\partial n} = \frac{u_{s,i}}{\varDelta C_w}, \end{equation}and finally compute the wall stress at any location using the following formula:
\begin{equation} \tau_{w,i}=\nu \left.\frac{\partial {u}_{h,i}}{\partial n}\right|_{w} -\tau_{i,n}^{S G S}|_{w}. \end{equation}
Hence, contribution of the mean wall stress is only present from the viscous and the sub-grid stresses. In addition to the 
$C_w$ obtained by post-processing the traditional approach solutions, the 
$C_w$ computed using the slip-velocity Reynolds number-based model are also used. The value of 
$\lambda$, required for implementing the slip-velocity Reynolds number-based approach is obtained from the traditional method. The slip-velocity Reynolds number-based model does not require the specification of different 
$C_w$ for each 
$Re_{\tau }$ case; however, it requires one 
$\lambda$ which remains constant across all the cases with different 
$Re_{\tau }$. Figure 11 shows the stream-wise mean velocity profiles, the root mean square of different velocity components and the Reynolds shear stress profiles at different friction Reynolds numbers. The vertical dashed lines show the location of the first, second and third off-wall grid points. For the traditional wall model, the wall-stress is computed using the velocity components at the third off-wall grid point. The slip-wall model does not require any such exchange location. It is clear from figure 11 that when the correct 
$\lambda$ is used, the Reynolds number dependence is captured accurately and the statistics obtained using the slip-wall model are identical to the traditional wall model.
Comparison of the first-order and second-order statistics using the traditional method (Trad.), by re-using the slip-wall model with 
$C_w$ computed from the traditional WMLES solution (Slip) and with the 
$C_w$ computed using the slip-wall Reynolds number formulation (
$Re_{slip}$) at different friction Reynolds number. The vertical dashed lines show the locations of the first, second and third grid points. For the traditional wall model, velocity is sampled at the third off-wall grid point.
The previous tests presented in figure 11 showed that the model is able to capture the Reynolds number dependence on a single grid. The next step is to change the grid resolution and check whether similar results also hold true for the new grid. Before performing numerical experiments with our proposed slip-wall model, an a priori study could be performed by using the results from the traditional WMLES solutions. Two different meshes are now used with 
$N_x \times N_y \times N_z: 12 \times 12 \times 12$ and 
$18 \times 12 \times 18$ elements, respectively. The resolution in the wall-normal direction is kept the same, whereas the resolution in the span-wise and stream-wise case are smaller in the case of the 
$18 \times 12 \times 18$ mesh because the size of the channel is kept constant.
Figure 12 shows the 
$\lambda$-normalised 
$C_w$ obtained for different resolutions for two different types of normalisations. Different normalisations are used because the effective 
$\varDelta$ is not known in the case where the element is not cubic. The plots marked by ‘WN’ and ‘VOL’ use the wall-normal grid spacing and the cube root of the cell volume as 
$\varDelta _e$, respectively. It can be observed that when the wall-normal grid spacing is used as 
$\varDelta _e$, the 
$\lambda$ values required for the two different resolutions are different. This suggests that if the wall-normal grid resolution is used for 
$\varDelta _e$, our proposed slip-wall model will not generalise to a different grid for the same 
$\lambda$ value. On the other hand, when the cell volume was used for 
$\varDelta _e$, the 
$\lambda$ values required to ensure that both the curves collapse was found to be same. This suggests that, for the resolutions considered here, the cube root of the cell volume is an ideal candidate for 
$\varDelta _e$ to ensure that the slip-wall model generalises to a new grid for the same value of 
$\lambda$. Hence, the proposed model will require the specification of the model constant 
$\lambda$ and it is expected to work on different grids and Reynolds numbers. Figure 12 also shows plot of 
$C_{w,\lambda }$ at two other resolutions of 
$12 \times 16 \times 12$ and 
$32 \times 12 \times 16$. While constructing the model fit these resolutions have not been considered. The plots for these specific resolutions are used later to explain the success of the slip-wall model on these unseen resolutions.
The 
$\lambda$-normalised 
$C_w$ versus the mean slip velocity 
$\langle u_s \rangle$ based Reynolds number 
${\langle u_s \rangle \varDelta }/{p \lambda \nu }$. The 
$\lambda$-normalised 
$C_w$ are computed by 1-D projection of the Reichardt profile and compared to the same obtained using the traditional WMLES solution. The plots marked by ‘WN’ and ‘VOL’ use the wall-normal grid spacing and the cube root of the cell volume for specification of 
$\varDelta _e$, respectively.
Figures 13 shows the stream-wise mean velocity profiles, the root mean square of different velocity components and the Reynolds shear stress profiles at different friction Reynolds numbers on two different meshes with 
$12 \times 16 \times 12$ and 
$32 \times 12 \times 16$ elements, respectively, that are not part of the data used for fitting the model for 
$C_{w,\lambda }$. Clearly, the model not only captures the effect of 
$Re_{\tau }$ but also generalises to a new resolution. The performance of the proposed slip-wall model is comparable to the traditional wall model which is a major improvement over the dynamic slip-wall model proposed by Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019). The results indicate that the proposed wall model is able to work reasonably well even at considerably different resolutions. The excellent performance of the slip-wall model can be also explained by computing 
$C_{w,\lambda }$ using the traditional WMLES solutions on these grids. The 
$C_{w,\lambda }$ values estimated using the traditional WMLES solutions from two different meshes with 
$12 \times 16 \times 12$ and 
$32 \times 12 \times 16$ elements, respectively, are plotted in figure 12. The accurate prediction of 
$C_{w,\lambda }$ by the model fit explains the excellent predictive performance of our slip-wall model.
Comparison of the first-order and second-order statistics obtained using the traditional method (Trad.) and the proposed slip-wall model at different friction Reynolds numbers. The solution is computed on two different meshes with 
$12 \times 16 \times 12$ and 
$32 \times 12 \times 16$ elements, respectively, that is not part of the data used for fitting the model. The vertical dash-dotted and dashed lines show the locations of first, second and third grid points for the meshes with 
$12 \times 16 \times 12$ and 
$32 \times 12 \times 16$ elements, respectively. For the traditional wall model, velocity is sampled at the third off-wall grid point.
6. Perspectives on improved slip-wall models
Slip-based wall models (Bose & Moin Reference Bose and Moin2014; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019; Whitmore et al. Reference Whitmore, Griffin, Bose and Moin2021) allow for easy implementation of the wall boundary condition without the need to sample velocity components at a few grid points away from the wall and also allows the possibility to model flow separation. The dynamic slip-wall model proposed by Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) shows good performance on the zero-pressure gradient flat plate case, albeit for a narrow range of 
$Re_\theta$. However, this model requires improvements when compared with traditional wall models for the equilibrium channel flow case at similar resolutions. In § 4, we provided a priori results on the model form for 
$C_w$ for equilibrium channel flows. In § 5, we used some of the insights obtained from § 4 to improve the performance of the existing slip-wall model to at least the traditional WMLES level. Although the performance of our proposed slip wall model was found to be acceptable, more insights from § 4 can be used to further improve the accuracy of both the proposed and the existing slip-wall models. To this end, we outline the following ingredients for the construction of a more generalisable slip-based wall model forms.
(i) The slip model coefficients can be different in the stream-wise, span-wise and the wall-normal directions as observed in figure 8(c). In a more complex 3-D case, the choice of stream-wise, span-wise and the wall-normal direction is a bit ambiguous. However, the mean flow can be used to identify these directions. However, this needs to be iteratively done because the mean flow can itself change when changing these directions. Another approach is to use the flow direction at the first off-wall grid point, similar to how the traditional wall models are implemented. In § 5, we used the same
$C_w$ for all directions. The effect of using different $C_w$ for different velocity components on our proposed model is a topic of further research.(ii) If a dynamic modelling procedure is performed to obtain
$C_w$, the value of $C_w$ cannot be assumed to be same at the original grid and the test filtered grid. Figure 8(c) shows that $C_w$ changes when the resolution is changed from $\varDelta ^+$ to the test filtered grid resolution $2\varDelta ^+$. In addition to $C_w$ being not constant across different grid levels, there is a dependence on the wall units. This dependence is generally not considered in the existing slip-wall model forms. However, this dependence is present in traditional wall models which are found to perform well for equilibrium wall-bounded flow cases. In § 5, we were able to improve the performance of the dynamic slip-wall model of Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) by just augmenting the model form without performing any dynamic procedure.(iii) The discrepancy in figures 9 and 10 between the
$C_{w,\lambda }$ curves obtained by the optimal projection of the Reichardt profile, and that obtained using the solutions of the traditional WMLES approach suggests that the current WMLES approaches are sub-optimal due to the presence of wall modelling and sub-grid modelling errors. This also suggests that there is a lot of scope for improvements, and our optimal projection framework can be used to assess the WMLES performance of different combination of sub-grid models and wall models.(iv) The final comment is on the choice of the parameter that should be used to performing the dynamic procedure. We saw in § 4 that the value of
$\lambda$ effectively captures the effect of the order of projection and hence the numerical method. The corresponding function $g_1$ is fairly universal for different orders. Hence, it is imperative that the dynamic modelling be performed on $\lambda$ rather than $C_w$. The model form for $g_1$ can be empirically obtained from DNS data, the Reichardt profile or from the solution of an existing model such as the wall-stress-based WMLES models. We further observed in § 5 that if the cube root of the cell volume is used for $\varDelta _e$, $\lambda$ remains fairly constant across different resolutions and Reynolds numbers for a given sub-grid model. A dynamic model that determines $\lambda$ without requiring the solutions from the traditional model is a topic of further research.
In this work, two different model forms for 
$C_{w,\lambda }$ were proposed as shown in figures 9 and 10. The first model uses 
$\varDelta ^+$ as a feature, whereas the second model uses the Reynolds number based on the slip velocity as a feature. The implementation of the first model is slightly more complex because the proposed expression for 
$C_w$ is a function of two parameters: 
$\lambda$ and 
$\varDelta ^+$. Assuming that the grid size 
$\varDelta$ is known, to compute 
$\varDelta ^+$ from 
$\varDelta$, an estimate of the average wall stress 
$\langle \tau _w\rangle$ is required:
\begin{equation} \langle\tau_{w}\rangle=v\left\langle\left.\frac{\partial {u}_{h,1}}{\partial y}\right|_{w}\right\rangle-\langle{u}_{h,1} {u}_{h,2}|_{w}\rangle- \langle\tau_{12}^{S G S}|_{w}\rangle. \end{equation}
This average quantity influences the slip velocities through 
$C_w$ which, in turn, affects the average itself. As discussed earlier, an alternate approach is to use the equilibrium wall profile to obtain 
$\tau _w$ as done in Whitmore et al. (Reference Whitmore, Griffin, Bose and Moin2021). However, this requires the sampling of the velocity fields from the off-wall grid points which makes the implementation of slip-wall models as cumbersome as the traditional wall model. An alternate approach is to use the Reynolds number based on the slip velocity as described in § 5. In this work, the optimal estimates of 
$C_w$ were obtained from the DNS solution by projecting on uniform elements of different sizes. However, our projection framework by using anisotropic elements also allows us to study the effect of the grid aspect ratios. As discussed in § 5, one approach to account for mesh anisotropy is to replace 
$\varDelta _e$ with an effective grid size such as the cube root of the cell volume. However, it is advisable to include the aspect ratio in the model form as well to ensure optimal performance across different types of meshes. Finally, the present model form has been derived from the channel flow data and its accuracy in the spatially developing flows such as the flat plate or in the separated flow regions has not been assessed. Given the excellent performance of the traditional wall models on the flat-plate cases we expect similar performance from our proposed model, however, this is a topic of further research.
7. Conclusion
The projection-based scale-separation approach is an essential part of the VMS and uses the grid effectively as a filter. It is applicable to cases where the filter length is anisotropic, varies in space or filtering needs to performed on an unstructured grid. These filter properties were found to be essential for a priori assessment of existing coarse-grained methods for wall-bounded turbulent flows, where the grids can be highly anisotropic and vary in size along a particular direction.
An a priori assessment of the optimal solutions at three different limits, the wall-resolved LES, the HRLES and the WMLES limit, was performed by projecting DNS on different grids suitable for these scale-resolving approaches. For each of these cases, while projecting the DNS on to the coarse space, weak imposition of the boundary condition was made by not enforcing no-slip boundary conditions at the boundary nodes. In the wall-resolved LES limit, the mean velocity was found to be well resolved, no-slip was naturally satisfied and the turbulent stresses were well represented. In the HRLES limit, which was obtained by coarsening the wall-resolved LES mesh in the span-wise and stream-wise directions, the mean velocity was well resolved and the no-slip boundary condition was naturally satisfied. However, the turbulent stresses were found to be well represented only at the centre of the channel and under-represented in the near-wall region where sufficient resolution was not present. In the WMLES limit, which is obtained by further coarsening the HRLES grid in the wall normal direction, the mean profile is no longer represented accurately near the wall and a slip velocity is obtained. The turbulent stresses in WMLES are relatively well represented at the centre of the channel compared to the near-wall region. In the near-wall region, the stream-wise and the span-wise velocity fluctuations were found to be non-zero at the wall, whereas the resolved wall-normal fluctuations and the turbulent shear stress were found to be under-represented. All these trends were found to be consistent with existing solutions in the literature suggesting that the present framework can be utilised to assess, augment and calibrate existing methods.
The ability to obtain the slip-velocity directly from 3-D projections of DNS on coarse near-wall meshes enabled further assessment of the existing slip-wall-based wall models. As a first step, estimates of the slip-wall model coefficient 
$C_w$ were obtained from the mean velocity profile in the inner layer through 1-D projections of the Reichardt profile. The 
$C_w$ estimates from the mean profile were found to be strongly dependent on the order of projection suggesting that the numerical method has considerable effect on the optimal value of 
$C_w$. In addition to this, the resolution for a given slip velocity and projection order was found to scale with the wall units. To make modelling more tractable, we introduced an extra resolution normaliser 
$\lambda$ to express the effect of projection order through a single coefficient, similar in scope to the Smagorinsky model coefficient 
$C_s$. When this analysis was extended to three dimensions, similar dependence on the polynomial order 
$p$ on 
$C_w$ was found for the stream-wise and the span-wise velocity components. However, on re-introduction of resolution normaliser 
$\lambda$ and reusing the 
$\lambda$ values corresponding to the 1-D projections, similar collapse in the 
$C_{w,\lambda }$ values was also observed for the 3-D case. The value of 
$C_{w,\lambda }$ was also found to be different for the stream-wise, span-wise and the wall-normal velocity components.
The ultimate goal of any a priori analysis is to improve the model performance in a posteriori calculations. As a first step towards better slip-based wall models, the performance of existing slip-based wall models was compared with traditional WMLES for channel flows. To establish an equivalence between the two methods, 
$C_{w,\lambda }$ curves were evaluated using the solution of the traditional WMLES approach and compared with the curves obtained for the Reichardt profile. The 
$C_{w,\lambda }$ curves for the traditional WMLES solutions were found to be identical to those obtained using the Reichardt profile at high 
${\varDelta _e^+}/{\lambda p}$. However, at low 
${\varDelta _e^+}/{\lambda p}$, the 
$C_{w,\lambda }$, the curves were found to differ suggesting the presence of sub-grid modelling and wall-modelling errors in the solution. To reduce the implementation challenges associated with using 
${\varDelta _e^+}/{\lambda p}$ as a feature, a slip-velocity Reynolds number-based (
$Re_{slip}$) feature was introduced. Finally, by choosing 
$\varDelta _e$ to be the cube root of the cell volume and re-using 
$\lambda$ from the traditional WMLES solution, a model form was constructed by fitting the 
$C_{w,\lambda }$ versus 
$Re_{slip}$ curve. The resulting model was shown to generalise to different resolutions, element aspect ratios and Reynolds numbers in a posteriori simulations.
Acknowledgements
We acknowledge Professor K. Fidkowski for valuable discussions.
Funding
This research was funded by NASA under the project ‘Scale-resolving turbulence simulations through adaptive high-order discretisations and data-enabled model refinements’, grant number 80NSSC18M0149 (Technical monitor: Dr Gary Coleman).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical computation of $L_2$-projection
The problem of finding an 
$L_2$ projection is equivalent to the problem of finding a 
$u_h\in {\mathcal {V}_h}$ such that
\begin{equation} (u_h,w_h) = (u,w_h) \quad \forall \, {w_h} \in {{\mathcal{V}_h}}.\end{equation}
The first step is to determine the coarse space 
$\mathcal {V}_h$. The coarse space should be low dimensional in comparison with the original solution to ensure that the projection operation acts as a filter. The low-dimensionality of the coarse space can be ensured by using fewer grid points or modes. There are many choices for the coarse space (e.g. the Fourier basis functions, the global Chebyshev polynomial basis functions and the piece-wise polynomial basis functions). Once the coarse space is fixed, the coarse solution can be written as a linear combination of the basis functions as follows:
\begin{equation} u_h = \boldsymbol{w}^{T}_h \boldsymbol{a}_h,\end{equation}
where 
$\boldsymbol {w}^{T}_h$ is a vector of coarse-scale basis functions spanning the coarse space and 
$\boldsymbol {a}_h$ is vector containing the corresponding basis coefficients. Substituting (A2) into (A1), we obtain
\begin{equation} \boldsymbol{M} \boldsymbol{a}_h = \boldsymbol{r},\end{equation}
where the mass matrix 
$\boldsymbol {M}$ and the right-hand side vector 
$\boldsymbol {r}$ is given by
\begin{equation} \boldsymbol{M} = \int \boldsymbol{w}_h \boldsymbol{w}^{T}_h \, {\rm d}\varOmega, \quad \textrm{and} \quad \boldsymbol{r} = (u,\boldsymbol{w}_h). \end{equation}
The coarse-scale basis coefficients are obtained as 
$\boldsymbol {a}_h = \boldsymbol {M}^{-1} \boldsymbol {r}$. The mass matrix 
$\boldsymbol {M}$ is local (block diagonal) when DG basis functions are used. In the case of CG basis functions, the mass matrix has to be assembled by adding contributions from individual element mass matrices. The computation of 
$\boldsymbol {M}$ is not as expensive as compared with the right-hand side vector 
$\boldsymbol {r}$, especially when 
$u$ is high dimensional. The elements of the matrix 
$\boldsymbol {M}$ can be precisely computed using a Gauss quadrature rule appropriate for the order of the polynomial used to define the coarse space. The computation of 
$(u,\boldsymbol {w}_h)$, however, needs special care because it requires the computation of the inner product of a high-dimensional solution 
$u$ with the coarse basis functions 
$w_h$ as shown in figure 14. The high-dimensional solution 
$u$ can come from a finite-difference, finite-volume, spectral or finite-element simulation.
The high-dimensional solution 
$u$ is multiplied with DG coarse-scale basis 
$w_{h,i}$ to obtain 
$uw_{h,i}$. The right-hand side 
$\boldsymbol {r}$ is finally computed by evaluating 
$\int uw_{h,i} \, {\rm d} \varOmega$ for all basis function 
$w_{h,i}$ spanning the coarse space.
A general approach to compute the elements of the right-hand side vector 
$\boldsymbol {r}$ is by using numerical integration. As can be observed in figure 14, the solution obtained after multiplication of the coarse basis functions 
$w_{h,i}$ with 
$u$ still contains high-dimensional features and requires a fine grid for numerical integration. The grid on which 
$u$ exists is assumed to be sufficiently fine for performing the numerical integration. In case the projected solution depends on the order of numerical integration or the size of the numerical integration grid, the solution can be injected on a more finer grid to perform the numerical integration. Once the integration grid is set, the trapezoidal rule or the Simpson formula can be applied to compute the integral over the 
$uw_{h,i}$ fields to obtain the right-hand side vector 
$\boldsymbol {r}$.
Open access
