Characterization of Flow Structures Induced by Highly Rough Surface Using Particle Image Velocimetry, Proper Orthogonal Decomposition and Velocity Correlations

L. R. Andersson¹, I. A. S. Larsson¹, J. G. I. Hellström¹, P. Andreasson^1,2, A. G. Andersson¹, T. S. Lundström¹

¹Division of Fluid and Experimental Mechanics, Luleå University of Technology, Luleå, Sweden.
²Vattenfall Research and Development, Älvkarleby, Sweden.

DOI: 10.4236/eng.2018.107028                

Abstract

High Reynolds number flow inside a channel of rectangular cross section is examined using Particle Image Velocimetry. One wall of the channel has been replaced with a surface of a roughness representative to that of real hydropower tunnels, i.e. a random terrain with roughness dimensions typically in the range of ≈10% - 20% of the channels hydraulic radius. The rest of the channel walls can be considered smooth. The rough surface was captured from an existing blasted rock tunnel using high resolution laser scanning and scaled to 1:10. For quantification of the size of the largest flow structures, integral length scales are derived from the auto-correlation functions of the temporally averaged velocity. Additionally, Proper Orthogonal Decomposition (POD) and higher-order statistics are applied to the instantaneous snapshots of the velocity fluctuations. The results show a high spatial heterogeneity of the velocity and other flow characteristics in vicinity of the rough surface, putting outer similarity treatment into jeopardy. Roughness effects are not confined to the vicinity of the rough surface but can be seen in the outer flow throughout the channel, indicating a different behavior than postulated by Townsend’s similarity hypothesis. The effects on the flow structures vary depending on the shape and size of the roughness elements leading to a high spatial dependence of the flow above the rough surface. Hence, any spatial averaging, e.g. assuming a characteristic sand grain roughness factor, for determining local flow parameters becomes less applicable in this case.

Keywords

CFD, Validation, Hydraulic Roughness, PIV, Hydropower

1. Introduction

Water tunnels are frequently used to convey water to and from hydropower turbines and in other sectors of infrastructure. The tunnels are often a key part of the design and their durability is vital for the continued operation. Tunnel excavation by rock-blasting is a relatively swift, and therefore popular, method compared to using tunnel boring machines [1] . One significant drawback of the method would be that the resulting walls of the tunnels have varying cross section and considerable surface roughness [2] . The scales of roughness in unlined hydropower tunnels range from a few millimeters to meters, which may be in the range of ≈10% - 20% of the hydraulic radius. A likely treatment of rough walls by today’s industrial standards is to estimate the roughness from the actual physical features of the roughness, e.g. estimates according to Manning, Chezy, grain size distribution [3] . These approaches have a long track record but only provide spatially averaged quantities and nothing about the actual dynamics within the tunnel. The physics of the flow in highly rough tunnels includes large variations and gradients of pressure [4] and velocity [5] , resulting in intermittent pressure forces and increased local shearing load acting on the walls of the tunnel. Such forces may very well jeopardize the structural integrity of the walls, causing events including erosion or even partial collapse of the tunnel [6] . Collapsing hydropower tunnels is a known problem and documented cases show that even after 30 - 40 years of usage some tunnels have experienced sudden severe failure [7] . The connection between the details of the flow and fatal failure of tunnels is not completely understood, but it has been theorized that intermittent pressure fluctuations directly coupled to surface roughness may have induced or facilitated the process. Similar results were presented in a recent study by [8] , showing that pressure-driven cyclic injection of fluid into rock material leads to larger damage, and at lower pressure than by static injection. Hydropower tunnels are a very hostile environment to measure within, hence, accurate results are almost non-existing. Capturing the geometry of existing tunnels also presents several problems: The tunnels themselves are dark and humid, making accurate measurements difficult [1] . Closing down the tunnels, and thereby any machinery operating downstream, is very expensive and puts the system in danger [6] . Therefore, there exist very few cases where experiments have been performed on models reflecting actual tunnels. This problem was highlighted by [9] but to a degree still remain today.

