On the modeling of rotors with rolling element bearings using bond graphs

Journal Pre-proof On the modeling of rotors with rolling element bearings using bond graphs Jone Torsvik, Eilif Pedersen

PII:

S0022-460X(20)30036-5

DOI:

https://doi.org/10.1016/j.jsv.2020.115205

Reference:

YJSVI 115205

To appear in:

Journal of Sound and Vibration

Received Date: 15 August 2019 Revised Date:

15 January 2020

Accepted Date: 20 January 2020

Please cite this article as: J. Torsvik, E. Pedersen, On the modeling of rotors with rolling element bearings using bond graphs, Journal of Sound and Vibration (2020), doi: https://doi.org/10.1016/ j.jsv.2020.115205. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

CRediT author statement Jone Torsvik: Conceptualization, Methodology, Software, Formal analysis, Investigation, Visualization, Writing - Original Draft. Eilif Pedersen: Conceptualization, Methodology, Software, Formal analysis, Supervision, Writing - Review & Editing.

On the Modeling of Rotors with Rolling Element Bearings Using Bond Graphs Jone Torsvik, Eilif Pedersen Norwegian University of Science and Technology (NTNU) Department of Marine Technology Otto Nielsens vei 10 7491 Trondheim, Norway [email protected] [email protected]

Abstract This paper describes the development of a generic rotordynamic model incorporating rolling element bearings. Aimed at examining model fidelity requirements for bearings, the work contributes to advancing knowledge on rotordynamic system models suitable for studying mechatronics- and controls-related applications, transient issues and fault conditions. The bond graph method is applied in a consistent manner, giving a clear model structure. Body-fixed equations of motion, derived from Lagrange’s method, are contained in field elements. An alternative cage implementation and a flexible outer ring are introduced, providing a more complete basis for bearing modeling. The flexible outer ring is modeled using the modal approach with a solution for moving loads. The derivation of the equations of motion and the bond graph formulation are thoroughly presented. Following the display of versatility and connectivity in the modeling approach, a simulation example is presented to demonstrate the capability and usefulness of the model for transient analyses. A classic time-series analysis of a run-up of an unbalanced rotor, with and without a damaged bearing, directly yields interesting results; the complex interaction between model elements following the passing of the first critical speed leads to a stationary vibratory condition taking significantly long time to develop in the damaged bearing case and involves persistent slip. Keywords: rotordynamic, bond graph, flexible shaft, rolling element bearing, flexible ring, cage model

Preprint submitted to Sound and Vibration

January 21, 2020

Nomenclature

u, u U v, v V w x, y, z X, Y, Z Y, Y αi- j αR β δ δ∗

Body displacement Strain energy General translational velocities Potential energy Circular beam cross-plane width Body-fixed coordinates Inertial frame coordinates Mode shapes Contact angle in X-Y plane, body i - body j Radial distance from beam reference surface Rotation angle, cross-sectional plane Hertzian contact displacement Dimensionless Hertzian contact displacement parameter  Imbalance eccentricity  Ellipticity parameter εθ , ε0θ Circumferential strain, 0 ≡ zeroth order ζ Damping ratio η Modal participation factors κ Hertzian contact stiffness λ Modal damping coefficient µ Contact slip dynamic friction coefficient ν Poisson’s ratio ρ Density % Curvature (inverse of radius) σθ , σiθ Circumferential stress, i ≡ from centrifugal forces τ, τ Generalized forces φ, ψ, θ Body rotations about x, y, z Φ Twist (torsion) in ring or beam ω Angular frequency ω General angular velocities vector Ω Constant angular velocity about Z-axis

Circular beam undeformed reference plane radius Beam cross-sectional area Modal amplitudes Hertzian contact damping coefficient Modal amplitude ratio Bending stiffness Damping matrix Unit vector Effort variable Modulus of elasticity (Young’s modulus) Hertzian contact equivalent modulus of elasticity Elliptical integral of the second kind Curvature difference Matrix of centrifugal, Coriolis and gyroscopic coupling terms h Circular beam in-plane thickness H(δ) Heaviside function I Mass moment of inertia IA Second moment of area k Beam distributed (unit) additional stiffness K Membrane stiffness K Stiffness matrix Kθ Change of ring curvature K Elliptical integral of the first kind m Mass, modal mass M Mass and inertia matrix m Transformation matrix (transformer modulus matrix) mi- j Transformation matrix from body i to contact in direction body j N Finite number (e.g. of modes) p, p Generalized momenta q, q Generalized coordinates r Radius r Position vector R, θ, Z Transverse (in-plane radial), circumferential (inplane tangential) and axial (cross-plane) ring coordinates s Normalized position along straight beam s Hertzian contact slip T Kinetic energy T Forcing function vector

a A A, B c C D D e e0 , e0 E E0 E F(%) G

Subscripts: b Pertaining to ball(roller) c Pertaining to cage d Pertaining to disc ir Pertaining to inner ring n Pertaining to a single mode or Pertaining to outer ring s Pertaining to shaft k Variable multiple per mode n  Pertaining to imbalance

1. Introduction The purpose of this paper is to investigate model fidelity requirements for rotordynamic systems incorporating rolling element bearings. Notwithstanding that a rolling element bearing itself constitutes a rotordynamic system, order of magnitude in some cases allows for simplified representation of bearing elements. For the same reason, classic rotordynamic theory often focus on the characteristics of fluid bearings, due to their lower stiffness, higher damping and usually more significant impact on overall system dynamics. Model considerations, however, depend not only on the basic properties of the system, but also the mode of operation and the analysis objectives. Industry-standard rotordynamic models, most often implemented using the FE approach, are readily available for the analysis of eigenvalues 2

and steady-state behavior of structurally complex systems. Recent years, however, have seen an increased interest in models representing a mechatronic or system perspective, taking active rotor control, rotor transient issues and fault conditions into consideration. Moreover, application areas are emerging that expand the envelope of traditional rolling element bearing usage and challenge common modeling approaches. Examples include main bearings in wind turbines and backup bearings in magnetic bearing systems. Rotordynamics is an inherently challenging branch of mechanics, and mathematical models are instrumental in analyzing and visualizing characteristic responses of even the most basic systems. The effects of imbalance and inertial coupling were explored in the classical works of Jeffcott [1], Stodola [2] and Green [3], giving names to the Jeffcott and Stodola-Green rotordynamic models. In modern engineering, complex structures, high speeds, extreme temperatures, fluid-structure interaction and electromagnetic coupling all add to the challenge of modeling, and, today, rotordynamic research has a large publication base. The fundamentals are accessible and have been well-presented over the last three decades through the works of Vance [4], Childs [5], Lalanne and Ferraris [6], Yamamoto and Ishida [7], Friswell [8], Matsushita et al. [9] and others. Broad application-related areas such as fluid-structure interaction and electromagnetic coupling have also seen substantial literary contributions. Applying a system perspective to rotordynamics, Xu et al. [10], in a recent review of the subject of electromagnetic coupling, pointed to multiphysics approaches, system-level analyses and model coupling as a future trend. Wagner et al. [11], in a review on model reduction methods applied to rotordynamics, acknowledged that system analyses based on direct application of complex FE models can yield prohibitively heavy computations. They concluded that more attention to the inherent characteristics of rotordynamics and also large-scale industrial application examples are needed. Regarding rolling element bearings, basic theory has been treated extensively in the literature, for example in Harris and Kotzalas [12, 13]. Gupta [14] notably contributed to rigid multi-body rolling element bearing modeling. Wensing and van Nijen [15] advocated for the importance of considering the effect of structural flexibility on system dynamic behavior. Cao et al. [16] recently reviewed model development of rolling element bearing-based rotordynamic systems, highlighting the tradeoff between accurate prediction and computational burden. They concluded that more research on advanced rolling element bearing models and model coupling is needed. Multidisciplinary analyses inevitably urge modelers and analysts to seek compromises between accurate and detailed representation of a system on one hand, and computational efficiency on the other. Finding the right balance in fidelity is an important task in its own right. Consequently, the overall modeling process will benefit greatly from modularity and scalability. The bond graph method is a modeling approach visualizing energy flow between elements representing inertia, compliance, dissipation and energy conversion in a dynamic system. Among its traits are connectivity and consistent model representation across energy domains, identification of proper state variables and system causality issues as well as straightforward generation of the system’s nonlinear state-space equations. Energy-based methods such as Lagrange’s method are useful for deriving the coupled equations of motion inherent in rotordynamic systems. These coupled equations are conveniently stored in energy-conserving fields in the bond graph, connected by scalable multidimensional bonds representing the number of coordinates employed. As equations are derived for different rotating bodies or structures, the associated fields serve as reusable objects. Supporting and encouraging multidisciplinarity and object-oriented modeling, the bond graph method provides a readily customizable basis for assembling and investigating models of varying complexity. Rotordynamic modeling based on the bond graph method has only received limited attention in the literature; Campos et al. [17] presented a bond graph model of the Jeffcott rotor, focusing on the differences between Lagrange’s method and the Hamiltonian method for obtaining the coupled equations of motion for a rotating body. The equations were contained in a multiport capacitive element (C-field), embodying storage of the sum of both kinetic and potential energy. Lagrange’s method was found to be the most straightforward method for obtaining the model equations. Pedersen [18] applied bond graph formulation to the classical Jeffcott and Stodola-Green rotors, using Lagrange’s method to derive the equations of motion. Successive approaches to structuring the bond graphs were demonstrated, leading to a distinct graphical adaptation of the IC-field known from Karnopp et al. [19]. The work was extended through the implementation of a flexible shaft using the finite mode approach (Pedersen [20]) and later the inclusion of shear correction (Pedersen and Polic [21]). As with rotordynamics, rolling element bearing modeling using the bond graph method has received only limited attention in the literature. Nakhaeinejad and Bryant [22] presented a multi-body rolling element bearing model used 3

