Optimising subgrid-scale closures for spectral energy transfer in turbulent flows

Abstract

Subgrid-scale (SGS) modelling is formulated using a local transport of spectral kinetic energy estimated by a wavelet multiresolution analysis. Using a spectrally and spatially local decomposition by wavelet, the unresolved inter-scale energy transfer and modelled SGS dissipation are evaluated to enforce explicitly and optimally their balance a priori over a range of large-eddy simulation (LES) filter widths. The formulation determines SGS model constants that optimally describe the spectral energy balance between the resolved and unresolved scales at a given cutoff scale. The formulation is tested for incompressible homogeneous isotropic turbulence (HIT). One-parameter Smagorinsky- and Vreman-type eddy-viscosity closures are optimised for their model constants. The algorithm discovers the theoretical prediction of Lilly (The representation of small-scale turbulence in numerical simulation experiments. In Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences, pp. 195–210) at a filter cutoff scale in the inertial subrange, whereas the discovered constants deviate from the theoretical value at other cutoff scales so that the spectral optimum is achieved. The dynamic Smagorinsky model used a posteriori shows a suboptimal behaviour at filter scales larger than those in the inertial subrange. A two-parameter Clark-type closure model is optimised. The optimised constants provide evidence that the nonlinear gradient model of Clark et al. (J. Fluid Mech., vol. 91, issue 1, 1979, pp. 1–16) is prone to numerical instability due to its model form, and combining the pure gradient model with a dissipative model such as the classic Smagorinsky model enhances numerical stability but the standard mixed model is not optimal in terms of spectral energy transfer. A posteriori analysis shows that the optimised SGS models produce accurate LES results.

1. Introduction

1.1. Background

One of the fundamental questions in large-eddy simulation (LES) concerns its closure model. With some exceptions such as King, Hamlington & Dahm (2016), most of the algebraic subgrid-scale (SGS) closures assume a priori their specific functional forms that depend on resolved-scale quantities. In addition, the closure models are parameterised by a length scale frequently chosen to be proportional to the LES grid width. See Piomelli, Rouhi & Geurts (2015) for a grid-independent choice of a model length scale and Bose, Moin & You (2010) for examining the grid dependence of LES solutions using an explicit filter. A presumed analytical closure model generally behaves well in the regimes where the assumption remains valid. Examples include the Smagorinsky model (Smagorinsky 1963) designed to work when a LES filter width is in the inertial subrange and the closure approach developed by Bardina, Ferziger & Reynolds (1980) requiring a scale similarity. For additional discussions, see Pope (2001). However, the validity of the modelling assumption is generally not warranted, and nontrivial modelling errors are introduced to LES prediction. This also applies to the commonly used eddy-viscosity models based on the Boussinesq hypothesis, since the SGS stress is not necessarily aligned with the resolved strain-rate tensor (Meneveau & Katz 2000). Furthermore, a presence of the LES grid spacing in a closure model (approximating the model filter width) makes it complicated to interpret the grid-convergence results of LES prediction.

