Nonlinear Hydraulic Vibration Modeling and Dynamic Analysis of Hydro-Turbine Generator Unit with Multiple Faults

Abstract

Hydraulic vibration is a critical factor affecting the stability and health of the hydro-turbine generator unit (HGU). This paper investigates the nonlinear dynamic characteristics of the HGU occurring from the hydraulic vibration. First, using the turbulence flow model, the lift and drag coefficients of the asymmetric airfoil blades are calculated at different Reynolds numbers and angles of attack. Based on the lift and drag coefficients of the asymmetric airfoil blade, a novel nonlinear mathematical model of asymmetric hydraulic forces is proposed for the asymmetric airfoil blade. Then, combining mechanic and electric asymmetric excitations, a coupled nonlinear mathematical model of the HGU is built using Lagrange equations. Moreover, the correctness of this model is verified. Finally, dynamic behaviors of the HGU changing with attack angle, deviation of the outlet blade angle, outlet guide vane angle, and Reynolds number are analyzed by bifurcation diagrams. It is found that increasing the Reynolds number can enhance the hydraulic vibration, which can lead to more complex nonlinear bifurcation and chaotic motions. In addition, some interesting physical phenomena of HGU under critical ranges and values are obtained. More importantly, these results could give theoretical guidance in the designing and stability research of hydropower units.

1. Introduction

Vibration is an inevitable phenomenon of the hydro-turbine generator unit (HGU) and is one of the important factors affecting the safe production of hydropower stations [1,2,3]. An in-depth analysis of the vibration mechanism of the HGU can better ensure the reliability of the unit’s operation, and effectively avoid or reduce the possible harm of vibration failure to the unit, which is of great significance to its safe and stable operation. In recent years, with the rapid development and wide application of hydropower [4,5,6,7,8,9,10,11], vibration problems of the HGU have received scholars’ great attention.

Due to the complex structure of HGUs, there are many reasons to induce vibration [12,13,14,15,16,17,18,19,20,21,22], and it is caused by the comprehensive action of many factors. Here, the vibration stability of the HGU is analyzed from hydraulic, mechanical, and electrical aspects. Mechanical vibration is caused by the interference of friction and other forces of the unit machinery. The main factors include misalignment of the rotor, unbalance of the rotor mass, collision, and so on. Electrical vibration is the vibration caused by the unbalanced magnetic pull force during the operation of the generator. The main factors are the short circuit of the rotor winding and the uneven air gap. Hydraulic vibration in the HGU is mainly caused by hydraulic pressure fluctuations, vortex forces, hydraulic unbalance, and other factors. Moreover, the main factor generating these is the turbulent flow regime, which was detailedly described by Moreno [23]. In order to study the influence of these vibration problems on the stability of the HGU, nonlinear dynamic theory and methods are introduced. Up to now, a great deal of work has been done on mechanical and electrical vibration, mainly focusing on the misalignment of the rotor, unbalance of the rotor mass, collision, and the uneven air gap. For example, Lees [24] analyzed the characteristics of the misalignment subsystem by establishing the misalignment constraint relationship between two rotors and taking dimensionless factors as parameters. Molka et al. [25] analyzed the response of the faulty rotor system under the condition of misalignment angle and put forward the characteristic frequency representing the misalignment fault in the vibration spectrum. Khanlo et al. [26] studied the performance of correlated chaotic motion caused by the collision of the correlated spindle-disk system. Khanlo. et al. [27] studied the nonlinear dynamic behaviors of the flexible disk system due to a lateral-torsional collision. Dirani et al. [28] established a nonlinear two-dimensional finite element model of the HGU with an uneven air gap. On the basis of the model analysis, the specific effects of an air gap on the magnetic force density of a hydro-turbine generator are obtained. Yan et al. [29] developed a nonlinear fractional-order model of the shafting system for a bending-torsional coupled hydro-generator set considering the coupling reaction of the asymmetric magnetic pull, unbalance of the rotor mass, and contact friction faults and the dynamic behaviors of the system were analyzed. Shi et al. [30] established a mathematical model of the bending-torsional coupling vibration for an unbalanced rotor, and investigated nonlinear vibration sensitivities of the HGU with some parameters by bifurcation diagrams. Xu et al. [31] proposed a new method to calculate the turbine torque. Then, considering the misalignment force, rubbing force, and unbalanced magnetic pull force, a nonlinear mathematical model of the hydro-generator set was built. Moreover, based on the model, the sensitivities of some electro-mechanical parameters to the operation of the unit were obtained. Su and Guo [32] presented a method to investigate the overall nonlinear dynamics of a coupled hydro-generator shaft-foundation system under unbalanced magnetic pull force and nonlinear seal forces, and The variable laws and inner mechanism of the coupled system were revealed further. The results of all these analyses have enriched the dynamical behaviors of a coupled hydro-generator shaft-foundation system. Xu et al. [33] introduced the nonlinear modal method to analyze the dynamic modal interactions between subsystems, and the results given by the different methods were compared to verify the method’s feasibility. The effect of the second-order modes was quantified to investigate its effect on the dynamic characteristics of francis hydro-turbine generator unit (FHTGU), and the vibration characteristics affected by the wind generation system were also investigated. As for hydraulic vibration, many scholars have performed studies on hydraulic variability because of the complex water movement [34,35,36,37,38,39]. For example, Dimitriadis et al. [34] identified the stochastic nature of the water movement and intense hydraulic variability in turbulent regimes. Finite systems of deterministic ordinary nonlinear differential equations were designed to represent forced dissipative hydrodynamic flow, and solutions of these equations can be identified with trajectories in phase space [35]. Frisch reviewed, discussed, and defined turbulence and water movement [36]. Shen et al. [37] analyzed pressure fluctuations in the turbine and identified possible exciting sources that may induce hydraulic vibration in hydropower systems by experimental tests and numerical calculations. Zhuang et al. [38] and Xu et al. [39] proposed nonlinear mathematical modeling of unbalanced hydraulic forces based on the fitting formula of the experimental data of lift and drag coefficients under two-dimensional plate flow and analyzed the effects of hydraulic variability on the stability of the HGU.

