A recursive neural-network-based subgrid-scale model for large eddy simulation: application to homogeneous isotropic turbulence

Abstract

We introduce a novel recursive procedure to a neural-network-based subgrid-scale (NN-based SGS) model for large eddy simulation (LES) of high-Reynolds-number turbulent flow. This process is designed to allow an SGS model to be applicable to a hierarchy of different grid sizes without requiring expensive filtered direct numerical simulation (DNS) data: (1) train an NN-based SGS model with filtered DNS data at a low Reynolds number; (2) apply the trained SGS model to LES at a higher Reynolds number; (3) update this SGS model with training data augmented with filtered LES (fLES) data, accommodating coarser filter size; (4) apply the updated NN to LES at a further higher Reynolds number; (5) go back to Step (3) until a target (very coarse) filter size divided by the Kolmogorov length scale is reached. We also construct an NN-based SGS model using a dual NN architecture whose outputs are the SGS normal stresses for one NN and the SGS shear stresses for the other NN. The input is composed of the velocity gradient tensor and grid size. Furthermore, for the application of an NN-based SGS model trained with one flow to another flow, we modify the NN by eliminating bias and introducing a leaky rectified linear unit function as an activation function. The present recursive SGS model is applied to forced homogeneous isotropic turbulence (FHIT) and successfully predicts FHIT at high Reynolds numbers. The present model trained from FHIT is also applied to decaying homogeneous isotropic turbulence and shows an excellent prediction performance.

1. Introduction

Large eddy simulation (LES) is an effective tool for predicting turbulent flow because it does not require significant computational resources towards resolution of smaller eddies. This resource reduction is realized through modelling of eddy motions at scales smaller than the grid size (called subgrid scales). The subgrid-scale (SGS) stress from these SGS motions should be modelled to close the governing equations for LES.

The SGS models have been derived for decades based on turbulence theory and statistical approximation. A subset of these models that favour simplicity and stability are represented as linear eddy-viscosity models. Examples include the Smagorinsky (Smagorinsky 1963) and dynamic Smagorinsky (Germano et al. 1991) models, based on the theory that the rate of energy transfer to smaller eddies is counterbalanced by the viscous dissipation within the inertial subrange, WALE model (Nicoud & Ducros 1999), using square of the velocity gradient tensor and analysing proper near-wall scaling for the eddy viscosity, and Vreman model (Vreman 2004), rooted in the principle that the SGS stress should be reduced in near-wall region or vanish in laminar flow. On the other hand, the scale-similarity model (SSM) (Bardina, Ferziger & Reynolds 1980) assumed that the interplay between resolved and modelled eddies can be accurately delineated by the difference between the filtered and doubly filtered velocities. The gradient model (GM) (Clark, Ferziger & Reynolds 1979) was derived from the Taylor expansion of the definition of SGS stress. Nevertheless, these models contain some limitations. The linear eddy-viscosity models, in particular, show very low correlations with true SGS stresses in a priori tests (Liu, Meneveau & Katz 1994; Salvetti & Banerjee 1995; Park, Yoo & Choi 2005), while they successfully predict turbulence statistics for various flows in a posteriori tests (i.e. actual simulations) (see, for example, Germano et al. (1991) for turbulent channel flow, Le Ribault, Sarkar & Stanley (1999) for turbulent jet, Kravchenko & Moin (2000) for flow over a circular cylinder). In contrast, both SSM and GM exhibit improved results in a priori tests, but they do not sufficiently dissipate turbulent kinetic energy in actual simulations, thereby leading to numerical instability (Vreman, Geurts & Kuerten 1996, 1997).

An alternative approach for the derivation of an SGS model involves a direct employment of direct numerical simulation (DNS) data. For example, the optimal LES (Langford & Moser 1999; Völker, Moser & Venugopal 2002) adopted stochastic estimation techniques (Adrian et al. 1989; Adrian 1990) to diminish the discrepancy between the ideal and simulated flow variables. This approach necessitated the use of multipoint correlation data as inputs, which in turn required corresponding DNS data. Moser et al. (2009) presented an optimal LES that did not require DNS data, but the SGS model was created based on the assumption of isotropic flow. Another example is an employment of deep learning techniques with a particular emphasis on artificial neural networks (NNs) (Sarghini, de Felice & Santini 2003; Gamahara & Hattori 2017; Maulik et al. 2018). Creating an SGS model through NNs requires the accumulation of extensive data through filtering DNS data. This approach is based on an assumption that there may be a complex but well-defined relationship between the SGS stress and resolved flow variables. The pioneering application of NNs in the calculation of SGS stress, to the best of our knowledge, was implemented by Sarghini et al. (2003). They employed an NN to find an optimal turbulent viscosity coefficient for a hybrid model (combined Smagorinsky and similarity models) applied to turbulent channel flow with nine velocity gradients and six resolved Reynolds stresses as inputs, and demonstrated its ability to accurately approximate the coefficient.

