2.2. Closure for the total stress at the wall-model height

The relaxation wall model requires specification of the total (molecular viscous $+$ turbulent) stress $\bar {{\boldsymbol {\tau }}}_{\varDelta }$ at $y=\varDelta$ as function of known LES quantities there, such as the LES velocity ${\boldsymbol U}_{LES}$ at the wall-model point. That is, we denote $\tilde {\boldsymbol u}(y=\varDelta ) = {\boldsymbol U}_{LES}$, where $\tilde {\boldsymbol u}({\boldsymbol x},t)$ is the velocity field being solved in LES. We argue that introducing a closure model for the total stress there ($\bar {{\boldsymbol {\tau }}}_{\varDelta }$) is more appropriate than closing the stress at the wall, since the wall is at a different position than the position at which the LES velocity is known. For now the simplest approach, which we shall adopt in this paper, is to use a standard ‘equilibrium’ RANS closure to relate the known velocity ${\boldsymbol U}_{LES}(y=\varDelta )$ to the total (viscous plus turbulent) shear stress $\bar {{\boldsymbol {\tau }}}_{\varDelta }$ at the same position. In this way, we can connect the new LaRTE wall model to the traditional equilibrium wall model: the latter is obtained simply by letting $T_s \to 0$, in which case the wall stress is set equal to the total stress at $y=\varDelta$. This equality of stresses would be justified under full equilibrium conditions, i.e. if the turbulence small-scale unresolved motions operated on much shorter time scales than the macroscopic variables (time-scale separation) and $T_s \to 0$ would be appropriate. However, for turbulence that lacks scale separation under quasi-equilibrium conditions with temporal and spatial variations and pressure gradient effects, the assumption that $T_s \to 0$ and equality of stresses are not formally justified. Then the relaxation equation can be solved instead.

Next, we describe the proposed closure for the stress $\bar {\boldsymbol {\tau }}_\varDelta$. We use the approach developed in Meneveau (Reference Meneveau2020) in which an equilibrium layer partial differential equation is numerically integrated and various fitting functions are developed for this solution. The model of Meneveau (Reference Meneveau2020) is expressed in the form of a friction Reynolds number that depends on a $U_{LES}$-based Reynolds number and a dimensionless pressure gradient parameter via a dimensionless fitting function

(2.19)\begin{equation} \frac{\langle u_\tau \rangle \varDelta }{\nu} = Re_{\tau\varDelta}^{pres}(Re_\varDelta, \psi_p),\quad {\rm where}\ Re_\varDelta = \frac{U_{LES} \varDelta}{\nu} \quad {\rm and} \quad \psi_p = \frac{1}{\rho} (\boldsymbol{\nabla}_h P \boldsymbol{\cdot} \hat{\boldsymbol e}_u) \frac{\varDelta^3}{\nu^2} \end{equation}

where the superscript ‘pres’ indicates that the fit includes pressure gradient dependence (unlike $Re_{\tau \varDelta }^{fit}(Re_\varDelta )$ also provided in Meneveau Reference Meneveau2020). This fit is repeated in Appendix C for completeness. Also, $\hat {\boldsymbol e}_u={\boldsymbol U}_{LES}/|{\boldsymbol U}_{LES}|$ is the unit vector in the ${\boldsymbol U}_{LES}$ direction. The pressure gradient $\boldsymbol {\nabla }_h P$ represents a steady or very low-frequency background pressure gradient, included in the full equilibrium part of the dynamics. Details of how $\boldsymbol {\nabla }_h P$ is determined in simulations are provided later, in § 4.1. The fit for $Re_{\tau \varDelta }^{pres}(Re_\varDelta,\psi _p)$ provided in Meneveau (Reference Meneveau2020) was obtained by numerically integrating the simple steady one-dimensional RANS equations that, unlike (2.4), did not include time dependence and can thus be characterized as a ‘fully equilibrium model’ (as opposed to quasi-equilibrium) assumption. The friction velocity $\langle u_\tau \rangle$ is the corresponding ‘full equilibrium’ friction velocity. Under these conditions, the full equilibrium vertically integrated momentum equation implies that $0 = -\varDelta \boldsymbol {\nabla }_h P/\rho + \bar {\tau }_\varDelta - \langle u_\tau \rangle ^2$, i.e. $\langle u_\tau \rangle ^2$ obtained from applying (2.19) represents a model for the combined $\bar {\tau }_\varDelta - \varDelta \boldsymbol {\nabla }_h P/\rho$.

Following usual practice of EQWMs, we assume that the total stress modelled by the fitted equilibrium expression is aligned with the LES velocity at the first grid point and write

(2.20)\begin{equation} \bar{\boldsymbol{\tau}}_\varDelta - \frac{1}{\rho} \varDelta {\boldsymbol \nabla}_h P = \frac{1}{2} c_{f}^{wm} U_{LES}^2 \hat{\boldsymbol e}_u = \left( Re_{\tau\varDelta}^{pres} \frac{\nu}{\varDelta}\right)^2 \hat{\boldsymbol e}_u ,\end{equation}

where $Re_{\tau \varDelta }^{pres}=Re_{\tau \varDelta }^{pres}(Re_\varDelta, \psi _p)$ is the fit. Note that the latter is related to the ‘equilibrium wall-model friction factor’ $c_{f}^{wm}$ according to $c_{f}^{wm} = 2 ({Re_{\tau \varDelta }^{pres}}/{Re_\varDelta })^2$. Finally then, the two-dimensional partial differential equation (PDE) governing the evolution of the friction–velocity vector in the LaRTE model reads as follows:

