G Coffignal

P Lorong

L Illoul

Pimm

A general method to accurately simulate material removal in virtual machining of flexible workpieces

Keywords: progressive Boolean intersection deformable bodies, virtual machining, non-rigid workpiece, vibrations, machined surface, Boolean chip, material removal simulation

Multi-axis milling and other computer numerical control machining processes allow us to create very complex geometries and thin parts. In this context, virtual machining is a powerful tool, but the simultaneous vibrations of the tool and workpiece are not easy to define and take into account. This paper presents a general method with which to simulate material removal when both the workpiece and tool are assumed to be non-rigid. We consider that they both vibrate when we define the Boolean chip. This is not usually considered with the aim of predicting the machined surface vibrations and the resulting geometric defects including roughness. By extending the material frame associated with the non-rigid workpiece, our method precisely defines the material removal for any tool or workpiece. It then allows us to establish a method of deriving efficient numerical approximations with which to simulate a succession of machining operations from roughening to finishing. Two kinds of finite element approximations are linked. One is a classical elastic finite element model including damping. The second, which is kinematically linked to the first, accurately describes the relative motion of each part of the tool with respect to the workpiece, and ensures the description of material removal and related forces. Two industrial examples show the potential of the method.

Introduction

Multi-axis milling and other computer numerical control (NC) machining processes allow the creation of very complex geometries and thin parts. In this context, the optimization of tool paths, workpieces, machines and tools must avoid the occurrence of unwanted vibrations. Before considering vibrations, it must be noted that even when considering purely geometric approaches, the works of Lee and Nesler [START_REF] Won | Complete swept volume generation, part i: Swept volume of a piecewise c1-continuous cutter at five-axis milling via gauss map[END_REF][START_REF] Won | Complete swept volume generation, part ii: Nc simulation of self-penetration via comprehensive analysis of envelope profiles[END_REF][START_REF] Won | Virtual workpiece: Workpiece representation for material removal process[END_REF] show that the need for accurate and robust descriptions of the volume swept by a rigid cutting tool in a rigid workpiece still produces research work. By including the dynamics of the whole workpiece-tool-machine system (WTMS), virtual machining (VM) aims to accurately predict phenomena occurring during any machining operation and thus to explain and master the phenomena.

Here, by VM we mean the use of a numerical process and dedicated software that should allow simulation in the time domain of a wide range of machining processes, machines, workpieces, tools and tool paths, taking into account numerical control of the machine, vibrations and material removal in an integrated, coupled and general approach based on physical models at a convenient scale. This modeling should be predictive. With this aim, VM software would be able to model and simulate any kind of machining operation, including the machining of a free-formed surface or thin-walled workpiece and five-axis milling.

Chatter and discontinuous cutting are among the important causes of vibrations, and since works by Tobias [START_REF] Andrew | A critical comparison of two theories of machine tool chatter[END_REF], many authors have been working on the subject. Since early works, it was clear that simulation programs must model regenerative forces. Their definition was recalled by Smith and Tlusty [START_REF] Smith | Efficient simulation programs for chatter in milling[END_REF]: "Regenerative Force means that the force on any tooth in the cut depends not only on the feed per tooth, and on the deflection of the cutter, but also on the surface which was left by the passage of previous teeth". In this context, many cutting force models have been described [START_REF] Smith | An overview of modeling and simulation of the milling process[END_REF]. They employ the concept of uncut chip thickness h, which corresponds to an idealized vision of the tool as a three-dimensional (3D) eraser (Fig. 1) and is valid in a wide range of relative tool/workpiece motions. This corresponds to simulations at a macroscopic level.

Altintas et al. gave an overview of research work dealing with the dynamics of milling and grinding [START_REF] Altintas | Chatter stability of metal cutting and grinding[END_REF] and a description of the concept of the virtual machine tool in [START_REF] Altintas | Virtual machine tool[END_REF]. All these topics were surveyed in [START_REF] Altintas | Manufacturing Automation[END_REF] by Altintas. Brecher et al. [START_REF] Brecher | Interaction of manufacturing process and machine tool[END_REF] presented the state of research on process-machine tool interaction in 2009. It must be noted that most of the works that deal with chatter and discontinuous cutting have the objective of fast stability analysis in the frequency domain or in the context of a periodic solution. This implies that some simplifications, linearizations and dedicated analytical developments are necessary to obtain delayed differential equations accounting for regenerative forces and then, obtain solutions to describe the stability limits. Time-domain simulations do not require as many simplifications in the models, and thus have the potential to follow more accurately the physical machining process. In return, they are more time consuming.

Research works, for the most part, have considered vibrations of the workpiece to be negligible and have ignored them for simplicity. When the dynamic behavior of the workpiece cannot be neglected, which is now often the case because of component mass minimization, the dynamic behavior must be accounted for in the modeling; this is a current research topic.

As the system has multiple degrees of freedom (DOF) the dynamics is governed by the natural frequencies and mode shapes. In addition, because of material removal, there are continuous changes in the distribution of mass and rigidity that modify the multi-DOF system. As a consequence, natural frequencies and mode shapes of the workpiece change along the tool path and this affects the frequency response function (FRF) at a given location. Altintas et al [START_REF] Altintas | Dynamic peripheral milling of flexible structures[END_REF] and Corduan et al. [START_REF] Corduan | Experimental approach of milling stability of thin walled parts, comparison with time domain simulation[END_REF] developed such models of a flexible workpiece in interaction with rigid tools. Coffignal et al. [START_REF] Assouline | Numerical simulation of machining at the macroscopic scale: dynamic models of the workpiece[END_REF][START_REF] Lorong | Contribution of computational mechanics in numerical simulation of machining and blanking: state of the art[END_REF][START_REF] Coffignal | Virtual machining: a general approach to deal with flexible workpieces[END_REF] also included the possibility of tool deformations in their time-domain approach, but for a simple geometry of the workpiece and three-axis milling. Bravo et al. [START_REF] Bravo | Stability limits of milling considering the flexibility of the workpiece and the machine[END_REF] accounted for the workpiece dynamics and introduced a 3D lobe diagram to describe the evolution of the stability limits along the tool path. Alan et al. [START_REF] Salih Alan | Analytical prediction of part dynamics for machining stability analysis[END_REF] studied the change in part dynamics and Budak et al. [START_REF] Budak | Prediction of workpiece dynamics and its effects on chatter stability in milling[END_REF] predicted its effects on chatter stability in milling. The successive FRF corresponding to both the change in position and the modification of the model is efficiently calculated using a method proposed by Özgüven [START_REF] Özgüven | Structural modifications using frequency response fonctions[END_REF]. Arnaud et al. [START_REF] Arnaud | Simulation of low rigidity part machining applied to thin-walled structures[END_REF] gave an example of a time-domain simulation using an evolutive finite element (FE) model of the workpiece. Kersting et al. [START_REF] Kersting | Simulation concept for predicting workpiece vibrations in five-axis milling[END_REF] presented a comparison of experimental results and simulations for a thin vibrating plate in five-axis milling. In [START_REF] Eksioglu | Discrete-time prediction of chatter stability, cutting forces, and surface location errors in flexible milling systems[END_REF] Eksioglu et al. presented a general discrete-time model of the milling system. This model was used to solve the chatter stability and take into account the dynamics of the tool and workpiece, including interrupted contact.

Except for [START_REF] Budak | Prediction of workpiece dynamics and its effects on chatter stability in milling[END_REF] and [START_REF] Lorong | Contribution of computational mechanics in numerical simulation of machining and blanking: state of the art[END_REF][START_REF] Coffignal | Virtual machining: a general approach to deal with flexible workpieces[END_REF][START_REF] Kersting | Simulation concept for predicting workpiece vibrations in five-axis milling[END_REF], it appears that these research works consider simple geometries of workpieces where the tool-workpiece interaction can be modeled in parallel planes. In [START_REF] Lorong | Contribution of computational mechanics in numerical simulation of machining and blanking: state of the art[END_REF][START_REF] Coffignal | Virtual machining: a general approach to deal with flexible workpieces[END_REF][START_REF] Kersting | Simulation concept for predicting workpiece vibrations in five-axis milling[END_REF], the interaction is not restricted to such planes but the workpiece is a thin plate.

Time-domain five-axis milling simulations of the machining of turbine blades, including tool and workpiece vibrations and regenerative forces, were given by Biermann et al. [START_REF] Biermann | A general approach to simulating workpiece vibrations during five-axis milling of turbine blades[END_REF]. Recently, Lorong et al. [START_REF] Lorong | Simulation of a finishing operation: Milling of a turbine blade and influence of damping[END_REF] presented such an example of simulation using the method described in the present paper. FE models of the tool and workpiece are used in both approaches but the tool-workpiece interactions are based on very different methods.

The method used in [START_REF] Kersting | Simulation concept for predicting workpiece vibrations in five-axis milling[END_REF] and [START_REF] Biermann | A general approach to simulating workpiece vibrations during five-axis milling of turbine blades[END_REF] is based on a method described by Wienert et al. [START_REF] Weinert | Modeling regenerative workpiece vibrations in five-axis milling[END_REF][START_REF] Weinert | Modeling of surface structures resulting from vibrating milling tools[END_REF] to model regenerative workpiece vibrations in milling. It allows the construction of a realistic geometric model of the surface structures by means of a discrete set of surface points with associated normals. The undeformed chip domain is obtained using a recursive constructive solid geometry (CSG) model. The corrugated shape of the chip and associated cutting forces due to vibrations are then obtained by means of corrections using an evaluation of the toolworkpiece relative displacements. An incremental numerical scheme applied to the delayed differential equations of motion gives the relative displacement.

The literature shows that it is possible to build predictive models that account for flexible workpieces and are able to well describe the actual physical observations. These include the structure and the location errors of the final surface. Because it is difficult to describe the geometric interaction between the tool and workpiece when neither of their vibrations can be neglected, the authors developed models that, even if they are often able to describe a class of problems, cannot deal with general situations. It appears that a definition of what we call the "Boolean chip" (Fig. 1b), which is considered in all simulations at the macroscopic level to define h, is not usually given for workpieces that undergo time-varying deformations. This does not allow a detailed understanding of all the assumptions that are made in describing the machined surface, and their consequences in terms of range of application of the methods.

It is not a straightforward task to build and mix dynamic models in efforts to simulate general machining operations while accounting for complex tool paths and flexible workpieces of any shape.

