2.1. Direct numerical simulations

Direct numerical simulations of turbulent forced convection over 3-D roughness have been performed in an open-channel configuration. The continuity, Navier–Stokes and passive scalar transport equations solved in this study are

(2.1)\begin{gather} \frac{\partial u_{i}}{\partial x_{i}}=0, \end{gather}

(2.2)\begin{gather}\frac{\partial u_{i}}{\partial t} + u_{j}\frac{\partial u_{i}}{\partial x_{j}}={-}\varPi\delta_{1i} -\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}+\nu\frac{\partial^{2} u_{i}}{\partial x_{j}\partial x_{j}} , \end{gather}

(2.3)\begin{gather}\frac{\partial \theta}{\partial t} + u_{j}\frac{\partial \theta}{\partial x_{j}}= \alpha \frac{\partial^{2} \theta}{\partial x_{j}\partial x_{j}} + Q , \end{gather}

where the velocity components in the streamwise ($x$), spanwise ($y$) and wall-normal ($z$) directions are $u$, $v$ and $w$, respectively, and $t$ is time. Buoyancy (due to gravity) and thermal expansion ($\beta$) properties of the fluid are neglected, i.e. ${g\beta \varTheta _\tau h/U_\tau ^2 = 0}$ (forced convection), and so temperature behaves as a passive scalar. The flow is driven through the open-channel domain by adding a spatially uniform body force to the right-hand side of the streamwise momentum equation, $\varPi \equiv -(1/\rho )\,\textrm {d}P/\textrm {d}\kern0.7pt x>0$ (where $\rho$ is density and $P(x)$ is the mean pressure). The hydrodynamic pressure that is solved for is the fluctuating (or periodic) component, $p$. To simulate forced convection, a body force was also added to the right-hand side of the transport equation for temperature, $Q=-u\,\textrm {d}T_w/\textrm {d}\kern0.7pt x>0$, where $u$ is the instantaneous streamwise velocity, $\textrm {d}T_w/\textrm {d}\kern0.7pt x$ is a prescribed mean wall-temperature gradient which drives the heat transfer and $\theta$ is the fluid temperature relative to the wall. This internal-heating technique can correspond to both a fluid moving along a cooled ($\mathrm {d}T_w/\mathrm {d}\kern0.7pt x<0$, $\varTheta >\varTheta _w$, $\varTheta _\tau >0$) or heated ($\mathrm {d}T_w/\mathrm {d}\kern0.7pt x>0$, $\varTheta <\varTheta _w$,$\varTheta _\tau <0$) wall and is a common forcing strategy amongst passive scalar channel-flow DNS (e.g. Kasagi, Tomita & Kuroda Reference Kasagi, Tomita and Kuroda1992; Alcántara-Ávila et al. Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021). Among both configurations, the dimensionless mean temperature $\varTheta ^+ \equiv (\varTheta - \varTheta _w)/\varTheta _\tau$ remains the same. Our present framework adopts a prescribed temperature fluctuation $\theta = 0$ (Dirichlet) boundary condition at the wall, which, with our present choice of temperature forcing $Q = -u\,\mathrm {d}T_w/\mathrm {d}\kern0.7pt x$, can be shown to correspond to a statistically uniform wall heat flux (e.g. Kasagi et al. Reference Kasagi, Tomita and Kuroda1992; Kays et al. Reference Kays, Crawford and Weigand2005) that mimics a Neumann boundary condition at the wall. For a statistically steady channel flow, volume integration shows that $U_\tau$ is related to the mean pressure gradient body forcing, $\varPi$, by $\varPi \equiv -(1/\rho )\,\mathrm {d}P/\mathrm {d}\kern0.7pt x = U_\tau ^2/h$ (e.g. Pope Reference Pope2000). This then enables a target value of $Re_\tau \equiv h U_\tau / \nu$ to be achieved by prescribing $\varPi$ to set $U_\tau$, for a given value of $h$ and $\nu$. For our present temperature forcing $Q = -u\,\mathrm {d}T_w/\mathrm {d}\kern0.7pt x$, volume integration shows that the average wall heat flux is given by $\langle q_w \rangle /(\rho c_p) = h U_b(\mathrm {d}T_w/\mathrm {d}\kern0.7pt x)$, where $U_b \equiv (1/h)\int ^h_0 U\, \mathrm {d}z$ is the bulk velocity (e.g. Alcántara-Ávila et al. Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021). This then sets the friction temperature for normalisations of temperature quantities, $\varTheta _\tau \equiv \langle q_w\rangle / (\rho c_p U_\tau )$. We imposed no-slip, impermeability conditions ($u=v=w=0$), a zero fluid–wall temperature difference ($\theta = 0$) at the bottom wall and free slip, impermeable and adiabatic boundary conditions at the top boundary. Other internal-heating techniques are possible, e.g. a spatially uniform body force (Kim & Moin Reference Kim and Moin1989; Pirozzoli, Bernardini & Orlandi Reference Pirozzoli, Bernardini and Orlandi2016; Peeters & Sandham Reference Peeters and Sandham2019). The average wall heat flux is defined here as $\langle q_{w} \rangle /(\rho c_{p})=\alpha \overline {\left \langle \partial \theta /\partial n|_{wall}\right \rangle }$, where $n$ is the local wall-normal direction, the overline denotes temporal averaging and angle brackets denote an integral over the surface wetted area, normalised on the total plan area. The thermal diffusivity $\alpha$ is a property of the working fluid and is related to the molecular Prandtl number, $Pr \equiv \nu /\alpha$. Presently, we fix the Prandtl number to match that of air at standard atmospheric conditions, $Pr = 0.7$. The simulations have been run at a constant, prescribed bulk velocity as opposed to setting the pressure gradient forcing, $\varPi$. For each simulation, the bulk velocity is set by a trial and error procedure such that the friction Reynolds number, $Re_{\tau }\equiv U_{\tau }h/\nu$, matches its target value to within $\pm 1\,\%$. This latter step is done to set the value of the roughness Reynolds number $k^+ \equiv kU_\tau /\nu$ which, for our simulations, will be determined by $k^+ \equiv (k/h)Re_\tau$, where $k/h$ is the inverse of the channel blockage ratio.

