David Tlalolini

Mathieu Ritou

Clément Rabréau

Sébastien Le Loch

Benoît Furet

Clément Rabreau
email: clement.rabreau@univ-nantes.fr

Benoit Furet

Modeling and characterization of an electromagnetic system for the estimation of Frequency Response Function of spindle

Keywords: Electromagnetic force, Frequency Response Function, Machine tool, spindle dynamics

Introduction

The dynamic behavior of spindles has a great influence on the quality of machined parts, notably in the aircraft industry where instabilities in milling can lead to nonquality and reduced lifetime of spindle bearings [START_REF] De Castelbajac | Monitoring of distributed defects on HSM spindle bearings[END_REF] and cutting tool [START_REF] Ritou | Angular approach combined to mechanical model for tool breakage detection by eddy current sensors[END_REF]. To improve the productivity of machine tools and to avoid the occurrence of chatter during milling, stability lobes diagrams are commonly used. The Frequency Response Function (FRF) at the tool tip is required and classically obtained by hammer impact test. It provides good results in classical milling. However, as shown in several studies [START_REF] Cao | A General Method for the Modeling of Spindle-Bearing Systems[END_REF][START_REF] Albrecht | High frequency bandwidth cutting force measurement in milling using capacitance displacement sensors[END_REF][START_REF] Rabréau | Phenomenological model of preloaded spindle behavior at high speed[END_REF] the dynamic behavior of High Speed Machining spindles varies at high speed due to dynamic effects. Schmitz et al. [START_REF] Schmitz | A method for predicting chatter stability for systems with speed-dependent spindle dynamics[END_REF] were the first to make measurements on a rotating spindle. Tap tests were performed at the tool tip, during spindle rotation. The displacement was measured with eddy-current sensors. In order to improve the measurement accuracy and repeatability, Albrecht et al. [START_REF] Albrecht | High frequency bandwidth cutting force measurement in milling using capacitance displacement sensors[END_REF] used a specific tool with a ball bearing to make hammer impact test measurement on a fixed surface. Active Magnetic Bearings (AMB) were then used to apply non contact excitation to the spindle [START_REF] Knospe | Active magnetic bearings for machining applications[END_REF][START_REF] Abele | Electromagnetic Actuator for in Process Non-Contact Identification of Spindle-Tool Frequency Response Functions[END_REF] and to measure the cutting forces [START_REF] Auchet | A new method of cutting force measurement based on command voltages of active electro-magnetic bearings[END_REF]. Rantatalo et al. [START_REF] Rantatalo | Milling machine spindle analysis using FEM and non-contact spindle excitation and response measurement[END_REF] have presented a contact-less dynamic spindle excitation device that can apply a dynamic force of ±15N from 400 to 2000Hz. Matsubara et al. [START_REF] Matsubara | Non-contact measurement of spindle stiffness by using magnetic loading device[END_REF] succeeded to apply a dynamic force of 110N from 0 to 1000Hz with their electromagnetic system, however it presents non-linearities. Strong requirement are required in order to analyze the dynamic behavior of high speed spindles.

This paper presents an excitation device capable of applying a dynamic load of ±80N at up to 5000Hz. A specific current amplifier has been designed for this purpose. The control is linearized with the use of pre-magnetisation and an experimental protocol to measure the flux density at a large range of frequency is proposed. Thus, the applied force can be precisely estimated over the entire frequency range. First of all, the system design is described, then electromagnetic force models are developed for the static and dynamic excitation modes. The specific experiments carried out to validate the model and to update the permeability of the magnetic core are described. The dynamic characteristics of the excitation device are then studied. Finally, FRF measurement are made with the excitation system on a HSM spindle at several speeds and are compared with a classic tap test experiment. The excitation system shown in figure 1, consists in a radial magnetic bearing with two perpendicular axes (x and y) and a dummy tool with an HSK 63A interface for the mounting into the spindle rotor. Each axis is composed of two identical electromagnets with magnetic field of opposing directions to produce a radial force. A pre-magnetisation is applied in order to improve the performance and linearize the control. The system is designed to provide a magnetic force with amplitudes until ±190N for the static mode and ±80N for the dynamic mode. Each electromagnet is made of two coils of wire wrapped around a core of ferromagnetic steel sheets. The coils are wired in parallel one to each other in order to improve dynamic efficiency. The magnetic core is made from laminated sheet of grain oriented electrical steel (3% Si-Fe) The core structure comprises 170 sheets with a thickness of 0.08mm which are assembled shifting each oriented sheet from the previous one by an angle of 90 • . A large cross section A e of the magnetic core was chosen. This configuration, studied in [START_REF] Lopez | Grain oriented steel assembly characterization for the development of high efficiency ac rotating electrical machines[END_REF], gives the best performances concerning the magnetizing current, the core loss amount, the relative permeability and the reactive power. A stack of ferromagnetic steel sheets is also integrated in the dummy tool to strengthen the magnetic attraction. In order to determine the analytic force model, a single two-pole electromagnet is considered, as presented in figure 2. The figure illustrates the geometrical parameters of the electromagnet, i.e. effective magnetic path length l e , air gap s and effective crosssection area A a in the air gap. Since the air gap lengths are sufficiently small compared with their cross-sectional dimensions, the fringing flux effect can be neglected. In this consideration, the magnetic flux φ is homogeneous in the core (c) and the air gaps (a), i.e. φ = B c A e = B a A a . Then, the flux density B is the same along the path l e in the magnetic core (B c ) and in the air gap (B a ), i.e. B a = B c = B. Thus, the cross-sectional areas are the same A a = A e . Under these assumptions, the magnetic system can be assumed lossless and the magnetic force can be derived from the energy W a stored in the two air gaps, by

W a = 1 2 BH a A e 2s (1)
where H a is the magnetic field intensity in one air gap, A e is the cross-sectional area and s is the air gap length. The flux density B is obtained through the permeability of free space µ 0 : B = µ 0 H a . Let's consider a radial displacement x of the shaft, along the radial force f . It involves a variation of the air gap of x.cos(α) (where α is the semi-angle between the two elementary forces, acting at each air gap). The force can be obtained from the principle of the virtual work dW = f dx. Assuming that B remains constant, the magnetic force f of the two-pole electromagnet derives from [START_REF] Schweitzer | Active magnetic bearings -chances and limitations[END_REF]:

f = dW dx = B 2 µ 0 A e cos(α) (2)

Static force model

Considering a single two-pole electromagnet (as shown in figure 2), the static force model is derived from the magnetomotive force equations (mmf) generated by a current i that is distributed in the two coils [START_REF] Balakrishnan | Air-gap reluctance and inductance calculations for magnetic circuits using a schwarz-christoffel transformation[END_REF]:

2N i = l e B µ 0 µ rs + 2sB µ 0 ( 3
)
where N is the number of wire turns of the electromagnet coil and µ rs is the relative permeability of the magnetic core for the static mode. The first term of right part of the equation (3) describes the mmf due to the magnetic core and the second term is the mmf due to the two air gaps. The flux density B is obtained by solving equation (3) in relation to the actual current i and the air gap s:

B = 2µ 0 N le µr s + 2s i (4)
As mentioned before, each axis is composed of two identical electromagnets in a differential configuration that generate the magnetic force. This configuration is depicted in figure 3. According to equation (4), the static magnetic force of the axis F s can then be obtained as the difference of the forces of both electromagnets as shown in the following equation:

F s = f + -f -= 4µ 0 N 2 A e cos(α)    i 2 + le µr s + 2s + 2 - i 2 - le µr s + 2s - 2    (5)
where i + and i -are the currents of each electromagnet. s + = s 0 + x and s -= s 0 -x are the respective instant air gaps for a radial displacement x of the dummy tool and a nominal air gap s 0 .

The magnetic force produced by each electromagnet shows the quadratic dependence on the current and the inversely quadratic dependence on the air gap. In order to linearize the control, a constant premagnetization current i 0 is applied. In addition, the control current i c , leads to i + = i 0 + i c and i -= i 0 -i c . Assuming that the system operates within its linear domain and that the displacement x of the rotor is much less than the nominal air gap s 0 , the static force can be simplified as:

F s = 16A e µ 0 N 2 cos(α) le µr s + 2s 0 2 i 0 .i c (6)
The value of µ rs is further obtained from experiments. The assumptions of linear domain and constant B during a small displacement x (for equation ( 2)) are verified in section 4.2.

f f x i + i + i - i -

Dynamic force model

For harmonic excitation, the control current is sinusoidal, i.e., i c = Î sin(ωt). The variation of flux density closely approximate a sinusoidal function throughout the length of the core. Assuming that displacement x is negligible (compared to s 0 ) and considering a premagnetization i 0 , the magnetic flux in equation ( 4) can be written as a function of the sum i 0 + i c and the difference i 0 -i c . Besides, the premagnetization flux B 0 and the dynamic flux Bd can be distinguished.

B = 2µ 0 N le µr + 2s 0 (i 0 ± i c ) = 2µ 0 N le µr s + 2s 0 i 0 ± 2µ 0 N le µr d + 2s 0 Î sin(2πf e t) (7)
B = B 0 + Bd sin(2πf e t) (8)
with

B 0 = 2µ 0 N le µr s + 2s 0 i 0 and Bd = 2µ 0 N le µr d + 2s 0 Î
For one axis, the magnetic force F d for dynamic excitation is then:

F d = 4A e cos(α) µ 0 B 0 Bd sin(2πf e t) (9)
µ r d is the dynamic relative permeability of the ferromagnetic core and depends on the pre-magnetisation current and on the frequency. The value is also different from the static µ rs .

The equation [START_REF] Knospe | Active magnetic bearings for machining applications[END_REF] enables the study of the behavior of the magnetic force under a sinusoidal excitation, but requires the flux density to be known. The determination of the flux density by means of an experimental test will be developed in the section 4.3. It will result in the identification of µ r d .

Experimental characterization of the electromagnetic system

In this section, the experimental procedure to characterize the magnetic flux density, both for the static and dynamic modes, is described. For the static model, the flux density is determined from equations ( 4) and ( 6) and from the measure of the magnetic force. The flux density for the dynamic model will be estimated with a specific experimental set-up and measurement of induced voltage.

Experimental set-up

The experimental set-up, illustrated in figure 4 and5, includes the excitation device mounted on a dynamometer table as well as the dummy tool clamped into the spindle rotor by an HSK 63A interface. The maximum spindle speed is 24000 rpm. The rotor is excited by the electromagnets with either Direct Current (DC) for quasi-static measurement or harmonic Alternating Currents (AC) for dynamic measurement. The excitation can be performed in both x and y radial directions. Due to the identical design of the x and y axes, the study is carried out in the x axis only.

The excitation current, for each axis, is provided by two Pulse Width Modulation (PWM) amplifiers, specifically designed for this application and can provide a maximum continuous current of 15A and 10A rms per electromagnet in a frequency range from 50 to 5000 Hz. In order to achieve this requirement, the amplifiers are designed with an internal current control loop. A non-contact displacement sensor is used to measure the instantaneous position of the rotor at the tip of the dummy tool. Other displacement sensors can be used to measure the displacement from another position of the dummy tool. For a DC excitation, the magnetic force is measured by a dynamometer table KISTLER 9255B. The National Instruments CompactDAQ system is used to interface the input and output analogic signals of the system at a sampling frequency f sp = 50kHz. A user interface was developed with LabVIEW software to manage the experiment. The measured datas are then converted and post-processed with Matlab.

Static characterization and model identification

A specific experiment was designed to characterize the static force behavior and to identify the relative permeability µ rs . The experiment is carried out for one electromagnet without pre-magnetisation. The control current i + is increased and the applied force is measured with the dynamometer table. The magnetic flux density B can then calculated from equation [START_REF] Albrecht | High frequency bandwidth cutting force measurement in milling using capacitance displacement sensors[END_REF].

From equation (4), B(i) of a single electromagnet is assumed to be linear, in a given domain. The magnetic force model, in the linear domain, of one electromagnet is determined from equations ( 2) and (4) : The value of the magnetic permeability is then updated from equation (10) using the parameters shown in Table 1 and the experimental results presented on figure 6a. The results presented are obtained with µ rs = 7432H/m.

f + = 4µ 0 N 2 A e cos(
The figure 6b present the experimental measurement of the magnetic flux density B and the linear approximation obtained with the updated model. Two domain can be identified corresponding to a linear evolution of the Magnetic Flux with the current and a saturated region when B > 0.94T . To stay in the linear domain, the applied current must be limited to i max = 4.6A in each electromagnet. 6b) and have opposed saturation domain, the current i 0 should be chosen in the middle of the electromagnets linear domain, i.e. i 0 = i max /2 = 3.5A

Dynamic characterization and model identification

One of the major issue to obtain FRF is to know the force level that is applied at a given frequency. Due to the low bandwidth frequency of the dynamometer table, this force level can't be measured at high frequencies (≤ 1kHz). To evaluate the magnetic force from equations ( 8) and ( 9), the amplitude of flux density as function of the excitation frequency f e must be obtained. An identical current is provided by a PWM amplifier, whatever the frequency. The link between the control current Î and the magnetic flux B is therefore required. To do so, a specific experimental setup is used [START_REF] Lamsont | A Method of Measuring the Magnetic Properties of Small Samples of Transformer Laminations[END_REF]. The operating flux density of one core can be measured through the induced voltage in the other coil. A digital oscilloscope is used to monitor and measure the RMS value of induced voltage as shown in figure 7. A sinusoidal current is supplied in the excitation coil while the induced voltage from the other coil is measured. Using Faraday's law of induction, the relation between the magnetic flux ϕ and the RMS of induced voltage V rms is given by:

ϕ = V rms √ 2πf e N ( 11
)
where ϕ is the magnetic flux produced by the control current supplied to one excitation coil. Since the two-pole electromagnet consist in two identical windings, the resultant magnetic flux throughout the length of the core is produced by the total mmf of the two windings. Assuming a linear flux-mmf relationship, the total magnetic flux is that is generated by the device in dynamic mode (see 3) is:

φ = 2ϕ = √ 2 πf e N V rms (12)
Under the assumption that the magnetic flux density is uniform across the cross section of the magnetic circuit, the magnetic flux is defined as φ = BA e . It results from equation ( 12) that the flux density can be estimated by:

B = √ 2 πf e N A e V rms (13)
The amplitude of flux in dynamic mode can hence be evaluated from equation ( 13) by the measure of induced voltage over the frequency range of 50 to 5000Hz. Figure 8 shows the amplitude of flux density versus excitation frequency for a given value of current Î. The optimal value of control current between the saturation boundaries of both electromagnets of an axis is Î = 1.5A, around the pre-magnetisation current i 0 = 3.5A (see figure 6b)

The eventual variation of th flux density in relation to the excitation is modeled through the dynamic relative permeability µ r d in equation ( 8). An analytical function of relative permeability as function of the frequency can be proposed to determine the amplitude of flux density from equation [START_REF] Abele | Electromagnetic Actuator for in Process Non-Contact Identification of Spindle-Tool Frequency Response Functions[END_REF]. A second order rational function is proposed, given by:

µ r d (f e ) = Σ 2 i=0 a i f i e Σ 2 i=0 b i f i e (14)
The coefficients a i and b i are obtained by curve fitting experimental data shown in figure 8a with equation [START_REF] Abele | Electromagnetic Actuator for in Process Non-Contact Identification of Spindle-Tool Frequency Response Functions[END_REF]. Simulation of flux density and relative permeability are performed with the identified coefficient and are presented in figure 8.

The identified model of the dynamic force, equation ( 9), is verified by comparing the measured and simulated frequency spectrum of the force over a frequency range of excitation from 50 to 1000Hz. The results illustrated in figure 9 validate the model update. It can be noticed that a dynamic force F d of 80N can be applied until 5kHz. The decrease of B after 7kHz is due to the bandwidth of the PWM.

Application of the excitation system

