4.1. Energy transfer in physical space

For isotropic turbulence and corresponding spherical, sharp spectral cutoff filter, the total SGS dissipation at cutoff $k_c$ is related to the physical space calculation as

(4.1) \begin{equation} \widehat {\varepsilon }_{\textit{SGS}}(k_c) = \frac {1}{N^3} \sum _{\boldsymbol {x}} \varepsilon _{\textit{SGS}}(\boldsymbol {x}), \quad \text{ where } \varepsilon _{\textit{SGS}}(\boldsymbol {x}) = -\tau _{ij}^{\textit{SGS}} S^\lt _{ij}, \end{equation}

with summation over all $N^3$ grid points in $\boldsymbol {x}$ . Unfortunately, there is no equivalent spectral space formula for the local SGS dissipation $\varepsilon _{{SGS}}(\boldsymbol {x})$ . Nevertheless, analogous to the scale separation of $k \leqslant a k_c$ in the previous section, we first restrict the estimation of eddy viscosity to subfiltered (SF) scales. Specifically, for general implementations in physical space, the kinetic energy transfer for a test-filtered field can be computed based on the subfilter scales (SFS)

(4.2) \begin{equation} \tau _{ij}^{{SFS}} = \tau _{ij}^{\text{sim}} = \overline {u^\lt _i u^\lt _j} - \overline {u}^\lt _i \overline {u}^\lt _j, \end{equation}

where ‘ $\bar {\enspace }$ ’ denotes a test-filtering, whose characteristic filtering width $\overline {\unicode{x1D6E5} }$ is usually larger than the LES mesh size $\unicode{x1D6E5}$ . Note that the current SFS stress is based on the resolved LES scales, which differs from the ‘unresolved’ SFS studied by e.g. Rouhi et al. (Reference Rouhi, Piomelli and Geurts2016). Specifically, the SFS stress distinguishes between the filter scale $\overline {\unicode{x1D6E5} }$ and the LES grid scale $\unicode{x1D6E5}$ ; the splitting of the former is performed only on resolved scales with either spectral or smooth filter. A detailed description of SFS stress and its applicability in LES is discussed in detail in the monograph of Berselli et al. (Reference Berselli, Iliescu and Layton2006), among others. Note that the subfilter stress is exactly the same as in Bardina’s similarity model (Bardina et al. Reference Bardina, Ferziger and Reynolds1980, Reference Bardina, Ferziger and Reynolds1983). In general, the similarity model is multiplied by a constant coefficient $C_{\text{sim}}$ , which is $O(1)$ based on a priori studies (Liu et al. Reference Liu, Meneveau and Katz1994) and theoretical analysis (Cook Reference Cook1997). Similarity-type models can achieve very high correlations in a priori analyses, but in a posteriori tests, due to the large backscatter, they are widely reported to be inadequately dissipative. In actual simulations with relatively low grid resolutions, it is normally necessary to apply ad hoc procedures to reduce the backscatter and stabilise the simulations (see e.g. Zang et al. Reference Zang, Street and Koseff1993; Liu et al. Reference Liu, Meneveau and Katz1994).

Similar to the SGS dissipation, the test-filtered energy transfer can also be computed using standard LES formulae based on the SFS tensor

(4.3) \begin{equation} \widehat {\varepsilon }_{{SFS}}(ak_c) \approx \frac {1}{N_x N_y N_z} \sum _{\boldsymbol {x}} \varepsilon _{{SFS}}(\boldsymbol {x}),\quad \text{ where } \varepsilon _{{SFS}}(\boldsymbol {x}) = -\tau _{ij}^{{SFS}} \overline {S}^\lt _{ij}, \end{equation}

in which the filtered velocity field is used to compute the strain rate. Following the convention in LES, a positive value of $\varepsilon$ indicates energy dissipation. Note that the concept of spectral eddy viscosity in the previous section is based on sharp spectral cutoff filtering, which is a natural choice for the Fourier spectral method. However, filtering in physical space leads to complex spectral support, thus, an exact mathematical relation between characteristic filtering width $\overline {\unicode{x1D6E5} }$ and the sharp spectral cutoff $a k_c$ does not exist. Nevertheless, the above expression suggests that $\varepsilon _{{SFS}}(\boldsymbol {x})$ is analogous to $-\widehat {T}^{res}(ak_c)$ in (3.9), implying that it serves a similar role in SGS modelling. Moreover, test-filtering is essential from the perspective of multi-scale energy transfer. Recently, Domaradzki (Reference Domaradzki2024) demonstrated that for HIT, the ITM based on a spectral cutoff filter and physical space filters with the same characteristic filtering width produces very similar results. This motivates further investigation into the applicability of scale separation in the context of physical space filtering. In particular, to establish a more direct connection with the spectral cutoff, it may be preferable to estimate the characteristic filtering width based on the filtering kernel in spectral space, which will be discussed in greater detail later in § 6.

