2.1 Theoretical motivation

Let us consider a wall-bounded flow. For the Navier–Stokes equations, the velocity at the wall is given by the no-slip boundary condition (Stokes et al. Reference Stokes1901)

(2.2) $$\begin{eqnarray}u_{i}|_{w}=0,\quad i=1,2,3.\end{eqnarray}$$

If we interpret LES as the solution of the filtered Navier–Stokes equations (Leonard Reference Leonard1975), the filtering operation in the wall-normal direction will result, in general, in non-zero velocities at the wall. For wall-resolved LES, where the effective filter size near the wall is small, $\bar{u}_{i}|_{w}$ can still be approximated by the no-slip boundary condition (Ghosal & Moin Reference Ghosal and Moin1995). However, when the filter size is large or the near-wall resolution is coarse, such as in wall-modelled LES, a modified wall boundary condition different from the usual no-slip is required for the three velocity components. Consider a one-dimensional symmetric filter kernel ${\mathcal{G}}(\unicode[STIX]{x1D712})$ with non-zero filter size $\unicode[STIX]{x1D6E5}_{{\mathcal{G}}}$ defined by its second moment

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{{\mathcal{G}}}^{2}=\int _{-\infty }^{\infty }{\mathcal{G}}(\unicode[STIX]{x1D712}^{\prime }){\unicode[STIX]{x1D712}^{\prime }}^{2}\,\text{d}\unicode[STIX]{x1D712}^{\prime }.\end{eqnarray}$$

Then, far from the wall ( $x_{2}\gg \unicode[STIX]{x1D6E5}_{{\mathcal{G}}}$ ), where $x_{2}=0$ is the location of the wall, the filtered velocity field can be computed as

(2.4) $$\begin{eqnarray}\bar{u}_{i}(\boldsymbol{x})=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }{\mathcal{G}}(-x_{1}+\unicode[STIX]{x1D712}_{1}^{\prime }){\mathcal{G}}(-x_{2}+\unicode[STIX]{x1D712}_{2}^{\prime }){\mathcal{G}}(-x_{3}+\unicode[STIX]{x1D712}_{3}^{\prime })u_{i}(\unicode[STIX]{x1D74C}^{\prime })\,\text{d}\unicode[STIX]{x1D74C}^{\prime },\end{eqnarray}$$

where $\boldsymbol{x}=(x_{1},x_{2},x_{3})$ and $\unicode[STIX]{x1D74C}^{\prime }=(\unicode[STIX]{x1D712}_{1}^{\prime },\unicode[STIX]{x1D712}_{2}^{\prime },\unicode[STIX]{x1D712}_{3}^{\prime })$ . However, in the near-wall region, the filter kernel in the wall-normal direction has a functional dependence on the wall-normal distance, which becomes prevalent for $x_{2}\rightarrow 0$ , as depicted in figure 2. The new filter operator restricted by the wall then becomes

(2.5) $$\begin{eqnarray}\bar{u}_{i}(\boldsymbol{x})=\int _{-\infty }^{\infty }\int _{0}^{\infty }\int _{-\infty }^{\infty }{\mathcal{G}}(-x_{1}+\unicode[STIX]{x1D712}_{1}^{\prime }){\mathcal{G}}^{\ast }(-x_{2}+\unicode[STIX]{x1D712}_{2}^{\prime };x_{2}){\mathcal{G}}(-x_{3}+\unicode[STIX]{x1D712}_{3}^{\prime })u_{i}(\unicode[STIX]{x1D74C}^{\prime })\,\text{d}\unicode[STIX]{x1D74C}^{\prime }\end{eqnarray}$$

where

(2.6) $$\begin{eqnarray}{\mathcal{G}}^{\ast }(\unicode[STIX]{x1D712};x_{2})=\frac{{\mathcal{G}}(\unicode[STIX]{x1D712})}{\displaystyle \int _{-x_{2}}^{\infty }{\mathcal{G}}(\unicode[STIX]{x1D711})\,\text{d}\unicode[STIX]{x1D711}}\end{eqnarray}$$

is the effective (rescaled) kernel at a wall-normal distance  $x_{2}$ .

Figure 2. Sketch of the change in effective kernel approaching the wall.

Assuming a no-slip boundary condition for the unfiltered velocity and approximating the velocity profile near the wall as a Taylor expansion $u_{i}(\boldsymbol{x})=\sum _{m=1}^{p}a_{m}x_{2}^{m}$ , the filtered velocity can be expressed as

(2.7) $$\begin{eqnarray}\bar{u}_{i}|_{w}=\mathop{\sum }_{m=1}^{p}a_{m}M_{{\mathcal{G}}_{0}^{\ast }}^{(m)},\end{eqnarray}$$

where $M_{{\mathcal{G}}_{0}^{\ast }}^{(m)}$ is the $m$ th moment of ${\mathcal{G}}^{\ast }(\unicode[STIX]{x1D712};0)$ , the effective filter at the wall. From (2.7), it can be shown that

(2.8) $$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}=\mathop{\sum }_{m=1}^{p}-2a_{m}{\mathcal{G}}(0)M_{{\mathcal{G}}_{0}^{\ast }}^{(m)}+ma_{m}M_{{\mathcal{G}}_{0}^{\ast }}^{(m-1)}.\end{eqnarray}$$

Since $M_{{\mathcal{G}}_{0}^{\ast }}^{(m)}\sim \unicode[STIX]{x1D6E5}_{{\mathcal{G}}}^{m}$ and ${\mathcal{G}}(0)\sim \unicode[STIX]{x1D6E5}_{{\mathcal{G}}}^{-1}$ , we can use (2.7) and (2.8) to obtain a second-order approximation of the boundary condition of the form

(2.9) $$\begin{eqnarray}\bar{u}_{i}|_{w}=l\left.\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}x_{2}}\right|_{w},\end{eqnarray}$$

where $l=M_{{\mathcal{G}}_{0}^{\ast }}^{(1)}/(1-2{\mathcal{G}}(0)M_{{\mathcal{G}}_{0}^{\ast }}^{(1)})$ .

This boundary condition is exact for a linear velocity profile but is expected to deteriorate as the linear approximation is no longer valid for $x_{2}<\unicode[STIX]{x1D6E5}_{{\mathcal{G}}}$ . In this case, the second- and higher-order terms excluded in (2.9) may result in different slip lengths and extra terms for each velocity component as in (2.1) to achieve an accurate representation of the flow at the wall. The particular expressions for $l_{i}$ and $v_{i}$ depend formally on the filter shape, size and instantaneous configuration of the filtered velocity vector at the wall.

Equation (2.9) motivates the use of the slip boundary condition for wall-modelled LES. Nevertheless, it is important to highlight a few remarks regarding the derivation and the consistency of the slip condition. The first observation is that in the case of explicitly filtered LES (Lund & Kaltenbach Reference Lund and Kaltenbach1995; Lund Reference Lund2003; Bose Reference Bose2012; Bae & Lozano-Durán Reference Bae and Lozano-Durán2017) with a well-defined filter operator, the filter size is a given function of the wall-normal distance, and the slip lengths and velocities can be computed explicitly. However, in the present study, we focus on traditional implicitly filtered LES, where the filter operator is not distinctly defined and, consequently, neither is the filter size, typically assumed to be proportional to the grid size. This supposed relation between the filter and grid sizes is not always valid (Lund Reference Lund2003; Silvis et al. Reference Silvis, Trias, Abkar, Bae, Lozano-Durán and Verstappen2016) and worsens close to the wall. Therefore, in the near-wall region, it is reasonable to assume that the effective filter size is an unknown function of the wall-normal distance, and the slip lengths and velocities must be modelled as they cannot be computed explicitly. As a final remark, note that commutation of the filter and derivative operators is necessary to formally derive the LES equations, which, in turn, entails a constant-in-space filter size or a filter operator that is constructed to be commutative (Marsden, Vasilyev & Moin Reference Marsden, Vasilyev and Moin2002). This condition is not met by (2.1), but given that the filter size for implicitly filtered LES is an unknown function of space, we also neglect terms arising from commutation errors.

2.2 A priori evaluation

A priori testing of the slip boundary condition is conducted to assess the validity of (2.1) in filtered DNS data of turbulent channel flow from Del Álamo et al. (Reference Del Álamo, Jiménez, Zandonade and Moser2004), Hoyas & Jiménez (Reference Hoyas and Jiménez2006) and Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014a ) at friction Reynolds numbers of $Re_{\unicode[STIX]{x1D70F}}\approx 950$ , 2000 and 4200, respectively.

In the following, $u_{\unicode[STIX]{x1D70F}}$ is the friction velocity, $\unicode[STIX]{x1D708}$ is the kinematic viscosity and the channel half-height is denoted by $\unicode[STIX]{x1D6FF}$ . Wall units are defined in terms of $\unicode[STIX]{x1D708}$ and $u_{\unicode[STIX]{x1D70F}}$ and denoted by the superscript $+$ and outer units in terms of $\unicode[STIX]{x1D6FF}$ and $u_{\unicode[STIX]{x1D70F}}$ . In each case, DNS velocity vector is filtered in the three spatial directions with a box-filter with filter size equal to $\unicode[STIX]{x1D6E5}_{1}\times \unicode[STIX]{x1D6E5}_{2}\times \unicode[STIX]{x1D6E5}_{3}$ . The resulting filtered data contain $\bar{u}_{i}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{i}/\unicode[STIX]{x2202}x_{2}|_{w}$ , which can be used to test the accuracy of (2.1) by computing their joint probability density function (PDF).

