2.1. Wall shear stress model

We empower the capability of the neural network model in estimating wall shear stress for different flow regimes using the data from separated flows and the law of the wall. Specifically, the data from the periodic hill flow and the logarithmic law for the mean streamwise velocity profile, which is introduced in Appendix A, are employed for training the wall shear stress model. The periodic hill case, even though its geometry is relatively simple, contains several flow regimes, i.e. flows with pressure gradients, flow separation and flow reattachment. The inclusion of the logarithmic law in model training then enables the trained model to react properly to flows in the equilibrium state.

A FNN is employed for building the connection between the near-wall flow and the wall shear stress. The neural network comprises an input layer, six hidden layers with 15 neurons in each layer and an output layer. The hyperbolic tangent (tanh) function serves as the activation function. The input features consist of six flow quantities at three wall-normal points, with the distance between two adjacent points denoted as $\Delta y_n/\delta _0 = 0.03$,

(2.1)\begin{equation} X_I = \left\{ \ln\left(\frac{y_n}{y^{*}}\right), \frac{u_{w,t}}{y_n} \frac{\delta_0}{u_b}, \frac{u_{w,n}}{y_n} \frac{\delta_0}{u_b}, \frac{u_s}{y_n} \frac{\delta_0}{u_b}, \frac{\partial p}{\partial w_t} \frac{y_n}{\delta_0} \frac{\delta_0}{u_b^2}, \frac{\partial p}{\partial w_n} \frac{y_n}{\delta_0} \frac{\delta_0}{u_b^2} \right\}, \end{equation}

where $y_n$ is the wall-normal distance, $u_{w,t}$, $u_{w,n}$ and $u_s$ are the velocity components in the wall-tangential, wall-normal and spanwise directions, ${\partial p}/{\partial w_t}$, ${\partial p}/{\partial w_n}$ are the pressure gradients in the wall-tangential and wall-normal directions, $u_b$ is the bulk velocity, $\delta _0$ is the global length scale. It is noted that the length scale $\delta _0$ represents the scale of the outer flow. It is set as the hill height $h$ in periodic hill flows, and the half-channel width $\delta$ in turbulent channel flows. For flows in a complex geometry with several geometric length scales, a systematic way for defining the length scale is yet to be developed. A case-by-case approach is probably needed. The wall-normal distance is normalized by a near-wall length scale $y^{*} = \nu /u_{\tau p}$ (Duprat et al. Reference Duprat, Balarac, Métais, Congedo and Brugiére2011), where $\nu = \mu /\rho$ is the kinematic viscosity, $u_{\tau p} = \sqrt {u_{v}^2+u_p^2}$, $u_{v} = \sqrt {| {\nu u_{w, t}}/{y_n} |}$, $u_p = | ({\nu }/{\rho } )({\partial p}/{\partial w_t}) |^{1/3}$. The output labels are the wall-tangential and spanwise wall shear stresses,

(2.2)\begin{equation} Y_O^{(w)} = \left\{ \frac{\tau_{w, t}}{u_b^2}, \frac{\tau_{w, s}}{u_b^2} \right\}. \end{equation}

The cost function is defined as the mean square error between the predicted output and the real output. The error backpropagation scheme (Rumelhart, Hinton & Williams Reference Rumelhart, Hinton and Williams1986) implemented with TensorFlow (Abadi et al. Reference Abadi2016) is employed to train the $\text {FEL}_{\tau _w}$ model by optimizing the weight and bias coefficients to minimize the cost function. A detailed description of the training procedures can be found in Zhou et al. (Reference Zhou, He and Yang2021a).

2.2. Eddy viscosity model

The SGS modelling error and discretization error are the two major causes for the suboptimal a posteriori performance of a data-driven wall model in a WMLES environment. In separated flows, a near-wall layer of constant shear stress does not exist. In WMLES, the near-wall layer, which is in a non-equilibrium state for separated flows, is not resolved by the grid. The wall shear stress boundary condition alone is not enough to model the way the unresolved near-wall flow affects the outer flow. In this work, we postulate that a proper eddy viscosity can approximate the effects of the near-wall flow not captured by the wall shear stress, and accounts for the SGS modelling error and discretization error in the meantime.

There are several options to establish an eddy viscosity model, e.g. an analytical approach, a data-driven approach or a hybrid approach. To compromise between generalizability and predictive capability of the model, a hybrid approach is employed. Specifically, we employ a modified version of the mixing length model with its coefficients learned in an embedded way, in the following form:

(2.3)\begin{equation} \nu_t = \exp\left( E_\nu \frac{u_p}{u_{\tau p}} \right) (\kappa y_n)^2 \left| \bar{S} \right| D \text{ at the first off-wall grid nodes,} \end{equation}

where the von Kármán constant $\kappa \approx 0.4$, $D = [ 1 - \exp (-(y^+ / A^+)^3) ]$, $A^+=25$, $y^+ = y_n u_\tau /\nu$, $u_\tau = \sqrt {\tau _w/\rho }$ is the friction velocity given by the wall shear stress model in § 2.1. With (2.3), the eddy viscosity at the first off-wall grid node is modelled using the modified mixing length model instead of the dynamic Smagorinsky model (DSM) employed in the outer flows. The exponential form coefficient, i.e. $\exp ( E_\nu ({u_p}/{u_{\tau p}}))$ is proposed to prevent negative eddy viscosity. While a clipping procedure for negative $\nu _t$ can be employed as well (Zhou et al. Reference Zhou, He, Wang and Jin2019; Park & Choi Reference Park and Choi2021), it may affect the predictive capability of the model.

To approximate the modified mixing length model using a neural network, the model coefficient $E_\nu$ is used as the output label. The neural network is composed of an input layer, six hidden layers with 15 neurons in each layer and an output layer. In the input layer, the input features are the same as those for the wall shear stress model, i.e. (2.1). The hyperbolic tangent (tanh) function is also employed as the activation function for this new neural network. Using the ensemble Kalman method, the $\text {FEL}_{\nu _t}$ model is trained in the WMLES environment in an embedded way.

2.3. Ensemble Kalman method

The ensemble Kalman method is a statistical inference method based on Monte Carlo sampling, and has been widely used in various applications (Chen et al. Reference Chen, Chang, Meng and Zhang2019; Zhang et al. Reference Zhang, Xiao, Gomez and Coutier-Delgosha2020; Schneider, Stuart & Wu Reference Schneider, Stuart and Wu2022; Zhang et al. Reference Zhang, Xiao, Wu and He2022b; Liu, Zhang & He Reference Liu, Zhang and He2023). In this work, the method is employed to learn the weights in the neural network model for eddy viscosity. It utilizes the statistics of the weights and model predictions to compute the gradient and Hessian of the cost function to update the model. The cost function is given as

