4.1. Determination of $C_{k}$

There are a number of ways in which the relation between the filter width and the local integral length scale can be specified. The model coefficient $C_{k}$ can, for example, be specified by the user based on the computational resources available. Alternatively, one could find an optimum value for $C_{k}$ that minimizes the total simulation error in the prediction of a set of selected quantities that may be known from experiments, theoretical considerations or DNS. For instance, Piomelli & Geurts (Reference Piomelli and Geurts2010, Reference Piomelli, Geurts, Kuerten, Geurts, Armenio and Fröhlich2011) required that the coefficient minimizes the error in the prediction of the skin-friction coefficient, $C_{f}=2{\it\tau}_{w}/{\it\rho}U_{b}^{2}$ (where ${\it\tau}_{w}$ is the wall stress, ${\it\rho}$ the fluid density and $U_{b}$ the average velocity in the channel), which is known from experimental data.

An alternative is to use an internal measure of the simulation resolution to determine $C_{k}$ . For this study, therefore, we define two measures of the contribution to the transport by the SFS. The first is the ratio of SFS dissipation to the total dissipation, defined as (Meyers et al. Reference Meyers, Geurts and Baelmans2003)

(4.1) $$\begin{eqnarray}s_{{\it\varepsilon}}=\frac{\langle {\it\varepsilon}_{sfs}\rangle }{\langle {\it\varepsilon}_{sfs}\rangle +\langle {\it\varepsilon}_{{\it\nu}}\rangle }=\frac{\langle 2{\it\nu}_{sfs}\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle }{\langle 2({\it\nu}+{\it\nu}_{sfs})\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle }.\end{eqnarray}$$

Geurts & Fröhlich (Reference Geurts and Fröhlich2002) and Meyers et al. (Reference Meyers, Geurts and Baelmans2003) observed that for large filter-width-to-grid ratios, $s_{{\it\varepsilon}}$ depends mostly on ${\it\Delta}$ . An appropriate choice of the averaging operator might be case-dependent; here, we have averaged over the entire domain and time.

A second measure is the SFS contribution to the Reynolds stresses. This measure should, however, take into account the fact that the SFS stress model, in the eddy-viscosity formulation, does not include the SFS energy. Therefore, the simplest coordinate-invariant measure of the Reynolds stresses is the second invariant of the anisotropic part of the Reynolds stress tensor. Let

(4.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D619}_{ij}^{a}=\overline{u}_{i}^{\prime }\overline{u}_{j}^{\prime }-{\it\delta}_{ij}\frac{\overline{u}_{k}^{\prime }\overline{u}_{k}^{\prime }}{3};\quad {\it\tau}_{ij}^{a}={\it\tau}_{ij}-{\it\delta}_{ij}\frac{{\it\tau}_{kk}}{3},\end{eqnarray}$$

where $\overline{u}_{i}^{\prime }=\overline{u}_{i}-\langle \overline{u}_{i}\rangle _{hd,t}$ is the resolved turbulent fluctuation, difference between resolved and mean velocity (averaged over homogeneous directions and time). We define

(4.3) $$\begin{eqnarray}s_{{\it\tau}}=\left[\frac{\langle {\it\tau}_{ij}^{a}{\it\tau}_{ij}^{a}\rangle }{\langle (\unicode[STIX]{x1D619}_{lm}^{a}+{\it\tau}_{lm}^{a})(\unicode[STIX]{x1D619}_{lm}^{a}+{\it\tau}_{lm}^{a})\rangle }\right]^{1/2}.\end{eqnarray}$$

Once a desired level of $s_{{\it\varepsilon}}$ or $s_{{\it\tau}}$ is assigned, the value of $C_{k}$ that ensures that the SFS model yields the desired subfilter activity can be found by performing a series of calculations on coarse grids; production calculations can subsequently be performed with higher reliability on finer meshes. With this procedure, $C_{k}$ becomes a user-defined resolution parameter, related to the cost of the calculation, and to the computational power available. Large values of $s_{{\it\varepsilon}}$ or $s_{{\it\tau}}$ result in large values of $C_{k}$ , corresponding to fairly coarse LES in which most of the energy transfer is modelled, and/or a large proportion of the Reynolds stresses are unresolved. Such calculations can be expected to converge on coarse grids. Small values of $s_{{\it\varepsilon}}$ or $s_{{\it\tau}}$ , conversely, will result in small values of $C_{k}$ , and finer meshes will be required. In the limit $s\rightarrow 0$ , the filter width vanishes and the DNS limit is recovered as the mesh is refined.

