Surface correlations of hydrodynamic drag for transitionally rough engineering surfaces

Rough surfaces are usually characterised by a single equivalent sand-grain roughness height scale that typically needs to be determined from laboratory experiments. Recently this method has been complemented by a direct numerical simulation approach, whereby representative surfaces can be scanned and the roughness eﬀects computed over a range of Reynolds number. This development raises the prospect over the coming years of having enough data for diﬀerent types of rough surfaces to be able to relate surface characteristics to roughness eﬀects, such as the roughness function that quantiﬁes the downward displacement of the logarithmic law of the wall. In the present contribution, we use simulation data for 17 irregular surfaces at the same friction Reynolds number, for which they are in the transitionally rough regime. All surfaces are scaled to the same physical roughness height. Mean streamwise velocity proﬁles show a wide range of roughness function values, while the velocity defect proﬁles show a good collapse. Proﬁle peaks of the turbulent kinetic energy also vary depending on the surface. We then consider which surface properties are important and how new properties can be incorporated into an empirical model, the accuracy of which can then be tested. Optimised models with several roughness parameters are systematically developed for the roughness function and proﬁle peak turbulent kinetic energy. In determining the roughness function, besides the known parameters of solidity (or frontal area ratio) and skewness, it is shown that the streamwise correlation length and the root-mean-square roughness height are also signiﬁcant. The peak turbulent kinetic energy is determined by the skewness and root-mean-square roughness height, along with the mean forward-facing surface angle and spanwise eﬀective slope. The results suggest feasibility of relating rough-wall ﬂow properties throughout the range from hydrodynamically smooth to fully-rough to surface parameters.


