Stator-Rotor Contact Force Estimation of Rotating Machine

: In turbomachines, dry friction resulting from stator–rotor contacts is a severe problem that may degrade lifetime of the machine or even lead to complete failure. Knowledge about the system states and contact forces is beneﬁcial for system monitoring or to prevent contacts through, e.g., active magnetic bearings. In this paper, a nonlinear model is derived that describes the lateral rotor vibrations in the case of contact and no contact. The elastic behavior of the shaft is modeled based on the ﬁnite-element method. The contact is described by a dry friction model. An augmented system description is formulated that allows estimation of rotor displacements and contact forces by means of nonlinear ﬁltering approaches like an extended Kalman ﬁlter. A simulation study was conducted that explicitly considered the hazardous backward whirl. The suggested approach shows suitable estimation performance related to both state and contact force estimation.


Introduction
In turbo-machines like compressors, turbines or pumps, minimization of the clearance between rotor and stator is desirable to maximize energy efficiency [1,2]. However, vibration of the rotor may occur due to unbalance, trapped fluids, shaft cracks, thermal bends. etc. [3]. As a consequence, rub between stator and rotor may appear. The rub can be specified as point rub, partial rub, or full annular rub [2]. Among these kinds of rubs, full annular rub is considered most dangerous. It may lead to effects known as forward and backward whirling. Forward whirl denotes continuous rubbing in direction of the shaft rotation, whereas during backward whirling, the rubbing occurs in the opposite direction of the shaft rotation. The stator-rotor contacts may lead to heating of the shaft and occurrence of hot spots [1,2]. Generally, the rubbing degrades the lifetime of the machine. It may lead to damage of the contact surface, large deformation, high frequency stress, and even to complete failure [4,5].
Several papers have been proposed to study the dynamics induced by stator-rotor contacts. A simple but still recent modeling approach is to consider the shaft to be linearelastic and to model the contact with the stator based on a linear dry friction model [6]. Although the model is simple, it becomes nonlinear in the case of contact, inducing effects such as chaos, bifurcation, whirling, and wiping [6]. Conditions of existence for forward and backward whirls depending on system parameters have been investigated [5]. For instance, an overview about different effects resulting from stator-rotor contacts can be found in [1,7]. Experimental validation of the described effects has been undertaken in [8]. A study of the rotordynamics of more complex rotating systems with multiple rotor-stator interface locations is conducted in [9].
Estimation of rotor-stator contact forces has not been considered in the literature so far. However, to predict the remaining lifetime and for the purpose of system monitoring, estimation of the contact forces could be beneficial.
A promising strategy to prevent whirl responses is the usage of active elements such as active bearings [3]. In [3], a state-feedback controller is applied to stabilize the rotor displacements and mitigate effects resulting from contacts. Due to the special structure of the feedback controller, the approach does not assume the system states to be measured, but it requires all interacting forces between rotor and the surroundings to be measured. To avoid additional sensors, an estimation strategy could be applied. Another active bearing control approach based on H ∞ control design is developed in [12]. It requires several rotations and displacements to be measured, which could be reduced by applying estimation approaches.
From the perspective of filtering, the estimation of unknown contact forces is an unknown input estimation problem. Related to linear stochastic systems, two main estimation approaches exist: the minimum-variance unbiased (MVU) filtering approach and the augmented state Kalman filtering (ASKF) approach. The MVU approaches provide unbiased state estimations with minimum error variance under the influence of unknown inputs. For this purpose, the Kitanides filter is the first MVU estimation approach that has been derived [13]. However, it does not explicitly estimate the unknown inputs itself. Later, an MVU estimator known as Gillijns-De-Moore filter was proposed that provides explicit estimation of the unknown inputs [14]. For this filter, it is shown that, in addition to the state estimations, the estimations of the unknown inputs are unbiased and have minimum error variance. The fundamental property of the MVU approaches is that they neither know nor assume any dynamics for the unknown inputs. This is in contrast to the ASKF approaches. In augmented state Kalman filtering, the system model is augmented by the dynamics of the unknown inputs. This was first described in [15] for the treatment of a constant bias. However, in the case of time-varying inputs, the related dynamics of the unknown inputs are typically not known and have to be approximated by, e.g., a piecewise constant model [16]. The variance of the process noise related to the augmented state influences the estimation performance of the filter. It is a design parameter that is typically tuned by trial and error. In [17], it was recently shown that a direct connection between the MVU and ASKF approaches exists. If the considered noise variance related to the augmented state is selected as infinity, then the ASKF algorithm equals the Gillijns-De-Moore filter. Consequently, the MVU approaches appear as a special case of augmented-state Kalman filtering.
In this paper, an estimation strategy to estimate states and contact forces in rotating machines is proposed. First, a linear-elastic model based on the finite-element method is derived to describe the elastic movement of the shaft. The contact with the stator is described based on a dry friction model. As a consequence, a nonlinear augmented state description is obtained that exactly describes the system dynamics in case of contact and no contact. No assumptions about the dynamics of the unknown input are required to be made. By applying a nonlinear filtering approach to the augmented state description, the displacements and bendings, as well as the contact forces can be estimated. As the augmented state description is an exact description, an optimal estimation approach like the particle filter would lead to optimal estimation results of the unknown contact force. However, due to the large number of states and the resulting computational costs, particle filtering is considered to be unapplicable. Instead, the extended Kalman filter is applied, and the related Jacobians are derived. In the conducted simulation example, contact force estimation under the hazardous backward whirl is explicitly considered.
The paper is organized as follows. The modeling of the elastic shaft and the contact force is described in Section 2. In Section 3 the estimation strategy to estimate the states and contact forces is proposed. The contact force estimation under forward and backward whirl is considered in the simulation study of Section 4.
Throughout the paper I n denotes an identity matrix of dimension n × n and 0 n×m denotes a zero matrix of dimension n × m.

