Alexandr Klimchik
email: alexandr.klimchik@mines-nantes.fr

Dmitry Bondarenko
email: dmitry.bondarenko@irccyn.ec-nantes.fr

Anatol Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Sébastien Briot

Benoît Furet
email: benoit.furet@irccyn.ec-nantes.fr

Sébastien Briot

Benoît Furet

Sebastien Briot
email: sebastien.briot@irccyn.ec-nantes.fr

Benôit Furet

Compensation of Tool Deflection in Robotic-Based Milling

Keywords: Industrial robot, milling, compliance error compensation, dynamic machining force model, non-linear stiffness model

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Currently, robots become more and more popular for a variety of technological processes, including high-speed precision machining. For this process, the robot is subjected by external loading which caused by the machining force. This force is generated by the interaction between the tool mounted on the robot end-effector and the workpiece during the material removal [START_REF] Dépincé | Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF]. It is a contact force and it is distributed along the affected area of the tool cutting part. To evaluate the influence and to analyze the robot behavior while machining, the cutting force should be defined either experimentally or using accurate mathematical model.

To evaluate the force caused by interaction between the tool and the workpiece, two approaches can be used. The static approach allows computing the average cutting force without any consideration of dynamic aspect in machining system. This force serves as an external loading of the robot. This approach is widely used in analysis of conventional machining processes using CNC machines [START_REF] Altintas | Manufacturing automation, metal cutting mechanics, machine tool vibrations and CNC design[END_REF], where the stiffness is high. In contrast, robots have relatively low structural stiffness. For this reason, in the case of robotic-based machining, an additional source of dynamic displacements of the endeffector with respect to the desired trajectory induced by robot compliance may arise. Such behavior leads to the variable contact between the machining tool and the workpiece. Thus, the generated contact force depends on the current position of the robot end-effector on the trajectory. Consequently, the cutting force cannot be evaluated correctly using the static approach. In this case, the dynamic approach, which will be used in the paper, is required. It is based on computing of the force at each instant of machining process that defines loading of the robot for the next instant of processing. As a result, the dynamic aspect of robot motion under such variable cutting force can be examined for whole process.

Usually, in the robot-based machining this force causes essential deflections that decrease the quality of the final product. The problem of the robot error compensation can be solved in two ways that differ in degree of modification of the robot control software:

(a) by modification of the manipulator model, which better suits to the real manipulator and is used by the robot controller (in simple case, it can be limited by tuning of the nominal manipulator model, but may also involve essential model enhancement by introducing additional parameters, if it is allowed by a robot manufacturer);

(b) by modification of the robot control program that defines the prescribed trajectory in Cartesian space (here, using relevant error model, the input trajectory is generated in such way that under the loading the output trajectory coincides with the desired one, while input trajectory differs from the target one).

Moreover, with regard to the robot-based machining, there is a solution that does not require force/torque measurements or computations [START_REF] Dépincé | Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF], where the target trajectory for the robot controller is modified by applying the "mirror" technique. An evident advantage of this technique is its applicability to the compensation of all types of the robot errors, including geometrical and compliance ones. However, this approach requires carrying out additional preliminary experiments which are quite expensive. So, it is suitable for the large-scale production only. Another compensation methodology has been proposed by [START_REF] Eastwood | A gravitational deflection compensation strategy for HPKMs[END_REF] that was used for gravitational deflection compensation for hybrid parallel kinematic machines.

This paper focuses on the modification of control program that is considered to be more realistic in practice. This approach requires also accurate stiffness model of the manipulator. From point of view of stiffness analysis, the external and forces directly influence on the manipulator equilibrium configuration and, accordingly, may modify the stiffness properties. So, they must be undoubtedly taken into account while developing the stiffness model. However, in most of the related works the Cartesian stiffness matrix has been computed for the nominal configuration [START_REF] Kao | Conservative Congruence Transformation for Joint and Cartesian Stiffness Matrices of Robotic Hands and Fingers[END_REF][START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]. Such approach is suitable for the case of small deflections only. For the opposite case, the most important results have been obtained in [START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF][START_REF] Tyapin | Kinematic and elastostatic design optimization of the 3-DOF Gantry-Tau parallel kinamatic manipulator[END_REF]Pashkevich, 2011), which deal with the stiffness analysis of manipulators under the endpoint loading.

Thus, to compensate errors caused by the machining process, it is required to have an accurate stiffness model and precise cutting force model. In contrast to the previous works, the compliance error compensation technique presented in this work is based on the non-linear stiffness model of the manipulator (Pashkevich, 2011) and dynamic model of technological process that generates the cutting force.

PROBLEM STATEMENT

For the compliance errors, the compensation technique must rely on two components. The first of them describes distribution of the stiffness properties throughout the workspace and is defined by the stiffness matrix as a function of the joint coordinates. The second component describes the forces/torques acting on the end-effector while the manipulator is performing its machining task (manipulator loading).

The stiffness matrix required for the compliance errors compensation highly depends on the robot configuration and essentially varies throughout the workspace. From general point of view, full-scale compensation of the compliance errors requires essential revision of the manipulator model embedded in the robot controller. In fact, instead of conventional geometrical model that provides inverse/direct coordinate transformations from the joint to Cartesian spaces and vice versa, here it is necessary to employ the so-called kinetostatic model [START_REF] Su | A Polynomial Homotopy Formulation of the Inverse Static Analyses of Planar Compliant Mechanisms[END_REF]. It is essentially more complicated than the geometrical model and requires rather intensive computations that are presented in Section 3.. The dynamic behavior of the robot under the loading F caused by technological process can be described as

C C C    M δt C δt K δt F (1) where C
M is 6 6  mass matrix that represents the global behavior of the robot in terms of natural frequencies, C

C is 6 6  damping matrix, C K is 6 6  Cartesian stiffness matrix of the robot under the external loading F , , δt δt and δt are dynamic dis- placement, velocity and acceleration of the tool endpoint in a current moment respectively [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF].

In general, the cutting force F c has a nonlinear nature and depends on many factors such as cutting conditions, properties of workpiece material and tool cutting part, etc [START_REF] Ritou | A new versatile in-process monitoring system for milling[END_REF]. But, for given tool/workpiece combination, the force F c could be approximated as a function of an uncut chip thickness h, which represents the desired thickness to cut at each instant of machining.

Hence, to reduce errors caused by cutting forces in the robotic-based machining it is required to obtain an accurate elasto-static model of robot and elasto-dynamic model of machining process. These problems are addressed in the following sections taking into account some particularities of the considered application (robotic-based milling).

MANIPULATOR MODEL

Elasto-static model

Elasto-static model of a serial robot is usually defined by its Cartesian stiffness matrix, which should be computed in the neighborhood of loaded configuration. Let us propose numerical technique for computing static equilibrium configuration for a general type of serial manipulator. Such manipulator may be approximated as a set of rigid links and virtual joints, which take into account elasto-static properties (Figure 1). Since the link weight of serial robots is not negligible, it is reasonable to decompose it into two parts (based on the link mass centre) and apply them to the both ends of the link. All this loadings will be aggregated in a vector

  1 ... n  G G G , where i
G is the loading applied to the i-th node- point. Besides, it is assumed that the external loading F (caused by the interaction of the tool and the workpiece) is applied to the robot end-effector.

...

  T T T θ 1 δ δ δ n jj j        G t F t τθ (2)
where the virtual displacements δ j t can be com- puted from the linearized geometrical model derived from

() θ δ δ , 1.. j j jn  tJθ , which includes the Jaco- bian matrices   () θ , j j    J
g q θθ with respect to the virtual joint coordinates.

So, expression (2) can be rewritten as

    T ( ) T ( ) T θ θ θ 1 δ δ δ n jn j j          GJ θ F J θ τ θ
(3) which has to be satisfied for any variation of δθ . It means that the terms regrouping the variables δθ have the coefficients equal to zero. Hence the force balance equations can be written as