The model errors are not necessarily reduced as the LES grid is refined or a more accurate numerical discretisation is employed. Thus, they could have leading-order effects on LES prediction in a grid-resolved limit or in a laminar regime, unless the SGS stress is designed to vanish under such conditions. Moreover, turbulent fluid motions interacting with multiphysical effects bring additional modelling challenges. For instance, turbulent premixed flames at moderate Karlovitz numbers (${Ka} > 1$) have a sufficient time for the thermal effects to interact with the Kolmogorov eddies. In such regimes, a SGS closure may potentially require to incorporate the effects of the combustion-induced energy backscatter (O'Brien et al. 2017; Sabelnikov et al. 2023), which is difficult to model via the resolved strain rate alone. Similarly, a priori study demonstrates that a modelled SGS stress for a two-way coupled LES of particle-laden turbulence of the mass loading of unity should include the inter-phase coupling effects if the particle inertia is significant (Nabavi, Di Renzo & Kim 2022). For such configurations, the specific functional form of the modelled SGS stress is not usually known a priori, and related theoretical studies are scarce, making the modelling study difficult.

A number of previous studies examined the model errors of the SGS closure and proposed alternative formulations different from the prevailing eddy-viscosity-type models, some of which are described in Vreman, Geurts & Kuerten (1997), Meneveau & Katz (2000) and Pope (2001). Lund & Novikov (1992) combined the symmetric and anti-symmetric components of the velocity gradient tensor and obtained a complete set of their products. Silvis, Remmerswaal & Verstappen (2017) investigated the desired mathematical and physical properties of SGS model and developed a framework for understanding and analysing existing SGS models, by which a new model can be developed. The autonomic closure (King et al. 2016) does not make any assumption about SGS model form and determines its general non-parametric relation dynamically. Langford & Moser (1999) developed a SGS modelling framework by which an approximation to an ‘ideal’ closure that accurately determines one-time, multipoint statistics is generated from two-point statistics using stochastic estimation (Adrian 1990), later tested a posteriori by the same authors (Langford & Moser 2004). Algorithms based on, for example, the neural networks could provide some breakthroughs to SGS modelling (Duraisamy, Iaccarino & Xiao 2019; Zhou et al. 2019; Sarghini, De Felice & Santini 2003; Xie, Yuan & Wang 2020; Vinuesa & Brunton 2022; Xu et al. 2023).

In addition, several previous studies demonstrated that spectral energy transfer can be leveraged in SGS modelling. The recent series of work done by Domaradzki (2021a,b, 2022) formulated LES in terms of the Fourier spectral energy transfer. An effective, spectral eddy viscosity is introduced with an assumption about the inertial subrange and a model spectrum. It is noteworthy that the formulation is capable of minimising modelling inputs (thus, nearly autonomous) (Domaradzki 2022). Although spectral energy flux is not incorporated directly, the minimum dissipation model (Rozema et al. 2015) is built upon an idea involving spectral kinetic energy. A similar argument is used to develop the dynamic global models (Park et al. 2006).

1.2. Proposed SGS modelling framework

This study develops a framework in which a SGS closure model that satisfies a certain energy criterion is extracted a priori from the corresponding direct numerical simulation (DNS) data. Similar to many previous studies, it is assumed that the SGS stress can be modelled algebraically using the resolved-scale quantities in the LES equations. However, its specific functional form is not presumed but determined using the DNS data a priori. This procedure implies that grid-independent LES is feasible. A second and more important assumption is that a balance of spectral energy exchange between the unresolved triadic interactions and modelled SGS dissipation is both necessary and sufficient for stable and accurate LES prediction. This is consistent with the fundamental concept of LES where a SGS closure should dissipate the correct amount of the resolved kinetic energy that should be transferred via the triadic interactions to SGS motions. It is noted that such energy transfer should hold not only statistically but also locally in both scale and space, since the modelled SGS dissipation is not just a statistical quantity if the LES equations are time integrated in the physical space.

The second assumption is justified for turbulence in equilibrium by examining the transport equation of the Fourier spectral energy in the LES context (Pope 2001)

(1.1)

where $\hat {E}$ is the Fourier energy spectrum, $\kappa$ is the Fourier wavenumber, $\kappa _{cut}$ is the filter cutoff wavenumber (thus, $\kappa < \kappa _{cut}$ in LES) and $\nu$ is the kinematic viscosity. In (

1.1

), the viscous dissipation is always negative and local in $\kappa$ (thus, no inter-scale interactions). Although $\hat {T}_{f}$ (second term on the right-hand side) involves the triadic interactions, they are resolved on the LES grid. Only $\hat {T}_{r}$ (third term on the right-hand side) is unclosed, requiring SGS modelling. Integrating (

1.1

) over the resolved scales $\kappa < \kappa _{cut}$,

(1.2)

it can be seen that a balance among the spectral energy fluxes must be achieved in LES by correctly modelling $\hat {T}_{r}$. In general, existing closure models are not designed to satisfy such energy balance independently of grid resolution. Without it, spectral energy transfer could be altered in a way dependent on model and grid, providing excessive or insufficient dissipation to the resolved scale motions.

Although spectral energy balance is elucidated using the Fourier theory, it is not straightforward how such theory can be utilised to enforce a condition similar to (1.1) or (1.2). Even though spectral kinetic energy and spectral energy fluxes can be computed using the Fourier transformation at an excellent spectral resolution (Batchelor 1953), the domain has to be periodic and all spatial information is lost in evaluating the transformation (that is, no spatial resolution). Spatially local interactions are often encountered in turbulence with multiphysics (and even isotropic turbulence) and require to combine a large number of Fourier modes. Thus, the Fourier formulation is not ideal for modelling the SGS stress in the physical space, although it is an excellent statistical analysis tool for ideal turbulent flows such as homogeneous turbulence.

The modelling challenges discussed previously can be addressed to certain extent using different techniques. Windowed Fourier transformation (in particular, the Gabor transformation) and coarse graining (Aluie 2013), among several others, can be used as alternatives for their finite spatial resolution. Using the Reynolds decomposition and the fluctuating turbulence kinetic energy (TKE), local energy transfer across a prescribed cutoff scale can be computed (Meneveau & Katz 2000). Spatial information is made available by formulation. However, the procedure is statistical and suited more for analysis than modelling.

Orthogonal wavelet transformation, in particular the multiresolution framework (Mallat 1989), is one of the techniques by which the known limitations of the Fourier analysis can be addressed for the SGS modelling purposes. The fundamental differences between the wavelet framework and the standard (or windowed) Fourier transformation are well documented in the literature (Mallat 1999). By design, wavelet provides a spectrally and spatially local decomposition of turbulent motions in a physically consistent and computationally efficient fashion. Moreover, it can be formulated in a way similar to the Fourier transformation to examine the spectral energy transport of turbulent flows (Meneveau 1991). A common practice is that the spatial distribution of spectral kinetic energy and spectral energy fluxes obtained by a wavelet transformation is interpreted as statistical variability or intermittency (Meneveau 1991).

This study proposes that the spectral and spatial information made available by wavelet can be exploited for optimising or discovering SGS closure terms a priori. For existing closure models such as the Smagorinsky model, the proposed optimisation procedure generates scale-dependent, optimal model constants in a sense of spectral energy balance between the resolved and unresolved-scale motions. In doing so, the behaviour of SGS models as a nominal LES grid is refined is studied systematically in a priori context. In addition, the procedure is useful to examine the limiting behaviour of closure models at DNS-like grid resolution. Consequently, it is argued that wavelet is more than a statistical tool and can be crafted as a useful technique for constructing a turbulence modelling framework including configurations where turbulent momentum transfer interacts in two ways with multiphysics.

This paper is structured as follows. The wavelet multiresolution analysis (WMRA) technique is formulated in the LES context in § 2. The algorithm for optimising a SGS closure is outlined in § 3. The proposed formulation is applied to the DNS data of isotropic turbulence, and the simulation results of DNS and the corresponding LES using the standard dynamic Smagorinsky model (DSM) are found in § 4. The optimisation results are reported, and related discussions are provided in § 6. A posteriori validations using the optimised constants follow in § 7. The paper concludes in § 9 proposing relevant future works.

2. WMRA framework for SGS modelling

The present study employs a WMRA framework originally developed by Meneveau (1991) for studying the spectral energy transport and later extended to turbulent channel flows (Dunn & Morrison 2003), turbulent combustion (Kim et al. 2018), droplet-laden turbulence (Freund & Ferrante 2019) and particle-laden turbulence (Nabavi et al. 2022). Compared with the Fourier analysis of spectral energy transfer, wavelet offers several distinct advantages (Meneveau 1991; Kim et al. 2018; Nabavi et al. 2022). 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 1991; Kim et al. 2018; Nabavi et al. 2022), 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. (2022).

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 1999; Meneveau 1991) is employed, and the complete details are found in Nabavi et al. (2022). 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. (2022) and Meneveau (1991), 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. (2022). 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 2001), 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 (

2.5

), the scale decomposition proposed by Meneveau (1991) 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 2004; De Stefano & Vasilyev 2012; De Stefano, Dymkoski & Vasilyev 2022), 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 1991). 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 2001). 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 2001), spectral enrichment (Bassenne et al. 2019; Ghate & Lele 2020) and autonomic closure (King et al. 2016). 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 2000). 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 2000), 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 (

2.14

) 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 1991; Nabavi et al. 2022), 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. (2017) 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 (1999). A primary conceptual difference is that Langford & Moser (1999) 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 (1999) 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 (1999) used stochastic estimation (Adrian 1990) 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 (1999), 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 1963), in particular, its dynamic variant (Germano et al. 1991) 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. 1991), the Smagorinsky constant $C$ is not prescribed but determined using the dynamic procedure. Computational work is done using the NGA solver (Desjardins et al. 2008) 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 (2018). To sustain turbulence, an external forcing (Bassenne et al. 2016) 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 1987). 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 (1967).

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. (2016). In the LES context, Lilly (1967) 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. (2018). For the $32^3$-point grid, the LES result of Bassenne et al. (2019) is also compared. Under the same flow condition and using the same SGS closure, the energy spectra of Bassenne et al. (2019) 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 1999). 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. (2018) (DNS on $512^3$ grid points) and Bassenne et al. (2019) (DSM on $32^3$ grid points) are shown. In (b), the Fourier enstrophy spectrum of Bassenne et al. (2018) (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 2000), 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. 2016), 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 1991). 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 2004). 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 2004), 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 1979), the unknown constant is expected to be close to $1/12$ at a DNS-like grid resolution (Pope 2001). 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 1967). 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 2001), 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 2010) 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. 1979) 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. 1979) 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 1996) or its dynamic variant is used (Vreman et al. 1996, 1997). 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. 1979).

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. (1979) and Vreman et al. (1996) 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. 1979). 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. (1979), the eddy viscosity is arbitrarily included to ensure numerical stability, which is however reported to cause excessive dissipation in transitional regimes (Vreman et al. 1996). 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. 1997) 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. 1996, 1997) 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. 1991), 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. (2019) 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. (2019) 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. 1979). 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. (2019). 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. (2019) 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 (2017) 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. 2017). 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 1967) 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 (1967).

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 2001). Reduced spectral energy at the largest scales appears to stem from using the different external forcing technique by Bassenne et al. (2016), 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.

