Full Dynamic Ball Bearing Model with Elastic Outer Ring for High Speed Applications

Ball bearings are commonly used in high speed turbomachinery and have a critical influence on the rotordynamic behavior. Therefore, a simulation model of the bearing to predict the dynamic influence is essential. The presented model is a further step to develop an accurate and efficient characterization of the ball bearing’s rotor dynamic parameters such as stiffness and deflections as well as vibrational excitations induced by the discrete rolling elements. To make it applicable to high speed turbomachinery, the model considers centrifugal forces, gyroscopic effects and ball spinning. The consideration of an elastic outer ring makes the bearing model suitable for integrated lightweight bearing constructions used in modern aircraft turbines. In order to include transient rotordynamic behavior, the model is built as a full dynamic multibody simulation with time integration. To investigate the influence of the elasticity of the outer ring, a comparison with a rigid formulation for several rotational speeds and loads is presented.


Introduction
Widely used in common machinery, rolling element bearings permit a rotational motion between two components.The advantages of ball bearings like low friction and maintenance make them suitable for most applications.In contrast to journal bearings, no rotor instability occurs, which makes them perfect for high speed applications.The trend towards more powerful and lightweight turbomachinery requires precise control and prediction of rotor vibrations and dynamic behavior.A precise bearing model is essential for the following two reasons: first, it allows the calculation of typical rotordynamic coefficients like speed and load-dependent stiffness.These parameters are important for the determination of the rotor's critical speeds as well as for the definition of the operational ranges.Second, a precise bearing model provides insight into the bearings excitation frequency.For instance, under radial load, there is always an excitation in the range of the ball passing frequency due to the discrete rolling elements.A damaged bearing with faults at the raceways or the balls generates additional system excitations in frequency ranges, which can be calculated in the dynamic simulation.In high speed lightweight applications such as turbopumps or compressors, angular contact ball bearings are used.They are axially preloaded to avoid clearance.The associated ball kinematics cause, besides the centrifugal forces, gyroscopic moments and ball spinning.This causes additional loads and wear that have to be considered in machinery design.In lightweight turbomachinery, mostly for highly integrated aircraft applications, the bearing raceways are directly integrated in the shaft and the housing itself.This hinders the direct prediction of the bearing effects on the rotordynamics.To follow this aspect in this contribution, the outer bearing ring is modeled elastically using the finite element method (FEM).This allows the determination of the interaction between the balls and the elastic structure resulting in overall shaft support stiffness.
The issue of bearing modeling has been dealt with for over 50 years.Among the first who dealt deeper with kinematic and kinetic analyzes were Jones [1] and Harris [2].Due to the lack of computing power, no extensive dynamic simulations could be carried out, which led to a quasi-static model with many simplifications.Firstly, the geometric relationships and quantities within the bearing were outlined by Harris.Secondly, contact forces and deformations (see Figure 1) were derived based on the Hertzian-contact theory.
Using the simplification stated by the race control hypothesis, which assumes that there is no relative motion between the ball and the bearing rings surface ( ω s = 0 and ∆v = 0) and the inclusion of centrifugal forces and gyroscopic moments M g , equilibrium relationships were derived.Further simplifications are made by neglecting the ball-cage interaction and assuming identical angular distribution of the balls around the rotation axis of the bearing.The model of Koch [3] and later Tüllmann [4] is an extension of Harris' study.The aim of this work was the investigation of increased maximum speed for axially loaded angular contact ball bearings.They extended the modeling by a degree of rotational freedom for the cage, rotating around the center of the outer ring.
To compute the angular velocity of the balls ω b , they first determine the angle α b by geometric relations, assuming pure rolling on both raceways.The different contact angles between the ball and the raceway (α i and α o ) cause a relative rotational movement between the surfaces (spinning motion) ω s .Using a friction model of the elastohydrodynamic (EHD) theory, a spinning friction loss results as a function of the ball rotating angle α b .By minimizing the power loss, the required rotating angle is calculated.The gyroscopic moment can then be determined by the ball rotation ω b and the angle α b .This is taken into account as a tangential force F t at the contact point in the force equilibrium of the inner and outer ring (Figure 1).
Another approach to investigate the dynamic bearing behavior is made by Oest [5] and later Fritz [6].They describe the rotational speed of the balls v b and the cage in the assembly purely from geometrical considerations in the static, unloaded position.This is possible assuming an equal ball contact angle α i = α o on the inner and the outer ring.Simplifications are also made by neglecting the gyroscopic moments, the centrifugal forces, and the drilling friction forces of the balls.These assumptions are made under the condition of pure rolling at low operation speed.
A very detailed model for rolling bearings is the ADORE program developed by Gupta [7] over decades.It describes the behavior of the bearing through a quasi-static as well as a dynamic model.With the static model, starting solutions for the dynamic simulation are determined with assumptions similar to Koch's (no cage interaction, calculation of the ball rotation via a race control assumption under kinematic relations).The dynamic model does not require any kinematic relationships or geometric assumptions and describes both roller skidding and skewing as well as the lubricating film behavior.For this purpose, six degrees of freedom (DOF) per bearing component are used, which are coupled through the forces.
Another model is the BEAST program developed by Stacke and Fritzson [8] for the company SKF (Gothenburg, Sweden).As the model of Gupta, each bearing component has six degrees of freedom, but the contact forces between the cage and the rolling elements are described in more detail.This allows for computing additional effects numerically, such as the power loss, the lubrication film behavior as well as the forces on the cage and its behavior.
To study the bearings' vibrational behavior, Sassi [9] presents a 1D model limited to 3 DOF (both rings and one ball) to study the dynamic response of a localized defect in the bearing.Another more detailed approach is made by Tkachuk [10], who uses a 3D dynamic model with 4 DOF for each part of the bearing.They analyze the effect of axial load and misalignment on the vibration signal.
Tadina [11] uses a bearing simulation model with finite element housing to simulate the vibrational response of the bearing to different local faults, modeled as ellipsoids on the races.Modeling the raceway faults or imperfections with sinusoidal functions is proposed in [12].
To determine the effects of a varying ball number and centrifugal forces on the dynamic response and excitation of the rotor system, Vakharia [13] developed a simulation model that includes specific geometric constraint assumptions for the cage speed and the ball positions.
The description of the lubricant's behavior, the friction model as well as the EHD damping effects can require time-consuming computations.However, Wijnant [14] presents an alternative approximation method for lubricated contacts in order to gain computing efficiency and keep an acceptable level of accuracy for the EHD contact loads and damping coefficients [15].For instance, in the frame of a rotor dynamic study, the outer ring as well as the housing and the shaft can be modeled as flexible bodies using finite element techniques including a modal reduction [14,15].
Another method consists of a multi-level model, which interpolates the solution from a set of precomputed values to provide faster results for EHD computations (see [8,16]).
The motivation for the authors' work is to develop a fast dynamic simulation model, which includes all relevant effects for high speed applications.Therefore, a minimal set of degrees of freedom, analytical and physical motivated solutions are used to avoid the most kinematic and kinetic simplifications.Besides the classic rigid bearing model, the elasticity of the outer ring is included.
The structure of this paper is as follows: in Section 2, the kinematics of the bearing components and the dominating forces are described.Section 3 outlines the modeling of the elastic outer bearing ring including a reduction of the full finite element model of the ring.The full dynamic model is summarized in Section 4 and the time integration scheme is outlined.Section 5 gives a case study of a bearing simulation and analyzes the influence of the elasticity of the outer ring on the rotordynamic coefficients and the vibrational behavior.In Section 6, conclusions are drawn.

