Cutting force and nonlinear chatter stability of ball-end milling cutter

Ball-end milling cutters are one of the most widely used cutters in the automotive, aerospace, die, and machine parts industries. Milling chatter reduced the surface quality and production efficiency, resulting in noise. It is particularly important to model the cutting force and analyze the flutter stability of ball-end milling cutters. In this study, a simplified milling force model of ball-end milling cutter with three degrees of freedom was established based on Merchant bevel cutting theory. The model simplified the milling force coefficient. The expressions of instantaneous milling area considering the vibration displacements in X, Y, and Z directions were derived. The nonlinear dynamic cutting force model of ball-end milling cutter with three degrees of freedom was proposed. The nonlinear chatter vibration mechanical model of ball-end milling cutter with three degrees of freedom was developed by introducing the time delay term, and the stability analysis is carried out by the stability lobe diagram. The proposed models were experimentally verified.


Introduction
Ball-end mills are widely used to produce complex surfaces. The spherical of cutter head induced the polishing of any concave surface with a radius not less than the radius of the ball [1]. However, the complex surface geometry of the workpiece required the high machining accuracy, posing a challenge to the working performance of the ball-end milling cutter [2]. In the process of milling, the geometry of the flute changes continuously along the complex tool path, causing difficulties to the calculation of the milling force. Milling force is one of the important physical parameters in CNC machining. The change of milling force directly or indirectly affects tool wear, machining error, surface topography, and cutting heat [3][4][5] and determines the power consumed in the cutting process of the CNC milling machine [6], and tool wear significantly affects milling efficiency and cost [7]. Milling force is an important basis for designing processing technology and optimizing cutting parameters [8][9][10], and it is of great significance to carry out research on the prediction of milling force.
Hossein et al. [11] analyzed and fitted the depth of cut, feed rate, and axial and tangential depth of cut and gave the relative input parameters. With the help of the statistical software package Minitab, the model analysis was performed to predict the milling force. Milfelner et al. [12] used the measured cutting force and genetic programming to establish the genetic equation of the cutting force for ballend milling. The system is used to conduct an experimental study on the cutting force monitoring in the process of ballend milling. The studies adopted experimental empirical methods to develop milling force models, and empirical formulas are subject to limited milling parameters and milling force experimental data. Yucesan et al. [13] proposed a semi-mechanical model that predicts the shear and friction load distribution of the rake face and flank face of the spiral groove of a ball-end mill. However, the determination of the cutting force coefficient required a large number of ball-end milling experiments. Moreover, the mathematical model was complex. Azeem et al. [14] proposed a simple and effective method to determine the cutting force coefficient of the ball-end milling model. The unique feature of this new method is that it can calibrate the effective empirical force coefficient under a wide range of cutting conditions by performing only one half groove cutting. Similarly, Jeong et al. [15] used instantaneous cutting force coefficients that have nothing to do with the cutting conditions, proposing a mechanical cutting force model for 3D ball-end milling. However, the cutting force coefficient was determined by experiment. Roushan et al. [16] combined the basic analytical method with orthogonal cutting finite element simulation to establish a micro-end milling force prediction model. The model takes into account the number of tool slots, milling volume, and the material properties of the workpiece. This basic analytical method does not require the support of milling parameters and milling force databases, and the developed model only needs to be verified by experiments. Merchant [17] proposed the oblique cutting theory based on orthogonal cutting and predicted the milling force. Based on the Merchant oblique cutting theory in the main section coordinate system and the generalized constant helix lead curve, Li et al. [18] proposed a general dynamic force model for a ball-end milling cutter. However, the cutting performance of the cutter edge curve studied by Li and coworkers is worse than that of the orthogonal edge curve. The orthogonal edge curve is commonly used in the actual machining process. Therefore, the cutting force model of a ball-end milling cutter with orthogonal edge curve is studied in this paper.
Arnold [19] used experiments to prove that the cutting force decreased with the increase of cutting speed. Since self-excited vibration is the difficult vibration to predict and suppress, it is very important for machining and industrial development. Based on the consideration of different processing types and the influencing factors of chatter, many scholars carried out theoretical analysis, experimental research, monitoring, and control for the problem of regenerative chatter [20][21][22][23][24][25][26]. Zhang et al. [27] established a multi-degree-of-freedom chatter model of the tool-workpiece system in modal coordinates to decouple the delay differential equation and solve the stable lobe diagram of ball-end milling by means of full discretization. The cutting tests conducted on the machining center verified the accuracy of the prediction model and observed the chatter phenomenon caused by different dominant modes. Kim et al. [28] proposed a dynamic ball-end milling force model, developed software for calculation of cutting shapes, and predicted dynamic cutting force in general CNC machining and single-track machining. The frequency spectrum of chatter during ball-end milling can be detected. Somkiat et al. [29] used Daubechies wavelet to develop a method in the ball-end milling process, which can detect chatter in the time domain and frequency domain at the same time. Experiments showed that under any cutting conditions, the chatter can be easily detected by using the proposed algorithm and the obtained parameters. Sun et al. [30] proposed a method to automatically adjust the direction of the tool axis to avoid chattering along the tool path. The above studies provide a foundation for solving the differential equations with time delay and a solution for the simulation of time domain and frequency spectrum.
Previous studies often identified milling coefficients through a large number of experiments. This study simplifies the complex milling force coefficients based on the Merchant oblique angle cutting theory and proposes a simplified three-degree-of-freedom ball-end milling cutter dynamic cutting force model, which can consider stability and instability cutting force. Based on the few studies of three-degree-of-freedom dynamic model of ball-end milling cutters, a three-degree-of-freedom ball-end milling cutter nonlinear dynamic model with time delay is established. First, according to Merchant's bevel cutting theory and the geometric model of the ball-end milling cutter with orthogonal cutting edge curves, the ball-end mill cutting force model is obtained through integration. Subsequently, the expressions of the instantaneous dynamic milling area generated by the vibration displacement in the X, Y, and Z directions are proposed. Furthermore, the dynamic cutting force model considering the vibration displacement is obtained. On this basis, a nonlinear dynamic model of the ball-end milling cutter is developed in combination with the dynamic theory. The accuracy of the proposed models is verified through experiments. Finally, the effects of the milling parameters on the milling force and the cutting stability rule of the ball-end milling cutter are discussed by the lobe diagram.

