Toward modeling of turbulent flow over surfaces of arbitrary roughness

The prediction of the effects of rough surfaces on the mean velocity field in turbulent flow is a long-standing challenge for turbulence modelers and has many practical applications. When contrasted with modeling of smooth flat-wall boundary layers, in which only the tangential viscous stress $μ∂u¯/∂y$ is exerted by the fluid on the wall, this shear stress is both modified by the flow and supplemented by the roughness-induced form drag force per unit area, caused by the flow-induced pressure distribution around each surface roughness element. The combination of these effects and of any blockage caused by roughness elements must, therefore, be modeled for rough-wall applications.

In most of the modeling approaches developed to date, the surface roughness is quantified by a single lengthscale within a closure relationship. An empirical calibration is then sought so that predictions made with the roughness model, together with any other required turbulence closures, approximate selected reference data—usually for sandgrain roughness—with some degree of fidelity. An early example of the single-lengthscale approach is Rotta's¹ proposal of $Δy$ as an additive shift to the surface-normal coordinate employed in models (⁠$y→y+Δy$⁠) to accommodate the effects of surface roughness, with the lengthscale $Δy$ modeled as an empirical function of the equivalent sandgrain roughness size k_s. Similarly, in their one-equation eddy-viscosity transport model, Aupoix and Spalart² modified the smooth-wall model of Spalart and Allmaras³ by introducing a roughness lengthscale as a wall-offset distance to increase the near-wall eddy viscosity in order to represent the effect of surface roughness.

For the two-equation k–$ϵ$ closure, Durbin et al.⁴ introduced a two-layer model, which described the effects of roughness through the additive shift of the wall-normal coordinate for turbulent quantities, which gave good predictions of flows over ramps, sand dunes, and flat walls at different pressure gradients. In a similar approach, Chen and Patel⁵ used a one-equation closure for k and a wall-function model for $ϵ$⁠, together with separate lengthscales in the near-wall eddy viscosity model and the $ϵ$ function, to incorporate the effect of surface roughness. Brereton and Yuan⁶ introduced an additive wall-roughness eddy viscosity across the roughness sublayer, which scaled on equivalent sandgrain size k_s and $Reτ$⁠, to represent surface roughness while, in a related wall-function approach, Suga et al.⁷ modeled roughness with the k–$ϵ$ closure using analytical functions to describe the near-surface velocity, dissipation, and production of k, in which a viscous sublayer thickness was scaled on k_s to represent effects of wall roughness. In the k–ω turbulence model,⁸ the effect of surface roughness was modeled by generalizing the wall-boundary value of the timescale ω from its single smooth-wall value to one determined from an empirically derived function of the equivalent sandgrain size k_s, to model rough-wall mean velocity profiles and friction coefficients. Yang et al.⁹ proposed a scheme based on the von Kármán–Pohlhausen polynomial velocity-profile method, in which the mean velocity in the roughness sublayer was modeled as an analytical function which grew exponentially with surface-normal distance, damped by a single prescribed attenuation coefficient which carried information on the surface roughness. In an analysis of measured velocity fields within roughness sublayers, Nikora et al.¹⁰ proposed physical scaling arguments for sublayer velocity profiles which were (i) uniform; (ii) linearly increasing with y; and (iii) exponentially increasing with y, each of which were consistent with different experimental data. At a more fundamental level, Forooghi et al.¹¹ have shown how the effect of sandgrain surface roughness can be represented in direct numerical simulations over smooth walls, by adding to the x-momentum equation a distributed body force modeled as two terms which depend linearly and quadratically on the local velocity, and on the y-dependence of the roughness geometry.

In most of these modeling approaches, roughness has been represented primarily by a single parameter, typically correlated empirically to yield a downward shift in a log-linear mean velocity profile or an increase in friction coefficient as a function of either the roughness height (for sandgrain roughness) or an equivalent sandgrain roughness height for surfaces of other textures. The main shortcoming of these single lengthscale models is that they do not provide a way of making an accurate estimate of a sandgrain roughness height from the true shape of a given rough surface. This point is of considerable importance because, for many realistic rough surfaces, the equivalent sandgrain roughness height k_s may be several times larger than the height of the roughness crest k_c and so an incorrect estimate of the equivalent sandgrain size can lead to a significant underestimation of the surface drag. Another shortcoming of the sandgrain roughness characterization is that the precise structure of the grains used by Nikuradse¹² in the experiments which formed the basis for this approach was not documented, though Thakkar et al.¹³ have replicated Nikuradse's results with direct numerical simulations of channel flow over a grit-blasted surface.

Flack and Schultz¹⁴ have examined proposed correlations between k_s and other measures of roughness topography and observed that (i) surfaces with an “effective slope” of less than 0.35 were “wavy” rather than rough, in which case k_s did not scale on roughness topography at all; and (ii) good correlation could be achieved with a simple algebraic function for k_s in terms of rms surface height and skewness, for seven different surfaces of predominantly positive skewness. While other correlations for k_s have been proposed by Musker¹⁵ and Medhurst,¹⁶ the development of a reliable correlation between general roughness topography and k_s over a wide range of roughness types and magnitudes is far from complete. Recently, Aghaei Jouybari et al.¹⁷ have presented a neural network model for determining k_s in the fully rough regime from a database of flows over many different kinds of surface roughness. They also showed that the downward shift in the log-linear mean velocity regime on account of roughness, observed in many sandgrain roughness experiments at high Reynolds numbers, is also found for all other kinds of roughness they explored at low Reynolds numbers, when written as $u+=(1/κ) ln (y+−d+)+B$⁠, where d⁺ is the zero plane displacement height.

