Design of Experiments Applied to Francis Turbine Draft Tube to Minimize Pressure Pulsations and Energy Losses in Off-Design Conditions

This paper proposes an original approach to investigate the influence of the geometry of Francis turbines draft tube on pressure fluctuations and energy losses in off-design conditions. It is based on Design of Experiments (DOE) of the draft tube geometry and steady/unsteady Computational Fluid Dynamics (CFD) simulations of the draft tube internal flow. The test case is a Francis turbine unit of specific speed Ns=120 m-kW which is required to operate continuously in off-design conditions, either with 45% (part-load) or 110% (full-load) of the design flow rate. Nine different draft tube geometries featuring a different set of geometrical parameters are first defined by an orthogonal array-based DOE approach. For each of them, unsteady and steady CFD simulations of the internal flow from guide vane to draft tube outlet are performed at part-load and full-load conditions, respectively. The influence of each geometrical parameter on both the flow instability and resulting pressure pulsations, as well as on energy losses in the draft tube, are investigated by applying an Analysis of Means (ANOM) to the numerical results. The whole methodology enables the identification of a set of geometrical parameters minimizing the pressure fluctuations occurring in part-load conditions as well as the energy losses in both full-load and part-load conditions while maintaining the requested pressure recovery. Finally, the results of the CFD simulations with the final draft tube geometry are compared with the results estimated by the ANOM, which demonstrates that the proposed methodology also enables a rough preliminary estimation of the draft tube losses and pressure fluctuations amplitude.


Introduction
Francis turbine, which roughly represents 60% of the installed hydropower capacity worldwide, is a mature technology featuring high efficiency up to 95% at its design point. In addition, Francis turbines are capable of load variations from no-load to full-load by adjusting the opening angle of the guide vanes. Francis turbine hydropower stations are traditionally used as base-load generating units with operations limited within a narrow range of operating conditions around their design point. In the recent years, the increasing electricity production by intermittent renewable energy sources, such as wind and solar, has increased the need of the electrical power system for flexibility. In this context, hydropower stations, including Francis turbine units, can play a crucial role by providing ancillary services to the grid and balancing the electricity production and consumption [1]. This requires an increasing number of start and stop sequences, load variations, as well as operations in off-design conditions to fulfill with the electrical power system requirements.
Francis turbines operating in off-design conditions are subject to the development of flow instabilities inducing vibrations and additional fatigue of the rotating and stationary components [2]. In the draft tube, the residual swirl leaving the runner gives rise to the development of the so-called cavitation vortex rope in the draft tube, whose characteristics depend on the operating conditions of the machine. At part-load (with a flow rate lower than the design value), the high swirling flow leads to the formation of a precessing vortex core (PVC) featuring a helical structure and a periodic motion of precession around the draft tube center axis with a frequency between 0.2 and 0.4 times the runner frequency [3][4][5]. It induces pressure fluctuations at the precession frequency in the draft tube and the upstream waterways. Cavitation may appear within the vortex core if the pressure in the draft tube is sufficiently low, which slightly affects the dynamics of the PVC and the resulting pressure fluctuations. In cavitation conditions at full-load (with a flow rate higher than the design value), cavitation surge may occur in the draft tube, characterized by self-excited oscillations of the cavitation volume and severe pressure pulsations [6,7]. An abundant literature focusing on the physical mechanisms involved in Francis turbine draft tube flow instabilities has been produced in the past few decades, including experimental, numerical and analytical investigations. Such off-design flow instabilities induce pressure pulsations in the entire hydraulic circuit, potential resonances, destabilizing power swings and additional vibrations, thus limiting the operating range of hydropower units and their capability to provide flexibility to the electrical power system. This also causes additional dynamic loadings on the machine components, increasing the fatigue, the risk of crack and finally the maintenance costs for hydropower utilities [8][9][10][11].
The prediction of Francis turbine draft tube flow in off-design conditions and their interaction with the surrounding system remains challenging despite the recent advancements in 3D Computational Fluid Dynamics (CFD) simulations [12,13], 1D hydro-acoustic modelling [14,15] or 1D-3D coupled simulations [16]. In addition to predictive numerical methods, methods of active and passive flow control have been developed and optimized to mitigate the development of the precessing vortex core at part-load conditions and the onset of cavitation surge at full load conditions-see, for instance, [17,18] for a review. Active flow control based on water and air injection has been extensively studied by both numerical and experimental approaches [19][20][21].
The draft tube of Francis turbines is a key component, transforming the remaining kinetic energy into pressure, thus increasing the available head for the turbine. In addition, it is subject to the development of the flow instabilities described previously, directly impacting the operations and safety of the machine in off-design conditions. Nevertheless, the influence of the draft tube geometry on the onset and intensity of off-design flow instabilities has not yet been studied in details and is usually not considered in the design of the machine but it may be crucial to ensure smoother operations in off-design conditions. The components of Francis turbine units are generally designed by performing steady CFD simulations coupled with optimization algorithms to reach the best efficiency at the design point or within a wider operating range [22,23]. However, the unsteady nature of the flow in the draft tube in off-design conditions requires unsteady CFD simulations to assess the dynamic behaviour of the draft tube flow and the resulting pressure fluctuations, making the optimization of the draft tube geometry by considering operations in off-design conditions time-consuming and computationally expensive and, therefore, unfeasible in standard hydropower projects.
In this context, this paper proposes an original approach to assess the impact of the draft tube geometry on the development of flow instabilities and energy losses in off-design conditions. It finally enables the determination of an ideal combination of geometrical parameters reducing pressure fluctuations and energy losses. The method is based on a state-of-the-art Design of Experiments (DOE) approach applied to the geometry of the draft tube and the results of unsteady and steady CFD simulations of the draft tube internal flow at part-load and full-load conditions, respectively. The test-case is a Francis turbine unit with a specific speed Ns = 120 m-kW located in Japan. The power station is required to operate with a constant flow rate of 0.6 m 3 s −1 during 6 months of the year (part-load conditions) and 1.5 m 3 s −1 the rest of the year (full-load conditions), and is, therefore, operating continuously in off-design conditions. A set of nine draft tubes' geometry is first determined by applying a DOE orthogonal array approach to the draft tube geometry that is reduced to a set of four geometrical parameters. Steady and unsteady CFD simulations of the turbine internal flow from guide vane inlet to draft tube outlet are performed for each draft tube geometry at full-load and part-load conditions, respectively. Based on the results of the complete set of CFD simulations, the influence of each geometrical parameter on both the energy losses and pressure fluctuations is then assessed by an Analysis of Means (ANOM). This finally enables the determination of a draft tube geometry minimizing both energy losses and pressure fluctuations in off-design conditions while maintaining sufficient pressure recovery to achieve the desired effective head.

