A Lagrangian relaxation towards equilibrium wall model for large eddy simulation

Abstract

A large eddy simulation wall model is developed based on a formal interpretation of quasi-equilibrium that governs the momentum balance integrated in the wall-normal direction. The model substitutes the law-of-the-wall velocity profile for a smooth surface into the wall-normal integrated momentum balance, leading to a Lagrangian relaxation towards equilibrium (LaRTE) transport equation for the friction–velocity vector ${\boldsymbol u}_\tau (x,z,t)$. This partial differential equation includes a relaxation time scale governing the rate at which the wall stress can respond to imposed fluctuations due to the inertia of the fluid layer from the wall to the wall-model height. A priori tests based on channel flow direct numerical simulation (DNS) data show that the identified relaxation time scale ensures self-consistency with assumed quasi-equilibrium conditions. The new approach enables us to formally distinguish quasi-equilibrium from additional, non-equilibrium contributions to the wall stress. A particular model for non-equilibrium contributions is derived, motivated by laminar Stokes layer dynamics in the viscous sublayer when applying fast-varying pressure gradients. The new wall modelling approach is first tested in standard equilibrium channel flow in order to document various properties of the approach. The model is then applied in large eddy simulation of channel flow with a suddenly applied spanwise pressure gradient (SSPG). The resulting mean wall-stress evolution is compared with DNS with good agreement. At the onset of the SSPG, the laminar Stokes layer develops rapidly while the LaRTE portion of the stress has a delayed response due to its inherent relaxation dynamics. Results also highlight open challenges such as modelling the response of near-wall turbulence occurring above the viscous sublayer and at time scales faster than quasi-equilibrium conditions.

1. Introduction

Large eddy simulation (LES) has become a popular prediction tool for unsteady turbulent flows due to its ability to resolve the less universal large-scale motions while modelling the more universal small scales. However, near the wall, the momentum-carrying structures scale with the viscous length. Thus, LES that resolve these structures near the wall (wall-resolved LES) incur computational costs that are not significantly lower than those of direct numerical simulation (DNS). For recent detailed analyses of computational costs of DNS and wall-resolved LES, see Choi & Moin (2012) and Yang & Griffin (2021). Wall-resolved LES becomes computationally intractable for high Reynolds number flows relevant for many engineering and geophysical applications. Therefore, in order for LES to be a practical flow prediction tool, the near-wall region must be modelled through wall models.

Numerous wall models have been proposed over the years (see Piomelli & Balaras (2002), Piomelli (2008), Larsson et al. (2016) and Bose & Park (2018), for reviews of LES wall modelling). Larsson et al. (2016) broadly identify different classes for wall models based on the coupling between the resolved and modelled regions. The wall models used and discussed in the present work are considered to be ‘wall-stress models’ in which the LES formally extends all the way to the wall and the wall stress is applied as a boundary condition to the LES. The equilibrium wall model (EQWM) is the simplest and most commonly used wall-stress model. The algebraic EQWM assumes that the velocity profile follows some known functional form. Typically, one assumes the velocity satisfies the ‘law of the wall’ such that $\langle u \rangle ^+ = f(y^+)$, where $y$ is the wall-normal coordinate and ‘$+$’ indicates inner units non-dimensionalization with the friction velocity $u_\tau$ and the kinematic viscosity $\nu$. If the wall-model height, $y_{wm}=\varDelta$, lies within the log layer then one assumes the log law. For rough walls, this is the approach most often followed in the geophysical literature (Moeng 1984; Bou-Zeid, Meneveau & Parlange 2005). Inversion of a composite profile including the viscous and buffer layers is more general for smooth wall applications, see e.g. Luchini (2018), Gonzalez, Adler & Gaitonde (2018) and Adler et al. (2020). Algebraic EQWMs usually assume the velocity profile to be valid locally and instantaneously such that the LES velocity at the wall-model height may be used to find the local friction velocity and thus the local wall stress. The log law for smooth walls is a nonlinear equation for the friction velocity and thus must be solved numerically, typically using iterations. As a more practical alternative, Meneveau (2020) developed fitting functions that directly compute the friction velocity given the velocity at the wall-model height, also including pressure gradient and roughness effects.

Differential forms of the EQWM are often also used through simplified forms of the thin boundary-layer equations. Over the past decade, the name ‘equilibrium wall model’ often refers to the numerical solution of the wall-normal diffusion equation $\partial _y [(\nu + \nu _T)\partial _y \langle u \rangle ] = 0$, where the functional form for the eddy viscosity, $\nu _T$, is assumed to be known, typically a mixing length model. Algebraic EQWMs are explicit solutions or approximations to the numerically solved differential EQWMs (Meneveau 2020).

Most wall models are based on the thin boundary-layer equations (TBLE) which may be considered Reynolds-averaged Navier–Stokes (RANS)-like in nature since the local momentum-carrying scales of motion are similar to or smaller than the computational LES grid. So-called ‘zonal’ methods, a type of hybrid method discussed in Larsson et al. (2016), solve these equations on a separate RANS mesh, refined in the wall-normal direction, below a user-defined RANS/LES interface height. The advantage of zonal methods is their ability to include all terms in the TBLE (such as unsteady, acceleration and pressure gradient terms), thereby allowing for a greater range of applicability. However, zonal methods have increased cost due to the wall-normal grid refinement which ultimately can lead to costs approaching wall-resolved LES. Therefore, there is continued interest towards methods that do not require a separate RANS mesh. The dynamic slip wall model (Bose & Moin 2014; Bae et al. 2019) is an alternative type of wall model that does not solve RANS-like equations or make assumptions about the flow (e.g. the law of the wall). The approach has the benefit of making connections to the dynamic model (Germano et al. 1991) using resolved flow information near the wall. However, at increasing Reynolds numbers, it is unclear that such modelling can accurately capture subtle Reynolds number dependencies of the predicted friction drag because information regarding the viscous sublayer is unavailable in such approaches.

