Investigation of Flow-Induced Instabilities in a Francis Turbine Operating in Non-Cavitating and Cavitating Part-Load Conditions

Abstract

The integration of intermittent renewable energy resources to the grid system requires that hydro turbines regularly operate at part-load conditions. Reliable operation of hydro turbines at these conditions is typically limited by the formation of a Rotating Vortex Rope (RVR) in the draft tube. In this paper, we investigate the formation of this vortex using the scale-resolving methods SST-SAS, wall-modeled LES (WMLES), and zonal WMLES. The numerical results are first validated against the available experimental data, and then analyzed to explain the effect of using different scale-resolving methods in detail. It is revealed that although all methods can capture the main features of the RVRs, the WMLES method provides the best quantitative agreement between the simulation results and experiment. Furthermore, cavitating simulations are performed using WMLES method to study the effect of cavitation on the flow in the turbine. These effects of cavitation are shown to be highly dependent on the amount of vapor in the RVR. If the amount of vapor is small, cavitation induces broadband high-frequency fluctuations in the pressure and forces exerted on the turbine. As the amount of cavitation increases, these fluctuations tend to have a distinct dominant frequency which is different from the frequency of the RVR.

1. Introduction

Due to their large storage capacity and quick response, hydro power plants are more and more being used to facilitate the integration of intermittent energy produced by other renewable sources of energy such as wind and solar to the grid system. This requires that hydro power plants expand their range of operation and also operate frequently at off-design conditions. In such off-design operation, the flow in different components of a hydro power plant is prone to complex behaviors and instabilities. One example is the formation of a Rotating Vortex Rope (RVR) in the draft tube which is associated with a high level of pressure fluctuations in the system. These pressure fluctuations can cause oscillation in the mechanical loading of the runner, causing a significant swing in the power produced by the hydro power plant [1,2]. Furthermore, the pressure drop in the core of the RVR can lead to the formation of vapor pockets in the water flow, known as cavitation. The presence of cavitation in hydro turbines is known to have several undesirable effects, such as cavitation erosion [3], performance degradation [4], and increased level of pressure fluctuations and vibration [5].

The detrimental effects of the RVR for hydro turbines operating at off-design conditions have motivated several experimental studies with the aim to investigate the formation of the RVR and its associated cavitating structures. Nishi et al. [6] experimentally examined non-cavitating and cavitating RVRs at different flow conditions. They observed that the pressure fluctuations induced by the RVR can be decomposed into synchronous and convective components. Iliescu et al. [7] studied the dynamics of non-cavitating and cavitating RVRs using High-Speed Visualization (HSV) and Particle Image Velocimetry (PIV) measurements. They concluded that the formation of vapor pockets in an RVR can lead to fluctuations in its geometrical features. Favrel et al. [8] and Arpe et al. [9] studied the effect of cavitation on the pressure fluctuations induced by the RVR. They both concluded that cavitation formation in the RVR can influence both the amplitude and the frequency of the pressure fluctuations. Further, Landry et al. [10] showed that the presence of cavitation in the RVR can alter the eigenfrequency of the hydraulic system by changing the sound propagation speed in the liquid. This can lead to severe vibrations and power swing if there is a match between this eigenfrequency and the frequency of induced pressure fluctuations by the RVR [1,2].

In addition to experimental investigations, CFD simulations have been used to investigate this flow. In the majority of these simulations, Reynolds-averaged Navier–Stokes (RANS) approaches are used as they have lower computational cost compared to scale-resolving methods. These RANS simulations have been shown to have deficiencies in capturing the correct behavior of RVRs. Ciocan et al. [11] studied the dynamics of the RVR using unsteady RANS and showed that their results under-predicted the level of swirl in the flow entering the draft tube and differed by 13% in RVR frequency compared to the experimental data. Similar deficiency has been observed by Liu et al. [12] who performed unsteady RANS simulations of the RVR in a draft tube and compared the level of pressure fluctuations with the experimental data. Furthermore, the simulations by Ruprecht et al. [13] have shown that using a RANS approach over-estimates the decay of the swirling flow around the RVR, which results in the under-prediction of the pressure fluctuations caused by the RVR. Despite these deficiencies, a few studies have used the RANS approach to study the effect of cavitation on the behavior of RVRs. Yu et al. [14] performed RANS simulations of a cavitating RVR and revealed that the presence of cavitation affects the vorticity production around the RVR as well as the pressure fluctuations generated in the draft tube. Jošt et al. [15] performed simulations of cavitating and non-cavitating flows in a Francis turbine and compared the results with experimental data. They showed that it is necessary to consider cavitation in the simulation to correctly reproduce the RVR-induced pressure fluctuations in the experiment.

To avoid the limitation of RANS approaches in capturing the correct behavior of RVRs, several numerical studies have used scale-resolving approaches to examine non-cavitating and cavitating RVRs. Salehi et al. [16] and Salehi and Nilsson [17] investigated the formation of RVRs during the startup and shutdown in a Francis turbine using a scale-resolving approach. By analyzing the results, they were able to explain the complex flow behavior leading to the formation of RVRs in the draft tube. Foroutan and Yavuzkurt [18] studied the RVR using DES and found a good agreement between their numerical results and experimental data. By analyzing the DES results, they postulated that the RVR is created in the shear layer between the region of low axial velocity in the center of the draft tube cone and the flow outside of this region, due to Kelvin–Helmholtz instability. Minakov et al. [19] used DES to study the RVR in the draft tube at different guide vane openings, and found that changing the guide vane opening leads to a different swirl number and discharge coefficient, which leads to different behavior of the RVR. Rajan and Cimbala [20] studied the flow in a simplified draft tube using the DES approach for different discharge coefficients. It was shown that the level of pressure fluctuations increases substantially when the RVR is formed in the draft tube. Guo et al. [21] used the LES approach to simulate cavitating flows in the draft tube. They highlighted the importance of including the runner to predict the low-pressure region in the center of the RVR. Pacot et al. [22] performed LES simulations for cavitating RVRs and showed that a large amount of cavitation in the RVR can prevent its precession motion, leading to a reduction of the RVR-induced pressure fluctuations.

