2.1. Large eddy simulation

In LES, the flow variables are filtered using the following operation (Pope Reference Pope2000):

(2.1)$$\begin{gather} \bar{\phi}(\boldsymbol{x},t)=\int {\bar{G}}(\boldsymbol{r})\phi(\boldsymbol{x-r},t)\,{\rm d}\boldsymbol{r}, \end{gather}$$

(2.2)$$\begin{gather}\int {\bar{G}}(\boldsymbol{r})\,{\rm d}\boldsymbol{r}=1, \end{gather}$$

where $\bar \phi$ is a filtered flow variable, $\bar {G}$ is a filter function, $\boldsymbol {x}$ and $\boldsymbol {r}$ denote position vectors, and $t$ is time. The spatially filtered continuity and Navier–Stokes equations for incompressible flows are

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

(2.4)$$\begin{gather}\frac{\partial {\bar u}_i}{\partial t}+\frac{\partial {\bar u}_i {\bar u}_j}{\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}, \end{gather}$$

where $x_i$ are the coordinates ($x,y,z$), $u_i$ are the corresponding velocities ($u,v,w$), and $\rho$ and $\nu$ are the fluid density and kinematic viscosity, respectively (Pope Reference Pope2000). The effect of the SGS eddy motions is modelled as the anisotropic part of SGS stress, ${\tau }^r_{ij}$:

(2.5)$$\begin{gather} {\tau}_{ij} = \overline{{u}_i{ u}_j}-{\bar u}_i{\bar u}_j, \end{gather}$$

(2.6)$$\begin{gather}{\tau}^r_{ij} = {\tau}_{ij}-\tfrac{1}{3}{\tau}_{kk}\delta_{ij}, \end{gather}$$

where $\delta _{ij}$ is the Kronecker delta. The filtered pressure $\bar p^{*}$ includes the isotropic components of SGS stress:

(2.7)\begin{equation} \bar p^{*} = \bar p + \tfrac{1}{3} \rho \tau_{kk}. \end{equation}

In the case of HIT, the filtered continuity and Navier–Stokes equations can be transformed in a spectral space. The pseudo-spectral method is used for spatial discretization, and the zero-padding method, augmented with the 3/2 rule, is applied to control aliasing errors. For temporal integration, a third-order Runge–Kutta method is used for the convection term, and the second-order Crank–Nicolson method is applied to the diffusion term.

To evaluate the performance of an NN-based SGS model relative to those of traditional SGS models, LESs with the constant Smagorinsky model (CSM, Smagorinsky Reference Smagorinsky1963), dynamic Smagorinsky model (DSM, Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992), gradient model (GM, Clark et al. Reference Clark, Ferziger and Reynolds1979) and dynamic mixed model (DMM, Anderson & Meneveau Reference Anderson and Meneveau1999) are carried out. The Smagorinsky models determine the anisotropic component of SGS stress from

(2.8)\begin{equation} {\tau}^r_{ij} ={-}2(C_s \bar \varDelta)^2 |\bar S| \bar S_{ij}, \end{equation}

where $C_s$ is the Smagorinsky model coefficient, $\bar \varDelta$ is the grid size, $\bar S_{ij}=(1/2)(\partial \bar u_i/\partial {x_j}+\partial \bar u_j/\partial {x_i})$ and $|\bar S|=\sqrt {2\bar S_{ij} \bar S_{ij}}$. For CSM, $C_{s}$ is a fixed value of 0.17 (Pope Reference Pope2000). For DSM, the Smagorinsky model coefficient is obtained as

(2.9)\begin{gather} C_s^2 = \max [\langle L_{ij}^{r} M_{ij} \rangle_h/ \langle M_{ij}M_{ij} \rangle_h, 0], \end{gather}

(2.10)$$\begin{gather} L_{ij}=\widetilde{\bar u_i \bar u_j} - {\widetilde{\bar u}_i} {\widetilde{\bar u}_j}, \end{gather}$$

(2.11)$$\begin{gather}{L}^r_{ij} = {L}_{ij}-\tfrac{1}{3}{L}_{kk}\delta_{ij}, \end{gather}$$

(2.12)$$\begin{gather}M_{ij}={-}2\tilde{\varDelta}^2 |\tilde{\bar S}| \tilde{\bar S}_{ij} + 2 {\bar \varDelta}^2 \widetilde{|\bar S|\bar S_{ij}}, \end{gather}$$

where $\overline {({\cdot })}$ and $\widetilde {({\cdot })}$ correspond to the grid- and test-filtering operations, respectively, and $\tilde {\varDelta } = 2\bar \varDelta$. The notation $\langle {{\cdot }} \rangle _h$ denotes an instantaneous averaging over homogeneous directions. The SGS stress in GM is derived through the application of Taylor expansions of the filtered velocity (Clark et al. Reference Clark, Ferziger and Reynolds1979; Vreman et al. Reference Vreman, Geurts and Kuerten1996) as

(2.13)$$\begin{gather} \tau_{ij}^r=\frac{1}{12} \left( \beta_{ij} - \frac{1}{3}\beta_{kk}\delta_{ij} \right), \end{gather}$$

(2.14)$$\begin{gather}\beta_{ij}=\sum_{l} {\bar \varDelta_l^2} \frac{\partial \bar u_i}{\partial{x}_l} \frac{\partial \bar u_j}{\partial{x}_l}, \end{gather}$$

where $\bar \varDelta _l$ is the grid size in the $x_l$ direction. The DMM (Anderson & Meneveau Reference Anderson and Meneveau1999) is a combination of CSM and GM with the model coefficients $C_s$ and $C_g$ that are dynamically determined:

(2.15)\begin{equation} {\tau}^r_{ij} ={-}2(C_s \bar \varDelta)^2 |\bar S| \bar S_{ij}+C_{g} ( \beta_{ij} - \tfrac{1}{3}\beta_{kk}\delta_{ij} ), \end{equation}

where

(2.16)$$\begin{gather} C_s^2 = \max \left[ \frac{\langle L_{ij}^r M_{ij}\rangle_h\langle N_{ij} N_{ij}\rangle_h-\langle L_{ij}^r N_{ij}\rangle_h\langle M_{ij} N_{ij}\rangle_h}{\langle M_{ij} M_{ij}\rangle_h\langle N_{ij} N_{ij}\rangle_h-\langle M_{ij} N_{ij}\rangle_h \langle M_{ij} N_{ij}\rangle_h}, 0 \right], \end{gather}$$

(2.17)$$\begin{gather}C_g = \frac{\langle L_{ij}^r N_{ij}\rangle_h\langle M_{ij} M_{ij}\rangle_h-\langle L_{ij}^r M_{ij}\rangle_h\langle M_{ij} N_{ij}\rangle_h}{\langle M_{ij} M_{ij}\rangle_h\langle N_{ij} N_{ij}\rangle_h-\langle M_{ij} N_{ij}\rangle_h \langle M_{ij} N_{ij}\rangle_h}, \end{gather}$$

(2.18)$$\begin{gather}N_{ij}=\left( B_{ij} - \frac{1}{3} B_{kk}\delta_{ij} \right)-\left( \widetilde \beta_{ij} - \frac{1}{3}\widetilde\beta_{kk}\delta_{ij} \right), \end{gather}$$

(2.19)$$\begin{gather}B_{ij}=\sum_{l} {\widetilde \varDelta_l^2} \frac{\partial {\widetilde{\bar u}}_i}{\partial{x}_l} \frac{\partial {\widetilde{\bar u}}_j}{\partial{x}_l}. \end{gather}$$

2.2. NN-based SGS model

So far, most NN-based SGS models have optimized the weights and biases of the NNs, despite using different deep learning techniques (see § 1). For example, Park & Choi (Reference Park and Choi2021) implemented an NN consisting of two hidden layers and 128 neurons per hidden layer to compute six components of the SGS stress:

(2.20)\begin{align} \left. \begin{array}{ll} h_i^{(1)}= \bar q_i & (i=1,2,\ldots,N_q), \\ \displaystyle h_j^{(2)}=\max \left[0,\gamma_{j}^{(2)} \left( \sum_{i=1}^{N_q}W^{(1)(2)}_{ij}h^{(1)}_i + b_{j}^{(2)} - \mu_{j}^{(2)} \right) / \sigma_{j}^{(2)} + \beta_{j}^{(2)} \right] & (j=1,2,\ldots ,128), \\ \displaystyle h_k^{(3)}=\max \left[0,\gamma_{k}^{(3)} \left( \sum_{j=1}^{128}W^{(2)(3)}_{jk}h^{(2)}_j + b_{k}^{(3)} - \mu_{k}^{(3)} \right) / \sigma_{k}^{(3)} + \beta_{k}^{(3)}\right] & (k=1,2,\ldots ,128),\\ \displaystyle h_l^{(4)}=s_l=\sum_{k=1}^{128}W^{(3)(4)}_{kl}h^{(3)}_k + b_{l}^{(4)} & (l=1,2,\ldots ,6). \end{array}\right\} \end{align}

