Hostname: page-component-7857688df4-xgplm Total loading time: 0 Render date: 2025-11-14T08:16:55.273Z Has data issue: false hasContentIssue false

Wall-model response to departure from equilibrium in turbulent flows over rough surfaces

Published online by Cambridge University Press:  14 November 2025

Teresa Salomone*
Affiliation:
Department of Mechanical and Materials Engineering, Queen’s University , Kingston, ON K7L3N6, Canada Engineering Department, University of Campania Luigi Vanvitelli , Aversa 81031, Italy
Charles Meneveau*
Affiliation:
Mechanical Engineering, Johns Hopkins University , Baltimore, MD 21218-2625, USA
Giuliano De Stefano*
Affiliation:
Engineering Department, University of Campania Luigi Vanvitelli , Aversa 81031, Italy
Ugo Piomelli*
Affiliation:
Department of Mechanical and Materials Engineering, Queen’s University , Kingston, ON K7L3N6, Canada
*
Corresponding authors: Teresa Salomone, tereslm@stanford.edu; Charles Meneveau, meneveau@jhu.edu; Giuliano De Stefano, giuliano.destefano@unicampania.it; Ugo Piomelli, ugo@queensu.ca
Corresponding authors: Teresa Salomone, tereslm@stanford.edu; Charles Meneveau, meneveau@jhu.edu; Giuliano De Stefano, giuliano.destefano@unicampania.it; Ugo Piomelli, ugo@queensu.ca
Corresponding authors: Teresa Salomone, tereslm@stanford.edu; Charles Meneveau, meneveau@jhu.edu; Giuliano De Stefano, giuliano.destefano@unicampania.it; Ugo Piomelli, ugo@queensu.ca
Corresponding authors: Teresa Salomone, tereslm@stanford.edu; Charles Meneveau, meneveau@jhu.edu; Giuliano De Stefano, giuliano.destefano@unicampania.it; Ugo Piomelli, ugo@queensu.ca

Abstract

Most turbulent boundary-layer flows in engineering and natural sciences are out of equilibrium. While direct numerical simulation and wall-resolved large-eddy simulation can accurately account for turbulence response under such conditions, lower-cost approaches like wall-modelled large-eddy simulation often assume equilibrium and struggle to reproduce non-equilibrium effects. The recent ‘Lagrangian relaxation-towards-equilibrium’ (LaRTE) wall model (Fowler et al. 2022 J. Fluid Mech. vol. 934, 137), formulated for smooth walls, applies equilibrium modelling only to the slow dynamics that are more likely to conform to the assumed flow state. In this work, we extend the LaRTE model to account for wall roughness (LaRTE-RW) and apply the new model to turbulent flow over heterogeneous roughness and in accelerating and decelerating flows over rough surfaces. We compare predictions from the new LaRTE-RW model with those from the standard log-law equilibrium wall model (EQWM) and with experimental data to elucidate the turbulence response mechanisms to non-equilibrium conditions. The extended model transitions seamlessly across smooth-wall and fully rough regimes and improves prediction of the skin-friction coefficient, especially in recovering trends at roughness transitions and in early stages of pressure-gradient-driven flow acceleration or deceleration. Results show that LaRTE-RW introduces response delays that are beneficial when EQWMs react too quickly to disturbances, but it is less effective in flows requiring rapid response, such as boundary layers subjected to accelerating–decelerating–accelerating free stream conditions. These findings emphasize the need for further model refinements that incorporate fast turbulent dynamics not currently captured by LaRTE-RW.

Information

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2025. Published by Cambridge University Press

1. Introduction

Surface roughness is present in many turbulent boundary-layer flows encountered in engineering and natural sciences. Sources of surface roughness can include corrosion, ice accumulation and manufacturing processes, among others. Essentially, its presence alters the development of turbulent boundary layers by enhancing drag and modifying the turbulence-generation cycle. At sufficiently high Reynolds numbers, every surface must be considered rough as the viscous length scale becomes comparable or smaller than to the roughness length scale. The practical importance of roughness effects on fluid flow has led to extensive research over the years, beginning with the pioneering works of Nikuradse (Reference Nikuradse1933) and Colebrook (Reference Colebrook1939). The key findings have since been reviewed by Raupach, Antonia & Rajagopalan (Reference Raupach, Antonia and Rajagopalan1991), Jiménez (Reference Jiménez2004) and more recently Chung et al. (Reference Chung, Hutchins, Schultz and Flack2021), and include (i) an increase in drag (Nikuradse Reference Nikuradse1933; Colebrook Reference Colebrook1939), (ii) a reduction in the anisotropy of turbulent Reynolds stresses (Shafi & Antonia Reference Shafi and Antonia1995), (iii) an alteration of the near-wall region by disrupting streaky structures and (iv) an additional mechanism for turbulent kinetic energy (TKE) production (Raupach et al. Reference Raupach, Antonia and Rajagopalan1991; Finnigan Reference Finnigan2000). In particular, the increase in TKE production is linked to a rise in local skin friction compared with that of an equivalent smooth-wall flow. Increased skin friction is reflected in a downward shift of the logarithmic region of the normalized velocity profile. The magnitude of the downward shift is referred to as the ‘roughness function’ $\Delta U^+$ , where $U$ is the mean streamwise velocity, and the superscript $+$ indicates normalization using viscous scales, specifically the friction velocity $u_{\tau }$ and the kinematic viscosity $\nu$ .

Quadrant analysis has indicated that turbulence structures near the wall in a rough surface environment are altered, leading to stronger sweeps in the vicinity of the roughness elements and a corresponding scarcity of ejections (Krogstad, Antonia & Browne Reference Krogstad, Antonia and Browne1992). Schultz & Flack (Reference Schultz and Flack2007) reported similar trends, emphasizing that strong sweeps dominate Reynolds shear-stress contributions near rough surfaces, particularly in high $k^+ = k \, u_\tau / \nu$ regimes. While roughness significantly influences the flow characteristics in the near-wall region, particularly within the roughness sublayer – where roughness effects are pronounced – it has been suggested that the flow above the roughness sublayer follows outer-layer similarity, a concept known as ‘Townsend’s similarity hypothesis’ (Townsend Reference Townsend1976), if a proper velocity scale is used. Several studies support this hypothesis (Jiménez Reference Jiménez2004; Flack, Schultz & Shapiro Reference Flack, Schultz and Shapiro2005; Flack & Schultz Reference Flack and Schultz2014), especially when there is a significant scale-separation between the boundary-layer thickness ( $\delta$ ), and the mean-roughness height ( $k$ ), with a ratio of at least $\delta / k \gt 40$ (Jiménez Reference Jiménez2004), for plane-channel, pipe and zero-pressure gradient (ZPG) flows. However, while the concept of outer-layer similarity has proven robust under equilibrium conditions, roughness heterogeneity and pressure gradients can significantly alter the turbulent boundary-layer structure.

1.1. Roughness heterogeneity

In most applications roughness is not spatially uniform. Adjacent surface regions often have different roughness heights, such as in boundary layers over plant canopies or mixed terrains with built areas, woods and fields. These scenarios involve abrupt changes in the wall boundary-condition, which significantly affect the flow. A straightforward example of heterogeneous roughness is the abrupt transition from a two-dimensional rough patch perpendicular to the flow to a smooth surface, or vice versa. This configuration is common in fluid mechanics, particularly in atmospheric boundary layers (Cheng & Castro Reference Cheng and Castro2002; Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2004; Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2009) and meteorological flows (Garratt Reference Garratt1990). Antonia & Luxton (Reference Antonia and Luxton1971, Reference Antonia and Luxton1972), in their pioneering experimental work, analysed this particular configuration, considering both smooth-to-rough (Antonia & Luxton Reference Antonia and Luxton1971) and rough-to-smooth (Antonia & Luxton Reference Antonia and Luxton1972) transitions. As the flow undergoes the transition from one surface condition to another, initially only the region close to the wall is affected by the new wall condition. However, this effect propagates into the interior of the flow, leading to the development of an internal boundary layer (IBL). Within the IBL, only a small fraction – approximately 5 % of the IBL thickness – is in equilibrium with the new wall condition, in what is known as the equilibrium layer (Garratt Reference Garratt1990), while the remaining region is in a non-equilibrium state (Shir Reference Shir1972; Rao, Wyngaard & Coté Reference Rao, Wyngaard and Coté1974).

The rate at which the IBL develops also determines the spatial extent needed for the flow to reach equilibrium with the new wall condition. A common observation is that, after the transition, the different flow variables return to equilibrium at different downstream locations. For instance, the skin-friction adapts more rapidly to the new condition than the mean velocity and Reynolds stresses. The growth rate of the IBL is actually dependent on whether the transition is from rough-to-smooth surface or vice versa. Antonia & Luxton (Reference Antonia and Luxton1971, Reference Antonia and Luxton1972) showed that in a smooth-to-rough transition the flow outside the IBL is not influenced by the new surface condition and the readjustment of the flow variables takes place fairly rapidly (in less than 20 boundary-layer thicknesses downstream of the transition); in a rough-to-smooth transition, however, reaching a self-sustained state is a much slower process, with both the inner and outer layers being affected by the new surface condition. In the latter case, they attributed the slower readjustment to the higher turbulence content carried by the outer-layer structures in the rough-wall boundary layer, which is hardly dissipated on the following smooth wall.

The recovery of the skin-friction coefficient to equilibrium conditions is a crucial quantity to be examined in these complex scenarios. In fact, non-equilibrium effects are particularly evident in the discontinuities observed in the skin-friction coefficient. Moving from smooth to rough strips, there is a sharp increase in the skin-friction coefficient, due to the initial blockage effect caused by the roughness elements. Conversely, at the rough-to-smooth transition, there is a steep drop in the skin-friction coefficient, resulting from the sheltering effect, which causes the flow to detach from the roughness elements, creating a small separation region. However, after both transitions, the skin-friction coefficient relaxes to the equilibrium values characteristic of fully developed rough or smooth surfaces (Antonia & Luxton Reference Antonia and Luxton1971, Reference Antonia and Luxton1972; Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2004; Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2009; Hanson & Ganapathisubramani Reference Hanson and Ganapathisubramani2016; Li et al. Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019).

1.2. Combined effects of roughness and pressure gradients

The majority of studies on roughness effects have focused on plane channel flows and ZPG boundary layers (Jiménez Reference Jiménez2004; Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). The effects of pressure gradients for rough walls have also been analysed, showing that they induce changes in the boundary-layer structure, altering the distribution of velocity and turbulence within the flow (Coleman, Moffat & Kays Reference Coleman, Moffat and Kays1977; Cal et al. Reference Cal, Brzek, Johansson and Castillo2009).

The influence of a favourable pressure gradient (FPG) on a rough wall upon near-wall turbulence and TKE actually depends on the FPG strength, with weak FPG leading to higher near-wall TKE energy. Also, the application of an adverse pressure gradient (APG) reduces the streamwise normal Reynolds stress near the surface, by enhancing the outward convection of vortices formed by roughness elements (Shin & Jin Song Reference Shin and Jin Song2014). Stronger FPG can cause quasilaminarization, leading to a reduction in both TKE and skin friction, particularly for smooth (Spalart Reference Spalart1986) and low-roughness (Piomelli & Yuan Reference Piomelli and Yuan2013) surfaces. When roughness heights are significant, on the other hand, relaminarization is inhibited (Yuan & Piomelli Reference Yuan and Piomelli2014, Reference Yuan and Piomelli2015), as roughness alters the mechanism driving the generation and growth of turbulent eddies. However, near-wall wall-normal and spanwise fluctuations continue to persist during flow acceleration. Further investigations by Tachie & Shah (Reference Tachie and Shah2008) and Cal et al. (Reference Cal, Brzek, Johansson and Castillo2009) showed that, under FPG conditions, acceleration can decrease outer-layer Reynolds stresses while enhancing fluctuations near the wall, a behaviour not observed in smooth-wall flows.

The combined effect of APG and roughness intensifies the drag over the surface, amplifies the wake region and reduces the length of the log region (Tay, Kuhn & Tachie Reference Tay, Kuhn and Tachie2009; Tsikata & Tachie Reference Tsikata and Tachie2013). Due to a larger momentum deficit in the near-wall region caused by roughness, flow separation can occur earlier, and the enhanced mixing in rough walls makes the normal Reynolds stresses less sensitive to APG, compared with a smooth wall (Song & Eaton Reference Song and Eaton2002; Aubertine, Eaton & Song Reference Aubertine, Eaton and Song2004).

Notably, the studies mentioned above have explored the impact of pressure gradients on rough walls, usually concentrating on a single type of pressure gradient. There has been limited research into understanding how the combined presence of different pressure gradients, both favourable and adverse, impacts rough surfaces. In fact, few studies have examined smooth-wall flows with varying pressure gradients (Spalart & Watmuff Reference Spalart and Watmuff1993; Piomelli, Balaras & Pascarelli Reference Piomelli, Balaras and Pascarelli2000; Bobke et al. Reference Bobke, Vinuesa, Örlü and Schlatter2017; Devenport & Lowe Reference Devenport and Lowe2022; Fritsch et al. Reference Fritsch, Vishwanathan, Todd Lowe and Devenport2022b ; Parthasarathy & Saxton-Fox Reference Parthasarathy and Saxton-Fox2023) uncovering their non-equilibrium and non-self-similar nature, as well as their dependency on flow history.

When analysing pressure-gradient flows, it is crucial to differentiate between equilibrium and non-equilibrium cases. An equilibrium flow is defined as one in which the flow characteristics depend solely on local conditions, and the convective terms in the turbulent transport equation are negligible (Bradshaw & Ferriss Reference Bradshaw and Ferriss1965). This description accurately represents fully developed channel or pipe flows where the mean velocity profile and flow statistics do not vary in the streamwise direction. However, the introduction of pressure gradients complicates the flow dynamics, as the sink flow is the only true equilibrium flow under such conditions (Rotta Reference Rotta1962; Jones, Marusic & Perry Reference Jones, Marusic and Perry2001). Achieving equilibrium in external flows proves more challenging due to factors such as flow growth, entrainment and three-dimensional effects (Fritsch et al. Reference Fritsch, Vishwanathan, Todd Lowe and Devenport2022b ). Therefore, a more adaptable definition of equilibrium can be used, requiring instead that the mean velocity profile exhibits self-similarity at sufficiently high Reynolds numbers (Marusic et al. Reference Marusic, McKeon, Monkewitz, Nagib, Smits and Sreenivasan2010).

Maintaining this self-similarity becomes challenging when the flow experiences, for example, a strong APG, especially when it approaches separation. Pressure gradients can be characterized using the Clauser parameter $\beta = (\delta ^*/\tau _w) \, \textrm{d}p/\textrm{d}x$ , where $\delta ^*$ is the displacement thickness, $\tau _w$ is the wall shear stress and $\textrm{d}p/\textrm{d}x$ is the streamwise pressure gradient (Clauser Reference Clauser1954). For a near-equilibrium flow, the parameter $\beta$ must remain constant.

Research on flows with varying $\beta$ on smooth walls is limited, but studies on rough walls are even scarcer. Recently, Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022a ) extended their investigation of flow with varying pressure-gradient distributions from smooth walls (Fritsch et al. Reference Fritsch, Vishwanathan, Todd Lowe and Devenport2022b ) to rough walls. In their experiment, a boundary layer is formed beneath a NACA0012 airfoil in a wind tunnel. The airfoil is designed to produce a bidirectional, continually varying pressure-gradient distribution on a rough wall. Modifying the angle of attack of the airfoil causes either an FPG followed by APG or vice versa. Their findings indicated that, while the general effect of a pressure gradient on the mean and integrated boundary-layer parameters is similar for both rough-and smooth-walls, the boundary-layer flow responds more rapidly over rough walls. This behaviour is attributed to the increased dominance of small-scale, near-wall motions in rough-wall flows, which are known to be less reliant on flow history.