(2.4)\begin{equation} J= \| w_m^{n+1} - w_m^n \|_{\mathsf{{P}}}+ \| \mathcal{H}[w_m^{n+1}] - \mathsf{{y}} \|_{\mathsf{{R}}}, \end{equation}

where $w_m$ is the weight of neural network, $m$ is the sample index, $n$ is the iteration index, $\mathsf{{P}}$ is the error covariance of neural network model, $\mathsf{{R}}$ is the error covariance of observation data, $\mathsf{{y}}$ is the observation data that obey the normal distribution with zero mean and variance of $\mathsf{{R}}$ and $\mathcal {H}$ is the model operator that maps the neural network weights to model predictions.

The method uses the ensemble of the neural network samples $W$ ($=[w_1, w_2, \ldots, w_M]$) to estimate the sample mean $\bar {W}$ and covariance $\mathsf{{P}}$ as

(2.5)\begin{equation} \left.\begin{aligned} \bar{W} & = \frac{1}{M} \sum_{m=1}^M w_{m} {,} \\ {}\mathsf{{P}} & = \dfrac{1}{M-1} (W - \bar{W})(W-\bar{W})^\top {,} \end{aligned} \right\} \end{equation}

where $M$ is the sample size. Based on the Gauss–Newton method, the first- and second-order derivatives of the cost function are required to update the weights. The ensemble Kalman method uses the statistics of these samples to estimate the derivative information (Luo et al. Reference Luo, Stordal, Lorentzen and Nævdal2015). At the $n{\rm th}$ iteration, each sample $w_m$ is updated based on

(2.6)\begin{equation} w_m^{n+1} = w_m^n + {\mathsf{{PH}}}^\top ({\mathsf{{HPH}}}^\top + {\mathsf{{R}}})^{{-}1} ({\mathsf{{y}}}_m^n - {\mathsf{{H}}}w_m^n) {,} \end{equation}

where ${\mathsf{{H}}}$ is the tangent linear model operator. The readers are referred to Zhang et al. (Reference Zhang, Xiao, Luo and He2022a) for details of the employed ensemble Kalman method.

2.4. Procedure for embedded training

The procedure for embedded training of the wall model is divided into four steps as depicted in figure 2, including (i) pretrain the model; (ii) obtain the predictions for an ensemble of neural network models; (iii) compute the error of the model predictions; (iv) update the neural network using the ensemble Kalman method; (v) iterate steps (ii) to (iv) until a certain threshold or the maximum number of iterations is reached. Specifics for steps (ii)–(iv) are listed as follows.

1. (i) Pretrain the model: the neural network model for the eddy viscosity is trained using the mixing length model with $E_{\nu }=0.0$ ($\exp (E_\nu ) = 1.0$) to obtain the initial guess on the neural network weights. This is important for accelerating the embedded training of the model. The wall shear stress model is trained using the high-fidelity data and the law of the wall, and remains unchanged during the training of the eddy viscosity model.

2. (ii) Obtain the WMLES predictions for an ensemble of neural network models: this is done by applying the FEL model to WMLES. In this step, it is important to consider the balance between the computational cost and training efficiency when determining the number of samples.

3. (iii) Compute the errors of the model predictions: the WMLES predictions and the observation data are employed to compute the error. The quantities and their locations essentially determine the capability of the learned model, which has to be selected with care.

4. (iv) Update the neural networks using the ensemble Kalman method: the errors obtained in the last step are employed to update the weights of the neural network model for $\nu _t$ based on the ensemble Kalman method.

Figure 2. Embedded training of the FEL wall model. EnKF denotes the ensemble Kalman filter, BCs denotes the boundary conditions.

3.1. The WMLES environment

The VFS-Wind (virtual flow simulator) code is employed for WRLES and WMLES (Yang et al. Reference Yang, Sotiropoulos, Conzemius, Wachtler and Strong2015b; Yang & Sotiropoulos Reference Yang and Sotiropoulos2018; Zhou, Wu & Yang Reference Zhou, Wu and Yang2021b; Zhou et al. Reference Zhou, Li, He and Yang2022). The governing equations are the spatially filtered incompressible Navier–Stokes equations in non-orthogonal, generalized curvilinear coordinates as follows:

(3.1)\begin{align} \left. \begin{aligned} J \dfrac{ \partial U^{j} }{\partial \xi^{j}} & = 0, \\ \dfrac{1}{J} \dfrac{\partial U^{i}}{\partial t} & = \dfrac{\xi^{i}_{l}}{J} \left( -\dfrac{\partial}{\partial \xi^{j}} (U^{j}u_{l}) - \dfrac{1}{\rho} \dfrac{\partial}{\partial \xi^{j}} \left( \dfrac{\xi^{j}_{l} p}{J} \right) + \dfrac{\mu}{\rho} \dfrac{\partial}{\partial \xi^{j}}\left(\dfrac{g^{jk}}{J} \dfrac{\partial u_{l}}{\partial \xi^{k}} \right) - \dfrac{1}{\rho} \dfrac{\partial \tau_{lj}}{\partial \xi^{j}} + f_l \right), \end{aligned} \right\} \end{align}

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 SGS stress tensor, which is modelled using the Smagorinsky model (Smagorinsky Reference Smagorinsky1963),

(3.2)\begin{equation} \tau_{ij} - \tfrac{1}{3} \tau_{kk} \delta_{ij} ={-} 2 \nu_{t} \overline{S_{ij}}, \end{equation}

where the overbar denotes the grid filtering, and $\overline {S_{ij}} = \frac {1}{2} ( {\partial \overline {u_i}}/{\partial x_j} + {\partial \overline {u_j}}/{\partial x_i} )$ is the filtered strain-rate tensor. The eddy viscosity $\nu _t$ is calculated by

(3.3)\begin{equation} \nu_{t} = (C_S \varDelta)^2 |\bar{S}|, \end{equation}

where $C_S$ is the Smagorinsky coefficient, $|\bar {S}| = \sqrt {2\overline {S_{ij}}\,\overline {S_{ij}}}$ and $\varDelta = J^{-1/3}$ is the filter size, where $J^{-1}$ is the cell volume.

