NONLINEAR DYNAMIC MODEL OF A TURBINE BLADE CONSIDERING VIBRATION AND CRACK COUPLING

In order to further investigate dynamic characteristics of turbine runner blades under the coupling eﬀect of vibration and crack, in this article, a runner blade with a crack was taken as the research object. The coupling eﬀect of vibration and crack was analyzed, a nonlinear dynamic model considering the coupling of the runner blade was developed, and the vibration and fatigue characteristics were investigated. First, based on the contact characteristics of the breathing crack surfaces in the runner blade under hydraulic excitation, a breathing crack surface contact model was established. Subsequently, a nonlinear dynamic model considering the coupling eﬀect of vibration and crack was obtained. The crack surface contact force and crack stiﬀness matrix were established and the vibration and fatigue characteristics were analyzed. Finally, the feasibility of the dynamic model was veriﬁed by a case study, and the dynamic and vibration fatigue characteristics under the coupling eﬀect of vibration and crack were revealed. The research results show that with propagation of the crack, the crack surface contact force increases and the dynamic stress amplitude at the crack tip increases as well. When the sum of the frequency of the hydraulic excitation and the crack surface contact force acting on the runner blade is close to the natural frequency of the runner blade, a combined resonance will occur. When the coupling eﬀect of vibration and crack is taken into consideration, the vibration fatigue crack propagation model is more accurate and can provide a basis for fatigue strength and life prediction of runner blades.


Introduction
Runner blades are the main force bearing and transmitting components of hydraulic turbine generator units.Under the action of long-term hydraulic excitation, runner blades often experience complex vibration phenomena, and long-term vibration can lead to generation of fatigue cracks in different degrees and different ways.The generation of fatigue cracks can further aggravate the vibration of runner blades, causing a significant decline in the dynamic performance of the unit, and even affect its safety.In order to assess the operation efficiency and reliability of hydraulic turbines, it is necessary to develop a dynamic model of runner blade with cracks.
Currently, the available research on dynamic models of runner blades with cracks has mainly considered the effect of fluid-solid coupled as well as the effects of hydraulic parameters, structural parameters, and crack parameters and has developed dynamic models of runner blades through simulation and theoretical analysis.For instance, Yang et al. (2014) used a singular element of the Ansys software to simulate the tip effect of cracks.A two--dimensional finite element analysis model of the crack was developed, and then, a dynamic model of a blade with crack was constructed.In another study, Fernandes et al. (2016) constructed a finite element model of a crack by using mixed-type finite elements, and developed a dynamic model of a cracked blade with different crack angles.However, under hydraulic excitation, the runner blade crack can expand changing the stiffness and resulting in a change of the dynamic response and the crack propagation rate.The dynamic response and crack propagation rate affect each other.This phenomenon is called the coupling effect between vibration and crack propagation.This coupling relationship can affect vibration characteristics and lead to crack arrest and instability propagation of the runner blade.In order to better reveal the vibration and fatigue characteristics of the runner blades, it is necessary to take the coupling relationship between vibration and fatigue cracks propagation into consideration.
At present, the research on the coupling of vibration and crack has been mainly focused on two aspects: crack stiffness and crack surface contact force.In terms of crack stiffness, scholars have considered changes in the stiffness during crack opening and closing, and have developed breathing crack models.For example, Chondros et al. (2001) used a nonlinear model to simulate crack stiffness and investigate vibration characteristics of a cracked beam.Liu and Chen (2010) studied the coupling problem between vibration and fatigue crack propagation of cracked beams using a single degree-of-freedom spring vibrator model and a nonlinear breathing crack stiffness model.Nevertheless, the study ignored the opening-closing phenomenon of the breathing crack, which is local contact behavior.Andreaus and Baragatti (2011) developed a contact finite element model to investigate forced vibration of structures containing cracks.The contact finite element breathing crack model considered elastic deformation when simulating the opening and closing of cracks, increasing the accuracy of the simulated stiffness changes.However, during the opening and closing process of the breathing crack, contact forces are generated on the crack surface, thus, it is necessary to study the contact forces on the crack surface in depth.For instance, Kucher et al. (2007) considered the contact on a blade breathing crack surface as a friction-free contact problem, introduced nonlinear contact stiffness at the crack as a penalty factor, and obtained an approximate expression of the contact force on the crack surface.To solve the convergence problem, Duan and Singh (2007) introduced a tangent smooth function into the Kucher equation, combining the Lagrange multiplier and penalty function methods in a contact problem.Bednarz (2017) considered the effect of crack parameters such as the crack gap on the crack surface contact force of a blade structure, and established an expression for the crack surface contact force using the penalty stiffness method.Although the crack stiffness and surface contact force models constructed in the above research included structural material parameters and crack parameters, the effect of vibration characteristics on the crack stiffness and surface contact force was not taken into consideration.The vibration characteristics affect the crack opening-closing cycle as well as the gap size of the crack during each cycle, which in turn affect the crack stiffness and contact force characteristics on the crack surface.Therefore, it is necessary to take into account the effect of vibration characteristics and develop a mathematical model able to reflect the internal relationship between crack stiffness and crack surface contact force characteristics.
Taking a turbine runner blade with cracks as the research object, this paper constructs a breathing crack contact model of a runner blade, develops a crack stiffness model including structural parameters, material parameters, crack parameters and vibration characteristics, formulates an analytical equation of the crack surface contact force, and obtains a nonlinear dynamic model which couples the vibration and crack propagation characteristics of the runner blade.The dynamic characteristics are obtained by decoupling the dynamic model, and then, the fatigue crack propagation model of the runner blade is obtained.