Based on the above analysis, mechanical and electrical vibration was the main research target in the early studies on vibration problems of the HGU. There are few nonlinear dynamic studies on hydraulic vibration, and the existing studies on hydraulic vibration also have some shortages. Firstly, due to the difficulty in calculating the lift and drag coefficients of the airfoil blade, the nonlinear mathematical model of hydraulic variability forces was derived by using the fitting formula of the experimental data of lift and drag coefficients under two-dimensional plate flow. Since the two-dimensional plates are symmetric, this method can only be applied to symmetric airfoil blades, but not asymmetric airfoil blades. Secondly, the experiment is troublesome and the working conditions are limited, so the experimental data can not be obtained under any working conditions, which means the obtained mathematical model of the hydraulic variability forces is not systematic and regular. In addition, in engineering practice, airfoil blades are asymmetrical. Therefore, to seek a numerical simulation method to calculate the lift and drag coefficients of the asymmetric airfoil blades is not only of engineering practice value but also of profound academic significance to establish the nonlinear mathematical model of hydraulic variability forces.

In light of the above analysis, we conclude the following innovations. Firstly, a method for calculating the lift and drag coefficients of the hydro-turbine blades at different Reynolds numbers by using the S-A turbulence flow model is proposed. The specific series expressions of the lift coefficient and lift-drag ratio with attack angle at different Reynolds numbers are obtained. Second, considering the forces acting on the airfoil blade, based on the mathematical expressions of the lift coefficient and lift-drag ratio, the nonlinear mathematical model of the asymmetric hydraulic forces for the asymmetric airfoil blade is established for the first time containing attack angle (α), deviation of the outlet blade angle (χ), outlet guide vane angle (α₁), and Reynolds number (Re). Then the nonlinear dynamic model of the HGU including the hydraulic vibration factors is further obtained by bringing it into the Lagrange equation. Third, based on the established dynamic model of the HGU including hydraulic vibration factors, the influences of attack angle (α), deviation of the outlet blade angle (χ), outlet guide vane angle (α₁), and Reynolds number (Re) on the nonlinear dynamic evolution characteristics of the HGU under different Reynolds numbers are analyzed in detail. Besides, the effects of hydraulic instability on the system stability of the HGU are evaluated, which has a certain theoretical guiding significance in the reduction and elimination of hydraulic vibration in practical engineering.

The rest of the paper is organized as follows: In Section 2, the mathematical modeling of asymmetric hydraulic forces for the asymmetric airfoil blade and the establishment of the nonlinear coupled mathematical model of the HGU is presented. Section 3 explains the numerical simulation and analysis. Section 4 presents the conclusions.

2. Mathematical Modeling of HGU

2.1. Modeling of Asymmetric Hydraulic Forces

The manufacturing errors of the outlet flow angle of the hydraulic turbine runner blade have been detected in some engineering observations [15,40], which can cause the inconsistency of the blade channel and lead to strong hydraulic vibration of the HGU. Moreover, when the attack angle of the airfoil blade is large, flow separation and vortex structure will occur that can further enhance the hydraulic vibration. Therefore, it is important and significant to find a mathematical model of asymmetric hydraulic forces for the runner’s airfoil blades to carry out a hydraulic vibration analysis of the HGU.

2.1.1. Radial Hydraulic Force of a Single Airfoil Blade

Force analysis of a single airfoil blade is presented in Figure 1. From Figure 1, when fluid flows over the blade, the pressure difference between the upper and lower surfaces causes lift L. Meanwhile, there is a backward drag D acting on the blade as the water flows. Then, the resultant force R can be estimated. Further, by breaking up the resultant force into radial and peripheral components, the radial force R_m can be obtained.

Based on Figure 1, the angle between the lift (L) and the resultant force (R) is λ, which can be written as

$$\lambda =\arctan \left(\frac{D}{L}\right)=\arctan \left(\frac{C_d}{C_L}\right)$$

(1)

where C_L and C_d are the lift and drag coefficients, separately.

Here, for the convenience of modeling, we define k as the drag-lift ratio, namely that

$$k=\frac{D}{L}=\frac{C_d}{C_L}$$

(2)

According to the Kutta-Joukowski theorem [40] and combining Equations (1) and (2), the radial hydraulic force R_m of a single blade can be estimated as

$$R_m=C_L\frac{\gamma }{g}\frac{{\left|{\mathbf{W}}_{\mathbf{m}}\right|}^2}{2}\frac{F}{\cos \left\{\arctan \left(k\right)\right\}}\cos \left\{{\beta }_m-\arctan \left(k\right)\right\}$$

(3)

where F is the maximum projected area of turbine blades; W_m is the average relative velocity of the blade, and β_m represents its direction angle; γ is the unit volume of liquid weight; g is the gravity coefficient.

From Equation (3), the lift coefficient (C_L), drag coefficient (C_d), and average relative velocity of the blade (W_m) are important parameters influencing the radial hydraulic force, whose expressions and values will be discussed in the following.

2.1.2. Lift and Drag Coefficients of the Asymmetric Airfoil Blade

Lift and drag coefficients (C_L and C_d) are dimensionless of lift L and drag D, and they can be defined as

$$C_L=L/\left(0.5\rho V_{\infty }{}^{2}c\right)$$

(4)

and

$$C_d=D/\left(0.5\rho V_{\infty }{}^{2}c\right)$$

(5)

With the attack angle and Reynolds number changing, lift coefficient C_L and drag-lift ratio k will change as a power series function, which can be expressed as follows:

$$C_L=f\left(\mathbf{Re},\alpha \right)={\displaystyle \sum _{i=0}^5\left\{p_i{\alpha }^i\right\}}={\displaystyle \sum _{i=0}^5\left\{w\left(\mathbf{Re}\right){\alpha }^i\right\}}$$

(6)

and

$$k=g\left(\mathbf{Re},\alpha \right)={\displaystyle \sum _{j=0}^7\left\{q_j{\alpha }^j\right\}}={\displaystyle \sum _{j=0}^7\left\{y\left(\mathbf{Re}\right){\alpha }^j\right\}}$$

(7)

where Re is the Reynolds number; α is the attack angle; p_i and q_j are coefficients of the power series of the lift coefficient C_L and drag-lift ratio k, respectively.

2.1.3. Average Relative Velocity of the Blade

Flow velocity in the blade entrance and exit indicated by subscripts 1 and 2, can be broken into relative velocity (W), convected velocity (U), and absolute velocity (V), which are presented in Figure 2 and Figure 3.