Figure 3. Joint PDF of (a) $\bar{u}_{1}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{1}/\unicode[STIX]{x2202}x_{2}|_{w}$ , (b) $\bar{u}_{2}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{2}/\unicode[STIX]{x2202}x_{2}|_{w}$ and (c) $\bar{u}_{3}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{3}/\unicode[STIX]{x2202}x_{2}|_{w}$ for box-filtered DNS with $\unicode[STIX]{x1D6E5}_{i}=0.01\unicode[STIX]{x1D6FF}$ (○), $0.02\unicode[STIX]{x1D6FF}$ (▿) and $0.03\unicode[STIX]{x1D6FF}$ (♢) at $Re_{\unicode[STIX]{x1D70F}}\approx 950$ . For each probability distribution the contours are 50 % and 95 %. The straight-dashed lines are obtained by the least-squares fitting to the joint PDF. (d) The $l_{1}$ dependence on $Re_{\unicode[STIX]{x1D70F}}$ with $\unicode[STIX]{x1D6E5}_{2}=0.05\unicode[STIX]{x1D6FF}$ (black, bottom) and $\unicode[STIX]{x1D6E5}_{2}=0.1\unicode[STIX]{x1D6FF}$ (blue, top) calculated from box-filtered DNS channel flow data (○), and estimation from box-filtered logarithmic layer approximation $l_{1}/\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6E5}_{2}/(2\unicode[STIX]{x1D6FF})[\log (Re_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6E5}_{2}/(2\unicode[STIX]{x1D6FF}))-1]+\unicode[STIX]{x1D705}B\unicode[STIX]{x1D6E5}_{2}/(2\unicode[STIX]{x1D6FF})$ with $\unicode[STIX]{x1D705}=0.41$ and $B=5.3$ (dashed line).

The joint PDF for $Re_{\unicode[STIX]{x1D70F}}\approx 950$ with $\unicode[STIX]{x1D6E5}_{1}=\unicode[STIX]{x1D6E5}_{2}=\unicode[STIX]{x1D6E5}_{3}$ are plotted in figure 3(a–c). The results show, on average, a linear correlation between $\bar{u}_{i}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{i}/\unicode[STIX]{x2202}x_{2}|_{w}$ , which supports the suitability of the slip boundary condition given in (2.1) with $v_{i}=0$ . However, the spread of the joint PDFs increases with increasing filter size. A trend similar to that of figure 3(a–c) appears when increasing the Reynolds number for a constant filter size $\unicode[STIX]{x1D6E5}_{i}/\unicode[STIX]{x1D6FF}$ (not shown), which implies that the second-order approximation deteriorates as the filter size increases in wall units. However, despite this scaling, a linear relationship between $\bar{u}_{i}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{i}/\unicode[STIX]{x2202}n|_{w}$ is still satisfied on average, and it will be shown in § 3 that this is enough to obtain accurate predictions of the mean velocity profile in actual LES.

Finally, figure 3(d) shows the streamwise slip length $l_{1}$ as a function of Reynolds number for a constant $\unicode[STIX]{x1D6E5}_{i}/\unicode[STIX]{x1D6FF}$ . The slip length was computed as the average ratio of $\bar{u}_{1}|_{w}$ and $\unicode[STIX]{x2202}\bar{u}_{1}/\unicode[STIX]{x2202}x_{2}|_{w}$ . The plot shows that the dependence of the streamwise slip length on Reynolds number is stronger for smaller $Re_{\unicode[STIX]{x1D70F}}$ . This behaviour can be explained by the relative thickness of the filter size and the buffer layer. When the ratio is of $O(1)$ , the filtered velocity at the wall takes into account contributions from the buffer layer, and $l_{1}$ is expected to be sensitive to changes in Reynolds number. However, when the buffer layer is a small fraction of the filter size, most of the contribution to $l_{1}$ comes from the logarithmic layer, which has a universal behaviour with $Re_{\unicode[STIX]{x1D70F}}$ . Neglecting the effect of the buffer layer, the approximate functional dependence of the slip length on Reynolds number can be estimated from the logarithmic velocity profile,

(2.10) $$\begin{eqnarray}\frac{\langle u_{1}\rangle }{u_{\unicode[STIX]{x1D70F}}}=\frac{1}{\unicode[STIX]{x1D705}}\log (x_{2}^{+})+B,\end{eqnarray}$$

where $\langle \cdot \rangle$ denotes ensemble average in the homogeneous directions and time, $\unicode[STIX]{x1D705}$ is the von Kármán constant and $B$ is the intercept constant. In the limit of high Reynolds numbers, the box-filtered streamwise velocity and its wall-normal derivative can be estimated by assuming a logarithmic law in the entire near-wall region and integrating (2.10). This gives an approximation for the average streamwise slip length:

(2.11) $$\begin{eqnarray}\frac{l_{1}}{\unicode[STIX]{x1D6FF}}\sim \frac{\langle \bar{u}_{1}/u_{\unicode[STIX]{x1D70F}}\rangle }{\unicode[STIX]{x2202}\langle \bar{u}_{1}/u_{\unicode[STIX]{x1D70F}}\rangle /\unicode[STIX]{x2202}(x_{2}/\unicode[STIX]{x1D6FF})}\sim \frac{\unicode[STIX]{x1D6E5}_{2}}{2\unicode[STIX]{x1D6FF}}\left[\log \left(Re_{\unicode[STIX]{x1D70F}}\frac{\unicode[STIX]{x1D6E5}_{2}}{2\unicode[STIX]{x1D6FF}}\right)-1\right]+\unicode[STIX]{x1D705}B\frac{\unicode[STIX]{x1D6E5}_{2}}{2\unicode[STIX]{x1D6FF}}\end{eqnarray}$$

(dashed lines in figure 3 d), which predicts a weak $Re_{\unicode[STIX]{x1D70F}}$ dependence for large Reynolds numbers. It is important to remark that $l_{1}$ from the figure is only an estimation from a priori testing and the particular values are not expected to work in an actual LES, although we expect the trends to be relevant.

2.3 Consistency constraints on the slip parameters

In an actual LES implementation, the choice of $l_{i}$ and $v_{i}$ must comply with the symmetries of the flow. Moreover, it is also necessary to satisfy on average the impermeability constraint of the wall to preserve the physics of the flow (more details are offered in § 3.6). Therefore, the slip boundary condition for a plane channel flow should fulfil

(2.12) $$\begin{eqnarray}\langle \bar{u}_{i}|_{w}\rangle =\left\langle \left.l_{i}\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}n}\right|_{w}\right\rangle +\langle v_{i}\rangle =0,\quad i=2,3.\end{eqnarray}$$

Equation (2.12) can be further simplified by assuming $l_{i}$ and $v_{i}$ to be constant in the homogeneous directions. Since $\langle \bar{u}_{i}|_{w}\rangle =0$ and $\langle \unicode[STIX]{x2202}\bar{u}_{i}/\unicode[STIX]{x2202}x_{2}|_{w}\rangle =0$ for $i=2,3$ , $\langle v_{2}\rangle$ and $\langle v_{3}\rangle$ must be set to zero. We can also set $\langle v_{1}\rangle =0$ without loss of generality, since its average effect can be absorbed by moving the frame of reference at constant uniform velocity. Then, the slip boundary condition consistent with the symmetries of the channel is of the form

(2.13) $$\begin{eqnarray}\bar{u}_{i}|_{w}=l_{i}\left.\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}n}\right|_{w}.\end{eqnarray}$$

When the flow is no longer homogeneous in $x_{1}$ , as in a spatially developing flat-plate boundary layer, the above arguments based on the symmetry of the channel do not hold. Then, the slip velocity $v_{i}$ can be used to impose zero mean mass flow through the walls and ensure that the boundary behaves, on average, as a non-permeable wall. Since the mass flow through a flat wall is only affected by $v_{2}$ , we can still set $v_{1}$ and $v_{3}$ to zero for simplicity.

3.1 Numerical experiments

Table 1. List of cases. The numerical experiments are labelled following the convention [SGS model]-[ $Re_{\unicode[STIX]{x1D70F}}$ ](-[other cases]). SGS models used are the dynamic Smagorinsky model (DSM), constant coefficient Smagorinsky model (SM), anisotropic minimum-dissipation model (AMD) and no model (NM). Grid resolutions different from the baseline case are denoted by c1 and c2. Three additional cases with different slip length than the baseline case are labelled s1, s2 and s3. See the text for details.

We perform a set of plane turbulent channel LES listed in table 1. The simulations are computed with a staggered second-order finite difference (Orlandi Reference Orlandi2000) and a fractional-step method (Kim & Moin Reference Kim and Moin1985) with a third-order Runge–Kutta time-advancing scheme (Wray Reference Wray1990). The code has been validated in previous studies in turbulent channel flows (Lozano-Durán & Bae Reference Lozano-Durán and Bae2016; Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2018) and flat-plate boundary layers (Lozano-Durán, Hack & Moin Reference Lozano-Durán, Hack and Moin2018). The size of the channel is $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}\times 2\unicode[STIX]{x1D6FF}\times \unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ in the streamwise, wall-normal and spanwise directions, respectively. It has been shown that this domain size is sufficient to accurately predict one-point statistics for $Re_{\unicode[STIX]{x1D70F}}$ up to 4200 (Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014a ). Periodic boundary conditions are imposed in the streamwise and spanwise directions. The eddy viscosity $\unicode[STIX]{x1D708}_{t}$ is computed at the cell centres and the values at the wall are obtained by assuming Neumann boundary conditions for the discretised $\unicode[STIX]{x1D708}_{t}$ , which is motivated by the fact that for coarse grid resolutions the SGS contribution at the wall must be non-zero. The channel flow is driven by imposing a constant mean pressure gradient, and all simulations were run for at least 100 eddy turnover times, defined as $\unicode[STIX]{x1D6FF}/u_{\unicode[STIX]{x1D70F}}$ , after transients.

The slip boundary condition from (2.13) is used on the top and bottom walls. We have tested the variability of $l_{i}$ in time by oscillating $l_{i}(t)$ with different amplitudes and frequencies around a given mean. The frequency of the oscillation considered were 0.5, 1 and 2 times the natural frequency given by the size of the grid and $u_{\unicode[STIX]{x1D70F}}$ , and the amplitudes imposed were up to 0.5 times the value of the mean. The different cases resulted in almost identical one-point statistics as those obtained with a constant $l_{i}$ of the same mean with the relative difference in the resulting wall stress less than 0.5 % for all cases. Thus, $l_{i}$ will be fixed to a constant value in both homogeneous directions and time for the remainder of the section.

We take as a baseline case the friction Reynolds number $Re_{\unicode[STIX]{x1D70F}}\approx 2000$ with a uniform grid resolution of $128\times 40\times 64$ in the $x_{1}$ , $x_{2}$ and $x_{3}$ directions, respectively. The grid size in outer units is $0.050\unicode[STIX]{x1D6FF}$ in the three directions, and follows the recommendations by Chapman (Reference Chapman1979) for resolving the large eddies in the outer portion of the boundary layer. The baseline SGS model used is the dynamic Smagorinsky model (DSM) (see Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992). The baseline slip lengths are $l_{i}=0.008\unicode[STIX]{x1D6FF}$ , $i=1,2,3$ for reasons given in § 3.2. Three additional cases with different slip lengths, which are given in the first group of table 1, are used to study the effect of $l_{i}$ on the one-point statistics.

