RANS computations for identification of 1-D cavitation model parameters: application to full load cavitation vortex rope

Due to the massive penetration of alternative renewable energies, hydropower is a key energy conversion technology for stabilizing the electrical power network by using hydraulic machines at off design operating conditions. At full load, the axisymmetric cavitation vortex rope developing in Francis turbines acts as an internal source of energy, leading to an instability commonly referred to as selfexcited surge. 1-D models are developed to predict this phenomenon and to define the range of safe operating points for a hydropower plant. These models involve several parameters that have to be calibrated using experimental and numerical data. The present work aims to identify these parameters with URANS computations with a particular focus on the fluid damping rising when the cavitation volume oscillates. Two test cases have been investigated: a cavitation flow in a Venturi geometry without inlet swirl and a reduced scale model of a Francis turbine operating at full load conditions. The cavitation volume oscillation is forced by imposing an unsteady outlet pressure conditions. By varying the frequency of the outlet pressure, the resonance frequency is determined. Then, the pressure amplitude and the resonance frequency are used as two objectives functions for the optimization process aiming to derive the 1-D model parameters.


Introduction
Due to the extensive development of new renewable energy sources during the last decade, the electrical transmission system undergoes large power flow fluctuations, since these energy sources are strongly dependent on the meteorological conditions. In order to manage these fluctuations, other sources of energy are used to help stabilizing the electrical grid. Among them, hydraulic power plants, taming also green energy, are used effectively due to their ability to respond quickly to a variation of the load. However, in order to inject the suitable amount of power in the grid, the hydraulic turbines have to operate far from their best efficient point. In the case of Francis turbines, the turbine off design operating conditions yield to the development of flow instabilities such as vortex ropes in the turbine draft tube. If the pressure inside the vortex rope reaches the vapour pressure, then cavitation appears which may induce pressure fluctuations that propagate in the whole hydraulic system. The interaction between this excitation source and the system may lead to resonance phenomena at part load [1]. At full load, self-excited pressure and discharge oscillations at the system's eigenfrequency are observed [2]. To investigate these pressure oscillations, a one-dimensional (1-D) unsteady mathematical model of the system is set up, including a model of the cavitation dynamics [3][4][5][6][7][8][9][10][11][12][13]. However, the aforementioned 1-D cavitation dynamics model is based on physical parameters which need to be known trough either experimental or numerical simulation approaches. In the framework of the European FP7 research project Hyperbole, a thorough set of experimental investigations are performed with the reduced scale physical model of a Francis turbine, including the identification of the wave speed and the second viscosity parameters of the cavitation dynamics model [14]. These parameters can also be derived from RANS numerical simulations of cavitation flows. Indeed, RANS models coupled with a cavitation model enable to investigate the dynamics of the cavitation vortex rope. Previous investigations with CFD simulations already identified the cavitation compliance and the mass flow gain factor parameters [15][16][17][18][19]. However, extracting from CFD simulations, the second viscosity parameter modelling the dissipation induced by phase change during cavitation volume fluctuations is less common. Analytical formula of the second viscosity combined with CFD simulation results allowed to derive this parameter [20]. However, the used formula was initially developed for uniform bubbly flow in pipes which flow configuration is far from the cavitation vortex rope. In the case of a 2-D Venturi geometry, an approach based on flow excitation through the outlet boundary condition of the CFD domain has been performed, allowing to compute the cavitation compliance and the second viscosity of this particular case study [21][22].
In the present paper, this approach is used to compute the aforementioned parameters by comparing the system dynamics response of both the CFD and 1-D unsteady models. The comparisons are achieved by using the CFD results as objective functions for the 1-D model. The main results of this approach applied to the Venturi geometry are summarized in this paper and the methodology has been applied to the full load cavitation vortex rope of a Francis turbine. The set of parameters obtained from the identification process are compared between the two test cases.

1-D cavitation model
The modelling of the draft tube cavitation flow is based on both continuity and momentum equations (1) and (2) including the convective terms and the divergent geometry [13].
This set of equations involve four cavitation vortex rope parameters to be found: -the local wave speed a in a control volume of length dx which yields the value of the cavitation compliance c C according to equation (3): -the mass flow gain factor  corresponding to a variation of the cavitation volume as function of the inlet discharge by keeping constant pressure, see equation (4) . In the case of the Francis turbine, this parameter can be interpreted as a variation of the cavitation volume as function of the inlet swirl  [13]. Indeed by varying the inlet discharge with a constant runner rotational speed, the swirl is modified. The cavitation volume as a function of the upstream (inlet) and downstream discharges was recently measured on a simplified test configuration with a microturbine [23]. These local continuity and momentum equations (1) and (2) are respectively integrated over continuity and momentum control volumes which are overlapped in space leading to a staggered grid. Depending on the cavitation length with respect to the investigated system length a lumped or a distributed model has to be selected. The distributed model is characterized by several continuity control volumes along the draft tube length where equation (1) is applied, whereas the lumped model features one continuity control volume. In the case of the distributed model, the mass flow gain factor of each continuity equation is applied to the inlet discharge state variable and is assumed to be constant along the draft tube, see equation (5). This assumption means that the cavitation volume in each control volume varies at the same rate as function of a variation of the inlet discharge.

