A Balancing Method for Multi-Disc Flexible Rotors without Trial Weights

: Rotor dynamic balancing is a classical problem. Traditional balancing methods such as the inﬂuence coefﬁcient method and the modal balancing method, have low balancing efﬁciency because they need to run many times to add trial weights. Although the model-based balancing method improves the balancing efﬁciency, it cannot accurately identify the position, amplitude and phase of each unbalance fault for rotors with multi-disc structures, so it is difﬁcult to apply it to actual balancing. To solve the above problems, based on the traditional modal balancing theory, this paper deduces that the continuous and isolated unbalance in the rotor-bearing system can be represented by isolated unbalance on several balancing planes approximately. The model-based method is used to identify the above-mentioned equivalent isolated unbalances, and then the corrected mass is added to the balancing planes so as to complete the balance of multiple ﬂexible rotor without trial weights. Considering the practical situation, the proposed balancing method includes two steps: low-speed balancing and high-speed balancing. The proposed balancing method is veriﬁed using three and four-disc rotors. The simulation results show that the balancing method can effectively reduce the vibration of the ﬂexible rotor after low-speed and high-speed balancing, and the amplitude at the measurement point is reduced by 79.74~97.60%, respectively.


Introduction
Rotating machinery is one of the most critical types of machinery in the industry, and its long-term reliable and safe operation is essential. Vibration caused by an imbalanced rotor is the main reason for the abnormal operation of rotating machinery. Rotor dynamic balancing is the main method to reduce the vibration caused by unbalancing fault, and it is an important topic of long-term concern in industry and academia. Numerous scholars have reviewed dynamic balancing methods [1,2]. The classical methods of rotor dynamic balancing mainly include the influence coefficient method and the modal balancing method, which are widely used because of their simplicity and stable performance. However, in the process of balancing, many trial runs are needed, which reduces the economic benefits of the unit. In addition, because of the complexity of the process, the application of the balancing methods based on artificial intelligence [3,4] and empirical technology [5] is limited. To improve the balancing efficiency, academics have proposed numerous trial-weight-free balancing methods [6].
Hundal [7] first proposed a trial-weight-free dynamic balancing method based on modal parameters, which can obtain the accurate correction mass through a series of calculations after obtaining accurate modal and experimental data. Subsequently, many scholars conducted a series of studies on the improvement of the trial-weight-free modal balancing method [8][9][10], but the above methods rely on the modal analysis of the rotor and have certain limitations in use. In order to develop an efficient and stable unbalance parameter identification method, scholars applied the finite element model to the research of the balancing method without trial weight and proposed a model-based unbalance parameter identification method. Zou et al. [11] proposed a novel time-domain unbalanced load identification method with the finite element method together with an augmented Kalman filter algorithm. Zhao et al. [12] proposed a balancing method based on transient characteristics, combined with dynamic load identification techniques, to identify the unbalance parameters of the rotor system. Wang et al. [13] proposed two algorithms called the single direction algorithm and the two orthogonal direction algorithm, in order to identify the unbalance of multi-disc and multi-span rotors from the unbalanced response. Bachschmid et al. [14] proposed a model-based method for identifying multiple faults including unbalance faults, which introduces residual plots as a visualization tool for fault identification and localization. Yao et al. [15] proposed a single-disk rotor unbalance fault identification method based on modal expansion and optimization algorithms, which can detect the location, amplitude, and phase of unbalance faults. In addition, the paper also proposed a two-disk rotor unbalance fault identification method integrating modal expansion, inverse problem, and optimization algorithm. Shrivastava et al. [16] proposed a model-based unbalance fault identification method using a Kalman filter and recursive least squares, which requires the location of the unbalance fault to be known in advance. Sudhakar et al. [17] proposed an improved equivalent load minimization method and vibration minimization method, and the above methods can identify the location, amplitude, and phase of a single unbalance fault with a few measurement points. Chatzisavvas et al. [18] proposed a technique for unbalance identification based on sparse vibration measurements using the equivalent load method, which is capable of identifying the location of multiple unbalance faults without requiring a priori knowledge.
It is usually assumed that the unbalance fault occurs at the rotor disc. The accurate identification of the unbalanced fault at the disc by the model-based method proposed above needs to be based on the accurate identification of positions. When there are more discs, none of the above methods can accurately identify the location, amplitude, and phase of the unbalance, which usually happens in gas turbines, aero-engines, and other rotating equipment with multistage blades. In addition, when balancing the rotor by the above method, the balancing planes must be consistent with the planes where the unbalance faults occurs. However, in practical application, due to many restrictions, the number of balancing planes may be less than that of unbalanced faults planes, thus limiting the application of the model-based dynamic balancing method.
In this paper, a model-based balancing method without trial weight is proposed for the unbalance faults of a multi-disc rotor-bearing system. Based on the modal balance theory, we first transform the rotor unbalance into isolated unbalance on several correction planes, and then identify the above unbalance by using the model-based method. Based on the above technology, a balancing method of flexible rotor with multi-disc is proposed in this paper, which includes two processes: low-speed balancing and high-speed balancing. Numerical examples show the effectiveness and efficiency of the balancing method in reducing the rotor vibration. This paper is organized as follows. In Section 2, an analytical expression is established based on the modal balance theory, which theoretically shows that the unbalance of the rotor distributed in continuous or isolated form can be approximated as the equivalent isolated unbalance on several balancing planes. In Section 3, a model-based unbalance parameter identification method is proposed which can identify the equivalent isolated unbalance on the balancing planes. The balancing process proposed is briefly described in Section 4. In Section 5, the balancing method proposed is simulated and verified, and the factors affecting the balancing performance are discussed. The conclusion is given in Section 6.