To study the effects of slip boundary condition on grid resolution, we define two meshes with $82\times 26\times 42$ and $64\times 20\times 32$ grid points distributed uniformly in each direction, which correspond to a uniform grid size of $0.077\unicode[STIX]{x1D6FF}$ and $0.100\unicode[STIX]{x1D6FF}$ , respectively as listed in the second group of table 1. The resolutions were chosen such that the first interior point lies in the logarithmic region and is far from the inner-wall peak of the streamwise r.m.s. velocity. The intention is to avoid capturing (even partially) the dynamic cycle in the buffer layer, since the wall-normal lengths of the near-wall vortices and streaks scale in viscous units, and that scaling is incompatible with the computational efficiency pursued in wall-modelled LES. The range of grid resolutions is limited due to the fact that the outer layer still needs to be resolved by the LES. However, it will be shown in § 3.5 that the selected range of grid resolutions is sufficient to show the sensitivity of the slip lengths to the grid resolution.

To investigate the effect of the Reynolds number, we consider three cases DSM-950, DSM-2000 and DSM-4200, which constitute the third group of table 1. The sensitivity to the SGS model will be assessed by comparing results from DSM, constant-coefficient Smagorinsky model (SM) (see Smagorinsky Reference Smagorinsky1963) without a damping function at the wall, the anisotropic minimum-dissipation model (AMD) (see Rozema et al. Reference Rozema, Bae, Moin and Verstappen2015) and cases without an SGS model (NM), given in the fourth group of table 1. Finally, LES results will be compared with DNS data from Del Álamo et al. (Reference Del Álamo, Jiménez, Zandonade and Moser2004), Hoyas & Jiménez (Reference Hoyas and Jiménez2006) and Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014a ).

3.2 Control of the wall stress and optimal slip lengths

In an LES of channel flow with the slip boundary condition (2.13), the wall stress is given by

(3.1) $$\begin{eqnarray}\langle \unicode[STIX]{x1D70F}_{w}\rangle =\unicode[STIX]{x1D708}\left.\left\langle \frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle -\langle \bar{u}_{1}\bar{u}_{2}|_{w}\rangle -\langle \unicode[STIX]{x1D70F}_{12}^{\mathit{SGS}}|_{w}\rangle ,\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{w}$ is the stress at the wall, $\unicode[STIX]{x1D70F}_{12}^{\mathit{SGS}}$ is the tangential SGS stress tensor and $\langle \bar{u}_{1}\bar{u}_{2}|_{w}\rangle$ is the result of the non-zero velocity provided by the slip condition. The slip lengths can be explicitly introduced by substituting $\bar{u}_{1}\bar{u}_{2}$ from (2.13) such that

(3.2) $$\begin{eqnarray}\langle \unicode[STIX]{x1D70F}_{w}\rangle =\unicode[STIX]{x1D708}\left.\left\langle \frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle -\left.\left\langle l_{1}l_{2}\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\frac{\unicode[STIX]{x2202}\bar{u}_{2}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle -\langle \unicode[STIX]{x1D70F}_{12}^{\mathit{SGS}}|_{w}\rangle ,\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{12}^{\mathit{SGS}}$ may also depend on the slip lengths. Therefore, the wall stress (and hence the mass flow) can be controlled by the proper choice of slip lengths. This is an important property of the slip boundary condition, and it is illustrated in figure 4. For coarse LES with no-slip boundary conditions, the near-wall region cannot be accurately computed owing to the inadequacy of the current SGS models and large numerical errors in the near-wall region, even if the resolution is sufficient to resolve the outer layer eddies. This mainly results in under- or over-predictions of the wall stress, among other effects, and the shift of mean velocity with respect to DNS. Figure 4 shows the mean streamwise velocity profile for cases DSM-2000 and DSM-2000-s[1–3]. The results reveal that increasing $l_{1}$ (at constant $l_{2}$ ) moves up the mean velocity profile by 8 %, while increasing $l_{2}$ (at constant $l_{1}$ ) have the opposite effect and decreases the mean by 15 %. Although not shown, it was tested that varying $l_{3}$ has a second-order effect on the mean velocity profile when compared to changes of the same order in $l_{1}$ and $l_{2}$ . For instance, the change in mean velocity profile is 1.2 % when $l_{3}$ is changed from $0.008\unicode[STIX]{x1D6FF}$ to $0.004\unicode[STIX]{x1D6FF}$ for case DSM-2000. The result is not totally unexpected since $\bar{u}_{1}$ and $\bar{u}_{2}$ are active components of the mean streamwise momentum balance in a channel flow (see (3.2)), while $\bar{u}_{3}$ enters only indirectly. All calculations in the present study have been performed with $l_{3}$ equal to  $l_{1}$ .

Figure 4. Mean streamwise velocity profile as a function of (a) outer units and (b) wall units for DSM-2000 $(l_{1},l_{2})=(0.008\unicode[STIX]{x1D6FF},0.008\unicode[STIX]{x1D6FF})$ , ○; DSM-2000-s1 $(l_{1},l_{2})=(0.008\unicode[STIX]{x1D6FF},0.004\unicode[STIX]{x1D6FF})$ , ▿; DSM-2000-s2 $(l_{1},l_{2})=(0.004\unicode[STIX]{x1D6FF},0.008\unicode[STIX]{x1D6FF})$ , ▫; DSM-2000-s3 $(l_{1},l_{2})=(0.097\unicode[STIX]{x1D6FF},0.045\unicode[STIX]{x1D6FF})$ , $\times$ ; and DNS (dashed line).

The observations in figure 4 may be explained in terms of the mean streamwise momentum balance at the wall and non-zero streamwise slip. With respect to the former, increasing $l_{2}$ enhances the $\langle \bar{u}_{1}\bar{u}_{2}\rangle$ contribution at the wall, which is translated into a lower mean velocity profile due to the higher momentum drain at the boundaries. The same argument applies when increasing $l_{1}$ . However, higher $l_{1}$ also implies larger slip in $x_{1}$ , which overcomes the previous momentum drain, and the resulting net effect is an increase of the mean mass flow. For the laminar Poiseuille flow with the slip boundary condition, the shift in the mean velocity profile can be computed analytically and is shown to be proportional to  $l_{1}$ .

The duality between the streamwise and wall-normal slip lengths makes it possible to always achieve the correct wall stress by an appropriate selection of $(l_{1},l_{2})$ , which we define as optimal slip lengths. The optimal slip lengths are effectively computed by running an LES with slip boundary condition using (3.2) with $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D70F}_{w}^{\mathit{DNS}}$ as a constraint, and then averaging in time the values obtained for $l_{1}$ and $l_{2}$ . Note that we have the freedom to impose the ratio $l_{1}/l_{2}$ , and the optimal slip lengths are not unique. Two examples are shown in figure 4. It is also important to remark that the control of the mean velocity profile is not possible in general without wall-normal transpiration. In particular, if an LES with $l_{1}/\unicode[STIX]{x1D6FF}=l_{2}/\unicode[STIX]{x1D6FF}=0$ (no-slip) already over-predicts the mean velocity profile with respect to DNS, the only possible outcome of increasing $l_{1}/\unicode[STIX]{x1D6FF}$ while maintaining $l_{2}/\unicode[STIX]{x1D6FF}=0$ is a positive shift of the mean velocity profile. Since our experience shows that negative values of $l_{1}$ will result in an unstable solution, the conclusion is that the correct mean velocity profile cannot be achieved in this case unless $l_{2}\neq 0$ , making the wall-normal slip length indispensable.

3.3 Prediction of the logarithmic layer

A second observation from figure 4 is that the shape of the mean velocity profile remains roughly constant for different slip lengths, and changes in $l_{1}$ and $l_{2}$ are mainly responsible for a shift along the mean velocity axis. We would like to connect the previous observation with the classic logarithmic profile for the mean streamwise velocity given in (2.10). Assuming that the filter operation does not alter the logarithmic shape of $\langle \bar{u}_{1}\rangle$ for the typical filter sizes (or grid resolutions), (2.10) should also hold for LES. However, it is not clear whether this would be the case for an actual LES. For example, Millikan’s asymptotic matching argument (Millikan Reference Millikan1938) requires a scale separation that tends to infinity as the Reynolds number increases, which is not the case in wall-modelled LES as the length scales are fixed in outer units. Other arguments, such as the Prandtl’s mixing length hypothesis (Prandtl Reference Prandtl1925) would suggest that the correct wall-normal mixing of the flow should be obtained in order to recover the logarithmic profile. Alternatively, from the point of view of Townsend’s attached eddy hypothesis (Townsend Reference Townsend1980), the flow from the LES should be populated by a self-similar hierarchy of eddies with sizes proportional to the wall distance and the proper number of eddies per unit area. In all cases, the SGS model plays a non-negligible role in fulfilling these conditions, especially at high Reynolds numbers and coarse grids. As a consequence, not all SGS models are expected to recover the correct shape of the mean velocity profile, and this is further discussed in § 3.5.

The prominent role of the SGS model in the correct representation of the logarithmic layer can be seen from the integrated mean streamwise momentum balance for filtered velocities

(3.3) $$\begin{eqnarray}\langle \bar{u}_{1}^{+}\rangle (x_{2}^{+})=\underbrace{\langle \bar{u}_{1}^{+}|_{w}\rangle }_{{\sim}B}+\underbrace{x_{2}^{+}\left(1-\frac{x_{2}}{2\unicode[STIX]{x1D6FF}}\right)+\displaystyle \int _{0}^{x_{2}^{+}}\langle \bar{u}_{1}^{+}\bar{u}_{2}^{+}+\unicode[STIX]{x1D70F}_{12}^{\mathit{SGS}+}\rangle \,\text{d}x_{2}^{\prime +}}_{{\sim}1/\unicode[STIX]{x1D705}\log (x_{2}^{+})}.\end{eqnarray}$$

By comparing the structure of (2.10) and (3.3), it is reasonable to hypothesise that the slip boundary condition mainly influences the intercept $B$ , which is independent of $x_{2}$ , while the SGS model controls the $x_{2}$ -dependent slope $1/\unicode[STIX]{x1D705}$ , related to the wall-normal mixing of the flow by the attached eddies. Note that this is not strictly the case, and some coupling is expected between all terms in (3.3). For example, the value of the integrand at the wall will depend on the slip lengths.