It is well established that rough walls modify the behavior of the flow [10] , [11] however to what extent has been thoroughly debated. For rough surfaces, the turbulence is associated with shear layers formed at the crests of the roughness elements where flow separation may occur [12] . For surfaces of sufficiently large roughness, individual surface aberrations frequently penetrate into the inertial sublayer [13] [14] leading to a breakdown of the logarithmic law of the wall. Yet to this day, a common way to model flow in hydropower tunnels in industry is to replace the natural roughness with numerical wall-roughness functions. This method relies on the conventional concept that roughness effects are confined to the inertial sublayer near the surface and have no direct effect on the outer flow [15] , a theory which for some flow cases have been questioned [16] [17] . Many studies have been directed at the understanding of flow heterogeneity over rough surfaces. Measurements by [18] and [19] provided proof of the localized velocity perturbations connected to the rough surface. This study generally has significantly higher Re and larger roughness elements (relative to the hydraulic radius) as compared to the mentioned studies, however, this is an applied study and the results serve well for comparison. A study by [4] employed flush mounted pressure sensors at points of interest, such as peaks and valleys, at the rough surface. Sensors placed on peaks revealed elevated frequencies of the pressure fluctuations as compared to those placed in valleys and there were a large variation in the pressure magnitude as a function of the spatial coordinates. The present study, which is carried out in the same channel as in [4] , will complement the previous study by further analysis using PIV. The scope of this article will be the following: 1) To further visualize the events surrounding singular roughness elements in order to assess the current industrial evaluation standards of uniform roughness treatment. 2) To visualize the spatial heterogeneity connected to the rough walls of the tunnel. 3) To bridge the gap between flow in hydropower tunnels and other fluvial flows. Due to the large-scale and randomness of the surface roughness, local flow patterns are unpredictable both spatially and temporally. Applying only temporal averaging for the analysis of the flow may, therefore, provide misleading results. Adding a spatial averaging to a plane parallel to the mean flow may filter away the smallest perturbations due to e.g. laser reflections or insufficient boundary layer resolution, while accounting for the largest events such as flow separation. This technique is called double averaging [20] . To identify the flow structures created from the rough surface interaction, the PIV-data is analyzed using proper orthogonal decomposition (POD).

2. Experimental Setup

The experimental setup consisted of a closed loop water system with a 10 m long rectangular Plexiglass (PMMA) channel having one rough surface, a pump, an electromagnetic flow meter, two tanks placed on different levels and a PIV-system. The function of the tank placed upstream of the channel is to provide a constant head on the system and to avoid air entrainment inside the channel. A schematic of the experimental setup can be seen in Figure 1 (not to scale).

A detailed description of the setup can be found in [4] . The flow rate was approximately 62 l/s which varied with approximately 4% throughout the campaign. This corresponds to a Reynolds number of $\mathrm{Re}=y{U}_0/2\nu \approx 200,000$ , where y/2 is the channel half-height, $U_0$ is the bulk-velocity and $\nu$ is the viscosity. The channel has a depth of 0.145 m and a width of 0.25 m. The measuring section, represented by the red box in Figure 1, is positioned about 6.8 m downstream of the tunnel inlet and has a length of 0.48 m. Figure 2 depicts the

Figure 1. The setup used in the experimental campaign. The channel along with the laser setup have been mirrored in this figure to provide a more apprehensible overview of the setup.

Figure 2. The rough surface in the measuring section, (a), the two coloured lines mark the position of the two measured planes also depicted in (b). The colour of the surface represents the relative slope of the surface topography. The flow is from left to right in the figure.

section of the rough surface over which the flow was measured during the experiments. One measuring plane was placed in the center of the channel (denoted middle), while another was placed closer to the camera (denoted upper). The color of the surface represents the slope of the surface topography and the two colored lines mark the position of the two measured planes. As can be seen in the figure, a ridge is passing through both measuring planes at x ≈ 7.06 m. The ridge is of interest for the measurements since the maximum height relative to y = 0 (the mean height) is similar for the two lines, additionally, the gradient of the surface is also nearly the same. However, the relative size of the ridge differs significantly between the lines. For the blue line defining the middle plane, the final slope of the ridge is very sharp but the roughness leading up to the element is relatively small, consequently, the relative size of the roughness element is small. On the red line (upper plane), the ridge is preceded by a “valley”, making the relative height larger. Ideally, four rough walls would have been used to keep the setup as realistic as possible. However, PIV requires optical access from at least 2 directions and therefore only one rough wall was used.

In this study, a right-handed coordinate system is employed with the x-coordinate (u-velocity component) originating from the tunnel entrance pointing in the flow direction. The y-coordinate (v-velocity component) is perpendicular to the rough surface with y = 0 defined as the average elevation of the rough surface, hence all presented heights are relative to the average height of the surface. Accordingly, the z-coordinate originates from the bottom wall in the flow direction.

2.1. The Rough Surface Model

The rough surface model used in the experiments is a 1:10 scale side wall of an existing rock tunnel whose topography has been captured by high resolution laser scanning, a method which has been proven efficient for determining surface roughness [1] . One of the main characteristics used to describe a rough surface is the height distribution function $p\left(k\right)$ , in this case the Gaussian distribution. The meaning of $p\left(k\right)$ is that the probability of any surface height between k and $k+dk$ is $p\left(k\right)dk$ [21] . Another important factor for characterizing a rough surface is the root mean square (RMS) roughness factor, which describes the average elevation of the roughness elements on the surface. In the current study the RMS roughness factor is denoted as the equivalent sand grain roughness factor $k_s$ for the surface and is defined as [1] [22]

$k_s^2={\displaystyle {\int }_{-\infty }^{+\infty }{h}^2}p\left(h\right)dh.$ (1)

As mentioned, $k_s$ is solely based on the height on the rough surface and does not take into consideration e.g. shape or aspect ratio of the roughness. To evaluate the spatial difference, the auto-correlation function R(r) over a specified length L is introduced in the stream wise direction (x-direction) of the rough surface [21] . The result can be seen in Figure 3, which is a z-direction average of the autocorrelation function.