9. Conclusions

LES has been an attractive approach simulating high-Reynolds-number turbulent flows in terms of computational efficiency and accuracy compared with the traditional, low-fidelity simulation techniques. Although it is still being debated whether LES will replace completely the Reynolds-averaged Navier–Stokes formulations in industry design processes in the foreseeable future, the technology is certainly useful and widely used in academia and some parts of the industry as a high-fidelity prediction tool. Nevertheless, many of the technological aspects of LES as well as its fundamental modelling concept are not entirely free from criticism (Pope 2004).

Numerical errors associated with LES have been studied extensively, including the effects of filter commutativity (Vasilyev, Lund & Moin 1998; Haselbacher & Vasilyev 2003; van der Bos & Geurts 2005; Klein & Germano 2020), discretisation error (Ghosal 1996; Geurts 2009; Viré & Knaepen 2009), truncation error (Kravchenko & Moin 1997), numerical dispersion and dissipation (Mittal & Moin 1997; Yalla, Oliver & Moser 2021), the averaging effects for the dynamic models (Lilly 1992; Meneveau, Lund & Cabot 1996), grid anisotropy (Rozema et al. 2015; Schumann, Toosi & Larsson 2020) and grid-independent LES (Bose et al. 2010). Some of those errors are intrinsically tied with SGS model errors, which have been also investigated to some extent (Vreman et al. 1997; Meneveau & Katz 2000; Pope 2001; Silvis et al. 2017). However, the modelling errors are much less straightforward to assess, and a systematic approach by which different SGS models are compared is lacking.

