Modelling spanwise heterogeneous roughness through a parametric forcing approach

Abstract

Inhomogeneous rough surfaces in which strips of roughness alternate with smooth-wall strips are known to generate large-scale secondary motions. Those secondary motions are strongest if the strip width is of the order of the half-channel height and they generate a spatial wall shear stress distribution whose mean value can significantly exceed the area-averaged mean value of a homogeneously smooth and rough surface. In the present paper it is shown that a parametric forcing approach (Busse & Sandham, J. Fluid Mech., vol. 712, 2012, pp. 169–202; Forooghi et al., Intl J. Heat Fluid Flow, vol. 71, 2018, pp. 200–209), calibrated with data from turbulent channel flows over homogeneous roughness, can capture the topological features of the secondary motion over protruding and recessed roughness strips (Stroh et al., J. Fluid Mech., vol. 885, 2020, R5). However, the results suggest that the parametric forcing approach roughness model induces a slightly larger wall offset when applied to the present heterogeneous rough-wall conditions. Contrary to roughness-resolving simulations, where a significantly higher resolution is required to capture roughness geometry, the parametric forcing approach can be applied with usual smooth-wall direct numerical simulation resolution resulting in less computationally expensive simulations for the study of localized roughness effects. Such roughness model simulations are employed to systematically investigate the effect of the relative roughness protrusion on the physical mechanism of secondary flow formation and the related drag increase. It is found that strong secondary motions present over spanwise heterogeneous roughness with geometrical height difference generally lead to a drag increase. However, the physical mechanism guiding the secondary flow formation, and the resulting secondary flow topology, is different for protruding roughness strips and recessed roughness strips separated by protruding smooth surface strips.

1. Introduction

Turbulent flows over spatially inhomogeneous rough surfaces are present in a variety of engineering and environmental applications, whereas smooth surfaces are rather the exception. In recent years, spanwise inhomogeneous rough surfaces have attracted increasing attention in research. Such surfaces are known to generate secondary flows of Prandtl's second kind (Prandtl 1931; Hinze 1967; Anderson et al. 2015) that appear as large-scale vortical structures in the time-averaged flow field perpendicular to the main flow direction (Hinze 1967; Barros & Christensen 2014; Vanderwel et al. 2019). These secondary flows amount to a few per cent of the main flow's energy only, but change the flow properties of the main flow, which in turn leads to a considerable modification of friction and heat transfer (Stroh et al. 2020a) depending on the width of the rough surface strip (Chung, Monty & Hutchins 2018; Wangsawijaya et al. 2020).

Within spanwise heterogeneous rough surfaces two types of surfaces are often distinguished in the literature, the so-called ridge-type and strip-type roughnesses (Wang & Cheng 2006). The former surface is characterized by significant differences in the spanwise wall elevation, e.g. streamwise aligned ridges or bars, which produce large secondary motions with upward motion above the elevated region. This type is studied in experiments (e.g. Vanderwel & Ganapathisubramani 2015; Medjnoun, Vanderwel & Ganapathisubramani 2018; Zampiron, Cameron & Nikora 2020) and through various numerical simulations (e.g. Wang & Cheng 2006; Hwang & Lee 2018; Stroh et al. 2020a). Strip-type roughness is characterized by spanwise alternating wall shear stress conditions, where spanwise wall elevations are not present or negligibly small. These surfaces produce large-scale secondary motions with upward motion above the lower stress patch and downward motion of the circulation above the high stress region if the spanwise wavelength of the surface structure is of the order of the boundary layer thickness (Chung et al. 2018). In numerical simulations idealized strip-type roughness conditions can be generated by alternating wall-normal velocity gradients at the wall (e.g. Willingham et al. 2014; Chung et al. 2018) while in experiments the smooth-wall region in between roughness strips can be elevated with the aim of avoiding offsets in the virtual origin between rough and smooth surfaces (Wangsawijaya et al. 2020).

While results documenting the behaviour of secondary motions are often described in a time-averaged and, in the case of direct numerical simulation (DNS), additionally often also in a phase-averaged sense, Wangsawijaya et al. (2020) point out that this averaging procedure actually masks a time-dependent meandering behaviour of the secondary motions that they found to be largest in the case of roughness strips with widths of the order of the boundary layer thickness.

Recently, Stroh et al. (2020b) demonstrated that the rotational direction of secondary motions, considered from an averaged point of view, above numerically resolved rough-wall strips depends on the relative surface elevation (protruding or recessed roughness) such that a secondary motion behaviour corresponding to either strip- or ridge-type roughness can be produced. Capturing the switch between the strip- and ridge-type behaviour was possible due to the fact that the roughness was resolved in the DNS of Stroh et al. (2020b). This cannot be achieved if one employs idealized boundary conditions or wall functions, e.g. those used in the DNS by Chung et al. (2018) or the large-eddy simulation (LES) approach by Willingham et al. (2014). However, roughness-resolving DNS can be highly demanding in terms of computational cost; what serves as a motivation to identify alternative numerical approaches that are able to capture this effect.

In the present study, we show that a parametric forcing approach (PFA) (Busse & Sandham 2012; Forooghi et al. 2018), in which the roughness effects are introduced into the momentum equations by a volume force, is able to capture the general drag increase for spanwise inhomogeneous roughness with spanwise wavelength of the order of the boundary layer thickness and, in particular, the different averaged secondary flow behaviours above protruding and recessed roughness strips as reported by Stroh et al. (2020b). While the direct representation of a rough surface (through boundary fitted grids or immersed boundary method) typically requires increased spatial resolution in a DNS, the application of a PFA roughness model provides the advantage that the resolution requirements for flow simulations over rough surfaces remain identical to the ones of smooth-wall DNS.

2. Methodology

The DNS are performed analogous to the flow configuration in Stroh et al. (2020b), where a fully developed turbulent open channel flow is driven by a prescribed pressure gradient adjusted to maintain a constant friction Reynolds number of $Re_\tau = 500$. The friction Reynolds number is $Re_\tau = u_\tau \delta _{eff} /\nu$, with the friction velocity $u_\tau$, the effective half-channel height $\delta _{eff}$, whose determination will be described in detail below, and the kinematic viscosity $\nu$. The streamwise, wall-normal and spanwise directions are denoted by $x$, $y$, $z$, respectively. The bottom wall consists of spanwise alternating strips of smooth and rough surface walls with varying relative smooth-wall elevations. The three different rough and smooth surface configurations from Stroh et al. (2020b) are depicted in figure 1($a$–$c$), the protruding roughness case $(a)$, the intermediate case $(b)$ and the recessed roughness case $(c)$.

