高速鉆骨過程中切削力的仿真與實驗研究

張立斌，王宇，趙宏強，趙世佳，胡穎

高速鉆骨過程中切削力的仿真與實驗研究

張立斌1, 2，王宇2, 3，趙宏強1，趙世佳2，胡穎2

(1. 中南大學 機電工程學院，湖南 長沙，410083；2. 中國科學院 深圳先進技術(shù)研究院，廣東 深圳，518055；3. 哈爾濱工業(yè)大學(深圳) 機電工程與自動化學院，廣東 深圳，518055)

基于彈塑性有限元模型，實現(xiàn)高速鉆骨時骨組織?器械交互狀態(tài)的仿真；根據(jù)動物標本鉆骨實驗結(jié)果，通過優(yōu)化方法實現(xiàn)模型損傷參數(shù)的快速辨識；結(jié)合有限元仿真和標本骨實驗，研究高速鉆骨過程中切削力隨不同鉆削參數(shù)的變化規(guī)律。研究結(jié)果表明：鉆削軸向力仿真結(jié)果與實驗結(jié)果較吻合，驗證了損傷參數(shù)的有效性及有限元模型的適用性；手術(shù)中應該選擇較小的進給量和更高的鉆速以實現(xiàn)安全鉆削。

有限元方法；高速鉆骨；損傷參數(shù)；優(yōu)化方法

在骨科手術(shù)中，一些骨科疾病如骨折經(jīng)常會用鋼釘?shù)戎踩塍w固定骨頭，釘?shù)赖馁|(zhì)量會影響植入體的固定效果及術(shù)后恢復，因此鉆釘?shù)朗顷P(guān)鍵的手術(shù)流程。在鉆釘?shù)肋^程中產(chǎn)生的過大鉆削力可能會對骨周圍組織造成一定損傷甚至引起裂紋[1]。為了優(yōu)化鉆削效果，有必要對鉆骨過程中的鉆頭與骨組織的交互作用進行研究。作為實驗方法的重要補充，計算模型特別是有限元模型已被廣泛應用于組織?器械交互的研究[2?6]。LUGHMANI等[1]建立了一個精確的有限元模型研究了低速條件下改變鉆削參數(shù)而引起的鉆削皮質(zhì)骨的軸向力和扭矩變化規(guī)律。SEZEK等[7]基于有限元方法研究了不同的轉(zhuǎn)速、進給量、鉆頭尺寸等參數(shù)和不同性別的骨在鉆削過程中所發(fā)生的溫度變化和熱壞死現(xiàn)象。ALAM等[8?9]基于有限元探究了在平面切削皮質(zhì)骨時切削深度、切削速度等參數(shù)對切削力和切削溫度的影響。JIN等[10]基于力信號提出了在鉆削過程中骨科手術(shù)機器人系統(tǒng)狀態(tài)感知的方法。張青云[11]通過有限元法和實驗進行鉆削皮質(zhì)骨性能的研究。國內(nèi)一些研究人員基于實驗分析了在不同鉆削參數(shù)下皮質(zhì)骨鉆削溫度和鉆削力的變化規(guī)律，如宋金榜等[12]研究了高速微細骨鉆過程中的切削力和切削溫度的變化規(guī)律。上述基于有限元模型模擬鉆骨過程的研究中，采用的鉆削轉(zhuǎn)速與實際臨床醫(yī)用鉆速相差過大，不具臨床參考價值。鉆削仿真中涉及材料的損傷，而其中損傷演化的位移參數(shù)設(shè)置尤為重要。本文作者基于優(yōu)化方法，通過實驗數(shù)據(jù)驗證，實現(xiàn)有限元在模擬高速鉆骨時模型的關(guān)鍵參數(shù)——損傷演化的等效塑性位移的快速辨識，進而完善皮質(zhì)骨鉆削的有限元模型，用于高速鉆削過程研究；通過鉆骨實驗驗證損傷參數(shù)的有效性，其中鉆削轉(zhuǎn)速設(shè)置為8 000~23 000 r/min，接近臨床中醫(yī)用骨鉆的轉(zhuǎn)速。

1 有限元模型

1.1 幾何模型

