Modelling roughness and acceleration effects with application to the flow in a hydraulic turbine

This study reports the numerical predictions of flows over turbine blades, which include flow acceleration and deceleration. Two issues are addressed: (1) accurately predicting roughness effects, and (2) evaluating the performance of Reynolds-Averaged Navier-Stokes (RANS) simulations on moderately accelerating flows. For the present turbine surfaces, it is found that roughness correlations based on roughness surface slope better predict the roughness effects than both the correlations based on the moments of roughness height statistics and the IEC standard approach. It is shown that RANS simulations reproduce the flow evolution over rough-wall accelerating turbulent boundary layers, although, on a smooth wall, they fail to capture strong non-equilibrium flow behaviours. Finally, a hydraulic turbine simulation is performed to show the significant roughness impact on the total losses.


Introduction
Hydraulic turbines have long lifespans; the oldest machines in operation at Hydro-Québec are over 100 years old. Over a long lifespan, surface degradation occurs and significant roughness contaminates the solid surfaces: corrosion, impacts from debris, cavitation, and cracks alter the initial smooth finish. While mitigating measures (e.g., surface smoothing) exist, they are costly and it is difficult to calculate the return on investment due to several technical problems from an engineering standpoint. This paper reports partial results of a research project aimed at developing accurate numerical methods capable of predicting the flow over realistic rough surfaces, using engineering models. This type of model would allow engineers to predict when it is financially beneficial to proceed to blade smoothing by quantifying the losses due to the roughness.

Roughness effects
The first study on roughness was performed by Nikuradse [1] in long pipes covered with sand grain of various heights. Colebrook [2] extended the results by fitting a formula appropriate for industrial applications, and Moody designed an easy-to-use chart still taught today in engineering classes. Schlichting [3] suggested using an "equivalent sand grain height", k s , to describe the effects of different rough surfaces by comparing them to the canonical sand grain. The main effect of roughness is an augmentation of drag, which is reflected in the mean velocity profile as an offset from the smooth-wall profile in the logarithmic region; this offset is called the "roughness function", ΔU + , where the subscript "+" indicates normalization using the kinematic viscosity, ν, and the friction velocity, u τ [4]. In solvers of the Reynolds-Averaged Navier-Stokes (RANS) equations, the single parameter used to describe the roughness is k s . What is usually measured, however, is the arithmetic average height, Ra. The ASME standard [5] provides some guidance to measure the roughness of machined surfaces. The drag, or k s , depends on both Ra and the surface texture. Many studies have been carried out to find a relation between geometrical characteristic of the surface and k s , but most were application specific [6]. Acknowledging this, the IEC code [7] recommends to use a simple linear relation, k s = 5Ra; but, as will be shown later on, the actual k s may differ from this prediction. Note that the linear relation between k s and Ra applies only in the fully rough regime, in which the viscous sublayer in the turbulent boundary layer is completely destroyed due to roughness.
Correlations have been proposed to predict the ratio k s /Ra solely based on the geometric characteristics of the surface. Reviews are found in Bons [8] and Flack and Schultz [9]. A set of commonly used correlation is based on the surface slope. For instance, Sigal and Danberg [10] proposed where and S, S f , and S s are the reference surface area before adding roughness, the total frontal area of the roughness, and the total windward wetted surface, respectively. Other forms of Λ s have been proposed for regular 2D or 3D roughness [11,12], but equation (2) is best suited to realistic surfaces. The constants a and b are empirical coefficients. Another slope-based correlation was proposed by Bons [8], where α is the local streamwise slope angle, which can be easily obtained from a 1D surface trace. A second class of correlations is based on the moments of height statistics, including the second (Rq), the third (skewness, Sk), and sometimes higher-order moments (kurtosis, Ku). For example, Flack and Schultz [9] proposed a correlation that was shown to work well on several types of realistic surfaces (gravel, honed surfaces, etc.): The "effective slope" (ES) of a surface is an additional parameter that can be used to categorize the roughness types with regard to the importance of surface slope. It is defined as the average magnitude of local surface slope, where k(x) is the surface height array (obtained from 1D surface trace), x the streamwise direction, and L is the sample length in x. It was shown that the surface slope affects ΔU + only when the surface is not sufficiently steep, i.e., in the waviness regime; when the surface is steep, it is said to be in the roughness regime [9,13].

