Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings David Noel, Mathieu Ritou, Benoît Furet, Sébastien Le Loch


Introduction
The dynamic behavior of a high speed rotor needs to be studied during its design.The dynamic stability is one of the most important criteria.Experimental tests can determine the instability, however, their major drawbacks are their cost, thus limiting the possible number of tests.With a numerical modeling approach, the stability can be analyzed in detail.As an example, numerical models contribute to choosing the optimum set of operating conditions in high speed machining (HSM) in order to maximize the productivity while ensuring workpiece quality and spindle health [1].Experimental steps are still necessary to update and validate the model.Updating would increase the accuracy, thus providing the values of parameters such as damping [2].
One main part of the model relies on the guidance of the rotor.The paper deals with the most widespread type: the angular contact ball bearing guided shaft.Its model, associated with a finite element model, enables the building of a global model of the rotor [3,4].In this process, the accuracy of the bearing stiffness values is crucial.Indeed, the frequency response function and, in particular, the critical eigenfrequencies are directly linked to the bearings stiffness [5].
Bearings manufacturers have developed hybrid bearings for high speed rotors.The balls are made of ceramic instead of steel, which decreases the dynamic effects on the balls and the subsequent load on the outer race at high speed due to centrifugal forces.Thus, these high precision components can reach more than 2:510 6 d m N (ball orbital diameter in millimeters times the shaft speed in rpm), which is characteristic of very critical applications [1].
The mechanical model of the bearing aims at obtaining not only the relation between the global loads f and global displacement d, but also the stiffness values.The models are characterized by the number of degrees of freedom (DOF).Since no explicit constitutive relation can be obtained for multi-DOF models, local variables (normal loads ðQ i ; Q o Þ and normal contact deformations ðd i ; d o Þ) need to be considered for each ball.Historically, two different approaches are possible.
The analytical approach is based on a hypothesis on the local load distribution.Sjova ¨ll [6] first proposed the use of the load distribution factor "e" and the Sjova ¨ll integrals.Palmgren [7] popularized this method with a 2-DOF model.Houpert [8] generalized the method with a 5-DOF model.Hernot et al. [9] computed the stiffness matrix.Unfortunately, the dynamics effects on the balls cannot be taken into account with this approach.
The numerical approach is based on a rigid body displacement hypothesis of the inner ring (see Fig. 1).The approach has been introduced by Jones [10] with a 5-DOF model.The main steps of the model's elaboration are exposed in Refs.[4,10,11] and the stiffness matrix computation is detailed in Refs.[4,10,12].With the numerical approach, the dynamic effects on the balls are taken into account during local and global equilibriums.Two types of dynamic effects are considered (see Fig. 2).The centrifugal force F c is induced by the ball orbiting at the diameter d m .A gyroscopic moment M g results from the shifting direction of the ball's spin axis.
Previous works proved that dynamic effects play a significant role in spindle behavior [2,4]: the dynamic effects on the shaft and on the balls of the bearings.Indeed the stiffness varies according to the shaft speed.Therefore, the numerical approach has been chosen in the paper to express the relation between global loads and displacements.
After presenting a 5-DOF bearing model and the different kinematic assumptions, the numerical computation of the bearing model is explained.A new complete analytical expression of the stiffness matrix, taking into account the dynamic effects on balls, is detailed.In a second section, the improvement of the accuracy and the Fig. 1 Principle of the bearing model related computational cost is then studied.The impact of the kinematic hypotheses is discussed.Finally, the relevance of the method is shown by studying the first modes of a high speed spindle.
2 Dynamic Model of Bearing 2.1 Elaboration of the Model.This subsection explains the main steps of the model elaboration schematized in Fig. 1.

Rigid Body Displacement and Geometric
Relations.The geometry at the contact is shown in Fig. 3, corresponding to an azimuth location w of a ball defined in Fig. 4. The relative displacement at the ball's location is established from a rigid body displacement hypothesis of the inner ring According to Fig. 3, the Pythagorean theorem gives two geometric equations In this formulation, the global deformations of the rings are not taken into account.Raceways are assumed circular and without radial expansion.Ring deformations can be computed with a finite element model or with a continuum mechanics model.The present method can also be refined with radial ring expansions [13].

