Characterizing the Roughness in Channel Flows Using Direct Numerical Simulations

We simulated a steady channel flow with direct numerical simulation (DNS) in which we solve the incompressible Navier–Stokes equationssubject to the continuity equationwhere $t$ = time; $u_i$ = velocity vector; $x_j$ = coordinate vector; ${\rho }_0$ = density of the fluid; $p$ = pressure; $\nu$ = kinematic viscosity of the fluid; ${\mathrm{\Pi }}_c$ = driving pressure gradient; ${\delta }_{ij}$ = Kronecker delta function; and $F_{\mathrm{IBM}}$ = immersed boundary force used to model the effect of the roughness elements. The coordinate axes $x_1$, $x_2$, and $x_3$ are aligned in the streamwise, spanwise, and vertical directions, respectively. The channel is periodic in the streamwise and spanwise directions, whereas a no-slip boundary condition is imposed at the bottom wall where the roughness elements are located. The top wall has boundary conditions given bywhere $H$ = channel depth.

The governing equations were solved with a direct-numerical simulation code that had been extensively validated in previous studies of turbulent channel flows (Lozano-Durán and Bae 2016, 2019; Patil and Fringer 2022). In the code, all terms in the momentum equation are discretized in space with second-order accurate finite-difference schemes. In time, the fractional-step method (Kim and Moin 1985) is used along with a third-order accurate Runge–Kutta time-advancement scheme (Rai and Moin 1991). The advection terms are discretized in time explicitly, whereas the viscous terms are discretized implicitly to eliminate the stability constraint associated with viscous diffusion where the vertical grid resolution is refined. The pressure-Poisson equation was solved using a direct-Fourier–based method, as detailed in Costa (2018). Fourier transforms are applied in the $x_1$ (streamwise) and $x_2$ (spanwise) directions, which have uniform grid spacing, whereas a finite-difference method is applied in the vertical, which can have variable grid spacing. These discretizations result in tridiagonal systems for the horizontal Fourier modes that are solved efficiently with a tridiagonal solver. A bumpy wall was introduced with a direct forcing immersed boundary method based on the method proposed by Scotti (2006). The bumpy wall was generated by placing randomly oriented ellipsoids with fixed semiaxes such that the mean roughness height could be estimated a priori through the roughness function (${\psi }_r$), which is the area fraction as a function of height, as shown in Fig. 2. Here, the area fraction at a given vertical coordinate is given by $A_f(x_3)=2{\pi }^2{H}^2[1-{\psi }_r(x_3)]$. Additional details of the computational framework can be found in Patil and Fringer (2022). As shown by Jiménez and Moin (1991) and Flores and Jiménez (2010), the near-wall, nonlinear turbulent kinetic energy production cycle maintains “healthy turbulence” for wall-bounded flows. The study by Jiménez and Moin (1991) provided crucial insights into the statistical flow properties and elucidated the geometric requirements for “healthy turbulence” in channel flow geometries, also called minimal-span channels. This concept of the minimal-span channel was subsequently used to understand the mean flow drag without resolving the entire velocity profile (Chung et al. 2015; MacDonald et al. 2017). As the name indicates, the minimal-span channels limit the domain size in the spanwise direction such that the flow domain resolves the minimal dynamics (i.e., the interaction of two streamwise streaks) required to correctly resolve the near-wall region that is responsible for most of the turbulence production (Jiménez and Moin 1991; Flores and Jiménez 2010; Pope 2000). As a result, by correctly tuning the domain size in the spanwise direction, the mean velocity profile can be resolved adequately up to $x_3^{+}\equiv x_3{u}_{*}/\nu \approx 160$ (Jiménez and Moin 1991; Flores and Jiménez 2010; Chung et al. 2015; MacDonald et al. 2017).

