Results on flows in a draft tube of a constant-head, constant-specific speed, model Francis turbine are presented based on computational fluid dynamics (CFD) simulations and theoretical analysis. A three-dimensional, unsteady, Navier–Stokes solver with the detached-eddy simulation (DES) model and the realizable k–*ϵ* (RKE) model is used to analyze the vortex rope formed at different discharge coefficients. The dominant amplitude of the pressure fluctuations at a fixed point in the draft tube increases by 13 times, and the length of the rope increases by 3.4 times when the operating point of the turbine shifts from a discharge coefficient of 0.37 to 0.34. A perturbation analysis based on a steady, axisymmetric, inviscid, incompressible model for the mean flow is performed to obtain a Sturm–Liouville (SL) system, the solutions of which are oscillatory if the discharge coefficient is greater than 0.3635, and nonoscillatory otherwise.

## Introduction

Prior studies [1–3] on draft tube surge have shown that turbines operating away from the design point (or the best efficiency point (BEP)) are laden with undesirable effects in the form of violent pressure fluctuations caused by a helical vortex formed in the draft tube. This helical vortex, known as the vortex rope, is a complicated fluid structure composed of spinning fluid particles. Besides, the vortex rope itself revolves around the axis of the draft tube and is hence described as a “precessing” fluid structure with high vorticity. Typically originating at the nose cone, the vortex rope may appear as a conical outward spiral in the draft tube, where the flow is extremely turbulent and corresponds to high Reynolds numbers.

Hydroturbines operating at off-design conditions are characterized by high swirl levels at the inlet of the draft tube. Furthermore, the draft tubes also comprise a shear layer formed between a centrally stalled region (directly underneath the nose cone of the runner) typical of diminished axial velocities and fluid flow in the surrounding vicinity with substantially high axial velocities. When the swirling main flow combines with the shear layer, a vortex rope is formed in the draft tube. The core of the vortex rope is a region of low pressure, which sometimes may be inherently lower than the vapor pressure of water, leading to cavitation. In such cases, the vortex rope appears as a cavitated core. In cases that are not intrinsically cavitating, the vortex rope can be visualized in experiments by sufficiently decreasing the pressure at the outlet of the draft tube such that the local pressure in the vortex rope core drops below the vapor pressure of water, leading to cavitation (see Ref. [4] for the development of a vapor core rope in the draft tube cone of a hydroturbine during an experiment).

The vortex rope manifests itself by means of high-amplitude, low-frequency pressure fluctuations that are nearly periodic in nature due to its precession. These fluctuations are extremely intense and induce severe vibrations and noise in the draft tube, thereupon disrupting the functioning of the turbine and reducing its efficiency, which calls for increased maintenance. For example, see Refs. [5–7] for an estimate of these pressure fluctuations obtained under different operating conditions. Recently, Müller et al. [8] also performed experiments in a small-scale test rig (in the form of a conical draft tube of a microturbine) to measure the amplitudes of the draft tube wall-pressure fluctuations caused due to the precessing vortex rope. They also pointed out that these fluctuations caused by the oscillating vortex rope can play a crucial role on the stability of hydraulic machines. Therefore, it is worthwhile to study and analyze the vortex rope extensively to gain information that may be used in its mitigation, and to minimize its harmful effects on the turbine. Of late, efforts have been taken (see Refs. [9] and [10], for example) to mitigate the vortex rope by injecting water axially into the discharge cone.

### Previous Work.

The last few decades have witnessed a great deal of work on the vortex rope and its mitigation. Several numerical models have been used by researchers to simulate swirling flows and the vortex rope in a draft tube, both at full-load (BEP) and at partial-load conditions. Shyy and Braaten [11] used the steady Reynolds-averaged Navier–Stokes (RANS) equations closed with the *k–ϵ* model to analyze swirling flow in a draft tube, focusing on numerical aspects and on the draft tube pressure recovery. Ruprecht et al. [12] discuss the importance of turbulence models in capturing the effects of the vortex rope accurately and conclude that the standard *k–ϵ* model is incapable of predicting the vortex rope due to excessive damping. Therefore, they use the extended model of Chen and Kim [13], which was also used by Ogor et al. [14]. The poor performance of the two-equation *k–ϵ* model at partial-load conditions was also pointed out by Vu et al. [15]. The two-equation shear stress transport (SST) model [16] and the one-equation eddy-viscosity model [17] were used to simulate swirling flows by Yaras and Grosvenor [18], who concluded that these two models are inappropriate for use in simulating strongly swirling confined flows.

