The Complex Dynamics of the Precessing Vortex Rope in a Straight Diffuser

The decelerated swirling flow in the discharge cone of Francis turbines operated at partial discharge develops a self-induced instability with a precessing helical vortex (vortex rope). In an axisymmetric geometry, this phenomenon is expected to generate asynchronous pressure fluctuations as a result of the precessing motion. However, numerical and experimental data indicate that synchronous (plunging) fluctuations, with a frequency lower than the precessing frequency, also develops as a result of helical vortex filament dynamics. This paper presents a quantitative approach to describe the precessing vortex rope by properly fitting a three-dimensional logarithmic spiral model with the vortex filament computed from the velocity gradient tensor. We show that the slope coefficient of either curvature or torsion radii of the helix is a good indicator for the vortex rope dynamics, and it supports the stretching - breaking up - bouncing back mechanism that may explain the plunging oscillations.


Introduction
When the hydraulic turbines are operated far from the best efficiency regime the residual swirl exiting the runner, and further decelerated in the discharge cone, leads to self-induced instabilities with associated pressure fluctuations that hinder the turbine operation. Dörfler et al. [1,Ch.2] review the physical mechanisms that produce such, sometimes severe, flow pulsations although the flow at runner outlet can be considered quasi-axisymmetric. The analysis of the wall pressure fluctuation [2] allows the representation of the pressure signal as a superposition of synchronous (plunging), asynchronous (rotating) and random components, with relative importance varying with the discharge. In particular, the operating regime at 70% from the best efficiency discharge displays strong asynchronous pulsations generated by the precessing helical vortex, also known as "vortex rope". However, a synchronous component is also present. The origin of this synchronous pressure fluctuation is usually related to the interaction between the precessing vortex rope and the draft tube elbow, and as a result Dörfler et al. [1, §2.2.11.2] advocate that such fluctuations should not be present in straight, axisymmetrical diffusers.
The precessing vortex rope in the discharge cone of hydraulic turbines was extensively investigated both experimentally and numerically [3]. The numerical experiment considered in this paper is actually reproducing the flow phenomenon studied in [3], and further explored when the vortex core is cavitating [4]. However, instead of using the full three-dimensional draft tube geometry we consider a simplified axisymmetric geometry with variable wall radius equal to the equivalent hydraulic radius of the real draft tube [5], as shown in Fig. 1. When considering an axisymmetric flow model [5] we obtain a steady flow field displaying a central stagnant region separated from the main swirling flow by a vortex sheet. We showed that the measured vortex rope is actually wrapped on this vortex sheet. This result was further exploited to find the main vortex rope parameters (precessing speed, helix geometry, wall pressure pulsation) without actually computing the unsteady 3D flow [6] [7]. The present paper is focused on describing the geometry of the helical vortex filament by considering a three-dimensional logarithmic spiral model. Such model was used in [8][9] to represent a segment of the vortex with the vortex core determined experimentally with 3D-PIV. A more rigorous approach is developed in this paper using the arc-length parameterization of the helix and a leastsquares fitting procedure to determine the curve parameters. The particular 3D spiral considered in this paper to approximate the vortex filament has the property that both curvature and torsion radii increase linearly with the arc-length, while their ratio remains constant. More general 3D spirals with independent linear curvature and torsion radii could also be considered as in [10]. The paper is organized as follows. Section 2 briefly describes the numerical experiment and the vortex filament identification, as detailed in [11]. Section 3 is devoted to the three-dimensional logarithmic spiral model and the least-squares fitting procedure. Section 4 examines the vortex rope dynamics using its geometrical parameters, and the last section summarizes the conclusions.

