2.1. Wall model formalism

We start by considering algebraic EWMs to provide some context for our proposed extension. In these models, we employ the LoW ansatz in (1.1) as an approximation of the near-wall flow. This is done by matching the LoW mode locally with the wall-parallel LES velocities such that the LoW mode is applied along the local streamwise direction. As (1.1) has only one mode, i.e. the LoW mode, we can determine the coefficients ${\boldsymbol c} = (c_{x},c_{z})$ from matching information at a single off-wall location. This gives the equations

(2.1a,b)\begin{equation} c_x {LoW}(\kern0.7pt y^+) = U , \quad c_z {LoW}(\kern0.7pt y^+) = W,\end{equation}

where $y^+$ is the wall-unit-scaled height of the matching location while $U$ and $W$ are the local wall-parallel LES velocities at this height. We can then calculate the wall-shear stress ${\boldsymbol \tau _w} = (\tau _{x},\tau _{z})$ from (1.1) as follows:

(2.2)\begin{equation} {\boldsymbol \tau_w} = \nu \left. \frac{{\rm d}{\boldsymbol U}}{{\rm d}y} \right|_{y=0} = \nu {\boldsymbol c} \left. \frac{{\rm d} {LoW}}{{{\rm d} y}^+} \right|_{y=0} \left. \frac{{{\rm d} y}^+}{{\rm d} y} \right|_{y=0} = {\boldsymbol c} u_\tau , \end{equation}

where we have used that ${{{\rm d}} y}^+/{{{\rm d}} y} = u_\tau /\nu$ and ${\rm d}{LoW}/{{{\rm d}} y}^+|_{y=0}=1$ by definition. To relate the friction velocity $u_\tau$ to the coefficients ${\boldsymbol c}$, we note that EWMs assume a local state of equilibrium. The wall-shear stress amplitude is therefore given by $\tau _w = |\boldsymbol \tau _w| = u_\tau ^2$. This, together with (2.2), gives the following relation:

(2.3)\begin{equation} u_\tau = (c_x^2 + c_z^2)^{1/2}.\end{equation}

Thus, the appearance of the friction velocity in (2.1a,b) (contained in $y^+$) results in this being a nonlinear set of equations that would require an iterative method to obtain the solution. However, this requirement can be removed by using $u_\tau$ from the previous time step or from long-time averaging. In this case, (2.1a,b) can be solved directly and (2.2) is then invoked to calculate the wall-shear stress.

The proposed extension in (1.2), on the other hand, has two modes, i.e. LoW and $g$. To determine the coefficients ${\boldsymbol c}_{1}=(c_{1x}, c_{1z})$ and ${\boldsymbol c}_2=(c_{2x},c_{2z})$, matching with the LES at two off-wall locations is used

(2.4) \begin{equation} \left.\begin{array}{c@{}} c_{1x}{LoW}(\kern0.7pt y_1^+)+ c_{2x}g(\kern0.7pt y_1^+)= U_1 , \quad c_{1z}{LoW}(\kern0.7pt y_1^+)+ c_{2z}g(\kern0.7pt y_1^+)= W_1 ,\\ c_{1x}{LoW}(\kern0.7pt y_2^+)+ c_{2x}g(\kern0.7pt y_2^+)= U_2 , \quad c_{1z}{LoW}(\kern0.7pt y_2^+)+ c_{2z}g(\kern0.7pt y_2^+)= W_2 . \end{array}\right\} \end{equation}

Here, $y_1^+$ and $y_2^+$ are the wall-unit-scaled heights of the two matching locations while $U_{1,2}$ and $W_{1,2}$ are the local wall-parallel LES velocities at these heights. Note that similar to EWMs, the augmented LoW ansatz is applied along the local streamwise direction. The wall-shear stress ${\boldsymbol \tau }_w = (\tau _x , \tau _z)$ can then be calculated similar to the EWM as

(2.5)\begin{equation} {\boldsymbol \tau}_w = \nu\left. \frac{{\rm d} {\boldsymbol U}}{{{\rm d}} y}\right|_{y=0} = \nu \left( {\boldsymbol c}_1 \left. \frac{{\rm d} {LoW}}{{{{\rm d}} y}^+} \right|_{y = 0} + {\boldsymbol c}_2 \left. \frac{{\rm d} g}{{{{\rm d}} y}^+} \right|_{y=0} \right) \left. \frac{{{{\rm d}} y}^+}{{{\rm d}} y} \right|_{y=0} = ({\boldsymbol c}_1 + {\boldsymbol c}_2) u_\tau , \end{equation}