The governing equations (2.1) and (2.3) are solved using a low-storage third-order Runge–Kutta scheme (Spalart, Moser & Rogers Reference Spalart, Moser and Rogers1991) in conjunction with a fractional step algorithm (Perot Reference Perot1993), which was implemented in past work by Chung et al. (Reference Chung, Chan, MacDonald, Hutchins and Ooi2015). Spatial discretisation is achieved using a fully conservative fourth-order symmetry preserving scheme (Verstappen & Veldman Reference Verstappen and Veldman2003). The passive scalar transport (2.3) is solved using a fourth-order implicit scheme for the wall-normal diffusive term, whereas the streamwise and spanwise terms are treated with an explicit scheme and advected using the Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme (Leonard Reference Leonard1979). We have validated these discretisation schemes extensively in prior work applied to thermal convection configurations, demonstrating that both thermal dissipation (to 2 %) and the global heat transfer (to 0.1 %) are accurately captured provided the characteristic mesh size is at most two times larger than the Kolmogorov length scale $\eta \equiv (\nu ^3/\varepsilon )^{1/4}$ for ${Pr = 0.7}$ where $\varepsilon$ is the turbulent dissipation rate (table 2 Rouhi et al. Reference Rouhi, Lohse, Marusic, Sun and Chung2021). For our present study, the most stringent resolution requirements are in the vicinity of the roughness and Zhong et al. (Reference Zhong, Hutchins and Chung2023) reported that $\eta ^+(z=k) \approx (\kappa k^+)^{1/4}$ in the fully rough regime for our present 3-D sinusoids at fixed $\varLambda = 0.18$. For our $k^+ \approx 33$–94, this yields an estimate $\eta ^+(z=k) \approx 1.9\unicode{x2013}2.5$, which we amply resolve in the wall-normal direction ($\Delta z^+_w \lesssim 1$, cf. table 2). This is relaxed in the streamwise grid resolution for instance, where we have used $\Delta x^+ \approx 7.3\unicode{x2013}9.9 > \eta ^+(z=k)$. The roughness geometry is resolved using a direct-forcing variant of the immersed boundary method (IBM) based on a volume-of-fluid interpolation detailed in Rouhi, Chung & Hutchins (Reference Rouhi, Chung and Hutchins2019). Our present IBM code has been validated in earlier works for both momentum and heat transfer (Rouhi et al. Reference Rouhi, Chung and Hutchins2019; Zhong et al. Reference Zhong, Hutchins and Chung2023). Periodic boundary conditions were enforced in the streamwise and spanwise directions. A free-slip boundary condition was enforced on the top boundary of the computational domain. A no-slip boundary condition was enforced on the bottom boundary using the IBM.

Table 2.

Simulation parameters for the nineteen rough-wall cases performed in this study. Sets L, M and H contain the low solidity $(\varLambda =0.09)$, medium solidity $(\varLambda =0.18)$ and high solidity $(\varLambda =0.36)$ cases, respectively. The simulation parameters include: the friction Reynolds number $(Re_{\tau })$; the ratios of the open-channel mean height to the roughness semi-amplitude $(h/k)$, the roughness wavelength to the roughness semi-amplitude $(\lambda /k)$ and the roughness wavelength to the open-channel mean height $(\lambda /h)$; viscous-scaled roughness height or roughness Reynolds number $(k^+)$ and roughness wavelength $(\lambda ^+)$; viscous-scaled domain length $(L_{x}^+)$ and width $(L_{y}^+)$; number of grid points in the streamwise $(N_{x})$, spanwise $(N_{y})$ and wall-normal $(N_{z})$ directions; viscous-scaled streamwise ($\Delta x^+$), spanwise ($\Delta y^+$), wall-normal spacing within the roughness canopy $(\Delta z_{w}^+)$ and at the open-channel mean height $(\Delta z^+_{h})$; number of grid points per roughness wavelength in the streamwise direction ($\lambda /\Delta x$), spanwise direction ($\lambda /\Delta y$) and number of grid points from the roughness peak-to-trough ($k_p/\Delta z_w$); the ratio of the virtual-origin shift to the roughness height ($d/k$); the non-dimensional sampling period $(TU_{\tau }/\min (z_{c},h))$. Viscous-scaled quantities are non-dimensionalised using $U_{\tau }$ and $\nu$.

2.2. Roughness geometry and simulation parameters

The bottom boundary of the computational domain is covered with 3-D sinusoidal roughness (Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2015), which follows a 2-D height distribution of the form

(2.4)\begin{equation} z_w(x,y) = k\cos(2{\rm \pi} x/\lambda)\cos(2{\rm \pi} y/\lambda), \end{equation}

where $k$ is the semi-amplitude (i.e. half of the peak-to-trough height, $k_p/2$) and $\lambda$ is the wavelength. The mean absolute height $(k_{a})$ and root-mean-squared height $(k_{rms})$ of the present 3-D sinusoids are constant multiples of its semi-amplitude, where ${k=2 k_{rms}=({\rm \pi} ^{2}/4)k_{a}}$, whilst its skewness $(Sk)$ is equal to zero (Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2015). In addition, the effective slope ($ES$) of the present roughness can be expressed as ${ES=(8/{\rm \pi} )(k/\lambda )}$, which is the mean absolute streamwise gradient of the height distribution (Napoli, Armenio & De Marchis Reference Napoli, Armenio and De Marchis2008), and is related to frontal solidity through the formula $\varLambda \equiv ES/2 =(4/{\rm \pi} )(k/\lambda )$ (MacDonald et al. Reference MacDonald, Chan, Chung, Hutchins and Ooi2016). For a fixed roughness amplitude, a greater frontal solidity or, equivalently, a steeper $ES$, can be achieved by shortening the roughness wavelength. Alternatively, frontal solidity can be held constant by varying $k$ and $\lambda$ in proportion such that the geometric scaling factor $k/\lambda$ remains fixed (Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2018; Jelly et al. Reference Jelly, Chin, Illingworth, Monty, Marusic and Ooi2020). The present study directs attention to both techniques with the principal aim of understanding how systematic changes of frontal solidity affect heat and momentum transfer in rough-wall turbulent flow.

