Autonomous large-eddy simulations of turbulence using eddy viscosity derived from the subgrid-scale similarity stress tensor

Abstract

A previously developed method for large-eddy simulations (LES), based on spectral eddy-viscosity models, is generalised to the physical space representation. The method estimates the subgrid-scale (SGS) energy transfer using a similarity-type model expression for the SGS tensor obtained using Gaussian filtering of velocity fields advanced in the simulations. Following steps for the spectral space representation, the SGS transfer in the physical space is used to obtain a spatially varying eddy viscosity at each time step in LES. The computed eddy viscosity is employed to model the SGS stress tensor in the familiar Boussinesq form for use in LES. The method is tested in LES of isotropic turbulence at high Reynolds numbers where the inertial range dynamics is expected and for lower Reynolds number decaying turbulence under conditions of the classical Comte-Bellot and Corrsin experiments. In both cases the agreement with reference data is very good and the SGS transfer computed for the proposed eddy-viscosity model is highly correlated with the transfer computed for the similarity-type stress tensor.

1. Introduction

The fundamental equations of large-eddy simulation (LES) are obtained by the filtering approach proposed by Leonard (1974). A filtering operation for a function $f(\boldsymbol {x})$ is a linear operation, denoted by an overbar, and is defined as a convolution integral with a given smoothing kernel $G$

(1.1)\begin{equation} \bar{f}(\boldsymbol{x})=\int G(\boldsymbol{x}-{\boldsymbol{x}}^{\prime};\varDelta) f({\boldsymbol{x}}^{\prime})\,{\rm d}\kern0.7pt{\boldsymbol{x}}^{\prime}, \end{equation}

where $\varDelta$ is a prescribed filter width. Applying a filtering operation to the incompressible Navier–Stokes equations for velocity $u_i$ results in Navier–Stokes equations for the filtered velocity $\bar {u}_i$ that contain subgrid-scale (SGS) stress tensor

(1.2)\begin{equation} \tau_{ij}=\overline{u_{i}u_{j}}-\bar{u}_{i}\bar{u}_{j}. \end{equation}

This term represents the effects of SGSs of turbulence which are unresolved by the LES mesh. A SGS model accounting for the effects of the small unresolved scales of turbulence is typically introduced directly into the governing equations in place of $\tau _{ij}$. Traditional SGS models fall into three general categories: eddy-viscosity models, similarity models and so-called mixed models which combine eddy viscosity and similarity expressions. There are a number of excellent reviews of theory and practice of SGS modelling in the traditional framework, e.g. by Lesieur & Metais (1996), Piomelli (1999), Meneveau & Katz (2000), Pope (2000) and Sagaut (2002).

Perhaps the most important role of a SGS model is its ability to account for the energy flux from large, resolved scales, to small, unresolved SGS scales, which physically appears as a dissipation of kinetic energy of the resolved scales. Because of that SGS models that have received the most attention in the literature on the subject are dissipative eddy-viscosity models. The eddy-viscosity concept goes back to observation by Boussinesq in late 19th century that there is a similarity between turbulent transport and molecular transport. Mathematically, that similarity leads to the turbulent stress tensor having the same functional representation as the viscous stress tensor, but with a different, flow-dependent eddy viscosity. Smagorinsky (1963) proposed the first SGS model for LES that assumed the SGS stresses obey such a gradient–diffusion relationship. Since that pioneering work, many SGS models developed in this framework have followed.

A different observation, concerning similarity of SGS stress tensors computed from fully resolved and filtered velocity fields, led to an alternative family of SGS models, self-similarity models. In proposing the first of such models, Bardina, Ferziger & Reynolds (1983) followed Leonard (1974) decomposition of the SGS stress into contributions from the SGS Reynolds stress $(\overline {u_{i}^{\prime }u_{j}^{\prime }})$, cross-stress $(\overline {\bar {u}_{i}u_{j}^{\prime }}+\overline {u_{i}^{\prime }\bar {u}_{j}})$ and the Leonard stress $(\overline {\bar {u}_{i}\bar {u}_{j}}-\bar {u}_{i}\bar {u}_{j})$ terms. The final form of the scale-similarity model that follows from the work of Bardina et al. (1983) is written as

(1.3)\begin{equation} \tau_{ij}^{sim}=\overline{\bar{u}_{i}\bar{u}_{j}}-\bar{\bar{u}}_{i}\bar{\bar{u}}_{j}. \end{equation}

Using a Gaussian-type filter, Bardina et al. (1983) demonstrated in a priori analyses of decaying homogeneous isotropic turbulence and in sheared homogeneous turbulence, that the similarity model predictions were highly correlated with the actual SGS quantities, computed from the fully resolved direct numerical simulation (DNS) fields, with correlation coefficients significantly larger than for the Smagorinsky model. Another noted benefit of the similarity model was its ability to predict the reverse energy transfer, so-called backscatter, which is common in turbulence. The similarity concept was generalised by Liu, Meneveau & Katz (1994). The generalised form of the similarity model is

(1.4)\begin{equation} \tau_{ij}= c_{L}(\widehat{\bar{u}_{i}\bar{u}_{j}}-\hat{\bar{u}}_{i}\hat{\bar{u}}_{j}), \end{equation}

where the hat represents an explicit test filter with the filter width larger than for the bar filter, and $c_L$ is a model constant. Liu et al. (1994) conducted model testing using a model constant of $c_{L}=1$ and addressed the importance of the filter shape in the application of the similarity model. A top-hat or box filter, Gaussian filter and a spectral cutoff filter were implemented and the results were compared in the study. Comparisons of the SGS energy flux were also made. Testing showed that the model could produce good correlations between predicted and measured stresses when the box filter and the Gaussian filter were implemented. Much weaker correlations were reported when a sharp spectral cutoff was used. Conclusions from this testing were that the similarity model, when implemented with a box filter or Gaussian filter, showed high correlations with actual stresses and captured the characteristics of backscatter as well.

Equation (1.3) is an exceedingly simple expression that amounts to approximating the full velocity in the exact expression for the SGS stress (1.2) with the known in LES, resolved velocity. Therefore, in that framework there is no need to develop different functional forms of the SGS models. There is no even need to refer to physics of turbulence since (1.3) is just a mathematical relation involving velocity $\bar {u}_{i}$, known in LES, and a prescribed filter. Consequently, for a given filter the LES equations with the similarity stress tensor are autonomous in the same sense as Navier–Stokes equations, where knowledge of velocity field $\bar {u}_{i}({\boldsymbol {x}},t_0)$ at time $t_0$, together with the boundary conditions and the molecular viscosity, is sufficient to obtain the solution at times $t>t_0$. That simplicity has been one of the reasons for substantial interest in this type of model. However, Bardina et al. (1983), as well as Liu et al. (1994), observed that the model in its proposed form was not sufficiently dissipative in a posteriori tests. This usually leads to similarity-type models failing in actual LES, especially when the viscous dissipation is small.

For this reason mixed models were proposed which combine the similarity model with an eddy-viscosity expression to benefit from the positive attributes of each class of model, namely the dissipation of eddy-viscosity models along with the similarity model's correlations with actual SGS stresses. Such models were proposed and implemented by Bardina et al. (1983), Vreman, Geurts & Kuerten (1994, 1996) and Winckelmans et al. (2001) among others. The need for an additional eddy-viscosity term underscores the failure of the pure similarity model. However, recently Johnson (2022) developed a new methodology for LES (physics-inspired coarsening) which provides a mixed model expression as a natural outcome of the method. It is a more refined view of the mixed model than a combination of two seemingly unrelated expressions in the original mixed models.

The purpose of this work is to show that while the similarity-type models are not suitable for use directly in LES, they contain sufficient information for designing an autonomous, dissipative, eddy-viscosity model. The eddy viscosity developed here is autonomous in a sense that its functional form is not postulated but is obtained from a simulated LES field itself at each time step in simulations. Because it is derived from the similarity model expression, the respective SGS dissipations are also highly correlated. It should be noted that, in general, similarity models are developed without any explicit connection to eddy-viscosity models. One exception is work of Brun, Friedrich & da Silva (2006) who used the structure function model of Metais & Lesieur (1992), derived from a spectral eddy viscosity and reformulated it as the scale similarity-type model.

The current research is guided by recent work (Domaradzki 2021a,b, 2022) that developed an autonomous eddy-viscosity model in spectral space based on the physics of interscale energy transfer among resolved scales. One of the goals here is to extend the methodology to the physical space representation because it is most widely used in LES practice where various finite-difference/volume codes are used. However, spectral (Fourier) representation benefits from a direct link to the wealth of physical results offered by theories of turbulence. This provides a very firm basis for turbulence modelling compared with formulations in the physical space using graded filters. This is because there is a wide variety of graded filters to choose from and, additionally, for the Cartesian grids, three-dimensional (3-D) filters are often written as a product of one-dimensional (1-D) filters to simplify a numerical implementation. This causes LES results to be somewhat ambiguous because their dependence on the filter type and the implementation, e.g. a 3-D filter or a sequence of 1-D filters, cannot be excluded.

Analytical theories of isotropic turbulence as originated by Kraichnan's direct interaction approximation (Kraichnan 1959) provide closure expressions for the energy transfer term $T(k)$ in the spectral kinetic energy equation in terms of the energy spectrum $E(k)$. Kraichnan (1976) was the first to employ such spectral closures for SGS modelling. A good overview of analytical theories of turbulence and their use in modelling can be found in Lesieur (1997), Lesieur, Metais & Comte (2005) and Zhou (2021). Compared with traditional SGS models that are based on more-or-less phenomenological arguments to develop functional forms of the modelling expressions, the distinguishing feature of the spectral SGS models is that the primary quantity used in modelling is the SGS energy transfer $T_{SGS}(k\,|\,k_c)$, which is presumed to be known accurately from the underlying turbulence theory. The notation $T_{SGS}(k\,|\,k_c)$ indicates the energy transfer from a range of resolved scales $k \leq k_c$ caused by nonlinear interactions involving SGSs $k > k_c$, where $k_c$ is a cutoff wavenumber of a sharp spectral filter. In spectral SGS models the computed transfer is used to derive an expression for the spectral eddy viscosity $\nu _{eddy}(k\,|\,k_c)$ by normalising $T_{SGS}(k\,|\,k_c)$ by $2 k^2 E(k)$ for the range of resolved wavenumbers $k < k_c$. The modelling is completed by replacing the unknown SGS term in spectral LES equations for the cutoff $k_c$ with a term $-\nu _{eddy}(k\,|\,k_c)k^2 u_n(\boldsymbol {k})$, linear in the velocity $u_n(\boldsymbol {k})$ in the spectral representation. This term has a dissipative character as long as $\nu _{eddy}(k\,|\,k_c)$ is a positive quantity. Kraichnan (1976) used a particular analytical theory, the test field model (TFM), whereas Chollet & Lesieur (1981) used another formulation, the eddy-damped quasi-normal Markovian (EDQNM) approximation, to compute the spectral viscosity for the infinite inertial range spectrum $E(k) \sim k^{-5/3}$. In both cases the computed eddy viscosity is positive definite and has a relatively simple form with a constant plateau for wavenumbers $k$ less than approximately $0.4 k_c$ and rising in a form of a cusp to the maximum value at $k=k_c$. In actual LES performed with such models at high Reynolds numbers the inertial range spectrum is correctly recovered. The main advantage of such an approach to SGS modelling is that the functional form of the eddy viscosity is not postulated but is derived from the SGS energy transfer $T_{SGS}(k\,|\,k_c)$ obtained from theory.

