Unsteady Numerical Simulation of Flow in Draft Tube of a Hydroturbine Operating Under Various Conditions Using a Partially Averaged Navier–Stokes Model

Hydropower is a renewable and sustainable energy source with many economic and environmental benefits. Storage capability and fast response make hydropower an ideal form of power generation. Hydropower plants provide 16.5% of the electricity worldwide, making hydropower the largest source of renewable electricity generation [1]. However, the variable energy demand requires hydroturbines to be operated over an extended range of conditions, sometimes far from their design point where the efficiency of the turbine is at its maximum. Several undesirable behaviors such as flow blockage, cavitation, and pressure and output power fluctuations may occur when a hydroturbine, especially a Francis turbine, operates under off-design conditions [2]. They result in efficiency reduction, structural vibration, and power swing, and reduce the optimal operating range of the machine.

One of the main components of a hydroturbine where unwanted flow phenomena appear under off-design conditions is the draft tube. The draft tube is the last component of a hydraulic turbine where the flow exiting the runner is decelerated, thereby converting the excess of kinetic energy into static pressure. It increases the efficiency of the plant at BEP and makes it possible to locate the turbine above the tail-water. However, as shown in Fig. 1, the draft tube displays an abrupt increase in hydraulic losses as the operating regime departs from the BEP, while other components (e.g., stay vanes, guide vanes, and runner) have a rather modest variation with the operation regime [3]. Therefore, it is admitted that for modern Francis turbines the shape of the performance chart (hill-chart) is dictated by the losses in the draft tube [4].

Theoretically, the flow leaving the runner and ingested by the draft tube is almost purely axial at the BEP. However, in practice, to minimize hydraulic losses associated with kinetic energy-to-static pressure conversion in the draft tube, a certain level of residual swirl is provided at the runner outlet. This swirling flow at the draft tube inlet is tuned for optimal performance at the best efficiency operating point. However, the swirl ingested by the draft tube departs significantly from the best configuration as the turbine discharge varies, and hydroturbines, especially Francis turbines, operating at off-design conditions often have a high residual swirl at the draft tube inlet.

The interaction between the decelerated swirling flow and the wake of the crown cone may lead to efficiency reduction and severe pressure fluctuations in the draft tube [5]. These pressure fluctuations become more dangerous, if their corresponding frequency approaches the natural frequency of the power plant structures. This can result in vibration of the whole installation [2]. Depending on the operating conditions, two forms of the flow instability can be observed. At part-load conditions, for which the flow rate is lower than the one at the BEP, the flow has a positive absolute circumferential velocity in the same sense as the runner revolution. In this case, a helical precessing vortex called the “vortex rope” develops in the draft tube with the frequency of precession of about 0.2–0.4 of the runner's rotation frequency [6]. This helical vortex is formed due to the rolling up of the shear layer between the wake of the crown cone (also called the stagnant region) and highly swirling outer flow [5]. At full-load conditions corresponding to higher flow rate than that of the BEP, the absolute circumferential velocity is negative inducing a swirling flow rotating in the opposite direction of the runner. In this case, the vortex rope takes a nearly axisymmetric shape, sometimes called the “torch.” In both cases, cavitation may occur in the vortex core.

Understanding of the physics of the flow in the draft tube allows engineers to avoid the above-mentioned undesirable phenomena. The knowledge gained can be used to elimination of these phenomena, and thereby enhancing the operating range of hydropower plants, adding significantly to the power generation economy. The low cost and availability of detailed data sets of make CFD an important tool for investigation of flow in hydroturbines, especially in draft tubes. The numerical simulation of unsteady turbulent flow in a draft tube, however, is a challenging task, and accurate numerical calculations of the flow parameters require a careful choice of turbulence closure method. Several previous draft tube flow predictions have been carried out employing the RANS equations with various turbulence closure models [7–10]. Although there are cases of successful draft tube simulation using two-equation turbulence RANS models, these models generally fail because of the isotropy assumption which is not suitable for such strongly swirling flows with highly curved streamlines [11]. For example, Vu et al. [12] performed simulations of flow in a draft tube using the k-ε turbulence model, and stated that as the operating conditions depart from the BEP, the performance of the k-ε model tends to deteriorate. Additionally, it is shown [5,13] that unsteady RANS (URANS) models overdamp the flow unsteadiness in the draft tube and give steady solutions without resolving unsteady coherent structures of the flow. As discussed in detail by Foroutan and Yavuzkurt [5], this is due to the fact that the time-scales of draft tube instabilities are close to that of large turbulent energy carrying eddies, and hence, they are being averaged-off by URANS models.

