Dynamic equation and response of turbine runner blade under flow excitation

Abstract

Considering the complex vibration phenomenon of turbine runner blades under the flow excitation, the nonlinear coupling dynamic equation, which contains hydraulic parameters and blade structure parameters, is established by the finite element method, based on geometrical nonlinearity and fluid–structure interaction of runner blades. According to the dynamic equation, the dynamic response of turbine runner blades under the flow excitation is simulated and analyzed through an example and is compared with the experimental results. The research shows that the proposed equation can well reflect the relation between the dynamic performance of turbine runner blades and its hydraulic parameters, structural parameters, and material parameters. Therefore, the equation provides a necessary theoretical basis for further study on deep-seated problems such as the dynamic characteristics of runner blade and the mechanism of vibration fatigue of runner blade.

Introduction

As the primary part of the hydraulic turbine generator unit, the runner blade directly interacts with the water flow with energy, which is the source of the flow excitation of the unit. The vibration problems constantly exist in the runner blade in the actual operation.^1,2 And, the fatigue crack is easily generated on the blade under the long time vibration. While the crack propagates to a certain extent, the shape of the blade will be greatly changed, which causes hydraulic unbalance and a significant decline in the dynamic performance of the unit and even affects the safety of the unit. ³ For example, since the cracks appeared in the runner blades of the hydraulic turbine generator units of Yantan hydropower station, the operational efficiency of the generator units decreased greatly and strong vibration took place in the generator units. As a consequence, this hydropower station was forced to shut down to deal with the cracks, which greatly affected the economic benefits of the hydropower station. ⁴ Three Gorges hydroelectric power station is the largest hydropower station in China, and this world-renowned hydroelectric power station is also under the risk of blade cracking. ⁵ Therefore, in order to ensure the safe, reliable, and efficient operation of hydropower station, it is necessary to study deep-level questions such as the vibration characteristics of turbine blades, the law of dynamic changes of fatigue cracks, and the safety of the operation.

In the beginning, the vibration problems of the runner blade were researched mainly through the experiment.^6,7 But, the cost of experiments was high, the amount of external disturbances factors was big, and some experiments were infeasible indeed. For the above reasons, the numerical method started to be used to study the vibration problems. With the development of numerical computing methods such as the finite elements and the boundary elements, some scholars began to try to analyze the vibration characteristics of blade structure in aqueous environment by the finite element method or the boundary element method. For instance, Dubas and Schuch ¹ analyzed the dynamic characteristics of the runner blades in both air and water by the finite element method. Using the nonlinear elastic transient analysis, Fu et al. ² studied the hydraulic vibration which was caused by the pressure fluctuation of flow on the runner blade. Zhang and colleagues^8,9 took the Francis turbine as the research object and presented a calculation method to solve the strong coupling flow-induced vibration of elastic structure of a small deformation based on the stability finite element method. RA Saeed et al. ¹⁰ studied the property of flow excitation acting on the runner blade with different operating conditions by the flow analysis software CFD (computational fluid dynamics), and so on. However, the operating environment of runner blades is very severe. The runner blade not only undergoes the start-up and the shutdown but also often runs at off-design conditions. And the state of the water flow is usually so unstable for this reason that intense vibration in the runner blade is generated by complex flow excitation. Therefore, in order to do the further study on the dynamic performance of turbine runner blades and reveal the vibration mechanism of runner blades, it is necessary to establish a dynamic model that can fully reflect the intrinsic relation between the dynamic performance of turbine runner blades and hydraulic parameters, structural parameters, and material parameters.

In this article, the nonlinear fluid–structure coupling dynamic equations of Francis turbine runner blades are formulated and the intrinsic relation between the dynamic response of runner blades and hydraulic parameters, structural parameters, and material parameters of blades are revealed by the finite element displacement method. These studies provide a theoretical basis for a further study on the dynamic characteristics of turbine runner blades.

Finite element model of runner blades

As shown in Figure 1, the shape of turbine runner blades is the surface with a finite length. For the simplicity and generality of the model, the surface with a finite length can be assumed as a combination of numerous flat shell elements, so the runner blades are replaced by the rectangular flat shell elements for finite element simulation, during the process of establishing the finite element model of the rectangular flat shell element. Meanwhile, in order to decrease the difficulty of solving equations, the impact of elastic rotation angles is also neglected, as shown in Figure 2.

Figure 1.

Diagram of structure of runner blade.

Figure 2.

Finite element model of rectangular flat shell element.

Displacement components

In the finite element model of rectangular flat shell element shown in Figure 2, the coordinate system o-xyz is the local coordinate system of elements; a is the width of rectangular flat shell element, b is the length of rectangular flat shell element; $u_i,v_i,w_i(i=1,2,3,4)$ are the displacement components for node i along x, y, and z axes, respectively; and the vector of the rectangular flat shell element in the generalized coordinates can be expressed as follows

