Flow-induced pulsations in Francis turbines during startup - A consequence of an intermittent energy system

Hydraulic turbines are increasingly responsible for regulating the electric grid, due to the rapid growth of the intermittent renewable energy resources. This involves a large increase in the number of starts and stops, which cause severe ﬂ ow-induced pulsations and ﬂ uctuating forces that deteriorate the machines. Better knowledge of the evolution of the ﬂ ow in the machines during transients makes it possible to avoid hazardous conditions, plan maintenance intervals, and estimate the costs of this new kind of operation. The present work provides an in-depth and comprehensive numerical study on the ﬂ ow-induced pulsations and evolution of the ﬂ ow ﬁ eld in a high-head model Francis turbine during a startup sequence. The ﬂ ow simulation is carried out using the OpenFOAM open-source CFD code. A thorough frequency analysis is conducted on the ﬂ uctuating part of different pressure probes and force components, utilizing Short-Time Fourier Transform (STFT) to extract the evolution of the frequency and amplitude of pulsations. Low-frequency oscillations are detected during the startup, which are induced by the complex ﬂ ow structure in the draft tube. A decomposition is performed on the draft tube pressure signals, and the variations of the synchronous (plunging) and asynchronous (rotating) modes are studied. The plunging mode is stronger at minimum and deep part load conditions, whereas the rotating mode is dominant during the presence of the Rotating Vortex Rope (RVR) at part load. The velocity ﬁ eld in the draft tube is validated against experimental data, and the complex ﬂ ow structures formed during the startup procedure are explained using the l 2 vortex identi ﬁ cation method. © 2022 The Author(s). Published by Elsevier Ltd. This is


Introduction
The use of renewable electric energy resources has been growing fast to respond to the increasing global electric energy consumption. Nowadays, the inevitable intermittency of electrical energy resources such as solar and wind power is compensated through hydropower systems [1e3]. The hydraulic turbines are not necessarily working at the steady Best Efficiency Point (BEP) condition anymore. They are being used in different transient operating sequences to stabilize the electrical grid, leading to entirely different engineering requirements for such machines.
Transient operations usually produce complex flow structures, such as flow separation, vortices, destructive pressure pulsations, cavitation, etc. Frequent occurrence of such undesirable flow structures could seriously deteriorate the turbine lifetime and cause fatigue stresses, wear and tear on different components [4]. Currently, Francis turbines may experience over 500 start-stop cycles per year [5], while they are usually designed to tolerate up to 10 cycles [6,7]. Undoubtedly, the accumulated damages from such abundant cycles degrade the machine's performance and may lead to its failure. Hence, it is crucially important to study and provide a profound understanding of the turbine flow field during transient operations such as startup.
Gagnon et al. [8] examined the influence of startup schemes on the fatigue-based life expectancy of a Francis turbine. It was explained that an optimization of the scheme could improve the turbine lifetime. Nicolle et al. [9] assessed the startup operations of a low-head Francis turbine using a reduced CFD model. Two different startup scenarios based on the guide vane opening scheme were investigated. Comparisons were made with limited experimental measurements and a general agreement was achieved.
The impact of the guide vane opening scheme on the startup procedure of a high-head Francis turbine has been experimentally assessed by Trivedi et al. [10]. The angular speed of the guide vanes was for one scheme almost twice as for the other scheme. Inappropriate rapid rotation of the guide vanes amplified the unsteadiness and developed undesirable pressure pulsations. Goyal et al. [7] performed an experimental study on the same high-head Francis turbine during startup. The startup sequence was split into two phases, namely, phase I, to synchronize the turbine with the generator, and phase II, to reach the steady-state condition. The second phase was accomplished using three different guide vane opening schemes which ended in Part Load (PL), BEP, and High Load (HL) conditions, respectively. The Rotating Vortex Rope (RVR) frequency was observed in both velocity and pressure data of the first scheme.
More recently, the startup of a prototype Francis turbine was experimentally and numerically investigated [11]. Two guide vane opening schemes, namely, conventional and reduced opening limit schemes were studied and it was shown that the reduced scheme decreased the fatigue damage. The draft tube vortices were shown to have a significantly higher impact on the dynamical stresses compared to the interblade vortices. It was concluded that disturbing the draft tube vortex could alleviate the damaging effects on the runner during startup.
Although the experimental investigations are trustworthy resources to assess the turbine flow field during startup, they are expensive and there are many limitations on accessibility and measured flow details. Numerical studies provide a reliable addition to assess and understand the details of the flow field during turbine startup. The startup is recognized as one of the most harmful operating conditions of hydraulic turbines [8]. Therefore, achieving a profound understanding of the complex flow field of a hydraulic turbine during startup is essential to reduce the damaging effects and improve the life expectancy of these machines.
The present article provides a comprehensive and detailed analysis of the transient flow field and its pulsations during a startup sequence of a Francis turbine. Such in-depth analyses are crucial for a better understanding of the hazardous pulsations to be able to ultimately reduce and avoid them. The simulation is performed utilizing the OpenFOAM open-source CFD code. The variation of the pressure field, velocity field, and forces are carefully assessed. One of the main focuses of the present study is to extract the flow-induced pulsations. No investigations are found in the literature on the draft tube pressure signal decomposition during a startup sequence. In the current study, for the first time, the variation of plunging and rotating modes of the fluctuating pressure during the startup operation is examined. An in-depth explanation of the complex flow structures downstream the runner, which play a crucial role in the generation of the pulsations, is presented. The paper is organized as follows. The investigated test case, including the geometrical and operational details, is introduced in Section 2. The mathematical formulations of the assessed problem are described in Section 3.1, while the details of the numerical framework are described in Section 3. Section 4 provides the numerical results and discussions, and finally, the concluding remarks of the paper are provided in Section 5.