Methodology
The proposed methodology to find the 1-D cavitation model parameters is sketched in Figure 1. The momentum excitation source, modelling the pressure source fluctuations induced by the cavitation volume dynamics, is not considered in this study. In the two investigated case studies, the cavitation volumes feature a steady behaviour. Hence, three parameters need to be found: the mass flow gain factor 1  , the wave speed a and the second viscosity ''  . First of all, the mass flow gain factor is determined by following a quasi-static approach: several CFD simulations are performed by keeping constant the pressure just downstream the cavitation volume for different inlet discharges. The variation of the cavitation volume as function of the inlet discharge yields the mass flow gain factor. Then, a global optimization process is performed to derive both the wave speed and the second viscosity parameters. In this optimization process, both CFD and 1-D system responses to a fluctuating outlet static pressure are compared. The set of parameters of the 1-D model is optimized to get the best fitting between the two system responses. By imposing outlet pressure fluctuations in a given range of frequency values to the CFD model, the resonance frequency p f of the system is identified when maximum of cavitation volume and pressure fluctuations are experienced. Two quantities measured in the CFD domain are considered as the objectives for the optimization process: the resonance frequency p f of the system and the pressure fluctuation amplitudes at a given location x h . The optimization process searching the set of 1-D cavitation parameters computes both the eigenfrequency and the forced response in the frequency domain of the 1-D model. By minimizing the errors for the two objectives between the 1-D and the CFD models, the set of 1-D cavitation parameters is found.

Cases studies
Two test cases are considered: a 2-D Venturi geometry and a simplified 3-D Francis turbine, see Figure 2. The flow in a Venturi geometry is featuring the development of a cavitation sheet in the diffuser of the Venturi downstream the throat. This Venturi test case was initially investigated in order to assess if 1-D cavitation model parameters, especially the second viscosity, can be derived from CFD computations.
The Francis turbine geometry is simplified by taking into consideration only the runner and the cone of the draft tube. The turbine is investigated at full load operating condition featuring the development of an axial cavitation vortex rope in the draft tube cone. To mitigate interactions between the cavitation vortex rope and the outlet boundary condition, an extension of the draft tube is added.
The advantage of the Venturi test case is the low computation effort required compared to the 3D Francis turbine and the fact that the flow does not feature a swirling component leading to numerically sensitive computations.

3D CFD Models
The CFD numerical simulations are performed by using Ansys CFX 15.0. The two-phase flow is modelled based on the homogeneous mixture assumption assuming that the two phases shared the same pressure and velocity fields.. For the 2-D Venturi, a structured mesh is set up with 15'000 nodes. For the 3D Francis turbine, a structured mesh is set up in the runner and the cone domains, whereas an unstructured mesh is used for the draft tube extension in order to enhance the dissipation of the vortex rope preventing numerical instabilities if the cavitation volume interacts with the boundary conditions. For the both cases, the mesh resolution has been assessed for a constant outlet pressure by comparisons with experimental measurement [21,24]. The total number of nodes is closed to 4.5 million with 2.5 million of nodes in the runner domain. In both cases, the velocity field is imposed at the inlet and a sinusoidal static pressure is specified at the outlet to excite the system. The inlet velocity field imposed at the runner is extracted from a steady cavitation computation on a full domain [24].

