Sensitivity analysis of temperature and force in robotic bone drilling process using Sobol statistical method

ABSTRACT The bone drilling process is indispensable in orthopaedic surgeries and treating bone breakages. It is also very important in dentistry and bone sampling operations. Bone is a very complex material and the process of drilling is very sensitive. Thus, bone drilling is one of the most important, common in the field of biomedical engineering. The bone drilling process can be promoted using automatic drilling machines and surgery-assisting robots. The problematic issue during operation is the high increase in drilling process temperature (higher than 47 °C) which leads to the so-called ‘thermal necrosis’ or cell death, and local burn in bone tissue. Furthermore, imposing higher forces to bone might yield to breaking or cracking, and consequently causes serious damages in bone. In this paper, the tool rotational speed, feed ratio and tool diameter were taken into account as process input parameters, and process temperature and thrust force were taken as output parameters. Design of experiments using response surface methodology was followed. Then, second linear governing equation was assigned to the model and its accuracy was evaluated. Later, Sobol statistical sensitivity analysis was used to ascertain the effect of process input parameters on process temperature and force. The results showed that among all effective input parameters, the tool rotational speed, feed rate and tool diameter have the highest influence respectively on process temperature and force. The behaviour of each output parameter with variation in each input parameter was further investigated.


Introduction
During bone breakage treatment, the goal is to help the broken bone relocate to its exact original place. Thus, broken bone segments should be placed exactly and accurately in a fixed location. Bolts are more often applied to hold broken parts solid. Knowing temperature and force during the bone drilling process is crucial for successful orthopaedic operations. Advancements in design of drilling tools and usage of automatic drills or surgery-assisting robots have already driven researchers. Optimisation of the bone drilling operation and avoidance of unacceptable errors and damages to the bone have been attempted [1]. Study and research on experimentation, analysis, modelling and simulation of the bone drilling process with automatic tools is remarkably ever-increasing [2]. Currently, using automatic drills and surgery-assisting robots, orthopaedic operations have been monumentally developed. Louredo et al. [3] achieved higher tool accuracy in layer removing and more accurate imposed force to the issue implementing a robotic system. Aziz et al. [4] increased the accuracy of tool positioning, control of imposed force and tool's forward/backward movements by introducing an algorithm. Boiadjiev et al. [5], for the purpose of analytical modelling of the process temperature when automatic drills were applied, investigated the mathematical temperature distribution modelling by thermal flow equation. They calculated the temperature of the tip of the drill bit and the specific coefficients for bones were taken into account [5]. Later, this study was further investigated by other researchers [6]. The force and temperature trend during bone drilling must be noted for desirable results [7]. Thermal necrosis and cell death will occur with temperature increase in bone [8]. A high temperature in bone drilling changes the state of the bone phosphates to alkaline. This leads to thermal necrosis and cell death, and causes the death of the bone tissue and decrease in material stiffness in the drilling affected zone [9]. This loosens the bolts applied in the operation [10] and also makes the recovery period lengthy, which is unacceptable. Thermal necrosis hinders the blood flow to the bone tissue and thus, it leads to cell death, local loss of bone tissue and weakening of its structure [11].
Thermal necrosis makes bolts loose, which leads to relative movements of the fixed parts and this is a failure in surgery [7].
The level of damage to the bone is connected to the increase in temperature and to the time of its exposure to heat [12]. According to the available literature, thermal necrosis is prone to happen for a wide range of temperature (44-100 C). When the temperature exceeds 70 C, thermal necrosis occurs immediately [13]. An increase in temperature from 47 to 50 C with exposure time of one minute influences bone tissue. However, when the temperature is less than 44 C, the thermal effect is not considerable if the exposure time is one minute at most [14]. Overall, most of the researches revealed that temperatures higher than 47 C within a minute cause thermal necrosis in bone tissue [15].
Presence of micro-cracks and injury to the bone tissue caused by high imposed forces aggravates surgery and hinders patient recovery [16]. Furthermore, the imposed force on the bone tissue corresponds to the temperature increase in the cortical bone [17]. The tool's rotational speed, feed rate, diameter and geometry highly affect the temperature and force evolution during bone drilling. Previously, many researchers have studied the effect of the tool's rotational speed and feed rate on the force and the temperature trend. Apart from the numerous investigations, studies on the process parameters to improve force and temperature reveal conflicts and contrasts [18]. Different studies about the effect of the cutting speed on the process force report conflicting results. According to Alam et al. [19], Basiaga et al. [20] and Jacob et al. [21], a raise in the tool's rotational speed reduces the force during bone drilling, whereas Lee et al. [22] interpreted that increasing the rotational speed boosts up the process force. Udiljak et al. found that tool rotational speed is not influential on axial force and it can be tagged as an ineffective parameter [23].
The effect of the tool's rotational speed and feed rate on temperature behaviour has also been contradictorily reported [13]. Vaughn et al. [24] stated that process temperature increases with a raise in rotational speed. Augustin et al. [25,26], Karaca et al. [27], Lee et al. [28], Udiljak et al. [23] and Pandey and Panda [29] reported that increasing the cutting speed and decreasing the feed rate leads to an increase in temperature. Matthews and Hirsch [30] observed that boosting up the rotational speed from 345 rpm to 2900 rpm did not have a discernible influence on process temperature when focused on human femur bone drilling process, and Sharawy et al. [31], that increasing the rotational speed from 1225 to 2500 rpm reduced the process temperature. Moreover, Shakouri et al. [32] concluded that drilling with high rotational speeds reduces the process temperature. Augustin et al. [25] reported that the maximum machining temperature is subsided with an increasing feed rate and, according to Pandey and Panda [18], a reduction in the feed rate decreases the process temperature. Alam [33] found that the process temperature with a feed rate of 20 mm/min was lower than that with a feed rate of 50 mm/min.
Considering the reported literature, up to now, statistical sensitivity analysis methods have not been implemented to quantitatively study the accurate effect of different parameters on temperature and force during bone drilling. Sensitivity analysis studies the uncertainty in the output of a given model and defines how the uncertainty in the output of a model can be related to the uncertainty in the inputs of the model [34]. This method is used to specify the effective and the ineffective parameters on the output. Sensitivity analysis models are classified into two local and general groups [35]. The E-Fast method was introduced by Cukier et al. [36] and later, Saltelli et al. [37] improved it.
In this paper, considering the tool rotational speed, feed rate and tool diameter as input parameters and force and temperature as output parameters of the bone drilling process, first design of experiments was planned using response surface methodology. Then, a secondorder regression governing equation was extracted for each output parameter and the accuracy of the model was evaluated. After reviewing the different sensitivity of the methods of analysis, Sobol statistical sensitivity analysis method was chosen and the effect of the input parameters on each output parametertemperature and forcewas scrutinized.