(2.21)\begin{equation} \frac{\partial {\boldsymbol{u}}_\tau}{\partial t} + {\boldsymbol V}_\tau \boldsymbol{\cdot} \boldsymbol{\nabla}_h {\boldsymbol{u}}_\tau = \frac{1}{T_s} \left[ \frac{1}{u_\tau}\left(-\frac{\varDelta}{\rho} \boldsymbol{\nabla}_h \bar{p}' + (Re_{\tau \varDelta}^{pres} \nu/\varDelta)^2 \hat{\boldsymbol e}_u \right) - {\boldsymbol u}_\tau\right] + u_\tau \frac{\delta^*_\varDelta}{\varDelta} \frac{\partial {\boldsymbol s}}{\partial t}, \end{equation}

where $\bar {p}'$ is the pressure fluctuation that excludes the very slow background pressure gradient based on $P$ and is further described in § 4.1. Also note that we have neglected the horizontal diffusion term in writing this equation (further discussion of horizontal diffusion is provided in § 4.2). The solution to (2.21) is then used to determine the quasi-equilibrium part of the wall stress $\bar {\boldsymbol {\tau }}_w=u_\tau {\boldsymbol u}_\tau$.

In practice, for flows such as channel flow or zero pressure gradient boundary layers that do not display major non-equilibrium effects, one would not expect to see noticeably different results in overall flow statistics, whether one applies the EQWM instantaneously at the wall as is usually done, or if one applies the proposed Lagrangian time relaxation, i.e. with some time delay and smoothing. However, the new formulation enables us to operationally separate the model self-consistently into a part that genuinely represents quasi-equilibrium dynamics, and a remainder for which additional modelling is needed. As introduced in § 3, the total modelled wall stress may also include a non-equilibrium component $\boldsymbol {\tau }_w^{\prime \prime }$, representing additional contributions to the wall stress that do not arise from the quasi-equilibrium dynamics encapsulated in (2.14) and (2.21). For example, in § 3 we introduce an additional term $\boldsymbol {\tau }_w^{\prime \prime }$ intended to capture the non-equilibrium laminar response to a rapidly changing pressure gradient in the viscous sublayer. But first, in the next section we present an a priori data analysis to study properties of several averaging time scales.

2.3. A priori analysis of quasi-equilibrium dynamics in channel flow

Here, we examine channel flow DNS data at $Re_\tau = 1000$ to show that the time scale $T_s$ identified in (2.15) provides a self-consistent decomposition of the flow into quasi-equilibrium and non-equilibrium (the remainder) components. DNS data were obtained from the Johns Hopkins Turbulence Database (JHTDB 2021) for the channel flow $Re_\tau =1000$ dataset (Graham et al. Reference Graham2016). The velocity data were Gaussian filtered horizontally at scale $\varDelta _x^+=\varDelta _z^+=196$, commensurate with the LES grid resolution for simulations considered later in this paper. The velocity was collected for all times available, $0\leq t u_\tau / h \leq 0.3245$, at a single point, $(x_0,z_0)$, over a wall-normal height $0\leq y^+ \leq \varDelta ^+ \approx 34$, close to the height of the first LES grid point that will serve as the wall-model height.

The DNS velocity is then temporally filtered below $y=\varDelta$ using a one-sided exponential time filter, computed as $\tilde {u}^n = \epsilon u^n + (1-\epsilon )\tilde {u}^{n-1}$, where $\epsilon = \delta t/T$, $\delta t$ is the time-step size, $n$ is the time step index and $T$ is the averaging time scale. Three different averaging time scales are considered to see which time scale is most consistent with quasi-equilibrium assumptions. Quasi-equilibrium is satisfied when the filtered velocity profile collapses to $\tilde {u} = u_{\tau,l} f(y u_{\tau,l}/\nu )$ for all time, where $u_{\tau,l}(x,z,t)$ is the local friction velocity. The local friction velocity is computed using $u_{\tau,l} = \nu /\varDelta Re_{\tau }^{fit}(U_{LES}\varDelta /\nu )$ where $Re_{\tau }^{fit}$ is the inverted law of the wall fit from Meneveau (Reference Meneveau2020) and $U_{LES} = \tilde {u}(x,y=\varDelta,z,t)$ is the filtered velocity at the wall-model height. In the limit $T \to \infty$, the velocity profile is static and thus full equilibrium is achieved. In this limit the local friction velocity tends towards the global friction velocity, computed as $u_{\tau,g} = \sqrt {(h/\rho )(-{\rm d}\langle p \rangle /{{\rm d}\kern0.7pt x})}$ where ${\rm d}\langle p \rangle /{{\rm d}\kern0.7pt x}$ is the bulk pressure gradient forcing for channel flow. Inner units normalization of the velocity profile using $u_{\tau,g}$ is shown in the top row of figure 2 whereas normalization with $u_{\tau,l}$ is shown in the bottom row of figure 2. In order, from left to right in figure 2, the averaging time scales considered are: panels (a,d) no temporal filtering, i.e. entirely local; panels (b,e) intermediate time scale $T_1=\varDelta /u_{\tau,g}$; and panels (c,f) LaRTE predicted time scale $T_s = {\rm \Delta} f({\rm \Delta} u_{\tau,g}/\nu )/u_{\tau,g}$. Note that $T_1$ is similar to what is used in Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2015), however, it is off by a factor $\theta /\kappa$ with $\theta =1$ used for their work. They also mentioned that a longer time scale may be needed, even suggesting a minimum of $\theta = 5$ which yields a coincidentally rather similar time scale as $T_s$ for $\varDelta ^+ = 34$. We stress that, here, $T_s$ is derived based on a momentum balance and does not require tuneable parameters as was the case in Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2015).

Figure 2. Velocity profiles at a single point $(x_0,z_0)$ at various times, from a priori tests from DNS data, Gaussian filtered in the horizontal directions at $\varDelta _x^+=\varDelta _z^+=196$. Different lines represent different times, lighter line colour corresponds to earlier time, separated by $t u_{\tau,g}/h=0.13$. Profiles are normalized using the global friction velocity $u_{\tau,g}$ in (a–c) while (d–f) use the local friction velocity $u_{\tau,l}$ as averaged over the same time filtering scale as used for the profile. (a,d) No time filter; (b,e) exponentially time filtered with $T_1=\varDelta /u_{\tau,g}$ time filter; (c,f) exponentially time filtered using $T_s = {\rm \Delta} f({\rm \Delta} u_{\tau,g}/\nu )/u_{\tau,g}$ consistent with the LaRTE approach. The $y^+$ dependence in the vertical axis matches that of the horizontal axis (i.e. $y^+ = y u_{\tau,g}/\nu$ for (a–c) and $y^+ = y u_{\tau,l}/\nu$ for (d–f)).

From figure 2(f) it is clear that filtering with the relaxation time scale $T_s$ most closely satisfies quasi-equilibrium assumptions, as it almost completely collapses to the law of the wall when normalized with the local friction velocity. This a priori test therefore provides justification that the relaxation towards EQWM, which responds within the relaxation time scale $T_s$ as discussed in § 2.1, is consistent with quasi-equilibrium assumptions.

4.1. Pressure gradient decomposition

In this section we discuss the various pressure gradient inputs to the model: $\boldsymbol {\nabla }_h P$, $\boldsymbol {\nabla }_h \bar {p}'$ and $\boldsymbol {\nabla }_h \tilde {p}''$. The first ($\boldsymbol {\nabla }_h P$) is used in evaluating the fully equilibrium fitted part to evaluate the turbulent stress as input to the LaRTE equation; $\boldsymbol {\nabla }_h \bar {p}'$ is the fluctuating pressure gradient input that directly affects the LaRTE dynamics. The last term $\boldsymbol {\nabla }_h \tilde {p}''$ is the forcing term for the non-equilibrium laminar response model described in § 3. We begin from the pressure gradient available from LES, which corresponds to the pressure gradient horizontally filtered to the size of the LES grid, denoted by $\boldsymbol {\nabla }_h \tilde {p}$, where $\tilde {p}=p_{LES}$ at $y=\varDelta$. We decompose it according to these three contributions

(4.1)\begin{equation} \boldsymbol{\nabla}_h \tilde{p} = \boldsymbol{\nabla}_h P + \boldsymbol{\nabla}_h \bar{p}' + \boldsymbol{\nabla}_h \tilde{p}'' = \boldsymbol{\nabla}_h \bar{p} + \boldsymbol{\nabla}_h \tilde{p}'', \end{equation}

where $\boldsymbol {\nabla }_h \bar {p}=\boldsymbol {\nabla }_h P + \boldsymbol {\nabla }_h \bar {p}'$.

