Proper Orthogonal Decomposition of Pressure Fields in a Draft Tube Cone of the Francis (Tokke) Turbine Model

The simulations of high head Francis turbine model (Tokke) are performed for three operating conditions - Part Load, Best Efficiency Point (BEP) and Full Load using software Ansys Fluent R15 and alternatively OpenFOAM 2.2.2. For both solvers the simulations employ Realizable k-e turbulence model. The unsteady pressure pulsations of pressure signal from two monitoring points situated in the draft tube cone and one behind the guide vanes are evaluated for all three operating conditions in order to compare frequencies and amplitudes with the experimental results. The computed velocity fields are compared with the experimental ones using LDA measurements in two locations situated in the draft tube cone. The proper orthogonal decomposition (POD) is applied on a longitudinal slice through the draft tube cone. The unsteady static pressure fields are decomposed and a spatio-temporal behavior of modes is correlated with amplitude-frequency results obtained from the pressure signal in monitoring points. The main application of POD is to describe which modes are related to an interaction between rotor (turbine runner) and stator (spiral casing and guide vanes) and cause dynamic flow behavior in the draft tube. The numerically computed efficiency is correlated with the experimental one in order to verify the simulation accuracy.


Introduction
The numerical simulations of hydraulic turbines are recently common for industrial companies. The main goal is to predict hydraulic efficiency of the machine for an entire range of operating conditions and thus to verify turbine hydraulic design. Nevertheless in the industry there are not usually enough computational resources and time to deal with complex calculations to verify turbulence models, discretisation schemes and perform time-consuming post-processing on large computational grids. Therefore the Francis-99 test case is a good option to deal with the post-processing methods e.g. proper orthogonal decomposition (POD) and compare it with the experimental results of the model Francis turbine [1]. POD offers complete spatio-temporal description of unsteady flows and it is especially suited to flows with strong coherent structures. The model turbine of specific speed

Computational study
In order to perform numerical calculations, both commercial and open source software were chosen with opportunity to test code performance. The turbine geometry and meshes are prepared in ICEM software and provided by the workshop organizers through the website [2]. Only small necessary corrections were made to enhance quality of the interfaces between the spiral case and turbine runner. The meshes contain two stationary (spiral casing and draft tube) and one rotational (runner) parts. The mesh for each operating point includes approximately 12.5 million of elements. All computations ran employing linux cluster partitioned into 16 cores in case of Ansys Fluent and into 64 cores in case of OpenFOAM software. The unsteady solutions were performed for all three operating points until the periodic flow state was achieved, then the time averaged flow fields and pressure monitors were recorded. The geometrical view of computational domain is presented in figure 2.

Computation using Ansys Fluent R15
Using commercial software brings relatively easy case set-up due to proved setting options and graphic user interface. On the other hand user is restricted by a number of available licences and inaccessible code modifications. The Ansys Fluent was chosen because it is not as widely spread in turbo-machinery computations as for example CFX and also due to authors' previous experiences.
The case set-up of Ansys Fluent calculation presented in table 1 is identical for all three operating regimes. The material properties and operating conditions (density, kinematic viscosity, runner speed and flow rate) were defined according to the measurements [2]. The time-step size was chosen to correspond with 1° of the runner rotation.

Computation using OpenFOAM
The open source software OpenFOAM has become arguably the best option for the license-free calculations. We have used OpenFOAM in version 2.2.2 and incompressible pimpleDyMFoam solver with implemented AMI interfaces between the mesh parts. The calculations are carried out using identical meshes as in case of Ansys Fluent which were transformed to the OpenFOAM format via fluent3DMeshToFoam utility.
The case set-up for calculations using OpenFOAM is similar with Fluent except particular differences in the spatio-temporal discretization, see table 2. At the outlet of the domain the value of pressure boundary condition was set to zero considering that the OpenFOAM use kinematic pressure in form .

