2. WMRA framework for SGS modelling

The present study employs a WMRA framework originally developed by Meneveau (Reference Meneveau1991) for studying the spectral energy transport and later extended to turbulent channel flows (Dunn & Morrison Reference Dunn and Morrison2003), turbulent combustion (Kim et al. Reference Kim, Bassenne, Towery, Hamlington, Poludnenko and Urzay2018), droplet-laden turbulence (Freund & Ferrante Reference Freund and Ferrante2019) and particle-laden turbulence (Nabavi et al. Reference Nabavi, Di Renzo and Kim2022). Compared with the Fourier analysis of spectral energy transfer, wavelet offers several distinct advantages (Meneveau Reference Meneveau1991; Kim et al. Reference Kim, Bassenne, Towery, Hamlington, Poludnenko and Urzay2018; Nabavi et al. Reference Nabavi, Di Renzo and Kim2022). With a sufficiently complex wavelet basis, a spatial localisation of the analysis at a good spectral resolution can be achieved. This study extends the analysis framework to SGS modelling in which an optimal model constant (or a set of constants) is estimated based on an argument that spectral energy balance is required between the resolved and unresolved scales. Moreover, it is suggested that the modelling framework can potentially extract an optimal functional form of the SGS stress satisfying the spectral energy balance. The model constant estimation is performed a priori using the corresponding DNS database, but a posteriori results validate the findings.

Following the wavelet formulations developed to analyse the spectral energy transfer of turbulence (Meneveau Reference Meneveau1991; Kim et al. Reference Kim, Bassenne, Towery, Hamlington, Poludnenko and Urzay2018; Nabavi et al. Reference Nabavi, Di Renzo and Kim2022), this study applies the WRMA framework to evaluate the spectrally and spatially local energy transfer due to the unresolved triadic interactions of momentum. Given a SGS closure, this study also uses WMRA to compute the modelled SGS dissipation per scale and location. The current modelling procedure enforces weakly a balance between the two spectral energy fluxes so that the SGS closure is optimised. In this section, the WMRA framework for spectral energy transfer is summarised first, and the specific procedure for estimating optimal model coefficients and discovering an optimal model form is described. The notation is consistent with that of Nabavi et al. (Reference Nabavi, Di Renzo and Kim2022).

For a scalar field $C[\boldsymbol {x}_{0},t]$ (such as a componentwise velocity $u_i$, where $i = 1, 2, 3$) defined at the cell centre locations $\boldsymbol {x}_{0}$ of a computational grid with uniform spacing $\varDelta$ discretising a triply periodic domain, the three-dimensional (3-D) discrete wavelet transformation is applied to obtain a wavelet series expansion of $C[\boldsymbol {x}_{0},t]$. In this study, $\varDelta$ is used to denote the DNS grid spacing. The 24-point Coifman basis is used to conduct wavelet transformation. The standard orthonormal wavelet series expansion (Mallat Reference Mallat1999; Meneveau Reference Meneveau1991) is employed, and the complete details are found in Nabavi et al. (Reference Nabavi, Di Renzo and Kim2022). For an integer scale index $s$, the corresponding wavenumber is defined as $\kappa _s = 2{\rm \pi} /\ell _{s}$ where $\ell _s = 2^{s}\varDelta$ is the dyadic wavelet grid spacing along a direction. The maximum scale level supported by a discrete data of size $N$ along a direction is $\mathcal {S} = \log _2 N$. Thus, $s = 1$ and $\mathcal {S}$ indicate the smallest resolved and largest scales, respectively, with the corresponding length scales being $\ell _1 = 2\varDelta$ and $\ell _\mathcal {S} = 2^\mathcal {S}\varDelta$ (that is, the side length of the domain). Such definitions show that the spectral and spatial resolutions of the wavelet analysis become coarser and finer, respectively, toward small scales (and vice versa for large scales). Using $\ell _s$, the wavelet colocation grid points $\boldsymbol {x}_{s}$ are defined. The wavelet transformation retains information on direction, characterised by the integer directionality index $d = 1, \ldots, 2^3-1$.

For the incompressible Navier–Stokes equations

(2.1)\begin{gather} \frac{\partial u_{i}}{\partial x_{i}} = 0, \end{gather}

(2.2)\begin{gather}\frac{\partial u_{i}}{\partial t}+ u_{j}\frac{\partial u_{i}}{\partial x_{j}} ={-}\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}+\nu\frac{\partial^2 u_{i}}{\partial x_{j}\partial x_{j}} + f_i, \end{gather}

where $u_{i}$ is a velocity component ($i = 1, 2, 3$), $p$ is the hydrodynamic pressure, $\rho$ is the mass density of fluid and $f_i$ is an external forcing, TKE is defined as $k = \tfrac {1}{2}u_i u_i$, and the wavelet energy spectrum at a given scale $s$ is given locally at $\boldsymbol {x}_{s}$ by

(2.3)\begin{equation} {\check{E}\vphantom{}^{(s)}}[\boldsymbol{x}_{s},t]= \frac{2^{{-}3s}}{\kappa_s \ln{2}}\sum_{d=1}^{7} \frac{1}{2}{\check{u}\vphantom{}^{(s,d)}_{i}}[\boldsymbol{x}_{s},t]{\check{u}\vphantom{}^{(s,d)}_{i}}[\boldsymbol{x}_{s},t], \end{equation}

in a triply periodic domain. In a non-periodic domain, the direction dependence is kept in (2.3) as ${\check{E}\vphantom {}^{(s,d)}}$. By formulation, wavelet-transformed quantities retain both spectral and spatial information (that is, dependence on $s$ and ${\boldsymbol {x}}_s$) at a cost of reducing spectral resolution at high wavenumbers. However, the reduced spectral resolution can be compensated by using a more complex basis such as the Coifman family. The directional information (that is, dependence on $d$) is useful to examine anisotropy in turbulence. Following Nabavi et al. (Reference Nabavi, Di Renzo and Kim2022) and Meneveau (Reference Meneveau1991), the transport equation of the wavelet energy spectrum is given as

(2.4)\begin{equation} \frac{\partial}{\partial t}{\check{E}\vphantom{}^{(s)}}[\boldsymbol{x}_{s},t] = {\check{\mathcal{T}}\vphantom{}^{(s)}_{V}}[\boldsymbol{x}_{s},t] + {\check{\mathcal{T}}\vphantom{}^{(s)}_{C}}[\boldsymbol{x}_{s},t] + {\check{\mathcal{T}}\vphantom{}^{(s)}_{P}}[\boldsymbol{x}_{s},t] + {\check{\mathcal{T}}\vphantom{}^{(s)}_{F}}[\boldsymbol{x}_{s},t], \end{equation}

where the right-hand side terms correspond to spectral kinetic energy fluxes due to the viscous diffusion, convective interactions, pressure transport and external forcing, respectively. The complete details of the spectral energy fluxes are referred to Nabavi et al. (Reference Nabavi, Di Renzo and Kim2022). In this study, the mean wavelet energy spectrum ${{\check{E}\vphantom {}^{(s)}}}[t]= \langle {{\check{E}\vphantom {}^{(s)}}}[\boldsymbol {x}_{s},t] \rangle _{\boldsymbol {x}_{s}}$ is examined where $\langle \ \rangle _{\boldsymbol {x}_{s}} = (2^{3s}/N^3 ) \sum _{\boldsymbol {x}_{s}}$ denotes average over the wavelet-colocation grid points at scale $s$. The mean energy fluxes are defined in the same way as ${{\check {\mathcal {T}}\vphantom {}^{(s)}_{X}}}[t]= \langle {{\check {\mathcal {T}}\vphantom {}^{(s)}_{X}}}[\boldsymbol {x}_{s},t] \rangle _{\boldsymbol {x}_{s}}$ where $X = V$, $C$, $P$ and $F$. Summing the mean wavelet energy spectrum and the mean spectral energy fluxes recovers TKE (as required by the Plancherel theorem) and the global energy balance in the physical space, respectively.

For SGS modelling study, the spectral energy fluxes are further formulated in the LES context. Similar to the formulation in the Fourier space (Pope Reference Pope2001), the transport equation of the resolved wavelet energy spectrum is written at $s\geqslant s_{cut}$ as

(2.5)

where the right-hand side terms correspond to the spectral energy fluxes due to different mechanisms. The directional dependence is retained for optimisation. In evaluating the triadic interaction terms in (

), the scale decomposition proposed by Meneveau (Reference Meneveau1991) is used, by which velocity field is decomposed into components smaller and larger than a prescribed filter cutoff scale $s_{cut}$ via

(2.6)\begin{equation} u_{i}[\boldsymbol{x}_{0},t] = u_{i}^{{<}s_{cut}}[\boldsymbol{x}_{0},t] + u_{i}^{\geqslant s_{cut}}[\boldsymbol{x}_{0},t] = \sum_{s=1}^{s_{cut}-1} u_{i}^{(s)}[\boldsymbol{x}_{0},t] + \sum_{s=s_{cut}}^{\mathcal{S}} u_{i}^{(s)}[\boldsymbol{x}_{0},t]. \end{equation}

Using (2.6), the nonlinear convection term is rewritten as