More direct ways of characterizing the effect of roughness on turbulent flow are now possible through direct numerical simulations of turbulent flow over rough surfaces, with no-slip/no-penetration boundary conditions enforced at fine-grained levels using immersed-boundary methods. With these techniques, detailed data on rough-wall boundary layers can be determined which resolve the instantaneous flow from the far field to within the roughness sublayer (Yuan and Piomelli,¹⁸ Busse et al.,¹⁹ and Chan et al.²⁰). These computational advances provide the capability to develop databases of turbulent flow over surfaces of arbitrary roughness size/texture, under the current restriction of the Reynolds number $Reτ$ less than a few thousand, which can then be used to construct surface-specific flow models.

In this paper, we present an approach to RANS (Reynolds-averaged Navier–Stokes) modeling of rough-wall boundary-layer flows, which is based on the effect of roughness on terms in the x-momentum equation. From analyses of the distributions of drag forces, and turbulent and dispersive stress gradients, it is shown that their complexity over different rough surfaces is such that the most practical technique for modeling flows over surfaces of arbitrary roughness is to model the sublayer velocity directly as a roughness-specific velocity profile or wall function, and to model the turbulent stresses beyond the roughness sublayer as their smooth-wall counterparts. Each kind of roughness is then described as a particular near-wall velocity profile, in wall units, which is extracted from a database of surrogate channel-flow direct numerical simulation results. In this paper, we examine the effectiveness of this approach as a means of modeling turbulent flow over surfaces with very different roughnesses.

The turbulent boundary layer in flows over rough walls may be described as an outer region above a roughness sublayer which extends from the deepest troughs of roughness elements to just above their peaks. When the fluid velocity and pressure are decomposed into (i) the time-space average (the spatial $⟨⟩$ average of a time $¯$ average); (ii) the “dispersive” spatial variation about the time-space average $̃$⁠; and (iii) the time-unsteady fluctuating component $′$⁠, as $ui(x¯,t)=⟨u¯i⟩(y)+ũi(x¯)+ui′(x¯,t)$⁠, the double-averaged Navier–Stokes equations (Raupach and Shaw²¹ and Nikora et al.²²) may be written as

Here, $⟨ ⟩s$ signifies superficial spatial area averaging, where the y-dependent fluid variable is averaged per unit total planar area, rather than per unit fluid planar area. It is assumed that spatial averaging is carried out over a thin “slab” of sufficient extent in the x–z plane to smooth out spatial heterogeneities of the surface and sufficiently thin to resolve the behavior of flow variables in the y direction (Finnegan²³). In Eq. (1), $f¯v,i$ and $f¯p,i$ are pseudo body forces per unit fluid mass, the streamwise components of which are the viscous and pressure drag forces. They are defined using integrals around the surface S of roughness elements in a given averaging plane as

where V₀ is the slab volume.

When an immersed boundary method is used in simulations to enforce a no-slip condition locally at the fine-grained fluid/solid boundary of a rough surface, the incompressible, constant-property, linear momentum equations are written as

where F_i is the local, instantaneous body force required to enforce the no-slip boundary condition. When area-averaged over all boundary cells in a given x–z plane, the x-component of the immersed body force F_i at a given y location is equal to the total drag force (Yuan and Piomelli²⁴), i.e.,

Thus, the sum $f¯v,x+f¯p,x$ in Eq. (1), which we rewrite as $⟨fx⟩s$⁠, can be found readily from immersed boundary method simulations of rough-wall flow. The distribution of this pseudo body force across the roughness sublayer—the region extending from the deepest roughness troughs at y = 0 to just beyond the roughness crests at y = k_c—is shown in Fig. 1, for fully developed flow in a channel with walls of synthesized sandgrain roughnesses of two different sizes at $Reτ=1000$ (Yuan and Piomelli¹⁸), where $Reτ$ is the Reynolds number referenced to friction velocity and channel half-width h.

Conventional turbulence modeling is concerned with developing single-point representations of the unclosed terms in the Reynolds-averaged Navier–Stokes equations, for solution on a grid that is usually very coarse relative to the fine-scale topography of a rough surface. It is, therefore, desirable to incorporate the roughness geometry as part of the rough-wall turbulence closure, so the model equations can be solved on a smooth-wall computational grid. This effect is achieved by describing flow variables at locations within the roughness sublayer as their Reynolds averages based on the superficial spatial area average $⟨ ⟩s$ defined earlier. The time-space averaged momentum equations with a body-force roughness representation and superficial planar-area averaging for solution on a smooth-wall grid are then

To close the streamwise component of this equation in boundary-layer applications, models might be sought for the y-distribution of the roughness body-force $⟨fx⟩s$ and the superficially averaged dispersive and Reynolds stresses $⟨ũṽ+u′v′¯⟩s$ or their gradients, either separately or in some combination.

The roughness body force is an important descriptor of the effects of roughness on flow and provides a useful reference scale for them. It is convenient to express this pseudo body force f_i in viscous units as $fi+=fiν/uτ3$⁠, where $uτ=τ0/ρ$ and is the friction velocity. For the case of fully developed channel flow, the streamwise component of Eq. (5) may be integrated with respect to y across a channel of half-width h, with y = 0 at the lower wall, under the assumption of y-invariant $∂p¯/∂x$⁠, to yield the stress balance in viscous units,

At $y+=0$⁠, the lowest point in the roughness sublayer, all components of velocity are zero and Eq. (6) becomes

If we note that $⟨fx⟩s→0$ at the edge of the roughness sublayer, just beyond the roughness crest where y = k_c, we can write this expression for the stress at y = 0 as

Thus, the shear stress at the lowest point in the roughness sublayer can be expressed as the sum of the roughness body-force integral over the sublayer and the corresponding “smooth-wall” viscous term—an interpretation which is important for modeling rough-wall flow. At sufficiently large Reynolds numbers, the contributions of both the explicit viscous term and the viscous part of $⟨fx⟩s$ become negligible and the flow is termed “fully rough.”

