A wall model based on neural networks for LES of turbulent flows over periodic hills

In this section, we describe the numerical method and the case setup for generating the data employed for developing a data-driven wall model, which can take into account the non-equilibrium effects, e.g. flow separation and reattachment, for turbulent flows over periodic hills.

We employ the Virtual Flow Simulator (VFS-Wind) Yang et al. (2015b); Yang and Sotiropoulos (2018) code for WRLES of turbulent flows over periodic hills. The VFS code has been successfully applied to industrial and environmental turbulent flows Kang and Sotiropoulos (2012); Kang et al. (2014); Khosronejad and Sotiropoulos (2014); Yang et al. (2017a); Yang and Sotiropoulos (2019a, b); Foti et al. (2019); Khosronejad et al. (2020). In VFS-Wind code, the governing equations are the three-dimensional unsteady spatially filtered incompressible Navier-Stokes equations in non-orthogonal, generalized curvilinear coordinates shown as follows:

where $x_{i}$ and $\xi^{i}$ are the Cartesian and curvilinear coordinates, respectively, $\xi^{i}_{l}=\partial\xi^{i}/\partial x_{l}$ are the transformation metrics, $J$ is the Jacobian of the geometric transformation, $u_{i}$ is the $i$-th component of the velocity vector in Cartesian coordinates, $U^{i}=(\xi^{i}_{m}/J)u_{m}$ is the contravariant volume flux, $g^{jk}=\xi^{j}_{l}\xi^{k}_{l}$ are the components of the contravariant metric tensor, $\rho$ is the fluid density, $\mu$ is the dynamic viscosity, and $p$ is the pressure. In the momentum equation, $\tau_{ij}$ represents the anisotropic part of the subgrid-scale stress tensor, which is modeled by the dynamic eddy viscosity subgrid-scale model,

where $\overline{S}_{ij}$ is the filtered strain-rate tensor and $\nu_{t}$ is the eddy viscosity calculated by

where $C$ is the model coefficient calculated dynamically using the procedure of Germano et al. Germano et al. (1991), $|\overline{S}|=\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}$ and $\Delta=J^{-1/3}$ is the filter size, where $J^{-1}$ is the cell volume.

The governing equations are spatially discretized using a second-order accurate central differencing scheme, and integrated in time using the fractional step method. An algebraic multigrid acceleration along with generalized minimal residual method (GMRES) solver is used to solve the pressure Poisson equation. A matrix-free Newton-Krylov method is used for solving the discretized momentum equation. More details about the flow solver can be found in the referencesGe and Sotiropoulos (2007); Kang et al. (2011); Yang et al. (2015b).

The geometry of the periodic hills is shown in Figure 1, which has been extensively employed in experiments Breuer et al. (2009); Rapp and Manhart (2011) and numerical simulations Temmerman et al. (2003); Fröhlich et al. (2005); Krank et al. (2018). As seen, the height of the hill is $h$, with a flat wall placed $2.036h$ above the crest of the hill and the distance between the crests of two hills is $L_{x}=9h$. In the spanwise direction, the size of the computational domain is $L_{z}=4.5h$. The Reynolds number based on the bulk velocity $U_{b}$, which is defined as $U_{b}=Q/(\rho L_{z}\times(L_{y}-h))$ (where $Q$ is the mass flux), and the height of the hill is $Re_{h}=\rho U_{b}h/\mu$. No slip boundary condition is applied at the top wall and the surface of the hills. In the streamwise and spanwise directions, periodic boundary condition is applied. The flow is driven by a pressure gradient uniformly applied to whole domain to maintain a constant mass flux.

In this work, the WRLES of flow over periodic hills at three Reynolds numbers are carried out for the training and testing of data-driven wall model, as shown in Table 1. The computational domain is discretized using a body-fitted curvilinear grid (as shown in Figure 1). The height of the first off-wall grid nodes in wall units, $\Delta y^{+}=\Delta y_{cc}u_{\tau}/\nu$, is in the range of 0.056 to 3.95 at $Re_{h}=10595$ as shown in Figure 12 in Appendix A. Here, $\Delta y_{cc}$ is half height of the first off-wall grid, $u_{\tau}=\sqrt{\tau_{w}/\rho}$ denotes the friction velocity and $\nu=\mu/\rho$.

The size of time step is $\Delta t=0.01h/U_{b}$. The simulation is first carried out for about $22T$ (flow-through time $T=L_{x}/U_{b}$) for the flow to achieve a fully developed state. Then the flow is further simulated for about $50T$ for time-averaged quantities and flowfield data on slices for training the data-driven wall model.

