Optimization of bone drilling parameters using Taguchi method based on finite element analysis

Thermal necrosis results fracture problems and implant failure if temperature exceeds 47 °C for one minute during bone drilling. To solve this problem, this work studied a new thermal model by using three drilling parameters: drill diameter, feed rate and spindle speed. Effects of those parameters to heat generation were studied. The drill diameters were 4 mm, 6 mm and 6 mm; the feed rates were 80 mm/min, 100 mm/min and 120 mm/min whereas the spindle speeds were 400 rpm, 500 rpm and 600 rpm then an optimization was done by Taguchi method to which combination parameter can be used to prevent thermal necrosis during bone drilling. The results showed that all the combination of parameters produce confidence results which were below 47 °C and finite element analysis combined with Taguchi method can be used for predicting temperature generation and optimizing bone drilling parameters prior to clinical bone drilling. All of the combination parameters can be used for surgeon to achieve sustainable orthopaedic surgery.


Introduction
In orthopaedics surgery, bone drilling is widely known to create hole for insertion of screw [1]. Before the insertion of screw, drilling then tapping are conducted [2]. About 95% of post-trauma therapies need drilling holes that are needed to install the screws for plate attachments or other prosthetics or for fixing and correcting fractures of bone [3]. Study on bone drilling has retrieved many attentions. One significant aspect of the study concerns the generation of heat during the drilling process [4]. A few surgical methods have been explained in terms of influences during osteotomy of implant site [5].
The heat generation in bone drilling could be partially dissipated by the fluids of tissue or blood and partially involved by the chips created [6]. Cutting action between drill bit and bone results in significant heat because of removal of material and friction coefficient between the drill bit and bone [7]. Drilling actions when preparation of site of dental implant could cause not only physical failure to the bone but also a rise of temperature in surrounding the implant site [8]. Temperatures over 47 °C leads to thermal osteonecrosis providing to loosening of screw and subsequently refractures and failures of implant [9], reduces protein regeneration capacity [10] and irreversible bone death [11]. Bone drilling factors that affect the temperature change when the preparation of implant site including drill bit geometry, applied force, spindle speed [12], drill wear, the effect of saline temperature applied for irrigation and guide of drill were studied in several investigations resulting substantial importance [13].
Temperature measurement and distribution at the interface of drill bit and bone is hard. Engineers and scientists have been depeloving numerical method for heat transfer analysis. Finite Element Analysis (FEA) shows effective proficiency in all engineering areas including in medical science. Bone In 1999, Davidson performed a parametric analysis to study the thermal impact on bone drilling. Parameter used were drill diameter, spindle speed, feed rate, drill helix and point angle, density, bone thermal conductivity and specific heat. Highest temperature increased by feed rate from 0.45 to 1.8 mm/s then decreased slightly by feed rate of 4.5 mm/s and lowest temperature occured at helix angle of 38 o and point angle of 130 o [14]. In 2004, Li et al. compared surface heating and volumetric sources in the laser melting modeling for ceramic materials concluding that the simulation result with volumetric heating source was accurate and stable [15]. In 2011, Lee et al. worked out a new thermal model for bone drilling by a sensitivity analysis and single parametric study [16] followed by Mokhtar and Fawad in 2012 which was to accomplish theory of bone drilling modeling to gain the heat flux [17].
In 2015, Gok et al. studied the drilling and thrust power including coefficients of heat transfer for bone drilling numerical simulations using a K-wire. DEFORM-3D software was used to conduct the bone drilling process resulting there was a great consistency between FEA and experimental results. This has shown the finite element capability [18].
Other finite element method functions are to investigate drilling mechanism and to predict temperature of drilling. In 2009, Tu et al. applied a three-dimensional finite element model showing that temperature decrease was caused by decreasing initial temperature of cutting tool and the thermal influenced zone size [19]. In 2013, Tu et al used an elastic-plastic finite element model for analyzing the behaviour of thermal contact between a drill bit and bone. The results showed the distribution of and rise of temperature rise surrounding the hole drilled could be successfully assessed [20]. In 2012, Sezek et al. analysed change of temperature during drilling of bone by application of Finite Elements Method (FEM) by deciding the points drilling in bone and found a secure zone or below temperature of 45 o C [21]. The numerical method in bone drilling simulation could be a foundation for surgical tools optimization including for determination of their properties [22] and provide tools for development of surgical operations [2].
Although many numerical simulations have been conducted by proposing new models, optimization of finite element analysis in bone drilling is still lacking. This paper presents an elaborated study from fundamental to optimization of three-dimensional model. A method is proposed to develop a thermal model to predict rise of temperature during bone drilling. In this study, effects of drill diameter, feed rate and spindle speed on generation of temperature in the bone drilling simulation were investigated then optimized using the Taguchi experimental design method. To optimize the effects of the process parameters, signal-to-noise (S/N) ratio and orthogonal array were calculated.