Integrating the auto-correlation function according to Equation (2) produces the integral length scale of the surface

${\tau }_r\equiv {\displaystyle {\int }_0^{\infty }R\left(x\right)dx}.$ (2)

Conclusively, $k_s$ is a quantity representative for the roughness height while

Figure 3. The auto-correlation function of the rough surface.

${\tau }_r$ represent the roughness length of the surface. Using Equation 1, $k_s$ is determined to 9.4 mm, while ${\tau }_r$ is found to be 39.1 mm using Equation 2. Hence, $k_s$ is about 6.4% of the hydraulic radius, which can be compared to the largest global roughness elements being about 20% of the hydraulic radius. The difference is a clear indicator of the spatial heterogeneity of the rough surface used in this study. For additional numerical comparison, the relative height of the ridge in the middle plane is about 9.44 mm, which is very close to $k_s$ . One can hereby conclude that the ridge studied is a fitting representation of the roughness of the entire surface. Additionally, the ridge in question is far from unique on the surface but appears at regular intervals, therefore, the sample size is deemed large enough to be spatially independent.

2.2. PIV-Setup and Error Estimation

The PIV-system used is a commercially available system from LaVision GmbH which has been applied in a number of studies, including [23] . It consists of a Litron Nano L PIV laser, i.e. a double pulsed Nd: YAG with a maximum repetition rate of 100 Hz and a pulse energy of 50 mJ. A 10-bit LaVision Flow Master Imager Pro CCD-camera with a spatial resolution of 1280 × 1024 pixels per frame is used for image acquisition. Sheet optics and mirrors produced a 1.5 mm thick laser sheet and directed it to the desired location, and a Nikon 50 mm f/1.8 D lens was fitted on the camera. The laser was mounted on a traverse allowing a simultaneous repositioning of the laser sheet and camera of up to 500 mm in the x-, y- and z-directions. The traverse was placed and operated independently of the experimental setup, a precaution preventing any large loads acting on the channel and to prevent vibrations from the rig to interfere with the camera. To cover the entire measuring section the traverse needed to be repositioned between image capturing. Accordingly, the planes had to be divided into 14 (middle) and 24 (upper) smaller subsets which were measured individually and then manually merged together. The spatial dimensions of each subset are about 100 mm (x-direction) by 80 mm (y-direction). To consider the laser sheet attenuation in the image periphery, subsequent positions are set to give a 30 mm overlap of the images. Stitching together the domain using the measured sets may create discontinuities in the domain, due to the different sets not matching perfectly. However by careful stitching, a good statistical convergence and choosing appropriate overlap the largest discrepancy between two sets which was no more than 1.3% in the streamwise direction. In the spanwise direction, the discrepancy is less than 1%. The tracer particles used were the previously proven feasible [24] Akzo Nobels Expancel 461 WU 20 hollow thermoplastic spheres with a diameter ranging from 2 μm to 30 μm, and a density of 1.2 g/cm³. The measurements were performed with a frequency of 75 Hz during 9.49 s, corresponding to a total of 712 image pairs for each recorded set. To account for the localized variations in velocity and pixel displacement, the time interval between the laser pulses ranged from 150 μs - 275 μs depending on the measuring position. The results were a typical mean displacement over the whole velocity field of 0.3 pixels in the y-direction and 7 pixels in the x-direction with a characteristic particle image diameter of 2 pixels. At approximately 490 samples the temporally averaged velocity converges towards a stable value. Beyond 560 samples the velocity differ no more than 1.5% from the converged value. A PIV experimental setup consists of several sub systems, and hence there are a number of potential error sources. The overall measurement accuracy in PIV is a combination of a variety of aspects extending from the recording process all the way to the methods of evaluation [25] . A cornerstone in all experimental design is to randomize the measuring procedure. By proper randomization, the effects of extraneous factors that may be present have less impact on the result [26] . The measurement uncertainties consist of those due to systematic biased errors and random precision errors (or due to erroneous measurements) [27] . The biased error associated with the scaling from pixels to meters is estimated to be 0.5% as derived from measurements over a known length scale. The primary source of random error is introduced by the sub-pixel estimator in the cross-correlation. This error is estimated to be 10% of the particle image diameter, which is the diameter in pixels of the particle as seen through the camera [28] . The mean particle image diameter in the present case is about 2 pixels, and a typical displacement between image pairs is 7 pixels in the main flow direction. Therefore the estimated random error of the measured velocity vector in each interrogation area is about 4% for the streamwise velocity component. The rough surface reflected light from the laser which in some images saturated the camera, inhibiting measurements of the near-wall flow. However, these effects where highly localized and within $0\le y/k_s<1$ of the rough surface, and hence, never affected the bulk-flow. Since the rough surface was placed on a side wall of the channel, the camera was placed above the channel facing downward, see Figure 1. The roughness elements closer to the camera sometimes covered parts of the plane intended for measuring, generating additional difficulties in measuring the near wall behavior of the flow. The laser sheet was initially positioned at the center of the channel (middle plane) to get a measurement where the effect of the side walls was as small as possible, and to get a measurement over a distinct roughness element. The second measurement section (upper) was placed at the same x-coordinate as the first one but closer to the camera, see Figure 2.