In this context, the present paper focuses on a method that we proposed and developed to model efficiently the relative motions and associated material removal in time-domain simulations. The originality of the method is that by following the definition of the Boolean chip in a workpiece that undergoes deformation, the method provides an easy way to take into ac-count the vibrations in a wide range of machining situations (e.g., milling, turning, reaming, and boring) including multi-axis scenarios. It also allows us to follow the complex history of the incremental removal of matter and accurately calculate the regenerative forces without any specific analytical development for each situation. All other components of the method have been demonstrated to be suitable in the literature and are present in the method proposed by Biermann et al. [START_REF] Biermann | A general approach to simulating workpiece vibrations during five-axis milling of turbine blades[END_REF]. The next section (section 2) briefly outlines our approach and introduces the domain swept by the rake faces in the material domain of the workpiece that defines the Boolean chip. The section describes FE models, step-by-step integration and modeling of cutting forces.

Section 3 defines the model of material removal, which usually is only clearly defined in the case of a perfectly rigid workpiece. A definition of the Boolean chip that is valid when the workpiece undergoes deformations is derived. With this aim, a material frame R m and its associated material coordinates are introduced. It is shown that an extension of R m outside the workpiece must be defined. In the present approach, a class of one-to-one mappings Φ + is introduced to make this extension. These mappings allow the definition and accurate extension of the usual model of material removal once a kinematic model is chosen for the workpiece deformation. Procedures with which to obtain the geometry of the Boolean chip are then given for any choice of Φ + . Section 4 gives an example of a practical implementation of the approach: the one-to-one mapping Φ + is built by means of an FE description. A dexel model [START_REF] Van Hook | Real-time shaded nc milling display[END_REF] of the workpiece is chosen and described. The model follows the vibrations and deformations of the FE model of the workpiece. The BREP model of the domain swept by the rake faces and its intersections with the dexel model of the workpiece, which gives the model of the Boolean chip, are described. This finally allows the calculation of the cutting forces and machined surface, time step after time step, and a description of this surface in terms of its roughness, undulations and form defects.

The final section gives two examples of the machining of industrial workpieces that demonstrate possibilities of the method. Conclusions are then presented.

Outline of our approach

Today, the FE method and associated algorithms, including reduction and time integration, have reached a high level of maturity and versatility when considering vibration analysis of elastic bodies. Commercial software can manage efficiently and accurately static and modal linear analyses of very large models. The method that we propose as a part of what could be a complete VM system was designed to take advantage of the power and robustness of FE solvers and their pre-and post-processors. As shown in Fig. 2, the problem is thus described by geometric and mechanical models as available in commercial FE software. In addition to these models of the tool and the workpiece, detailed geometric models of the active parts of the tool and their associated cutting force models are provided and the NC program must be given. These data, and the choice of the time-domain integration scheme, are sufficient to be able to simulate the operation.

Our method is implemented in Nessy, which is research software that was developed in our laboratory. There is no need for any analytical preprocessing depending upon the tool or the tool path. The small time step ∆t of the integration scheme allows the simulation of the process. This is done time increment after time increment without any requirement other than to follow the relative motion imposed by the NC program (APT) and maintain the dynamic equilibrium of the WTMS as described by Eq. ( 1). Additional models are automatically built and linked within Nessy to a reduced model described by Eq. ( 5), and the interactions and material removal are modeled independently of the mesh used in the FE software to build Eq. ( 1).

In the following, we only consider the modeling of the mechanics of the interaction between the tool and workpiece and the resulting vibrations. To accurately describe the relative motions of the deformable workpiece and cutting edges of the tool, it is necessary to introduce several frames. A coordinate system is assumed to be associated with each frame R i we introduce. To simplify the presentation, we use a unique notation R i for the frame and its associated coordinate system. R g denotes the inertial frame, R s is the frame in which the workpiece is supported and undergoes small deformations and rotations, and R c denotes the frame in which the cutter undergoes small deformations and rotations. The coordinate system associated with each of these three frames is assumed to be Cartesian.

The motions of R s and R c are supposed to be known, and we denote by R s (τ ) and R s (τ ) the known positions of these frames and associated origins and Cartesian axes in R g at time τ .

In milling, the model of the spindle is included in the FE model, which is described in R c by the mesh M c . In turning, the model of the spindle is included in the FE model, which is described in R s by the mesh M s . In both cases, there is a large rotation that can be described by means of the rotation vector θ c/s of R c with respect to R s . Fig. 2 describes the main steps of our approach. Several models are used and some are linked with others. Fig. 3 gives an outline of these links. Meshes M c and M s (section 2.1) are used to build all the FE matrices for the dynamics. Three-node triangles (T3) are used to describe the boundary displacements and positions of the active parts of the tool by means of M Σ (section 4.2) and the swept domain during a time step by means of M (section 4.2). Four-node tetrahedrons (T4) are used in M + m (section 4.1), which allows us to easily create a one-to-one mapping Φ + . The displacements of the nodes of M Σ , M and M + m are constrained to follow the displacements of their respective parent models: these displacements are given by the interpolation functions and nodal displacements of M c for the active parts of the tool (meshes M Σ and M) and M s for the one-to-one mapping (mesh M + m ). A dexel model D (section 4.3) allows us to follow, in a material frame and in detail, the geometrical changes of the workpiece, as machining develops.

FE model used to simulate the dynamics of the WTMS

Our approach to virtual machining merges in a unique solver-the FE simulation of the dynamics of the WTMS and the modeling of material removal. This implies the modeling of the geometry of the tool and the transformation of the domain Ω W describing the workpiece geometry. For efficiency, the FE matrices and the reductions are built in external commercial FE software and then imported. The FE method employs a matrix system Eq. ( 1). The construction of the lhs in Eq. ( 1) follows the usual FE approach ( [START_REF] Gmür | Dynamique des Structures. Analyse modale numérique[END_REF]) for the workpiece and its supports (index s) and the cutter (index c).

Eq. ( 1) is obtained using two FE meshes: M s in R s for the workpiece and M c in R c for the tool. Any type of element can be used to set up Eq. ( 1). If there are any non-rigid supports, they are included in the model. q s and q c are columns having degrees of freedom:

M s • •• q s + Ds • • q s + Ks • q s = Q sT + Qs , (1a) M c • •• q c + Dc • • q c + Kc • q c = Q cW + Qc . ( 1b
)
In this paper, for simplification we assume that material removal does not significantly affect the matrices on the lhs of Eq. (1a). We thus consider only one FE model of the workpiece and one associated mesh M s . Elements are not removed from the FE element model of the workpiece to continuously account for material removal in our method, but several meshes of the workpiece could be considered and combined. This was proposed in [START_REF] Assouline | Numerical simulation of machining at the macroscopic scale: dynamic models of the workpiece[END_REF][START_REF] Lorong | Contribution of computational mechanics in numerical simulation of machining and blanking: state of the art[END_REF], where a possible solution to account for these changes of mass and stiffness of the workpiece was given: macro time steps are defined to calculate and store FE matrices corresponding to the theoretical geometry of the workpiece (perfect machining) and to the current configuration of the WTMS at theses macro steps. Then, between two macro-time steps T a and T b , two solutions q a and q b are calculated using each set of matrices. The solution for τ ∈ [T a T b ] is approximated by q(τ ) = (1 -α)q a (τ ) + αq b (τ ) where α = (τ -T a )/(T b -T a ).

By means of Eq. (2), M s provides at any time τ the displacement u s and coordinates P s (τ ) of any point P of the workpiece in R s by means of H s , the displacement interpolation matrix for the workpiece. M c plays the same role for the tool, by means of H c , and leads to u c and coordinates C c (τ ) in R c of any point C of the tool. Associated virtual displacements are denoted δu s , resp. δu c . Initial positions at τ = t I are denoted C = C c (t I ) and P = P s (t I ) and are taken as reference positions; i.e., u s ( P, t I ) = 0 and u c ( C, t I ) = 0. Thus,

P s (τ ) = P + u s ( P, τ ), C c (τ ) = C + u c ( C, τ ), (2a)
u s ( P, τ ) = H s ( P) • q s (τ ), u c ( C, τ ) = H c ( C) • q c (τ ), (2b
)

δu s ( P) = H s ( P) • δq s , δu c ( C) = H c ( C) • δq c . ( 2c
)
On the rhs of Eq. (1), Q cW and Q sT link Eq. (1b) and Eq. (1a). Q cW , the generalized action of the workpiece on the tool, and Q sT , the generalized action of the tool on the workpiece, are constructed from the local cutting forces that model the interaction between the tool and workpiece across their common area of contact Σ C . The surface Σ C is not a simple one, particularly in milling (Fig. 4). It is a subset of Σ, the union of all rake faces of the tool, and very often, several teeth and/or inserts are involved in the contact at the same time τ . This leads us to consider a set of elementary tools T (k) in describing any insert or tool as an assembly of such elementary parts. We denote by f T /W the surface density of forces exerted on the workpiece by the tool at time τ , everywhere on Σ C ⊂ Σ. At the coincident point, the action on the workpiece is thus f W/T = -f T /W . Taking virtual works δu s T • f T /W and δu c T • f W/T associated with these densities, we have

Q sT = Σ C H sT • f T /W • dΣ, ( 3a
)
Q cW = Σ C H cT • f W/T • dΣ, ( 3b
)
where

f W/T , resp. f T /W , are the components of f W/T in R c , resp. of f T /W in R s . rake face Σ tooth cutting edge Γ tooth instantaneous surface of contact Σ c tooth
elementary tools zoom (a) A tooth, its rake face Σ tooth ⊂ Σ and an example of elementary tools.

Γ c tooth Σ c tooth F W/T (k) rake face cutting edge Γ tooth ⊂ Γ h (k) Γ (k) cutting edge of T (k) Σ (n) Σ (m) Σ (k) Σ tooth T (n) T (k) T (m) C F (k) Σ (k) : rake face of the elementary tool T (k) C 4 (k) C 1 (k) C 2 (k) C 3 (k) (b)
Close-up of rake faces Σ (i) for three elementary tools T (i) , i ∈ {m, k, n}. Details for T (k) .

Figure 4: End mill example. An elementary tool T (k) associated with a tooth. Rake face Σ tooth ⊂ Σ and rake faces Σ (i) of elementary tools;

Σ (i) ⊂ Σ tooth ⊂ Σ.
The dimensions of matrices and the number of unknowns in Eq. ( 1) can be very large if the workpiece or the tool is a real-life industrial object (of the order of 10 6 degrees of freedom in the examples given in the last section). When simulating machining, we only have to consider lowfrequency phenomena in the sense that only the few first modes of vibration make a significant contribution to the dynamic response of the WTMS. Typically, the number of modes kept in the reduced model is of the order of 10 for simple workpieces to hundreds for complicated shapes with cyclic symmetry.