Theory
The deformation y(z) of a rotor with bending stiffness EI(z) and continuous mass m(z) under an unbalanced force F can be described by the differential equation [19]: where ω is the angular velocity of the rotor. Unbalance faults are usually caused by the inconsistency between the center of mass and the center of rotation of the rotor. Rotating machinery is important equipment in modern industry. Their unbalanced distribution has two forms: isolated distribution and continuous distribution [20], as shown in Figure 1.

Theory
The deformation ( ) of a rotor with bending stiffness ( ) and continuous mass ( ) under an unbalanced force can be described by the differential equation [19]: where is the angular velocity of the rotor. Unbalance faults are usually caused by the inconsistency between the center of mass and the center of rotation of the rotor. Rotating machinery is important equipment in modern industry. Their unbalanced distribution has two forms: isolated distribution and continuous distribution [20], as shown in Figure 1. When the rotor rotates at an angular speed , the unbalance distributed in isolated and continuous form generates an unbalance force , which is expressed as [19] = ( ) + where the continuous unbalance ( ) can be expanded by the principal mode in the following form where ( ) is the principal mode, which is the solution of Equation (1) when = 0.
is the nth modal component of the unbalanced distribution, and its expression is where is the modal mass of the nth mode and can be expressed as The response of the rotor under the unbalanced force can be expressed in the following form [19] ( , ) = − + 1 ( ) According to Equation (6), we can conclude that a single isolated unbalance can excite all modes ( = 1,2, … ) unless it happens to be located at a node of a certain mode When the rotor rotates at an angular speed ω, the unbalance distributed in isolated and continuous form generates an unbalance force F, which is expressed as [19] where the continuous unbalance u(z) can be expanded by the principal mode in the following form u(z) = ∑ n c n m(z)ϕ n (z) (n = 1, 2, . . .) where ϕ n (z) is the principal mode, which is the solution of Equation (1) when F = 0. c n is the nth modal component of the unbalanced distribution, and its expression is (4) where N n is the modal mass of the nth mode and can be expressed as The response y of the rotor under the unbalanced force F can be expressed in the following form [19] According to Equation (6), we can conclude that a single isolated unbalance u p can excite all modes (n = 1, 2, . . .) unless it happens to be located at a node of a certain mode shape ϕ n z p = 0 without excitation of that mode. According to Equation (3), the continuous unbalance u(z) only excites the modes it contains, e.g., when u(z) = c 2 m(z)ϕ 2 (z), only the second order principal mode is excited.
Rotor unbalance is usually caused by manufacturing errors, assembly errors and material non-uniformity, and its distribution is complex, so it is difficult to accurately measure its distribution characteristics in practice. To solve the above problems, the purpose of this paper is to find a set of unbalance with a simpler distribution which excites a vibration response equal to the vibration response excited by the original unbalance of the system. The isolated unbalances in accordance with the characteristics of the simpler distribution form and the following discussion focuses on whether there is such a set of isolated unbalances.
Suppose there are several isolated unbalances, and the vibration caused by these unbalances is the same as that caused by the original unbalances in the rotor. The unbalance U K at multiple axial positions z k is used to represent the above-mentioned isolated unbalance. According to Equation (6), to make the vibration the same, the unbalance U K should meet the following conditions: When the angular velocity of rotor is lower than ω N , it is generally considered that the influence of modes above the N order can be ignored, i.e., it is considered that ω 2 / ω 2 n − ω 2 in Equation (6) tends to 0 when n > N, at which time Equation (7) becomes Further expanding Equation (8), the relationship between the equivalent unbalance U K and the original unbalance of the system can be expressed by Equation (9) Obviously, if and only if K = N, Equation (9) has a solution. Therefore, we can conclude that when the rotor runs below ω N , the original unbalance in the rotor can be approximately equivalent to the isolated unbalances in K balancing planes. It should be noted that any of these K planes should not be located at the node of any vibration mode, which will lead to the number K of unknown parameters in Equation (9) being less than the number N of equations, and thus the equivalent unbalance U K can not be solved.