It is also known that the SGS energy transfer as well as its wavenumber distribution can be obtained from DNS of isotropic turbulence, as shown for the first time by Domaradzki et al. (1987). The method introduced in that paper was later extended and used in numerous investigations as a diagnostic tool to elucidate and understand physics of nonlinear interactions acting in isotropic and wall-bounded turbulent flows, simulated using DNS and LES (Domaradzki & Rogallo 1990; Zhou 1993; Brasseur & Wei 1994). More recently, it was shown by Domaradzki (2021a,b, 2022) that such a numerical analysis method, in combination with the spectral eddy viscosity methodology of Kraichnan (1976), can be used to develop a self-contained SGS model without reference to explicit expressions of the analytical theories or any other classical SGS models. Specifically, the SGS energy transfer among resolved scales and its wavenumber distribution is computed from LES fields at each time step and cast in the form of a spectral eddy viscosity. Such a computed eddy viscosity is then modified to make it consistent with two known asymptotic properties of energy flux in the inertial range and used in the eddy-viscosity term added to the Navier–Stokes spectral solver as a SGS modelling term. Effectively, the procedure allows self-contained LES without use of extraneous SGS models or, equivalently, at each time step the model is obtained from a simulated field itself and asymptotic properties of the energy flux in the inertial range. The method has been tested in LES of isotropic turbulence at high Reynolds numbers and at lower Reynolds numbers decaying turbulence under conditions of the classical Comte-Bellot & Corrsin (1971) experiments. In both cases the agreement with reference data was excellent, providing numerical justification for the proposed approach.

The purpose of the current work is to extend this methodology to the physical space representation. While the exact, formal equivalence between modelling in the spectral and the physical representation does not exist, the model development in the physical space can productively follow the steps used in developing the model in the spectral space. Those steps are described and contrasted for both representations in the next section.

2. Description of the method

To contrast the physical and spectral space representation for the purpose of model development we summarise in the following the relevant equations and the corresponding goals of LES in both representations. The Navier–Stokes equations in the physical space for incompressible flow in Cartesian coordinates are

(2.1)\begin{gather} \frac{\partial {u}_n}{\partial x_n} = 0, \end{gather}

(2.2)\begin{gather}\frac{\partial {u}_n}{\partial t} ={-} \frac{\partial}{\partial x_j}({u}_n {u}_j)- \frac{\partial {p}}{\partial x_n} + \nu \frac{\partial^2 {u}_n}{\partial x_j \partial x_j}, \end{gather}

where ${\boldsymbol {u}}({\boldsymbol {x}},t) = (u_1,u_2,u_3)$ is the velocity field, and $x_1,x_2,x_3$ represent coordinates, respectively. The constant density $\rho$ is incorporated in pressure $p$.

2.1. The spectral space representation of Navier–Stokes equations

In this work we consider homogeneous, isotropic turbulence in a triply periodic box with a linear dimension $L$. Such conditions allow a straightforward implementation of numerical pseudo-spectral (Fourier) methods where the physical and the Fourier representation are used interchangeably to maximise the numerical accuracy and efficiency. Specifically, the domain is discretised using $N$ uniformly spaced grid points in each direction resulting in a mesh size $\Delta x = L/N$ and a total of $N^3$ grid points. Independent variables, here ${\boldsymbol {u}}({\boldsymbol {x}},t)$, are transformed between physical and spectral space using the discrete Fourier transform

(2.3)\begin{equation} \boldsymbol{u}(\boldsymbol{k}) = \frac{1}{N^3}\sum_{\boldsymbol{x}}\boldsymbol{u}(\boldsymbol{x}) \exp(-{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}), \end{equation}

and the inverse transform

(2.4)\begin{equation} \boldsymbol{u}(\boldsymbol{x}) =\sum_{\boldsymbol{k}}\boldsymbol{u}(\boldsymbol{k}) \exp({\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}), \end{equation}

where $\boldsymbol {x}$ are the mesh points in physical space and $\boldsymbol {k}$ are the discrete wavenumbers with components $k_i = \pm n_i \Delta k$, $n_i = 0,1,2,\ldots, N/2$, $i = 1,2,3$, and $\Delta {k} = 2{\rm \pi} /L = 1$. Equations (2.1)–(2.2) transformed into spectral (Fourier) space have the following form (see, e.g., Lesieur 1997 and Pope 2000):

(2.5)\begin{gather} {\rm i}k_{n}u_{n}(\boldsymbol{k},t)=0, \end{gather}

(2.6)\begin{gather}\frac{\partial u_{n}(\boldsymbol{k},t)}{\partial t} = N_n(\boldsymbol{k},t) - {\rm i}k_n p(\boldsymbol{k},t) - \nu k^{2} u_{n}(\boldsymbol{k},t), \end{gather}

where wavenumbers $\boldsymbol {k}$ are associated with scales of turbulent motions and $N_n$ is the Fourier transform of the nonlinear term

(2.7)\begin{equation} N_n(\boldsymbol{k},t) ={-}{\rm i}k_j \int {\rm d}\boldsymbol{p} u_{j}(\kern0.7pt \boldsymbol{p},t)u_{n}(\boldsymbol{k}-\boldsymbol{p},t). \end{equation}

Note that the same notation is used for dependent variables, with the physical space indicated by an independent variable ${\boldsymbol {x}}$ and the spectral space indicated by ${\boldsymbol {k}}$. The equation for the energy amplitudes $\frac {1}{2}|u(\boldsymbol {k},t)|^{2}= \frac {1}{2}u_{n}(\boldsymbol {k},t)u_{n}^{*}(\boldsymbol {k},t)$, where the asterisk denotes a complex conjugate, follows from (2.6)

(2.8)\begin{equation} \frac{\partial}{\partial t} \frac{1}{2}|u(\boldsymbol{k},t)|^{2} = T(\boldsymbol{k},t)- 2\nu k^{2}\frac{1}{2}|u(\boldsymbol{k},t)|^{2}, \end{equation}

where $T(\boldsymbol {k},t)$ is the nonlinear energy transfer

(2.9)\begin{equation} T(\boldsymbol{k},t)={\rm Real}(u_{n}^{*}(\boldsymbol{k}) {N_n}(\boldsymbol{k},t)), \end{equation}

and ${\rm Real}$ is the real part of a complex expression. Physical quantities of interest for isotropic turbulence are described in terms of the scalar wavenumber $k = |\boldsymbol {k}|$ by averaging over thin spherical shells defined for an arbitrary quantity $f(\boldsymbol {k})$ as

(2.10)\begin{equation} \langle\, f(\boldsymbol{k})\rangle=\frac{1}{N_{\boldsymbol{k}}}\sum_{\boldsymbol{k}}f(\boldsymbol{k}), \end{equation}

where $\langle \cdots \rangle$ denotes the shell average and the summation extends over all ${N_{\boldsymbol {k}}}$ modes in the shell of thickness $\Delta k$ centred at $k = |\boldsymbol {k}|$. For instance, the energy spectrum is defined as follows:

(2.11)\begin{equation} E(k,t)=4 {\rm \pi}k^2 \langle \tfrac{1}{2}u_n(\boldsymbol{k},t)u_n^*(\boldsymbol{k},t)\rangle. \end{equation}

2.2. SGS modelling in the spectral space representation

Because of the prescribed geometry the largest turbulent length scale is on the order $O(L)$ but there is no limit on the size of small scales for continuous velocity field $\boldsymbol {u}(\boldsymbol {x})$. Physics of isotropic turbulence, however, allows to introduce the smallest length scale, the Kolmogoroff length $\eta$, such that influence of scales smaller than $\eta$ on scales larger than $\eta$ is negligible for most practical purposes. This is the basis for DNS of isotropic turbulence where (2.2) are solved numerically in a box with a side $L$ on a mesh with the grid size $\Delta x = O(\eta )$. It is believed that simulations performed with a high-order numerical method on a grid satisfying the above condition will accurately represent turbulence dynamics, in particular the energy decay rates and the velocity correlations, represented either in a form of the structure functions or the energy spectra (Pope 2000). Equivalent DNS in the spectral space using (2.6) involves discretisation of wavenumbers with the smallest wavenumber $\Delta k=2{\rm \pi} /L$ and the largest wavenumber $k_{max} = O(1/\eta )$.

In LES one attempts to obtain turbulence properties for the same physical flow but using mesh much coarser than for well-resolved DNS, i.e. with $\Delta x \gg \eta$ or, equivalently, $k_{max} \ll (1/\eta )$. It is convenient to think about a LES velocity field as a projection of a high-dimensional DNS field on a lower-dimensional space, equivalent to sampling of a DNS field on a coarse mesh in the physical space, or to truncation of modes with high wavenumbers in the spectral space. Such a projection will be denoted by a ‘less than’ symbol, e.g. $\boldsymbol {u}^<(\boldsymbol {x})$ or $\boldsymbol {u}^<(\boldsymbol {k})$. The goal of LES is to obtain the velocity $\boldsymbol {u}^<$ that would have the same statistics as the actual DNS field projected/truncated to the LES subspace. The most straightforward implementation of such a program is in the spectral space using a sharp spectral filter with the cutoff $k_c \ll O(1/\eta )$, which removes all modes with $k > k_c$ but leaves all modes with $k < k_c$ unaffected. The resulting LES equations have the standard LES form in the physical space

(2.12)\begin{equation} \frac{\partial }{\partial x_n} u^<_n(\boldsymbol{x},t) = 0 , \end{equation}

(2.13)\begin{align} \frac{\partial }{\partial t} u^<_n(\boldsymbol{x},t) &={-} \frac{\partial }{\partial x_j}(u^<_n(\boldsymbol{x},t)\;u^<_j(\boldsymbol{x},t))^<- \frac{\partial} {\partial x_n} {p^<(\boldsymbol{x},t)}\nonumber\\ &\quad +\nu \frac{\partial^2}{\partial x_j \partial x_j} u^<_n(\boldsymbol{x},t) - \frac{\partial } {\partial x_j} \tau^<_{nj}(\boldsymbol{x},t), \end{align}

