Transient two-phase CFD simulation of overload pressure pulsation in a prototype sized Francis turbine considering the waterway dynamics

At high load operation points, Francis turbines generally produce large cavitation volumes of central vortex character in the draft tube. In order to gain a deeper understanding of the flow behaviour at high load conditions a combined 1D-3D transient two-phase numerical investigation at prototype size was carried out and these results were compared with measured site data. A one-dimensional model to capture hydroacoustic effects along a pipeline will be presented. The corresponding PDEs were solved using an implicit finite difference scheme on a staggered grid. In contrast to previous studies this model is coupled to the commercial software ANSYS CFX through an interface which exchanges pressure and discharge data within every time step until convergence. Results of the one-dimensional approach as well as the coupled solution were validated with commercial one-dimensional software (SIMSEN) and a full threedimensional calculation for hydroacoustic test cases. Unlike former investigations the described 1D-3D approach is used to compare site data with a numerical analysis at prototype size focused on the amplitude and frequency of the pressure pulsation at overload condition. The combined model is able to capture the occurring phase change in the draft tube as well as the propagating pressure oscillation through the hydraulic system without solving for the whole penstock in a 3D manner, thus saving time and computational resources.


Introduction
In Francis turbines at operating conditions such as full load and overload the prevailing cavitation rope in the draft tube may act as an energy source, which leads to self exciting pressure oscillations propagating along the whole water conduit. The corresponding pressure pulsation frequency depends on load and tail water level. Typical frequencies are about 20%-100% of the runner rotation frequency [1]. The state of the art approach is to investigate the dynamic behaviour for the whole hydraulic circuit in a one-dimensional (1D) hydroacoustic manner. Therefore, all circuit elements must be replaced by their 1D representatives. The cavitation rope effect is represented by two lumped parameters, the cavitation compliance C = ∂V C /∂H and mass flow gain factor χ = ∂V C /∂Q as suggested by Brennen and Acosta [2] and adopted by several authors [3,4,5,6]. These parameters can be derived either from experiments or CFD calculations. However, the biggest drawback of this approach, beside some uncertainties in the definition of the lumped parameters, are missing insights in complex three-dimensional (3D) flow features. Doerfler et al. [1] used ANSYS CFX to investigate the transient behaviour depending on different numerical variations and boundary conditions, whereas Chirkov et al. [7]  Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under licence by IOP Publishing Ltd 1 or Cherny et al. [8] coupled a 1D approach representing the waterway with a three-dimensional in-house CFD code. The present work combines the advantages of a well validated, efficient and maintained commercial code to simulate distributor, turbine and draft tube in an unsteady 3D manner and attaches a 1D model to represent the penstock and spiral case. Prototype size model was considered because of existing site measurement data and to avoid uncertainties in model scale up. Benefits resulting from this work are full insights into the unsteady 3D flow through the turbine and a physical representation of the penstock to enable transient boundary conditions for the 3D computational model. A variation of numerical parameters allows to improve the consistency with the available test data regarding amplitude and frequency.

Hydroacoustic Modeling
The underlying mathematical model describing the dynamic behavior of a pipe element with length dx is derived from mass and momentum equation. Because of high wave speeds at low flow velocities, convective terms can be neglected with respect to propagative terms, which leads to eq (1).

discrete system
To solve the hyperbolic partial differential equation a discrete system is used. The pipe can be discretized in space using a centered scheme on a staggered grid as illustrated in Figure 1 [9]. Following Nicolet [10] an analogy to electric circuits can be drawn, where the model constants Figure 1: spatial discretisation are replaced by their electric representative.
Using these values and the centered scheme at location i + 1/2 yields to the implicit matrix form  The remaining temporal differential is integrated using a second order implicit Euler scheme: The resulting implicit matrix equation is solved with a Gauss-Seidel method and by applying an additional under-relaxation factor.

case study I: valve
The first test case considers a pipe with a fast closing valve at the downstream location. Due to the additional head loss across the valve, a loss term is added to the hydroacoustic resistance at the valve position.
The loss function is dependent on the type of the valve and in this case represented by the following function: where y is a power function according to eq (7).
The loss function has an asymptotic behavior and needs to be limited to an arbitrary high value. Figure 2 shows a comparison between the developed code and the commercial software SIMSEN for the test case under discussion. The time-dependent pressure at pipe midpoint is normalized with the head. Numerical setup regarding time step and grid size is identical. Time   discretization differs between the explicit fourth order Runge Kutta scheme (SIMSEN) and the implicit second order Euler backward scheme 1 . Besides an offset in closing time of the valve, which leads to a slight pressure overestimation, a good agreement can be stated.