Test-Case and Motivations
The test-case is a hydropower station located in Japan and composed of one Francis turbine unit with a rated available head (difference of elevation between the free-surface of the headwater and tailwater reservoirs) of approximately H = 125 m. The specificity of this power station is that it is required to operate with a constant flow rate of Q = 0.6 m 3 s −1 during 6 months of the year (typically from May to September-summer flow rate) and with a higher constant flow rate Q = 1.5 m 3 s −1 during the rest of the year (typically from October to April-winter flow rate) due to restrictions from the downstream water network. Therefore, two different runners are used, each one designed for a specific flow rate, and are exchanged every six months to fit with the required flow rate.
In the framework of a collaborative project between Mitsubishi Electric Plant Engineering Corporation (MPEC), Tohoku Small Hydropower Co., Ltd., and Waseda University, the complete machine from spiral casing to draft tube is redesigned and replaced. The project aims at designing a single machine (with a specific speed of Ns = 120 m-kW) capable of operating at both flow rate values to reduce the manufacturing cost on one hand, and to reduce the maintenance and runner exchange costs on the other hand. The geometry of all the components of the new machine, including Spiral Casing (SC), Stay Vanes (SV), Guide Vanes (GV), Runner Vanes (RV) and Draft Tube (DT), has been designed by Waseda University by targeting maximum efficiency and power output at both flow rate values. A 3D view of the final machine is provided in Figure 1, together with its specifications. The design procedure for the complete machine is out of the scope of this paper and will not be discussed in the following.  The efficiency and flow rate as a function of the guide vane opening (GVO) angle obtained by performing steady CFD simulations of the internal flow of the final machine from SC to DT at the prototype scale are given in Figure 2. According to this figure, the machine will operate continuously in off-design conditions, either in part-load conditions (summer flow rate, with GVO ≈ 40%) or full-load conditions (winter flow rate, with GVO ≈ 100%). Therefore, the machine will be subject to the development of flow instabilities in the draft tube, notably at part-load conditions for which a precessing vortex core (PVC) may develop in both cavitation and cavitation-free conditions. In the first step of the project, the machine components including the draft tube are designed by targeting minimum energy losses without considering the development of draft tube flow instabilities and pressure fluctuations, which would require time-consuming unsteady CFD simulations. In a second step, an original approach is proposed to define a draft tube geometry enabling us to: • Minimize the energy losses in the draft tube at both flow rate values; • Mitigate the flow instability and minimize pressure fluctuations at part-load conditions; • Ensure the necessary pressure recovery to achieve the desired net head.
In this second step, only the draft tube geometry is modified and the components from SC to RV designed during the first step are used as a result and input of the 3D model of the machine. The complete design procedure is summarized in Figure 3. This paper focuses only on the second step, i.e., the procedure to design the final draft tube geometry by considering its impact on the pressure fluctuations at part-load conditions and the energy losses at both load conditions. For this purpose, an original approach based on Design of Experiments (DOE) and an Analysis of Means (ANOM) of the results of unsteady and steady CFD simulations is proposed in the following.
It must be mentioned that the phenomenon of cavitation surge potentially occurring at full-load conditions is not considered in the present paper. First, cavitation surge instabilities occur only in certain cases and may be mitigated by air injection on the prototype machine if the present test-case is subject to unstable full-load conditions. In addition, this phenomenon occurs in cavitation conditions at one unstable eigenfrequency of the complete hydraulic circuit, and its prediction by CFD simulation would require twophase flow unsteady simulations coupled with an one-dimensional model of the complete hydraulic circuit [16], which is out of the scope of this paper. CFD aided-design (steady CFD simulations)

Approach:
Design of Experiments (orthogonal array) Figure 3. Diagram of the design procedure for the present test-case, together with the limit of the scope of this paper.