Investigated test case
A high-head Francis turbine model is used as the investigated test case. The Francis-99 turbine model, provided by the Francis-99 workshop series [12], is a 1:5.1 scale model of a prototype Francis turbine [13]. The runner consists of 15 full-length and 15 splitter blades. The prototype and model net heads are about H prototype z 377 m and H model z 12 m, respectively. Fig. 1a and b show two cross-sections of the Francis-99 model.
The axial and horizontal in-plane velocity components have been experimentally measured at a PIV plane. The PIV plane is shown as a red line and a gray shaded area in the z-normal and y-normal sections, respectively. The velocity measurements are reported on three PIV lines, two horizontal lines (Lines 1 and 2), and one axial line (Line 3). Moreover, the static pressure is reported for three sensor locations, namely, VL2, DT5, and DT6. In the experiments, the draft tube pressure sensors were piezoelectric, and only instantaneous fluctuation of pressure was measured [14]. Two additional numerical probes (RP1 and RP2) are defined in the rotating zone (Runner) to sample the pressure field throughout the sequence. The numerical probes are placed in the middle of one runner passage (in between two neighboring main and splitter blades) at different axial positions. The current work concerns a startup sequence that commences from the minimum load operating condition. The guide vanes are nearly closed with an opening angle of a ¼ 0.8 and the flow rate is

Computational framework and numerical aspects
The CFD simulation is carried out with OpenFOAM-v1912 [15,16]. The governing equations are discretized using the finitevolume approach on a collocated mesh. The current section briefly describes the governing equations and the employed numerical methods and schemes. More detailed information about the numerical aspects of the performed CFD simulation is provided by Salehi et al. [17], who used the same approach for a shutdown sequence of the same case.

Mathematical formulation
A transient incompressible turbulent flow can be modelled by the Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations, given by where Àru i u j represents the unknown Reynolds stress tensor. The Shear Stress Transport (SST) based Scale-Adaptive Simulation URANS model (i.e., SST-SAS) [18,19] is here employed for the calculation of the Reynolds stress tensor. SST-SAS is a turbulenceresolving URANS model, used for simulations of industrial transient flows. Its formulation decreases the local eddy viscosity to resolve the turbulent spectrum and break-up of large eddies, providing LES-like solutions. Several research studies verified the performance of the SST-SAS model in the simulation of hydraulic machinery flows [11,20e25].

Discretization schemes
The second-order backward implicit scheme is employed for the discretization of the temporal derivative terms. The time step of the simulation is chosen as Dt ¼ 1.25 Â 10 À4 s, corresponding to runner and guide vane rotations of 0.25 and 1.625 Â 10 À4 in each time step. The average and maximum CFL numbers at the highest flow rate (BEP) are 0.025 and 55. It should be noted that the CFL number is less than 2 for 99.4% of the cells.
The convective terms in the momentum equation are discretized using the Linear-Upwind Stabilised Transport (LUST) scheme [26], which blends the central and second-order upwind schemes with a blending factor of 0.75. In other words, the face values are calculated blending 75% second-order central and 25% second-order upwind schemes, balancing accuracy and numerical stability. The second-order upwind scheme approximates other convective terms (i.e. in k and u equations).
The Laplacian terms in the transport equations are estimated using the second-order central scheme. An explicit non-orthogonal correction due to the high skewness of the cells at some locations is inevitable because of the complex geometry.

Pressure-velocity coupling
The PIMPLE pressure correction algorithm is employed for the pressure-velocity coupling. It combines two pressure correction algorithms, namely, SIMPLE [27] and PISO [28] as outer and inner correction loops, respectively. A maximum of 10 outer correction loops is performed in each time step, controlled by a residual criterion. At most time steps the flow solution is converged after four outer correction loops. Each outer loop conducts two inner correction loops. After each inner loop, one additional nonorthogonal correction loop is performed to assure convergence of explicit terms. It has been shown that the OpenFOAM implementation of the pressure correction algorithm is in line with the Rhie-Chow interpolation technique [29,30].

