2.1. Equations of motion

We begin by introducing the equations of motion for the quasi-linear approximation of Hwang & Eckhardt (Reference Hwang and Eckhardt2020) in fully-developed incompressible turbulent channel flow, bounded by two plates that are infinitely long and wide. Here, the streamwise, wall-normal and spanwise coordinates are denoted as $x$ , $y$ and $z$ , respectively, and the corresponding velocity vector is represented as $\boldsymbol{u} = (u, v, w)$ . The channel half-height is set as $h$ , and the two plates are located at $y = 0, 2h$ , respectively. To implement the quasi-linear approximation, the velocity field is decomposed into a time-averaged mean $\boldsymbol{U} = (U(y), 0, 0)$ , and fluctuations about this mean $\boldsymbol{u}' = (u', v', w')$ : i.e. the Reynolds decomposition. This yields the following coupled equations for the mean and fluctuations:

(2.1a) \begin{align}&\qquad\qquad\qquad \nu \frac {{\rm d}U}{{\rm d}y}-\overline {u'v'} = \frac {\tau _w}{\rho }\left(1-\frac {y}{h}\right), \end{align}

(2.1b) \begin{align}& \frac {\partial \boldsymbol{u}'}{\partial t}+ (\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} )\boldsymbol{u}' + (\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{U} =-\frac {1}{\rho }\,\boldsymbol{\nabla} {p'}+\nu\, {\nabla} ^2 {\boldsymbol{u}'}+ \boldsymbol{N}, \end{align}

where

(2.1c) \begin{equation} \boldsymbol{N} = -\boldsymbol{\nabla} \boldsymbol{\cdot} \left(\boldsymbol{u}'\boldsymbol{u}'-\overline {\boldsymbol{u}'\boldsymbol{u}'}\right). \end{equation}

Here, $t$ denotes time and $\overline {({\cdot} )}$ represents the average in time and two homogeneous directions, with $\tau _w$ the time-averaged wall shear stress, $p'$ the fluctuating pressure, $\rho$ the fluid density, and $\nu$ the kinematic viscosity. In the present study, the effects of roughness on near-wall turbulence are modelled with the TRM boundary conditions (e.g. Lācis et al. Reference Lācis, Sudhakar, Pasche and Bagheri2020). The boundary conditions for the turbulent mean velocity are written as

(2.2) \begin{equation} {U} =l_x\left .\frac {\partial {U}}{\partial y}\right \vert _{y=0,2h}, \end{equation}

and the fluctuating velocities at the two boundaries $y=0,2h$ are defined as

(2.3a) \begin{align}&\qquad\qquad\,\, {u'} =l_x\left .\frac {\partial {u'}}{\partial y}\right \vert _{y=0,2h}, \end{align}

(2.3b) \begin{align}&\qquad\qquad {w'} =l_z\left .\frac {\partial {w'}}{\partial y}\right \vert _{y=0,2h}, \end{align}

(2.3c) \begin{align}& {v'} =-m_x\left .\frac {\partial {u'}}{\partial x}\right \vert _{y=0,2h}-m_z \left .\frac {\partial {w'}}{\partial z}\right \vert _{y=0,2h}, \end{align}

where the coefficients $l_x$ , $l_z$ , $m_x$ and $m_z$ are slip and transpiration lengths, respectively. The TRM boundary conditions not only incorporate the Navier-slip boundary conditions for the streamwise and spanwise velocities, but also provide a transpiration condition for the wall-normal velocity that accounts for the wall-normal momentum transport associated with the surface texture through the transpiration lengths. It has been shown that these boundary conditions model well for the rough surfaces in the transitionally rough regime, which do not significantly modify the near-wall turbulence. From this perspective, it is ideal to test the modelling capacity of the quasi-linear approximations developed for smooth walls (Hwang & Eckhardt Reference Hwang and Eckhardt2020; Holford et al. Reference Holford, Lee and Hwang2024a ).

The mean equation (2.1a ) includes the Reynolds shear stress term that feeds back from the fluctuating velocity field. The equation for the fluctuating velocity (2.1b ) is of the form linearised about the mean velocity if the nonlinear term $\boldsymbol{N}$ is dropped. Following Hwang & Eckhardt (Reference Hwang and Eckhardt2020) and Holford et al. (Reference Holford, Lee and Hwang2024a ), the nonlinear term $\boldsymbol{N}$ is modelled as

(2.4a) \begin{equation} \boldsymbol{N}_{\nu _t,f} = \boldsymbol{\nabla} \boldsymbol{\cdot} \big(\nu _t \big(\boldsymbol{\nabla} \boldsymbol{u}'+ \boldsymbol{\nabla} {\boldsymbol{u}'}^{\rm T}\big)\big)+\boldsymbol{f}', \end{equation}

where the eddy viscosity $\nu _t$ is adopted from the empirical expression of Cess (Reference Cess1958),

(2.4b) \begin{equation} \nu _t(\eta ) = \frac {\nu }{2} \left\{1+\frac {\kappa ^2\,Re_{\tau }^2}{9}\left(1-\eta ^2\right)^2\left(1+2\eta ^2\right)^2(1-\exp[(|\eta |-1)\,Re_{\tau }/A])^2\right\}-\frac {\nu }{2}, \end{equation}

with $\eta = (y-h)/h$ , and $\boldsymbol{f}'=(f'_u,f'_v,f'_w)$ is a stochastic forcing term whose amplitude and colour are to be determined. Here, the eddy viscosity parameters are set to $\kappa = 0.426$ and $A = 25.4$ , with which the best least squares fitting of the mean velocity profile of (2.1a ) is obtained with that of DNS at $Re_{\tau }\approx 2000$ (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006). The inclusion of the eddy viscosity term in (2.4a ) has been understood to significantly simplify the modelling process, and consolidates the key physical elements, as recently analysed in detail by Holford et al. (Reference Holford, Lee and Hwang2024b ) through a comparison of spectral energy budgets between DNS and quasi-linear approximations. Furthermore, the eddy viscosity term has been consistently shown to enhance the capacity of (2.1b ) for the predictions of turbulent statistics and energy spectra, as well as for the performance of state estimations (Hwang Reference Hwang2016; Illingworth, Monty & Marusic Reference Illingworth, Monty and Marusic2018; Madhusudanan, Illingworth & Marusic Reference Madhusudanan, Illingworth and Marusic2019; Towne, Lozano-Durán & Yang Reference Towne, Lozano-Durán and Yang2020; Holford et al. Reference Holford, Lee and Hwang2023).

As discussed earlier, the formulation of the quasi-linear approximations in this study requires a known mean velocity. Indeed, like RANS, (2.1) has a closure problem due to unknown $\boldsymbol{N}$ – if $\boldsymbol{N}$ is known, then (2.1) is solvable to obtain the mean and fluctuating velocities. Conversely, if part of the solution, say the mean velocity, is assumed to be known, then it is possible to construct $\boldsymbol{N}$ inversely, although $\boldsymbol{N}$ is not necessarily unique (Skouloudis & Hwang Reference Skouloudis and Hwang2022; see also § 2.2). This is the central theoretical element of the present quasi-linear approximations. In the present study, the construction of $\boldsymbol{N}$ is performed using $\boldsymbol{N}_{\nu _t,f}$ in (2.4), where the stochastic forcing term is determined with the known mean velocity. The mean velocity here is obtained by solving (2.1a ), with the well-known eddy-viscosity-based closure given by $-\overline {u'v'}=\nu _t\, {{\rm d}U}/{{\rm d}y}$ with $\nu _t$ in (2.4b ).

2.2. Forcing models

To apply the framework of the stochastic linear dynamical system for (2.1b ) (Farrell & Ioannou Reference Farrell and Ioannou1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; Hwang & Cossu Reference Hwang and Cossu2010), we initially consider that the forcing is temporally white and spatially and componentwisely uncorrelated. Given the homogeneous nature in the wall-parallel directions, we perform the Fourier transform along the $x$ and $z$ directions,

(2.5) \begin{equation} \,\hat {\!\boldsymbol{f}}'(t,y;k_x,k_z)=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty } \boldsymbol{f}'(t,x,y,z)\,{\rm e}^{-{\rm i}(k_xx+k_zz)}\,{\rm d}x\,{\rm d}z, \end{equation}