Chung & Pullin (2009) derived a wall model based on vertically integrating the unsteady term in the TBLE. A ‘law-of-the-wall’ assumption for the velocity profile is invoked such that $\langle u \rangle = u_\tau f(y^+)$ and thus, upon differentiation with time, the chain rule gives a $\partial u_\tau /\partial t$ term leading to an ordinary differential equation (ODE) in time for the wall stress. This allows us to obtain an evolution equation for the wall stress as a function of known LES quantities at the model height, where advection terms in the TBLEs are approximated by the resolved LES terms, i.e. assuming plug flow profile below the wall-model height. Inspired by this work but allowing for deviations from plug flow, the integral wall model (iWMLES) introduced by Yang et al. (2015) aims to achieve an algebraic closure based on an approach similar to the von Kármán Pohlhausen method. The assumed velocity profile is linear in the viscous sublayer up to the buffer layer and then switches to a log profile plus a linear correction to represent possible pressure gradient effects. A series of coefficients are determined using matching conditions consisting of the boundary conditions and the full vertically integrated TBLEs. Since iWMLES is based on the full TBLEs it can in principle capture deviations from full equilibrium conditions while still maintaining affordability due to its algebraic nature. However, while simpler than solving additional ODEs on finer meshes, the iWMLES approach still involves solving several coupled nonlinear equations, requires specifying an empirically chosen time scale for exponential time filtering of the velocity input at the wall-model height and involves a discontinuity in the slope of the assumed velocity profile in the buffer region.

Additionally, in practice the so-called log-layer mismatch is a known problem in wall modelling. Kawai & Larsson (2012) argue that the near-wall region of wall-modelled LES (WMLES) is inherently under-resolved numerically, so the wall-model height should ideally be chosen further away from the wall in order to reduce the log-layer mismatch. Bou-Zeid et al. (2005) showed through the Schwartz inequality that a local law of the wall formulation inherently over-predicts the wall stress due to LES velocity fluctuations. They then proposed reducing this error by filtering the velocity using a $2\varDelta$ spatial filter. Similarly, Yang, Park & Moin (2017) suggest local temporal filtering and/or wall-parallel spatial filtering to reduce log-layer mismatch. Recently, Hosseinzade & Bergstrom (2021) tested various exponential filtering time scales of the input velocity and pressure gradient in the context of solving unsteady RANS equations with a finely resolved wall-normal grid. Unsteady and horizontal advection terms were also included. Moreover, and consistent with the recommendations by Kawai & Larsson (2012), Hosseinzade & Bergstrom (2021) found that placing the wall-model grid point inside the LES region (they used the fifth LES grid point) improved results. They also find that the use of time filtering is more important when placing the wall-model height at the first grid point as opposed to inside the LES region.

In the present work we aim to extend the iWMLES and Chung & Pullin (2009) approaches to arrive at a unified formulation that is both based on formal derivations starting from the underlying governing equations and is robust for applications. It will be shown that the model takes the form of a Lagrangian relaxation, hence it is termed the Lagrangian relaxation towards equilibrium (LaRTE) model. As will be seen, a number of arbitrary choices that had to be made in the context of iWMLES, such as specifying an empirically chosen exponential filtering time scale (also needed in the approach by Hosseinzade & Bergstrom 2021), will no longer be required. The derivation of the LaRTE wall model is presented in § 2 and is based on a self-consistent interpretation of ‘quasi-equilibrium’ in the assumed velocity profile below the wall-model height. This new interpretation then enables us to formally distinguish quasi-equilibrium from additional, non-equilibrium contributions to the wall stress and the latter can then be modelled separately.

A second part of this work then concerns a model for the non-equilibrium flow and stress response in the viscous sublayer to temporally changing pressure gradients. For applications involving time-varying applied pressure gradients, it is important to supply the new quasi-equilibrium LaRTE model with additional components that reflect the deviations from the assumed simple velocity profile in the near-wall region. Numerous prior efforts have been made to understand unsteady effects on wall-bounded turbulent flows (e.g. the works by Jung, Mangiavacchi & Akhavan (1992), Coleman, Kim & Le (1996), Scotti & Piomelli (2001), He & Seddighi (2015), Weng, Boij & Hanifi (2016), Jung & Kim (2017), Sundstrom & Cervantes (2018c) and Lozano-Durán et al. (2020) and more detailed background is provided in § 3). These studies, in which a laminar Stokes layer is observed near the wall, motivate the non-equilibrium model introduced in § 3. These studies also discuss how turbulence structure itself is affected by the unsteadiness above the viscous sublayer, however, we leave modelling of such effects for future work.

After introducing the two basic ingredients for the new wall model, in § 4 we address practical implementation issues. Applications to equilibrium and non-equilibrium channel flows are described in § 5. Various features of the model are presented via time series at individual points as well as contour plots of modelled wall-stress components. The tests in non-equilibrium flow are performed for channel flow upon which a very strong sudden spanwise pressure gradient (SSPG) is applied. Comparisons with DNS results by Lozano-Durán et al. (2020) are also included. Summary and conclusions are presented in § 6.

2. LaRTE wall model