Here, $\bar {\boldsymbol {q}}$ is the grid-filtered input, $N_q$ is the number of the input components, $\boldsymbol {W}^{(m)(m+1)}$ is the weight matrix between $m$th and $(m+1)$th layers, $\boldsymbol {b}^{(m)}$ is the bias vector of the $m$th layer, $\boldsymbol {s}$ is the output (six components of the SGS stress), and $\boldsymbol {\gamma }^{(m)}$, $\boldsymbol {\mu }^{(m)}$, $\boldsymbol {\sigma }^{(m)}$ and $\boldsymbol {\beta }^{(m)}$ are the parameters for batch normalization (Ioffe & Szegedy Reference Ioffe and Szegedy2015). During an NN training process, the parameters ($\boldsymbol {W}^{(m)(m+1)}$, $\boldsymbol {b}^{(m)}$, $\boldsymbol {\gamma }^{(m)}$, $\boldsymbol {\mu }^{(m)}$, $\boldsymbol {\sigma }^{(m)}$ and $\boldsymbol {\beta }^{(m)}$) are optimized to minimize the training error.

In the present study, we construct an NN-based SGS model using a dual NN architecture (figure 1), where the output of one NN (NN$_{N}$) is the SGS normal stresses $(\tau _{11}^r,\tau _{22}^r,\tau _{33}^r)$ and that of the other (NN$_{S}$) is the SGS shear stresses $(\tau _{12}^r,\tau _{13}^r,\tau _{23}^r)$. The reason for using two NNs is that the ranges of the normal and shear SGS stresses are quite different from each other, and thus separate treatments of these SGS stresses increase the prediction capability of the present NN-based SGS model. Each of the present NNs consists of two hidden layers and 64 neurons per hidden layer, and the output of the $m$th layer, $\boldsymbol {h}^{(m)}$, is computed as

(2.21)\begin{equation} \left.\begin{array}{ll} h_i^{(1)}= \bar q_i & (i=1,2,\ldots,9), \\ \displaystyle h_j^{(2)}=\max[0.02r_j^{(2)},r_j^{(2)}], \quad r_j^{(2)}=\sum_{i=1}^{9}W^{(1)(2)}_{ij}h^{(1)}_i & (j=1,2,\ldots ,64), \\ \displaystyle h_k^{(3)}=\max [0.02r_k^{(3)},r_k^{(3)} ], \quad r_k^{(3)}=\sum_{j=1}^{64}W^{(2)(3)}_{jk}h^{(2)}_j & (k=1,2,\ldots ,64), \\ \displaystyle h_l^{(4)}=s_l=\sum_{k=1}^{64}W^{(3)(4)}_{kl}h^{(3)}_k & (l=1,2,3). \end{array} \right\} \end{equation}

Note that we remove the bias ($\boldsymbol {b}$) and batch normalization parameters ($\boldsymbol {\gamma }$, $\boldsymbol {\mu }$, $\boldsymbol {\sigma }$ and $\boldsymbol {\beta }$). Also, we use the leaky rectified linear unit function (leaky ReLU) as an activation function: $f(x) = \max [ax, x]$, which was proposed by Maas, Hannun & Ng (Reference Maas, Hannun and Ng2013) to overcome the vanishing gradient problem (Bengio, Simard & Frasconi Reference Bengio, Simard and Frasconi1994). The negative slope is determined to be $0.01 \leq a \leq 0.2$ (Xu et al. Reference Xu, Wang, Chen and Li2015), and we choose $a=0.02$. Note that one of the reasons for introducing the bias $\boldsymbol {b}$ is to produce non-zero output even if the input is zero (Goodfellow, Bengio & Courville Reference Goodfellow, Bengio and Courville2016). For the present flow problem, the bias is not necessarily needed because the output should be zero when the input is zero. The effects of removing the bias and batch normalization parameters and replacing ReLU with leaky ReLU on the results of a priori and a posteriori tests are discussed in Appendix A, where we show that setting $\boldsymbol {b} = \boldsymbol {\mu } = \boldsymbol {\beta } = 0, ~\boldsymbol {\gamma } = \boldsymbol {\sigma } =1$ and replacing ReLU with leaky ReLU do not deteriorate the numerical solutions for the present flow.

Figure 1. Schematic diagram of the present NN-based SGS model that consists of a dual NN: NN$_{N}$ and NN$_{S}$ predict the SGS normal and shear stresses, respectively.

Without the bias and batch normalization parameters and with the use of leaky ReLU, the present NN complies with the following conditions:

(2.22)$$\begin{gather} \boldsymbol{{\mathrm{NN}}} (c\boldsymbol{A}) = c \boldsymbol{{\mathrm{NN}}} (\boldsymbol{A}) \quad \text{only if} \ c \geq 0, \end{gather}$$

(2.23)$$\begin{gather}\boldsymbol{{\mathrm{NN}}} (\boldsymbol{A}+\boldsymbol{B}) \ne \boldsymbol{{\mathrm{NN}}} (\boldsymbol{A}) + \boldsymbol{{\mathrm{NN}}} (\boldsymbol{B}), \end{gather}$$

where $c$ is a positive scalar, and $\boldsymbol {A}$ and $\boldsymbol {B}$ are arbitrary tensors ($\boldsymbol {A} \ne \boldsymbol {B}$). For example, for $c=2$ and $\boldsymbol {B} = 2 \boldsymbol {A}$, $\boldsymbol {{\mathrm {NN}}} (c\boldsymbol {A}) = c \boldsymbol {{\mathrm {NN}}} (\boldsymbol {A})$ and $\boldsymbol {{\mathrm {NN}}} (\boldsymbol {A}+\boldsymbol {B}) = \boldsymbol {{\mathrm {NN}}} (\boldsymbol {A}) + \boldsymbol {{\mathrm {NN}}} (\boldsymbol {B})$, but for $c=2$ and $\boldsymbol {B} = - 2 \boldsymbol {A}$, $\boldsymbol {{\mathrm {NN}}} (c\boldsymbol {A}) = c \boldsymbol {{\mathrm {NN}}} (\boldsymbol {A})$ and $\boldsymbol {{\mathrm {NN}}} (\boldsymbol {A}+\boldsymbol {B}) \ne \boldsymbol {{\mathrm {NN}}} (\boldsymbol {A}) + \boldsymbol {{\mathrm {NN}}} (\boldsymbol {B})$ due to the leaky ReLU used. Thus, the present NN still retains its nonlinearity, making it suitable for nonlinear regression between the SGS stresses and local flow variables. Equation (2.22) allows us to apply an NN trained from one flow to another flow. Let us show an example for turbulent channel flow and flow over a circular cylinder. Assume that the input and output variables are $\bar \varDelta ^2 \vert \bar \alpha \vert \bar \alpha _{ij}$ and $\tau ^r_{ij}$ for these two flows, where $\bar \varDelta = (\bar \varDelta _x \bar \varDelta _y \bar \varDelta _z)^{1/3}$, $\bar \alpha _{ij}=\partial \bar u_i / \partial x_j$ and $|\bar \alpha |=\sqrt {\bar \alpha _{ij}\bar \alpha _{ij}}$. An NN is constructed from turbulent channel flow with input and output variables normalized by the wall-shear velocity $u_\tau$: i.e. $\tau ^r_{ij} / u_\tau ^2 = \boldsymbol {{\mathrm {NN}}} ( \bar \varDelta ^2 \vert \bar \alpha \vert \bar \alpha _{ij} / u_\tau ^2 )$. Then, this trained NN can be used for flow over a circular cylinder with input and output variables normalized by the free-stream velocity $u_\infty$ from the following relation:

