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 [69]. 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 [1012].

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 [1720]. 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. [2931], 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 Lx and Ly, 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 Ly. This causes the separation distance between cubes with the same spanwise coordinate to have a separation distance of 2 Lx. 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 Lz, and the cube side length is h, which can be seen in the isometric view of the topology in Fig. 3.

Fig. 1
Top view of aligned geometry with flow from left to right. Dashed line encloses one period of the array.
Fig. 1
Top view of aligned geometry with flow from left to right. Dashed line encloses one period of the array.
Close modal
Fig. 2
Top view of staggered geometry with flow from left to right. Dashed line encloses one period of the array.
Fig. 2
Top view of staggered geometry with flow from left to right. Dashed line encloses one period of the array.
Close modal
Fig. 3
Isometric sketch of topology
Fig. 3
Isometric sketch of topology
Close modal

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, λ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, Lz, 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.

While there are multiple methods of defining the friction velocity in a rough wall setting, this paper uses the definition
uτ=fbLzρ
(1)
Here, fb is the volumetric body force that drives the flow. The attendant friction Reynolds number is defined as
Reτ=uτLzν
(2)

For the sparse packed cubes, Reτ4700, consistent with Yang et al. [30]. For the high packing density cases, Reτ=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.

Table 1

Cube configurations

Case IDλpLx×Ly×Lz
A70 and S7070%1.195×1.195×4.0
A50 and S5050%1.41×1.41×4.0
A25 and S2525%2.0×2.0×4.0
A440 and S4404.4%4.8×4.8×3.5
A100 and S1001.0%10×10×3.5
A008 and S0080.081%35×35×3.5
Case IDλpLx×Ly×Lz
A70 and S7070%1.195×1.195×4.0
A50 and S5050%1.41×1.41×4.0
A25 and S2525%2.0×2.0×4.0
A440 and S4404.4%4.8×4.8×3.5
A100 and S1001.0%10×10×3.5
A008 and S0080.081%35×35×3.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 fb, 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 [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.

Fig. 4
Streamwise cross section of A440 grid
Fig. 4
Streamwise cross section of A440 grid
Close modal
Fig. 5
Mean velocity profiles for A25
Fig. 5
Mean velocity profiles for A25
Close modal

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. Uxy 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.

Fig. 6
Mean velocity profiles for LPD aligned cases
Fig. 6
Mean velocity profiles for LPD aligned cases
Close modal
Fig. 7
Mean velocity profiles for LPD staggered cases
Fig. 7
Mean velocity profiles for LPD staggered cases
Close modal

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.

Fig. 8
Mean velocity profiles for S25
Fig. 8
Mean velocity profiles for S25
Close modal

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].

Fig. 9
Local velocity profiles for S25
Fig. 9
Local velocity profiles for S25
Close modal

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.

Fig. 10
Mean velocity profiles for HPD cases
Fig. 10
Mean velocity profiles for HPD cases
Close modal

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 1113 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.

Fig. 11
Averaged streamwise velocity with in plane velocity vectors for A100 k–ϵ
Fig. 11
Averaged streamwise velocity with in plane velocity vectors for A100 k–ϵ
Close modal
Fig. 12
Averaged streamwise velocity with in plane velocity vectors for A100 SST
Fig. 12
Averaged streamwise velocity with in plane velocity vectors for A100 SST
Close modal
Fig. 13
Averaged streamwise velocity with in plane velocity vectors for A100 RST
Fig. 13
Averaged streamwise velocity with in plane velocity vectors for A100 RST
Close modal

Further evidence of the difference in the flow fields can be seen in Figs. 1416, 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.

Fig. 14
Velocity contours at Z/h = 3 A100 k–ϵ
Fig. 14
Velocity contours at Z/h = 3 A100 k–ϵ
Close modal
Fig. 15
Velocity contours at Z/h = 3 A100 SST
Fig. 15
Velocity contours at Z/h = 3 A100 SST
Close modal
Fig. 16
Velocity contours at Z/h = 3 A100 RST
Fig. 16
Velocity contours at Z/h = 3 A100 RST
Close modal

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-ϵ), 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. 1719. 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.

Fig. 17
Averaged streamwise velocity with in plane velocity vectors for S100 k–ϵ
Fig. 17
Averaged streamwise velocity with in plane velocity vectors for S100 k–ϵ
Close modal
Fig. 18
Averaged streamwise velocity with in plane velocity vectors for S100 SST
Fig. 18
Averaged streamwise velocity with in plane velocity vectors for S100 SST
Close modal
Fig. 19
Averaged streamwise velocity with in plane velocity vectors for S100 RST
Fig. 19
Averaged streamwise velocity with in plane velocity vectors for S100 RST
Close modal

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.