Local Equilibrium.
For each ball, Newton's second law of motion leads to (see Fig. 4) Equation ( 3) and Fig. 4 suppose that contact surfaces can provide sufficient reaction forces to the gyroscopic moment M g , i.e., k o M g =D lQ o and k i M g =D lQ i where l is the friction coefficient at the contacts.This assumption is reasonable for spindles since a high preload is applied on the bearings.The coefficients k i and k o express the gyroscopic moment distribution (see Sec. 2.2.3).The local variables of Eqs.(2) and (3) are x ¼ ðX 1 ; X 2 ; d o ; d i Þ.The contact angles ða o ; a i Þ and local loads ðQ o ; Q i Þ are expressed as a function of x, respectively, with the geometric relations of Fig. 3 and with the Hertz constitutive relation [14]: d ¼ KQ 2=3 .

Global Equilibrium.
Once the local equilibrium has been solved, all of the local loads and contact angles are known.The global loads f are obtained by summing each ball's quantities, according to Newton's second law applied on the inner ring  2.2 Details on Bearing Dynamics.The ball motion is driven by all forces applied on the ball, notably the distribution of friction forces on both raceways.It results from complex physical phenomena.Moreover, these phenomena rely on ball equilibrium, i.e., local loads ðQ i ; Q o Þ and contact angles ða i ; a o Þ.The variables are still to be determined at this stage.To lower the complexity of the bearing model, several assumptions are generally made [4,10,11].A constant shaft speed and no gross slippage are assumed.The ball rotational speed about the tangent of the pitch circle is neglected.More importantly, a hypothesis is made concerning the pitch angle b, as detailed in Sec.2.2.2.[11] determined the centrifugal force and gyroscopic moment expressions as

Expressions of Dynamic Effects. Harris and Kotzalas
For a rotating inner ring, Harris and Katzalas demonstrated the following speed ratios [11] x The pitch angle b corresponds to the angle between the bearing axis and the rolling axis of the ball; see Fig. 2. Its expression is missing for a complete determination of the kinematic model.A hypothesis must be made in order to be able to express this angle.The following paragraph summarizes the different hypotheses found in the literature and another one introduced from a geometrical relation.

Kinematic Hypotheses and Pitch Angle b
Control theory.The Jones's theory, called the "control theory," supposes that there is pure rolling on one raceway and rolling and spinning on the other one [10].On the controlling raceway, the spinning motion x s is, therefore, assumed nonexistent (see Fig. 5).For instance, concerning the outer-race control, the spinning motion x so is null.The pitch angle b is expressed according to each assumption (see Table 1).An additional criterion validates if the chosen control race is adapted to a given configuration (see Table 1).The criterion is based on a comparison of the frictional moments required to spin the ball about the normal directions to the contact surfaces.The controlling raceway is the raceway for which the moment is higher.
Figure 6 illustrates the control regions in relation to the axial load F x and shaft speed x.The bearing studied in Sec. 3 has been adopted for this uniaxial example.The small region between the two boundaries B irc and B orc corresponds to a region where both inner and outer-race control criteria are met.A drawback of this theory is that the chosen control theory must be verified afterwards, according to the criterion.A new computation eventually needs to be done if the initial choice was not the proper one.For the uniaxial case, the boundaries have been fitted with quadratic equations.Once the analytic expressions of the limits have been established for the considered bearing, they can be used to first choose the adequate theory for a uniaxial loading case.
The control theory has been mostly used for the past few decades, even if it is not always experimentally verified.Nevertheless, the outer-race control is adapted to an oil-air lubricated bearing under light loading and high speed, as discussed in Ref. [11].
Hypothesis based on d'Alembert's principle.Changan et al. [15] established a new formula to determine the pitch angle b by applying d'Alembert's principle on balls.The ratio C of the frictional spinning moment is considered with Contrary to the control theory, a unique expression continuously varying is obtained, regardless of the values of the speed and load.The results seem to be more accurate with this pitch angle calculation [15,16].In the following, this theory will be referred to as the "hybrid" theory.
This hypothesis corresponds to the assumption that the spin-toroll ratios at each contact area are equal The hypothesis leads to simpler equations for x m and x R .Even if the hypothesis seems, at first, limited compared to the other ones, it is retained in the paper for further investigation.In the following, this theory will be referred to as the "geometric" theory.