Following the ideas underlying the iWMLES (Yang et al. 2015) we assume that between the wall and the wall-model height an LES grid point at a distance $\varDelta$ away from the wall there exists a quasi-equilibrium mean velocity profile (see figure 1a). An overline denotes the corresponding averaging operation, which may be interpreted as a horizontally grid-filtered quantity at the LES scale in the wall-parallel plane, and additional (implicit) temporal averaging whose properties will become apparent from the derivation itself. The key assumption underlying the proposed wall model is that in the horizontal (wall parallel, $x$–$z$) plane, the mean velocity $\bar {\boldsymbol u}_s = \bar {u} \hat {\boldsymbol \imath } + \bar {w} \hat {\boldsymbol k}$ ($\hat {\boldsymbol \imath }$ and $\hat {\boldsymbol k}$ are the two unit vectors on the wall) can be written according to

(2.1)\begin{equation} \bar{\boldsymbol u}_s(x,y,z,t) = {\boldsymbol u}_{\tau}(x,z,t) f(y^+), \end{equation}

where ${\boldsymbol u}_\tau (x,z,t)$ is the friction–velocity vector and is a slowly varying function of the horizontal positions $x,z$ and time $t$. The characteristic time scale characterizing what is termed ‘slow’ evolution is not prescribed a priori but will be shown to arise directly from the assumption of quasi-equilibrium. The inner similarity function $f(y^+)$ with $y^+ = y u_\tau /\nu$ is the assumed velocity profile in inner units, and $u_{\tau } = |{\boldsymbol u}_{\tau }|$. Typically, $f(y^+)$ includes a linear region near the wall merging with a logarithmic portion above the buffer layer but the precise shape of $f(y^+)$ is not important at initial stages of development. We remark that in the present work we deal exclusively with smooth planar walls.

Figure 1. (a) Sketch of assumed inner velocity profile (in blue) representing a quasi-equilibrium RANS solution in the inner layer, responding to an outer ‘applied’ total shear stress $\bar {\boldsymbol {\tau }}_\varDelta$ at the wall-model height at $y=\varDelta$. (b) Sketch of stresses acting on the fluid layer between $y=0$ and $y=\varDelta$, leading to inertia term with response time scale $T_s$ proportional to $\varDelta$.

The full quasi-steady velocity is then given by

(2.2)\begin{equation} \bar{\boldsymbol u} = \bar{\boldsymbol u}_s + \bar{v} \hat{\boldsymbol \jmath}, \end{equation}

where $\bar {v}$ is the wall-normal velocity and $\hat {\boldsymbol \jmath }$ the unit vector in the $y$-direction. The friction–velocity vector ${\boldsymbol u}_{\tau } = u_{\tau x}{\hat {\boldsymbol \imath }} + u_{\tau z}{\hat {\boldsymbol k}}$ is defined such that the (kinematic) wall stress vector $\bar {\boldsymbol \tau }_w$ (its two components in the wall plane) is given by

(2.3)\begin{equation} \bar{\boldsymbol\tau}_w = u_{\tau} {\boldsymbol u}_{\tau}, \end{equation}

i.e. $\bar {\boldsymbol u}_s$ and ${\boldsymbol u}_{\tau }$ are in the same direction as $\bar {\boldsymbol \tau }_w$. This direction will be represented by unit vector ${\boldsymbol s}$ (that also can depend on $x,z,t$), i.e. ${\boldsymbol u}_\tau = u_\tau {\boldsymbol s}$ and $\bar {\boldsymbol u}_s = u_\tau f(y^+) {\boldsymbol s}$ (figure 1a).

Next, we aim to derive an evolution equation for the friction velocity vector ${\boldsymbol u}_{\tau }(x,z,t)$ that is consistent with the RANS evolution for $\bar {\boldsymbol u}$

(2.4)\begin{equation} \frac{\partial \bar{\boldsymbol{u}}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \bar{\boldsymbol{u}} \bar{\boldsymbol{u}} \right) ={-} \frac{1}{\rho} \boldsymbol{\nabla}\bar{p}+ \boldsymbol{\nabla} \boldsymbol{\cdot} \left[(\nu+\nu_T) (\boldsymbol{\nabla} \bar{\boldsymbol{u}}+ \boldsymbol{\nabla} \bar{\boldsymbol{u}}^\top) \right], \end{equation}

where $\nu _T(x,y,z,t)$ is the position-dependent eddy viscosity associated with the RANS model being considered and $\bar {p}(x,z,t)$ is the quasi-equilibrium pressure with no wall-normal dependence, consistent with boundary-layer approximations. The momentum equation for the wall-parallel velocity $\bar {\boldsymbol {u}}_s$ (2 components) reads

(2.5)\begin{align} \frac{\partial \bar{\boldsymbol{u}}_s}{\partial t} + \boldsymbol{\nabla}_h \boldsymbol{\cdot} \left( \bar{\boldsymbol{u}}_s \bar{\boldsymbol{u}}_s \right) + \partial_y ( \bar{v} \bar{\boldsymbol{u}}_s ) &={-} \frac{1}{\rho} \boldsymbol{\nabla}_h\bar{p}+ \frac{\partial}{\partial y} \left[ (\nu+\nu_T) \frac{\partial \bar{\boldsymbol{u}}_s}{\partial y} \right] \nonumber\\ &\quad + \boldsymbol{\nabla}_h \boldsymbol{\cdot} \left[(\nu+\nu_T) \left( \boldsymbol{\nabla}_h \bar{\boldsymbol{u}}_s+\boldsymbol{\nabla}_h \bar{\boldsymbol{u}}_s^\top \right) \right], \end{align}