Introduction
Rough surfaces are encountered in a large number of applications; from roughness in conjunction with industrial heat exchangers [1], turbomachinery [2,3], ship propellers and hulls [4][5][6] to roughness induced by plant canopies and vertical structures in an urban environment pertaining to atmospheric flows [7,8]. The drag associated with the transport of goods by ships is of particular interest, given the associated emissions. According to [6], any solid surface exposed to the marine environment will be affected by fouling. Marine fouling, which is caused by the accumulation of organic molecules, microorganisms, plants and animals on a body submerged in the water [4], leads to an increase in roughness of the hull and hence its hydrodynamic drag. The drag penalty causes a decrease in ship speed and maneuverability and an increase in fuel consumption. Propeller fouling, al-though a small part of the fouling on the marine vehicle, is also important from the point of view of increased friction and fuel consumption which in turn hampers performance. Due to extended periods of service, turbines, compressors and other turbomachinery components are adversely affected by roughness since their surface quality degrades due to phenomena such as erosion, corrosion and deposition [2]. Heat exchangers utilise roughness to improve their efficiency [1] as the increase in wall friction causes an increase in the wall shear stress which enhances the heat transfer rate.
Many industrial surface finishing processes produce materials that are classified as rough. Examples of such processes include grinding, shotblasting, spark-erosion, casting etc. Understandably, the connection between surface topology and drag is a fundamental topic in fluid dynamics. Previous work has been mostly limited to numerical and experimental studies on regular rough surfaces made from systematic arrangements of cubes, bars, cylinders, rods, spheres etc. possessing a small number of characteristic length scales and whose surface properties could be easily evaluated. The objective of the current study is to conduct direct numerical simulations (DNS) of a range of well-characterised, scanned irregular rough surfaces seen in practical applications and methodically relate their surface parameters to various flow properties.
The primary effect of roughness is an increase in the surface friction compared to a smooth wall, which is seen as a downward shift in the mean streamwise velocity profile when plotted in wall-units. This shift is quantified by the roughness function, ∆U + , also known as the roughness effect. Based on the smooth-wall log-law profile, this velocity deficit can be represented as where '+' superscripts indicate wall-units, z + and k + are the wall-normal distance and roughness height respectively in wall-units, κ ≈ 0.4 is the von Kármán constant and A = 5.5 is the additive constant. Within the fully-rough regime, ∆U + can be empirically related to the equivalent sand-grain roughness, k s,eq using an equation of the form which can be derived from equations (3) and (4) in Jiménez [9]. It must be noted that k s,eq is not known a priori and must be determined experimentally or numerically using DNS by performing a Reynolds number sweep from the transitionally rough into the fully-rough regime. The data is then matched in the fully-rough regime with a standard reference curve such as the Colebrook universal interpolation formula (refer Jiménez [9]), given by ∆U + = 1 κ ln(1 + 0.26k + s,eq ).
(3) Figure 1 displays ∆U + against k + s,eq on semilogarithmic axes for the Colebrook relation as well as the fully-rough asymptote. The piecewise linear curve obtained from the uniform sand pipe-flow experiments of Nikurade [10], which is regarded as one of the benchmark studies in roughness, is also shown.
Schlichting [11] performed a series of experiments on regular rough surfaces that included staggered arrangements of spheres, spherical segments, cones and angular plates. His experiments were conducted in the fully rough regime i.e. the regime where the friction factor is independent of Reynolds number. The experiments were carried out in rectangular channels at Re = ud/ν = 4.3 × 10 5 where u is the mean velocity of the flow, d is the channel hydraulic diameter and ν is the fluid kinematic viscosity. One of the most important objectives of Schlichting's study was to develop a model to predict the surface friction for rough surfaces similar to those used in his experiments but at other Reynolds numbers and roughness ratios, k/r h , where k is the absolute height of the roughness elements from the plate on which they were mounted and r h is the hydraulic radius. This involved determining an equivalent sand-grain roughness, k s,eq which was the equivalent size of sand grains as used in the experiments of Nikuradse [10] and which had the same resistance as the geometry under consideration. Schlichting proposed that the surface resistance depended not only on the relative roughness, r h /k but also on the roughness density, S f /S, where S f is the total projected area of the roughness elements on a plane normal to the direction of the flow (or the frontal area of the roughness elements) and S is the surface area of the plate on which the roughness elements are mounted. He also proposed a resistance coefficient for rough surfaces as C f = 2W r /(ρu 2 k S f ), where W r = W − W g is the resistance due to the roughness elements alone, W is the total resistance of the rough plate, W g is the resistance of the smooth areas between the roughness elements and u k is the velocity at a distance from the wall y = k. It was found that C f was independent of S f /S for small values of roughness density and decreased rapidly for large values of roughness density.
The equivalent sand-grain roughness of Nikuradse [10], k s,eq , has subsequently become the universal currency of exchange, as mentioned by Bradshaw [12], in the study of rough surfaces and many researchers have aimed at its prediction using correlations to surface parameters. Sigal & Danberg [13,14] conducted a study to determine a suitable geometric correlation relating to the roughness density effect. Their relation was based on a database of results obtained by Schlichting's experiments [11] and twelve other regular roughness studies (see [13,14] and references therein). Their new roughness density parameter, Λ s was given as where S is the planform area of the corresponding smooth surface, S f is the total frontal area of all roughness elements (S and S f are equivalent to those used in Schlichting's [11] studies), A f is the frontal area of a single roughness element and A s is the wetted area of a single roughness element. where k is the absolute height of the roughness elements from the surface on which they are mounted (which is equivalent to that used in Schlichting's [11] studies). The scarcity of data for three-dimensional roughness prevented the authors from developing a completely general model and as such their parameter is known to be better suited for two-dimensional roughness [13,14]. Through a series of channel flow experiments on smooth, patterned rough and completely rough surfaces van Rij et al. [15] proposed a more generalized form of the Sigal-Danberg parameter which could be applied to three-dimensional irregular roughness. In case of irregular three-dimensional roughness, A f /A s is replaced by S f /S w , the ratio of the total frontal area to the total wetted area for all the roughness elements. Hence, the modified version of the parameter was given as where S w is the total area of all roughness elements wetted by the flow. A modified equation for the equivalent sand-grain roughness was also proposed as k s,eq k = where k is the average roughness element height of the surface (S a in the notation of the present paper, refer Appendix A for definition).
The work of Bons [16] is important in the field of turbomachinery roughness correlations as he defined a streamwise forward-facing surface angle, α i of roughness elements (refer Appendix A for definition). Based on experimental data for six types of turbine blade roughness, he also proposed an associated correlation as where k is again the average roughness element height of the surface (S a in the notation of the present paper) and α is the average streamwise forward-facing surface angle (in degrees). The review by Flack & Schultz [17] on previously proposed roughness correlations in various roughness regimes covered many experiments on different types of regular roughness including mesh, spheres, pyramids and square bars etc. and irregular roughness including different types of sandpaper, honed surfaces, uniform sand and turbine blades subject to pitting and corrosion etc. They suggested that the correlations proposed in the past were useful only for a subset of rough surfaces and could not be applied to roughness in general, especially irregular roughness. Hence, their aim was to propose a suitable new correlation that could be used more generally and that could be applied to a wider selection of irregular and three-dimensional rough surfaces and hence provide a method to enable drag prediction based solely on the surface topography. They considered flows mainly in the fully rough regime due to the availability of a large quantity of experimental results. A statistical analysis conducted by the authors on various roughness scaling parameters indicated that the root-mean-square (rms) roughness height, S q and the skewness, S sk of the surface elevation probability density function correlated strongly with k s,eq . The proposed correlation was given by It accurately predicted k s,eq values for most of the surfaces considered by the authors, although complete generality was not achieved with it.
There have been other notable contributions in the study of hydrodynamic drag prediction using surface property correlations. Musker [18] proposed a new relation for an effective roughness height and correlated it with the roughness function using seven surfaces representative of a variety of ship-hull roughness. The surface geometric properties included in the relation were the rms roughness height, surface skewness, kurtosis and the average slope of roughness elements. Waigh & Kind [19] formulated relations for the roughness effect, based on 16 experimental studies comprising of various types of regular roughness geometries of differing shape and distribution in the fully-rough regime. Their relations included a roughness spacing parameter, ratio of the roughness height to spanwise length for a single element and ratio of the wetted area to frontal area of a single roughness element. In the field of urban roughness and the atmospheric boundary layer, Wieringa [20] and Grimmond & Oke [21] have provided a number of empirical correlations.
It must be noted that all studies mentioned above were experimental. More recently, it has been possible to conduct numerical simulations that complement the experimental database. Yuan & Piomelli [22] conducted studies to estimate and predict the roughness function and k s,eq on realistic surfaces. They carried out large-eddy simulations in turbulent open-channel flows over sand-grain roughness and realistic roughness replicated from hydraulic turbine blades. Both transitionally rough and fully rough regimes were covered by considering different roughness heights at two Re τ = 400 and 1000. They evaluated the performance of three existing correlations, proposed by van Rij et al. [15], Bons [16], and Flack & Schultz [17], to predict k s,eq . These correlations have already been displayed in equations (6), (7) and (8) respectively above. Data collapse was obtained for the first two correlations, which are slope-based whereas the third correlation, which is moment based, showed data scatter. The reason for this scatter was that moments did not contain slope information and thus did not scale with k s,eq in cases where surface slope was an important parameter. Surface slope was influential in their studies because their surfaces were in the waviness regime (as described by [23]), which meant that there was a dependence of ∆U + on the effective slope, ES (refer Appendix A for definition).
The objectives of the above-mentioned studies were to characterise irregular roughness purely on the basis of geometrical considerations. The present study has similar aims but based on DNS data. Direct numerical simulations of 17 industrially relevant irregular rough surfaces scaled to the same physical roughness height at the same friction Reynolds number are conducted. The relatively large size of the surface database means a wide range of rough surfaces seen in practical applications with a broad spectrum of topographical properties has been considered. Based on the simulation data, surface topographical properties influencing the hydrodynamic drag and turbulent kinetic energy are systematically determined. A description of the numerical methodology, computational geometry, boundary conditions and meshing criteria is given in Section 2. Section 3 gives a brief description of the rough surface samples studied and their corresponding surface properties. Section 4 shows statistical results for the mean streamwise velocity, velocity defect profiles and turbulent kinetic energy. Section 5 provides a methodical approach to determine which surface properties are influential in determining the hydrodynamic drag and turbulent kinetic energy. Finally, Section 6 describes the conclusions from this study.

Numerical methodology
A three-step methodology, as described by Busse et al. [24], is used to conduct the simulations; surface data acquisition, data filtering and direct numerical simulations.

Surface data acquisition and data filtering
The surface data for all samples has been obtained using an Alicona Infinite Focus microscope which measures the surface height by focus variation. The data is obtained in the form of a height map of surface z coordinates (roughness heights from the mean reference plane) as a function of its x and y (streamwise and spanwise) coordinates.
For all rough surfaces, the surface scan is obtained as a discrete height map on a regular cartesian grid in x and y; x = 0, ∆s, 2∆s, . . . , (M − 1)∆s and y = 0, ∆s, 2∆s, . . . , (N − 1)∆s where ∆s is the spacing of the measurement points as obtained during the scan and M and N are the number of data points in the streamwise and spanwise directions respectively. The sample for simulation is obtained as a smaller sub-section of the scan. In order to select a sample which is representative of the rough surface under consideration, the physical size of the subsection is initially determined based on a visual inspection of the surface scan. The sub-section, which maintains a fixed 2 : 1 (streamwise to spanwise domain width) aspect ratio, must be chosen to retain sufficient roughness features, but taking into account the computational cost. After selection, the sub-section is checked, and if necessary re-selected, so that it maintains adequate roughness correlation lengths within the streamwise computational domain and (since the simulations are conducted in channels) adequate domain size in terms of channel half-heights. The smallest streamwise domain lengths have an extent of approximately 5 times the mean channel half-height. As was shown by Busse et al. [24], this is sufficient to obtain domain size independent rough-wall mean flow and Reynolds stress statistics. The described technique can be adopted as most surfaces in the current study exhibit a homogenous distribution of roughness features. The location of the subsection on the scan is determined based on rms errors in roughness heights at the lateral boundaries. In order to minimize non-periodicity in the lateral boundaries, the sub-section with least rms errors in roughness heights between its streamwise boundaries and between its spanwise boundaries must be selected. A combined error for the streamwise and spanwise boundaries is used to determine the best sub-section. The sample in its raw form is unsuitable for simulation and the surface data needs to be filtered. Filtering is done in Fourier space and is essentially a smoothing step. It needs to be done for the following reasons: firstly, the surface scan usually contains a finite amount of measurement noise which is typically on small spatial scales [25]. It is essential to remove this noise. Secondly, due to computational constraints, it is not possible to resolve all the length scales of roughness. From an aerodynamic perspective, the smallest roughness scales are usually not relevant [26] and according to Jiménez [9], the effect of roughness is known to be dominated by the largest features of a rough surface. Filtering removes the smallest scales which are below a user-defined threshold. Thirdly, periodic boundary conditions are used in the streamwise and spanwise directions to reduce computational cost and perform efficient simulation in reasonably small computational domains. Filtering makes the rough surface sample periodic. In the case of non-periodic boundary conditions, very large domains would be required in order to ensure independence of the flow parameters from the inlet and outlet boundary conditions, which would significantly increase computational cost. The surface data is hence filtered using a low-pass filter to obtain an approximate model of the 3D surface topography. Details of the filtering process are described in Busse et al. [24]. Figure 2 shows an example of a surface sample before and after filtering.
It is essential to choose an appropriate value for the cut-off wavenumber k c , all wavenumbers above which are filtered out. If k c is too low, the filtered data will be too smooth and will not be an accurate representation of the original data. If it is too high, a lot of small and aerodynamically irrelevant scales would be present in the data which would significantly increase the computational cost. The value of k c depends to a great extent on the topography of the rough surface and hence no general recommendations can be given. However, studies conducted by Busse et al. [24] on one of the samples considered in the present work showed that a difference of 8% between the filtered and unfiltered values of the average and rms roughness height, S a and S q , retained most of the large scale surface topography. The same criterion is used in the present work to determine k c whose value is adjusted accordingly depending on the sample. Once k c and the domain size are specified, the periodic sample is a precisely-defined representation of the original surface and can be used together with DNS in a rigorous manner.

Geometry, boundary conditions and meshing criteria for DNS
The rough surface samples are used as no-slip wall boundaries in incompressible turbulent channel flow. The streamwise, spanwise and wall-normal directions in the computational domain are denoted by x, y and z respectively, with corresponding domain lengths L x , L y and L z . Considering the cut-off wavenumber criterion mentioned in the previous section in conjunction with the streamwise domain length, the maximum streamwise wavenumber is given by k c L x . The samples have an aspect ratio of 2 : 1, which means L x = 2L y . The rough surface on the upper boundary corresponds to a mirror image of that on the lower boundary but translated by L x /2 and L y /2 in the streamwise and spanwise directions respectively. The mean surface height is set as the mean reference plane, z = 0 at the bottom boundary and z = 2δ at the top boundary, where δ is the channel half-height. The channel height of 2δ is measured as the distance between the bottom and top mean reference planes. The domain length in the wall-normal direction, L z , is slightly larger than 2δ to take into account the height of the roughness features. Figure 3 shows a schematic representation of the computational domain. The friction Reynolds number for this study is Re τ = δu τ /ν = 180, where u τ is the friction velocity of the fluid, for which the flow is in the transitionally rough regime and where DNS is feasible for large number of samples. All rough surface samples are scaled to the same roughness height, k, defined for this study by the mean peak-to-valley height, S z,5×5 (refer Appendix A for definition), such that k = δ/6. It has been recommended that in order to study universal roughness behaviour, k should be small compared to δ. Jiménez [9] recommends δ/k in excess of 40. In order to achieve a significant roughness effect for δ/k > 40, very high Re τ , in excess of 1000, would be required. This in turn would lead to extremely dense meshes as the small scales of motion, especially close to the rough walls would need to be resolved. These factors lead to a prohibitively high computational cost. Hence, δ/k = 6 is used, which leads to a clear roughness effect at Re τ = 180. This is discussed in Section 4. Also discussed in Section 4 is the effect of the relatively low δ/k within the context of outer-layer similarity.
Uniform grid spacing is used in the streamwise and spanwise directions, taking the minimum of the following two criteria; ∆x + = ∆y + ≈ 5 ('+' superscripts indicate wall-units) and ∆x = ∆y ≈ λ min /12, where ∆x and ∆y are the streamwise and spanwise grid spacings and λ min is the smallest wavelength of the rough surface after filtering, defined as the inverse of the filter cut-off wavenumber, k c . A stretched grid is used in the wall-normal direction. In the region of the roughness features, min(h(x, y)) < z < max(h(x, y)), uniform grid spacing is used with ∆z + min < 1 and gradual stretching is applied towards the channel centre with ∆z + max ≤ 5. The grid resolution has been validated in Busse et al. [24] on the basis of a grid refinement study.
The three-dimensional, incompressible Navier-Stokes equations, nondimensionalised by δ and u τ , are discretized using a standard 2 nd order central difference scheme in the spatial domain which operates on a staggered cartesian grid and use a 2 nd order Adams-Bashforth method for temporal discretization. The flow is driven by a constant mean streamwise pressure gradient, which fixes the value of the friction velocity, u 2 τ = − δ ρ dP dx = 1, where ρ is the fluid density and dP/dx is the mean streamwise pressure gradient. An immersed boundary method, described in detail in Busse et al. [24], is used to resolve the rough walls.

Surface samples and topographical properties
Flow over a total of 17 rough surface samples has been simulated. The database includes two carbon-carbon composite surfaces, a concrete surface, a graphite surface, as well as surfaces subject to the processes of casting, hand filing (2 cases), gritblasting, grinding, shotblasting, spark-erosion (5 cases) and replicas of two ship propeller surfaces eroded by periods of service. In order to simplify naming, samples are assigned names as given in Table 1. These are the names used henceforth. The composite and graphite samples were exposed to arc-heating in order to simulate the environment experienced by space vehicles while re-entering the atmosphere. The cast, filed, gritblasted, ground, ship propeller, shotblasted and one out of the five spark-eroded samples (spark-eroded 5 from Table 1) were taken from standard roughness comparators. The remaining four spark-eroded samples were taken from a spark-eroded surface provided by an industrial third-party. These four samples along with the ship propeller samples were selected as different sub-sections from the same respective larger rough surface scan. The concrete sample was taken from a larger block of concrete. Surface plots of all 17 rough surfaces are shown in Appendix B. The composite and filed samples have strong directional alignment of their roughness features. In order to study this phenomenon, two samples of each are evaluated; one having features aligned in the streamwise direction and the other spark-eroded 4 s17 spark-eroded 5 having features aligned in the spanwise direction. In Appendix B, Figures B1 (b) and (c) show the composite samples, s2 and s3 with features aligned in the streamwise and spanwise directions respectively and Figures B1 (e) and (f) show the filed samples, s5 and s6 with features aligned in the spanwise and streamwise directions respectively. The ground sample, s9 also shows strong directional alignment of features in the spanwise direction, as shown in Figure B1 (i). A large number of parameters [27] can be used to characterise rough surfaces. Table 2 displays a broad list of parameters for the current dataset of 17 samples, whose description and computation is given in Appendix A. These parameters are computed for the filtered surface samples. Flack & Schultz [17] mention that skewness, S sk , is a quantitative way of describing whether a rough surface has more pronounced peaks or valleys. A negative value of skewness indicates that the surface is pitted, for example, due to corrosion or surface wear, whereas a positive value indicates roughness due to isolated large peaks, for example, due to deposition of foreign materials (as in biological fouling). A surface skewness value close to zero indicates a more or less homogenous distribution of peaks and valleys. The s1, s2, s7, s14 and s15 samples have a positive value for skewness whereas all other samples have a negative value. Also, s15 has a skewness value close to zero.
The largest correlation lengths are exhibited by s5 and s6 in their respective spanwise and streamwise directions. This is attributed to the strong anisotropy of their roughness features. The directionality of the features of a given rough surface sample can be obtained from the surface texture aspect ratio, S tr which is given by the ratio of the shortest to longest correlation lengths of the sample. If S tr > 0.5 then the sample is regarded as statistically isotropic whereas anisotropic samples have S tr < 0.3 (refer [27]). All samples, with the exception of s2, s3, s5, s6, s9 and s10 have S tr > 0.5 and hence are statistically isotropic. The s10 sample has S tr = 0.41 and can be considered weakly anisotropic. Both the composite samples, s2 and s3 are anisotropic with S tr = 0.28 for s2 and its dominant features oriented in the streamwise direction and S tr = 0.21 for s3 with its dominant features oriented in the spanwise direction.
A parameter called the flow texture ratio, S flow tr , has been defined as the ratio of sample spanwise to streamwise correlation lengths, S flow tr = L cor y /L cor x (refer Appendix A for definitions). This parameter is another indicator of the anisotropy of the roughness features. If S flow tr 1, for example, S flow tr = 29.9664 for the s5 sample, its roughness features have strong directional preference in the spanwise direction (refer Figure B1 (e)) and if S flow tr 1, for example, S flow tr = 0.0345 for the s6 sample, its roughness features have strong directional preference in the streamwise direction (refer Figure B1 (f)).
The effective slope, ES, as introduced by Napoli et al. [28], represents the overall gradient of the roughness elements of an irregular rough surface. Higher values of ES indicate more dense roughness whereas lower values are obtained for relatively sparse roughness. In the case of three-dimensional roughness, the effective slope is computed in the streamwise and spanwise directions and denoted by ES x and ES y respectively. Most samples have similar values of ES x and similar values of ES y . Based on these values, s4, s7 and s8 can be considered relatively more densely rough whereas s9, s11 and s12 can be considered relatively sparsely rough. A closer look at the values of S f /S and ES x from Table 2 shows that 2 × S f /S ≈ ES x for the current dataset. This relation was also pointed out by Napoli et al. [28] in their studies. Figure 4 (a) shows a plot of the two quantities and clearly establishes this relationship as all points fall on the straight line given by 2 × S f /S = ES x . Figure 4 (b) shows the variation of Λ s with ES x and a clear dependence is seen for these two quantities as well. This dependence seems sensible as S f /S is an integral Table 2. Topographical properties of the rough surface samples in this study. Sz,5×5 = δ/6, kcLx = maximum streamwise wavenumber. All length scales are non-dimensionalised by the mean channel half-height, δ. α and αrms are computed in degrees. Refer Appendix A for a description of each property and its calculation. Surface sample naming convention described in Table 1.   part of Λ s . Bons [16] mentions that the mean streamwise forward facing surface angle, α is geometrically related to the Sigal-Danberg parameter. This can be seen from Figure 5 (a), which shows a semilogarithmic plot of the 2 quantities for the samples in this study. From this, it is logical to follow that S f /S is also related to α, Figure 5 (b). This serves as motivation to also look at the variation of the root-mean-square of the streamwise surface angle, α rms with Λ s , Figure 5 (c) and with S f /S, Figure 5 (d). These relationships also confirm that the streamwise forward facing surface angle is approximately linearly related to the frontal area of the roughness elements. α rms was proposed by [2] as an important parameter characterising real roughness in the context of turbine blades. It is also worth noting that tan α ≈ ES x , the tangent of the mean streamwise forward facing surface angle approximately represents the streamwise effective slope.
The above observations establish that ES x , α and α rms are all closely linked to the solidity, S f /S for the current set of samples and as such cannot be regarded as independent parameters. This is particularly important for the parametric fitting studies conducted in Section 5 as, if one of 4 properties enters the fit at a certain stage, none of the other 3 would provide any more useful information at a later stage. Refer Section 5 for further details.

Simulation parameters and results
The streamwise, spanwise and wall-normal domain lengths for all samples along with number of grid cells in each direction and the grid spacing are displayed in Table 3. ∆z + max is the maximum wall-normal grid spacing at the channel centre. The large variation in the streamwise and spanwise domain extents is seen due to the fact that the roughness height for all samples is maintained the same, at k = δ/6 (where k = S z,5×5 ). All simulations are performed at Re τ = 180, for which the roughness Reynolds number, k + = ku τ /ν = 30. Corresponding parameters for the smooth-wall reference case are as follows: L x /δ = 12, L y /δ = 6, L z /δ = 2, n x = 256, n y = 256, n z = 224, ∆x + = 8.4375, ∆y + = 4.2188, ∆z + max = 4.70, mean streamwise velocity, U /u τ = 15.77 and centreline velocity, U c /u τ = 18.44.
Time-averaged mean streamwise velocity profiles in wall-units against the wallnormal distance are shown in Figure 6 (a)-(d) on semilogarithmic axes. Smoothwall profiles have also been shown for comparison. The roughness function, ∆U + is generally measured as the downward shift of the log region for a given roughwall profile from the corresponding smooth-wall profile. However, due to the low Reynolds number in this study, no clearly defined log region is present. Thus ∆U + is computed by subtracting the centreline velocity for a given rough surface simulation from the corresponding smooth-wall centreline velocity [29]. ∆U + values for all samples are given in Table 3. A significant roughness effect is seen for all samples, Table 3. Rough surface sample domain sizes, non-dimensionalised by δ, number of grid cells and grid spacings, for the samples in this study. nx, ny and nz are the number of grid cells in the streamwise, spanwise and wall-normal directions respectively. Also shown are the values of ∆U + and peak profile TKE.  from the downward shift in the mean velocity profiles. There is a wide range, from ∆U + = 1.28 (s6 sample) to ∆U + = 5.02 (s7 sample), despite all samples being scaled to the same roughness height. This is a clear indication that the roughness function depends not only on the roughness height for a given sample but also on its detailed roughness topography. The s6 sample has the smallest roughness function value at ∆U + = 1.28. This is because the roughness features of this sample are strongly aligned in the streamwise direction (refer Figure B1 (f)) and this anisotropic topography gives less resistance to the flow. This leads to a comparatively lower increase in surface friction and hence a smaller value of ∆U + compared to other samples. The s7, s4 and s8 samples show some of the largest values of ∆U + , at 5.02, 4.95 and 4.36 respectively. This closeness in ∆U + values is possibly due to similar values of some of their surface properties; for example, L cor x , S f and Λ s (refer Table 2). The ship-propeller samples, s10 and s11 also exhibit similar values of ∆U + , despite their surface properties showing numerous differences.
The variation of ∆U + with the effective slope is shown in Figure 7, where it is seen that ∆U + increases with ES x . This general dependence of ∆U + on the effective slope indicates that the current set of samples lies in the waviness regime (described by [23]). The curve appears to be levelling out at higher values of ES x implying that the critical ES x separating the waviness and roughness regimes may be close to these values. The variation of ∆U + with Λ s is shown in Figure 8 on a semilogarithmic plot. It can be seen that ∆U + decreases as Λ s increases. The plot also shows a good collapse of the data, with the coefficient of determination,  R 2 = 0.8836 and rms error of the fit, σ = 0.3753. Several parameters such as the roughness density, shape and direction with respect to the mean flow are taken into account while computing Λ s and the data already scales very well with the roughness function. This serves as an initial motivation to perform a more methodical study of surface properties that influence ∆U + , as described in Section 5.
The time-averaged mean streamwise velocity defect profiles for the samples are shown in Figure 9 (a)-(d). U (z) + denotes the temporally-and spatially-averaged streamwise velocity profiles in the wall-normal direction in wall-units. The centreline velocity, U + c , is obtained at z = δ. The profiles are then computed by taking the difference between the two, U + c − U (z) + , for each wall-normal location. They are plotted against the wall-normal distance normalised by the channel half-height, z/δ. The region where these profiles widely differ from each other for the different samples is close to the roughness features. This is seen from the plots for z/δ 0.1. Beyond this region and closer to the channel centre, the profiles follow the smoothwall data to a good degree, thus obtaining a good collapse. This indicates that outer layer similarity is preserved and Townsend's wall similarity hypothesis [30] is satisfied. An important observation from this study is that outer layer similarity for the irregular rough samples is achieved despite the relatively low Re τ and ratio of channel half-height to roughness height of δ/k = 6.
Profiles of the turbulent kinetic energy (TKE) in wall-units, given by TKE = [(u/u τ ) 2 + (v/u τ ) 2 + (w/u τ ) 2 ]/2 are shown in Figure 10 (a)-(d), plotted against the wall-normal distance, along with corresponding smooth-wall profiles. Additionally, Figure 11 shows profiles of streamwise, (u/u τ ) 2 , spanwise, (v/u τ ) 2 and wall-normal, (w/u τ ) 2 fluctuations for the smooth-wall case and the s8 sample as an example. This figure enables us to see graphically the contribution of each to the TKE. It is clear from the figure that peak streamwise fluctuations are higher in magnitude than both the peak spanwise and peak wall-normal fluctuations for the smooth-wall as well as the s8 sample. Also, the smooth-wall peak streamwise fluctuations are higher than the s8 peak streamwise fluctuations. The above two observations were found to be true for all rough surface samples in this study. To quantify the anisotropy of the velocity fluctuations and hence TKE, the diagonal components of the normalised Reynolds stress anisotropy tensor, where δ i,j is the Kronecker delta, are computed at the peak value of TKE for all samples. Since it is known that In general, an increased amount of roughness is accompanied by a decrease in the peak streamwise fluctuations [29,31,32], which means samples with higher peak magnitudes have lower ∆U + compared to samples with lower peak magnitudes.
This influences TKE as well, as samples with higher ∆U + show a lower value of peak TKE (refer Table 3). Within the region of the roughness features and close to z/δ ≈ 0, all samples show TKE values greater than the smooth-wall value. The reason is that rough-wall velocity fluctuations can occur very close to the roughness features, including at and below the mean wall location, z/δ = 0. This is not possible in the case of smooth walls. In general, an increase in TKE is observed close to the rough walls whereas a collapse with the smooth-wall data is seen away from the walls. For the smooth-wall case, the TKE peak is located at z/δ ≈ 0.1, whereas rough-wall peaks for all samples are located slightly above that.

Parametrisation of topographical properties for ∆U + and TKE
Although a general solution to the roughness problem must be delayed until a more complete dataset is available for fully rough cases and including a wider Reynolds number range, it seems sensible to try to make as much progress as possible with the present set of restricted data, first of all to find out the issues that arise in formulating a general empirical model and, secondly, to guide the next set of simulations to best exploit the available computational resources. As can be seen from Table 3, there is a large variation in the roughness function ∆U + that must be due to other parameters, besides the height, that govern the surface topography. To obtain surface properties that possibly influence the roughness function, a fitting process is employed whereby ∆U + is plotted against a combination of surface properties and the quality of the fit is improved by successively adding other properties. Additions are made based on a systematic testing of all available properties using specific mathematical forms (algebraic, exponential, logarithmic or power) and selecting the property and form that gives the best possible fit. In Table 4, all forms tried for this process are listed for an example surface property, p. The particular form for a property may not necessarily be optimal, since only the above-mentioned four mathematical forms are tested and there may be other forms which might influence the fitting process. Combinations of surface parameters are denoted by λ n , where n = 0 for a baseline model, n = 1 for a 1-parameter model and so on. An example for n = 2 could be where c 1 and c 2 are fitting constants.
To measure the success of the method we use the root-mean-square error, σ between the data and a straight-line curve fit using the derived parameter, as well as the value of the coefficient of determination, R 2 of the fit. In order to maximise the number of property combinations tested, fits giving the 3 lowest values of σ (which means the 3 highest values of R 2 ) are retained for further improvement. This means, for example, in case of n = 0 (or the baseline) model, where ∆U + is fitted with a single surface property, fits of those properties which obtain the 3 lowest values of σ are selected for further improvement by addition of more properties. However, in the following description, only the best fits are reported. The reason for selecting multiple property combinations is also because a given property combination that gives the lowest value of σ for n = 1, for example, may not necessarily give the lowest σ value for n = 2 because of the interactions between different surface properties. Parameters are continued to be added until no significant improvement of the fit is obtained and the fit with the final lowest value of σ is selected as the best.
For n = 0 to fit ∆U + , we consider the performance of a solidity parameter, expressed here in logarithmic form, Table 4. Mathematical forms of properties tested during the fitting process. c is a fitting constant and p is the value of a given surface property.  Figure 12. Linear fits to the DNS data, correlating the roughness function, ∆U + with different parameters λ 0 , λ 1 ,λ 2 and λ 3 , corresponding to equations (9), (11), (12), (13) in the text.
Using just this parameter, the best fit to the data is ∆U + = aλ 0 +b, with a = 2.0438 and b = 8.9035 and with σ = 0.3807 and R 2 = 0.8802. The resulting straight line is plotted in Figure 12 (a) and already a reasonable fit to the data can be seen. Given the success of the simplest measure, our strategy is to introduce modifications to the definition of λ 0 based on additional surface properties as shown in Table 2.
Extensions to the solidity parameter have already appeared in the literature and, for irregular surfaces, van Rij et al. [15] redefined the parameter introduced by Sigal & Danberg [14], previously given in equation (4). A generalisation of this approach is to set where S w is the wetted area of the forward-facing elements of the surface and Sigal & Danberg used β = 1.6 (note that Sigal & Danberg used the inverse of this parameter whereas we prefer a definition where λ can be interpreted as the solidity or density of the roughness). However, using this value of β in the present study led to no improvement in the standard error. A separate exercise was undertaken to optimise the value of the exponent, giving β = 0.18, but with a barely measurable increase in R 2 . The reasons for the failure of this additional term are clear from Table 2, since for the types of roughnesses considered the wetted area parameter is always about half of the planform area i.e. there is an approximate symmetry (of forward-facing and rearward-facing roughness elements) in the roughness samples. Thus the term S f /S w introduces no additional useful information. One could continue using (10) as the reference parameter and get the same results, but we prefer only to use parameters that are justified by the data and instead revert to the simple solidity, λ 0 as our baseline property. The next step is to test each of the potential surface parameters as modifications to λ 0 . The best two, with almost identical performance, were the streamwise correlation length parameter L cor x /S z,5×5 and the flow texture ratio S flow tr . The success of both suggests that the spanwise correlation length is less important and we retain the best-performing parameter, with a single optimised coefficient for n = 1 to give The improved fit to the data is shown in Figure 12 (b), with σ = 0.3073 and R 2 = 0.9220. It is interesting that a streamwise correlation length enters as the nextmost important parameter after the solidity since this type of parameter doesn't appear in many correlations. The parameter is additionally intriguing since dense roughness cases will have low values of L cor x /S z,5×5 and from the correlation this would lead to lower λ 1 and hence lower ∆U + , which is indeed what is observed (Jiménez [9]). The absence of dense roughness cases (S f /S > 0.15) from the current sample set means that we cannot test this fully, and addressing this point would be a priority for future simulations.
We continue the process to define the best models for n = 2 and 3. The best model for n = 2 is found to include the relative rms roughness height parameter, S q as with σ = 0.1806 and R 2 = 0.9731. The best model for n = 3 includes the skewness, S sk as with σ = 0.1383 and R 2 = 0.9842. Figure 12 (c) and (d) show the continued improvements seen with the λ 2 and λ 3 representations. The largest remaining errors in the fit to the data are less than 0.1u τ . Additional parameters were tested, but with no significant further improvements found. Fit parameters for λ 0 , λ 1 , λ 2 and λ 3 are summarised in Table 5. Tests were also run by removing parameters individually, confirming that a ranking in order of importance is (i) solidity, (ii) streamwise correlation length non-dimensionalised by the mean peak-to-valley height, (iii) rms roughness height non-dimensionalised by the mean peak-to-valley height and (iv) skewness. Note that the roughness height is not one of these parameters since all the cases were run for the same S z,5×5 . Had the simulations been in the fullyrough regime, the equivalent sand-grain roughness k + s,eq would be determined as a constant (dependent on all the above parameters) multiplied by some suitable measure of the roughness height e.g. S + q or S + z,5×5 . Both rms roughness height and skewness are part of the Flack & Schultz model [17] so it is no surprise to see them here. Also, as we have seen, the effective slope, mean and rms streamwise forward facing surface angles are proportional to the solidity for these surfaces, so cannot Table 5. Best fit parameters for λ0, λ1, λ2 and λ3, corresponding to equations (9), (11), (12), (13) in the text. ∆U + = aλn + b, where a = slope of the fit and b = y-axis intercept. σ = rms error of the fit and R 2 = coefficient of determination. be considered as independent parameters. We should caution that the above analysis is only the first step. As more samples are added covering different types of roughnesses (dense, for example) we might expect that additional parameters would be required. What is important is that we now have a systematic method to incorporate additional parameters. We caution again that the models in equations (9) to (13) should not be used for k + s,eq since the current data were all taken in the transitionally rough regime. What we have been able to do is identify a number of parameters that contribute significantly to the roughness function in this regime and it is likely that the same parameters contribute to the determination of k + s,eq . The same numerical coefficients would only be found if all the samples followed the same path through the transitionally rough regime, which is unlikely.
A similar approach as above was also utilised to fit surface property data to the value of peak TKE from Table 3. Different parameters are seen to appear in the model as the fluctuations behave differently with property combinations as compared to ∆U + . The various models obtained for this process are given below. n = 0: n = 1: n = 2: n = 3: The fits for the above models are shown in Figure 13 (a)-(d) and fit parameters are summarised in Table 6. Although fits are reported only up to n = 3, further improvements, following the same systematic fitting approach as described for ∆U + , are seen up to n = 5. Values up to σ = 0.0244 and R 2 = 0.9820 are obtained when the average roughness height, S a , in its algebraic form and shortest correlation length, S al , in its power form (both properties non-dimensionalised by S z,5×5 ) are included in the model. However, due to the relatively small size of the sample database for fitting purposes, the influence of these properties is probably not as significant as the ones seen up to n = 3. For n < 3, property combinations other  Figure 13. Linear fits to the DNS data, correlating the peak TKE, with different parameters λ 0 , λ 1 ,λ 2 and λ 3 corresponding to equations (14), (15), (16) and (17) in the text.  (15) and (16) finally lead to equation (17), which ultimately gives the lowest σ value of all final fits tested, it is selected as the best fit and is discussed here. The baseline parameter is the mean forward-facing surface angle, α, which is an angle parameter as opposed to S f /S, which is an area parameter, seen in the case of ∆U + . However, it has been shown in Figure 5(b) that S f /S and α are approximately linearly related. Also, it is understood that higher values of α correspond to a higher roughness effect (from higher values of ∆U + , refer Tables 2 and 3) and hence its influence on the fluctuations would be significant. Other baseline parameters that gave σ values comparable to α include α rms , S f /S, and ES x , all in logarithmic form. This is not too surprising as the four parameters are interrelated. Surface skewness, S sk is the next important parameter and after that comes the rms roughness height non-dimensionalised by the mean peak-to-valley-height, S q /S z,5×5 , both of which also appear in λ 1 and λ 2 respectively for the other baseline parameters, albeit in different forms. Spanwise effective slope, ES y is the next parameter to enter the fit, the appearance of which could relate to how the streamwise flow navigates around roughness features. It would be interesting in the future to understand why certain parameters enter the fit as opposed to others, which is not considered in this study. The baseline TKE fit is not as good as the fit for ∆U + but significant improvement is seen until λ 3 . A separate fitting study conducted for the profile peak streamwise fluctuations, (u/u τ ) 2 alone also gave the same properties influencing the fit, although in a slightly different order.

Conclusions
A direct numerical simulation study of 17 industrially relevant rough surface samples with the same physical roughness height and at the same friction Reynolds number, in the transitionally rough regime has been presented. Mean streamwise velocity profiles show a clear roughness effect for all samples. A wide range of the roughness function is obtained, from ∆U + = 1.28 to 5.02, despite all samples being scaled to the same roughness height of k = S z,5×5 = δ/6. It is thus clear that the roughness effect depends not only on the roughness height but also on the detailed roughness topography.
A process to determine which surface properties influence ∆U + is then formulated. The process involves fitting the roughness function with one or more surface properties. Addition of properties to the model is based on a systematic testing of all available properties and selecting the combination that provides the best fit as measured by the value of root-mean-square error between the data and the fit as well as the coefficient of determination. A similar procedure is also applied to fit the profile peak turbulent kinetic energy. Optimized fits are developed for ∆U + and peak TKE. Properties influencing ∆U + include the solidity (S f /S) and surface skewness (S sk ), which are known from literature, and additionally, the streamwise correlation length (L cor x ) and the rms roughness height (S q ). Properties influencing peak TKE are slightly different. The surface skewness and rms roughness height are still seen in the TKE fits together with the mean forward-facing surface angle (α) and spanwise effective slope (ES y ). The identification of key surface parameters that determine critical flow properties is an important conclusion of the current work. The final fits obtained for both ∆U + and peak TKE are of high quality (with the value of R 2 = 0.9842 for ∆U + and 0.9288 for TKE). Since all the simulation data was taken in the transitionally rough regime, this process is used only to determine the surface properties that influence ∆U + and not predict an equivalent sand-grain roughness, which would require data in the fully-rough regime. The process establishes that possibly the same properties would influence flow parameters in the fully-rough regime as well. An extension of the procedure to the fully-rough regime is straight forward but requires larger computational resources. Increasing the size of the data set by introducing more surfaces having different properties (for example, surfaces in the roughness regime with respect to their effective slope) would serve to increase the generality of the fitting process and could be considered a priority for future work.

Acknowledgements
This work has been supported by the Engineering and Physical Sciences Research Council grant number EP/I032576/1. We would also like to thank the A large number of surface parameters are used to characterise the rough surfaces used in the current study. [27] gives a very extensive range of metrological parameters that may be used to describe rough surfaces in general. It is important to note that since the mean surface height is taken as the mean reference plane, z = 0, we get a relation for the mean surface height, h as The shortest correlation length is defined as S al = min (l∆s) 2 + (m∆s) 2 |R h (l, m) ≤ 0.2 and the longest correlation length is defined as S sl = max (l∆s) 2 + (m∆s) 2 |R h (l, m) ≥ 0.2 ∩ (l, m) ∈ central lobe .
The central lobe of the areal autocorrelation function is the simply connected area where R h > 0.2 that contains (0, 0). The surface texture aspect ratio, S tr is given by the ratio of the shortest to longest correlation lengths, The surface correlation lengths are given as, streamwise correlation length: L cor x = min{l∆s|R h (l, 0) ≤ 0.2}, spanwise correlation length: L cor y = min{m∆s|R h (0, m) ≤ 0.2}.
A parameter called the flow texture ratio, S flow tr , which depends on the streamwise and spanwise correlation lengths, has been defined as

A.3. Aerodynamic parameters
In the context of aerodynamics, several other geometric parameters for the characterisation of rough surfaces have been defined. In the context of two-dimensional roughness, Napoli et al. [28] introduced the streamwise effective slope, ES. The solidity or frontal area ratio has been used extensively in literature and is an indication of the roughness density. It is given by the ratio of the total frontal area of all roughness elements to the planform area of the sample, S f /S. The generalised Sigal-Danberg parameter as defined by van Rij et al. [15] is given as where S = L x L y is the planform area of the corresponding smooth surface. S f is the total frontal area of all roughness elements and is given as Ly 0 r x ∂h ∂x + r y ∂h ∂y W (x, y)dxdy.  Table 1 for naming convention. All plots have the same colourbar, shown at the bottom. Plots coloured by roughness height. Refer Table 1 for naming convention. All plots have the same colourbar, shown at the bottom.