(2.7)\begin{equation} \frac{\partial}{\partial x_{j}} (u_{i}u_{j}) = \frac{\partial}{\partial x_{j}} \left(u_{i}u_{j} - u_{i}^{\geqslant s_{cut}} u_{j}^{\geqslant s_{cut}} \right) + \frac{\partial}{\partial x_{j}} \left(u_{i}^{\geqslant s_{cut}} u_{j}^{\geqslant s_{cut}} \right). \end{equation}

In the current a priori modelling framework, the true LES field is not available and should be estimated. One possible approach is to use a threshold filtering (Goldstein & Vasilyev Reference Goldstein and Vasilyev2004; De Stefano & Vasilyev Reference De Stefano and Vasilyev2012; De Stefano, Dymkoski & Vasilyev Reference De Stefano, Dymkoski and Vasilyev2022), which is not followed by the current study. In the Fourier context, the resolved-scale motions can be estimated using scale-sharp filtering. In general, such practice does not lead to an accurate representation of the true LES solution. Furthermore, numerical errors associated with performing LES a posteriori are not taken into account. This is presumably the case for wavelet scale-cutoff filter (Meneveau Reference Meneveau1991). However, in order to develop the optimisation framework, this study assumes that the superfilter scale velocity field approximates the resolved-scale motions, namely

(2.8)\begin{equation} \bar{u}_i \approx u_i^{\geqslant s_{cut}}, \end{equation}

which solves the filtered momentum equations

(2.9)\begin{equation} \frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j\frac{\partial \bar{u}_i}{\partial x_j} ={-}\frac{1}{\rho}\frac{\partial \bar{p}}{\partial x_i} + \nu\frac{\partial^2 \bar{u}_i}{\partial x_j \partial x_j} - \frac{\partial \tau^{r}_{ij}}{\partial x_j} + \bar{f}_i, \end{equation}

where the overbar denotes a LES filtering operation and $\tau ^{r}_{ij}$ is the deviatoric part of the SGS stress tensor, namely, $\tau ^{R}_{ij} = \tau ^{r}_{ij} + \tfrac {1}{3}\tau ^{R}_{kk}\delta _{ij}$. The isotropic part of the SGS stress is included in the filtered pressure (Pope Reference Pope2001). Similar to the Fourier formulation, the last term in (2.7) represents the resolved triadic interactions, and the corresponding spectral energy flux is estimated a priori as

(2.10)\begin{equation} {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C,f}}[\boldsymbol{x}_{s},t] ={-}\dfrac{2^{{-}3s}}{\kappa_s \ln 2} {\check{u}\vphantom{}^{(s,d)}_{i}}[\boldsymbol{x}_{s},t] \check{u_{j}^{\geqslant s_{cut}}\dfrac{\partial u_{i}^{\geqslant s_{cut}}}{\partial x_{j}}}^{(s,d)}{}[\boldsymbol{x}_{s},t]. \end{equation}

The spectral energy flux due to the resolved triadic interactions can be also written for $s > s_{cut}$ as

(2.11)\begin{equation} {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C,f}}[\boldsymbol{x}_{s},t] = \sum_{s^\prime,s^{\prime\prime} = s_{cut}}^{\mathcal{S}} \sum_{{\boldsymbol{x}}_{s^\prime},{\boldsymbol{x}}_{s^{\prime\prime}}} \sum_{d^\prime,d^{\prime\prime} = 1}^{7} {\check{u}\vphantom{}^{(s,d)}_{i}}[\boldsymbol{x}_{s},t] \ \check{u}\vphantom{}^{(s^\prime,d^\prime)}_{j}[\boldsymbol{x}_{s},t] \check{\dfrac{\partial u_i}{\partial x_j}}^{(s^{\prime\prime},d^{\prime\prime})}[\boldsymbol{x}_{s},t], \end{equation}

where a scale triad involves three resolved scales $(s, s^\prime, s^{\prime \prime })$, all larger than $s_{cut}$. The unresolved component of the triadic energy transfer ${\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}$ is not available a posteriori. It can be estimated a posteriori using, for instance, presumed probability density function (Pope Reference Pope2001), spectral enrichment (Bassenne et al. Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019; Ghate & Lele Reference Ghate and Lele2020) and autonomic closure (King et al. Reference King, Hamlington and Dahm2016). However, those approaches require extra modelling efforts or assumptions. For modelling, it can be evaluated a priori using the spectral energy flux due to the full triadic interactions (represented by the DNS data) via

(2.12)\begin{equation} {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C,r}}[\boldsymbol{x}_{s},t] = {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C}}[\boldsymbol{x}_{s},t] - {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C,f}}[\boldsymbol{x}_{s},t]. \end{equation}

In (2.5), the physical viscous diffusion is primarily an energy sink, that is, ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{V}}} \lesssim 0$. Similar to the LES formulation in the Fourier context, the triadic energy transfer involving the resolved-scale motions, ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,f}}}$, is represented on the LES grid and does not require any modelling, as indicated by the resolved scale triad in (2.11). It is the unresolved triadic interaction term ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}}$ that requires a SGS closure. The pressure transport ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{P}}}$ does not involve any inter-scale transfer. Summing (2.5) over all resolved scales ($s > s_{cut}$), wavelet colocation points, and directions yields a spectral energy balance for turbulence in equilibrium in the LES context:

(2.13)

This study seeks to find a SGS closure that enforces a spectral energy balance constraint between the resolved and unresolved scales. Such energy balance is checked a priori or a posteriori (Meneveau & Katz Reference Meneveau and Katz2000). However, this study enforces the energy balance explicitly so that an analytical form of the optimal (in a sense of spectral energy balance) SGS closure is discovered a priori. Although the modelling procedure uses a corresponding DNS database, the optimisation uncovers an analytical form of the SGS closure terms, by which further generalisation can be made towards different flow conditions. Considering that the convective energy flux is responsible for redistributing energy over a range of scales (Meneveau & Katz Reference Meneveau and Katz2000), it is proposed that the spectral energy transfer due to the modelled SGS stress should balance with ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}}$. In order to do so, (2.9) is formulated similar to (2.5) to obtain a transport equation of the resolved spectral kinetic energy equation as

(2.14)

where the overbar denotes the wavelet energy spectrum and spectral energy fluxes evaluated using the LES velocity and pressure. With an assumption that (

) converges to (2.5) by using an ideal SGS closure, the following condition is enforced as

(2.15)\begin{equation} {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C,r}}[\boldsymbol{x}_{s},t] \approx {\check{\bar{\mathcal{T}}}\vphantom{}^{(s,d)}_{SGS}}[\boldsymbol{x}_{s},t] \quad \mathrm{where} \ s>s_{cut}, \end{equation}

where $s_{cut}= 1, \ldots, \mathcal {S}$, and ${{\check {\bar {\mathcal {T}}}\vphantom {}^{(s,d)}_{SGS}}}$ is the modelled SGS dissipation evaluated a posteriori. One may presume a closure model and conduct LES to obtain ${{\check {\bar {\mathcal {T}}}\vphantom {}^{(s,d)}_{SGS}}}$. If so, however, enforcing (2.15) between the filtered DNS and a posteriori LES solutions locally in space and in time becomes questionable. Instead, the modelled SGS energy flux ${{\check {\bar {\mathcal {T}}}\vphantom {}^{(s,d)}_{SGS}}}$ is approximated by its a priori subfilter-scale (SFS) flux ${{\check {\bar {\mathcal {T}}}\vphantom {}^{(s,d)}_{SFS}}}$ as

(2.16)\begin{equation} {\check{\bar{\mathcal{T}}}\vphantom{}^{(s,d)}_{SGS}}[\boldsymbol{x}_{s},t] \approx {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{SFS}}[\boldsymbol{x}_{s},t] ={-}\dfrac{2^{{-}3s}}{\kappa_s \ln 2} {\check{u}\vphantom{}^{(s,d)}_{i}}[\boldsymbol{x}_{s},t] {\check{\dfrac{\partial \tau^{r}_{ij}}{\partial x_{j}}}\vphantom{}^{(s,d)}_{}}[\boldsymbol{x}_{s},t]. \end{equation}

This formulation is based on an observation that ${\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}$ is the only term in which the unresolved scales are allowed to alter energy contained in the resolved scales (in a priori context). The other mechanisms are either scale local or dependent on the resolved scales only, and thus no closure model is needed. In this study, (2.15) and (2.16) are used as a constraint for which a functional form of a SGS closure $\tau ^{r}$ is optimised in terms of spectral energy transfer. However, due to the stochastic nature of turbulence and a limitation set by a priori approach, the constraint is enforced weakly.

3. Algorithm for discovering an optimal SGS closure

This section describes the algorithm for optimising an existing SGS closure and discovering a new model, both from a perspective of optimal energy balance between the resolved and unresolved scales. It is emphasised that the proposed formulation determines an analytical form of the SGS closure $\tau ^{r}_{ij}$ in (2.9) in such a way that it achieves optimally the spectral energy balance (2.15) with respect to the corresponding DNS data. Although the modelling procedure is a priori, it is not limited to approximating the resolved-scale velocity field. The formulation educes the SGS stress by enforcing weakly the spectral energy balance existing in the analytical LES transport equation (2.5). Thus, it is believed that uncertainties associated with a priori estimation could be reduced to a certain extent.

Figure 1 illustrates schematically the proposed modelling procedure for turbulent flows in a triply periodic box. The overall procedure remains the same for other configurations.