Figure 1. Schematic of the open channel domain with roughness strips at the wall. Variation of the smooth-wall elevation with fully resolved roughness (cases $a$–$c$) and modelled roughness (cases $d$–$f$). Cases $(a$,$d)$ correspond to protruding roughness, $(b$,$e)$ to an intermediate roughness and $(c$,$f)$ to a recessed roughness configuration.

The structured surfaces in Stroh et al. (2020b), as well as the smooth-wall elevation, are numerically resolved by the immersed boundary method (IBM) of Goldstein, Handler & Sirovich (1993). The spanwise wavelength, L, of one pair of smooth- and rough-wall patches is fixed to $L/\delta = 1$ and the ratio of the smooth-wall width $W$ to the wavelength is set to $W/L = 0.5$. The striped surface texture is generated by distributing randomly several discrete roughness elements and deleting those elements whose roughness centre position is placed in a smooth-wall strip, while fulfilling the prescribed roughness statistics. Thus, the individual roughness strips are not identical in their detailed topography. Moreover, foothills of roughness elements placed closely to the rough–smooth border can protrude slightly into the smooth-wall strip. For the homogeneous rough surface the mean elevation of the rough surface is $\bar {k} = 0.043\delta$, the maximum roughness height $k_{max} = 0.1 \delta$ and the root mean square of the roughness height distribution is $k_{rms}= 0.024\delta$. The generation procedure of the rough-wall strips is based on the same statistical roughness properties as Stroh et al. (2020b) and further details on the properties of the resolved roughness can be found in this reference. For the present study the roughness patch is modelled by the PFA roughness model (Forooghi et al. 2018), while the smooth-wall elevation is represented by the same IBM. Since the PFA, as explained below, does not resolve local roughness effects and is horizontally homogeneous within the rough region, the rough strips are all identical to each other in this case.

The governing equations for the incompressible turbulent flow

(2.1a,b)\begin{equation} \frac{\partial u_i}{\partial x_i} = 0, \quad \frac{\partial u_i}{\partial t} + \frac{\partial u_i u_j}{\partial x_j} ={-} \frac{1}{\rho}\frac{\partial p}{\partial x_i} + \nu \frac{\partial^2 u_i}{\partial x_j \partial x_j} + \varPi \delta_{i1} + F_{{IBM,i}} + F_{{r,i}} \end{equation}

are solved using the spectral solver SIMSON (Chevalier et al. 2007), where $\rho$ is the density, p is the pressure and $\varPi$ is the forcing term required for constant pressure gradient simulations. The three velocity components are denoted by $(u_1, u_2, u_3) = (u,v,w)$. The term $F_{{IBM,i}}$ represents the IBM volume force while $F_{{r,i}}$ is the PFA forcing term for modelling the rough surface.

The introduction of the rough surface strips and the elevation of the smooth-wall strips on top of the bottom wall ($y=0$) leads to a local surface elevation and in consequence to a reduction of the effective cross-section seen by the fluid, while the numerical domain size is kept constant. The effective half-channel height $\delta _{eff}$ is obtained by subtracting the global melt-down height $h_{eff}$ from the constant half-channel height $\delta$, such that $\delta _{eff} = \delta - h_{eff}$. Here, $h_{eff}$ is obtained by averaging the surface elevation over the full channel length and width, including roughness and elevated smooth-wall regions. In order to maintain the same friction Reynolds number among all configurations, the effective half-channel height is taken into account for the adjustment of the pressure gradient $\varPi$ (Stroh et al. 2020b).

The PFA forcing term $F_{{r,i}}$ consists of the sum of a linear and a quadratic contribution of the form

(2.2)\begin{equation} F_{r,i} ={-} A(y) u_i - B(y) u_i |u_i|. \end{equation}

The general idea behind the derivation of the two model functions – $A(y)$ and $B(y)$ – is briefly presented in the following while full details can be found in Forooghi et al. (2018). The first and second terms on the right-hand side of (2.2) aim to reproduce the viscous drag and form drag per unit volume caused by roughness elements at a certain wall distance $y$, respectively. An analogy between roughness and porous media is employed to derive an expression for the first function, based on the Kozney–Carman porous medium permeability model

(2.3)\begin{equation} A(y) = k_K\frac{\nu s(y)^2}{\epsilon(y)^3}. \end{equation}

In (2.3) $\epsilon$ (porosity) is the fluid volume per unit total volume, $s$ is the total surface area of the roughness per unit total volume and $k_k$ is an empirical constant. In order to find the values of $s$ and $\epsilon$ as functions of $y$, the roughness perimeter (area) resulting from intersection of $y$-planes with the roughness surface (volume) is used. The function $B(y)$ is derived such that the corresponding term represents the form drag due to all roughness elements. That is

(2.4)\begin{equation} B(y) = c_D\frac{ s_f(y)}{2}, \end{equation}

where $s_f$ denotes the total ‘windward-projected’ surface area of roughness per total volume and $c_D$ is the effective drag coefficient of the roughness. One should keep in mind that the three functions $s_f(y)$, $s(y)$ and $\epsilon (y)$ are uniquely determined based on the specific roughness geometry. They can also be considered as statistical representations of a roughness topography. The two constants $k_K$ and $c_D$ serve as model constants, which enable tuning of the model. The values of these constants have been tuned by Forooghi et al. (2018) based on a number of DNS cases with homogeneous roughness with systematically varied topographies. In the present work, we slightly readjust the constants in order to reproduce the mean velocity profile for the specific roughness topography under consideration as closely as possible the homogeneous rough case. Note that the PFA forcing is not applied in the wall-normal direction, i.e. $F_{r,2} = 0$ following the suggestion by Busse & Sandham (2012).