To validate the employed numerical method and case setup, we compare the profiles of the mean velocity, Reynolds shear stress, turbulence kinetic energy, and the skin friction and pressure coefficients computed in this work with the results from measurements Temmerman et al. (2003) and DNS by Krank et al. Krank et al. (2018) and demonstrate an overall good agreement as shown in Appendix A for validating the employed numerical method and case setup.

The WRLES data are further processed to prepare the data for training the data-driven wall model. In WMLES, the wall shear stress is often computed using the velocity at the first off-wall grid node. If the data-driven model is developed using the velocity at a specific location, it may only be applicable to grids of grid spacings. Moreover, if the wall model is developed using the data at a certain streamwise location, e.g. the location where the flow is attached or the location where flow separation occurs, it may only be valid for a certain flow condition. To avoid these two issues, we are devoted to develop a data-driven wall model applicable to different spatial resolutions and not limited to certain flow conditions using the data in the near-wall region of the periodic hills at all streamwise locations. Specifically in this work, the flow data in the near-wall region with wall-normal distance in the range of $0.006\leq\eta/h\leq 0.1$ are employed, where $\eta$ denotes the wall-normal coordinate. The top boundary of the region at $\eta/h=0.1$ is determined considering that the flow field above is less correlated with the wall shear stress and is usually well-resolved by WMLES. It is noticed that the region with $\eta/h<0.1$ is defined as the inner layer for a turbulent channel flow, where $h$ is the half width of the channel. The bottom boundary at $\eta/h=0.006$ is defined to preclude the effects of viscous sublayer and with the consideration that no wall model is needed if the viscous scale is resolved.

A step-by-step diagram for preparing the training data is shown in Figure 2. Saving the three-dimensional flow fields at every time step, which requires a significant amount of disk space, is not feasible. Instead, we save the WRLES data at four spanwise ($x-y$) slices located at $z/h=0.0$, $1.125$, $2.25$ and $3.375$, respectively. To make the most of the WRLES data and meanwhile keep the cost for training the model at a reasonable level (in other words, avoid using the flowfields close in time, which can be very similar), 9 snapshots of the instantaneous flowfields are extracted on each slice for one flow-through time. And in total we obatain 450 snapshots for the whole $50$ flow-through times at each Reynolds number. For each snapshot, the flow field data at different wall-normal locations are further extracted at ${N_{x}}-1$ streamwise locations due to periodic boundary conditon to include different flow features along the lower wall of periodic hills. At each streamwise location of the lower wall, the flowfield data at 95 nodes, which are uniformly distributed in $\eta/h\in[0.006,0.1]$, are then extracted to form 95 pairs of input-output data along the wall-normal direction. Finally, the flowfield data along the wall-normal direction at different streamwise locations for all considered spanwise slices and snapshots form the complete training and testing datasets, which contain approximately $1.8\times 10^{8}$ input-output pairs. It is noticed that the grid nodes in the wall-normal ($\eta$) direction in general do not coincide with the curvilinear grid nodes employed for solving the flow. The linear interpolation approach based on triangulation is employed to obtain the flowfield data at the 95 points along the wall-normal coordinate.

Input features and output labels are critical for the success of data-driven models. The wall shear stress including the streamwise and spanwise stresses $\tau_{w,t}$, $\tau_{w,s}$, which are often applied as boundary conditions for outer flow simulations in WMLES, are employed as the output labels. To prepare the dataset for model training, the wall shear stresses are directly calculated using the WRLES data. As for the input features, the wall-normal distance $\eta$, the three velocity components $u_{w,t}$, $u_{w,n}$ and $u_{s}$ in the wall-tangential, -normal and spanwise directions and the pressure gradients $\frac{\partial p}{\partial w_{t}}$, $\frac{\partial p}{\partial w_{n}}$ in the wall-tangential and -normal directions are employed. It has been shown by Duprat et al. Duprat et al. (2011) that using a near-wall scaling defined with the classic friction velocity and the streamwise pressure gradient can improve the performance of wall models for separated flows. To take into account such knowledge when constructing the neural networks for data-driven wall models, the wall-normal distance normalized using the near-wall scale and written in the logarithmic form, i.e. $\ln\left(\eta/y^{*}\right)$, where $y^{*}=\nu/u_{\tau p}$, $u_{\tau p}=\sqrt{u_{\tau}^{2}+u_{p}^{2}}$, $u_{\tau}=\sqrt{\left|\tau_{w}\right|/\rho}$, $u_{p}=\left|(\mu/\rho^{2})(\partial P/\partial x)\right|^{1/3}$, is employed as the input feature for training the data-driven wall model. Here the wall shear stress is computed as $\tau_{w}=\mu u_{t}/\eta$, which implicitly accounts for no-slip boundary conditions at the wall and approximates the wall shear stress resolved in WMLES. To further improve the generality of the trained model, the pressure gradients are multiplied by $\frac{\eta}{h}$ before taken as input features as suggested by Yang, Bose and Moin Yang et al. (2017b).