1. (i) Wavelet transformation. Orthonormal wavelet transformation is applied to an instantaneous snapshot obtained from DNS, namely,

   (3.1)\begin{equation} {\check{u}\vphantom{}^{(s,d)}_{i}}[\boldsymbol{x}_{s},t] = \sum_{\boldsymbol{x}_{0}} u_i[\boldsymbol{x}_{s},t] \mathcal{G}^{(s,d)}[\boldsymbol{x}_{0}-\boldsymbol{x}_{s}], \end{equation}
   where $\mathcal {G}^{(s,d)}[\boldsymbol {x}_{0} - \boldsymbol {x}_{s}]$ is a 3-D discrete wavelet basis obtained by a tensor product of one-dimensional wavelet basis.

2. (ii) Estimate the resolved-scale motions. The filter cutoff scale $s_{cut}$ is chosen such that $s_{cut} = \log _2(\bar {\varDelta }_{f}/\varDelta )$ where $\bar {\varDelta }_{f}$ is a priori LES filter width. The resolved-scale motions are estimated using (2.8). Even if the resolved fluid motions are approximated by the filtered velocity field, $\bar {\varDelta } \approx \bar {\varDelta }_{f}$ is not generally true. Instead, most of the literature assume that the LES filter width is several times larger than the LES grid spacing. Wavelet coefficients at $s < s_{cut}$ are set to be zero, which is equivalent to suppressing fluid motions smaller than the a priori LES filter width $\bar {\varDelta }_{f} = 2^{s_{cut}}\varDelta$. Although such operation does not behave as a scale-sharp filter in the same way as the Fourier transformation does, scales smaller than $s = s_{cut}$ are suppressed and only larger (more precisely, superfilter) scales are retained:

   (3.2)\begin{equation} \bar{u}_i[\boldsymbol{x}_{0},t] \approx u_i^{\geqslant s_{cut}}[\boldsymbol{x}_{0},t] = \sum_{s = s_{cut}}^{\mathcal{S}}\sum_{\boldsymbol{x}_{s}}\sum_{d = 1}^7 \check{u}_i[\boldsymbol{x}_{s},t] \mathcal{G}^{(s,d)}[\boldsymbol{x}_{0}-\boldsymbol{x}_{s}]. \end{equation}

3. (iii) Evaluate the targeted spectral energy flux across the cutoff scale. Using the WMRA framework (Meneveau Reference Meneveau1991; Nabavi et al. Reference Nabavi, Di Renzo and Kim2022), the inter-scale energy flux ${\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}$ involving the unresolved-scale motions for a given $s_{cut}$ is computed using (2.10) and (2.12). The inter-scale energy flux is defined at the wavelet colocation points, and its minimum spatial resolution is the same as the corresponding LES grid (that is, $\bar {\varDelta }$). Thus, it can serve as a targeted spectral energy flux for a selected filter cutoff scale, which a ‘correct’ SGS closure should transfer to scales smaller than $s_{cut}$ (that is, subgrid scales).

4. (iv) Inverse wavelet transformation. Wavelet coefficients at $s < s_{cut}$ are discarded and only those at $s \geqslant s_{cut}$ are transformed back to the physical space. Thus, the inverse transformation is obtained on a grid having the same resolution as the LES grid, not the original DNS grid. This step completes performing LES a priori.

5. (v) Evaluate the modelled spectral energy flux on a nominal LES grid. Using the estimated resolved-scale velocity field $u_i^{\geqslant s_{cut}} \approx \bar {u}_i$ defined on the physical-space grid of spacing $\bar {\varDelta }$, the modelled SGS energy dissipation is calculated. An existing SGS closure such as the Smagorinsky model is employed so that the SGS dissipation is calculated using (2.16) parameterised by unknown coefficients $\eta _k$, where $k = 1, 2, \ldots, \mathcal {K}$ with $\mathcal {K}$ being the number of the undetermined model constants (for instance, $\mathcal {K} = 1$ for the Smagorinsky model). The unknown constants are determined at the next optimisation step.

   Although this study evaluates the SGS energy dissipation using an existing closure model, the following generalisation is possible without any reformulation. In other words, a functional form of the SGS stress may remain undetermined, and a general functional relation combining velocity and its gradients is assumed as

   (3.3)\begin{equation} \tau^{r}_{ij} \approx f({\boldsymbol{u}},{\boldsymbol{\nabla}}{\boldsymbol{u}},{\boldsymbol{\nabla}}^2{\boldsymbol{u}}, \ldots), \end{equation}
   with the corresponding unknown coefficients $\eta _k$. Indeed, this generates a combinatorially large number of the candidate terms (that is, $\mathcal {K} \gg 1$) since the nonlinear products and the magnitude of the individual terms are also taken into account. This practical difficulty can be addressed partially by performing optimisation a priori and enforcing the energy constraint weakly. In addition, dimensional, physical and tensorial requirements can eliminate a large number of terms in (3.3) in a way similar to what Silvis et al. (Reference Silvis, Remmerswaal and Verstappen2017) proposed. Using standard algorithms in multivariate optimisation, it is possible to determine the set of unknowns for many candidate terms in (3.3). Still, (3.3) is not most convenient to incorporate for the current study where the feasibility of the optimisation framework is assessed. Thus, (3.3) is not used in this study, and to demonstrate the feasibility of the formulation, four different models for $\tau ^{r}_{ij}$ containing one or two undetermined constants are examined in § 6.

6. (vi) Optimise for the model constants. For a general SGS closure for $\tau ^{r}_{ij}$, the unknowns $\eta _k$ where $k = 1, 2, \ldots, \mathcal {K}$ are determined a priori by enforcing (2.15) weakly in the least-squares sense. Although the number of the unknowns (that is, $\mathcal {K}$) can be very large depending on the choice of the candidate terms in $\tau ^{r}_{ij}$, the number of data points on the filtered DNS database is usually much larger than $\mathcal {K}$, making (2.15) a strongly overdetermined system. By doing so, it is anticipated that uncertainties associated with the present a priori modelling approach and numerical errors associated with wavelet transformation could be reduced.

   For a given filter cutoff scale $s_{cut}$ and a flow snapshot at $t = t_m$, where $m = 1, 2, \ldots, \mathcal {M}$, the unknown model coefficients are determined by minimising the scale-normalised, squared error between the unresolved triadic and the modelled SGS energy fluxes via

   (3.4)\begin{equation} \eta^{(s_{cut})}_k[t_m] = \arg \min_{\hspace{-0.6cm}\eta_k} \sum_{s = s_{cut}}^\mathcal{S}\sum_{\boldsymbol{x}_{s}}\sum_{d = 1}^7 \frac{2^{{-}3s}}{\kappa_s \ln 2} \lvert {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{C,r}}[\boldsymbol{x}_{s},t_m] - {\check{\mathcal{T}}\vphantom{}^{(s,d)}_{SFS}}[\boldsymbol{x}_{s},t_m] \rvert^2, \end{equation}
   where $k = 1, 2, \ldots, \mathcal {K}$. It should be noted that ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}}$ remains the same for different closure models. The minimisation problem is solved using a standard multivariate optimisation algorithm. In addition, the minimisation is repeated for a series of flow snapshots, and $\eta _k$ is averaged over time if turbulence is stationary. The optimisation procedure is repeated over a range of the filter cutoff scale $s_{cut}$. In addition, the minimisation problem can be penalised by including a regularisation parameter to promote the sparsity of the unknown constants, useful to make the discovered SGS closure compact. However, the present study does not include any regularisation term. The modelling algorithm is summarised as follows.

The discovered SGS closure is examined for its mathematical and physical properties. In addition, the closure is compared with the existing LES models to understand the advantages and drawbacks, useful in improving the algorithm outlined in this section. The discovered model is implemented in a computational fluid dynamics code and tested a posteriori for canonical flows.

Figure 1. Proposed SGS modelling using a wavelet multiresolution framework.

This section concludes by summarising some of the differences and similarities between the current wavelet framework and the optimal-LES formulation of Langford & Moser (Reference Langford and Moser1999). A primary conceptual difference is that Langford & Moser (Reference Langford and Moser1999) adopted a concept of an ‘ideal’ LES in a context of one-time, multipoint statistical accuracy, whereas the current formulation focuses on numerical stability as a consequence of matched spectral energy transfer across a nominal filter cutoff scale. In formulating optimisation, Langford & Moser (Reference Langford and Moser1999) minimised a least-squares difference between the ‘true’ SGS force and modelled SGS force, both found on the right-hand side of the filtered Navier–Stokes equations. The wavelet framework uses the unresolved triadic energy flux and modelled SGS dissipation, both found on the right-hand side of the filtered kinetic energy equation. In performing optimisation, Langford & Moser (Reference Langford and Moser1999) used stochastic estimation (Adrian Reference Adrian1990) with $N$-point correlations as model inputs so that the unknown estimation kernels are determined. The current modelling is formulated in a wavelet multiresolution context, and available DNS fields are provided so that model constants are determined. Both formulations allow to keep an arbitrary number of terms in the modelled SGS stress, although the number is not known a priori. In Langford & Moser (Reference Langford and Moser1999), a correct inter-scale energy transfer towards the SGS motions was reported as a consequence of using an optimal LES formulation, whereas such energy transfer was evaluated a priori and used by the wavelet framework as a part of optimisation criteria. Using wavelet offers additional benefits by further constraining optimisation. For instance, the unresolved triadic energy flux can be conditioned to be negative so that energy backscatter effects on SGS closure can be studied. Directional information available in wavelet statistics allows to study anisotropy effects on SGS modelling. Overall, although both formulations are rooted in the concept of optimisation for determining SGS models, technical differences seem to be significant.

