Thermo-mechanical stresses distribution on bone drilling: Numerical and experimental procedures

In bone drilling, the temperature and the level of stresses at the bone tissue are function of the drilling parameters. If certain thresholds are exceeded, irreversible damages may occur on the bone tissue. One of the main challenges in the drilling process is to control the associated parameters and even more important, to avoid the surrounding tissue damage. In this study, a dynamic numerical model is developed to determine the thermo-mechanical stresses generated during the bone drilling, using the finite element method. The numerical model incorporates the geometric and dynamic characteristics involved in the drilling processes, as well the developed temperature inside the material. The numerical analysis has been validated by experimental tests using polyurethane foam materials with similar mechanical properties to the human bone. Results suggest that a drill bit with lower drill speed and higher feed rate can reduce the strains and stresses in bone during the drilling process. The proposed numerical model reflected adequately the experimental results and could be useful in determination of optimal drilling conditions that minimise the bone injuries.


Introduction
The drilling process is one of the most frequent mechanical procedures in machining operations.The industrial concepts of productivity and surface integrity in material removal processes can be translated to medical applications.3][4][5][6][7] Currently, it is well known that bone drilling process is a typical coupled thermo-mechanical problem.During the drilling, the heat generation occurs as a result of plastic deformation and friction along the cutting tool and the workpiece interface. 8,91][12] Besides the negative thermal impact, high drill speed with higher cutting forces and tool vibrations can also cause damages on the bone microstructure, resulting in microcracks formation and bone fracture.The mechanical damage to bone cells could delay the healing process after the surgery, reducing the fixation strength and in some cases even the implant failure. 9,13,14Piattelli et al. 15 reported eight cases of failed implants due to suspected thermally induced bone necrosis.Recently, Augustin et al. 2 reported that the implant failure rate for leg osteosynthesis ranges between 2.1% and 7.1% and it is higher than the failure rate of upper limbs due to physiological pressure during locomotion.
The growing increase of bone drilling as a medical intervention and the importance of the associated problems have motivated the development of methodologies to find out effects of drilling parameters that minimise the bone injuries.However, there is a still lack of information with regard to the strain and thermal stresses distribution in bone tissue during drilling.Although several methodologies have been proposed for estimated values for bone temperature and cutting forces, none of them includes the thermal stresses distributions on the bone.The large number of significant drilling parameters makes the study of bone drilling, through experiments, impracticable.The development of accurate numerical models is necessary in order to reproduce the factors that influence the output variables (mainly temperature, drilling forces and surface integrity).In the literature, only a few numerical simulations have attempted to model the bone drilling.Recently, Marcos et al. 1 have made a bibliographic review of the main contributions on the field of numerical modelling of bone cutting, including the bone drilling effect.They found that the numerical models are still far from clinical applications and there is a clear requirement for improvement in the models available for drilling simulations.Most of the models presented in the literature were performed with limited parameters and use two-dimensional (2D) simulations; furthermore, the analysis was not validated comparatively with experimental tests. 6,13,16Experimentally, several techniques have been developed to measure the strain level in body surfaces during the drilling processes for many industrial applications, but have not been widely used for medical situations.Usually, the strain level is measured using strain gauges and a data acquisition system from which the stresses can be calculated.
The present study was designed for the evaluation of thermo-mechanical stresses generated during the bone drilling, as function of different drilling parameters (drill speed, feed rate, drill diameter and hole depth) and temperature changes during the process.A three-dimensional (3D) dynamic finite element (FE) model with the element removal scheme was developed to simulate the thermo-mechanical behaviour of the contact region between the drill bit and the bone tissue.The simulations were performed using the LS-DYNA program, as a nonlinear transient dynamic finite element.This program incorporates the dynamic characteristics involved in the process with the accurate geometrical considerations, as well the developed temperature inside the material.Physical problems involving severe loading conditions applied in a short time are best simulated with explicit dynamics solutions, such as problems with nonlinearities, changing contact, and failure or material separation.In addition, an experimental methodology was developed using solid rigid polyurethane foam materials (Sawbones, Pacific Research Laboratories, Inc., Vashon.WA, USA) as an alternative for human bone.The solid polyurethane foam of 0.80 g/cm 3 density was selected for the cortical bone, which was approved by the American Society for Testing and Materials for a standard material for testing orthopaedic devices and instruments.The foams were instrumented with strain gauges to measure the level of strains on the surface during the drilling.Simultaneous, the temperature distribution on the foams surface and in the drill bit were measured with the thermography equipment.Inside the foams' material, at different positions from the drilling area, thermocouples were used to measure the temperature.The results of the numerical simulations are discussed and compared with the experimental results.
The used methodologies in this study are an important contribution in the bone drilling modelling.The incorporation of temperature data as a prescribed temperature on the FE model provides not only the mechanical behaviour of the bone, but also the thermal-mechanical behaviour with different feed rates and drill speeds.