To minimize the number of coarse-grid calculations required to determine $C_{k}$ , the error-minimization technique based on the successive inverse polynomial interpolation (SIPI) procedure, as proposed by Geurts & Meyers (Reference Geurts and Meyers2006) to optimize the Smagorinsky constant, is used. The SIPI method is based on the definition of a measure ${\it\phi}$ of the difference between some target quantity and the LES prediction of this quantity. Geurts & Meyers (Reference Geurts and Meyers2006), for instance, performed simulations of HIT, and based the definition of ${\it\phi}$ on the TKE. The method consists of performing three LES on a coarse mesh, for values of $C_{k}$ that bracket the optimum (one evaluation is always performed with $C_{k}=0$ , one with a fairly large value of $C_{k}$ , the third with an intermediate value). We then construct an interpolating parabola, and find its minimum, which is used for the next evaluation of ${\it\phi}$ . The interpolation is then repeated using three pairs $[C_{k},{\it\phi}(C_{k})]$ (including the newest one) that result in a parabola with an upwards concavity. A new minimum is found and a new error evaluation carried out. This procedure resulted in an accurate identification of the minimum with typically less than six error evaluations. The cost of the whole procedure is small, given that the SIPI test runs are performed on a mesh that is much coarser than that used for the production runs. The use of the coarse mesh is made possible by the fact that the integral length scale (3.3) (and, hence, $\ell$ and ${\it\Delta}$ ) is representative of the largest turbulence scales, and is expected to be predicted reasonably well even by coarse calculations. Obviously, the latter requirement does implicitly set some lower bound on the coarsest resolutions that actually can be adopted in such a set of SIPI test runs. In the applications considered here, the overhead created by these SIPI runs was found to be less than 5 %.

Figure 3 demonstrates the application of SIPI on a channel flow using the proposed model at $\mathit{Re}_{{\it\tau}}=950$ on a coarse grid ( $48\times 65\times 48$ ). The parameters of the simulation are discussed later. For this case ${\it\phi}$ is based on subfilter activity $s_{{\it\varepsilon}}$ :

(4.4) $$\begin{eqnarray}{\it\phi}_{s}\hspace{2.22198pt}\%=100\times \left|s_{{\it\varepsilon}}-s_{tgt}\right|^{2}\end{eqnarray}$$

where $s_{tgt}$ is set to 0.23; this value would be typical of a well-resolved LES of a wall-bounded flow; the sensitivity of the results to the target value will be examined later. The first three calculations were performed with $C_{k}=0.000$ (point (i)), $C_{k}=0.009$ (point (ii)) and $C_{k}=0.0045$ (point (iii)). The parabola constructed from these three points yields a minimum for $C_{k}=0.007$ . A fourth calculation was carried out with this value, and the next parabola created by points (ii), (iii) and (iv); it also has a minimum at $C_{k}=0.007$ , indicating this value as the target value. In total four simulations were performed to optimize $C_{k}$ . The optimum value obtained directly by sampling ${\it\phi}_{s}$ in a number of cases was also verified by performing nine simulations (figure 3 b), and is the same as that found using SIPI.

Figure 3. SIPI optimization for plane channel flow at $\mathit{Re}_{{\it\tau}}=950$ with subfilter activity $s_{{\it\varepsilon}}$ as quantity to be optimized with target value $s_{tgt}=0.23$ and grid $48\times 65\times 48$ . (a) Error optimization based on SIPI; the first polynomial has been constructed using points (i)–(iii) (——); the second using (ii)–(iv) (- - - -). (b) Error obtained from nine simulations with varying values of $C_{k}$ .

As mentioned, all of the averages required to evaluate the SFS ratios, $s_{{\it\varepsilon}}$ and $s_{{\it\tau}}$ , are carried out over the entire computational domain and time. A single value of the coefficient $C_{k}$ obtains; the model length scale $\ell$ , however, may vary if the TKE and dissipation are only averaged in homogenous directions (as will be the case in the plane channel flow). This choice is by no means unique, and is made here to simplify the initial applications of the model. Its effects and limitations will be discussed in § 4.3. More local averages can also be used, and are the subject of ongoing work.

4.2. Forced HIT

We first performed simulations of HIT, with Reynolds number (based on $\mathscr{K}^{1/2}$ and the integral scale $L$ ) $\mathit{Re}_{L}=250$ . To achieve a steady state, a force was added to the right-hand side of the momentum equation. Following Rosales & Meneveau (Reference Rosales and Meneveau2005) $A=0.1333$ in (2.3). The domain size was $2{\rm\pi}\times 2{\rm\pi}\times 2{\rm\pi}$ , and the initial condition was a solenoidal isotropic velocity field with random phases and the energy spectrum:

(4.5) $$\begin{eqnarray}E_{T}(k)=\frac{16}{\sqrt{{\rm\pi}/2}}\frac{u_{o}^{2}k^{4}}{k_{o}^{5}}\exp \left(-\frac{2k^{2}}{k_{o}^{2}}\right),\end{eqnarray}$$

where $u_{o}=1$ is the initial root-mean-square (r.m.s.) velocity, and $k_{o}=20$ the wavenumber where the spectrum is maximum. These input parameters correspond to one of the three DNS cases performed by Rosales & Meneveau (Reference Rosales and Meneveau2005) using a pseudo-spectral method. Note that, with the forcing defined in (2.3), at steady state $A=\langle \overline{{\it\varepsilon}}_{tot}\rangle _{{\it\Omega},t}/2\langle \mathscr{\overline{K}}\rangle _{{\it\Omega},t}$ , yielding $L=\langle \mathscr{\overline{K}}\rangle _{{\it\Omega},t}^{1/2}/2A$ . The DNS simulation by Rosales & Meneveau (Reference Rosales and Meneveau2005) showed that the TKE reaches an asymptotic state around 0.3; therefore, we expect $L\simeq 2.05$ .