More recently, Volino & Schultz (Reference Volino and Schultz2023) conducted a similar study where the upper wall of the wind tunnel consisted of four flat plates that could be independently adjusted to create the desired pressure gradient. In their set-up, the transition between the FPG and APG regions was separated by a ZPG zone, allowing the flow to readjust to equilibrium conditions. They found that, in both smooth and rough walls, the FPG reduced the turbulence in the boundary layer and increased the number of turbulent structures, while lowering their inclination angle with respect to the wall. In the subsequent ZPG region the flow returned to its equilibrium state, with this adjustment occurring sooner in the rough wall than in the smooth one. In the APG region, opposite trends were observed, with changes due to the adverse gradient being significantly more pronounced in the rough wall compared with the smooth one. They also observed a general trend for the equivalent sand-grain roughness-height, $k_s$ , which represents a measure of the drag induced by the roughness, showing an increase of this parameter in FPG regions, and a decrease in APG ones. This observation is crucial for modelling purposes since $k_s$ is often an input in turbulence models for rough walls.

1.3. Wall models

Numerical simulation represents an essential tool to study rough-wall flows, but direct numerical simulation (DNS) and wall-resolved large-eddy simulation (LES) are limited to low or moderate-Reynolds-number flows (Jiménez Reference Jiménez2004). Practically, modelling the near-wall region is the only option to perform eddy-resolving simulations of high-Reynolds-number flows over rough walls at a reasonable computational cost.

Several wall models for LES have been developed and reviewed in the literature (Piomelli & Balaras Reference Piomelli and Balaras2002; Piomelli Reference Piomelli2008; Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016; Bose & Park Reference Bose and Park2018). The most common type of wall models are those known as ‘wall-stress models’, where the inner-layer solution is replaced by a single quantity, the wall shear stress $\tau _w$ , which is used as a wall-boundary condition for the LES field. In this context, the classical log-law equilibrium wall model (EQWM) (Deardorff Reference Deardorff1970; Schumann Reference Schumann1975; Piomelli et al. Reference Piomelli, Ferziger, Moin and Kim1989) remains the most commonly used, but it can result in errors for non-equilibrium flows (Bose & Park Reference Bose and Park2018).

The separation of near-wall effects from the outer flow allows wall models to focus on capturing the roughness impact near the surface, while applying universal turbulent scaling laws in the outer region. This approach is generally valid for equilibrium flows at high Reynolds numbers. However, under non-equilibrium conditions – such as varying pressure gradients or roughness heterogeneity – this assumption may no longer hold. Both near-wall and outer-layer dynamics can be influenced, raising questions about the accuracy of these wall models in such complex scenarios. While most studies have focused on equilibrium turbulent flows (e.g. flat-plate boundary layers, channels and pipes), real-world applications often involve non-equilibrium effects, such as mean flow three-dimensionality, unsteadiness, varying pressure gradients or heterogeneous roughness, which alter flow dynamics significantly. Addressing these gaps requires further exploration into the applicability of wall models under non-equilibrium conditions.

In fact, simulating non-equilibrium flows motivated research on more elaborate wall models. For example, a generalized Moody diagram (GMD) method (Meneveau Reference Meneveau2020) for wall-modelled LES (WMLES) directly expresses the wall stress as function of input parameters, all properly expressed in non-dimensional forms. The modelling procedure, which is referred here to as the ‘GMD model’, was constructed by introducing various fitting functions to the solution of the one-dimensional Reynolds-averaged Navier–Stokes (RANS) equation in its simple near-wall boundary-layer approximation with a mixing-length closure. The use of fitting functions allows the calculation of the wall shear stress without iterative procedures not only in equilibrium conditions, but also in the presence of mild and strong pressure-gradients, as well as spatially uniform surface roughness.

Based on integral modelling ideas pioneered in LES by Chung & Pullin (Reference Chung and Pullin2009), Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2015) proposed the idea of integral WMLES, where the turbulent boundary layer equations are integrated in the wall-normal direction, while an algebraic closure for the assumed velocity profile is utilized following the von Kàrmàn–Polhausen method. Specifically, the near-wall region is replaced by a velocity profile that is linear in the viscous sublayer and merges with a logarithmic layer in the core of the flow, with the inclusion of a linear correction to include pressure-gradient effects. The velocity profile is described by a series of coefficients that are determined using matching conditions at the inner–outer layer interface, which include boundary conditions, and the full integrated turbulent boundary layer equations. Although this methodology can be more accurate than EQWM, due to the inclusion of all terms in the momentum equation, it entails solving multiple nonlinear equations. In addition, it necessitates the specification of an empirically determined time scale for the exponential filtering of velocity inputs at the wall-model height.

The issues with the method proposed by Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2015) were addressed by Fowler, Zaki & Meneveau (Reference Fowler, Zaki and Meneveau2022) through the introduction of the ‘Lagrangian relaxation-towards-equilibrium’ (LaRTE) model. This model employs an evolution equation for the friction velocity, akin to the previous approach, but adopts a Lagrangian formulation instead of an Eulerian one. Importantly, it does not necessitate the a priori specification of a relaxation time scale, which instead arises as a direct consequence of the underlying transport equations. The LaRTE model allows the identification of quasiequilibrium (qEQ) dynamics that can be supplemented by additional non-equilibrium terms (Fowler, Zaki & Meneveau Reference Fowler, Zaki and Meneveau2023). However, the existing formulation of LaRTE does not explicitly account for roughness at the surface, since it is based on an assumed universal velocity profile in inner units that includes the linear viscous sublayer near the presumed smooth wall.

1.4. Objectives

The primary objective of this paper is to extend the LaRTE model (Fowler et al. Reference Fowler, Zaki and Meneveau2022, Reference Fowler, Zaki and Meneveau2023) to rough surfaces, enabling it to account for either homogeneous or heterogeneous roughness. The present extension for rough walls (LaRTE-RW) allows the model to capture the complex flow behaviour associated with transition between adjacent smooth and rough surfaces, which constitutes a challenge for existing models. A second objective is to examine and compare the responses of various wall models, including the standard log-law EQWM, the GMD model and the newly developed LaRTE-RW model, under non-equilibrium conditions. This study focuses on the accuracy of these models in estimating key flow quantities such as skin-friction coefficient and mean-velocity profile in non-equilibrium flows over rough surfaces. Third, applying several wall models to two different non-equilibrium classes of flows which are heterogeneous roughness and rough-wall boundary layers with FPG and APG, we aim to pinpoint some shortcomings of existing models, and indicate directions for future improvements.

The present heterogeneous roughness configuration consists of alternating rough and smooth strips oriented perpendicular to the mean flow direction. Specifically, our numerical simulations are conducted under flow conditions and domain characteristics comparable to the experimental work of Li et al. (Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019). Additionally, we compare present results with the experimental data from Antonia & Luxton (Reference Antonia and Luxton1971) to further evaluate the accuracy of the proposed wall models in predicting the readjustment of flow variables after a surface transition. In fact, this scenario represents a flow configuration where strong non-equilibrium effects originate at the wall. Moreover, for the streamwise pressure gradient cases, we perform simulations replicating the experiments of Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022a ). This configuration enables the simultaneous examination of roughness and pressure-gradient effects. Indeed, the variation in pressure gradient distributions results in flows exhibiting qEQ behaviour, particularly in the downstream regions of the test section.

The rest of this manuscript is structured as follows. The underlying LES and reference wall modelling methodologies are outlined in § 2, while § 3 introduces the proposed extension of the LaRTE model to rough walls. The results of preliminary simulations for homogeneous smooth and rough channel flows, used as validation benchmarks, are provided in § 4. The predictive capabilities of the new model for non-equilibrium flows are demonstrated in § 5, where the results of WMLESs are presented and discussed. The heterogeneous roughness configuration is considered in § 5.1, highlighting the strong non-equilibrium behaviour associated with roughness transitions. Moreover, in § 5.2, the simulations of turbulent boundary layers under streamwise pressure gradients, representing qEQ conditions, are shown. Finally, in § 6, we draw the main conclusions of this study, and offer recommendations for future research.

2. Methodology

2.1. Governing equations and numerical method

The governing equations for LES of incompressible flow are as follows:

(2.1) \begin{equation} \frac {\partial \hat {u}_i}{\partial x_i} = 0, \quad \frac {\partial \hat {u}_i}{\partial t} + \frac {\partial }{\partial x_{\!j}} (\hat {u}_i \hat {u}_{\!j}) = -\frac {\partial \hat {p}}{\partial x_i} + \nu {\nabla} ^2 \hat {u}_i - \frac {\partial \tau _{\textit{ij}}}{\partial x_{\!j}}, \end{equation}

where the hat ( $\hat {\boldsymbol{\cdot }}$ ) represents LES-filtered quantities. In this system, $\hat {u}_i$ are the filtered velocity components in the three spatial directions ( $x_1 \equiv x$ streamwise, $x_2 \equiv y$ wall-normal and $x_3 \equiv z$ spanwise), corresponding to the velocity components $u$ , $v$ and $w$ , respectively. The term $\hat {p} = \hat {P}/\rho$ represents the modified pressure, $\nu$ is the kinematic viscosity and $\tau _{\textit{ij}} = \widehat {u_i u_{\!j}} - \hat {u}_i \hat {u}_{\!j}$ stand for the unresolved subfilter-scale (SFS) stresses. These unknown terms are modelled using an eddy-viscosity assumption for the deviatoric part of the stress tensor,

(2.2) \begin{equation} \tau _{\textit{ij}}^d = \tau _{\textit{ij}} - \frac {1}{3} \delta _{\textit{ij}} \tau _{kk} = -2\nu _{SFS} \hat {S}_{\textit{ij}}, \end{equation}

where the resolved strain-rate tensor is given by

(2.3) \begin{equation} \hat {S}_{\textit{ij}} = \frac {1}{2} \left ( \frac {\partial \hat {u}_i}{\partial x_{\!j}} + \frac {\partial \hat {u}_{\!j}}{\partial x_i} \right ). \end{equation}

In this work, we employ the SFS model proposed by Vreman (Reference Vreman2004) with its standard model constant.

The LES governing equations are integrated in time using a fractional step method (Chorin Reference Chorin1968; Kim & Moin Reference Kim and Moin1985). The wall-normal diffusion is advanced using the Crank–Nicolson scheme, while a third-order Runge–Kutta scheme is applied to the remaining terms. Spatial discretization is performed with conservative second-order finite differences on a staggered grid. The code is parallelized using the message passing interface method, and has been rigorously validated through previous applications to similar boundary-layer flows (Keating et al. Reference Keating, Piomelli, Bremhorst and Nešić2004; Yuan & Piomelli Reference Yuan and Piomelli2015; Wu & Piomelli Reference Wu and Piomelli2018; Ambrogi, Piomelli & Rival Reference Ambrogi, Piomelli and Rival2022). Details on particular boundary conditions, placement of the wall model interface and roughness parameter $k_s$ , are provided in §§ 4, 5.1 and 5.2, for open-channel (or half-channel), heterogeneous roughness and streamwise-varying pressure gradient flow configurations, respectively.

2.2. Wall models

Wall-stress models used as boundary conditions for LES relate the wall shear stress to the resolved velocity in the outer region. Practically, a wall-model interface is set at a predefined wall-normal location, $y = \varDelta$ . In the subsequent applications, we will use $\varDelta = 5\,\%$ of the boundary-layer height for channel flow simulations. For boundary-layer applications with varying pressure gradients, the boundary-layer thickness varies in the streamwise direction (due to the developing nature of the boundary layer and the applied pressure gradients), so the relative position of the wall-model interface with respect to the boundary layer thickness will vary. At the strongest pressure gradient imposed, the interface is located between $13\,\%$ and $6\,\%$ of the boundary layer height, with $13\,\%$ being the location at the reference plane where there is ZPG and $6\,\%$ at the end of the domain.

2.2.1. The log-law EQWM

One of the most commonly used models of this type is the standard log-law EQWM, in which the wall stress is derived by numerically solving the logarithmic law-of-the-wall for smooth and rough surfaces, respectively,

(2.4) \begin{align} U^+ &= \kappa ^{-1} \log (\varDelta ^+) + B , \end{align}
(2.5) \begin{align} U^+ &= \kappa ^{-1} \log (\varDelta /k_s) + 8.5, \end{align}

where $\kappa =0.4$ is the von Kàrmàn constant, and $B=5$ . Equation (2.4) must be solved with iterative methods, whereas (2.5) can be explicitly solved once $k_s$ has been fixed. Consequently, the model for the wall stress magnitude $\tau _w(x,z,t) = \rho u_\tau ^2$ is

(2.6) \begin{align} \tau _w &= \rho \left [ \frac {\kappa U_{\textit{LES}}}{\log (\varDelta ^+) + \kappa B} \right ]^2, \end{align}
(2.7) \begin{align} \tau _w &= \rho \left [ \frac {\kappa U_{\textit{LES}}}{\log (\varDelta /k_s) + \kappa \, 8.5} \right ]^2, \end{align}

for smooth and rough walls, respectively, where $U_{\textit{LES}} = \| U_{\textit{LES}}^i \| = [\hat {u}(x,\varDelta ,z,t)^2 + \hat {w}(x,\varDelta ,z,t)^2]^{1/2}$ is the magnitude of the wall-parallel component of the LES velocity at the wall-model interface. As mentioned before, models of this type calculate the wall stress $\tau _w$ from outer-flow quantities; the computed stress is then used in the numerical model to close the momentum balance in the first grid cell.

2.2.2. The GMD model

The GMD model (Meneveau Reference Meneveau2020) relies on the full equilibrium momentum balance, where unsteady and advection terms are neglected. Very briefly, the streamwise RANS momentum equation in its boundary-layer approximation is integrated, while fitting functions are developed for its solution. The equation is integrated in a dimensionless framework, where two different Reynolds numbers are defined as

(2.8) \begin{equation} {\textit{Re}}_{\varDelta } = \frac {U_{\textit{LES}} \, \varDelta }{\nu },\, {\textit{Re}}_{\tau \varDelta }^{\textit{uf}} = \frac {u_\tau \,\varDelta }{\nu }. \end{equation}

Here, ${\textit{Re}}_{\varDelta }$ is the main input parameter of the wall model, and ${\textit{Re}}_{\tau \varDelta }^{\textit{uf}}$ is the fitted solution from which the friction velocity is calculated. Other input parameters include a dimensionless pressure gradient term $\psi _p = \rho ^{-1} \varDelta ^3 / \nu ^2 ( \partial \hat {p}_{\textit{LES}} / \partial x_i ) s_i$ , where the index $i$ spans the horizontal directions in the wall-parallel plane (i.e. $i = 1, 3$ ), and the repeated indices are intended to imply summation. Furthermore, $s_i = U_{\textit{LES}}^i / \|U_{\textit{LES}}^i \|$ is the $i$ th component of the unit vector in the direction of the resolved velocity at the wall-model interface. For rough surfaces, the model also requires the dimensionless roughness length, $k_s/\varDelta$ . Unlike most other EQWM formulations, and similarly to the standard Moody diagram approach, no iterations are required to determine the wall stress even for smooth-wall applications. Note that the GMD method can also be written in a more familiar form according to

(2.9) \begin{equation} \tau _w = \rho \frac {1}{2} \, c_{\!f}^{w\textit{m}}({\textit{Re}}_\varDelta ,\psi _p, k_s/\varDelta ) \, U_{\textit{LES}}^2, \end{equation}

where $c_{\!f}^{w\textit{m}}({\textit{Re}}_\varDelta , \psi _p, k_s/\varDelta ) = 2({\textit{Re}}_{\tau \varDelta }^{\textit{uf}}/{\textit{Re}}_{\varDelta })^2$ . Interested readers are referred to Meneveau (Reference Meneveau2020) for further details.