本文采用標準醫(yī)療鉆頭，其結(jié)構(gòu)參數(shù)如下：直徑為2.5 mm，頂角為118°，螺旋角為25°，長度為2.5 mm。首先，用SolidWorks軟件繪制鉆頭三維幾何模型，然后導入ABAQUS中，進行材料屬性賦予、網(wǎng)格劃分以及邊界條件加載等操作。

研究使用微細的標準醫(yī)療鉆頭，在鉆骨過程中鉆頭對骨頭的作用區(qū)域僅是很小的一部分，因此，在皮質(zhì)骨的建模中將骨組織簡化為圓柱體，直徑為6 mm，高為2 mm。

1.2 材料屬性

1.2.1 骨組織力學特性

骨組織的力學行為在彈性范圍內(nèi)遵循Hill的各向異性材料潛在理論和與速率相關(guān)的彈性準則[1]。由于皮質(zhì)骨在力學上具有較低的各向異性相關(guān)性，為了便于分析，本文把皮質(zhì)骨定義為各向同性的彈塑性材料。鉆頭與皮質(zhì)骨材料屬性如表1所示[1, 13]。

表1 材料參數(shù)[1, 13]

1.2.2 骨組織本構(gòu)模型

在鉆削加工中，力的變化是一個非線性過程，材料處在高溫、大應變和大應變率情況下。Johnson-Cook模型很適合描述這種依賴應變?應變率的材料屬 性[14]。Johnson-Cook材料模型本構(gòu)方程如下：

表2 Johnson-Cook模型參數(shù)

1.2.3 材料失效準則

1.3 網(wǎng)格劃分

鉆頭部分采用自由網(wǎng)格劃分技術(shù)，劃分四面體單元，單元類型為C3D4。骨模型采用細分網(wǎng)格技術(shù)，這對模擬大變形和非線性材料行為有很重要的意義。首先分割模型，在鉆孔附近細化網(wǎng)格，其他部位使用粗網(wǎng)格，這樣劃分更加科學有效還可以提高計算效率，節(jié)省計算資源。減縮積分單元能更好地模擬網(wǎng)格扭曲十分嚴重的問題，因此，骨模型單元類型為C3D8R。減縮積分單元容易變形扭曲，產(chǎn)生零能模式，這里使用增強型沙漏控制。本文骨組織有限元模型共80 124個單元，鉆頭模型共33 615個單元，網(wǎng)格模型如圖1所示。

圖1 網(wǎng)格模型

1.4 邊界條件和接觸

在有限元仿真中，皮質(zhì)骨模型側(cè)面被固定，以限制6個自由度。鉆頭被約束為剛體，同時設(shè)置參考點，并對參考點施加轉(zhuǎn)速和進給量以控制鉆頭的運動。

接觸使用通用接觸運算法則，該法則產(chǎn)生的接觸力基于罰摩擦公式。鉆頭和皮質(zhì)骨接觸時的摩擦模型用庫侖摩擦來描述，該模型用摩擦因數(shù)來表征2個表面之間的摩擦行為，摩擦因數(shù)為0.3[15]。

2 鉆骨實驗

2.1 實驗方法

為了驗證模型的有效性，探究高速鉆骨時軸向力的變化規(guī)律，進行不同鉆削參數(shù)下的鉆骨實驗。鉆骨實驗在三軸直線機器人上進行，實驗裝置如圖2所示。骨鉆的動力由Maxon電機提供，鉆頭通過骨鉆帶動。實驗材料選用新鮮的豬肩胛骨，鉆頭為直徑2.5 mm的標準醫(yī)療不銹鋼鉆，使用ATI(Nano43)力傳感器采集數(shù)據(jù)，傳感器通過采集卡連接電腦。

圖2 鉆骨實驗裝置

實驗開始時，鉆頭和骨皮質(zhì)處于接觸狀態(tài)，待鉆頭轉(zhuǎn)速達到平穩(wěn)狀態(tài)，三軸直線機器人帶動骨鉆以預先設(shè)置的進給量穿透皮質(zhì)骨，同時，ATI力傳感器實時記錄并保存鉆削軸向力數(shù)據(jù)。

鉆骨實驗中采用的鉆削參數(shù)如表3所示。鉆頭與皮質(zhì)骨接觸時的力有一定波動，取峰值時的力在MATLAB軟件中進行均值方差計算，將其計算結(jié)果作為實驗結(jié)果。在有限元仿真中，計算后的力波動較大，使用MATLAB對數(shù)據(jù)濾波后取均值。

