Direct numerical simulation of turbulent channel flow over a surrogate for Nikuradse-type roughness

A tiled approach to rough surface simulation is used to explore the full range of roughness Reynolds numbers, from the limiting case of hydrodynamic smoothness up to fully rough conditions. The surface is based on a scan of a standard grit-blasted comparator, subsequently low-pass filtered and made spatially periodic. High roughness Reynolds numbers are obtained by increasing the friction Reynolds number of the direct numerical simulations, whereas low roughness Reynolds numbers are obtained by scaling the surface down and tiling to maintain a constant domain size. In both cases, computational requirements on box size, resolution in wall units and resolution per minimum wavelength of the rough surface are maintained. The resulting roughness function behaviour replicates to good accuracy the experiments of Nikuradse (1933 VDI-Forschungsheft, vol. 361), suggesting that the processed grit-blasted surface can serve as a surrogate for his sand-grain roughness, the precise structure of which is undocumented. The present simulations also document a monotonic departure from hydrodynamic smooth-wall results, which is fitted with a geometric relation, the exponent of which is found to be inconsistent with both the Colebrook formula and an earlier theoretical argument based on low-Reynolds-number drag relations.


Introduction
The equivalent sand-grain roughness height, k s , was first defined by Schlichting (1936) as the size of the sand grains from the experiments of Nikuradse (1933) that gave effectively the same frictional resistance as the roughness under consideration. This parameter has subsequently served as a 'universal currency of exchange' (Bradshaw 2000) in roughness studies, and many researchers have aimed at its prediction. Nikuradse (1933) was the first to classify the flow into three regimes, which, depending on his k + s , had the following ranges: k + s < 4 for the hydraulically smooth regime, 4 k + s 70 for the transitionally rough regime and k + s > 70 for the fully rough regime (where k + s = k s u τ /ν is the equivalent sand-grain roughness height in wall units, with u τ the friction velocity and ν the kinematic viscosity). According to Jiménez (2004), the k s of a given rough surface guarantees collapse of all roughness functions in the fully rough regime. Outside this regime, the problem becomes more difficult. Nikuradse (1933) remains the most complete study of the transitionally rough regime, while subsequent work (Ligrani & Moffat 1986;Shockling, Allen & Smits 2006;Schultz & Flack 2007) has produced evidence for non-universal behaviour. Other recent roughness studies covering a wide Reynolds number range include Flack et al. (2016), where the skin-friction behaviour of surfaces blasted with different grades of grit was studied, and Barros, Schultz & Flack (2017), where skin-friction studies were conducted on mathematically generated irregular roughness.
Nikuradse's original surfaces are long-since lost, and his method of coating sieved sand grains with lacquer makes the surfaces difficult to reproduce. In addition, there have been comments on his idiosyncratic experimental technique (Hager & Liiv 2008) which also helped to motivate the present study based on a direct numerical simulation (DNS) approach. Relatively few DNS studies over a wide range of k + s of roughness exist, mainly due to limitations on the computational resources that are required to conduct high-Reynolds-number simulations necessary to achieve fully rough flow. For sinusoidal pipe roughness, Chan et al. (2015) used body conforming grids to study the effect of height and wavelength. An extensive DNS study was published recently by , where a comparison of the near-wall flow for two realistic irregular rough surfaces, a graphite surface and a grit-blasted surface (which is also being studied in the current paper), was made for a wide Reynolds number range. The fully rough regime was reached for both surfaces. This study, however, did not provide DNS data for the lower part of the transitionally rough regime. This regime is also interesting due to previous speculations about the behaviour as k + s → 0. In particular, Bradshaw (2000) objected to the lower limit of k + s = 4 given by Nikuradse (1933) and suggested that the change must be gradual, for example U + ∼ (k + s ) α (where U + is the roughness function), with potentially α = 1 from the Colebrook curve or α = 2 from arguments put forward by Bradshaw (2000) based on the Oseen solutions for drag on an isolated sphere.
The goals of the present study are to investigate the rough-wall flow over an irregular industrial grit-blasted surface, obtained from a standard comparator sample, in the entire Reynolds number range using DNS, and make comparisons with the benchmark experimental data of Nikuradse (1933). A three-step methodology, consisting of surface scanning, pre-processing (which includes surface scaling and Fourier filtering) and DNS of incompressible rough channel flow, is used to conduct the simulations. The methodology is described in detail by Busse, Lützner & Sandham (2015) and . The surface is scanned using a variable focus microscope, after which it is scaled and then filtered using a low-pass Fourier filter. Filtering removes roughness length scales below a user-defined threshold to obtain a smoothly varying topography (figure 1) and makes the surface periodic in the streamwise and spanwise directions. Figure 1(c) shows an example of a weighted streamwise roughness elevation spectrum, k x G(k x ), before filtering, as a function of the streamwise wavenumber, k x , showing the extent of the streamwise filter cutoff Direct numerical simulation of a surrogate for Nikuradse-type roughness wavenumber, k c , and inverse of the streamwise grid spacing, 1/ x. An exhaustive list of surface topographical properties characterising the filtered sample is also provided by Thakkar et al. (2017), where the grit-blasted sample is denoted as s8.