4. Simulations of homogeneous isotropic turbulence

The proposed modelling strategy is tested for forced homogeneous isotropic turbulence (HIT) simulated by DNS and a series of LES, respectively. The Smagorinsky closure (Smagorinsky Reference Smagorinsky1963), in particular, its dynamic variant (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) is employed for validation. In the classic Smagorinsky model (CSM), the Boussinesq approximation is invoked, and the deviatoric part of the SGS stress tensor $\tau ^{r}_{ij}$ is modelled as

(4.1)\begin{equation} \tau^{r}_{ij} \approx{-} 2 C^2 \bar{\varDelta}^2 \lvert\bar{S}\rvert \bar{S}_{ij}, \end{equation}

where $C$ is the Smagorinsky constant, $\bar {S}_{ij}$ is the filtered strain-rate tensor, and $\lvert \bar {S}\rvert$ is the magnitude of the filtered strain-rate tensor, namely,

(4.2)\begin{equation} \lvert\bar{S}\rvert = \sqrt{2 \bar{S}_{ij} \bar{S}_{ij}}. \end{equation}

Thus, the eddy viscosity is defined as ${\nu _{t}} = C^2 \bar {\varDelta }^2 \lvert \bar {S}\rvert$ so that $\tau ^{r}_{ij} = - 2{\nu _{t}} \bar {S}_{ij}$. In the DSM (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991), the Smagorinsky constant $C$ is not prescribed but determined using the dynamic procedure. Computational work is done using the NGA solver (Desjardins et al. Reference Desjardins, Blanquart, Balarac and Pitsch2008) that supports scalable, massively parallel simulations of incompressible turbulent flows. A kinetic-energy-conserving, second-order, centred finite-difference discretisation is employed, which is fully validated in the past for the LES of various turbulent flows. Additional tests are performed to confirm that the NGA code is conservative and its numerical errors do not affect the conclusion of this study.

Both DNS and LES runs are performed in the same computational domain of a triply periodic box with the side length of $L = 2{\rm \pi}$. A total of $512^3$ control volumes discretise the computational box of DNS (thus, $\varDelta = 2{\rm \pi} /512$), and the LES grid resolution varies from $\bar {\varDelta } = 2{\rm \pi} /256$ ($256^3$ grid points) to $\bar {\varDelta } = 2{\rm \pi} /32$ ($32^3$ grid points) by a factor of 2. Thus, 4 different LES grid resolutions ($32^3$, $64^3$, $128^3$ and $256^3$ grid points) are examined in this study. Variables are made dimensionless in the same way as Bassenne, Moin & Urzay (Reference Bassenne, Moin and Urzay2018). To sustain turbulence, an external forcing (Bassenne et al. Reference Bassenne, Urzay, Park and Moin2016) is applied with TKE kept constant over time, that is $k \approx k_\infty$, in which variables having a subscript $\infty$ are measured after the DNS prediction becomes statistically stationary. The Reynolds number based on the Taylor microscale $\lambda$ is ${Re}_\lambda = u_{\ell }\lambda /\nu =85$, where $u_{\ell }$ is the integral velocity (also, the root mean square of the fluctuating velocity). For DNS, spatial resolution is $\kappa _{max}{\ell _{k,\infty }}=1.6$, where $\kappa _{max}={\rm \pi} /\varDelta$ is the maximum wavenumber supported by the DNS grid and ${\ell _{k,\infty }}$ is the Kolmogorov length scale. The nominal eddy-turnover time $t_{\ell,\infty } = \tfrac {1}{3}(u_i^\prime u_i^\prime )/\epsilon$ estimated using the DNS result is used to normalise time. The velocity field is initialised to be solenoidal and isotropic using a prescribed Passot–Pouquet spectrum (Passot & Pouquet Reference Passot and Pouquet1987). The computational time-step size is determined so that the resulting Courant number is approximately $0.5$. For analysis and modelling, 12 instantaneous snapshots are sampled every $1.6 t_{\ell,\infty }$ after the initial transient effects become insignificant at $t \geqslant 14.4 t_{\ell,\infty }$ (see figure 2a).

Figure 2. Time evolution of (a) the normalised viscous dissipation rate and (b) the Smagorinsky constant determined by the dynamic process. In (b), the horizontal dotted line corresponds to the theoretical estimation of Lilly (Reference Lilly1967).

Figure 2(a) shows the time evolution of the dissipation rate. As the LES grid is refined, the dissipation rate converges to its true value predicted by DNS. Following the initial overshoot, the dissipation rate remains more or less constant for all cases (even though it is not forced explicitly to do so), qualitatively similar to Bassenne et al. (Reference Bassenne, Urzay, Park and Moin2016). In the LES context, Lilly (Reference Lilly1967) estimated the Smagorinsky constant via $C_{S}^2 = {1}/{{\rm \pi} ^2} ({2}/{3 C_{K}})^{3/2} \approx 0.1732^2$, where $C_{K} = 1.5$ is the Kolmogorov constant and the LES filter width is assumed in the inertial subrange. In figure 2(b), the theoretical value is compared with the dynamically determined Smagorinsky constants. As the LES grid is refined, the model constant decreases monotonically but does not converge to zero (as expected). With $64^3$ grid points, the dynamic model predicts the Smagorinsky constant very close to the theoretical estimation, implying that the LES filter width of the $64^3$-point grid LES using the standard dynamic model lies in the inertial subrange.

For validation, the Fourier energy spectrum is plotted in figure 3. The current DNS result shows good agreement with that of Bassenne et al. (Reference Bassenne, Moin and Urzay2018). For the $32^3$-point grid, the LES result of Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019) is also compared. Under the same flow condition and using the same SGS closure, the energy spectra of Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019) and the present study agree well with each other. On $128^3$-point grid, it appears that the LES grid width is in the dissipation range, and LES achieves the DNS-like resolution with $256^3$ points as the energy spectra are indistinguishable between DNS and LES.

Figure 3. Time-averaged Fourier energy spectra for LES results using (a) the DSM and (b) no explicit SGS model.

Figure 3(b) shows results when no explicit SGS model is used. On coarse grids (${\lesssim }64^3$ points), LES prediction is largely inaccurate due to insufficient spatial resolution and missing SGS dissipation. The energy spectrum becomes similar to that of DNS on $256^3$-point grid. However, mild energy pileup is observed at high wavenumbers on $128^3$-point grid. Figure 3(b) shows that on ${\lesssim }64^3$ grids, SGS models are needed for accurate turbulence prediction. It also shows that numerical errors associated with the computational code and its discretisation are not large enough to affect the conclusion of the present modelling and simulation.

5. Mean wavelet statistics for the LES database

Wavelet analysis is performed on the LES-generated snapshots of the forced HIT to examine the LES-grid convergence of the mean wavelet energy fluxes. Thus, the analysis results reported in this section are a posteriori. The standard DSM is used. Figure 4(a) shows the wavelet energy spectra for the current DNS and LES predictions. Due to the good spectral accuracy of the Coifman basis coif4 used in this analysis, the wavelet spectra follow closely the Fourier energy spectrum. LES results using ${\geqslant }64^3$ points show insignificant differences with each other. Energy pileup at high wavenumbers is expected and attributed to the spectral leakage of wavelet transformation (Mallat Reference Mallat1999). In figure 4(b), the wavelet enstrophy spectra are shown. Unlike the wavelet energy spectrum, the LES-grid dependence is more pronounced, showing a relevance of using the enstrophy spectrum to show the LES-grid effects. In both figures, LES achieves a DNS-like resolution with ${\gtrsim }256^3$ grid points.

Figure 4. Time-averaged wavelet spectra of (a) TKE and (b) enstrophy. In (a), the Fourier energy spectra of Bassenne et al. (Reference Bassenne, Moin and Urzay2018) (DNS on $512^3$ grid points) and Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019) (DSM on $32^3$ grid points) are shown. In (b), the Fourier enstrophy spectrum of Bassenne et al. (Reference Bassenne, Moin and Urzay2018) (DNS on $512^3$ grid points) is shown.

Figures 5(a–e) show the spectral energy fluxes due to different mechanisms of the momentum transport in (2.5). For a stationary HIT, the weighted sum of all the energy fluxes is equal to zero (less than $0.01\,\%$ of $u_{k,\infty }^3$ for all the cases) so that the instantaneous kinetic energy remains the same, namely, ${\rm d}\bar {k}/{\rm d}t \approx 0$, as required by (2.13). For its nature of redistributing energy across scales (Meneveau & Katz Reference Meneveau and Katz2000), the wavenumber-weighted sum of $\langle {{\check {\mathcal {T}}\vphantom {}^{(s)}_{C,f}}}\rangle _{\boldsymbol {x}_{s}}$ is negligible. The physical viscous and modelled SGS energy transfer mechanisms are purely dissipative on average across all scales (see figure 5b,c) and thus negative if integrated with respect to $\kappa _s$. For the current incompressible flow, the pressure transport mechanism is much less important than the others except for the largest scale (see figure 5d). The linear forcing that sustains turbulence is a net energy source by design (Bassenne et al. Reference Bassenne, Urzay, Park and Moin2016), and thus its mean energy input is positive (see figure 5e).