where we have again used that ${{\textrm {d}} y}^+/{{\textrm {d}} y} = u_\tau /\nu$ and $\textrm {d}{LoW}/{{\textrm {d}} y}^+|_{y=0}=1$, while $\textrm {d}g/{{\textrm {d}} y}^+|_{y=0} = 1$ is enforced through normalization of $g$. Further, as for EWMs, we assume that the wall-shear stress magnitude can be expressed through a local friction velocity such that $\tau _w = | {\boldsymbol { \tau _w}}| = u_\tau ^2$ still holds for the proposed model. Using this assumption, together with (2.5), we get

(2.6)\begin{equation} u_\tau = [(c_{1x}+c_{2x})^2+(c_{1z}+c_{2z})^2]^{1/2}.\end{equation}

Comparing with the EWM result from (2.3), we see that the structure is very similar. Further, the friction velocity $u_\tau$ in (2.4) (contained in $y_1^+,y_2^+$) results in a set of nonlinear equations, which can be handled in the same way as for the EWM. With the coefficients determined, (2.5) is then invoked to compute the wall-shear stress. It is worth noting that the model reduces to an algebraic EWM when ${\boldsymbol c}_2=0$.

Following Yang et al. (Reference Yang, Park and Moin2017), we apply additional filtration to the velocities at both matching locations. This filtration allows us to place the matching location(s) close to the wall without a log-layer mismatch, which is preferred to placing the matching location(s) farther away; since the latter incurs penalties on parallel computing, among other disadvantages. Yang et al. (Reference Yang, Park and Moin2017) also found similar behaviours using temporal and spatial filters in 2-D channel flow cases. We have chosen to only consider spatial filtering in this work since the model will be validated against highly unsteady flows, for which spatial filtering is more appropriate.

2.2. Wall model modes

We now describe the two modes, i.e. ${LoW}$ and $g$, in greater detail. The construction of the $g$ mode is covered separately in § 2.3 below, and we discuss the effect of the choice of the $g$ mode in Appendix B. Figure 3(a) shows the LoW mode. This is straightforward. The mode contains the viscous sublayer, the buffer layer and the logarithmic layer, but not the wake layer. The behaviour of the LoW mode conforms to $U^+=\log (\kern0.7pt y^+)/\kappa +B$ at sufficiently large $y^+$. Analytical expressions for the viscous sublayer and buffer layer can be found in Reichardt (Reference Reichardt1951) and Spalding (Reference Spalding1961). Figure 3(b) shows the POD-based mode $g$. The mode is approximately constant away from the wall and conforms to the no-slip condition at the wall.

Figure 3.

(a) The $LoW$ mode: the two dashed lines correspond to ${LoW} = y^+$ and ${LoW} = \log (\kern0.7pt y^+)/\kappa +B$. (b) The $g$ mode: the two dashed lines correspond to $g=y^+$ and $g=9.7$.

A basic requirement of a wall model is that it degenerates to the no-slip condition at the wall-resolved limit. Here, we discuss the behaviour of the model when the viscous sublayer is resolved. When the viscous layer is resolved, we have

(2.7)\begin{equation} {\boldsymbol U}=\left.\frac{{\rm d} {\boldsymbol U}}{{{\rm d}} y}\right|_{y=0} y,\end{equation}

at first off-wall grid points. As seen in figure 3, both the LoW mode and the mode $g$ are linear functions of $y^+$ in the viscous sublayer, so in this case, we have

(2.8)\begin{equation} {LoW} = g = y^+.\end{equation}

It follows from (2.7) and (2.8) that (2.4) reduces to

(2.9a,b)\begin{equation} ({ {\boldsymbol c}_{1}+{\boldsymbol c}_{2}})y_1^+{=}\left.\frac{{\rm d} {\boldsymbol U}}{{{\rm d}} y}\right|_{y=0} y_1, \quad ({ {\boldsymbol c}_{1}+{\boldsymbol c}_{2}})y_2^+{=}\left.\frac{{\rm d} {\boldsymbol U}}{{{\rm d}} y}\right|_{y=0} y_2 . \end{equation}