From Figure 2, the average relative velocity of the blade W_m can be expressed as

$$W_{\mathbf{m}}=\left(W_{\mathbf{1}}+W_{\mathbf{2}}\right)/2$$

(8)

where |W_m|, |W₁| and |W₂| represent their values, which can be written as

$$\left|W_{\mathbf{1}}\right|=V_{m1}/\sin {\beta }_1=\mathrm{Q} /\left(F_1\sin {\beta }_1\right)$$

(9)

and

$$\left|W_{\mathbf{2}}\right|=V_{m2}/\sin {\beta }_2=\mathrm{Q} /\left(F_2\sin {\beta }_2\right)$$

(10)

where β₁ represents inlet blade angle; β₂ represents outlet blade angle; Q represents the flow of the hydraulic turbine; F₁, F₂ are the area of the blade entrance and exit section, respectively, which can be expressed as

$$F_1=\left(\pi D_1-Z{S}_u\right)b_0$$

(11)

and

$$F_2=\left(\pi D_2-Z{S}_u\right)D_2$$

(12)

where Z represents the number of runner blades; S_u is the blade thickness; D₁ and D₂ are the diameters of the blade inlet and outlet, respectively; b₀ is the guide vane height.

Based on Equations (6)–(8), |W_m| and β_m can be written as

$$\left|{\mathbf{W}}_{\mathbf{m}}\right|=Q\sqrt{1/{(F_1\sin {\beta }_1)}^2+1/{(F_2\sin {\beta }_2)}^2+2\cos ({\beta }_2-{\beta }_1)/\left(F_1{F}_2\sin {\beta }_1\sin {\beta }_2\right)}/2$$

(13)

and

$${\beta }_m=\arcsin \arcsin \left[2(F_1+F_2)/\sqrt{{\left\{F_2/\left(\sin {\beta }_1\right)\right\}}^2+{\left\{F_1/\left(\sin {\beta }_2\right)\right\}}^2+\left(2F_1{F}_2\cos ({\beta }_2-{\beta }_1)\right)/\left(\sin {\beta }_1\sin {\beta }_2\right)}\right]$$

(14)

2.1.4. Asymmetric Hydraulic Forces of the Turbine Runner

In practical engineering, because of geometrical manufacturing deviations, the outlet blade angles of the runner blades are always inconsistent, which can result in synchronous vibration occurring from an asymmetric hydraulic force rotating with the runner. Here, for the convenience of modeling, we assume that the deviation of the outlet flow angle (χ) only exists in only a pair of blades, χ = β_m2 − β_m1. The asymmetric hydraulic forces acting on the runner are described in Figure 4.

Assumption 1: Since this study provides a method for calculating the lift and drag coefficients of the airfoil blade, as a preliminary generalized study, not specifically for a particular runner, numerical simulations of the flow characteristics based on a two-dimensional airfoil blade are carried out.

Assumption 2: The turbine blade grille density is small, and the interaction between the runner airfoil blades is little. Thus, the liquid flow around the airfoil blade in the blade grille is close to the flow of a single airfoil blade.

From Figure 4, R_m1 and R_m2 are the radial hydraulic forces of the pair of blades with inconsistent outlet flow angles, respectively, and based on Equations (3)–(5), (11) and (12), they can be written as

$$\{\begin{array}{l}{R}_{m1}=\frac{C_L\gamma F}{2g\cos \left\{\arctan \left(k\right)\right\}}{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{1}}\right|}^2\cos \left\{{\beta }_{m1}-\arctan \left(k\right)\right\}\\ R_{m2}=\frac{C_L\gamma F}{2g\cos \left\{\arctan \left(k\right)\right\}}{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{2}}\right|}^2\cos \left\{{\beta }_{m2}-\arctan \left(k\right)\right\}\end{array}$$

(15)

Based on this, radial asymmetric hydraulic forces of the runner considering the asymmetry of the airfoil blade H_m can be expressed as

$$H_m=R_{m1}-R_{m2}=\frac{C_L\gamma F}{2g\cos \left\{\arctan \left(k\right)\right\}}\left\{{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{1}}\right|}^2\cos \left\{{\beta }_{m1}-\arctan \left(k\right)\right\}-{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{2}}\right|}^2\cos \left\{{\beta }_{m2}-\arctan \left(k\right)\right\}\right\}$$

(16)

where |W_m1| and |W_m2| are values of the average relative velocity for the pair of blades with inconsistent outlet flow angle, β_m1 and β_m2 are their direction angles, which are written as

$$\{\begin{array}{l}\left|{\mathbf{W}}_{\mathbf{m}\mathbf{1}}\right|=Q\sqrt{1/{(F_1\sin {\beta }_1)}^2+1/{(F_2\sin {\beta }_{21})}^2+2\cos ({\beta }_{21}-{\beta }_1)/\left(F_1{F}_2\sin {\beta }_1\sin {\beta }_{21}\right)}/2\\ \left|{\mathbf{W}}_{\mathbf{m}\mathbf{2}}\right|=Q\sqrt{1/{(F_1\sin {\beta }_1)}^2+1/{(F_2\sin {\beta }_{22})}^2+2\cos ({\beta }_{22}-{\beta }_1)/\left(F_1{F}_2\sin {\beta }_1\sin {\beta }_{22}\right)}/2\\ {\beta }_{m1}=\arcsin \left[2(F_1+F_2)/\sqrt{{\left\{F_2/\left(\sin {\beta }_1\right)\right\}}^2+{\left\{F_1/\left(\sin {\beta }_{21}\right)\right\}}^2+\left(2F_1{F}_2\cos ({\beta }_{21}-{\beta }_1)\right)/\left(\sin {\beta }_1\sin {\beta }_2{}_1\right)}\right]\\ {\beta }_{m2}=\arcsin \left[2(F_1+F_2)/\sqrt{{\left\{F_2/\left(\sin {\beta }_1\right)\right\}}^2+{\left\{F_1/\left(\sin {\beta }_{22}\right)\right\}}^2+\left(2F_1{F}_2\cos ({\beta }_{21}-{\beta }_1)\right)/\left(\sin {\beta }_1\sin {\beta }_2{}_2\right)}\right]\end{array}$$

(17)

Assuming that the position angle of the pair of blades with inconsistent outlet flow angle is δ, H_m can be decomposed into x and y components as

$$\left\{\begin{array}{l}{H}_x=\left|\cos \delta \right|H_m\\ H_y=\left|\sin \delta \right|H_m\end{array}\right\}$$

(18)

Combining Equations (14) and (16), the asymmetric hydraulic forces considering the asymmetry of the runner airfoil blade in x and y components can be written as

$$\{\begin{array}{c}{H}_x=\left|\cos \delta \right|\frac{C_L\gamma F}{2g\cos \left\{\arctan \left(k\right)\right\}}\left\{{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{1}}\right|}^2\cos \left\{{\beta }_{m1}-\arctan \left(k\right)\right\}-{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{2}}\right|}^2\cos \left\{{\beta }_{m2}-\arctan \left(k\right)\right\}\right\}\\ H_y=\left|\sin \delta \right|\frac{C_L\gamma F}{2g\cos \left\{\arctan \left(k\right)\right\}}\left\{{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{1}}\right|}^2\cos \left\{{\beta }_{m1}-\arctan \left(k\right)\right\}-{\left|{\mathbf{W}}_{\mathbf{m}\mathbf{2}}\right|}^2\cos \left\{{\beta }_{m2}-\arctan \left(k\right)\right\}\right\}\end{array}$$