(2.24)\begin{align} \frac{\tau^r_{ij}}{u_\infty^2} = \frac{u_\tau^2}{u_\infty^2} \frac{\tau^r_{ij}}{u_\tau^2} = \frac{u_\tau^2}{u_\infty^2} \textbf{NN} \left( \frac{\bar \varDelta^2|\bar \alpha| \bar \alpha_{ij}}{u_\tau^2} \right) = \textbf{NN} \left( \frac{u_\tau^2}{u_\infty^2} \frac{\bar \varDelta^2|\bar \alpha| \bar \alpha_{ij}}{u_\tau^2} \right) = \textbf{NN} \left( \frac{\bar \varDelta^2|\bar \alpha| \bar \alpha_{ij}}{u_\infty^2} \right). \end{align}

During the training process, the weight parameters $\boldsymbol {W}^{(m)(m+1)}$ are optimized to minimize the mean square error (loss) given by

(2.25)\begin{equation} L=\frac{1}{3}\frac{1}{N_{b}}\sum_{l=1}^{3} \sum_{n=1}^{N_b}(s_{l,n}^{{fDNS}}-s_{l,n})^2, \end{equation}

where $N_{b}$ denotes the minibatch size of 256, and $s_{l,n}^{{fDNS}}$ corresponds to the anisotropic components $\tau ^r_{ij}$ obtained from fDNS. Here, the loss of each minibatch is calculated and then the weights are updated for every minibatch.

We apply the Adam algorithm (a type of gradient descent, Kingma & Ba Reference Kingma and Ba2014) and learning rate annealing method (Simonyan & Zisserman Reference Simonyan and Zisserman2014; He et al. Reference He, Zhang, Ren and Sun2016) to optimize the weights and enhance the training speed, respectively. Early stopping (Goodfellow et al. Reference Goodfellow, Bengio and Courville2016) is employed in the training process to avoid overfitting. The learning rate is initially 0.025, and is subsequently reduced by a factor of 10 if no improvement in the training error is observed over five epochs. If the learning rate is decreased three times but there is no reduction in the training error over five epochs, training is stopped. The training requires 150 to 300 epochs. The entire process to train and execute the present NN is carried out with the PyTorch open-source library in the Python programming environment.

We develop an NN-based SGS model whose input and output are $\bar \varDelta ^2|\bar \alpha |\bar \alpha _{ij}$ and $\tau ^r_{ij}$, respectively (NN-based velocity gradient model, called NNVGM hereafter). Note that the dimensions of input and output are matched to be the same. The velocity gradient tensor as an input has been used by many previous studies (Gamahara & Hattori Reference Gamahara and Hattori2017; Wang et al. Reference Wang, Luo, Li, Tan and Fan2018; Xie et al. Reference Xie, Wang, Li and Ma2019b; Zhou et al. Reference Zhou, He, Wang and Jin2019; Pawar et al. Reference Pawar, San, Rasheed and Vedula2020; Prat, Sautory & Navarro-Martinez Reference Prat, Sautory and Navarro-Martinez2020; Xie et al. Reference Xie, Wang, Li, Wan and Chen2020a,Reference Xie, Wang, Li, Wan and Chenb,Reference Xie, Wang and Weinanc; Park & Choi Reference Park and Choi2021; Kim et al. Reference Kim, Kim, Kim and Lee2022; Abekawa et al. Reference Abekawa, Minamoto, Osawa, Shimamoto and Tanahashi2023; Kim et al. Reference Kim, Park and Choi2024), and the same form of input was considered by Jamaat & Hattori (Reference Jamaat and Hattori2022) for one-dimensional Burgers turbulence. The NNVGM obtains inputs from each grid point and computes corresponding SGS stresses.

2.3. Forced homogeneous isotropic turbulence (FHIT) and filtering method

FHIT involves an additional forcing term in the Navier–Stokes equations at low wavenumbers (Ghosal et al. Reference Ghosal, Lund, Moin and Akselvoll1995; Rosales & Meneveau Reference Rosales and Meneveau2005; Park et al. Reference Park, Lee, Lee and Choi2006):

(2.26)\begin{equation} f_{i} = \epsilon_{t}~ \frac{u_{i}}{\sum_{0<\vert \boldsymbol{k} \vert <2} |\hat{u}_{j}(\boldsymbol{k})|^{2}}, \end{equation}

where $f_{i}$ is the forcing term, $\hat u_i$ is the Fourier coefficient of $u_i$, $\epsilon _{t}$ is the prescribed mean total dissipation rate determining the turbulent energy injection rate, $\vert \mathbf {k} \vert =\sqrt {k_x^2 + k_y^2 + k_z^2}$, and $k_i$ is the wavenumber in the $i$ direction. Note that $\epsilon _t$ encompasses both the resolved and modelled dissipation.

In FHIT, two distinct Reynolds numbers, $Re_{\lambda }=u_{{rms}}{\lambda }/\nu$ and $Re_L = UL/\nu$, can be defined. Here, $u_{{rms}}$ is the r.m.s. velocity fluctuations, $\lambda (= \sqrt {15 u_{{rms}}^2 \nu / \epsilon _t})$ is the Taylor microscale, $L$ is associated with the computational domain size of $2{\rm \pi} L \times 2{\rm \pi} L \times 2{\rm \pi} L$, and $U$ is the characteristics velocity such that $U=(\epsilon _t L)^{1/3}$. The Taylor microscale Reynolds number $Re_{\lambda }$ has been widely used in turbulence studies. Nevertheless, the extraction of $u_{rms}$ requires intricate experimental data or DNS. In contrast, $Re_L$ does not require such preliminary outcomes, and can be deduced from $\epsilon _t L/U^{3} = 1$, $\eta k_{{max,DNS}}$ and $N_{{DNS}}$, where $\eta$ is the Kolmogorov length scale, and $k_{{max,DNS}}$ and $N_{{DNS}}$ are the largest wavenumber and number of grid points in each direction in DNS, respectively. According to Pope (Reference Pope2000), $\eta k_{{max,DNS}}=3/2$ or $\varDelta _{{DNS}}/\eta = 2{\rm \pi} / 3$ is set to achieve sufficiently high resolution in DNS, where $\varDelta _{{DNS}}$ is the grid size in DNS. Then, $L = \varDelta _{{DNS}} N_{{DNS}} / (2{\rm \pi} ) = \eta N_{{DNS}}/3$. With these relations together with $\eta = (\nu ^3 / \epsilon _t )^{1/4}$, one can obtain $Re_L = (N_{{DNS}} / 3)^{4/3}$. Table 1 lists the cases of DNS at various Reynolds numbers by changing the number of grid points for FHIT, where each case is named as DNS$N_{{DNS}}$. Note that the cases listed in table 1 do not necessarily require actual DNS except for the case of DNS128 (the recursive procedure introduced below requires DNS only at a low Reynolds number), and LESs are conducted for higher Reynolds number cases through the recursive procedure.

Table 1. Cases of DNS at various Reynolds numbers for FHIT. Here, $Re_L = (N_{DNS}/3)^{4/3}$ and $\varDelta _{DNS}/\eta = 2{\rm \pi} / 3\ (N_{DNS} \varDelta _{DNS} = 2{\rm \pi} L)$. DNSs are performed for DNS128 and DNS256, and the corresponding $Re_\lambda$ values are given in this table. Here, $\varDelta _{{DNS}}$ and $N_{DNS}$ are the size and number of grid points in each direction in DNS, respectively.

Table 2 shows the cases considered for fDNS and LES of FHIT. The fDNS data are obtained by applying a transfer function $\hat {G}(\mathbf {k})$ to the Fourier-transformed DNS data $\hat u_i(\mathbf {k})$:

(2.27)$$\begin{gather} \hat{\bar u}_{i}(\boldsymbol{k}) = \hat{u}_{i}(\boldsymbol{k}) \hat{\bar{G}}(\boldsymbol{k}), \end{gather}$$

(2.28)$$\begin{gather}\hat{\bar{G}}(\boldsymbol{k}) = \int \exp({\rm i} \boldsymbol{k} {\cdot} \boldsymbol{r}) \bar{G} (\boldsymbol{r})\,{\rm d}\boldsymbol{r}. \end{gather}$$

The fDNS data can be obtained by applying a spectral cutoff filter or Gaussian filter in the spectral space. However, these filterings have limitations when the filtered variables are compared with LES results. The spectral cutoff filtering removes the Fourier coefficients $\hat u_i$ at $k > k_c$, and allows filtered variables to be placed on coarser grids of LES, where $k_c (={\rm \pi} / \bar \varDelta )$ is the cutoff wavenumber. However, this filtering fails to satisfy the realizability condition due to negative weights in the physical space (Vreman, Geurts & Kuerten Reference Vreman, Geurts and Kuerten1994). Additionally, the spectral cutoff filtering may exhibit non-local oscillations (Meneveau & Katz Reference Meneveau and Katz2000; Pope Reference Pope2000). On the other hand, the Gaussian filtering is free from these drawbacks, and requires filtering even when the wavenumber exceeds the LES limit (i.e. $k > {\rm \pi}/\bar \varDelta _{LES}$). Although the filter size may coincide with the LES grid size, the Gaussian filtered data are allocated at DNS grids, not at LES grids.