Problem formulation
The streamwise, spanwise and wall-normal directions are denoted by x, y and z respectively. The friction Reynolds number, Re τ = δu τ /ν, where δ is the mean channel half-height, varies in the range 180 Re τ 720, whereas the roughness Reynolds number, k + = ku τ /ν = (k/δ)Re τ , varies in the range 3.75 k + 120 (where k is the mean peak-to-valley height and a '+' superscript indicates wall units), as shown in table 1. All cases for k + 30 are scaled to a roughness height of k/δ = 1/6, which has been previously shown to be small enough to collapse the mean flow in the outer layer (Busse et al. 2015), as well as being directly relevant to flow over gravel beds (Mohajeri et al. 2015) and urban environments (Orlandi & Leonardi 2006). The computational domain size scales with the roughness height, and k/δ = 1/6 leads to streamwise, L x , and spanwise, L y , domain extents of 5.63δ and 2.815δ respectively, where L y = L x /2. For a typical industrial rough surface, Busse et al. (2015) showed that a minimum streamwise domain extent of approximately 5δ was sufficient to obtain domain-size-independent mean flow and Reynolds stress statistics. A maximum streamwise wavenumber of 24 is obtained after Fourier filtering.  Here, n x , n y and n z are the numbers of grid cells in the streamwise, spanwise and wall-normal directions respectively. height is used to set a datum for the surface at z = 0 for the lower wall, and similarly for the upper wall at z = 2.
In principle, there are several ways to achieve a particular k + by varying k/δ and Re τ . The U + results from a study at three different values of Re τ , namely 120, 180 and 360, with corresponding k/δ = 1/4, 1/6 and 1/12, and hence the same k + = 30, are shown in figure 2(a). The corresponding domain sizes are 8.445δ × 4.223δ, 5.63δ × 2.815δ and 2.815δ × 1.408δ respectively, and all three cases have uniform streamwise and spanwise grid spacings, with x + = y + 5, and a stretched grid in the wallnormal direction, with z + min ≈ 0.67 (in the roughness layer, i.e. from the lowest valley to the highest peak of the sample) and z + max < 5 (at the channel centre). It is known that different flows having the same k + should produce the same U + (Krogstad & Antonia 1999). The results show deviations from a constant U + only for the Re τ = 120 case. A similar study was conducted by Chan et al. (2015), and they chose Re τ = 180 as their lowest Reynolds number for DNS of sinusoidal pipe roughness.
To increase k + above 30, Re τ is simply increased while maintaining the same k/δ. To obtain k + < 30, however, decreasing the Reynolds number below Re τ = 180 may lead to low-Reynolds-number effects. Hence, instead of reducing the Reynolds number below Re τ = 180, samples for k + = 15, 7.5 and 3.75 are obtained by shrinking and tiling the Re τ = 180 (k + = 30) case. For example, to obtain the k + = 15 case, the k + = 30 case is scaled down by a factor of two, which reduces the roughness height from k/δ = 1/6 to k/δ = 1/12, giving k + = 15 at Re τ = 180. The topography within each tile is identical to the k + = 30 reference case. The shrunk sample is then tiled (or replicated) in a 2 × 2 fashion. Tiling maintains the required minimum streamwise domain size of approximately 5δ (as mentioned in the first paragraph of this section) and automatically satisfies the periodic boundary conditions without discontinuities. The same procedure is applied for the k + = 7.5 and 3.75 cases, considering the appropriate scaling and tiling of the k + = 30 case, i.e. one-fourth and 4 × 4 tiles for k + = 7.5, and one-eighth and 8 × 8 tiles for k + = 3.75. The sample for k + = 15 is shown in plan view in figure 2(b). Thakkar (2017) shows how tile-to-tile correlations within the flow are confined to the immediate roughness layer.
The overall simulation parameters are provided in table 1. All cases have a domain size of 5.63δ × 2.815δ, a uniform streamwise and spanwise grid, with x + = y + 5, and a stretched wall-normal grid, with z + min ≈ 0.67 (in the roughness layer, i.e. from the lowest valley to the highest peak of the sample) and z + max < 5 (at the channel centre). Moreover, the grid contains a minimum of 12 points per smallest wavelength of the roughness.