The simulations of Yuan and Piomelli,¹⁸ for flow over surfaces which are rough in their entirety, indicate that the viscous stress $∂⟨u¯⟩s+/∂y+=0$ at $y+=0$⁠. It follows that the integral of the roughness body force for such surfaces is constrained as

and that after multiplication by $kc+$ the roughness-related terms in Eq. (6), or Eq. (5) in wall units, are of order unity. This rescaling provides a rational way of comparing roughness-related quantities from flows over different, entirely rough surfaces.

Profiles of the roughness body force $⟨fx+⟩s$ across the roughness sublayer are plotted in Fig. 1 for fully developed channel flow at $Reτ=1000$⁠, over walls of synthesized sandgrain roughness of heights $kc+=28$ and $kc+=86$⁠, or $kc/h=0.028$ and $kc/h=0.086$⁠, respectively. The abscissa scaling chosen for this figure is $kc+⟨fx+⟩s$⁠, consistent with the integral constraint of Eq. (9). The surface with the smaller roughness had an equivalent sandgrain roughness of $ks+=22$ and so was considered to be in the “transitionally rough” region, whereas the rougher surface had an equivalent sandgrain value of $ks+=72$⁠, which is in the “fully rough” and so Reynolds-number-independent regime (Yuan and Piomelli¹⁸).

In Reynolds-averaged turbulence closures, it is conventional to model the turbulent stresses rather than their gradients. It is, therefore, important to assess the relative sizes of the turbulent and dispersive stresses, their distribution within the flow, and their similarity to smooth-wall counterparts, for which the modeling literature is extensive. Profiles of Reynolds and dispersive stresses from the channel-flow simulations of Yuan and Piomelli¹⁸ over sandgrain roughnesses at $Reτ=1000$ are plotted in Fig. 2 in viscous units, denoted by the $+$ superscript, for the two roughnesses shown in Fig. 1. The corresponding profile for $u′v′¯+$ in a smooth-walled channel at the same Reynolds number is also plotted for the purpose of comparison. Beyond each roughness sublayer (⁠$y/h>0.026$ and 0.086), the dispersive stress $⟨ũṽ⟩s+$ is significantly smaller than the rough-wall Reynolds stress $⟨u′v′¯⟩s+$⁠, although it can be seen by inspection that the sum of the rough-wall stresses $⟨ũṽ⟩s++⟨u′v′¯⟩s+$ is almost equal to the smooth-wall Reynolds stress $u′v′¯+$ for the two roughnesses shown. This comparability between outer-flow stress profiles in smooth- and rough-wall flow is consistent with the outer-layer similarity hypothesis of Townsend,²⁵ The effect of roughness elements within the sublayer is to interfere with the near-wall turbulence regeneration processes that create $⟨u′v′¯⟩s$ in such a way that the region of growth and the location of the peak are shifted outward, relative to the Reynolds stress distribution above a smooth wall.

One way to close the x-momentum equation of Eq. (5) is to seek a model- or simulation-based representation of the distribution of the roughness body force $⟨f¯x⟩s$ across the roughness sublayer. Since the stress gradients within this sublayer must also be modeled, it is useful to compare $⟨f¯x⟩s$ with the y-derivatives of $⟨u′v′¯⟩s$ and $⟨ũṽ⟩s$ in this region to develop a sublayer modeling strategy. Profiles of $⟨f¯x+⟩s,−∂⟨u′v′¯⟩s+/∂y+$⁠, and $−∂⟨ũṽ⟩s+/∂y+$ are plotted in Fig. 3 for the case of sandgrain roughness of height $kc+=86$⁠, with each term multiplied by $kc+$ for ease of comparison with the integral of the scaled roughness body force, which is unity [Eq. (9)]. The sum of these terms must be balanced by the axial pressure gradient and the viscous stress gradient according to Eq. (5). It is clear from Fig. 3 that while the roughness body force $⟨f¯x⟩s$ is the largest of these three terms, the y-derivatives of $⟨u′v′¯⟩s$ and $⟨ũṽ⟩s$ are of the same order of magnitude within the roughness sublayer. These gradients of Reynolds and dispersive stresses are also much larger within the roughness sublayer than in the outer flow. The behavior of stress gradients for the corresponding flow over sandgrain roughness of $kc+=28$ was very similar. The relative sizes of the y-gradients of $⟨u′v′¯⟩s$ and $⟨ũṽ⟩s$ in Fig. 3 indicate that these terms cannot be ignored within the roughness sublayer, but must also be included in a rough-wall closure. The prospect of modeling each of them in terms of local flow variables and surface topography for arbitrary rough surfaces seems quite remote.

An additional modeling problem becomes apparent when the form of the DNS-resolved mean velocity $⟨u¯⟩s+$ within the roughness sublayer is considered. It is plotted in Fig. 4 for flow at $Reτ=1000$ over sandgrain roughnesses of two different heights. The mean velocity is almost zero over the 20% of the sublayer closest to the wall for both roughnesses, and almost linear over the outer 50% for $kc/h=0.028$⁠. When the velocity profile is linear in y, the viscous diffusion term in Eq. (5) is zero and the modeled body force must balance the pressure gradient and the modeled turbulent and dispersive stress gradients in a fully developed flow. Where viscous diffusion is not zero, the small difference between these large modeled terms is required to take a particular form to yield the correct y-behavior in $⟨u¯⟩s+$⁠. The difficulties inherent in such a modeling approach are such that it appears to be more prudent to employ directly a roughness-sublayer model for $⟨u¯⟩s+$ that is surface specific, as deduced from DNS data, thereby removing the need to solve the x-momentum equation in the roughness sublayer altogether. For these reasons, this approach is explored and developed in Secs. VII and VIII of this paper.

