# Prediction Modeling and Sensitivity Analysis of Robot Bone Milling Temperature Operated by a Doctor

Authors: Heqiang Tian; Jingbo Pan; Yu Gao; Xiaoqing Dang; Bin Tian; Debao Meng; Yanan Yao
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Bone milling is a common method in robot orthopedic surgery. However, excessive milling temperature will cause thermal necrosis of bone cells and tissues. It is necessary to carry out further research and analysis on the robot bone milling process considering the lamina milling skills of spinal surgeons and clinical practice to reduce the damage to bone cells and nearby tissues and obtain good cutting surface quality. Considering the randomness of milling parameters during operation, a prediction method of milling temperature model for ball milling cutter considering the doctor’s surgical skills was proposed based on response surface method. Because of material anisotropy and microstructure difference between the cortical bone and cancellous bone, this paper would analyze the influencing factors in different bone layers to establish the prediction model of milling temperature in the segments of cortical bone and cancellous bone. Also, the influence and distribution of milling parameters on temperature in three cutting modes such as parallel cutting mode, cross cutting mode, and vertical cutting mode in the cortical bone region were analyzed. The parameter sensitivity of the milling temperature prediction model was analyzed by the Sobol method, and the influence of the input parameters on the output milling temperature was analyzed quantitatively.


---

## Body

## 1. Introduction

In recent years, medical robot has become the research hotspot of medicine and robot. Because the bone has the characteristics of rigid structure and not easy to deform and the bone image obtained before surgery has high coincidence rate and repeatability with the actual anatomical structure during surgery, bone surgery has become the main application field of robot technology [1]. Robot orthopedic surgery commonly uses ball milling cutter for bone removal, because ball milling with sharp edge, high strength, and impact resistance can disperse the cutting force and control the cutting depth with higher cutting efficiency and better cutting quality [2]. Bone tissue is an anisotropic and heterogeneous biomass composite material composed of bioactive bone cells. In the process of bone milling, excessive milling temperature will cause thermal necrosis of bone cells and tissues [3]. The surface quality of bone milling directly affects bone regeneration and implant fixation and also plays a crucial role in the surgery success [4].It is the key to avoid excessive temperature in the process of bone cutting in modern surgery. Through theoretical analysis of the mechanism of bone cutting, the mathematical prediction model of cutting temperature can be established, and the optimal cutting parameters can be found by optimizing the prediction model, which can minimize the bone tissue damage in the process of bone processing. The bone cutting temperature model can also be used for bone processing state recognition and optimization of cutting parameters setting of each layer of bone. The related work has been paid more attention by scholars all over the world. Federspil et al. obtained optimized parameters for drilling the calvarian or mastoid bone procedure and established the first setup for robotic milling of the lateral skull base by conducting experimental milling of the skull base on two human specimens using a hexapod robot [5]. Zhang et al. deduced the thermal generation rate and other parameters in reverse by a thermal model and then used the finite element model to simulate the temperature distribution during the skull grinding, which established the basis for obtaining the optimal parameters [6]. Al‐Abdullah et al. developed an artificial neural network force model to estimate the milling force of different bone densities as a function of the milling feed rate and spindle speed based on a series of bone milling experiments with commercial artificial bones [7]. Kianmajd et al. used CAM to calculate the position of the milling cutter and the material removal rate to estimate the robot milling force along the entire path of the cutter [8]. Feldmann et al. studied the basic principle of orthogonal cutting processing of cortical bone and cutting force, temperature rise, and chip formation under two different anterior angles and six different cutting depths to establish a mechanical model [3]. Shin and Yoon established a prediction model to estimate the maximum grinding temperature and heat distribution by conducting the grinding on the fresh bovine femur [9]. Sugita et al. pointed out that the cutting temperature increased with the increase of feed speed, estimated the heat distribution inside the cortical bone, regarded the grinding head as a linear heat source moving on a semi-infinite plane (bone surface), and theoretically calculated the cutting temperature distribution [10]. Bruce et al. predicted the bone milling temperature in neurosurgery based on the thermal energy conversion equation and the finite element heat transfer model of the linear correlation between the bone milling heat and the electric energy input of the motor [11]. Al-Abdullah et al. studied the robot milling force and temperature of cancellated bone and proposed an artificial neural network model to estimate the average value of milling force and temperature [12]. The above studies explored the mechanism of bone tissue device interaction, discussed the thermal damage phenomenon of bone tissue device interaction, and proposed the prediction modeling method and parameter optimization method of bone cutting parameters.However, most of the above studies focus on the prediction model of cortical bone layer, without modeling the prediction model of cancellous bone and transition area, without considering the influence of different bone density, microstructure, and specific milling trajectory, and lack quantitative analysis of the influence of cutting parameters. At the same time, the parameters obtained and the model established do not consider the lamina milling skills of spinal surgeons and lack the basis of clinical practice.This paper studies the problems related to robot bone milling operated by a doctor with combining theoretical analysis and experiment. The spinal surgery cooperative robot for lamina milling was operated by a doctor to measure the parameters of lamina milling and analyze their laws and effects when the medical ball end milling cutter was used. Based on the experiment, the prediction model considering the doctor’s surgical skills was established and the sensitivity analysis of the effective parameters was carried out. The main study in this study is summarized as follows:(1)
In order to obtain the optimal milling parameters of bone and avoid the bone damage caused by robot milling, the influence of individual differences of bone density, microstructure characteristics, and milling methods should be considered(2)
The empirical prediction model of robot bone milling temperature was established in different bone regions (cortical and cancellated bone regions)(3)
The effect of effective robot bone milling parameters on the output milling temperature was quantitatively studied to discuss the parameter sensitivity(4)
From exploring lamina milling skills of spinal surgery cooperative robot operated by a doctor, the measurement and law research of lamina milling parameters are carried out in the process of doctor operating robot surgery and a reasonable physical description and expression among robot lamina milling skill learning parameters are established to establish a theoretical basis for safety interaction and dynamic control of robot lamina milling


## 2. Experimental Method of Bone Milling

### 2.1. Robot Bone Milling System

The experimental system of bone milling was established. The experimental device consists of a six-axis cooperative robot, a force sensor, a thermocouple thermometer, and a milling tool. The robot was used in the bone milling experiment. The robot has a payload of 6 kg, a maximum working radius of 914 mm, and a maximum moving speed of 2800 mm/s. The K-type thermocouple was used to complete bone milling temperature measurement with embedded in the bone samples, whose range is −50∼1300°C and resolution is 0.1°C. The milling tool adopts a professional spherical milling cutter with a cemented carbide material ball head diameter of 4 mm for bone resection. The experimental system of robot bone milling is shown in Figure1.Figure 1
Robot bone milling experiment system.


### 2.2. Preparation of Bone Specimen

According to the density and structure, the bone is mainly divided into cancellous bone and cortical bone [13], as shown in Figure 2. Cancellous bone structure is composed of trabeculae and pores, in which trabeculae play a supporting role and pores are filled with bone marrow. There is a very thin transition between the cortical bone and cancellous bone, which is called inner circumferential lamellae. The cortical bone is located in the outer layer of bone, which is composed of osteons and bone matrix. Osteons are parallelly arranged and buried in the matrix, and osteons are arranged on the surface of bone as the outer circumferential lamellae. Therefore, the structure of cortical bone makes the bone material anisotropic, that is, it has different thermodynamic properties in different directions. Because the interface strength between the bone and bone plate is relatively weak, it may lead to cracking during cutting. In bone cutting, the formation of bone milling debris is similar to that of unidirectional composites, which mainly depends on the cutting depth and osteon orientation [14], as shown in Figure 3.Figure 2
Microstructure of the bone.Figure 3
Three cutting modes: (a) cross cutting mode, (b) parallel cutting mode, and (c) vertical cutting mode.
(a)(b)(c)Pig femur was used as the experimental material in this paper. The pig femur was processed into a number of cortical and cancellous bone blocks with a length of 20 mm, a width of 20 mm, and a height of 6 mm. According to the bone structural characteristics, the cortical bone part can be divided into three directions such as parallel direction, cross direction, and vertical direction relative to the pig femur axis in Figure4. A thermocouple was placed inside the cutting area perpendicular to the feed direction of a bone block so as to obtain the real-time temperature of bone milling and take the average value of the highest temperature as the final temperature data.Figure 4
Cortical bone part and three cutting directions.Milling experiments were carried out according to preset experimental sequence. Under different cutting parameters (rotation speed, feed speed and, milling depth), a series of milling experiments were carried out in three cutting modes (vertical cutting mode, parallel cutting mode and cross cutting mode) in the cortical bone area and cancellous bone area. The temperature during the process of bone milling was measured by using a thermocouple; each group of experiment values was filled in Table1 to explore the influence of rotation speed, feed speed, and milling depth on milling temperature.