An alternative solution is to apply hybrid URANS/large eddy simulation (LES) models. Hybrid URANS/LES methods aim at combining the best of URANS and LES turbulence models. In this approach, the URANS models are mainly employed in the near-wall region, while LES treatment is applied to regions far away from the wall. Hybrid URANS/LES has become one of the main modeling framework for quantitatively accurate predictions of complex unsteady flows at high Reynolds numbers. Various hybrid URANS/LES models have been used for unsteady draft tube flow simulations including very large eddy simulation (VLES) [13], the detached eddy simulation [14], and the scale-adaptive simulation [15]. Most of these models have shown relative improvement over the URANS models. However, detailed evaluation and investigation of these models for various operating conditions of the hydroturbines are seldom addressed.

Foroutan and Yavuzkurt [16] developed a new hybrid URANS/LES turbulence model in the PANS framework, which can be used in unsteady simulations of turbulent swirling flows. This paper presents results obtained using this PANS model in simulations of flow in the draft tube of a Francis turbine operating under various conditions. Therefore, the present paper tries not only to validate the newly developed PANS model [16] for a wide range of operating conditions but also to discuss and compare physical phenomena occurring under each condition. In addition, the present PANS results are compared to the computational results obtained using RANS models and capability/incapability of these models for various operating conditions are discussed.

The PANS [17] method is a bridging, seamless model between RANS and direct numerical simulation. It aims to provide the best possible closure at any given level of computational grid size, and improve accuracy with increasing the resolution. The PANS is advantageous in comparison with the LES, since the grid does not need to be fine enough to fulfill the requirement of LES; therefore, less number of grid points and less computational resources would be needed. Furthermore, in PANS modeling, a RANS turbulence model is used for the near-wall region, which would be much more effective than a simple subgrid scale model used in LES simulations. The seamless, continuous nature of the PANS model is a results of the so-called “partial averaging” concept [17], which corresponds to averaging (filtering) a portion of the fluctuating scales. This portion is based on the ratios of unresolved-to-total turbulent kinetic energy and dissipation rate.

$νu$ denotes the PANS eddy viscosity (eddy viscosity of unresolved scales) and obtained as a function of the unresolved turbulent kinetic energy k_u and unresolved turbulence dissipation rate ε_u [17]. Therefore, in order to close the system of PANS equations, two transport equations for the unresolved turbulent kinetic energy k_u and its dissipation rate ε_u are required.

The model coefficients are given in Table 1. The two key parameters in PANS modeling are f_k = k_u/k and f_ε = ε_u/ε, which represents the unresolved-to-total ratios of turbulent kinetic energy and turbulence dissipation rate, respectively. The unresolved-to-total kinetic energy ratio f_k was prescribed as a constant in the original PANS model [17]. However, using partial integration of the complete turbulence energy spectrum, Foroutan and Yavuzkurt [16] developed a new relation for the f_k parameter as follows:

where $Δ=(ΔxΔyΔz)1/3$ is the grid length scale, and $Λ=k3/2/ε$ is the turbulence integral length scale. When the grid is coarse (⁠$Λ≪Δ$⁠), Eq. (9) gives an f_k value close to 1, and the present PANS model returns to its background RANS model. When the grid size is much smaller than the local turbulence length scale, more and more turbulent kinetic energy can be resolved as f_k decreases. The parameter f_k in Eq. (9) is a dynamic parameter, changing at each location at the end of every time-step, and then it is used as a fixed value at the same location during the next time-step, as in Han et al. [19]. The parameter f_ε represents the unresolved-to-total ratio of the turbulence dissipation rate. In a high Reynolds number turbulent flow, most of the dissipation is associated with small scale structures, which are much smaller than the grid size used in practical applications. Therefore, it is unlikely that the dissipative small scales to be resolved; hence, the parameter f_ε is taken to be equal to 1, as in Ref. [19], implying that ε_u = ε.