Two strategies for modeling $⟨u¯⟩s+$ (henceforth referred to as $u+$ in modeling contexts) across the sublayer for particular rough surfaces are presented: a coarse-mesh wall-function model intended to predict the correct wall shear stress, without resolving the mean velocity in the sublayer; and a fine-mesh inner velocity-profile model to resolve both the wall shear stress and the complete sublayer velocity profile. In both strategies, the roughness-specific near-wall model is coupled with a standard k-$ϵ$ outer-flow closure

In the conventional approach to wall-function modeling, an algebraic description of $u+(y+)$ is provided as a function called by a RANS equation solver. Its purpose is to provide the velocity in wall units $up+$ at the first nodal point beyond the wall $yp+$⁠, in a standard finite-volume discretization. The mesh size must be such that this point is within an assumed log-linear region in which the velocity is described by $u+=(1/κ) ln y++B$ or $u+=(1/κ) ln Ey+$⁠, where $E=eκB$ and κ is the von Kármán constant. The effective viscosity ν_w to be applied at the wall is chosen to equate the wall shear stress there, where the velocity is modeled as linear, to that at y_p in the log-linear region, i.e.,

so that

Thus, the wall shear stress is made to take the value required by the assumed form of the log-linear velocity at y_p. Consequently, the requirement that y_p is no closer to the wall than the inner extent of the log-linear region places a lower limit on the mesh size, or an upper limit on the number of uniformly distributed cells across a domain. If the near-wall velocity field (in wall units) deduced from one flow is universally applicable to other near-wall flows over the same roughness (in wall units), one needs only to determine the value of E and the inner extent of the log-linear region from simulation data to construct such a wall function for a particular rough surface. With a suitable mesh and a model relationship to specify k⁺ and $ϵ+$ at $yp+$ as additional boundary conditions, a coarse-mesh RANS calculation can then be made for flow over the rough surface of interest.

In DNS of fully developed flow over rough surfaces, the value of τ_w and thus of $uτ$ is determined as that required to balance the pressure gradient. When this value of $uτ$ is used to attempt a log-linear fit over part of the mean velocity profile, an offset to the $y+$ coordinate of d⁺—the zero-plane displacement defined as $d=∫yFxdy/∫Fxdy$ and rescaled in local viscous units—is necessary to achieve a uniform slope of $1/κ$⁠. Therefore, any wall function intended to describe a log-linear region in flow over a rough wall must also include such an offset if it is to yield the desired values of $uτ$ and τ_w, and so is written as $u+=(1/κ) ln (y+−d+)+B=(1/κ) ln (E(y+−d+))$⁠. Either d⁺ must be included in the wall function, or the mesh node at y = 0 is set to y = d. The inclusion of this offset is important for flows at low Reynolds numbers, but much less so at high ones. While this coarse-mesh approach to rough-wall modeling enforces the wall-function-provided value of τ_w, it can place significant limits on the spatial resolution of the solution, it does not represent the average velocity across the cell adjacent to the wall correctly, and it provides no information on the velocity profile within the roughness sublayer, which would be essential to treating coupled heat- and mass-transfer problems. This wall-function modeling approach is evaluated in fully developed channel flows with walls comprising four very different kinds of roughness, in Sec. VIII A.

The second proposed approach is to specify the profile of $u+$ within the roughness sublayer and combine it with a Reynolds-averaged closure over the outer part of the flow. The sublayer profile of $u+$ is selected as a roughness-specific function from a library of simulation results for fully rough flow over surfaces of many different kinds of roughness. Such surface-specific sublayer profiles are specified as $u+(y+)$ and typically require a companion determination of $uτ$ at each step in a numerical solution to determine u(y). Since the mean velocity profile is effectively a scaled library function, the turbulence model equations can be solved within the roughness sublayer by imposing the sublayer velocity field and applying standard wall conditions at $y+=0$⁠. Alternatively, one can impose the sublayer velocity field and solve only beyond the roughness sublayer, in a manner similar to that of the wall-function method, but without the lower limit on mesh size. In that case, boundary values for, say, k and $ϵ$ are applied at y = k_c as those values at $y+$ in the corresponding smooth-wall flow, consistent with the observations of Yuan and Piomelli,¹⁸ or from standard model relationships.

The outer flow is modeled by the Reynolds-averaged linear momentum equations and a standard two-equation closure from roughly $y+=kc+$ outward, with far-field boundary conditions appropriate to the particular flow of interest. The value of $uτ$⁠, which couples the solutions in the two regions, is determined by balancing the wall shear stress τ₀ with the pressure gradient in fully developed duct flows, or from the momentum integral equation in a developing boundary layer. Whereas the inferred value of $uτ$ enables u_p to be determined from Eq. (10) in the wall-function model, it allows the local sublayer u(y) profile to be found from the given $u+(y+)$ profile in this model. Furthermore, there is no restriction on the number or size of cells in the roughness sublayer or beyond it.

Roughness-sublayer velocity profiles were determined from the library of 42 fully rough flows over different surfaces, compiled by Aghaei Jouybari et al.¹⁷ These flows were computed using direct numerical simulation of a half-channel flow with an upper “slip wall,” at $Reτ=1000$⁠, and appear to be in very good agreement with more expensive full-channel simulations as far from the wall as $y/h=0.4$⁠, which is well beyond the roughness sublayer for all these flows. Although a discrete set of velocity-profile data can be extracted from the library for each kind of roughness, it is more convenient to represent each sublayer profile with a composite model equation in which only a small number of surface-dependent parameters are specified.

In this composite profile model, illustrated in Fig. 5 and described in detail in the Appendix, a region of zero or negligible mean velocity begins at the deepest trough within the roughness elements, at $y+=0$⁠, and is blended to an outer log-linear region with a cubic polynomial. By introducing the variable $y∗=y+−d+$⁠, and identifying $y0∗$ and $yl∗$ (or $y0+$ and $yl+$⁠) as the point at which $u+$ first exceeds 0 and the innermost extent of the log-linear region, the composite model for $u+$ in rough-wall boundary layers may be specified as:

where $f(y∗)$ is the polynomial and the parameters d, B, y₀, and $yl∗$ are provided from the library/database for the chosen rough surface. For many rough-wall flows, the roughness sublayer does not reach $yl∗$ and so only the $u+=0$ region and a part of the $u+=f(y∗)$ section are needed.

Preliminary tests of this modeling approach were made against the simulation data of Yuan and Piomelli¹⁸ for flow in channels with sandgrain roughnesses of two different crest heights $kc+=28$ and 86, at $Reτ=1000$⁠. Of these test cases, the former with the smaller roughness was in the transitionally rough regime, while the latter was considered fully rough. The coefficients of the composite profile model were first determined from results of each of the simulations, and the sublayer velocity-profile models applied from $y+=0$ to $y+=kc+$⁠. The outer region of the two-layer model was described by a k–$ϵ$ closure and the combined model solved with finite volume discretization and a pseudo-unsteady Crank–Nicolson scheme, with iteration on the value of $uτ$ until convergence was reached at $Reτ=1000$⁠. The computation was carried out on a half-channel domain, with symmetry boundary conditions on $u¯$⁠, k, and $ϵ$ at the channel centerline. The cell spacing was chosen to be $Δy+=1$ at $y+=0$⁠, with each successive cell increasing in size by 8%.

The results of these test computations are shown in Fig. 6, together with a computed smooth-wall profile as a reference. The excellent agreement between the model and DNS results across the roughness sublayer indicates only that the composite model describes accurately the simulated velocity profile, with correctly chosen coefficients. The smooth blending of the outer profile with its sublayer counterpart and the general agreement over the log-linear region implies that the outer-region k–$ϵ$ model and its inner matching with the prescribed sublayer profile is satisfactory for this application. It is interesting to note that in neither case is the edge of the roughness sublayer within the log-linear region. Therefore, $kc+<yl+$ and so the outer-flow k-$ϵ$ model is used to determine $u+$ in the buffer layer from $kc+$ to $yl+$⁠, and beyond. The model friction coefficients (defined as $τ0/(ρU2/2))$ of $Cf=0.007 1$ and 0.014 0 are about 9% lower than the simulation values of $Cf=0.007 8$ and 0.015 3, primarily as a consequence of predicting a more pronounced wake than the simulations. Slight improvements in the prediction of C_f were achieved when the sublayer profile model was applied as far from the wall as $1.1kc+$ and $1.2kc+$⁠. However, they provide less of a test for this modeling approach when the profile model is applied as far from the wall as $yl+$⁠, since the values of the k-$ϵ$ model coefficients are chosen to ensure consistency with some features of a log-linear velocity profile.

To test these roughness sublayer modeling approaches, we consider as target data fully rough turbulent channel-flow simulations at $Reτ=1000$ of flow over the following rough surfaces: (i) sandgrain roughness with $kc+=90$⁠; (ii) a uniformly spaced array of half-ellipsoids of height $kc+=90$ on a smooth surface; (iii) a two-dimensional sinusoidally wavy surface of height $kc+=60$⁠; and (iv) random (white noise) low-order spatial Fourier modes of height $kc+=120$⁠. These surfaces correspond to cases C37, C17, C30, and C40, respectively, in the DNS database of 42 rough-surface half-channel flows compiled by Aghaei Jouybari et al.¹⁷ These particular surfaces—each with large roughness heights—were chosen because their different topographies result in quite different sublayer velocity profiles, and the prospect of predicting them by other modeling approaches appears to be quite remote.

To evaluate the wall-function modeling approach, k-$ϵ$ calculations of fully developed turbulent channel flow at $Reτ=1000$ were carried out, with a user-defined wall function $u+=(1/κ) ln (E(y+−d+))$ in the finite-volume code Fluent, in which E and d⁺ were taken from the log-linear region of the mean velocity in DNS of channel flow over the same roughness. In these calculations, the area-averaged velocity U was varied systematically until the target value of $Reτ$ was reached, after which the modeled friction coefficient C_f was calculated and compared with the corresponding DNS value. This test was performed for each of the four very different rough surfaces with large roughness crest heights k_c shown in Figs. 7–10 and described in more detail in Subsection VIII B. The results are shown in Table I, together with the values of B (or E) and d⁺ provided for each roughness from the DNS database. For each calculation, 12 cells of uniform width $Δy+≃80$ were distributed between the wall and the channel center, with $yp+≃40$⁠. In each case, the predicted friction coefficient was between 3% and 6% lower than the DNS value, notwithstanding the use of so coarse a mesh. In the context of this test, these results imply that the coarse-mesh k-$ϵ$ outer-flow calculation beyond the specified wall-function velocity at y_p still yields an accurate representation of the computed bulk flow rate, despite the misrepresentation of the average velocity in the wall-adjacent cell.

The first surface considered was a large sandgrain roughness with $kc+=90$ and is shown together with its near-wall velocity profile in Fig. 7. It was chosen as a reference against which sublayer velocity profiles for other surfaces might be compared because of the familiarity of this surface to researchers. As a surface that is rough in its entirety, $d⟨u¯⟩s+/dy+=0$ at $y+=0$ and the superficially averaged velocity $⟨u¯⟩s+$ first exceeds 0 at $y+≃25$⁠. The parameters B, d⁺, $y0+$⁠, and $yl+$ provided from the DNS database for the profile model equation for this surface are given in Table II and are consistent with those in a generalized velocity-profile model for sandgrain roughness proposed in Sec. IX B. The friction coefficient computed with this sublayer velocity-profile model and a standard k-$ϵ$ model in the outer flow was 0.016 0 and matches, within modeling error, the value of 0.015 6 obtained from the DNS computation.