Identification of Equivalent Unbalance
Due to the unknown of the right side in Equation (9), the traditional modal balancing method requires several runs for adding trial weights to complete the balancing of the rotor, which greatly reduces the efficiency. To improve the efficiency, the model-based unbalance parameter identification method has been greatly developed in recent years. It can estimate unbalance parameters according to the response signals collected in one run. This section introduces the identification of the equivalent unbalance U K by using the model-based method.
The equation of motion of a flexible rotor can be written as [21]: where M, C, K, and G are the mass, damping, stiffness, and gyroscopic matrices, respectively, q is the displacement vector of the nodes, and f is the unbalanced force vector. To improve solution efficiency and even the unsolvable solution due to the large dimensionality of the finite element model, the reduced-order model is used to establish the equivalent unbalance solution model. The Iterated Improved Reduced System (IIRS) method proposed by Choi et al. [22] is an effective technique for obtaining condensed matrices. We can obtain the exact characteristic properties of the system from the condensed matrices without consuming expensive calculation costs. However, the condensed matrices of the gyroscopic matrices cannot be obtained by using the above method, so the improved IIRS method [23] is used to reduce the order of the model of the rotor-bearing system. This method must first select the degrees of freedom (dof) to be retained. In this paper, the translational dof of the node where the equivalent unbalances and the measuring sensor is located are selected as the retained dof. The reduced-order equation of Equation (10) is [23].
The state space form of Equation (11) is [24] .
is the input matrix, where S f is the selection matrix of unbalanced forces, the elements of the translational dof of the nodes where the equivalent unbalances are located are usually set to 1 and the others are set to 0. Y(t) is the output vector. The state output matrix H c is the selection matrix of the displacement response. D c is the direct transfer matrix. In this paper, it is the zero matrix because only the displacement is measured. The superscript c denotes the continuous form of the matrix.
The discrete form (denoted by the superscript d) of Equation (12) can be expressed as [25] .
where ∆t is the sampling time for discretization. For the zero initial condition, Equation (13) can be written as the following moving average model [26] where R is called the Markov parameter [27], where Y = [Y(0) Y(1) · · · Y(n − 1)] T is the response vector at the measurement point, · · · f (n − 1)] T is the unbalanced forec vector at the disk, n is the total number of sampling points. The expression of R is as follows: Equation (15) is established to solve for the equivalent unbalanced forces. We can obtain the matrix R when the parameters of the rotor-bearing system are confirmed. We can measure the displacement response Y(j), j = 0, 1, . . . , n − 1, and arrange it to obtain the vector Y. By solving Equation (15), we can obtain the vector F, which can be decomposed to obtain the equivalent unbalanced force at each moment. When a displacement transducer is used to measure the rotor vibration, all elements of the matrix R 0 are zero. According to the Picard condition [28], the solution of F in Equation (15) is an ill-posed inverse problem.
According to the literature [24], the solution of ill-posed inverse problem usually has the following form where α is the regularization parameter selected in the solution process, and in this paper the L-curve criterion [29] is used to determine the parameter α.
After the vector F is decomposed to obtain the equivalent unbalanced forces at each moment, the amplitude and phase of the equivalent unbalanced force are extracted by using the spectrum correction method [30]. The spectral correction method is an improved FFT method, which can improve the extraction accuracy of signal frequency, amplitude and phase.

Balancing Process
Considering the actual working condition, the unbalance may cause excessive vibration at the critical speed, which will lead to the problem that the rotor cannot safely pass the critical speed, so the low speed balancing is needed in advance. Therefore, the balancing process is set up as two processes: low-speed and high-speed balancing. The process is shown in Figure 2. Since the processes of low-speed balance is the same as that of high-speed balance, and only the rotational speed is different, this paper takes low-speed balance as an example to illustrate the process of this balancing method.
Equation (15) is established to solve for the equivalent unbalanced forces. We can obtain the matrix when the parameters of the rotor-bearing system are confirmed. We can measure the displacement response ( ), = 0,1, … , − 1, and arrange it to obtain the vector . By solving Equation (15), we can obtain the vector , which can be decomposed to obtain the equivalent unbalanced force at each moment. When a displacement transducer is used to measure the rotor vibration, all elements of the matrix are zero. According to the Picard condition [28], the solution of in Equation (15) is an ill-posed inverse problem.
According to the literature [24], the solution of ill-posed inverse problem usually has the following form where is the regularization parameter selected in the solution process, and in this paper the L-curve criterion [29] is used to determine the parameter .
After the vector is decomposed to obtain the equivalent unbalanced forces at each moment, the amplitude and phase of the equivalent unbalanced force are extracted by using the spectrum correction method [30]. The spectral correction method is an improved FFT method, which can improve the extraction accuracy of signal frequency, amplitude and phase.

Balancing Process
Considering the actual working condition, the unbalance may cause excessive vibration at the critical speed, which will lead to the problem that the rotor cannot safely pass the critical speed, so the low speed balancing is needed in advance. Therefore, the balancing process is set up as two processes: low-speed and high-speed balancing. The process is shown in Figure 2. Since the processes of low-speed balance is the same as that of highspeed balance, and only the rotational speed is different, this paper takes low-speed balance as an example to illustrate the process of this balancing method.   Firstly, the balancing plane and the sensor position are selected, and then the low-speed balance shaft speed is selected to establish a reduced-order model. Secondly, construct the matrix R of Equation (15) based on the reduced-order model. Thirdly, the vector Y of Equation (15) can be obtained by collecting and assembling the vibration data at the balanced rotating speed, and the equivalent unbalance force of the rotor-bearing system can be obtained by solving Equation (15) by the regularization method. Fourth, the amplitude and phase of the equivalent unbalance force at the balancing planes can be obtained by spectral correction. Finally, the low-speed balance of the rotor can be achieved by adding the same amount of unbalance at the opposite position of the balancing planes.