The two equations have the following solution:

(2.10)\begin{equation} {{\boldsymbol c}_{1}+{\boldsymbol c}_{2}}=\frac{\boldsymbol U}{|{\boldsymbol U}|}u_\tau. \end{equation}

Although ${\boldsymbol c}_1$ and ${\boldsymbol c}_2$ are not uniquely determined, ${{\boldsymbol c}_{1}+{\boldsymbol c}_{2}}$ is. As a result, the wall-shear stress is uniquely determined. We also notice that (2.10) is non-singular even when the wall-shear stress is 0. Hence, $\tau _w=0$ is a removable singularity ($\tau _w=0$ is a singularity because $y^+$ is not defined when $\tau _w=0$, it is removable because it does not create any mathematical singularity to the calculation of the wall-shear stress). In terms of numerical implementation, one could choose to do a minimum norm solution, which would be quite elegant. We have chosen to add a simple ‘if’ statement that allows the code to treat wall-resolved grids separately.

Finally, when implementing the wall model, we tabulate the $LoW$ mode and the $g$ mode in a lookup table. An inquiry about the $LoW$ and $g$ values at a given $y^+$ is then interpolated from the lookup table. Here, we prefer lookup tables over analytical formulation because ‘addressing’ a lookup table is more efficient than computing a closed-form expression. In computer science terms, a lookup table is an array that replaces runtime computation with a simpler array indexing operation. The tables may be precalculated and stored in static program storage, prefetched as part of a program's initialization, or stored in hardware in application-specific platforms. The process is called ‘direct addressing’. Because retrieving a value from memory is often faster than carrying out a more expensive computation or input/output operation, lookup tables can save processing time.

2.3. Construction of the g mode

Define ${\boldsymbol r} = {\boldsymbol U}-{\boldsymbol c}_1 {LoW}$ such that ${\boldsymbol r} = (r_x,r_z)$ describes the deviation from the LoW mode. The $g$ mode, or more precisely ${\boldsymbol c}_2 g$, is intended to capture ${\boldsymbol r}$. Any function, as long as it is not exactly the LoW mode, should capture part of that deviation. Here, we want a mode that captures as much energy in $\boldsymbol r$ as possible. This directly leads to POD. Specifically, we consider a 1-D scalar variant of POD. It is the POD of the fluctuating streamwise velocity component along the wall-normal direction in a 2-D turbulent channel. Further details on the numerical implementation are given in Appendix A. As POD can be performed in multiple spatial dimensions, we provide a few observations in this regard. First, if the POD analysis were to be performed in three-dimensional space, a local coupling between the wall model and LES would no longer be possible. Instead, a more global coupling over some wall-parallel plan(s) would be necessary which could be limiting for flows in complex geometries. Second, a global coupling between LES and the wall model could be problematic in the parallelization of the WMLES code.

We perform this 1-D POD for wall-normal intervals ranging from the wall to the log-layer region in a 2-D turbulent channel with $Re_\tau =5200$ (Lee & Moser Reference Lee and Moser2015). Three different heights within the log-layer are considered. The result for the first POD mode is shown in figure 4. We see that the near-wall part of the POD mode does not change and that the modes asymptote to a constant in the logarithmic layer. This gives rise to the $g$ mode in figure 3(b). Further, the POD spectra for the three different analyses are given in figure 5. In figure 5(a), it can be observed that the first and second POD modes carry roughly 50 % and 15 % of the total energy, respectively, for all three cases. This shows that the first POD mode is significantly more important than the subsequent modes. It is curious to note that the $g$ mode used here is almost identical to the laminar mode in Fowler et al. (Reference Fowler, Zaki and Meneveau2022), where the mode describes the response of the flow to an alternating pressure gradient. This suggests that the instantaneous flow in an equilibrium channel is also subjected to large instantaneous pressure gradients.

Figure 4.

(a) First POD mode for three wall-normal intervals which start from the wall and end at $y/\delta = 0.1$, $0.2$ and $0.3$. (b) Zoom in on the near-wall region.

Figure 5.

(a) The POD eigenvalue spectra for the three cases in figure 4. (b) Same as in (a) but in a double logarithmic plot.