Velocity profiles in the draft tube
The experimental velocity profiles measured in two lines situated downstream in the draft tube cone were provided in order to verify numerically computed flow exiting the turbine runner [2]. The locations of LDA measuring lines are shown in figure 3. The locations are 64mm (LDA-1) and 382mm (LDA-2) below the draft tube inlet. On the other hand tangential velocity profiles computed using OpenFOAM have better agreement with the experimental one than in case of Ansys Fluent calculation. Analysis shows that axial velocity profiles are much better predicted than the tangential ones for both Fluent and OpenFOAM. While in the wall region the inaccuracy between measurements and computations is probably caused by the near wall modelling, the significant variation in the central region (especially at BEP and FullLoad) is attributed to both inaccuracy of LDA measurements (see figure 6b, where the measured tangential velocity in the axis suppose to be zero) and inaccuracy in turbulence models (large discrepancy between Fluent and OpenFOAM). It has to be mentioned, that slight difference between measured and computed profiles at BEP and Full Load can arise due to different absolute values of the head, flow rate and runner speed used for the measurements in order to suppress vibrations induced in the test rig [2].

Pressure pulsations
The numerical unsteady pressure pulsations were evaluated in a three different probe locations. The first probe (VL01) is situated in the gap between the spiral distributor and the turbine runner, the second probe (DT11) and the third probe (DT21) are situated in the cone of the draft tube. The probe locations correspond with the locations of pressure transducers during the experiment, see figure 7.  First, the average values of computed pressure are estimated and compared with the experiment, see table 3. The numerical results are presented in a form of percentage deviation from the experimental measurements. The largest discrepancy between computed and measured pressure was found at High Load for DT11 probe where the value of numerical pressure underestimates the experimental one by 7.7 % (Fluent) and 8.3 % (OpneFOAM). The best agreement between the computed and measured pressure is at BEP for OpenFOAM where the numerically computed pressure in DT21 probe underestimate the experimental one just by 0.3%. The amplitude-frequency spectra are carried out using general FFT function in MATLAB software for all three operating points. The frequencies are normalized assuming the runner angular speed, see table 4. The runner has 15 full length and 15 splitter blades, thus the normalized frequency, corresponding with the blade passing frequency is = 30 with its higher harmonics = 60, 90, … Due to influence of splitter blades the normalized frequency = 15 also appears.
The frequency spectra of DT11 and VL01 probes are presented in figures 8, 10 and 12 and figures 9, 11 and 13 respectively. The spectra of measured signals can be compared with the spectra of computed signal of both Ansys Fluent and OpenFOAM (please note different scales of vertical axes in all plots). The best agreement for VL01 probe was found at High Load (figure 13) where the amplitude magnitudes and normalized frequencies of both numerical spectra well correspond with the measured one. On the other hand, the largest disagreements in amplitude magnitudes were found for DT11 probe (figures 8b, 10b and 12b) in case of numerical simulations using Ansys Fluent where the high frequency noise can be observed. It is clear that the source of this noise arises from the numerical calculation and we think that the reason can be influence of the mesh interfaces interaction produced downstream to the draft tube.
Furthermore, the particular discrepancies in the amplitude magnitudes can arise from the different lengths of pressure signal and different sampling frequencies used for the measurements and computations.
The and is pressure and friction torque respectively, n is the runner angular speed, Q is the flow rate and ∆ and ∆ is differential static pressure and differential dynamic pressure respectively. In case of OpenFOAM the integral values of pressure and velocity were calculated using patchAverage utility on corresponding patches. The best agreement between computed and measured hydraulic efficiency is at BEP. While the Fluent calculated hydraulic efficiency with 0.64% error compare to the experiment the OpenFOAM was more accurate with error only 0.22%. However the significant discrepancies were obtained for off-design conditions at Part Load and High Load. At Part Load the Fluent calculated hydraulic efficiency with 11.87% error compare to the experiment and at High Load with 1.75% error. Moreover the hydraulic efficiency calculated using OpenFOAM shows irrational results of hydraulic efficiency = 99.5% for Part Load and 100.3% for High Load. At High Load the head and differential pressure computed using OpenFOAM are in much better agreement with the experiment than the Fluent calculation. Nevertheless the runner torque is considerably overestimated which is here the main reason of nonsensical hydraulic efficiency. At Part Load the OpenFOAM wrongly computes both the runner torque and net head which is influenced by the inaccuracy of tangential velocity profile at the runner outlet, see figure 4. Considering the Euler equation ( = − ) the computed head respectively specific energy is higher than in case of measurements.
Due to discrepancy in results at Part Load and High Load the PRESTO! scheme for pressure discretization was tested for Ansys Fluent calculations. The PRESTO! scheme is generally