where $\hat {(\cdot)}$ denotes the Fourier transformed state, and $(k_x, k_z)$ is the considered wavenumber pair with the streamwise and spanwise wavelengths denoted as $\lambda _x = 2\pi /k_x$ and $\lambda _z = 2\pi /k_z$ , respectively. The same definitions of the Fourier transform are used for the other flow variables. Then the spectral covariance matrix of the stochastic forcing is set to be given by

(2.6) \begin{equation} \begin{aligned} &E\left[\,\hat {\!\boldsymbol{f}}'(t,y;k_x,k_z)\,\hat {\!\boldsymbol{f}}'^H(t',y';k_x,k_z)\right]\\&= \begin{bmatrix} W_{u}(k_x,k_z) &0 &0\\[4pt] 0& W_{v}(k_x,k_z) &0\\[4pt] 0 &0 & W_{w}(k_x,k_z) \end{bmatrix} \delta (y-y')\,\delta (t-t'), \end{aligned} \end{equation}

where $(\cdot)^H$ denotes the complex conjugate transpose, $E[\cdot]$ is the expectation operator on different stochastic realisations, $\delta (\cdot)$ is the Dirac delta function, and $W_{l}$ with $l=\{u,v,w\}$ are component weights as functions of $k_x$ and $k_z$ , through which the amplitude of the stochastic forcing is determined in the streamwise and spanwise directions.

The following spectral (spatial and componentwise) covariance matrix of velocity fluctuations can be computed by solving (2.1b ) in the Fourier space, if the stochastic forcing is specified:

(2.7) \begin{equation} \phi _{\boldsymbol{\it uu}}\left(y,y';k_x,k_z\right)=E\left[\hat {\boldsymbol{u}}'(t,y;k_x,k_z)\,\hat {\boldsymbol{u}}'^H(t,y';k_x,k_z)\right]. \end{equation}

Given the existence of a linear relationship between the velocity and forcing spectral covariance matrices, $\phi _{\boldsymbol{\it uu}}(y,y'; k_x,k_z)$ is written as

(2.8) \begin{equation} \phi _{\boldsymbol{\it uu}}(y,y';k_x,k_z)=\sum _{{l}=u,v,w} W_{{l}}(k_x,k_z)\, \phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z), \end{equation}

where $\phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z)$ is the velocity spectral covariance matrix corresponding to each component of the forcing with the unit amplitude (Holford et al. Reference Holford, Lee and Hwang2024a ). In particular, using the linear nature of (2.4a ), $\phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z)$ in (2.8) is obtained from the linear response to the white-in-time forcing with the unit amplitude (i.e. setting $W_{{l}}(k_x,k_z)=1$ for each $l\ (=u,v,w)$ ). For instance, $\phi _{\boldsymbol{\it uu},u}(y,y';k_x,k_z)$ in (2.8) is the velocity spectral covariance matrix computed with the forcing

(2.9) \begin{equation} \begin{aligned} E[&\,\hat {\!\boldsymbol{f}}'(t,y;k_x,k_z)\,\hat {\!\boldsymbol{f}}'^H(t',y';k_x,k_z)]= \begin{bmatrix} \delta (y-y') &0&0\\ 0& 0 &0\\ 0 &0 & 0 \end{bmatrix} \delta (t-t'). \end{aligned} \end{equation}

The stochastic forcing in (2.6) serves as a starting point for the quasi-linear approximation approach here. With the pre-computed $\phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z)$ , the remaining problem becomes the determination of $W_l$ in such a way that the resulting Reynolds shear stress corresponds to that in (2.1) using the prescribed mean velocity (for the details, see Holford et al. Reference Holford, Lee and Hwang2024a ). However, a forcing that is white in time and decorrelated in the wall-normal direction for each $k_x$ and $k_z$ is physically unrealistic. Indeed, this kind of forcing covariance has previously been shown to generate a non-physical large energy response in the region close to the channel centre (Hwang & Eckhardt Reference Hwang and Eckhardt2020; Holford et al. Reference Holford, Lee and Hwang2024a ). It has been proposed that a simple ad hoc solution to this is considering some of the leading proper orthogonal decomposition (POD) modes from the response to the stochastic forcing, which can capture the primary energy-containing motions at integral length scales – the reader may refer to Hwang & Eckhardt (Reference Hwang and Eckhardt2020) and Holford et al. (Reference Holford, Lee and Hwang2024b ), where the origin of the non-physical stochastic response, the effect of the number of the POD modes, and the related energetics are fully discussed for the quasi-linear approximations considered here. Employing this approach, the velocity spectral covariance matrix is further approximated by considering some of the leading POD modes from $\phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z)$ in (2.8):

(2.10) \begin{equation} \phi _{\boldsymbol{\it uu},l}^{\textit{POD}}(y,y';k_x,k_z) \approx \sum _{i=1}^{N_{\textit{POD}}} \sigma _i\, \hat {\boldsymbol{u}}_{i,POD}(y;k_x,k_z)\,\hat {\boldsymbol{u}}_{i,POD}^H(y';k_x,k_z), \end{equation}

where $N_{\textit{POD}}$ is the number of the considered leading POD modes, and $\sigma _i$ and $\hat {\boldsymbol{u}}_{i,POD}$ are the eigenvalues and eigenfunctions of the velocity spectral covariance matrix obtained from the white-in-time and spatially decorrelated forcing $\phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z)$ in (2.8). Based on previous findings (Hwang & Eckhardt Reference Hwang and Eckhardt2020; Holford et al. Reference Holford, Lee and Hwang2024a ), $N_{\textit{POD}} = 2$ is considered in this study.

Following Holford et al. (Reference Holford, Lee and Hwang2024a ), the forcing weight $W_{l}(k_x,k_z)$ for each component is further prescribed in the form

(2.11) \begin{equation} W_{l}(k_x,k_z)=W_{l,k_x}(k_x/k_z)\,W_{k_z}(k_z). \end{equation}

Here, the self-similar streamwise weight $W_{l,k_x}(k_x/k_z)$ is used to account for the self-similar nature of the energy-containing part in the velocity and forcing spectra (i.e. the attached eddy hypothesis; for further details, see Holford et al. Reference Holford, Lee and Hwang2023), and $W_{k_z}(k_z)$ is used to determine the forcing amplitude for each spanwise scale that varies from the near-wall ( $\lambda _z^+ \approx O(100)$ ) to the outer one ( $\lambda _z/h \approx O(1)$ ) (Hwang & Eckhardt Reference Hwang and Eckhardt2020).

We note that the streamwise weight $W_{l,k_x}(k_x/k_z)$ is essentially a modelling choice for the quasi-linear approximations in the present study. This is because, for any arbitrary choice of $W_{l,k_x}(k_x/k_z)$ , there may well exist $W_{k_z}(k_z)$ that provides the Reynolds shear stress identical to the one given by the prescribed mean velocity profile (for a further discussion, see Skouloudis & Hwang Reference Skouloudis and Hwang2022). In the present study, therefore, three models for $W_{l,k_x}(k_x/k_z)$ are considered.

1. (i) The MQLA model considers the form of the self-similar streamwise weight as

   (2.12) \begin{equation} W_{l,k_x}(k_x/k_z)=\tfrac{1}{3}\delta (k_x/k_z) \end{equation}
   for all $l$ (Hwang & Eckhardt Reference Hwang and Eckhardt2020). This is the simplest form of the quasi-linear approximations, since it completely ignores the statistical structure of the velocity spectra in the streamwise direction. For the same reason, it is computationally the cheapest option, since the covariance matrix $\phi _{\boldsymbol{\it uu},l}$ needs to be pre-computed only for $k_x=0$ .

2. (ii) The DQLA model considers the streamwise weight $W_{l,k_x}(k_x/k_z)$ obtained by assimilating DNS data in the logarithmic region at $Re_\tau =5200$ (Holford et al. Reference Holford, Lee and Hwang2024a ). In the present study, $W_{l,k_x}(k_x/k_z)$ obtained for $k_z h=126$ in Holford et al. (Reference Holford, Lee and Hwang2024a ) is used (see their figure 1). This is a spanwise wavenumber associated with the logarithmic region (i.e. $100/Re_{\tau }\lesssim \lambda _z/h \lesssim 1.5$ ). We note that in Holford et al. (Reference Holford, Lee and Hwang2024a ), $k_z h=14,30,50,76,126$ were tested for $Re_\tau =5200$ . The chosen $k_z h=126$ corresponds to $\lambda _z^+=259$ at $Re_\tau =5200$ and is most relevant to the region near the wall among all the weights computed. By construction, DQLA offers the full 2-D spectra as a result of the non-zero streamwise weight. However, this improvement comes at a largely increased computational cost, as it requires pre-computations of the covariance matrix $\phi _{\boldsymbol{\it uu},l}$ in the full range of $k_x \ne 0$ for all velocity components – indeed, the computational cost of DQLA is typically a few hundred times greater than that of MQLA (see table 1).