Table 2. Cases of fDNS and LES at various Reynolds numbers for FHIT. Here, $N_c$ is the number of grid points in each direction for fDNS or LES, and $\bar \varDelta$ denotes the filter size for fDNS ($\bar \varDelta _{fDNS}$) or the grid size for LES ($\bar \varDelta _{LES}$). The name of each case indicates fDNS$N_c / N_{DNS}$ or LES$N_c / N_{DNS}$. Note that $\varDelta _{DNS}/\eta = 2{\rm \pi} / 3$.

In the present study, we use the cut-Gaussian filtering (Zanna & Bolton Reference Zanna and Bolton2020; Guan et al. Reference Guan, Chattopadhyay, Subel and Hassanzadeh2022, Reference Guan, Subel, Chattopadhyay and Hassanzadeh2023; Pawar et al. Reference Pawar, San, Rasheed and Vedula2023) which retains the transfer function of Gaussian filtering but truncates the filtered variables at the wavenumbers exceeding the cutoff wavenumber $k_{c}$. Table 3 shows three different filters in the spectral space, corresponding transfer functions and largest wavenumbers of filtered variables, respectively. Figure 2 shows the energy spectra from DNS128, fDNS32/128 using cut-Gaussian, Gaussian and spectral cutoff filters, and LES32/128 with DSM and GM, respectively. The energy spectrum with GM is in excellent agreement with that of fDNS using the spectral cutoff filter. In contrast, the energy spectrum with DSM is closer to that of fDNS using the cut-Gaussian filter. Given the increased adaptability of DSM over GM, we measure the prediction capability of the SGS models based on the cut-Gaussian filtered DNS data.

Table 3. Comparison of three filters. Here, $H$ is the Heaviside step function.

Figure 2. Three-dimensional energy spectra: black line, DNS128 ($Re_L =149.09; Re_\lambda = 93.03$); red line, fDNS32/128 with cut-Gaussian filtering; blue line, fDNS32/128 with Gaussian filtering; green line, fDNS32/128 with spectral cutoff filtering; blue $\circ$, LES32/128 with DSM; green $\circ$, LES32/128 with GM.

3.1. A recursive algorithm

The present study employs a recursive procedure to construct an NNVGM adapted to various grid sizes normalized by the Kolmogorov length scale from the following steps:

1. (i) obtain training data (fDNS data) by filtering DNS data at a low Reynolds number ($Re_L$);

2. (ii) train an NNVGM with the fDNS data;

3. (iii) apply the NNVGM to LES at a higher Reynolds number, where the ratio of LES grid size to the Kolmogorov length scale ($\bar \varDelta _{LES}/\eta$) is equal to that of filter size to the Kolmogorov length scale ($\bar \varDelta _{fDNS}/\eta$) (see Appendix B for the effect of different filter to grid ratios);

4. (iv) filter the LES data and include the filtered LES (fLES) data in the training dataset, where the new filter size is twice the LES grid size;

5. (v) train the NNVGM using the augmented training data;

6. (vi) apply the updated NNVGM to LES at even higher Reynolds number, where the LES grid size is equal to the filter size defined in Step (iv);

7. (vii) repeat steps (iv)–(vi) for LES at higher Reynolds numbers.

Steps (i) and (ii) are similar to those used in constructing previous NN-based SGS models. The present recursive procedure aims at improving the performance of an NNVGM at a wider range of non-dimensional grid sizes by recursively performing LES and training it using fLES data.

3.2. NNVGM trained with fDNS data: a priori and a posteriori tests

First, fDNS data are used to train an NN. The Reynolds number and number of grid points in each direction for DNS are $Re_{L} = 149.09$ and $N_{DNS}=128$, respectively (table 1). The filter sizes considered are $\bar \varDelta _{fDNS}/\eta = 8.39$ and $4.19$, corresponding to the cases of fDNS32/128 and fDNS64/128, respectively (table 2). Including two different filter sizes in constructing fDNS data helps to improve the prediction capability of NNVGMs (see, for example, Park & Choi Reference Park and Choi2021) by broadening the ranges of the input and output, when they are appropriately normalized. In the present study, we train a dual NN (NN$_{N}$ and NN$_{S}$ in figure 1) with the training data (input and output) normalized by r.m.s. SGS normal ($\tau _{11, rms}^r$) and shear ($\tau _{12, rms}^r$) stress fluctuations, respectively. Note that $\tau _{11, rms}^r (=\tau _{22, rms}^r = \tau _{33, rms}^r)$ and $\tau _{12, rms}^r (=\tau _{13, rms}^r = \tau _{23, rms}^r)$ are obtained from (2.5) and (2.6) with DNS data. During actual LES, however, the SGS normal and shear stress fluctuations are a priori unknown, and thus providing normalized input data to the NN is not possible. Hence, in actual LES, we can only provide the input normalized by the characteristic velocity scale ($\bar {\varDelta }^2|\bar {\alpha }|\bar {\alpha }_{ij}/U_{LES}^2$). Thanks to the important property of the present NN, (2.22), we obtain the following relation (no summation on $i$ and $j$):

(3.1)\begin{equation} \frac{\tau^r_{ij}}{U_{LES}^2} = \frac{\tau_{ij, rms}^r}{U^2_{LES}} \frac{\tau^r_{ij}}{\tau_{ij, rms}^r} = \frac{\tau_{ij, rms}^r}{U^2_{LES}} \textbf{NN} \left( \frac{\bar \varDelta^2|\bar \alpha| \bar \alpha_{ij}}{\tau_{ij, rms}^r} \right) = \textbf{NN} \left( \frac{\bar \varDelta^2|\bar \alpha| \bar \alpha_{ij}}{U^2_{LES}} \right). \end{equation}

This relation indicates that, during actual LES, one can provide the input normalized by $U_{LES}$ to the NN trained with the input normalized by the r.m.s. SGS stresses.

However, the SGS stresses obtained from DNS data have a problem of imbalanced distribution because they tend to be mainly distributed around zero. So, we choose the normalized output and corresponding input through the process known as undersampling (Drummond & Holte Reference Drummond and Holte2003; Liu, Wu & Zhou Reference Liu, Wu and Zhou2008). During this process, we do not use all the filtered data as the training one, but choose data such that the possibility of choosing data as the training one decreases when its magnitude is near to zero: i.e. the probability ($\mathcal {P}$) to choose an SGS shear stress and corresponding inputs as training data is

(3.2)$$\begin{gather} \mathcal{P} = \begin{cases} \sin^{2}{\theta} & \text{if } 0 \le \theta < {\rm \pi}/{2} \\ 1 & \text{if } \theta \geq {\rm \pi}/{2} \end{cases}, \end{gather}$$

(3.3)$$\begin{gather}\theta = \frac{\rm \pi}{8} \frac{\sqrt{\{(\tau_{12}^{r})^{2}+(\tau_{23}^{r})^{2}+(\tau_{13}^{r})^{2}\}/3}}{\tau_{12,rms}^{r}}. \end{gather}$$

Figure 3 illustrates the effect of undersampling on the probability density function (p.d.f.) of training data. This undersampling reduces the probability of zero SGS stresses and enhances the occurrence of strong SGS stresses, leading to generate high SGS dissipation and thus improving the performance of an NNVGM. The SGS normal stresses are also similarly processed. For each fDNS dataset, approximately 600 000 pairs of $\bar \varDelta ^2 |\bar \alpha | \bar \alpha _{ij}$ and $\tau _{ij}^r$ are used for training. Other types of the probability $\mathcal {P}$ for undersampling are also considered and the results from different choices of $\mathcal {P}$ are discussed in Appendix C.

Figure 3. Probability density function (p.d.f.) of normalized $\tau _{12}^r$ before (blue) and after (red) undersampling.