Proper orthogonal decomposition
The proper orthogonal decomposition (POD) was introduced into the field of fluid mechanics by Lumley [4] as a technique which enables to bridge time and frequency domain and obtains spatiotemporal information about the flow in a form of spatio-temporal eigenfunctions. Comprehensive theoretical description of POD, including classical approach and its properties, can be found in [5]. The POD provides a basis for the modal decomposition of a set of functions and most efficient way of capturing the dominant components, i.e. coherent structures. The POD can be used to analyze experimental as well as numerical data being applied to scalar or vectorial functions. The scalar function ( ; ), e.g. pressure field, is considered in following theoretical explanation. For instance, a set of data ( , ) as a function of physical space and time t, can be expressed using POD as a set of orthogonal spatial basis functions ( ) (i.e. spatial modes), and temporal functions ( ) (i.e. temporal modes), respectively. Where = 1,2, … , , N is the number of grid points and = 1,2, … , , M the number of snapshots. Accordingly, the approximation of the data set onto the first k snapshots can be written in terms of the spatial and temporal functions as follow where < has the largest mean square projection [6].
In the method of snapshots suggested by Sirovich [7], the general N x N eigenvalue problem is reduced to M x M eigenvalue problem. Especially, this solution brings a substantial reduction of computational effort if the number of grid points N significantly exceeds number of data set snapshots. Numerical computation is an appropriate example to use snapshot method. The inner product between every pair of pressure fields, i.e. snapshots, is the temporal correlation matrix ( , ) used as the kernel and defined by eq.  The spatial modes ϕ (x) are calculated based on orthogonal property: Particularly, the flow structure is fully described using spatio-temporal representation of modes ( ), ( ). The correlation matrix ( , ) is symmetric and positive definite. Thus the eigenvalues are real and define the contribution of the coherent structure into complex flow [8]. Correspondingly, the eigenvalue magnitude of one single mode represents an amount of this contribution from the total sum ∫ 〈p (x)p (x)〉dx Ω [7]. Finally, the computed modes are sorted out taking into account the largest value of the eigenvalue magnitude in order to identify the most significant structures. As a result, the temporal modes ( ) are calculated by projecting the data function ( ; ) on the eigenfunction ( ).
Throughout, the decomposed function ( ; ) is considered in dimensionless form in order to compare the values associated to different investigated cases. The dimensionless form of the pressure field is done as follows where ∆ ( , ) = ( , ) − , ( , ) is the static pressure field obtained in numerical simulation at each time step t, p static pressure mean value on the outlet surface of the domain (draft tube outlet) computed as a mass-weighted surface integral in form: and corresponds to bulk velocity at the draft tube inlet.