3.4 Turbulence intensities and Reynolds stress contribution

The sensitivity of $\langle u_{i}^{\prime 2}\rangle ^{1/2}$ to the choice of slip lengths is examined in figure 5(a). Three cases DSM-2000, DSM-2000-s1 and DSM-2000-s3 are considered, two of them supplying the correct mean velocity profile. The r.m.s. velocities are insensitive, especially away from the wall, even when the slip lengths are such that the mean velocity profile does not match that of DNS. The most noticeable difference is observed at the wall, where larger $l_{1}$ results in smaller r.m.s. values. The consequence is that even if the mean velocity profile matches that of DNS for an optimal $(l_{1},l_{2})$ , some pairs are preferred in order to avoid the near-wall under- and over-shoot in the r.m.s. velocity fluctuations. A more comprehensive study of the near-wall turbulent intensities with the slip boundary condition can be found in Bae et al. (Reference Bae, Lozano-Durán, Bose and Moin2018).

Figure 5. (a) Streamwise, spanwise and wall-normal r.m.s. velocity fluctuations (from top to bottom), and (b) mean Reynolds stress contribution for DSM-2000 $(l_{1},l_{2})=(0.008\unicode[STIX]{x1D6FF},0.008\unicode[STIX]{x1D6FF})$ , ○; DSM-2000-s1 $(l_{1},l_{2})=(0.008\unicode[STIX]{x1D6FF},0.004\unicode[STIX]{x1D6FF})$ , ▿; DSM-2000-s3 $(l_{1},l_{2})=(0.097\unicode[STIX]{x1D6FF},0.045\unicode[STIX]{x1D6FF})$ , $\times$ ; and DNS (dashed line).

Figure 5(b) stresses another important property of the slip condition. As was the case for r.m.s. velocities, the changes in the mean tangential Reynolds stresses are negligible to varying values of $l_{1}$ and $l_{2}$ , except at the walls. For the case with larger $l_{1}$ , the wall-stress contribution from the $\langle \bar{u}_{1}\bar{u}_{2}\rangle$ is roughly 50 %, and the remaining stress is then carried by the SGS and viscous terms. This can be considered an advantage compared to the classic no-transpiration condition since the SGS model, usually known to under-predict the wall stress (Jiménez & Moser Reference Jiménez and Moser2000), is not constrained to account for the resolved non-zero $\bar{u}_{i}\bar{u}_{j}$ at the wall.

Finally, the structure of the streamwise velocity at the wall for filtered DNS and wall-modelled LES is shown in figure 6. The filtered DNS data was obtained by box-filtering the streamwise velocity with filter size $\unicode[STIX]{x1D6E5}_{i}=0.050\unicode[STIX]{x1D6FF}$ , $i=1,2,3$ , that coincides with the LES grid resolution. Although our analysis is qualitative, the figures show that despite the comparable intensities, the filtered DNS is organised into more elongated streaks. Also note that for a constant $l_{1}$ , the slip boundary condition forces the velocity and its wall-normal derivative to have the same structure close to the wall that is inconsistent with box-filtered DNS. This suggests that an accurate representation of the flow structure at the wall is neither expected nor necessary in order to obtain accurate predictions of the low-order flow statistics far from the wall. This is consistent with previous studies indicating that the outer layer dynamics are relatively independent of the near-wall cycle (Del Álamo et al. Reference Del Álamo, Jiménez, Zandonade and Moser2006; Flores & Jiménez Reference Flores and Jiménez2006; Jiménez Reference Jiménez2013; Mizuno & Jiménez Reference Mizuno and Jiménez2013; Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014b ; Dong et al. Reference Dong, Lozano-Durán, Sekimoto and Jiménez2017).

Figure 6. Instantaneous snapshot of the streamwise velocity at the wall for (a) box-filtered DNS ( $Re_{\unicode[STIX]{x1D70F}}\approx 2000$ ), and (b) wall-modelled LES (DSM-2000) of channel flow for the $x_{1}{-}x_{3}$ plane. For both cases, the filter or grid size is $\unicode[STIX]{x1D6E5}_{i}/\unicode[STIX]{x1D6FF}=0.050$ , $i=1,2,3$ . Colours indicates velocity in wall units.

Figure 7. (a) Effect of SGS models on the mean velocity profile for $l_{i}=0.008\unicode[STIX]{x1D6FF}$ . (b) The mean velocity profiles have been shifted to compare the shapes of the mean velocity profile, where the shift is given by $\unicode[STIX]{x0394}u=u_{1}^{+}(\unicode[STIX]{x1D6FF})-u_{1}^{+\mathit{DNS}}(\unicode[STIX]{x1D6FF})$ . Dynamic Smagorinsky model (○), constant coefficient Smagorinsky model (▫), anisotropic minimum-dissipation model (♢) and no model (▿) are given for the turbulent channel with $Re_{\unicode[STIX]{x1D70F}}\approx 2000$ . DNS (dashed line).

3.5 Sensitivity to SGS model, Reynolds number and grid resolution

In this section, we study the effect of the SGS model, Reynolds number and grid resolution on the mean velocity profile for the slip boundary condition. The discussion is necessary for understanding the most relevant sensitivities of wall models based on the slip boundary condition.

Figure 7 shows the sensitivity of the mean velocity to different SGS models for DSM-2000, SM-2000, AMD-2000 and NM-2000. In all of the cases, the slip lengths are fixed and equal to $0.008\unicode[STIX]{x1D6FF}$ such that the velocity profile at the centre of the channel for DSM-2000 matches the DNS data. Note that this particular choice is arbitrary, and that alternative values of slip lengths could be selected to find the best match between SM-2000, AMD-2000 or NM-2000 and DNS. However, the results below are independent of this choice, since the relative shift between cases is barely affected. The results in figure 7 reveal that not only the shape but also the mean mass flow, and thus the optimal slip lengths for each SGS model, are impacted by the SGS model at grid resolutions typical of wall-modelled LES. Regarding the shape of $\langle \bar{u}_{1}\rangle$ (figure 7 b), for low-dissipation SGS models (e.g. NM), the flow becomes more turbulent, causing the mean velocity profile to flatten due to the enhanced mixing. On the other hand, for highly dissipative SGS models (e.g. SM), the shape approaches a parabolic profile, that is closer to the laminar solution.

The effect of each SGS model on the mass flow rate in wall units can be understood by considering the definition of the friction velocity,

(3.4) $$\begin{eqnarray}u_{\unicode[STIX]{x1D70F}}^{2}=-\langle \bar{u}_{1}\bar{u}_{2}|_{w}\rangle +\left\langle \left.\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle +\left\langle \left.\unicode[STIX]{x1D708}_{t}\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle .\end{eqnarray}$$

For a channel flow driven by a constant mass flow rate, the last term in (3.4) is zero for LES without an SGS model, which results in lower $u_{\unicode[STIX]{x1D70F}}$ and therefore a positive shift of the mean velocity profile scaled in wall units. For non-zero eddy viscosity, the mean SGS stress at the wall will contribute to increase $u_{\unicode[STIX]{x1D70F}}$ , creating a negative shift in the mean velocity profile in wall units. The actual impact of the SGS model on $u_{\unicode[STIX]{x1D70F}}$ is more intricate owing to the coupling between $\unicode[STIX]{x1D708}_{t}$ and the flow velocities. However, the qualitative behaviour of $u_{\unicode[STIX]{x1D70F}}$ described above still holds.

Conversely, the effects on the mean velocity profile can also be explained for a channel flow driven by a constant pressure gradient. In this case, the left-hand side of (3.4) is fixed. Hence, variations in the SGS stress at the wall must be compensated by variations in the Reynolds and viscous stress terms. We have observed that these changes are balanced by the viscous stress, $\unicode[STIX]{x1D708}\unicode[STIX]{x2202}\bar{u}_{1}/\unicode[STIX]{x2202}x_{2}|_{w}$ , rather than by the Reynolds stress term. The variation in the mass flow can then be understood through the slip boundary condition, where larger $\unicode[STIX]{x2202}\bar{u}_{1}/\unicode[STIX]{x2202}x_{2}|_{w}$ implies a larger slip at the wall and, hence, higher mass flow.

Figure 8. Eddy viscosity, $\unicode[STIX]{x1D708}_{t}$ , as a function of wall-normal height with $l_{i}=0.008\unicode[STIX]{x1D6FF}$ for: (a) DSM-2000 (○) and AMD-2000 (▿); and (b) DSM-2000 (○), DSM-4200 (▿) and DMS-2000-c2 (▫).

Figure 8(a) shows the eddy viscosity as a function of wall-normal height for various SGS models. In particular, the eddy viscosity of DSM and AMD model are comparable far from the wall, consistent with the shape of the mean velocity profile in the outer region shown in figure 7(b). However, $\unicode[STIX]{x1D708}_{t}$ differs notably for DSM and AMD model in the first two grid points off the wall, leading to the differences in mean velocity profile observed in figure 7(a). Hence, the results above are indicative of the fact that different SGS models demand different optimal slip lengths.

The grid resolution and Reynolds number sensitivity are studied in figure 9, again maintaining constant slip lengths. Regarding the resolution, coarsening the grid increases the mass flow. This phenomenon is also observed in LES with no-slip boundary condition and, in the present case, is probably related to an inconsistency between the choice of slip lengths and the wall-normal momentum flux provided by the SGS model. Figure 9(b) shows a weak dependence of the mean velocity profile on the Reynolds number. The most notable observation is the under-estimation of the mass flow for the lowest Reynolds number, but overall the optimal slip lengths are quite insensitive to $Re_{\unicode[STIX]{x1D70F}}$ , consistent with our analysis in § 2.2 regarding the effect of $Re_{\unicode[STIX]{x1D70F}}$ . Finally, increasing the Reynolds number or coarsening the grid resolution augments the eddy viscosity for DSM (figure 8 b), which is consistent with the expected behaviour from SGS models.

Figure 9. Effect of (a) the grid resolution and (b) Reynolds number on the mean velocity profile for $l_{i}=0.008\unicode[STIX]{x1D6FF}$ : (a) $\unicode[STIX]{x1D6E5}_{i}/\unicode[STIX]{x1D6FF}=0.050$ (○), 0.077 (▿), and 0.100 (▫); (b) $Re_{\unicode[STIX]{x1D70F}}\approx 950$ (▿), $Re_{\unicode[STIX]{x1D70F}}\approx 2000$ (○), $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ , (▫). DNS (dashed line).