For complex three-dimensional unsteady flows, the use of direct numerical simulation (DNS) is preferred since it offers the highest accuracy. However, it is not practical for industrial applications because of its demand for extensive computational resources. The next preference is large-eddy simulation (LES) or very large-eddy simulation (VLES), also known as detached-eddy simulation (DES). The DES is a promising model for use in unsteady turbulent flows and is reasonable in terms of computational time and cost. The DES has been preferred over RANS-based models by several researchers including Ruprecht et al. [12], Ogor et al. [14], Paik et al. [19], Foroutan and Yavuzkurt [20], and Foroutan and Yavuzkurt [21]. Due to the inherent limitation of RANS simulations in capturing unsteady features of swirling flows, Guo et al. [22] used LES to analyze the cavitating and noncavitating flows in the draft tube of a Francis turbine operating at partial-load conditions. To analyze the vortex rope formed in the draft tube of a model Francis turbine at two specific flow rates, which correspond to 70% and 91% of the flow rate at the BEP, Foroutan and Yavuzkurt [21] used four different models, namely, the standard *k–ϵ* model [23], the realizable *k–ϵ* (RKE) model [24], the SST *k–ω* model [25], and the DES model based on the SST *k–ω* model [26]. They concluded that the DES model captures the structure of the vortex rope much better than the other models considered. Foroutan and Yavuzkurt [27] also tested the performances of a partially averaged Navier–Stokes (PANS) model and the traditional RANS model in simulating the flow in the draft tube of a Francis turbine for four different operating conditions. Based on their results, they concluded that though both the RANS and PANS models performed well in predicting the flow when the operating point of the turbine was close to the BEP, the results of the RANS model did not agree well with the experimental measurements, for conditions in which the operating point of the turbine was not close to the BEP. Recently, Javadi et al. [28] studied the precessing vortex rope in a conical diffuser, with rotor–stator interactions, experimentally and numerically, for different runner rotational speeds in the range 400–920 rpm, and a constant flow rate of 30 l/s. For their numerical simulations, they used the delayed DES method [29], conjugated with the Spalart–Allmaras turbulence model [17], which they demonstrate to be effective in predicting the flow features.

Although considerable results on the vortex rope have been obtained using several computational models, the analysis of swirling flow in the draft tube was simplified by Susan-Resiga et al. [30], who used an axisymmetric model to compute circumferentially averaged pressure and velocity fields. The axisymmetric approximation amounts to treating the flow field as two-dimensional, and restricting the analysis to the symmetry plane, which reduces the computational effort. However, due to its fundamental simplicity, this model cannot capture the unsteady features of the flow. The axisymmetric model renders a steady, averaged flow picture of an extremely unsteady phenomenon and is therefore incapable of predicting the vortex rope location and its structure. Nonetheless, in spite of being incompetent in capturing the vortex rope, the axisymmetric model is capable of recovering the circumferentially averaged flow field associated with the full three-dimensional unsteady swirling flow. Therefore, the axisymmetric model could also be used in parametric studies to provide a background for theoretical analysis. Susan-Resiga et al. [31] used the axisymmetric model to perform an eigenvalue analysis of the linearized governing equations and showed that the swirl reaches a critical state where a sudden variation in draft tube pressure recovery is observed. To analyze swirling flows in the draft tube and to obtain experimental data, a popular project dealing with the flow investigation in draft tubes (FLINDT) was developed [32]. Susan-Resiga et al. [33] modified the draft tube used in the FLINDT project to obtain a geometry in the form of an axisymmetrical straight diffuser.

*n*= 0.56. They conducted a parametric study to analyze the effect of the runner angular speed,

_{s}*ω*, on the velocity field at the runner outlet, which is also the inlet of the draft tube. They carried out this analysis by changing the runner angular speed,

*ω*, which in turn changes the characteristic Reynolds number, Re, defined as

*ν*is the kinematic viscosity of water and

*R*

_{ref}is a reference radius, which is typically the outlet radius of the runner. Consequently, a runner of outlet diameter 0.4 m, rotating at an angular speed of 1000 rpm, would correspond to the flow of water ($\nu =10\u22126m2/s$) at a Reynolds number of about 8.4 × 10

^{6}. To determine the parameters which influence the draft tube inlet velocity profiles, Susan-Resiga et al. [31] also varied the turbine operating point, which is specified by the discharge coefficient, $\varphi $, and the specific energy coefficient,

*γ*, defined as

*H*is the turbine head,

*ω*is the runner angular speed,

*g*is the acceleration due to gravity, and

*γ*< 1.3. Hence, they concluded that within this range, the discharge coefficient, $\varphi $, is the only relevant parameter. Based on their experimental measurements, they also modeled the inlet velocity field as the superposition of two Batchelor vortices on a base flow in the form of a solid body rotation. They further identified that one of the vortices is counter-rotating (and coflowing) with respect to the base flow while the other is corotating (and counterflowing) relative to the solid body rotation. Thus, the velocity field at the inlet of the draft tube may be described as a sum of three components so that the circumferential component of velocity,

*U*, and the axial component of velocity,

_{θ}*U*, are

_{Z}*R*∈ [0,

*R*) is the radial coordinate,

_{i}*R*=

_{i}*R*

_{ref}= 0.2 m is the radius of the draft tube inlet, and

*ω*,

_{i}*U*, and

_{i}*R*are coefficients that are functions of $\varphi $ (see Refs. [31] and [33]). Further, Susan-Resiga et al. [33] also noted that the discharge coefficient corresponding to the BEP of the Francis turbine is $\varphi BEP\u22480.37$. A dimensionless parameter

_{i}*ξ*may be defined as the ratio of the flow rates at the operating point of the turbine and the BEP. Mathematically

where the flow rate, $Q(\varphi )$, is defined in Eq. (3).

### Scope and Outline of the Present Work.