Contact Hypothesis: Distribution of the Gyroscopic
Moment.The contact hypothesis concerns the gyroscopic moment distribution between the inner and outer raceways.The distribution is taken into account through the coefficients k i and k o (see Fig. 4).
In previous works, tangential loads on the ball, corresponding to the reaction to the gyroscopic moment, are assumed equal on both raceways [4,10,17].Others associate a distribution according to the kinematic hypothesis [11,15,16,18].The associated distribution parameters k i and k o are given in Table 2.
In this work, tangential loads on balls are considered equal on both raceways k i ¼ k o ¼ 1 in order to enable comparisons of the stiffness computations with Ref. [4].

Model
Solving Methods.The model has been established in the previous subsection, depending on the application case.Either the load f is known and the displacement d is required and vice versa.Thus, two types of global solving methods are needed and exposed in this subsection.For both, local solutions, i.e., for each ball, are obtained first.It is numerically solved.

Local Solution.
In this first step, the global displacement d is known.The distance between the groove curvature centers is computed with Eq. (1).Each ball equilibrium needs to be solved.The four equations (see Eqs. ( 2) and ( 3)) expressed with the variables x are simultaneously solved.Since there is no explicit solution for the multi-DOF model, a numerical algorithm is required.The Newton-Raphson method is commonly chosen due to its relative simplicity.Here, e n is the residue vector of the four equations at iteration number n The matrix A is built from the partial derivatives of the four equations with respect to the local variables x.The matrix A is discussed in detail in Refs.[4,10].
It should be noted that the convergence of the algorithm depends on the starting point.In the present case, all functions are not defined on the real field.To avoid any problem, the solution of the static uniaxial case analytically obtained should be chosen as a starting point.

Global Load Determination.
The load determination consists of computing the load f that needs to be applied on the inner ring to obtain a given displacement d.Once each ball equilibrium is solved (see Sec. 2.3.1), the local loads ðQ i ; Q o Þ are summed according to Eq. ( 4) in order to obtain f.

Global Displacement Determination.
The displacement determination is the inverse of the load determination.It consists of computing the displacement d resulting from the application of a given load f.It is required to study the behavior of the spindle bearings.A new Newton-Raphson algorithm is proposed.It is based on the load determination method (see Fig. 7) where K n is the stiffness matrix and f n the result of the load determination at iteration number n.The implementation of the algorithm is simple once the stiffness matrix is made available.Its computational cost is discussed in Sec.3.1.

Computation of the Stiffness Matrix.
The stiffness matrix K is needed for the method of global displacement determination.Moreover, it is required in order to investigate the dynamic behavior of rotors, as discussed in Sec.3.4.The K matrix is a linearized relation between the loads Df and displacements Dd, for a given loading state.This is a 5 Â 5 Jacobian matrix built from the partial derivatives of loads with respect to the displacements: K ¼ ½@f=@d.

Numerical
Computation by Finite Differences.The computation by finite differences is the easiest method to implement because no analytic expression is needed.The stiffness matrix is calculated from a numerical gradient obtained after five additional load determinations (for the five DOF) close to the imposed displacement configuration.The method does not work for a displacement determination and is time-consuming (see Sec. 3.1).However, the numerical computation has the benefit of providing reference values for the validation of analytical computations of the stiffness matrix.