Thermo-mechanical model of drilling
In this paper, a 3D FE model was developed using the explicit dynamic finite element code, LS-DYNA (LSTC, Livermore, CA, USA).The model aims to simulate the drilling process and to calculate the thermo-mechanical stress distribution in the bone material.The FE model consists of a set made up of a twist drill bit (Ø4 mm, point angle of 118 and a helix angle of 30 ) and a bone block with appropriate boundary conditions.The drill bit geometry was designed through the CAD software SolidWorks Õ (Dassault Systems S.A., France) and imported into LS-DYNA program for numerical simulation.The bone block was modelled with a rectangular shape with the following dimensions: 14 Â 10 Â 4 mm.For mesh generation, dual mesh density technology was used in order to improve the calculation accuracy and to reduce the processing time.The element size is critical for the computational cost of the simulation, solved with an explicit integration scheme.The selection of the element size should be stated ensuring equilibrium between computational efficiency and precision of the model. 17Since the main part in this type of analysis is the bone region located in the immediate vicinity of the drilled hole, a circular disc with 6 mm of diameter was formed in this part and a mesh discretisation was applied with an element size of 0.5 mm to capture high stress gradients during the drilling process.In the remaining block was used a coarse mesh with an element size of 1 mm.To define those two densities, separate simulations were run to obtain better convergence and smooth contours, and varying mesh densities were tried to find the optimal values that ensure the time of calculation and the stability of results.Previous mesh sensibility studies showed that the use of dual mesh density is a common practice in this kind of simulation to achieve both accuracy and time efficiency of the calculation. 12,18The element type chosen for the numerical model was the 3D Solid 164 (eight nodes with three degrees of freedom at each node in X, Y, Z directions), only used for explicit dynamic analysis.The final numerical model of the block and the drill bit were meshed with 20477 elements.In the explicit analysis, it is important to define nodal groups called components to apply the loads.The whole FE model of the bone drilling along with the names of all its components is displayed in Figure 1(a).
Following the dynamics characteristics involved in these processes, the block was fixed on the bottom, while the drill bit movement was explicitly simulated via a dynamic FE approach to rotate and move only about its own longitudinal axis (Y axis) with a specific drill speed (o) and feed rate (V) vertically downwards into the bone block, as shown in Figure 1(b).To explore the effects of feed rate and drill speed on the generated thermal stresses, the numerical simulations were performed using three different drill speeds (600, 800 and 1200 r/min) with a constant feed rate of 50 mm/min, and three different feed rates (25, 50 and 75 mm/min) with a constant drill speed of 800 r/min.
Thermal stress analysis is driven by strains generated in the structure by a temperature load.Through the temperature data obtained with experimental tests, the time temperature history was applied to the bone block component as a prescribed temperature load, producing the heating source effect inside the block.The initial temperature of the bone block component was defined as 19 C.