Though several studies have claimed that stronger vortex ropes are obtained when turbines operate at off-design conditions, it is not known how the structure of the vortex rope changes with varying operating conditions. Additionally, it is also of interest to know the conditions under which a vortex rope may be classified as weak and strong, and to estimate the critical condition at which a transition may be expected. Such a condition may be obtained from experiments, CFD simulations, and/or mathematical analyses. An estimate of such a critical condition helps in deciding the range within which the turbine may be operated, so that the undesirable effects of the vortex rope are minimal. We note that Susan-Resiga et al. [31] performed an eigenvalue analysis using the axisymmetric approximation, to estimate such a critical condition. One goal of this work is to show that a reduced-order mathematical analysis (here, we use a SL formulation) is capable of predicting the critical condition, whose validity may be verified by results obtained from the CFD simulations, which reveal the structure of the vortex rope and the range of pressure fluctuations corresponding to different operating conditions.

*γ*(for $1.0\u2264\gamma \u22641.3$), we present results only for varying discharge coefficients, $\varphi $. Therefore, here we consider swirling flows in the discharge cone of a Francis turbine whose runner outlet has a radius of

*R*

_{ref}= 0.2 m, operating at a constant speed of

*ω*= 1000 rpm, and a constant head,

*H*, such that

*γ*= 1.18. Note that fixing these two parameters fixes the specific speed of the turbine,

*n*, defined as

_{s}*b*), we must specify the functions $Ui(\varphi )$ and $Ri(\varphi )$. Using regression analysis on the data collected by digitizing Fig. 2 of Ref. [33], we obtain

*U*'s are in units of (m/s) and the

_{i}*R*'s are in units of (m). Similarly, we also obtain $\omega i(\varphi )$ using regression analysis on the data collected by digitizing Fig. 2 of Ref. [33]

_{i}where the *ω _{i}*'s are in units of (rad/s). For a runner speed of

*ω*= 1000 rpm = 104.72 rad/s, and a reference radius of

*R*

_{ref}= 0.2 m, a turbine operating at a specific energy coefficient of

*γ*= 1.18 has a constant head of

*H*≈ 26.4 m (according to Eq. (2

*a*)). From Eqs. (3), (7), and (8), we obtain the flow rate at the BEP, $QBEP=Q(\varphi \u22480.37)\u22480.143m3/s$. Using these values in Eq. (6), we obtain

*n*≈ 0.6.

_{s}Here, we analyze swirling flows operating at partial load, the flow rates for which correspond to $0.919\u2264\xi \u22641.0(0.34\u2264\varphi \u22640.37)$. We are interested in investigating the effects of the vortex rope on the draft tube, using the CFD simulations and a simple theoretical analysis motivated by Susan-Resiga et al. [31]. To this end, (i) in Sec. 2, we present the results from the CFD simulations, which include the pictorial depictions of the vortex rope, estimates of the pressure fluctuations, and the strength of the rope determined by its geometrical parameters, and (ii) in Sec. 3, we follow Susan-Resiga et al. [31] and present a theoretical analysis which uses standard perturbation theory to formulate an SL system, based on which the flow may be classified into regions of oscillatory (and bounded) solutions and regions of nonoscillatory (and unbounded) solutions.

## Numerical Analysis

Here, we present the results from the CFD simulations. We are interested in the effect of the discharge coefficient on the structure of the vortex rope and on the pressure fluctuations caused by the vortex rope. The geometry used, the details of the mesh, the boundary conditions, and other computational details are presented in Sec. 2.1. The validity of the computational scheme used is established in Sec. 2.2, in which the numerical results are compared with the experimental data available from Ref. [33]. The analysis of the pressure fluctuations is presented in Sec. 2.3, and the structures of the vortex ropes obtained for different discharge coefficients are presented in Sec. 2.4.

### Computational Details.

For the computational analysis, we consider the simplified draft tube used by Foroutan and Yavuzkurt [21]. The simplified draft tube, as shown in Fig. 1, is typically a conical diffuser of angle 17 deg and a length-to-inlet diameter ratio of 1.685, followed by a discharge section of constant diameter. Foroutan and Yavuzkurt [21] also performed a detailed grid sensitivity analysis for this problem of analyzing swirling flow in the simplified draft tube and observed that grids of sizes smaller than that used for their analysis led to errors as high as 10%. Therefore, for the numerical simulations, we use the mesh of Foroutan and Yavuzkurt [21], which is shown in Fig. 1. The computational mesh of Foroutan and Yavuzkurt [21], used here, has about 2 × 10^{6} cells and is finer in the central core region and near the walls to capture high gradients due to the vortex rope and the walls, respectively.

^{®}[34]. We used a three-dimensional, transient, pressure-based Navier–Stokes solver, with the DES model, along with the RKE model. Following Susan-Resiga et al. [33], we specified inlet profiles for the three velocity components, the turbulent kinetic energy,

*k*, and the turbulent dissipation rate,

_{t}*ϵ*. For the circumferential and axial velocity profiles at the inlet, we used Eq. (4), along with Eqs. (7)–(9). Susan-Resiga et al. [33] pointed out that their measurements of the radial velocity profiles are inaccurate due to the low magnitudes; therefore, they recommend a simple linear fit (with a value of zero at the geometric center of the inlet and a value of 0.15 ×

_{t}*U*

_{Z}_{,wall}close to the wall) for the radial component of velocity. So, we specified an inlet radial velocity that varies according to

*R*∈ [0,

*R*) and

_{i}*R*=

_{i}*R*

_{ref}= 0.2 m is the radius of the draft tube inlet. The number 0.15 in Eq. (10) is due to the geometry of the draft tube: Since the half-angle of the draft tube cone is 8.5 deg, the ratio