Response surface method
Response surface method is a mathematical-statistical method used to model and analyse problems as complex functions of some variables. The goal of 'response surface methodology' (RSM) is to statistically model and optimize the problem [38]. The basics of the RSM are design of experiments and statistical optimization. The design of experiments is a suitable tool for engineers in developing experiments with less time and expense. Applying this method requires less process time and costs [39]. Evaluation of the accuracy of experiments, governing the mathematical model of the experiments, developing interaction diagrams of input variables, experiments optimization and assuring of the exact reliance of the developed model are some of the advantages of RSM [40]. Considering the factors and effective interactions, the general form of the equation is Equation (1) [39]: Furthermore, RSM is able to model the relation between the inputs and outputs and represent it as a secondorder linear integration equation [41].

Sensitivity analysis methods
Sensitivity analysis is a tool to study systems in order to reveal the effect of input parameters on output parameters. Sensitivity analysis is categorized from its functionality to deterministic and probabilistic analyses. If a model form is focused, sensitivity analysis can be classified into graphic, mathematical and statistical. Graphic method: in this method, sensitivity is depicted in the form of diagrams, tables or surfaces. The graphic method is used to show output variation under the effect of input parameters.
Mathematical method: in this method, sensitivity can be found from output variation based on input variation. This method employs formulations and computations to study the output with a small change in inputs.
Statistical method: This sensitivity analysis method tries to simulate the system with probabilistic distribution of input variables. Then, the effect of this distribution on output variables is analysed. In this method, the interactive effect of several input variables on the behaviour of an output variable can be studied. The Sobol method is a statistical sensitivity analysis which is independent of the model and it is based on variance degradation. This method can be used for nonlinear and nonuniform functions [42].
In this method, for the defined model, Y = f(X), where Y is the output of the model, X (x 1 ; x 2 ; . . . ; x n ) is the input parameter vector and the output variance (V) is the summation of the variances of each degraded term (Equation (2)): where V i is the first effect for each input factor ...;n being the interaction of n factors. Sensitivity factor is the ratio of each order variance to the total variance (S i ¼ V i V is the first-order sensitivity factor and S ij ¼ V ij V is the second-order sensitivity factor). The total sensitivity factor which equals the overall effect of the parameter is the summation of all orders of the sensitivity factors of these parameters, which can be computed using Equation (3):

