Modal balancing of flexible rotors with distributed unbalance

 Abstract: A discrete dynamic balancing method for flexible rotors is proposed in this paper. The polynomial functions are adopted to characterize the mass eccentricity curve of rotor, and the relationship between the parameters of eccentricity curve and vibration response is established by finite element substructure method. Based on the vibration response at the bearing location, the parameters of mass eccentricity curve are identified, which are then decomposed into various modes to achieve the modal balance of the flexible rotor. Upon the proposed method, the simulation of a hollow rotor under cantilever support is performed, which indicates the effectiveness of the discrete dynamic balancing method.


Introduction
Rotating equipment is the throat component of power machinery. In industrial production, rotor unbalance is considered to be the main factor causing vibration of rotating machinery. Severe vibration will not only bring noise, reduce the operating efficiency and service life of machinery and equipment, and even cause dangerous problems. According to the different dynamic balance principles, the flexible rotor dynamic balance can be classified as: influence coefficient method [1][2], which has the disadvantages of starting and stopping the equipment multiple times, modal balance method [3][4], and comprehensive balance method [5], which is developed to integrate the two aforementioned balancing haptic methods. A wide range of scholars have conducted quantities of research on rotor dynamic balance, but most of them are based on lumped mass models [6]. Tiwari et al [7] adopt the concentrated lumped mass model to simultaneously identify the unbalance and the bearing dynamic stiffness parameters. However, very little attention was devoted to the topics of distributed unbalance, support installation position and counterweight plane limitation. Lee and Shih [8] used general transfer matrix method to discretize the unbalanced eccentric curve into Fourier series and analyzed the balance effect. Lee et al [9] proposed that the unbalance of the rotor should be continuously distributed and the lumped mass model is only applicable to thin disks, not to shafts. Due to the complexity of the method, Yang and Lin [10] proposed a finite element method that uses a polynomial curve to fit the unbalanced rotor distribution for the identification of eccentricity of long continuous rotors. Hundal et al [11] proposes a method to balance all modes order by order by establishing the relationship between unbalanced modal components and corrected modal components. Deepthikumar et al [12] established the finite element rotor model, carried out the dynamic balance experiment, and introduced the concept of distributed unbalance quantization for the first time.
In the present study, based on the finite element method, a polynomial curve which is used to fit the distribution unbalance is obtained from the relationship between vibration and unbalance force, then, based on the mode orthogonal theory and the fitted mode function, the cantilever supported rotor is balanced in one single plane. The rotor bearing system is usually consisting of rotating shaft, rigid disks, and bearings. Assume that the system is linear and the damping of the shaft is not considered in the present study.

Shaft model
The rotating shaft is generally discretized into a finite number of elements and the node of an element expressed at the geometric center on the left cross section of the i-th element is shown in Figure 1, where O-XYZ is the fixed coordinate system, and o-xyz the rotating coordinate system attached of the rotor. The origin of the fixed coordinate system and the rotating coordinate system coincide, as shown in Figure 1. The rotor rotates around the Z axis with an angular velocity Ω, and u, v is the node displacement of the rotor in X, Y axis, and θx, θy the node deflection angle of the rotor around the X, Y axis. The mass center Me, geometric center Ge, eccentric component and phase of any cross-section are depicted in Figure 2. The appropriate number of elements depends on the modal order expected to be analyzed and the geometry of the rotor. It is assumed that the cross-sectional shape, dimension, and material constant are uniform in each shaft element. In this paper, the rotor is modelled by using two-nodes Timoshenko beam elements, and two dofs in translation while another two dofs in rotation of each node are considered.  For a shaft element as is shown in Fig.1, the equation of motion can be expressed as: where [M] (e) , [G] (e) , and [K] (e) are the mass, gyroscopic, and stiffness matrix of shaft element, {q} (e) and {f} (e) are the displacement and force vector of the element node, respectively. (A detail of the matrix is given in [13]) Assuming that the unbalance of the rotor is continuous piecewise along the axis of the rotor, and expressed in m-degree polynomial, the unbalance of the rotor in the rotating coordinate system is expressed as follows: a jb x s y s l The gravity factor is not considered in this paper, which is given by [13] for the opposite circumstance. The generalized unbalanced force can be derived from the beam element shape function and virtual work principle.