We define the drag coefficient aswhere $u_{*}\equiv \sqrt{\tau /{\rho }_0}$ is the friction velocity; $\tau$ = time-averaged bottom stress averaged over the bottom boundary; and $U$ = domain-integrated and time-averaged streamwise velocity. The friction velocity is fixed by choosing the driving pressure gradient ${\mathrm{\Pi }}_c=u_{*}^2/H$, which in turn ensures that the bottom stress is given by $\tau ={\rho }_0{u}_{*}^2$. As a result, because $u_{*}$ is fixed, the drag coefficient changes due to changes in $U$. In this problem, there are seven relevant parameters: the bottom stress ($\tau /{\rho }_0$), velocity ($U$), channel depth ($H$), kinematic viscosity of the fluid ($\nu$), and the three semiaxis lengths of the ellipsoid ($\alpha {\overline{k}}_s$, $\beta {\overline{k}}_s$, and $\gamma {\overline{k}}_s$). Additionally, there are two rotation angles (uniformly distributed) that are required to define the Euler angle rotations of the ellipsoids. However, because there are a relatively large number of roughness elements, it is assumed that the effect of sampling the rotation angles is not substantial, as shown in Yuan and Piomelli (2014). Thus, using the Buckingham-pi theorem, the drag coefficient can be shown to have the functional dependencewhere $\mathrm{Re}\equiv \mathrm{UH}/\nu$ is the Reynolds number; $C_o$ = Corey shape factor defined in Eq. (3); $H/{\overline{k}}_s$ = blocking factor; and $S_p\equiv {(\alpha \beta /{\gamma }^2)}^{1/3}$ is the sphericity of the ellipsoids. In Eq. (8), the last three terms on the right-hand side of the equation depend on the statistical properties of the bed. Therefore, rather than conducting an exhaustive study of the effects of the bed parameters, we focused on the effects of $C_o$ while holding $H/{\overline{k}}_s$ fixed and then related the effective bottom roughness to the statistical properties of the bed. The sphericity ($S_p$) is not held constant, although we minimized its effects by ensuring $S_p\ge 0.84$.

To understand the effect of $C_o$, we chose friction Reynolds numbers 350 and 700, fixed $H/{\overline{k}}_s=13.15$, and ran a set of full-span simulations along with a series of minimal-span channel flows to extend the range of $C_o$ while limiting the computational cost. To directly compare the flat wall cases to the bumpy wall cases, Jiménez (2004) recommended $H/{\overline{k}}_s>40$, such that the change in the effective depth does not significantly affect the comparison. However, because all the bumpy wall cases have identical values of $H/{\overline{k}}_s$, they can be compared directly. The relation ${⟨·⟩}^{+}\equiv ⟨·⟩u_{*}/\nu$ indicates nondimensionalization using wall units unless specified otherwise. As shown in Table 1, we ran a set of 13 simulations to understand the effect of changing $C_o$ and ${\mathrm{Re}}_{*}$ on the mean flow drag.