( )T ( )T θ θ θ 1 n jn j j       τ J G J F (4)
These equations can be re-written in blockmatrix form as

(G)T (F)T θ θ θ     τ J G J F (5)
where

(F) ( ) θθ n  JJ , (G) (1) ( ) θ θ θ T TT ... n    J J J , T TT 1 ... n    G G G
. Finally, taking into account the virtual spring reaction θθ  τ K θ , where

  1n θ θ θ ,..., diag  K K K
, the desired static equilib- rium equations can be presented as

(G)T (F)T θ θ θ      J G J F K θ (6)
To obtain a relation between the external loading F and internal coordinates of the kinematic chain θ corresponding to the static equilibrium, equations ( 6) should be solved either for different given values of F or for different given values of t . Let us solve the static equilibrium equations with respect to the manipulator configuration θ and the external load- ing F for given end-effector position

   tg θ and the function of auxiliary-loadings   G θ     (G)T (F)T θ θ θ ;;      K θ J G J F t g θ G G θ (7)
where the unknown variables are   , θF . Since usually this system has no analytical solution, iterative numerical technique can be applied. So, the kinematic equations may be linearized in the neighborhood of the current configuration i θ

      (F) θ 11 ; i i i i i      tg θ J θ θ θ (8)
where the subscript 'i' indicates the iteration number and the changes in Jacobians (G) (F) θθ , JJ and the auxil- iary loadings G are assumed to be negligible from iteration to iteration. Correspondingly, the static equilibrium equations in the neighborhood of i θ may be rewritten as

(G)T (F)T θ 1 θ 1 θ ii       J G J F K θ . (9)
Thus, combining (8), ( 9) and analytical expression for

1 (G)T (F)T θ θ θ ()      θ K J G J
F , the unknown variables F and θ can be computed using follow- ing iterative scheme

        1 (F) 1 (F)T θθ (F) (F) 1 (G)T θ θ θ 1 (G)T (F)T θθ 1 1 1 θ 1 • i i i i i i i i                     F J K J tg θ J θ J K J G θ K J G J F (10)
The proposed algorithm allows us to compute the static equilibrium configuration for the serial robot under external loadings applied to any point of the manipulator and the loading from the technological process.

Stiffness matrix

In order to obtain the Cartesian stiffness matrix, let us linearize the force-deflection relation in the neighborhood of the equilibrium. Following this approach, two equilibriums that correspond to the manipulator state variables ( , , ) F θt and ( δ , δ , δ )    FF θ θ t t should be considered simul- taneously. Here, notations δF , δt define small in- crements of the external loading and relevant displacement of the end-point. Finally, the static equilibrium equations may be written as

  (G)T (F)T θ θ θ ;       tg θ K θ J G J F (11) and             T (G) (G) θ θ θ T (F) (F) θθ δδ δ δ δ δδ              t t g θθ K θ θ J J G G J J F F (12)
where θ , , , , t F G K θ are assumed to be known. After linearization of the function () g θ in the neighborhood of the loaded equilibrium, the system (11), ( 12) is reduced to equations

(F) θ (G) (G) (F) (F) θ θ θ θ θ δδ δ δ δ δ δ       tJ θ K θ J G J G J F J F (13)
which defines the desired linear relations between δt and δF . In this system, small variations of Jaco- bians may be expressed via the second order derivatives

(F) (F) θ θθ δδ  JHθ , (G) (G) θ θθ δδ  JHθ , where (G) 2 θθ 1 2 T j j j n      H g G θ , (F θθ 2 )2 T    H g F θ .
Also, the auxiliary loading G may be computed via the first order derivatives as δδ     GG θθ Further, let us introduce additional notation

(F) (G) (G)T θθ θθ θθ θ      H H H J G θ
, which allows us to present system (13) in the form

(F) θ (F)T θ θ θθ δδ δ                  0J tF 0 θ J K H (14)
So, the desired Cartesian stiffness matrices C K can be computed as

  1 (F) 1 (F)T C θ θ θθ θ ()    K J K H J (15)
Below, this expression will be used for computing of the elasto-static deflections of the robotic manipulator.

Mass matrix

To evaluate dynamic behaviour of the robot under the loading, in addition to the Cartesian stiffness matrix C K it is required to define the mass matrix C M . Comprehensive analysis and definition of this matrix have been proposed in [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF]. Here, let us summarise the main results that will be used further in the error compensation technique.

Similar to the stiffness matrix, here physical properties defined by the mass matrix

z i ii a I L   , i m is physical mass of i-th link, i  is density of i-th link, i L is link length, p i I is the polar moment, , yz ii
II are the second moments of the area. Since the mass matrix θ M is defined in the joint coordinates it can be recomputed into Cartesian coordinates associated with the tool end-point using the Jacobian matrix θ J (which depend on the robot configuration q and computed with respect to virtual joint coordi- nates θ ) using following expression θ C θθ

T  M J M J (17)
Thus, using expressions ( 16) and ( 17) it is possible to compute mass matrix C M for a given ro- bot configuration q .

MACHINING PROCESS

Let us obtain the model of the cutting force which depends on the relative position of the tool with respect to the workpiece at each instant of machining. As follows from previous works [START_REF] Brissaud | The Fractional Model for the Determination of the Cutting Forces[END_REF], for the known chip thickness h, the cutting force F c can be expresses as

    2 0 ,0 1 ss c s p h h r h h F h k a h hh    ( 18
)
where p a is a depth of cut, 0 1 r k k   depends on the parameters k ∞ , k 0 that define the so called stiffness of the cutting process for large and small chip thickness h respectively (Figure 2) and h s is a specific chip thickness, which depends on the current state of the tool cutting edge. The parameters k 0 , h s , r are evaluated experimentally for a given combination of tool/working material. To take into account the possible loss of contact between the tool and the workpiece, expression (18) should be supplement by the case of 0 h  as h s ~k∞ ~k0

F c h Range of h in case of conventional CNC machining

Figure 2

Fractional cutting force model

F c (h)   0, if 0 c F h h  ( 19
)
For the multi-edge tool the machining surface is formed by means of several edges simultaneously. The number of working edges varies during machining and depends on the width of cut. For this reason, the total force F c of such interaction is a superposition of forces F c,i generated by each tool edge i, which are currently in the contact with the workpiece. Besides, the contact force F c,i can be decomposed by its radial F r,i and tangential F t,i components (Figure 3). In accordance with Merchant's model [START_REF] Merchant | Mechanics of metal cutting process. I-Orthogonal cutting and type 2 chip[END_REF], the t-component of cutting force F t,i can be computed with the equation ( 18). The r-component F r,i is related with F t,i by following expression [START_REF] Laporte | Phenomenological model for drilling operation[END_REF] ,,

r i r t i F k F  ( 20
)
where the ratio factor k r depends on the given tool/workpiece characteristics.

It should be mentioned that in robotic machining it is more suitable to operate with forces expressed in the robot tool frame {x,y,z}. Then, the corresponding components F x , F y (Figure 3) of the cutting force F c can be expressed as follows ,, 11 ,, 11 cos sin sin cos

zz zz nn x r i i t i i ii nn y r i i t i i ii F F F F F F           (21)
where n z is the number of currently working cutting edges, φ i is the angular position of the i-th cutting edge (the cutting force in z direction F z is negligible here). So, the vector of external loading of the robot to the machining process can be composed in the frame {x,y,z} using the defined components F x , F y as F=[F x ,F y ,0,0,0,0] T . It should be stressed that the cutting force components F r,i , F t,i mentioned in equation ( 18),(20) are computed for the given chip thickness h i , which should be also evaluated. Let us define model for h i using mechanical approach. Then the chip thickness h i removed by i-th tooth depends on the angular position φ i of this tooth and it can be evaluated using to the geometrical distance between the position of the given tooth i and the current machining profile (Figure 3). It should be mentioned, that the main issue here is to follow the current relative position between the i-th tooth and the working material or to define whether the i-th tooth is involved in cutting for given instant of process. Because of the robot dynamic behavior and the regenerative mechanism of surface formation [START_REF] Tlusty | Basic non-linearity in machining chatter[END_REF] this problem cannot be solved directly using kinematic relations.

 i i F r,i F t,i F x F y v f W TCP Workpiece Cutter y x Cutting direction