We employ the DSM (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992) for modelling the subgrid scales, in which the model coefficient $C_S$ is calculated to minimize the mean square error between the resolved stress at the grid filter and the test filter,

(3.4)\begin{equation} C_S^2 = \frac{\left\langle L_{ij}M_{ij} \right\rangle}{\left\langle M_{kl}M_{kl} \right\rangle}, \end{equation}

where $L_{ij} = \widehat {\overline {u_i u_j}} - \widehat {\overline {u_i}}\,\widehat {\overline {u_j}}$, $M_{ij} = 2\varDelta ^2 \widehat {\overline {S_{ij}}|\bar {S}|} - 2\hat {\varDelta }^2 \overline {S_{ij}} |\hat {\bar {S}}|$; the hat denotes the test filtering which involves the 27 grid nodes surrounding a given grid node in three dimensions. Here $\hat {\varDelta }$ is the size of test filter. For the generalized curvilinear coordinates employed in this work, the following formulation is employed (Armenio & Piomelli Reference Armenio and Piomelli2000):

(3.5)\begin{equation} C_S^2 = \frac{\left\langle L_{il}M_{im}G_{lm} \right\rangle}{\left\langle M_{kp}M_{kq}G_{pq} \right\rangle}, \end{equation}

where $G_{lm}$ is the covariant metric tensor. The homogeneous direction can be further employed for computing $C_S$. In all the simulated cases, only the local averaging is employed. No noticeable differences were observed when employing additional averaging in the spanwise direction.

The governing equations are spatially discretized using a second-order accurate central difference scheme, and integrated in time using the fractional step method. An algebraic multigrid acceleration along with a GMRES (generalized minimal residual method) solver is used to solve the pressure Poisson equation. A matrix-free Newton–Krylov method is used for solving the discretized momentum equation.

3.2. Embedded training

The ensemble Kalman method is an efficient method for embedded learning. To learn the model, an ensemble of samples are required to estimate the gradient and Hessian based on the statistics derived from these samples (Luo et al. Reference Luo, Stordal, Lorentzen and Nævdal2015). A sample corresponds to a WMLES case with the FEL wall model, which runs over a time period long enough to obtain the statistics to be compared with the reference data (the WRLES results in this work). The ensemble of WMLES samples are executed in parallel without communication during running the cases. With the error computed based on the discrepancy between the WMLES predictions and reference data, the FEL model is updated for the next iteration.

The periodic hill flow case employed for embedded learning of the FEL wall model is illustrated in figure 3. The hill geometry is given by analytical expressions. The different slopes shown in figure 3(a) are obtained by multiplying a factor $\alpha$ to the hill width (Xiao et al. Reference Xiao, Wu, Laizet and Duan2020; Zhou et al. Reference Zhou, He and Yang2021a). The computational domain ($L_x=9.0h, L_y=3.036h, L_z=4.5h$) and the employed curvilinear mesh on an $x$–$y$ plane are illustrated in figure 3(b) for the case with $\alpha = 1.0$, with the contour of time-averaged streamwise velocity with streamlines from the simulation at $Re_h=\rho u_b h/\mu =10\,595$ (where $h$ represents the hill height, $u_b = Q/(\rho L_z (L_y-h))$ denotes the bulk velocity and $Q$ is the mass flux). A periodic boundary condition is applied in the streamwise and spanwise directions. A uniform pressure gradient maintaining a constant mass flux is applied over the entire domain to drive the flow. On the top wall and the surface of the hills, a no-slip boundary condition is applied in WRLES, while in WMLES the FEL wall model is employed.

Figure 3. Schematic of (a) the periodic hill geometry with various slopes (parameterized by multiplying a factor $\alpha$ to the hill width), and (b) the computational domain ($L_x=9.0h, L_y=3.036h, L_z=4.5h$), the employed curvilinear mesh on an $x\unicode{x2013}y$ plane (every fifth grid) and the contour of time-averaged streamwise velocity with streamlines for the baseline case ($\alpha = 1.0$) at $Re_h=10\,595$.

Table 1 lists the parameters of the WRLES (‘WR’) and WMLES (‘WM’) cases for the flow over periodic hill with different slopes, grid resolutions and Reynolds numbers. In the table, ‘H0.5’, ‘H1.0’ and ‘H1.5’ represent the periodic hill configurations with $\alpha =0.5$, 1.0 and 1.5, respectively. The parameters $N_x$, $N_y$ and $N_z$ denote the grid number in the streamwise, vertical and spanwise directions, respectively. Here $\Delta x_f$ corresponds to the grid size in the streamwise direction at the crest of the hill, $\Delta y_f$ represents the height of the first off-wall grid cell and $\Delta y_c^+ = (\Delta y_f/2) u_\tau /\nu$ is the dimensionless wall-normal distance of the volume centre of the first off-wall grid. As the friction velocity $u_\tau$ varies with the streamwise locations, the $\Delta y_c^+$ and $\Delta y_{c,{max}}^+$ listed in table 1 are, respectively, the approximate average value and the maximum value over various streamwise locations. For the grids with $\Delta y_f/h = 0.03$, 0.06 and 0.09, the grid numbers are the same. The grid nodes are non-uniformly distributed in the vertical direction to adjust the sizes of the first off-wall grid cells. The maximum grid spacings $\Delta y/h$ are approximately 0.142, 0.107 and 0.105 for $\Delta y_f/h = 0.03$, 0.06 and 0.09, respectively. It is noted that the resolution of the finest WMLES grid employed for the periodic hill case with $Re_h=10\,595$ is close to that for a WRLES. For the case at $Re_h=37\,000$, on the other hand, the grid resolutions with $\delta y_f = 0.06$, 0.09 are typical for WMLES. The objective is to examine whether the proposed model can work when varying the grid resolution from wall-modelled to that close to wall-resolved.

Table 1. Parameters for the WRLES and WMLES periodic hill with different slopes and Reynolds numbers.