One of the challenges any CFD study of draft tube faces is the lack of a test case for which the geometry and boundary conditions, as well as detailed experimental data are available. These experimental data are needed for assessment and validation of numerical simulations. The measurements should include not only pressure fluctuations but also profiles of velocity components as well as turbulent quantities. Among the several experimental studies reviewed for this study, those related to the FLINDT project [20] contain these details. The FLINDT project experiments were carried out on a scaled model (1:10) of a Francis turbine (see Fig. 2) with specific speed of 0.56. The turbine model has a spiral case, stay ring of ten stay vanes, a distributor made of 20 wicket gates (guide vanes), a 17-blade runner of 0.4 m outlet diameter, and an elbow draft tube with one pier [20]. Many experimental and numerical studies have been carried out under the FLINDT project [6,7,9,12]; however, the exact details of the draft tube geometry is not available in open literature. Therefore, the simplified FLINDT draft tube has been investigated by some researchers [21,22]. A comprehensive investigation of the previously published articles within the FLINDT project is performed in this study, in order to build a complete database on details of the draft tube geometry (as much as possible) and available experimental data. Using this database, the three-dimensional FLINDT draft tube geometry was rebuilt as shown in Fig. 3. Therefore, this study can provide a publicly available test case with all details required for a CFD simulation.

The draft tube geometry includes the runner crown cone, the draft tube cone of 17 deg angle followed by a 90 deg curved elbow, and a rectangular section diffuser with a pier. Section S0 in Fig. 4 is located just downstream of the trailing edge of runner blades. Sections S1 and S2 are the locations within the draft tube cone where radial profiles of velocity components were measured in experimental studies [7,9]. They are located 84 mm and 242 mm downstream of the runner outlet, respectively. The static pressure was measured at the draft tube wall at sections S1 to S5 using piezoresistive pressure transducers [7,23].

According to the objective of this study, four cavitation-free operating points (case A–D) are selected covering a wide range of operating conditions. Table 2 summarizes the characteristics of these four operating points. The points of interest are selected for the same head (energy) coefficient of ψ = 2gH/($ω2R2$) = 1.18 and different flow rate coefficients of φ = Q/($πωR3$) = 0.41, 0.368, 0.34, and 0.26 corresponding to about 110%, 99%, 91%, and 70% of the flow rate at BEP.

For case D, the circumferentially averaged mean velocity components and turbulent quantities (k and ε) just downstream of the trailing edge of runner blades are available [24]. Therefore, the inlet section for this case is chosen to be the surface swept by the trailing edge of the runner (section S0 in Fig. 4). These data, however, are not available for the cases A–C; hence, a downstream section in the draft tube cone (section S1 in Fig. 4) is considered as the inlet section where the experimental data, including axial and circumferential velocity components and turbulent kinetic energy profiles, are available. Because of the relatively high uncertainty and low magnitude of the radial velocity component in measurements at section S1, a linear variation of this component is considered [22] at inlet for cases A–C. Turbulence dissipation rate is not available at section S1 from experiments. Therefore, the inlet profiles for the dissipation rate are computed from the turbulent kinetic energy profiles as $ε=k3/2/l$ [25] for cases A–C. The turbulence length scale $l$ = 0.01R (R is the runner radius) has been chosen according to previous studies [22], which are based on the analysis of laser Doppler velocimetry measurements.