where $\boldsymbol {\nabla }_h=\partial _x \hat {\boldsymbol \imath } +\partial _z \hat {\boldsymbol k}$ represents the horizontal gradients on the $x-z$ wall plane, and diffusion cross-terms involving the (small) vertical velocity $\bar {v}$ have been neglected. Into this equation we replace the main ansatz (2.1). And, following the logic by Chung & Pullin (2009) and the iWMLES by Yang et al. (2015), we integrate from $y=0$ to the wall-model height at $y=\varDelta$. The procedure will be illustrated via the first term, the Eulerian time derivative $\partial _t \bar {\boldsymbol u}_s$, but similar steps can be applied to the advective derivative terms also invoking the continuity equation as shown in Appendix A. The key steps for the time derivative read as follows:

(2.6)\begin{align} \frac{\partial}{\partial t} \int_0^\varDelta \bar{\boldsymbol{u}}_s \,{{\rm d}y} &= \frac{\partial}{\partial t}\left[{\boldsymbol s} \int_0^\varDelta u_\tau f \left(\frac{y u_\tau}{\nu} \right){{\rm d}y}\right] \nonumber\\ &= {\boldsymbol s} \int_0^\varDelta \frac{\partial u_\tau}{\partial t} \left( f(y^+) + u_\tau f^\prime(y^+) \frac{y}{\nu} \right) {{\rm d}y} + \frac{\partial {\boldsymbol s}}{\partial t} \int_0^\varDelta u_\tau f(y^+)\,{{\rm d}y}. \end{align}

As recognized by Chung & Pullin (2009) in their derivation of an integral boundary-layer equation-based wall model, the first integral on the right-hand side can be rewritten with ${\rm d}/{{\rm d}y}^+[y^+ f(y^+)]$ as integrand, resulting in

(2.7)\begin{align} \frac{\partial}{\partial t} \int_0^\varDelta \bar{\boldsymbol{u}}_s \,{{\rm d}y} &= {\boldsymbol s} \frac{\partial u_\tau}{\partial t} \int_0^{\varDelta^+} \frac{{\rm d}}{{{\rm d}y}^+} \left[ y^+f(y^+) \right] {{\rm d}y}^+ \frac{\nu}{u_\tau} + u_\tau \frac{\partial {\boldsymbol s}}{\partial t} \int_0^\varDelta f(y^+) \, {{\rm d}y} \nonumber\\ &= {\boldsymbol s} \frac{\partial u_\tau}{\partial t} {\rm \Delta} f(\varDelta^+) + u_\tau \frac{\partial {\boldsymbol s}}{\partial t} \int_0^\varDelta f(y^+) \, {{\rm d}y} \nonumber\\ &= \frac{\partial (u_\tau{\boldsymbol s}) }{\partial t} {\rm \Delta} f(\varDelta^+) + u_\tau \frac{\partial {\boldsymbol s}}{\partial t} \left(\int_0^\varDelta [f(y^+)-f(\varDelta^+)]\, {{\rm d}y}\right). \end{align}

The last term motivates definition of a ‘cell displacement thickness’

(2.8)\begin{equation} \delta^*_\varDelta = \int_0^\varDelta \left(1-\frac{\bar{u}_s(y)}{\bar{u}_s(\varDelta)} \right){{\rm d}y} \quad\to \quad \frac{\delta^*_\varDelta}{\varDelta} =\frac{1}{\varDelta^+} \int_0^{\varDelta^+} \left(1-\frac{f(y^+)}{f(\varDelta^+)} \right) {{\rm d}y}^+, \end{equation}

analogous to the boundary-layer displacement thickness but integrated only up to $y=\varDelta$. Finally, the Eulerian time derivative term can be written according to

(2.9)\begin{equation} \frac{\partial}{\partial t} \int_0^\varDelta \bar{\boldsymbol{u}}_s \,{{\rm d}y} = {\rm \Delta} f(\varDelta^+) \frac{\partial {\boldsymbol u}_\tau}{\partial t} - u_\tau f(\varDelta^+) \delta^*_\varDelta \frac{\partial {\boldsymbol s}}{\partial t}. \end{equation}

Integration of the advective term, i.e. $\int _0^\varDelta \boldsymbol {\nabla }_h \boldsymbol {\cdot } (\bar {\boldsymbol {u}}_s \bar {\boldsymbol {u}}_s ) + \partial _y ( \bar {v} \bar {\boldsymbol {u}}_s )\, {{\rm d}y}$ requires the mean vertical velocity at $y=\varDelta$, since $\int _0^\varDelta \partial _y ( \bar {v} \bar {\boldsymbol {u}}_s) \, {{\rm d}y} = \bar {v}(\varDelta ) \bar {\boldsymbol {u}}_s(\varDelta ) {\boldsymbol s}$ (and ${\boldsymbol s}$ does not depend on $y$). Using the continuity equation $\partial _s \bar {u}_s + \partial _y \bar {v}=0$ we obtain