This allows a huge reduction of the problem employing usual methods of reduction. In the frame of Rayleigh-Ritz methods [START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Géradin | Théorie des Vibrations[END_REF], this can be defined as

q i ≈ q R i (i = s or i = c) and q R s = Ψ s • y s , and δq R s = Ψ s • δy s , ( 4a
) q R c = Ψ c • y c , and δq R c = Ψ c • δy c , ( 4b
)
where rectangular matrices Ψ i have far fewer columns than rows. The dimension of y i is thus much smaller than that of q i . The introduction of the approximations Eq. ( 4) in Eq. ( 1) leads to the reduced system:

M s • •• y s + D s • • y s + K s • y s = Y s , (5a) M c • •• y c + D c • • y c + K c • y c = Y c . ( 5b
)
All the matrices in Eq. ( 5) are obtained with the same procedure used for M i in the case of the lhs and the same procedure used for Y i in the case of the rhs:

M i = Ψ T i • M i • Ψ i , Y i = Ψ T i • Q i . ( 6
)
For instance, the generalized cutting forces Y cW and Y sT are

Y cW = Ψ T c • Q cW and Y sT = Ψ T s • Q sT . ( 7
)
Matrices Ψ i usually gather the few first natural modes of vibration from Eq. ( 1) and may be improved [START_REF] Géradin | Théorie des Vibrations[END_REF] using well-chosen static deformations.

Any convenient FE software is used to get the required matrices, natural frequencies and mode shapes as input data, when building Eq. ( 5) in Nessy.

The main steps to get the constant matrices of the reduced model in Eq. ( 5) are shown in the upper part of Fig. 3: on the left is depicted the flowchart for the tool, taking the example of a milling cutter. The flowchart for the workpiece is shown on the right. Rectangles and ovals respectively represent processes and data.

Additional comments about the cutting forces and the Boolean chip

At the macroscopic level that we consider, the details of contact forces cannot be modeled because the non-linear and coupled thermomechanics are not fully considered during the formation of the chip; instead, a crude model, the accuracy of which is sufficient at this scale, is used. This well-established model when dealing with vibrations and chatter [START_REF] Andrew | A critical comparison of two theories of machine tool chatter[END_REF][START_REF] Altintas | Chatter stability of metal cutting and grinding[END_REF][START_REF] Altintas | Manufacturing Automation[END_REF] is based on the assumption that, to predict the cutting forces and the final surface, it is sufficient to consider the tool as a 3D eraser.

Considering the mathematical operation corresponding to this assumption, we define what we call a "Boolean chip" as opposed to the actual chip. Once this model is adopted for a vibrating workpiece, the last ingredient is an associated model that provides a sufficiently accurate estimation of the cutting forces, mainly from the local geometry of the Boolean chip. Mechanistic models [START_REF] Smith | An overview of modeling and simulation of the milling process[END_REF][START_REF] Altintas | Chatter stability of metal cutting and grinding[END_REF][START_REF] Altintas | Manufacturing Automation[END_REF] or any other kind of representative model at this scale can be used. In these models, the surface density f W/T is approximated by a force per unit length φ W/T , which is defined for each point C belonging to Γ and in contact with the workpiece.

For instance and for the sake of simplicity, Eq. ( 8) recalls a simple Kienzle model [START_REF] Kienzle | Spezifische schnittkrafte bei der metallbearbeitung[END_REF]. Like most models at this scale, it is based on the local depth of cut h, also called the uncut chip thickness, which is nothing else than the Boolean chip thickness near Σ. At point C, in a local frame linked to Σ and Γ , each component φ j of φ W/T can be written as

φ j = K j h h 0 a j , ( 8
)
where φ j is a force for the unit length of Γ and, by definition, h ≥ 0. K j and a j are parameters related to a given edge-material pair and h 0 is a reference length. Their values are obtained from experiments. Usually, this approach leads us to consider a discretization of Γ and Σ. In our work, we consider that the whole rake face is an assembly of rake faces of a set of elementary tools (Fig. 4). This allows us to take into account a great variety of cutting force models and tool geometries. Once the discretization of elementary tools is chosen, a model must be described to calculate F (k)

W/T and its point of application C (k)

F for each elementary tool T (k) . This will not be discussed here in detail, but φ W/T must be integrated along Γ (k) . At each point C (k) F (Fig. 4), in a local frame linked to Σ (k) and Γ (k) , each component F j of F W/T can be written as

F j = b 0 φ j dΓ, ( 9
)
where dΓ is the element of length along Γ (k) and b the length of Γ (k) . φ j can be obtained by any model. Usually, C

F is supposed to be on Γ (k) . Once such a discretization has been carried out for n F elementary T (k) , and

F (k) W/T are cal- culated, we have the forces F (k) T /W = -F (k)
W/T applied to the workpiece at the coincident points

P (k) F = C (k)
F . These points play the role of sampling points in a numerical integration of the rhs in Eq. ( 3). The equation can be written as

Q cW = n F k=1 H cT ( C(k) F ) • F (k) c W/T , (10)
Q sT = n F k=1 H sT ( P(k) F ) • F (k) s T /W . ( 11
)
Y cW and Y sT are obtained from Q cW and Q sT using Eq. ( 7), and since they are the main sources of excitation in Eq. ( 5), the quality of the simulation greatly relies on the accuracy of h, or other relevant geometric parameters of the Boolean chip and cutting conditions, and the accuracy of their evolution as a function of τ , everywhere on Σ C . The Boolean chip is of course very different to an actual chip, but the literature shows that it allows very satisfactory simulation at the macroscopic level for the two types of expected results: sufficiently accurate models of the cutting forces on the one hand (e.g., [START_REF] Jeong Hoon Ko | Development of a virtual machining system, part 1 approximation of the size effect for cutting force prediction[END_REF][START_REF] Yun | Development of a virtual machining system, part 3: cutting process simulation in transient cuts[END_REF][START_REF] Roth | Mechanistic modelling of the milling process using an adaptive depth buffer[END_REF][START_REF] Roth | Mechanistic modelling of 5-axis milling using an adaptive and local depth buffer[END_REF][START_REF] Jeong | Time domain model of plunge milling operation[END_REF][START_REF] Surmann | Simulation of milling tool vibration trajectories along changing engagement conditions[END_REF][START_REF] Merdol | Virtual cutting and optimization of three-axis milling processes[END_REF][START_REF] Damir | Force prediction and stability analysis of plunge milling of systems with rigid and flexible workpiece[END_REF][START_REF] Wan | Phase width analysis of cutting forces considering bottom edge cutting and cutter runout calibration in flat end milling of titanium alloy[END_REF][START_REF] Zhang | Process modeling and toolpath optimization for five-axis ball-end milling based on tool motion analysis[END_REF][START_REF] Murat | Stability of peripheral milling operations with long end mills[END_REF]) and of the machined surface S M (t E ) on the other hand (e.g., [START_REF] Yun | Development of a virtual machining system, part 2 prediction and analysis of a machined surface error[END_REF][START_REF] Surmann | Simulation of milling tool vibration trajectories along changing engagement conditions[END_REF][START_REF] Weinert | Modeling regenerative workpiece vibrations in five-axis milling[END_REF][START_REF] Kersting | Simulation concept for predicting workpiece vibrations in five-axis milling[END_REF][START_REF] Biermann | A general approach to simulating workpiece vibrations during five-axis milling of turbine blades[END_REF][START_REF] Quinsat | Characterization of 3d surface topography in 5-axis milling[END_REF][START_REF] Desai | On cutter deflection surface errors in peripheral milling[END_REF]).

Step-by-step numerical solution

A simplified flowchart of the incremental numerical process that is used to solve Eq. ( 5) is given in the lower part of Fig. 3. Among the data used to establish Eq. ( 5), the NC definition of the tool/workpiece relative motion allows us to obtain at any time τ the positions of the frames R c and R s in R g . This relative motion is described by means of a known matrix Rc s (τ ). This relative motion is responsible for the occurrence of machining, and thus gives rise to Q cW and Q sT in Eq. [START_REF] Won | Complete swept volume generation, part i: Swept volume of a piecewise c1-continuous cutter at five-axis milling via gauss map[END_REF], which transform into Y cW and Y sT in Eq. ( 5), the system we solve. This indicates that both Eq. (1) and Eq. ( 5) are non-linear systems.

Eq. ( 5) is solved numerically from the time t 0 to the final time t E using a step-by-step process that gives the solution for t

∈ {t 1 , • • • , t i , • • • , t E },
where, for instance, t i = t i-1 +∆t if a constant time step ∆t is used. The numerical scheme [START_REF] Bathe | Finite Element Procedures[END_REF] may be implicit (e.g., the Newmark scheme), which implies that iterations are carried out to obtain the solution at the end of each time step t i , or explicit (e.g., the central difference).

In any case, the process assumes that a solution is known at time t = t i and gives the new solution at time t i+1 = t i + ∆t. In the method that we proposed ( [START_REF] Lorong | Contribution of computational mechanics in numerical simulation of machining and blanking: state of the art[END_REF][START_REF] Lorong | Simulation of a finishing operation: Milling of a turbine blade and influence of damping[END_REF]), it is implied that we build an accurate model of the Boolean chip Ω B (t i , t i + ∆t) for each time step, and for each iteration if needed, to compute Y cW and Y T s on the rhs of Eq. ( 5). As shown in the following sections, this is not straightforward when dealing with non-rigid workpieces.

We introduce t I and [t I , t 0 ], which is an artificial interval of time that describes the change between the initial configuration of the workpiece I C and its clamped configuration 0 C at the beginning of the machining operation. This change is a purely static process, which means that during [t I , t 0 ], the time τ is simply a means with which to describe continuously the change from I C to 0 C (kinematic time).

Modeling of material removal

This section presents the main aspects of material removal modeling in the frame of simulations at the macroscopic scale. We consider that y s and y c are known step by step, by means of the numerical scheme, on any [t i , t i+1 ] ⊂ [t 0 , t E ]. See for instance [START_REF] Lorong | Simulation of a finishing operation: Milling of a turbine blade and influence of damping[END_REF] for details. Practical implementations are detailed in the next section.

As explained in the previous section, the tool/workpiece relative motion must be accurately described. This is a requirement first to obtain an accurate coupling by means of Y cW and Y cW , and second because this is the only way to obtain a reliable model of the machined surface S M (τ ), whose final state S M (t E ) is usually one of the main expected results.

In other words, we must be able to calculate where any point C of the tool (particularly of Σ) is, with respect to any point P of the workpiece, at any time τ ∈ [t I , t E ]. If these two points are coincident at time τ , they belong to the Boolean chip and matter in P is erased.

To describe the tool/workpiece relative motion and define the Boolean chip geometry, usual continuum mechanics concepts and Lagrangian formulation are used: a material frame R m and its associated material coordinates have to be used as far as the material removal of a non-rigid workpiece is concerned. We have to pay attention to different configurations of the workpiece that must be clearly specified. I C is the initial, undeformed, reference configuration. 0 C is the configuration once the workpiece is clamped on its supports. τ C is the deformed and machined configuration at time τ . In these configurations, the domains occupied by the workpiece material are denoted by Ω W (t I ), Ω W (t 0 ) and Ω W (τ ). The boundaries of these domains, which model the surface of the workpiece, are denoted ∂Ω W (t I ), ∂Ω W (t 0 ) and ∂Ω W (τ ).

Relative motion of R c with respect to R s

To consider the relative motion of the tool with respect to the workpiece, calculations are made using the coordinates of reference C of a point C in R c , and then obtaining its deformed position

C c (τ ), also in R c . The transformation of coordinates from R c to R s then gives C s (τ ) in R s .
The relative motion of C with respect to the points P of the workpiece can then be defined. This is the main object of Section 3.

The relative motion of R c with respect to R s is defined per the NC program and by the rotation of a spindle on which either the tool or workpiece is linked.

Using homogeneous coordinates X, the transformation of coordinates from Xg in R g to Xs in R s can be written

Xs = Rg c • Xg and Xg = Rc g • Xs = Rg c -1 • Xs , ( 12
)
where Rg c is a 4 × 4 matrix that describes both the translation and large rotation of R s with respect to R g and

Xg ∆ = [X T g 1] T (13)
is a 4 × 1 column matrix. The change of coordinates from R c to R s is simply

Xs = Rc s • Xc where Rc s = Rg s • Rc g . ( 14
)
The assumption that the motions of R c and R s in R g are known and described by the NC program implies that Rc s (τ ), Rg s (τ ) and Rc g (τ ) are known for any τ ∈ [t I , t E ].

Kinematic models of the active parts of the tool

As previously mentioned, to treat a wide variety of tools, the set Σ of all points corresponding to the active areas of the tool is partitioned into elementary rake faces Σ (k) associated with elementary tools T (k) . Each insert or each tooth is modeled by one or several T (k) and we have

Σ = k Σ (k) (15)
Fig. 4 and Fig. 5 show examples of T (k) and their Σ (k) . An example of a rake face To calculate the material removal, the relative motion of each Σ (k) with respect to the workpiece matter must be found. The first step is to obtain C s (τ ), the coordinates of any point C of the tool, at time τ , in R s .

Σ (k) = C 1 C 2 C 3 C 4 is shown in
Σ (k 1 ) (t i ) m Σ (k 1 ) (t i+1 ) m Σ (k 4 ) (t i ) m Γ (k 2 ) (t i+1 ) m Γ (k 2 ) (t i ) (d)
The column C c (τ ) that comprises the coordinates of a point C of the tool, at instant τ , as a function of its coordinates C ∆ = C c (t I ) in the initial configuration I C of the WTMS is a consequence of Eq. (2b) and Eq. (4b):

C c (τ ) = C + H c ( C) • Ψ c • y c (τ ). ( 16
)
The position of any point C of the tool, in R s , at any time is obtained from its reference position C using the procedure c-to-s described in Proc. 1.

Procedure c-to-s: for any required τ ∈ [τ 1 , τ 2 ] and for any C, get C s . Steps: 1) get C c from Eq. ( 16), 2) set up homogeneous coordinates Ĉc from C c and Eq. ( 13), 3) use Eq. ( 14) to get Ĉs from Ĉc , 4) extract C s from homogeneous coordinates Ĉs using Eq. ( 13).