$$\mathit{u}=[\begin{array}{cccccccccccc}{u}_1\hfill & v_1\hfill & w_1\hfill & u_2\hfill & v_2\hfill & w_2\hfill & u_3\hfill & v_3\hfill & w_3\hfill & u_4\hfill & v_4\hfill & w_4\hfill \end{array}{]}^{\mathrm{T}}$$

According to the finite element method, the mode of neutral plane displacements of the rectangular flat shell element can be written as follows

$$\{\begin{array}{c}u(x,y,t)=\sum _{i=1}^4{N}_{\mathrm{ui}}(x,y)u_i(t)\hfill \\ \multicolumn{1}{c}{w(x,y,t)=\sum _{i=1}^4{N}_{\mathrm{wi}}(x,y)w_i(t)}\\ \multicolumn{1}{c}{v(x,y,t)=\sum _{i=1}^4{N}_{\mathrm{vi}}(x,y)v_i(t)}\end{array}$$

(1)

In above equations, $N_{\mathrm{ui}}(x,y)$ , $N_{\mathrm{vi}}(x,y)$ , $N_{\mathrm{wi}}(x,y)$ are the shape functions, and

$$\{\begin{array}{c}{N}_{u1}=N_{v1}=N_{w1}=\frac{1}{4}(1-\frac{x}{a})(1+\frac{y}{b})\hfill \\ N_{u2}=N_{v2}=N_{w2}=\frac{1}{4}(1-\frac{x}{a})(1-\frac{y}{b})\hfill \\ N_{u3}=N_{v3}=N_{w3}=\frac{1}{4}(1+\frac{x}{a})(1-\frac{y}{b})\hfill \\ N_{u4}=N_{v4}=N_{w4}=\frac{1}{4}(1+\frac{x}{a})(1+\frac{y}{b})\hfill \end{array}$$

(2)

Kinetic energy of elements

According to the finite element method, the kinetic energy of rectangular flat shell element can be represented as follows

$$\begin{array}{l}T=\frac{1}{2}{\displaystyle \underset{s}{\int }\rho (x,y){[{\stackrel{.}{u}}_a(x,y,t)]}^{2}dS}\\ +\frac{1}{2}{\displaystyle \underset{s}{\int }\rho (x,y){[{\stackrel{.}{v}}_a(x,y,t)]}^{2}dS}\\ +\frac{1}{2}{\displaystyle \underset{s}{\int }\rho (x,y){[{\stackrel{.}{w}}_a(x,y,t)]}^{2}dS}\end{array}$$

(3)

where $\rho (x,y)$ is the mass distribution function of elements, ${\stackrel{·}{u}}_a(x,y,t)$ is the absolute velocity of elements along the x-axis direction at the position $(x,y)$ , ${\stackrel{·}{v}}_a(x,y,t)$ is the absolute velocity of elements along the y-axis direction at the position $(x,y)$ , ${\stackrel{·}{w}}_a(x,y,t)$ is the absolute velocity of elements along the z-axis direction at the position $(x,y)$ , and $\mathrm{dS}$ is the area of elements.

Substituting equation (1) into equation (3), and rearranging, the formula is

$$T=\frac{1}{2}{\stackrel{·}{\mathit{u}}}_z^{\mathrm{T}}{\mathit{m}}_t{\stackrel{·}{\mathit{u}}}_a$$

(4)

where ${\stackrel{·}{\mathit{u}}}_a$ is the vector of the absolute velocity of nodes of element, ${\mathit{m}}_t$ is the mass matrix of elements, and

$${\mathit{m}}_t=\underset{s}{\int }\rho (x,y){\mathit{N}}^{\mathrm{T}}\mathit{N}\mathrm{dS}$$

(5)

In the above equation, $\mathit{N}$ is

$$\mathit{N}=\left[\begin{array}{ccc}{N}_{u1}& 0& 0\\ 0& N_{v1}& 0\\ 0& 0& N_{w1}\end{array}\begin{array}{ccc}{N}_{u2}& 0& 0\\ 0& N_{v2}& 0\\ 0& 0& N_{w2}\end{array}\begin{array}{ccc}{N}_{u3}& 0& 0\\ 0& N_{v3}& 0\\ 0& 0& N_{w3}\end{array}\begin{array}{ccc}{N}_{u4}& 0& 0\\ 0& N_{v4}& 0\\ 0& 0& N_{w4}\end{array}\right]$$

Potential energy of elements

In the case of considering the large geometric deformation of runner blades, the strain–displacement relation of the rectangular flat shell element is ¹¹

$$\{\begin{array}{c}{\epsilon }_x=\frac{\partial u}{\partial x}+z\frac{{\partial }^{2}w}{\partial x^2}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\hfill \\ {\epsilon }_y=\frac{\partial v}{\partial y}+z\frac{{\partial }^{2}w}{\partial y^2}+\frac{1}{2}{\left(\frac{\partial w}{\partial y}\right)}^2\hfill \\ {\gamma }_{\mathrm{xy}}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+2z\frac{{\partial }^{2}w}{\partial x\partial y}+\frac{\partial w}{\partial x}\frac{\partial w}{\partial y}\hfill \end{array}$$

(6)

Substituting equations (1) and (2) into equation (6), and rearranging, the following formula can be obtained

$$\mathbf{\epsilon }={\mathit{S}}_0\mathit{u}+\frac{1}{2}\sum _{j=1}^3{\mathit{u}}^{\mathrm{T}}{\mathit{k}}_j\mathit{u}{\mathit{S}}_j$$

(7)

where $\mathbf{\epsilon }$ is

$$\begin{array}{c}\mathbf{\epsilon }={\left[\begin{array}{ccc}{\epsilon }_x\hfill & {\epsilon }_y\hfill & {\gamma }_{\mathrm{xy}}\hfill \end{array}\right]}^{\mathrm{T}}\hfill \\ {\mathit{S}}_0=\frac{1}{4\mathrm{ab}}\left[\begin{array}{ccc}-b+y\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -a+x\hfill & 0\hfill \\ \multicolumn{1}{c}{-a+x}& -b+y\hfill & 0\hfill \end{array}\begin{array}{ccc}b-y\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -a-x\hfill & 0\hfill \\ \multicolumn{1}{c}{-a-x}& b-y\hfill & 0\hfill \end{array}\begin{array}{ccc}b+y\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& a+x\hfill & 0\hfill \\ \multicolumn{1}{c}{a+x}& b+y\hfill & 0\hfill \end{array}\begin{array}{ccc}-b-y\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& a-x\hfill & 0\hfill \\ \multicolumn{1}{c}{a-x}& -b-y\hfill & 0\hfill \end{array}\right],\hfill \\ {\mathit{S}}_1={\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}},{\mathit{S}}_2={\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}},{\mathit{S}}_3={\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]}^{\mathrm{T}}\hfill \\ {\mathit{k}}_1={\mathit{N}}_1^{\mathrm{T}}{\mathit{N}}_1,{\mathit{k}}_2={\mathit{N}}_2^{\mathrm{T}}{\mathit{N}}_2,{\mathit{k}}_3=2{\mathit{N}}_1^{\mathrm{T}}{\mathit{N}}_2\hfill \end{array}$$