(19)

where δ = δ₀ + ωt.

2.2. Modeling of the Mechanic and Electric Asymmetric Forces

Except for the asymmetric hydraulic forces of the runner, the dynamic response of the HGU is also significantly affected by some other asymmetric forces. They are summarized in

Table 1. Forces acting on different components of the hydro-turbine generator unit (HGU).

Component.	Damping Force	Oil-Film Force	Asymmetric Magnetic Pull	Hydraulic Asymmetric Force	Rub-Impact Force
Generator rotor	√, Ref. [17]		√, Ref. [41]		√, Ref. [31]
Turbine runner	√, Ref. [17]			√
Generator bearing		√, Ref. [42]
Water guide bearing		√, Ref. [42]

The “√” represents that forces acting on the components of the hydro-turbine generator unit (HGU) can be seen in corresponding references.

, and the specific expressions are shown in the listed references.

2.3. Modeling of the HGU

The shaft system of the HGU is installed in the vertical position. Here, only transverse vibrations are considered and both the rotor and coupling are assumed rigid. The structure of the HGU is presented in Figure 5.

From Figure 5, the coordinate relationships [42] of the HGU can be written as

$$\{\begin{array}{l}{x}_{o2}=x_{o1}+r\cos \theta ,y_{o2}=y_{o1}+r\sin \theta \\ x_{c1}=x_{o1}+e_1\cos \phi ,y_{c1}=y_{o1}+e_1\sin \phi \\ x_{c2}=x_{o2}+e_2\cos \phi ,y_{c2}=y_{o2}+e_2\sin \phi \end{array}$$

(20)

where (x_o1, y_o1) and (x_o2, y_o2) are axis coordinates for the rotor and runner, separately; θ₀ is the initial phase of the HGU, and θ = ωt + θ₀; (x_c1, y_c1) and e₁ are mass coordinates and eccentricity for the rotor, respectively; (x_c2, y_c2) and e₂ represent mass coordinates and eccentricity for the runner, respectively; φ= ωt + φ₀.

The Lagrange function [42] of the HGU is

$$\begin{array}{cc}\hfill L& =T-U=\left(T_r+T_G\right)-\left(\frac{1}{2}{k}_1{\left|o{o}_1\right|}^2+\frac{1}{2}{k}_2{\left|o{o}_2\right|}^2\right)\hfill \\ & =\frac{1}{2}{m}_{o1}\left({\dot{x}}_{o1}{}^2+{\dot{y}}_{o1}{}^2+e_1^2{\dot{\phi }}^2-2{\dot{x}}_{o1}{e}_1\dot{\phi }\sin \phi +2{\dot{y}}_{o1}{e}_1\dot{\phi }\cos \phi \right)+\frac{1}{2}\left(J_1+m_{o1}{e}_1^2\right){\dot{\phi }}^2\hfill \\ & -\frac{1}{2}{k}_{o1}\left(x_{o1}^2+y_{o1}^2\right)-\frac{1}{2}{k}_{o2}\left(x_{o1}^2+y_{o1}^2+2r{x}_{o1}\cos \theta +2r{y}_{o1}\sin \theta +r^2\right)\hfill \\ & +\frac{1}{2}{m}_{o2}\left({\dot{x}}_{o2}{}^2+{\dot{y}}_{o2}{}^2+e_2^2{\dot{\phi }}^2-2{\dot{x}}_{o2}{e}_2\dot{\phi }\sin \phi +2{\dot{y}}_{o2}{e}_2\dot{\phi }\cos \phi \right)+\frac{1}{2}\left[J_1+m_{o2}\left(r^2+e_2^2\right)\right]{\dot{\theta }}^2\hfill \end{array}$$

(21)

where J₁, m₁, k₁ are the rotation inertia, the mass, and the bearing stiffness of the generator rotor, respectively; J₂, m₂, k₂ are the rotation inertia, the mass, and the bearing stiffness for the turbine runner, respectively.

The Lagrange equations [43] of the HGU are

$$\{\begin{array}{l}\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{x}}_{o1}}\right)-\frac{\partial L}{\partial x_{o1}}=\Sigma F_x\\ \frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{y}}_{o1}}\right)-\frac{\partial L}{\partial y_{o1}}=\Sigma F_y\end{array}$$

(22)

From the above, Equations (19) and (20) can be combined into

