Modal Balancing of Warped Rotors without Trial Runs Using the Numerical Assembly Technique

Abstract

The increasing use of high-speed machinery leads to a growing demand for efficient balancing methods for flexible rotors. Conventional balancing methods are costly and time-consuming since they require multiple trial runs. For this reason, recent research focuses on model-based balancing methods, which substitute measurements with simulations. This work presents and examines a model-based modal balancing method, which utilizes the Numerical Assembly Technique (NAT) for the in situ balancing of warped rotors with flexible behaviour. NAT is a successive modification of discrete–continuous modelling that leads to analytical harmonic solutions and is very computationally efficient. In this version of NAT, internal damping is also included with a viscoelastic material model using fractional time derivatives. The modal balancing procedure is adapted to handle measurements outside of the critical speeds and the effect of the pre-bend on the rotor. The accuracy of the simulations is shown by comparing measured mode shapes and eigenvalues with values calculated with NAT. Furthermore, the first two modes of a rotor test bed are successfully balanced without trial runs.

1. Introduction

Imbalance and bow are the most common causes of excessive rotor vibrations, which can be harmful to the machinery and its surroundings. Techniques for the balancing of rotors with flexible behaviour, like modal balancing [1] or the influence coefficient method [2], are necessary to reduce imbalance if the operational speed exceeds 70% of the first critical speed [3]. Additional consideration needs to be taken if the rotor is warped due to thermal distortion, gravity sag or mechanical bow. The balancing procedure can be either performed independently or used to compensate for the effects of the bow at certain spin speeds [4]. While the modal balancing method is efficient and accurate up to higher modes, it is necessary to know the mode shapes of the system to calculate the orthogonal weight sets [5]. Orthogonal weight sets can be found either intuitive [6], experimental [7] or with model-based methods [8].

Modern balancing techniques for rotors with flexible behaviour aim to improve upon both traditional methods in various ways. The unified balancing approach combines the advantages of influence coefficient and modal balancing technique [9]. Transient rotor balancing techniques use run-up data to speed up and improve the balancing process [10]. Another approach, often referred to as homologous information fusion technology, is to use multiple types of information besides vibration measurement, for example, flow, pressure and temperature for balancing [11]. Advances were also made in the balancing of non-standard rotors, e.g., warped rotors [12] or rotor-bearing systems with non-linear behaviour [13,14]. The focus of many modern balancing methods is to substitute measurements with simulations. Recent examples of model-based balancing techniques were proposed by Yun et al. [15], Carvalho et al. [16] and Nordmann et al. [17,18].

This paper investigates the application of a model-based modal balancing technique for warped rotors utilising the Numerical Assembly Technique (NAT). NAT is an analytic method used for determining the steady-state harmonic response, mode shapes and eigenvalues of rotor-bearing systems. Upon its introduction in 1999 by Wu and Chou [19,20], it was only capable of solving very simple beam models. Since then, NAT has been extended significantly, especially in the field of rotor dynamics. Chen and Wu [21,22] applied NAT to calculate the natural frequencies of uniform Euler–Bernoulli beams. Lin [23,24] included multiple lumped masses, multiple-pinned supports, and rotary inertias. In their study [25], Wang et al. incorporated the Timoshenko beam theory into NAT and investigated the effects of shear coefficients and slenderness ratios. Yesilce [26,27] studied the effects of axial forces on Timoshenko and Reddy–Bickford beams. NAT was extended to rotating systems in 2014, allowing the calculation of forward and backward whirling speeds [28]. Farghaly [29] utilized NAT to investigate Timoshenko beams on elastic supports and Vaz [30] examined the effects of various material and geometric discontinuities using NAT. Klanner et al. [31] expanded NAT to encompass distributed loading, imbalance, and fractional derivative material damping. Within recent years, Quinz et al. presented a modal balancing technique [32] and an influence coefficient balancing technique [33], where measurements were substituted with NAT simulations.

The main aim of this work is to adjust the modal balancing method based on NAT first presented in [32] to be used under real-world conditions and verify its utility experimentally. The novelty of this work is:

• Investigation of the optimal adaption of the model-based modal balancing procedure for warped rotors.
• Adaption of the procedure to measurements outside of critical speeds.
• Experimental verification of the proposed procedure.

Numerical analyses demonstrate the functionality, experiments are performed on a test bed to compare measured and calculated values and the first two modes of a test setup are balanced with the proposed method.

2. Materials and Methods

In this section, the proposed model-based modal balancing method—where the imbalance responses, eigenvalues and mode shapes are calculated with a simulation of the rotor-bearing system using NAT—is described. Only a vibration measurement of the initial rotor without trial weights is necessary for balancing. As long as the bow of the rotor is significantly smaller than its diameter, it does not influence the accuracy of the NAT simulation but is considered in the subsequent modal balancing process.

2.1. Numerical Assembly Technique

NAT employs a continuum model to derive analytic solutions when imbalance is introduced at discrete points. In the case of distributed imbalances, it also leverages the Fourier extension method to generate semi-analytical results. NAT can be used to accurately balance linear rotor-bearing systems with stepped shafts and multiple disks that are supported on roller bearings [32] or fluid film bearings [33]. It has been shown that the results of state-of-the-art finite element method (FEM) models converge to the results of NAT when the number of nodes is increased [33]. Furthermore, NAT reduces the computational time by a factor of ten compared to FEM [31]. This version of NAT has been improved from the ones described in [32] as Timoshenko beam theory is used instead of the Rayleigh theory and external damping is considered. Additionally, material damping is taken into account by a fractional derivative damping model. For further details on this version of NAT, readers are referred to [34], where it was used for an influence coefficient method.

The space fixed coordinate system, denoted as $Oxyz$, is established in accordance with the setup illustrated in Figure 1. In this configuration, x and y are mutually orthogonal and lie in the transverse plane to z, which traverses the undeflected axis of the rotor within its bearings.