The second test surface was a smooth wall upon which resided a uniformly spaced array of semi-ellipsoids of height $kc+=90$⁠, shown together with its near-wall velocity profile in Fig. 8. The roughness elements cover only part of the otherwise-smooth surface and $⟨u¯⟩s+$ appears to grow almost linearly from its no-slip condition at $y+=0$⁠, in which case $y0+=0$⁠. The values of the other parameters which fitted the model velocity profile to the simulation result are given in Table II. It is particularly difficult to propose rational models for turbulent flow over this kind of surface, which combines elements of smooth- and rough-wall flow. Although of the same roughness crest height as the previous surface, the downward shift of the log-linear region is greater and the friction coefficient of 0.021 6 is 40% higher. The corresponding sublayer velocity-profile model prediction was 0.022 8, which is about 5% larger than the half-channel DNS value 0f 0.021 6.

The third test surface modeled was a two-dimensional sinusoidally wavy wall with a peak-to-peak amplitude of $kc+=60$⁠, shown in Fig. 9. Short-wavelength walls of this kind generate some of the flow features of other kinds of roughness and the ability to describe the surface analytically makes them useful for developing and testing rough-surface-flow codes (e.g., Yuan et al.²⁶). An interesting feature of this case is that the superficially averaged streamwise velocity is slightly negative over about half the roughness sublayer, with very steep surface-normal gradients in $⟨u¯⟩s+$ near the roughness crest. The sublayer velocity-profile model, with a sublayer velocity profile which approximated the reverse-flow region as one of zero velocity, predicted $Cf=0.012 9$—slightly higher than the simulation result of 0.012 8.

The final test surface considered is shown in Fig. 10 and was constructed as low-order x- and z-direction Fourier modes of white noise with height $kc+=120$⁠. It was chosen as a randomly textured surface to contrast with the first surface, for which each ellipsoidal element had the same orientation and semi-axis length to mimic a uniform sandgrain roughness. For this surface, the sublayer velocity-profile model prediction of C_f was 0.012 1, which is slightly higher than the simulation result of 0.011 6. It is interesting to note that the effect of randomness in surface texture and a roughness crest of $kc+=120$ yields a value of C_f that is 25% lower than that of the regular sandgrain roughness with $kc+=90$ and emphasizes the importance of orderliness and other subtle effects of texture on the resulting surface drag.

In summary, each sublayer velocity-profile model prediction was between 1% and 6% higher than the corresponding half-channel DNS result. The uncertainty that arises in fitting the coefficients of the model to the roughness sublayer profile is thought to lead to errors of a few percent in C_f , to which other modeling and numerical uncertainties are added. It is also worth pointing out that the slip-wall upper boundary condition used in half-channel DNS has the effect of increasing slightly the outer-flow velocity, while reducing that of the inner flow, so that a half-channel DNS tends to predict values of C_f that can be up to 10% lower than those from full-channel DNS (MacDonald et al.²⁷). This observation may explain why all the sublayer-profile-model C_f predictions, with symmetry boundary conditions applied to $u¯$⁠, k, and $ϵ$ at the channel centerline, are consistently slightly higher than the half-channel DNS results. It bears an interesting contrast with the wall-function results of Subsection VIII A, for which values of C_f were 3%–6% lower than the DNS results, because of mesh coarseness and the overestimation of $u¯$ in the cell adjacent to the wall.

In the case of turbulent, fully developed pipe flow, Colebrook's²⁸ formula can be used as an algebraic model of target data against which the wall shear stresses computed by the wall function and sublayer velocity-profile models can be tested, over a wide range of roughness heights and flow Reynolds numbers.

For the case of the wall-function model of flow in pipes with sandgrain roughness, only the quantities d and B need be specified. Since this model is used for coarse-mesh flow computations, little additional error is usually introduced if it is also assumed that the log-linear region of the velocity extends from $y+=d+$ to the pipe center, since wakes in fully developed pipe flows are small. The velocity-profile equation $u+=(1/κ) ln (y+−d+)+B$ can be integrated from r = 0 to $r=h−d$ (or y = d to y = h), where h is the pipe radius, to yield the mean or area-averaged velocity,

The Colebrook formula for the friction factor in pipe flow with rough walls of average roughness height ϵ (⁠$≃kc/2$⁠) is

after setting f = 4C_f, Re_D = 2 Re_τ U∕u_τ, and $U/uτ=2/Cf$⁠, and replacing the logarithm with a natural one. If we note that the empirical coefficient 0.407 in Eq. (13) is, perhaps fortuitously, very nearly equal to the von Kármán constant κ, and if we replace it with κ, we can equate Eqs. (12) and (13) to yield a simple expression for B (and its downward shift relative to smooth-wall flow $ΔB$⁠) as

Since setting $d=ϵ=0$ yields B = 5.77, rather than its usual smooth-wall value of 5.5, the value of B in Eq. (14) was reduced proportionally to compensate for assumptions made in its derivation. With this revised expression for B, pipe-flow computations with the wall-function model and a zero-plane displacement of $d=0.5kc$⁠, with 12 cells of uniform width across the pipe radius, reproduced the value of C_f from the Colebrook formula for flow in pipes with sandgrain roughness to within about 2%, for $0<ϵ/D<0.04$ and 10⁴ < Re_D < 10⁸, as shown in Fig. 11.

In the case of pipe flow, it is possible to generalize the sandgrain roughness velocity-profile model for arbitrary roughness size k_c and Reynolds number by calibrating its parameters to yield the friction coefficients of the formula of Colebrook.²⁸ Coefficient B, which describes the downward shift of the log-linear region on account of roughness, has been shown by Schlichting²⁹ to follow the equation,