We also comment on the choice of $Re_\tau =5200$ for the POD analysis. As WMLES is primarily aimed at high-Reynolds-number flows, it is preferable to base the POD analysis on such a case. Currently, $Re_\tau = 5200$ is the largest Reynolds number for which the instantaneous velocity field, needed for the POD analysis, is publicly available (Graham et al. Reference Graham2016). Still, an ideal wall model should be able to perform well also at low Reynolds numbers. Therefore, we have performed a simple test at different $Re_\tau$ for both the EWM (1.1) and POD-augmented WM (1.2), with the $g$ mode from $Re_\tau = 5200$ in figure 3(b), to compare them in terms of capturing the turbulent kinetic energy (TKE). We do this using DNS data from four different Reynolds numbers: $Re_\tau = 180$, 540, 1000 and 5200. The data for the two cases $Re_\tau = 180$ and $540$ are generated using a code with the same numerics as in Kim, Moin & Moser (Reference Kim, Moin and Moser1987) while the cases $Re_\tau = 1000$ and 5200 are from Lee & Moser (Reference Lee and Moser2015), Graham et al. (Reference Graham2016). We then fit the streamwise coefficients of both wall models in (1.1) and (1.2) by a least squares regression using the streamwise velocity from DNS data over the interval $0 \leq y / \delta \leq 0.1$. This interval is chosen to be consistent with the subsequent a priori and a posteriori testing for a 2-D turbulent channel in §§ 3 and 4, respectively. The results are shown in figure 6. We observe that the POD-augmented WM captures significantly more of the TKE than the EWM at all Reynolds numbers considered. In regards to the decrease in the captured TKE with increasing Reynolds number, which is observed for both wall models, this can be explained by analogy with the momentum cascade in homogeneous isotropic turbulence. Specifically, as the Reynolds number increases, smaller and smaller scales are activated in the flow which then start contributing to the total TKE. As these small scales cannot be captured by a few large-scale modes, the relative fraction of the TKE captured by such a finite set of large-scale modes will inevitably decrease as the Reynolds number increases.

Figure 6.

Normalized streamwise TKE captured by the EWM in (1.1) and POD-augmented WM in (1.2) for different $Re_\tau$. The normalization is by the TKE from DNS.

Admittedly, the POD modes in another flow will not be the same as the ones in a 2-D channel, but we argue that this is not a huge concern here given the purposes of this work. Specifically, considering POD modes from other wall-bounded turbulent flows will only be meaningful if the present approach yields more accurate results than the EWM. If this turns out to be the case, a more detailed analysis involving POD modes from other flows should be considered to determine the potential benefits of using particular modes in different flow types. Still, we note that the flow physics that governs the instantaneous flow in a channel is also present in other wall-bounded flows. Therefore, wall modelling improvements based on a mode from this case would be an encouraging sign for application in other flow types as well. Last, we comment on the choice of the flow configuration. Considering that the inner layer of a boundary layer is no different from that of a channel, POD analysis of a boundary-layer flow should give us similar results. Nonetheless, DNS data of channel flows is more extensively available than DNS data of boundary layers.

3.1. Equilibrium channel