Here $C_{k}$ was determined by choosing a target value of the subfilter activity $s_{{\it\varepsilon},tgt}=0.15$ . Simulations were then performed on a coarse grid with $32^{3}$ points, and the optimum value of $C_{k}$ was found using SIPI; only 5 simulations were required to find that $C_{k}\simeq 0.006$ gave $s_{{\it\varepsilon},tgt}=0.15$ . To evaluate the sensitivity of the results to increasing the range of eddies that are resolved, we also performed the optimization using $s_{{\it\varepsilon},tgt}=0.3$ , a value expected to give a more significant contribution of the SFS to the momentum transport. The coefficient increased to 0.010.

Figure 4(a,b) show the energy spectra $E_{T}(k)$ calculated, for $C_{k}=0.006$ and $C_{k}=0.010$ , using different grids. Several observations can be made; in both figures the energy spectra collapse onto each other for $128^{3}$ grid points when $C_{k}=0.006$ , and $32^{3}$ when $C_{k}=0.010$ . This is consistent with the observations by Geurts & Fröhlich (Reference Geurts and Fröhlich2002) and Meyers et al. (Reference Meyers, Geurts and Baelmans2003) in which, for higher subfilter activity, grid convergence was achieved at coarser resolution. The value $L\simeq 2.05$ is also achieved to good approximation using both $C_{k}$ values when $64^{3}$ grid points are used (not shown). Finally, we observe that when $C_{k}=0.006$ , smaller eddies (higher wavenumbers) are resolved while the use of a higher $C_{k}$ value tends to increase the contribution of the subfilter model, resulting in a strong attenuation of the predicted spectra at high wavenumbers (figure 4 b). For the higher value of $C_{k}$ , the spectrum decays more rapidly at high wavenumbers. Since the DNS data are not filtered, this results in a larger difference with the reference data.

Figure 4. Total energy spectra at grid resolutions of $---$   $32^{3}$ , — ⋅ —  $64^{3}$ , - - - -  $128^{3}$ and ——  $192^{3}$ with: (a) $s_{{\it\varepsilon},tgt}=0.15$ and $C_{k}=0.006$ ; (b) $s_{{\it\varepsilon},tgt}=0.3$ and $C_{k}=0.010$ . The ● shows DNS (Rosales & Meneveau Reference Rosales and Meneveau2005).

In figure 5 we compare the flow structures obtained, for different grid resolutions, for the $C_{k}=0.010$ case. The eddies are visualized using isosurfaces of the second invariant of the velocity gradient tensor, $\unicode[STIX]{x1D64C}$ (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988; Dubief & Delcayre Reference Dubief and Delcayre2000):

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D64C}=-\frac{1}{2}\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{j}}{\partial x_{i}}.\end{eqnarray}$$

We observe progressively smaller structures being resolved as the mesh is refined, and also more finely detailed appearance of the resolved eddies. At coarse grids the contribution of numerical discretization errors is quite significant. Conversely, beyond a resolution of $128^{3}$ the global visual impression of the flow appears quite similar, suggesting the start of convergence to a grid-independent LES. For the lower value of $C_{k}$ (not shown) smaller scales appear, as expected.

Figure 5. Isosurfaces of $\unicode[STIX]{x1D64C}=5.00$ for $C_{k}=0.010$ : (a) $32^{3}$ points; (b) $64^{3}$ points; (c)  $128^{3}$ points; (d) $192^{3}$ points.

We also performed simulations at a higher Reynolds number, $\mathit{Re}_{L}=4490$ . The spectra are shown in figure 6, for $s_{{\it\varepsilon}}=0.84$ (resulting in $C_{k}=0.006$ ) and $s_{{\it\varepsilon}}=0.44$ (which gives $C_{k}=0.0015$ ). Grid-convergence appears to set in beyond $128^{3}$ and $256^{3}$ grid points, respectively. The correspondence between the spectrum obtained here, with an undefined filter function, and the model spectrum coupled with a sharp Fourier cutoff discussed before is only approximate; however, we can remark, with reference to figure 2(b), that, at this Reynolds number, $s_{{\it\varepsilon}}=0.84$ corresponds to $C_{{\it\Delta}}\simeq 0.05$ ; the filter width is approximately 20 times smaller than the length scale $L_{est}$ . The corresponding wavenumber is approximately $k_{{\it\Delta}}\simeq 31$ (approximately 77 for $s_{{\it\varepsilon}}=0.44$ ), and is indicated by a dashed line in the figure. The LES results match the DNS data well up to wavelengths $\simeq 2{\it\Delta}$ . Near the cutoff the spectra decay considerably more rapidly, reflecting the high value of $s_{{\it\varepsilon}}$ , which causes the model to be very dissipative. Compared with a Smagorinsky model, in which the eddy viscosity is given by

(4.7) $$\begin{eqnarray}{\it\nu}_{smag}=(C_{s}h)^{2}|\overline{S}|,\end{eqnarray}$$

