Swirling Flow Computation at the Trailing Edge of Radial-Axial Hydraulic Turbines

Modern hydraulic turbines require optimized runners within a range of operating points with respect to minimum weighted average draft tube losses and/or flow instabilities. Tractable optimization methodologies must include realistic estimations of the swirling flow exiting the runner and further ingested by the draft tube, prior to runner design. The paper presents a new mathematical model and the associated numerical algorithm for computing the swirling flow at the trailing edge of Francis turbine runner, operated at arbitrary discharge. The general turbomachinery throughflow theory is particularized for an arbitrary hub-to-shroud line in the meridian half-plane and the resulting boundary value problem is solved with the finite element method. The results obtained with the present model are validated against full 3D runner flow computations within a range of discharge value. The mathematical model incorporates the full information for the relative flow direction, as well as the curvatures of the hub-to-shroud line and meridian streamlines, respectively. It is shown that the flow direction can be frozen within a range of operating points in the neighborhood of the best efficiency regime.


Introduction
In the recent years the development and integration of renewable energy resources emphasized the essential role of hydropower to compensate for the stochastic nature or variable energy sources (e.g. wind and solar) thus preserving the electrical grid stability 1 . As a result, the hydraulic turbines are often required to operate from overload to deep part load, while traditionally their design was focused on maximizing the peak efficiency. On the other hand, particularly for the wide-spread Francis turbines of medium-low head, the shape of the hill chart is driven by the steep increase in draft tube hydraulic losses at off-design operating points, [1]. In addition, far from the best efficiency point, the decelerated swirling flow downstream the runner developes various self-induced instabilities with associated pressure pulsations [2]. It becomes quite clear that the turbine runner has to be designed in tandem with the draft tube in order to insure efficient and safe turbine operation within a wide operating range.
Simultaneous shape optimization of both hydraulic turbine runner and draft tube, with the objective of increasing turbine's efficiency within a wide range of operating points, Lyutov et al. [3] can be approached by considering a geometrical parameterization of runner crown, band and blades, as well as of the draft tube. However, the rather large number of parameters (e.g. 37 1 The HYPERBOLE project -HYdropower plants PERformance Figure 1. Meridian half-plane, (z, r), for a Francis turbine, and the intrinsic coordinates (ξ, ζ) on the hub-to-shroud line corresponding to the meridian projection of the blade trailing edge, and (s, n) on the meridian streamline, respectively.
in [3]) combined with three-dimensional turbulent flow computation for several operating points hinder the tractability of such approach. Another approach is to first optimize the swirling flow at the draft tube inlet, then design the runner accordingly. Using a swirling flow parameterization introduced in [4], Galván et al. [5] optimized the axial and circumferential velocity profiles at the draft tube inlet for minimum draft tube hydraulic losses, at a given operating point. It is obvious that this approach can be used for optimization within a range of operating points only if the swirl ingested by the draft tube can be reliably computed at arbitrary operating regimes. Such a model has been introduced by Susan-Resiga et al. [6] and validated against experimental data on a cross-section at the draft tube inlet, in the discharge cone. This model was further used to examine the influence of the velocity field at the draft tube inlet on the draft tube hydraulic losses within an operating range, [7], or to compute the swirling flow downstream a Kaplan runner at two operaing points [8]. A full optimization procedure for the swirling flow at the draft tube inlet is exemplified in [9]. Since the runner blades are to be designed for a given swirling flow at the runner outlet, it is obviously more convenient for design purposes to compute the flow as close as possible to the blades trailing edge. A first step in this direction is taken in [10], where the swirling flow exiting a Francis runner is computed on a conical surface (straight line in the meridian half plane, from crown to band). The next obvious step was to consider an arbitrary hub-to-shroud line close to the blade trailing edge, with the model extension presented in [11] and further examined in [12]. It became clear that although the model for the swirling flow exiting the runner captures the main swirling flow features as the turbine discharge is varying within a range of operating points, some discrepancies with respect to reference data obtained from runner flow simulation suggest that not all significant terms in the mathematical formulation are accounted for.
In this paper, we present the complete and rigorously derived mathematical model for computing the swirling flow exiting the runner of radial-axial hydraulic turbines. Practically, we consider the most general case where the hub-to-shroud line, on which the axisymmetric swirling flow is computed, coincides (or it is in the close neighborhood) to the meridian projection of the runner blade trailing edge, Fig. 1. As a result, when computing the draft tube flow one does not need anymore to consider more or less arbitrary draft tube inlet sections. The draft tube is directly fed with the swirl emerging from the runner. Moreover, the model allows the swirling flow computation for arbitrary operating points, as long as there are no severe detachments or secondary flows in the runner interblade channels.
The axial symmetry assumption is quite reasonable for runners with large number of blades, as it is the case for Francis turbines. Computing axisymmetric swirling flows further downstream into the draft tube provides useful information on the occurence and development of the vortex breakdown, [13,14], with further estimation of the flow unsteadiness (e.g. precessing vortex rope [15,16]) without actually computing the unsteady three-dimensional flow. As a result, the model further detalied and validated in this paper could be quite useful in the early design and optimization stages of a new or refurbished runner, allowing a quite accurate assessment of the draft tube hydrodynamics within a range of operating points before actually designing the runner.