We consider a 2-D turbulent channel with $Re_\tau = 1000$ (Graham et al. Reference Graham2016). We filter the DNS data following Graham et al. (Reference Graham2016) and Yang et al. (Reference Yang, Zafar, Wang and Xiao2019) so that the study is more faithful to WMLES than DNS. A Gaussian filter is used with a filtration length scale of $300 (\nu / u_\tau )$. This length scale is chosen to match the grid resolution used in the subsequent a posteriori investigations in § 4.2 to make the results more comparable. The filtered DNS data is then used to fit for the coefficients in (1.1) and (1.2). Note that we use the same DNS data for fitting both the EWM and POD-WM to ensure a fair comparison. Specifically, the matching locations are at $y/\delta =0.5/12$ and 1.5/12, i.e. the first and the second off-wall grid point locations in the subsequent a posteriori study in § 4.2. The wall-shear stress is computed according to (2.2) for the EWM and (2.5) for the POD-augmented WM. The wall-shear stress computed from the EWM and POD-augmented WM is compared with the wall-shear stress from both DNS and filtered DNS in figure 7. We observe the following. First, both the EWM and POD-augmented WM predict the correct mean stress. Second, the EWM captures the large-scale variations of the wall-shear stress, but it does not capture the more extreme wall-shear stress events. On the other hand, the POD-augmented WM captures the large-scale variation and some of the extreme events. For a more quantitative comparison, we note that the root-mean-square (r.m.s.) of the wall-shear stress fluctuations $\tau _{w,{rms}}^+$ is 0.42, 0.13, 0.11 and 0.20 in DNS, filtered DNS, predicted by the EWM and predicted by the POD-augmented WM, respectively. Further, we also give the r.m.s. of the error between the filtered DNS and the EWM and POD-WM predictions which are 0.09 and 0.16, respectively. Thus, we see that the r.m.s. of the error is larger for the POD-augmented WM than the EWM when the models are compared with filtered DNS data. As the POD-augmented WM contains two modes, i.e. the LoW mode and POD-based mode $g$, it is interesting to investigate the behaviour of the coefficients in front of these two modes. Figure 8 shows the contours of these coefficients and their correlation. First, we observe that the two coefficients are of the same magnitude. Second, there is a strong anticorrelation between the $LoW$ mode and the $g$ mode coefficients. This suggests that when the equilibrium LoW term deviates from the mean value, the non-equilibrium term pulls it back. We have repeated the above a priori test at other Reynolds numbers and matching locations. We see that $c_{1x}$ and $c_{2x}$ are always anticorrelated, but the correlation depends on the Reynolds number and the matching location. The exact physical mechanism behind this observation is not entirely clear to us, but this phenomenon, if it translates to WMLES, will help to stabilize the numerical simulation.

Figure 7.

(a) Wall-shear stress from DNS. (b) Wall-shear stress from filtered DNS. (c) Wall-shear stress predicted by the EWM. (d) Wall-shear stress predicted by the POD-augmented WM. The WM calculations are discussed in the text. Note that we show only part of the computational domain for presentation purposes.

Figure 8.

(a) Contours of $c_{1x}$. (b) Contours of $c_{2x}$. (c) Instantaneous realizations of $c_{1x}'$ and $c_{2x}$.

3.2. Non-equilibrium channel

The main purpose of the POD-augmented WM is to capture non-equilibrium effects that are not captured by the EWM. In this subsection, we consider the time evolution of a 2-D channel flow subjected to suddenly applied APGs. Using the nomenclature $R[Re_\tau ] A[APG/(\tau _{w,0}/\delta )]$, we consider the two cases: R1000A10 and R1000A100 with moderate and strong adverse pressure gradients. Again, we fit for the coefficients in (1.1) and (1.2) using the same matching information for both models with the matching locations being the same as those used in the a posteriori study in § 4.3. The wall-shear stress is computed according to (2.2) and (2.5). We then compare the computed wall-shear stresses with the DNS.

Figure 9(a,b) shows the R1000A10 results. We see from figure 9(a) that there is a good agreement between the LoW coefficients $c_x$ and $c_{1x}$ from the EWM and POD-augmented WM, respectively, and that the $g$ mode in the POD-augmented WM is almost entirely ‘turned off’ as $c_{2x} \approx 0$. Hence, the POD-augmented WM formulation reduces to the EWM in near-equilibrium conditions. In regards to the accuracy of the wall-shear stress predictions, we see that both the EWM and POD-augmented WM produce accurate results as seen in figure 9(b). Still, it should be observed that the POD-augmented WM performs slightly worse than the EWM for this near-equilibrium case. Next, the R1000A100 results are shown in figure 9(c,d). From figure 9(c), we see that the LoW coefficients for the EWM and POD-augmented WM no longer agree and that the $g$ mode in the POD-augmented WM is active. Specifically, for the POD-augmented WM, we observe that the LoW coefficient $c_{1x}$ remains almost constant while the $g$ coefficient $c_{2x}$ continuously decreases. Further, the wall-shear stress predictions in figure 9(d) show that the POD-augmented WM does a much better job of capturing the strong decrease in the wall-shear stress than the EWM.

Figure 9.

(a) The WM coefficients $c_x$ from the EWM and $c_{1x}$ and $c_{2x}$ from the POD-augmented WM in the case R1000A10. (b) Wall-shear stress from DNS, the EWM and the POD-augmented WM in case R1000A10. (c,d) Same as (a,b) but for case R1000A100.