The ‘fully equilibrium pressure gradient’, $\boldsymbol {\nabla }_h P$ is to be used in the fitting function to model the turbulent stress at $y=\varDelta$. It is obtained by temporal filtering $\boldsymbol {\nabla }_h \tilde {p}$ at a long time scale $n \,T_s$ where $n$ is some constant sufficiently greater than one and $T_s$ is the relaxation time scale of the LaRTE model. We thus write

(4.2)\begin{equation} \boldsymbol{\nabla}_h P = \boldsymbol{\nabla}_h \langle \tilde{p} \rangle_{nT_s}, \end{equation}

where the brackets indicate one-sided exponential time filtering and the subscript denotes the corresponding filtering time scale. The rationale for this choice is that the equilibrium time scale should be greater than the quasi-equilibrium time scale, $T_s$, such that only very slow pressure changes are included in the fitted full equilibrium model. We chose $n=3$ as a practical compromise that works well in applications to be shown later, and results appear to be quite insensitive to this choice.

The laminar Stokes layer that develops near the wall is caused by the high-frequency component (fastest changing) pressure gradient fluctuations. Therefore we define the non-equilibrium pressure gradient input to be a high-pass temporally filtered version of the pressure gradient. This is achieved in practice by subtracting from the LES pressure gradient another low-pass-filtered signal, but low-pass filtered at a high frequency. Specifically, we write

(4.3)\begin{equation} \boldsymbol{\nabla}_h \tilde{p}'' = \boldsymbol{\nabla}_h \tilde{p} - \boldsymbol{\nabla}_h \langle \tilde{p} \rangle_{t_\nu}, \end{equation}

where $\langle \cdot \rangle _{t_\nu }$ represents a temporal low-pass filter at time scale $t_\nu$. For $t_\nu$, the appropriate filtering time scale should be the diffusion time from the wall to the edge of the Stokes layer ($y=l_s$). We define this time scale to be $t_\nu \equiv l_s^2/\nu$. Then rewriting the Stokes layer thickness in inner units, the time scale becomes

(4.4)\begin{equation} t_\nu = \frac{(l_s^+)^2 \nu}{u_\tau^2}. \end{equation}

The Stokes layer is assumed to be confined to the viscous sublayer, therefore as an approximation we let $l_s^+ \approx 12$ and $u_\tau$ is obtained from the LaRTE model.

With $\boldsymbol {\nabla }_h P$ and $\boldsymbol {\nabla }_h \tilde {p}''$ so determined, the input to the LaRTE transport equation is the ‘band-pass filtered’ version of the pressure gradient equal to

(4.5)\begin{equation} \boldsymbol{\nabla}_h \bar{p}' = \boldsymbol{\nabla}_h \langle \tilde{p} \rangle_{t_\nu} - \boldsymbol{\nabla}_h \langle \tilde{p} \rangle_{3T_s}, \end{equation}

recalling that $\tilde {p} = p_{LES}$ is the pressure available from LES at the first wall-model point away from the wall.

Figure 3 shows the different time scales and wall distances considered for the wall modelling region beneath $y=\varDelta$. The laminar Stokes layer is confined to the fastest time scale, $t_\nu$, and smallest wall distance, $l_s$, considered. The quasi-equilibrium and full equilibrium regions, on the other hand, correspond to the largest time scales considered ($T_s$ and $nT_s$, respectively). Note that the LaRTE model has some high-frequency content coming from $\boldsymbol {\nabla }_h \bar {p}'$ and $\bar {\boldsymbol {\tau }}_\varDelta$, thus the blue region extends somewhat further left than $T_s$. The remaining region left in white corresponds to the turbulent portion in the wall-modelled region, below $y=\varDelta$, at scales faster than $T_s$ but slower than the viscous time scale $t_\nu$. The response of turbulence in this region (e.g. reduction of turbulent stresses due to scrambling) requires separate modelling not yet included in the present work. Note that interesting phenomena such as drag reduction due to spanwise wall oscillations or applied pressure gradients (Jung et al. Reference Jung, Mangiavacchi and Akhavan1992) would depend to a large extent on such modelling.

Figure 3. Schematic of the various time scales and wall-normal distances considered for modelling. Different coloured regions identify the corresponding wall modelling components.

4.2. Discretization of the LaRTE evolution equation

The full model embodied by (2.14) is a nonlinear PDE for ${\boldsymbol u}_\tau (x,y,t)$, with an elliptic diffusion term and advective term. Following the logic of the Lagrangian dynamic model implementation (Meneveau, Lund & Cabot Reference Meneveau, Lund and Cabot1996) and acknowledging the approximate nature of various modelling assumptions to be made, we opt for efficiency over high-order numerical accuracy in the proposed numerical implementation, while aiming to maintain the main features of the model. We discretize the LaRTE evolution equation using a forward Euler method such that the friction velocity vector may be solved for explicitly. This requires evaluating all terms, including the Lagrangian time derivative, at the previous time step $t_{n-1} = t_n - \delta t$ (with $n$ the time-step index) where all terms are known. We then propose to discretize the Lagrangian derivative at time $t_{n-1}$ using a semi-Lagrangian scheme (Staniforth & Côté Reference Staniforth and Côté1991)

(4.6)\begin{align} \left[\frac{\partial {\boldsymbol u}{_\tau}}{\partial t} + {\boldsymbol V}_\tau \boldsymbol{\cdot} \boldsymbol{\nabla}_h {\boldsymbol u}{_\tau} \right](x^\prime_i,z^\prime_k,t_{n-1}) &= \left[\frac{{{\rm d}}_s {\boldsymbol u}_{\tau}}{{\rm d} t}\right] (x^\prime_i,z^\prime_k,t_{n-1}) \nonumber\\ &\approx \frac{1}{\delta t} \left[ {\boldsymbol u}{_\tau}(x_i,z_k,t_n) - {\boldsymbol u}{_\tau}(x^\prime_i,z^\prime_k,t_{n-1}) \right], \end{align}

where $x^\prime _i = x_i-V_{\tau x} \delta t$ and $z^\prime _k = z_k-V_{\tau z} \delta t$ (with $i$ and $k$ position indices), and ${\boldsymbol V}_\tau$ is evaluated from (2.12) using ${\boldsymbol u}_\tau (x_i,z_k,t_{n-1})$. Replacing into the equation $[ {\boldsymbol u}{_\tau }(x_i,z_k,t_n) - {\boldsymbol u}{_\tau }(x^\prime _i,z^\prime _k,t_{n-1}) ] /\delta t = {\boldsymbol{RHS}}(x^\prime _i,z^\prime _k,t_{n-1}),$ where ${\boldsymbol {RHS}}$ is the entire right-hand side of (2.21) yields

(4.7)\begin{equation} {\boldsymbol u}{_\tau}(x_i,z_k,t_n) = {\boldsymbol u}{_\tau}(x^\prime_i,z^\prime_k,t_{n-1}) + \delta t {\boldsymbol{RHS}}(x^\prime_i,z^\prime_k,t_{n-1}). \end{equation}