(2.10)\begin{align} \bar{v}(\varDelta) ={-} \int_0^\varDelta \partial_s \bar{u}_s(y)\,{{\rm d}y} &={-} \frac{\partial}{\partial s} \left[u_\tau \int_0^\varDelta f(y^+)\,{{\rm d}y}\right] \nonumber\\ &={-} \frac{\partial u_\tau}{\partial s} \int_0^\varDelta \left[f(y^+) + y^+ f'(y^+)\right] {{\rm d}y} ={-} \frac{\partial u_\tau}{\partial s} {\rm \Delta} f(\varDelta^+). \end{align}

As further shown in detail in Appendix A the entire integral of the advective term can then be written as

(2.11)\begin{equation} \int_0^\varDelta \left[\boldsymbol{\nabla}_h \boldsymbol{\cdot} \left( \bar{\boldsymbol{u}}_s \bar{\boldsymbol{u}}_s \right) + \partial_y ( \bar{v} \bar{\boldsymbol{u}}_s ) \right]{{\rm d}y} ={\rm \Delta} f(\varDelta^+)\,{\boldsymbol V}_\tau \boldsymbol{\cdot} \boldsymbol{\nabla}_h {\boldsymbol u}_\tau, \end{equation}

where

(2.12)\begin{equation} {\boldsymbol V}_\tau = \left(1-\frac{\delta^*_\varDelta}{\varDelta} - \frac{\theta_\varDelta}{\varDelta}\right) f(\varDelta^+) {\boldsymbol u}_\tau \end{equation}

is the advective velocity and a ‘cell momentum thickness’

(2.13)\begin{equation} \theta_\varDelta = \int_0^\varDelta \frac{\bar{u}_s(y)}{\bar{u}_s(\varDelta)} \left(1-\frac{\bar{u}_s(y)}{\bar{u}_s(\varDelta)} \right){{\rm d}y} \quad \to \quad \frac{\theta_\varDelta}{\varDelta} = \frac{1}{\varDelta^+} \int_0^{\varDelta^+} \frac{f(y^+)}{f(\varDelta^+)} \left(1-\frac{f(y^+)}{f(\varDelta^+)} \right) {{\rm d}y}^+, \end{equation}

has been introduced arising from the integrals of quadratic advection terms.

Integrating each term in (2.5) between $y=0$ and $y=\varDelta$, it is convenient to define the total (molecular viscous $+$ turbulent viscous) shear stress at $y=\varDelta$ according to $\bar {\boldsymbol \tau }_\varDelta = (\nu +\nu _T)\partial \bar {\boldsymbol u}_s/\partial y$. Collecting terms, replacing ${\boldsymbol u}_\tau = u_\tau {\boldsymbol s}$ and dividing the entire equation by ${\rm \Delta} f(\varDelta ^+)$, the evolution equation for the friction–velocity vector can now be written according to

(2.14)\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} + \bar{{\boldsymbol{\tau}}}_\varDelta \right) - {\boldsymbol u}_\tau\right] + u_\tau \frac{\delta^*_\varDelta}{\varDelta} \frac{\partial {\boldsymbol s}}{\partial t} + \frac{1}{{\rm \Delta} f(\varDelta^+)} {\boldsymbol \nabla}_h \boldsymbol{\cdot} \boldsymbol{\mathsf{D}}_\tau, \end{equation}

where $T_s$ is given by

(2.15)\begin{equation} T_s = f(\varDelta^+) \frac{\varDelta }{u_\tau}. \end{equation}

It represents a time scale that arises from the derivation of (2.14) and does not require additional ad hoc assumptions.

The horizontal diffusion flux tensor integrated in the vertical direction, $\boldsymbol{\mathsf{D}}_\tau$, is defined according to

(2.16)\begin{equation} \boldsymbol{\mathsf{D}}_\tau = \int_0^\varDelta (\nu+\nu_T) \left( \boldsymbol{\nabla}_h \bar{\boldsymbol{u}}_s+\boldsymbol{\nabla}_h \bar{\boldsymbol{u}}_s^\top \right) {{\rm d}y}, \end{equation}

and is further detailed in Appendix B. When coupled with appropriate models for $\bar {{\boldsymbol {\tau }}}_\varDelta$, $f(\varDelta ^+)$, $\delta ^*_\varDelta$, $\theta _\varDelta$ and $\boldsymbol{\mathsf{D}}_\tau$, we refer to (2.14) as the evolution equation underlying the LaRTE wall model.

2.1. Discussion

For the sake of initial discussion, it is instructive to consider a simplified form for the (e.g.) streamwise $x$-component of (2.14), for now neglecting the pressure gradient and diffusive terms, as well as the direction change term $\partial {\boldsymbol s}/\partial t$. Under these simplifying conditions, (2.14) can be written as

(2.17)\begin{equation} \frac{{{\rm d}}_s u_{\tau x}}{{\rm d} t} = \frac{1}{T_s} \left( \bar{{\tau}}_{{\rm \Delta} x}/u_\tau - {u}_{\tau x} \right), \end{equation}

with ${{\rm d}}_s/{\rm d} t=\partial _t + {\boldsymbol V}_\tau \boldsymbol {\cdot } {\boldsymbol \nabla }_h$ representing a Lagrangian time derivative on the surface. In this form, it becomes apparent that the model represents a Lagrangian relaxation dynamics, with $T_s$ serving as the relaxation time scale for how the friction–velocity component ${u}_{\tau x}$ approaches the stress at the wall-model grid point in LES ($\bar {{\tau }}_{{\rm \Delta} x}$) (the latter divided by the friction–velocity magnitude $u_\tau$). For the present discussion, neglecting $\bar {\tau }_{w z}$, i.e. with $\bar {\tau }_{wx}=u_\tau u_{\tau x} = u_{\tau x}^2$, we see by multiplying (2.17) by $u_{\tau x}$ that in terms of the wall stress the Lagrangian relaxation equation can equivalently be written as

(2.18)\begin{equation} \frac{{\rm d}_s \bar{\tau}_{w x}}{{\rm d} t} = \frac{2}{T_s} \left(\bar{\tau}_{{\rm \Delta} x} - \bar{\tau}_{w x}\right), \end{equation}