Analytical
Computation.This paragraph presents a new method for the analytical computation of the stiffness matrix.The method is based on Refs.[4,10] and is extended to accurately take into account the dynamic effects.In this way, the complete analytical expression of the stiffness matrix is achieved.This is particularly suitable for the computation of critical bearing applications (with high d m N), such as high speed spindles.
The stiffness matrix is K ¼ ½@f=@d.Since the load f results from the sum of each ball load (see Eq. ( 4)), the bearing stiffness results from the stiffness contribution of each ball.Thus, the matrix is built from a sum of terms corresponding to each ball.Therefore, each entry of K is expressed by a linear combination of elementary partial derivatives of the dynamic variables Here, i, j, k, and z are, respectively, the indices of the global loads, global displacements, dynamic local variables, and balls.The @f i =@x k d , @f i =@A 1 , and @f i =@A 2 terms are analytically obtained by differentiating the equations of the global equilibrium (see Eq. ( 4)).Trigonometric manipulations related to Fig. 3 are necessary.The @A 1 =@d j and @A 2 =@d j terms are easily obtained from Eq. ( 1).
The difficulty consists of obtaining the remaining terms @x k d =@d j .For this purpose, the six other Eqs.( 2), (3), and ( 5) are differentiated with respect to the global displacements d j .For instance the equation of the centrifugal force F c gives The six differentiated equations lead to a linear combination of @x d =@d j .This can be formulated as a matrix product B:@x d =@d j ¼ s j 0 0 0 0 with The upper left terms B ij with (i, j) from 1 to 4 correspond to the matrix A evoked in Eq. (10) and detailed in Refs.[4,10].The vector s j contains known terms concerning A 1 and A 2 .For instance Here, B is a 6 Â 6 matrix that is, in fact, independent of the DOF j.The elementary partial derivatives @x d =@d j are simultaneously computed by inverting the matrix (see Eq. ( 14)).The numerical computation is carried out for each ball and for each DOF j.
For each kinematic hypothesis seen in Sec.2.2.2, the analytical expressions corresponding to the dynamic effects are different (the two last rows of B and s j ).A symbolic software is useful in order to establish these expressions.For example, B 64 for the outer-race control is given in the Appendix.
This method differs from the literature.The most accomplished method [4] is equivalent to neglecting the terms B 35 , B 36 , and B 46 .In that work, @x k d =@d j for k 2 f1; 2; 3; 4g are computed apart and @x k d =@d j for k 2 f5; 6g are computed afterwards.This simplification decreases the complexity of the implementation because F c and M g do not need to be analytically expressed in the function of The relevance of the new method, based on the complete analytical expression of the stiffness matrix, is investigated in the following section.

Results and Comparison
The five computational methods studied in this section are presented in Table 3.Each method will be referred by its number.All configurations are covered.
The following high precision bearing has been retained for the study: SNFA VEX70/NS9CE3.It corresponds to a hybrid angular contact ball bearing that can be found in high speed and high power spindles.Table 4 presents its characteristics.
3.1 Computational Cost.In this subsection, computation costs are compared for the methods presented in Table 3.The computational costs for each method are given in Table 5 in percentages, relative to the reference case number 1.
The displacement determination is more costly because of the additional loop of the algorithm (see Fig. 7).In the example, the algorithm expressed by Eq. ( 11) converges in six steps.The differences between the analytical methods are minor for both the load and displacement determination methods (respectively, methods (1) to ( 3) and ( 2) to ( 4)).In fact, only some additional elementary operations are added.The new simultaneous computation of the partial derivatives of the dynamic effects, presented in Sec.2.4.2, does not significantly increase the computational cost.
The finite difference method 5 is more costly than the analytical computations 1 and 3. Indeed, six load determinations have to be carried out to compute the K matrix.The implementation of the analytical stiffness matrix is, therefore, beneficial because the computational cost is almost divided by 6 (the difference between methods 5 to 1) and also because it enables us to use the displacement determination computation.

Stiffness Matrix
Accuracy.The accuracy of the stiffness matrix is studied through the axial stiffness K 11 evaluated in the uniaxial case and considering the outer-race control hypothesis.The results from the load and displacement determination are identical because the stop criterion "tol" of Fig. 7 is small enough.  Figure 8 presents the relative error on the axial stiffness of the previous analytical method 1, compared with the new one 3 The relative error increases with the shaft speed x and decreases with the axial load F x .This maximum stiffness error corresponds to the configuration for which dynamics effects are significant compared to local loads.The error reaches 16% in the example.The stiffness error is significantly higher for steel ball bearings due to the ball density (the magnitude of dynamic effects is 2.5 higher).
The simplification adopted in previous work [4] was correct because it dealt with moderate d m N applications.Nevertheless, here, the stiffness error is too important for accurately modeling the dynamic behavior of bearings in the case of critical applications (high d m N).Indeed, Fig. 8 attests that the new stiffness computation is particularly interesting beyond 15,000 rpm.
In the example, the maximum error between methods 3 and 5 is less than 0.01%, regardless of the loads and speeds.Besides, the finite difference method is known to provide the reference results after a proper convergence study.Thus, we can conclude that the new analytical expression is complete because it leads to the exact results, contrarily to the previous ones.