and

$$\begin{array}{l}{\mathit{N}}_1=\frac{1}{4ab}\left[\begin{array}{cccccccccccc}0& 0& -b+y& 0& 0& b-y& 0& 0& b+y& 0& 0& -b-y\end{array}\right],\\ {\mathit{N}}_2=\frac{1}{4ab}\left[\begin{array}{cccccccccccc}0& 0& -a+x& 0& 0& -a-x& 0& 0& a+x& 0& 0& a-x\end{array}\right]\end{array}$$

From the stress–strain relation, it is known that

$$\mathbf{\sigma }=\mathit{D}\mathbf{\epsilon }$$

(8)

where $\mathbf{\sigma }$ is $\mathbf{\sigma }={\left[{\sigma }_x{\sigma }_y{\tau }_{\mathrm{xy}}\right]}^{\mathrm{T}}$ and $\mathit{D}$ is the stress–strain transformation matrix, and

$$\mathit{D}=\frac{E}{1-{\mu }^2}\left[\begin{array}{ccc}1\hfill & \mu \hfill & 0\hfill \\ \multicolumn{1}{c}{\mu }& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1-\mu }{2}\hfill \end{array}\right]$$

where $\mu$ is the Poisson ratio of materials, and E is the elasticity modulus.

Then, elastic potential energy of the rectangular flat shell element can be written as follows

$$P_e=\frac{1}{2}\underset{V}{\int }{\mathbf{\sigma }}^{\mathrm{T}}\mathbf{\epsilon }\mathrm{dV}$$

(9)

Substituting equations (7) and (8) into equation (9), respectively, and rearranging, the equation can be written as follows

$$\begin{array}{c}\multicolumn{1}{c}{P_e=\frac{1}{2}{\mathit{u}}^{\mathrm{T}}\mathit{ku}+\frac{1}{2}\sum _{j=1}^3\underset{V}{\int }{\mathit{u}}^{\mathrm{T}}{\mathit{g}}_j\mathit{u}{\mathit{u}}^{\mathrm{T}}{\mathit{k}}_j\mathit{u}\mathrm{d}V}\\ \multicolumn{1}{c}{+\frac{1}{8}\sum _{j=1}^3\sum _{k=1}^3\underset{V}{\int }{\mathit{h}}_{\mathrm{kj}}({\mathit{u}}^{\mathrm{T}}{\mathit{k}}_k\mathit{u}{\mathit{u}}^{\mathrm{T}}{\mathit{k}}_j\mathit{u})\mathrm{d}V}\end{array}$$