Rigid disc model
The thin disc is still modeled by adopting the concentrated and lumped mass method, the motion of the rigid disk can be expressed as: where where md, ed, and φd are the mass, disc eccentricity, and eccentric phase of the disc, respectively.

Bearing model
The bearing is modeled with eight stiffness and damping coefficients, and the force on each bearing can be expressed as: where {q}B and {f}B are the displacement and force vector of the bearing, respectively

Equation of motion of the rotor substructure
The equation of motion of the rotor substructure can be obtained by assembling the equation of motion of the shaft element and the disk, expressed as: where are the rotor mass, gyroscopic, and stiffness matrix, respectively.
The overall system is composed of the rotor substructure and the bearing substructure, and the equation of motion can be expressed as: where {q} and {f} are called the overall generalized nodal displacement and force vectors, respectively. Matrices Following the procedure adopted by Yang and Lin [10], the trajectory of the axis in the steady-state response of the lateral vibration of the shaft can be expressed as: In order to simplify the expression, the node displacement of the rotor is expressed as: Correspondingly, the synchronous external force on the rotor can be expressed as: Substituting equation (11) and equation (12) into equation (9) , we obtain:

Balancing methodology
Based on the finite element model established in the previous section, assuming a set of global and local discrete unbalanced polynomial functions, the relationship between the global and local unbalanced polynomial coefficients is established through the node position information, and the rotor unbalanced force is obtained by the local polynomial coefficients, furthermore, the relationship between the vibration response and the global unbalanced polynomial is determined. The mass imbalance is expressed in the form of modal components, and the first-order bending speed dynamic balance is performed on the flexible rotor based on the mode orthogonal theory. The general process of discrete dynamic balancing method is shown in the flow chart below, Figure 3.

