Force model of freeform surface multi-axis machining with fillet end mill based on analytical contact analysis

In the multi-axis machining of freeform surface, compared with ball end mill, the fillet end mill has higher machining efficiency under the same residual height and has been widely used. As the most important physical quantity in machining process, milling force has always been the focus of research. In this paper, the contact geometry between fillet end mill and freeform surface is analyzed by analytical method, and then the milling force prediction model of multi-axis machining is established. Based on differential discretization, the cutter location point of multi-axis machining of freeform surface is approximate to multi-axis machining of oblique plane, which simplifies the research object. The inclination angle is defined to describe the relationship among cutter axis, feed, and workpiece in cutter coordinate system. The space range of the cutting edge element participating in material cutting is constructed by the swept surface of previous tool path, the to-be machined surface and the feed direction surface, and the in-cut cutting edge is determined by judging the cutting edge element one by one. Considering cutter run-out, the element cutting forces on the cylindrical and fillet surfaces of fillet end mill are derived, and all the element forces within in-cut cutting edge are summed by vector to obtain the overall milling force of fillet end mill. Simulation results show that, compared with the solid method, this contact analysis method between cutter and workpiece can take both efficiency and accuracy into account. In the machining experiment, the measured force and predicted force are consistent in trend and amplitude, which verifies the effectiveness of the milling force prediction model.