Bearing Kinematics and Forces
The presented angular contact ball bearing model is based on a multibody simulation, using 6 DOF for the inner ring, 2 DOF for each ball (their axial and radial displacements) and 250 DOF for the elastic outer ring.The coordinate frames used to describe the bearing configuration are shown in Figure 2. The origin of the inertial coordinate system R I is fixed at the center of the outer ring, so that r o = 0.In order to describe the balls' positions, the local coordinate frame (e ax , e rad ) is used.Each ball's local system is established according to the angle Φ b,k and determined by the cage rotation ω c .
To describe the ball-inner race and ball-outer race contacts and interactions, local coordinate frames (R k,i ) and (R k,o ) are used (see Figure 3).e ax and e rad represent the balls' degrees of freedom.

Deflections, Normal Forces and Cage Speed
The main step to describe the bearing behavior is to calculate the contact deflections δ i and δ o as well as the contact angles α i and α o between the ball and the inner and outer raceway (Figure 3).Therefore, a detailed description of the inner bearing geometry and geometric projections of the degrees of freedom are used.In the general bearing configuration, the contact angles are different between inner and outer ring because of the acting centrifugal forces.This is taken into account by using the 2 DOF of each ball and can be calculated with the dynamic equilibrium.
Assuming a pure rolling condition, meaning no slip in circumferential direction, the velocities of the ball centers can be calculated (see Figure 4): and the ball rotation is given by with the ball contact positions r bo and r bi , respectively, and the velocity parts in circumferential direction v cd bo and v cd bi .Since each ball can be positioned differently in the rings, depending on the contact geometry and the ring positions, the velocity of the ball contact point can differ in addition to the velocities v b of the center of the balls.However, the cage constrains them all to a prescribed velocity determined by the rotational velocity ω c of the cage.To calculate the resulting rotational velocity ω c , different methods exist.Here, in this contribution, a method based on a weighting approach is used and presented in the following.
The first step is to calculate the ball normal Hertzian contact forces F n,i and F n,o using the deflections δ i and δ o (see, for instance, [6]): with the Hertzian stiffness coefficient K H (see [5]).Note that, in this present study, the effect of lubrication is neglected.Its influence may require further investigations and could be considered, for instance, by a contact force law as derived by Wijnant [14].
The improved method to calculate the cage speed ω c is to take a weighted mean of the ball center velocities with the normal contact forces as weighting factors (see Figure 4): with using k for the ball number.This weighting approach ensures that the velocities of balls, which are highly loaded, are taken more into account than the velocities of balls, which are not loaded.Therefore, the higher the force F n , the more a ball is considered to satisfy the pure rolling condition.If a ball has no contact, an unloaded clearance zone occurs and no forces act on the cage.As the gyroscopic moments counterbalancing tangential forces are acting perpendicular to the cage rotation speed, they are not taken into account in this weighted averaging.
Having the cage speed ω c , the centrifugal forces F c acting on the balls are calculated with the ball mass m b and the difference of r b and r c : They are acting in the direction of e rad .