Heat generation
Coupled temperature-displacement method is applied to calculate and solve concurrently temperature and stress/displacement problems. A coupled analysis is applied as thermal and mechanical solutions influence strongly one another. Abaqus/Explicit software allows coupled temperature-displacement analysis modes. Fully coupled thermal-stress analysis in the software is always transient. In the software, equations of heat transfer are calculated by integration of explicit forward-difference time below [23].
as the temperature generated at node N whereas the subscript i as the number of increment in a step of explicit dynamic. The current temperatures are referred by values of ̇ from the foregoing increment. The values of ̇a re calculated at the start of the increment as follow [23]. Where ( ) as the lumped capacitance matrix, as the vector of applied nodal source and as the vector of internal flux. The explicit method integrates by means of time by applying a lot of small time increments. The forward-difference and central-difference operators are stable. The limit of stability for the two operators is obtained by [23]: Where as the maximum frequency in the equations system of the response of mechanical solution and as the largest eigenvalue in the equations system of the response of thermal solution. An approximation to the limit of stability for the forward-difference operator in the response of thermal solution is calculated by [23]:

Where
, as the minimum dimension of element in the mesh and = as the thermal diffusivity of the material. The letters k, and c are the thermal conductivity, density, and specific heat of material, serially. Boundary conditions may be applied to determine displacements/rotations and temperatures at nodes in fully coupled thermal-stress analysis. Boundary conditions at a step of dynamic coupled temperaturedisplacement response should apply suitable references of amplitude if boundary conditions are determined for the step without references of amplitude, they are applied instantly at the step start. Due to the software does not recognize jumps in displacement, a boundary condition of non-zero displacement value that is determined without reference of amplitude will be neglected and a boundary condition of zero velocity will be implemented [23]. The heat generated in the primary zone of shear deformation is calculated by [24]: Where is the force of shear along the plane of shear, is the velocity of shear, is the yield stress of shear of bone, and As is the area of shear. They all are functions of drill bit radius r.

Thermal model and simulation
A three-dimensional dynamic coupled temperature-displacement finite element analysis was used to simulate cortical bone drilling process. Simulations were carried out using Abaqus 6.14 software (Dassault Systemes Simulia Corp, USA). The bone drilling process was modeled in Abaqus/Explicit mode. Due to the simulation was dynamic, damage failure criteria was applied to manage the removal of element. Drill bit and cortical bone were modelled. The cortical bone was modeled as a rectangular work piece by a length of 30 mm, a width of 30 mm and a thickness of 6 mm. The drill bit was modeled with a helix angle of 23° and a point angle of 120°. Since drill bit stiffness is much higher than that of the bone and in order to reduce the analysis time, rigid body was assigned to the drill bit. For appropriate prediction of temperature during drilling, fine mesh was used at and in the surrounding of the cortical bone to be drilled, whereas coarse mesh was used to discretize the cortical bone away from the drill bit path. The cortical bone was meshed with C3D8T 8-node 110070 elements consisting of 61830 hexahedral elements and 48240 wedge elements whereas the drill bit was meshed with 6227 C3D4 4node tetrahedron elements. Due to cortical bone has relatively small mechanical anisotropy, isotropic equivalent homogeneous material was applied instead of anisotropic and heterogeneous structure [25]. The interest region was the surrounding of the drilled hole where the temperature was highest. The  To guarantee stability and convergence of the simulation, a time period step of 0.0005 s was used. Thermal properties of the drill bit and the bone are listed in Table 1. Boundary condition was fixed on both vertical faces of the bone and for the drill bit was set free for rotation along the bone axis and vertical movement. Initial temperatures for the bone and the drill bit were set at 23 o C following room temperature. These initial and boundary temperatures were the same for all simulations. Contact interaction type of the drill bit and the bone was surface-to-surface contact with friction coefficient of 0.7. First surface (master) was assigned to the drill bit, whereas second surface (slave) was assigned to the cortical bone. Johnson-cook damage model was used for the chip -workpiece removal in the drilling simulations. [31].
Where ̅ is the equivalent flow stress, A is the material initial yield strength, the parameters B, n, C (as listed in table 2) and m are material model parameters, − is the equivalent plastic strain and ̇− is the rate of equivalent plastic strain which is normalized with a rate of reference strain ̇0 and ̂ is the homologus temperature calculated in equation 7 [31]. Where is the workpiece instantaneous temperature, is the room temperature and is the material melting temperature and is the material melting temperature. Element failure has disadvantage which is material deletion. The main disadvantage with element failure is the deletion of material. It is mostly non-physical that mass is deleted from the drilling process and this deletion affects the forces (pressure) between the drill bit and the bone. Therefore, it affects all simulation result. To minimize this effect, the density of mesh must be very fine [31].