The excitation system has been designed to analyze the quasi-static and dynamic behavior of the spindle rotor, under different spindle speeds. In this section, firstly, the quasi-static excitation leads to the estimation of radial stiffness of the spindle (at the dummy tool tip). Particularly, the evolution of spindle stiffness with rotation speed can be observed. Then, swept sine excitation is applied for the dynamic mode and the FRF are estimated, through the dynamic force model. Lastly, the FRF evolution with spindle speed is observed.

Quasi-static results

The experimental set-up is the one presented in figures 3 and 4. The rotor is initially centered in the excitation system and the electromagnets are supplied by a programmed control current in order to apply a square signal load, as shown in figure 10a. A radial force F s is applied to the dummy tool at different rotating speed. The resulting displacement and the applied load are measured. In order to obtain a maximum force amplitude, the pre-magnetization current is i 0 = 3.5A and maximal control current is i cmax = 3.5A.

In this experiment, the saturated zone presented in figure 6 is reached. Therefore, the linear control of the radial force with the control current is not verified. The static force model was developed for the linear domain and cannot be used here, out of the linear domain. Consequently, the force is directly measured with the dynamometer table and not estimated from the current measurement.

The force and displacement measurements data are filtered to remove the signals noise. The selected cutoff frequency is low and thus eliminates the tool runout noise and high frequency noise from the displacement signal. The filtered signal is shown for four different speeds in figure 10b. The assumption of a negligible displacement of the dummy-tool compared to the initial air-gap s 0 = 250µm is verified from the quasi-static measurement.

As a result, the evolution of the radial stiffness of the spindle is obtained from the quasi-static test presented (see figure 11). It is mainly due to the centrifugal forces on bearing balls that modifies the contact angle and radial stiffness.

Dynamic validation

The dynamic validation of the excitation system consists in accessing the FRF at the dummy tool tip by exciting the rotor over a frequency range of 50 to 5000Hz. The experimental setup corresponds to figures 3 and 4. The FRFs are estimated for different spindle speeds. The magnetic force applied to the rotor in the x-direction is evaluated from equation ( 9) and from the measure of the control current i c .

Swept sine control current

In order to obtain a FRF over a specified frequency range, the excitation current is carried out by a sine sweep with a linearly varying frequency. The control current is described by: where f min and f max are respectively the minimum and maximum frequencies and T is the duration of the dynamic test in seconds. The time series of a swept sine between 0 and 10Hz are shown in figure 12a. The frequency spectrum of the signal is given in figure 12b. As shown in figure 12b, a swept sine solicitation excites all frequencies between f min and f max . The experimental dynamic test of spindles are performed with a sweep duration of 30 seconds in the frequency range 50-4000Hz.

i c = Î sin 2πt f max -f min 2T t + f min (15)

Signal processing

A specific signal processing is necessary for measured displacements. After a low-pass filter to suppress the high frequency noise, it is necessary to suppress the displacements corresponding to the tool runout. This is observed in the frequency domain by contributions at harmonics of the shaft frequency (see raw signal in figure 13). The filtered signal is obtained after the use of an H1 formulation to compute the FRF (i.e. cross spectrum of the displacement and load/autospectrum of the load). This formulation is used to reduce the noise on the displacement signal. The Matlab smooth function with the "rloess" option is then used to suppress the undesired harmonic content of the frequency domain.

Experimental dynamic results

Dynamic measurements were performed at different spindle speeds. Figure 14 shows the experimental results after post-processing. The modes frequencies and amplitudes varies with spindle speeds, which justify the necessity to use such a device. It is due to the evolution of the dynamic behavior of the spindle at high speed. A validation tap test was performed to confirm the FRF prediction made with the electromagnetic excitation device. It is plotted in dashed line in the figure and it is in good agreement with the FRF measurement made with the electromagnetic excitation device at 0 rpm. Discrepancies can be seen at high frequencies (higher than 2kHz) and are due to limitations of tap test configuration.

Conclusion

A specific device has been developed to analyze the quasi-static and dynamic behavior of HSM spindles. Force models are developed both for static and dynamic modes. Besides, a method to characterize the electromagnetic system behavior in static and dynamic mode is presented in this paper. The limit between the linear and saturation domain of the electromagnets has been experimentally identified; as well as the static and dynamic relative permeabilities that are required for the force models. Maximum forces of ±190N for the static mode and ±80N for the dynamic mode up to a frequency of 5kHz were obtained. Quasi-static and dynamic measurement were performed at up to 24000rpm on an industrial spindle to illustrate the system performance. The approach was also validated by a tap test measurement of FRF. Consequently, the electromagnetic system, remains a highly efficient tool to analyze the dynamic behavior of a machine-tool spindle.

13

Figure 1 :

 1 Figure 1: (a) The arrangement of the two pairs of electromagnets. (b) Sectional view of the excitation system and dummy tool.

3 .Figure 2 :

 32 Figure 2: A two-coils electromagnet with magnetic cores made of stack of ferromagnetic steel sheets.

Figure 3 :

 3 Figure 3: Specific arrangement of two electromagnets in order to achieve a linear control.

Figure 5 :

 5 Figure 5: Experimental setup diagram