Spinning, Skidding and Gyroscopic Moments
To calculate the gyroscopic effects of the rolling elements, the balls' rotational velocity ω b is needed.The ball rotation axis ω b cannot be normal for both contact ellipses at the outer and inner ring contacts simultaneously because of the different angles α o and α i .This causes relative rotation velocity at the contact points, called here spinning motion and denoted by ω s,o and ω s,i at the raceways (see Figure 4).It generates the power loss P loss and wear on the surfaces due to the friction forces at the contact ellipse.The determination of the ball rotation ω b is not possible using only geometrical constraints.
A modern approach for this calculation of ω b is made by [3,4].It is applied to every ball in this dynamic bearing model.The methodology is to minimize the ball spinning power loss P loss to get the ball's rotation angle α b .Therefore, the ball spinning motion for each contact is calculated depending on the ball rotation vector: ω s = f ( ω b ).The power loss is then calculated using the spin vector ω s and the relative skidding velocity ∆ v to: where v(P) is the slip velocity and τ(P) = p(P) • µ is the friction shear stress in the contact ellipse at point P, the integration taking place over the entire ellipse (see Figure 5).The chosen dry coulomb friction coefficient for steel is µ = 0.1 and p(P) is the local hertz pressure.The local slip velocity v(P) at point P is: The skidding velocity ∆ v is calculated from the cage rotation ω c and the contact point velocities assuming that the skidding velocity at both contact points is equal The calculated power loss is then added for both contacts for each ball and is a function of the ball rotation vector P loss = f ( ω b ) and its angle α b and magnitude.A numerical minimization of the power loss leads to ω b and must be performed in each timestep.The reached power loss and the spinning motion are a part of the whole bearing losses and can be used for frictional, thermal or wear calculations.
The gyroscopic moments caused by the misalignment between the balls' rotation ω b and the bearing axis are shown in Figure 1.It is defined with the balls' moment of inertia J b as: We assume that these gyroscopic moments are in equilibrium with the moments created on the ball by the tangential forces F t (Figure 1), acting at the ball-race contact points.