Figure 5 illustrates profiles of velocity components applied at the inlet boundary, i.e., section S1 for cases A–C and section S0 for case D. The evolution of the axial velocity component at the inlet for cases A–C (Fig. 5(a)) is characterized by increasing the low velocity region at the center as turbine's flow rate decreases. Also, the circumferential velocity increases as the operating condition deviates from the BEP. However, the direction of the circumferential velocity (i.e., the direction of the swirl) is inverted for the operating points at the higher flow rates (Fig. 5(b)). This can be explained by considering the velocity triangles at runner exit (see Ref. [26] for details). Figure 5(c) shows the axial, radial, and circumferential velocity profiles at section S0 (the inlet of for case D). The circumferential velocity is quite high in this case resulting in a high level of swirl (swirl number of 0.63) at the inlet to the draft tube. The increase in velocity components near the band (see Fig. 5(c)) in this case is related to the secondary flows in the blade channel at partial load with formation of interblade vortices [27].

In order to overcome back-flows at the outlet, which could introduce numerical instabilities, the outflow boundary is located farther downstream after the draft tube discharge. Nevertheless, it has been shown [28] that the outlet boundary plays a relatively less important role in this case. No-slip conditions are applied on all walls, while the crown cone is defined as a rotating wall with an angular velocity equal to that of the runner (for case D with S0 as the inlet section).

The main computational grid for cases A–C (where the computational domain starts from section S1) consists of 1,324,220 cells and for case D (where the computational domain starts from section S0) consists of 2,284,220 cells, which are quite moderate grid sizes. This, in fact, is one of the advantages of the PANS model where the requirement of the computational grid is not as rigid as for an LES model and a relatively coarser mesh can be used. The grid is refined in areas of large variable gradients, i.e., near-wall and centerline region. To make sure that a grid-independent solution is obtained, simulations are also performed for case A with a finer grid of 3,716,064 cells (results are discussed in Sec. 5). The first cell center normal to the wall is placed at y⁺$≈$ 2 everywhere in both grids.

The governing equations (Eqs. (1), (2), (5), and (6)) are numerically solved using finite volume method in ansys-fluent 14.0 [29], where the newly developed PANS model is implemented by means of a user-defined function. All simulations are performed using a bounded central-differencing scheme for the momentum equation, while a second-order upwind scheme is used for turbulence quantities. The temporal differencing is carried out using the second-order implicit scheme. The SIMPLE algorithm is used to couple the velocity and pressure fields. Unsteady PANS simulations are initialized by a k-ε RANS steady solution. The time-step size is taken corresponding to one degree rotation of the runner (2.2$×$10⁻⁴ s) which is known to be sufficient for hydroturbine applications [9]. The convergence criterion of the residuals for each time-step is set to three orders of magnitude drop or maximum 30 subiterations. In order to make sure that a quasi-steady condition is reached, unsteady simulations are continued for a very long time. Specifically, simulations are performed for 73,000 iterations corresponding to about 200 rotations of the runner, and about 37 through-flow time (the time required by the mean flow to pass through the domain once). All time-averaged parameters are computed by temporal averaging of the results of 60,000 time-steps (13.2 s of operation of the runner corresponding to 167 rotations) after setting the simulations to run for an initial 13,000 time-steps (2.86 s of operation of the runner corresponding to 36 rotations), i.e., initial transient. On average, each time-step takes about 27 s (48 s for the case D) for the main grid and 82 s for the fine grid on a 16-processor Linux cluster. Simulations are performed on a high performance computing cluster within the Pennsylvania State University's Research Computing and Cyber-infrastructure (RCC) unit.

Two global parameters are considered for overall evaluation of the present numerical simulations:

1. (1)
   The pressure recovery coefficient defined as [7]

   $χ≡p6−p212ρ(QA2)2$

   (10)

   and calculated between sections S2 and S6 in Fig. 4. The pressure recovery coefficient quantifies the overall hydraulic performance of a turbine draft tube in converting the excess of kinetic energy at the runner outlet into static pressure.

2. (2)

   The portion of the flow exiting the “left” channel (see Fig. 3 for the definition of left and “right”).

