Inserting machine-learned virtual wall velocity for large-eddy simulation
of turbulent channel flows

The concept of this study is mainly composed of two parts — a priori test and a posteriori test — as illustrated in Fig. 1. As described in introduction part, we aim at reducing the number of grid points in the near-wall region in an LES using a staggered grid system as much as possible, while still resolving the near-wall flow structure. However, when we apply a no-slip wall boundary condition in the wall-normal direction, a velocity at the first point from the wall may be overestimated as mentioned in the introduction, which causes non-negligible error in the mean velocity distribution. Similarly, although not illustrated in Fig. 1, all the fluctuating velocity components are also subjected to the similar discretization error. Therefore, we here consider the artificial wall slip velocity ${\bm{u}}_{w}$ to fix all the velocity components at the first point from the wall ${\bm{u}}_{y_{1}}$.

In a priori test, a machine-learning model ${\cal F}$ is constructed to estimate the artificial slip velocity ${\bm{u}}_{w}$, from the velocity on two $x-z$ cross sections near wall region $\{{\bm{u}}_{y_{1}},{\bm{u}}_{y_{2}}\}$ such that ${\bm{u}}_{w}\approx{\cal F}({\bm{u}}_{y_{1}},{\bm{u}}_{y_{2}})$, where $y_{1}$ and $y_{2}$ denote the locations of the first and second grid points used in LES. Hereafter, this artificial wall velocity is referred to as virtual wall surface velocity. For the training of machine-learning model, we use subsampled $x-z$ cross-sectional velocity data $\tilde{\bm{u}}$ obtained by a direct numerical simulation (DNS) at the same friction Reynolds number ${\rm Re}_{\tau}=180$ as that used in the target LES. The subsampling operation for the DNS data enables us to match the streamwise and spanwise grid resolutions to those in the target LES. As illustrated in Fig. 1$(b)$, the virtual wall surface velocity ${\bm{u}}_{w}$ used for a training process is prepared using the linear extrapolation from the DNS data at the first point from wall in the target LES, $y^{+}_{1}$ (where “+” denotes the wall units), and the next point in DNS, $y^{+}_{1}+\delta y^{+}$, as

Here, the subsampling operation is denoted as $\tilde{(\cdot)}$. In a priori test, a machine-learning model $\cal F$ estimates the virtual wall surface velocity $\tilde{\bm{u}}_{w}^{\rm{ML}}$ from the two cross-sectional DNS data at the first and the second points from the wall in a target LES, $\tilde{\bm{u}}_{y^{+}_{1}}^{\rm DNS}$ and $\tilde{\bm{u}}_{y^{+}_{2}}^{\rm DNS}$, as

and it is compared with the reference, $\tilde{\bm{u}}_{w}^{\rm{DNS}}$. The constructed model ${\cal F}$ is then applied to an LES in a posteriori test by using it as the boundary condition, i.e.,

As a machine-learning model, we capitalize on convolutional neural network (CNN) LBBH1998 originally developed in image recognition. The filters, trainable parameters inside the CNN, are able to handle high-dimensional data efficiently and extract key features. Thanks to its unique capability in handling high-dimensional data, the use of CNN has also been spread in the fluid dynamics field in recent years FFT2020 ; FNKF2019 ; SP2019 ; MFF2019 ; FNF2020 ; FHNMF2020 ; fukami2020sparse ; morimoto2021convolutional ; FukamiVoronoi ; matsuo2021supervised .

As shown in Fig. 2$(a)$, the basic operation of CNN is to take a summation of a Hadamard product between a designated region of input data and a trainable filter $h$. Usually, a CNN consists of several convolutional layers to build a certain relationship between inputs and outputs. In this study, velocity fields of $x-z$ cross-sections at two designated $y^{+}$ are fed into the first convolutional layer and then $q^{(1)}$ will be obtained as an output of the first layer. The procedure in obtaining $q^{(l)}$ from $q^{(l-1)}$ is repeated until $l<l_{\rm max}$, where the final output $q^{(l_{\rm max})}$ corresponds to a virtual wall surface velocity in our case. The operation inside the convolutional layer can be expressed as,

where $C={\rm floor}(H/2)$, $b_{k}^{(l)}$ is a bias, $\phi$ is an activation function, $M$ is the number of input data channels and $k$ is the number of filters (equals to number of output data channels), respectively.