1-D SIMSEN Models
The Venturi and Francis turbine test cases 1-D are modelled with the SIMSEN software which includes the cavitation model described previously, see figure 4. The cavitation zone, where cavitation parameters must be identified, is spatially discretized with a lumped model in the case of the Venturi geometry, whereas a distributed model is used for the Francis turbine draft tube. In the cavitation free region, the wave speed is set to 1'000 m/s. Finally, the boundary conditions of the 3-D model are reproduced in the 1-D model: a constant discharge source at the inlet and a sinusoidal static pressure at the outlet. To derive the mass flow gain factor, several computations are performed by keeping the pressure, the closest to the cavitation volume, constant for different inlet discharges. Hence, the pressure is kept constant by adjusting the outlet pressure. For the Venturi test case, the pressure is kept constant at the plane 3 and for the Francis turbine, the pressure is kept constant at the interface between the cone and the extension of the draft tube, see Figure 2. The evolution of the dimensionless cavitation volume as function of the dimensionless inlet discharge is shown in Figure 5. The variation of the volume is lower in the Venturi test case compared to the Francis turbine. The increase of the inlet discharge at a constant pressure level near the cavitation volume induces changes of the flow field allowing a cavitation volume development. In the case of the Francis turbine, the increase of the inlet discharge changes the inlet swirling flow which increases the vortex core region where cavitation is developed. However, although the flow does not feature swirling flow in the case of the Venturi geometry, the increase of the inlet discharge induces a higher flow separation at the throat where cavitation is developed. By comparing the two curves' slopes it is shown that the cavitation volume increase due to the swirl at the runner outlet is higher than the flow separation at the throat. A linear regression is used to determine the value of the mass flow gain factor for the two test cases. The values are summarized in Table 1. It is shown that the mass flow gain factor is about thirty times lower for the Venturi test case than the Francis turbine test case.   ).
For out of resonance conditions, the amplitude of pressure fluctuations decreases from the outlet to the cavitation sheet, which behaves like a pressure node for the system. However, for resonance conditions the amplitude near the cavitation sheet is higher than the outlet excitation source. Moreover, pressure fluctuations do not feature sinusoidal shape anymore. These typical shapes can be observed in Francis turbine draft tubes at resonance conditions due to system non-linearities. This different behaviour between resonance and non-resonance conditions are also observed for the Francis turbine test case. In Figure 8 and Figure    For the Venturi, the pressure upstream the cavitation volume at Plane 2 is considered whereas for the Francis turbine, it is the pressure at the interface. The goals for the two objective functions are given in Table 2. Finally, by performing the optimization process, the set of 1-D cavitation model parameters leading to a matching of the eigenfrequency and the amplitude of pressure fluctuations are given in Table 3 for the two test cases. For each test case, two identifications have been performed with and without mass flow gain factor to assess its influence on the wave speed and the second viscosity values. In the case of the Venturi, the low value of the mass flow gain factor barely influences the parameters. The effect of the flow separation on the cavitation volume as function of the inlet discharge is negligible for the modelling of the cavitation dynamics. In the case of the Francis turbine, the parameters are slightly changed. Indeed, the mass flow gain factor modelling the effect of the inlet swirl has a stabilizing effect at full load operating point [13], [17]. Hence, damping is added to the system which induces a decrease of the eigenfrequency. To compensate this stabilizing effect and to keep the same eigenfrequency, the wave speed parameter is increased with the second viscosity.  In Figure 10 and in Figure 11, the forced response of the 1-D model is plotted along the system abscissa. The used cavitation parameters consider the mass flow gain factor parameter, see Table 3. The dashed grey line corresponds to the outlet excitation frequency set at the first eigenfrequency of the 1-D system. The two remaining series represent an excitation frequency out of resonance conditions. The CFD pressure fluctuation amplitudes at planes P2 and P4 for the Venturi and at the interface for the Francis turbine are compared for the two investigated frequencies of each case study. The CFD objectives for the identification of the cavitation parameters are surrounded. For the excitation set to the first eigenfrequency the pressure fluctuations in the 1-D model are in good agreement with the amplitudes simulated with the CFD model. Indeed, the cavitation parameters have been identified for this frequency. However, by applying the same set of parameters for a higher excitation frequency, the difference with the CFD model is increased. Although the cavitation volume behaves like a pressure node for the system like in the CFD model, the predicted amplitudes with the 1-D model are lower than the CFD model for both case studies. This difference observed for out of resonance condition could be explained by either frequency dependent cavitation parameters or nonlinearities of the system at resonance conditions.

Conclusion
One dimensional cavitation vortex rope models used for stability analysis of power plants involve parameters to be identified either experimentally or numerically. The second viscosity parameter modelling the dissipation induced by the phase change during cavitation volume fluctuations is decisive to predict the stability limit and the system response to the vortex rope excitation. This study is aiming to show that two phase flow RANS computations are able to derive this parameter. Hence, a methodology has been setup and applied to two test cases. First of all, the methodology has been applied to a 2-D Venturi geometry. With this simplified test case it has been proven that the second viscosity can be identified from RANS computations. However, the values obtained for the wave speed and the second viscosity are very low compared to experimentally obtained parameters for Francis turbines found in the relevant literature. The methodology is then applied to a Francis turbine.