## Abstract

Flow over arrays of cubes is an extensively studied model problem for rough wall turbulent boundary layers. While considerable research has been performed in computationally investigating these topologies using direct numerical simulation (DNS) and large eddy simulation (LES), the ability of sublayer-resolved Reynolds-averaged Navier–Stokes (RANS) to predict the bulk flow phenomena of these systems is relatively unexplored, especially at low and high packing densities. Here, RANS simulations are conducted on six different packing densities of cubes in aligned and staggered configurations. The packing densities investigated span from what would classically be defined as isolated, up to those in the d-type roughness regime, filling in the gap in the present literature. Three different sublayer-resolved turbulence closure models were tested for each case: a low Reynolds number k–*ϵ* model, the Menter k–*ω* SST model, and a full Reynolds stress model. Comparisons of the velocity fields, secondary flow features, and drag coefficients are made between the RANS results and existing LES and DNS results. There is a significant degree of variability in the performance of the various RANS models across all comparison metrics. However, the Reynolds stress model demonstrated the best accuracy in terms of the mean velocity profile as well as drag partition across the range of packing densities.

## Introduction

Rough wall turbulent boundary layers are a common feature in many real engineering applications. For example, at large scales, atmospheric flows interact with buildings to create highly complex turbulent structures around urban canopies [1,2]. At smaller scales, the performance of internal cooling channels in turbine blades and heat exchangers is directly impacted by the surface roughness, which depends on material properties, wear, and other environmental factors as described in Bons et al. [3]. In order to accurately predict pressure drop and heat transfer in open channel and internal flows, understanding the effects that surface roughness has on the underlying physics is essential.

Investigation of this class of turbulent boundary layer problems was long restricted to experimental observations. Early studies of friction loss in sand-grained roughened pipes by Nikuradse [4] and Moody [5] spawned seventy years of experiments covering a wide range of surface roughness morphologies [6–9]. These experiments have led to the development of a host of correlations that attempt to predict the effects of deterministic and irregular roughness on the mean flow as a function of the surface statistics [10–12].

In addition, with the growth of additive manufacturing as a viable alternative to conventional machining, there is active research in the characterization of how the variable and complex roughness of an additive manufacturing surface affects friction and heat transfer; see Stimpson et al. [13], Kirsch and Thole [14], and McClain et al. [15].

However, due to the high overhead costs and limited bandwidth of experimental investigation, using computational fluid dynamics (CFD) as a predictive method for rough wall turbulent boundary layers has been turned to as an alternative, although there are still several shortcomings with this approach. Scale resolving methods such as direct numerical simulation (DNS) or large eddy simulation (LES) are computationally expensive at the Reynolds numbers encountered in many engineering applications, and become prohibitively so for real rough surfaces [16]. That being the case, simpler geometries are often relied upon.

Cube arrays are a common surrogate for more complex roughness morphologies. Not only are the elements notionally representative of the surface roughness encountered in many real applications, but a great deal of variation in drag and heat transfer can be realized through simple measures such as changing intercube spacing or orientation. Cubes also have the added benefit of having been extensively studied in their isolated form, in both experimental and numerical contexts [17–20]. The single-cube data are a benchmark for the lower limit of packing density. In addition, due to their comparative geometric simplicity and periodicity, these cube arrays have been extensively investigated using DNS and LES.

Experimental results of flow over cubic roughness elements can be found in Cheng et al. [21], and Castro et al. [22], and span low to moderate surface packing densities. Coceal et al. [23] performed DNS on both staggered and aligned cube arrays at 25% packing density, and those results served as a benchmark for many of the following computational investigations. DNS data for a wider variety of packing densities was further reported in Leonardi and Castro [24], and Xu et al. [25]. LES results can be found in Inagaki et al. [26], Kono et al. [27], Kanda et al. [28], Yang et al. [29–31], and Li et al. [32,33].