where $\tau ^<_{nj}$ is the SGS stress tensor

(2.14)\begin{equation} \tau^<_{nj}(\boldsymbol{x},t)=(u_n(\boldsymbol{x},t) u_j(\boldsymbol{x},t))^<- (u^<_n(\boldsymbol{x},t) \; u^<_j(\boldsymbol{x},t))^<. \end{equation}

The above equations are usually written using overbar to denote truncation and the SGS stress is written simply as $\tau _{nj}$ but that notation is reserved for later use to denote spatial filtering, see (2.31)–(2.33). The use of spectral truncations for computing SGS stress $\tau ^<_{nj}$ is indicated by an explicit superscript <. In addition, instead of index $n$ usually $i$ is used; we have changed the notation to avoid conflict with the imaginary unit $i$ appearing in the Fourier representation.

In spectral representation (2.5)–(2.6) lead to LES equations

(2.15)\begin{gather} {\rm i}k_{n}u^<_n(\boldsymbol{k},t)=0, \end{gather}

(2.16)\begin{gather}\frac{\partial u^<_n(\boldsymbol{k},t) }{\partial t} = N^<_n(\boldsymbol{k},t)- {\rm i}k_n p^<(\boldsymbol{k},t) - \nu k^{2} u^<_n(\boldsymbol{k},t)+ N^{SGS}_n(\boldsymbol{k},t),\quad k \leq k_c, \end{gather}

where the SGS term is

(2.17)\begin{equation} N^{SGS}_n(\boldsymbol{k},t)= [(N_n(\boldsymbol{k},t)- {\rm i}k_n p(\boldsymbol{k},t)) - (N^<_n(\boldsymbol{k},t)- {\rm i}k_n p^<(\boldsymbol{k},t))], \quad k \leq k_c. \end{equation}

Note that in (2.16) and (2.17) truncation symbol < applied to the nonlinear and pressure terms implies that these terms are computed using only resolved modes $u^<_n(\boldsymbol {k})$ and the result is further truncated to modes with $k \leq k_c$, indicated by notation $k\leq k_c$ following each equation. Terms $N_n$ and $p$ are computed using all dynamically relevant modes (resolved and subgrid modes) and the results are then truncated to the resolved modes. Following analytical theories of turbulence the SGS term in (2.16) is represented in the same form as the viscous term

(2.18)\begin{align} \frac{\partial u^<_n(\boldsymbol{k},t) }{\partial t} = N^<_n(\boldsymbol{k},t)- {\rm i}k_n p^<(\boldsymbol{k},t) - \nu k^{2} u^<_n(\boldsymbol{k},t) - \nu_{eddy}(k,t)k^2u^<_n(\boldsymbol{k},t),\quad k \leq k_c, \end{align}

with the eddy viscosity to be determined. Details of the method to determine $\nu _{eddy}(k,t)$ are described in Domaradzki (2021a,b, 2022). In this section we summarise main features of the method for the purpose of modifying it for LES in the physical space.

The starting point is the spectral LES energy equation for scales $k\leq k_c$, obtained from (2.16) by following derivation of (2.8). The evolution equation for energy of individual modes $(1/2)|u^<(\boldsymbol {k},t)|^{2}=(1/2)u^<_n(\boldsymbol {k})u^<_n(\boldsymbol {k})^*$ is

(2.19)\begin{equation} \frac{\partial }{\partial t} \frac{1}{2}|u^<(\boldsymbol{k},t)|^{2} = T^<(\boldsymbol{k},t)- 2\nu k^{2}\frac{1}{2}|u^<(\boldsymbol{k},t)|^{2}+T_{SGS}(\boldsymbol{k},t), \end{equation}

where $T_{SGS}(\boldsymbol {k},t)$ is the SGS transfer for an individual wavenumber mode $\boldsymbol {k}$. Because of isotropy the energy equation can be averaged over wavenumber shells centred at $k=|\boldsymbol {k}|$, giving

(2.20)\begin{equation} \frac{\partial}{\partial t} E^<(k\,|\,k_c)=T^<(k\,|\,k_c)-2 \nu k^2E^<(k\,|\,k_c)+T_{SGS}(k\,|\,k_c), \quad k \leq k_c, \end{equation}

where $E^<(k\,|\,k_c)$ is the energy spectrum of the resolved modes $k< k_c$, $T^<(k\,|\,k_c)$ is the energy transfer among resolved modes, and the SGS energy transfer term is

(2.21)\begin{equation} T_{SGS}(k\,|\,k_c)=T(k)-T^<(k\,|\,k_c), \quad k \leq k_c, \end{equation}

where $T(k)$ is the full nonlinear energy transfer computed using all modes, resolved and SGS. Total SGS energy transfer across $k_c$ is

(2.22)\begin{equation} T_{SGS}(k_c)= \int_{0}^{k_c} T_{SGS}(k\,|\,k_c)\,{\rm d}k. \end{equation}

For simplicity the dependence of the above quantities on time variable is omitted.

Following Kraichnan (1976), the SGS spectral energy equation can be formally rewritten as

(2.23)\begin{equation} \frac{\partial}{\partial t} E^<(k\,|\,k_c)=T^<(k\,|\,k_c)-2 \nu k^2 E^<(k\,|\,k_c)- 2 \nu_{eddy}(k\,|\,k_c) k^2 E^<(k\,|\,k_c), \end{equation}

where the SGS energy transfer is expressed in the same functional form as the molecular dissipation term by introducing the theoretical, effective eddy viscosity

(2.24)\begin{equation} \nu_{eddy}(k\,|\,k_c) ={-}\frac{T_{SGS}(k\,|\,k_c)}{2k^2E^<(k\,|\,k_c)}. \end{equation}

This is the quantity sought to be used in solving (2.18). As an example, assuming infinite inertial range spectrum $k^{-5/3}$, theoretical formulae for $T_{SGS}(k\,|\,k_c)$ can be computed numerically (Kraichnan 1976; Chollet & Lesieur 1981; Lesieur 1997) and the eddy viscosity (2.24) is well fitted by the expression derived by Chollet (1984)

(2.25)\begin{equation} \nu_{eddy}(k\,|\,k_c)=\left(\frac{E(k_c)}{k_c}\right)^{1/2}{C_K}^{{-}3/2}(0.441 +15.2\, {\rm e}^{{-}3.03k_c/k}) \equiv C_m\,f_1(k\,|\,k_c). \end{equation}

Here $C_K$ is the Kolmogorov constant, taken usually as $1.4$, $E(k_c)$ is the energy spectrum at the cutoff $k_c$ and $f_1=0.441 +15.2 \exp (-3.03k_c/k)$ is a spectral model shape function dependent only on $k/k_c$, its form here specific to the Chollet–Lesieur model. The numerical level of the eddy viscosity is determined by the model coefficient $C_m$ which depends on time through $E(k_c)$. It is worth repeating here that the eddy viscosity is not postulated but derived from the primary physical quantity which is the energy transfer across a wavenumber cutoff $k_c$ between resolved scales ($k< k_c$) and SGSs ($k>k_c$).

It was shown previously (Domaradzki 2021a,b, 2022) that the task of modelling $T_{SGS}(k\,|\,k_c)$ can be approached with limited reliance on the analytical theories by exploiting information already contained in numerical LES fields being simulated. That approach splits the task of modelling $T_{SGS}(k\,|\,k_c)$ into finding the total SGS transfer/dissipation, integrated over $0< k< k_c$ and, separately, its distribution in wavenumbers $k$. Following (2.25) the assumed form of the eddy viscosity is $C_m\, f(k\,|\,k_c)$, where the coefficient $C_m$ and the shape function $f(k\,|\,k_c)$ are determined from the resolved LES data. We enumerate in the following steps in that approach to be used subsequently as a guidance for model development in the physical space.

1. (i) The eddy viscosity $\nu _{eddy}(k\,|\,k_c)$ in (2.18) is obtained from the SGS energy transfer averaged over shells, $T_{SGS}(k\,|\,k_c)$, i.e. it accounts for an aggregate effect of transfer on a set of many modes, not for an effect on individual modes described by $T_{SGS}(\boldsymbol {k}\,|\,k_c)$.

2. (ii) For a resolved LES field with the cutoff $k_c$, a sharp test filter with the cutoff $ak_c, \ a<1$, is introduced to analyse SGS energy transfer within resolved range of scales. While different values of $a$ can be considered (see Domaradzki 2022) $a=1/2$ is a recommended choice.

3. (iii) Using the test cutoff $(1/2)k_c$ the resolved, $k$-dependent, SGS transfer $T^{res}_{SGS}(k\,|\,\frac {1}{2}k_c)$ can be computed from LES data in the simulations.

4. (iv) Using $T^{res}_{SGS}(k\,|\,\frac {1}{2}k_c)$ the resolved eddy viscosity $\nu _{eddy}^{res} (k\,|\,\frac {1}{2}k_c)$ can be computed from the definition (2.24).

5. (v) The resolved eddy viscosity $\nu _{eddy}^{res} (k\,|\,\frac {1}{2}k_c)$ can be rescaled to the actual LES cutoff $k_c$ using the similarity variable $k/k_c$ (see (2.25)).

6. (vi) The plateau level for rescaled $\nu _{eddy}(k\,|\,k_c)$ can be determined from the theoretical limit $k \rightarrow 0$.

7. (vii) Using $T^{res}_{SGS}(k\,|\,\frac {1}{2}k_c)$ the total resolved transfer $T^{res}_{SGS}(\frac {1}{2}k_c)$ across test cutoff $(1/2)k_c$ can be computed using integration (2.22).

8. (viii) The value of the model coefficient $C_m$ can be determined from the condition that the eddy viscosity dissipation must be equal to the total SGS transfer $T_{SGS}(k_c)$, (2.26).

9. (ix) The total SGS transfer $T_{SGS}(k_c)$ across the actual LES cutoff $k_c$ is not known directly from the LES data but can be determined from the formula

   (2.26)\begin{equation} T_{SGS}(k_c)=\frac{1}{1-a^{4/3}}T_{SGS}^{res}(ak_c). \end{equation}
   For $a=1/2$ the constant factor $D=1/(1-a^{4/3})=1.66$ and $T_{SGS}^{res}(ak_c)$ is found through step (vii).

Formula (2.26) was derived in Domaradzki (2021b) using the ultraviolet scaling of the energy flux for the inertial range dynamic established by Kraichnan (1971a,b)