We perform a priori tests with the present NNVGMs and traditional models using fDNS32/128 ($Re_L = 149.09$) and fDNS64/256 ($Re_L = 375.69$), respectively. Table 4 provides the correlations of the SGS normal and shear stresses and SGS dissipation, and the magnitude of the SGS dissipation from various SGS models. Here, NNVGM(32+64)/128 is trained using fDNS32/128 and fDNS64/128 datasets, and NNVGM(64+128)/256 is done using fDNS64/256 and fDNS128/256 datasets. For both Reynolds numbers, NNVGMs provide the highest correlations of the SGS stresses and SGS dissipation, and the magnitudes of the SGS dissipation closest to those of fDNS data. It is notable that NNVGM(32+64)/128 trained at $Re_L = 149.09$ successfully predicts the correlations and SGS dissipation even when it is applied to a higher Reynolds number of $Re_L = 375.69$.

Table 4. Statistics from a priori tests at $Re_L=149.09$ and 375.69 with $\bar \varDelta _{fDNS}/\eta =8.38$: Pearson correlation coefficients $(R)$ of $\tau _{ij}^{r}$ and SGS dissipation $\epsilon _{SGS} (=-\tau _{ij}^{r}\bar S_{ij})$, and magnitude of $\epsilon _{SGS}$.

Figure 4(a,b) shows the p.d.f.s of $\epsilon _{SGS}$, $\tau _{11}^{r}$ and $\tau _{12}^{r}$ from various SGS models for $Re_L = 149.09$ and 375.69, respectively. As shown, the p.d.f.s of $\epsilon _{SGS}$ from NNVGM and GM agree very well with fDNS data. The p.d.f.s from both CSM and DSM appear only at positive SGS dissipation, but do not agree well with fDNS data. However, DMM provides large positive SGS dissipation (its curves nearly fall on those of CSM), leading to overestimated SGS dissipation (see table 4). Similarly, NNVGM and GM predict SGS stresses better than other models. It should be noted that NNVGM(32+64)/128 trained at $Re_L=149.09$ predicts the p.d.f.s accurately at $Re_L=375.69$ as much as NNVGM(64+128)/256 does.

Figure 4. Probability density functions of the SGS dissipation and SGS normal and shear stresses (a priori test with $\bar \varDelta _{fDNS}/\eta =8.38$): (a) $Re_L=149.09$; (b) $Re_L=375.69$. In panel (a), $\blacksquare$, fDNS32/128; red line, NNVGM(32+64)/128; blue line, CSM32/128; blue dashed line, DSM32/128; blue dotted line, GM32/128; blue dash-dotted line, DMM32/128. In panel (b), $\blacksquare$, fDNS64/256; red $\circ$, NNVGM(64+128)/256; red line, NNVGM(32+64)/128; blue line, CSM64/256; blue dashed line, DSM64/256; blue dotted line, GM64/256; blue dash-dotted line, DMM64/256.

Now, we perform a posteriori tests (actual LES) with the present NNVGMs and traditional models for $Re_L=149.09$ and 375.69, and compare the results with fDNS32/128 and fDNS64/256 data, respectively. Figures 5 and 6 show the three-dimensional energy spectra, and p.d.f.s of the SGS dissipation, SGS stresses, strain rates ($\bar S_{11}$ and $\bar S_{12}$) and rotation rate ($\bar \varOmega _{23}=0.5 ( \partial \bar u_3 / \partial x_2 - \partial \bar u_2 / \partial x_3 )$), respectively. The energy spectra from NNVGM and GM agree well with fDNS data, whereas CSM, DSM and DMM provide rapid fall-off at high wavenumbers (figure 5). For p.d.f.s of the SGS dissipation and stresses, NNVGM, GM and DMM provide accurate predictions, whereas CSM and DSM do not. For the strain and rotation rates, NNVGM and GM provide the most accurate results (figure 6). These results indicate that LES with NNVGM trained at a Reynolds number predicts turbulence statistics very well for a different (higher) Reynolds number flow if $\bar \varDelta _{LES}/\eta = \bar \varDelta _{fDNS}/\eta$. A similar conclusion was also made by Park & Choi (Reference Park and Choi2021), where the grid size was non-dimensionalized in wall units rather than by the Kolmogorov length scale.

Figure 5. Three-dimensional energy spectra (a posteriori test; LES32/128 and LES64/256 for $Re_L=149.09$ and 375.69, respectively, with $\bar \varDelta _{LES}/\eta =8.38$): (a) $Re_L=149.09$; (b) $Re_L=375.69$. In panel (a), $\blacktriangle$, DNS128; $\blacksquare$, fDNS32/128; red line, NNVGM(32+64)/128; blue line, CSM32/128; blue dashed line, DSM32/128; blue dotted line, GM32/128; blue dash-dotted line, DMM32/128. In panel (b), $\blacktriangle$, DNS256; $\blacksquare$, fDNS64/256; red $\circ$, NNVGM(64+128)/256; red line, NNVGM(32+64)/128; blue line, CSM64/256; blue dashed line, DSM64/256; blue dotted line, GM64/256; blue dash-dotted line, DMM64/256. Here, $u_\eta$ is the Kolmogorov velocity scale.

Figure 6. Probability density functions of the SGS dissipation, SGS stresses, strain rates and rotation rate (a posteriori test; LES32/128 and LES64/256 for $Re_L=149.09$ and 375.69, respectively, with $\bar \varDelta _{LES}/\eta =8.38$): (a) $Re_L=149.09$; (b) $Re_L=375.69$. In panel (a), $\blacksquare$, fDNS32/128; red line, NNVGM(32+64)/128; blue line, CSM32/128; blue dashed line, DSM32/128; blue dotted line, GM32/128; blue dash-dotted line, DMM32/128. In panel (b), $\blacksquare$, fDNS64/256; red $\circ$, NNVGM(64+128)/256; red line, NNVGM(32+64)/128; blue line, CSM64/256; blue dashed line, DSM64/256; blue dotted line, GM64/256; blue dash-dotted line, DMM64/256.

3.3. NNVGM trained with fDNS and fLES data: a priori and a posteriori tests

As mentioned before, we train the NN with fLES data as well as fDNS data. The theoretical basis for filtering LES data can be attributed to the test filter introduced by Germano et al. (Reference Germano, Piomelli, Moin and Cabot1991):

(3.4)\begin{equation} \tilde{\phi}(\boldsymbol{x},t)=\int \phi(\boldsymbol{x-r},t)\tilde{G}(\boldsymbol{r}) \,{\rm d}\boldsymbol{r}, \end{equation}

where $\tilde G$ is a filter function and $\tilde {\bar G} = \tilde G \bar G$. The test-filtered SGS stresses are written as

(3.5)$$\begin{gather} T_{ij} = \widetilde{\overline{{u}_i{ u}_j}}-\tilde{\bar{u}}_i\tilde{\bar{u}}_j, \end{gather}$$

(3.6)$$\begin{gather}{T}^r_{ij} = {T}_{ij}-\tfrac{1}{3}{T}_{kk}\delta_{ij}, \end{gather}$$

and ${T}^r_{ij}$ is obtained by the following relation:

(3.7)\begin{equation} {T}^r_{ij} = \tilde{\tau}^r_{ij}+{L}^r_{ij}, \end{equation}

where

(3.8)$$\begin{gather} L_{ij} = \widetilde{\bar{u}_{i} \bar{u}_{j}}-\tilde{\bar{u}}_i\tilde{\bar{u}}_j, \end{gather}$$

(3.9)$$\begin{gather}{L}^r_{ij} = {L}_{ij}-\tfrac{1}{3}{L}_{kk}\delta_{ij}. \end{gather}$$

The present NNVGM is updated through the accumulation of new datasets with the input of $\tilde {\varDelta }^2|\tilde {\bar \alpha }|\tilde {\bar \alpha }_{ij}$ and the output of ${T}^r_{ij}$. The test filter size is twice the LES grid size, i.e. $\tilde {\varDelta }=2\bar \varDelta _{LES}$. The filtered LES data, referred to as fLES32/256, are obtained by filtering LES64/256 data. The fLES data are selected using the same normalization and undersampling techniques described in § 3.2. The fLES32/512, fLES32/1024, $\ldots$, fLES32/32768 data can be created by filtering LES64/512, LES64/1024, $\ldots$, LES64/32768, respectively, during the recursive procedure. Table 5 summarizes the present NNVGMs considered and corresponding training data.