2.2.3. The LaRTE model

In the LaRTE model, the governing equation for the friction velocity vector $\boldsymbol{u}_\tau$ is derived by integrating the RANS equations for the wall-parallel momentum components. The integration takes place within the layer beneath $\varDelta$ , producing a governing equation that describes the qEQ evolution of the wall stress. In fact, the estimation of this qEQ wall stress is governed by a relaxation time scale, $T_s$ , that arises from the derivation without adjustable parameters (Fowler et al. Reference Fowler, Zaki and Meneveau2022). The temporal relaxation process takes place in a Lagrangian frame, moving at a velocity $\boldsymbol{V}_\tau$ that also arises from the fundamental integral momentum equation. The LaRTE formulation enables applying the equilibrium modelling only to the dynamics occurring at time scales smaller than $T_s$ , i.e. self-consistently to where equilibrium modelling is justified.

Additional contributions to the total wall stress must be independently modelled. For smooth surfaces, the response of laminar viscous sublayers developing in response to rapidly (non-equilibrium) applied pressure gradients was modelled using analytical solutions to near-wall Stokes-layer dynamics (Fowler et al. Reference Fowler, Zaki and Meneveau2022). The additional non-equilibrium turbulence wall stress term subsequently proposed by Fowler et al. (Reference Fowler, Zaki and Meneveau2023) describes wall stress evolution at intermediate time scales, namely, between the viscous time scale ( $t_{\nu } \sim \nu /u_{\tau }^2$ ) and $T_s$ .

The turbulent non-equilibrium (turbNEQ) term is modelled using the attached eddy hypothesis of Townsend (Reference Townsend1976), which describes the wall-attached eddies as large-scale structures that influence the flow near the wall. These eddies are assumed to produce a velocity fluctuation $u^{\prime \,i}_\varDelta = U_{\textit{LES}}^i - u_{\tau , i} F(y^+)$ , where $F(y^+)$ is the inner-similarity function, which will be discussed in the next section. In our study, we have adopted the turbNEQ model as presented in the literature, which assumes a linear velocity profile in the vicinity of the wall. This assumption, however, may not fully account for the effects of roughness, which are critical to our analysis. As such, future developments may need to refine this model to better capture the influence of rough surfaces on the turbNEQ contribution. The linear profile assumption is valid within a layer characterized by a thickness $l_s^+$ , where $l_s^+ = 12$ is the typical value for the viscous sublayer. In this layer, the velocity profile is assumed to transition linearly from the wall to the flow. The stress contribution from the turbNEQ component is then computed as $\tau _{w,i}^\prime (x,z,t) = u_{\tau } u^{\prime \,i}_\varDelta /l_s^+$ (Fowler et al. Reference Fowler, Zaki and Meneveau2023).

For cases where laminar non-equilibrium contributions can be neglected, the total wall stress is evaluated as the superposition of the qEQ wall stress, $\overline {\tau }_{w,i}(x,z,t) = u_\tau (x,z,t) u_{\tau ,i}(x,z,t)$ ( $i=1,3$ ) and the turbNEQ component, $\tau _w^\prime (x,z,t)$ , according to

(2.10) \begin{equation} \tau _{w,i}(x,z,t) = \overline {\tau }_{w,i}(x,z,t) + \tau _{w,i}^\prime (x,z,t), \end{equation}

where, generally, $(\overline {\boldsymbol{\cdot }})$ represents qEQ model terms, while a single prime denotes turbNEQ terms. Importantly, this formulation enables a separation of the wall stress into qEQ and non-equilibrium components, where the latter accounts for additional stress contributions not arising from the qEQ dynamics that are described through the LaRTE modelling procedure. The interested reader is referred to Fowler et al. (Reference Fowler, Zaki and Meneveau2023) for a detailed description of the turbNEQ component. In the next section, the proposed extension of the qEQ component, $\overline {\tau }_w(x,z,t)$ , of the LaRTE model for rough walls is presented.

3. The LaRTE model: roughness extension

Following the derivation of the LaRTE model for smooth walls (Fowler et al. Reference Fowler, Zaki and Meneveau2022), we assume that, in the flow region between the rough wall and the wall-model interface ( $0 \lt y \lt \varDelta$ ), a horizontal qEQ mean velocity profile, $\overline {u}_{\textit{si}}$ , exists, and can be written as $\overline {\boldsymbol{u}}=\overline {u}_{s1}\boldsymbol{ i}+\overline {v}{\boldsymbol{j}}+\overline {u}_{s3}\boldsymbol{k}$ , where $\overline {u}_{\textit{si}}$ ( $i=1,3$ ) is the two-dimensional velocity in the wall-parallel plane while $\overline {v}$ is the wall normal velocity component, and $\boldsymbol{i}$ , $\boldsymbol{j}$ and $\boldsymbol{ k}$ are the usual Cartesian unit vectors. The former is expressed in terms of inner scaling according to

(3.1) \begin{equation} \overline {u}_{\textit{si}} (x,y,z,k_s^+,t) = u_{\tau i}(x,z,k_s^+,t) F(y^+,k_s^+), \end{equation}

where $u_{\tau i}(x, z, k_s^+, t) = u_{\tau } s_i$ , with $u_{\tau } = \| u_{\tau i} \|$ . The friction velocity vector $u_{\tau } s_i$ is a function of spatial coordinates along the wall, $x$ , $z$ and time $t$ , as well as surface roughness, which is characterized by the equivalent sand-grain roughness height $k_s^+ = k_s u_{\tau }/\nu$ . The inner-similarity function, $F(y^+,k_s^+)$ , depends on $y^+ = y u_{\tau }/\nu$ and contains a log-law region that merges with a linear region in the viscous sublayer. Additionally, for roughness, a dependency on $k_s^+$ must be added to the original formulation, so as to yield

(3.2) \begin{equation} F(y^+, k_s^+) = f(y^+) - f_{\textit{NW}}(y^+) \Delta U^+(k_s^+), \end{equation}

where $f(y^+)$ is the smooth inner-similarity function as invoked in Fowler et al. (Reference Fowler, Zaki and Meneveau2022), while $f_{\textit{NW}}(y^+) \Delta U^+(k_s^+)$ accounts for roughness effects through the roughness function $\Delta U^+$ . The function $f_{\textit{NW}}(y^+)$ is active in the near-wall region to drive the roughness contribution to zero in this region. When $k_s^+ = 0$ , $F(y^+, k_s^+) = f(y^+)$ and the smooth inner-similarity function is recovered. On the contrary, when $k_s^+ \gt 0$ , the velocity profile shifts downward as an effect of the increasing drag. In this work, the amount by which the profile is shifted is obtained from the empirical Colebrook formula (Colebrook Reference Colebrook1939),

(3.3) \begin{equation} \Delta U^+ = \kappa ^{-1} \log (1+0.26 \, k_s^+). \end{equation}

Equation (3.3) was preferred to other possible relations because it allows analytical evaluation of the derivatives. The dynamical evolution equation for $u_{\tau i}(x,z,k_s^+,t)$ is derived from the integration of the wall-parallel RANS momentum equation for the velocity $\overline {u}_{\textit{si}}$ , which reads

(3.4) \begin{align} \frac {\partial \overline {u}_{\textit{si}}}{\partial t} + \partial _{\!j} (\overline {u}_{\textit{si}} \overline {u}_{\textit{sj}}) + \partial _y (\overline {v} \, \overline {u}_{\textit{si}}) &= - \frac {1}{\rho } \partial _i \overline {p} + \frac {\partial }{\partial y} \left [ (\nu + \nu _T) \frac {\partial \overline {u}_{\textit{si}}}{\partial y} \right ] \nonumber \\ &\quad + \partial _{\!j} [(\nu + \nu _T) (\partial _i\overline {u}_{\textit{sj}} + \partial _{\!j} \overline {u}_{\textit{si}})] ,\end{align}

where $i = 1,3$ are the wall-parallel directions, $\nu _T$ is the eddy viscosity and $\overline {p}(x,z,t)$ is the qEQ pressure that, within the boundary layer-approximation, does not depend on $y$ . As outlined in Fowler et al. (Reference Fowler, Zaki and Meneveau2022), the incorporation of the diffusion term $\partial _{\!j} [(\nu + \nu _T) (\partial _i\overline {u}_{\textit{sj}} + \partial _{\!j} \overline {u}_{\textit{si}})$ can be avoided in practical LES, as diffusion is artificially introduced by the time integration scheme. Therefore, the diffusion term is neglected in the derivation of the equation for the friction-velocity vector.

Equation (3.4), without the diffusion term, is integrated from the wall, $y=0$ , to the wall-model interface, $y=\varDelta$ . Here, only the treatment of the Eulerian time derivative term is described, whereas the advective term will be treated in Appendix A.

Following the approach used by Fowler et al. (Reference Fowler, Zaki and Meneveau2022), we begin by defining the following integrals:

(3.5) \begin{eqnarray} \delta _\varDelta ^* = \int _{0}^\varDelta \left [1-\frac {\overline {u}_{\textit{si}} (y)}{\overline {u}_{\textit{si}}(\varDelta )}\right ]\textrm{d}y \Rightarrow \int _{0}^{\varDelta } F(y^+, k_s^+)\textrm{d}y = F_\varDelta (\varDelta -\delta _\varDelta ^*), \end{eqnarray}

where $F_\varDelta = F(\varDelta ^+, k_s^+)$ , and

(3.6) \begin{align} \theta _\varDelta = \int _{0}^\varDelta \frac {\overline {u}_{\textit{si}}(y)}{\overline {u}_{\textit{si}}(\varDelta )} \left [1-\frac {\overline {u}_{\textit{si}}(y)}{\overline {u}_{\textit{si}}(\varDelta )}\right ] \textrm{d}y \Rightarrow \int _{0}^\varDelta F^2(y^+, k_s^+)\textrm{d}y = F^2_\varDelta (\varDelta -\delta _\varDelta ^*-\theta _\varDelta ). \end{align}

The above two equations determine the ‘cell displacement thickness’ $\delta _\varDelta ^*$ , and the ‘cell momentum thickness’ $\theta _\varDelta$ , respectively.

Next, consider the integration of the Eulerian time derivative,

(3.7) \begin{eqnarray} \frac {\partial }{\partial t} \int _{0}^\varDelta \overline {u}_{\textit{si}} \,\textrm{d}y &=& \frac {\partial }{\partial t} \left [ \int _{0}^\varDelta u_{\tau i} F(y^+, k_s^+) \,\textrm{d}y \right ] \nonumber \\ &=& \frac {\partial u_{\tau i}}{\partial t} \int _{0}^\varDelta F(y^+, k_s^+) \,\textrm{d}y + u_{\tau i} \int _{0}^\varDelta \frac {\partial }{\partial t} F(y^+, k_s^+) \,\textrm{d}y. \end{eqnarray}

Substituting (3.5) in the expression for the first integral, we obtain

(3.8) \begin{eqnarray} \frac {\partial }{\partial t} \int _{0}^\varDelta \overline {u}_{\textit{si}} \,\textrm{d}y &=& \frac {\partial u_{\tau i}}{\partial t} F_\varDelta (\varDelta -\delta _\varDelta ^*) + u_{\tau i} \int _{0}^\varDelta \frac {\partial }{\partial t} F(y^+, k_s^+) \,\textrm{d}y. \end{eqnarray}

The term $\int _{0}^\varDelta \partial _t F(y^+,k_s^+) \,\textrm{d}y$ can be separated into the smooth and rough contributions, as follows:

(3.9) \begin{eqnarray} \int _{0}^\varDelta \partial _t F(y^+,k_s^+) \,\textrm{d}y &=& \int _{0}^\varDelta \partial _t f(y^+) \,\textrm{d}y - \int _{0}^\varDelta \partial _t \left [f_{\textit{NW}}(y^+) \Delta U^+(k_s^+)\right ] \,\textrm{d}y. \end{eqnarray}

The first integral on the right-hand side can be handled similarly to Fowler et al. (Reference Fowler, Zaki and Meneveau2022), using the notation $f^\prime (y^+) = \textrm{d} f(y^+)/\textrm{d} y^+$ and $f_\varDelta = f(\varDelta ^+)$ ,

(3.10) \begin{eqnarray} \int _{0}^\varDelta \partial _t f(y^+) \,\textrm{d}y &=& \int _{0}^{\varDelta ^+} f^\prime (y^+) \frac {\partial y^+}{\partial t} \,\textrm{d}y^+ \frac {\nu }{u_{\tau }} \nonumber \\ &=& \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} f_\varDelta \varDelta - \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} \int _{0}^\varDelta f(y^+) \,\textrm{d}y. \end{eqnarray}

Moreover, applying the product rule to the second integral in (3.9), one obtains

(3.11) \begin{align} \int _{0}^\varDelta \partial _t \left [f_{\textit{NW}}(y^+) \Delta U^+(k_s^+)\right ] \,\textrm{d}y &= \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} f_{\textit{NW},\varDelta } \Delta U^+ \varDelta - \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} \int _{0}^\varDelta f_{\textit{NW}} \Delta U^+ \,\textrm{d}y \nonumber \\ &\quad + \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} k_s^+ \left (\Delta U^+\right )^\prime \int _{0}^\varDelta f_{\textit{NW}} \,\textrm{d}y, \end{align}

Table 1. Definitions of the momentum and displacement thicknesses required by the new LaRTE-RW model.

where $ (\Delta U^+ )^\prime = d \Delta U^+(k_s^+)/\textrm{d} k_s^+$ represents the sensitivity of the roughness function $\Delta U^+$ to the value of $k_s^+$ . Expression (3.9) then becomes

(3.12) \begin{align} \int _{0}^\varDelta \partial _t F(y^+,k_s^+) \,\textrm{d}y = \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} F_{\varDelta } \delta ^* - \frac {1}{u_{\tau }} \frac {\partial u_{\tau }}{\partial t} k_s^+ \left (\Delta U^+\right )^\prime f_{\textit{NW},\varDelta } (\varDelta - \delta ^*_{\varDelta , \textit{NW}}), \end{align}

where $\delta ^*_{\varDelta , \textit{NW}}$ is the cell displacement thickness defined using the function $f_{\textit{NW}}(y^+)$ (see table 1). Consequently, (3.7) leads to

(3.13) \begin{align} \frac {\partial }{\partial t} \int _{0}^\varDelta \overline {u}_{\textit{si}} \,\textrm{d}y& = \frac {\partial u_{\tau i}}{\partial t} F_\varDelta \left (\varDelta -\delta _\varDelta ^* \right ) \nonumber \\ &\quad + s_i \frac {\partial u_\tau }{\partial t} \left [F_\varDelta \delta _\varDelta ^* - k_s^+ (\Delta U^+)^\prime f_{\textit{NW},\varDelta } \left (\varDelta -\delta ^*_{\varDelta , \textit{NW}}\right )\right ]\!, \end{align}

that, after some manipulation, can be expressed as

(3.14) \begin{align} \frac {\partial }{\partial t} \int _{0}^\varDelta \overline {u}_{\textit{si}} \,\textrm{d}y& = \Delta F_{\varDelta , \textit{HR}} \frac {\partial u_{\tau i}}{\partial t} \nonumber \\ &\quad - u_{\tau } \left [ F_\varDelta \delta _\varDelta ^* - k_s^+ (\Delta U^+)^\prime f_{\textit{NW},\varDelta } \left (\varDelta -\delta ^*_{\varDelta , \textit{NW}}\right ) \right ] \frac {\partial s_i}{\partial t}, \end{align}

where $F_{\varDelta , \textit{HR}} = F_\varDelta - k_s^+ (\Delta U^+)^\prime f_{\textit{NW},\varDelta } (1 - \delta ^*_{\varDelta , \textit{NW}}/\varDelta )$ . Notice that when the wall is smooth, i.e. $k_s^+ = 0$ and $F(y^+, k_s^+) = f(y^+)$ , the previous equation reduces to