Results
Values of U + are computed as the difference between the mean centreline velocities, U c /u τ , of the rough-and smooth-wall cases at the same Re τ (Busse et al. 2015). Reference data for the smooth-wall simulations are shown in table 2, where all cases have a domain size of 12δ × 6δ and z + min ≈ 0.67. As the roughness decreases (lower k + ) and the flow approaches hydrodynamically smooth conditions, the statistics may become sensitive to the domain size. To assess this issue, for the three smallest values of k + , namely 15, 7.5 and 3.75, U + was also computed based on a smooth-wall simulation with the same domain size as the rough sample, i.e. 5.63δ × 2.815δ, and Re τ = 180. Other parameters for this simulation were n x × n y × n z = 128 × 128 × 224, with x + = 7.92, y + = 3.96, z + min = 0.67, z + max = 4.70, U c /u τ = 18.46 and the streamwise bulk velocity, U b /u τ = 15.83. Compared with the smooth-wall data in table 2, U c /u τ is only slightly different, and hence the value in table 2 was used to calculate U + .   is observed. The profile for k + = 3.75 is quite close to the smooth-wall profile and has a very small value of U + = 0.35. The velocity defect profiles, shown in figure 3(b), collapse for z/δ 0.1, indicating that they have already achieved a Reynolds-number-independent shape in the outer regions of the flow, despite the relatively high k/δ for most cases. Differences are observed only in the immediate vicinity of the rough wall.