Proc. 1: Pseudo-code of procedure c-to-s to get the position of a point of the tool in R s from its reference coordinates in R c .

This procedure c-to-s must be applied to any point of Σ but it is not sufficient to compute the Boolean chip. This is because the workpiece undergoes deformations in R s . Thus, the definition of m Ω S (τ 1 , τ 2 ), the domain swept in the workpiece by Σ, for τ ∈ [τ 1 , τ 2 ] requires additional effort. The following sections describe the method that we propose.

y s x s O s Ω W (t I ) P(t I ) Σ target (a) Initial configuration I C at t I < t 0 y s O s x s Ω W (t 0 ) P(t 0 ) Σ target (b) Clamped configuration 0 C at t 0

Material frame and the kinematic model of the workpiece

To illustrate the geometric aspects that we consider, a model workpiece, the simple geometry of which is depicted in Fig. 6 and Fig. 7, was chosen to propose 2D illustrations of the concepts by means of cross sections. The position of any point P of the workpiece is obtained at any instant τ in R s from its reference position P and displacement u s , which are given by Eq. (2b) and Eq. (4b):

P s (τ ) = P + H s ( P) • Ψ s • y s (τ ) (17)
Eq. ( 17) defines a one-to-one mapping Φ that depends on the time τ by means of y s . Φ |τ describes the transformation of the initial domain Ω W (t I ) occupied by the workpiece in I C to the current non-machined one Ω W (τ ) at instant τ . Note that Ω W (τ ) is a virtual domain, because it is larger than the actual domain Ω W (τ ) if material has been removed from the workpiece:

Ω W (τ ) ⊆ Ω W (τ ). Φ |τ :    Ω W (t I ) → Ω W (τ ), P → P s (τ ), P s (τ ) = P + H s ( P) • Ψ s • y s (τ ) . ( 18
)
For any time τ ∈ [t I , t E ] the one-to-one mapping Φ |τ defines a curvilinear coordinate system (α m , β m , γ m ) |τ on the domain Ω W (τ ). In this coordinate system (Fig. 6, Fig. 9) that follows the workpiece deformation, a point P of the workpiece has constant coordinates α P , β P and γ P . This leads to the notion of the material frame.

We denote by R m the material frame defined by the one-to-one mapping Φ and its associated coordinate system (α m , β m , γ m ). We denote by Rm the initial configuration of R m :

Rm ∆ = R m (t I ). We thus write (α m , β m , γ m ) |τ = Φ |τ (α m , β m , γ m ) |t I (19a) R m (τ ) ∆ = Φ |τ Rm . ( 19b
)
There is interest in considering R m only when the relative motion of the tool with respect to the workpiece is considered. Only the relative motion of the tool with respect to R m has meaning in the modeling of material removal.

Fig. 6 illustrates the link of the material coordinate system to the matter; there is a transform from the initial configuration (x m , y m , z m ) = (α m , β m , γ m ) |t I to the configuration (α m , β m , γ m ) |t 0 once the workpiece is clamped on its supports.

y m x m z m (a) Initial workpiece: m Ω W (t I ) y m x m z m plane Π (b) m Ω W (t 3 )
Figure 7: 3D views in Rm , at t I and t 3 , of the upper part of the workpiece shown in Fig. 6 by means of its section in the plane Π.

In the initial configuration I C, it is simpler to start with a Cartesian coordinate system. To emphasize this choice, we denote position by (x m , y m , z m ) in this particular initial instance (α m , β m , γ m ) |t I . Although not necessary, to simplify all expressions, it is natural to choose this initial coordinate system such that (x m , y m , z m ) = (x s , y s , z s ). Fig. 6 and Fig. 9 illustrate this choice. This leads to P = [α P β P γ P ] T .

Provided that the domain affected by the transformation Φ is extended, R m defines a material frame in which it is possible to define and observe simply the material removal and the relative motions between the rake faces with respect to the matter. This is described in the next section.

To simplify the problem, we denote by m Ξ(τ ), with a left exponent m, the image in the reference configuration I C, by means of the one-to-one mapping Φ -1 |τ of a geometric entity Ξ, which is in motion and the domain of which is Ξ(τ

) at time τ in R s . m Ξ W (τ ) ∆ = Φ -1 |τ Ξ W (τ ) . ( 20
)
At time t I , we thus have

m Ξ W (t I ) = Ξ W (t I ). (21)
For instance, m Ω W (τ ) is the image, for the initial configuration, of the domain Ω W (τ ) filled by the workpiece at time τ in R s :

m Ω W (τ ) ∆ = Φ -1 |τ Ω W (τ ) , m Ω W (t I ) = Ω W (t I ).
Note that for a point C in motion with respect to R s , we could use the notation defined by Eq. ( 20). However, for the sake of simplicity, we prefer to use the notation C ≡ m C for a point and C ≡ m C for its coordinates, which leads us to write C(τ )

∆ = Φ -1 |τ C(τ ) and C(τ ) ∆ = Φ -1
|τ C s (τ ) . This notation is coherent with that chosen for a point P of the workpiece in Eq. ( 17) and for R m in Eq. ( 19), but C(τ ) is a function of time. In contrast, note that neither P nor Rm is a function of time.

m Σ (k) (t i+1 ) m Γ (k) (t i+1 ) C 3 (t i ) C 2 (t i ) C 4 (t i ) C 1 (t i ) C 3 (t i+1 ) C 2 (t i+1 ) C 4 (t i+1 ) C 1 (t i+1 )
(a) Domain swept by the image in Rm of the rake face

Σ (k) = C1 C2 C3 C4 of the elementary tool T (k) (3D-view) C 3 (t i ) C 2 (t i ) C 2 (t i+1 ) C 3 (t i+1 ) workpiece boolean chip (b) Domain swept in Rm by the image C2 (τ ) C3 (τ ), for τ ∈ [t i , t i+1 ] (2D pro- jection)
Figure 8: Image of the rake face Σ (k) of the elementary tool T (k) , and its swept domain in Rm .

Until now, R m has only been defined for the points belonging to the workpiece. This is not sufficient because T (k) are not permanently inside the workpiece domain Ω W . In milling, for instance, the elementary tools are always outside the material or partly outside. Fig. 8 depicts an example of such a situation for a subset C 1 C 2 C 3 C 4 of an elementary rake face that is described by its four vertices. At time t i , all the points are outside the material. At time t i+1 , C 1 and C 2 are inside, as seen in Fig. 8b. Fig. 10 and Fig. 11 illustrate this situation in Rm for most of the rake faces.

It is essential to detect the time intervals during which contact occurs and removes material. Positions of the points in Σ must be carefully followed, even when they are not in contact with the workpiece. This is particularly true for the vertices used to describe the geometry of any Σ (k) because their positions with respect to the workpiece are required for several calculations. It is thus necessary to build a one-to-one mapping on a domain larger than Ω W (t I ).

To extend the domain relating to Φ, we define a domain m B such that Ω W (t I ) ⊂ m B. The one-to-one mapping Φ + is built for any point P in B to enrich Eq. ( 18), using an additional function A P, τ such that, for all τ ,

