1. Introduction
Wall-bounded turbulent flows are of particular relevance to many engineering applications. Computational costs of large eddy simulations (LES) increase significantly with the increase in Reynolds number. This is especially true for flows at a friction-velocity-based (
$u_\tau$) Reynolds number 
$Re_\tau \gt 10^3$, which is the range of Reynolds numbers relevant to industrial applications (Smits & Marusic Reference Smits and Marusic2013). A direct numerical simulation (DNS) resolves all the relevant scales of motion and offers the highest possible fidelity (Moin & Mahesh Reference Moin and Mahesh1998). However, substantial grid requirements along with time-step limitations at high Reynolds numbers make DNS infeasible for computing flows of practical relevance. On the other hand, Reynolds-averaged Navier–Stokes (RANS) equations model all the relevant scales of motion and places less restrictive demands on computational costs but offers a lower fidelity (Wilcox Reference Wilcox1998). It may not be a reliable tool for computing flows for which the turbulence models are not calibrated.
A wall-resolved LES (WRLES) resolves dynamically important energy-carrying eddies and models the nearly universal and nearly isotropic small-scale structures, i.e. subgrid scales (SGS) (Sagaut Reference Sagaut2005). For a WRLES of a turbulent boundary layer at high Reynolds number, however, a vast majority of the computational resources have to be spent on the viscous and logarithmic layers since the grid point requirement for each of these layers scale as 
$O(Re_{\tau }^{2})$ (Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016). To alleviate this ‘near-wall problem of LES’, wall-modelled LES (WMLES) offers a practical solution, which aims to bypass the resolution of the inner layer of the turbulent boundary layers. In WMLES, turbulent motions in the inner layer are modelled, whereas outer layer turbulent motions are resolved as in a conventional LES (Piomelli & Balaras Reference Piomelli and Balaras2002). Wall-stress models and hybrid LES/RANS are the two different approaches to model the inner layer and perform WMLES. As these approaches still resolve the outer layer of the turbulent boundary layer, they can -- in principle -- offer better fidelity than RANS techniques.
A hybrid RANS/LES technique, including the detached eddy simulation in a WMLES set-up, uses RANS equations in the inner layer to estimate the wall stress and switches to the LES mode in the outer layer (Heinz Reference Heinz2020). The LES solution is used to feed information to a RANS model at some distance away from the wall. The predictions, however, depend on the choice of the RANS model and the modelling of the RANS/LES interface (Piomelli et al. Reference Piomelli, Balaras, Pasinato, Squires and Spalart2003; Davidson & Dahlström Reference Davidson and Dahlström2005; Davidson & Billson Reference Davidson and Billson2006; Keating & Piomelli Reference Keating and Piomelli2006; Shur et al. Reference Shur, Spalart, Strelets and Travin2008; Choi, Edwards & Baurle Reference Choi, Edwards and Baurle2009). On the other hand, a wall-stress model computes wall shear stress using a log law or the solution of some form of thin boundary layer equations (TBLE) on an embedded grid between the first grid point and the wall (Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016; Bose & Park Reference Bose and Park2018). In this approach, the filtered wall shear stress is estimated at each time step using the LES information from an off-wall grid point on an LES mesh. The wall shear stress is then passed onto the LES grid as a Neumann boundary condition. Several wall-stress modelling strategies with varying complexities have been developed and studied over the years (Deardorff Reference Deardorff1970; Schumann Reference Schumann1975; Piomelli et al. Reference Piomelli, Ferziger, Moin and Kim1989; Piomelli Reference Piomelli1999; Cabot & Moin Reference Cabot and Moin2000; Piomelli & Balaras Reference Piomelli and Balaras2002; Sagaut Reference Sagaut2005; Piomelli Reference Piomelli2008; Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016; Bose & Park Reference Bose and Park2018).
A traditional algebraic equilibrium wall-stress model (EQWM) using a log law, e.g. Reichardt’s profile (Reichardt Reference Reichardt1951), has the advantage of low computational cost. However, it generally performs poorly in non-equilibrium conditions and complex geometries, in particular, involving flows with boundary layer separation (Park Reference Park2017; Goc, Bose & Moin Reference Goc, Bose and Moin2020; Whitmore, Bose & Moin Reference Whitmore, Bose and Moin2024). Moreover, non-monotonic grid convergence in the prediction of the size of turbulent separation bubbles is also observed for the EQWMs (Goc et al. Reference Goc, Bose and Moin2020; Whitmore et al. Reference Whitmore, Griffin, Bose and Moin2021; Agrawal et al. Reference Agrawal, Whitmore, Griffin, Bose and Moin2022). On the other hand, several non-equilibrium wall models have shown promise in predicting separated flows (Balaras, Benocci & Piomelli Reference Balaras, Benocci and Piomelli1996; Wang & Moin Reference Wang and Moin2002; Hickel et al. Reference Hickel, Touber, Bodart and Larsson2013; Park & Moin Reference Park and Moin2014, Reference Park and Moin2016a). However, the two-layer zonal wall models employing embedded grids require some effort in grid generation and domain decomposition in the pre-processing step, especially on unstructured meshes (Bodart & Larsson Reference Bodart and Larsson2011; Park & Moin Reference Park and Moin2016a,Reference Park and Moinb). The model predictions obtained using a wall-stress model depend on the choice of the exchange location (EL), i.e. the location from which the instantaneous LES solution is sampled to feed into the wall model, even for the simple geometry of a turbulent flow in a channel (Kawai & Larsson Reference Kawai and Larsson2012; Frère et al. Reference Frère, Carton de Wiart, Hillewaert, Chatelain and Winckelmans2017). Consequently, the EL becomes a parameter of the simulation, which needs to be adjusted according to the flow field characteristics. Setting the EL requires knowledge of the boundary layer thickness, which is a property of the solution. Moreover, many of the wall-stress models introduce complexities with accompanying empiricism to treat complex flows, e.g. sensors to turn off the wall-stress models at the separation point (Bodart & Larsson Reference Bodart and Larsson2011; Bose & Park Reference Bose and Park2018; Agrawal, Bose & Moin Reference Agrawal, Bose and Moin2024).
Accurate and reliable prediction of separated flows at high Reynolds numbers remains a pacing research issue within the computational fluid dynamics (CFD) community. Several efforts to validate the state-of-the-art WMLES techniques in predicting separated flows at appropriate Reynolds numbers in a realistic external aerodynamics configuration have been undertaken recently. NASA CFD Vision 2030 report (Slotnick et al. Reference Slotnick, Khodadoust, Alonso, Darmofal, Gropp, Lurie and Mavriplis2014)) has identified WMLES for complex three-dimensional flows of practical relevance as one of the key milestones along the CFD technology development roadmap. Park & Moin (Reference Park and Moin2016b), Lehmkuhl et al. (Reference Lehmkuhl, Park, Bose and Moin2018), Goc et al. (Reference Goc, Bose and Moin2020, Reference Goc, Lehmkuhl, Park, Bose and Moin2021) have investigated predictive capabilities of equilibrium and non-equilibrium models within the WMLES framework in the characterization of the flow around an aircraft by considering the JAXA standard model and NASA common research model with wing/body/tail configuration, showing promise in practical applications, yet identifying several areas of improvement.
Since the solution may not be accurately computed near the wall on coarser near-wall LES grids due to the presence of steep wall-normal gradients, one promising alternative to the wall-stress models is the virtual-wall model in which the LES domain is terminated at some finite distance above the wall (Chung & Pullin Reference Chung and Pullin2009; Inoue & Pullin Reference Inoue and Pullin2011; Cheng, Pullin & Samtaney Reference Cheng, Pullin and Samtaney2015). Instantaneous slip velocities obtained using a reduced form of the TBLE are then provided at this virtual boundary or the ‘virtual wall’, which corresponds to the location of the bottom boundary in the LES. The offset of the virtual wall is set to be proportional to the mesh size. This treatment of the wall slip boundary condition has been shown to capture the quantitative features of a separation--reattachment turbulent boundary layer flow at low to moderately large Reynolds numbers. However, identification of the lifted virtual wall can be challenging for complex practical engineering geometric configurations.
Bose & Moin (Reference Bose and Moin2014) propose a wall modelling strategy in which formal boundary conditions for the filtered Navier–Stokes equations are derived instead of relying on the true (unfiltered) boundary conditions for the filtered fields. Unlike traditional wall-stress models and hybrid RANS/LES approaches, the method does not use a wall-stress model or a RANS model in the inner layer to estimate the wall stress. As a result, sampling of the LES solution at the off-wall grid points is not required. The model is derived using the properties of a modified form of the differential filter (Germano Reference Germano1986), and it does not make any assumptions about the local state of the boundary layer or any RANS/LES hybridization. Instantaneous wall slip velocities can be estimated using this slip-wall model when the near-wall solution is under-resolved in the case of coarse grid resolutions. The model is compatible with an arbitrary LES filter and can be motivated using the RANS-type momentum equation (Yang & Bose Reference Yang and Bose2017). It offers a promising alternative to the wall-stress models to predict high-Reynolds-number flows with complex geometries involving separation and reattachment.
The slip-wall model relates the velocity field at the wall to the wall-normal derivative of the velocity field via a wall-adjacent length scale called slip length. The slip length depends on a model coefficient 
$C_w$ and the near-wall grid resolution 
$\unicode{x1D6E5}$. The model recovers the no-slip condition as the near-wall grid is refined and in the limit 
$\unicode{x1D6E5} \rightarrow 0$, and smoothly admits a wall slip velocity as the near-wall grid resolution is coarsened and the flow is no longer fully resolved. The slip-wall model is a general boundary condition applicable to any geometrically complex surface, including two orthogonally or non-orthogonally intersecting walls. Moreover, it is naturally suited to handle boundary layer separation as it will smoothly revert to a no-slip condition at the separation point without additional sensors or damping functions.
The value of the model parameter, i.e. the slip length, is found to depend on the Reynolds number 
$Re_\tau$ of the flow, grid resolution, SGS model and the numerical discretization (Carton & Murman Reference Carton de Wiart and Murman2017; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019; Pradhan & Duraisamy Reference Pradhan and Duraisamy2023). Bose & Moin (Reference Bose and Moin2014) proposed a dynamic procedure to calculate the slip length based on the Germano identity. However, attempts to reproduce the results for a high 
$Re_\tau$ channel flow were unsuccessful (Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019). The wall-stress invariant dynamic wall model (WSIM) of Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) provides an alternate dynamic procedure to estimate the slip length. The model predictions for the channel flow at the high 
$Re_\tau$ cases are found to depend on the grid resolution and grid convergence studies were not carried out. Numerical experiments with arbitrary constant values of slip length using NASA’s discontinuous Galerkin (DG) solver eddy in the implicit LES set-up failed to yield stable computations when applied to a channel flow at 
$Re_\tau \gtrsim 1000$ with a high-order polynomial basis (
$p = 3$ and 
$p = 7$) (Carton & Murman Reference Carton de Wiart and Murman2017).
Existing dynamic versions of the slip-wall model of Bose & Moin (Reference Bose and Moin2014), Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) are found to be sensitive to the numerical implementation details, including the numerical discretization and the choice of SGS model in a finite-volume framework. Moreover, some models (Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019) show a significant log-layer mismatch with the DNS. Given the limitations of the existing dynamic slip-wall models, recent works use a Prandtl mixing length-based model to estimate slip lengths (Whitmore et al. Reference Whitmore, Bose and Moin2024). Other strategies to estimate optimal slip lengths include an optimization procedure to reproduce a known wall shear stress distribution for turbulent channel flows at a range of Reynolds numbers and grid resolutions and model the behaviour of the slip lengths using a curve fit (Whitmore & Bose Reference Whitmore and Bose2023). Ad hoc sensor-based modelling strategies to change the model forms for the slip lengths that switch between the mixing-length-based and parameterized forms in the separation regions have also been studied (Whitmore et al. Reference Whitmore, Bose and Moin2024). Application of the static slip-length models to separated flows in the Boeing speed bump and the JAXA standard model configurations suggests that robust separation predictions require the development of an accurate method for computing slip lengths.
Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) employed an optimal finite-element projection framework to obtain a priori estimates of the wall slip velocity for a typical WMLES using DNS data for a channel flow (Lee & Moser Reference Lee and Moser2015) and proposed improvements to the slip-wall model of Bose & Moin (Reference Bose and Moin2014). The optimal projection framework is used to modify the slip length, and it is represented as a function of the Reynolds number based on local slip-velocity magnitude and near-wall local grid resolution 
$Re_{slip}$. A new model parameter 
$\lambda$ is introduced to represent the effect of the numerical method or the order of projection 
$p$ in the DG set-up and SGS model. Using an a priori estimate for 
$\lambda$, the 
$Re_{slip}$ model for the modified slip length is shown to give good predictions for a range of high 
$Re_\tau$ channel flow cases with the constant coefficient Smagorinsky SGS model using a DG solver with orders of projection up to 
$p = 3$.
The present study begins with the modified form of the slip-wall model proposed by Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023), which uses the 
$Re_{slip}$ model for the modified model coefficient. The main objective is to establish a dynamic modelling procedure for the model parameters. We use the dynamic Smagorinsky model (DSM) as the SGS model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991). The choice of the SGS model is found to be critical to obtain the correct slope of the velocity profile in the log layer (Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019). On the other hand, values of the slip length are found to be responsible for a shift in the mean velocity profile relative to the DNS data, and they do not affect the shape of the mean velocity profile.
The paper is organized as follows. Section 2 presents the DG discretization framework used in the present work. Section 3 provides an overview of the original slip-wall model formulation by Bose & Moin (Reference Bose and Moin2014) along with the modification introduced by Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023). The dynamic modelling procedures of Bose & Moin (Reference Bose and Moin2014) and Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) to calculate the slip length are discussed briefly in § 4. The proposed dynamic modelling strategy to estimate the model parameter 
$\lambda$ using a modified form of the Germano identity and the 
$Re_{slip}$ model is presented in § 5. The key assumptions to arrive at the final form of the dynamic model are also discussed. The proposed dynamic slip-wall model is tested on a range of channel and periodic hill flows in § 6 and results are compared with the available DNS and experimental data along with an EQWM. Finally, conclusions are drawn in § 7.
2. Discontinuous Galerkin discretization
The governing equations in this work are the compressible Navier–Stokes equations in their conservative form written as
\begin{equation} \frac {\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{U}) - \nabla \cdot \mathbf{G}(\mathbf{U}, \nabla \mathbf{U}) = \mathbf{0}, \end{equation}
where 
$\mathbf{U} \in \mathbb{R}^s$ is the conservative state vector of rank 
$s$, consisting of density, momentum and total energy components, 
$\mathbf{F}$ is the inviscid flux and 
$\mathbf{G}$ is the viscous flux. We note that boldface denotes a state vector. We use the DG method for the spatial discretization. The DG method combines the concepts of finite-element and finite-volume methods and allows for high-order approximations, geometric flexibility and natural parallelization. The computational domain 
$\Omega$ is divided into non-overlapping elements 
$K$, each having a subdomain 
$\Omega _K$ and boundary 
$\partial \Omega _K$. These elements can have arbitrary shapes and sizes, allowing for efficient representation of complex geometries. A polynomial approximation is typically used to represent the solution using a 
${L}_2$ projection within each element. The degree of the polynomial 
$p$ can vary, and higher-degree polynomials enable higher-order accuracy. The DG space 
$\mathcal{V}_h$ is defined as
\begin{equation} {\mathcal{V}_h} \triangleq \left \{\phi \in {L_2}(\Omega ):\phi _h\equiv \phi |_{\Omega _K}\in P^p,\forall \Omega _K\in \Omega \right \}, \end{equation}
where the space of polynomials up to degree 
$p$ is denoted as 
$P^p$ and 
$\phi _h$ is the basis function defined on 
$\Omega _K$. Defining 
$\mathcal{V}_h$ in this manner allows for discontinuities in the solution across element boundaries. The elementwise solution 
$\mathbf{U}_h$ that approximates 
$\mathbf{U}$ in 
$\Omega _k$ takes the form
\begin{equation} \mathbf{U}(\mathbf{x},t) \approx \mathbf{U}_h(\mathbf{x},t) = \sum _{j=1}^{n_p} \mathbf{W}_{k,j}(t) \phi _{k,j}(\mathbf{x}), \hspace {0.1in} \mathbf{x} \in \Omega _k, \end{equation}
where 
$\mathbf{U}_{k,j}$ represents the coefficients associated with the 
$j$th basis function 
$\phi _{k,j}$ and 
$n_p$ represents the total number of degrees of freedom within the element 
$k$ of order 
$p$.
The DG method employs a weak formulation of the governing equations that is obtained by multiplying (2.1) by test functions, which are the same as the basis functions, integrating by parts and coupling the elements via numerical fluxes:
\begin{align} & \int _{\Omega _K} \mathbf{\phi }_h^T \frac {\partial \mathbf{U}_h}{\partial t} \, \textrm{d}\Omega \; - \int _{\Omega _K} \nabla \mathbf{\phi }_h^T \cdot \left [ \mathbf{F} \left (\mathbf{U}_h \right ) - \mathbf{G} \left ( \mathbf{U}_h, \nabla \mathbf{U}_h \right ) \right ] \, \textrm{d}\Omega \nonumber \\&
\quad + \int _{\partial \Omega _K} \mathbf{\phi }_h^T \left [\, \widehat {\mathbf{F}} \left ( \mathbf{U}_h^+, \mathbf{U}_h^- \right ) - \widehat {\mathbf{G}} \left ( \mathbf{U}_h^+,\mathbf{U}_h^-,\nabla \mathbf{U}_h^+, \nabla \mathbf{U}_h^- \right ) \right ] \cdot \mathbf{n} \, \textrm{d}S \nonumber \\&\quad - \int _{\partial \Omega _K} \left ( \mathbf{U}_h^+- \{ \mathbf{U}_h \} \right )^T \mathbf{G} \left ( \mathbf{U}_h^+, \nabla \mathbf{\phi }_h^+ \right ) \cdot \mathbf{n} \, \textrm{d}S = \mathbf{0} \quad\text{for all } \mathbf{\phi }_h \in \boldsymbol{\mathcal{V}}_h. \end{align}
Here 
$\partial \Omega _K$ represents the element boundary and, on that boundary, 
$(\cdot )^+$ and 
$(\cdot )^-$ represent quantities taken from the current and neighbouring element, respectively. Approximate numerical fluxes are denoted by 
$\widehat {(\cdot )}$, 
$\{\cdot \}$ represents a face average or boundary value and 
$\mathbf{n}$ is the outward pointing normal vector. The boundary conditions are set through the numerical fluxes.
Substituting (2.3) into (2.4), we get the final update equation that can be written as
\begin{equation} \mathbf{M} \dfrac {{\textrm d} \mathbf{W}}{{\textrm d}t} = -\mathbf{RHS}, \end{equation}
where 
$\mathbf{M}$ is the spatial mass matrix and 
$\mathbf{RHS}$ consists of the volume and surface integrals. Then, the spatial residual vector can be defined as
\begin{equation} \mathbf{R} \equiv \dfrac {{\textrm d} \mathbf{W}}{{\textrm d}t} = -\mathbf{M}^{-1} \mathbf{RHS}. \end{equation}
We solve for the expansion coefficients 
$\mathbf{W}$ that then provide an approximation of the solution to the governing equations over the entire computational domain. The solver used in the present study is discussed in § 6 and Appendix A.
3. Slip-wall modelling
The slip-wall model is essentially a wall boundary condition. The main idea is that the slip velocity is a natural consequence of the near-wall under-resolution of the LES mesh. This has also been shown by Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) using the optimal finite-element projection framework wherein the 
$L_2$ projection of channel flow DNS data onto grids suitable for a WMLES results in slip velocities at the wall. Also, the magnitude of the slip velocity is shown to increase with an increase in near-wall grid under-resolution. This shows that the near equivalence in the boundary conditions for the unfiltered and filtered variables does not hold in the case of a coarse LES when wall modelling becomes necessary. A slip-wall model is an alternative to the traditional wall-stress modelling approach wherein the wall stress is not estimated directly but is indirectly affected through the non-vanishing filtered velocities at the wall. It provides estimates of the slip velocities at the wall when the LES grid resolution is insufficient to accurately resolve the near-wall region and the no-slip condition is not satisfied.
Bose & Moin (Reference Bose and Moin2014) use the properties of a modified differential filter to derive a slip-velocity boundary condition given as
\begin{equation} \overline {u}_i = \overline {C}_w \overline {\unicode{x1D6E5} }_w \dfrac {\partial \overline {u}_i}{\partial n}, \end{equation}
where 
$n$ is the wall-normal direction, 
$\overline {C}_w$ is a tunable model coefficient, whereas 
$\overline {\unicode{x1D6E5} }_w$ is related to the near-wall grid resolution. In (3.1) the slip velocity only depends on the wall-normal derivative of the velocity field and is a direct consequence of the constraint placed on the differential filter that the slip length vanishes at the boundaries. The magnitude of the slip length, i.e. 
$\overline {C}_w \overline {\unicode{x1D6E5} }_w$, imposes a filter length scale at the wall; if it vanishes at the wall then the filtered velocity field will exactly satisfy a no-slip boundary condition. The slip-wall boundary condition smoothly admits a wall slip velocity as the near-wall LES resolution is coarsened and the flow is no longer fully resolved. It is pertinent to note that while (3.1) is derived from a specific choice of the form of the filter kernel, previous studies (Carton & Murman Reference Carton de Wiart and Murman2017; Pradhan & Duraisamy Reference Pradhan and Duraisamy2023) show that the slip-wall model can still perform well even without using the specified filter explicitly.
Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) characterize the wall slip velocity in a WMLES in terms of a Reynolds number based on slip-velocity magnitude and near-wall under-resolution using the optimal finite-element projection framework and propose a modified form of the slip-wall model given by
\begin{equation} \overline {u}_i = \dfrac {\overline {C}_{w,\lambda }}{\lambda } \overline {\unicode{x1D6E5} }_w \dfrac {\partial \overline {u}_i}{\partial n} , \quad \text{where } \overline {C}_w = \dfrac {\overline {C}_{w,\lambda }}{\lambda }. \end{equation}
The model coefficient 
$ \overline {C}_{w,\lambda }$ is a function of the slip-velocity-based Reynolds number 
$\overline {Re}_{slip}$ and 
$ \lambda$, where 
$\overline {Re}_{slip} = \overline {u}_s (\overline {\unicode{x1D6E5} }^{e}_{w} / p)/ \overline {\nu }$. Here, 
$\overline {u}_s$ is the magnitude of the wall slip velocity and 
$\overline {\nu }$ is the kinematic viscosity of the fluid. Note that 
$p$ denotes the order of the polynomial basis used in the DG solver with 
$\overline {\unicode{x1D6E5} }^{e}_{w}$ being the element size adjacent to the wall, and their ratio represents the effective grid size. The model parameter 
$\lambda$ contains the effect of the order of projection 
$p$ and, hence, the numerical method along with the SGS model. Using the above form of 
$\overline {C}_{w,\lambda }$, it is found that given an SGS model 
$ \overline {C}_{w,\lambda } / \lambda$ admits a universal scaling relationship for a particular value of 
$\lambda$ for a wide range of the parameter space. As a result, the model incorporates the effect of Reynolds number, near-wall grid under-resolution, SGS model and numerical discretization.
4. Previous dynamic slip-wall models
Bose & Moin (Reference Bose and Moin2014) presented a dynamic procedure to compute the slip length 
$(\overline {C}_w \overline {\unicode{x1D6E5} }_w)$ in the slip-wall model given by (3.1). It uses a modified form of the Germano identity, which represents the invariance of the total Reynolds stress at the test-filtered level. The model coefficient 
$(\overline {C}_w \overline {\unicode{x1D6E5} }_w)$ is computed as
\begin{equation} \left ( \overline {C}_w \overline {\unicode{x1D6E5} }_w \right )^2 \unicode{x1D6E5} _{R}^{2} \dfrac {\partial \widehat {\overline {u}}_i}{\partial n} \dfrac {\partial \widehat {\overline {u}}_j}{\partial n} + T_{ij} - \widehat {\tau }_{ij} = \widehat {\overline {u}_i \overline {u}_j } - \widehat {\overline {u}}_i \widehat {\overline {u}}_j, \end{equation}
where 
$ \unicode{x1D6E5} _R = ( \widehat {\overline {\unicode{x1D6E5} }}_{w} / \overline {\unicode{x1D6E5} }_w )$ is the ratio of the test filter width to the grid filter width at the wall, and a value of 
$\unicode{x1D6E5} _R = 1.4$ is recommended. Here, 
$\overline {(\cdot )}$ represents a grid-filtered quantity, a hat, i.e. 
$\widehat {(\cdot )}$, denotes the test filtering operation, 
$T_{ij}$ and 
$\tau _{ij}$ depict the SGS stress tensors at the test and grid filter levels, respectively. The slip length is assumed to be equal for the three spatial directions. Equation (4.1) is solved for 
$ ( \overline {C}_w \overline {\unicode{x1D6E5} }_w )$ using a least squares method. The model was tested on a series of high-Reynolds-number channel flows and NACA 
$4412$ airfoil at near-stall conditions.
Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) proposed an alternate dynamic modelling strategy for the slip length 
$ ( \overline {C}_w \overline {\unicode{x1D6E5} }_w )$ as an improvement over the Bose & Moin (Reference Bose and Moin2014) dynamic model. The dynamic model is based on a combination of the invariance of the wall-stress condition under test filtering and a modified form of the Germano identity and is referred to as the wall-stress invariant model (WSIM). The proposed dynamic modelling approach, however, is not unique, and different modelling choices are possible. The dynamic model is given by
\begin{equation} \left ( \overline {C}_w \overline {\unicode{x1D6E5} }_w \right )^2 = \dfrac {L_{ij} M_{ij} + F_{ij} M_{ij}}{M_{kl} M_{kl}}, \end{equation}
where
\begin{align}&\qquad\quad L_{ij} = \overline {u}_i \overline {u}_j - \widehat {\overline {u}}_i \widehat {\overline {u}}_j, \end{align}
\begin{align}& M_{ij} = \left [ \dfrac {\partial \overline {u}_i}{\partial n} \dfrac {\partial \overline {u}_j}{\partial n} - \unicode{x1D6E5} _{R}^{2} \dfrac {\partial \widehat {\overline {u}}_i}{\partial n} \dfrac {\partial \widehat {\overline {u}}_j}{\partial n} \right ], \end{align}
and 
$F_{ij}$ contains different wall stresses, namely Reynolds stress, subgrid stress, viscous stress and pressure tensors computed from the specified velocity field. The model was tested on a statistically stationary plane turbulent channel, a non-equilibrium three-dimensional transient channel and a zero-pressure-gradient flat-plate turbulent boundary layer.
5. A new dynamic slip-wall model
We propose a dynamic procedure to compute the model coefficient 
$\lambda$ in the modified slip-wall model given by (3.2) rather than the slip length in (3.1) following insights from Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023). We start with the Germano identity (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991), which can be written as
\begin{equation} T_{ij} - \widehat {\tau }_{ij} = \widehat {\overline {u}_i \overline {u}}_j - \widehat {\overline {u}}_i \widehat {\overline {u}}_j , \end{equation}
where the SGS stresses at the grid and test-filtered levels are given by
\begin{equation} \tau _{ij} = \overline {u_i u_j} - \overline {u}_i \overline {u}_j \quad \text{and} \quad T_{ij} = \widehat {\overline {u_i u_j}} - \widehat {\overline {u}}_i \widehat {\overline {u}}_j. \end{equation}
Equation (5.1) represents an exact identity and does not involve any assumptions. Subtracting (
$\widehat {\overline {u}_i \overline {u}_j } - \overline {u}_i \overline {u}_j$) from both sides of (5.1), we get
\begin{equation} T_{ij} - \widehat {\tau }_{ij} - ( \widehat {\overline {u}_i \overline {u}_j } - \overline {u}_i \overline {u}_j ) = \overline {u}_i \overline {u}_j - \widehat {\overline {u}}_i \widehat {\overline {u}}_j .\end{equation}
We assume that the slip velocity at the test-filtered level takes a form similar to that for the grid-filtered level, and it is given by
\begin{equation} \widehat {\overline {u}}_i = \widehat {\overline {C}}_w \widehat {\overline {\unicode{x1D6E5} }}_{w} \dfrac {\partial \widehat {\overline {u}}_i}{\partial n}, \end{equation}
where the model coefficient 
$ \widehat {\overline {C}}_w = \widehat {\overline {C}}_{w,\lambda } / \lambda$ has a form similar to that of the coefficient at the grid-filtered level 
$\overline {C}_{w} = \overline {C}_{w,\lambda }/ \lambda$. Next, we assume that 
$\lambda$ is constant between the grid and test-filtered levels. This assumption is based on the findings of Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) using the optimal finite-element projection framework. The universality of the model coefficient 
$\lambda$ depends on the choice of the length scale 
$\unicode{x1D6E5} _w$ used in the slip-wall model. If the cube root of the cell volume is used as the length scale, 
$\lambda$ is found to remain fairly constant across different resolutions and Reynolds numbers for a given SGS model. Next, we assume the slip length to be equal in the streamwise and spanwise directions, while zero in the wall-normal direction, i.e. no transpiration based on the findings of Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023). In general, it can be different in the streamwise, spanwise and wall-normal directions as observed in the a priori studies on the slip-wall model using the optimal finite-element projection framework by Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023). However, numerical experimentation using arbitrary constant values of the slip lengths in the streamwise and spanwise directions for turbulent channel flows did not significantly affect the model predictions. Similar observations are made in the numerical experiments of Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019). Notably, in the previous works on dynamic slip-wall models (Bose & Moin Reference Bose and Moin2014; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019), slip length is assumed to be the same in the three spatial directions.
The model coefficient 
$ \widehat {\overline {C}}_{w,\lambda }$ is assumed to be a function of Reynolds number based on slip-velocity magnitude and the near-wall grid resolution at the test-filtered level along with 
$\lambda$. Substituting for the slip velocities at the grid-filtered level and test-filtered level using (3.1) and (5.4) on the right-hand-side of (5.3), we get
\begin{equation} \left ( \overline {C}_w \overline {\unicode{x1D6E5} }_w \right )^2 \dfrac {\partial \overline {u}_i}{\partial n} \dfrac {\partial \overline {u}_j}{\partial n} - \left ( \widehat {\overline {C}}_w \widehat {\overline {\unicode{x1D6E5} }}_{w} \right )^2 \dfrac {\partial \widehat {\overline {u}}_i}{\partial n} \dfrac {\partial \widehat {\overline {u}}_j}{\partial n} = \overline {u}_i \overline {u}_j - \widehat {\overline {u}}_i \widehat {\overline {u}}_j \hspace {0.02in} . \end{equation}
Now, 
$\overline {\unicode{x1D6E5} }_w$ and 
$\widehat {\overline {\unicode{x1D6E5} }}_{w}$ depend on the grid resolution, 
$p$, and the filter used. On the other hand, the model coefficients 
$ \overline {C}_w$ and 
$ \widehat {\overline {C}}_w$ depend on 
$Re_{slip}$, 
$p$ at the grid- and test-filtered levels, respectively, and the model coefficient 
$\lambda$. In principle, we can use the above equation to find 
$\lambda$ for a given model for 
$ \overline {C}_w$ and, hence, 
$ \widehat {\overline {C}}_w$. However, this would result in a significantly complex nonlinear equation in 
$\lambda$. We choose an alternate approach to simplify the process with an aim to keep a balanced mixture of physical content and mathematical simplicity and rewrite (5.5) as
\begin{equation} \left ( \overline {C}_w \overline {\unicode{x1D6E5} }_w \right )^2 \left [ \dfrac {\partial \overline {u}_i}{\partial n} \dfrac {\partial \overline {u}_j}{\partial n} - C_{wR}^{2} \unicode{x1D6E5} _{R}^{2} \dfrac {\partial \widehat {\overline {u}}_i}{\partial n} \dfrac {\partial \widehat {\overline {u}}_j}{\partial n} \right ] = \overline {u}_i \overline {u}_j - \widehat {\overline {u}}_i \widehat {\overline {u}}_j, \end{equation}
where 
$ C_{wR} = \widehat {\overline {C}}_w / \overline {C}_w$. In this work, we use the value for 
$ \unicode{x1D6E5} _{w}$ as per Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) and 
$\unicode{x1D6E5} _{R}$ is given by
\begin{equation} \unicode{x1D6E5} _{R} = \dfrac {\widehat {\overline {\unicode{x1D6E5} }}_{w}}{\overline {\unicode{x1D6E5} }_w} = \dfrac {p}{p^\star }, \end{equation}
where 
$p^\star$ is the sharp modal cutoff filter order, as discussed in Appendix A. On the other hand, Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) show that the model coefficient 
$C_{w, \lambda }$ is a function of grid resolution, and its value increases when the grid resolution is changed from 
$\unicode{x1D6E5} ^+$ to 
$2 \unicode{x1D6E5} ^+$. In other words, given that the test filter width is coarser than the grid filter width, 
$ \widehat {\overline {C}}_w$ can be expected to be greater than 
$\overline {C}_w$, thereby resulting in the ratio 
$ C_{wR}$ to be greater than one. We use a value of 
$ C_{wR} = 2$ in this work. Sensitivity studies using different plausible values of 
$ C_{wR}$ are shown in Appendix C.
Let
\begin{equation} M_{ij} = \left [ \dfrac {\partial \overline {u}_i}{\partial n} \dfrac {\partial \overline {u}_j}{\partial n} - C_{wR}^{2} \unicode{x1D6E5} _{R}^{2} \dfrac {\partial \widehat {\overline {u}}_i}{\partial n} \dfrac {\partial \widehat {\overline {u}}_j}{\partial n} \right ] \hspace {0.1in} \text{and} \hspace {0.1in} L_{ij} = \overline {u}_i \overline {u}_j - \widehat {\overline {u}}_i \widehat {\overline {u}}_j \end{equation}
for notational convenience. Equation (5.6) can then be equivalently written as
\begin{equation} \left ( \dfrac {\overline {C}_{w,\lambda }}{\lambda } \overline {\unicode{x1D6E5} }_w \right )^2 M_{ij} = L_{ij} \hspace {0.02in} . \end{equation}
Equation (5.9) represents six independent equations in space for a single unknown 
$\lambda$ given the model for 
$\overline {C}_{w,\lambda }$. Thus, the system is overdetermined and we use the method of least squares to obtain 
$\lambda$, which is then given by
\begin{equation} \left ( \dfrac {\overline {C}_{w,\lambda }}{\lambda } \overline {\unicode{x1D6E5} }_w \right )^2 = max \left ( \dfrac {\langle L_{ij} M_{ij}\rangle }{\langle M_{kl} M_{kl}\rangle }, 0 \right ), \end{equation}
where 
$\langle \rangle$ indicates that the numerator and denominator are first averaged over an element followed by an averaging over the directions of homogeneity, i.e. streamwise and spanwise in the case of channel flows and spanwise in the case of periodic flows, and the ratio is clipped to have a maximum value of zero. We can rewrite the above equation as
\begin{equation} \dfrac {\overline {C}_{w,\lambda }}{\lambda } \overline {\unicode{x1D6E5} }_w - l_s = 0, \end{equation}
where
\begin{equation} l_s = \sqrt {max \left ( \dfrac {\langle L_{ij} M_{ij} \rangle }{\langle M_{kl} M_{kl} \rangle }, 0 \right )} \hspace {0.02in} . \end{equation}
Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) provide a model fit for 
$\overline {C}_{w,\lambda }$ based on 
$L_2$ projected channel flow DNS data, which is given by
\begin{equation} \overline {C}_{w,\lambda } = \begin{cases} 0.725 \,\, log_{10}\,\, (\overline {Re}_{slip}/\lambda ) - 0.925 & \text{if} \,\, \log _{10} (\overline {Re}_{slip}/\lambda ) \gt 4.25, \\ 0.6 \,\, log_{10}\,\, (\overline {Re}_{slip}/\lambda ) - 0.41 & \text{if} \,\, \log _{10} (\overline {Re}_{slip}/\lambda ) \gt 3.18 , \\ 0.475 \,\, log_{10}\,\, (\overline {Re}_{slip}/\lambda ) - 0.003 & \text{if} \,\, \log _{10} (\overline {Re}_{slip}/\lambda ) \gt 0.7 , \\ 0.33 & \text{otherwise.} \end{cases} \end{equation}
Substituting for 
$\overline {C}_{w,\lambda }$ in (5.11), we get a nonlinear equation with 
$\lambda$ as the only unknown, which can be found dynamically using a numerical method. The Secant method is used to find 
$\lambda$ using (5.11) and (5.13). The parameter 
$\lambda$, thus obtained, can reach unrealistically high values, especially at high Reynolds numbers on coarse near-wall LES meshes. We prescribe an upper limit to 
$\lambda$ as
\begin{equation} \lambda _f = max(\lambda ,3 \lambda _{CCSM}), \quad \text{where } \lambda _{CCSM} = 1.125. \end{equation}
Here, 
$\lambda _{CCSM}$ is the value of 
$\lambda$ for the constant coefficient Smagorinsky model obtained in the a priori analysis of channel flow DNS data and traditional wall model predictions using the optimal finite-element projection framework in Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023). The limiter value is based on numerical experimentation across the range of Reynolds numbers and flow configurations considered in this work. Finally, we use 
$\lambda _f$ in (3.2).
The dynamic model in (5.11) is essentially identical to the dynamic slip-wall model of Bose & Moin (Reference Bose and Moin2014). The model given by (5.11) can be transformed to (4.1) using the Germano identity. However, the model form in (4.1) is found to be sensitive to the implementation details including the numerical discretization and the choice of the SGS model in a finite-volume framework in previous studies (Park & Moin Reference Park and Moin2016b; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019). Our attempts to use this model form in the DG framework with 
$p \geqslant 2$ resulted in unstable simulations with the constant coefficient and the dynamic Smagorinsky SGS models.
We remark that the dynamic modelling procedure to obtain 
$\lambda$ as discussed above is not unique. The model coefficient can be obtained using a number of modelling choices, e.g. modified form of the Germano identity used in Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019). Equation (5.10) has a form similar to that of the dynamic model of Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019), but it does not contain the additional wall-stress terms in the numerator, which originates from the invariance of the wall-stress condition under test filtering. Those additional wall-stress terms are expected to predict the same wall stress regardless of the grid resolution (or filter) and act as an effective self-regulating mechanism to control the changes in slip length to predict the correct wall stress. In our case, a similar effect is obtained by enforcing the Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) model for 
$\overline {C}_{w,\lambda }$, which is a function of the slip-wall velocity and near-wall grid resolution via the slip-based Reynolds number. The parameter 
$\overline {C}_{w,\lambda }$ and, hence, the slip length appropriately changes when the near-wall grid is refined or coarsened and/or when the Reynolds number is increased or decreased, resulting in an appropriate change in the wall stress. Note that the wall stress can be related to the slip length 
$ (\overline {C}_{w} \overline {\unicode{x1D6E5} }_w)$ and slip velocity 
$U_{slip}$ using the slip-wall model as 
$\tau _w = (\nu + \nu _{SGS}) ( U_{slip} /(\overline {C}_{w} \overline {\unicode{x1D6E5} }_w) )$, where 
$\nu _{SGS}$ is the SGS or eddy viscosity.
5.1. Implementation of the wall boundary condition
We assume that there is no transpiration and slip is only allowed in the wall-parallel directions. It is useful to note that, the slip-wall model allows for transpiration as considered in the previous studies of Bose & Moin (Reference Bose and Moin2014), Carton & Murman (Reference Carton de Wiart and Murman2017) and Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019). However, using the optimal finite-element projection framework, Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) show that the slip length associated with the wall-normal velocity is approximately zero. Hence, it can be set to zero without significant loss of generalizability. In the current implementation of the dynamic slip-wall model, we compute the wall-normal derivatives of the slip-velocity components at the wall using (3.2). The wall-parallel slip-velocity components at the wall are computed using the solution inside the element adjacent to the wall. This is then used to compute the wall stress and it is applied as a Neumann boundary condition complemented by a slip boundary condition for the velocity. The numerical implementation is done using the following steps.
(i) At every integration point, a ghost value is created, where the wall-parallel slip-velocity components at the wall are obtained from the element interior state
$\mathbf{U}^{+}_{h}$ as
(5.15)Wall-normal velocity gradients are also calculated using the interior solutions.
\begin{equation} \overline {u}^{b}_{h,i} = \overline {u}^{+}_{h,i} - \overline {u}^{+}_{h,j} n_j n_i. \end{equation}
-
(ii) Slip-wall parameters
$C_{w,\lambda }$ and $\lambda$ are computed using the dynamic slip-wall model. -
(iii) The wall-normal derivatives of the slip-velocity components at the wall are then computed using the slip-wall model given by (3.2) as
(5.16)
\begin{equation} \left .\frac {\partial \overline {u}_{h,i}}{\partial n}\right |_{w} = \frac {\overline {u}^{b}_{h,i}}{\unicode{x1D6E5} (C_{w,\lambda }/\lambda )}. \end{equation}
-
(iv) Finally, wall-stress components at the quadrature points of the boundary faces are computed using
(5.17)where we consider the contribution of the mean wall stress from the viscous and subgrid stresses for the wall stress.
\begin{equation} \tau '_{w}= \left ( \nu + \nu _{SGS} \right ) \left .\frac {\partial \overline {u}_{h,i}}{\partial n}\right |_{w}, \end{equation}
-
(v) The projected wall stress
$\tau _{w,i}= \tau '_{w} n_i$ is applied as a Neumann boundary condition.
6. Numerical experiments
In this work we use CaslabDG, an in-house DG solver for the computations. The governing equations are the filtered compressible Navier–Stokes equations in conservative form. The solver was successfully used previously to compute statistically stationary channel flows at high 
$Re_\tau$ (Pradhan & Duraisamy Reference Pradhan and Duraisamy2023) using a constant coefficient Smagorinsky model with up to 
$3$ orders of the polynomial basis 
$p$. The solver is parallelized using the message passing interface (MPI). Inviscid fluxes are approximated using the Roe approximate Riemann solver (Roe Reference Roe1981). An SGS model is used for the unresolved SGS stresses in the filtered Navier–Stokes equations. The SGS viscosity is added to the molecular viscosity and the viscous flux contains both molecular and turbulence contributions. The second form of Bassi and Rebay (Bassi & Rebay Reference Bassi and Rebay2000), popularly known as the BR-2 scheme, is used for the viscous fluxes. The governing equations are marched in time using an explicit third-order Runge--Kutta total variation diminishing (RK3-TVD) scheme.
The solver uses the Lagrange nodal basis evaluated at the Gauss–Legendre quadrature points and the number of quadrature points 
$ngp$ in each of the three directions is related to the polynomial degree of approximation by 
$ngp = (p + 2)$. The integrals are approximated using the Gauss quadrature rule. The basis and test functions are created using a tensor product of the one-dimensional Lagrange interpolating polynomials that forms a non-hierarchical nodal basis. The corresponding number of degrees of freedom in each element is 
$(p + 1)^3$. We use a Lagrange polynomial basis of degree 
$p = 2$ for all the WMLES computations. It is to be noted that we do not use an explicit filter for the spatial filtering operation, but rely on implicit filtering through the numerical discretization and grid resolution. Also, the polynomial basis degree of 
$p = 2$ used in this work does not warrant for polynomial dealiasing, which can be achieved by explicitly filtering the solution at every time step (Diosady & Murman Reference Diosady and Murman2013; Gassner & Beck Reference Gassner and Beck2013; Brazell et al. Reference Brazell, Brazell, Stoellinger and Mavriplis2015).
The results obtained using the dynamic slip-wall model are compared with those obtained using an EQWM. For the EQWM, we compute the wall friction 
$\tau _w$ from the instantaneous velocity taken at the furthest distance from the wall inside the first element. The computed wall friction is then used as the Neumann boundary condition applied at the quadrature points of the boundary faces. Our implementation is similar to the work of Carton & Murman (Reference Carton de Wiart and Murman2017). We use the Reichardt function of the form
\begin{equation} u^+ = \dfrac {1}{\kappa } ln(1 + \kappa z^+) + 6.646 \left [ 1 - exp\left ( \dfrac {-z^+}{11} \right ) - \dfrac {z^+}{11} exp\left ( \dfrac {-z^+}{3} \right ) \right ], \end{equation}
where 
$ \kappa = 0.38$, as the EQWM, and use the Newton–Raphson method to iterate on the values of 
$u^+$ and 
$y^+$. This wall function supports the theoretical velocity profile down to the wall. We have chosen this approach for its simplicity and efficiency, and it is shown to give excellent results for statistically steady channel flows at high Reynolds numbers in Carton & Murman (Reference Carton de Wiart and Murman2017). Note that, for the simulations of separated flows in periodic hill configuration, we do not use any adjustments to the implementation of the EQWM like turning off the model in the separation region.
6.1. Sharp modal cutoff filter as a test filter
The dynamic modelling procedure requires filtering at two different levels, i.e. grid filter and test filter, to calculate the value of the model coefficient 
$\lambda$. In a DG framework this is equivalent to using two different orders of polynomial basis for approximating the solution. The Lagrange interpolation polynomials, which are used as the basis functions within our work, are not hierarchical, i.e. every basis function contains high-order solution content. As a result, unlike a spectral method, we cannot directly use a sharp cutoff filter to remove the higher-order modes. To reduce the order of projection that would result in a coarser filtering operation, the solution coefficients need to be transformed to a modal representation, the hierarchical form of which allows for a classification of solution modes based on polynomial degree. The solution can then be coarsely filtered by setting the higher-order modes to zero or by scaling the higher-order coefficients by a factor 
$ \alpha \in [0, 1]$. In this work, we use a cutoff filter order 
$p^\star = 1$ and set the modes of degree greater than one to zero for the test filter. This is equivalent to assuming a test filter to be about twice the width of the grid filter, which is generally followed in finite difference or finite-volume methods (Pope Reference Pope2000). Once the filtered forms of modal solution coefficients are obtained, an inverse transformation is performed to get the filtered nodal solution coefficients, thereby obtaining the coarse-filtered solution. In this work, we follow the procedure outlined by Brazell et al. (Reference Brazell, Brazell, Stoellinger and Mavriplis2015) to implement the sharp modal cutoff filter in our solver, and it is discussed in detail in Appendix A.
6.2. Dynamic Smagorinsky model
The DSM is a simple eddy viscosity model that relates the unresolved SGS stresses to the resolved strain rate 
$\overline {S}_{ij}$ via a turbulent viscosity 
$\nu _{SGS}$ as
\begin{align} \tau _{ij}^{SGS} = 2 \overline {\rho } \nu _{SGS} \overline {S}_{ij}, \quad \text{where } \overline {S}_{ij} = \dfrac {1}{2} \left ( \dfrac {\partial \overline {u}_i }{\partial x_j } + \dfrac {\partial \overline {u}_j }{\partial x_i } \right ) . \end{align}
The SGS eddy viscosity is related to a characteristic velocity and a length scale on dimensional grounds, and it is given by
\begin{equation} \nu _{SGS} = ( C_S \unicode{x1D6E5} )^2 | \overline {S} | , \end{equation}
where 
$C_S$ is the Smagorinsky coefficient, 
$ | \overline {S} | = \sqrt {2 \overline {S}_{ij} \overline {S}_{ij} }$ is the strain-rate magnitude and 
$\unicode{x1D6E5}$ is the filter width or a representative grid size. The DSM improves upon the original Smagorinsky model by dynamically adjusting the model coefficient 
$(C_S \unicode{x1D6E5} )$ based on local flow properties. The idea is to seek a more accurate representation of the turbulence, especially in regions with varying flow conditions. The DSM also provides a near-wall correction that can lead to proper near-wall behaviour of the SGS viscosity without the use of wall-damping functions. The dynamic calculation of the coefficient is based on an explicitly performed second-level filter operation called the test filter that is applied to the grid-filtered variables. As mentioned before, we denote the test filter operation by a hat and we use the sharp modal cutoff filter as the test filter, as discussed in (6.1). The model coefficient 
$(C_S \unicode{x1D6E5} )$ is calculated as
\begin{equation} (C_s \unicode{x1D6E5} )^2 = \dfrac {1}{2} \dfrac { L_{ij}^{d} M_{ij} }{M_{kl} M_{kl}}, \end{equation}
where the Leonard stress tensor 
$L_{ij}$ and its deviatoric part 
$ L_{ij}^{d}$ are given by
\begin{align} L_{ij} = \widehat { \overline {u}_i \overline {u}_j} - \widehat { \overline {u}}_i \widehat { \overline {u}}_j , \quad L_{ij}^{d} = L_{ij} - \dfrac {1}{3} L_{kk} \delta _{ij} \end{align}
and
\begin{equation} M_{ij} = \widehat { | \overline {S} | \overline {S}_{ij} } - \unicode{x1D6E5} ^{2}_{R} | \widehat {\overline {S}} | \widehat {\overline {S}}_{ij} \hspace {0.02in} . \end{equation}
The derivative and the test filter operations do not commute for the sharp modal cutoff filter. We follow Brazell et al. (Reference Brazell, Brazell, Stoellinger and Mavriplis2015) to determine the second term of 
$M_{ij}$ by computing the test-filtered velocity followed by the derivatives of the test-filtered velocity to form the strain-rate tensor. Another possible choice is to use the test filter operation on the grid-filtered strain rate, but this approach does not have any advantages over the method used here as shown by Brazell et al. (Reference Brazell, Brazell, Stoellinger and Mavriplis2015). The parameter 
$\unicode{x1D6E5} _R$ is calculated as per the recommendation of Brazell et al. (Reference Brazell, Brazell, Stoellinger and Mavriplis2015) and it is given by
\begin{equation} \unicode{x1D6E5} _{R} = \dfrac {\widehat {\overline {\unicode{x1D6E5} }}_{w}}{\overline {\unicode{x1D6E5} }_w} = \dfrac {p + 1}{p^\star + 1}. \end{equation}
The numerator in (6.4) can assume local negative values and this is physically consistent as it corresponds to energy backscatter, i.e. energy from the SGS scales is transferred back to the resolved scales. However, negative SGS viscosity values can numerically destabilize the simulation, especially when the sum 
$(\nu + \nu _{SGS})$ becomes negative. Thus, it is customary to perform some type of averaging of the numerator and denominator, generally in the directions of homogeneity. In this work, we perform two-step averaging. First, the numerator and denominator are averaged over an element to get their representative single values in each element. After this, the numerator and denominator are averaged over the homogeneous directions to get the final averaged numerator and denominator as 
$ \langle L_{ij}^{d} M_{ij} \rangle$ and 
$\langle M_{kl} M_{kl}\rangle$, respectively. Finally, the ratio 
$[0.5 (\langle L_{ij}^{d} M_{ij}\rangle / \langle M_{kl} M_{kl} \rangle )]$ is clipped to get non-negative values.
6.3. Application to statistically stationary channel flows
The new dynamic slip-wall model is applied to a series of statistically stationary turbulent channel flows that are homogeneous in directions parallel to the wall. The fully developed turbulent flow between the two parallel walls is separated by a distance 
$2 \delta$ in the 
$z$ direction, where 
$\delta$ is the half-channel height. The flow is assumed to be periodic in the streamwise (
$x$) and spanwise (
$y$) directions. The simplicity of geometry and boundary conditions makes this canonical flow configuration an appealing test case, and it has been used to validate the performance of previous dynamic slip-wall models (Bose & Moin Reference Bose and Moin2014; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019).
The friction Reynolds number is set as 
$Re_\tau = \overline {\rho } u_\tau \delta / \mu$, with 
$u_\tau = \sqrt {\tau _w/ \overline {\rho }}$, 
$u_\tau$ being the friction velocity based on the wall shear stress 
$\tau _w$ taken from the available DNS. The friction Reynolds number is imposed through a constant forcing in the 
$x$-momentum equation using the pressure gradient 
$ \langle {\textrm d}p/{\textrm d}x \rangle = - \tau _w / \delta$. A number of WMLES studies use the constant pressure gradient as the forcing strategy in finite-volume (Yang & Bose Reference Yang and Bose2017; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019) and DG frameworks (Carton & Murman Reference Carton de Wiart and Murman2017; Frère et al. Reference Frère, Carton de Wiart, Hillewaert, Chatelain and Winckelmans2017; Lv et al. Reference Lv, Yang, Park and Ihme2021). For a DNS, Quadrio, Frohnapfel & Hasegawa (Reference Quadrio, Frohnapfel and Hasegawa2016) show that the specific choice of the forcing term, i.e. constant mass flow or constant pressure gradient, does not produce important statistical consequences and one-point statistics do not show any appreciable difference between the simulations.
The size of the computational domain is 
$2 \pi \delta$ in the 
$x$ direction and 
$\pi \delta$ in the 
$y$ direction. The degree of polynomial 
$p$ used for all the simulations presented here is 
$2$ and the sharp modal cutoff order 
$p^\star$ is 
$1$. For all the cases considered, the flow is initially evolved for at least time 
$20 \delta / u_\tau$ units and statistics are sampled for an additional 
$10 \delta / u_\tau$ time units. One-point statistics including the mean velocity and Reynolds shear and normal stresses are compared with the DNS of Lee & Moser (Reference Lee and Moser2015) and Hoyas, Oberlack & Laux (Reference Hoyas, Oberlack, Alcántara-Ávila, Kraheberger and Laux2022).
The performance of the proposed dynamic slip-wall model is validated using the cases listed in table 1, which shows the simulated Reynolds numbers and the grid resolutions in inner and outer layer units. The meshes are uniform in the streamwise, spanwise and wall-normal directions. The first element size in the wall-normal direction for all the considered cases is significantly coarser than a conventional LES mesh and the resolution is insufficient to resolve near-wall turbulent structures. As a result, none of the simulated Reynolds numbers with the grids given in table 1 are wall resolved.
Summary of mesh parameters for the different simulated Reynolds numbers. Here, 
$\unicode{x1D6E5} _x$, 
$\unicode{x1D6E5} _y$ and 
$\unicode{x1D6E5} _z$ are the effective grid sizes in the streamwise (
$x$), spanwise (
$y$) and wall-normal (
$z$) directions, respectively, 
$\delta$ is the half-channel height and 
$\unicode{x1D6E5} _x^+$, 
$\unicode{x1D6E5} _y^+$ and 
$\unicode{x1D6E5} _z^+$ are normalized with wall units. Here
$N_x$, 
$N_y$ and 
$N_z$ represent the number of elements in the streamwise, spanwise and wall-normal directions, respectively. The number of degrees of freedom in each direction is given by 
$(p+1) N_x$, 
$(p+1) N_y$ and 
$(p+1) N_z$, where 
$p$ is the degree of the polynomial basis. Note that the numerical experiments are labelled following the convention: [dynamic slip-wall model (DSW)]-[
$Re_\tau$]-[grid resolution].
A grid sensitivity study for channel flow at 
$Re_\tau \approx 10\,000$ is shown in figure 1. Starting with a coarse mesh with 
$8 \times 8 \times 8$ elements, the number of elements in each of the three directions is doubled at each level of refinement. This corresponds to three different coarse near-wall resolutions of 
$\unicode{x1D6E5} z = 0.125 \delta , 0.0625 \delta ,$ and 
$0.03125 \delta$, i.e. 
$\unicode{x1D6E5} _{z_w}^{+} \approx 1253$, 
$626$ and 
$313$. The details of the mesh parameters are given in table 1. For all the cases, the first off-wall grid point lies in the log layer. The results are plotted starting from the second off-wall element at 
$(p+1)$ quadrature points in each element. It can be seen that the mean velocity for the coarsest mesh G1 has a slight positive log layer mismatch that reduces upon grid refinement. The difference between the model predictions at each of the successive grid refinement levels is less than 
$1 \,\%$. A grid refinement study for the other two Reynolds number cases shows a similar trend.
Grid refinement study for the proposed dynamic wall model showing comparisons between model predictions and DNS for the streamwise mean velocity at 
$Re_\tau \approx 10\,000$. (a) Classical visualization and (b) visualization focused on the bulk profile and the top interface of the first element.
A comparison between the dynamic slip-wall model and the EQWM predictions with the DNS is shown in figure 2. The model predictions are obtained on grid G2 with 
$16 \times 16 \times 16$ elements. The slope of the mean velocity profiles obtained using the two models is similar, however, there is a slight shift between them. Both the mean velocity profile predictions match well with the DNS. The two model predictions for the Reynolds shear stress profiles also match the DNS well. The streamwise Reynolds stress predicted by the two models is also similar but there is a slight mismatch with the DNS. On the other hand, the spanwise and wall-normal Reynolds stress profiles obtained using the two models closely agree with the DNS.
Comparison between the proposed dynamic slip-wall model and EQWM predictions using grid G2 with the DNS for the (a) full mean velocity profile, (b) mean velocity in the log region, (c) Reynolds shear stress, (d) root-mean-square (r.m.s.) velocity fluctuations at 
$Re_\tau \approx 10\,000$.
One-point statistics on grid G2 at 
$Re_\tau \approx 2000, 5200$ and 
$10\,000$ are presented in figure 3 and compared with the available DNS. The first and second moments agree well with the DNS at the three Reynolds numbers. The model parameter 
$C_w$ and mean streamwise slip velocity 
$U_{slip}$ at the three Reynolds numbers and grid resolutions are plotted in figure 4. It is important to note that the slip-wall model is sensitive to the numerical implementation details, including the numerical discretization and the choice of the SGS model. Consequently, different optimal slip lengths are required for the correct prediction of the wall stress depending on the numerical set-up. As a result, the quantitative assessment of model parameters like 
$C_w$ or 
$\lambda$ with DNS is difficult. However, 
$C_w$ variation with the near-wall grid resolution 
$\unicode{x1D6E5} _{w}^{+}$ follows an expected trend that is qualitatively similar to that observed by Whitmore & Bose (Reference Whitmore and Bose2023) for the optimal slip length estimates, i.e. larger slip lengths on coarser near-wall grid resolutions and decay of the slip length in the limit 
$\unicode{x1D6E5} _{w}^{+} \to 0$. Moreover, the slip-wall velocity trend is also consistent with the a priori filtering tests of Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) using the optimal finite-element projection framework, i.e. an increase in the slip velocity with an increase in the Reynolds number on identical grids, and an increase in slip velocity at the same Reynolds number upon coarsening the grid.
Comparison between DNS and proposed dynamic slip-wall model predictions using grid G2 for the (a) full mean velocity, (b) mean velocity profile in the log region, (c) Reynolds shear stress, (d) r.m.s. velocity fluctuations at 
$Re_\tau \approx 2000, 5200$ and 
$10\,000$.
Slip parameter 
$C_w$ (filled symbols) and streamwise slip velocity 
$U_{slip}$ (unfilled symbols) as a function of near-wall grid resolution 
$\unicode{x1D6E5} _{w}^{+}$ at 
$Re_\tau \approx 2000$ (squares), 
$5200$ (circles) and 
$10\,000$ (deltas). Colour code: black for grid G1, red for grid G2, blue for grid G3.
Change in the mass flow rate relative to the nominal value plotted as a function of the flow through time obtained at 
$Re_\tau \approx 2000$ on grid G2.
The model parameter 
$C_w$ and, hence, the slip length are found to increase as the near-wall grid is coarsened at the considered Reynolds numbers (see figure 4). The slip length modifies how momentum is transferred between the wall and the fluid. It effectively tunes the wall stress to account for near-wall turbulence effects without fully resolving the boundary layer. A larger slip length increases the velocity slip at the wall as seen in figure 4, leading to a higher velocity near the boundaries. For the considered cases, this results in the reduction of wall shear stress. For the same driving pressure gradient, this implies that more of the pressure gradient contributes to accelerating the bulk flow instead of being dissipated near the walls. As a result, the velocity profile is modified such that the mass flow rate 
$U_b = (1/2 \delta ) \int _{-\delta }^{\delta } \overline {u(z)} {\textrm d}z$ and, hence, the bulk Reynolds number 
$Re_b = U_b \delta / \nu$ are found to increase for these cases as evident from the shift in mean velocity profiles relative to the DNS in figures 1 and 3. The maximum increase in the mass flow rate relative to the nominal value is about 
$6\,\%$ on the coarsest grid G1 and the increase reduces upon grid refinement. Figure 5 shows the change in the mass flow rate variation relative to its nominal value with the flow through time (
$L_x/U_b$) at 
$Re_\tau \approx 2000$ on grid G2, and similar trends are observed for the other cases.
It is important to note that, the underlying assumptions for the slip condition become invalid for significantly larger near-wall grid resolutions commensurate to larger filter sizes. The shift in the mean velocity relative to the DNS is expected to increase as the near-wall grid is coarsened, as seen in figure 1. Similar observations are made in the study of Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) using the constant pressure gradient forcing term. Note that the increase in mass flow upon coarsening the grid is also observed in LES with the no-slip boundary condition albeit on relatively finer grids than the WMLES grids.
A quantitative assessment of the dynamic slip-wall model is performed in terms of the normalized 
$L_2$ error in the streamwise mean velocity 
$U^+$ predictions with respect to the DNS for all the cases presented in table 1. The calculations exclude the first near-wall element. The normalized 
$L_2$ error is determined between the second off-wall element 
$\unicode{x1D6E5} ^{+}_{2e}$ and the half-channel height 
$\delta ^+$ as
\begin{equation} \mathcal{E} = \left [ \dfrac {\int _{\unicode{x1D6E5} ^{+}_{2e}}^{\delta ^{+}} (U^{+}_{DSW} - U^{+}_{DNS})^2 {\textrm d}z^+}{\int _{\unicode{x1D6E5} ^{+}_{2e}}^{\delta ^{+}} (U^{+}_{DNS})^2 {\textrm d}z^+} \right ]^{1/2} . \end{equation}
Here, 
$U^{+}_{DSW}$ and 
$U^{+}_{DNS}$ represent the mean velocity obtained using the proposed dynamic slip-wall model and DNS, respectively. The error is evaluated at 
$(p+1)$ quadrature points within each element in the wall-normal direction 
$z^+$, and the integration for each element is performed using quadrature. The error 
$\mathcal{E}$ is plotted as a function of the representative grid size 
$\unicode{x1D6E5}$ in figure 6. We consider 
$\unicode{x1D6E5}$ based on element volume, 
$\unicode{x1D6E5} = (\unicode{x1D6E5} _x \unicode{x1D6E5} _y \unicode{x1D6E5} _z)^{1/3}$ with 
$\unicode{x1D6E5} _x$, 
$\unicode{x1D6E5} _y$ and 
$\unicode{x1D6E5} _z$ taken as the effective grid sizes in the streamwise, spanwise and wall-normal directions, respectively. The 
$L_2$ error slightly increases with an increase in Reynolds number on an identical grid. However, the maximum error is less than 
$3 \,\%$ for all the cases considered here, demonstrating the performance of the model at practically relevant Reynolds numbers on significantly under-resolved near-wall LES mesh resolutions.
Normalized 
$L_2$ error, 
$\mathcal{E}$, in streamwise mean velocity 
$U^+$ as a function of grid resolution 
$\unicode{x1D6E5}$ at 
$Re_\tau \approx 2000, 5200$ and 
$10\,000$.
Snapshots of the (a) normalized streamwise slip velocity and (b) vorticity magnitude on the bottom wall obtained using the new dynamic slip-wall model at 
$Re_\tau \approx 10\,000$ using grid G2.
Instantaneous snapshots of the streamwise slip velocity normalized by the friction velocity 
$u_\tau$ on the bottom wall at 
$Re_\tau \approx 10\,000$ employing the G2 grid is shown in figure 7(a). The mean slip velocity at the wall increases as the Reynolds number increases and the simulation resolves a smaller fraction of the inner layer of the boundary layer. The mean streamwise slip velocities at the wall are approximately 
$10.9 u_\tau$, 
$13.9 u_\tau$ and 
$15.35 u_\tau$ for 
$Re_\tau \approx 2000, 5200$ and 
$10\,000$, respectively, on grid G2; the centreline velocity is approximately 
$28 u_\tau$. This behaviour is consistent with the a priori filtering tests of Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) using the optimal finite-element projection framework. A snapshot of the vorticity magnitude levels on the bottom wall is also shown. The visualization of the near-wall eddies is shown in figure 8 using the Q criterion.
Isosurafces of the Q criterion coloured with normalized streamwise velocity 
$\overline {u}/u_\tau$ obtained using the new dynamic slip-wall model at 
$Re_\tau \approx 10\,000$ using grid G2.
6.4. Application to separated flows
We next apply the new dynamic slip-wall model to periodic hill flows at different Reynolds numbers. The flow configuration consists of a channel flow with constrictions and forms a generic case of an internal flow separating from a curved surface. Periodic boundary conditions are applied in the streamwise (
$x$) and spanwise (
$y$) directions. The flow separates at the hill crest resulting in a large recirculation bubble, and it reattaches further downstream. A Reynolds number based on the bulk velocity at the crest 
$Re_b = \overline {\rho } U_b h/\mu$ determines the flow conditions for this case with 
$U_b$ being the bulk velocity and 
$h$ the hill height. The constant mass flow rate is ensured by adding a source term in the 
$x$-momentum equation. This forcing term is dynamically adjusted to provide the correct mass flow rate at the hill crest, and therefore, the correct bulk Reynolds number. The mass flux is constant to five digits of accuracy for the present computations after an initial transient.
The periodic hill case has been extensively studied over the past 15 years, both experimentally and numerically. Rapp & Manhart (Reference Rapp and Manhart2011) performed experiments in a water channel at 
$Re_b$ ranging from 
$5600$ to 
$37\,000$. Several DNS and LES (Breuer et al. Reference Breuer, Peller, Rapp and Manhart2009; Balakumar, Park & Pierce Reference Balakumar, Park and Pierce2014; Diosady & Murman Reference Diosady and Murman2014; Gloerfelt & Cinnella Reference Gloerfelt and Cinnella2015; Krank, Kronbichler & Wall Reference Krank, Kronbichler and Wall2018) studies have also been conducted. Many studies have also been performed to test the performance of WMLES (Balakumar et al. Reference Balakumar, Park and Pierce2014; Carton & Murman Reference Carton de Wiart and Murman2017). The availability of high-quality data from experiments, DNS and LES makes this a good benchmark test case to evaluate the performance of the slip-wall model in the presence of separation and reattachment processes. Carton et al. (Reference Carton de Wiart, Larsson and Murman2018) place two-dimensional periodic hill cases at level 4 complexity in their list of benchmark test cases, which grow in complexity from level 1 to level 5. This test case is well documented and well posed with consistent DNS/LES predictions between different codes that match well with the experiments. As pointed out in Gloerfelt & Cinnella (Reference Gloerfelt and Cinnella2019), this benchmark test case has been a choice in several European projects and workshops to investigate the reliability of RANS/LES strategies.
The periodic hill flows involve massive separation on the hill’s leeward sides, the length of which is about 
$50 \,\%$ of that of the periodic segment. The principal challenge of this flow arises from the separation on the curved hill surface and the fact that the reattachment point, and hence, the whole flow, are highly sensitive to the separation process. The flow exhibits complex dynamics, including separation, reattachment, an unsteady shear layer, a large recirculation bubble and strong acceleration on the windward wall. Resolving these delicate flow details, especially on significantly coarser grids is a challenging task. For example, in a recent work, Zhou, He & Yang (Reference Zhou, He and Yang2021) developed a data-driven wall model using the feedforward neural network and training data from the WRLES of the periodic hill flows with different Reynolds numbers and hill geometries. The model shows a good performance for the turbulent channel flows at 
$Re_\tau$ as high as 
$5200$. However, significant discrepancies in the mean velocity and Reynolds stress predictions are observed between their data-driven wall model and the WRLES for the periodic hill case with validation studies limited to 
$Re_b \approx 10\,600$.
The size of the computational domain is 
$L_x = 9h$, 
$L_y = 4.5h$ and 
$L_z = 3.035h$ in the streamwise (
$x$), spanwise (
$y$) and wall-normal (
$z$) directions, respectively. Piecewise third-order polynomial functions give the coordinates of the curved hill, and the second hill geometry is described by the same equations with a horizontal translation (Rapp & Manhart Reference Rapp and Manhart2011). We use two grids; a coarse grid with 
$50 \times 24 \times 9$ elements, i.e. 
$150 \times 72 \times 27$ (
$ = 0.2916$ million) degrees of freedom, and a fine grid with 
$75 \times 36 \times 15$ elements, i.e. 
$225 \times 108 \times 45$ (
$ \approx 1.1$ million) degrees of freedom. In comparison to our grids, a DNS of 
$Re_b \approx 10\,600$ performed by Krank et al. (Reference Krank, Kronbichler and Wall2018) using a seventth-order DG solver, used 
$128 \times 64 \times 64$ elements, i.e. 
$896 \times 448 \times 448$ (
$ \approx 180$ million) degrees of freedom, whereas to perform an implicit LES, a mesh with 
$448 \times 224 \times 224$ (
$ \approx 22.5$ million) degrees of freedom was used.
The grids are approximately uniform in the streamwise and spanwise directions and a mild stretching is used in the wall-normal direction. The mesh is perpendicular to the wall in the first cell away from the wall. The effective element sizes at the hill crest, a key region for the periodic hill flow, are 
$\unicode{x1D6E5} _{x} (= \unicode{x1D6E5} ^{e}_{x}/p) \approx 0.105h$ and 
$\unicode{x1D6E5} _{z} (= \unicode{x1D6E5} ^{e}_{z}/p) \approx 0.093h$ for the coarse grid and 
$\unicode{x1D6E5} _{x} \approx 0.065h$ and 
$\unicode{x1D6E5} _{z} \approx 0.064h$ for the fine grid. Figure 9 shows the two grids used in the computations. We consider two high-Reynolds-numbers cases of 
$ Re_b \approx 10\,600$ and 
$37\,000$ for which high-quality experimental data are available.
Coarse and fine grids used to compute the periodic hill flows.
Effect of grid refinement on the mean velocity profiles in the streamwise (U) and vertical (W) directions at different stations for the 
$Re_b \approx 10\,600$ case. Red solid lines, grid G1; blue solid lines, grid G2; unfilled circles, Rapp & Manhart (Reference Rapp and Manhart2011) experiment; filled circles, WRLES by Breuer et al. (Reference Breuer, Peller, Rapp and Manhart2009).
We first study the effect of mesh resolution on the dynamic slip-wall model predictions at 
$Re_b \approx 10\,600$. Wall-normal variation of the mean streamwise and wall-normal velocity profiles obtained on the coarse and fine grids is shown in figure 10, while Reynolds streamwise, wall-normal and shear stress profiles are shown in figure 11. The model predictions are compared with the experimental data and the WRLES results of Breuer et al. (Reference Breuer, Peller, Rapp and Manhart2009) at four streamwise locations of 
$x = 1h, 2h, 4h$ and 
$8h$ that cover the separated as well as post-reattachment regions. The mean streamwise and wall-normal velocity profile predictions at these locations on the two grids closely match each other and they compare well with the experiments and WRLES. Note that the wall-normal velocity prediction at 
$x/h = 8$ is particularly sensitive to the grid refinement. This has also been observed in previous LES studies of Gloerfelt & Cinnella (Reference Gloerfelt and Cinnella2019). Notably, the WRLES results for wall-normal velocity, which has a lower amplitude than the mean streamwise velocity, show discrepancies with the experiments at 
$x/h = 8$. Interestingly, WRLES results match well with the dynamic slip-wall predictions at this location.
The reverse flow velocities are captured well on the two grids. The Reynolds shear stress profiles on the two grids are also similar to each other and they show a good match with the experiment at the four locations. On the other hand, streamwise and wall-normal Reynolds stress profile predictions on the coarse grid follow the qualitative trend well and the predictions improve on the fine grid and get closer to the experimental data. It is to be noted that, despite the very coarse grid resolution, the agreement between the resolved part of the Reynolds shear and normal stresses with the measurement is reasonably good.
The periodic hill flow is then computed at a higher Reynolds number of 
$Re_b = 37\,000$ using the fine grid and comparisons between the new dynamic slip-wall model and the EQWM predictions for the mean velocities and Reynolds stresses are shown in figures 12 and 13, respectively. The EQWM significantly underpredicts the separation and shows a faster recovery and the mean streamwise and wall-normal velocity profiles in the separated regions show a significant mismatch with the experiments. On the other hand, the dynamic slip-wall model accurately captures the separation and shows an excellent match with the experiments for the mean velocities at the four locations in figure 12. A discrepancy is observed for the wall-normal velocity at 
$x/h = 8$ similar to the 
$Re_b \approx 10\,600$ case, which is consistent with the study of Gloerfelt & Cinnella (Reference Gloerfelt and Cinnella2019). The Reynolds shear stress profiles predicted by the dynamic slip-wall model also closely agree with the experiments whereas the EQWM model predictions show a considerable mismatch. This is also the case for the streamwise Reynolds stress profiles. However, the dynamic slip-wall model predictions for the wall-normal Reynolds stress show some discrepancies and overpredict the levels found in the experiments. This is again consistent with the LES studies of Gloerfelt & Cinnella (Reference Gloerfelt and Cinnella2019) that showed dramatic overprediction for the wall-normal Reynolds stresses while getting a good match for the other Reynolds stress components. Overall, the dynamic slip-wall model predictions are considerably better than the EQWM.
Effect of grid refinement on the profiles of Reynolds stresses in the streamwise (U) and vertical (W) directions at different stations for the 
$Re_b \approx 10\,600$ case. Red solid lines, EQWM; red solid lines, grid G1; blue solid lines, grid G2; unfilled circles, Rapp & Manhart (Reference Rapp and Manhart2011); experiment; filled circles, WRLES by Breuer et al. (Reference Breuer, Peller, Rapp and Manhart2009).
Mean velocity profiles in the streamwise (U) and vertical (W) directions at different stations for the 
$Re_b \approx 37\,000$ case. Red solid lines, EQWM; blue solid lines, dynamic slip-wall model; unfilled circles, Rapp & Manhart (Reference Rapp and Manhart2011) experiment.
Profiles of Reynolds stresses in the streamwise (U) and vertical (W) directions at different stations for the 
$Re_b \approx 37\,000$ case. Red solid lines, EQWM; blue solid lines, dynamic slip-wall model; unfilled circles, Rapp & Manhart (Reference Rapp and Manhart2011) experiment.
Plot of 
$\lambda$-normalized slip parameter 
$C_w$, i.e. 
$C_{w,\lambda }$, as a function of 
$Re_{slip}$ and 
$\lambda$ at 
$Re_b \approx 10\,600$ (squares) and 
$Re_b \approx 37\,000$ (circles). CColour code: black for grid G1; red for grid G2.
Figure 14 plots the variation of 
$\lambda$-normalized slip parameter 
$C_w$, i.e. 
$C_{w,\lambda }$, as a function of the slip-velocity-based Reynolds number 
$Re_{slip}$ and 
$\lambda$ for the two Reynolds numbers of 
$Re_b \approx 10\,600$ and 
$37\,000$ obtained on grids G1 and G2. The curves for 
$C_{w,\lambda }$ are found to collapse as a function of 
$Re_{slip}$ for the two grid resolutions, and the universal scaling relationship discovered in Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) for turbulent channel flows using an a priori analysis is found to remain valid for these separated flows as well in the posteriori computations. The parameter 
$C_{w,\lambda }$ follows the expected trend; assumes lower values at the lower Reynolds number and coarser near-wall grid resolutions and decays as the near-wall grid is coarsened or the Reynolds number is decreased. The streamwise variation of the slip parameter 
$C_w$ is shown in figure 15 and it also follows the expected trend at the two Reynolds numbers and grid resolutions. Variations in 
$C_w$ along the streamwise direction are also observed for these cases. The sudden peaks at the start of the hill and at 
$x/h \approx 7$ for the high-Reynolds-number case of 
$Re_b \approx 37\,000$ are indicative of the unphysical effects due to the coarse grid resolution and these peaks vanish on the fine grid G2.
Streamwise distribution of the slip parameter 
$C_{w}$ at 
$Re_b \approx 10\,600$ (solid lines) and 
$Re_b \approx 37\,000$ (dashed lines) obtained on grid G1 (black lines) and grid G2 (red lines).
Streamwise distribution of the slip velocity 
$U_{slip}$ at 
$Re_b \approx 10\,600$ (black line) and 
$Re_b \approx 37\,000$ (red line) obtained on grid G2.
The streamwise distributions of the mean streamwise slip-wall velocity at 
$Re_b \approx 10\,600$ and 
$37\,000$ obtained on the fine grid G2 are shown in figure 16. The slip velocity 
$U_{slip}$ for the two cases changes the sign near the separation point, remains negative in the separation bubble and reverses the sign near the reattachment location. For the 
$Re_b \approx 10\,600$ case, the locations where the slip velocity changes the sign closely match the experimental measurements and WRLES results for the separation and reattachment points at 
$x/h \approx 0.2$ and 
$x/h \approx 4.6$, respectively. As indicated by the negative values of the slip velocity, the separation region slightly reduces in size at the higher Reynolds case, which is physically consistent.
7. Conclusion
Several strategies have been proposed to bypass the stringent near-wall grid resolution requirement for performing LES of high-Reynolds-number flows in the presence of solid walls. In this work, we focus on the slip-wall modelling approach -- originally proposed by Bose & Moin (Reference Bose and Moin2014) -- and replace the conventional no-slip velocity boundary condition with slip velocities at the wall. The major objective is to accurately capture the mean flow characteristics at Reynolds numbers of practical relevance using a significantly coarse near-wall LES mesh, and do so in a robust manner.
We present a new formulation of a dynamic slip-wall model that is consistent with the DG framework and is tightly integrated with DG operators. The model coefficients of the modified slip-wall model of Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) are based on a priori estimates obtained using an optimal finite-element projection framework. Here, we propose a dynamic modelling procedure to compute the scaling parameter 
$\lambda$ for the slip-wall model coefficient 
$C_w$. The dynamic part of the model is based on a modified form of the Germano identity and coupled with the DSM. The level of under-resolution is represented by a slip Reynolds number and the proposed model attempts to also incorporate the effects of the numerical discretization and the SGS model.
The canonical case of a statistically stationary turbulent channel flow is first used to validate the new dynamic slip-wall model. The model predictions are compared with the available DNS data at three Reynolds numbers of 
$Re_\tau \approx 2000, 5200$ and 
$10\,000$. Grid independence studies are performed at these Reynolds numbers by considering significantly under-resolved LES meshes with streamwise, spanwise and wall-normal grid resolutions corresponding to 
$\unicode{x1D6E5} _{x} \simeq 0.1-0.4\delta$, 
$\unicode{x1D6E5} _{y} \simeq 0.05-0.2\delta$ and 
$\unicode{x1D6E5} _{z} \simeq 0.03-0.125\delta$, respectively. These mesh resolutions are significantly coarser than the WMLES mesh recommendations of Larsson et al. (Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016) corresponding to 
$\unicode{x1D6E5} _{x} \simeq 0.08 \delta$, 
$\unicode{x1D6E5} _{y} \simeq 0.05 \delta$ and 
$\unicode{x1D6E5} _{z} \simeq 0.01-0.05 \delta$. Mean velocity profiles show an excellent match with the DNS at the considered Reynolds numbers on all the grids with an 
$L_2$ error less than 
$3 \,\%$ for all the cases. Reynolds shear and normal stress profiles resolved on the significantly coarse grids also show excellent agreement with the DNS. The model performance is shown to be similar to that of the EQWM, which is known to predict the equilibrium wall-bounded flows without separation accurately. This is a considerable improvement over the dynamic slip-wall model of Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) that shows a significant log layer mismatch at similar Reynolds numbers but on comparatively finer grid resolutions.
The model performance is evaluated in flow separation and reattachment over periodic hills at Reynolds numbers of 
$Re_b \approx 10\,600$ and 
$37\,000$ using two different grid resolutions. The meshes used for the computations are significantly coarser than the conventional LES meshes, e.g. the fine mesh used here has about 
$20$ times fewer degrees of freedom than the implicit LES performed by Krank et al. (Reference Krank, Kronbichler and Wall2018). The streamwise and wall-normal mean velocity profile predictions obtained using the dynamic slip-wall model on the two grids compare well with the experimental data in the separated and post-reattachment flow regions at 
$Re_b \approx 10\,600$. Reynolds shear stress predictions obtained using the two grids also match very well with experiments at different streamwise locations. However, the Reynolds normal stresses are better predicted on the fine grid. Computations at 
$Re_b \approx 37\,000$ using the fine grid show that the dynamic slip-wall model predictions for the mean velocity profiles agree well with the experiments. The Reynolds shear stress profiles are also in excellent agreement with the experiments, with some discrepancies in the Reynolds normal stress predictions. On the other hand, EQWM for this case shows significant discrepancies with the experimental data for the mean velocities as well as Reynolds shear and normal stresses.
The new model contains parameters like 
$C_{w,\lambda }$, 
$ C_{wR}$ and 
$\unicode{x1D6E5} _R$ that are empirically established, but the model does not assume the state of a boundary layer. It is important to note that the main purpose of a slip-wall model is similar to that of a traditional wall-stress model, i.e. accurate estimation of wall shear stress. Achieving this goal without prior assumptions regarding the state of the boundary layer or embedded empirical parameters is an outstanding challenge. Moreover, the instantaneous velocity field is intertwined with the effects of the LES grid resolution and Reynolds number for a given numerical discretization and SGS model as discussed in Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019). The modelling choices made in this work are consistent with the observations of Pradhan & Duraisamy (Reference Pradhan and Duraisamy2023) and the works of Bose & Moin (Reference Bose and Moin2014), Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019). The empirical parameters in the proposed model, especially 
$C_{w,\lambda }$, provide an explicit reference to how the near-wall flow should behave at different near-wall grid resolutions and Reynolds numbers. The model is found to be somewhat insensitive to the parameters 
$ C_{wR}$ and 
$\unicode{x1D6E5} _R$ in their plausible range for the considered cases, which includes the smooth body-separated flows with separation and reattachment.
The new model can consistently predict mean velocity and Reynolds shear and normal stress profiles for the equilibrium as well as separated flows at high Reynolds numbers using significantly coarse near-wall LES meshes. The model performs at a computational cost similar to the EQWM, which is the cheapest state-of-the-art WMLES strategy. In the authors’ opinion, the excellent performance of the model may be attributed to the integration of the optimal finite-element projection framework used to obtain the slip-wall parameters with the consistent dynamic procedures for the SGS and slip-wall modelling coupled with the DG framework. This work is a step towards making the slip-wall model a viable computing tool for predicting complex engineering flows, and further evaluations are required.
While the present work demonstrates the implementation and validation of the dynamic slip-wall model within a DG framework, we emphasize that the fundamental approach is not inherently tied to DG methods. The essential requirement is a rigorous coarse graining or scale separation operator, which can be provided by various numerical frameworks including finite-element methods, variational multiscale approaches or projection-based methods. For finite-volume and finite difference techniques, while the implementation is less straightforward, agglomeration-based techniques (see, e.g. Gravemeier Reference Gravemeier2006) can provide viable pathways. Changing the numerical framework (for instance, to continuous finite elements) would necessitate recalibration of the base parameters, along with a redefinition of the length scales. This requirement stems from our fundamental observation that subgrid and wall models cannot be decoupled from the underlying numerical method due to the strong interactions between unresolved and coarsely resolved dynamics.
Acknowledgement
This research was funded by NASA under the project `Scale-resolving turbulence simulations through adaptive high-order discretizations and data-enabled model refinements’, grant number 80NSSC18M0149 (technical monitor: Dr Gary Coleman). We acknowledge Professor Krzysztof Fidkowski, Dr Gary Coleman and Dr Aniruddhe Pradhan for their valuable discussions.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Sharp modal cutoff filter implementation
The implementation of the sharp modal cutoff filter is a three-step process: transforming the nodal solution coefficients to a hierarchical modal representation, applying a filter on the modal coefficients and then transforming back into the nodal representation. Let us denote the nodal solution coefficients as 
$u_j$, the nodal basis functions as 
$\phi _j$, the modal solution coefficients as 
$b_j$ and the modal basis functions as 
$\psi _j$. Then, we can write the approximation 
$u_h$ to any flow variable 
$u$ in an element as
\begin{equation} u \approx u_h = \sum _{j=1}^{p+1} u_j \phi _j = \sum _{j=1}^{p+1} b_j \psi _j \hspace {0.02in} . \end{equation}
Multiplying the above equation by 
$ \psi _i$ and integrating over the standard element, we get
\begin{equation} \sum _{j=1}^{p+1} C_{ij} u_j = \sum _{j=1}^{p+1} M_{ij} b_j \end{equation}
or in matrix and vector form, we can write
\begin{equation} [C] \boldsymbol{u} = [M] \boldsymbol{b} \hspace {0.02in} . \end{equation}
Here, the modal mass matrix [M] is given by
\begin{equation} [M] = M_{ij} = \int _{\Omega _k} \psi _i \psi _j \,{\textrm d}x . \end{equation}
The mixed mass matrix [C] is given by
\begin{equation} [C] = C_{ij} = \int _{\Omega _k} \psi _i \phi _j \,{\textrm d}x . \end{equation}
Using (A3), we can obtain the modal solution coefficients, 
$\boldsymbol{b}$, from the nodal solution coefficients 
$\boldsymbol{u}$ by inverting the modal mass matrix [M], i.e.
\begin{equation} \boldsymbol{b} = [M]^{-1} [C] \boldsymbol{u}. \end{equation}
Now that the hierarchical modal basis coefficients have been obtained, a square filter matrix, [F], can be applied as a matrix-vector product, i.e.
\begin{equation} \widehat {\boldsymbol{b}} = [F] \boldsymbol{b}, \end{equation}
where 
$\widehat {\boldsymbol{b}}$ are the filtered modal solution coefficients. The sharp cutoff filter matrix 
$[F]$ is diagonal with its entries being 
$0$ or 
$1$. If all entries are 
$1$ giving the identity matrix, the filtering operation returns the original solution. To obtain a cutoff filter of order 
$(p^\star + 1)$, all diagonal entries of the filter matrix are 
$1$ up to and including the (
$p^\star + 1$) diagonal entry with the rest of the entries 
$0$. The last step in the modal decomposition filtering procedure is to transform the filtered modal coefficients back to nodal basis coefficients to give the filtered nodal solution. This reverse transformation can be performed as follows :
\begin{align} \widehat {\boldsymbol{u}} & = [C]^{-1} [M] \widehat {\boldsymbol{b}}, \\[-10pt] \nonumber \end{align}
Let
\begin{equation} [B] = [M]^{-1} [C]\quad \text{and} \quad \widehat {[F]} = [B]^{-1} [F] [B] \hspace {0.02in} . \end{equation}
We can then write the final filtered nodal solution coefficients as
\begin{align} \widehat {\boldsymbol{u}} = \widehat {[F]} \boldsymbol{u}. \end{align}
The final filter matrix 
$ \widehat {[F]}$ can be assembled as a pre-processing step as it does not have a dependence on the solution.
Comparisons of analytical solution for a function 
$u(x) = cos(2x) + 0.3sin(8x)$ in the range 
$ [-\pi ,\pi ]$ with the best fit obtained using one-dimensional nodal DG employing 10 elements and 
$p = 3$ along with filtered solutions for 
$p^\star = 0,1,2$ and 
$3$. Unfilled green circles indicate the quadrature points within each element.
The sharp modal filter is first tested in a one-dimensional DG set-up. The objective is to find a best fit for the function 
$ u(x) = cos(2x) + 0.3sin(8x)$. The degree of the Lagrange basis function is set to 
$p = 3$, and the sharp modal cutoff filter is tested for 
$p^{\star } = 0,1,2$ and 
$3$. The number of elements used is 
$10$ with four quadrature points on each element, and the results are presented in figure 17. The analytical solution is plotted for the domain 
$[-\pi ,\pi ]$. The nodal DG solution with 
$p = 3$ matches closely with the analytical solution, and the jumps at the element approximation denote the discontinuous nature of the approximation. The sharp modal cutoff filter with 
$p^{\star } = 3$ does not affect the solution and the results are identical to the original solution. On the other hand, lower cutoff orders of 
$p^{\star } = 2,1$ and 
$0$ result in a piecewise quadratic, linear and constant solution, respectively.
The modal sharp cutoff filter in its one-dimensional form discussed above is extended to three dimensions in a tensor product fashion and applied before every RK3-TVD step in our in-house DG code. The effect of the test filtering operation using the modal sharp cutoff filter on the normalized instantaneous streamwise velocity 
$\overline {u}/ u_\delta$ is shown in figure 18. The degree of polynomial used is 
$p = 3$ and results are shown for the 
$Re_\tau \approx 544$ case with filter orders of 
$p^\star = 3,2,1$ and 
$0$. The snapshots of 
$\overline {u}/ u_\tau$ show the loss of information and decrease in resolution of the flow field as the filter cutoff order is reduced.
Snapshots of normalized streamwise velocity in the 
$xz$ plane passing midway through the spanwise domain dimension with 
$p^\star = 3, 2, 1$ and 
$0$ showing the effects of the test filtering operation at 
$Re_\tau \approx 544$.
Appendix B. Wall-resolved LES at $Re_\tau \approx 544$
A WRLES of the turbulent channel flow is performed at 
$Re_\tau \approx 544$ using the DSM as the SGS model to verify our in-house DG solver. The mesh size is 
$36 \times 30 \times 24$ elements in the streamwise, spanwise and wall-normal directions, respectively. The grid is uniform in the streamwise and spanwise directions and it is geometrically stretched in the wall-normal direction with a stretching ratio of 
$1.2$. The effective grid sizes in each direction in wall units are 
$\unicode{x1D6E5} _x^+ \approx 38$, 
$\unicode{x1D6E5} _y^+ \approx 19$ and 
$\unicode{x1D6E5} _{z}^+$ at the wall is 
$ \unicode{x1D6E5} _{z_w}^+ \approx 4.5$ and at the channel centre is 
$ \unicode{x1D6E5} _{z_c}^+ \approx 32.2$. The grid resolution is based on the recommendation of Bose & Moin (Reference Bose and Moin2014) for a WRLES, i.e. 
$\unicode{x1D6E5} _x^+ \lesssim 50$, 
$\unicode{x1D6E5} _y^+ \lesssim 30$ and 
$\unicode{x1D6E5} _{z_w}^+ \sim O(1)$. In comparison, the grid resolution of the available DNS (Lee & Moser Reference Lee and Moser2015) is 
$\unicode{x1D6E5} _x^+ \approx 8.9$, 
$\unicode{x1D6E5} _y^+ \approx 5$, 
$ \unicode{x1D6E5} _{z_w}^+ \approx 0.019$ and 
$ \unicode{x1D6E5} _{z_c}^+ \approx 4.5$. Please note that, the effective grid sizes 
$\unicode{x1D6E5} _x$, 
$\unicode{x1D6E5} _y$ and 
$\unicode{x1D6E5} _z$ for the finite-element grid are defined as 
$\unicode{x1D6E5} _x = \unicode{x1D6E5} ^e_x/p$, 
$\unicode{x1D6E5} ^e_y/p$ and 
$\unicode{x1D6E5} ^e_z/p$, respectively. The quantities 
$\unicode{x1D6E5} ^e_x$, 
$\unicode{x1D6E5} ^e_y$ and 
$\unicode{x1D6E5} ^e_z$ represent the actual element sizes in the finite-element mesh.
(a) Snapshot of normalized streamwise velocity in the 
$xz$ plane passing midway through the spanwise dimension, (b) snapshot of the vorticity magnitude on the bottom wall and (c) isosurface of the Q criterion coloured with normalized streamwise velocity 
$\overline {u}/u_\tau$ for a WRLES at 
$Re_\tau \approx 544$ obtained using the DSM SGS model.
Wall-normal variation of the (a) mean Smagorinsky coefficient 
$\langle C_s \rangle$, (b) mean velocity, (c) Reynolds shear stress and (d) r.m.s. velocity fluctuations for a WRLES employing DSM as the SGS models compared with the DNS at 
$Re_\tau \approx 544$.
The instantaneous streamwise velocities 
$\overline {u}$ normalized with 
$u_\tau$ for the WRLES in the 
$xz$ plane are shown in figure 19(a). The solution is reasonably resolved. Snapshots of the vorticity magnitude on the bottom wall are shown in figure 19(b). On the other hand, figure 19(c) shows the isometric view of the isosurfaces of the Q criterion to visualize the near-wall eddies.
The wall-normal variation of the mean Smagorinsky coefficient 
$\langle C_s \rangle$ along with mean velocity and Reynolds shear and normal stresses for the WRLES compared with the DNS is shown in figure 20. The Smagorinsky coefficient 
$C_s$ for the DSM assumes a value of zero at the wall as 
$L_{ij}^{d}$ in (6.4) is equivalently zero at the wall owing to the no-slip velocity boundary condition. It gradually increases in the viscous sublayer before reaching a value of about 
$0.12$ in the log layer at about 
$ z^+ \approx 60$ after which it remains close to 
$0.1-0.12$ till the half-channel height 
$\delta ^+$. The mean velocity and Reynolds shear and normal stress profiles obtained using the WRLES are nearly identical to the DNS data.
Appendix C. Sensitivity analysis
The proposed dynamic slip-wall model involves two model parameters, namely 
$ C_{wR}$ and 
$\unicode{x1D6E5} _R$. The definition of the filter operation fixes the value of 
$\unicode{x1D6E5} _R$, which is the ratio of the test filter width to the grid filter width. We used a value for 
$\unicode{x1D6E5} _R$ as recommended by Brazell et al. (Reference Brazell, Brazell, Stoellinger and Mavriplis2015). Numerical experiments using different values for 
$\unicode{x1D6E5} _R$ in the plausible range of 
$ \unicode{x1D6E5} _{R} = [1,2]$ for 
$p = 2$ and 
$p^\star = 1$ resulted in negligible differences in the results and these observations are similar to those made by Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2019) for their dynamic slip-wall model. On the other hand, the parameter 
$ C_{wR}$ comes into the picture because of the use of different values of the model coefficient 
$C_{w}$ at the test filter and grid filter levels. The sensitivity to 
$ C_{wR}$ is tested for values in the plausible range 
$ C_{wR} = [1,2]$. Results for the two extreme values in this range, i.e. 
$ C_{wR} = 1$ and 
$2$, are shown in figure 21 for the 
$Re_\tau \approx 10\,000$ case obtained using grid G2. The effect of 
$ C_{wR}$ on the mean velocity and Reynolds stress predictions is also found to be negligible and the results are almost identical. This suggests that the model coefficient 
$ C_{w}$ can be taken to be the same at the test- and grid-filtered levels, which is the general practice (Bose & Moin Reference Bose and Moin2014; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019).
Effect of the parameter 
$ C_{wR}$ values on the proposed dynamic wall model predictions for the case DSW-10000-G2 along with DNS comparisons for the (a) mean velocity, (b) Reynolds shear stress and (c) r.m.s. velocity fluctuations.
Appendix D. Computational cost
The simulations were performed on NASA’s Pleiades Supercomputer on the Broadwell compute nodes consisting of E5-2680v4 Intel Xeon processors at 
$2.4$ GHz. For the channel flow computations on the finest mesh G3 consisting of 
$32 \times 32 \times 32$ elements with about 
$0.885$ million degrees of freedom, the dynamic slip-wall model takes about 
$0.135$ s of wall time per time step on 
$512$ processors. For the 
$Re_\tau \approx 10\,000$ case, the dynamic slip-wall model requires a wall time of about 
$11.8$ mins for a single flow through (
$= L_x /U_b$) on the grid G3. On the other hand, for the periodic hill cases using the fine grid consisting of 
$ 75 \times 36 \times 15$ elements with about 
$1.1$ million degrees of freedom, the wall time required by the dynamic slip-wall model per time step is approximately 
$ 0.088$s on 
$ 3330$ processors. For the 
$Re_b = 37\,000$ case, the wall time required for a single flow through is about 
$22$ mins. The EQWM requires a similar time per time step as that of the dynamic slip-wall model for the channel flow and periodic hill cases on identical grids. On the other hand, a static slip-wall model using an arbitrary constant value of the slip length takes about 
$1\,\%$ less time per time step for the channel flow on the G3 mesh and about 
$7\,\%$ less time per time step for the periodic hill case on the G2 grid in comparison to the dynamic slip-wall model on an identical number of processors.
Open access