表3 鉆削參數(shù)

2.2 實驗結(jié)果

根據(jù)上述方法進行實驗，實驗結(jié)果分別如表4和表5所示。

表4 不同進給量下的軸向力(N=8 000 r/min)

表5 不同轉(zhuǎn)速下的軸向力(f=1 mm/s)

3 損傷演化參數(shù)辨識

3.1 梯度下降法

梯度下降法是一種最早也最簡單的優(yōu)化方法，其迭代方式簡單，易于理解，因此，本文采用梯度下降法對損傷參數(shù)進行優(yōu)化。梯度下降法是沿負梯度方向搜索，求解表達式最小值(或沿上升方向搜索求最大值)的一種方法，過程中越接近目標值，步長越小。

鉆削時材料在損傷演化階段，函數(shù)關(guān)系如下[1]：

3.2 損傷演化的位移參數(shù)辨識

ABAQUS通過損傷演化模擬損傷的發(fā)生，損傷演化的類型使用位移模式，在保持模型其他參數(shù)不變的情況下，損傷演化的等效塑性位移決定了軸向力，通過調(diào)整等效塑性位移參數(shù)優(yōu)化仿真結(jié)果。本文參數(shù)調(diào)節(jié)具體流程如圖3所示。圖3中，為仿真與實驗結(jié)果的誤差，其表達式為

式中：FEA為仿真軸向力；EXP為實驗軸向力。

實驗時保持轉(zhuǎn)速為8 000 r/min，進給量為 1.5 mm/s，實驗得到的軸向力為10.3 N。建立有限元模型，設(shè)置初始損傷演化的等效塑性位移參數(shù)并計算仿真結(jié)果。

圖3 損傷演化調(diào)節(jié)流程圖

將此結(jié)果與實驗結(jié)果進行對比，若相對誤差超出10%，則說明此位移參數(shù)不合理，然后以一定步長調(diào)整此參數(shù)，優(yōu)化仿真結(jié)果，最終確定1個合理的位移參數(shù)。

仿真實驗中等效塑性位移調(diào)整了4次，結(jié)果如圖4所示。保持轉(zhuǎn)速為8 000 r/min，進給量為1.5 mm/s，設(shè)置初始位移參數(shù)為0.10 mm，初始步長為0.1 mm，仿真得到的軸向力為7.5 N，與實驗結(jié)果相比誤差為27%；隨后設(shè)置位移參數(shù)為0.20 mm，仿真結(jié)果仍在誤差范圍外，相對誤差為13%；然后設(shè)置位移參數(shù)為0.30 mm，計算后相對誤差為21%，且仿真與實驗結(jié)果之差大于0 N，此時步長減半；位移參數(shù)取0.25 mm，結(jié)果相對誤差為2%，在允許的范圍內(nèi)。因此，最終本文等效塑性位移設(shè)定為0.25 mm。

1—實驗結(jié)果；2—仿真結(jié)果。

4 實驗驗證

4.1 切削應力

鉆骨Mises應力云圖如圖5所示。在仿真研究中材料分離時的最大應力約為150 MPa，這與文獻[6]中切削骨仿真和拉伸實驗斷裂時的應力基本一致。

圖5 鉆骨Mises應力云圖

4.2 不同鉆削參數(shù)下軸向力變化

保持等效塑性位移為0.25 mm不變，改變進給量和鉆速，計算仿真結(jié)果并與實驗結(jié)果對比，結(jié)果分別如圖6和圖7所示。

從圖6可以看出：保持鉆速為8 000 r/min不變，隨著進給量增加，軸向力也逐漸增大；當進給量從 1.0 mm/s增加到2.0 mm/s時，軸向力在8~15 N范圍內(nèi)波動。實驗軸向力從最低到最高增加了7.12 N，仿真軸向力增加了6.64 N，尤其是在實驗過程中，軸向力幾乎呈線性增長。保持進給量為1.0 mm/s不變，鉆速從8 000 r/min起每次增加5 000 r/min，直至 23 000 r/min(見圖7)，隨著鉆速增加，軸向力減小 2.13 N。因此，在手術(shù)過程中，在其他條件不變的情況下，應盡量減小進給量，增加鉆速，從而產(chǎn)生更小的切削力。