The details of the proposed model are summarsized in Table 1. The weights on the filters ${\bm{w}}$ are optimized through the back propagation Kingma2014 by minimizing a cost function computed from output ${\bm{q}}^{(l_{\rm max})}$ and reference data ${\bm{q}}_{\rm ref}$, such that

We use the $L_{2}$ error as the cost function.

To clarify the advantage of nonlinear CNN, we also perform a linear regression analysis for virtual wall-surface velocity estimation nakamura2021comparison ; SH2017 . The linear regression is able to express the output as a linear map of input,

where ${\bm{P}}$, ${\bm{Q}}$, ${\bm{\beta}}$ is an input matrix, an output matrix, and an weight matrix, respectively. Since we use two velocity sectional fields $\{\tilde{\bm{u}}^{\rm DNS}_{y_{1}^{+}}$, $\tilde{\bm{u}}^{\rm DNS}_{y_{2}^{+}}\}$ as the input data, the virtual wall surface velocity estimated by the linear regression $\tilde{\bm{u}}^{\rm LR}_{w}$ is expressed as the linear sum of input velocity and weight, such that

where a single snapshot of velocity data is reshaped into a one-dimensional vector, ${\bm{\beta}_{1}}$ and ${\bm{\beta}_{2}}$ are the weight matrices. The weight matrices ${\bm{\beta}_{1}}$ and ${\bm{\beta}_{2}}$ are optimized so that the $L_{2}$ norm between the left-hand side and the right-hand side of Eq. 8 can be minimized over the trained snapshots. It can mathematically be formed as

We do not use the $L_{1}$ or $L_{2}$ penalization terms for the fair comparison to the CNN (i.e., Eq. 6) nakamura2021comparison .

We construct a machine-learning model to estimate the virtual wall surface velocity, which will be applied for LES. Let us visualize the estimated virtual wall surface velocities for each case in Fig. 3. The values underneath each contour represent the $L_{2}$ error norm normalized by the velocity fluctuation $\epsilon={||\tilde{{\bm{u}}}_{\rm{DNS}}-\tilde{{\bm{u}}}_{\rm{ML}}||_{2}}/{||\tilde{{\bm{u}^{\prime}}}_{\rm{DNS}}||_{2}}$. The linear regression cannot estimate the virtual wall surface velocity accurately in terms of both the contours and the $L_{2}$ error norm. In contrast, the flow fields estimated by the machine-learning model are in qualitative agreement with the reference DNS for all three cases.

We also investigate the ability of models using probability density function (PDF), as presented in Fig. 4$(a)$. The velocity distributions of DNS and machine-learning-based estimation are generally consistent, but the low probability events do not show good agreement with each other. This is because the present model is trained to minimize the $L_{2}$ error as stated in Eq. 6, thereby leading to output the average value of fluctuation components.

We then evaluate the estimation performance of the machine-learning model on wave space using the energy spectrum, as presented in Fig. 5$(b)$. For each grid width, the machine-learning model is able to estimate well especially the low wavenumber components. However, the overestimation can be found at the high wavenumber counterparts, which implies that machine-learning model preferentially estimates the low wave-number components. This observation is consistent with previous studies of turbulence analysis using supervised machine learning methods FFT2020b ; SSS2020 .

The dependence of the estimation ability on the amount of the training snapshots is also examined in Fig. 6. The estimation accuracy improves with increasing the number of snapshots used for training of machine-learning model in all cases. Notably, the decreasing rate of the error becomes smaller with the cases of more than 3000 snapshots. Therefore, we hereafter use machine-learning models trained with 3000 snapshots in a posteriori test.

In a posteriori test, the machine-learning model trained in a priori test is applied to the LES. We use f2py F2PY to combine a FORTRAN-based simulation codes with a python-based machine-learning module. In LES, the governing equations are the filtered and coarse-grained continuity equation and the Navier–Stokes equation,