3.6 The role of slip velocity in imposing zero mean mass flow through the walls

As discussed in § 2.3, we require the slip velocities to guarantee no net mass flow through the walls for flows that are inhomogeneous in the wall-parallel direction. The requirement is demonstrated here in an LES of a zero-pressure-gradient flat-plate turbulent boundary layer.

The numerical method is similar to that of the channel flow presented in § 3.1 with the exception of the boundary conditions and the Poisson solver, which was modified to take into account the non-periodic boundary conditions in the streamwise direction. The simulation ranges from $Re_{\unicode[STIX]{x1D703}}\approx 1000$ to 10 000, where $Re_{\unicode[STIX]{x1D703}}$ is the Reynolds number based on the momentum thickness. This range is comparable to the boundary layer simulation by Sillero, Jiménez & Moser (Reference Sillero, Jiménez and Moser2013) that will be used for comparisons.

The slip boundary condition from (2.1) is used at the wall, located at $x_{2}=0$ . In the top plane, we impose $u_{1}=U_{\infty }$ (free-stream velocity), $u_{3}=0$ and $u_{2}$ estimated from the known experimental growth of the displacement thickness for the corresponding range of Reynolds numbers as in Jiménez et al. (Reference Jiménez, Hoyas, Simens and Mizuno2010). This controls the average streamwise pressure gradient, whose nominal value is set to zero. The turbulent inflow is generated by the recycling scheme of Lund, Wu & Squires (Reference Lund, Wu and Squires1998), in which the velocities from a reference downstream plane, $x_{ref}$ , are used to synthesise the incoming turbulence. The reference plane is located well beyond the end of the inflow region to avoid spurious feedback (Nikitin Reference Nikitin2007; Simens et al. Reference Simens, Jiménez, Hoyas and Mizuno2009). In our case, $x_{ref}/\unicode[STIX]{x1D703}_{0}=890$ , where $\unicode[STIX]{x1D703}_{0}$ is the momentum thickness at the inlet. A convective boundary condition is applied at the outlet with convective velocity $U_{\infty }$ (Pauley, Moin & Reynolds Reference Pauley, Moin and Reynolds1990) and small corrections to enforce global mass conservation (Simens et al. Reference Simens, Jiménez, Hoyas and Mizuno2009). The spanwise direction is periodic. The length, height and width of the simulated box are $L_{x}=1060\unicode[STIX]{x1D703}_{avg}$ , $L_{y}=18\unicode[STIX]{x1D703}_{avg}$ and $L_{z}=35\unicode[STIX]{x1D703}_{avg}$ , where $\unicode[STIX]{x1D703}_{avg}\approx 2.12\unicode[STIX]{x1D703}_{0}$ denotes the momentum thickness averaged along the streamwise coordinate. This domain size is similar to those used in previous studies (Jiménez et al. Reference Jiménez, Hoyas, Simens and Mizuno2010; Schlatter & Örlü Reference Schlatter and Örlü2010; Sillero et al. Reference Sillero, Jiménez and Moser2013). The streamwise and spanwise resolutions are $\unicode[STIX]{x1D6E5}_{1}/\unicode[STIX]{x1D6FF}=0.05$ ( $\unicode[STIX]{x1D6E5}_{1}^{+}=118$ ) and $\unicode[STIX]{x1D6E5}_{3}/\unicode[STIX]{x1D6FF}=0.04$ ( $\unicode[STIX]{x1D6E5}_{3}^{+}=84.3$ ) at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ . The number of wall-normal grid points per boundary layer thickness is chosen to be ${\sim}20$ at the inlet, which is in line with the channel flow simulations in the previous sections. The grid is slightly stretched in the wall-normal direction with minimum $\unicode[STIX]{x1D6E5}_{2}/\unicode[STIX]{x1D6FF}=0.01$ ( $\unicode[STIX]{x1D6E5}_{2}^{+}=20.8$ ). All computations were run with $CFL=0.5$ and for 50 washouts after transients.

The slip lengths are computed to match the empirical friction coefficient, $C_{f}$ , from White & Corfield (Reference White and Corfield2006). The connection between the slip parameters and the friction coefficient is

(3.5) $$\begin{eqnarray}\frac{1}{2}U_{\infty }^{2}\langle C_{f}\rangle =\unicode[STIX]{x1D708}\left.\left\langle \frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle -\left.\left\langle l_{1}l_{2}\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\frac{\unicode[STIX]{x2202}\bar{u}_{2}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle +\left.\left\langle l_{1}v_{2}\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right|_{w}\right\rangle -\langle \unicode[STIX]{x1D70F}_{12}^{\mathit{SGS}}|_{w}\rangle ,\end{eqnarray}$$

which is equivalent to (3.1) with the slip boundary condition applied to the $u_{1}u_{2}$ term. The slip lengths are now a function of the streamwise coordinate to take into account the inhomogeneity of the flow in $x_{1}$ . To ensure numerical stability, exponential filtering in time with filter size $0.2\unicode[STIX]{x1D6FF}/u_{\unicode[STIX]{x1D70F}}$ was applied to the slip lengths in addition to averaging in the homogeneous direction. Equation (3.5) is key to guaranteeing the correct wall stress, but we have the freedom to impose two more conditions to fully determine $l_{1}$ , $l_{2}$ and $l_{3}$ . For simplicity, we set $l_{1}=l_{2}=l_{3}=l$ , and compute $l$ from (3.5) with the correct $C_{f}$ prescribed. The value of $v_{2}$ is computed at each time step to ensure global zero mean mass flow through the wall such that

(3.6) $$\begin{eqnarray}v_{2}(t+\unicode[STIX]{x0394}t)=-\langle u_{2}(x_{1},0,x_{3},t)\rangle _{w},\end{eqnarray}$$

where $\unicode[STIX]{x0394}t$ is the time step and $\langle \cdot \rangle _{w}$ denotes average over the entire wall. The average $v_{2}$ obtained is of the order of $10^{-3}U_{\infty }$ . The slip velocities $v_{1}$ and $v_{3}$ in (2.1) are set to zero.

Figure 10(a) shows the resulting mean slip lengths computed to produce the target $C_{f}$ , which is successfully achieved as shown in figure 10(b). It is important to stress again that one of the main differences of the boundary layer case with respect to the channel flow is the necessity of a non-zero $v_{2}$ term from (2.1) to guarantee that the wall behaves as a no-transpiration boundary on average. We have implemented a global condition (constant-in-space $v_{2}$ , equation (3.6)) that does not prevent instantaneous local mass flow at a particular streamwise location as seen in figure 10(c). However, the mean mass flow remains locally close to zero for all streamwise locations.

Figure 10. (a) Mean slip lengths $l$ normalised by $\unicode[STIX]{x1D703}_{avg}$ , the average momentum thickness, (b) the friction coefficient from the LES with slip boundary condition (black ——) and the empirical friction coefficient from White & Corfield (Reference White and Corfield2006) (red ○), and (c) the instantaneous (black solid line) and time-averaged (red dashed line) wall-normal velocities at the wall as a function of  $Re_{\unicode[STIX]{x1D703}}$ .

Figure 11. (a) Mean streamwise velocity profile and (b) r.m.s. streamwise ( $\times$ ), spanwise (○) and wall-normal (▿) fluctuation profiles at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ . Symbols are for LES. DNS from Sillero et al. (Reference Sillero, Jiménez and Moser2013) (dashed line).

The mean streamwise velocity and the three r.m.s. velocity fluctuations at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ ( $Re_{\unicode[STIX]{x1D70F}}\approx 1989.5$ ) are shown in figure 11 and compared with Sillero et al. (Reference Sillero, Jiménez and Moser2013). As expected, the mean DNS and LES velocities match in the wake region, as the correct $C_{f}$ in the LES is imposed. The shape of the profile is also well predicted. The r.m.s. velocities are reasonably well reproduced at this Reynolds number, with no over-prediction of the streamwise r.m.s. velocity and under-prediction of the other two components close to the wall, consistent with the analysis in § 3.4. Overall, these results along with those from LES of channel flow in the previous sections show that the slip boundary condition successfully reproduces the one-point statistics of the flow as long as the slip lengths accurately reflect the correct mean wall stress.

Two more cases (not shown) were run to test the effect of the slip boundary condition on the net mass flow through the wall. In the first case, a slip boundary condition was imposed such that the net mass flow through the wall is positive (incoming flow through the wall) such that $\langle \bar{u}_{2}|_{w}\rangle _{w}\approx 0.01U_{\infty }$ . In this case, the boundary layer thickness grew five times faster than the reference DNS. In contrast, when the simulation was run with net negative mass flow through the walls ( $\langle \bar{u}_{2}|_{w}\rangle _{w}\approx -0.01U_{\infty }$ ), the flow remained laminar. The results are consistent with observations in previous studies on blowing and suction of boundary layers (Simpson, Moffat & Kays Reference Simpson, Moffat and Kays1969; Antonia et al. Reference Antonia, Fulachier, Krishnamoorthy, Benabid and Anselmet1988; Chung & Sung Reference Chung and Sung2001; Yoshioka & Alfredsson Reference Yoshioka and Alfredsson2006) and highlight the relevance of imposing a correct zero net mass flow through the walls to faithfully predict the boundary layer growth.

4.1 Previous dynamic models

Bose & Moin (Reference Bose and Moin2014) introduced a dynamic wall model based on the slip boundary condition free of any a priori parameters. The slip length, assumed to be equal for the three spatial directions, is computed via a modified form of the Germano’s identity (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991),

(4.1) $$\begin{eqnarray}l^{2}\left(\unicode[STIX]{x1D6E5}_{R}^{2}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{j}}{\unicode[STIX]{x2202}n}-\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\bar{u}_{j}}{\unicode[STIX]{x2202}n}\right)+\unicode[STIX]{x1D61B}_{ij}^{\mathit{SGS}}-\widehat{\unicode[STIX]{x1D70F}_{ij}}^{\mathit{SGS}}=\widehat{\bar{u}_{i}\bar{u}}_{j}-\bar{u}_{i}\bar{u}_{j},\end{eqnarray}$$