Numerical Simulation Verification
In this section, the simulation examples of two rotor-bearing systems (with three and four discs, respectively) are used to verify the effectiveness of the proposed balancing method.

Balancing of 3-Disc Rotor
The rotor bearing system with three discs shown in Figure 3a is selected to verify the proposed balancing method. The density of the rotor is 7850 kg/m 3 , the elastic modulus is 2.06 × 10 11 N/m 2 and the Poisson ratio is 0.3. The diameter of rotor is 0.025 m and the total length is 0.72 m. The rotor-bearing system has three discs D1, D2, and D3 with the same parameters, whose mass are 8.9 kg, diameter inertia is 0.03196 kg·m 2 , and pole inertia is 0.06392 kg·m 2 . The distances of the three discs from the left end of the rotor are 0.24 m, 0.32 m and 0.48 m, respectively. The two bearings of the system have the same parameters and are set to be anisotropic to simulate the actual bearings. The stiffness in x and y directions is 7.5 × 10 8 N/m and 8 × 10 8 N/m and damping in x and y directions is 2 × 10 3 N·s/m and 1.8 × 10 3 N·s/m, respectively. The rotor is simulated as 18 identical Timoshenko beams, each beam has two nodes, and each node has 4 dofs: translation and rotation in X and Y directions, respectively. The finite element model has 19 nodes and 76 dofs. Three discs are located at nodes 6, 9 and 13, and two bearings are located at nodes 1 and 19. Displacement sensors are located in the x and y directions at nodes 3 and 17, respectively, to collect the vibration displacement response of the rotor. Firstly, the balancing plane and the sensor position are selected, and then the lowspeed balance shaft speed is selected to establish a reduced-order model. Secondly, construct the matrix of Equation (15) based on the reduced-order model. Thirdly, the vector of Equation (15) can be obtained by collecting and assembling the vibration data at the balanced rotating speed, and the equivalent unbalance force of the rotor-bearing system can be obtained by solving Equation (15) by the regularization method. Fourth, the amplitude and phase of the equivalent unbalance force at the balancing planes can be obtained by spectral correction. Finally, the low-speed balance of the rotor can be achieved by adding the same amount of unbalance at the opposite position of the balancing planes.

Numerical Simulation Verification
In this section, the simulation examples of two rotor-bearing systems (with three and four discs, respectively) are used to verify the effectiveness of the proposed balancing method.