While turbulent forced convection over the present 3-D sinusoids has been investigated from smooth up to the fully rough regime in a minimal channel configuration by MacDonald et al. (Reference MacDonald, Hutchins and Chung2019), only a single value of frontal solidity of $\varLambda =0.18$ was considered in their study. The present study builds upon the past work of MacDonald et al. (Reference MacDonald, Hutchins and Chung2019) by extending the investigation of turbulent forced convection over 3-D sinusoids to two new values of frontal solidity, namely, $\varLambda =0.09$ and $\varLambda =0.36$, in addition to the intermediate value of $\varLambda =0.18$. Three-dimensional sketches of the sinusoidal roughness for each value of frontal solidity considered here are shown in figure 2. The three solidity cases $\varLambda = 0.09, 0.18$ and 0.36 are obtained by varying the roughness wavelength as $\lambda \approx 14.2k$, $7.1k$ and $3.6k$ shown in figures 2(a), 2(b) and 2(c), respectively.

Figure 2.

Open-channel-flow configuration with (a) low solidity $(\varLambda =0.09)$, (b) medium solidity $(\varLambda =0.18)$ and (c) high solidity $(\varLambda =0.36)$ sinusoidal roughness on the bottom boundary. The roughness mean plane is located at $z/h=0$ and is shown as the dashed black line (- - -) on the side views.

Table 2 details the simulation parameters for the nineteen rough-wall cases considered in this study. The simulations have been grouped into three separate sets according to their frontal solidity as low (L), medium (M) and high (H) corresponding to $\varLambda = 0.09$, 0.18 and 0.36. Set L contains the low solidity cases $\varLambda =0.09$, which falls into the sparse regime of $\varLambda < 0.18$, as shown by MacDonald et al. (Reference MacDonald, Chan, Chung, Hutchins and Ooi2016) for similar roughness geometry. This set has a blockage ratio of $h/k=36$, except for one case where $h/k=72$ to ensure that $Re_{\tau }\gtrsim 395$ for all cases to minimise low-Reynolds-number pressure gradient effects (Nickels Reference Nickels2004). The roughness Reynolds number, $k^+ \equiv kU_\tau /\nu$, for this set covers $5\lesssim k^+ \lesssim 67$. Set M contains the medium solidity cases $(\varLambda =0.18)$, which has a blockage ratio of $h/k=18$, except for one case where $h/k=36$, and covers the roughness Reynolds-number range ${11\lesssim k^+ \lesssim 94}$. Set H contains the high solidity cases ${\varLambda =0.36}$, which falls into the dense regime of $\varLambda > 0.18$ (MacDonald et al. Reference MacDonald, Chan, Chung, Hutchins and Ooi2016), which also has a blockage ratio of ${h/k=18}$, except for one case where ${h/k=36}$, and covers the roughness Reynolds-number range ${11 \lesssim k^+ \lesssim 94}$. Past work by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015) and Chung et al. (Reference Chung, Chan, MacDonald, Hutchins and Ooi2015) fitted the data of the present 3-D sinusoidal roughness with $\varLambda =0.18$ in the fully rough regime, obtaining an equivalent sand-grain roughness of $k_{s} \approx 4.1k$. Based on this previously obtained value, the group M cases span the range of $21 \lesssim k_{s}^+ \lesssim 268$ and encompass both the transitionally and fully rough regimes. The ratio of $k_s/k$ will be further investigated in § 3.3 after applying the virtual-origin shift. To reduce the computational cost of achieving fully rough conditions (which is needed to compute an equivalent sand-grain roughness for groups L and H), the spanwise domain width, $L_{y}$, was kept narrow in the spirit of the minimal channel for rough-wall flows proposed by Chung et al. (Reference Chung, Chan, MacDonald, Hutchins and Ooi2015). Although restricting the spanwise channel width produces unphysical results in the outer layer (e.g. the formation of an artificial wake), the near-wall flow remains unaffected and closely resembles a full-span simulation up to a wall-normal critical height of $z_{c}\approx 0.4 L_{y}$ (Flores & Jiménez Reference Flores and Jiménez2010). For all cases considered here, the spanwise domain width accommodates a minimum of three roughness wavelengths ($L_{y}/\lambda =3$), which is a factor of three wider than that used by MacDonald et al. (Reference MacDonald, Hutchins and Chung2019). The spanwise extent of our computational domain therefore satisfies the recommendations made by Chung et al. (Reference Chung, Chan, MacDonald, Hutchins and Ooi2015), namely $L_{y}\gtrsim \max (100\nu /U_{\tau },k/0.4,\lambda )$. In addition, our present choice of $L_y$ also ensures the minimal channel critical height, $z_c = 0.4L_y$ is greater than the roughness sublayer thickness estimated as $z_r \approx 0.5\lambda$ for our present roughness based on the wall-normal extent of roughness-induced dispersive motions (Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2018). The $z_c > z_r$ condition then ensures that the near-wall roughness sublayer flow is correctly resolved by our minimal channels. The critical height ${z_c}$ also applies to the temperature field, as shown by MacDonald et al. (Reference MacDonald, Hutchins and Chung2019) and Zhong et al. (Reference Zhong, Hutchins and Chung2023). The streamwise domain lengths used throughout this study also satisfy the condition $L_{x}\gtrsim \max (3L_{y},1000\nu /U_{\tau },\lambda )$, as recommended by MacDonald et al. (Reference MacDonald, Chung, Hutchins, Chan, Ooi and García-Mayoral2017). For the wall-normal grid, a uniform grid spacing of $\Delta z_{w}^+$ is used from the roughness trough $(z/k=-1)$ up to the roughness peak $(z/k=1)$, before being gently stretched to a maximum value of $\Delta z_{h}^+$ at the top boundary of the open-channel domain $(z/h=1)$ using a hyperbolic tangent function. The peak-to-trough height of the roughness canopy $(k_{p} \equiv 2k)$ is resolved using a minimum of eighty-five points ${(k_{p}/\Delta z_{w} \gtrsim 85)}$. Uniform grid spacing is used in the streamwise and spanwise directions, which ranges from ${3.3<\Delta x^+ < 9.9}$ and ${0.8<\Delta y^+ <5.4}$ for all cases considered here. In line with past recommendations made by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015), the roughness topography was resolved using a minimum of twenty-four points per wavelength in the streamwise direction ($\lambda /\Delta x=24$, where $\Delta x$ is the streamwise grid spacing) and a minimum of forty-eight points per wavelength in the spanwise direction ($\lambda /\Delta y=48$, where $\Delta y$ is the spanwise grid spacing). Further details are given in table 2. Validation of the current simulations against the past work of MacDonald et al. (Reference MacDonald, Hutchins and Chung2019) using a different (finite volume) solver at matched medium solidity (${\varLambda = 0.18}$) is included in Appendix A. Statistical quantities were computed using a space-then-time averaging methodology. Specifically, an instantaneous field variable, say, $\phi (x,y,z,t)$, was first spatially averaged across the $xy$-plane at each wall-normal height and then time averaged. Only in-fluid cells contribute to spatially averaged quantities within the roughness canopy $(-1< z/k<1)$, which corresponds to an intrinsic average (Gray & Lee Reference Gray and Lee1977). Throughout this manuscript, $U(z)$ and $\varTheta (z)$ denote the mean streamwise velocity and mean fluid temperature relative to the wall respectively. Statistical quantities reported herein were collected for a minimum of fifteen large-eddy turn over times, depending on the spanwise domain width and the friction Reynolds number.