where $(\bar{\;\cdot\;})$ represents a filter operation, and $\bar{\bm{\tau}}$ denotes the sub-grid scale (SGS) stress tensor. The size of the computational domain is the same as that of DNS, i.e., $(L_{x}\times L_{y}\times L_{z})=(4\pi\delta\times 2\delta\times 2\pi\delta)$, and the number of grid points is $(N_{x}\times N_{y}\times N_{z})=(64\times N^{{\dagger}}_{y}\times 64)$, where $N^{{\dagger}}_{y}$ represents the number of grid points in the $y$ direction. The uniform grid is used in the $y$ direction. As the grid width in the $y$ direction, we consider three cases as summarized in Table 2. The time step in the present LES is $\Delta t^{+}=6.30\times 10^{-2}$. We here use the constant Smagorinsky model Smago1963 as the baseline SGS model. In the present demonstration, four cases of LES are compared as follows;

1. 1.

   LES without wall models (case 1),

2. 2.

   LES with van Driest’s damping function Driest1957 (case 2),

3. 3.

   LES assisted with machine-learned model, but without van Driest’s damping function (case 3),

4. 4.

   LES assisted with machine-learned model and van Driest’s damping function (case 4).

The flow fields obtained by each LES are visualized using the second invariant of the velocity gradient tensor $(Q^{+}=0.005)$ CPC1990 in Fig. 7$(a)$. Note that we only compare among cases 1, 2, and 4 with $\Delta y^{+}=20.0$ because case 3 has shown an unstable behavior of the simulation due to the low accuracy of machine-learned model trained in a priori test. These visualized fields exhibit the vortex structures in a reasonable manner. The time history of bulk Reynolds number ${\rm{Re}}_{b}$ is also evaluated to investigate the correction of the simulation itself as shown in Fig. 7$(b)$. Note that the time history of case 3 with $\Delta y^{+}=20.0$ is only shown until around $t^{+}=5.00\times 10^{3}$ due to the unstable simulation as mentioned above. The averaged bulk Reynolds numbers provided by DNS are also shown as the gray line in each case for comparison. For each grid width, the bulk Reynolds number is not sufficiently corrected with case 3. On the other hand, case 4 is able to correct it with $\Delta y^{+}=\{5.00,10.0\}$. Note that such correction cannot be observed with $\Delta y^{+}=20.0$. This is likely caused by the low estimation ability of the model trained in a priori test. Moreover, mean streamwise velocity profiles are also compared in Fig. 7$(c)$. Analogous to the observation in Fig. 7$(b)$, case 4 shows its reasonable ability with $\Delta y^{+}=\{5.00,10.0\}$. Note again that the overestimation with $\Delta y^{+}=20.0$ is likely caused by the lack of estimation ability as stated above. Hence, a reasonable performance of a priori test, at least, is required for the present correction method.

To further examine the physical validity of the present LES in each case, the root-mean square (RMS) of velocity and vorticity fluctuations are summarized in Fig. 8. The results in case 4 with $\Delta y^{+}=\{5.00,10.0\}$ show closer distributions to that of the DNS compared to the other cases especially near the wall. This is likely because the SGS viscosity is corrected by applying the van Driest’s damping function, thereby leading to the correction of the RMS values near the wall. Therefore, the utilization of both the van Driest model (i.e., the physical correction for the SGS model) and the machine-learned model (i.e., the correction for the error due to discretization of boundary condition) employ well to capture near wall behavior correctly.

As mentioned above, we have used the Smagorinsky model for the present analyses in a posteriori test. We here discuss the generalizability of the proposed machine-learning-based wall model with regard to SGS models. In addition to Smagorinsky model, let us consider two SGS models; Wall-Adapting Local Eddy-viscosity model (WALE) ND1999 , and Coherent Structure Model (CSM) Kobayashi2005 . The results of time history of bulk Reynolds number and the mean streamwise velocity profile are shown in Fig. 9. The augmentation of LES can be seen with $\Delta y^{+}=\{5.00,10.0\}$ while failing with $\Delta y^{+}=20.0$, which is the same trend as the case with the constant Smagorinsky model. Therefore, the proposed machine-learning-based wall model is robust against the choice of SGS model.