Force model of ball-end milling cutter
This study is based on the Merchant oblique cutting theory which correlates the main section coordinate system with the orthogonal cutting edge curve. The cutting force of a ball-end milling cutter is obtained from the perspective of theoretical analysis.

Differential milling force model
The geometry of the cutting edge of the ball-end mill is a complicated spatial spiral. When studying the cutting force of a ball-end milling cutter, it is usually based on the assumption of differential dispersion. The cutting edge is divided into infinitely small units called the differential cutting edge segment. Then, the complex milling process can be considered many differential oblique cutting along the cutting edge, and the force acting on the differential cutting edge segment is called a differential cutting force. Based on the Merchant oblique cutting theory, this paper establishes a differential cutting force model firstly and then develops the total cutting force through integration. Figure 1 shows the Merchant oblique cutting model. Moreover, the geometric relationship between cutting speed V, shear speed V s , and chip flow speed V c on the effective plane (PBCDEF) is displayed. In single-edge oblique cutting, the tangential force F t , the axial force F a , and radial force F r on the effective plane can be determined by [17] where ϕ e is the shear angle on the effective plane, γ e is the rake angle on the effective plane, β e is the effective friction F a = F f sin cos s − cos sin s + sin e tan s − F t tan s (3) F r = F f cos cos n + sin e tan n cos s − F t tan n cos s angle, τ s is the shear stress, φ λ is the chip flow angle on the rake face, λ s is the inclination angle, γ n is the rake angle on the normal surface, and F f is the frictional force. The frictional force F f can be expressed as where a c and a w are the chip thickness and the chip width, respectively. Because the milling process of the ball-end mill cutter can be considered as many differential oblique cutting along the cutting edge, the differential tangential force, differential axial force, and differential radial force of the differential cutting edge segment can be expressed as dF t , dF a , and dF r . The chip width a w can be expressed as the differential chip width da w . Then, the expressions of a c and da w are derived from the geometric model of the ball-end milling cutter in the following section.