Introduction
Mechanical manufacturing technology is serving aerospace, automobile, petrochemical, and other industries, and cutting is currently the most important machining method in the field of mechanical manufacturing. In the multi-axis machining of freeform surface, fillet end mill is more and more widely used as a kind of high-efficiency machining cutter. As shown in Fig. 1, the fillet end mill has a larger cutting bandwidth L than the ball end mill under the same residual height, which can improve the machining efficiency. Fillet end mill is the product of the combination of flat end mill and ball end mill. On the basis of flat end mill, the corner is rounded, which avoids the phenomenon of cutting edge chipping in the machining process and improves cutter life. Milling force, as the most important basic physical quantity, is of great significance to machining parameter optimization, tool path planning, machining vibration, and deformation. It has been a hot topic to realize the efficient and accurate prediction of milling force in machining process. There are many factors affecting milling force, including cutter geometry and material, workpiece geometry and material, machining parameters, and tool path and cutter axis vector. For the multi-axis machining of freeform surface with fillet end mill, the contact analysis between cutter and workpiece is a difficult point in the modeling of milling force due to the continuous changes of cutter axis, feed direction, and workpiece geometry.
The main purpose of contact geometry analysis between cutter and workpiece is to determine the real-time contact between cutting edge and workpiece material under the cutter feed motion and self-rotation motion. At present, contact geometry analysis mainly includes solid modeling method, Zmap discrete method, and analytical method. The solid modeling method is the earliest application and has high simulation accuracy, but its time and space complexity are too high. Although many scholars have proposed improvements, they have not fundamentally solved the problem [1][2][3]. The Zmap method overcomes the problem of low efficiency but suffers a large loss in simulation accuracy [4][5][6]. The analytical method uses mathematical expressions to describe cutter workpiece engagement and in-cut cutting edge. This method combines efficiency and accuracy and is more suitable for complex freeform surface multi-axis machining [7][8][9]. Spence et al. [1] used solid modeling method in the 2 1/2 axis milling and creatively established the contact model between cutter and workpiece through Boolean operations. Ju et al. [2] bypassed the cutter swept volume and used solid modeling method to extract the contact area in complex surface machining, which reduced the data amount of workpiece updating and improved efficiency. Feery et al. [3] used a semidiscrete solid modeling technology to discretize cutter and cutter swept surface into countless thin slices and calculated the intersection curve between cutter and thin slices. For multi-axis machining, Aras et al. [4] proposed to represent workpiece surface as a series of discrete vectors in virtual environment, which is approximate. Combined with the identification function of logic matrix, Wei et al. [5] proposed a Zmap method to efficiently obtain the contact area between cutter and workpiece. Taking ring mill as the research object, Kiswanto et al. [6] proposed an algorithm for the contact area between cutter and workpiece by combining analytical method with discrete method. For five-axis semi-finishing, Kiswanto et al. [7] also established an analytical contact area model for flat end mill and ring mill. Combined with the geometric contour of workpiece and other conditions, Sun et al. [8] determined the space rang of the cutting edge participating in material cutting by analytical method in multi-axis machining of freeform surface. Wei et al. [9] discretized cutting depth and analytically deduced the contact area of fiveaxis machining of freeform surface with flat end mill. In a word, the above contact geometry analysis between cutter and workpiece is mainly focused on flat end mill, ball end mill, and ring mill. The analytical method of contact analysis for multi-axis machining of freeform surface with fillet end mill needs further study.
The cutting force coefficient represents the force per unit cutting area, which has a great influence on the accuracy of milling force prediction. The cutting force coefficient has many nonlinear factors, including cutter diameter, rake angle, helix angle, friction angle, cutting thickness, shear stress of workpiece material, and so on. At present, it is the most commonly used method to obtain the cutting force coefficient by trial cutting experiment. Zhang et al. [10] established the Lumped-mechanism cutting force model of flat end mill, calculated the cutting force coefficient under different cutting thicknesses, and fitted the expression of cutting force coefficient. Wang et al. [11] compared the difference of cutting force coefficients under different polynomials and pointed out that an appropriate polynomial of cutting force coefficient should be selected for force prediction. Based on the average cutting force of groove cutting with ball end mill, Wei et al. [12] expressed the shear coefficient as a polynomial function of axial position, the plowing coefficient as a constant, and fitted the cutting force coefficient by multiple linear regression. Based on the linear mechanical model, Dikshit et al. [13] emphatically analyzed the influence of rotational speed on the cutting force coefficient. In short, using instantaneous cutting force to identify cutting force coefficient is a trend, which can reflect the size effect. In five-axis machining with cutter axis constantly changing, in addition to the size effect, the difference of the cutting force coefficients of the cutting edge elements with different axial heights should be considered.
Aiming at the multi-axis machining of freeform surface with fillet end mill, this paper analyzes the contact geometry between cutter and workpiece based on analytical method and then establishes a milling force prediction model. The main contents are as follows: (1) the multi-axis machining of freeform surface is divided into multi-axis machining of oblique plane by differentiation. Based on the parametric definition of the relationship among cutter axis, feed, and workpiece position, the algorithm flow of in-cut cutting edge involved in material cutting is deduced by using space-limited method. (2) Taking into account the influence of cutter run- Fig. 1 Bandwidth comparison of fillet end mill and ball end mill out, the cutting thickness and width of cutting edge element on the cylindrical and fillet surfaces of fillet end mill are deduced, respectively, and the overall milling force prediction model of fillet end mill is established by summing all the element force vectors.
2 Cutter-workpiece contact analysis As shown in Fig. 2, the curve tool path in CNC machining process is linear interpolation, which is composed of a series of straight line segments. Based on the idea of differential discretization, freeform surface of multi-axis machining can be regarded as micro-inclined plane of multi-axis machining at each cutter location point. During the CNC machining of freeform surface, the tool path, cutter axis vector, machining parameters, and the normal direction of freeform surface at cutter location point can all be derived by CAM software post-processing, which is the basis of analyzing the contact between cutter and workpiece with analytical method.
After the above analysis, the research object can be transformed from freeform surface to inclined plane. As shown in Fig. 1, in order to realize maximum chip removing volume, the tilt angle in the five-axis machining of fillet end mill is often set to zero. In other words, the feed direction f, cutter axis p, and normal direction n of workpiece are adjusted into a plane, which is formulated as follows: The establishment of coordinate systems is the premise of milling force modeling. The cutter coordinate system O c -X c Y c Z c and workpiece coordinate system O w -X w Y w Z w are defined, respectively. The center of end plane of fillet end mill is set as the origin O c ; Z c axis is along the direction of cutter axis; X c axis is the cross production of feed direction f and normal direction n; Y c axis is determined by the rule of right-hand coordinate system. The programming coordinate system is mostly set as the workpiece coordinate system O w -X w Y w Z w , which is statically related to workpiece.
In order to describe the space position relationship among the feed direction, the cutter axis and the workpiece in the multi-axis machining of inclined plane, the inclination angle ε that the angle between normal direction n and cutter axis p is defined, and its expression is as follows:

Cutter profile and cutting edge
The geometry of fillet end mill is shown in Fig. 3. In the coordinate system, the analytical expressions of cylindrical surface and fillet surface of the cutter profile are shown in formula (3).
where R is the overall radius of cutter, r is the fillet radius, and the γ is intermediate variable.
Take the equal lead right-handed spiral cutting edge as an example, the vector expressions of cutting edge in cylindrical surface and fillet surface are as follows, respectively.
where θ is the radial position angle of cutting edge element. k is the axial position angle of cutting edge element between the line determined by cutting edge element and fillet center point O r and the line passing through point O r and parallel to cutter axis. Radial position angle θ is the difference between position angle ϕ of cutting edge and lag angle φ as shown in formula (5), and β is the normal helix angle.