Breathing crack contact model
The blade can be considered as a shell structure.In order to develop an analytical model for the surface contact of breathing cracks in the blade excited by hydraulic excitation, the following assumptions need to be introduced: (1) It is assumed that the surface contact behavior of the breathing crack is based on a gradual opening and closing process, the change trend of which is linear.(2) When the crack surfaces contact with each other, the effect of axial displacement is relatively small and can be ignored, and only the effect of the radial displacement along the contact surfaces of the crack on the contact force is taken into consideration.(3) According to the opening closing behavior characteristics of breathing cracks, the contact between corresponding point pairs on the crack surfaces can be used to simulate this behavior.
Based on the above assumptions, a breathing crack surface contact model of a runner blade was developed, as shown in

Crack surface contact force
Based on the shell-like mechanical model of the breathing crack surface of the blade shown in Fig. 1.When the crack opens, the surface contact force of the breathing crack in the local coordinate system can be expressed as (Luo et al., 2019) and Equation (3.1) is the crack surface contact force in the first contact form, f ′ a and f ′ b , represent the crack surface contact force acting on points a and b, respectively.Similarly, Eqs.(3.2) is the crack surface contact force in the second contact form, f a ′ and f b ′ represent the crack surface contact force acting on points a ′ and b ′ , respectively.In addition, δ(u a − u b − b c ) is the unit step function, which expressed as Moreover, u a and u b represent the radial displacement of points a and b, respectively, along the x-axis in the local coordinate system, u a ′ and u b ′ represent the radial displacement of points a ′ and b ′ , respectively, along the x-axis in the local coordinate system, b c is the crack gap width, and K * is the contact stiffness of the shell element, also called penalty stiffness, which can be expressed as where V is the volume of the crack contact area of a plane shell element, A c is the crack contact area of a plane shell element, ν is Poisson's ratio of the runner blade material, and E is its elastic modulus.Subsequently, u a and u b in Eqs.(3.1) and u a ′ , u b ′ in Eqs.(3.2), can be calculated based on the neutral surface of the crack, and the distance between the contact-point pair can be transformed into the displacement of the corresponding points A and A ′ on the neutral surface.When the first contact form is expressed in the local coordinate system of the crack in Fig. 1 (II), u a and u b in Eqs.(3.1) can be expressed as follows where u A and u A ′ represent the lateral displacement of points A and A ′ along the x-axis, θ yA and θ yA ′ represent the angular displacement of points A and A ′ around the y-axis, and h denotes the thickness of the shell element.When the second contact form is expressed in the local coordinate system of the crack in Fig. 1 (III), u a ′ and u b ′ in Eqs.(3.2) can be expressed as follows For the first contact form, the contact force acting on the crack surface of points a and b can be obtained by substituting Eq. (3.4) into Eqs.(3.1) For the second contact form, the contact force acting on the crack surface of points a ′ and b ′ can be obtained by substituting Eqs.(3.5) into Eqs.(3.2) The radial and angular displacement of points A and A ′ in the local coordinate system of the crack in Eqs.(3.6) and (3.7) can be expressed by a displacement interpolation function.Therefore, the contact force acting on the crack surface of points a and b in the first type of contact can be expressed as where ψ 1 is 1 × 20 row vector, and (ψ 1 , and the rest are zero, which is related to the structural parameters and crack parameters of the shell element.In addition, l c is the length of the crack from the left end of the plane shell element.
Similarly, the crack surface contact force in the second type of contact can be expressed as where ψ 2 is 1 × 20 row vector, and (ψ 2 , and the other elements are zero, which is related to the structural parameters and crack parameters of the shell element. According to Eqs. (3.8) and (3.9), the contact force on the crack surface is affected by vibration characteristics, which are related to material parameters, such as the elastic modulus and Poisson's ratio, structural parameters such as length and thickness of the shell element, crack parameters such as crack length, gap width, crack angle, contact area and contact area volume, vibration characteristics including radial displacement along x-axis and the angular displacement around the y-axis.
According to Eqs. (3.8) and (3.9), the crack surface contact force of a shell element includes two step functions, (ψ ), which are nonlinear terms, thus, the crack surface contact force is a nonlinear force.

Crack stiffness matrix
According to deformation characteristics of plane shell elements with cracks under the action of crack propagation, strain energy is released.In addition, according to the stress characteristics of runner blades with cracks, considering only the strain energy released by type I (open type) and type II (sliding type) cracks, the strain energy P c released due to crack propagation can be expressed as where K I and K II correspond to mode I and II stress intensity factors at different crack inclination angles, which can be expressed as follows where α is the crack propagation angle, σ α is the normal stress along the crack propagation angle direction at any time, τ α is the shear stress along the crack propagation angle direction at any time.F I (α) and F II (α) are the propagation angle correction coefficients of mode I and II cracks, respectively, which can be determined based on the crack propagation angle α, λ I (l) and λ II (l) represent the stress intensity factor correction functions of mode I and II cracks, respectively, which can be determined based on relative length of the crack.
Considering large geometric deformation of the runner blade and ignoring the effect of bending deformation of the y-axis, the displacement-strain relationship for the crack section of the plane shell element can be expressed as Equation (3.12) can be transformed into a matrix form as follows where ) is 20 × 20 matrix, where and The relationship between stress and strain is where σ c = [σ lx , 0, τ lxy ] T , D is the elastic matrix of the shell element, which can be expressed as By substituting Equations (3.11) 1 to (3.14) into Eq.(3.10), the following expression can be obtained where Subsequently, the shell element stiffness matrix k c caused by crack propagation can be obtained According to Eq. (3.17), the crack-stiffness-matrix is affected by the vibration characteristics, which are related to material parameters, such as the elastic modulus and Poisson's ratio, structural parameters such as length and thickness of the shell element, crack parameters such as crack length, crack propagation angle, crack contact area and contact area volume, and vibration characteristics such as normal and shear stress amplitudes.