with $C_{s}\simeq 0.17{-}0.23$ , the present model gives higher dissipation by a factor $(C_{k}L_{est}/C_{s}h)^{2}$ . For the grid-converged simulation, with $h=2{\rm\pi}/192\simeq 0.016L_{est}$ , this implies that the dissipation predicted by the present model is over three times larger than that of the Smagorinsky model. Choosing a smaller value of the SFS activity, $s_{{\it\varepsilon}}=0.44$ results in a wider range of scales in agreement with the DNS data; however, a finer mesh ( $256^{3}$ grid points) is required to obtain a similar level of grid-converged data; again, eddies of size ${>}2{\it\Delta}$ appear accurately captured.

Figure 6. Total energy spectra at grid resolutions $---$   $64^{3}$ , — ⋅ —  $128^{3}$ , - - - -  $192^{3}$  and ——  $256^{3}$ with $\mathit{Re}_{L}=4490$ and: (a)  $s_{{\it\varepsilon},tgt}=0.84$ and $C_{k}=0.006$ . (b)  $s_{{\it\varepsilon},tgt}=0.44$ and $C_{k}=0.0015$ . The ● shows DNS (Rosales & Meneveau Reference Rosales and Meneveau2005).

This study highlights an important feature of this SFS model: the fact that the coefficient $C_{k}$ is really a user-defined parameter, determined by the SFS activity level chosen. The larger $s_{{\it\varepsilon}}$ (or any other measure of the SFSs used), the larger the coefficient and the more significant the contribution of the unresolved scales. As $C_{k}$ (and, consequently, $\ell$ ) are increased, a coarser grid is sufficient to resolve eddies of size comparable to the filter width, but modelling errors are expected to increase as well. An optimal value of the SFS activity may be case-dependent; in the present case, the LES was found to be accurate for scales larger than $2{\it\Delta}$ , independent of the value of ${\it\Delta}$ that is implied by the definition of $\ell$ chosen by the user.

It should be observed that, for $C_{k}=0.006$ , i.e.  $C_{{\it\Delta}}\approx 0.05$ , nearly 10 % of the TKE resides in the unresolved scales (figure 2 a). Approximating $\mathscr{K}$ in (3.10) with $\mathscr{K}_{res}$ may, therefore, result in errors in ${\it\nu}_{sfs}$ , particularly since this viscosity scales with $\mathscr{K}^{3}$ . To verify the significance of these errors, three additional simulations were performed. First, the exact value of $L$ (obtained from the DNS data) was used; second, we considered a case in which the TKE from the DNS was used in the numerator of (3.10); finally, $\mathscr{K}_{sfs}$ was estimated as (Yoshizawa Reference Yoshizawa1986)

(4.8) $$\begin{eqnarray}\mathscr{K}_{sfs}=0.0886(C_{{\it\Delta}}L_{est})^{2}|\overline{S}|,\end{eqnarray}$$

and was added to $\mathscr{K}_{res}$ in the numerator of (3.10). The results of the four simulations were within 2 % of each other, indicating that neglecting the SFS contribution to the TKE should not give significant errors when the SFSs responsible for less than 10 % of the TKE itself.

4.3. Plane channel flow

Channel flow simulations were performed at $\mathit{Re}_{{\it\tau}}=950$ and 2000. The domain size is $6{\it\delta}\times 2{\it\delta}\times 3{\it\delta}$ in the streamwise, wall-normal and spanwise directions, respectively. While Piomelli & Geurts (Reference Piomelli and Geurts2010, Reference Piomelli, Geurts, Kuerten, Geurts, Armenio and Fröhlich2011) determined $C_{k}$ by minimizing the error in the prediction of skin friction coefficient $C_{f}$ , which is known from correlation of experimental data (Dean Reference Dean1978), we chose $C_{k}$ to achieve a desired level of the subfilter activity $s_{{\it\varepsilon}}$ . The relationship of the chosen measure to others (the amount of Reynolds shear stresses or TKE provided by the model, for instance) will be discussed later.

Figure 3 showed the error minimization for channel flow at $\mathit{Re}_{{\it\tau}}=950$ and grid size $48\times 65\times 48$ illustrating the effectiveness of the SIPI procedure. Figure 7 shows the mean velocity profile (in wall units) and the turbulence intensities obtained at various grid resolutions. Approximate convergence appears with $48\times 65\times 48$ grid points for the mean velocity, while $128\times 129\times 128$ grid points are required to obtain near grid-independent Reynolds stresses. The fact that the mean-flow statistics can be predicted reasonably accurately with only $48\times 65\times 48$ grid points (corresponding to ${\rm\Delta}x^{+}\simeq 119$ and ${\rm\Delta}z^{+}\simeq 60$ ) is rather surprising: with this resolution the quasi-streamwise vortices in the wall layer cannot be resolved adequately, and LES is not expected to be accurate. This issue will be discussed further later.

Figure 7. Turbulence statistics for $\mathit{Re}_{{\it\tau}}=950$ and $C_{k}=0.007$ : (a) mean velocity; (b) $u_{rms}$ ; (c) $v_{rms}$ ; (d) $w_{rms}$ ; $---$   $48\times 65\times 48$ grid points; — ⋅ —  $64\times 97\times 64$ grid points; - - - -  $128\times 129\times 128$ grid points; ——  $192\times 193\times 192$ grid points; $+$ DNS (Hoyas & Jiménez Reference Hoyas and Jiménez2006).