2.3. PIV Post-Processing

Post-processing of the PIV-data was done using the commercial software DaVis by LaVision [29] . To calculate the particle displacement, a min/max filter for particle intensity normalization followed by a multi-pass scheme with decreasing window size and offset was used. The interrogation window size was 32 × 32 pixels for the first pass and 16 × 16 pixels for the second pass with adaptive window shift, both with an overlap of 25%. The cross-correlation was performed using the standard cyclic FFT-algorithm with a three-point Gaussian peak fit to estimate the sub-pixel displacement, followed by vector post-processing by applying a median filter to reject spurious vectors (less than 2%) and to interpolate from surrounding interrogation windows [30] . The processed data were imported into Matlab using PIV mat where further analysis was performed. The double averaging process is performed through two decompositions, the first part is the Reynolds decomposition where $\theta =\bar{\theta }+{\theta }^{\prime }$ . $\bar{\theta }$ is the temporally averaged quantity and ${\theta }^{\prime }$ is the quantity fluctuating in time. The second part is the spatial decomposition, $\bar{\theta }=\langle \bar{\theta }\rangle +\tilde{\theta }$ where the angle brackets denote spatial averaging over the desired plane and tilde denotes the temporal deviation from the double averaged component $\langle \bar{\theta }\rangle$ . The Reynolds decomposition is applied when post-processing the raw PIV-images, while the spatial decomposition is applied a posteriori on the processed PIV images. POD is an algorithm which determines and hierarchically ranks the dominant structures in the flow with respect to their energy content, allowing the statistical capture of flow structures despite eventual shortcomings such as insufficient temporal resolution and/or sample size. The POD modes stem from calculation of the singular eigenvalue decomposition

$K{v}_i={\lambda }_i{v}_i$ (3)

of the auto covariance matrix K, given by

$K=U^{T}U$ . (4)

i is the total number of eigenvalues, U is a matrix where the columns consist of the instantaneous fluctuating velocity snapshots, according to

$U=\left[\begin{array}{cccc}|& |& ⋮& |\\ {u^{\prime }}_0& {u^{\prime }}_1& ⋮& {u^{\prime }}_{Nt}\\ |& |& ⋮& |\end{array}\right]$ , (5)

where ${u^{\prime }}_j=u_j-\bar{u}$ for $j=\{1,\mathrm{2...}{N}_t\}$ . $N_t$ is the number of snapshots, 712 in the current case. The modes are then calculated and normalized by

${\varphi }_i=\frac{{\displaystyle {\sum }_{n=1}^{N_t}{v}_n^i{u}_n}}{\Vert {\displaystyle {\sum }_{n=1}^{N_t}{v}_n^i{u}_n}\Vert }$ . (6)

The method stems from [31] and a more recent introduction to POD can be found in [32] .

3. Results and Discussion

The middle plane is placed at z = 125 mm, and is represented by a diamond symbol in the figures. The upper plane was placed at z = 165 mm and will be represented by an x symbol in the figures. To avoid cluttering only a portion of the data have been plotted, typically every fourth point. This does not affect the results and is solely for the purpose of making the data easier to distinguish. The velocity components of the flow (u, v) are denoted as the vector u. To evaluate the flow the u- and v-components of the velocity were averaged over time for one measurement (see Sec. 2.) to produce the temporally averaged velocity components $\bar{u}$ . Some of the results are then spatially averaged in the streamwise direction, denoted by $\langle \bar{u}\rangle$ . In the first section below, temporally averaged and Quadrant analysis of the velocity to discern the spatial heterogeneity of the flow are presented. In the second section, integral length scales applied to the temporally averaged velocity are discussed and in the third section POD is applied on the instantaneous velocity field. The instantaneous contribution to the Reynolds stresses are calculated according to

$\bar{u}{\bar{v}}_Q=\frac{1}{T}{\displaystyle {\sum }_{i=1}^{n}S}\left(u_i-\bar{u}\right)\left(v_i-\bar{v}\right)$ , (7)

where S is a sorting term. If $\bar{u}\bar{v}$ falls into quadrant Q then S = 1, otherwise S = 0. T is the total measuring time for each sample. An introduction to Quadrant analysis can be found in [33] or [11] .

3.1. Average Velocity and Quadrant Analysis

Figure 3 shows the u-component of the time averaged velocity field of the middle measurement plane. The flow field above the rough surface exhibits a highly localized behavior induced by roughness elements, exemplified by the zone of high velocity formed at the crest of the roughness element (the ridge) positioned at x $\simeq$ 7.06 m.

The average bulk velocity in the channel is U₀ = 1.562 m/s. A zone of negative velocity (recirculation zone) is present behind the crest of the ridge, an effect similarly visualized by [18] . Consistent with the previously mentioned study, the negative horizontal flow in the separation zone is about $-0.15U_0$ .