Table 5. NNVGMs and corresponding training data. Here, NNVGM128D and NNVGM256DD are the same as NNVGM(32+64)/128 and NNVGM(64+128)/256 discussed in § 3.2, respectively. A directory with the NNs (NNVGM128D, NNVGM(128D+256L), NNVMG(128D+512L), $\ldots$, NNVMG(128D+32768L)) and a customizable Jupyter notebook can be accessed at https://www.cambridge.org/S0022112024009923/JFM-Notebooks/files/Table_5/Table_5.ipynb. Here, the weights $W_{ij}^{(1)(2)}, W_{jk}^{(2)(3)}$ and $W_{kl}^{(3)(4)}$ in (2.21) are also accessible for NNVGM128D to NNVGM(128D+32768L). The SGS stresses $\tau _{ij}^r/U^2$ from NNVGMs are calculated with user's input values of $\bar \alpha _{ij} L/U$.

Let us conduct a priori and a posteriori tests for $Re_L=375.69$ with $\bar \varDelta /\eta =16.76$ (twice that of LES64/256 and LES32/128). During the training process, NNVGM(128D+256D) employs fDNS data only (fDNS32/128, fDNS64/128 and fDNS32/256), but NNVGM(128D+256L) uses a combination of fDNS and fLES data (fDNS32/128, fDNS64/128 and fLES32/256). Table 6 and figure 7 show the results of a priori tests. For the correlation coefficients, the performances of NNVGMs are similar to those of GM and DMM, and are much better than those of CSM and DSM, whereas the SGS dissipation is best predicted by NNVGMs (table 6). For p.d.f.s (figure 7), similar predictions as discussed in figure 4 are observed, and thus we do not repeat here. A notable observation is that the results of NNVGM(128D+256L) are nearly identical to those obtained from NNVGM(128D+256D), implying that fLES data can replace fDNS data in constructing NNVGMs. Currently, we reach this conclusion only from numerical simulation, but a more mathematical analysis may be required to confirm the usefulness of fLES data in constructing the NN.

Table 6. Statistics from a priori test at $Re_L=375.69$ with $\bar \varDelta /\eta =16.76$: Pearson correlation coefficients $(R)$ of $\tau _{ij}^{r}$ and $\epsilon _{SGS}$, and magnitude of $\epsilon _{SGS}$.

Figure 7. Probability density functions of the SGS dissipation and stresses (a priori test at $Re_L=375.69$ with $\bar \varDelta /\eta =16.76$): $\blacksquare$, fDNS32/256; red $\circ$, NNVGM(128D+256D); red line, NNVGM(128D+256L); blue line, CSM32/256; blue dashed line, DSM32/256; blue dotted line, GM32/256; blue dash-dotted line, DMM32/256.

Figures 8 and 9 show the results of energy spectra and p.d.f.s from a posteriori tests, where the result from GM32/256 is not shown because the simulation diverged. Again, NNVGM(128D+256D) and NNVGM(128D+256L) perform well among the SGS models considered, but show underpredictions at intermediate wavenumbers. The DMM predicts p.d.f.s very well, but its energy spectrum rapidly decreases at high wavenumbers. Thus, we perform additional LES with NNVGM(128D+1024L) whose trained filter sizes are $4.19 \le \bar \varDelta _{fDNS}/\eta$ and $\bar \varDelta _{fLES}/\eta \le 67.02$ (note that the filter sizes for NNVGM(128D+256L) are $4.19 \le \bar \varDelta _{fDNS}/\eta$ and $\bar \varDelta _{fLES}/\eta \le 16.76$). The energy spectrum with NNVGM(128D+1024L) is in excellent agreement with fDNS32/256, suggesting that one should expect successful LES when $\bar \varDelta _{LES} / \eta$ is within the range of trained filter sizes (see § 4.2 for further discussion). It is also important to note that the present NNVGMs do not show the inconsistency between a priori and a posteriori tests observed from the traditional SGS models (Park et al. Reference Park, Yoo and Choi2005).

Figure 8. Three-dimensional energy spectra (a posteriori test at $Re_L=375.69$ with $\bar \varDelta _{LES}/\eta =16.76$; $N^3 = 32^3$): $\blacktriangle$, DNS256; $\blacksquare$, fDNS32/256; red $\circ$, NNVGM(128D+256D); red line, NNVGM(128D+256L); red dotted line, NNVGM(128D+1024L); blue line, CSM32/256; blue dashed line, DSM32/256; blue dash-dotted line, DMM32/256.

Figure 9. Probability density functions (a posteriori test at $Re_L=375.69$ with $\bar \varDelta _{LES}/\eta =16.76$; $N^3 = 32^3$): (a) SGS dissipation and stresses; (b) strain and rotation rates. $\blacksquare$, fDNS32/256; red $\circ$, NNVGM(128D+256D); red line, NNVGM(128D+256L); red dotted line, NNVGM(128D+1024L); blue line, CSM32/256; blue dashed line, DSM32/256; blue dash-dotted line, DMM32/256.

3.4. A posteriori test using a finite difference method with a box filter

In this subsection, we apply a trained NNVGM (from a spectral method with the cut-Gaussian filter) to FHIT using the second-order central difference method (CD2) with a box filter, to see if the NNVGM trained from one numerical method can be successfully applied to LES using a different numerical method. All the spatial derivative terms in the filtered continuity and Navier–Stokes equations (2.3) and (2.4) are discretized using CD2 in staggered grids, and the filtered flow variable is obtained using a box filter in the physical space as

(3.10)$$\begin{gather} \bar \phi(\boldsymbol{x},t)= \int_{-\infty}^{\infty} \bar{G}(\boldsymbol{r}) \phi(\boldsymbol{x-r},t) \,{\rm d}\boldsymbol{r}, \end{gather}$$

(3.11)$$\begin{gather}\bar{G}(\boldsymbol{r})=\begin{cases} 1/\bar \varDelta^3 & \text{if } |\boldsymbol{r}| \le \bar \varDelta /2 \\ 0 & \text{if } |\boldsymbol{r}| > \bar \varDelta /2 \end{cases}. \end{gather}$$

We perform DNS of FHIT at $Re_L=375.69$ with $N^3 = 256^3$ and LESs with $N^3 = 32^3$ and $64^3$ using NNVGM(128D+32768L) and traditional SGS models, respectively. The results are given in figure 10. The energy spectra from DNS and fDNS using the CD2 and box filter are in excellent agreements with those using the spectral method and cut-Gaussian filter (figure 10a). The predictions of the energy spectrum by SGS models with $N^3 = 32^3$ are quite similar among themselves but show some deviations from the fDNS data (figure 10b). On the other hand, with $N^3 = 64^3$, the predictions are much better than those with $N^3 = 32^3$, and reasonably agree with the fDNS data (figure 10c). Hence, the prediction using the NNVGM trained from the spectral method is quite accurate, even though LES is performed with the CD2 and box filter, which may indicate the potential of the present recursive procedure to flow over/inside a complex geometry.

Figure 10. Three-dimensional energy spectra (a posteriori test at $Re_L=375.69$ using CD2 and box filter): (a) DNS256 (black line, spectral method; $\square$, CD2) and fDNS64/256 (dashed line, spectral method and cut-Gaussian filter; $\blacksquare$, CD2 and box filter); (b) a posteriori test with $N^3 = 32^3$; (c) a posteriori test with $N^3 = 64^3$. In panels (b) and (c), $\blacksquare$, fDNS; red line, NNVGM(128D+32765L); blue line, CSM; blue dashed line, DSM; blue dash-dotted line, DMM; dashed line, no SGS model (diverge with $N^3 = 32^3$).

4.1. Forced homogeneous isotropic turbulence at high Reynolds numbers

LESs of FHIT with the number of grid points of $64^3$ (LES64/256, LES64/512, LES64/1024, LES64/2048, LES64/4096, LES64/8192, LES64/16384 and LES64/32768) are conducted at eight different Reynolds numbers ($Re_{L} = 375.69\unicode{x2013} 242347.33$) using nine different SGS models (NNVGM128D, NNVGM(128D+512L), NNVGM(128D+2048L), NNVGM(128D+8192L), NNVGM(128D+32768L), CSM, DSM, GM and DMM). The computational time step is set at $\Delta t U/L=0.0025$, at which the Courant–Friedrichs–Lewy (CFL) numbers are below 0.3.