The full-span channels had dimensions $2\pi H\times \pi H\times H$ and were discretized using $768\times 512\times 256$ grid points in the streamwise, spanwise, and vertical directions, respectively. Uniform grid spacing was used to resolve the roughness region, beyond which hyperbolic tangent grid stretching was used. For ${\mathrm{Re}}_{*}=700$, this gave $\mathrm{\Delta }{x}_1^{+}=5.73$, $\mathrm{\Delta }{x}_2^{+}=4.30$, and $\mathrm{\Delta }{x}_{3,\min }^{+}=0.66$ over the roughness region and $\mathrm{\Delta }{x}_{3,\max }^{+}=10.0$ at the top of the channel. The dimensions of the minimal-span channels were $2\pi H\times 200\mathrm{\Delta }{x}_2^{+}\times H$ and discretized using $768\times 64\times 256$ grid points in the streamwise, spanwise, and vertical directions, respectively. For all simulations, a time-step size of $\mathrm{\Delta }{t}^{+}\equiv u_{*}^2\mathrm{\Delta }t/\nu =0.045$ was used, based on a maximum Courant number of 0.4. All bumpy wall simulations corresponded to hydraulically transitional flow conditions, that is, $4\le u_{*}{\overline{k}}_s/\nu \le 70$, with $u_{*}{\overline{k}}_s/\nu =26.6$ and $u_{*}{\overline{k}}_s/\nu =53.2$ for ${\mathrm{Re}}_{*}=350$ and ${\mathrm{Re}}_{*}=700$, respectively. Ideally, for DNS, higher-order methods are preferred, although the deficiencies inherent to the second-order accurate method can be alleviated with grid refinement for the requisite scales of interest, as discussed in Moin and Verzicco (2016). For the simulations carried out in this paper, the region that encapsulates the roughness features was uniformly discretized with a vertical grid resolution of $\mathrm{\Delta }{x}_3^{+}=0.6<1.0$. This grid resolution accurately resolved the flow features that correctly predict the Reynolds-stress dynamics such that the relative residual in computing the turbulent kinetic energy budget was less than 5%. Additionally, implementing boundary conditions and the immersed boundary method are relatively straightforward for second-order numerical stencils when compared to higher-order stencils or spectral methods. As a result, despite the second-order accurate spatial discretization, the present computational framework ensured accurate simulation of turbulent flows without introducing additional computational costs while reducing the model complexity.

The channel flow simulations were initialized with a flow field from precursor simulations interpolated and scaled to match the friction Reynolds number. The simulations were run for a total of 15 eddy-turnover times ($T_{ϵ}=H/u_{*}$) with an initial transient of $10T_{ϵ}$. Time-averaged statistics discussed in this paper were averaged over the last $5T_{ϵ}$ unless otherwise specified. The flow is said to be statistically converged when the total stress profile above the roughness elements follows the linear stress profile as discussed in Patil and Fringer (2022). The full-span channels were run on 256 central processing units (CPU) and required about 276,500 CPU hours to simulate a total of 15 eddy-turnover times for both values of the Reynolds numbers. The minimal-span channels were run on 32 CPUs and required 7,700 CPU hours to simulate a total of 15 eddy-turnover times, reflecting savings in computational cost by a factor of 36 when compared to the full-span simulations.

Changing the parameters and their impact on the mean flow drag can be inferred by observing the location of the log-law region in the time-averaged and planform-averaged velocity profiles. Fig. 3(a) shows a comparison of the time-averaged and planform-averaged streamwise velocity for the full-span channel flow cases. Comparing the velocity profiles to the canonical flat-wall channel case (CF), the presence of roughness decreased the mean flow $U$ and thus increased the bottom drag coefficient for the bumpy wall cases. The bumpy wall log-law is given by Eq. (2), where Nikuradse (1933) found that $z_0={\overline{k}}_s/30$ for sand-grain type roughness (i.e., $C_o=1.0$). As shown in Fig. 3, the red dashed line corresponds to the log-law estimate given by Eq. (2) and the Nikuradse (1933) estimate for $z_0$. Case C350C1 exactly matched this prediction, and, more importantly, $z_0$ was not regressed, unlike in the other cases (i.e., $C_o=0.6$). The mean flow drag for case C700C1 was larger when compared to case C350C1 because the drag increased with increasing ${\mathrm{Re}}_{*}$. A similar observation can be made when case C700C06 is compared to case C350C06. The full-span channel results suggest that there was a consistent increase in the mean flow drag when decreasing $C_o$ for the two ${\mathrm{Re}}_{*}$ considered in this study.

To further understand the effect of $C_o$ on the mean flow drag, we validated the use of minimal-span channels (Jiménez and Moin 1991; Chung et al. 2015; MacDonald et al. 2017). As seen in Fig. 3(a), as the spanwise domain was restricted to incorporate the interaction of just two streamwise streaks, the velocity profiles in the minimal-span channels matched the full-span counterparts for $x_3^{+}\lesssim 160$, beyond which the profiles deviated from the log law. Therefore, because minimal-span channels accurately captured the near-wall velocity profiles, we used the minimal-span channels to extend the range of $C_o$ without imposing the large computational cost that would be required for the full-span channel cases.