Obviously, $A(y)$ and $B(y)$ are zero for $y>k_{max} = 0.1\delta$ and their specific distributions, as shown in figure 2, are restricted to $y\leq k_{max} = 0.1\delta = 2.3 \bar {k}$ as there is no roughness above this height. Close to the wall, porosity approaches zero, leading for $A(y)$ to assume very large values. Therefore, as visible in figure 2$(a)$, $A(y)$ is bounded near the wall to ensure numerical stability. This has a negligible effect on the flow since the mean velocity is very small for $y<\bar {k}/2$.

Figure 2. The PFA model functions $A$ and $B$ over the wall-normal distance $y$ normalized using the kinematic viscosity $\nu$ and the mean roughness height $\bar {k}$.

The domain length in the streamwise, wall-normal and spanwise directions is $(L_x, L_y, L_z) = (8\delta, \delta, 4\delta )$ and the same boundary conditions as in Stroh et al. (2020b) are employed. The modelled roughness cases employ two sets of grid resolutions, one being the same as for the resolved roughness cases with a grid resolution $(N_x, N_y, N_z) = (768, 301, 384)$ and a second set with lower resolution in the streamwise and wall-normal directions with $(N_x, N_y, N_z) = (384, 201, 384)$. This allows us to confirm that the grid resolution for the case of the PFA roughness model can be reduced to standard DNS resolution without significant impact on the results.

In order to obtain statistically converged results, statistics were time integrated over a period of at least 50 flow-through time units. Since the smooth- and rough-wall patches cause a spanwise heterogeneity, the velocity field is decomposed into its streamwise- and time-averaged mean and related random fluctuation defined by $u_i(x,y,z,t) = \bar {u}_i(y,z) + u_i'(x,y,z,t)$, where the overbar $\overline {(\cdot )}$ indicates temporal and streamwise averaging. This mean velocity can be decomposed into its global mean and the related spanwise deviation $\bar {u}_i(y,z) = \langle \bar {u}_i \rangle (y) + \tilde {u}(y,z)$, where angular brackets $\langle \cdot \rangle$ indicate spanwise averaging. For the present data sets the spatial averages in wall-parallel planes are obtained through extrinsic averaging, which includes the solid region with zero velocity values.

3. Results

3.1. Global flow properties

The PFA model was originally developed to represent the effect of homogeneous rough surfaces on a turbulent flow field and is applied to spanwise heterogeneous roughness in the present investigation. In order to match the skin friction coefficient of a particular homogeneous rough surface with high accuracy, the coefficients in the model functions $A$ and $B$ of (2.2) require a fine tuning which yields the particular distribution of $A$ and $B$ shown in figure 2 for the present homogeneous rough reference surface. All results obtained in terms of global flow properties are summarized in table 1. These include the smooth-wall and homogeneous rough-wall (resolved and modelled) reference cases, the three resolved roughness strips of Stroh et al. (2020b) with smooth-wall elevations of $h=0$, $h=0.97 \bar {k}$ and $h=1.70 \bar {k}$ and a number of PFA model cases with additional $h$-values in the range $0.25\bar {k} \leq h \leq 2.50\bar {k}$. Note that the definition of $\bar {k}$ is based on the homogeneous rough reference in the present study while $\bar {k}$ is defined separately for each type of inhomogeneous roughness in Stroh et al. (2020b). This different definition induces slight differences of the order of $\Delta \bar {k} / \delta < 0.002$ due to the limited computational box size.

Table 1. Global flow properties from DNS of the homogeneous smooth- and rough-wall cases and the heterogeneous smooth–rough cases with varying smooth-wall distance $h$.

The reduced cross-sectional area of the channel – through the introduction of IBM-based roughness elements or the PFA model – is taken into account through the effective half-channel height $\delta _{eff}$ introduced in § 2. Therefore, the bulk velocity $U_b$ is evaluated as

(3.1)\begin{equation} U_b= \frac{1}{\delta_{eff}}\int_0^\delta \langle \bar{u} \rangle (y)\,\mathrm{d}y. \end{equation}

Since the simulations are run at constant $Re_\tau$ the introduction of roughness leads to a reduced bulk Reynolds number $Re_b = U_b \delta _{eff} / \nu$ and a decrease of the normalized bulk velocity $U_b^+$ (where the superscript + represents viscous units obtained through a normalization with the friction velocity $u_\tau$). Note that $u_\tau$ is defined through the effective wall shear stress which is obtained through an extrapolation of the linear total shear stress distribution to the location $y_0=\delta - \delta _{eff}$ (Chan-Braun, Garcća-Villalba & Uhlmann 2011). The reduced flow rate in the rough channel at constant $Re_\tau$ translates into an increased friction coefficient $c_f = 2u_\tau ^2/U_b^2$ compared with the smooth-wall reference. As reported in Forooghi et al. (2018) the tuned PFA model captures the drag increase in terms of $c_f/c_{f,s}$ (where $c_{f,s}$ represents the skin friction drag coefficient of a smooth wall at the same $Re_\tau$) very well for the homogeneous rough surface. In addition, it can be seen that the reduction of the spatial resolution to standard DNS dimensions for the PFA model does not influence the global flow parameters.

For inhomogeneous roughness the effective (i.e. spanwise-averaged) skin friction coefficient can strongly differ from the area-averaged mean value of homogeneous rough and smooth walls (Türk et al. 2014; Chung et al. 2018; Stroh et al. 2020b; Wangsawijaya et al. 2020). For the present heterogeneous rough case – in which the area-averaged mean value of the smooth wall and the homogeneous rough case corresponds to $c_f/c_{f,s}=2.00$ – the results summarized in table 1 clearly indicate a significant relative drag increase for all considered types of roughness strips. For the resolved roughness the largest relative drag increase is predicted for $h=0$ ($23.5\,\%$) along with the strongest reduction of $Re_b$.