Identification of unbalanced coefficient of global distribution
The local eccentricity curve is defined for each shaft element unit, while the global eccentricity curve is throughout the rotor. The local eccentric curve is expressed in the x-z and y-z planes as: Correspondingly, the global eccentric curve can be expressed as: where ai and bi (Ai and Bi) are local (global) eccentric polynomial coefficients in the x-z and y-z planes, respectively. X(z) is the projection of the global eccentricity curve on the x-z plane, and Y(z) the projection of the global eccentricity curve on the y-z plane. s is the local coordinate in a shaft, l the length of the shaft element and z the global axial coordinate. In order to ensure that the eccentric curve coefficient can be solved, the following condition must be satisfied: 8(n+1)>=2(m+ 1), where n is the number of axis element and m is the degree of the polynomial. For ease of expression, the following uses a five polynomial curve as an example to illustrate, the polynomial curve of any order can be obtained in the same way, so the overall eccentric curve and the local eccentric curve are: For the same axis element, in order to match the local eccentric curve and the global eccentric curve, the values and derivatives of the local and global eccentric curves at the node of the element are required to be equal. The following uses the boundary conditions of the x-z plane as an example to illustrate, and the y-z plane is the same: where l=l2-l1 , l1 , l2 are the z coordinate at the node of the element. ''' '' denotes the derivative with respective to z. In order to simplify the expression, the above formula is organized into a matrix form: In order not to lose the generality, the above expression can be extended to a general expression of order m and integrate all axis element, we obtain: where m is an odd number because if m is even, the derivatives up to order of m/2 of the local and global eccentric curves are matched at one node of the element, while only the derivatives up to m/2 are matched at another node. Similarly, the function value and derivative value of the global eccentricity coefficient at the node can be expressed as: where {U5} is the global eccentric curve coefficient vector and {TL} e the transfer matrix between the two vectors.
In order not to lose the generality, the above expression can be extended to a general expression of order m and integrate all axis element: So far, we have determined the coefficient relationship between the global eccentric curve and the local eccentric curve: At the same time, the relationship between the unbalanced force of the shaft element and the local eccentricity coefficient can be determined by Equation (2 Assembling the element equations, we obtain: If the thin disc is considered to have mass eccentricity, the identification coefficient matrix can be expressed as: Correspondingly, the unbalanced force of the disc should be embedded in the position of the node where the disc is located. so, the unbalanced force matrix is modified to: Substituting Equation (25) and (27) into Equation (11), we obtain:  [10]. At the same time, in order to ensure the existence of solution to the equation, the total degree of freedom of the measurement should be greater than 2(m+1)+2d.

Single plane modal balancing
This paper is more interested in modal dynamic balance with limited weight plane, so single-plane modal dynamic balancing method is used for first-order bending valancing, with emphasis to limit the vibration amplitude at location of the bearings.
The rotor mass imbalance distribution can be expressed by the following two methods according to the centroid eccentricity and modal components (the following uses the Y plane as an example to illustrate): where ηi ' is the mass unbalanced modal component coefficient Based on the mode orthogonal theory, the above formula can be expressed as: where Mi is the i-th modal generalized mass. Similarly, the counterweight is expressed as a modal component: where βi is the counterweight modal component coefficient and c is the location of the counterweight. In order to balance a certain mode, the coefficients of mass unbalance and counterweight in this mode should be opposite to each other, so the equation required to balance the i-th modal component is as follows Therefore, the balance weight and phase required to balance the first-order bending speed at the rotor axial position c are given by:

Numerical simulation
In order to explain the above-mentioned distribution unbalance prediction and modal dynamic balance process, this paper establishes a horizontally installed rotor system, which is composed of a steel stepped shaft and rolling bearings as is shown in Figure 4. Among them, the shaft is discretized into ten two-node Timoshenko beam elements and the cross stiffness and damping terms are not considered in rolling bearings. The total length of the shaft is 2m and the counterweight plate to balance the first-order critical speed is in the cantilever state at 11 nodes. The material properties and other parameters of the rotor system with equivalent bearing stiffness is shown in Table 1. The numerical simulation process is completed in Matlab (Ver. 2017b)

Figure 4 Schematic diagram of rotor bearing system
Since the rotor is divided into ten units, the mass matrix [M], generalized damping matrix [C], and stiffness matrix [K] of order 44×44 can be determined, we can then, fit the modal shape data calculated by the finite element model. The mode functions of the x-z and y-z planes can be obtained, which is shown in Figure 6. The first/second-order bending mode function is: Furthermore, the natural frequency of the system can be obtained through Equation (13). which is given in Table 2. Under the condition of cantilever support, the natural frequency of the rotor system will be lower than that under the simple support condition.
The bearing stiffness is exactly the same in the x and y directions, which can be account for the similarity of the modal shape in the x-z and y-z panes. The discrete mode shape data graph and the fitted mode shape function graph are shown in Figure 6.   The vibration at the bearing before balancing is shown in Figure 7. As can be seen from the diagram, the rotor reaches resonance at the first bending frequency. After crossing the first-order critical speed, the vibration increases and the rotor gradually approaches the second-order resonance region. Due to the limitation of the position of the counterweight plane, through equation (25), the amplitude and phase of the counterweight required to balance the first-order bending which can be calculated by balancing methodology are: 23932.625gmm / 77.8449 degree.  From table 3, we can draw that after the single-plane modal balancing, the vibration of the bearing at the first-order bending frequency is reduced from the initial 189/233.7 to 100.9/120.5, which is reduced by about 50%. The effect of dynamic balance is remarkable, and by measuring the vibration at the support, the equation of the eccentricity curve after balancing can be obtained.

Conclusions
The finite element substructure method combined with the m-degree polynomial curve is proposed to establish the discrete unbalanced model of the rotor system in the present study and the relationship between the vibration at the bearing locations and the polynomial curve coefficients, the mass unbalance and the counterweight modal components is analyzed. A polynomial curve of degree 5 is adopted to simulate the distributed unbalance of the hollow rotor in cantilever state and the modal method to dynamically balance the rotor. The effectiveness of the discrete dynamic balancing method and the validity of the model are verified by simulation. The vibration at the bearing locations decreased to about 50% by adding the counterweight obtained from the procedure. The developed programs of discrete dynamic balancing in this paper based on finite element and modal balance theory can be used for on-site dynamic balance analysis of large flexible rotors.