Modeling and Coefficient Identification of Cortical Bone Milling Forces of Ball-End Milling Cutter for Orthopaedic Robot

In orthopedic robot surgery, milling is the most commonly used method for cutting bone tissue. Medical ball-end milling cutters, as specialized cutting tools, interact with both cortical and cancellous bone layers, directly generating milling forces. Particularly when cutting the hard cortical bone layer, significant cutting chatter, large cutting forces, and heat are often produced. Since milling forces are typically a function of cutting parameters and bone tissue characteristics when the tool remains unchanged, accurately establishing a cortical bone milling force model and assessing the milling force in terms of suppressing cortical bone cutting chatter, reducing cutting force heat damage, and optimizing process parameters are significantly important.

In the realm of modeling milling forces, the application of multivariate regression and theoretical analysis methods is pervasive. Multivariate regression hinges on the orchestration of cutting experiments and the measurement of cutting parameters, while theoretical analysis necessitates the identification of coefficients. Presently, scholars have employed diverse mathematical modeling and analytical techniques to formulate milling force models. Within the domain of multivariate regression modeling, Noordin et al. [1] scrutinized the interrelation between cutting parameters and responses through response surface methodology. They established predictive models for cutting forces and surface roughness, scrutinizing the impact of process parameters. Wu et al. [2] utilized the least squares multivariate regression method to construct a predictive model for cutting forces based on experimental data, appraising the model's precision. Tahmasbi et al. [3] introduced a second-order linear regression model to prognosticate drilling forces in bone drilling processes. Abdul-lateef Al-Abdullah et al. [4] employed artificial neural network (ANN) methods to formulate a milling force model for cancellous bone tissue milling, with an emphasis on real-time optimization and control. Zheng et al. [5] established a milling force model using support vector regression machine learning methods, offering predictions under varying processing parameters. Rabiee et al. [6] leveraged support vector regression (SVR) in machining force prediction, allowing surgeons to anticipate forces under diverse machining parameters prior to surgery.In terms of theoretical modeling of milling forces, Plaskos et al. [7] posited a model of cortical bone milling force based on the orthogonal cutting distribution method. This model represents bone milling force as the equivalent effect of instantaneous cutting force during the rotational cycle. Shang et al. [8] developed a predictive model for estimating cutting forces during two-dimensional vibration-assisted micro-milling of bone materials, incorporating anisotropic structural characteristics and kinematic properties of the milling tool. Zheng et al. [9] proposed a theoretical model linking milling forces to parameters such as spindle speed, feed rate, and bone thickness. This model subdivides the cutting edge into cutting microelements and calculates the total cutting force using the specific cutting energy method. Liao et al. [10] devised a cutting stress model for bone materials, factoring in anisotropic characteristics and predicting cutting force coefficients by considering bone orientation, tool geometry, and edge effects. Arbabtafti et al. [11] introduced a milling force modeling method using voxel principles to simulate virtual bone and chip removal, specifically applicable to spherical cutting burrs in robot-assisted joint replacement surgeries.Concerning milling force coefficient identification, Guzel and Lazoglu [12] employed a layered milling method for coefficient identification. Wojciechowski [13] identified milling force coefficients using the slot cutting method. Aydın and Köklü [14] proposed a method utilizing orthogonal milling experiments to identify coefficients of ball-end milling cutters, acknowledging the constant shear coefficients along the Z-axis while recognizing the associated limitations in accuracy. Weng et al. [15] suggested a method for identifying the coefficients of ball-end milling cutters based on the zero frequency terms of Fourier series, though without considering the influence of the changing position of milling microelements on the coefficients. Li et al. [16] introduced a novel method amalgamating the average milling force approach with particle swarm optimization (PSO). They constructed milling force and identification models based on the average milling force method, providing a PSO algorithm to enhance the search for precise cutting force coefficients.

Due to the greater heterogeneity and non-uniformity of bone tissue, its cutting behavior may be more complex. Traditional multivariate regression methods usually rely on extensive experimental data, but when dealing with bone tissue, this may be limited as obtaining a sufficient number of samples may not be easy. Additionally, the properties of bone tissue may vary among patients, including variations in bone density and strength. In such cases, traditional multivariate regression methods may not adapt well to this variability, while theoretical analysis and parameter identification methods can better adapt to different situations. Furthermore, multivariate regression models are only related to cutting parameters. Although they can be used for parameter optimization to avoid force-heat damage, they cannot be used to analyze robot milling chatter. In contrast, theoretical analysis and coefficient identification methods can establish theoretical models related to the milling state, which can be used for the analysis of robot milling chatter stability. Therefore, this study uses theoretical analysis and employs coefficient identification methods to establish a cortical bone milling force model for the ball-end milling cutter of orthopedic surgical robots. Although this method has certain limitations, this study aims to explore its applicability within a certain range, providing a theoretical basis for studying robot milling chatter problems. First, a microelement milling model of the ball-end milling cutter was established, and from it, a milling force model was derived. Then, a milling force coefficient identification model was established based on the average milling force model. Finally, this model was used to calculate the milling forces, which were then compared and analyzed with actual values to verify the accuracy of the established milling force model and identified coefficients.