Balancing of 3-Disc Rotor
The rotor bearing system with three discs shown in Figure 3a is selected to verify the proposed balancing method. The density of the rotor is 7850 kg/m , the elastic modulus is 2.06 × 10 11 N/m and the Poisson ratio is 0.3. The diameter of rotor is 0.025 m and the total length is 0.72 m. The rotor-bearing system has three discs D1, D2, and D3 with the same parameters, whose mass are 8.9 kg, diameter inertia is 0.03196 kg • m , and pole inertia is 0.06392 kg • m . The distances of the three discs from the left end of the rotor are 0.24 m, 0.32 m and 0.48 m, respectively. The two bearings of the system have the same parameters and are set to be anisotropic to simulate the actual bearings. The stiffness in x and y directions is 7.5 × 10 8 N/m and 8 × 10 8 N/m and damping in x and y directions is 2 × 10 3 N•s/m and 1.8 × 10 3 N•s/m, respectively. The rotor is simulated as 18 identical Timoshenko beams, each beam has two nodes, and each node has 4 dofs: translation and rotation in X and Y directions, respectively. The finite element model has 19 nodes and 76 dofs. Three discs are located at nodes 6, 9 and 13, and two bearings are located at nodes 1 and 19. Displacement sensors are located in the x and y directions at nodes 3 and 17, respectively, to collect the vibration displacement response of the rotor.  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Figure 3b shows the mode shapes of the first two orders mode of the system shown in Figure 3a. Figure 3c shows the Campbell diagram of this system. Table 1 shows the first and second order synchronous forward whirl critical speed (FWCS) and synchronous back whirl critical speed (BWCS). The rotor balancing method proposed in this paper is verified by numerical simulation with the rotor operating lower than the critical speed as an example. In the current work, two discs located at nodes 6 and 13 are selected as balancing planes and nodes 3 and 17 are selected as the measurement position of the sensors. And the displacement responses in x and y direction are used to calculate the equivalent unbalance. According to the Section 3, before identifying the equivalent unbalance parameters, the reduced-order model of the finite element model must be established first, and the translational dof of nodes [3,6,13,17] are chosen as the dof of the reduced-order model. Two unbalanced distribution are used to verify the reliability of the proposed method. The amplitude and phase of the unbalance are shown in Table 2. The unbalanced masses of all three discs were placed at 0.1 m from the center of the disc. The effectiveness of the balancing method is verified when the unbalanced distribution is in case 1 in Table 2. Because the unbalance may cause the rotor to vibrate too much at the critical speed, and then the rotor cannot safely cross the critical speed, it is necessary to balance at low speed. In this paper, the speed for low-speed balancing is selected as 40 Hz. The equivalent unbalance is identified by the identification method proposed in Section 3. The length of the displacement data needs to be specified to obtain the vector in Equation (15). In this paper, the equivalent unbalanced force is reconstructed by using the vibration displacement signal of 30 cycles and 30 sampling points in each cycle. The iden-  Figure 3b shows the mode shapes of the first two orders mode of the system shown in Figure 3a. Figure 3c shows the Campbell diagram of this system. Table 1 shows the first and second order synchronous forward whirl critical speed (FWCS) and synchronous back whirl critical speed (BWCS). The rotor balancing method proposed in this paper is verified by numerical simulation with the rotor operating lower than the critical speed ω 2 as an example. In the current work, two discs located at nodes 6 and 13 are selected as balancing planes and nodes 3 and 17 are selected as the measurement position of the sensors. And the displacement responses in x and y direction are used to calculate the equivalent unbalance. According to the Section 3, before identifying the equivalent unbalance parameters, the reduced-order model of the finite element model must be established first, and the translational dof of nodes [3,6,13,17] are chosen as the dof of the reduced-order model. Two unbalanced distribution are used to verify the reliability of the proposed method. The amplitude and phase of the unbalance are shown in Table 2. The unbalanced masses of all three discs were placed at 0.1 m from the center of the disc. The effectiveness of the balancing method is verified when the unbalanced distribution is in case 1 in Table 2. Because the unbalance may cause the rotor to vibrate too much at the critical speed, and then the rotor cannot safely cross the critical speed, it is necessary to balance at low speed. In this paper, the speed for low-speed balancing is selected as 40 Hz. The equivalent unbalance is identified by the identification method proposed in Section 3. The length of the displacement data needs to be specified to obtain the vector Y in Equation (15). In this paper, the equivalent unbalanced force is reconstructed by using the vibration displacement signal of 30 cycles and 30 sampling points in each cycle. The identification results of equivalent unbalance at discs D1 and D3 are shown in Table 3.  To test the effect of low-speed balancing, the steady-state unbalance response before and after low-speed balancing is calculated using the full-order model. The response is calculated using the full-order model in the subsequent part of this paper and will not be emphasized in the future. Figure 4 shows the spatial shape of the rotor before and after balancing. The amplitude of the displacement in the x-direction at nodes 3 and 17 before and after balancing at 40 Hz is shown in Figure 5. The amplitude in the x-direction of node 3 is reduced from 2.218 × 10 −5 m to 4.139 × 10 −6 m, and the amplitude in the x-direction of node 17 is reduced from 2.364 × 10 −5 m to 5.532 × 10 −6 m, with the obvious effect of vibration reduction.  To test the effect of low-speed balancing, the steady-state unbalance response before and after low-speed balancing is calculated using the full-order model. The response is calculated using the full-order model in the subsequent part of this paper and will not be emphasized in the future. Figure 4 shows the spatial shape of the rotor before and after balancing. The amplitude of the displacement in the x-direction at nodes 3 and 17 before and after balancing at 40 Hz is shown in Figure 5. The amplitude in the x-direction of node 3 is reduced from 2.218 × 10 −5 m to 4.139 × 10 −6 m, and the amplitude in the x-direction of node 17 is reduced from 2.364 × 10 −5 m to 5.532 × 10 −6 m, with the obvious effect of vibration reduction.
Before low-speed balancing After low-speed balancing  In addition, as shown in Figure 6, the spatial shape of the rotor at 120 Hz before and after low-speed balance is also calculated. Figure 7 shows the amplitude of the rotor at nodes 3 and 17 in the x-direction before and after low-speed balancing at 120 Hz. The amplitude of node 3 x-direction is reduced from 3.744 × 10 −5 m to 2.402 × 10 −5 m, and the amplitude of node 17 x-direction is reduced from 3.099 × 10 −5 m to 2.116 × 10 −5 m. It can be found that after low-speed balancing, the vibration response of the rotor decreases at high speed, but it is not as significant as that at low speed. In addition, as shown in Figure 6, the spatial shape of the rotor at 120 Hz before and after low-speed balance is also calculated. Figure 7 shows the amplitude of the rotor at nodes 3 and 17 in the x-direction before and after low-speed balancing at 120 Hz. The amplitude of node 3 x-direction is reduced from 3.744 × 10 −5 m to 2.402 × 10 −5 m, and the amplitude of node 17 x-direction is reduced from 3.099 × 10 −5 m to 2.116 × 10 −5 m. It can be found that after low-speed balancing, the vibration response of the rotor decreases at high speed, but it is not as significant as that at low speed.
Before low-speed balancing After low-speed balancing  In addition, as shown in Figure 6, the spatial shape of the rotor at 120 Hz before and after low-speed balance is also calculated. Figure 7 shows the amplitude of the rotor at nodes 3 and 17 in the x-direction before and after low-speed balancing at 120 Hz. The amplitude of node 3 x-direction is reduced from 3.744 × 10 −5 m to 2.402 × 10 −5 m, and the amplitude of node 17 x-direction is reduced from 3.099 × 10 −5 m to 2.116 × 10 −5 m. It can be found that after low-speed balancing, the vibration response of the rotor decreases at high speed, but it is not as significant as that at low speed.
Before low-speed balancing After low-speed balancing