Simplified Parametrization of the Draft Tube Geometry
Elbow draft-tube of Francis turbines can be modelled by 20 planar sections, each section being described by one or several parameters. This may result in a parametrization of the Francis turbine draft tube by up to hundreds of parameters. Usually, the number of parameters is reduced by using a simplified parametrization with 9 to 20 parameters, such as in [23][24][25] among others.
In the studies mentioned previously, the optimization of the draft tube geometry targets the maximization of the performance of the draft tube, which can be achieved by performing steady CFD simulations of the draft tube domain only [25] or the coupled domains runner-draft tube [23]. In the present study, the minimization of the pressure fluctuations in part-load conditions is also targeted, which requires unsteady CFD simulations of the draft tube domain together with the runner domain, increasing, therefore, the computational cost. To decrease the latter, the number of parameters modelling the draft tube geometry is reduced to the minimum, as illustrated in Figure 4. This is made possible by the fact that some parameters are imposed by the implementation of the existing machine, including the difference of elevation between the runner outlet and the center of the draft tube outlet section (∆Z in Figure 4), and the distance between the machine center axis and the draft tube outlet section (∆X in Figure 4). In addition, the size of the considered hydropower unit allows to make additional simplifications. First, only circular planar sections are considered along the draft tube, whereas more complex section shapes, including oval and rectangular ones, are usually considered beyond the draft tube cone in the design of Francis turbines of medium and large size. Finally, a linear area diffusion curve is considered in the elbow, whereas it may feature a peak point and then a trough point in the elbow [18]. Based on these simplifications, the geometry of the draft tube can be modelled and fully defined by only 4 free parameters (listed below and defined in Figure 4):  The parameter 3, i.e., the elbow section ratio, is defined as the difference between the outlet diameter and inlet diameter of the elbow normalized by the length of the elbow centerline L elb. = πR elb. /2. It can be interpreted as the variation of the section diameter per unit of length along the elbow. In addition, the centerline of the elbow is defined as a quarter circle with a curvature radius R elb. , which depends on the length of the cone L cone .

Design of Experiments Applied to the Draft Tube Geometry
In the literature and industrial projects, the optimization of Francis turbine draft tube geometry usually targets maximum efficiency, i.e., minimization of the energy losses and maximization of pressure recovery, the latter being specially important for large hydropower stations with large flow rate values. Optimization algorithms such as stochastic evolutionary algorithms [25], multi-objective genetic algorithms [23], Particle Swarm Optimization algorithms [26], have been widely applied for the design of Francis turbine components, including runner and draft tube. In this case, steady CFD simulations are performed to estimate the losses and efficiency for each design candidate. The application of optimization algorithms may require a total number of design candidates up to several hundreds. If the minimization of pressure fluctuations in off-design conditions is also targeted as in the present test-case, unsteady CFD simulations of the draft tube internal flow are required to capture the development of unsteady vortical structures and corresponding pressure fluctuations. Therefore, the use of optimization algorithms is irrelevant and unfeasible in this case due to the computational cost it would require.
In the following, a method based on a state-of-the-art Design of Experiments (DOE) approach is proposed for the design of the draft tube component. This approach enables a drastic reduction of the number of required simulations compared with optimization algorithms, and will provide additional insights regarding the influence of each geometrical parameter (called factor in the DOE terminology) on the characteristic function, or objective function, to be minimized (the pressure fluctuations and the energy losses in the present test-case). DOE has been widely applied to turbomachines, including Francis turbine runners [27,28]. However, from the authors' knowledge, this approach has never been applied to unsteady CFD simulation results of hydraulic turbines for the evaluation of the influence of the design parameters on draft tube flow instabilities and resulting pressure fluctuations.
The factorial design approach, which considers all possible combinations of levels, is an expensive approach due to the number of required experiments if the number of factors is higher than 2. Therefore, a fractional design approach based on orthogonal arrays is applied in the following. For each factor, three different levels are considered, see Table 1. The method of Taguchi [29] is used to express the orthogonal array, as follows: with N the number of rows for the orthogonal array, n the number of factors, s i the number of levels and k i the number of factors with s i levels. Given the number of factors (n = 4) and the number of levels for each factor (s i = 3), the orthogonal array L 9 (3 4 ) with 9 rows is used.
Each row corresponds to a specific draft tube design (called experiment), characterized by a certain combination of levels defined by the theory of the orthogonal array (not described in the following, the reader can refer to [30]). The orthogonal array is provided in Table 2, and the 9 resulting draft tube geometries are illustrated in Figure 5. Table 1. Geometrical parameters (factors) and corresponding levels used for the DOE analysis. Table 2. Level of the factors for the 9 different draft tube geometries.

Level of the Factors
For each draft tube design, the internal flow is computed by CFD simulations at two different operating points: one at part-load conditions with a guide vane opening angle corresponding to 40% of the maximum opening and one at full-load conditions with a guide vane opening angle corresponding to 100% of the maximum opening. These two operating points correspond to the required flow conditions on the prototype machine, i.e., Q = 0.6 m 3 s −1 for GVO 40% (summer conditions) and Q = 1.5 m 3 s −1 for GVO 100% (winter conditions). It must be mentioned that steady CFD simulations are performed at full-load conditions whereas unsteady CFD simulations are performed at part-load conditions to capture the unsteady behaviour of the draft tube flow and the resulting pressure fluctuations. The numerical set-up is presented in the following section.