Kinematic Hypotheses.
The kinematic hypothesis concerns the computation of the pitch angle b (see Sec. 2.2.2).Its impact on the stiffness is presented here.The axial stiffness K 11 is evaluated with the new analytical method for a medium preload of 1220 N on the bearing (see Fig. 9).
The loss of stiffness at high speed is due to the dynamic effects on the balls.As seen in Fig. 9, the stiffness drop is significant, which is characteristic of elastically preloaded bearings.Note that the kinematic hypothesis has a notable impact on the estimated stiffness at high speed.The estimated stiffness with the hybrid hypothesis is located between the values obtained with the inner and outer-race control, which is logical because the control theory supposes extreme situations concerning ball spinning motions.The hybrid theory is closer to the outer-race control at high speed, which is consistent with the boundaries presented in Fig. 5.The stiffness obtained with the geometric theory is almost in the middle of the values obtained with the inner and outer-race control because of the identical spin-to-roll ratios on both raceways.
In the example, the drop in stiffness in relation to the shaft speed is 75.7% and 78.3%, respectively, for the hybrid theory and outer-race control.The drops would be smaller for a more loaded bearing.The relative error for the axial stiffness K 11 between the outer-race control and hybrid theory is 10.4% at 30,000 rpm.
The main difference in the model comes from the estimation of the pitch angle b in Fig. 10.This angle barely changes the values of the speeds x R and x m or centrifugal forces F c .However, the impact on the gyroscopic moment M g is important (see Fig. 11).At high speed, the difference between the hybrid and control theory reaches 30% for M g .Hence, with a neglected gyroscopic moment, the impact of the kinematic hypothesis is quasi null.
Let us compare kinematic theories at high speed.Harris [11] determined that the outer-race control theory is valid at high speed, but only with light loads on the bearing.In that case, stiffness values obtained from the hybrid theory tend towards the outer-race theory values (even more than in Fig. 9).Indeed, x so then tends towards zero; therefore, the hybrid theory is suitable.In the same way, results from the hybrid theory tend towards the geometric theory ones with large loads, due to an increase of x so , which is also observed in Harris's book.Thus, the hybrid theory is also adapted with large loads.
In conclusion, the hybrid theory seems to be suitable for modeling the bearing kinematic because the behavior varies in relation to the bearing load.However, for a simplified implementation of the analytical stiffness matrix, the geometric hypothesis could provide a good estimation of the stiffness values for a significantly loaded bearing.
3.4 Impact on Spindle Dynamic Behavior.In Sec.3.2, a study of the stiffness accuracy is presented.It highlights the major error made at high speeds if the analytical expressions of the stiffness are not properly derived.The present subsection establishes the impact of this error on the computation of the eigenfrequencies of the spindle rotor.The impact of the kinematic hypothesis choice is also investigated.
The rotor volumetric model is presented in Fig. 12.The bearings are taken into account as linearized springs applied to the center of the inner ring (corresponding to each bearing DOF).The spring stiffness values are the diagonal terms of the stiffness matrix K.They are computed with the new analytical method proposed in the paper.The off diagonal terms of K have been neglected for the purposes of this paper.The effect of the off diagonal terms is discussed in Ref. [19].
The four bearings are identical (their characteristics are shown in Table 4).They are in a back-to-back arrangement of tandems and elastically preloaded by springs.Thus, the axial load F x of each bearing corresponds to half of the preload and is considered constant.For the example, the bearing preload is 1220 N, corresponding to the stiff preload recommended by the bearing manufacturer.The bearing inner rings, spacers, and motor parts are modeled as cylinder elements with mass but without stiffness.For the ease of comparison, the dynamic effects on the shaft are not considered in this simplified model of the rotor.
Catia V R finite element simulations have been carried out for the kinematic hypothesis of the outer-race control and hybrid theory.The simulations have also been made with results from methods 1 and 3 in order to observe the impact of the new computation method.The eigenfrequencies of the first two bending modes are presented in Table 6.
As for a 1-DOF mass spring-system, a drop of stiffness means lower eigenfrequencies, which explains the significant frequency evolution between 0 and 30,000 rpm.
The uncertainty on the eigenfrequencies resulting from the kinematic hypotheses (outer-race control or hybrid theory) is about 7% in this example.Additional work focusing on refining the bearing kinematics would probably enable us to reduce this uncertainty.
The new stiffness computation of method 3 avoids errors of, respectively, 6% and 4% for the outer-race control and the hybrid theory compared with the previous analytical method.The bending mode frequency errors vary from 16 Hz to 27 Hz between the different methods.This is crucial in order to predict stable conditions for high speed milling.The new complete analytical method is, therefore, of great interest in the elaboration of high speed spindle models.