4.1. Code

The POD-augmented WM is implemented in our in-house pseudospectral code LESGO. The code solves the incompressible Navier–Stokes equations in a half-channel with a slip top boundary and periodicity in both the streamwise and the spanwise directions. The flow is driven by a constant streamwise pressure gradient. The code employs a pseudospectral method in the streamwise and the transverse directions, a second-order finite difference method in the wall-normal direction and the second-order Adam–Bashforth method for time stepping. The grid is uniform in each direction. The time step size is such that the Courant–Friedrichs–Lewy number is always smaller than 0.06 (due to the Adams–Bashforth time-stepping scheme). Further, a Gaussian filter is applied to the wall-parallel velocities before they are used as input to the wall model. Last, the subgrid-scale stresses are modelled via the Lagrangian scale-dependent model (Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005). The code is well-validated (Abkar & Moin Reference Abkar and Moin2017; Yang et al. Reference Yang, Chen, Hu and Abkar2022) and further details of the code's numerics can be found in Yang & Abkar (Reference Yang and Abkar2018) and the references cited therein.

4.2. Fully developed channel

We first apply the POD-augmented WM in 2-D equilibrium channels. This exercise will serve as validation: the WM must be able to predict the LoW at all Reynolds numbers on both wall-resolved and wall-modelled grids. Table 1 shows the details of the WMLES. The nomenclature of the cases is $R[Re_\tau ]N[N_y]$. The Reynolds number is from $Re_\tau =180$ to $10^5$, the grid is from DNS-like (in case R180N48) to typical-WMLES-like (R1000N12, R1e5N12) and the domain is of the size as in Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014). The LES/WM matching locations are at the first two off-wall grid points, which are used by both wall models. The matching locations are in the viscous sublayer (in case R180N48) or the logarithmic layer (in cases R1000N12, R1e5N12, R1e5N24, R1e5N48). Following the established best practice (Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016), the matching locations are not placed in the buffer layer. Placing the matching location within the buffer layer incurs log-layer mismatch for WMLES in general – even if the matching location is not the first off-wall grid point and the velocity at the matching location is filtered. Nonetheless, placing the matching location is rarely a concern for flows at high Reynolds numbers since the buffer layer will not be resolved. This issue comes up here only because we validate against DNS data, which is limited to moderate Reynolds numbers. Figure 10(a) shows the mean flow. The profiles follow the LoW irrespective of the grid and the Reynolds number. We see some difference at high Reynolds numbers when the grid resolution is low, which may be due to the interplay between the subgrid-scale model and the wall model. Figure 10(b,c) shows the instantaneous contours of $c_{1x}$ and $c_{2x}$. The mean flow is $U=u_\tau {LoW}(\kern0.7pt y^+)$ in a channel. Hence, $c_{1x}$'s mean should be $u_\tau$ and $c_{2x}$'s mean should be $0$. The above expectation bears out in figure 10(b,c) where $c_{1x} / u_\tau = 0.96$ and $c_{2x} / u_\tau = 0.05$ when averaged over the wall-parallel directions. We note that the observed deviation from $c_{1x} = u_\tau$ and $c_{2x} = 0$ is due to instantaneous effects as they disappear upon time averaging. Further, the local instantaneous deviations in $c_{1x}$ from $u_\tau$ is a result of non-equilibrium effects, where we should see the POD-based $g$ mode. This expectation also bears out in figure 10(b,c).

Table 1.

The WMLES details. The half-channel height is used to normalize the numbers in the ‘Domain’ column. The first off-wall grid point is at $\Delta y/2$ (the half is because of the staggered grid).

Figure 10.

(a) Mean velocity profiles. The square symbols indicate the location of the first off-wall grid point. The solid lines are POD-WM results, and the dashed lines are EWM results. (b) Contours of $c_{1x}$. (c) Contours of $c_{2x}$.

4.3. Adverse pressure gradient