The FEL model is trained using the H0.5-WM-0.06 and H1.5-WM-0.06 cases at $Re_h=10\,595$. The vertical profiles of the mean streamwise velocity at $x/h=[0.05, 0.5, 1, 2, \ldots, 6]$ for the H0.5-WR case and $x/h=[0.05, 0.5, 1, 2, \ldots, 10]$ for the H1.5-WR case are employed as the observation data to learn the weights of the embedded neural network. Here we set the sample size $M=20$, and the maximum iteration index $N=20$. The training employs 40 samples, including 20 samples of H0.5-WM cases and 20 samples H1.5-WM cases, respectively. Each case requires approximately 27 h to run on 64 central processing unit (CPU) cores for the totally 20 iterations. Thus, the computational cost for training the model amounts to 68 266 core hours, which is primarily attributed to the WMLES solver. In comparison, the off-line FNN wall models in our previous work (Zhou et al. Reference Zhou, Yang, Zhang and Yang2023b) were trained using a graphic processing unit of NVIDIA RTX2080 for approximately 3 h, or the Intel-i7 CPU for approximately 11 h. The embedded training exhibits a higher demand on the computing environment. As for the large-scale problems which require a much finer grid resolution for WMLES, it is a great challenge to train the wall model using the embedded learning approach because of the enormous computational resources, especially the large number of CPU cores in parallel.

The obtained FEL wall model is then applied to the cases shown in table 1. The computational cost using the proposed model is found comparable to the algebraic wall model, such as the Werner–Wengle (WW) model (Werner & Wengle Reference Werner and Wengle1993). Results from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are shown in the below and compared with those from the WW-DSM model for the H0.5-WM-0.06 and H1.5-WM-0.06 cases employed for training. Here, ‘$\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$’ and ‘WW-DSM’ are named as the combination of two submodels, i.e. wall shear stress model and eddy viscosity model at the first off-wall grid nodes. Figure 4 displays the contours of time-averaged streamwise velocity with streamlines. How the hill slope affects the separation and reattachment points, and the shape of the separation bubble is demonstrated in the wall-resolved results as shown in figures 4(a) and 4(b). The $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model captures such variations of the separation bubble with the hill slope as shown in figures 4(c) and 4(d). The WW-DSM model, on the other hand, underestimates the size of the separation bubble (figure 4e,f), which was also demonstrated in the literature (Temmerman et al. Reference Temmerman, Leschziner, Mellen and Fröhlich2003; Breuer et al. Reference Breuer, Kniazev and Abel2007; Duprat et al. Reference Duprat, Balarac, Métais, Congedo and Brugiére2011). Quantitative evaluations of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are shown in figure 5 for the training cases. For both cases, an overall good agreement with the wall-resolved data is observed for the time-averaged streamwise velocity ($\langle u \rangle$), vertical velocity ($\langle v \rangle$) and the primary Reynolds shear stress ($\langle u'v' \rangle$) for the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. The TKE ($k = \frac {1}{2}\langle u'u'+v'v'+w'w' \rangle$), on the other hand, is somewhat underpredicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. This is reasonable because of the limited range of scale resolved by the employed coarse grid. We will show later in Appendix B that the TKE predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model agrees well with the filtered WRLES results.

Figure 4. Contours of time-averaged streamwise velocity with streamlines obtained from the two training cases at $Re_h = 10\,595$, (a) H0.5-WR; (b) H1.5-WR ($\Delta y_f/h=0.003$); (c) H0.5-WM, $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$; (d) H1.5-WM, $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ ($\Delta y_f/h=0.06$); (e) H0.5-WM, WW-DSM; (f) H1.5-WM, WW-DSM ($\Delta y_f/h=0.06$).

Figure 5. Vertical profiles of (a,e) the time-averaged streamwise velocity $\langle u \rangle$ and (b,f) vertical velocity $\langle v \rangle$, (c,g) primary Reynolds shear stress $\langle u'v' \rangle$ and (d,h) turbulence kinetic energy (TKE) $k = \frac {1}{2}\langle u'u'+v'v'+w'w' \rangle$ from the WRLES and WMLES with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-DSM models for the H0.5 case ($\textit {a}$–$\textit {d}$) and H1.5 case ($\textit {e}$–$\textit {h}$) with $\Delta y_f/h=0.06$ at $Re_h = 10\,595$.

To further quantify the prediction accuracy of the proposed model, we introduce the relative error for the flow statistics as follows:

(3.6)\begin{equation} {\rm err}_f = \frac{\sum |\,f_{WR}-f_{WM}| \Delta y}{\sum |\,f_{WR}| \Delta y}, \end{equation}

where $f$ represents the flow statistics, such as $\langle u \rangle$ and $k$, $\sum$ denotes the integral along the vertical direction. The obtained relative errors are shown in figure 6 for the time-averaged streamwise velocity $\langle u \rangle$ and TKE $k$ for the training cases. For the H0.5-WM-0.06 case, the errors of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are approximately 5 % for $\langle u \rangle$ and 25 % for $k$. As for the H1.5-WM-0.06 case, the errors of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are approximately 10 % for $\langle u \rangle$ and 20 % for $k$. The errors of the WW-DSM model predictions are much larger than those from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model.

Figure 6. The relative errors of (a) time-averaged streamwise velocity $\langle u \rangle$ and (b) TKE $k$ between the WMLES and WRLES for the H0.5 and H1.5 cases with $\Delta y_f/h=0.06$ at $Re_h = 10\,595$.

4.1. Cases with different grid resolutions

The proposed model was trained on a specific grid resolution. In this section, we examine the performance of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model for different grid resolutions with $\Delta y_f/h=0.01$, $0.03$ and $0.09$. The contours of the time-averaged streamwise velocity with streamlines obtained from the three grid resolutions different from the training cases are presented in figures 7 and 8. It is seen that with the refining or coarsening of the grids, consistent results are obtained for both grids. Surprisingly, with the coarsest grid (with only around 11 grid cells over the hill height), not bad recirculation bubbles are still predicted using the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model for both cases, which are significantly underpredicted or even not captured using the WW-DSM model.

Figure 7. Contours of time-averaged streamwise velocity with streamlines obtained from the H0.5-WR and H0.5-WM cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-DSM models at $Re_h = 10\,595$: (a) H0.5-WR, $\Delta y_f/h=0.003$; (b) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.01$; (c) WW-DSM, $\Delta y_f/h=0.01$; (d) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.03$; (e) WW-DSM, $\Delta y_f/h=0.03$; (f) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.09$; (g) WW-DSM, $\Delta y_f/h=0.09$.

Figure 8. Contours of time-averaged streamwise velocity with streamlines obtained from the H1.5-WR and H1.5-WM cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-DSM models at $Re_h = 10\,595$: (a) H1.5-WR, $\Delta y_f/h=0.003$; (b) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.01$; (c) WW-DSM, $\Delta y_f/h=0.01$; (d) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.03$; (e) WW-DSM, $\Delta y_f/h=0.03$; (f) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.09$; (g) WW-DSM, $\Delta y_f/h=0.09$.