Figure 11 shows the three-dimensional energy spectra at $Re_L = 946.67$, $6010.98$, $38167.31$ and 242347.33 using NNVGM128D, NNVGM(128D+512L), NNVGM(128D+2048L), NNVGM(128D+8192L) and NNVGM(128D+32768L), respectively. The energy spectra predicted from NNVGMs trained with low-Reynolds-number data such as NNVGM128D and NNVGM(128D+512L) show energy pile-up at high wavenumbers for high-Reynolds-number simulations. This energy pile-up at high wavenumbers decreases significantly by NNVGMs trained with higher-Reynolds-number data. The results from NNVGM128D indicate an inherent limitation of the non-recursive NN-based SGS model, when the LES grid size is bigger than the filter sizes of fDNS employed for training the NN. It is interesting to note that the energy spectrum for the case of $Re_L = 242347.33$ can be reasonably predicted by NNVGM(128D+2048L) and NNVGM(128D+8192L) that includes the training data up to $Re_L=6010.98$ and 38167.31, respectively. This result indicates that, when the training data contain the length scales in the inertial range large enough compared with the dissipation scale, the present SGS model should predict isotropic turbulence even at the Reynolds number higher than that trained.

Figure 11. Three-dimensional energy spectra from various NNVGMs ($N^3 = 64^3$): (a) $Re_L = 946.67$; (b) 6010.98; (c) 38167.31; (d) 242347.33. $\blacksquare$, filtered Kolmogorov energy spectrum, $E(k)=\frac {3}{2}{\epsilon }^{2/3} {k}^{-5/3} \exp (-k^2 \bar \varDelta ^2/12)$; red dotted line, NNVGM128D; red dash-dot-dotted line, NNVGM(128D+512L); red dash-dotted line, NNVGM(128D+2048L); red dashed line, NNVGM(128D+8192L); red line, NNVGM(128D+32768L).

Figure 12 shows the three-dimensional energy spectra at eight Reynolds numbers from various SGS models including traditional ones. First, LESs with GM reveal instability at high Reynolds numbers (see also Vollant, Balarac & Corre Reference Vollant, Balarac and Corre2016). LESs without SGS model exhibit non-physical energy pile up at high wavenumbers (figure 12f) due to insufficient dissipation, as observed from LESs with NNVGM128D and NNVGM(128D+512L) in figure 11. NNVGM(128D+32768L) suitably calculates the energy spectra that are quite similar to those of CSM, DSM and DMM. The predicted energy spectra from NNVGM, CSM, DSM and DMM follow the Kolmogorov energy spectrum at low and intermediate wavenumbers.

Figure 12. Three-dimensional energy spectra at $Re_L = 375.69\unicode{x2013}242347.33$ from various SGS models ($N^3 = 64^3$): (a) NNVGM(128D+32768L); (b) CSM; (c) DSM; (d) GM; (e) DMM; ( f) no SGS model. The black dashed line denotes the Kolmogorov energy spectrum, $E(k)=\frac {3}{2}{\epsilon }^{2/3} {k}^{-5/3}$.

In §§ 3.2 and 3.3, we observed that GM and NNVGM performed similarly in a priori tests for small filter sizes, whereas, in actual LES, GM diverged but NNVGM performed stably for large grid sizes. One of the reasons for these different results is that the two models compute the SGS dissipation differently. As shown in table 6, for a large filter size, GM predicts the SGS dissipation significantly smaller than fDNS data, while NNVGM computes it accurately. Another reason may be that the error of the SGS stress prediction by GM grows exponentially with increasing grid size (Vreman et al. Reference Vreman, Geurts and Kuerten1996), since GM is derived from Taylor expansion. However, NNVGM is trained with the SGS stresses from small to large filter sizes, and thus the error of the SGS stress prediction by NNVGM does not increase with increasing grid size.

4.2. Decaying homogeneous isotropic turbulence

The present approach is based on a hypothesis that the NN-based SGS models in LES successfully predict the turbulence statistics when the LES grid size normalized by the Kolmogorov length scale is within the range of filter sizes used for training (see also Park & Choi Reference Park and Choi2021). In figure 13, we indicate the filter sizes (normalized by the Kolmogorov length scale) used for fDNSs and fLESs of FHIT with NNVGMs and the normalized grid sizes used for LESs of two DHITs. This figure should be understood as follows: for DHIT of Comte-Bellot & Corrsin (Reference Comte-Bellot and Corrsin1971), NNVGM(128D+1024L) or NNVGM(128D+32768L) can be used for successful LES, but not NNVGM128D and NNVGM(128D+256L); for DHIT of Kang, Chester & Meneveau (Reference Kang, Chester and Meneveau2003), only NNVGM(128D+32768L) can be used for successful LES. To examine this hypothesis, LESs are performed for two DHITs (untrained flows) of Comte-Bellot & Corrsin (Reference Comte-Bellot and Corrsin1971) and Kang et al. (Reference Kang, Chester and Meneveau2003) with the present NNVGMs trained only with FHIT data. This task should be a notable challenge due to the transient nature of DHIT as opposed to FHIT.

Figure 13. Illustration of the filter-size ranges of fDNS and fLES used to generate each NNVGM (FHIT), and the grid-size ranges required for the simulations of DHITs (Comte-Bellot & Corrsin Reference Comte-Bellot and Corrsin1971; Kang et al. Reference Kang, Chester and Meneveau2003). On the top of this figure (FHIT), each solid circle ($\bullet$) denotes the filter size of fDNS or fLES used for generating training data. For each NNVGM, the line connecting the first and last solid circles indicates the range of grid sizes for successful LES. On the bottom of this figure (DHIT), each solid square ($\blacksquare$) indicates the grid size normalized by the initial Kolmogorov length scale. As time goes by, the Kolmogorov length scale increases and thus $\bar {\varDelta }_{LES}/\eta$ decreases denoted as a thick dashed line.