4.1. Sheltered–exposed area decomposition

In this section, we present an analysis of techniques for estimating the exposed and sheltered regions, both through direct analysis of the DNS data and through the use of a ray-tracing model derived from the roughness geometry. Additionally, we propose a simplified fitting expression for predicting the area fractions ($A_e/A_w$ and $A_s/A_w$), based on the ray-tracing model. First, we estimate the sheltered regions directly from the DNS by whether the time-averaged ${u^+ =0}$ isosurface that lies above the roughness covers the underlying surface. If not, the flow is considered exposed. The red isosurface in figure 8(a) is at $\overline {u^+} = 0$ for $\varLambda = 0.18$ and $k^+ \approx 33$, showing the area covered within the recirculation region. The area partition produced from the isosurface is illustrated in the $xy-$plane in figure 8(d) where the sheltered and exposed regions are shown by the white and black contours respectively (labelled ‘S’ and ‘E’, respectively). A further possible method to estimate the sheltered and exposed regions is from the local, time-averaged viscous stresses ${\tau _\nu /\rho \equiv \nu (\partial u/\partial n|_{wall})}$. Here, we note that the local $\tau _\nu$ differs from the total wall shear stress, $\langle \tau _w \rangle$, from which the global friction velocity $U_\tau \equiv (\langle \tau _w \rangle / \rho )^{1/2}$ has been defined in § 1, as $\langle \tau _w\rangle$ encompasses both the viscous and pressure drag contributions for rough walls. From this $\tau _\nu$ method, the sheltered regions are determined as the area enclosed by the ${\tau _\nu = 0}$ contour (red line in figure 8b). The area partition produced from the former method in figure 8(d) and from the time-averaged contours of the ${\tau _\nu = 0}$ in figure 8(e) produce comparable results with no more than 7 % difference of the produced area partition for all cases.

Figure 8.

Methods employed for determining sheltered and exposed areas using the ${k^+ = 33}$ and ${\varLambda = 0.18}$ case as an example. Throughout, roughness peaks and valleys are marked by black circles and white squares, respectively. (a) Sheltered regions are determined by whether the ${\overline {u^+} = 0}$ isosurface (red) covers the underlying flow. If not, the flow is considered exposed. (b) Time-averaged contours of the viscous stress $\tau _\nu /\rho \equiv \nu \partial \bar {u}/\partial n|_{wall}$ (darker for increasing magnitude). Sheltered regions are determined as the area enclosed by the ${\tau _\nu = 0}$ contour (red line) and exposed otherwise. (c) The sheltering is modelled as a straight line starting from the crest and angled with sheltering angle ($\theta _s$). The exposed and sheltered areas are shown with black and white lines in (c), respectively. Area partition produced from (d) the isosurface $\overline {u^+} = 0$; (e) the ${\tau _\nu = 0}$ isocontour and; (f) the sheltering model. Sheltered regions are shown by white filled contours (labelled ‘S’) and exposed regions are shown in black (labelled ‘E’).

The predictions of the exposed and sheltered regions in figures 8(a,d) and 8(b,e) require the flow fields. To enable prediction of the area fraction without need for computing flow fields, we obtained the exposed and sheltered regions directly from the roughness geometry. For this technique (called ray tracing), we assume that the sheltering extends as a straight line from the crest to the consecutive roughness element with an angle $\theta _s$ measured from the horizontal line, as previously proposed by Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) and shown here in figure 8(c). For simplicity, the sheltering angle for the ray tracing, is assumed constant (i.e. independent of $k^+$ and $\varLambda$) at $\theta _s = 15^\circ$, although $\theta _s$ obtained from DNS data using the ${\overline {u^+} = 0}$ isosurface method (figure 8a) varied in the range ${\theta _s = 8^\circ }$–$15^\circ$ over the three solidity cases and approached a constant value of $\theta _s \approx 10^\circ$, $13^\circ$ and $15^\circ$ for $\varLambda = 0.09$, $0.18$ and 0.36, respectively, in the fully rough regime. The estimation of $\theta _s$ through the ${\overline {u^+} = 0}$ isosurface can be exemplified, for instance, by figure 7(a). In this case, $\theta _s$ corresponds to the angle of a straight-line approximation fitted to the sheltering contour (white line). Angles in the vicinity of this range (i.e. ${\theta _s \approx 15^\circ }$) are typical of those found for separated flows (Raupach Reference Raupach1992; Prasad & Williamson Reference Prasad and Williamson1997; Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016, Reference Yang, Stroh, Chung and Forooghi2022). Figure 8(f) shows that the area partition from the straight-line sheltering assumption with ${\theta _s = 15^\circ }$ can provide a reasonable prediction of the area fractions compared with the other methods. Figure 9 shows the area fraction from the time-averaged contours of ${\tau _\nu = 0}$ and the ray tracing of the straight-line sheltering for the three solidity cases at fixed $k^+ \approx 33$. The validity of the straight sheltering assumption (as shown in figure 9e–f) for $\varLambda = 0.18$ and 0.36 is demonstrated by the good approximation of the area partition compared with figure 9(b,c). However, for the case of low solidity with $\varLambda = 0.09$, the straight-line sheltering assumption appears to underestimate the sheltered area (as illustrated in figure 9d) in contrast to the area partition in figure 9(a). This may be attributed to either the inaccurate estimation of $\theta _s$ for the low solidity, or the assumption of no lateral spreading, as will be discussed below.