In practice, a SGS closure is selected a priori considering the state of turbulent flow to be simulated. Such choice is often informed by corresponding a priori studies. Since a SGS model is designed typically with certain assumptions about turbulence (such as flow regime, scale similarity and filter scale in the inertial subrange), such practice becomes less reliable if turbulence is in a strongly non-equilibrium state or non-turbulent regions are found in the LES domain. One of the consequences is that the model errors become significant, and may not disappear (or may even increase) if the LES grid is refined. In such cases, it is not straightforward to distinguish model errors from those originating from different sources. Other approaches such as the autonomic closure (King et al. 2016) or machine learning (Duraisamy et al. 2019; Vinuesa & Brunton 2022) can be incorporated to address the closure problem dynamically and adaptively. However, fundamental and technical questions regarding the physical consistency of their LES prediction and computational efficiency remain unaddressed. In particular, a dynamic determination of correct SGS stresses via, for instance, a data-driven procedure seems technically feasible and useful, but its predictability, uniqueness and consistency to the filtered governing equations are not well understood as of yet.

This study proposes a framework by which DNS data are leveraged to inform the selection of an optimal SGS closure model a priori. It is assumed that an algebraic SGS model describes the behaviour of LES-predicted solutions over a range of filter cutoff scales. An optimal SGS model is selected by using a WMRA and by weakly enforcing a spectral energy balance between the filtered and unfiltered scales. Using WMRA, spectral energy fluxes for inter-scale energy transfer and modelled SGS energy dissipation are computed locally in scale and space. Constraining a SGS model with the energy balance is rooted in the fundamental idea of LES that a SGS model should dissipate a correct amount of energy transferred by the resolved scale of fluid motions towards the unresolved motions.

This article has presented the formulation and validated it by demonstrating that optimal model constants are found for prescribed SGS closure models for forced HIT. Two eddy-viscosity models (static Smagorinsky and static Vreman models), the nonlinear gradient closure (Clark et al. 1979) and the Clark model (Clark et al. 1979) have been tested. The optimisation results, even if they are obtained a priori, are consistent both quantitatively and qualitatively with the theoretical estimation, the analytical limiting behaviour and a posteriori LES prediction using the DSM (Germano et al. 1991). A posteriori tests performed using the optimised model constants confirm such findings. The formulation can be used for SGS models having more than one undetermined constant.