(10)

where $\mathit{k}$ is the stiffness matrix of elements, and $\mathit{k}={\int }_V{\mathit{S}}_0^{\mathrm{T}}\mathit{D}{\mathit{S}}_0\mathrm{d}V$ , ${\mathit{g}}_i={\mathit{S}}_0^{\mathrm{T}}\mathit{D}{\mathit{S}}_i$ , ${\mathit{h}}_{\mathrm{ij}}={\mathit{S}}_i^{\mathrm{T}}\mathit{D}{\mathit{S}}_j$ .

Assuming that the fluid between two adjacent blades is steady and the water flow is a small-disturbance flow that is incompressible and non-viscous, the pressure energy of flow of the runner blade on the surface of the rectangular flat shell element can be expressed as follows ¹²

$$P_w={\mathit{u}}^{\mathrm{T}}\mathit{p}$$

(11)

where $\mathit{p}$ is the vector of hydrodynamic pressure of the node of elements.

For non-viscous and incompressible fluid with the perturbation, based on the theory of hydromechanics, the vector of hydrodynamic pressure of the node $\mathit{p}$ can be expressed as follows ¹²

$$\mathit{p}={\mathit{m}}_p\stackrel{··}{\mathit{u}}$$

(12)

where ${\mathit{m}}_p$ is the additional mass matrix and related to hydraulic parameters.

Substituting equation (12) into equation (11), and rearranging, the following formula of fluid pressure energy of the element, $P_w$ , can be obtained

$$P_w={\mathit{u}}^{\mathrm{T}}{\mathit{m}}_p\stackrel{··}{\mathit{u}}$$

(13)

The total potential energy of rectangular flat shell element of runner blade is made up of the elastic potential energy of the element $P_e$ and the fluid pressure energy of the element $P_w$ , and it can be expressed as follows

$$P=P_e+P_w$$

(14)

Substituting equations (10) and (13) into equation (14), and rearranging, the above equation (14) can be written as follows

$$P=\frac{1}{2}{\mathit{u}}^{\mathrm{T}}\mathit{k}\mathit{u}+\frac{1}{2}{\displaystyle \sum _{j=1}^3{\displaystyle \underset{V}{\int }{\mathit{u}}^T{\mathit{g}}_j\mathit{u}{\mathit{u}}^T{\mathit{k}}_j\mathit{u}\mathrm{d}V}}+\frac{1}{8}{\displaystyle \sum _{j=1}^3{\displaystyle \sum _{k=1}^3{\displaystyle \underset{V}{\int }{\mathit{h}}_{kj}({\mathit{u}}^T{\mathit{k}}_k\mathit{u}{\mathit{u}}^T{\mathit{k}}_j\mathit{u})\mathrm{d}V}}}+{\mathit{u}}^T{\mathit{m}}_p\stackrel{··}{\mathit{u}}$$

(15)

Differential equations of motions of elements

Based on the Lagrange equation, differential equations of motion of rectangular flat shell element of runner blades can be obtained

$$\frac{d}{\mathrm{dt}}\left(\frac{\partial T}{\partial \stackrel{·}{\mathit{u}}}\right)-\frac{\partial T}{\partial \mathit{u}}+\frac{\partial P}{\partial \mathit{u}}=\mathit{f}+\mathit{q}$$

(16)

where $\stackrel{·}{\mathit{u}}$ is the vector of the elastic velocity of elements, $\mathit{f}$ is the vector of the external generalized force, and $\mathit{q}$ is the force vector of the element exerted by other elements connecting to this element.

Substituting the formula of kinetic energy equation (4) and the formula of potential energy equation (15) into equation (11) and rearranging, taking no account of the effect of the damping, the differential equations of motions of the rectangular flat shell element of runner blades can be represented as follows

$$\begin{array}{c}({\mathit{m}}_t+{\mathit{m}}_p)\stackrel{··}{\mathit{u}}+\mathit{ku}=\mathit{f}+\mathit{q}-{\mathit{m}}_t{\stackrel{··}{\mathit{u}}}_r\hfill \\ -\frac{1}{2}\sum _{j=1}^3\underset{V}{\int }({\mathit{g}}_j{\mathit{u}}^{\mathrm{T}}{\mathit{k}}_j\mathit{u}+2{\mathit{u}}^{\mathrm{T}}{\mathit{g}}_j{\mathit{k}}_j\mathit{u})\mathrm{dV}\hfill \\ -\frac{1}{4}\sum _{j=1}^3\sum _{k=1}^3\underset{V}{\int }{\mathit{h}}_{\mathrm{jk}}({\mathit{k}}_k\mathit{u}{\mathit{u}}^{\mathrm{T}}{\mathit{k}}_j\mathit{u}+{\mathit{k}}_j\mathit{u}{\mathit{u}}^{\mathrm{T}}{\mathit{k}}_k\mathit{u})\mathrm{dV}\hfill \end{array}$$