Roughness function
From table 1, it is clear that the two highest values of k + , namely 90 and 120, are in the fully rough regime as their U b /u τ values are approximately constant (the difference being approximately 0.1 %), which is an indication that the effective friction factor, 2(u τ /U b ) 2 , is independent of k + and hence Re τ . Moreover, the fact that their U + values are greater than 7 further confirms their presence in the fully rough regime .
Figure 4(a) shows the plot of U + against k + on semilogarithmic axes. Figure 4(b) shows the piecewise linear curve obtained from the data of Nikuradse (1933), the Colebrook relation, given by U + = (1/κ) ln(1 + 0.3k + s ), and the fully rough asymptote, given by Direct numerical simulation of a surrogate for Nikuradse-type roughness that in the fully rough regime, k + s collapses the roughness functions for all types of roughness. Thus, U + values for the grit-blasted sample at higher k + are matched with the expected fully rough behaviour, to obtain k + s ≈ 0.87k + , and plotted in figure 4(b).
From figure 4(b), it is also observed that the current data follow the Nikuradse (1933) curve quite closely throughout the Reynolds number range. This shows that, for the first time, roughness closely resembling the sand-grain roughness of Nikuradse (1933) has been simulated over the entire Reynolds number range, from the hydrodynamically smooth to the fully rough regime. Experiments replicating those of Nikuradse (1933) have appeared only in a conference paper by Nakato et al. (1985), while the agreement of the current DNS data means that we have a surrogate of the type of surface that Nikuradse used. The similarities may extend to the smallest roughness scales. Nikuradse (1933) applied his lacquer twice, once to stick the sand grains to the pipe surface and then again 'in order to obtain a better adherence of the grains'. Nikuradse (1933) also mentioned that 'this lacquer formed a direct coating on the wall and also a covering on the grains no thicker than the penetration of these grains into the lacquer coating of the wall'. This statement means that effectively the sand grains were completely coated in lacquer, which would have led to elimination of the extremely small length scales, which is exactly what Fourier filtering of the grit-blasted surface does. To gain a visual understanding of the preceding statement, the reader is referred to figure 4 in Nikuradse (1933), which shows a photograph of his sand grains.
Experiments on honed roughness carried out by Shockling et al. (2006) and Schultz & Flack (2007), and grit-blasted surfaces studied by Flack et al. (2016) also showed similar trends to those of Nikuradse (1933) over the Reynolds number range, consistent with the present results. However, this tends to contradict the observation of Bradshaw (2000), who suggested that realistic roughness should exhibit more Colebrook-type than Nikuradse-type behaviour due to the large range of length scales present compared with Nikuradse's sand grains. Rather, the absence of modern studies that reproduce Colebrook-type behaviour would appear to be a more pressing concern.
3.3. Nikuradse's scaling Further comparisons are made between the two data sets by plotting the DNS data in Nikuradse's scaling, described by his log-region velocity profile parameter, A, computed by Nikuradse (1933) as For a channel, the pipe radius, r, is replaced by δ and the results are plotted in figure 5. The range of k + s is divided into three regimes (denoted by black dashed lines) by Nikuradse: the hydraulically smooth regime for log 10 (k + s ) 0.55, the transitionally rough regime for 0.55 < log 10 (k + s ) 1.83 and the fully rough regime for log 10 (k + s ) > 1.83. The transitionally rough regime is further divided into three regions (denoted by grey dashed lines in figure 5), depending on whether the value of A increases (0.55 < log 10 (k + s ) 0.85), remains mostly constant (0.85 < log 10 (k + s ) 1.15) or decreases (1.15 < log 10 (k + s ) 1.83). Excellent agreement between the current DNS and the experimental data is observed. Similar comparisons were made by Ligrani & Moffat (1986)  as with the current data. In particular, the honed pipe roughness of Shockling et al. (2006) did not exhibit as pronounced an inflectional behaviour as shown by the data of Nikuradse (1933) or Ligrani & Moffat (1986), and reached the fully rough regime at comparatively lower k + s ≈ 40 (or log 10 (k + s ) ≈ 1.6).
3.4. Data characterisation Least-square curve fits to the DNS data allow further comparisons with experiment and also allow the asymptotic limit of small roughness to be studied in more detail. Although there are fewer data points in the DNS, these points are known to be high precision, due to the extensive validations and the long run times to collect good statistics, as described by Busse et al. (2015) and Thakkar et al. (2017). For fitting purposes, the data are classified into three sections: a lower transitionally rough regime (using the first three data points, i.e. 0 k + s 13.05), an upper transitionally rough regime (13.05 < k + s 78.3) and the fully rough regime (k + s > 78.3). We also use the fact that U + = 0 at k + s = 0 (smooth-wall condition) provides a reference point for hydrodynamically smooth flow. Fit quality is determined by the root mean square error, σ , and the coefficient of determination, R 2 . Curves are fitted using polynomial, exponential or power laws, and those fits with the lowest values of σ , and hence highest values of R 2 , are discussed in the following. Figure 6 shows the curve fits and table 3 shows the fit quality parameters in the various regimes. The equations of the three fits are given as 0 k + s 13.05 : Data points in the lower transitionally rough regime exhibit a power-law behaviour, and the curve fit, equation (3.2), satisfies the smooth-wall condition, U + = 0 at k + s = 0. In the upper transitionally rough regime, a log-law behaviour is observed. The fit quality in both of the above regimes is excellent, as observed from table 3. The data of Nikuradse (1933) in the range 14.13 k + s 67.61 gave U + = 1.79 ln(k + s ) − 6. Equation (3.3) lies in a similar range and, although the data agree with Nikuradse (1933) visually (figure 4b), the difference in the fit is probably due to the limited number of points in the current data set. In the fully rough regime, the last two data points are simply connected and hence the fit shows no root mean square error and R 2 = 1 (table 3). The slope of the fully rough asymptote is 1/κ = 1/0.4 = 2.5 837 R1-8    > 70), the slope of the curve being 2.50. Compared with these slope values, the difference in slope of the fully rough regime fit, equation (3.4), is approximately 11 %. This difference could be reduced by having more data points in the fully rough regime.
The curve fit in the lower transitionally rough regime, equation (3.2), is of particular interest following an earlier speculation of Bradshaw (2000). He noted that there was no lower limit to k + s below which the flow would be hydrodynamically smooth, as might be misinterpreted from Nikuradse's curve fits. Interestingly, he also speculated on the form the departure from hydrodynamically smooth conditions might take, suggesting a quadratic law based on an argument involving the Oseen correction to Stokes flow, predicting an exponent of 2 in equation (3.2). The numerical value of 1.37 is somewhat sensitive to the fit methodology, given the limited number of points available for analysis, but is clearly not in agreement with the quadratic law, which may be recovered only much closer to the wall. In this context, it would be interesting to repeat the present DNS study for other samples, since it is possible that the initial departure from smooth-wall flow depends on the particular surface topology, given that the roughness height at which roughness effects are first observed and the length of the transitionally rough regime have been proved to depend on the surface topology (Flack, Schultz & Rose 2012).

Conclusions
The behaviour of an irregular realistic grit-blasted surface, in the entire Reynolds number range from the transitionally rough to the fully rough regime, has been studied 837 R1-9 Downloaded from https://www.cambridge.org/core. using DNS, with a tiling approach employed to simulate cases with low roughness Reynolds number, which presents a different set of challenges to cases with high Reynolds number. Excellent agreement of all cases with the behaviour observed by Nikuradse is found, both in terms of the roughness function, U + , and Nikuradse's own scaling parameter, A. One would expect other sand-grain-like surfaces to also follow Nikuradse's trend, so the current surface is not unique. Nevertheless, the results are useful in providing a digital representation of a rough surface that shows 'Nikuradse' hydrodynamic roughness behaviour. In addition, the simulations identify theoretical issues relating to the approach to hydrodynamic smoothness, since neither the Colebrook curve nor Bradshaw's quadratic law is close to the DNS results.