Numerical Set-Up
The flow is modeled based on the homogeneous Reynolds-Averaged Navier Stokes (RANS) and Unsteady Reynolds-Averaged Navier Stokes (U-RANS) at GVO 100% (fullload conditions) and GVO 40% (part-load conditions), respectively. The Shear Stress Transport (SST) and SST-Scale-Adaptive Simulation (SST-SAS) turbulence models [31] are adopted at full-load and part-load conditions, respectively. Both models are widely used for the simulation of steady and unsteady flows in hydraulic machines [32][33][34][35]. The flow is simulated by using the commercial CFD code ANSYS CFX 19.0. The computational model goes from the guide vanes inlet to the draft tube outlet, and includes the full runner domain and only one guide vane to reduce the computational cost. Neglecting the spiral casing and stay vanes domain, and considering one single guide vane have already provided results in good agreement with experimental results in the literature [36]. The parameters of the CFD model are summarized in Table 3. The high resolution advection scheme is used for both operating points. For the unsteady simulation at partload, the second order backward Euler scheme is used for the unsteady scheme. No slip condition is applied to the solid walls. Interface GV-RV Stage mixing plane with constant total pressure Interface RV-DT Transient rotor-stator Frozen rotor Total pressure and flow direction are imposed as boundary condition at the inlet section, while the static pressure is imposed at the outlet section. The flow conditions at the guide vane inlet, including the flow direction and the total pressure value, are determined based on the results of a steady simulation of the full domain (including the spiral casing and the stay vanes) and the total available head of the power station. In terms of interface between the domains, the interface between the guide vane and the runner domain is set as stage mixing plane with constant total pressure at both operating points, and the interface between the runner and draft tube domains is set as frozen rotor and transient rotor stator at full-load and part-load, respectively. The convergence criteria are defined as RMS below 10 −5 for both simulation types.
For the unsteady simulation at part-load, the time step is set at 3.6 deg. of the runner rotation, for a total simulation time of 120 runner rotations. The runner frequency n (s −1 ) is set at the actual synchronous frequency on the prototype machine.
Since CFD simulations are performed at the design stage of the machine, no experimental data is available for comparison and validation of the CFD model at both operating points. However, previous numerical studies using a similar CFD modelling for the prediction of hydraulic turbines internal flow have shown a good agreement with experimental results measured during model test, demonstrating the capability of this numerical set-up to reproduce the main features of the draft tube flow behaviour in off-design conditions [37,38].
The flow domains are meshed by hexahedral elements. The structured mesh in the guide vane and runner domains are generated by using the turbine blade meshing commercial software Ansys Turbogrid 19.0 while the structured mesh in the draft tube domain is generated by Ansys ICEM CFD 19.0. The influence of the mesh density on the simulated results is investigated by RANS simulations performed at one operating point close to the design point with a GVO equal to 90% of the maximum value, for which the draft tube flow is stable and does not feature unsteady coherent structures. The draft tube geometry 1 is used for the mesh dependency study.
Five different meshes are generated, with a number of elements ranging from 3,385,000 (mesh 1) to 30,200,000 (mesh 5), as given in Table 4, while the average y + -value is kept constant. The influence of the elements number on the value of the flow rate, runner losses and draft tube losses is given in Figure 6a-c, respectively. Based on these results, the mesh number 3 is used in the following study to reduce the computational cost while minimizing the influence of the mesh density on the results. The characteristics of the final mesh, including number of elements in each domain and average value of y + at both part-load and full-load conditions are given in Table 5. For the draft tube domain, a range of values is given since the number of elements and the average y + slightly vary from one draft tube geometry to another.

Preliminary Computational Results
The operating parameters of the machine obtained at both operating points for each draft tube geometry are presented in Figure 7. It includes the flow rate, the net head of the machine H e (from GV to DT), the hydraulic efficiency η and the output power P. The parameters are defined as follows: where Pt inlet and Pt outlet are the total pressure at the GV inlet and DT outlet, respectively, ω = 2π · n is the angular velocity of the runner, T is the total torque applied to the runner, ρ is the water density and g is the gravity. It must be mentioned that the net head is usually based on the total pressure difference between spiral casing inlet and draft tube outlet. The simulated flow rate is not significantly affected by a modification in the draft tube geometry at both operating points, the maximum difference being 2.4% at GVO 40% between the results in draft tube 1 and draft tube 4, and 4% at GVO 100% between the results in draft tube 1 and draft tube 2. The effective net head only slightly fluctuates for the draft tube geometries 2 to 9, but a deficit up to 5 m is observed in the case of draft tube 1 compared with the values in the other cases. This is explained by the absence of pressure recovery in this case since the inlet and outlet sections of the draft tube 1 have equal diameter. The geometry of the draft tube also affects the efficiency of the machine, with values ranging from 91.9% to 94.2% at GVO 100% and from 90.7% to 92.3% at GVO 40%.  Although the efficiency is maximum in the case of the draft tube 1 at both operating points, this geometry does not provide the maximum output power at both conditions (1.49 MW for the draft tube 1 against 1.63 MW for the draft tube 2 at full-load, and 0.66 MW for the draft tube 1 against 0.69 MW for the draft tube 4 at part-load), which is explained by the deficit of effective head due to the absence of pressure recovery.
In the following section, the influence of the draft tube geometry on the behaviour of the precessing vortex core (PVC) and corresponding pressure fluctuations at part-load conditions will be first illustrated in two cases. Then, an ANOM will be applied to the numerical results in terms of pressure fluctuations and energy losses in the draft tube. For all draft tube geometries, the pressure fluctuations at the wall are monitored in two different sections of the draft tube cone during unsteady CFD simulations at GVO 40%, see Figure 8. In each section, the pressure is monitored at four different locations regularly spaced by an angle of 90 • . Section 1 is located 0.18 × D 1 downstream of the runner outlet in all draft tube geometries, whereas the z-coordinate of Section 2 changes from one draft tube to another due to the variation of the draft tube cone length. Section 2 is, however, placed 0.03 × D 1 upstream of the elbow inlet in all draft tube geometries.