(17)

where ${\stackrel{··}{\mathit{u}}}_r$ is the vector of the acceleration of the rigid body of this rectangular flat shell element.

Dynamic equations of the runner blade

If the runner blade is divided into rectangular flat shell elements, then according to equation (17), the dynamic equation of the $i\mathrm{th}$ element of the runner blade in the local coordinate system can be expressed as follows

$$\begin{array}{c}({\mathit{m}}_{\mathrm{ti}}+{\mathit{m}}_{\mathrm{pi}}){\stackrel{··}{\mathit{u}}}_i+{\mathit{k}}_i{\mathit{u}}_i={\mathit{f}}_i+{\mathit{q}}_i-{\mathit{m}}_{\mathrm{ti}}{\stackrel{··}{\mathit{u}}}_{\mathrm{ri}}\hfill \\ -\frac{1}{2}\sum _{j=1}^3\underset{V}{\int }({\mathit{g}}_{\mathrm{ij}}{{\mathit{u}}_i}^{\mathrm{T}}{\mathit{k}}_{\mathrm{ij}}{\mathit{u}}_i+2{{\mathit{u}}_i}^{\mathrm{T}}{\mathit{g}}_{\mathrm{ij}}{\mathit{k}}_{\mathrm{ij}}{\mathit{u}}_i)\mathrm{dV}\hfill \\ -\frac{1}{4}\sum _{j=1}^3\sum _{k=1}^3\underset{V}{\int }{\mathit{h}}_{\mathrm{ijk}}({\mathit{k}}_{\mathrm{ik}}{\mathit{u}}_i{{\mathit{u}}_i}^{\mathrm{T}}{\mathit{k}}_{\mathrm{ij}}{\mathit{u}}_i+{\mathit{k}}_{\mathrm{ij}}{\mathit{u}}_i{{\mathit{u}}_i}^{\mathrm{T}}{\mathit{k}}_{\mathrm{ik}}{\mathit{u}}_i)\mathrm{dV}\hfill \end{array}$$

(18)

where ${\mathit{u}}_i$ is the generalized coordinate vector of the $i\mathrm{th}$ element; ${\stackrel{··}{\mathit{u}}}_{\mathrm{ri}}$ is the vector of the stiffness acceleration of the $i\mathrm{th}$ element; ${\mathit{m}}_{\mathrm{ti}}$ and ${\mathit{m}}_{\mathrm{pi}}$ are the mass matrix and the additional mass matrix of the $i\mathrm{th}$ element, respectively; ${\mathit{k}}_i$ is the stiffness matrix of the $i\mathrm{th}$ element; ${\mathit{f}}_{\mathrm{wi}}$ is the vector of external generalized force exerted on the $i\mathrm{th}$ element; ${\mathit{q}}_i$ is the vector of the applied force; ${\mathit{k}}_{\mathrm{ij}}$ , ${\mathit{g}}_{\mathrm{ij}}$ , and ${\mathit{h}}_{\mathrm{ijk}}$ are the matrix of the $i\mathrm{th}$ element corresponding to ${\mathit{k}}_j$ , ${\mathit{g}}_j$ , and ${\mathit{h}}_{\mathrm{jk}}$ of equation (10), respectively.

The setup is that ${\mathit{R}}_i$ is the transformation matrix between the $i\mathrm{th}$ element coordinates and the global coordinates and ${\mathit{B}}_i$ is the coordinate compatible matrix between the local number and the global number of the $i\mathrm{th}$ element. Hence, the differential equations of motion of the $i\mathrm{th}$ element which is numbered in the global coordinate can be written as follows

$$\begin{array}{c}({\mathit{M}}_{\mathrm{ti}}^e+{\mathit{M}}_{\mathrm{pi}}^e)\stackrel{··}{\mathit{U}}+{\mathit{K}}_i^e\mathit{U}={\mathit{F}}_i^e+{\mathit{Q}}_i^e-{\mathit{M}}_{\mathrm{ti}}^e{\stackrel{··}{\mathit{U}}}_r\hfill \\ -\frac{1}{2}\sum _{j=1}^3\underset{V}{\int }({\mathit{G}}_{1\mathrm{ij}}{\mathit{U}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{ij}}\mathit{U}+2{\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{U}}^{\mathrm{T}}{\mathit{G}}_{2\mathrm{ij}}\mathit{U})\mathrm{dV}\hfill \\ -\frac{1}{4}\sum _{j=1}^3\sum _{k=1}^3\underset{V}{\int }({\mathit{H}}_{1\mathrm{ikj}}\mathit{U}{\mathit{U}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{ij}}\mathit{U}+{\mathit{H}}_{2\mathrm{ijk}}\mathit{U}{\mathit{U}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{ik}}\mathit{U})\mathrm{dV}\hfill \end{array}$$