The throughflow turbomachinery theory
The simplest model for turbomachinery flows, which captures most of the swirling flow physics, assumes a perfect fluid (incompressible and inviscid) and a steady axisymmetric flow. While the perfect fluid assumption is reasonably valid for large hydraulic machines, the steady and axialsymmetric flow implies an infinite number of blades, each one infinitely thin. This assumption might be acceptable for Francis turbines, as a first approximation, thanks to the relatively large number of blades. Nevertheless, as crude as it may seem, the above simplified model became the backbone of turbomachinery design and it is known as throughflow model. We briefly revisit it, and then we particularize it for the runner outlet.
Throughout this paper all quantities are considered dimensionless with respect to a reference length (R ref , runner outlet radius) and a reference velocity (ΩR ref , where Ω is the runner angular speed). Due to this choice, the dimensionless radius is r ≡ R/R ref , and the dimensionless transport velocity is (ΩR)/(ΩR ref ) = r.
Let us consider the cylindrical coordinates (z, r, θ) with the corresponding unit vectors ẑ,r,θ . The velocity vector can be expressed as, v =ẑv z +rv r +θv θ = ∇ψ × ∇θ + κ∇θ, where, v m ≡ẑv z +rv r is the meridian projection of the velocity, ψ is the Stokes' streamfunction and κ ≡ rv θ is the circulation function, respectively, while in cylindrical coordinates we have ∇θ =θ/r. The relative velocity, w, with respect to the runner rotating frame, has the meridian projection w m = v m and the circumferential component Within the runner bladed region we consider that the relative flow follows a prescribed streamsurface given by the equation θ − ϕ(z, r) = constant, thus the flow tangency condition can be written as w · ∇ (θ − ϕ(z, r)) = 0, or (∇ψ × ∇θ) · ∇ϕ = κ r 2 − 1. (1a) On the tangent plane to the relative flow streamsurface, one can project the steady Euler equation of motion along the streamline direction. The result is a Bernoulli-like theorem stating that the rothalpy i ≡ p + v 2 /2 − κ remains constant along a streamline, i.e. we have a function i(ψ). This function is given by the flow downstream the guide vanes of the hydraulic turbine, and it remains unaltered throughout the runner. The second projection is on the direction normal to the streamline, but still in the plane tangent to the relative flow streamsurface, leading to the so-called principal equation [17], Since Eqs. (1) involve three functions, one can consider two kind of problems: design, when κ(z, r) is given while ϕ(z, r) and ψ(z, r) are computed, or analysis, when ϕ(z, r) is given while κ(z, r) and ψ(z, r) are computed, where (z, r) are the coordinates in the meridian half-plane, Fig. 1.

Natural coordinates
Since the purpose of this paper is computing the swirling flow at the runner outlet, on a hubto-shroud line in the meridian half-plane, it is convenient to employ natural coordinates on the streamline (s, n), Fig. 1, instead of (z, r). The unit vectorŝ is tangent to the streamline and downstream oriented, while the normal to the streamlinen is chosen such thatŝ ×n =θ. The meridian velocity v m will be expressed using its magnitude v m and angle with respect to the axial direction, φ, defined by tan φ = v r /v z . As a result, one can re-write the principal equation According to the Kutta-Joukowski condition, the runner blade loading vanishes at the trailing edge, thus ∂κ/∂s = 0 and the last term in the left-hand side vanishes. As a matter of fact this remains true further downstream in the bladeless region of the discharge cone, and the principal equation reduces to the Bragg-Hawthorne equation. Moreover, for the first two terms in the l.h.s. of Eq. (2) correspond to the shear, −∂v m /∂n, and curvature, v m ∂φ/∂s, components of the circumferential vorticity ω θ . The left-hand side of Eq. (1a) is simply v m ∂ϕ/∂s.