$$\{\begin{array}{cc}\hfill \frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{x}}_{o1}}\right)-\frac{\partial L}{\partial x_{o1}}& =\left(m_{o1}+m_{o2}\right){\ddot{x}}_{o1}-\left[\left(m_{o1}{e}_1+m_{o2}{e}_2\right)\sin \phi +m_{o2}r\sin \theta \right]\hfill \\ & -\left[\left(m_{o1}{e}_1+m_{o2}{e}_2\right)\cos \phi +m_{o2}r\cos \theta \right]{\omega }^2-\left(k_{o1}+k_{o2}\right)x_{o1}+k_{o2}r\cos \theta \hfill \\ & =F_{x-ump}+F_{x-oil}-F_{x-f}+H_x+F_{x-rub}\hfill \\ \hfill \frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{y}}_{o1}}\right)-\frac{\partial L}{\partial y_{o1}}& =\left(m_{o1}+m_{o2}\right){\ddot{y}}_{o1}+\left[\left(m_{o1}{e}_1+m_{o2}{e}_2\right)\cos \phi +m_{o2}r\cos \theta \right]\hfill \\ & -\left[\left(m_{o1}{e}_1+m_{o2}{e}_2\right)\sin \phi +m_{o2}r\sin \theta \right]{\omega }^2+\left(k_{o1}+k_{o2}\right)y_{o1}+k_{o2}r\sin \theta \hfill \\ & =F_{y-ump}+F_{y-oil}-F_{y-f}+H_y+F_{y-rub}\hfill \end{array}$$

(23)

Furthermore, Equation (21) can be rewritten as