Within the framework of NAT, the rotor-bearing system is characterized by N stations connected by $M=N-1$ segments. These stations correspond to components such as disks, bearings, steps, and the ends of the rotor, whereas the segments correspond to the cylindrical elements interconnecting these stations. Bearings are modelled as spin-speed-independent spring and damper elements, although the spin-speed-dependent behaviour of, e.g., fluid film bearings can also be modelled [33]. Disks are characterized by their mass $m^{\left(i\right)}$, mass moments of inertia about the x- and y-axis ${\Theta }_t^{\left(i\right)}$, mass moment of inertia about the z-axis ${\Theta }_p^{\left(i\right)}$, amount of eccentricity ${ϵ}^{\left(i\right)}$, and angular position of eccentricity ${\beta }^{\left(i\right)}$. Segments are defined by the frequency-dependent complex Young’s modulus $E_{\mathsf{\ell }}^{★}$ and shear modulus $G_{\mathsf{\ell }}^{★}$, density ${\rho }_{\mathsf{\ell }}$, area of the cross-section $A_{\mathsf{\ell }}$, its diametric moment of area about the x- and y-axis $I_{\mathsf{\ell }}$, external damping coefficient $d_{e\mathsf{\ell }}$ and the shear correction factor $k_{S\mathsf{\ell }}$. $\Omega$ denotes the spin speed of the system and $\omega$ its whirling speed. When the system is only subjected to imbalance excitation, $\Omega$ and $\omega$ are equal. The state of the rotor is defined by the vector $\overrightarrow{x}$ consisting of the displacements $r_x$ and $r_y$, the rotations of the cross-section ${\phi }_x$ and ${\phi }_y$, the bending moments $M_x$ and $M_y$, and the shear forces $Q_x$ and $Q_y$. The $\widetilde{•}$ symbol denotes the Fourier transform of the value •. Due to the fact that the solutions are complex conjugated, only four of the eight resulting differential equations are necessary to describe the system. Therefore only the solutions with the + index are considered.

2.1.1. Characteristic Equations

The Timoshenko beam theory [35] is used together with a fractional derivative damping model developed by Pritz et al. [36] to obtain the characteristic equations of a segment ℓ

$$\begin{array}{c}\begin{array}{cc}\hfill & \frac{{\partial }^2{\tilde{r}}_{x\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}^2}+\left({\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}-{\overline{d}}_{e\mathsf{\ell }}\right){\tilde{r}}_{x\mathsf{\ell }}^{+}-\frac{\partial {\tilde{\phi }}_{y\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}}=-\frac{{\overline{\Omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}}{2}{\mathrm{e}}^{\mathrm{i}{\beta }_{\mathsf{\ell }}},\hfill \end{array}\hfill \end{array}$$

(1a)

$$\begin{array}{c}\begin{array}{cc}\hfill & \frac{{\partial }^2{\tilde{r}}_{y\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}^2}+\left({\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}-{\overline{d}}_{e\mathsf{\ell }}\right){\tilde{r}}_{y\mathsf{\ell }}^{+}+\frac{\partial {\tilde{\phi }}_{x\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}}=\frac{\mathrm{i}{\overline{\Omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}}{2}{\mathrm{e}}^{\mathrm{i}{\beta }_{\mathsf{\ell }}},\hfill \end{array}\hfill \end{array}$$

(1b)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{c}\begin{array}{cc}\hfill & \frac{{\partial }^2{\tilde{\phi }}_{y\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}^2}+\left({\overline{\omega }}_{\mathsf{\ell }}^2-\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}\right){\tilde{\phi }}_{y\mathsf{\ell }}^{+}-\mathrm{i}\frac{{\overline{E}}_{1\mathsf{\ell }}}{{\overline{E}}_{0\mathsf{\ell }}}\frac{{\partial }^2{\tilde{\phi }}_{x\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}^2}+\mathrm{i}{\overline{g}}_{\mathsf{\ell }}{\tilde{\phi }}_{x\mathsf{\ell }}^{+}+\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}\frac{\partial {\tilde{r}}_{x\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}}=0,\hfill \end{array}\hfill \end{array}$$

(1c)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{c}\begin{array}{cc}\hfill & \frac{{\partial }^2{\tilde{\phi }}_{x\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}^2}+\left({\overline{\omega }}_{\mathsf{\ell }}^2-\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}\right){\tilde{\phi }}_{x\mathsf{\ell }}^{+}+\mathrm{i}\frac{{\overline{E}}_{1\mathsf{\ell }}}{{\overline{E}}_{0\mathsf{\ell }}}\frac{{\partial }^2{\tilde{\phi }}_{y\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}^2}-\mathrm{i}{\overline{g}}_{\mathsf{\ell }}{\tilde{\phi }}_{y\mathsf{\ell }}^{+}-\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}\frac{\partial {\tilde{r}}_{y\mathsf{\ell }}^{+}}{\partial z_{\mathsf{\ell }}}=0,\hfill \end{array}\hfill \end{array}$$

(1d)

with

$$\begin{array}{cc}{\overline{\omega }}_{\mathsf{\ell }}^2\hfill & =\frac{{\omega }^2{\rho }_{\mathsf{\ell }}}{E_{0\mathsf{\ell }}^{★}},{\overline{\Omega }}_{\mathsf{\ell }}^2=\frac{{\Omega }^2{\rho }_{\mathsf{\ell }}}{E_{0\mathsf{\ell }}^{★}},{\overline{E}}_{0\mathsf{\ell }}=\frac{E_{0\mathsf{\ell }}^{★}}{G_{\mathsf{\ell }}^{★}{k}_{S_{\mathsf{\ell }}}},{\overline{E}}_{1\mathsf{\ell }}=\frac{E_{1\mathsf{\ell }}^{★}}{G_{\mathsf{\ell }}^{★}{k}_{S_{\mathsf{\ell }}}},\hfill \\ {\overline{d}}_{e\mathsf{\ell }}\hfill & =\frac{\mathrm{i}\omega d_{a\mathsf{\ell }}}{G_{\mathsf{\ell }}^{★}{k}_{S_{\mathsf{\ell }}}{A}_{\mathsf{\ell }}},r_{\mathsf{\ell }}^2=\frac{I_{\mathsf{\ell }}}{A_{\mathsf{\ell }}},{\overline{g}}_{\mathsf{\ell }}=\frac{2\omega {\rho }_{\mathsf{\ell }}\Omega }{E_{0\mathsf{\ell }}^{★}}.\hfill \\ E_{0\mathsf{\ell }}^{★}\hfill & =a_{0\mathsf{\ell }}^E+\frac{a_{1\mathsf{\ell }}^E}{2}[{\left(\mathrm{i}(\omega -\Omega )\right)}^{{\alpha }_{\mathsf{\ell }}^E}+{\left(\mathrm{i}(\omega +\Omega )\right)}^{{\alpha }_{\mathsf{\ell }}^E}],\hfill \\ E_{1\mathsf{\ell }}^{★}\hfill & =\frac{a_{1\mathsf{\ell }}^E}{2}[{\left(\mathrm{i}(\omega -\Omega )\right)}^{{\alpha }_{\mathsf{\ell }}^E}-{\left(\mathrm{i}(\omega +\Omega )\right)}^{{\alpha }_{\mathsf{\ell }}^E}],\hfill \\ G_{\mathsf{\ell }}^{★}\hfill & =\frac{a_{0\mathsf{\ell }}^E+{\left(\mathrm{i}\omega \right)}^{{\alpha }_{\mathsf{\ell }}^E}{a}_{1\mathsf{\ell }}^E}{2(1+{\nu }_{\mathsf{\ell }})}.\hfill \end{array}$$