3. (iii) The minimal two-mode quasi-linear approximation (M2QLA) model in the present study considers the following simple form of the streamwise weight $W_{l,k_x}(k_x/k_z)$ :

   (2.13) \begin{equation} W_{l,k_x}(k_x/k_z)=\tfrac{1}{3}\left [\delta (k_x/k_z)+ \delta (k_x/k_z \pm 0.5)\right ] \end{equation}
   for all $l$ . By construction, this model resembles MQLA, but includes a minimal variation in the streamwise direction. The first term on the right-hand side of (2.13) is designed to model long streaky structures in the near-wall, logarithmic and outer regions ( $\lambda _x\gt 10 \lambda _z$ ), while the second term is for the related streamwise wavy motions involving quasi-streamwise vortical strictures ( $\lambda _x \approx 2 \lambda _z$ ) (for a detailed discussion, see Hwang Reference Hwang2015). The computational cost of M2QLA is only a couple of times greater than that of MQLA, as it only requires pre-computations of $\phi _{\boldsymbol{\it uu},l}$ for $k_x=0$ and $k_x/k_z=\pm 0.5$ .

Table 1.

Numerical and optimisation parameters used in the present study at $Re_{\tau }=180$ : $N_y$ , the number of wall-normal collocation points; $N_{k_x}$ , the number of streamwise wavenumbers; $N_{k_z}$ , the number of spanwise wavenumbers. Here, $\|\cdot\|_Q^2 \equiv ({u_\tau ^2/h})\int _0^{2h}(\cdot)^2\, Q(y)\, {\rm d}y$ and $\|\cdot\|_{L_2}^2\equiv ({u_\tau ^2/h})\int _0^{2h}(\cdot)^2\, {\rm d}y$ .

2.3. Quasi-linear approximations

Thus far, we have discussed how we construct the forcing for each $k_z$ , i.e. the self-similarity-based decomposition of $W_l(k_x,k_z)=W_{l,k_x}(k_x/k_z)\,W_{k_z}(k_z)$ , the form of $W_{l,k_x}(k_x/k_z)$ , and the use of POD modes. Now, we discuss how the forcing structure for each $k_z$ is used to determine $W_{k_z}(k_z)$ , which will make the fluctuations from (2.1b ) consistent with the mean equation (2.1a ).

Given the discussion in § 2.2, the final form of the velocity covariance matrix in the Fourier space is written as

(2.14) \begin{equation} \phi _{\boldsymbol{\it uu}}(y,y';k_x,k_z)=W_{k_z}(k_z) \sum _{l=u,v,w} W_{{l},k_x}(k_x/k_z)\, \phi _{\boldsymbol{\it uu},l}^{\textit{POD}}(y,y';k_x,k_z), \end{equation}

with the dimensions $3N_{y}\times 3N_{y}$ , for each $k_x$ and $k_z$ . Before applying the quasi-linear approximation for the TRM boundary condition in (2.2) and (2.3), we establish the solution procedure of the quasi-linear approximations for the no-slip boundary condition. Following the procedure documented in Willis, Hwang & Cossu (Reference Willis, Hwang and Cossu2010), the spectral covariance operator $\phi _{\boldsymbol{\it uu},l}(y,y';k_x,k_z)$ in (2.8) is first obtained from the response of the linearised Navier–Stokes equations to white-in-time noise provided componentwise with the unit amplitude: for example, $\phi _{\boldsymbol{\it uu},u}^{\textit{POD}}(y,y';k_x,k_z)$ is obtained by setting $W_u=1$ and $W_v=W_w=0$ . For this purpose, the eigenvalue problem from (2.1b ) with the nonlinear term model $\boldsymbol{N}_{\nu _t,f}$ in (2.4a ) is first solved in terms of the primitive variables $(u',v',w',p')$ using a Chebyshev collocation method (Weideman & Reddy Reference Weideman and Reddy2000). The Lyapunov equation for the velocity spectral covariance operator $\phi _{\boldsymbol{\it uu},l}$ is subsequently formulated in the resulting eigenfunction space and is solved using the lyap function in MATLAB (see Appendix A for further details). Finally, the velocity spectral covariance operator $\phi _{\boldsymbol{\it uu},l}$ is further approximated using the two leading POD modes (i.e. $N_{\textit{POD}}=2$ ), as discussed with (2.10).

Now, $\phi _{\boldsymbol{\it uu},l}^{\textit{POD}}(y,y';k_x,k_z)$ was determined as detailed above, and $W_{l,k_x}(k_x/k_z)$ was prescribed through § 2.2. Therefore, only the determination of $W_{k_z}(k_z)$ remains. For the determination of the spanwise weight $W_{k_z}(k_z)$ , an optimisation problem is formulated to minimise the difference between the Reynolds shear stress required for the mean profile $\overline {u'v'}(y)$ and that generated by the fluctuation equation $E[{u}'{v}'](y)$ . In particular, the following optimisation problem for $W_{k_z}(k_z)$ is considered following Holford et al. (Reference Holford, Lee and Hwang2024a ):

(2.15a) \begin{align} \min _{W_{k_z}} \left [\!\frac {\int _0^{2h}(\overline {u'v'}(y)-E[u'v'](y))^2\,Q(y)\,{\rm d}y}{\int _0^{2h}(\overline {u'v'}(y))^2\,Q(y)\,{\rm d}y} \!\right ]^{0.5}\!+\gamma \left [\! \int _0^{\infty }\!\left (\frac {{\rm d}^2W_{k_z}(k_z)}{{\rm d}(\ln k_z)^2}\right )^2\! R_{uv}(k_z)\,{\rm d}k_z\!\right ]^{0.5}\nonumber\\[4pt] \end{align}

subject to

(2.15b) \begin{equation} W_{k_z}(k_z)\geqslant 0, \end{equation}

where $Q(y)=(1-|\eta |)^{-1}$ to place equal emphasis on wall-normal grid points via a logarithmic scaling with distance from the wall. The first term of (2.15a ) is a relative error in the Reynolds shear stress with the term introduced in its denominator, and the second term is a regularisation term for the generation of a smooth $W_{k_z}(k_z)$ in the logarithmic coordinate of $k_z$ (see below for a further explanation). The combination of the two square-rooted quadratic functions of $W_{k_z}(k_z)$ enables us to use the second-order cone programming (e.g. Holford et al. Reference Holford, Lee and Hwang2023, Reference Holford, Lee and Hwang2024a ). From this point of view, different forms of the optimisation problem can be formulated to achieve the same goal, as in Hwang & Eckhardt (Reference Hwang and Eckhardt2020) and Skouloudis & Hwang (Reference Skouloudis and Hwang2022).

In (2.15), the weight function $W_{k_z}(k_z)$ is constrained to be in the range between $ \lambda _z^+=10$ and $ \lambda _z=10h$ , with the values of $W_{k_z}(k_z)$ at the smallest and largest spanwise wavelengths set to zero, i.e. $W_{k_z}(\lambda ^+_z= 10) =W_{k_z} (\lambda _z = 10h) = 0$ . Given the spanwise scale typically varying from $\lambda _z^+ \approx 100$ in the near-wall region to $\lambda _z \approx 1.5h$ in the outer region (Hwang Reference Hwang2015), such a constraint on $W_{k_z}(k_z)$ is able to cover a full range of spanwise length scales, from near-wall motions to very large-scale motions in the outer layer. The Reynolds shear stress, denoted as $\overline {u'v'}(y)$ , is solved by considering the classical mixing length model with the eddy viscosity in (2.4b ) for the mean momentum equation in (2.1a ): $-\overline {u'v'}(y)=\nu _t {\rm d}U/{\rm d}y$ . The Reynolds shear stress $E[u'v'](y)$ , generated from the fluctuation equation (2.1b ) with the incorporation of the nonlinear term model (2.4a ), is formulated using

(2.16) \begin{equation} \begin{aligned} E[u'v'](y)=&\frac {1}{4\pi ^2}\int _0^{2h}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty } \phi _{uv}(y,y';k_x,k_z)\,\delta (y'-y)\,{\rm d}k_x\,{\rm d}k_z\,{\rm d}y', \end{aligned} \end{equation}