for simulating vibration response caused by bearing surface defects and compared the results with experimental data. Each rolling element was represented by a rigid body, whose translational coordinates in the inertial frame and rotations of the body-fixed coordinate system were separated onto two three-dimensional multibond 1-junctions. Inertial coupling and gyroscopic effects were implemented on the body-fixed rotations in the form of an Eulerian Junction Structure (EJS), known from Breedveld [23]. Mishra et al. [24] presented three rolling element bearing models of different complexity, also aimed at simulating vibrations caused by defects and thus offsetting the need for physical experiments. The middle alternative was a bond graph-based multi-body model, which, unlike the model in Nakhaeinejad and Bryant [22], included a cage. They concluded that the most advanced model, developed using commercial multi-body software, better replicated the experimental data. However, they also pointed out that the bond graph model has significant potential in terms of integration into larger multi-domain systems as well as implementation of more advanced contact models in the bearing itself. Considering flexible bodies, Karnopp et al. [19] pointed to the finite mode approach as the most accurate low-order representation for distributed parameter systems and showed how readily they are implemented using bond graphs. One advantage of finite mode representation is the way in which it handles moving loads, although with some special measures in place, as briefly described in Margolis [25]. Huang and Soedel [26] derived the equations of motion for a rotating ring, including strain stiffening from centrifugal forces. By use of the equations, it was shown that a rotating ring under a stationary point load does not exhibit resonant behavior in the same way as a stationary ring under a moving point load. Resonance phenomena in an elastic ring under a moving load has in itself attracted interest, for example, as discussed by Forbes and Randall [27]. In this paper, we use the bond graph method to develop a rotordynamic model incorporating rolling element bearings. The multi-body system is consistently represented by body-fixed equations of motion, derived using Lagrange’s method and formulated as graphically recognizable I- or IC- fields, depending on the coupling terms. We show how the modeling approach results in a clear and accessible bond graph structure. We also include roller cages in the bearings, connected differently than in Mishra et al. [28], demonstrating how the versatile approach readily supports alternative configurations. Both the shaft and the bearing outer ring are modeled as flexible bodies, using the finite mode approach. We take advantage of the fact that the equations of motion for a rotating ring with strain stiffening, as derived by Huang and Soedel [26], provide for an elegant, general bond graph implementation of a flexible ring with a moving load—including when the ring is not rotating. In the last part of the paper, we apply a rather synthetic approach. The model is used to simulate a rotor run-up with damage introduced to one of the bearings, which shows the effect a given condition has on the transient response of the system as a whole. 2. Rotor model 2.1. Coordinates and equations of motion Consider the rotor with a flexible shaft and associated coordinate systems shown in Fig. 1. The Jeffcott rotor is only modeled in the transverse coordinates uX and uY , and introduces us to inertial coupling effects from mass imbalance in planar motion, known as centrifugal and Coriolis terms. It provides for understanding the concepts of synchronous motion and critical speeds. The Stodola-Green rotor adds the rotational coordinates φ and ψ about the x and y axes, and, with that, the effect of gyroscopic coupling. Since motion in the X-Z and Y-Z planes uncouple for zero rotor speed, the gyroscopic effects are understandably speed dependent, as are the natural frequencies of the rotor. In the literature, the Stodola-Green rotor is seen modeled with both rigid and flexible shafts—the difference being that the bearings provide the only compliance in the former case. The work presented in this paper is based on the Stodola-Green rotor model from [20]. A rigid disc is attached to a rotating Euler-Bernoulli shaft, modeled using the finite mode approach. Shear correction is not considered here, but readily implemented as shown in Pedersen and Polic [21].

4

Y



y

i

Y

y

m ,i

x

i

Y





dz

X

 x

X

z

Z 



Z

uY

X

uX

z

Z Figure 1: Rotor assembly and coordinate systems

The kinetic energy of a rigid disc i attached to the shaft at an axial position z is straightforwardly expressed as   1  1 T d, i = md, i u˙ 2X, i + u˙ 2Y, i + I x, i ω2x, i + Iy, i ω2y, i + Iz, i ω2z, i 2 2 (1)     1   1 1 = md, i u˙ 2X, i + u˙ 2Y, i + I x, i φ˙ 2i + ψ˙ 2i + Iz, i θ˙ 2 − 2φi ψ˙ i θ˙ . 2 2 2 The kinetic energy of the flexible shaft is obtained by integrating the kinetic energy of the differential shaft element dz over the full length of the shaft 1 Ts = 2

ZL

ρA



u˙ 2X

+

u˙ 2Y



1 dz + 2

0

ZL





ρIA φ + ψ dz − 2θ˙ ˙2

˙2

0

ZL

˙ + ρIA Lθ˙ 2 . ρIA φψdz

(2)

0

The shaft coordinates are continuous functions in z. Separation of variables is employed to solve the shaft equations of motion, assuming the availability of a solution where the coordinates are sums of products of kinematically admissible mode shape functions YX (z), YY (z), and time-dependent, generalized coordinates qX (t), qY (t) uX (z, t) =

∞ X

uY (z, t) =

YXn (z)qXn (t),

n=1

φ(z, t) =

YYn (z)qYn (t),

n=1

∞ X dYXn (z) n=1

∞ X

dz

ψ(z, t) = −

qXn (t),

∞ X dYYn (z) n=1

dz

(3) qYn (t).

The shaft sees no bending moment or shear force at z = 0 and z = L, that is, the second- and third-order partial derivatives of uX and uY with respect to z equal zero in these positions, and thus Y 00 (0) = Y 00 (L) = Y 000 (0) = Y 000 (L) = 0. The force-free boundary conditions give a general solution to the spatial equation for lateral bending of a straight Euler-Bernoulli beam, identical for the individual mode shapes in each of the two orthogonal planes Yn (s) = sin kn Ls + sinh kn Ls − Cn (cos kn Ls + cosh kn Ls) ,

Cn =

sin kn L − sinh kn L cos kn L − cosh kn L

(4)

where s = z/L is the normalized position along the shaft length L, and k˜ n is given by ρA 2 k˜ n4 = ω, EIA n 5

(5)

and obtained by solving the free-free beam model equation. Selecting a finite number of modes N, including the translational and rotational rigid body contributions in the X- and Y- directions, gives two coordinate arrays qX (t) and qY (t), and the mode shape array Y(s) qX (t) =

h

qX00

qX0

qX1

···

qXN

qY (t) =

h

qY00

qY0

qY1

···

qY N

Y(s) =

h

1

iT iT

L(s − 1/2) Y1 (s) · · ·

,

,

(6)

YN (s)

i

.

By combination of Eq. (6) and Eq. (3), and insertion into the shaft kinetic energy equation (Eq. [2]), the N + 2 by N + 2 inertia matrices M s1 , M s2 and M s3 for the first three terms of Eq. (2) can be found through integration 1 1 1 T q˙ X M s1 q˙ X + q˙ TY M s1 q˙ Y = q˙ TX ρAL 2 2 2

Z1

1 Y (s)Y(s)ds q˙ X + q˙ TY ρAL 2

0

1 T 1 1 q˙ M s2 q˙ X + q˙ TY M s2 q˙ Y = q˙ TX ρIA L 2 X 2 2

Z1

T

0

Z1

Y0T (s)Y0 (s) 1 ds q˙ X + q˙ TY ρIA L 2 L2

0

˙ TX M s3 q˙ Y = −2θq ˙ TX ρIA L −2θq

YT (s)Y(s)ds q˙ Y , Z1

Y0T (s)Y0 (s) ds q˙ Y , L2

(7)

0

Z1

Y0T (s)Y0 (s) ds q˙ Y . L2

0

Similarly, the N by N stiffness matrix K associated with the shaft flexible modes can be obtained from the bending strain energy expression U=

1 T 1 q K1 qX + qTY K1 qY 2 X 2

1 = qTX EIA L 2

Z1

Y00 T (s)Y00 (s) 1 ds qX + qTY EIA L 2 L4

0

Z1

Y00 T (s)Y00 (s) ds qY . L4

(8)

0

Revisiting the disc kinetic energy equation (Eq. [1]), it can now also be seen that inertial matrices for each of the first three equation terms, for a total number of discrete discs ND, are obtained as the respective inertial parameters md, i , I xd, i and Iyd, i multiplied by sums of the mode shapes evaluated at the individual positions zi for each disc i Md1 =

ND X

md, i YT (zi )Y(zi ),

i=1

Md2 =

ND 1 X I xd, i Y0 T (zi )Y0 (zi ), L2 i=1

Md3 =

ND 1 X Iyd, i Y0 T (zi )Y0 (zi ). L2 i=1

(9)

2.2. Imbalance excitation The kinetic energy of a single imbalance, represented by a point mass m,i , located at a distance i from the center of the shaft with an orientation ∆θi in the shaft body-fixed coordinate system, can be expressed as h i  1 ˙ i u˙ X cos θ˜i − u˙ Y sin θ˜i m, i u˙ 2X + u˙ 2Y + i2 θ˙ 2 + 2θ 2   1 = θ˙ 2 m, i i2 + m, i i2 θ˙ Y(zi ) q˙ X cos θ˜i − Y(zi ) q˙ Y sin θ˜i , 2

T , i =

where θ˜i = θ + ∆θi , and state variables are evaluated at the position of each individual imbalance along the shaft.

6

(10)

Neglecting the contribution from imbalance mass to the total rotor inertia, the forcing function resulting from the centrifugal forces from a number N of discrete imbalance masses can be expressed as        T X   −θ˙ 2 a˜ sin θ + b˜ cos θ   =     , T =  T Y   −θ˙ 2 a˜ cos θ − b˜ sin θ  (11) N N X X     m, i i2 YT (zi ) sin ∆θi . m, i i2 YT (zi ) cos ∆θi , b˜ = a˜ = i=1

i=1

2.3. Bond graph representation and the IC-field Energy methods are used to derive the coupled equations of motion for IC-field representation as introduced by Karnopp et al. [19]. Given that the inertia states will be represented by generalized momenta, it shall be considered here that they are already contained in Lagrange’s equations as the partial derivative of the kinetic energy with respect to the generalized velocities   M1+2 ∂T  p= =  0 ∂q˙  2qTY MT3 M1+2 = M s1 + M s2 + Md1 + Md2 ,

0 M1+2 0

2M3 qY 0 Iz, tot

   q˙ X     q˙ Y  ˙ θ

M3 = M s3 + Md3 ,

    ,  Iz, tot = Izs +

(12) ND X

Izd, i .

i=1

The second contribution to Lagrange’s equation is the partial derivative of the kinetic energy with respect to generalized displacements, obtained as       q˙ X  0    ∂T = θ˙  −2MT3   q˙ Y  , (13)   ˙  ∂q 0 θ Lastly, the term containing the potential energy U from the flexible shaft modes becomes   K1 ∂U  =  0 ∂q  0