In addition to the mean velocity profiles, there were distinct changes observed in the root-mean-squared (RMS) velocity profiles, as shown in Fig. 3(b). First, for cases with $C_o=0.6$, there was a strong decrease in the RMS velocity components when compared to the case with $C_o=1.0$ in the near-wall region. Further away from the wall, cases with $C_o=0.6$ were identical to those with $C_o=1.0$, further confirming the decoupled nature of the near-wall and outer regions of the flow (Townsend 1976). Similar mean flow response for the two Reynolds numbers suggested that the effect of $C_o$ was to modify the wall boundary condition such that decreasing values of $C_o$ resulted in a larger effective $z_0$. As for the Reynolds and viscous stress profiles, most of the variations occurred in the vicinity of the roughness elements. For the viscous stress, there was a prominent decrease in the maximum value with decreasing $C_o$. Additionally, as the friction Reynolds number increased, the relative contribution of the viscous stress decreased. These changes in the viscous stress profiles support the hypothesis that the stress contribution due to flow separation is expected to increase at the expense of the viscous stress for decreasing $C_o$ (discussed later). The Reynolds stress profiles, on the other hand, were independent of $C_o$ and followed the linear stress (blue dashed line) profile as expected.

To compare the drag coefficient for the full-span and minimal-span channels, we define the drag coefficientwhere $U_r$ = time-averaged and planform-averaged streamwise velocity evaluated at $x_3^{+}=120$. This definition of the drag coefficient allows for comparison of the full-span and minimal-span channels for cases where the full log-law velocity profile may not be available. Such a definition of the drag coefficient is quite common, especially in field experiments where only point measurements may be available (Egan et al. 2019). Fig. 4 shows a comparison of the drag coefficient for all cases detailed in Table 1. The overall trend was that with decreasing values of $C_o$, there was an increase in the drag coefficient for the two friction Reynolds numbers considered. Additionally, the drag coefficient predicted using the minimal-span channels was consistent with the full-span channels, providing further impetus to use the concept of minimal-span channels to predict the mean flow drag (Chung et al. 2015; MacDonald et al. 2017). For ${\mathrm{Re}}_{*}=350$, the drag coefficient increased as $C_o$ decreased until about $C_o=0.6$, beyond which $C_d^r$ saturated and started to decrease slightly. However, for ${\mathrm{Re}}_{*}=700$, a monotonic increase in the drag coefficient was observed with decreasing $C_o$. Because the estimate of $C_d^r$ is sensitive to $z_{\mathrm{ref}}$, $C_d^r$ seemed to saturate for ${\mathrm{Re}}_{*}=350$ but increased monotonically for ${\mathrm{Re}}_{*}=700$. Although it is unclear why case MC350C04 appears to be an outlier in the overall trend, it is omitted in the following analysis. The value of $C_d^r$ was more sensitive to $C_o$ for higher ${\mathrm{Re}}_{*}$, which is a result of the higher contribution of the form drag at higher ${\mathrm{Re}}_{*}$, as discussed later. It is expected that for increasing values of ${\mathrm{Re}}_{*}$, the drag coefficient becomes independent of ${\mathrm{Re}}_{*}$, similar to the behavior of pipe friction in the Moody chart. The independent behavior of the transition to ${\mathrm{Re}}_{*}$ is expected to occur at lower ${\mathrm{Re}}_{*}$ for increasing relative roughness such that once the flow is sufficiently hydraulically rough (i.e., ${\mathrm{Re}}_k\gg 70$), the drag coefficient will asymptote to a constant value. Because of the relatively low values of ${\mathrm{Re}}_{*}$ in this paper, due to computational constraints, there was a nontrivial dependence of the mean flow drag on ${\mathrm{Re}}_{*}$ for identical values of the Corey shape factor. As a result, for identical values of $C_o$ and ${\overline{k}}_s$, the mean flow drag changed nontrivially.