Since our method relies on the training data provided by DNS as expressed in Eq. 6, of particular interest here is its capability at higher Reynolds numbers than that in its training process. Let us apply the machine-learned model trained at ${\rm{Re}}_{\tau}=180$ to the LES at ${\rm{Re}}_{\tau}=360$. The size of computational domain and grid points in the LES at ${\rm{Re}}_{\tau}=360$ are $(L_{x}\times L_{y}\times L_{z})=(2\pi\delta\times 2\delta\times\pi\delta),(N_{x}\times N_{y}\times N_{z})=(64\times N^{*}_{y}\times 64)$). The time step in the present simulation is $\Delta t^{+}=6.30\times 10^{-2}$. A non-uniform grid is applied in the $y$ direction at $y^{+}>40$ multiplied by a given stretch rate, while the uniform grid is considered at $y^{+}<40$.

The number of grid points in the $y$ direction $N^{*}_{y}$ depends on the grid width, as summarized in Table 3. We here use the constant Smagorinsky model as the SGS model.

The performance of the present model (case 4) is assessed using the time history of bulk Reynolds number and the mean streamwise velocity profile in Fig. 10. For comparison, we also present the results of the DNS at ${\rm Re}_{\tau}=360$ and case 2 which applies the van Driest function. As can be expected, the model does not work with $\Delta y^{+}=20.0$ due to the lack of the estimation ability. However, what is notable here is that the reasonable simulations can be achieved with $\Delta y^{+}=\{5.00,10.0\}$ despite that the test Reynolds number is considered an extrapolation from the training range. The present investigation suggests that we can expect a reasonable performance of a machine-learned model even at a higher Reynolds number for the grid width where the model employs well in a training Reynolds number range.

Towards practical applications of the proposed method, the generalizability of the machine learning model in terms of the domain size of the training data is preferable. Hence, our interest here is the use of a locally trained ML model for both a priori and a posteriori tests. One of the techniques for local domain training is to zoom-in and/or -out target images MG2018 . Morimoto et al. morimoto2020generalization has recently investigated the possibility of this concept and reported the effectiveness for various fluid flow data.

Let us use this zoom-in/out concept and discuss the dependence of the LES performance on the domain size in training data. We consider four cases in terms of the domain size of the training data; half $(L_{x}\times L_{z})=(2\pi\delta\times\pi\delta)$, one quarter $(\pi\delta\times{\color[rgb]{0,0,0}(}1/2{\color[rgb]{0,0,0})}\pi\delta)$, and one eighth $({\color[rgb]{0,0,0}(}1/2{\color[rgb]{0,0,0})}\pi\delta\times{\color[rgb]{0,0,0}(}1/4{\color[rgb]{0,0,0})}\pi\delta)$ compared to the global domain size. The machine-learning model is trained with each local domain. The constructed model is then applied to the global domain (i.e., $(L_{x}\times L_{z})=(4\pi\delta\times 2\pi\delta)$) in the test cases as summarized in Fig. 11. We use ${\rm{Re}}_{\tau}=180$ as the training and test Reynolds number for this demonstration.

The virtual wall surface velocities with $\Delta y^{+}=5.00$ are shown in Fig. 12. The flow fields estimated by the machine-learning models are in reasonable agreement with DNS data for all cases. On the other hand, the higher $L_{2}$ errors are shown with especially with the case of one eighth. This is likely because the training data used as the input of the machine-learning model loses the low-wave number components by using zooming-in technique, which leads to decrease the estimation accuracy, since the present model dominantly estimates the low wavenumber components as discussed in Section 3.1.

Moreover, the comparison of $L_{2}$ error in the other $y$ combinations are summarized in Fig. 13. The similar trends can be found with $\Delta y^{+}=\{10.0,20.0\}$ as well as with $\Delta y^{+}=5.00$. Therefore, we hereafter use the machine-learned models trained with the domain size of one quarter in a posteriori test.

The results of time history of bulk Reynolds number and the mean streamwise velocity profile are shown in Fig. 14. Note that the time history with the machine-learning model trained by the local domain with $\Delta y^{+}=20.0$ is shown until around $t^{+}=1.00\times 10^{4}$ due to the unstable simulation as the same reason mentioned in Section 3.2. Although the fluctuations can be found with time history of bulk Reynolds number compared to the case with the global training, the LES can be augmented with $\Delta y^{+}=\{5.00,10.0\}$, although not with $\Delta y^{+}=20.0$. Summarizing above, the present model shows the generalizability in terms of the domain size used in a training pipeline.