Pressure-gradient effects
In hydro turbines, large areas of mild to moderate flow acceleration and deceleration occur frequently ( Figure 1), characterized by favourable (for acceleration) or adverse (for deceleration) pressure gradients (FPG or APG). The acceleration parameter, K, is usually used to quantify the level of acceleration/deceleration; here, K is defined as K=(ν/U 2 )dU/dx along streamlines, where x is the flow direction and U is the flow velocity. RANS solvers and wall functions are widely used in turbine simulation to control the calculation cost. Roughness effects on the turbulence are introduced by setting the boundary value of eddy viscosity or by setting u τ in the wall models (to produce ΔU + ). The existing rough-wall modifications are usually calibrated in zero-pressure-gradient (ZPG) flows. It is not clear how these modifications behave under combined effects of roughness and pressure gradients.

Objectives
In this study, we investigate the accuracy of multiple approaches to determine k s and the performance of RANS in flows with moderately accelerating flows with or without roughness. First, we carried out large-eddy simulations (LES) in a turbulent open-channel flow to obtain the actual k s , to investigate the accuracy of the IEC standard and existing k s -correlations (section 2). Then, we investigate the performance of RANS models by comparing them with the LES on FPG flat-plate turbulent boundary layers over smooth and rough walls (section 3). Finally, the fully rough flow in the distributor and the runner of a Francis turbine is simulated using RANS, and the sensitivity of the results on the value of k s is studied (section 4).