The entire right-hand side of (4.7) at the upstream position $(x_i-V_{\tau x} \delta t,z_k-V_{\tau z} \delta t)$ at the time $t_{n-1}$ is obtained using first-order bilinear spatial interpolation of the grid values on the plane, as was done in three dimensions in Meneveau et al. (Reference Meneveau, Lund and Cabot1996). The additional numerical diffusion associated with the low-order interpolation reduces the need to include the horizontal diffusion term ${\boldsymbol \nabla }_h \boldsymbol {\cdot } \boldsymbol{\mathsf{D}}_\tau$ which would require additional modelling and numerical cost associated with solving an elliptic problem. Thus, we neglect the term ${\boldsymbol \nabla }_h \boldsymbol {\cdot } \boldsymbol{\mathsf{D}}_\tau$ altogether in practical implementations in our LES (however, see Appendix B for explicit form of one part of this term). Finally, evaluation of ${\boldsymbol{RHS}}$ requires the $\partial _t \boldsymbol s$ term. It is discretized using backward differencing, as $\partial _t \boldsymbol s|_{n-1} = (\boldsymbol s|_{n-1}- \boldsymbol s|_{n-2})/\delta t$ all evaluated at the interpolated position $(x_i-V_{\tau x} \delta t,z_k-V_{\tau z} \delta t)$.

4.3. Evaluating the temporal convolution integral

Jiang et al. (Reference Jiang, Zhang, Zhang and Zhang2017) developed a method for ‘fast evaluation of the Caputo fractional derivative’ which significantly reduces storage and computational cost requirements thus making the computation of the convolution integral practical. To summarize, their method decomposes the integral into a local and history parts where the history contribution is evaluated efficiently by making a sum-of-exponentials approximation to the kernel. An exponential kernel has the advantage that the current value of the convolution depends only on the previous time-step value of the convolution and a local term, as exploited in many applications where time filtering is needed (as e.g. in Meneveau et al. (Reference Meneveau, Lund and Cabot1996) and in other instances of exponential time filtering applied in this paper).

Since the sum-of-exponentials approximation algorithm is critical for the model, we will describe here the basic details of it pertaining to our application with $\alpha = 1/2$. Our task is to find an efficient way of computing the convolution integral in (3.5). To simplify notation we let $\boldsymbol {G}(t) \equiv -\rho ^{-1}\boldsymbol {\nabla }_h \tilde {p}''(t)$. Then the non-equilibrium wall stress is given by

(4.8)\begin{equation} \boldsymbol{\tau}_w''(t_n) = \sqrt{\nu/{\rm \pi}} \int_{t_0}^{t_n} \boldsymbol{G}(t')(t_n-t')^{{-}1/2} \,{\rm d} t', \end{equation}

where $t_n$ is the current time and $n$ is, as before, the time step index. The sum-of-exponentials (SOE) approximation for the kernel reads

(4.9)\begin{equation} (t_n - t')^{{-}1/2} \approx \sum_{m=1}^{N_{exp}} \omega_m \exp({-s_m (t_n - t')}), \end{equation}

where the constants $\omega _m$ and $s_m$ are determined a priori as a function of the time-step size for the SOE approximation, $\delta t$, the time duration considered, $T$, and the desired maximum error for the SOE approximation of the kernel, $\epsilon$. According to Jiang et al. (Reference Jiang, Zhang, Zhang and Zhang2017), the number of exponential terms, $N_{exp}$, is also a function of these parameters and can be estimated by the expression

(4.10)\begin{equation} N_{\exp }={O}\left(\log \frac{1}{\epsilon}\left(\log \log \frac{1}{\epsilon}+\log \frac{T}{\delta t}\right)+\log \frac{1}{\delta t}\left(\log \log \frac{1}{\epsilon}+\log \frac{1}{\delta t}\right)\right). \end{equation}

Here, $\epsilon$ is the error associated with the approximation in (4.9) (not to be confused with the error from discretizing (4.8)). The simulations in this paper will mostly use $\delta t \sim 4\times 10^{-4}$ and so the constants were computed using $\delta t=4\times 10^{-4}$. Also, it was found that, in order to guarantee good accuracy in all cases considered, we required $\epsilon =10^{-9}$. We used $T=1$ although this parameter was seen to affect the coefficients very little as long as $T \gg \delta t$. The optimization approach by Jiang et al. (Reference Jiang, Zhang, Zhang and Zhang2017) yields $N_{exp}=48$ although fewer terms (obtained by using larger $\epsilon$) could be used while still yielding reasonable accuracy. Appendix D provides information about the computed constants as well as a more detailed verification of the numerical method.

The integral is divided into local and history parts (in time)

(4.11)\begin{equation} \boldsymbol{\tau}_w'' = \boldsymbol{\tau}_{w,l}'' + \boldsymbol{\tau}_{w,h}'', \end{equation}

where

(4.12)\begin{equation} \left.\begin{gathered} \boldsymbol{\tau}_{w,l}''(t_n) \equiv \sqrt{\nu/{\rm \pi}} \int_{t_{n-1}}^{t_n} \boldsymbol{G}(t')(t_n-t')^{{-}1/2} \,{\rm d} t' \\ \boldsymbol{\tau}_{w,h}''(t_n) \equiv \sqrt{\nu/{\rm \pi}} \int_{t_0}^{t_{n-1}} \boldsymbol{G}(t')(t_n-t')^{{-}1/2} \,{\rm d} t'. \end{gathered}\right\} \end{equation}

The local part is evaluated using the ‘L1 method’ (Li & Zeng Reference Li and Zeng2015)

(4.13)\begin{equation} \boldsymbol{\tau}_{w,l}'' \approx 2 \boldsymbol{G}(t_{n-1/2}) \sqrt{\frac{\nu {\rm \Delta} t_n}{\rm \pi}}, \end{equation}

where ${\rm \Delta} t_n = t_n - t_{n-1}$ and $\boldsymbol {G}(t_{n-1/2}) = 0.5(\boldsymbol {G}(t_{n})+\boldsymbol {G}(t_{n-1}))$. The history part is evaluated by replacing the kernel with a SOE approximation from (4.9). This SOE approximation is useful because it allows the integral to be computed recursively. The history term can then be computed using

(4.14)\begin{equation} \boldsymbol{\tau}_{w,h}''\approx \sqrt{\nu/{\rm \pi}} \sum_{m=1}^{N_{exp}} \omega_m \boldsymbol{I}_m (t_n), \end{equation}

where