P (t I ) ∈ Ω W (t I ) ∪ ∂Ω W (t I ) ⇒ A P, τ = 0, (22)
Φ + |τ :    m B → B(τ ), P → P s (τ ) P s (τ ) = Φ |τ + A P, τ . ( 23
)
β m α m m Ω W (t I ) = Ω W (t I ) Ω W (t) τ = t I Φ + (t) α m β m τ = t y m = y s x m = x s O m =O s y m = y s x m = x s O m =O s ω m ω m Φ(t) m B = B(t I ) = B B(t) m Ω W (t) P C(t) Φ(t)

C(t) P(t)

Φ + -1 (t) We keep the notation Eq. ( 20) for the entities transformed by (Φ + |τ ) -1 . A construction of Φ + Eq. ( 23) that respects Eq. ( 22) is easy to achieve making an FE approximation. This approach ensures that Φ and Φ + are identical to each other for any point of the workpiece. Note that when making an FE approximation, there is still a wide choice of possible extensions of Φ into Φ + . This gives a wide range of possible A, the only requirement being that Φ + must be a one-to-one mapping. An example of practical implementation is proposed Section 4. Fig. 9, Fig. 10 and Fig. 13 illustrate Φ, Φ + and B in a plane Π orthogonal to z s . In these figures, the lines iso-α m and iso-β m , which are the images of x s = i.∆l and y s = j.∆l by Φ + are shown inside B.

Definition of the Boolean chip and material removal

By means of Eq. ( 25), we define m Ω B (τ 1 , τ 2 ), which we call the Boolean chip. The Boolean chip is generated during the time interval [τ 1 , τ 2 ]. m Ω B is illustrated in Fig. 1, Fig. 8 We consider m Ω S (τ 1 , τ 2 ), the 3D domain swept in Rm by the 2D domain m Σ between times τ 1 and τ 2 :

m Ω S (τ 1 , τ 2 ) ∆ = τ ∈ [τ 1 ,τ 2 ] m Σ(τ ) . ( 24
)
Fig. 1, Fig. 5, Fig. 10, Fig. 11 and Fig. 12 give several examples of domains swept by rake faces in Rm . The plane Π of these figures is depicted in Fig. 7.

By definition, the Boolean chip is the domain m Ω B (τ 1 , τ 2 ) erased inside the workpiece by the motion of m Σ between τ 1 and τ 2 :

m Ω B (τ 1 , τ 2 ) ∆ = m Ω W (τ 1 ) ∩ m ΩS (τ 1 , τ 2 ), (25)
where m Ω W (τ 1 ) is the open domain occupied by the workpiece material at time τ

1 in Rm . At τ 2 , this domain is m Ω W (τ 2 ) = m Ω W (τ 1 ) -m Ω B (τ 1 , τ 2 )
, where for any couple (τ 1 , τ 2 ) such that

τ 2 > τ 1 > t 0 , we have m Ω W (τ 2 ) ⊆ m Ω W (τ 1 ) ⊆ m Ω W (t 0 ), (26)
and

m Ω W (τ 2 ) = m Ω W (t 0 ) -m Ω B (t 0 , τ 2 ). (27)
The contact zone Σ C between the tool and workpiece at time τ is easily obtained because, as a consequence of the definition of the Boolean chip, it is the common boundary between the workpiece and tool at τ : m Ω S tooth 1 (t 0 , t 5 )

