唐勇濤，魏思奇，程剛，趙建華，蔡繁繁，歐陽漢斌

?

基于兩種不同建模方式的脛骨有限元對比分析

唐勇濤，魏思奇，程剛，趙建華，蔡繁繁，歐陽漢斌

【摘要】目的探討兩種不同建模方式對脛骨有限元模型生物力學(xué)特性的影響。方法采集1例正常志愿者脛骨CT斷層掃描數(shù)據(jù)，對其脛骨模型進(jìn)行三維重建，分別重建出單純脛骨外輪廓模型以及帶松質(zhì)骨和髓腔結(jié)構(gòu)的腔體模型，采用相應(yīng)的灰度公式賦值法和均一賦值法完成模型材料參數(shù)賦值，同時對兩組設(shè)定相同的載荷和邊界條件，模擬站立位狀態(tài)下脛骨上段受到的軸向壓縮和扭轉(zhuǎn)載荷，對比兩組模型的脛骨應(yīng)力和位移分布云圖、軸向壓縮剛度以及扭轉(zhuǎn)剛度。結(jié)果兩組不同方式構(gòu)建的脛骨模型總體位移分布基本一致；灰度賦值組的脛骨有限元模型在兩種載荷作用下產(chǎn)生的應(yīng)力峰值和位移幅值均小于均一賦值組，其軸向壓縮剛度和扭轉(zhuǎn)剛度均高于均一賦值組；均一賦值組模型在脛骨中下段位置的應(yīng)力分布相較于灰度賦值組更為集中，相對符合臨床實際情況。結(jié)論兩種建模方式均可滿足脛骨有限元分析需求，均一賦值法建模具有更高的準(zhǔn)確性，模型應(yīng)用的靈活性更高。

【關(guān)鍵詞】脛骨；骨重建；模型，解剖學(xué)；有限元分析；材料賦值；應(yīng)力，物理；生物力學(xué)；計算機(jī)模擬

作者單位：518106廣東，深圳市光明新區(qū)人民醫(yī)院骨一科（唐勇濤，魏思奇，程剛，趙建華，蔡繁繁）；510515廣州，南方醫(yī)科大學(xué)人體解剖學(xué)教研室（歐陽漢斌）

E-mail：1241195475@qq.com

近年來，以生物力學(xué)為基礎(chǔ)的有限元分析方法在人體骨骼結(jié)構(gòu)力學(xué)分析中發(fā)揮著巨大作用。相對于傳統(tǒng)生物力學(xué)實驗手段，有限元仿真方法能夠更為有效地預(yù)測骨骼在生理性和異常沖擊受力狀態(tài)下的應(yīng)力分布以及多個方向的應(yīng)變情況[1]，預(yù)測準(zhǔn)確性高[2-3]。隨著螺旋CT的廣泛應(yīng)用和醫(yī)學(xué)影像三維重建技術(shù)的日趨成熟，有限元分析方法也在不斷完善，學(xué)者們提出各種類型的建模方法，以獲得更為準(zhǔn)確的預(yù)測模型[4-5]。脛骨作為人體下肢骨的一個重要承載結(jié)構(gòu)，有關(guān)其尸體標(biāo)本結(jié)構(gòu)特性的生物力學(xué)研究已有文獻(xiàn)報道，但運(yùn)用有限元分析方法對脛骨進(jìn)行生物力學(xué)研究，各家報道不一，特別是在脛骨模型材料屬性賦值方法的選用上，尚未得出一致的結(jié)論[6-8]。本文就脛骨有限元分析建模的材料屬性賦值問題展開初步研究，旨在闡明不同建模方法對模型最終分析結(jié)果的影響，為脛骨生物力學(xué)研究提供可靠的參考依據(jù)。

1　材料與方法

1.1脛骨模型數(shù)據(jù)采集與三維重建