(19)

where $\mathit{U}$ , $\stackrel{·}{\mathit{U}}$ , $\stackrel{··}{\mathit{U}}$ , ${\stackrel{··}{\mathit{U}}}_r$ are the generalized coordinate vector, the generalized velocity vector, the generalized acceleration vector, and the rigid-body acceleration vector of the runner blade in the global coordinate system, respectively, and

$$\begin{array}{c}{\mathit{M}}_{\mathrm{ti}}^e={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{m}}_{\mathrm{ti}}{\mathit{R}}_i{\mathit{B}}_i\hfill \\ {\mathit{M}}_{\mathrm{pi}}^e={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{m}}_{\mathrm{pi}}{\mathit{R}}_i{\mathit{B}}_i\hfill \\ {\mathit{K}}_i^e={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{k}}_i{\mathit{R}}_i{\mathit{B}}_i\hfill \\ {\mathit{F}}_i={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{f}}_i\hfill \\ {\mathit{Q}}_i={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{q}}_i\hfill \\ {\mathit{K}}_{\mathrm{ij}}={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{k}}_{\mathrm{ij}}{\mathit{R}}_i{\mathit{B}}_i\hfill \\ {\mathit{G}}_{2\mathrm{ij}}={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{g}}_{\mathrm{ij}}{\mathit{k}}_{\mathrm{ij}}{\mathit{R}}_i{\mathit{B}}_i\hfill \\ {\mathit{G}}_{1\mathrm{ij}}={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{g}}_{\mathrm{ij}}\hfill \\ {\mathit{H}}_{1\mathrm{ijk}}={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{h}}_{\mathrm{ijk}}{\mathit{k}}_{\mathrm{ik}}{\mathit{R}}_i{\mathit{B}}_i\hfill \\ {\mathit{H}}_{2\mathrm{ijk}}={\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{h}}_{\mathrm{ijk}}{\mathit{k}}_{\mathrm{ij}}{\mathit{R}}_i{\mathit{B}}_i\hfill \end{array}$$

As a result, the relation between the generalized coordinate vector ${\mathit{u}}_i$ of the $i\mathrm{th}$ element of the runner blade and its generalized coordinate vector $\mathit{U}$ in the global coordinate system is

$${\mathit{u}}_i={\mathit{R}}_i{\mathit{B}}_i\mathit{U}$$

(20)

All differential equations of the motion of elements are added together, and the effect of damping is estimated approximately by the theory of viscous damping. Consequently, a damping term is introduced into the equation. Then, the dynamic equation of the runner blade with the influence of damping is obtained as follows ¹³

$$\begin{array}{c}\mathit{M}\stackrel{··}{\mathit{U}}+\mathit{C}\stackrel{·}{\mathit{U}}+\mathit{KU}=\mathit{F}-{\mathit{M}}_t{\stackrel{··}{\mathit{U}}}_r\hfill \\ -\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^3\underset{V}{\int }({\mathit{G}}_{1\mathrm{ij}}{\mathit{U}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{ij}}\mathit{U}+2{\mathit{B}}_i^{\mathrm{T}}{\mathit{R}}_i^{\mathrm{T}}{\mathit{U}}^{\mathrm{T}}{\mathit{G}}_{2\mathrm{ij}}\mathit{U})\mathrm{dV}\hfill \\ -\frac{1}{4}\sum _{i=1}^n\sum _{j=1}^3\sum _{k=1}^3\underset{V}{\int }({\mathit{H}}_{1\mathrm{ijk}}\mathit{U}{\mathit{U}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{ij}}\mathit{U}+{\mathit{H}}_{2\mathrm{ijk}}\mathit{U}{\mathit{U}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{ik}}\mathit{U})\mathrm{dV}\hfill \end{array}$$

(21)

where $\mathit{C}$ is the system damping matrix, and $\stackrel{·}{\mathit{U}}$ is the generalized velocity vector of runner blades in the global coordinate system, and

$$\begin{array}{l}{\mathit{M}}_t={\displaystyle \sum _{i=1}^n{\mathit{M}}_{ti}^e},{\mathit{M}}_p={\displaystyle \sum _{i=1}^n{\mathit{M}}_{pi}^e},\\ \mathit{M}={\mathit{M}}_t+{\mathit{M}}_p,\mathit{K}={\displaystyle \sum _{i=1}^n{\mathit{K}}_i^e},\mathit{F}={\displaystyle \sum _{i=1}^n{\mathit{F}}_i^e}\end{array}$$

It is worth noting that the interaction forces ${\mathit{Q}}_i^e$ between the elements are offset by each other after adding the differential equations of the motion of elements together. Since the lower band of the runner is located at the lower end of the runner blade and connected with the blade, when the dynamic equation of the turbine runner blade is established using the parameters above, the mass matrix ${\mathit{M}}_t$ of the lower band is superposed to the mass matrix as the lumped mass.