In addition, when the subfilter stress or the similarity model (4.2) is used as an SGS model, it is multiplied by the rate-of-stress tensor $S^\lt _{ij}$ without test-filtering when computing the SGS dissipation (4.3). Since the SF scales are primarily defined for $k \geqslant ak_c$ , the dissipation tends to be more dominant around $ak_c$ rather than at the true LES cutoff $k_c$ . The resulting backscatter at moderate wavenumbers near the filter cutoff $ak_c$ is speculated to be a key factor contributing to numerical instabilities and inaccurate predictions for the similarity model (Domaradzki Reference Domaradzki2021a ). Therefore, to avoid the scale mismatch, we compute the total energy transfer only for the test-filtered scales. Following the idea of the classic test field model (Kraichnan Reference Kraichnan1976) and the spectral eddy viscosity in § 3, here, we convert the SF dissipation to the eddy viscosity form while maintaining the same total energy transfer

(4.4) \begin{equation} \overline {\nu }_t(\boldsymbol {x}) = C \frac {\varepsilon _{{SFS}}(\boldsymbol {x})}{ |\overline {S}|^2(\boldsymbol {x})},\quad \text{ where } |\overline {S}| = \sqrt {2\overline {S}^\lt _{ij} \overline {S}^\lt _{ij}}, \end{equation}

which leads to an effective eddy viscosity on an LES grid. From the perspective of energy transfer in LES, the coefficient $C$ should be chosen such that the SGS dissipation computed using (4.4) in (2.5) correctly approximates the actual total SGS energy transfer in (2.4). It was found by Domaradzki (Reference Domaradzki2024) that for common Gaussian and top-hat filters in the simulations of HIT, the coefficient $C$ was of the order of unity.

In practice, a convolution filter used for scale separation in physical space is defined as

(4.5) \begin{equation} \overline {u}^\lt _i(x) = \int ^\infty _{-\infty } G(x - x') u^\lt _i(x') {\rm d}x', \qquad \overline{u}^{\prime}_{i} = u^\lt _{i} - \overline {u}^\lt _i, \end{equation}

where $\overline{u}^{\prime}_{i}$ is the subfilter velocity. For non-uniform meshes in wall-normal direction, we apply locally a quadratic interpolation to uniformly spaced points before applying the above formulae (Loh & Domaradzki Reference Loh and Domaradzki1999), see the descriptions of numerical methods in § 5. Then, the evolution equation for filtered velocity becomes

(4.6) \begin{equation} \partial _t \overline {u}^\lt _i + \partial _j(\overline {u}^\lt _i \overline {u}^\lt _j + \overline {\tau _{ij}^{{SGS}}} + \tau _{ij}^{{SFS}}) = - (1/\rho ) \partial _i \overline {p}^\lt + 2 \nu \partial _j \overline {S}^\lt _{ij}, \end{equation}

in which the commutation error between filtering and differentiation is absorbed into $\tau _{ij}^{{SGS}}$ and $\tau _{ij}^{{SFS}}$ . The latter can be computed directly based on (4.2) or reconstructed (Bull & Jameson Reference Bull and Jameson2016) in simulations. However, in many works, the SF stress (or equivalently, the similarity model) is the only explicit term added to the LES equation; the real SGS effects ( $\tau _{ij}^{{SGS}}$ or $\overline {\tau _{ij}^{{SGS}}}$ ) are approximated by numerical dissipation, which is challenging to control for low-order schemes and almost negligible for higher-order ones such as the spectral methods. Furthermore, the various implementations of filtering lead to additional uncertainties. In fact, it is widely reported that when only the SF stress (i.e. the similarity model) is used, it fails to provide adequate energy dissipation (Meneveau & Katz Reference Meneveau and Katz2000). To overcome these issues, our recent approach (Domaradzki Reference Domaradzki2024) splits the task of modelling energy transfer at resolved scales into the following steps:

1. (i) select a physical space filtering procedure, denoted as ‘ $\bar {\enspace }$ ’;

2. (ii) based on a LES mesh, compute the subfilter energy transfer/dissipation $\varepsilon _{{SFS}}$ at resolved scales, (4.3);

3. (iii) compute the eddy viscosity $\overline {\nu }_t$ based on the resolved energy transfer from (4.4) and constant $C=1$ ;

4. (iv) advance the LES equation (2.2) in time employing SGS stress (2.5) with the eddy viscosity $\overline {\nu }_t$ .

Note that unlike the spectral space implementations outlined in § 3, we have not yet taken care of the shape of eddy viscosity among scales. In the following, we show that the above modelling strategy, despite its simplicity, is capable of recovering major features of spectral interscale energy transfer model.

4.2. Scale distribution of energy transfer

Based on the concept of spectral eddy viscosity in § 3, ideally, an eddy viscosity model should exhibit a plateau at low wavenumbers and a cusp near the LES cutoff $k_c$ . The former helps to maintain the desired Kolmogorov −5/3 law in the inertial range, while a cusp helps to damp the non-physical energy accumulation at small scales and stabilise the simulations. However, it is not possible to recover the precise shape of spectral eddy viscosity for inhomogeneous flow. In a series of recent papers (Domaradzki Reference Domaradzki2021a , Reference Domaradzkib , Reference Domaradzki2022a , Reference Domaradzkib , Reference Domaradzki2024), it was shown that as long as the total SGS transfer was enforced, the details of its wavenumber distribution (i.e. its $k$ -dependence) were of secondary importance. In fact, it is observed that the shape functions with different plateau levels and cusps produced largely the same results for energy spectra when the total transfer was computed using (3.9) (Domaradzki Reference Domaradzki2022b ). Therefore, in this subsection, rather than finding a direct proxy for the spectral shape function $f$ , we show that the eddy viscosity (4.4) can produce a similar plateau–cusp shape at desired scales.

The original interscale model is developed based on an additional spectral cutoff in (3.9), in which the evolution of spectral energy $T_{{SGS}}$ at LES scales $k \leqslant ak_c$ is considered. Nevertheless, filtering in physical space normally generates broad spectral support (see e.g. Sagaut Reference Sagaut2006), indicating that a complete scale separation cannot be achieved in the same way; (4.6) still includes an aggregate effect of transfer due to non-local interactions between small and large scales. A common strategy in LES modelling is to create a cusp at small scales near the LES cutoff, see figure 1. From the perspective of multi-scale energy transfer, a cusp near $k_c$ is equivalent to introducing large dissipation to reduce energy at small resolved scales. Thus, instead of focusing on the transport equation of large scale kinetic energy $\overline {E} =({1}/{2}) \overline {u}^\lt _i \overline {u}^\lt _i$ , we first study the subfiltered scales. The evolution equation for subfilter energy $\overline {E}' = ({1}/{2}) \overline{u}^{\prime}_{i} \overline{u}^{\prime}_{i}$ is

(4.7) \begin{equation} \partial _t \left (\frac {1}{2} \overline{u}^{\prime}_{i} \overline{u}^{\prime}_{i} \right ) = T(\overline {u}^\lt _i, \overline{u}^{\prime}_{i}) - \overline {\nu }_t \partial _j \overline{u}^{\prime}_{i} \partial _j \overline{u}^{\prime}_{i} - \nu \partial _j \overline{u}^{\prime}_{i} \partial _j \overline{u}^{\prime}_{i}, \end{equation}

where $T(\overline {u}^\lt _i, \overline{u}^{\prime}_{i})$ represents total energy transfer. In spectral space implementation for HIT (Domaradzki Reference Domaradzki2021a ), a cusp of spectral eddy viscosity observed in the well-known Chollet–Lesieur model (Chollet & Lesieur Reference Chollet and Lesieur1981) is capable of controlling energy accumulation at small scales, thus stabilising the flow at large $Re$ and reproducing DNS results accurately. However, subsequent investigations (Domaradzki Reference Domaradzki2022b ) have shown that it may not be necessary to reproduce the exact shape of spectral eddy viscosity very accurately, as long as the plateau–cusp shape is roughly maintained, see figure 1. Additionally, the concept of reducing energy accumulation near $k_c$ has also been used in our previous works (Sun & Domaradzki Reference Sun and Domaradzki2018; Sun & Yang Reference Sun and Yang2023), leading to accurate results for wall-bounded flows. This suggests that an exact reproduction of the spectral eddy viscosity (3.6) in physical space may not be necessary, we simply need to retain the essential effects of the plateau–cusp shape.