The geometric model of the ball-end milling cutter
Due to the excellent cutting performance of tools with orthogonal cutting edge curves, the cutting edge curves of most spiral tools in actual production are orthogonal cutting edge curves. The cutting edge curve is expressed as where P s is the lead, z is the axial height, and φ (z) is the lag angle between the axial position z and z = 0 on the flute. Figure 2 shows the detailed geometry of the flute of the ball-end mill. As shown in Fig. 2a and b, the P (z) on the flute is the point at the axial location z, and the angular position θ i is the angle between the positive direction of the Y axis and the P (z) point in X-Y plane. In Fig. 2c, the geometry of the rake surface of the flute is obtained, and the radial rake angle and relief angle are defined in X-Y plane.
As shown in Fig. 2d, the height angle q is the angle between the positive direction of the X axis and the point P (z) in X-Z plane, and q can be expressed as where R 0 is the sphere radius of the tool measured at O. The local radius R z in X-Y plane at the point P (z) is given by According to Eq. (7), the angular position of the point P (z) can be expressed as where χ i is the rotation angle of the flute i at height z = 0, ω is the angular velocity of the spindle, φ (z)i is the lag angle of the point P (z) on the flute I, and n g is the number of teeth. Figure 3 shows the geometry of the differential cutting flute segment and chip in a two-dimensional coordinate system. The normal rake angle γ n and the normal shear angle ϕ n on the rake face of the tool are shown in Fig. 3d. From the geometry of the chip, the expressions of the chip thickness a c and the differential chip width a w are given by where a f is the feed per tooth, and the expressions of a f is The geometry of the ball-end mill and the chip thickness a c and the differential chip width da w are then used in the following section to calculate the milling forces.

Milling force model of ball-end milling cutter
Because the differential cutting force expression is the product of the cutting force coefficient and the cutting area on the effective plane, the differential cutting force can be expressed as In the main section coordinate system, the cutting characteristics of differential cutting edge segment are similar to oblique cutting. The cutting coefficient K needs to be measured by multiple experiments. Here, the simplified cutting force coefficient B is substituted, and the expression is  (13), the simplified expressions of the cutting force coefficients B r , B a , and B t are given by (11) dF t = K t a c a w , dF r = K r a c a w , dF a = K a a c a w (12)  Transforming the main profile coordinate system (X 2 , Y 2 , Z 2 ) into the global coordinate system (X, Y, Z) to obtain the equation of the differential cutting force of the ball-end milling cutter, they can be expressed as Substituting Eqs. (6) and (8) into Eq. (17), taking z as the limit of integration, the model of the ball-end milling force is obtained by integration. The expressions of the cutting forces F xi , F yi , and F zi of the flute i in the X, Y, and Z directions are given by The milling force is the sum of the milling forces of multiple teeth, and the overall milling force of the ball-end milling cutter can be expressed as Since the milling force model contains integral limits z 1 and z 2 , it is necessary to solve for integral limits z 1 and z 2 in order to solve the milling force.

Dynamic integration interval of ball-end milling
According to the position change relationship of the ballend milling cutter during the machining process, the relationship between the in-cut angle, the out-cut angle, and the axial position are established. Consequently, the integration interval under the milling state is obtained. The integral limits z 1 and z 2 can be expressed by the lag angle, and the overall milling force model of the tool is derived.
When the axial depth of cut is h a , the expressions of the axial height z a and z b are Then, the lag angles φ a and φ b of points P (a) and P (b) can be expressed as Figure 4 shows that the length of the flute involved in milling is dynamically changing when the tool is in the machining state after cutting into the workpiece. In this case, the cutting area on the X-Y plane is a semicircular sector, and the cutting edge involved in milling gradually increases from the point P (b) to P (a) . The whole action can be divided into four stages. In the first stage, the cutting edge of the tool does not participate in the meshing, and the milling force on the cutter is zero. In the second stage, the contact length of cutter-workpiece increases from 0 to the maximum. Thus, the force gradually increases. In the third stage, the cutting process is relatively stable because the contact length of cutter-workpiece reaches its maximum and remains unchanged. The fourth stage is opposite to the second stage. The contact length of the cutter-workpiece is reduced from the maximum to 0. Therefore, the cutting force is gradually reduced. The integration interval of the ball-end milling force in the milling state can be derived. As shown in Table 1, the equations of the upper and lower limits of the integral are obtained.