(4.15)\begin{align} \boldsymbol{I}_m(t_n) &= \int_{t_0}^{t_{n-1}} \boldsymbol{G}(t')\exp({-s_m(t_n-t')})\,{\rm d} t' \nonumber\\ &= \exp({-s_m(t_n-t_{n-1})}) \boldsymbol{I}_m(t_{n-1}) + \int_{t_{n-2}}^{t_{n-1}} \boldsymbol{G}(t')\exp({-s_m(t_n-t')})\,{\rm d} t' \nonumber\\ &\approx \exp({-s_m {\rm \Delta} t_n}) \left[ \boldsymbol{I}_m(t_{n-1}) + \frac{\boldsymbol{G}(t_{n-3/2})}{s_m} (1-\exp({-s_m {\rm \Delta} t_{n-1}})) \right]. \end{align}

The total non-equilibrium wall stress can then be computed using (4.11) together with (4.13), (4.14) and (4.15). The advantage of this method is that it requires $O(N_{exp})$ storage and $O(N_T N_{exp})$ computational work whereas a direct method requires storing the entire time evolution, i.e. $O(N_T)$ storage and $O(N_T^2)$ work which becomes unwieldy for long simulations.

5.1. Statistically stationary channel flow

First, the LaRTE wall model together with the non-equilibrium part is implemented in a simulation of statistically steady-state channel flow at various Reynolds numbers. This is a flow in which the traditional EQWM typically provides good results. The objective is thus mainly to ensure that similarly good results are obtained using the new model as well as to document its various features, such as typical orders of magnitudes of the terms appearing in the Lagrangian relaxation transport equation for the friction–velocity vector. Simulations use LESGO, an open-source, parallel, mixed pseudo-spectral and centred finite difference LES code available on Github (LESGO 2021). The Lagrangian scale-dependent dynamic subgrid stress model (Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005) is used in the bulk of the flow. The near-wall region is modelled using the new wall models proposed here: the LaRTE model governed by (2.14) (with closure and simplifications according to (2.21)) and the laminar non-equilibrium model governed by (3.5). A wall-stress boundary condition is applied consisting of the superposition between the two models (i.e. $\tilde {\boldsymbol \tau }_w = \bar {\boldsymbol \tau }_w + \boldsymbol {\tau }_w''$). Further notes regarding implementation are discussed in § 4 and fits needed for the LaRTE model are provided in Appendix C.

First, simulations are performed with friction Reynolds numbers based on the half-channel height of $Re_\tau = 1000$ and 5200. The domain size, number of grid points, and grid size are $(L_x, L_y, L_z)/h = (8{\rm \pi},2,3{\rm \pi} )$, $(N_x, N_y, N_z) = (128,30,48)$, and $(\varDelta _x, \varDelta _y, \varDelta _z)/h = (0.196,0.067,0.196)$, respectively. In inner units the grid sizes for $Re_\tau = 1000$ and 5200 are $(\varDelta _x^+, \varDelta _y^+, \varDelta _z^+) = (196,67,196)$ and $(\varDelta _x^+, \varDelta _y^+, \varDelta _z^+) = (1021,347,1021)$, respectively. Several additional simulations are performed at even higher Reynolds numbers ($Re_\tau = \{0.2,\ 1,\ 5\}\times 10^5$) using the same number of grid points in order to ensure applicability at arbitrarily high Reynolds numbers. As can be seen these are very coarse WMLES, very different from the much finer resolutions required for WRLES.

In LESGO the wall model takes information from the first grid point away from the wall (i.e. $\varDelta = \varDelta _y/2$). The wall-model heights for all friction Reynolds numbers considered are summarized in table 1. These wall-model heights lie within the log layer. The proposed new wall model is applied using the LES data at $y=\varDelta$. A ‘$2\varDelta$ spatial filter’, like that used in Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005), is applied to the LES velocity at $y=\varDelta$ which is provided as the velocity input to (2.19) to model the turbulent stress at that position. This is primarily done to reduce log-layer mismatch (Yang et al. Reference Yang, Park and Moin2017) without causing an excessively sluggish response in the wall stress which would occur if the velocity was time filtered instead. The pressure gradient, on the other hand, is not spatially filtered but instead is temporally filtered with the single-sided exponential filter with the decomposition and filtering time scales described in § 4.1.

Table 1. Wall-model height in inner units $\varDelta ^+$ for all friction Reynolds numbers simulated; $Re_\tau = 3170$ corresponds to the steady-state friction Reynolds number long after the application of the spanwise pressure gradient presented in § 5.2.

Additional simulations are carried out using the traditional equilibrium wall model (without pressure gradient effects) in order to separate wall modelling dependencies from other dependencies such as grid resolution, subgrid-scale (SGS) modelling or the numerical discretizations used in the code. The EQWM used here computes the wall stress using the fitting function $Re_{\tau \varDelta }^{fit}(Re_\varDelta )$ from Meneveau (Reference Meneveau2020) which is also summarized in algorithm 1 presented in Appendix C. Then, the wall-stress vector is computed using $\boldsymbol {\tau }_w = (\nu \varDelta ^{-1} Re_{\tau \varDelta }^{fit})^2 \hat {\boldsymbol {e}}_u$ with $\hat {\boldsymbol e}_u={\boldsymbol U}_{LES}/|{\boldsymbol U}_{LES}|$.

First, we show mean velocity profiles for the five Reynolds numbers tested (figure 4) and mean Reynolds stresses for $Re_\tau =1000$ and $Re_\tau =5200$ (figure 5). For the two lowest Reynolds numbers, the WMLES using the new wall model, with combined LaRTE and laminar non-equilibrium parts, is compared with the DNS of Lee & Moser (Reference Lee and Moser2015) and the WMLES using the EQWM. The EQWM results are nearly indistinguishable from the new wall-model results. All velocity profiles follow the expected law of the wall but with a slight log-layer undershoot for $Re_\tau = 1000$ and an overshoot of the profile in the wake region at the centre of the channel. The LES Reynolds stresses generally follow the same trends as the DNS but with an overpredicted spanwise variance and underpredicted streamwise variance in the near wall region. Similar trends have been obtained in WMLES using different codes and SGS models (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2015), and are likely attributable to simulation details (i.e. grid resolution or SGS modelling) other than the wall model. The agreement between the two WMLES with completely different wall models further supports this claim. Wang, Hu & Zheng (Reference Wang, Hu and Zheng2020) showed the effect of different choices in SGS modelling, wall modelling and grid resolution on various turbulence statistics. Their primary finding is that various one-point statistics in the outer region are not significantly affected by the wall model, but are sensitive to the SGS model. As shown in the recent wall-model-independent analysis by Lozano-Durán & Bae (Reference Lozano-Durán and Bae2019), LES accuracy in the outer region of wall-bounded flows is highly sensitive to details of the ratio of grid resolution compared with the outer length scale (rather than Reynolds number). We therefore conclude that the wall model is unlikely to be the cause for the observed level of differences between the LES and DNS statistics and find that the new wall model, when applied to a standard equilibrium channel flow, performs similarly well as the classic equilibrium wall model.

Figure 4. Open circles: mean velocity profiles from WMLES using the LaRTE and non-equilibrium wall model for $Re_\tau =1000$ (red), 5200 (blue), 20 000 (magenta), $10^5$ (green) and $5\times 10^5$ (black). Plus signs: mean velocity profiles for WMLES using the EQWM at $Re_\tau = 1000$ (red) and 5200 (blue). Lines: DNS from Lee & Moser (Reference Lee and Moser2015) at $Re_\tau = 1000$ (red line), 5200 (blue line) and log law $\langle u \rangle ^+=\ln (y^+)/0.4+5.0$ (dashed line).

Figure 5. Mean Reynolds stresses for (a) $Re_\tau =1000$ and (b) $Re_\tau =5200$. Colours correspond to $\langle u'u' \rangle$ (blue), $\langle v'v' \rangle$ (red), $\langle w'w' \rangle$ (green) and $\langle u'v' \rangle$ (black). Open circles: WMLES using the new wall model with LaRTE and laminar non-equilibrium parts. Plus signs: WMLES using the EQWM. Lines: DNS from Lee & Moser (Reference Lee and Moser2015).

Next, we illustrate by means of time signals at a representative point on the wall the various terms in (2.14) that are being solved at each point following the implementation described in § 4.2. In figure 6 we show, in cyan, signals of the input stress vector $\bar {\boldsymbol {\tau }}_\varDelta$ at $y=\varDelta$, evaluated using the fitted equilibrium model in (2.19). The stress is divided by $u_\tau$, the magnitude of the obtained friction velocity. The grey line shows the same, but with the horizontal pressure gradient added, the quantity towards which the friction–velocity vector ${\boldsymbol u}_\tau$ relaxes, with relaxation time scale $T_s$. The cyan and grey lines are generally close, showing that the effect of the pressure gradient is smaller than but not negligible compared with the imposed turbulent stress. The blue line in figure 6 shows the friction velocity resulting from the LaRTE solution. In this flow, the characteristic mean value of $T_s$ can be estimated as $T_s \langle u_\tau \rangle / h = (\varDelta /h) f(\varDelta ^+) = (1/30) f(1000/30) \approx 0.45$. As is evident, major fluctuations of ${\boldsymbol u}_\tau$ occurring at time scales smaller than $T_s$ have been filtered out almost entirely. Only low-frequency variability is left, internally consistent with the notion of quasi-equilibrium that underlies the assumption of the profile scaling in inner units. Note that if an equilibrium model were used the wall stress would fluctuate at levels comparable to the cyan signal.

Figure 6. Time signals of relevant terms in the LaRTE model at some arbitrary representative point at the wall from LES or channel flow at $Re_\tau =1000$ and $\varDelta /h=1/30$. Time signals shown are for the terms $\bar {\boldsymbol \tau }_\varDelta /u_\tau$ (cyan), $(-\varDelta {\boldsymbol \nabla }_h \bar {p}/\rho + \bar {\boldsymbol \tau }_\varDelta )/u_\tau$ (grey) and ${\boldsymbol u}_\tau$ (blue). (a) Shows the $x$-component and (b) the $z$-component terms. The vertical dashed line shows the relaxation time scale $T_s \langle u_\tau \rangle /h \approx 0.45$.

Next, signals of the individual terms in (2.14) are presented in figure 7. The Eulerian time derivative shown in black displays some anticorrelated trend with the advective term shown in red. This is expected for transported quantities, and leads to smaller magnitudes of the Lagrangian time derivative as compared with the Eulerian time derivative. The blue line shows the entire relaxation towards equilibrium term which essentially drives the rate of change of the friction–velocity vector. The non-standard term with the Eulerian time derivative of the orientation vector ${\boldsymbol s}$ is negligible in the $x$-direction while it shows some contribution in the spanwise direction.

Figure 7. Time signals of terms in the evolution equation for ${\boldsymbol u}_\tau$, (2.14) at some arbitrary representative point at the wall from LES or channel flow at $Re_\tau =1000$ and $\varDelta /h=1/30$. Time signals shown are for the terms $\partial \boldsymbol {u}_\tau /\partial t$ (black), ${\boldsymbol V}_\tau \boldsymbol {\cdot } {\boldsymbol \nabla }_h {\boldsymbol u}_\tau$ (red), $-T_s^{-1}[u_\tau ^{-1}(-\varDelta {\boldsymbol \nabla }_h \bar {p}/\rho + \bar {\boldsymbol \tau }_\varDelta )-{\boldsymbol u}_\tau ]$ (blue) and $-u_\tau (\delta ^*_\varDelta /\varDelta ) \partial \bar {\boldsymbol s}/\partial t$ (green).

Wall-stress contours for $Re_\tau = 1000$ are presented in figure 8 for a single snapshot of one of the LES realizations. As can be seen from (a,b), the LaRTE quasi-equilibrium stress shows elongated structures that extend over relatively long distances downstream. The fluctuations occur, as expected around a value of $\bar {\tau }_{wx} \approx 1$. The spanwise stress component $\bar {\tau }_{wz}$ has zero mean and fluctuations that appear to occur at smaller scales, generally consistent with elongated structures that display larger variability in the transverse direction than in the streamwise direction. Panels (c,d) show the contribution from the laminar non-equilibrium portion of the model. In spite of the backward time integration that should smooth signals to some degree, these fields display much smaller-scale fluctuations. These reflect fluctuations in pressure gradients in both streamwise and spanwise directions that tend to occur at scales similar to the LES grid scale. Panels (e,f) show contours of the sum of both contributions, combining the streamwise elongated structure and the smaller-scale fluctuations from the laminar non-equilibrium part of the model.

Figure 8. Snapshots of the wall stress for $Re_\tau =1000$. (a,b) Show the quasi-equilibrium stress from the LaRTE model for both streamwise and spanwise components, (c,d) show the laminar layer non-equilibrium portion and (e,f) show the total stress ($\boldsymbol {\tau }_w = \bar {\boldsymbol {\tau }}_w+\boldsymbol {\tau }_w''$).

It is of interest to explore further the qualitative differences between an Eulerian and a Lagrangian time derivative in applying the LaRTE model. To this effect, we select some time during the LES using the LaRTE approach and denote that time as $t=0$. Then, we continue the LES using the Lagrangian version of LaRTE and perform another simulation that continues using the Eulerian version, i.e. simply omitting the advective derivative ${\boldsymbol V}_\tau \boldsymbol {\cdot } {\boldsymbol \nabla }_h {\boldsymbol u}_\tau$ from the evolution equation. Figure 9 shows the results in the form of contour plots of the $x$-component of the modelled wall stress, $\bar {\tau }_{wx}=u_\tau u_{\tau x}$. By construction, they both agree at $t=0$ but begin to differ at later times, significantly. As confirmed by examining animations, the Eulerian version ‘pins’ fluctuations at the wall while perturbations from imposed stress at $y=\varDelta$ travel downstream. The time filtering implicit in the relaxation equation then ‘smears’ and elongates the structures excessively in the streamwise direction. In the Lagrangian version shown to the left, perturbations are allowed to travel downstream, including the time-filtered versions that therefore maintain their more compact integrity as time progresses. We conclude that the Lagrangian version appears more physically reasonable. We remark, however, that we have no ‘true’ distribution (e.g. from DNS) to compare with, since we would need to evaluate either Eulerian or Lagrangian time averaging from the DNS, and similar differences would be obtained, without necessarily indicating which one is ‘better’ or ‘true’. Having seen significant differences in predicted stress distribution between Eulerian and Lagrangian versions of the model, and the latter being directly motivated by the underlying integral momentum equation, we continue using the Lagrangian version for the rest of this paper.

Figure 9. Streamwise wall stress, $\bar {\tau }_{wx}$, contours for several time instances for $Re_\tau = 1000$. Comparing Lagrangian relaxation towards equilibrium (RTE) (a,c,e) with Eulerian RTE (b,d,f) where the advection term is excluded. Both models are initialized with the same data as shown in the top row. Time is non-dimensionalized with $\langle V_{\tau x} \rangle \approx 7.93 \langle u_\tau \rangle$ and $L_x=8{\rm \pi} h$.

More quantitative characterization of the stress fluctuations is provided by the probability density function (PDF) of each component of the wall stress. The PDFs obtained from $Re_\tau =1000$ and 5200 could be compared with filtered DNS data at the same Reynolds numbers. LES data were collected over five separate uncorrelated simulations to obtain better convergence of statistics. DNS data were obtained from a public database (JHTDB 2021) and the instantaneous local wall stress was spatially filtered horizontally using a Gaussian filter at the same scale as the LES grid. The PDFs are shown in figure 10. The PDF from the filtered DNS (dashed line) peaks around $\bar {\tau }_{wx}=1$ and $\bar {\tau }_{wz}=0$. The quasi-equilibrium (LaRTE) part of the model (blue lines) peaks at the same expected values, but display significantly narrower distributions owing to the time filtering that reduces the fluctuations consistent with the notion of quasi-equilibrium. The laminar Stokes layer model that only models the fast laminar response in the viscous sublayer provides additional fluctuations. However, for the streamwise directions, these fluctuations are of smaller magnitude than those for the filtered DNS. This shows that the model is still missing significant parts of the streamwise stress fluctuations. Additional modelling is likely needed to account for these additional fluctuations that belong neither to the quasi-equilibrium nor the rapid laminar sublayer response parts of the dynamics. We note that in the spanwise direction, the PDFs agree better, in fact slightly overestimating the fluctuations for the $Re_\tau = 1000$ case but predicting the spanwise fluctuations PDF for the $Re_\tau = 5200$ case very well. From figure 10 we can also see that as the Reynolds number increases, the PDFs of the non-equilibrium components (red curves) narrow. As can be expected from (3.5), that shows the laminar non-equilibrium portion of the stress to be proportional to $\nu ^{1/2}$, the stress contribution from the laminar Stokes layer near the wall in fact vanishes in the limit of infinite Reynolds number, unlike fluctuations expected to occur due to turbulence in the wall layer. These contributions are not included in the current model and must await further developments outside of the scope of the present paper.

Figure 10. PDFs for $\tau _{wx}$ and $\tau _{wz}$; (a,b) $Re_\tau = 1000$; (c,d) $Re_\tau = 5200$. The PDF curves correspond to the filtered DNS (dashed black), the LaRTE model (blue), the non-equilibrium model (red) and the composite model (LaRTE $+$ non-equilibrium) (dashed black). DNS data obtained from the Johns Hopkins Turbulence Database JHTDB (2021) and Graham et al. (Reference Graham2016). The DNS PDFs are obtained from the Gaussian-filtered wall stress where the filtering size is the same as the LES mesh size in the horizontal directions.

We also note that when using the single-sided exponential filter with a fluctuating and short filtering time scale such as $t_\nu$, some undesirable trends can be obtained such as that the mean (in space or time) of a variable may not be exactly equal to the mean of the filtered variable. Because of this last detail, the PDF of the non-equilibrium model has a non-zero mean as seen in figure 10. This should be kept in mind whenever using the temporal exponential filter with a time-dependent averaging time scale.

5.2. Channel flow with SSPG

Next we discuss a highly non-equilibrium test case following the work of Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020). A large spanwise pressure gradient, $\partial p_\infty /\partial z$, is applied to a statistically steady turbulent channel base flow after $t=0$. Particularly, we follow the case presented in their wall modelling results section in which the initial flow ($t=0$) is standard channel flow with $Re_\tau =$1,000 after which ($t>0$) a spanwise pressure gradient is suddenly applied with strength $\partial p_\infty /\partial z=10 \partial p_\infty /\partial x = 10 \rho u_{\tau 0}^2/h$ (where $u_{\tau 0}$ is the mean friction velocity of the initial condition and $h$ the channel half-height). The flow is initialized with the results from § 5.1. The results presented in this section use the same code with the same mesh, subgrid scale and wall model, etc. We should note that dynamic time stepping is used in order to maintain a constant Courant–Friedrichs–Lewy (CFL). The time step size stays within the range $1\times 10^{-4} \leq \delta t u_{\tau 0}/h \leq 4 \times 10^{-4}$ from steady state to long after the application of the SSPG.

First, in figure 11(a,b) we show pressure gradient signals at an arbitrary point corresponding to the LES pressure gradient input $\partial \tilde {p}/\partial z$ (black line), and its three constituent parts consistent with the discussion of § 4.1: the long-time average pressure gradient $\partial P/\partial z$ (green line) entering into the full equilibrium fitted model, the band-pass-filtered fluctuating pressure gradient $\partial \bar {p}'/\partial z$ (blue line) that enters the quasi-equilibrium LaRTE equation and the rapid non-equilibrium $\partial \tilde {p}''/\partial z$ (red line) that affects mostly the viscous sublayer if sufficiently fast. As is evident in figure 11(a), $\partial \tilde {p}''/\partial z$ captures the majority of the LES pressure gradient fluctuations, $\partial P/\partial z$ captures only the ‘equilibrium’ or very slowly varying pressure gradient, and $\partial \bar {p}'/\partial z$ captures any remaining fluctuations. Figure 11(b) shows more clearly that at the onset of the SSPG ($t=0$) the equilibrium pressure gradient slowly relaxes to its new steady-state value and that the strength of the quasi-equilibrium pressure gradient fluctuations grows. Both of these pressure gradient signals are inputs to the LaRTE model whose wall stress and relevant relaxation terms are shown in figure 11(c). Here, we can see the importance of the quasi-equilibrium pressure gradient in the LaRTE model grows upon application of the SSPG.

Figure 11. Time signals at some arbitrary representative horizontal point for (a,b) spanwise pressure gradient components and (c) spanwise quantities in the LaRTE model relevant for $\bar {\tau }_{wz}$. (a,b) LES pressure gradient $\partial \tilde {p}/\partial z$ (black), non-equilibrium pressure gradient $\partial p'' /\partial z$ (red), band-pass-filtered pressure gradient $\partial \bar {p}'/\partial z$ (blue) and equilibrium pressure gradient $\partial P/\partial z$ (green) all normalized with $h/(\rho u_{\tau 0}^2)$. (c) Quasi-equilibrium spanwise wall stress $\bar {\tau }_{wz}$ (blue), $-\varDelta (\partial \bar {p}/\partial z)/\rho + \bar {\tau }_{{\rm \Delta} z}$ (grey) and $\bar {\tau }_{{\rm \Delta} z}$ (cyan) all normalized with $u_{\tau 0}^2$. Thin dashed horizontal lines indicate steady-state values before and after the SSPG.

Next we present the plane-averaged wall-stress response to the SSPG. Figure 12 shows the spanwise wall stress after the SSPG has been applied compared with the DNS of Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020). Panel (a) shows the wall-stress decomposition after the initial transient and (b) shows the wall-stress behaviour long after the SSPG was applied. The trends are in agreement with expectations. For a brief time ($0 \leq t u_{\tau 0}/h \leq 0.05$) the wall stress follows the laminar solution closely, during which the non-equilibrium component is dominant compared with the quasi-equilibrium component. Afterwards the balance is reversed. Note that the LaRTE model responds faster than the relaxation time scale which is $T_s u_{\tau 0}/h \approx 0.45$. This is due to the inclusion of the band-pass-filtered pressure gradient, $\boldsymbol {\nabla }_h \bar {p}'$, in the LaRTE model. Without this pressure gradient, $\bar {\tau }_{wz}$ is delayed by a time of order $T_s$. On the contrary, if no high-pass filtering is done, $\boldsymbol {\nabla }_h p''= 0$, $\bar {\tau }_{wz}$ is nearly linear initially and unable to capture the $\sqrt {t}$ trend corresponding to the laminar Stokes layer. Therefore, low-pass filtering is needed to prevent the overly sluggish behaviour of the quasi-equilibrium model and high-pass filtering plus the inclusion of the laminar non-equilibrium model is needed to get the correct $\sqrt {t}$ behaviour initially. Finally, note the new wall model has a closer agreement to the DNS than the EQWM which gives a linear response to the SSPG. We can attribute the faster wall-stress response to the laminar non-equilibrium model whereas the quasi-equilibrium model has a slow initial response which approaches the EQWM curve long after the SSPG. We should also note that the EQWM performance shown here is closer to the DNS than that reported in Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) for their implementation of the EQWM. We have verified that the difference is because the wall-model height used here is smaller than that of Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) which leads to a faster wall-stress response.