Material properties
It is well known that the FE simulation results are very sensitive to the material properties, and thus it is important to choose the most suitable material models in accordance with the analysis type.Materials subject to drilling are highly affected by large and high strain rates, which finally lead to failure.To define the bone block submitted to high impact deformation and the thermal expansion due to the temperature load, different material models were implemented considering different components of the model.
The material model for the bone block component was *MAT_ELASTIC_PLASTIC_THERMAL, which is the *MAT type 4 in LS_DYNA material library.This model allows the definition of temperature-dependent material coefficients in a thermoelastic-plastic behaviour and takes into account the thermal coefficient of expansion.The coefficient of thermal expansion allows the calculation of thermal strain and stress using the prescribed temperature data.
In the selection of material model for the hole component, it should be taken into consideration just the material models which consider the strain-rate with dependency of the materials plastic curve.According to Skrlec and Klemenc, 19 the three most common applied material models that consider the strain-rate effects are: Cowper-Symonds, Johnson-Cook and Zerilli-Armstrong.1][22] In the FE model, the Cowper-Symonds law is implemented through the *MAT_PLASTIC_KINEMATIC (*MAT type 3) in LS-DYNA material library.Yield stress according to the Cowper-Symonds material model is defined with equation (1) 23,24 where s0 is the reference yield stress in MPa, _ " the strain rate in s À1 , b is the hardening parameter (between 0 for kinematic hardening and 1 for isotropic hardening), C and P are the Cowper-Symonds strain rate parameters, " eff p the effective plastic strain and Ep the plastic hardening modulus which

Fernandes et al.
is dependent on the Young's modulus E, and the E tan tangent modulus is given by The mechanical properties of the bone block component and hole component were obtained from the uniaxial tensile tests and have been comprehensively defined in our previous studies. 25][28] The material parameters that consider the strain rate effects are not simply measured and determined, some of them are empirically determined through special experimental and optimisation process. 19he drill bit was assumed to be a rigid body, since its stiffness is much higher than the bone.This is a practice way to reduce computational cost and resources involved in the drilling simulation.All the material properties used in the numerical analysis are listed in Table 1.
In equation ( 2), the strain rate influences only the yield stress, which means that the plastic curves (the flow stress as a function of strain) are parallel. 19he larger the strain rate, the higher the flow stress, as represented in Figure 2, according to the experimental data published in the literature for other materials. 19ntact interaction and failure criteria The contact algorithm was *CONTACT_ ERODING_SURFACE_TO_SURFACE, which is accessible in the LS-DYNA program. 24This type of contact is used when a body surface penetrates the surface of another body.Eroding contact methods are suggested when solid elements involved in the contact definition are subject to erosion (element deletion) due to material failure criteria.The surface of the drill bit was classified as the contact surface.The target surfaces include all the block surfaces that will be contacted with drill bit.The friction behaviour between the drill bit and the block is assumed to be governed by Coulomb's friction law.In this case, for the frictional contact between the drill bit and the block, a constant coefficient of friction of 0.3 was used. 29,30 failure criteria was applied to simulate the element removal and hole generation.As mentioned above, the approach of element erosion can be used to simulate the perforation when any element exceeds a specified plastic strain.These elements are deleted from the solution after they fail.Based on the bone properties, the failure strain reaching 0.05 is adopted as the criterion in the erosion algorithm implementation for the numerical simulation.Dynamic analysis  was used with the simulation range subdivided into 15,000 time increments of 8.0 Â 10 À4 s.ANSYS/ LS-DYNA requires very small time steps with many iterations to ensure stability of solution.Each simulation took about 68 h to execute on a workstation with a quad-core Intel I7-4790k with 16 GB RAM.

Experimental validation
An experimental methodology was executed to validate the numerical thermo-mechanical model of drilling.Solid rigid polyurethane foams were used as an alternative material to cortical bones, due its consistent and homogeneous structural properties. 31,32The samples were supplied in rectangular shape with the dimension of 130 Â 180 Â 40 mm and the material has a closed cell with density of 0.80 g/cm 3 (Figure 3(a)).
In order to validate the numerical model, the drilling tests were divided in two groups.For the first group, the holes were performed at drill speeds of 600, 800 and 1200 r/min with a constant feed rate equal to 50 mm/min.In the second group, three different feed rates (25, 50 and 75 mm/min) were applied throughout the drilling procedures, with a constant drill speed of 800 r/min.Drillings were performed in the Mechanical Laboratory at the Polytechnic Institute of Braganc¸a.In total, 30 holes with 30 mm of depth were made at room temperature (without cooling) using a twist drill bit of 4 mm diameter, point angle equal to 118 and helix angle of 30 .A control of the drilling parameters was provided by a computer numerically controlled (CNC) system.All drills were used no more than 15 times before being replaced with a new one.For each combination of parameters, the average of six drillings was used to present the results.