Despite some encouraging results, the proposed framework requires non-trivial improvements in both fundamental and technical aspects. The presumed availability of a DNS database makes the framework less straightforward to extend to higher-Reynolds-number flows in which DNS is not always feasible. One possible solution is to leverage the high-resolution 3-D experimental data of high-Reynolds-number turbulent flows. By using experimental data, potential artifacts stemming from physical models often used for multiphysical turbulence simulation can be circumvented.

These difficulties can also be addressed by optimising SGS model constants for several low- and intermediate-Reynolds-number turbulent flows and extrapolating the models statistically or with an aid of machine learning algorithms to flow regimes where DNS is not tractable. Such extrapolation is common in conventional flow modelling. Early pioneering works in LES modelling focused on (primarily due to limited computing resources) low-Reynolds-number canonical configurations of turbulent flows such as HIT, free shear flow and channel flows (Lesieur & Metais 1996; Meneveau & Katz 2000). Insight gained from the analysis and modelling of lower-Reynolds-number turbulence was extended and applied to realistic turbulent flows at high Reynolds numbers. It is expected that the current formulation based on spectral energy transfer is presumably as effective at high Reynolds numbers where inertial subrange is well developed (as demonstrated for a Smagorinsky-type model in § 8) and the optimisation criterion (2.15) is justified provided that LES filter width is not excessively small to resolve the dissipation range.

A specific extrapolation procedure is presumably dependent on optimisation results (such as their statistical fluctuations and asymptotic behaviour) at lower Reynolds numbers and not known a priori. However, it seems feasible that a practice of extrapolating wind-tunnel or subscale flight test data to full-scale flight conditions in terms of Reynolds number can be one of the candidates (Pettersson & Rizzi 2008; Kulkarni et al. 2022)

The extrapolation efforts to higher Reynolds numbers can be streamlined by employing recent advances in data-driven modelling. For instance, Lozano-Durán & Bae (2023) combined canonical, low-to-intermediate Reynolds number, building-block flows to train artificial neural networks and improve LES wall-modelling in realistic turbulent flows around complex geometry. Examples of other recent works leveraging machine learning algorithms for novel LES modelling are discussed by, for example, Vinuesa & Brunton (2022).

In addition, the optimisation procedure constrains a SGS model with the spectral energy balance. Although the energy criterion seems consistent to the fundamental concept of LES, it concerns primarily the numerical stability of LES prediction (rather than prediction accuracy). This suggests that there could be another criterion that should be enforced simultaneously on a SGS model, and the energy criterion employed in this study is one of the necessary conditions for the true SGS closure as pointed out by several previous studies (Meneveau 1994; Silvis et al. 2017).

Although this study uses incompressible HIT at a low Reynolds number to demonstrate the feasibility of the modelling framework, the framework is not limited to such condition. The current formulation is developed based on the spectral transport of TKE as shown in (2.14) and an energy balance across a grid-cutoff scale (2.15), neither of which depend on flow regimes. For instance, the unresolved inter-scale energy flux (2.12) can be formulated and evaluated for compressible flows. An existing algebraic SGS model for compressible turbulence (for instance, Moin et al. (1991)) can be used to compute ${{\check {\bar {\mathcal {T}}}\vphantom {}^{(s,d)}_{SGS}}}$.

Technical improvements are desired regarding the use of wavelet transformation. It is well known that the analysis results using wavelet transformation depends on the choice of its basis (Mallat 1999). Thus, it is likely that the optimisation results would depend on the wavelet basis used by this study. For instance, a factor of $2^{-1}$ rescaling is made on the wavenumber so that the inertial subrange behaviour is matched between the current optimisation and the dynamic Smagorinsky LES. The factor depends presumably on wavelet basis or the spectral accuracy of numerical discretisation. Thus, several commonly used wavelet bases should be tested, and a priori optimisation results should be accompanied by a posteriori validation.

In (2.8), the resolved-scale motions are estimated by performing a scale-cutoff filtering in the wavelet space (Meneveau 1991). Although this choice provides a priori estimation of the filtered velocity, it is well known in the Fourier context that a scale-sharp filter does not provide an accurate estimation of the grid-resolved fluid motions, which is presumably the case for wavelet as well. It is believed that a better estimation of the resolved-scale motions improves the optimisation results, which is left for a future work.

Although wavelet is used as a key technique for constructing the proposed modelling framework, the choice is not unique. In other words, a decomposition technique that is not based on wavelet transformation can be incorporated in the current framework to evaluate the spectral energy fluxes. However, it is still required that such technique provides a finite resolution in scale and in space simultaneously.