*U*/

_{R}*U*< tan(8.5 deg) ≈ 0.15. Susan-Resiga et al. [33] considered a four-parameter expression for the turbulent kinetic energy,

_{Z}*k*, so that

_{t}*R*,

_{B}*R*,

_{C}*k*, and

_{B}*k*are functions of the discharge coefficient, $\varphi $. Susan-Resiga et al. [33] also provide experimental data based on which one can obtain these functions. We determined these functions by digitizing Figs. 2 and 3 of Ref. [33], and using regression analysis, so that

_{C}*R*and

_{C}*R*are in units of (m), and

_{B}*k*and

_{C}*k*are in units of (m/s)

_{B}^{2}. Finally, for the turbulent dissipation rate,

*ϵ*, at the inlet, we used the profile recommended by Susan-Resiga et al. [33]

_{t}with *C _{μ}* = 0.09 and

*l*= 0.02

_{r}*R*.

_{i}At the outlet of the draft tube, we specified a radial pressure condition, $\u2202P/\u2202R=\rho (U\theta )2/R$. At the walls, we used a no-slip condition. Using the boundary conditions described above, we simulated the unsteady swirling flow (using the DES model available in ansys-fluent^{®} [34], along with the RKE scheme for the RANS model) in the simplified draft tube for four different discharge coefficients in the range $0.34\u2264\varphi \u22640.37$. The discretization scheme is second-order accurate in both time (implicit) and space. For numerical simulation of flow in hydraulic turbines, Ciocan et al. [35] recommend using a time step equivalent to a single degree of rotation of the runner. In our case, since the speed of the runner is *ω* = 1000 rpm, it takes about 1.667 × 10^{−4} s to rotate through 1 deg. Therefore, in the transient numerical simulations, we used a time step of Δ*t* = 10^{−4} s, which is less that the time taken by the runner to rotate through 1 deg. A similar value for the time step was used by Foroutan and Yavuzkurt [21] in their numerical simulations with the DES model. Here, we are interested in analyzing the swirling flow, which is periodic due to the precession of the vortex rope. To analyze the flow in the periodic, quasi-steady state, we first obtained a steady-state solution and used this solution as the initial condition for the transient (periodic, quasi-steady) simulations, as recommended by Ciocan et al. [35]. When using such a steady-state solution as the initial condition for a transient simulation, Ciocan et al. [35] also pointed out that “the transient solution should converge to the desired periodic behavior after several runner rotations.” We observed from the numerical simulations that the periodic, quasi-steady state was achieved in less than 12,000 time steps, which corresponds to a physical time of 1.2 s or 20 runner rotations. Therefore, we present the results obtained during the time period 2–10 s. Within each time step, we allowed for 20 iterations to obtain a converged solution for that particular time step. For every time step, we observed that the residuals dropped to values less than 10^{−8} at the end of 20 iterations. We used ansys-cfx^{®} [36] to postprocess the results obtained.

### Validation of the CFD Results.

To validate the results obtained using the CFD simulations, we used the experimental data of Susan-Resiga et al. [33], who provided the steady (in the mean), axisymmetric (circumferentially averaged) values of the axial and circumferential velocities at a horizontal section 158 mm (see Fig. 1) downstream of the draft tube inlet. We consider two lines, *L*_{1} and *L*_{2} (see Fig. 1), along which we obtain the axial and circumferential velocities, at two different instants of time, *t* = 6 s and *t* = 10 s. In Fig. 2, we present the comparisons of the results of the numerical simulations with the experimental data of Susan-Resiga et al. [33]. We observe that the numerical results obtained for the four different cases (two different time instants and two different directions along the *X*- and *Y*-axes) considered are similar to one another. We also note that the agreement between the experimental and numerical results is reasonably good.

### Pressure Fluctuations.

*P*is an appropriate measure of the range of pressure fluctuations (and

*a*is the dominant amplitude of the pressure fluctuations),

_{P}*P*

_{ref}is a suitable reference pressure, which is constant and could be assumed to be unity without loss of generality,

*L*is the length of the draft tube,

*D*is the diameter of the draft tube inlet,

*n*is the specific speed of the (scaled model) of the turbine, Fr is the Froude number (ratio of the inertial force on the fluid to the gravitational force on it), and We is the Weber number (ratio of the inertial force on the fluid to the force on it due to surface tension). For the range of discharge coefficients ($0.34\u2264\varphi \u22640.37$) considered here, there is no cavitation, and hence, the effect of

_{s}*σ*can be neglected. The effect of Fr and We is assumed to be minor (see Ref. [1]). Similarly, the effect of Re is disregarded, since the flow in a conventional turbine is highly turbulent. In the present analyses, we consider fluid flow in a fixed draft tube of a scaled model of a Francis turbine operating at a constant specific speed of 0.6, and the runner of which rotates at a fixed rate of 1000 rpm. This further eliminates the dependence of Δ

*P*on

*L*/

*D*and

*n*. The functional form, under these conditions, then reduces to

_{s}Therefore, the range of pressure fluctuations varies only with the discharge coefficient. For the purpose of analyzing the pressure fluctuations, we collect data at a specific point, *P* (see Fig. 1), where the effects of the vortex rope are expected to be maximum. After eliminating the initial unsteady variations occurring within the first 2 s, we show the results after the flow becomes stationary periodic.