(3.15) \begin{eqnarray} \frac {\partial }{\partial t} \int _{0}^\varDelta \overline {u}_{\textit{si}} \,\textrm{d}y &=& \varDelta f_\varDelta \frac {\partial u_{\tau i}}{\partial t} - u_{\tau } f_\varDelta \delta ^*_{\varDelta S} \frac {\partial s_i}{\partial t}. \end{eqnarray}

In fact, (3.15) corresponds to the Eulerian time-derivative originally derived for smooth walls in the LaRTE model (Fowler et al. Reference Fowler, Zaki and Meneveau2022).

Using the same procedure for the advection term in (3.4), by defining the advection velocity of the friction-velocity vector $V_{\tau j}$ , and the relaxation time scale $T_s$ as

(3.16) \begin{align} V_{\tau j} &= \left ( 1- \frac {\delta _\varDelta ^*}{\varDelta }-\frac {\theta _\varDelta }{\varDelta } \right ) \frac {F_\varDelta ^2}{F_{\varDelta , \textit{HR}}} u_{\tau j} , \end{align}
(3.17) \begin{align} T_s &= F_{\varDelta , \textit{HR}} \frac {\varDelta }{u_\tau }, \end{align}

the new LaRTE-RW equation, which accounts for roughness, is given by

(3.18) \begin{eqnarray} \frac {\partial u_{\tau i}}{\partial t} &+& V_{\tau j} \partial _{\!j} ( u_{\tau i}) = \frac {1}{T_s} \!\left [ \frac {1}{u_\tau } \left (\! -\frac {\varDelta }{\rho } \partial _i \overline {p} + \overline {\tau }_{\Delta i} \!\right ) -u_{\tau i}\!\right ] + u_{\tau } \frac {F_\varDelta }{F_{\varDelta , \textit{HR}}} \frac {\delta _\varDelta ^* }{\varDelta } \frac {\partial s_i}{\partial t} + \mathcal{F}_{1 i} \frac {F_\varDelta }{F_{\varDelta , \textit{HR}}} \frac {\theta _\varDelta }{\varDelta } \nonumber \\ &+& 2 \frac {f_\varDelta }{F_{\varDelta , \textit{HR}}} \mathcal{F}_{2 i} \left ( 1 - \frac {\delta ^*_{\Delta S}}{\varDelta } - \frac {\theta _{\Delta f \, f_{\textit{NW}}}}{\varDelta }\right ) -2 \frac {f_{\textit{NW},\varDelta } \Delta U^+}{F_{\varDelta , \textit{HR}}} \mathcal{F}_{2 i} \left ( 1 - \frac {\delta ^*_{\varDelta , \textit{NW}}}{\varDelta } - \frac {\theta _{\varDelta , \textit{NW}}}{\varDelta }\right ) \nonumber \\ &-& \frac {F_\varDelta }{F_{\varDelta , \textit{HR}}} \mathcal{F}_{2 i} \left (1 - \frac {\delta ^*_{\varDelta , \textit{NW}}}{\varDelta }\right ) - u_{\tau } k_s^+ \frac { d (\Delta U^+)}{\textrm{d} k_s^+} \frac {f_{\textit{NW},\varDelta }}{F_{\varDelta , \textit{HR}}} \left (1-\frac {\delta ^*_{\varDelta , \textit{NW}}}{\varDelta } \right ) \frac {\partial s_i}{\partial t} . \end{eqnarray}

In this equation, the functions $\mathcal{F}_{1 i}$ and $\mathcal{F}_{2 i}$ are defined as

(3.19) \begin{align} \mathcal{F}_{1 i} &= u_\tau u_{\tau i} \, (\partial _{\!j} s_{\!j}) \, F_\varDelta , \end{align}
(3.20) \begin{align} \mathcal{F}_{2 i} &= u_{\tau i} u_{\tau j} \, (\partial _{\!j} \Delta U^+) \, f_{\textit{NW},\varDelta }. \end{align}

The first three terms on the right-hand side of (3.18) were already present in the previous smooth formulation (Fowler et al. Reference Fowler, Zaki and Meneveau2022), where the term $\mathcal{F}_{1 i} \theta _\varDelta / \varDelta$ was actually added by Fowler et al. (Reference Fowler, Zaki and Meneveau2023), whereas the additional four terms appear from the current roughness extension. Here, the displacement and momentum thicknesses, $\delta ^*_{\varDelta }$ and $\theta _\varDelta$ , are evaluated using the new similarity function $F(y^+, k_s^+)$ . The terms $\delta ^*_{\Delta S}$ and $\delta ^*_{\varDelta , \textit{NW}}$ are determined using the smooth function $f(y^+)$ and the near-wall function $f_{\textit{NW}}(y^+)$ , respectively. The momentum thickness $\theta _{\Delta ff_{\textit{NW}}}$ involves the product $f(y^+)f_{\textit{NW}}(y^+)$ , while $\theta _{\varDelta , \textit{NW}}$ utilizes $f_{\textit{NW}}(y^+)$ . A summary of the displacement and momentum thicknesses appearing in (3.18) is shown in table 1.

Equation (3.18) describes the qEQ development of the wall stress, $\overline {\tau }_{w i} = u_{\tau } u_{\tau i}$ , over the time scale $T_s$ , under the influence of $\overline {\tau }_{\varDelta i}$ . This equation incorporates roughness effects through a modification of $V_{\tau j}$ and $T_s$ (compared with the previous smooth formulation), and through the functions $\mathcal{F}_{2 i}$ and $k_s^+ d (\Delta U^+)/\textrm{d} k_s^+$ in the last term (accounting for roughness heterogeneity). In fact, when roughness heterogeneity is present, the gradient $\partial \Delta U^{+} / \partial x_{\!j}$ appearing in the function $\mathcal{F}_{2 i}$ (3.20) activates the fourth, fifth and sixth terms on the right-hand side of (3.18). The final (seventh) term involves the quantity $k_s^+ \, d (\Delta U^+)/\textrm{d} k_s^+$ , which represents the sensitivity of the roughness function $\Delta U^+$ to variations in $k_s^+$ in rough-wall regions (i.e. where $k_s^+ \gt 0$ ). Under the Colebrook formulation, this term is proportional to $1/k_s^+$ in the fully rough regime.

For homogeneous smooth walls, where the terms $k_s^+ \, d (\Delta U^+)/\textrm{d} k_s^+$ and $\partial _{\!j} \Delta U^+$ vanish, (3.18) simplifies to

(3.21) \begin{equation} \frac {\partial u_{\tau i}}{\partial t} + V_{\tau j} \partial _{\!j} ( u_{\tau i}) = \frac {1}{T_s} \left [ \frac {1}{u_\tau } \left ( -\frac {\varDelta }{\rho } \partial _i \overline {p} + \overline {\tau }_{\Delta i} \right ) -u_{\tau i}\right ] + u_{\tau } \frac {\partial s_i}{\partial t} \frac {\delta _\varDelta ^* }{\varDelta } + \mathcal{F}_{1 i} \frac {\theta _\varDelta }{\varDelta }. \end{equation}

This equation, which is analogous to the corresponding equation in Fowler et al. (Reference Fowler, Zaki and Meneveau2022), excluding the diffusion term and including the $\mathcal{F}_{1 i} \theta _\varDelta /\varDelta$ term, differs primarily in that $V_{\tau j}$ and $T_s$ are defined here by using the rough-wall inner-similarity function $F(\varDelta ^+,k_s^+)$ .

The friction velocity relaxes to $ -\varDelta /\rho \, \partial _i \overline {p} + \overline {\tau }_{\varDelta i}$ on the right-hand side of (3.18), which requires a closure that accounts for surface roughness. We follow the approach proposed by Fowler et al. (Reference Fowler, Zaki and Meneveau2023), where the stress at the interface is computed using (2.9), and the pressure is taken directly from the LES field, i.e. $ \overline {p} = \hat {p}_{\textit{LES}}$ . In contrast to Fowler et al. (Reference Fowler, Zaki and Meneveau2023), where the contribution of the fast laminar Stokes layer was subtracted from the LES pressure input, we retain the full LES pressure. This choice is justified by the rough-wall configuration considered here, in which a laminar Stokes layer is not expected to exist.

Moreover, while Fowler et al. (Reference Fowler, Zaki and Meneveau2022) use a filtered pressure gradient – obtained by applying a temporal filter over scales of order $ 3T_s$ – as input to the wall stress model, in order to simplify implementation, we use the instantaneous LES pressure gradient. We have verified that the simulation results are not sensitive to this modelling choice. Consequently, using the mild pressure gradients fit (Meneveau Reference Meneveau2020), the closure for the stress at the interface reads

(3.22) \begin{equation} -\frac {\varDelta }{\rho } \, \partial _i \overline {p} + \overline {\tau }_{\varDelta i} = \frac {1}{2} \, c_{\!f}^{w\textit{m}}({\textit{Re}}_\varDelta , \chi , k_s/\varDelta ) \, U_{\textit{LES}}^2 \, s_i, \end{equation}

where $ c_{\!f}^{w\textit{m}}$ is the modelled skin-friction coefficient evaluated with the mild-pressure gradients fit utilizing the non-dimensional pressure gradient parameter $\chi$ that is defined in Appendix C, where the algorithm used to compute $ c_{\!f}^{w\textit{m}}$ is detailed.

It is worth noting that the equilibrium model $ U^+ = \kappa ^{-1} \log (\varDelta /k_s) + 8.5$ assumes fully rough conditions. While this model is well-suited for fully rough regimes, the GMD model offers greater flexibility as it can be applied to transitionally rough cases. The GMD model includes a transition between smooth wall and fully rough behaviour, with the associated ‘Moody diagram’ displaying $ C_{\!f} \sim {\textit{Re}}^{-1}$ behaviour at low Reynolds numbers and becoming horizontal in the fully rough regime (see figure 10 of Meneveau Reference Meneveau2020). The LaRTE-RW model uses the GMD formulation to compute stress at the interface, making it applicable to transitionally rough conditions as well. However, it is important to note that the relationship between $ \Delta U^+$ and $ k_s^+$ is based on the Colebrook formula, which is universal in the fully rough regime, but it is not in transitionally rough cases.

The LaRTE-RW evolution (3.18) is discretized in time using a forward Euler method, allowing for the explicit update of $u_{\tau i}$ . The inherent numerical diffusion from this discretization scheme led to the omission of the wall-parallel diffusion term. Spatial gradients resulting from the advection term in the Lagrangian time derivative are computed using a semi-Lagrangian approach, where the right-hand side of the equation is evaluated at the upstream position with bilinear spatial interpolation. Some more details on the numerical solution of (3.18) can be found in the discussion provided in Fowler et al. (Reference Fowler, Zaki and Meneveau2022) for the corresponding evolution equation for smooth walls.

4. Validation

To provide a preliminary assessment of the LaRTE-RW model’s accuracy, WMLESs were carried out for statistically steady open-channel (or half-channel) flow, considering both smooth and homogeneous rough-wall conditions. The computational domain was $L_x \times L_y \times L_z$ = $6.4 \delta \times \delta \times 3.2 \delta$ , with a uniform grid resolution of $192 \times 60 \times 94$ points. Periodic boundary conditions were applied in the streamwise and spanwise directions, while the top boundary was treated with symmetry conditions. The bottom-wall boundary condition depended on the specific wall stress model employed. In particular, the LaRTE-RW model, which computes wall stress according to (3.18), was compared with both the GMD model and log-law EQWM.

The wall-model interface was located at $ \varDelta /\delta = 0.05$ , ensuring that three grid points were below this interface, following the recommendations of Kawai & Larsson (Reference Kawai and Larsson2012). Although the flow below the interface is unphysical, due to the limited resolution, having points in this region helps facilitate eddy generation and allows the mixing length to exceed the grid size.

In both GMD and LaRTE-RW models, the input velocity from the LES field, $ U_{\textit{LES}}$ , was filtered in time using a single-sided exponential filter of the form

(4.1) \begin{equation} U_{\textit{LES}} = \alpha \, U_{\textit{LES}}^{n} + (1 - \alpha ) \, U_{\textit{LES}}^{n-1}, \end{equation}

where $ U_{\textit{LES}}^{n}$ is the instantaneous horizontal velocity at the wall-model height at time step $ n$ , and $ U_{\textit{LES}}^{n-1}$ is the temporally filtered velocity at previous time step. The smoothing coefficient is defined as

(4.2) \begin{align} \alpha = \frac {\Delta t}{T_{\textit{avg}}}, \end{align}

where $ \Delta t$ is the LES time step, and the averaging window is $ T_{\textit{avg}} = 3 T_s$ . This exponential filter stabilizes the input to the wall model by attenuating under-resolved, high-frequency fluctuations, while still preserving physically relevant unsteadiness in the large-scale motions. In contrast, the log-law model utilized the instantaneous velocity at the wall-model interface. Simulations with a sand-grain roughness height of $k_s^+ \simeq 130$ were performed, along with smooth-wall simulations, for comparison with the results from Lee & Moser (Reference Lee and Moser2015).

Figure 1. Inner-scaled mean velocity profiles for smooth and rough homogeneous surfaces, with the inset showing the corresponding roughness function. The DNS data from Lee & Moser (Reference Lee and Moser2015) are used for the smooth wall at ${\textit{Re}}_{\tau } = 5200$ , while rough wall simulations are conducted at a roughness Reynolds number of $k_s^+ = 130$ .

The simulated inner-scaled mean-velocity profiles are shown in figure 1. Under equilibrium conditions, both the LaRTE-RW and GMD models produced similar results, with no notable deviations. All velocity profiles closely follow the smooth log-law solution, represented by the solid black line. The smooth-wall simulations show a slight overprediction of the wake region amplitude compared with reference DNS data (Lee & Moser Reference Lee and Moser2015), a discrepancy also observed by Fowler et al. (Reference Fowler, Zaki and Meneveau2022) for the smooth LaRTE model, likely due to numerical factors such as grid resolution and SFS model approximation (Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2019; Wang, Hu & Zheng Reference Wang, Hu and Zheng2020), rather than to the wall models themselves.

For the homogeneous rough surface, the velocity profiles again display a log-law region. The measured roughness function ( $\Delta U^+ \approx 8.5$ ) is consistent with the Colebrook formula, as well as Nikuradse’s data, as shown in the inset of figure 1. It is noted that this agreement with known results is a verification of the model’s consistency and correct implementation.

Figure 2. Sketch of the computational domain for heterogeneous-wall simulations (not to scale).

5. Results

To evaluate the predictive capabilities of the extended LaRTE model for rough walls in non-equilibrium flows, WMLESs are performed for channels with heterogeneous roughness, and turbulent boundary layers subjected to varying streamwise pressure gradient distributions.

5.1. Heterogeneous roughness

Simulations of flows with heterogeneous roughness were performed to evaluate the LaRTE-RW model’s ability to handle abrupt transitions between rough and smooth surfaces. Specifically, the present simulations were designed to replicate the experimental work conducted by Li et al. (Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019). However, since experimental data were available only for the smooth strip, numerical results for the rough strip were validated against the earlier experimental findings of Antonia & Luxton (Reference Antonia and Luxton1971), which, however, where at a slightly different Reynolds number.