### 2.3. Robot Bone Milling Progress

The cooperative robot emphasizes the dominant position of doctors in the surgical process and uses force (tactile) feedback means to enable doctors and robots to jointly control surgical tools to complete the surgical process, which maintains doctors’ working habits, and enable doctors to work according to conventional surgical skills. The process of lamina milling of spinal surgeons has strict specifications, and spinal surgeons have received strict surgical skill training. In order to make the spinal surgery robot reach the skill level of doctors and improve the flexibility and safety of lamina milling of spinal surgery robot, the spinal surgery cooperative robot for lamina milling is operated by a doctor to measure the parameters of lamina milling and analyze their laws and effects, and the reasonable physical description and expression among the parameters of robot lamina milling skills are established, which can establish a theoretical basis for the safe interaction and dynamic control of robot lamina milling.In this study, doctors operate the spinal surgery cooperative robot for lamina milling, mainly use the cooperative motion of hand and arm to control the mechanical arm, and perceive the surgical operation with the help of human vision, hearing, touch, and motion sense. At the same time the tissue instrument action parameters of the end surgical tools during the doctor’s operation are continuously recorded, and the milling temperature and other parameters at different speeds, feed depth, and milling speed are extracted by stages to establish the empirical prediction model of milling skill parameters by stages. On this basis, the influence of input milling skill parameters on output parameters is clarified and quantitatively analyzed, and the best milling skill parameters are found through multiobjective optimization.


### 2.4. Measurement of Bone Milling Temperature

The working principle diagram of the thermocouple thermometer used in this paper is shown in Figure5. In order to ensure the accuracy of temperature measurement data, the temperature values of T1 and T2 of the working end 1 and reference end 2, respectively, need to be first calibrated with a constant temperature bath. The temperature of the constant temperature bath can be set to 45°C, 55°C, 65°C, 75°C, 85°C, and 98°C, and reach the set temperature after a period of time. Insert the two ends of the thermocouple thermometer into the same position in the water at the same time, read the temperature values of T1 and T2, and calculate the difference between T1 and T2 and the maximum error between the set value and the measured value. The calibration experiment data sheet is shown in Table 2. It can be seen that the displayed values of T1 and T2 are close to the set temperature, the maximum error is 0.2°C, and the maximum difference between T1 and T2 is 0.3°C. The difference is less than 1°C, indicating that the accuracy of the thermocouple thermometer is high and meets the experimental requirements.Figure 5
Working principle of the thermocouple thermometer.Table 1
Experiment sequence and data results.
Serial numberOrderA (r/min)B (mm/s)C (mm)T1 (°C)T2 (°C)T3 (°C)T4 (°C)8230001200.343.24647.64112320001200.235.438.64033.11942000750.23740.842.135653000300.34951.75346.91162000300.241.24445.7391073000750.24245.246.639.73810001200.12326.127.921293000300.145.14849.34317102000750.237.240.341.935.113112000750.136.1394133.89121000750.225.82930.423.620132000750.236.839.841.534.618142000750.237.14041.6355151000300.332.936.137.630.714162000750.339.84344.237.571710001200.3273031.624.915182000750.2374041.334.816192000750.237.140.241.7351201000300.12931.332.827In order to obtain the real bone milling temperature required to establish the prediction model, the robot uses a cemented carbide ball end milling cutter to carry out orthogonal cutting experiments and thermocouples to collect the experimental data of bone milling temperature under different cutting modes and different cutting parameters. First, the surface of the bone sample completes one cutting to obtain a flat cutting plane, so as to control the cutting depth in the bone milling experiment, and the thermocouple is inserted and fixed on the surface of the bone sample, so that the thermocouple is cut off when the tool cuts through the temperature measuring point in the actual cutting process, so as to obtain the real instantaneous cutting temperature. As shown in Figure6, the installation holes are drilled at three different positions in the milling direction on the bone block workpiece to prepare for the subsequent installation of thermocouples. Second, fix the prepared bone samples on the test bench and install thermocouples, accurately adjust the position relationship between the tool and the bone, observe the whole milling process and the surface state of bone milling, measure the change of cutting temperature with time, carry out milling experiments with different milling parameters in different areas, and collect and record the milling temperature results in real time. Figure 7 shows the thermocouple temperature measurement. The ball end milling cutter blade cuts the thermocouple during a single pass through the bone milling surface. The peak value of the pulse in Figure 8 represents the maximum temperature of the bone milling surface measured by the thermocouple. This method can estimate the actual internal temperature of the cutting point.Figure 6
Bone drilling for thermocouple installation.Figure 7
Thermocouple temperature measurement.Figure 8
Temperature signal measured by using a thermocouple.


### 2.5. Response Surface Method

The response surface method consists of two parts: experimental design and functional fitting, which determine the fitting function accuracy. When using the response surface method to discuss the mathematical relationship between design variables and responses, the factors and levels need to determined first. The number of levels is 2, while the number of factors generally does not exceed 4. All these factors are measured data. Next, select the experimental design method, create the experimental design, and determine the test operation sequence. Then, the experiment results are obtained and the experimental data are analyzed. Finally, the selected model is fitted and the model validity is analyzed. The undetermined coefficients of the linear polynomials are determined by the experimental design to fit the response relationship. Let the variabley be a function of x1,x2,…,xn, which is expressed as follows:(1)y=fx1,x2,…,xn+ε,where fx1,x2,…,xn is the response coefficient of factor x1,x2,…,xn and ε is the random error. It can be seen from the mathematical method that(2)Ey=fx1,x2,…,xn.If the relation betweenEy=fx1,x2,…,xn is known, the relation between y and x1,x2,…,xn can be found. In order to find the optimal combination of research objectives more conveniently, a reasonable model is needed to reflect the research objectives. In general, a second-order model can meet the research requirement, and the expression is shown as follows:(3)y=β0+∑i=1kβixi+∑i=1kβiixi2+∑i∑jβijxixj+ε,where y is the response function; β is the regression coefficient matrix; βi, βii, and βij are the parameter coefficients of linear, square, and interaction factor effects, respectively; xi and xj represent all the influencing factors; and ε is the error term. The second-order response model considers the influence of single factor and its secondary effects and interaction effects on the response, so it can predict the response more effectively. Finally, ANOVA can be used to verify the model.


### 2.6. Central Composite Design Method