To generate a plateau and a cusp across different spectral regions, we follow the minimum-dissipation criterion proposed by Verstappen (Reference Verstappen2011) to evalute the transport equation of $\overline {E}'$ . By equating the first two terms on the right-hand side of (4.7), the actions of small scales will be purely dissipative due to the molecular dissipation. In other words, the energy of the subfilter scales would decay exponentially, akin to the effect of a cusp in the spectral eddy viscosity (figure 1). This also ensures that the subfiltered energy is independent of large, resolved scales, which helps to form a plateau. However, due to the lack of widely acknowledged physical space results analogous to the Chollet–Lesieur model (Chollet & Lesieur Reference Chollet and Lesieur1981), we do not require $({1}/{2}) \overline{u}^{\prime}_{i} \overline{u}^{\prime}_{i}$ to decay at precisely a certain rate, which is also not crucial based on our previous observations (Domaradzki Reference Domaradzki2021a,b ).

To circumvent the difficulty in directly analysing small-scale energy, following Verstappen (Reference Verstappen2011), the Poincaré–Wirtinger inequality is applied,

(4.8) \begin{equation} \int _{\Omega _\delta } \frac {1}{2} \overline{u}^{\prime}_{i} \overline{u}^{\prime}_{i} {\rm d}x \leqslant \int _{\Omega _\delta } \frac {1}{2} (\overline {\partial }_i \overline {u}^\lt _j) (\overline {\partial }_i \overline {u}^\lt _j) {\rm d}x, \end{equation}

where $\overline {\partial }_i = \sqrt {C_i} \overline {\unicode{x1D6E5} }_i \partial _i$ is the scaled gradient operator and $C_i$ is a modified Poincaré constant. Therefore, an exponential decay of the above gradient term serves as a necessary condition for the subfilter energy $\overline {E}' = ({1}/{2}) \overline{u}^{\prime}_{i} \overline{u}^{\prime}_{i}$ in (4.7) to decay exponentially. To derive the transport equation for the right-hand side, assume the subfiltered quantities ( $\overline {u}^\lt _i$ , $\overline {p}^\lt$ ) obey a periodic boundary condition on a box filter domain $\Omega _\delta$ , which is mostly homogeneous near the smallest scales of motion (Verstappen Reference Verstappen2011; Verstappen et al. Reference Verstappen, Rozema and Bae2014). Similarly, based on the local isotropic assumption of subgrid scales (Moser et al. Reference Moser, Haering and Yalla2021), the eddy viscosity $\overline {\nu }_t$ that models small-scale actions is also assumed to be locally constant. Then, multiply (4.6) by $\overline {\partial }_j \overline {\partial }_j \overline {u}^\lt _i$ ,

(4.9) \begin{align} \partial _t \left (\frac {1}{2} (\overline {\partial }_i \overline {u}^\lt _j)(\overline {\partial }_i \overline {u}^\lt _j) \right ) = & -(\overline {\partial }_k \overline {u}^\lt _i) (\overline {\partial }_k \overline {u}^\lt _j) \overline {S}^\lt _{ij} - (\overline {\partial }_k \overline {S}^\lt _{ij}) \overline {\partial }_k(2\overline {\nu }_t \overline {S}^\lt _{ij}) - 2\nu (\overline {\partial }_k \overline {S}^\lt _{ij}) \overline {\partial }_k(\overline {S}^\lt _{ij})\nonumber\\ & +\, \partial _i f_i, \end{align}

where the subfiltered terms are modelled using the eddy viscosity. Here, $\partial _i f_i$ is a combination of flux terms that vanish due to the periodic boundary condition on $\Omega _\delta$ . Note that the term on the left-hand side serves as a surrogate to the subfiltered energy based on (4.8). To further simplify the eddy viscosity term in (4.9), we again apply the Poincaré inequality