With the rapid development of deep learning technology, more advanced techniques have been employed to develop SGS models for diverse applications. NN-based SGS models may be classified into three categories: closed-form, reassembling and direct NN-based SGS models, respectively. The closed-form NN-based SGS models, analogous to the study conducted by Sarghini et al. (2003), postulate that SGS models may be derived in a closed form from physical, empirical or analytical tools. The coefficients of these models are adjustable, and are computed using NNs with resolved flow variables as inputs. Wollblad & Davidson (2008) identified the coefficients pertinent to proper orthogonal decomposition for turbulent channel flow through an NN. Beck, Flad & Munz (2019) and Pawar et al. (2020) suggested that calculating eddy viscosity could be employed as a surrogate to the direct computation of SGS stress and contribute to enhanced stability. Moreover, Xie et al. (2019a), Xie, Yuan & Wang (2020d) and Wang et al. (2021) suggested that the SGS stress can be evaluated through a polynomial function including the strain rate, rotation rate, velocity gradient and grid size, with NNs used to determine the coefficients of the polynomial. Yu et al. (2022) and Liu et al. (2023) obtained the coefficients of Smagorinsky and helicity SGS models using NN, respectively.

Reassembling NN-based SGS models extract unfiltered variables from filtered ones, and subsequently use these extracted variables for the computation of SGS stress. In this approach, NNs play a prominent role in the unfiltering or filtering process. Maulik et al. (2018) engaged NNs in the training for unfiltering (Maulik & San 2017) and filtering for two-dimensional decaying homogeneous isotropic turbulence (DHIT) in LES. The first NN collected filtered vorticity and stream function data across multiple grids and gauged the values of unfiltered vorticity and stream function at a single grid. The second NN was subsequently deployed to approximate the filtered nonlinear term, a crucial element for the determination of the SGS stress. Aligning with this methodology, Yuan, Xie & Wang (2020) employed an NN for the estimation of SGS stresses in three-dimensional forced homogeneous isotropic turbulence (FHIT), specifically using the NN to unfilter the velocities.

Direct NN-based SGS models derive the SGS stress or force directly from resolved flow variables through the application of NNs. Gamahara & Hattori (2017) employed NNs to compute the SGS stress as immediate outputs for turbulent channel flow from four input groups such as the strain rate, rotation rate, wall distance and velocity gradient. Previous direct NN-based SGS models have used various techniques and methods to enhance model performance and stability. These include an ad hoc method such as wall-damping function or clipping (Gamahara & Hattori 2017; Maulik et al. 2019; Zhou et al. 2019), derivation of inputs from multiple grids for single or multiple outputs (Beck et al. 2019; Maulik et al. 2019; Xie et al. 2019b; Zhou et al. 2019; Pawar et al. 2020; Sirignano, MacArt & Freund 2020; Xie et al. 2020a,b; Xie, Wang & Weinan 2020c; MacArt, Sirignano & Freund 2021; Stoffer et al. 2021; Cheng et al. 2022; Guan et al. 2022; Liu et al. 2022; Guan et al. 2023), incorporation of second derivatives of velocities (Wang et al. 2018; Xie et al. 2019b; Pawar et al. 2020; Sirignano et al. 2020; Xie et al. 2020a; MacArt et al. 2021) and consideration of wall distance or filter size as inputs (Gamahara & Hattori 2017; Zhou et al. 2019; Abekawa et al. 2023). Others (Sirignano et al. 2020; MacArt et al. 2021; Guan et al. 2023; Sirignano & MacArt 2023) defined the training error (or objective function) using flow variables other than the SGS stresses. Sirignano et al. (2020), MacArt et al. (2021) and Sirignano & MacArt (2023) trained NNs with adjoint-based, PDE (partial differential equation)-constrained optimization methods. Other improvements involve retraining NNs through transfer learning for higher Reynolds number flows (Subel et al. 2021; Guan et al. 2022) and training NNs in complex flows (MacArt et al. 2021; Sirignano & MacArt 2023; Kim, Park & Choi 2024). On the other hand, other types of techniques to improve the performance of NN-based SGS models have also been explored, such as the use of a convolutional neural network (Beck et al. 2019; Pawar et al. 2020; Guan et al. 2022; Liu et al. 2022; Guan et al. 2023), reinforcement learning (Novati, de Laroussilhe & Koumoutsakos 2021; Kim et al. 2022; Kurz, Offenhäuser & Beck 2023) and a graph neural network (Abekawa et al. 2023).

