Balancing of Flexible Rotors Supported on Fluid Film Bearings by Means of Influence Coefficients Calculated by the Numerical Assembly Technique

Abstract

In this paper, a new method for the balancing of rotor-bearing systems supported on fluid film bearings is proposed. The influence coefficients necessary for balancing are calculated using a novel simulation method called the Numerical Assembly Technique. The advantages of this approach are quasi-analytical solutions for the equations of motion of complex rotor-bearing systems and very low computation times. The Numerical Assembly Technique is extended by speed-dependent stiffness and damping coefficients approximated by the short-bearing theory to model the behavior of rotor systems supported on fluid film bearings. The rotating circular shaft is modeled according to the Rayleigh beam theory. The Numerical Assembly Technique is used to calculate the steady-state harmonic response, influence coefficients, eigenvalues, and the Campbell diagram of the rotor. These values are compared to simulations with the Finite Element Method to show the accuracy of the procedure. Two numerical examples of rotor-bearing systems are successfully balanced by the proposed balancing method.

1. Introduction

Rotary machines are always subject to imbalances, leading to an excitation of the rotor and the surrounding structure. This excitation causes audible noises and reduces the lifespan of the system. At spin speeds above 70% of the first flexural critical speed, the displacement of the shaft cannot be neglected, and the system is considered flexible [1]. In this case, flexible balancing is used to minimize the elastic displacement of the rotor throughout the whole length of the shaft and all operating speeds. There are two main groups of methods for the balancing of flexible rotors: modal balancing methods and influence coefficient methods. Combined methods are also frequently applied [2]. The modal balancing method requires only a small set of test runs and is accurate up to higher modes [3]. However, it can only be properly applied to rotor-bearing systems with roller bearings, where modal analysis is possible [4]. The influence coefficient method, on the other hand, requires only the assumption of linear relations between the unbalance forces and the vibration responses [5]. Usually, many test runs are necessary to determine the influence coefficients, which may be very time consuming and expensive. This can be avoided with calculated influence coefficients, as long as a very accurate model for the rotor-bearing system is used, including all rotor dynamic effects [6]. In this paper, the Numerical Assembly Technique (NAT) is used to generate a suitable model.

NAT is an efficient, quasianalytic method to calculate the steady-state harmonic response, influence coefficients, eigenvalues, and the Campbell diagram of rotor-bearing systems. NAT was introduced by Wu and Chou [7]. In its original form, it is applicable to harmonic vibrations of one-dimensional structures. In subsequent years, the method has been extended to additional applications in mechanics, especially in the field of rotor dynamics. Wu and Chen extended the method to the Timoshenko beam theory [8] and analyzed the natural frequencies of various nonuniform beams using continuum models [9,10,11]. The effects of multiple lumped masses, multiple-pinned supports, and rotary inertias on the calculation of Euler–Bernoulli beams were discussed in [12,13,14]. In [15], these approaches were applied to the Timoshenko beam theory, and the effect of the slenderness ratio and the shear coefficient were investigated. Lin combined various concentrated elements including point masses, rotary inertias, linear springs, rotational springs, and spring-mass systems for the calculation of Euler–Bernoulli and Timoshenko beams [16,17]. The effects of axial forces on Timoshenko [18,19] and Reddy-Bickford beams [20] were studied by Yesilce. In [21], NAT was used to calculate forward and backward whirling speeds of rotor systems, and the gyroscopic effects were analyzed. The method was extended to include Euler–Bernoulli beams with variable geometry or material discontinuities in [22]. Timoshenko beams with elastic supports were examined in [23]. In recent years, Klanner et al. introduced distributed loading [24], unbalance [25], and fractional derivative damping [26] in the NAT scheme. The advantages of NAT are that it leads to quasi-analytical solutions, and it is computationally more efficient than the Finite Element Method (FEM). In [25], numerical comparisons showed a reduction in computational time by a factor of ten compared to FEM. This research builds upon the previous work of Quinz et al. [27], where a balancing method using NAT was proposed, which was limited to isotropic rotor-bearing systems.

The main aim of this work is to find an analytic method that allows for a significant reduction in the vibration of rotor-bearing systems supported on fluid bearings, without the need for expensive and time-consuming test runs, and to prove its accuracy and numerical efficiency. The novelty of the introduced research is the extension of NAT to perform the balancing of rotor-bearing systems supported on fluid film bearings, the implementation of an improved search algorithm to calculate Campbell diagrams with NAT, and the validation of these techniques with numerical tests. The original in-house software was created based on the presented methodology.

The remainder of this paper is organized as follows: Section 2 describes the calculation and the balancing method. Section 3 presents the numeric tests and the obtained results. Section 4 discusses the advantages of the proposed method and other state-of-the-art techniques. Finally, Section 5 concludes the paper.

2. Material and Methods

In Figure 1, the general rotor problem is shown.

The space fixed coordinate system Oxyz is chosen so that Oz passes the undeflected axis of the rotor in its bearings. Ox and Oy are perpendicular to each other and transverse to Oz. The rotor is supported on fluid film bearings modeled as anisotropic springs and dampers including cross-coupling terms and rotates at a constant spin speed $\Omega$. Several disks are mounted on the shaft. Each disk has a mass $m^{\left(i\right)}$, a mass moment of inertia about the x- and y-axis ${\Theta }_t^{\left(i\right)}$, a mass moment of inertia about the z-axis ${\Theta }_p^{\left(i\right)}$, an amount of eccentricity ${\epsilon }^{\left(i\right)}$, and an angular position of eccentricity ${\beta }^{\left(i\right)}$. The rotor may also have a distributed unbalance defined by $\epsilon \left(z\right)$ and $\beta \left(z\right)$.

The simulation of the rotor-bearing system is carried out with NAT. The method gathers all the parameters necessary to describe the system in two arrays: one for the stations of the rotor bearings system and the other one for its segments. Stations represent disks, bearings, steps, and the ends of the rotor. Segments describe the cylindrical elements between the stations.

The rotor segments are modeled according to the Rayleigh beam theory. This theory is based on the same assumption as the Euler–Bernoulli beam theory but includes rotary inertias and gyroscopic effects [28]. In the following, the formulation of Klanner et al. [25] is extended to include the behavior of fluid film bearings.