Although the numerical studies reviewed above have highlighted the importance of using scale-resolving methods to capture the correct behavior of the RVR, only a few of them have performed a detailed comparison between the performance of different scale-resolving approaches. Furthermore, detailed analyses of how cavitation affects the flow features of the RVR are also scarce in the literature. To address this knowledge gap, in the present paper, the formation of non-cavitating and cavitating RVRs in the Francis-99 turbine is investigated using three scale-resolving methods. The simulation results are compared with the experimental data provided through the Francis-99 workshop, and the differences between the results obtained by different scale-resolving methods are highlighted. Then, cavitating simulations are performed for three cavitating conditions, covering the condition near the inception of cavitation in the RVR to the fully cavitating RVR. Using these simulations, the effect of cavitation on the behavior of the RVR, the generated pressure fluctuations and the forces exerted on different component of the turbine are discussed in detail. This paper is organized into five sections. After this introduction, the next section briefly summarizes the employed numerical set-up. Then, the studied turbine and flow conditions are described in the following section. The results are presented in the fourth section, where first the effects of different scale-resolving methods on capturing the flow features are discussed, and then, the effects of cavitation on the flow features of the RVR are explained in detail.

2. Numerical Set-Up

In this paper, the modified interPhaseChangeFoam solver [23,24,25] from the OpenFOAM-2.2.x framework [26] is used to perform the simulations. The governing equations in this solver are the incompressible Navier–Stokes equations for two-phase (liquid-vapor) isothermal flows. Using the homogeneous mixture approach, only one set of equations is solved for the two-phase mixture. Similar to turbulent single-phase flows, turbulent two-phase mixture flows include a wide range of scales. For engineering flow, RANS, LES, or zonal hybrid RANS/LES approaches can be employed, where the governing equations of a two-phase mixture with homogeneous assumption read

$$\frac{\partial }{\partial t}\left(\rho \right)+\nabla ·\left(\rho \tilde{\mathbf{u}}\right)=0,$$

(1)

$$\frac{\partial }{\partial t}\left(\rho \tilde{\mathbf{u}}\right)+\nabla ·(\rho \tilde{\mathbf{u}}\otimes \tilde{\mathbf{u}})+\nabla ·\left([p\mathbf{I}-\tau ]\right)+\nabla ·\left({\tau }_{unres}\right)=0,$$

(2)

where $\rho$ is the mixture density, $\mathbf{u}$ is the velocity vector, p is the pressure, $\mathbf{I}$ is the identity tensor, $\tau$ is the viscous stress tensor, and ${\tau }_{unres}$ is the stress tensor due to unresolved turbulence in the flow. The tilde operation on the velocity vector in the above equations, $\tilde{\mathbf{u}}$, represents a time-averaging operation in the RANS approach or spatial filtering in the LES approach. Assuming that the mixture of liquid and vapor is homogeneous and that the dynamic viscosity and the density in each phase are constant, the mixture viscous stress tensor, $\tau$, and the mixture density $\rho$ can be obtained from

$$\begin{array}{c}\hfill \tau =({\alpha }^l{\mu }^l+{\alpha }^v{\mu }^v)S,\end{array}$$

(3)

$$\begin{array}{c}\hfill \rho ={\alpha }^l{\rho }^l+{\alpha }^v{\rho }^v,\end{array}$$

(4)

where S is the mixture strain tensor, and ${\alpha }^l$, ${\alpha }^v$, ${\rho }^l$, and ${\rho }^v$ are, respectively, the liquid volume fraction, the vapor volume fraction, the density of the liquid and the density of the vapor. In order to close the governing equations, ${\alpha }^l$, ${\alpha }^v$ and ${\tau }_{unres}$ should be determined. In non-cavitating simulations, the liquid volume fraction, ${\alpha }^l$, is equal to 1 while the volume fraction of the vapor phase, ${\alpha }^v$, is equal to zero in the entire domain. For two-phase cavitating flows, ${\alpha }^l$ is obtained using a transport equation, reading

$$\frac{\partial }{\partial t}\left({\alpha }^l{\rho }^l\right)+\nabla ·\left({\alpha }^l{\rho }^l\tilde{\mathbf{u}}\right)=\dot{m},$$

(5)

where $\dot{m}$ is the mass transfer term due to vaporization and condensation. The vapor volume fraction, ${\alpha }^v$, is then determined using ${\alpha }^v+{\alpha }^l=1$. To model the mass transfer term, $\dot{m}$, the Schnerr–Sauer model [27] is used in this paper. This term can written as

$$\dot{m}={\alpha }^l\left({\dot{m}}_{{\alpha }_v^l}-{\dot{m}}_{{\alpha }_c^l}\right)+{\dot{m}}_{{\alpha }_c^l},$$

(6)

where ${\dot{m}}_{{\alpha }_c^l}$ and ${\dot{m}}_{{\alpha }_v^l}$ are, respectively, condensation and vaporization terms and can be defined as

$${\dot{m}}_{{\alpha }_c^l}=C_c{\alpha }^l\frac{3{\rho }^l{\rho }^v}{\rho R_B}\sqrt{\frac{2}{3{\rho }^l}}\sqrt{\frac{1}{\left|p-p_v\right|}}\max \left(p-p_v,0\right),$$

(7)

$${\dot{m}}_{{\alpha }_v^l}=C_v\left(1+{\alpha }_{Nuc}-{\alpha }^l\right)\frac{3{\rho }^l{\rho }^v}{\rho R_B}\sqrt{\frac{2}{3{\rho }^l}}\sqrt{\frac{1}{\left|p-p_v\right|}}\min \left(p-p_v,0\right).$$

(8)

In the above equations, $C_c$ and $C_v$, the condensation and vaporization constants, are 1, $p_v$ is the vapor pressure, ${\alpha }_{Nuc}$ is the initial volume fraction of nuclei, and $R_B$ is the radius of the nuclei which is obtained from

$$R_B=\sqrt[3]{\frac{3}{4\pi n_0}\frac{1+{\alpha }_{Nuc}-{\alpha }^l}{{\alpha }^l}}.$$

(9)

The initial volume fraction of nuclei is calculated from

$${\alpha }_{Nuc}=\frac{\frac{\pi n_0{d}_{Nuc}^3}{6}}{1+\frac{\pi n_0{d}_{Nuc}^3}{6}},$$