In orthopedic surgery, the medical ball-end milling cutter is a specialized tool used for milling bone tissue. Its cutting edges are evenly distributed over the surface of the ball end and are concentrated at the top of the ball end. When milling, the medical ball-end milling cutter usually chooses a down-milling method and needs to be tilted at a certain angle to overcome the limitation of zero milling force. To analyze the milling force of the medical ball-end milling cutter, it is necessary to differentially process each of its milling edges, dividing each milling edge of the ball-end milling cutter into many small milling elements, and then adding up the milling force of each element to get the instantaneous milling force of the ball-end milling cutter. As shown in Figure 1(a), this study established a rectangular coordinate system $O-XYZ$ with the top point of the ball-end milling cutter as the origin. As shown in Figure 1(b), the milling edge of the ball-end milling cutter is projected onto the $XOY$ plane to obtain the projection of the milling edge. The angle between the tangent at the top point of the milling edge and the $Y$-axis represents the rotation angle of the ball-end milling cutter, denoted as $\theta$. $\varphi$ represents the lag angle of the milling edge, which can describe the position of the milling element through the value of $\varphi$. $\psi$ represents the angle between the line connecting the milling element and the center of the ball-end milling cutter and the positive direction of the $Y$-axis.

The ball-end milling cutter can be divided into two types: constant helix lead and constant helix angle. This study chooses to use a ball-end milling cutter with a constant helix lead. Since the helix lead of the ball-end milling cutter remains constant, the spiral line at any point on the maximum milling radius and the milling edge of the tool can be projected onto the cylindrical surface of the circle where the center of the ball is located, and unfolded on the $XOZ$ plane, as shown in Figure 2. Wherein, $\beta$ represents the helix angle of the milling edge element, and $R$ represents the milling radius of the milling edge. The helix angle on the plane where the center of the ball is located is defined as the maximum helix angle, denoted as $\beta_0$. The radius of the arc on this plane represents the maximum milling radius during the milling process of the tool, denoted as $R_0$.

Since the milling force of each milling edge of a ball-end milling cutter can be considered as the vector sum of the milling force of each milling element, it is necessary to determine the equation expression of the milling edge of the ball-end milling cutter and use this expression to integrate the elemental milling forces to obtain the overall milling force. In the coordinate system of the ball-end milling cutter, assume the coordinates of any milling element are ($x$, $y$, $z$). The parameters of the milling element correspond to each other, and the position of the milling element of the ball-end milling cutter can be represented by the axial position angle, whose expression is shown in Eq. (1), wherein $k$ represents the axial position angle of the milling element.

The lag angle of the milling edge of the ball-end milling cutter is represented as:

The vector equation of the ball-end milling element is represented as:

The arc length of the milling edge is represented as:

The surface of the bone is very complex. To ensure the accuracy of the surgery, the milling path cannot be limited to a single plane, but the milling cutter must be analyzed in space. At the same time, due to the complex shape of the milling edge of the medical ball-shaped milling cutter, it is not possible to directly calculate the force on the milling edge, so it is necessary to discretize the milling edge into multiple elements and analyze the milling force of each element. As shown in Figure 3, each element on the milling edge can be decomposed into radial force, normal force, and axial force, respectively, along the radial, normal, and axial directions of the tool. The instantaneous milling force of the milling element can be represented as:

where, $K_{t e}$ represents the tangential plowing force coefficient of the milling element, $K_{r e}$ represents the radial plowing force coefficient of the milling element, $K_{a e}$ represents the axial plowing force coefficient of the milling element, $K_{t c}$ represents the tangential shearing force coefficient of the milling element, $K_{r c}$ represents the radial shearing force coefficient of the milling element, $K_{a c}$ represents the axial shearing force coefficient of the milling element, $h_i(\psi, \theta, \alpha)$ represents the deformation milling thickness, and $\mathrm{d} b$ represents the milling width of the milling element of the ball-end milling cutter.

In above equations, it can be observed that in order to obtain the milling force of the ball-end milling cutter, it is also necessary to obtain the expression of undeformed chip thickness. In the milling process of the ball-end milling cutter, since the milling speed of the cutter is much greater than the feed speed of the cutter, the milling path of the ball-end milling cutter can be approximately regarded as an arc. The undeformed chip thickness can be represented by Eq. (6), where $f_z$ represents the feed per tooth of the ball-end milling cutter.