where $\phi _{uv}$ is from (2.1). The constraint $W_{k_z}(k_z)\geqslant 0$ is given to ensure that the forcing covariance remains positive definite. Additionally, $W_{k_z}$ must exhibit sufficient smoothness to prevent the existence of any non-physical flow features. To ensure this smoothness in the logarithmic $k_z$ coordinate, a global regularisation term with the related second-order derivative of $W_{k_z}$ is included in (2.15a ), with $\gamma$ controlling its relative importance, such that a set of physically plausible smooth velocity spectra can be achieved. The smoothness regularisation is further weighted by $R_{uv}(k_z)$ , which is defined as

(2.17) \begin{equation} R_{uv}(k_z)=\frac {1}{2\pi }\int _h^{2h}\int _{-\infty }^{\infty } \sum _{{l}=u,v,w} W_{{l},k_x}(k_x/k_z)\, \phi _{uv,{l}}^{\textit{POD}} (y,y;k_x,k_z)\,{\rm d}k_x\,{\rm d}y, \end{equation}

and it tends to have large values at large spanwise wavelengths. Since spectra of the Reynolds shear stress contain large energy at large spanwise wavelengths, the regularisation term is more heavily weighted at such wavelengths. This weighting accounts for the rapid decay in Reynolds shear stress spectra observed in previous studies (Skouloudis & Hwang Reference Skouloudis and Hwang2022; Holford et al. Reference Holford, Lee and Hwang2024a ). It prevents erroneous behaviours in velocity spectra at larger scales, encouraging a smoothly attached compact support at the large spanwise length scales.

The optimisation problem in (2.15) is discretised along the streamwise and spanwise wavenumber axis using a uniform logarithmic spacing, with $\varDelta (\ln k_xh) \leqslant 0.15$ for $\lambda _x/h \in [10/Re_\tau ,100]$ (DQLA) and $\varDelta (\ln k_zh) \leqslant 0.15$ for $\lambda _z/h \in [10/Re_\tau ,10]$ ( $Re_\tau \equiv u_\tau h/\nu$ , where $u_\tau$ is the friction velocity). The integration with respect to the streamwise and spanwise wavenumbers is conducted using the trapezoidal method in the logarithmic coordinates. The optimisation problem in (2.15) is formulated in a standard form of a second-order cone program, such that it can be efficiently solved with the MOSEK solver in MATLAB; for further details, the reader may refer to Holford et al. (Reference Holford, Lee and Hwang2024a ). Table 1 shows the details on the number of streamwise and spanwise wavenumbers, collocation points in the wall-normal direction, the value of the penalty $\gamma$ in (2.15), and the associated errors in the $Q$ and $L_2$ norms of the optimisation problem at $Re_{\tau }=180$ . Figure 1(a) shows the weight $W_{k_z}(k_z)$ obtained by solving (2.15) for M2QLA. Using this weight, an almost exact agreement can be achieved between the Reynolds shear stress profiles from the mean equation (2.1a ) and the fluctuation equations (2.1b ), as shown in figure 1(b). The same level of agreement is also obtained for $W_{k_z}(k_z)$ for MQLA and DQLA (table 1).

Figure 1.

Outputs of M2QLA optimisation: (a) the spanwise weight of Fourier modes; (b) the Reynolds shear stress profiles from the mean equation (blue dashed line) and the fluctuating equations (red solid line). Here, $\gamma =0.001$ , $Re_{\tau }=180$ and $N_y=98$ .

Figure 2.

Outputs of optimisation: (a) $u^{\prime }_{rms}/ u_{\tau }$ , (b) $v^{\prime }_{rms}/u_{\tau }$ , (c) $w^{\prime }_{rms}/u_{\tau }$ , with MQLA (red), DQLA (blue) and M2QLA (green) at $Re_{\tau }=180$ . The optimisation parameters are listed in table 1. The DNS results from Lee & Moser (Reference Lee and Moser2015) are also plotted for comparison (black).

Figure 2 compares the velocity fluctuations from the present quasi-linear approximations with those of DNS (Lee & Moser Reference Lee and Moser2015) at $Re_{\tau }=180$ . The MQLA shows the most anisotropic turbulent velocity fluctuations: the streamwise fluctuations are much greater than those of DNS, while the cross-streamwise ones are much smaller. This is because streamwise-dependent Fourier modes are excluded, preventing any energy transfer from the streamwise velocity to the rest through the streamwise pressure strain Holford et al. (Reference Holford, Lee and Hwang2024a ,Reference Holford, Lee and Hwang b ). Indeed, inclusion of an additional streamwise-varying Fourier mode improves this undesirable anisotropy, as seen in the velocity fluctuations from M2QLA. Finally, the anisotropy of DQLA is comparable to that of DNS, as expected. Despite the encouraging improvement in DQLA, its turbulence statistics still show some non-negligible differences from those of DNS, especially in the wall-normal and spanwise velocity fluctuations. This is because the simple nonlinear term model in (2.4a ) and the inclusion of only two leading POD modes do not correctly model the nonlinear processes involved in the pressure-strain transport, especially by slow pressure. For a detailed discussion, the reader may refer to Holford et al. (Reference Holford, Lee and Hwang2024b ).

2.4. The TRM boundary conditions

Having established the quasi-linear approximations for no-slip boundary conditions, we now extend this framework for the case with TRM boundary conditions. The key task here is to develop a quasi-linear framework that predicts the change in the mean and fluctuation velocities by TRM boundary conditions subject to the equations of motion. We start by considering the steady mean momentum equation for a rough surface modelled by TRM boundary conditions:

(2.18) \begin{equation} -\frac {1}{\rho }\frac {{\rm d} P_0}{{\rm d} x} + \nu \frac {{\rm d}^2 U_r}{{\rm d} y^2}-\frac {{\rm d} \overline {u'v'}_r}{{\rm d} y}=0, \end{equation}

where the subscript $r$ denotes the flow quantities over the rough surface, and ${\rm d}P_0/{\rm d}x$ is the applied streamwise pressure gradient. The mean velocity and Reynolds shear stress are first decomposed as $U_r=U_s+\delta U$ and $\overline {u'v'}_{\!\!r}= \overline {u'v'}_{\!\!s} + \delta \overline {u'v'}$ , with the subscript $s$ representing the quantities on the smooth wall. Here, the applied streamwise pressure gradient is chosen to be unchanged in both rough- and smooth-wall cases. Substituting the decompositions into (2.18), and extracting the mean equation without roughness for the given streamwise pressure gradient, yields

(2.19a) \begin{equation} \nu \frac {{\rm d}^2 \delta U}{{\rm d} y^2} - \frac {{\rm d} \delta \overline {u'v'}}{{\rm d} y}=0. \end{equation}

Similarly, substituting the decomposition for mean velocity into (2.2) and extracting the boundary condition over the smooth wall lead to the boundary condition for $\delta U$ :

(2.19b) \begin{equation} \left [\delta U-{l_x}\frac {{\rm d} \delta U}{{\rm d} y}\right ]_{y=0,2h}={l_x}\left .\frac {{\rm d} U_s}{{\rm d} y}\right |_{y=0,2h}. \end{equation}

In the absence of a streamwise pressure gradient term, (2.19) can be solved if $\delta \overline {u'v'}$ is given. Within the present quasi-linear approximation framework, $\delta \overline {u'v'}$ will have to be obtained by solving the turbulent fluctuation (2.1b ) with $\boldsymbol{N}_{\nu _t,f}$ in (2.4a ) subject to the TRM boundary conditions (2.3). However, (2.1b ) is also coupled with (2.19a ) through the mean velocity, indicating that (2.19a ) must be solved simultaneously with (2.1b ).