Figure 3 shows the time series of pressure data (at point *P*) and the corresponding frequency spectra for four different discharge coefficients. The dominant amplitude, *a _{P}*, and frequency,

*f*, obtained from the frequency spectra are listed in Table 1, and the range, Δ

_{P}*P*, of pressure fluctuations and the average pressure,

*P*

_{avg}, are shown as functions of the discharge coefficient in Fig. 4. From Table 1 and Fig. 4, we see that lower discharge coefficients indicate a higher mean pressure and range of pressure fluctuations, suggesting that the vortex rope is much stronger at low discharge coefficients, and that it decreases in strength as the operating point of the turbine moves closer to the BEP. From Table 1, we also see that the amplitude of pressure fluctuations increases by an order of magnitude (approximately 13 times) as the operating point of the turbine shifts from the BEP ($\varphi =0.37$) to $\varphi =0.34$, and likewise, the mean pressure increases roughly by a factor of 2. These numbers indicate the strength of the vortex rope at $\varphi =0.34$ relative to that at the BEP. The dominant frequencies of the pressure fluctuations obtained for all four values of the discharge coefficients are small (less than 20%) when compared to the frequency of rotation of the runner, which is 16.67 Hz. These values are consistent with the condition that the vortex rope is a high-amplitude, low-frequency phenomenon. An actual estimate of the strength of the vortex rope may be obtained by visualizing the vortex rope via isocontours of a few derived variables, and by measuring the size of the rope, as discussed in Sec. 2.4.

$\varphi $ | P_{avg} | a_{P} | f_{P} | ΔP = 2a_{P} | $\Delta P(\Delta P)BEP$ |
---|---|---|---|---|---|

(Pa) | (Pa) | (Hz) | (Pa) | ||

0.34 | −6596.5 | 126.10 | 2.50 | 252.20 | 12.9 |

0.35 | −8721.3 | 77.78 | 1.25 | 155.56 | 7.9 |

0.36 | −10,834.2 | 33.50 | 3.25 | 67.00 | 3.4 |

0.37 | −11,571.2 | 9.76 | 0.25 | 19.51 | 1.0 |

$\varphi $ | P_{avg} | a_{P} | f_{P} | ΔP = 2a_{P} | $\Delta P(\Delta P)BEP$ |
---|---|---|---|---|---|

(Pa) | (Pa) | (Hz) | (Pa) | ||

0.34 | −6596.5 | 126.10 | 2.50 | 252.20 | 12.9 |

0.35 | −8721.3 | 77.78 | 1.25 | 155.56 | 7.9 |

0.36 | −10,834.2 | 33.50 | 3.25 | 67.00 | 3.4 |

0.37 | −11,571.2 | 9.76 | 0.25 | 19.51 | 1.0 |

### Structure of the Vortex Rope.

As discussed in Sec. 1, the vortex rope is formed in swirling flows due to high shear and appears in discharge cones as a conical helix. In principle, the vortex rope is a region of low pressure, high vorticity, and high swirl strength. Since it is a zone of low pressure, it is common to visualize the vortex rope via an isocontour of pressure, at an appropriate contour level. Additionally, there are also other ways to visualize the vortex rope. These are essentially by virtue of derived variables (from the velocity and vorticity field), which are characteristic of the vortex core region and effective in illustrating the rope. We consider five such variables here (also see Ref. [36]). They are

- (i)
parameter

*q*: the second invariant of the velocity gradient tensor, defined in Eq. (21) - (ii)
swirling discriminant, Δ: the discriminant of velocity gradient tensor for complex eigenvalues, defined in Eq. (23)

- (iii)
swirling strength,

*s*: the imaginary part of complex eigenvalues of velocity gradient tensor, defined in Eq. (24) - (iv)
parameter

*λ*_{2}: the negative value of the second eigenvalue of the symmetry square of velocity gradient tensor, defined in Eq. (26) - (v)
real eigen helicity,

*h*: the dot product of vorticity and the real eigenvector of velocity gradient tensor, defined in Eq. (27)

**R**= [

*r*] are the velocity gradient tensor, the strain-rate, and rotation-rate tensors, respectively, and if

_{ij}*s*may be defined as

*λ*

_{2}is defined as the second eigenvalue of the sum of squares of the strain-rate and rotation-rate tensors. Mathematically

**u**and

**e**are the velocity vector and the real eigenvector of

**D**, then the real eigen helicity is defined as

Figure 5 shows the vortex rope formed in the draft tube at *t* = 10 s, visualized by all the methods described previously, for a discharge coefficient of 0.34. The contour levels of all the variables are chosen carefully after several trials, so that they all correspond to approximately the same rope length and angle (as defined in Fig. 7). It is observed that the isopressure contour does not provide an adequate estimate of the vortex rope structure and results in a premature cutoff of the rope tail. In other words, even though the rope is defined as a region of low pressure, the structure of the rope cannot be fully pictured by means of an isopressure surface, as shown in Fig. 5. The derived variables, on the other hand, are more effective in depicting the vortex rope. It is also worth mentioning that both large- and small-scale features of the vortex rope are ideally illustrated when it is visualized via these variables, and the quantitative structure of the rope is nearly the same in all cases. Therefore, presenting results for $0.34\u2264\varphi \u22640.37$ using all of these methods would be repetitive. So, in the parametric study, we present visualizations of the vortex rope obtained using the *q*-criterion, because it is mathematically easier to calculate and is one of the commonly used methods for the analysis and visualization of vorticity-dominated flows.