Fig. 20
Streamwise velocity contours at X = 0. Left: A50 RST. Right: DNS Xu et al. [25].
Fig. 20
Streamwise velocity contours at X = 0. Left: A50 RST. Right: DNS Xu et al. [25].
Close modal
Fig. 21
Top-down view of streamwise velocity contours and in plane streamlines at Z/h = 0.8. Left: A50 RST. Right: DNS Xu et al. [25].
Fig. 21
Top-down view of streamwise velocity contours and in plane streamlines at Z/h = 0.8. Left: A50 RST. Right: DNS Xu et al. [25].
Close modal

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 and dispersive stresses. While Reynolds stresses are a result of time averaging, and represent the product of temporally fluctuating velocities, dispersive stresses analogously arise from spatial averaging, and represent the product of spatial fluctuations in velocity. The dispersive stresses are calculated using
Dxz=ρ(UWxyUxyWxy)
(3)

where UWxy 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 ρuτ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 Rxz,Dxz 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 (τT+) very closely matches the theoretical value τT+=1+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.

Fig. 22
Reynolds and dispersive stress for A100 RST
Fig. 22
Reynolds and dispersive stress for A100 RST
Close modal
Fig. 23
Reynolds and dispersive stress for S100 RST
Fig. 23
Reynolds and dispersive stress for S100 RST
Close modal
Fig. 24
Reynolds stress components for A50 RST
Fig. 24
Reynolds stress components for A50 RST
Close modal

Drag Coefficients and Drag Partition

The drag coefficients and drag partition are also examined for comparison, as these are often quantities of interest in engineering applications. The drag coefficient is evaluated using
Cd=FρUh2h2
(4)

where F is the drag force on one wall-mounted cube, which includes both the pressure and viscous drag, and Uh 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 Uh. Second, the RST model provides the better prediction for the A/S008 and A/S100 cases, but not for A/S440.

Fig. 25
Top: drag coefficient for aligned cubes and bottom: drag coefficient for staggered cubes
Fig. 25
Top: drag coefficient for aligned cubes and bottom: drag coefficient for staggered cubes
Close modal

Figure 26 shows the ratio of drag force on the channel substrate FS to the drag force on the cubes FD 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.

Fig. 26
Drag partition for aligned HPD cubes
Fig. 26
Drag partition for aligned HPD cubes
Close modal

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

     
  • CD =

    drag coefficient

  •  
  • Dij =

    dispersive stress tensor

  •  
  • fb =

    volumetric body force

  •  
  • Lx =

    streamwise domain length

  •  
  • Ly =

    spanwise domain length

  •  
  • Lz =

    wall-normal domain height

  •  
  • Rij =

    Reynolds stress tensor

  •  
  • uτ =

    friction velocity

  •  
  • Uh =

    mean velocity at cube height

  •  
  • Ux =

    comprehensive spatial average in the x-direction

  •  
  • Uxy =

    comprehensive spatial average in the xy plane

  •  
  • λp =

    packing density

  •  
  • ν =

    kinematic viscosity

  •  
  • ρ =

    density

  •  
  • τT+ =

    sum of Reynolds and dispersive stress

References