Figure 8 shows contours of streamwise velocity fluctuations obtained with four grids at Reynolds number $\mathit{Re}_{{\it\tau}}=950$ . One can observe how refining the mesh results in improved resolution of finer scales, as expected. This grid-refinement study shows that there is very little difference between the results obtained at the two finer meshes. Apparently, for these grids, $\ell$ does not change much with resolution, as the integral scale is well-approximated throughout the domain. The smallest resolved eddies that result from the choice $C_{k}=0.007$ can be represented reasonably well using $128\times 129\times 128$ grid points.

Figure 8. Contours of $u$ -velocity fluctuations at $y^{+}=11$ superposed on the underlying grid with $\mathit{Re}_{{\it\tau}}=950$ and $C_{k}=0.007$ : (a) $48\times 65\times 48$ , (b) $64\times 97\times 64$ , (c) $128\times 129\times 128$ and (d) $192\times 193\times 192$ grid points.

We also conducted simulations at $\mathit{Re}_{{\it\tau}}\simeq 2000$ . The SIPI procedure was carried out on a $64\times 97\times 64$ grid using the estimate for subfilter activity based on the dissipation, and a target value $s_{tgt}=0.23$ . More refined calculations were subsequently carried out with the value of $C_{k}=0.005$ that was obtained from the SIPI algorithm. The turbulence statistics are shown in figure 9, and again demonstrate improved agreement with the DNS data as the mesh is refined while keeping $C_{k}$ fixed. This indicates that with the filter-width definition based on flow physics, as proposed here, the convergence to a grid-independent LES can be pursued effectively. At proper (sufficiently low) $C_{k}$ the grid-independent LES results will show close agreement with DNS data. As $C_{k}\rightarrow 0$ , the DNS limit would be achieved when the grid is refined.

Figure 9. Turbulence statistics for $\mathit{Re}_{{\it\tau}}=2000$ and $C_{k}=0.005$ : (a) mean velocity; (b) $u_{rms}$ ; (c) $v_{rms}$ ; (d) $w_{rms}$ ; — ⋅ —  $64\times 97\times 64$ grid points; $---$   $128\times 129\times 128$ grid points; - - - -  $192\times 193\times 192$ grid points; ——  $256\times 257\times 256$ grid points; $+$ DNS (Hoyas & Jiménez Reference Hoyas and Jiménez2006).

Profiles of several quantities related to the SFS model are shown in figure 10. Figure 10(a) shows the ratio of the eddy viscosity and the molecular viscosity. Some dependence on the grid resolution can be observed on coarse grids as the eddy viscosity converges with grid refinement. This dependence on resolution, however, is much weaker than what would be seen in a classical model of the Smagorinsky type, in which the viscosity depends on $h^{2}$ ; in the present case the grid resolution is changed by a factor of 4 in each direction between coarsest and finest simulations, which would result in a decrease of the Smagorinsky eddy viscosity by approximately a factor of 16. The eddy viscosity as defined in this paper only decreases by 33 % over the entire range of under-resolved to approximately grid-independent resolutions. Also, the two finest grids give well-converged results for ${\it\nu}_{sfs}$ despite the fact that the grid size is different by 50 %. The near-wall behaviour, as expected, is proportional to $y^{6}$ , steeper than the theoretical prediction of $y^{3}$ .

Figure 10. Profiles of: (a) SFS eddy viscosity; (b) estimated integral scale $L_{est}$ ; (c) ratio $L_{est}/h$ ; and (d) ratio $C_{k}L_{est}/C_{S}h$ ; channel flow, $\mathit{Re}_{{\it\tau}}=950$ and $C_{k}=0.007$ ; $---$   $48\times 65\times 48$ grid points; — ⋅ —  $64\times 97\times 64$ ; - - - -  $128\times 129\times 128$ ; ——  $192\times 193\times 192$ grid points.

Figure 10(b,c) show the integral length scale $L_{est}$ and its ratio to the grid size. First, we observe a rapid rise of $L_{est}$ near the wall followed by a more gradual increase in the channel core. As the mesh is refined the prediction of $\mathscr{K}$ and ${\it\varepsilon}_{tot}$ becomes more accurate, which results in changes in $L$ and ${\it\nu}_{sfs}$ ; a nearly grid-converged value of $L$ is reached for the $128\times 129\times 128$ grid. The length scale $L_{est}$ is of the order of the channel half-width in the core of the flow, reflecting the expected size of the eddies there. The mixing length $\ell =-\langle u^{\prime }v^{\prime }\rangle /|\text{d}U/\text{d}y|^{2}$ , in the core of the channel, is approximately equal to $0.15{\it\delta}$ , of the same order as $L_{est}$ . At this resolution, 20 grid points are required to resolve $L_{est}$ (as shown by the ratio $L_{est}/h$ in figure 10 c), at least in the channel core.