0 K1 0

0 0 0

   qX     qY  θ

    . 

(14)

As is evident, the equations of motion are only coupled in the inertial terms. Thus, an I-field representation is appropriate, as shown in Fig. 2. Shown in the same figure, the transformer (TF) element for connecting a bearing in position z has the modulus    Y(z) 0 0    m(z) =  0 (15) Y(z) 0    0 0 1 The choices made so far prompt a comment on a practical aspect of the resulting bond graph; extending the equations of motion to include additional coordinates will not alter the graphical structure, just increase the number of velocities on the multibonds. 3. Rolling element bearing model 3.1. Coordinate systems and transformations Only planar bearing motion is considered in this paper. The bearing inner ring simply represents the undeformed circular section of the flexible shaft at a given axial position. The equations of motion for the rollers, cage and outer 7

I

C

Se



1

1

TF

u X  Velocities, bearing   u connection point  Y   

Shaft modal  q X ( z )  velocities   qY ( z )    Rotor spin   

TF

1 u X   uY  Velocities, bearing   connection point  

R

Figure 2: Rotor bond graph

ring are described in body-fixed coordinates and derived using Lagrange’s method. The rollers and cage are treated as rigid bodies in general motion, and their equations of motion are coupled in the inertial terms if any imbalance is present. The overall bearing bond graph is essentially a 3 degree of freedom multibond structure representing the ˙ T. body-fixed velocities [ x˙ y˙ θ] A coordinate system is aligned with the line of Hertzian contact between two circular-section bodies i and j. It is rotated by an angle αi- j given by the relative positions of the origins of the body-fixed coordinate systems in the inertial frame, that is ! Y j − Yi αi- j = atan2 . (16) X j − Xi Thus, the multibonds at the contact between the two bodies reduce to two degrees of freedom in the planar case. Fig. 3 gives an overview of the employed coordinate systems. The transformation of forces and displacements to and from the contact between the two circular-section bodies i and j is modulated by the difference between contact coordinate system angle αi- j and the angular position of the adjoining body-fixed coordinate systems, θi or θ j . Except for the flexible outer ring introduced in Sec. 4, the defining body geometries are considered undeformed by forces. For the rollers and the inner ring, this means that any contact force acts at a constant radius ri from the body center. The roller side transformation of a contact towards the inner ring will, for example, take the form:   cos (αb-ir − θb ) − sin (αb-ir − θb )  mb-ir =  sin (αb-ir − θb ) cos (αb-ir − θb )  0 rb

    . 

(17)

A similar transformation is applicable for contact between neighboring rollers in a full complement bearing. 3.2. Why and how to include a cage A rolling element bearing model, without other roller supports than the Hertzian contacts towards the inner and outer rings, can operate effectively only as long as slip does not occur. The equidistant spacing between rollers, as set by the initial conditions, will be disturbed at the first loss of roller traction. Thus, additional contacts between the cage and the rollers, or directly between the rollers in a full complement bearing, are required if situations involving slip are considered. Moreover, cage kinematics should arguably be considered to the same extent as roller kinematics. Assigning mass and inertia to a cage will, for example, enable the capturing of cage run-out effects. There are several ways to include a cage. One way is to connect the rollers to a cage body by orthogonal spring damper elements in the inertial frame, as shown by Mishra et al. [28]. Another option is to use dedicated coordinate systems rotated to the lines of contact between the cage and the rollers, or directly between neighboring rollers. The latter choice facilitates any heterogeneous connection in the radial and circumferential directions, and provides freedom in representing clearance and pocket geometry. For example, cage pocket sides need not be taken as parallel, and can be inclined to provide Hertzian contact force components in the radial direction. With only one contact point 8

Contact rotation gives FR component in θ-direction and vice versa

Outer ring

Contact rotation from ring curvature change

Cage

yb-or

yb-c rc

xb-ir

Y yc yor

Deformed reference plane

xb

xb-c

yb-ir

Inner ring

xb-or

yb

yir

θb

Undeformed reference plane

xb-c

θc θir

αor-b

yb-c

rc

a

αb-c αb-ir

xc

X

X

xor xir Figure 3: Bearing coordinate systems

on each side of the pocket, however, the risk of the cage pinching the roller should be considered. Depending on the actual pocket geometry, an alternative is to include multiple cage-roller contacts. Here, a cage with one Hertzian contact facing each side of the rollers is implemented. The cage pockets are taken as having planar, parallel sides. Positive clearance is given by the contact displacement, initial value δ(0). The model maintains linear damping with constant coefficient for the Hertzian contact—regardless of clearance. However, a different damping model is easily implemented, for example, reducing or switching off damping under positive clearance. To limit cage run-out and the corresponding change in spacing between the rollers, the cage is also modeled as being guided by the inner ring. An orthogonal pair of linear spring damper elements connects the cage and the inner ring in the inertial frame. The roller side transformation of the roller-cage contact is  
π  cos ϕ − sin ϕ  ϕ = αb-c − − θb − αb-c − θc − θc∗   2 mb-c =  sin ϕ cos ϕ  , (18) π   ∗ + θ − θb , = θ c− c 0 ±rb 2 where the bracketed quantity on the right side of Eq. (18) is the contribution to the transformation angle from roller movement within the cage pocket clearance. Locked to a fixed pocket center angle θc∗ in the cage coordinate system, αc-b disappears from the roller side transformation, and the relative rotation between the roller and the cage is determined only by the respective body rotations θb and θc . The cage side transformation of a cage-roller contact is       q  cos θc∗ − π − sin θc∗ − π 
2 2 (Xb − Xc )2 + (Yb − Yc )2 , rcx = cos αb-c − θc − θc∗       π π ∗ ∗   mc-b =  sin θc − 2 (19) q cos θc − 2  ,   ∗
2 2 r = sin α − θ − θ (X − X ) + (Y − Y ) , cy b c c b c b c c −r ±r cx

cy

9

where the individual pocket rotation angles θc∗ are constant and the torque arm components rcx and rcy are variable in the cage body-fixed coordinate system. 3.3. Hertzian contact model A rolling element bearing is essentially characterized as a non-linear spring. As described by Hertzian theory, the restoring forces from elastic deformation between, and in direction normal to, two surfaces more or less singly or doubly curved, follows the law F = κδn , where κ is the stiffness and δ is the contact displacement or mutual approach of the contacting bodies. The exponent n depends on the type and geometry of the rolling element bearing, and is usually taken as n = 3/2 for a regular deep groove ball bearing. In this paper the focus is solely on solid mechanics. It is, however, important to note that omitting the effects of operation in an Elasto-Hydrodynamic Lubrication (EHL) regime represents a simplification with respect to restoring forces. EHL will both increase the stiffness, due to viscosity effects in the lubricant, and alter the effective contact angle. The situation is different with respect to damping. Material hysteresis in the elastic surface deformation is negligible compared to the viscous damping from EHL. The damping force in direction normal to the surfaces can be ˙ showing the fact that the damping coefficient c is itself a function of displacement. The damping written as F = c(δ)δ, coefficient is obtained on an empirical basis, and a constant value for c is used here. Based on the works of Wijnant et al. [29] and Wensing and van Nijen [15], values from 25 to 100 Ns/m have been considered for the simulations presented in Sec 5. Harris and Kotzalas [12] presents the following formula for obtaining the contact stiffness κ: s 2 2E 0 P , (20) κ= % 3δ∗ 3/2 where % is curvature, defined as the inverse of a radius, and X

%=

P

% is the curvature sum

1 1 1 1 + + + ri1 ri2 r j1 r j2

,

(21)

where indexes i and j refer to body i and j, and indexes 1 and 2 refer to orthogonal section planes 1 and 2 normal to the contact surfaces. E 0 is the equivalent modulus of elasticity (Young’s modulus) for the contact between two bodies i and j  2 −1  1 − νi2 1 − ν j  0 E =  + (22)  , Ei Ej  and δ∗ is a dimensionless displacement found tabulated or graphed against an auxiliary quantity F(%) called curvature difference ! ! 1 1 1 1 − + − ri1 ri2 r j1 r j2 P F(%) = . (23) % δ∗ is defined as

2K  π 1/3 , (24) π 2 2 E where  is the ellipticity ratio and K and E are the elliptical integrals of the first and second kind. Wensing [30] demystifies the normalized quantity δ∗ by reminding us that  = 1 and K = E = π/2 for purely circular contact profiles. Thus, δ∗ vanishes from Eq. (20). This assumption simplifies the task of approximating the stiffness for any bearing for which the exact geometry is not known, although at the expense of a slightly lowered stiffness. The friction model for slip s at the contact is based on Kragelsky et al. [31], and the friction coefficient µ is expressed as n o µ = (−0.1 + 22.28|s|) e(−181.46|s|) + 0.1 sign (s) (25) δ∗ =

The friction coefficient comprises stick-slip, Coulomb and Stribeck characteristics. Viscous friction is not present, which is consistent with a dry contact situation disregarding EHL. 10

R

IC

1 q 0

Next roller

q n 

T

T

T

OUTER RING



MTF C R R

C 1

MTF

uR u 

1

 s 

 or−b−b

   s 

 c − b

T

 xb  y   b b 

 c−b , rc−b , c

ROLLER

R

1

MTF

 s 

T

0

MTF

   , r ,   s  c−b c−b c

0

MTF

T

 xc yc  c  b

MTF

mg:Se

ir

1

I

CAGE T

c

1

Global

Next roller

1

Se:mg

0 1

 or −b

MTF

1

MTF

αi-j , ri-j

MTF  or−b , ror−b MTF

INNER RING u X uY  

R

 c − b

b−ir−ir

Next roller

1

1

I

R

C

 or −b

C

1

b−ir−b

Sf : 

MTF

0

MTF

MTF

0

T

1

C R

Inertial frame connections Figure 4: Bearing bond graph

The complete contact model is expressed as     F   κδ3/2 H(δ) + c(δ)δ˙ x   =   F =    F   µ κδ3/2 H(δ) + c(δ)δ˙ y

   ˙ ,  = f (q, q)

   δ    q =   ,  0 

   δ˙    q˙ =   .  s 