High-Speed Balancing
In order to further reduce the vibration of the rotor-bearing system, the rotor is balanced at high speed, and the speed of high-speed balancing is set to 120 Hz. After the correction mass for the low-speed balancing is added, the original unbalanced distribu-

High-Speed Balancing
In order to further reduce the vibration of the rotor-bearing system, the rotor is balanced at high speed, and the speed of high-speed balancing is set to 120 Hz. After the correction mass for the low-speed balancing is added, the original unbalanced distribution becomes as shown in the previous column of Table 4. The equivalent unbalances at discs D1 and D3 are identified by the identification method introduced in Section 3. The identification results are shown in the last column of Table 4. The system is balanced at high speed by adding correction masses of 18.27 g and 18.13 g at −96.05 • and −60.15 • of discs D1 and D3, respectively. The steady-state response of the rotor before and after high-speed balancing is calculated. The spatial shape of the rotor is shown in Figure 8, and the amplitude of the rotor in the x-direction at nodes 3 and 17 before and after highspeed balancing are shown in Figure 9. The amplitude in x-direction of node 3 is reduced from 2.402 × 10 −5 m to 7.584 × 10 −6 m, and the amplitude in x-direction of node 17 is reduced from 2.116 × 10 −5 m to 5.786 × 10 −6 m, with an obvious effect of vibration reduction. Figure 10 shows the amplitude of the rotor's speed-up response at node 3 and 17 x-direction after low-speed and high-speed balancing when the unbalanced distribution is case 1 in Table 2. By comparing the amplitudes, it can be found that the balancing method proposed can effectively reduce the vibration amplitude of the rotor. Before high-speed balancing After high-speed balancing      To verify the effectiveness of the proposed method for different unbalanced distributions, the unbalanced distribution is set to case 2 in Table 2. The proposed method is used to balance the rotor-bearing system at low and high speeds. To save the space of the paper, the spatial shape of the rotor and the response at nodes 3 and 17 at the balancing speed are not listed. Figure 11 shows the amplitude of the rotor's speed-up response at node 3 and 17 x-direction after low-speed and high-speed balancing when the unbalanced distribution is case 2 in Table 2. By comparing the amplitudes, it can be found that the balancing method proposed can effectively reduce the vibration amplitude of the rotor. Table 5 summarizes the vibration amplitudes of nodes 3 and 17 in the x-direction at 120 Hz before and after balancing for different unbalance excitations (Case 1 and Case 2). The results show that the vibration amplitude of the rotor decreases by 79.74% to 95.65%, which proves the effectiveness of the balancing method.
to balance the rotor-bearing system at low and high speeds. To save the space of the paper, the spatial shape of the rotor and the response at nodes 3 and 17 at the balancing speed are not listed. Figure 11 shows the amplitude of the rotor's speed-up response at node 3 and 17 x-direction after low-speed and high-speed balancing when the unbalanced distribution is case 2 in Table 2. By comparing the amplitudes, it can be found that the balancing method proposed can effectively reduce the vibration amplitude of the rotor.  Table 5 summarizes the vibration amplitudes of nodes 3 and 17 in the x-direction at 120 Hz before and after balancing for different unbalance excitations (Case 1 and Case 2). The results show that the vibration amplitude of the rotor decreases by 79.74% to 95.65%, which proves the effectiveness of the balancing method. Table 5. Amplitude comparison at 120 Hz before and after low speed and high speed balancing.

Balancing of 4-Disc Rotor
The 4-disc rotor used for the simulation in this section is based on the 3-disc rotor in Section 5.1, with an additional disc at node 16, and its parameters are consistent with those of other discs. Figure 12 shows the finite element model, the modal shapes of the first two modes, and Campbell diagrams of the four-disc rotor bearing system. The first and second order FWCS and BWCS are listed in Table 6. (a) Figure 11. Comparison of the steady-state response of the rotor in the node 3 and 17 x-directions before and after balancing (case 2). Table 5. Amplitude comparison at 120 Hz before and after low speed and high speed balancing.

Balancing of 4-Disc Rotor
The 4-disc rotor used for the simulation in this section is based on the 3-disc rotor in Section 5.1, with an additional disc at node 16, and its parameters are consistent with those of other discs. Figure 12 shows the finite element model, the modal shapes of the first two modes, and Campbell diagrams of the four-disc rotor bearing system. The first and second order FWCS and BWCS are listed in Table 6.
bution is case 2 in Table 2. By comparing the amplitudes, it can be found that the balancing method proposed can effectively reduce the vibration amplitude of the rotor. Figure 11. Comparison of the steady-state response of the rotor in the node 3 and 17 x-directions before and after balancing (case 2). Table 5 summarizes the vibration amplitudes of nodes 3 and 17 in the x-direction at 120 Hz before and after balancing for different unbalance excitations (Case 1 and Case 2). The results show that the vibration amplitude of the rotor decreases by 79.74% to 95.65%, which proves the effectiveness of the balancing method. Table 5. Amplitude comparison at 120 Hz before and after low speed and high speed balancing.