m Σ C (τ ) = m Σ(τ ) ∩ m ∂Ω W (τ ). (28
m Ω S tooth 2 (t 0 , t 5 )

m Ω S tooth 3 (t 0 , t 5 )

Figure 11: Cross section in plane Π of domains m Ω S tooth i (t 0 , t 5 ) swept in Rm by the rake faces of teeth 1, 2 and 3.

m Ω W (t 5 )

m Ω W (t 4 ) m Ω W (t 3 ) m Ω W (t 2 ) m Ω W (t 0 ) m Ω S tooth 1 (t 1 , t 2 )
m Ω S tooth 1 (t 0 , t 1 )

m Ω S tooth 2 (t 2 , t 3 )

m Ω S tooth 2 (t 3 , t 4 )

m Ω S tooth 2 (t 4 , t 5 ) { part of m ∂Ω W (t E ) { part of m ∂Ω W (t E ) part of m ∂Ω W (t E ) m Ω W (t 1 ) { m Ω B (t 0 , t 1 ) m Ω B (t 1 , t 2 ) m Ω B (t 2 , t 3 ) m Ω B (t 3 , t 4 ) m Ω B (t 4 , t 5 ) { part of m ∂Ω W (t E ) { part of m ∂Ω W (t E )
Ω W (t 8 ) Ω W (t 8
) The position C of any contacting point C of Σ(τ ) (i.e., such that C(τ ) ∈ m Σ C (τ )) can then be obtained from its image C(τ ) in Rm . Let us recall that this is necessary to calculate Y cW and Y sT on the rhs of Eq. ( 5).

Ω W (t 2 ) Ω W (t 2 ) Σ(t 2 ) Σ(t 2 ) Σ(t 8 ) Σ(t 8 ) Φ R (t 8 ) Φ + (t 8 ) Φ R (t 2 ) Φ + (t 2 ) Φ + (t 2 ) Φ R (t 2 ) Φ R (t 8 ) Φ + (t 8 ) α m β m t = t I rigid workpiece: Φ R non-rigid workpiece: Φ + t = t 2 t = t 8
To build m Ω B , we must be able to calculate at least for each point C of Σ and for any time τ its image C(τ ) in Rm . This is always possible once Φ has been extended to a domain B that is large enough. As a matter of fact, it is only necessary that B contains the whole Boolean chip and the corresponding swept domains that intersect the workpiece during the machining operation. The points C of the tool that are not in this last defined B are of no utility.

Using the procedure c-to-m given in Proc. 2, the coordinates in R s are obtained from the coordinates of C in R c . We then obtain the coordinates C of the image m C of C in Rm by means of (Φ + |τ ) -1 :

C(τ ) = (Φ + |τ ) -1 C s (τ ) . ( 29
)
Eq. ( 29) is used to calculate m Ω B (τ 1 , τ 2 ) for any τ 1 ≤ τ 2 . All Boolean operations can then be done in Rm . On a given time interval [τ 1 , τ 2 ], the relative motion between (R s ) and (R c ) is known; thus, Rc s (τ ) is known in this interval. As soon as y s (τ ) and y c (τ ) are known on the whole interval, the relative motion of the tool with respect to the workpiece is calculated using the procedure c-to-m described in Proc. When needed, for instance to plot the trajectory of any point C in R m (t) like in Fig. 10, Fig. 11, Fig. 12 or Fig. 13, it is also possible, for any t ≥ τ , to obtain at t the image of C(τ ) in R m by means of its coordinates

Φ + |t C(τ ) . Φ +
|t is also useful because it permits us to calculate the position, on the tool, of the integration points C (k) F used in Eq. ( 10) to obtain Y cW . This implies the calculation, for each

T (k) involved, C (k) F c (τ ) in R c from C (k) F s (τ ) in R s .
Note that the coordinates of these points are functions of time because the shape of the section of the Boolean chip is a function of τ . This leads to the procedure m-to-c, which is depicted in Proc. 3.

At this stage, it is important to note that in the case that it is possible to assume rigid-body motion for the workpiece, which is the usual assumption made in most other works, procedures c-to-m and m-to-c become straightforward because Φ, and thus Φ + , becomes a simple isometry Φ R that is equivalent to a change of Cartesian coordinates. This is illustrated on the left of Fig. [START_REF] Assouline | Numerical simulation of machining at the macroscopic scale: dynamic models of the workpiece[END_REF].

In addition to the comparison of the two families of one-to-one mappings, with or without the rigid-body assumption (at τ = t 2 and τ = t 8 ), Fig. 13 also shows the position of some rake faces and the trajectories of their cutting edges with respect to the material at τ = t 2 and τ = t 8 in coherence with Fig. 10 and Fig. 11. Note than in R m , and in the case of the non-rigid Procedure m-to-c: for one point C of the tool, for any required τ ∈ [τ 1 , τ 2 ] (mainly for points C ∈ Σ C required to compute Y cW ), get C from C(τ ), y c (τ ), y s (τ ) and Rs c (τ ) = Rc s (τ ) -1 . Steps: 1) get C s (τ ) from its coordinates C(τ ) using Φ + |τ , 2) build homogeneous coordinates Ĉs (τ ) from C s (τ ) using Eq. ( 13), 3) use Eq. ( 14) to get Ĉc (τ ) from Ĉs (τ ), 4) extract C c (τ ) from homogeneous coordinates Ĉc (τ ) using Eq. ( 13), 5) get C from C c (τ ) by inversion of Eq. ( 16).

Proc. 3: Pseudo-code of procedure m-to-c, used to obtain the position of reference C of a point of the tool in R c from the coordinates of its image C(τ ) in Rm . workpiece, the shape does not remain the same as time passes. This is not the case for the rigid workpiece. It must be noted in Fig. 13 that the positions of the parts of Σ on the plotted trajectories in R m are the positions of the rake faces in R s at time t i only for Σ(t i ). This is why the definition of the Boolean chip makes sense only in R m .

Concluding remarks

We note the following.

1. The transformation of the domain Ω W of the workpiece as material removal develops can only be observed in a material frame such as R m where all the points of the workpiece have constant coordinates. Among these material frames, Rm simplifies the calculation because the points of the workpiece do not move with respect to each other or because they have disappeared as a consequence of material removal.

The domain of the workpiece at time

τ in Rm is m Ω W (τ ) and m Ω W (τ ) ⊆ m Ω W (t I ).
Before machining, at time τ = t 0 , the workpiece domain is m Ω W (t 0 ) = m Ω W (t I ). After machining, at time t E , the domain reduces

to m Ω W (t E ), with m Ω W (t E ) ⊂ m Ω W (t 0 ).
This change of m Ω W and Ω W is depicted in Fig. 7, Fig. 12 and Fig. 13 for different values of

t i ∈ [t 0 , t E ].
3. Before machining, the "workpiece surface" is the boundary m ∂Ω W (t 0 ) of m Ω W (t 0 ). After machining, the "workpiece surface" is the boundary

m ∂Ω W (t E ) of m Ω W (t E ).
The "machined surface" S M (t E ), which is one of the simulation aims, corresponds to the part of the boundary m ∂Ω W (t E ) that does not belong to the initial surface m ∂Ω W (t 0 ) of the workpiece.

Material removal modeling is an important aspect of the approach because it defines

m Ω W (t E ), and consequently the geometry of the machined surface at the end of machining. More generally, simulation must allow a fine description of the geometry of m Ω W at the end of each simulation step. Indeed, this is necessary because the calculation of the fluctuations of the cutting forces depends directly on the quality of the description of m Ω W between two successive actions of the tool. This fluctuation induces vibrations that, in turn, modify the swept domain m Ω S and thus modify the fluctuations of cutting forces. It is at this level that there is coupling between the dynamic behaviors of the elementary tools and the workpiece, by means of Q cW and Q sT that appear on the rhs of Eq. (1b) andEq. (1a).

5. Note that it is not easy to accurately prescribe a tool trajectory to obtain a given expected machined surface. For instance, if the aim is to obtain a machined surface S M = Σ target parallel to a plane y s constant, the envelop of the cutting edges of the tool must follow the geometry of Σ target . Such a surface is illustrated by its trace in Fig. 6a and Fig. 6b. The trajectory of the tool must take into account the change in Σ target , as depicted in Fig. 6b, as soon as the workpiece undergoes significant deformation. One source of such deformation is clamping, but it is also encountered in thin-walled pieces. In Fig. 6, Fig. 10, Fig. 11, Fig. 12 and Fig. 13, trajectories of cutting edges that satisfy the requirement of generating Σ target are observed.

6. Finally, as depicted in Fig. 12, the surface is generated progressively by several successive actions of the cutting edges. These edges usually belong to different elementary tools (part of an insert or part of a tooth in milling), which does not simplify the modeling.

The figure shows that a very small part of m ∂Ω W (t E ) is generated between t 0 and t 5 . This reinforces the need to have an accurate representation of each intermediate state of

m ∂Ω W .

Practical implementation of the theoretical approach

To solve Eq. ( 5), m Ω B (t I , t) and m ∂Ω W (t) must be described with a good level of accuracy owing to their close link with the variation in cutting forces. In practice, this leads to a huge calculation time step after time step, and appropriate modeling must be carried out to obtain the desired accuracy in a sufficiently short computation time. To focus on such modeling, we assume in this section that each kinematic quantity is known on [t I , t i ] and also on [t i , t i+1 ] by means of an appropriate solving process. Rc s , y c and y s are thus known quantities on these time intervals.

The theoretical approach used to define and build the Boolean chip has been described in Section 3. Some approximations have to be made for a practical implementation of the method and simulations.

In this section, we first propose an efficient description of Φ + and the box B that we introduced to define the extension of the mapping Φ on any time interval [t I , t]. We then explain how we set up the model that we use to describe m Ω S (t i , t i+1 ), the domain swept by Σ during [t i , t i+1 ]. In this context, Ω W is described by means of a model using dexels ( [START_REF] Van Hook | Real-time shaded nc milling display[END_REF]). This allows us to describe easily and efficiently the Boolean chip m Ω B (t i , t i+1 ) and the evolution of m Ω W from t i to t i+1 .

Finite element model of R m

The kinematics of the material domain R m is approximated by an FE description. For the sake of efficiency, the mesh is restricted to a portion of the workpiece and its outside that are sufficient to describe the interaction between the tool and the workpiece; it is not required that m B completely contain m Ω W .

A mesh M + m is thus built and used. It must include the whole portion of the workpiece that is subjected to machining, and a sufficient subset of the outside of the workpiece. To obtain an explicit one-to-one mapping in each element, and have the smallest number of numerical operations, four-node tetrahedrons elements (T4) are chosen. This is done whatever type of element is used in M s to model the workpiece dynamics and obtain systems Eq. ( 1) and Eq. ( 5).

A simplified flowchart is given on the right of Fig. 3. An example of a choice of m B is given in Fig. 15c for the workpiece whose partial mesh M s is depicted in Fig. 14. Fig. 15a shows the trajectory of the cutter center that has been taken into account to choose m B and the part of the mesh M s (T10; isoparametric 10-node elements) that was transformed into the mesh M m (T4) is shown in Fig. 15d.

Introducing the matrix H + of interpolation functions, and q + m , the column that contains the nodal displacements (DOF) of M + m , we make the approximation

P s (τ ) = P + H + ( P) • q + m (τ ). ( 30
)
If M s contains types of elements other than T4, or if the meshes are not the same, there is a slight additional approximation to the initial one given by Eq. ( 23). This can be controlled by the size of the elements of M + m . As a matter of fact, this approximation is not of great importance because we only use an approximation of displacements, not of their gradients.

Because of its purpose, the mesh M + m is partitioned into two sets of T4 that respect the boundary of the workpiece, as shown in Fig. 15b. We call M m the mesh of the workpiece and M m the mesh of its outside. The mesh M m is constructed with the constraint that first it fills the domain inside m B, but outside M m , and second the faces of the generated T4 are the same as those of M m on their common boundary.

The nodes of M m are fastened to the displacements of the workpiece given by y s using Eq. (2b) for each node and Eq. ( 4). This leads to the following partition of q + m into q m and q m , where q m comprises the displacements of all the nodes of M m . This results in a partition of H + into H m and H m , and Eq. (30) becomes

P s (τ ) = P + H m ( P) H m ( P) • q m (τ ) q m (τ ) . ( 31
)
q m contains all the nodal displacements of M m that are not associated with the nodes belonging to m ∂Ω W (t I ). It is seen that this corresponds to the choice

A P, τ = H m ( P) • q m (τ ), (32)
and a sufficient condition that Eq. ( 32) satisfies the requirements of Eq. ( 22) is given by

P (t I ) ∈ Ω W (t I ) ∪ ∂Ω W (t I ) ⇒ H m ( P) = 0. (33)
As stated before, the only requirement on q m is that it should lead to a one-to-one mapping.

Once the displacements of the nodes of M m are obtained by means of Eq. ( 17), the constant matrix G s m such that q m = G s m • q s is known and we obtain the relationship Eq. ( 34) between q m and y s :

q m = Ψ m • y s with Ψ m ∆ = G s m • Ψ s . ( 34
)
Ψ m is thus independent of time and it can be set up once and for all before starting the simulations.

To define the mapping Φ + outside the workpiece domain, we have to decide which procedure to choose. Eq. ( 34) enforces the displacements of the nodes that belong to the boundary of the workpiece because the displacements are a subset of those in q m . The displacement of the other nodes of M m (i.e., q m ) must be defined.

Among possible methods, q m can be obtained by an FE calculation using the model described by M m and an arbitrary material (e.g., the workpiece material). This allows us to obtain a matrix Ψ m with the same number of columns as Ψ m by means of a simple process. Each column of Ψ m is obtained by solving a static problem where the (small) displacements of the nodes belonging to m ∂Ω W are equal to their corresponding value in Ψ m . Fig. 16 illustrates the deformation of B using displacements prescribed by Ψ 2|m , the second column of Ψ m , which describes the mode shape of the second natural mode in the machined region. Outside the workpiece, the deformation is given by Ψ 2|m , the second column of Ψ m , which extends the mode shape to the whole domain B. The mode shape is only shown inside B.

Ψ m gathers deformed shapes associated with unit values of each component of y s and q m = Ψ m • y s with that choice. As far as the workpiece deformation remains small, this ensures a one-to-one mapping.

Once Ψ m is known, Ψ m is calculated once and for all before the simulation. All the components of q + m are linked to y s and we have

q m q m = Ψ m Ψ m • y s . ( 35
)
Eq. ( 31) can be written P s (τ ) = P + H+ ( P) • y s ,

where

H+ ( P) = H m ( P) • Ψ m H m ( P) • Ψ m .

FE model of the rake faces and BREP description of the swept domain

This section describes the method proposed to model the swept domain in the interval

[t i , t i+1 ].
The aim is to model the 3D domain m Ω S (t i , t i+1 ) that is swept by the 2D domain m Σ in Rm .

As stated in Section 3.2, to be able to model a wide range of tools, a first partition of the active parts of the tool in T (k) is carried out. This leads to a partition of Σ into n Σ elementary rake faces Σ (k) . Examples of such first decomposition are depicted in Fig. 4 and Fig. 5. This decomposition of Σ by means of a mesh M Σ is then described using three-node plane triangles T3, which describe each Σ (k) with a mesh M (k) Σ . Examples of mesh patterns to model Σ (k) and its cutting edge Γ (k) are shown in Fig. 17a.

T3 elements are supposed to be small enough that their images obtained using Φ + or its inverse can also be approximated by T3 without significant loss of accuracy. This is not a strong requirement because we assume small displacements of the tool in R c and of the workpiece in R s , and thus, Φ + is close to an isometry.

The 3D domain m Ω S = m Ω S + m ∂Ω S is modeled by the union of m Ω (k)

S swept by m Σ (k) in Rm : m Ω S = k m Ω (k) S . ( 36
) Each swept domain m Ω (k) S is modeled using a BREP model M (k) S ; i.e., M (k)
S is a model of

m ∂Ω (k) S . M (k)
S is a set of T3. Faces at t i and t i+1 are given by the two positions of M

Σ at these time steps as shown in Fig. 17b.

As a consequence, m Ω S is indirectly modeled by a mesh M that is the union of all M (k)

S : M = k M (k) S . ( 37
)
Note that M is not exactly a BREP model of m Ω S as it also contains the meshes of internal boundaries in addition to its description of m ∂Ω S . Moreover, provided that a convenient sequential process is used to make the Boolean operations, the domains Σ (k) may overlap, mainly near small radii of the cutting edge Γ , thus leading to the overlapping of M (k) S . However, Γ (k) never overlap in our approach.

The topology of M

S is chosen, once and for all. An example of the mesh of side parts that are generated by the motion of the edges C 3 C 4 and C 4 C 1 of the boundary of Σ (k) is given in Fig. 17b. In this case, M S in any frame R provided that the coordinates of the nodes are given and lead to a non-degenerated representation (i.e., the Jacobian of the transformation from the position of the reference where the topology is described to that in R is positive anywhere). One or more intermediate positions must be used, for instance, when the time step ∆t = t i+1 -t i is too large to accurately follow large rotations, which would introduce large chord error.

Γ (k) C 1 C 2 C 3 C 4 C 1 C 5 C 2 C 3 C 6 C 4 C 1 C 4 C 5 C 6 C 7 C 2 C 3 (a) Example of types of mesh patterns M (k)
Σ for Σ (k) and associated location of Γ (k) . positions of the rake face Σ (k) of T (k)

Γ (k) cutting edge of T (k) 0.
Σ (k) = C 1 C 2 C 3 C 4 .
Figure 17: Examples of generic mesh patterns that approximate the elementary rake faces Σ (k) and related domains Ω (k)

S (t 1 , t 2 ) swept on [t 1 , t 2 ].
To complete the definition of M Σ and M (k) Σ , we denote by CΣ the column that gathers the coordinates in R c of the n N nodes of M Σ in the reference configuration of the tool:

CΣ ∆ = [ CT 1 CT 2 • • • CT n N ] T . ( 38
)
As in the case of the workpiece, M Σ is completely different to M c , but its nodes are fastened to the tool and follow its motion, as modeled by Eq. ( 16). The procedure c-to-m can thus be applied to each node C p : its coordinates Cp (τ ) in Rm are obtained from Cp . This gives CΣ (τ ), which comprises the images of all the nodes at τ in CΣ . We denote by C(k) Σ the subset of CΣ that only comprises the nodes describing the mesh M (k) Σ of Σ (k)

Geometric model of the workpiece obtained using dexels

In the proposed approach, the workpiece domain Ω W and its evolution are modeled using dexels ( [START_REF] Van Hook | Real-time shaded nc milling display[END_REF][START_REF] Won | Virtual workpiece: Workpiece representation for material removal process[END_REF]) that we briefly describe in this section. In practice, only the part of the workpiece domain inside the box B requires to be represented by dexels. The approximation of Ω W using dexels is interesting because the positions of the end points of a dexel in one direction are very accurate (they are defined by real numbers), which is important when modeling machined surfaces. Furthermore, in the context of Boolean operations with BREP domains that are described by T3 meshes-which is our choice for M-this leads to simple, efficient and robust algorithms. Robustness is essential for industrial problems.

Let R D denote a frame and its associated Cartesian coordinate system (u, v, w). First, a set of supports that are straight lines parallel to w is defined. The intersections of the supports with the plane (u, v) correspond to the centers of a regular grid of ∆u × ∆v rectangles. Dexels can then be seen as parallelepipedic domains of section ∆u × ∆v that are defined along the supports in the interval [w min , w max ]. The faces of dexels are parallel to the three orthogonal planes defined by (u, v, w).

A dexel is thus defined by (i, j, w min , w max ) where (i, j) are the integer coordinates of the support. The actual coordinates of the dexel are (u, v) = (i • ∆u, j • ∆v) in R D . A dexel can share the boundary or the support of other dexels, but in any case, different dexels must correspond to disjoint domains. This allows us to approximate any 3D domain using a set of dexels ( [START_REF] Van Hook | Real-time shaded nc milling display[END_REF]) associated with one or more directions. This is what we do to model Ω W (t i ) using a set that we denote by D(t i ) or simply D before . Usually, three orthogonal directions are used (triple dexels) ( [START_REF] Won | Virtual workpiece: Workpiece representation for material removal process[END_REF]). As it is shown in Fig. 2 and Fig. 3, it is the geometric model of the workpiece that is used to build the initial dexel model D(t I ) of m Ω W (t I ) (in practice, only the part of m Ω W (t I ) inside m B is built).

Even if the work is always carried out in Rm for efficiency, the dexels follow the workpiece deformation. Their supports are defined in Rm and are placed along isolines of α m , β m and γ m . Examples of their deformation are thus given by the shape of the isolines of α m and β m , inside the workpiece domain, in the different figures where they appear (Fig. 6, Fig. 9 and Fig. 13).

In the examples given in the last section, only one direction of dexels, orthogonal to the machined surface, was used. This is also the case in Fig. 18 where we illustrate the principle of the virtual material removal and its associated set of dexels.

The procedure to calculate the intersection between two domains, one being a BREP M modeled with T3 and the other being a set D of dexels, is based on the intersection, when it occurs, between each T3 and each dexel. The set of dexels belonging to D before , and that are involved in the intersection, is denoted D cut . The remaining dexels are denoted D uncut . D before = D uncut ∪ D cut and D uncut ∩ D cut = ∅ During the intersection process, which is depicted in Fig. 18, two other sets of dexels are considered: D kept and D erased .

D erased is the set of old or new dexels that have been found inside the BREP M. D kept comprises all new dexels that are not inside M and that contribute with D uncut to the description of D after . For the sake of simplicity, we denote by M ∩ D before the result of the intersection, and by D before the new set of dexels that describes the domain after removing M ∩ D before :

D after = D before -D cut + D kept , D erased = M ∩ D before , D after = D uncut ∪ D kept , D uncut ∩ D kept = ∅.
The Boolean chip m Ω B (t i , t i+1 ) is modelled using the set D erased .

Calculation of material removal and related cutting forces

The calculations of the Boolean chip and the related cutting forces are made elementary tool by elementary tool. B denote, for an elementary tool T (k) , the domain it sweeps and the Boolean chip it generates.

Let I τ ∆ = {t i , t i+1 }. I τ also contains all intermediate time steps, if there are any. For instance, in the case that there are two intermediate steps, τ ∈ I τ ∆ = {t i , 1 3 t i + 2 3 t i+1 ), 2 3 t i + 1 3 t i+1 , t i+1 }. These intermediate steps may be useful in improving the quality of the trajectories of the rake faces when ∆t = t i+1 -t i is too large for this aim (see Section 4.2). As a matter of fact, ∆t is mainly driven by the numerical scheme used to solve Eq. ( 5) and it is not always suitable to reduce it, because each time step requires a number of important numerical operations.