In-cut cutting edge
In-cut cutting edge is the cutting edge section which is really involved in material removal. The in-cut cutting edge used to define cutting width is an important parameter for force prediction in the multi-axis machining of freeform surface, because it changes in real time with the change of tool path and cutter axis. As shown in Fig. 4, the contact between cutter and workpiece is a space triangle area on the cutter profile, which is composed of three spatial curves. Obviously, the cutting edge that is in the cutter workpiece engagement is involved in cutting. Wei et al. [14] obtained the boundary curves of cutter workpiece engagement based on the surface intersection and further solved the in-cut cutting edge, but the method is complex. In this paper, based on the space-limited method, a simpler analytical algorithm for in-cut cutting edge is proposed, which does not depend on the known cutter workpiece engagement.
The discrete cutting edge element is regarded as a spatial point (x y z) in cutter coordinate system, and participating in material cutting needs to meet the following three spatial conditions at the same time: (1) Out of the swept surface W of previous tool path. During the machining process, the change of workpiece geometry is reflected in the previous tool path. Taking the tool path, cutter axis vector, machining parameters, and cutter geometric profile as known variables, the swept surface of previous tool path of each cutter location point can be calculated analytically. When the cutting edge element is on the cylindrical surface and the fillet surface, the corresponding expressions are formula (6) and (7), respectively.
where s is the cutting step. When down milling, the symbol of s is "+"; when up milling, the symbol is "-." The γ 1 is intermediate variable, expressed as formula (8).
(2) Under to-be machined surface Q, the expression is: where d n is the cutting depth.
(3) In front of feed direction surface S. Feed direction surface S separates the cutter profile into two parts. Taking the feed direction f as front, the part of cutter behind S does not participate in cutting, while the front part may Cylinder surface  (10) and formula (11), respectively.
In-cut cutting edge is an integrated range of cutting edge; that is, the key to its solution is to determine the upper boundary point k u and the lower boundary point k d . If the upper and lower boundary points exist, there is only one set. As shown in Fig. 5, based on the above cutting conditions, the discrete cutting edge elements will be judged one by one whether they are involved in cutting from the tip of cutting edge. In the process of judging one by one, if an element does not participate in cutting and the next element participates in cutting, the k d is between the two elements. Then the dichotomy can be used to find k d accurately. At this time, the accurate k u can be approached step by step by using the dichotomy between the element in cutting and the element farthest from the tip of cutting edge. If no cutting edge element satisfies the cutting conditions, the cutting edge does not participate in cutting. Based on the algorithm flow, the in-cut cutting edge corresponding different cutter rotation positions can be obtained by changing the radial position angle of cutting edge.

Milling force model
In the process of force modeling, the discrete cutting edge element is approximated as straight-edge oblique cutting. According to Altintas's mechanical liner model, the radial, tangential, and axial forces of the cutting edge element on cutting edge j are as follows [15]: In formula (12), K r , K t and K a are radial, tangential, and axial cutting force coefficients, expressed as follows [16]: where τ s is the shear stress. β n is the normal friction angle. α n is the normal rake angle. ϕ n is the shear angle. η i is the cutting edge inclination angle. η f is the chip outflow angle. In formula (12), db is the cutting width, the calculation in the cylindrical and fillet surfaces are as shown in the following formula.
In formula (12), h(θ j , k) is the undeformed cutting thickness, which is comprehensively affected by the cutter axis vector and feed direction in multi-axis machining. Define the projection component of feed f on the normal n c of cutting edge element as the nominal cutting thickness h 1 (θ j , k) [17], and its vector form in the cutter coordinate system is as formula (15).
where f z is the feedrate per tooth; the normal vector n c of cutting edge element on cylindrical surface and fillet surface are as follows: Different from the normal cutting thickness, cutter run-out is another factor affecting the cutting thickness. As shown in Fig. 6, cutter run-out is a common phenomenon, which means that the rotary center O' does not coincide with the geometrical center O'. According to the run-out theory of flat end mill [18], Yes Yes Condition: The cutting edge element is in the space limited by surfaces W, Q and S.