選取一名正常志愿者右側(cè)脛骨作為數(shù)據(jù)采集來源，排除受檢者受試肢體相關(guān)骨骼結(jié)構(gòu)先天異常、損傷及退變等病理變化。采用64排螺旋CT機(jī)（西門子公司，德國）沿脛骨干長軸方向進(jìn)行掃描，范圍覆蓋膝關(guān)節(jié)及踝關(guān)節(jié)，掃描層距1 mm，共獲得連續(xù)橫斷面1 000張512×512像素的CT圖像，并將獲得的CT圖像以通用DICOM標(biāo)準(zhǔn)格式刻錄至光盤導(dǎo)出影像工作站。進(jìn)一步將脛骨影像數(shù)據(jù)導(dǎo)入Mimics 14.0軟件（Materialise公司，比利時）進(jìn)行脛骨皮質(zhì)外輪廓、脛骨髓腔、脛骨遠(yuǎn)近端的圖像分割處理，將分割的各個結(jié)構(gòu)蒙板提取處理并計算出脛骨皮質(zhì)骨、松質(zhì)骨以及髓腔的三維模型，導(dǎo)出后以STL文件形式保存（圖1）。

1.2模型表面處理與脛骨腔體結(jié)構(gòu)的實體建模

圖1　脛骨髓腔內(nèi)部復(fù)合結(jié)構(gòu)的二維圖像分割及三維重建

首先將上述模型（皮質(zhì)骨、近端松質(zhì)骨、遠(yuǎn)端松質(zhì)骨、髓腔）的STL文件分別導(dǎo)入Geomagic 12軟件（Geomagic公司，美國）中進(jìn)行模型表面平滑處理和精確曲面構(gòu)建；然后在三維建模軟件UG NX 8.5（西門子公司，德國）中進(jìn)行含髓腔結(jié)構(gòu)的脛骨建模相關(guān)布爾操作，最后得到完整實體模型（圖2），并以STP文件格式導(dǎo)出皮質(zhì)骨和遠(yuǎn)近端松質(zhì)骨部件。

1.3灰度賦值法與均一賦值法脛骨有限元模型建立

將脛骨實體模型導(dǎo)入有限元分析軟件Abaqus 6.10（達(dá)索公司，法國）中劃分四面體網(wǎng)格，并以inp文件形式導(dǎo)入Mimics賦予材料屬性，Mimics將根據(jù)材料屬性公式[9]，在原位二維影像的相應(yīng)位置自動將灰度轉(zhuǎn)換為局部有限元網(wǎng)格單元的彈性模量和泊松比（圖3）；對均一賦值法模型采用均一材料屬性在Abaqus軟件中完成賦值，將兩種骨組織結(jié)構(gòu)均視為各向同性材料，其中皮質(zhì)骨彈性模量為16 700 MPa，松質(zhì)骨彈性模量為155 MPa，泊松比均取值0.3[10]。

1.4軸向壓縮及扭轉(zhuǎn)載荷與邊界條件設(shè)定

圖2　脛骨三維模型的表面處理與實體建模

圖3　灰度賦值法構(gòu)建的脛骨模型材料屬性分布圖

根據(jù)文獻(xiàn)報道，模擬人體站立位狀態(tài)下脛骨平臺關(guān)節(jié)面受到的軸向壓縮載荷為700 N、扭轉(zhuǎn)載荷為20 Nm，載荷作用區(qū)域范圍位于脛骨平臺上方以及脛骨的干骺端上部[11]。另外，將脛骨下端踝關(guān)節(jié)關(guān)節(jié)面上方2.5 cm以下范圍設(shè)置為6個自由度完全約束的邊界條件，然后建立2個分析步以先后施加軸向壓縮和扭轉(zhuǎn)載荷。完成設(shè)置后，利用Abaqus軟件的網(wǎng)格劃分功能對模型進(jìn)行四面體網(wǎng)格劃分，采用C3D4實體單元，獲得兩組網(wǎng)格化模型。其中灰度賦值模型單元數(shù)為217 006個，均一賦值模型單元數(shù)為549 535個。最后將模型提交到有限元求解器中進(jìn)行運(yùn)算。

2　結(jié)果

2.1基于灰度賦值法和均一賦值法的脛骨模型應(yīng)力、位移分布云圖