Another way to understand the effect of changing the geometry of the bumps is through changes in the roughness parameter $z_0$, as shown in Fig. 5. For the full-span cases, $z_0$ can be regressed to best fit the log law [Eq. (2)] because the log-law region is well resolved. This is evident in Fig. 3, which shows that the minimal-span channels reproduced the full-span velocity profiles. However, for cases without companion full-span channels ($C_o=0.4$, 0.8), the lack of a significant log law region did not allow such a regression to compute $z_0$. Consequently, to enable consistent comparison between the full-span and minimal-span channels, $z_0^r$ was inferred from the log law at a reference height of $z_{\mathrm{ref}}=120\nu /u_{*}$, such thatwhere $C_D^r$ = drag coefficient defined in Eq. (9). Fig. 5(a) suggests a similar overall trend as observed from the drag coefficient for changing $C_o$ and ${\mathrm{Re}}_{*}$. Additionally, it is clear that for increasing values of $C_o$, ${\alpha }_k$ increased, suggesting that the roughness height $z_0^r$ decreased for increasing $C_o$.

Although these observations provide a consistent way to relate the roughness characteristics, it is often more practical to relate the roughness height $z_0^r$ to statistical properties of the bed height in addition to the mean roughness (bed) height. Indeed, the results in the paper show very clearly that the bottom drag varied significantly through changes in $C_o$ even though ${\overline{k}}_s$ was constant for all simulations. To relate $z_0$ to the statistical properties of the bed, we calculated the standard deviation of the bed heightwhere $k_s^i$ = height of each of the ellipsoids; and $N\approx 360$ is the number of ellipsoids ($N$ is approximate because the number of ellipsoids that can fit on the bottom wall is subject to the random orientation angles). We then regressed the standard deviation to the roughness height with $z_0^r=\chi k_s^{\sigma }$, where $\chi$ is the regression parameter. As shown in Fig. 5(b), good correlation can be observed between $k_s^{\sigma }$ and the roughness parameter; that is, $z_0^r$, where $\chi ({\mathrm{Re}}_{*}=700)=0.188$ and $\chi ({\mathrm{Re}}_{*}=350)=0.112$. For ${\mathrm{Re}}_{*}=350$, the data point for $C_o=0.4$ was not included in the regression because it appears to be an outlier. Although it is unclear why this point was an outlier, we expect this to be a consequence of the minimal-span nature of the channel. Using linear regression, we observed $R^2=0.891$ for ${\mathrm{Re}}_{*}=700$, and $R^2=0.765$ for ${\mathrm{Re}}_{*}=350$, suggesting a strong correlation between the roughness characteristics and the expected roughness parameter. In addition to the ${\mathrm{Re}}_{*}$ dependence, the definition of $C_o$ does not account for varying values of the sphericity ($S_p$) for identical $C_o$, as suggested by Julien (2010). It is clear to see that $z_0^r$ is sensitive to ${\mathrm{Re}}_{*}$, which can be inferred from the regression parameter. Some of the scatter observed in the data presented here can be attributed to the transitional roughness Reynolds number regime because the flow separation was localized to some roughness elements that penetrated beyond the viscous sublayer, as shown in Fig. 6 and discussed in Schultz and Myers (2003) and Flack et al. (2012). These observations may help explain the variability observed in conventional methods to estimate the mean flow drag (e.g., Moody 1944) or the roughness function ($\mathrm{\Delta }{U}^{+}$) defined in Schultz and Myers (2003).