Damping in Ball Race Contact
In our dry bearing case, there is only little damping compared to the lubricated one, where the elastohydrodynamic oil film generates the damping forces.However, the cyclic deformation of a linear-elastic material causes energy losses that correspond to the hysteresis on a load-displacement diagram.For a specific material, the loss factor Ψ d corresponds to the dissipated energy E D and the strain energy E S over a deformation cycle (see [5]): For this model, Ψ d = 1.5% is chosen for each bearing contact.To propose a viscous damping coefficient c, the energy dissipation is compared for a 1D oscillator system.The resulting coefficient c is a function of Ψ d , the Hertzian contact stiffness k h and the deformation frequency ω bp = || ω c ||/n k , the ball passing frequency (n k is the number of balls), to [5]:

Modeling of the Elastic Outer Ring
So far, the local penetration between ball and outer raceway is taken into account due to the Hertzian contact model (see Section 2.1).In order to cover the global deflection of the outer raceway in addition to the local one, an elastic model for the bearing outer ring is considered.The inner ring elasticity can also be modeled.As an example, only the elasticity of the outer ring is considered.In this work, it is represented by a reduced finite element model with linear elastic material behavior.
The following steps describe representatively the process of the outer ring modeling.First, the physical dimensions define the geometry of the elastic ring.Second, a discretization by the finite element method gives the governing equations for the full outer ring model.In the last step, an adequate reduction method decreases the size of the full model for a later efficient consideration in the dynamic simulation.Here, the Craig-Bampton approach is applied.

Geometry
The geometry of the bearing outer ring is defined by the physical dimensions of the rigid bearing.As an example, without loss of generality, the outer ring is approximated in the following by a cuboid of dimensions D o × D o × W o with the side length D o and the width W o of the outer ring.The cuboid has a hollow cylinder with diameter d o representing the outer raceway inside the bearing ring.Despite the simple FE geometry, the detailed geometric curvature of the outer raceway is still considered by an analytical superposition.Thus, it is possible to choose a very crude or highly reduced FE mesh and to represent a real bearing raceway geometry (represented by the curvature radii and its centers) for calculating the deflections and contact forces.

Finite Element Approach
In the next step, a finite element approach is used in order to discretize the structural equations of the bearing ring.This step is performed in a finite element software tool.The elastic volume is discretized by using three-dimensional incomplete quadrilateral elements (each with 20 nodes) with serendipity shape functions.A structured mesh is applied to the interface of the hollow cylinder by using 40 elements in circumferential and four elements in axial directions, respectively.Bi-quadratic finite elements (each with eight nodes) are used for the interface.The degrees of freedom of the nodes of the cuboid's bottom are fixed.Figure 6 shows the finite element model of the bearing outer ring structure.The homogeneous space-discretized equation of the full finite element model is stated as

Reduction Methodology
Due to the large number of elastic degrees of freedom, the full finite element model can hardly be handled in a dynamic multibody simulation.Therefore, in the following, the full model with its N degrees of freedom is reduced by the classical Craig-Bampton method.
For the Craig-Bampton reduction, the vector x of the full model is decomposed into a vector x b containing n b displacements of boundary nodes and a vector x i containing n i displacements of inner nodes.The degrees of freedom are thus partitioned as with the total number of nodes N = n b + n i .As the boundary nodes are later retained for the reduced model, they are chosen as the nodes that belong to the circle of the outer raceway, indicated in blue in Figure 6.The residual nodes correspond to the inner nodes of the model.Then, a reduction x = T q el by the Craig-Bampton approach is defined as: The first columns of the reduction basis T are n b interface constraint modes based on unit displacements of the boundary coordinates x i .The remaining columns contained in the matrix Φ v a set of n v n i fixed interface normal modes.They are determined by restraining all boundary nodes and solving the obtained eigenvalue problem for the first n v eigenvalues.Using the reduction of Equation ( 14), the following reduced mass and stiffness matrices are obtained: In the later example, 80 boundary nodes and hence n b = 3 × 80 constraint modes and n v = 10 vibration modes are used.
It is noteworthy that other approaches could also be used for the reduction of the full outer ring structure model.For instance, the already reduced model could be further reduced by a modal truncation approach, as followed by Novotny [17] for elastic journal bearings.In addition, a load dependent approach by using free interface normal modes and attachment modes could have been followed (see [18][19][20]).