We use a multi-hidden-layer feedforward nueral network (FNN) to establish the relation between the near-wall flow and the wall shear stress on the surface of periodic hills. As shown in Figure 3, the employed FNN consists of an input layer, multi hidden layers and an output layer. Each layer has a number of neurons, which are computational units that propagate the weighted sums of the inputs to an activation function and calculate the output. The detailed procedures for calculating the output based on the input in the FNN is shown in Appendix B, which includes the linear matrix manipulation of the weight and bias coefficients and the nonlinear mapping using the activation function.

The activation function used in this paper is the rectified linear unit (ReLU) (Goodfellow et al., 2016),

The prepared input and output data are normalized using the Min-Max scaling,

The loss function is defined as

where $N_{s}$ is the number of training samples, $w$ is the weight coefficient, and $\lambda_{0}$ is the regularization rate, which is set to 0.001. The first term in Eq. (6) denotes the mean square error (MSE) between the FNN output ${\mathbf{Y}}^{\text{*}}$ and the labeled output $\mathbf{Y}$ from the WRLES. The second term is an L2 regularization term included to avoid overfitting.

We use the error backpropagation (BP) scheme [56] implemented with TensorFlow [57] to train the FNN by optimizing the weight and bias coefficients to minimize the loss function. The key steps for FNN training are as follows:

(1) Provide training data to the input layer, propagate data signal forward layer by layer, and compute the result in the output layer. Details on this step can be found in Appendix B.

(2) Compute the loss function according to Eq. (6) using FNN output and the labeled output.

(3) Adjust the weight and bias coefficients using the gradient descent algorithm,

where $v$ denotes the weight and bias coefficients in the FNN, $\eta\in\left(0,1\right)$ denotes the learning rate, which is dynamically adjusted using the Adam optimizer Kingma and Ba (2014).

(4) Repeat the above steps until the stop criterion is met.

In this section, three sets of cases, i.e. one with different number of input features, the other one with different number of neurons in each hidden layer, and the third one with different number of output labels, are carried out to evaluate the performance of different setups for training wall model using FNN and the flow data at $Re_{h}=5600$ and $10595$. The key requirement for a wall model is to accurately predict the wall shear stress. In conventional wall models, the wall shear stress is determined by an empirical relation (e.g. the power law or the logarithmic law) or simplified Navier-Stokes equations (e.g. the thin boundary layer equation) using the velocity at one off-wall grid node (usually the first or the second off-wall node). To compensate the lack of governing equations in data-driven wall models, flowfield data at more than one off-wall grid nodes are probably needed. In the first set of cases, we test the effects of input features obtained using different numbers of wall-normal points (five inputs per point) with the distance between two adjacent points $0.03h$. How well a data-driven model represents near-wall turbulence with flow separations and reattachments depends on the complexity of the employed neural network, i.e. the number of hidden layers and the number of neurons in each layer. In the second set of cases, we examine the effects of the number of neurons ranging from 3 to 100 on the training performance with fixed number of hidden layers and input features. Details on these two sets of cases can be found in Table 2.

Figure 4 plots the variations of loss with the training epochs for both training and validation datasets for the two sets of cases. The number of training samples is $1.24\times 10^{7}$, of which $90\%$ are used as the training dataset and the rest $10\%$ are used as validation dataset, and the batch size is $2\times 10^{5}$. Intially, the loss is large because the weight coefficients are randomly set and the bias coefficients are set to zero. Then, the weight and bias coefficients are adjusted and the loss rapidly decreases during the first few epochs. After the initial stage, the loss tends to approch a steady small value after approximately 1000 training epochs.