The quasi-linear approximations in Hwang & Eckhardt (Reference Hwang and Eckhardt2020) and Holford et al. (Reference Holford, Lee and Hwang2024a ) were about the self-consistent determination of the forcing term in (2.4a ), with the assumption that mean velocity is known prior. In contrast, the problem here is to predict the changes in mean and fluctuating velocities due to the TRM boundary condition. In this case, the form of $\boldsymbol{N}_{\nu _t,f}$ should be known prior to the determination of the changes in mean and fluctuations – otherwise, the problem will have a closure issue. A simple closure for $\boldsymbol{N}_{\nu _t,f}$ that we have chosen in the present study is to fix $\nu _t$ in (2.4b ), and $W_l(k_x,k_z)$ in (2.6), for $\boldsymbol{N}_{\nu _t,f}$ in (2.4a ), and this is seen to perform well (see § 3) as long as the flow modification by the TRM boundary condition remains modest (i.e. the smooth-wall-like regime). This choice certainly deserves further discussion, which will be given at the end of this section.

Given the form of $\boldsymbol{N}_{\nu _t,f}$ , we now introduce the procedure to obtain the mean and fluctuating velocities. To solve (2.19a ), we introduce an artificial unsteady term, such that

(2.20) \begin{equation} \frac {\partial \delta U}{\partial \tau }=\nu \frac {\partial ^2 \delta U}{\partial y^2} - \frac {\partial \delta \overline {u'v'}}{\partial y} \end{equation}

with an artificial time $\tau$ , while obtaining $\delta \overline {u'v'}$ by solving the algebraic (or steady) Lyapunov equation related to (2.1b ) with the boundary conditions in (2.3). Equation (2.20) is subsequently discretised with a mixed first-order Euler method, where the viscous term is treated implicitly, while the Reynolds shear stress term is done explicitly – we note that given the purpose of obtaining the steady solutions, the accuracy of the time integration here is not very important. At each time step, the Reynolds shear stress is obtained by solving the Lyapunov equation for (2.1b ) with the mean velocity $U_r$ , updated with $\delta U$ in (2.20), and the boundary conditions (2.3). The initial condition for (2.20) is obtained by solving (2.19a ) with the mixing length model, where $-\overline {u'v'}=\nu _t\, {\rm d}U/{\rm d}y$ , with $\nu _t$ from (2.4b ). The time integration is performed with $\Delta \tau\, u_\tau /h =0.003$ until it provides a well-converged solution for (2.19a ) with the given quasi-linear approximation framework in § 2.3. Further details on the algorithm described here are summarised in Algorithm1.

Algorithm 1

A quasi-linear approximation for the TRM boundary condition

Finally, some discussion about the form of $\boldsymbol{N}_{\nu _t,f}$ is necessary. At first glance, it appears that the form of $\nu _t$ in (2.4b ) would need to be updated as the mean velocity changes during the solution procedure in Algorithm1 – indeed, the form of $\nu _t$ originates from the relation $-\overline {u'v'}=\nu _t\, {\rm d}U/{\rm d}y$ . However, it is important to note that $\boldsymbol{N}_{\nu _t,f}$ is not for the mean velocity but for turbulent fluctuations. As such, in principle, there is no strict reason why $\nu _t$ should be changed with the updated $U_r(y)$ . In fact, the most important role of the eddy viscosity in $\boldsymbol{N}_{\nu _t,f}$ is the removal of the fluctuation energy from the integral length scales, mimicking the energy transferred from the integral length scales to the Kolmogorov length scales in real flow (for a detailed discussion, see Holford et al. Reference Holford, Lee and Hwang2024b ). Therefore, the update of $\nu _t$ according to the changes in $U(y)$ through $-\overline {u'v'}=\nu _t\, {\rm d}U/{\rm d}y$ is expected to significantly alter the amount of fluctuation energy removed from the integral length scales; in fact, our attempt to update $\nu _t$ through $-\overline {u'v'}=\nu _t\, {\rm d}U/{\rm d}y$ yields blown-up solutions, which provides additional evidence that the eddy viscosity only acts as a sink term for dissipating fluctuation energy. To keep the amount of fluctuation energy removed from the integral length scales approximately at the same level, it is therefore decided to fix $\nu _t$ even with the TRM boundary conditions. If $\nu _t$ is fixed, then the forcing weight $W_l(k_x,k_z)$ in (2.6) should be kept at the same level. Therefore, $W_l(k_x,k_z)$ is also set to be fixed with the values obtained from the quasi-linear approximations with no-slip boundary conditions. These modelling choices appear to impose the implicit assumption that the flow operates in the smooth-wall-like regime so that the near-wall turbulent dynamics is not modified structurally. Therefore, they are expected to limit the quasi-linear approximations with TRM boundary conditions to be applicable when the slip and transpiration lengths are modest. In practice, the closure appears to work fairly well in a relatively wide range of transpiration and slip lengths, given the comparisons of the roughness functions modelled with the quasi-linear approximations with those from DNS (see § 3).

Table 2.

Summary of the prediction results with various slip ( $l_{x}^+$ , $l_{z}^+$ ) and transpiration ( $m_{x}^+$ , $m_{z}^+$ ) lengths at $Re_{\tau }=180$ . The subscripts $m$ , $2$ and $d$ denote the quantities predicted with MQLA, M2QLA and DQLA, respectively. Here, $Re_{\tau }$ is based on $u_{\tau }$ at $y=0$ , and $Re_{\tau }$ with a prime is based on $u_{\tau }$ at $y=-l_{uv}$ with $h'=h+l_{uv}$ . The virtual origins of the mean flow and the Reynolds shear stress are computed a posteriori. The roughness function here is computed based on the virtual origin framework (i.e. (3.1)). The first seven cases are with isotropic transpiration lengths $m_{x}^+=m_{z}^+$ . The last eight cases are with anisotropic transpiration lengths, with L $\langle \cdot \rangle$ MX $\langle \cdot \rangle$ denoting only streamwise transpiration imposed, whereas L $\langle \cdot \rangle$ MZ $\langle \cdot \rangle$ represents only spanwise transpiration imposed. For direct comparison, the case names here are chosen to be identical to those in Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022).

2.5. Summary of modelling assumptions and their assessment

Thus far, we have developed a quasi-linear framework that can predict turbulence statistics and spectra subject to the TRM boundary conditions. Given that some non-trivial assumptions are made for the formulation, here we provide a summary of these assumptions and their expected effects on the development of the quasi-linear models with the TRM boundary conditions in the present study.

1. (i) Nonlinear term $\boldsymbol{N}_{\nu _t,f}$ in the form of (2.4). We have assumed that the nonlinear terms in the fluctuation equations are modelled as $\boldsymbol{N}_{\nu _t,f} = \boldsymbol{\nabla} \boldsymbol{\cdot} (\nu _t (\boldsymbol{\nabla} \boldsymbol{u}'+ \boldsymbol{\nabla} {\boldsymbol{u}'}^{\rm T}))+\boldsymbol{f}'$ as in (2.4). Given the purpose of the quasi-linear approximations modelling ‘energy-containing scales’, the eddy viscosity model is introduced to mimic the energy transported from the large integral scale to the dissipative small scale (i.e. the energy removal mechanism at large integral scale), and the forcing is given to generate a non-zero solution for the framework. This is a minimal model to represent the energy balance on the integral length scale. For the same purpose of modelling the integral length scale, the leading POD modes from the fluctuation equations are also used. This assumption allows us to incorporate the energy cascade and the driving mechanism of non-trivial velocity fluctuations in a minimal manner. It is a crucial modelling assumption for the given quasi-linear approximation framework, but its physical relevance and validity have been extensively discussed and tested in Holford et al. (Reference Holford, Lee and Hwang2024a ).

2. (ii) Forcing weight $W_{l}(k_x,k_z)$ in (2.11). The streamwise weight $W_{l,k_x}(k_x/k_z)$ is fundamental to account for the self-similar nature of the energy-containing motions, i.e. the attached eddy hypothesis (Townsend Reference Townsend1976; Hwang Reference Hwang2015). In the present study, three different choices for $W_{l,k_x}(k_x/k_z)$ are made to test their influence on the quasi-linear approximations made in the present study. Indeed, the detailed form of $W_{l,k_x}(k_x/k_z)$ does not appear to significantly influence the predictions of the roughness functions, as we will see in § 3 (see table 2, in particular).

3. (iii) Fixed eddy viscosity $\nu _t$ for fluctuation equations. This assumption has been found to be essential for the development of quasi-linear models that are numerically stable, as discussed in detail in § 2.4. From a physical viewpoint, this assumption implies that $\boldsymbol{N}_{\nu _t,f}$ behaves as in the smooth-wall-like regime. As such, this assumption will make the present quasi-linear approximations work presumably in the transitionally rough regime.