Figure 5. Time-averaged mean wavelet spectral energy fluxes due to (a) the resolved convective momentum transport, (b) physical viscosity, (c) modelled SGS energy transfer, (d) pressure gradient and (e) external forcing. In (a,d), the horizontal dotted line denotes the neutral energy transfer.

Figure 5(a) shows the mean wavelet energy flux due to the resolved triadic interactions. For DSM ($128^3$) and DSM ($256^3$), $\langle {{\check {\mathcal {T}}\vphantom {}^{(s)}_{C,f}}} \rangle _{\boldsymbol {x}_{s}}$ nearly coincides with the DNS flux, whereas the lower-resolution LES results show clear deviation. Thus, $N^3 \gtrsim 128^3$ points is sufficient to directly resolve the mean triadic interactions of this forced HIT. As figure 5(b) shows, the mean dissipation range is nearly resolved using $N^3 = 256^3$ grid points, whereas it is not completely the case for DSM ($128^3$). As a result, the SGS energy transfer for DSM ($128^3$) is small but not negligible as shown in figure 5(c). On the other hand, the mean spectral energy fluxes for DSM ($256^3$) are nearly indistinguishable with the DNS fluxes (except for the pressure mechanism, as discussed in the following), and its SGS energy transport remains virtually neutral.

The lower-resolution LES ($32^3$ and $64^3$ grid points) demonstrates a much stronger dependence on the SGS energy transport. For the triadic energy transfer (figure 5a), the crossover wavenumber at which scales switch between a net energy source and sink moves toward large scales as LES grid becomes coarse. This indicates an artificial reduction of the inertial subrange and extension of the dissipation range (see also figure 5b). The modelled SGS dissipation has $O(1)$ impacts on spectral kinetic energy transfer, as demonstrated in figure 5(c).

In figure 6, the mean spectral energy fluxes due to the resolved triadic interactions, physical viscous transport and modelled SGS energy transfer are plotted to show the relative importance of the SGS energy transfer over the physical viscous dissipation mechanism at several LES-grid resolutions. As the LES-grid resolution improves (figure 6a–d), contributions from the mean SGS energy transport decay to near-zero for the dynamic Smagorinsky closure, leaving the physical viscous dissipation as the only energy sink. Even if the model constant $C$ (shown in figure 2b) does not converge to zero in the DNS-resolution limit (for instance, $256^3$-point grid), the mean SGS energy flux converges to zero in the wavelet context. For DSM ($64^3$), the viscous and SGS energy fluxes are nearly identical with less than $2.5\,\%$ errors. For DSM ($32^3$), the modelled SGS dissipation is a dominant energy sink across all scales. Similarity between the viscous and SGS energy fluxes is expected as the Smagorinsky closure is built upon the eddy-viscosity hypothesis. As a consequence, the modelled SGS stress takes the same functional form as the physical viscous stress and increases the effective viscosity of turbulence (or decreases the effective Reynolds number), explaining the extended dissipation range observed in figure 5(b). Figure 6 shows the LES-grid convergence of the modelled SGS dissipation but also emphasises in the wavelet context the well-known limitation of the Bousinessq hypothesis.

Figure 6. Time-averaged mean wavelet spectral energy fluxes due to the resolved triadic interactions, physical viscous transport and modelled SGS energy transfer for (a) DSM ($32^3$), (a) DSM ($64^3$), (a) DSM ($128^3$) and (a) DSM ($256^3$). The horizontal dotted line denotes the neutral energy transfer.

Figure 7 shows the mean cumulative energy flux of convective transfer involving the unresolved (or subfilter-)scale motions. This energy flux describes the amount of energy flow rate across a grid (or filter) cutoff scale $s_{cut}$ involving the unresolved (or subfilter) scale motions (that is, no interactions between resolved scales are included). Thus, it is evaluated as a wavenumber-weighted sum of the detailed convective spectral energy flux (2.12) as

(5.1)\begin{equation} \check{\varPi}\vphantom{}^{(s_{cut},d)}_{{SFS}}[{\boldsymbol{x}}_{s_{cut}}] = \sum_{m = s_{cut}}^\mathcal{S} (\kappa_m \ln 2) \check{\mathcal{T}}\vphantom{}^{(m,d)}_{{C,r}}[{\boldsymbol{x}}_m\vert{\boldsymbol{x}}_{s_{cut}}], \end{equation}

where the superfilter-scale detailed fluxes are evaluated on the cutoff scale grid so that there is no redundancy in information and the Plancherel formula is satisfied between the physical and wavelet spaces. Note that $s_{cut} >1$ since $s_{cut} = 1$ does not give any SFS information. A positive value of the cumulative energy flux denotes net energy transfer toward scales smaller than $s_{cut}$ (that is, forwardscatter or downscale energy transfer), whereas its negative mean value indicates net energy backscatter towards superfilter scales (Meneveau Reference Meneveau1991). For all simulation results, net forward energy scatter is observed, as expected for the current HIT setup. For DNS, a peak of the mean cumulative energy flux is found at the dimensionless wavenumber of $0.2$ within the inertial subrange (see figure 4a), wherein the mean cumulative energy flux nearly coincides with the mean energy dissipation rate. Figure 7 shows that the physical inertial subrange is predicted using ${\gtrsim }128^3$ LES-grid points. This can be compared with figure 2(b) where $64^3$-point LES predicts the theoretical Smagorinsky constant most accurately. Such comparison implies that estimating a correct model constant in LES using, for instance, the dynamic procedure is not always sufficient to predict a physically correct turbulent flow and also shows an importance of model error.

Figure 7. Time-averaged spectra of the cumulative triadic energy flux. Only $s_{cut} = 2, \ldots, \mathcal {S}$ are shown for each case.

It should be noted that all the analyses made in this section concern wavelet statistics averaged in space per scale. Thus, the spatially local characteristics of wavelet statistics is averaged out. Such spatial locality is important for high-Reynolds-number turbulent flows, but even for homogeneous turbulence, spatial locality is known to be dynamically significant at small scales. For instance, figure 6 shows a posteriori balance involving several different mechanisms of spectral energy transport. Locally in space, however, the stochastic nature of turbulence causes a dynamic energy imbalance, some of which should be taken care of by the SGS closure model. However, existing SGS models are not designed to consider such spatial locality, partially due to the fundamental modelling assumption. This study is motivated by an argument that the local wavelet statistics (that is, before taking the averaging operation $\langle \phantom {\;} \rangle _{\boldsymbol {x}_{s}}$) can be made useful to overcome the limitation set by inherent modelling assumptions and guide the SGS modelling study by discovering an optimal closure model in terms of spectral energy balance such as (2.15).

6. SGS model optimisation

6.1. Discovering spectrally optimal model constants

An analytical form of a SGS closure that optimally describes the spectral energy balance (2.15) is discovered using the DNS database. Since the true SGS stress is not available a posteriori, this study performs optimisation a priori using the triadic interactions involving the SFS motions (2.12) and the modelled SGS stress (2.16), both evaluated a priori. This approach seems to allow only a priori interpretation for the discovered SGS terms. However, since the spectral energy flux requirement is enforced weakly, the discovered model is expected to be affected indirectly by the discrepancies between a priori and a posteriori formulations. Instead of assuming a most general form of SGS closure such as (3.3), it is possible to compute the SGS stress using existing closures such as the Smagorinsky model (4.1). If so, however, any error associated with the model form persists, fundamentally limiting the modelling efforts and its extension to a range of turbulent flows. The algorithm outlined in § 3 is not limited to a particular closure model. However, for demonstration purposes, this study focuses on existing closure models and the optimisation results are presented for the forced HIT described in § 4 and analysed in § 5. Specific SGS closure models examined in this study are given by

(6.1)\begin{gather} \tau^{r}_{ij} \approx{-}2 \big(\eta_1 \bar{\varDelta}^2 \lvert\bar{S}\rvert\big) \bar{S}_{ij} \quad \text{(one-parameter Smagorinsky-type closure)}, \end{gather}

(6.2)\begin{gather}\tau^{r}_{ij} \approx{-}2 \left(2.5 \eta_1 \sqrt{\frac{B_\beta}{\alpha_{ij}\alpha_{ij}}}\right) \bar{S}_{ij} \quad \text{(one-parameter Vreman-type closure)}, \end{gather}

(6.3)\begin{gather}\tau^{r}_{ij} \approx \eta_1 \bar{\varDelta}^2 \frac{\partial \bar{u}_i}{\partial x_k} \frac{\partial \bar{u}_j}{ \partial x_k} \quad \text{(one-parameter nonlinear gradient closure)}, \end{gather}

(6.4)\begin{gather}\tau^{r}_{ij} \approx \eta_1 \bar{\varDelta}^2 \frac{\partial \bar{u}_i}{\partial x_k} \frac{\partial \bar{u}_j}{ \partial x_k} -2 \big(\eta_2 \bar{\varDelta}^2 \lvert\bar{S}\rvert\big) \bar{S}_{ij} \quad \text{(two-parameter Clark-type closure)}. \end{gather}