Numerical experiment and vortex filament computation
The present paper develops a general methodology for quantifying the shape of the precessing helical vortex and uses the numerical experiment detailed in [11], further summarized in this section. The computational domain corresponds to the surrogate draft tube shown in Fig. 1 right, which includes the discharge cone of an actual draft tube followed by a variable radius pipe. The radius corresponding to the elbow region is considered to be the equivalent hydraulic radius of the actual non-circular crosssections shown in Fig. 1 left [5]. The inlet surface corresponds to the outlet of the Francis runner, [11,Fig. 2], and an average constant pressure condition is considered on the outlet section. A block structured hexahedral grid of 1.5 million cells is built for the axisymmetric flow domain.
Although the computational domain is axisymmetric, and the inlet boundary condition corresponds to a steady axisymmetric velocity field, the decelerated swirling flow develops a self-induced instability and becomes unsteady and three-dimensional. The FLUENT 15 expert CFD software with Scale-Adaptive Simulation (SAS) method for turbulence modeling is employed for the numerical flow simulation. The SAS model, introduced in [12], has the useful feature that for unstable flows it changes smoothly from a Large Eddy Simulation (LES) through various stages of eddy-resolution back to a Reynolds-Averaged Navier-Stokes (RANS) model based on the specified time step. Menter and Egorov claim that SAS is an advanced unsteady RANS model which can produce the spectral content for unstable flows. This claim is supported by a series of examples presented in [13]. The FLUENT setup uses SIMPLEC for pressure-velocity coupling method, Least Squares Cell Based approximation for gradients, second-order for pressure, bounded central differencing for momentum, second order upwind for turbulent kinetic energy and specific dissipation rate, respectively. The transient formulation was set to bounded second order implicit. A time step of 0.167 ms was considered for a total flow time more than 10 seconds. The results analyzed in this paper correspond to the last two seconds of this interval.
Our particular interest is the analysis of the precessing helical vortex, which is a characteristic feature of the decelerated swirling flows downstream hydraulic turbine runners operated at part-load. Identifying the vortex filament is not a trivial task. Although it is related to the velocity field, most of the time the so-called vortex rope is visualized as a constant pressure surface. This approach has the main drawback that in decelerated flows the average pressure level rises downstream and the constant pressure surface approximately surrounds an upstream segment of the vortex filament. Other approaches directly address the velocity field. When considering the classical splitting of the velocity gradient tensor in symmetrical (strain rate tensor) and anti-symmetrical (with associated vorticity vector) components, one can consider the ratio of these two tensor magnitudes to build the Q-number indicator. Such approaches rely on an arbitrary value for iso-surface visualizations.
A more rigorous approach for vortex visualizations considers the velocity gradient tensor ∇v computed on each grid cell. The eigenvalues of this tensor λ 1 ,λ 2 ,λ 3 ( ) are the fundamental quantities which determine the flow pattern. When we have a real eigenvalue λ r and a complex conjugate pair, the flow is spiraling around the direction associated with the eigenvector associated with λ r . The intersection of the support line of this vector with the cell faces generates a vortex filament segment. This method, introduced in [14] is implemented by the routine FX_VORTEXCORE from the Fluid Feature Extraction Tool-kit [15] which returns a set of vortex core segments with associated vortex core strength (vorticity magnitude) values. Although it has a solid theoretical foundation, due to various numerical approximations this approach does not produces contiguous lines, sometimes it locates features that are not vortices, and is sensitive to other non-local vector field features [14].   Figure 2 shows several samples of vortex filaments identified with the above algorithm implemented in the TECPLOT software. One can easily identify the characteristic helical shape, with the corresponding precessing motion. However, we noticed in [11] that besides rotating this helical vortex is also changing its shape in a quasi-periodical low-frequency cycle. This paper attempts to quantify this evolution, as we conjectured that this might be the cause of the associate plunging pulsations.

Three-dimensional logarithmic spiral model
The 3D logarithmic spiral was shown in [8][9] to adequately represent the experimental data for the vortex rope filament measured with 3D Particle Image Velocimetry. Instead of using the polar angle as parameter, we consider the 3D logarithmic spiral parameterized by the arc-length [16] [17], x(s), y(s), z(s) as: In Eqs. (1) the shifted arc-length coordinate s * is considered with respect to the cone apex. The coefficient c 0 is the phase shift, corresponding to the precession angle of the helix. Note that both x(s) and r(s) are linear functions of the arc-length coordinate.
Since we consider an arc-length parameterization of the curve, the coefficients must satisfy the following constraint: As a result the model (1) has only five independent parameters, a 0 ,a 1 ,b 0 ,b 1 , and c 0 , which provide the required flexibility to fit the numerical data. Let us summarize geometrical meaning of the parameters describing the three-dimensional helix: • b 0 b 1 gives the location of the cone apex corresponding to s * = 0 ; • b 1 a 1 gives the cone half-angle γ , i.e. tanγ = b 1 a 1 ; • a 0 accounts for the arbitrary origin of the axial coordinate; • c 0 corresponds to the phase, i.e. the precession angle; • c 1 is related to the helix pitch.
Using the general formula for the curvature κ of a smooth 3D curve we obtain the curvature radius, R c , for the 3D logarithmic spiral (1) as, Using the general formula for the torsion τ of a smooth 3D curve we obtain the torsion radius R t , for the 3D logarithmic spiral (1) as, We note from Eqs. (3) and (4) that the ratio of curvature and torsion radii is R c R t = a 1 An alternative to this coupled fitting procedure would be to separately fit x(s) = a 0 + a 1 s and r(s) = b 0 + b 1 s , then find the precessing angle c 0 . However, we found that the segregated fit is useful to find a first guess for the fit parameters, then one can proceed with the rigorous coupled approach as in Eq. (5). The numerical implementation uses the IMSL subroutine UNLSF which solves a nonlinear least-squares problem using a modified Levenberg-Marquardt algorithm and a finite-difference Jacobian [18].