最終提取出兩組脛骨模型的應(yīng)力分布云圖見圖4。兩種受力狀態(tài)下灰度賦值組脛骨的應(yīng)力幅值要顯著低于均一賦值組：灰度賦值組軸向壓縮的最大應(yīng)力峰值為70.57 MPa，扭矩作用下最大應(yīng)力峰值為17.21 MPa；而均一賦值組軸向壓縮的最大應(yīng)力峰值為108.40 MPa，扭矩作用下最大應(yīng)力峰值為27.42 MPa。此外，灰度賦值組脛骨干整體應(yīng)力分布較為均勻，高應(yīng)力區(qū)主要位于骨干中段以下；均一賦值組脛骨干的整體應(yīng)力分布主要集中在脛骨中下段交界處，呈現(xiàn)出明顯的應(yīng)力集中現(xiàn)象。

兩組脛骨模型的位移分布如圖5所示。兩種受力狀態(tài)下灰度賦值組脛骨發(fā)生的位移幅值略低于均一賦值組：灰度賦值組軸向壓縮最大位移為13.97 mm，扭矩作用下最大位移為3.25 mm；而均一賦值組軸向壓縮最大位移為15.75 mm，扭矩作用下最大位移為3.87 mm。此外，兩組模型無論是在軸向壓縮還是在扭轉(zhuǎn)作用下，脛骨總體產(chǎn)生的位移分布均較一致，總體趨勢符合生理狀態(tài)下脛骨受力后發(fā)生變形的規(guī)律。

2.2兩組脛骨模型軸向壓縮剛度和扭轉(zhuǎn)剛度對比結(jié)果

通過計算載荷-位移曲線的線性斜率，得到兩組模型的整體剛度，灰度賦值組脛骨模型的軸向壓縮剛度和扭轉(zhuǎn)剛度（50.11 N/mm，13.36 Nm/°）均高于均一賦值組（44.44 N/mm，11.23 Nm/°）。

圖4　兩組脛骨模型在軸向壓縮與扭轉(zhuǎn)作用下的應(yīng)力分布云圖（單位：MPa）

3　討論

一直以來，有限元分析方法被廣泛應(yīng)用于航天、土木和機(jī)械結(jié)構(gòu)等工程領(lǐng)域，被公認(rèn)為是一種有效、準(zhǔn)確、低成本的力學(xué)結(jié)構(gòu)分析方法[12-13]。在骨科生物力學(xué)領(lǐng)域，有限元方法可以對人體骨骼、肌肉和韌帶等組織進(jìn)行仿真分析，除了為骨科植入物的結(jié)構(gòu)力學(xué)分析和改良設(shè)計提供可靠依據(jù)外，在人體結(jié)構(gòu)生物力學(xué)基礎(chǔ)研究方面同樣具有廣闊的應(yīng)用前景，有效彌補(bǔ)了傳統(tǒng)生物力學(xué)試驗手段的不足。近年來，國內(nèi)外學(xué)者圍繞人體骨骼的有限元分析方法提出了不少新的理論和方法，其中材料參數(shù)模擬方法一直是領(lǐng)域內(nèi)爭論的焦點(diǎn)[14-15]。作為下肢承重結(jié)構(gòu)的重要組成部分，脛骨結(jié)構(gòu)的生物力學(xué)分析研究有了長足的發(fā)展，但對脛骨有限元建模方法方面的相關(guān)研究較少，特別是脛骨建模過程中采用何種材質(zhì)賦值方法更符合臨床實際情況，目前尚無定論[16]。本研究分別通過灰度賦值法和均一賦值法建立兩種內(nèi)部幾何結(jié)構(gòu)及自身材料屬性參數(shù)均截然不同的脛骨有限元模型，模擬軸向壓縮載荷和扭轉(zhuǎn)載荷作用下脛骨結(jié)構(gòu)的生物力學(xué)響應(yīng)，并據(jù)此探索出一種更符合臨床實際情況的建模方案。

圖5　兩組脛骨模型在軸向壓縮與扭轉(zhuǎn)作用下的位移分布云圖（單位：mm）