由圖6和圖7可知：改變進給量時軸向力變化幅度為6~7 N，波動較大，而改變鉆速時軸向力變化幅度小于3 N，波動相對較小。通過對比可知進給量對軸向力的影響比鉆速的影響大，這也符合鉆削軸向力經(jīng)驗公式所體現(xiàn)的規(guī)律。

1—實驗結(jié)果；2—仿真結(jié)果。

1—實驗結(jié)果；2—仿真結(jié)果。

本文作者使用相關(guān)系數(shù)對數(shù)據(jù)進行量化，相關(guān)系數(shù)表達式為

式中：和為變量；Cov(,)為和的協(xié)方差，Var[]和Var[]分別為和的方差。通過求解可得圖6和圖7中2組數(shù)據(jù)相關(guān)系數(shù)分別為1=0.997 2，2=0.982 0。從2個相關(guān)系數(shù)可以得到在改變進給量和鉆速時，仿真所得到的鉆削軸向力結(jié)果與實驗結(jié)果有很高的吻合度。這證明了本文作者所提的基于優(yōu)化方法，通過實驗驗證，實現(xiàn)位移參數(shù)快速辨識方法的正確性，同時也證明了仿真所設(shè)置的等效塑性位移參數(shù)設(shè)置是合理、有效的，表明本文作者建立的有限元模型適用于模擬高速鉆骨過程。

5 結(jié)論

1)建立了1個彈塑性有限元模型用來模擬高速鉆骨過程，基于優(yōu)化方法，通過實驗驗證，實現(xiàn)了有限元在模擬高速鉆骨時模型的關(guān)鍵參數(shù)即等效塑性位移參數(shù)的快速辨識。

2)在手術(shù)中，為了避免過大的切削力對骨周圍組織的損傷，應該選擇較小的進給量和更高的鉆速；在臨床手術(shù)中，為了不產(chǎn)生過大的切削力，進給量應小于1.0 mm/s，轉(zhuǎn)速應高于23 000 r/min。

[1] LUGHMANI W A, BOUAZZA-MAROUF K, ASHCROFT I. Drilling in cortical bone: a finite element model and experimental investigations[J]. Journal of the Mechanical Behavior of Biomedical Materials, 2015, 42(8): 32?42.

[2] ZHAO Shijia, GU Linxia, FROEMMING S R. On the importance of modeling stent procedure for predicting arterial mechanics[J]. Journal of Biomechanical Engineering, 2012, 134(12): 121005-1?121005-6.

[3] 聶涌, 馬俊, 胡欽勝, 等. 髖臼假體放置角度對髖臼周圍應力分布的影響[J]. 醫(yī)用生物力學, 2014, 29(4): 1?7..NIE Yong, MA Jun, HU Qinsheng, et al. Influence from acetabular component orientation on stress distribution of periacetabulum[J]. Journal of Medical Biomechanics, 2014, 29(4): 1?7.

[4] Karaca F, Aksakal B, Kom M. Influence of orthopaedic drilling parameters on temperature and histopathology of bovine tibia: an in vitro study[J]. Medical Engineering and Physics, 2011, 33(10): 1221?1227.

[5] Wang Yu, Cao Meng, Zhao Xiangrui, et al. Experimental investigations and finite element simulation of cutting heat in vibrational and conventional drilling of cortical bone[J]. Medical Engineering and Physics, 2014, 36(11): 1408?1415.

[6] Qi Lin, Wang Xiaona, Meng M Q. 3D finite element modeling and analysis of dynamic force in bone drilling for orthopedic surgery[J]. International Journal for Numerical Methods in Biomedical Engineering, 2014, 30(9): 845?856.

[7] SEZEK S, AKSAKAL B, KARACA F. Influence of drill parameters on bone temperature and necrosis: a FEM modelling and in vitro experiments[J]. Computational Materials Science, 2012, 60(5): 13?18.

[8] ALAM K, MITROFANOV A V, SILBERSCHMIDT V V. Finite element analysis of forces of plane cutting of cortical bone[J]. Computational Materials Science, 2009, 46(3): 738?743.

