- Email: [email protected]
An inverse method to reconstruct the heat flux produced by bone grinding tools
International Journal of Thermal Sciences 101 (2016) 85e92
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
An inverse method to reconstruct the heat flux produced by bone grinding tools Guangjun Wang a, b, Lihui Zhang c, Xudong Wang a, Bruce L. Tai d, * a
School of Power Engineering, Chongqing University, Chongqing 400044, PR China Key Laboratory of Low-grade Energy Utilization Technologies and Systems, Ministry of Education, Chongqing University, Chongqing 400044, PR China c Key Laboratory of Special Purpose Equipment and Advanced Processing Technology, Zhejiang University of Technology, Ministry of Education, Hangzhou, 310014, PR China d Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843, USA b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 20 August 2015 Received in revised form 29 September 2015 Accepted 15 October 2015 Available online xxx
Surgical bone grinding using high-speed, small spherical abrasive tools could cause thermal damage to the surrounding neurovascular structures. Our prior work developed a thermal model based on the twodimensional grinding theory to calculate the temperature distribution. To verify the grinding-theorybased model, this study used experimental data incorporated with an inverse heat transfer method to mathematically estimate the heat flux distribution. This inverse method also considered a time-varying system to reflect a temporal change of the heat flux. Specifically, a coupled approach combining sequential function specification method (SFSM) and sequential quadratic programming (SQP) was employed to calculate the temporal and spatial variables simultaneously. Numerical tests were performed to determine the effectiveness and the limitations of this method. Then, the experimental data of prior work was applied to reconstruct the heat flux. The results verified a nearly time-invariant heat flux and demonstrated a consistent trend in the spatial distribution. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Bone grinding Inverse heat transfer method Finite element method
1. Introduction Bone grinding with high-speed abrasive tools is often used in orthopedic surgery, neurosurgery, otolaryngology, and plastic surgery. Heat is of concern in these procedures since grinding is known to produce more energy compared to other cutting processes, consequently leading to a greater risk to adjacent living tissues. Bone necrosis can occur at about 50 C [1,2]; neural tissues are more vulnerable to the elevated temperature with a threshold level of 43 C [3]. In a minimally invasive skull base surgery, for example, neurosurgeons use a small grinding bit to remove cranial bone surrounded by optic nerves and internal carotid arteries to provide exposure for tumor removal. Excessive heat can cause catastrophic consequences such as stroke and vision damage [4]. Thermal issues in bone cutting processes have been studied extensively in the fields of dentistry, orthopedics and trauma [5e7]. Their results provided the knowledge in tool selection, proper operating parameters, and the level of heat production in the
* Corresponding author. E-mail address: [email protected] (B.L. Tai). http://dx.doi.org/10.1016/j.ijthermalsci.2015.10.021 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.
process. However, there is limited literature investigating bone grinding with abrasive tools. Kondo et al. [8] performed grinding experiments on human cranial bones and recorded temperature with the thermography. Our prior work in collaboration with Shih et al. [4,9,10] conducted grinding experiments on bovine bones and used finite element method (FEM) for estimating temperature distribution. The heat flux distribution was assumed time-invariant in FEM and the spatial distribution was built based on the 2D grinding theory [11]. Despite a decent agreement between FEM and experiments, the heat flux distribution was not completely verified. This study, therefore, is aimed to develop an inverse heat transfer method (IHTM) to solve the heat flux distribution based on experimental inputs and compare to the prior work [4]. In essence, IHTM estimates the unknown thermal boundary conditions based on measurable temperature information [12e16]. The current algorithms for solving IHTM include regularization method, sequential function specification method (SFSM), gradient-based optimization methods and so on. Among them, SFSM is useful in the transient heat transfer problem [17], provided that it can obtain stable inverse solution when a proper number of the future time steps is selected. Samadi et al. [18] applied SFSM to estimate transient heat flux imposed on the rake face of a cutting
86
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
tool during a cutting operation. Behbahani-nia et al. [19] combined dual reciprocity boundary element and SFSM to estimate the time and space-varying boundary heat flux in a two-dimensional configuration. However, when applying SFSM to a moving heat source scenario with unknown magnitude and spatial distribution, the sensitivity coefficients of measured temperature with respect to each elemental heat flux component need to be re-calculated, because they change over time due to the change of sensor position relative to the moving heat source. To overcome this issue, this paper presents a hybrid method using both SFSM and sequential quadratic programming (SQP). The total heat flux distribution is considered to consist of two parts f(x0 , y0 ) and qðtÞ. The former describes the spatial distribution of heat flux; the latter reflects the variation of average heat flux over time. SFSM and SQP are then applied to solve f(x0 , y0 ) and qðtÞ, respectively. This paper is organized as follows. The grinding-theory-based model from our prior work is first reviewed and then followed by the mathematical model used for IHTM. A series of numerical tests are conducted to verify the feasibility and identify limitations of the proposed method, including the effects of initial values and the number of future time steps in SFSM. Finally, this method is applied to the experimental data to estimate the heat flux distribution and compare to the grinding-theory-based model. 2. Review of the grinding theory-based model Fig. 1 shows a schematic view of the grinding process on a bone surface. The small spherical grinding tool with diameter D (typically 2e4 mm) and rotational speed u moves across the bone surface along the feeding direction (x-axis) at the speed v and the cutting depth ac. A tilt angle a is defined as the angle between the tangent to the surface and the axis of the grinding tool. In the grinding theory-based model, the spherical tool is divided into a finite amount of elemental grinding wheels (EGW). Each wheel performs grinding individually with its own tangential speed (re u) and depth of cut (ae). The averaged heat flux under an EGW (qe) then can be calculated by
qe ¼
Kðre uÞae Ae
(1)
Fig. 2 shows the finite element model (FE-model) to which the grinding model is applied. The ground groove is located in the middle of the top surface. The heat flux (shaded region) moves along this groove during the grinding process. TCi represents the ith temperature sensor placed on both sides of the grinding groove. Since this is a pure heat transfer model without considering material removal, the groove is built in the model and simplified as an existing slot with the width of lb. This width is defined by the geometrical intersection between the spherical tool and bone surface at point p, as shown in Fig. 1. The number of EGWs in contact with the surface is determined by the EGW size and la, which is the horizontal distance between points p and q. As the example shown in Fig. 2, there are four active EGWs and each one has two elements underneath. The relative strength of heat flux (i.e., with an undefined K) along the feed direction is calculated based on Eq. (1), and across the feed direction is fitted by a triangular distribution of ratio 1 to 0.33. Finally, the experimentally measured temperatures at sensors are used to adjust the FE-model to find an optimal K in Eq. (1), and complete constructing the heat flux. In this model, the spatial heat flux distribution is determined based on the EGW concept and the assumption of triangular distribution. Also, the scale factor K is assumed a constant throughout the process for a time-invariant system. To further justify these assumptions, a mathematical model, along with a time-varying setting, is used to reconstruct the spatial and temporal distributions, as detailed in the following sections. 3. Mathematical heat flux model Rather than using EGWs to map the spatial heat flux distribution to the FE-model, the mathematical model describes the heat flux by two geometrical variables, kx and ky. Assuming a moving coordinate system x0 ey0, the coordinate conversion can be described by
x0 ¼ x þ vt and y0 ¼ y þ ly 2 þ lb =2
(2)
where v is a constant feed rate and t is time. Then, the heat flux field on the thermal load surface is defined using kx and ky by
qðx0 ; y0 ; tÞ ¼ kx x0 þ ky y0 þ cðtÞ
(3)
where Ae is the area where the heat is applied under the EGW and K is a scale factor. The averaged heat flux qe is distributed in a triangular shape on the workpiece surface as shown in Fig. 1. The shape of triangular distribution depends on the grinding direction [11]; the feed-in edge has a thicker chip load and thus produces more heat.
where c(t) is the magnitude of heat flux at time t and location (x0 , y0 ) ¼ (0,0). The geometrical variables are assumed constant and time invariant in this setting. The temporal term c(t) alters the heat flux magnitude during the grinding process to reflect the environmental noise, inhomogeneous material, or the heat sink effect of the tool. At an instantaneous time, the heat flux field q(x0 , y0, t) can be illustrated by a plane as shown in Fig. 3. Although a high-order
Fig. 1. Schematic view of the bone grinding process with spherical grinding tool.
Fig. 2. Finite element thermal model of bone grinding.
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
87
Fig. 4. Temperature response to a discretized heat flux input. Fig. 3. Schematic view of the spatial distribution of heat flux over rectangular grinding zone.
polynomial or power law may be used to form a more complex shape, they will increase the number of unknown variables to solve in the following inverse method. The total grinding heat Q(t) can be obtained by integrating the heat flux field over the thermal load surface:
Zla Zlb Q ðtÞ ¼ 0
qðx0 ; y0 ; tÞdy0 dx0
(4)
0
By substituting Eq. (3) into Eq. (4), c(t) is obtained as follows:
cðtÞ ¼ qðtÞ 0:5$kx la 0:5$ky lb
(5)
where qðtÞ ¼ Q ðtÞ=la lb is the average value of heat flux in the grinding zone at time t. By substituting Eq. (5) into Eq. (3), Eq. (3) can be modified as follows: 0
0
0
0
qðx ; y ; tÞ ¼ f ðx ; y Þ þ qðtÞ
(6)
where f ðx0 ; y0 Þ ¼ kx x0 þ ky y0 0:5$kx la 0:5$ky lb . Assume k ¼ [kx ky], k and qðtÞ can uniquely determine the heat flux distribution at any given time for a direct heat transfer problem. To map this mathematical model to the FE-model (Fig. 2), the heat flux in the geometric center of each element is taken as an average to be uniformly applied to the corresponding element. The temperature Ti(j) at a given location TCi and time step j can be obtained by
Ti ðjÞ ¼ Ti ðjÞjq¼0 þ Ti ðjÞjf ¼0 ði ¼ 1; 2; …; G;
j ¼ 1; 2; …; MÞ
(7)
where G is the number of temperature sensors and M is the total time steps; Ti ðjÞjq¼0 represents the response temperature Ti(j) when only f(x0 , y0 ) is applied during the time interval [0, j]; Ti ðjÞjf ¼0 represents the response temperature Ti(j) when only qðtÞ is applied during the time interval [0, j]. Note that time step j stands for a total time of t ¼ j Dt, where Dt is the time interval between steps in FEM. The Ti ðjÞjq¼0 is numerically solved by FEM with the boundary condition q(x0 , y0, t) ¼ f(x0 , y0 ) and a pre-defined initial temperature T0 and k. The Ti ðjÞjf ¼0 is described by a linear time-invariant system. The discretized heat flux qðtÞ is expressed by qðnÞ. Each component in qðnÞ represents the average heat flux within a time interval Dt at step n, as shown in Fig. 4. The temperature response can be seen as an accumulated result of all heat flux components qðnÞ at different time steps and positions, such that
Ti ðjÞjf ¼0 ¼
j X n¼1
Tin ðjÞ ¼
j X n¼1
Ii ðn; j nÞqðnÞ
(8a)
where Tin ðjÞ is the temperature response to an individual heat flux component qðnÞ, and Ii(n, j n) is the sensitivity coefficient. Because of a moving heat source, the first term of Ii(n, j n) means the position of the heat source at time step n; the second term is the time interval between the occurrence of heat input and the temperature. In other words, Ii(n, j n) can be seen as a temperature response at time step j, as shown in Fig. 4, when a heat flux q0 ¼ 1 [W/mm2] is applied at time step n. Note that time steps n and j are equivalent but used in different time domains. Ii(n, j n) can be mathematically expressed by Eq. (8b), but it is calculated as a discrete function using FEM by placing q0 at different times and locations.
Ii ðn; j nÞ ¼
vTin ðjÞ vqðnÞ
(8b)
4. Inverse heat transfer method The average heat flux qðnÞ and the geometrical variables k are the unknowns in this inverse problem. An inverse method coupling SFSM and SQP is introduced in this section to solve this problem. 4.1. Estimation of qðnÞ using SFSM SFSM [16] is applied to estimate the time-varying qðnÞ with a given k. At any time step N, the components from qð1Þ to qðN 1Þ are known. To estimate qðNÞ for the step interval [N 1, N], the temperature measurements of different sensors within the time period [N 1, N þ r 1] are used, where r is a pre-selected number of future time steps. For stability, it is temporarily assumed that qðnÞ is constant during the period of [N 1, N þ r 1]. Then, the objective function Jq for estimating qðNÞ is defined as the sum of squared error function,
Jq ðqðNÞÞ ¼
G Nþr1 X X i¼1
Timea ðjÞ TiFEM ðjÞ
2
(9)
j¼N
where Timea and TiFEM are the measured and FEM-calculated temperatures at sensor TCi. For j ¼ N, N þ 1, …, N þ r 1, the future temperature TiFEM ðjÞ is
TiFEM ðjÞ ¼ TiFEM ðjÞjq¼0 þ
N1 X n¼1
Ii ðn; j nÞqðnÞ þ
j X
! Ii ðn; j nÞqðNÞ
n¼N
(10) Then, by substituting Eq. (10) to Eq. (9), differentiating Eq. (9) with respect to qðNÞ, and setting it equal zero, the heat flux at the current time step N is given by:
88
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
PG
PNþr1
i¼1
j¼N
Timea ðjÞ TiFEM ðjÞ
qðNÞ ¼
q¼0
PN1 n¼1
! Ii ðn; j nÞqðnÞ $
Pj
n¼N Ii ðn; j
nÞ
!2
(11)
Since the SFSM method requires no iteration and involvement of FEM, it has high computational efficiency; however, selection of a proper future time step r could influence the result. The r is selected based on the sensitivity of the heat transfer system (i.e., thermal diffusivity and sensor location) and the noise level of measured data. In general, a rapid response system requires a smaller r; a high noise system needs a larger r. The r could also be a function of time if the relative position between sensors and heat source varies. For simplicity, a constant r is used here for analysis.
(m K), and initial temperature 20 C. The total number of discrete elements along grinding groove was 52 2 with the element dimensions of 0.38 mm 1.2 mm. Locations of the temperature sensors in the experiment are listed in Table 1 based on the configuration in Fig. 2. A heat flux qðtÞ with a rapid temporal change was artificially created for testing an extreme case:
PG
i¼1
PNþr1 n¼N
Pj
n¼N Ii ðn; j
nÞ
4.2. Estimation of f(x0 , y0 ) using SQP method The basic concept of SQP is to obtain an approximate solution by solving a sequence of quadratic programming sub-problems. This algorithm has global and superliner convergence [20]. Using k ¼ [kx ky] as independent variables, the objective function F(k) is defined as the sum of squared errors between the measured temperature Timea and calculated temperature TiFEM of all temperature sensors over the whole time domain, such that
FðkÞ ¼
G X M X
Timea ðjÞ TiFEM ðjÞ
i¼1 j¼1
2 k;qðnÞ
qðt Þ ¼
(12)
0:2; 0:3;
t2ð0; M=4Þ∪ð3M=4; M Þ t2½M=4; 3M=4
h
i . W mm2
where M is the total amount of time steps. Taking kx ¼ 0.05 W/mm3 and ky ¼ 0.03 W/mm3, the corresponding spatial distribution function of heat source was given as:
f ðx0 ; y0 Þ ¼ 0:05x0 0:03y0 þ 0:006
h
W=mm2
i
(13b)
According to the heat flux artificially generated in Eqs. (13a) and (13b), this “exact” temperature Tiexa at given sensors was obtained by solving a direct heat transfer problem. Then, Tiexa was used to
where TiFEM is calculated with pre-defined qðnÞ (from Section 4.1) and k using FEM. The commercial software, MATLAB, has fmincon function in its Optimization Toolbox, which can exploit SQP for identifying optimal k to minimize Eq. (12). 4.3. Coupling of SFSM and SQP Coupling of SFSM and SQP is formulated for simultaneous estimation of k and qðnÞ, as shown in the flowchart in Fig. 5. In the coupled algorithm, the optimization of k using SQP is considered to be mainline, and the inverse estimation of qðnÞ is treated as the subroutine. During an iteration, k is fixed through the process while qðnÞ is optimized using SFSM to minimize the objection function, Jq (Eq. (9)). Then, with an optimized qðnÞ, SQP is applied to the objective function F(k) (Eq. (12)) to determine a new set of k for the next iteration. The iteration continues until k or F(k) are converged to the step interval less than d (¼104). These stopping criteria are tuned to be two magnitudes higher than the default settings (106), based on trial and error, considering the noise level and the number of inputs. 5. Numerical tests of IHTM A series of numerical tests were performed to measure the robustness of the inverse method and to identify the effects of initial values of k and future time step r. The following parameters were used to setup the FE-model in accordance with our prior work [4]: D ¼ 4 mm, u ¼ 60,000 rpm, a ¼ 30 , v ¼ 20 mm/min, ac ¼ 0.4 mm, lx ¼ 20 mm, ly ¼ 20 mm, lz ¼ 6 mm, the bone density 2050 kg/m3, specific heat 516 J/(kg K), thermal conductivity 0.55 W/
(13a)
Fig. 5. The flowchart of SFSM and SQP coupling strategy.
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
89
Table 1 Locations of the temperature sensors in experiment (unit: mm). Sensors
TC1
TC2
TC3
TC4
(x, y, z)
(8.5, 6.7, 6.0)
(14.2, 6.9, 6.0)
(8.1, 13.5, 6.0)
(13.5, 13.5, 6.0)
generate the simulated temperature measurements using Eq. (14) with randomized errors:
Timea ¼ Tiexa þ ms
(14)
where m is a random variable over the range 2.58 m 2.58 following a normal distribution with zero mean and s is the standard deviation. The thermocouples in the experiments (36 gauge Ktype, Omega Engineering) had a maximum variation of 1.1 C according to the data sheet, and thus s was set as 0.55 C. With a normal distribution setting, 95% of the randomized errors were within 1.1 C. 5.1. Effects of initial values of k The first task was to study the effects of different initial k values on the final results. The r was set as 3 and s was 0 (no noise) for a controlled condition. The calculated results for k and qðnÞ are shown in Table 2 and Fig. 6, respectively. It is found that both inverse solutions are not sensitive to the initial values of k. Errors of kx range from 2% to 8% and ky remains 0.0%. The results of qðnÞ are fairly consistent. However, qðnÞ can hardly reproduce a rapid profile change due to the use of future time steps for stability, which will be discussed in Section 5.2. Compared to ky, a relatively large error of kx is due to its smaller sensitivity. This sensitivity coefficient (S) can be defined by the temperature difference caused by a small change in geometrical variables, such that
Sx;i ðjÞ ¼
Sy;i ðjÞ ¼
Ti ðjÞjkx Ti ðjÞjkx Dkx Dkx
6. Heat flux reconstruction based on experimental data
Ti ðjÞjky Ti ðjÞjky Dky
(15b)
Dky
5.2. Effects of selection of r The second task was to understand the impact of r selection on the final results. The noise level was set at s ¼ 0.55 C and initial values were k0 ¼ [0.1 0.1]. Various r's (3, 6, and 9) were taken for Table 2 Inverse solution of k with different initial values (s ¼ 0 C).
1 2 3
Initial
this numerical test. Table 3 and Fig. 8 show the results of k and qðnÞ, respectively. Again, the relative errors for kx are larger than those of ky because of the sensitivities as mentioned in Section 5.1. Also, kx exhibits a significant error when r ¼ 9, while it may be less affected when r < 6. On the other hand, the future time step has noticeable effects on the results of qðnÞ. As seen from Fig. 8, the case of r ¼ 3 fluctuates vastly at the two ends of the estimated heat flux when s ¼ 0.55 C is included, compared to the first task in Fig. 6 with s ¼ 0. This issue can be mitigated by larger r's (6 and 9); however, the rapid changes at 14 s and 42 s cannot be captured. The fluctuation occurs due to low sensitivity when the moving heat source begins and ends because it is far from the sensors. A larger r behaves like a moving average function to smooth the noisy profile but also creates more bias errors to the result.
(15a)
where Ti ðjÞjkx is the temperature response using kx and Dkx is a small disturbance. The S values are changing over time due to the moving heat source, as the results shown in Fig. 7. Since TC1eTC4 are quite symmetric along the grinding groove, Sx,i is similar at the measurement sensors, and Sy,i is nearly mirrored. The absolute values of Sx,i are significantly lower than those of Sy,i. This indicates the uncertainty of kx, which leads to a greater error when different initial values are given.
#
Fig. 6. Inverse solution of heat flux qðnÞ with different initial k values.
Estimated
Errors
kx
ky
kx
ky
kx
ky
0.1 2.0 10.0
0.1 5.0 20.0
0.052 0.054 0.051
0.030 0.030 0.030
4.0% 8.0% 2.0%
0.0% 0.0% 0.0%
With the FE-model setup as described in Section 5, the experimentally measured temperature data were applied to IHTM to calculate the heat flux distribution with initial k0 ¼ [0.1 0.1]. Two r's, 3 and 9, were used to compare the results. The experimental setup is briefly shown in Fig. 9, which is detailed in Ref. [4]. TC1eTC4 were all K-type thermocouples at designated positions relative to the grinding path. The actual positions were measured under a microscope and used in the IHTM. The grinding bit was tilted by 30 from the sample surface and moved horizontally at a velocity of 20 mm/min. The bone sample was flattened to ensure a constant cutting depth of 0.4 mm. Fig. 10 shows the temperature fitting results at thermocouples TC1eTC4 for both r ¼ 3 and 9, respectively. A good agreement between FEM and actual temperature can be seen in both cases with maximum temperature discrepancy less than 3 C. This indicates an adequate amount of degrees of freedom in the heat flux distribution that can be tuned to match with the temperatures. It is also found that selection of r has no significant effects on the temporal temperature fitting, likely due to the fact that the temporal change of qðnÞ is small in this experiment. Fig. 11 shows the results of qðnÞ for r ¼ 3 and r ¼ 9, respectively. Both cases exhibit high heat flux value in the beginning due to the fact that the sensitivities are low when the moving heat source is far from the sensors, similar to the numerical case in Fig. 8. Therefore, this region (<15 s) may not be taken as a real heat flux change. In comparison between the two cases, data of r ¼ 3 has
90
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
Fig. 7. Sensitivity coefficients for (a) kx and (b) ky.
more fluctuating profile than that of r ¼ 9. This is because a smaller r is more sensitive to the change of heat flux while also easily affected by the noise. However, despite some differences in the profiles, both cases show a convergent heat flux around the level of 0.1 W/mm2 after 15 s. Fig. 12 shows the results of the heat flux spatial distribution calculated by r ¼ 3 and r ¼ 9 at 30 s (about half-way through the process) in comparison to the baseline, the grinding theory-based model in the prior study. The calculated geometrical variables [kx ky] of r ¼ 3 and r ¼ 9 are [0.106 0.073] and [0.009 0.081], respectively. Note that the geometrical variables (k) are constant over time. Considering the fluctuating profile of case r ¼ 3, the total heat generation (using Eq. (4)) ranges from 0.29 to 0.47 W. The total heat generations for r ¼ 9 and baseline are 0.35 W and 0.38 W,
respectively. Although the resultant heat is similar among three, significant variation exists for the case of r ¼ 3, which has been discussed in Section 5.2. It is also important to point out that r ¼ 3 and r ¼ 9 generate different spatial distributions as shown by Fig. 12(a) and (b). This is because of the uncertainty of kx as a result of the current thermocouple arrangements, as described in Section 5.1. The r plays a role like moving average of the temperature profile. Since the heat source moves in the x-direction, a larger r can further reduce the sensitivity of kx, resulting in a flat heat flux distribution. However, it can still be seen a similar ascending trend toward the feed direction among all cases. Apart from kx, ky and qðnÞ are of good agreement between two cases. The baseline model adopts a triangle distribution with a ratio of 0.33 across two elements in y0 direction. For the case r ¼ 3 shown in Fig. 12(a), the ratio ranges from 0.23 to 0.55 (the case of 0 heat flux is not counted). The case of r ¼ 9 with qðnÞ about 0.1 W/mm2 has the heat flux ratio between 0.27 and 0.37 (Fig. 12(b)). In summary, IHTM depends heavily on the selection of input positions since they determine the sensitivities of the geometrical variables. Fig. 10(a) and (b) shows almost identical temperature fitting results regardless the spatial differences between cases r ¼ 3 and r ¼ 9 (Fig. 12) due to the low sensitivity of kx to the input data. The parameter r does not directly determine the heat flux spatial distribution but can influence the searching of kx. One solution
Fig. 8. Inverse solution of heat flux qðnÞ with different future time steps.
Fig. 9. Experimental setup for the temperature measurements during bone grinding. Actual thermocouple positions were measured under a microscope for IHTM.
Table 3 Inverse solution of k with different r (s ¼ 0.55 C). #
4 5 6
r
3 6 9
Initial
Estimated
Errors
kx
ky
kx
ky
kx
ky
0.1 0.1 0.1
0.1 0.1 0.1
0.073 0.064 0.136
0.031 0.032 0.030
46.0% 28.0% 172.0%
3.33% 6.67% 0.00%
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
91
Fig. 10. Temperature fitting results in IHTM using (a) r ¼ 3 and (b) r ¼ 9.
Fig. 11. Results of qðnÞ using (a) r ¼ 3 and (b) r ¼ 9.
Fig. 12. Results of the spatial heat flux using (a) r ¼ 3, (b) r ¼ 9 at 30 s grinding time, and (c) the grinding-theory-based model (prior work).
could be adding thermocouples on the grinding path, lined up in xdirection, to increase the sensitivity of kx. 7. Conclusions In this study, a methodology named IHTM was developed to estimate the temporal and spatial distributions of the bone grinding heat flux. This SFSM-SQP-coupled IHTM could perform simultaneous estimation of spatial and temporal heat flux distributions. Overall, the results have proven that i) a nearly triangular distribution along the traverse direction of the grinding burr, ii) a
higher heat flux magnitude toward the feed direction, and iii) an almost time-invariant heat flux, indicating zero the heat sink effect. Moreover, it was found that the inverse results were influenced by the future time step (r) and the initial values of geometrical variables (k), both of which depended on the experimental setup. Ideally, a low-noise, high-accuracy data is preferred when using smaller r to capture rapid heat flux changes over time. A setup that has high-sensitivity of geometrical variables to the sensor temperatures should be used to avoid the disturbance of initial values. Although this study is limited by the experimental setup in the
92
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
prior work, it is possible to improve the setup based on the IHTM for a better outcome. For examples, thermocouples may be embedded under the grinding path to increase the sensitivity of kx; more thermocouples may be used along the feed direction to increase the sensitivity of qðnÞ at all positions. To further improve the robustness of this method, the future work will focus on the selection of r and sensor arrangements based on the expected heat source profile. These two factors were pre-defined in this study, but could significantly influence the inverse solutions with improper selections. A well-developed IHTM can be used to characterize and compare the heat fluxes under different tool types and operating conditions, such as feed, spindle speed, and tilt angle, for tool optimization and thermal damage prediction. Acknowledgments The authors are grateful to the support of the National Natural Science Foundation of China (Nos. 51176211 and 51505427) and Zhejiang Provincial Natural Science Foundation of China (LQ16E050011). References [1] S. Karmani, The thermal properties of bone and the effects of surgical intervention, Curr. Orthop. 20 (2006) 52e58. [2] T. Udiljak, D. Ciglar, S. Skoric, Investigation into bone drilling and thermal bone necrosis, Adv. Manuf. Eng. Manag. 2 (2007) 103e112. [3] N. McDannold, N. Vykhodtseva, F. Jolesz, K. Hynynen, MRI investigation of the threshold for thermally induced blood-brain barrier disruption and brain tissue damage in the rabbit brain, Magn. Reson. Med. 51 (2004) 913e923. [4] L. Zhang, B.L. Tai, G. Wang, K. Zhang, S. Sullivan, A.J. Shih, Thermal model to investigate the temperature in bone grinding for skull base neurosurgery, Med. Eng. Phys. 35 (2013) 1391e1398.
[5] G. Augustin, T. Zigman, S. Davila, T. Udilljak, T. Staroveski, D. Brezak, S. Babic, Cortical bone drilling and thermal osteonecrosis, Clin. Biomech. 27 (2012) 313e325. [6] M.T. Hillery, I. Shuaib, Temperature effects in the drilling of human and bovine bone, J. Mater. Process. Technol. 92e93 (1999) 302e308. [7] S. Sezek, B. Aksakal, F. Karaca, Influence of drill parameters on bone temperature and necrosis: a FEM modelling and in vitro experiments, Comput. Mater. Sci. 60 (2012) 13e18. [8] S. Kondo, Y. Okada, H. Iseki, T. Hori, K. Takakura, A. Kobayashi, H. Nagata, Thermological study of drilling bone tissue with a high-speed drill, Neurosurgery 46 (2000) 1162e1168. [9] L. Zhang, B.L. Tai, A.C. Wang, A.J. Shih, Mist cooling in neurosurgical bone grinding, CIRP Ann.-Manuf. Technol. 62 (2013) 367e370. [10] A.J. Shih, B.L. Tai, L. Zhang, S. Sullivan, S. Malkin, Prediction of bone grinding temperature in skull base neurosurgery, CIRP Ann.-Manuf. Technol. 61 (2012) 307e310. [11] S. Malkin, C. Guo, Grinding Technology: Theory and Applications of Machining with Abrasives, Industrial Press, New York, 2008. [12] A. Brosse, P. Naisson, H. Hamdi, J.M. Bergheau, J. Mater. Process. Technol. 201 (2008) 590e595. [13] J. Shi, J. Wang, Inverse problem of estimating space and time dependent hot surface heat flux in transient transpiration cooling process, Int. J. Therm. Sci. 48 (2009) 1398e1404. [14] J. Zhou, Y. Zhang, J.K. Chen, Z.C. Feng, Inverse estimation of spatially and temporally varying heating boundary conditions of a two-dimensional object, Int. J. Therm. Sci. 49 (2010) 1669e1679. [15] H.-J. Kim, N.-K. Kim, J.-S. Kwak, Heat flux distribution model by sequential algorithm of inverse heat transfer for determining workpiece temperature in creep feed grinding, Int. J. Mach. Tools Manuf. 46 (2006) 2086e2093. [16] V.M. Luchesi, R.T. Coelho, An inverse method to estimate the moving heat source in machining process, Appl. Therm. Eng. 45e46 (2012) 64e78. [17] J.V. Beck, B. Blackwell, C.R. StClair, Inverse Heat Conduction: Ill-posed Problems, Wiley Press, New York, 1985. [18] F. Samadi, F. Kowsary, A. Sarchami, Estimation of heat flux imposed on the rake face of a cutting tool: a nonlinear, complex geometry inverse heat conduction case study, Int. Commun. Heat Mass Transf. 39 (2012) 298e303. [19] A. Behbahani-nia, F. Kowsary, A dual reciprocity BE-based sequential function specification solution method for inverse heat conduction problems, Int. J. Heat Mass Transf. 47 (2004) 1247e1255. [20] P.T. Boggs, J.W. Tolle, Sequential quadratic programming, Acta Numer. 4 (1996) 1e51.
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
An inverse method to reconstruct the heat flux produced by bone grinding tools Guangjun Wang a, b, Lihui Zhang c, Xudong Wang a, Bruce L. Tai d, * a
School of Power Engineering, Chongqing University, Chongqing 400044, PR China Key Laboratory of Low-grade Energy Utilization Technologies and Systems, Ministry of Education, Chongqing University, Chongqing 400044, PR China c Key Laboratory of Special Purpose Equipment and Advanced Processing Technology, Zhejiang University of Technology, Ministry of Education, Hangzhou, 310014, PR China d Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843, USA b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 20 August 2015 Received in revised form 29 September 2015 Accepted 15 October 2015 Available online xxx
Surgical bone grinding using high-speed, small spherical abrasive tools could cause thermal damage to the surrounding neurovascular structures. Our prior work developed a thermal model based on the twodimensional grinding theory to calculate the temperature distribution. To verify the grinding-theorybased model, this study used experimental data incorporated with an inverse heat transfer method to mathematically estimate the heat flux distribution. This inverse method also considered a time-varying system to reflect a temporal change of the heat flux. Specifically, a coupled approach combining sequential function specification method (SFSM) and sequential quadratic programming (SQP) was employed to calculate the temporal and spatial variables simultaneously. Numerical tests were performed to determine the effectiveness and the limitations of this method. Then, the experimental data of prior work was applied to reconstruct the heat flux. The results verified a nearly time-invariant heat flux and demonstrated a consistent trend in the spatial distribution. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Bone grinding Inverse heat transfer method Finite element method
1. Introduction Bone grinding with high-speed abrasive tools is often used in orthopedic surgery, neurosurgery, otolaryngology, and plastic surgery. Heat is of concern in these procedures since grinding is known to produce more energy compared to other cutting processes, consequently leading to a greater risk to adjacent living tissues. Bone necrosis can occur at about 50 C [1,2]; neural tissues are more vulnerable to the elevated temperature with a threshold level of 43 C [3]. In a minimally invasive skull base surgery, for example, neurosurgeons use a small grinding bit to remove cranial bone surrounded by optic nerves and internal carotid arteries to provide exposure for tumor removal. Excessive heat can cause catastrophic consequences such as stroke and vision damage [4]. Thermal issues in bone cutting processes have been studied extensively in the fields of dentistry, orthopedics and trauma [5e7]. Their results provided the knowledge in tool selection, proper operating parameters, and the level of heat production in the
* Corresponding author. E-mail address: [email protected] (B.L. Tai). http://dx.doi.org/10.1016/j.ijthermalsci.2015.10.021 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.
process. However, there is limited literature investigating bone grinding with abrasive tools. Kondo et al. [8] performed grinding experiments on human cranial bones and recorded temperature with the thermography. Our prior work in collaboration with Shih et al. [4,9,10] conducted grinding experiments on bovine bones and used finite element method (FEM) for estimating temperature distribution. The heat flux distribution was assumed time-invariant in FEM and the spatial distribution was built based on the 2D grinding theory [11]. Despite a decent agreement between FEM and experiments, the heat flux distribution was not completely verified. This study, therefore, is aimed to develop an inverse heat transfer method (IHTM) to solve the heat flux distribution based on experimental inputs and compare to the prior work [4]. In essence, IHTM estimates the unknown thermal boundary conditions based on measurable temperature information [12e16]. The current algorithms for solving IHTM include regularization method, sequential function specification method (SFSM), gradient-based optimization methods and so on. Among them, SFSM is useful in the transient heat transfer problem [17], provided that it can obtain stable inverse solution when a proper number of the future time steps is selected. Samadi et al. [18] applied SFSM to estimate transient heat flux imposed on the rake face of a cutting
86
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
tool during a cutting operation. Behbahani-nia et al. [19] combined dual reciprocity boundary element and SFSM to estimate the time and space-varying boundary heat flux in a two-dimensional configuration. However, when applying SFSM to a moving heat source scenario with unknown magnitude and spatial distribution, the sensitivity coefficients of measured temperature with respect to each elemental heat flux component need to be re-calculated, because they change over time due to the change of sensor position relative to the moving heat source. To overcome this issue, this paper presents a hybrid method using both SFSM and sequential quadratic programming (SQP). The total heat flux distribution is considered to consist of two parts f(x0 , y0 ) and qðtÞ. The former describes the spatial distribution of heat flux; the latter reflects the variation of average heat flux over time. SFSM and SQP are then applied to solve f(x0 , y0 ) and qðtÞ, respectively. This paper is organized as follows. The grinding-theory-based model from our prior work is first reviewed and then followed by the mathematical model used for IHTM. A series of numerical tests are conducted to verify the feasibility and identify limitations of the proposed method, including the effects of initial values and the number of future time steps in SFSM. Finally, this method is applied to the experimental data to estimate the heat flux distribution and compare to the grinding-theory-based model. 2. Review of the grinding theory-based model Fig. 1 shows a schematic view of the grinding process on a bone surface. The small spherical grinding tool with diameter D (typically 2e4 mm) and rotational speed u moves across the bone surface along the feeding direction (x-axis) at the speed v and the cutting depth ac. A tilt angle a is defined as the angle between the tangent to the surface and the axis of the grinding tool. In the grinding theory-based model, the spherical tool is divided into a finite amount of elemental grinding wheels (EGW). Each wheel performs grinding individually with its own tangential speed (re u) and depth of cut (ae). The averaged heat flux under an EGW (qe) then can be calculated by
qe ¼
Kðre uÞae Ae
(1)
Fig. 2 shows the finite element model (FE-model) to which the grinding model is applied. The ground groove is located in the middle of the top surface. The heat flux (shaded region) moves along this groove during the grinding process. TCi represents the ith temperature sensor placed on both sides of the grinding groove. Since this is a pure heat transfer model without considering material removal, the groove is built in the model and simplified as an existing slot with the width of lb. This width is defined by the geometrical intersection between the spherical tool and bone surface at point p, as shown in Fig. 1. The number of EGWs in contact with the surface is determined by the EGW size and la, which is the horizontal distance between points p and q. As the example shown in Fig. 2, there are four active EGWs and each one has two elements underneath. The relative strength of heat flux (i.e., with an undefined K) along the feed direction is calculated based on Eq. (1), and across the feed direction is fitted by a triangular distribution of ratio 1 to 0.33. Finally, the experimentally measured temperatures at sensors are used to adjust the FE-model to find an optimal K in Eq. (1), and complete constructing the heat flux. In this model, the spatial heat flux distribution is determined based on the EGW concept and the assumption of triangular distribution. Also, the scale factor K is assumed a constant throughout the process for a time-invariant system. To further justify these assumptions, a mathematical model, along with a time-varying setting, is used to reconstruct the spatial and temporal distributions, as detailed in the following sections. 3. Mathematical heat flux model Rather than using EGWs to map the spatial heat flux distribution to the FE-model, the mathematical model describes the heat flux by two geometrical variables, kx and ky. Assuming a moving coordinate system x0 ey0, the coordinate conversion can be described by
x0 ¼ x þ vt and y0 ¼ y þ ly 2 þ lb =2
(2)
where v is a constant feed rate and t is time. Then, the heat flux field on the thermal load surface is defined using kx and ky by
qðx0 ; y0 ; tÞ ¼ kx x0 þ ky y0 þ cðtÞ
(3)
where Ae is the area where the heat is applied under the EGW and K is a scale factor. The averaged heat flux qe is distributed in a triangular shape on the workpiece surface as shown in Fig. 1. The shape of triangular distribution depends on the grinding direction [11]; the feed-in edge has a thicker chip load and thus produces more heat.
where c(t) is the magnitude of heat flux at time t and location (x0 , y0 ) ¼ (0,0). The geometrical variables are assumed constant and time invariant in this setting. The temporal term c(t) alters the heat flux magnitude during the grinding process to reflect the environmental noise, inhomogeneous material, or the heat sink effect of the tool. At an instantaneous time, the heat flux field q(x0 , y0, t) can be illustrated by a plane as shown in Fig. 3. Although a high-order
Fig. 1. Schematic view of the bone grinding process with spherical grinding tool.
Fig. 2. Finite element thermal model of bone grinding.
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
87
Fig. 4. Temperature response to a discretized heat flux input. Fig. 3. Schematic view of the spatial distribution of heat flux over rectangular grinding zone.
polynomial or power law may be used to form a more complex shape, they will increase the number of unknown variables to solve in the following inverse method. The total grinding heat Q(t) can be obtained by integrating the heat flux field over the thermal load surface:
Zla Zlb Q ðtÞ ¼ 0
qðx0 ; y0 ; tÞdy0 dx0
(4)
0
By substituting Eq. (3) into Eq. (4), c(t) is obtained as follows:
cðtÞ ¼ qðtÞ 0:5$kx la 0:5$ky lb
(5)
where qðtÞ ¼ Q ðtÞ=la lb is the average value of heat flux in the grinding zone at time t. By substituting Eq. (5) into Eq. (3), Eq. (3) can be modified as follows: 0
0
0
0
qðx ; y ; tÞ ¼ f ðx ; y Þ þ qðtÞ
(6)
where f ðx0 ; y0 Þ ¼ kx x0 þ ky y0 0:5$kx la 0:5$ky lb . Assume k ¼ [kx ky], k and qðtÞ can uniquely determine the heat flux distribution at any given time for a direct heat transfer problem. To map this mathematical model to the FE-model (Fig. 2), the heat flux in the geometric center of each element is taken as an average to be uniformly applied to the corresponding element. The temperature Ti(j) at a given location TCi and time step j can be obtained by
Ti ðjÞ ¼ Ti ðjÞjq¼0 þ Ti ðjÞjf ¼0 ði ¼ 1; 2; …; G;
j ¼ 1; 2; …; MÞ
(7)
where G is the number of temperature sensors and M is the total time steps; Ti ðjÞjq¼0 represents the response temperature Ti(j) when only f(x0 , y0 ) is applied during the time interval [0, j]; Ti ðjÞjf ¼0 represents the response temperature Ti(j) when only qðtÞ is applied during the time interval [0, j]. Note that time step j stands for a total time of t ¼ j Dt, where Dt is the time interval between steps in FEM. The Ti ðjÞjq¼0 is numerically solved by FEM with the boundary condition q(x0 , y0, t) ¼ f(x0 , y0 ) and a pre-defined initial temperature T0 and k. The Ti ðjÞjf ¼0 is described by a linear time-invariant system. The discretized heat flux qðtÞ is expressed by qðnÞ. Each component in qðnÞ represents the average heat flux within a time interval Dt at step n, as shown in Fig. 4. The temperature response can be seen as an accumulated result of all heat flux components qðnÞ at different time steps and positions, such that
Ti ðjÞjf ¼0 ¼
j X n¼1
Tin ðjÞ ¼
j X n¼1
Ii ðn; j nÞqðnÞ
(8a)
where Tin ðjÞ is the temperature response to an individual heat flux component qðnÞ, and Ii(n, j n) is the sensitivity coefficient. Because of a moving heat source, the first term of Ii(n, j n) means the position of the heat source at time step n; the second term is the time interval between the occurrence of heat input and the temperature. In other words, Ii(n, j n) can be seen as a temperature response at time step j, as shown in Fig. 4, when a heat flux q0 ¼ 1 [W/mm2] is applied at time step n. Note that time steps n and j are equivalent but used in different time domains. Ii(n, j n) can be mathematically expressed by Eq. (8b), but it is calculated as a discrete function using FEM by placing q0 at different times and locations.
Ii ðn; j nÞ ¼
vTin ðjÞ vqðnÞ
(8b)
4. Inverse heat transfer method The average heat flux qðnÞ and the geometrical variables k are the unknowns in this inverse problem. An inverse method coupling SFSM and SQP is introduced in this section to solve this problem. 4.1. Estimation of qðnÞ using SFSM SFSM [16] is applied to estimate the time-varying qðnÞ with a given k. At any time step N, the components from qð1Þ to qðN 1Þ are known. To estimate qðNÞ for the step interval [N 1, N], the temperature measurements of different sensors within the time period [N 1, N þ r 1] are used, where r is a pre-selected number of future time steps. For stability, it is temporarily assumed that qðnÞ is constant during the period of [N 1, N þ r 1]. Then, the objective function Jq for estimating qðNÞ is defined as the sum of squared error function,
Jq ðqðNÞÞ ¼
G Nþr1 X X i¼1
Timea ðjÞ TiFEM ðjÞ
2
(9)
j¼N
where Timea and TiFEM are the measured and FEM-calculated temperatures at sensor TCi. For j ¼ N, N þ 1, …, N þ r 1, the future temperature TiFEM ðjÞ is
TiFEM ðjÞ ¼ TiFEM ðjÞjq¼0 þ
N1 X n¼1
Ii ðn; j nÞqðnÞ þ
j X
! Ii ðn; j nÞqðNÞ
n¼N
(10) Then, by substituting Eq. (10) to Eq. (9), differentiating Eq. (9) with respect to qðNÞ, and setting it equal zero, the heat flux at the current time step N is given by:
88
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
PG
PNþr1
i¼1
j¼N
Timea ðjÞ TiFEM ðjÞ
qðNÞ ¼
q¼0
PN1 n¼1
! Ii ðn; j nÞqðnÞ $
Pj
n¼N Ii ðn; j
nÞ
!2
(11)
Since the SFSM method requires no iteration and involvement of FEM, it has high computational efficiency; however, selection of a proper future time step r could influence the result. The r is selected based on the sensitivity of the heat transfer system (i.e., thermal diffusivity and sensor location) and the noise level of measured data. In general, a rapid response system requires a smaller r; a high noise system needs a larger r. The r could also be a function of time if the relative position between sensors and heat source varies. For simplicity, a constant r is used here for analysis.
(m K), and initial temperature 20 C. The total number of discrete elements along grinding groove was 52 2 with the element dimensions of 0.38 mm 1.2 mm. Locations of the temperature sensors in the experiment are listed in Table 1 based on the configuration in Fig. 2. A heat flux qðtÞ with a rapid temporal change was artificially created for testing an extreme case:
PG
i¼1
PNþr1 n¼N
Pj
n¼N Ii ðn; j
nÞ
4.2. Estimation of f(x0 , y0 ) using SQP method The basic concept of SQP is to obtain an approximate solution by solving a sequence of quadratic programming sub-problems. This algorithm has global and superliner convergence [20]. Using k ¼ [kx ky] as independent variables, the objective function F(k) is defined as the sum of squared errors between the measured temperature Timea and calculated temperature TiFEM of all temperature sensors over the whole time domain, such that
FðkÞ ¼
G X M X
Timea ðjÞ TiFEM ðjÞ
i¼1 j¼1
2 k;qðnÞ
qðt Þ ¼
(12)
0:2; 0:3;
t2ð0; M=4Þ∪ð3M=4; M Þ t2½M=4; 3M=4
h
i . W mm2
where M is the total amount of time steps. Taking kx ¼ 0.05 W/mm3 and ky ¼ 0.03 W/mm3, the corresponding spatial distribution function of heat source was given as:
f ðx0 ; y0 Þ ¼ 0:05x0 0:03y0 þ 0:006
h
W=mm2
i
(13b)
According to the heat flux artificially generated in Eqs. (13a) and (13b), this “exact” temperature Tiexa at given sensors was obtained by solving a direct heat transfer problem. Then, Tiexa was used to
where TiFEM is calculated with pre-defined qðnÞ (from Section 4.1) and k using FEM. The commercial software, MATLAB, has fmincon function in its Optimization Toolbox, which can exploit SQP for identifying optimal k to minimize Eq. (12). 4.3. Coupling of SFSM and SQP Coupling of SFSM and SQP is formulated for simultaneous estimation of k and qðnÞ, as shown in the flowchart in Fig. 5. In the coupled algorithm, the optimization of k using SQP is considered to be mainline, and the inverse estimation of qðnÞ is treated as the subroutine. During an iteration, k is fixed through the process while qðnÞ is optimized using SFSM to minimize the objection function, Jq (Eq. (9)). Then, with an optimized qðnÞ, SQP is applied to the objective function F(k) (Eq. (12)) to determine a new set of k for the next iteration. The iteration continues until k or F(k) are converged to the step interval less than d (¼104). These stopping criteria are tuned to be two magnitudes higher than the default settings (106), based on trial and error, considering the noise level and the number of inputs. 5. Numerical tests of IHTM A series of numerical tests were performed to measure the robustness of the inverse method and to identify the effects of initial values of k and future time step r. The following parameters were used to setup the FE-model in accordance with our prior work [4]: D ¼ 4 mm, u ¼ 60,000 rpm, a ¼ 30 , v ¼ 20 mm/min, ac ¼ 0.4 mm, lx ¼ 20 mm, ly ¼ 20 mm, lz ¼ 6 mm, the bone density 2050 kg/m3, specific heat 516 J/(kg K), thermal conductivity 0.55 W/
(13a)
Fig. 5. The flowchart of SFSM and SQP coupling strategy.
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
89
Table 1 Locations of the temperature sensors in experiment (unit: mm). Sensors
TC1
TC2
TC3
TC4
(x, y, z)
(8.5, 6.7, 6.0)
(14.2, 6.9, 6.0)
(8.1, 13.5, 6.0)
(13.5, 13.5, 6.0)
generate the simulated temperature measurements using Eq. (14) with randomized errors:
Timea ¼ Tiexa þ ms
(14)
where m is a random variable over the range 2.58 m 2.58 following a normal distribution with zero mean and s is the standard deviation. The thermocouples in the experiments (36 gauge Ktype, Omega Engineering) had a maximum variation of 1.1 C according to the data sheet, and thus s was set as 0.55 C. With a normal distribution setting, 95% of the randomized errors were within 1.1 C. 5.1. Effects of initial values of k The first task was to study the effects of different initial k values on the final results. The r was set as 3 and s was 0 (no noise) for a controlled condition. The calculated results for k and qðnÞ are shown in Table 2 and Fig. 6, respectively. It is found that both inverse solutions are not sensitive to the initial values of k. Errors of kx range from 2% to 8% and ky remains 0.0%. The results of qðnÞ are fairly consistent. However, qðnÞ can hardly reproduce a rapid profile change due to the use of future time steps for stability, which will be discussed in Section 5.2. Compared to ky, a relatively large error of kx is due to its smaller sensitivity. This sensitivity coefficient (S) can be defined by the temperature difference caused by a small change in geometrical variables, such that
Sx;i ðjÞ ¼
Sy;i ðjÞ ¼
Ti ðjÞjkx Ti ðjÞjkx Dkx Dkx
6. Heat flux reconstruction based on experimental data
Ti ðjÞjky Ti ðjÞjky Dky
(15b)
Dky
5.2. Effects of selection of r The second task was to understand the impact of r selection on the final results. The noise level was set at s ¼ 0.55 C and initial values were k0 ¼ [0.1 0.1]. Various r's (3, 6, and 9) were taken for Table 2 Inverse solution of k with different initial values (s ¼ 0 C).
1 2 3
Initial
this numerical test. Table 3 and Fig. 8 show the results of k and qðnÞ, respectively. Again, the relative errors for kx are larger than those of ky because of the sensitivities as mentioned in Section 5.1. Also, kx exhibits a significant error when r ¼ 9, while it may be less affected when r < 6. On the other hand, the future time step has noticeable effects on the results of qðnÞ. As seen from Fig. 8, the case of r ¼ 3 fluctuates vastly at the two ends of the estimated heat flux when s ¼ 0.55 C is included, compared to the first task in Fig. 6 with s ¼ 0. This issue can be mitigated by larger r's (6 and 9); however, the rapid changes at 14 s and 42 s cannot be captured. The fluctuation occurs due to low sensitivity when the moving heat source begins and ends because it is far from the sensors. A larger r behaves like a moving average function to smooth the noisy profile but also creates more bias errors to the result.
(15a)
where Ti ðjÞjkx is the temperature response using kx and Dkx is a small disturbance. The S values are changing over time due to the moving heat source, as the results shown in Fig. 7. Since TC1eTC4 are quite symmetric along the grinding groove, Sx,i is similar at the measurement sensors, and Sy,i is nearly mirrored. The absolute values of Sx,i are significantly lower than those of Sy,i. This indicates the uncertainty of kx, which leads to a greater error when different initial values are given.
#
Fig. 6. Inverse solution of heat flux qðnÞ with different initial k values.
Estimated
Errors
kx
ky
kx
ky
kx
ky
0.1 2.0 10.0
0.1 5.0 20.0
0.052 0.054 0.051
0.030 0.030 0.030
4.0% 8.0% 2.0%
0.0% 0.0% 0.0%
With the FE-model setup as described in Section 5, the experimentally measured temperature data were applied to IHTM to calculate the heat flux distribution with initial k0 ¼ [0.1 0.1]. Two r's, 3 and 9, were used to compare the results. The experimental setup is briefly shown in Fig. 9, which is detailed in Ref. [4]. TC1eTC4 were all K-type thermocouples at designated positions relative to the grinding path. The actual positions were measured under a microscope and used in the IHTM. The grinding bit was tilted by 30 from the sample surface and moved horizontally at a velocity of 20 mm/min. The bone sample was flattened to ensure a constant cutting depth of 0.4 mm. Fig. 10 shows the temperature fitting results at thermocouples TC1eTC4 for both r ¼ 3 and 9, respectively. A good agreement between FEM and actual temperature can be seen in both cases with maximum temperature discrepancy less than 3 C. This indicates an adequate amount of degrees of freedom in the heat flux distribution that can be tuned to match with the temperatures. It is also found that selection of r has no significant effects on the temporal temperature fitting, likely due to the fact that the temporal change of qðnÞ is small in this experiment. Fig. 11 shows the results of qðnÞ for r ¼ 3 and r ¼ 9, respectively. Both cases exhibit high heat flux value in the beginning due to the fact that the sensitivities are low when the moving heat source is far from the sensors, similar to the numerical case in Fig. 8. Therefore, this region (<15 s) may not be taken as a real heat flux change. In comparison between the two cases, data of r ¼ 3 has
90
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
Fig. 7. Sensitivity coefficients for (a) kx and (b) ky.
more fluctuating profile than that of r ¼ 9. This is because a smaller r is more sensitive to the change of heat flux while also easily affected by the noise. However, despite some differences in the profiles, both cases show a convergent heat flux around the level of 0.1 W/mm2 after 15 s. Fig. 12 shows the results of the heat flux spatial distribution calculated by r ¼ 3 and r ¼ 9 at 30 s (about half-way through the process) in comparison to the baseline, the grinding theory-based model in the prior study. The calculated geometrical variables [kx ky] of r ¼ 3 and r ¼ 9 are [0.106 0.073] and [0.009 0.081], respectively. Note that the geometrical variables (k) are constant over time. Considering the fluctuating profile of case r ¼ 3, the total heat generation (using Eq. (4)) ranges from 0.29 to 0.47 W. The total heat generations for r ¼ 9 and baseline are 0.35 W and 0.38 W,
respectively. Although the resultant heat is similar among three, significant variation exists for the case of r ¼ 3, which has been discussed in Section 5.2. It is also important to point out that r ¼ 3 and r ¼ 9 generate different spatial distributions as shown by Fig. 12(a) and (b). This is because of the uncertainty of kx as a result of the current thermocouple arrangements, as described in Section 5.1. The r plays a role like moving average of the temperature profile. Since the heat source moves in the x-direction, a larger r can further reduce the sensitivity of kx, resulting in a flat heat flux distribution. However, it can still be seen a similar ascending trend toward the feed direction among all cases. Apart from kx, ky and qðnÞ are of good agreement between two cases. The baseline model adopts a triangle distribution with a ratio of 0.33 across two elements in y0 direction. For the case r ¼ 3 shown in Fig. 12(a), the ratio ranges from 0.23 to 0.55 (the case of 0 heat flux is not counted). The case of r ¼ 9 with qðnÞ about 0.1 W/mm2 has the heat flux ratio between 0.27 and 0.37 (Fig. 12(b)). In summary, IHTM depends heavily on the selection of input positions since they determine the sensitivities of the geometrical variables. Fig. 10(a) and (b) shows almost identical temperature fitting results regardless the spatial differences between cases r ¼ 3 and r ¼ 9 (Fig. 12) due to the low sensitivity of kx to the input data. The parameter r does not directly determine the heat flux spatial distribution but can influence the searching of kx. One solution
Fig. 8. Inverse solution of heat flux qðnÞ with different future time steps.
Fig. 9. Experimental setup for the temperature measurements during bone grinding. Actual thermocouple positions were measured under a microscope for IHTM.
Table 3 Inverse solution of k with different r (s ¼ 0.55 C). #
4 5 6
r
3 6 9
Initial
Estimated
Errors
kx
ky
kx
ky
kx
ky
0.1 0.1 0.1
0.1 0.1 0.1
0.073 0.064 0.136
0.031 0.032 0.030
46.0% 28.0% 172.0%
3.33% 6.67% 0.00%
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
91
Fig. 10. Temperature fitting results in IHTM using (a) r ¼ 3 and (b) r ¼ 9.
Fig. 11. Results of qðnÞ using (a) r ¼ 3 and (b) r ¼ 9.
Fig. 12. Results of the spatial heat flux using (a) r ¼ 3, (b) r ¼ 9 at 30 s grinding time, and (c) the grinding-theory-based model (prior work).
could be adding thermocouples on the grinding path, lined up in xdirection, to increase the sensitivity of kx. 7. Conclusions In this study, a methodology named IHTM was developed to estimate the temporal and spatial distributions of the bone grinding heat flux. This SFSM-SQP-coupled IHTM could perform simultaneous estimation of spatial and temporal heat flux distributions. Overall, the results have proven that i) a nearly triangular distribution along the traverse direction of the grinding burr, ii) a
higher heat flux magnitude toward the feed direction, and iii) an almost time-invariant heat flux, indicating zero the heat sink effect. Moreover, it was found that the inverse results were influenced by the future time step (r) and the initial values of geometrical variables (k), both of which depended on the experimental setup. Ideally, a low-noise, high-accuracy data is preferred when using smaller r to capture rapid heat flux changes over time. A setup that has high-sensitivity of geometrical variables to the sensor temperatures should be used to avoid the disturbance of initial values. Although this study is limited by the experimental setup in the
92
G. Wang et al. / International Journal of Thermal Sciences 101 (2016) 85e92
prior work, it is possible to improve the setup based on the IHTM for a better outcome. For examples, thermocouples may be embedded under the grinding path to increase the sensitivity of kx; more thermocouples may be used along the feed direction to increase the sensitivity of qðnÞ at all positions. To further improve the robustness of this method, the future work will focus on the selection of r and sensor arrangements based on the expected heat source profile. These two factors were pre-defined in this study, but could significantly influence the inverse solutions with improper selections. A well-developed IHTM can be used to characterize and compare the heat fluxes under different tool types and operating conditions, such as feed, spindle speed, and tilt angle, for tool optimization and thermal damage prediction. Acknowledgments The authors are grateful to the support of the National Natural Science Foundation of China (Nos. 51176211 and 51505427) and Zhejiang Provincial Natural Science Foundation of China (LQ16E050011). References [1] S. Karmani, The thermal properties of bone and the effects of surgical intervention, Curr. Orthop. 20 (2006) 52e58. [2] T. Udiljak, D. Ciglar, S. Skoric, Investigation into bone drilling and thermal bone necrosis, Adv. Manuf. Eng. Manag. 2 (2007) 103e112. [3] N. McDannold, N. Vykhodtseva, F. Jolesz, K. Hynynen, MRI investigation of the threshold for thermally induced blood-brain barrier disruption and brain tissue damage in the rabbit brain, Magn. Reson. Med. 51 (2004) 913e923. [4] L. Zhang, B.L. Tai, G. Wang, K. Zhang, S. Sullivan, A.J. Shih, Thermal model to investigate the temperature in bone grinding for skull base neurosurgery, Med. Eng. Phys. 35 (2013) 1391e1398.
[5] G. Augustin, T. Zigman, S. Davila, T. Udilljak, T. Staroveski, D. Brezak, S. Babic, Cortical bone drilling and thermal osteonecrosis, Clin. Biomech. 27 (2012) 313e325. [6] M.T. Hillery, I. Shuaib, Temperature effects in the drilling of human and bovine bone, J. Mater. Process. Technol. 92e93 (1999) 302e308. [7] S. Sezek, B. Aksakal, F. Karaca, Influence of drill parameters on bone temperature and necrosis: a FEM modelling and in vitro experiments, Comput. Mater. Sci. 60 (2012) 13e18. [8] S. Kondo, Y. Okada, H. Iseki, T. Hori, K. Takakura, A. Kobayashi, H. Nagata, Thermological study of drilling bone tissue with a high-speed drill, Neurosurgery 46 (2000) 1162e1168. [9] L. Zhang, B.L. Tai, A.C. Wang, A.J. Shih, Mist cooling in neurosurgical bone grinding, CIRP Ann.-Manuf. Technol. 62 (2013) 367e370. [10] A.J. Shih, B.L. Tai, L. Zhang, S. Sullivan, S. Malkin, Prediction of bone grinding temperature in skull base neurosurgery, CIRP Ann.-Manuf. Technol. 61 (2012) 307e310. [11] S. Malkin, C. Guo, Grinding Technology: Theory and Applications of Machining with Abrasives, Industrial Press, New York, 2008. [12] A. Brosse, P. Naisson, H. Hamdi, J.M. Bergheau, J. Mater. Process. Technol. 201 (2008) 590e595. [13] J. Shi, J. Wang, Inverse problem of estimating space and time dependent hot surface heat flux in transient transpiration cooling process, Int. J. Therm. Sci. 48 (2009) 1398e1404. [14] J. Zhou, Y. Zhang, J.K. Chen, Z.C. Feng, Inverse estimation of spatially and temporally varying heating boundary conditions of a two-dimensional object, Int. J. Therm. Sci. 49 (2010) 1669e1679. [15] H.-J. Kim, N.-K. Kim, J.-S. Kwak, Heat flux distribution model by sequential algorithm of inverse heat transfer for determining workpiece temperature in creep feed grinding, Int. J. Mach. Tools Manuf. 46 (2006) 2086e2093. [16] V.M. Luchesi, R.T. Coelho, An inverse method to estimate the moving heat source in machining process, Appl. Therm. Eng. 45e46 (2012) 64e78. [17] J.V. Beck, B. Blackwell, C.R. StClair, Inverse Heat Conduction: Ill-posed Problems, Wiley Press, New York, 1985. [18] F. Samadi, F. Kowsary, A. Sarchami, Estimation of heat flux imposed on the rake face of a cutting tool: a nonlinear, complex geometry inverse heat conduction case study, Int. Commun. Heat Mass Transf. 39 (2012) 298e303. [19] A. Behbahani-nia, F. Kowsary, A dual reciprocity BE-based sequential function specification solution method for inverse heat conduction problems, Int. J. Heat Mass Transf. 47 (2004) 1247e1255. [20] P.T. Boggs, J.W. Tolle, Sequential quadratic programming, Acta Numer. 4 (1996) 1e51.