Table 3 shows the draft tube pressure recovery coefficient obtained with the present simulations in comparison with the experimental data [23]. Numerical results obtained by Mauri [7], who used the standard k-ε turbulence model, are also presented for comparisons (where available). As expected, the highest pressure recovery is achieved at the BEP, while the pressure recovery coefficient drops dramatically as conditions moves away from the BEP. The k-ε RANS model gives good results at the BEP, i.e., 3.2% difference with experimental data. The draft tube flow at the BEP has very moderate level of swirl and is quite “well-behaved.” No strong flow instability and vortex rope exists in the flow; Therefore, RANS closure models are able to give reasonable predictions. Even in case A, with 10% more flow rate than the BEP, the flow instability is not high enough to form the vortex rope and RANS simulations are able to predict the pressure recovery coefficient by only 3.7% deviation from the experimental data values. It should be noted, however, that in both cases A and B, predictions of the present PANS model are in very good agreement with experimental data (only 2.6% and 0.6% difference in case A and B, respectively), while improvements are seen comparing to the k-ε RANS predictions. Nevertheless, as shown previously in Ref. [12], the clear drawback of the RANS model occurs under partial load conditions, where flow is highly unsteady and a strong helical vortex rope forms inside the draft tube. The RANS closure models are unable to predict the flow behavior in these conditions, while there still exists a good agreement between the predictions of the PANS model, which resolves important vortical structures of the flow, and experimental data. Specifically, deviation more than 13% is seen between experimental data and RANS simulations for case C, while PANS predictions show less than 4% deviation. For case D with 70% of the BEP flow rate, no pressure recovery calculation is found in Ref. [7], therefore, RANS simulations are performed which show more than 58% overprediction. Result of the present PANS simulations deviates by only 6% from experimental data for this case. Using the fine grid (with 3,716,064 cells), results of the present PANS model show a slight improvement and only 1.4% difference with experimental data is seen for the case A. However, the quite considerable increase in the computation time (more than three times) cannot be justified. Therefore, results obtained using the main grid (with 1,324,220 cells) is considered to be “grid-independent” within a reasonable error.