For the convenience of analyzing the force condition of the milling element of the ball-end milling cutter during the milling process, combining with the schematic diagram of the milling force element of the ball-end milling cutter (Figure 3), the forces $\mathrm{d} F_t$, $\mathrm{d} F_r$, and $\mathrm{d} F_a$ on the milling element are decomposed in the $x$, $y$, $z$ directions, respectively, as $\mathrm{d} F_x$, $\mathrm{d} F_y$, and $\mathrm{d} F_z$, and the force conditions are:

The instantaneous resultant milling force of the ball-end milling cutter refers to the resultant force of all milling elements participating in milling when the tool is at the turning angle $\theta$, that is, the vector sum of the milling forces of the milling elements.

Assume $j$ is the $j$-th milling edge on the ball-end milling cutter, and the contact boundary of the $j$-th milling edge is $\left(\kappa_j^l, \kappa_j^u\right)$. When the turning angle of the ball-end milling cutter is $\theta$, the instantaneous milling component forces acting on the $x$, $y$, and $z$ axes of the ball-end milling cutter are $F_x$, $F_y$, and $F_z$, respectively, then the instantaneous milling component forces of the ball-end milling cutter are:

where,

In the equation, $N_f$ represents the number of milling edges of the ball-end milling cutter.

The milling force coefficient of a ball-end milling cutter refers to the milling force required per unit area of chips removed by each milling element during the milling process. Since different materials use different tools, the milling force coefficients often vary. This study proposes a method based on the average milling force model coefficient identification for medical ball-end milling cutters. This model transforms the milling force model into a functional relationship between milling coefficients and milling forces and uses the least squares linear regression method to obtain the average milling force through slot milling experiments, thus calculating the milling coefficients. In this model, the shear force coefficient is represented by the position angle of the milling element of the ball-end milling cutter, while the plowing force coefficient is a constant. Then the milling force coefficient identification model can be expressed as Eq. (11), where$K_{t c 0}, K_{t c 1}, K_{t c 2} \ldots K_{t c m}, K_{r c 0}, K_{r c 1}, K_{r c 2} \ldots K_{r c m}, K_{a c 0}, K_{a c 1}, K_{a c 2} \ldots K_{a c m}$ are all constants.

According to Eq. (8), the milling force of the $j$-th milling edge of the ball-end milling cutter can be expressed as:

Let:

Then in Eq. (12):

To facilitate the calculation of milling force coefficients, it is necessary to solve the average milling force per revolution of the ball-end milling cutter. An important characteristic of milling force is the average milling force per revolution, which can comprehensively reflect the effect of all milling elements of the ball-end milling cutter in one rotation and can also eliminate the influence of factors such as tool vibration and eccentricity. The average milling force of the ball-end milling cutter can be regarded as each milling edge element experiencing a milling region once. Therefore, the average milling force model can be expressed as:

Extracting the milling force coefficients in Eq. (15), it can be expressed as:

Eq. (16) represents the equation for calculating the average milling force of the ball-end milling cutter and can be used for the milling coefficient identification experiment. However, to obtain the specific values of the milling force coefficients, it is necessary to determine the unknowns in the equation, so a slot milling experiment is required to obtain related results. For easy data collection, the experiment was carried out using the experimental setup shown in Figure 4, which mainly consists of an orthopedic surgical robot, a six-dimensional force sensor, a bone grinding drill, and a force sensor data acquisition system. The model of the orthopedic surgical robot is EC66, with a maximum effective load of 6 Kg, a maximum end movement speed of 2.8 m/s, a working radius of 914 mm, and a joint movement range of -360° to +360°. The six-dimensional force sensor model is KWR75B, with a resolution of up to 0.03 %FS, capable of measuring a maximum force of 200 N and a maximum torque of 8 Nm in the $x$, $y$, $z$ directions. During the experiment, the data collection frequency of the six-dimensional force sensor is set to 1000 Hz, and the milling force data is collected in real-time through the data acquisition system. For accurate data acquisition, the six-dimensional force sensor is installed at the tool mounting end of the robot, and the bone grinding drill is mounted on the six-dimensional force sensor along the $Z$-axis of the robot tool, as shown in Figure 5.

In order to quickly measure the milling forces in the $x$, $y$, $z$ directions and complete the milling coefficient identification, as shown in Figure 4, the robotic arm performs a pig bone milling experiment in a vertical posture. In the experiment, only the error generated by the vertical direction gravity needs to be eliminated, and the milling forces collected by the sensor in the three directions should be the same as the milling forces at the end of the milling tool in the three directions. Due to the action of gravity, the moment arm of the force with the sensor coordinate system produces additional torque, causing the torques in the three directions collected by the sensor to differ from the actual contact action torques. The milling coefficient identification experiment only uses the milling forces in the $x$, $y$, $z$ directions and does not require torque, so the method of this study is feasible.