(26)

3.4. Complete bearing model Fig. 4 shows the complete bearing bond graph, including the inner and the outer ring, the cage, and, for clarity, a single roller. Each body is represented by a multibond 1-junction, representing the body-fixed velocities. The forces from the Hertzian contact laws are summed in local 1-junctions and applied to the 0-junctions representing the relative velocities between the bodies. The signal arrows indicate that the state of the restoring C-element modulates the friction in the R-element. Modulated transformers (MTFs) provide the required coordinate transformations between the body-fixed velocities and all contacts. There are two potential contact paths between the roller and the cage, since contact on either side of the pocket requires its own unique transformation. The inner ring 1-junction is highlighted 11

 or −b+

Flexible ring angle orrections β and γ

C R R

C b−c−c

MTF

1

OUTER RING

MTF

uR u 

T

1

 s 

 or−b +  −b

T

0

MTF

I

b−c−b

   s 

R

C

ROLLER

1

  b−c−c  s 

CAGE

MTF

0

I

Roller 1  xb yb b  velocities

b−ir−b

0.015

T

0.01

MTF

0.005

C

0

1 R

 s 

0

MSf

MTF

b−ir−ir

T

Inner ring T velocities 1 u X uY ir 

-0.005 -0.01 -0.015 - /32

D

(a) Damage flow source on bond graph

- /64 (

b_ir

-

ir

-

* ) D

/64

/32

(b) Superimposed velocity profile

Figure 5: Introducing damage to inner ring

with a white dot, indicating that it is identical to the velocities of a shaft connection point in Fig. 2. The dashed rectangle in the lower right part of the figure encompasses the transformations to the inertial frame necessary for calculating the angles αi- j and radii ri- j , and for including the linear cage to inner ring contact. The outer ring contact 0-junction sees an additional bond representing a velocity correction from the deflection of the ring and given by a transformation of the ring’s states. The details of the outer ring part of the bond graph are explained in Sec. 4.6 and Sec. 4.7. Since applying damage to the bearing is considered a special case, the inclusion in the bond graph is shown separately in Fig. 5, to not make Fig. 4 more busy than necessary. 3.5. Introducing bearing damage The local coordinate systems at the body contacts provide for straightforward implementation of a Hertzian contact velocity disturbance as caused by geometrical imperfections of the bodies. The modulated flow source (MSf) as shown ˙ representing a passing spall cavity—in in Fig. 5 superimposes a slightly trapezoidal velocity profile on the velocity δ, this case on the contact between the roller and the inner ring. Belonging to the inner ring, the damage is located at a fixed angle θ∗D in the inner ring body-fixed coordinate system, and its passing of each roller is timed by the angle αb ir − θir − θ∗D . In this simplified method, the contact coordinate system is not rotated to follow the actual geometrical disturbance, and tangential components from the Hertzian contact force are ignored. The only change in tangential forcing is from any resulting roller slip. 4. Flexible bearing outer ring model 4.1. Coordinates and assumptions The target is a finite mode representation for a vibrating, closed circular beam. Considering the relatively thin rings, Euler-Bernoulli theory is employed, assuming that plane sections remain plane and normal to the neutral axis after deformation. Neglecting transverse shear deformation leaves only circumferential strain εθ to be considered. For the same reason, it is ignored that differential bending strain in a curved beam is not proportional to the distance from the neutral axis, and that the neutral axis does not coincide with the centroid axis, as shown in Fig. 6(a). Seidel and Erdelyi [32] refer to a ratio of 1/10 or lower between ring thickness and ring neutral axis radius as rendering the Euler-Bernoulli assumption valid. For cases requiring more careful consideration, Rao and Sundararajan [33] provide insight on the gradual effects of omitting rotational inertia and shear deformation in rings. Rao [34] describes how a total of ten equations in ten unknowns represents the three-dimensional flexural vibrations of a circular ring. If the undistorted neutral axis of the ring is plane, that is, of single curvature, the equations can be shown to fall into two sets—each expressible by an internally coupled displacement variable pair. As shown 12

Y

uR

uθ Neutral axis θ

h

kR

a

θ+β

P

uZ

Reference axis, e.g. centroid dA

M

Φ

w

O'

θ

M αR

Z

(a) Bending strain

kθ

Ω

X

(b) Coordinate system Figure 6: Circular beam

in Fig. 6(b), these are the in-plane displacement pair uR and uθ , and the pair made up by the out-of-plane (axial) displacement uZ and the cross section torsional (twist) displacement Φ. Only in-plane motion is considered in this paper, and uR and uθ are referred to as transverse and circumferential displacement, respectively. The equations of motion are developed for a circumferential strip with unit cross-plane width (Z-direction). Thus, the restraining effect expressed by Poisson’s . .h ratio ν is i neglected, and, as a simplification 2 3 of general plate- and shell theory, the bending stiffness, D = Eh 12(1 − ν ) , reduces to D = Eh3 12, and the . membrane stiffness, K = Eh (1 − ν2 ), reduces to K = Eh. Rao and Sundararajan [33] give guidance on deriving the equations of motion for rings with boundary conditions imposed in discrete locations. In this paper, the ring is assumed to be free, although, as shown later in Eq. (41), provisions will be made for including an additional, distributed stiffness in transverse and circumferential directions. Moreover, the ring is assumed to have a constant cross section around the circle and a neutral axis with a constant radius. The Fourier series is used for the mode shapes. Wensing and van Nijen [15] also use the Fourier series for in-plane vibration, but use Chebyshev polynomials for out-of-plane vibration. 4.2. Equations of motion for a vibrating, thin ring Soedel [35] shows that the equations of motion in both the transverse and circumferential coordinates uR and uθ can be written as ! ! ! D ∂4 uR K ∂3 uθ ∂uθ ∂2 uR − + 2 uR + + ρh = τR , a4 ∂θ4 a ∂θ3 ∂θ ∂t2 ! ! ! (27) D ∂3 uR K ∂uR ∂2 uθ ∂2 uθ ∂2 uθ − 2 + + ρh = τθ . − a4 ∂θ3 a ∂θ2 ∂t2 ∂θ2 ∂θ As is generally known from Euler-Bernoulli theory, longitudinal vibrations in a beam typically have natural frequencies one order of magnitude higher than transverse vibrations. It is not uncommon to assume the ring neutral axis to  be inextensible, and, by introducing the relation uR = − ∂uθ ∂θ, decouple the two equations of motion, obtaining a single equation in one or the other coordinate. The choice is a matter of the relative importance of circumferential forcing and circumferential inertial terms. Both coordinates are retained throughout the model development in this paper. It is, however, shown later under Sec. 4.8 how inextensional theory can be used to reduce the forced solution.

13

As in Eq. (3), separation of variables is used, assuming the availability of solutions in the following form uR (θ, t) =

∞ X

YRn (θ) qRn (t) ,

uθ (θ, t) =

n=0

∞ X

Yθn (θ) qθn (t).

(28)

n=0

Substituting the time-dependent functions qRn (t) and qθn (t) with eiωn t , and setting the forcing to zero, the equations of motion can be written as ! ! D ∂4 YRn K ∂3 Yθn ∂Yθn − + 2 YRn + + ρhω2n YRn = 0, a4 a ∂θ4 ∂θ3 ∂θ ! ! (29) D ∂3 YRn K ∂YRn ∂2 Yθn ∂2 Yθn 2 − − 2 + + ρhωn Yθn = 0. a4 a ∂θ3 ∂θ2 ∂θ ∂θ2 Two intuitive aspects simplify the consideration of boundary conditions and mode shapes for the force-free ring. First, since the ring is continuous, and θ and θ + 2π are the same physical point (sharing all physical quantities), the spatial solution must be periodic. Thus, a possible solution is YRn (θ) = An ei n(θ−ϕ) ,

Yθn (θ) = Bn ei n(θ−ϕ) .

(30)

The angle ϕ is appropriately included in the functions to illustrate that the force free ring has no preference for the orientation of its modes. Secondly, YRn and Yθn must be orthogonal; this is illustrated in Fig. 7(b), which shows the n = 2 mode, or first bending mode, for ϕ = 0. As point P moves outward on the horizontal axis, such that YR2 (0) = A2 and Yθ2 (0) = 0, point R moves inward on the vertical axis such that YR2 (π/2) = −A2 and Yθ2 (π/2) = 0. Point Q moves mainly circumferentially instead of transversely, such that YR2 (π/4) = θ and Yθ2 (π/4) = −B2 . By the insertion of Eq. 30 in Eq. 29, the equations of motion can be written as n o K  n o  D 4 3 n A − i n B + i nB + A + ρhω2n An = 0, n n n n a4 a2 o  Kn  D n 3 o 2 2 i n A − n B − i nA + n B + ρhω2n Bn = 0, n n n n a4 a2

(31)

and arranged in matrix form    K  n4 D + K + ρhω2   A  3D i −n + n n    n  4 2 4 2 a a a a     (32)     = 0.       2D 2K 2   −i −n3 D + n K  Bn −n 4 − n 2 + ρhωn a4 a2 a a   Since the trivial solution An Bn T = 0 is not interesting, the determinant of the left matrix is required to be zero. This yields a fourth order frequency equation, ω4n + a2 ω2n + a0 = 0,

(33)

which has two roots in each n, here named ωnk , k = 1, 2. As previously stated, ωn2  ωn1 . The amplitude ratio between the circumferential and transverse displacements can be expressed as a function of mode number, Cnk = Ank /Bnk = f (n). Ank and Bnk are, as such, arbitrary, since only the ratio Cnk defines the mode shapes. Looking at lower n numbers, it can be shown that 1 Cn2 ≈ − . n

Cn1 ≈ −n,

(34)

This is an important result, as it means that uR dominates at the k = 1 frequencies and uθ dominates at the k = 2 frequencies. For n = 0, ω01 disappears, analogously to rigid body rotation, and ω02 represents the ’breathing mode’ of the ring (see Fig. 7[a]). For n = 1, ω11 = 0; thus, bending still does not occur, analogously to rigid body translation. 14

R

Q P

‘Breathing Mode’

n=0

n=2

n=1

(a) Ring n = 0 and n = 1 modes, none of which include bending

(b) n = 2 mode, first bending mode

Figure 7: Ring modes