We have chosen, in this paper, to consider a global average of the SFS activity measures. This choice implies that, in inhomogeneous flows, the actual SFS contribution may be higher or lower than the target in different regions of the flow. To understand better the consequences of this choice, in figure 11 we show profiles of a surrogate measure of the SFS activity, plane by plane. This measure is defined as

(4.9) $$\begin{eqnarray}\widetilde{s}_{{\it\varepsilon}}=\frac{\langle {\it\varepsilon}_{sfs}\rangle _{xzt}}{\langle {\it\varepsilon}_{tot}\rangle _{xzt}}\frac{s_{{\it\varepsilon}}}{\displaystyle \int _{0}^{h}\frac{\langle {\it\varepsilon}_{sfs}\rangle _{xzt}}{\langle {\it\varepsilon}_{tot}\rangle _{xzt}}\,\text{d}y}\end{eqnarray}$$

where $\langle \cdot \rangle _{xzt}$ indicates averaging over wall-parallel planes and time, whereas $\langle \cdot \rangle$ indicates averaging over the entire volume (and time). The second term in (4.9) ensures that the average of $\widetilde{s}_{{\it\varepsilon}}$ is $s_{{\it\varepsilon}}$ . A similar expression defines $\widetilde{s}_{{\it\tau}}$ . The figure shows that throughout most of the channel the SFS activity is nearly constant. On the coarse grid (on which the SIPI procedure was carried out) its average value is very close to the target one; on finer grids it is somewhat lower: the grid-converged value of the integral scale is slightly higher than that predicted on the coarse grid (figure 10 d), and on the fine mesh the resolved eddies have a higher-than-expected contribution to the momentum and energy transport. The local SFS activity measures show significant variations only near the wall, for $y<0.1{-}0.2~(y^{+}<100{-}200)$ . This result indicates that in calculations of attached boundary layers (especially if wall models are used, in which the first grid point is located in the logarithmic layer) the use of global averaging may be adequate. On the other hand, in flows including separation or other non-equilibrium effects, use of a more local measure of SFS activity may be more suitable. Within that contest (and, in particular, in wall-modelled LES, in which the viscous dissipation is negligible compared with that supplied by the SFSs), we point out that $s_{{\it\tau}}$ is a better indicator than $s_{{\it\varepsilon}}$ . First, as $\mathit{Re}\rightarrow \infty$ one loses much sensitivity when $s_{{\it\varepsilon}}$ is used to determine $C_{k}$ . Second, requiring that the SFS contribution to the dissipation be equal to $s_{{\it\varepsilon}}$ at each point implies that ${\it\nu}_{sfs}={\it\nu}/(1-s_{{\it\varepsilon}})$ everywhere as well (whereas with the present approach, this is true only on average). Such a model cannot be expected to be accurate in most flows.

Figure 11. Profiles of plane-averaged SFS activity measures, defined in (4.9): (a) $\widetilde{s}_{{\it\varepsilon}}$ ; (b)  $\widetilde{s}_{{\it\tau}}$ ; $---$   $48\times 65\times 48$ grid points; — ⋅ —  $64\times 97\times 64$ ; - - - -  $128\times 129\times 128$ ; ——  $192\times 193\times 192$ grid points. The thin dotted line shows the target value.

Note that the filter width ${\it\Delta}=C_{{\it\Delta}}L$ is unknown: the constant $C_{{\it\Delta}}$ does not need to be determined as part of the simulation algorithm. The model coefficient $C_{k}$ also accounts for the (unknown) constant of proportionality for the time scale. One can estimate $C_{{\it\Delta}}$ , however, from figure 2, assuming that the turbulence is isotropic near the centre of the channel. The Reynolds number based on $L$ and the turbulence intensity $\mathscr{K}^{1/2}$ is approximately 800 in this region, and the subfilter activity calculated at grid convergence is approximately 0.17; figure 2 indicates that these values correspond to $C_{{\it\Delta}}\simeq 0.06$ . The filter width is, therefore, smaller than the integral scale by a factor of 15–20 (similar to the values observed for the simulations of HIT discussed before). For the finest mesh, in which $L/h\simeq 30$ (figure 10 c), this gives a ratio ${\it\Delta}/h\simeq 1.5{-}2$ , which is consistent with earlier observations of subfilter resolutions ${\it\Delta}/h$ that would be required to achieve nearly grid-independent numerical solutions (Vreman, Geurts & Kuerten Reference Vreman, Geurts and Kuerten1994; Lund Reference Lund2003). Another estimate of the ratio between filter width and grid size can be obtained by comparing the quantity $C_{k}L$ with the analogous quantity used in the Smagorinsky model, $C_{s}h$ . As is the case for $C_{k}$ , $C_{s}$ also includes the constant of proportionality for the time scale. The ratio between the two quantities is shown in figure 10(d). In the channel core this ratio is approximately equal to 1.5 for the $128\times 129\times 128$ grid, increasing to 2 for the more refined grid.

Figure 12 shows the variation of the two measures of subfilter activity defined above in (4.1) and (4.3) with $C_{k}$ . We also include a third quantity that measures the SFS contribution to the shear stress,