Milling parameters of cutting force model
It is necessary to determine the milling parameters included in the ball-end milling force model before the simulation. In the ball-end milling force model, the basic milling parameters mainly consist of friction angle β, normal shear angle ϕ n , and shear stress τ, the geometric parameters, and cutting conditions of the cutter. According to the literature [31] and the cutting data, we can get the shear stress τ and the friction angle β when the cemented carbide tool cuts the 45 steel workpiece. The friction angle β is expressed as According to Merchant's first equation, the normal shear angle ϕ n can be expressed as The milling parameters of the ball-end milling force model selected in this paper are shown in Table 2.

Dynamic milling force model of ball-end milling mutter
The chatter in the machining process is mainly caused by the change of the cutting thickness with time, and this chatter is also called regenerative chatter. The main mechanism is shown in Fig. 5. In the machining process, small disturbances of external factors cause small changes of the dynamic cutting force. The small changes act on the processing system, causing the system to produce the vibration displacement. Then, the flute i leaves chatter marks on the surface of the workpiece. After a time lag τ, the milling edge i + 1 will be milled again on the workpiece's surface with chatter marks. However, due to the presence of chatter marks, the instantaneous milling thickness of the flute i + 1 will change. Then, the dynamic cutting force will change. Finally, the milling surface will also have corrugated chatter marks.
According to the mechanism of the regeneration effect, when chatter occurs in the milling process, the vibration displacement of cutter promotes the change in the instantaneous dynamic milling thickness. Figure 6 shows the instantaneous chip shape during stable cutting. The instantaneous dynamic milling thickness is expressed as Figure 7 shows the changes in the instantaneous dynamic milling area when there is vibration displacement in the X directions. Then, the expression of the changed instantaneous dynamic milling thickness is Figure 8 shows the changes of the instantaneous dynamic milling area when there is vibration displacement in the Y directions. Then, the expression of the changed instantaneous dynamic milling thickness is    Figure 9 represents the changes of the instantaneous dynamic milling area when there is vibration displacement in the Z direction. Then, the expression of the changed instantaneous dynamic milling thickness is The expression of the cut-in angle φ a ' of the changed milling force model is Figure 10 shows the change in the instantaneous dynamic milling area when there are vibration displacements in the X, Y, and Z directions at the same time.
Then, the changed instantaneous dynamic chip thickness is expressed as Due to the occurrence of vibration displacement, both of the dynamic milling thickness and milling width change, and then the instantaneous dynamic cutting force changes. Substituting the instantaneous dynamic chip thickness (Eq. (31)) into the ball-end milling force model (Eq. (21)), nonlinear dynamic cutting force model of ball-end milling cutter with three degrees of freedom is obtained.
The dynamic cutting force of the ball-end milling cutter F x (t) in the X direction is expressed as The dynamic cutting force of the ball-end milling cutter F y (t) in the Y direction is expressed as The dynamic cutting force of the ball-end milling cutter F z (t) is expressed as