Table 4 presents the portion of the flow rate through the left channel, where the present PANS predictions are compared with experimental data as well as those obtained by the k-ε RANS simulations [7]. Again, a fairly good agreement is seen between present PANS predictions and experimental data, with PANS results for case A with the fine grid being in excellent agreement with data. Another interesting point is the uneven partition of the flow rate in the left and the right channels, which become more distinct at lower flow rates. It is seen that at 70% flow rate, only 20% of flow exits through the right channel. This is physically due to the interactions between swirling flow entering the draft tube and secondary flows due to the curvature of the draft tube elbow. A numerical test performed in this study indicates that reversing the swirl direction at inlet (i.e., reversing the runner's rotation direction) switches the portion of flow going through each channel.

Figures 6–9 shows the time-averaged axial and circumferential velocity components in the draft tube obtained from the present PANS simulations. For cases A–C, with section S1 as the inlet, results are plotted on section S2 (Figs. 6–8). Case D has the benefit of having the inlet at section S0; therefore, results can be plotted at two sections, i.e., S1 and S2 (Fig. 9). Results are compared to the experimental data [7,9,22] as well as those obtained from the k-ε RANS simulations. As shown by previous studies [5,13,14] and discussed here in Sec. 1, URANS models damped out the unsteadiness of the flow and give steady results. Hence, time-averaged k-ε URANS results found to be the same as the k-ε RANS steady results.

Comparing Figs. 6(a), 7(a), 8(a), 9(a), and 9(c) reveals that the axial velocity component is characterized by increasing the low velocity region (stagnant region) at the center of the draft tube as the turbine's flow rate decreases. Furthermore, comparing Figs. 9(a) and 9(c) show that this region radially expands towards downstream in the draft tube. Comparing the circumferential velocity shows that the minimum swirl is associated with the BEP flow rate. Also, the swirl is damped moving downstream in the draft tube (compare Figs. 9(b) and 9(d)).

As shown in Fig. 6, both RANS and PANS models correctly predict the axial and circumferential velocities for case A. PANS, however, seems to capture more details of the circumferential component at the center of the draft tube. For case B (Fig. 7), the RANS k-ε model overpredicts the axial and circumferential velocities by as much as 19% and 35%, respectively, while PANS predictions show excellent agreement with data. The deficiency of the RANS model, as also shown by Vu et al. [12], appears mainly for the partial load conditions. As illustrated in Fig. 8(a), RANS underpredicts the axial velocity by as much as 64% near the centerline. Predictions can be improved by about 47% using the PANS model, showing quite good agreement with experimental data. Predictions of both models for the circumferential velocity are nearly the same, showing a relatively good agreement with data, namely 7% average deviation (see Fig. 8(b)). The deviation between RANS simulation results and data becomes larger as turbine's discharge decreases (farther away from the BEP). As shown in Figs. 9(a) and 9(c), the RANS k-ε model completely fails in correctly predicting the level of the axial velocity near the centerline for case D. It predicts a large backflow region, and as a result, overpredicts the level of the velocity outside of the shear flow region, while no reverse flow is observed within experimental data. This class of models has been developed for and calibrated by use of data from simple, steady flows. Therefore, they are insufficient for predicting strong, unstable shear layers as in the case of the partial load draft tube flow where a strong vortex is present. Predictions of the present PANS simulations show considerable improvements. The level of the axial velocity is well calculated in the shear layer, while better agreement is seen with experimental data near the centerline, although there still exists a notable underprediction. Based on these results, it is concluded that RANS models cannot correctly predict the flow behavior at partial load, where the low velocity inner region interacts with outer flow and the vortex rope forms. The precessing vortex rope enhances the mixing and turbulence production and diffusion that cannot be modeled using the RANS turbulence models. Using the PANS turbulence model, the level of mixing and flow unsteadiness is better predicted and results are improved.

The wall static pressure is shown at four sections (Sections S2–S5 in Figs. 4 and 10) in Fig. 11 for three operating conditions, i.e., cases A–C. In the experimental measurements [23], several pressure transducers were located peripherally at each section. Figure 10 shows the top view of these four sections together with the pressure monitoring points. The “inner” and “outer” sides of the draft tube's elbow are specified to help understanding the directions in Fig. 10. As seen in Fig. 11, wall pressure is lower near the inner side of the elbow (e.g., P5–P8 at section S3) and is higher near the outer side (e.g., P1, P2, P11, and P12 at section S3) for all operating conditions. This is expected due to the flow streamline curvature forced by the elbow. Results obtained from the PANS simulations are compared to the experimental data and the k-ε RANS results of Mauri [7]. The global trend of the experimental data is predicted fairly well by both models. Locally, however, the difference between RANS predictions and experimental data reaches 90% of the measured value. The RANS calculations clearly overestimate the influence of elbow at the section S2, but the differences decrease in the following sections. Predictions are significantly improved by using the PANS model, specifically at the outer side of the elbow for section S2 (points P1, P2, P11, and P12), and the inner side of the elbow for sections S3 and S4 (points P4–P8).

In order to compare simulation results and experimental data more quantitatively, deviations between RANS and PANS predictions and the experimental data are calculated point by point and are summarized in Table 5. For each case of A–C and each section of S2–S5, the maximum and the average percentage difference between calculated results and experimental data is shown in Table 5. In case A, the difference between the k-ε RANS predictions and data can reach 49% (Section S4) with the average difference for four sections being 15%. The average difference between PANS predictions and data is 9%, while improvements are seen for all sections. Similar results are seen for cases B and C with 7% and 8% improvements of predictions on average using the PANS model. Locally, a difference of 92% is seen between RANS results and data for case C while the maximum difference of PANS results and data never exceeds 30%.

Using the present unsteady PANS simulations, transient features of the draft tube flow can be investigated. Figure 12 shows contours of instantaneous (for an arbitrary instance in time) and time-averaged axial velocity in the draft tube for all operating conditions. Contours are plotted in the symmetry plane, showing the draft tube cone and elbow up to the pier. The evolution of the axial velocity from over-load to partial load can be clearly seen in Fig. 12. In case A, corresponding to 110% of the BEP flow rate, axial velocity is quite high especially at the center of the draft tube. By decreasing the flow rate towards the cases B and C, a stagnant region starts to grow at the center of the draft tube (the dark region in Fig. 12). This stagnant region is the result of the wake of the crown cone as well as the swirling nature of the flow which tends to decrease the momentum near the center and increase it near the wall. In case B, corresponding to 99% of the BEP flow rate, the stagnant region is small and is limited to a very narrow region close to the centerline of the draft tube. However, as the flow rate further decreases, the stagnant region grows in size and expands in the draft tube. The shear layer forming between this stagnant region and highly swirling outer flow becomes unstable and rolls up resulting in formation of a helical vortex which is known as the vortex rope. The vortex rope rotates around the vortex core as well as it precesses around the draft tube centerline axis. As shown for case D, by further decreasing of the flow rate, a very strong precessing vortex rope forms in the draft tube. Contours of the time-averaged axial velocity in Fig. 12 show the mean extent of the stagnant region in the draft tube, increasing by decreasing the flow rate.

As shown in Fig. 13, the present PANS model captures a strong precessing vortex rope in draft tube for the case D. The helical vortex is visualized in Fig. 12 at three arbitrary instants by the isopressure surfaces corresponding to P = −16,000 Pa. The vortex rope has a very unsteady nature, and its shape can dramatically change over time. The tail of the vortex rope may impact the inner side of the elbow wall as seen in Fig. 13(a). This impact, called the “shock phenomenon,” induces strong acoustic noise, pressure fluctuations, and even structural vibrations [30]. The unsteady pressure forces arising from this will be discussed in Sec. 5.4.

To demonstrate fluctuations in the draft tube due to the formation of a strong vortex rope in case D, wall pressure is monitored on the draft tube wall for this case. Figure 14 shows the evolution of wall pressure with time for eight points in the draft tube. The top row in Fig. 14 labeled with “inner path” shows pressure fluctuations for four points on sections S1, S2, S3, and S5 located in the inner side of the draft tube elbow. The bottom row labeled with “outer path” shows pressure fluctuations at the same sections but located on the outer side of the elbow. These eight points are shown in Fig. 4 with four inner points marked with red squares and four outerpoints marked with green circles. Results are shown for about 12 s after letting the simulations to run for an initial 4 s to make sure that a periodic unsteady (quasi-steady) state is reached. It is seen that pressure fluctuations show similar behavior for the inner and the outer points at section S1, which is due to the symmetry of the draft tube cone. It should be noted that, as expected, pressure fluctuations for these two points show a 180 deg phase lag. However, this 180 deg phase difference is not perfect for the all instants due to the effect of the elbow, which can be felt even in the cone, as well as the vortex break-apart causing additional fluctuations. The peak-to-peak amplitude of pressure fluctuations can reach 4,000 Pa in this section. The asynchronous (with a 180 deg phase lag) nature of these fluctuations exerts severe oscillatory forces on the draft tube wall. The long-term operation under this condition may results in structural damages to the draft tube components [2] as shown in Fig. 15. In cases where the draft tube is hung below the turbine and is not integrated with the powerhouse structure, these forces may cause failure of the lateral supports that brace the draft tube to the powerhouse [31]. Moving further downstream toward sections S2 and S3, pressure fluctuations show a very different behavior for the inner and outer points. Specifically, the inner side of the elbow demonstrates very severe pressure fluctuations with the peak-to-peak amplitude as high as 14,000 Pa. This is mainly due to the interactions between the vortex rope, the secondary flows, and the low pressure region in the inner side of the elbow. This is the region where vortex rope may impact the draft tube wall causing the “shock phenomenon” [30] as shown in Fig. 13(a). The concrete erosion in the draft tube elbow [31] can be related to the strong forces exerted by these phenomena.

Figure 16 shows the root-mean-square (RMS) of pressure (P_rms) for the eight monitored points in Fig. 14. It is seen that P_rms in the inner side of the elbow increases to as much as 2.4 times the one in the cone (from 1,397 Pa to 3,343 Pa) and then decreases further downstream to 2,800 Pa, while the outer path shows a relatively consistent behavior in terms of changing P_rms (around 1,500 Pa). Therefore, it is concluded that for a draft tube of a hydroturbine operating under partial load, where a strong vortex rope forms in the draft tube, the most critical region where severe pressure fluctuations are felt is the inner side of the draft tube elbow.

Unsteady numerical simulations of flow in the draft tube of a model Francis turbine is carried out in this study. The test case is the FLINDT draft tube reconstructed in this study compiling data from several previously published papers. To the best of authors' knowledge, this is the first attempt on regenerating the FLINDT draft tube and making it available in open literature. A newly developed PANS model is used as the turbulence closure model. Four operating conditions ranging from 110% to 70% of the BEP flow rate is considered and several parameters including the pressure recovery coefficient, mean velocity, and wall pressure obtained from the present PANS simulations are compared with those from the experimental measurements as well as those obtained by RANS k-ε simulations. It is shown that RANS and PANS both can predict the flow behavior close to the BEP operating condition. However, RANS results deviate considerably from the experimental data as the operating condition moves away from the BEP. The pressure recovery factor predicted by the k-ε model shows more than 13% overprediction for case C with 91% of BEP flow rate, while axial velocity is underpredicted by more than 64% near the centerline. In case D, with even lower flow rate (70% of BEP), the k-ε model fails substantially in predicting mean axial velocity while more than 58% difference is seen in prediction of the recovery factor. Predictions can be improved dramatically using the present unsteady PANS simulations. Specifically, the pressure recovery factor is predicted to less than 4% deviation and mean axial velocity is calculated to less than 17% difference with experimental data for case C. Even for the case D, with a strong vortex rope forming in the draft tube and flow being very unstable, the pressure recovery factor is predicted to only 6% difference with data, while considerable improvement is seen in compassion to the RANS results, although still notable underprediction is seen near the center of the draft tube. Furthermore, details of the unsteady flow field, including the unstable shear layer and the stagnant region, and the precessing helical vortex rope are captured using the unsteady PANS simulations. It is shown that the formation of the vortex rope in partial load conditions results in severe pressure fluctuations with the peak-to-peak amplitude as high as 14,000 Pa. The oscillatory forces exerted on the draft tube by these fluctuations may result in structural damages as have been observed in real–life. Investigations of the wall unsteady pressure reveals that the critical region in the draft tube with highest pressure fluctuations is located at the inner side of the elbow.

The work presented herein was funded in part by the Office of Energy Efficiency and Renewable Energy (EERE), U.S. Department of Energy, under Award No. DE-EE0002667 (the DOE/PSU Graduate Student Fellowship Program for Hydropower Research) and DE-EE0002668 (The HRF Fellowship), and the Hydro Research Foundation.

f_k =

unresolved-to-total ratio of turbulent kinetic energy

f_ε =

unresolved-to-total ratio of turbulence dissipation rate

H =

turbine head

k =

turbulent kinetic energy

P =

pressure

Q =

volume flow rate

r =

radial coordinate

R =

runner radius

$u¯i$ =

partially averaged velocity

x_i, x, y, z =

coordinates

$y+$ =

dimensionless wall distance

Greek Symbols

$Δ$ =

grid length scale

$ε$ =

dissipation rate of turbulent kinetic energy

$νu$ =

pans eddy viscosity

$τij$ =

subfilter stress tensor

$ϕ$ =

flow rate (discharge) coefficient (= Q/($πωR3$))

$ψ$ =

head (energy) coefficient (= 2gH/($ω2R2$))

ω =

angular velocity of the runner

$Λ$ =

turbulence length scale

$χ$ =

pressure recovery factor

Abbreviations

BEP =

best efficiency point

LES =

large eddy simulation

PANS =

partially averaged Navier–Stokes

RANS =

Reynolds-Averaged Navier–Stokes

URANS =

unsteady RANS