1. Introduction
Scientists have, for years, been documenting the relationship between surface roughness and its hydraulic resistance, the former pertaining to geometry of the surface while the latter is dictated by the dynamics of flow over the roughness (see reviews by Jiménez Reference Jiménez2004, and Flack & Schultz Reference Flack and Schultz2010). The cataloguing process can be a never-ending exercise because the characteristics of each rough surface are unique. In order to make predictions on full-scale operations, it is necessary to establish the equivalent sand-grain roughness
$k_{s}$
of a given surface, which relates the drag increment of the given surface to an equivalent surface composed of monodisperse sand grains. Once
$k_{s}$
has been determined, it is possible to predict the drag penalty at application Reynolds numbers using either the Moody chart (Moody Reference Moody1944), for pipe and channel flows, or developments of this for zero-pressure-gradient boundary layers (Prandtl & Schlichting Reference Prandtl and Schlichting1955; Granville Reference Granville1958). Generally speaking, the approach to date has been to first identify a particular rough surface of scientific or engineering interest, and then to characterise its hydraulic resistance through well-controlled laboratory experiments. By exposing the rough surface to an increasing range of flow speeds, until such point beyond which the resistance coefficient
$C_{f}$
becomes constant, referred to as the ‘fully rough’ asymptote, it is possible to ascertain
$k_{s}$
. A convenient dimensionless group here is the equivalent roughness Reynolds number
$k_{s}^{+}\equiv k_{s}U_{{\it\tau}}/{\it\nu}$
, where the friction velocity
$U_{{\it\tau}}\equiv \sqrt{{\it\tau}_{0}/{\it\rho}}\equiv U_{\infty }\sqrt{C_{f}/2}$
;
${\it\tau}_{0}$
is the wall drag force per plan area;
${\it\rho}$
is the density;
$U_{\infty }$
is the free-stream velocity or centreline velocity for internal geometries; and
${\it\nu}$
is the kinematic viscosity. For accurate predictions, it is not sufficient to merely establish
$k_{s}$
; rather the increment of
$C_{f}$
caused by surface roughness must be mapped with respect to
$k_{s}^{+}$
at all conditions from the dynamically smooth up to the fully rough regimes, which generally covers the approximate range
$5\lesssim k_{s}^{+}\lesssim 100$
. Over the last century or so, this painstaking and time-consuming procedure has been repeated for many surfaces of interest, and as a result there is now a rich database of roughness amassed in the published literature over time.
One problem here is that there is no widely applicable function that relates
$k_{s}$
, a dynamic parameter, to readily observable or measurable geometric properties of the surface, such as the root-mean-square (r.m.s.) roughness height
$k_{rms}$
or average roughness height
$k_{a}$
. Certain classes of surfaces, say sand-grain roughness, may exhibit an approximate proportionality between
$k_{s}$
and some physical surface length scale, but as a general rule this proportionality will not hold across different surface topologies. There are numerous attempts in the literature to formulate more complicated functions, often involving some measure of mean surface height such as
$k_{rms}$
or
$k_{a}$
and other properties relating to the shape or arrangement of roughness elements such as skewness, effective slope, solidity and so on (see review by Flack & Schultz Reference Flack and Schultz2010). Though many of these parameters have some success in describing the particular class of surfaces for which they were formulated (e.g. painted or sanded surfaces), none are widely applicable across the almost limitless range of surface topologies that are encountered in engineering and meteorological applications.
The present method to directly evaluate
$k_{s}$
represents a paradigm shift from such parametrisations that are based on a handful of geometrical factors. Recognising that all roughness geometries are unique, and that a one-size-fits-all formulaic solution is proving elusive, we have sought an approach to minimise the expense involved in experimentally determining
$k_{s}$
. Specifically, the present method relates raw profilometry data, described by many degrees of freedom, directly to
$k_{s}$
. In contrast, previous indirect approaches first reduce raw profilometry data to a limited subset of geometrical factors such as
$k_{rms}$
and skewness before relating these factors to
$k_{s}$
. The inevitable loss of information in reducing raw profilometry data to geometrical factors is avoided by the present method and, in this sense, the present method represents a direct evaluation of the roughness function.
The approach relies on the direct computation of hydraulic resistance by direct numerical simulation (DNS), which we presently show can be made substantially cheaper, and therefore faster, than previously thought. We circumvent the otherwise prohibitive cost of conventional DNS by employing minimal-span channels (Jiménez & Moin Reference Jiménez and Moin1991; Hwang Reference Hwang2013), the rationale of which is discussed in the following. Normally, a straightforward and direct computation of roughness drag using DNS employing full-span channels is extremely expensive as it entails simultaneously capturing both the bulk flow, which scales with the half-channel height,
$h$
, and the near-wall flow around the roughness elements, which scales with a characteristic roughness height,
$k$
. Given a rough surface of fixed blockage ratio
$k_{s}/h\lesssim 1/40$
(Jiménez Reference Jiménez2004; Flack, Schultz & Shapiro Reference Flack, Schultz and Shapiro2005), a complete characterisation of hydraulic resistance requires parametric simulations that sweep through the equivalent roughness Reynolds numbers,
$k_{s}^{+}\approx 5$
–100, corresponding to the hydraulically smooth and the fully rough regimes, respectively. For the blockage ratio,
$h/k_{s}=40$
, this means performing parametric simulations at the friction Reynolds numbers,
$\mathit{Re}_{{\it\tau}}\equiv hU_{{\it\tau}}/{\it\nu}=(h/k_{s})(k_{s}U_{{\it\tau}}/{\it\nu})\approx 200$
–4000, which are currently unfeasible. Recall that the cost of DNS, counting the number of spatial and temporal degrees of freedom, scales unfavourably as
$\mathit{Re}_{{\it\tau}}^{3}$
(Pope Reference Pope2000, § 9.1.2). For grids conforming to the surface of the roughness elements, this cost is further exacerbated by the need for increased mesh density, and reduced time steps.
However, the extreme cost associated with conventional DNS employing full-span channels seems unnecessary. The quantity of interest from an engineering point of view is the retardation in the mean flow over the roughness relative to the smooth-wall flow. This relative flow retardation or downward velocity shift,
${\rm\Delta}U$
, occurs mostly in the vicinity of the roughness layer, but holds constant above a few roughness heights, well into the log layer (if it exists) and the wake region, a behaviour described by Townsend’s outer-layer similarity hypothesis (Townsend Reference Townsend1976). This suggests that a simulation of only the near-wall region and its interaction with the roughness geometry is required in order to extract
${\rm\Delta}U^{+}\equiv {\rm\Delta}U/U_{{\it\tau}}$
, which is known as the (Hama) roughness function. This distillation of the problem is consistent with the observation that
${\rm\Delta}U^{+}$
does not depend on the bulk flow but only on
$k_{s}^{+}$
and other details of the roughness geometry.
A framework for simulating only the near-wall dynamics is the minimal channel. The concept is first described by Jiménez & Moin (Reference Jiménez and Moin1991) and is currently receiving renewed attention in various contexts of understanding wall-bounded turbulence (Flores & Jiménez Reference Flores and Jiménez2010; Hwang Reference Hwang2013; Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014). Here we exploit this framework for measuring
${\rm\Delta}U^{+}$
by fully resolving the near-wall Navier–Stokes dynamics and its interaction with the roughness geometry. The prohibitive cost of conventional DNS is alleviated by use of these minimal-span channels, which are designed to preclude the bulk flow that scales with
$h$
. Without the bulk flow, the cost of DNS with roughness now only scales as
$k_{s}^{+3}$
, which is quite feasible for the engineering task at
$k_{s}^{+}\approx 5$
–100. In principle, the computational cost is potentially
$(h/k_{s})^{3}$
times less than that of a conventional DNS in a full-span channel.
In the remainder, we demonstrate the efficacy of this approach and develop guidelines for its use. Beginning first with the parametric forcing model of Busse & Sandham (Reference Busse and Sandham2012), we carefully confirm that the minimal channels return the same estimate for the roughness function as that given by full-span simulations (§ 3.1). An important subtlety here lies in the choice of the minimal spanwise unit, which we demonstrate is related not only to the usual span of the near-wall cycle (approximately 100 viscous wall units), but also to some physical roughness height
$k$
and, presumably, the roughness texture under investigation. Having established the accuracy with the simple forcing model, we then test this procedure with the more realistic test case of a three-dimensional roughness, explicitly described by grids conforming to the rough surface, and compare results with full-domain simulations (§ 3.2).
Computational domain of (a) the full channel and (b) the minimal channel. The magnified view of the channel surface illustrates the rough wall represented by (c) the roughness-forcing region,
$0<z<k$
, or (d) the explicit sinusoidal roughness with semi-amplitude
$k$
and mean at
$z=0$
. Both top and bottom channel walls are rough.
Simulation cases listed with nominal
$\mathit{Re}_{{\it\tau}}$
. Each of these 17 cases is run using both smooth and rough walls, that is, a total of 34 separate simulations are run. In Set 1 (modelled roughness)
${\it\alpha}k=1/40$
and in Set 2 (explicit sinusoidal roughness)
${\it\lambda}_{x}/k={\it\lambda}_{y}/k\approx 7.1$
. The grid spacing is uniform in the streamwise and spanwise direction, while the grid spacing in the wall-normal direction is stretched, with the finer wall grid spacing
${\rm\Delta}z_{w}$
and the coarser centreline grid spacing
${\rm\Delta}z_{c}$
.
2. Direct numerical simulations
We present two sets of direct numerical simulations (table 1). In Set 1, roughness is modelled by body forces (figure 1
c) in a fourth-order staggered-grid code as described by Morinishi et al. (Reference Morinishi, Lund, Vasilyev and Moin1998). The code is written by the first author and has been used in other studies (e.g. Chung, Monty & Ooi Reference Chung, Monty and Ooi2014). In Set 2, roughness is explicitly represented by a roughness-conforming grid (figure 1
d) in the code CDP described by Mahesh, Constantinescu & Moin (Reference Mahesh, Constantinescu and Moin2004). The purpose of the Set 1 simulations with modelled roughness is to test the present minimal-channel method for a hypothetical roughness in the absence of any imposed wall-parallel length scales. Such a configuration, which we shall call homogeneous roughness, is ideal for investigating how the characteristic roughness height (as opposed to wall-parallel length scales) sets the minimum allowable span
$L_{y}$
of the channel such that
${\rm\Delta}U^{+}$
is still accurately captured. The purpose of the Set 2 simulations with grids conforming to the rough surface is to demonstrate that the present minimal-channel method works as expected for a specific roughness geometry, namely the sinusoidal ‘egg-carton’ roughness that has previously been characterised (in a pipe) by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015).
In both sets of simulations, we solve the following Navier–Stokes equations of motion:
where
$u_{i}$
is the velocity,
$t$
is time,
$x_{j}$
is the spatial coordinate,
$p$
is the pressure, and
$f(t)$
is the spatially uniform, time-varying, mean pressure gradient that drives the flow at constant mass flux. The streamwise, spanwise and wall-normal directions are referred to as either
$x_{1}$
,
$x_{2}$
and
$x_{3}$
or
$x$
,
$y$
and
$z$
respectively. Periodic boundary conditions are imposed in the streamwise and spanwise directions with the respective domain sizes,
$L_{x}$
and
$L_{y}$
(figure 1
a,b). In the reference full-span simulations,
$L_{y}$
measures
${\rm\pi}h$
for Set 1 and
$2{\rm\pi}h$
for Set 2 (at several selected
$\mathit{Re}_{{\it\tau}}$
), while in the minimal-span simulations,
$L_{y}$
is constrained by the consideration of three issues:
-
(i) the span of the near-wall cycle (assumed to be approximately
$100\,{\it\nu}/U_{{\it\tau}}$
), -
(ii) the span of the roughness elements,
${\it\lambda}_{y}$
(absent for Set 1 simulations), and -
(iii) the physical height of the roughness elements,
$k$
.
Regarding point (i), previous studies (Jiménez & Moin Reference Jiménez and Moin1991; Hwang Reference Hwang2013) have shown that
$L_{y}^{+}>100$
or else the near-wall cycle, which is the signature of wall turbulence, cannot be properly captured. Regarding point (ii), the span of the channel must be wide enough to accommodate the widest scale of the roughness elements in order to properly represent the nature of the roughness in question. For the Set 1 simulations, which have a homogeneous forcing, this point is not an issue. Regarding point (iii), it turns out that the height of the roughness elements also plays a role in setting the minimum
$L_{y}$
, an issue that will be discussed in § 3.1. The chosen resolutions (table 1) are comparable to other channel-flow direct numerical simulations (e.g. Moser, Kim & Mansour Reference Moser, Kim and Mansour1999; Bernardini, Pirozzoli & Orlandi Reference Bernardini, Pirozzoli and Orlandi2014) and the streamwise domain lengths are long enough to accommodate several near-wall streaks, which are approximately
$1000\,{\it\nu}/U_{{\it\tau}}$
. Although the present study focuses on minimising the spanwise dimension of the channel, which is the critical dimension in setting the dynamics of the log layer (Flores & Jiménez Reference Flores and Jiménez2010), one presumes that the streamwise dimension will also be subject to constraints. For a smooth wall, Chin et al. (Reference Chin, Ooi, Marusic and Blackburn2010) show that the mean velocity profile is slightly overestimated when
$L_{x}^{+}\lesssim 1000$
(at
$\mathit{Re}_{{\it\tau}}\approx 180$
for pipe flow), suggesting that the near-wall streaks need to be properly captured. It is also clear that
$L_{x}$
should be greater than
${\it\lambda}_{x}$
, where
${\it\lambda}_{x}$
is the largest characteristic streamwise scale of the roughness elements. The characteristic roughness sizes
$k$
are selected so that they are a fixed and small fraction of the half-channel height
$h$
. The roughness Reynolds number
$k^{+}$
is increased by increasing the overall friction Reynolds number
$\mathit{Re}_{{\it\tau}}$
, which is the typical way roughness is encountered in practice (experimentally, and in practice, the roughness Reynolds number is usually increased by increasing the flow speed).
In the Set 1 simulations (modelled roughness), the equations of motion, (2.1), are numerically solved between two no-slip, impermeable walls at
$z=0$
,
$2h$
and the effect of roughness is represented by the parametric-forcing model of Busse & Sandham (Reference Busse and Sandham2012), whereby a forcing term is added to the right-hand side of (2.1
a) of the form
In principle, the effect of any roughness geometry, which includes both pressure and viscous drag, can be formally written as a body force on the right-hand side of the Navier–Stokes equation, but for the present purposes, the form (2.2) is meant to represent a hypothetical homogeneous roughness. The roughness forcing is only active in the streamwise direction, opposes the mean flow and always dissipates kinetic energy. The roughness factor,
${\it\alpha}$
, is thought to scale with the roughness density, that is, the frontal area per unit volume (Nikora et al.
Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007; Busse & Sandham Reference Busse and Sandham2012), measured in inverse-length (area per unit volume) units. Here,
$k=h/40$
and
${\it\alpha}=1/(40k)$
for all simulations in Set 1 (refer to table 1 for details). A cosine mapping is used to stretch the grid in the wall-normal direction.
Grid for the Set 2
$\mathit{Re}_{{\it\tau}}\approx 180$
minimal-span simulation.
In the Set 2 simulations, the sinusoidal roughness previously studied by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015) is explicitly represented by a roughness-conforming grid composed of hexahedral cells (figure 2). The grid is stretched in the wall-normal direction, with the stretching ratio held below 1.05. Grid information for all cases in Set 2 is detailed in table 1. The roughness surface described by the sinusoidal function,
locates the no-slip, impermeable wall relative to the reference plane at
$z=0$
. Under this definition, we consider the characteristic roughness height
$k$
to be the semi-amplitude of this sinusoid. For the present sinusoidal roughness, the peak-to-valley roughness height
$k_{t}~(=2k)$
, the r.m.s. roughness height
$k_{rms}~(=k/2)$
or other physical roughness heights could instead be used. Ultimately, the present method seeks to determine
$k_{s}/k$
directly and the specific choice of the physical roughness height is therefore irrelevant, as all one requires in practice is the ratio between the equivalent sand-grain roughness and some readily measurable geometric property of the rough surface. For the present sinusoidal roughness,
$k_{s}/k=2k_{s}/k_{rms}=k_{s}/k_{t}/2$
. The roughness surface at the top wall is the mirror image (across the centreline) of the roughness surface at the bottom wall. Here,
$k=h/18$
and
${\it\lambda}_{x}={\it\lambda}_{y}\approx 7.1k$
for all simulations in Set 2 (refer to table 1 for details).
Mean velocity profiles of simulated turbulent channel flow at (a)
$\mathit{Re}_{{\it\tau}}\approx 590$
and (b)
$\mathit{Re}_{{\it\tau}}\approx 4000$
: solid, smooth; dashed, modelled roughness; grey, minimal span; black, full span. The vertical dashed line marks the top of the roughness-forcing region. The inset shows that the velocity shift stays constant for both minimal- and full-span channels above the roughness-forcing region. The various grey profiles in (b) correspond to the various minimal spans,
$L_{y}^{+}\approx 116$
, 232 and 463. The inset in (b) shows that
$L_{y}^{+}\approx 116$
is too narrow to accurately capture the roughness function for roughness heights above
$k^{+}=100$
. The full-span
$\mathit{Re}_{{\it\tau}}\approx 4000$
data in (b) are from Bernardini et al. (Reference Bernardini, Pirozzoli and Orlandi2014).
3. Results and discussion
3.1. Modelled roughness
In figure 3(a), we show the mean velocity profile from four simulations with modelled roughness (Set 1 simulations) at
$\mathit{Re}_{{\it\tau}}\approx 590$
, corresponding to the smooth full-span, the rough full-span, the smooth minimal-span and the rough minimal-span simulations. Consistent with the findings of Hwang (Reference Hwang2013), figure 3(a) shows that the minimal span of
$L_{y}^{+}\approx 116$
is sufficient for capturing the mean smooth-wall velocity profile given by the full-span simulations up to the height
$z_{c}^{+}\approx 40$
. We will use the notation
$z_{c}^{+}$
to refer to the point above which the minimal profile departs from the full-span profile. The present simulations show that this is also the case for the profile above modelled rough walls (given by the dashed lines in figure 3
a). Above the roughness height of
$k^{+}=\mathit{Re}_{{\it\tau}}/(h/k)\approx 590/40\approx 15$
, the inset of figure 3(a) shows that
${\rm\Delta}U^{+}$
for the full-span channels is reproduced by the minimal-span channels, even though
$U^{+}$
itself above
$z_{c}^{+}\approx 40$
for the minimal-span channels fails to capture
$U^{+}$
of the full-span channels. We will use the notation
${\rm\Delta}U^{+}$
to refer to the nominally constant value of the downward velocity shift at matched
$\mathit{Re}_{{\it\tau}}$
,
$U_{s}^{+}-U_{r}^{+}$
, evaluated above the roughness height. The similarity of the mean relative velocity is, of course, a manifestation of Townsend’s outer-layer similarity: the downwards shift in velocity is purely determined by the near-wall flow over the roughness, which sets the wall drag. This shift will remain the same irrespective of the state of the outer profile, which is vastly different for minimal- or full-span channels. These results demonstrate the efficacy of the present method. It is important to emphasise that the aforementioned method to evaluate
${\rm\Delta}U^{+}$
requires matched
$\mathit{Re}_{{\it\tau}}$
, which is simple to achieve computationally since the outer length scale and the driving pressure gradient (and therefore
$U_{{\it\tau}}$
) can be fixed between smooth- and rough-wall simulations. Experimentally, where
$\mathit{Re}_{{\it\tau}}$
between smooth and rough walls are seldom well-matched,
${\rm\Delta}U^{+}$
is typically evaluated from a shift in the log region (provided it exists).
When the top of the roughness-forcing region is higher than
$z_{c}^{+}$
, a wider minimal channel is required. A systematic study with various channel spans at
$\mathit{Re}_{{\it\tau}}\approx 4000$
, as listed in table 1, investigates this effect. These results are shown in figure 3(b). For the smooth walls, the data suggest that the minimal channel is able to reproduce the full channel profiles up to
$z_{c}^{+}\approx 0.4\,L_{y}^{+}$
. That is, the minimal profile from
$L_{y}^{+}\approx 116$
reproduces the full profiles up to
$z_{c}^{+}\approx 0.4\,(116)\approx 46$
, the minimal profile from
$L_{y}^{+}\approx 232$
reproduces the full profiles up to
$z_{c}^{+}\approx 0.4\,(232)\approx 93$
, the minimal profile from
$L_{y}^{+}\approx 463$
reproduces the full profiles up to
$z_{c}^{+}\approx 0.4\,(463)\approx 185$
, and the pattern presumably goes on as long as
$z_{c}^{+}\approx 0.4\,L_{y}^{+}$
remains in the log layer
$z^{+}<0.15\,\mathit{Re}_{{\it\tau}}$
, where the size of eddies is thought to scale with distance from the wall. This behaviour is first reported by Flores & Jiménez (Reference Flores and Jiménez2010), who suggested
$z_{c}^{+}\approx 0.3\,L_{y}^{+}$
. We refer to the region
$z<z_{c}$
as being ‘unconfined’ since the turbulence in this region is faithfully captured, showing no signs of being constrained by the minimal span. A difficulty arises when the top of the roughness rises above the unconfined region (when
$k>z_{c}$
). For the narrowest channel, we observe that the roughness function is over-predicted, giving
${\rm\Delta}U^{+}\approx 11$
instead of the
${\rm\Delta}U^{+}\approx 9$
obtained for the two wider (yet still minimal) channels. This sensitivity leads to an additional constraint for the minimum channel span, both of which must be satisfied:
Physically, this means that the channel must be wide enough to fully immerse the roughness elements (
$z_{c}>k$
) in natural, unconfined wall turbulence (which only occurs when
$L_{y}^{+}>100$
). Although the constraint that
$L_{y}>k/0.4$
is developed from the present modelled homogeneous roughness, where
$k$
is the modelled roughness height, we expect that, in general,
$L_{y}>O(k)/0.4$
for any physical roughness height
$k$
such as
$k_{rms}$
or
$k_{a}$
. In practice, a convergence study with wider (still minimal) channels can be performed (see figure 3
b) to manage this effect. A conservative approach is to take
$L_{y}>k_{t}/0.4$
, where
$k_{t}$
is the maximum peak-to-valley roughness height. The third constraint, requiring that
$L_{y}>{\it\lambda}_{y}$
(where
${\it\lambda}_{y}$
is some characteristic spanwise length scale of the roughness texture, or repeating unit of the roughness) is not active for this case of modelled roughness, where the homogeneous forcing suggests
${\it\lambda}_{y}\rightarrow 0$
.
A sweep in roughness Reynolds numbers, from
$k^{+}\approx 5$
to 100, corresponding to
$\mathit{Re}_{{\it\tau}}\approx 180$
to 4000 (table 1) is needed in order to fully characterise the roughness transition from the hydraulically smooth to the fully rough regime. This in turn enables determination of the equivalent sand-grain roughness
$k_{s}/k=C$
, where
$C$
is a constant that best scales the highest roughness data onto the fully rough asymptote. Through the use of minimal-span channels, even the simulations at
$\mathit{Re}_{{\it\tau}}\approx 4000$
are feasible. The Hama roughness function for this sweep is presented in figure 4. The data from the full-span simulations at
$\mathit{Re}_{{\it\tau}}\approx 590$
and 950 are shown by the filled circles. It is noted that the roughness function profile obtained from the minimal-span simulations matches very well with the full-span simulations. Fitting to the fully rough regime, shown by the solid curve on figure 4, the present data show that
$k_{s}/k\approx 1.6$
for
${\it\alpha}=1/(40k)$
and so we have characterised the present modelled roughness from the transitionally rough regime,
$k_{s}^{+}\approx 7$
, to the fully rough regime,
$k_{s}^{+}\approx 160$
. The obtained value for
$k_{s}$
could now, in theory, be used to predict full-scale performance under this roughness condition (Granville Reference Granville1958; Schultz Reference Schultz2007). Comparing the minimal data with the uniform sand-grain roughness data of Nikuradse (Reference Nikuradse1933), shown in light grey markers, reveals that the parametric-forcing model with the box-profile shape function closely mimics the roughness transition of uniform-sand-grain roughness (as previously noted by Busse & Sandham Reference Busse and Sandham2012).
The Hama roughness function at fixed
$h/k=40$
for the modelled roughness, where
$k_{s}/k\approx 1.6$
, simulated in channels. The superimposed uniform sand-grain roughness data of Nikuradse (Reference Nikuradse1933) are extracted from Jiménez (Reference Jiménez2004).
Mean velocity profiles of turbulent channel flow: (a)
$\mathit{Re}_{{\it\tau}}\approx 180$
, smooth and explicit roughness (
$h/k=18$
) in minimal- and full-span channels; (b)
$\mathit{Re}_{{\it\tau}}\approx 180$
, 360, 540, 1080, smooth and sinusoidal roughness (
$h/k=18$
) in minimal-span channels. The vertical dashed line marks the top of the roughness element. The inset shows that the velocity shift stays the same for both minimal- and full-span channels above the roughness elements.
The Hama roughness function at fixed
$h/k=18$
for the sinusoidal roughness, where
$k_{s}/k\approx 4.1$
, simulated in channels and pipes. The pipe data are from Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015) and have the maximum blockage ratio,
$R/k=6.75$
. The superimposed uniform sand-grain roughness data of Nikuradse (Reference Nikuradse1933) are extracted from Jiménez (Reference Jiménez2004).
3.2. Sinusoidal roughness
Figure 5(a) is a plot of the mean velocity profile from both the full and the minimal
$L_{y}^{+}\approx 141$
channel with explicitly represented sinusoidal roughness elements (Set 2 simulations) at
$\mathit{Re}_{{\it\tau}}\approx 180$
. Consistent with the simulations with modelled roughness, the mean profiles from the minimal channel collapse onto the mean profiles from the full channel below about
$z_{c}^{+}\approx 0.4\,L_{y}^{+}\approx 56$
. For
$z^{+}>56$
, the minimal profiles diverge from the full channel, but the inset in figure 5(a) shows that the roughness function,
${\rm\Delta}U^{+}$
, remains the same for both full and minimal channels. In other words, the minimal channel is sufficient for predicting the roughness function. The method works in this case because the top of the roughness elements at
$z^{+}=k^{+}\approx 10$
is well within the unconfined region
$z^{+}\approx 0.4\,L_{y}^{+}\approx 56$
, and the roughness wavelength
${\it\lambda}_{y}^{+}~({\approx}71)$
is smaller than
$L_{y}^{+}$
. An issue that often arises is the appropriate wall-normal location to evaluate the roughness function. In high-Reynolds-number flows, this location is frequently taken to be somewhere in the log layer. Here we observe that
$U_{s}^{+}-U_{r}^{+}$
is relatively insensitive to this location as long as
${\rm\Delta}U^{+}$
is evaluated above the roughness height. The insets in figures 3 and 5 demonstrate that this observation is consistent for all our results, for both explicit and modelled roughness.
Figure 5(b) shows the mean velocity profiles for the minimal channel at a fixed blockage ratio
$h/k=18$
, where
$k^{+}$
is increased from 10 to 60 by increasing the friction Reynolds number from
$\mathit{Re}_{{\it\tau}}\approx 180$
to 1080. For the lowest
$k^{+}\approx 10$
, the span of the minimal channel is chosen to be
$L_{y}^{+}\approx 141$
in order to properly capture the near-wall cycle and accommodate two periods of the sinusoidal roughness. In the cases with the higher roughness Reynolds numbers,
$k^{+}\approx 20$
, 30 and 60, the minimal channel only captures one period of sinusoidal roughness, that is,
${\it\lambda}_{y}=L_{y}$
. The corresponding predicted roughness functions
${\rm\Delta}U^{+}$
are shown in the inset of figure 5(b).
In order to assess the efficacy of the present method for an explicitly represented sinusoidal roughness, we plot the Hama roughness function versus the equivalent sand-grain roughness
$k_{s}^{+}$
in figure 6. Fitting the data in the fully rough regime, we are able to characterise this particular sinusoidal roughness with
${\it\lambda}_{x}={\it\lambda}_{y}\approx 7.1\,k$
to find that
$k_{s}/k\approx 4.1$
. The minimal-channel predictions, shown by the open markers, collapse onto full pipe-flow simulations of the same roughness previously studied by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015), albeit at various pipe blockage ratios
$R/k$
, as shown by the filled circles in figure 6. This collapse further validates the ability of the minimal channel simulations to recover the same estimate for
$k_{s}$
as given by full-size simulations.
It is worth pointing out that recent direct numerical simulations have also been simulated in the fully rough regime (Busse & Sandham Reference Busse and Sandham2012; Yuan & Piomelli Reference Yuan and Piomelli2014). An important distinction here is that these previous cases used domain sizes that are deemed sufficiently large to capture the full range of turbulent motions. In contrast, the present minimal simulations have a severely confined outer flow with much lower intensities and a radically altered wake profile. The potential saving in computational cost for the minimal simulations could be redeployed to explicitly represent the roughness (see the present Set 2 simulations), without the need for roughness models.
4. Conclusions
We have presented a novel, fast and direct method for characterising the hydraulic resistance of any given surface roughness. The way in which a particular roughness transitions from the hydraulically smooth to the fully rough regime is, to the first approximation, described by how the roughness geometry interacts with the near-wall flow. The method presented herein shows that this interaction, as far as the mean drag is concerned, is accurately captured using minimal-span channels. There are other roughness effects, such as the change in turbulent kinetic energy, and it remains to be seen how these effects are represented by the minimal-channel technique. We have validated the present method with a specific sinusoidal roughness and showed that the method is able to fully characterise the hydraulic (drag) behaviour of roughness, not only in the transitionally rough regime but also in the fully rough regime, yielding an estimate of
$k_{s}$
which closely matches that given by full-span simulations.
The dynamic drag characterisation of a rough surface is encapsulated in the Hama roughness function, which we show can be accurately determined using minimal-span channels. The fact that
${\rm\Delta}U^{+}$
is plotted against
$k^{+}$
(or
$k_{s}^{+}$
) in the literature and not against
$\mathit{Re}_{{\it\tau}}$
acknowledges that
${\rm\Delta}U^{+}$
, to a large extent, depends only on the roughness-affected near-wall flow. The minimal-span channel is a method for simulating this near-wall flow of thickness
$O(k)$
without resolving the outer scale
$h\gg k$
, thereby breaking the curse of the (outer) Reynolds number. The savings in computational cost for the present method are possible because capturing only the near-wall flow requires far fewer grid points than capturing the full flow. For example, the full-span DNS of Bernardini et al. (Reference Bernardini, Pirozzoli and Orlandi2014) at
$\mathit{Re}_{{\it\tau}}\approx 4000$
requires
$N_{x}\times N_{y}\times N_{z}=8192\times 4096\times 1024=3.4\times 10^{10}$
grid points. In contrast, the minimal-span (
$L_{y}^{+}\approx 463$
) DNS at the same
$\mathit{Re}_{{\it\tau}}$
(table 1), which is sufficient to reach the fully rough regime, requires only
$N_{x}\times N_{y}\times N_{z}=384\times 96\times 1024=3.8\times 10^{7}$
grid points, or a three-orders-of-magnitude reduction in number of grid points. The actual wall-clock time will, of course, depend on the code and parallelisation, as well as the machine. As an indicative comparison using the same code and machine, the full-span (
$L_{y}/h={\rm\pi}$
)
$\mathit{Re}_{{\it\tau}}\approx 950$
simulation (table 1) requires
$2600\,\text{CPU hours}/(h/U_{{\it\tau}})$
, while the minimal-span (
$L_{y}^{+}\approx 112$
)
$\mathit{Re}_{{\it\tau}}\approx 950$
simulation (table 1) requires only
$42\,\text{CPU hours}/(h/U_{{\it\tau}})$
. These figures amount to a two-orders-of-magnitude reduction in CPU hours.
The present method can be used to characterise the drag characteristics of many surfaces very quickly. Such a powerful tool enables the researcher to focus on how the geometry of roughness affects the turbulent flow, instead of focusing on setting up expensive and time-consuming experiments, physical or numerical. Indeed, the present method even enables researchers to reassess the physics underlying the success and failure of previously proposed geometrical factors. Recall that the minimal-span method retains the same benefit from a full-span DNS in that the friction velocity formed from both viscous and pressure contributions is directly measured.
In some ways, the present idea is not entirely new. In the large-eddy simulation subgrid-scale modelling methodology, the geometry-dependent large eddies are directly simulated whilst the subgrid small eddies are assumed to be universal and are modelled or understood through the energy-cascade phenomenology of Richardson, Taylor and Kolmogorov. Here, this reasoning is reversed for the roughness problem, in which the geometry-dependent small eddies are directly simulated whilst the universal ‘supergrid’ large eddies are assumed to be universal and are modelled or understood through the outer-layer similarity hypothesis of Townsend, a kind of small-eddy simulation, a term coined by Jiménez (Reference Jiménez2003).
We have also provided guidelines on selecting the minimal channel. To obtain an accurate prediction of
${\rm\Delta}U^{+}$
, the roughness element must be submerged within unconfined near-wall flow, which yields the first two constraints:
$L_{y}>k/0.4$
and
$L_{y}^{+}>100$
respectively. Potentially, if
$k$
were very large, the first of these constraints could lead to increasingly large ‘minimal’ channels. However, in practice it seems unlikely that these simulations would need to be conducted for
$k_{s}^{+}\gg 100$
in order to establish the fully rough regime, and hence
$L_{y}$
defined by this constraint should remain manageable. More problematic is the constraint that
$L_{y}>{\it\lambda}_{y}$
. For heterogeneous or sparse rough surfaces, the spanwise length scale or repeating unit
${\it\lambda}_{y}$
could potentially be very large in terms of viscous scaling, presenting an obvious Achilles’ heel for this proposed method, although this would be equally true for full-span simulations. Even for surfaces which have a small
${\it\lambda}_{y}$
yet where the roughness is randomly arranged (such as a sanded or painted surface), the minimal span would need to be sufficiently large to contain a statistically representative sample of the rough surface. Even so, this
${\it\lambda}_{y}$
would in most cases be significantly less than the
${\rm\pi}h$
box width often required of full-span simulations. Within these constraints, this method could potentially be extremely useful for all rough surfaces that exhibit homogeneity over a relatively small scale. Examples could include surface finishes obtained from machining, painting or spray coating. One could also potentially use minimal channels scaled in this way to investigate drag reducing rough surfaces such as riblets. Indeed, other passive or active surfaces where the primary effect of the surface is an alteration of the near-wall structure, or near-wall slip, could also be investigated using minimal channels (e.g. super-hydrophobic surfaces, compliant walls, porosity, acoustic liners).
We can also point to certain possible improvements or extensions of this technique in the future. One issue with the minimal channels is that the centreline velocity becomes very high (relative to
$U_{{\it\tau}}$
), which adds computational expense through required reductions in the time step to satisfy the Courant–Friedrichs–Lewy (CFL) condition. A possible solution here would be to add a body-forcing term at the channel centreline to appropriately manage the magnitude of this velocity. An alternative is to use a simple eddy-viscosity term that is active sufficiently far above the roughness elements to account for absent eddies that are wider than the minimal span, which would reduce the centreline velocity. It is also possible that the number of grid points could be further halved by using an open channel (Scotti Reference Scotti2006). From a different approach, it could be the case that wider minimal channels may be more efficient in terms of obtaining statistically converged statistics for given CPU hours (due to reductions in the ‘burstiness’ associated with minimal channels where
$L_{y}^{+}$
is close to 100). Finally, and perhaps most promising, we are keen to explore the capability of a single temporal simulation, starting from the highest
$\mathit{Re}_{{\it\tau}}$
, and reducing the bulk velocity to map
${\rm\Delta}U^{+}$
versus
$k^{+}$
(and hence obtain
$k_{s}^{+}$
) in a single numerical experimental sweep.
Acknowledgements
This research was undertaken on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government, and also on the Victorian Life Science Computational Institute (VLSCI). The authors gratefully acknowledge the financial support of the Australian Research Council (ARC).