Dynamic model of milling with three degrees of freedom
Regarding the cutting regeneration effect, the general machine tool cutting system is often expressed as a delay differential equation (DDE). As shown in Fig. 8, the milling area is determined by the Q' point and the Q point. The Q' point is the position point of the current milling edge i.  where τ x , τ y , and τ z are the time lags in the X, Y, and Z directions, respectively. According to the expressions of chatter displacements Δx, Δy, and Δz, the dynamic milling forces F x (t), F y (t), and F z (t) are also functions of time t and time lag τ. According to the characteristics of milling chatter, the expressions of the time lag τ x , τ y , and τ z are given by where L is the length of the workpiece in the direction of the feed speed. Machining errors, installation eccentricity, and manufacturing errors can cause the eccentricity e of the spindle of the machine tool. The spindle generates eccentric load in the X and Y directions during rotation. The eccentric loads F 0x and F 0y can be expressed as where m is the mass. The expression m·e can be identified through experiments. Figure 11 shows a schematic diagram of the dynamics of the milling process. During the milling process of the machine tool, the cutting tool generates cutting vibration in the X, Y, and Z directions. Assuming that the damping forces of the tool subsystem and the workpiece subsystem are proportional to the vibration velocity v x , v y , and v z in their respective directions, the following coupled three degrees of freedom dynamic equations are established. The expressions are where m x , m y , and m z are the equivalent masses of the tool in the X, Y, and Z directions, respectively. c x , c y , and c z represent the equivalent damping of the tool in the X, Y, and Z directions, respectively. k x , k y , and k z represent the equivalent stiffness of the tool in the X, Y, and Z directions, respectively.

Experiment of the ball-end milling force
The experiment is shown in Fig. 12. The milling experiment in this study is carried out on the XK5032A CNC milling machine. A three-phase force sensor (CL-YD-3210) is used to measure the milling force of the three perpendicular directions (X, Y, and Z). The basic parameters of the cutter for milling force simulation and experiment are listed in Table 2. Three sets of milling test are carried out under the milling conditions, as shown in Table 3

Experiment of dynamic milling force of the ball-end cutting
In order to verify the accuracy of the dynamic milling force model, the displacement sensor is used to measure the vibration displacement of the cutter during the machining process. The measured vibration displacement is substituted into the dynamic cutting force model (Eqs. (32)-(34)), and the dynamic cutting force is obtained through simulation. Finally, the simulated dynamic cutting force is compared with the experimentally measured dynamic milling force. Figure 16 shows the experiment of the dynamic milling force measurement. A laser displacement sensor (optoNCDT2300-20) is used to measure the vibration displacement of the cutter, and a three-phase force sensor (CL-YD-3210) measured the dynamic milling force during milling. It is worth noting that due to the limitation of the experimental equipment, it is impossible to measure the vibration displacement in the Z direction. The milling process is also carried out on the XK5032A CNC milling machine.
The experiments are carried out under the same conditions (spindle speed n s of 600 r/min, feed speed V f of 90 mm/ min, and depth of cut h a of 1 mm). The results are shown in Fig. 17. Figure 17a, b, and c are the curves of the dynamic milling force F x (t), F y (t), and F z (t), respectively. The vibration displacements in the X and Y directions measured by the laser displacement sensor are substituted into the dynamic milling forces (Eqs. (32)-(34)). At this time, the ball-end milling force changes from static to dynamic. As shown in Fig. 17, the dynamic milling force curve becomes no longer smooth, and the curves obtained from the test data is basically the same as the simulation curves of the dynamic milling force model. Considering the errors of installation and clamping, it is within the error range.

Identification of the dynamic characteristic parameters
Before simulating and analyzing the dynamic model of the milling system, it is necessary to obtain the dynamic characteristic parameters of the milling system. In this paper, the modal hammer test is designed to identify the natural frequency and damping ratio of the milling system in the X, Y, and Z directions. Figure 18a and b are the detailed diagrams of the modal hammer test in the X and Z directions, respectively, which mainly describe the placement of the impact hammer and acceleration sensor. In addition, the test methods in the Y direction and the X direction are basically the same. The obtained frequency response curves in the X, Y, and Z directions are expressed in Fig. 19. Through the modal hammer test in the X, Y, and Z directions designed in this paper, the results of the natural frequency ω n and damping ratio ξ in the X, Y, and Z directions obtained by Fourier analysis are shown in Table 4.
In this paper, a stiffness test is carried out to determine the stiffness of the system. The principle is that the tool is subjected to squeezing force by the hydraulic jack, and the deformation of the tool is measured by the laser displacement sensor. As shown in Fig. 20, when the static stiffness is measured in the X and Z directions, the details of the test equipment such as hydraulic jacks, tools, force sensors, and displacement sensors are shown. Figure 21a and b show the calculation results of the stiffness in the X direction and Z direction, respectively. In addition, the stiffness in the X direction and Y direction is considered to be the same here. According to the fitting curves in Fig. 21a and b, the stiffness k x is 1.3257 × 10 4 N/mm in the X direction, and stiffness k z in the Z directions is 1.9987 × 10 5 N/mm.