(4.10) \begin{equation} \int _{\Omega _\delta } (\overline {\partial }_k \overline {S}^\lt _{ij}) (\overline {\partial }_k \overline {S}^\lt _{ij}) {\rm d}x \geqslant \int _{\Omega _\delta } \overline {S}^\lt _{ij} \overline {S}^\lt _{ij} {\rm d}x = \frac {1}{2} \int _{\Omega _\delta } |\overline {S}|^2 {\rm d}x. \end{equation}

Consequently, the production of small-scale energy $\overline {E}' = ({1}/{2}) \overline{u}^{\prime}_{i} \overline{u}^{\prime}_{i}$ can be prevented by comparing the transfer of small-scale energy and the eddy viscosity term (i.e. the first two terms on the right-hand side of (4.9)) in combination with (4.10),

(4.11) \begin{equation} \overline {\nu }_t \int _{\Omega _\delta } |\overline {S}|^2 {\rm d}x \overset {?}{\geqslant } - \int _{\Omega _\delta } (\overline {\partial }_k \overline {u}^\lt _i) (\overline {\partial }_k \overline {u}^\lt _j) \overline {S}^\lt _{ij} {\rm d}x. \end{equation}

In pratice, it is natural to take the minimum eddy dissipation that satisfies the above condition, which is analogous to generating a constant plateau for filtered scales at relatively low wavenumbers,

(4.12) \begin{equation} \overline {\nu }_t = \frac {\int _{\Omega _\delta } (\overline {\partial }_k \overline {u}^\lt _i) (\overline {\partial }_k \overline {u}^\lt _j) \overline {S}^\lt _{ij} {\rm d}x}{\int _{\Omega _\delta } -|\overline {S}|^2 {\rm d}x}. \end{equation}

The above equation is similar to the AMD model proposed by Rozema et al. (Reference Rozema, Bae, Moin and Verstappen2015). A notable difference being the denominator, in which the nominal dissipation $\varepsilon _\nu$ is used instead of the pseudo-dissipation,

(4.13) \begin{equation} \varepsilon _\nu = 2 \nu \overline {S}^\lt _{ij} \overline {S}^\lt _{ij} = \nu \partial _j \overline {u}^\lt _i \partial _j \overline {u}^\lt _i + \nu \frac {\partial ^2 \overline {u}^\lt _i \overline {u}^\lt _j}{\partial x_i \partial x_j} = \varepsilon _\nu ^{\text{pseudo}} + \nu \frac {\partial ^2 \overline {u}^\lt _i \overline {u}^\lt _j}{\partial x_i \partial x_j}, \end{equation}

where the last term is small in most circumstances (Pope Reference Pope2000), especially after filtering.

Furthermore, we relate (4.12) to the SF energy transfer in the previous subsection. As discussed by e.g. Vreman et al. (Reference Vreman, Geurts and Kuerten1996), Thiry & Winckelmans (Reference Thiry and Winckelmans2016) and Abkar & Moin (Reference Abkar and Moin2017), a series expansion of the similarity model based on the box or Gaussian filtering leads to

(4.14) \begin{equation} \tau _{ij}^{{SFS}} = \tau _{ij}^{\text{sim}} = \overline {u^\lt _i u^\lt _j} - \overline {u}^\lt _i \overline {u}^\lt _j = \frac {1}{12} (\overline {\unicode{x1D6E5} }_k \partial _k \overline {u}^\lt _i) (\overline {\unicode{x1D6E5} }_k \partial _k \overline {u}^\lt _j) + O(\overline {\unicode{x1D6E5} }_k^4). \end{equation}