Conclusion
After presenting a 5-DOF model for angular contact ball bearings, a new method was proposed to compute a stiffness matrix that properly takes into account the dynamic effects on the balls.Complete analytical expressions of stiffness matrix have been obtained.The numerical results have revealed an enhanced accuracy for the stiffness values.The new method is particularly relevant for lightly loaded bearings for critical applications (both the high values of shaft speed and ball orbital diameter).The computational cost has been analyzed for several solving methods and has not revealed any significant cost increase for the new stiffness computation method.Kinematic hypotheses concerning the ball spinning motion have been addressed.The hybrid theory seems more suitable because it takes into account loads using a friction ratio.The impact of the new method has been illustrated through a high speed spindle model.The eigenfrequencies have confirmed that the new method is crucial in order to accurately predict the dynamic behavior of the spindle.

Nomenclature Capital Letters
A 1 , A 2 ¼ distances between groove curvature centers (loaded) A ¼ Jacobian matrix for the computation of the local variables x B ¼ f i þ f o À1 gap between groove curvature centers (unloaded)

Appendix: Example of Analytical Expression in B Matrix
Considering the outer-race control hypothesis

Fig. 2 Fig. 3
Fig. 2 Dynamic effects on the balls

Fig. 5
Fig. 5 Spinning and rolling motion of the ball

Fig. 6
Fig.6Example of the inner and outer-race control regions for the uni-axial case, according to Jones's theory

Fig. 8 Fig. 9
Fig. 8 error on the axial stiffness DK 11 considering the uni-axial case under different loads F x and speed x

Fig. 12 CAD
Fig. 12 CAD model of the spindle rotor

Table 1
Criteria for the control theory and b expressions

Table 2
Gyroscopic moment distributions

Table 3
Solving methods

Table 5
Computational costs depending on the solving method

Table 6
Influence of the stiffness error on the eigenfrequencies values Jacobian matrix for the computation of the stiffness matrix D ¼ ball diameter E ¼ modulus of elasticity F c ¼ centrifugal force on ball J ¼ mass moment of inertia of ball K ¼ contact constant according to Hertz's theory K ¼ bearing stiffness matrix L ¼ complete elliptic integral of the second kind calculated in accordance to the Hertz theory M g ¼ ball gyroscopic moment Q ¼ ball-raceway normal load < ¼ radius to locus of raceway groove curvature center< i ¼ 0:5d m þ ðf i À 0:5ÞD cos a 0 X 1 , X 2 ¼ distancesbetween the inner groove curvature center and the ball center (loaded) Lowercase Letters d ¼ global displacement of the inner ring; d j corresponds to d z for j ¼ 3 d ¼ ðd x ; d y ; d z ; h y ; h z Þ d m ¼ ball orbital diameter f ¼ r/D f ¼ global load on inner ring (f and d are expressed at O h , the center of the outer ring; see Fig.4).f ¼ ðF x ; F y ; F z ; M y ; M z Þ m ¼ ball mass r ¼ raceway groove curvature radius x ¼ local variables x ¼ ðX 1 ; X 2 ; d o ; d i Þ x d ¼ local dynamic variables x d ¼ ðX 1 ; X 2 ; d o ; d i ; F c ; M g Þ z ¼ number of ballsGreek Symbolsa ¼ contact angles b ¼ ball pitch angle c ¼ D=d m d ¼ ball-raceway normal displacement k ¼ distribution parameter for the gyroscopic moment l ¼ friction coefficient at the contact ¼ Poisson's ratio q ¼ mass density w ¼ ball angular position on the pitch circle, azimuth position x ¼ shaft speed x m ¼ orbital speed of ball x R ¼ speed of ball about its own axis x roll ¼ rolling motion x s ¼ spinning motion