where $l$ is the slip length, $(\hat{\cdot })$ is the test filter, $\unicode[STIX]{x1D6E5}_{R}$ is the filter size ratio between the test and grid filters, $\unicode[STIX]{x1D70F}_{ij}^{\mathit{SGS}}$ and $\unicode[STIX]{x1D61B}_{ij}^{\mathit{SGS}}$ represent the grid and test filter SGS tensors, respectively. Equation (4.1) is then solved for $l$ by using least squares.

In Bose & Moin (Reference Bose and Moin2014), the model was tested for a series of LES of turbulent channel flow and NACA 4412 airfoil. However, our attempts to reproduce the channel flow results did not perform as expected with our current implementation, which uses a different SGS model and numerical discretisation. The discrepancies motivated a deeper study of the slip boundary condition and investigation of alternative dynamic wall models as that presented in the next section.

4.2 Wall-stress invariant dynamic wall model

We present a dynamic wall model formulated for the slip boundary condition based on the invariance of wall stress under test filtering. We will refer to this new model as the wall-stress invariant model (WSIM). We will assume a single slip length $l_{1}=l_{2}=l_{3}=l$ and neglect any potential contribution from the terms $v_{i}$ . The goal is to define a dynamic procedure to obtain $l$ at each time step such that the mean wall stress obtained from (3.1) is an accurate representation of the stress resulting from solving the Navier–Stokes equations with DNS resolution.

Let us consider the wall-stress operator ${\mathcal{T}}_{ij}$ given by

(4.2) $$\begin{eqnarray}{\mathcal{T}}_{ij}(\boldsymbol{u})=-\unicode[STIX]{x1D619}_{ij}(\boldsymbol{u})|_{w}-\unicode[STIX]{x1D70F}_{ij}^{\mathit{SGS}}(\boldsymbol{u})|_{w}+2\unicode[STIX]{x1D708}\unicode[STIX]{x1D61A}_{ij}(\boldsymbol{u})|_{w}-p(\boldsymbol{u})|_{w}\unicode[STIX]{x1D6FF}_{ij},\end{eqnarray}$$

where terms $\unicode[STIX]{x1D619}_{ij}$ , $\unicode[STIX]{x1D70F}_{ij}^{\mathit{SGS}}$ , $\unicode[STIX]{x1D61A}_{ij}$ and $p\unicode[STIX]{x1D6FF}_{ij}$ are the Reynolds stress, subgrid stress, strain-rate and pressure tensors, respectively, computed from the specified velocity field $\boldsymbol{u}$ . In our formulation of ${\mathcal{T}}_{ij}$ , the wall is assumed to be smooth, but the wall stress can be extended to encode information about the type of wall (smooth, rough, hydrophobic, etc.) by adding the appropriate drag terms into the right-hand side of (4.2).

The modelling choice for the dynamic wall model is

(4.3) $$\begin{eqnarray}\displaystyle & {\mathcal{T}}_{ij}(\bar{\boldsymbol{u}})-{\mathcal{T}}_{ij}(\hat{\bar{\boldsymbol{u}}})=0, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & {\mathcal{T}}_{ij}(\hat{\hat{\bar{\boldsymbol{u}}}})-\widehat{{\mathcal{T}}_{ij}}(\hat{\bar{\boldsymbol{u}}})=0, & \displaystyle\end{eqnarray}$$

where the different wall stresses, ${\mathcal{T}}_{ij}$ , are obtained by either computing the stress of the filtered velocity field or by filtering the total stress. Condition (4.3) enforces the invariance of wall stress under test filtering and allows the wall model to predict the same wall stress regardless of the grid resolution (or filter). A similar approach was also adopted by Anderson & Meneveau (Reference Anderson and Meneveau2011) for rough walls. Condition (4.4) is analogous to the Germano’s identity for the total stress of the filtered velocity and can be interpreted as a consistency condition between the wall stress and filter operator. The proposed model is given by combining (4.3) and (4.4) such that

(4.5) $$\begin{eqnarray}{\mathcal{F}}_{ij}={\mathcal{T}}_{ij}(\bar{\boldsymbol{u}})-{\mathcal{T}}_{ij}(\hat{\bar{\boldsymbol{u}}})+{\mathcal{T}}_{ij}(\hat{\hat{\bar{\boldsymbol{u}}}})-\widehat{{\mathcal{T}}}_{ij}(\hat{\bar{\boldsymbol{u}}})=0.\end{eqnarray}$$

The rationale for this particular combination of (4.3) and (4.4) is given by considering ${\mathcal{F}}_{12}$ , the dominant shear stress component in a boundary layer. This term can be simplified by test filtering only in the wall-normal direction with a box filter with filter size $\hat{\unicode[STIX]{x1D6E5}}+x_{2}$ at $x_{2}<\hat{\unicode[STIX]{x1D6E5}}$ . Assuming that the subgrid stress tensor is given by an eddy viscosity, i.e.  $(\unicode[STIX]{x1D70F}_{ij}^{\mathit{SGS}}-1/3\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D70F}_{kk}^{\mathit{SGS}})(\bar{\boldsymbol{u}})=-2\unicode[STIX]{x1D708}_{t}\unicode[STIX]{x1D61A}_{ij}(\bar{\boldsymbol{u}})$ and that $\unicode[STIX]{x1D708}_{t}$ is constant in the near-wall region, the first-order approximation of ${\mathcal{F}}_{12}$ can be shown to be

(4.6) $$\begin{eqnarray}{\mathcal{F}}_{12}=\frac{\hat{\unicode[STIX]{x1D6E5}}}{2}\left(\frac{\unicode[STIX]{x2202}\bar{u}_{1}\bar{u}_{2}}{\unicode[STIX]{x2202}x_{2}}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{2}}\left[(\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{t})|_{w}\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x_{2}}\right]\right)+\mathit{O}(\hat{\unicode[STIX]{x1D6E5}}^{2}).\end{eqnarray}$$

Hence, ${\mathcal{F}}_{12}=0$ implicitly enforces the well-established constant stress layer across the wall-normal direction. Note that (4.3) and (4.4) may be combined differently to produce alternative versions of (4.5), but not all these groupings lead to a first-order approximation consistent with the form reported in (4.6). For example, enforcing only (4.3) would not be consistent with (4.6). Nevertheless, wall models using different variants of (4.5) may be constructed based on similar principles; a broader family of dynamic models postulated on stress-invariant principles was investigated in Lozano-Durán et al. (Reference Lozano-Durán, Bae, Bose and Moin2017).

Next we introduce an explicit dependence on the slip length condition in (4.5). We will assume that the slip boundary condition also applies for the test-filtered velocity field,

(4.7) $$\begin{eqnarray}\hat{\bar{u}}_{i}=\hat{l}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{i}}{\unicode[STIX]{x2202}n},\end{eqnarray}$$

where $\hat{l}$ is the slip length at the test filter level. We further suppose that a linear functional dependence of the slip length with the filter size of the form $\hat{l}=\unicode[STIX]{x1D6E5}_{R}l$ holds. This assumption will be shown to be reasonably well satisfied in figure 13(c), which contains a posteriori values for the optimal slip length (3.2) as a function of grid resolution. Introducing the slip boundary condition at grid- and test-filter levels for the Reynolds stress terms in ${\mathcal{T}}_{ij}(\bar{\boldsymbol{u}})$ and ${\mathcal{T}}_{ij}(\hat{\bar{\boldsymbol{u}}})$ , such that

(4.8a,b ) $$\begin{eqnarray}\bar{u}_{i}\bar{u}_{j}=l^{2}\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\bar{u}_{j}}{\unicode[STIX]{x2202}n},\quad \hat{\bar{u}}_{i}\hat{\bar{u}}_{j}=l^{2}\unicode[STIX]{x1D6E5}_{R}^{2}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{j}}{\unicode[STIX]{x2202}n},\end{eqnarray}$$

the WSIM becomes

(4.9) $$\begin{eqnarray}l^{2}\left(\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\bar{u}_{j}}{\unicode[STIX]{x2202}n}-\unicode[STIX]{x1D6E5}_{R}^{2}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{j}}{\unicode[STIX]{x2202}n}\right)=\bar{u}_{i}\bar{u}_{j}-\hat{\bar{u}}_{i}\hat{\bar{u}}_{j}+{\mathcal{T}}_{ij}(\bar{\boldsymbol{u}})-{\mathcal{T}}_{ij}(\hat{\bar{\boldsymbol{u}}})+{\mathcal{T}}_{ij}(\hat{\hat{\bar{\boldsymbol{u}}}})-\widehat{{\mathcal{T}}}_{ij}(\hat{\bar{\boldsymbol{u}}}),\end{eqnarray}$$

which can be rewritten as

(4.10) $$\begin{eqnarray}l^{2}{\mathcal{M}}_{ij}={\mathcal{L}}_{ij}+{\mathcal{F}}_{ij},\end{eqnarray}$$

with

(4.11a,b ) $$\begin{eqnarray}{\mathcal{M}}_{ij}=\frac{\unicode[STIX]{x2202}\bar{u}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\bar{u}_{j}}{\unicode[STIX]{x2202}n}-\unicode[STIX]{x1D6E5}_{R}^{2}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{i}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\hat{\bar{u}}_{j}}{\unicode[STIX]{x2202}n},\quad {\mathcal{L}}_{ij}=\bar{u}_{i}\bar{u}_{j}-\hat{\bar{u}}_{i}\hat{\bar{u}}_{j}.\end{eqnarray}$$

The system (4.10) is over-determined and $l$ is computed via least squares as

(4.12) $$\begin{eqnarray}l^{2}=\frac{({\mathcal{L}}_{ij}+{\mathcal{F}}_{ij}){\mathcal{M}}_{ij}}{{\mathcal{M}}_{ij}{\mathcal{M}}_{ij}}=\frac{{\mathcal{L}}+{\mathcal{F}}}{{\mathcal{M}}},\end{eqnarray}$$

where repeated indices imply summation and the compact notation ${\mathcal{L}}={\mathcal{L}}_{ij}{\mathcal{M}}_{ij}$ , ${\mathcal{F}}={\mathcal{F}}_{ij}{\mathcal{M}}_{ij}$ and ${\mathcal{M}}={\mathcal{M}}_{ij}{\mathcal{M}}_{ij}$ is used. For incompressible flows, the isotropic part of $\unicode[STIX]{x1D70F}_{ij}^{\mathit{SGS}}$ is usually not defined by the SGS models and, since the system is already over-determined, we exclude the $i=j$ components of (4.12).