Experiment and simulation of dynamic model
Due to the eccentricity e of the spindle of the machine tool, the eccentric force F 0 is produced. By collecting the displacement data during milling for analysis and identification, the product of m and e in Eq. (37) is obtained as 0.21 kg·m.
During the milling process of the ball-end milling cutter, the laser displacement sensor measures the displacement of the cutter and compares the test results with the vibration response obtained from the time domain simulation to verify the dynamic model. Because the laser displacement sensor cannot collect the vibration displacement in the Z direction, this paper uses the two sets of milling  Table 5 to verify the X and Y directions. The milling parameters in Table 2 and modal parameters determined by tests are selected. Figures 22 and 23 show the time domain diagram and the chatter displacement data measured by the displacement sensor in the X and Y directions. The response curves obtained from the simulation in the X and Y directions are basically consistent with the curves of the displacement data obtained in the experiment. Therefore, the accuracy of the dynamic model is verified based on the time-domain diagrams of the X and Y directions and the vibration displacement fitting results collected from the experiment. In the cutting process, the effects of milling parameters on the cutting force are different. Prior knowledge can help process planners choose more suitable cutting conditions to reduce excessive tool wear or damage and even ensure the machining accuracy of parts. This requires simulation and analysis of cutting force models with different milling parameters. Since this paper mainly establishes the model of the ball-end milling force during the slot milling, the axial depth of cut h a and the feed rate V f are taken as the research objects, and MATLAB simulation is performed to analyze their influence on the milling force.

The influence of the axial depth of cut on milling force
A simulation analysis of the influence of the axial depth of cut h a on the milling force is carried out. When the spindle speed n s is 600 r/min and the feed speed V f is 30 mm/min, the axial depth of cut h a is 0.3 mm, 0.5 mm, and 0.7 mm. As shown in Fig. 24, the force of the workpiece under different h a is simulated. When the axial depth of cut h a increases, the values of the milling force in the X, Y, and Z directions all increase, and the curve characteristics do not change basically. According to Eqs. (9) and (22), the increase in the axial depth of cut h a increases the milling width, which ultimately leads to the increase in the milling forces.

The influence of the feed speed on milling force
hen the spindle speed n s is 600 r/min and the axial depth of cut h a is 0.3 mm, the ball-end milling forces under different conditions of the feed speed V f of 30 mm/min, 45 mm/ min, and 60 mm/min are compared. As shown in Fig. 25a, b, and c, the changes in the ball-end milling forces F x , F y , and F z are displayed. It can be seen that when the feed rate increases, the values of the milling force increase; however, the characteristics of the curves have no change basically. According to Eqs. (9)- (11), the increase in the feed rate V f causes the feed per tooth a f to increase, so that the milling thickness a c increases. Consequently, the milling forces F x , F y , and F z increase.