Full Dynamic Bearing Model
The bearing assembly can be seen as two different parts: The bearing elements (the balls, cage and inner ring) and the elastic outer ring structure.The movements of the bearing elements lead to a multibody simulation, in which the elastic outer ring is a structural dynamics element.
The two parts are coupled using projections of the displacements and forces with interpolations to the constraint modes.This is discussed in detail next.

Multibody Simulation
Modeling the bearing components, the kinematics and kinetics lead to a dynamic equilibrium for each bearing component.For the inner ring, the dynamic equation writes: where h i is the vector of external forces on the inner ring.For the dynamic equation of motion of each ball k follows, using r k = (e ax e rad ) T : where the matrix A k transforms the ball forces into the local ball coordinate system.The centrifugal forces are summarized in the vector h k .The outer ring follows with r o = x, and the equation: where the Jacobi matrix J T k projects the k-th ball forces into the direction of the coordinates of the finite element model of the outer ring.
These three equations are summarized to a set of dynamic equations of motion of the whole bearing.A time integration scheme (MATLAB ode15s, R2016b, 7 September 2016, The MathWorks Inc., Natick, MA, USA) applied to the state space formulation of the system is used in order to get the component velocities and displacements over time.

Deflection Calculation
The contact deflection δ i and δ o of each ball on the inner and outer raceway, respectively, is calculated from the current state variables.The deflection δ i,k depends on the coordinate vector r i of the inner raceway and the coordinate vector r k of the k-th ball.The deflection δ o,k depends on the coordinate vector r k of the k-th ball and the current elastic deformation of the bearing housing at the ball angle Φ b,k (Figure 4), which is in the following denoted by d o,k .The latter deformation d o,k is interpolated from the deformation of the neighboring nodes by using the bi-quadratic shape functions used for the finite element discretization.

Force Projections
The ball forces F n,i , F t,i and F n,o , F t,o on the inner and outer raceway, respectively, are calculated on the basis of the local deflections (see Equation ( 3)).While they can be directly applied to the dynamic equations of the inner ring (see Equation ( 16)), they have to be transformed into the local coordinate system by a matrix A k , when applying them to the dynamic equation of the k-th ball (see Equation ( 18)).At the outer raceway, they have to be projected onto the elastic coordinates at the ball angle Φ b,k (Figure 4).Therefore, the Jacobi matrix J T k is used (see Equation ( 18)).It contains the bi-quadratic shape functions used for the finite element discretization and evaluated at the current position Φ b,k .

Bearing Model Implementation
The coupling of the two parts is done directly by projecting the elastic housing's deformation to the contact deflections of the outer raceway.The normal and tangential forces are then projected back from the raceway to the housing.In this way, it is possible to add the state space formulation of the housing directly to the bearing state vector.Thus, the overall dynamic equilibrium is solved for each timestep.

Outer Ring Structural Damping
To add damping to the elastic housing, two approaches are combined.The elastic housing modes are damped by 1% modal damping using the modal expansion to build the viscous damping matrix C el (see [21]): with µ s = x (s) M el x T (s) for the generalized mass of the state space mode x (s) .M el is the mass matrix of the elastic housing with stiffness matrix K el .
The static degrees of freedom of the housing (viscous damping matrix C static ) are damped using the Rayleigh approach: with α = β = 1%.The whole damping matrix C el is then assembled: To speed up the time integration scheme and to choose realistic initial conditions, the static equilibrium is solved beforehand.