(10)

where the average number of nuclei per cubic meter of liquid volume, $n_0$, and the initial nuclei diameter, $d_{Nuc}$, are assumed to be ${10}^{11}$ and ${10}^{-5}$, respectively.

To model the ${\tau }_{unres}$ term, the eddy-viscosity hypothesis is used as

$${\tau }_{unres}=-2{\mu }_{unres}\overline{S}+\frac{2}{3}{k}_{unres}\mathbf{I},$$

(11)

where ${\mu }_{unres}$ is the eddy viscosity and $k_{unres}$ is the unresolved turbulent kinetic energy. In this paper, we use three models to approximate ${\mu }_{unres}$ and $k_{unres}$ which are explained in the following.

2.1. SST-SAS RANS Approach

The eddy viscosity, ${\mu }_{unres}$, and unresolved turbulent kinetic energy, $k_{unres}$, represent, respectively, the turbulent viscosity, ${\mu }_t$, and the turbulent kinetic energy, $k_t$ within the context of the RANS approach. Here, these terms are modeled using the Shear Stress Transport-based Scale-Adaptive Simulation model (SST-SAS) [28,29]. In this turbulence model, the turbulent kinetic energy, $k_t$, and the turbulent specific rate of dissipation, ${\omega }_t$, are calculated using transport equations. These two terms are then used to calculate ${\mu }_t$ as

$${\mu }_t=\rho a_1\frac{k_t}{\max (a_1{\omega }_t,b_1{F}_{23}S)},$$

(12)

where S is the invariant measure of the strain rate, $a_1$ and $b_1$ are constants, and $F_{23}$ is a blending function. The transport equation for $k_t$ in $\mathrm{SST}-\mathrm{SAS}$ is the same as the equation in the $\mathrm{k}-\omega \mathrm{SST}$ turbulence model [30], and can be written as

$$\frac{\partial k_t}{\partial t}+\nabla ·\left(\tilde{u}{k}_t\right)=P_k+{\nabla }^2\left[\frac{1}{\rho }(\mu +{\mu }_t)k_t\right]-{\beta }^{*}{\omega }_t{k}_t,$$

(13)

where $P_k$ and ${\beta }^{*}$ are, respectively, the production term and a constant. The transport equation for ${\omega }_t$ can be read as