In (6.2), $\alpha _{ij} = \partial \bar {u}_i/\partial x_j$ is the resolved velocity gradient tensor, and $B_\beta$ is the second invariant of $\beta _{ij} = \bar {\varDelta }^2 \alpha _{mi} \alpha _{mj}$ (Vreman Reference Vreman2004). In the two eddy-viscosity models (6.1) and (6.2), the terms inside the parenthesis correspond to the eddy viscosity ${\nu _{t}}$ of the models, respectively. In (6.1), the unknown constant $\eta _1$ is expected to converge to the Lilly's theoretical prediction $C_{S}^2 = 0.1732^2$ if the nominal LES filter width lies in the inertial subrange where the CSM (4.1) is designed. Since the theoretical value of the Vreman model constant for HIT is $c = 2.5C_{S}^2 = 2.5(0.17)^2 = 0.07$ (Vreman Reference Vreman2004), the unknown value $\eta _1$ in (6.2) is also expected to converge to $0.17^2$. For the nonlinear gradient model (6.3) (Clark, Ferziger & Reynolds Reference Clark, Ferziger and Reynolds1979), the unknown constant is expected to be close to $1/12$ at a DNS-like grid resolution (Pope Reference Pope2001). However, such condition is unlikely to be obtained in the current a posteriori grid resolution (see § 4). It should be noted that no explicit constraint is made on the sign of the unknown $\eta _k$, where $k = 1, \ldots, \mathcal {K}$. Employing a constrained optimisation algorithm can find optimal constants requiring to have a certain sign.

In the presence of strong downscale energy transfer (that is, forward energy scatter) across a cutoff scale, ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{C,r}}}[\boldsymbol {x}_{s},t]$ becomes strongly positive locally in space. The spectral energy balance (2.14) dictates that ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{SGS}}}[\boldsymbol {x}_{s},t]$ should be sufficiently large. Thus, the current a priori optimisation is expected to predict a large model constant to dissipate large spectral energy transfer toward the subfilter scales. A large model constant is also expected where the nominal LES grid resolution is coarse (equivalently, low cutoff wavenumber). As the LES grid width $\bar {\varDelta }$ becomes large, more energy should be dissipated by the SGS closure is supposed, thus making the optimised model constant larger. On the other hand, the optimised value of $\eta _k$ is expected to be very small if a priori LES grid resolution is sufficiently high and, thus, the unresolved triadic interaction is weak in transferring energy.

6.2. Eddy-viscosity SGS closures

Optimisation is performed to determine an unknown constant $\eta _1$ in (6.1) and (6.2). Figure 8 compares the optimal model constants for the Smagorinsky-like and Vreman-like models for a range of the filter cutoff scales. In the inertial subrange (see also figure 4a), the optimal constants for both models recover the Lilly's theoretical prediction $C_{S}^2 = 0.1732^2$, which validates that the current optimisation produces a result consistent to the theoretical estimation (Lilly Reference Lilly1967). Since the Smagorinsky closure is designed specifically for the inertial subrange, the optimal model constant is not expected to remain the same as the theoretical value $C_{S}^2$ across all cutoff scales. Still, the current optimisation framework is useful to estimate (at least, a priori) if a given closure provides a sufficient amount of dissipation to balance the inter-scale energy transfer towards the subfilter scales.

Figure 8. The optimised model constants for the one-parameter Smagorinsky-like (6.1) and Vreman-like (6.2) closure models. A posteriori estimation of the Smagorinsky constants is shown as $\ast$ (green) as a function of the grid-based wavenumber scaled by $2^{-1}$.

If the LES filter width is reduced, figure 8 shows that optimisation predicts progressively smaller constants for both Smagorinsky- and Vreman-like closures, presumably toward zero in the limit of $\kappa _{cut} \rightarrow \infty$. This is a desired property of SGS closure model if LES grid is refined and, thus, more turbulent motions are resolved. Figure 8 suggests that such property is closely tied to the optimal spectral energy balance the current optimisation enforces. It is also suggested that a SGS model designed by the current optimisation formulation has a potential to ensure such property (at least a priori and to be checked a posteriori in a future work). Indeed, the physical viscous mechanism dominates on average over the SGS dissipation at such grid resolution regardless of the value of model constant (Pope Reference Pope2001), which can be also seen a posteriori in figure 6(c,d). Still, spatially local SGS energy transfer is not taken into consideration in such argument, which could be significant where small-scale turbulent motions are known to play a crucial role such as those in chemically reacting turbulent flows. Furthermore, the SGS model constant is often used to compute the SGS scalar transport, which may amplify the model error in two-way interactions involving fluid momentum. In addition, figure 8 demonstrates that the standard practice of using the same Smagorinsky constant $C_{S}$ at refined LES-grid resolution cannot avoid excessive dissipation as many previous studies report. In addition, figure 8 provides quantitative information as to which value of the model constant should be used to achieve the spectral energy balance when LES grid resolution is varied. As $\kappa _{cut} \rightarrow 0$, the estimated model constants become almost several times larger than $C_{S}$ for both models. In such a limit, the SGS model is expected to converge to the Reynolds stress. Thus, the observed trend in figure 8 is qualitatively consistent since the SGS model is responsible for nearly all energy dissipation.

In figure 8, the dynamically determined Smagorinsky constants reported in figure 2(b) are plotted together. It should be noted that there is a fundamental ambiguity in directly comparing the filter cutoff wavenumber $\kappa _{cut}$ defined in the context of a priori wavelet analysis (see Schneider & Vasilyev Reference Schneider and Vasilyev2010) and the grid cutoff wavenumber $2{\rm \pi} /\bar {\varDelta }$ in a posteriori sense. In addition, it is not straightforward to estimate a quantitative relation between the wavelet filtering (3.2) used to evaluate ${{\check {\mathcal {T}}\vphantom {}^{(s,d)}_{SGS}}}$ and the grid filtering implicitly applied to the current DSM formulation. In this study, the wavenumber used to plot the model constants from the individual DSM simulations is scaled by $1/2$. Although the factor of $1/2$ is arbitrary, this choice is justified by the observation that in figure 2(b), DSM ($64^3$) predicts the model constant $C^2 \approx C_{S}^2$ and the Smagorinsky-like closure (6.1) has a spectrally optimal constant of $\eta _1 \approx C_{S}^2$ at $\kappa _{cut}{\ell _{k,\infty }} \approx 0.39$. Assuming that the inertial subrange behaviour of the static Smagorinsky model and DSM are the same, the filter and grid cutoff wavenumbers are matched by introducing the artificial factor of $1/2$. It should be also noted that in the literature of the implicitly filtered LES it is usually assumed that the LES filter width is several times larger than the LES grid width. This study suggests that as far as the optimal spectral energy transfer is concerned, the constant is $2$. However, this number is a consequence of using the wavelet filtering for a priori analysis and likely to be specific to the current study only. Figure 8 shows that the DSM becomes suboptimal in the context of the spectral energy balance at $\kappa _{cut} \rightarrow 0$. Thus, the eddy viscosity estimated by the dynamic procedure is expected to be smaller than the optimal value at a coarse LES grid resolution. Indeed, the current LES performed using a finite difference cannot rule out possibilities that numerical errors contribute to the reduced eddy viscosity, which will be examined in a future work. Regardless, the dynamically determined Smagorinsky constant behaves in a way that is spectrally optimal as the LES grid is refined so that the inertial subrange is resolved.

Figure 8 demonstrates that the proposed optimisation formulation generates SGS model constants in a way consistent with the theoretical and a posteriori numerical studies. It also shows qualitatively correct grid convergence of the model coefficients for both Smagorinsky- and Vreman-type eddy viscosity models. This is useful to understand the existing SGS models, to assess and compare their performance over a range of filter cutoff scales, to design new SGS closure models and potentially to predict SGS scalar fluxes or SGS interactions with different physics. Table 1 summarises the model constants optimised for each SGS closure model considered in this study. At the filter cutoff wavenumber $\kappa _{cut}{\ell _{k,\infty }} = 0.2$, both Smagorinsky- and Vreman-type eddy viscosity models select the constants close to $C_{S}^2 = 0.1732^2 = 0.03$, which is also obtained for DSM ($64^3$).

Table 1. The optimised model constants for a range of the filter cutoff wavenumbers. The equivalent LES grid resolution is determined so that the inertial subrange behaviour is matched between a priori and a posteriori results. The theoretical estimation of the Smagorinsky constant in the inertial subrange is $C_{S}^2 = 0.1732^2 = 0.03$, and the gradient model (Clark et al. Reference Clark, Ferziger and Reynolds1979) has its coefficient equal to $1/12 = 0.083$ as $\kappa _{cut} \rightarrow \infty$.

The above discussions are obtained for the time-averaged model constants found by the optimisation algorithm. Although the forced HIT is statistically stationary, the instantaneous variation of the model constants $\eta _k[t_m]$, where $m = 1, 2, \ldots, 12$, is examined. In figure 9, each line corresponds to the result from each snapshot, and the lines with filled symbols correspond to the time-averaged model constants (reported in figure 8). Although the optimised constants show similar trends at different time instants, they also show a significant temporal variation, in particular, for large filter widths. This is attributed not only to the relatively smaller sample size of large-scale wavelet statistics but also partially to the fundamentally inhomogeneous and anisotropic nature of large-scale fluid motions.