Figure 9.

Comparison between sheltered–exposed area partitions at varying $\varLambda$ for (a–c) the present DNS at $k^+ \approx 33$ computed obtained from the $\tau _\nu =0$ partition and (d–f) the present model adopting a fixed sheltering angle, $\theta _s \approx 15^\circ$. White patches correspond to the sheltered (‘S’) area, while black patches are for the exposed (‘E’) area. Black circles and white squares mark roughness peaks and valleys, respectively.

The straight-line sheltering assumption clearly does not account for the lateral spreading that is evident in the flow field data (e.g. figure 9a–c). This lateral spreading could also be included into the current model following the same approach as Millward-Hopkins et al. (Reference Millward-Hopkins, Tomlin, Ma, Ingham and Pourkashanian2011), who developed an empirical relation to relate the lateral spreading angle to the roughness height and width for uniform cube array roughness. However, the current technique appears to provide reasonable area fraction prediction with minimal assumptions.

The previous results demonstrate that, while the ray-tracing technique is a reliable tool for predicting the area fraction of 3-D sinusoidal roughness, it lacks a closed-form analytical expression, which may make its implementation in a predictive model cumbersome. For our present 3-D sinusoids, the area ratio $A_e/A_w$ from ray-tracing may instead be well approximated by the following empirical relation:

(4.3)\begin{equation} \frac{A_e}{A_w} = \mathrm{min} \left(1.0,\ 0.5\left(\frac{\tan(\theta_s)}{\varLambda}\right)^{1/2} \right), \end{equation}

where $\theta _s$ is the sheltering angle, as illustrated in figure 8(c) and assumed throughout to be $15^\circ$. A limiting condition when all wetted areas are exposed ${A_e = A_w}$ is incorporated into (4.3). The sheltered area fraction will follow from $A_s/A_w \equiv 1-A_e/A_w$. Figure 10 shows the good agreement of the area fraction expression (4.3) with the DNS data. The DNS data in figure 10 $A_e/A_w$ are computed from the $\tau _\nu = 0$ contour method (figure 8b,e). Figure 10 also shows the ray-tracing representation (grey lines) of the exposed to wetted area ratio for validation of the fit (4.3). The comparison with our present DNS also validates one crucial simplification used in our area fraction expression of (4.3), namely, that this area fraction is independent of $k^+$, which appears to approximately hold true for larger $k^+$. Thus, hereafter in this section, the area fraction expression (4.3) will be used to determine ($A_e/A_w$ and $A_s/A_w$).

Figure 10.

Exposed wetted area fraction $A_e/A_w$. The markers show area fractions determined from the DNS, based on the $\tau _v = 0$ contour method, with the same marker styles as figure 6. These data are compared against the ray-tracing algorithm (grey lines) and the area ratio expression (4.3) (black lines) for varying choices of the sheltering angle $\theta _s$, with a best fit of $\theta _s\approx 15^\circ$.

4.2. Sheltered–exposed heat-transfer models

Here, we will provide details on the exposed and sheltered heat-transfer models for $C_{h,e}$, $C_{h,s}$ required in (4.2) for prediction. For the present DNS, $C_{h,e}$ and $C_{h,s}$ are computed using the exposed–sheltered decomposition obtained from the $\tau _\nu = 0$ contour method (e.g. figure 8b,e), consistent with figures 9 and 10. Figure 11 summarises the content we will present in this section. Namely, models for $C_{h,e} = f(k^+,\varLambda )$ and $C_{h,s} = f(k^+,\varLambda )$ and a comparison with the current DNS and that of Zhong et al. (Reference Zhong, Hutchins and Chung2023) (markers).

Figure 11.

(a) The exposed heat-transfer coefficient $C_{h,e}$ for the present DNS (markers with lines). The marker styles are the same as in figure 6. Direct numerical simulations at fixed $\varLambda = 0.18$ for varying $\textit {Pr}$ from Zhong et al. (Reference Zhong, Hutchins and Chung2023) have also been included. The exposed heat-transfer model $C_{h,e}\textit {Pr}^{2/3} = 0.5(k^+)^{-1/2}\varLambda ^{1/2}$ is shown by the coloured lines matching each solidity in the DNS. (b) Sheltered heat-transfer coefficient compensated by an empirical Prandtl number scaling $C_{h,s}Pr^{0.45}$. The present model approximates (4.5) this behaviour as both a $C_{h,s}\approx 0.012 \textit {Pr}^{-0.45}$ (horizontal black line) and zero.