During the experiment, the pig bone samples were fixed on a fixture, and bone milling was performed using the robot cutting system. The milling process adopted down-milling method, performing slot milling experiments on the cortical bone layer of the pre-fabricated pig bone samples along the $Y$ direction. The spindle speed was kept constant at 5000 r/min. The milling depth and feed speed of the ball-end milling cutter were continuously varied. To ensure accurate resolution of the milling force coefficients, this study assumed that the milling force coefficients of the ball-end milling cutter are a cubic function of the position angle of the milling element and set the number of experiments to 10 to improve accuracy. During the milling process, the milling force signals were collected simultaneously, and Matlab was used for data reading and processing.

In the processing of one set of collected signals, the $X$-direction milling force data was read as shown in Figure 6(a). A low-pass filter was used to remove noise from the signal, as shown in Figure 6(b). The Envelope function was used to obtain the envelope diagram of bone milling, as shown in Figure 6(c). After accumulation, its average value was taken to obtain the instantaneous average milling force, as shown in Figure 6(d). Then, the average milling force during the milling process of the ball-end milling cutter was obtained through univariate linear regression. The processing of $Y$ and $Z$ direction signals is shown in Figure 7 and Figure 8. The milling force signals were processed using the above method, and the slot milling experimental data obtained were placed in Table 1.

Based on the data measured in the pig cortical bone slot milling experiment, 30 milling force coefficient equations can be obtained. These equations can be fitted using the least squares method to obtain the milling force coefficient calculation equation as follows:

The milling force coefficients can be obtained from Eq. (17), and by substituting the obtained milling force coefficients into Eq. (8), the milling force of any milling element can be obtained.

Through the milling force coefficient identification experiment, a complete milling force model of the ball-end milling cutter was obtained. To verify the accuracy of the established milling force model and the experimentally obtained milling force coefficients, it is necessary to calculate the trend of the milling force model and compare it with the average results of 5 actual measured milling forces. When calculating, the commonly used milling parameters in orthopedic surgery were selected: spindle speed $n$=5000 r/min, milling depth $a_p$=0.5 mm, feed speed $V_f$=60 mm/min, the axial position angle during milling was 1, and the corresponding milling coefficients were: $K_{t c}$=375.61 $\mathrm{N} / \mathrm{mm}^2$, $K_{r c}$=380.12 $\mathrm{N} / \mathrm{mm}^2$, $K_{a c}$=315.72 $\mathrm{N} / \mathrm{mm}^2$. The comparison diagrams of milling forces in each direction obtained through simulation and experiment are shown in Figure 9, Figure 10, and Figure 11.

From the above figures, it can be seen that the milling force curve obtained through calculation is basically consistent with the milling force curve obtained experimentally in terms of values and changing trends. However, there is still a certain deviation between the two curves, with the maximum deviations of the milling forces in the $X$, $Y$, and $Z$ directions being 0.5 N, 0.8 N, and 0.6 N, respectively. This is because errors in data acquisition are caused by factors such as vibrations during the milling process and the tilt of the milling cutter not reaching or exceeding the preset value. Despite these deviations, the average values and peak values of the two curves are quite stable. The established milling force model and the obtained milling force coefficients are reliable and accurate within an acceptable range of error.

This study has established a milling force model for the ball-end milling cutter of an orthopedic robot during cortical bone milling through theoretical derivation and experiments, achieving accurate prediction and in-depth analysis of the milling forces during the process. This achievement is of great significance for improving the stability of cortical bone cutting, reducing cutting damage, and optimizing cutting process parameters. The following are the main contributions and conclusions of this study:

(1) Through theoretical derivation using the microelement method, the milling force of the medical ball-end milling cutter element was obtained, and the relationship between the individual milling element force and milling parameters was established, thus constructing a milling force model for the ball-end milling cutter of the orthopedic surgical robot. This model provides a solid theoretical foundation for the analysis of robot bone tissue cutting chatter stability and parameter optimization.

(2) Based on the average milling force model, a milling coefficient identification model was proposed and experiments were conducted to identify the coefficients in the milling force model, using the least squares method for identification, resulting in milling force coefficients consistent with actual conditions. This lays an important foundation for the establishment of the milling force model of the orthopedic surgical robot.

(3) The obtained milling force model was experimentally verified. A comparative analysis of the milling force images obtained from simulation with those obtained from actual milling confirmed the high accuracy of the identified milling force model within an acceptable range of error. This verification further confirms the reliability and accuracy of the established milling force model of the ball-end milling cutter of the orthopedic robot for cortical bone milling.