The minimal-span channel for rough-wall turbulent flows

TL;DR

This paper investigates the minimal-span channel approach for simulating rough-wall turbulent flows, optimizing domain size, boundary conditions, and computational cost to accurately capture near-wall dynamics with reduced resources.

Contribution

It introduces a refined minimal-span channel framework, including domain length criteria, slip wall boundary conditions, and a cost estimation method for efficient simulations of rough-wall turbulence.

Findings

Optimal streamwise domain length is three times the spanwise width or 1000 viscous units.

Half-height channels with slip walls replicate near-wall behavior with fewer grid points.

A cost model relates simulation duration to statistical uncertainty in roughness function.

Abstract

Roughness predominantly alters the near-wall region of turbulent flow while the outer layer remains similar. This makes it a prime candidate for the minimal-span channel, which only captures the near-wall flow by restricting the spanwise channel width to be of the order of a few hundred viscous units. Recently, Chung et al. (J. Fluid Mech., vol. 773, 2015, pp. 418-431) showed that a minimal-span channel can accurately characterise the hydraulic behaviour of roughness. Following this, we aim to investigate the fundamental dynamics of the minimal-span channel framework with an eye towards further improving performance. The streamwise domain length of the channel is investigated with the minimum length found to be three times the spanwise width or 1000 viscous units, whichever is longer. A half-height (open) channel with slip wall is shown to reproduce the near-wall behaviour seen in a standard channel, but with half the number of grid points. Next, a forcing model is introduced into the outer layer of a half-height channel to reduce the high streamwise velocity. Finally, an investigation is conducted to see if varying the roughness Reynolds number with time is a feasible method for obtaining the full hydraulic behaviour of a rough surface. An empirical costing argument is developed to determine the cost in terms of CPU hours of minimal-span channel simulations a priori. This argument involves counting the number of eddy lifespans in the channel, which is then related to the statistical uncertainty of the streamwise velocity. For a given statistical uncertainty in the roughness function, this can then be used to determine the simulation run time. Following this, a finite-volume code with a body-fitted grid is used to determine the roughness function for square-based pyramids using the above insights. Good agreement with the literature for the same roughness geometry is observed.

Tables

Table 1. Table 1: Description of the simulations performed with the finite difference code to investigate streamwise domain length. All simulations done as both smooth wall and modelled rough wall. Entries: R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} , nominal friction Reynolds numbers; L 𝐿 L , domain length; N 𝑁 N , number of grid points; Δ Δ \Delta , grid-spacing in viscous units, subscript w 𝑤 w and h ℎ h refers to the wall-normal grid spacing at the wall and channel centre; L z / h = 2 subscript 𝐿 𝑧 ℎ 2 L_{z}/h=2 denotes standard-height (two no-slip walls), Δ ​ t S + Δ superscript subscript 𝑡 𝑆 \Delta t_{S}^{+} and Δ ​ t R + Δ superscript subscript 𝑡 𝑅 \Delta t_{R}^{+} is the smooth- and rough-wall time step. Roughness height k = h / 40 𝑘 ℎ 40 k=h/40 so k + ≈ 15 superscript 𝑘 15 k^{+}\approx 15 at R ​ e τ = 590 𝑅 subscript 𝑒 𝜏 590 Re_{\tau}=590 .

ID	$R{e}_{\tau }$	$\frac{L_x}{h}$	$L_x^{+}$	$L_y^{+}$	$\frac{L_z}{h}$	$N_x$	$N_y$	$N_z$	$\Delta x^{+}$	$\Delta y^{+}$	$\Delta z_w^{+}$	$\Delta z_h^{+}$	$\Delta t_S^{+}$	$\Delta t_R^{+}$
Full-span channel
FS	590	$2\pi$	3707	1854	2	384	384	256	9.7	4.8	0.04	7.2	0.26	0.29
Minimal-span channel with varying streamwise length
MS1	590	$0.1\pi$	190	118	2	24	24	256	7.7	4.9	0.04	7.2	0.28	0.32
MS2	590	$0.16\pi$	300	118	2	48	24	256	6.2	4.9	0.04	7.2	0.24	0.27
MS3	590	$0.25\pi$	460	118	2	48	24	256	9.7	4.9	0.04	7.2	0.36	0.42
MS4	590	$0.5\pi$	930	118	2	96	24	256	9.7	4.9	0.04	7.2	0.36	0.42
MS5	590	$\pi$	1850	118	2	192	24	256	9.7	4.9	0.04	7.2	0.36	0.41
Minimal-span channel with varying streamwise length (wider span)
MSW1	590	$0.1\pi$	190	354	2	24	72	256	7.7	4.9	0.04	7.2	0.37	0.45
MSW2	590	$0.16\pi$	300	354	2	48	72	256	6.2	4.9	0.04	7.2	0.31	0.40
MSW3	590	$0.25\pi$	460	354	2	48	72	256	9.7	4.9	0.04	7.2	0.46	0.52
MSW4	590	$0.5\pi$	930	354	2	96	72	256	9.7	4.9	0.04	7.2	0.43	0.46
MSW5	590	$\pi$	1850	354	2	192	72	256	9.7	4.9	0.04	7.2	0.38	0.40

Table 2. Table 2: Full-span and minimal-span channel simulations, for standard-height (two no-slip walls, L z / h = 2 subscript 𝐿 𝑧 ℎ 2 L_{z}/h=2 ) and half-height (one no-slip and one slip wall, L z / h = 1 subscript 𝐿 𝑧 ℎ 1 L_{z}/h=1 ) channels. Entries are same as table 1 .

ID	$R{e}_{\tau }$	$\frac{L_x}{h}$	$L_x^{+}$	$L_y^{+}$	$\frac{L_z}{h}$	$N_x$	$N_y$	$N_z$	$\Delta x^{+}$	$\Delta y^{+}$	$\Delta z_w^{+}$	$\Delta z_h^{+}$	$\Delta t_S^{+}$	$\Delta t_R^{+}$
Full-span channel
FS	590	$2\pi$	3707	1854	2	384	384	256	9.7	4.8	0.04	7.2	0.26	0.29
FSH	590	$2\pi$	3707	1854	1	384	384	128	9.7	4.8	0.04	7.2	0.29	0.32
Minimal-span channel
MS6	590	$2\pi$	3707	118	2	384	24	256	9.7	4.9	0.04	7.2	0.35	0.40
MSH6	590	$2\pi$	3707	118	1	384	24	128	9.7	4.9	0.04	7.2	0.36	0.41

Table 3. Table 3: Minimal-span simulations with outer-layer damping. Refer to table 1 for definitions. z d + superscript subscript 𝑧 𝑑 z_{d}^{+} indicates location where damping starts (figure 1 ), cases with no damping indicated by a hyphen.

ID	$R{e}_{\tau }$	$\frac{L_x}{h}$	$L_x^{+}$	$L_y^{+}$	$\frac{L_z}{h}$	$N_x$	$N_y$	$N_z$	$\Delta x^{+}$	$\Delta y^{+}$	$\Delta z_w^{+}$	$\Delta z_h^{+}$	$z_d^{+}$	$\Delta t_S^{+}$	$\Delta t_R^{+}$
Minimal-span, half-height channel with outer-layer damping
MSHD1	590	$2\pi$	3707	118	1	384	24	128	9.7	4.9	0.04	7.2	100	0.48	0.59
MSHD2	590	$2\pi$	3707	118	1	384	24	128	9.7	4.9	0.04	7.2	150	0.46	0.54
MSHD3	590	$2\pi$	3707	118	1	384	24	128	9.7	4.9	0.04	7.2	200	0.44	0.51
MSHD4	590	$2\pi$	3707	118	1	384	24	128	9.7	4.9	0.04	7.2	300	0.41	0.47
MSH6	590	$2\pi$	3707	118	1	384	24	128	9.7	4.9	0.04	7.2	-	0.36	0.41
Minimal-span, half-height channel with outer-layer damping (wider span)
MSWHD1	590	$2\pi$	3707	236	1	384	48	128	9.7	4.9	0.04	7.2	200	0.43	0.48
MSWHD2	590	$2\pi$	3707	236	1	384	48	128	9.7	4.9	0.04	7.2	250	0.43	0.47
MSWHD3	590	$2\pi$	3707	236	1	384	48	128	9.7	4.9	0.04	7.2	300	0.42	0.46
MSWH6	590	$2\pi$	3707	236	1	384	48	128	9.7	4.9	0.04	7.2	-	0.40	0.43
Minimal-span, half-height channel with outer-layer damping (increased $R{e}_{\tau }$)
MSHRD1	2000	$0.59\pi$	3707	300	1	384	60	384	9.7	5.0	0.02	8.2	300	0.25	0.37
MSHRD2	2000	$0.59\pi$	3707	300	1	384	60	384	9.7	5.0	0.02	8.2	400	0.24	0.35
MSHR6	2000	$0.59\pi$	3707	300	1	384	60	384	9.7	5.0	0.02	8.2	-	0.24	0.31

Table 4. Table 4: Initial conditions for temporal sweep in a minimal-span channel. Refer to table 1 for definitions.

ID	$R{e}_{\tau }$	$\frac{L_x}{h}$	$L_x^{+}$	$L_y^{+}$	$\frac{L_z}{h}$	$N_x$	$N_y$	$N_z$	$\Delta x^{+}$	$\Delta y^{+}$	$\Delta z_w^{+}$	$\Delta z_h^{+}$	$\Delta t_S^{+}$	$\Delta t_R^{+}$
Temporal sweep (at initial conditions)
MS6	590	$2\pi$	3707	354	2	384	80	256	9.6	4.4	0.04	7.2	0.38	0.43

Table 5. Table 5: Description of the simulations performed with the finite volume code for the square-based pyramids. k t + superscript subscript 𝑘 𝑡 k_{t}^{+} is the trough-to-peak pyramid height, λ + superscript 𝜆 \lambda^{+} the pyramid wavelength, α 𝛼 \alpha the pyramid edge angle, and C f , m ​ o ​ d = 2 ​ U τ 2 / U b , f ​ u ​ l ​ l 2 subscript 𝐶 𝑓 𝑚 𝑜 𝑑 2 superscript subscript 𝑈 𝜏 2 superscript subscript 𝑈 𝑏 𝑓 𝑢 𝑙 𝑙 2 C_{f,mod}=2U_{\tau}^{2}/U_{b,full}^{2} the predicted full-span skin friction coefficient (see text). Other entries are same as table 1 . All simulations conducted using a half-height channel with no outer-layer damping. SF09 refers to Schultz & Flack ( 2009 ) , HKS11 to Hong et al. ( 2011 ) , and DCO16 to Di Cicca & Onorato ( 2016 ) .

Case	$R{e}_{\tau }$	$k_t^{+}$	${\lambda }^{+}$	$\alpha$	$L_x^{+}$	$L_y^{+}$	$N_x$	$N_y$	$N_z$	$\Delta x^{+}$	$\Delta y^{+}$	$\Delta z_w^{+}$	$\Delta z_h^{+}$	$C_{f,mod}$	$\Delta U^{+}$
Sm_198	834	-	-	-	1180	197	288	48	320	4.1	4.1	0.14	6.2	0.0053	-
Sm_313	850	-	-	-	1201	316	288	72	320	4.2	4.4	0.14	6.3	0.0051	-
Sm_396	851	-	-	-	1204	401	288	96	320	4.2	4.2	0.14	6.3	0.0050	-
40_198	848	40.4	200	22	1200	200	288	48	320	4.2	4.2	0.14	6.3	0.015	8.0
40_198_f	837	39.2	197	22	1184	197	504	84	560	2.3	2.3	0.24	3.6	0.015	7.9
40_198_2	879	41.7	207	22	1243	414	288	96	320	4.3	4.3	0.15	6.5	0.015	8.3
63_313	883	66.3	329	22	1249	329	288	72	320	4.3	4.6	0.24	6.5	0.021	10.2
80_396	868	82.4	409	22	1227	409	288	96	320	4.3	4.3	0.24	6.4	0.024	11.0
SF09	4020	39.9	198	16	-	-	-	-	-	-	-	-	-	-	6.4
SF09	3900	62.7	310	16	-	-	-	-	-	-	-	-	-	-	7.0
SF09	4290	87.5	433	16	-	-	-	-	-	-	-	-	-	-	8.7
DCO16	670	39	147	28	-	-	-	-	-	-	-	-	-	-	8.9
HKS11	3520	63.3	442	16	-	-	-	-	-	-	-	-	-	-	8.2