The computational domain dimensions were $L_x \times L_y \times L_z$ = $56 \delta \times \delta \times 7 \delta$ , where $\delta$ represents the channel half-height, corresponding to the experimental data in Li et al. (Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019), recorded at the interface between rough and smooth surfaces. A schematic of the domain is presented in figure 2, reproducing the experimental set-up, where the rough section of length $L_r$ covers $66 \,\%$ of the domain. Here, using the same convention as in Li et al. (Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019), $x_0$ represents the roughness-transition location, and $\hat {x} = x - x_0$ stands for the relative streamwise position. The Reynolds number, defined using the bulk velocity, was ${\textit{Re}}_b = U_b \delta / \nu = 121\,000$ . The simulation employed a grid resolution of $1792 \times 60 \times 224$ points, based on a grid convergence study previously conducted for the same flow configuration by Salomone, Piomelli & De Stefano (Reference Salomone, Piomelli and De Stefano2023). Therein, it was observed that the difference in the mean velocity profiles (at locations $95\,\%$ along the lengths of smooth and rough strips) between the grid used here and a finer grid was less than $3\,\%$ . The resolution in plus units corresponds to $ \Delta x_{s}^+ = \Delta z_{s}^+ = 153$ and $ \Delta y_{s}^+ = 81$ over the smooth strip, and $ \Delta x_{r}^+ = \Delta z_{r}^+ = 218$ and $ \Delta y_{r}^+ = 116$ over the rough strip, where the friction velocity used is evaluated at the smooth and rough reference locations, respectively. The $\Delta U^+$ is imposed through (3.3) where $k_s^+=130$ matches the experiment of Li et al. (Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019) at the reference location over the rough strip, and is $\Delta U^+ = 0$ over the smooth strip. Spatial variations of $\Delta U^+$ are accounted for in the LaRTE-RW evolution (3.18) through the inclusion of the additional terms that involve the spatial derivatives of $\Delta U^+$ .

The profiles of the skin-friction coefficient, normalized by a reference value taken at 95 % of the rough-strip extent, as predicted by the log-law, GMD and LaRTE-RW models – with (LaRTE-RW turbNEQ) and without (LaRTE-RW qEQ) the turbNEQ component – are shown in figure 3. The skin-friction coefficient is defined as

(5.1) \begin{equation} C_{\!f}(x) = \frac {2 \, \tau _w }{\rho _{\infty } U_{cl}^2}, \end{equation}

where $ U_{cl}(x) = U(x, y = \delta )$ is the centreline velocity, obtained by averaging in time and in the spanwise (homogeneous) direction. At the upstream edge of the rough strip, the different models predict distinct values of $ C_{\!f}$ , but they converge to the same value farther downstream. Over the smooth strip, all models predict nearly identical $ C_{\!f}$ profiles, except for the LaRTE-RW model without the turbNEQ contribution, which exhibits a slower relaxation towards the equilibrium $ C_{\!f}$ value.

In figure 4, the skin-friction coefficient $ C_{\!f}$ is normalized by the reference value specific to each surface type (rough or smooth) and compared with experimental data. Both the log-law model (blue solid line) and the GMD model (dashed red line) effectively capture the reestablishment of equilibrium skin friction after the rough patch.

Figure 3. Normalized skin-friction coefficient, $C_{\!f}/C_{fr}$ , where $C_{fr}$ is the skin-friction coefficient measured at $95\,\%$ of the rough strip’s length.

Figure 4. Normalized skin-friction coefficient on the rough strip (a) and smooth strip (b). $C_{fr}$ and $C_{fs}$ are the skin-friction coefficient measured at $95\,\%$ of the rough- and smooth-strip lengths, respectively. Experimental uncertainty from Antonia & Luxton (Reference Antonia and Luxton1971) (a) and Li et al. (Reference Li, de Silva, Rouhi, Baidya, Chung, Marusic and Hutchins2019) (b) is depicted with grey shading.

Over the rough surface following the transition, the log-law and GMD equilibrium models somewhat underpredict the friction coefficient, although the GMD model predicts a slightly higher (more realistic) initial $C_{\!f}$ over the rough region, due to the influence of a weak pressure gradient. The LaRTE-RW (dotted line) and LaRTE-RW including the turbNEQ term (orange dot–dashed line) provide predictions in good agreement with the reference data (black triangles with grey shading illustrating experimental uncertainty).

Notably, in the LaRTE-RW model, roughness-induced discontinuities at the transition are captured, whereas these discontinuities cannot be resolved by the original LaRTE formulation (not shown) that would rely only on an equation like (3.21) that lacks information about spatially changing surface conditions. As the proposed LaRTE-RW model shows better performance on the rough patch compared with the log-law and GMD approaches, the relevance of terms like those in (3.20), accounting for roughness variability, is clearly demonstrated.

Interestingly, in the smooth section, while the log-law, GMD and LaRTE-RW models with the turbNEQ term align well with the experimental data, the LaRTE-RW model without the turbNEQ term predicts an excessively slow recovery. This is attributed to the effects of the relaxation time scale, $T_s$ , which delays the alignment between the local wall stress and the stress at the interface. This discrepancy is resolved by incorporating the turbNEQ stress term $\tau _{w i}^\prime$ , which re-establishes the connection between wall and interface stresses. Indeed, the solution labelled as LaRTE-RW turbNEQ shows a skin friction behaviour that is much closer to that observed in the log-law and GMD models.

The key distinction between the new LaRTE-RW and traditional models lies in the fact that the wall stress derived from solving (3.18) differs from the stress computed at the interface (as in the GMD model, for instance). This is demonstrated in figure 5, where the different components of the stress in the streamwise direction, $\tau _{wx}$ , are illustrated through instantaneous contour maps. Similar to (2.10), the total wall-stress consists of the LaRTE-RW qEQ and the turbNEQ components, $\tau _{wx} = u_\tau u_{\tau x} + \tau _{wx}^\prime$ . The LaRTE-RW component leads to extended regions of higher wall stress at the upstream edge of the rough strip and reduced stress at the leading edge of the smooth strip, an effect that can be attributed to the relaxation time scale $ T_s$ . The turbNEQ contribution, $ \tau _{wx}^\prime$ , primarily affects the wall stress near the transition zones, where it accelerates the local stress adjustment. In the interior of each strip, its influence is minimal. The stress at the interface, $ \tau _{\varDelta x}$ , presents more noise, due to its direct correlation with the unfiltered LES pressure gradient. In contrast, the time scale $ T_s$ acts as a smoothing operator, which explains why both $ u_\tau u_{\tau x}$ and the total wall stress $ \tau _{wx}$ actually appear smoother.

Figure 5. Instantaneous contours of the wall-stress (a–c) and the stress at the wall-model interface (d) with LaRTE-RW model for the streamwise component.

5.2. Streamwise-varying pressure gradient

The WMLESs of rough-wall turbulent boundary layers have been performed, using the experimental dataset of Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022a ) as a reference. In the experiments, two Reynolds numbers were considered, namely ${\textit{Re}} = 2 \times 10^6$ and $3.5 \times 10^6$ , based on the inflow velocity and model chord length, while the angle of attack of the airfoil used for producing the streamwise-varying pressure gradient was varied from $\alpha = -10^\circ$ to $\alpha = 12^\circ$ , with increments of $2^\circ$ .

For the present numerical simulations, the lower Reynolds number and four angles of attack among the ones provided have been selected. The angles of attack chosen are representative of a spectrum of pressure gradients, ranging from weak to moderate pressure-gradient conditions. According to the experimental set-up, the pressure gradient development always begins with an FPG region. This is followed by either a single APG region or a sequence of an APG region followed by another FPG region. Specifically, for the FPG–APG configuration, angles of $\alpha = -10^\circ$ (moderate) and $\alpha = 0^\circ$ (weak) have been chosen, while, for the FPG–APG–FPG configuration, the angles considered are $\alpha = 6^\circ$ and $\alpha = 12^\circ$ . All quantities shown hereafter are averaged in both time and spanwise direction, unless otherwise noted.

5.2.1. Simulations set-up

The computational domain is shown in figure 6, with a spatial extension of $L_x \times L_y \times L_z = 600 \delta _{0}^* \times 45 \delta _{0}^* \times 75 \delta _{0}^*$ , where $\delta _{0}^*$ is the displacement thickness at the inlet plane. The Reynolds number is ${\textit{Re}}_{\delta ^*_{0}} = 17\,000$ , based on the inlet free stream velocity, $U_{0}$ , and $\delta _{0}^*$ . A uniform grid is employed, with $\Delta x/\delta ^*_{0} = \Delta z/\delta ^*_{0} = 0.3472$ , and $\Delta y /\delta ^*_{0} = 0.3461$ , This corresponds to a resolution in plus units of $\Delta x^+ = \Delta z^+ \approx \Delta y^+ = 310$ , with the friction velocity calculated at the reference plane. The wall-model interface is located at $\varDelta /\delta _{0}^* = 1.3$ , allowing four points below the interface, which is within the range recommended by Kawai & Larsson (Reference Kawai and Larsson2012), and the equivalent sand-grain roughness height that is used as input to the models is $k_s/\delta _{0}^* = 0.31$ . The turbulent inflow at the inlet is generated through the recycling and rescaling method of Lund, Wu & Squires (Reference Lund, Wu and Squires1998), while a convective boundary condition (Orlanski Reference Orlanski1976) is used at the outlet.

The reference transverse plane corresponds to the initial point in the ZPG region of the experimental dataset. In our analysis, we use the reference length, $L_{\textit{ref}}$ , to normalize the streamwise coordinates; this length represents the full streamwise extent of the computational domain in the simulations, and corresponds to the rough-wall length of 5.43 m used in the experiment. For consistency, when discussing the results, the origin for the streamwise coordinate corresponds to the reference plane location. Practically, FPG and APG regions are achieved by imposing a spatially varying free stream velocity in the streamwise direction, $U_{\infty }(x)$ , with the wall-normal component, $V_{\infty }(x)$ , being calculated through mass balance. The free stream velocity profile, $ U_{\infty }(x)$ , is prescribed to match experimental data at seven streamwise locations for each of the four angles of attack considered. These profiles are shown in figure 7, along with the corresponding acceleration parameter, defined as

(5.2) \begin{equation} K(x) = \frac {\nu }{U_{\infty }^2} \frac {\textrm{d}U_{\infty }}{\textrm{d}x}, \end{equation}

Figure 6. Computational domain for the varying pressure gradients simulations. The green area is a visual representation of the applied FPG–APG condition.

Figure 7. Prescribed free-stream velocity profiles $ U_{\infty } / U_{\textit{ref}}$ , fitted to match experimental data at seven streamwise locations for each angle of attack. The corresponding acceleration parameter is shown in the insets. Light blue highlights FPG regions, light yellow highlights APG regions. Lines represent simulations; black circles denote experimental data from the Virginia Tech Stability Wind Tunnel (VPI), as reported by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022a ).

Figure 8. Streamwise skin-friction coefficient development for different wall models at each angle of attack, compared with experimental data from the VPI reported by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022a ). The experimental uncertainty is depicted with grey shading.

which is displayed in the insets of each panel. At angles of attack corresponding to FPG–APG conditions ( $ \alpha = -10^\circ$ and $ \alpha = 0^\circ$ ), a region of constant free stream velocity is maintained up to $ x = 0$ , after which $ U_{\infty }$ increases until approximately $ x/L_{\textit{ref}} \simeq 0.35$ , and then decreases. Accordingly, the acceleration parameter is positive in the FPG region, zero at the point of peak velocity and negative in the subsequent APG region.

For FPG–APG–FPG conditions ( $ \alpha = 6^\circ$ and $ \alpha = 12^\circ$ ), the free stream velocity exhibits an initial increase, followed by a decrease and then another increase. The APG region is more extended for $ \alpha = 12^\circ$ , beginning at $ x/L_{\textit{ref}} \simeq 0.23$ and ending at $ x/L_{\textit{ref}} \simeq 0.45$ , compared with $ \alpha = 6^\circ$ , where the APG starts at $ x/L_{\textit{ref}} \simeq 0.31$ and ends at the same downstream location. In both cases, the second FPG region begins at $ x/L_{\textit{ref}} \simeq 0.45$ , and the acceleration parameter changes its sign according to these transitions.

5.2.2. Simulations results

The streamwise development of the skin-friction coefficient, defined here as

(5.3) \begin{equation} C_{\!f}(x) = \frac {2 \, \tau _w }{\rho _{\infty } U_{\infty }^2} , \end{equation}

is shown in figure 8. The overall trend in the wall-model predictions of the skin-friction coefficient is a slight overestimation in the initial FPG regions, but still remaining within the uncertainty range of the experimental data (except for the GMD model that falls above the data range). This is followed by a subsequent underestimation of $C_{\!f}$ from any of the three wall models considered, regardless of whether the downstream region is characterized by a single APG zone (figures 8 a and 8 b), or an opposite APG–FPG sequence (figures 8 c and 8 d).

Making a direct comparison between the present skin-friction coefficient and free stream velocity reveals important history effects in the current flow configurations. As observed in the experimental data, the locations of maximum $ C_{\!f}$ align well with the peak free stream velocity in the FPG–APG cases. For the FPG–APG–FPG configurations, the dip in $ C_{\!f}$ occurs upstream of the corresponding decrease in free stream velocity, consistent with the behaviour of the acceleration parameter shown in the insets of figure 7. In contrast, looking at the WMLES results, we see that the peak $C_{\!f}$ occurs before the maximum velocity in FPG–APG configurations, and after the velocity dip in FPG–APG–FPG cases. This fact indicates that, during FPG phases, the wall models over-accelerate the near-wall boundary layer flow (or causes too weak deceleration) while, during APG phases, the deceleration is more persistent and the flow does not respond fast enough to the subsequent FPG condition that would tend to accelerate the near-wall flow (and thus the wall stress). Moreover, while the $C_{\!f}$ profile in FPG–APG cases follows the expected pattern – rising in FPG regions and falling in APG ones – the behaviour is different in FPG–APG–FPG cases. After the first experimental data point in the APG region, the $C_{\!f}$ does not rise in the subsequent FPG zone, but instead remains mostly flat, as shown in figure 8(c,d), leading to a noticeable underestimation of the skin-friction coefficient.

The discrepancies observed in these comparisons were also present in smooth-wall simulations (not shown) performed using the same computational domain and grid resolution as the rough-wall cases. This confirms that the inaccuracies in the rough-wall results are not due to the inclusion of roughness effects in the wall models, but are instead intrinsic to the wall-modelling methodologies themselves.

The LaRTE-RW model, overall, performs marginally better than other wall models, avoiding the overprediction of the skin-friction coefficient in the initial stages of pressure-gradient development. The GMD model tends to either overestimate or underestimate the skin-friction coefficient in comparison with other models. This discrepancy is most pronounced at $\alpha =-10^\circ$ (see figure 8 a), where the GMD model’s prediction falls outside the uncertainty bounds, even during the early phase of pressure gradient development. This error can be attributed to the exclusion of the convective term in the model development, when explicitly accounting for the pressure gradient. In fact, a non-zero streamwise pressure gradient is invariably accompanied by flow acceleration or deceleration, which should be captured by the convective term. This aligns with the findings of Hickel et al. (Reference Hickel, Touber, Bodart and Larsson2013), who argued that incorporating the pressure gradient effects without including the convective (unsteady) terms creates an inconsistency with the momentum balance, particularly outside the boundary layer, where the Bernoulli principle applies. Notably, when the LaRTE-RW model is coupled with the GMD model at the interface, these overpredictions and underpredictions are mitigated. This is because the LaRTE-RW model solves a partial differential equation for the wall stress that explicitly includes the convective term.

It must be noted that the wall models considered in this work are all strongly coupled to the flow dynamics at the wall-model interface, which may result in the overlooking of important behaviour occurring below this interface. For example, Parthasarathy & Saxton-Fox (Reference Parthasarathy and Saxton-Fox2023) showed that when an APG region follows an upstream FPG one, a bimodal structure in the Reynolds stresses can emerge, with an intensification of the near-wall peak across the APG region. Here, the near-wall resolution is likely to be too coarse to capture the increase in Reynolds stresses close to the wall, which would otherwise mitigate the drop in skin-friction coefficient. Instead, the wall models interpret the reduced velocity in the APG region as a decrease in wall stress, without accounting for the additional Reynolds stress that would be present in a real flow. As a result, the predicted $ C_{\!f}$ values in the APG region of figures 8(a) and 8(b) are lower than those observed in experiments.