Boundary conditions
The guide vanes open up by rotating with a constant rotational speed of 1.3 /s. The time-variation of the guide vane opening angle is plotted in Fig. 2a. As seen in the figure, a smooth transition is implemented at the start and stop of the rotation (t ¼ 2 and t ¼ 9 s) to minimize the numerical instability caused by the sudden movement of the guide vanes. The total time of the sequence, t ¼ 12 s, corresponds to 66.52 runner revolutions.
The guide vane movement is imposed through an ad-hoc developed boundary condition that requires the guide vane rotational speed as input. Therefore, the rotational speed u of the guide vanes is shown in Fig. 2b.
It is assumed that the inlet volume flow rate of the turbine varies linearly with respect to the guide vanes angle. This assumption is according to the Francis-99 workshop series recommendation due to inaccurate measurements of the flow rate during transient operation [12]. Hence, a time-varying spatially uniform velocity, according to the flow rate, is imposed at the inlet of the spiral casing. A fixed turbulence intensity (I ¼ 7%) and viscosity ratio (n t / n ¼ 100) is considered for the inflow condition. The inlet pressure is extrapolated from the inside domain using a zero-gradient assumption. All quantities at the outlet boundary are computed using the zero-gradient condition, except the pressure which is set by a fixed value.
As previously described, there are four different mesh regions in the simulation (spiral casing, guide vanes, runner, and draft tube). The Cyclic Arbitrary Mesh Interface (cyclicAMI) [31,32] was utilized to transfer the information between the different domains.
In order to reach a statistically stationary state at minimum load condition, the flow is solved for 4 s flow time corresponding to over 22 runner rotations, and then the startup sequence presented in Fig. 2 is initiated.

Dynamic mesh framework
CFD analysis of the transient operation of Francis turbines includes two types of simultaneous mesh motion, i.e, mesh deformation of the guide vane domain due to the rotation of each guide vane and solid body rotation of the runner domain. Therefore, a Laplacian displacement mesh morphing solver is employed to deform the guide vane domain mesh while the solid-body rotation function handled the runner rotations. In each time step, the mesh is updated at the beginning of the first PIMPLE outer correction loop. Then, the face fluxes are calculated based on the face swept volumes and relative fluid velocity [33,34].
The mesh morphing is governed by a Laplace equation, given by where G is the motion diffusivity and d cell is the displacement vector of the cell centers. The Laplace equation is solved for the cellcentered displacement (d cell ) and then the solution is interpolated to get the point displacements (d points ). Finally, the new point locations (at time t þ Dt) are simply computed as The motion diffusivity (G) is obtained using a quadratic inverse distance scheme with respect to the guide vane surfaces.
Severe mesh deformation of the guide vane region due to the large rotation of the guide vanes during the startup sequence could potentially result in low-quality mesh cells and consequently deterioration of convergence and accuracy of the numerical results. Therefore, in this study, the mesh quality parameters were monitored during the mesh deformation. The guide vane region was remeshed two times at guide vane openings of a ¼ 3.47 and a ¼ 6.85 to maintain an acceptable mesh quality. More information on the numerical aspects and mesh deformation, as well as the open-source case and codes of the current study, is provided by Salehi and Nilsson [35], as the same case and codes are employed in the present work.
A block-structured mesh is created for the CFD simulation. The mesh at BEP contains a total of 16 million cells (for more information please see our previous studies [17,35]).

Parallel processing
The scotch [36] domain decomposition approach is used to split the computation domain and distribute roughly equal loads to the processors while minimizing their interconnections. The job is submitted to a Linux cluster using 320 CPU cores. The full startup sequence consumed a computational cost of 170,000 core hours.

Results and discussion
This section presents the results of the transient startup sequence of the Francis-99 model turbine.