Increased mean flow drag with a simultaneous decrease in the viscous stress for decreased $C_o$ suggests that there was a net increase in the form drag components as the flow separated at the crest of the roughness elements. As shown in Fig. 6, there was relatively more flow separation for cases C350C06 and C700C06 when compared to cases C350C1 and C700C1. These observations, along with the attenuated viscous stress profiles and increased drag coefficient, suggest that the form drag for smaller $C_o$ increased. In addition to Fig. 6, a vorticity comparison can be seen in Fig. 7, which shows instantaneous vorticity for the four full-span cases. Comparing cases with ${\mathrm{Re}}_{*}=350$ against ${\mathrm{Re}}_{*}=700$ clearly shows the effect of increased Reynolds number on the vorticity distribution close to the wall. The effect of $C_o$ was more pronounced for lower ${\mathrm{Re}}_{*}$, as indicated by the slightly more negative values for case C350C06 when compared to case C350C1. For higher ${\mathrm{Re}}_{*}$, there were no obvious differences in the vorticity between cases C700C1 and C700C06.

Because the immersed boundary force ($F_{\mathrm{IBM}}$) is imposed at every substep in the Runge–Kutta time-integration scheme, only the divergence-free velocity at the end of each time step is available (Yuan and Piomelli 2014, 2015). Therefore, although the IBM force can be used to compare the total drag force on the bed, it does not give the relative contributions of viscous and form drag. To separately compute these components of the drag, we began with the streamwise momentum equationwhere $P=p/{\rho }_0$ = modified or reduced pressure. Defining ${\mathcal{V}}_f$ as the volume occupied by the fluid above the roughness elements and integrating Eq. (12) over ${\mathcal{V}}_f$ gives, after using Gauss’s theorem and noting that $F_{\mathrm{IBM}}=0$ in ${\mathcal{V}}_f$ and imposing periodicity in the $x_1$ and $x_2$ directions and free-slip condition at $x_3=H$,where ${\mathcal{A}}_{\mathrm{B}}$ = control surface at the bumpy wall that corresponds to the top of the roughness elements; and $e_j$ = unit normal vector pointing outward relative to ${\mathcal{V}}_f$. After time-averaging Eq. (13), the unsteady term vanishes, thus giving the momentum partitioningwhere the overbar represents the time average. The advective term on the left-hand side is nonzero because the IBM method imposes a force that produces a vanishing cell-centered velocity component and a small nonzero face-centered component where ${\mathcal{A}}_{\mathrm{B}}$ is defined. As a result, the advective term does not vanish, although it is much smaller than the other terms. Table 2 lists the normalized contribution of the three terms in the mean momentum equations for the full-span channel flow cases. The results indicate an enhanced form drag component for $C_o=0.6$ that can be interpreted as a consequence of increased flow separation because the roughness elements are taller when compared to the case with $C_o=1.0$. Additionally, this increase in the form drag occurs with a simultaneous decrease in the viscous drag for the two Reynolds numbers. Because the driving pressure gradient (i.e., $u_{*}^2/H$) is constant for varying $C_o$ for a given ${\mathrm{Re}}_{*}$, only the bottom boundary conditions are responsible for an increase in the drag coefficient and the mean momentum partition. It is clear from Table 2 that there was a definitive increase in the form drag for decreasing $C_o$ for both ${\mathrm{Re}}_{*}$.

The relative importance of the viscous and form drag (relative drag) is depicted for varying $C_o$ in Fig. 8. The minimal span channels capture the overall trend as predicted by the full span channels, with a consistent overprediction for all cases when compared to the full span data. As shown in Fig. 3, the time-averaged and planform-averaged velocity profiles for the minimal-span channels led to a flow velocity that was comparatively large in magnitude away from the wall. Because the minimal-span channels did not effectively mix momentum in the vertical direction due to the spanwise domain constraint, we anticipated larger form drag as a result of this increased mean flow velocity away from the wall when compared to the full-span channel cases (Chung et al. 2015; Yuan and Piomelli 2015). For both values of ${\mathrm{Re}}_{*}$, the relative drag increased for decreasing values of $C_o$. Additionally, it is clear from these data that the relative drag increased more quickly for ${\mathrm{Re}}_{*}=700$ when compared to ${\mathrm{Re}}_{*}=350$, further validating the changes observed in $C_d^r$ (see Fig. 4).