(4.10) $$\begin{eqnarray}s_{uv}=\frac{\left\langle {\it\nu}_{sfs}\unicode[STIX]{x1D61A}_{12}\right\rangle }{\left\langle -\overline{u}^{\prime }\overline{v}^{\prime }+{\it\nu}_{sfs}\unicode[STIX]{x1D61A}_{12}\right\rangle }.\end{eqnarray}$$

While this measure of SFS activity is not generally applicable, being coordinate-dependent, and difficult to generalize to complex geometries, it is useful to relate the measures introduced earlier in the literature to $s_{uv}$ , which gives a direct measure of the contribution of the SFSs to momentum transport. The figure shows, first, that the range of $s_{{\it\varepsilon}}$ considered (between 0 % and 40 %) corresponds to an SFS contribution to the total shear stress of less than 10 %; the corresponding values of $s_{{\it\tau}}$ are under 5 %. As we already observed, on fine grids all of the metrics considered have lower values than those estimated on the coarse mesh for the SIPI procedure used to assign $C_{k}$ . Thus, the grid-converged calculations resolve more eddies than expected; indicating that the procedure used to evaluate $C_{k}$ is conservative.

Figure 12. Measures of SFS activity; channel flow, $\mathit{Re}_{{\it\tau}}=950$ : (a) $s_{{\it\varepsilon}}$ ; (b) $s_{{\it\tau}}$ (c) $s_{uv}$ ; line, $48\times 65\times 48$ grid points; line with symbols, $128\times 129\times 128$ grid points.

It is possible also to choose the value of $C_{k}$ that optimizes the prediction of some quantity known experimentally. This is predicated on the availability of prior knowledge, and is effective only if the quantity of interest is not strongly dependent on the small scales, since the SIPI procedure is performed on a coarse grid. Piomelli & Geurts (Reference Piomelli and Geurts2010, Reference Piomelli, Geurts, Kuerten, Geurts, Armenio and Fröhlich2011), for instance chose the value of $C_{k}$ that minimized the error in the prediction of the skin-friction coefficient, $C_{f}$ . While we believe that the use of an internal measure such as $s_{{\it\tau}}$ is preferable, we also observe that the values of $C_{k}$ chosen here are close to the value that optimizes the prediction of $C_{f}$ . On the coarse mesh, in fact, the minimum $C_{f}$ is obtained for $C_{k}=0.0061$ , close to the value chosen for fine simulations ( $C_{k}=0.007$ ).

Next, we examined the sensitivity of the numerical results to the value of $C_{k}$ chosen. We performed simulations on the $128\times 129\times 128$ grid, using coefficients ranging from 0.005 to 0.012, corresponding to SFS contributions to the dissipation ranging between 10 % and 33 % (as evaluated on the fine grid) and to the shear stress between 1 % and 5 %. For SFS dissipation between 10 % and 25 % of the total we observe no difference in the mean velocity profiles (figure 13 a) and little difference in the r.m.s. turbulence intensities (figure 13 b; the other components have similar behaviours). When the filter width is increased so that the subfilter eddies contribute over 30 % of the dissipation, we begin to observe errors in the buffer region (see the inset in the figure). In particular, the point of maximum $u_{rms}$ moves away from the wall (from $y^{+}\simeq 30$ for $C_{k}=0.012$ to 13 for $C_{k}=0.005$ ), and the region of the peak becomes wider, indicating an overestimated diffusion due to a slightly overestimated eddy viscosity. From figure 2, one can estimate that this case corresponds to a filter width ${\it\Delta}\simeq 0.07{-}0.1L$ , substantially larger than that corresponding to $C_{k}=0.007$ . Even with such a large filter width, however, the differences with the DNS data are surprisingly small.

Figure 13. Profiles of (a) mean velocity and (b) $u_{rms}$ ; for $\mathit{Re}_{{\it\tau}}=950$ ; ——  $C_{k}=0.005$ ; $\cdots \cdots$   $C_{k}=0.007$ ; - - - -  $C_{k}=0.009$ ; — ⋅ —  $C_{k}=0.012$ ; $128\times 129\times 128$ grid points.

4.4. Comparison with other models and cost considerations

We next compare the present model with the dynamic eddy-viscosity model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992), which is in widespread use. We used the plane-averaged formulation with both grid and test filter widths proportional to the grid size, and performed calculations of the channel flow at $\mathit{Re}_{{\it\tau}}=950$ using both a coarse grid (with $48\times 65\times 48$ points) and a refined one using $128\times 129\times 128$ points. Figure 14 compares the mean velocity profiles and the streamwise r.m.s. turbulence intensity, $u_{rms}$ . First, the present model, on the coarsest mesh, gives results in surprisingly good agreement with the DNS, while the dynamic model results in an incorrect prediction of the near-wall flow, with an overestimation of the intercept of the logarithmic layer and an excessively high peak of $u_{rms}$ , both due to a 12 % underestimation of the wall stress. On a finer grid, on the other hand, both models give similar results. With $C_{k}=0.007$ , the SFS eddies are expected to contribute between 16 % and 23 % of the dissipation, and 2 %–4 % of the Reynolds shear stresses (figure 12).