In this study, central composite design was used to design the experimental scheme. As shown in Figure9, three types of test points (cubic points, axial points, and center points) are usually selected in the central composite design in the experiment space. The specific test points are composed as follows:(4)N=mc+mr+mo,where mc represents the cubic point. If there are k impact factors, each factor is taken as two levels (−1, 1), then there are 2k test points in mc. mr represents the axial point, which is distributed in the axial direction. If there are k impact factors, there are a total of 2k test points in mr. mo is the center point, which means the number of times to repeatedly take the zero level. In general, mo≥3. The distance between the axial point and the center point (α) is called the asterisk arm, and α is an undetermined parameter. Multiple starting points can be used in the selection. If there are k impact factors, α=2k/4 should be set, and the expected experimental function can be obtained by adjusting α.Figure 9
Central composite design diagram.In this study, the maximum cutting force and the highest surface temperature of bone milling area were taken as indexes, with the milling parameters (rotation speedA, feed speed B, and milling depth C) for the influencing factors, the second-order response surface central composite design was chosen, and the experiment design factor level table was established, where there are three factors, k = 3, so the cubic point mc=8 and the axial point mr=6. mo=6 was chosen as the center point mainly to estimate the experimental error and where α≈1. Table 3 shows the arrangement of experiment points in the center of composite design.Table 2
Calibration results of the thermocouple thermometers.
No.Temperature setting (°C)T1 (°C)T2 (°C)T2−T1 (°C)Maximum error (°C)14545.144.80.30.225544.955.00.10.136565.264.80.40.247575.175.20.10.258585.085.10.10.169898.198.00.10.1In this paper, three factors (rotation speedA, feed speed B, and milling depth C) affecting the milling force and temperature were mainly studied. Table 4 is the variation level table of milling experiment parameters, and Table 5 is the coding table of the variation range of experiment factors and the level of quadratic regression design factors.Table 3
Arrangement of experiment points for central composite design.
ItemNumberNumber of factorsP3Cubic pointmc8Axial pointmr6Center pointmo6Center point number (orthogonal)9Total experiment number (uniform precision)20Total experiment number (orthogonal)23Table 4
Level changes of milling parameters.
xUnder levelAbove levelZero levelRotation speedA (r/min)100030002000Feed speedB (mm/min)3012075Milling depthC (mm)0.10.50.3


### 2.7. Sobol Method

The Sobol method is a variance-based sensitivity analysis [15]. The variance-based method depends on the following number of estimates:(5)S=VxEy|xVy,where y is the output variable, x is the input variable, E(y | x) is the expected value of y at a fixed value of x, the variance takes all possible values of x, and V(y) is the total variance of the output y.Let the model bey=fxx=x1,x2,…,xk, where xi is uniformly distributed, and f2x is integrable in the domain, which can be decomposed into(6)fx=f0+∑i=1nfixi+∑i=1nfi≤j≤nxi,xj+⋯+f1,…kxi1,xi2,…,xik.Therefore, based on the Sobol method, the total variance of modely in this study is decomposed into the variance of a single parameter and the variance generated by the interaction between parameters. The total variance V of model y can be decomposed as follows:(7)V=∑i=1nVi+∑i≤j≤nnVij+⋯+V1,…,k,where Vi represents the folk square difference generated by parameter xi, Vij represents the variance of the interaction between two parameters xi and xj, and Vk represents the variance of k parameters acting together.The specific steps of Sobol global sensitivity analysis method are as follows: first, the parameters that need sensitivity analysis are selected. In this paper, milling parameters (rotation speedA, feed speed B, and milling depth C) of ball milling cutter are presented. Then, the variation range of each parameter is determined. The range of rotation speed A is 1000–3000 r/min, the range of feed speed B is 30–120 mm/min, and the range of milling depth C is 0.1–0.3 mm. The distribution of each input milling parameter is uniformly distributed. The mean value f0 and the total variance V of each parameter are calculated by the Monte Carlo method as follows:(8)f0=∫fxdx≈1N∑k=1Nfxk,(9)V=∫f2xdx−f02≈1N∑k=1Nf2xk−f02,where the average value f0 must be a constant, N is the size of Monte Carlo sampling, and the number of samples in this study was selected as 3000. By fixing the value of parameters and changing other parameters, the folk square difference or first-order effect Vi of each parameter was calculated, that is, the influence degree of a single parameter on the model results was defined as follows:(10)Vi=∫fidxi.The influence degree of the interaction between parameters on the model results can be expressed by the following equation:(11)Vi1,i2,…,ik=∫fi1,i2,…,ik2dxi1,dxi2,…,dxik.Let the single parameter of the objective function model and the sensitivity indexS of the interaction between the parameters be expressed as:(12)Si1,…,is=Vi1,…,isV,1≤i1<…<is≤k.Therefore, for parameterxi, the first-order sensitivity index can be obtained as Si=Vi/V, which measures the main influence of parameter xi on the output model result y, that is, the partial contribution of parameter xi to the variance Vy. Similarly, for i≠j, the second-order sensitivity index of Vy is Sij=Vij/V, which measures the influence of the interaction between parameters xi and xj on the output model result y.After normalization of equation (12), the following equation can be obtained:(13)1=∑i=1NSi+∑i≤j≤nNSij+⋯+S1,2,…,k.The global sensitivity of each input parameter refers to the sensitivity of all the parameters. The global sensitivity definition index can be expressed as follows:(14)STi=Si+∑i≠jSij+⋯.


## 2.1. Robot Bone Milling System

The experimental system of bone milling was established. The experimental device consists of a six-axis cooperative robot, a force sensor, a thermocouple thermometer, and a milling tool. The robot was used in the bone milling experiment. The robot has a payload of 6 kg, a maximum working radius of 914 mm, and a maximum moving speed of 2800 mm/s. The K-type thermocouple was used to complete bone milling temperature measurement with embedded in the bone samples, whose range is −50∼1300°C and resolution is 0.1°C. The milling tool adopts a professional spherical milling cutter with a cemented carbide material ball head diameter of 4 mm for bone resection. The experimental system of robot bone milling is shown in Figure1.Figure 1
Robot bone milling experiment system.


## 2.2. Preparation of Bone Specimen

According to the density and structure, the bone is mainly divided into cancellous bone and cortical bone [13], as shown in Figure 2. Cancellous bone structure is composed of trabeculae and pores, in which trabeculae play a supporting role and pores are filled with bone marrow. There is a very thin transition between the cortical bone and cancellous bone, which is called inner circumferential lamellae. The cortical bone is located in the outer layer of bone, which is composed of osteons and bone matrix. Osteons are parallelly arranged and buried in the matrix, and osteons are arranged on the surface of bone as the outer circumferential lamellae. Therefore, the structure of cortical bone makes the bone material anisotropic, that is, it has different thermodynamic properties in different directions. Because the interface strength between the bone and bone plate is relatively weak, it may lead to cracking during cutting. In bone cutting, the formation of bone milling debris is similar to that of unidirectional composites, which mainly depends on the cutting depth and osteon orientation [14], as shown in Figure 3.Figure 2
Microstructure of the bone.Figure 3
Three cutting modes: (a) cross cutting mode, (b) parallel cutting mode, and (c) vertical cutting mode.
(a)(b)(c)Pig femur was used as the experimental material in this paper. The pig femur was processed into a number of cortical and cancellous bone blocks with a length of 20 mm, a width of 20 mm, and a height of 6 mm. According to the bone structural characteristics, the cortical bone part can be divided into three directions such as parallel direction, cross direction, and vertical direction relative to the pig femur axis in Figure4. A thermocouple was placed inside the cutting area perpendicular to the feed direction of a bone block so as to obtain the real-time temperature of bone milling and take the average value of the highest temperature as the final temperature data.Figure 4
Cortical bone part and three cutting directions.Milling experiments were carried out according to preset experimental sequence. Under different cutting parameters (rotation speed, feed speed and, milling depth), a series of milling experiments were carried out in three cutting modes (vertical cutting mode, parallel cutting mode and cross cutting mode) in the cortical bone area and cancellous bone area. The temperature during the process of bone milling was measured by using a thermocouple; each group of experiment values was filled in Table1 to explore the influence of rotation speed, feed speed, and milling depth on milling temperature.


## 2.3. Robot Bone Milling Progress