for roughness heights of $ks+>70$ and high Reynolds numbers, and so applies in the “fully rough” regime. A general model for B therefore requires low Reynolds number and low roughness height corrections to Eq. (15). A simple low-Reynolds-number correction that led to good agreement between C_f when determined by the Colebrook formula and by this model computation at large roughness heights was to multiply $(1/κ) ln [ks+/3.5]$ by $(1−6/Reτ)$⁠. A compact correction for smaller roughness heights was to multiply it by the factor exp{(max[70∕k_c⁺, 1] − 1) $Reτ$∕110}, where the apparent singularity at low values of $kc+$ is of little concern because thin roughness sublayers do not extend to the log-linear region. Better agreement with the Colebrook formula was found when a further factor of 1.7 was introduced after $ks+$ was replaced by $kc+$ so that a general formula for B is

The powers in this empirical equation were approximated as square roots for convenience and the coefficients were rounded for compactness, since the Colebrook formula for C_f has an uncertainty of up to 10% and a more precise form of Eq. (16) would not necessarily be more accurate. The value of 4 in the logarithmic term was added to model the observation of hydraulic smoothness for $ks+ ≲3.5$⁠, and the beginning of log-linear region was set as y_l⁺ = max[30, 1.6 k_c⁺]. The location of $y0+$—the point at which the mean velocity first exceeds zero—was matched to simulation data and model computations consistent with Colebrook's formula by setting $y0+=10 ln [1+(kc+−4)/10]$⁠, where the offset of 4 enforces hydraulically smooth flow for small sandgrain roughness heights. The zero-plane displacement d⁺ was set to $0.5kc+$ as an approximate match to DNS results.

Computations of the friction coefficient C_f were made with this sublayer velocity-profile model, generalized to describe sandgrain roughness, with the computational scheme described earlier and a standard k-$ϵ$ outer-flow closure, over a wide range of values of $kc+$ and Reynolds number, and are shown in Fig. 12 together with C_f from Colebrook's formula. This calibration of the roughness-sublayer velocity profile is clearly sufficient for the computational model to match the Colebrook friction coefficient, to at least the accuracy of the Colebrook formula itself. The excellent agreement between the model calculation and its calibration equation illustrates that the composite model of the velocity profile within the roughness sublayer, determined from two low-$Reτ$ direct numerical simulations, and generalized for other values of roughness height and Reynolds number, can indeed be combined with a standard outer-layer turbulence closure in a fine-mesh computation to match target values of the integral flow parameter C_f.

The problem of modeling turbulent flow over rough surfaces is complicated by the infinite variety of surface topographies. Each topography corresponds to a particular two- or three-dimensional surface at which a no-slip, and often a no-penetration boundary condition is applied. The different geometries of surfaces over which boundary conditions are applied lead in turn to a variety of flow phenomena within the roughness sublayer. Examples include recirculating flow patterns, isolated from the outer flow, in d-type roughness (Perry et al.³⁰); recirculation directly behind, and interaction with the outer flow further downstream of roughness elements in k-type roughness; modifications of these flow patterns through the mutual sheltering of roughness elements (Jiminéz³¹); enhanced streamwise vortices near roughness crests (Rotta¹); enhanced turbulence production in the wake of certain roughness elements and enhanced redistribution/isotropization of component turbulent kinetic energy in the sandgrain roughness sublayer (Yuan and Piomelli¹⁸); amplification of turbulent motions by interaction between roughness sublayer and outer eddies (Anderson et al.³²); and reduction of turbulence scale by sublayer roughness (Aghaei Jouybari et al.³³).

The phenomena cited above have each been identified in flow over some regular, prescribed roughness. In practice, the texture of surfaces may also be irregular, in which case it is even more difficult to anticipate which flow phenomena or combinations thereof might be present in the sublayer over some arbitrary rough surface, let alone propose how the turbulent and roughness-induced features of such a flow might be modeled accurately. The solution described in this paper is to use the results of direct numerical simulations to prescribe directly either the coefficients B and d⁺ in a roughness-specific wall function or B, d⁺, $y0+$⁠, and $yl+$ in an algebraic model of the mean velocity profile in the roughness sublayer that arises from flow over the rough surface of interest—an approach made possible by advances in direct numerical simulation with immersed-boundary methods—rather than computing it by trying to develop roughness-dependent models for the body forces and stress gradients in the unclosed momentum equations. Extensive libraries of velocity profiles within and beyond the roughness sublayer can now be compiled for flows over many different kinds of surface roughness to serve this purpose, such as, the database of 42 different rough-wall flows described by Aghaei Jouybari et al.¹⁷ 

The application of reference results from one near-wall flow to another relies on their assumed “universality” when expressed in wall units. It also relies on a consistent way of representing velocities in wall units, which is provided by the recognition that when the log-linear part of the boundary-layer velocity profile is shifted by the zero-plane displacement d⁺, as $u+=(1/κ) ln (y+−d+)+B$⁠, it appears to describe all rough-wall flows. Under these assumptions, wall-function models and sublayer velocity-profile models in wall units devised from simulations of rough-wall channel flow may then be presumed to describe other kinds of boundary-layer flows over the same rough surface. A roughness-specific wall function or sublayer velocity-profile model can then be coupled with an outer-flow closure and its particular far-field boundary conditions. These approaches are well suited to k–$ϵ$⁠, k-ω, and Spalart–Allmaras models for attached boundary-layer flows over locally homogeneous roughness, when the mean velocity field is expected to be log-linear beyond about y⁺ of 40. While they appear to perform well in steady, fully developed duct flows, they have still to be tested in developing and transient flows.