Figure 12. Spanwise wall stress after SSPG (a) after a short period and (b) after a long period. (a) DNS from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) (thick dashed black), composite wall stress $\langle \bar {\tau }_{wz} \rangle + \langle \tau _{wz}'' \rangle$ (solid black), quasi-equilibrium wall stress $\langle \bar {\tau }_{wz} \rangle$ (blue), non-equilibrium wall stress $\langle \tau _{wz}'' \rangle$ (red), EQWM (grey) and laminar solution for Stokes’ first problem (thin dashed black). Dotted line in (b) indicates the steady-state value after the SSPG. Angled brackets indicate ensemble averaging over the horizontal plane and five separate simulations for (a) and ensemble averaging over the horizontal plane for (b). Here, $Re_{\tau 0}=1000$ and $\partial p_\infty /\partial z = 10 \partial p_\infty /\partial x$ for $t>0$.

The streamwise wall-stress response to the SSPG is shown in figure 13. As can be seen, the wall model is unable to capture the slight initial decrease in $\tau _{wx}$ due to a complex three-dimensional mechanism discussed in Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020). The reason is that the scrambling of momentum transporting turbulent structures due to the sudden spanwise pressure gradient is not included in any part of the present model. Significantly more sophisticated modelling of the eddy viscosity in the RANS model used to derive the LaRTE equation would be required. In this case, however, the difference is less than 1 %–2 % of $u_{\tau 0}$. The model correctly relaxes towards the DNS trend for $tu_{\tau 0}/h>1$. The increase in $\tau _{wx}$ may be attributed to the increase in Reynolds number as the mean pressure gradient increases in magnitude even as its alignment rotates away from the $x$-axis. This behaviour is ‘slow’ and thus is expected to be captured well by the LaRTE model. This also means that the EQWM is able to capture this behaviour well, as evident from the similarity between the equilibrium and quasi-equilibrium curves in figure 13. After a long time, the new equilibrium condition is reached in which the $x$-component of the pressure gradient must be balanced by the wall stress and thus the wall stress reduces back to unity as shown in figure 13(b).