To model the exposed heat transfer, we will adopt a Reynolds-analogy-type model, $C_{h,e}\propto (k^+)^{-1/2}\varLambda ^{1/2}\textit {Pr}^{-2/3}$. The model hinges on the phenomenological ideas of Owen & Thomson (Reference Owen and Thomson1963) and Yaglom & Kader (Reference Yaglom and Kader1974), where they postulated that the local viscous–conductive regions close to the roughness elements behave analogously to a flat-plate laminar boundary layer. The intuition is restricted to the fully rough regime (high-$k^+$) such that a scale separation between the roughness height and viscous scale emerges $k \gg \nu /U_\tau$. Owing to this scale separation, they envision the flow locally as having attached thermal and viscous boundary layers, developing over the roughness height $k$. The consequence in envisioning this phenomenology, as argued by Owen & Thomson (Reference Owen and Thomson1963) and Yaglom & Kader (Reference Yaglom and Kader1974), is that we may then generalise the Reynolds-analogy prediction for heat transfer $C_{h,x}\propto \textit {Re}_x^{-1/2}\textit {Pr}^{-2/3}$ (Kays et al. Reference Kays, Crawford and Weigand2005) to the case of a rough wall. Here, $\textit {Re}_{x}\equiv xU_\infty /\nu$ is a Reynolds number based on the streamwise development length, $x$ and free-stream velocity, $U_\infty$, while $C_{h,x} \equiv q_{w,x}/(\rho c_p U_\infty \varTheta _\infty )$ is the dimensionless heat transfer (Stanton number) defined on the free-stream temperature relative to the wall, $\varTheta _\infty$. In generalising to a rough wall, Owen & Thomson (Reference Owen and Thomson1963) and Yaglom & Kader (Reference Yaglom and Kader1974) set the streamwise fetch as $x \propto k$ and the velocity and temperature scales, $U_\infty \propto U_\tau \propto U_k$, $\varTheta _\infty \propto \varTheta _k$. Adopting these relations, this generalises $\textit {Re}_x \equiv xU_\infty /\nu$ to the roughness Reynolds number $k^+ \equiv kU_\tau /\nu$ and the heat transfer $C_{h,x} \equiv q_{w,x}/(\rho c_p U_\infty \varTheta _\infty )$ to a rough-wall heat-transfer coefficient, say, $C_{h,k}$, yielding the predictions: $C_{h,k} \propto (k^+)^{-1/2}\textit {Pr}^{-2/3}$.

Our present work will adopt two additional modifications to the original form put forward by Owen & Thomson (Reference Owen and Thomson1963) and Yaglom & Kader (Reference Yaglom and Kader1974), such that the prediction may be sensitised to the solidity, $\varLambda$. First, rather than take the Reynolds-analogy model to describe the total heat transfer, $C_{h,k}$, we instead only take it to describe the exposed heat transfer, $C_{h,e}$. The intuition underpinning this decision is based on the observations from figure 7, where we isolate that it is only in the windward, exposed regions where the flow remains attached and tends to experience high-shear flow, correlating to higher heat transfer. Consequently, we may expect ideas from a canonical wall-attached boundary-layer flow to be applicable in these regions. Notably, this behaviour is distinct when compared with sheltered regions which appear to follow a different heat-transfer scaling entirely. Second, rather than taking $x\propto k$ as the streamwise fetch, we instead take $x \propto \lambda$. The intuition underpinning this choice comes from recognising that $x$ is understood as the wall-parallel development length of the boundary layer. Thus, rather than taking $x$ to be analogous to the roughness height $k$ in the case of a rough wall, the more appropriate analogue to take for $x$ would be the local arc length along the roughness topography over which the local boundary layer may develop. In the interest of simplicity, we may approximate as the in-plane roughness length scale $\lambda$. Adopting these changes, the model changes from $C_{h,k}\propto (k^+)^{-1/2}\textit {Pr}^{-2/3}$ to $C_{h,e}\propto (\lambda ^+)^{-1/2}\textit {Pr}^{-2/3}$, since $\lambda \propto k/\varLambda$ by definition

(4.4)\begin{equation} C_{h,e} = 0.5(k^+)^{{-}1/2}\varLambda^{1/2}\textit{Pr}^{{-}2/3}, \end{equation}

where the $0.5$ pre-factor is a tuning constant. As seen in figure 11(a), (4.4) is able to capture the scaling trends of our present DNS reasonably well, best seen for our high-$\varLambda$ cases (${\gtrsim }0.18$). Although the model may appear to work well for transitionally rough values of $k^+$ (i.e. non-fully rough), this is not by design. Strictly speaking, our $C_{h,e}$ model is best applied to the fully rough regime, as the greater $k \gg \nu /U_\tau$ scale separation helps to ensure that the postulated phenomenology of a local laminar boundary layer is upheld. From figure 11(a), however, we see that our $C_{h,e}$ model still performs well in the transitionally rough regime at $k^+ \approx O(10)$. Over a smooth wall, the viscous sublayer (i.e. the local laminar boundary layer) has a characteristic thickness of $5\nu /U_\tau$ (Pope Reference Pope2000). Consequently, moderate transitionally rough values of $k^+ \approx O(10)$ may not provide significant scale separation between the roughness size $k$ and the boundary-layer thickness $5\nu /U_\tau$. It is moreover expected that our model should fail if low-$k^+$ values are considered such that $k$ is smaller than or of comparable size to $5 \nu /U_\tau$. The roughness elements in this instance are submerged below or of comparable size to the viscous sublayer thickness (e.g. Luchini Reference Luchini1996). Consequently, this would invalidate the phenomenology of a local laminar boundary layer developing over the roughness topography which (4.4) relies upon, owing to insufficient scale separation between $k$ and $5 \nu /U_\tau$. To verify the validity of the Reynolds-analogy $\textit {Pr}^{-2/3}$ scaling, we have also included in figure 11(a) data from the DNS of Zhong et al. (Reference Zhong, Hutchins and Chung2023). The study considered an identical 3-D sinusoidal roughness with fixed $\varLambda = 0.18$, while varying $k^+$ and $\textit {Pr}$. As evidenced by figure 11(a), the $Pr^{-2/3}$ scaling yields an excellent collapse for the exposed heat transfer, $C_{h,e}$, within the limited investigated range of $0.5 \lesssim Pr \lesssim 2.0$.