Figure 9. The optimised model constants for the (a) Smagorinsky-type and (b) Vreman-type closures for all 12 DNS snapshots. The filled symbols indicate the time-averaged model constants reported in figure 8.

6.3. Gradient SGS model

The nonlinear gradient model (Clark et al. Reference Clark, Ferziger and Reynolds1979) is developed based on an argument that in the DNS-resolution limit $\bar {\varDelta }/\ell _k \ll 1$ the exact SGS stress behaves proportional to the nonlinear product of velocity gradients via $\tau ^{R}_{ij} \sim ({\bar {\varDelta }^2}/{12})({\partial \bar {u}_i}/{\partial x_k}) ({\partial \bar {u}_j}/{\partial x_k}) + O(\bar {\varDelta }^4)$. The model is, however, known to be unstable and is thus combined with the static Smagorinsky closure for numerical stability (Vreman, Geurts & Kuerten Reference Vreman, Geurts and Kuerten1996) or its dynamic variant is used (Vreman et al. Reference Vreman, Geurts and Kuerten1996, Reference Vreman, Geurts and Kuerten1997). Using the pure gradient-type closure given in (6.3) (thus, the original Clark model without the eddy-viscosity part), the SGS energy flux is estimated a priori using (2.16), and the optimisation is performed to determine the unknown constant $\eta _1$ via (3.4). In the DNS-like resolution (that is, $\kappa _{cut} \rightarrow \infty$), $\eta _1$ is expected to converge to $1/12$.

Figure 10(a) shows the model constant $\eta _1$ optimised for the pure gradient-type model (6.3). Unlike the eddy-viscosity models discussed in § 6.2, the discovered model constant does not show monotonic variation over scales. In addition, the optimised constants are consistently larger than the theoretical estimation $1/12$. However, this is expected because the theoretical value is justified only when the Kolmogorov scale is well resolved (Clark et al. Reference Clark, Ferziger and Reynolds1979).

Figure 10. The optimised model constants for the pure gradient-type closure (6.3): (a) time-averaged and (b) instantaneous model constants where the line with filled symbols indicates the time average.

Similar to figure 9, the optimised constants for all 12 snapshots are plotted in figure 10(b). It is observed that at coarse LES-grid resolution, the optimised constant sometimes takes a negative value so that the spectral energy balance is achieved. Since the current HIT has a net energy transfer towards subfilter scales (see figure 7), figure 10(b) shows that the nonlinear gradient product in (6.3) increases (rather than dissipates) the resolved kinetic energy, showing its inadequacy as a SGS closure model if LES grid resolution is coarse. This observation is consistent to an anti-dissipative nature of the original Clark model when $\eta _1 = 1/12$ is used by Clark et al. (Reference Clark, Ferziger and Reynolds1979) and Vreman et al. (Reference Vreman, Geurts and Kuerten1996) where numerical instability is observed unless an additional dissipation mechanism by, for instance, the static Smagorinsky model is included.

The present optimisation procedure determines not only scale-dependent model constants but also points out if a model form of a certain SGS closure produces physically consistent effects on LES. Figure 10(a) suggests that due to its model form based on velocity gradient product, the gradient model constant at very coarse LES grid resolution should be sometimes negative so that the triadic energy transfer towards the unresolved scales can be balanced by the modelled SGS stress. The constant is scale dependent as indicated by table 1, but should converge to $1/12$ in the DNS-like grid resolution (Clark et al. Reference Clark, Ferziger and Reynolds1979). The constant appears to decrease towards the theoretical value $1/12$ as $\bar {\varDelta }$ decreases further. However, with the current DNS grid resolution of $512^3$ grid points where $\varDelta /\ell _k = O(1)$, it is difficult to confirm such convergence.

6.4. Clark model

The proposed SGS modelling framework can be used to optimise a SGS closure having more than a single unknown constant so that the spectral energy criterion (2.15) is satisfied. The optimisation is applied to the Clark-like SGS closure (6.4). The scale-dependent, spectrally optimal model constants are determined by the optimisation algorithm and plotted in figure 11 (also see table 1). Note that they are determined to satisfy (3.4) simultaneously but reported separately in figures 11(a) and 11(b). In figure 11(a), the gradient part of the Clark-like closure is compared with the optimised gradient model (6.3). Similarly, figure 11(b) compares the optimised eddy-viscosity component of (6.4) with the optimised Smagorinsky closure (6.1). In Clark et al. (Reference Clark, Ferziger and Reynolds1979), the eddy viscosity is arbitrarily included to ensure numerical stability, which is however reported to cause excessive dissipation in transitional regimes (Vreman et al. Reference Vreman, Geurts and Kuerten1996). As shown in figure 11(b), the eddy-viscosity contribution remains similar to the standalone Smagorinsky closure at $\kappa _{cut}{\ell _{k,\infty }} < 0.39$. As a result of the additional dissipation, the spectrally optimal condition (2.15) requires the other constant $\eta _1$ associated with the gradient part of the model to be close to the theoretical value $1/12$ (see figure 11a). Thus, if the standard Clark model is used (that is, $\eta _1 = 12$ and $\eta _2 = 0.1732$), the current optimisation suggests that the model is spectrally suboptimal. This seems to explain why the dynamic Clark model (Vreman et al. Reference Vreman, Geurts and Kuerten1997) is successful at coarse LES grid since the dynamic procedure chooses a correct SGS dissipation dynamically.

Figure 11. The optimised model constants for the Clark-type closure (6.4) for (a) the gradient and (b) the Smagorinsky parts of the closure.

If the filter cutoff scales become small, the gradient part of the Clark model should dissipate three or four times more energy for spectral energy balance, represented by higher $\eta _1$ in figure 11(a) at $\kappa _{cut}{\ell _{k,\infty }} \gtrsim 0.39$. If such $\eta _1$ is used, the optimised Smagorinsky constant is required to be negative at $\kappa _{cut}{\ell _{k,\infty }} \gtrsim 0.8$. Thus, if a Smagorinsky-type model is used with the standard or dynamically determined coefficient, the LES prediction is expected to be too dissipative since the Smagorinsky constant should be positive. Thus, across all filter cutoff wavenumbers, combining the gradient model and the static Smagorinsky model is not the best choice, creating artificial dissipation originating from the model itself. Using the dynamic procedure is useful to alleviate the dissipation, but cannot address the issue completely. This conclusion appears to be consistent to the previous studies where the dynamic Clark model produces much better accurate LES prediction (Vreman et al. Reference Vreman, Geurts and Kuerten1996, Reference Vreman, Geurts and Kuerten1997) since the Smagorinsky constant is adjusted to better achieve a spectral energy balance. Additional a posteriori studies using the optimised (using the constants in table 1) and dynamic Clark models are needed to support the conclusion.

7. A posteriori validation

The SGS model constants optimised a priori to satisfy the spectral energy criterion (2.15) are used to conduct LES a posteriori. A posteriori study is performed in the following way. For a given equivalent grid resolution reported in table 1, the corresponding optimal constant for a model is prescribed instead of its theoretical value (for instance, $C_{S}^2 = 0.1732^2$ for the Smagorinsky-type model). The LES runs are performed at a series of equivalent grid resolution listed in the table. Under the same flow and simulation conditions, additional LES is performed using the standard models, including the static Smagorinsky closure with $\eta _1 = C_{S}^2 = 0.1732^2$ in (6.1), the DSM (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991), the standard gradient with $\eta _1 = 1/12$ in (6.3) and the standard Clark model with $\eta _1 = 1/12$ and $\eta _2 = C_{S}^2 = 0.1732^2$ in (6.4). In addition, LES without an explicit SGS model is performed for comparison. The Vreman-type model is not tested since its optimised constants are similar to those of the Smagorinsky-type closure as shown in figure 8.

Figure 12 compares the Fourier energy spectra obtained a posteriori for the Smagorinsky-type models. On a coarse grid ($32^3$ points), model dependence is observed, whereas finer-grid ($128^3$-point) results are similar among different models since a posteriori SGS energy flux is already very small on finer-resolution grids (see figure 5c). At all grid resolution examined in this study, the optimised Smagorinsky constants give a posteriori prediction close to those using the dynamic model. In particular, the result in figure 12(a) follows the dynamic Smagorinsky prediction of Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019) who used the same flow condition and a numerical scheme of the same spatial accuracy. The energy spectrum obtained for the static Smagorinsky model remains slightly higher (thus, closer to the DNS prediction). This is expected because its model constant $C_{S}^2 = 0.1732^2$ is lower than the dynamically determined constant and the optimal constant found by the proposed algorithm (see figure 8), thus reducing the SGS energy dissipation. However, the $32^3$ point grid does not sufficiently resolve the inertial subrange, making the choice of $C_{S}^2 = 0.1732^2$ questionable. Thus, the close agreement of the static Smagorinsky model is presumably a consequence of model error. On the finer $128^3$-point grid, energy spectra are similar among different models. The optimised model performs nearly the same as the DSM, which is expected because the optimal constant is nearly the same as the dynamically determined model constant (see figure 8). The simulation using the static Smagorinsky model slightly underpredicts spectral energy since the optimal dissipation is obtained $\eta _1 < C_{S}^2 = 0.1732^2$ as reported in figure 8. As pointed out in figure 3(b), a mild energy pileup still exists at small scales if a SGS model is not used.