The strategy of using the results of fine-grained simulations to describe high-shear roughness sublayers, where they are essential because the processes there are most complex and most difficult to model, and a low-order turbulence model in the low-shear outer flow, is similar to that of several wall-layer models for smooth-wall flow (Piomelli and Balaras³⁴). In the context of rough-wall flows, it has been applied as DNS-derived, roughness-specific, wall-function and sublayer velocity-profile models together with an outer-flow k-$ϵ$ model. The strategy was tested by first checking that the DNS-informed models adequately reproduced the mean velocity profiles and friction coefficients of the simulations from which they were deduced: channel flows at $Reτ=1000$ over surfaces of four very different kinds of roughness: sandgrains, wavy walls, walls with spaced arrays of semi-ellipsoids, and walls with random roughness. The largest discrepancy between the corresponding model computations of skin friction C_f and simulation results was less than 6%. The sandgrain roughness models were then used to predict successfully Colebrook's equation describing fits to experimental pipe-flow data of C_f at much higher Reynolds numbers.

The database used in this study included many kinds of surfaces, at a small number of different roughness heights, with flow at two different values of $Reτ$—many in the “fully rough” regime. For flow over surfaces of smaller roughness heights or at lower Reynolds numbers in the transitionally rough regime, the roughness-specific sublayer velocity profile must be modified or recalculated from a surrogate simulation at those conditions. Potential modifications are to rescale the model sublayer velocity-profile parameters with multiplicative corrections for Reynolds number and roughness height, or to re-express B as a function of the roughness size in wall units, in the wall-function model. For the case of sandgrain roughness, a four-parameter model was calibrated in this way so that the modeled sublayer velocity profile, together with an outer-flow k–$ϵ$ closure, reproduced Colebrook's description of C_f as a function of $Red$ for pipe flow with good accuracy. An algebraic equation deduced for the value of B in the wall-function model, for transitionally and fully rough flow over sandgrain roughness, was equally accurate.

For surfaces with roughnesses for which no simulation result is available, one can either carry out the required DNS or attempt to infer the likely values of model parameters by interpolation within a database. In the related work of Aghaei Jouybari et al.,¹⁷ it has been shown how the machine-learning technique of Gaussian process regression can be used to predict the equivalent sandgrain roughness height as a function of as many as nine different surface texture parameters, after model training with an extensive DNS database of rough surfaces. The same approach can be used to determine estimates of B, d⁺ and other model parameters for some new rough surface, together with a companion prediction of their likely accuracy.

Given the complexity of modeling near-surface phenomena in rough-wall flow, the capabilities of direct numerical simulation, and the adequacy of low-order turbulence models in the low-shear outer flow, these DNS-informed wall-function and sublayer velocity-profile modeling approaches seem to be most promising for the general problem of accurate prediction of boundary-layer flow over surfaces of arbitrary roughness.

The authors gratefully acknowledge the support of the Office of Naval Research under Grant No. N00014-17-1-2102.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

A simple composite description is presented of the superficially averaged mean velocity profile from the lowest point in the roughness sublayer at $y+=0$ to just beyond the roughness crest at $y+=kc+$⁠. It appears to describe quite accurately the mean velocity in the roughness sublayer for a wide range of different rough surfaces, using four parameters. For surfaces of large roughness, it is not unusual for $kc+$ to exceed 100, which would be within the log-linear part of the corresponding smooth-wall boundary-layer velocity profile, and it is convenient to begin development of the velocity-profile model there.

In smooth-wall flow, functional analysis of the equivalence within an overlap region of velocity profiles which scale on (i) inner variables y, $uτ$⁠, and ν; and (ii) outer variables y, δ, $uτ$⁠, and U_e, requires that the velocity in the overlap region is logarithmic in $y+$ so that

where U_e is the far-field or centerline velocity and the constants κ, B, and A take values of 0.4, 5.5, and between 0.65 and 2.35, respectively, depending on the kind of flow.

In rough-wall flows, mean velocity profiles determined by both experiment and simulation also indicate the presence of such a log-linear region, though shifted downward to lower values of $u¯/uτ$ at a given $y+$⁠, and with slope $1/κ$ only when $y+$ is replaced by $y+−d+$⁠, where d⁺ is the zero-plane displacement defined in Sec. VII A, or some comparable y-offset lengthscale. Since the geometry of the rough wall introduces additional dimensionless variables to the functional analysis problem, which may describe the roughness through k_c, k_rms, skewness, etc., the observation of an overlap region logarithmic in $y+$ implies that

where C, D, and d⁺ are constants for any given roughness but, in general, functions of dimensionless roughness variables. If this observation of a logarithmic overlap region applies to all rough-wall boundary layers, as appears to be the case from the simulations of Aghaei Jouybari et al.¹⁷ over many different rough surfaces, one can write a general form for D as

where the function $ΔU+$ is zero for smooth walls but is otherwise non-zero. Under these assumptions, a smooth-wall velocity-profile model for the log-linear region can describe a rough-wall flow equally well when the value of B in Eq. (A1) is reduced below 5.5 and the $y+$ offset of d⁺ is included, as in Eq. (A2).

It appears to be simplest to construct a general analytical model which spans the near-wall and log-linear regions by blending a cubic polynomial with Eq. (A2). If we introduce the variable $y∗=y+−d+$ and identify the points $y0∗$ and $yl∗$ (or $y0+$ and $yl+$⁠) as the locations at which $u+$ first exceeds zero and the innermost extent of the log-linear region, a simple composite model for $u+$ in rough-wall boundary layers is then

where the coefficients of the cubic polynomial $f(y∗)$ are chosen as conditions on $u+$ and $du+/dy+$ to ensure smooth blending at $y∗=y0∗$ and $y∗=yl∗$⁠. The cubic polynomial may be constructed in either $y+$ or $ln y+$ space, depending on which provides a better match to reference data for any particular rough surface. This composite model reverts to its smooth-wall equivalent when $d=y0=0$⁠, B = 5.5, and $y∗=y+$⁠.