The method is quasi-analytical in the sense that the governing equations are fulfilled exactly. Minor errors are introduced only at the interface conditions due to double-precision floating-point arithmetic when solving the system of linear equations. The application of the Fourier extension method in the case of distributed unbalances results also in small numerical errors. This is not the case for the proposed balancing technique, since the influence coefficients are calculated as reaction to a point force, which is calculated analytically.

2.1. Equations of Motion

The state of the rotor is represented by the displacements $u_x$ and $u_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 equilibrium conditions lead to the equations of motion for each segment (index ℓ)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & E_{\ell }{I}_{\ell }\frac{{\partial }^4{u}_{x\ell }(z_{\ell },t)}{\partial z_{\ell }^4}-{\rho }_{\ell }{I}_{\ell }\frac{{\partial }^4{u}_{x\ell }(z_{\ell },t)}{\partial z_{\ell }^2\partial t^2}+{\rho }_{\ell }{A}_{\ell }\frac{{\partial }^2{u}_{x\ell }(z_{\ell },t)}{\partial t^2}-2{\rho }_{\ell }{I}_{\ell }\Omega \frac{{\partial }^3{u}_{y\ell }(z_{\ell },t)}{\partial t\partial z_{\ell }^2}=\hfill \\ \hfill & q_x(z_{\ell },t)+{\rho }_{\ell }{A}_{\ell }\epsilon \left(z_{\ell }\right){\Omega }^{2}cos(\Omega t+\beta \left(z_{\ell }\right)),\hfill \end{array}\end{array}$$