According to the derivation process of equation (21), the mass matrix $\mathit{M}$ consists of ${\mathit{M}}_t$ and ${\mathit{M}}_p$ . ${\mathit{M}}_t$ is related to the structural parameters and material parameters of the blade, and ${\mathit{M}}_p$ is introduced in order to consider the influence of the pressure energy of the flow. Therefore, the mass matrix $\mathit{M}$ is related to not only the structural parameters of the system but also the hydraulic parameters, which is the fluid–structure coupling term. $-{\mathit{M}}_t{\stackrel{··}{\mathit{U}}}_r$ is the self-excited inertia force of the system, which is not only related to motion parameters of the system but also related to the structural parameters and hydraulic parameters of the system; thus, it is the fluid–structure coupling term. $\mathit{F}$ is the flow excitation exerted on the runner blade. It mainly includes (1) the pressure fluctuation because of the asymmetry of the force of the flow when the water flows into or out of the runner, (2) the pressure fluctuation caused by the low-frequency vortex in the draft tube, and (3) the pressure fluctuation caused by the Karman vortex at the exit of the guide vane or the runner blade. Meanwhile, the nonlinear term which is introduced because of the geometrical nonlinearity of the runner blade is also included in equation (21), which includes not only the square term of the system generalized coordinate vector $\mathit{U}$ but also the cubic term of the system generalized coordinate vector $\mathit{U}$ . Obviously, this equation (21) is a nonlinear fluid–structure coupling dynamic equation. Based on equation (21), the first approximate solution of the dynamic response of the runner blades of the turbine is obtained using the multi-scale method.

Simulation and analysis of dynamic response

Simulation object

Taking a 302MW Francis turbine whose type is HL286-LJ-800 of a hydropower plant in Guangxi as the research object, ¹⁴ the parameters of the runner are as follows: the maximum diameter of the runner is $D_0=8.33\mathrm{m}$ ; the height is $h=5.19\mathrm{m}$ ; 13 blades are uniformly distributed around the runner; the mass of the single blade is $m_t=7.79\times {10}^3\mathrm{kg}$ ; the mass of the upper canopy is $m_s=101.4\times {10}^3\mathrm{kg}$ ; the mass of the lower band is $m_r=68\times {10}^3\mathrm{kg}$ ; the moment of inertia of the upper canopy and the lower band is $J_s=8.52\times {10}^5\mathrm{kg}{\mathrm{m}}^2$ and $J_x=1.15\times {10}^6\mathrm{kg}{\mathrm{m}}^2$ , respectively; the material density of the blade is $\rho =7.85\times {10}^3\mathrm{kg}/{\mathrm{m}}^3$ ; the modulus of elasticity of the blade is $E=210\mathrm{GPa}$ ; and the Poisson ratio of the blade is $\mu =0.3$ .

In the literature, ¹⁴ the dynamic stress on the runner blade has been measured. The layout mode of measuring points and strain gauges is shown in Figure 3, and the lead of the test system is shown in Figure 4. Strain signals are collected by the telemeter, recorded by TEAC310C tape recorder, and handled by the data acquisition and analysis system HP38528. Figure 5 is the experimental time–domain dynamic response characteristics curve of measuring point 2 with the operating condition 120 MW. ¹⁴ According to Figure 5, the maximum dynamic stress, the minimum dynamic stress, and the root mean square values of the dynamic stress of this measuring point are 26.82, −19.05, and 7.51 MPa, respectively.

Figure 3.

Distribution of the measuring point of dynamic stress on the blade.

Figure 4.

Schematic representation of lead of dynamic stress test.

Figure 5.

Experimental time–domain dynamic response characteristics curve of measuring point 2.

Simulation and analysis

In order to compare with the experimental results in the previous work, ¹⁴ the dynamic characteristics of turbine runner blades are analyzed when the operating condition is 120 MW. Because the operating condition 120 MW is the low-load condition of the Francis hydraulic turbine generator unit, the value of pressure fluctuation in the draft tube is large in this condition. Thus, the main consideration is the pressure fluctuation caused by the asymmetry of the force of the flow when the water flows into or out of the runner and the low-frequency vortex in the draft tube. According to the work by Yang, ¹⁵ when the operating condition is 120 MW, the excitation amplitude of the pressure fluctuation is ${\mathit{F}}_1=150\mathrm{kN}$ due to the asymmetry of the force of the flow when the water flows into or out of the runner, and the excitation frequency of ${\mathit{F}}_1$ is ${\Omega }_1=2.5\mathrm{Hz}$ ; the amplitude of the pressure fluctuation is ${\mathit{F}}_2=1100\mathrm{kN}$ due to the low-frequency vortex in the draft tube, and the excitation frequency of ${\mathit{F}}_2$ is ${\Omega }_2=0.3\mathrm{Hz}$ .