Park & Choi (2021) developed a direct NN-based SGS model for turbulent channel flow whose input and output were the velocity gradient or strain rate and the SGS stress, respectively. They performed a comprehensive analysis of the salient characteristics inherent to direct NN-based SGS models, and demonstrated from a posteriori test without introducing ad hoc techniques such as wall-damping function or clipping that SGS models relying on a single grid input can compute turbulence statistics with a reasonably decent level of accuracy. In addition, two specific results from this study are noteworthy: first, when the LES grid size falls within the range of filter sizes for training NN, the SGS model successfully predicts the flow; second, the SGS model is also successful in predicting the flow at higher Reynolds number when the grid size in wall units is the same as that trained. Hence, these observations indicate that the applicability of direct NN-based SGS models depends on the dimensionless LES grid size relative to the trained grid (or filter) size. For the development of a general NN-based SGS model for complex flow, it is essential that it should function properly over a wide range of non-dimensional filter or grid sizes.

Despite their promising characteristics, previous direct NN-based SGS models have confronted substantial challenges when aimed for application to untrained flows. These challenges originate from the difficulty in developing universal non-dimensional input and output variables for various flows and the nonlinear property (and thus occurrence of potentially large errors in extrapolation) of NN. To achieve regularized input and output values, various techniques have been explored, including the min-max method (Xie et al. 2019b; Pawar et al. 2020; Wang et al. 2021), Gaussian normalization method (Stoffer et al. 2021; Cheng et al. 2022; Liu et al. 2022), normalization with root-mean-square (r.m.s.) values (Xie et al. 2020a,b,c; Guan et al. 2022, 2023), non-dimensionalization in wall units (Park & Choi 2021; Kim et al. 2022), and prescribed velocity and length scales (MacArt et al. 2021). However, each method has its own limitations: e.g. it requires an equilibrium state or homogeneous direction(s), does not show generalizability when switched to other flows, or proves to be infeasible as it requires DNS results. Consequently, a plausible approach towards general direct NN-based SGS models would be to employ only local variables for non-dimensionalization, like those described by Prakash, Jansen & Evans (2022) and Abekawa et al. (2023). On the other hand, two strategies may be proposed for the issue of extrapolation errors of NN. The first strategy is to discover normalized flow variables and SGS stresses that are bounded within a certain range for various flows. The second strategy is to incorporate all accessible data to cover an exhaustive range of flow fields. The former presents considerable challenges, because it requires a formidable task of pinpointing normalized flow variables that remain invariant with respect to factors like the Reynolds number, filter size and flow topology. The latter may demand DNS data of various flows at fairly high Reynolds numbers and thus imposes a huge computational burden.

Therefore, the objective of the present study is to devise a new NN-based SGS model designed to overcome the limitations of existing direct NN-based SGS models. We modify the structure of an NN and formulate a recursive procedure to generate an SGS model valid for a wide range of grid sizes. Here, an NN is trained using fDNS data at low Reynolds number and then updated by data at higher Reynolds number collected through filtered LES data. To validate the performance of the present recursive NN-based SGS model, several SGS models are employed to simulate FHIT at various Reynolds numbers. In addition, two DHIT cases are simulated using the NN trained with FHIT. Section 2 provides the LES framework, NN-based SGS models under consideration, FHIT, and filtering method employed. Section 3 describes the recursive algorithm for constructing recursive NN-based SGS models, and provides the results of a priori and a posteriori tests. In § 4, the present SGS models are applied to FHIT at high Reynolds numbers and DHIT, and the results are discussed, followed by conclusions in § 5.

2. Numerical details

2.1. Large eddy simulation

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

(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 2000). 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 1963), dynamic Smagorinsky model (DSM, Germano et al. 1991; Lilly 1992), gradient model (GM, Clark et al. 1979) and dynamic mixed model (DMM, Anderson & Meneveau 1999) 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 2000). 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. 1979; Vreman et al. 1996) 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 1999) 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 (2021) 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 2015). 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 (2013) to overcome the vanishing gradient problem (Bengio, Simard & Frasconi 1994). The negative slope is determined to be $0.01 \leq a \leq 0.2$ (Xu et al. 2015), 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 2016). 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 2014) and learning rate annealing method (Simonyan & Zisserman 2014; He et al. 2016) to optimize the weights and enhance the training speed, respectively. Early stopping (Goodfellow et al. 2016) 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 2017; Wang et al. 2018; Xie et al. 2019b; Zhou et al. 2019; Pawar et al. 2020; Prat, Sautory & Navarro-Martinez 2020; Xie et al. 2020a,b,c; Park & Choi 2021; Kim et al. 2022; Abekawa et al. 2023; Kim et al. 2024), and the same form of input was considered by Jamaat & Hattori (2022) 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. 1995; Rosales & Meneveau 2005; Park et al. 2006):