Results and Analysis
A cross-spectral analysis between the pressure signals simulated at the monitoring points p 1 and p 3 is performed in the case of DT3 and DT4. The results for DT3 and DT4 are given in Figure 9a   In both cases, a clear frequency peak is observed at f = 0.46 × n and f = 0.38 × n for DT3 and DT4, respectively. This frequency corresponds to the vortex precession frequency f PVC that locally induces pressure fluctuations on the draft tube wall. The amplitude of these pressure fluctuations in the case of DT3 is slightly higher than that in DT4. The phase shift at f PVC between p 1 and p 3 (bottom figures) reveals the nature of the pressure fluctuations in the corresponding cross-section: in the case of DT3, the phase shift is almost equal to 0, which corresponds to synchronous (or in-phase) pressure pulsations, whereas it is equal to ≈π/2 in the case of DT4. The precession of a vortex core in the draft tube of Francis turbine is known to induce two different types of pressure fluctuations, both occurring at the precession frequency f PVC . The convective pressure component, or asynchronous, is locally induced by the rotation of the inhomogeneous pressure field with the vortex core. If the pressure fluctuations measured by two sensors located in the same cross-section and spaced by an angle θ are purely convective, the phase shift at f PVC between both pressure signals is equal to θ. In the case of an elbowed draft tube, the interaction between the PVC and the elbow dissymmetric geometry induces an additional type of pressure fluctuations at f PVC , the so-called synchronous component [39]. It has equal phase and amplitude in a given cross-section, and propagates through the complete hydraulic circuit, including the pipe upstream of the machine. Therefore, the pressure fluctuations measured in a Francis turbine elbowed draft tube correspond to the superimposition of both components. The phase shift between two pressure sensors in the same cross-section is not equal to neither 0 (pure synchronous nature) nor the angular difference between both sensors (pure convective nature) but to a value depending on the amplitude and phase of each component. For the present case, it can be assumed that the synchronous component in the case of DT3 is dominant in the corresponding section of interest, leading to a phase shift close to 0.
Several mathematical tools have been developed to decompose the pressure fluctuations induced by a PVC in Francis turbine draft tube into synchronous and convective components [4,40,41]. If N sensors are equally distributed along the wall in one crosssection, the pressure components can be decomposed by using the following equations [42]: where p syn corresponds to the synchronous pressure component common to all pressure sensors, p conv,i is the convective pressure component relative to the pressure sensor i and p i is the individual raw pressure signal i. In the present case, four monitoring pressure points located in Section 1 (N = 4) are used to calculate the pressure components. A Fast Fourier Transform (FFT) is then applied to each pressure component. The FFT obtained in Section 1 are given in Figure 10a,b for the cases DT3 and DT4, respectively. As concluded by analysing the phase of the cross-spectral density function, the synchronous component is dominant in Section 1 of DT3, whereas the convective component amplitude is equal to twice the amplitude of the synchronous component in the DT4. However, the amplitude of the convective components in Section 1 is similar in both draft tubes while the synchronous component in DT3 is almost five time higher than the one in DT4. The results obtained in Section 2 are given in Figure 11. While the amplitude of the synchronous pressure components is close to the one observed in Section 1 in both draft tubes, the amplitude of the convective components in DT3 remarkably increases compared with the one in Section 1 (see Figure 10).
The amplitude of the convective component in a given section is directly related to the PVC trajectory radius in this section, as already demonstrated in simplified test cases generating swirling flows [43]. In addition, the synchronous component results from the interaction between the precession motion of the vortex core and the dissymmetric geometry of the elbow, and it has been assumed that a wider PVC trajectory in the elbow may induce higher synchronous pressure pulsations [3].
To understand the differences observed between both draft tubes in terms of pressure components amplitude, the instantaneous position of the vortex centre is estimated in 4 cross-sections along the draft tube, including Sections 1 and 2, and 2 additional sections in the elbow. A total of 30 time steps are considered, corresponding to approximately three complete PVC rotations. To determine the PVC center position at each time step, the vortex centre identification algorithm proposed by Graftieaux et al. [44] is applied. It relies on the computation of the criterion Γ 1 , which is only based on the flow topology and does not involve any velocity gradients. It has been successfully used for the identification of PVC in hydraulic turbines [37] and simplified swirling flow applications [45], and tip vortices in hydrofoil applications [46].
The instantaneous positions of the PVC centre are given at 30 time steps in Figures 12 and 13 for DT3 and DT4, respectively. The black dots indicate the position of the consecutive instantaneous vortex centres while the red line corresponds to a mean PVC trajectory based on the 30 instantaneous vortex centres by considering a pure circular trajectory.
In DT3, the vortex center trajectory is very close to the section center in Section 1 and then strongly widens along the draft tube cone, with a vortex trajectory radius r v in Section 2 equal to six times the radius in Section 1. It should be noted that the dispersion of the instantaneous vortex centre around the mean trajectory is remarkably low, highlighting a high level of periodicity in the precession motion. These features are also observed in the elbow, even if the centre of the vortex trajectory is slightly shifted from the geometrical centre of the section. This strong periodical precession of the vortex in the elbow may be responsible for the high-amplitude synchronous pressure pulsations observed in the draft tube cone. Therefore, the results in Figure 10 can be explained by a weak vortex precession in Section 1 (low convective pressure fluctuations) combined with a strong, periodical vortex precession in the elbow producing high synchronous pressure fluctuations that propagate upstream.  In DT4, the evolution of the vortex precession trajectory is remarkably different from the one in DT3. First, the instantaneous vortex centres feature a higher dispersion around the mean trajectory compared with the previous case. In addition, the radius r v of the vortex trajectory remains similar in the draft tube cone and then slightly widens along the elbow, but with a radius much lower than the one in DT3. These features, combined with the weak periodicity of the vortex precession, explain the notable difference observed between the amplitude of the pressure components in DT3 and DT4. In DT4, the precession of the vortex in the elbow only produces low-amplitude synchronous fluctuations, making the convective pressure fluctuations at the runner outlet dominant. This preliminary analysis, based on the results obtained in only two draft tubes featuring radically different geometries, clearly highlights the influence of the draft tube geometry on the structure of the PVC and the amplitude of the associated pressure fluctuations. Moreover, it also reveals the complexity of the phenomenon: in a given geometry, the structure of the PVC evolves along the draft tube, producing convective pressure fluctuations whose amplitude depends on the local vortex trajectory and changes along the draft tube. The synchronous pressure component is directly related to the precession of the vortex in the elbow, and the study of the impact of the draft tube geometry on the flow stability should ideally consider the vortex structure along the complete draft tube.
However, as this is not feasible from a practical point of view in the case of an analysis by DOE and ANOM, the following will only focus on the pressure fluctuations (including both convective and synchronous components) at the runner outlet in the monitoring Section 1 of the draft tube cone. This simplification is supported by the two following arguments. First, the runner will be impacted by the local convective pressure fluctuations, i.e., the precession of the vortex just at the runner outlet. In addition to that, the synchronous pressure pulsations measured at the runner outlet also propagate upstream, through the runner blade channels, and affect directly the runner dynamic stresses and fatigue. Therefore, minimizing the pressure fluctuations at the runner outlet will result in a reduction of the dynamic stresses on the runner blades and an improvement of the runner lifetime. One should, however, keep in mind that the development of convective pressure fluctuations farther in the draft tube may induce additional vibrations of the stationary components and additional fatigue. This aspect is, however, out of the scope of the present paper.