Modeling of Rotating Shaft
In this section, a rotating shaft is modeled as an elastic beam. As visualized in Figure 1a, the shaft is assumed to rotate around the z-axis. The mass of the shaft is considered to be unbalanced, so that the resulting unbalance force leads to radial displacement of the beam in x and y-directions. Based on finite element method [18], the dynamics of displacements of the shaft can be modeled as where M ∈ R 36×36 , D ∈ R 36×36 , K ∈ R 36×36 describe the mass, damping, and stiffness of the shaft; d ∈ R 36 comprises the coordinates of the nodes; F u ∈ R 2 with its input matrix B u ∈ R 36×2 describes excitation resulting from the unbalance force; and G = mg ∈ R with input matrixB g ∈ R 36 describes excitation resulting from gravity, where m ∈ R is the mass of the shaft without the unbalance mass. The vector d ∈ R 36 is defined as where θ y i , θ x i ∈ R are angles and y i , x i ∈ R are displacements around the y and x axis at node i ( Figure 1b). The derivation of the matrices M, D, K can be found in Appendix A. The unbalance force and its input matrix are given as where m u ∈ R is the unbalance mass, γ ∈ R is the distance between the unbalance mass and the geometric center, and ϕ ∈ R and Ω ∈ R are the rotating angle and angular velocity of the shaft. The gravity is distributed over all nodes leading to B g = 1 9 1 0 0 0 1 0 0 0 · · · 1 0 0 0 . (4)