(2.27)\begin{equation} \varPi_{uv}(k_c\,|\,k) \sim (k_c/k)^{4/3}\varPi(k_c), \quad k \gg k_c, \end{equation}

where $\varPi (k_c)$ is the total energy flux across $k_c$ (equal to the total SGS transfer $T_{SGS}(k_c)$) and $\varPi _{uv}(k_c\,|\,k), \; k>k_c,$ measures the amount of energy flux across $k_c$ caused by interactions involving at least one wavenumber mode with a wavenumber greater than $k$.

2.3. SGS modelling in the physical space representation

The procedure summarised previously, using sharp spectral filters, is a natural choice in LES using Fourier modes. Moreover, such a spectral formalism is strongly favoured by theories of isotropic turbulence, providing a wealth of well-established, unique concepts (e.g. the spectral energy flux) and results (e.g. the inertial range spectral form) that provide a conceptual guidance and unambiguous benchmarks for DNS and LES. However, in LES performed in the physical space representation with finite-volume or finite-difference numerical codes, sharp spectral filters are not feasible and graded filters are normally employed. Although it is natural to ask if and how the modelling procedure in the spectral representation can be extended to a more typical LES framework in the physical space representation, the answer is not obvious. Some preliminary results were discussed in Domaradzki (2021a) where in order to address this question Gaussian and box filters were used; see, e.g., Pope (2000). Specifically, the 1-D Gaussian filter kernel is

(2.28)\begin{equation} G(x,x^\prime)=\sqrt{\frac{6}{{\rm \pi} \Delta^2}} \exp\left( -\frac{6 |x-x^\prime|^2}{\Delta^2} \right) \end{equation}

and the 1-D box filter kernel