(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 (2000), $\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 1994). Additionally, the spectral cutoff filtering may exhibit non-local oscillations (Meneveau & Katz 2000; Pope 2000). 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 2020; Guan et al. 2022, 2023; Pawar et al. 2023) 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. A recursive NN-based SGS model (NNVGM)

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 2021) 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 2003; Liu, Wu & Zhou 2008). 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 (2021), 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. (1991):

(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. 2005).

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. Applications to FHIT at high Reynolds numbers and DHIT

To assess the applicability of the present NNVGMs (table 5), we conduct LESs of FHIT at high Reynolds numbers and of DHIT. During simulation, a clipping is applied to the model coefficients of DSM and DMM when they are negative (see (2.9) and (2.16)), whereas NNVGMs do not use any clipping.

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 2016). 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. 1996), 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 2021). 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 (1971), 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 (2003), 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 (1971) and Kang et al. (2003) 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 1971; Kang et al. 2003). 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 (1971) (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 2010). 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. 1998) and rescaling method (Kang et al. 2003). 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. (1996). 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. (2003). 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. (2003). 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. (2003).

The LES results on DHIT by Kang et al. (2003) 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. (2003) 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. 2003) 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. (2003) 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. 2003). 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. 2003) 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 (2023) 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. 2024 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.

5. Conclusions

In the present study, we developed an NN-based SGS model (NN-based velocity gradient model; NNVGM) using a dual NN architecture (the output of one NN is the SGS normal stresses and that of the other is the SGS shear stresses) with the input of $\bar {\varDelta }^2|\bar {\alpha }|\bar \alpha _{ij}$, where $\alpha _{ij}$ is the velocity gradient tensor. We aimed at overcoming the difficulties encountered during the SGS model development, i.e. how to obtain high-Reynolds-number data for training and how to apply a trained NN to untrained flow. By eliminating bias and batch normalization, and employing the leaky ReLU function as an activation function within the NNs, the present NN retains its nonlinearity but satisfies $\boldsymbol {{\mathrm {NN}}} (c\boldsymbol {A}) = c \boldsymbol {{\mathrm {NN}}} (\boldsymbol {A})$ for a positive scalar $c$. This important property of NN allows us to apply the NN to a flow different from the trained one, because the parameters used for non-dimensionalization can be included in the scalar $c$. To generate training data, we adopted a cut-Gaussian filter rather than the Gaussian or spectral cutoff filter, considering its application to LES with coarser grids and realizability condition. We also designed a recursive procedure which consisted of the following steps: (1) conducting DNS; (2) training an NN with fDNS data; (3) conducting LES at a higher Reynolds number; (4) training the NN with augmented data including fLES data; (5) going to Step (3) for higher Reynolds number flows. For the present FHIT, the grid and filter sizes were normalized by the Kolmogorov length scale, and these normalized sizes became double at every recursive procedure. The NN trained through this recursive procedure contained a wide range of filter sizes normalized by the Kolmogorov length scale such that the LES grid size required could be located within this range of filter sizes. Testing NNVGMs on FHIT showed that fLES data can be a practical alternative to fDNS data for training an NN, and one can avoid using costly DNS to extract training data.

We conducted LESs of forced and decaying homogeneous isotropic turbulence (FHIT and DHIT) with NNVGMs and traditional SGS models, respectively. In FHIT, the same number of grid points was used for different Reynolds numbers. Among the SGS models considered, NNVGM constructed through the recursive procedure performed very well for FHIT. In contrast, NNVGM trained only with fDNS data at a low Reynolds number showed energy pile up at high wavenumbers. We considered two different DHIT, one by Comte-Bellot & Corrsin (1971) and the other at much higher Reynolds number by Kang et al. (2003). In both DHIT, the recursive NNVGM predicted the decay of resolved turbulent kinetic energy and its energy spectra in time most accurately among the SGS models considered, while most other SGS models provided excessive energy spectra at intermediate or high wavenumbers.

In the present study, we applied the NNVGM to FHIT and DHIT, but this procedure may not be applicable to laminar and inhomogeneous turbulent flows. Hence, the next step is to extend the present approach to those flows by employing the dynamic approach used in DSM or training NN with more flows having different topology. This work is being carried out in our group (Cho & Choi 2023; Kim & Choi 2023).