Dynamic coupling model of vibration and crack
The relationship between the generalized coordinate vector i of element u i and the generalized coordinate vector U of the shell structure in the global coordinate system can be expressed as where B i is the coordinate matrix between the unit number and system number, and R i is the conversion matrix between the coordinates unit i and the global coordinates.
For large geometric deformation and the crack surface contact force, the dynamic equation of a runner blade with a crack can be written as follows where In addition, U, U and Ü correspond to the generalized coordinate vector, velocity vector, and acceleration vector in the global coordinate system of the blade with crack, Üd is the rigid body acceleration vector, M d Üd is the self-excited inertial force of the system, M t and M p are the mass matrix and additional mass matrix, respectively, C is the damping matrix, K is the stiffness matrix, F and F ′ are the hydraulic excitation and crack surface contact force, respectively.
According to Eq. (4.2), the dynamic equation of the blade with the crack includes the nonlinear crack surface contact force, a nonlinear term representing the large geometric deformation of the blade n i=1 3 dV , and a nonlinear term representing the crack strain energy of the blade therefore the dynamic equation of the blade with the crack is a nonlinear equation.In the dynamic equation of the blade, the F ′ (nonlinear crack surface contact force) and K c U (elastic restoring force) terms caused by the crack as well as the nonlinear terms generated by crack strain energy of the blade include not only crack parameters, but also vibration parameters.Consequently, this is a coupling relationship between vibration and the crack.