showing that for the stress the relaxation time scale is $T_s/2$. This time scale was originally derived by Chung & Pullin (2009) where they also showed the wall stress tends (in an Eulerian sense) towards its steady-state value at a rate corresponding to this time scale (in their work $1/(\varLambda \tilde {\eta }_0)$ is equivalent to $T_s/2$ shown here).

It can be seen from (2.9) that the relaxation time scale $T_s$ arises from integrating the assumed velocity profile between $y=0$ and $y=\varDelta$ (this integral is similar to the term $\partial L/\partial t$ that arises in the iWMLES approach by Yang et al. 2015). As a result of the analysis, $T_s$ is proportional to the volume per unit area ($\varDelta$) in the fluid layer under consideration, between the wall and the wall-model height $y=\varDelta$. It represents the inertia of the fluid in that layer and can be seen to cause a time delay between the stress at $y=\varDelta$ and at the wall, under unsteady conditions. The thicker the layer (large $\varDelta$), the more the time delay due to added fluid inertia. Conversely, the stronger the turbulence (large $u_\tau$), the faster the relaxation leading to a smaller time delay. We note that, in Yang et al. (2015) for the iWMLES approach, an explicit time filtering at a time scale $\sim \varDelta /\kappa u_\tau$ was introduced, operationally similar to but shorter than $T_s$ by a factor $\kappa f(\varDelta ^+) = \kappa \bar {u}(\varDelta )/u_\tau$. Here, such temporal relaxation behaviour has been derived formally from the momentum equation (unsteady RANS) and the assumed validity of a quasi-equilibrium velocity profile (2.1).

Also, for the full vector problem with two stress components, we remark that, when attempting to write the relaxation equation in terms of wall-parallel components $\bar {\tau }_{wx}$ and $\bar {\tau }_{wz}$ instead of friction–velocity vector components $u_{\tau x}$ and $u_{\tau z}$ (or ${\boldsymbol u}_{\tau }$), the resulting equation is far less intuitive compared with the relatively simple form of (2.14). The latter resembles a standard transported vector field equation with a relaxation source term and a diffusion term (only the $\partial {\boldsymbol s}/\partial t$ term is non-standard) and is therefore much preferable.

Finally, we note that when including the pressure gradient term $\boldsymbol {\nabla }_h \bar {p}$, (2.14) shows that the implied relaxation dynamics is how ${\boldsymbol u}_\tau$ approaches the total cell forcing (vector) term $(-\rho ^{-1}{\boldsymbol \nabla }_s \bar {p}\varDelta + \bar {{\boldsymbol \tau }}_\varDelta )/u_\tau$.

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 (2020) in which an equilibrium layer partial differential equation is numerically integrated and various fitting functions are developed for this solution. The model of Meneveau (2020) 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 2020). 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 (2020) 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. 2016). 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 (2020) 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. (2015), 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. (2015).

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.

3. Non-equilibrium laminar fast response model

So far we have developed a model for $\bar {\boldsymbol {\tau }}_w$, the quasi-equilibrium part of the wall stress, which responds to external conditions (changes in velocity and pressure gradients at the wall-model height) at a time scale $T_s$ consistent with the assumption of quasi-equilibrium. In order to supplement the quasi-equilibrium Lagrangian relaxation model with additional physics, we now focus on the rapid response of the inner-most part of the boundary layer, the viscous sublayer. The response of near-wall structures of turbulence to rapidly changing pressure gradients has been studied extensively in the past. Jung et al. (1992), Karniadakis & Choi (2003), Quadrio & Ricco (2003), Ricco et al. (2012) and Yao, Chen & Hussain (2019) studied spanwise wall oscillations due to their drag-reducing capabilities. Vardy & Brown (2003) and Vardy et al. (2015) attempt to understand the wall shear stress for water hammer pipe flows. Experimental and numerical studies of pulsatile flows have been performed by Scotti & Piomelli (2001), Tardu & da Costa (2005), Tardu & Maestri (2010), Weng et al. (2016), Sundstrom & Cervantes (2018a) and Cheng et al. (2020). Streamwise accelerating flows were considered by He & Jackson (2000), Greenblatt & Moss (2004), He, Ariyaratne & Vardy (2008), He, Ariyaratne & Vardy (2011), He & Ariyaratne (2011), He & Seddighi (2013), He & Seddighi (2015), Jung & Chung (2012), Jung & Kim (2017), Sundstrom & Cervantes (2017) and Sundstrom & Cervantes (2018c). Then there are studies with step changes in either the wall boundary condition or in the pressure forcing such as the sudden spanwise wall movement of Coleman et al. (1996), Tang & Akhavan (2016) and Abe (2020), the sudden spanwise pressure gradient of Moin et al. (1990) and Lozano-Durán et al. (2020) or a change in the direction of the pressure gradient forcing (for a variety of directions) as done in de Wiart, Larsson & Murman (2018). One of the common observations amongst all of these flows is the existence of a laminar Stokes layer near the wall. For pulsatile flows, the wall-stress deviation from its steady-state value follows the solution to Stokes’ second problem exactly for high-frequency oscillations (Weng et al. 2016). For streamwise accelerating flows, the wall stress deviation from its initial value follows the solution to Stokes’ first problem during the first stage of the acceleration (He & Seddighi 2015; Jung & Kim 2017; Sundstrom & Cervantes 2018c). Sundstrom & Cervantes (2018b) even showed that the wall stress for low-frequency pulsations follows the Stokes solution during the acceleration phase of the pulse. They further showed the similarity of the wall stress during the acceleration phase of pulsatile flows with the initial phase of streamwise accelerating flows. For the sudden spanwise wall movement and SSPG, the spanwise velocity and wall stress follows Stokes’ first problem during the early response (Coleman et al. 1996; Abe 2020; Lozano-Durán et al. 2020).