F

  max = 130N in linear and 190N in the saturation domain. Since the signals of the two electromagnets are in phase opposition around a pre-magnetisation current i 0 (see figure

Figure 6 :

 6 Figure 6: Experimental and simulated behavior of only one electromagnet, in relation to the coil current i + , and identification of a limited linear domain of magnetic flux.

Figure 7 :

 7 Figure 7: Experimental set-up for the determination of the magnetic flux in dynamic mode.

Figure 8 :

 8 Figure 8: (a) Estimation of the flux density over frequency range of 50 to 5000Hz for a control current Î = 1.5A. (b) Model of the dynamic relative permeability µr d .

Figure 9 :

 9 Figure 9: Comparison of measured and model-based estimated dynamic forces.

  Excitation currents and radial force.

  Displacements at different spindle speeds N [rpm].

Figure 10 :

 10 Figure 10: Measurements in quasi-static mode.

Figure 11 :

 11 Figure 11: Variation of the spindle radial stiffness (at tool tip) with rotation speed.

Figure 12 :

 12 Figure 12: Example of a swept sine excitation from 1 to 10Hz for Î = 1A.

Figure 13 :

 13 Figure 13: Estimated FRFs in frequency domain with and without post-processing.

Figure 14 :

 14 Figure 14: FRF at the dummy tool tip in relation to spindle speed N [rpm].

Table 1 :

 1 Parameter values of the two-pole electromagnet

	Term	Symbol	Value	Unit
	Effective cross-sectional area	A e	182.2	mm 2
	Effective magnetic path length	l e	171.4	mm
	Permeability of free space	µ 0	4π × 10 -7 H/m
	Number of turns on a coil	N	45
	Nominal air gap	s 0	0.25	mm
	Pole angle	α	π/8	rad
	As a consequence of this identification, the maximal reachable forces are obtained.

Acknowledgments

The research was conducted within the context of the French FUI project "QuaUsi". The authors would like to thank the company Europe Technologies and Precise where some experiments were made.


A Boudelier

M M Ritou

S Garnier

B Furet

Cutting force model for machining of CFRP laminate with diamond abrasive cutter

Keywords: cutting force model, composite, trimming, multidirectional CFRP laminate, diamond abrasive cutter, cutting mechanism

The article presents a cutting force model for trimming operations of CFRP laminate with diamond abrasive cutters. Those tools are more and more encountered on industrial applications of CFRP trimming, due to their abrasion resistance and their low cost. Contrary to endmills, they consist of a large number of cutting grits, randomly distributed around the tool. To tackle the issue, a continuous model of tool engagement is proposed. Validity of the approach is verified. A mechanical model of cutting forces, adapted to CFRP laminate, is then presented. The evolution of specific cutting coefficient in relation to fibres orientation is investigated through a piecewise constant model. It leads to the proposal of a sine model for the specific cutting coefficients. The simulated forces are in good agreement with the experimental results of cutting tests, carried out in multidirectional CFRP laminate for different fibres orientation and widths of cut. Cutting mechanisms are finally discussed depending on fibres orientation.

INTRODUCTION

Machining of composite laminates is often difficult because of their mechanical and thermal properties: heterogeneity, anisotropy and low thermal conductivity [START_REF] Teti | Machining of Composite Materials[END_REF].

Moreover, their high abrasiveness, notably mentioned by [START_REF] Ramkumar | Effect of workpiece vibration on drilling of GFRP laminates[END_REF] in drilling operations and [START_REF] Dumas | Development of a robotic cell for trimming of composite parts[END_REF] in milling operations, conducts to premature tool wear.

That is why diamond abrasive cutters, constituted of small diamond grits (Figure 1), have been specifically developed for CFRP trimming operations. In comparison with carbide or polycrystalline diamond tools (PCD), tool life and productivity are improved while material integrity is respected. Moreover, costs are reduced.

Very few studies concern diamond abrasive cutters. Colligan and Ramulu (1999) characterised trimmed surfaces depending on tool grits size and feedrate. They also performed preliminary tests to study influence of these parameters on cutting forces. [START_REF] Boudelier | Investigation of CFRP machining with diamond abrasive cutters[END_REF] studied the CFRP machinability with diamond abrasive cutters. However, no cutting force model has been proposed for this tool technology whereas it is necessary in order to optimize cutting conditions and avoid material integrity defects. [START_REF] Slamani | A study of the combined effects of machining parameters on cutting force components during high speed robotic trimming of CFRPs[END_REF] studied the cutting force during the robotized trimming of CFRP.

Force models have been proposed for orthogonal cutting with carbide [START_REF] Bhatnagar | On the machining of fiber reinforced plastic (FRP) composite laminates[END_REF] and PCD tools (Wang et al., 1995a;[START_REF] Ramulu | Machining and surface integrity of fibre-reinforced plastic composites[END_REF]. Influence of fibres orientation, feed rate and cutting speed has essentially been studied. Empirical (Wang, 1995b), numerical [START_REF] Arola | Finite element modeling of edge trimming fiber reinforced plastics[END_REF][START_REF] Nayak | Machining studies of UD-GFRP composites part 2: Finite element analysis[END_REF][START_REF] Venu Gopala Rao | Three-Dimensional Macro-Mechanical Finite Element Model for Machining of Unidirectional-Fiber Reinforced Polymer Composites[END_REF] and analytical models have been developed.

The latter are based on a tangential -radial cutting force approach. At a macroscopic scale, the influence of fibres orientation can be integrated into the specific cutting coefficients through sine [START_REF] Sheikh-Ahmad | Multiple regression and committee neural network force predition models in milling FRP[END_REF][START_REF] Bonnet | CFRP drilling model: fiber orientation influence on mechanical load and delamination[END_REF][START_REF] Karpat | Mechanistic force modelling for milling of unidirectionnal carbon fiber reinforced polymer laminates[END_REF] or polynomial models [START_REF] Sheikh-Ahmad | Machining of Polymer Composites[END_REF]. At a microscopic scale, approaches based on the chip formation were developed, considering fibres breaking mechanisms (Bhatnagar et al., 3 1995;[START_REF] Puw | Chip formation model of cutting fiber-reinforced plastics perpendicular to fiber axis[END_REF][START_REF] Zhang | A force prediction model for cutting unidirectional fibre-reinforced plastics[END_REF]. Also, [START_REF] Yuan | Development of a cutting force prediction model based on brittle fracture for carbon fiber reinforced polymers for rotary ultrasonic drilling[END_REF] proposed a force model for the ultrasonic machining of CFRP.

In this paper, an analytical model of cutting force for trimming operations of multidirectionnal CFRP laminate with diamond abrasive cutters is presented. Contrary to endmills, diamond cutters consist of a multitude of grits with irregular shapes that are randomly distributed around the tool [START_REF] Malkin | Grinding Technology: Theory and Applications of Machining with Abrasives[END_REF]. That is why a continuous model of tool engagement is proposed. Secondly, as the influence of fibres orientation on the specific cutting coefficients is unknown, a first investigation is carried out based on a piecewise ( ) ( ). . sin( ).

t t z p r r z p F K f a F K f a         (1)
where K t and K r are the tangential and radial specific cutting coefficients (MPa), f z is the feed per tooth (mm/rev/tooth), a p the depth of cut (mm), φ the tooth angular position.

In the case of machining of CFRP laminates, the specific cutting coefficients vary with fibre orientation θ that corresponds to the angle between cutting speed direction and fibre orientation of the considered unidirectional ply in the laminate (Figure 2). During tool revolutions, it evolves with the tool angular position.

Continuous model of tool engagement

Considering trimming operations with diamond abrasive cutters, cutting models must be adapted. Indeed, a diamond abrasive cutter is close to a grinding tool but with a larger grit size. Contrary to tools with cutting edges, it relies on a large number of grits, with a random distribution of height and position around the tool body (Figure 1). From one tool to another, the distribution is different. Contrary to endmills, multiplicity and randomness of cutting grits lead to constant cutting forces during machining. Consequently, the classical approach with a localized engagement of each cutting edge is not appropriated.

5 N a e V f ψ θ x y χ F re F te χ* dχ f F n F f

Figure 3. Elementary cutting forces applied on a dχ sector

The assumption is made that the grits distribution around the tool is homogenous and uniform. This is the reason why a continuous tool engagement model is proposed in the paper.

At each instant, it is considered that the diamond abrasive cutter generates cutting forces all along the tool engagement domain. The domain is discretized into elementary sectors dχ (Figure 3). In order to obtain the elementary tangential F te (χ)

K F f a K F f a         ( 2
)
where f is the feed per revolution (mm/rev).

Coherence between the two approaches

The objective of the paragraph is to demonstrate that K t (θ), which has been introduced in the expression of elementary tangential force with diamond abrasive cutter approach (Equation 2), is similar to the one considered with endmills (Equation 1). To do so, the mechanical works of the spindle during a revolution are compared for the two tool engagement approaches.

6

In case of machining with endmills, the spindle work W c during one tool revolution is the sum of the work of each cutting edge j. It can be obtained by summing the instantaneous spindle power during one tool revolution period:

1 0 .. T Z c c t j W V F dt     (3)
where T is the tool revolution period, V c the tool cutting speed and R the tool radius.

Considering that V c =R.dφ/dt, and introducing Equation ( 1) in (3), the spindle work can be expressed as:

* . . . ( ). sin( ). c p t W R f a K d        (4)
In case of machining with diamond abrasive cutter, the instantaneous spindle power P d resulting from the elementary tangential forces F te (χ) can be expressed as: * . ( ).

d c t e P V F d      (5)
From Equation (2), it follows that: * () . . .

. sin( ). 2.

t d c p K P V f a d        (6)
As the tool is considered as continuously engaged (Figure 3), the spindle power is constant. Thus, the spindle work for abrasive cutter W d is simply given by: 0 ..

T d d d W P dt P T   (7)
As the period T=2πR/V c , the work can be written as: * . . ( ). sin( ).

d p t W R f a K d        (8) 7
The analytical expressions of the spindle work (Equations 4 and 8) are identical for the two tool engagement models. Consequently, the specific cutting coefficients that will be identified in the paper can be compared with ones obtained for other cutting tool technologies.

CUTTING FORCE MODEL FOR MULTIDIRECTIONAL LAMINATE

The CFRP multidirectional laminates of aeronautic structure consist in a stack of plies with, generally, unidirectional carbon fibres (see Figure 4); according to the laminate lay-up that defines the sequence of ply orientations (see lay-up in Table 1). ψ represents the fibre orientation of a given ply, in relation to feed direction.

Let's consider the machining of a given ply with a diamond abrasive cutter, the local depth of cut corresponds to the ply thickness e. The feed force F f (ψ) and the normal force F n (ψ) can be obtained by summing the elementary forces along the tool engagement domain.

After projection in the feed reference system (x,y), they can be expressed by the following expressions:

* ( ) cos sin ( ) . . sin( ). . . ( ) sin cos ( ) 2. ft nr FK f ed FK                                  (9)
where ψ represents the fibre orientation in relation to feed direction. θ is the relative orientation between fibre and cutting direction:

  2; 2           (10)
In order to determine the average cutting forces generated during the trimming of a multidirectional laminate, a mixture approach is used (Wang, 1995b;[START_REF] Sheikh-Ahmad | Multiple regression and committee neural network force predition models in milling FRP[END_REF][START_REF] Sheikh-Ahmad | Machining of Polymer Composites[END_REF]. It assumes that the interactions between plies are neglected.

Indeed, the interface between two plies is only composed of polymer matrix. Its mechanical behaviour during trimming is negligible in comparison with plies one's. This hypothesis will be validated in a following section. Consequently, the cutting forces can be simply obtained 8 by summing the forces necessary to machine each unidirectional ply of the laminate (Figure 4):

* 1 cos sin ( ) . . . sin( ). . . sin cos ( ) 2. i m ft i nr i FK f n e d FK                                   ( 11
)
where m is the number of different ply orientations, n ψi the number of ply oriented at an angle of ψ i degrees and e i the ply thickness, for a given orientation i.

MODEL OF SPECIFIC CUTTING COEFFICIENTS

At this step of the study, the evolution of the specific cutting coefficients in relation to the relative fibre orientation θ is unknown. This is the reason why the expression of K t (θ) and K r (θ) are, in a first approach, modelled by a piecewise constant function in order to observe their evolution.

Investigation of the influence of fibres orientation on specific cutting coefficients

The relative fibre orientation θ varies from -π/2 to + π/2. It is divided into n z different sectors of β z degrees. The specific cutting coefficients K t (θ) and K r (θ) are assumed constant for each sector. As the feed direction and the ply orientation ψ remain identical during the

. . . ( ). sin( ). cos( ). 2. (

). sin ( ).

z z z ze z z z ze n m j f i i t i j j i j n n j r i j j jn f F n e K d Kd                                (12)
where θ ij is the relative fibre orientation of an angular sector j, for a ply orientation i; and n ze is the index of the first cutting sector of the tool, evaluated by the floor function of n z .χ*/π+1. In this way, the Feed and Normal force models lead to a linear combination of the specific cutting coefficients, which will simplify their identification. Since the radial and tangential specific cutting coefficients are represented by a piecewise constant function, Equation ( 12) can be expressed as a matrix product of the specific cutting coefficients K:

2 ,2 2 ,1 2 ,1 . , cut z cut z f t n n n r n n F K F A K wit h F K and A F K                (13) 10
where n cut is the number of cutting tests. Therefore, the specific cutting coefficients can be identified by matrix inversion: K=A -1 .F.

Experimental setup

An aeronautic CFRP laminate of T800-M21 was trimmed to carry out the experiments. The laminate was composed of 37 plies of unidirectional graphite fibre fabric, made by alternating four fibres orientations (Table 1). The matrix used was a thermosetting epoxy resin M21.

The diamond abrasive cutter used during the experiments was composed of natural diamond grits fixed with a nickel bond on a cylindrical body by electro deposition (Figure 1).

A 16 mm tool diameter with 852 μm grits size was used. Widths of cut of 16 mm (slot) and 5.94 mm (in down-milling) were machined. For a discretization of the tool engagement in 12 sectors (n z =12) in the piecewise model (Figure 5), these two widths of cut have solicited respectively 12 (χ*=0°) or 5 tool sectors (χ*=105°). For each width of cut, 12 machining directions from ψ=0° to 180° were studied. F n are evaluated and resultant of cutting forces F is calculated.

All the experiments were carried out with external and internal water-based coolant, and industrial HSM cutting conditions were used, as defined in Table 2.

Results of the investigation

The specific cutting coefficients of the piecewise model have been identified by matrix inversion (Equation 12) and based on the feed forces F f and normal forces F n of the n cut =24 cuts of the experiments. Results are illustrated in Figure 6.

ψ n z =12 β z =15°a verage a b n z =12 β z =15°a verage K t (θ) (Mpa) K r (θ) (Mpa) θ (deg) θ (deg)

Specific cutting coefficients model

Due to the previous results (Figure 6), an analytical model of the specific cutting coefficients in relation to fibres orientation is proposed:

( ) . sin(2. ) ( ) . sin(2. ) t t m t a t r rm ra r K K K K K K             ( 14
)
where K tm and K rm are the mean values, K ta and K ra their amplitude and ψ t and ψ r their phase angle. In accordance with previous investigation, a 2θ appears since θ is π-periodic. Finally, the expressions of cutting forces (11) can be expressed as:

F K K f n e d F K K                                        (15)
In order to simplify the identification of model parameters, the model can be linearized:

0 1 2 0 1 2 ( ) . sin(2. ) . cos(2. ) ( ) . sin(2. ) . cos(2. ) t r K a a a K b b b             (
m f i i o i o f F n e a d a d a d b d b d b d                                             (18)
It can be expressed as a matrix product:

2 ,1 2 ,6 6,1 . cut cut n n F A K  ( 19
)
with F is the matrix containing the experimental cutting forces of the n cut experiments, A the matrix representing chip sections, and K the matrix containing the specific cutting coefficients

(a 0 ,a 1 ,a 2 ,b 0 ,b 1 ,b 2 )
. These parameters of the analytical model can hence be identified by matrix inversion (Equation 19). 2, i.e. 24 couples of feed/normal forces, have been considered in order to identify the parameters of the specific cutting coefficients model 14 (Equation 14) through matrix system inversion (Equation 19). Results are presented in Table 3 and illustrated in Figure 7.

Table 3. Identified parameters of the specific cutting coefficient model

Maximum tangential and radial specific cutting coefficients are obtained respectively for a fibre orientation of -54° and -40° whereas their minimum corresponds to an orientation of 36° and 50°. Considering the identified model of specific cutting coefficient, feed and normal forces have been calculated (Equation 15). Simulated cutting forces are compared with experiments in Figure 8.

Ψ (deg) F f (N) F n (N) F n (N) Ψ (deg) Ψ (deg) Ψ (deg)
F f (N) F n (N) F n (N) Ψ (deg) Ψ (deg) Ψ (deg)
The prediction of feed and normal forces is in accordance with the experimental data.

The relative errors are respectively estimated at ± 7.8% and ± 5.8% (or ± 39 N and ± 24 N)

for machining operation with a low width of cut (5.94 mm) whereas they are evaluated at ± 3.2% and ± 6.1% (± 46 N and ± 39 N) for a 16 mm tool engagement.

DISCUSSION

First point of discussion concerns the levels of specific cutting coefficients. During machining with diamond abrasive cutters, and contrary to machining with endmills, the radial component is more important than the tangential one (Figure 7). Tool geometry explains this trend. Indeed, rake angles are generally positive for milling with carbide or PCD tools (from 0 to 18°) (Wang, 1995a;[START_REF] Bhatnagar | On the machining of fiber reinforced plastic (FRP) composite laminates[END_REF]1995). Consequently, cutting configuration is favourable and specific cutting coefficients associated to the couple Tool/Material are low.

During machining operation with diamond abrasive cutters, cutting is less efficient. Indeed, the rake angles are highly negatives because grits have rather a spherical shape. By tool scanning, [START_REF] Boudelier | Optimisation of Process Parameters in CFRP Machining with Diamond Abrasive Cutters[END_REF] established that the average rake angle is -65° (Figure 9).

Consequently, level of specific cutting coefficients is higher. In order to explain the influence of fibres orientation on the variations of specific cutting coefficients, results have been transposed to an orthogonal cutting configuration (Figure 10). The grit is represented with the average rake angle. A chip section of 1 mm 2 has been considered (f=0.1 mm/rev, a p =10 mm). Tangential F t , radial F r and resultant F res forces have been simulated depending on fibres orientation, using the previously identified parameters (Table 3). Then, the resultant force F res have been projected onto the fibres reference system (x fib ; y fib ). Two components are obtained: F axi oriented in the fibres direction and F trans orthogonal to the fibres (Figure 10). Results are shown in the Figure 11. 17

F (N) F axi F t F res F trans F r

Figure 11. Cutting forces in orthogonal cutting configuration

Firstly, the lowest resultant force is obtained for a fibers orientation of about 60°. In that case, the grit cutting face is nearly perpendicular to fibres axis (Figure 10a). Transverse force F trans is negligible and only an axial force F axi is generated (Figure 11). compressive load is applied along fibres direction and a significant shear at the fibre-matrix interface appears. It is very favourable to crack initiation and propagation. Indeed, the fibre-matrix interface is considered as the weak point of polymer matrix composites due to their low interlaminar shear strength [START_REF] Sheikh-Ahmad | Machining of Polymer Composites[END_REF]. As resultant force is minimal, this configuration of fibres orientation is optimum for cutting operation.

The highest resultant force is obtained for fibres orientation around -30°. In that case, the grit cutting face is nearly parallel to fibres axis (Figure 10b). In this configuration, a significant transverse force F trans is generated, perpendicularly to the grit cutting face. Fibers are compressed the one on the other. Besides, axial force F axi applied on fibres is very low. That leads to very low shear at the fibre-matrix interface. As the weakness of the material is lowly solicited, this configuration is unfavorable for the cut.

CONCLUSIONS

A cutting force model dedicated to the machining of multidirectional CFRP laminate with diamond abrasive cutters was proposed. The analytical approach is based on the 18 variation of specific cutting coefficients depending on fibres orientation. The following conclusions can be drawn:

1. The model allows accurate prediction of cutting forces during trimming of multidirectional CFRP with diamond abrasive tools, over the entire range of fibres orientation (-90° to +90°) and for different widths of cut. Indeed, the simulated forces were in good agreement with the experimental results. Model error was evaluated at ± 7.8% for feed force and ± 5.8% for normal force.

2. An original approach was proposed in order to consider the large number of grits that are randomly distributed around the tool. Conversely to the local engagement of endmill cutting edge, a continuous engagement of the tool was introduced. Coherence with endmill model was verified. This approach represents a significant improvement over previous models and could be extended to grinding applications.

3. Variations of specific cutting coefficients depending on fibres orientation were investigated using a piecewise constant model. It led to the proposal of a sine model for the specific cutting coefficients. Combined with a mixture approach, accurate prediction of cutting forces during trimming operation of multidirectional laminate was obtained.

4. Results of the model were transposed to an orthogonal cutting operation and cutting forces were expressed in the fibres reference system. Cutting mechanisms were then discussed depending on fibres orientation in order to explain the variation of specific cutting coefficients.

Figure 1 .

 1 constant model. The results lead to the proposition of a sine model of specific cutting coefficients. Then, accuracy of the model is evaluated and compared to experiments carried out over the entire range of fibres orientation (-90° to +90°). Finally, results of the model are expressed in the fibres reference system in order to discuss cutting mechanisms depending on fibres orientation. Diamond abrasive cutter with grits size of 852 μm TOOL ENGAGEMENT MODELS Cutting edge model Classical milling cutters rely on a given number of cutting edges. Tool engagement is modelled by localized engagements of cutting edges. Then, cutting forces generated by each 4 cutting edge during material removal are proportional to chip thickness. Tangential and radial forces, applied by tool on material, are respectively defined by the following expressions:

Figure 2 .

 2 Figure 2. Cutting forces at the cutting edge

Figure 4 .

 4 Figure 4. Example of configuration of a multidirectional laminate,

  the relative fibre orientation θ also remains identical considering a given sector of the tool (Figure5).

Figure 5 .

 5 Figure 5. Piecewise constant model

Figure 6 .

 6 Figure 6. Evolution of the specific cutting coefficients (a) K t (θ) and (b) K r (θ) identified through a piecewise constant model

Figure 7 .

 7 Figure 7. Analytical model of the specific cutting coefficients in relation to fibre orientation, (a) K t (θ) and (b) K r (θ), identified from the experiments. All the experiments presented inTable 2, i.e. 24 couples of feed/normal forces, have

Figure 8 .

 8 Figure 8. Experimental and simulated feed forces F f (a-c) and normal forces F n (b-d), a e = 16 mm (a-b) and a e = 5.94 mm (c-d)

Figure 9 .Figure 10 .

 910 Figure 9. (a) Definition of grits rake angle, (b) Distribution of grits rake angle on a diamond

Table 1 .

 1 Mechanical properties of carbon/epoxy material T800-M21(Hexcel)

	Number of plies	37
	Ply thickness	0.27 mm
	Fibre volume fraction	56.6%
	Laminate Density	1580 kg/m 3
	Compression modulus	136 GPa
	Tensile modulus	172 GPa
	Lay up	(90/-45/0/45/0/45/0/-45/0/0/0/0/-45/90/90/45/0/90/90/0/45/90/
		-45/90/-45/0/0/0/0/-45/0/45/0/45/0/-45/90)

Table 2 .

 2 Cutting conditions

	Cutting speed	1400 m/min
	Spindle speed	27850 RPM
	Feed per revolution	0.25 mm/rev
	Depth of cut	9.96 mm
	Widths of cut	16 and 5.94 mm

ACKNOWLEDGMENTS:

This work was carried out within the context of the working group Manufacturing'21 which gathers 16 French research laboratories. The topics approached are: modelling of the manufacturing processes, virtual machining, and emergence of new manufacturing methods.


Mathieu Ritou

Clément Rabréau

Sébastien Le Loch

Benoît Furet

Didier Dumur

Benoit Furet

Influence of spindle condition on the dynamic behavior

Keywords: Spindle, Wear, Dynamics

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Cutting dynamics is of major importance for the optimization of High Speed Machining (HSM). Moreover, the efficiency of the spindle is also a key point of the machine-tool performance [START_REF] Abele | Machine tool spindle units[END_REF]. In particular, spindle dynamics at high speed is of major interest and only partly known. Several electromagnetic devices were developed to investigate experimentally the dynamic behavior [START_REF] Tlalolini | Modeling and characterization of an electromagnetic system for the estimation of Frequency Response Function of spindle[END_REF][START_REF] Matsubara | Non-contact measurement of spindle stiffness by using magnetic loading device[END_REF][START_REF] Rantatalo | Milling machine spindle analysis using FEM and non-contact spindle excitation and response measurement[END_REF]. Numerical mechanical models were proposed, based on accurate bearing models that consider the dynamic effects [START_REF] Cao | Mechanical model development of rolling bearing-rotor systems: A review[END_REF]. Classical updatings (identification of model parameters) of rotor dynamics model fail, due to too few possible measuring points (compared to classical structural dynamics) and too complex evolutions at high speed. More detailed and phenomenological approaches can provide a better understanding [START_REF] Cao | A General Method for the Modeling of Spindle-Bearing Systems[END_REF][START_REF] Rabréau | Phenomenological model of preloaded spindle behavior at high speed[END_REF][START_REF] Jedrzejewski | Modelling of angular contact ball bearings and axial displacements for high-speed spindles[END_REF]. Although the impact of, notably, spindle speed or preload on the Frequency Response Function (FRF) were studied [START_REF] Cao | A comparative study on the dynamics of high speed spindles with respect to different preload mechanisms[END_REF][START_REF] Zivkovic | Mathematical modeling and experimental testing of high-speed spindle behavior[END_REF][START_REF] Ozturk | Investigation of spindle bearing preload on dynamics and stability limit in milling[END_REF]; the impact of spindle condition on dynamics has been very rarely examined up to now. HSM spindles are very efficient but also present a high risk of spindle failure [START_REF] De Castelbajac | Monitoring of distributed defects on HSM spindle bearings[END_REF]. Due to high cost of curative maintenance, spindles are used as long as possible, before their repairs. In this way, partly damaged spindle can produce for months in industry, whereas it has an unknown impact on the cutting dynamics. Besides, it is arduous to perceive since very long periods must be considered. Therefore, there is a strong interest in the investigation of the impact of spindle condition on the dynamic behavior and machining performances.

The paper proposes an original experimental approach based on accurate characterizations of the dynamic behavior for different spindle conditions (repaired or damaged). The FRFs of the spindle were obtained for different spindle conditions and speeds, through an electromagnetic excitation device. The axial behavior was also investigated. The main failure modes of HSM spindle are identified and modeled. Their effects on FRF are simulated and conduct to general conclusions concerning their impacts on spindle dynamics. Lastly, cutting tests reveals the impact of spindle condition on stability lobe diagrams.

Experiments

Procedure

Seven experimental campaigns have been carried out in order to characterize accurately the evolutions of the dynamic behavior of a spindle, in relation to the spindle condition. Four repaired conditions (noted Rep) and three damaged conditions (noted Dam, just before unmounting for maintenance) of a unique spindle have been studied, over a period of more than three years (Fig. 1). There were industrial productions between them. The damaged conditions were determined by the company, notably through vibration measurement. The experiments consist in the measurement of the spindle FRF with an electromagnetic excitation device, at several speeds [START_REF] Tlalolini | Modeling and characterization of an electromagnetic system for the estimation of Frequency Response Function of spindle[END_REF]. A non-contact swept sine excitation of 80 N is applied over a range of 50-5000Hz to a dummy tool. The radial displacement at tool tip is measured by an Eddy current sensor. The radial force is estimated based on a measure of the output current and an identified model of the electromagnetic excitation device [START_REF] Tlalolini | Modeling and characterization of an electromagnetic system for the estimation of Frequency Response Function of spindle[END_REF]. Then, the FRF is computed by H1 formulation.

Besides, the spindle axial behavior is evaluated over a large range of load (-1800 N to +1500 N), during spindle rotations and with an axial loading device [START_REF] Rabréau | Phenomenological model of preloaded spindle behavior at high speed[END_REF]. The axial applied force F and the axial shaft deflection (in relation to spindle housing) u are measured and enable the study of the preload and of the axial stiffness evolution with spindle speed. All measures have been performed at the same steady-state operating temperature, at several spindle speeds between 4000 and 24 000 RPM (only results for these extreme values will be presented and discussed).

The studied spindle is a Fischer MFW2310, 70 kW, 24 000 RPM, that consists of 5 hybrid angular contact ball bearings (three SKF VEX70 at the front and two VEX60 at the rear) in back-to-back arrangements with two spring preload systems (Fig. 2). Here, it equips a Huron KX30 machine-tool.

Results

Figure 3 presents the results of the axial measurements (axial shaft deflection u and axial force F) at low and high spindle speeds. The different plots refer to the different spindle conditions introduced in Figure 1. At low speed, the slope at low force level shows that there are low variations of the axial stiffness. Note that, contrary to the other ones, Rep0 experiments were conducted over +/-2100 N. More deviations can be observed under -1500 N, where there is a significant drop of the stiffness. It is due to a contact loss in the front bearings. The behavior in this zone is determined by the preload systems and the rear bearings. At high speed, there is no contact loss in the bearings. The axial stiffness (the slope) only differs from one reassembly of the spindle to another one (between colors). The reference of shaft deflection u=0µm corresponds to F=0N at low speed, for every experimental campaign (see Fig. 3). At high speed, the shaft deflects out of the spindle housing (negative vertical displacement u), mainly due to centrifugal forces on balls. At F=0N, small variations can be observed between repaired spindles. Lower and more dispersed values are obtained for damaged spindles. At low speed, eigenfrequencies are visible at around 500 and 1100Hz. The frequency of the mode at 500Hz is constant with speed, but the other frequencies are evolving due do dynamic effects. Some modes appear at high speed (at 600 and 2100Hz).

It has been verified that the observed effects are greater than the variability of the experiment. According to the spindle condition, it can be noted that relatively similar FRFs are obtained. It is quite logical since the major design characteristics of the spindle, such as the shaft dimensions and the bearing types and positions, remain identical over time. However, interesting and nonnegligible variations of the dynamic behavior can be observed.

Larger deviations are obtained between reassemblies of the spindle, rather than between repaired and damaged conditions. Notably for Rep2 and Rep3, at low speed, the dominant mode (which is at 1150Hz for Rep0 and Rep1) has switched with the 1000Hz mode due to evolutions of mode coupling. This phenomenon is also visible at high speed, for Rep0. The frequency of the 500Hz mode is almost constant, but its amplitude decreases for Dam1 damaged spindles. The frequency can also vary, as for Dam2 at 900Hz at high speed.

Therefore, spindle dynamics can be affected by the spindle condition and, paradoxically, by the reassemblies.

Failure modes of HSM spindle

In order to explain the observed evolutions of FRF, the main failure of HSM spindle need to be identified and modeled. In a FMEA conducted in UsinAE project, around 250 reports of repairer's diagnostic at disassemblies of this spindle (used in aeronautic industry) had been analyzed. It revealed that the main failure modes of HSM spindles are, in order of importance: faulty ball bearings, modified geometry of HSK or shaft, defective bearing lubrication, clamping system or cooling system [START_REF] Maleki | Tailored ontology supporting sensor implementation for the maintenance of industrial machines[END_REF]. Note that, during any spindle repair, all bearings are changed (because they are damaged before or during disassembly) and preload is adjusted, to meet the manufacturer's specifications.

Concerning the failures of Dam0, Dam1 and Dam2, several ball bearings were systematically damaged (at front or rear bearings) and presented a white strip on the contact zone of raceways, which is common for aluminum HSM spindles in aeronautics [START_REF] De Castelbajac | Monitoring of distributed defects on HSM spindle bearings[END_REF]. Bearing lubrication was defective for Dam0. Components of the tool clamping system and of the rotary union were worn for Dam0 and Dam1. Dimensions were out of tolerance (of several microns) concerning the HSK interface (Dam0&1) and the shaft (Dam0). Thus, the observed damaged spindles are representative of the main failures of HSM spindle.

We have investigated the failure of a Dam1 rear bearing, due to a tool breakage (Fig. 5). The raceways were measured with Alicona Infinite Focus. The worn white strip is a distributed defect on the whole circumference of the contact zone on the raceways with an areal roughness Sa=936µm (versus 643nm for unworn zone). It consists of a number of micro-indentations, of an average depth of -5µm. This bearing failure mode is a falsebrinelling and it is due to the intensive vibrations [START_REF] Skf | Railway technical handbook[END_REF].

The main failure modes of HSM spindle are henceforth known and models can be proposed to predict their effects on dynamics.

Modeling of the spindle condition

Spindle model

The main components of the numerical simulations are presented in Figure 2. The shaft and the dummy-tool are modeled by 3D finite-elements. Mass and stiffness have been added on the location of the motor with 2D finite elements. The spindle housing is considered flexible and modeled by 3D finite-elements. The ball bushing between the rear sleeve and the spindle housing is modeled by a 5DoF spring joint: the preload springs in the axial direction and a 2x2DoF for rotations and translations in each radial direction. The machine-tool dynamics is modeled by a 2x2DoF joint, of which modal parameters are identified for the two principal modes of the machine-tool in each radial direction. The HSK interface is modelled by 2x2DoF stiffness. LMS Virtual.Lab software using linearized temporal computations were used for the multi-body dynamics simulations.

The bearings behavior is taken into account through 5x5 stiffness matrix that depends on the spindle speed. These matrices are obtained by the updating of the analytical model of the spindle axial behavior. It is based on a 5DoF model of the angular contact ball bearings [START_REF] Noel | Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings[END_REF]. Based on axial experiments (Fig. 3), a model update of the axial equilibrium of the preloaded system is performed. It provides the operating preload (force and stiffness) values and other bearing parameters with regards to the speed and the axial load [START_REF] Rabréau | Phenomenological model of preloaded spindle behavior at high speed[END_REF]. The preload systems and arrangement of bearings (with identified values) are integrated in the 3D model. The phenomenological model takes into account the dynamic effects in the bearings: centrifugal forces and gyroscopic moments on balls, and also macroscopic deformations of the shaft due to dynamic and thermal effects. The presence of a stroke limit and of friction between the rear sleeve and the housing is considered.

The transparent grey shaft in Figure 2 illustrates the spindle dominant mode, which is a shaft bending mode. Simulations of the 3D spindle model are in good agreement with the electromagnetic device FRF, particularly at high speeds.

Phenomenological modeling of spindle wear

The aim is to propose a model of the main failure modes of HSM spindle. It is used to investigate the impact of these failures on the dynamic behavior of HSM spindle. In order to understand the complete variability of spindle dynamics, not only failures are considered (evolution from repaired to damaged condition), but also the spindle repairs (variability between reassemblies). Variation ranges of parameters were chosen below based on experimental observations and literature.

Modified preload.

The axial experiments have shown variation of the preload, particularly between reassemblies (Fig. 3). Indeed, bearings are changed at each spindle repair. They are mounted with interference fit on the shaft and the preload is adjusted again, to meet the specification for the considered spindle. However, unwitting variation of preload force can be observed, due to tolerance stack that slightly modifies the springs compression. Preload springs are also sometimes changed (e.g. for Rep1). Consequently, this variability is modeled through modified values of preload force and preload spring stiffness. It modifies the axial equilibrium of the bearings arrangement and the new bearing stiffness matrixes are computed for the spindle 3D model. The impact of a preload variation has been studied. Figure 6 present the results for preloads of 1700N +/-200N . It is shown that variations of the preload force have a noticeable effect on all the spindle modes, modifying eigenfrequency and amplitude. Besides, lower spring stiffness decreases the main eigenfrequency and also modifies amplitudes. Deformed raceway of bearing. Due to wear or false-brinelling, the geometry of the bearing raceway can have been modified. Indeed, false-brinelling causes distributed defects on the entire raceway circumference, by vibrations-induced indentation. It modifies the ball-raceway contact. It can be modeled with a different value of raceway groove curvature radius in the bearing model. Two hypotheses have been formulated: the groove radius can have been increased around the contact zone by plastic deformations at a macroscopic scale; or micro-indentation has conducted to a smaller local groove radius, close to ball radius. Both hypotheses modify the bearing stiffness matrixes. Figure 7 presents the simulation results for +3% variation of the outer raceway groove radius, for front or rear bearings. It mainly modifies the amplitudes of the modes. Defective front bearings mainly affect the main mode, at 1000Hz; whereas defective rear bearings affect lower frequency mode, at 600Hz.

Modified bearing lubrication.

The distributed microindentation induced by false-brinelling degrades the roughness at the ball-raceway contact, as observed in Figure 5. It has an impact on the bearing lubrication, particularly the oil film thickness, which modifies the bearing damping [START_REF] Sarangi | Stiffness and damping characteristics of lubricated ball bearings considering the surface roughness effect. Part 2: numerical results and application[END_REF][START_REF] Dietl | Rolling Bearing Damping for Dynamic Analysis of Multi-body Systems-experimental and Theoretical Results[END_REF]. The consequence of an increased damping (from 12% to 22%) in the front or rear bearings has been studied. Figure 8 shows that defective front and rear bearings significantly modify the amplitude of the dominant mode, which is a shaft bending mode the but front and rear bearings also contribute. A defective rear bearing significantly decreases the amplitude of the 600Hz mode (rear sleeve mode). Defective tool clamping. Components of the tool clamping system can be worn or a lack of grease can modify the friction. As a consequence, the clamping force can be reduced, which decreases the stiffness at the HSK interface. Wang et al. [START_REF] Wang | Stiffness of HSK Tooling System in High Speed Machining[END_REF] have observed experimentally a stiffness drop of 67% for a clamping of 5kN (instead of 18kN). However, these variations are visible only if the angular deflections are large enough. Introducing Wang's values in the spindle model, Figure 7 shows that a 67% drop of the 2x2DoF HSK stiffness has a limited impact on the principal modes, whereas an extremely low clamping force was used.

Unbalanced rotor

After shocks, the balance of the rotor components assembly can be modified. It can be modeled by an additional mass on the shaft. The spindle temporal simulations have shown that small unbalance have no effect on the FRF. Only an increased contribution at the spindle frequency can be measured in the vibration spectrum.

Discussion

The experimental results can henceforth be analyzed and understood through the dynamics simulations with models of the main failure modes of HSM spindle. Paradoxically, the most significant changes of spindle dynamics that have been observed are not due to a worsened spindle condition, as expected, but they are due to preload variations at spindle reassembly. Indeed, a +100N increase (from 1700 to 1800N) of the preload force has been identified with the axial experiments for Rep2. It modifies the bearing stiffness matrices and the couplings around the dominant bending mode at 1000Hz. However, the spindle condition has also a noticeable impact on the spindle dynamics. The variations are mainly explained by bearing failures and the deterioration of their raceways. Compared to Rep1, Dam1 spindle dynamics presented lower amplitude for the rear sleeve mode at 500Hz (Fig. 4). Thus, the distributed false-brinelling that is common and was observed on bearing seems to increase the damping, as simulated in Figure 8.

As a synthesis, the major effect is due to the preload variation between spindle reassemblies that modifies the frequency and amplitude of the dominant mode, through bearing stiffness. Faulty bearings mainly affect the modes amplitude, through increased damping. The rear bearing has a great impact on the rear sleeve mode (at 600Hz); both front and rear bearings have a lower effect on the dominant bending mode amplitudes (at 1000Hz); and front bearing on the higher frequency bending modes. The frequency shift between repaired and damage conditions is mainly explained by the degraded HSK that modifies the tool-shaft bending.

Impact on cutting dynamics

Experiments were performed in order to investigate the impact of the spindle condition (Dam1 VS Rep2) on cutting dynamics. A shrink-fitted 20mm diameter, 176mm long and 2 teeth carbide tool was chosen for milling 7050 aluminum alloy. Figure 9 presents the stability lobes diagrams (SLD) that were obtained. Note that around 15 and 22kRPM, instability could not be reached due to excessively high vibration levels (asterisks). Elsewhere, a significant effect of spindle condition on cutting dynamics is revealed. At 10 000 RPM, the stability limit increases of 33% for Dam1 damaged spindle condition. Conversely, the stability limit drops of 26% at 18 000 RPM. Consequently, spindle condition can significantly degrade the quality of industrial production.

Conclusion

An original approach has been proposed for the study of the impact of spindle condition on the spindle dynamics. 7 experimental campaigns were performed over 3 years, with accurate characterization of the dynamic behavior with an electromagnetic device. The main failures modes of HSM spindle were identified (bearing false-brinelling notably) and modeled. Their effects on FRF were simulated and general conclusions were obtained concerning their impacts on spindle dynamics. Paradoxically, the major effect is due to preload variation between spindle reassemblies (for which additional attention should be paid). Faulty bearings mainly modify mode amplitudes and a degraded HSK modifies the frequencies of bending modes. The significant impact of spindle condition on cutting stability was also observed.

Figure 1 .

 1 Figure 1. Chronology of the experimental campaigns.

Figure 2 .

 2 Figure 2. Numerical model of the spindle.

Figure 3 .

 3 Figure 3. Spindle axial behavior, in relation to the spindle conditions.The results of radial experiment at low and high speeds with the electromagnetic excitation device are presented in Figure4. It reveals the dynamic behavior at the tip of the dummy tool.

Figure 4 .

 4 Figure 4. Spindle FRF, in relation to the spindle conditions.

Figure 5 .

 5 Figure 5. Bearing failure of Dam1 spindle condition.

Figure 6 .

 6 Figure 6. FRF simulations with modified preload.

Figure 7 .

 7 Figure 7. FRF simulations for faulty bearings or HSK interface.

Figure 8 .

 8 Figure 8. FRF simulations with increased bearing damping.

Figure 9 .

 9 Figure 9. Impact of spindle condition on stability lobes diagrams.

Journal homepage: www.elsevier.com/locate/cirp

Acknowledgment

This research was supported by the French FUI QuaUsi project. The authors thank Europe Technologies and Precise companies for their contribution.