While much of the investigation of these roughness arrays has been carried out using turbulence-resolving CFD (DNS and LES), the significance of the interactions between turbulent fluctuations and the mean flow on bulk performance predictions remains unresolved. In order to explore this, Reynolds-averaged Navier–Stokes (RANS) simulations, which do not resolve any turbulence, can be used. Xie and Castro [34], Cheng et al. [35], and Santiago et al. [36] have explored using RANS on cube geometries with varying success depending upon the packing density and RANS turbulence closure model used. However, existing RANS studies have relied almost exclusively upon wall-modeled forms of the turbulence closure equations. Additionally, little attention has been paid to very low packing densities, and very high packing densities.

Here, sublayer-resolved RANS simulations are carried out on wall mounted cubes for a range of six packing densities in two different configurations; aligned and staggered, using three well-established turbulence models. By comparing steady RANS results to high pedigree LES and DNS, the accuracy and drawbacks of this tool can be assessed, and it can be determined whether the prediction of certain mean or secondary flow features requires scale-resolving tools. If not, the significantly reduced computational cost of RANS could be of benefit for fluid dynamics practitioners, allowing more computationally efficient investigation of internal rough wall flows.

The remainder of the paper is organized as follows: First, the different geometric configurations investigated are described, and the computational procedures are reported. Then, results from a series of numerical studies on multiple arrays of cubic roughness elements at various packing densities are reported, and comparisons are made to existing data. Finally, an assessment of the capabilities and drawbacks of RANS is given.

## Approach

All of the geometric configurations considered in this study are of wall-mounted cubes in a half-channel. Two different configurations are considered: aligned and staggered. Figures 1 and 2 show top-down views of these arrangements, respectively. *X* is the streamwise coordinate, *Y* is the spanwise coordinate, and *Z* is the wall normal coordinate. In the aligned configuration of Fig. 1, cubes are each separated in the streamwise and spanwise directions by equivalent distances *L _{x}* and

*L*, with the front face perpendicular to the flow direction, which is from left to right. For cubes in the staggered arrangement in Fig. 2, every other column of cubes is shifted in the spanwise direction by half of the separation distance

_{y}*L*. This causes the separation distance between cubes with the same spanwise coordinate to have a separation distance of 2

_{y}*L*. While there are multiple combinations of streamwise and spanwise separation distances that can yield the same packing density, this pairing is the one most commonly used in the literature. The half-channel height is

_{x}*L*, and the cube side length is

_{z}*h*, which can be seen in the isometric view of the topology in Fig. 3.

By varying the magnitude of the cube separation distances, different packing or surface coverage densities, *λ _{p}*, are achieved. Surface coverage density is defined as the ratio of area obstructed by cubes to the total ground area of the channel. For each configuration, six different packing densities are investigated, 08%, 1.0%, 4.4%, 25%, 50%, and 70%. For the sake of brevity, this paper refers to the.08%, 1.0%, and 4.4% cases collectively as low packing density (LPD), and the 25%, 50%, and 70% cases collectively as high packing density (HPD).

At very low packing density, $\lambda p<O(3%)$, wall-mounted cubes can be classified as isolated roughness [30]. In this regime, the cubes are separated enough that they have limited effect on one another, and thus act as isolated or single roughness elements. As packing density increases to more moderate values, cube arrays begin to exhibit *k*-type roughness behavior. In this regime, the wakes of the upstream roughness elements directly affect those downstream, a phenomena, which has been termed sheltering [31,37,38]. In this study, the staggered 25% case is in the *k*-type regime. The aligned 25% case, and the rest of the HPD configurations studied, are d-type in nature. There, the packing density is high enough that there is little active momentum exchange between the outer flow above the roughness height, and the flow in the region below and between the elements [7]. In this regime, the total drag on the substrate is small compared to the total drag on the roughness element.

The height of the wall-mounted cube, *h*, is set equal to 1.0 for all configurations, and the half channel height, *L _{z}*, is 3.5

*h*for the three low packing densities, and 4.0

*h*for the three higher packing densities, consistent with Yang et al. [30], Coceal et al. [23], and Xu et al. [25]. These studies have shown that the roughness sublayer, or the region of the flow affected by the presence of the wall-mounted roughness, is relatively thin in comparison to the roughness height.

*f*is the volumetric body force that drives the flow. The attendant friction Reynolds number is defined as