Results of Representative Example
Showing the differences in the bearing behavior between a rigid and elastic outer ring formulation is the task of this chapter.
The examined ball bearing is a standard deep groove ball bearing of type 6404.Its outer diameter is 72 mm, the inner is 20 mm, and it includes six balls with a diameter of 15.1 mm.The material is steel and no lubrication is used.In the elastic housing case, the outer ring is replaced by the flexible structure.For a realistic rotor case simulation, some shaft mass (1 kg) and stiffness are added to the bearing's inner ring, which leads to neglecting the rotation of the inner ring system around e y and e z .
The results should not necessarily represent a realistic bearing.Instead, they are presented to show the differences of the model methodology and the advantage of detailed elastic formulations.
Figure 7 shows a snapshot of the dynamic simulation (time-dependent) of the bearing.It is loaded by a constant radial force of 1000 N in a positive y-direction (upwards in the figure).The shaft rotates at 1047 rpm.The green points represent the position of the elastic housing.It can be seen that the balls have a Hertzian deflection δ additionally to the housing's deformation.The deflections and deformations are amplified for better visualization.In contrast to a homogeneous black box formulation, the inhomogeneous load distribution can be seen clearly.

Bearing Displacements
The radial inner ring displacements (shaft movement) are shown in Figure 8.The bearing is loaded with radial force, in the y-direction, at several rotational speeds (3000, 6000, 12,000, 24,000 rpm) for a rigid and elastic housing formulation.
Within the rigid configurations, the speed dependency is only minor.The centrifugal forces are small with regard to the Hertzian stiffness.In the elastic housing case, a larger influence of the rotational speed can be seen.The higher the shaft speed, the higher the radial displacement.The absolute inner ring displacement of the elastic formulation is higher compared to the rigid one.This shows the weakening effect of a soft outer ring structure.Figure 9 shows the radial bearing stiffness calculated from the displacements.

Bearing Excitation Frequencies
The consideration of the bearing as an inhomogeneous, discrete parted element leads to a rotating angle-dependent behavior.Depicted in Figure 10 is the orbit plot of the inner ring motion of the constant radial loaded bearing (F rad =1000 N) at a constant rotational speed of 1000 rpm.It indicates a movement of the shaft in an elliptic orbit.The center location is not constant unlike the load, so the dynamic simulation is useful to predict the amplitudes.The frequency of this orbit is described mainly by the ball passing frequency (see Figure 11).It has to be considered to a rotor resonance at that frequency.In this case, the inner ring rotates at 1000 rpm, corresponding to 16.7 Hz.
The dominating excitation frequency is about 33.9 Hz.It excites higher frequencies of the elements and the housing.The ball passing frequency is determined by the cage speed ω c and the number of balls.An accurate calculation, presented in this paper, is useful for the machinery design process.

Conclusions
In this paper, a full dynamic multibody bearing simulation including an elastic outer ring is presented.The discrete elements of the bearing, the inner ring, the balls and the cage are described.To calculate the acting forces, the Hertzian theory and deflections between the bearing elements are used.To model a bearing integrated a lightweight structure, the bearing outer ring is modeled as an elastic structure using the FEM approach and reduction techniques.Projections of forces and displacements are used to couple the multibody simulation to the structural part.This formulation makes it possible to examine the ball-structure interaction and the whole dynamic bearing behavior.A time integration leads to the shaft and bearing elements dynamic motion.To compare the elastic and the rigid outer ring approach, some stiffness and displacement results are shown.The orbit movement can lead to a vibrational excitation of the rotor.The excitation frequency is calculated using a detailed physically motivated bearing cage approach.The model still contains simplifications like the assumption of constant radii for inner and outer raceway shapes, the neglect of cage inertia and clearance effects, the simplified FE geometry or the absence of lubrication.These simplifications give space for future investigations; for instance, the validation of the model by experimental data could provide deeper knowledge on the underlying effects.

Figure 2 .
Figure 2. Degrees of freedom and coordinate systems.

Figure 4 .
Figure 4. Velocities in the ball bearing.

Figure 5 .
Figure 5. Power loss, velocity and pressure in the contact ellipse.

Figure 6 .
Figure 6.Finite element model of an elastic bearing outer ring.

Figure 8 .
Figure 8. Inner ring displacement according to radial force at different rotational speeds.

Figure 9 .
Figure 9. Radial bearing stiffness according to radial force at different rotational speeds.

Figure 10 .
Figure 10.Inner ring orbit over time at constant radial load and rotational speed.

Figure 11 .
Figure 11.Inner ring displacement over time at constant radial load and rotational speed.