HSM Spindle Model Updating With Physical Phenomena Refinements

The modeling of High Speed Machining (HSM) spindles is a complex task due to the numerous physical phenomena involved in the dynamic behavior. Modeling is still rarely used in the industry, although sophisticated research work has been achieved. The boundary conditions of rotor models, which correspond to the ball bearings, are crucial and difficult to define. Indeed, they affect the dynamic behavior of the rotor in a non-linear and sometimes in an unpredictable way. The aim of the paper is to determine a relevant spindle model, i.e. the adequate level of complexity. To do so, a dynamic bearing model is introduced and the axial model of a spindle is established in relation to the preloaded bearing arrangement. Then, the operating stiffness of the spindle has been obtained experimentally with a new specific device that applies axial load and measures the resulting displacement, whatever the spindle speed. The model updating with the experimental data combined to sensibility analysis have led to the model refinement with additional physical phenomena, in order to account for non-linearities observed experimentally. The parameters of the model are also identified experimentally. As a result, a relevant spindle model is obtained and validated by the good agreement between simulations and experiments.Copyright © 2013 by ASME


Introduction
The numerical modeling of HSM spindle is a challenge for many scientists and many industrial user. The use of numerical models decreases the costs for choosing the optimal operating conditions in machining. It also avoids the expensive prototypes during the design stage of new spindles.
Angular contact ball bearings are mainly used for HSM spindles. Multi-Degrees Of Freedom (DOF) dynamic models of bearings are generally based on the work of Jones [1]. Stiffness matrix computation is also decisive for including the bearing model into the rotor model. Computation techniques are detailed in [2,3,4]. The structural model of spindle is mostly based on beam or 3D Finite Elements Models [5,3]. State of the art also include the thermal behavior [6].
Nevertheless, numerical models are not sufficiently widespread in the industry. Are models too complex or not enough refined? Do they not reflect the reality? The aim of the presented work is to combine experiments and simulations to build a model reflecting the real behavior. The analytical model is based on the structural setup. The updated parameters are related to the preload system and the ball bearings.
First, the 5DOF numerical model of the ball bearing is shortly presented. The axial model of rotor is detailed along with the solving method. Then, the experimental protocol is explained. A new experimental device for applying axial loads is introduced. Finally, two updating steps including sensitivity analysis are performed. As a result, the model has been enriched and crucial parameters such as preloads have been identified. The prediction of the behavior of the HSM spindle is particularly ac-

Numerical Modeling 2.1 Industrial System
This paper deals with a Fischer Electrospindle MFW 2310 (Speed 24 000 rpm, Power 70 kW ) for HSM milling applications. Its structure is presented in Fig. 1. Bearings are numbered 1 to 3 and ball bearings a to e. a,b and c ball bearings are SNFA VEX70/NS9CE3 and ball bearings d and e are SNFA VEX60/NS9CE3. The spindle is elastically preloaded with two sets of springs. The corresponding preloads will be referred as front and rear preload. The bearing 1 is preloaded by both systems. The motor is located between bearings 2 and 3. Ball bearings d and e are mounted in a sleeve that can translate in the spindle housing, in order to compensate thermal elongation of the rotor due to motor heat.

Ball Bearing Model
For a complete rotor model, a dynamic 5DOF ball bearing model is required. It consists in giving the constitutive relation between global displacements d = (δ x , δ y , δ z , θ y , θ z ) t and the global loads f = (F x , F y , F z , M y , M z ) t . Because no direct expression exists for a multi DOF model, the local behavior is expressed for each ball (see Fig. 2).
The first method, called analytical approach, solves the global displacements d from the global loads f. It is based on an hypothesis on the local load repartition [7,8] and can not take into account dynamic effects on balls (centrifugal force F c and gyroscopic moment M g ). This method is rejected for high speed application (i.e. high values of the d m N criterion).
The second method, called the numerical approach, gives the global loads f from the global displacements d. It is based on the three steps represented on Fig. 2 [1,3,9]. The linearized relation given by the stiffness matrix K is computed analytically as detailed in previous personal work [4]. The stiffness matrix is necessary for determining the global displacements d from the global loads f. Furthermore K is needed to compute the shaft equilibrium.