On the other hand, in the transition from APG to FPG, the flow accelerates and the near-wall velocity gradient steepens. In cases of strong FPG, this could lead to relaminarization, reduced Reynolds stress and consequently lower $ C_{\!f}$ . However, in the present case, the FPG is not strong enough to relaminarize the flow. Instead, it compresses the turbulent structures towards the wall and increases the near-wall shear. This leads to a rapid rise in skin friction, which is often underpredicted by wall models that lack the resolution or physics to capture the turbulent memory and re-energization process associated with such transitions. Moreover, the lag in the response of $ C_{\!f}$ in figures 8(c) and 8(d) is due to the wall models being coupled to the outer-layer dynamics through the wall-model interface, which adjusts more slowly to the imposed forcing than the near-wall turbulence would in a real-world scenario.

Figure 9. Outer-units velocity profiles for each angle of attack. From bottom to top, the seven profiles correspond to experimental locations at $x / L_{\textit{ref}}=0, 0.11, 0.22, 0.34, 0.45, 0.56$ and $0.67$ . Each velocity profile, except the first one, is shifted upward by 0.5 units for clarity. Black circles represent experimental data from the VPI, as reported by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022 a). The dashed line represents the location of the wall-model interface.

Figure 9 shows the velocity profiles in outer units, for the four different flow configurations. At the first experimental location ( $x/L_{\textit{ref}}=0$ ), numerical and experimental results agree well, indicating that the incoming boundary layer flow is consistent between the two approaches. However, moving downstream, some differences between experiment and simulations become apparent. All the models overestimate the velocity deficit following the reversal of the pressure gradient at the location of the wall-model interface. This velocity deficit, introduced at the wall-model interface, is reflected in the underestimation of $C_{\!f}$ at the corresponding downstream locations. Additionally, in areas where the pressure gradient distorts the velocity profile (near the peak/dip of the free stream velocity), the simulations predict a profile that maintains a log-law shape. The distortion of the profiles begins to appear towards the end of the domain, and is particularly evident at the last location ( $x/L_{\textit{ref}}=0.67$ ) for the strongest FPG–APG case ( $\alpha = -10^\circ$ ).

5.2.3. Effect of turbNEQ portion on the LaRTE-RW model

Simulations using the LaRTE-RW model without the turbNEQ component, referred to as LaRTE-RW qEQ, were also conducted. The skin-friction coefficient results, shown in figures 10(a) and 10(b) for angles of attack $\alpha = -10^\circ$ and $\alpha = 12^\circ$ , reveal that in the FPG–APG case, the LaRTE-RW qEQ provides improved predictions, although the peak $C_{\!f}$ is still under-predicted. However, in the FPG–APG–FPG configuration, the same model struggles to accurately capture the $C_{\!f}$ development.

This discrepancy is due to the time scale $T_s$ in the LaRTE model, which introduces a delay in how the wall shear-stress, $\tau _{w i}$ , adjusts to the qEQ stress at the wall-interface $\overline {\tau }_{\varDelta i}$ . Reasonably, this delay is advantageous in cases where models respond too quickly to free stream disturbances (as in FPG–APG), yet it proves detrimental in situations requiring a faster response (as in FPG–APG–FPG). In essence, an effective wall model should have the flexibility to either maintain a delayed response or, alternatively, adapt more rapidly based on free stream conditions.

The LaRTE-RW model is constrained by its qEQ framework, even when the turbNEQ portion is added, limiting it to approach, but not exceed, the response of the EQWM. This limitation is especially apparent in situations where a faster-than-equilibrium response would be preferable. Indeed, developing an approach that allows the wall model to ‘overreact’ (responding more aggressively than the qEQ solution) could significantly enhance prediction accuracy but it is far from clear how to implement such a response mechanism in a wall model.

5.2.4. Effects of grid refinement and SFS model

To investigate whether grid refinement could enhance the accuracy of skin-friction coefficient predictions downstream of the APG region, further computations were conducted using anisotropically refined grids. Specifically, three different simulations were performed by doubling the grid points in one space direction, while leaving the resolution in the other two directions unchanged. The WMLES solutions for the flow configuration FPG–APG–FPG corresponding to $\alpha =12^\circ$ were obtained using the log-law EQWM. By inspection of figure 11(a), showing $C_{\!f}$ profiles, as minor differences appear at downstream locations, grid refinement is demonstrated not to significantly impact the solution.

Figure 10. Streamwise development of the skin-friction coefficient $ C_{\!f}$ for two representative angles of attack: (a) FPG/APG case $ \alpha = -10^\circ$ , and (b) FPG/APG/FPG case $ \alpha = 12^\circ$ . Results are shown with the turbNEQ contribution included (dot–dashed line) and excluded (dotted line). Black circles represent experimental data from the VPI, as reported by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022 a). Experimental uncertainty is depicted with grey shading.

Figure 11. (a) Streamwise development of the skin-friction coefficient $C_{\!f}$ , with anisotropic grid refinement. (b) Comparison of Vreman SFS and LDSM for log-law EQWM, and LaRTE-RW at angle of attack $\alpha = 12^\circ$ . Black circles represent experimental data from the VPI, as reported by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Lowe and Devenport2022a ). Experimental uncertainty is depicted with grey shading.

Furthermore, the influence of the SFS model was assessed performing additional simulations employing the Lagrangian dynamic SFS model (LDSM) proposed by Meneveau, Lund & Cabot (Reference Meneveau, Lund and Cabot1996). These solutions were obtained for the same FPG–APG–FPG case, following either log-law EQWM or LaRTE-RW approach. The predicted $C_{\!f}$ profiles are shown in figure 11(b). The use of LDSM mitigates the initial overprediction of $C_{\!f}$ in the FPG region. However, after the pressure gradient reversal, the solution underpredicts this parameter, and none of the the different SFS models enhance the prediction of $C_{\!f}$ after the pressure gradient reversal, where the predicted values remain markedly lower than the experimental ones.

Very interestingly, despite the top-wall boundary condition being unmodified and, thus, identical to that shown in figure 7(d) – indicating that the flow locally experiences the same imposed pressure gradient – different SFS models result in distinct $C_{\!f}$ developments, particularly at downstream regions. This fact highlights the influence of upstream flow history on the evolution of the turbulent boundary layer, suggesting that it is not solely governed by local levels of free stream velocity and pressure gradient. Consequently, differences in predicting the boundary layer behaviour at upstream locations (through different eddy viscosity models, for instance) can propagate downstream, leading to observed deviations in WMLES predictions.

6. Summary and conclusions

This paper introduces a significant extension of the LaRTE model for rough walls, enabling seamless prediction of flow behaviour across smooth and fully rough regimes. The extended model’s performance was evaluated against traditional models, including the log-law EQWM and the GMD model, under various conditions, including homogeneous and heterogeneous roughness, and varying pressure gradients. The model was successfully applied to WMLES of flow over patches with smooth-to-rough and rough-to-smooth transitions, showing improved accuracy compared with the traditional log-law wall model. In particular, the case of rough to smooth transition highlighted the importance of including a turbNEQ component into the wall model to accelerate convergence towards the smooth wall asymptote.

Furthermore, the study of pressure gradient cases revealed important insights into wall-model responses to free stream forcing. While wall models accurately capture the general trends of increasing friction in FPG regions, and decreasing friction in APG regions, inaccuracies emerge downstream of pressure-gradient reversals. Notably, in FPG–APG cases, wall models tend to react too quickly to free stream disturbances, whereas in FPG–APG–FPG cases, the response is overly delayed. Transitions between FPG and APG regions further revealed limitations in the current modelling framework. In particular, when the flow transitions from FPG to APG, wall models tend to underpredict the skin-friction coefficient due to their inability to capture the re-energization of near-wall turbulence. Conversely, in APG to FPG transitions, experimental observations show a rapid recovery in wall shear stress, attributed to the realignment of turbulent structures with the wall. However, the LaRTE-RW model, like other WMLES approaches, underpredicts this response due to its dependence on the slowly evolving outer-layer signal. These findings reinforce the need for wall models that incorporate both non-equilibrium dynamics and a more faithful representation of near-wall turbulence memory and response times. Therefore, future work will focus on further refining the LaRTE-RW model with enhanced non-equilibrium components to address this particular need.

Findings from smooth-wall simulations add further support to these observations, suggesting that the present inaccuracies observed in rough-wall configurations are not caused by the incorporation of roughness effects, but are instead rooted in limitations of the underlying wall-modelling methodologies themselves. This study also highlights the impact of SFS models on $C_{\!f}$ predictions, particularly highlighting that inaccuracies in upstream boundary layer predictions can propagate downstream, underscoring the importance of accurately resolving boundary layer states at early stages.

However, it is worth noting that a key advantage of the newly proposed formulation is in the proven ability to separate the qEQ component of wall stress from its non-equilibrium counterpart. This modular structure facilitates targeted refinement and tuning. In contrast, models like the log-law EQWM and GMD integrate these dynamics inseparably, making such improvements more difficult to implement. We note that although the assumption of the mixing length for deriving the mean velocity profile from the simplified boundary layer momentum equation may not adequately apply in the vicinity of rough-to-smooth transitions or under varying pressure gradients and it represents a drawback in the log-law and GMD models, in the LARTE-RW context is logical, as it aligns with the qEQ conditions the LaRTE-RW equation is designed to represent. Departures from equilibrium effects must be considered through the non-equilibrium portions of the model.

Overall, the current results emphasize the benefits of the new LaRTE-RW model in WMLES for rough-wall flows under non-equilibrium conditions. The model demonstrates superior predictive capabilities, particularly during the initial stages of FPG development, also when paired with the LDSM. The over-predictions of wall stress observed in the GMD model are effectively mitigated by the LaRTE-RW model, leading to slightly improved performance. The LaRTE-RW introduces advantageous time delays where standard EQWMs provide excessively fast responses to imposed flow conditions. However, the time delays prove detrimental in situations requiring more swift responses, such as for boundary layers developing below accelerating–decelerating–accelerating free stream conditions. Finally, the LaRTE-RW model surpasses traditional models like the log-law EQWM in handling complex configurations such as heterogeneous roughness transitions. The present research findings demonstrate the LaRTE-RW model’s potential for advancing the predictive accuracy of WMLES in complex turbulent flow configurations, but that additional research is needed to fine-tune response times and mechanisms to further improve wall modelling in LES.

Funding

The authors acknowledge the financial support by the Natural Sciences and Engineering Research Council of Canada (NSERC-CRNSG) under the Discovery Grant programme. This research was enabled in part by computational support provided by the Digital Research Alliance of Canada. C.M.’s wall modelling research activities are supported by the Office of Scientific Research (grant no. N00014-2312185)

Declaration of interests

The authors report no conflict of interest.

Appendix A. Analytical derivation of the advection integral

In this appendix, we provide a summary of the integration of the advection term for the proposed LaRTE-RW model. Starting with the continuity equation,

(A1) \begin{equation} \partial _i \overline {u}_{\textit{si}} + \partial _y \overline {v} = 0, \end{equation}

we derive the wall-normal velocity component, $\overline {v}(\varDelta )$ , required for the advection term,

(A2) \begin{equation} \overline {v}(\varDelta ) = - \int _0^\varDelta \partial _i \overline {u}_{\textit{si}} \, \textrm{d}y = -\partial _i u_{\tau i} \int _0^\varDelta F \, \textrm{d}y - u_{\tau i} \int _0^\varDelta \partial _i F \, \textrm{d}y. \end{equation}

Following a procedure similar to that used by Fowler et al. (Reference Fowler, Zaki and Meneveau2023), the previous expression becomes

(A3) \begin{eqnarray} \overline {v}(\varDelta ) &=& - (\partial _s u_\tau ) \, \varDelta F_\varDelta - u_\tau (\partial _i s_i) \int _0^\varDelta F \, \textrm{d}y + u_{\tau i} \, (\partial _i \Delta U^+) \int _0^\varDelta f_{\textit{NW}} \, \textrm{d}y, \end{eqnarray}

where $\partial _s = s_i \partial _i$ represents differentiation in the direction of $s_i$ . This expression follows the form of the corresponding equation in Fowler et al. (Reference Fowler, Zaki and Meneveau2023), except with the function $F(y^+, k_s^+)$ instead of $f(y^+)$ . The last term stems from the spatial derivatives of $f_{\textit{NW}} \Delta U^+$ .

The advection term is given by

(A4) \begin{equation} \int _0^\varDelta \left [ \partial _{\!j} (\overline {u}_{\textit{si}} \overline {u}_{\textit{sj}}) + \partial _y (\overline {v} \, \overline {u}_{\textit{si}}) \right ] \textrm{d}y = \int _0^\varDelta \partial _{\!j} (\overline {u}_{\textit{si}} \overline {u}_{\textit{sj}}) \textrm{d}y + \overline {v}(\varDelta ) u_{\tau } s_i F_\varDelta . \end{equation}

Focusing on the first integral, it holds that

(A5) \begin{eqnarray} \int _0^\varDelta \partial _{\!j} (\overline {u}_{\textit{si}} \overline {u}_{\textit{sj}}) \textrm{d}y &=& u_{\tau i} u_{\tau j} \int _0^\varDelta \partial _{\!j} F^2 \, \textrm{d}y + \partial _{\!j} (u_{\tau i} u_{\tau j}) \int _0^\varDelta F^2 \, \textrm{d}y, \end{eqnarray}

where the function $F^2(y^+,k_s^+)$ is defined as follows:

(A6) \begin{eqnarray} F^2(y^+, k_s^+) &=& f^2(y^+) + f^2_{\textit{NW}}(y^+) \Delta U^{+2}(k_s^+) - 2 f(y^+) f_{\textit{NW}}(y^+) \Delta U^+(k_s^+). \end{eqnarray}

Equation (A5) can be expanded with the inclusion of the derivatives of the two functions $f^2_{\textit{NW}}(y^+) \Delta U^{+2}(k_s^+)$ and $- 2 f(y^+) f_{\textit{NW}}(y^+) \Delta U^+(k_s^+)$ , resulting in

(A7) \begin{align} \int _0^\varDelta \partial _{\!j} (\overline {u}_{\textit{si}} \overline {u}_{\textit{sj}}) \textrm{d}y &= u_{\tau i} (\partial _s u_\tau ) \Delta F^2_\varDelta + [u_\tau (\partial _s u_{\tau i})+ u_{\tau }^2 s_i (\partial _{\!j} s_{\!j})] \int _{0}^\varDelta F^2 \, \textrm{d}y \nonumber \\ &\quad - 2 \, u_{\tau i} u_{\tau j} \, (\partial _{\!j} \Delta U^+) f_{\textit{NW},\varDelta } \, f_\varDelta (\varDelta - \delta ^*_{\varDelta S} - \theta _{\Delta f \, f_{\textit{NW}}})\nonumber \\ &\quad + 2 \, u_{\tau i} u_{\tau j} \, (\partial _{\!j} \Delta U^+) f^2_{\textit{NW},\varDelta } \, \Delta U^+ (\varDelta - \delta ^*_{\varDelta , \textit{NW}} - \theta _{\varDelta , \textit{NW}}). \end{align}

On the right-hand side of this equation, the first two terms correspond to analogous terms defined by Fowler et al. (Reference Fowler, Zaki and Meneveau2023), while the last two terms arise from considering the additional functions $f^2_{\textit{NW}}(y^+)\Delta U^{+2}(k_s^+)$ and $-2f(y^+)\, f_{\textit{NW}}(y^+) \, \Delta U^{+}(k_s^+)$ , respectively.