_{b}For the sparse packed cubes, $Re\tau \u22484700$, consistent with Yang et al. [30]. For the high packing density cases, $Re\tau =500$, consistent with Xu et al. [25]. At these Reynolds numbers, the flow can be considered to be fully rough at most intermediate packing densities, and becomes transitionally rough at either end of the packing density distribution. Table 1 lists spacing information for all runs. Case identifiers start with A for aligned and S for staggered, followed by integers indicating packing density.

Case ID | λ_{p} | $Lx\xd7Ly\xd7Lz$ |
---|---|---|

A70 and S70 | 70% | $1.195\xd71.195\xd74.0$ |

A50 and S50 | 50% | $1.41\xd71.41\xd74.0$ |

A25 and S25 | 25% | $2.0\xd72.0\xd74.0$ |

A440 and S440 | 4.4% | $4.8\xd74.8\xd73.5$ |

A100 and S100 | 1.0% | $10\xd710\xd73.5$ |

A008 and S008 | 0.081% | $35\xd735\xd73.5$ |

Case ID | λ_{p} | $Lx\xd7Ly\xd7Lz$ |
---|---|---|

A70 and S70 | 70% | $1.195\xd71.195\xd74.0$ |

A50 and S50 | 50% | $1.41\xd71.41\xd74.0$ |

A25 and S25 | 25% | $2.0\xd72.0\xd74.0$ |

A440 and S440 | 4.4% | $4.8\xd74.8\xd73.5$ |

A100 and S100 | 1.0% | $10\xd710\xd73.5$ |

A008 and S008 | 0.081% | $35\xd735\xd73.5$ |

The RANS simulations were conducted using starccm+ (v2019.2) [39], a finite volume CFD package that is widely used in industry and research settings. All cases employed second-order spatial discretization to solve the incompressible RANS equations. Flow is driven using a constant volumetric body force *f _{b}*, chosen to achieve the desired $Re\tau $. In this study, we focus on the effects of employing sublayer resolved forms of the RANS turbulence closure equations. All cases are simulated using three different models. We apply two different widely used two-equation eddy viscosity models; the Menter-

*k*–

*ω*SST model [40], and the Lien

*k*–

*ϵ*model [41]. In addition, we apply a sublayer resolved Reynolds stress transport (RST) model based on Launder and Shima [42]. RSTs are more computationally expensive than eddy viscosity models as they solve modeled transport equations for the six individual Reynolds stresses themselves (plus a length scale providing dissipation scalar). However, they are also more accurate for many flows, in particular where stress anisotropy significantly impacts the mean-flow evolution. For the pressure strain term, we rely on a two-layer formulation based on Gibson and Launder [43], and Rodi [44]. All three turbulence models are available in starccm+, and no changes were made to any of the standard model constants.

One period of the cube array is simulated for each configuration, where each of the four sides of the computational domain is a cyclic boundary condition, and the top of the domain at the half-channel height is a symmetry boundary. The dashed lines in Figs. 1 and 2 enclose the periodic domain. The cube surface and substrate are treated as no-slip walls. Because the domain is periodic, no specific inflow or far-field boundary conditions are needed for the RANS and turbulence transport equations, which reduces the variability in turbulence model performance uncertainty. A structured wall-resolved mesh is used for all cases, with a *y*^{+} less than unity in conformance with the low Reynolds number forms of the turbulence models used. Across all 36 cases, 0.93 was the maximum value of *y*^{+} in the cell adjacent to the wall. Figure 4 shows a streamwise cross section of the A440 grid for reference. The cell count ranges from 100,000 for the A70 case where the packing density is highest, to 6.0 million cells for the S008 case where the packing density is lowest. In order to assess the resolution of the grid, a mesh refinement study was performed on the A25 RST case. The baseline grid of 200,000 cells was refined by a factor of 2 in every direction, leading to a fine grid of 1.6 million cells. A coarse grid of 25,000 cells was also developed. The mean velocity profile and drag partition results were compared for the fine, coarse, and baseline meshes. The velocity profile results from all three of these meshes are included in Fig. 5. The maximum difference in the mean velocity profiles is less than 2%. This demonstrates that the baseline mesh is sufficiently fine to consider the results grid-independent.