Quantitative assessments of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model for different grid resolutions are demonstrated in figure 9(a–d), it is seen that the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model performs well in predicting the turbulence statistics at grid resolutions of $\Delta y_f/h=0.01$ and 0.03 for the H0.5 case, but shows some discrepancies at the coarsest grid with $\Delta y_f/h=0.09$. Regarding the assessment for the H1.5 case shown in figure 9(e–g), the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model generally predicts the vertical profiles of the turbulence statistics at the coarser grid resolutions with $\Delta y_f/h=0.03$ and 0.09. A certain degree of discrepancy is observed at the finest grid with $\Delta y_f/h=0.01$ for the second-order turbulence statistics (i.e. $\langle u'v' \rangle$ and $k$).

Figure 9. Vertical profiles of (a,e) time-averaged streamwise velocity $\langle u \rangle$ and (b,f) vertical velocity $\langle v \rangle$, (c,g) primary Reynolds shear stress $\langle u'v' \rangle$ and (d,h) TKE $k$ from the WRLES and WMLES with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model for the H0.5 case ($\textit {a}$–$\textit {d}$) and H1.5 case ($\textit {e}$–$\textit {h}$) at $Re_h = 10\,595$.

The relative errors of vertical profiles between the WMLES and WRLES are shown in figure 10 for the time-averaged streamwise velocity $\langle u \rangle$ and TKE $k$ for the H0.5 and H1.5 cases at $Re_h = 10\,595$. As shown in figures 10(a) and 10(c) for the H0.5 case, the errors for the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are only 2 % for $\langle u \rangle$ and 10 % for $k$ at grid resolution with $\Delta y_f/h=0.01$, and increase to approximately 3 % for $\langle u \rangle$ and 17 % for $k$ at grid resolution with $\Delta y_f/h=0.03$. For the coarsest grid resolution with $\Delta y_f/h=0.09$, the errors further increase to approximately 10 % and 20 % for $\langle u \rangle$ and $k$, respectively. The errors of the WW-DSM model predictions are much larger than those from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. For the assessment using the H1.5 case as shown in figures 10(b) and 10(d), an overall better performance is observed for the proposed model as well in comparison with the WW-DSM model. For the two coarser grids, the errors for $\langle u \rangle$ and $k$ are approximately 5 % and 15 %, respectively, for the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. For the fine grid with $\Delta y_f/h=0.01$, the errors for $\langle u \rangle$ are less than 5 %, but the errors for $k$ increase anomalously to 24 %.

Figure 10. The relative errors of ($\textit {a}$–$\textit {b}$) time-averaged streamwise velocity $\langle u \rangle$ and ($\textit {c}$–$\textit {d}$) TKE $k$ between the WMLES and WRLES for the H0.5 case (a,c) and H1.5 case (b,d) at $Re_h = 10\,595$.

The mean skin friction coefficients ($C_f = \tau _w / \frac {1}{2} \rho u_b^2$) predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are compared with WRLES results for different grid resolutions in figure 11. It can be observed that the skin friction coefficients predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are in good agreement with the WRLES results at most streamwise locations for the H0.5 case with $\Delta y_f/h=0.01$ and 0.03 and the H1.5 case with $\Delta y_f/h=0.03$ and 0.09, respectively. For the H0.5 case with $\Delta y_f/h=0.09$ and the H1.5 case with $\Delta y_f/h=0.01$, the peak values of the skin friction coefficient near the hill crest are somewhat overestimated.

Figure 11. Comparison of the time-averaged skin friction coefficient from the WMLES with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model and the WRLES of the (a) H0.5 and (b) H1.5 case at $Re_h = 10\,595$.

It is noted that the WMLES prediction on a coarser grid is more accurate for the H1.5 case (with the error comparison shown in figure 10d). As several factors can affect the prediction accuracy of WMLES, such as grid resolution, discretization errors, errors from the SGS models and the employed wall models, it is hard to identify the cause for the observed grid inconsistency. Such inconsistency was also reported in the literature. For instance, Zhou & Bae (Reference Zhou and Bae2024a) observed non-monotonic convergence in capturing the separation bubble of a two-dimensional (2-D) Gaussian-shaped bump for the anisotropic minimum dissipation model (Rozema et al. Reference Rozema, Bae, Moin and Verstappen2015) and Vreman model (Vreman Reference Vreman2004). It is acceptable considering that no similarity property was incorporated into the design of the input and output quantities of the model to ensure its applicability to different grid resolutions, and cases with different grid resolutions were not included in the embedded model training, either.

4.2. Cases with different hill slopes

In this section, the performance of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model is evaluated using the H1.0 case, which has a slope different from the training cases.

Figure 12 shows the contours of time-averaged streamwise velocity with streamlines obtained from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model for the H1.0 case with different grid resolutions. It is seen that the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-DSM models have similar predictions of the separation bubble for the grid with $\Delta y_f/h=0.01$. For the coarser grids with $\Delta y_f/h=0.03$ and 0.09, the separation bubble predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model is close to those from the WRLES, with large discrepancies observed for the WW-DSM model.

Figure 12. Contours of time-averaged streamwise velocity with streamlines obtained from the H1.0-WR and H1.0-WM cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-DSM models at $Re_h = 10\,595$: (a) H1.0-WR, $\Delta y_f/h=0.003$; (b) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.01$; (c) WW-DSM, $\Delta y_f/h=0.01$; (d) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.03$; (e) WW-DSM, $\Delta y_f/h=0.03$; (f) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.09$; (g) WW-DSM, $\Delta y_f/h=0.09$.

Figure 13 compares the vertical profiles of the flow statistics obtained from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model with the WRLES results for the H1.0-WM-0.01/0.03/0.09 cases with $Re_h = 10\,595$ for various streamwise locations. Good agreements with the WRLES results are observed for all the three grid resolutions, although somewhat discrepancies are observed for TKE $k$ in the lower half of the vertical profiles.

Figure 13. Vertical profiles of (a) time-averaged streamwise velocity $\langle u \rangle$ and (b) vertical velocity $\langle v \rangle$, (c) primary Reynolds shear stress $\langle u'v' \rangle$ and (d) TKE $k$ from the H1.0-WR case, the H1.0-WM cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model at $Re_h = 10\,595$.

The relative errors in $\langle u \rangle$ and $k$ predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model and the WW-DSM model are shown in figure 14. As observed, the error is less than 5 % (sometimes even 2 %) for $\langle u \rangle$, and the maximum error is around 20 % for $k$, respectively, for the cases with different grid resolutions for the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. The errors for the WW-DSM model, on the other hand, are significantly higher in comparison with the proposed model.