The vortex rope visualized via isocontours of “*q*” is shown in Fig. 6, for four different discharge coefficients. All the results are shown for the same contour level (*q* = 1345 1/s^{2}). A qualitative observation reveals that the size of the vortex rope increases as the discharge coefficient decreases. While the rope at $\varphi =0.37$ (BEP) is very thin and restricted to a relatively small region around the axis of the draft tube, the vortex rope at $\varphi =0.34$ spreads out in the radial and axial directions and occupies a larger portion of the draft tube. These results are consistent with the variation in the mean pressure and the range of pressure fluctuations shown in Fig. 4. To obtain a quantitative estimate of the rope strength, we measure the rope length and the angle as described in Fig. 7. Since the geometry of the rope enables it to be inscribed in a cone of a specified length and angle, the strength of the rope can be defined in terms of the dimensions of the cone. From Fig. 7, the half-angle of the rope, denoted by *α*, is $\u2220BAC$. While the angle can be clearly defined for the rope, there are different ways to define its length. An appropriate definition of the length of the rope would be the distance measured axially from the inlet to the midpoint of the last visible coil of the rope. Typically, the length is measured such that portions of the coils (if any) occurring beyond the measured length are not taken into account for calculating the angle of the rope. From Fig. 7, the length of the rope, denoted by *l*, is the distance *OC*. To obtain accurate estimates of the angle, *α*, and the length, *l*, we use the results from two of the variables, the swirling strength, *s*, and *q*, denoted by subscripts *s* and *q*, respectively. The rope angle and length are then obtained as the mean values of these two results. The values obtained for the rope angle and length are shown in Table 2 for the four discharge coefficients considered. The variation in the rope length is similar to that in the average pressure and the range of pressure fluctuations, as observed in Fig. 8.

$\varphi $ | l (m)_{s} | l (m)_{q} | l (m) | α (deg)_{s} | α (deg)_{q} | α (deg) |
---|---|---|---|---|---|---|

0.34 | 0.612 | 0.609 | 0.611 | 17.76 | 17.99 | 17.88 |

0.35 | 0.432 | 0.432 | 0.432 | 15.71 | 15.47 | 15.59 |

0.36 | 0.268 | 0.265 | 0.267 | 12.16 | 12.79 | 12.48 |

0.37 | 0.183 | 0.179 | 0.181 | 09.97 | 10.45 | 10.21 |

$\varphi $ | l (m)_{s} | l (m)_{q} | l (m) | α (deg)_{s} | α (deg)_{q} | α (deg) |
---|---|---|---|---|---|---|

0.34 | 0.612 | 0.609 | 0.611 | 17.76 | 17.99 | 17.88 |

0.35 | 0.432 | 0.432 | 0.432 | 15.71 | 15.47 | 15.59 |

0.36 | 0.268 | 0.265 | 0.267 | 12.16 | 12.79 | 12.48 |

0.37 | 0.183 | 0.179 | 0.181 | 09.97 | 10.45 | 10.21 |

The average values (bold) are graphed in Fig. 8.

It is observed from Table 2 that the length of the rope increases by a factor of 3.35, as the discharge coefficient decreases from 0.37 to 0.34. Likewise, the angle of the rope increases by a factor of 1.75. These numbers confirm the fact that the rope gets much stronger and leads to severe pressure fluctuations as the operating point of the turbine shifts farther from the BEP. The farther the point of operation, the stronger the rope and the greater its effects.

## Perturbation Analysis

*R*,

*θ*,

*Z*), neglecting the body forces, are

*U*,

_{R}*U*,

_{θ}*U*) are the dimensional radial, tangential, and axial velocities, respectively,

_{Z}*P*is the dimensional pressure, and

*ρ*is the density of the fluid. Since axisymmetric flows enable analysis in a two-dimensional plane, a stream function,

*ψ*, may be defined such that

are strictly functions of *ψ*, by Bernoulli's theorem and Kelvin's theorem, respectively.

The SLE has been used extensively in the analysis of swirling flows and flows involving vortex tubes. Leibovich [45] dealt with the stability and breakdown of incompressible vortices at high Reynolds numbers using the SLE and estimated the point of breakdown for a given set of conditions. Lopez [46] used the SLE to obtain bifurcation diagrams for axisymmetric flow in a constricted pipe in a two-dimensional parameter space comprising the Reynolds number and the relative swirl of the incoming flow. Wang and Rusak [47] studied the linear stability of solid body rotation and general swirling flows in a constant cross section, finite-length pipe, based on an eigenvalue analysis, and used the SLE to show that stability characteristics of the flow change at the critical swirl number. Wang and Rusak [48] also studied the stability characteristics and the time-asymptotic behavior of such flows using the unsteady SLE. Gallaire and Chomaz [49] used the unsteady SLE to study the stability properties of swirling flows in a pipe of finite length and found that the boundary conditions drive the instability of the flow in the limit of long, but finite pipes. Susan-Resiga et al. [31] used the SLE to perform an eigenvalue analysis of the swirling flow in draft tubes and used the concept of critical swirl to explain the drop in the pressure recovery coefficient observed close to the BEP in the scaled model of the Francis hydroturbine.