In this case it is reasonable to introduce a special rectangular grid, which decomposes the workpiece area into segments and allows tracking the tool/workpiece interaction and the formation of the machining profile (Figure 4).

Here, Steps Δs x , Δs y between grid nodes are constant and depend on the tool geometry, cutting condition and time discretization Δτ. Each node j ( 1, w jN  , N w is the number of nodes) of the grid can be marked as "1" or "0": "1" corresponds to nodes situated in the workpiece area with material (rose nodes in Figure 5), "0" corresponds to nodes situated in workpiece area that was cut away (white nodes in Figure 5).

In order to define the number of currently cut nodes by the i-th tooth, the previous instant of machining process should be considered. Let us define A i as an amount of working material that is currently cut away by the i-th tooth (Figure 5). So, if node j marked as "1" is located inside the marked sector (green nodes in Figure 5), it changes to "0" and A i is increasing by xy ss DD . Analyzing all potential nodes and computing A i , the chip thickness h i , removed at given instant of the process by the i-th tooth, can be estimated by ,

i i i h A R  D 1, z iN 
. The angle Δφ i determines the current angular position of the ith tooth regarding to its position at the instant τ-Δτ and referred to the position of TCP at τ-Δτ. Evaluating the tool/workpiece intersection A i and computing the corresponding chip thickness h i Described mechanism of chip formation and the machining force model ( 18) allow computing the dynamic behavior of the robotic machining process where models of robot inertia and stiffness are discussed in the section 3 of the paper. The detailed algorithm that is used in numerical analysis is presented in Figure 6, where the analysis of the robot dynamics is performed in the tool frame with respect to the dynamic displacement of the tool δt dyn fixed on the robot end-effector around its position on the trajectory.. Computing the radial and tangential cutting force components Frj, Ftj

Computing the machining forces Fx, Fy (tool frame)

Analysis of dynamics of the robotic machining process (tool frame)

Next time step τ :

τ=τ0:Dτ:τMAX Z j , 1  1 1 1 dyn c c dyn c c dyn c        δt M C δt M K δt M F Figure 6
Algorithm for numerical simulation of robotic machining process dynamics

COMPLIANCE ERROR COMPENSATION TECHNIQUE

In industrial robotic controllers, the manipulator motions are usually generated using the inverse kinematic model that allows us to compute the input signals for actuators 0 ρ corresponding to the desired end-effector location 0 t , which is assigned assum- ing that the compliance errors are negligible. However, if the external loading F is essential, the ki- nematic control becomes non-applicable because of changes in the end-effector location. It can be computed from the non-linear compliance model as

  1 F 0 | f   t F t ( 22
)
where the subscripts 'F' and '0' refer to the loaded and unloaded modes respectively, and ' | ' separates arguments and parameters of the function   f . Some details concerning this function are given in our previous publication (Pashkevich, 2011).

To compensate this undeterred end-effector displacement from 0 t to F t , the target point should be modified in such a way that, under the loading F , the end-platform is located in the desired point 0 t . This requirement can be expressed using the stiffness model in the following way

  (F) 00 | f  F t t (23)
where (F) 0 t denotes the modified target location. Hence, the problem is reduced to the solution of the nonlinear equation ( 23) for (F) 0 t , while F and 0 t are assumed to be given. It is worth mentioning that this equation completely differs from the equation

0 ( | ) f  F t t ,
where the unknown variable is t . It means that here the compliance model does not allow us to compute the modified target point (F) 0 t straightforwardly, while the linear compensation technique directly operates with Cartesian compliance matrix [START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF].

To solve equation ( 23) for (F) 0 t , similar numerical technique can be applied. It yields the following iterative scheme

  (F) (F) 1 (F) 0 0 0 0 ( | ) • f       t t t F t (24)
where the prime corresponds to the next iteration, (

is the scalar parameter ensuring the convergence. More detailed presentation of the developed iterative routines is given in Figure 7.

  F t ) F ( 1 0 | ε f    F t t Compute end-platform location No Yes   1 (F) 0 | f   F t F t Next iteration     0 ; ( ) j j    t t F 0 F F t Obtain trajectory     j j  t t F t Repeat for different t0   0 j  t t Stopping criterion is satisfied?   (F) (F) 0 0 0       F t t t t 0 (F) 0  t t
Compute corresponding equilibrium configuration and stiffness matrix

( ) C
, , ,

i i i i q F K θ Figure 7
Procedure for compensation of compliance errors Hence, using the proposed computational techniques, it is possible to compensate a main part compliance errors by proper adjusting the reference trajectory that is used as an input for robotic controller. In this case, the control is based on the inverse kinetostatic model (instead of kinematic one) that takes into account both the manipulator geometry and elastic properties of its links and joints. Implementation of developed compliance error compensation technique presented in Figure 8.

EXPERIMENTAL VERIFICATION

The developed compliance error compensation technique has been verified experimentally for robotic milling with the KUKA KR270 robot along a simple trajectory in aluminum workpiece. It is assumed that at the beginning of the technological process the robot is in the configuration q (see Ta- ble 1 Figure 9). The parameters of the stiffness model for the considered robot have been identified in [START_REF] Dumas | Joint stiffness identification of six-revolute industrial serial robots[END_REF] and are presented in Table 1. Link masses required for the mass matrix of the robot are presented also in Table 1.

Table 1

Initial data for robotic-based milling Joint coordinates, [deg] q 1 q 2 q 3 q 4 q 5 q 6 90 -50 120 Taking into account that the workpiece has a straight borders let us assume that at the instant t=0 one of the teeth of the tool is in contact with the workpiece material as it is shown in the Figure 10. It is also assumed that the machining process is performing with the constant feed rate v f =4 m/min (applied in x-direction of the robot tool frame) and the constant spindle rotation Ω=8000 rpm along the straight line of 80 mm. Experimental verification and numerical simulation of the described case of the milling process with KUKA KR-270 robot using the algorithm shown in Figure 6 allows us to trace the evolution of machining force x,y-components for the whole process (Figure 11). The corresponding dynamic displacement of the tool around its current position on the trajectory is shown in Figure 12 In accordance with the obtained results the system robot/machining process realize complex vibratory motion. The high frequency component of this motion (about 700 Hz, Figure 11) is related to the spindle rotation and the number of tool teeth N z . In certain cases such behavior can excites the dynamics of the robot (natural modes) but this study remains out the frame of the presented paper. On the contrary, the low frequency component of robot/tool motion (about 7 Hz, Figure 12), especially in the y-direction (that is perpendicular to the applied feed) influences directly the quality of final product. Such motion is related to the robot compliance and it can be compensated using the error compensation technique described in the paper. Hence, let us form the modified trajectory based on the dynamic displacement of the robot end-effector in the y-direction (Figure 13):

It should be stressed that the time step between referenced points of this modified trajectory is limited with the characteristics of the controller used in the robot (in the presented case this step is chosen 0.05 sec). The corresponding feed rate v fy for the modified trajectory has been computed. So, this new data (feed f y and feed rate v fy ) with the data defined in the beginning of this section allow us to compensate the trajectory error during machining caused by the robot compliance. The resulted compensated trajectory in the y-direction (in time domain) is presented in Figure 14. It should be noted that the part of the trajectory while machining tool is engaging into the workpiece does not have effect on the quality of final product (surface). During this stage the contact area between the tool and the workpiece is increasing progressively. Hence, at each instant of processing the cutter corrects the machining profile and eliminates trajectory errors produced during all previous instants. On the contrary, during the stage of machining with the fully engaged tool the trajectory in x,y-directions define directly the final machining profile and this part of trajectory is analyzed here (Figure 14). Comparison results presented in Figure 12 and Figure 14 are summarized in Table 2. So after applying error compensation technique the static deviation in y direction has been reduced from 0.058 mm to 0.00014 mm (99.8%). Maximum defilation in the machining profile has been reduced from 0.063 mm to 0.0047 mm (92.6%). Low frequency remained the same for both cases.

Table 2

Milling trajectory accuracy before and after compliance error compensation Hence, obtained results show that the developed compliance error compensation allows us significantly increase the accuracy of the robotic-based machining.

CONCLUSION

In robotic-based machining, an interaction between the workpiece and technological tool causes essential deflections that significantly decrease the manufacturing accuracy. Relevant compliance errors highly depend on the manipulator configuration and essentially differ throughout the workspace. Their influence is especially important for heavy serial robots. To overcome this difficulty this paper presents a new technique for compensation of the compliance errors caused by technological process. In contrast to previous works, this technique is based on the nonlinear stiffness model and the reduced elastodynamic model of the robotic based milling process.

The advantages and practical significance of the proposed approach are illustrated by milling with of KUKA KR270. It is shown that after error compensation technique significantly increase the accuracy of milling. In future the proposed technique will be integrated in a software toolbox.

  Figure 3Forces of tool/workpiece interaction

  Figure 5Evaluating the tool/workpiece intersection A i and computing the corresponding chip thickness h i

  Computing the of TCP and position/ orientation of all teeth (tool frame) Analysis of interaction between j-th tooth and grid nodes, computing Aj Updating the grid, computing the chip thickness hj Definition of the external loading caused by technological process (tool frame) : F=[Fx Fy 0

  Figure 8Implementation of compliance error compensation technique

Figure 9 Figure 11

 911 Figure 9 Starting pose of the KUKA KR270 robot to perform the operation of milling For the milling, the cutter with the external diameter D=20 mm and four teeth (N z =4) distributed uniformly over the tool is used. For the given combination of the tool and the workpiece material the following parameters correspond to the cutting force model defined in (18): k 0 = 6 5 10  N/m, h s = 5 1.8 10   m, r=0.1, k r =0.3.

Figure 14

 14 Figure 14Evolution of the dynamic displacement obtained after involving the error compensation technique into the analysis of robotic milling process

  . Figure 12 Evolution of the tool dynamic displacement δt dyn that is composed from x TCP and y TCP components Figure 13 Modified trajectory f y and corresponding feed rate v fy in y-direction, computed based on the original dynamic displacement of the tool δt dyn

		0.08
		0.06
	δt dyn , [mm]	0 0.02 0.04				x TCP y TCP
	-0.02
	0 -0.04	0.2	0.4	0.6	0.8	1	1.2 τ, [s]
		-5
		x 10 10 -2
	f y , y TCP , [mm]	-5 0 5		y TCP Referenced points Modified trajectory in y-direction
		0	0.2	0.4	0.6	0.8	1	1.2 τ, [s]
		-5
		x 10
		0
	v fy , [m/s]	-2		Referenced points
		-4		Modified feed rate in y-direction
		0	0.2	0.4	0.6	0.8	1	1.2 τ, [s]

ACKNOWLEDGEMENTS

The authors would like to acknowledge the financial support of the ANR, France (Project ANR-2010-SEGI-003-02-COROUSSO) and the Region "Pays de la Loire", France.


Dmitry Bondarenko

Anatol Pashkevich

Sébastien Briot

Mathieu Ritou

Benoît Furet

Sébastien

) A Pashkevich

) S Briot