Problem formulation
x 1 , x 2 , and x 3 (or x, y, and z) are, respectively, the streamwise, wall-normal and spanwise directions, and u i (or u, v, and w) are the filtered velocity components in those directions; P = p/ρ is the modified pressure, ρ the density. The sub-grid stress τ ij is modeled using the Lagrangian-averaged eddyviscosity model [14]. The overbar denoting the filtering operation is hereafter omitted for simplicity. The simulations are performed using a well-validated code that solves (6) on a staggered grid using second-order, central differences for all terms, and a second-order semi-implicit time advancement scheme. An open-channel flow with half-channel height h is simulated with no-slip and symmetric boundary conditions applied to the bottom wall and to the top boundary, respectively. Periodic conditions are used in both the x-and z-directions. An immersed-boundary method [15] is used to impose no-slip boundary condition on the rough surfaces.  Roughness samples were obtained from real turbines. Figure 2 shows the samples S1 to S4 used in the present work; these surfaces are representative of those taken on various sites. Digitalization and measurement were carried out using an optical profilometer CHR150 from Microphotonics. The surfaces are different in the characteristics ( texture [5] is tentatively applied for the heavily rusted surfaces; a Gaussian filter was used to separate the surface oscillations from roughness. Table 1 shows that the filtering length significantly affects the measured roughness height. In this study, we considered the unfiltered surfaces. Since the size of surface sample is smaller than the horizontal simulation domain, the sample is called a "tile"; the final surface is obtained from duplicating the tile of the same type with random rotations, to form a horizontal domain of size at least 6h×3h. Visual inspection of multiple rough samples confirms that no noticeable orientation of the roughness elements is present. A linear smoothing function is applied to the interface between two patches to ensure a seamless match between tiles. A numerical sand-grain (SG) roughness model [16] is also considered. The kurtosis values for final surfaces S1 to S4 and the SG are close to 3 (table 2), indicating that randomness is achieved.
The domain size in x and z ranges from 6h×3h to 12h ×12h, sufficient to accommodate the largest turbulent structures with the current Reynolds number. To cover both transitionally and fully rough regimes, two Reynolds numbers and three Ra values are used: the Reynolds number based on u τ and h are Re τ =400 and 1000, and Ra values are 0.020h (K1), 0.040h (K2), and 0.067h (K3). The grid sizes Δx + and Δz + range between 3 and 40, similar to the values used in previous rough-wall simulations using the same numerical techniques [15]; the resolution in y satisfies Δy + <1 in the region y<Rt and Δy + <25 at the top boundary. The total number of grid points ranges from 128×128 to 640×512 in x and z, and from 82 to 300 in y.

Results
The profiles of U + for cases with K3 and Re τ =1000 are shown in figure 3(a). The zero-plane displacement, d ≈ Ra, is subtracted from the wall-normal location to collapse the logarithmic regions of the smooth and rough cases [15]. The smooth-wall profile collapses with the experimental data by Schultz and Flack [17]. In the logarithmic region, the rough-wall profiles display an offset from the smooth-wall one, with the amount of the offset, ΔU + , depending on the surface topography.
Since the linear relation between k s and Ra is valid only inside the fully rough regime, the beginning of this regime needs to be first identified: it is determined by varying Ra + and found as the start of the logarithmic relation between ΔU + and Ra + . The fully rough regime starts at a lower roughness height for the turbine surfaces (especially S1 and S4) than for SG, where it is established for k + s,crit ≈80.
In the fully rough regime, k s is obtained from ΔU + for each simulation [1], where κ is the Von Kármán constant. The results are tabulated in table 3; k s /Ra varies significantly (from 0.3 to 2.7), and the IEC standard (k s /Ra=5) over-predicts k s for the current surfaces. Figure 3(b) shows the relation between ΔU + and ES with a fixed Ra. An overall increase of ΔU + is observed as ES increases, indicating that the current surfaces are in the waviness regime. Also plotted Experiment SG S1 S2 S3 S4 Figure 3. Results in cases with K3 and Re τ =1000: (a) streamwise mean velocity; + Schultz and Flack [17]; (b) dependence of the roughness function on the effective slope; + Napoli et al. [13]. The fully rough data (cases with K3 and Re τ =1000) are used to evaluate the existing k s -correlations specified in equations (1) (slope/shape method), (3) (slope-rms method), and (4) (moment method). For each correlation, the constants a and b are fitted with the actual k s values obtained from LES, and the resultant correlation is plotted in figure 4. The forms of the two slope-based correlations (figure 4(a) and 4(b)) represent reasonably well the variation of the current data. Larger scatter is found when the moment method is used (figure 4(c)). In particular, surfaces S3 and S4 have similar Rq/Ra (within 22%) and Sk (within 20%) as shown in table 2, but differ by a factor of 4 in k s /Ra. Since moments of height statistics do not contain slope information directly, the moment-based correlation is not suitable in the waviness regime. It is suggested that the slope-based correlations be used to predict k s in hydraulic-turbine applications.

Problem formulation
To study the performance of turbulence models for the RANS equations, both 3D LES and 2D RANS are carried out for the same flat-plate accelerating boundary layer configuration: a streamwise freestream velocity U ∞ (x) varying in the streamwise direction is imposed on the top boundary of the domain; the vertical freestream velocity is obtained by ensuring continuity. Figure 5(a) shows that the acceleration parameter K, matches the one used in earlier experimental studies [18]. The reference quantities are the displacement thickness δ o * and the freestream velocity U ∞,o at a reference location x=0, where the Reynolds number based on U ∞ (x) and δ * (x) first reaches 1240; in most part of the domain, this Reynolds number ranges between 1000 and 2000 for the smooth case, and between 2000 and 6000 for the rough case; it is lowered by acceleration, and increased by roughness and boundarylayer growth. For LES, an unsteady inflow boundary condition was obtained from an inlet region of 60δ o * using the recycling/rescaling method [19], and a convective outflow boundary condition was used at the domain outlet. For RANS, an inlet region of around 400δ o * is used for the boundary layer to develop from uniform velocity inlet (U=U ∞,o , V=0). Around 100 millions grid points are used for the LES, which gives 20<Δx + <50, 10<Δz + <20, and Δy + <1 for y<Rt; for RANS, around 100 thousands grid points are sufficient to resolve the flow, with the first vertical grid point located below y + =1. For the rough case, fully rough regime is achieved throughout the domain; k s + ≈70 at x = 0. To study the RANS behaviors in non-equilibrium flow, the k-ω model [20] with roughness modification proposed by Fuhrman et al. [21] is compared to LES. Also considered are two models in CFX, a commercial code widely used for turbine analysis: the shear stress transport (SST) k-ω model (with or without wall function) and k-ε model (with wall function).

Results
On the smooth wall, the friction coefficient C f from the LES is validated against experimental results ( figure 5(b)). C f increases at the beginning of acceleration due to a faster increase of u τ compared to U ∞ ; as K strengthens, the turbulence starts to relaminarize [22,23]; during this process, C f decreases and the mean velocity profile becomes laminar-like. However, the RANS model does not capture this strongly non-equilibrium process. Figure 5(b) shows that, on a smooth wall, C f from both RANS simulations varies uniformly with K, without showing a dip. The velocity profiles in figure 6(a) shows that the LES and RANS results match well in the equilibrium region (x/δ o * =50); in the strong-K region (x/δ o * =220), both results show a decrease of logarithmic slope, an effect of FPG [22]; in the relaminarization region (x/δ o * =300), however, LES gives a laminar-like mean flow, while RANS predicts an instantaneously recovery from the acceleration. The results from the SST model in CFX are consistent with the k-ω model.
In the rough-wall case, relaminarization does not occur, due to the breakup of the viscous sublayer in the fully rough regime [23]; thus, compared to a smooth-wall flow, the rough-wall FPG flows are expected to be predicted more accurately by RANS. Figure 5(c) shows that the C f from RANS matches the trend of the LES results: an increase of C f in the non-equilibrium region (100<x/ δ o * <300). The k-ω model without wall function better captures this effect of FPG, while the models with wall functions (both SST and k-ε) under-predict it. The velocity profiles in figure 6(b) illustrate this difference: the models with wall function collapse with LES and the k-ω model in the ZPG region (x/δ o * =50 and 300), while they slightly under-predict the offset in the high-FPG region (x/δ o * =220). Nevertheless, such effect of the wall function is acceptable for mildly or moderately accelerating flows. x

Roughness effects on flows around realistic turbine geometry
This section evaluates what impact the roughness produces on turbine performance using commercial software. Three simulations were performed with CFX in the hydraulic turbine presented in figure 2; the turbine geometry induces flow acceleration and deceleration, as shown in figure 1. Different approaches are used to characterize roughness. Table 4 shows the values of k s obtained from the LES simulations and from ASME/IEC standard; k s varies significantly, up to one order of magnitude for the guide vane and runner. For the distributor lower and upper surfaces, S1 is assumed since no measurement is available.
The numerical setup has been validated to obtain a performance hill chart for this turbine. The head is imposed at the inlet, and a static pressure is placed at the outlet. The simulation includes a periodic sector of the distributor, the runner, and a Moody draft tube. The k-ε model with the scalable wall function is used. k s is the only variable in the simulations.
where C is the velocity, μ and μ T the molecular and eddy viscosities, Ω the volume, and A is the inlet area. Figure 7 shows that the losses clearly correlate with the roughness height. Results presented in table 5 show that the roughness impact on the efficiency is between 1 and 2.5% (sum of runner and distributor losses); the values of k s obtained from LES give 1.5% higher efficiency compared to those obtained from ASME/IEC standards. The flowrate and power are also affected. Therefore, to obtain an accurate estimate of the absolute efficiency for existing machines, roughness should be taken into account, and accurate k s values must be used.

Conclusion
This study investigates the approaches currently applied to predict turbine performance using RANS simulation with the value of k s predicted from the surface geometry. The actual values of k s on multiple turbine-blade surfaces are obtained using LES in an open-channel flow. For the current turbine roughness, which is in the waviness regime, it is found that the IEC standard over-predicts k s ; the slope-based k s -correlations best predict k s , while the moment-based method shows larger scatter.
To validate the RANS models on non-equilibrium boundary layers over roughness, a flat-plate turbulent boundary layer under moderate freestream acceleration is simulated using both LES and RANS; RANS captures the flow evolution in the fully rough flows; however, it does not predict the strongly non-equilibrium flow reversion in the smooth case. The wall function is shown to slightly under-predict the drag increase in the high-FPG region.
The importance of k s is shown using a RANS simulation of flow over an actual turbine geometry, which results in moderate FPG and APG levels. It is shown that different k s -prediction approaches provide significantly different losses. This demonstrates the importance of the roughness effects and properly selecting k s -prediction method; it also provides further validation for the use of RANS model in hydraulic turbine simulations.