Axial dynamic load identification of hydraulic turbine based on Chebyshev orthogonal polynomial approximation

An analysis method is proposed to identify axial dynamic loads acting on the Francis turbine based on Chebyshev orthogonal polynomial expansion theory. Dynamic loads are expressed as functions of time and polynomial coefficients. The dynamic load identification model is constructed through discretized integral convolution of the loads, such as the Duhamel integral. However, the discretized numerical integral has a time-cumulative error problem that decreases the recognition accuracy of the dynamic load. Compared with the traditional method, the algorithm proposed in this paper constructs the relationship between the modal displacement and force using polynomial orthogonality and derivative relation between displacement and velocity or acceleration. The new method could avoid the Duhamel integral and time-cumulative error problem. This algorithm not only requires less measuring point information, but is also highly efficient. Compared with genetic algorithm identification, orthogonal polynomial algorithm is not easy falling into local convergence, and does not require multiple repetitions positive analysis trial to evaluate individual fitness value. Numerical simulations demonstrate that the identification and assessment of dynamic loads are effective and consistent when the proposed method is used.


Introduction
Nowadays, dynamic loads play an important role in many practical engineering problems [1][2][3][4][5][6].The dynamic load must be known when various methods are used in dynamics analysis to ensure the reliability and safety of engineering structures.However, in some cases, such as axial dynamic loads of hydro-generator turbines, it is difficult to directly measure the dynamic loads due to limits of technical or economic conditions.However, the dynamic response is easily measured, which leads to the development of the theory of load identification.
There are two main identification methods: frequency-domain technique and time-domain technique [7,8].The frequency-domain technique determines the load spectrums by the frequency response functions and the response spectrums measured, or calculates the modal forces in modal space through coordinate changes [9][10][11][12].Compared with frequency-domain identification, the time-domain technique has good accuracy and clear physical meaning.H. Ocry, H. Glaser [13] determined the load firstly using the time-domain technique through modal coordinates changed in 1985.Recently, the time-domain technique has been greatly developed.
In practical engineering problems, the physical, geometrical, and boundary characteristics are generally uncertain due to errors in manufacture and installation, as well as complex working conditions.This prevents traditional identification methods from achieving good accuracy.Therefore, developing an efficient method to estimate the load of the rotating machinery has theoretical value and wide engineering significance.R. Tiwari and V. C. simultaneously identified the residual unbalanced load and bearing dynamic parameters [14].Chu and Zhang proposed an improved time-domain identification method through the Duhamel integral; Zhang also identified a two-dimensional distributed dynamic load through the Duhamel integral, and pointed out that the Duhamel integral has a time-cumulative error effect [15][16][17].Some researchers determined a hydro-electric unit dynamic load through genetic algorithm or heuristic search algorithm according to the measured displacement [18][19][20][21].The two algorithms identification processes require a number of trials to obtain better population reproduction.Furthermore, the calculation efficiency is poor.Liu et al. identified structures for dynamic loads through various methods, such as the Galerkin method which weighed the total least squares method; Gegenbauer polynomial approximation method; second-order Taylor-Series expansion method, SVM regression approach; and virtual work principle [22][23][24][25][26][27].Some scholars are committed to solve the impact of noise during the identification by proposing a method that integrates denoising, correction, regularization preconditioned iteration method, and Kalman filter or Monte Carlo filter method [28][29][30].
For a water turbine generator set, the distances of upper and lower guide bearings are asymmetrical relative to the rotor.The rotor rotates forming arcuate cyclotron, has not only horizontal displacement, but also vertical displacement (Angular displacement).That is the gyroscopic effect.Usually, modal uncoupling of gyroscopic systems is difficult to carry out without additional assumptions and simplifications.However, this study aims to identify the unknown desired axial dynamic load of a rotor bearing system based on Chebyshev orthogonal polynomial theory, instead of radial dynamic loads.This study constructs the unit shaft system and pier coupling vibration model consider vertical freedom only.The gyroscopic effect is not considered during the modal decoupling.
This paper is organized as follows.Section 2 introduces the basic theories of the Chebyshev Polynomials and determines the algorithm of the dynamic load on the basis of orthogonal polynomial expansion.Compared with other research that used the traditional method, this study constructs the relationship between the modal displacement and force using the Chebyshev Polynomial orthogonality instead of the Duhamel integral.This method could avoid a Duhamel integral time-cumulative error effect.Moreover, Section 3 adopts the proposed method to solve the determined equation and identify the unknown desired dynamic load through numerical examples.Finally, Section 4 provides the conclusion.