Strain measurement
Fifteen flexible linear strain gauges (1-LY18-6/120, 120 AE 0.35% from HBM) were attached on the blocks' surface to measure the strain during the drilling process (Figure 3(b)).The blocks were properly clean to create a uniform surface and the locations of the holes were marked to always keep the same distance of 3.5 mm, between the edge of the hole and the strain gauge.All strain gauges were identified for each channel and attached with bonding adhesive following the instructions of the manufacturer.Each lead wire was connected to the quarter bridge in a data acquisition system (Vishay Micro Measurements P3 Strain Indicator and Recorder).Strain measurements were taken continuously during each step of the drilling procedure.The corresponding profiles of stresses versus drilling depth from the block surface were calculated.

Temperature control
Temperature measurement was carried out using two methods.In the first method, thermography (ThermaCAM 365, FLIR Systems) was used with the lens located at distance of 1.5 m from the drilling area.This method allowed to obtain thermal images of the blocks and the drill bit surfaces, before and immediately after the drilling.The measured temperature is function of the surface conditions, represented by their emissivity.The imposed parameters to the camera during image acquisition are listed in Table 2.
In the second method, a set of chromel-alumel Ktype thermocouples (Omega Engineering Inc., Stamford, USA) installed close to the hole surface  were used for measuring the temperature inside of the block.The thermocouples with wire diameter of 3.5 mm were installed in two opposite sides of the block and at the same distance from the hole for all tests, as shown in Figure 4.After the thermocouples have been placed, the entrance of each hole was appropriately sealed with glue to fix the thermocouples and isolate them from external heat.All thermocouples were connected to a data acquisition system (HBM MGCplus) for the acquisition of temperature along the drilling time.
The temperature rise inside of the material during drilling was used as prescribed temperature in the FE model.An average of the registered temperature values in different holes through the drill time was considered.These curves were imposed as prescribed temperature on the bone block component.Figure 5 shows the evolution of the mean temperatures recorded inside of the material considering the different drill parameters.As the numerical model was modelled with only 4 mm in height, the temperature data were collected only for the correspondent drilling times.For each feed rate were considered different drilling times, ensuring the full block drilling.
Each block was prepared to accommodate eight holes with a distance of 20 mm between them.The complete experimental setup is shown in Figure 6.In this setup, the position of the drill is fixed and a CNC machine was used to perform the holes.All experiments started from room temperature (approximately 19 C) and the drilling direction was perpendicular to the block axis.Between successive experiments, sufficient time was allowed for the block and the drill bit to return to initial conditions.
Previous publications from numerical and experimental tests showed that temperature measurements, on the block top surface at 3.5 mm from the drilled zone, represent small values with no interference in the strain measurements. 7,25sults and discussion

FE results
The ability to measure the thermo-mechanical stress is crucial to predict the behaviour of bone tissue during drilling and also to aid in developing predictive capabilities by verifying models.The overall  objective of this study was to analyse the thermomechanical stress as a function of the different drill parameters.
The results obtained by using the FE model are given under this section.A drilling depth of 3 mm was chosen to show the von-Mises stress distribution at different feed rates and drill speeds in the numerical simulations.In order to evaluate each drill parameter individually, the drill speed was fixed at the value of 800 r/min to measure the effect of feed rate and fixed at the value of 50 mm/min to measure the effect of drill speed.Figure 7 presents the contours of von-Mises stress distribution in the bone drilling for each drill situation.
For the same time instant was observed through the numerical results that the von-Mises stresses  increase with increase of feed rate and drill speed.However, this trend becomes more evident with regard to the feed rate parameter.This may be explained in part by the differences between drilling depths at different feed rates.For the same time instant, the holes made with a feed rate equal to 75 mm/min reach a higher depth than feed rates equal to 25 and 50 mm/min.It is also important to emphasise that the maximum values of von-Mises stresses in all numerical simulations were verified in the drilled zone and its immediate vicinity, exceeding the 5 MPa.During the numerical simulation, the elements in the drilled hole were removed when the plastic strain failure in the material was reached.

