MODELING OF CUBE ARRAY ROUGHNESS; RANS, LES AND DNS

Flow over arrays of cubes is an extensively studied model problem for rough wall turbulent boundary layers. Not only are cube elements notionally representative of the surface roughness encountered in real engineering applications such as urban canopies and heat exchangers, but a great deal of variation in drag and heat transfer can be realized through simple measures such as changing inter-cube spacing and orientation. While considerable research has been performed in computationally investigating these topologies using DNS and LES, the ability of sublayer-resolved RANS to predict the bulk ﬂow phenomena of these systems is relatively unexplored, especially at high packing densities. Here, RANS simulations are conducted on six different packing densities of cubes in aligned and staggered conﬁgurations. The packing densities investigated span from what would classically be deﬁned as isolated, up to those in the d-type roughness regime. Three different sublayer-resolved turbulence closure models were tested for each case; a low Reynolds number k-epsilon model, the Menter k-omega SST model, and a full Reynolds stress model. Comparisons of the velocity ﬁelds, secondary ﬂow features, and drag coefﬁcients are made between the RANS results and existing LES and DNS results from the literature. There is a high degree of variability in the performance of the various RANS models across all metrics of comparison. The Reynolds stress model however demonstrated the best predictive value in terms of the mean velocity proﬁle 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. At large scales, atmospheric flows interact with buildings to create highly complex turbulent structures around urban canopies [1]. 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. [2]. 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 problem was long restricted to experimental observations. Early studies of friction loss in sand-grained roughened pipes by Nikuradse [3] and Moody [4] spawned seventy years of experiments covering a wide rage of surface roughness morphologies [5,6]. These experimental results have allowed for 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 [7,8].
In addition, with the growth of additive manufacturing as a viable alternative to conventional metallurgy, there is active research in the characterization of how the variable and complex roughness of an AM surface affects friction and heat transfer; see Stimpson et al. [9], Kirsch and Thole [10], and Hanson et al. [11].
However, due to the high overhead costs and limited bandwidth of experimental investigation, using CFD as a predictive method for rough wall turbulent boundary layer problems 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 very computationally expensive at the Reynolds numbers encountered in many real engineering applications, and become more prohibitively so for real rough surfaces. 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 [12][13][14][15]. The single-cube data provides 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. [16] , Castro et al. [17] and Perret et al. [18] and span low to moderate surface packing densities. Coceal et al. [19] 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 [20]. LES results can be found in Inagaki et al. [21], Kono et al. [22], Kanda et al. [23] and Yang et al. [24][25][26].
While much of the investigation of these roughness arrays has been through the use of turbulence resolving CFD (DNS and LES), the question of whether accurate prediction of the interactions between turbulent fluctuations and the mean flow is necessary to capture bulk performance remains. In order to explore this, Reynolds Averaged Navier-Stokes (RANS) simulations, which do not resolve any turbulence, can be used.
Xie and Castro [27], Cheng et al. [28] and Santiago et al. [29] 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, and little attention has been paid to very low packing densities, and very high packing densities.
Here, sublayer-resolved RANS simulations are carried out on a range of 6 packing densities with 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 prediction of certain mean or secondary flow features require scale resolving tools. If not, the significantly reduced computational cost of RANS compared to LES and DNS could be of benefit for fluid dynamics practitioners, and allow more computationally efficient investigation of internal rough wall flow without sacrificing fidelity.
This paper is organized as follows. First, the different geometric configurations investigated are delineated, 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 and recommendations for further investigation are made.

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 y , with the front face perpendicular to the flow direction, which is from left to right. For cubes in the staggered arrangement, as seen in Fig. 2, every other column of cubes is shifted in the spanwise direction by half of the separation distance L y . This causes the separation distance between cubes with the same spanwise coordinate to have a separation distance of 2L x . 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 L z , and can be seen from 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 often refers to the .08%, 1.0%, and 4.4% cases collectively as Low Packing Density cases (LPD), and 25%, 50% and 70% together as the High Packing Density cases (HPD).
At very low packing density, λ p < O(3%), wall mounted cubes can be classified as isolated roughness [25]. In this regime, the cubes are separated enough that they have a limited effects 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 wake of the upstream roughness elements directly affects those downstream, a phenomena which has been termed sheltering [26,30,31]. 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. Here, 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 [6]. In effect, the flow tends to skim over the top of the cubes without significant interaction with the bottom wall.
The height of the wall mounted cube is h is set equal to 1.0 for all configurations, and the half channel height, L z , is 3.5h for the three low packing density cases, and 4.0h for the three higher packing densities, consistent with Yang et al. [25], Coceal et al. [19], and Xu et al. [32]. 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, and that the vertical height of the domain in comparison to the roughness height is sufficient to capture the flow.
While there are multiple methods of defining the friction velocity in a rough wall setting, this paper uses the definition given by Eqn. 1.  Here, f b is the volumetric body force. The attendant friction Reynolds number is defined as For the sparse packed cubes, Re τ =4735, consistent with Yang et al. [25], and Re τ =500 for all of the high packing density cases, consistent with Coceal et al. [19] and Xu et al. [32]. 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 cases. Cases are identified by their case ID, with A for Aligned and S for Staggered, followed by integers indicating packing density. The RANS simulations were conducted using STARCCM+ (v2019.2) [33], 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 τ . 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 [34], and the Lien k-ε model [35]. In addition, we apply a sublayer resolved Reynolds Stress Transport (RST) model based on Launder and Shima [36]. RSTs are more computationally expensive than eddy viscosity models as they solve modelled 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 [37] and Rodi [38]. 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 Fig. 1 and Fig. 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. Figure 4 shows a cross-section of the A440 grid for reference. The cell count ranges from 200,000 for the A70 case where the packing density is highest, to 6.0 million cells for the S008 case at the other end of the test matrix. A grid study was performed leading to certainty of grid independence in the presented 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 (ISA) can be used to represent averaged flow quantities [39]. All averages used herin will be 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 5 and 6 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. [25] for these cases is 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 performers the best for both A100 and S100, while k-ε performer the best for both A440 and S440.  In the LES of 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 secondary flow phenomena will be discussed further in the next section, these features have a pronounced impact upon the average profiles, especially above the cube height, with the difference especially striking for A100 case. For both the aligned and staggered .08% cases, which approach the limit of a smooth channel flow, the SST model results provide the best agreement with the LES data, while both k-ε and RST predict that the velocity deficit caused by the roughness to be greater than what is expected, although the differences between the three turbulence models is the smallest among the low packing density cases. The Log-Law line in Fig. 5 and Fig. 6 is that of a smooth channel.
We now move to the HPD cases. Figure 7 shows the CSA for the S25 RANS results, with DNS from Coceal et al. [19] for comparison. In addition, wall modeled RANS data produced by Santiago et al. [29] is included. Here again, we observe that SST produces a significant center-line velocity overestimation, while both k-ε and RST provide a much more accurate prediction, and outperform the existing wall-modeled results.   Figure 8 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 [27] are shown as well. Locally, both the wall resolved k-ε and RST velocity profiles conform well to 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. 9, 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. [32].  Figure 10 details the velocity profiles for the other d-type cases. For case A50, the RST performs significantly better than its other RANS counterparts when compared to the DNS of Xu et al. [32]. 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 is 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 accuracy decreasing as it moves to higher or 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. Here U x , V x and W x are the the three components of velocity comprehensively averaged in the streamwise direction. Figures 11, 12, and 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 transfer low momentum fluid from below the cube height to the upper regions of the boundary layer, causing spanwise variation in the streamwise velocity. This also leads to the comparatively lower center-line mean velocity of the RST case observed earlier when compared to the 2-equation models. The maximum in-plane velocity for the RST case is on the order of ten percent of the center-line velocity, which is stronger than either the k-ε or SST results, and is consistent with the results of Yang et al. [25]. Of added interest here is that while the SST model better reflects the correct secondary flow behavior, this does not reflect in an accurate mean velocity profile.
Further evidence of the difference in the flow fields can be seen in Figure 14, 15, and 16, which depict contours of streamwise velocity on a horizontal plane 3h from the wall for A100. Both the SST and RST models predict the presence of a high momentum pathway above the roughness strip (y=0), and low mo-  mentum pathways between the roughness strips, with the magnitude of the pathways being comparatively larger in the RST case. This is the appropriate trend as evidenced by the LES data (see Yang et al. [25] for corresponding figures), but the opposite is observed for the k-ε model.
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. 17, with RANS results on the left and DNS from Xu et al. [32] 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 show more homogeneity in the spanwise direction above the cube. Focusing in on the region below the height of the cube, Fig. 18 depicts contours of streamwise velocity as well as in-plane streamlines at Z/h=.8 for A50 RST. RANS is able to correctly predict the presence of the arch-vortex in the wake of the cubes, but marginally over-predict the streamwise velocity in the trench to either side of the cubes compared to the DNS.