Analysis of Means of Pressure Fluctuations in Part-Load Conditions
In this section, the influence of each geometrical parameter on the amplitude of the pressure fluctuations simulated in the monitoring Section 1 is assessed. The pressure fluctuations are first made dimensionless by defining the pressure coefficient C p as follows: where p is the pressure,p is the time-averaged pressure and H is the available head of the hydropower unit. Pressure fluctuations in Francis turbines are commonly normalized by the net head H e of the machine. However, in the present study, the net head varies from one case to another (see Figure 7), making the normalization by the net head irrelevant for further comparisons between the nine different cases. The Root-Mean-Square (RMS) value of the pressure coefficient C p, i , where the index i = 1 . . . 9 is the number of the experiment, is the characteristic function of the experiments. The mean value m of this characteristic function is calculated by considering its value over the nine experiments as follows: The pressure coefficient RMS value C p, i is given for each DT geometry in Figure 14, together with the mean value m. The DT geometry has a great influence on the amplitude of the pressure coefficient C p , its value ranging from 0.05% (DT1) to 0.57% (DT8) while the mean value m is equal to 0.36%. The effect of each level of the geometrical parameters on the pressure coefficient amplitude is assessed by applying an ANOM approach [30] to the results given in Figure 14. Considering a level j (j ranging from 1 to 3) of a geometrical parameter a k (k ranging from 1 to 4), the effect of the geometrical parameter a k on the characteristic function when equal to its level j can be estimated by the difference between the average value of the characteristic function over the experiments featuring the same level j of the parameter a k (noted m a k =j ) and the mean value m. For example, the first level of the factor a 1 corresponding to the DT cone angle, i.e., a 1 = α cone = 0 deg. (see Table 1), is used in three different DT geometries among the nine DT geometries (experiments DT1, DT2 and DT3). Therefore, the effect of the level 1 of the factor a 1 on the pressure coefficient amplitude can be estimated as follows [30]: For each geometrical parameter, the effect of the three different levels on the pressure coefficient amplitude is estimated by using Equation (9). The results are reported in Figure 15. In this figure, the mean characteristic function m is not subtracted to the mean value of each level m a k =j to avoid negative value, but the mean value is indicated by the horizontal dashed line. Therefore, values below/above this line correspond to a mitigating (green zone in Figure 15)/enhancing (red zone in Figure 15) effect with respect to the mean value m, respectively. Based on the ANOM of the pressure coefficient amplitude presented in Figure 15, the following conclusions can be drawn. First, the parameters with the most influence on the pressure coefficient amplitude are the DT cone angle and the elbow section ratio. For both parameters, as their value is increased, the amplitude of the pressure coefficient is amplified. This is particularly significant for the elbow section ratio: when this parameter is modified from its lowest value (SR elb. = 0) to its maximum value (SR elb. = 0.2), its effect on the pressure coefficient amplitude is multiplied by a factor of about 2.5. At its lowest level, this parameter decreases the amplitude of the pressure coefficient with respect to the mean value m (mitigating effect of −0.16%), while the amplitude is increased with respect to the mean value m when this parameter is set at its highest level (enhancing effect of +0.12%). In the case of the DT cone angle α cone , its effect on the pressure coefficient amplitude varies from −0.07% (mitigating effect) with α cone = 0 deg to +0.07% (enhancing effect) with α cone = 4 deg.
The influence of these two parameters can be explained by the fact that increasing the diffuser effect in the cone results in a higher adverse pressure gradient, increasing the potential flow recirculation and the PVC trajectory radius. This, in turn, induces local convective pressure fluctuations with a higher amplitude. Concerning the influence of the elbow section ratio, it was shown in Section 4.1 that an increase in the section area in the elbow may induce a wide PVC trajectory in the cone on one hand, and high synchronous pressure pulsations on another hand. Reducing the diffuser effect in both DT cone and elbow seems therefore crucial to reduce the amplitude of the pressure fluctuations at the runner outlet. The DT cone length has also a slight influence on the amplitude of the pressure pulsations, while the influence of the diffuser angle does not appear clearly: the minimum value of the pressure coefficient amplitude is observed at α di f f . = 0 deg, then its value increases and decreases as the value of α di f f . is increased.
Based on the ANOM, it can be concluded that the combination of levels that produce the minimum amplitude of pressure fluctuations in the monitoring Section 1 corresponds to the draft tube geometry 1, i.e., the value of each geometrical parameter is set at its first level. This is in agreement with the comparison between the results obtained with the different draft tube geometries, see Figure 14.
However, using this geometry as a final design is not a reliable solution, since the outlet and inlet sections of the draft tube have equal diameter: there is no pressure recovery in the draft tube, which leads to a drop of the net head H e of the turbine and of the output power P compared with the other geometries, as shown in Section 3.4. This effect is particularly important at full-load conditions (GVO 100%) with a higher flow rate and remaining kinetic energy in the draft tube. In the following, the energy losses and pressure recovery in the DT are considered.