Swirling flow equation at the runner blade trailing edge
The governing equation for the swirling flow exiting a radial-axial hydroturbine runner is a particular form of Eq. (2) projected on a hub-to-shroud line that corresponds to the meridian projection of the runner blade trailing edge, Fig. 1. Let us start with the flow-tangency condition, Eq. (1a), and define the relative flow angle β at runner outlet as, According to the definition (3) we have β < 0 since w θ is in oposite direction with respect to the transport velocity.
The hub-to-shroud line is given in parametric form as z(ξ) and r(ξ), with ξ the arc-length coordinate running from 0 at the hub to at the shroud, Fig. 1. On this line we consider the normal unit vectorζ, downstream oriented, and the tangent unit vectorξ, such thatζ ×ξ =θ. The angle between the unit vectorζ and the axial directionẑ is denoted by γ, with cos γ = dr/dξ and sin γ = −dz/dξ.
By introducing the angle χ ≡ φ − γ, between the streamline tangent and the normal to the hub-to-shroud line, one can rewrite the principal equation (2) as a second order differential equation for ψ(ξ), with variable coefficients a(ξ), b(ξ) and c(ξ), where the right-hand side term introduces a mild non-linearity.  (4) and (5) are the main result of the present theoretical development. The boundary conditions are ψ(0) = 0 and ψ( ) = q/2, respectively. Once the flow direction is frozen, i.e. the coefficients a(ξ), b(ξ) and c(ξ) do not change, the only information about the operating point comes from the discharge coefficient q ≡ Q πΩR 3 ref . In order to solve numerically the differential equation (4) with essential boundary condition we consider its weak form where weighting functionψ(ξ) is an arbitrary perturbation of the solution ψ(ξ), which satisfies homogeneous essential boundary conditions. The left-hand side is a linear operator, while the right-hand side has to be iteratively updated because of the nonlinear rothalpy function i(ψ).
The weak formulation (6) is the basis of the finite element method. In this paper we have employed cubic Hermite 1D finite elements.

Information required for computing the swirling flow
In order to solve the Eq. (6) one needs the information for computing the coefficients (5), as well as the rothalpy versus the streamfunction i(ψ). This section shows how this information can be inferred using numerical data from the circumferentially averaged flow field at runner outlet.

The hub-to-shroud line
The hub-to-shroud line is parameterically represented as z(ξ), r(ξ), where ξ is the arc-length coordinate running from zero at the hub up to the maximum value at the shroud. Initially, the line is given as a set of points (ξ i , z i , r i ), ordered for increasing ξ. A cubic spline interpolation allows the geometrical manipulation in order to compute the curvature and to specify a set of grid points. The grid points, as the ones shown with circles in Fig. 2, are chosen such that the    area of the revolution surface corresponding to any segment is the same.
The curvature is computed using the rotation-invariant formula: Equation (7) reveals the need for cubic spline representation of the arc-length parameterization z(ξ) and r(ξ), with both first and second order derivatives continuous. The plot from Fig. 3 shows that in our case the hub-to-shroud line is a cubic spline since the curvature is continuous and practically piecewise linear. Moreover, dγ/dξ is positive, thus the angle γ between the axial direction and the normal pointing downstream is increasing monotonically from hub to shroud.  The simplified model considered in this paper assumes that the flow direction at runner outlet is known. Given the velocity components, one can compute the meridian flow direction, tan φ = v r /v z , and the circumferential relative flow direction, tan β = (v θ − r) v 2 z + v 2 r . One can see that the flow direction is not significantly altered by the turbine discharge, as long as the operating points is not too far from the best efficiency one. However, the dashed red lines in both Fig. 4 and Fig. 5, corresponding to the part load operating regime suggest that the flow begins to detach from the runner blade the hypothesis of a "frozen" flow direction is no longer valid. On the other hand, one can use the extent of the operating range with unmodified flow direction as a quality indicator of the runner blade: the wider this discharge range, the better the flow guidance within the interblade channels.