The cooperative robot emphasizes the dominant position of doctors in the surgical process and uses force (tactile) feedback means to enable doctors and robots to jointly control surgical tools to complete the surgical process, which maintains doctors’ working habits, and enable doctors to work according to conventional surgical skills. The process of lamina milling of spinal surgeons has strict specifications, and spinal surgeons have received strict surgical skill training. In order to make the spinal surgery robot reach the skill level of doctors and improve the flexibility and safety of lamina milling of spinal surgery robot, the spinal surgery cooperative robot for lamina milling is operated by a doctor to measure the parameters of lamina milling and analyze their laws and effects, and the reasonable physical description and expression among the parameters of robot lamina milling skills are established, which can establish a theoretical basis for the safe interaction and dynamic control of robot lamina milling.In this study, doctors operate the spinal surgery cooperative robot for lamina milling, mainly use the cooperative motion of hand and arm to control the mechanical arm, and perceive the surgical operation with the help of human vision, hearing, touch, and motion sense. At the same time the tissue instrument action parameters of the end surgical tools during the doctor’s operation are continuously recorded, and the milling temperature and other parameters at different speeds, feed depth, and milling speed are extracted by stages to establish the empirical prediction model of milling skill parameters by stages. On this basis, the influence of input milling skill parameters on output parameters is clarified and quantitatively analyzed, and the best milling skill parameters are found through multiobjective optimization.


## 2.4. Measurement of Bone Milling Temperature

The working principle diagram of the thermocouple thermometer used in this paper is shown in Figure5. In order to ensure the accuracy of temperature measurement data, the temperature values of T1 and T2 of the working end 1 and reference end 2, respectively, need to be first calibrated with a constant temperature bath. The temperature of the constant temperature bath can be set to 45°C, 55°C, 65°C, 75°C, 85°C, and 98°C, and reach the set temperature after a period of time. Insert the two ends of the thermocouple thermometer into the same position in the water at the same time, read the temperature values of T1 and T2, and calculate the difference between T1 and T2 and the maximum error between the set value and the measured value. The calibration experiment data sheet is shown in Table 2. It can be seen that the displayed values of T1 and T2 are close to the set temperature, the maximum error is 0.2°C, and the maximum difference between T1 and T2 is 0.3°C. The difference is less than 1°C, indicating that the accuracy of the thermocouple thermometer is high and meets the experimental requirements.Figure 5
Working principle of the thermocouple thermometer.Table 1
Experiment sequence and data results.
Serial numberOrderA (r/min)B (mm/s)C (mm)T1 (°C)T2 (°C)T3 (°C)T4 (°C)8230001200.343.24647.64112320001200.235.438.64033.11942000750.23740.842.135653000300.34951.75346.91162000300.241.24445.7391073000750.24245.246.639.73810001200.12326.127.921293000300.145.14849.34317102000750.237.240.341.935.113112000750.136.1394133.89121000750.225.82930.423.620132000750.236.839.841.534.618142000750.237.14041.6355151000300.332.936.137.630.714162000750.339.84344.237.571710001200.3273031.624.915182000750.2374041.334.816192000750.237.140.241.7351201000300.12931.332.827In order to obtain the real bone milling temperature required to establish the prediction model, the robot uses a cemented carbide ball end milling cutter to carry out orthogonal cutting experiments and thermocouples to collect the experimental data of bone milling temperature under different cutting modes and different cutting parameters. First, the surface of the bone sample completes one cutting to obtain a flat cutting plane, so as to control the cutting depth in the bone milling experiment, and the thermocouple is inserted and fixed on the surface of the bone sample, so that the thermocouple is cut off when the tool cuts through the temperature measuring point in the actual cutting process, so as to obtain the real instantaneous cutting temperature. As shown in Figure6, the installation holes are drilled at three different positions in the milling direction on the bone block workpiece to prepare for the subsequent installation of thermocouples. Second, fix the prepared bone samples on the test bench and install thermocouples, accurately adjust the position relationship between the tool and the bone, observe the whole milling process and the surface state of bone milling, measure the change of cutting temperature with time, carry out milling experiments with different milling parameters in different areas, and collect and record the milling temperature results in real time. Figure 7 shows the thermocouple temperature measurement. The ball end milling cutter blade cuts the thermocouple during a single pass through the bone milling surface. The peak value of the pulse in Figure 8 represents the maximum temperature of the bone milling surface measured by the thermocouple. This method can estimate the actual internal temperature of the cutting point.Figure 6
Bone drilling for thermocouple installation.Figure 7
Thermocouple temperature measurement.Figure 8
Temperature signal measured by using a thermocouple.


## 2.5. Response Surface Method

The response surface method consists of two parts: experimental design and functional fitting, which determine the fitting function accuracy. When using the response surface method to discuss the mathematical relationship between design variables and responses, the factors and levels need to determined first. The number of levels is 2, while the number of factors generally does not exceed 4. All these factors are measured data. Next, select the experimental design method, create the experimental design, and determine the test operation sequence. Then, the experiment results are obtained and the experimental data are analyzed. Finally, the selected model is fitted and the model validity is analyzed. The undetermined coefficients of the linear polynomials are determined by the experimental design to fit the response relationship. Let the variabley be a function of x1,x2,…,xn, which is expressed as follows:(1)y=fx1,x2,…,xn+ε,where fx1,x2,…,xn is the response coefficient of factor x1,x2,…,xn and ε is the random error. It can be seen from the mathematical method that(2)Ey=fx1,x2,…,xn.If the relation betweenEy=fx1,x2,…,xn is known, the relation between y and x1,x2,…,xn can be found. In order to find the optimal combination of research objectives more conveniently, a reasonable model is needed to reflect the research objectives. In general, a second-order model can meet the research requirement, and the expression is shown as follows:(3)y=β0+∑i=1kβixi+∑i=1kβiixi2+∑i∑jβijxixj+ε,where y is the response function; β is the regression coefficient matrix; βi, βii, and βij are the parameter coefficients of linear, square, and interaction factor effects, respectively; xi and xj represent all the influencing factors; and ε is the error term. The second-order response model considers the influence of single factor and its secondary effects and interaction effects on the response, so it can predict the response more effectively. Finally, ANOVA can be used to verify the model.


## 2.6. Central Composite Design Method

In this study, central composite design was used to design the experimental scheme. As shown in Figure9, three types of test points (cubic points, axial points, and center points) are usually selected in the central composite design in the experiment space. The specific test points are composed as follows:(4)N=mc+mr+mo,where mc represents the cubic point. If there are k impact factors, each factor is taken as two levels (−1, 1), then there are 2k test points in mc. mr represents the axial point, which is distributed in the axial direction. If there are k impact factors, there are a total of 2k test points in mr. mo is the center point, which means the number of times to repeatedly take the zero level. In general, mo≥3. The distance between the axial point and the center point (α) is called the asterisk arm, and α is an undetermined parameter. Multiple starting points can be used in the selection. If there are k impact factors, α=2k/4 should be set, and the expected experimental function can be obtained by adjusting α.Figure 9
Central composite design diagram.In this study, the maximum cutting force and the highest surface temperature of bone milling area were taken as indexes, with the milling parameters (rotation speedA, feed speed B, and milling depth C) for the influencing factors, the second-order response surface central composite design was chosen, and the experiment design factor level table was established, where there are three factors, k = 3, so the cubic point mc=8 and the axial point mr=6. mo=6 was chosen as the center point mainly to estimate the experimental error and where α≈1. Table 3 shows the arrangement of experiment points in the center of composite design.Table 2
Calibration results of the thermocouple thermometers.
No.Temperature setting (°C)T1 (°C)T2 (°C)T2−T1 (°C)Maximum error (°C)14545.144.80.30.225544.955.00.10.136565.264.80.40.247575.175.20.10.258585.085.10.10.169898.198.00.10.1In this paper, three factors (rotation speedA, feed speed B, and milling depth C) affecting the milling force and temperature were mainly studied. Table 4 is the variation level table of milling experiment parameters, and Table 5 is the coding table of the variation range of experiment factors and the level of quadratic regression design factors.Table 3
Arrangement of experiment points for central composite design.
ItemNumberNumber of factorsP3Cubic pointmc8Axial pointmr6Center pointmo6Center point number (orthogonal)9Total experiment number (uniform precision)20Total experiment number (orthogonal)23Table 4
Level changes of milling parameters.
xUnder levelAbove levelZero levelRotation speedA (r/min)100030002000Feed speedB (mm/min)3012075Milling depthC (mm)0.10.50.3