Axial Spindle Model
The numerical model of the rotor is based on Fig. 3. Because the aim is to identify the characteristics of preload setup and ball bearings, a rigid rotor with a quasi-static behavior is considered. The non-linear dynamic model of ball bearing described in section 2.2 is considered. For tandem setup of ball bearings, the load and the stiffness are doubled to take into account the two ball bearings (for instance for ball bearings a and b : The input of the model is the axial load on the shaft F and the preload parameters: preload (P 1 , P 2 ) and preload springs stiffness (K p1 , K p2 ). The results are the axial loads on bearings (F x1 , F x2 , F x3 ) and axial displacements u = (u, u p2 , u p1 ) of respectively the shaft, the outer ring of bearing 2 and the sleeve.
The model is based on three equations respectively corresponding to equilibriums of the shaft, of the third bearing's outer rings and of the outer ring of bearing 2: Free state Preloaded, no external load With axial load The numerical solving is done by a Newton-Raphson algorithm. The Jacobian matrix J is found by differentiating Eq. (1) with respect to the axial displacements u = (u, u p2 , u p1 ). For this purpose, the expressions for axial loads are linearized: with (δ x1,0 , δ x2,0 , δ x3,0 ) corresponding to the axial displacement of the ball bearings at the preload state. K xx is the linearized axial stiffness of the bearings. Because the ball bearing model is non-linear, its value depend on the axial load F x . It need to be computed at each step of the Nexton-Raphson algorithm.
At last, the solving algorithm is: with ξ the residual vector of equations Eq. (1) and: The algorithm is initialized by the preloaded state. At each step n of the algorithm, the ball bearing model is solved 3 times (i.e. for each bearings 1, 2 and 3).

Experimental Setup
An experimental device has been developed for applying axial loads in both directions on the spindle, whatever the shaft speed (see Fig. 4). Two inductive sensors are placed axially (see on   of the rotor under working conditions, i.e. at a given speed and considering a warm and thermically stable spindle. Indeed the aim is not to develop the complementary thermal model. That is why, before each test, a long warm-up has been carried out in order to always work with the same steady thermal conditions given by the CNC controller. Also, tests are done quickly to avoid major temperature change. For that mater, a relatively high sampling frequency is chosen for the sensors (12, 5 kHz). Repeatability has been verified with excellent results.
The data processing consists in a moving average to suppress the runout of the front face of the spindle. The data are calibrated considering (u, u p2 ) = (0, 0) for a warm state of the spindle, without shaft speed or axial load F.

Experimental Results
The experimental data after post processing are displayed in Fig. 5 and Fig. 6. The circle markers are the selected points for the updating.
With an increasing shaft speed, the axial displacements (u, u p2 ) can be observed on Fig. 6. They are mainly due to centrifugal forces on balls that tend to change the balls location. The displacement of the front of the spindle u is directly axial deflexion of ball bearings (a,b). The displacement of the rotor sleeve u p2 is the sum of deflexions of ball bearings (a,b) and (d,e).
On to the contact loss of the front bearings 1. The corresponding load F A is approximately the preload on the front bearing P 1 + P 2 = 1 200 N because the preload springs are soft (see discussion in section 6.1). The axial stiffness for F < F A is given by Hence the behavior of the spindle is quasi linear for F < F A because the ball bearings are a lot stiffer than the preload springs. The slope of the line is approximately the equivalent stiffness K eq = K p1 + K p2 = 12.0 N/µm. With nonzero shaft speeds, there is no more loss of contact due to the separation of contact angles in relation to dynamic effects on balls: α i on the inner ring and α o on the outer ring.

Model Updating (step 1)
In order to update the spindle model, a cost function is built. Two experimental data types are chosen. First, some displacements u in relation to axial loads F have been selected for N equals 4 000, 16 000 and 24 000 rpm (see blue circles on Fig. 5). Second, the deflexion of the rear sleeve u p2 during its rotation at 16 000 and 24 000 rpm have also been chosen (see blue circles on Fig. 6).
The cost function is: The number of experimental points (i + j) has been reduced to a minimum to ensure a reasonable computational cost.
To determined the parameters that need to be updated, a One-Factor-At-a-Time (OFAT) sensitivity analysis is realized. The parameters for the sensitivity analysis are given in Tab  of parameters are fixed: ±5% for material values, ±0.1% for geometrical values and α 0 = 25 •+3 −2 for the free state contact angle as provided by the bearing manufacturer. The preloads and preload stiffness are unknown at this stage. Nevertheless minimal preload is fixed as the recommended light preload: single row for the front preload P 1 = 100 N and tandem for the rear preload P 2 = 150 N. The maximal values are approximated from the experimental data at 4 000 rpm: P 1 + P 2 = 1 200 N and K p1 + K p2 = 12.0 N/µm, as explained in section 4.
The sensitivity values are expressed in percentages and have been computed using a 1% magnitude perturbation of the parameter's range. Results of the first sensitivity analysis are presented in Tab. 1. The five first parameters are retained for the updating: (P 1 , P 2 , K p1 , K p2 , α 0 ). The model updating has been conducted with an optimization Matlab algorithm: fmincon. The values minimizing the cost functions are presented in Tab. 1: (P 1 , P 2 , K p1 , K p2 , α 0 ) = (978, 150, 11.4, 5.15, 28). The simulations of the updated model after step 1 are displayed on Fig. 7 and compared to the experimental results.
Results are approximately in agreement at 4 000 and 16 000 rpm in spite of stiffness errors that can be observed by the slopes at null force F. But at high speed and large negative force, if the model forecasts a drop of stiffness, the experiments reveals an increase that is not compatible with the current model.

Stroke limit on Front Bearing
Model On the experimental results of Fig. 7, inflexion points B, B' and B" are observed. This behavior can not yet be described by the numerical model. The stiffness measured is higher before these points (F < (F B , F B , F B )) as would be a rigidly preloaded rotor. According to the authors, a physical limit could restrict the motion of the outer ring of the number 2 bearing. A stop can be considered and limits the motion of the outer ring in the left direction. A new constraint is established for numerical simulations: u p1 ≥ u p1,l with u p1,l < 0.
The solving algorithm is modified: if the numerical solution u p1 < u p1,l , the displacement u p1 is fixed at u p1,lim and the solution is recomputed again using the following algorithm: with I = J(1 : 2, 1 : 2).

Impact of the stroke limit on the behavior
To understand the effect of the stop, it is necessary to look over the behavior of the rigidly and elastically preloaded systems   FIGURE 9. Axial behavior of a rotor with elastic preload at various spindle speed. Simulation results are presented on Fig.  8 and Fig. 9 considering the spindle of Fig. 1. The solid lines correspond to the front bearing 1 and the dashed lines to the rears bearings 2 and 3. On Fig. 8, the curves of a rigidly preloaded system without shaft speed are well known in the literature (see examples in [10]). The load F can be found for a given displacement u by the equilibrium of the rotor: F = F 1 − (F 2 + F 3 ). As an example, the Fig. 8 displays loads corresponding to the shaft displacement u = −20 µm at 24 000 rpm. The preload value P 1 + P 2 is directly found in the intersection point of the curves F 1 and F 2 + F 3 . Thus, P 0 ,P 16 and P 24 are the preloads respectively for 0, 16 000 and 24 000 rpm. The preload increases significantly with shaft speed (see Fig. 8). In the example, the preload increases from 1 205 N to 3 770 N with spindle speed. This is one of the reasons why the rigid preload is put aside for high speed rotors. The points A corresponds to the loss of contact of bearing 1.
The simulations concerning the rotor with front and rear elastic preloads are presented in Fig. 9. Loads on the front bearing F 1 in relation to the rotor axial displacement u are equal to the rigidly preloaded arrangement. Besides, loads F 2 + F 3 evolve linearly. Indeed, the equivalent stiffness of bearings 2 and 3 is k eq = K xx2 K p1 K xx2 +K p1 + K xx3 K p2 K xx3 +K p2 . Because the stiffness of springs are much lower than the stiffness of ball bearings, the equivalent stiffness becomes k eq ∼ k p1 + k p2 . In that case, the equivalent stiffness does not evolves with axial load F. Hence, the curves F 2 + F 3 are lines. As a result, the elastically preloaded spindle stiffness K sp is far below the rigidly preloaded one. Contrary to the rigidly preloaded system, no major increase of the   FIGURE 10. Axial behavior of a rotor with elastic preload in presenceof a stop limit on bearing 2 preload P 1 + P 2 with shaft speed is observed. In the example, the preload at 0 rpm is P 0 = 1 205 N and evolves to P 24 = 1 650 N at 24 000 rpm.
If a stop is considered on bearing 2 of the arrangement, other inflexion points are observed: B,B' and B". Passed the boundary, the behavior of the front preload is rigid while the rear preload system (bearing 3) is still working elastically. The displacement corresponding to the boundary depends on the speed because of the axial displacement of each ball bearing due to the spindle speed.
The temperature difference between rotor and stator has a major influence on this boundary. Nevertheless, heat is not dangerous for the lifespan because it tends to reduce the rigid preload. The presence of this stop can be of a good advantage for the stiffness of the spindle. If it is well chosen, the spindle can work as a rigidly preloaded system at high speed.

Radial Ring Expansion
The first update was unsuccessful. The load F A corresponding to the lost of contact of the front bearing 1 and the stiffness after contact loss k eq were correct but the error was too important on the spindle stiffness K sp for F > 0 (see Fig. 7). A physical phenomenon is still missing in the model. In fact, the rings expansions due to shaft speed are not yet taken into account in the model. The rotating inner ring expands due to centrifugal forces and the immobile outer ring expand due to the excessive load of balls at high speed (i.e. centrifugal force on balls). Even if expansions are low, they are of major influence on the behavior, because ring curvature radius r i and r o are close to the ball A method for taking into account rings expansion is explained in [11]. The change of contact angle due to the rings radial expansion is given by: where ∆u = u i − u o is the difference of radial expansion between the inner and outer rings. One solution is to model accurately the ring expansion either with a Finite Element or continuum mechanics model. The other solution is to identify the evolution of ∆u with shaft speed from experimental data. Thus additional parameters are added to the update algorithm: (∆u 04 , ∆u 16 , ∆u 24 ).

Model Updating (step 2)
A new updating step has been conducted with the refined model. Four parameters are added: the end stop clearance on bearing 2 u p1,l , and the radial expansions of rings (∆u 04 , ∆u 16 , ∆u 24 ). The initial values are the values resulting from the first update in section 5.
A new sensitivity analysis has been conducted (see step 2 in Tab. 1 in Appendix A). It has led to the selection of the parameters of major significance for the update: (P 1 , P 2 , K p1 , K p2 , u p1,l , ∆u 04 , ∆u 16 , ∆u 24 ). The parameter α 0 has been ruled out because it is redundant with the expansion parameters ∆u (due to Eq. (8)). The selected experimental data are identical to step 1.
The result of the second update is displayed in Fig. 11. The mean deviation between simulation and experimental displace-ments decrease from 4, 37µm to 1, 03 µm, which proves the very good quality of the model. It is particularly accurate for medium shaft speed (4 000 to 16 000 rpm). The slight difference at high speed is probably due to the required kinematic hypothesis chose for computing the dynamic effects in ball bearings. This reason also explains the inaccuracy for F < −1 200 N. Fortunately the spindle is never used under such extreme loading conditions. In conclusion, model's parameters such as preloads and preload spring stiffness have been identified from experimental data. New refinement has been made with success on the final updated model. The axial model is finally validated.

Conclusion
For ensuring the accuracy of the HSM spindle model, an updating stage is undertaken. After presenting the ball bearings model and the rotor model, the experimental protocol has been exposed. Experimental tests have been done on a CNC machine with a new specific device to load the shaft in both directions whatever the spindle speed. Sensitivity analyses have enabled to select the relevant parameters for the updating process. The first step revealed unexpected physical phenomena in the model: an end stop limiting the motion of one of the bearings and speeddependent radial expansions of ball bearing rings. The impact of the end stop has been detailed and compared to the behavior of rotors with rigid and elastic preloads. The end stop is of great advantage regarding the spindle stiffness at high speed. A second updating step has been realized considering the rifined model. Preload forces, preload stiffness, ball bearing rings radial expansions and the end stop clearance have been identified. They correspond to the real operating conditions. The behavior obtained with the final updated numerical model describes with good agreement the reality. This step of model updating is crucial in the process of building a 3D dynamic and non-linear model of a high speed rotor.

Capital letters D
ball diameter E modulus of elasticity K bearing stiffness matrix K p preload spring stiffness K sp equivalent axial stiffness of the spindle K xx = K(1, 1) axial stiffness of the ball bearing N shaft speed P preload Q ball-raceway normal load Lowercase letters d = (δ x , δ y , δ z , θ y , θ z ) global displacement of inner ring. d m ball orbital diameter f = (F x , F y , F z , M y , M z ) t global load on inner ring (f and d are expressed at the center of the outer ring) f r/D r raceway groove curvature radius u = (u, u p2 , u p1 ) t axial displacement vector u axial displacement of the shaft u p1 axial displacement of the outer ring of bearing 2 u p1,l stop clearance for u p1 u p2 axial displacement of the rear sleeve (bearing 3) x axial direction y, z radial direction Greek letters α contact angle δ ball-raceway normal displacement ∆u N = u i − u o difference between radial ring expansions at speed N ε mean deviation between experiments and simulations (in microns) ν Poisson's ratio ξ residual vector of the equilibrium equations Eq. (1) ρ mass density Subscripts b ball i inner ring o outer ring r ring