Note that the first part of the right-hand-side of (4.9), $\bar{u}_{i}\bar{u}_{j}-\hat{\bar{u}}_{i}\hat{\bar{u}}_{j}$ ( ${\mathcal{L}}_{ij}$ ), is the result of applying the boundary condition at the grid and test filter levels. The remaining terms, ${\mathcal{T}}_{ij}(\bar{\boldsymbol{u}})-{\mathcal{T}}_{ij}(\hat{\bar{\boldsymbol{u}}})+{\mathcal{T}}_{ij}(\hat{\hat{\bar{\boldsymbol{u}}}})-\widehat{{\mathcal{T}}}_{ij}(\hat{\bar{\boldsymbol{u}}})$ ( ${\mathcal{F}}_{ij}$ ), then act as an effective control such that the slip length increases if the current wall stress is under-predicted, and decreases if the wall stress is over-predicted. This self-regulating mechanism can be examined by analysing the terms ${\mathcal{M}}$ , ${\mathcal{L}}$ and ${\mathcal{F}}$ for three test cases of a channel flow at $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ with DSM and grid $\unicode[STIX]{x1D6E5}_{1}=\unicode[STIX]{x1D6E5}_{2}=\unicode[STIX]{x1D6E5}_{3}=0.05\unicode[STIX]{x1D6FF}$ . The first case is computed by imposing the optimal slip length, $l=l_{opt}=0.009\unicode[STIX]{x1D6FF}$ . The second and third cases are analogous to the first case but with $l=1.70l_{opt}$ and $l=0.35l_{opt}$ , respectively. The terms ${\mathcal{M}}$ , ${\mathcal{L}}$ and ${\mathcal{F}}$ were evaluated after the cases were run with their corresponding slip lengths fixed in time until the statistically steady state was reached. Note that ${\mathcal{L}}$ can be interpreted as the model prior to applying the control mechanism, and this allows us to define two slip lengths, namely, $l_{{\mathcal{L}}}={\mathcal{L}}/{\mathcal{M}}$ and $l_{{\mathcal{L}}+{\mathcal{F}}}=({\mathcal{L}}+{\mathcal{F}})/{\mathcal{M}}$ .

The terms ${\mathcal{M}}$ , ${\mathcal{L}}$ and ${\mathcal{L}}+{\mathcal{F}}$ evaluated from the three test cases are plotted in figure 12(a), and the corresponding slip lengths $l_{{\mathcal{L}}}$ and $l_{{\mathcal{L}}+{\mathcal{F}}}$ are presented in figure 12(b). The results show that application of WSIM recovers the optimal slip length, and thus the correct wall stress, through the control mechanism ${\mathcal{F}}_{ij}$ . The analysis provided here is performed a priori, that is, the wall model was used to predict $l$ at time $t+\unicode[STIX]{x0394}t$ for a given flow field at time $t$ , but without an actual dynamic coupling between WSIM and LES. The remainder of the paper is devoted to testing the WSIM in LES under different test scenarios.

Figure 12. (a) Plots of ${\mathcal{L}}$ (●), ${\mathcal{L}}+{\mathcal{F}}$ (○) and ${\mathcal{M}}$ (▿) computed from LES of channel flow using the slip boundary condition with fixed $l$ equal to $l=0.35l_{opt}=0.003\unicode[STIX]{x1D6FF}$ , $l=l_{opt}=0.009\unicode[STIX]{x1D6FF}$ , and $l=1.70l_{opt}=0.015\unicode[STIX]{x1D6FF}$ . (b) The slip lengths $l_{{\mathcal{L}}}$ (●), $l_{{\mathcal{L}}+{\mathcal{F}}}$ (○) normalised by the optimal slip length. Here $\unicode[STIX]{x1D6E5}_{R}$ was assigned to be $1.6$ (see § 5.1). The vertical dotted lines are $l=l_{opt}$ . Red arrows highlight the improvement achieved by including the control term ${\mathcal{F}}_{ij}$ . See the text for further details.

Figure 13. Error in the streamwise mean velocity profile, ${\mathcal{E}}$ , as a function of (a) grid size (for $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ ) and (b) Reynolds number (for grid G1). WSIM (○) and EQWM (♢). The predicted slip lengths $l/\unicode[STIX]{x1D6FF}$ for WSIM (solid lines) and optimal slip lengths (dashed lines) as a function of (c) grid resolution for $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ and (d) Reynolds number for grid G1.

5.1 Numerical experiments

To test the performance of the WSIM, three flow configurations are considered: a statistically steady plane turbulent channel (2-D channel), a non-equilibrium three-dimensional transient channel (3-D channel) and a zero-pressure-gradient flat-plate turbulent boundary layer. The numerical methods of the simulations are the same as those given in §§ 3.1 and 3.6. The 2-D channel flow and the turbulent boundary layer were discussed in § 3. The three-dimensional transient channel flow is a temporally developing turbulent boundary layer in a planar channel subjected to a sudden spanwise forcing as in Moin et al. (Reference Moin, Shih, Driver and Mansour1990).

The size of the 2-D and 3-D channel domain is $8\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}\times 2\unicode[STIX]{x1D6FF}\times 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ in the streamwise, wall-normal and spanwise directions, respectively. For both 2-D and 3-D channel flows, periodic boundary conditions are applied in the streamwise and spanwise directions. For the top and bottom walls, we impose either the no-slip (NS), slip boundary condition for WSIM, or Neumann boundary condition for cases with the equilibrium wall model (EQWM). The formulation for the EQWM follows Kawai & Larsson (Reference Kawai and Larsson2013) with a matching location at the third grid cell for the streamwise velocity, although recent studies have shown that the first grid cell may be used for the EQWM when the velocities are filtered using a spatial or temporal filter (Yang, Park & Moin Reference Yang, Park and Moin2017) following the methodology first introduced for algebraic wall models (Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2004).

For the 2-D channel, the flow is driven by imposing a constant mean pressure gradient and the simulations are started from a random initial condition run for at least $100\unicode[STIX]{x1D6FF}/u_{\unicode[STIX]{x1D70F}}$ after transients. In the case of the 3-D channel, the calculations were started from a 2-D fully developed plane channel flow driven by a streamwise mean pressure gradient. The subsequent calculations were performed with a transverse (spanwise) mean pressure gradient of $\unicode[STIX]{x2202}P/\unicode[STIX]{x2202}x_{3}=10\unicode[STIX]{x1D70F}_{w0}/\unicode[STIX]{x1D6FF}$ , where $\unicode[STIX]{x1D70F}_{w0}$ is the mean wall shear stress of the unperturbed channel. The 3-D channel simulations were run for $10u_{\unicode[STIX]{x1D70F}0}/\unicode[STIX]{x1D6FF}$ and averaged over seven realisations, where $u_{\unicode[STIX]{x1D70F}0}$ is the friction velocity of the 2-D initial condition.

For the boundary layer, the set-up is identical to that in § 3.6. The range of $Re_{\unicode[STIX]{x1D703}}$ is from 1000 to 10 000. The length, height and width of the simulated box are $L_{x}=1060\unicode[STIX]{x1D703}_{avg}$ , $L_{y}=18\unicode[STIX]{x1D703}_{avg}$ and $L_{z}=35\unicode[STIX]{x1D703}_{avg}$ with the streamwise and spanwise resolutions of $\unicode[STIX]{x1D6E5}_{1}/\unicode[STIX]{x1D6FF}=0.05$ ( $\unicode[STIX]{x1D6E5}_{1}^{+}=118$ ) and $\unicode[STIX]{x1D6E5}_{3}/\unicode[STIX]{x1D6FF}=0.04$ ( $\unicode[STIX]{x1D6E5}_{3}^{+}=84.3$ ) at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ . The grid is slightly stretched in the wall-normal direction with minimum $\unicode[STIX]{x1D6E5}_{2}/\unicode[STIX]{x1D6FF}=0.01$ ( $\unicode[STIX]{x1D6E5}_{2}^{+}=20.8$ ). The inlet, outlet and top boundary conditions are as in § 3.6 with $x_{ref}/\unicode[STIX]{x1D703}_{0}=890$ . For the wall, we impose the slip boundary condition for the WSIM.

It is important to note the details of the filter operation, as dynamic wall models are particularly sensitive to this choice (see Appendix). Test filtering a variable $f$ in a given spatial direction at point $i$ is computed as $1/6f(i-1)+2/3f(i)+1/6f(i+1)$ (Simpson’s rule). The operation is repeated for all three directions away from the wall. This corresponds to a discrete fourth-order quadrature over a cell of size $2\unicode[STIX]{x1D6E5}_{1}\times 2\unicode[STIX]{x1D6E5}_{2}\times 2\unicode[STIX]{x1D6E5}_{3}$ for a uniform grid. At the wall, the same filtering operation is used in the horizontal directions while the wall-normal filter is one-sided and given by $2/3f(1)+1/3f(2)$ , with $f(1)$ and $f(2)$ denoting values at the first and second wall-normal grid points. This is an integration over a cell of size $2\unicode[STIX]{x1D6E5}_{1}\times \unicode[STIX]{x1D6E5}_{2}\times 2\unicode[STIX]{x1D6E5}_{3}$ . In addition, the definition of the filter operation fixes the value of $\unicode[STIX]{x1D6E5}_{R}$ , which is the ratio between the grid and test filter sizes at the wall. In this case, the $\unicode[STIX]{x1D6E5}_{R}$ based on the cell volume is given by $\sqrt[3]{2\times 1\times 2}\approx 1.6$ .

The cases for the 2-D and 3-D channel are labelled following the convention ([Channel type]-[Wall model]-[Reynolds number]-[Grid]), where the grid labels G0, G1, and G2 given in table 2 correspond to $320\times 25\times 120$ ( $\unicode[STIX]{x1D6E5}_{1}=\unicode[STIX]{x1D6E5}_{2}=\unicode[STIX]{x1D6E5}_{3}=0.080\unicode[STIX]{x1D6FF}$ ), $512\times 40\times 192$ ( $\unicode[STIX]{x1D6E5}_{1}=\unicode[STIX]{x1D6E5}_{2}=\unicode[STIX]{x1D6E5}_{3}=0.050\unicode[STIX]{x1D6FF}$ ), and $1024\times 80\times 384$ ( $\unicode[STIX]{x1D6E5}_{1}=\unicode[STIX]{x1D6E5}_{2}=\unicode[STIX]{x1D6E5}_{3}=0.025\unicode[STIX]{x1D6FF}$ ), respectively. The wall models applied are labelled NS, EQWM or WSIM. Additional cases with anisotropic grids, different values of $\unicode[STIX]{x1D6E5}_{R}$ , test-filtering operations or SGS models were run to study the sensitivity of the model to these choices. They are discussed in the Appendix but not included in the table.