The procedure to obtain the Boolean chip is described in Proc. 4. For usual mechanistic models or for the Kienzle model Eq. ( 8), step 3.e allows us to calculate cutting forces from the geometry of D

(k) erased ≈ m Ω (k) B .
Once phase 3.d is complete for Σ (k) , all geometric and kinematic calculations required to feed the chosen model for the cutting forces can be made. This allows us to obtain

P(k) F = C(k) F and F (k) c W/T (Section 2.2).

C(k)

F is then calculated using the procedure m-to-c given in Proc. 3. Finally, Y cW and Y sT are obtained from Q cW and Q sT using Eq. ( 7). This last phase is completed quickly because the displacements of the nodes of M + m and M Σ are linked once and for all, and by means of contant matrices, to y s and y c , respectively. Moreover, only the nodes of elements containing the integration points P(k) F and C(k) F are affected by this operation. Because we assume that the transformation of a T3 is a T3, the procedure m-to-c is simplified: the position of

C(k) F in Σ (k)
(k) erased = M (k) S ∩ D (k-1) ; e) knowing D (k) erased , calculate C(k) F and F (k) c W/T ; f) D (k) := D (k-1) -D (k) cut + D (k) kept ; g) D (k) boolean = D (k-1) boolean + D (k) erased . D after = D (n Σ ) is the model of the new workpiece domain m Ω W (t i+1 ). D erased = D (n Σ ) boolean is the model of the Boolean chip m Ω B (t i , t i+1 ).
Proc. 4: Pseudo-code of the procedure to calculate the material removal

Examples of simulations of industrial workpieces

The examples presented are industrial examples, and we only intend to illustrate the ability of the proposed method to deal with complex non-rigid workpieces and complex tools such as a milling cutter. Kienzle models were used for the cutting force in both examples. The total time to complete the simulation, once the system Eq. ( 5) has been built is given for the two examples. Simulations were done on a laptop computer (MacBook Pro, OS X, Intel Core i7, 2.66 GHz, RAM 8 GB). Currently, algorithms have been coded with care but only use one core and do not take advantage of the GPU. An example of five-axis milling of a turbine blade was given in [START_REF] Lorong | Simulation of a finishing operation: Milling of a turbine blade and influence of damping[END_REF].

Example of a turning operation

In this first example, an axisymmetric workpiece (height of 1m, clamped on its lower face where the displacement is prescribed and equal to zero) was modeled to understand the evolution of machining stability during a turning operation on its upper part: facing of a circular ring (mean diameter 1.088 m, width 26.8 mm). Several constant rotational speeds were tested. The machined surface shown in Fig. 20 and Fig. 21 corresponds to a rotational speed leading to chatter during the operation. The radius of the circular rake face of the active part of the tool was 3 mm. Twenty-two elementary tools were disposed to model a sector of 100 deg of the insert to have an accurate description of the cutting edge.

The FE mesh M s is partially shown in Fig. 19. The upper part is modeled with 42,120 T10 that are then easily transformed into T4 to construct M m . The rest of the model is built using 32,040 Q8 shell elements. The two sets of elements are linked using rigid elements. The model has more than 765,463 DOF. The machined surface is orthogonal to the axis of symmetry and the machining goes from the outside to the axis of symmetry as depicted in Fig. 20 (radial feed = 0.35mm/revolution). The tool was assumed to be rigid and only 25 modes were kept to build Eq. (4a). The machined surface was described by only one set of 19,500,000 (390×50,000) dexels parallel to the axis of rotation. This corresponds to 14.57 dexels/mm (212.3 dexels/mm 2 ). A small part is depicted in Fig. 21 to show the degree of detail that can be obtained. For this surface, the maximum defect is small near the external diameter of the machined zone (z maxz min ≈ 5 µm) and important near its internal diameter where chatter occurs (z max -z min ≈ 300 µm) at the end of the operation. The total simulation was completed in 78 minutes and corresponds to 77 revolutions of the workpiece (∆t = 2.5×10 -4 seconds; 615,000 increments).