灰度賦值建模手段比較簡單，通過參考既往文獻(xiàn)中的“灰度-彈性模量”經(jīng)驗公式即可完成材料屬性的指派，整個建模自動化程度相對較高，理論上比較符合個體化有限元分析的需求[7]。然而，這一途徑也存在局限性：首先，CT圖像質(zhì)量決定了材料賦值的分布，其中包括數(shù)據(jù)來源因素，如尸體標(biāo)本與活體骨骼的成像差異以及成像設(shè)備的參數(shù)差異等；其次，灰度賦值是根據(jù)模型單元體素和二維影像灰度像素的位置關(guān)系進(jìn)行聯(lián)系的，臨床應(yīng)用中的一些骨骼虛擬操作如截骨、內(nèi)固定等，將使骨骼相對原有的幾何形態(tài)和位置發(fā)生改變，從而影響灰度賦值的準(zhǔn)確性。

對于均一賦值法建模，既往文獻(xiàn)也報道了支持的觀點(diǎn)和研究結(jié)果，Vulovic等[17]利用均一賦值法分別對正常志愿者股骨三維有限元模型中的骨皮質(zhì)、骨松質(zhì)、ward三角區(qū)松質(zhì)骨等結(jié)構(gòu)進(jìn)行模擬生理載荷下的有限元分析，結(jié)果顯示，股骨結(jié)構(gòu)的生物力學(xué)響應(yīng)接近于尸體標(biāo)本測試數(shù)據(jù)。本研究中采用的均一賦值法，在建模階段和材料參數(shù)指派方面與灰度賦值法存在著不小的差異。均一賦值法要求分別建立皮質(zhì)骨、松質(zhì)骨和髓腔的實體模型，通過布爾運(yùn)算獲得僅由皮質(zhì)骨和松質(zhì)骨兩部分組成的具有腔體結(jié)構(gòu)的實體模型，在幾何結(jié)構(gòu)上更接近真實脛骨，并可分別對松質(zhì)骨和皮質(zhì)骨賦予相對恒定的彈性模量和泊松比。然而，均一賦值法建模過程同樣存在一些不足，包括幾何建模的準(zhǔn)確性存在人為差異、建模路線復(fù)雜以及類似股骨近端骨小梁特征性的形態(tài)分布在均一化處理中被忽略等等。

本研究結(jié)果顯示，在應(yīng)力分布方面，灰度賦值組脛骨模型的應(yīng)力分布相對均一賦值組更為均勻，主要原因之一可能是因為髓腔內(nèi)脂肪組織同樣存在灰度。盡管其賦予的彈性模量比較小，但低彈性模量單元填充整個髓腔后與骨性結(jié)構(gòu)構(gòu)成一個整體并共同承載傳遞的應(yīng)力效應(yīng)，這一機(jī)制可能在一定程度上造成了脛骨干應(yīng)力均勻分布的現(xiàn)象。在應(yīng)力、位移和剛度方面，灰度賦值法模型也表現(xiàn)出更高的抗壓和抗扭轉(zhuǎn)力學(xué)性能，這一結(jié)果也與大量髓腔單元填充產(chǎn)生的應(yīng)力分載效應(yīng)有關(guān)。而均一賦值組則呈現(xiàn)出脛骨中下段部位發(fā)生的應(yīng)力集中現(xiàn)象，因為此處為脛骨直徑由寬變窄的過渡部位，工程材料截面發(fā)生變化的位置恰恰是應(yīng)力集中的位置，所以這一應(yīng)力分布現(xiàn)象更加符合實際情況。位移分布方面，兩組模型之間呈現(xiàn)的分布云圖具有較高的一致性，通過分析位移分布方向和脛骨受力后變形的形態(tài)，均可觀察到這樣的一致性，表明兩種建模方法獲得的模型變形趨勢，無論是在軸向壓縮載荷還是扭轉(zhuǎn)載荷，均符合正常人體脛骨受力變形規(guī)律，這也驗證了本研究建模方法的有效性。而對于兩組模型脛骨結(jié)構(gòu)剛度的分析，如同其應(yīng)力、位移幅值差異一樣，也存在類似的差異，灰度賦值組脛骨模型的軸向壓縮剛度和扭轉(zhuǎn)剛度均大于均一賦值組。