where $a_{0\mathsf{\ell }}^E$, $a_{1\mathsf{\ell }}^E$ and ${\alpha }_{\mathsf{\ell }}^E$ are the parameters of the fractional derivative material damping model. These are derived by measurements and are available for many commonly used materials [37]. A detailed derivation of the characteristic equations is found in [34].

2.1.2. Homogeneous Solution

To find the homogeneous solution the imbalance excitation is set to zero and solutions of the form of

$$\begin{array}{cc}{\tilde{r}}_{hx\mathsf{\ell }}^{+}\hfill & =C_{r_x\mathsf{\ell }}{e}^{\mathrm{i}k{z}_{\mathsf{\ell }}},\hfill \end{array}$$

(2a)

$$\begin{array}{cc}{\tilde{r}}_{hy\mathsf{\ell }}^{+}\hfill & =C_{r_y\mathsf{\ell }}{e}^{\mathrm{i}k{z}_{\mathsf{\ell }}},\hfill \end{array}$$

(2b)

$$\begin{array}{cc}{\tilde{\phi }}_{xh\mathsf{\ell }}^{+}\hfill & =C_{{\phi }_x\mathsf{\ell }}{e}^{\mathrm{i}k{z}_{\mathsf{\ell }}},\hfill \end{array}$$

(2c)

$$\begin{array}{cc}{\tilde{\phi }}_{yh\mathsf{\ell }}^{+}\hfill & =C_{{\phi }_y\mathsf{\ell }}{e}^{\mathrm{i}k{z}_{\mathsf{\ell }}},\hfill \end{array}$$

(2d)

are assumed. This results in a system of equations