Flow direction
In addition to the flow direction on the hub-to-shroud line, we need to evaluate the meridian streamlines curvature. This is performed by considering streamline segments between the two cross-sections downstream the runner, reprezented as cubic polynomials r(z), as shown in Fig. 6. The streamline curvature will be: where r(z) is represented as piecewise cubic polynomial using the end point values for the function and its first derivative. One can see from Fig. 7 that the numerical data for the streamline curvature dφ/ds are clustered in the neighborhood of a master curve represented as a cubic polynomial. Numerical uncertainties are present near the shroud, most likely due to the small length of the streamline segments in this region, as seen in Fig. 6. As a general conclusion for this subsection, the flow direction as well as the streamline curvature can be well represented with simple analytical expressions, which are valid within a certain range of discharge values.  Dimensionless rothalpy versus streamfunction, i(ψ). Solid lines are reconstructed from the master fit in Fig. 9. The rothalpy i ≡ p + v 2 /2 − κ is constant along a meridian streamline, thus we are looking for the function i(ψ) corresponding to each operating point. The numerical data, shown with circles in Fig. 8, are obtained from circumferentially averaged velocity and pressure fields at runner outlet.

Rothalpy
In order to obtain a simple analytical representation we collapse all data as shown in Fig. 9. First, the abscissa is normalized at each operating point by the streamfunction value at the shroud,ψ ≡ ψ/( q/2 ) . Second, the ordinate for each data set is shifted for zero average. Note that the function i(ψ) is defined up to an arbitrary additive constant, thus shifting the ordinate does not change the value of di/dψ . The clustered data points are then fitted with a cubic polynomial, shown with solid line in Fig. 9. The reversed transformation of this master curve generates the lines for each operating point, shown in Fig. 8. Note that the rothalpy decreases from hub to shroud because the circulation κ varies, while the total pressure is practically constant. This is due to the columnar guide vanes geometry within the region with radial-axial meridian flow, as it is the case for low-medium head Francis turbines.

Numerical results for the velocity components at runner blade trailing edge
The numerical solution of Eq. (6) is obtained with our finite element code which employs 1D Hermite elements. Both the streamfunction ψ values and its derivative dψ/dξ are obtained at the nodal points as part of the solution. Then, the velocity components are recovered using the meridian and circumferential flow direction information from Figs. 4 and 5. The second term in the right-hand-side of Eq. (6) contains a mild non-linearity due to the rothalpy derivative di/dψ . As a result, this term should be updated iteratively, starting with an initial guess for ψ(ξ) which corresponds to a constant discharge velocity from hub to shroud. Three iterations are found to be sufficient for the solution to converge. Figure 10 shows the numerical results obtained with the present model, solid lines, and the circumferentially averaged velocity components obtained from 3D turbulent flow simulation in the runner, circles. One can observe a very good agreement of the model predictions with the numerical data, particularly for the meridian (axial and radial) velocity components. On the other hand, the circumferential velocity component is quite sensitive to the approximation from Fig. 5. Nevertheless, the present results are a significant improvement with respect to our previous attempts [11] [12] which implicitly discarded part of the information used for computing the coefficients (5). One can conclude from Fig. 5 that the hypothesis of frozen flow direction remains valid within a certain range of operating points, as long as there are no severe flow detachments in the runner.
On the other hand, Fig. 11 shows that as we move far away from the best effciency regime, for example at part load, the viscous effects lead to an alteration of the relative flow direction and as a result the model can capture only qualitatively the velocity components profiles. Nevertheless, in spite of this drawback, one can use the present model for a-priory assessment of various configurations of the runner blades in the neighborhood of the trailing edge. We conjecture that as the range of the "frozen flow direction" hypothesis increases the runner blade desing improves leading to a better flow guidance within the interblade channel.

Conclusions
The paper introduces a rigorously derived mathematical model for computing the swirling flow at the trailing edge of runner blades when the turbine is operated at various discharge values.
The model accounts for the meridian and tangential flow direction, the curvatures of the hubto-shroud line, the curvature of the meridian streamline, and the variation of the rothalpy from hub to shroud. As a result, one needs a full 3D runner flow simulation at the design operating point in order to assess the above information. Then, the swirling flow exiting the runner at various discharge values can be quickly and robustly obtained with the present model, in order to assess the draft tube performance (hydraulic losses and pressure recovery) within a range of operating points. Moreover, this model is well suited for optimizing the runner blade outlet geometry by significantly reducing the computational effort when considering the optimization within a range of operating points.    Figure 10. Axial, radial and circumferential velocity profiles on the hub-to-shroud line near the runner blade trailing edge, for operating points in the neighborhood of the best efficiency regime. Data from the runner 3D numerical simulation, circumferentially averaged, are shown with circles while the solid lines are obtained with the present model.  Figure 11. Velocity components at part-load regime, q = 0.237993 (guide vane opening 20 • ).