Taguchi method
The Taguchi method is the most significant technique for effective and systematic approach for better optimization solution proposed by Taguchi (1990) [33] and to provide a significant reduction in experiment number for same number of process control factors [34]. Smaller-the-better criterion was used to calculate signal-to-noise (S/N) ratio where T is temperature generated, then response was calculated using the obtained S/N ratio.

S/N Ratio= -10 log T 2
In Taguchi method, for conducting the experiment, factors and levels were defined. Factors are controllable parameters which affect the experiment significantly. Three principal drilling parameters were chosen for numerical calculations same as the situation of real drilling processes. They were drill bit diameter, feed rate and spindle speed,. The bone temperature was the target during the drilling process, so the maximum temperature of the entire bone was extracted. The drill diameter, feed rate and spindle speed were the controllable parameters for this experiment since they can influence temperature generated during bone drilling. The controllable parameters along with the levels used are shown in Table 3.  After defining the factors and their levels, orthogonal array was defined. To select an appropriate orthogonal array, a degree of freedom (DOF) was calculated. DOF is needed to determine which level is better and is defined as the number of comparisons. DOF was 1 for mean value and 8 for two each for the remaining factors, so the total DOF was 9. Therefore L9 orthogonal array or a number of nine experiments were selected for this study as shown in Table 4. Figure 2 shows heat generation during drilling simulation. Maximum temperature resulted from each experiment was used for the optimization. For calculating S/N ratio, the objective function of the temperature generated was smaller-the-better criterion due to the lower the temperature, the better the bone drilling process. The S/N ratios of all the experiments were calculated and tabulated as shown in Table 5. The average of all the factors at each level or response was then calculated from the obtained S/N ratio as shown in Table 6.      This section discusses the result obtained by using the methodology explained in the previous section. The variation in temperature generated are shown in Figure 3. Temperature increases as the feed rate increases. The feed rate of 120 mm/min generates highest temperature followed by the feed rate of 100 mm/min and 80 mm/min. The change effect in the spindle speed for the drill bit of 4 mm indicates a significant increase in the temperature generated due to with the increasing of the spindle speed, bone chips may erode the drilled surface with more centrifugal force. The change effect in the diameter indicates a decrease in the temperature generated due to by increasing the diameter, contact area between the drill bit and the bone is higher with lesser pressure. Temperature decreases when the diameter changes from the diameter of 4 mm to 6 mm then slightly increases when the diameter changes from the diameter of 6 mm to 8 mm. However, the results show that the all temperatures generated are below 47 o C.

Results and discussion
According to the calculation results from the Taguchi method that are shown in Table 6, the combination of drilling parameters that produces the best temperature generated is by the drill diameter of 6 mm, the feed rate of 80 mm/min and the spindle speed of 400 rpm. However all combination can be used for drilling process since those parameters produce temperatures which are below 47 o C.

Conclusions
The present paper used finite element analysis for temperature generation then the results were optimized by Taguchi method. The above results and discussions are summarized as follows: -From the three linear graphs, it is clear that the optimum values of the factors and their levels are the diameter of 6 mm, the feed rate of 80 mm/min and the spindle speed of 400 rpm. -The results show that the drill diameter, the feed rate and the spindle speed give significant effects on temperature generation in bone drilling. It is found that the feed rate has the highest influence on temperature generated followed by the drill diameter then the spindle speed respectively. -All the factors and levels used can be reference for parameters used by surgeon in bone drilling to prevent thermal necrosis. -Combination of finite element analysis with Taguchi method is suitable for analysis prior to clinical bone drilling since experimental bone drilling needs complex tools and preparation.