The stability analysis based on the stability lobe diagram
In this study, the influence of the axial depth of cut h a on dynamic model is studied by drawing the stability lobe diagram based on the dynamic model of the ball-end milling cutter, which shows the relationship between the spindle speed n s and the axial depth of cut h a . This study analyzes the milling stability of the system based on the stability lobe diagram.
The dynamic characteristic parameters of the machine tool system are the key to drawing the stability lobe  diagram. In this section, the natural frequency, damping ratio, and stiffness determined by the experiment are selected. Figure 26 shows the drawn stability lobe diagram. X axis is the spindle speed n s , and Y axis is the axial limit depth of cut h a . The horizontal line drawn by the dashed line is the minimum limit depth of cut. According to the different properties, the stability lobe diagram in Fig. 26 can be divided into three areas: absolute stability area, unstable area, and conditional stable area. The area below the minimum cutting depth is called the absolute stable area, and the spindle speed and axial depth of cut in this area are in stable cutting. The area above the lobes drawn by the solid red line is called the unstable area, and the vibration occurs in this area. Moreover, the machining process is in unstable cutting. The area between the absolute stable zone and the unstable zone is the conditional stable zone, which is located between the horizontal line corresponding to the minimum limit depth of cut and the lobe. Then, dynamic simulations are performed on different regions of the stability lobe diagram. Figure 27 shows the time domain diagram and phase diagram of a point in the absolute stability region. Figure 27a, c, and e are the time domain diagrams in the X, Y, and Z directions, respectively. Figure 27b, d, and f are the phase diagrams in the three directions, respectively. The spindle speed n s is 800 r/min, the feed speed V f is 60 mm/min, and the axial depth of cut h a is 0.15 mm, and it is less than the limit axial depth of cut at critical stability. As shown in Fig. 27, the vibration amplitude in the three directions is very small, and the vibration amplitude changes periodically. The phase diagram is a closed circle. At this time, the tool is in the stable milling state. Figure 28 shows the dynamic simulations of the unstable zone, which are the time domain diagram and phase diagram when the spindle speed n s is 800 r/min, the feed speed V f is 60 mm/min, and the axial depth of can be seen that the axial depth of cut h a is negatively related to stable milling. In order to reduce tool vibration and improve the stability of milling, it is necessary to select a smaller axial depth of cut h a . Figure 29 shows the waterfall diagram in the X, Y, and Z directions when the spindle speed n s is 600 r/min. It can be seen that the amplitude increases with the axial depth of cut h a . The amplitude peaks at a frequency of 20 Hz, which is twice the harmonic of the spindle speed n s . This is because the ball-end milling cutter selected in this paper is a two-tooth. There are multiple frequency multiplications in the X, Y, and Z directions. Figure 30 shows the waterfall diagram in the X, Y, and Z directions when the spindle speed n s is 800r/min. It can be seen that there are frequency multiplications in the X, Y, and Z directions, and the amplitude increases with the axial depth of cut h a . The amplitude reaches the peak value

Conclusions
In this study, the theoretical modeling, simulation analysis, and experimental verification of the milling force, dynamic cutting force, and dynamics of ball-end milling were investigated. The main conclusions are as follows: (1) A three-degree-of-freedom cutting force model of ballend milling cutters is established, and the simulation results based on the theoretical model are basically consistent with the experimental results. According to the curve, it can be found that the cutting force of the ball-end milling reaches the peak value in the Y and Z directions at the same time, and the phase of the X direction lags behind the Y and Z directions. According to the test results, the influence of depth of cut and feed speed on the milling force of ball-end milling cutter is obtained. (2) The nonlinear dynamic cutting force model of the ballend mill and the nonlinear dynamic model of milling are established. The test results show that the dynamic cutting force model can accurately predict the cutting force of the ball-end milling cutter during stable and unstable milling. Moreover, the proposed nonlinear dynamic model of milling and its solution method can accurately predict the stability of milling. (3) The cutting stability rule of the ball-end milling cutter is obtained. The drawn stability lobe diagram can accurately predict the minimum limit depth of cut, and the vibration response in different regions on the stability lobe diagram has obvious differences. From the time domain diagram and phase diagram of different regions on the stability lobe diagram, it can be seen that as the axial depth of cut increases, the more unstable the system. The drawn waterfall chart shows the frequency characteristics of the ball-end milling cutter, and the waterfall charts indicate that the peak value is reached at the frequency twice the harmonic of the spindle speed and the existence of frequency multiplications.
Funding The work was supported by the Fundamental Research Funds for the Central Universities (Grant Nos. N2103003 and N2003006) and National Natural Science Foundation of China (Grant Nos. 52075087 and U1708254).

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Competing interests
The authors declare no competing interests.