Figure 13. Streamwise wall stress after SSPG (a) after a short period and (b) after a long period. (a) DNS from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) (thick dashed black), composite wall stress $\langle \bar {\tau }_{wx} \rangle + \langle \tau _{wx}'' \rangle$ (solid black) and EQWM (grey). Dotted line in (b) indicates the steady-state value after the SSPG. Angled brackets indicate ensemble averaging over the horizontal plane and five separate simulations for (a) and ensemble averaging over the horizontal plane for (b). Here, $Re_{\tau 0}=1000$ and $\partial p_\infty /\partial z = 10 \partial p_\infty /\partial x$ for $t>0$.

Figure 14 shows mean spanwise velocity profiles after the SSPG for several different time instances. Generally, the LES agrees with the DNS for both wall models quite well. This is consistent with the results reported in Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) where all wall models considered produced good spanwise velocity profiles upon application of the SSPG. This is also consistent with the notion that one-point statistics are less sensitive to the wall model relative to other simulation parameters, as argued earlier in this section.

Figure 14. Mean spanwise velocity profiles at $t u_{\tau 0}/h = 0.21$, 0.405, 0.615, 0.825 and 1.02 from bottom to top. Open circles: WMLES with composite LaRTE and laminar non-equilibrium components. Plus signs: WMLES with the EQWM. Lines: DNS from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020).