The asymmetric channel flow case (one rough wall opposite of a smooth one) has been well documented by [34] . Traditionally, the rough surface acts as a sink for momentum for the flow, due to the outer similarity treatment the velocity close to the surface can then be approximated using the logarithmic law of the wall. Accordingly, the maximum of the double averaged streamwise velocity component is shifted away from the rough surface [35] .

For the current case, the maximum velocity ${\bar{u}}_{\max }=1.45U_0$ is shifted towards the rough surface ( $y/k_s=4.1$ ) (see Figure 5), similar to the results of [24] . The localized roughness aberrations produce flow alterations of significant magnitude, which becomes representative for the flow close to the rough surface when applying spatial averaging. Similar local velocity alterations, for flow over a gravel bed, were reported by [19] . From Figure 5, a note can also be dedicated to the spatial heterogeneity of the flow. The measured planes are positioned with a lateral distance of only 40 mm between them yet the position of the maximum velocity differs substantially, as the maximum of the double averaged velocity is positioned closer to the center of the channel for the upper plane ( $y/k_s\simeq 7$ ). The size of the roughness elements differs no more than 3% between the planes but the standard deviation of the roughness elements is significantly higher in the upper plane. Quadrant analysis is a method to disclose the instantaneous point-contribution of the velocity fluctuations to the Reynolds stress in relation to a defined quantity H, defined as

$H=\frac{u^{\prime }{v}^{\prime }}{\overline{u^{\prime }}\overline{v^{\prime }}}$ . (8)

In short, the contribution to the Reynolds stresses is divided into one of four possible quadrants depending on the sign of the instantaneous velocity fluctuations. As H increases, low magnitudes of the instantaneous velocity products are sorted out, thus only the significant contributions are left for comparison. Quadrant 2 events generally, but not always, represent ejection and similarly, quadrant 4 events represent sweeps [11] . Figure 6 present the quadrant analysis applied in at the points presented in Figure 4. Point a) is placed at ${\bar{u}}_{\max }$ showing an overwhelming dominance of Q2 events for all H, consistent with [18] , indicating ejections of low velocity fluid away from the rough surface. Point b) and c) are placed relatively close to each other, both on either side of the shear layer separating the accelerating and recirculating flow. Both points show the highest measured $\bar{u}\bar{v}$ magnitudes of all points, 8.41 and 7.67 m²/s² respectively. Point b) show a similar magnitude for Q2 and Q4 for $H=0$ , however, the most

(a) (b)

Figure 4. The left figure (a) denote the $\bar{u}$ -component of the flow field from the middle measuring plane, the colour scale represents velocity and is given in m/s, the cyan box denotes the subset in which the POD was performed. The right figure (b) clarify the points in which the quadrant analysis was performed.

Figure 5. The double averaged u-component of the velocity for both measured planes, is the average free stream velocity. The red line denotes the height of the ridge, below the line the available amount of spatial samples decline.

dominant event for $H>2$ is Q1. Point c), which is placed in the separation zone, show an overall dominance of Q4 events for all H. This is an indicator of sweeps of high velocity fluid moving towards the wall. Points c) and e) are placed at similar heights leeward of the roughness element, thus, displaying similar behavior. The main difference being the displayed $\bar{u}\bar{v}$ magnitudes which have dissipated significantly by point e). Presumably, if no other roughness element would occur downstream, Q2 events would again become dominant as similarly theorized by [18] . However, within the selected section there is no visible zone of reattachment. Surprisingly, for $H<3$ , the major contribution in point d) is Q3 followed by Q1. d) also has the lowest recorded measured $\bar{u}\bar{v}$ magnitude for $H=0$ , which is −0.0747 m²/s². For larger H Q2 events become dominant. One possible reason for this might be traces of decelerating high velocity flow from the ridge still present where the point is located.

3.2. Spatial Velocity Correlation

To characterize the size of the flow structures above the rough surface a correlation length approach is utilized. The streamwise spatial velocity correlation is calculated by [36]

$R_u=\frac{1}{{\widetilde{\stackrel{¯}{u}}}^2\left(x_2-x_1\right)}{\displaystyle {\int }_{x_1}^{x_2}\tilde{u}}\left(x,y\right)\tilde{u}\left(x+r,y\right)dx$ (9)

where u = (u, v), r is the streamwise incremental coordinate and

$\tilde{u}=\bar{u}\left(x,y\right)-\frac{1}{x_2-x_1}{\displaystyle {\int }_{x_1}^{x_2}\bar{u}}\left(x,y\right)dx=\bar{u}-\langle \bar{u}\rangle$ . (10)

Using Equation (2), a characteristic length scale can be derived for each velocity component. This operation continues from the crest (not $y/k_s=0$ ) to $y/k_s\simeq 8$ , and an integral length scale is calculated for each acquired auto-correlation function, see Figure 7. The values labeled smooth are measurements near the Plexiglas wall opposite of the rough surface in the center plane

(a)(b)(c)(d)(e)