By combining (A7) with the term $ \overline {v}(\varDelta ) u_{\tau i} F_\varDelta$ , the full expression for the advection term becomes

(A8) \begin{align} \int _0^\varDelta \left [ \partial _{\!j} (\overline {u}_{\textit{si}} \overline {u}_{\textit{sj}}) + \partial _y (\overline {v} \, \overline {u}_{\textit{si}}) \right ] \textrm{d}y &= u_{\tau j} (\partial _{\!j} u_{\tau i}) F^2_\varDelta (\varDelta - \delta ^*_\varDelta - \theta _\varDelta ) - \mathcal{F}_{1 i} F_\varDelta \theta _\varDelta \nonumber \\ &\quad - 2 \, \mathcal{F}_{2 i} \, f_\varDelta (\varDelta - \delta ^*_{\varDelta S} - \theta _{\varDelta f \, f_{\textit{NW}}}) \nonumber \\ &\quad + 2 \, \mathcal{F}_{2 i} \, f_{\textit{NW},\varDelta } \Delta U^+ (\varDelta - \delta ^*_{\varDelta , \textit{NW}} - \theta _{\varDelta , \textit{NW}}) \nonumber \\ &\quad + \mathcal{F}_{2 i} \, F_\varDelta (\varDelta - \delta ^*_{\varDelta , \textit{NW}}), \end{align}

where $ \mathcal{F}_{1 i}$ and $ \mathcal{F}_{2 i}$ are defined as

(A9) \begin{align} \mathcal{F}_{1 i} &= u_\tau u_{\tau i} \, (\partial _{\!j} s_{\!j}) \, F_\varDelta , \end{align}
(A10) \begin{align} \mathcal{F}_{2 i} &= u_{\tau i} u_{\tau j} \, (\partial _{\!j} \Delta U^+) \, f_{\textit{NW},\varDelta }. \end{align}

This leads to the final expression for the advective velocity $ V_{\tau j}$ as defined in (3.16) of the main text, and consequently to (3.18).

Appendix B. Inner-similarity and near-wall functions in the LaRTE-RW model

In this appendix, we present the functions used to describe the velocity profile in the proposed LaRTE-RW model. The function $ F(y^+, k_s^+)$ represents the inner-similarity function that accounts for the combined effect of wall-normal distance and surface roughness. It is constructed using two different components, specifically, the smooth-wall inner-similarity function $ f(y^+)$ , and the near-wall roughness correction function $ f_{\textit{NW}}(y^+)$ , combined with the roughness function $ \Delta U^+(k_s^+)$ . The function $ f(y^+)$ follows the form of the equivalent function introduced in Fowler et al. (Reference Fowler, Zaki and Meneveau2022),

(B1) \begin{equation} f(y^+) = \left [ \frac {1}{\kappa } \log (\kappa _2 + y^+) + B \right ] \left [ 1 + \left ( \frac {y^+}{\kappa _1} \right )^{-\beta } \right ]^{-\frac {1}{\beta }}, \end{equation}

while the full expression for $F(y^+, k_s^+)$ is given by

(B2) \begin{equation} F(y^+, k_s^+) = f(y^+) - f_{\textit{NW}}(y^+) \Delta U^+(k_s^+), \end{equation}

where

(B3) \begin{equation} f_{\textit{NW}}(y^+) = \left [ 1 + \left ( \frac {y^+}{\kappa _1} \right )^\beta \right ]^{-\frac {1}{\beta }}. \end{equation}

In the above equations, the constants are defined as

(B4) \begin{align} B = 4.95, \quad \kappa _2 = 9.753, \quad \beta = 1.903, \quad \kappa _1 = \frac {1}{\kappa } \log (\kappa _2) + B, \end{align}

while $ \kappa = 0.4$ is the von Kármán constant. Basically, the function $ f(y^+)$ describes the inner-similarity for smooth walls, while $ f_{\textit{NW}}(y^+)$ introduces the correction for the near-wall behaviour in rough-wall conditions, ensuring that the roughness contribution, represented by $ \Delta U^+(k_s^+)$ , diminishes near the wall.

The resulting function $F(y^+, k_s^+)$ shifts downward as $ k_s^+$ increases, reflecting the increased drag imposed by the rough surface. In the logarithmic region, the term $ f_{\textit{NW}}(y^+) \Delta U^+(k_s^+)$ reaches the values of $ \Delta U^+(k_s^+)$ predicted by the Colebrook formula, and approaches zero in the viscous sublayer. Figure 12 shows the plotted functions for the inner-similarity and near-wall corrections across different values of $ k_s^+$ . Figure 12(a) depicts the overall behaviour of the inner-similarity function, while figure 12(b) highlights the near-wall correction.

Figure 12. Inner-similarity function (a) and near-wall function (b) in the LaRTE-RW model.

Appendix C. Evaluation of the friction factor $\boldsymbol{c}_{\!\boldsymbol{f}}^{\boldsymbol{w}\boldsymbol{m}}$ for the interface stress

The baseline wall stress model used in this work follows the empirical fitting procedure proposed by Meneveau (Reference Meneveau2020), which provides a unified framework to evaluate the friction Reynolds number ${\textit{Re}}_{\tau \varDelta }^{\textit{uf}}$ . The model avoids numerical iteration by providing an explicit expression based on dimensionless quantities available from LES input.

We begin by evaluating the non-dimensional parameter

(C1) \begin{align} \chi = \frac {N \varDelta }{U^2_{\textit{LES}}}\min \left [ \frac {{\textit{Re}}_{\varDelta }}{{\textit{Re}}_{\tau \varDelta }^{\textit{fit}}}, \; \frac {1}{\kappa } \log \left (\frac {\exp (\kappa \, 8.5) \varDelta }{k_s}\right ) \right ]^2, \end{align}

with $N = \rho ^{-1} \partial \hat {p}_{\textit{LES}} / \partial x$ . To ensure the validity of the fits and derivations, $\chi$ should remain within the range $-0.2 \leqslant \chi \leqslant 0.2$ . In practice, $\chi = \textrm{sign}(\chi ) \boldsymbol{\cdot }\min (|\chi |, 0.2)$ . The procedure consists of the following three stages.

  1. (i) To compute the baseline fit for ZPG,

    (C2) \begin{align} \begin{aligned} \beta _1 &= \left (1 + 0.155 / {\textit{Re}}_\varDelta ^{-0.03} \right )^{-1}, \quad \beta _2 = 1.7 - \left (1 + 36 / {\textit{Re}}_\varDelta ^{0.75} \right )^{-1}, \\ \quad \kappa &= 0.4, \quad \kappa _3 = 0.005, \quad \kappa _4 = \kappa _3^{\beta _1 - 1/2}, \\ {\textit{Re}}_{\tau \varDelta }^{\textit{fit}} &= \kappa _4 {\textit{Re}}_\varDelta ^{\beta _1} \left [ 1 + (\kappa _3 {\textit{Re}}_\varDelta )^{-\beta _2} \right ]^{\frac {\beta _1 - 1/2}{\beta _2}}. \end{aligned} \end{align}
  2. (ii) To include mild pressure gradients, evaluating the composite friction Reynolds number using the following set of equations:

    (C3) \begin{align} {\textit{Re}}_{\tau \varDelta ,\nu } &= \left ( 1 + 0.5 \chi \right )^{-1/2} {\textit{Re}}_{\tau \varDelta }^{\textit{fit}}, \nonumber \\ {\textit{Re}}_{\varDelta }^* &= {\textit{Re}}_{\varDelta } - \frac {\chi }{2 \kappa } {\textit{Re}}_{\tau \varDelta }^{\textit{fit}} \left ( 1 - \frac {11}{{\textit{Re}}_{\tau \varDelta }^{\textit{fit}}} \right ) \left [ 1 + \left ( \frac {50}{{\textit{Re}}_{\tau \varDelta }^{\textit{fit}}} \right )^2 \right ]^{-1/2}, \nonumber \\ \beta _1^* &= \left (1 + 0.155 / {\textit{Re}}_\varDelta ^{* \,0.03} \right )^{-1}, \quad \beta _2^* = 1.7 - \left (1 + 36 / {\textit{Re}}_\varDelta ^{* \, 0.75} \right )^{-1}, \nonumber \\ \kappa _4^* &= \kappa _3^{\beta _1^* - 1/2}, \nonumber \\ {\textit{Re}}_{\tau \varDelta ,{in}} &= \kappa _4^* ({\textit{Re}}_{\varDelta }^*)^{\beta _1^*} \left [ 1 + (\kappa _3 {\textit{Re}}_{\varDelta }^*)^{-\beta _2^*} \right ]^{(\beta _1^* - 1/2)/\beta _2^*}, \nonumber \\ \theta &= \left ( 1 + {\textit{Re}}_\varDelta /400 \right )^{-1}, \nonumber \\ {\textit{Re}}_{\tau \varDelta }^{\textit{com}} &= \theta \, {\textit{Re}}_{\tau \varDelta ,\nu } + (1 - \theta ) \, {\textit{Re}}_{\tau \varDelta ,{in}}. \end{align}
  3. (iii) To merge with rough-wall representations, defining the asymptotic fully rough friction Reynolds number as

    (C4) \begin{equation} {\textit{Re}}_{\tau \varDelta }^{\infty } = {\textit{Re}}_\varDelta \left [ {\dfrac {1}{\kappa } \log \left (\frac {\exp (\kappa \, 8.5) \varDelta }{k_s} \right ) + \frac {\chi }{2 \kappa }\left ( 1 - \frac {k_s}{\exp (\kappa \, 8.5) \varDelta } \right )} \right ]^{-1}. \end{equation}
    The unified friction Reynolds number, accounting for viscous, inertial, pressure-gradient and roughness effects, is then obtained as
    (C5) \begin{equation} {\textit{Re}}_{\tau \varDelta }^{\textit{uf}} = \left [ \left ( {\textit{Re}}_{\tau \varDelta }^{\textit{com}} \right )^6 + \left ( {\textit{Re}}_{\tau \varDelta }^{\infty } \right )^6 \right ]^{1/6}. \end{equation}

Based on this equation, both the modelled friction velocity $u_\tau ^{w\textit{m}} = {\textit{Re}}_{\tau \varDelta }^{\textit{uf}} \nu / \varDelta$ and the EQWM friction factor $ c_{\!f}^{w\textit{m}} = 2 ( {{\textit{Re}}_{\tau \varDelta }^{\textit{uf}}}/{{\textit{Re}}_{\varDelta }} )^2$ can be evaluated. Finally, the modelled total stress at the model interface is obtained from (3.22).

References

