Modeling of cutting forces in five-axis milling of sculptured surfaces directly from a CAM tool path

The present work is motivated by the fact that for predicting cutting forces, which arise under various cutting conditions, workpiece-tool pairs and machining depths, the numerical methods are slower and less efficient than the analytical methods. In addition, recent developments in Computer Assisted Machining (CAM) techniques have enabled analytical methods to be applied even with a complex workpiece geometry. The present paper presents a practical and powerful analytical method which is based on the use of the toolpath file as the main information source for the machined surface. This method takes into account tool position and orientation in a five-axis milling process with tool ball-end modeled as a sphere called bull-nose. Also assumptions are made to get a good approximation in the calculation of global and local cutting forces, in the aim of developing an analytical model able to predict the cutting forces for five-axis milling process, easy to apply for any practical case.


Nomenclature
(O, X, Y, Z) Global reference system (er, eκ, eψ) Local cutting edge base vectors (E, x, y, z)  Introduction 5-axis machining strategies are widely used in the manufacturing of sculptured free-form surfaces such as propellers, turbines, and dies/molds. Modeling the cutting forces is the basis for the prediction of tool wear and breakage, machine-tool vibration, cutting parameter optimization and surface quality. To generate G-Code for a specific machine tool, the user have to provide all of the required input parameters such as cutter geometry and material, desired depth of cut, spindle speed, and feed rate. These inputs are often based on handbooks and years of experience. The prediction of cutting forces allows the enhancement of both cutting conditions and tool path but requires precise knowledge of the cutting process. A general model for predicting cutting forces which is based on analytical calculations, is developed for 5-axis machining programs to plan and optimize the machining process. A literature review shows that there have been many models developed to predict the cutting forces in milling processes. Modeling of cutting can be based on three different and complementary approaches: an analytical approach devoted to describe and use the thermomechanical phenomena occurring during cutting, a mechanistic approach attempting to calculate the cutting forces simply and effectively by using reference measurements transposed to other situations, or a numerical approach usually based on Finite Elements. Fontaine et al. (2006Fontaine et al. ( , 2007 described a complete geometrical model for 3-axis milling and used a thermomechanical oblique cutting model taking into account the part material behaviour and the friction conditions at the tool-chip interface. The capability of this thermomechanical approach for the modelling of other machining processes (turning, drilling and milling) was also demonstrated by Molinari et al (2007).
Vector clipping is a mechanistic method which was proposed in Chappel (1983). In this method, the surface is approximated with a cloud of points and the surface normal vectors are calculated at various points. To simulate the machining process, the intersection of each vector with the tool swept envelope must be calculated.
In the Z-mapping method, the corresponding height of each grid point is stored as a Zvalue. For each movement of the tool, Z-values are compared with the tool position and then updated. Anderson (1978) proposed Z-mapping in collision detection and elimination in NC machining. Kim et al. (2003) used the Z-mapping technique to determine the cutter contact area in 3-axis sculptured surface ball-end milling. In Lazoglu's study (2003), Z-mapping is used to determine the chip load. Lee et al. (2002) developed an enhanced Z-mapping method by applying the supersampling technique. As a result, both efficiency and accuracy of the Zmapping method were improved. One of the applications of Z-mapping in 5-axis milling is the one presented by Fussel et al. (2003). In that study, tool-workpiece engagement was identified using an extended Z-buffer method. They used the swept envelope of the cutter to determine the intersection between the tool and Z-buffer elements. Also, the 5-axis movement of the tool was approximated as a 3-axis motion by keeping the desired accuracy.
Voxel and Octree are 3D solid discrete methods. Voxel has a simple database. As a result, it is easy and fast to apply. Also, the workpiece stock update can be performed quickly. Karunakaran et al. (2010) applied this method to simulate material removal from the stock in milling operations. Kim et al. (2006) used the Octree method to identify cutting regions in milling operations.
In analytical methods, information on the workpiece geometry is used together with the tool position to analyze and simulate the process geometry. A Book by Shah et al. (1995) presents detailed theory and applications of surface modeling techniques for CAD/CAM of complex surface geometries. One of the analytical studies in 3-axis sculptured surface machining is the one presented by Meng Lim et al. (1995). Another book by Choi et al. (2012) presents detailed theory and applications for 5-axis sculptured surface machining. Du et al. (2005) and Bailey et al. (2002) developed analytical methods to simulate 5-axis tool motion and process geometry. Ozturk et al. (2006) analytically modeled the chip load in ballend milling of free form surfaces. Lee et al. (2000) estimated the depth of cuts by positioning the tool axis coincident with the surface normal at that point. Imani et al. (1998) modeled the cutter-workpiece engagement boundaries using a geometric simulation system which uses a commercial solid modeler as the geometric engine.
In Zeroudi et al (2010), an analytical model was presented to predict the cutting forces for 3-axis milling. In the present paper, we extend this model to 5-axis milling by using the Cutter Localisation data (CL data) coming from a CAM application to calculate the tool path, orientation and engagement in the workpiece material.

1.
Description of the tool trajectory Fig. 1 shows that the z'-axis is the main axis of the spindle, while θeq is the angle between the z-axis and the z'-axis which can be either clock-wise or counter clock-wise. The workpiece coordinate system is referenced using capital letters and is denoted (X, Y, Z) with O as the origin. The z-axis and the Z-axis have the same direction but different center point.
In the present study, we consider the displacement of the tool from point E to E' which is only due to the rotational movements around X, Y, Z respectively. As a result, phenomena such as axial eccentricity, deflection and vibrations are not taken into account.
In practice, the origin O of the workpiece is palpated by the 5-axis milling machine while the origin E of the local coordinate system is obtained by measuring the length of the tool fixed in its holder using a measuring tool apparatus. Fig. 1 shows that the purchased theoretical surface is formed by a succession of a tilted straight line segments which are created from CL points resulting from CAM application. The tilt angle of the tool δ is due to the inclination of the segments forming the trajectory Tk. It can be calculated from CLX (i,k), CLY (i,k) and CLZ (i,k) as follows: Where:

Five-axis milling machine configuration
A 5-axis milling machine refers to a machine's ability to move a tool or a part in five different axes simultaneously. Fig. 2 represents a schematic representation of the tool orientation angles.

Fig. 2
Schematic representation of the swiveling angles θ n , θ t , the equivalent angle θ eq and the elliptic deformation angle θ z'y To calculate the tool position in the global coordinate system, the transformation relationship between the tool coordinate system and the global coordinate system is expressed as follows: Then, the components (I, J, K) of the tool carrying axis are obtained via :

The swiveling and the elliptic deformation angle
The swiveling angles indicate the instantaneous tool position in the global coordinate system. In the present study, these angles are known and given by a CAM software in APT format. However, this consideration imposes to express the swiveling angles as a function of I, J and K which are directly extracted from the CL data.
To calculate the tool position in the global coordinate system, the contact between the tool and the machined surface is assumed to be a point. When the spherical part of the tool swivels around its center, it does not generate changes in the global coordinate system. Nevertheless, the change in the orientation of the tool may be observed in its cylindrical part. The two tool orientation angles, which are represented in Fig. 2, describe the tool movement. The first angle θn corresponds to a spin around the x-axis, while the second angle θt corresponding to a spin around the y-axis.
The transformation matrix developed in Section 2 is used to express the tool's new position in the global coordinate system. The components of the axis carrying the tool are used to express θn and θt as a function of I, J and K as follows: The equivalent angle θeq, which is due to the swiveling angles θn and θt, can be calculated by using the scalar product and assuming unit vector in the z and z' axes. The expression of θeq is given as follows: As function of I, J and K, the expression of θeq so obtained is given as follows: As shown in Fig. 2, the elliptical deformation angle θz'y, which is obtained by projecting the z'and y axes in the (x, y) plane, may be determined as follows : The final expression of θz'y as a function of I, J and K, is given as follows: • If 0 I  and 0 J  , then : The influence of swiveling angles The geometry of the tool is made of two parts. The first part, called the body, is a cylinder of radius R0 equal to the nominal radius of the second part which is hemispherical. The center of the hemispherical part is denoted C and the tip of the tool is denoted E. When the tool moves in 3, 4 or 5 axes configuration, the swiveling angles strongly affect the modeling variables such as the current point P and the height z along the rotation angle (see Fig. 3).
The effect of the swiveling angles can be described via the following points: • For the 3-axis configuration, the speed is zero at point E, whereas for the 5-axis configuration, the speed is zero at point E'. • On the ED1 segment, there is no difference in the volume of removed matter observed for both 3 and 5-axis configurations. The dimensions D1C, CD2, ED1, and ED2 are calculated from the triangles CD1F1 and CD2F2 with a similar method: In the interval 1 2

ED z ED 
, the dimensions G1H1 and H2G2 are deduced from the triangles F1G1H1 and F2G2H2, respectively: The coordinates (XP, YP, ZP) of point P in the global coordinate system are given by :

Tool penetration
Considering a 5-axis machining process, the penetration of the tool may be calculated for any value of the swiveling angle or the direction of the motion relative to the global coordinate system. As shown in Fig. 3, the z * may be calculated with the following formula in the interval EB z EC

Tool effective radius R*(z)
For a point P located on the envelope of the tool at a machining depth z * and along the axis of rotation, the corresponding radius * () Rz of the tool measured in the (x, y) plane is as follows:

Position of the current point P along the cutting edge
The current point P is located on a cutting edge in the local frame (x', y', z') with coordinates (x'P, y'P, z'P) for a given tooth j and an altitude z'. The following expression gives the coordinates of P in the global coordinate system :  When the tool is oriented in a given direction, it forms an ellipse or a half-ellipse with the machined surface, depending on the equivalent angle θeq. Fig. 4 illustrates different marks in the penetration process between the tool and a plane surface.
are expressed as follows : Here,

Tool engagement
The engagement of the tool in the matter must be determined in order to know the quantity of matter machined by each tooth. For each cutting edge, it is necessary to determine the position of the current point P and project it onto the uncut surface attacked by the tooth. This surface is delimited by three types of boundary conditions: • The relative position of the workpiece uncut surfaces along the three axes and specifically along the Z-axis (upper initial or pre-form surface). • The previous tool path considered without any tool deflection and with perfect surface finish (path interval Δp defined in (X, Y) plane, see Fig. 1). • The path of the previous tooth.

Description of the tool contact
As shown in Fig. 6 (a), from the point P we draw a normal unit vector denoted n , which is perpendicular to the vector f defined as the instantaneous advance vector. As shown in Fig. 6 (b), the vector f can be calculated using the x  and z  angles or using the advance by tooth between two successive CL points. the vector ( ) , , x y z ff f f may be also expressed in the global coordinate system using the following expression: (x, y, z) ( 1, ) ( 1, ) ( 1, ) ( , )

Engagement compared to the trace of the preceding tooth
In the case of non-deformed chip, the trace left by the preceding tooth is taken into account in the calculation of the chip thickness * where: • r e is the normal vector at the tool envelope at point P.
• j e is the vector of eccentricity for tooth j.
• 1 j e − is the vector of eccentricity for tooth j-1.

Cutter geometry discretization
As shown in Fig. 7, a spherical coordinate system (

41
Based on the assumption that the helix stepover is constant along the tool, both expressions (40) and (41) are valid along the tool. This assumption is an approximation, and amounts to considering the inclination angle of the edge constant along the helix. In practice, cutters of constant stepover are more widespread, and this is why this approximation is acceptable.

Total rotation of the tool *
 The rotation of the tool is locally propagated and the angle of rotation between two successive CL points can be easily determined from the local interpolation segment, the feed per tooth ft and the local feed angle φz. The total rotation of the tool between two successive points CL is noted Δθ * (i) and may be calculated for each step as follows: Here, *  is equal to the summation of the Δθ * (i)'s, which are calculated for each elementary displacement along the tool trajectory.

6.
Tool engagement in the previous machined surface 7.1 Height of previous machined surface (z sp ) To consider, as simply as possible, the engagement of the tool in the previous machined surface at distinct height z values during a series of machining passes, it is enough to compare the position of the current point P with those of the preceding trajectories. The previous surface can be described in a precise way using the peak error method or an approximate way which considers the previous surface as perfect.
The simplest idea is to consider a plane surface but in the case of a more complex surface, the upper surface can be described starting from the preceding nominal surface while considering no defects in the resulting surface. Point P is in a position of cut if zP≤ zsp, where zsp is the height of the preceding surface corresponding to x = xP and y = yP. Fig. 9 Parameters defining the trajectory of the tool for a rectilinear advance in the (X, Y, Z) plane, (case of upright machining) The height of the previous machined surface zsp can be calculated by the transversal inclination angle φΔp, which represents the slope of the machined surface in the transversal direction of the advance, as shown in Fig. 9 and Fig. 10.
The calculation of the transverse inclination angle φΔp and the height of preceding surface zsp are the same as for the case of 3-axis milling. The expressions are as follows:

Calculation of point C P
The transverse offset between two successive tool paths in the plane (X, Y) corresponds to the pitch Δp. As shown in Fig. 12, the point Cp is the intersection between the axis of the tool and the preceding trajectory of the tool in the (X, Y) plane, independently from the current point P and the height (Z = ZP). The preceding trajectory of tool is considered perfect without any influence of the tool eccentricity. At point Cp, the radius of the tool envelope is denoted RP.
The current point P is in position of cut if  PP PC R .

Fig. 11
Comparison between point P on tooth j and the previous tool path for a straight path in the plane (X, Y) at Z = Z P with Δp> 0 The distance PCP is calculated by the same method as in 3-axis milling case. The distance in absolute value is given by the following exression: Calculation of the radius R P The position of P point in the (X, Y) plane is denoted CL, and its height P CL Z .The calculation of CLP on the Tk-1 trajectory remains unchanged compared to the case of 3-axis milling. When machining up or down, the height P CL Z on the Tk-1 trajectory can be calculated as follows: ( , 1) ( , 1) ( , 1) 1 ( 1, 1) ( 1, 1) ( , 1) ( 1, 1) ( , 1) In 5-axis milling, it is necessary to calculate the position of P and Cp points along the preceding trajectory Tk-1 in order to calculate the radius RP. RP can then be calculated by distinguishing the cases illustrated in Fig. 12    ( )

Radial tool eccentricity
Tool eccentricity is a shift between the axis of rotation of the spindle and the central axis of the tool. There are two types of eccentricity: (a) axial eccentricity which is the angular shift between the axis of the tool and the axis of the spindle, and (b) radial eccentricity which is the the distance between the axis of the tool and the axis of the spindle. Radial eccentricity is modeled by a value of eccentricity denoted e. As shown in  where Δψ* is the shift angle defined in formula (41).
Since the eccentricity does not change with tool orientation, * ee →→ = , and the local frame (x', y', z') is not centered anymore at the point E' defined above, the coordinates of P point are calculated as follows: Due to the radial eccentricity of the tool, the equivalent radius Re * (zP) at P point must be calculated in order to deduce the local cutting speed Vce. In the (x', y') plane, P' is the perpendicular projection of point P on the fR CC vector, as shown in Fig. 14  To be more rigourous, the cutting speed direction is also variable, and that also affects the angles α * n et λ * s angles along the considered elementary edge. This does not affect significantly the cutting forces. Therefore, there is no real interest in taking into account such a consideration.
In the (x', y', z') frame, the normal vector r e and the vectors of eccentricity j e and 1

Global forces on the tool
Starting from the local geometry (shear angle  * , edge tilt angle λ * s, elementary cutting width dw * and thickness of the non-deformed chip * 0 t ) for each elementary cutting edge, an estimate of the local cutting forces ( r dF , dF  and dF  ) can be calculated using either a mechanistic (based on the empirical coefficients of cut) or a thermomechanical method (based on a constitutive law and friction parameters).
In order to simplify and optimize the computing time in the calculation of the cutting forces along the trajectory of the tool, the following considerations were taken during the calculations: • Calculations are made for each CL point during a rotation of the tool of 2π/Nt (Nt is the number of teeth of the tool) ; • The result obtained in a CL point can be reproduced throughout the distance between the point considered CL(i, k) and the following point CL(i+1, k) ;

14
• The elementary cutting forces r dF , dF  and dF  are calculated for each edge engaged in the workpiece ( * 0 0 t  ) for a given angular position *  and expressed in the ( * r e , * e  , * e  ) frame.