4.3. Ring rotation and strain stiffening Including the stiffening effect from centrifugal forces in a rotating ring begins with modifying the kinetic energy and strain energy terms. As shown in Ginsberg [36], the velocity vector of a point P on a body can be expressed as
v = v00 + vP/00 rel + ω × rP/00 , (35)
where v00 is the velocity vector of the origin O0 of a body-fixed coordinate system, vP/00 rel is the velocity vector of the point P relative to the body-fixed coordinate system, ω is the angular velocity vector of the body fixed coordinatesystem, and rP/00 is the position vector of the point P in the body-fixed coordinate system. By revisiting Fig. 6(b), it can be seen that a point on the reference surface of an undistorted flexible ring at radius a from the ring center, given that the ring rotates about the z-axis with velocity Ω, has velocity vectors and position vector v00 = 0eR + aΩeθ + 0eZ ,
vP/00 rel = u˙ R eR + u˙ θ eθ + 0eZ ,

ω = 0eR + 0eθ + ΩeZ ,

(36)

rP/00 = uR eR + uθ eθ + 0eZ .

The total velocity vector is v = (˙uR − uθ Ω) eR + (˙uθ + uR Ω + aΩ) eθ + 0eZ ,

(37)

As shown in Huang and Soedel [26], the kinetic energy can now be expressed as ρha T= 2

Z2π h i (˙uR − uθ Ω)2 + (˙uθ + uR Ω + aΩ)2 dθ.

(38)

0

Considering the strain-displacement relations for the ring, these are εθ =

ε0θ

1 + αR Kθ + 2a2

 !2 !2    u u ∂ ∂ R θ  uR + + − uθ  , ∂θ ∂θ

ε0θ =

1 a

! ∂uθ + uR , ∂θ

Kθ =

1 ∂β 1 = a ∂θ a2

∂uθ ∂2 uR − ∂θ ∂θ2

! (39)

ε0θ represents the membrane strain and αR Kθ the bending strain. Kθ is the change of curvature from bending, β the corresponding rotation angle, and αR the distance from the reference surface at which bending strain is considered. The boxed term in Eq. (39) is a non-linear term, which has to be included to consider the stiffening effect from tension due to the centrifugal forces, which in turn will contribute to the restoring mechanism at high angular velocities. A guitar string serves as a relevant yet extreme analogy. The strain energy can now be expressed as U=a

Zh/2 Z2π "

# 1 σθ εθ + σiθ εθ dθ, 2

−h/2 0

where σiθ is the unit stress from the centrifugal forces, and can be approximated as ρa2 Ω2 . 15

(40)

4.4. Additional, distributed stiffness As shown in Huang and Soedel [26], an additional restoring mechanism, not directly related to the physical properties of the ring itself, is straightforwardly introduced in the form of a distributed stiffness acting in the transverse and/or circumferential directions. The principle is shown in Fig. 6(b), where the additional stiffness is illustrated as a pair of orthogonal springs with unit constants kR and kθ acting on the differential ring element dθ. Thus, the remaining contribution to the total energy can be expressed as Z2π "    # 1 2 2 kR uR + kθ uθ − τR uR + τθ uθ dθ, V=a 2

(41)

0

where τR and τθ are the generalized unit forces in transverse and circumferential direction, respectively. 4.5. Modified equations of motionR Using Hamilton’s principle, δ (U + V − T )dt = 0, the coupled equations of motion are now obtained as D a4

! ! " !# σiθ h K 1 ∂4 uR ∂3 uθ ∂uθ ∂uθ ∂2 uR − − + 2 uR + + 1+ uR + 2 + kR uR + ρh a a a ∂θ4 ∂θ3 ∂θ ∂θ ∂θ2

∂2 uR ∂uθ − 2Ω − Ω2 uR − aΩ2 ∂t2 ∂t

!

∂ uθ ∂uR + 2Ω − Ω2 uθ ∂t2 ∂t

!

= τR (42)

D a4

K ∂ uR ∂ uθ − − 2 a ∂θ3 ∂θ2 3

2

!

∂uR ∂ uθ + + ∂θ ∂θ2 2

!

σiθ h a2

∂uR ∂ uθ − ∂θ ∂θ2 2

uθ − 2

! + kθ uθ + ρh

2

= τθ ,

The boxed terms in Eq. (42) are the result of including ring rotation. Note that contributions from the straindisplacement relations in Eq. (39) to anything but tension from centrifugal forces is neglected, meaning that the equations of motion are effectively linearized about an operating point given by the rotational speed of the ring. By substituting σiθ with ρa2 Ω2 in the boxed terms, and redoing the steps in Eq. (28) to Eq. (31), the equations of motion can be written in matrix form as     K K     4D 2 2 2 3D 2 n + − n ρhΩ + k + ρhω i −n + n + 2nρhΩ + 2ρhΩω R n n 4 2 4 2    An  a a a a     (43)     = 0        D K D K   −i −n3 + n + 2nρhΩ2 + 2ρhΩωn −n2 4 − 2 − n2 ρhΩ2 + kθ + ρhω2n  Bn a4 a2 a a Requiring that the determinant of the left matrix is zero, a new fourth order frequency equation is obtained, ω4n + a2 ω2n + a1 ωn + a0 = 0,

(44)

with coefficients a2 , a1 and a0 in n and Ω. The presence of the first-order term a1 ωn is the result of including ring rotation, and leads to that solving for ωn now gives four distinct frequencies for each mode, beginning from n = 1, except for Ω = 0. The additional two frequencies are effectively bifurcations due to Coriolis acceleration. The bifurcation frequencies are unaffected by the support stiffness represented by kR and kθ , despite their contribution to the overall stiffness. The solutions to the mode shapes can be written as YRnk (θ, t) = Ank ei (nθ+ωnk t) ,

Yθnk (θ, t) = Bnk ei (nθ+ωnk t) ,

(45)

where n = 0, ±1, ±2, ..., k = 1, 2, 3, 4, and the amplitude ratio between the circumferential and transverse displacements is now a function of frequencies and angular velocity, Cnk = i (Bnk /Ank ) = f (ωnk , Ω). The time-dependent nature of these mode shapes, as opposed to the mode shapes in Eq. (30), is related to the way in which they rotate around the ring, as they are neither stationary nor synchronous with the imposed rotational speed Ω. As in Eq. (34), it can be shown for the amplitude ratio that Cnk ≈ −n,

1 Cnk ≈ − , n

k = 1, 2, 16

k = 3, 4.

(46)

4.6. Forced equations and bond graph implementation For continuous systems with an infinite number of degrees of freedom, solutions can be represented as complex Fourier series, uR (θ, t) =

∞ X 4 X

ηnk (t) cos (nθ + ωnk t) ,

uθ (θ, t) =

∞ X 4 X

n=0 k=1

Cnk ηnk (t) sin (nθ + ωnk t) .

(47)

n=0 k=1

By collecting the sin(nθ) and cos(nθ) terms, Eq. (47) can be recast to express the displacements uR and uθ as functions of a set of four coupled generalized modal coordinates per mode, uR (θ, t) =

∞ X   q1n (t) cos (nθ) + q2n (t) sin (nθ) ,

uθ (θ, t) =

n=0

∞ X 

 q3n (t) cos (nθ) + q4n (t) sin (nθ) ,

(48)

n=0

where the generalized modal coordinates are given as        qn =      

q1n (t) q2n (t) q3n (t) q4n (t)

  4   P ηnk (t) cos (nθ + ωnk t)     k=1 4   P   − ηnk (t) sin (nθ + ωnk t)   k=1  =  P   4 C η (t) cos (nθ + ω t) nk nk nk     k=1 4   P Cnk ηnk (t) sin (nθ + ωnk t)

        .     

(49)

k=1

Conversely, Eq. (48) is inserted into the forced equation of motion, Eq. (42), to express the generalized modal coordinates qn as a function of the forces τR and τθ ∞ X

∞ X

[Rn1 (t) cos (nθ) + Rn2 (t) sin (nθ)] = τR (θ, t),

[Rn3 (t) cos (nθ) + Rn4 (t) sin (nθ)] = τθ (θ, t),

(50)

n=0

n=0

where the functions R1n , ..., R4n are given as R1n = ρh (q¨ 1n − 2Ωq˙ 3n ) + c˜ q1n + n˜aq4n ,

˜ 3n − n˜aq2n , R3n = ρh (q¨ 3n + 2Ωq˙ 1n ) + bq

R2n = ρh (q¨ 2n − 2Ωq˙ 4n ) + c˜ q2n − n˜aq3n ,

˜ 4n + n˜aq1n , R4n = ρh (q¨ 4n + 2Ωq˙ 2n ) + bq

(51)

˜ c˜ are shorthands given as and a˜ , b, a˜ = n2

D K + + 2ρhΩ2 , a4 a2

D K  b˜ = n2 4 + 2 + ρhΩ2 + kθ , a a

c˜ = n4

D K + + ρhn2 Ω2 + kR . a4 a2

(52)

A set of linear second-order differential equations in the generalized modal coordinates can now be assembled; mq¨ n + mGq˙ n + Kn qn + Dn q˙ n = T n (θ, t) ,   0  0 G =   2Ω 0

0 0 0 2Ω

−2Ω 0 0 0

 0   −2Ω  , 0  0

  λn1  0 Dn =   0 0

0 λn2 0 0

0 0 λn3 0

0 0 0 λn4

    ,  

  c˜  0 Kn =   0 n˜a

0 c˜ −n˜a 0

0 −n˜a b˜ 0

n˜a 0 0 b˜

    ,  

(53)

where m is the modal mass, equal to the circumferentially distributed mass ρh per unit width of the ring, independently of the mode number. G is the skew-symmetric matrix of inertial coupling terms, also independent of mode number. Dn is the diagonal modal damping matrix, with the frequency-dependent modal damping coefficient given as λn = 2ρhζn ωn . Kn is the modal stiffness matrix. T n is the modal forcing function. 17

By introducing generalized modal momenta, pn = mq˙ n ,

(54)

the four first order differential equations per mode are obtained as e0 n (q˙ n ,qn )