Ambrogi, F., Piomelli, U. & Rival, D.E. 2022 Characterization of unsteady separation in a turbulent boundary layer: mean and phase-averaged flow. J. Fluid Mech. 945, A10.10.1017/jfm.2022.561CrossRefGoogle Scholar
Antonia, R.A. & Luxton, R.E. 1971 The response of a turbulent boundary layer to a step change in surface roughness Part 1. Smooth to rough. J. Fluid Mech. 48 (4), 721761.10.1017/S0022112071001824CrossRefGoogle Scholar
Antonia, R.A. & Luxton, R.E. 1972 The response of a turbulent boundary layer to a step change in surface roughness. Part 2. Rough-to-smooth. J. Fluid Mech. 53 (4), 737757.10.1017/S002211207200045XCrossRefGoogle Scholar
Aubertine, C.D., Eaton, J.K. & Song, S. 2004 Parameters controlling roughness effects in a separating boundary layer. Intl J. Heat Fluid Flow 25 (3), 444450.10.1016/j.ijheatfluidflow.2004.02.007CrossRefGoogle Scholar
Bae, H.J., Lozano-Durán, A., Bose, S.T. & Moin, P. 2019 Dynamic slip wall model for large-eddy simulation. J. Fluid Mech. 859, 400432.10.1017/jfm.2018.838CrossRefGoogle ScholarPubMed
Bobke, A., Vinuesa, R., Örlü, R. & Schlatter, P. 2017 History effects and near equilibrium in adverse-pressure-gradient turbulent boundary layers. J. Fluid Mech. 820, 667692.10.1017/jfm.2017.236CrossRefGoogle Scholar
Bose, S.T. & Park, G.I. 2018 Wall-modeled large-eddy simulation for complex turbulent flows. Annu. Rev. Fluid Mech. 50 (1), 535561.10.1146/annurev-fluid-122316-045241CrossRefGoogle ScholarPubMed
Bou-Zeid, E., Meneveau, C. & Parlange, M.B. 2004 Large-eddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: blending height and effective surface roughness. Water Resour. Res. 40, W02505.10.1029/2003WR002475CrossRefGoogle Scholar
Bradshaw, P. & Ferriss, D.H. 1965 The response of a retarded equilibrium turbulent boundary layer to the sudden removal of pressure gradient. National Physical Laboratory, Aerodynamics Division.Google Scholar
Cal, R.B., Brzek, B., Johansson, T.G. & Castillo, L. 2009 The rough favourable pressure gradient turbulent boundary layer. J. Fluid Mech. 641, 129155.10.1017/S0022112009991352CrossRefGoogle Scholar
Chamorro, L.P. & Porté-Agel, F. 2009 Velocity and surface shear stress distributions behind a rough-to-smooth surface transition: a simple new model. Boundary-Layer Meteorol. 130, 2941.10.1007/s10546-008-9330-xCrossRefGoogle Scholar
Cheng, H. & Castro, I.P. 2002 Near-wall flow development after a step change in surface roughness. Boundary-Layer Meteorol. 105 (3), 411432.10.1023/A:1020355306788CrossRefGoogle Scholar
Chorin, A.J. 1968 Numerical solution of the Navier–Stokes equations. Maths Comput. 22 (104), 745762.10.1090/S0025-5718-1968-0242392-2CrossRefGoogle Scholar
Chung, D., Hutchins, N., Schultz, M.P. & Flack, K.A. 2021 Predicting the drag of rough surfaces. Annu. Rev. Fluid Mech. 53 (2021), 439471.10.1146/annurev-fluid-062520-115127CrossRefGoogle Scholar
Chung, D. & Pullin, D.I. 2009 Large-eddy simulation and wall modelling of turbulent channel flow. J. Fluid Mech. 631, 281309.10.1017/S0022112009006867CrossRefGoogle Scholar
Clauser, F.H. 1954 Turbulent boundary layers in adverse pressure gradients. J. Aeronaut. Sci. 21 (2), 91108.10.2514/8.2938CrossRefGoogle Scholar
Colebrook, C.F. 1939 Turbulent flow in pipes with particular reference to the transition region between smooth- and rough-pipe laws. J. Inst. Civil Engrs Lond. 11, 133156.10.1680/ijoti.1939.13150CrossRefGoogle Scholar
Coleman, H.W., Moffat, R.J. & Kays, W.M. 1977 The accelerated fully rough turbulent boundary layer. J. Fluid Mech. 82 (3), 507528.10.1017/S0022112077000810CrossRefGoogle Scholar
Deardorff, J.W. 1970 A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J. Fluid Mech. 41 (2), 453480.10.1017/S0022112070000691CrossRefGoogle Scholar
Devenport, W.J. & Lowe, K.T. 2022 Equilibrium and non-equilibrium turbulent boundary layers. Prog. Aerosp. Sci. 131, 100807.10.1016/j.paerosci.2022.100807CrossRefGoogle Scholar
Finnigan, J. 2000 Turbulence in plant canopies. Annu. Rev. Fluid Mech. 32 (1), 519571.10.1146/annurev.fluid.32.1.519CrossRefGoogle Scholar
Flack, K.A. & Schultz, M.P. 2014 Roughness effects on wall-bounded turbulent flows. Phys. Fluids 26 (10), 101305.10.1063/1.4896280CrossRefGoogle Scholar
Flack, K.A., Schultz, M.P. & Shapiro, T.A. 2005 Experimental support for Townsend’s Reynolds number similarity hypothesis on rough walls. Phys. Fluids 17 (3), 035102.10.1063/1.1843135CrossRefGoogle Scholar
Fowler, M., Zaki, T.A. & Meneveau, C. 2022 A Lagrangian relaxation towards equilibrium wall model for large eddy simulation. J. Fluid Mech. 934 (A44), 137.10.1017/jfm.2021.1156CrossRefGoogle Scholar
Fowler, M., Zaki, T.A. & Meneveau, C. 2023 A multi-time-scale wall model for large-eddy simulations and applications to non-equilibrium channel flows. J. Fluid Mech. 974, A51.10.1017/jfm.2023.585CrossRefGoogle Scholar
Fritsch, D.J., Vishwanathan, V., Roy, C.J., Lowe, K.T. & Devenport, W.J. 2022 a Turbulence and pressure fluctuations in rough wall boundary layers in pressure gradients. Exp. Fluids 63, 140.10.1007/s00348-022-03476-9CrossRefGoogle Scholar
Fritsch, D.J., Vishwanathan, V., Todd Lowe, K. & Devenport, W.J. 2022 b Fluctuating pressure beneath smooth wall boundary layers in nonequilibrium pressure gradients. AIAA J. 60 (8), 47254743.10.2514/1.J061431CrossRefGoogle Scholar
Garratt, J.R. 1990 The internal boundary layer — a review. Boundary-Layer Meteorol. 50, 171203.10.1007/BF00120524CrossRefGoogle Scholar
Hanson, R.E. & Ganapathisubramani, B. 2016 Development of turbulent boundary layers past a step change in wall roughness. J. Fluid Mech. 795, 494523.10.1017/jfm.2016.213CrossRefGoogle Scholar
Hickel, S., Touber, E., Bodart, J. & Larsson, J. 2013 A parametrized non-equilibrium wall-model for large-eddy simulations. In TSFP-8. Begel House Inc.Google Scholar
Jiménez, J. 2004 Turbulent flows over rough walls. Annu. Rev. Fluid Mech. 36, 173196.10.1146/annurev.fluid.36.050802.122103CrossRefGoogle Scholar
Jones, M.B., Marusic, I. & Perry, A.E. 2001 Evolution and structure of sink-flow turbulent boundary layers. J. Fluid Mech. 428, 127.10.1017/S0022112000002597CrossRefGoogle Scholar
Kawai, S. & Larsson, J. 2012 Wall-modeling in large eddy simulation: length scales, grid resolution, and accuracy. Phys. Fluids 24 (1), 015105.10.1063/1.3678331CrossRefGoogle Scholar
Keating, A., Piomelli, U., Bremhorst, K. & Nešić, S. 2004 Large-eddy simulation of heat transfer downstream of a backward-facing step. J. Turbul. 5 (1), 020.10.1088/1468-5248/5/1/020CrossRefGoogle Scholar
Kim, J. & Moin, P. 1985 Application of a fractional-step method to incompressible Navier–Stokes equations. J. Comput. Phys. 59 (2), 308323.10.1016/0021-9991(85)90148-2CrossRefGoogle Scholar
Krogstad, P.-Å., Antonia, R.A. & Browne, L.W.B. 1992 Comparison between rough-and smooth-wall turbulent boundary layers. J. Fluid Mech. 245, 599617.10.1017/S0022112092000594CrossRefGoogle Scholar
Larsson, J., Kawai, S., Bodart, J. & Bermejo-Moreno, I. 2016 Large eddy simulation with modeled wall-stress: recent progress and future directions. Mech. Engng Rev. 3 (1), 123.Google Scholar
Lee, M. & Moser, R.D. 2015 Direct numerica simulation of turbulent channel flow up to ${\textit{Re}}_\tau \approx 5200$ . J. Fluid Mech. 774, 395415.10.1017/jfm.2015.268CrossRefGoogle Scholar
Li, M., de Silva, C.M., Rouhi, A., Baidya, R., Chung, D., Marusic, I. & Hutchins, N. 2019 Recovery of wall-shear stress to equilibrium flow conditions after a rough-to-smooth step change in turbulent boundary layers. J. Fluid Mech. 872, 472491.10.1017/jfm.2019.351CrossRefGoogle Scholar
Lund, T.S., Wu, X. & Squires, K.D. 1998 Generation of inflow data for spatially-developing boundary layer simulations. J. Comput. Phys. 140, 233258.10.1006/jcph.1998.5882CrossRefGoogle Scholar
Marusic, I., McKeon, B.J., Monkewitz, P.A., Nagib, H.M., Smits, A.J. & Sreenivasan, K.R. 2010 Wall-bounded turbulent flows at high Reynolds numbers: recent advances and key issues. Phys. Fluids 22, 065103.10.1063/1.3453711CrossRefGoogle Scholar
Meneveau, C. 2020 A note on fitting a generalised Moody diagram for wall modelled large-eddy simulations. J. Turbul. 21 (11), 650673.10.1080/14685248.2020.1840573CrossRefGoogle Scholar
Meneveau, C., Lund, T.S. & Cabot, W.H. 1996 A Lagrangian dynamic subgrid-scale model of turbulence. J. Fluid Mech. 319, 353385.10.1017/S0022112096007379CrossRefGoogle Scholar
Nikuradse, J. 1933 Laws of flow in rough pipes. NACA TM 1292. National Advisory Committee for Aeronautics.Google Scholar
Orlanski, I. 1976 A simple boundary condition for unbounded hyperbolic flows. J. Comput. Phys. 21, 251269.10.1016/0021-9991(76)90023-1CrossRefGoogle Scholar
Parthasarathy, A. & Saxton-Fox, T. 2023 A family of adverse pressure gradient turbulent boundary layers with upstream favourable pressure gradients. J. Fluid Mech. 966, A11.10.1017/jfm.2023.429CrossRefGoogle Scholar
Piomelli, U. 2008 Wall-layer models for large-eddy simulations. Prog. Aerosp. Sci. 44 (6), 437446.10.1016/j.paerosci.2008.06.001CrossRefGoogle Scholar
Piomelli, U. & Balaras, E. 2002 Wall-layer models for large-eddy simulations. Annu. Rev. Fluid Mech. 34 (1), 349374.10.1146/annurev.fluid.34.082901.144919CrossRefGoogle Scholar
Piomelli, U., Balaras, E. & Pascarelli, A. 2000 Turbulent structures in accelerating boundary layers. J. Turbul. 1 (1), 001.10.1088/1468-5248/1/1/001CrossRefGoogle Scholar
Piomelli, U., Ferziger, J.H., Moin, P. & Kim, J. 1989 New approximate boundary conditions for large eddy simulations of wall-bounded flows. Phys. Fluids A 1, 10611068.10.1063/1.857397CrossRefGoogle Scholar
Piomelli, U. & Yuan, J. 2013 Numerical simulations of spatially developing, accelerating boundary layers. Phys. Fluids 25, 101304121.10.1063/1.4825033CrossRefGoogle Scholar
Rao, K.S., Wyngaard, J.C. & Coté, O.R. 1974 The structure of the two-dimensional internal boundary layer over a sudden change of surface roughness. J. Atmos. Sci. 31 (3), 738746.10.1175/1520-0469(1974)031<0738:TSOTTD>2.0.CO;22.0.CO;2>CrossRefGoogle Scholar
Raupach, M.R., Antonia, R.A. & Rajagopalan, S. 1991 Rough-wall boundary layers. Appl. Mech. Rev. 44 (1), 125.10.1115/1.3119492CrossRefGoogle Scholar
Rotta, J.C. 1962 Turbulent boundary layers in incompressible flows. Prog. Aerosp. Sci. 2 (1), 195.10.1016/0376-0421(62)90014-3CrossRefGoogle Scholar
Salomone, T., Piomelli, U. & De Stefano, G. 2023 Wall-modeled and hybrid large-eddy simulations of the flow over roughness strips. Fluids 8, 10.10.3390/fluids8010010CrossRefGoogle Scholar
Schultz, M.P. & Flack, K.A. 2007 The rough-wall turbulent boundary layer from the hydraulically smooth to the fully rough regime. J. Fluid Mech. 580, 381405.10.1017/S0022112007005502CrossRefGoogle Scholar
Schumann, U. 1975 Subgrid-scale model for finite difference simulation of turbulent flows in plane channels and annuli. J. Comput. Phys. 18, 376404.10.1016/0021-9991(75)90093-5CrossRefGoogle Scholar
Shafi, H.S. & Antonia, R.A. 1995 Anisotropy of the Reynolds stresses in a turbulent boundary layer on a rough wall. Exp. Fluids 18, 213215.10.1007/BF00230269CrossRefGoogle Scholar
Shin, J.H. & Jin Song, S. 2014 Pressure gradient effects on smooth and rough surface turbulent boundary layers—Part I: favorable pressure gradient. J. Fluids Engng 137 (1), 011203.10.1115/1.4027474CrossRefGoogle Scholar
Shir, C.C. 1972 A numerical computation of air flow over a sudden change of surface roughness. J. Atmos. Sci. 29 (2), 304310.10.1175/1520-0469(1972)029<0304:ANCOAF>2.0.CO;22.0.CO;2>CrossRefGoogle Scholar
Song, S. & Eaton, J.K. 2002 The effects of wall roughness on the separated flow over a smoothly contoured ramp. Exp. Fluids 33, 3846.10.1007/s00348-002-0411-1CrossRefGoogle Scholar
Spalart, P.R. 1986 Numerical study of sink-flow boundary layers. J. Fluid Mech. 172, 307328.10.1017/S0022112086001751CrossRefGoogle Scholar
Spalart, P.R. & Watmuff, J.H. 1993 Experimental and numerical study of a turbulent boundary layer with pressure gradients. J. Fluid Mech. 249, 337371.10.1017/S002211209300120XCrossRefGoogle Scholar
Tachie, M.F. & Shah, M.K. 2008 Favorable pressure gradient turbulent flow over straight and inclined ribs on both channel walls. Phys. Fluids 20 (9), 095103.10.1063/1.2973205CrossRefGoogle Scholar
Tay, G.F.K., Kuhn, D.C.S. & Tachie, M.F. 2009 Influence of adverse pressure gradient on rough-wall turbulent flows. Intl J. Heat Fluid Flow 30, 249265.10.1016/j.ijheatfluidflow.2009.01.001CrossRefGoogle Scholar
Townsend, A.A. 1976 The Structure of Turbulent Shear Flow. Cambridge University Press.Google Scholar
Tsikata, J.M. & Tachie, M.F. 2013 Adverse pressure gradient turbulent flows over rough walls. Intl J. Heat Fluid Flow 39, 127145.10.1016/j.ijheatfluidflow.2012.11.001CrossRefGoogle Scholar
Volino, R.J. & Schultz, M.P. 2023 Comparison of smooth- and rough-wall non-equilibrium boundary layers with favourable and adverse pressure gradients. J. Fluid Mech. 959, A35.10.1017/jfm.2023.160CrossRefGoogle Scholar
Vreman, A. W. 2004 An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications. Phys. Fluids 16 (10), 36703681.10.1063/1.1785131CrossRefGoogle Scholar
Wang, L., Hu, R. & Zheng, X. 2020 A comparative study on the large-scale-resolving capability of wall-modeled large-eddy simulation. Phys. Fluids 32 (3), 035102.10.1063/1.5144669CrossRefGoogle Scholar
Wu, W. & Piomelli, U. 2018 Effects of surface roughness on a separating turbulent boundary layer. J. Fluid Mech. 841, 552580.10.1017/jfm.2018.101CrossRefGoogle Scholar
Yang, X.I.A., Sadique, J., Mittal, R. & Meneveau, C. 2015 Integral wall model for large eddy simulations of wall-bounded turbulent flows. Phys. Fluids 27, 025112132.10.1063/1.4908072CrossRefGoogle Scholar
Yuan, J. & Piomelli, U. 2014 Numerical simulations of sink-flow boundary layers over rough surfaces. Phys. Fluids 26. 015113.10.1063/1.4862672CrossRefGoogle Scholar
Yuan, J. & Piomelli, U. 2015 Numerical simulation of a spatially developing accelerating boundary layer over roughness. J. Fluid Mech. 780, 192214.10.1017/jfm.2015.437CrossRefGoogle Scholar
Figure 0

Table 1. Definitions of the momentum and displacement thicknesses required by the new LaRTE-RW model.

Figure 1

Figure 1. Inner-scaled mean velocity profiles for smooth and rough homogeneous surfaces, with the inset showing the corresponding roughness function. The DNS data from Lee & Moser (2015) are used for the smooth wall at ${\textit{Re}}_{\tau } = 5200$, while rough wall simulations are conducted at a roughness Reynolds number of $k_s^+ = 130$.

Figure 2

Figure 2. Sketch of the computational domain for heterogeneous-wall simulations (not to scale).

Figure 3

Figure 3. Normalized skin-friction coefficient, $C_{\!f}/C_{fr}$, where $C_{fr}$ is the skin-friction coefficient measured at $95\,\%$ of the rough strip’s length.

Figure 4

Figure 4. Normalized skin-friction coefficient on the rough strip (a) and smooth strip (b). $C_{fr}$ and $C_{fs}$ are the skin-friction coefficient measured at $95\,\%$ of the rough- and smooth-strip lengths, respectively. Experimental uncertainty from Antonia & Luxton (1971) (a) and Li et al. (2019) (b) is depicted with grey shading.

Figure 5

Figure 5. Instantaneous contours of the wall-stress (a–c) and the stress at the wall-model interface (d) with LaRTE-RW model for the streamwise component.

Figure 6

Figure 6. Computational domain for the varying pressure gradients simulations. The green area is a visual representation of the applied FPG–APG condition.

Figure 7

Figure 7. Prescribed free-stream velocity profiles $ U_{\infty } / U_{\textit{ref}}$, fitted to match experimental data at seven streamwise locations for each angle of attack. The corresponding acceleration parameter is shown in the insets. Light blue highlights FPG regions, light yellow highlights APG regions. Lines represent simulations; black circles denote experimental data from the Virginia Tech Stability Wind Tunnel (VPI), as reported by Fritsch et al. (2022a).

Figure 8

Figure 8. Streamwise skin-friction coefficient development for different wall models at each angle of attack, compared with experimental data from the VPI reported by Fritsch et al. (2022a). The experimental uncertainty is depicted with grey shading.

Figure 9

Figure 9. Outer-units velocity profiles for each angle of attack. From bottom to top, the seven profiles correspond to experimental locations at $x / L_{\textit{ref}}=0, 0.11, 0.22, 0.34, 0.45, 0.56$ and $0.67$. Each velocity profile, except the first one, is shifted upward by 0.5 units for clarity. Black circles represent experimental data from the VPI, as reported by Fritsch et al. (2022a). The dashed line represents the location of the wall-model interface.

Figure 10

Figure 10. Streamwise development of the skin-friction coefficient $ C_{\!f}$ for two representative angles of attack: (a) FPG/APG case $ \alpha = -10^\circ$, and (b) FPG/APG/FPG case $ \alpha = 12^\circ$. Results are shown with the turbNEQ contribution included (dot–dashed line) and excluded (dotted line). Black circles represent experimental data from the VPI, as reported by Fritsch et al. (2022a). Experimental uncertainty is depicted with grey shading.

Figure 11

Figure 11. (a) Streamwise development of the skin-friction coefficient $C_{\!f}$, with anisotropic grid refinement. (b) Comparison of Vreman SFS and LDSM for log-law EQWM, and LaRTE-RW at angle of attack $\alpha = 12^\circ$. Black circles represent experimental data from the VPI, as reported by Fritsch et al. (2022a). Experimental uncertainty is depicted with grey shading.

Figure 12

Figure 12. Inner-similarity function (a) and near-wall function (b) in the LaRTE-RW model.