Dynamic load identification based on orthogonal polynomial decomposition
Generalized Chebyshev orthogonal polynomials can be used to approximate any continuous single-valued function.As long as the order is reasonable and sufficient, it is better to identify the time-varying dynamic load [17].
The fitting form of the generalized time interval [0, ] is as follows: The recursive formula is shown below: The Chebyshev orthogonal polynomial coefficients are determined by: where, ℎ( ) is weight function.
The system motion equation for an arbitrary order multi-degrees is: where [ ], [ ], and [ ] are the mass matrix, damping matrix, and stiffness matrix of the system; ( ), ( ), ( ) are acceleration, velocity and displacement of the measuring point.( ) is the load on node , which is waited to be identified.
Assuming that damping is proportional damping, modal coordinate transformation is introduced: where, [Φ] is principal modal matrix, ( ) is the th modal displacement: Left-multiplied [Φ] , the system can be decoupled into single degree of freedom systems according to modal shape orthogonality: * ( ) + * ( ) where * , * , * and ( ) are the modal mass, modal damping, modal stiffness, and modal force in the th order; and ( ), ( ), ( ) are modal acceleration, velocity, and displacement in the th order.Chebyshev orthogonal polynomial expansion of the th order modal force ( ) is: According to Duhamel integral: where Substitute Eq. ( 10) into Eq.( 6): where, [ ] = [ ⋯ ], = ⋯ .It is assumed that each order modal force needs terms expression to satisfy accuracy requirement.The system has order modals.The number of coefficients to be identified is * .Dispersing the measuring point displacement at different * times in the calculation interval could obtain * equations.All the orthogonal polynomials coefficients could be obtained by solving the equations.Subsequently, the dynamic load could be identified according to Eqs. ( 8) and ( 9).
The method above constructs the relationship between modal displacement and modal force by Duhamel integral in modal space.Eq. (10) shows that integral upper limit of the Duhamel integral is calculated instantaneously.In others words, the Duhamel integral calculation must be performed at each discrete time; the integral results are different at different times.Thus, the method above not only requires a large amount of computation, but also has a time-accumulation error problem.
The problem could be improved by decreasing time step and increasing orthogonal polynomials order number.However, it needs more calculation and require more measuring point information.
In this study, a new algorithm is proposed to solve the problem.The algorithm constructs the relationship between modal displacement and modal force by polynomial orthogonality and derivative relation between the displacement and velocity or acceleration.This method could avoid Duhamel integral and integral time-cumulative error problem as well as increase identification accuracy.
The orthogonal polynomial expansions of the th order modal force, displacement, velocity and acceleration are as follows: Substitute Eq. ( 12) into Eq.( 8): Assuming that the calculation of time interval is [0, ], multiply ( )ℎ( ) on both sides of Eq. ( 13) and integral in interval [0, ].According to polynomials weighted orthogonality: where , and are determined by: According to the derivative relation between the displacement and acceleration or velocity: Substitute Eq. ( 16) into Eq.( 15): Substitute Eq. ( 17) into Eq.( 14).The new relationship is constructed between displacement and force in mode space: where [ ] is the unit matrix, and: , where , are the coefficient vectors of the modal displacement and force in th order.The measuring point time domain displacement could be expressed by principal modal matrix and modal displacement.The orthogonal polynomial coefficients could be solved by constructing equations through dispersing the measuring point displacement at different times.Consequently, dynamic load identification could be performed.
The algorithm proposed in this paper constructs the relationship between modal displacement and force according to polynomial orthogonality and derivative relation between the displacement and velocity or acceleration.This method could avoid Duhamel integral and increase identification accuracy.Although each element in [ ] or [ ] matrices must be calculated during the integration, the integral upper limit is the same.That is, the integral result remains constant at different times.Implementing integral operation at each discrete time becomes unnecessary.The algorithm increases recognition accuracy and improves calculation efficiency.
The modal mass, stiffness, frequency and shape must be known during load identification.The frequency and mode shape could be easily obtained through modal analysis.The modal mass could be obtained through Eq. ( 20) as demonstrated in an earlier study [31]: where is the th order modal kinetic energy.Furthermore, the modal stiffness and damping could be deduced.

Numerical examples
This study used a practical hydro-generator unit as example.A hydro-generator shaft system and pier coupling vertical vibration model was constructed.The pier was regarded as a rigid foundation.The thrust bearing and frame centrosome was regarded as a spring element.The large shaft was simplified as a non-mass elastic beam.Its mass was distributed into three nodes, which represent the exciter, rotor and runner.The rotor arm was simplified as a non-mass elastic rod.The average quality was assigned to the edge of the rotor bracket and frame centrosome.The mass spring system model is shown in Fig. 1.Moreover, the values and meanings of all parameters are shown in Table 1 [5][6].The total stiffness matrix of the system is: Magnetic yoke and magnetic pole mass at rotor edge 1.2×10 5