Analysis of Means of Draft Tube Energy Losses
The ANOM is applied to the energy losses in the draft tube at both part-load (GVO 40%) and full-load (GVO 100%) conditions. The losses H r,DT in the draft tube are estimated based on steady-state and unsteady-state CFD simulations at GVO 40% and GVO 100%, respectively, by considering the difference of total pressure between the inlet and outlet of the draft tube: The results of the ANOM of the draft tube losses are given in Figure 16. The elbow section ratio SR elb. and the diffuser angle α di f f . have a clear influence on the energy losses in the draft tube at both operating points: the effect of both parameters continuously increase as the value of their level is increased. This is particularly significant for the elbow section ratio at full-load conditions: it has a mitigating effect of −0.64% with respect to the mean value m as the parameter is set at its first level, whereas it has an enhancing effect of +0.87% with respect to the mean value m as the parameter is set at its third level.
The influence of the DT cone angle and length on the losses is more ambiguous. At both operating points, the losses decrease and then increase as the value of the cone angle is increased. However, this effect is weak and the minimum losses value is only 9% lower than the maximum value. Concerning the DT cone length, the losses continuously increase as the value of this parameter is increased at GVO 40%, whereas they decrease and then increase in the case of GVO 100%. The combination of geometrical parameters inducing the minimum energy losses in the draft tube is given in Table 6. Concerning the level of the DT cone length, the one minimizing the losses at GVO 100% is selected, since its effect at GVO 40% is negligible. The combination minimizing the pressure fluctuations amplitude (previous section) is also included in Table 6.
The minimization of either the energy losses or the pressure fluctuations amplitude does not lead to the same set of geometrical parameters and are therefore conflictive. In addition, the pressure recovery must be considered to ensure a sufficient effective head for the turbine. In the following section, the geometry featuring the best compromise between sufficient pressure recovery, minimization of both draft tube losses and pressure fluctuations is discussed based on additional analysis. The ANOM approach is finally validated by comparing the results provided by CFD simulation of the final design with the values estimated by ANOM.

Determination of the Final Draft Tube Geometry
The value of a characteristic function (i.e., the draft tube losses or the pressure coefficient amplitude in the present case) for any combination of level a k = j can be estimated based on the results provided by the ANOM as follows [30]: For instance, for the combination of parameters minimizing the draft tube energy losses (combination 1 in Table 6), the value of the draft tube losses H r,DT at GVO 40% can be estimated as follows: Equation (11) is applied to estimate the value of the draft tube energy losses at both GVO 40% and GVO 100% and the pressure coefficient amplitude at GVO 40% for both combinations of geometrical parameters given in Table 6. The results are given in Table 7. As shown in Table 7, the energy losses in the draft tube and the amplitude of the pressure coefficient cannot be minimized at the same time. The amplitude of the pressure coefficient is for instance multiplied by 3.5 if the combination 1 minimizing the losses is used instead of the combination 2. A compromise between mitigation of energy losses and pressure fluctuations amplitude may be reached by setting the geometrical parameters 1 to 3 to their first level, while the diffuser angle (parameter 4) will be set at a value enabling a sufficient pressure recovery, as described in the following.
The pressure recovery in a Francis turbine draft tube corresponds to the increase in the static pressure due to the deceleration of the flow along the draft tube induced by the increase in the cross-section area. This enables the reduction of the losses at the draft tube outlet by minimizing the remaining kinetic energy, thereby increasing the net head of the turbine.
The ideal pressure recovery H recov in the draft tube is defined as the difference of static pressure between the inlet and outlet of the draft tube without considering the energy losses, i.e.,: where R 1 and R 4 correspond to the radius of the draft tube section at the inlet and outlet, respectively, as defined in Figure 4. The influence of the ratio R 4 /R 1 on the ideal pressure recovery at both operating conditions is given in Figure 17. A maximum of 6% and 1% of the total available head H can be recovered and transformed into pressure by the draft tube at GVO 100% and 40%, respectively. According to Figure 17, the head recovered by the draft tube does not significantly further increase beyond a ratio R 4 /R 1 higher than 2 at both operating conditions, i.e., for a diffuser angle α di f f . > 4 deg. The difference between the recovered head with α di f f . = 4 deg. and the maximum value of recovered head does not exceed 0.38% and 0.06% at full-load and part-load conditions, respectively. However, this analysis does not consider the energy losses in the draft tube: the latter induce a decrease in the recovered head as they increase when the diffuser angle is increased (see results of the ANOM in Figure 16).
The balance between the ideal recovered head and the head losses in the draft tube as a function of the diffuser angle is, therefore, considered in the following to determine a suitable value of the diffuser angle. To do so, three additional draft tube geometries are investigated (experiments 10 to 12). They feature different values of diffuser angle, ranging from α di f f . = 1.5 deg. to α di f f . = 4 deg., while the other geometrical parameters are fixed and set at their first level as defined previously. The parameters for each draft tube geometry are summarized in Table 8. For each geometry, steady and unsteady CFD simulations are performed at full-load (GVO 100%) and part-load (GVO 40%) conditions, respectively. The results of the simulation with α di f f . = 3 deg. (experiments 11) are also used to validate the methodology and the consistency of the estimation provided by the ANOM in terms of losses and pressure fluctuations amplitude, since the level α di f f . = 3 deg. is included in the experiments 1 to 9.  The influence of the diffuser angle on the draft tube losses H r,DT , the ideal recovered head H recov. and the actual recovered head ∆H recov. is presented in Figure 18. The actual recovered head corresponds to the difference between the ideal recovered head and the draft tube losses. First, it can be noticed that the pressure recovery is not sufficient at partload conditions and does not compensate for the draft tube losses, leading to a negative value of actual recovered head. However, the pressure recovery plays a crucial role in full-load conditions where maximum 4% of the head H can be recovered.
At both part-load and full-load conditions, the draft tube losses increase as the diffuser angle is increased, which is in agreement with the results obtained by the ANOM in Section 4.3. In addition, the ideal pressure recovery does not significantly increase beyond α di f f . = 4 deg., as highlighted previously. Since the losses increase in the meantime as α di f f . is further increased beyond 4 deg., the actual recovered head ∆H recov remains almost constant from α di f f . = 3 deg. to α di f f . = 4 deg. at GVO 100% and even slightly decreases in the case of GVO 40%. It can be expected than increasing further the diffuser angle may lead to a decrease in both the actual recovered head and the net head of the turbine.
The influence of the diffuser angle on the pressure fluctuations amplitude in Section 1 at part-load conditions is provided in Figure 19. The tendency is similar to the one obtained by the ANOM in Section 4.3. The amplitude first increases as the diffuser angle is increased and it then decreases beyond α di f f . = 3 deg.