Numerical results and analysis
As mentioned in §2, the algorithm based on the velocity gradient tensor analysis provides a vortex segment for each cell where the set of three eigenvalues contains a real values and a pair of complex conjugates. As a result, the following post-processing procedure is applied in order to obtain a helical vortex filament. First, we select the vortex segments with the vortex core strength within a proper interval. This step removes spurious segments which are not actually related to the helical vortex filament. Second, we order the points for increasing axial coordinate. Third, we compute the arclength coordinate by successively adding the length of each cell segment, thus obtaining the numerical data s i , x i , y i , z i ( ) , shown with circles in Fig. 3 left. Figure 3 shows a typical fit for the arc-length parameterization of the vortex filament, using the approach from §3. One can see that the reconstructed three-dimensional logarithmic spiral (green line) closely matches the numerical vortex filament (black line).  The fit parameters obtained for the helical vortex filament at various instants in time are shown in Tab. 1. Examples of the resulting fits are further shown in Fig. 4, to illustrate the robustness of the approach introduced in this paper. One can see that even when there is a gap in the numerical data, our leastsquares approach produced good results. The analysis of the precession angle c 0 is shown in Fig. 5, where the values from Tab. 1 have been corrected with multiples of 2π in order to obtain a monotonically increasing function versus time.
The slope of the linear regression (solid red line in Fig. 5) provides the average angular speed, with a corresponding precessing frequency of 3.582 Hz. The corresponding experimental value is 3.75 Hz [3], thus we can say that our numerical experiment closely match the measurements with respect to this parameter. However, one can notice in Fig. 5 that the precessing angular speed fluctuates around the average value, indicating that in spite of the simple problem setup the rope has a more complicated motion.   Fig. 6 reveals low frequency changes in both curvature and torsion slopes. As a result, one can imagine the helical vortex filament as a spring that is successively elongating and compressing, thus producing the low frequency plunging fluctuations. We conjecture that this behavior is related to a cycle of vortex stretchingbreaking up -bouncing back which is going to be further investigated.

Conclusions
The paper introduces a rigorous approach to quantify the shape of the precessing helical vortex and to describe its complex dynamics. First, we employ a method for vortex filament extraction from a three-dimensional velocity field. This method has a solid theoretical foundation and it is based on the eigenvalue analysis of the velocity gradient tensor. A further post-processing of cell vortex segment information provides a numerical representation of the vortex filament as a collection of points s i , x i , y i , z i ( ) ordered for increasing arc-length coordinate s .
Second, we consider the three-dimensional logarithmic spiral model, with the arc-length natural parameterization, and we show that a set of five parameters provides the required flexibility to represent the numerical data. These parameters are determined at successive instants in time by a coupled least squares fit that simultaneously considers the functions x(s), y(s), and z(s) . It is shown that the reconstructed analytical spiral closely matches the vortex filament found from the numerical simulation.  Third, from the variation in time of the precessing angle we determine an average precessing frequency in good agreement with the measured value.
Fourth, by observing that both curvature and torsion radii increase linearly with the shifted arclength coordinate (measured from the cone apex), we consider the corresponding proportionality coefficients as the relevant parameters to describe the dynamics of the helix shape. Low values of these coefficients indicate an elongated spiral, while large values would correspond to a compressed spiral, thus prompting the mechanical analogy with a spring that is oscillating axially with a frequency lower than the precessing one. We conjecture that these axial pulsations generate the plunging fluctuations even in axisymmetric geometries, and correspond to an intrinsic elasticity of the inherently unstable decelerated swirling flow.