We now consider a model for the sheltered heat transfer, $C_{h,s}$. Figure 11(b) shows $C_{h,s}$ for the present DNS as well as the $\varLambda = 0.18$, varying $\textit {Pr}$ data from Zhong et al. (Reference Zhong, Hutchins and Chung2023). Here, we observe a behaviour where $C_{h,s}$ is held approximately constant with $k^+$ in the fully rough regime (high-$k^+$), with a mild dependence on the constant value with respect to $\varLambda$ as well as a dependency on $\textit {Pr}$. A similar behaviour can be observed when examining the local variations in figure 7(d–f) or (p–r) where the heat transfer is held approximately uniform in sheltered regions. Of notable difference, however, is how the exposed heat transfer tends to provide the dominant contribution to the total when examining figure 7(p–r), but when taking integrated measures across the entire surface area (i.e. $C_{h,e}$ and $C_{h,s}$), the trends of figure 11(a,b) suggest both $C_{h,e}$ and $C_{h,s}$ will approach similar orders of magnitude for larger $k^+$. Slight qualitative differences will also be present in comparing the sheltered regions of figure 7 and the sheltered heat transfer $C_{h,s}$ in figure 11(b), as $C_{h,s}$ is computed based on sheltered regions from the $\tau _\nu = 0$ contour (figure 10b,e) whereas, in figure 7, sheltered regions are identified as reversed-flow regions ($\overline {u^+}=0$). In light of these observations and in the interest of a model of maximal simplicity, we will test two separate models for the sheltered heat transfer

(4.5a,b)\begin{equation} C_{h,s}= {0.012\textit{Pr}^{{-}0.45}}, \quad C_{h,s} = 0. \end{equation}

The first option $C_{h,s} = 0.0012\textit {Pr}^{-0.45}$ is an empirical fit to the medium–high solidity data (figure 11b, black line). The second option of $C_{h,s} = 0$ coincides with the physical assumption that the total heat transfer, $C_{h,k}$, may be predicted solely through the exposed contribution in (4.2) alone and is one of the simplest approximations one may have arrived at based on the $q_w$ observations of figure 7(p–r). It is worth noting that this approximation may not align with the integrated sheltered heat transfer $C_{h,s}$ (figure 11b), as $C_{h,e}$ and $C_{h,s}$ can approach similar orders of magnitude at higher $k^+$. Despite this, we will still examine the case where $C_{h,s}= 0$ in order to assess its impact on the total rough-wall heat-transfer coefficient, $C_{h,k}$.

4.3. Model comparison data

In this section, we compare our present model with the $C_{h,k}$ measured from the DNS, having now given various expressions in (4.3), (4.4) and (4.5) as required in (4.2) for prediction. Figure 12 presents the comparison against our present DNS. Adopting the choice of $C_{h,s}=0.012\textit {Pr}^{-0.45}$ (solid lines), our model captures the trend of the data reasonably well, best seen for the higher solidity cases. Worth noting too is the performance of the model prediction when neglecting any sheltered contribution, $C_{h,s}=0$ (dashed lines). Here, we see that the dominant scaling behaviour of $C_{h,k}$ with respect to $k^+$ can still be maintained despite this, and does not suffer too harshly in performance for $k^+ \approx 30\unicode{x2013}100$. This result then showcases that one may obtain a good estimate for the total heat transfer, $C_{h,k}$, by solely considering contributions from exposed regions of the flow for this $k^+$ range. We emphasise the final point in particular regarding the $k^+$ range, because as seen in figure 12(a) the $C_{h,s} = 0.012\textit {Pr}^{-0.45}$ and $C_{h,s}=0$ choices eventually produce disparate asymptotic behaviours for $C_{h,k}$ for higher $k^+$, indicating that sheltered heat transfer may become prominent further into the fully rough regime.

Figure 12.

(a) The rough-wall heat-transfer coefficient $C_{h,k}$ for our present DNS (coloured markers) compared against our present model (coloured lines) from (4.2) using the various proposals of (4.3)–(4.5) and a fixed sheltering angle, $\theta _s \approx 15^\circ$. The solid lines adopt a constant sheltered heat transfer $C_{h,s} = 0.012\textit {Pr}^{-0.45}$ while the dashed lines neglect the sheltered heat-transfer contribution ($C_{h,s} = 0$). The solid lines adopt a constant sheltered heat transfer $C_{h,s} = 0.012 Pr^{-0.45}$ while the dashed lines neglect the sheltered heat-transfer contribution ($C_{h,s} = 0$), resulting in less accurate predictions of $C_{h,k}$. ($b$) Value of $C_{h,k}$ plotted with respect to $\varLambda$ (darker coloured lines correspond to increasing $k^+$). The black lines connect DNS data at matched values of $k^+$. Only model lines with constant $C_{h,s}=0.012\textit {Pr}^{-0.45}$ are shown for clarity.

A further comment can be made regarding the range of validity with respect to $\varLambda$ of our present model. The model framework relies on a clear distinction between exposed and sheltered regions which for particular limits of $\varLambda$ may be ill defined. For sparse–wavy regimes ($\varLambda \to 0$), smooth-wall conditions potentially absent of any flow separation will be approached (Raupach et al. Reference Raupach, Antonia and Rajagopalan1991; MacDonald et al. Reference MacDonald, Chan, Chung, Hutchins and Ooi2016), which would invalidate the concept of sheltered flow. Although for our $\varLambda = 0.09$ case, flow separation can be observed locally (cf. figure 7a,j, $\overline {u^+} \leq 0$), the differences in the sheltered and exposed heat-transfer magnitudes are not as pronounced compared with the higher solidity cases (compare figure 7p with q,r). Consequently, the distinction between sheltered and exposed behaviour is perhaps not as evident for our $\varLambda = 0.09$ case, which could explain why our present model is less successful for $\varLambda = 0.09$.