## 2.7. Sobol Method

The Sobol method is a variance-based sensitivity analysis [15]. The variance-based method depends on the following number of estimates:(5)S=VxEy|xVy,where y is the output variable, x is the input variable, E(y | x) is the expected value of y at a fixed value of x, the variance takes all possible values of x, and V(y) is the total variance of the output y.Let the model bey=fxx=x1,x2,…,xk, where xi is uniformly distributed, and f2x is integrable in the domain, which can be decomposed into(6)fx=f0+∑i=1nfixi+∑i=1nfi≤j≤nxi,xj+⋯+f1,…kxi1,xi2,…,xik.Therefore, based on the Sobol method, the total variance of modely in this study is decomposed into the variance of a single parameter and the variance generated by the interaction between parameters. The total variance V of model y can be decomposed as follows:(7)V=∑i=1nVi+∑i≤j≤nnVij+⋯+V1,…,k,where Vi represents the folk square difference generated by parameter xi, Vij represents the variance of the interaction between two parameters xi and xj, and Vk represents the variance of k parameters acting together.The specific steps of Sobol global sensitivity analysis method are as follows: first, the parameters that need sensitivity analysis are selected. In this paper, milling parameters (rotation speedA, feed speed B, and milling depth C) of ball milling cutter are presented. Then, the variation range of each parameter is determined. The range of rotation speed A is 1000–3000 r/min, the range of feed speed B is 30–120 mm/min, and the range of milling depth C is 0.1–0.3 mm. The distribution of each input milling parameter is uniformly distributed. The mean value f0 and the total variance V of each parameter are calculated by the Monte Carlo method as follows:(8)f0=∫fxdx≈1N∑k=1Nfxk,(9)V=∫f2xdx−f02≈1N∑k=1Nf2xk−f02,where the average value f0 must be a constant, N is the size of Monte Carlo sampling, and the number of samples in this study was selected as 3000. By fixing the value of parameters and changing other parameters, the folk square difference or first-order effect Vi of each parameter was calculated, that is, the influence degree of a single parameter on the model results was defined as follows:(10)Vi=∫fidxi.The influence degree of the interaction between parameters on the model results can be expressed by the following equation:(11)Vi1,i2,…,ik=∫fi1,i2,…,ik2dxi1,dxi2,…,dxik.Let the single parameter of the objective function model and the sensitivity indexS of the interaction between the parameters be expressed as:(12)Si1,…,is=Vi1,…,isV,1≤i1<…<is≤k.Therefore, for parameterxi, the first-order sensitivity index can be obtained as Si=Vi/V, which measures the main influence of parameter xi on the output model result y, that is, the partial contribution of parameter xi to the variance Vy. Similarly, for i≠j, the second-order sensitivity index of Vy is Sij=Vij/V, which measures the influence of the interaction between parameters xi and xj on the output model result y.After normalization of equation (12), the following equation can be obtained:(13)1=∑i=1NSi+∑i≤j≤nNSij+⋯+S1,2,…,k.The global sensitivity of each input parameter refers to the sensitivity of all the parameters. The global sensitivity definition index can be expressed as follows:(14)STi=Si+∑i≠jSij+⋯.


## 3. Results and Discussion

### 3.1. Prediction Model of Milling Temperature

Based on the central composite experimental design method, the central composite in response surface model was selected for experiment design. The rotation speedA, feed speed B, and milling depth C are filled in Table 3, and 20 experiment design points were automatically analyzed and constructed by Design-Expert software. In the cortical bone area, milling experiments were carried out in three cutting directions: parallel cutting direction, cross cutting direction, and vertical cutting direction. The milling temperatures are T1 in parallel cutting mode, T2 in cross cutting mode, T3 in vertical cutting mode, and T4 in the cancellous area, respectively. The sequence of experiment points and experiment data results are shown in Table 1.Table 5
Range of test factors and coding table of quadratic regression design.
Range of variablexActual variablexStandardization ofzVariablex1x2x3z1z2z3Above levelx12x22x32111Zero levelx10x20x3000Under levelx11x21x31−1−1Distance fromr horizontalr×Δ1+x10r×Δ2+x20r×Δ3+x30rrrDistance from −r horizontalr×Δ1+x10r×Δ2+x20r×Δ3+x30−r−r−rStandard deviationΔ1Δ2Δ3Considering the influence of bone microstructure and milling trajectory, the empirical mathematical models of milling temperature were established in cortical and cancellous bone segments. The reliability was set as 95%, andP<0.05 was the level of statistical significance. The variance analysis of effective parameters was used to compare the differences of different parameters and combinations to verify the validity of different model items of milling parameters. The prediction model of the maximum milling temperatures T1, T2, and T3 under parallel cutting mode, cross cutting mode, and vertical cutting mode in the cortical bone region and the prediction model of milling temperature T4 in cancellous bone regions are shown as follows:(15)T1=16.605+0.021A−0.16B−17.452C−8.471×10−2AB+2.5×10−4AC+0.011BC−3.155×10−6A2+6.15B2+89.545C2,(16)T2=17.333+0.021A−0.137B−6.224C−5.555×10−7AB−1.5×10−3AC−0.022BC−3.123×10−6A2+5.32×10−4B2+77.727C2,(17)T3=19.2+0.021A+0.026B−0.138C−2.777×10−7AB−1.125×10−3AC−0.025BC−3.236×10−6A2+5.499×10−4B2+86.364C2,(18)T4=16.605+0.021A−0.16B−17.452C−8.471×10−2AB+2.5×10−4AC+0.011BC−3.155×10−6A2+6.15B2+89.545C2,where A represents rotation speed, B represents feed speed, and C represents milling depth.


### 3.2. Prediction Result Analysis