Next, we apply the POD-augmented WM in a channel subjected to a suddenly imposed APG. A schematic of the flow was shown in figure 2(a) already, and figure 2(c) shows the time evolution of the mean flow. Table 2 shows the WMLES details. The nomenclature is $R[Re_\tau ]A[APG/(\tau _{w,0}/\delta )]$, where the subscript 0 indicates a variable evaluated prior to the application of the APG. The channel is initially at a Reynolds number $Re_\tau =1000$. The APGs are $10(\tau _{w,0}/\delta )$ and $100(\tau _{w,0}/\delta )$ in cases R1000A10 and R1000A100, respectively. The flow decelerates, and the Reynolds number decreases. LES information at the second and the fourth off-wall grid points are fed to the POD-augmented WM so that we do not match in the buffer layer. For a fair comparison, the EWM makes use of the LES information at these two off-wall locations as well. Specifically, the EWM predicts wall-shear stress that yields the best fit of the velocity through the second and fourth off-wall grid points.

Table 2.

The WMLES details. The half-channel height is used to normalize the numbers in the ‘Domain’ column. The APG column shows the magnitude of the suddenly imposed APG. The numbers are normalized using $\tau _{w,0}/\delta$.

Figure 11(a,b) shows the flow rate as a function of time for R1000A10 and R1000A100, respectively, and we compare WMLES and DNS. The DNS data are by the present authors as well. The set-ups are similar to those in He & Seddighi (Reference He and Seddighi2015) and are not detailed here for brevity. We see from figure 11(a,b) that the flow rate is accurately predicted irrespective of the wall model.

Figure 11.

(a,b) Flow rate as a function of time: (a) R1000A10; (b) R1000A100. Here, ‘Inviscid’ corresponds to (4.1). (c,d) Wall-shear stress as a function of time: (c) R1000A10; (d) R1000A100. Here DNS, values given by the DNS, which we take as the truth; POD-WM, predictions by the WMLES that employs the POD-augmented WM for near-wall turbulence modelling; EWM, predictions by the WMLES that employs the EWM for near-wall turbulence modelling.

We explain why wall modelling is not critical to the prediction of the flow rate and the velocity itself. The flow rate is governed by the following volume-integrated $x$ momentum equation: ${\textrm {d}(\rho U_b)}/{\textrm {d}t}=-{\tau _w}/{\delta }-{\textrm {d}P}/{{\textrm {d}}\kern0.06em x}$. If the flow is inviscid, this equation becomes ${\textrm {d}(\rho U_b)}/{\textrm {d}t}=-{\textrm {d}P}/{{\textrm {d}}\kern0.06em x}$, which directly leads to

(4.1)\begin{equation} U_b=U_{b,0}-\frac{1}{\rho}\frac{{\rm d} P}{{{\rm d}}\kern 0.06em x}t .\end{equation}

The above inviscid estimate is plotted in figure 11(a,b) and is fairly accurate. In fact, the imposed adverse pressure gradient overwhelms the wall-shear stress and dominates the evolution of the velocity field. This makes wall modelling less critical to the prediction of the flow rate.

A more difficult test is the wall-shear stress. Figure 11(c,d) shows the (horizontally averaged) wall-shear stress as a function of time for R1000A10 and R1000A100, respectively. The EWM captures the wall-shear stress in case R1000A10 but grossly over-predicts the wall-shear stress in case R1000A100. The POD-augmented WM, on the other hand, captures the wall-shear stress well in both cases – although its prediction becomes less accurate as the flow approaches separation in case R1000A100. We also observe that the POD-augmented WM is slightly less accurate than the EWM for the R1000A10 case which was also observed in the a priori tests.

We have shown in § 3 that the POD-augmented wall model captures non-equilibrium effects. Here, we explain why the POD-augmented WM is more accurate than the EWM. We note that the velocity and the wall-shear stress are coupled, and the accurate prediction of one depends on the accurate prediction of the other. Nonetheless, for the two problems here, the imposed pressure gradient overwhelms the wall-shear stress and dominates the evolution of the velocity. Consequently, the task of wall modelling here simplifies to the prediction of the wall-shear stress when given the velocity. Under this simplification, a WM is accurate if it gives an accurate reconstruction of the near-wall velocity field. Figure 12(a) shows the DNS velocity profile at $tu_{\tau,0}/\delta =0.035$ in case R1000A100 as well as the reconstructions of the mean velocity profiles according to the POD-augmented WM and the EWM. We see from figure 12(a) that the POD-augmented WM provides a more accurate reconstruction of the mean flow than the EWM.

Figure 12.