Figure 14. Profiles of (a,c) mean velocity and (b,d) $u_{rms}$ for $\mathit{Re}_{{\it\tau}}=950$ ; —— present model, $C_{k}=0.007$ ; — ⋅ — dynamic eddy-viscosity model; $+$ DNS (Hoyas & Jiménez Reference Hoyas and Jiménez2006); (a,b) $48\times 65\times 48$ grid points; (c,d) $128\times 129\times 128$ grid points.

Figure 15 shows that, on the coarse mesh, the present model yields an eddy viscosity that is much larger than that predicted by the dynamic model. This balances the inability of the resolved scales to transport momentum at the appropriate rate: the fluctuation contours shown in figure 8 show streaks of much larger scale than physical, and a more quiescent flow. Thus, the contribution of the SFS to the Reynolds stresses is enhanced (figure 15 b), and exceeds that of the resolved scales in the viscous and buffer layers. In fact, the total Reynolds stress (sum of resolved and SFS components) may be excessively high in the buffer region; the mean velocity gradient must, in that region, be lower, to satisfy the momentum balance; this phenomenon is noticeable in figure 14(a). With the dynamic model, on the other hand, the subfilter stresses are insufficient to balance the unresolved momentum, and the mean velocity gradient must increase to close the momentum balance, leading to the high intercept of the logarithmic layer. On the fine mesh, on the other hand, the dynamic model gives accurate results indicative of the high fidelity achieved at such resolutions.

Figure 15. Profiles of (a,c) eddy viscosity and (b,d) shear stress; channel flow, $\mathit{Re}_{{\it\tau}}=950$ . —— Present model, $C_{k}=0.007$ ; - - - - dynamic eddy-viscosity model; ●  $\mathscr{K}{-}{\it\omega}$ turbulence model in (a). (a,c): $48\times 65\times 48$ grid points. (b,d): $128\times 129\times 128$ grid points. For the shear stress figures: lines: resolved stress; lines with squares: SFS stress; lines with circles: total Reynolds stress.

SFS models usually require the energy- and momentum-transporting eddies to be resolved to give accurate results. On the coarse mesh, such is not the case, and the dynamic model cannot be expected to be accurate. The model proposed here, on the other hand, has some features of a RANS turbulence model, through the definition of the integral scale. Figure 15(a) shows the eddy-viscosity obtained for this case by solving the RANS equations with the $\mathscr{K}{-}{\it\omega}$ turbulence model (Wilcox Reference Wilcox1993). In the viscous and buffer layers the present model yields virtually the same eddy viscosity as that predicted by the RANS turbulence model. At finer resolution the present model behaves like a standard SFS model, yielding a much lower viscosity than the $\mathscr{K}{-}{\it\omega}$ model. It is not to be expected, however, that this behaviour could extrapolate to cases in which $s_{{\it\varepsilon}}\rightarrow 1$ : unlike the PANS approach, which solves explicitly an equation for the unresolved TKE, the present model requires that $\mathscr{K}_{tot}\simeq \mathscr{K}_{res}$ . Even on the coarse mesh, the resolved eddies support 95 % of the shear stress, integrated over the entire channel (figure 12 c). Presumably, they will support roughly the same percentage of TKE, leading to an accurate prediction of the dissipation length. For $s_{{\it\varepsilon}}=0.3$ and $s_{{\it\tau}}=0.03$ errors in the prediction of the r.m.s. fluctuations become more significant, even on the fine mesh; these values may approach the limit of applicability of the present model.

An important issue to be considered when developing SFS models is their cost. The model proposed here, in terms of computational cost, is comparable to the Smagorinsky model, and significantly cheaper than dynamic models that require 60 % more CPU time per time step and grid point due to the repeated filtering operations required. The proposed model requires, of course, five or six precursor simulations to determine a reliable estimate of the optimum value of the coefficient $C_{k}$ . Since they are performed on a grid that is much coarser than that used for production runs, however, this overhead is not excessive. Tables 1 and 2 show the details of the simulation cost using the proposed SFS model and other models, respectively, for the $\mathit{Re}_{{\it\tau}}\simeq 1000$ case; all of the simulations were run for $10{\it\delta}/u_{{\it\tau}}$ on eight 2.52 GHz Sparc64 VII processors.

Table 1. Simulation cost for channel flow at $\mathit{Re}_{{\it\tau}}\simeq 1000$ using the proposed SFS model on 8 processors.

Table 2. Simulation cost for channel flow at $\mathit{Re}_{{\it\tau}}\simeq 1000$ using 8 processors and other commonly used models.

The cost of the simulations on the fine mesh for the present model was only 6 % higher compared with a simulation using the Smagorinsky model, and nearly 20 % lower than the dynamic model. The precursor simulations, which were performed on a much coarser grid added only 4 % to the total cost. Including the cost of the precursor simulations, the present model is only marginally more expensive than the Smagorinsky model, and less costly than the dynamic model. These considerations depend to a very large extent on the algorithm used: the present code has a very efficient direct solver for the Poisson equation, which requires a fraction of the total CPU time per step; in codes in which less efficient solvers are used, especially in complex geometries, the Poisson solution may require a much larger fraction of the CPU time, and the cost of the SFS model evaluation may not be a significant consideration.