Figure 14. Relative errors of (a) time-averaged streamwise velocity $\langle u \rangle$ and (b) TKE $k$ between the WMLES and WRLES for the H1.0 case at $Re_h = 10\,595$.

The mean skin friction coefficients $C_f$ predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are compared with WRLES predictions in figure 15 for the three grid resolutions. It is evident that the predictions from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model closely align with the WRLES results across the three grid resolutions.

Figure 15. Comparison of the time-averaged skin friction coefficient from the WMLES with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model and the WRLES of H1.0 case at $Re_h = 10\,595$.

4.3. Cases with different Reynolds numbers

In this section, the generalization ability of the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model is tested for a higher Reynolds number (i.e. $Re_h = 37\,000$) using the H1.0 case.

The contours of time-averaged streamwise velocity obtained from the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are compared with WRLES in figure 16. It is seen that the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model successfully predicts the separation bubble for the three grid resolutions for this higher Reynolds number case. On the other hand, the WW-DSM model significantly underpredicts the size of the separation bubble.

Figure 16. Contours of time-averaged streamwise velocity with streamlines obtained from the H1.0-WR and H1.0-WM cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-DSM models at $Re_h = 37\,000$: (a) H1.0-WR, $\Delta y_f/h=0.003$; (b) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.01$; (c) WW-DSM, $\Delta y_f/h=0.01$; (d) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.03$; (e) WW-DSM, $\Delta y_f/h=0.03$; (f) $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\Delta y_f/h=0.09$; (g) WW-DSM, $\Delta y_f/h=0.09$.

The vertical profiles of the flow statistics predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are shown in figure 17 for the high Reynolds number case. Good agreement with the WRLES results is observed for all the three grid resolutions. The relative errors in $\langle u \rangle$ and $k$ are shown in figure 18 for the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and the WW-DSM models. As seen, the errors in $\langle u \rangle$ are less than 10 %, and the errors in $k$ are less than 20 %, respectively, for the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. In comparison, the errors of the WW-DSM model predictions are relatively larger.

Figure 17. Vertical profiles of the (a) time-averaged streamwise velocity $\langle u \rangle$ and (b) vertical velocity $\langle v \rangle$, (c) primary Reynolds shear stress $\langle u'v' \rangle$ and (d) TKE $k$ from the H1.0-WR and H1.0-WM cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model at $Re_h = 37\,000$.

Figure 18. Relative errors of (a) time-averaged streamwise velocity $\langle u \rangle$ and (b) TKE $k$ between the WMLES and WRLES for H1.0 case at $Re_h = 37\,000$.

The mean skin friction coefficients $C_f$ predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model are compared with the WRLES predictions in figure 19, with the difference $\Delta C_f = C_f (Re_h=10\,595) - C_f (Re_h=37\,000)$ calculated to show the Reynolds number effect. It is seen that the magnitudes of both the maximum $C_f$ at the windward side of the hill ($x/h \approx 8.5$) and the minimum $C_f$ in the recirculation zone ($x/h \approx 2\sim 5$) decrease while the magnitude of the second maximum $C_f$ close to the separation point increases with the increase of Reynolds number. The $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model succeeds in predicting the variation of the maximum $C_f$ with the Reynolds number, but does not predict the variations with Reynolds numbers at other locations.

Figure 19. (a) The time-averaged skin friction coefficient from the WMLES with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model and the WRLES of H1.0 case at $Re_h = 37\,000$ and (b) the difference $\Delta C_f = C_f (Re_h=10\,595) - C_f (Re_h=37\,000)$.

5.1. Comparisons of predictions from different submodel combinations

To understand the roles of the proposed wall shear stress model and eddy viscosity model, we perform the WMLES cases with different combinations of the two submodels, which include the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model, the $\text {FEL}_{\tau _w}$-DSM model, the $\text {FEL}_{\tau _w}$-$\text {zero}_{\nu _t}$ model with the zero eddy viscosity ($\nu _t=0$) at the first off-wall grid nodes, the $\text {FEL}_{\tau _w}$-SM model, the $\text {FEL}_{\tau _w}$-LDM model and the WW-$\text {FEL}_{\nu _t}$ model. Here, ‘SM’ denotes the Smagorinsky model with constant coefficient $C_S^2 = 0.01$, ‘LDM’ denotes the Lagrangian dynamic model proposed by Meneveau, Lund & Cabot (Reference Meneveau, Lund and Cabot1996),

(5.1)\begin{equation} C_S^2 = \frac{\mathcal{J}_{LM}}{\mathcal{J}_{MM}}. \end{equation}

The values of $\mathcal {J}_{LM}$ and $\mathcal {J}_{MM}$ are updated from time step ‘$n$’ to ‘$n+1$’ at the grid point $\boldsymbol {x}$,

(5.2) \begin{equation} \left.\begin{aligned} & \mathcal{J}_{MM}^{n+1}(\boldsymbol{x}) = \varepsilon [M_{ij}M_{ij}]^{n+1}(\boldsymbol{x}) + (1-\varepsilon) \mathcal{J}_{MM}^n (\boldsymbol{x} - \tilde{\boldsymbol{u}}^n \Delta t), \\ & \mathcal{J}_{LM}^{n+1}(\boldsymbol{x}) = H\{\varepsilon [L_{ij}M_{ij}]^{n+1}(\boldsymbol{x}) + (1-\varepsilon) \mathcal{J}_{LM}^n (\boldsymbol{x} - \tilde{\boldsymbol{u}}^n \Delta t)\},\\ & \varepsilon = \frac{\Delta t/T^n}{1 + \Delta t/T^n} \ (T^n = 1.5\varDelta (\mathcal{J}_{LM}^n \mathcal{J}_{MM}^n)^{{-}1/8}), \end{aligned} \right\}\end{equation}

where $\Delta t$ is the time step, $H\{x\}$ is the ramp function ($H\{x\}=x$ if $x\geqslant 0$, and zero otherwise).

