Multibody modeling of non-planar ball bearings

– This work presents the dynamic modeling of ball bearing which uses multibody dynamic formalism. Such formalism allows immediate integration of the model in dynamic simulations of helicopter main gear boxes. Ball bearing is considered non-lubricated in order to predict its behavior in case of lubrication system failure. Rolling contacts are treated with the method proposed by Kalker. This approach is based on polynomial approximation of relative displacement on the contact ellipse. For low computational cost and without any spatial discretization, it gives a good estimation of tangential traction and creep. Also, a regularization of the Kalker linear creep theory is proposed. It is used here to facilitate the global convergence of the Newton iterative process. It is well suited for multibody dynamic simulations which do not need a very ﬁne treatment of rolling contact. A numerical example of a ball bearing under thrust load is presented. cas de d´efaillance du syst`eme de lubriﬁcation. Le mod`ele mis en place pour les contacts roulants pseudo-ponctuels, issu des travaux de J.J. Kalker, se fonde sur une approximation po-lynomiale du d´eplacement relatif sur l’ellipse de contact. Ce dernier fournit, pour un temps de calcul r´eduit et sans discr´etisation spatiale, une bonne estimation des eﬀorts et micro-glissements au contact. Aussi, une r´egularisation de la th´eorie lin´eaire de Kalker est propos´ee. Elle est utilis´ee pour faciliter la convergence globale de l’algorithme de Newton. Elle est ´egalement bien adapt´ee aux simulations dynamiques multicorps qui ne n´ecessitent pas une mod´elisation tr`es ﬁne du contact roulant. Un roulement `a billes soumis `a un eﬀort axial est pr´esent´e comme exemple num´erique.


Introduction
In helicopters, the mechanical transmission of power is achieved by the Main Gear Boxes (M.G.B). Rolling bearings are the most critical component in such mechanical system since it permits the relative rotation motion between shaft and housing, and it carries high load and speed. The understanding of such components is essential. Thus the local behavior of rolling element bearings has been widely studied. Notable works include quasi-static analysis [1] and dynamic analysis with the well-known environment cannot be neglected. For instance, misalignment induced by housing deformation may cause damages to rolling bearings and a premature failure [8].
Dynamic modeling of multibody systems allows analysis, design, optimization and control of complex mechanical systems. Usually, joints are assumed to be ideal or simplified [9,10]. Studies dealing with the use of physical joint in multibody dynamic simulation are of growing interest. In multibody simulation, real cylindrical or spherical joint have been used [11,12]. A recent work proposes a model of planar deep groove ball bearing that can be used in planar multibody system [13].
The purpose of this study is to introduce a non-planar dynamic model of non-lubricated ball bearing with clearance which can be used in dynamic simulation of M.G.B. This model is directly integrable in multibody dynamic simulations. Its utilization is not limited to simplified shaft architecture as could be specialized codes for bearing simulation [14]. Also, specialized codes for gearboxes simulation are often limited to quasi-static analysis. Some of them, as RomaxDESIGNER c , propose a coupling with an external multibody dynamic software. In such simulations, rolling bearings are treated as nonlinear stiffness matrices. It allows efficient vibration analysis of gearboxes but it is unable to predict internal behaviour of the rolling bearing. The present model is dedicated to direct time integration of the equations of motions. Then, it furnishes the internal dynamic behaviour of the ball bearing.
Non lubricated model is useful to predict behavior of rolling bearings in case of lubrication system failure. In such model, the description of rolling contact is of prior importance. The number of simulated contacts becomes very important, which could lead to long CPU time. Kalker's approach [15] has been chosen since for low computational cost and without any spatial discretization, it gives an accurate estimation of tangential traction and creep. A simplified model, based on a regularization of Kalker linear theory has been developed in order to facilitate the convergence of the global nonlinear dynamic resolution. Rolling contact model deriving from Kalker linear theory has been studied [16].

Description of multibody formalism
Several formalisms are available in multibody dynamic. The present work uses the finite element approach [10]. A fundamental aspect in spatial multibody dynamic is the representation of large finite rotation. Most approaches use Euler angles or Euler parameters. The first method provides intuitive description of rotational motion but has the drawback to exhibit singularity. The second method avoids this singularity but needs one redundant variable. Their use requires imposing additional nonlinear constraint of normality. The present work employs Conformal Rotation Vector. This parameterization is free of any singularity over one full rotation and uses only three independent parameters [17]. A correction procedure of the Conformal Rotation Vector is used to allow rotation of any magnitude. Thus, each body needs the use of six parameters. All generalized coordinates are collected in a column vector q.
Bodies are usually interconnected with joints. These connections are expressed as algebraic equations. Holonomic constraints Φ are formulated as implicit functions of the generalized coordinates and time (1). Nonholonomic constraints Φ nh encountered in this work are always linear in velocities and can be expressed in the form (2) whereq notation denotes for time derivative of q.
The constrained dynamic problem is treated with the augmented Lagrangian method. Non-holonomic constraints are considered using a dissipation function D, where p is a penalty coefficient and k is a scaling factor. λ is a set of Lagrange multipliers. The superscript T is the transpose operator.
Dynamic equilibrium equations of constrained system are then obtained. M , B are respectively the mass matrix and the matrix of holonomic constraint gradient. The force vector g is the sum of external, internal and complementary inertia force. r is the dynamic residual vector.
This nonlinear second order differential-algebraic system (4) is solved using HHT scheme [18]. As Newmark method, it allows a direct resolution of second order differential system. It is unconditionally stable, second order accurate and induces less numerical dissipation than damped Newmark algorithm. This scheme needs the linearized form of the dynamic equilibrium (4).
Knowing an approximate solution (q ,q ,q , λ ,λ ) at time t and considering a correction of the solution (Δq, Δq, Δq, Δλ, Δλ), the differentiation of (4) provides: where K t and C t are the tangent stiffness and the tangent damping matrices. Finally, the nonlinear equation of motion is solved iteratively by Newton method.

Rolling contact model
The description of rolling contact is of prior importance because net traction force defines the time evolution of the ball bearing. Boussinesq has solved elastic contact between non-conforming surfaces. Relative elastic displacement u = (u, v, w) T , contact pressure p z and tangential traction f = (f x , f y ) T are linked through integral equations. The local Coulomb model is used, let μ be the friction coefficient. In this work, the coupling between normal and tangential contact problems is neglected. It is exact if both bodies have the same elastic properties. Thus, in case of contact between two non-conforming bodies, pressure distribution and contact area E are given by the Hertz theory of frictionless contact. The contact area is an ellipse, a and b are its semi-minor and semi-major axes respectively. The pressure distribution is of the form (7).
The normal load F n is linked with the contact deflection δ by (8). k h is the contact rigidity which depends on curvatures of contacting surfaces and material properties of contacting bodies. More information about Hertz contact theory can be found in [19]. Once the contact area and the contact pressure are obtained, the tangential problem of contact is solved.

Tangential contact problem
Kinematic description of contacting surfaces is needed to solve the tangential problem of contact. Let R be a fixed Cartesian coordinate system. Consider a particle M of body i which lies at time t in x in the undeformed state, and at x + u i in the deformed state. For each contacting body, the velocity of the particle is found by differentiating particle position with respect to time.
Assuming small strains and neglecting second order terms, the relative velocity of micro-slip between contacting points S(M, t) is given by: In case of stationary rolling contact and dividing by the rolling speed V r , (10) becomes: where, V r is the rolling speed. s = (s x , s y ) T , c = (c x , c y ) and ϕ z are respectively micro-slip, creep and spin velocities adimensionalized by the rolling speed. A classical assumption is to assume that the rolling direction follows the x axis of the contact ellipse which means that V rx = V r and V ry = 0. For ball bearings, such condition cannot be asserted.

Kalker nonlinear creep theory
Kalker's approach makes use of the fact that in case of material elastic symmetry, the load displacement equation is valid [15]. If we assume a traction of the form (12), then the relative displacements inside the elliptical area E are polynomials in x and y of the same degree. Polynomial coefficients (a mn , b mn ) and (d pq , e pq ) are linked linearly.
Classical methods [19] assume a priori the stick and slip zones and often lead to a closed form for the traction and creep. Nevertheless, it breaks down when the spin motion becomes large. The difficulty of such a problem lies in the different boundary conditions which have to be satisfied in slip and stick zone. Those areas are not known in advance. Kalker's approach avoids this difficulty by combining boundary conditions. Let T and S be defined by (13). In the stick area, S is null, whereas in the slip area tangential traction is equal to the maximal Coulomb value and T is null. Finally on the whole contact area E, the product (T S) must be null. Thus, the tangential contact problem is solved by minimizing the integral of (T S) over E. Since the value of the integrand of (T S) over E has no physical interest, the minimization process is made on an elliptical grid. Such an approach blurs the distinction between stick and slip zones, which are now identified a posteriori. A polynomial degree of M = 3 is satisfactory for most situations. This method gives net traction forces and power dissipation of the rolling contact.
Figure 1a compares net traction between rolling contact models and experimental data. Except the present model curve, other curves arise from [19]. The test case is a ball rolling on a plane with increasing longitudinal creep. Strip theory [19] gives also good results but it cannot be used when spin motion is present. Classical traction distribution induced by longitudinal creep is shown in Figure 1b, the stick zone is close to the leading edge (right) whereas slip zone is close to the trailing edge (left). Even if this rolling contact model gives good estimation of net traction for a slight computational cost, the underlying minimization process can break down for very unnatural situations. Unfortunately, such unnatural situations can arise from the Newton-Raphson process of dynamic resolution. Indeed, at each time step, first iterations of the Newton-Raphson process frequently lead to excessive unbalance, which produces convergence difficulties of nonlinear rolling contact models. To solve this problem, a simplified and more robust rolling contact model has been developed. This model is used until the dynamic residual has decreased sufficiently to avoid unnatural situations. These cases correspond to inconsistent values of (c x , c y , ϕ z ).

Kalker linear creep theory
Linear creep theory is limited to infinitesimal creep and also derived from load displacement equation [15]. It assumes complete stick of the contact area (11). The minimization of the surface traction on the leading edge of the contact area leads to a linear system. Finally, tangential forces are linear in creep and spin velocities. Creep coefficients C ij depend only on ellipse ratio a b and friction coefficient.
A simplified approach based on Kalker linear creep theory is compared to an exact approach by Chevalier [16]. This model needs very low computational cost. It is only valid for an infinitesimal creep since net traction can exceed the force of limiting friction. If contact curvatures are not changing, as in the case of contact between balls and cage pockets, creep coefficients C ij are calculated only one time. In ball bearings, the range of variation of creep and spin velocities is large. Such model is not applicable as it is, however it is a physical starting point of a regularized rolling contact model.

Regularized Kalker linear creep theory
It was shown in 3.1.1 that a simplified rolling contact model is needed to facilitate the global convergence of the Newton-Raphson procedure. This simplified model will be used until the absolute dynamic residual r is below a certain value ε switch . Then, Kalker nonlinear model is used in order to get accurate values of contact forces. To avoid a huge increase of the global dynamic residual after model switching, the simplified model has to be close enough to the nonlinear model. Also, the threshold has to be well chosen. A too high value of ε switch may lead to convergence difficulties of the contact model. It increases also the computational cost since more contacts will be treated with the nonlinear model. A too small value may cause a huge rise of the dynamic residual after model switching. Usually, the dynamic convergence threshold is 10 −6 . From numerical experiments, it appears that, a threshold ε switch of the order of magnitude of 10 −1 is sufficient to avoid convergence problems of the nonlinear contact model.
In the example presented in Section 5, the dynamic convergence threshold was 10 −6 . The threshold ε switch was fixed at 10 −3 , in order to avoid lack of convergence of the nonlinear contact model and save computational time. This choice is based on empirical considerations. Naturally, this value depends slightly on the ball bearing studied.
Contact models used in multibody dynamic are often a regularization of the Coulomb model. The regularization function (15) proposed by [20] is especially interesting since it guaranties function and derivative continuity.
A key point of such model is to choose the parameter . Very often, this choice is only based on numerical considerations because of its direct influence on the numerical behaviour of the model. Since it defines the tangent at zero creep, this regularization parameter has also a real physical meaning. The approach proposed here is to compute the regularization parameter which gives the slope of the Kalker linear theory for infinitesimal creep. Then, the model developed combines linear creep theory and regularization functions in order to get a fast and accurate rolling contact model. It is summed up in (16).
where (c z , c 0 , c ) are equivalent creep parameters (17), and φ ij are regularization functions. If the coupling between transverse creep and spin is neglected φ yz ≡ φ zy ≡ 0, the model never exceeds maximal friction force if φ x ≤ 1 and φ y ≤ 1. In our case, it is more interesting to take into account such coupling. Then, an artificial renormalization, with preservation of the force direction, is used when resulting force exceeds maximal friction force. The maximal spin moment is computed assuming complete slip of the contact area.
For φ x and φ y , it appears that regularization functions of the form (15) were sufficient. For both functions, regularization function parameters are chosen to equal Kalker linear model at zero creep. Other regularization functions are rational fractions. Some of them need extra parameters which have been identified by minimizing the leastsquares error between regularized model and Kalker nonlinear model. Each identification process is made with only one creep parameter not null. For instance, parameters identification of φ yz is done by minimizing transverse tangential force error for rolling contacts with spin but without transverse creep. Fifty uniformly spaced spin values were used for this identification. A similar procedure was followed for the identification of φ zy and φ z . It seems that these extra parameters are sensitive to friction coefficient but not to ellipse ratio. Then, for a ball bearing dynamic simulation, only one identification process of these extra parameters is needed. A detailed expression of regularization functions can be found in Appendix A. Figure 2 compares the nonlinear creep theory and regularized contact model for ball rolling on a plane with only longitudinal creep. A good agreement is found between models even if it appears that the regularized model slightly overestimates the friction force. Figure 3 shows the case of combined creep. The longitudinal creep remains constant while the transverse creep increases. A good agreement is also found. One can remark that Kalker linear model grossly overestimates net traction. Figure 4 shows an arbitrary test case which combines longitudinal creep, transverse creep and spin. One can note that both Kalker nonlinear model and regularized model reproduce well camber thrust phenomenon. This effect is not predicted by complete slip theory [19] since it is due to tangential elastic compliance of the surface. The increase of creep case complexity leads to a higher discrepancy between regularized model and Kalker nonlinear model.
Even if an estimation of the dissipated power can be defined as (18), this estimation is inconsistent since linear model assumes complete stick of the contact zone. Moreover, due to the coupling between transverse creep and spin, this dissipated power could become positive which enforces its non-physical meaning. Anyway, the main objective of such a simplified model is to facilitate global convergence of the Newton-Raphson procedure. However, this rolling contact model can be used in multibody dynamic simulations which do not need a very fine treatment of rolling contact.

Rolling contact in multibody dynamic context
The rolling contact model presented previously is replaced in the multibody dynamic formalism. Contact force is still separated in normal and tangential components. The normal force only depends on generalized coordinates whereas F T depends also on generalized coordinates derivatives. The normal load F n is expressed in the form (19). Again, δ is the contact deflection. The contact rigidity k h depends on curvatures of contacting surfaces and material properties. Finally, contact force F is expressed in the form (20) where n is the normal vector.
The resolution of equilibrium equations (4) needs the computation of tangent operators. If the computation of such operators is exact, the convergence of the Newton-Raphson procedure is quadratic. However, such exact linearized form is often too costly to obtain. An approximate linearized form is frequently more efficient in terms of computational cost. The contribution of the normal load in the tangent stiffness operator is computed considering a small variation of the normal load (21). In this expression, the influence of contact rigidity variation is neglected. Since geometrical interactions are computed analytically, explicit expressions of δ and n are available.
Thus, an approximate explicit form of this contribution is obtained.
ΔF n (q) = ΔF n (q) n(q) + F n (q) Δn(q) The contribution of the tangential load in tangent operators is more cumbersome to derive. The tangential force can be expressed in the form (22). The variation of tangential load is then given by (23). The first term in (23) is already computed from (21). Because of the underlying minimization process in the nonlinear rolling contact model, an explicit expression of the contribution of ΔF T is really hard to achieve. This contribution can be approximated by finite differences. From numerical experiments, it appears that a numerical evaluation of this contribution is costly and does not improve so much the convergence of the Newton-Raphson process. Thus, this contribution is neglected.

Dynamic ball bearing model description
Ball bearing consists of three parts: a number of rolling elements, the outer and the inner races and the cage. The interaction between the bearing elements constitutes the basic formulation of the dynamic model. In case of local computation, specific boundary conditions are applied on races. For global simulation, the races are linked with the shaft and the housing. In case of incompatible model e.g. rigid races and finite element shaft and housing, rigid link can be used.

Boundary conditions
For a local study, the outer race is linked to a reference frame using hinge joint [10]. The reference frame could be an absolute frame or a fixed frame which consider housing misalignment due to housing flexibility. Such reference frame is derived from static FEM analysis of M.G.B. The inner race is submitted to external load and is free to move to obtain dynamic equilibrium. Any force or moment can be applied to the inner race. Both races have prescribed rotation velocities around their own rotation axis via non holonomic constraint of type: where Ω, z and ω d (t) are respectively the vector of material angular velocity, the axis of rotation of the race and the prescribed angular velocity. According to rotation parametrization, a matrix T links angular velocity  of the body with rotation parameters and its derivatives. Identification with (2) gives:

Internal interactions
Internal interactions are classified into two main parts, interaction between rolling elements and races or cage, and interaction between cage and race. These interactions are treated with compliant contact models. With such models, contact forces are expressed as continuous functions of contact parameters and contact geometry. Knowing positions and velocities of contacting bodies, the calculation of geometric interactions gives these contact parameters as for example, indentation, slip and rolling velocities. The geometric model of the ball bearing is shown in Figure 6. All geometric calculations are done analytically. Then interaction model provides contact forces and moments.

Rolling elements and races or cage interactions
Normal load F n of point contacts between balls and races or cage is obtained with the classical Hertz theory of elastic contact (19). The computation of tangential forces F T and moments is done using rolling contact models based on the work of Kalker. The approach followed is detailed in Section 3.

Cage and race interactions
Unlike punctual contact, Hertz theory does not provide a uniquely defined load deflection relation for cage race interaction. The treatment of the interaction between cage and race follows the one proposed in [2]. It employs Lundberg empirical load deflection relation to determine normal load. Tangential forces are obtained assuming complete slip of the contact area.

Results and discussion
This section presents a local study of a ball bearing which dimensions are summarized in Figure 6. Outer race is fixed rigidly with the housing. A thrust load of 1500 N is applied on the inner race. The rotation velocity of the inner race is 6000 rev.min −1 . Estimated friction coefficient between balls and races is 0.08. Quasi-static simulation provides initial values of generalized coordinates and initial velocities are obtained assuming pure rolling at contact points. Because of approximate initial values, the ball bearing behavior exhibits a stabilization period (see Fig. 7). In order to facilitate convergence, a simplified contact model is used during this period. This simplified model is the regularized rolling contact model presented in 3.1.3 in which coupling between transverse creep and spin was neglected. The HHT method is used to integrate the equation of motion. The damping parameter α f , which has to be in (0, 1 3 ] is chosen equal to 0.05. The time step is selected automatically. The automatic step size control is based on the computation of an estimate of the local truncation error (see Chap. 11 in [10]). For this simulation the time step was initialized at 10 −6 s , it lied between 10 −8 s and 10 −5 s. Since the present study focused on contact power dissipation, only absolute dynamic residual was adapted. Especially here, power dissipation of rolling contacts was small quantities since contact conditions were close to pure rolling. The absolute unbalance allowed is 10 −6 N. For constraints, the convergence threshold is 10 −8 .
Only results of the first ball are shown since each ball gets the same behavior. Figure 7a shows time evolution of adimensionalized creep velocities and net traction. Switching between simplified model and Kalker model produces, essentially, transverse creep. It was expected since in this model, spin motion and transverse tangential force are not uncoupled anymore. Even when creep velocity tends to zero, net traction is not null because spin motion produces transverse tangential force. This is coherent with experimental observations [19]. As can be seen in Figure 7b, the simplified model overestimates spin moment. The system exhibits oscillating behavior which amplitude slightly decreases over time because of numerical damping and friction. The oscillating behavior is mainly due to Hertz contact theory. This contact force model is equivalent to a nonlinear spring which rigidity depends on geometry and material properties. Balls are compressed between two nonlinear springs. Therefore, the system is prone to oscillating behavior. In order to limit this phenomenon, a possibility is to use dissipative contact force models. A study of some of relevant compliant contact force models for multibody systems dynamics can be found in [21].
Mean power dissipation of the rolling contact between ball and inner ring is 74 mW; outer contact 3 mW. Greater dissipation of inner contact was expected since inertial force tends to unload inner contact. Thus it is more disposed to get relative motion.

Conclusions
This work presented the dynamic modeling of nonplanar ball bearing. It combined multibody dynamic formalism and realistic rolling contact model. Ball bearing was considered non-lubricated in order to predict its behavior in case of lubrication system failure. Normal contact forces on rolling elements were computed using Hertz contact theory. Kalker nonlinear creep theory was used to solve tangential rolling contact problems between balls and races. Without spatial discretization, it gave accurate estimation of net traction, moment, creep and power dissipation of rolling contacts. In order to facilitate global convergence of the Newton-Raphson algorithm used to solve dynamic equilibrium, a new rolling contact model has been developed. This model combines Kalker linear creep theory and regularization function. Such model is used until the dynamic residual is below a certain value. It enforces robustness of the simulation and saves CPU time. Ball bearing dynamic behavior was integrated using HHT scheme and boundary conditions were imposed with augmented Lagrangian method. A numerical example of a ball bearing under thrust load was used to show its dynamic behavior and rolling contact characteristics.
The use of realistic joint in multibody dynamic simulation becomes essential especially for flexible mechanisms. The model presented is immediately usable in multibody dynamic simulations. It is the starting point for integration of dynamic behavior of helicopters main gear boxes.