Dynamic characteristics of runner blades
By solving Eq. (4.2) using the multi-scale method, the dynamic response of the blade with the crack can be obtained where ϕ is the regular modal matrix and η is the corresponding modal coordinate array of the blade.Subsequently, the approximate solution of the displacement response U of the blade with the crack can be expressed as Based on the developed finite element model of the blade, the relationship between displacement and strain is as follows where ε = [ε x , ε y , γ xy ] T , can be expressed as where S j (j = 1, 2, 3) is 20 × 20 matrix, and The relationship between stress and strain is σ = Dε (5.5) Refer to Eq. (3.15) for the D matrix above.By combining Eqs.(4.1), (5.3) and (5.5), the dynamic stress at any position of the blade with the crack can be obtained as follows

Vibration fatigue characteristics of blades with cracks
According to the Pairs fatigue crack propagation model, the fatigue crack propagation rate of a runner blade with a crack can be expressed as where l is the crack length of the member with the crack, N is the number of cycles of dynamic stress, m and C 0 are material constants of the members with the crack, and ∆K is the stress intensity factor amplitude.Considering the contact and vibration characteristics of the crack surface, the amplitude of the stress intensity factor can be expressed as where F (l) is the crack shape correction factor, which is related to the crack size and which can be obtained from the stress intensity factor manual, ∆K I and ∆K II are the amplitudes of K I and K II , respectively, and ∆f (U) is the amplitude of the contact force acting on the crack surface, which can be determined by combining Equations (3.8), (3.9) and (5.2).By substituting Eq. (5.8) into Eq.(5.7), the crack propagation rate model of the runner blade with the crack can be expressed as According to Equations (5.6) and (5.9), due to the coupling effect between vibration and the crack, the contact force on the crack surface and the crack stiffness change continuously with propagation of the crack, resulting in a continuous change of the dynamic stress characteristics.Consequently, the crack propagation model of the runner blade is related not only to the amplitude of dynamic stress, but also to the amplitude of the contact force on the crack surface.

Case study
In this Section, a large Francis turbine in Guangxi is taken as an example.The main parameters of the runner are as follows: rated speed of 75 RPM, rated head of 59.4 m, rated flow of Q 0 = 580 m 3 /s, rated working condition of 302.In actual operation, due to the thin water outlet edge of the blade and the stress concentration phenomenon at the "T-shaped" connection between the blade and upper crown, penetrating cracks often appear at the water outlet edge of the blade.Table 1 lists the statistical values of through cracks found on the outlet edge of the blade of a hydropower station.The runner blade material is ZG0Cr13Ni5Mo, the material constant is m = 3.65, C 0 = 3.46 • 10 −10 and its allowable stress is σ = 43.27• 10 7 Pa.
In order to verify the correctness of the developed model, No. 3 blade under rated working conditions was taken as the research object.Based on the characteristics of hydraulic excitation acting on the blade, a crack 100 mm away from the upper crown was analyzed.It was assumed that the blade contained only a single crack.Two different crack lengths of 100 mm and 150 mm  and crack gap widths of 1.6 mm and 2 mm were investigated.In order to simplify the simulation model, the runner blade was divided into 9 rectangular elements, each blade had 60 degrees-of--freedom, the element number was represented by 1 ○, 2 ○, A 1 , A 2 ,..., indicated the crack number, the node numbers were represented by 1, 2, . . .(Table 2).The finite element model of the blade with the crack is illustrated in Fig. 2. Time-domain simulation diagrams of the contact force on the crack surface and dynamic stress at the crack tip were obtained according to Eqs. (3.8), (3.9) and (5.6), for a crack length of 100 mm and crack gap width of 1.6 mm.The results are shown in Figs.3a and 3b, respectively.For a crack length of 150 mm and crack gap width of 2 mm, the results are shown in Figs.4a  and 4b, respectively.According to Figs. 3a and 4a, the contact force on the crack surface increases with propagation of the crack.In addition, as the crack gap width increases, the period of the contact force on the crack surface increases as well.According to Figs. 3b and 4b, the dynamic stress amplitude at the crack tip increases with propagation of the crack.
As shown in Figs.5a and 5b, the correctness of the developed model is verified.The fracture module in Ansys workbench was used to establish the crack, the surface unit surface body was inserted at the crack location and the grid was divided, dynamic load was applied and transient dynamic analysis was carried out.The results show that the closer the crack position is, the greater is the dynamic stress.Among them, the dynamic stress at the crack is 32.7 MPa and the mean dynamic stress in Fig. 3b is 30.1 MPa, the error is 2.6 MPa, and the error rate is 7.9%.
According to Eq. (3.17), the crack stiffness increases due to propagation of the crack, which leads to a continuous reduction of the natural frequency of the runner blade with the crack.Moreover, width of the crack gap increases due to propagation of the crack, and the period of the contact force on the crack surface increases as well.When the sum of the frequency of hydraulic excitation and the crack surface contact force acting on the blade is close to the natural frequency of the blade, a combined resonance of the runner blade will be caused.Based on the characteristics of hydraulic excitation acting on the blade, when the crack extends to 132 mm, the combined resonance will occur.By comparing Figs.3b and 4b, it can be observed that when the crack length is 150 mm, the dynamic stress at the crack tip was more irregular than that when the crack length was 100 mm, which also provides a theoretical basis for unstable propagation of the crack under resonance.In Fig. 6, when only dynamic stress is considered without considering the coupling effect of vibration and crack, the crack propagation life curve is the blue line; and when the contact force and dynamic stress on the crack surface are considered, taking the coupling effect of vibration and crack into consideration, the crack propagation life curve is the one red.It can be observed that when the working condition was only the rated load, the crack expanded to the length of 150 mm after 4.53 • 10 7 s (524.3 days).When the coupling effect of vibration and crack was taken into consideration, the crack propagation to 150 mm took 4.42 • 10 7 s (511.6 days).According to (Yang et al., 2014)], the turbine blades are mainly operated near the rated working conditions and are stopped for inspection every 4.25 • 10 7 s (491.7 days).When considering the coupling effect of vibration and crack, the error between the simulation and test results is 4.04%, while, without considering the coupling effect of vibration and crack, the error is 6.63%.This indicates that the vibration fatigue crack propagation model considering the coupling effect of vibration and crack is more accurate, and provides a basis for fatigue strength and life prediction of blades.