REYNOLDS AND DISPERSIVE STRESSES
Another area of interest in this work is how accurately the RANS RST model predicts the various components of the Reynolds stress. The XZ components of the comprehensively   [32] averaged Reynolds and dispersive stress tensors for A100 RST are plotted in Fig. 19, along with the corresponding LES data from Yang et al. [25]. The RANS data shows excellent agreement with the LES data above the height of the roughness, but less so within the roughness layer itself. Indeed, the RANS dispersive stress is more consistent with that of a higher packing density than one percent, which have positive values below the cube height. Additionally, for much of the domain, the sum of the Reynolds and dispersive stresses (τ + T ) very closely matches the theoretical line τ + T = −1 + Z L z , which is plotted as grey dashes in Fig. 19.
In the case of a higher packing density arrangement, such as A50, the RST model provides a somewhat different picture. In Fig. 20, we see that RANS over-predicts the normal components of the Reynolds stress compared to the DNS data in Xu et al. [32], although there is excellent agreement with the XZ component. The error in the normal components is especially evident at the height of the cube.

DRAG PARTITION
The drag partition and drag coefficients are also examined for comparison, as these are often quantities of interest in engineering applications. The drag coefficient is evaluated using Equation 3, 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.
(3) Figure 21 shows the drag coefficient for all of the sparse cube cases, once again with LES data from Yang et al. [25] 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 over-prediction of the velocity profiles, leading to a larger than expected value of U h . The second is that the RST model pro-vides the better prediction for the A/S008 and A/S100 cases, but not for A/S440.  Figure 22 details the ratio of drag force on the channel substrate F S to the drag force on the cubes F D as a function of packing density for the aligned HPD cases, with Xu et al. [32] DNS data 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.

CONCLUSION
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, all of the models were incapable of capturing the wide swath of physics across the test matrix. The RANS RST model had the most consistent performance across the roughness configurations, especially in terms of mean velocity and drag partition. The periodic nature of these geometries

ACKNOWLEDGMENT
This material is based upon work supported by the Department of Energy under Award Number(s) DE-FE0001730. Disclaimer: 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.