Radial position angle ith cutting edge element uth
Find accurate k u (i) between element u and max(u) using Dichotomy i=i+1 Find accurate k d (i) between element u-1 and u using dichotomy run-out distance ρ and run-out angle λ are used to measure the degree of run-out. The run-out affects the effective cutting radius R' of each cutting edge, and then changes the cutting thickness of each cutting edge element. The change is expressed as the perturbed cutting thickness h 2 (θ j , k).
where the effective radius difference ΔR' between current and front cutting edge elements is as follows: Finally, the actual cutting thickness is the sum of nominal value and disturbed value.
The transformation matrix of three-way force of cutting edge element on fillet surface to cutter coordinate system is shown in formula (20). When k=π/2, the formula can be serviced for the cutting edge element on cylindrical surface.
In the in-cut cutting edge, all the element forces are summed up. Furthermore, the overall milling force of fillet end mill is obtained by further summing N cutting edges: For the multi-axis machining of oblique plane with cutter location points as the research object, the force model is established based on the cutter coordinate system O c -X c Y c Z c , and each cutter location corresponds to a cutter coordinate system. In order to measure uniformly, the milling force under different cutter location points on the tool path needs to be converted to the static workpiece coordinate system O w -X w Y w Z w . The relationship between them is: where η is the angles between the axes of the two coordinate systems. Because the feed direction f, cutter axis vector p and workpiece normal n are all known in the programming coordinate system, according to the relevant definitions in the "Cutter-workpiece contact analysis" section of this paper, the axis vectors of cutter coordinate system can be obtained in the workpiece coordinate system. Taking η xwxc as an example, the calculation formula is as follows: Now the force model of multi-axis machining of freeform surface with fillet end mill, which can be used to predict milling force under any conditions of cutter axis, feed direction, feed rate, cutting depth, cutting step, and rotating speed, is established.

Experiment and simulation
In order to verify the force prediction model of freeform surface multi-axis machining with fillet end mill and the analytical method of cutter workpiece contact established in this paper, the following simulations and experiments are arranged.

Simulation of in-cut cutting edge
Aiming at multi-axis machining of freeform surface, the contact geometry between fillet end mill and workpiece is analyzed. In order to verify the effectiveness of the analytical method, a simulation example is designed to compare with the solid modeling method. As shown in Fig. 7, based on the modeling and machining module of Unigraphics NX software, the design of workpiece geometry and the planning of tool path are executed. Parameters including cutter diameter, Fig. 6 Cutter run-out fillet radius, number of cutting edge, helix angle, inclination angle, cutting step and cutting depth are 10mm, 3mm, 4, 30°, 30°, 4mm, and 2.5mm, respectively. Six cutter location points are selected equally along the tool path shown in Fig. 7. Based on the analytical algorithm of incut cutting edge in oblique plane machining, the in-cut cutting edge of sampled cutter location points can be calculated by using MATLAB numerical analysis software. As shown in Fig. 8, the four cutting edges will sequentially cut into and out of the cutter workpiece engagement with cutter rotation, and the intersection between the cutting edge and the cutter workpiece engagement, that is, the upper and lower boundary points, are constantly changing. The change of the upper and lower boundary points for one cutter rotation cycle is shown in Fig. 8. Under the simulation condition, the cutting edges of cylindrical and fillet surfaces are involved in cutting at the same time, as shown by the blue dotted line and the black dotted line, respectively. In addition, this analytical method is efficient in calculating the in cut cutting edge, and it takes less than one second by using a personal office laptop.
The analytical method proposed in this paper is based on oblique plane machining, which will produce some errors when applied to freeform surface machining. Although the solid modeling method is inefficient in analyzing the contact between cutter and workpiece, its accuracy is high enough. In order to explore the accuracy of the analytical method, the following comparison is made. Based on the simulation machining parameters, the fillet end mill, cutting edge, and freeform surface are modeled by using Unigraphics NX software, and Boolean operation is performed to measure the intersection points of cutting edge and cutter workpiece engagement. The upper and lower boundary points and the length of in-cut cutting edge obtained by UG simulation are shown in Table 1 under the cutter rotation angle of 306°when the in-cut cutting edge is the longest. The length l is defined as l=k u −k d .
Similarly, at the radial position angle of 306°, the k u , k d and the length l analytical of in-cut cutting edge calculated by the analytical method are 110.5129°, 34.3602°, and 76.1527°, respectively. If the error of the analytical method relative to solid modeling method is defined as formula (24), the error results corresponding the six sampling points are shown in column 5 of Table 1.
The simulation results show that the analytical method of contact analysis between cutter and workpiece has high efficiency and slight accuracy loss when applied to freeform surface, so it is completely acceptable for force modeling.