$$\begin{array}{c}\frac{\partial {\omega }_t}{\partial t}+\nabla ·\left(\tilde{u}{\omega }_t\right)=\alpha \frac{P_k}{{\mu }_t}-{\nabla }^2\left[\frac{1}{\rho }(\mu +{\mu }_t){\omega }_t\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill -\beta {\omega }_t^2-2(F_1-1)\rho {\sigma }_{w2}\frac{1}{\omega }\nabla k_t\nabla {\omega }_t+Q_{SAS},\end{array}$$

(14)

where $\alpha$, $\beta$, and ${\sigma }_{w2}$ are constants, $F_1$ is a blending function, and $Q_{SAS}$ is an extra term compared to the equation in the $\mathrm{k}-\omega \mathrm{SST}$ turbulence model. This extra term can be obtained as

$$Q_{SAS}=\max \left[\rho {\zeta }_2\kappa S^2{\left(\frac{L}{L_{vK}}\right)}^2-C\frac{2\rho k_t}{{\sigma }_{\mathsf{\Phi }}}\max \left(\frac{{\left|\nabla {\omega }_t\right|}^2}{{\omega }_t^2},\frac{{\left|\nabla k_t\right|}^2}{k_t^2}\right),0\right],$$

(15)

where ${\zeta }_2$, $\kappa$, ${\sigma }_{\mathsf{\Phi }}$, and C are constants. In Equation (15), L and $L_{vK}$ are the modeled turbulence length scale and the von Karman length scale, respectively, which are defined as

$$L=\frac{\sqrt{k}}{c_{\mu }^{1/4}\omega },L_{vK}=\max \left(\frac{\kappa S}{\left|{\nabla }^2\tilde{\mathbf{u}}\right|},C_s\sqrt{\frac{\kappa {\zeta }_2}{\beta /c_{\mu }-\alpha }}\right),$$

(16)

where $c_{\mu }$, $\beta$, $C_s$ and $\alpha$ are constants.

2.2. WALE LES Approach

Within the LES context, ${\mu }_{unres}$ and $k_{unres}$ are the sub-grid scale viscosity, ${\mu }_{sgs}$, and the sub-grid turbulent kinetic energy, $k_{sgs}$, respectively. Here, the wall-adapting local eddy-viscosity (WALE) LES model proposed by Nicoud and Ducros [31] is used to model these terms as

$${\mu }_{sgs}=C_k\Delta \sqrt{k_{sgs}},$$

(17)

$$k_{sgs}={\left(\frac{C_w^2\Delta }{C_k}\right)}^2\frac{{\left({\tilde{S}}^d{\tilde{S}}^d\right)}^3}{{\left({\left(\tilde{S}\tilde{S}\right)}^{5/2}+{\left({\tilde{S}}^d{\tilde{S}}^d\right)}^{5/4}\right)}^2},$$

(18)

where $\Delta$ is the cell length scale calculated based on the cubic root of the cell volume and $C_k=1.6$ and $C_w=0.325$ are the model constants. $\tilde{S}$ and ${\tilde{S}}^d$ are also, respectively, the resolved-scale strain rate tensor and traceless symmetric part of the square of the velocity gradient tensor.

2.3. Zonal RANS/LES Approach

In the zonal approach, the $\mathrm{SST}-\mathrm{SAS}$ RANS and WALE LES approaches are used in predefined regions of the domain. To mark these predefined regions, a scalar field $r_{LES}$ is set to 1 in the LES regions and zero in the RANS regions. Using this scalar field, the ${\mu }_{unres}$ and $k_{unres}$ can be obtained as

$${\mu }_{unres}=r_{LES}{\mu }_{sgs}+(1-r_{LES}){\mu }_t,$$

(19)

$$k_{unres}=r_{LES}{k}_{sgs}+(1-r_{LES})k_t,$$

(20)

where ${\mu }_{sgs}$ and $k_{sgs}$ are obtained using the WALE LES approach from Equations (17) and (18), while ${\mu }_t$ and $k_t$ are determined using the $\mathrm{SST}-\mathrm{SAS}$ RANS approach from Equations (12) and (13).

2.4. Discretization Schemes and Solution Algorithm

To discretize convective terms in the momentum equations, the Linear-Upwind Stabilized Transport (LUST) convection scheme is used. This scheme blends 2nd order upwind (25%) and central differencing (75%) schemes [32]. The discretization of diffusion terms in the momentum equations are done using the linear scheme. For the convective term in the liquid volume fraction transport equation, the first-order upwind scheme is used. Using this scheme is consistent with the mixture assumption and is recommended in the original publication [33] for the Schnerr–Sauer model. A second-order backward implicit scheme is used for time discretization and a pressure-based PIMPLE approach, a combination of SIMPLE and PISO algorithms, is employed to solve the discretized equations. To ensure the convergence of the algorithm, the residual target for the pressure and velocities are set to ${10}^{-13}$ for the simulations with cavitation and ${10}^{-7}$ for the simulations without cavitation. For more details about the solution procedure, the reader can refer to Asnaghi et al. [24] and Bensow and Bark [23].

3. The Francis-99 Turbine and Computational Mesh

In this paper, the Francis-99 turbine model is used. This turbine is a scaled-down model of the prototype Francis turbine installed at the Tokke power plant in Norway [34]. The computational domain shown in Figure 1 includes the spiral casing, the stay vanes, the guide vanes, the runner, and the draft tube. There are 14 stay vanes and 28 guide vanes. The runner includes 15 splitter blades and 15 full-length blades. The inlet and outlet diameters of the runner are 0.63 m and 0.347 m, respectively. The mesh in different components of the turbine are produced using the Pointwise V18.3 mesh generation software. The mesh specifications for each component are shown in

Table 1. Mesh resolution.

Components	# of Cells	Average ${\mathit{y}}^{+}$
SC+SV+GV	10.12 M	20
RU	5.13 M	12
DT	5.63 M	14
Total	20.88 M	17

. The mesh has 20.88 M cells, which according to the mesh dependency studies performed for simulations of the Francis-99 turbine found in the literature [35,36,37,38,39] is enough for mesh-independent results of the flow in the draft tube. The table also shows that the average $y^{+}$ value is larger than 1 for all components; therefore, the wall function based on Spalding’s law is used at the walls [25,40,41]. It should be mentioned that using the wall function for the simulations of the flow in the draft tube of Francis turbines is justified as the work by Wilhelm et al. [42] showed that resolving near-wall regions instead of using the wall function has insignificant effects on the captured flow in the draft tube.

As mentioned in Section 2, we use a zonal WMLES approach where the RANS and WMLES approaches are used in different regions of the domain. Since the focus of the paper is to investigate the RVR in the draft tube, we use the WMLES approach in the draft tube (region labeled by DT in Figure 1) and in the rest of the domain (regions labeled by RU, SC, SV, and GV in Figure 1), the $\mathrm{SST}-\mathrm{SAS}$ RANS approach is employed.

The experimental data used for validation in the present study are from the second Francis-99 workshop. These data include both pressure and velocity measurements at different locations. The velocity measurements includes the axial and horizontal components over three lines in the draft tube, Line 1, Line 2, and Line 3 in Figure 2a. The data also include the static pressure fluctuations at two probes in the draft tube, Probe 2 and Probe 3 in Figure 2a and the static pressure in the vaneless region between the guide vane blades and runner blades, is marked by Probe 1 in Figure 2b.

3.1. Studied Flow Conditions

In this paper, the part-load (PL) condition in the second Francis-99 workshop is studied. In this condition, the guide vane opening is $\alpha ={6.72}^{\circ }$, the runner angular speed is $n=332.84\mathrm{r}/\min$, and the discharge is $Q=0.13962{\mathrm{m}}^3/\mathrm{s}$. In order to study the effect of cavitation in this PL condition, the simulations are performed for both non-cavitating and cavitating conditions. The cavitation number, $\sigma$, in these simulations is defined as

$$\sigma =\frac{(p_2-p_v)/{\rho }_l+\frac{1}{2}{(Q/A_2)}^2}{gH},$$

(21)

where $p_2$ and $A_2$ are, respectively, the pressure and the cross-section area at the draft tube outlet, $p_v$ is the vapour pressure, ${\rho }_l$ is liquid density, and H is the turbine head. Here, we study the flows at $\sigma =\left\{0.07,0.06,0.05\right\}$ which includes the near-inception cavitation in the RVR ($\sigma =0.07$) as well as the condition corresponding to the fully cavitating RVR ($\sigma =0.05$).

4. Results

The results are divided into two parts. The first part presents the effect of turbulence modeling on capturing the global quantities, the velocity profiles and the pressure fluctuations in the draft tube. The second part is devoted to studying the effects of cavitation on the same flow features.

4.1. Effect of Turbulence Modeling

In

Table 2. Effect of turbulence modeling technique on predicted global quantities, torque, $M_z$, head, H, and efficiency $\eta$.

Quantities	SAS	Zonal WMLES	WMLES	Exp.
$M_z$ (Nm)	440.57	440.51	439.68	420.79
Relative error for $M_z$	4.7%	4.7%	4.5%	-
H (m)	12.71	12.70	12.48	11.87
Relative error for H	7.1%	7.0%	5.1%	-
$\eta$	88.14	88.19	89.58	90.13
Relative error for $\eta$	2.2%	2.2%	0.6%	-

, the global quantities captured with different turbulence modeling techniques are compared with the experimental values. The relative errors in this table are obtained by dividing the difference between the experimental and numerical values by the experimental values. This comparison shows that all of the turbulence modeling techniques can capture these global quantities with a relatively small error, although a slightly lower relative error can be seen in the WMLES results. According to Čelič and Ondráčka [43], this error can be due to neglecting the labyrinth seal and disk friction losses.

Figure 3 compares the time-averaged experimental and numerical velocity profiles along the three lines shown in Figure 2a. The experimental profiles are taken from the data provided by the second Francis-99 workshop [44]. The experimental axial velocity over Line 1 and Line 2 indicates that a region with low values of absolute axial velocity exists near the center of the draft tube. The comparison between the numerical results in Figure 3a,c indicates that this region is captured by all simulations, although there are some quantitative differences between the different results. As it can be seen in the axial velocity profiles on Lines 1 and 2 (Figure 3a), the regions with low values of absolute axial velocity in the SAS and zonal WMLES results are more confined to the center of the draft tube cone compared to the WMLES results and the experimental data. The axial velocity along the centerline (Figure 3c) shows that while the WMLES simulation predicts a negative averaged axial velocity along the entire Line 3, similar to the experimental data, positive values for the time-averaged axial velocity can be seen in the SAS and zonal WMLES results. Considering the definition of the axial direction (z), shown in Figure 2a, this positive axial velocity means that there is a reversed flow along the centerline of the draft tube cone in these two simulations. The cause of this difference is explained later in this paper. The horizontal velocity profiles (Figure 3b,c) show that the simulations can capture the trends similar to the experiment. However, the values of horizontal velocity in these simulations are shifted compared with the experimental values. It should be noted that the absolute value of the horizontal velocity is very close to zero over the measurement lines, the relative uncertainty for these velocities is higher compared with that for the axial velocities according to Salehi and Nilsson [45]. This higher uncertainty can be one reason for the difference between the results for these velocity profiles.

It is well-established that the region with low values of absolute axial velocity shown in Figure 3 is caused by the formation of the RVR in the draft tube at part-load condition. Figure 4 shows snapshots of the RVRs of the different numerical results. In this figure, the RVRs are visualized using an iso-surface of Q, which is the second invariant of the velocity gradient tensor. The level of the iso-surface is $5000{\mathrm{s}}^{-2}$, which is chosen for an optimal visualization of the RVRs. In the zoom-in view, it can be seen that the RVR consists of many vortices wrapping around each other. At the location of these vorticies, the positive value of axial velocity (red color in the figure) indicates the presence of a reversed flow. The figure also presents the time history of the axial velocity on Line 1 for one period of the RVR rotation. The time instance corresponding to the snapshot of the iso-surfaces is shown by black dashed lines. These time histories show that as the RVR passes Line 1, the axial velocity on these lines becomes positive (red regions). The sweeping motion of the vortices causes the tilted streaks seen in the time history. A comparison between the results from the different simulations indicates that the number and size of the vortices in the RVR are strongly influenced by the selection of the turbulence modeling technique. In the SAS simulation, the RVR consists of a small number of large vortices, as shown in the zoom-in views, and this leads to fewer and large streaks in the time history of the axial velocity. In the zonal WMLES and WMLES simulations, however, the RVR has a large number of smaller vortices, which creates thinner streaks in the time history.

As mentioned earlier and shown in Figure 3, the region with a low value of absolute axial velocity in the center of the draft tube cone is affected by the selection of the turbulence modeling. To provide a reason for this effect, the phase-averaged axial velocity over Lines 1 and 2 in the experiment and simulations is shown in Figure 5. To obtain the phase-averaged values in this paper, the data signal corresponding to one cycle of the RVR rotation, $T_{RVR}$, is divided into 30 windows. The data corresponding to each window are averaged together. Two instances A and B are marked in these figures. At instance A, the axial velocity has a high positive value which is due to the passage of the RVR over the measurement lines. At instance B, the RVR has rotated further and left the measurement lines, which leads to the observed decrease in the axial velocity toward the negative values. The comparison in Figure 5 shows that when the vortex passes the measurement lines at instance A, all simulations predict a distribution of axial velocity which is quite similar to the experimental data. However, this is not case for the instance when the RVR leaves the measurement lines (instance B). At this instance, the axial velocity over Line 1 near the center of the draft tube ($y=0$ mm) has positive values in the SAS and zonal WMLES simulations while the values in the experiment and the WMLES simulations are negative. The reason for this difference is that the RVR in the SAS and zonal WMLES simulations rotates on a path which is closer to the center of draft tube as compared with the experiment and the WMLES simulation. To clearly show this, the figure shows vertical black lines passing through the maximum values of axial velocity when the RVR is on the measurement line. The distance between these lines is also shown. It can be seen that distance between these lines is larger in the experiment and the WMLES results as compared with the other two simulations. This difference means that a portion of the RVR in the SAS and zonal WMLES simulations has overlaps with Line 1 near the center of the draft tube all through the rotational cycle of the RVR. Since the axial velocity in the RVR is positive according to Figure 4, this overlap would lead to the positive values of axial velocity over Line 1 at instance B in the SAS and zonal WMLES simulations. It also leads to the positive values of the averaged axial velocity in these two simulations which is shown in Figure 3.

The RVR in the draft tube creates a large amount of pressure fluctuations which can lead to vibrations. To evaluate how well these pressure fluctuations can be captured by the different turbulence modeling techniques, Figure 6 compares the experimental and numerical pressure fluctuation signals in Probe 2 in the draft tube. For each signal, the Root Mean Square (RMS) of the fluctuations is noted in the plot. It can be seen that in all of the numerical results, the RMS values are lower compared with the values in the experiment. The numerical RMS values however increase towards the experimental one as the resolution of the turbulence modeling technique increases. It can be seen that, as expected, it is mainly the smaller scales of the fluctuations that differ between the results from the different turbulence modeling techniques, while the amplitude and frequency of the large-scale RVR motion are similar.

To investigate the reasons why the RMS values of the pressure fluctuations are lower in the simulations compared with the experiment, Figure 7 shows a frequency analysis of the pressure fluctuations. For the experiment, the analysis (Figure 7a) includes the power spectrum analysis of the pressure signal from Probe 2 (left plot), the coherence (middle plot), and the phase difference (right plot) between the pressure signals from Probes 2 and 3. The coherence and phase difference between the two signals are calculated from the cross-spectral density, $P_{xy}$, which is obtained using Welch’s method [46]. The frequency analysis of the experimental pressure signal (left plot) indicates the existence of two dominant frequencies, $f_1$ and $f_2$. For these dominant frequencies, the coherence between the signals at Probes 2 and 3 (middle plot) is almost one, which means that these dominant frequencies also exist in the signal from Probe 3. The phase difference between the two signals (right plot) is $\pi$ for the dominant frequency $f_1$ and its harmonic $2f_1$. Considering that the locations of Probes 2 and 3 are exactly at opposite sides of the cone region in the draft tube, this phase difference suggests that the dominant frequency, $f_1$, is due to the precession of the RVR in the draft tube cone. It should be mentioned that the precession frequency $0.29f_{runner}$ is in the range 0.2–0.4 $f_{runner}$, which has been found in previous studies [7,47]. For the dominant frequency $f_2$, the phase difference is almost zero indicating that the corresponding pressure fluctuations are synchronous meaning that they have the same phase and amplitude for the pressure sensors located in the same cross section of the draft tube. The frequency analyses of the pressure signals in the simulations with different turbulence modeling techniques show that the dominant frequency of the RVR, $f_1$, and its harmonic, $2f_1$, is captured by all turbulence modeling techniques. The experimental dominant frequency $f_2$ can however not be seen in any of the numerical results, which causes a reduced RMS value.

To investigate the origin of the dominant frequency $f_2$ in the experiment, Figure 8 shows a spectral analysis of the pressure fluctuations in the draft tube at the Best Efficiency Point (BEP) and High Load (HL) conditions for which there is no RVR. Similar to the PL condition, a dominant frequency at $f_2$ can be seen also for BEP and HL (left plots). As for the PL condition, the coherence (middle plots) is close to one and the phase difference (right plots) is zero for this frequency. This indicates a synchronous nature of these pressure fluctuations, at a frequency that is rather independent of the operating condition. Based on this observation, we can conclude that the synchronous pressure fluctuations seen in the experimental results is related to a component in the system rather than the flow features in the components studied here. It should be mentioned that similar synchronous pressure fluctuations have been observed by Favrel et al. [8], Arpe et al. [9], with frequencies in the range of 2–4 $f_{runner}$. These studies, however, have shown that this type of pressure fluctuations occurs in cavitating conditions and can be attributed to the interaction between the cavitating RVR and the elbow in the draft tube.

To further analyze the high-frequency content of the signals and its effects on the RMS values, the fluctuations are decomposed into two components as

$$p^{\prime }=p_{RVR}^{\prime }+p_{Other}^{\prime },$$

(22)

where $p_{RVR}^{\prime }$ denotes the pressure fluctuations due to the precession of the RVR and $p_{Other}^{\prime }$ denotes the pressure fluctuations due to other sources. To perform this decomposition, we assume that the pressure fluctuations due to RVR at Probe 2 have a phase difference of $\pi$ with the corresponding fluctuations at Probe 3. This assumption is shown to be true in Figure 7a. Based on this assumption, the frequencies of pressure fluctuations due to RVR are determined. These frequencies are then filtered from the signal of pressure fluctuations to obtain the pressure fluctuations due to other sources, $p_{Other}^{\prime }$. The fluctuations due to RVR, $p_{RVR}^{\prime }$, then can be obtained by subtracting $p_{Other}^{\prime }$ from the original pressure fluctuation signal according to Equation (22).

Figure 9 shows the different components of the pressure fluctuations at Probe 2 according to Equation (22) and their corresponding RMS values from the results of the simulations and experiment. For the pressure fluctuations due to the RVR, the difference between the predicted RMS values and the experimental RMS value correspond to 16%, 14% and 4% for the SAS, zonal WMLES and WMLES simulations, respectively. The reason for these differences will be explained later. The corresponding differences for the pressure fluctuations due to other sources are 77%, 59%, and 39%. The reason for these large differences is mainly that (as shown before) the numerical fluctuations do not include the synchronous pressure fluctuations with the experimental dominant frequency of $f_2$ shown in Figure 7a. The comparison between the results of the different simulations also shows that the RMS values of $p_{Other}^{\prime }$ are very sensitive to the selected turbulence modeling technique. This type of pressure fluctuations includes the pressure fluctuations due to wakes of the guide vanes, runner blades, and runner crown which are captured to a larger extent in the WMLES simulation compared with the zonal WMLES and SAS simulations, for which the region upstream the draft tube is resolved using a RANS approach.

The comparison between the RMS values of the pressure fluctuations due to the RVR in Figure 9a showed that the predicted value from the WMLES simulation is closer to the experimental value than the values using the other turbulence modeling techniques. In order to investigate the reason for this and also study the effect of turbulence modeling on the pressure field in the draft tube, Figure 10 shows the phase-averaged pressure over Line 1 with the different turbulence modeling techniques. The core of the RVR, where the pressure is low, passes Line 1 at instances A and B. It can be seen that the RVR core pressure drops more in the WMLES results compared with the results of the other turbulence modeling techniques. A more quantitative comparison is shown in Figure 10d, at time instance B only. It can be seen that the pressure field far from the vortex core ($y>50$ mm) is almost the same in all simulations. However, the pressure in the near-field of the RVR ($y<0$ mm) is affected by the choice of turbulence modeling technique. In the WMLES results, the minimum pressure in the near-field is lower and the extent of the low pressure region is larger and slightly closer to the nearest wall compared with the other two results. This can explain the higher RMS level of the pressure fluctuations due to the RVR in the WMLES simulation as shown in Figure 9a. It should particularly be stressed that capturing the correct pressure drop in the RVR is very important for cavitating simulations as this pressure drop is the driving force of the cavitation formation. In WMLES simulations, where a larger pressure drop can be captured, a larger volume of cavitation should be expected in cavitating simulations.

To explain the reason for the larger pressure drop in the RVR of the WMLES simulation, Figure 11 shows the phase-averaged normal velocity over Line 1 at time instance B that was shown in Figure 10. The center of the RVR approximately corresponds to $U_n=0$ due to the rotating flow around the center of the RVR. The comparison between the numerical results shows that the gradient of the normal velocity $\frac{\partial U_n}{\partial y}$ is larger around the center of the RVR in the WMLES result than with the other turbulence model techniques. This indicates that the swirling motion around the RVR in the WMLES simulation is stronger, which leads to the larger pressure drop shown in Figure 10.

4.2. Effect of Cavitation

As mentioned in Section 3.1, cavitation simulations using WMLES are performed for three different cavitation numbers to study the effect of cavitation on the global quantities, the velocity profiles and the pressure fluctuations in the draft tube. It should be mentioned that since there are no experimental data for cavitating conditions in a Francis-99 turbine, no comparison is made between the simulation results and experimental ones. To show the extent of the cavitating region in these simulations, Figure 12 presents the cavitating part of the RVR using a blue iso-surface of ${\alpha }_v=0.99$. This figure also shows the variation of the total volume of vapor in the RVR as well as a spectral analysis of this variation. It can be seen that the cavitation inception at $\sigma =0.07$ happens at the root of the RVR near the runner crown. The total vapor content at this condition exhibits significant fluctuations, indicating that the cavitation is highly unstable. The spectral analysis of the vapor volume variation shows that this instability in the cavity volume does not have any dominant frequency. By decreasing the cavitation number to $\sigma =0.06$, the cavitation starts to incept in the small vortices further downstream the runner exit. Similar to the previous condition, the cavitation is unstable, although with a lower frequency. The spectral analysis indicates that although there is an increase in the PSD level of frequencies lower than $f=3f_{runner}$, this increase does not lead to a dominant peak in the PSD level. By further decreasing the cavitation number to $\sigma =0.05$, the cavitating region covers almost the entire root of the RVR near the runner exit. The variation of vapor volume indicates that there are fluctuations in the size of the cavitating region of the vortex. Unlike the other two cavitating conditions, the spectral analysis shows that these fluctuations have a dominant frequency at $f_3/f_{runner}=0.43$.

Table 3. Effect of cavitation number on the torque, $M_z$, head, H, and efficiency $\eta$.

$\mathit{\sigma }$	∞	$0.07$	$0.06$	$0.05$
$M_z$ (Nm)	439.68	439.02	439.00	438.99
H (m)	12.48	12.46	12.46	12.46
$\eta$	89.58	89.61	89.60	89.60

presents the global quantities in the form of torque, $M_z$, head, H, and efficiency $\eta$, for the different cavitation numbers and for the non-cavitating condition ($\sigma =\infty$). It can be seen that the cavitation number does not have any significant effect on these quantities, as the maximum variation in these quantities with respect to the cavitation number is less than 0.2 percent. This is expected, as the studied cavitation numbers are far from the cavitation breakdown for the studied turbine [48].

To study the effect of cavitation on the velocity field, Figure 13 shows the time-averaged velocity profiles on Lines 1–3 (shown in Figure 2) for different cavitation numbers. It can be seen that cavitation does not have any effect on these velocity profiles for $\sigma =0.06$ and $\sigma =0.07$. For these conditions, the size of the region with low values of absolute axial velocity is almost identical to the size of this region in the non-cavitating condition. In the fully cavitating RVR in the simulation with $\sigma =0.05$, however, the velocity profiles are slightly affected by the presence of cavitation. This effect is more dominant in the profiles for Line 1 as this line is closer to the cavitating part of the RVR.

In order to investigate the effects of cavitation on the structure of the RVR and the instantaneous velocity field, Figure 14 presents the iso-surface of the Q criterion (left plots) and the history of the axial velocity over Line 1 for one period of vortex rotation (right plots) for different cavitation numbers. The iso-surface of the Q-criterion shows that similar to the non-cavitating condition, shown in Figure 4c, the RVR of the cavitating conditions consists of small vortices and there is a reverse flow at the location of these small vortices. The plots of history of the axial velocity on Line 1 show that streaks are formed as these vorticies and their reverse flows pass Line 1. A comparison between the results for the different cavitation numbers shows that the reversed flow in these streaks is highly affected by the presence of cavitation. At $\sigma =0.07$, for which the amount of cavitation in the RVR is small, the reverse flow in the streaks is quite similar to the non-cavitating condition (see Figure 4c). As the amount of cavitation increases, for $\sigma =0.06$ and $\sigma =0.05$, the reverse flow in the streaks becomes weaker.

Figure 15a shows the pressure fluctuations at Probe 2 for different cavitation numbers, including their RMS values. It can be seen that the inception of cavitation (at $\sigma =0.07$) leads to spikes in the pressure fluctuations, which is due to the collapse of the cavitation region. This can be seen in Figure 15b, where the total vapor volume decreases to near-zero values at the time of the pressure spikes. Due to these spikes, the RMS of the pressure fluctuations increases by 54% compared to the non-cavitating case. It should be noted that the spikes are truncated in the plot in order to keep a scale that still shows the variations due to the RVR. By slightly increasing the amount of cavitation (for $\sigma =0.06$), the spikes are not as frequent as those at $\sigma =0.07$, which indicates that the cavitation region is less frequently entirely collapsing. This is confirmed in Figure 12b, where the total volume fraction goes to zero less frequently than in Figure 12a. There is however a further increase in the RMS value as the cavitation number is decreased from $\sigma =0.07$ to $\sigma =0.06$, indicating that the collapses of the larger cavitation regions give higher pressure pulses. Again, it should be noted that the spikes are truncated in the plot. In the case of the fully cavitating RVR, at $\sigma =0.05$, most of the spikes are gone. This indicates that the cavitation region never collapses entirely (confirmed in Figure 12c), and that the collapses of smaller cavitation regions in the freestream give much smaller pressure spikes. This leads to smaller RMS values. On the other hand, it can clearly be seen that the variations due to the RVR is much less periodic at $\sigma =0.05$ than for the other cavitation numbers, indicating that the general flow features are influenced to a larger extent. In accordance with the increase in the RMS value compared to the non-cavitating condition, the cavitation increases the amplitude of the variations due to the RVR as a major contributor to the RMS value.

Figure 16 shows spectral analyses of the pressure fluctuations at Probe 2 for the different cavitation numbers, as well as the phase difference between the pressure fluctuations at Probes 2 and 3. The spectral analyses show that the dominant frequency, $f_1/f_{runner}=0.3$, is not affected by the presence of cavitation. As mentioned earlier, this frequency is related to the frequency of the RVR rotation. A comparison between the non-cavitating and cavitating conditions indicates that cavitation mostly affects the PSD level of the higher frequencies rather than that of the relatively low RVR frequency. For $\sigma =0.07$, at cavitation inception, an increase can be seen in the PSD level of frequencies larger than $2.0f_{runner}$. By further increasing the amount of cavitation, at $\sigma =0.06$, the increase in the PSD level of the pressure fluctuations happens already at $f>f_{runner}$. At $\sigma =0.05$, the increase in the PSD level approaches the $f_1$ frequency, with a local peak at $f_3/f_{runner}=0.43$, and the PSD level of the higher frequencies again decreases. The frequency of the additional peak is the same as the dominant frequency of the vapor volume fluctuations for this cavitation number, as shown in Figure 12c. The phase difference between the pressure fluctuations at Probes 2 and 3, shown in Figure 16b, shows that for the frequencies where there is an increase in the PSD level due to cavitation, the phase difference is highly reduced (approaching zero). This means that the increased pressure fluctuations in these frequencies are synchronous.

Figure 17 shows the effects of cavitation on the different components of the pressure fluctuations and their RMS values, decomposed according to Equation (22). It can be seen that the cavitation has insignificant effects on the RMS values of the pressure fluctuations due to the RVR, as the maximum difference between the RMS values for the different cavitation numbers is around 7%. However, the RMS values of the other sources are highly affected by the presence of cavitation. Similar to the trends shown in Figure 15a, the RMS of the synchronous pressure fluctuations first increases as the amount of cavitation increases ($\sigma =0.07$ and $\sigma =0.06$), and then decreases when the RVR is fully cavitating ($\sigma =0.05$).

Figure 18 shows the effects of cavitation on the forces exerted on the runner and draft tube in the frequency domain. For the runner, only the z-component of the force (the blue curve) is affected by the cavitation, while for the draft tube, both the x- and z-components of the forces (red and blue curves, respectively) are affected by the cavitation. The trends of the changes due to cavitation, however, are the same for these affected force components, and they are quite similar to the trends for the pressure fluctuations shown in Figure 16a. At $\sigma =0.07$, where the amount of cavitation is small, there is an increase in the PSD level of the high-frequency fluctuations of the affected force components. The same increase can be seen for $\sigma =0.06$, although the increase in the PSD level starts to appear already at lower frequencies. At $\sigma =0.05$, this increase in PSD level leads to the dominant frequency $f_3/f_{runner}=0.43$, which is the same as the dominant frequency of the vapor volume fluctuations as shown in Figure 12c. It should be mentioned that the changes in the forces discussed here are caused by cavitation-induced pressure fluctuations, which are shown to be synchronous in Figure 16. Due to the synchronous nature of these pressure fluctuations, they affect only the forces in the directions where the geometry is asymmetrical. In the symmetrical directions, the changes in the forces due to these pressure fluctuations cancel each other out. For the runner, the geometry is almost symmetrical with respect to the x- and y-directions, and therefore, the cavitation-induced pressure fluctuations can affect only the forces in the z-direction. In the draft tube, however, the geometry is symmetric only with respect to the y-direction, and the effects of cavitation can be seen both in the x- and z-components of the forces.

5. Conclusions

In this work, we examine non-cavitating and cavitating RVR (Rotating Vortex Rope) in the Francis-99 turbine model using scale-resolving approaches. The non-cavitating simulations are performed using the SST-SAS, zonal wall-modeled LES and wall-modeled LES (WM) approaches, and the results are compared with the experimental data made available by the Francis-99 workshop. Furthermore, cavitating simulations are conducted for three cavitation numbers using the WMLES approach. The results from these simulations are used to study the effects of cavitation on the flow features, such as the velocity distribution and pressure fluctuations in the draft tube, and the forces acting on the draft tube and the runner.

The comparison between the results of the non-cavitating simulations and the experimental data reveals that the averaged velocity profiles and the pressure fluctuations predicted by the WMLES approach are in better agreement with the experimental data compared with the SST-SAS and zonal WMLES approaches. It is shown for the first time that the better velocity prediction is related to a correctly predicted rotating path of the RVR in the WMLES simulation. In the SST-SAS and zonal WMLES simulations, the RVR rotates on a path which is closer to the center of the draft tube as compared to the WMLES simulation and experiment. This yields an over-prediction of the influence of the RVR on the time-averaged velocity profile near the center of the draft tube. The comparison between the numerical results and the experimental data also shows that all methods under-predict the RMS values of the pressure fluctuations in the experiment, although the predicted RMS of the pressure fluctuations using the WMLES approach is closest to the experimental value. Using a detailed analysis of the pressure fluctuations in the simulations and experiment, the reason for the difference between the numerical and experimental values is shown to be related to a synchronous type of fluctuations only appearing in the experimental data, with a dominant frequency of around 2.8 times the frequency of the runner. An analysis of the pressure fluctuations in the experimental data shows that this synchronous pressure fluctuations can be seen also at the BEP and HL conditions, which suggests that these synchronous pressure fluctuations are created by a component in the system rather than the RVR. Such a finding has not been reported in the previous studies on the Francis-99 turbine. It is also shown that the WMLES simulation can capture a larger pressure drop in the center of the RVR compared to the other two approaches, which is shown to be due to the larger swirling velocity around the RVR in the WMLES simulation.

The results from the cavitating simulations reveal that cavitation slightly affects the average velocity profiles in the draft only if there is a large amount of cavitation in the RVR. They also show that the presence of cavitation damps out the instantaneous reverse flow in the small vortices in the RVR. The pressure fluctuations are also shown to be significantly affected by the presence of cavitation. Cavitation induces synchronous pressure fluctuations with a frequency larger than the frequency of the RVR. Capturing these synchronous pressure fluctuations in simulations has previously not been reported in the literature. When the amount of cavitation in the RVR is small, these fluctuations have a broadband high-frequency spectrum, while in the case of a fully cavitating RVR, they have a dominant frequency close to the dominant frequency of the RVR. The analysis of the forces reveals that the cavitation-induced pressure fluctuations have different effects on the forces exerted on the runner and draft tube. In the runner, the presence of cavitation induces significant force fluctuations only in the direction aligned with the rotational axis of the runner, while the force fluctuations in the draft tube are additionally affected in the direction of the bend. This finding can be used to design a method to detect cavitation in the turbine based on the direction of the vibrations in the draft tube and runner.