Balancing of 4-Disc Rotor
The 4-disc rotor used for the simulation in this section is based on the 3-disc rotor in Section 5.1, with an additional disc at node 16, and its parameters are consistent with those of other discs. Figure 12 shows the finite element model, the modal shapes of the first two modes, and Campbell diagrams of the four-disc rotor bearing system. The first and second order FWCS and BWCS are listed in Table 6.  The comparison between Tables 1 and 6 shows that the addition of disc 4 leads to the decreases of the critical speed of the system. Therefore, the balancing speeds are set to 35 and 90 Hz in this section. The balancing goal of this section is the same as that of the threedisc rotor, that is, balancing the rotor-bearing system running below the second order critical speed . The parameters selected in this section (including the nodes of the reduced-order model, the length of the vibration signal, and the number of sampling points per cycle) are the same as those in Section 5.1.
In this section, the unbalanced distribution is set to two cases, case 3 and case 4, as shown in Table 7, to verify the balancing ability of the balancing method for rotors with different unbalanced distribution. The balancing process is the same as that of the 3-disc rotor described in Section 5.1, and the detailed balancing process will not be described again. After low-speed and high-speed balancing, Figures 13 and 14 show the steady-state displacement of the rotor at nodes 3 and 17 x-direction at each speed, respectively. Table  8 shows the comparison of the displacement of the rotor at nodes 3 and 17 in the x-direction at the balancing speed of 90 Hz before and after balancing. It can be found that the vibration amplitude of the rotor is significantly reduced after the low and high-speed balancing.   The comparison between Tables 1 and 6 shows that the addition of disc 4 leads to the decreases of the critical speed of the system. Therefore, the balancing speeds are set to 35 and 90 Hz in this section. The balancing goal of this section is the same as that of the three-disc rotor, that is, balancing the rotor-bearing system running below the second order critical speed ω 2 . The parameters selected in this section (including the nodes of the reduced-order model, the length of the vibration signal, and the number of sampling points per cycle) are the same as those in Section 5.1.
In this section, the unbalanced distribution is set to two cases, case 3 and case 4, as shown in Table 7, to verify the balancing ability of the balancing method for rotors with different unbalanced distribution. The balancing process is the same as that of the 3-disc rotor described in Section 5.1, and the detailed balancing process will not be described again. After low-speed and high-speed balancing, Figures 13 and 14 show the steady-state displacement of the rotor at nodes 3 and 17 x-direction at each speed, respectively. Table 8 shows the comparison of the displacement of the rotor at nodes 3 and 17 in the x-direction at the balancing speed of 90 Hz before and after balancing. It can be found that the vibration amplitude of the rotor is significantly reduced after the low and high-speed balancing.

Discussion
In this paper, an efficient balancing method for multi-disc rotors is proposed. Based on the modal balance theory, it is first demonstrated that the original unbalance in the rotor can be approximated as isolated unbalance on multiple balancing planes, then the model-based identification method is adopted to identify the above equivalent unbalance and, finally, the rotor balancing is completed by adding the same unbalance mass on the reverse position of the balancing planes. Considering that the actual rotor may not exceed the critical speed due to excessive vibration, the balancing process is set to low-speed balancing and high-speed balancing. From the comparison of the amplitudes of each speed before and after balancing shown in Figures 10, 11 13 and 14, the proposed balancing

Discussion
In this paper, an efficient balancing method for multi-disc rotors is proposed. Based on the modal balance theory, it is first demonstrated that the original unbalance in the rotor can be approximated as isolated unbalance on multiple balancing planes, then the model-based identification method is adopted to identify the above equivalent unbalance and, finally, the rotor balancing is completed by adding the same unbalance mass on the reverse position of the balancing planes. Considering that the actual rotor may not exceed the critical speed due to excessive vibration, the balancing process is set to low-speed balancing and high-speed balancing. From the comparison of the amplitudes of each speed before and after balancing shown in Figures 10, 11 13 and 14, the proposed balancing