Experiment of milling force prediction
The multi-axis machining experiment of freeform surface with fillet end mill, as shown in Fig. 9, has universally representative significance. The workpiece material is 45 steel. The tool path is set to Z-constant type and sine type, respectively. Unigraphics NX software is used to realize workpiece modeling and tool path planning. In this experiment, the cemented carbide fillet end mill with 4 cutting edges in the study of Wei et al. was used for cutting, and the cutter diameter, fillet radius, and helix angle are 10 mm, 2 mm, and 30°, respectively. The run-out distance is 0.0052 mm and run-out angle is 1.23rad [19]. The cutting step and cutting depth are 1 mm. The inclination angle of cutter axis is set to 30°along tool path. The experiment was completed in the five-axis  Fig. 9 Experiment of freeform surface multi-axis machining with fillet end mill Table 2 The polynomial coefficients of K t , K r , and K a  (13), the cutting force coefficient is an important parameter in the element force model, but its influence factors are numerous and the calculation is complicated. It is a common method in engineering practice to obtain cutting force coefficient based on trial cutting experiment [20]. The coefficient database corresponding to axial position angle k and cutting thickness h is obtained by changing cutter axis inclination angle and using instantaneous cutting force, as shown in Fig. 10 [19]. When the axial position angle is greater than 90°, the cutting force coefficient corresponding to the cylindrical surface has nothing to do with the axial position angle, only related to the cutting thickness. In the process of milling force prediction, the cutting force coefficients of different cutting edge elements can be extracted from the database with high efficiency and precision. The analytical expression of the three-way cutting force coefficient database is  Table 2.
Under the dynamometer coordinate system O m -X m Y m Z m , the measured force along the Z-constant tool path and the sine tool path are shown in Fig. 12 and Fig. 13, and the horizontal axis is time. Along the two tool paths, the numbers of cutter location points are 80 and 98, respectively. The tool path postprocessing file of Unigraphics NX software directly contains cutter location point coordinates and cutter axis vector p information, and the feed direction f can be calculated by connecting two adjacent cutter location points. Further, the cutter axis vector is set perpendicular to the workpiece surface, and the cutter axis vector obtained by the post-processing is exactly the normal direction n of cutter location point. For each cutter location corresponding to the multi-axis machining of oblique plane with fillet end mill, based on the theory established above, the milling force of one cutter rotation cycle can be predicted. Figure 11 shows the comparison of the cycle milling force at cutter location point synchronized with cutter rotation. The horizontal axis is cutter rotation period, and the vertical axis is cutting force amplitude. In Fig. 11, the measured results are in good agreement with the predicted results in trend, and the error in amplitude is also acceptable. The error in the direction of the maximum milling force F my is about 10%.
Based on the predicted cycle milling force of each cutter location point, the maximum and lowest cutting forces within the cycle are extracted as the maximum and minimum predictive force values of the cutter location point. The predicted milling force of the whole tool path can be obtained by connecting the predicted force of each cutter location point. The predicted milling force along the Z-constant tool path and sine tool path are shown in the lines on the right in Fig. 12 and Fig. 13, respectively. It can be seen that the predicted force and the measured force along the tool path are in good agreement in trend, the amplitude is basically consistent, and the local maximum error is also about 20%. In a word, considering various machining errors and model errors, the experiment strongly proves that the milling force prediction model of freeform surface multi-axis machining with fillet end mill established in this paper is correct.

Conclusion
Based on the idea of differential discretization, the multi-axis machining of freeform surface is transformed into multi-axis machining of oblique plane, which simplifies the research object. The inclination angle is defined to describe the space position relationship among the cutter axis, workpiece, feed direction, and cutter coordination system. Based on the spacelimited method, the judgment conditions of the cutting edge elements involved in material cutting of fillet end mill are given, and the cutting edge elements are judged one by one to obtain the complete in cut cutting edge. The analytical method for contact analysis between the freeform surface and fillet end mill can take into account both efficiency and accuracy and has a broad application prospect. Starting from straight-edge oblique cutting and considering the influence of cutter run-out on cutting thickness, the element cutting force model in the cylindrical and fillet surfaces of fillet end mill is derived. Combined with the cutting force coefficient database, the milling force prediction of multi-axis machining of freeform surface with fillet end mill can be realized by summing all the element force vectors within in-cut cutting edge.
Aiming at the multi-axis machining of freeform surface with fillet end mill, a series of simulations and experiments about the contact geometry and the milling force are set up. Compared with the solid modeling method, the analytical algorithm of in-cut cutting edge has high efficiency and slight loss of accuracy on freeform surface, and it is completely acceptable for force modeling. In the machining experiment with arbitrary tool path, the measured and predicted forces agree well in trend and amplitude, and the local maximum error is less than 20%. In a word, these simulations and experiments prove the effectiveness of the contact geometry analysis and milling force prediction model.  Data availability The datasets used in our study are available from the corresponding author on reasonable request.

Declarations
Ethical approval Not applicable.

Consent to participate Not applicable.
Consent to publish All the authors have reached agreement for publication.
Competing interests The authors declare no competing interests.