We here use concepts inspired by these prior works to complement the quasi-equilibrium model presented in § 2. We decompose the velocity $\tilde {\boldsymbol {u}}$ (this velocity is spatially filtered in the two-dimensional horizontal plane but not time filtered, and so it may still contain strong time and $y$-dependent deviations from the quasi-equilibrium profile $\bar {\boldsymbol {u}}$) according to $\tilde {\boldsymbol {u}} = \bar {\boldsymbol {u}} + \tilde {\boldsymbol {u}}''$, where $\tilde {\boldsymbol {u}}''$ is the non-equilibrium velocity to be modelled here. The deviations from the quasi-equilibrium velocity distribution $\bar {\boldsymbol {u}}$, such as the deviations visible in figures 2(a) and 2(d), can arise from a variety of sources such as time-dependent turbulent fluctuations and pressure gradients. The fastest changing pressure gradient fluctuations can induce oscillatory flow conditions even in the viscous sublayer, not unlike those involved in the Stokes first and second problems. Since for the quasi-laminar part of the flow in the viscous sublayer an analytical solution can be developed, we aim now to model that part of the non-equilibrium wall stress arising directly from the response of the laminar sublayer to rapid pressure gradient fluctuations. We denote the corresponding laminar velocity response as $\tilde {\boldsymbol {u}}_{l}''$ where subscript ‘$l$’ stands for laminar component.

In the viscous sublayer the linear terms of the Navier–Stokes equation dominate and hence we argue that $\tilde {\boldsymbol {u}}_{l}''$ obeys

(3.1)\begin{equation} \frac{\partial \tilde{\boldsymbol{u}}_l''}{\partial t} ={-} \frac{1}{\rho} \boldsymbol{\nabla}_h \tilde{p}'' + \nu \frac{\partial^2 \tilde{\boldsymbol{u}}_l''}{\partial y^2}, \end{equation}

where, as before, subscript $h$ represents the horizontal directions $x$ and $z$. The boundary conditions are $\tilde {\boldsymbol {u}}_l''(y=0,t) = 0$ and $\partial \tilde {\boldsymbol {u}}_l''/\partial y(y\to \infty,t)=0$ with the initial condition $\tilde {\boldsymbol {u}}_l''(y,t_0) = 0$. To simplify the problem, it is useful to define a ‘non-equilibrium free stream velocity’, $\tilde {\boldsymbol {u}}''_\infty$, defined as the velocity that would exist as an inviscid response to the non-equilibrium pressure gradient

(3.2)\begin{equation} \frac{\partial \tilde{\boldsymbol{u}}_\infty''}{\partial t} ={-}\frac{1}{\rho}\boldsymbol{\nabla}_h \tilde{p}'', \quad \to \quad \tilde{\boldsymbol{u}}_\infty''(t) = \int_{t_0}^t -\frac{1}{\rho} \boldsymbol{\nabla}_h \tilde{p}'' \,{\rm d} t'. \end{equation}

We can then use this velocity to eliminate the pressure gradient using the variable transformation $\hat {\boldsymbol {u}}(y,t)=\tilde {\boldsymbol u}_\infty ''(t) - \tilde {\boldsymbol u}_l''(y,t)$. The problem then reduces to the generalized Stokes problem where the wall velocity is $\tilde {\boldsymbol u}_\infty ''(t)$. From Schlichting & Gersten (2017) this has the solution