Discussion
In this paper, an efficient balancing method for multi-disc rotors is proposed. Based on the modal balance theory, it is first demonstrated that the original unbalance in the rotor can be approximated as isolated unbalance on multiple balancing planes, then the model-based identification method is adopted to identify the above equivalent unbalance and, finally, the rotor balancing is completed by adding the same unbalance mass on the reverse position of the balancing planes. Considering that the actual rotor may not exceed the critical speed due to excessive vibration, the balancing process is set to low-speed balancing and high-speed balancing. From the comparison of the amplitudes of each speed before and after balancing shown in Figures 10, 11, 13 and 14, the proposed balancing method effectively reduces the rotor vibration. From the results in Tables 5 and 8, it can be seen that the amplitude at the measuring point decreases by 79.74% to 97.60%.
However, it should be noted that, in the derivation of Equations (7) to (8), we ignore the influence of modes above the Nth order, but in practice, higher order modes always affect the rotor vibration response at low speeds, although this influence is small. The above reasons will lead to errors of the equivalent unbalances calculated by Equation (9), which will reduce the balancing effect. On the other hand, the accuracy of the equivalent unbalance identification method introduced in Section 3 will also affect the balancing effect of the rotor to some extent. In extension, Equation (15) is constructed based on the reducedorder model, and the prediction accuracy of the unbalance response of the reduced-order model determines the identification accuracy of the equivalent unbalance, which in turn affects the balancing effect of this method.
To illustrate the effect of the reduced-order model, the unbalance response of the 3-disk rotor in Section 5.1 is predicted by the reduced-order model and the full-order model, respectively. Figure 15a shows the comparison of the unbalance responses under the excitation of unbalance (1 × 10 −3 kg-m ∠ 135 • at node 13), and Figure 15b shows the relative error of the unbalance response predicted by the reduced-order model. It can be seen that the unbalance response predicted by the reduced-order model at different speeds is different from that predicted by the full-order model. This deviation can lead to the identification error of the equivalent unbalanced forces, which will affect the balancing effect of the proposed method. It is also to be noted that the maximum deviation occurs near the first-order critical velocity, which should not be chosen as balancing speed, otherwise, it will cause a large distortion in the equivalent unbalance identification. This implies that a more accurate equivalent unbalances identification method will help to improve the balancing effect of the rotor.  Tables 5 and 8, it can be seen that the amplitude at the measuring point decreases by 79.74% to 97.60%. However, it should be noted that, in the derivation of Equations (7) to (8), we ignore the influence of modes above the th order, but in practice, higher order modes always affect the rotor vibration response at low speeds, although this influence is small. The above reasons will lead to errors of the equivalent unbalances calculated by Equation (9), which will reduce the balancing effect. On the other hand, the accuracy of the equivalent unbalance identification method introduced in Section 3 will also affect the balancing effect of the rotor to some extent. In extension, Equation (15) is constructed based on the reduced-order model, and the prediction accuracy of the unbalance response of the reduced-order model determines the identification accuracy of the equivalent unbalance, which in turn affects the balancing effect of this method.
To illustrate the effect of the reduced-order model, the unbalance response of the 3disk rotor in Section 5.1 is predicted by the reduced-order model and the full-order model, respectively. Figure 15a shows the comparison of the unbalance responses under the excitation of unbalance (1 × 10 −3 kg-m ∠ 135°at node 13), and Figure 15b shows the relative error of the unbalance response predicted by the reduced-order model. It can be seen that the unbalance response predicted by the reduced-order model at different speeds is different from that predicted by the full-order model. This deviation can lead to the identification error of the equivalent unbalanced forces, which will affect the balancing effect of the proposed method. It is also to be noted that the maximum deviation occurs near the first-order critical velocity, which should not be chosen as balancing speed, otherwise, it will cause a large distortion in the equivalent unbalance identification. This implies that a more accurate equivalent unbalances identification method will help to improve the balancing effect of the rotor.

Conclusions
This paper proposes a multi-disc flexible rotor balancing method without trial weights. The paper first explains, based on the modal balance theory, that when the rotor operates below the th order critical speed , the unbalance, no matter how complex it is, can be approximated as isolated unbalances on balancing planes. The paper also suggests a model-based identification method to identify the equivalent isolated unbalances on balancing planes. Based on the above theory and method, the paper proposes a balancing method that includes two steps of low-speed balancing and high-speed balancing to reduce the rotor vibration caused by the unbalance. Two rotor bearing systems with three and four disks, respectively, are used to verify the proposed balancing method. The simulation results show that the amplitude in the x-direction at the measurement

Conclusions
This paper proposes a multi-disc flexible rotor balancing method without trial weights. The paper first explains, based on the modal balance theory, that when the rotor operates below the Nth order critical speed ω N , the unbalance, no matter how complex it is, can be approximated as isolated unbalances on N balancing planes. The paper also suggests a model-based identification method to identify the equivalent isolated unbalances on N balancing planes. Based on the above theory and method, the paper proposes a balancing method that includes two steps of low-speed balancing and high-speed balancing to reduce the rotor vibration caused by the unbalance. Two rotor bearing systems with three and four disks, respectively, are used to verify the proposed balancing method. The simulation results show that the amplitude in the x-direction at the measurement point decreases by 79.74~97.60% after balancing. Meanwhile, the paper analyzes the factors affecting the performance of the balancing method and concludes that the balancing effect can be further enhanced by improving the accuracy of the identification method.
The proposed balancing method can be easily applied to rotor bearing systems operating at higher speeds, although this paper verifies the balancing method by using rotors operating below second order critical speed. In the future, we will carry out further research on the improvement of the accuracy of the unbalanced parameter identification method and the experimental confirmation of the balancing method.