In Figure 4(a), the loss in the FNN model I5-O1 (only using input features at one wall-normal point) is significantly worse than the other models, that the loss at 6000 epochs is about 1.5 times larger than others, while the values of loss from the models using input features at more than 1 points are similar with each other. This indicates that only using the input features at one point is not surfficient to accurately model the complex near-wall turbulence encountered in this periodic hill case, while adding just one point can significantly improve the training performance. To ensure the training performance without increasing the computational cost in the meantime, we choose 15 input features from 3 wall-normal points for case set 2 and other cases in this work. Figure 4(b) compares the loss of FNN models with different number of neurons. If no overfitting occurs, the more neurons are employed, the smaller the loss can be obtained. In this work we use 20 neurons for the proposed FNN wall model.

To consider the influence of output labels on the FNN training, the model ‘FNN-1’ with only one output label ($\tau_{w,t}$) and ‘FNN-2’ with two output labels ($\tau_{w,t}$, $\tau_{w,s}$) are trained, validated and tested, respectively. The settings of the two FNN wall models are shown in Table 3. Figure 5 shows the variations of loss with training epochs for FNN-1 and FNN-2. It can be observed that including the spanwise wall shear stress as the output label has some effects on the model training process for epochs less than 100 but little influence on the final loss of tangential wall shear stress.

To evaluate the prediction accuracy and generalization capacity of the trained FNN wall model, we apply it to predict the wall shear stress for different snapshots and spanwise slices for both training datasets and testing datasets obtained from the cases with $Re_{h}=5600$ and $10595$.

We first evaluate the FNN wall models for predicting wall shear stresses using the training dataset. In Figure 6 we show the comparison of the instantaneous friction coefficient (which is defined as $C_{f}=\tau_{w}/\frac{1}{2}\rho U_{b}^{2}$) at an instant, the correlation coefficient $\rho_{\text{FNN-LES}}$ between the instantaneous wall shear stress predicted by the FNN model and the WRLES, and the error $\varepsilon_{\tau}$ of the instantaneous wall shear stress predicted by the FNN model relative to the WRLES predictions. Here, the correlation coefficient $\rho_{\tau}$ is defined as

where $\left\langle{\;}\right\rangle$ denotes the average over snapshots.

To further assess the prediction accuracy of FNN model on the fluctuations of wall shear stress over time, we define the instantaneous relative error $\varepsilon_{\tau}$ as follows:

where $\left|\left\langle\tau_{w}^{\text{LES}}\right\rangle\right|_{\max}$ denotes the peak value of averaged wall shear stress among all the streamwise locations. The relative error of the time-averaged wall shear stress, which will be shown in Figure 7(d), is defined as follows:

As seen in Figure 6(a)$\sim$(b), the instantaneous skin friction coefficient $C_{f}$ predicted by the FNN model in general agrees with that from WRLES at most streamwise locations. Many sharp peaks are observed in streamwise variation of the instantaneous wall shear stress. The trained FNN wall model is observed being able to predict these abrupt variations at most streamwise locations, although the peak amplitude is underpredicted at some locations, e.g. at $x/h\approx 0.18$ for this instant for the $Re_{h}=10595$ case.. It is also noticed that the wall tangential component of wall shear stress predicted by the FNN-1 and FNN-2 almost collapse with each other. Moreover, the spanwise component of the wall stress predicted by the FNN-2 model is also in good agreement with the WRLES predictions. In Figure 6(c)$\sim$(d), it is observed that the correlation coefficients are larger than 0.7 and the instantaneous relative errors are smaller than 0.1 at almost all streamwise locations except near the crest of the hill at $x/h=0.2$ and $x/h=8.5$ (where the correlation coefficient is around 0.6, the instantaneous relative error is around 0.15), indicating that the large temporal fluctuations there are not well captured by the FNN model. It is noticed that correlation coefficient for the spanwise wall shear stress is similar with that of the tangential wall shear stress although the magnitude of the instantaneous spanwise wall shear stress is one order of magnitude smaller than the tangential component, which makes it difficult to train the corresponding FNN model for both components.