Pressure fluctuations
As previously described, a number of pressure probes, namely VL2, DT5, DT6, RP1, and RP2, were defined in the computational domain (see Fig. 1), and the variation of the static pressure is recorded throughout the entire startup sequence. The experimental results of the static pressure are available for the VL2 probe, while only the pressure fluctuations were monitored at DT5 and DT6. In transient (time-varying) turbulent flows, as in the current case, the obtained signals (for instance, pressure) consists of two different parts, the mean and fluctuating parts ðp 0 ¼ p À p̄Þ. The mean signal changes through time due to the variation of the operating condition. Therefore, in order to extract the fluctuating pressure, the instantaneous mean should be calculated. The present study employs the Savitzky-Golay finite impulse response filter [37] for smoothing the obtained signals and calculating instantaneous mean and fluctuations. A variable window size is chosen to capture the fluctuations more accurately. The window is much smaller at the start and end of the transient sequence, where the variation of the pressure level due to the change in operating condition is sharper. Fig. 3 shows the time-variation of the static pressure, its instantaneous mean, and the fluctuations of the static pressure from its instantaneous mean. In general, the numerical prediction of the VL2 pressure ( Fig. 3a) sufficiently matches the experimental data, although with slightly lower values at the BEP condition at the end of the sequence. The maximum relative error, calculated as | p num À p exp |/p exp Â 100, is 4.25%. Each plot contains a zoomed view that covers a 90 rotation of the runner in either the stationary minimum load or BEP condition. The VL2 zoomed views show clear smooth pressure pulsations due to the Rotor-Stator Interaction (RSI) in the vaneless space (between the runner blades and the guide vanes). Since the runner consists of 30 full and splitter blades, 7.5 pressure pulsations can be seen in these zoomed views. The vaneless space static pressure fluctuates around a nearly constant mean pressure at the minimum load condition. Some lowfrequency oscillations are also visible at the minimum load condition in the VL2 pressure, which could be due to large unsteady flow structures in the massively separated flow in the draft tube. When the startup sequence commences at t ¼ 2 s, the guide vanes start opening up, and consequently, the pressure increases with the turbine flow rate growth. The rate of the pressure increment is initially higher and then it reduces and reaches a constant level until the end of the sequence. The numerical results suggest an overall pressure rise from 160 kPa (at minimum load) to 174 kPa (at BEP). The numerical results reach a stationary condition at t ¼ 9 s when the sequence finishes. In contrast, the experimental pressure results show that the flow still needs some time to reach the steady condition, due to dynamics in the experimental open-loop hydraulic system. Some low-frequency oscillations are also visible in the numerical pressure results after the initiation of the sequence. These oscillations are most likely produced by large flow structures formed in the draft tube in the low load conditions and will be discussed in detail later. One can see such pulsations more apparently in the fluctuating pressure shown in Fig. 3b. Distinct periodic oscillatory patterns are seen between t ¼ 4.5 s and t ¼ 6.5 s that are probably caused by the formation and diminish of the RVR. The static and fluctuating pressure in one of the draft tube probes (DT6) is also shown in Fig. 3c and d. There is not a clear sign of the RSI fluctuations in the presented zoomed views at BEP. Here again, There are two statistically stationary phases in the whole simulated sequence, namely, the initial minimum load and the final BEP conditions. The Fast Fourier Transform (FFT) analysis technique enables us to identify the excited frequencies and their amplitude of the obtained signals. Therefore, FFT was applied on the fluctuating part of the VL2 and DT6 pressures at both stationary conditions, and the results are plotted in Fig. 4. It should be noted that in the present study, all the frequencies are normalized by the runner rotational frequency f n ¼ 5.543 Hz. The runner blade passing frequency f b ¼ 30f n (15 full-length blades and 15 splitter blades) is the dominant frequency in the VL2 probe. The amplitude of f b is much larger at BEP compared to the minimum load condition. Peaks are visible at the harmonics of the runner passing frequency (15f n and 60f n ). A low frequency of approximately 0.3f n is also excited at the minimum load condition of the VL2 pressure, which could be explained by the large separated flow region in the draft tube at such conditions. The draft tube fluctuating pressure (DT6) seems to be most excited at the above-mentioned low frequency. A moderate peak can also be seen at the frequency of 15f n , corresponding to the full-length blade passing frequency, as only full-length blades are elongated to the draft tube. In other words, the DT6 draft tube probe can sense the rotation of the full-length blades much more than the splitters.
Due to the time-varying nature of the obtained signals in transient sequences, such as turbine startup, both the excited frequencies and their amplitudes change throughout the sequence. Hence, a Short Time Fourier Transform (STFT) analysis is required for time-frequency analysis. STFT divides the full-time domain into small subdomains and performs the Fourier transform on each subdomain. The time-variations of the amplitudes of different frequencies of the VL2 and DT6 fluctuating pressures are illustrated as spectrograms in Fig. 5. The runner blade passing frequency (f b ¼ 30f n ) is the dominant frequency of the vaneless space pressure throughout the whole sequence (Fig. 5a). The harmonic frequencies (i.e., 15f n , 45f n , 60f n , 75f n , etc.) are also clearly excited. A wide range of excited stochastic frequencies are visible in the minimum load condition (t < 2 s), indicating a complex flow field including large separations and vortex breakup. When the guide vanes start to open up and the flow rate increases, such frequencies diminish slightly after t ¼ 2 s. The zoomed view of the VL2 spectrogram suggests the existence of low-frequency high-amplitude oscillations during the transient sequence. The RVR phenomenon is most likely responsible for such types of pulsations. The DT6 spectrogram denotes a deterministic frequency of 15f n , corresponding to the passing of the full-length blades. The RVR low-frequency oscillations are also clearly visible here.
The time-variation of the amplitude of different excited frequencies is extracted from the STFT calculations and presented in Fig. 6. For the VL2 sensor, the runner blade passing frequency (30f n ) is dominant throughout the whole sequence. The amplitude is nearly constant and slightly increases when the turbine reaches the BEP condition. The amplitude of the RVR frequency (0.3f n ) is increased with the initiation of the transient sequence and then decreases as the large RVR structures diminish when the turbine approaches the BEP condition. On the other hand, for the DT6 probe, the amplitude of 0.3f n is dominant in the entire sequence, except for the BEP condition where the large draft tube vortical structures are washed away. A sudden rise and then decrease is observed in the amplitude of 0.3f n in the middle of the sequence due to the formation and collapse of the RVR.
Hydraulic turbine draft tube cones generally experience two different types of pressure pulsations at low load conditions [38,39]. The pressure signals can be decomposed into synchronous and asynchronous modes. The synchronous mode (also known as the plunging mode) is somehow similar to the water hammer pressure waves which travel throughout the whole hydraulic system. The asynchronous mode (rotating mode), produced by the local instabilities such as the RVR, is only active in the crosssections. The pressure signal decomposition can be performed using the unsteady signals of two different pressure probes which are positioned at the opposite sides of the draft tube cone with the same height, through ðSynchronous component or plunging modeÞ; ðAsynchronous component or rotating modeÞ: A few researchers have studied the draft tube pressure signal decomposition to identify the appearance of plunging and rotating modes at low load conditions of hydraulic turbines (e.g., Refs. [38e42]). However, no investigation can be found in the literature on the decomposition of plunging and rotating modes of a hydraulic draft tube during a startup sequence. Extracting such modes from pressure signals can be particularly helpful for explaining the appearance and collapse of the RVR in transient sequences like shutdown and startup.
In the present test case, the DT5 and DT6 sensors are placed on opposite sides (180 apart) of the conical part of the draft tube and could be used for signal decomposition. First, the synchronous and observed. The sudden rise of rotating mode after t ¼ 4 s could be a sign of the formation of rotating vortical structures (i.e., RVR) which decays with further increasing turbine load after t ¼ 6 s. After t > 7 s the turbine approaches the BEP condition and the large vortical flow structures inside the draft tube cone vanish, the pressure fluctuations predominantly contribute to the synchronous mode, and the asynchronous mode is rather negligible at the design condition.
In order to better understand the formation and collapse of the RVR and its impact on the decomposed pressure modes, a bandpass filter with a narrow frequency range of 0.1f n , centered at the fundamental frequency of RVR (0.3f n ), was applied to the decomposed signals to isolate the RVR effects in the plunging and rotating modes. As previously seen in Figs. 4 and 5, the frequency of 0.3f n is  the dominant frequency inside the draft in a wide range of low-load conditions. The filtered signals displayed in Fig. 7b reveal that the plunging effect is mostly the dominant mode at minimum load and deep part load conditions (before t ¼ 4 s), suggesting that disintegrated stochastic flow structures at such conditions primarily cause axial pulsations that are sensed throughout the whole system at the same time. Nonetheless, when the startup sequence of the turbine initiates the rotating effects gradually increase in time while the plunging mode weakens. The fact that the rotating mode is dominant between t ¼ 4.5 s and t ¼ 6.5 s could be a clear sign of the formation and collapse of the RVR. As expected, no large vortical structures should exist at the BEP condition. Therefore, both the plunging and rotating modes decay after t ¼ 7 s. A time-dependent frequency analysis was performed on the decomposed signals and the results are shown as spectrograms in Fig. 8. Here again, synchronous fluctuations are observed at minimum load and deep part load conditions (t < 4 s), while asynchronous pressure pulsations can be detected between t ¼ 4.5 s and t ¼ 6.5 s. More specially, the fundamental frequency of the RVR is much more pronounced in the rotating mode than the plunging mode during the presence of RVR (between t ¼ 4.5 s and t ¼ 6.5 s), as also pointed out by Goyal et al. [42].
As explained in Section 2, two probes are defined in the rotating domain of the runner, namely, RP1 and RP2 (see Fig. 1b) and their pressure variation throughout the startup sequence is demonstrated in Fig. 9. Predictably, the RP1 pressure is generally higher than that at RP2, as it is closer to the runner inlet. Both pressure probes exhibit a gradual rise during the transient sequence. The RP1 pressure increases by 13.3 kPa, whereas the RP2 pressure grows by 7.2 kPa. High-frequency RSI fluctuations are visible through the provided zoomed views. Here the probes are rotating with the runner and therefore the pressure is expected to show a peak whenever the probe is passing a guide vane trailing edge. The fluctuating part of the pressure shows stronger high-frequency RSI fluctuations for RP1 as it is closer to the guide vanes. Both probes contain low-frequency oscillations which are slightly amplified between t ¼ 4.5 s and t ¼ 6.5 s. Fig. 10 plots the FFT of the fluctuating pressure of the rotating probes at the stationary conditions (minimum load and BEP). As expected, the fluctuations have a dominant frequency at the guide vane passing frequency (f gv ¼ 28f n ) which is stronger at BEP. The first harmonic of this frequency (f gv ¼ 56f n ) also shows a small peak. Additionally, some low-frequency peaks, due to the formation and breakup of vortical flow structures are detected by the FFT analysis at minimum load conditions. Since the probes are rotating, the runner rotation frequency (f n ) and its first few harmonics (2f n , 3f n , 4f n , etc.) are also excited at both conditions. An STFT analysis can further explain the variation of the amplitudes of the excited frequency during the transient sequence. Fig. 11 presents rather similar trends for the time-variation of the amplitudes for both rotating probes. The guide vane passing frequency is a deterministic and dominant frequency during the whole sequence. The zoomed views (Fig. 11b and d) denote that at minimum load condition, a vast range of stochastic frequencies is excited which decay after a short while into the transient sequence. Although an excited low frequency is observed during the formation of the RVR (between t ¼ 4.5 s and t ¼ 6.5 s), the value of that frequency is larger than the RVR fundamental frequency. When the turbine reaches the design condition, the excitation of the runner rotation frequency (f n ) and its harmonics are clearly visible in the zoomed view of both probes.