Figure 20 compares the vertical profiles of the flow statistics obtained from the H1.0-WR case and H1.0-WM-0.09 cases with $Re_h = 10\,595$. It is seen that the $\text {FEL}_{\tau _w}$ model with the SM, DSM and LDM SGS models without additional treatments at the first off-wall grid nodes cannot accurately predict the separation bubble and flow statistics. With the learned eddy viscosity submodel, the predictions of flow statistics at various streamwise locations are significantly improved for the WW-$\text {FEL}_{\nu _t}$ model, although the $\text {FEL}_{\nu _t}$ submodel is not trained with the WW model. As for the comparison between the $\text {FEL}_{\tau _w}$-$\text {zero}_{\nu _t}$ and $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ models, the results are similar with the predicted mean velocities in the lower region of the separation bubble somewhat better predicted by the $\text {FEL}_{\nu _t}$ model.

Figure 20. Vertical profiles of (a) time-averaged streamwise velocity $\langle u \rangle$ and (b) vertical velocity $\langle v \rangle$, (c) primary Reynolds shear stress $\langle u'v' \rangle$ and (d) TKE $k$ from the H1.0-WR case and the H1.0-WM-0.09 cases for different submodel combinations. The periodic hill case with $Re_h = 10\,595$ is employed for testing.