1.
Anderson
,
W.
,
Li
,
Q.
, and
Bou-Zeid
,
E.
,
2015
, “
Numerical Simulation of Flow Over Urban-Like Topographies and Evaluation of Turbulence Temporal Attributes
,”
J. Turbul.
,
16
(
9
), pp.
809
831
.10.1080/14685248.2015.1031241
2.
Giometto
,
M.
,
Christen
,
A.
,
Meneveau
,
C.
,
Fang
,
J.
,
Krafczyk
,
M.
, and
Parlange
,
M.
,
2016
, “
Spatial Characteristics of Roughness Sublayer Mean Flow and Turbulence Over a Realistic Urban Surface
,”
Boundary-Layer Meteorol.
,
160
(
3
), pp.
425
452
.10.1007/s10546-016-0157-6
3.
Bons
,
J.
,
Taylor
,
R.
,
McClain
,
S.
, and
Rivir
,
R.
,
2001
, “
The Many Faces of Turbine Surface Roughness
,”
ASME J. Turbomach.
,
123
(
4
), pp.
739
748
.10.1115/1.1400115
4.
Nikuradse
,
J.
,
1937
, “
Law of Flow in Rough Pipes
,”
NACA
,
Washington, DC
, Report No. 1292.
5.
Moody
,
L.
,
1944
, “
Friction Factors for Pipe Flow
,”
Trans. ASME
,
66
, pp.
671
681
.http://www.ipt.ntnu.no/~asheim/TPG4135/Moody.pdf
6.
Schultz
,
M.
, and
Flack
,
K.
,
2003
, “
Turbulent Boundary Layers Over Surfaces Smoothed by Sanding
,”
ASME J. Fluids Eng.
,
125
(
5
), pp.
863
870
.10.1115/1.1598992
7.
Jimenez
,
J.
,
2004
, “
Turbulent Flows Over Rough Walls
,”
Annu. Rev. Fluid Mech.
,
36
, pp.
173
196
.10.1146/annurev.fluid.36.050802.122103
8.
Coleman
,
H.
,
Hodge
,
B.
, and
Taylor
,
R.
,
1984
, “
A Re-Evaluation of Schlichting's Surface Roughness Experiment
,”
ASME J. Fluids Eng.
,
106
(
1
), pp.
60
65
.10.1115/1.3242406
9.
van Rij
,
J.
,
Belnap
,
B.
, and
Ligrani
,
P.
,
2002
, “
Analysis and Experiments on Three-Dimensional, Irregular Surface Roughness
,”
ASME J. Fluids Eng.
,
124
(
3
), pp.
671
677
.10.1115/1.1486222
10.
Bons
,
J.
,
2010
, “
A Review of Surface Roughness Effects in Gas Turbines
,”
ASME J. Turbomach.
,
132
(
2
), p.
021004
.10.1115/1.3066315
11.
Flack
,
K.
, and
Schultz
,
M.
,
2010
, “
Review of Hydraulic Roughness Scales in the Fully Rough Regime
,”
ASME J. Fluids Eng.
,
132
(
4
), p.
041203
.10.1115/1.4001492
12.
Forooghi
,
P.
,
Stroh
,
A.
,
Magagnato
,
F.
,
Jakirlic
,
S.
, and
Frohnapfel
,
B.
,
2017
, “
Toward a Universal Roughness Correlation
,”
ASME J. Fluids Eng.
,
139
(
12
), p.
121201
.10.1115/1.4037280
13.
Stimpson
,
C.
,
Snyder
,
J.
,
Thole
,
K.
, and
Mongillo
,
D.
,
2016
, “
Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels
,”
ASME J. Turbomach.
,
138
(
5
), p.
051008
.10.1115/1.4032167
14.
Kirsch
,
K.
, and
Thole
,
K.
,
2018
, “
Numerical Optimization, Characterization, and Experimental Investigation of Additively Manufactured Communicating Microchannels
,”
ASME J. Turbomach.
,
140
(
11
), p.
111003
.10.1115/1.4041494
15.
McClain
,
S.
,
Hanson
,
D.
,
Cinnamon
,
E.
,
Snyder
,
J.
,
Kunz
,
R.
, and
Thole
,
K.
,
2021
, “
Flow in a Simulated Turbine Blade Cooling Channel With Spatially Varying Roughness Caused by Additive Manufacturing Orientation
,”
ASME J. Turbomach.
,
143
(
7
), p.
071013
.10.1115/1.4050389
16.
Yang
,
X.
, and
Griffin
,
K.
,
2021
, “
Grid-Point and Time-Step Requirements for Direct Numerical Simulation and Large-Eddy Simulation
,”
Phys. Fluids
,
33
(
1
), p.
015108
.10.1063/5.0036515
17.
Akins
,
R.
,
Peterka
,
J.
, and
Cermak
,
K.
,
1977
, “
Mean Force and Moment Coefficients for Buildings in Turbulent Boundary Layers
,”
J. Wind Eng. Ind. Aerodyn.
,
2
(
3
), pp.
195
209
.10.1016/0167-6105(77)90022-8
18.
Lim
,
H.
,
Thomas
,
T.
, and
Castro
,
I.
,
2009
, “
Flow Around a Cube in a Turbulent Boundary Layer: LES and Experiment
,”
J. Wind Eng. Ind. Aerodyn.
,
97
(
2
), pp.
96
109
.10.1016/j.jweia.2009.01.001
19.
Rodi
,
W.
,
1997
, “
Comparison of LES and RANS Calculations of the Flow Around Bluff Bodies
,”
J. Wind Eng. Ind. Aerodyn.
,
69–71
, pp.
55
75
.10.1016/S0167-6105(97)00147-5
20.
Tominaga
,
Y.
, and
Stathopoulos
,
T.
,
2010
, “
Numerical Simulation of Dispersion Around an Isolated Cubic Building: Model Evaluation of RANS and LES
,”
Build. Environ.
,
45
(
10
), pp.
2231
2239
.10.1016/j.buildenv.2010.04.004
21.
Cheng
,
H.
,
Hayden
,
P.
,
Robins
,
A.
, and
Castro
,
I.
,
2007
, “
Flow Over Cube Arrays of Different Packing Densities
,”
J. Wind Eng. Ind. Aerodyn.
,
95
(
8
), pp.
715
740
.10.1016/j.jweia.2007.01.004
22.
Castro
,
I.
,
Cheng
,
H.
, and
Reynolds
,
R.
,
2006
, “
Turbulence Over Urban-Type Roughness: Deductions From Wind Tunnel Measurements
,”
Boundary-Layer Meteorol.
,
118
(
1
), pp.
109
131
.10.1007/s10546-005-5747-7
23.
Coceal
,
O.
,
Thomas
,
T.
,
Castro
,
I.
, and
Belcher
,
S.
,
2006
, “
Mean Flow and Turbulence Statistics Over Groups of Urban-Like Cubical Obstacles
,”
Boundary-Layer Meteorol.
,
121
(
3
), pp.
491
519
.10.1007/s10546-006-9076-2
24.
Leonardi
,
S.
, and
Castro
,
I.
,
2010
, “
Channel Flow Over Large Cube Roughness: A Direct Numerical Simulation Study
,”
J. Fluid Mech.
,
651
, pp.
519
539
.10.1017/S002211200999423X
25.
Xu
,
H.
,
Altland
,
S.
,
Yang
,
X.
, and
Kunz
,
K.
,
2021
, “
Flow Over Closely Packed Cubical Roughness
,”
J. Fluid Mech.
,
920
, pp.
1
24
.10.1017/jfm.2021.456
26.
Inagaki
,
A.
,
Castillo
,
M.
,
Yamashita
,
Y.
,
Kanda
,
M.
, and
Takimoto
,
H.
,
2012
, “
Large-Eddy Simulation of Coherent Flow Structures Within a Cubical Canopy
,”
Boundary-Layer Meteorol.
,
142
(
2
), pp.
207
222
.10.1007/s10546-011-9671-8
27.
Kono
,
T.
,
Tamura
,
T.
, and
Ashie
,
Y.
,
2010
, “
Numerical Investigations of Mean Winds Within Canopies of Regularly Arrayed Cubical Buildings Under Neutral Stability Conditions
,”
Boundary-Layer Meteorol.
,
134
(
1
), pp.
131
155
.10.1007/s10546-009-9434-y
28.
Kanda
,
M.
,
Moriwaki
,
R.
, and
Kasamatsu
,
F.
,
2004
, “
Large-Eddy Simulation of Turbulent Organized Structure Within and Above Explicitly Resolved Cube Arrays
,”
Boundary-Layer Meteorol.
,
112
(
2
), pp.
343
368
.10.1023/B:BOUN.0000027909.40439.7c
29.
Yang
,
X. I. A.
, and
Meneveau
,
C.
,
2016
, “
Large Eddy Simulations and Parameterisation of Roughness Element Orientation and Flow Direction Effects in Rough Wall Boundary Layers
,”
J. Turbul.
,
17
(
11
), pp.
1072
1085
.10.1080/14685248.2016.1215604
30.
Yang
,
X.
,
Xu
,
H.
,
Huang
,
X.
, and
Ge
,
M.
,
2019
, “
Drag Force on Sparsely Packed Cube Arrays
,”
J. Fluid Mech.
,
880
, pp.
992
1019
.10.1017/jfm.2019.726
31.
Yang
,
X.
,
Sadique
,
J.
,
Mittal
,
M.
, and
Meneveau
,
C.
,
2016
, “
Exponential Roughness Layer and Analytical Model for Turbulent Boundary Layer Flow Over Rectangular-Prism Roughness Elements
,”
J. Fluid Mech.
,
789
, pp.
127
165
.10.1017/jfm.2015.687
32.
Li
,
Q.
,
Bou-Zeid
,
E.
,
Anderson
,
W.
,
Grimmond
,
S.
, and
Hultmark
,
M.
,
2016
, “
Quality and Reliability of Les of Convective Scalar Transfer at High Reynolds Numbers
,”
Int. J. Heat Mass Transfer
,
102
, pp.
959
970
.10.1016/j.ijheatmasstransfer.2016.06.093
33.
Li
,
Q.
, and
Bou-Zeid
,
E.
,
2019
, “
Contrasts Between Momentum and Scalar Transport Over Very Rough Surfaces
,”
J. Fluid Mech.
,
880
, pp.
32
58
.10.1017/jfm.2019.687
34.
Xie
,
Z.
, and
Castro
,
I.
,
2006
, “
LES and RANS for Turbulent Flow Over Arrays of Wall-Mounted Obstacles
,”
Flow Turbul. Combust.
,
76
(
3
), pp.
291
312
.10.1007/s10494-006-9018-6
35.
Cheng
,
Y.
,
Lien
,
F.
,
Yee
,
E.
, and
Sinclair
,
R.
,
2003
, “
A Comparison of Large Eddy Simulations With a Standard k–ϵ Reynolds-Averaged Navier–Stokes Model for the Prediction of a Fully Developed Turbulent Flow Over a Matrix of Cubes
,”
J. Wind Eng. Ind. Aerodyn
,.,
91
(
11
), pp.
1301
1328
.10.1016/j.jweia.2003.08.001
36.
Santiago
,
J. L.
,
Coceal
,
O.
,
Martilli
,
A.
, and
Belcher
,
S.
,
2008
, “
Variation of the Sectional Drag Coefficient of a Group of Buildings With Packing Density
,”
Boundary-Layer Meteorol.
,
128
(
3
), pp.
445
457
.10.1007/s10546-008-9294-x
37.
Raupach
,
M.
,
1992
, “
Drag and Drag Partition on Rough Surfaces
,”
Boundary-Layer Meteorol.
,
60
(
4
), pp.
375
325
.10.1007/BF00155203
38.
Yang
,
X.
, and
Ge
,
M.
,
2021
, “
Revisiting Raupach's Flow-Sheltering Paradigm
,”
Boundary-Layer Meteorol.
,
179
(
2
), pp.
313
323
.10.1007/s10546-020-00597-8
39.
Siemens Industries Digital Software
, 2019, “
Star CCM+ v 2019.2 User Manual
,” Siemens.
40.
Menter
,
F.
,
1994
, “
Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications
,”
AIAA J.
,
32
(
8
), pp.
1598
1605
.10.2514/3.12149
41.
Lien
,
F.
,
Chen
,
W.
, and
Leschziner
,
M.
,
1996
, “
Low Reynolds-Number Eddy-Viscosity Modelling Based on Non-Linear Stress-Strain/Vorticity Relations
,”
Eng. Turbul. Modell. Exp.
,
3
, pp.
91
100
.10.1016/B978-0-444-82463-9.50015-0
42.
Launder
,
B.
, and
Shima
,
N.
,
1989
, “
Second-Moment Closure for the Near-Wall Sublayer: Development and Application
,”
AIAA J.
,
27
(
10
), pp.
1319
1325
.10.2514/3.10267
43.
Gibson
,
M.
, and
Launder
,
B.
,
1978
, “
Ground Effects on Pressure Fluctuations in the Atmospheric Boundary Layer
,”
J. Fluid Mech.
,
86
(
3
), pp.
491
511
.10.1017/S0022112078001251
44.
Rodi
,
W.
,
1991
, “
Experience With Two-Layer Models Combining the k-ϵ Model With a One-Equation Model Near the Wall
,”
AIAA
Paper No. 91-0216.10.2514/6.1991-216
45.
Xie
,
Z. T.
, and
Fuka
,
V.
,
2018
, “
A Note on Spatial Averaging and Shear Stresses Within Urban Canopies
,”
Boundary-Layer Meteorol.
,
167
(
1
), pp.
171
179
.10.1007/s10546-017-0321-7
46.
Kuwata
,
Y.
, and
Suga
,
K.
,
2013
, “
Modelling Turbulence Around and Inside Porous Media Based on the Second Moment Closure
,”
Int. J. Heat Fluid Flows
,
43
, pp.
35
51
.10.1016/j.ijheatfluidflow.2013.03.001
47.
Aupoix
,
B.
,
2016
, “
Revisiting the Discrete Element Method for Predictions of Flows Over Rough Surfaces
,”
ASME J. Fluids Eng.
,
138
(
3
), p.
031205
.10.1115/1.4031558
48.
Taylor
,
R.
,
Coleman
,
H.
, and
Hodge
,
B.
,
1985
, “
Prediction of Turbulent Rough-Wall Skin Friction Using a Discrete Element Approach
,”
ASME J. Fluids Eng.
,
107
(
2
), pp.
251
257
.10.1115/1.3242469
49.
Kato
,
M.
, and
Launder
,
B.
,
1993
, “
The Modelling of Turbulent Flow Around Stationary and Vibrating Square Cylinders
,”
Ninth Symposium on Turbulent Shear Flows
, Kyoto, Japan, Aug. 16–18, pp.
1
6
.