(3.3)\begin{equation} \hat{\boldsymbol u}(y,t) = \int_{t_0}^t \frac{\partial \tilde{\boldsymbol u}_\infty''}{\partial t} \textrm{erfc} \left( \frac{y}{2\sqrt{\nu(t-t')}} \right) {\rm d} t'. \end{equation}

Rewriting in terms of $\tilde {\boldsymbol u}_l''$ and $\boldsymbol {\nabla }_h \tilde {p}''$ gives

(3.4)\begin{equation} \tilde{\boldsymbol u}_l''(y,t)= \int_{t_0}^t \left(-\frac{1}{\rho}\boldsymbol{\nabla}_h \tilde{p}''\right) \textrm{erf} \left( \frac{y}{2\sqrt{\nu(t-t')}} \right) {\rm d} t', \end{equation}

from which the stress contribution can be obtained by differentiation, evaluation at $y=0$ and multiplication by $\nu$, and reads as follows:

(3.5)\begin{equation} \boldsymbol{\tau}_w''(t)= \sqrt{\nu/{\rm \pi}} \int_{t_0}^t - \frac{1}{\rho}\boldsymbol{\nabla}_h \tilde{p}^{\prime \prime}(t') (t - t')^{{-}1/2}\,{\rm d} t'. \end{equation}

Interestingly, we can use $\tilde {\boldsymbol u}_\infty ''$ to relate the non-equilibrium wall stress with the Caputo fractional derivative

(3.6)\begin{equation} \boldsymbol{\tau}_w'' ={-}\sqrt{\nu/{\rm \pi}} \varGamma(1/2) D_t^{1/2} (\tilde{\boldsymbol{u}}_\infty''), \end{equation}

where the Caputo fractional derivative of order $\alpha$ of a signal $v(t)$ is defined as (Samko 1993)

(3.7)\begin{equation} D_t^\alpha v(t) = \frac{1}{\varGamma(1-\alpha)} \int_{t_0}^t \frac{v^{(1)}(t')}{(t-t')^{\alpha}}{\rm d} t'. \end{equation}

In the equation above $0 < \alpha < 1$ is the order of the fractional derivative. The rapid wall-stress model uses $\alpha = 1/2$. Relating the wall stress with the Caputo fractional derivative is useful because an efficient numerical evaluation of this type of non-local integral operator is possible, as described in § 4.3.

4. Pressure decomposition and numerical implementations

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. 1992) 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 1996) 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é 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. (1996). 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. (2017) 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. (1996) 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. (2017), 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. (2017) 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 2015)

(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. Tests in equilibrium and non-equilibrium channel flow

To test the new wall model (with both quasi-equilibrium LaRTE and non-equilibrium components), LESs are conducted for statistically stationary channel flow as well as for channel flow with a large step change in the spanwise pressure gradient (referred to as SSPG).

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. 2005) 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. (2005), 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. 2017) 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 (2020) 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 (2015) 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. 2015), 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 (2020) 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 (2019), 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 (2015) 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 (2015).

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. (2016). 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. (2020). 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. (2020). 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. (2020) 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. (2020) 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. (2020) (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. (2020). 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. (2020) (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. (2020) 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. (2020).

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.

6. Summary and conclusions

We have introduced the LaRTE wall model, representing a quasi-equilibrium dynamics. The LaRTE model consists of an evolution equation for the friction velocity (and thus the wall stress) using a method similar to that introduced in Chung & Pullin (2009) where the law of the wall is utilized to rewrite the unsteady term as $\partial \boldsymbol {u}_\tau /\partial t$. Also similar to Chung & Pullin (2009) and Yang et al. (2015), the LaRTE model is based on the vertically integrated RANS-like thin boundary-layer equations. Remarkably, it is found that generalization to include the advection terms leads to a Lagrangian form rather than an Eulerian one. Moreover, the right-hand side of this equation is in the form of a term describing relaxation towards the stress at the wall-model height, with a relaxation time scale $T_s$. A priori testing based on DNS channel flow data shows that the relaxation time scale is consistent, whereas use of shorter time scales would be inconsistent, with the assumption of quasi-equilibrium.

The proposed formalism allows for separate modelling of non-equilibrium effects not captured by the LaRTE model. In § 3 we introduce a non-equilibrium wall model to capture quick transient pressure gradient effects. The approach is well suited for modelling the laminar Stokes layer observed in the literature for flows with a rapidly changing pressure gradient.

The LaRTE plus laminar non-equilibrium model is tested for (a) simple channel flow with a constant pressure gradient and (b) the sudden spanwise pressure gradient test case introduced in Lozano-Durán et al. (2020). These two cases were selected in order to verify that the LaRTE model elements indeed exhibit quasi-equilibrium behaviour, develop intuition regarding the model performance and examine its response in a case where the non-equilibrium effects are known to follow the Stokes laminar solution. Time signals, wall-stress contour plots and wall-stress PDFs reveal interesting and useful physical insight as to how the model operates. For example, time signals of the LaRTE-predicted wall stress are consistent with the idea that the friction velocity relaxes towards its equilibrium value at the relaxation time scale. Wall-stress contour plots show that the structures in the LaRTE model are rather large, of length of the order of or larger than the channel half-width. LaRTE thus implicitly averages out the turbulence leaving only large-scale fluctuations that are internally consistent with ‘quasi-equilibrium’, i.e. under conditions that are sufficiently averaged so that using the locally determined friction velocity one may assume the law of the wall to hold.

A comparison between Lagrangian and Eulerian versions of the relaxation towards equilibrium model show the Lagrangian version tends to advect wall-stress structures whereas the Eulerian version tends to thin out and elongate structures in the streamwise direction. Wall-stress PDFs show the distribution of wall-stress fluctuations of each component and how the non-equilibrium model diminishes in importance as Reynolds number increases. Differences between PDFs from the wall model and filtered DNS show that additional ingredients will be required to fully capture the turbulent fluctuations in wall stress, especially in the streamwise direction. For the SSPG test case, time signals of the plane-averaged mean spanwise and streamwise wall stresses show the model is in good agreement with the DNS data of Lozano-Durán et al. (2020) and that long after the application of the SSPG the model approaches the correct steady-state values.

In the LaRTE approach, closure for the total stress $\bar {\boldsymbol \tau }_\varDelta$ at the wall-model height is required. Also, the pressure gradient decomposition into various time scales involves some modelling choices. While a number of other options exist and could be explored, the corresponding choices used in this paper were justified through physical interpretation and the model's ability to yield good results for the test cases considered here.

In summary, the proposed LaRTE model represents a new framework for wall modelling. Specifically, the formal identification of quasi-equilibrium dynamics enables us to model the remainder non-equilibrium parts more rigorously. Moreover, the quasi-equilibrium portion can be further expanded beyond the case of smooth flat-plate surfaces to include, e.g. roughness effects and additional dependencies on other physical parameters such as streamline curvature, thermal and compressibility effects, etc. The choice of model for the turbulence stress at the wall-model height and the pressure gradient decompositions could also be further improved. We leave such extensions for future work. Also, the new wall model must be tested in other flows such as flows with streamwise pressure gradients, flow over a cylinder, steps and others. Finally, for flows with rapidly oscillating walls for e.g. drag reduction (Jung et al. 1992) in which the turbulent structure below $y=\varDelta$ is heavily affected by non-equilibrium scrambling effects, additional modelling will be required, corresponding to the ‘white’ region in figure 3.