Experimental design
The input parameters in bone drilling analysis are rotational speed of the tool (V), feed rate (f) and tool diameter (D). Crucial output variables are the maximum process temperature (T) and maximum thrust force (F). An automatic drilling machine was used, where the drill bit used in the experiments was made of 'high-speed steel' (HSS). The holes' depth in the analysis was 8 mm.
The thrust force was measured with a dynamometer. Ktype thermocouples were implemented to measure the temperature at a depth of 3 mm and a distance of 0.5 mm from the hole wall [43]. Figure 1 illustrates the typical experimental set-up. Bovine femur cortical bone, which is similar to the human cortical bone, was used. To have experiments more similar to a real surgery, it was ensured that the bone tissue was alive a few hours before experimentation. The location of the thermocouples is shown in Figure 2.

Mathematical modelling
Rotational speed, feed rate and tool diameter were three input variables and 3 3 full factorial experiments were done. RSM and central composite design were employed to develop the model. In Table 1, the input variables and   Table 2. Minitab v.16 software package was used to analyse the results and also to obtain the coefficients of the governing empirical equation. Using RSM and data analysis, a second-order linear regression equation was derived for each output variable based on the input variables. The results are further discussed and model optimization was also performed.

Data analysis and process modelling interpretation
Analysis of variance (ANOVA) based on temperature and force data analysis is presented in Tables 3 and 4. Assuming average reliability of 95% in a precise engineering experimentation, a P-value less than 0.05 is considered to ascertain the effectiveness of different model terms [39].
By the second-order linear model, the prediction error sum of squares (PRESS) value for temperature is 136.470 and for force it is 435.473. The second-order linear regression equations governing the temperature and force behaviours are according to Equations (4) and (5), respectively.
Force ¼ 145:607 À 0:0432796V þ 0:150566f The values for temperature are R 2 = 96.31%, R 2 (pred) = 90.12% and R 2 (adj) = 94.35% and for force they are R 2 =     Figure 3, it can be inferred that the accuracy of the developed models for temperature and force is acceptable.

Parameters sensitivity analysis
In this section, we tried to clarify the effect of the tool's rotational speed, feed rate and tool diameter on the process temperature and force in the bone drilling process.
The plots in Figures 4 and 5 show the behaviour of the temperature and the force with different input variables using a developed model with the RSM method. Then, the scatter points which are temperature and force   response with input parameters were extracted using SimLab software based on the Sobol method. In this method, opposite to graphic methods in which all inputs but one is constant, using a special algorithm, all parameters were varied at the same time. For details, see [42].