3.1. Quantification of roughness effects: the virtual origins

We first introduce how we quantify the roughness effects for the purpose of assessing the performance of the proposed quasi-linear approximations compared to DNS of Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022). A well-known metric for assessing the effects of rough surfaces is the roughness function, denoted as $\Delta U$ . The concept of the roughness function is founded on the existence of the logarithmic law in regions sufficiently distant from the direct influence of surface roughness, analogous to the behaviour observed over smooth walls. In the presence of roughness on a surface, a vertical shift in the mean velocity profile $\Delta U$ within the logarithmic region can be induced (Clauser Reference Clauser1954; Hama Reference Hama1954). This shift, represented by the roughness function, describes the momentum deficit caused by surface roughness. Consequently, the roughness function serves as a measure of drag modification relative to a smooth surface, with $\Delta U \lt 0$ indicating an increase in drag, and $\Delta U \gt 0$ a reduction in drag.

The roughness function can be quantified using the virtual origin framework (e.g. Ibrahim et al. Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021), and this approach is employed here for direct comparison with the previous DNS (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022), which used the same approach to quantify the rougness function. Specifically, it is determined as the difference between the virtual origin of the mean flow ( $l_{U}$ ), representing the point where the mean flow perceives an idealised smooth wall, and the virtual origin of near-wall turbulence ( $l_{T}$ ), where turbulence in the near-wall region interacts with the surface. Luchini (Reference Luchini1996) was the first to reveal a distinction between the virtual origin of the mean flow ( $l_{U}$ ) and that of near-wall turbulence ( $l_{T}$ ). This difference arises because the turbulence in the near-wall region is primarily driven by the quasi-streamwise vortices, which would experience a displacement due to the weakening of the cross-flow shear and impermeability with the existence of rough surfaces. The virtual origin of the mean flow can be established from the slip velocity at the origin of the smooth wall. The virtual origin may be defined for various flow quantities, such as streamwise ( $l_u$ ), wall-normal ( $l_v$ ) and spanwise ( $l_w$ ) velocities, or Reynolds shear stress ( $l_{uv}$ ). Ibrahim et al. (Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021) proposed that the most appropriate choice to define the virtual origin of near-wall turbulence is based on the Reynolds shear stress (i.e. $l_{T}=l_{uv}$ ). Following previous studies (Ibrahim et al. Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021; Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022), the virtual origin of near-wall turbulence here is determined a posteriori by computing the wall-normal shift in the Reynolds shear stress profile. In particular, based on trial and error, the wall-normal shift is sought by aligning the rough-wall solution with the smooth-wall counterpart for the best fit in the near-wall region. With the given definitions of the virtual origins, we note that the TRM model is expected to model only small to moderate roughness effects on the flow.

The roughness function can be predicted based on the virtual origin framework with well-converged values of $U_r^+$ and $\overline {u'v'}_r^+$ (here, the superscript $+$ denotes the viscous inner scaling). In particular, the roughness function is quantitatively determined by the difference between the virtual origin of mean velocity ( $l_U^+$ ) and the virtual origin of Reynolds shear stress ( $l_{uv}^+$ ) (Ibrahim et al. Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021):

(3.1) \begin{equation} \Delta U^+ = l_U^+ - l_{uv}^+ , \end{equation}

where $l_U^+$ is the well-converged slip velocity of $U_r^+$ at the boundary plane ( $y^+ = 0$ ), and $l_{uv}^+$ corresponds to the wall-normal shift of $\overline {u'v'}_{\!\!r}^+$ relative to $\overline {u'v'}_{\!\!s}^+$ , which can give the best fit in the near-wall region (i.e. $10\lt y^+\lt 25$ ) (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022). We note that $\Delta U^+$ produced with this approach shows very little difference from that computed by the conventional method based on the mean velocity difference in the log or outer region.

3.2. Roughness functions

In this subsection, we will mainly focus on analysing the effects of isotropic roughness by computing the virtual origins alongside the roughness functions with various slip and transpiration lengths, i.e. $m_x=m_z$ and $l_x=l_z$ . The effect of anisotropic roughness ( $m_x\neq m_z$ , in particular) is also studied using the present quasi-linear frameworks. However, it has been found that the modelling capability of the quasi-linear frameworks remains essentially the same, and as such, it is presented in Appendix C for brevity.

The virtual origins of mean flow and near-wall turbulence, along with the resulting roughness functions, predicted using the MQLA, M2QLA and DQLA frameworks, are presented in table 2. The predictions from MQLA, M2QLA and DQLA are generally very similar to each other, especially in the cases with mild roughness. In the cases with zero transpiration lengths (i.e. L2M0 and L5M0), drag reduction is achieved; here, all the predicted roughness functions remain positive, with the virtual origin of the mean flow positioned above that of near-wall turbulence, indicating increased momentum relative to the smooth-wall case. This also agrees with many previous numerical findings under purely slip boundary conditions throughout the literature (Min & Kim Reference Min and Kim2004; Busse & Sandham Reference Busse and Sandham2012; Gómez-de-Segura & Garcia-Mayoral Reference Gómez-de-Segura and García-Mayoral2020). In contrast, drag increase is observed, accompanied by a negative roughness function, in all the cases with non-zero transpiration lengths, except L10M10 based on MQLA and M2QLA predictions. In case L10M10 predicted by MQLA and M2QLA, the virtual origin of mean velocity, determined by the streamwise slip length via (2.2), is sufficiently deep that it remains below the virtual origin of near-wall turbulence, leading to a positive $\Delta U$ . On the other hand, the virtual origin predicted by DQLA of mean velocity $l_{U_d}^+$ is slightly higher than values predicted by MQLA and M2QLA, while the DQLA-predicted virtual origin of near-wall turbulence is much deeper, resulting in an increase in drag with a negative roughness function. We note that L10M10 is the case most poorly predicted by the present quasi-linear approximations due to the largest values of the slip and transpiration lengths (see discussion below).

In comparison to the cases with zero transpiration length, transpiration length has a non-trivial impact on the position of the virtual origin of near-wall turbulence, which pushes it deeper. When the virtual origin of near-wall turbulence is positioned deeper than that of the mean velocity, an increase in drag is observed. This phenomenon is a result of the non-zero transpiration lengths enabling the displacement of near-wall quasi-streamwise vortices towards the boundary plane, which can result in the generation of wall-normal velocity fluctuations (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022). However, the slip lengths applied to the wall-parallel velocity components affect not only the virtual origin of mean velocity but also that of near-wall turbulence. Comparisons among the cases with the same transpiration lengths (e.g. L2M0 and L5M0, L2M5 and L5M5, as well as L5M10 and L10M10) clearly reveal that the increased slip lengths lead to deeper virtual origins for both mean velocity and near-wall turbulence. The effect of slip lengths on the virtual origin of mean velocity is straightforward (see (2.2)), while their impact on virtual origin of near-wall turbulence can be understood through a reformulation of the TRM boundary condition for the wall-normal velocity (2.3c ) into