$$\begin{array}{c}\begin{array}{cc}\hfill & \left[\begin{array}{cccc}-k^2+{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}-{\overline{d}}_{e\mathsf{\ell }}& 0& -\mathrm{i}k& 0\\ 0& -k^2+{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}-{\overline{d}}_{e\mathsf{\ell }}& 0& \mathrm{i}k\\ \mathrm{i}k\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}& 0& -k^2+{\overline{\omega }}_{\mathsf{\ell }}^2-\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}& \mathrm{i}(k^2\frac{{\overline{E}}_{1\mathsf{\ell }}}{{\overline{E}}_{0\mathsf{\ell }}}+{\overline{g}}_{\mathsf{\ell }})\\ 0& -\mathrm{i}k\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}& -\mathrm{i}(k^2\frac{{\overline{E}}_{1\mathsf{\ell }}}{{\overline{E}}_{0\mathsf{\ell }}}+{\overline{g}}_{\mathsf{\ell }})& -k^2+{\overline{\omega }}_{\mathsf{\ell }}^2-\frac{1}{{\overline{E}}_{0\mathsf{\ell }}{r}_{\mathsf{\ell }}^2}\end{array}\right]\left[\begin{array}{c}{C}_{r_x\mathsf{\ell }}\\ C_{r_y\mathsf{\ell }}\\ C_{{\phi }_y\mathsf{\ell }}\\ C_{{\phi }_x\mathsf{\ell }}\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\hfill \end{array}\hfill \end{array}$$

(3)

with the non-trivial solutions of k

$$\begin{array}{cc}{k}_{1,2}\hfill & =\pm \sqrt{\frac{1}{2}\left(-{\overline{d}}_{e\mathsf{\ell }}-{\overline{g}}_{\mathsf{\ell }}^{+}+{\overline{\omega }}_{+\mathsf{\ell }}^2+{\overline{E}}_{0\mathsf{\ell }}{\overline{\omega }}_{\mathsf{\ell }}^2-R_{+\mathsf{\ell }}\right)}=\pm {\alpha }_{1\mathsf{\ell }}\hfill \end{array}$$

(4a)

$$\begin{array}{cc}{k}_{3,4}\hfill & =\pm \sqrt{\frac{1}{2}\left(-{\overline{d}}_{e\mathsf{\ell }}-{\overline{g}}_{\mathsf{\ell }}^{+}+{\overline{\omega }}_{+\mathsf{\ell }}^2+{\overline{E}}_{0\mathsf{\ell }}{\overline{\omega }}_{\mathsf{\ell }}^2+R_{+\mathsf{\ell }}\right)}=\pm {\alpha }_{2\mathsf{\ell }}\hfill \end{array}$$

(4b)

$$\begin{array}{cc}{k}_{5,6}\hfill & =\pm \sqrt{\frac{1}{2}\left(-{\overline{d}}_{e\mathsf{\ell }}+{\overline{g}}_{\mathsf{\ell }}^{-}+{\overline{\omega }}_{-\mathsf{\ell }}^2+{\overline{E}}_{0\mathsf{\ell }}{\overline{\omega }}_{\mathsf{\ell }}^2-R_{-\mathsf{\ell }}\right)}=\pm {\alpha }_{3\mathsf{\ell }}\hfill \end{array}$$

(4c)

$$\begin{array}{cc}{k}_{7,8}\hfill & =\pm \sqrt{\frac{1}{2}\left(-{\overline{d}}_{e\mathsf{\ell }}+{\overline{g}}_{\mathsf{\ell }}^{-}+{\overline{\omega }}_{-\mathsf{\ell }}^2+{\overline{E}}_{0\mathsf{\ell }}{\overline{\omega }}_{\mathsf{\ell }}^2+R_{-\mathsf{\ell }}\right)}=\pm {\alpha }_{4\mathsf{\ell }}\hfill \end{array}$$

(4d)

with

$$\begin{array}{c}\begin{array}{cc}{R}_{\pm \mathsf{\ell }}\hfill & =\sqrt{-\frac{4{\overline{d}}_{e\mathsf{\ell }}}{{\overline{E}}_{\pm \mathsf{\ell }}{r}_{\mathsf{\ell }}^2}+\frac{4{\overline{\omega }}_{\pm \mathsf{\ell }}^2}{r_{\mathsf{\ell }}^2}+{\left({\overline{d}}_{e\mathsf{\ell }}\mp {\overline{g}}_{\mathsf{\ell }}^{\pm }+{\overline{\omega }}_{\pm \mathsf{\ell }}^2-{\overline{E}}_{0\mathsf{\ell }}{\overline{\omega }}_{\mathsf{\ell }}^2\right)}^2},\hfill \end{array}\hfill \\ \begin{array}{cc}{\overline{\omega }}_{\pm \mathsf{\ell }}^2\hfill & =\frac{{\omega }^2{\rho }_{\mathsf{\ell }}}{E_{\pm \mathsf{\ell }}^{★}},{\overline{E}}_{\pm \mathsf{\ell }}=\frac{E_{\pm \mathsf{\ell }}^{★}}{G_{\mathsf{\ell }}^{★}{k}_{S\mathsf{\ell }}},\hfill \\ {\overline{g}}_{\pm \mathsf{\ell }}\hfill & =\frac{2\omega {\rho }_{\mathsf{\ell }}\Omega }{E_{\pm \mathsf{\ell }}^{★}},E_{\pm \mathsf{\ell }}^{★}=a_{0\mathsf{\ell }}^E+a_{1\mathsf{\ell }}^E{(\mathrm{i}(\omega \mp \Omega )}^{{\alpha }_{\mathsf{\ell }}^E}.\hfill \end{array}\hfill \end{array}$$

Inserting either $k_1,k_2,k_3$ or $k_4$ into Equation (3) leads to the relations

$$\begin{array}{cc}\hfill & C_{r_y\mathsf{\ell }}=-\mathrm{i}{C}_{r_x\mathsf{\ell }},\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill & C_{{\phi }_x\mathsf{\ell }}=\mathrm{i}{C}_{{\phi }_y\mathsf{\ell }},\hfill \end{array}$$

(5b)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill & C_{{\phi }_y\mathsf{\ell }}=\mathrm{i}\frac{k_{•}^2-{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}+{\overline{d}}_{e\mathsf{\ell }}}{k_{•}}{C}_{r_x\mathsf{\ell }},\hfill \end{array}$$

(5c)

while inserting either $k_5,k_6,k_7$ or $k_8$ yields

$$\begin{array}{cc}\hfill & C_{r_y\mathsf{\ell }}=\mathrm{i}{C}_{r_x\mathsf{\ell }},\hfill \end{array}$$

(6a)

$$\hspace{1em}\begin{array}{cc}\hfill & C_{{\phi }_x\mathsf{\ell }}=-\mathrm{i}{C}_{{\phi }_y\mathsf{\ell }},\hfill \end{array}$$

(6b)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill & C_{{\phi }_y\mathsf{\ell }}=\mathrm{i}\frac{k_{•}^2-{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}+{\overline{d}}_{e\mathsf{\ell }}}{k_{•}}{C}_{r_x\mathsf{\ell }}.\hfill \end{array}$$

(6c)

The state of each segment is defined by the matrix equation

$${\overrightarrow{\tilde{x}}}_{h\mathsf{\ell }}^{+}={\underline{\underline{B}}}_{\mathsf{\ell }}\overrightarrow{c_{\mathsf{\ell }}},$$

(7)

where ${\overrightarrow{\tilde{x}}}_{h\mathsf{\ell }}^{+}$ is the homogeneous solution vector, ${\underline{\underline{B}}}_{\mathsf{\ell }}$ is the state variable matrix, and $\overrightarrow{c_{\mathsf{\ell }}}$ is the vector of arbitrary constants, which are defined as

$${\overrightarrow{\tilde{x}}}_{h\mathsf{\ell }}^{+}={[{\tilde{r}}_{hx\mathsf{\ell }}^{+},{\tilde{r}}_{hy\mathsf{\ell }}^{+},{\tilde{\phi }}_{hx\mathsf{\ell }}^{+},{\tilde{\phi }}_{hy\mathsf{\ell }}^{+},{\tilde{M}}_{hx\mathsf{\ell }}^{+},{\tilde{M}}_{hy\mathsf{\ell }}^{+},{\tilde{Q}}_{hx\mathsf{\ell }}^{+},{\tilde{Q}}_{hy\mathsf{\ell }}^{+}]}^T,$$

$${\overrightarrow{c}}_{\mathsf{\ell }}^{+}={[c_{1\mathsf{\ell }},c_{2\mathsf{\ell }},c_{3\mathsf{\ell }},c_{4\mathsf{\ell }},c_{5\mathsf{\ell }},c_{6\mathsf{\ell }},c_{7\mathsf{\ell }},c_{8\mathsf{\ell }}]}^T$$

and

$$\begin{array}{cc}\hfill & {\underline{\underline{B}}}_{\mathsf{\ell }}^T=\hfill \\ \hfill & \left[\begin{array}{cccccccc}{E}_1^{+}& E_1^{-}& E_2^{+}& E_2^{-}& E_3^{+}& E_3^{-}& E_4^{+}& E_4^{-}\\ -\mathrm{i}{E}_1^{+}& -\mathrm{i}{E}_1^{-}& -\mathrm{i}{E}_2^{+}& -\mathrm{i}{E}_2^{-}& \mathrm{i}{E}_3^{+}& \mathrm{i}{E}_3^{-}& \mathrm{i}{E}_4^{+}& \mathrm{i}{E}_4^{-}\\ k_{\phi 1}{E}_1^{+}& -k_{\phi 1}{E}_1^{-}& k_{\phi 2}{E}_2^{+}& -k_{\phi 2}{E}_2^{-}& k_{\phi 3}{E}_3^{+}& -k_{\phi 3}{E}_3^{-}& k_{\phi 4}{E}_4^{+}& -k_{\phi 4}{E}_4^{-}\\ \mathrm{i}{k}_{\phi 1}{E}_1^{+}& -\mathrm{i}{k}_{\phi 1}{E}_1^{-}& \mathrm{i}{k}_{\phi 2}{E}_2^{+}& -\mathrm{i}{k}_{\phi 2}{E}_2^{-}& -\mathrm{i}{k}_{\phi 3}{E}_3^{+}& \mathrm{i}{k}_{\phi 3}{E}_3^{-}& -\mathrm{i}{k}_{\phi 4}{E}_4^{+}& \mathrm{i}{k}_{\phi 4}{E}_4^{-}\\ -k_{M1}^{+}{E}_1^{+}& -k_{M1}^{+}{E}_1^{-}& -k_{M2}^{+}{E}_2^{+}& -k_{M2}^{+}{E}_2^{-}& -k_{M3}^{-}{E}_3^{+}& -k_{M3}^{-}{E}_3^{-}& -k_{M4}^{-}{E}_4^{+}& -k_{M4}^{-}{E}_4^{-}\\ -\mathrm{i}{k}_{M1}^{+}{E}_1^{+}& -\mathrm{i}{k}_{M1}^{+}{E}_1^{-}& -\mathrm{i}{k}_{M2}^{+}{E}_2^{+}& -\mathrm{i}{k}_{M2}^{+}{E}_2^{-}& \mathrm{i}{k}_{M3}^{-}{E}_3^{+}& \mathrm{i}{k}_{M3}^{-}{E}_3^{-}& \mathrm{i}{k}_{M4}^{-}{E}_4^{+}& \mathrm{i}{k}_{M4}^{-}{E}_4^{-}\\ -k_{Q1}{E}_1^{+}& k_{Q1}{E}_1^{-}& -k_{Q2}{E}_2^{+}& k_{Q2}{E}_2^{-}& -k_{Q3}{E}_3^{+}& k_{Q3}{E}_3^{-}& -k_{Q4}{E}_4^{+}& k_{Q4}{E}_4^{-}\\ \mathrm{i}{k}_{Q1}{E}_1^{+}& -\mathrm{i}{k}_{Q1}{E}_1^{-}& \mathrm{i}{k}_{Q2}{E}_2^{+}& -\mathrm{i}{k}_{Q2}{E}_2^{-}& -\mathrm{i}{k}_{Q3}{E}_3^{+}& \mathrm{i}{k}_{Q3}{E}_3^{-}& -\mathrm{i}{k}_{Q4}{E}_4^{+}& \mathrm{i}{k}_{Q4}{E}_4^{-}\end{array}\right]\hfill \end{array}$$

with

$$\begin{array}{cc}{E}_{•}^{+}\hfill & =e^{\mathrm{i}{\alpha }_{•\mathsf{\ell }}{z}_{\mathsf{\ell }}},E_{•}^{-}=e^{\mathrm{i}{\alpha }_{•\mathsf{\ell }}(L_{\mathsf{\ell }}-z_{\mathsf{\ell }})},k_{\phi •}=\mathrm{i}\frac{{\alpha }_{•\mathsf{\ell }}^2-{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}+{\overline{d}}_{e\mathsf{\ell }}}{{\alpha }_{•\mathsf{\ell }}},\hfill \\ k_{M•}^{\pm }\hfill & =E_{\pm \mathsf{\ell }}^{★}{I}_{\mathsf{\ell }}\left({\alpha }_{•\mathsf{\ell }}^2-{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}+{\overline{d}}_{e\mathsf{\ell }}\right),k_{Q•}=\mathrm{i}{G}_{\mathsf{\ell }}^{★}{k}_{S_{\mathsf{\ell }}}{A}_{\mathsf{\ell }}\left(\frac{-{\overline{\omega }}_{\mathsf{\ell }}^2{\overline{E}}_{0\mathsf{\ell }}+{\overline{d}}_{e\mathsf{\ell }}}{{\alpha }_{•\mathsf{\ell }}}\right).\hfill \end{array}$$

2.1.3. Boundary and Interface Conditions

To yield a unique solution, the governing equations require boundary conditions on the left and right end of the rotor and interface conditions on every station within the rotor. Since solving Equation (1) leads to fourth-order equations, eight boundary or interface conditions are necessary for each segment. The process in which the boundary- and interface conditions are found for different bearing types is described in detail in [33].

2.1.4. Particular Solution

The particular solutions for concentrated imbalances and arbitrarily distributed imbalances are derived utilizing the Fourier extension method [38], Greens function method [39], the residue theorem and Jordan’s lemma [40]. The complete derivation is omitted for the sake of brevity.

2.1.5. Assembly and Solution Procedure

The total solution vector of each rotor segment consists of the homogeneous and the particular solution ${\overrightarrow{x}}_{\mathsf{\ell }}={\overrightarrow{x}}_{h\mathsf{\ell }}+{\overrightarrow{x}}_{p\mathsf{\ell }}$. A system of linear equations is found by assembling the characteristic equations of each segment with the boundary and interface conditions

$$\underline{\underline{A}}\overrightarrow{c}=\overrightarrow{b},$$

(8)

where $\underline{\underline{A}}$ is the system matrix, $\overrightarrow{b}$ the right-hand side vector and $\overrightarrow{c}$ a vector of unknown constants. The unknown constants $\overrightarrow{c}$ can be transformed into the total solution vector $\overrightarrow{x}$, which defines the behaviour of the rotor, using Equation (7). A detailed explanation of the assembly procedure of NAT is found in [31].

2.2. Modal Balancing Method

The proposed method builds on the N plane method, which corrects modal components of unbalance in a progressive way near critical speeds without affecting already balanced modes [1]. To avoid confusion with the index N for the rotor stations, the index P is used for the number of balanced modes in this paper. Modal balancing using NAT was first presented in [32]. The proposed version adds a distinction between balancing speeds and critical speeds, to make it better suited to balance systems under real-world conditions, where critical speeds often can not be reached safely before balancing.

The proposed model-based balancing method is built on three matrices: the measured vibrations of the initial system at the balancing speeds ${\underline{\underline{r}}}_I$, the calculated vibrations of the system under unit unbalance $U_U$ at the balancing speeds ${\underline{\underline{r}}}_S$ and the calculated vibrations of the system under unit unbalance at the critical speeds ${\underline{\underline{r}}}_{crit}$. Since ${\underline{\underline{r}}}_S$ and ${\underline{\underline{r}}}_{crit}$ are simulated using NAT, trial runs are avoided and the rotor-bearing system is balanced with a single measurement run of the initial system without additional weights. With this approach, the assembly needs to be dismantled only once to mount the balancing weights, as long as suitable measurement equipment is in place.

Orthogonal balancing weights that only influence a single mode shape are found by multiplying generalised unbalances with a modal matrix $\underline{\underline{\varphi }}$. This matrix consists of the relations between rotor displacements at each critical speed and measurement position. The modal matrix $\underline{\underline{\varphi }}$ is found by making the simulated matrix ${\underline{\underline{r}}}_{crit}$ non-dimensional. This procedure allows us to determine accurate mode shapes, without measuring the unbalanced system at critical speeds, which is often infeasible. A balanced rotor is excited by the unbalance mounted for compensation ${\overrightarrow{U}}_C$ and its initial unbalance ${\overrightarrow{U}}_I$. To achieve perfect balance, both excitations need to cancel each other out and all displacements must vanish

$${\underline{\underline{\varphi }}}^T{\overrightarrow{U}}_C=-{\underline{\underline{\varphi }}}^T{\overrightarrow{U}}_I.$$

(9)

Through modal decomposition, Equation (9) can be divided into one equation for each critical speed

$${\overrightarrow{\varphi }}_p^T{\overrightarrow{U}}_C=-{\overrightarrow{\varphi }}_p^T{\overrightarrow{U}}_{I}p=1\dots P.$$

(10)

By comparing the measured excitation caused by the unknown initial unbalance to the excitation caused by a unit unbalance calculated with NAT at the same position and spin speed, a general unbalance $U_{gen}$ can be calculated, which quantifies how much each mode shape is excited by the initial unbalance

$$U_{gen}^p=\frac{1}{Q}{\sum }_{q=1}^Q\frac{r_I^{pq}}{r_S^{pq}}{U}_U$$

(11)

where p is the index of the critical speed, q is the index of the measurement position and Q is the number of measurement positions. Since there are usually multiple measurement positions available, the surplus information is used to optimise the calculation of the generalised unbalance $U_{gen}$ using the least squares method. Simply mounting the generalised unbalance at the opposite side of the rotor would minimize the vibration of the according mode but at the cost of influencing the vibration of all other modes. To avoid this problem, the generalised unbalances are modified with the modal matrix to generate orthogonal balancing weights

$${\overrightarrow{U}}_C^p=-{\left({\overrightarrow{\varphi }}_p^T\right)}^{-1}\left({\overrightarrow{\varphi }}_p^T{\overrightarrow{U}}_{gen}\right).$$

(12)

A set of balancing weights is called orthogonal if it balances a single mode without influencing all other considered modes. In theory, if a set of orthogonal balancing weights is mounted for each mode of the system—which means $P\to \infty$ for a continuous system—the system would be perfectly balanced and thus

$$\overrightarrow{U_I}=-{\sum }_{p=1}^P{\overrightarrow{U_C}}^p.$$

(13)

In practice, it is sufficient to only balance all modes up to the balancing speed and sometimes one or two modes above the critical speed to reach acceptable vibration levels [1].

2.3. Theory of Bow Compensation

There are three classical ways of bow compensation. Nicholas et al. [4] first described these three methods for the Laval rotor. Rao [41], Deepthikumar et al. [12,42], Sanches and Pederiva [43] among others expanded upon these approaches to identify bow and unbalance and balance various rotor-bearing systems. The total excitation of a warped rotor consists of the unbalance excitation and the excitation due to the pre-bend of the shaft. The resulting displacement amplitude of the rotor is referred to as the total whirl, while the displacement caused only by the unbalance is called the net whirl.

Method 1 aims to reduce the total whirl at the balancing speed to zero. The correction unbalance is therefore

$$U_C=-\frac{r_I}{r_S}{U}_U.$$

(14)

This method corresponds to the balancing of a straight rotor and usually takes place when the operator is not aware of the bow of the shaft. Unfortunately, this method leads to the worst results, since the balancing weights not only counteract the unbalance but also the effect of the bow at the balancing speed. This “overbalancing” leads to residual unbalance at every other speed, except the balancing speed, and is especially harmful near the critical speed. The only way to consistently reduce unbalance over the whole frequency range with this approach is to balance at the critical speed, which is in practice often not feasible.

The second method only balances the net whirl and takes the effect of the warped shaft out of the balancing process. This is achieved by deducting the slow roll bow from the vibration measurements at the balancing speed

$$U_C=-\frac{r_I-r_0}{r_S}{U}_U$$

(15)

where $r_0$ is the complex residual bow. This term describes the vector from the journal axis to the geometrical centre of the shaft at a standstill in a specific plane in complex coordinates. For a completely straight rotor, the complex residual bow is 0 at every axial position. The interaction between residual bow $r_0$, displacement due to unbalance excitation $r_u$ and total displacement $r_I$ is demonstrated in Figure 2.

The second balancing procedure leads to a very good overall reduction in vibration. Since the vibration amplitude due to the warped shaft is not constant and only approximately corresponds to the slow-roll bow, a certain displacement at the critical speed remains, even with faultless balancing.

The third approach—which is recommended by Nicholas et al. [4]—aims to minimize the total whirl at the critical speed, by measuring the subcritical balancing speed. This is performed by modifying the vibration due to the shaft bow with a factor

$$U_C=-\frac{r_I}{r_S}{U}_U+\left(\frac{{\Omega }_{cr}^2-{\Omega }_b^2}{{\Omega }_b^2}\right)\frac{r_0}{r_S}{U}_U$$

(16)

where ${\Omega }_b$ is the balancing speed and ${\Omega }_{cr}$ is the first critical speed. In theory, this leads to the best reduction in vibration for a warped Laval rotor, where the vibration amplitude never exceeds the slow-roll bow and is even minimized at the first critical speed.

The ideal results of the three balancing methods for warped Laval rotors are summarized in Figure 3.

Which of these approaches is optimal for the presented model-based modal balancing method using NAT on a multi-disk rotor-bearing system is investigated in Section 3.4.

3. Results and Discussion

This section provides a detailed description of the rotor-dynamic test bed, its corresponding NAT model, and the experiments conducted. First, the proposed method is demonstrated for a numerical example. After that, the optimal method of bow compensation for the model-based balancing approach is investigated on a one-disk rotor. Then, the proposed balancing technique is applied to balance the first mode and later, the first two modes of the test setup. Furthermore, the eigenvalues and mode shapes are calculated and compared to measurements.

3.1. Test Bed

Figure 4 shows the rotor dynamic test bed of the Institute of Mechanics of the Graz University of Technology.

The elastic shaft, with a maximum diameter of 15 mm, can be mounted with any number of disks and is supported by two roller bearings. Due to the considerably greater stiffness of the supports compared to the isotropic shaft, the anisotropy of the bearing housing is neglected during the balancing process. The shaft is connected to a Kemmerich HM90D-140 P electric motor through a magnetic coupling, with a maximum spin speed of 400 Hz and a power of 6 kW. Due to safety requirements, the whole setup is limited to spin speeds up to 150 Hz, which makes measurements above the second eigenfrequency infeasible. The unbalance response is measured by four Keyence LK-H055 laser displacement measurement systems, with an accuracy of 0.025 $\mathsf{\mu }\mathrm{m}$ and the phase is monitored by an IE5349 inductive sensor manufactured by IFM electronics. The test bed is safeguarded by two safety bearings and blast protection during operation. Five disks of varying weights were mounted onto the shaft, as seen in Figure 5. To maximize the effect of material damping on the rotor-bearing system, external damping was minimized by using roller bearings and no external damping elements.

The initial rotor is excited by the unbalance caused by manufacturing errors and the bow of the shaft.

Table 1. Bow of the pre-bent shaft.

	Left Disk	Right Disk
Axial position	88 mm	448 mm
Bow	0.08694 mm	0.11108 mm
Angular position	${173.91}^{\circ }$	${196.53}^{\circ }$

lists the slow-roll bow of the shaft.

Additionally, two unbalance weights are mounted on the disks: 5.57 g at 42.5 mm radius at the 0${}^{\circ }$ position of the first balancing disk at $z=0.088$ m and the 180${}^{\circ }$ position of the second balancing disk at $z=0.448$ m. These unbalances excite especially the second mode of the system.

3.2. Rotor Model

First, an NAT simulation of the test bed is generated. The parameters of the stations of this mathematical model are shown in

Table 2. Stations of the NAT model.

z	m	${\mathbf{\Theta }}_{\mathit{xy}}$	${\mathbf{\Theta }}_{\mathit{zz}}$	${\mathit{c}}_{\mathit{x}}$	${\mathit{c}}_{\mathit{y}}$	d
m	kg	kg m2	kg m2	N/m	N/m	Ns/m
0	0	0	0	0	0	0
0.008	0	0	0	282,540,000	129,968,400	80
0.012	0	0	0	0	0	0
0.088	2.0200	0.00150	0.0028	0	0	0
0.143	0.1300	$2.139\times {10}^{-5}$	$3.4589\times {10}^{-5}$	0	0	0
0.268	0.0926	$1.297\times {10}^{-5}$	$1.9986\times {10}^{-5}$	0	0	0
0.448	2.0200	0.00150	0.0028	0	0	0
0.524	0	0	0	0	0	0
0.528	0	0	0	282,540,000	129,968,400	80
0.574	0.6471	0.0020322	0.00034134	0	0	0
0.596	0	0	0	0	0	0

, where z is the axial position, m the mass of the disks, ${\Theta }_{xy}$ the angular mass of the disk about the x- and y-axis, ${\Theta }_{zz}$ the angular mass about the z-axis, $c_x$ the bearing stiffness in x-direction, $c_y$ the bearing stiffness in y-direction and d the damping coefficient of the bearing in both directions.

The cylindrical segments of the rotor are represented in the model with the following properties:

• Density: 7700 $\mathrm{kg}/{\mathrm{m}}^3$;
• Shear modulus: $8.1\times {10}^{11}\mathrm{N}/{\mathrm{m}}^2$;
• Shear correction factor: 0.89;
• External damping coefficient: 45 $\mathrm{Ns}/\mathrm{m}$.

The fractional derivative damping parameters are as follows:

• $a_0^E$: $2.2\times {10}^{11}$;
• $a_1^E$: $6.887\times {10}^9$;
• ${\alpha }^E$: $0.3$.

The Zener material model is also capable of modelling the very low viscoelastic properties of steel [37].

3.3. Numerical Analysis

Before the experimental investigation is carried out, a numerical analysis of the proposed balancing method is presented to show how the procedure works in principle. All calculations for the experiments were conducted on a 3.2 GHz Intel^® Core^TM i7-8700 CPU running Windows 10 using MATLAB^TM R2019a. For the numeric analysis, the rotor bow is neglected. The eigenvalues and Campbell diagram are found with a recursive search algorithm described in [32,33]. As is shown in Figure 6, the first two forward bending modes are calculated to be at 59.107 Hz and 145.48 Hz. These values will later be compared to measurements.

For the numeric balancing test, a concentrated unbalance at the disk and a distributed eccentricity of the shaft are assumed. The amount of eccentricity of the shaft, shown in Figure 7, is arbitrarily chosen follow to the function

$$ϵ=6.6\times {10}^{-4}z-5.4\frac{1}{\mathrm{m}}\times {10}^{-4}{z}^2+1\mathrm{m}\times {10}^{-4}log(\frac{z}{\mathrm{m}}+1)+8\mathrm{m}\times {10}^{-5}sin\left(2\pi \frac{z}{\mathrm{m}}\right).$$

(17)

Its direction changes in a corkscrew manner according to

$$\beta =2.5\pi z$$

(18)

and the eccentricities of the disks are listed in

Table 3. Eccentricity of the disks of the numeric example.

Axial Position	Amount	Direction
0.088 m	3.7233 $\times {10}^{-4}$ m	4.6251
0.143 m	1.600 $\times {10}^{-3}$ m	3.1400
0.268 m	1.600 $\times {10}^{-3}$ m	4.800
0.448 m	3.7233 $\times {10}^{-4}$ m	5.0789

.

The balancing planes are at position $z=0.088$ m and $z=0.448$ m. 50 Hz and 128 Hz are chosen as balancing speeds, which are in the information zone, but not so close to the critical speeds that measurements of an unbalanced system would be unrealistic. The displacements of the unbalanced example are compared to displacements of a system with unit unbalance, mode shapes are calculated and the correction unbalances, shown in

Table 4. Balancing weights of the numeric example.

	Left Disk	Right Disk
Amount	948 g mm	1297 g mm
Position	${249.54}^{\circ }$	${283.97}^{\circ }$

, are found within 0.17 s.

The unbalance responses of the balanced and unbalanced system are shown in Figure 8. The amplitude of vibration is measured at the three disks at positions $z=$ 0.088 m, 0.268 m and 0.448 m. In this paper, the comparisons of vibrational amplitude show always the average amplitude of these three measurement positions. The average reduction in vibration across all balancing planes and measured speeds is 93.43%.

In the next steps, the balancing procedure is replicated with an experimental setup. The aim of the experimental investigation is to find out how well the model-based balancing procedure works under the additional challenges of real-world systems like:

• Unknown parameters, especially the influence of the foundation and the coupling.
• Warped shafts.
• Measurement and mounting errors.

3.4. Investigation of Optimal Bow Compensation

Before the testing of the modal balancing procedure on the two disk rotor, the optimal way of compensating for the rotor bow with the model-based balancing method is investigated. The three classical compensation methods by Nicholas et al. [4] are considered. Method 1 minimizes the total displacement at the balancing speed, method 2 minimizes only the unbalance vibration while ignoring the residual bow and method 3 minimizes the total vibration at the first critical speed through measurements at the balancing speed as described in Section 2.3. For this investigation, the test bed described in Section 3.1 is modified so that only a single 2 kg disk is mounted on the middle of the shaft at $z=0.268$ to make the results comparable to the literature. The balancing weights, mounted at a radius of 42.5 mm, for all three approaches, are shown in

Table 5. Balancing weights of the one disk rotor.

	Method 1	Method 2	Method 3
Amount	8.71 g	5.80 g	7.63 g
Position	${5.720}^{\circ }$	${16.97}^{\circ }$	${9.224}^{\circ }$

.

The results of mounting these balancing weights are shown in Figure 9. The measured eigenfrequency matches the eigenfrequency calculated with an NAT of 43.96 Hz very well.

Method 1 (and the unbalanced case) did not allow for a safe passage through the first critical speed, which is in line with [4], where the authors advise against the use of the method. Although Nicholas et al. recommend method 3 for single disk systems, method 2 is chosen for the following experiments. As can be seen in Figure 9, small modelling and mounting errors can reduce the efficiency of method 3 significantly. Since the balancing weights calculated with NAT are slightly too big, the minimum total vibration is at 42 Hz instead of 44 Hz. This leads to a rapid increase in vibration amplitude towards the resonance frequency comparable to method 1. Method 3 also only minimizes the total vibration at a specific spin speed and plane at the rotor, which is less useful if multiple disks and modes are balanced, unlike method 2, which circumvents the problem [42].

Therefore, the bow of the shaft is considered as follows for the experimental investigation: the vector of the slow roll bow is deducted from the subcritical vibration measurements according to method 2. The calculated unbalance response to a unit unbalance using NAT remains unchanged. This way, the balancing is performed corresponding to the net whirl instead of the total whirl. At supercritical speeds, the effect of unbalance dominates the influence of the shaft bow [44] and, therefore, the total whirl is used for balancing.

3.5. Balancing of the First Mode of the Multi-Disk System

In the first experimental test of the multi-disk rotor, only the first mode of the system is balanced. The left disk at $z=$0.088 m is chosen as the balancing plane. The unbalanced system cannot reach the first critical speed safely; therefore, 50 Hz is chosen as the balancing speed. As described before, the complex residual bow is deducted from the measured total whirl at 50 Hz to obtain the net whirl. The net whirl is compared to the response to a unit unbalance calculated with NAT at 50 Hz. This leads to a balancing weight of 9.5458 g at a radius of 42.5 mm and the ${340.77}^{\circ }$ position within a total calculation time of 0.071 s. Mounting this balancing weight reduces the vibration of the first mode significantly as is shown in Figure 10.

The average amplitude of vibration is reduced by 61.39% at 55 Hz. Since only the first mode is balanced, the amplitude of the second mode is not reduced but even increased. Furthermore, it is observable that the calculated eigenfrequency of 59.107 Hz matches the measured frequency of the highest amplitude well. After balancing, the first critical speed can be transversed safely and the calculated and measured first mode shape is compared. The calculated dimensionless mode shape matches the measurement accurately, as is shown in Figure 11.

3.6. Balancing of the First and Second Mode of the Multi-Disk System

Even though the balancing weight of the previous example reduces the amplitude of the first mode significantly, it does increase the excitation of the second mode. For this reason, the first two modes are balanced simultaneously in the second example. The slow roll bow and the total whirl of the test bed at 50 Hz and 128 Hz are measured. Comparing these values with values calculated with NAT following the presented model-based modal balancing leads to the balancing weights shown in

Table 6. Balancing weights of the second experimental example.

	Left Disk	Right Disk
Amount	1.4808 g	9.0608 g
Position	${12.136}^{\circ }$	${359.82}^{\circ }$

. The total calculation time is 0.1042 s.

The calculated balancing weights reduce the vibration of the first two flexural modes of the system significantly, as is shown in Figure 12.

The average amplitude of vibration is reduced by 71.95% at 55 Hz and 45.96% at 128 Hz. The calculated and measured dimensionless second mode shape is compared after balancing, as is shown in Figure 13. As for the first mode shape, the calculation matches the experimental data well.

3.7. Discussion

The experiments show that the proposed simulation and balancing method works for the investigated examples and opens up interesting directions for further research:

• Despite the unknown foundation parameters, the eigenfrequencies were accurately predicted. In this test case, the shaft was substantially more flexible than the bearings and bearing foundations, reducing their influence. Further research could extend NAT with detailed foundation behaviour and investigate systems on elastic supports. Also, the balancing on rotors supported on fluid film bearings using NAT has only been theoretically presented [33], but has not yet been experimentally verified.
• The fractional time derivative damping model depicted the behaviour of the shaft correctly. Since the shaft is made out of steel, which shows very low viscoelastic properties, a less sophisticated material model would suffice. The applied material model in combination with NAT has already been used to calculate eigenfrequencies of strongly viscoelastic materials like Polyvinylchlorid [37], suggesting that the balancing of rotors made out of these materials could also be possible. This should be verified experimentally.
• Although NAT assumes a straight shaft, the vibration amplitudes of the first two modes of a warped rotor were successfully reduced using the presented procedure. Since the pre-bend of the test bed is significantly higher than the usually accepted values for high-speed machinery, it is inferred that a slight bow of the shaft poses no problem for the method.

4. Conclusions

The proposed model-based balancing procedure was numerically and experimentally validated to be effective for an isotropic rotor-bearing system with a warped shaft and low damping. The method is first demonstrated on a numerical example and then verified on a rotor test bed, where the first two modes are balanced without trial runs. The simulation shows high numerical efficiency as the balancing weights of all test cases are calculated in less than a second. The numerical test was conducted under idealized conditions and led to a total reduction in vibration of 93.43%. The experimental tests included common challenges for model-based balancing methods, like unknown foundation parameters and warped shafts, and led to a reduction in vibration amplitude between 45.96% and 71.95%, depending on the spin speed and test setup. Calculated mode shapes and eigenfrequencies match the measured values well, indicating that the proposed simulation method is sufficiently accurate.