Force pulsations
Sharp variations and oscillations of forces and moments exerted on different parts of hydraulic turbines during transient operations could cause serious damages and negatively affect the lifetime of the turbine. Therefore, performing force analysis during transient sequences like the startup is essential for mitigating such damaging effects. Forces and moments acting on the runner surfaces (i.e., hub, shroud, main blades, and splitters) as well as a single guide vane are monitored during the startup sequence and the results are presented in this section. Although the runner force analysis is performed on the whole runner in the present work, investigating the fluctuating forces on one individual runner blade is suggested for future studies as it could be beneficial for assessing the fatigue effects and lifetime. Fig. 12 shows the x and z components of the force acting on the runner, as well as the runner torque (z component of moment  vector). F x shows strong low-frequency oscillations between t ¼ 4.5 s and t ¼ 6.5 s, in which the vortex rope is formed and rotates around the turbine axis (Fig. 12a). At the BEP condition, the x-force is fluctuating around a non-zero value, indicating that the flow distribution around the runner is not perfectly axisymmetric. The zoomed view of the fluctuating part of F x (Fig. 12d) denotes both high and low-frequency oscillations at the BEP condition. The axial force (z-force) is initially oscillating at negative (downward) values at the minimum load condition and shortly after the commencement of the transient sequence it increases and becomes positive (upward) and continues its growth until the BEP condition (Fig. 12b). Comparing the fluctuating parts of F z and F x suggests that the low-frequency oscillations during the formation and collapse of the RVR are much weaker for F z (compared to their corresponding instantaneous mean). In other words, the RVR mostly affects the horizontal (radial) forces rather than the axial force. This is compatible with the signal decomposition analysis presented in Section 4.1. The axial forces are expected to oscillate with the plunging mode of the RVR, while the radial forces vary with the rotating mode. As elaborated in Fig. 7b, the rotating mode of the RVR is the dominant mode during t ¼ 4.5 s and t ¼ 6.5 s, and thus the radial force oscillations are greater. The variation of the runner axial torque through time exhibits a smooth linear growth in absolute value of the torque with turbine load increase, from 29.8 N m at minimum load to 630.3 N m at BEP. It is also seen that the fluctuating part of the torque signal is negligible with respect to its instantaneous mean. Fig. 12f reveals that the formation of the RVR barely affects the fluctuating torque and the maximum fluctuating  torque that occurs at the minimum load condition is less than 1.5 N m. The spectrograms in Fig. 13 exhibit STFT analysis of the runner forces. The runner blade passing frequency is a deterministic dominant frequency throughout the entire sequence for both the horizontal and axial forces (Fig. 13a and c). The f b frequency is less isolated in the f z force and it is more affected by the wide range of stochastic frequencies. The zoomed view of the F x spectrogram (Fig. 13b) shows mainly stochastic frequencies at minimum load that vanish with the initiation of the sequence. Then the impacts of the RVR on the horizontal forces are clearly seen as low-frequency oscillations. However, the axial forces present a wide range of excited low-frequencies that are not limited to the formation or the collapse of the RVR (Fig. 13d). Here again, we can deduce that complex flow structures at deep part load have strong plunging effects and result in fluctuations of the axial force. The runner rotation frequency (f n ) has an important role in the variation of the horizontal forces at the BEP condition.
To further assess the variation of forces in the Francis-99 startup sequence, the forces and moments acting on a single guide vane are studied. Fig. 14 depicts the time-variation of the radial force and torque (axial moment around the guide vane rotational axis) of the guide vane nearest to the volute tongue. The negative torque acts as to open the guide vane and vice versa. Both plots display smooth variations during the transient sequence with nearly constant ranges of RSI fluctuations, which are larger for the radial force. The F r is maximum at minimum load and reduces with load increase. More importantly, during the formation and collapse of the RVR (between t ¼ 4.5 s and t ¼ 6.5 s), F r oscillates with some lowfrequency oscillation but M z does not show any impact from the RVR. The spectrogram of F 0 r , illustrated in Fig. 15, demonstrates a broad span of stochastic frequencies at the minimum load condition. This could be the impact of complex separated flow structures formed behind the trailing edge of the guide vanes at minimum load condition. Fig. 16 employs an iso-surface of l 2 ¼ 7500 s À2 to reveal these structures. As expected, the 0.3f n frequency is distinctly evident during the existence of the RVR in Fig. 15.