(3.2) \begin{equation} {v'} =-m_xl_x\frac {\partial ^2 {u'}}{\partial x\,\partial y} \left \vert _{y=0,2h} -m_zl_z\frac {\partial ^2 {w'}}{\partial z\,\partial y}\right \vert _{y=0,2h}, \end{equation}

via substituting the slip boundary conditions of wall-parallel velocity components (2.3a ) and (2.3b ) into (2.3c ). Evidently, the slip lengths are capable of modifying the wall-normal velocity, thereby exerting further non-trivial impact on the Reynolds shear stress, which leads to the displacement of the virtual origin of near-wall turbulence.

Table 3.

Comparison between the M2QLA predictions and the corresponding DNS data from Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022) with various slip ( $l_{x}^+$ , $l_{z}^+$ ) and transpiration ( $m_{x}^+$ , $m_{z}^+$ ) lengths. Here, the subscript ${D}$ denotes the DNS data. The virtual origins of the velocity fluctuations $l_{u}^+$ , $l_{v}^+$ and $l_{w}^+$ are also computed a posteriori, via extrapolation of their root mean square profiles.

Given the general similarity among the predictions with MQLA, M2QLA and DQLA (see also Appendix B for a further comparison), we select the prediction results of M2QLA to compare with the DNS data from Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022). In Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022), the DNS were performed in a computational domain with $L_x=2\pi h$ , $L_y=2 h$ and $L_z=\pi h$ , with the numbers of grid points $N_x=196$ , $N_y=144$ and $N_z=160$ at $Re_{\tau }=180$ . The flow is driven by an imposed constant mean pressure gradient. The periodic boundary conditions are applied along the streamwise and spanwise directions, whereas the TRM boundary conditions are imposed on the wall-normal domain boundaries at $y = [0, 2h]$ to model the surface roughness. As shown in table 3, indeed, the M2QLA predictions align very well with the DNS results, especially for virtual origins and roughness functions in the mild isotropic roughness cases. The discrepancy with the DNS data is gradually enlarged with increasing slip and transpiration lengths, as expected. Specifically, in the cases with large isotropic roughness (e.g. L10M10), the M2QLA predictions provide deeper virtual origins of streamwise and spanwise velocities than those from DNS. In contrast, the M2QLA-predicted virtual origin of wall-normal velocity is shallower. This appears to be the result of the overestimation of the streamwise and spanwise turbulence intensities, while underestimating the wall-normal intensity by M2QLA under the large TRM boundary conditions. Furthermore, the virtual origins of the mean velocity from DNS exhibit a pronounced reduction from the corresponding streamwise slip lengths in the cases with large roughness (i.e. $l_{U_D}^+\lt l_{x}^+$ ; see L10M10, for example). However, this trend is not well captured by the M2QLA predictions. Indeed, this discrepancy persists in the predictions from all of MQLA, M2QLA and DQLA. Similarly, the Reynolds shear stress profiles predicted by M2QLA are greater than those from DNS, leading to larger mean velocities predictions, and contributing to the discrepancy of the roughness functions.

3.3. One-point turbulence statistics

Figure 3.

(a,c) Mean velocity and (b,d) Reynolds shear stress profiles of the cases L2M2 (red), L2M5 (green), L5M5 (blue), L5M10 (magenta), L10M10 (cyan) and smooth wall (black), with (a,b) M2QLA and (c,d) DNS (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022) at $Re_{\tau }=180$ .

Considering M2QLA providing better statistical predictions than MQLA with a computational cost much smaller than DQLA, we present the predicted one-point statistics with M2QLA in this subsection – note also that MQLA cannot describe the effect of $m_x$ by construction, as seen from (2.3). Figure 3(a) illustrates the mean velocity profiles obtained with various slip and transpiration lengths for the isotropic roughness cases. In the region very close to the wall, the mean velocities are seen to be heavily influenced by the imposed streamwise slip boundary conditions. In this region, indeed, the mean velocity profiles exhibit consistent behaviours across the cases with the same streamwise slip lengths (i.e. the cases L2M2 and L2M5, as well as L5M5 and L5M10). Additionally, increasing the streamwise slip length leads to the greater mean velocity in this region. In the region further away from the wall, most mean velocity profiles align closely with that of a smooth wall, with the exception of the case L5M10, whose mean velocity profile exhibits a clear downward shift in the channel centre, which agrees well with DNS observations (see figure 3 c). However, there exists a discrepancy between the mean velocity predicted by M2QLA and that of DNS for the case L10M10. In this case, the mean velocity predicted by M2QLA in the region very close to the wall is much higher than the DNS solution, indicating that the virtual origin of the mean velocity from the M2QLA prediction is much deeper than that from DNS. This also leads to an evident inconsistency in the region further away from the wall, where the predicted mean velocity of M2QLA is always greater than the smooth-wall solution resulting in a positive roughness function, while the DNS data exhibit a downward shift from the smooth-wall case in the centre of the channel, with a negative roughness function suggesting an increase in drag (see figures 3(a,c) and table 3).

As presented in figure 3(b), the non-zero transpiration length can produce finite-amplitude Reynolds shear stress at the boundary plane. Larger transpiration lengths, in particular, can induce stronger Reynolds shear stresses in the region close to the wall, which effectively shifts the virtual origin of near-wall turbulence downwards. This effect may counteract the slip of mean velocity, leading to the increase in drag relative to the smooth-wall scenario (e.g. $l_{uv}^+\approx 3.5$ and $\Delta {U}^+\approx -1.5$ for L2M5). In the mild roughness cases, the predicted Reynolds shear stresses align well with the DNS data, whose peaks move towards the boundary plane and become higher with increasing slip and transpiration lengths (see figure 3 b,d). However, this is not the case when the transpiration and slip lengths are large. For example, the M2QLA-predicted Reynolds shear stress of the L10M10 case exhibits a clear increase from that of the L5M10 case, which is inconsistent with the DNS observations. Furthermore, the Reynolds shear stress is overestimated compared to its DNS counterpart in the near-wall region, owing to the existence of larger mean velocity there. In the region farther away from the wall, all the Reynolds shear stress profiles collapse to the smooth-wall profile, which agree qualitatively with the previous DNS findings in Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022).

Figure 4.

Cross-streamwise view of the leading POD mode at $(\lambda _x^+,\lambda _z^+)\approx (\infty ,100)$ in the cases (a) smooth wall, (b) L5M5, (c) L10M10 with M2QLA at $Re_{\tau }=180$ . The contours denote streamwise velocity, and the vectors represent cross-streamwise velocities. The velocity field of each POD mode is normalised by its maximum streamwise velocity.

To understand the observations in figure 3 from a physical viewpoint, the related leading POD mode with $(\lambda _x^+,\lambda _z^+)\approx (\infty ,100)$ ( $(k_x^+,k_z^+)=(0,0.063)$ ) is visualised in figure 4. In the smooth-wall case, the leading POD mode appears in the form of streamwise streaks accompanied by (weak) counter-rotating streamwise vortices (Hwang & Cossu Reference Hwang and Cossu2010). In case L5M5, where the predictions of the roughness functions are reasonably good, the streaks are only slightly shifted towards the wall, and the structures around the peak do not seem to be largely affected. However, for the case L10M10, a significant part of these streaks is ingested into the boundary plane, indicating that the flow structures modelled by M2QLA are heavily affected. This explains why the predictions of the quasi-linear approximations for this become deteriorated – the assumption that turbulence remains smooth-wall-like made for the nonlinear term model (assumption (iii) in § 2.5) is presumably no more valid for the large slip and transpiration lengths, resulting in the relatively poor performance of the quasi-linear approximations.

Figure 5.

(a,c) Mean velocity and (b,d) Reynolds shear stress profiles of the cases L2M2 (red), L2M5 (green), L5M5 (blue), L5M10 (magenta), L10M10 (cyan) and smooth wall (black) from (a,b) M2QLA and (c,d) DNS data (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022).

Then the impact of surface roughness on friction velocity $u_{\tau }$ is analysed. Considering that the near-wall turbulence perceives a virtual smooth wall at $y=-l_{uv}$ , the friction velocity $u_{\tau }$ is expected to correspond to the shear stress at this position. Furthermore, in a channel flow, the total stress varies linearly with the wall-normal coordinate. As such, the friction velocity at $y=-l_{uv}$ can be determined by extrapolating the total stress curve from $y=0$ (Gómez-de-Segura & Garcia-Mayoral Reference Gómez-de-Segura, Fairhall, MacDonald, Chung and Garcia-Mayoral2020), which leads to

(3.3) \begin{equation} u_{\tau }|_{y=-l_{uv}} =u_{\tau }|_{y=0}\sqrt {\frac {h+l_{uv}}{h}}. \end{equation}

This friction velocity is subsequently used to rescale flow quantities, and the resulting mean velocity and Reynolds shear stress profiles with the origin set at $y^+=-l_{uv}^+$ are shown in figures 5(a) and 5(b). Indeed, the mean velocity and Reynolds shear stress from the majority of the cases considered so far align closely with the smooth-wall solution. On the other hand, noticeable deviations from the smooth-wall solution can be observed in cases L5M10 and L10M10 with large roughness. Specifically, the cases L5M10 and L10M10 exhibit higher mean velocities with stronger Reynolds shear stresses relative to the smooth-wall profile in the near-wall region. Compared to DNS observations (figure 5 c,d), the deviations from the smooth-wall solution in figures 5(a) and 5(b) are more pronounced, confirming that M2QLA overestimates the alteration of near-wall turbulence for these two cases.

Figure 6.

(a,d) Streamwise, (b,e) wall-normal and (c,f) spanwise turbulent intensity profiles of L2M2 (red), L2M5 (green), L5M5 (blue), L5M10 (magenta), L10M10 (cyan) and smooth wall (black) cases, with (a,b,c) M2QLA and (d,e,f) DNS (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022) at $Re_{\tau }=180$ .