## Results

### Velocity Profiles.

First, we examine the velocity profiles for the low packing density arrays. For rough wall turbulent boundary layers, either comprehensive spatial averaging (CSA) or intrinsic spatial averaging can be used to represent averaged flow quantities [45]. All averages in this paper use CSA, i.e., with the roughness occupied space included in the spatial average with zero velocity. *U _{xy}* denotes the streamwise velocity, comprehensively averaged in both the streamwise and spanwise directions. Figures 6 and 7 show the CSA of streamwise velocity, normalized by the friction velocity for the three sparse packing densities, in aligned and staggered configurations, respectively. The LES data from Yang et al. [30] for these cases are also plotted for comparison.

For four of the cases A440, S440, A100, and S100, all three RANS models overestimate the velocity at the channel centerline, although better agreement is observed below the height of the roughness elements (*Z* = *h* = 1). The SST model for these four cases diverges the most from the LES data, with both *k*–*ϵ* and RST more closely reflecting the LES profile. RST performs the best for both A100 and S100, while *k*–*ϵ* performs the best for both A440 and S440.

In the LES results of the aforementioned four cases, a flattening of the velocity profile above the roughness element is observed. This is caused by strong secondary flow, which transports low velocity fluid from below the roughness height to the center of the channel, causing a nearly flat averaged velocity profile above the cube height. While this flow phenomena will be discussed further in the section Secondary Flow Features, these features have a pronounced impact upon the average profiles, especially above the cube height. This is especially striking for A100.

For both the aligned and staggered 0.08% cases, which approach the limit of a smooth channel flow, the SST model results provide the best agreement with the LES data. Both *k*–*ϵ* and RST predict that the velocity deficit caused by the roughness is greater than expected. The Log-Law line in Figs. 6 and 7 is that of a smooth channel.

We now move to the HPD cases. Figure 8 shows the CSA for the S25 RANS results, with DNS from Coceal et al. [23], for comparison. In addition, wall-modeled RANS data from Santiago et al. [36] is included. Here again, we observe that SST produces a significant centerline velocity overestimation, while both *k*–*ϵ* and RST provide a much more accurate prediction, and outperform the existing wall-modeled results.

Figure 9 depicts the local velocity profiles at a point directly between the cubes (location depicted by the red X included in the sketch). For comparison, both LES and wall-modeled RANS data from Xie and Castro [34] are shown as well. Locally, both the wall-resolved *k*–*ϵ* and RST velocity profiles agree well with the LES data, and are almost identical to their wall-modeled counterparts. Once again, SST performs poorly both above and below the cube height. Turning to the aligned case and Fig. 5, there is a greater spread in the performance of the three models for A25 compared to S25. While SST remains inaccurate and *k*–*ϵ* produces a small over-prediction, RST demonstrates excellent agreement with the DNS results reported in Xu et al. [25].

Figure 10 details the velocity profiles for the other d-type cases. For case A50, the RST performs significantly better than the other RANS models when compared to the DNS of Xu et al. [25]. For both A70 and S70, all three RANS models give very similar predictions for the profiles above the cubes, although all overestimate the velocity compared to the DNS in the case of A70.

Based on these results, it is clear that the predictive capabilities of RANS are highly sensitive to the turbulence model used. The RST performs the best on average across the case matrix, with SST performing the worst. This tends to support the hypothesis that capturing the effects of turbulence anisotropy is necessary to predict the flow fields produced by these cubes. For the low packing density cases, RANS provides stronger predictive capability below the height of the roughness, with performance deteriorating with increasing height above the roughness. Performance is similar to that of some of the wall-modeled results previously published. In the low packing density region, the LES results show that there is very little difference between the aligned and staggered profiles for a given packing density. This trend is also generally reflected in the RANS results.

In addition, there is no clear correlation between the packing density and the performance of a given turbulence model. The RST, for example, produces the best results for A25, with less accuracy at higher and lower packing densities.

## Secondary Flow Features