Note that the whole machined surface, which corresponds to a circular strip, was modeled. In this particular case, to keep the axisymmetry of the geometric model and avoid spurious effects, an additional one-to-one transformation Φ a from a cylindrical coordinate system (r m , θ m , z m ) to (x m , y m , z m ) and dexels are implemented in the z direction on a regular grid in a domain a B = [r min , r max ] × [0, 2π[×[z bottom , z top ]. In this case, m B = Φ a ( a B), and for all operations associated with intersections with dexels, Φ + must be replaced by the composition Φ + • Φ a .

Example of a face milling operation

In this second example, the influence of the workpiece deformation can be seen first because clamping generates a global defect, and second because the induced vibrations are not the same on the five parts of the machined surface, thus leading to different surface defects. Fig. 23 illustrates these two points.

The machined surface is orthogonal to the z axis (enclosing rectangle with dimensions of 375 mm × 76 mm, axial depth of cut 3 mm). The dynamics of the tool includes the fact that the spindle undergoes deformations, while the milling cutter is assumed to be perfectly rigid. To build the reduced model, only two bending modes modes were kept for the tool and 10 modes were used for the workpiece. The total simulation was completed in 74 minutes and corresponds to 245 revolutions of the spindle (∆t = 7.0 × 10 -5 seconds; 151,000 increments).

The maximum defect is z max -z min ≈ 150 µm and is mainly due to the static deformation induced by clamping (simulation of an error of location of 0.1 mm in the Z direction for one of the supports). Nevertheless, the simulation shows very different distributions of the defects, depending on which part of the machined surface is considered. This clearly shows the effect of the deformation of the workpiece on the geometry of the final surface; this deformation allows vibrations and the response is dependent on the mode shapes. For a given harmonic excitation, the dynamic stiffness is not the same, depending on which of the five parts of the workpiece the excitation is located; thus, the action of the tool, which is more complex, has no reason to be the same. The different parts of the generic mesh M are depicted in Fig. 24. It is observed on the swept domains that one of the teeth has a special shape that reduces undulation defaults on the machined surface. This is easily integrated in our approach by its geometric description; we do not make any assumption concerning the way the matter is erased by successive teeth. Seventeen elementary tools where used to describe the active parts of each of the 12 teeth.

machined surface

Conclusion

This paper presented a general method with which to simulate material removal in the case that both the workpiece and tool, including their holding devices, are not assumed to be rigid. It was not intended here to demonstrate that this method can give meaningful physical results.

The literature shows that this kind of approach at the macroscopic level is mature and gives satisfactory results, including the prediction of the surface structure and location errors, as far as the models of the dynamics of the WTMS and the cutting force models are sufficiently refined. The examples show that the machining of significative zones of industrial workpieces can be simulated on a laptop computer.

The difficulty when dealing with flexible workpieces is to imagine general methods to account for the relative motion of the active parts of the tool with respect to the workpiece. Our new method proposes an alternative to the methods available in the literature and has this aim. First, our method precisely defines the material removal in this general situation. It then allows us to establish a method of deriving efficient numerical approximations to simulate a whole machining operation. This is done in the context of a reduced numerical model based on a set of natural modes that describes the relative motions in dynamics. This relatively small model (having 10 to a few hundreds of DOF if required) is obtained using the FE models of the whole WTMS that may be very large. Additionally, geometric models of the swept domains and the associated dexel model of the workpiece are automatically linked to this reduced model. FE meshes that are linked to the reduced model allow us to describe the relative kinematics of the active parts of the tool with respect to an extended material frame, and thus to the dexel model of the workpiece it contains. The material removal can thus be modeled as accurately as desired by convenient choice of the degree of refinement of each model. The combination of dexels and finite elements allows us to control independently the accuracy of the geometric model of the surface by means of dexel density and the accuracy of the description of the relative motions, as usual in FE analyses, by sufficiently refined meshes and selection of the set of modes kept in the reduced model.

The evolution with time of mechanical fields and also models of the final and intermediate machined surfaces are calculated in the simulation, and thus are numerical outputs for the user that may be plotted, even during simulation. The geometrical model of the surface allows us to access its roughness, undulations and form defects, as for a real machined surface, under the assumption that the main material removal phenomenon can be described by a Boolean chip. Regenerative effects are inherently taken into account by the method without any a priori assumption about the motion of elementary tools.

As a time-domain simulation, this approach may deal with any non-linear and non-periodic or transient phenomena, as long as physical models, including cutting force models, have been established and validated. This allows us to keep the simulation as close as possible to the physical phenomena at the macroscopic level for general situations. Relaxation of initial residual stress and thermal deformations are not discussed here but they may be taken into account.

As shown in the literature, such time-domain simulation tools and stability-oriented analyses complement each other; inclusion of the time domain leads to simulations that are more time consuming than stability analyses, but stability analyses require approximations that are sometimes difficult to assess without conducting physical tests. Time-domain simulations allow us to check these approximations before conducting physical tests when expensive workpieces are to be machined.

The proposed method may be applied to complex industrial tools and workpieces, in almost any cutting process; e.g., five-axis milling [START_REF] Lorong | Simulation of a finishing operation: Milling of a turbine blade and influence of damping[END_REF], turning (Section 5.1) and face milling (Section 5.2). A workpiece may have holes, groves and any kind of particular geometric shape, and be in any position. Moreover, it is possible to simulate successive cuts, typically in studying the interaction between successive roughing operations and finishing. Under the condition that it would lead to efficient and robust algorithms, any other type of geometric model of the workpiece could be used in the material frame Rm instead of the dexel model. Our approach makes it possible to build a general mechanical core of a full VM simulator. At present, three limitations can be noted. First, the way we describe the incremental BREP model of the swept domain could degenerate if some points are near the axis of rotation (e.g. in drilling), which could lead to self-intersection. Second, we do not build a model of the flank faces of elementary tools. This prevents us from accounting for certain cutting force models including process damping that necessitate us to consider the interaction of the flank face and the workpiece. The next step of our developments will be a natural incorporation of that interaction in the method. Third, when the variation in shape of the workpiece induces significative changes in the dynamics of the workpiece, the method must account for these changes. We proposed such a method in [START_REF] Assouline | Numerical simulation of machining at the macroscopic scale: dynamic models of the workpiece[END_REF][START_REF] Lorong | Contribution of computational mechanics in numerical simulation of machining and blanking: state of the art[END_REF]. It must be validated and automated. We are now investigating the validation of such an approach for thin-walled cylinders.

Another perspective of the proposed method, and at least the concept of Φ + and the related extended material frame, is that it could be used in contact problems, when looking for neighboring points between two flexible structures undergoing unilateral contact.

Figure 1 :

 1 Figure 1: Material removal and Boolean chip (cross sections).

Figure 2 :

 2 Figure 2: Main steps of our approach.

Figure 3 :

 3 Figure 3: Flowchart of our approach. 8

Fig

  Fig.4and Fig.5show examples of T(k) and their Σ(k) . An example of a rake faceΣ (k) = C 1 C 2 C 3 C 4 is shown in Fig. 4b.

Figure 5 :

 5 Figure 5: Example of a milling cutter. Elementary tools T (k1) , T (k2) . . . T (k6) defined by a slice of the geometric model.

Figure 6 :

 6 Figure 6: Configurations I C and 0 C of the workpiece and workpiece-attached isolines of α m and β m .

Figure 9 :

 9 Figure 9: Configurations I C and t C of the (upper part of the) workpiece before material removal and workpiece-attached isolines of α m and β m .

Figure 10 :

 10 Figure 10: Cross section of (the upper part of) the workpiece in plane Π of Rm . First traces of rake faces at t i = n∆t and trajectories of cutting edges in the plane Π are shown for teeth 1, 2 and 3.

Figure 12 :

 12 Figure 12: Modeling of material removal and Boolean chips in Rm (sections in plane Π).

Figure 13 :

 13 Figure 13: Comparison of the two families of one-to-one mappings: with the rigid-body assumption (Φ R |ti is an isometry; left images) and without (Φ + |ti ; right images).

2 .

 2 Procedure c-to-m: for one point C of the tool and for any required τ ∈ [τ 1 , τ 2 ] (mainly for points C required to compute the Boolean chip), get C(τ ) from C, y c (τ ), y s (τ ) and Rc s (τ ). Steps: 1) get C s (τ ) from C (procedure c-to-s), 2) get C(τ ) from Eq. (29). Proc. 2: Pseudo-code of procedure c-to-m, used to obtain in Rm the coordinates of the image C(τ ) of a point of the tool from its position of reference C in R c .

Figure 14 :

 14 Figure 14: Partial view of the mesh M s used to construct Eq. (1a): isoparametric 10-node tetrahedron (T10) elements. The motion of mesh M + m (T4) is linked to the motion of M s .

  (a) Trajectory in Rm of the cutter center and location of the mesh M m (T4) chosen to build the one-to-one mapping Φ + . (b) A portion of M + m : on the left, faces of T4 elements modeling the boundaries of m B and, inside it, the boundaries of M m ; on the right, M m the mesh of T4 elements of the outside of the workpiece. (c) The domain m B that surrounds the machined part of the workpiece. The mesh M + m (tetrahedrons T4), shown in Fig. 15b, describes the motion of any point in m B. (d) The mesh M m (T4) of the machined part of the workpiece. The motion of M m is linked to the motion of the mesh M s of the whole workpiece shown in Fig. 14.

Figure 15 :

 15 Figure 15: Ingredients for the construction of Φ + and its inverse.

Figure 16 :

 16 Figure 16: Example of the deformation of B, as it is described by the nodes of the boundary of M + m , and of M m inside it. This deformation extends the second mode shape as given by M s outside the workpiece. It corresponds to Ψ 2|m , the second column of Ψ m .

S

  uses one intermediate position of the edges at 1 2 (t i + t i+1 ) in each time step. M (k) S describes the boundary of a closed domain that models Ω (k)

  D before = D uncut ∪ D cut and its partitions D uncut and D cut . Sets of dexels D uncut , D kept and D erased . D uncut M t i+1 (c) The set of dexels D after : D after = D uncut ∪ D kept .

Figure 18 :

 18 Figure 18: 2D illustration of the sets of dexels D before ≈ Ω W (t i ), D uncut and D cut (a), D kept and D erased (b) and D after ≈ Ω W (t i+1 ) (c).

  boolean and D (0) by setting:D (0) boolean := ∅ and D (0) := D before ; 3)for each Σ (k) (for k = 1 to k = n Σ ) a) extract C(k) Σ (τ ) for any τ ∈ I τ ; b) calculate the position of M (k) S in Rm from C(k) Σ (τ ∈ I τ ), where M (k) S is the BREP model of m Ω (k)S , the domain swept by Σ (k) during I τ ; c) cut dexels with T3 elements of M

Figure 19 :

 19 Figure 19: Example 1 -Turning of an industrial axisymmetric workpiece: partial view of the mesh M s .

Figure 20 :

 20 Figure 20: Example 1 -Position of the machined zone and close-up view of a small part of the zone.

Figure 21 :

 21 Figure 21: Example 1 -View of the machined surface S M (t E ) (interpolation of dexel ends) shown in Fig. 20.

Figure 22 :Figure 23 :

 2223 Figure 22: Example 2 -Industrial workpiece in face milling. 143,414 T4 elements, 107,900 DOF.

Figure 24 :

 24 Figure 24: Example 2 -Different parts of the generic mesh M used to model the swept domain, for each time step and each tooth, and a magnified view of the mesh for an active part of a tooth (17 elementary tools).

Figure 25 :

 25 Figure 25: Example 2 -Central part of the machined surface.

  Domain in Rm between t i and t i+1 swept by the elementary rake faces of the slice: a subset of

m Ω S (t i , t i+1 )

  , and Fig.12.

	y m	rake face, tooth 1	t 8	t 9
		rake face, tooth 2	t 9				t 0
			t 8				t 1
		rake face, tooth 3 workpiece m Ω(t 0 )	t 9	t 8		t 1	t 0	t 1	t 0
		rake face, tooth 4 first action in interval t 8 ≤ t <t 9	x m
	O m

  Boolean procedure on [t i , t i+1 ]: get D erased and D after = D(t i+1 ) from D before = D(t i ). Steps: 1) for each τ ∈ I τ calculate CΣ (τ ) from CΣ using the procedure c-to-m; 2) initialize D

is known by means of the T3 in which it is, and by its reduced coordinates in it. This gives directly the reference position

C(k) F in R c .

Acknowledgment

The examples given in this paper were extracted from real problems in the context of research collaborations with SNECMA (SAFRAN Group), PSA (Peugeot Citroën Group) and Process Conception Ingénierie S.A.