[9] ALAM K, MITROFANOV A V, SILBERSCHMIDT V V. Thermal analysis of orthogonal cutting of cortical bone using finite element simulations[J]. International Journal of Experimental and Computational Biomechanics, 2010, 1(3): 236?251.

[10] JIN Haiyang, HU Ying, DENG Zhen, et al. Model-based state recognition of bone drilling with robotic orthopedic surgery system[C]// The 2014 International Conference on Robotics and Automation (ICRA). Chicago, USA: IEEE,2014: 3538?3543.

[11] 張青云. 基于ABAQUS的皮質(zhì)骨鉆削性能的仿真和試驗研究[D]. 天津: 天津理工大學機電工程學院, 2014: 46.ZHANG Qingyun. Finite element analysis and experimental research of cortical bone drilling performation based on ABAQUS[D]. Tianjin: Tianjin University of Technology. School of Mechanical Engineering, 2014: 46.

[12] 宋金榜, 陳明, 趙梓涵, 等. 高速微細骨鉆過程中切削力與切削溫度試驗研究[J]. 機械設(shè)計與制造, 2015(4): 137?139. 143. SONG Jinbang, CHEN Ming, ZHAO Zihan, et al. An investigation on cutting force and temperature during high speed micro-drilling on bone[J]. Machinery Design & Manufacture, 2015(4): 137?139, 143.

[13] TU Yuankun, CHEN Liwen, CIOU J S, et al. Finite element simulations of bone temperature rise during bone drilling based on a bone analog[J]. Journal of Medical & Biological Engineering, 2013, 33(3): 269?274.

[14] 劉森. 骨銑削過程中切削力和溫度的仿真與實驗研究[D]. 哈爾濱: 哈爾濱工業(yè)大學機電工程與自動化學院, 2014: 58.LIU Sen. Research on simulation and experimentation of cutting force and temperature for bone milling[D]. Harbin: Harbin Institute of Technolog.School of Mechanical Engineering and Autmation, 2014: 58.

[15] MELLAL A, WISKOTT H W A, BOTSIS J, et al. Stimulating effect of implant loading on surrounding bone[J]. Clinical Oral Implants Research, 2004, 15(2): 239?248.

[16] 陳海波, 呂咸青, 喬彥松. 梯度下降法在沉積物粒度分布擬合中的應用[J]. 數(shù)學的實踐與認識, 2011, 41(13): 78?87. CHEN Haibo, Lü Xianqing, QIAO Yansong. Application of gradient descent method to the sedimentary grain-size distribution fitting[J]. Mathematics in Practice and Theory, 2011, 41(13): 78?87.

(編輯 伍錦花)

Computational simulation and experimental research on cutting force during high-speed bone drilling

ZHANG Libin1, 2, WANG Yu2, 3, ZHAO Hongqiang1, ZHAO Shijia2, HU Ying2

(1. College of Mechanical and Electrical Engineering, Central South University, Changsha 410083, China;2. Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China;3. School of Mechanical Engineering and Automation, Harbin Institute of Technology, Shenzhen 518055, China)

A finite element model was constructed to simulate the interaction between bone tissue and device during high-speed bone drilling. Based on the result fromdrilling experiment and the optimization method, the damage parameter used in the finite element model was identified. Combining the finite element simulation and theexperiment, variations of cutting force along with the changes in drilling parameters during the process of high-speed bone drilling were studied. The results show that simulation result of drilling axial force coincides with that ofexperiment, which validates the effectiveness of damage parameter and the applicability of finite element model; lower feeding rate and higher drilling speed should be adopted to ensure drilling safety in orthopedic surgery.

finite element method; high-speed bone drilling; damage parameter; optimization method

10.11817/j.issn.1672-7207.2018.09.011

TH781

A

1672?7207(2018)09?2191?06

2017?10?31；

2017?12?19

國家自然科學基金資助項目(61573336)；深圳市基礎(chǔ)研究重點項目(JCYJ20160608153218487) (Project(61573336, U1613224) supported by National Natural Science Foundation of China; Project(JCYJ20160608153218487) supported by Key Fundamental Research Program of Shenzhen)

胡穎，博士，研究員，從事醫(yī)療機器人研究；E-mail: ying.hu@siat.ac.cn