Modeling of Stator-Rotor Contact
The unbalance force leads to radial displacement of the nodes. If the displacement exceeds the clearance, contact between rotor and stator occurs. It is assumed that the contact occurs at node five.
As visualized in Figure 2a, the distance d c ∈ R between the stator and rotor at node five is given by where r r ∈ R is the radius of the rotor, r s ∈ R is the inner radius of the stator, r 0 ∈ R is the clearance, and d 5 ∈ R is the displacement of node five, known as It follows from (5) that d c ≥ 0 characterizes contact and d c < 0 characterizes no contact. In the literature, it is common to model rotor-stator, rotor-bearing, or rotor-casing contact forces by spring and damper elements [5,6,19,20]. In this contribution, the bearing forces are modeled by springs and dampers but without clearance and without rub. The focus is on the contact of rotor and stator with rub. Following [5,6], this contact is modeled by a set of radial springs as visualized in Figure 2. The springs generate the normal force F N =kd c ∈ R, wherek = k s 10 9 > 0 ∈ R denotes the constant stiffness of the springs and k s is the prefactor of the spring stiffness. The factor 10 9 results from high values of the Young modulus E Y > 10 9 for the materials in contact. The force of friction is given by the dry friction model F f = µF N ∈ R with friction coefficient µ > 0 ∈ R. In addition, lubrication of the contact surfaces can be considered by proper choice of the friction coefficient. The contact force F c ∈ R 2 in x − y coordinates of the inertial coordinate system is where v c,x (y 5 , x 5 ,ẏ 5 ,ẋ 5 ) =ẋ 5 cos(ψ) −ẏ 5 sin(ψ) + Ωr r , is the velocity of the shaft at the contact point projected on axis x of the body-fixed coordinate system (Figure 2b). The final model that comprises occurance of rotor stator contact is given by is the input matrix of the contact force and Θ ∈ R is one in the case of contact (d c ≥ 0) and zero otherwise (d c < 0). By rearranging (9) and assuming full rank of M, the state-space representationη with η ∈ R 72 , u ∈ R 3 , ∈ R 72 , A ∈ R 72×72 , B ∈ R 72×3 , and N ∈ R 72×72 can be achieved. The state vector η is defined as η = d ḋ , and the input vector u is The system matrix A, the input matrix B of known respectively measured inputs, and the input matrix N related to the contact force are Let η(i) denote the i-th element of η. The nonlinear term (η, k s , µ) is given by In the following, discretization of the model is discussed. From the well-known discretization rules of linear time-invariant systems (see, e.g., [21]) it follows that the system matrix and input matrix can be discretized as where A is assumed to have full rank and t s denotes the sampling time. The contact force related term (η, k s , µ) is assumed to stay constant between two time steps (k + 1)t s and kt s which is a suitable approximation if the sampling time is sufficient small. For the piece-wise constant input (η, k s , µ) the corresponding input matrix N can be discretized as Finally, the nonlinear discrete-time model is achieved.
The Kalman filter is applied to estimate the states of the augmented system. Quantitieŝ ζ k ,d 5,k ,Θ k ,F cy,k+1 ,F cx,k+1 denote the estimation of ζ k , d 5,k , Θ k , F cy,k+1 , F cx,k+1 . First, the estimated displacementd of node five is determined. It is checked if the known clearance r 0 is exceeded or not, i.e., The a priori state estimationζ is obtained from the nonlinear augmented model f . The input vector u k is known and equals u from (12) at time step k. The a priori error covariance is determined as where Λ k denotes the Jacobian known from (20). The a posteriori state estimation is calculated aŝ where y k+1 are the measured node displacements and K k+1 denotes Kalman filter gain according to The a posteriori error covariance is obtained from According to (7), the estimation of the contact forces at time step k + 1 is given by