Milling temperatureT1 in the parallel cutting mode, temperature T2 in the cross cutting mode, milling temperature T3 in the vertical cutting mode were collected, respectively, and then milling temperature T4 was collected in cancellous bone region. Table 6 shows the results of variance analysis of quadratic regression equation prediction model of milling temperature under the parallel cutting mode of cortical bone area by ANOVA module, where the unimportant model items were eliminated. It can be seen from Table 6 that, in the influence factors of the model, rotation speed A and feed speed B are significant, milling depth C is not significant, and quadratic terms such as AB, AC, BC, A2, B2, and C2 are not significant. The misfit term in the table is used to indicate the degree of difference between the model and the experiment. In this fitting, P>0.05, the misfit term is not significant, which is beneficial to the model, and there is no misfit factors, indicating that the prediction model can replace the simulation point to analyze the simulation results. It can be seen from Table 6 that the fitting effect of quadratic multinomial model is significant and rotation speed A and feed speed B are the significant factors influencing milling temperature, and rotation speed A has the most significant influence.Table 6
Analysis of variance for quadratic regression equationT1.
FactorsSum of squares of deviationsDegrees of freedomMean square errorFP (Prob > F)Model803.78989.316001.180.0201SignificantRotation speedA649.641649.6443653.000.0098SignificantFeed speedB87.62187.625887.450.0054SignificantMilling depthC38.81138.812607.810.0136SignificantAB0.00010.0000.0001.0000AC5.000E − 00315.000E − 0030.340.5750BC0.02010.0201.340.2733A227.37127.371838.87<0.0001B24.2714.27286.64<0.0001C22.2112.21148.17<0.0001Residual0.15100.015Misfit term0.05550.0110.590.7090Not significantPure error0.09350.019Figure10(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T1, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 10(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Figure 10
Residual analysis of prediction equationT1. (a) Residual normal probability distribution and (b) relationship between the predicted value and the actual value.
(a)(b)Figure11 shows the three-dimensional response curves of the influence of milling parameters on temperature T1. It can be clearly seen from Figure 11(a) that, within the range of study parameters, the milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. It can be concluded that the highest region of milling temperature is concentrated in the region with high rotational speed and low feed speed. It can be seen from Figure 11(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 11(c) that, within the range of study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 11
Response surface curves of the influence of milling parameters onT1. (a) Effect of rotation speed A and feed rate B on T1, (b) effect of rotation speed A and milling depth C on T1, and (c) effect of feed speed B and milling depth C on T1.
(a)(b)(c)Table7 shows the results of analysis of variance (ANOVA) for the prediction model of milling temperature T2 in the cross cutting mode of cortical bone region, respectively.Table 7
Analysis of variance for quadratic regression equationT2.
FactorsSum of squares of deviationsDegrees of freedomMean square errorFP (Prob > F)Model798.09988.681137.670.0361SignificantRotation speedA649.641649.648334.490.0081SignificantFeed speedB79.52179.521020.250.0324SignificantMilling depthC40.80140.80523.490.0101SignificantAB5.000E − 00315.000E − 0030.0640.4933AC0.1810.182.310.2976BC0.08010.0801.030.5736A226.82126.82344.040.0004B23.1913.1940.940.0007C21.6611.6621.320.1472Residual0.78100.078Misfit term0.1750.0340.280.4084SignificantPure error0.6150.12Figure12(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T2, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 12(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Figure 12
Residual analysis of prediction equationT2. (a) Residual normal probability distribution and (b) relationship between the predicted value and actual value.
(a)(b)Figure13 shows the three-dimensional response curves of the influence of milling parameters on temperature T2. It can be clearly seen from Figure 12(a) that, within the range of study parameters, the effects of rotational speed and feed speed on milling temperature T2 in the cross mode are similar to those in the parallel cutting mode. The milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. The highest area of milling temperature is concentrated in the case of high rotating speed and low feeding speed. It can be seen from Figure 13(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 13(c) that, within the range of the study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 13
Response surface curves of the influence of milling parameters onT2. (a) Effect of rotation speed A and feed speed B on T2, (b) effect of rotation speed A and milling depth C on T2, and (c) effect of feed speed B and milling depth C on T2.
(a)(b)(c)Table8 shows the results of variance analysis of the prediction model of the quadratic regression equation of temperature T3 in the vertical cutting mode of cortical bone region by using ANOVA module. Figure 14(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T3, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 14(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Table 8
Analysis of variance for quadratic regression equationT3.
FactorsSum of squares of deviationsDegrees of freedomMean square errorFP (Prob > F)Model781.71986.86912.730.0351SignificantRotation speedA638.401638.406708.620.0106SignificantFeed speedB76.18176.18800.490.0091SignificantMilling depthC37.25137.25391.430.0263SignificantAB1.250E − 00311.250E − 0030.0130.9110AC0.1010.101.060.3266BC0.1010.101.060.3266A228.80128.80302.68<0.0001B23.4113.4135.840.0001C22.0512.0521.550.0009Residual0.95100.095Misfit term0.5450.111.330.3808Not significantPure error0.4150.082Figure 14
Residual analysis of prediction equationT3. (a) Residual normal probability distribution and (b) relationship between the predicted value and the actual value.
(a)(b)Figure15 shows the three-dimensional response curves of the influence of milling parameters on temperature T3. It can be clearly seen from Figure 15(a) that, within the range of study parameters, the milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. The highest area of milling temperature is concentrated in the case of high rotating speed and low feeding speed. It can be seen from Figure 15(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 15(c), within the range of the study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 15
Response surface curves of the influence of milling parameters onT3. (a) Effect of rotation speed A and feed speed B on T3, (b) effect of rotation speed A and milling depth C on T3, and (c) effect of feed speed B and milling depth C on T3.
(a)(b)(c)Table9 shows the results of variance analysis of the prediction model of the quadratic regression equation of milling temperature T4 in the cancellous bone region by using ANOVA module. Figure 16(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T4, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 16(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Table 9
Analysis of variance for quadratic regression equationT4.
FactorsDeviation sum of squaresDegrees of freedomMean square errorFP (Prob > F)Model796.83988.543370.220.0341SignificantRotation speedA643.201643.2024483.930.0086SignificantFeed speedB88.80188.803380.380.0243SignificantMilling depthC37.64137.641432.640.0101SignificantAB0.01110.0110.430.5276AC0.03110.0311.190.3010BC0.03110.0311.190.3010A226.82126.821020.780.0015B24.4914.49170.780.0002C22.1212.1280.560.6750Residual0.26100.026Misfit term0.09450.0190.560.7296Not significantPure error0.1750.034Figure 16
Residual analysis of prediction equationT4. (a) Residual normal probability distribution and (b) relationship between the predicted value and the actual value.
(a)(b)Figure17 shows the three-dimensional response curves of the influence of rotation speed and feed speed on milling parameters on T4. It can be clearly seen from Figure 17(a) that, within the range of study parameters, the milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. From the above, it can be concluded that the highest region of milling temperature is concentrated in the region with high rotational speed and low feed speed. It can be seen from Figure 17(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 17(c) that, within the range of study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 17
Response surface curves of the influence of milling parameters onT4. (a) Effect of rotation speed A and feed speed B on T4, (b) effect of rotation speed A and milling depth C on T4, and (c) effect of feed speed B and milling depth C on T4.
(a)(b)(c)According to the prediction model of the above cortical bone regions under different cutting modes, the anisotropy of bone material has a significant influence on the milling force and milling temperature, which cannot be ignored. Considering the influence of bone units in cortical bone structure, three cutting directions such as parallel to the bone unit, cross with the bone unit, and perpendicular to the bone unit should be used in the cortical bone area. Under the vertical cutting mode, the blade needs to accumulate more energy through the bone unit due to the enhancement of bone unit fiber compared to the other two cutting directions, so the milling temperature is higher than the other two cutting directions. Under the cross cutting mode, the milling crack grows along the bonding line, the matrix, and then the bone unit “peeling.” The milling force required for the bone material crack growth or fracture is larger than that in the parallel mode, and the heat generated is correspondingly less. The failure of bone tissue in the parallel cutting mode occurs mainly along the weak bonding line, so the milling temperature is minimal under the parallel cutting mode.


### 3.3. Global Sensitivity Analysis

The relative contribution rate of each parameter to the output parameter can be obtained by taking the model of milling force and milling temperature in different areas obtained by the response surface method as the objective function and calculating the global sensitivitySTi of the milling parameters of each variable. The first-order sensitivity index and global sensitivity index of each input milling parameter can be obtained, and the influence degree of input milling parameters on output milling force and milling temperature can be obtained by sorting them. The first-order sensitivity and global sensitivity of the prediction model obtained by sensitivity analysis are shown in Table 10. In order to more intuitively represent the relative contributions of different milling parameters in different areas in the above table about milling force and milling temperature, the first-order sensitivity and global sensitivity of the prediction model are expressed in Figure 12.Table 10
First-order sensitivity and global sensitivity of the prediction model.
Cutting modesCutting parametersT first-order sensitivityT global sensitivityParallel cutting mode of cortical bone regionRotation speedA0.660.68Feed speedB0.210.23Milling depthC0.050.07Cross cutting mode of cortical bone regionRotation speedA0.680.71Feed speedB0.190.22Milling depthC0.050.07Vertical cutting mode of cortical bone regionRotation speedA0.690.72Feed speedB0.210.23Milling depthC0.040.05Cancellous bone areaRotation speedA0.790.84Feed speedB0.110.12Milling depthC0.040.05As can be seen from Figure18, for a given parameter distribution, the global sensitivity of the milling temperature prediction model is not much different from that of the first-order sensitivity numerically, so the interaction between the parameters of the respective variables is not particularly obvious. This also verifies the results obtained from the analysis of variance for the prediction model above. It can also be seen that, for the milling temperature model established according to the experiment, the importance and influence of each input parameter on the milling temperature are basically similar under the consideration of different regions and different milling modes. The parameter that has the greatest influence on the milling temperature is the rotation speed, followed by the feed speed, and the influence of the milling depth is not obvious.Figure 18
First-order sensitivity and global sensitivity of the milling temperature prediction model in different regions.


## 3.1. Prediction Model of Milling Temperature

Based on the central composite experimental design method, the central composite in response surface model was selected for experiment design. The rotation speedA, feed speed B, and milling depth C are filled in Table 3, and 20 experiment design points were automatically analyzed and constructed by Design-Expert software. In the cortical bone area, milling experiments were carried out in three cutting directions: parallel cutting direction, cross cutting direction, and vertical cutting direction. The milling temperatures are T1 in parallel cutting mode, T2 in cross cutting mode, T3 in vertical cutting mode, and T4 in the cancellous area, respectively. The sequence of experiment points and experiment data results are shown in Table 1.Table 5
Range of test factors and coding table of quadratic regression design.
Range of variablexActual variablexStandardization ofzVariablex1x2x3z1z2z3Above levelx12x22x32111Zero levelx10x20x3000Under levelx11x21x31−1−1Distance fromr horizontalr×Δ1+x10r×Δ2+x20r×Δ3+x30rrrDistance from −r horizontalr×Δ1+x10r×Δ2+x20r×Δ3+x30−r−r−rStandard deviationΔ1Δ2Δ3Considering the influence of bone microstructure and milling trajectory, the empirical mathematical models of milling temperature were established in cortical and cancellous bone segments. The reliability was set as 95%, andP<0.05 was the level of statistical significance. The variance analysis of effective parameters was used to compare the differences of different parameters and combinations to verify the validity of different model items of milling parameters. The prediction model of the maximum milling temperatures T1, T2, and T3 under parallel cutting mode, cross cutting mode, and vertical cutting mode in the cortical bone region and the prediction model of milling temperature T4 in cancellous bone regions are shown as follows:(15)T1=16.605+0.021A−0.16B−17.452C−8.471×10−2AB+2.5×10−4AC+0.011BC−3.155×10−6A2+6.15B2+89.545C2,(16)T2=17.333+0.021A−0.137B−6.224C−5.555×10−7AB−1.5×10−3AC−0.022BC−3.123×10−6A2+5.32×10−4B2+77.727C2,(17)T3=19.2+0.021A+0.026B−0.138C−2.777×10−7AB−1.125×10−3AC−0.025BC−3.236×10−6A2+5.499×10−4B2+86.364C2,(18)T4=16.605+0.021A−0.16B−17.452C−8.471×10−2AB+2.5×10−4AC+0.011BC−3.155×10−6A2+6.15B2+89.545C2,where A represents rotation speed, B represents feed speed, and C represents milling depth.


## 3.2. Prediction Result Analysis

Milling temperatureT1 in the parallel cutting mode, temperature T2 in the cross cutting mode, milling temperature T3 in the vertical cutting mode were collected, respectively, and then milling temperature T4 was collected in cancellous bone region. Table 6 shows the results of variance analysis of quadratic regression equation prediction model of milling temperature under the parallel cutting mode of cortical bone area by ANOVA module, where the unimportant model items were eliminated. It can be seen from Table 6 that, in the influence factors of the model, rotation speed A and feed speed B are significant, milling depth C is not significant, and quadratic terms such as AB, AC, BC, A2, B2, and C2 are not significant. The misfit term in the table is used to indicate the degree of difference between the model and the experiment. In this fitting, P>0.05, the misfit term is not significant, which is beneficial to the model, and there is no misfit factors, indicating that the prediction model can replace the simulation point to analyze the simulation results. It can be seen from Table 6 that the fitting effect of quadratic multinomial model is significant and rotation speed A and feed speed B are the significant factors influencing milling temperature, and rotation speed A has the most significant influence.Table 6
Analysis of variance for quadratic regression equationT1.
FactorsSum of squares of deviationsDegrees of freedomMean square errorFP (Prob > F)Model803.78989.316001.180.0201SignificantRotation speedA649.641649.6443653.000.0098SignificantFeed speedB87.62187.625887.450.0054SignificantMilling depthC38.81138.812607.810.0136SignificantAB0.00010.0000.0001.0000AC5.000E − 00315.000E − 0030.340.5750BC0.02010.0201.340.2733A227.37127.371838.87<0.0001B24.2714.27286.64<0.0001C22.2112.21148.17<0.0001Residual0.15100.015Misfit term0.05550.0110.590.7090Not significantPure error0.09350.019Figure10(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T1, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 10(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Figure 10
Residual analysis of prediction equationT1. (a) Residual normal probability distribution and (b) relationship between the predicted value and the actual value.
(a)(b)Figure11 shows the three-dimensional response curves of the influence of milling parameters on temperature T1. It can be clearly seen from Figure 11(a) that, within the range of study parameters, the milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. It can be concluded that the highest region of milling temperature is concentrated in the region with high rotational speed and low feed speed. It can be seen from Figure 11(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 11(c) that, within the range of study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 11
Response surface curves of the influence of milling parameters onT1. (a) Effect of rotation speed A and feed rate B on T1, (b) effect of rotation speed A and milling depth C on T1, and (c) effect of feed speed B and milling depth C on T1.
(a)(b)(c)Table7 shows the results of analysis of variance (ANOVA) for the prediction model of milling temperature T2 in the cross cutting mode of cortical bone region, respectively.Table 7
Analysis of variance for quadratic regression equationT2.
FactorsSum of squares of deviationsDegrees of freedomMean square errorFP (Prob > F)Model798.09988.681137.670.0361SignificantRotation speedA649.641649.648334.490.0081SignificantFeed speedB79.52179.521020.250.0324SignificantMilling depthC40.80140.80523.490.0101SignificantAB5.000E − 00315.000E − 0030.0640.4933AC0.1810.182.310.2976BC0.08010.0801.030.5736A226.82126.82344.040.0004B23.1913.1940.940.0007C21.6611.6621.320.1472Residual0.78100.078Misfit term0.1750.0340.280.4084SignificantPure error0.6150.12Figure12(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T2, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 12(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Figure 12
Residual analysis of prediction equationT2. (a) Residual normal probability distribution and (b) relationship between the predicted value and actual value.
(a)(b)Figure13 shows the three-dimensional response curves of the influence of milling parameters on temperature T2. It can be clearly seen from Figure 12(a) that, within the range of study parameters, the effects of rotational speed and feed speed on milling temperature T2 in the cross mode are similar to those in the parallel cutting mode. The milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. The highest area of milling temperature is concentrated in the case of high rotating speed and low feeding speed. It can be seen from Figure 13(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 13(c) that, within the range of the study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 13
Response surface curves of the influence of milling parameters onT2. (a) Effect of rotation speed A and feed speed B on T2, (b) effect of rotation speed A and milling depth C on T2, and (c) effect of feed speed B and milling depth C on T2.
(a)(b)(c)Table8 shows the results of variance analysis of the prediction model of the quadratic regression equation of temperature T3 in the vertical cutting mode of cortical bone region by using ANOVA module. Figure 14(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T3, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 14(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Table 8
Analysis of variance for quadratic regression equationT3.
FactorsSum of squares of deviationsDegrees of freedomMean square errorFP (Prob > F)Model781.71986.86912.730.0351SignificantRotation speedA638.401638.406708.620.0106SignificantFeed speedB76.18176.18800.490.0091SignificantMilling depthC37.25137.25391.430.0263SignificantAB1.250E − 00311.250E − 0030.0130.9110AC0.1010.101.060.3266BC0.1010.101.060.3266A228.80128.80302.68<0.0001B23.4113.4135.840.0001C22.0512.0521.550.0009Residual0.95100.095Misfit term0.5450.111.330.3808Not significantPure error0.4150.082Figure 14
Residual analysis of prediction equationT3. (a) Residual normal probability distribution and (b) relationship between the predicted value and the actual value.
(a)(b)Figure15 shows the three-dimensional response curves of the influence of milling parameters on temperature T3. It can be clearly seen from Figure 15(a) that, within the range of study parameters, the milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. The highest area of milling temperature is concentrated in the case of high rotating speed and low feeding speed. It can be seen from Figure 15(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 15(c), within the range of the study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 15
Response surface curves of the influence of milling parameters onT3. (a) Effect of rotation speed A and feed speed B on T3, (b) effect of rotation speed A and milling depth C on T3, and (c) effect of feed speed B and milling depth C on T3.
(a)(b)(c)Table9 shows the results of variance analysis of the prediction model of the quadratic regression equation of milling temperature T4 in the cancellous bone region by using ANOVA module. Figure 16(a) shows the normal probability distribution of residual errors in the prediction model of bone milling temperature T4, which is most in line with the ideal situation. As can be seen from the figure, the points are all distributed near the straight line. Figure 16(b) shows the relationship between the predicted value and the actual value of the model. The closer the point is to the straight line, the smaller the difference between the actual simulation value and the predicted value of the model. According to the above analysis, the prediction models constructed are all reasonable.Table 9
Analysis of variance for quadratic regression equationT4.
FactorsDeviation sum of squaresDegrees of freedomMean square errorFP (Prob > F)Model796.83988.543370.220.0341SignificantRotation speedA643.201643.2024483.930.0086SignificantFeed speedB88.80188.803380.380.0243SignificantMilling depthC37.64137.641432.640.0101SignificantAB0.01110.0110.430.5276AC0.03110.0311.190.3010BC0.03110.0311.190.3010A226.82126.821020.780.0015B24.4914.49170.780.0002C22.1212.1280.560.6750Residual0.26100.026Misfit term0.09450.0190.560.7296Not significantPure error0.1750.034Figure 16
Residual analysis of prediction equationT4. (a) Residual normal probability distribution and (b) relationship between the predicted value and the actual value.
(a)(b)Figure17 shows the three-dimensional response curves of the influence of rotation speed and feed speed on milling parameters on T4. It can be clearly seen from Figure 17(a) that, within the range of study parameters, the milling temperature will increase with the increase of rotating speed and decrease with the increase of feeding speed. From the above, it can be concluded that the highest region of milling temperature is concentrated in the region with high rotational speed and low feed speed. It can be seen from Figure 17(b) that, within the range of study parameters, the milling temperature increases with the increase of rotating speed and milling depth. Therefore, when the milling temperature is large, both the rotation speed and the milling depth are large. It can be seen from Figure 17(c) that, within the range of study parameters, the milling temperature decreases with the increase of feed speed and increases with the increase of milling depth. Therefore, the larger area of milling temperature is concentrated when the feed speed is smaller and the milling depth is larger.Figure 17
Response surface curves of the influence of milling parameters onT4. (a) Effect of rotation speed A and feed speed B on T4, (b) effect of rotation speed A and milling depth C on T4, and (c) effect of feed speed B and milling depth C on T4.
(a)(b)(c)According to the prediction model of the above cortical bone regions under different cutting modes, the anisotropy of bone material has a significant influence on the milling force and milling temperature, which cannot be ignored. Considering the influence of bone units in cortical bone structure, three cutting directions such as parallel to the bone unit, cross with the bone unit, and perpendicular to the bone unit should be used in the cortical bone area. Under the vertical cutting mode, the blade needs to accumulate more energy through the bone unit due to the enhancement of bone unit fiber compared to the other two cutting directions, so the milling temperature is higher than the other two cutting directions. Under the cross cutting mode, the milling crack grows along the bonding line, the matrix, and then the bone unit “peeling.” The milling force required for the bone material crack growth or fracture is larger than that in the parallel mode, and the heat generated is correspondingly less. The failure of bone tissue in the parallel cutting mode occurs mainly along the weak bonding line, so the milling temperature is minimal under the parallel cutting mode.


## 3.3. Global Sensitivity Analysis

The relative contribution rate of each parameter to the output parameter can be obtained by taking the model of milling force and milling temperature in different areas obtained by the response surface method as the objective function and calculating the global sensitivitySTi of the milling parameters of each variable. The first-order sensitivity index and global sensitivity index of each input milling parameter can be obtained, and the influence degree of input milling parameters on output milling force and milling temperature can be obtained by sorting them. The first-order sensitivity and global sensitivity of the prediction model obtained by sensitivity analysis are shown in Table 10. In order to more intuitively represent the relative contributions of different milling parameters in different areas in the above table about milling force and milling temperature, the first-order sensitivity and global sensitivity of the prediction model are expressed in Figure 12.Table 10
First-order sensitivity and global sensitivity of the prediction model.
Cutting modesCutting parametersT first-order sensitivityT global sensitivityParallel cutting mode of cortical bone regionRotation speedA0.660.68Feed speedB0.210.23Milling depthC0.050.07Cross cutting mode of cortical bone regionRotation speedA0.680.71Feed speedB0.190.22Milling depthC0.050.07Vertical cutting mode of cortical bone regionRotation speedA0.690.72Feed speedB0.210.23Milling depthC0.040.05Cancellous bone areaRotation speedA0.790.84Feed speedB0.110.12Milling depthC0.040.05As can be seen from Figure18, for a given parameter distribution, the global sensitivity of the milling temperature prediction model is not much different from that of the first-order sensitivity numerically, so the interaction between the parameters of the respective variables is not particularly obvious. This also verifies the results obtained from the analysis of variance for the prediction model above. It can also be seen that, for the milling temperature model established according to the experiment, the importance and influence of each input parameter on the milling temperature are basically similar under the consideration of different regions and different milling modes. The parameter that has the greatest influence on the milling temperature is the rotation speed, followed by the feed speed, and the influence of the milling depth is not obvious.Figure 18
First-order sensitivity and global sensitivity of the milling temperature prediction model in different regions.


## 4. Conclusions and Limitations

(1)
Considering the influence of bone microstructure and milling trajectory, the prediction model considering the doctor’s surgical skills of the maximum milling temperature under parallel cutting mode, cross cutting mode, and vertical cutting mode in the cortical bone region and the prediction model of milling temperature in cancellous bone regions are established.(2)
By observing the influence curve of milling parameters on milling temperature using the ball milling cutter, it can be seen that, under the same cutting parameters, the milling temperature under the vertical cutting mode is the largest, under the cross cutting mode is the second, and under the parallel cutting mode is the smallest.(3)
It can be seen from the influence curve of milling parameters on milling temperature in cancellous bone region that the highest area of milling temperature is concentrated in the area with high rotating speed and low feeding speed, while the area with large milling temperature is concentrated in the area with large rotating speed and milling depth, small feeding speed, and large milling depth.(4)
By analyzing the parameter sensitivity of the temperature prediction model and the influence of the input parameters on the output milling temperature, it is concluded that rotation speed has the greatest influence on milling temperature, followed by feed speed and with the least influence for milling depth.(5)
In this experiment, doctors operated the spinal surgery cooperative robot for lamina milling, and the fresh pig femur was used instead of human bone, and the thermocouple was used to measure the milling temperature. Although multiple measurements were used to reduce the measurement error, there is still a certain gap between the measurement method and the actual milling temperature. In the future research, we will use human bone as the experiment sample, increase the sample size, and use the infrared thermal imager for temperature measurement, so as to improve the accuracy of the research.



---

*Source: 6806689-2021-10-14.xml*