In Figure 7, we evaluate the capability the FNN model in predicting the mean wall shear stress and the standard deviation of wall shear stress. As seen in Figure 7(a), the mean skin friction coefficients at both $Re_{h}=5600$ and $10595$ predicted by the FNN models (i.e. FNN-1 and FNN-2) are in perfect agreement with WRLES results at all streamwise locations. The mean spanwise wall shear stress component, on the other hand, is close to zero for both FNN-2 and WRLES predictions (not shown in the figure). In Figure 7(b)$\sim$(c), we compare the normalized standard deviations of the fluctuations of wall shear stresses predicted by the FNN models with those from WRLES. It is observed that the FNN predictions are smaller than the WRLES predictions for both tangential and spanwise components. Interestingly, it is observed that the standard deviation of the spanwise shear stress is similar with the tangential component although its instantaneous value is one order of magnitude smaller than the tangential component. Figure 7(d) shows the error of time-averaged wall shear stress predicted by the FNN model relative to that from WRLES. It is seen that the absolute values of the errors are smaller than 0.05 at most streamwise locations, except at locations close to the crest of the hill.

Overall, we have shown that both FNN-1 and FNN-2 wall models can accurately predict the instantaneous and time-averaged wall shear stresses for the training dataset. Next, we will evaluate the performance of the FNN wall models using the testing dataset, which is obtained from a different $x-y$ slice (located at $z/h=2.25$) from the cases with $Re_{h}=5600$ and $10595$ totally different from the training dataset.

In Figure 8, we first evaluate the capability of the FNN models in predicting the instantaneous wall shear stress using the testing dataset. As seen in Figure 8(a)$\sim$(b), both tangential and spanwise instantaneous skin friction coefficients predicted by the FNN models are in an overall good agreement with the WRLES predictions except for some sharp peaks. Figure 8(c)$\sim$(d) show the correlation coefficient and relative error between the FNN and LES predictions. As seen in the range of $x/h=1$ to $8.5$, the correlation coeffcients are in general larger than 0.7 and the relative errors are smaller than 0.1 for both tangential and spanwise components. At locations close to the crest of the hill, lower coefficients and larger errors are observed especially close the separation point for the spanwise component.

In Figure 9, we compare the mean skin friction coefficients (averaged over snapshots), the standard deviations of the fluctuations and the time-averaged relative error of wall shear stresses predicted by different FNN wall models. As seen in Figure 9(a), good agreements between the FNN and WRLES predictions at $Re_{h}=5600$ and $10595$ are obtained for the mean friction coefficients for both FNN models. For the standard deviations of the wall shear stresses, the predictions by the FNN models are significantly smaller than those from WRLES for both tangential and spanwise components. In Figure 9(d), the relative errors are smaller than 0.05 at most streamwise locations, which are close to the results in the training dataset.

To further evaluate the generalization capacity of the FNN wall model, we apply the trained FNN wall models to the testing dataset at $Re_{h}=19000$ in this section.

Figures 10 and 11 evaluate the predictive capability of the FNN models on the instantaneous and mean wall shear stresses, respectively. We first examine the predictons of the instantaneous skin friction coefficient. As shown in Figure 10(a), the FNN predictions are consistent with those from WRLES at most streamwise locations. The correlation coeffcients (Figure 10(b)) of instantaneous wall shear stresses between the FNN and LES predictions are also observed in general larger than 0.6 in the range of $x/h=1$ to $8.5$, although they are somewhat smaller than those computed using the training dataset. The instantaneous relative errors are then examined in Figure 10(d). It is observed that the magnitude of the errors are slightly larger than those in training dataset, but still smaller than 0.1 at most streamwise locations.

Here we examine the performance of the FNN models for predicting the mean skin friction coefficient. As seen in Figure 11(a), the mean skin friction coefficients at $Re_{h}=19000$ predicted by the FNN-2 model are in good agreement with WRLES results at all streamwise locations. Some discrepancies are observed between the FNN-1 predictions and the WRLES predictions highlighting the importance of having the spanwise wall shear stress as one of the output labels. As expected, the mean spanwise wall shear stress component is close to zero for both FNN model predictions and WRLES predictions. Similar with those observed in the cases with $Re_{h}=5600$ and $10595$, the normalized standard deviations of wall shear stress fluctuations predicted by the FNN models are smaller than those from WRLES, as shown in Figure 11(b). In Figure 11(d), the relative errors of the time-averaged wall shear stress are smaller than 0.1 at most streamwise locations, which are slightly larger than those computed using the training dataset. It is also noticed that the magnitude of the relative error for the predictions from the FNN-2 model is smaller than that from the FNN-1 model especially at locations around $x/h=8.4$.

Although further improvement is still needed for the FNN wall model to accurately predict the standard deviations of wall shear stress fluctuations, the evaluations against the testing datasets at different Reynolds numbers demonstrate the excellent generalization capacity of the trained FNN wall models.