Figure 6. Instantaneous $\bar{u}\bar{v}$ contributions from the four quadrants (Q_1-4) at the points (a)-(e) visible in Figure 4(b). The sub-figures share a common x-axis and the y-axis is given by Equation (7).

Figure 7. The integral length scale $\tau$ for the u-component.

and are meant to represent a smooth wall case. These values are similar for both the upper and middle case, therefore only one the middle one is displayed.

Near the surface, the flow exhibit very different behavior between the two planes up until about half the channel $\left(y/k_s\simeq 8\right)$ . The middle plane show a sharp increase in the length scale at $y/k_s\simeq 2.7$ , about $0.37k_s$ above the crest of the roughness element with a magnitude of $\tau /{\tau }_r\simeq 1.6$ . The peak, to a degree, indicates the shear layer forming between the bulk and recirculation zone following the roughness peak visible in Figure 4. Scaling the results with the integral length scale of the rough surface, ${\tau }_r$ , as done in Figure 7 provides some interesting insight into the size of the flow structures. There is a strong correlation between the length scales of the rough surface and the flow above the rough surface since $\tau /{\tau }_r\simeq 1$ in the bulk flow, see Figure 7. The obvious explanation for this would be that the rough surface creates similar length scales of the flow above the rough surface. Which would also suggest that the effects of the surface roughness is visible in the entire channel and not just in the vicinity of the wall, a notion proven for similar flow applications [19] [37] . This is a likely hypothesis and the idea is quite intriguing, however, the extent of this phenomenon has to be further investigated before any definite conclusions can be drawn.

3.3. Proper Orthogonal Decomposition

While the integral length scale was applied to the mean velocity, POD was applied to the instantaneous velocity fluctuations. The subsets from the middle plane could not be measured at the same time, hence, any perturbations in the flow cannot be tracked from one subset to another and thereby limiting the visualization of the data. Additionally, the temporal resolution made it difficult to capture sufficient snapshots of the same structure. However, POD is a statistical tool which provides an opportunity to visualize the distribution of energy within each subset, thereby avoiding the problem of synchronized pictures and large temporal resolution. As mentioned in Sec. 3.1, the effects on the mean velocity and higher-order statistics suggest vortex shedding behind the ridge in the

Figure 8. The first three modes from the POD, the left figure is the first mode, the textboxes denote the amount of captured kinetic energy within the mode. The location of the POD is given by the cyan box in Figure 4.

middle plane. To further investigate this, POD was applied to a subset downstream of the roughness element positioned close to the center of the middle field as defined by the cyan box in Figure 4. In Figure 8 the first three modes of the POD are visualized. These three modes capture in total 33.2% of the fluctuating velocity kinetic energy associated with the modes. The modes are dimensionless and the energy content comes from the eigenvalues in Equation (3). Assuming that the POD-mode ends at 20% of the maximum, a longitudinal length scale of about ${\tau }_{pod}/{\tau }_r=1.55$ can be determined, which is similar to the peak in Figure 7. It should be noted that POD-modes cannot be automatically assumed to represent vortices. None the less, it can be applied to zones where shedding is known to appear and identify the dominant structures of the flow, as in this case.

The scale in Figure 8 is between −1 and 1, where 1 corresponds to the maximum energy captured within the mode. The first mode represents the flow structures containing the most energy, which also represents the flow structures rendering the peak in the integral length scale (Figure 7).

4. Conclusion

Results from PIV measurements of flow over a rough hydraulic surface are presented. The surface is produced from laser scanning an existing rock surface and the dimensions of the experimental tunnel are made to reflect real conditions for hydropower tunnels. A likely treatment of such surfaces in the industry is to assume uniform (and thereby small scale) roughness which, according to these results, would lead to erroneous estimations of flow parameters. The results include profiles of double averaged velocity, higher-order statistics, quadrant analysis, correlation length scales and POD. The presented measurements reveal a highly localized behavior of the flow connected to the rough surface. Even small deviations from the local mean height in the surface roughness produce perturbations in the flow which will be visible in the results. In contrast to classical results in asymmetric channel flow the maximum velocity is shifted towards the surface. This shows that the effects from the rough surface are large enough to manifest even in the double averaged velocity. It should be noted that the roughness element studied is not unique, but similar ones occur regularly on the rough surface. Therefore, using large spatial samples would produce the same distortions when averaging. The higher order statistics indicate that the flow above and behind the ridge is characterized by ejection and intermittent bursts of velocity, this is also where the highest point-contributions to the Reynolds stresses where recorded. Similarly, earlier measurements showed higher frequencies of fluctuating pressure at the same ridge. Research has shown that these are unfavorable conditions for rock surfaces from a durability point of view and may hasten or induce an eventual process of tunnel breakdown. As theorized, the problem of tunnel breakdown is likely connected to the flow-roughness effects. Evaluation of the correlation lengths of the flow reveals a significant difference between how the roughness elements interact with the flow. Both the integral length scales and POD approximately predicted the position of the largest flow structures formed in vicinity of the rough surface at $y/k_s\simeq 3$ . Consequently, similar data can be obtained from both the time-averaged velocity and the instantaneous velocity fluctuations. The streamwise length scales of the flow holds close resemblance to the length scales of the rough surface, since $\tau /{\tau }_r\simeq 1$ for $y/k_s>7$ . In contrast to Townsends’s similarity hypothesis where the effects of the rough surface is visible beyond the range of the rough surface, and throughout the channel. Similar effects have been shown in studies concerning flow over riverbeds, dunes or in rivers. It should be noted that such cases usually employ lower Re and larger roughness height to surface ratio and would not regularly be associated with the current application. It does, however, highlight that for hydropower applications, rough surfaces cannot be treated as uniform and only friction-inducing. This is particularly important when modelling flow using CFD, whose role has grown vastly in many industrial applications the past two decades. The results provided here show that the resolution of the roughness is very important, which could have large implications on the future evaluation of hydropower tunnels. A different proposed approach would involve a modification of the uniform wall functions currently used by today’s standards. A thorough mapping of the turbulent kinetic energy in correlation with the rough surface would yield the necessary data, as shown by the quadrant analysis in this study. Thereby the law of the wall could be modified with a spatially stochastic localized increase in wall-near velocity gradients. The roughness length-scales employed in this study might suffice for such an endeavoring as implied by the integral length-scales, however, different methods of deriving these should nevertheless be explored. One limitation of this study is the usage of only one rough wall. Hence, the flow in an actual hydropower tunnel cannot be claimed to be fully understood yet, albeit further understood. There has been no indication of eventual scaling effects between the experiments and the actual case. But if one would consider the problems within hydropower tunnels today and the agreement with other studies the authors assume that the scaling effects, if any, would be insignificant.