case study II: tail water
A fast closing valve leads to large gradients and occurring back flow. Furthermore, the representation of the closing function implicates additional challenges for the 3D code. Therefore, a test case with a time-dependent tail water level is investigated. Figure 3 shows the model setup and the 1D-3D coupling. The commercial character of the 3D code impedes an implicit At every inner iteration loop the resulting averaged total pressure at inlet boundary is exchanged with the one-dimensional approach and a new discharge based on the 1D flow field is determined. Figure 4 shows the resulting pressure history for three different cases: 1D approach, compressible 3D calculation and coupled solution. Besides the frequency there are good agreements of all  three calculations regarding amplitude, period and damping. Deviations of the 3D approach are caused by a variable density based on absolute pressure distribution. This results in a slightly varying wave speed. For the coupled simulation the 3D pipe segment was short and modeled as being incompressible, therefore no deviations in the wave speed occurred and a good agreement could be achieved.
3. Francis turbine at overload condition 3.1. numerical setup The previously described and validated coupling is applied to a hydropower project. A Francis turbine of the upper power range with unbranched bare steel penstock and a draft tube with direct outflow to the tail water was investigated. During commissioning various measurements were carried out and lead to a good size database. For comparison with the numerical results an operating point according to OP#1 in Table 1 was chosen. The simulations were performed in prototype size taking the buoyancy into account to enable direct comparisons with prototype measurements and to avoid scale up uncertainties. The reference level for the buoyancy model was set to the measured tail water level. Considering the axisymmetric appearance of the vortex rope, it is assumed sufficient to model only one passage of stay vane, wicket gate and runner and applying periodical boundary conditions. The draft tube was modeled completely including elbow and piers. Transfer between stationary and rotating frame is realized by circumferentially averaging the flow variables (stage interface) or by transient rotor-stator coupling, accounting for all transient flow characteristics. To capture the dynamic interaction between the cavitation volume and discharge variations together with pressure oscillations, a multiphase simulation was carried out. Thus, unphysical vorticity peaks in the vortex center can be avoided. The mass transfer between the phases representing evaporation and condensation is described by the Rayleigh-Plesset cavitation model in terms oḟ where the model constant are as follows: Turbulent scales were modeled by the SST k-ω-model with automatic wall function. In general, ANSYS CFX 15.0 uses a finite volume based discretization scheme up to second order accuracy for convective fluxes and truly second order accuracy for diffusive fluxes. Time dependent computations were performed with second order accurate time scheme (eq. (4)).
The computational grid consists of three domains (tandem cascade, runner and draft tube) with block-structured hexahedral cells. The mesh size of SVWG (780,000) and RU (640,000) was kept constant for all investigations.

boundary conditions
Similar to the case study in chapter 2.3 an upstream total energy for the 1D penstock inlet boundary condition is assumed. Based on a measured pressure at the spiral case, pipe losses and an expected discharge according to the hill chart, the fixed total specific energy follows as: At the 1D-3D interface an averaged pressure extracted from the 3D stay vane inlet is set as a 1D pressure boundary condition. A discharge return value of the 1D model in addition with an appropriate inflow angle serves as the 3D domain inlet condition. As buoyancy was considered the hydrostatic pressure profile at the draft tube outlet was set in accordance to the measured tail water level. Figure 5a shows the computational domain in meridional view with the resulting vortex rope visualized by means of an isosurface corresponding to a volume fraction of 0.5. The yellow colored area represents the interface position between runner and draft tube and serves as evaluation plane to monitor the mass flow Q RU DT into the draft tube. Another evaluation plane is located at the end of the elbow to monitor the mass flow Q DT . Figure 5b shows contour plots of absolute pressure and circumferential velocity which are separated by the cavitation rope. The absolute pressure rise due to the hydrostatic effect is observable. Approaching the vortex rotational axis, the pressure drops in accordance with an increase in circumferential velocity which is plotted on the right-hand side. Locations where the cavitation volume exhibits contractions coincide with spots of maximum circumferential velocity. Inside the rope, the velocity field can be investigated through vectors colored according to their axial velocity (Figure 5c). In the first 'bubble', reverse flow occurs. High velocities at the phase boundary force the fluid inside the bubble downstream and enforces back flow to satisfy continuity. The reverse flow is causing the rope to contract.