For the simplicity and the generality, the turbine runner blade is divided into six rectangular elements, 12 nodes, and 27 displacements, so there are 27 degrees of freedom for each turbine runner blade. The element number is presented by ①, ②, ③, ④, ⑤, ⑥, and the node number is presented by 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, where $U_1$ , $U_4$ , $U_7$ , $U_{10}$ , $U_{13}$ , $U_{16}$ , $U_{19}$ , $U_{22}$ , and $U_{25}$ are the elastic displacements along the X-axis, respectively; $U_2$ , $U_5$ , $U_8$ , $U_{11}$ , $U_{14}$ , $U_{17}$ , $U_{20}$ , $U_{23}$ , and $U_{26}$ are the elastic displacements along the Y-axis, respectively; $U_3$ , $U_6$ , $U_9$ , $U_{12}$ , $U_{15}$ , $U_{18}$ , $U_{21}$ , $U_{24}$ , and $U_{27}$ are the elastic displacements along the Z-axis, respectively; and O-XYZ is the global coordinate system, as shown in Figure 6. The generalized coordinate vector of the turbine runner blade can be expressed as follows

$$\mathit{U}={\begin{array}{ccccccc}{U}_1\hfill & U_2\hfill & \cdots \hfill & U_j\hfill & U_{26}\hfill & \cdots \hfill & U_{27}\hfill \end{array}}^{\mathrm{T}},(j=1,2,3,\dots ,27)$$

Figure 6.

Schematic representation of the dynamical model of the turbine runner.

According to the dynamic model of the runner shown in equation (22), the intrinsic properties of turbine runner blades are solved through programming by MATLAB. The first five-order natural frequencies of turbine runner blades are shown in Table 1.

Table 1.

First five-order natural frequencies of runner.

Order	1	2	3	4	5
Frequency (Hz)	36.25	63.07	80.97	106.98	158.26

The first approximate solution of the dynamic response of the runner blades of the turbine is obtained using the multi-scale method, based on the dynamic equation of the runner blade shown in equation (21). The dynamic stress along the x and y directions in the local coordinate system in the X and Y directions along the local coordinate system at any location on the blade can be calculated according to equations (7), (8), and (20).

The dynamic stress perpendicular to the outlet edge is mainly simulated and calculated because the cracks on the outlet edge of runner blades are generally perpendicular to the outlet edge. ¹⁴ Figures 7–9 show the simulation time–domain dynamic response characteristics curves of nodes 6, 9, and 12, respectively. Node 6, 9, and 12 are all the nodes on the outlet edge of the blade, and node 6 is the same as the measuring point shown in Figure 3. From to the simulation time–domain dynamic response characteristics curves of nodes 6, 9, and 12, it can be seen that the dynamic stress of node 6 near the upper canopy and node 12 near the lower band is larger than the dynamic stress of the node 6 at the middle of the blade’s outlet edge. According to the theory of fracture mechanics, the node 6 near the upper canopy and the node 12 near the lower band are more prone to generating fatigue cracks than the node 9 at the middle position. In the work by Yang et al., ¹³ the locations of 24 cracks that produced on the runner blades of the turbine runner blades were counted. The statistical results were that 12 cracks appeared near the canopy of the blade, 9 cracks appeared near the lower band of the blade, and the other 3 cracks were on other parts. In other words, the fatigue cracks of the runner blades mainly appear near the upper canopy and the lower band. It is obvious that results of the simulation analysis are basically consistent with the experimental results.

Figure 7.

Simulation time–domain dynamic response characteristics curve of node 6.

Figure 8.

Simulation time–domain dynamic response characteristics curve of node 9.

Figure 9.

Simulation time–domain dynamic response characteristics curve of node 12.

In the meantime, according to the simulation time–domain dynamic response characteristics curve of node 9 shown in Figure 8, the maximum dynamic stress, the minimum dynamic stress, and the root mean square values of the dynamic stress of node 9 are 22.88, −16.89, and 11.02 MPa, respectively, which are basically consistent with the experimental results shown in Figure 5.

Conclusion

1. Considering the flow excitation, the dynamic equations of turbine blades are established by the finite element displacement method. The proposed equations not only contain the fluid–solid coupling term but also contain the nonlinear term due to the large geometric deformation of the blades. The dynamic equations are the nonlinear fluid–structure coupling dynamic equations.

2. When the runner blades are forced by the flow excitation, the dynamic responses near the upper canopy and the lower band of the turbine runner are larger than that at the middle position. As a result, the fatigue cracks are more likely to appear in the upper canopy and the lower band of the runner blade.

3. The established coupling dynamic model of runner blade not only reflects the intrinsic relations between the dynamic responses and hydraulic parameters, structural parameters, and material parameters but also provides the necessary theoretical basis to investigate into the deep-seated problems about the coupling mechanism of vibration and fatigue of runner blade.