Therefore, the gradients in the numerator of (4.12) can be replaced by the similarity model if higher-order terms are neglected. Compared with the scale separation based on sharp spectral cutoff, the similarity model computed using smooth filters actually exhibits higher correlations in a priori analysis (Liu et al. Reference Liu, Meneveau and Katz1994), supporting the implementation in physical space. Notably, the leading-order term is the well-known gradient model (or the nonlinear part of the tensor-diffusivity model, sometimes also referred to as the Clark model (Clark et al. Reference Clark, Ferziger and Reynolds1979). Akin to the similarity model, the gradient model is known to correlate well with the SGS stress, but tends to cause instability in practical simulations (Vreman et al. Reference Vreman, Geurts and Kuerten1996). For this reason, it is often applied with zero-clipping (Vreman et al. Reference Vreman, Geurts and Kuerten1997; Balarac et al. Reference Balarac, Sommer, Meunier and Vollant2013) or combined with the eddy viscosity model to generate mixed models (Clark et al. Reference Clark, Ferziger and Reynolds1979; Fabre & Balarac Reference Fabre and Balarac2011). As we will discuss later in § 6, recasting the similarity model to the eddy viscosity form generates a similar effect. Moreover, it is worthwhile to note that the factor ${1}/{12}$ in the gradient model can be changed if a different filter is chosen, so that the modelling may become suboptimal. Following the notation in the previous section, for generality, an additional modelling coefficient $C$ should be added to (4.12) to account for various filters.

Recall that the objective of this subsection is to demonstrate that the eddy viscosity (4.4) can generate the desired plateau–cusp characteristics. To reproduce the expression, we combine (4.12) and (4.14) to obtain an eddy viscosity model that prevents the production of small-scale energy $\overline {E}'$ , achieving an effect similar to a spectral eddy viscosity cusp at high wavenumbers. In actual implementations, it may not be necessary or intuitive to apply the additional box filter or cell averaging. Following a previous study (Rozema et al. Reference Rozema, Bae, Moin and Verstappen2015), one can apply the mid-point rule for integration,

(4.15) \begin{equation} \overline {\nu }_t \approx \frac {\int _{\Omega _\delta } - C \tau _{ij}^{{SFS}} \overline {S}^\lt _{ij} {\rm d}x}{\int _{\Omega _\delta } |\overline {S}|^2 {\rm d}x} \approx \frac { - C \tau _{ij}^{{SFS}} \overline {S}^\lt _{ij}}{|\overline {S}|^2} = C \frac {\varepsilon _{{SFS}}}{ |\overline {S}|^2}. \end{equation}

Thus, the eddy viscosity based on the interscale energy transfer concept (4.4) is recovered, which further confirms that concepts behind the original spectral ITM can be transfered directly to the physical space. In practice, the box filter can be replaced by an averaging procedure to enhance numerical stability.

Finally, as we discussed before, the remaining modelling coefficient $C$ in (4.15) depends on the choice of filtering. In the following large eddy simulations in § 6, by selecting smooth filters with characteristic filtering width $\overline {\unicode{x1D6E5} }$ close to the mesh size, the cusp region is limited to the high wavenumbers and of less importance to the turbulent statistics of interest in LES. In this way, the remaining subfiltered scales are almost negligible and there is no need to consider the rescaling coefficient. In addition, note that in our previous works (Sun & Domaradzki Reference Sun and Domaradzki2018), it is shown that the exact shape (and even the time interval) of the dissipation is not crucial, as long as the total SGS dissipation effectively prevents the non-physical energy pileup near $k_c$ . The same observation also holds for the simulations of HIT based on the interscale model (Domaradzki Reference Domaradzki2021a , Reference Domaradzkib , Reference Domaradzki2024). Therefore, we set $C=1$ and do not impose additional constraints.

4.3. Summary of the present approach

In summary, the new interscale model offers several clear advantages over traditional eddy viscosity models or similarity-type models.

1. (i) The current method is conceptually inherited from the classic spectral eddy viscosity (Kraichnan Reference Kraichnan1976), providing a more solid theoretical foundation for constructing SGS dissipation based on the test-filtered velocity field.

2. (ii) Recasting the model into an eddy viscosity formulation (4.15) imposes an additional constraint, serving as a straightforward way to enhance the numerical robustness of the similarity model.

3. (iii) The model reproduces effects similar to the plateau–cusp features predicted by EDQNM theory (see figure 1), mitigating the scale mismatch issue commonly observed in similarity-type models.

In addition, the theoretical analyses and satisfactory simulation results for homogeneous, isotropic turbulence (Domaradzki Reference Domaradzki2021a , Reference Domaradzkib , Reference Domaradzki2024) at high Reynolds numbers give us more confidence to extend the modelling strategy to more complex flows. As we will see in the a posteriori analysis in § 6, the ability to generate an appropriate amount of SGS dissipation among different scales is essential for the effectiveness of an SGS model.