$$\{\begin{array}{l}{\dot{x}}_{o1}=v_x\\ {\dot{v}}_x=1/\left(m_1+m_2\right)\left\{\begin{array}{l}\left[\left(m_{o1}{e}_1+m_{o2}{e}_2\right)\cos \phi +m_{o2}r\cos \theta \right]{\omega }^2+\left(k_{o1}+k_{o2}\right)x_{o1}-k_{o2}r\cos \theta \\ +F_{x-ump}+F_{x-oil}-F_{x-f}+H_x+F_{x-rub}\end{array}\right\}\\ {\dot{y}}_{o1}=v_y\\ {\dot{v}}_y=1/\left(m_1+m_2\right)\left\{\begin{array}{l}\left[\left(m_{o1}{e}_1+m_{o2}{e}_2\right)\sin \phi +m_{o2}r\sin \theta \right]{\omega }^2-\left(k_{o1}+k_{o2}\right)y_{o1}-k_{o2}r\sin \theta \\ +F_{y-ump}+F_{y-oil}-F_{y-f}+H_y+F_{y-rub}\end{array}\right\}\end{array}$$

(24)

3. Numerical Simulation and Analysis

3.1. Calculation and Verification of the Lift and Drag Coefficients

Functions of the lift coefficient C_L and drag-lift ratio k of the NACA2415 asymmetric airfoil blade are established in Section 2.1.2. Here, to get the coefficients of the power series of the lift coefficient C_L and drag-lift ratio k, the calculation and verification are carried out, and the grid for the single airfoil blade is generated by ICEM CFD, which is presented in Figure 6. The grid extends from −9 chords upstream to 14 chords downstream. The upper and lower boundary extends 10 chords from the profile. The grid presented is a “C”-like structured grid. The independent grid can be obtained by increasing grid nodes until it reaches the stage where the solution will not change with the number of nodes. Here, the number of nodes is set as 28,400.

FLUENT is used for the numerical calculation of the lift and drag coefficients (C_L and C_d), and the Spalart-Allmaras model is used to model the turbulence. The finite volume-based method has been used to convert the governing equations of flow into algebraic equations that can be solved numerically. The pressure-velocity coupling has been achieved by the SIMPLE algorithm. The boundary conditions are pressure far-field boundary conditions. The time step of unsteady flow calculation is 0.001 s, and the number of iterations is 2000. To verify the numerical simulation set-up and ensure the accuracy of the results, a comparison of the calculated lift and drag coefficients with the experimental data [44,45] and results from the thin airfoil theory [46] at Reynolds numbers of 2.0 × 10⁵ are presented in Figure 7.

From Figure 7, it can be seen that the lift and drag coefficients obtained by the thin airfoil theory are close to the experimental data [44,45] when the attack angle is smaller than 10 degrees. For the numerical results obtained by the FLUENT S-A turbulence model, the changing curves of lift and drag coefficients are basically consistent with that of the experiment [44,45], and the error is large only around the critical stall attack angle, which is because at high attack angles, there is a separation vortex in the flow field of the airfoil blade, and the accuracy of the turbulent flow model decreases leading to a decrease in the accuracy of the numerical simulation of the lift and drag coefficients of the airfoil blade, which in turn leads to the numerical calculation deviating from the experimental values at high attack angles. Therefore, in general, the study of the unsteady flow around the asymmetric airfoil blade using the FLUNENT Spalart-Allmaras turbulence model in this paper is basically accurate and reliable. Based on this, by numerical simulation, the lift and drag coefficients (C_L and C_d) at different Reynolds numbers are obtained, which are presented in Figure 8. Then, using MATLAB, the fitting of lift coefficient C_L and drag-lift ratio k by power series are carried out and the fitting results are presented in Figure 9, from which orders 5 and 7 have been chosen for the polynomials of lift coefficient C_L and drag-lift ratio k, respectively. The fitting coefficients of the power series of the lift coefficient C_L and drag-lift ratio k are shown in

Table 2. Fitting coefficients of the power series of the lift coefficient C_L and drag-lift ratio k.

Conditions		Coefficients of the Power Series
i or j		7	6	5	4	3	2	1	0
(1)	pi	-	-	0.0003	−0.0027	0.0141	−0.0365	0.1278	0.1541
	qj	−0.0004	0.0035	−0.0174	0.0489	−0.0608	0	0	0.1578
(2)	pi	-	-	0.0001	−0.0009	0.0039	−0.0083	0.1049	0.1824
	qj	−0.0002	0.0018	−0.009	0.0254	−0.0317	0	0	0.0889
(3)	pi	-	-	0.0001	−0.0008	0.0043	−0.0115	0.112	0.1839
	qj	−0.0002	0.0016	−0.0079	0.0225	−0.0281	0	0	0.0793

Note: Conditions are (1): Re = 0.5 × 10⁵; (2): Re = 2.2 × 10⁵; (3): Re = 0.5 × 10⁵.

.

3.2. Nonlinear Dynamic Analysis of the HGU

The nonlinear dynamics of the HGU are solved by the numerical calculation method, and by analyzing the nonlinear dynamic behaviors of the HGU system with parameters, the essential characteristics of the unit system can be analyzed. Based on the fitting coefficients of the power series of the lift coefficient C_L and drag-lift ratio k obtained in Section 3.1, using MATLAB the effects of the attack angle (α), deviation of outlet blade angle (χ), outlet guide vane angle (α₁), and Reynolds number (Re) on the dynamic behaviors of the HGU are investigated. Here, the initial values [42] of the simulation are (x_o1, v_x, y_o1, v_y)^T = (0.0001, 0.0001, 0.0001, 0.0001). On the basis of references [42,47], detailed parameter values conducted in the simulation are: m_o1 = 1.6 × 10⁴ kg, m_o2 = 1.0 × 10⁴ kg, k_o1 = 8.4 × 10⁷ N/m, k_o2 = 6.4 × 10⁷ N/m, e₁ = 0.4 mm, e₂ = 0.4 mm, ω = 14.28 rad/s, Q = 160.3 m³/s, D₁ = 4.34 m, Z = 14, δ₀ = 0.6 rad.

3.2.1. Comparison and Verification

According to the solution grid and calculation conditions in Section 3.1, the velocity streamlines and pressure distribution of the asymmetric airfoil blade under different attack angles (α) at Re = 3.5 × 10⁵ are obtained in Figure 10. Meanwhile, the time-domain diagrams of axis coordinates (x_o1 and y_o1) are presented in Figure 11.

Comparing Figure 10 and Figure 11, when the attack angle (α) is smaller than 10°, there is no separation in the flow field of the airfoil blade. Meanwhile, the vibration amplitude of the HGU is small. When the attack angle (α) reaches 11°, a small clockwise separation vortex appears at the trailing edge of the airfoil, and accordingly, the vibration amplitude shows a slight increase. Then, when the attack angle (α) increases to 13°, the separation vortex and separation region gradually expand and move along the upper airfoil toward the leading edge. At this point, the axis coordinates (x_o1 and y_o1) increase towards a larger vibration amplitude. Especially at α = 19°, the separation vortex increases further and the separation region reaches the maximum thickness of the airfoil. Meanwhile, the axis coordinates (x_o1 and y_o1) enter an oscillation state with a very large vibration amplitude. Based on the above analyses, changes in the separation zone and separation vortex in the flow field of the airfoil blade are basically consistent with the vibration characteristics of the axis coordinates (x_o1 and y_o1). This is because flow separation can lead to a decrease in the lift-drag ratio and flow performance of the airfoil blade, which can further intensify the hydraulic vibration. In fact, the effects of the separation of the flow field and the separation vortex on the performance of the airfoil blade and the HGU vibration have been detected in engineering observations.

Above all, the proposed nonlinear mathematical model of the HGU can reflect the effects of the flow separation and separation vortex on the vibration of the HGU by attack angle (α) well, which also verifies the accuracy of the nonlinear mathematical model of the HGU.

3.2.2. Effects of Attack Angle (α)

Attack angle (α) and the Reynolds number are two critical parameters that can determine the operating conditions of the airfoil bade. Here, to study the specific effects of these two parameters, dynamic evolutions of asymmetric hydraulic forces (H_x and H_y) and axis coordinates (x_o1 and y_o1) of the HGU with increasing attack angle (α) at different Reynolds numbers are shown in Figure 12 and Figure 13, respectively.

As shown in Figure 12 and Figure 13, from the perspective of vibration amplitude, with the attack angle (α) increasing from 0° to 19°, the asymmetric hydraulic forces (H_x and H_y) increase by approximately 2.0 × 10⁵ N, 2.03 × 10⁵ N and 2.08 × 10⁵ N respectively, at Reynolds numbers of 0.5 × 10⁵, 2.2 × 10⁵ and 3.5 × 10⁵. Meanwhile, under the action of asymmetric hydraulic forces (H_x and H_y), the axis coordinates (x_o1 and y_o1) of the HGU have an approximate increase of 0.26 mm, 0.43 mm, and 0.46 mm respectively, at Re = 0.5 × 10⁵, 2.2 × 10⁵ and 3.5 × 10⁵. This is because as the attack angle (α) increases, the flow separation zone and separation vortex of the airfoil blade also increase, which can decrease the flow performance of the airfoil blade and further enhance the hydraulic vibration. On the other hand, from the changing process of dynamic evolutions of the HGU, with the attack angle (α) changing from 0° to 19°, the axis coordinates (x_o1 and y_o1) of the HGU always maintain a two-period oscillation when α < 15.78° at Reynolds numbers of 0.5 × 10⁵. Then, the axis coordinates (x_o1 and y_o1) remain in an unstable state with high variability. With the Reynolds number increasing to 2.2 × 10⁵, the dynamic behaviors of the HGU become more complex. In detail, the axis coordinates (x_o1 and y_o1) show a HOPF bifurcation at α = 8.311° and then enter into an irregular and high vibration state at α = 9.766°. Interestingly, when 12.7° < α < 13.11°, the axis coordinates (x_o1 and y_o1) enter into a short stable state with a two-period oscillation, which changes into an unstable state with high vibration when α > 13.11°. The dynamic evolutions of the axis coordinates (x_o1 and y_o1) at the Reynolds number of 3.5 × 10⁵ are similar to that with a Reynolds number of 2.2 × 10⁵. Detailedly, when the attack angle is 6.957°, a HOPF bifurcation appears that changes into an unstable state with high variability when α = 8.685°. Then the system breaks out of the high variability state and transforms into a periodic oscillation when 11.71° < α < 12.16°. Especially at α = 12.16°, the shaft coordinates transform into an unstable state with high variability again.

From the above, several conclusions can be drawn. First, with the attack angle (α) increasing, changing trends of axis coordinates (x_o1 and y_o1) are consistent with those of asymmetric hydraulic forces (H_x and H_y), which means that the hydraulic instability has decisive effects on the vibration amplitude of the HGU. Second, at Re = 0.5 × 10⁵, the shaft system enters into high variability at α = 17.11°. However, at Re = 2.2 × 10⁵ and Re = 3.5 × 10⁵, a phenomenon of high variability appears at α = 9.766° and α = 8.685°, respectively, which indicates that the dynamic evolutions of the HGU with attack angle (α) are very sensitive to the Reynolds number. In other words, with the increase of the Reynolds number, the bifurcation and high variability points of the axis coordinates (x_o1 and y_o1) move forward. Third, at Re = 2.2 × 10⁵, when α > 8.311°, the hydraulic instability has a large effect on the stability of the HGU. Out of this range, effects are small. Besides, when 12.7° < α < 13.11° at Re = 2.2 × 10⁵ and 11.71° < α < 12.16° at Re = 3.5 × 10⁵, the axis coordinates (x_o1 and y_o1) are in periodic oscillation rather than a high variability state, indicating that within this range, the hydraulic instability has a weak sensitivity to the operation state of the HGU, while it has a strong sensitivity in other ranges. Finally, in actual engineering, the best safety ranges of the attack angle (α) are from 0° to 15° at Re = 0.5 × 10⁵, from 0° to 8° at Re = 2.2 × 10⁵, and from 0° to 7° at Re = 3.5 × 10⁵.

3.2.3. Effects of the Deviation of Outlet Blade Angle (χ)

The opening value of the blade outlet is an important parameter of the blade design, which is a characterization of the overflow. Additionally, when the opening value of the blade outlet for the turbine runner is inconsistent, the outflow of the turbine runner will be not symmetrical, and then it can lead to the generation of asymmetric hydraulic forces (H_x and H_y), which in turn causes hydraulic vibration of the HGU. Particularly, deviation of the outlet blade angle (χ) is an important parameter characterizing the uneven opening value of the blade outlet. For studying the effects of the deviation of the outlet blade angle (χ), dynamic evolutions of asymmetric hydraulic forces (H_x and H_y) and axis coordinates (x_o1 and y_o1) at different Reynolds numbers are described in Figure 14 and Figure 15, respectively.

As shown in Figure 14 and Figure 15, with the increase of Reynolds number and the deviation of the outlet blade angle (χ), vibration amplitudes of asymmetric hydraulic forces (H_x and H_y) and axis coordinates (x_o1 and y_o1) both increase, and their changing trends are similar. In detail, with the deviation of outlet blade angle (χ) increasing from 0 rad to 0.091 rad, the asymmetric hydraulic forces (H_x and H_y) increase by approximately 8.97 × 10⁶ N, 1.01 × 10⁷ N, and 1.037 × 10⁷ N respectively, at Reynolds numbers of 0.5 × 10⁵, 2.2 × 10⁵ and 3.5 × 10⁵. Under the action of asymmetric hydraulic forces (H_x and H_y), the axis coordinates (x_o1) of the HGU have an approximate increase of 1.75 mm, 2.04 mm, and 2.18 mm respectively, at Re = 0.5 × 10⁵, 2.2 × 10⁵, and 3.5 × 10⁵. Meanwhile, the axis coordinate (y_o1) of the HGU increases by 2.2 mm, 2.4 mm, and 2.46 mm respectively, at Re = 0.5 × 10⁵, 2.2 × 10⁵, and 3.5 × 10⁵. These phenomena show that the asymmetric hydraulic forces (H_x and H_y) have determined effects on the stability of the HGU during this changing process. In detail, at Re = 0.5 × 10⁵, 2.2 × 10⁵, and 3.5 × 10⁵, with the deviation of outlet blade angle (χ) increasing from 0 rad to 0.091 rad, the vibration amplitudes of the asymmetric hydraulic forces (H_x and H_y) also increase.

In addition, taking a further study of Figure 15a–f, at Re = 0.5 × 10⁵, 2.2 × 10⁵, and 3.5 × 10⁵, when 0 rad < χ < 0.091 rad, the axis coordinates (x_o1 and y_o1) of the HGU go through two-period oscillation—high variability state—two-period oscillation—multi-period oscillation—HOPF bifurcation—high variability state. However, it is worth noting that the occurrence points of the high variability state are different when the Reynolds numbers are different. To be specific, at Re= 0.5 × 10⁵, the HGU enter a high variability state when 0.0087 rad < χ < 0.0155 rad and χ > 0.07624 rad. Then, at Re = 2.2 × 10⁵, they come into the high variability state when 0.0075 rad < χ < 0.013 rad and χ > 0.0641 rad. Finally, with the Reynolds numbers increasing to 3.5 × 10⁵, the axis coordinates enter into the high variability state when 0.00598 rad < χ <0.01042 rad, and χ > 0.04929 rad.

The above analyses can lead to several meaningful results. First, the derivation of the outlet blade angle (χ) has decisive effects on the dynamic evolutions of the HGU, meaning that the operating state of the HGU is sensitive to the hydraulic instability. Second, when 0.0087 rad < χ< 0.0155 rad at Re = 0.5 × 10⁵, 0.0075 rad < χ < 0.013 rad at Re = 2.2 × 10⁵ and 0.00598 rad < χ < 0.01042 rad at Re = 3.5 × 10⁵, the value of the asymmetric hydraulic forces (H_x and H_y) is small, but the axis coordinate s (x_o1 and y_o1) are in the state of high variability, showing that hydraulic instability could intensify mechanic and electric instabilities among these ranges. In addition, with the increase of Reynolds number, occurrence points of the high variability state move forward, and the ranges of instability increase. Namely, the stability of HGU with the derivation of outlet blade angle (χ) is very sensitive to the Reynolds number. Finally, the derivation of outlet blade angle (χ) should avoid these unstable ranges: 0.0087 rad < χ < 0.0155 rad and χ > 0.07624 rad at Re = 0.5 × 10⁵, 0.0075 rad < χ < 0.013 rad and χ > 0.0641 rad at Re = 2.2 × 10⁵, 0.00598 rad < χ < 0.01042 rad and χ > 0.04929 rad Re = 3.5 × 10⁵.

3.2.4. Effects of Outlet Guide Vane Angle (α₁)

Outlet guide vane angle (α₁) is an important parameter determining the flow condition before the hydraulic turbine runner, which can further affect the asymmetric hydraulic forces (H_x and H_y) and the stability of the HGU.To investigate the influences of the outlet guide vane angle (α₁), dynamic evolutions of asymmetric hydraulic forces (H_x and H_y), and axis coordinates (x_o1 and y_o1) with increasing outlet guide vane angle (α₁) at different Reynolds numbers are presented in Figure 16 and Figure 17.

From Figure 16, at Re = 0.5 × 10⁵, 2.2 × 10⁵, and 3.5 × 10⁵, qualitatively speaking, the dynamic changing process of the asymmetric hydraulic forces (H_x and H_y) are similar. When the outlet guide vane angle (α₁) ranges from 0 rad to 0.11 rad, the value of the asymmetric hydraulic forces (H_x and H_y) decreases rapidly. Especially at α₁ = 0.11 rad, the value of the asymmetric hydraulic forces (H_x and H_y) decreases to about 0. Moreover, as the outlet guide vane angle (α₁) continues to range from 0.11 rad to 0.425 rad, the value of the asymmetric hydraulic forces (H_x and H_y) fluctuates around 0 all the time. Then, with the outlet guide vane angle (α₁) increasing further, asymmetric hydraulic forces (H_x and H_y) begin to increase and show a slight jump when it increases to a certain extent. Finally, the asymmetric hydraulic forces (H_x and H_y) get rid of the jump and go into a periodic oscillation. At the same time, it is worth noting that at Reynolds numbers of 2.2 × 10⁵ and 3.5 × 10⁵, the value of the asymmetric hydraulic forces (H_x and H_y) is basically the same, which is larger than that at Reynolds numbers of 0.5 × 10⁵. Moreover, when 0 rad < α₁ < 0.11 rad and 0.425 rad < α₁ < 1.4 rad, larger differences of the asymmetric hydraulic forces (H_x and H_y) appear at different Reynolds numbers. In addition, quantitatively speaking, with the outlet guide vane angle (α₁) increasing from 0 rad to 1.4 rad, the asymmetric hydraulic forces (H_x and H_y) decrease by approximately 5.5 × 10⁶ N at Reynolds numbers of 0.5 × 10⁵. Under the action of asymmetric hydraulic forces (H_x and H_y), the axis coordinates (x_o1) and (y_o1) of the HGU has an approximate decrease of 8.5 mm and 11 mm at Re = 0.5 × 10⁵. With the Reynolds number increasing to 2.2 × 10⁵ and 3.5 × 10⁵, the asymmetric hydraulic forces (H_x and H_y) decrease by 6.2 × 10⁶ N and 6.4 × 10⁶ N respectively, which have small differences. Therefore, the changes in the vibration amplitudes of the axis coordinates (x_o1) and (y_o1) are almost the same, decreasing by 10.5 mm and 13 mm.

Based on the above, comparing Figure 16 with Figure 17, some important results can be concluded. First, the variation of the dynamic evolutions of the HGU is determined by the asymmetric hydraulic forces (H_x and H_y). Second, the Reynolds number is less sensitive to the outlet guide vane angle (α₁) than to the attack angle (α) and deviation of the outlet blade angle (χ). Dynamic characteristics of the HGU can be affected by large changes in Reynolds number. Third, at Reynolds numbers of 0.5 × 10⁵, 2.2 × 10⁵, and 3.5 × 10⁵, when 0.11 rad< α₁ < 0.425 rad, the vibration amplitude of the system maintains at about 0, meaning that this range is close to the optimal working condition. Out of this range, whether the outlet guide vane angle (α₁) increases or decreases, the vibration amplitude of the HGU will increase significantly. Because a small outlet guide vane angle (α₁) can lead to the water rushing to the front of the guide blade profile, the attack angle is positive, while the outlet guide vane angle (α₁) is large, the water flow will rush to the back of the guide blade profile with a negative attack angle. These two conditions will both aggravate the hydraulic instability and affect the efficiency of the HGU. In fact, in engineering practice, when the HGU is in operation, a larger or smaller guide vane opening will both affect the stability and performance of the unit. Third, when the outlet guide vane angle (α₁) ranges from 0.1008 rad to 0.103 rad and from 0.9925 rad to 1.32 rad, the coordinates of the shaft system enter into high variability at Reynolds numbers of 0.5 × 10⁵, while at Reynolds numbers of 2.2 × 10⁵, they go into high vibration in the ranges from 0.1006 rad to 0.1042 rad and from 0.917 rad to 1.364 rad, which illustrates that larger Reynolds numbers can enhance the hydraulic instability in these ranges. In other words, the axis coordinates (x_o1 and y_o1) have a strong sensitivity to the Reynolds number in these ranges. Finally, in order to guarantee the efficient operation of the HGU, the outlet guide vane angle (α₁) should be controlled within the range of 0.11 rad to 0.425 rad as far as possible.

4. Conclusions

To investigate the nonlinear dynamics of the HGU under hydraulic vibration, the lift and drag coefficients of the asymmetric airfoil blades are calculated at different Reynolds numbers and angles of attack using the turbulence flow model. Based on the lift and drag coefficients of the asymmetric airfoil blade, a nonlinear mathematical model of the asymmetric hydraulic forces (H_x and H_y) is firstly proposed for the asymmetric airfoil blade. Then, taking the mechanic and electric asymmetric excitations into consideration, a nonlinear coupled mathematical model of the HGU is built, and the correctness of this model is verified. Using the latter, the nonlinear vibration characteristics of the HGU are investigated with the changing of attack angle (α), deviation of the outlet blade angle (χ), outlet guide vane angle (α₁), and Reynolds number (Re), and some important physical phenomena are obtained as follows:

(1)

In the changing process of attack angle (α), deviation of the outlet blade angle (χ), and outlet guide vane angle (α₁), asymmetric hydraulic forces have decisive effects on the vibration amplitude of the HGU.

(2)

The Reynolds number is less sensitive to the outlet guide vane angle (α₁) than to the attack angle (α) and deviation of the outlet blade angle (χ), and the dynamic characteristics of the HGU with increasing outlet guide vane angle (α₁) can only be effected by large changes in the Reynolds number.

(3)

Increasing the Reynolds number can enhance the effects of the attack angle (α) and the deviation of outlet blade angle (χ) on dynamic evolutions of the HGU, which leads to the points of the high variability state moving forward. Meanwhile, the stable ranges are respectively analyzed.

(4)

A larger or smaller outlet guide vane angle (α₁) can both enhance the nonlinear hydraulic vibration, and the best range is from 0.11 rad to 0.425 rad. In addition, in the ranges of outlet guide vane angle (α₁) from 0.917 rad to 1.4 rad, the Reynolds number can intensify the nonlinear hydraulic vibration.

The conclusions obtained are a guiding significance to studying the nonlinear dynamics of hydraulic vibration and give theoretical guidance in the designing and stability research of the HGU. Besides, due to the limitations of this study, in future studies, the complex three-dimensional flow characteristics of the hydraulic machinery will be studied in-depth, taking into account the interaction between the turbine blades.