z }| { d p˙ n = m q˙ n = mq¨ n = −mGq˙ − Kn qn −Dn q˙ n + T n (θ, t) . dt

(55)

To make a complete set of eight first order differential equations, the remaining four are simply taken as q˙ n =

pn . m

(56)

Each set of modal velocities is linearly coupled in the non-diagonal inertial and stiffness matrices in Eq. (53) and u =  Yn ( ) q ( t ) can be represented by an IC field, as shown in Fig. 8(a). The diagonal damping matrix in Eq. (53) is appropriately n Yn ( ) q ( t ) u= n included as a separate multibond R-element.

 

τ τ u u

R R

Yn ( ) Yn ( )

n

n

TF TF

0 0

en (q n , q n ) en (q n , q n ) pn pn

1 1

qn qn

IC IC

 

IC IC

 

(a) Constant evaluation coordinate θ

 

τ τ u u

R R

Yn ( (t ) ) Yn ( (t ) )

n n

MTF1 MTF1

0 0

en (q n , q n ) en (q n , q n ) pn pn

1 1

qn qn

Yn ( (t ) ) q n ( t ) MTF2 Yn ( (t ) ) q n ( t ) MTF2

1 1 

qn qn

+ − + −



+ +

+ +

* *

(b) Time-dependent evaluation coordinate θ(t)

u =  Yn ( (t ) ) q ( t ) +  Yn ( )  ( t ) q ( t  u= n Yn ( (t ) ) q ( t ) +  n Yn ( )  ( t ) q ( t

Figure 8: IC-field for flexible outer ring

n

1 Tn = ν

Z2π

( Yn (θ) · γ · τ (t) dθ,

n , 0, ν = π n = 0, ν = 2π

n

) ,

0

 0 0 0  cos (nθ)  (nθ) 0 − sin 0 0 Yn (θ) =  0 0 cos (nθ) 0  0 0 0 − sin (nθ)

    , 

  1  1 γ =   0 0

0 0 1 1

    , 

(57) " τ (t) =

τR τθ

# ,

Eq. (57) shows how the modal forcing T n couples to the generalized transverse and circumferential forces τR and τθ , through the mode shape function Yn (θ). γ is a connectivity matrix. The Hertzian contact forces FR and Fθ are applied here as point loads, defined by use of the Dirac delta function 
as τ[R,θ] = F[R,θ] (t) a δ (θ). Since the equations of motion are developed for a ring of unit width, the forces FR and 18

Fθ must be divided by the ring cross-plane width w. The modal forcing as function of the Hertzian contact forces becomes TF z }| { " # ( ) 1 FR /w n , 0, ν = π Tn = · Y (θ) · γ · , , (58) Fθ /w n = 0, ν = 2π aν defining the modulus of the transformer TF in Fig. 8. 4.7. Moving load Analogously to Eq. (28), the relation between the generalized velocities u˙ and the generalized modal displacements q˙ is given as TF z }| { ( ) ∞ X 1 n , 0, ν = π T T u˙ (θ, t) = · γ · Yn (θ) ·q˙ n (t) , , (59) n = 0, ν = 2π aν n=0 where γ is the same connectivity matrix as in Eq. (57). In our model, the evaluation coordinate θ of the modal shape functions Yn (θ) is not constant, but changes with the roller-to-outer-ring contact. A time-dependent θ leads to the ˙ inclusion of a second term in the expression for displacements u; MTF1

MTF2

z }| { z }| { ∞ ∞ X X 1 1 T T T 0T · γ · Yn (θ(t)) ·q˙ n (t) + · γ · Yn (θ(t)) · qn (t) ·θ˙ (t) , u˙ (θ, t) = aν aν n=0 n=0

(

n , 0, ν = π n = 0, ν = 2π

) .

(60)

This can be explained as follows: The changing ring curvature rotates the contact angle, introducing a circumferential component of the transverse forcing and vice versa, as shown in Fig. 3. The rate of change of θ(t) represents a flow that multiplies with these new force components to produce power. This is easily overlooked for forcing with zero mean value, but will otherwise manifest as an integrating power imbalance that will cause displacements uR and uθ to drift off. Margolis [25] describes how the energy imbalance caused by a time-dependent evaluation coordinate can be compensated for by the inclusion of an ancillary flow source Sf in the bond graph. Here, however, the rate of change of θ is the result of the system response to other external excitation, and the reciprocal modal force components to θ˙ are transformed, summed, and applied to the roller. θ˙ corresponds to the velocity α˙ b−or , accessible through a transformation of roller position in the inertial frame from Cartesian to polar coordinates, and not through derivation of Eq. (16). Fig. 8(b) shows the bond graph for the outer ring IC-field with a time-dependent evaluation coordinate. θ∗ is the roller to outer ring contact initial position, which is identical to the initial value of the global angle αb-or . 4.8. Inextensional deformation The inextensible neutral axis assumption mentioned under Sec. 4.2 is briefly revisited. Huang and Soedel [26] propose a method for neglecting factors associated with extensional modes in the equations for a rotating ring. Eq. (46) and Eq. (49) give us the following approximate relations 1 ηn4 (t) ≈ − ηn1 (t), n Consequently, Eq. (48) can be rewritten as ∞   X uR (θ, t) = q1n (t) cos (nθ) + q2n (t) sin (nθ) ,

ηn3 (t) ≈

uθ (θ, t) =

n=0

∞  X 1 n=0

1 ηn2 (t). n

 1 q2n (t) sin (nθ) − q1n (t) cos (nθ) , n n

(61)

(62)

reducing the number of generalized modal coordinates from four to two. The matrices for Eq. (53) then becomes         ˜ c − n2 a˜ 2    b˜ b˜ + n2 a˜         0 0 −2Ω        ˜  ˜   λ 0 n a ˜ + b a ˜ + b n         Gn =  , Dn =  (63) , Kn =       ,      ˜ c − n2 a˜ 2    b˜ + n2 a˜ b˜ 0 λn    2Ω        0 0    n a˜ + b˜ a˜ + b˜ 19

and the forcing function takes the form      T n =    



         . Tn2 b˜ + n˜aTn3      a˜ + b˜ Tn1 b˜ − n˜aTn4   a˜ + b˜

(64)

These equations are only coupled in the Gn matrix, and can be represented in a bond graph by a two-dimensional multibond I-field, together with separate C and R elements. 4.9. Stationary ring forced equation Lastly under the bond graph implementation, Eq. (27) is revisited, and a model development ignoring ring rotation in the first place is shown. The forced solution can then be written as uR (θ, t) =

∞ X 2 X 

 ηn1k (t) cos (nθ) + ηn2k (t) sin (nθ) ,

n=0 k=1

(65)

∞ X 2 X   uθ (θ, t) = Cnk ηn1k (t) sin (nθ) − Cnk ηn2k (t) cos (nθ) , n=0 k=1

involving four modal coordinates h iT h iT qn = q1n (t) q2n (t) q3n (t) q4n (t) = ηn11 (t) ηn12 (t) ηn21 (t) ηn22 (t) .

(66)

The amplitude ratio can be expressed as Cnk =

Bnk Ank

  n   D  2 n + K a2 a2 = ! . 2  D n 2 +K ρhωnk − 2 a a2

The coefficient matrices for Eq. (53) now becomes   0 0   λn1 0   0 λn1 0 0  , Dn =  0 λn2 0   0 0 0 0 λn2

 2  ωn1  0 Kn =   0 0

(67)

0 ω2n1 0 0

0 0 ω2n2 0

0 0 0 ω2n2

    ,  

(68)

with no matrix of Coriolis terms or centripetal terms present. Thus, the modal equations are internally decoupled and can be represented by ordinary four-dimensional multibond I-, C-, and R-elements. The resulting forcing function is slightly more involved, with FR and Fθ connecting through two different modulated transformer elements (MTFs)     0 0 0  cos(nθ)   1/Nn1      1  0 sin(nθ) 0 0   1/Nn1  FR T nn = ·  ·   1/Nn2  · w + 0 0 cos(nθ) 0 aν      0 0 0 sin(nθ) 1/Nn2  0 0 0  Cn1 sin (nθ)  1  0 −Cn1 cos (nθ) 0 0 ·  0 0 Cn2 sin(nθ) 0 aν   0 0 0 −Cn2 cos(nθ) (

n , 0, ν = π n = 0, ν = 2π

)

    1/Nn1     1/Nn1  ·  1/N n2   1/Nn2

  2 Nnk = π 1 + Cnk .

, 20

    Fθ  · ,  w

(69)

In terms of bond graph implementation, it can be seen by comparison that the modified equations of motion in Sec. 4.6 in fact yields an overall elegant solution. 5. Simulation The model is implemented in 20-sim and simulated using the parameters shown in Tab. 1. The presented simulation results focus on the effects of bearing damage on system behavior during run-up. The rotor is accelerated from standstill to rated speed in approximately one second, during which time it passes through the first critical speed, and kept at a constant speed thereafter. A mass imbalance is present in the same axial position as the disc. Two cases are compared: one with undamaged bearings and one with a spall cavity introduced to the inner ring in one of the bearings. All state variables are set to zero before the start of the simulation, except the Hertzian contact displacements. The roller-cage contact displacements are set to an initial value equal to half the total roller-pocket clearance, and the roller-inner ring and roller-outer ring contact displacements are set to correspond to a preload of 10 N in each contact. The preload is applied to provide some margin towards slip, which still will occur due to dynamic loads in the damaged bearing case. The contact preload also effectively gives an initial membrane stress in the outer ring. All external forces are onset at the start of the simulation. The rotor is accelerated by the application of a torque with proportional control, giving a final speed of 3585 rpm in the damaged bearing case versus a rated speed of 3600 rpm.

Bearing inner ring

Disc (mid shaft)

5.1. Rotor behavior

Figure 9: Lateral displacement, x-direction, disc (top) and inner ring (bottom), comparing cases with undamaged bearings and one damaged bearing

Fig. 9 shows the lateral displacements of the disc and the bearing inner ring for both the damaged and undamaged bearing cases. The disc and the bearing inner ring lateral displacements peak during the quick passing through the first critical speed and decays exponentially in the undamaged bearing case. In the damaged bearing case, the lateral displacement of the disc tapers off faster and in a seemingly more linear fashion. During this period, the displacement of the inner ring remains approximately at the level reached during the passing of the critical speed. After the displacement of the disc from the passing of the critical speed subsides, the damaged bearing sees a steadily increasing vibration level, stabilizing first after approximately 25 seconds. 21

Young’s modulus Poisson’s ratio Density

210 x 109 Pa 0.25 7810 kg/m3

Rotor: Shaft diameter Shaft length Shaft flexible modes Shaft rigid modes Modal damping coefficient Bearing positions on shaft Disc diameter Disc thickness Disc position on shaft Imbalance mass Imbalance radius Imbalance phase angle Imbalance position on shaft

15 mm 600 mm 5 2 0.1 10 mm, 590 mm 200 mm 10 mm 300 mm 50 g 15 mm 0◦ 300 mm

Bearing: Type Outer diameter Width No of balls Ball diameter Inner race contact diameter Outer race contact diameter Inner race groove radius Outer race groove radius Cage mass Cage moment of inerta

6202 35 mm 11 mm 8 5.953 mm 19.347 mm 31.253 mm 3.0465 mm 3.2165 mm 14 g 2.24 x 10−6 kgm2

Cage pocket clearance

0.2 mm

Cage guiding on inner ring: Translational stiffness Translational damping Rotational damping

1 x 106 N/m 24 Ns/m 5 x 10−4 Nms

Hertzian contact parameters: Ball preload Stiffness inner race Stiffness outer race Stiffness cage Hertzian contact damping

10 N 5.096 x 109 N/m3/2 6.252 x 109 N/m3/2 2.751 x 109 N/m3/2 36 Ns/m

Bearing surface damage: Depth Width Damage phase angle

0.25 mm 3.0 mm 50◦

Bearing outer ring: Thickness Reference plane radius Angular velocity Additional circumferential stiffness Additional transverse stiffness Modal damping coefficient Flexible modes (incl. bending) Included parts of non-bending n=0 and n=1 order modes

1.8735 mm 16.56325 mm 0 s−1 0 N/m 0 N/m 0.02 4 n=0 transverse, ’breathing mode’

Table 1: Model parameters

Fig. 10 shows the magnitude of the bearing force vector together with the force vector rotational speed and whirl speed normalized to rotor speed. During the passing of the critical speed, the rotational frequency of the bearing force vector drops from following the shaft rotational frequency and converges towards the first rotor natural frequency, which is approximately 33.6 Hz. This speed ratio persists as the lateral force from natural vibration decays, until the minimum value of the superposed force fluctuation from the rotor imbalance hits zero. From that point, the rotational speed of the force vector again starts gaining on the rotor speed, while the force vector fluctuation continues to settle around the value determined by the rotor imbalance. The time from passing critical speed to the force hitting zero is more than halved in the damaged bearing case. After hitting zero, the force vector, averaging around the imbalance force, is very much influenced by the impulses from bearing damage, quickly eliminating the trace of natural vibration decay. The convergence of the force vector rotational speed after the transition point also happens faster in the damaged bearing case. The whirling pattern of the bearing is influenced by the rate at which the inner ring passes the rollers. Otherwise similar in behavior to the force vector rotational speed, the bearing whirl speed is super-synchronous—close to two times the rotor speed. 5.2. Cage behavior Fig. 11 shows the orbits of the inner ring and cage mass centers through two full shaft rotations at four different times for the damaged bearing case. Scatter plots are used, with marker color intensity fading in a backwards direction. 22

200

2.0

150

1.5

100

1.0

50

0.5

0

5

10

15

20

25

30

Figure 10: Bearing rotating force vector - force magnitude and force rotational speed relative to rotor speed (with and without damage).

The density of markers gives an indication of speed variation. In Fig. 11(a), at one second, just after reaching full speed, the inner ring orbit shows a marked polygonal pattern, as well as speed variation from passing of the rollers. Here, the large rotating force from rotor resonance causes the inner ring to wedge the rollers apart as it orbits. The between five and six vertices correspond to the total number of rollers minus roller position shift (cage rotation) in one shaft rotation. The super-synchronous orbit speed of just above 1.5 times the rotor speed is evident in that the inner ring orbit completes a little more than three revolutions per two shaft rotations. As the inner ring forces the rollers apart, the cage pocket geometry forces the cage in the radially opposite direction from the inner ring, resulting in a smaller orbit radius for the cage. The bearing sees little slip at this point. In Fig. 11(b), at 5.4 seconds, just after the force vector hits zero, the inner ring and cage orbital patterns are irregular and in the process of changing. Interestingly, the inner ring clearly traces out four complete revolutions, of which one collapses. The cage vaguely traces out a stationary triangular pattern, switching to an inverted threepointed star pattern corresponding to the collapsing orbit of the inner ring. This coincides with periodic slipping events reaching the highest peak slip velocities seen through the simulation run, as also shown in Sec. 5.3. In Fig. 11(c), at ten seconds, a larger plot scale is used. The change in contact force from the damage now causes the inner ring to orbit at a larger radius. The inner ring traces out four circles, producing a basic 2X deformation pattern. However, due to the damage, the inner ring orbit also causes a 1X cyclic variation in the roller passing, which is slightly more pronounced than the 2X pattern at this point. The combined motion is more visible in the cage orbit, now exhibiting a wavelike pattern. Upon closer examination, the orbit can be seen to gradually assume a rotating, triangular shape. Only moderate slip occurs at this point. As the run proceeds, relative motion between cage and inner ring increases until it reaches a peak at around 20 seconds. In Fig. 11(d), at 25.3 seconds, plotted at an even larger scale, the inner ring orbital pattern has developed into a triangular shape in which it remains steady. Cyclic bearing deformation is now synchronous at 1X. The train has broken free and is slipping steadily, rotating at 0.5X instead of the geometrical 0.38X. The cage orbit now follows the inner ring orbit, with only moderate relative motion between the two.

23

Inner ring Cage

(a) @1.0s

(b) @5.4s

(c) @10.0s

(d) @25.3s

Figure 11: Damaged bearing - inner ring and cage orbits through two full inner ring rotations at different times. Note increasing radial scale from left to right.

5.3. Flexible outer ring behavior The rigid body motion of the bearing outer rings is not considered in these simulations. Of the n = 0 and n = 1 modes (see Fig. 7[a]), only the n = 0 contribution to transverse displacement is included. Sometimes called the ’breathing mode’, it has a natural frequency of 50019 Hz in the given example. The lowest ring natural frequency is 4375 Hz, belonging to the n = 2 first bending mode. The added stiffness as described in Sec. 4.4 is set to zero, leaving only the inertial and restoring forces of the outer ring to balance the forces from the Hertzian contacts. A somewhat unrealistic case for a 6202 type bearing specifically, this will exaggerate the outer ring deflection in the simulation. Moreover, the ring is taken as not spinning, that is, Ω is zero. Figs. 12 to 15 illustrate the displacements of the inner ring and roller mass centers, and the outer ring deflection at different points during the 30-second run for the damaged bearing case. Each figure captures a full rotation of the shaft, pictured through nine frames, (a) to (i), taken at equidistant times. Fading shadows trace the changes from one frame to the next. The figures are augmented by multiplying outer ring transverse displacement from base radius a by a factor of 100, and positioning the inner ring and roller mass centers at 67 times the actual radial distance from mean radial positions. The figures also show the angular positions of the rotor imbalance, the inner ring damage and a tracer roller. The inner ring rotating force vector and the roller to outer ring Hertzian contact force vectors are indicated by red arrows, having the same scale in all figures. At one second in Fig. 12, there appears a singly symmetric oval or egg-shaped outer ring deformation rotating at the inner ring whirling frequency of just above 1.5 times rotor speed. The largest deformation radius follows the inner ring orbit, as the inner ring forces the rollers evenly against the outer ring. At 5.4 seconds in Fig. 13, the symmetry is less pronounced. In the last frame of Fig. 13(i), the ring assumes close to a circular shape during a slipping event as described in Sec 5.2. The trailing shadows of the rollers show how the rotational speed of the train increases as the rollers break free from the outer ring. However, at ten seconds in Fig. 14, the outer ring again assumes a symmetric deformation shape. Not affected by the forces from the lateral vibration of the rotor, the deformation is less egg-shaped than at 1 second, if not fully doubly symmetric. Both 1X and 2X deformation are present. At 25.3 seconds in Fig. 15, the fully developed pattern shows a clear 1X behavior, with the vertices of the triangular orbit shown in Fig. 11, corresponding to a pronounced egg-shaped deformation. The deformation is not rotating but 24

(a) @0

(b) @π/4

(c) @π/2

(d) @3π/4

(e) @π

Disc imbalance angle Bearing damage angle

Force vector Tracer roller

(f) @5π/4

(g) @3π/2

(h) @7π/4

(i) @2π

(j) Legend

Figure 12: Damaged bearing - One full shaft rotation @1.0s

(a) @0

(b) @π/4

(c) @π/2

(d) @3π/4

(e) @π

Disc imbalance angle Bearing damage angle

Force vector Tracer roller

(f) @5π/4

(g) @3π/2

(h) @7π/4

(i) @2π

(j) Legend

Figure 13: Damaged bearing - One full shaft rotation @5.4s

rather reshaping itself in the different directions, alternating between the egg-shape and a more circular shape. 6. Conclusion In summary, we have modeled a rotordynamic system incorporating rolling element bearings using the bond graph method. Our model repository has been expanded with an alternative cage support and a flexible outer ring, forming a more complete basis and set of building blocks for representing different bearing configurations. We have also used the model to simulate a run-up of an unbalanced rotor, with and without a damaged bearing. The results show that the lateral displacement at the bearing inner ring or shaft end reaches the same amplitude when passing the rotor’s first critical speed in both the undamaged and damaged bearing cases. In the damaged bearing, however, the amplitude does not decay but remains virtually constant, while energy from the rotor’s normal mode vibration is dissipated. As 25

(a) @0

(b) @π/4

(c) @π/2

(d) @3π/4

(e) @π

Disc imbalance angle Bearing damage angle

Force vector Tracer roller

(f) @5π/4

(g) @3π/2

(h) @7π/4

(i) @2π

(j) Legend

Figure 14: Damaged bearing - One full shaft rotation @10.0s

(a) @0

(b) @π/4

(c) @π/2

(d) @3π/4

(e) @π

Disc imbalance angle Bearing damage angle

Force vector Tracer roller

(f) @5π/4

(g) @3π/2

(h) @7π/4

(i) @2π

(j) Legend

Figure 15: Damaged bearing - One full shaft rotation @25.3s

the force from the rotor approaches that of the imbalance alone, the bearing sees an increase in slip, occurring in periodic bursts with high slip velocities. Then follows a transitional period with gradually increasing amplitude until a severe, steady inner ring or shaft end vibration pattern is developed where the bearing train slips continuously. The lateral vibration amplitude in the disc position is dominated by the imbalance and the first bending mode, and the most apparent effect from the bearing damage on the rotor is the increased rate of decay after passing the first critical speed, due to the increased bearing damping. Throughout the modeling process, we have demonstrated the versatility of the chosen approach by elaborating on the development of the equations of motion and the structuring of the bond graph. The storage of body-fixed equations of motion in field elements provides for an object-oriented approach and reuse of elements. The multibonds’ dimensions conveniently size with the number of coordinates employed per body. The bond graph structure with transformers provides for changing and multiplying mechanical bodies and the connections between them, and for 26

interfacing other objects and domains, and brings clarity to the overall model. From the simulation, we have learnt that the system transient response is not readily intuitive or inferred from the characteristics and parameters of the individual components and chosen constraints. Particularly striking is the complicated sequence playing out in the damaged bearing, and the amount of time it takes for a stationary condition to develop. Although affected by the low bearing stiffness due to the unsupported outer ring in this example, the total system behavior demonstrates some of the capability and usefulness that we want to achieve with the model. Besides this qualitative aspect, any potential use of the results for the given physical rotor-bearing system in continuous operation at a moderate constant speed is arguably limited to producing a vibration spectrum and determining whether loads are so severe that bearing deterioration will accelerate. Further time-series analysis would, however, be required if our study were to be expanded to include operation involving time-varying speed or load, irregular loads from, for example, rubbing, or even time-varying external forces to the bearing outer ring, such as from a squeeze film damper arrangement—all of which are possible uses for the model. Acknowledgments Funding: This work was supported by Equinor ASA; and the Research Council of Norway [project number 263819]. [1] H. H. Jeffcott, XXVII. The lateral vibration of loaded shafts in the neighbourhood of a whirling speed.–The effect of want of balance, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 37 (219) (1919) 304–314, doi:\bibinfo{doi}{10.1080/ 14786440308635889}. [2] A. Stodola, Dampf- und Gasturbinen: mit einem Anhang u¨ ber die Aussichten der W¨armekraftmaschinen [Steam and Gas Turbines: with an Appendix on the Future of Heat Engines], Springer, Berlin, 6. aufl. edn., 1924. [3] R. B. R. Green, Gyroscopic effects on the critical speeds of flexible rotors, Journal of Applied Mechanics–Transactions of The ASME 15 (4) (1948) 369–376. [4] J. M. Vance, Rotordynamics of turbomachinery, Wiley, New York, ISBN 0471802581, 1988. [5] D. Childs, Turbomachinery rotordynamics : phenomena, modeling, and analysis, Wiley, New York, ISBN 047153840X, 1993. [6] M. Lalanne, G. Ferraris, Rotordynamics prediction in engineering, Wiley, Chichester, 2nd ed. edn., ISBN 0471972886, 1998. [7] T. Yamamoto, Y. Ishida, Linear and nonlinear rotordynamics : a modern treatment with applications, Wiley, New York, ISBN 0471181757, 2001. [8] M. I. Friswell, Dynamics of rotating machines, vol. 28, Cambridge University Press, Cambridge, ISBN 9780521850162, 2010. [9] O. Matsushita, M. Tanaka, H. Kanki, M. Kobayashi, P. Keogh, Vibrations of Rotating Machinery : Volume 1. Basic Rotordynamics: Introduction to Practical Vibration Analysis, vol. 16, Springer Japan : Imprint: Springer, Tokyo, ISBN 4431554564, doi:\bibinfo{doi}{10.1007/ 978-4-431-55456-1}, 2017. [10] X. Xu, Q. Han, F. Chu, Review of Electromagnetic Vibration in Electrical Machines, Energies 11 (7) (2018) 1779, doi:\bibinfo{doi}{10.3390/ en11071779}. [11] M. B. Wagner, A. Younan, P. Allaire, R. Cogill, Model Reduction Methods for Rotor Dynamic Analysis: A Survey and Review, International Journal of Rotating Machinery 2010 (273716), doi:\bibinfo{doi}{10.1155/2010/273716}. [12] T. A. Harris, M. N. Kotzalas, Rolling bearing analysis : [1] : Essential concepts of bearing technology, vol. [1], Taylor & Francis, Boca Raton, 5th ed. edn., ISBN 084937183X, 2007. [13] T. A. Harris, M. N. Kotzalas, Advanced concepts of bearing technology, CRC/Taylor & Francis, Boca Raton, Fla, 5th ed. edn., ISBN 0849371821, 2007. [14] P. K. Gupta, Advanced Dynamics of Rolling Elements, Springer New York, New York, NY, ISBN 9781461297673, doi:\bibinfo{doi}{10. 1007/978-1-4612-5276-4}, 1984. [15] J. A. Wensing, G. C. van Nijen, The dynamic behaviour of a system that includes a rolling bearing, Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 215 (6) (2001) 509–518, doi:\bibinfo{doi}{10.1243/1350650011543763}. [16] H. Cao, L. Niu, S. Xi, X. Chen, Mechanical model development of rolling bearing-rotor systems: A review, Mechanical Systems and Signal Processing 102 (2018) 37–58, doi:\bibinfo{doi}{10.1016/j.ymssp.2017.09.023}. [17] J. Campos, M. Crawford, R. Longoria, Rotordynamic modeling using bond graphs: modeling the Jeffcott rotor, Magnetics, IEEE Transactions on 41 (1) (2005) 274–280, doi:\bibinfo{doi}{10.1109/TMAG.2004.838924}. [18] E. Pedersen, Rotordynamics and bond graphs: basic models, Mathematical and Computer Modelling of Dynamical Systems 15 (4) (2009) 337–352, ISSN 13873954, doi:\bibinfo{doi}{10.1080/13873950903063116}. [19] D. C. Karnopp, D. L. Margolis, R. C. Rosenberg, System dynamics : a unified approach, Wiley, New York, 2nd ed. edn., ISBN 0471621714, 1990. [20] E. Pedersen, A bond graph approach to rotor dynamics, in: Proceedings of the 2012 - 10th International Conference on Bond Graph Modeling and Simulation, ICBGM’12, Part of SummerSim 2012 Multiconference, vol. 44, ISBN 9781618399854, 44–49, 2012. [21] E. Pedersen, D. Polic, Bond graph modeling of rotordynamic systems with a flexible shaft including shear correction, in: ICBGM’2014, vol. 46, The Society for Modeling and Simulation International, ISBN 07359276, ISSN 07359276, 205–212, 2014. [22] M. Nakhaeinejad, M. D. Bryant, Dynamic modeling of rolling element bearings with surface contact defects using bond graphs, Journal of Tribology 133 (1), doi:\bibinfo{doi}{10.1115/1.4003088}.

27

[23] P. C. Breedveld, Proposition For An Unambiguous Vector Bond Graph Notation, J Dyn Syst Meas Control Trans ASME V 104 (3) (1982) 267–270, ISSN 15289028, doi:\bibinfo{doi}{10.1115/1.3139707}. [24] C. Mishra, A. K. Samantaray, G. Chakraborty, Ball bearing defect models: A study of simulated and experimental fault signatures, Journal of Sound and Vibration 400 (2017) 86–112, doi:\bibinfo{doi}{10.1016/j.jsv.2017.04.010}. [25] D. Margolis, Finite mode bond graph representation of vehicle-guideway interaction problems, Journal of the Franklin Institute 302 (1) (1976) 1–17, doi:\bibinfo{doi}{10.1016/0016-0032(76)90099-5}. [26] S. C. Huang, W. Soedel, Effects of coriolis acceleration on the free and forced in-plane vibrations of rotating rings on elastic foundation, Journal of Sound and Vibration 115 (2) (1987) 253–274, ISSN 10958568, doi:\bibinfo{doi}{10.1016/0022-460X(87)90471-8}. [27] G. Forbes, R. Randall, Resonance phenomena of an elastic ring under a moving load, Journal of Sound and Vibration 318 (4–5) (2008) 991–1004, ISSN 0022460X, doi:\bibinfo{doi}{10.1016/J.JSV.2008.05.021}. [28] C. Mishra, A. K. Samantaray, G. Chakraborty, Bond graph modeling and experimental verification of a novel scheme for fault diagnosis of rolling element bearings in special operating conditions, Journal of Sound and Vibration 377 (2016) 302–330, doi:\bibinfo{doi}{10.1016/j. jsv.2016.05.021}. [29] Y. H. Wijnant, J. A. Wensing, G. C. Nijen, The Influence of Lubrication on the Dynamic Behaviour of Ball Bearings, Journal of Sound and Vibration 222 (4) (1999) 579–596, doi:\bibinfo{doi}{10.1006/jsvi.1998.2068}. [30] J. A. Wensing, On the dynamics of ball bearings, Ph.D. thesis, University of Twente, Enschede, 1998. [31] I. V. Kragelsky, M. N. lrm;Dobychin, V. S. lrm;Kombalov, Friction and Wear: Calculation Methods, Elsevier Science & Technology, ISBN 0080273203, 1982. [32] B. S. Seidel, E. A. Erdelyi, On the Vibration of a Thick Ring in Its Own Plane, Journal of Engineering for Industry 86 (3) (1964) 240, doi:\bibinfo{doi}{10.1115/1.3670524}. [33] S. S. Rao, V. Sundararajan, In-Plane Flexural Vibrations of Circular Rings, Journal of Applied Mechanics 36 (3) (1969) 620, doi:\bibinfo{doi} {10.1115/1.3564726}. [34] S. S. Rao, Vibration of continuous systems, John Wiley & Sons, Hoboken, N.J., ISBN 0471771716, 2007. [35] W. Soedel, Vibrations of shells and plates, vol. 177, Marcel Dekker, New York, 3rd ed., r edn., ISBN 0824756290, 2004. [36] J. H. Ginsberg, Advanced engineering dynamics, Cambridge University Press, Cambridge, 2nd ed. edn., ISBN 0521470218, 1995.

28