Elasto-Dynamic Model of Robotic Milling Process Considering Interaction between Tool and Workpiece

Keywords: Robot elasto-dynamic model, cutting force model, vibratory analysis, milling

published or not. The documents may come

L'archive ouverte

INTRODUCTION

Parallel robots have increasingly been used in industry in the last few years, mainly for pick-and-place applications or high-speed machining [START_REF] Brogardh | Present and future robot control developmentan industrial perspective[END_REF], [START_REF] Chanal | A study of the impact of machine tool structure on machining processes[END_REF]. This interest is due to their main properties, i.e. their higher rigidity and dynamic capacities compared with their serial manipulator counterparts.

Having a good knowledge of the elasto-dynamic behavior of a manipulator plus its interactions is a crucial point. In this sense, accurate elasto-dynamic models are necessary at both the control stage [START_REF] Singer | Residual vibration reduction in computer controlled machines[END_REF]- [START_REF] Pelaez | Input shaping reference commands for trajectory following Cartesian machines[END_REF] and design stage [START_REF] Voglewede | Overatching framework for measuring the closeness to singularities of parallel manipulators[END_REF]- [START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF], in order to optimize the geometry, as well as the shape of the elements of the manipulator. This will lead to the creation of a mechanism in which vibrations will be minimized.

Several elasto-dynamic models have been proposed and used in the literature. Three main general methods can be distinguished:  Finite element analysis (FEA); the FEA method is proved to be the most accurate and reliable, since the links/joints are modeled with its true dimension and shape [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot uncoupled kinematics[END_REF], [START_REF] Shabana | Dynamics of multibody systems[END_REF]- [START_REF] El-Khasawneh | Computation of stiffness and stiffness bounds for parallel link manipulators[END_REF], but is highly time-consuming. This method is usually applied at the final design stage for the verification and component dimensioning.  Matrix structural analysis (MSA) method is a common technique in mechanical engineering [START_REF] Deblaise | Systematic analytical method for PKM stiffness matrix calculation[END_REF], [START_REF] Imbert | Analyse des structures par éléments finis[END_REF]; it incorporates the main ideas of the FEA but operates with rather large flexible elements (beams, arcs, cables, etc.). This obviously yields reduction of the computational expense and, in some cases, allows analytical stiffness matrix to be obtained. However, this method can only be applied to links with simple shapes and requires improved skills in FEA.  Virtual joint methods (VJM) [START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF], [START_REF] Wittbrodt | Dynamics of flexible multibody systems[END_REF], which is also referred to as ''lumped modeling", is based on the expansion of the traditional rigid model by adding virtual joints (localized springs), which describe the elastic deformations of the manipulator components (links, joints and actuators).

Generally, lumped modeling is simpler to use than MSA and provides acceptable accuracy in reduced computational time. The main limitations of these models are:  Because of the large number of elements they use, they are still highly time-consuming and can lead to few accurate results because of the problem dimension;  They do not incorporate accurate modeling of the robot interactions, especially regarding the milling process. The purpose of this paper is to propose a reduced elastodynamic modeling approach on parallel robots, based on VJM, combined with an accurate modeling of the milling process efforts. This approach is validated on a KUKA KR240 robot used for milling of an aluminum workpiece for automobile industry.

REDUCED ELASTO-DYNAMIC MODELING OF ROBOTS 2.1 Problem Statement

Let us consider a general robot made of n elements and m actuators. Considering that the robot is in an equilibrium position, its elastic potential energy is given by [START_REF] Shabana | Dynamics of multibody systems[END_REF] 12 T e tot V  q K q , where K tot is the global stiffness matrix of size p×p, p being the number of elastic coordinates taken into account, and q the vector of elastic generalized coordinates. The robot kinetic energy is given by 12 T tot T  q M q , where M tot is the global mass matrix of size p×p. Starting from these definitions, and considering that external forces F are applied, it can be demonstrated that the system is governed by the relation [START_REF] Shabana | Dynamics of multibody systems[END_REF], [START_REF] Wittbrodt | Dynamics of flexible multibody systems[END_REF]:

tot tot  M q K q F (1)
Solving this problem involves finding the p generalized of the problem and also inverting the p×p matrix M tot , which is highly time-consuming and can lead to few accurate results because of the problem dimension. In order to avoid such kind of drawbacks [START_REF] Pashkevich | Stiffness analysis of over constrained parallel manipulators[END_REF] has recently proposed a new procedure to compute the deformations of the robot when a force is applied at the end-effector. This procedure computes a stiffness (and also compliance) matrix K r (and S r , resp.) of dimension six that represents the behavior of the robot in terms of deformations. Moreover, during the procedure, only inversions of 6 dimensional matrices are involved, which considerably reduces computational time and avoids accuracy problems due to illconditioning of the large global stiffness matrix. Thus, the global p dimensional problem defined with respect to all variables q has been reduced to a 6 dimensional problem defined with respect to platform deflections t only. As a result, the entire robot can locally be seen as a virtual spring of dimension six that deforms when applying a wrench on the end-effector.

Starting from these considerations, it would be interesting to reduce the dimension of the problem by expressing equation (1) as a function of the reduced stiffness matrix K r , of the platform deflections t and, also, as a function of a reduced mass matrix that will be denoted as M r , which could represent the global behavior of the robot in terms of natural frequencies. The first section of the paper will be focused on this problem.

Computation of the Robot Natural Frequencies -Rayleigh-Ritz Approximation

Considering the robot with free oscillations, equation ( 1) becomes tot tot  M q K q 0 . A solution q l of this equation can be found by solving the system:

  2 l tot tot l   M K q 0 , with 2 ll f   (2)
where q l represents the vectors of the shape of free vibrations of the system for the l-th natural mode, f l and  l are their corresponding natural frequency and pulsation, respectively. If the matrix 2 l tot tot   MK is singular (which is always the case when  l is one of the modal pulsation of the robot), q l becomes non-null and is the eigenvector corresponding to the pulsation  l . There is an infinity of vectors q l validating (2) (for a given  l ), but all are proportional to the others. Vector q l is not necessarily dependent of time, but almost represents the amplitude of the vibrations. Therefore, when only the l-th mode of the system is excited, the displacements of all springs may be written under the form:

  * sin l l l l t   qq (3)
where  l is a phase difference corresponding to the mode l.

When all modes are excited, the displacement q of all spring centers may be written in the form:

  * 11 sin pp l l l l ll t       q q q (4)
The most common way to find the values of the pulsation  l is to solve the eigenvalue problem

  21 det 0 l tot tot    I M K ( 5
)
where I is an identity matrix of dimension p.

Another way to find the natural pulsation  l of the mode l would be to know the exact amplitude  l of the displacements q l for this mode. For the natural mode l, the potential and kinetic energies of the system are given by:

    * * 2 * * 2 2 1 2 1 2 sin 1 2 1 2 cos TT e l tot l l tot l l l TT l tot l l l tot l l l Vt Tt           q K q q K q q M q q M q (6)
From the principle of energy conservation, it follows that

    max max 0 e VT  , i.e.   2 0 T l l tot tot l   q M K q (7)
It is obvious that the exact knowledge of the amplitude q l is an impossible task without a direct measure on the system of all the displacements of the robot nodes. However, this vector may be approximated by another denoted as ˆl q that is close to the exact amplitude q l . Introducing this approximated vector ˆl q into (7) will allow us to find a corresponding value of ˆl  and, as a result, ˆl f which is the approximated natural frequency of the system. Such kind of elasto-dynamic problem resolution is called the Rayleigh-Ritz method [START_REF] Birkhoff | Rayleigh-Ritz approximation by piecewise cubic polynomials[END_REF].

The better the approximation, the more accurate the value of ˆl  . The designer's skills in terms of understanding and analysis of robot physical behavior here are of the utmost importance. In this sense, let us recall that the first natural frequencies are associated with the highest level of energy due to vibrations, and represent the highest displacements of the structure.

Using the Rayleigh-Ritz approximation in order to compute the first natural frequencies, the stresses for which the maximal displacements appear have to be found. From our experience in elastic behavior of robots, it is assumed that a good approximation of these maximal displacements will be the deformation of the robot with a load applied at the end-effector, and it can be shown in the following that this hypothesis is valid. Using this assumption, the displacements q l of all springs can be computed as a function of the end-effector displacements t, i.e. it is possible to define a matrix J q such as:

l  p qJ δt (8)
As a result, introducing (8) into (7) will lead to

  2 T T T l tot tot   q p q q δt J M J J K J δt 0 (9)
where the matrices T tot qp J M J and T tot qq J K J are now of dimension 6. It should be mentioned here that, in the case where an external load is applied at the end-effector only, the term T tot qq J K J is equal to the reduced stiffness matrix K r of the robot. Therefore, (9) can be rewritten as:

  2 T l   rr δt M K δt 0 , with T tot  r q q
M J M J [START_REF] Cammarata | The elastodynamics of a 3-CRU spherical robot[END_REF] In the following sections it is explained how to obtain expressions [START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF] and [START_REF] Cammarata | The elastodynamics of a 3-CRU spherical robot[END_REF].

Reduction of the Link Mass Matrix

It is possible to decompose the previously cited task into two sub-problems. First, the displacements of a beam j can be expressed as a function of the displacement of its extremities. Then, one can express the beam extremity displacements as a function of the platform displacements t. Using this approach will allow for a reduction in the size of the link mass matrices, and thus avoiding creating global mass matrices M tot with very large dimension.

Two main ways can be followed to reduce the size of the link mass matrices. The first one consists in discretizing the beam j into p j rigid links and springs and to express their displacements as a function of the beam extremity displacements. However, such numerical method must be repeated for each link and, thus, increases the size of the algorithm and decreases its efficiency. As a result, it is preferred to use the following procedure which allows analytical expressions to be obtained for the reduced link mass matrices. Let us consider the link j, modeled as a beam (Figure 1). At this beam is attached a local frame represented by the vectors x j , y j and z j . Before any deformation of the system, the beam j is linked (rigidly or by a passive joint) to beams (j-1) and (j+1) at points O j and O (j+1) , respectively (Figure 1). After deformation of the robot, the beam extremity located at O j is displaced from

1 1 1 1 2 6
, , , T j j j q q q    1 j q and the one located at

O (j+1) is displaced from 2 2 2 2 1 2 6
, , ,

T j j j q q q    j q
, where the three first components of each vector correspond to the translational displacements along local x j , y j and z j axes, respectively, and the three last components to the rotational displacements along the same axes. It should be mentioned that in the remainder of the paper, the left superscript "0" will stand for the coordinates expressed in the global frame. If no left superscript is mentioned, the vector is expressed in the local frame attached to the link j.

The general formula for the kinetic energy of an elastic Bernoulli beam is equal to [START_REF] Imbert | Analyse des structures par éléments finis[END_REF]:

  0 12 diag , , , , , j L T ij j p y z j j j j j j T dx A A A I I I     j j j j q Q q Q (11)
In this expression, j q represents the velocity of the beam cross-section located at position x from the local reference frame (Figure 1), L j is the length of the beam j,  j the mass density at cross-section x, A j its area, p j I its torsional constant and y j I , z j I , the quadratic momentums along y j and z j , respectively.

For the l-th natural mode, and from ( 4), the kinetic energy can be rewritten as:

  22 0 1 2 cos j L T j l l l j T t dx       j j j q Q q (12)
q j being the amplitude of the displacement of the beam crosssection located at position x from the local reference frame (Figure 1). In the Rayleigh-Ritz approximation, considering that the deformations due to the natural vibrations are similar to those obtained when an external load is applied at the robot endeffector only, each link of the structure will deform due to the stresses transmitted through the robot joints at points O j . As a result, the deformations  j of the beam cross-section (Figure 1) can be approximated by the deformations of a tip-loaded beam, given by [START_REF] Imbert | Analyse des structures par éléments finis[END_REF]   diag , , , , ,

j j j j j j f g g f h h  2 jj δδ ( 13
)
where

  2 j xL  jj δδ
represents the deformation of the beam at its tip and

          23 2 , 0.5 3 , 0.5 j j j j j j j j f x x L g x x L x L h x x L x L     ( 14
)
As a result, the global displacement q j of the beam crosssection at x can be expressed as a sum of two terms:

  3 33       1 j j j ID qq δ 0I , with   0 0 0 00 00 Dx x        (15)
In this sum, the left terms corresponds to the displacement of the undeformed beam due to the displacement of the node located at O j . Introducing ( 13) to ( 15) into ( 12) leads to the following equation:

    2 2 0 0 0 0 0 1 2 cos T red j l l l j Tt             1 2 1 2 j j j j q q M q q (16)
where the expressions of each components of matrix 0 red j M are given in [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF].

Reduction of the Robot Mass Matrix

Using the results of the previous section, the total kinetic energy of the system is given by:

  2 1 2 cos T j l l l tot j T T t        q M q (17) with   0 0 0 12 diag , , , red red red tot p  M M M M and             0 1 0 2 0 1 0 2 0 1 0 2 1 1 2 2 , , , , , ,
T T T T T T T pp q     q q q q q q
It is necessary to express the relationship linking all vectors 0 v j q (v = 1, 2) to the end-effector displacements t. For a robot composed of k legs and using the VJM [START_REF] Pashkevich | Stiffness analysis of over constrained parallel manipulators[END_REF], this displacement is equal to:

i i pi i   δt J θ J p ( 18
)
where  i represents the deformations of all virtual springs of the leg i and p i the displacements of its passive joints, and J i and J pi are Jacobian matrices relying these displacements to the displacement t. These matrices can be obtained by the differentiation of the global transformation matrix T i of the chain i including rigid, passive and elastic coordinates [START_REF] Pashkevich | Stiffness analysis of over constrained parallel manipulators[END_REF].

After several mathematical derivations extracted from [START_REF] Pashkevich | Stiffness analysis of over constrained parallel manipulators[END_REF], the displacement p i of the passive joints of leg i can be expres- sed as a function of t [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF] and it can finally be shown that the relation between q and t can be written as:

 q qJ δt [START_REF] Altintas | Manufacturing automation, metal cutting mechanics, machine tool vibrations and CNC design[END_REF] where the expression of J q is detailed in [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF]. Introducing this relation into (17) will lead to:

  22 1 2 cos T l l l Tt     r δt M δt ( 20
)
where the expression of M r is given at [START_REF] Cammarata | The elastodynamics of a 3-CRU spherical robot[END_REF]. As a result, from [START_REF] Cammarata | The elastodynamics of a 3-CRU spherical robot[END_REF], finding the robot natural frequencies relies on solving the 6 dimensional eigenvalues problem

  21 6 det 0 l    rr I M K (21)

Reduced Elasto-Dynamic Model

To conclude this section, it is necessary to mention that, from the previous obtained results, the robot free dynamic behavior can be modeled using the following expression:  rr M δt K δt 0 [START_REF] Merrit | Theory of self-excited machine-tool chatter[END_REF] As a result, if external milling forces F are applied on the tool, [START_REF] Merrit | Theory of self-excited machine-tool chatter[END_REF] becomes:  rr M δt K δt F [START_REF] Tlusty | Basic non-linearity in machining chatter[END_REF] This expression will be used in the following sections in order to analyze the robot behavior during milling process. It should be mentioned that in order this model to be valid, the displacements of the robots should be sufficiently slow so that the vibratory phenomenon are mainly due to the milling process excitations. Let us now define the expressions for the computation of the milling forces.

CUTTING FORCE MODEL

In milling process, the cutting force F that appears in equation ( 23) is caused by the interaction between the tool and the workpiece. It is a contact force and it is distributed along the affected area of the tool cutting part [START_REF] Altintas | Manufacturing automation, metal cutting mechanics, machine tool vibrations and CNC design[END_REF]- [START_REF] Jayaram | Estimation of the specific cutting pressures for mechanistic cutting force models[END_REF]. In robotic machining, the tool is mounted on the robot end-effector. In that case, the cutting force serves as a source of an additional external loading for robot and influences its motion. To evaluate this influence and to correctly analyze the robot behavior while machining, model of cutting force should be defined.

Basic Expressions

Most of existing works in the area of machining are based on Merchant's model of cutting, where he assumed that the tool/workpiece contact force is distributed uniformly along the tool cutting edge [START_REF] Merchant | Mechanics of metal cutting process. I-Orthogonal cutting and type 2 chip[END_REF]. Based on this hypothesis, a complex three-dimensional process of cutting can be analyzed referred to the one separated tool cross-section. In this cross-section the distributed cutting force can be presented as an equivalent point force F c applied to the cutting edge.

In general, the cutting force F c has a nonlinear nature and depends on many factors such as cutting conditions, properties of workpiece material and tool cutting part, etc. [START_REF] Merrit | Theory of self-excited machine-tool chatter[END_REF], [START_REF] Tlusty | Basic non-linearity in machining chatter[END_REF]. But, for given tool/workpiece combination, the force F c could be approximated as a function of an uncut chip thickness h, which represents the desired thickness to cut.

The typical diagram F c (h) is presented in Fig.

3 [START_REF] Armarego | The machining of metals[END_REF]. To describe it analytically, different models (linear, exponential, fractional, etc.) can be found in literature. The linear model perfectly suits to conventional CNC based machining (especially with single-tip tools). Here, usually the high rate removal of working material produces large chip thickness and the relation F c (h) could be correctly approximated with a linear function [START_REF] Altintas | Manufacturing automation, metal cutting mechanics, machine tool vibrations and CNC design[END_REF], [START_REF] Armarego | The machining of metals[END_REF] (Fig. 3-a):

    ,0 c T E F h k h k L h    (24)
where k T is the so called tool/workpiece cutting energy, k E is the tool/workpiece edge force, L is the width of cut. The model parameters k T , k E are estimated experimentally for a given combination of tool/working material. However, in robotic based machining the robot compliance could be the source of considerable relative tool/workpiece displacement. Even, the loss of tool/workpiece contact can be observed. As a result of such behavior, the material removal rate and the corresponding cutting force depend on current position of the tool on its path. So, the linear model F c (h) cannot be applied correctly and some nonlinear approximations are required (exponential, fractional, etc.).

The exponential model [START_REF] Jayaram | Estimation of the specific cutting pressures for mechanistic cutting force models[END_REF] corresponding to Fig. 3-b

is based on the equation   ,0 m cF F h k h L h  ( 25
)
where k F is the so called specific cutting force and the power m<1 depends on the properties of the workpiece and the tool cutting part. The parameters k F , m are also estimated experimentally for a given combination of tool/working material. This model is quite popular in industrial applications with machines CNC to analyze different machining operations with multi-edge tools [START_REF] Peigne | Impact of the cutting dynamics of small radial immersion milling operations on machined surface roughness[END_REF]. But the exponential approximation does not fit well the physical phenomena of cutting process in case of small removal of material (i.e. for small h). The main reason for this is that the function [START_REF] Peigne | Impact of the cutting dynamics of small radial immersion milling operations on machined surface roughness[END_REF] has infinite derivative when h tends to zero, which does not correspond to the reality.

In order to avoid this drawback, the fractional model [START_REF] Brissaud | The Fractional Model for the Determination of the Cutting Forces[END_REF] was proposed, that corresponds to Fig. 3-c and is based on the expression

    2 0 ,0 1 ss c s h h r h h F h k L h hh    (26)
where 0 1 r k k   depends on the parameters k ∞ , k 0 that define the so called stiffness of the cutting process for large and small h respectively (see Fig. 3-c) and h s is a specific chip thickness, which depends on the current state of the tool cutting edge. The parameters k 0 , h s , r are evaluated experimentally for a given combination of tool/working material.

In this work, the cutting force will be computed using the fractional model [START_REF] Brissaud | The Fractional Model for the Determination of the Cutting Forces[END_REF]. To take into account the latter phenomena, it is allowed for h to be either positive or negative, assuming that

  0, if 0 c F h h  (27)
For multi-edge tool the machining surface could be formed by means of several edges simultaneously. The number of working edges varies during machining and depends on the relative tool/workpiece position. Thus, the total force F c of such interaction is a superposition of forces F c,i generated on each tool edge i, which is currently in contact with the workpiece. Due to the presence of two different types of tool motion (spindle rotation, feed) the contact force F c,i can be described with its radial F r,i and tangential F t,i components. In accordance with Merchant's model [START_REF] Merchant | Mechanics of metal cutting process. I-Orthogonal cutting and type 2 chip[END_REF], the t-component of cutting force F t,i can be computed with the equation [START_REF] Brissaud | The Fractional Model for the Determination of the Cutting Forces[END_REF]. The r-component F r,i is related with F t,i by following expression [START_REF] Laporte | Phenomenological model for drilling operation[END_REF] ,,

r i r t i F k F  ( 28
)
where the ratio factor k r depends on the tool/workpiece characteristics.

In robotic machining it is more suitable to operate with forces expressed in the robot tool frame {x,y,z}. Then, the corresponding components F x , F y of the machining force F c are expressed as follows ,, 11 ,, 11 cos sin sin cos

zz zz nn x r i i t i i ii nn y r i i t i i ii F F F F F F           (29)
where n z is the number of currently working cutting edges, φ i is the angular position of the i-th cutting edge (the cutting force in z direction F z is negligible here) It should be stressed that the cutting force components F r,i , F t,i mentioned in equation (29) are computed regarding the corresponded chip thickness h i , which should be also evaluated. Using mechanical approach of analysis of machining operation, the parameter h i is computed in general as the geometrical distance between the position of the given tooth i and the current machining profile. Let us suppose that the feed is applied only in x direction (robot base frame). Below, the instant t=0 corresponds to the tool position when one of the teeth is in contact with the working material. Also, it will be supposed that the feed rate v f and the spindle rotational speed Ω are constant during whole machining process.

h i v i a i i  i i F r,i F t,i F x F y v f W TCP Workpiece Tool a i v i v f v Wx v Wy

Case A: Rigid Tool Fixation

The objective of this case study is to understand the mechanics of the tool/workpiece interaction, while any dynamic aspect related to the robot compliance is excluded from the analysis. In that case the applied feed rate and spindle rotational speed totally determine the position {x i , y i }, 1, z iN  of all tool cutting edges referred to robotic base frame at each instant of machining t

   cos , sin 2 1 , 1, i f i i i z iz x v R y R N i i N            W    ( 30
)
The orientation α i of the i-th tooth velocity v i can be defined by the feed rate v f and the spindle rotational speed Ω at each instant as to the i-th tooth, when it is situated inside the working material.

    atan2 i f x y v +v v a WW  with sin , cos x i y i v R v R  WW  W  W . The
If the given tooth is located outside the working material, the corresponding chip thickness h i is equal to zero. The following expressions allow evaluating h i for all possible positions of the i-th tooth on its path while machining 0, sin , 2 , 1, 2 sin , 2

i z i i i i f z f z i i f z xR h x R x v N i N v N x v N a  a        W    W  W  (31)
The advantage of the presented algorithm of computing the chip thickness is that different phases of tool/workpiece interaction illustrated in Fig. 5 can be identified.

It should be mentioned that the phase of tool approaching to the workpiece corresponding to the zero machining force is not considered here. For the remaining phases a detailed analysis is presented below:  The phase of tool engaging into the working material (phase I) corresponds to the variable contact area between the tool and the workpiece. The TCP during this phase is located always outside the workpiece. In fact, the phase I can be divided into two sub phases:  Phase Ia: In the beginning of milling operation a small area of workpiece is affected by the machining process. This fact and presence of two types of motion (feed, spindle rotation) form the case, when only one tooth can participate in cutting at the same time (Fig. 5-a). Such behavior produces intermittent machining forces F x and F y with the frequency 2 z N  W Hz. The sub phase 1 is very limited in time and its duration depends on the feed rate, the spindle rotational speed and the number of teeth N z . For example, if v f =4m/min, Ω=10 4 rpm, N z =4 the duration of phase Ia is only 0.04sec.  Phase Ib: It corresponds to a case, when several teeth can participate in cutting at the same time, but the TCP does not reach the workpiece border (Fig. 5-b). As a result an oscillatory periodic behavior in machining forces F x and F y is observed. But, because of different number of currently working teeth, the force patterns are not homogenous.  The phase of machining with fully engaged tool (Phase II) starts when TCP reaches the workpiece border (Fig. 5-c). In that case always the same number of teeth (n c =2 for the tool with N z =4) is working at every instant of cutting process. It produces harmonic periodic machining forces F x and F y with the frequency 2 Z  W Hz. So, even with the analysis of tool motion during machining without any dynamic aspect, the oscillatory behavior of machining forces can be detected for whole process. The high frequency of such oscillation (for example, Ω=10 4 rpm, N z =4 give frequency of 667Hz) does not affect the motion of robot but it can be crucial for robot control system and should be considered in design of robotic machining process.

Case B: Tool Fixation with Compliance in x Direction

In contrast to Case A, the dynamic aspect of tool motion associated with the robot compliance is considered here. Thus, at each instant of machining, the position x TCP is defined as a superposition of tool displacement x f due to feed and a dynamic displacement δx due to compliance of the fixation:

TCP f x x x  . The first component ff xv 
 is known at each instant of process while the second one depends on the current position of the tool regarding to the machining profile. In this case, the dynamic displacement δx can be obtained by reducing the equation ( 23) to a one-dimensional problem and by introducing the damping related to the machining process and robot control algorithms

x x x M x C x K x F       ( 32
)
where M is the equivalent mass of the tool fixation, K x and C x are its stiffness and damping respectively and F x is the machining force in x direction. The damping 2 xx c k m   is related to the damping factor ζ, which can be estimated experimentally.

{x i ,y i } τ TCP τ TCP τ-Dτ  i τ-Dτ  i τ Da R {x i ,y i } τ-Dτ W A i Workpiece area h i Ds x Ds y Nodes x f + x Dτ τ-Dτ
z i TCP i i i x x R y R i N      .
Comparative analysis of this position with respect to the current machining profile defines the chip thickness h i removing by i-th tooth. But, the main issue here is to define whether i-th tooth participates in cutting for given instant of process. For this reason, it is proposed to create a mesh on the workpiece, where each node j ( 1, w jN  , N w is the number of nodes) can be filled with "1" or "0": "1" corresponds to nodes situated in the workpiece area with material, "0" corresponds to nodes situated in workpiece area that was cut away.

In order to define the number of currently cut nodes by the i-th tooth, the previous instant of machining process should be considered. Let us define A i as an amount of working material that is currently cut away by the i-th tooth. So, if node j filled with "1" is located inside the sector, it changes to "0" and A i is increasing by xy ss DD (Δs x , Δs y are node steps in x and y directions respectively). Analyzing all potential nodes and computing A i , the chip thickness h i , removed at given instant of process by the i-th tooth, can be estimated by ,

i i i h A R a D 1, z iN 
. The angle Δα i determines the current angular position of the i-th tooth regarding to its position at the instant τ-Δτ (Δτ is the time step) and referred to the position of TCP at τ-Δτ.

Here, in contrast to the previous case, the dynamic aspect of the tool motion allows  estimating of deviation in tool motion from the desired one because of the robot compliance in the feed direction. It should be mentioned that this deviation affects the Cartesian stiffness of robot but does not influence the machining profile quality. For example, following parameters M=100 kg, K x = 5 3 10  N/m, ζ=0.05, N z =4, Ω=10 4 rpm, v f =4 m/min provide deviation in the feed direction of 0.13mm (Fig. 8). But it is not essential for this application.  detecting of vibratory behavior of tool motion while it is engaging into the workpiece. In some cases the low frequency of such motion can excite robot natural frequency, destabilize machining operation and even damage the tool or/and workpiece. For example, the milling process with following parameters M=100 kg, K x = 5 3 10  N/m, ζ=0.05, N z =4, Ω=10 4 rpm, v f =4 m/min generates vibration of f 1 =8.8Hz from the beginning of cutting process in addition to the frequency 667Hz related to the spindle rotation. More details on influence of the tool fixation parameters (M, K

x ) on the dynamics of its motion during machining are presented in Table 1, which covers the range of values for M, K x computed for different configurations of the robot KUKA KR 240 using equations presented in Section 2.

As it can be observed from the table, changing the fixation parameters (i.e. the robot configuration) influences low frequencies (about 10 -20 Hz) of the tool motion. Suitable robot configurations to perform given machining operation could be defined. But, it should be mentioned that this one dimensional equivalent model presented here allows analysis of machining process dynamics in the feed direction only. In order to evaluate behavior of the tool motion more closely to the real machining operation, this model should be extended.

Case C: Tool Fixation with Compliance in x and y Directions

In this case, a dynamic aspect of tool motion in feed direction (x) and orthogonal to it (y) is considered. Then, similarly to the Case B, at each instant of machining process, the position of TCP is determined by ,

TCP f TCP f x x x y y y      , where ,, , f f x f f y x v y v   .
The dynamic displacements δx, δy could be obtained from the equation ( 23) which, in this case, is reduced to

x y F x x x F y y y                               rr Μ C Κ (33)
where M r (2×2) is the equivalent mass matrix of the tool fixation, the matrices K r (2×2) and C (2×2), where , , , , 2 i j i j i j i j   C K M characterize the fixation stiffness and damping respectively (which can be estimated experimentally),

F

x and F y are the machining forces in x and y directions. In contrast to the previous case, the position of the i-th tooth at each process instant t includes dynamic components in both directions: sin ,

i TCP i x x R   cos , i TCP i y y R   1, z iN 
. The algorithm of computing the chip thickness h i for given position of tooth {x i , y i } is similar to Case B. In order to illustrate the advantages of this two dimensional model and its ability to detect some phenomena (that are not visible in Cases A, B) the robotic milling process is simulated for KUKA KR 240 robot with the following parameters

N kg m                     r r M K C
.

Here, the equivalent mass matrix M r is computed in accordance with the method presented in the Section 2 of this paper. The stiffness K r is the structural stiffness of the robot, referred to its end-effector. It should be noted that in practice the non diagonal elements in these matrices are non zero. But, with this simplified case it is possible to identify qualitatively the dynamic nature of the tool behavior in the direction (y) orthogonal to the feed. Simulation results corresponding to this case study are presented in Fig. 9.

It should be also stressed that the tool dynamic behavior in the feed direction (x) is similar to the results, which were obtained in the Case B. The displacement in y-direction has an essential dynamic component during the phase of tool engagement (Phase I) into the workpiece and becomes constant while machining with fully engaged tool. The corresponding frequency f 1y =7.3Hz is comparable with the frequency f 1x =8.2Hz of the tool dynamic displacement in feed direction. More details on the tool motion during milling process regarding to the parameters of the tool fixation are presented in It should be noted that changing the robot configuration affects the tool dynamics in y direction which is crucial regarding to the quality of the final product. Hence, considering the tool displacement orthogonal to the feed direction is essential and the Case C gives more realistic results comparing Cases A and B. As it is shown in Fig. 9, the deviation (0.31 -0.56mm) in machining profile from the desired one has a vibratory behavior during the phase of tool engagement into the workpiece. Thus it cannot be suppressed by straightforward compensation methods and other compensation techniques should be proposed.

CONCLUSION

In this paper, a reduced elasto-dynamic model of the robotic based milling process has been presented. In contrast to previous works, it takes into account the interaction between the milling tool and the workpiece that depends on the endeffector position, process parameters and cutting conditions (spindle rotation, feed rate, geometry of the tool, etc.). To reduce the dimension of the problem, the robot dynamics have been described as an equivalent mass-spring-damper system with six dimensions. This approach aims at decreasing computational cost and at avoiding inaccuracy due to illconditioning in the full size model. To achieve a realistic modeling of the milling process, the machining efforts due to the interaction between robot tool and working material have been introduced into the robot model and calculated at each time instant. This model has allowed selecting the best process parameters and avoiding the vibratory behavior of this machining system which can dramatically affect the milling quality.

The developed model has been applied to the elastodynamic behavior analysis of KUKA KR240 robot used for milling of an aluminum workpiece for automobile industry. This allowed us estimating the deviation in motion of the robot end-effector during machining caused by the flexibilities in robot links and joints, which essentially influence the quality of the final product.

Figure 1 :

 1 Figure 1: Displacements and elastic deformations of a beam.

Figure 2 :

 2 Figure 2 : Architectural representation of a general parallel manipulator

Figure 3 :

 3 Figure 3: Linear (a), exponential (b) and fractional (c) cutting force models Fc(h).

Figure 4 :

 4 Figure 4: The force interaction between the i-th tooth and workpiece and corresponding kinematics.Let us present algorithms for the computation of h i during process of robotic milling for three different types of tool fixation, assuming that machining is performed in aluminum workpiece with straight borders. The cutter of the external radius R=10mm with N z =4 teeth distributed uniformly over the

Figure 5 :

 5 Figure 5: Different phases of tool/workpiece interaction (case Nz=4) and corresponding machining forces; (a)process of tool engaging into the workpiece, only one tooth can be in contact with workpiece at the same time, (b)process of tool engaging into the workpiece, several teeth could be in contact with workpiece at the same time, (c)machining with fully engaged tool.

Figure 6 :

 6 Figure 6: Machining force patterns with average forces referenced to the frame {x,y}.

Figure 7 :

 7 Figure 7: Evaluating the tool/workpiece intersection Ai and computing the corresponding chip thickness hi.

Figure 8 : 5 3 × 10 N

 8510 Figure 8: Machining force patterns and TCP dynamic displacement in case of 1DOF model (M=100 kg, Kx=

Figure 9 :

 9 Figure 9: Machining force patterns and TCP dynamic displacements in case of 2DOF model (Mr,xx=Mr,yy=100 kg, Kr,xx=Kr,yy=

Table 1 : Influence of tool fixation parameters on the tool motion; τs is the settling time, PO is the overshoot, f1 is the first frequency of the tool dynamic displacement in feed direction; ζ=0.05, Nz =4, Ω=10 4 rpm, vf=4 m/min

 1

	M, kg K x , N/m					τ	s , sec |x	s |, mm PO, %	f 1 , Hz
	100	0.05 10 	6	1.2	0.80	52	3.5
	100	0.30 10 	6	0.6	0.13	30	8.8
	100	0.60 10 	6	0.5	0.07	23	12.4
	100	1.00 10 	6	0.4	0.04	17	15.6
	100	2.00 10 	6	0.4	0.02	11	22.6
	150	2.00 10 	6	0.5	0.02	14	18.3
	200	2.00 10 	6	0.5	0.02	18	15.6

Table 2 : Influence of the tool fixation parameters on tool dynamic behavior; f1x, f1y are first frequencies of the tool dynamic displacement in x and y directions respectively; Nz =4, Ω=10 4 rpm, vf=4 m/min M r,xx =M r,yy , kg

 2

		K r,xx =K r,yy ,	|x	s |,	|y	s |,	f 1x ,	f 1y ,
		N/m						mm	mm	Hz	Hz
	100	0.05 10 	6	0.80	2.61	3.2	2.2
	100	0.30 10 	6	0.13	0.43	8.2	7.3
	100	0.60 10 	6	0.07	0.22	12.4	11.4
	100	1.00 10 	6	0.04	0.13	15.6	14.6
	100	2.00 10 	6	0.02	0.06	22.4	21.5
	150	2.00 10 	6	0.02	0.06	18.3	17.4
	200	2.00 10 	6	0.02	0.06	15.6	15.6

Table 2 (

 2 it covers the range of values for M

	r,xx , M	r,yy , K r,xx , K r,yy

ACKNOWLEDGMENTS

The authors would like to acknowledge the financial support of the ANR, France (Project ANR-2010-SEGI-003-02-COROUSSO) and the Region "Pays de la Loire", France.