從應(yīng)用角度來看，均一賦值方法有限元模型應(yīng)用的靈活性相對灰度賦值法建模更高，可用于復(fù)雜矯形外科手術(shù)的仿真分析及應(yīng)力應(yīng)變計算。李長有等[18]報道Chiari骨盆截骨術(shù)中外展肌力負(fù)荷對股骨頭軟骨生物力學(xué)影響的有限元分析，其采用的均一賦值骨盆有限元建模方法可實現(xiàn)對各種骨盆截骨術(shù)式的模擬分析，便于設(shè)置不同外展肌力等載荷條件，以獲得髖關(guān)節(jié)軟骨接觸壓應(yīng)力的分布結(jié)果。戴金龍等[19]應(yīng)用有限元方法比較分析傳統(tǒng)股骨轉(zhuǎn)子下封閉楔形外翻截骨術(shù)與改良“Z”形轉(zhuǎn)子間外翻截骨術(shù)的術(shù)后力學(xué)差異，其采用的也是均一賦值的材料屬性設(shè)定方法，在分析過程中可靈活設(shè)置股骨模型的截骨角度，以便于組間模型的對比。相比之下，灰度賦值法的有限元建模要求賦值前的網(wǎng)格模型空間坐標(biāo)位置需和原有CT圖像所覆蓋的范圍保持一致[20]，一旦骨骼結(jié)構(gòu)的網(wǎng)格模型空間位置發(fā)生變化（如旋轉(zhuǎn)截骨、內(nèi)外翻截骨等），通過經(jīng)驗公式賦予材料屬性將無法實現(xiàn)或形成錯誤的屬性賦值，采用均一賦值法可有效解決這一問題。

綜上，盡管兩種脛骨有限元模型在建模和賦值操作上有明顯差異，但最終通過仿真分析得到的結(jié)果在生物力學(xué)變化趨勢上具有一致性，從而驗證了有限元模型的有效性。而作為一種生物力學(xué)分析預(yù)測手段，均一賦值法構(gòu)建的脛骨有限元模型較之灰度賦值模型能更準(zhǔn)確地反映人體脛骨的生物力學(xué)特性，在影像數(shù)據(jù)要求和后續(xù)模型的修改應(yīng)用方面也更具優(yōu)勢。

參考文獻(xiàn)

[1]何旅洋,張志強(qiáng),鄭百林.外骨骼框架生物力學(xué)設(shè)計及強(qiáng)度分析[J].醫(yī)用生物力學(xué), 2014, 29(6): 504-510.

[2]Keyak JH, Rossi SA, Jones KA, et al. Prediction of femoral fracture load using automated finite element modeling [J]. J Biomech, 1998, 31(2): 125-133.

[3]Wang X, Sanyal A, Cawthon PM, et al. Prediction of new clinical vertebral fractures in elderly men using finite element analysis of CT scans [J]. J Bone Miner Res, 2012, 27(4): 808-816.

[4]汪方,石杜芳,王秋根,等.基于冰凍切片的人體骨盆有限元模型的建立與初步驗證[J].復(fù)旦學(xué)報(醫(yī)學(xué)版), 2010, 37(4): 384-390.

[5]Dall'Ara E, Luisier B, Schmidt R, et al. A nonlinear QCT-based finite element model validation study for the human femur tested in two configurations in vitro [J]. Bone, 2013, 52(1): 27-38.

[6]徐志才,胡廣洪,黃振宇,等.脛骨模型對膝關(guān)節(jié)有限元分析結(jié)果影響的探討[J].中國數(shù)字醫(yī)學(xué), 2014, 9(4): 69-72.

[7]張國棟,廖維靖,陶圣祥,等.股骨頸有限元分析的賦材料屬性方法探討及有效性驗證[J].中國組織工程研究與臨床康復(fù), 2009, 13(52): 10263-10268.

[8]Chen G, Wu FY, Liu ZC, et al. Comparisons of node-based and element-based approaches of assigning bone material properties onto subject-specific finite element models [J]. Med Eng Phys, 2015, 37(8): 808-812.

[9]Lengsfeld M, Schmitt J, Alter P, et al. Comparison of geometry-based and CT voxel-based finite element modeling and experimental validation [J]. Med Eng Phys, 1998, 20(7): 515-522.

[10]Completo A, Duarte R, Fonseca F, et al. Biomechanical evaluation of different reconstructive techniques of proximaltibia in revision total knee arthroplasty: an in-vitro and finite element analysis [J]. Clin Biomech, 2013, 28(3): 291-298.