Validation and Discussion
To validate the methodology of DOE-ANOM, the results obtained by CFD simulations of the draft tube geometry with α di f f . = 3 deg. (experiment 11 in Table 8) are compared with the results estimated by using the ANOM. The value of the draft tube losses and pressure fluctuations are estimated by applying Equation (11) to the results of the ANOM given in Sections 4.2 and 4.3. The geometrical parameters are fixed at their first level, except the diffuser angle (level 2, α di f f . = 3 deg.). The comparison between simulated and estimated results is given in Table 9. The simulated and estimated results are in good agreement at part-load conditions, with an overestimation by the ANOM approach of +12% for the draft tube losses and +9% for the pressure fluctuations amplitude. The deviation is, however, much higher for the draft tube losses at full-load conditions. The difference between both operating conditions may be explained by the fact that only steady CFD simulations have been performed at full-load conditions. In addition, the interactions between the different parameters have not been considered and may explain this deviation between estimated and simulated results at GVO 100% [30]. Further investigations may focus on the influence of such interactions between the parameters on the results in full-load conditions. However, the approach proposed in this paper provides a reliable estimation of the influence of each draft tube geometrical parameter on the vortex-induced pressure fluctuation in part-load conditions and the energy losses in both full-load and part-load conditions. This makes possible the determination of a suitable combination of geometrical parameters minimizing both pressure fluctuations amplitude and energy losses, while keeping the required pressure recovery in the draft tube. This methodology based on Design of Experiments and Analysis of Means of CFD simulations results in the draft tube flow represents an original alternative to state-of-the-art design approaches based on optimization algorithms. The latter remain unfeasible from a practical point of view if the minimisation of pressure fluctuations in off-design conditions is targeted, since this would require to perform unsteady CFD simulations for several hundred of design candidates. By applying the proposed approach, the computational time and cost are drastically reduced in the case of simple draft tube geometries. In addition, it also enables a rough preliminary estimation of the characteristic function of interest (the pressure fluctuations amplitude and energy losses in the present case) for several combination of geometrical parameters by using the results of the Analysis of Means.

Conclusions
An original method is proposed in this paper to characterize the influence of the draft tube geometry on the pressure fluctuations amplitude and energy losses during off-design conditions. The method relies on the application of state-of-the-art Design of Experiments (DOE) and Analysis of Means (ANOM) approach to the results of unsteady and steady CFD simulations of the draft tube internal flow. The proposed methodology is applied to a Francis turbine unit that is required to operate continuously in off-design conditions, either with 45% (part-load) or 110% (full-load) of the design flow rate.
A set of nine different draft tube geometries featuring different sets of geometrical parameters is first defined by applying an orthogonal array-based DOE approach. For each geometry, unsteady and steady CFD simulations are performed at part-load conditions (with a guide vane opening angle equal to 40% of the maximum value) and full load conditions (maximum guide vane opening), respectively. By applying an ANOM of the results obtained in the nine draft tube geometries, the complexity of the problem is reduced and the influence of each geometrical parameter on the pressure fluctuations amplitude and energy losses is evaluated. Among the different findings, it is confirmed that the diffuser effect along the draft tube has a great impact on the development of the PVC and resulting pressure fluctuations. By considering the ANOM results and the required pressure recovery in the draft tube, an ideal set of geometrical parameters can be finally deduced. The CFD simulation results of a final draft tube geometry are finally compared with the estimation provided by the ANOM. The good agreement, especially in part-load conditions, demonstrates that the method also enables a rough preliminary estimation of the draft tube losses and pressure fluctuations amplitude.
The methodology proposed in this paper represents an original approach to estimate the influence of the draft tube design on vortex-induced pressure fluctuations and energy losses during off-design conditions with reasonable computational costs and time. In addition, it is an alternative to state-of-the-art design approaches based on optimization algorithms for the design of draft tube geometry considering flow instabilities in off-design conditions. Such approaches remain unfeasible from a practical point of view if the minimization of pressure fluctuations is targeted, since this would require us to perform unsteady CFD simulations for several hundred of design candidates. By applying the proposed approach, the computational time and cost are drastically reduced in the case of simple draft tube geometries. Therefore, this may pave the way to new design procedures for the development of hydraulic machines draft tube capable of smooth operations in off-design conditions to fulfill with the new requirements of the electrical power system in terms of flexibility.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to confidentiality agreement with industrial partners.