Acknowledgements

The research presented was carried out as a part of “Swedish Hydropower Centre-SVC”. SVC has been established by the Swedish Energy Agency, Energiforsk and Svenska Kraftnät together with Luleå University of Technology, KTH Royal Institute of Technology, Chalmers University of Technology and Uppsala University.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Bråtveit, K., Lia, L. and Olsen, N.R.B. (2012) An Efficient Method to Describe the Geometry and the Roughness of an Existing Unlined Hydro Power Tunnelin. Energy Procedia, 20, 200-206. https://doi.org/10.1016/j.egypro.2012.03.020
[2]	Perfect, E. (1997) Fractal Models for the Fragmentation of Rocks and Soils: A Review. Engineering Geology, 48, 185-198. https://doi.org/10.1016/S0013-7952(97)00040-9
[3]	Chanson, H. (2004) The Hydraulics of Open Channel Flow: An Introduction, 2nd Edition, Elsevier, Amsterdam.
[4]	Andersson, L.R., Larsson, I.A.S., Hellström, J.G.I., Andreasson, P. and Andersson, A.G. (2016) Experimental Study of Head Loss over Laser Scanned Rock Tunnelin. 6th International Symposium on Hydraulic Structures, Portland, 27-30 June 2016, 22-29.
[5]	Andersson, L.R., Hellström, J.G.I., Andreasson, P. and Andersson, A.G. (2015) Numerical Simulation of Artificial and Natural Rough Surfacesin. APS 2015 Proceedings, Orlando, 20-25 June 2015.
[6]	Bråtveit, K., Bruland, A. and Brevik, O. (2016) Rock Falls in Selected Norwegian Hydropower Tunnels Subjected to Hydropeaking. Tunnelling and Underground Space Technology, 52, 202-207. https://doi.org/10.1016/j.tust.2015.10.003
[7]	Reinius, E. (1986) Rock Erosion. International Water Power & Dam Construction, 38, 43-48.
[8]	Patel, S.M., Sondergeld, C.H. and Rai, C.S. (2017) Laboratory Studies of Hydraulic Fracturing by Cyclic Injection. International Journal of Rock Mechanics and Mining Sciences, 95, 8-15. https://doi.org/10.1016/j.ijrmms.2017.03.008
[9]	Grass, A.J. (1971) Structural Features of Turbulent Flow over Smooth and Rough Boundaries. Journal of Fluid Mechanics, 50, 233-255. https://doi.org/10.1017/S0022112071002556
[10]	Jiménez, J. (2004) Turbulent Flows over Rough Walls. Annual Review of Fluid Mechanics, 36, 173-196. https://doi.org/10.1146/annurev.fluid.36.050802.122103
[11]	Pope, S.B. (2001) Turbulent Flows. IOP Publishing, Bristol.
[12]	Kruse, N., Kuhn, S. and Rudolf von Rohr, P. (2006) Wavy Wall Effects on Turbulence Production and Large-Scale Modes. Journal of Turbulence, 7, No. 31.
[13]	Cheng, H. and Castro, I.P. (2002) Near Wall Flow over Urban-Like Roughness. Boundary-Layer Meteorology, 104, 229-259. https://doi.org/10.1023/A:1016060103448
[14]	Schlichting, H. and Gersten, K. (2003) Boundary-Layer Theory. 8th Edition, Springer Science, Berlin.
[15]	Seddighi, M., He, S., Pokrajac, D., O’donoghue, T. and Vardy, A.E. (2015) Turbulence in a Transient Channel Flow with a Wall of Pyramid Roughness Journal of Fluid Mechanics, 781, 226-260.
[16]	Krogstad, P. and Antonia, R. (1999) Surface Roughness Effects in Turbulent Boundary Layers. Experiments in Fluids, 27, 450-460. https://doi.org/10.1007/s003480050370
[17]	Patel, V. (1998) Perspective: Flow at High Reynolds Number and over Rough Surfaces—Achilles Heel of CFD. Journal of Fluids Engineering, 120, 434-444. https://doi.org/10.1115/1.2820682
[18]	Bennett, S.J. and Best, J.L. (1995) Mean Flow and Turbulence Structure over Fixed, Two-Dimensional Dunes: Implications for Sediment Transport and Bedform Stability. Sedimentology, 42, 491-513. https://doi.org/10.1111/j.1365-3091.1995.tb00386.x
[19]	Buffin-Bélanger, T., Rice, S., Reid, I. and Lancaster, J. (2006) Spatial Heterogeneity of Near-Bed Hydraulics above a Patch of River Gravel. Water Resources Research, 42, W04413. https://doi.org/10.1029/2005WR004070
[20]	Nikora, V., McEwan, I., McLean, S., Coleman, S., Pokrajac, D. and Walters, R. (2007) Double-Averaging Concept for Rough-Bed Open-Channel and Overland Flows: Theoretical Background. Journal of Hydraulic Engineering, 133, 873-883. https://doi.org/10.1061/(ASCE)0733-9429(2007)133:8(873)
[21]	Zhao, Y., Wang, G.-C. and Lu, T.-M. (2006) Characterization of Amorphous and Crystalline Rough Surface: Principles and Applications. Academic Press, Cambridge.
[22]	Sarkar, S. and Dey, S. (2010) Double-Averaging Turbulence Characteristics in Flows over a Gravel Bed. Journal of Hydraulic Research, 48, 801-809. https://doi.org/10.1080/00221686.2010.526764
[23]	Larsson, I.A.S., Granström, B.R., Lundström, T.S. and Marjavaara, D. (2012) PIV Analysis of Merging Flow in a Rotary Kiln. Experiments in Fluids, 53, 545-560.
[24]	Andersson, A.G., Andreasson, P., Hellström, J.G.I. and Lundström, T.S. (2012) Modelling and Validation of Flow over a Wall with Large Surface Roughness.
[25]	Raffel, M., Willert, C.E., Wereley, S.T. and Kompenhans, J. (2007) Particle Image Velocimetry: A Practical Guide.
[26]	Montgomery, D.C. (2012) Design and Analysis of Experiments.
[27]	Coleman, H.W. and Steele, W.G. (2009) Experimentation, Validation, and Uncertainty Analysis for Engineers. 3rd Edition.
[28]	Balakumar, B.J., et al. (2009) High Resolution Experimental Measurements of Richtmyer-Meshkov Turbulence in Fluid Layers after Reshock Using Simultaneous PIV-PLIF. Proceedings of the APS Topical Group on Shock Compression of Condensed Matter, Nashville, December 2009, Volume 1195. https://doi.org/10.1063/1.3295225
[29]	(2007) LaVision Gmbh Product Manual for Davis 7.2.
[30]	Westerweel, J. and Scarano, F. (2005) Universal Outlier Detection for PIV Data. Experiments in Fluids, 39, 1096-1100.
[31]	Bakewell, H.P. (1967) Viscous Sublayer and Adjacent Wall Region in Turbulent Pipe Flow. The Physics of Fluids, 10, 1880.
[32]	Erik Meyer, K., Cavar, D. and Pedersen, J.M. (2007) POD as Tool for Comparison of PIV and LES Data.
[33]	Yue, W., Meneveau, C., Parlange, M.B., Zhu, W., Van Hout, R. and Katz, J. (2007) A Comparative Quadrant Analysis of Turbulence in a Plant Canopy. Water Resources Research, 43, W05422. https://doi.org/10.1029/2006WR005583
[34]	Hanjalic, K. and Launder, B.E. (1972) Fully Developed Asymmetric Flow in a Plane Channel. Journal of Fluid Mechanics, 51, 301-335. https://doi.org/10.1017/S0022112072001211
[35]	Nakagawa, S., Na, Y. and Hanratty, T. (2003) Influence of a Wavy Boundary on Turbulence. I. Highly Rough Surface. Experiments in Fluids, 35, 422-436. https://doi.org/10.1007/s00348-003-0681-2
[36]	Tennekes, H. and Lumley, J.L. (1972) A First Course in Turbulence.
[37]	Roy, A.G., Buffin-Bélanger, T., Lamarre, H. and Kirkbride, A.D. (2004) Size, Shape and Dynamics of Large-Scale Turbulent Flow Structures in a Gravel-Bed River. Journal of Fluid Mechanics, 500, 1-27. https://doi.org/10.1017/S0022112003006396