After showing the improvements in the vertical profiles of flow statistics for the WW-$\text {FEL}_{\nu _t}$ model, it is interesting to further examine how well the skin friction and pressure coefficients are predicted. It should be noted that the skin friction coefficient is calculated from the wall shear stress models (i.e. $\text {FEL}_{\tau _w}$ and WW). The pressure coefficient is defined as $C_p = (p - p_{ref}) / \frac {1}{2} \rho u_b^2$, where $p_{ref}$ is the reference pressure at $x/h=0$ on the bottom wall (Zhou et al. Reference Zhou, Wu and Yang2021b). As seen in figure 21(a), the peak value of $C_f$ is overestimated at the leeward hill face while underestimated at the windward hill face for the WW-$\text {FEL}_{\nu _t}$ model, which, on the other hand, is accurately predicted by the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model. As for the flow separation, it is delayed by the WW-$\text {FEL}_{\nu _t}$ model ($x_{sep}/h = 0.57$ (0.22 for WRLES)), resulting an early occurrence of the reattachment ($x_{ret}/h = 3.78$ (4.35 for WRLES)). Without the treatment of eddy viscosity at the first off-wall grid nodes, the predictions of $C_f$ from the $\text {FEL}_{\tau _w}$-DSM and WW-DSM models are further from the WRLES results. In figure 21(b), it is shown that both the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ and WW-$\text {FEL}_{\nu _t}$ models accurately predict the variations of $C_p$ along the streamwise direction, especially the rise immediately downstream of the hill crest. As for the $\text {FEL}_{\tau _w}$-DSM and WW-DSM models, discrepancies are observed in the $C_p$ curves especially around the hill top and at the leeward side. Moreover, it is noticed that the favourable pressure gradients at the windward side ($x/h \in [7.5, 8.8]$) are higher when using the DSM at the first off-wall grid nodes in comparison with the $\text {FEL}_{\nu _t}$ model. At last, the vortex structures, which are visualized using the isosurfaces of instantaneous pressure fluctuations ($p'/\rho U_b^2=-0.02$) are examined in figure 22 for different model combinations. It is seen that the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model can generate large-scale coherent flow structures similar to the WRLES results. The $\text {FEL}_{\tau _w}$-DSM and WW-DSM models, on the other hand, fail to generate such flow structures.

Figure 21. Comparison of the time-averaged (a) skin friction and (b) pressure coefficients between the WRLES, the WMLES with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$, $\text {FEL}_{\tau _w}$-DSM, WW-DSM and WW-$\text {FEL}_{\nu _t}$ models for the H1.0 case at $Re_h = 10\,595$.

Figure 22. Vortex structures identified using isosurfaces of instantaneous pressure fluctuation ($p'/\rho U_b^2=-0.02$) from the H1.0-WR case and the H1.0-WM-0.09 cases for different submodel combinations. The instantaneous streamwise velocity is used for colouring the isosurfaces: (a) H1.0-WR; (b) H1.0-WM, $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$; (c) H1.0-WM, $\text {FEL}_{\tau _w}$-DSM; (d) H1.0-WM, WW-DSM.

The above analyses have shown that the eddy viscosity model at the first off-wall grid nodes, affecting the overall flow pattern, plays an important role on the improvement of WMLES predictions. In the next section, the SGS stresses are examined to probe into the underlying physics for different SGS models.

5.2. Examination of SGS stresses and energy transfer

The proposed eddy viscosity model, i.e. the $\text {FEL}_{\nu _t}$ submodel is assessed by comparing the computed SGS stresses and energy transfer with those from the WRLES results and the DSM. Eddy viscosity in the first off-wall grid cell models the interaction of the unresolved near-wall flow structures with the outer flow. To compute the SGS quantities on the WMLES grid using the WRLES results, the flow fields obtained from the H1.0-WR case (with grid number $296 \times 192 \times 186$) are spatially filtered onto a grid with the resolution of $148 \times 16 \times 31$ in the streamwise, spanwise and vertical directions, respectively. By filtering the WRLES flow field, the contravariant counterpart of the SGS stress tensor in the curvilinear coordinates is obtained in the following form (Armenio & Piomelli Reference Armenio and Piomelli2000):

(5.3)\begin{equation} \sigma_i^k = \overline{J^{{-}1} \xi_j^k u_j u_i} - \overline{J^{{-}1} \xi_j^k u_j}\, \overline{u_i} = \overline{U_k u_i} - \overline{U_k}\, \overline{u_i}. \end{equation}

As for the eddy viscosity model, the SGS stress tensor is computed by

(5.4)\begin{equation} \tau_{ij} ={-}2 \nu_t \overline{S_{ij}}. \end{equation}

In figure 23, the scatter distribution of the SGS stresses at the first off-wall grids computed from the $\text {FEL}_{\nu _t}$ model and DSM are compared with those from the filtered WRLES. To compute values shown by the scatter points, the moving average of 10 continuous snapshots of flow fields, which cover $1/90$ flow-through time, is employed. It is seen that the magnitudes of the different component SGS stresses are confined in a narrow region for the $\text {FEL}_{\nu _t}$ model being in agreement with the filtered WRLES results, which, on the other hand, exhibits over a much larger range for the DSM.

Figure 23. The SGS stresses at the first off-wall grids in the vertical from the moving average of 10 continuous snapshots of flow fields: (a) $\tau _{xx}$; (b) $\tau _{xy}$; (c) $\tau _{xz}$; (d) $\tau _{yy}$; (e) $\tau _{yz}$; (f) $\tau _{zz}$. Here, $U_b^2$ is used for non-dimensionalization.

The energy transfer rate to the residual motions is another important quantity to examine, which is given by

(5.5)\begin{equation} P_r ={-} \tau_{ij} \overline{S_{ij}}. \end{equation}

Figure 24 shows the scatter distribution of energy transfer rate $P_r$ at the first off-wall grid nodes. For the filtered WRLES flow field, the energy transfer shows both positive and negative values, but it is always positive (sometimes zero) for the $\text {FEL}_{\nu _t}$ model and the DSM. It is observed that the energy dissipation (positive value of $P_r$) predicted by the $\text {FEL}_{\nu _t}$ model is closer to that from the filtered WRLES when compared with the DSM (which is much higher). Figure 25 plots the average energy transfer rate $P_r$ for the filtered WRLES, the $\text {FEL}_{\nu _t}$ model and the DSM, with 95 % confidence intervals annotated as semitransparent bands. The average is carried out along the spanwise direction and 10 continuous snapshots of flow fields. The overprediction of $P_r$ by the DSM is further confirmed.

Figure 24. Energy transfer rate $P_r$ at the first off-wall grids from the moving average of 10 continuous snapshots of flow fields, where $U_b^3/h$ is used for non-dimensionalization.

Figure 25. The average value of energy transfer rate $P_r$ at the first off-wall grids, with 95 % confidence intervals annotated as semitransparent bands. The average is carried out along the spanwise direction and 10 continuous snapshots of flow fields.

To further compare different eddy viscosity models at the first off-wall grid nodes, the probability distribution functions (p.d.f.s) of the SGS stresses and energy transfer rate are examined. The p.d.f.s are computed using 200 snapshots with the time interval between two adjacent snapshots $2T/9$. Figures 26 and 27 show the normalized p.d.f.s of $\tau _{ij}$ and $P_r$ at the first off-wall grids located at $x/h \approx 0.01$ and $y/h \approx 1.021$. It is seen that the p.d.f.s of the SGS stresses and the energy transfer rate predicted by the $\text {FEL}_{\nu _t}$ model are close to the filtered WRLES results, especially for the shear stresses $\tau _{xy}$, $\tau _{xz}$ and $\tau _{yz}$. Large discrepancies, on the other hand, are observed for the DSM.

Figure 26. Normalized p.d.f. of the SGS stresses at the first grid in both the streamwise and vertical directions for the filtered WRLES, the $\text {FEL}_{\nu _t}$ model and the DSM: (a) $\tau _{xx}$; (b) $\tau _{xy}$; (c) $\tau _{xz}$; (d) $\tau _{yy}$; (e) $\tau _{yz}$; (f) $\tau _{zz}$. The results are calculated from 200 snapshots of flow fields that cover a total simulation time $22T$.

Figure 27. Normalized p.d.f. of the energy transfer rate $P_r$ at the first grid in both the streamwise and vertical directions for the filtered WRLES, the $\text {FEL}_{\nu _t}$ model and the DSM.

The eddy viscosity and the model coefficients at the first off-wall grids are then examined in figure 28 for the proposed $\text {FEL}_{\nu _t}$ and DSM models (figure 28a,c) and varying grid spacings (figure 28b,d). It is seen in figure 28(a) that the eddy viscosity predicted by the $\text {FEL}_{\nu _t}$ submodel is much small explaining the not bad performance of the zero $\nu _t$ model. In figure 28(b), it is observed that the eddy viscosity monotonously decreases with the grid spacing, demonstrating a proper grid-sensitive property of the $\text {FEL}_{\nu _t}$ submodel. With the eddy viscosity, the model coefficient $C_S = \sqrt {\nu _t/| \bar {S} |} / \varDelta$ is computed and plotted in figure 28(c,d). In the DSM, the threshold value of $(C_S^2)_{max} = 0.1$ is used to avoid excessive damping of the resolvable turbulence. Behaviours similar to the eddy viscosity are observed for the model coefficients for different models and grid spacings.

Figure 28. Comparison of (a,b) the eddy viscosity $\nu _t$ and (c,d) its coefficient $C_S$ ($=\sqrt {\nu _t/| \bar {S} |} / \varDelta$) at the first off-wall grids: (a,c) the DSM and the $\text {FEL}_{\nu _t}$ model (2.3) simultaneously calculated from the H1.0-WM-0.09 case with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model; (b,d) the $\text {FEL}_{\nu _t}$ model calculated from the H1.0-WM-0.01/0.03/0.06/0.09 cases at $Re_h = 10\,595$. These results are obtained by the moving average of 10 continuous snapshots of flow fields.

Figure 29 compares the eddy viscosity coefficient $D$ (2.3) at the first off-wall grids calculated from the time-averaged flow fields for H1.0-WM-0.01/0.03/0.06/0.09 cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model at $Re_h = 10\,595$. Similar to figures 26 and 27, these results are averaged by 200 snapshots of flow fields that cover a total time of $22T$. It is observed that the magnitudes of $D$ in general increase with the first off-wall grid spacing $\Delta y_f/h$, with the streamwise-averaged value of $D=$ 0.322, 0.178, 0.0634 and 0.0041 for $\Delta y_f/h = 0.09$, 0.06, 0.03 and 0.01, respectively.

Figure 29. Comparison of the eddy viscosity coefficient $D$ (in (2.3)) at the first off-wall grids calculated from the time-averaged flow fields for H1.0-WM-0.01/0.03/0.06/0.09 cases with the $\text {FEL}_{\tau _w}$-$\text {FEL}_{\nu _t}$ model at $Re_h = 10\,595$.

Overall, we have observed that the proposed $\text {FEL}_{\nu _t}$ model significantly lowers the magnitude of $\nu _t$ at the first off-wall grid nodes, to an extent that simply setting the eddy viscosity to zero at the first off-wall grid nodes can also improve the WMLES predictions. Decreases in $\nu _t$, which are associated with the decreases in SGS stresses, decrease the mixing of the near-wall region with the outer flow and lower the available momentum therein, favouring flow separation, and thus solve the problem of the DSM with no flow separation or small separation bubbles. The detailed mechanism on how the eddy viscosity at the first-off grid nodes affects the flow evolution involves the momentum transport via turbulence, pressure gradient and others. Its investigations require a systematic study using carefully designed cases, and will be carried out in future work. A final note in this section is on the sound property of the proposed $\text {FEL}_{\nu _t}$ model, that it varies with the grid resolution as expected, although only one grid resolution is employed in its training.