Results
In the following, the lateral behavior of a rotating shaft under the influence of rotorstator contact forces is simulated. The system parameters used in the simulation are summarized in Table 1.
The behavior of the shaft during forward and backward whirling is considered with the aim to estimate the appearing contact forces. During forward whirling, the angular velocities Ω and Ω of Figure 2a have the same algebraic signs, whereas during backward whirling, Ω and Ω have opposite signs. The lowest eigenfrequency of the shaft is 23.71 [1/s]. The rotational speed Ω of the shaft is increased according to A simulation duration of 5 [s] is considered. During that time, the shaft runs through its first eigenmode. The resulting resonance peak leads to stator-rotor contacts. Dependent on the values µ, k s of the dry friction model, forward or backward whirling occurs. The simulated measurement noise has standard deviation σ r = 1 × 10 −3 [mm] and the sampling time of the filter is chosen as T s = 0.001 [s]. The initial states and error covariance arê ζ = 0 74×1 , P 0 = I 74 × 0.0001. The design parameters σ k s , σ µ are chosen as σ k s = 100, σ µ = 100.
Selecting µ, k s as µ = 0.2, k s = 1 establishes a forward whirl as shown in Figure 3. The unknown radial displacement of node five is estimated well by the filter. The exceedance of the clearance can be correctly detected from the estimated displacements. The estimated contact forces are shown in Figure 4. The forward whirl remains stable until the end of the simulation. A strong vibrational effect can be detected. For both x− and y−components of the contact force, shape as well as amplitude of the true and estimated signal closely fit to each other.
Selecting µ, k s as µ = 0.5, k s = 1 establishes a backward whirl as visualized by Figure 5. It can be seen that the displacement of node five is correctly estimated so that the exceedance of the clearance can be detected. The backward whirl remains stable for the remaining simulation time. In comparison to the forward whirl, more than four times higher amplitude values can be observed, which also vary in time. The estimations of the contact forces are visualized in Figure 6. Again, shape and amplitude of true and estimated signals closely fit to each other.   In the following, the contact force estimation performance is quantified based on the measures ∀F cx,k = 0 : ∀F cy,k = 0 : For both scenarios, the performance measures are stated in Table 2. The performance measures are sensitive to relatively high estimation errors occuring when |F cx,k | or |F cy,k | are small. However, from Figures 4 and 6, it is known that the amplitudes of the contact forces are large most of the time. For that reason, the performance values (44), (45) are evaluated for all |F cy,k | > 5000 and |F cx,k | > 5000 additionally.

Conclusions
In this paper, an estimation strategy suitable to estimated forces resulting from statorrotor contacts is proposed. An augmented system description is developed that exactly describes the system dynamics for the considered type of contact force. The extended Kalman filter (EKF) is applied to an augmented system description to estimate the system states including contact forces. The signals of the estimated and true contact forces closely fit to each other regarding shape and amplitude. Apart from the EKF, other nonlinear filtering approaches may be applied to the proposed model in the future. In addition, plastic deformation could be considered.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available because the data can be easily generated by the reader based on the descriptions of the paper. considered beam elements can be found in [18] for instance. The lateral displacement of the rotor shown in Figure 1a is described by the following model where the forces resulting from the bearings are not yet considered. The eight beam elements form the mass M 0 ∈ R 36×36 , stiffness K 0 ∈ R 36×36 , and damping matrices D 0 ∈ R 36×36 according to and M i , K i ∈ R 8×8 given as The quantities used are: the material density ρ ∈ R, the cross section of the i-th beam element Σ i ∈ R, the length of the i-th beam element l i ∈ R, the Young modulus E Y ∈ R, the second moment of area of the i-th beam element J i ∈ R, and the damping factors α, β ∈ R. The transformation matrix T A ∈ R 64×36 describes the relation between the generalized coordinates of the eight beam elements v ∈ R 64 v = y 1l x 1l θ y 1l θ x 1l y 1r x 1r θ y 1r θ x 1r y 2l x 2l θ y 2l θ x 2l y 2r x 2r θ y 2r θ x 2r · · · y 8l x 8l θ y 8l θ x 8l y 8r x 8r θ y 8r θ x 8r T , and the displacements and bending angles of the nine nodes d ∈ R 36 are known from (2).
In the following, the presently missing forces resulting from the bearings are considered. The excitation p b ∈ R 36 of system (A1) resulting from the bearings can be described by with spring stiffness k b ∈ R and damper coefficient d b ∈ R. The location of the forces are defined by Π ∈ R 36×36