Figure 12. Time-averaged Fourier energy spectra for a posteriori LES runs using Smagorinsky-type closures on (a) $32^3$ and (b) $128^3$ grid points. In (a), the LES result of Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019) is shown for comparison.

In figure 13, a posteriori LES results on $32^3$-point grid obtained using the gradient-type models are plotted. The energy spectrum of the pure gradient model (6.3) is similar to that of no-model LES, confirming that its SGS dissipation is insufficient and thus the model cannot be used standalone (Clark et al. Reference Clark, Ferziger and Reynolds1979). The standard Clark model with $\eta _1 = 1/12$ and $\eta _2 = C_{S}^2 = 0.1732^2$ in (6.4) predicts an energy spectrum similar to that of the dynamic Smagorinsky. The optimised Clark model using $\eta _1 = 0.0953$ and $\eta _2 = 0.0503$ found by the present optimisation (see table 1) predicts an energy spectrum close to the dynamic model result of Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019). A posteriori LES using the optimised two-parameter Clark models on ${\gtrsim }128^3$-point grid are not stable, and thus the results are not reported. However, this is expected because the corresponding optimal constants for the Smagorinsky part of the mixed model are negative (see table 1 and figure 11b).

Figure 13. Time-averaged Fourier energy spectra for a posteriori LES runs using gradient-type closures on $32^3$ grid points. Bassenne et al. (Reference Bassenne, Esmaily, Livescu, Moin and Urzay2019) used the DSM.

Comparison of the Fourier energy spectra shows that the optimised constants found a priori by the current wavelet framework result in reasonable a posteriori prediction accuracy. For both Smagorinsky- and Clark-type closures, the optimised SGS models behave a posteriori similar to the DSM. Table 1 indicates that the optimised model constants of the Smagorinsky- and Clark-type models (the eddy-viscosity component of the latter, in particular) are close to each other on $32^3$ point grid ($0.0529$ vs $0.0503$). The gradient part of the optimised Clark model takes $\eta _1 = 0.0953$ (see table 1), slightly larger than the standard choice of $\eta _1 = 1/12 = 0.083$ which does not produce sufficient SGS dissipation as figure 13 shows. Thus, the optimised two-parameter Clark model is essentially the optimised Smagorinsky model, explaining observed similarities in figures 12 and 13. Both a priori and a posteriori results suggest that the dynamic procedure works presumably in a way consistent to the spectrally optimal state that this study employs, that is, the unresolved triadic energy transfer is balanced by the modelled SGS dissipation.

8. Optimisation at higher Reynolds numbers

Results so far are obtained for low-Reynolds-number HIT at ${Re}_\lambda = 85$. Although the formulation and algorithm do not make any assumption regarding Reynolds number or inertial subrange, it is useful to demonstrate their applicability to higher Reynolds numbers. In this section, the optimisation framework is applied to the space–time-resolved DNS database of Cardesa, Vela-Martín & Jiménez (Reference Cardesa, Vela-Martín and Jiménez2017) at ${Re}_\lambda = 315$, in which a pseudospectral method is employed to simulate forced incompressible HIT in a triply periodic box using $1024$ Fourier modes per direction. The ratio of the domain width $L = 2{\rm \pi}$ to the nominal Kolmogorov length is $L/{\ell _{k,\infty }} = 1516$, whereas $L/{\ell _{k,\infty }} = 1005$ for HIT at ${Re}_\lambda = 85$ described in § 4. A total of 10 snapshots taken approximately one integral time scale apart are used to optimise the Smagorinsky-like closure (6.1). Figure 14(a) compares the Fourier energy spectra between ${Re}_\lambda = 85$ (corresponding to the DNS result on figure 3) and $315$ (Cardesa et al. Reference Cardesa, Vela-Martín and Jiménez2017). At ${Re}_\lambda = 315$, a substantially extended inertial subrange is observed up to $\kappa {\ell _{k,\infty }} \approx 0.1$. At wavenumbers smaller than $\kappa {\ell _{k,\infty }} \approx 0.1$, spectral energy density is amplified compared with the result for ${Re}_\lambda = 85$ by an order of magnitude at the largest scale, partially attributed to different external forcing techniques used by the two DNS setups. Agreement at $\kappa {\ell _{k,\infty }} \gtrsim 0.1$ is reasonably good.

Figure 14. (a) Time-averaged Fourier energy spectra. (b) Optimised model constants for the one-parameter Smagorinsky-like (6.1) model. Dynamic estimation of the Smagorinsky constants is shown as $\ast$ (green) for ${Re}_\lambda = 85$ (or ${\times }$ (blue) for ${Re}_\lambda = 315$) as a function of the grid-based wavenumber scaled by $2^{-1}$ for ${Re}_\lambda = 85$ (or $2^{-3}$ for ${Re}_\lambda = 315$).

The model optimisation is performed using the DNS data for ${Re}_\lambda = 315$, and spectrally optimal Smagorinsky constants are shown in figure 14(b) and also tabulated in table 2 as a function of the reduced filter cutoff wavenumber. In table 2, the equivalent LES grid resolution is obtained in the same way as table 1. Also shown in figure 14(b) is the optimisation result for ${Re}_\lambda = 85$, the same as the one for the Smagorinsky closure in figure 8. At small filter cutoff scales resolving the inertial subrange, both sets of constants decrease in magnitude at similar rates. For ${Re}_\lambda = 315$, the constant is reduced to $0.0046$ (see table 2), whereas the constant is $0.0184$ for ${Re}_\lambda = 85$ (see table 1). The filter cutoff wavenumbers at which the optimal Smagorinsky constants recover the theoretical prediction (Lilly Reference Lilly1967) are different between the low- and high-Reynolds-number flows at $\kappa _{cut} {\ell _{k,\infty }} = 0.39$ and $0.088$, respectively. This is expected from figure 14(a) since the inertial subranges for the two flows extend up to different wavenumbers, respectively. For ${Re}_\lambda = 315$, the inertial subrange is resolved if the filter cutoff wavenumber is $\kappa _{cut} {\ell _{k,\infty }} \gtrsim 0.1$ (see figure 14a), whereas for ${Re}_\lambda = 85$, $\kappa _{cut} {\ell _{k,\infty }} \gtrsim 0.4$ is expected to resolve inertial subrange motions. If the filter cutoff scales do not resolve the inertial subrange, the optimised Smagorinsky constants of the two flows differ significantly from each other. However, this is expected because figure 14(a) shows that spectral energy contents at those scales are considerably different.

Table 2. The optimised model constants for a range of the filter cutoff wavenumbers. The equivalent LES grid resolution is determined so that the inertial subrange behaviour is matched between a priori and a posteriori results. The theoretical estimation of the Smagorinsky constant in the inertial subrange is $C_{S}^2 = 0.1732^2 = 0.03$.

Figure 14(b) shows that if the nominal LES filter width resolves the inertial subrange, spectrally optimal Smagorinsky constants show little sensitivity to Reynolds number. This is encouraging for extending the optimisation framework to more practical Reynolds numbers. However, since the current comparison is affected by Reynolds number difference, external forcing that sustains HIT, and potential numerical errors associated with finite difference used for HIT at ${Re}_\lambda = 85$, a subsequent study comparing optimisation results using the two numerical methods at the same Reynolds number of HIT forced using the same external forcing technique is necessary.

Using the equivalent grid resolution in table 2, a posteriori LES runs are conducted for the standard and optimised Smagorinsky models, respectively. In addition, LES runs using the standard DSMs are performed on the same equivalent grid resolution. Figure 15 shows the dynamic Smagorinsky constants on the equivalent grid resolution. Due to higher Reynolds number, the dynamic procedure predicts the Smagorinsky constant the same as the theoretical estimation on $128^3$-point grid (as opposed to $64^3$ for ${Re}_\lambda = 85$). The dynamically determined Smagorinsky constants are plotted in figure 14(b) as $\times$ (blue). Agreement between the spectrally optimal (a priori) and dynamically determined (a posteriori) coefficients is better at ${Re}_\lambda = 315$ than ${Re}_\lambda = 85$, presumably due to the extended inertial subrange shown in figure 14(a).

Figure 15. Time evolution of the Smagorinsky constant determined by the dynamic process. The horizontal dotted line corresponds to the theoretical estimation of Lilly (Reference Lilly1967).

Further comparison is made in figure 16 by showing the Fourier energy spectra obtained a posteriori. On $32^3$-point grid (figure 16a), the optimised Smagorinsky constant results in energy spectrum nearly indistinguishable from that of the DSM, consistent with figure 12 obtained for ${Re}_\lambda = 85$. The standard Smagorinsky model using $C_{S}^2 = 0.03$ is less dissipative on the equivalent LES grid resolution of $32^3$ (see table 2), also observed in figure 16(a). At a refined LES grid of $128^3$ points (figure 16b), all a posteriori energy spectra collapse as designed for Smagorinsky-type closures (Pope Reference Pope2001). Reduced spectral energy at the largest scales appears to stem from using the different external forcing technique by Bassenne et al. (Reference Bassenne, Urzay, Park and Moin2016), rather than numerical errors or Reynolds number effects.

Figure 16. Time-averaged Fourier energy spectra for a posteriori LES runs using Smagorinsky-type closures on (a) $32^3$ and (b) $128^3$ grid points.