Analysis and investigation of the effect of input variables on bone drilling temperature
Analysis of drill rotational speed. The interaction of feed rate and rotational speed on temperature behaviour is illustrated in Figure 4 for different tool diameters. This plot shows a second-order surface with a saddle point. This saddle point possesses a maximum regarding the rotational speed and a minimum regarding the feed rate. As a result, the behaviour of the temperature near this point is different. Some contrary results [27,29,31,33] reported previously might be attributed to the temperature variation in the vicinity of this saddle point. For example, at a constant feed rate, the temperature initially increases with the raise in rotational speed and after this saddle point it decreases. Therefore, both contrary reports of better conditions of machining at low and high rotational speeds can be true. Noticing this critical point is a must which was not adequately considered due to lack of precise modelling. Considering the reduction in the slope of the rotational speed diagram and also the saddle-like area in the feed rate and rotational speed interactive diagram, it can be predicted that high speed machining causes a lower rise in temperature and thus, lower probability for the presence of thermal necrosis.
According to Figures 4 and 6, it can be observed that with an increase in rotational speed, the maximum temperature increases. Therefore, the least damage to the bone, from thermal necrosis view point, can be achieved at low rotational speeds. Since a saddle point is present (Figure 4), it is predicted that, in high speed machining (after saddle point), temperatures less than 47 C are achievable. As shown in Figure 6, the slope of the rotational speed-temperature diagram decreases gradually. This is attributed to easy exit of chip from the hole. This matter was studied recently and current results show that high-speed drilling decreases the possibility of thermal necrosis and the rate of temperature increase [32].
Seldom are tools with diameter of 6 mm and higher used in surgeries. Therefore, in this study, similar to some reports, we chose to use more common tools with diameters of 2.5-5 mm [6]. In preliminary experiments, our attempts to use tools with a diameter of 6, 8 and 10 mm in the bone drilling process showed that the process temperature rises very quickly with the tools' diameter (data not shown). It seemed that thermal necrosis was inevitable in operations with tool diameter equal to or greater than 8 mm. As can be seen from the above diagrams, with rotational speed of 2500 rpm and tool diameter in the range of 6-10 mm, the process temperature becomes critical (third diagram) and even for lower rotational speeds with tool diameter of 8-10 mm, thermal necrosis is probable (Figure 7).
Effect of drill diameter. Considering the F-value and the coefficient of the regression equation, the tool diameter in the range of or experiments was very influential on the generated heat and the development of thermal necrosis. It can be implied from the plots of Figures 4  and 8 that with the increase in tool diameter, the process temperature increases.
As it can be seen in the plots of Figures 4 and 8, an increase in drill diameter leads to a rise in the process maximum temperature. With an increase in tool diameter, the contact surface between the bone and the tool increases, which in turn raises the frictional forces and consequently, the generated heat increases. Moreover, increasing the tool diameter causes higher drilling forces to be imposed on the bone and the drill itself. This leads to an increase in the temperature (Figures 4 and 8).
With increasing the tool diameter, there is a remarkable rise in the rate of temperature growth. As can be observed in the plots of Figure 4, increasing the tool diameter leads to an increase in the temperature. The allowable range with no possibility for thermal necrosis   Figure 7. Effective interaction diagrams of feed rate and diameter in experimental study with 6-10 mm tool diameters. for a tool diameter of 5 mm is much lower than that for tool diameter of 4 or 2.5 mm. It is noteworthy from the biomedical point of view that a smaller tool diameter corresponds to a shorter recovery period [15].
Feed rate. The effect of feed rate on process temperature behaviour is complex. In lower feed rates, the applied force on the bone friction of tool-bone and thickness of the deformed chip are initially lower, so the chip easily finds its way out of the hole. This leads to less generated heat and lower process temperature. Later, the time of tool-bone contact increases, which leads to an increase in the heat transferred from the tool to the bone. This consequently causes a rise in bone temperature. Then, the applied force to the bone considerably increases and finally, the non-deformed chip becomes thicker, while the tool-bone contact duration decreases. Depending on the situation, it is possible that the minimum generated heat occurs at lower, or higher, or even intermediate values in the feed rate range. Therefore, the effect of feed rate on process temperature depends on the experimental conditions and the range of the variables.
As it can be seen from the effect of feed rate (Figure 9) within the presented range of experiments, the minimum rate of temperature increase occurs in the middle of the feed rate range. In this state, two improvements are obtained instantaneously: a) the tool-bone contact time is shorter than that with lower feed rates and b) the imposed and frictional forces between the tool and the bone are less than those at higher feed rates. These benefits are observable in Figure 4. The figure shows that with increasing the feed rate at a constant rotational speed, at first, the maximum temperature decreases, and then it increases again. Initially, increasing the feed rate facilitates the escape of the bone's brittle chip and develops less friction due to shorter bone-tool contact time but then, due to the thicker non-deformed chip and also a rise in the imposed and the frictional forces, the process temperature increases.