(a) The DNS profiles and a priori WM reconstructions at $tu_{\tau,0}/\delta =0.035$ for case R1000A100. The two vertical lines are at $1.5/12\delta$ and $3.5/12\delta$, which correspond to the second and fourth grid points in our WMLES. The $LoW$ mode and the $g$ mode in the POD-augmented WM are also plotted. This corresponds to $c_{1x} LoW$ and $c_{2x}g$. Their contributions to the wall-shear stress are of opposite signs, which is consistent with figure 10. (b,c) Mean velocity profiles from DNS, WMLES with EWM and POD-augmented WM, (b) normalization by inner units, (c) normalization by outer units.

We can then quantify the error in the predicted wall-shear stress by computing the wall-shear stress due to the $LoW$ mode and the $g$ mode. Table 3 tabulates results. The $LoW$ mode in the EWM gives $\tau _w =0.50 u_{\tau,0}^2$, leading to a 61 % error. The LoW mode and the $g$ mode in the POD-WM give $0.62 u_{\tau,0}^2$ and $-0.24 u_{\tau,0}^2$, leading to an overall error of 23 %. Table 3 also tabulates the error in figure 11 at the same time instant. There, the errors are 116 % and 29 %. The larger error in the a posteriori test is because the error has accumulated in the WMLES.

Table 3.

Contributions of each mode to the wall-shear stress at time $tu_{\tau,0}/\delta =0.035$ in case R1000A100. The EWM does not contain a $g$ mode, and therefore its contribution to EWM's prediction is 0. The reconstruction result corresponds to figure 12, and the WMLES result corresponds to figure 11. The wall-shear stress values are normalized with $\rho u^2_{\tau,0}$. The DNS value is 0.31, which we take as the truth.

Before we conclude this section, we comment on the presentation of the results. Figure 12(b,c) show the DNS and WMLES profiles at $tu_{\tau,0}/\delta =0.035$ in case R1000A100. At this time instant, the DNS predicts $\tau _w=0.31\tau _{w,0}$. The POD-augmented WM predicts $\tau _w=0.40\tau _{w,0}$, and the EWM predicts $\tau _w=0.67\tau _{w,0}$, leading to a 116 % error in EWM's prediction and a 29 % error in POD-augmented WM's prediction. The EWM is not at all accurate. This, however, is only clear in figure 12(b), where normalization is by the inner units, not in figure 12(c), where normalization is by the outer units. Again, this is because the imposed adverse pressure gradient overwhelms the wall-shear stress, and the evolution of the velocity profile is dominated by the imposed pressure gradient for the initial period. We would also like to comment on the flow statistics we show. WMLES and LES, in general, resolve only part of the energy spectrum, and therefore they are not meant to capture turbulence quantities like the velocity r.m.s. and the energy spectrum, which is why we have chosen to focus on mean flow quantities.

4.4. Transverse pressure gradient

Last, we apply the POD-augmented WM in a channel subjected to a suddenly imposed TPG. A schematic of the flow was shown in figure 2(b) already. The channel is initially at a Reynolds number $Re_\tau =1000$. The magnitude of the TPG is $10(\tau _{w,0}/\delta )$. Table 4 shows the WMLES details. Figure 13(a,b) shows the WMLES result. We observe the following. Firstly, the POD-augmented wall model captures the initial decrease in the streamwise wall-shear stress and gives wall-shear stress predictions that agree more closely with the DNS than the other models (comparing figures 13a and 2d). Secondly, the streamwise wall-shear stress behaves differently from one realization to another, whereas the spanwise wall-shear stress behaves similarly between the different realizations. Thirdly, the POD-augmented WM gives more accurate spanwise wall-shear stress predictions than the EWM.

Table 4.

The WMLES details. The normalization is the same as in table 2. Here $N$ is the number of independent realizations.

Figure 13.

Time evolution of (a) the streamwise wall-shear stress and (b) the spanwise wall-shear stress. The DNS is by Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020). The POD-augmented WM results are ensemble averages from 10 statistically independent realizations.

Repeating the exercise in § 4.3, we can show that wall modelling is not critical to the prediction of the flow rate. We can also reconstruct the velocity profiles and explain why the POD-augmented WM outperforms the EWM. These results, however, are repetitive and therefore are not shown here for brevity.