*ψ*

_{0}= 0 on the axis of the draft tube. Equation (36) is the function corresponding to the base flow. For their analysis, Susan-Resiga et al. [31] added a perturbation to the base flow

*K*is the axial wavenumber of this perturbation,

*ϵ*≪ 1 is a small parameter. They substituted Eq. (37) in Eq. (35) and linearized the system to obtain a second-order ordinary differential equation (ODE) for $\psi \u0303(Y\u0302)$

*E*and

*I*and used Euler's equation, $dP/dR=\rho U\theta 2/R$, eventually transforming Eq. (39) to

where *U _{θ}* and

*U*are defined in Eqs. (4

_{Z}*a*) and (4

*b*). Then, they performed an eigenvalue analysis to obtain the discharge coefficient at which the swirl reaches a critical stage. Here, we use Sturmian theory to obtain similar results.

Equation (38) represents an SL system for $\psi \u0303$ in the domain $0\u2264Y\u0302\u2264Y\u0302max$ ($Y\u0302max$ is chosen here to correspond to the edge of the wall boundary layer, in order to avoid viscous effects). The solutions of the SL system depend on the sign of the term in brackets, namely, $C(Y\u0302)+K2/2Y\u0302$. Using Sturmian theory [50], we obtain two cases. The first case is $C(Y\u0302)+K2/2Y\u0302>0$. In this case, the solutions of Eq. (38) are exponential, and $Y\u0302\u22600$ at any point in the domain. The amplitude of the perturbation grows; therefore, these are unstable solutions of the SL system. The second case is $C(Y\u0302)+K2/2Y\u0302<0$. In this case, there are one or more zeroes for $Y\u0302\u2208[0,Y\u0302max]$, and the solutions of Eq. (38) are oscillatory. Thus, these are stable solutions of the SL system. The condition $C(Y\u0302)+K2/2Y\u0302=0$ represents an inflection point. However, it is very rare that this condition holds for the entire domain; it may occur locally at a finite number of points. The coefficients $(Ui,Ri,\omega i)$ are functions of $\varphi $ (see Eqs. (7)–(9)), which implies that the function $C(Y\u0302)+K2/2Y\u0302$ depends on $\varphi $. To obtain the regions of stability and instability, we plot the function $C(Y\u0302)+K2/2Y\u0302$ for various discharge coefficients, $\varphi $.

The parameter *K* is the eigenvalue of the system and may be obtained numerically, after discretizing Eq. (38) through a suitable scheme, such as finite differences. Thus, we obtain as many eigenvalues as the size of the system. The largest of these eigenvalues is the most influential and is hence used in the computations. Figure 9 shows the dimensionless quantity $c(y\u0302)+k2/2y\u0302$ as a function of $y\u0302$ over a range of discharge coefficients, where we have used the length scale *R*_{ref} and the time scale $\omega \u22121$ to nondimensionalize *C*, *K*, *R*, and $Y\u0302$. (Lowercase letters *c*, *k*, *r*, and $y\u0302$ denote dimensionless parameters corresponding to the respective uppercase letters.)

Figure 9 indicates that $c(y\u0302)+k2/2y\u0302<0\u2200y\u0302\u2208[0,y\u0302max]$ if $\varphi \u22730.3635$. Thus, for $\varphi \u22730.3635$, one may expect unconditional stability. On the other hand, if $\varphi \u22720.3635$, one may expect instability (growing perturbations). At $\varphi =0.37$, which may be classified under the stable regime, the vortex rope is thin (see Fig. 6), whereas for $\varphi \u22640.36$, which may be classified under the unstable regime, the vortex rope is relatively larger in size and much stronger.

It is also observed from Fig. 9 that there are no inflection points for $\varphi \u22730.3635$. On the contrary, there is at least one inflection point for $\varphi \u22720.3629$, and there are at least two for $0.3629\u2272\varphi \u22720.3635$. Thus, the function $c(y\u0302)+k2/2y\u0302$ changes sign in $y\u0302\u2208[0,y\u0302max]$ once or twice accordingly, and there are subregions where it is positive or negative. As a result, for $\varphi \u22640.3635$, we obtain conditional instability. For $\varphi \u22640.3629$, most of the domain falls in the unstable zone, except for a small finite region around the axis. The magnitude, $|c(y\u0302)+k2/2y\u0302|$, gives an estimate of how large $\psi \u0303$ is, since $\psi \u0303\u223c\u2009exp(mY\u0302)$, where $m\u221d|c(y\u0302)+k2/2y\u0302|$. The inflection points are locally neutral points, and slight deviations at these points may result in the system being pushed into a stable or unstable zone. Based on their eigenvalue analysis, Susan-Resiga et al. [31] characterize the flow as subcritical and supercritical for discharge coefficients less than and greater than 0.365, respectively. The present analysis using Sturmian theory yields a critical discharge coefficient that is within 0.4% of the value obtained by Susan-Resiga et al. [31], thus reaffirming their result.

The results of the theoretical analysis are consistent with those of the CFD simulations: The CFD simulations show that both the size of the vortex rope and the dominant amplitude of the pressure fluctuations increase as the operating point of the turbine moves away from the BEP; the theoretical analysis indicates that there is a transition from stability to instability at the critical discharge coefficient, where the stability occurs close to the BEP and the instability occurs at discharge coefficients away from the BEP. Both these trends show that the undesirable effects of the vortex rope on the draft tube increase as the operating point of the turbine shifts away from the BEP.