The first DHIT to simulate is the experiment by Comte-Bellot & Corrsin (Reference Comte-Bellot and Corrsin1971) (called CBC hereafter). In CBC, turbulence data were experimentally generated through grid turbulence using a mesh size of $M=5.08$ cm and free-stream velocity of $U_0=10$ m s$^{-1}$. The Reynolds number based on the Taylor microscale was $Re_{\lambda }=71.6$ at $tU_0/M=42$. For LES, flow variables are non-dimensionalized by $L=11M /(2{\rm \pi} )$ and $U=\sqrt {3/2} u_{rms} \vert _{tU_0/M=42}$ (Lee, Choi & Park Reference Lee, Choi and Park2010). The numbers of grid points tested are $N^3=32^3$ and $64^3$, respectively. A divergence-free initial field is obtained using the logistic polynomial approximation method (Knight et al. Reference Knight, Zhou, Okong'o and Shukla1998) and rescaling method (Kang et al. Reference Kang, Chester and Meneveau2003). LESs are performed 50 times by varying initial fields to obtain ensemble averages. The computational time step is fixed to be ${\Delta } t U / L = 0.01$. As turbulence decays in time, the Kolmogorov length scale $\eta$ increases in time, and thus $\bar \varDelta _{LES}/\eta = 60.2$ and 26.5 ($N^3 = 32^3$), and $30.1$ and $13.2$ ($N^3 = 64^3$) at $t U_0 / M=42$ and $171$, respectively. Each LES starts from large $\bar {\varDelta }_{LES}/\eta$ (denoted as a solid square in figure 13) and moves to smaller $\bar {\varDelta }_{LES}/\eta$ due to increased $\eta$ in time (denoted as a thick dashed line).

The LES results on CBC DHIT with $N^3 = 32^{3}$ (starting from $\bar \varDelta _{LES} / \eta = 60.2$) are given in figure 14. Without SGS model, the turbulent kinetic energy decays very slowly and energy pile up occurs at high wavenumbers. With GM, resolved kinetic energy first decreases and increases later in time, and simulation finally diverges, as shown by Vreman et al. (Reference Vreman, Geurts and Kuerten1996). DSM and CSM are relatively good at predicting turbulence decay and energy spectra, but they overestimate the energy spectra at intermediate wavenumbers within the inertial range. However, NNVGM(128D+32768L) performs best in predicting energy spectra, while NNVGM(128D+1024L) and DMM also perform better than other types of SGS models (DMM shows energy fall-off at high wavenumbers), whereas NNVGM(128D+1024L) shows energy pile up at high wavenumbers, possibly because $\bar \varDelta _{LES}/ \eta \ (\le 60.2)$ is only slightly smaller than $\bar \varDelta _{fLES} / \eta =67.02$ (fLES32/1024). Note also that the prediction with NNVGM128D ($\bar \varDelta _{LES} / \eta \gg \bar \varDelta _{fDNS} / \eta =8.38$) is not good at all and worse than CSM and DSM, validating our conjecture on the choice of NNVGM for successful LES (figure 13). With $N^3 = 64^3$ (starting from $\bar \varDelta _{LES} / \eta = 30.1$), NNVGM(128D+1024L) also performs very good as much as NNVGM(128D+32768L) does (figure 15), because $\bar \varDelta _{LES} / \eta \ (\le 30.1)$ is sufficiently smaller than $\bar \varDelta _{fLES} / \eta =67.02$, while other SGS models (except DMM) still show some deviations from experimental data.

Figure 14. LES of CBC DHIT with $N^3 = 32^3$: (a) resolved turbulent kinetic energy; (b) $E(k)$ at $tU_{0}/M=42$, $98$ and $171$. $\blacksquare$, filtered CBC data; red dotted line, NNVGM128D; red dashed line, NNVGM(128D+1024L); red line, NNVGM(128D+32768L); blue line, CSM; blue dashed line, DSM; blue dotted line, GM; blue dash-dotted line, DMM; dashed line, no SGS model.

Figure 15. LES of CBC DHIT with $N^3 = 64^3$: (a) resolved turbulent kinetic energy; (b) $E(k)$ at $tU_{0}/M=42$, $98$ and $171$. $\blacksquare$, filtered CBC data; red dotted line, NNVGM128D; red dashed line, NNVGM(128D+1024L); red line, NNVGM(128D+32768L); blue line, CSM; blue dashed line, DSM; blue dotted line, GM; blue dash-dotted line, DMM; dashed line, no SGS model.

The second DHIT to simulate is the experimental one by Kang et al. (Reference Kang, Chester and Meneveau2003). The Reynolds number is much higher ($Re_{\lambda }=716$ at $tU_{0}/M=20$) than that of CBC DHIT, where $U_{0}=11.2$ m s$^{-1}$ and $M=0.152$ m. We consider the present NNVGMs and traditional SGS models. LESs are performed 50 times by varying initial fields to obtain ensemble averages. The flow variables are non-dimensionalized by $L = 33.68M/(2{\rm \pi} )$ and $U=\sqrt {3/2} u_{rms} \vert _{tU_0/M=20}$. The numbers of grid points tested are $N^3=32^3$ and $64^3$. The flows are initiated using the energy spectrum and rescaling method as done by Kang et al. (Reference Kang, Chester and Meneveau2003). For $N^3 = 32^3$, $\bar \varDelta _{LES}/\eta =1454.5$ and 888.9 at $tU_{0}/M=20$ and $48$, respectively. The trained grid sizes of NNVGM(128D+32768L) ($4.19 \le \bar \varDelta _{fDNS} / \eta$ and $\bar \varDelta _{fLES} / \eta \le 2144.66$) include this range of $\bar \varDelta _{LES}/\eta$, whereas those of NNVGM128D ($4.19 \le \bar \varDelta _{fDNS} / \eta \le 8.38$) do not (see figure 13), suggesting that NNVGM(128D+32768L) should properly predict turbulence decay and variation of the energy spectra of DHIT by Kang et al. (Reference Kang, Chester and Meneveau2003).

The LES results on DHIT by Kang et al. (Reference Kang, Chester and Meneveau2003) with $N^3 = 32^3$ are shown in figure 16. While CSM and DSM overpredict the energy spectra even at intermediate wavenumbers and thus overestimate turbulent kinetic energy, DMM predicts the resolved kinetic energy and energy spectra most accurately. NNVGM(128D+32768L) slightly overestimates the resolved turbulent kinetic energy due to its slight overprediction of energy at high wavenumbers. This may be because $\bar \varDelta _{fLES} / \eta =2144.66$ (fLES32/32768) is only slightly larger than $\bar \varDelta _{LES}/ \eta ~(\le 1454.5)$. By increasing the number of grid points ($N^3=64^3$), the predictions by NNVGM(128D+32768L) become much better (figure 17). It is interesting to see that the level of kinetic energy without SGS model is almost constant in time due to the energy pile up at high wavenumbers. This is because the LES grid size ($\bar \varDelta _{LES}/\eta \approx 10^3$) is so large that the corresponding eddies do not properly dissipate energy at these length scales, and thus total kinetic energy is nearly conserved without SGS model.

Figure 16. LES of DHIT by Kang et al. (Reference Kang, Chester and Meneveau2003) with $N^3 = 32^3$: (a) resolved turbulent kinetic energy; (b) $E(k)$ at $tU_{0}/M= 20$, 30, 40 and $48$. $\blacksquare$, filtered experimental data; red dotted line, NNVGM128D; red line, NNVGM(128D+32768L); blue line, CSM; blue dashed line, DSM; blue dotted line, GM; blue dash-dotted line, DMM; dashed line, no SGS model.

Figure 17. LES of DHIT (Kang et al. Reference Kang, Chester and Meneveau2003) with $N^3 = 64^3$: (a) resolved turbulent kinetic energy; (b) $E(k)$ at $tU_{0}/M=20$, 30, 40 and $48$. $\blacksquare$, filtered experimental data; red dotted line, NNVGM128D; red line, NNVGM(128D+32768L); blue line, CSM; blue dashed line, DSM; blue dotted line, GM; blue dash-dotted line, DMM; dashed line, no SGS model.

Finally, the effects of the present NNVGMs on the prediction of the resolved kinetic energy and spectrum of DHIT of Kang et al. (Reference Kang, Chester and Meneveau2003) are shown in figure 18. As shown before, NNVGM128D does not provide an accurate energy spectrum, because it is trained with the data from $4.19 \le \bar \varDelta _{fDNS}/\eta \le 8.38$ and does not have sufficient information on $8.38 < \bar \varDelta _{LES} / \eta \le 727.3$ (see figure 13 for the grid size range suggested for LES of DHIT of Kang et al. Reference Kang, Chester and Meneveau2003). As the NNVGM is recursively updated (from NNVGM128D to NNVGM(128D+32768L)), its prediction becomes increasingly accurate. NNVGM(128D+2048L) provides quite an accurate energy spectrum albeit showing energy pile up at high wavenumbers. This is because NNVGM(128D+2048L) is trained with the data from $4.19 \le \bar \varDelta _{fDNS}/\eta$ and $\bar \varDelta _{fLES}/\eta \le 134.04$ and contains large length scales in the inertial range. As expected, NNVGM(128D+32768L) provides accurate energy spectra at all wavenumbers.

Figure 18. Effects of NNVGMs on the prediction of DHIT (Kang et al. Reference Kang, Chester and Meneveau2003) with $N^3 = 64^3$: (a) resolved turbulent kinetic energy; (b) $E(k)$ at $tU_{0}/M= 48$. $\blacksquare$, filtered experimental data; red dotted line, NNVGM128D; red dash-dot-dotted line, NNVGM(128D+512L); red dash-dotted line, NNVGM(128D+2048L); red dashed line, NNVGM(128D+8192L); red line, NNVGM(128D+32768L).

4.3. Computational cost

The amounts of CPU time required for estimating the SGS stresses and advancing one computational time step for the traditional and NN-based SGS models are given in table 7. With the same number of grids ($N^3=64^3$), the amounts of CPU time required to obtain the SGS stresses by DSM, DMM and NNVGM are approximately 3.3, 5.1 and 2.8 times that by CSM, respectively. Therefore, the computational cost of NNVGM is higher than CSM but lower than those of DSM and DMM, which is consistent with the report of Kang, Jeon & You (Reference Kang, Jeon and You2023) studying a similar HIT. Note, however, that the total amounts of CPU time for advancing one computational time step are relatively similar among themselves. The reason that the CPU time of DSM or DMM for determining the SGS stresses is higher than that of NNVGM is because DSM and DMM require Fourier and inverse Fourier transforms to determine the dynamic coefficient(s). However, for flow over a complex geometry, such transforms are not required based on a finite difference method. In that case, an NN-based SGS model may require more CPU time than DSM or DMM (see Kim et al. Reference Kim, Park and Choi2024 as an example).

Table 7. Amounts of CPU time (second) required for estimating the SGS stresses and advancing one computational time step, respectively. Simulations are performed using eight CPU cores (Intel(R) Core(TM) i9-11900KF CPU @ 3.50 GHz), and the amounts of CPU time are obtained by averaging over 50 computational time steps. In the present simulation, $\tau _{ij}$ was obtained using one GPU (GeForce RTX 3060) for NNVGM(128D+32768L), but, for comparison with traditional models, the amounts of CPU time using eight CPU cores are provided in this table.