Table 2. Grid resolutions in outer units. The first column contains the label used to name LES cases for the 2-D and 3-D channel flow simulations computed with different grids. The second, third and fourth columns are the streamwise, wall-normal and spanwise grid resolutions, respectively.

The 2-D channel results are compared with DNS data from Hoyas & Jiménez (Reference Hoyas and Jiménez2006) and Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014a ) for $Re_{\unicode[STIX]{x1D70F}}\approx 2000$ and 4200, Yamamoto & Tsuji (Reference Yamamoto and Tsuji2018) for $Re_{\unicode[STIX]{x1D70F}}\approx 8000$ and with the law-of-the wall for $Re_{\unicode[STIX]{x1D70F}}>8000$ . For the boundary layer, the resulting friction coefficient is compared with the empirical $C_{f}$ from White & Corfield (Reference White and Corfield2006), and the mean velocity profiles are compared with the DNS data from Sillero et al. (Reference Sillero, Jiménez and Moser2013) at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ and the experimental data from Österlund (Reference Österlund1999) at $Re_{\unicode[STIX]{x1D703}}\approx 8000$ .

The performance of the WSIM in laminar flows has not been studied. However, in the limit of fine grids, the dynamic procedure of WSIM should produce zero slip lengths and revert to the no-slip boundary condition. Thus, for laminar cases with enough grid resolution to resolve the near-wall structures, we expect the dynamic wall model to naturally switch off.

5.2 Statistically steady two-dimensional channel flow

We assess the performance of the WSIM compared with EQWM and NS. The results are discussed in terms of the error in the streamwise mean velocity profile across the logarithmic region. This choice was necessary to include higher-Reynolds-number cases where the corresponding DNS was not available and the law of the wall is used instead. Restricting the error to be evaluated only in the logarithmic layer is justified as wall models mainly impact the solution by vertically shifting the mean velocity profile and do not alter its shape for the range of grid resolutions tested as shown in § 3.3. In particular, the error is measured as the normalised $L_{2}$ error of the streamwise mean velocity between the second grid point and $0.2\unicode[STIX]{x1D6FF}$ ,

(5.1) $$\begin{eqnarray}{\mathcal{E}}=\left[\frac{\displaystyle \int _{\unicode[STIX]{x1D6E5}_{2}}^{0.2\unicode[STIX]{x1D6FF}}(\langle \bar{u}_{1}\rangle -\langle u_{1}^{\mathit{DNS}}\rangle )^{2}\,\text{d}x_{2}}{\displaystyle \int _{\unicode[STIX]{x1D6E5}_{2}}^{0.2\unicode[STIX]{x1D6FF}}(\langle u_{1}^{\mathit{DNS}}\rangle )^{2}\,\text{d}x_{2}}\right]^{1/2}.\end{eqnarray}$$

In the case where the corresponding DNS does not exist, $\langle u_{1}^{\mathit{DNS}}\rangle$ is replaced by the law of the wall,

(5.2) $$\begin{eqnarray}\langle u_{1}^{\mathit{DNS}+}\rangle =\frac{1}{\unicode[STIX]{x1D705}}\log x_{2}^{+}+B,\end{eqnarray}$$

with $\unicode[STIX]{x1D705}=0.392$ and $B=4.48$ (Luchini Reference Luchini2017).

Figure 13(a,b) show ${\mathcal{E}}$ as a function of grid resolution and Reynolds number. At moderate Reynolds numbers ( $Re_{\unicode[STIX]{x1D70F}}<8000$ ) and all grid resolutions, the error for WSIM ( ${\mathcal{E}}\approx 2\,\%{-}6\,\%$ ) is similar to that of the EQWM ( ${\mathcal{E}}\approx 2\,\%{-}3\,\%$ ). With increasing Reynolds number, the performance degrades (up to ${\mathcal{E}}\approx 15\,\%$ at $Re_{\unicode[STIX]{x1D70F}}\approx 20\,000$ ), while the EQWM does not. The accurate results for EQWM are not surprising as its modelling assumptions are well satisfied for channel flow settings. The reason for the declining performance of WSIM at very high Reynolds number can be found in §§ 2.1 and 2.2, where it was argued that the underlying assumptions for the slip condition are invalidated for large filter sizes. However, it is worth mentioning that the errors for an LES with no wall model ( ${\mathcal{E}}\approx 100\,\%$ for 2D-NS-4200-G1) are an order of magnitude larger than the errors of WSIM for all cases. The mean velocity profiles and streamwise r.m.s. velocity fluctuations for WSIM for various cases are shown in figure 14. Additional sensitivities to anisotropic grids, different values of $\unicode[STIX]{x1D6E5}_{R}$ , test-filtering operations or SGS model are discussed in the Appendix.

The slip lengths predicted by the WSIM are shown in figure 13(c,d) as a function of grid resolution and Reynolds number and compared to the optimal slip lengths (3.2). It is remarkable that the WSIM captures the overall behaviour of the optimal slip lengths, that is, a strong dependence on grid resolution and a weak variation with Reynolds number.

Figure 14. (a) Mean velocity profiles and (b) streamwise r.m.s. velocity fluctuations for WSIM at $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ for grid G0 (▿), G1 (○) and G2 (▫). DNS for $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ (dashed line). (c) Mean velocity profiles and (d) streamwise r.m.s. velocity fluctuations for WSIM at $Re_{\unicode[STIX]{x1D70F}}\approx 2000$ (▫), 4200 (○), 8000 (▿) and 20 000 (♢) for grid G1. DNS for $Re_{\unicode[STIX]{x1D70F}}\approx 4200$ (dashed line).

5.3 Three-dimensional transient channel flow

The performance of WSIM in non-equilibrium scenarios is assessed in a three-dimensional transient channel flow (Moin et al. Reference Moin, Shih, Driver and Mansour1990). Note that in general wall models cannot be assumed to be effective at transferring information from the inner to the outer layer in non-equilibrium flows. Hence, the current flow set up, characterised by a spanwise boundary layer growing from the wall owing to viscous effects, is expected to be problematic for wall-modelled LES. A plane channel flow was modified to incorporate a lateral (transverse) pressure gradient 10 times that of the streamwise pressure gradient. The details of the simulations were given in § 5.1.

The wall models explored are the WSIM and EQWM. A case with the no-slip boundary condition is used for control, and the figure of merit is the evolution of the streamwise and spanwise wall stress as a function of time (figure 15). Note that the temporal increment of the wall stress magnitude involves an increase of the Reynolds number from $Re_{\unicode[STIX]{x1D70F}}\approx 932$ at $t=0$ to $Re_{\unicode[STIX]{x1D70F}}\approx 2600$ at $t=10\unicode[STIX]{x1D6FF}/u_{\unicode[STIX]{x1D70F}0}$ . The results show that the performance of the WSIM is similar to the EQWM despite its parameter-free nature. The streamwise and spanwise mean velocity profiles at various time instances are given in figure 16, which also shows that both the WSIM and EQWM predict similar time evolutions. Although there is no reference DNS available for the full time span of our simulations, in the limited time range from $t=0$ to $1\unicode[STIX]{x1D6FF}/u_{\unicode[STIX]{x1D70F}0}$ , Giometto et al. (Reference Giometto, Lozano-Durán, Park and Moin2017) showed that the EQWM predicts the evolution of the wall stresses with less than 10 % deviation in the spanwise wall-stress prediction from DNS, and thus the results from the WSIM are also expected to exhibit a similar error. Both the EQWM and WSIM entail a quantitative correction to the prediction provide by the no-slip boundary condition. Consequently, the computational simplicity and absence of a secondary mesh makes the WSIM an appealing approach at the cost of a moderate attenuation of the predictive capabilities compared with more sophisticated wall models.

Figure 15. Wall stress in (a) streamwise and (b) spanwise directions as a function of time for the WSIM (○), EQWM (♢) and NS ( $\times$ ).

Figure 16. Mean (a) streamwise and (b) spanwise velocity profiles as a function of $x_{2}/\unicode[STIX]{x1D6FF}$ at $tu_{\unicode[STIX]{x1D70F}0}/\unicode[STIX]{x1D6FF}=0$ (solid lines), 4.5 (dashed lines) and 9 (dotted lines) for the WSIM (○) and EQWM (♢).

5.4 Zero-pressure-gradient flat-plate turbulent boundary layer

Finally, the performance of WSIM is assessed in a flat-plate turbulent boundary layer. The details of the numerical set-up were discussed in § 5.1. The friction coefficient is shown in figure 17 from $Re_{\unicode[STIX]{x1D703}}\approx 1000$ to 10 000. Note that the recycling scheme of Lund et al. (Reference Lund, Wu and Squires1998) imposes an artificial boundary condition at the inlet, requiring an initial development region for the flow to fully adapt to the slip boundary condition, which is the reason for the discrepancy in $C_{f}$ near the inlet. Consistent with previous test cases, the results show that WSIM predicts the friction coefficient well within 4 % error for $Re_{\unicode[STIX]{x1D703}}>6000$ . The mean streamwise velocity profile at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ and 8000 and the r.m.s. velocity fluctuations at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ are also well predicted as reported in figure 18.

Figure 17. Empirical friction coefficient from White & Corfield (Reference White and Corfield2006) (red ○) and the friction coefficient from WSIM (black ——).

Figure 18. Mean streamwise velocity profile for (a) $Re_{\unicode[STIX]{x1D703}}\approx 6500$ and (b) $Re_{\unicode[STIX]{x1D703}}\approx 8000$ , and the r.m.s. (c) streamwise ( $\times$ ), (d) spanwise (○) and wall-normal (▿) fluctuation profiles at $Re_{\unicode[STIX]{x1D703}}\approx 6500$ . WSIM (symbols) and DNS (Sillero et al. Reference Sillero, Jiménez and Moser2013) or experiment (Österlund Reference Österlund1999) (dashed line).