Conclusions
• With propagation of the crack, the contact force on the crack surface increases, its period increases, and the dynamic stress amplitude at the crack tip increases as well.• The propagation of the crack increases the crack stiffness, which leads to a continuous reduction in the natural frequency of the blade.When the sum of the frequency of hydraulic excitation and the crack surface contact force acting on the blade is close to the natural frequency of the blade, amd a combined resonance occurs, which also provides a theoretical basis for unstable propagation of cracks under the resonance.• When the coupling effect of vibration and crack is taken into consideration, the developed vibration fatigue crack propagation model is more accurate, and can provide a basis for prediction of the fatigue strength and life of the blades.

Fig. 1 .
According to the third assumption of the breathing crack contact model, two groups of corresponding crack contact point-pairs were: point a to point b (upper surface), and point a ′ to point b ′ (lower surface).During the actual movement of the runner blade, two forms of contact behavior on the crack surface of the runner blade can occur.The first form is displayed in Fig 1 (II) and the second in Fig. 1 (III).

Fig. 1 .
Fig. 1.The breathing crack surfaces contact model of the turbine runner blade 5 MW, maximum runner diameter of 8.33 m, height of 5.19 m, upper crown and lower ring mass of 101.4 • 10 3 kg and 68 • 10 3 kg, respectively, and upper crown and lower ring moment of inertia of 8.52 • 10 5 kg•m 2 and 1.15 • 10 6 kg•m 2 , respectively.In total, 13 blades are evenly distributed in the runner.A single blade has mass of 7.7 • 10 3 kg, width of 1.65 m, thickness of 0.62 m, elastic modulus of E = 210 GPa, material density of 7.85 • 10 3 kg/m 3 , and Poisson's ratio of ν = 0.3.

Fig. 2 .
Fig. 2. Finite element model of the runner blade with cracks

Fig. 3 .
Fig. 3. (a) Simulation curve of the crack surfaces contact force and (b) time domain simulation curve crack-tip dynamic stress when the crack length was 100 mm and the crack gap width was 1.6 mm

Table 1 .
Statistical table of through cracks on the water outlet edge of blades of a hydropower station