The general drag increasing impact of the inhomogeneous roughness on $c_f$ above the area-averaged value is well captured by the PFA model. However, we find an underprediction of $c_f/c_{f,s}$ for $h=0$ and an overprediction for $h=1.70\bar {k}$, while there is a very good match for $h \approx \bar {k}$. A reduction of the resolution for the modelled roughness to the one of a standard smooth-wall DNS has a negligible to small effect on the obtained results, which is largest for $h=1.70\bar {k}$ with a difference of less than $2\,\%$ for $c_f/c_{f,s}$. Therefore, the additional variations of $h$ investigated with the PFA model only are simulated with low resolution. All results presented in the following refer to the low resolution configuration. Figure 3 shows the streamwise mean velocity profiles (averaged in time and the two wall-parallel spatial directions) in logarithmic scaling. The zero wall location is placed at a distance of $\delta$ from the channel centre for all cases. Figure 3$(a)$ contains the velocity profiles for the resolved and modelled homogeneous rough surface for reference. Slight deviations in the region in the range $\bar {k}< y< k_{max}$ are visible that can be traced back to a similar but not identical wall-normal force distribution in case of IBM and PFA. Overall, very good agreement between these two approaches is obtained. This is in agreement with the very similar integral $U_b^+$ values reported in table 1. The roughness function $\Delta U^+$ is $8.002$ and $8.046$ for the resolved and modelled roughness, respectively.

Figure 3. Comparison of the streamwise mean velocity profiles between IBM and PFA model cases for three smooth-wall elevations. In $(a)$ the homogeneous rough PFA case is presented by the light blue line.

As will be shown later, the streamwise mean velocity exhibits spanwise variations up to the channel centre in some of the heterogeneous rough cases. Nevertheless, these spanwise-averaged velocity profiles can provide some insight into the differences for modelled and resolved roughness reported in table 1. In the case of $h=0$ it can be seen that, in contrast to the homogeneous rough case, differences between IBM and PFA already emerge below $y=\bar {k}$. The PFA model leads to higher average velocities in this region. In this case the roughness strips are exposed to the surrounding flow (compare figure 1). As outlined in § 2, individual roughness elements whose centre is in the rough region can extend with their foothills into the otherwise smooth region (see figure 1$a$) for the present set-up. These roughness elements, which slightly stick out of the rough region, are not modelled in the PFA approach, which is restricted to a force distribution in the rough region (see figure 1$d$). Therefore, the additional drag exerted by the spanwise protruding roughness elements leads to larger global drag for the IBM case (see table 1) and a reduced average streamwise velocity (see figure 3). With increasing $h$ this effect is reduced and eventually flipped since spanwise protruding roughness elements are merged with the surrounding elevated smooth-wall area. In consequence, IBM and PFA results agree much better for $h>0$ (see figure 3$b{,}c$). In the case of $h=1.70 \bar {k}$, $U_b^+$ for the resolved roughness exceeds the one of the corresponding modelled case by 3 %–4 % (depending on the resolution). Translated into a drag coefficient, the PFA model thus produces larger drag than the IBM approach for this case. This difference can be related to an effectively narrower rough-wall region in case of IBM (since roughness elements at the edges are partially merged with the elevated smooth wall) which allows us to generate higher flow rates for the same pressure drop, opposite to the effectively wider rough-wall region for $h=0$. However, the comparison between IBM and PFA for this case is more complex, as can be seen from the streamwise velocity profiles at different spanwise locations discussed in the following.

The data displayed in figure 4 are obtained based on averaging in time, streamwise direction and exploiting the spanwise periodicity. The resulting spanwise coordinate $\bar {z} = z/\delta$ is in the range $0 \leq \bar {z} \leq 1$ with its origin $\bar {z}=0$ placed at the centre of the smooth-wall patch. The centre of the roughness patch is located at $\bar {z}=0.50$ and the transition between rough to smooth occurs at $\bar {z}\approx 0.25$. Therefore, the blue shaded colours correspond to spanwise locations over the smooth surface patch while red shaded colours represent locations over the rough surface part.

Figure 4. Mean streamwise velocity profiles at different spanwise locations $\bar {z}=z/\delta$. The dark blue line shows the centre of the smooth strip at $\bar {z}=0.0$ and the dark red line the centre of the rough strip at $\bar {z}=0.5$. The dashed and dotted black lines represent the streamwise mean velocity of the smooth and homogeneous rough case (IBM).

The velocity profiles at different spanwise locations collapse in the outer flow region for some of the investigated cases only (see figure 4$b{,}f$). The other cases reveal spanwise variations far into the bulk flow, indicating the presence of strong secondary motions, which are addressed in detail in the following section. From figure 4 it is apparent that the influence of the secondary flow on the streamwise mean flow among IBM and PFA differs, especially in the case of $h =1.70 \bar {k}$. A strong spanwise inhomogeneity is present for IBM in the range $0.18 \le y/\delta \le 0.85$ (figure 4$c$) in contrast to a more homogeneous streamwise flow field for PFA (figure 4$f$). The increased velocities for IBM discussed above can be seen to originate from the flow above the rough surface part. In this case the PFA forcing is located below the surrounding smooth-wall strips (see figure 1$f$) while individual IBM roughness peaks reach beyond $y=h$ (see figure 1$c$). At the same time, streamwise velocity can establish in between the IBM roughness elements. We will turn again to the discussion of larger drag for PFA in the case of $h =1.70 \bar {k}$ after the discussion of secondary motions (§ 3.2) and turbulent flow properties (§ 3.3).