(1a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & E_{\ell }{I}_{\ell }\frac{{\partial }^4{u}_{y\ell }(z_{\ell },t)}{\partial z_{\ell }^4}-{\rho }_{\ell }{I}_{\ell }\frac{{\partial }^4{u}_{y\ell }(z_{\ell },t)}{\partial z_{\ell }^2\partial t^2}+{\rho }_{\ell }{A}_{\ell }\frac{{\partial }^2{u}_{y\ell }(z_{\ell },t)}{\partial t^2}+2{\rho }_{\ell }{I}_{\ell }\Omega \frac{{\partial }^3{u}_{x\ell }(z_{\ell },t)}{\partial t\partial z_{\ell }^2}=\hfill \\ \hfill & q_y(z_{\ell },t)+{\rho }_{\ell }{A}_{\ell }\epsilon \left(z_{\ell }\right){\Omega }^{2}sin(\Omega t+\beta \left(z_{\ell }\right)),\hfill \end{array}\end{array}$$

(1b)

where E is the Young’s modulus, I is the diametric moment of area of the cross-section about the x- and y-axis, $\rho$ is the density, $q_{•}$ are the external distributed forces, A is the area of the cross-section, $\epsilon$ is the amount of eccentricity, and $\beta$ is its angular position. The rotations of the cross-section are defined by

$$\begin{array}{cc}\hfill {\phi }_{y\ell }(z_{\ell },t)& =\frac{\partial u_{x\ell }(z_{\ell },t)}{\partial z_{\ell }},\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill {\phi }_{x\ell }(z_{\ell },t)& =-\frac{\partial u_{y\ell }(z_{\ell },t)}{\partial z_{\ell }},\hfill \end{array}$$

(2b)

the bending moments are defined by

$$\begin{array}{cc}\hfill M_{y\ell }(z_{\ell },t)& =E_{\ell }{I}_{\ell }\frac{{\partial }^2{u}_{x\ell }(z_{\ell },t)}{\partial z_{\ell }^2},\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill M_{x\ell }(z_{\ell },t)& =-E_{\ell }{I}_{\ell }\frac{{\partial }^2{u}_{y\ell }(z_{\ell },t)}{\partial z_{\ell }^2},\hfill \end{array}$$

(3b)

and the shear forces can be computed by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Q_{x\ell }(z_{\ell },t)=& {\rho }_{\ell }{I}_{\ell }\frac{{\partial }^3{u}_{x\ell }(z_{\ell },t)}{\partial t^2\partial z_{\ell }}+2{\rho }_{\ell }{I}_{\ell }\Omega \frac{{\partial }^2{u}_{y\ell }(z_{\ell },t)}{\partial t\partial z_{\ell }}-E_{\ell }{I}_{\ell }\frac{{\partial }^3{u}_{x\ell }(z_{\ell },t)}{\partial z_{\ell }^3},\hfill \end{array}\end{array}$$

(4a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Q_{y\ell }(z_{\ell },t)=& {\rho }_{\ell }{I}_{\ell }\frac{{\partial }^3{u}_{y\ell }(z_{\ell },t)}{\partial t^2\partial z_{\ell }}-2{\rho }_{\ell }{I}_{\ell }\Omega \frac{{\partial }^2{u}_{x\ell }(z_{\ell },t)}{\partial t\partial z_{\ell }}-E_{\ell }{I}_{\ell }\frac{{\partial }^3{u}_{y\ell }(z_{\ell },t)}{\partial z_{\ell }^3}.\hfill \end{array}\end{array}$$

(4b)

The proposed method assumes a solution of the form

$$\begin{array}{cc}\hfill u_{x\ell }(z_{\ell },t)& ={\tilde{u}}_{x\ell }^{+}\left(z_{\ell }\right){\mathrm{e}}^{\mathrm{i}\omega t}+{\tilde{u}}_{x\ell }^{-}\left(z_{\ell }\right){\mathrm{e}}^{-\mathrm{i}\omega t},\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill u_{y\ell }(z_{\ell },t)& ={\tilde{u}}_{y\ell }^{+}\left(z_{\ell }\right){\mathrm{e}}^{\mathrm{i}\omega t}+{\tilde{u}}_{y\ell }^{-}\left(z_{\ell }\right){\mathrm{e}}^{-\mathrm{i}\omega t},\hfill \end{array}$$

(5b)

which leads to four ordinary differential equations, where $\omega$ is the complex eigenvalue. The solutions of $u_{•\ell }^{+}\left(z_{\ell }\right)$ and $u_{•\ell }^{-}\left(z_{\ell }\right)$ have to be complex conjugated to yield the real solutions for $u_{•\ell }(z_{\ell },t)$. Therefore, the differential equations for $u_{•\ell }^{+}\left(z_{\ell }\right)$ and $u_{•\ell }^{-}\left(z_{\ell }\right)$ are completely decoupled, and only two of these equations

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{{\partial }^4{\tilde{u}}_{x\ell }^{+}\left(z_{\ell }\right)}{\partial z_{\ell }^4}+{\rho }_{\ell }{\omega }^2\frac{{\partial }^2{\tilde{u}}_{x\ell }^{+}\left(z_{\ell }\right)}{E_{\ell }\partial z_{\ell }^2}\hfill \\ \hfill & -\frac{{\rho }_{\ell }{A}_{\ell }{\omega }^2}{E_{\ell }{I}_{\ell }}{\tilde{u}}_{x\ell }^{+}\left(z_{\ell }\right)-2\mathrm{i}{\rho }_{\ell }\Omega \omega \frac{{\partial }^2{\tilde{u}}_{y\ell }^{+}\left(z_{\ell }\right)}{E_{\ell }\partial z_{\ell }^2}=\frac{{\rho }_{\ell }{A}_{\ell }{\Omega }^2}{2E_{\ell }{I}_{\ell }}\epsilon \left(z_{\ell }\right){\mathrm{e}}^{\mathrm{i}\beta \left(z_{\ell }\right)},\hfill \end{array}\end{array}$$

(6a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{{\partial }^4{\tilde{u}}_{y\ell }^{+}\left(z_{\ell }\right)}{\partial z_{\ell }^4}+{\rho }_{\ell }{\omega }^2\frac{{\partial }^2{\tilde{u}}_{y\ell }^{+}\left(z\right)}{E_{\ell }\partial z_{\ell }^2}\hfill \\ \hfill & -\frac{{\rho }_{\ell }{A}_{\ell }{\omega }^2}{E_{\ell }{I}_{\ell }}{\tilde{u}}_{y\ell }^{+}\left(z_{\ell }\right)+2\mathrm{i}{\rho }_{\ell }\Omega \omega \frac{{\partial }^2{\tilde{u}}_{x\ell }^{+}\left(z_{\ell }\right)}{E_{\ell }\partial z_{\ell }^2}=-\frac{i{\rho }_{\ell }{A}_{\ell }{\Omega }^2}{2E_{\ell }{I}_{\ell }}\epsilon \left(z_{\ell }\right){\mathrm{e}}^{\mathrm{i}\beta \left(z_{\ell }\right)},\hfill \end{array}\end{array}$$

(6b)

are sufficient to solve the system of equations. If steady-state vibrations due to unbalance are analyzed, it is apparent that $\Omega =\omega$. External distributed loadings $q_{•}=0$ can also be neglected in this case.

2.2. Speed Dependent Parameters of Fluid Film Bearings

The stiffness and damping parameters of fluid film bearings depend on the geometry of the bearing, the viscosity of the lubricant, the static load of the bearing, and the spin speed of the rotor. For this paper, laminar flow, incompressible fluid, and short journal bearings are assumed. To find the equations to define the behavior of fluid-film bearings, an additional reference frame is introduced. The rotating reference frame $O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$ is coincident with the center of the bearing. The directions of the axes of this reference frame, whose axis $x^{\prime }$ contains the center of the bearing C, are not fixed in space if point C moves around point $O^{\prime }$, as is visualized in Figure 2.

First, the stationary operation of the bearing is described, where the center of the bearing is independent of time t and, therefore, the reference frame $O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$ is fixed in space. The pressure is linked to the thickness of the fluid by the Reynolds Equation [30]

$$\frac{1}{6}\left(\frac{1}{R_j^2}\frac{\partial }{\partial \theta }\right(\frac{h^3}{\mu }\frac{\partial p}{\partial \theta })+\frac{\partial }{\partial z}(\frac{h^3}{\mu }\frac{\partial p}{\partial z}\left)\right)=\Omega \frac{\partial h}{\partial \theta }+2\frac{\partial h}{\partial t},$$

(7)

where R is the Reynolds number, h the film thickness, p the pressure, and $\mu$ the viscosity of the lubricant. The film thickness is described as a function of the fixed coordinates and the clearance c

$$h=c(1\nabla \frac{x}{c}cos\left(\theta \right)\nabla \frac{y}{c}sin\left(\theta \right)),$$

(8)

where ∇ is the Nabla operator. In the case of short bearings, the flow in the circumferential direction may be neglected [31]. In combination with Equation (8), this leads to

$$\frac{1}{6}\left(\frac{\partial }{\partial z}\right(\frac{h^3}{\mu }\frac{\partial p}{\partial z}\left)\right)=(\Omega x\nabla 2\dot{y})sin\left(\theta \right)\nabla (\Omega y+2\dot{x})cos\left(\theta \right).$$

(9)

From Equation (9), the pressure distribution and the forces in static conditions ${\overrightarrow{f}}_{st}$ are obtained

$${\overrightarrow{f}}_{st}=\frac{\mu R_j{b}^3\Omega }{2c^2}\left[\begin{array}{c}{\int }_{{\theta }_1}^{{\theta }_2}\Lambda cos\left(\theta \right)\mathrm{d}\theta \\ \\ {\int }_{{\theta }_1}^{{\theta }_2}\Lambda sin\left(\theta \right)\mathrm{d}\theta \end{array}\right],$$

(10)

where

$$\begin{array}{ccc}\hfill & \Lambda =\frac{x_{st}sin\left(\theta \right)-y_{st}cos\left(\theta \right)}{c{h_{st}^{*}}^3},\hfill & \hfill h_{st}^{*}=1-\frac{x_{st}}{c}cos\left(\theta \right)-\frac{y_{st}}{c}sin\left(\theta \right),\end{array}$$

and b is the length of the bearing. After the stationary conditions have been defined, the linearized dynamic system is computed. Assuming that the bearing is displaced by only a small quantity $\Delta x,\Delta y$ from the static equilibrium position and moves with speed $\dot{x},\dot{y}$, the force of the fluid film on the journal is described as

$$\overrightarrow{f}=\overrightarrow{f_{st}}-{\left[\begin{array}{cc}\frac{F_x}{\dot{x}}& \frac{F_x}{\dot{y}}\\ \frac{F_y}{\dot{x}}& \frac{F_y}{\dot{y}}\end{array}\right]}_{x_{st},y_{st}}\left[\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right]-{\left[\begin{array}{cc}\frac{F_x}{x}& \frac{F_x}{y}\\ \frac{F_y}{x}& \frac{F_y}{y}\end{array}\right]}_{x_{st},y_{st}}\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right],$$

(11)

as long as the inertia of the oil film is neglected. The linearized behavior of the bearings for a given spin speed is expressed by eight coefficients (four stiffness coefficients $k_{•}$ and four damping coefficients $d_{•}$), which are used to describe the boundary and interface conditions of NAT. These values are obtained by

$$\begin{array}{cc}\hfill & {\left[\begin{array}{cc}{k}_{xx}& k_{xy}\\ k_{yx}& k_{yy}\end{array}\right]}_{x_{st},y_{st}}=-\frac{\mu R_j{b}^3\Omega }{2c^3}\left[\begin{array}{cc}{I}_1+I_2& -I_4+I_5\\ I_6+I_5& -I_1+I_3\end{array}\right],\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill & {\left[\begin{array}{cc}{d}_{xx}& d_{xy}\\ d_{yx}& d_{yy}\end{array}\right]}_{x_{st},y_{st}}=-\frac{\mu R_j{b}^3}{c^3}\left[\begin{array}{cc}{I}_4& I_1\\ I_1& I_6\end{array}\right],\hfill \end{array}$$

(12b)

where

$$\begin{array}{cccc}\hfill I_1& ={\int }_{{\theta }_1}^{{\theta }_2}\frac{sin\left(\theta \right)cos\left(\theta \right)}{{h_{st}^{*}}^3}\mathrm{d}\theta ,\hfill & \hfill I_2& =3{\int }_{{\theta }_1}^{{\theta }_2}\frac{\Lambda {cos\left(\theta \right)}^2}{h_{st}^{*}}\mathrm{d}\theta ,\hfill \\ \hfill I_3& =3{\int }_{{\theta }_1}^{{\theta }_2}\frac{\Lambda {sin\left(\theta \right)}^2}{h_{st}^{*}}\mathrm{d}\theta ,\hfill & \hfill I_4& ={\int }_{{\theta }_1}^{{\theta }_2}\frac{{cos\left(\theta \right)}^2}{{h_{st}^{*}}^3}\mathrm{d}\theta ,\hfill \\ \hfill I_5& =3{\int }_{{\theta }_1}^{{\theta }_2}\frac{\Lambda sin\left(\theta \right)cos\left(\theta \right)}{h_{st}^{*}}\mathrm{d}\theta ,\hfill & \hfill I_6& ={\int }_{{\theta }_1}^{{\theta }_2}\frac{{sin\left(\theta \right)}^2}{{h_{st}^{*}}^3}\mathrm{d}\theta .\hfill \end{array}$$

For a detailed description of these steps, the reader is referred to Genta [29] and Someya [32].

2.3. Boundary and Interface Conditions

To yield a unique solution, the governing equations require certain boundary conditions. Since Equation (6a,b) are fourth-order equations, eight boundary or interface conditions are necessary for each segment. The boundary conditions on the left end (station $\left(1\right)$) are defined by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & (-{\omega }^2{m}^{(1)}+j\omega d_{xx}^{(1)}(\Omega )+k_{xx}^{(1)}(\Omega )){\tilde{u}}_{x1}^{+}(0)+\hfill \\ \hfill & (+j\omega d_{xy}^{(1)}(\Omega )+k_{xy}^{(1)}(\Omega )){\tilde{u}}_{y1}^{+}(0)-{\tilde{Q}}_{x1}^{+}(0)=\frac{m^{(1)}{\tilde{\epsilon }}_{+}^{(1)}{\Omega }^2}{2},\hfill \end{array}\end{array}$$

(13a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & (-{\omega }^2{m}^{(1)}+j\omega d_{yy}^{(1)}(\Omega )+k_{yy}^{(1)}(\Omega )){\tilde{u}}_{y1}^{+}(0)+\hfill \\ \hfill & (+j\omega d_{yx}^{(1)}(\Omega )+k_{yx}^{(1)}(\Omega )){\tilde{u}}_{x1}^{+}(0)-{\tilde{Q}}_{x1}^{+}(0)=-j\frac{m^{(1)}{\tilde{\epsilon }}_{+}^{(1)}{\Omega }^2}{2},\hfill \end{array}\end{array}$$

(13b)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & -{\omega }^2{\Theta }_t^{\left(1\right)}{\tilde{\phi }}_{y1}^{+}\left(0\right)-j\Omega \omega {\Theta }_p^{\left(1\right)}{\tilde{\phi }}_{x1}^{+}\left(0\right)-{\tilde{M}}_{y1}^{+}\left(0\right)=0,\hfill \end{array}\end{array}$$

(14a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & -{\omega }^2{\Theta }_t^{\left(1\right)}{\tilde{\phi }}_{x1}^{+}\left(0\right)-j\Omega \omega {\Theta }_p^{\left(1\right)}{\tilde{\phi }}_{y1}^{+}\left(0\right)-{\tilde{M}}_{x1}^{+}\left(0\right)=0.\hfill \end{array}\end{array}$$

(14b)

The boundary conditions on the right end (station $\left(N\right)$) are considered analogously.

The interface conditions at station $\left(i\right)$ are defined by

$$\begin{array}{ccc}\hfill & {\tilde{u}}_{x\ell }^{+}\left(Z_i^{+}\right)-{\tilde{u}}_{x\ell -1}^{+}\left(Z_i^{-}\right)=0,\hfill & \hfill {\tilde{u}}_{y\ell }^{+}\left(Z_i^{+}\right)-{\tilde{u}}_{y\ell -1}^{+}\left(Z_i^{-}\right)=0,\end{array}$$

(15)

$$\begin{array}{ccc}\hfill & {\tilde{\phi }}_{y\ell }^{+}\left(Z_i^{+}\right)-{\tilde{\phi }}_{y\ell -1}^{+}\left(Z_i^{-}\right)=0,\hfill & \hfill {\tilde{\phi }}_{x\ell }^{+}\left(Z_i^{+}\right)-{\tilde{\phi }}_{x\ell -1}^{+}\left(Z_i^{-}\right)=0,\end{array}$$

(16)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & -{\omega }^2{\Theta }_t^{\left(i\right)}{\tilde{\phi }}_{y\ell }^{+}\left(Z_i^{+}\right)-j\Omega \omega {\Theta }_p^{\left(i\right)}{\tilde{\phi }}_{x\ell }^{+}\left(Z_i^{+}\right)-{\tilde{M}}_{y\ell }^{+}\left(Z_i^{+}\right)+{\tilde{M}}_{y\ell -1}^{+}\left(Z_i^{-}\right)=0,\hfill \end{array}\end{array}$$

(17a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & -{\omega }^2{\Theta }_t^{\left(i\right)}{\tilde{\phi }}_{x\ell }^{+}\left(Z_i^{+}\right)-j\Omega \omega {\Theta }_p^{\left(i\right)}{\tilde{\phi }}_{y\ell }^{+}\left(Z_i^{+}\right)-{\tilde{M}}_{x\ell }^{+}\left(Z_i^{+}\right)+{\tilde{M}}_{x\ell -1}^{+}\left(Z_i^{-}\right)=0,\hfill \end{array}\end{array}$$

(17b)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & (-{\omega }^2{m}^{(i)}+j\omega d_{xx}^{(i)}(\Omega )+k_{xx}^{(i)}(\Omega )){\tilde{u}}_{x\ell }^{+}(Z_i^{+})+(j\omega d_{xy}^{(i)}(\Omega )\hfill \\ \hfill & +k_{xy}^{(i)}(\Omega )){\tilde{u}}_{y\ell }^{+}(Z_i^{+})-{\tilde{Q}}_{x\ell }^{+}(Z_i^{+})+{\tilde{Q}}_{x\ell -1}^{+}(Z_i^{-})=-j\frac{m^{(i)}\epsilon (z_{\ell }){\mathrm{e}}^{\mathrm{i}\beta (z_{\ell })}{\Omega }^2}{2},\hfill \end{array}\end{array}$$

(18a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & (-{\omega }^2{m}^{(i)}+j\omega d_{yy}^{(i)}(\Omega )+k_{yy}^{(i)}(\Omega )){\tilde{u}}_{y\ell }^{+}(Z_i^{+})+(j\omega d_{yx}^{(i)}(\Omega )\hfill \\ \hfill & +k_{yx}^{(i)}(\Omega )){\tilde{u}}_{x\ell }^{+}(Z_i^{+})-{\tilde{Q}}_{y\ell }^{+}(Z_i^{+})+{\tilde{Q}}_{y\ell -1}^{+}(Z_i^{-})=-j\frac{m^{(i)}\epsilon (z_{\ell }){\mathrm{e}}^{\mathrm{i}\beta (z_{\ell })}{\Omega }^2}{2}.\hfill \end{array}\end{array}$$

(18b)

The locations $Z_i^{-}$ and $Z_i^{+}$ are infinitesimal to the left and right of station (i).

2.4. Homogeneous Solution

The general homogeneous solution of the governing Equation (6a,b) is found by setting the complex unbalance $\epsilon \left(z_{\ell }\right){\mathrm{e}}^{\mathrm{i}\beta \left(z_{\ell }\right)}$ to zero. Thus, assuming a solution in local coordinates $z_{\ell }$ of the form

$$\begin{array}{ccc}\hfill & {\tilde{u}}_{hx\ell }^{+}\left(z_{\ell }\right)=c_{ux\ell }{\mathrm{e}}^{ik{z}_{\ell }},\hfill & \hfill {\tilde{u}}_{hy\ell }^{+}\left(z_{\ell }\right)=c_{uy\ell }{\mathrm{e}}^{ik{z}_{\ell }}\end{array}$$

(19)

leads to a system of equations

$$\left[\begin{array}{cc}{k}^4-\frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}{k}^2-\frac{{\rho }_{\ell }{A}_{\ell }{\omega }^2}{E_{\ell }{I}_{\ell }}& 2\frac{\mathrm{i}{\rho }_{\ell }\Omega \omega }{E_{\ell }}{k}^2\\ -2\frac{\mathrm{i}{\rho }_{\ell }\Omega \omega }{E_{\ell }}{k}^2& k^4-\frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}{k}^2-\frac{{\rho }_{\ell }{A}_{\ell }{\omega }^2}{E_{\ell }{I}_{\ell }}\end{array}\right]\left[\begin{array}{c}{c}_{ux\ell }\\ c_{uy\ell }\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$

(20)

The non trivial solutions are found by setting the determinant to zero

$$(k^4-\frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}{k}^2-\frac{{\rho }_{\ell }{A}_{\ell }{\omega }^2}{E_{\ell }{I}_{\ell }}{)}^2-4\frac{{\rho }_{\ell }^2{\Omega }^2{\omega }^2}{E_{\ell }^2}{k}^4=0.$$

(21)

In the case of $\Omega =\omega$, the eight solutions of the characteristic equation are given by

$$\begin{array}{cc}\hfill & k_{1,2}=\pm \frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}\sqrt{\frac{1}{2}(\sqrt{1+\frac{4E_{\ell }{A}_{\ell }}{{\rho }_{\ell }{\omega }^2{I}_{\ell }}}-1)},\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill & k_{3,4}=\pm \frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}\sqrt{\frac{1}{2}(\sqrt{1+\frac{4E_{\ell }{A}_{\ell }}{{\rho }_{\ell }{\omega }^2{I}_{\ell }}}+1)},\hfill \end{array}$$

(22b)

$$\begin{array}{cc}\hfill & k_{5,6}=\pm \frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}\sqrt{\frac{3}{2}(\sqrt{1+\frac{4E_{\ell }{A}_{\ell }}{9{\rho }_{\ell }{\omega }^2{I}_{\ell }}}+1)},\hfill \end{array}$$

(22c)

$$\begin{array}{cc}\hfill & k_{7,8}=\pm \frac{{\rho }_{\ell }{\omega }^2}{E_{\ell }}\sqrt{\frac{3}{2}(\sqrt{1+\frac{4E_{\ell }{A}_{\ell }}{9{\rho }_{\ell }{\omega }^2{I}_{\ell }}}-1)}.\hfill \end{array}$$

(22d)

For $k_1-k_4$, the relation $c_{uy\ell }=-\mathrm{i}{c}_{ux\ell }$, and for $k_5-k_8$ the relation $c_{uy\ell }=\mathrm{i}{c}_{ux\ell }$ is valid. Therefore, the constants $c_{ux\ell }$ ans $c_{uy\ell }$ are not independent, and the general homogeneous solution is given by

$$\begin{array}{c}\hfill {\tilde{u}}_{hx\ell }^{+}\left(z_{\ell }\right)=c_{1\ell }cos\left({\alpha }_{1\ell }{z}_{\ell }\right)+c_{2\ell }sin\left({\alpha }_{1\ell }{z}_{\ell }\right)+c_{3\ell }{\mathrm{e}}^{{\alpha }_{2\ell }(z_{\ell }-L_{\ell })}+c_{4\ell }{\mathrm{e}}^{-{\alpha }_{2\ell }{z}_{\ell }}\\ \hfill +c_{5\ell }cos\left({\alpha }_{3\ell }{z}_{\ell }\right)+c_{6\ell }sin\left({\alpha }_{3\ell }{z}_{\ell }\right)+c_{7\ell }{\mathrm{e}}^{{\alpha }_{4\ell }(z_{\ell }-L_{\ell })}+c_{8\ell }{\mathrm{e}}^{-{\alpha }_{4\ell }{z}_{\ell }},\end{array}$$

(23a)

$$\begin{array}{c}\hfill {\tilde{u}}_{hy\ell }^{+}\left(z_{\ell }\right)=-\mathrm{i}(c_{1\ell }cos\left({\alpha }_{1\ell }{z}_{\ell }\right)+c_{2\ell }sin\left({\alpha }_{1\ell }{z}_{\ell }\right)+c_{3\ell }{\mathrm{e}}^{{\alpha }_{2\ell }(z_{\ell }-L_{\ell })}+c_{4\ell }{\mathrm{e}}^{-{\alpha }_{2\ell }{z}_{\ell }}\\ \hfill +c_{5\ell }cos\left({\alpha }_{3\ell }{z}_{\ell }\right)+c_{6\ell }sin\left({\alpha }_{3\ell }{z}_{\ell }\right)+c_{7\ell }{\mathrm{e}}^{{\alpha }_{4\ell }(z_{\ell }-L_{\ell })}+c_{8\ell }{\mathrm{e}}^{-{\alpha }_{4\ell }{z}_{\ell }}),\end{array}$$

(23b)

where $L_{\ell }$ is the length of the segment, $c_{1\ell }-c_{8\ell }$ are arbitrary constants, and ${\alpha }_{•}$ are the real-valued wave numbers ${\alpha }_{1\ell }=k_1,{\alpha }_{2\ell }=\mathrm{Im}\left(k_3\right),{\alpha }_{3\ell }=k_5$, and ${\alpha }_{4\ell }=\mathrm{Im}\left(k_7\right)$, with $\mathrm{Im}(•)$ as the imaginary part. To increase the stability of the numeric implementation, exponential terms are implemented instead of the commonly used $sinh(•)$ and $cosh(•)$ terms. With this approach, the amplitude of each function remains smaller or equal to 1 within each segment. According to Equations (2a,b)–(4a,b) the rotation, moments, and shear forces of the beam are found. The state within a rotor segment is defined by the matrix equation

$${\overrightarrow{x}}_{h\ell }\left(z_{\ell }\right)={\underline{\underline{B}}}_{\ell }\left(z_{\ell }\right)\overrightarrow{c_{\ell }},$$

(24)

where ${\overrightarrow{x}}_{h\ell }$ is the homogeneous solution vector, ${\underline{\underline{B}}}_{\ell }$ is the state variable matrix, and $\overrightarrow{c_{\ell }}$ is the vector of arbitrary constants.

2.5. Particular Solution

Unbalances are either concentrated on singular points or distributed along the whole length of the rotor described by any arbitrary distribution function. In the case of arbitrarily distributed unbalance functions, an approximation of the load is given by the Fourier extension method [33]. This method has the advantage of not suffering from the Gibbs phenomenon, unlike the classical Fourier series expansion. Therefore, the convergence rate of the method is very high [33] and very efficient implementations are also available [34]. For a detailed description, the reader is referred to Klanner et al. [25].

2.6. Assembly and Solution Procedure

The total solution vector of a rotor segment is defined as ${\overrightarrow{x}}_{\ell }\left(z_{\ell }\right)={\overrightarrow{x}}_{h\ell }\left(z_{\ell }\right)+{\overrightarrow{x}}_{p\ell }\left(z_{\ell }\right)$. The total solutions of the segments together with the boundary and interface conditions lead to a system of linear equations

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

(25)

where $\underline{\underline{A}}$ is the system matrix, which depends only on the concentrated elements at the stations. The right-hand side vector $\overrightarrow{b}$ also depends on the unbalance. The unknown constants $\overrightarrow{c}$ are the solution of the system of linear equations, which defines the state variables of the whole rotor.

2.7. Recursive Eigenvalue Search

In this section, a new approach to find the eigenvalues and the Campbell diagram of a NAT system matrix is implemented. At every eigenvalue, the coefficient matrix of the NAT model is singular. The prevalent approach in classical root search strategies is based on a sign change detection of the determinant of the coefficient matrix. If the step size in this approach is too large, eigenvalues, which are close to each other, or double roots might be missed.

The recursive algorithm of Bestle et al. [35] solves this problem by reducing the zero search to a minimization problem. At each local minimum of the absolute value of the determinant, the algorithm is repeated in an area of one step size in both directions until the desired accuracy is met. The algorithm searches in logarithmic steps, which is better suited for wide spread frequencies. In the case of a damped rotor-bearing system, the zero search problem depends on two variables, the damped critical speed ${\omega }_d$ and the modal damping coefficient d, which are the imaginary and the real part of the complex eigenvalues. In this case, the algorithm scans a matrix of those two variables for values, which are a local minimum in both directions [35].

For the generation of Campbell diagrams, the algorithm is adapted in a way that the eigenvalue search is performed at a defined spin speed. To increase numeric efficiency, the search area for the eigenvalues at a specific excitation frequency is chosen to be the value of this eigenfrequency at the previous excitation frequency plus or minus a small percentage, depending on the step size. Since Campbell diagrams consist of continuous lines, the efficiency can be increased significantly with this approach, without missing eigenvalues.

2.8. Influence Coefficient Method

The influence coefficients ${\alpha }_{ik}(\Omega )$ are defined as

$${\alpha }_{ik}(\Omega )=\frac{x_i(\Omega )}{U_k},$$

(26)

where $x_i(\Omega )$ is the displacement at position i caused by an unbalance $U_k$ at position k. Displacement, influence coefficients, and unbalances are usually complex, whereas the ratio of the imaginary and the real part contains the information of the angular position. If the influence coefficients of many combinations of unbalance planes k and measurement planes i are calculated, this input-output relationship can be described as a system of linear equations

$$\overrightarrow{x}=\underline{\underline{\alpha }}\overrightarrow{U},$$

(27)

where $\overrightarrow{x}$ are the displacements of the rotor bearing system caused by the addition of test weights, $\underline{\underline{\alpha }}$ is the matrix of influence coefficients, and $\overrightarrow{U}$ is the unbalance vector. If the vibrations of the system are measured, and the matrix of the influence coefficients is calculated, the unbalance of the system can easily be computed by inverting the influence coefficient matrix

$$\overrightarrow{U_c}=-{\underline{\underline{\alpha }}}^{-1}\overrightarrow{x_0},$$

(28)

where $\overrightarrow{U_c}$ are the correction unbalances, and $\overrightarrow{x_0}$ are the displacements of the rotor bearings system without test weights. Thus, the system is balanced by mounting correction masses proportional to the calculated unbalance on the opposite side of the rotor in the according balancing planes. Surplus information, on additional measurement planes or spin speeds, may be used to find the ideal balancing positions or reduce measurement errors with the least squares method [36].

3. Results

To show the viability of the balancing method, a rotor-bearing system was simulated in two different states of unbalance. At first, the critical speeds, Campbell diagram, and frequency response function (FRF) were compared with the well-known FEM code by Friswell et al. [37] to show the accuracy of the NAT simulation. Then, the influence coefficients of the system were calculated using the proposed method. Usually, the displacements of the rotor-bearing system are measured on site. As proof of concept, the displacements of the rotor-bearing system were calculated with the Friswell code in this paper. With these displacements and the influence coefficients calculated by NAT, the balancing weights were found, and the balancing success was analyzed.

In the first test, the unbalance was limited to the disks mounted on the shaft. Since all disks are used as balancing planes, the proposed quasi-analytical method was expected to find the exact unbalance of the system. The second test also included the arbitrarily distributed unbalance of the shaft. With a finite number of balancing planes, the vibrations of the system cannot be eliminated entirely for every position and every spin speed but are substantially reduced. The average amplitude of the vibrations before and after balancing were compared to determine the balancing success. All calculations were performed on an Intel^® Core^TM i7-8700 CPU with 3.2 GHz running on a Windows 10 operating system using MATLAB^TM version R2019a.

3.1. Rotor Bearing System

The rotor-bearing system consisted of a stepped shaft with three disks of different weights, which was supported on two fluid film bearings, as shown in Figure 3. This geometry was chosen because it included all relevant aspects of the rotor problem including gyroscopic effects, changing shaft diameter, and the behavior of fluid film bearings. The parameters of the stations are listed in

Table 1. Stations of the numerical example.

z	m	${\mathbf{\Theta }}_{\mathit{t}}$	${\mathbf{\Theta }}_{\mathit{p}}$	F	b	${\mathit{d}}_{\mathit{i}}$	$\mathit{\mu }$	c
m	kg	${\mathrm{kgm}}^2$	${\mathrm{kgm}}^2$	N	m	m	$\frac{\mathrm{Ns}}{{\mathrm{m}}^2}$	m
0	0	0	0	0	0	0	0	0
0.13	0	0	0	431	0.04	0.02	0.032	$8\times {10}^{-5}$
0.31	14.7	0.0412	0.0735	0	0	0	0	0
0.47	23	0.0966	0.0180	0	0	0	0	0
0.71	33	0.1960	0.3722	0	0	0	0	0
0.91	0	0	0	380	0.04	0.02	0.032	$8\times {10}^{-5}$
1.0	0	0	0	0	0	0	0	0

, where F is the static load on the fluid film bearing, b the length of the fluid film bearing, $d_i$ its inner diameter, c its clearance, and $\mu$ the viscosity of the lubricant. The disks were considered stiff. The stepped shaft had a density of 7800 $\frac{\mathrm{k}\mathrm{g}}{{\mathrm{m}}^3}$ and a Young’s modulus of $2.1\times {10}^{11}\frac{\mathrm{N}}{{\mathrm{m}}^2}$.

The fluid film bearing was simulated with cross-coupling effects and speed-dependent stiffness and damping. The characteristics of the stiffness and damping of the right fluid film bearing are shown in Figure 4. Dashed lines represent negative values.

The NAT simulation was compared with an FEM model built using the Friswell code, where the shaft was modeled with 100 elements á 0.01 m. Due to the discretization, FEM simulations with few elements tend to calculate higher natural frequencies than analytic calculations. This effect is reduced and nearly eliminated with increasing element count [29]. Therefore, the critical speeds calculated with the FEM model were slightly higher than those found by the quasi-analytical NAT simulation. As is seen in the Campbell diagram (Figure 5) and

Table 2. Flexural critical speeds in Hz.

NAT	FEM	Relative Error
43.014	43.015	$0.0031\%$
55.059	55.060	$0.0018\%$
154.28	154.29	$0.0031\%$
160.75	161.76	$0.6275\%$

, the differences between the calculations were negligible.

3.2. Influence Coefficient Matrix

The influence coefficients were calculated with an NAT model of the rotor-bearing system at 2000 rpm, without any information about the actual unbalance of the system.

$$\underline{\underline{\alpha }}=\left[\begin{array}{ccc}0.0126-0.0078i& 0.0135-0.0082i& 0.0102-0.0073i\\ 0.0135-0.0082i& 0.0161-0.0091i& 0.0134-0.0088i\\ 0.0102-0.0073i& 0.0134-0.0088i& 0.0136-0.0092i\end{array}\right]$$

The measurement positions and also the balancing positions were the planes of the disks at the positions $z=$ 0.31 m, 0.47 m, and 0.71 m. Since the matrix is symmetric, the influence coefficients fulfilled the theorem of Maxwell and Betti [38].

3.3. Concentrated Unbalance

In case A, the unbalance was only concentrated at the disks of the rotor bearing system, according to

Table 3. Unbalance of case A.

z	U	$\mathit{\beta }$
m	kg m	rad
0.31	$2.8229\times {10}^{-4}$	0
0.47	$5.5135\times {10}^{-4}$	$\frac{\pi }{2}$
0.71	$2.5994\times {10}^{-3}$	$\pi$

.

The Frequency Response Function (FRF) generated by the proposed method was congruent with the FRF calculated by the Friswell code, as shown in Figure 6. Displacement $u_x$ and $u_y$ represented the average amplitude at all balancing positions.

Figure 7 and

Table 4. Average relative error in %.

Number of Elements	x-Direction	y-Direction
6	$1.2965\times {10}^{-2}$	$1.3422\times {10}^{-2}$
24	$5.6796\times {10}^{-5}$	$5.6677\times {10}^{-5}$
100	$1.1306\times {10}^{-7}$	$1.3352\times {10}^{-7}$

show the relative error of the displacement between NAT and three FEM models, the aforementioned model with 100 elements, one with 24 elements, and one with six elements. It is clearly visible that the numeric FEM solutions converged to the quasianalytic NAT solution with an increasing amount of elements. The relative error was calculated as the absolute value of the difference in amplitude between the two simulation methods, divided by the value of the NAT model.

The proposed balancing procedure was used on the rotor system at 2000 rpm. This led to the balancing weights listed in

Table 5. Balancing weights of case A.

z	U	$\mathit{\beta }$
m	kg m	rad
0.31	$2.8229\times {10}^{-4}$	$\pi$
0.47	$5.5135\times {10}^{-4}$	−$\frac{\pi }{2}$
0.71	$2.5994\times {10}^{-3}$	0

, which were equal to the unbalance of the disks rotated 180°. Since the method was quasianalytic and the balancing occurred at all unbalance positions, mounting these weights completely eliminated the vibration of the numeric example. The computation time for the calculation of the influence coefficients and balancing weights was 0.08752 s.

3.4. Distributed Unbalance

In case B, the unbalance of the shaft was also taken into account. It was distributed along the whole length of the shaft, and the function of the amount of the eccentricity was arbitrarily chosen to be

$$\epsilon \left(z\right)=-7\times {10}^{-6}\times z+{10}^{-5}\times log(z+1)+8\times {10}^{-6}\times sin\left(4\pi z\right).$$

(29)

The direction of the unbalance also changed in a corkscrew manner

$$\beta \left(z\right)=2\pi \times z,$$

(30)

which is considered to be the most difficult case to balance [39]. Figure 8 shows a visualization of the eccentricity of the shaft.

For a better balancing result, the system was now balanced at 2000, 5000, 10,000, and 20,000 rpm, with a total computation time of 0.24256 s. At each speed, an influence coefficient matrix and a set of balancing weights were calculated. The arithmetic mean of these balancing weights was used as the final set of balancing weights, as shown in

Table 6. Balancing weights of case B.

z	U	$\mathit{\beta }$
m	kg m	rad
0.31	$3.0180\times {10}^{-4}$	−3.0885
0.47	$5.2940\times {10}^{-4}$	−1.5469
0.71	$2.6122\times {10}^{-3}$	0.0060

.

Mounting these balancing weights resulted in an average reduction in the amplitude of vibration of 99.07%, measured across all balancing positions and at every spin speed between 1 rpm and 20,000 rpm, as seen in Figure 9.

4. Discussion

The experimental procedure to determine influence coefficients requires several test runs with stops, restarts, and cooldowns of the machine, which can be very time-consuming and expensive [36]. Modal balancing methods, on the other hand, which require fewer measurements, are not properly applicable on rotor systems supported on fluid film bearings [4]. For this reason several approaches have been made in recent years to substitute the measurements of influence coefficients with an accurate calculation.

In this paper, NAT was extended to calculate the influence coefficients of rotor systems supported on fluid film bearings without test runs. Whereas the commonly used FEM has a tradeoff between the numerical speed and accuracy of the simulation [29], NAT led to quasianalytic solutions and was very computationally efficient. This makes NAT especially suited for applications with high rotational speed, where FEM requires a fine discretization to model higher modes accurately. However, a substantial advantage of FEM is its versatility: It is also used to estimate the behavior of the foundation [36], which is not yet feasible with NAT. In further research, this problem will be addressed, by supplementing the proposed method with other state-of-the-art techniques. The method of the complex dynamic matrix of compliances is well suited to calculate the foundation and support behavior of rotating systems [40]. An alternative is the calculation of these parameters with the polynomial chaos expansion, which needs additional measurements, but has already been shown to work with input data from NAT [41]. The investigation of the proposed balancing method in practice, including the influence of the foundation, is the subject of ongoing research.

Artificial neural networks have been used in recent years, not only to monitor dynamic systems and make maintenance decisions [42] but also to improve the balancing process itself in combination with FEM [43]. The high efficiency of NAT makes it well suited for the high computational demand of these techniques and might lead to even better balancing results.

5. Conclusions

In this paper, NAT was extended to include speed-dependent stiffness and damping coefficients of rotor-bearing systems supported on fluid film bearings. Furthermore, a new modification of a recursive search algorithm to efficiently generate Campbell diagrams with NAT was implemented. The results obtained by FEM converged with an increasing amount of nodes to the results obtained by the presented quasi-analytical method and were in good agreement. An FEM calculation of 100 elements showed only an average difference of $1.3352\times {10}^{-7}$ of the amplitude of vibration. An implementation of the proposed simulation technique for the balancing of flexible rotors by means of calculated influence coefficients was presented. This balancing method was successfully applied in two test cases. The influence of multiple disks, stepped shafts, distributed unbalance, and the characteristics of fluid film bearings, such as cross-coupling effects and speed-dependent parameters, were taken into account. The unbalance was identified exactly in the first test as expected for a quasi-analytical method. In the second test, the eccentricity of the shaft was included, and the balancing procedure led to an average reduction in the amplitude of vibration of 99.07% measured across all balancing positions and relevant spin speeds. The simulation via NAT was very computationally efficient, with computation times of 0.08752 s and 0.24256 s for the test cases. The results indicate that the proposed balancing method can be applied effectively.