transient simulation
Preliminary studies support the results of Doerfler et al. [1], showing strong numerical dependencies for transient two-phase flow with oscillating cavitation volume, in case of a coupled solution of momentum and pressure equations against the segregated solution of volume fraction on the continuity conservation between both phases. To quantify the mass conservation the ratio δ QV as introduced by Doerfler et al. [1] is used: with For very fine meshes and small time steps along with a good convergence of inner coefficient loops, this ratio can asymptotically reach values close to unity. In Figure 6, the mass flow  Figure 6: continuity error difference ∆Q, the cavitation volume V C and the temporal derivation of the cavitation volume dV C /dt are plotted over time. At the marked position, the draft tube grid #v01 was replaced by #v02 and the time step reduced accordingly. The influence of the spatial resolution in the draft tube and of the time step size on the ratio δ QV becomes obvious. After half a period at t = 0 s the coupling with the 1D-model was realized, corresponding to #v02 in Table 2. The ratio δ QV for the coarser grid and larger time step equals to 0.8 and could be increased to a value of 0.96 after the mesh refinement and a reduction of the time step. To further improve the mass conservation, Chirkov et al. [7] suggested a coupled solution algorithm for pressure and momentum equation together with the volume fraction, allowing coarser meshes and larger time steps (cf. #v03). In Figure 7 three different simulation settings according to Table 2 are compared by their pressure fluctuations in the draft tube. For variant #v01, a constant mass flow was prescribed at the inlet. The developed 1D hydroacoustic model was employed for variant #v02 and #v03. The model represents the penstock and the spiral case as a straight pipe with an overall length of 220 m. No measuring data of the wave speed in the penstock was available. From #v02 to #v03 the wave speed was reduced from initially 1368 m/s to 1000 m/s to account for penstock elasticity [11]. An initialization of the pipe is done with a linear pressure profile according to Bernoulli's formula.  Variant #v01 shows an increase in amplitude over time corresponding to a forced oscillation with negative damping. Therefore no stable condition can be achieved. Variant #v02 and #v03 are damped by the oscillating inlet mass flow and associated losses, resulting from the penstock modeling. At the current simulation state, the damping still exceeds the selfexcitation of the hydraulic machine and leads to a decay in amplitude. To quantify the decay behavior an exponential fitting function p DT (t) =p DT,∞ + A 0 e −R 0 t was applied for both variants. Hence, the stable state can be estimated by extrapolation. As a result of the changes applied in the setting for variant #v03 the decay is significantly decreased (R 0,#v02 = 0, 121 → R 0,#v03 = 0, 068). The oscillation for #v02 is expected to be damped out almost entirely, whereas the amplitudes of #v03 will reach a constant value of approximately 0.52 · 10 4 Pa. Thus #v03 is used for a comparison with the measured data below. In addition, frequency increase due to the 1D-3D coupling for both variants is observed.

comparison with measurement data
The comparison with site measured static pressure is done at two positions. Figure 8 shows the pressure oscillation in the draft tube measured at the man door together with the simulation results obtained in the 3D simulation domain. Figure 9 compares the pressure at the spiral case intake with the simulated pressure in the corresponding 1D pipe segment. The diagrams show only the dynamic portion for better comparability. The pressure measured at both positions oscillates synchronously and is reproduced well by the 1D-3D simulation. However, the frequency of 0.664 Hz is underestimated compared to the measured frequency of 0.8125 Hz. The amplitudes in the draft tube and at spiral case intake decay below the measurements.  Figure 8: dynamic pressure at draft tube measurement location

Conclusion
A 1D hydroacoustic model was presented and in contrast to existing approaches the coupling with a commercial CFD code is realized. The hydroacoustic model and the coupling were validated with two test cases. It was shown that the 1D model is able to represent the hydroacoustic behavior of a pipe and the coupling captures the interaction between the 3D flow field and the 1D hydroacoustic model. Furthermore, the hybrid 1D-3D approach was used for the investigation of an overload instability in a Francis turbine. Distributor, turbine and draft tube were simulated by means of two-phase CFD, whereas penstock and spiral case were represented with the 1D hydroacoustic model. In addition to previous studies, simulations were carried out in prototype size to enable direct comparison with site measuring data and avoid scale-up uncertainties. Different numerical settings were described, showing strong dependencies on solution strategy, mesh density and time step. In comparison with a constant mass flow boundary condition, using a 1D penstock model generates a variable mass flow into the 3D domain, resulting in a damped oscillation and an increasing frequency. An extrapolation of the simulation results indicates that the oscillation will reach a stable state. However, deviations to the measurement still remain with respect to frequency and amplitude. It is expected that the deviations may be further reduced by a 360 • simulation, inhomogeneous two-phase calculation or considering the spiral case geometry.