Figure 15 shows contours of fluctuations of the quasi-equilibrium stress $\bar {\boldsymbol \tau }_w$. Specifically, we show contours of $\bar {\tau }_{ws}'=\bar {\boldsymbol \tau }_w \boldsymbol {\cdot } \langle {\boldsymbol s} \rangle - \langle \bar {\boldsymbol \tau }_w \rangle \boldsymbol{\cdot} \langle {\boldsymbol s} \rangle$, where $\langle {\boldsymbol s} \rangle$ is the plane-averaged unit vector, i.e. in the direction of the mean LaRTE wall stress. The contours represent the wall-stress fluctuations aligned with the plane-averaged mean quasi-equilibrium wall-stress direction. As a reference, the dashed lines shown are aligned with the total wall stress and thus include the contributions from the laminar boundary layer developing due to the application of the SSPG. The application of the SSPG disrupts the orientation and shape of the structures as the mean flow rotates towards the $z$ direction. At later times, after the SSPG is applied, panels $(b,d,f)$ show that the structures have had enough time to orient and advect themselves with the mean wall-stress direction, albeit with a reduced size.

Figure 15. Contours of the quasi-equilibrium $s$-component wall stress ($\boldsymbol {s}$ introduced in § 2) with the plane-averaged mean subtracted, $\bar {\tau }_{ws}'$, for various times after the SSPG. (a,c,e) Immediately after SSPG; (b,d,f) later times after SSPG. Dashed lines are aligned with the plane-averaged total wall-stress angle (includes both quasi-equilibrium and non-equilibrium components).

Figure 16 shows the evolution of the angle (with respect to the $x$ axis) of the pressure gradients (dashed lines) and resulting plane-averaged total (black line) and quasi-equilibrium (blue line) wall stresses. As expected, the quasi-equilibrium component of the wall stress has a delayed and smoothed response to the SSPG, while at the initial transient the non-equilibrium component dominates the plane-averaged wall-stress response. After some adjustment time the quasi-equilibrium component becomes more dominant, again in establishing the trends of the plane-averaged wall-stress components such as its direction. At large times, the angle tends to the same angle as the net applied pressure gradient.

Figure 16. Plane-averaged pressure gradient and wall-stress angles after the SSPG; $\langle \boldsymbol {\nabla }_h \bar {p} \rangle$ (dashed blue), $\langle \boldsymbol {\nabla }_h p_\infty \rangle$ (dashed black), $\langle \bar {\boldsymbol {\tau }}_w \rangle$ (solid blue), $\langle \bar {\boldsymbol {\tau }}_w \rangle + \langle \boldsymbol {\tau }_w'' \rangle$ (solid black) where angled brackets indicate plane averaging.