(2.29)\begin{equation} G(x,x^\prime)= \left\{\begin{array}{ll} 1/\varDelta & \text{if } {|x-x^\prime|} \leq \varDelta/2, \\ 0 & \text{otherwise}, \end{array}\right. \end{equation}

where in both cases the filter width was set to $\varDelta = 2 \Delta x$. The filtering operation in Cartesian coordinates for an arbitrary function $f(x,y,z)$ is then given by the formula

(2.30)\begin{equation} \bar{f}(x,y,z)=\int G(x,x^\prime) G(y,y^\prime) G(z,z^\prime) f(x^\prime,y^\prime,z^\prime) \,{{\rm d}\kern0.7pt x}^\prime \,{{\rm d}y}^\prime \,{\rm d}z^\prime. \end{equation}

When filtering, denoted by an overbar, is applied to full Navier–Stokes (2.2), the LES equations in a standard notation for an incompressible flow are

(2.31)\begin{gather} \frac{\partial }{\partial x_i}\bar{u}_i = 0, \end{gather}

(2.32)\begin{gather}\frac{\partial }{\partial t} \bar{u}_i + \frac{\partial }{\partial x_j}\bar{u}_i\;\bar{u}_j ={-} \frac{\partial} {\partial x_i} \bar{p} +\nu \frac{\partial^2}{\partial x_j \partial x_j} \bar{u}_i - \frac{\partial } {\partial x_j} \tau_{ij}, \end{gather}

where $u_i=(u_1,u_2,u_3)=(u,v,w)$, $p$ and $\nu$ are the velocity, pressure and the kinematic viscosity, respectively, and $\tau _{ij}$ is the SGS stress tensor

(2.33)\begin{equation} \tau_{ij}=\overline{{u}_i {u}_j} - \bar{u}_i \bar{u}_j. \end{equation}

The form of (2.31) and (2.32) requires that the filtering and differentiation commute (Ghosal & Moin 1995; Vasilyev, Lund & Moin 1998). In practice, however, these equations are frequently the starting point in SGS modelling without regard to formal requirements for their derivation in the filtering framework.

There is a notable difference in interpretation of quantities computed using the sharp spectral filter and graded physical space filters. For both cases all quantities in LES are represented on a mesh with a finite mesh size $\Delta x$, implying spectral support with a cutoff wavenumber $k_c={\rm \pi} /{\Delta x}$. The explicit notation, e.g. for velocity $\boldsymbol {u}^<$, was used to signify such quantities (see (2.13)–(2.14)). In spectral space $\boldsymbol {u}^<(\boldsymbol {k})$ is obtained from a full velocity field $\boldsymbol {u}(\boldsymbol {k})$ by truncation with the sharp spectral filter with the cutoff $k_c \ll 1/\eta$. Its complement $\boldsymbol {u}^>(\boldsymbol {k})=\boldsymbol {u}(\boldsymbol {k})-\boldsymbol {u}^<(\boldsymbol {k})$ represents true SGSs that are not resolved by a coarse LES mesh. In LES equations in the physical space (2.31)–(2.33), a projection on the LES mesh is only implied. The velocity $\bar {\boldsymbol {u}}(\boldsymbol {x})$ is the full velocity $\boldsymbol {u}(\boldsymbol {x})$ first filtered to $\bar {\boldsymbol {u}}(\boldsymbol {x})$ and then projected on the LES mesh, i.e. for complete clarity it should be really denoted as $\bar {\boldsymbol {u}}(\boldsymbol {x})^<$. When a sharp spectral filter is applied to the full velocity $\boldsymbol {u}(\boldsymbol {x})$ at $k_c$, the resulting filtered quantities lose completely information about scales $k>k_c$; those scales are truly unknown, SGSs $\boldsymbol {u}^>(\boldsymbol {x})$ in the physical space representation. In case of the graded physical space filters a filtered field $\bar {\boldsymbol {u}}(\boldsymbol {x})$ still has the same spectral support as the unfiltered field $\boldsymbol {u}(\boldsymbol {x})$, i.e. all scales $k>0$ are modified by filtering but none is removed. In such a case the velocity field $u^\prime _i=u_i - \bar {u}_i$ is often called a subfilter velocity. In particular, for invertible filters there is no information loss caused by filtering because all scales in $\boldsymbol {u}(\boldsymbol {x})$ can be recovered from $\bar {\boldsymbol {u}}(\boldsymbol {x})$ by a filter inverse. A good discussion of how differences between sharp spectral filter and graded physical space filters affect interpretation of LES results is given by Langford & Moser (1999), Domaradzki, Loh & Yee (2002) and Domaradzki & Adams (2002). In particular, (2.32) with the explicit notation for LES truncation is obtained directly by filtering of (2.13), indicated by an overbar

(2.34)\begin{align} \frac{\partial }{\partial t} \bar{u}^<_i(\boldsymbol{x},t)& ={-} \frac{\partial }{\partial x_j}(\bar{u}^<_i(\boldsymbol{x},t)\;\bar{u}^<_j(\boldsymbol{x},t))^<- \frac{\partial} {\partial x_i} {\bar{p}^<(\boldsymbol{x},t)} \nonumber\\ &\quad +\nu \frac{\partial^2}{\partial x_j \partial x_j} \bar{u}^<_i(\boldsymbol{x},t) - \frac{\partial}{\partial x_j} \tau^{full}_{ij}(\boldsymbol{x},t), \end{align}

where the full SGS stress (2.33) is explicitly written as

(2.35)\begin{equation} \tau^{full}_{ij}=(\overline{{u}_i {u}_j})^<- (\bar{u}^<_i \; \bar{u}^<_j)^<= \overline{{u}_i {u}_j} - \bar{u}^<_i \bar{u}^<_j . \end{equation}

The second equality in (2.35) applies if one considers only numerical LES equations and velocity fields, when the most outside superscripts <, implying projection on the LES mesh, may be ignored in the notation. In addition, if the full velocity decomposition into resolved scales and SGSs is made explicit, $\boldsymbol {u}(\boldsymbol {x})={\boldsymbol {u}}(\boldsymbol {x})^< + \boldsymbol {u}^>(\boldsymbol {x})$, the full SGS stress tensor can be split into two components, $\tau ^{full}_{ij}=\tau ^{res}_{ij}+ \tau ^{phy}_{ij}$ (see, e.g., Domaradzki et al. 2002). The resolved SGS stress tensor

(2.36)\begin{equation} \tau^{res}_{ij}=(\overline{{u}^<_i {u}^<_j})^<- (\bar{u}^<_i \bar{u}^<_j)^<= \overline{{u}^<_i {u}^<_j} - \bar{u}^<_i \bar{u}^<_j \end{equation}

is a SGS similarity-like stress, computed using the velocity projected on the LES mesh ${\boldsymbol {u}}(\boldsymbol {x})^<$. The second equality in (2.36) indicates that the most outside superscripts < may be ignored in the notation if only LES resolution is considered. The remaining term $\tau ^{phy}_{ij}$ has a form

(2.37)\begin{equation} \tau^{phy}_{ij}=(\overline{{u}^<_i {u}^>_j})^<+ (\overline{{u}^>_i {u}^<_j})^<+ (\overline{{u}^>_i {u}^>_j})^<, \end{equation}

which accounts for the physics of the nonlinear interactions that involve unknown SGSs ${\boldsymbol {u}}(\boldsymbol {x})^>$.

While differences between sharp spectral filters and physical space graded filters suggest that quantities computed with different filters in our modelling procedure could lead to incongruent outcomes, we were surprised to find that the energy transfer quantities computed in the physical space were useful in the spectral LES method. In Domaradzki (2021a), LES were performed using a slight modification of the spectral space procedure. The modification consisted of replacing the total resolved SGS transfer $T^{res}_{SGS}(\frac {1}{2}k_c)$ computed in the spectral space by the resolved SGS energy transfer computed in the physical space. The SGS energy transfer associated with the full SGS stress tensor (2.35) is

(2.38)\begin{equation} \epsilon_{SGS}(\boldsymbol{x})= \tau^{full}_{ij}(\boldsymbol{x}) \bar{S}^<_{ij}(\boldsymbol{x}), \end{equation}

where $\bar {S}^<_{ij}$ is the resolved rate-of-strain tensor

(2.39)\begin{equation} \bar{S}^<_{ij}= \frac{1}{2} \left( \frac{\partial \bar{u}^<_i}{\partial x_j} + \frac{\partial \bar{u}^<_j}{\partial x_i} \right). \end{equation}

The contribution to the total SGS dissipation provided by the resolved SGS stress in the physical space representation is

(2.40)\begin{equation} \epsilon^{res}_{SGS}(\boldsymbol{x})= \tau^{res}_{ij}(\boldsymbol{x}) \bar{S}^<_{ij}(\boldsymbol{x}). \end{equation}

Note that if LES are performed with only resolved SGS stress retained in (2.34), in general they will fail because they lack information contained in $\tau ^{phy}_{ij}$ about dynamically important actual SGSs for modes with $k>k_c$. This deficiency leads to insufficient SGS dissipation in actual LES performed with such pure similarity models.

In the above formulae the overbar denotes a spatial filtering procedure, e.g. using (2.28) or (2.29), and the space-dependent resolved SGS dissipation represents the resolved energy transfer component for the LES (2.40). To make a connection with the spectral SGS transfer the above formulae can also be implemented for a sharp spectral filter with cutoff $\frac {1}{2}k_c$ (i.e. the overbar implies sharp spectral filtering, the filtering and differentiation are performed in the spectral space, and multiplication of $\tau ^{res}_{ij}$ and $\bar {S}_{ij}$ in the physical space). In that case

(2.41)\begin{equation} \langle \epsilon^{res}_{SGS} \rangle=\frac{1}{N^3}\sum_{\boldsymbol{x}}\epsilon^{res}_{SGS}(\boldsymbol{x})=T^{res}_{SGS} \left(\frac{1}{2}k_c\right), \end{equation}

where summation is over all $N^3$ mesh points $\boldsymbol {x}$ and $T^{res}_{SGS}(\frac {1}{2}k_c)$ is the integrated resolved SGS transfer computed using (2.22) with the upper limit of integration $(1/2)k_c$. Note that while integrated SGS transfers in the physical and spectral space are the same, there is no equivalent spectral space formula corresponding to the local SGS dissipation in the physical space representation (2.40). Note also that for typical turbulent fields these integrated transfers are negative, consistent with the energy transfer from large to small scales in the mean.

In Domaradzki (2021a) the eddy-viscosity shape function $f(k\,|\,k_c)$ was prescribed and the model coefficient $C_m$ was computed according to the spectral implementation where the average value of (2.40), given by formula (2.41), served as a proxy for the resolved SGS transfer $T^{res}_{SGS}(\frac {1}{2}k_c)$. The results were in a close agreement with the corresponding results of LES performed entirely with the spectral filtering (see figure 9 in Domaradzki 2021a), indicating that at least some spectral quantities required in the enumerated steps for the modelling in the spectral space can be obtained using physical space quantities. This also suggests that the implementation of the procedure entirely in the physical space may be possible. Steps for such a procedure, based on the method in the spectral space, are described in the following.

2.4. The physical space eddy-viscosity model

The procedure seeks an expression for the eddy viscosity in the physical space $\nu _{eddy}(\boldsymbol {x},t)$ that is appropriate for advancing in time the unknown LES velocity field ${\boldsymbol {u}}^<(\boldsymbol {x},t)$ using LES (2.12) and (2.13). Note that the equations and the unknown velocity in the procedure are the same for the spectral and the physical space method. The only difference is how the eddy viscosity is determined. In spectral space the task is to solve LES (2.16) for the velocity $\boldsymbol {u}^<(\boldsymbol {k})$ in spectral representation. Using the eddy-viscosity model (2.16) is replaced by (2.18) and the spectral eddy viscosity $\nu _{eddy}(k)$ is determined through steps (i)–(ix) enumerated at the end of § 2.2. In physical space the task is to solve (2.13) for the velocity $\boldsymbol {u}^<(\boldsymbol {x})$ in real space representation. The eddy-viscosity model for the SGS stress tensor in (2.13) is

(2.42)\begin{equation} \tau^<_{ij}={-}2\nu_{eddy}(\boldsymbol{x}){S}^<_{ij}(\boldsymbol{x}), \end{equation}

where

(2.43)\begin{equation} {S}^<_{ij}= \frac{1}{2} \left( \frac{\partial {u}^<_i}{\partial x_j} + \frac{\partial {u}^<_j}{\partial x_i} \right). \end{equation}

Note that (2.42) ignores the usual term with the trace of the SGS stress because the procedure uses only energy equations where that term does not contribute.

The physical space eddy viscosity $\nu _{eddy}(\boldsymbol {x})$ is determined following steps used to determine the spectral eddy viscosity. In spectral space one uses the resolved LES field $\boldsymbol {u}^<(\boldsymbol {k})$ that is then filtered with a test sharp spectral filter with the cutoff $(1/2)k_c$ to compute the resolved SGS energy transfer $T^{res}_{SGS}(k\,|\,\frac {1}{2}k_c)$. The physical space equivalent is application of a graded filter to the resolved field $\boldsymbol {u}^<(\boldsymbol {x})$. Using graded filters the resolved energy transfer is computed using (2.40), which serves as a quantity equivalent in the physical space to $T^{res}_{SGS}(\boldsymbol {k}\,|\,\frac {1}{2}k_c)$ in the spectral space, i.e. the resolved SGS energy transfer for an individual wavenumber mode $\boldsymbol {k}$.

The eddy-viscosity closure model for the resolved SGS stress tensor (2.36) is

(2.44)\begin{equation} \tau^{mod}_{ij}={-}2\nu_{eddy}^{mod}(\boldsymbol{x})\bar{S}^<_{ij}(\boldsymbol{x}), \end{equation}

where $\bar {S}^<_{ij}$ is the resolved rate-of-strain tensor given by (2.39). Equation (2.44) implies a formal relation for $\nu _{eddy}^{mod}(\boldsymbol {x})$ that involves computed resolved SGS transfer

(2.45)\begin{equation} \epsilon^{res}_{SGS}(\boldsymbol{x})= \tau^{mod}_{ij}(\boldsymbol{x}) \bar{S}^<_{ij}(\boldsymbol{x})={-}2\nu_{eddy}^{mod}(\boldsymbol{x})\bar{S}^<_{ij}(\boldsymbol{x})\bar{S}^<_{ij}(\boldsymbol{x}). \end{equation}

It is well known that a naive application of this formula on a pointwise basis

(2.46)\begin{equation} \nu_{eddy}^{mod}(\boldsymbol{x})={-} \frac{\epsilon^{res}_{SGS}(\boldsymbol{x})}{2(\bar{S}^<_{ij}(\boldsymbol{x}))^2}, \end{equation}

will cause numerical instabilities when implemented in actual LES. This is because the SGS dissipation $\epsilon ^{res}_{SGS}(\boldsymbol {x})$ will contain spatial regions of forward and inverse transfer and the computed eddy viscosity will contain negative values that may become a source of instabilities. The key observation from the spectral procedure (item (i)) is that in computing the eddy viscosity the SGS transfer averaged over shells is used. While different types of averaging in the physical space can be considered, a simple filtering, the same as applied to derive (2.31)–(2.33), was found to be sufficient, i.e. the final expression for the eddy viscosity becomes

(2.47)\begin{equation} \nu_{eddy}(\boldsymbol{x})={-} \frac{\bar{\epsilon}^{res}_{SGS}(\boldsymbol{x})}{2(\bar{S}^<_{ij}(\boldsymbol{x}))^2}. \end{equation}

An important observation from the spectral procedure is that the eddy viscosity is computed first for the filtered field $k<\frac {1}{2}k_c$ ($\nu _{eddy}^{res} (k\,|\,\frac {1}{2}k_c)$). Subsequently, it is rescaled and applied in simulations to the full LES field $k< k_c$. We follow this sequence in the physical space as well. The eddy viscosity is computed first for the filtered field, signified by the resolved rate of strain in the denominator in (2.47). In the physical space the filtered field $\bar {u}^<_i$ already has the same spectral support as unfiltered (but truncated) LES field $u^<_i$ so no rescaling to the mesh cutoff is required. However, the eddy viscosity computed through (2.47) is applied in LES to the full, unfiltered LES field $u^<_i$, i.e. the SGS stress tensor in (2.14) is modelled through (2.42)–(2.43) with the eddy viscosity (2.47).

A simple, physical interpretation is that for two similar fields, here $\bar {u}^<_i$ and $u^<_i$, the eddy viscosities are similar, here assumed to be the same.

In summary, the operational prescription for applying the method is as follows.

• • Select a physical space filtering procedure, denoted by an overbar (e.g. here (2.28) or (2.29)).

• • Initialise a velocity $\boldsymbol {u}(\boldsymbol {x})^<$ on a LES mesh.

• • Compute resolved SGS stress $\tau ^{res}_{ij}$, (2.36), resolved SGS dissipation $\epsilon ^{res}_{SGS}(\boldsymbol {x})$, (2.40), and filtered resolved SGS dissipation $\bar {\epsilon }^{res}_{SGS}(\boldsymbol {x})$.

• • Compute eddy viscosity $\nu _{eddy}(\boldsymbol {x})$ from (2.47).

• • Form the model SGS stress tensor $\tau ^<_{ij}(\boldsymbol {x})$, (2.42), and use it to advance $\boldsymbol {u}(\boldsymbol {x})^<$ in time by solving LES (2.12)–(2.13).

3. Results

3.1. A fully autonomous method

To test the proposed implementation of the method in the physical space we have repeated LESs for several forced and decaying cases simulated previously using the spectral implementation of the method (Domaradzki 2021a,b, 2022). Details of the numerical method are provided in those papers and are reprised here for the sake of completeness.

The LES equations are solved using a pseudo-spectral numerical method of Rogallo (1981) in the implementation of Yeung & Pope (1988). For steady-state cases we use the forcing scheme of Sullivan, Mahalingam & Kerr (1994) in which the sum of squared amplitudes of velocity modes in a sphere of prescribed radius $K_f=3$ is kept constant in time. This is accomplished by multiplying all modes in the forced sphere by the same constant factor at the end of each time step to restore the energy in the sphere to the value at the beginning of the time step. Turbulence parameters were obtained from spectral quantities. The integral of $E(k)$ over $k$ gives turbulent kinetic energy per unit mass $\frac {3}{2}{u^\prime }^2$, where $u^\prime$ is the root-mean-square (r.m.s.) turbulent velocity. The integrated dissipation spectrum gives the dissipation rate of the turbulent kinetic energy, $\varepsilon$. The Taylor microscale is computed as $\lambda =(15{u^\prime }^2 \nu /\varepsilon )^{1/2}$ and the microscale Reynolds number is $Re_\lambda =u^\prime \lambda /\nu$. A time scale for the evolution of turbulence is the large-eddy turnover time $T_e = L_{int}/u^\prime$, where $L_{int}$ is the integral length scale

(3.1)\begin{equation} L_{int}=\frac{\rm \pi}{2{u^\prime}^2}\int_{0}^{\infty} k^{{-}1}E(k)\,{\rm d}k. \end{equation}

The macroscale Reynolds number is defined as $Re = u^\prime L_{int}/\nu$. Note that some parameters defined previously depend on the viscous dissipation $\varepsilon$ which is dominated by the high-wavenumber part of the energy spectrum, not available in LES. For high-Reynolds-number cases the viscous dissipation $\varepsilon$ is assumed to be equal to the energy flux $\varPi$ across the spectrum. In forced LES we estimated $\varPi =\varepsilon$ in a steady state as a difference between the measured energy input rate by forcing and the known viscous dissipation in the resolved range. The viscous dissipation in the resolved range was found to be four orders of the magnitude less than estimated $\varepsilon$. Using this estimate the initial Taylor microscale Reynolds number $Re_{\lambda }$ in all forced cases exceeded $10^4$, indicating that the inertial range theory should apply. Because of that, if the modelling procedure is correct, LES should recover known features of the inertial range dynamics. In what follows, we show that, indeed, the method develops and maintains in time the inertial range spectrum with a correct value of the Kolmogoroff constant outside the forcing wavenumber band.

The computational domain is a cubic box with side $L=2{\rm \pi}$, discretised using $N$ mesh points in each Cartesian direction. In the majority of cases $N=64$ and a sequence of 1-D Gaussians (2.28) was used for filtering. One case with $N=32$ was run to see the effects of resolution on the model proposed. In another case 1-D Gaussian filters were replaced by 1-D box filters (2.29) to see the effects of filtering on the model. The simulations were initialised with a velocity field consistent with a prescribed energy spectrum and random phases. That was followed by a short DNS run for 50–100 time steps (with only the molecular viscosity active) to build physically meaningful phase relationships. The fields at the end of the precursor DNS run were used as initial conditions for LES with the eddy viscosity turned on.

Parameters in the simulations for high-Reynolds-number cases are reproduced in table 1 for entries considered in this paper. In all cases simulations were run for an order of 10 large-eddy turnover times, reaching statistically steady state in less than half of the total run time. This is indicated by a close agreement between energy spectra averaged over the second half of each run ($n_{step}$ from 1000 to 2000) and the instantaneous spectrum at the end of a run ($n_{step}=2000$).

Table 1. Simulated LES cases. In all cases the molecular viscosity is $\nu =2.5 \times 10^{-7}$, the time step is $\Delta t = 0.005$, and the total simulation time is $T=10$. Turbulence parameters are shown at the initial time (IC) for two different initial spectra (eKolm and ePulse) and at the final time step in the simulations ($n_{step}=2000$).

In figure 1 we plot energy spectra obtained in simulations initialised with the $k^{-5/3}$ function with no prefactors. The spectral energy slopes at late times are in good agreement with the $-5/3$ exponent, with minor departures in the vicinity of the LES cutoff $k_c$. The compensated spectra in a form of a $k$-dependent Kolmogoroff function

(3.2)\begin{equation} C_K(k)=\frac{E(k)}{\varepsilon^{2/3}k^{{-}5/3}}, \end{equation}

fall within the expected range 1.4–2.1 outside the forcing wavenumbers.

Figure 1. Results for forced LESs initialised with the inertial range spectral form (eKolm). Lines with symbols $\circ$, initial conditions; broken line, spectrum after $N_t=2000$ time steps (around 10 large eddy turnover times); solid line, spectrum averaged over last $1000$ steps. In this and all subsequent figures thin straight lines show, as appropriate, $-5/3$ slope, and a boundary of the forcing band at $k=3$. For compensated spectra (b) horizontal lines mark expected range of values for the Kolmogoroff constant. (a) Energy spectra and (b) compensated energy spectra.

Spectra shown in figure 2 illustrate dependence of the method on the numerical resolution and a filter type. Despite significant decrease in the resolution from $64^3$ to $32^3$ mesh points for the case eKolm32 as compared with the case eKolm, there is no visible deterioration in the quality of the prediction of the inertial range. However, replacing the Gaussian filter by the tophat filter leads to some deterioration of the inertial range prediction for the case eKolmBox: high-wavenumber modes appear slightly overdamped whereas the intermediate-wavenumber modes have slightly higher energy level than for the case eKolm. We may speculate that the performance of the Gaussian filter is better because of its more non-local character, i.e. involving information from more mesh points around a given node than for the tophat filter that involves only information from the nearest neighbours.

Figure 2. Results for forced LESs initialised with the inertial range spectral form (eKolm IC) and run with: (a) numerical resolution of $32^3$ modes and a Gaussian filter; (b) numerical resolution of $64^3$ modes and a box filter. The meaning of different line types is the same as in figure 1(a). (a) Energy spectra for LES with $N=32$ and (b) energy spectra for LES with a tophat filter.

Simulations were repeated with a pulse initial condition for which $E(k) = 0$ for $k > 4$. The energy spectra and the compensated spectra are shown in figure 3. Note that the initial condition for run ePulse shown in figure 3(a) is obtained from a short precursor DNS, run for 50 time steps with only molecular viscosity active. The results lead to similar conclusions as for the case eKolm: a very good agreement with the inertial range predictions outside the forcing band and minor rise near $k_c$. It should be noted that for the pulse initial condition the steady state spectrum settles at a higher energy level. Despite that the Kolmogoroff function is very little changed compared with the case shown in figure 1, indicating that the method provides robust prediction of inertial range dynamics. In addition, the method as implemented does not rely on explicit or implicit assumptions about the inertial range exponent nor about the value of the Kolmogoroff constant. Both are predictions of the method.

Figure 3. Results for forced LESs initialised with the pulse initial condition (ePulse). The meaning of different line types is the same as in figure 1. (a) Energy spectra and (b) compensated energy spectra.

For testing the spectral method at lower Reynolds number in Domaradzki (2021b), decaying turbulence results of the classical experiments of Comte-Bellot & Corrsin (1971) were used. We repeated LES for the same test case using the physical space method. Comte-Bellot & Corrsin (1971) provided data for energy spectra measured at three times $tU_0/M=42$, 98 and 171, where $U_0=1000$ cm s$^{-1}$ is a free-stream velocity in the wind tunnel and $M=5.08$ cm is the grid size. In numerical tests the velocity field is initialised to be consistent with the energy spectrum at $tU_0/M=42$ and the task is to advance the fields in simulations and compare computed turbulence parameters and spectra with the experimental spectra at $tU_0/M = 98$ and 171. The experimental results are provided in the dimensional form using cgs units, with the lowest reported wavenumber $k_e=0.15$ cm$^{-1}$ and the highest $k_e=20$ cm$^{-1}$ (subscript $e$ indicates experimental wavenumber). Since in pseudo-spectral numerical simulations the minimum wavenumber is by default equal to unity, the experimental data are converted into simulation data by changing the unit of length from centimetres to a new unit $L$: $1 \ {\rm cm} = k_e^{min} L$, where $k_e^{min}$ is the numerical value of the minimum experimental wavenumber to be represented in the simulations. We extrapolated the Comte-Bellot and Corrsin data to the minimum wavenumber $k_e^{min}=0.1$ cm$^{-1}$, with the maximum wavenumber kept at $k_e^{max}=20$ cm$^{-1}$. In new units the simulation parameters are $\Delta k= 10 \Delta k_e [L^{-1}]=1 [L^{-1}]$, $\nu =0.1^2 \nu _e [{L^2}/{s}]=0.0015 [{L^2}/{s}]$ (for air, $\nu _e=0.15$ cm$^2$ s$^{-1}$), the energy spectrum $E=0.1^3 E_e [{L^3}/{s^2}]$ and time remains dimensional in seconds. All turbulent quantities can be converted in the same way for a comparison between simulation and experimental data.

In figure 4(a) we plot experimental energy and dissipation spectra from Comte-Bellot & Corrsin (1971). The thin vertical line indicates the cutoff wavenumber selected for simulations ($k_c=30$, equivalent to $k_e=3$ cm$^{-1}$ in cgs units). For all three time instants in the experiments the viscous dissipation peak is located outside resolved range. Because the numerical resolution is unable to properly capture the dissipation, attempts at DNS with such a resolution fail as seen in figure 4(b) for a no-model case. Time evolution of the energy spectra in LES obtained using the physical space procedure is shown in figure 5 and can be compared with corresponding results for the spectral procedure (shown in figure 7 in Domaradzki 2021b). Time evolution of the energy spectrum is predicted quite well for both time intervals $U_0 t/{M}=[42,98]$ and $U_0 t/{M}=[98,171]$.

Figure 4. (a) Experimental energy and viscous dissipation spectra at different times. Lines of the same type show energy spectra, $E(k)$ (with symbols), and dissipation spectra, $D(k)$ (without symbols), respectively. Solid lines and circles, $U_0 t/{M}=42$; broken lines and squares, $U_0 t/{M}=98$; dotted lines and triangles, $U_0 t/{M}=171$. (b) Time evolution of energy spectra in under resolved DNS. Markers correspond to experimental data and lines show progression in time from the initial energy spectrum at $U_0 t/{M}=42$ to the final time at $U_0 t/{M}=171$ (dotted line).

Figure 5. Time evolution of energy spectra in an autonomous LES run: (a) for time interval $U_0 t/{M}=[42,98]$; (b) run continued for time interval $U_0 t/{M}=[98,171]$.

Some turbulence quantities computed in LES are collected in table 2 and compared with experimental results. They are in a very good agreement with the same data obtained in Domaradzki (2021b) (see table III in that paper) using fully spectral eddy viscosity. That reference contains also more extensive discussion of the comparison between LES data and experiments. Overall, these results demonstrate convincingly that the proposed physical space method is able to predict time evolution of turbulence for this classical case.

Table 2. A comparison of turbulence parameters between LES and experiments (Comte-Bellot & Corrsin 1971). LES performed with $64^3$ resolution, initialised with the energy spectrum at $U_0 t/{M}=42$ shown in figure 4: $u^\prime$, r.m.s. turbulent velocity; $\varepsilon$, dissipation rate; $\eta$, Kolmogoroff microscale; $\lambda$, Taylor microscale; $L_{int}$, integral length scale; $R_\lambda$, microscale Reynolds number. Experimental values are shown in parentheses, after conversion using a length unit $[L]=0.1$ cm. In the first row the experimental data are at the initial time ${U_0 t}/{M}=42$ whereas LES data at a somewhat later time ${U_0 t}/{M}=51$ to allow build-up of nonlinear interactions from random phases in the initial condition.

It may be noted that the fundamental quantity in the method development in the physical space is the resolved SGS dissipation $\epsilon ^{res}_{SGS}(\boldsymbol {x})$ (2.40). The resolved SGS stress tensor $\tau ^{res}_{ij}$ (2.36) can be considered as a generalised similarity model in a sense that it is computed using the velocity $\boldsymbol {u}(\boldsymbol {x})^<$ with the same spectral support as used for LES equations, i.e. no true SGSs with $k>k_c$ enter into its computation. Alternatively, it can be considered as a deconvolution model, in a sense that the velocity $\boldsymbol {u}(\boldsymbol {x})^<$ can be recovered from the filtered velocity $\bar {\boldsymbol {u}}(\boldsymbol {x})^<$ by inversion of the filtering operation, which for a Gaussian filter can be performed exactly. The SGS stress tensor for similarity or deconvolution models is known to be highly correlated with the exact SGS stress computed from full velocity data containing SGSs $k>k_c$ (Liu et al. 1994; Meneveau & Katz 2000). Despite that, similarity and deconvolution models fail in actual LES because they cannot maintain adequate SGS energy dissipation as simulation time progresses. On the other hand, eddy-viscosity-based models show very low correlations with exact SGS quantities but perform well in actual LES because of their good dissipative properties. In contrast to such common observations, for the method proposed here we find very high correlations between the eddy-viscosity results and the results obtained with the similarity model. Specifically, in figure 6 we plot SGS energy transfer for both models. Note that forward transfer is signified by negative values, i.e. acting as an energy sink in the LES dynamics, and backscatter is signified by positive values, i.e. acting as an energy source in the dynamics. Visual inspection of colour plots indicates that the energy transfer for both cases appears quite correlated. The computed correlation coefficient for these planar data is 0.81. We have also computed a correlation coefficient for the full 3-D data, getting values also in excess of $0.8$. Note these are much higher values than the value around $0.4$ found for the standard Smagorinsky model. High correlations are not surprising because the eddy viscosity is derived from the SGS dissipation of the similarity model. Yet the presence of backscatter for the SGS dissipation of the eddy-viscosity model is surprising, as it is commonly believed that it would lead to unstable simulations. However, none of the cases simulated with this approach showed any hints of instability. The analysis of the computed fields showed that the total forward energy transfer was at least an order magnitude greater than the backscatter. We suspect that its overall dominance in the total energy transfer may explain why relatively small negative values of eddy viscosity, fluctuating in space and time, do not lead to catastrophic instabilities. It must be noted that the above conclusions were reached for graded filters, here specifically for the Gaussian filter. Correlations between the actual SGS stresses and similarity-type stresses computed using sharp spectral filters are known to be significantly smaller (Liu et al. 1994). In addition, for sharp spectral filtering the forward and inverse SGS transfers are of the same order of magnitude, with the magnitude of the net forward transfer being much smaller than each individual forward/inverse component as shown in, e.g., Piomelli et al. (1991) and Domaradzki, Liu & Brachet (1993).

Figure 6. SGS energy transfer in a cross-sectional plane, with x and y axis labels corresponding to mesh point numbers along x and y direction: (a) for the similarity model; (b) for the eddy-viscosity model (2.42)–(2.47). Positive values are signified by red and orange colours, and negative values by yellow, green and blue.

3.2. A method with the total energy flux constraint

Simulation results using the SGS stress (2.42) with the eddy viscosity (2.47) are very encouraging but the method does not implement fully all steps involved in the spectral space representation. In Domaradzki (2021a) it was found that the last step (ix), with (2.26) that predicts the total energy flux at the cutoff wavenumber $k_c$, was critical for the accuracy of spectral results. In particular, under prediction of the total transfer led to gradual but significant loss of accuracy, in a way reminiscent of the classical similarity model. One may thus speculate that the current procedure could be further improved if that last step is implemented.

In accordance with the form of the eddy viscosity in spectral space (see (2.25)) it would require introducing a model constant in the definition of the eddy viscosity (2.47), i.e.

(3.3)\begin{equation} \nu_{eddy}^{(C)}(\boldsymbol{x})={-} C\frac{\bar{\epsilon}^{res}_{SGS}(\boldsymbol{x})}{2(\bar{S}^<_{ij}(\boldsymbol{x}))^2}. \end{equation}

The role of the model constant $C$ is to enforce the global constraint for the total SGS energy transfer (items (viii) and (ix) for the spectral space method). It is similar to the model coefficient $C_m$ in the spectral space method. However, because of the postulated form of the spectral eddy viscosity (see (2.25)) where function $f_1(k\,|\,k_c)$ is non-dimensional, $C_m$ has dimensions of viscosity whereas $C$ in (3.3) is non-dimensional. For the spectral representation the model coefficient $C_m$ was determined uniquely from the ultraviolet locality scaling of the energy flux for a sharp spectral filter and the inertial range spectrum (Domaradzki 2022). However, there are no such scaling results available for arbitrary, graded filters and the model coefficient $C$ in (3.3) must be treated as an adjustable constant that depends on the filter employed and the energy spectrum of a simulated field. For spatial filtering in the present work we have used a product of 1-D Gaussian filters with $\varDelta =2\Delta x$, (2.28) and (2.30). For this specific filter it was possible to compare resolved SGS dissipation $\langle \epsilon ^{res}_{SGS} \rangle$ with the total transfer $T_{SGS}(k_c)$ obtained for the same velocity fields $\boldsymbol {u}(\boldsymbol {x})^<$ using the spectral formulae. In the previous papers (Domaradzki 2021b, 2022) the sharp spectral filter was used to analyse two distinct velocity fields: one consistent with the inertial range spectrum and another one consistent with the experimental spectrum from Comte-Bellot & Corrsin (1971). Inspection of the numerical data showed that the model constant is $C \approx 0.8$ for both flows, thus not too far from $C=1.0$ corresponding to eddy viscosity (2.47). While this might be a fortuitous agreement for the particular filter type with a particular filter width used in this work, the obtained results encourage further work to quantify the total SGS dissipation for different filter types.

The majority of simulations reported in § 3.1 have been repeated using eddy viscosity (3.3) with $C=0.8$ and are shown in figures 7–9. The results are quite acceptable. For high Reynolds numbers the spectral energy slopes at late times are in a good agreement with the $-5/3$ exponent, except in the vicinity of the LES cutoff $k_c$. The $k$-dependent Kolmogoroff functions fall within the expected range 1.4–2.1 outside the forcing wavenumbers. However, as the cutoff $k_c$ is approached a steep increase in $C_K$ is notable. Such a behaviour of spectra in the vicinity of $k_c$ is consistent with insufficient SGS dissipation in that range. It could be explained by a reduced total dissipation ($C=0.8$ vs $C=1.0$) but more likely it is caused by inadequacies in representing distribution of SGS transfer in scale (wavenumber). The spectral procedure imposes scale distribution through an eddy viscosity shape function that is composed of a plateau and a cusp (see, e.g., (2.25)). The presence of a cusp at $k_c$ in the spectral eddy viscosity considered in previous papers (Domaradzki 2021a,b, 2022) led to better agreement with the inertial range form near $k_c$, suggesting that the physical space eddy viscosity is less active near $k_c$ than the spectral eddy viscosity.

Figure 7. Results for forced LESs, using the eddy viscosity (3.3) with $C=0.8$, initialised with the inertial range spectral form. The meaning of different line types is the same as in figure 1. (a) Energy spectra and (b) compensated energy spectra.

Figure 8. Results for forced LESs, using the eddy viscosity (3.3) with $C=0.8$, initialised with the pulse initial condition. The meaning of different line types is the same as in figure 1. (a) Energy spectra and (b) compensated energy spectra.

Figure 9. Time evolution of energy spectra using the eddy viscosity (3.3) with $C=0.8$: (a) for time interval $U_0 t/{M}=[42,98]$; (b) run continued for time interval $U_0 t/{M}=[98,171]$.

Time evolution of the energy spectra in LES for the Comte-Bellot & Corrsin (1971) case, obtained using the physical space procedure with $C=0.8$ is shown in figure 9 and can be compared with corresponding results for the spectral procedure (shown in figure 7 in Domaradzki 2021b). Time evolution of the energy spectrum is predicted reasonably well for the first time interval $U_0 t/{M}=[42,98]$, though closer inspection of figure 4(b) shows that LES results at $U_0 t/{M}=98$ overpredict experimental data at higher $k$. Corresponding results obtained using the spectral procedure (figure 7 in Domaradzki 2021b) showed better match in that range of wavenumbers. These observations are consistent with conclusions reached for high-Reynolds-number cases: the physical space method produces less dissipation near mesh cutoff than the spectral method. For the second time interval $U_0 t/{M}=[98,171]$ the energy peak is under predicted at the final instant, similar to the spectral case (figure 7 in Domaradzki 2021b), but at time $U_0 t/{M}=171$ the LES curve matches experimental spectra as accurately as LES with $C=1$.

3.3. Comparisons with other SGS modelling procedures

The primary goal of the current effort was to determine if and how the spectral space methodology can be extended to the physical space representation. To judge the extent to which this program has been successful it is most appropriate to compare performance of the proposed procedure with performance of equivalent spectral eddy-viscosity models. In figure 10 a line representing averaged energy spectrum from figure 7 is reproduced (solid black line) where it can be compared with corresponding spectra obtained using three different spectral eddy viscosities for the same initial condition eKolm in table 1. A broken black line has been obtained using the classical Chollet–Lesieur eddy viscosity given by (2.25) derived from analytical theories of turbulence. A dotted line is for the eddy viscosity obtained from numerical LES data with the plateau level set to a constant fraction of the cusp value at $k/k_c=1$ (Domaradzki 2021b, 2022). Finally, the broken dotted line is obtained with a $k$-independent eddy viscosity, i.e. constant in the domain, computed at each time step by requiring that its energy dissipation rate is equal to the value given by (2.26). That last curve practically overlaps with the solid black line for the eddy-viscosity case (3.3) for $C=0.8$. It suggests that in practice the eddy viscosity (3.3) acts as a time-dependent but constant-in-space viscosity. A similar conclusion was reached by Thiry, Winckelmans & Duponcheel (2019) for the dynamic Smagorinsky model in LES at high Reynolds numbers: ‘The derived spectral SGS viscosity was confirmed to be essentially uniform over all wavenumbers’. Even earlier, it was observed in LES of stratified turbulence by Siegel & Domaradzki (1994) that spectra of a standard Smagorinsky eddy viscosity were approximately uniform in wavenumbers. Nevertheless, all four eddy-viscosity models considered in figure 10, despite individual differences, produce energy spectra in general agreement among themselves and with the inertial range spectrum. The differences are visible only in the immediate vicinity of cutoff $k_c$ and can be attributed to differences in the wavenumber dependence of the eddy viscosities. In particular, the eddy viscosity that combines a fixed plateau and a cusp obtained from LES data (dotted line) shows the best agreement with the inertial range spectrum, extending to the cutoff wavenumber.

Figure 10. Results for forced LESs initialised with the inertial range spectral form and different eddy-viscosity models: solid line, eddy viscosity (3.3) with $C=0.8$; broken line, the classical Chollet–Lesieur eddy viscosity, (2.25); broken-dotted line, wavenumber-independent eddy viscosity; dotted line, eddy viscosity with a numerically determined cusp and a constant plateau level $p=0.37$ of the cusp at $k_c$ (Domaradzki 2021b, 2022).

The form of the eddy viscosity with the multiplicative constant $C$ (3.3) may suggest that the method could be amenable to the dynamic procedure of Germano et al. (1991). However, the dynamic procedure is not useful for the present method as discussed in the previous paper (Domaradzki 2021a). This is because both methods set different objectives. In the current approach the objective is to find the best estimate of the SGS transfer using a resolved LES field. In particular, the constant $C$ is selected by requiring that the total SGS transfer computed using graded filters approximates the total SGS transfer for a field truncated using a sharp spectral filter with the cutoff wavenumber $k_c$. In the dynamic procedure the objective is to find the best estimate of model coefficients using a resolved LES field, i.e. a model is selected first and then the Germano identity is used to determine model constants. As noted in Domaradzki (2021a) this implies that the SGS transfer predicted by the dynamic procedure will depend not only on a given velocity field but also on a selected model; i.e. for the same velocity field but different models the dynamic procedure can give different values of the SGS dissipation. It was also shown in that paper that the actual SGS transfer at $k_c$ is at most by a factor $1.66$ greater than the resolved SGS transfer at the test cutoff $(1/2)k_c$, i.e. values of a SGS dissipation predicted using a SGS model cannot be substantially different from the resolved SGS transfer.

Although it was shown by Ghosal et al. (1995) that it is possible to formulate a general integral equation for space- and time-dependent model coefficients, in practice a limited number of coefficients is determined by the least squares method as proposed by Lilly (1992). In most cases just one coefficient is sought, in the current context $C$, and the simplest application of the dynamic procedure is based on the energy form used in the original paper by Germano et al. (1991)

(3.4)\begin{equation} \tau^{res}_{ij}(\boldsymbol{x}) \bar{S}^<_{ij}(\boldsymbol{x})= (\tau^{full}_{ij}(\boldsymbol{x}) - \overline{\tau^{<}_{ij}}(\boldsymbol{x}))\bar{S}^<_{ij}(\boldsymbol{x}). \end{equation}

In the above equation $\tau ^{res}_{ij}$ is the resolved SGS tensor (2.36), known in LES. Here $\overline {\tau ^{<}_{ij}}$ is the unknown SGS tensor (2.14) obtained using spectral truncations, and then filtered with a graded filter, here (2.28) and (2.30). Finally $\tau ^{full}_{ij}$ is the SGS tensor (2.35) obtained using the same graded filter, projected on LES mesh (truncated spectrally at $k_c$). Modelling terms on the right-hand side of (3.4) by expression (2.42) with the eddy viscosity (3.3) leads to the following equation:

(3.5)\begin{equation} \epsilon^{res}_{SGS}(\boldsymbol{x})={-}2(C(\boldsymbol{x})\nu_{eddy}(\boldsymbol{x})\bar{S}^<_{ij}(\boldsymbol{x})- \overline{C(\boldsymbol{x})\nu_{eddy}(\boldsymbol{x}){S}^<_{ij}(\boldsymbol{x})})\bar{S}^<_{ij}(\boldsymbol{x}), \end{equation}

with the form of the eddy viscosity $\nu _{eddy}$ given by (2.47).

If $C(\boldsymbol {x})$ can be taken out of the filtered expression on the right-hand side the calculation of $C$ at each point $\boldsymbol {x}$ involves division by a term that fluctuates around zero causing the entire expression to be ill-conditioned, as could have been anticipated on the basis of the original work of Germano et al. (1991). Indeed, when (3.5) was used for a particular velocity field in the simulations, $C(\boldsymbol {x})$ fluctuated wildly. Averaging over the computational domain produced a value $C=-3.9$, completely unrelated to the value $C \approx 0.8$ determined from the global SGS transfer constraint, and with the large standard deviation on the order $O(10^3)$. A potentially better approach is to consider (3.5) averaged over computational domain

(3.6)\begin{equation} \langle \epsilon^{res}_{SGS}(\boldsymbol{x})\rangle ={-}2(\langle C(\boldsymbol{x})\nu_{eddy}(\boldsymbol{x})\bar{S}^<_{ij}(\boldsymbol{x})\bar{S}^<_{ij}(\boldsymbol{x})\rangle - \langle \overline{C(\boldsymbol{x})\nu_{eddy}(\boldsymbol{x}){S}^<_{ij}(\boldsymbol{x})}\;\bar{S}^<_{ij}(\boldsymbol{x})\rangle), \end{equation}

where $\langle \cdots \rangle$ denotes averaging over all mesh points (see (2.41)). This form also helps to understand better the source of difficulties in using the dynamic procedure for the model proposed here. Assuming constant $C$ and using the definition of the eddy viscosity (2.47) gives

(3.7)\begin{equation} \langle \epsilon^{res}_{SGS}(\boldsymbol{x})\rangle = C\left(\langle \bar{\epsilon}^{res}_{SGS}(\boldsymbol{x})\rangle - \left\langle \overline{\frac{\bar{\epsilon}^{res}_{SGS}(\boldsymbol{x})}{(\bar{S}^<_{ij}(\boldsymbol{x}))^2}{S}^<_{ij}(\boldsymbol{x})} \;\bar{S}^<_{ij}(\boldsymbol{x})\right\rangle \right). \end{equation}

The form of individual terms with $\epsilon ^{res}_{SGS}(\boldsymbol {x})$ in the above equation suggests that their values are comparable, which was confirmed by a direct calculation. Therefore, in order to enforce that equality large values of $C$ are needed, leading in turn to excessive SGS dissipation. This is a direct consequence of the definition of the eddy viscosity in terms of the resolved SGS transfer $\epsilon ^{res}_{SGS}(\boldsymbol {x})$. We thus conclude that the dynamic procedure is unlikely to offer improvements to the method proposed here because its difficulties in predicting the actual, physical SGS energy transfer.

The ultimate goal in developing eddy viscosity in the physical space is for applications to inhomogeneous flows and in the near future we plan to report such results for turbulent channel flow at Reynolds numbers up to $Re_{\tau }=2000$. However, there are several valid reasons for choosing homogeneous isotropic turbulence and Fourier spectral methods as a first step in the method development. For any SGS model to be viable it must be able to predict properties of isotropic turbulence at high Reynolds numbers. The current model is indeed able to predict the Kolmogoroff inertial range form. Moreover, without any additional input, simulations predict the Kolmogoroff constant in an accepted range of experimental and numerical values. This should be contrasted with some classical SGS models for isotropic turbulence (e.g. Smagorinsky or Chollet–Lesieur) that must assume the inertial range form and a value of the Kolmogoroff constant in deriving the model expressions. SGS models should also demonstrate predictive capabilities for flows at lower Reynolds numbers. This is often a simpler test case because the presence of the molecular dissipation assists the model SGS dissipation in LES. The current method, without any adjustments to the method used for high-Reynolds-number cases, was successful for the test case based on the experimental data of Comte-Bellot and Corrsin. In addition, it is well known that many standard finite-difference and finite-volume methods are plagued by a numerical dissipation, which in LES may be comparable to the model SGS dissipation. The numerical dissipation for actual LES is difficult to quantify and its presence may cloud conclusions regarding tests of various SGS models. The use of a de-aliased, spectral code in the present work assures that the numerical dissipation is negligible, i.e. effects observed in LES can be attributed entirely to the modelling procedure. Finally, there is also fundamental value in changing the focus of modelling from the SGS stress tensor form to the energy transfer. This change of focus, originally introduced by Kraichnan (1976), has been largely ignored in the modelling community, because the traditional theoretical formalism is limited to the spectral space representation. Yet focusing on the energy transfer allows one to use a wealth of theoretical knowledge about isotropic turbulence for SGS modelling as shown in recent papers by the current author. The method presented here is an attempt to extend the spectral SGS modelling formalism based on the energy transfer to the physical space representation. The model development uses the dissipation of the similarity model as a fundamental energy transfer quantity and follows steps in the model development for the spectral space representation. Although there is no one-to-one correspondence between these two representations, the eddy-viscosity model in the physical space based on the dissipation of the similarity stress tensor appears to be a successful implementation of the proposed program, at least for isotropic turbulence for the time being.

4. Conclusions

A previously proposed SGS modelling procedure based on the interscale energy transfer among resolved scales in LES, described for the spectral space implementation in Domaradzki (2021a,b, 2022), has been extended to the physical space representation. As the fundamental quantity the method employs the SGS energy transfer computed using a similarity-type model expression for the SGS tensor obtained using Gaussian filtering of velocity fields advanced in the simulations. The method development has followed several steps, enumerated for the spectral space representation, which provide guidance for designing a spatially varying eddy viscosity at each time step in LES. The computed eddy viscosity has then been employed to model the SGS stress tensor in the familiar Boussinesq form for use in LES. The method is autonomous in a sense that the form of the eddy viscosity (2.47) is not postulated but is extracted from the LES data without any adjustable constants. However, the method is unlikely to be universal because it is expected that it must depend on the filter type and the filter width. This should be contrasted with the spectral method utilising the sharp spectral filter where all model constants can be determined uniquely from the LES data and the analytical theories of turbulence (Domaradzki 2021a,b, 2022). Therefore, the physical space method is generalised by introducing in (3.3) a model constant $C$ to account for its dependence on a specific filter used and its potential influence on the computed total SGS transfer. The method was tested in LES of isotropic turbulence at high Reynolds numbers where the inertial range dynamics is expected and for lower Reynolds number decaying turbulence under conditions of the classical Comte-Bellot and Corrsin experiments. First, the fully autonomous method was considered where the eddy viscosity (2.47) does not refer to any adjustable constants, i.e. the model constant is implicitly $C=1$. For both flows the agreement with reference data is very good and the SGS transfer computed for the proposed eddy-viscosity model is highly correlated with the transfer computed for the similarity stress tensor. Subsequently, $C$ was determined by comparing resolved energy transfer computed in the physical space with the corresponding transfer computed in the spectral method. Inspection of the numerical data showed that the model constant for the filter type chosen was $C \approx 0.8$ for both flows, not too far from the value $C=1$ implied by the fully autonomous method. Because of that the spectral results were quite comparable in both cases. However, a nominally higher value of $C$ in the autonomous method suggests overall larger SGS dissipation, explaining minor differences in spectral results observed in figures in this paper for both cases.

In summary, the current procedure, based on the SGS energy transfer of the similarity model, can produce the eddy-viscosity expressions in the physical space representation that are not only as globally dissipative as standard eddy-viscosity models, but also that they predict modelled SGS dissipation which approximates the SGS dissipation of the similarity model well and is highly correlated with it.