To a weaker extent differences between IBM and PFA are also present for $h=0$ and $h\approx \bar {k}$. The high momentum pathways are located above the smooth surface parts for these cases (blue shaded colours) which corresponds to the flow distribution one would expect from a laminar flow (i.e. a flow without secondary motions of Prandtl's second kind) above surfaces with varying friction drag. This spanwise inhomogeneity of $\bar {u}$ is enhanced by the secondary motions typically found above ridge-type roughness, inducing a downwelling motion above the recessed area. This effect appears to be slightly stronger for the modelled roughness. Figure 4($a$,$d$) confirms that the increase of $U_b^+$ in case of PFA for $h=0$ is related to the velocity difference at the transition between smooth- and rough-wall areas ($\bar {z}=0.25$).

In case of increased elevation of the smooth-wall area, the spanwise distribution of the streamwise velocity changes in the sense that high momentum pathways are located above the rough surface parts (red shaded colours). This can be observed in the case of $h=1.70 \bar {k}$ for the resolved roughness and for $h=2.00 - 2.50 \bar {k}$ for the modelled roughness. This different behaviour is related to the reversed rotational direction of the secondary motions addressed in the following section. It is interesting to note that the spanwise inhomogeneity is most pronounced in the outer flow region in these cases, while the classical ridge-type behaviour is dominated by large spanwise variations in the near-wall region.

3.2. Secondary motion

As discussed before, protruding and recessed roughness were previously found to induce secondary motions that exhibit opposite rotational directions if evaluated in a time- and phase-averaged framework. In this context, the protruding roughness strips bear large similarities with ridge-type roughness while recessed roughness induces secondary motions similar to strip-type roughness (Stroh et al. 2020b).

Figure 5 shows the cross-sectional mean flow obtained for resolved and modelled inhomogeneous roughness. The secondary motion is extracted from a phase average (over $L/2$) of the mean flow field obtained through space averaging (along the streamwise direction) and temporal averaging over a total simulation time of at least 50 flow-through time units.

Figure 5. Contours of streamwise mean velocity and the induced secondary motion for resolved roughness cases $(a$–$c)$ and modelled roughness cases $(d$–$k)$.

The white lines with arrows in figure 5 represent the in-plane secondary motion, while the colour code corresponds to the streamwise mean velocity. Isolines of the streamwise mean velocity are plotted as grey lines to indicate the spanwise inhomogeneity of the mean flow. It can clearly be seen that the PFA model is able to predict large-scale secondary motions. These are in very good agreement with the results for the resolved roughness in case of the protruding rough surface, $h=0$ (see figure 5$a$,$d$), revealing an upward motion above the roughness strip that is strong enough to significantly deflect the isolines of streamwise mean velocity. In the near-wall region these deflections are slightly stronger for the modelled roughness. In the case of $h \approx \bar {k}$ (see figure 5$b$,$f$) the modelled roughness induces a secondary motion similar to the one at $h=0$ along with the corresponding bulging of the streamwise flow while the resolved roughness does not reveal any bulging of the streamwise velocity isolines, indicating a weaker secondary motion. The PFA case thus behaves more like a ridge-type roughness. For the recessed roughness ($h=1.70\bar {k}$ in figure 5$c{,}i$) the secondary flow topology is similar between IBM and PFA in the sense that the upward motion is located above the smooth-wall strip for resolved and modelled roughness. However, the secondary motion induces a deflection of the streamwise mean flow in case of the resolved roughness (figure 5$c$), especially at larger wall distance, while this cannot be observed for the modelled roughness (figure 5$i$), indicating weaker secondary currents.

The secondary flow topology in the case of recessed roughness, $h=1.70 \bar {k}$, encompasses only one large-scale vortex pair, while a more complex topology exists for smaller $h$. For $h=0$, in addition to the dominating large-scale vortex pair that reflects in the mean flow bulging, two small vortex pairs can be observed for the resolved and the modelled roughness: one on top of the roughness strip and another one above the smooth-wall area. Both small vortex pairs have an opposite rotational direction compared with the dominating large-scale vortex pair such that a downward motion above the centre of the roughness and an upward motion above the centre of the smooth patch is found in the near-wall region.

For the intermediate rough case, $h \approx \bar {k}$, the topology of the weak secondary motion for the resolved roughness (figure 5$b$) can be understood as a transition between the two other cases. At the edge of the partially protruding roughness strip (figure 5$b$) a first indication of an upward motion can be seen that starts to form the two small-scale vortex pairs still present for $h=0$.

The additional PFA simulations with varying $h$ reveal that this transition of the secondary flow topology can also be captured for the modelled roughness. In this case, a topology similar to (figure 5$b$) is realized with an increased $h$ of $h=1.50 \bar {k}$ (figure 5$h$) which also does not induce any bulging of the streamwise mean flow isolines. The fact that the modelled case requires larger $h$ values in order to match the secondary flow topology of the resolved roughness is also present for the more recessed roughness. The secondary flow and related isoline curvature generated for a resolved roughness with $h=1.70 \bar {k}$ (figure 5$c$) is captured for a modelled roughness with $h= 2.00 - 2.50 \bar {k}$ ($j,k$) which is in agreement with the discussion in § 3.1.

As a measure for the spatial extent and the strength of the secondary motion the spanwise-averaged wall-normal mean velocity magnitude $\langle \tilde {v}\tilde {v} \rangle$ is shown in figure 6. Since spanwise averages integrate a number of different features, especially in cases of complex flow topology, a colour map of the wall-normal spanwise distribution of $\tilde {v}$ is provided in the Appendix in figure 9.

Figure 6. Spanwise-averaged wall-normal dispersive stress profiles for IBM-resolved roughness cases in $(a)$ and PFA-modelled roughness cases in $(b)$.

For the IBM case with protruding roughness ($h=0$), $\langle \tilde {v}\tilde {v} \rangle$ spans a wide wall-normal range with a rather constant value with a local near-wall maximum around $y \approx k_{max}$. This can directly be related to the strong deflection at the smooth–rough interface. While small values for $\langle \tilde {v}\tilde {v} \rangle$ are found for $h \approx \bar {k}$, the recessed roughness case ($h=1.70\bar {k}$) features a single strong peak around $y\approx 0.3\delta$. Interestingly, the streamwise mean velocity isolines at this wall-normal location exhibit a weak bulging only (see figure 5), indicating that the secondary motion acts to homogenize the spanwise distribution of streamwise mean velocity in this case. Comparing $\langle \tilde {v}\tilde {v} \rangle$ for the modelled roughness (figure 6$b$) with the resolved one (figure 6$a$) reveals a number of differences. For $h=0$ the PFA model depicts a stronger peak around $y \approx k_{max}$ and lower values in the bulk of the flow. The stronger peak can be related to a stronger and more localized upward wall-normal velocity at the transition of the smooth to rough wall (see figure 9), while for the IBM model these upward deflections are less intense in the transition region. In terms of $\langle \tilde {v}\tilde {v} \rangle$ the modelled roughness with $h=0.50\bar {k}$ bears the largest resemblance to the resolved $h=0$ case. For $h=1.70\bar {k}$ the modelled roughness induces significantly lower values for $\langle \tilde {v}\tilde {v} \rangle$ as expected from the absence of isoline curvature in figure 5($i$). Actually, the wall nearest isoline in this figure shows a slight outward bulging above the roughness which a stronger secondary motion would annihilate and eventually reverse. For PFA this increase in secondary motion strength can be realized with an increase of $h$ such that not only the secondary flow topology but also the corresponding intensity for $h=1.70\bar {k}$ with IBM is reasonably well matched for $h = 2 - 2.5 \bar {k}$ with PFA.

Within the PFA modelled cases (figure 6$b$) the sole influence of $h$ on $\langle \tilde {v}\tilde {v} \rangle$ (employed here as a simple measure of the secondary flow strength) can nicely be analysed. On the one hand, protruding roughness induces very strong upward motions around $y \approx k_{max}$ which reduce in their intensity for a reduced protrusion of the roughness. On the other hand, recessed roughness also induces strong $\langle \tilde {v}\tilde {v} \rangle$; again with reduced intensity for decreased wall offsets. However, for recessed roughness, the peak of $\langle \tilde {v}\tilde {v} \rangle$ is located at a larger wall-normal distance. This is in agreement with the observation made in respect to figure 4 (§ 3.1), that spanwise inhomogeneity of the streamwise velocity profile is more pronounced in the outer flow region for recessed roughness and more pronounced in the near-wall region for protruding roughness. In the case of recessed roughness, the turbulence generated secondary motions homogenize the surface elevation induced inhomogeneities near the wall. Their influence on the streamwise mean flow is therefore more perceptible at larger wall distances up to which surface-induced inhomogeneities alone do not reach. In contrast, turbulence generated secondary motions over protruding roughness enhance surface-induced inhomogeneities of the flow field such that their combined influence is strongest near the wall.

We note that the observed tendency of the decrease and increase of the wall-normal velocity magnitude with increasing $h$ correlates well with the development of the skin friction coefficient of the PFA model (cf. table 1) and suggests that the appearance of strong secondary motions contributes to drag increase.

3.3. Turbulent flow properties

The flow above rough walls is typically characterized by increased turbulent kinetic energy and shear stress. Forooghi et al. (2018) show that the PFA model captures these features relatively well for homogeneous rough surfaces. Figure 7 shows the spanwise variation of the turbulent kinetic energy $K = (\overline {u'u'} + \overline {v'v'} + \overline {w'w'})/2$ slightly above the maximum roughness elevation at a wall-normal position of $y=0.115\delta$ for the investigated inhomogeneous roughness cases. The increased turbulent kinetic energy above the rough patches and also its spanwise distribution are well captured by the model. Comparing the magnitude of $K$ between resolved and modelled roughness for the same smooth-wall elevation $h$ indicates that the modelled roughness induces similar turbulent kinetic energy to the resolved one. Above the smooth-wall patches $K$ is similar for resolved and modelled roughness in the case of $h = 0.97 \bar {k}$ and $h = 1.70 \bar {k}$ but differs for $h=0$. The turbulent kinetic energy above the smooth area located in between resolved roughness elements is larger in this case. The wall-normal location at which $K$ is extracted corresponds to the one where $\langle \tilde {v}\tilde {v} \rangle$ is largest for the protruding roughness case and where a downwash is present above the smooth wall. The observed difference in $K$ is thus probably related to the different strength of the secondary motion in these two cases.

Figure 7. Spanwise variation of turbulent kinetic energy for IBM-resolved roughness cases in $(a)$ and PFA-modelled roughness cases in $(b)$ extracted at $y=0.115\delta$. Line colours same as in figure 6.

The local peaks of $K$ at the edges of the protruding ($h=0$) or partially protruding ($h=0.5\bar {k}$) roughness are known for ridge-type roughness (Hwang & Lee 2018; Schäfer et al. 2019) and can be related to the upwash that occurs on the corners, while the distribution of $K$ above the recessed roughness is similar to the one above strip-type roughness (Anderson et al. 2015); i.e. with the largest values of $K$ in the centre of the roughness strip. The PFA model captures these qualitative differences well. In particular it can be seen that the distribution of $K$ resembles the one for strip-type roughness for $h=1.25\bar {k}$ and $h=1.50\bar {k}$. Starting from $h=1.70\bar {k}$ local peaks of $K$ at the transition between rough and smooth strips start to emerge again. These peaks are located further towards the smooth region than for smaller $h$, indicating that the elevated smooth surface starts to influence the turbulent flow in a ridge-type manner, despite the fact that $K$ is generally larger above the rough-wall region at the considered wall-normal distance (see cross-sectional distribution of $K$ shown in the Appendix in figure 10).

A change towards ridge-type behaviour for the elevated smooth-wall surface is also visible in the distribution of $\overline {v'w'}$ shown in figure 8 for the $y$–$z$ cross-section of the flows. For protruding edges of roughness strips this component of the Reynolds stress tensor, which represents the wall-normal deflection of spanwise velocity fluctuations, is known to be an important quantity for the formation of secondary motion (Hwang & Lee 2018; Stroh et al. 2020b). The spatial extent of the $\overline {v'w'}$ contours for the protruding roughness cases in figure 8$(a$,$d)$ are similar; however, the modelled case (figure 8$d$) exhibits a localized maximum region at the upper edge of the roughness region, while the region of large $\overline {v'w'}$ is more spread out for the resolved roughness (figure 8$a)$. This can directly be related to the highly localized transition from smooth to rough for the PFA model, which is obviously more gradual for the resolved case if viewed from the streamwise-averaged perspective employed in the present plots. Along the same line the smaller penetration of the $\overline {v'w'}$-contours into the modelled roughness region can also be explained with the spatial homogenization of the model in the rough surface region. The maximum intensity of $\overline {v'w'}$ is larger for figure 8$(a)$, indicating stronger deflection events at the transition from smooth to rough. This is likely to be caused by individual larger roughness elements which are not present in the roughness model.

Figure 8. Contours of Reynolds shear stress component $\overline {v'w'}$ for resolved roughness cases $(a$–$c)$ and modelled roughness cases with low resolution $(d$–$k)$. The coloured dotted lines correspond to $\overline {v'w'}/U_b^2$ values of $( \pm 0.0001$, $\pm 0.0005)$. Isolines of the streamwise mean velocity are shown in grey.

For the intermediate roughness case, $h\approx \bar {k}$, the $\overline {v'w'}$ contours resemble those of the protruding roughness ($h=0$) but with slightly smaller wall-normal extent and lower magnitudes. In this scenario the resolved roughness case (figure 8$b)$ already comprises very small regions of oppositely signed $\overline {v'w'}$ at the smooth-wall edge inside the roughness region, indicating that deflections in the opposite direction are already present. This can only occur for flow located inside individual roughness valleys in which spanwise fluctuations towards the elevated smooth region are deflected upward. This effect is not present in the corresponding homogenized roughness model (figure 8$\,f)$. However, when $h$ is elevated further for the modelled roughness the first indications of oppositely signed $\overline {v'w'}$ at the transition from a smooth to a rough surface occur for $h=1.5 \bar {k}$ (figure 8$h$). This similarity in the qualitative distribution of $\overline {v'w'}$ between figures 8($b$) and ($h$) coincides with a similar secondary flow topology (see figure 5).

The fact that larger values of $h$ are required for the modelled roughness in order to achieve a similar distribution of $\overline {v'w'}$ can also be seen for the recessed roughness. The resolved roughness $h=1.70 \bar {k}$ (figure 8$c$) bears large similarities to the modelled roughness $h=2.00 \bar {k}$ (figure 8$j$). Again, this case exhibits a highly similar secondary motion which also induces similar deflections to the streamwise mean velocity profile. In general, an increase in $h$ induces a sign reversal of $\overline {v'w'}$, indicating that the smooth-wall area starts to act as the protruding surface part as already noted in respect to the spatial distribution of $K$. This is most pronounced for $h=2.50 \bar {k}$ (figure 8$k$). The sign reversal does not occur at identical $h$-values for resolved and modelled roughness but requires a larger smooth-wall elevation for the modelled roughness. This shift can be related to the homogenizing nature of the PFA model which does not capture the turbulence present in individual roughness valleys and its interaction with the elevated smooth wall. Overall, the present results confirm that $\overline {v'w'}$ has a strong impact on the formation of secondary flows and their respective strength.

The limiting case of strip-type roughness, in which a smooth wall with different shear stress is considered, cannot be realized with the present scale separation where deflections at the roughness edges appear to be present in all cases. Considering the $\overline {v'w'}$ distribution, the case of $h=\bar {k}$ is still clearly dominated by the ridge-type behaviour of the rough surface part. For the PFA cases the minimum influence of $\overline {v'w'}$ appears to exist around $h=1.5\bar {k}$ for the present roughness strips which is also the case with the lowest wall-normal mean velocities (cf. figure 6). Assuming $\overline {v'w'}$ to be an indicator for protrusion triggered secondary motion (Hwang & Lee 2018) suggests that $h \approx 1.5\bar {k}$ most closely generates a strip-type behaviour. For the present roughness configuration this case does not induce secondary motions that significantly alter the streamwise mean flow.

4. Discussion and conclusions

The PFA roughness model developed for the modelling of homogeneous rough surfaces (Busse & Sandham 2012; Forooghi et al. 2018) is used to predict the turbulent flow over spanwise heterogeneous roughness for different relative roughness elevations, i.e. protruding vs recessed streamwise roughness strips. The roughness elevation is varied indirectly through a stepwise elevation of the intermediate smooth-wall strip height $h$. For three cases of $h$ a direct comparison with roughness-resolving DNS based on IBM is available from Stroh et al. (2020b).

While the roughness model is applied in a clearly defined spanwise region of width $W$, the roughness-resolving IBM based DNS is not set up with a sharp cut at the edges of the rough region since this would have altered the roughness statistics of the rough strip. Instead, all roughness elements, whose centre position is located outside the rough region, are removed to create smooth strips. In consequence, foothills of roughness elements located close to the roughness edges can reach either onto the neighbouring smooth strips in the case of protruding rough strips or merge with the elevated neighbouring smooth strips in the case of recessed roughness. Therefore, the effective width of the rough strip is slightly larger for the protruding roughness ($h=0$) and slightly smaller for the recessed roughness ($h=1.70\bar {k}$). This difference directly reflects in the global flow properties in the sense that relatively larger drag is generated for resolved protruding roughness (due to effectively wider roughness strips) and relatively smaller drag is present for resolved recessed roughness (due to effectively narrower roughness strips).

The comparison of the global effective friction coefficient achieved with the resolved and modelled roughness approaches reflects the discussion above. Compared with the IBM based DNS, the PFA roughness model predicts a smaller friction coefficient for $h=0$ and a larger one for $h=1.70\bar {k}$ while good agreement is obtained for $h \approx \bar {k}$. In general, the relative drag increase of the heterogeneous roughness compared with the area-averaged mean value for smooth and rough contributions is well captured by the PFA model. Its simpler structure compared with IBM resolved roughness not only allows us to reduce the computational effort to the one of smooth-wall turbulent flow DNS but also to systematically investigate the effect of single parameters of heterogeneous rough surfaces, such as e.g. the relative roughness elevation in the present study, in a more straightforward manner.

In this context, the present results show that the PFA roughness model is able to capture all salient features of heterogeneous rough surfaces. In particular, the PFA model captures the transition of secondary motion from ridge- to strip-type behaviour through variation of $h$ and the related reversed rotational direction (Stroh et al. 2020b) which cannot be resolved with other numerical roughness modelling approaches that rely on the prescription of an effective wall shear stress (Anderson et al. 2015; Chung et al. 2018).

The relative drag increase of spanwise inhomogeneous surfaces in turbulent channel flows is to a first approximation qualitatively similar to the one obtained in corresponding laminar flow (Daschiel et al. 2012); i.e. relative drag increase is obtained for spanwise wavelengths of the order of the channel height (as in the present case) while relative drag decrease can be obtained with much larger wavelengths. The present PFA based results indicate the presence of a relatively constant drag coefficient in the range of $0.5\bar {k}< h<1.5\bar {k}$ which coincides with weak secondary motions. An increasing strength of secondary motions induces an increase of the drag coefficient. In the case of protruding roughness ($h=0$) this drag increase occurs because the secondary motion enhances the inhomogeneity of the streamwise velocity distribution. In case of recessed roughness, when the smooth-wall region starts to emerge as a protruding region, the secondary motion is such that it enhances low momentum pathways over the smooth regions and high momentum pathways in the rough regions. While this phenomenon occurs first in the bulk of the flow ($h=1.70\bar {k}$) a further increase of $h$ also induces clearly visible high momentum pathways in the near-wall region of the rough surface part ($h=2.50\bar {k}$). Therefore, the relative influence of the rough region on the total drag is increased. Overall, the presence of secondary motions (which are triggered through protruding rough or smooth strips in the present study) generally leads to larger drag increase; however, the underlying physical mechanisms differ for protruding and recessed roughness strips.

A direct comparison of the secondary flow topology between modelled and resolved roughness reveals that a shift towards larger values of $h$ is required for the PFA model to reproduce the IBM results. This suggests that the PFA roughness model induces a larger wall offset when applied in heterogeneous instead of homogeneous rough-wall flow conditions. While the PFA forcing (2.2) is applied in streamwise and spanwise directions only, the local wall-normal mean velocity shown in figure 9 shows that the PFA model nevertheless prevents the occurrence of downwash into the roughness region (negative values of $\tilde {v}$) which is visible for the resolved roughness. The secondary motions present for heterogeneous roughness generally induce a net downward motion above the roughness centres which penetrates further into the resolved rough region. In consequence, the fluid in resolved roughness valleys can interact with neighbouring smooth regions. This is particularly important in the case of recessed roughness strips where wall-normal deflections of spanwise velocity fluctuations can occur at the edges of the protruding smooth-wall region. The reduced wall-normal penetration in case of the PFA model thus yields smaller values of $\overline {v'w'}$ which is known to be the primary driving force for the secondary motions above ridge-type roughness (Hwang & Lee 2018). The reduced secondary motion strength predicted with the PFA roughness model compared with IBM for $h=1.70\bar {k}$ is thus caused by the reduced intensity of $\overline {v'w'}$. The comparison of the global drag for these two cases (see table 1) is governed by two factors: PFA induces lower drag due to weaker secondary motions, but higher drag due to the effectively narrower roughness strip width for IBM. The latter phenomenon dominates in the present set-up, which indicates a possible advantage of using a simplified roughness model instead of IBM if the sensitivity of global flow properties towards individual parameters is investigated. For such investigations the PFA roughness model can be used without any particular fine tuning.

Figure 9. Contours of local wall-normal mean velocity $\tilde {v}$ for resolved roughness cases $(a$–$c)$ and modelled roughness cases $(d$–$k)$. Isolines of the streamwise mean velocity are shown in grey.

While the secondary flow generated above protruding rough strips ($h=0$) resembles that of protruding smooth ridges with smooth valleys (Hwang & Lee 2018; Medjnoun et al. 2018; Stroh et al. 2020a) or ridges and valleys with identical roughness (Vanderwel & Ganapathisubramani 2015; Vanderwel et al. 2019), the present results indicate an interesting difference for a protruding smooth ridge with rough valleys. The latter encompasses one large-scale vortex pair only, while the other ridge-type cases typically contain three vortex pairs, as discussed in relation to figure 5($a$,$d$). This difference probably originates from the fact that the turbulent kinetic energy is consistently higher above the rough surface region than above the smooth one. In agreement with the suggestion of Hinze (1967, 1973), the secondary motion in the near-wall region is always directed from the region of high $K$ towards the region of lower $K$ (located above the smooth surface region). Deflections ($\overline {v'w'}$) on a protruding rough ridge counteract the secondary motion induced by the spanwise gradient of $K$, leading to the formation of multiple vortex pairs. However, the opposite occurs over surfaces with rough valleys. In this case deflections on the edges of the protruding smooth surface region enhance the secondary motion induced by the gradient of $K$ such that the secondary motions do not indicate a difference between protruding ridge-type surface structures and non-protruding strip-type ones.

We note that a pure strip-type behaviour can probably not be obtained with the scale separation ($\bar {k}/\delta$) of the present DNS, since $\overline {v'w'}$ generated at any protruding surface part appears to have a strong influence on the secondary flow formation. However, the corresponding influence is weak for $h \approx 1.5 \bar {k}$, indicating a small difference in the effective virtual wall location between smooth and rough surface parts for the present configuration. Since protruding smooth surface parts will enhance the ridge-type secondary flow, it is advisable to work with such a configuration if the effect of strip-type roughness is to be investigated numerically (with resolved roughness) or experimentally.

It is likely that a further increase of $h$ beyond the parameter space of the present study will reduce the relevance of high $K$ above the rough surface such that the peak of $K$ in the near-wall region of the elevated smooth surface will eventually dominate along a wall-parallel line. In this case, a recovery of the classical ridge-type secondary flow topology is expected.

We conclude that the PFA roughness model constitutes an attractive alternative for roughness-resolving DNS when investigating not only the effect of homogeneously but also of heterogeneously rough surfaces in turbulent channel flows. The PFA roughness model allows one to systematically study the effect of individual roughness parameters which can often be strongly interlinked for the resolved case. We are of the opinion that such parameter studies do not necessarily require a fine tuning of the PFA model to a particular rough surface (which would require additional roughness-resolving DNS) but can be used in an a priori manner, at least for turbulent channel and equilibrium boundary layer flows. However, for significantly different flow conditions (especially including flow impingement on rough surfaces) the model should not be applied without further checks.

Finally, we note that roughness models are also required in LES, which have been used extensively to study roughness-induced secondary motions at large Reynolds numbers (Willingham et al. 2014; Yang & Anderson 2018; Forooghi, Yang & Abkar 2020). In these simulations the small-scale height fluctuations of the rough surface typically lie below the first grid point, which necessitates the use of so-called wall models. At the present stage, these models have not been shown to predict the different flow behaviour for protruding and recessed roughness strips. Future work in this field should be dedicated to the question of how to improve these model approaches to also enable the prediction of the more complex secondary flow topology present over protruding roughness strips.