Abstract
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 timecumulative 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 timecumulative 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.
1. Introduction
Nowadays, dynamic loads play an important role in many practical engineering problems [16]. 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 hydrogenerator 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: frequencydomain technique and timedomain technique [7, 8]. The frequencydomain 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 [912]. Compared with frequencydomain identification, the timedomain technique has good accuracy and clear physical meaning. H. Ocry, H. Glaser [13] determined the load firstly using the timedomain technique through modal coordinates changed in 1985. Recently, the timedomain 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 timedomain identification method through the Duhamel integral; Zhang also identified a twodimensional distributed dynamic load through the Duhamel integral, and pointed out that the Duhamel integral has a timecumulative error effect [1517]. Some researchers determined a hydroelectric unit dynamic load through genetic algorithm or heuristic search algorithm according to the measured displacement [1821]. 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; secondorder TaylorSeries expansion method, SVM regression approach; and virtual work principle [2227]. 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 [2830].
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 timecumulative 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.
2. Dynamic load identification based on orthogonal polynomial decomposition
Generalized Chebyshev orthogonal polynomials can be used to approximate any continuous singlevalued function. As long as the order is reasonable and sufficient, it is better to identify the timevarying dynamic load [17].
The fitting form of the generalized time interval [0, $S$] is as follows:
The recursive formula is shown below:
The Chebyshev orthogonal polynomial coefficients ${a}_{j}$ are determined by:
where, $h\left(t\right)$ is weight function.
The system motion equation for an arbitrary $n$ order multidegrees is:
where [$M$], [$C$], and [$K$] are the mass matrix, damping matrix, and stiffness matrix of the system; ${\ddot{u}}_{k}\left(t\right)$, ${\dot{u}}_{k}\left(t\right)$, ${u}_{k}\left(t\right)$ are acceleration, velocity and displacement of the measuring point. ${F}_{m}\left(t\right)$ is the load on node $m$, which is waited to be identified.
Assuming that damping is proportional damping, modal coordinate transformation is introduced:
where, $\left[\mathrm{\Phi}\right]$ is principal modal matrix, $y\left(t\right)$ is the $i$th modal displacement:
Leftmultiplied ${\left[\mathrm{\Phi}\right]}^{T}$, the system can be decoupled into $n$ single degree of freedom systems according to modal shape orthogonality:
where ${M}_{i}^{*}$, ${C}_{i}^{*}$, ${K}_{i}^{*}$ and ${f}_{i}\left(t\right)$ are the modal mass, modal damping, modal stiffness, and modal force in the $i$th order; and ${\ddot{y}}_{i}\left(t\right)$, ${\dot{y}}_{i}\left(t\right)$, ${y}_{i}\left(t\right)$ are modal acceleration, velocity, and displacement in the $i$th order.
Chebyshev orthogonal polynomial expansion of the $i$th order modal force ${f}_{i}\left(t\right)$ is:
According to Duhamel integral:
$=\frac{{\sum}_{j=0}^{\infty}{a}_{ij}}{{M}_{i}^{*}{\omega}_{Di}}{\int}_{0}^{t}{T}_{j}\left(\tau \right){e}^{{\xi}_{i}{\omega}_{i}\left(t\tau \right)}\mathrm{s}\mathrm{i}\mathrm{n}\left[{\omega}_{Di}\left(t\tau \right)\right]d\tau =\left[\begin{array}{llll}{d}_{i0}& {d}_{ij}& \cdots & {d}_{i\infty}\end{array}\right]\left\{\begin{array}{c}{a}_{i0}\\ {a}_{ij}\\ \vdots \\ {a}_{i\infty}\end{array}\right\},$
where, ${C}_{i}^{*}=2{\xi}_{i}{\omega}_{i}{{M}_{i}}^{*}$, ${\omega}_{Di}={\omega}_{i}\sqrt{1{\xi}_{i}^{2}}$, ${d}_{ij}=\frac{1}{{M}_{i}^{*}{\omega}_{Di}}{\int}_{0}^{t}{T}_{j}\left(\tau \right){e}^{{\xi}_{i}{\omega}_{i}\left(t\tau \right)}\mathrm{s}\mathrm{i}\mathrm{n}\left[{\omega}_{Di}\left(t\tau \right)\right]d\tau $.
Substitute Eq. (10) into Eq. (6):
where, $\left[{D}_{i}\right]=\left[\begin{array}{llll}{d}_{i0}& {d}_{ij}& \cdots & {d}_{i\infty}\end{array}\right]$, $\left\{{A}_{i}\right\}={\left\{\begin{array}{llll}{a}_{i0}& {a}_{ij}& \cdots & {a}_{i\infty}\end{array}\right\}}^{T}$.
It is assumed that each order modal force needs $J$ terms expression to satisfy accuracy requirement. The system has $n$ order modals. The number of coefficients to be identified is ${J}^{*}n$. Dispersing the measuring point displacement at different ${J}^{*}n$ times in the calculation interval could obtain ${J}^{*}n$ 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 timeaccumulation 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 timecumulative error problem as well as increase identification accuracy.
The orthogonal polynomial expansions of the $i$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, $T$], multiply ${T}_{g}\left(t\right)h\left(t\right)$ on both sides of Eq. (13) and integral in interval [0, $T$]. According to polynomials weighted orthogonality:
where ${b}_{ig}^{\text{'}\text{'}}$, ${b}_{ig}^{\text{'}}$ and ${b}_{ig}$ 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 [$I$] is the unit matrix, and:
where $\left\{{b}_{i}\right\}$, $\left\{{a}_{i}\right\}$ are the coefficient vectors of the modal displacement and force in $i$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 [$P$] or [$Q$] 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 ${E}_{i}$ is the $i$th order modal kinetic energy. Furthermore, the modal stiffness and damping could be deduced.
3. Numerical examples
This study used a practical hydrogenerator unit as example. A hydrogenerator 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 nonmass elastic beam. Its mass was distributed into three nodes, which represent the exciter, rotor and runner. The rotor arm was simplified as a nonmass 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 [56].
Fig. 1Vertical vibration simplified model of the hydrogenerator
The total stiffness matrix of the system is:
Table 1Values and meanings of the parameters of the simplified vertical vibration model
Equivalent value  Equivalent meaning  
${m}_{1}$  8.28×10^{4 }kg  Exciter mass+1/2 shaft segment mass from top to rotor support frame 
${m}_{2}$  1.042×10^{6} kg  Rotor support frame centrosome mass +1/2 total arms mass +1/2 total shaft mass 
${m}_{3}$  3.29×10^{5} kg  Runner mass+ water addition quality+1/2 shaft segment mass from rotor support frame to runner 
${m}_{4}$  9×10^{5} kg  Magnetic yoke and magnetic pole mass at rotor edge 
${m}_{5}$  1.2×10^{5} kg  Under frame centrosome mass+1/2 total under frame arms mass 
${k}_{1}$  7.25×10^{10 }(N/m)  tensile stiffness of the shaft segment 
${k}_{3}$  5.72×10^{10} (N/m)  tensile stiffness of the shaft segment 
${k}_{4}$  2.32×10^{10} (N/m)  Total vertical stiffness of rotor frame arms 
${k}_{51}$  2.20×10^{12} (N/m)  thrust bearing equivalent vertical stiffness between m_{2} and_{}m_{5} 
${k}_{52}$  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 $\xi $ 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 $\alpha =2\xi {\omega}_{1}{\omega}_{2}/({\omega}_{1}+{\omega}_{2})$, $\beta =2\xi /({\omega}_{1}+{\omega}_{2})$.
The modal parameters obtained through model analysis are shown in Table 2.
Table 2Modal parameters of the vertical vibration model
Modal order  Parameter  
Circular frequency (Hz)  Modal mass (Kg)  Modal stiffness (N/m)  Modal damping (N·s/m)  Modal shape (−)  
1  59.75  2.08×10^{6}  7.42×10^{9}  1.24×10^{7}  [0.87,0.86,0.88,1,0.86]^{T} 
2  204.24  1.59×10^{6}  6.66×10^{10}  3.26×10^{7}  [−0.65,−0.62,−0.81,1,−0.62]^{T} 
3  474.47  4.45×10^{5}  1.00×10^{11}  3.99×10^{7}  [−0.40,−0.29,1,0.04,−0.30]^{T} 
4  972.32  9.01×10^{4}  8.51×10^{10}  3.27×10^{7}  [1,−0.08,0.02,0.002,−0.08]^{T} 
5  4531.06  1.34×10^{5}  2.73×10^{12}  1.04×10^{9}  [0.005,−0.12,0.001,0.0001,1]^{T} 
Construct axial hydrothrust for Francis turbine according to the literature [19]:
$287.4\times 1{0}^{3}\mathrm{c}\mathrm{o}\mathrm{s}\left(2\mathrm{\pi}\times 2.5t+5.07\right)+489.7\times 1{0}^{3}\mathrm{c}\mathrm{o}\mathrm{s}\left(2\mathrm{\pi}\times 8.75t+1.43\right).$
The dynamic load above was applied to the runner node, which is node 3. The underframe centrosome displacement could be obtained through calculation. The response and additional 10 % random noise was regarded as the known measuring point displacement. The axial hydrothrust 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 timecumulative error effect. The integral result of the [$P$] and [$Q$] 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. 2Identification result of Case 1
Fig. 3Identification result of Case 2
Fig. 4 shows relative instantaneous error and timecumulative 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.
Fig. 4Identification error of the two cases
Fig. 5Cumulative error due to discretion
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 twodimensional orthogonal polynomials for distributed dynamic load identification. For onedimensional 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 ${k}_{52}$, the dynamic load transferred from the under frame and the pier could be obtained, that is, the vertical exciting load of the pier.
Fig. 6Cumulative error due to load handing method
Fig. 7Cumulative error with polynomials order
4. 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 timecumulative 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 timecumulative 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.
References

Kashani Mir Tahmaseb, Jayasinghe Supun, Hashemi Seyed M. Dynamic finite element analysis of bearingtorsion coupled beams subjected to combined axial load and end moment. Shock and Vibration, Vol. 2015, 2015, p. 471270.

Hadi Jalali Mohammad, Ghayour Mostafa, ZiaeiRad Saeed, et al. Dynamic analysis of a high speed rotorbearing system. Measurement, Vol. 54, 2014, p. 19.

Si Xiaohui, Lu Wen Xiu, Chu Fuller Lateral vibration of hydroelectric generating set with different supporting condition of thrust pad. Shock and Vibration, Vol. 18, 2011, p. 317331.

Sun Wanquan, Huang Xionghui Identification of vibration transfer path for the coupled system of water turbine generator set and power house. Journal of Vibration and Shock, Vol. 33, Issue 6, 2014, p. 2328.

Ma Zhenyue, Dong Yuxin Dynamics of Water Turbine Generator Set. Dalian University of Technology Press, Dalian, 2003.

Zhi Baoping, Ma Zhenyue, Wu Qianqian Study on transfer paths of vertical vibrations in the head cover system of turbine. Journal of Hydroelectric Engineering, Vol. 32, Issue 3, 2013, p. 241245.

Xu Feng, Chen Huaihai, et al. Force identification for mechanical vibration stateoftheart and prospect. China Mechanical Engineering, Vol. 13, Issue 6, 2002, p. 526531.

Uhl T. The inverse identification problem and its technical application. Archive of Applied Mechanics, Vol. 77, Issue 5, 2007, p. 325337.

Hillary B., Ewins D. J. The use of strain Gauges in force determination and frequency response function measurements. Proceedings of the 2nd International Modal Analysis Conference, Florida, USA, 1984, p. 627634.

Zhu Dechun, Zhang Fang, Jiang Jinhui Identification technology of dynamic load location. Journal of Vibration and Shock, Vol. 31, Issue 1, 2012, p. 2023.

Thite A. N., Thompson D. J. The quantification of structureborne transmission paths by inverse methods. Part1: Improved singular value rejection methods. Journal of Sound and Vibration, Vol. 264, Issue 2, 2003, p. 411431.

Liu Y., Shepared Jr. W. S. Dynamic force identification based on enhanced least squares and total leastsquares schemes in the frequency domain. Journal of Sound and Vibration, Vol. 282, Issue 1, 2005, p. 3760.

Ory H., Glaser H., Holzdeppe D. The reconstruction of forcing function based on measured structural responses. Proceeding of 2nd International Symposium on Aeroelasticity and Structural Dynamics, Aachen, 1985.

Tiwari R., Chakravarthy V. Simultaneous identification of residual unbalances and bearing dynamic parameters from impulse responses of rotorbearing systems. Mechanical Systems and Signal Processing, Vol. 20, 2006, p. 15901614.

Chu Liangcheng, Qu Naisi, Wu Ruifeng Dynamic load identification in time domain. Chain Ocean Engineering, Vol. 5, Issue 3, 1991, p. 279286.

Zhang Yunliang, Lin Gao, Wang Yongxue, et al. An improved method of dynamic load identification in time domain. Chinese Journal of Computational Mechanics, Vol. 21, Issue 2, 2004, p. 209214.

Zhang Yongcheng A Time Domain Method on Two Dimensional Distributed Dynamic Load Identification. Nanjing University of Aeronautics and Astronautics, 2007.

Li Shouju, Liu Yingxi, Ren Mingfa, et al. Identification procedures of vibration loading parameters of hydraulic generator based on improved genetic algorithm. Engineering Mechanics, Vol. 20, Issue 5, 2003, p. 163169.

Wang Haijun, Lian Jijian, Yang Min, et al. Identification of axial dynamic load of a Francis Turbine. Journal of Vibration and Shock, Vol. 26, Issue 4, 2007, p. 123125.

Helio Fiori de Castro, Katia Lucchesi Cavalca, Lucas Ward Franco de Camargo, et al. Identification of unbalance forces by metaheuristic search algorithm. Mechanical Systems and Signal Processing, Vol. 24, 2010, p. 17851789.

Chisari Corrado, Bedon Chiara, Amadio Claudio Dynamic and static identification of baseisolated bridges using genetic algorithms. Engineering Structures, Vol. 102, 2015, p. 8092.

Liu Jie, Meng Xianghua, Jiang Chao Timedomain Galerkin method for dynamic load identification. International Journal for Numerical Methods in Engineering, Vol. 105, Issue 6, 2016, p. 620640.

Jia You, Yang Zhichun, Guo Ning Random dynamic load identification based on error analysis and weighted total least squares method. Journal of Sound and Vibration, Vol. 358, Issue 12, 2015, p. 111123.

Liu Jie, Sun Xingsheng, Han Xu, et al. Dynamic load identification for stochastic structures based on Gegenbauer polynomial approximation and regularization method. Mechanical Systems and Signal Processing, Vol. 56, 2015, p. 3554.

Li Xiaowang, Deng Zhongmin Identification of dynamic loads based on secondorder TaylorSeries expansion method. Shock and Vibration, Vol. 2016, 2016, p. 6461427.

Mao Wentao, Hu Dike, Yan Guirong A new SVM regression approach for mechanical load identification. International Journal of Applied Electromagnetics and Mechamics, Vol. 33, Issue 3, 2010, p. 10011008.

Xu Xun, Ou Jinping Force identification of dynamic systems using virtual work principle. Journal of Sound and Vibration, Vol. 337, 2015, p. 7194.

Xiao Yue, Chen Jian, Li Jiazhu, et al. A joint method of denoising correction and regularization preconditioned iteration for dynamic load identification in time domain. Journal of Vibration Engineering, Vol. 26, Issue 6, 2013, p. 854862.

Lourens E., Reynders E., De Roeck G., et al. An augmented Kalman filter for force identification in structural dynamics. Mechanical Systems and Signal Processing, Vol. 27, 2012, p. 446460.

Nasrellah H. A., Manohar C. S. Finite element method based Monte Carlo filters for structural system identification. Probabilistic Engineering Mechanics, Vol. 26, Issue 2, 2011, p. 294307.

Shang Lin Modal mass computation based on ANSYS finite element analysis. Missiles and Space Vehicles, Vol. 313, Issue 3, 2011, p. 5557.

Zhang Fang Research on Dynamic Load Identification Technology in Time Domain. Nanjing University of Aeronautics and Astronautics,1994.
About this article
This research was supported by the Natural Science Foundation of China (51479165), State Key Laboratory Base of Ecohydraulic Engineering in Arid Area Independent Research Project (2016ZZKT1) and Program 2013 KCT15 of the Shanxi Provincial Key Innovative Research team. Sincere gratitude is extended to the editor and anonymous reviewers for their professional comments and corrections, which significantly improved the presentation of the paper.