[11]Gray HA, Taddei F, Zavatsky AB, et al. Experimental validation of a finite element model of a human cadaveric tibia [J]. J Biomech Eng, 2008, 130(3): 031016.

[12]楊鋒平,羅金恒,張華,等.金屬延性斷裂準(zhǔn)則精度的評價[J].塑性工程學(xué)報, 2011, 18(2): 103-106.

[13]談梅蘭,武國玉,梁福祥.基于ABAQUS的連桿疲勞分析[J].中國機(jī)械工程, 2013, 24(5): 634-638.

[14]蘇佳燦,張春才,陳學(xué)強(qiáng),等.骨盆及髖臼三維有限元模型材料屬性設(shè)定及其生物力學(xué)意義[J].中國臨床康復(fù), 2005, 9 (2): 71-73.

[15]Taddei F, Pancanti A, Viceconti M. An improved method for the automatic mapping of computed tomography numbers onto finite element models [J]. Med Eng Phys, 2004, 26(1): 61-69.

[16]Austman RL, Milner JS, Holdsworth DW, et al. The effect of the density-modulus relationship selected to apply material properties in a finite element model of long bone [J]. J Biomech, 2008, 41(15): 3171-3176.

[17]Vulovic S, Korunovic N, Trajanovic M, et al. Finite element analysis of CT based femur model using finite element program PAK [J]. J Serb Soc Computat Mech, 2011, 5(2): 160-166.

[18]李長有,王禹祥,范廣宇. Chiari骨盆截骨術(shù)中外展肌張力對股骨頭軟骨生物力學(xué)影響的三維有限元分析[J].醫(yī)用生物力學(xué), 2007, 22(1): 79-83.

[19]戴金龍,李炫升,梁蘭麗,等.改良式股骨轉(zhuǎn)子間外翻截骨術(shù)的生物力學(xué)研究:三維有限元素分析[J].中華創(chuàng)傷骨科雜志, 2007, 9(6): 510-514.

[20]Synek A, Chevalier Y, Baumbach SF, et al. The influence of bone density and anisotropy in finite element models of distal radius fracture osteosynthesis: evaluations and comparison to experiments [J]. J Biomech, 2015, 48(15): 4116-4123.

（本文編輯：白朝暉）

基礎(chǔ)研究

A comparison between two finite element analysis methods of human tibia via different modeling patterns

TANG Yongtao*, WEI Siqi, CHENG Gang, ZHAO Jianhua, CAI Fanfan, OUYANG Hanbin. *Orthopaedic Department I, Shenzhen Guangming New District People's Hospital, Shenzhen, Guangdong 518106, China

【Abstract】Objective To investigate the effects of two modeling approaches on biomechanical response of finite element models of human tibia. Methods A group of data of a volunteer's CT scan was obtained for 3D reconstruction of tibia, and the solid tibial model with cavity structure was then reconstructed by isolating medullary canal and cancellous bone. Material properties in 2 models were assigned via homogeneous and inhomogeneous modulus, identical load case and boundary conditions were applied to simulate the axial compression and torsion loading on the tibia during the mid-stance phase of gait circle simultaneously. Distribution of stress and displacement, stiffness of axial compression and torsion were compared between two models. Results Regarding to the displacement, two models presented similar distribution on the bone surface. Under the two load cases, inhomogeneous model resulted in lower maximum stress and maximum displacementcomparing to homogeneous model. Similarly, inhomogeneous group presented higher axial compression stiffness and torsional stiffness than homogeneous group under the same load case. A physiological stress concentration was observed on the distal part of tibia diaphysis in the homogeneous model, which relatively matched the clinical situation. Conclusions Both of the two modeling patterns could be available for finite element analysis of tibia models. However, homogeneous modeling was associated with higher accuracy for prediction of biomechanical response and flexibility in application compared to inhomogeneous modeling.

【Key words】Tibia; Bone remodeling; Models, anatomic; Finite element analysis; Material assignment; Stress, mechanical; Biomechanics; Computer simulation

中圖分類號：R356.1

文獻(xiàn)標(biāo)識碼：A

文章編號：1674-666X(2016)01-026-07

DOI：10.3969/j.issn.1674-666X.2016.01.005

收稿日期：（2015-12-03；修回日期：2016-01-12）