Velocity variation
The velocity field is sampled through the entire startup sequence along the three lines shown in Fig. 1 and the numerical results are compared to the experimental data for validation. The variation of the flow field is carefully examined to understand and explain the draft tube flow field during the turbine startup. It should be mentioned that in this work, the horizontal velocity (U) represents the velocity component parallel to Line 1 and 2 (similar but not identical to radial velocity), while the normal velocity (V) is the velocity component normal to the PIV plane (similar but not identical to tangential velocity).
Figs. 17e19 present the time-variation of the numerical velocity components along the three PIV lines (previously shown in Fig. 1). The axial and horizontal velocity components are compared to the experimental measurements. The variable s represents the curve length of each line, which is normalized by its maximum in all plots. The comparison reveals that the numerical axial velocity (W) trend is quite similar to the experiment and thus is adequately well predicted by the simulation. At minimum load condition, the axial velocity direction is upward all over both Lines 1 and 2, varying with low-frequency oscillations. This indicates a massive reversed flow region that covers the entire extent of both lines, while the small mass flow through the draft tube cone passes outside those lines. Then, when the guide vanes start to open up at t ¼ 2 s, the reversed flow region gradually gets smaller. The low-frequency oscillations amplify with the establishment of the RVR. After reaching the BEP condition, the reversed flow region completely vanishes and the flow is entirely in the downward direction. At the design condition, the magnitude of W increases with the distance     to the draft tube walls, whereas it is slightly decreased at the center (s/s max ¼ 0.5) due to the runner cone wake. The fluctuations around the center are induced by the vortex shedding created behind the runner cone. The W contours on Line 3, which is located at the center of the draft tube, denote the existence of the same reversed flow region which is diminished just before reaching the design condition. Both the numerical and experimental data show a rather sharp change in velocity direction around t ¼ 8 s. Fig. 18 indicates that the horizontal velocity mainly fluctuates around zero in both the numerical and experimental data. The fluctuations are initially moderate at minimum load condition and clearly magnify during the formation and collapse of the RVR (between t ¼ 4.5 s and t ¼ 6.5 s). Then, at the design conditions, both the numerical and experimental results show near-zero U values with small fluctuations.
The normal velocity component (V) was not measured in the experimental study. Therefore, Fig. 19 displays only the numerical results of the time-variation of the normal velocity. A strong swirling flow exists at minimum load condition. The commencement of the startup sequence and load increase temporarily reduces the normal velocity. However, the V component remarkably increases and oscillates with the creation of the RVR. After the collapse of the RVR, the V velocity smoothly reduces. Suddenly before the steady BEP condition, the direction of the V velocity and consequently the swirl orientation changes and a weak counterrotating flow exists at the design condition.
The time-variation of the velocity field is further assessed using two points, namely Point 1 and 2 (see Fig. 1 To examine the swirling flow in the draft tube, the normal velocity (V) is presented in Fig. 20 as well, although no experimental data is available. Positive values of V at Points 1 and 2 indicate a swirling flow in the same direction as the runner rotation, and vice versa.
Water turbines are designed such that a nearly non-swirling flow leaves the runner at BEP. The presence of a weak swirling flow at the design condition could help the flow to stay attached to the draft tube walls. A residual positive (in the same direction of the runner) swirling flow exists at partial load condition, while a higher flow rate than design condition (high load) forms a negative (counter-rotating) residual swirl.
As expected, at minimum load condition a considerable positive tangential component exists, especially on Point 2 which is further from the draft tube center, indicating a large remaining positive swirl in the draft tube. Similar to U and W, the V velocity at Point 2 experiences the large oscillations of the RVR sooner than at Point 1. When the large rotational RVR structures in the draft tube are diminished, the V velocity decreases as the load increases. It settles at an insignificant negative value, indicating a weak counterrotating swirl at the design condition which could be intentional to keep the flow attached to the draft tube walls.

Flow structures in the draft tube
The formation and breakup of the vortical flow structures inside the draft tube during the startup sequence is analyzed in this section. Based on a previous study [17], the l 2 -criterion is employed to identify and visualize the vortical flow structures. It assumes a vortex to be a region with two negative eigenvalues of the S 2 þ U 2 tensor [43], where S and U are the strain and rotation tensors, given by S ¼ Therefore, a vortex can be identified as a region with a negative second largest eigenvalue, l 2 . The OpenFOAM function object Lambda2 changes the sign of the S 2 þ U 2 tensor eigenvalues, and thus a positive value should be used for the creation of the l 2 isosurfaces.  A video is supplied with the article for the readers to see the time-evolution of the vortical structures during the startup operation. Fig. 21 utilizes an iso-surface of l 2 ¼ 750 s À2 to unveil the evolution of the draft tube vortical structures by 12 snapshots during the turbine startup sequence. The corresponding times (t), guide vane opening angles (a), and the turbine load (normalized flow rate Q/Q BEP ) are denoted below each figure. The transient sequence starts from the minimum load condition. At this condition, a massively separated flow field with a significant residual positive swirl exists downstream of the runner (Fig. 21a). As a result, large persistent vortical structures are visible upstream of the draft tube elbow that produces low-frequency pulsations in the flow field and turbine forces. When the transient sequence initiates at t ¼ 2 s, the guide vanes start to open up and the turbine load increases. Consequently, the growing flow rate washes down the large vortical structures (Fig. 21b and c).
At t ¼ 4.2 s (Fig. 21d) the aforementioned large vortices are completely vanished and instead elongated vortical structures are formed downstream the runner. These are formed due to the instability of the shear layer between the swirling downward and separated upward flow regions. The separated region is still quite large and thus the shear layer is close to the draft tube wall. Four distinct draft tube vortices are formed in this region. Continuing the startup sequence, the vortical flow structures develop and expand (Fig. 21e). Thereafter, further opening the guide vanes, the stagnant (reversed flow) region shrinks. Accordingly, the unstable vortical structures gradually integrate and form a large unstable coherent structure that is helically wrapped around the stagnant region ( Fig. 21f and g). An integrated rotating vortex rope is clearly distinguishable at time t ¼ 6.0 s (Fig. 21g) due to KelvineHelmholtz instability of the sharp shear layer. This is in accordance with the results presented in Sections 4.1e4.3, where distinct low-frequency high-amplitude oscillations were observed between t ¼ 4.5 s and t ¼ 6.5 s.
The additional augmentation of the flow rate decreases the runner residual swirl and squeezes the stagnant region. Inevitably, the integrated central vortex becomes more stable and moves toward the center of the draft tube ( Fig. 21h and i). One can see small vortices that rotate around the central axis and merge into a stable slender vortex that is attached to the runner cone ( Fig. 21j and k). Finally, the turbine reaches the BEP condition where a small negative (counter-rotating) swirl leaves the runner and forms a stable and nearly stationary vortex at the center of the draft tube.

Conclusion
The present paper provides a detailed numerical study on the pulsations originated by the transient flow features during the startup of a high-head model Francis turbine. Our results contribute to a better knowledge of the evolution of the flow in the hydraulic turbines during startup operation that could provide the possibility to avoid harmful conditions and have a better estimate of maintenance intervals and costs.
The high-frequency pulsations generated by the blade passing rotor-stator interaction (30f n ) were the dominant excited frequency in the vaneless space throughout the entire startup sequence, while the guide vane passing frequency (28f n ) was the dominant mode for the pressure probe inside the runner rotating domain. Lowfrequency high-amplitude oscillations were observed in the middle of the sequence, suggesting the formation of the RVR. A signal decomposition of the draft tube pressure indicated that the complex flow structures formed at minimum and deep part load conditions have strong plunging (synchronous) effects. Increasing turbine load gave a sudden rise in rotating (asynchronous) mode during the formation of the RVR.
Low-frequency oscillations of the RVR affect the radial forces acting on the runner more than the axial force. The blade passing frequency is less isolated in the axial force compared to the radial component. The axial force is greatly affected by the plunging effects at the minimum load condition and thereby its STFT shows mainly stochastic frequencies. A frequency analysis of the fluctuating radial force exerted on the guide vanes revealed a broad span of stochastic frequencies at the minimum load condition due to the massively separated flow field behind the nearly closed guide vanes. The velocity field in the draft tube revealed the presence of a Fig. 21. Illustration of draft tube vortical structures using an iso-surface of l 2 ¼ 750 s À2 at different times corresponding to different guide vane openings (a).
large quasi-stagnant region with a large positive residual swirl that reduces during the sequence. Large persistent vortical structures were observed inside the draft tube at the minimum load condition. They are responsible for the low-frequency oscillations in such conditions. Gradually increasing the turbine load results in an integration of the unstable vortical structures and the formation of the RVR. At BEP a slender stable vortical structure was observed near the center of the draft tube.
CRediT authorship contribution statement

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.