kg
Under frame centrosome mass+1/2 total under frame arms mass 7.25×10 10 (N/m) tensile stiffness of the shaft segment 5.72×10 10 (N/m) tensile stiffness of the shaft segment 2.32×10 10 (N/m) Total vertical stiffness of rotor frame arms 2.20×10 12 (N/m) thrust bearing equivalent vertical stiffness between m2 and m5 9.41×10 9 (N/m) Total vertical stiffness of under frame arms The lumped mass matrix is: The damping matrix was acquired using the Rayleigh proportional damping model.The damping ratio was 0.05.The first and second order circular frequencies were used to calculate the coefficients of mass matrix and stiffness matrix in Eq. ( 23): where = 2 ( + ) ⁄ , = 2 ( + ) ⁄ .The modal parameters obtained through model analysis are shown in Table 2.
The dynamic load above was applied to the runner node, which is node 3. The under-frame centrosome displacement could be obtained through calculation.The response and additional 10 % random noise was regarded as the known measuring point displacement.The axial hydro-thrust was identified using the orthogonal polynomial decomposition algorithm.
The identification used two cases.Case 1 constructed the relationship between modal displacement and force through the Duhamel integral.Case 2 constructed the relationship according to the polynomial orthogonality and derivative relationship between displacement and velocity or acceleration as mentioned above.The identified dynamic loads of the two cases are shown in Fig. 2 and Fig. 3.The figures show that the algorithm proposed in this paper avoids Duhamel integral and improves time-cumulative error effect.The integral result of the [ ] and [ ] matrices remain constant at different times.Implementing integral operation at each discrete time has become unnecessary.The algorithm increases recognition accuracy and improves calculation efficiency.Fig. 4 shows relative instantaneous error and time-cumulative error of the two algorithms.The cumulative errors are sum of absolute value of the relative error in each calculation time.The cumulative error of the proposed method could converge at a certain time step.However, the cumulative error of the Duhamel integral method would diverge with time step.There are two kinds of factor which lead to the cumulative error.First, the discrete step could not be infinitely small because of the calculation efficiency.For example, in Fig. 5, the cumulative error is difference between the discrete summation and continuous integral.Second, the finite element load applying method could not be an ideal state step mutation.For instance, in Fig. 6, difference of the load applying methods also lead to a cumulative error.An empirical method in literature [32] could be used to determine the order number of the orthogonal polynomial.However, the formula mainly be applied to two-dimensional orthogonal polynomials for distributed dynamic load identification.For one-dimensional orthogonal polynomials, the cumulative error convergence trend chart with orthogonal polynomial order number could be used to determine the order number.Fig. 7 shows the curves of identification cumulative error with the orthogonal polynomials order number.The cumulative error decreases as the order increases.When the order number is up to 16, the error is convergent, and the result is credible.
The identified dynamic load was applied to the structure.The whole structure dynamic response can be obtained.Using the displacement response of node 5 multiplied by the stiffness of , the dynamic load transferred from the under frame and the pier could be obtained, that is, the vertical exciting load of the pier.

Conclusions
An analysis method is proposed to identify dynamic loads based on Chebyshev orthogonal polynomial expansion theory.This algorithm not only requires less measuring point information, but is also highly efficient.Compared with genetic algorithm identification, orthogonal polynomial algorithm is not easy falling into local convergence, and does not require multiple repetitions positive analysis trial to evaluate individual fitness value.
The traditional dynamic load identification model is constructed through discretized integral convolution of the loads, such as the Duhamel integral.The discretized numerical integral has a time-cumulative error problem that decreases the recognition accuracy of the dynamic load.Compared with the traditional method, the algorithm proposed in this paper constructs the relationship between the modal displacement and force using polynomial orthogonality and derivative relation between displacement and velocity or acceleration.The new method could avoid the Duhamel integral and time-cumulative error problem.Performing an integral calculation for a coefficient matrix at each discrete time becomes unnecessary.The proposed algorithm has higher identification accuracy and efficiency than the traditional Duhamel integral algorithm.

Fig. 4 . 5 .
Fig. 4. Identification error of the two cases Fig. 5. Cumulative error due to discretion

Cumulative error of case 1 Cumulative error of case 2 Instantaneous error of case 1 Instantaneous error of case 2 Fig. 6 . 7 .
Fig. 6.Cumulative error due to load handing method Fig. 7. Cumulative error with polynomials order

Table 1 .
Values and meanings of the parameters of the simplified vertical vibration model

Table 2 .
Modal parameters of the vertical vibration model