Experimental results and FE model validation
Experimental validation of numerical model was performed by comparing the results obtained in both models.The average of stress component (s zz ), which represents the normal stress on face perpendicular to the z axis, was calculated experimentally and compared with the numerical results at different drilling times.The distance between the edge of the drilled hole and the linear strain gauge, as well as the direction of the gauge was considered on the selection of the nodal points in the FE model.Figure 8(a) and (b) shows the mean and standard deviation of normal stresses in one direction (s zz ) obtained in both methods at different drilling times.In addition to the numerical simulations, the drill speed was fixed at 800 r/min to measure the effect of feed rate and fixed at 50 mm/min to measure the effect of drill speed.
Experimental and numerical results indicate a strong relationship between either feed rate, drill speed or drilling depth and the stress levels.Overall, it was observed that the normal stresses in drilling process increase with the increasing of feed rate and drill speed, for the same time instant.Regarding the drilling depth, experimental and numerical results showed a clear trend for an increase in level of stresses with tool penetration along the time, reaching a maximum value when the drill bit penetrated completely the block.The greater the drilled hole depth, the greater is the normal stresses in the bone block.
After analysing the bone behaviour at different drilling times, it is also important to analyse at the end of the process.The effect of feed rate and drill  speed for the complete penetration of the block was examined.Figure 9(a) and (b) presents the average of maximum normal stresses (s zz ) for both models analysed.
Analysing the entire process, it was observed that the maximum normal stresses retained the same behaviour with respect to the drill speed.Regarding the feed rate, the normal stresses decreased with the increase of feed rate.Although at the start of drilling, the generated normal stresses are higher for higher feed rates (75 mm/min), in the entire process are found higher stresses for lowest feed rate (25 mm/min) because the drilling time also increases.In generally, the numerical results were validated and are in good agreement with the experimental results.

Conclusions
In this paper, a three-dimensional thermo-elasticplastic dynamic model based on finite element method was built to simulate the drill bit penetration into the bone tissue.Thermo-mechanical stresses involved in the drilling processes were obtained and compared with the experimental tests.Based on numerical and experimental results, the following conclusions can be drawn: . The thermo-mechanical stresses generated in the material increase with tool penetration and, consequently, with increase of hole depth. .The normal stresses (s zz ) tend to decrease when the feed rate is higher and increase when the drill speed is higher. .The maximum values of von-Mises stresses are found in the drilled zone and its immediate vicinity for all situations. .The results presented here demonstrate that the experimental methodology coupled to the developed numerical model is an effective procedure for evaluating the thermo-mechanical behaviour of the contact region between the drill bit and the bone.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Figure 1 .
Figure 1.(a) FE model of drilling and (b) boundary conditions.

Figure 2 .
Figure 2. Flow stress for the Cowper-Symonds model.

Figure 3 .
Figure 3. (a) Solid rigid polyurethane foam and (b) installation of the linear strain gauges.

Figure 5 .
Figure 5. Temperature data applied in the Fe model.

Figure 7 .
Figure 7. Distribution of von-Mises stress (MPa) at different feed rates and drill speeds (3 mm depth).

Figure 9 .
Figure 9. Maximum normal stresses s zz (MPa) with (a) feed rate and (b) drill speed.

Figure 8 .
Figure 8. Normal stresses s zz (MPa) in the experimental and FE models with (a) feed rate and (b) drill speed.

Table 1 .
Material properties and Cowper-Symonds parameters used in numerical simulations.

Table 2 .
Parameters used for thermal images acquisition.