## Closure

In this work, we have investigated swirling flows in the simplified draft tube of a scaled model of a Francis turbine operating at a constant head of 26.4 m, whose runner rotates at a constant speed of 1000 rpm. Four discharge coefficients are considered for the numerical analysis. The dominant amplitude of the pressure fluctuations at a fixed point, *P*, on the draft tube wall, and the size (length and half-angle) of the vortex rope are obtained for each of these cases. As the operating point of the turbine shifts from $\varphi =0.37$ (BEP) to $\varphi =0.34$: (i) the dominant amplitude of the pressure fluctuations increases by approximately 13 times, (ii) the length of the rope increases by about 3.4 times, and (iii) the half-angle of the rope increases from a value of about 10.2 deg to about 17.9 deg. From these numbers, one may understand how strong a vortex rope gets as the operating point of the turbine shifts from the BEP. It is also found that pressure is not a suitable parameter for visualizing the rope. Nevertheless, the other chosen parameters depict the rope completely, and the structures of the vortex rope obtained using each of these variables are nearly identical to one another when comparing both large-scale and small-scale features.

The theoretical analysis indicates that the structure of the solutions to the SL system changes at $\varphi \u22480.3635$. This discharge coefficient, $\varphi \u22480.3635$, may be treated as a critical condition above which the solutions are oscillatory. Correspondingly, the swirling flow in the draft tube may be expected to yield a weaker rope. For discharge coefficients less than 0.3635, the solutions of the SL system grow exponentially. Correspondingly, for smaller discharge coefficients (away from the BEP), one may expect that the vortex rope gets strong enough, leading to large pressure fluctuations. It is therefore recommended that this turbine be operated at a discharge coefficient closer to the BEP, in order to avoid the undesirable effects due to a stronger vortex rope.

## Acknowledgment

The authors acknowledge the Research Computing and Cyberinfrastructure Unit of Information Technology Services at The Pennsylvania State University for providing advanced computing resources and services that have contributed to the research results reported in this paper. This work is affiliated with the hydropower research program at Penn State and is funded by a grant from DOE (Award No. DE-EE0002667).

## Nomenclature

*a*=_{P}dominant amplitude of the pressure fluctuations (Pa)

*D*=diameter of the draft tube inlet (m)

**D**=velocity gradient tensor

**e**=real eigenvector of

**D***E*=total specific energy (m

^{2}/s^{2})*f*=_{P}dominant frequency of the pressure fluctuations (Hz)

- Fr =
Froude number (dimensionless)

*g*=acceleration due to gravity (m/s

^{2})*h*=real eigen helicity (1/s)

*H*=turbine head (m)

*I*=circulation function (m

^{2}/s)*k*=axial wavenumber,

*KR*_{ref}(dimensionless)*K*=axial wavenumber (1/m)

*k*=_{t}turbulent kinetic energy (m

^{2}/s^{2})*l*=length of the vortex rope (m)

*L*=length of the draft tube (m)

*n*=_{s}specific speed of the turbine (dimensionless)

*P*=pressure (Pa)

*P*_{avg}=average value of the pressure fluctuations (Pa)

*P*_{ref}=reference pressure (Pa)

*q*=second invariant of the velocity gradient tensor (1/s

^{2})*Q*=volumetric flow rate (m

^{3}/s)*r*=radial coordinate,

*R*/*R*_{ref}(dimensionless)*R*=radial coordinate (m)

**R**=rotation-rate tensor

*R*=_{i}radius of the draft tube inlet (m)

*R*_{ref}=reference radius (m)

- Re =
Reynolds number (dimensionless)

- RKE =
realizable

*k–ϵ**s*=swirling strength (1/s)

**S**=strain-rate tensor

*t*=time (s)

**u**=velocity vector

*U*=_{R}radial component of velocity (m/s)

*U*=_{Z}axial component of velocity (m/s)

*U*=_{θ}circumferential component of velocity (m/s)

*X*=*X*-coordinate (m)*Y*=*Y*-coordinate (m)- $y\u0302$ =
modified radial coordinate, $Y\u0302/Rref2$ (dimensionless)

- $Y\u0302$ =
modified radial coordinate,

*R*^{2}/2 (m^{2})*Z*=*Z*-coordinate (m)*α*=half-angle of the vortex rope (deg)

*γ*=specific energy coefficient (dimensionless)

- Δ =
swirling discriminant (1/s

^{6})- Δ
*P*= dominant range of pressure fluctuations (Pa)

*ϵ*=ordering parameter (dimensionless)

*ϵ*=_{t}turbulent dissipation rate (m

^{2}/s^{3})*θ*=tangential coordinate (rad)

*λ*_{2}=second eigenvalue of the symmetry square of the velocity gradient tensor (1/s

^{2})*ν*=kinematic viscosity of water (m

^{2}/s)*ξ*=ratio of flow rates,

*Q*/*Q*_{BEP}(dimensionless)*ρ*=density of water (kg/m

^{3})- $\varphi $ =
discharge coefficient (dimensionless)

*ψ*=stream function (m

^{3}/s)*ω*=angular speed of the runner (rad/s)