Analysis and investigation of the effect of input variables on bone drilling force
In this section, knowing the interactive effect of rotational speed and feed rate at different tool diameter ( Figure 5) and also diagrams related to major effective factors (Figures 10-12), the influences on the force during the drilling process are studied.
Analysis of drill rotational speed. As can be seen from Figures 5 and 10, increasing the tool rotational speed reduces the imposed force to the bone. Moreover, a raise of the rotational speed eases the movement of the chips out of the hole and avoids blockage of chips on the wall, which in return leads to a decrease in the friction coefficient between the hole wall and the drilling tool. This yields a reduction in process force. Considering the saddle-like shape of the process temperature model, one can infer that in case high speed machining becomes practical in bone surgery, high rotational speed is beneficial for both the process temperature and the force.
Analysis of feed rate. Lower feed rates produce lower imposed forces to the bone. The minimum imposed force to the bone occurs with higher rotational speeds and lower feed rates. Increasing feed rate gives a rise in the thickness of the formed chip and consequently, the imposed force to the bone boosts up. With lower feed rates and less tool penetration into the bone, tissue deformation of the chip and its flee out of the hole become easier. Moreover, the formed chip is thinner and consequently lower force is imposed to the bone. This point can be clearly observed in Figures 5 and 11. Higher process force will increase the possibility of bone breakage or damage. However, it should be noted that, in orthopaedic surgery, the duration of operation is very important and when using surgery assisting robots, operation duration is directly proportional to the feed rate.
Effect of drill diameter. It can be implied from Figures 5 and 12 that with an increase in tool diameter, the contact area between the tip of the tool and the bone extends with the second order of the tool diameter and hence, the generated forces increase with tool diameter. Inappropriate application of tools with big diameters in orthopaedic operations increases the temperature and force in bone tissue and results in a lengthy recovery period. Using a micro-drilling process or drills with small diameter in bone surgery is crucial in the performance of operation.

Sobol analysis in sensitivity investigation of effective parameters on process temperature and force
Regarding the available results shown in Figures 6, 8, 9 and 10-12, sensitivity analysis for temperature and force in bone drilling process, it was observed that in the considered ranges for input parameters, the tool's rotational speed, feed rate and tool diameter have the most pronounced influence. The slope of the variation curve for each item dictates its effectiveness. Figure 13 shows sensitivity analysis using the Sobol method and it confirms the previous results. Among the input parameters, noticing computations using the Sobol method and SimLab software that with a special algorithm varies all input parameters at the same time to detect their effect on output variable [42], the tool's rotational speed with 61% is the parameter that most strongly affects the process temperature. This is followed by feed rate and tool diameter, with 26% and 13%, respectively. Moreover, regarding the process force, the tool's rotational speed, feed rate and tool diameter account for 59%, 24% and 17% of the observed effect, respectively.
The Sobol method is an accurate method but it involves high computation costs. A special algorithm that provides the possibility to change all input parameters simultaneously is used. This feature is different from other graphic methods in which all inputs but one is   Figure 11. Effect of feed rate in force.
constant. The algorithm used in the Sobol method can find how each input parameter and their interactions can affect the output variable more accurately.
The advantage of the Sobol method over ANOVA is the ability to quantitatively compute the exact effect of input parameters. Here, effective and ineffective parameters are clearly distinguished. However, in the analysis of interactive influence of parameters and complex relation, the response surface method yields better results.

Conclusions
In this paper, in addition to developing a model for bone drilling process using response surface methodology, where tool rotational speed, feed rate and tool diameter were included as input parameters and process temperature and force as output parameters, mathematical second-order linear regression equations for the process temperature and the force were derived. The accuracy of modelling as well as the importance and effect of each of the input parameters on the process force and temperature were studied by Sobol sensitivity analysis. The results from this study demonstrated that the interactive diagram of tool rotational speed and feed rate yields a second-order surface with a saddle point. The presence of this saddle point explains the conflicting results given by other researchers about temperature behaviour. With the increase in tool diameter, temperature increases. In the selected range of input parameters, the temperature increases with the rotational speed. The effect of the feed rate is rather complicated: the tool-bone contact time is lower than that with lower feed rates, and the applied and frictional forces between the tool and the bone are lower than at higher feed rates. By increasing the rotational speed and decreasing two parameters, 'feed rate' and 'tool diameter', the force imposed onto the bone will decrease. Considering a simultaneous change of all input parameters, the tool's rotational speed, feed rate and tool diameter are the most influential because, in the variation range for each of these parameters, the slope of the variation curve dictates their effectiveness. Generally, the parameters with greatest effect on process temperature are the tool's rotational speed, feed rate and tool diameter with 61%, 26% and 13%, respectively. Moreover, effective parameter on process force are tool's rotational speed, feed rate and tool diameter with 59%, 24% and 17%, respectively.

Disclosure statement
No potential conflict of interest was reported by the authors.