Also of interest in this work is the question of the predictive capability of RANS in the context of secondary flows. Figures 11–13 show streamwise averaged streamwise velocity contours, with in-plane streamwise averaged velocity vectors for the A100 geometry (velocity vectors are shown on only half the domain for clarity of visualization). The scaling of the in-plane velocity is shown on the left half of the figure. Here again, we observe more evidence of the wide variation in performance between the turbulence models. While the *k*–*ϵ* model predicts very little in the way of secondary flow, the RST model predicts the presence of strong counter-rotating vortices on either side of the cube. These vortices transport low momentum fluid from below the cube height to the upper regions of the boundary layer, leading to spanwise variation in the streamwise velocity. This also leads to the comparatively lower centerline mean velocity of the RST case observed earlier when compared to the two-equation models. The maximum in-plane velocity for the RST case is of the order of ten percent of the centerline velocity, which is higher than either the *k*–*ϵ* or SST results, and is consistent with the LES results of Yang et al. [30]. Of added interest here is that while the SST model better captures the secondary flow behavior, this does not correspond to an accurate mean velocity profile.

Further evidence of the difference in the flow fields can be seen in Figs. 14–16, which depict contours of streamwise velocity on a horizontal plane 3 h from the wall for A100. Both the SST and RST models predict the presence of a high momentum pathway above the roughness strip (i.e., near *Y* = 0), and low momentum pathways between the roughness strips. The spatial variation of velocity that characterizes these pathways is comparatively larger in the RST case. This is the appropriate trend as evidenced by the LES data (see Ref. [30] for corresponding figures), but the opposite trend is observed for the *k*–*ϵ* model.

The location of the high momentum pathway is directly influenced by the strength of the previously described counter-rotating vortices. When this secondary flow feature is strong, as in the A100 RST case, the low momentum fluid that is transported by these vortices reduces the local velocity between the cubes, and the high momentum pathway appears along the cube centerline. When the vortices are weak (k-$\u03f5)$, the flow is comparatively faster in the unobstructed channels on either side of the cube, and a low momentum strip is observed above the cube centerline.

Similar trends are observed for staggered cube arrangements as well. The S100 flow fields for all three turbulence models are plotted in Figs. 17–19. There is a large variation between the models in the strength of the secondary flow features, and again the RST model returns the most accurate result when compared with the LES data.

To further examine the relationship between a strong mean flow prediction and the accuracy of secondary flows, we turn to the HPD case A50 RST. Streamwise velocity contours at *X* = 0 (the cube centerline) are shown in Fig. 20, with RANS results on the left and DNS from Xu et al. [25], on the right. The black line represents the contour line passing through *Z/h* = 1.2 at *Y* = 0. There is a weak concavity to this line, suggesting the presence of a weak high momentum pathway above the cube. RANS here slightly overestimates the size of this feature, as the DNS shows more homogeneity in the spanwise direction above the cube. In Fig. 21, we turn to a top-down view of the A50 cubes at a height of *Z/h* = 0.8 (slightly below the roughness height). In the axial gap between the cubes, in-plane streamlines depict discrete counter-rotating vortices, a feature which is present in both the RANS and DNS results.

## Reynolds and Dispersive Stresses

where $\u27e8UW\u27e9xy$ is the product of the streamwise and wall normal velocities, comprehensively averaged in the streamwise and spanwise directions. These dispersive stresses are of interest in applications where spatial averaging is used, such as porous media [46], and discrete element roughness modeling [47,48]. The XZ components of the comprehensively averaged Reynolds and dispersive stress tensors for A100 RST are plotted in Fig. 22, along with the corresponding LES data from Yang et al. [30]. All stresses are normalized by $\rho u\tau 2$. The RANS data show excellent agreement with the LES data above the height of the roughness, but less so within the roughness layer itself. Unlike $\u2212Rxz,\u2212Dxz$ is not a positive definite quantity. Positive values of dispersive stress below the roughness height, as predicted by the RANS for this case, are generally a feature of packing density greater than 1%. Additionally, the sum of the Reynolds and dispersive stresses ($\tau T+$) very closely matches the theoretical value $\tau T+=\u22121+ZLz$, which is plotted as a gray dashed line in Fig. 22. This shows that the simulation is well converged, and that in the outer region, Reynolds and dispersive stresses balance the pressure gradient. While the dispersive stress is small compared to the Reynolds stress in the outer region, it is not zero, and contributes to the wall-normal momentum balance. The stress budget for the S100 RST case, as shown in Fig. 23, tells a similar story to its aligned counterpart. In the case of a higher packing density arrangement, such as A50, the RST model provides a somewhat different picture. In Fig. 24, we see that RANS underpredicts the normal components of the Reynolds stress compared to the DNS data in Xu et al. [25], although there is excellent agreement with the XZ component. The error in the normal components is especially evident at the height of the cube. The XZ component of the dispersive stress is confined almost exclusively to the roughness occupied region in this case, although RANS underpredicts its magnitude as well.

