- Email: [email protected]
Optimization of bone drilling process based on finite element analysis
Applied Thermal Engineering 108 (2016) 211–220
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Optimization of bone drilling process based on finite element analysis Xiashuang Li a, Wei Zhu a, Junqiang Wang b,⇑, Yuan Deng a,⇑ a b
Beijing Key Laboratory for Advanced Functional Materials and Thin Film Technology, School of Materials Science and Engineering, Beihang University, Beijing 100191, China Department of Orthopaedics and Traumatology, Jishuitan Hospital, Beijing 100035, China
h i g h l i g h t s
g r a p h i c a l a b s t r a c t
Volumetric heat source is loaded to
obtain more accurate simulation results. Multiple parameters’ synergistic effects on the drilling temperature are studied. Empirical formula with sensitivity analysis is proposed to predict the temperature. Intermittent feed drilling is demonstrated to reduce the drilling temperature.
a r t i c l e
i n f o
Article history: Received 1 October 2015 Revised 25 May 2016 Accepted 18 July 2016 Available online 19 July 2016 Keywords: Bone drilling Feed rate Spindle speed Drill diameter Intermittent feed drilling Finite element method
a b s t r a c t Bone drilling is frequently used in modern surgery operations, especially in fracture treatment. Drilling temperature should be strictly controlled, since when the temperature is 47 °C higher over 60 s, irreversible damages happen to bone. The study investigated the heat transfer of bone drilling by a threedimensional model based on finite element method. Three principal drilling parameters were studied: feed rate, spindle speed and drill diameter. A parametric study proved that drilling temperature increases once any of three parameters rises. Parameters’ effects on drilling temperature are proved synergistic. An empirical formula of maximum drilling temperature depending on three drilling parameters is given with a sensitivity analysis. To a drilling process of surgical operations, this formula could be used for temperature prediction and optimization of parameters respecting drilling temperature. Furthermore, the intermittent feed drilling is proposed and demonstrated to be very efficient to reduce thermal necrosis to bone tissue. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction Bone drilling is a significant part of the internal and external fixation processes in orthopedic surgery. The heat generated during drilling may be partially dissipated by the blood and tissue fluids and partially carried away by the chips formed. However, due to the poor thermal conductivity of fresh cortical bone, which is
⇑ Corresponding authors. E-mail addresses: [email protected] (J. Wang), [email protected] (Y. Deng). http://dx.doi.org/10.1016/j.applthermaleng.2016.07.125 1359-4311/Ó 2016 Elsevier Ltd. All rights reserved.
among 0.38–2.3 Wm1 K1, it’s difficult for bone to conduct the rest heat away from the cutting edge. So, the temperature rise at the cutting edge in a deep cortical hole would be extremely high. The phenomenon of bone injury, resulting from the rise of heat during drilling operations, is defined as thermal necrosis. This damage to bone cells would delay the healing process after the surgery and reduce the strength of the fixation. Eriksson and Albrektsson [1] observed that the cortical necrosis and delayed healing occur in living animal bone heated at 47 °C for one minute. Krause [2] states that osteoclasts begin to die when the drilling temperature is higher than 50 °C. According to the study of Augustin [3],
212
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Nomenclature {b} a c c1 c2 c3 c4 d [D] Dep Din Dout e f h Kxx Kyy Kzz l
unit outward normal vector intensity coefficient specific heat [J kg1 K1] first regression coefficient second regression coefficient third regression coefficient fourth regression coefficient diameter of drill [mm] thermal conductivity matrix drilling depth [mm] inside diameter of bone model [mm] outside diameter of bone model [mm] constant feed rate of drilling [mm min1] heat convection coefficient [W m2 K1] elemental thermal conductivity in the direction of x [W m1 K1] elemental thermal conductivity in the direction of y [W m1 K1] elemental thermal conductivity in the direction of z [W m1 K1] length of bone model [mm]
thermal damage was caused to the bone while temperature was between 47 °C and 50 °C. The thermal damage on bone tissue during orthopedic sawing was also investigated and the threshold temperature found was about 55 °C [4]. In order to minimize the damage caused by the high temperature (above 47 °C) in bone drilling operations, it is necessary to optimize the drilling parameters. Many researches have been conducted to find out effects of different drilling parameters on multiple quality characteristics of bone drilling. Drilling parameters such as spindle speed, feed rate, drilling depth, drill tip angle, drilling force and drill diameter were experimentally studied [5–7]. Pandey and Panda [8] used the analysis of variance (ANOVA) to determine the effect of each drilling parameter on drilling temperature. Beside, drilling experiments on bovine bone with spindle speeds from 800 rpm to 3800 rpm were performed by Lee et al. [9], and it is concluded that the drilling temperature rises with increase in spindle speed. Other methods like drilling experiments on bone-like material Poly Methyl Methacrylate (PMMA) were also used for prediction [10]. Previous researches were mainly aimed to investigate the effect of individual drilling parameter on the bone drilling temperature. This paper studies multiple parameters’ synergistic effects on drilling temperature with a three-dimensional model. The objective is to build an instruction for classic bone drilling operation in form of a general empirical formula that is concise and simple for using. Recently the numerical simulation based on finite element method is a possible substitute for high-cost and complex experimental work, especially as a useful tool for verification of analytical results and prediction. ANSYS is powerful and practical FEM-based numerical simulation software for analyzing complicated structure. ANSYS Mechanical is capable to analyze the thermal response of structures to heat transfer effects, involving conduction, convection, and radiation heat transfer [11]. It is widely used for modeling heat transfer problem owing to its excellent performance in solving problem of different modes including steady, transient, linear and nonlinear. Thiagarajan, King et al. [12] investigated a steadystate heat transfer simulation in ANSYS to quantify thermal losses of the pool boiling. Zhou et al. [13] simulated oxygen cutting, similar with procedure of drilling, using a composites heat source model to predict the cutting temperature in ANSYS.
L vector operator n spindle speed of drill [rpm] q_ drilling (x, y, z, t) volumetric heat generation rate [W m3] t time [s] T(x, y, z, t) current temperature of the bone [K] T1 ambition temperature [°C] empirical value of TMAXm issued from the regression Temp [°C] Texp experimental value of TMAXm issued from the model [°C] Thuman human body’s temperature [°C] TS temperature on the surface S [°C] TMAX(t) maximum temperature on the bone at moment t (drilling temperature) [°C] TMAXm maximum value of TMAX during a drilling process [°C] V velocity vector for mass transport of heat [m s1] vc cutting speed of drill [m s1] vf feed rate of drilling [m s1] q density [kg m3]
The Finite Element Method (FEM) presents effective performance not only in virtually all areas of engineering but also in applied science like medical science. The FEM could avoid the high-cost in term of expense and time of experimental study in medical science. For example, ANSYS is used for the simulation of biomaterials like bone and dent to investigate their mechanical and thermal properties. This paper analyzes the temperature’s change during cortical bone drilling with different parameters, including drill spindle speed n, feed rate f and drill diameter d via ANSYS. Specifically, with customizing ANSYS by APDL (ANSYS Parametric Design Language) script, a three-dimensional heat transfer model of bone was built and verified. In addition, a volumetric heating source was used in this model. Li et al. [14] compared volumetric and surface heating sources in the modeling of laser melting of ceramic materials and concluded that the simulation result with volumetric heating source is more accurate and stable. Other applications of finite element method were also used to investigate the mechanism of cutting process and to predict drilling temperature. Tu et al. [15] used a three-dimensional finite element model to show that lowering the initial temperature of Kirschner pin can decrease the temperature rise as well as the size of the thermally damaged zone. Mokhtar and Fawad [16] worked out a complete theory for bone drilling modeling to obtain the heat flux. A new thermal model has been worked out by Lee et al. [17] for bone drilling with a sensitivity analysis and a single parametric study was also investigated for this model. Sezek et al. [18] worked with a FEM-based model for bone drilling and found a safety zone for drilling parameters. To date, most studies of bone drilling modeling proposed a new model or a new theory without further practical research of the drilling temperature by exploring the model or the theory. This paper presents a detailed investigation from the foundation to the validation of the three-dimensional model, followed by results of different levels. Hundreds of simulations were completed in ANSYS with values that correspond with the real drilling procedure. The influence on drilling temperature is studied separately with single parameter and double parameters. The effects of multiple parameters on drilling temperature were proved synergistic. In order to generalize the numerical results, an empirical formula
213
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
is given to predict the drilling temperature in function of these three drilling parameters. At last but not the least, an innovative intermittent feed drilling method is also proposed and proved by this model to effectively reduce the drilling temperature. 2. Materials and methods 2.1. Thermal model ANSYS uses the finite-element method to solve the underlying governing equations and the associated problem-specific boundary conditions. With the presence of internal heat source due to the drilling process, the balance relation of the heat flow of a volume bounded by an arbitrary surface S is as follows [19]:
qc
@Tðx; y; z; tÞ þ fVgT fLgT þ fLgT fqg ¼ q_ drilling ðx; y; z; tÞ @t
ð1Þ
where q is the density, c is the specific heat, and T is the locationtime dependent temperature. {L} is the vector operator, {L} 8@9 < @x = @ = @y ; {V} is the velocity vector for mass transport of heat, {V} :@; 8 @z 9 < vx = = v y . q_ drilling (x, y, z, t) is the volumetric heat generation rate, : ;
vz
caused by the drilling movement and it’s given as follows [20,21]:
v 0:85 q_ drilling ðx; y; z; tÞ ¼ av 0:75 f c where a = 1E10, the coefficient of the intensity;
ð2Þ
vf
f ¼ 60000 , the feed
rate and v c ¼ p the cutting speed. {q} is the heat flux vector, and Fourier’s law is used to relate the heat flux vector to the thermal gradient: n 2 60
d , 2
fqg ¼ ½DfLgT
ð3Þ
where the conductivity matrix [D] is given by: [D] 2 3 K xx 0 0 = 4 0 K yy 0 5. Kxx, Kyy, Kzz are the elemental conductivity 0 0 K zz respectively in the direction of x, y, and z. {V} is the velocity vector for mass transport of heat, based on the proposed model we ignore the effect of the drilling chips carrying away heat on the temperature field of the drilling site, i.e. letting {V} = 0 in (1) and substitute (2) into it, the thermal governing equation of this quasi-steady heat conduction problem becomes:
qc
@Tðx; y; z; tÞ ¼ q_ drilling ðx; y; z; tÞ þ fLgT ½DfLgT @t
ð4Þ
The boundary condition is considered. Specified convection heat loss acting over non-drilled surface S (Newton’s law of cooling) is given as:
fqgT fbg ¼ hðT S T1 Þ
ð5Þ
where {b} is the unit outward normal vector, h is the heat convection coefficient acting on the surfaces of non-drilled bone (h = 3 Wm2 K1), T1 is the ambition temperature considered as 20 °C, TS is the temperature on the surface S, the radiation heat loss is negligible compared to the convection heat loss. The initiative condition is:
Tðx; y; z; t ¼ 0Þ ¼ T human
ð6Þ
Thuman is the average temperature of a human body, meaning 37 °C. In many simulations [15,17] and even experimental investigations [23] of the bone drilling process, 37 °C was set as the initial bone temperature.
2.2. Numerical modeling and loading In ANSYS, the process of the bone drilling was simulated by a model including a hollow cylinder and a mobile thermal source, respectively representing the bone and the cutting edge of drill. ANSYS is indeed proficient in modeling complex structure by introducing the physical model built in Solidworks. However, the drilling happens on very small part of bone so that what the problem is in fact local. Hence, the complex form of bone has very little change the simulation’s result that is ignorable. In many previous investigation of bone drilling problem, a simplified geometry like a hollow cylinder is often used for the bone and they all get satisfied results [12,22]. The APDL mode is used due to its flexibility over the Graphics User Interface (GUI) mode for the simulation.SOLID90 is chosen as the element type that has 20 nodes with a single degree of freedom, temperature, at each node. The 20-node elements have compatible temperature shapes and are well suited to model curved boundaries. It’s possible to set up several materials properties for the element such as, K xx ; K yy ; K zz ; q and c. Besides, this element could carry surface loads like convections and body loads like heat generations. Material properties of the bone used in the model are listed in Table 1. The physical model is a hollow cylinder as shown in Fig. 1(a) and (b) and obeying to a homogenous and anisotropic material whose geometry parameters are given in Table 2. As for the mesh generation, dual mesh density technology was used in order to improve the calculation accuracy and to reduce the processing time. Since the temperature gradient would rise in the drilling region, a higher mesh density was needed for the area near of drilling site, noted S2. While for the rest of model, a low mesh density was enough, noted S1. To define those two densities, an approximation method was applied to find the optimal values that ensure the time of calculate and the stability of results. Four pairs density were tested as follows: A1: A2: A3: A4:
S1 = 0.01 mm, S1 = 0.02 mm, S1 = 0.01 mm, S1 = 0.01 mm,
S2 = 0.001 mm; S2 = 0.001 mm; S2 = 0.0005 mm; S2 = 0.0008 mm.
The maximum drilling temperature calculated for each four groups are given in Table 3. The difference in temperature of these four groups is small enough to be ignored so the difference in calculation time is decisive to choose the density pair. Group A2 has the shortest calculation time but singularity problem appears, as the mesh is not smooth enough due to the large difference of densities in this pair. Finally, the length of meshed elements is chosen as 0.001 mm for high mesh density and 0.01 mm for low mesh density. The mesh of the transition part on the boundary of two regions were refined and smoothed as necessary to create the appearance of a natural transition avoiding singular points. The final meshed model is presented in Fig. 1(c). The convection and the heating source are two loads to be applied. The air convection between the bone and the environment was loaded with the heat convection coefficient of 3 W m2 K1 on all the surfaces of bone. The heat generation due to the drilling was modeled by a volumetric heating source that moves from the exterior to the interior of the bone following the drilling path until the end of drilling depth. To apply this mobile load, the to-be-drilled part is divided in six sections (six cylinder), numbered from 1 to 6, as shown in the Fig. 1(d). So that the volumetric heating source could move from section number 1 to section number 6 successively. At each time increment, this stepping load of heat generation rate is applied on only
214
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
one piece of these six sections to simulate the local heating resulting from the drilling. Then in the next time increment, the volumetric heat source was only applied on the next section. 2.3. Numerical simulation Three principal drilling parameters have been chosen for numerical calculations, same as the situation of the real drilling processes. They are f the feed rate, n the spindle speed, and d the diameter of the drill. The bone temperature is the target during the drilling process, so the maximum temperature of the entire bone is extracted. Numerical simulations were investigated for single parameter study at first and the typical values have been used for each parameter, as shown in the Table 4. Then numerical simulations were done for double parameter study and the test values are listed in Table 5. In the case of the single parameter study, the more values are tested, the more accurate is the result. However, for a double parameter study, it’s not the quantity but the representativeness and equal interval of test values that determine the accuracy of result. Therefore, values for the double parameter study were well selected, as shown in Table 5.
Table 1 Material properties of the bone. Item
Value
Density (kg m3) Specific heat capacity (J kg1 K1) Thermal conductivity (W m1 K1) at 20 °C Initial temperature (°C)
1700 1260 0.38 37
Table 2 Geometry parameters of the physical model. Item
Value
Outside diameter (mm) Inside diameter (mm) Length (mm) Drilling depth (mm)
100 70 100 10
Table 3 Calculation results of different density pairs. Group
TMAXm (°C)
A1 A2 A3 A4
44.012 44.011 44.002 43.985
Table 4 Tested values for single parameter study. Drilling parameters f (mm min n (rpm) d (mm)
1
)
Test values 30–50–60–75–100–120–150 150–300–600–900–1200–1500–1800 2.0–3.2–3.5–4.0–4.5–5.0
Table 5 Test values for double parameter study. Drilling parameters
Test values
f (mm min1) n (rpm) d (mm)
25–50–75–100 300–600–900–1200–1500–1800 2.0–2.5–3.0–3.5–4.0–4.5–5.0
Fig. 1. (a) Cross section of the hollow cylinder; (b) Physical model; (c) Meshed model; (d) Zoom of the part of bone to be drilled.
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Fig. 2. Time dependent drilling temperature with different feeding rates (n = 900 rpm, d = 3.2 mm).
3. Results and discussions 3.1. Single parameter results Firstly, the feed rate f has been chosen as the only variable for numerical calculations. The values of f used are in Table 4. The spindle speed n and the drill diameter d was kept constant: n = 900 rpm, d = 3.2 mm. Fig.2 displays the change of the timedependent maximum drilling temperature corresponding to different values of f. Initially the temperature was set at 37 °C. At each moment t, the maximum temperature of the entire bone is noted as TMAX, which is the drilling temperature, represented by one point in the Fig. 2. The time of drilling is shorter when the feed rate f increases for the same drilling depth. In general, the TMAX increases rapidly at first and becomes almost stable until the end of drilling process. In fact, the maximum temperature of the bone appears at the cutting edge and rises rapidly in the first several seconds. Gradually the drilled part cools down and the maximum temperature of bone always appears at the point where the drilling is going on. The value of the maximum temperature of bone in each drilling repeated keeps almost stable. For each trial, the TMAX reaches its maximum value, noted as TMAXm, in the first two or three seconds. In Fig. 2, the drilling temperature TMAX increases with f increasing. Keeping n = 900 rpm, d = 3.2 mm, TMAXm goes over
215
47 °C when f is higher than 75 mm min1. Fig. 2 indicates that f has a considerable influence on drilling temperature. When f goes from 30 mm min1 to 150 mm min1, the TMAXm increases from 40.63 °C to 51.55 °C, which means a relative increase of 26.85%. In a real case, when the drill is fed more rapidly, the drill force is stronger and so is the friction force. This phenomenon will generate more heat on the friction area in the same period of time. Same kind of calculation was conducted for different d and different n whose results are given in Fig. 3. As same as the situation for f, the TMAX of each trial reached its maximum value in the first two or three seconds of drilling process. The larger the parameter’s (n or d) value is, the longer it takes to reach the TMAXm. TMAXm rises with n and d in different speed. The relative increase is about respectively 14.45% for d going from 2.0 mm to 5.0 mm and 28.94% for n going from 150 rpm to 1800 rpm. In practice, when the drill’s diameter becomes larger, so does the friction area between the drill and the bone tissue, so more heat is generated. The influence of the increment of n or f on the increase of TMAXm increases is bigger than that of d. In reality, more heat is generated during the drilling process due to the increase of each of three parameters. However, the larger the drill diameter (d) is, the larger the diameter of the hole is, so that there is more surface for heat to transmit. In this way, the drilling temperature will not rise as much as that in the other two cases. 3D diagrams of TMAX in function of one single drilling parameter and time show the growth of drilling temperature more directly. Fig. 4(a) shows the change of TMAX with n and time. It manifests that for n higher than 1700 rpm, 1.5 s after the beginning of drilling, TMAX exceeds the security temperature (47 °C). It tells the maximum continuous time the drilling before causing thermal necrosis for a value of n given. It could also be used to find out possible values of n when the drilling time is set. Fig. 4(b) is the projection of this 3D diagram, from which statistics can be gathered more easily. Same kind of 3D diagram was also completed for f and time or d and time. 3.2. Double parameters results Furthermore, in order to evaluate three parameters’ different level of influence and their synergistic effects on drilling temperature, two groups of double parameter were considered as covariables: d-f and n-f. Tested values for each parameter are listed in Table 4. Keeping n constant (n = 900 rpm), for each possible value of d, numerical experiments were conducted with different f. In this part, 4 7 simulations were investigated and the results are presented in Fig. 5.
Fig. 3. Change of drilling temperature with (a) different d, f = 75 mm min1, n = 900 rpm; (b) different n, f = 75 mm min1, d = 3.2 mm.
216
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Fig. 4. (a) 3D diagram of TMAX in function of n and time; (b) Projection of 3D diagram of TMAX in function of n and time.
Fig. 5. (a) Different effects of d and f on drilling temperature; (b) Different effects of f and d on drilling temperature; (c) 3D diagram of TMAXm in function of d and f; (d) Projection of 3D diagram.
In Fig.5(a), each curve represents the change of TMAXm in function of f with different given values d. As mentioned in 3.1, for each bone drilling simulation, TMAXm represents the maximum value of TMAX during the whole drilling process. The larger the d is, the
more TMAXm grows with f. For each f given, the change of TMAXm in function of d is shown vertically: The higher the feed rate f is, the faster TMAXm increases with d. This means that the effects of d and f on the drilling temperature are synergistic even though
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
217
Fig. 6. (a) n and f on drilling temperature; (b) Different effects of f and n on drilling temperature; (c) 3D diagram of TMAXm in function of n and f; (d) Projection of 3D diagram.
they are independent from one to another. In practice, when the drill feeding is faster, more heat will be generated in the drilling process with a large drill diameter than that of a small drill diameter so that the rise of drilling temperature is bigger. When the drill diameter is larger, a larger feed rate will generate much more heat on the increased friction area than a small feed rate and consequently the drilling temperature increases more. In Fig.5(b), each curve indicates the change of TMAXm along with d for f disposed. It’s a different way to interpret the same data as in Fig. 5(a). These two diagrams prove that d = 3.5 mm is a safe drill diameter for n = 900 rpm and f 6 100 mm min1 and f = 50 mm min1 is a safe feed rate for n = 900 rpm and d 6 5 mm. Fig. 5(c) is the 3D diagram of TMAXm in function of d and f. This kind of diagram helps visualize the security zone (below 47 °C) bounded by different values of f and n, for a drill diameter d given. In this situation, non-security zone (over 47 °C) is colored in red1 and orange. In Fig. 5(d), the 3D diagram is projected to the plan of d and f, which makes it less intuitionistic but more convenient to read values. Same type of study was implemented keeping d constant (d = 3.2 mm) and variable n and f. 4 6 simulations were investigated in total for this part and the results are shown in Fig. 6. In a similar way, Fig. 6(a) and (b) confirm the synergistic effects of n and f on the drilling temperature: TMAXm increases more rapidly 1
For interpretation of color in Figs. 5 and 7, the reader is referred to the web version of this article.
with n when the feed rate f is larger, and vice versa. From these two diagrams it can be told that n = 1200 rpm is a safe drill diameter for d = 3.2 mm and f 6 100 mm min1 and f = 50 mm min1 is a safe feed rate for d = 3.2 mm and n 6 1800 rpm. Fig. 6(c) is the 3D diagram of TMAXm in function of n and f. The boundary between the security zone and the non-security zone is the one between the green band and the yellow band. Fig. 6(d) shows a projection of the 3D diagram. 3.3. Empirical formula In order to quantify the effects of parameters on the drilling temperature, an empirical formula of drilling temperature was built in MATLAB. The model was linearized at first considering these three drilling parameters as predictor variables. Then a multiple linear regression was realized with the function regress. Supposing the formula of TMAX in the following type: c2 c3
TMAX m ¼ ec1 f d nc4
ð7Þ
To linearize the formula, logarithm natural is taken at both sides, which makes:
ln ðTMAX m Þ ¼ c1 þ c2 ln ðf Þ þ c3 ln ðdÞ þ c4 ln ðnÞ
ð8Þ
52 groups of numerical calculations results were entered in MATLAB to determine the constants in the formula. Finally, the empirical formula obtained is:
218
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Fig. 7. (a) Residual case order plots; (b) Comparison between TMAXm experimental (Texp) and TMAXm empirical (Temp).
TMAX m ¼ e2:4812 f
0:1092 0:1387 0:0995
d
n
ð9Þ
A local sensitivity analysis is done to figure out how the uncertainty of the TMAXm can be distributed to different sources of uncertainty of the three parameters. The sensitivity of this empirical formula with changes of each parameter is equal to the partial derivative of TMAXm respectively to each. Obviously, the results correspond to the exponent number of each parameter, meaning 0.1092 for f, 0.1387 for d and 0.0995 for n. The empirical formula of maximum drilling temperature shows that the TMAXm increases non-linearly with the change of those parameters. This analysis also illustrates that in these current units of parameters, TMAXm is more sensitive to the drill diameter than to the feed rate and least sensitive to the drill spindle speed. During bone drilling, the operator need to be more conscious with a large drill diameter since the drilling temperature will arise faster in this situation with the augmentation of other drilling parameters. For a disposed bone drilling, the feed rate should be decreased firstly if the predicted drilling temperature is critical. Fig.7(a) displays an error bar plot of the confidence intervals on the residuals from this regression. The interval around the 29th, 51th and 52th residual (shown in red) does not contain zero. This indicates that the residual is larger than expected in 95% of new observations, and suggests the data point is an outlier. A comparison between TMAXm experimental (Texp) and TMAXm empirical (Temp) demonstrates that the empirical formula is reliable enough, as shown in Fig.7(b). As a concise and compact representation of the numerical results, this reliable empirical formula quantifies the influence of drilling parameters on maximum drilling temperature. Before an orthopedic surgery with a bone drilling procedure, by entering the values of these three drilling parameters, a prediction of maximum drilling temperature could be made. If this temperature exceeds 47 °C, the safe threshold, the empirical formula could tell which parameter has the most significant influence on drilling temperature. And the value of this parameter is mathematically be diminished in advance.
drilling actions. Simultaneously, the convection with air in the environment helps the bone to cool down. To complete the same work of drilling, same quantity of heat is generated but attributing to the intermittent feed, heat is better conducted or scattered in this case. The drilling temperature is effectively reduced as the result of simulation showed with the FEM based model in ANSYS. In this paper, the number of passages means how many steps the drilling is divided into and rest time is the time between two successive steps of drilling. Four groups of numerical experiments were carried out with same drilling parameters f = 75 mm min1, d = 3.2 mm and n = 1500 rpm. In the first group, the group of reference, the drilling process was continuous. In the three rest groups, the drilling process was respectively divided into 2, 6 and 12 passages with different rest time. The change of TMAXm in time of each simulation is given in Fig.8. In all the diagrams, the curves with squares represent the same one continuous process. The curve with disks and the curve with triangles in the same diagram are processes of intermittent feed with same number of passages but different rest time. For example, in Fig. 8(a), the curve with disks represents a process in 2 passages with a rest time of 2 s and the curve with triangles represents a process in 2 passages with a rest time of 5 s. When the drilling process is divided into 2 or 6 passages, the drilling temperatures did not decreased, regardless of the length of rest time. When the number of passages reaches 12 as presented in Fig. 8(c), the drilling temperature is obviously reduced: from 48.5 °C that is over the security temperature to below the security temperature 47 °C. By all means, even when the passages number is not large enough to decrease the drilling temperature as shown in Fig. 8 (b), the effective duration in the case of the drilling temperature is over 47 °C is reduced. The bone remained above 47 °C for 12 s with a continuous feed while for just 1 or 2 s, 6 times with an intermittent feed in 6 passages. Hence, the intermittent feed drilling is very effective to lower the irreversible damage on human body resulted by the high temperature during surgical operations.
4. Conclusions 3.4. Intermittent feed drilling An innovative method of feed, intermittent feed is discussed in this part to decrease the drilling temperature. Intermittent feed reduces the continuous drilling time as well as the continuous heating time on bone. This leaves the bone more time to conduct the generated heat away from the cutting edge between two
In this paper, a three-dimensional heat transfer model based on finite element method was built with a volumetric heating source to investigate the thermal process of bone drilling. The study of this model showed that maximum drilling temperature during the drilling process increases with each of these three parameters including feed rate f, drill spindle speed n and drill diameter d.
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
219
Fig. 8. Change of TMAXm of drilling in (a) 2 passages; (b) 6 passages; (c) 12 passages with 2 s or 5 s of rest between passages.
Moreover, their effects on the drilling temperature are synergistic. For the same increase of f, the increase of maximum drilling temperature is higher when n or d is higher. For the same increase of n or d, the increase of maximum drilling temperature is higher when f is higher. Based on numerical results from this model, an empirical formula of maximum drilling temperature that depends on three different drilling parameters was provided to quantify the effects of parameters and to predict the maximum temperature for these three parameters given. The simulations with this model proved that the innovative intermittent feed method could effectively minimize the thermal necrosis to bone tissue caused by the high temperature during the bone drilling. The intermittent feed method could on one hand reduce the drilling temperature and on the other hand shorten the duration of high temperature on the bone. The further study could be to integrate more parameters in the model such as bone density, bone sex or drilling force to complete it and to test the efficiency of intermittent feed method to reduce the thermal necrosis experimentally.
Acknowledgements The work was supported by the State Key Development Program for Basic Research of China (Grant No. 2012CB933200), the State Key Program of National Natural Science Foundation of
China (Grant No. 61534001 and 61372178), National Natural Science Fund for Innovation Group (Grant No. 51221163), and the Fundamental Research Funds for the Central Universities. References [1] A.R. Eriksson, T. Albrektsson, Temperature threshold levels for heat-induced bone tissue injury: a vital-microscopic study in the rabbit, J. Prosthet. Dent. 50 (1) (1983) 101–107. [2] W.R. Krause, Orthogonal bone cutting: saw design and operating characteristics, J. Biomech. Eng. 109 (3) (1987) 263–271. [3] G. Augustin, S. Davila, K. Mihoci, T. Udiljak, D.S. Vedrina, A. Antabak, Thermal osteonecrosis and bone drilling parameters revisited, Arch. Orthop. Trauma Surg. 71 (2008) 77–128. [4] H.C. Shin, Y.S. Yoon, Bone temperature estimation during orthopaedic round bur milling operations, J. Biomech. 33 (2006) 39–139. [5] J.E. Lee, O.B. Ozdoganlar, Y. Rabin, An experimental investigation on thermal exposure during bone drilling, Med. Eng. Phys. 34 (10) (2012) 1510–1520. [6] F. Karaca, B. Aksakal, M. Kom, Influence of orthopaedic drilling parameters on temperature and histopathology of bovine tibia: an in vitro study, Med. Eng. Phys. 33 (10) (2011) 1221–1227. [7] G. Augustin, T. Zigman, S. Davila, T. Udilljak, T. Staroveski, D. Brezak, S. Babic, Cortical bone drilling and thermal osteonecrosis, Clin. Biomech. 27 (4) (2012) 313–325. [8] R.K. Pandey, S.S. Panda, Optimization of multiple quality characteristics in bone drilling using grey relational analysis, J. Orthop. 12 (1) (2015) 39–45. [9] J. Lee, O.B. Ozdoganlar, Y. Rabin, An experimental investigation on thermal exposure during bone drilling, Med. Eng. Phys. 34 (2012) 1510–1520. [10] M.I. Hassan, V. Kalidindi, A. Carner, N. Lemmerman, M.V. Thomas, I.S. Jawahir, K. Saito, Scaling human bone properties with PMMA to optimize drilling conditions during dental implant surgery, in: Progress in Scale Modeling, Springer, Netherlands, 2008, pp. 497–507.
220
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
[11] D. Webster, F.A. Schulte, F.M. Lambers, G. Kuhn, R. Müller, Strain energy density gradients in bone marrow predict osteoblast and osteoclast activity: a finite element study, J. Biomech. 48 (5) (2015) 866–874. [12] S.J. Thiagarajan, R. Yang, C. King, S. Narumanchi, Bubble dynamics and nucleate pool boiling heat transfer on microporous copper surfaces, Int. J. Heat Mass Transf. 89 (2015) 1297–1315. [13] B. Zhou, Y. Liu, S.K. Tan, Efficient simulation of oxygen cutting using a composite heat source model, Int. J. Heat Mass Transf. 57 (1) (2013) 304–311. [14] J.F. Li, L. Li, F.H. Stott, Comparison of volumetric and surface heating sources in the modeling of laser melting of ceramic materials, Int. J. Heat Mass Transf. 47 (6) (2004) 1159–1174. [15] Y.K. Tu, Y.C. Liu, W.J. Yang, L.W. Chen, Y.Y. Hong, Y.C. Chen, L.C. Lin, Temperature rise simulation during a kirschner pin drilling in bone, in: 3rd IEEE International Conference on Bioinformatics and Biomedical Engineering, 2009, pp. 1–4. [16] M.M. Mokhtar, H. Fawad, Bone drilling modeling and simulation techniques, in: IEEE Symposium on Business, Engineering and Industrial Applications (ISBEIA), 2012, pp. 357–361.
[17] J.E. Lee, Y. Rabin, O.B. Ozdoganlar, A new thermal model for bone drilling with applications to orthopaedic surgery, Med. Eng. Phys. 33 (10) (2011) 1234– 1244. [18] S. Sezek, B. Aksaka, F. Karaca, Influence of drill parameters on bone temperature and necrosis: a FEM modeling and in vitro experiments, Comput. Mater. Sci. 60 (2012) 13–18. [19] F. Kong, R. Kovacevic, 3D finite element modeling of the thermally induced residual stress in the hybrid laser/arc welding of lap joint, J. Mater. Process. Technol. 210 (6) (2010) 941–950. [20] Y.U. Huairen, Theory of Metal Cutting (in Chinese) (1990) 100–105. [21] G.T. Smith, Cutting Tool Technology: Industrial Handbook, Springer Science & Business Media, 2008, pp. 96–102. [22] L. Qi, X. Wang, M.Q. Meng, 3D finite element modeling and analysis of dynamic force in bone drilling for orthopedic surgery, Int. J. Numer. Meth. biomed. Eng. 30 (9) (2014) 845–856. [23] B.L. Tai, A.C. Palmisano, B. Belmont, T.A. Irwin, J. Holmes, A.J. Shih, Numerical evaluation of sequential bone drilling strategies based on thermal damage, Med. Eng. Phys. 37 (9) (2015) 855–861.
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Optimization of bone drilling process based on finite element analysis Xiashuang Li a, Wei Zhu a, Junqiang Wang b,⇑, Yuan Deng a,⇑ a b
Beijing Key Laboratory for Advanced Functional Materials and Thin Film Technology, School of Materials Science and Engineering, Beihang University, Beijing 100191, China Department of Orthopaedics and Traumatology, Jishuitan Hospital, Beijing 100035, China
h i g h l i g h t s
g r a p h i c a l a b s t r a c t
Volumetric heat source is loaded to
obtain more accurate simulation results. Multiple parameters’ synergistic effects on the drilling temperature are studied. Empirical formula with sensitivity analysis is proposed to predict the temperature. Intermittent feed drilling is demonstrated to reduce the drilling temperature.
a r t i c l e
i n f o
Article history: Received 1 October 2015 Revised 25 May 2016 Accepted 18 July 2016 Available online 19 July 2016 Keywords: Bone drilling Feed rate Spindle speed Drill diameter Intermittent feed drilling Finite element method
a b s t r a c t Bone drilling is frequently used in modern surgery operations, especially in fracture treatment. Drilling temperature should be strictly controlled, since when the temperature is 47 °C higher over 60 s, irreversible damages happen to bone. The study investigated the heat transfer of bone drilling by a threedimensional model based on finite element method. Three principal drilling parameters were studied: feed rate, spindle speed and drill diameter. A parametric study proved that drilling temperature increases once any of three parameters rises. Parameters’ effects on drilling temperature are proved synergistic. An empirical formula of maximum drilling temperature depending on three drilling parameters is given with a sensitivity analysis. To a drilling process of surgical operations, this formula could be used for temperature prediction and optimization of parameters respecting drilling temperature. Furthermore, the intermittent feed drilling is proposed and demonstrated to be very efficient to reduce thermal necrosis to bone tissue. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction Bone drilling is a significant part of the internal and external fixation processes in orthopedic surgery. The heat generated during drilling may be partially dissipated by the blood and tissue fluids and partially carried away by the chips formed. However, due to the poor thermal conductivity of fresh cortical bone, which is
⇑ Corresponding authors. E-mail addresses: [email protected] (J. Wang), [email protected] (Y. Deng). http://dx.doi.org/10.1016/j.applthermaleng.2016.07.125 1359-4311/Ó 2016 Elsevier Ltd. All rights reserved.
among 0.38–2.3 Wm1 K1, it’s difficult for bone to conduct the rest heat away from the cutting edge. So, the temperature rise at the cutting edge in a deep cortical hole would be extremely high. The phenomenon of bone injury, resulting from the rise of heat during drilling operations, is defined as thermal necrosis. This damage to bone cells would delay the healing process after the surgery and reduce the strength of the fixation. Eriksson and Albrektsson [1] observed that the cortical necrosis and delayed healing occur in living animal bone heated at 47 °C for one minute. Krause [2] states that osteoclasts begin to die when the drilling temperature is higher than 50 °C. According to the study of Augustin [3],
212
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Nomenclature {b} a c c1 c2 c3 c4 d [D] Dep Din Dout e f h Kxx Kyy Kzz l
unit outward normal vector intensity coefficient specific heat [J kg1 K1] first regression coefficient second regression coefficient third regression coefficient fourth regression coefficient diameter of drill [mm] thermal conductivity matrix drilling depth [mm] inside diameter of bone model [mm] outside diameter of bone model [mm] constant feed rate of drilling [mm min1] heat convection coefficient [W m2 K1] elemental thermal conductivity in the direction of x [W m1 K1] elemental thermal conductivity in the direction of y [W m1 K1] elemental thermal conductivity in the direction of z [W m1 K1] length of bone model [mm]
thermal damage was caused to the bone while temperature was between 47 °C and 50 °C. The thermal damage on bone tissue during orthopedic sawing was also investigated and the threshold temperature found was about 55 °C [4]. In order to minimize the damage caused by the high temperature (above 47 °C) in bone drilling operations, it is necessary to optimize the drilling parameters. Many researches have been conducted to find out effects of different drilling parameters on multiple quality characteristics of bone drilling. Drilling parameters such as spindle speed, feed rate, drilling depth, drill tip angle, drilling force and drill diameter were experimentally studied [5–7]. Pandey and Panda [8] used the analysis of variance (ANOVA) to determine the effect of each drilling parameter on drilling temperature. Beside, drilling experiments on bovine bone with spindle speeds from 800 rpm to 3800 rpm were performed by Lee et al. [9], and it is concluded that the drilling temperature rises with increase in spindle speed. Other methods like drilling experiments on bone-like material Poly Methyl Methacrylate (PMMA) were also used for prediction [10]. Previous researches were mainly aimed to investigate the effect of individual drilling parameter on the bone drilling temperature. This paper studies multiple parameters’ synergistic effects on drilling temperature with a three-dimensional model. The objective is to build an instruction for classic bone drilling operation in form of a general empirical formula that is concise and simple for using. Recently the numerical simulation based on finite element method is a possible substitute for high-cost and complex experimental work, especially as a useful tool for verification of analytical results and prediction. ANSYS is powerful and practical FEM-based numerical simulation software for analyzing complicated structure. ANSYS Mechanical is capable to analyze the thermal response of structures to heat transfer effects, involving conduction, convection, and radiation heat transfer [11]. It is widely used for modeling heat transfer problem owing to its excellent performance in solving problem of different modes including steady, transient, linear and nonlinear. Thiagarajan, King et al. [12] investigated a steadystate heat transfer simulation in ANSYS to quantify thermal losses of the pool boiling. Zhou et al. [13] simulated oxygen cutting, similar with procedure of drilling, using a composites heat source model to predict the cutting temperature in ANSYS.
L vector operator n spindle speed of drill [rpm] q_ drilling (x, y, z, t) volumetric heat generation rate [W m3] t time [s] T(x, y, z, t) current temperature of the bone [K] T1 ambition temperature [°C] empirical value of TMAXm issued from the regression Temp [°C] Texp experimental value of TMAXm issued from the model [°C] Thuman human body’s temperature [°C] TS temperature on the surface S [°C] TMAX(t) maximum temperature on the bone at moment t (drilling temperature) [°C] TMAXm maximum value of TMAX during a drilling process [°C] V velocity vector for mass transport of heat [m s1] vc cutting speed of drill [m s1] vf feed rate of drilling [m s1] q density [kg m3]
The Finite Element Method (FEM) presents effective performance not only in virtually all areas of engineering but also in applied science like medical science. The FEM could avoid the high-cost in term of expense and time of experimental study in medical science. For example, ANSYS is used for the simulation of biomaterials like bone and dent to investigate their mechanical and thermal properties. This paper analyzes the temperature’s change during cortical bone drilling with different parameters, including drill spindle speed n, feed rate f and drill diameter d via ANSYS. Specifically, with customizing ANSYS by APDL (ANSYS Parametric Design Language) script, a three-dimensional heat transfer model of bone was built and verified. In addition, a volumetric heating source was used in this model. Li et al. [14] compared volumetric and surface heating sources in the modeling of laser melting of ceramic materials and concluded that the simulation result with volumetric heating source is more accurate and stable. Other applications of finite element method were also used to investigate the mechanism of cutting process and to predict drilling temperature. Tu et al. [15] used a three-dimensional finite element model to show that lowering the initial temperature of Kirschner pin can decrease the temperature rise as well as the size of the thermally damaged zone. Mokhtar and Fawad [16] worked out a complete theory for bone drilling modeling to obtain the heat flux. A new thermal model has been worked out by Lee et al. [17] for bone drilling with a sensitivity analysis and a single parametric study was also investigated for this model. Sezek et al. [18] worked with a FEM-based model for bone drilling and found a safety zone for drilling parameters. To date, most studies of bone drilling modeling proposed a new model or a new theory without further practical research of the drilling temperature by exploring the model or the theory. This paper presents a detailed investigation from the foundation to the validation of the three-dimensional model, followed by results of different levels. Hundreds of simulations were completed in ANSYS with values that correspond with the real drilling procedure. The influence on drilling temperature is studied separately with single parameter and double parameters. The effects of multiple parameters on drilling temperature were proved synergistic. In order to generalize the numerical results, an empirical formula
213
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
is given to predict the drilling temperature in function of these three drilling parameters. At last but not the least, an innovative intermittent feed drilling method is also proposed and proved by this model to effectively reduce the drilling temperature. 2. Materials and methods 2.1. Thermal model ANSYS uses the finite-element method to solve the underlying governing equations and the associated problem-specific boundary conditions. With the presence of internal heat source due to the drilling process, the balance relation of the heat flow of a volume bounded by an arbitrary surface S is as follows [19]:
qc
@Tðx; y; z; tÞ þ fVgT fLgT þ fLgT fqg ¼ q_ drilling ðx; y; z; tÞ @t
ð1Þ
where q is the density, c is the specific heat, and T is the locationtime dependent temperature. {L} is the vector operator, {L} 8@9 < @x = @ = @y ; {V} is the velocity vector for mass transport of heat, {V} :@; 8 @z 9 < vx = = v y . q_ drilling (x, y, z, t) is the volumetric heat generation rate, : ;
vz
caused by the drilling movement and it’s given as follows [20,21]:
v 0:85 q_ drilling ðx; y; z; tÞ ¼ av 0:75 f c where a = 1E10, the coefficient of the intensity;
ð2Þ
vf
f ¼ 60000 , the feed
rate and v c ¼ p the cutting speed. {q} is the heat flux vector, and Fourier’s law is used to relate the heat flux vector to the thermal gradient: n 2 60
d , 2
fqg ¼ ½DfLgT
ð3Þ
where the conductivity matrix [D] is given by: [D] 2 3 K xx 0 0 = 4 0 K yy 0 5. Kxx, Kyy, Kzz are the elemental conductivity 0 0 K zz respectively in the direction of x, y, and z. {V} is the velocity vector for mass transport of heat, based on the proposed model we ignore the effect of the drilling chips carrying away heat on the temperature field of the drilling site, i.e. letting {V} = 0 in (1) and substitute (2) into it, the thermal governing equation of this quasi-steady heat conduction problem becomes:
qc
@Tðx; y; z; tÞ ¼ q_ drilling ðx; y; z; tÞ þ fLgT ½DfLgT @t
ð4Þ
The boundary condition is considered. Specified convection heat loss acting over non-drilled surface S (Newton’s law of cooling) is given as:
fqgT fbg ¼ hðT S T1 Þ
ð5Þ
where {b} is the unit outward normal vector, h is the heat convection coefficient acting on the surfaces of non-drilled bone (h = 3 Wm2 K1), T1 is the ambition temperature considered as 20 °C, TS is the temperature on the surface S, the radiation heat loss is negligible compared to the convection heat loss. The initiative condition is:
Tðx; y; z; t ¼ 0Þ ¼ T human
ð6Þ
Thuman is the average temperature of a human body, meaning 37 °C. In many simulations [15,17] and even experimental investigations [23] of the bone drilling process, 37 °C was set as the initial bone temperature.
2.2. Numerical modeling and loading In ANSYS, the process of the bone drilling was simulated by a model including a hollow cylinder and a mobile thermal source, respectively representing the bone and the cutting edge of drill. ANSYS is indeed proficient in modeling complex structure by introducing the physical model built in Solidworks. However, the drilling happens on very small part of bone so that what the problem is in fact local. Hence, the complex form of bone has very little change the simulation’s result that is ignorable. In many previous investigation of bone drilling problem, a simplified geometry like a hollow cylinder is often used for the bone and they all get satisfied results [12,22]. The APDL mode is used due to its flexibility over the Graphics User Interface (GUI) mode for the simulation.SOLID90 is chosen as the element type that has 20 nodes with a single degree of freedom, temperature, at each node. The 20-node elements have compatible temperature shapes and are well suited to model curved boundaries. It’s possible to set up several materials properties for the element such as, K xx ; K yy ; K zz ; q and c. Besides, this element could carry surface loads like convections and body loads like heat generations. Material properties of the bone used in the model are listed in Table 1. The physical model is a hollow cylinder as shown in Fig. 1(a) and (b) and obeying to a homogenous and anisotropic material whose geometry parameters are given in Table 2. As for the mesh generation, dual mesh density technology was used in order to improve the calculation accuracy and to reduce the processing time. Since the temperature gradient would rise in the drilling region, a higher mesh density was needed for the area near of drilling site, noted S2. While for the rest of model, a low mesh density was enough, noted S1. To define those two densities, an approximation method was applied to find the optimal values that ensure the time of calculate and the stability of results. Four pairs density were tested as follows: A1: A2: A3: A4:
S1 = 0.01 mm, S1 = 0.02 mm, S1 = 0.01 mm, S1 = 0.01 mm,
S2 = 0.001 mm; S2 = 0.001 mm; S2 = 0.0005 mm; S2 = 0.0008 mm.
The maximum drilling temperature calculated for each four groups are given in Table 3. The difference in temperature of these four groups is small enough to be ignored so the difference in calculation time is decisive to choose the density pair. Group A2 has the shortest calculation time but singularity problem appears, as the mesh is not smooth enough due to the large difference of densities in this pair. Finally, the length of meshed elements is chosen as 0.001 mm for high mesh density and 0.01 mm for low mesh density. The mesh of the transition part on the boundary of two regions were refined and smoothed as necessary to create the appearance of a natural transition avoiding singular points. The final meshed model is presented in Fig. 1(c). The convection and the heating source are two loads to be applied. The air convection between the bone and the environment was loaded with the heat convection coefficient of 3 W m2 K1 on all the surfaces of bone. The heat generation due to the drilling was modeled by a volumetric heating source that moves from the exterior to the interior of the bone following the drilling path until the end of drilling depth. To apply this mobile load, the to-be-drilled part is divided in six sections (six cylinder), numbered from 1 to 6, as shown in the Fig. 1(d). So that the volumetric heating source could move from section number 1 to section number 6 successively. At each time increment, this stepping load of heat generation rate is applied on only
214
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
one piece of these six sections to simulate the local heating resulting from the drilling. Then in the next time increment, the volumetric heat source was only applied on the next section. 2.3. Numerical simulation Three principal drilling parameters have been chosen for numerical calculations, same as the situation of the real drilling processes. They are f the feed rate, n the spindle speed, and d the diameter of the drill. The bone temperature is the target during the drilling process, so the maximum temperature of the entire bone is extracted. Numerical simulations were investigated for single parameter study at first and the typical values have been used for each parameter, as shown in the Table 4. Then numerical simulations were done for double parameter study and the test values are listed in Table 5. In the case of the single parameter study, the more values are tested, the more accurate is the result. However, for a double parameter study, it’s not the quantity but the representativeness and equal interval of test values that determine the accuracy of result. Therefore, values for the double parameter study were well selected, as shown in Table 5.
Table 1 Material properties of the bone. Item
Value
Density (kg m3) Specific heat capacity (J kg1 K1) Thermal conductivity (W m1 K1) at 20 °C Initial temperature (°C)
1700 1260 0.38 37
Table 2 Geometry parameters of the physical model. Item
Value
Outside diameter (mm) Inside diameter (mm) Length (mm) Drilling depth (mm)
100 70 100 10
Table 3 Calculation results of different density pairs. Group
TMAXm (°C)
A1 A2 A3 A4
44.012 44.011 44.002 43.985
Table 4 Tested values for single parameter study. Drilling parameters f (mm min n (rpm) d (mm)
1
)
Test values 30–50–60–75–100–120–150 150–300–600–900–1200–1500–1800 2.0–3.2–3.5–4.0–4.5–5.0
Table 5 Test values for double parameter study. Drilling parameters
Test values
f (mm min1) n (rpm) d (mm)
25–50–75–100 300–600–900–1200–1500–1800 2.0–2.5–3.0–3.5–4.0–4.5–5.0
Fig. 1. (a) Cross section of the hollow cylinder; (b) Physical model; (c) Meshed model; (d) Zoom of the part of bone to be drilled.
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Fig. 2. Time dependent drilling temperature with different feeding rates (n = 900 rpm, d = 3.2 mm).
3. Results and discussions 3.1. Single parameter results Firstly, the feed rate f has been chosen as the only variable for numerical calculations. The values of f used are in Table 4. The spindle speed n and the drill diameter d was kept constant: n = 900 rpm, d = 3.2 mm. Fig.2 displays the change of the timedependent maximum drilling temperature corresponding to different values of f. Initially the temperature was set at 37 °C. At each moment t, the maximum temperature of the entire bone is noted as TMAX, which is the drilling temperature, represented by one point in the Fig. 2. The time of drilling is shorter when the feed rate f increases for the same drilling depth. In general, the TMAX increases rapidly at first and becomes almost stable until the end of drilling process. In fact, the maximum temperature of the bone appears at the cutting edge and rises rapidly in the first several seconds. Gradually the drilled part cools down and the maximum temperature of bone always appears at the point where the drilling is going on. The value of the maximum temperature of bone in each drilling repeated keeps almost stable. For each trial, the TMAX reaches its maximum value, noted as TMAXm, in the first two or three seconds. In Fig. 2, the drilling temperature TMAX increases with f increasing. Keeping n = 900 rpm, d = 3.2 mm, TMAXm goes over
215
47 °C when f is higher than 75 mm min1. Fig. 2 indicates that f has a considerable influence on drilling temperature. When f goes from 30 mm min1 to 150 mm min1, the TMAXm increases from 40.63 °C to 51.55 °C, which means a relative increase of 26.85%. In a real case, when the drill is fed more rapidly, the drill force is stronger and so is the friction force. This phenomenon will generate more heat on the friction area in the same period of time. Same kind of calculation was conducted for different d and different n whose results are given in Fig. 3. As same as the situation for f, the TMAX of each trial reached its maximum value in the first two or three seconds of drilling process. The larger the parameter’s (n or d) value is, the longer it takes to reach the TMAXm. TMAXm rises with n and d in different speed. The relative increase is about respectively 14.45% for d going from 2.0 mm to 5.0 mm and 28.94% for n going from 150 rpm to 1800 rpm. In practice, when the drill’s diameter becomes larger, so does the friction area between the drill and the bone tissue, so more heat is generated. The influence of the increment of n or f on the increase of TMAXm increases is bigger than that of d. In reality, more heat is generated during the drilling process due to the increase of each of three parameters. However, the larger the drill diameter (d) is, the larger the diameter of the hole is, so that there is more surface for heat to transmit. In this way, the drilling temperature will not rise as much as that in the other two cases. 3D diagrams of TMAX in function of one single drilling parameter and time show the growth of drilling temperature more directly. Fig. 4(a) shows the change of TMAX with n and time. It manifests that for n higher than 1700 rpm, 1.5 s after the beginning of drilling, TMAX exceeds the security temperature (47 °C). It tells the maximum continuous time the drilling before causing thermal necrosis for a value of n given. It could also be used to find out possible values of n when the drilling time is set. Fig. 4(b) is the projection of this 3D diagram, from which statistics can be gathered more easily. Same kind of 3D diagram was also completed for f and time or d and time. 3.2. Double parameters results Furthermore, in order to evaluate three parameters’ different level of influence and their synergistic effects on drilling temperature, two groups of double parameter were considered as covariables: d-f and n-f. Tested values for each parameter are listed in Table 4. Keeping n constant (n = 900 rpm), for each possible value of d, numerical experiments were conducted with different f. In this part, 4 7 simulations were investigated and the results are presented in Fig. 5.
Fig. 3. Change of drilling temperature with (a) different d, f = 75 mm min1, n = 900 rpm; (b) different n, f = 75 mm min1, d = 3.2 mm.
216
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Fig. 4. (a) 3D diagram of TMAX in function of n and time; (b) Projection of 3D diagram of TMAX in function of n and time.
Fig. 5. (a) Different effects of d and f on drilling temperature; (b) Different effects of f and d on drilling temperature; (c) 3D diagram of TMAXm in function of d and f; (d) Projection of 3D diagram.
In Fig.5(a), each curve represents the change of TMAXm in function of f with different given values d. As mentioned in 3.1, for each bone drilling simulation, TMAXm represents the maximum value of TMAX during the whole drilling process. The larger the d is, the
more TMAXm grows with f. For each f given, the change of TMAXm in function of d is shown vertically: The higher the feed rate f is, the faster TMAXm increases with d. This means that the effects of d and f on the drilling temperature are synergistic even though
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
217
Fig. 6. (a) n and f on drilling temperature; (b) Different effects of f and n on drilling temperature; (c) 3D diagram of TMAXm in function of n and f; (d) Projection of 3D diagram.
they are independent from one to another. In practice, when the drill feeding is faster, more heat will be generated in the drilling process with a large drill diameter than that of a small drill diameter so that the rise of drilling temperature is bigger. When the drill diameter is larger, a larger feed rate will generate much more heat on the increased friction area than a small feed rate and consequently the drilling temperature increases more. In Fig.5(b), each curve indicates the change of TMAXm along with d for f disposed. It’s a different way to interpret the same data as in Fig. 5(a). These two diagrams prove that d = 3.5 mm is a safe drill diameter for n = 900 rpm and f 6 100 mm min1 and f = 50 mm min1 is a safe feed rate for n = 900 rpm and d 6 5 mm. Fig. 5(c) is the 3D diagram of TMAXm in function of d and f. This kind of diagram helps visualize the security zone (below 47 °C) bounded by different values of f and n, for a drill diameter d given. In this situation, non-security zone (over 47 °C) is colored in red1 and orange. In Fig. 5(d), the 3D diagram is projected to the plan of d and f, which makes it less intuitionistic but more convenient to read values. Same type of study was implemented keeping d constant (d = 3.2 mm) and variable n and f. 4 6 simulations were investigated in total for this part and the results are shown in Fig. 6. In a similar way, Fig. 6(a) and (b) confirm the synergistic effects of n and f on the drilling temperature: TMAXm increases more rapidly 1
For interpretation of color in Figs. 5 and 7, the reader is referred to the web version of this article.
with n when the feed rate f is larger, and vice versa. From these two diagrams it can be told that n = 1200 rpm is a safe drill diameter for d = 3.2 mm and f 6 100 mm min1 and f = 50 mm min1 is a safe feed rate for d = 3.2 mm and n 6 1800 rpm. Fig. 6(c) is the 3D diagram of TMAXm in function of n and f. The boundary between the security zone and the non-security zone is the one between the green band and the yellow band. Fig. 6(d) shows a projection of the 3D diagram. 3.3. Empirical formula In order to quantify the effects of parameters on the drilling temperature, an empirical formula of drilling temperature was built in MATLAB. The model was linearized at first considering these three drilling parameters as predictor variables. Then a multiple linear regression was realized with the function regress. Supposing the formula of TMAX in the following type: c2 c3
TMAX m ¼ ec1 f d nc4
ð7Þ
To linearize the formula, logarithm natural is taken at both sides, which makes:
ln ðTMAX m Þ ¼ c1 þ c2 ln ðf Þ þ c3 ln ðdÞ þ c4 ln ðnÞ
ð8Þ
52 groups of numerical calculations results were entered in MATLAB to determine the constants in the formula. Finally, the empirical formula obtained is:
218
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
Fig. 7. (a) Residual case order plots; (b) Comparison between TMAXm experimental (Texp) and TMAXm empirical (Temp).
TMAX m ¼ e2:4812 f
0:1092 0:1387 0:0995
d
n
ð9Þ
A local sensitivity analysis is done to figure out how the uncertainty of the TMAXm can be distributed to different sources of uncertainty of the three parameters. The sensitivity of this empirical formula with changes of each parameter is equal to the partial derivative of TMAXm respectively to each. Obviously, the results correspond to the exponent number of each parameter, meaning 0.1092 for f, 0.1387 for d and 0.0995 for n. The empirical formula of maximum drilling temperature shows that the TMAXm increases non-linearly with the change of those parameters. This analysis also illustrates that in these current units of parameters, TMAXm is more sensitive to the drill diameter than to the feed rate and least sensitive to the drill spindle speed. During bone drilling, the operator need to be more conscious with a large drill diameter since the drilling temperature will arise faster in this situation with the augmentation of other drilling parameters. For a disposed bone drilling, the feed rate should be decreased firstly if the predicted drilling temperature is critical. Fig.7(a) displays an error bar plot of the confidence intervals on the residuals from this regression. The interval around the 29th, 51th and 52th residual (shown in red) does not contain zero. This indicates that the residual is larger than expected in 95% of new observations, and suggests the data point is an outlier. A comparison between TMAXm experimental (Texp) and TMAXm empirical (Temp) demonstrates that the empirical formula is reliable enough, as shown in Fig.7(b). As a concise and compact representation of the numerical results, this reliable empirical formula quantifies the influence of drilling parameters on maximum drilling temperature. Before an orthopedic surgery with a bone drilling procedure, by entering the values of these three drilling parameters, a prediction of maximum drilling temperature could be made. If this temperature exceeds 47 °C, the safe threshold, the empirical formula could tell which parameter has the most significant influence on drilling temperature. And the value of this parameter is mathematically be diminished in advance.
drilling actions. Simultaneously, the convection with air in the environment helps the bone to cool down. To complete the same work of drilling, same quantity of heat is generated but attributing to the intermittent feed, heat is better conducted or scattered in this case. The drilling temperature is effectively reduced as the result of simulation showed with the FEM based model in ANSYS. In this paper, the number of passages means how many steps the drilling is divided into and rest time is the time between two successive steps of drilling. Four groups of numerical experiments were carried out with same drilling parameters f = 75 mm min1, d = 3.2 mm and n = 1500 rpm. In the first group, the group of reference, the drilling process was continuous. In the three rest groups, the drilling process was respectively divided into 2, 6 and 12 passages with different rest time. The change of TMAXm in time of each simulation is given in Fig.8. In all the diagrams, the curves with squares represent the same one continuous process. The curve with disks and the curve with triangles in the same diagram are processes of intermittent feed with same number of passages but different rest time. For example, in Fig. 8(a), the curve with disks represents a process in 2 passages with a rest time of 2 s and the curve with triangles represents a process in 2 passages with a rest time of 5 s. When the drilling process is divided into 2 or 6 passages, the drilling temperatures did not decreased, regardless of the length of rest time. When the number of passages reaches 12 as presented in Fig. 8(c), the drilling temperature is obviously reduced: from 48.5 °C that is over the security temperature to below the security temperature 47 °C. By all means, even when the passages number is not large enough to decrease the drilling temperature as shown in Fig. 8 (b), the effective duration in the case of the drilling temperature is over 47 °C is reduced. The bone remained above 47 °C for 12 s with a continuous feed while for just 1 or 2 s, 6 times with an intermittent feed in 6 passages. Hence, the intermittent feed drilling is very effective to lower the irreversible damage on human body resulted by the high temperature during surgical operations.
4. Conclusions 3.4. Intermittent feed drilling An innovative method of feed, intermittent feed is discussed in this part to decrease the drilling temperature. Intermittent feed reduces the continuous drilling time as well as the continuous heating time on bone. This leaves the bone more time to conduct the generated heat away from the cutting edge between two
In this paper, a three-dimensional heat transfer model based on finite element method was built with a volumetric heating source to investigate the thermal process of bone drilling. The study of this model showed that maximum drilling temperature during the drilling process increases with each of these three parameters including feed rate f, drill spindle speed n and drill diameter d.
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
219
Fig. 8. Change of TMAXm of drilling in (a) 2 passages; (b) 6 passages; (c) 12 passages with 2 s or 5 s of rest between passages.
Moreover, their effects on the drilling temperature are synergistic. For the same increase of f, the increase of maximum drilling temperature is higher when n or d is higher. For the same increase of n or d, the increase of maximum drilling temperature is higher when f is higher. Based on numerical results from this model, an empirical formula of maximum drilling temperature that depends on three different drilling parameters was provided to quantify the effects of parameters and to predict the maximum temperature for these three parameters given. The simulations with this model proved that the innovative intermittent feed method could effectively minimize the thermal necrosis to bone tissue caused by the high temperature during the bone drilling. The intermittent feed method could on one hand reduce the drilling temperature and on the other hand shorten the duration of high temperature on the bone. The further study could be to integrate more parameters in the model such as bone density, bone sex or drilling force to complete it and to test the efficiency of intermittent feed method to reduce the thermal necrosis experimentally.
Acknowledgements The work was supported by the State Key Development Program for Basic Research of China (Grant No. 2012CB933200), the State Key Program of National Natural Science Foundation of
China (Grant No. 61534001 and 61372178), National Natural Science Fund for Innovation Group (Grant No. 51221163), and the Fundamental Research Funds for the Central Universities. References [1] A.R. Eriksson, T. Albrektsson, Temperature threshold levels for heat-induced bone tissue injury: a vital-microscopic study in the rabbit, J. Prosthet. Dent. 50 (1) (1983) 101–107. [2] W.R. Krause, Orthogonal bone cutting: saw design and operating characteristics, J. Biomech. Eng. 109 (3) (1987) 263–271. [3] G. Augustin, S. Davila, K. Mihoci, T. Udiljak, D.S. Vedrina, A. Antabak, Thermal osteonecrosis and bone drilling parameters revisited, Arch. Orthop. Trauma Surg. 71 (2008) 77–128. [4] H.C. Shin, Y.S. Yoon, Bone temperature estimation during orthopaedic round bur milling operations, J. Biomech. 33 (2006) 39–139. [5] J.E. Lee, O.B. Ozdoganlar, Y. Rabin, An experimental investigation on thermal exposure during bone drilling, Med. Eng. Phys. 34 (10) (2012) 1510–1520. [6] F. Karaca, B. Aksakal, M. Kom, Influence of orthopaedic drilling parameters on temperature and histopathology of bovine tibia: an in vitro study, Med. Eng. Phys. 33 (10) (2011) 1221–1227. [7] G. Augustin, T. Zigman, S. Davila, T. Udilljak, T. Staroveski, D. Brezak, S. Babic, Cortical bone drilling and thermal osteonecrosis, Clin. Biomech. 27 (4) (2012) 313–325. [8] R.K. Pandey, S.S. Panda, Optimization of multiple quality characteristics in bone drilling using grey relational analysis, J. Orthop. 12 (1) (2015) 39–45. [9] J. Lee, O.B. Ozdoganlar, Y. Rabin, An experimental investigation on thermal exposure during bone drilling, Med. Eng. Phys. 34 (2012) 1510–1520. [10] M.I. Hassan, V. Kalidindi, A. Carner, N. Lemmerman, M.V. Thomas, I.S. Jawahir, K. Saito, Scaling human bone properties with PMMA to optimize drilling conditions during dental implant surgery, in: Progress in Scale Modeling, Springer, Netherlands, 2008, pp. 497–507.
220
X. Li et al. / Applied Thermal Engineering 108 (2016) 211–220
[11] D. Webster, F.A. Schulte, F.M. Lambers, G. Kuhn, R. Müller, Strain energy density gradients in bone marrow predict osteoblast and osteoclast activity: a finite element study, J. Biomech. 48 (5) (2015) 866–874. [12] S.J. Thiagarajan, R. Yang, C. King, S. Narumanchi, Bubble dynamics and nucleate pool boiling heat transfer on microporous copper surfaces, Int. J. Heat Mass Transf. 89 (2015) 1297–1315. [13] B. Zhou, Y. Liu, S.K. Tan, Efficient simulation of oxygen cutting using a composite heat source model, Int. J. Heat Mass Transf. 57 (1) (2013) 304–311. [14] J.F. Li, L. Li, F.H. Stott, Comparison of volumetric and surface heating sources in the modeling of laser melting of ceramic materials, Int. J. Heat Mass Transf. 47 (6) (2004) 1159–1174. [15] Y.K. Tu, Y.C. Liu, W.J. Yang, L.W. Chen, Y.Y. Hong, Y.C. Chen, L.C. Lin, Temperature rise simulation during a kirschner pin drilling in bone, in: 3rd IEEE International Conference on Bioinformatics and Biomedical Engineering, 2009, pp. 1–4. [16] M.M. Mokhtar, H. Fawad, Bone drilling modeling and simulation techniques, in: IEEE Symposium on Business, Engineering and Industrial Applications (ISBEIA), 2012, pp. 357–361.
[17] J.E. Lee, Y. Rabin, O.B. Ozdoganlar, A new thermal model for bone drilling with applications to orthopaedic surgery, Med. Eng. Phys. 33 (10) (2011) 1234– 1244. [18] S. Sezek, B. Aksaka, F. Karaca, Influence of drill parameters on bone temperature and necrosis: a FEM modeling and in vitro experiments, Comput. Mater. Sci. 60 (2012) 13–18. [19] F. Kong, R. Kovacevic, 3D finite element modeling of the thermally induced residual stress in the hybrid laser/arc welding of lap joint, J. Mater. Process. Technol. 210 (6) (2010) 941–950. [20] Y.U. Huairen, Theory of Metal Cutting (in Chinese) (1990) 100–105. [21] G.T. Smith, Cutting Tool Technology: Industrial Handbook, Springer Science & Business Media, 2008, pp. 96–102. [22] L. Qi, X. Wang, M.Q. Meng, 3D finite element modeling and analysis of dynamic force in bone drilling for orthopedic surgery, Int. J. Numer. Meth. biomed. Eng. 30 (9) (2014) 845–856. [23] B.L. Tai, A.C. Palmisano, B. Belmont, T.A. Irwin, J. Holmes, A.J. Shih, Numerical evaluation of sequential bone drilling strategies based on thermal damage, Med. Eng. Phys. 37 (9) (2015) 855–861.