It is instructive to demonstrate that our framework for predicting the rough-wall heat-transfer coefficient, $C_{h,k}$, is also capable of handling different roughness topographies. Here, we have applied our framework to the experimental 2-D rib data of Webb et al. (Reference Webb, Eckert and Goldstein1971). Owing to the different geometry, applying the ray-tracing exercise as was done for our present sinusoids in § 4.1 will yield different expressions for the exposed–sheltered area fractions, $A_e/A_w$, $A_s/A_w$ required for prediction of $C_{h,k}$ from (4.2). These details have been left to Appendix B, and, identically to our present sinusoids, yield expressions solely determined by $\varLambda$ and an assumed sheltering angle $\theta _s$. Recall that the constant prefactors for the exposed and sheltered heat-transfer coefficient ${C_{h,e} = 0.5(k^+)^{-1/2}\varLambda ^{1/2}\textit {Pr}^{-2/3}}$, $C_{h,s} = 0.012\textit {Pr}^{-0.45}$ and sheltering angle $\theta _s = 15^\circ$ were obtained by tuning to our present DNS (figure 11), where the partitions for $C_{h,e}$ and $C_{h,s}$ were readily accessible. Presumably, these constants are likely dependent on the specific roughness topography and need to be tuned, but for the experimental data of Webb et al. (Reference Webb, Eckert and Goldstein1971), only the total coefficient $C_{h,k}$ is available. Consequently, we make no attempt to tune these coefficients, instead keeping them the same as the ones used for our present sinusoids. In figure 13, we show a comparison between our present model and the $C_{h,k}$ of both our present DNS and data of Webb et al. (Reference Webb, Eckert and Goldstein1971). Our model shows excellent agreement in capturing the fully rough (high-$k^+$) data of Webb et al. (Reference Webb, Eckert and Goldstein1971), despite not having made any changes to the $C_{h,e}$, $C_{h,s}$ constant prefactors and sheltering angle, $\theta _s = 15^\circ$. The trends of $C_{h,k}$ with respect to $\varLambda$ are also well captured. These results then demonstrate that our sheltered–exposed framework has the potential to generalise well to arbitrary roughness geometries. It is worth mentioning that the current model yields accurate predictions for the 2-D rib data of Webb et al. (Reference Webb, Eckert and Goldstein1971), even when low solidities of $\varLambda = 0.025$ are considered. This contrasts with the discrepancy noted between the model predictions and the DNS case of low solidity ($\varLambda = 0.09$) for the present 3-D sinusoids. This distinction might be attributed to the presence of sharp corners in the 2-D rib geometry, which forces the flow to separate and helps to ensure that there exist distinct sheltered and exposed regions: a key assumption of our present model.

Figure 13.

The rough-wall heat-transfer coefficient $C_{h,k}$ at $\textit {Pr} = 0.7$ for the present DNS (markers with lines) and the 2-D rib experiments of Webb et al. (Reference Webb, Eckert and Goldstein1971) (square markers). The model predictions at each $\varLambda$ (solid lines for present 3-D sinusoids, dashed lines for 2-D ribs) all employ the same model constants, $C_{h,e} = 0.5(k^+)^{-1/2}\varLambda ^{1/2}\textit {Pr}^{-2/3}$, $C_{h,s}=0.012\textit {Pr}^{-0.45}$, and sheltering angle $\theta _s = 15^\circ$. The area models employed for the 2-D ribs are given by (B1)–(B2).

Next, we discuss the effective scaling behaviour of our model $C_{h,k} \sim (k^+)^{-p}$, where the correct value to take for the exponent $p$ has remained a topic of debate (e.g. Li et al. Reference Li, Rigden, Salvucci and Liu2017, Reference Li, Bou-Zeid, Grimmond, Zilitinkevich and Katul2020; Chung et al. Reference Chung, Hutchins, Schultz and Flack2021; Zhong et al. Reference Zhong, Hutchins and Chung2023). Amongst theoretical proposals of note are the $p=1/2$ Reynolds-analogy-type proposal of Owen & Thomson (Reference Owen and Thomson1963), Yaglom & Kader (Reference Yaglom and Kader1974) and the $p=1/4$ surface-renewal scaling of Brutsaert (Reference Brutsaert1975). It has been argued that this scaling exponent may perhaps depend on the roughness topography (Li et al. Reference Li, Bou-Zeid, Grimmond, Zilitinkevich and Katul2020) and in some cases it has been hinted that this exponent varies with spatial location (Meinders, Hanjalic & Martinuzzi Reference Meinders, Hanjalic and Martinuzzi1999; Nakamura, Igarashi & Tsutsui Reference Nakamura, Igarashi and Tsutsui2001; Li et al. Reference Li, Bou-Zeid, Anderson, Grimmond and Hultmark2016). This latter intuition is one that directly coincides with our proposed model, where, namely, we have proposed that in exposed regions, $C_{h,e}\sim (k^+)^{-1/2}$ and for sheltered, $C_{h,s} \sim (k^+)^0$. Recall from (4.2) that the total heat-transfer coefficient follows from a weighted sum (based on area fractions dependent on $\varLambda$) between $C_{h,e}$ and $C_{h,s}$. In this sense, our $C_{h,k}$ model can be interpreted as a ‘blending’ function between a $C_{h,k}\sim (k^+)^{-1/2}$ and $C_{h,k} \sim (k^+)^{0}$ scaling behaviour, where the aggregated behaviour between these two scalings can potentially emerge as the surface-renewal scaling $C_{h,k}\sim (k^+)^{-1/4}$ depending on the particular value of $\varLambda$ and $k^+$ (figure 13). Thus, in contrast to past rough-wall heat-transfer theories, our present model framework proposes that rough-wall heat transfer is determined not by any singular mechanism (e.g. Reynolds analogy or surface renewal) but rather by a blending of distinct behaviours, where presently we identify these distinct behaviours as a sheltered and exposed behaviour.

Finally, we end this section by highlighting that our rough-wall heat-transfer coefficient has been defined on crest values, $C_{h,k}\equiv 1/(U^+_k \varTheta ^+_k)$. Although definitions based on crest values may be common in the roughness literature (e.g. Dipprey & Sabersky Reference Dipprey and Sabersky1963; Webb et al. Reference Webb, Eckert and Goldstein1971; Macdonald et al. Reference Macdonald, Griffiths and Hall1998), a more relevant quantity for practitioners to use in prediction may be an integrated or full-scale heat-transfer coefficient, e.g. the bulk Stanton number $C_h$. Therefore, in Appendix C we have provided comprehensive details on how our proposed $C_{h,k}$ model can be extended to enable full-scale prediction of the integrated heat transfer, $C_h$.