## Drag Coefficients and Drag Partition

where *F* is the drag force on one wall-mounted cube, which includes both the pressure and viscous drag, and *U _{h}* is the comprehensive spatial average of

*U*evaluated at the cube height.

Figure 25 shows the drag coefficient for all of the sparse cube cases, once again with LES data from Yang et al. [30] used for comparison. Here, there are two clear patterns. First, across all cases and turbulence models, RANS generally gives an underestimation of the drag coefficient. This is an artifact of the overprediction of the velocity profiles, leading to a larger than expected value of *U _{h}*. Second, the RST model provides the better prediction for the A/S008 and A/S100 cases, but not for A/S440.

Figure 26 shows the ratio of drag force on the channel substrate *F _{S}* to the drag force on the cubes

*F*as a function of packing density for the aligned HPD cases. DNS data from Xu et al. [25] are used for comparison. With increasing packing density, this ratio decreases as the substrate area shrinks and is increasingly sheltered from the outer flow. Here, RST exhibits good agreement with the DNS.

_{D}## Conclusions

RANS studies have been performed for a suite of twelve cube roughness array configurations for which DNS and LES data are available. The purpose of these studies was to assess the performance of various wall-resolved closure models, and determine whether their reduced computational cost could be leveraged without sacrificing fidelity. Overall, there was a high degree of performance variability between models and across packing densities for the various performance metrics. While the trends in most cases were reasonable, none of the models completely captured the wide variety of physics across the test matrix. The *k*–*ϵ* model displayed some utility in predicting the mean velocity and drag coefficient, but displayed limited capacity for predicting the expected secondary flow features. The standard SST model failed to capture even the mean velocity profile, so we would not recommend this model for this class of roughness. The RST model exhibited the most consistent performance across the roughness configurations, especially in terms of mean velocity, drag partition, and secondary flow features. These results indicate that, for certain quantities of interest, RANS RST models are sufficiently accurate that they can be used in place of more computationally expensive turbulence resolving CFD. We note that there are numerous variants of the *k*–*ϵ*, SST and RSM models available, some of which might improve the presented predictions. For example, it is likely worth exploring the Kato–Launder correction [49] to the production term in the two-equation models, since it is known to improve the prediction of some bluff body flows.

## Acknowledgment

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information disclosed, or represents that its use would not infringe privately owned rights. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

## Funding Data

Department of Energy (Award Number(s) DE-FE0001730; Funder ID: 10.13039/100000015).

## Nomenclature

*C*=_{D}drag coefficient

*D*=_{ij}dispersive stress tensor

*f*=_{b}volumetric body force

*L*=_{x}streamwise domain length

*L*=_{y}spanwise domain length

*L*=_{z}wall-normal domain height

*R*=_{ij}Reynolds stress tensor

- $u\tau $ =
friction velocity

*U*=_{h}mean velocity at cube height

*U*=_{x}comprehensive spatial average in the

*x*-direction*U*=_{xy}comprehensive spatial average in the

*x*–*y*plane*λ*=_{p}packing density

*ν*=kinematic viscosity

*ρ*=density

- $\tau T+$ =
sum of Reynolds and dispersive stress

## References

*k*-ϵ Model With a One-Equation Model Near the Wall