Equations

$$\nabla \cdot \mathbf{u}=0,\frac{\partial \mathbf{u}}{\partial t}+\nabla \cdot \left(\mathbf{uu}\right)=-\frac{1}{\rho }\nabla p+\nabla \cdot \left(\nu \nabla \mathbf{u}\right)+\mathbf{F}+\mathbf{G}+\mathbf{K},$$

$$F(z,k)=\{\begin{array}{ll}1,& \text{if z<k or 2h−k<z}.\\ 0,& \text{otherwise}.\end{array}$$

$$K_{i}=-\gamma\Gamma(z,z_{d})\Big{(}u_{i}(x,y,z,t)-\langle u_{i}\rangle(z_{d},t)\Big{)}\Big{|}u_{i}(x,y,z,t)-\langle u_{i}\rangle(z_{d},t)\Big{|},$$

$$\Gamma (z,z_d)=\{\begin{array}{ll}1,& \text{if z>zd and z<2h−zd}.\\ 0,& \text{otherwise}.\end{array}$$

$$\frac{\mathrm{d}{U}_b}{\mathrm{d}t}=A=\frac{U_{b,end}-U_{b,start}}{\Delta T}=const.,$$

$${\Delta }_p=\frac{\nu }{\rho U_{\tau }^3}\frac{\mathrm{d}p}{\mathrm{d}x}.$$

$$N_{\text{instant}}=\frac{L_x}{{\lambda }_{x,z_c}}\frac{L_y}{{\lambda }_{y,z_c}}\cdot \frac{L_z}{h}=\frac{L_x}{7.5z_c}\frac{L_y}{2.5z_c}\cdot \frac{L_z}{h}.$$

$$C^{\star }=\frac{T_{sim}}{T_b}{N}_{\text{instant}}=\frac{T_{sim}}{T_b}\cdot \frac{L_z}{h}\cdot \frac{L_x}{7.5z_c}\frac{L_y}{2.5z_c}.$$

$$U(z,T_{sim})=\frac{1}{N}\sum _{i=0}^{N-1}⟨u(z,t_i)⟩$$

$$T_0=1+2\sum _{k=1}^{N-1}\left(1-\frac{k}{N}\right)\rho (k)$$

$$s^2=\frac{1}{N-T_0}\sum _{i=0}^{N-1}{\left(⟨u(z,t_i)⟩-U(z,T_{sim})\right)}^2.$$

$$\Delta U^{+}\pm 2.77{ϵ}^{+}\approx \Delta U^{+}\pm 2.77K(C^{\star }{)}^{-1/2}\approx \Delta U^{+}\pm 91.4(C^{\star }{)}^{-1/2}/z_c^{+}.$$

$$\text{CPU hours}=\frac{T_{sim}^{+}}{\Delta t^{+}}{\mathcal{T}}_1{N}_{CPU},$$

Full text

1 Introduction

Conventional Direct Numerical Simulations (DNS) of wall-bounded turbulent flows represent a challenging computational problem, as both small and large scales need to be represented. The former requires a fine grid to resolve the small viscous scales, while the latter requires a large domain to capture the large outer-layer motions. However, pioneering work by Jiménez & Moin (1991) and Hamilton et al. (1995) into minimal-flow units showed that a small computational domain can be used to exclusively capture the turbulent near-wall cycle, independent of the large outer scales. This is achieved by simply restricting the domain of the channel to a small size, where the spanwise and streamwise lengths are prescribed in terms of viscous units. Jiménez & Moin (1991) showed that turbulence could be maintained in the form of the near-wall cycle of the buffer layer when the spanwise domain width was only on the order of $100\nu/U_{\tau}$; here $\nu$ is the kinematic viscosity and $U_{\tau}=\sqrt{\tau_{w}/\rho}$ is the friction velocity, defined using the wall shear stress $\tau_{w}$ and the fluid density $\rho$ . This was further supported by future studies into minimal-flow units (Jiménez & Pinelli, 1999; Flores & Jiménez, 2010; Hwang, 2013; Chung et al., 2015; Hwang & Bengana, 2016). In particular the work of Flores & Jiménez (2010) showed that the minimal-span channel can also capture the logarithmic layer of turbulent flows.

An important aspect of minimal-span channels is that the near-wall flow is accurately captured up to a critical wall-normal location, $z_{c}$. Above this point, the streamwise velocity increases compared to a full-span channel. This unphysical increase occurs as the narrow spanwise domain width of the minimal-span channel acts as a filter which limits the largest spanwise scale of energy-containing eddies. The flow does, however, retain turbulent scales smaller than the spanwise domain width, so it is not laminar (Jiménez & Pinelli, 1999). The critical value where the minimal-span channel departs from the full-span channel has been shown to scale with the spanwise domain width, $z_{c}\approx 0.4L_{y}$ (Chung et al., 2015), although a constant of 0.3 is suggested in Flores & Jiménez (2010) and 0.3–0.4 in Hwang (2013). Following Flores & Jiménez (2010), we will refer to the flow below $z_{c}$ as ‘healthy’ turbulence, as it is in a full-span channel.

In the context of roughness, the central question is how the geometry of a rough surface is related to its hydraulic behaviour; namely, what value the (Hama) roughness function, $\Delta U^{+}$, takes (Hama, 1954). This quantity reflects the flow retardation, or velocity shift, that the roughness imposes on the flow when compared to a smooth wall, for matched friction Reynolds numbers, $Re_{\tau}=U_{\tau}h/\nu$ (here $h$ is the half-channel height, boundary layer thickness, or pipe radius). For a given surface, $\Delta U^{+}$ can be obtained from various semi-empirical models and approximations, or directly from full-scale experiments and numerical simulations. The former are of varying accuracy and depend on the model selected and the rough surface in question, while the latter can be prohibitively expensive. In particular, full-scale numerical simulations suffer from the drawbacks mentioned above, in which both small- and large-scale motions need to be captured. However, roughness is thought to primarily alter only the near-wall flow, in a region called the roughness sublayer which typically extends 3–5$k$ from the wall, where $k$ is some characteristic height of the roughness (Raupach et al., 1991). Well outside the roughness sublayer the flow is only changed insofar as it depends on the friction velocity $U_{\tau}$. This is the basis of Townsend’s outer-layer similarity hypothesis (Townsend, 1976) which has received significant attention and has been supported by several rough-wall studies (Flack et al., 2005; Leonardi & Castro, 2010; Chan et al., 2015) and also numerical simulations with modified boundary conditions (Flores & Jiménez, 2006; Mizuno & Jiménez, 2013; Chung et al., 2014). If we assume Townsend’s hypothesis holds, then it follows that the roughness only alters the near-wall region of the flow which therefore makes it a prime candidate for use in the minimal-span channel framework.

The equivalent sandgrain roughness, $k_{s}$, is a single dynamic parameter that is used to describe the hydraulic behaviour of roughness. It is defined so that the roughness function collapses for all data in the fully rough regime, when the friction factor no longer depends on the bulk-velocity Reynolds number (Jiménez 2004, who called this $k_{s\infty}$). However, each rough surface will have a unique behaviour in the transitionally rough regime. It is an expensive process to determine the full hydraulic behaviour, as a range of simulations need to be conducted for the different roughness Reynolds numbers, each requiring a different body-fitted grid and its own initialisation period. It would therefore be desirable to conduct a simulation in which only a single computational grid is used, and the bulk velocity is changed over time to sweep through a range of roughness Reynolds numbers. This is a similar approach to how experimental studies are performed in that a single rough surface with fixed $k/h$ is tested at multiple flow speeds, so that the roughness Reynolds number varies. In order to obtain statistics at a desired friction Reynolds number in a temporally evolving flow, statistics are averaged over a small window in which the instantaneous friction Reynolds number is close to the desired one. Within this window, the flow is assumed to be quasi steady. This would therefore generate a near-continuous curve of $\Delta U^{+}$ versus $k^{+}$, rather than the conventional (steady) approach which only generates a few data points. An important consideration which will be investigated here is whether acceleration effects from the changing mass flux become significant and distort the estimation of $\Delta U^{+}$, that is, how quickly can the bulk velocity be varied such that the quasi-steady assumption remains approximately valid.

Recently, we have applied the idea of minimal-span channels to the roughness problem in a proof-of-concept study (Chung et al., 2015). This study demonstrated the feasibility of using the minimal-span channel to accurately compute the roughness function, when compared to full-span channels. However, only the spanwise width was investigated, with the roughness function being the primary quantity of interest. In this paper, we aim to further investigate the fundamental dynamics of the minimal-span channel framework for roughness. This will include analysing higher order flow statistics to identify and understand the essential features of the minimal-span channel. Firstly, the streamwise domain length of the minimal channel is investigated in §3.1 to identify the minimum streamwise length required to maintain healthy turbulence. Second, given that the outer-layer of minimal channels are inherently unphysical, then alterations to this region should not alter the healthy near-wall flow. To assess this claim, two alterations consisting of a half-height (open) channel and outer-layer damping are investigated in §3.2 and §3.3, respectively. Finally, an investigation in §3.4 is conducted in which the bulk velocity (and hence roughness Reynolds number) is varied with time, to see if the full hydraulic behaviour of the rough surface can be obtained in one unsteady simulation. A costing argument is then developed in §4 using the characteristic time and length scales of eddies, so that the CPU hours can be estimated a priori. The insights gained in this paper are then used in §5 to simulate the flow over square-based pyramids using a finite-volume code with a body-fitted grid, which has been experimentally studied in the literature by Schultz & Flack (2009) and others.

2 Numerical Method

The majority of Direct Numerical Simulations in this paper are conducted using the same finite difference code as in Chung et al. (2015), which uses a fully conservative fourth-order staggered-grid scheme. Time integration is performed using the third-order low-storage Runge–Kutta scheme of Spalart et al. (1991), and the fractional-step method of Kim & Moin (1985) is used. The Navier–Stokes equations are

[TABLE]

where $\mathbf{u}=(u,v,w)$ is the velocity in the streamwise ($x$), spanwise ($y$) and wall-normal ($z$) directions, $t$ is time, and $p$ is pressure. $\mathbf{G}=G_{x}(t)\mathbf{i}$ is the spatially invariant, time-varying streamwise forcing term which drives the flow at constant mass flux. The flow is solved in a reference frame in which the mass flux is zero at all times, although all equations in this paper are given in the stationary-wall reference frame. Solving in a zero-mass-flux reference frame permits a larger computational time step (Lundbladh et al., 1999), as well as reducing high-wavenumber convective disturbances produced by finite difference schemes (Bernardini et al., 2013). The roughness model, $\mathbf{F}=-\alpha F(z,k)u|u|\mathbf{i}$, is based on the work of Busse & Sandham (2012), which applies a forcing in the streamwise direction that opposes the flow. The roughness factor $\alpha=1/(40k)$ is kept constant for all simulations, with the roughness height $k=h/40$. The function $F(z,k)$ is a simple step function which applies the roughness forcing adjacent to the top and bottom no-slip channel walls at $z=0$ and $z=2h$,

[TABLE]

This roughness model removes any spanwise or streamwise roughness length scales from the problem and allows for simpler computations as the same smooth-wall grid can be used. This is a similar model to Borrell (2015), although there the author used an $\mathbf{F}\propto-u\mathbf{i}$ scaling. The final term in (1), $\mathbf{K}$, is a forcing function which damps the velocity in the outer layer of the minimal channel. It has the form

[TABLE]

(no summation over $i$) where angled brackets denote the spatial average of the instantaneous velocity over a wall-parallel plane. The factor $\gamma$ has units of inverse time and in this study is set according to $\gamma h/U_{\tau}\approx 1$. The parameter $\Gamma(z,z_{d})$ is similar to the function $F(z,k)$ of the roughness forcing model in that it indicates where the damping is applied, namely in the outer layer of the channel,

[TABLE]

The value of $z_{d}$ should be greater than the critical wall-normal location, $z_{c}$, of the minimal-span channel, so that the forcing does not contaminate the healthy turbulence of the near-wall flow. A step function for $\Gamma$ is used to minimise the parameter space, as well as to ensure the location of $z_{d}$ is unambiguous. Several different values of $z_{d}$ are tested and will be presented in the following section. In (3), the term $\langle u_{i}\rangle(z_{d},t)$ is the wall-parallel spatially averaged velocity at $z_{d}$. This is present to reduce the magnitude of the streamwise velocity; the minimal channel has a very high streamwise velocity which presents a time step restriction from the CFL number. The forcing is in this form so that the velocity is forced to remain at the same level as at $z=z_{d}$. This outer-layer damping is similar to the masks employed in Jiménez & Pinelli (1999), which damp fluctuations in the outer-layer region of minimal channels. The current forcing model simply sets the streamwise velocity to be approximately that at $z=z_{d}$, however it will be shown to work reasonably well for the purpose of obtaining $\Delta U^{+}$. Note that the average velocity in the spanwise and wall-normal directions is zero, so that the forcing in these directions can simply be $K_{i}=-\gamma\Gamma(z,z_{d})u_{i}|u_{i}|$, for $i=2,3$. Figure 1 shows the two forcing regions that are employed in the current study.

The streamwise and spanwise grid is evenly spaced, while the wall-normal grid is stretched with a cosine mapping. Periodic boundary conditions are applied in the spanwise and streamwise directions. In the case of standard-height channels, no-slip walls are located at $z=0$ and $z=2h$. For clarity, the word ‘full’ will be used to refer to the span (full span), while the case with two no-slip walls will be referred to as ‘standard’ height. For half-height (open) channels, a no-slip wall is still positioned at $z=0$, however now the top domain surface is a free-slip wall with $\partial u/\partial z=\partial v/\partial z=0$, positioned at $z=h$. This slip wall maintains the impermeability constraint, $w=0$. The spanwise width, $L_{y}$, of the minimal-span channel satisfies the guidelines of Chung et al. (2015), namely that $L_{y}\gtrsim\max(100\nu/U_{\tau},k/0.4,\lambda_{r,y})$, where $\lambda_{r,y}$ is the spanwise length scale of roughness elements. Because the homogeneous roughness forcing model has no spanwise length scale then the final constraint can be ignored. Simulations are conducted at a friction Reynolds number of $Re_{\tau}\approx 590$, with a few additional simulations conducted at $Re_{\tau}\approx 2000$. Relevant simulations are introduced at the start of each section for each of the investigations conducted in this study.

2.1 Temporal sweep

As discussed in the introduction, instead of conducting multiple steady simulations of different roughness Reynolds numbers in order to determine the equivalent sandgrain roughness, a single unsteady simulation is to be performed in which the bulk velocity is varied. In particular, we will investigate different rates of change of the bulk velocity and the effects that this acceleration has on the flow. The sweep will start with the highest friction number to be tested, $Re_{\tau,start}$ and will then be decelerated via an adverse pressure gradient. This ensures an adequate grid resolution at the start of simulation. At the final friction Reynolds number, $Re_{\tau,end}<Re_{\tau,start}$, the grid resolution would be finer than necessary, which for the current simulations would have 4–5 times more cells than if a conventional steady simulation were conducted at this final Reynolds number.

Presently, we vary the bulk velocity $U_{b}(t)$ with time,

[TABLE]

where $U_{b,end}$ is the desired end point of the simulation (the bulk velocity corresponding to $Re_{\tau}=180$), and $\Delta T$ is a time scale which should be sufficiently long enough that acceleration effects are not significant to the flow. The initial bulk velocity, $U_{b,start}$, is set corresponding to a friction Reynolds number of $Re_{\tau,start}=590$ and different values of $A$ are tested. An initial value of $A$ is selected such that over the course of the sweep $\Delta TU_{\tau 0}/h=30\Rightarrow U_{\tau 0}^{2}/hA\approx 30$, with subsequent runs using $2A$ and $4A$. In these cases of a higher rate of change of $U_{b}$, multiple sweeps are run from different initial conditions with the results ensemble averaged, to ensure that statistics are obtained over the same amount of simulation time. For example, the case where the gradient is quadrupled, four sweeps are conducted with four different initial conditions. This is visualised in figure 2(a), which shows the change in the friction Reynolds number as a function of time. An important consideration with this technique is that the friction velocity and hence the normalised spanwise domain width $L_{y}^{+}$ will vary with time. The domain width must remain larger than 100 viscous units at all times to ensure the turbulent flow is sustained, which may necessitate multiple stages in the sweep if $L_{y}^{+}$ becomes too small. Each stage in the sweep should be set up such that $L_{y}^{+}\gtrsim 100$.

A dimensionless measure of the flow acceleration is the pressure gradient parameter,

[TABLE]

Various experimental studies of decelerating boundary layers show that mild decelerations of $\varDelta_{p}<0.01$ have small ($<3\%$) errors in calculating $U_{\tau}$ from methods based on the assumption of zero pressure gradients (Patel, 1965; Jones et al., 2001). While we are interested in channel flow as opposed to boundary layers, this value should still be useful in providing an indication as to whether acceleration effects will be significant. Seddighi et al. (2014), meanwhile, compared a step acceleration with a slow ramp up in which the bulk velocity was linearly varied in a channel such that $-\varDelta_{p}\approx 0.73$ (favourable pressure gradient). Acceleration effects were still seen in this slow ramp up. Note that for steady channel flow, the pressure gradient parameter is $\varDelta_{p}=-1/Re_{\tau}<0$, which for the current simulations at $Re_{\tau}=590$ takes a value of $-0.0017$. In the current study, three different sweep rates are studied. As shown in figure 2(b), the linear variation in $U_{b}$ results in the pressure gradient parameter having a maximum value at the end of the simulation, when $Re_{\tau}\approx 180$. The different sweeps will therefore be referred to by their maximum pressure gradient parameter, which are $\varDelta_{p,Re_{\tau}=180}\approx(0.03,0.07,0.15)$. These values place the sweeps considered in this study into the range between negligible ($\varDelta_{p}<0.01$) and noticeable ($-\varDelta_{p}>0.73$) pressure gradient effects, and are substantially larger than the steady value ($-\varDelta_{p}=0.0017$).

3 Results

3.1 Varying streamwise domain length (table 1)

Wall-bounded turbulent flows are often characterised by very long structures on the order of tens of channel half heights (Kim & Adrian, 1999; Lozano-Durán & Jiménez, 2014). This becomes very expensive to simulate for conventional full-span channel simulations, so shorter domain lengths of $L_{x}/h=2\pi$ are often used (Chin et al., 2010; Munters et al., 2016). These seem to capture the majority of the outer-layer dynamics despite their relatively short length. However, these lengths are only necessary due to the large outer-layer structures which are present in the full-span simulations, which have streamwise lengths on the order of the channel half-height, $h$.

The minimal-span channel only captures the inner-layer flow which is not dependent on $h$, suggesting a shorter domain length can be used. The near-wall cycle of the buffer layer produces streaky structures with streamwise lengths of 1000 viscous units (Kline et al., 1967). However, these low- and high-speed streaks are nearly two-dimensional, which suggests they could be represented in the infinite ($k_{x}=0$) mode and hence the domain length doesn’t have to be this long. These streaks are accompanied by quasi-streamwise vortices with streamwise lengths of 200–300 viscous units (Jeong et al., 1997). The narrowest minimal-span channels with $L_{y}^{+}\approx 100$ therefore require $L_{x}^{+}\approx 200$–300 to capture the near-wall cycle in the buffer layer (Jiménez & Moin, 1991). The inertial logarithmic layer also have similar vortex clusters (del Álamo et al., 2006), which scale as 2–3 times their spanwise width (Flores & Jiménez, 2010; Hwang, 2015). Statistically stationary and homogenous shear turbulence, which shows similarities to the logarithmic layer in wall-bounded turbulence, also suggest a scaling of 2 times the spanwise length (Sekimoto et al., 2016). Simulations with larger spanwise domain widths would produce vortices with longer streamwise lengths, so that the largest captured eddy would have spanwise and streamwise lengths of $\lambda_{y}=L_{y}$ and $\lambda_{x}=2$–$3L_{y}$.

Two different spanwise domain widths are investigated for varying the streamwise length (table 1). Firstly, the smallest width of $L_{y}^{+}\approx 120$ would have the largest eddies having streamwise lengths of $\lambda_{x}^{+}=240$–350. Various streamwise domain lengths of $L_{x}^{+}=(190,300,460,930,1850)$ are chosen to see how this affects the largest eddies. The second spanwise domain width of $L_{y}^{+}\approx 350$ would be able to capture larger eddies, the largest of which would have a streamwise length of approximately $\lambda_{x}^{+}=710$–$1100$. The same set of streamwise domain lengths are simulated for this wider domain. A standard-height channel is used, with no outer-layer damping, $K_{i}=0$.

Figure 3 shows the mean velocity profile for the two spanwise widths, for all the streamwise lengths tested. The minimal channels (grey lines) tend to agree with the full-span channel (black line) up until the critical wall-normal location $z_{c}=0.4L_{y}$ (vertical dotted line), at which point the velocity profile of the minimal channels increases. For the smaller domain width (figure 3a,b), the shortest streamwise lengths of $L_{x}^{+}\approx 190$ and $300$ result in a reduction in the smooth-wall velocity above $z^{+}\gtrsim 20$. There is also a slight increase in the centreline velocity. As a result, the difference in smooth- and rough-wall velocity (inset of figure 3b) decreases relative to the full-span and longer minimal-span cases. This effect is not due to the differences in grid resolution of the two shortest channels (table 1), as the wider channel has the same differences in grid resolution but will be seen to not have this effect. Little discernible difference can be seen between the longer streamwise lengths of $L_{x}^{+}\geq 460$, in agreement with the scalings discussed above. This also indicates that the low- and high-speed streaks, with lengths of 1000 viscous units, do not need to be explicitly captured as they are aliased into the $k_{x}=0$ mode. The case of $L_{x}^{+}\approx 300=2.5L_{y}^{+}$ produces an incorrect profile and this suggests that a minimum streamwise length of $L_{x}^{+}=3L_{y}^{+}$ is required, especially for narrow channels in which $z_{c}$ is closer to the buffer layer than log layer of the flow.

The smallest streamwise length of $L_{x}^{+}\approx 190$ was unable to sustain turbulence for a prolonged period. It would decay to a laminar state after $T^{+}=TU_{\tau}^{2}/\nu\approx 14.8\times 10^{4}$ in the smooth-wall channel, and after $T^{+}\approx 9.7\times 10^{4}$ in the rough-wall channel. The data shown in the previous figures for this streamwise length is averaged only when the flow is turbulent. This behaviour is similar to what was observed in Jiménez & Moin (1991), who showed that turbulence could not be maintained in channels with streamwise lengths of around 200 viscous units. As discussed above, this is because it is shorter than the quasi-streamwise vortices which have streamwise lengths of $\lambda_{x}^{+}\approx 200$–$300$.

The wider spanwise domain width of $L_{y}^{+}\approx 354$ (figure 3c,d) results in a centreline velocity that is closer to the full-span channel, as these wider channels capture a wider range of turbulent motions. However, the case with the shortest streamwise length of $L_{x}^{+}\approx 190\approx 0.54L_{y}^{+}$ has a larger centreline velocity, with the critical wall-normal height appearing lower than $z_{c}=0.4L_{y}$. This somewhat resembles the channels with a narrower spanwise width of $L_{y}^{+}\approx 120$ (figure 3a,b), suggesting that the very short domain length restricts the size of the largest eddies. Their spanwise width would now be smaller than the width of the domain, i.e. the streamwise length is now the limiting length scale. Interestingly, the streamwise length of $L_{x}^{+}\approx 300=0.84L_{y}^{+}$ appears to agree quite well with the cases with longer lengths. This is in contrast to the narrow domain (figure 3a), which shows a clear difference for $L_{x}^{+}\approx 300\approx 2.5L_{y}^{+}$. Even the case with $L_{x}^{+}\approx 460\approx 1.3L_{y}^{+}$ is producing a velocity profile and roughness function comparable to that of channels with longer streamwise lengths, despite not having the requisite scaling of $L_{x}^{+}\gtrsim 2$–$3L_{y}^{+}$ discussed above. This is similar to the results of Toh & Itano (2005) who looked at wide spanwise channels with very short streamwise lengths. The velocity profiles they obtained look similar to a full-span, full-length channel, with no apparent increase in streamwise velocity in the outer layer that is characteristic of minimal-span channels. A possible explanation for the results seen here is that the wall-normal critical location $z_{c}^{+}\approx 140=0.24h^{+}$ is outside the log layer, which is generally believed to end at $z\approx 0.15h$ (Marusic et al., 2013). As such, the expected streamwise length scale of 2–3$L_{y}$ is no longer appropriate in the outer layer. In this case, the largest eddy at the edge of the logarithmic layer would have a streamwise length of approximately $\lambda_{x}^{+}\approx$440–660, which could explain why the simulation with $L_{x}^{+}\approx 460$ produces a reasonable velocity profile.

Premultiplied one-dimensional streamwise energy spectra of streamwise velocity, $k_{x}E_{uu}^{+}$, are shown in figure 4, where ${{u}_{rms}^{\prime 2}}=\int_{0}^{\infty}k_{x}E_{uu}\hphantom{.}\mathrm{d}\log(k_{x})$. Only the smooth-wall flow is shown as the rough-wall data exhibits the same qualitative trends as the smooth wall and so is omitted. This occurs because the rough-wall flow is in the transitionally rough regime, which leads to a weakening of the near-wall cycle (MacDonald et al., 2016). As it hasn’t been destroyed, the buffer-layer streaks and quasi-streamwise vortices are active and so the same trends as the smooth-wall flow are observed. Note that in the fully rough regime, this cycle is believed to be destroyed with the inertial logarithmic region started at the roughness crests (Jiménez, 2004). However, the scaling of $\lambda_{x}\approx 2$–$3\lambda_{y}$ would likely still hold as this is generated by a self-sustaining process in the logarithmic layer (Flores & Jiménez, 2010; Hwang, 2015). The vertical dashed lines correspond to the longest streamwise scale based on the log-layer scaling of $3L_{y}^{+}$ discussed above, while the termination of the line contours of the minimal channels show the channel streamwise length. The dashed square in the wider channels (figure 4b,d,f,h) shows the extra captured region due to having a larger $L_{y}^{+}$, compared to the narrower minimal channels (figure 4a,c,e,g).

It is clear that for the narrowest channel with $L_{y}^{+}\approx 120$, the shortest streamwise length of $L_{x}^{+}\approx 190=1.6L_{y}^{+}$ (figure 4a) is too short. This is seen to result in an increase in turbulent energy above $z_{c}^{+}$ (horizontal dashed line), relative to the full-span channel (shaded contour). A similar effect is observed for $L_{x}^{+}\approx 300=2.5L_{y}^{+}$ (not shown). This increase in energy at smaller streamwise length scales above $z_{c}^{+}$ is in agreement with previous minimal-channel simulations (Hwang, 2013). When the streamwise domain length exceeds the $3L_{y}$ scaling with $L_{x}^{+}\approx 460=3.9L_{y}^{+}$ (figure 4c) then reasonable agreement with the full-span channel is observed, despite the narrow range of scales captured by the minimal channel. Further increases to the streamwise length for this domain width (figure 4e,g) do not improve the collapse with the full-span channel, especially near $z_{c}^{+}$, as no additional turbulent motions are captured according to the $3L_{y}^{+}$ scaling discussed above. The increased length is however able to capture the near-wall cycle with peak at $z^{+}\approx 15$ and $\lambda_{x}^{+}\approx 1000$ (figure 4g).

A similar picture emerges for the wider minimal channel with $L_{y}^{+}\approx 354$. The shortest domain length tested, $L_{x}^{+}\approx 190\approx 0.54L_{y}^{+}$ (figure 4b), is unable to capture much of the turbulent motions, similar to that observed in the narrower minimal channel of figure 4(a). A streamwise domain length of $L_{x}^{+}\approx 463\approx 1.3L_{y}^{+}$ (figure 4d) has better agreement with the full-span channel, although there is still some discrepancy. The $3L_{y}$ scaling is almost reached in figure 4(f) with $L_{x}+\approx 927\approx 2.6L_{y}$, and excellent agreement is observed with the full-span channel, with only a slight increase in energy above $z_{c}^{+}$. Further increasing $L_{x}^{+}$ above $3L_{y}$ (figure 4h) does not provide any improvement other than capturing more of the near-wall cycle, further supporting this scaling argument.

The premultiplied spanwise energy spectra of streamwise velocity are shown in figure 5. Here, the vertical dashed line now shows the widest spanwise length scale that can exist based on the streamwise domain length using the $\lambda_{y}=3L_{x}$ scaling. For the narrower channel with $L_{y}^{+}\approx 120$, the scaling is reached when $L_{x}^{+}\approx 460=3.9L_{y}^{+}$ (figure 5c), and further increases to $L_{x}^{+}$ (figure 5e,g) do not result in any change to the turbulent energy. Similarly for the wider channel with $L_{y}^{+}\approx 350$, the scaling is approximately reached in figure 5(f) with $L_{x}^{+}\approx 930=2.6L_{y}^{+}$ and increasing $L_{x}^{+}$ to 1850 viscous units (figure 5h) does not result in any substantial change to the spectra. The cases with $L_{x}\lesssim 2$–$3L_{y}$ (figure 5a,b,d) result in enhanced turbulence energy, particularly in the near-wall peak at $z^{+}\approx 15$ and $\lambda_{y}^{+}\approx 100$.

Given the above results, it is possible to restrict the streamwise length to $L_{x}=3L_{y}$. This is due to quasi-streamwise vortices that exist in the logarithmic layer (where neither viscosity or the channel half height are relevant length scales), suggesting this scaling is independent of $Re_{\tau}$. Smaller streamwise lengths than this scaling lead to poor agreement between the minimal and full-span channel, as seen in figures 4 and 5. However, for very narrow channels this can result in a streamwise domain length less than 1000 viscous units. While it appears the buffer-layer streaks do not need to be captured in the domain, having such a small streamwise domain exacerbates the bursty nature of the minimal-span channel (Jiménez, 2015). Only one or maybe two quasi-streamwise vortices are present in the domain, which makes obtaining converged statistics difficult as the simulation needs to be run for a significantly long time, an issue discussed in detail in §4. We believe that a minimum streamwise length of approximately 1000 viscous units seems a reasonable length in these cases, so that several of the smallest quasi-streamwise vortices are present. As such, we recommend that $L_{x}\gtrsim\max(\hphantom{.}3L_{y},1000\nu/U_{\tau},\lambda_{r,x})$, where the last requirement pertains to the streamwise length scale, $\lambda_{r,x}$, of the roughness which is absent in the current roughness forcing model.

3.2 Half-height channel (table 2)

The previous section looked at the effect of the streamwise domain. Now we will consider the effect of the wall-normal domain, particularly in terms of the outer-layer flow which is unphysical in minimal channels. First, we will consider a half-height (open) channel which consists of a slip wall positioned at the channel centreline. Intuitively, changing the boundary condition at $z=h$ is unlikely to effect the flow at $z=z_{c}$, given the distance between these two heights ($z_{c}/h\approx 0.1$ in this study). For efficient roughness simulations, this has the benefit of reducing the size of the grid by a half when compared to a conventional standard-height channel with two no-slip walls. The simulations are detailed in table 2 and both full-span and minimal-span channels are simulated. Here, the streamwise length is fixed at $L_{x}^{+}\approx 3700$ and there is no outer-layer damping, $K_{i}=0$.

Figure 6 shows the effect of the half-height channel in terms of the mean velocity profile, when compared to the standard-height channel. This is done for both full-span (figure 6a) and minimal-span (figure 6b) channels, for smooth and rough walls. For clarity, the full-span and minimal-span channels are shown in different figures, however the near-wall behaviour is identical up until $z_{c}=0.4L_{y}\Rightarrow z_{c}^{+}\approx 50$, as observed previously (figure 3). Both sets of figures show that the use of the half-height channel with slip wall has a negligible effect on the flow. The main difference is seen in the wake region, where the half-height channel restricts the outer-layer structures, resulting in a slight decrease in the mean velocity. Crucially, this difference is the same for both smooth and rough walls, meaning that the difference between them, $U_{S}^{+}-U_{R}^{+}$, is the same. Moreover, the change in the half-height channel occurs above the critical wall-normal location, where the minimal-span flow is already altered compared to the full-span flow.

The streamwise, spanwise, and wall-normal root-mean-squared velocity fluctuations are shown in figure 7 for both full-span and minimal-span channels, comparing standard-height and half-height channels. For the full-span channel (figure 7a), these velocity fluctuations show very good agreement between the standard-height (black lines) and half-height (grey lines) channels in the near-wall region, in agreement with previous full-span half-height channel studies (Handler et al., 1999; Leonardi & Castro, 2010). The half-height channel has zero wall-normal velocity fluctuations at the channel centre due to the impermeability constraint of the slip wall. The streamwise velocity fluctuations of the half-height channels are slightly enhanced above $z^{+}\gtrsim 40$ relative to the standard-height channels, however the difference is relatively small. Moreover, outer-layer similarity is still maintained between the smooth- and rough-wall channels, suggesting the effect on the difference between the two flows will be minor.

For the minimal-span channel (figure 7b), a slightly different picture emerges. Firstly, it should be noted that a standard-height minimal-span channel (black lines of figure 7b) have enhanced streamwise and wall-normal velocity fluctuations compared to its full-span counterpart (black lines of figure 7a) in the outer layer. This has been noted in other minimal-span channel studies (e.g. Hwang 2013; MacDonald et al. 2016), so will not be discussed in depth here. However, when a half-height channel is used in a minimal-span channel, the streamwise velocity fluctuations (grey lines of figure 7b) are damped from $z^{+}\gtrsim 20$ compared to the standard-height minimal-span channel. This is the opposite of what occurred in a full-span channel. Interestingly, Hwang (2013) found that when the spanwise uniform eddies (that is, the $k_{y}=0$ mode) were filtered out of a minimal-span channel, a similar behaviour was observed in that the wall-normal and streamwise velocity fluctuations were damped compared to the unfiltered channel. It is unclear whether the half-height channel is performing a similar operation to this $k_{y}=0$ filtering, as this would imply the impermeability constraint is limiting the infinite spanwise motions in some way. However, it is clear that the near-wall flow is preserved, despite the imposition of an outer-layer boundary condition. This shows that a half-height channel flow can be simulated without significant near-wall detriment.

The pressure fluctuations, meanwhile, show a substantial increase in the minimal-span channel compared to the full-span channel. The premultiplied streamwise spectra of pressure is shown in figure 8(a), and it is clear that the fluctuations are larger at longer wavelengths ($\lambda_{x}^{+}\approx 1000$–2000) and at heights above $z_{c}^{+}$. The fluctuations at wavelengths shorter than $3L_{y}^{+}$ and below $z_{c}^{+}$ are, however, in reasonable agreement with the full-span channel. Despite this, because the root-mean-squared pressure fluctuations are the integral of this spectrum ($p_{rms}^{\prime 2}=\int_{0}^{\infty}k_{x}E_{pp}d\log(k_{x})$), then the profile is nearly an order of magnitude larger in the minimal-span channel (grey line, figure 8b). Note that the mean pressure is still correct, as this must obey the averaged wall-normal momentum equation, $\overline{w^{2}}+\overline{p}=C$, where $C$ is a constant.

A possible solution to this behaviour is to consider decomposing the pressure into its rapid and slow parts (Kim, 1989). Here, the rapid pressure emerges due to the mean streamwise velocity gradient while the slow pressure is due to the velocity fluctuation gradients. The mean velocity gradient $dU/dz$ in the minimal-span channel is enhanced above $z_{c}$, which means the rapid pressure will also be larger. If we instead define an altered streamwise velocity field whereby $u^{a}(x,y,z)=u(x,y,z)-(U(z)-U(z_{c}))$ for $z>z_{c}$ where $U(z_{c})$ is the mean streamwise velocity at $z=z_{c}$, this will then set the mean $dU^{a}/dz=0$. If the other two velocity components are unchanged and $u^{a}=u$ for $z\leq z_{c}$, then the rapid pressure would be zero above $z_{c}$. The pressure field arising from this altered velocity field, $\nabla^{2}p^{a}=-u^{a}_{i,j}u^{a}_{j,i}$ can then be computed and is seen to be in much better agreement with the full-span channel (figure 8b). There is a slight difference in the near-wall peak near $z^{+}\approx 30$, however this may be due to the discontinuity in $dU^{a}/dz$ at $z=z_{c}$. This correction to pressure is an additional step that needs to be performed whenever pressure statistics are to be outputted. The unaltered (discrete) pressure is still necessary at each time step to ensure the flow remains divergence free.

3.3 Outer-layer damping (table 3)

The previous section showed that a half-height (open) channel did not significantly alter the healthy near-wall flow. A more aggressive alteration is now attempted in which the outer-layer flow is damped through use of the $K_{i}$ forcing term (3) in a half-height channel. This will reduce the large centreline velocity (figure 6b) of the minimal-span channel, which places a restriction on the time step due to the CFL number. The streamwise domain length is fixed at $L_{x}^{+}\approx 3700$, and other relevant parameters are detailed in table 3. The height where the damping begins, $z_{d}$, is varied in channels of two different spanwise widths; $L_{y}^{+}\approx 120$ and $L_{y}^{+}\approx 240$. These minimal channels will have a healthy turbulence region up to $z_{c}^{+}\approx 47$ and $z_{c}^{+}\approx 94$, respectively.

Figure 9 shows the mean streamwise velocity profile for smooth- and rough-wall minimal-span channels, for both of the spanwise widths at $Re_{\tau}=590$. The smallest spanwise width, $L_{y}^{+}\approx 120\Rightarrow z_{c}^{+}\approx 47$ (figure 9a,b) has four different positions for $z_{d}$ examined, with values of $z_{d}^{+}=(100,150,200,250)$. In all cases, the damping tends to limit the streamwise velocity above $z>z_{d}$ to the value of $U(z=z_{d})$. It is clear in the smooth-wall flow (figure 9a) that the two smallest values of $z_{d}^{+}=100$ and 150 are too close to the wall, contaminating the healthy near-wall flow and resulting in an increase in the streamwise velocity below $z_{d}$. In particular, $z_{d}^{+}=100$ results in an obvious overestimation of $\Delta U^{+}$, as seen in the inset of figure 9(b). A value of $z_{d}^{+}=150$ has a similar, although not as strong effect. This is clearly not desirable and suggests that values of $z_{d}^{+}\geq 200$ are necessary to avoid contamination of the healthy near-wall flow.

To investigate how the location of $z_{d}$ is related to $z_{c}$ will require cases with a larger spanwise domain width, as $z_{c}\approx 0.4L_{y}$. This is done in figure 9(c,d), where the spanwise width is now $L_{y}^{+}\approx 240\Rightarrow z_{c}^{+}\approx 94$. Here, three different positions for $z_{d}$ are analysed, with values of $z_{d}^{+}=(200,250,300)$. There is little difference in the near-wall profiles of the three different $z_{d}$ positions. This suggests that, even though $z_{c}$ has doubled, we still require the damping to begin at least 200 viscous-units away from the wall. Having the outer-layer damping begin closer to the wall (i.e. reducing $z_{d}$) leads to its effects permeating the entire near-wall region, increasing the viscous stress in the log-layer and even buffer layer. This results in an overestimation of the mean streamwise velocity. The strong effect of the damping when $z_{d}$ is small could also be due to the strength of the damping, $\gamma h/U_{\tau}\approx 1$ in (3). Reducing the strength would likely reduce the strong near-wall effects that the damping has, although this is not pursued here.

Next, the friction Reynolds number is increased to $Re_{\tau}=2000$. Since the roughness height is fixed as a ratio of the channel half-height, $k=h/40$, then the roughness Reynolds number increases to $k^{+}=50$. The spanwise width also needs to be increased to submerge the roughness sublayer in healthy turbulence (Chung et al., 2015); here it takes a value of $L_{y}^{+}=300$ so that the critical wall-normal location is $z_{c}^{+}\approx 120=2.4k$. Figure 10 shows the mean velocity profile for these higher $Re_{\tau}$ simulations in which $z_{d}^{+}=(200,300)$. Only smooth-wall full-span channel data is available from Hoyas & Jiménez (2006) to compare with the minimal-span channel, and good agreement is seen up until $z^{+}\approx 120=0.4L_{y}^{+}$ (figure 10a). Here, both positions of $z_{d}$ produce velocity profiles which compare well with the minimal-span channel with no outer-layer damping, up to the forcing location $z_{d}$.

It is also of interest to see what effect the outer-layer damping has on second-order turbulence statistics. This is shown in figure 11 for the streamwise and wall-normal velocity fluctuations (figure 11a) and the Reynolds shear stress (figure 11b), for the case of outer forcing with $z_{d}^{+}=200$ and $L_{y}^{+}=120$. This is compared to the minimal-span channel with the same spanwise width, but with no outer-layer damping. The turbulence intensity up until $z^{+}\approx 30$ agrees reasonably well between the case with forcing and without, however because the forcing term scales as $(u_{i}-\langle u_{i}\rangle)^{2}$, fluctuations are almost zero above $z>z_{d}$. The fact that they are not exactly zero is likely because only the spatially averaged velocity $\langle u_{i}\rangle$ is used in the forcing term, and so will vary slightly with time given the small spatial averaging domain of the minimal channel.

The Reynolds shear stress (figure 11b) is almost zero above $z^{+}>300$, meaning the damping is not dissimilar to having a lower friction Reynolds number while keeping the roughness Reynolds number the same. However, flow structures can enter into the damping region and remain coherent for a period of time which likely depends on the damping strength $\gamma$. This is different to a half-height channel at a lower friction Reynolds number, in which the impermeability constraint confines the structures at the slip wall. An advantage of higher $Re_{\tau}$ with damping compared to lower $Re_{\tau}$ without damping is that the higher friction Reynolds number of the flow with damping will reduce the effect of low friction Reynolds number flows. Chan et al. (2015) showed that the roughness function $\Delta U^{+}$ was overestimated for $Re_{\tau}\approx 180$ compared to $Re_{\tau}\approx 360$ and $Re_{\tau}\approx 540$ for matched roughness Reynolds numbers. This is due to there being a non-negligible pressure gradient effect, which as discussed in §2.1 is $\varDelta_{p}=-1/Re_{\tau}$. It would therefore be preferable to conduct higher friction Reynolds number simulations with damping, rather than a lower friction Reynolds number simulations. While not investigated here, it also seems reasonable that the grid spacings above $z_{d}$ could be relaxed without altering the near-wall flow.

This damping is similar to the masks, or sponge regions, employed in Jiménez & Pinelli (1999). The authors explicitly filtered out disturbances in the outer layer, resulting in laminar-like flow in the outer layer but still maintaining a healthy near-wall cycle. The present damping is similar, although some fluctuations are still present in the outer layer (figure 11a) and the velocity is fixed rather than forming a laminar velocity profile. A benefit of this damping is that it reduces the computational time step, which is limited by the CFL number, $\textit{CFL}=\max_{1\leq i\leq 3}\left(|u_{i}|\Delta t/\Delta x_{i}\right)$. Here $|u_{i}|$ is the maximum instantaneous velocity in the $i$th direction and $\Delta x_{i}$ is the corresponding grid spacing at that location. Note that the viscous time step restriction does not become significant in any of the simulations performed here. For minimal-span channels at $Re_{\tau}=590$ with no outer-layer damping, the maximum CFL number occurs in the outer layer due to the large streamwise velocity. We see the time step improves by approximately 20–24% with the use of the outer-layer damping (table 3, $L_{y}^{+}\approx 118$) than without. This means the simulations can be completed quicker and hence use less computational resources. However, for the present simulations at higher friction Reynolds numbers of $Re_{\tau}=2000$, the time step for the smooth wall with outer-layer damping is unaltered compared to the case without damping. This is due to the use of the cosine mapping to define the wall-normal grid, which produces an excessively fine grid near the wall especially for high Reynolds number cases. As a result of such a small $\Delta z$, the wall-normal CFL number is now maximum and this occurs in the near-wall region, which is independent of the outer-layer damping. This clustering of points near the wall is a known issue of the cosine mapping for high $Re_{\tau}$ simulations (Lenaers et al., 2014), but is used here to ensure a systematic, controlled study. If a more appropriate mapping was used for these high friction Reynolds numbers (e.g. that in Lee & Moser 2015), we would expect a similar or better improvement in time step as at $Re_{\tau}=590$.

The most appropriate position for $z_{d}$ is difficult to establish from the above data. A minimum of $z_{d}^{+}\approx 200$ seems necessary to ensure that the damping does not interact with the near-wall flow. It also seems reasonable that $z_{d}$ should scale in some way on $z_{c}$, for larger domain widths. The three sets of simulations suggest that $z_{d}$ should be approximately two times $z_{c}$, so a tentative rule of thumb would be $z_{d}^{+}\gtrsim\max(200,\hphantom{.}2z_{c}^{+})$. For an appropriate wall-normal mapping, this allows a 20–24% improvement in the computational time step.

3.4 Temporal sweep (table 4)

In this section, the temporal sweep (§2.1) is performed in a standard-height channel, with no outer-layer damping ($K_{i}=0$) and a streamwise length of $L_{x}/h=2\pi$ (table 4). The initial friction Reynolds number is $Re_{\tau}=590$ which reduces to $Re_{\tau}=180$ via an adverse pressure gradient with different rates investigated. The roughness height is fixed at $k=h/40$, so the roughness Reynolds number varies as $5\lesssim k^{+}\lesssim 15$. Three different sweeps were conducted in which the pressure gradient parameter at the end of the sweep is $\varDelta_{p,Re_{\tau}=180}=(0.03,0.07,0.15)$.

Figure 12 shows the mean velocity profiles for the fastest and slowest sweeps, comparing the sweep data with full-span steady-flow data at the highest and lowest friction Reynolds numbers tested. Figure 12(a,b) shows the mean velocity for the slowest sweep, when the pressure gradient parameter (6) at $Re_{\tau}\approx 180$ was $\varDelta_{p,Re_{\tau}=180}=0.03$. Figure 12(c,d) shows the mean velocity for the fastest sweep, with $\varDelta_{p,Re_{\tau}=180}=0.15$. To obtain statistics for a desired $Re_{\tau}$, time-averaging is performed over a window of $\pm 10\%$ of that value. At the initial friction Reynolds number of the simulation, $Re_{\tau}=590$, both the sweeps (figure 12a,c) show good agreement with the full-span channel in the near-wall region. As the spanwise width is $L_{y}/h=0.6$, then the critical wall-normal location is approximately $z_{c}^{+}\approx 0.4L_{y}^{+}=142$ which agrees well with the figure. As $\varDelta_{p}$ is not significant yet (figure 2b) then there is little difference between the two sweeps.

Differences start to emerge at the end of the sweep when the friction Reynolds number is close to $Re_{\tau}=180$ and the $\varDelta_{p}$ is maximum. The slowest sweep (figure 12b) has reasonable agreement very close to the wall below $z^{+}\lesssim 15$, but the pressure gradient has reduced the streamwise velocity compared to the full-span steady flow. However, the difference between the smooth and rough-wall sweep flows (inset of figure 12b) agrees quite well with the full-span steady data. Taking the difference between two flows with the same weak pressure gradient ostensibly still produces a correct estimate of $\Delta U^{+}$. However, when the pressure gradient is increased to to $\varDelta_{p,Re_{\tau}=180}=0.15$ (figure 12d) then larger differences are seen. The near-wall region has been noticeably changed, and the difference in smooth and rough wall velocity (inset) is seen to overestimate the full-span data. Note that the data have been ensemble averaged over four runs, so that the amount of time the flow is averaged over is the same as in the slowest sweep. This therefore indicates that the overestimation of $\Delta U^{+}$ is a direct result of the stronger (adverse) pressure gradient.

The roughness function $\Delta U^{+}$ is now computed from the difference in smooth and rough-wall flows at matched $U_{\tau}$ values. As seen in the insets of figure 12, the smooth- and rough-wall velocity difference is not constant with $z^{+}$ as in the steady data. This indicates that the flow has not been averaged for long enough at the desired $Re_{\tau}$ value. As such, multiple sweeps from different initial conditions would need to be conducted to obtain a converged outer region. However, for the purposes of this investigatory study, the velocity difference is averaged from the crest of the roughness to the channel centre to obtain an estimate of $\Delta U^{+}$. This is done for a range of friction Reynolds numbers, with the resulting data plotted in figure 13 against the equivalent sandgrain roughness, $k_{s}^{+}$. Ideally, the temporal sweep data (lines) should agree with the steady flow data (black symbols). Indeed, the temporal sweep does show promise in this regard, especially for the weaker pressure gradient cases. If we interpolate the steady roughness function data from Chung et al. (2015), then an overall relative error term can be approximated as $\varepsilon=\int|1-\Delta U^{+}_{steady}/\Delta U^{+}_{sweep}|\hphantom{.}\mathrm{d}k_{s}^{+}$. The two slowest sweeps, with $\varDelta_{p}=0.03$ and 0.07 have similar relative errors of $\varepsilon\approx 12\%$ and 13%, respectively, while the fastest sweep is noticeably larger with an error of $\varepsilon\approx 26\%$.

It appears that the sweep approach could work for determining the hydraulic behaviour of roughness. A pressure gradient parameter of $\varDelta_{p}\leq 0.03$–0.07 is required, which is not too far from the recommendation in Patel (1965) that zero pressure gradient conditions can be assumed when $\varDelta_{p}\leq 0.01$. A value of $\varDelta_{p}\gtrsim 0.15$ results in pressure gradient effects that lead to inaccuracies in $\Delta U^{+}$. Future studies could prescribe a constant value of $\varDelta_{p}$ rather than linearly varying the bulk velocity, as this would ensure pressure gradient effects remain negligible, while efficiently traversing the range of roughness Reynolds numbers.

4 Cost and Convergence

Quantifying the cost of the simulations is extremely important when it comes to evaluating the benefit of minimal-span channels. There are two important costs to consider in these simulations; the cost of memory (or the number of processors required) and the cost of simulation time. The memory cost can be readily described by the size of the grid: the number of spatial degrees of freedom for full-span simulations scale as $Re_{\tau}^{9/4}$ (Pope, 2000, §9.1.2) which for $Re_{\tau}\gtrsim 10^{3}$ requires national-level high-performance computing facilities. The minimal channel, meanwhile, scales as $k_{s}^{+9/4}$ (Chung et al., 2015), making it feasible for smaller facilities at the university and industry level, or potentially high-end desktop computers.

The cost of the simulation runtime is conventionally given in terms of CPU hours, however an obvious question emerges; does a minimal channel with a domain volume that is, say, 10 times smaller than a full-span channel require the simulation runtime to be 10 times longer? This would imply a similar cost in terms of CPU hours, but to answer this question definitively we need to know how long it takes to obtain converged statistics, particularly in the first-order statistics necessary to obtain the roughness function. A benefit of this analysis is that it will enable us to predict how many CPU hours are required for a minimal-span simulation before actually performing the simulation.

Generally, full-scale turbulence simulations are run for approximately 10 large-eddy turnover times, where the largest eddy captured in the domain are of characteristic size $h$ and velocity $U_{\tau}$. This implies a simulation run time of $T_{sim}U_{\tau}/h\approx 10$ (see, for example, Hoyas & Jiménez 2006; Lozano-Durán & Jiménez 2014; Lee & Moser 2015) and this long runtime is why DNS of full-span channels are so expensive. However, the largest captured eddies in the minimal-span channel will be of characteristic size $z_{c}$, implying we need $T_{sim}U_{\tau}/z_{c}\approx 10$ eddy turnover times of this $z_{c}$-sized eddy. In this study, these are generally of the size $h/z_{c}\approx 10$, implying that $T_{sim}U_{\tau}/z_{c}\approx 10$ would be only 1 turnover time of the $h$-sized eddy (if it existed), representing an order of magnitude saving in time.

This argument, however, does not make reference to the wall-parallel size of the domain. Lozano-Durán & Jiménez (2014) suggested that channels with smaller domain sizes need to be run for proportionately longer, and showed that the statistical uncertainties in $u^{\prime}(z=h$) of full-span channels decrease in a manner inversely proportional to the square root of the wall area. However, this result was for two channels that, from the point of view of the current study, were full-span channels. The largest characteristic eddies in both channels were of size $h$, which explains the simple relationship between channel domain size and simulation runtime. In contrast, the largest characteristic eddies in minimal channels are of size $z_{c}$ which varies depending on the width of the channel. This means that such a simple relationship between channel domain size and run time is unlikely to hold. In the context of these $z_{c}$-sized eddies, a full-span simulation could potentially have hundreds of eddies distributed throughout the domain, while a minimal-span simulation will only have a few. For this minimal-span channel to obtain the same level of converged statistics as a full-span channel would presumably require the same number of these $z_{c}$-sized eddies to pass through the domain. It therefore becomes useful to count the number of $z_{c}$-sized eddies present in the domain. These eddies have a spanwise width of $\lambda_{y,z_{c}}=L_{y}\approx z_{c}/0.4=2.5z_{c}$ and, as discussed in the previous section, a streamwise length of approximately three times the span, $\lambda_{x,z_{c}}\approx 3\lambda_{y,z_{c}}=7.5z_{c}$. Hence, the number of $z_{c}$-sized eddies at any instant is

[TABLE]

The first factor is an approximation of the number of $z_{c}$-sized eddies present on one wall of the channel. The other factor, $L_{z}/h$, is present as a half-height channel ($L_{z}/h=1$) will have half as many $z_{c}$-sized eddies as a standard-height channel ($L_{z}/h=2$).

Over the course of the simulation, these $z_{c}$-sized eddies will grow and decay in a process termed ‘bursting’, so it becomes necessary to consider the time scale of these eddies. The bursting process in the buffer layer was investigated by Jiménez et al. (2005), who showed that the time scale depended on the friction Reynolds number. Low Reynolds number flows have a long bursting period of $T_{b}^{+}\approx 1000$ for $Re_{\tau}\approx 100$, however this saturates to $T_{b}^{+}\approx 400$ for $Re_{\tau}\gtrsim 1000$. This was determined for minimal channels with small spanwise widths such that they only captured the buffer layer, i.e. $L_{y}^{+}\lesssim 125$–$200\Rightarrow z_{c}^{+}\lesssim 50$–$80$. Figure 14(a) shows this buffer-layer bursting period, along with the current data for buffer-layer minimal channels ($L_{y}^{+}\lesssim 125$–$200$) at $Re_{\tau}=590$. Good agreement is seen, with the bursting period at this friction Reynolds number being approximately $T_{b}^{+}\approx 540$. The bursting period was determined in the same manner as Jiménez et al. (2005), defined as the weighted logarithmic average of the frequency spectra of the wall-shear stress.

In the logarithmic layer, Flores & Jiménez (2010) determined that the bursting period scales with distance from the wall according to $T_{b}U_{\tau}/z_{c}=6$, which was supported by Hwang & Bengana (2016). The bursting period of the largest eddies in the minimal channel, $T_{b}^{+}$, therefore depends on both $z_{c}$ and $Re_{\tau}$. Narrow minimal channels which only capture the buffer layer ($z_{c}^{+}\lesssim 50$–$80$) have an $Re_{\tau}$-dependence on the bursting period (figure 14a), while wider minimal channels which capture part of the logarithmic layer ($z_{c}^{+}\gtrsim 50$–$80$) have a bursting period that depends on $z_{c}$. This is summarised in figure 14(b), which also shows the present data. The behaviour for intermediate $z_{c}^{+}$ values between buffer and logarithmic flows (grey band) is unknown and likely depends on both $z_{c}^{+}$ and $Re_{\tau}$. Here, the inertial motions of the logarithmic layer are scarce and are likely competing with the $Re_{\tau}$-dependent motions of the buffer layer, with neither completely dominating. This distinction becomes irrelevant for higher $Re_{\tau}$ channels as the buffer-layer bursting period saturates to $T_{b}^{+}\approx 400$.

The product of the bursting timescale, $T_{b}$, and the instantaneous number of $z_{c}$-sized eddies, $N_{\text{instant}}$, gives the total number of $z_{c}$-sized eddies that exist over the simulation runtime $T_{sim}$,

[TABLE]

$C^{\star}$ is simply a running count of how many $z_{c}$-sized eddies have been sampled, and increases with simulation runtime. To investigate whether this is the correct measure for statistical convergence, the statistical uncertainty of streamwise velocity in different minimal channels is analysed. There are various ways of quantifying the statistical uncertainty; here we will estimate it using the standard error of the mean. Instantaneous snapshots of $u(x,y,z,t_{i})$ at time $t_{i}\leq T_{sim}$ are spatially averaged in the wall-parallel plane to get $\langle u(z,t_{i})\rangle$. The temporal average of these $N$ snapshots gives the overall mean,

[TABLE]

The standard error of this mean for statistically independent samples can be computed as $\epsilon=s/\sqrt{N}$, where $s$ is the unbiased standard deviation of the signal $\langle u(z,t_{i})\rangle$. The quantity $\epsilon$ has the same units as $u$, so can be non-dimensionalised to $\epsilon^{+}$. A 95% confidence interval can then be given in which we are 95% certain the true mean falls within $U^{+}\pm 1.96\epsilon^{+}$ (Benedict & Gould, 1996). The roughness function, being the difference in two estimated velocities, therefore has a 95% probability of falling within the range $\Delta U^{+}\pm 1.96\sqrt{\epsilon_{s}^{+2}+\epsilon_{r}^{+2}}$, where subscripts $s$ and $r$ refers to the uncertainties in the smooth-wall and rough-wall velocities.

In turbulence simulations the samples may not necessarily be statistically independent, so the assumption that $\epsilon=s/\sqrt{N}$ would not hold. Following Trenberth (1984) and Oliver et al. (2014), an effective number of statistically independent samples $N_{eff}=N/T_{0}$ can be defined, where

[TABLE]

is the decorrelation separation distance, ${\rho}$ being the autocorrelation function of $\langle u(z,t)\rangle$ normalised on its variance. The standard error can then be estimated as $\epsilon=s/\sqrt{N_{eff}}$, where

[TABLE]

If the samples were truly independent, then the autocorrelation $\rho$ should be near zero, resulting in $T_{0}=1\Rightarrow N_{eff}=N$ and the familiar unbiased standard deviation for $s$ is obtained (with divisor $N-1$) . The autocorrelation function can be noisy, especially for small $N$, prompting Oliver et al. (2014) to develop an open-source code to fit an autoregressive time series model to a given signal. However we have found that for minimal channel simulations $N$ is sufficiently large such that the autocorrelation definition appears to be adequate, with $T_{0}$ usually being $O(1)$. The coarse-graining approach for estimating the standard error for correlated samples, detailed in the appendix of Hoyas & Jiménez (2008), was also attempted although no appreciable differences were detected, again because the samples are nearly independent.

Increasing the simulation runtime $T_{sim}$ reduces the statistical uncertainties as evidenced in figure 15(a), which shows the standard error of the velocity at a representative wall location $z^{+}=47$. The very short minimal-span channels (light grey lines) only have one or two $z_{c}$-sized eddies present, and the bursting nature of only a few eddies greatly increases the statistical uncertainties, requiring a very long runtime. The dashed line shows $(T_{sim}U_{\tau}/h)^{-1/2}$, indicating that the uncertainties decrease inversely proportional to the square root of time. This is essentially the same result as Hoyas & Jiménez (2008) and Lozano-Durán & Jiménez (2014) who argued the decrease was inversely proportional to the square-root of wall area, the difference here being that we have considered the snapshots as a sampling of time rather than multiples of wall area.

When the standard error is instead considered as a function of $C^{\star}$ (figure 15b) we obtain a better collapse, indicating that counting $z_{c}$-sized eddies as in the $C^{\star}$ formulation is indeed an appropriate measure to use. The data show a collapse with $\epsilon^{+}=K(C^{\star})^{-1/2}$, where here it was determined $K=0.75$ for $z_{c}^{+}\approx 47$. This relationship shows that obtaining a desired error tolerance requires capturing a certain number of $z_{c}$-sized eddies ($C^{\star}$). Increasing the streamwise domain length does not affect this measure, other than to simulate more eddies at once.

Figure 16 shows that there is not a universal number of $z_{c}$-sized eddies that need to be simulated to obtain a desired $\epsilon^{+}$, that is, the coefficient $K$ in $\epsilon^{+}=K(C^{\star})^{-1/2}$ is a function of $z_{c}^{+}$. Figure 16(a) shows that for $z_{c}^{+}\approx 94$, $K=0.37$, while figure 16(b) with $z_{c}^{+}\approx 142$ has $K=0.23$. Figure 16(d) shows that this coefficient decreases with $z_{c}^{+}$, which implies that a wider minimal channel will have a smaller statistical uncertainty than a narrower one, for the same $C^{\star}$ value. This is possibly due to the way in which eddies aggregate in wider channels. As the channel widens and $z_{c}^{+}$ increases, the largest $z_{c}$-sized eddy increases, however there remains a hierarchy of smaller eddies below this largest eddy. The uncertainty $\epsilon^{+}$ of the mean velocity would have contributions from all of these eddies and so would depend on $z_{c}^{+}$. The reason this effect is not already captured in $C^{\star}$ is that this quantity only considers the largest $z_{c}$-sized eddy, whereas $\epsilon^{+}$ ostensibly considers all eddies up to and including those of $z_{c}$ size. The data in figure 16(d) show what appears to be a $-1$ power law, in which case a fit gives $K(z_{c})=33/z_{c}^{+}$.

Recall that the 95% confidence interval of the roughness function is $\Delta U^{+}\pm 1.96\sqrt{\epsilon_{s}^{+2}+\epsilon_{r}^{+2}}$. If both smooth- and rough-wall simulations were to have the same uncertainty, then this interval can be rewritten as

[TABLE]

The advantage of this formulation is that setting a desired error tolerance in $\Delta U^{+}$ enables the simulation run time (and hence CPU hours) to be determined a priori. To see this, the spanwise channel width, $L_{y}\gtrsim\max(100\nu/U_{\tau},k/0.4,\lambda_{r,y})$ (Chung et al., 2015), and streamwise channel length, $L_{x}\gtrsim\max(3L_{y},1000\nu/U_{\tau},\lambda_{r,x})$ (§3.1), must first be determined, which sets $z_{c}=0.4L_{y}$. The user must then determine a desired error tolerance, $\zeta$, for $\Delta U^{+}\pm\zeta$. This error tolerance will depend on the application; riblets have a roughness function $-\Delta U^{+}\lesssim 1$ (Spalart & McLean, 2011), which would require a smaller level of statistical uncertainty, than say an engineering roughness problem where $\Delta U^{+}\approx 10$. Once set however, this can be related to $C^{\star}$ through $\Delta U^{+}\pm 91.4(C^{\star})^{-1/2}/z_{c}^{+}$. From here, (8) can be solved for the simulation time $T_{sim}$ using the minimal channel dimensions and bursting period $T_{b}$. The number of CPU hours is then

[TABLE]

where $\Delta t^{+}$ is the computational time step, so $T^{+}/\Delta t^{+}$ is the number of time steps performed during the simulation. $\mathcal{T}_{1}$ is the average wall clock time taken to perform a single step and $N_{CPU}$ is the number of processors employed in the computation. The time step $\Delta t^{+}$ can be estimated if the streamwise CFL number at the channel centre is unity, i.e. $U_{h}\Delta t/\Delta x=CFL=1$, which gives $\Delta t^{+}=\Delta tU_{\tau}^{2}/\nu=\Delta x^{+}U_{\tau}/U_{h}$. If $\Delta x^{+}\approx 10$ and, for a full-span channel, the centreline velocity is $U_{h}/U_{\tau}\approx 25$ then we have $\Delta t^{+}\approx 0.4$, which agrees with that recommended in Choi & Moin (1994). Body-fitted roughness grids may have a smaller time step due to having a finer grid to resolve the roughness geometry. $\mathcal{T}_{1}$ and $N_{CPU}$ will depend on the code being used. The full-span channel used in this study, run for $TU_{\tau}/h=10$ large-eddy turnover times with $\Delta t^{+}\approx 0.18$, $\mathcal{T}_{1}\approx 26.4$ seconds and $N_{cpu}=64$ requires approximately 15800 CPU hours. For a minimal channel with $h/z_{c}\approx 10$ and 95% confidence interval of $\Delta U^{+}\pm 0.1$, the above analysis leads to a requirement of only 1420 CPU hours (for $N_{CPU}=64$ and $\mathcal{T}_{1}=2.4$ seconds, based on case MS6 of table 2), which is over ten times less than the full-span channel.

5 Minimal channel applied to pyramids

5.1 Computational setup

Following the insights discussed in this paper, a finite-volume code is used to determine the roughness function of square-based pyramids. This is the same code as in Chan et al. (2015) and Chung et al. (2015), based on the work of Mahesh et al. (2004) and Ham & Iaccarino (2004). No outer-layer damping will be used ($\mathbf{K}=\mathbf{0}$) and since a body-fitted grid is used, no roughness forcing is required ($\mathbf{F}=\mathbf{0}$). Three pyramids are simulated with different heights, but all with nominal $Re_{\tau}=840$ and angles between the pyramid edge and the horizontal of $\alpha=22^{\circ}$. Each pyramid will be referred to by its nominal trough-to-peak height, $k_{t}^{+}$ and wavelength $\lambda^{+}$, as a unique identifier given by the notation $k_{t}^{+}$  $\lambda^{+}$. Hence the case 40_198 has a height $k_{t}^{+}\approx 40$ and wavelength $\lambda^{+}\approx 198$. The other two cases are 63_313 and 80_396. These three cases were designed to match the pyramids studied in Schultz & Flack (2009), herein referred to as SF09. The angle of SF09 was reported as the edge angle, however this should be the angle between the pyramid face and the horizontal (M. Schultz, personal communication). The correct edge angles of $16^{\circ}$ are given in table 5. This makes the present pyramids steeper than SF09, however their study included pyramids with $\alpha=35^{\circ}$ and observed no appreciable difference in $\Delta U^{+}$, suggesting this discrepancy is insignificant. The smallest pyramid also has a similar $k_{t}^{+}$ to Di Cicca & Onorato (2016), herein referred to as DCO16, who had $k_{t}^{+}=39$. Their pyramid slope angle was $\alpha=28^{\circ}$ and their pyramids were separated from each other by a small gap of 32 viscous units. However, given that the roughness height is matched, these difference shouldn’t result in a substantially different roughness function. Additionally, the case 63_313 matches the roughness height of Hong et al. (2011) (herein HKS11). The cases that are simulated here, and those available in the literature, are summarised in table 5. The expected full-span skin-friction coefficient, $C_{f,mod}=2U_{\tau}^{2}/U_{b,full}^{2}$, is also given in this table. Due to the altered outer-layer velocity profile of the minimal-span channel, the composite profile of Nagib & Chauhan (2008) is fitted to estimate the full-span velocity profile for $z>z_{c}$. This therefore allows for an estimate of the full-span bulk velocity $U_{b,full}$ to be obtained. The estimated smooth-wall coefficient value of $C_{f,mod}\approx 0.005$ is in good agreement with the empirical fit of Dean (1978), in which $C_{f}=0.073(Re_{b})^{-1/4}=0.0054$ where $Re_{b}=2hU_{b}/\nu$ is the bulk Reynolds number. This shows that the minimal-span channel framework can still be used to estimate the full-span skin-friction coefficient, which is commonly used by engineers.

Now that the dimensions of the roughness are known, the channel domain size can be determined. Following the guidelines of Chung et al. (2015), the spanwise domain width must satisfy $L_{y}\gtrsim\max(100\nu/U_{\tau},k_{t}/0.4,\lambda_{r,y})$. For the smallest pyramids (case 40_198), we have $L_{y}^{+}\gtrsim\max(100,100,198)$ so it is the large wavelength which necessitates $L_{y}=\lambda$. An additional case with the same channel dimensions as 40_198 is simulated (case 40_198_f), but with a finer computational grid to ensure that all the turbulent scales are resolved. This spanwise width of $L_{y}=\lambda$ means that the critical wall-normal location $z_{c}=0.4L_{y}$ for these pyramids is set at $z_{c}^{+}=80\approx 2k_{t}^{+}$. The roughness sublayer may be larger than $2k_{t}$, so a channel with a larger width of $L_{y}=2\lambda$ is also simulated (case 40_198_2), so that the critical location is now $z_{c}=4k_{t}$. For the larger pyramids (cases 63_313 and 80_396), the pyramid wavelength is again the limiting constraint, so $L_{y}=\lambda$. The streamwise domain length was investigated in §3.1, where the recommendation was $L_{x}\gtrsim\max(3L_{y},1000\nu/U_{\tau},\lambda_{r,x})$. For the smaller pyramids (40_198), we have $L_{x}^{+}\gtrsim(594,1000,198)$ so it is the second constraint that is the limiting one. To ensure we have complete a number of pyramids in the domain, a streamwise length of $L_{x}=6\lambda$ is selected which in viscous units is $L_{x}^{+}=1188$. Case 63_313 is also limited by the second constraint. For the largest pyramids ($L_{y}^{+}=396$), the first constraint ($L_{x}\gtrsim 3L_{y}$) requires a streamwise length of $L_{x}=3L_{y}\Rightarrow L_{x}^{+}=1188$.

The length of time to achieve a converged result can now be estimated since the spanwise width $L_{y}$ and hence critical wall-normal location $z_{c}=0.4L_{y}$ have been set. For the smallest pyramid case (40_198), the spanwise width is $L_{y}^{+}=198$ so $z_{c}^{+}=79$. If we wish to be 95% confident that the roughness function is in the range $\Delta U^{+}\pm 0.1$, then from (12) the number of $z_{c}$-sized eddies that must be observed are $C^{\star}=(91.4/(0.1z_{c}^{+}))^{2}\approx 134$. We want to use this to determine the simulation runtime $T_{sim}$ through (8), which requires estimating the bursting period $T_{b}$. Since the channel is relatively wide and $z_{c}^{+}=79$, it’s likely to have the beginnings of a logarithmic layer (figure 14b), which implies $T_{b}U_{\tau}/z_{c}=6\Rightarrow T_{b}^{+}=6z_{c}^{+}\approx 475$. All that remains is to substitute this value and the channel dimensions into (8) to determine that $T_{sim}^{+}\approx 3.1\times 10^{4}$. As the nominal streamwise grid spacing is known to be $\Delta x^{+}\approx 4.2$, and if we assume a centreline velocity of $U_{h}^{+}\approx 30$, then the time step $\Delta t^{+}=\Delta x^{+}/U_{h}^{+}\approx 0.14$ for a CFL number of unity at the centreline. Finally, using (13) with an assumed wall clock time per step of $\mathcal{T}_{1}=6.4$ seconds and $N_{CPU}=128$, we expect to use around 52,000 CPU hours. A similar computation for cases 63_313 ($z_{c}^{+}\approx 125$, $C^{\star}\approx 53$) and 80_395 ($z_{c}^{+}\approx 158$, $C^{\star}\approx 33$) predicts 53,000 and 100,000 CPU hours, respectively. For comparison, a single full-span channel ($L_{x}=2\pi$, $L_{y}=\pi$, $L_{z}=2h$ and $T_{sim}U_{\tau}/h=10$) with pyramids would likely require at least $2\times 10^{6}$ CPU hours using the current code, more than a full order of magnitude greater than the minimal channel.

The pyramids are aligned so that the trough between neighbouring pyramids is 45∘ to the mean flow direction. However, because hexahedral cells are used which are aligned with the mean flow direction, the trough line falls across opposite vertices of the cell. This would mean the cell needs to be split into two pentahedrals to fit the roughness geometry. To mitigate this issue, the trough of the pyramids are rounded off and faceted. The resulting trough has a radius of curvature of approximately half the pyramid crest-to-peak height, $R=0.53k_{t}$. The crest of the pyramid retains its sharp edge, as evidenced in figure 17. The total volume of fluid in the channel with pyramids matches that of a smooth-wall channel, which would ensure a collapse of flow statistics in the outer layer, if it existed (Chan et al., 2015). In other words, the hydraulic half height (the channel equivalent to the hydraulic radius of a pipe) matches the smooth wall, which implies the virtual origin is at one third the height of the pyramid.

5.2 Pyramid results

Figure 18(a–c) shows the mean velocity profile for the pyramid roughness data. The insets show the difference in velocity between smooth-wall and pyramid-roughness channels, and it can be seen to be relatively constant from the pyramid crest to the channel centre. This indicates that the velocity profile for the smooth-wall and rough-wall are matched from the crest (apart from the offset $\Delta U^{+}$), or that the roughness sublayer for pyramids is small. The roughness function is obtained by averaging the difference in smooth- and rough-wall velocities over $k_{t}^{+}+5<z<z_{c}^{+}$.

For the smaller pyramids ($k_{t}^{+}=40$) in a narrow channel, we obtained a roughness function of $\Delta U^{+}=8.0$ which agrees well with the refined mesh ($\Delta U^{+}=7.9$). This indicates that both meshes are resolving all the relevant scales in the roughness sublayer. When the channel is widened to include two repeating pyramids in the spanwise direction (40_198_2), we achieve a similar result of $\Delta U^{+}=8.3$. While slightly larger than the previous two values, this difference is predominantly due to the increase in $Re_{\tau}$ for this channel. The roughness Reynolds number, $k_{t}^{+}\approx 41.8$, is almost 5% over the target value of $k_{t}^{+}=40$, which means the roughness function would also be larger. The otherwise good agreement between this case and the narrower case 40_198 supports the conclusion that the roughness sublayer is less than $2k_{t}$ in height above the pyramid crest. If the roughness sublayer was larger than $2k_{t}$, then the narrower minimal channel with $z_{c}=2k_{t}$ could not fully capture the roughness effects and hence would have a different roughness function to the wider minimal channel. Moreover, this agreement with case 40_198, when $L_{y}=\lambda_{y}$, suggests that the interaction of flow structures generated between neighbouring roughness elements is small, as this case only has a single repeating element in the spanwise direction. This agrees with our previous study (MacDonald et al., 2016), which compared a three-dimensional sinusoidal roughness in a full-span channel and a minimal-span channel in which there is only one repeating roughness element in the spanwise direction. It would only be when the spanwise roughness wavelength was very large, or infinite, as in the case of two-dimensional spanwise bar roughness and $d$-type surfaces that these spanwise effects would become significant. However, this would require further study and is outside the scope of this paper. The above three values of the roughness function can be compared with values of $6.4\pm 10\%$ from SF09 and $8.9$ from DCO16. The current data therefore falls between these two sources (figure 18d).

For the case with $k_{t}^{+}\approx 63$, the roughness function was found to be $\Delta U^{+}=10.2$, which is larger than $7.0\pm 10\%$ from SF09 and 8.2 from HKS11. For the largest pyramids, where $k_{t}^{+}\approx 80$, the roughness function was found to be $\Delta U^{+}=11.0$, which is again larger than the value of $8.7\pm 10\%$ from SF09. As a result, the equivalent sandgrain roughness appears to be $k_{s}\approx 3.5k_{t}$, compared to $k_{s}\approx 1.5k_{t}$ from SF09. Note that DCO11 would predict an even larger equivalent sandgrain roughness than $k_{s}\approx 3.5k_{t}$, given their larger $\Delta U^{+}$ at $k_{t}^{+}\approx 40$.

To show some of the higher order statistics that can be obtained from these simulations, the difference between sweep and ejection Reynolds stress contributions are shown in figure 19. This difference, $\Delta S_{0}=(\{u^{\prime}w^{\prime}\}_{4,0}-\{u^{\prime}w^{\prime}\}_{2,0})/\langle\overline{u^{\prime}w^{\prime}}\rangle$, was developed in Raupach (1981) in which the Reynolds shear stress is conditionally averaged on sweeps ($u^{\prime}>0$, $w^{\prime}<0$, quadrant 4) and ejections ($u^{\prime}<0$, $w^{\prime}>0$, quadrant 2), where $\{\cdot\}$ denotes the conditional average. The hyperbolic hole region $H=0$. From this statistic, it can be seen in figure 19(a) that sweeps are the dominant means of momentum transfer within the roughness canopy, in agreement with Raupach (1981). These sweeps are much more dominant than in a smooth wall, which has a slight preference for sweeps very close to the the wall ($z^{+}\lesssim 12$), with ejections dominant above. Above the roughness crest, both rough-wall and smooth-wall flows collapse with ejections being slightly larger than sweeps. The outer-layer is associated with ejections being dominant, with $\Delta S_{0}$ tending to $-1$ at $z=h$. However, in the minimal-span channel this statistic stays approximately constant above $z>z_{c}$, with $\Delta S_{0}$ being close to zero. If instead the wall-normal coordinate is normalised on $k_{t}$ (figure 19b) then we see that all three rough-wall flows nearly collapse, although the case with the smallest height ($k_{t}^{+}\approx 40$, red dashed line) there is a slight difference, with sweeps remaining dominant at the roughness crest.

Figure 20 shows the statistical uncertainty of $U^{+}(z=z_{c})$ as a function of $C^{\star}$, the number of $z_{c}$-sized eddies. This is the same format as figure 16(a–c). Importantly, the predicted number of $z_{c}$-sized eddies that need to be observed to obtain an error tolerance of $\Delta U^{+}\pm 0.1$ using (12) is indicated by the vertical dashed lines. The current data for both smooth and rough walls show reasonable agreement with this prediction, albeit the prediction tending to overestimate the necessary $C^{\star}$ value. While this means an overestimated number of CPU hours from the number given in §5.1, this is more desirable than underestimating the required number of CPU hours. It is worthwhile mentioning that the current data is for a different $Re_{\tau}$ and roughness geometry then what the predicted $C^{\star}$ empirical relation was developed with, showing a certain robustness to the method.

The roughness function values of $8\lesssim\Delta U^{+}\lesssim 11$ measured here place these pyramid cases in the fully rough regime, while the results and recommendations from the previous sections (§3–4) were conducted in the transitionally rough regime. Given that the pyramid cases were set up using these recommendations without any modifications and achieved reasonable agreement with the literature, this suggests that the results of this study can apply to both transitionally and fully rough flows.

6 Conclusions

A series of numerical experiments were conducted to investigate the fundamental dynamics of the minimal channel, with an emphasis on more efficiently simulating rough-wall flows in order to obtain the roughness function $\Delta U^{+}$. These experiments were performed using a finite difference code.

The streamwise length was investigated in §3.1 for a minimal channel. It appears that the largest log-layer eddy scales as $L_{x}=3L_{y}$, where the spanwise width is obtained from the guidelines in Chung et al. (2015), $L_{y}\gtrsim\max(100\nu/U_{\tau},k/0.4,\lambda_{r,y})$. However, for very narrow channels this implies a very short streamwise domain length, which may not be able to sustain the turbulent motions. The small spatial domain also makes obtaining converged statistics a very time-consuming process. As such, minimal channels would benefit from a minimum streamwise domain length of around 1000 viscous units and so a guideline for setting $L_{x}$ is $L_{x}\gtrsim\max(3L_{y},1000\nu/U_{\tau},\lambda_{r,x})$, where $\lambda_{r,x}$ would be some characteristic streamwise length scale of the roughness.

Two alterations to the outer-layer flow were investigated to see whether they affected the healthy near-wall flow. A half-height (open) channel was used in §3.2 which a slip-wall is positioned at the channel centre. This resulted in only a slight change to the mean velocity profile and turbulence intensities in both minimal and full-span channels, primarily affecting the wake region. The pressure fluctuations in the minimal channel were shown to be nearly an order of magnitude larger than in a full-span channel, however this was suggested to be as a result of the rapid pressure term which is related to the mean velocity gradient, $dU/dz$. When $dU/dz$ was instead artificially set to zero in the unphysical region above $z_{c}$, the pressure fluctuations from this altered velocity field were then in much better agreement with the full-span channel. A benefit of using a half-height channel is that only half the number of gridpoints are required, at the cost of running the simulation for twice as long to reach the same level of statistical convergence. It is up to the user to balance this trade off between memory and time.

A forcing model was applied to the outer layer in §3.3 to damp the fluctuations and reduce the centreline velocity. The near-wall flow was unaltered even though the Reynolds shear stress is zero above the location where the damping starts, $z_{d}$. This is similar to the filtering employed by Jiménez & Pinelli (1999), although the present damping model forces $U(z>z_{d})=U(z=z_{d})$. This allowed for an improvement in the computational time step of around 20–24%, although this did not extend to higher friction Reynolds numbers in the present simulations due to the grid using a cosine mapping in the wall-normal direction. The present data suggest that the location where the forcing starts should be $z_{d}\gtrsim\max(200\nu/U_{\tau},2z_{c})$. The first constraint is present for very narrow channels where $z_{c}^{+}$ is close to the wall, as the forcing will interfere with the near-wall flow otherwise.

A temporal sweep was conducted in §3.4 to see if the full roughness behaviour could be obtained by varying the roughness Reynolds number with time. A simple linear variation in the bulk velocity was attempted, with three different rates of change investigated. It was shown that the two slowest rates of change, with final pressure gradient parameter values of $\varDelta_{p,Re_{\tau}=180}=0.03$ and 0.07 could reasonably predict some of the $\Delta U^{+}$ vs $k_{s}^{+}$ curve. The largest value of $\varDelta_{p}=0.15$ resulted in substantial flow changes compared to the steady flow and so are not feasible. The early work of Perry & Joubert (1963) showed that $\Delta U^{+}$ was independent of an applied pressure gradient, and the present work goes on to provide a bound on $\varDelta_{p}$ for when this holds. This also suggests that pressure gradients, normally applied by a spatial acceleration of the flow in experimental facilities, can equally be simulated by a temporal acceleration. This is similar to the approach of Kozul et al. (2016) who represented a spatially developing boundary layer as temporally developing instead.

An eddy-counting argument was developed in §4 to analyse the statistical uncertainties of numerical simulations in minimal channels. The total number of $z_{c}$-sized eddies, $C^{\star}$, are counted over the duration of the simulation, which involves estimating their characteristic length and time scales. For a desired level of uncertainty, $\zeta$, in the roughness function, $\Delta U^{+}\pm\zeta$, the expected number of $z_{c}$-sized eddies that need to exist over the course of the simulation was estimated as $C^{\star}=(91.4/(\zeta\cdot z_{c}^{+}))^{2}$ as in (12). This means that for a known $z_{c}^{+}$, the user can determine how many $z_{c}$-sized eddies need to be observed through this empirical relation. The simulation run time, and hence number of CPU hours required, can then be estimated a priori. This would enable researchers to determine beforehand how best to allocate a limited number of CPU hours when studying rough-wall flows.

A case study of square-based pyramids (§5) was used to illustrate the minimal-span channel framework and the insights gained in this paper. The viscous dimensions of the pyramids were set to match those of Schultz & Flack (2009). The roughness functions from these simulations were overestimated compared to those of Schultz & Flack (2009), although underestimated compared to that of Di Cicca & Onorato (2016) who had a similar roughness geometry. The predicted $C^{\star}$ value required to obtain a desired level of statistical uncertainty was shown to be in reasonable agreement, if not somewhat conservative, with the data from the pyramid simulations. As such, the estimated CPU hours required for the simulation can be predicted reasonably accurately before performing the simulation. For this pyramid roughness, the minimal-span channel used nearly 20 times less CPU hours than the estimated cost of a full-span channel, as well as using substantially fewer CPUs.

Acknowledgements

This work was partly funded through the Multiflow program by the European Research Council. Computational time was granted under the Victoria Life Sciences Computational Initiative, which is supported by the Victorian Government, Australia.