Turbulent intensities predicted by M2QLA are also shown in figure 6(a–c). At the boundary plane, all the turbulence intensities exhibit finite amplitudes with non-zero transpiration lengths. The amplitudes at the boundary plane also increase with larger slip and transpiration lengths. In the near-wall region, compared to the DNS results in figure 6(d–f), qualitative agreement in the streamwise turbulent intensities is evident for the cases with mild roughness. In particular, the peaks of the streamwise turbulent intensities in these cases are all higher than that of the smooth-wall case, and are achieved in the region closer to the boundary plane. However, in cases L5M10 and L10M10, the DNS results exhibit a damping effect that reduces the streamwise turbulent intensities below the smooth-wall solution, while the M2QLA predictions are still higher than the smooth-wall profile (figure 6 a,d). On the other hand, the wall-normal and spanwise turbulent intensities predicted by M2QLA show qualitative agreement with the DNS results for all the considered cases, as they increase with larger slip and transpiration lengths in the near-wall region (figures 6(b,c) and 6(e,f)). Given the overestimation of streamwise intensities in the cases L5M10 and L10M10, the turbulence modification towards isotropisation in M2QLA appears not to be as obvious as that in DNS (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022).

The one-point turbulent statistics are subsequently rescaled using the friction velocity at $y^+=-l_{uv}^+$ , and the results are presented in figure 7, with the origin set at $y^+=-l_{uv}^+$ . Notably, in cases L2M2, L2M5 and L5M5, the predicted streamwise intensities are slightly higher relative to the smooth-wall case (figure 7 a), which aligns with the findings of Ibrahim et al. (Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021) and Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022) (see also figure 7 d). The wall-normal intensities in these cases remain approximately the same as those of the smooth-wall solution (figure 7 b). These behaviours are consistent with DNS results (figure 7 e). However, there are also some discrepancies between the spanwise intensities predicted by M2QLA and those from DNS. Specifically, the M2QLA predictions of spanwise intensities in these cases show slight increases (figure 7 c). However, DNS data reveal that the spanwise intensities collapse well with the smooth-wall solution (figure 7 f). In general, the discrepancies between the M2QLA predictions for turbulence intensities and the DNS data tend to enlarge with increasing slip and transpiration lengths. Notable discrepancies with DNS results are observed in case L10M10: for example, a significant decrease in streamwise intensity and a considerable increase in wall-normal intensity from those of the smooth-wall case were observed in DNS (figure 7 d,e), but these are not captured in the M2QLA.

Figure 7.

(a,d) Streamwise, (b,e) wall-normal and (c,f) spanwise turbulent intensity profiles of L2M2 (red), L2M5 (green), L5M5 (blue), L5M10 (magenta), L10M10 (cyan) and smooth wall (black) cases, based on (a,b,c) M2QLA and (d,e,f) DNS data (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022) with the origin set at $y^+=-l_{uv}^+$ and rescaled with the $u_{\tau }$ value at that plane.

3.4. The 2-D spectra

Figure 8.

Pre-multiplied 2-D spectral densities of (a,e,i,m) $u^2$ , (b,f,j,n) $v^2$ , (c,g,k,o) $w^2$ and (d,h,l,p) $uv$ of cases (a–d) L2M2, (e–h) L2M5, (i–l) L5M5 and (m–p) L10M10. The shaded regions correspond to the smooth-wall solution at $y^+\approx 15$ , while the solid lines are the TRM cases at $y^++l_{uv}^+\approx 15$ scaled using $u_{\tau }$ at $y^+=-l_{uv}^+$ . The contour increments for each column are 0.3, 0.01, 0.06, 0.03, respectively.

The 2-D premultiplied velocity spectra are also examined at a fixed wall-normal location for selective cases. The MQLA and M2QLA frameworks, which include zero and only one streamwise dependent mode, respectively, lack the capability to capture the structural characteristics through 2-D energy spectra. As such, the DQLA framework is employed here, including the non-zero streamwise weight. Following previous work (Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022), we present the 2-D spectra at $y^++l_{uv}^+\approx 15$ for the cases L2M2, L2M5, L5M5 and L10M10 (from top to bottom) in figure 8. The spectra for $u^2$ and $uv$ over smooth walls illustrate a strong qualitative agreement with the DNS results, where DQLA effectively replicates all essential flow features by its construction: for example, the peak locations in the $u^2$ and $uv$ spectra of DNS appear in $(\lambda _x^+, \lambda _z^+)\approx (1000,100)$ and $(\lambda _x^+, \lambda _z^+)\approx (600,80)$ , respectively (figure 5 in Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022), well-captured by the DQLA predictions. However, as in figure 2(b,c), some discrepancies with DNS data are observed in the $v^2$ and $w^2$ spectra. In particular, compared to DNS results from Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022) (see their figure 5), the DQLA-predicted primary peaks in the spectra of $v^2$ and $w^2$ over smooth walls are achieved with shorter streamwise wavelengths and longer spanwise wavelengths. These differences in the smooth-wall case mainly originate from the nonlinear term model in the DQLA that do not capture the nonlinear physical processes associated with the pressure-strain transport (i.e. streak instability and transient growth), as discussed in detail in Holford et al. (Reference Holford, Lee and Hwang2024a ,Reference Holford, Lee and Hwang b ).

Figure 9.

Flow visualisation of wall-normal velocity fluctuation of the Fourier mode associated with the Kelvin–Helmholtz instability ( $\lambda _x^+=119$ and $\lambda _z^+=182$ ) in (a,b) the $x{-}y$ plane at $z^+=0$ , and (c,d) the $x{-}z$ plane at $y^++l_{uv}^+\approx 15$ , with (a,c) the smooth-wall case, (b,d) the L10M10 case at $Re_{\tau }=180$ . Here, the contour levels are directly from the spectra without amplitude modifications of the mode.

When the TRM boundary conditions are applied, like the DNS results, the 2-D spectra for the L2M2 and L2M5 cases differ little from the smooth-wall spectra. On the contrary, the spectra for the L5M5 case exhibit some notable deviations in $u^2$ , $v^2$ , $w^2$ from the smooth-wall case (figure 8 i–k), although the Reynolds shear stress co-spectra (figure 8 l) remain largely smooth-wall-like. The L10M10 case presents strong deviations from the smooth-wall case in all spectra. Indeed, the near-wall spectra at small streamwise and large spanwise wavelengths become more energetic under the TRM boundary conditions with large slip and transpiration lengths, which agrees with DNS observations (figure 5 in Khorasani et al. Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022). Although weak (see also a discussion in the last paragraph of § 5), the deviations observed in the spectra of the wall-normal velocity fluctuations for $60 \leqslant \lambda ^+_x \leqslant 200$ , and for the spanwise wavelengths $\lambda ^+_z \geqslant 130$ , under the TRM boundary conditions appear to be reminiscent of those associated with the initiation of a Kelvin–Helmholtz instability, which has been observed in flows over large riblets with certain shapes (Garcia-Mayoral & Jiménez Reference García-Mayoral, Gómez-de-Segura and Fairhall2011; Endrikat et al. Reference Endrikat, Modesti, García-Mayoral, Hutchins and Chung2021a ,Reference Endrikat, Modesti, MacDonald, García-Mayoral, Hutchins and Chung b ) and highly permeable porous structures (Jimenez et al. Reference Jimenez, Uhlmann, Pinelli and Kawahara2001; Kuwata & Suga Reference Kuwata and Suga2017; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). For providing additional evidence of initiating a Kelvin–Helmholtz instability in the L10M10 case, wall-normal velocity fluctuations of the Fourier mode with $\lambda _x^+=119$ and $\lambda _z^+=182$ are visualised. Indeed, the wall-normal velocity fluctuations in the L10M10 case peak in the near-wall region and extend to the wall, whereas those in the smooth-wall case achieve their peaks at $y^+\approx 80$ and extend to the channel centre (figure 9 a,b). Furthermore, the wall-normal velocity fluctuations at $y^++l_{uv}^+\approx 15$ in the L10M10 case become largely stronger than those over a smooth wall (note that the amplitude of this mode is very weak in the smooth-wall case; see figure 8), suggesting that the fluctuations with small streamwise and large spanwise wavelengths are indeed initiated to be more energetic with large slip and transpiration lengths.