POD of pressure fields
The successful application of the POD for investigation of the vortex rope dynamic in the draft tube of the hydraulic turbine can be found in [9] and [10]. The previous investigations were carried out for a constrained domain of the simplified draft tube. In this paper the entire turbine domain ensure full dynamical influences, primarily the rotor-stator interaction. The POD analysis is applied on the numerical data from the Ansys Fluent calculations for all three operating points. The node values of the static pressure are extracted in the longitudinal slice along the draft tube cone, starts180 mm and ends 638 mm below the dividing plane of the spiral casing. The pressure data are stored with the specified number of constant time-step intervals. Then the POD is executed using script written in Matlab software where the eigevalues and eigenvectors are calculated employing the eig function. The resulting spatio-temporal representations of the static pressure modes can be correlated with the dynamic behaviour extracted from the pressure monitoring points. Thus we are able to decide which mode is consequence of the rotor-stator dynamic and which is caused by the self-induced instability in the draft tube. It is necessary to say that when the significant coherent structure (e.g. vortex rope) does not appear in the draft tube (as is case of these operating points), the resulting modes have very small eigenvalue magnitude compare to the mean flow, see table 5. In the following subsections the spatio-temporal representations of the first five modes for each operating point are presented and their relation to the turbine regime is discussed. The value of static pressure is normalized according to the equation (6). For all three operating points no significant vortex rope leading to the strong modes of self-induced instability appears [9]. Therefore in the present analysis we will distinguish two kinds of the static pressure modes. The modes related to the forced instability with the frequency at ratio of the runner rotational frequency and modes related to the rotorstator interaction at the blade passing frequency.    (figure 18c, d) with normalized frequency = 2 and the fifth mode (figure 18f) with = 0.2 are consequence of the forced instability forming unstable swirl rotating in the draft tube cone two times faster (the second and third mode) respectively five times slower (the fifth mode) than the turbine runner. The fourth mode (figure 18e) has dominant normalized frequencies = 15 and = 30 and spatial representation shows correlation to the rotor-stator interaction combining full length blades with the splitter blades.

Result discussion
At the BEP operating point one can see that the computed head using both Ansys Fluent and OpenFOAM is considerably overestimated. The cause lays in a wrongly computed tangential velocity profiles at the turbine outlet where the velocity magnitudes overestimate experimental measurements. Nevertheless the computed runner torques are also overestimated, thus the resulting numerically computed hydraulic efficiencies adequately correspond with the experimental ones.
Different situation is for off-design conditions. While the Fluent is able to compute hydraulic efficiency with a relatively small error at High Load and with larger error at Part Load, the results obtained using the OpenFOAM are not credible. The main reason of OpenFOAM failure was found in a computation of the runner torque which is for both Part Load and High Load significantly overestimated. This supposes to be further investigation for the future reliable usage of the OpenFOAM solver. The main difference between Fluent and OpenFOAM can be found in the different numerical schemes used in each solver and in the estimation of the runner torque. The unsteady pressure analysis shows large disagreement between computed and measured amplitudes in a case of DT11 probes. The spectra of numerical signal from Ansys Fluent include the high frequency noise which is not presented in the experimental results. We suppose that this noise is a main cause of the amplitudes disagreement. Nevertheless [2] informs about vibrations induced in the test rig at BEP and High Load, which interfered with the LDA measurements at these conditions. Regarding to the pressure measurements we think that the high amplitudes at the measured frequency = 54 (BEP) and = 49 (High Load) in DT11 spectra could be potentially cause of these vibrations as well.
From the POD analysis at Part Load and High Load operating points it was found that some temporal modes related to the forced instability have more than only one dominant frequency of the vortex rotation, the normalized frequencies related to the rotor-stator interaction are also visible. That means the rotor-stator interaction is significant source of the pressure pulsations in the draft tube. Moreover in a case of the rotor-stator interaction some modes with normalized frequency = 15 appear. These modes are influence of fifteen splitter blades. From the eigenvalue magnitude presented in table 4 it is clear that the Part Load is the most dynamic operating regime. Nevertheless when the strong coherent structure, in a form of the vortex rope, does not appears as it is case of these investigated turbine regimes, the dynamic behaviour of the pressure field in the draft tube is very low as can be seen in the percentage amount of the eigenvalue presented in table 5. The most dominant pressure modes were found modes related to the rotor-stator interaction and the modes caused by the forced instability in a relation to the runner rotation.
Due to the large discrepancy in the OpenFOAM results at Part Load and High Load, different numerical schemes were tested for the disretization of the time (backward) and the vector fields (Guass limitedLinear). It was found that the numerical schemes have large influence on the results accuracy especially on the runner torque calculation. Nonetheless using these schemes did not lead to more realistic value of the hydraulic efficiency at Part Load and High Load respectively.