人體足踝部有限元模型的建立及有效性分析

吳愷楊茂偉*都承斐

（1.中國醫(yī)科大學附屬第一醫(yī)院骨科，沈陽110001；2.北京航空航天大學生物與醫(yī)學工程學院，北京100191）

隨著工業(yè)、礦業(yè)及交通運輸業(yè)的發(fā)展，足踝部的損傷不斷增加，越來越受到骨科醫(yī)師及運動醫(yī)學專家的重視，人們對足踝部的生物力學研究也不斷深入?；铙w組織研究能得到較為精確的數(shù)據(jù)，但因倫理學的原因受到限制，而對足踝部的傳統(tǒng)研究方法是進行尸體標本的生物力學研究。尸體標本生物力學研究目前被認為是金標準，但也有其局限性，首先很難獲得尸體標本；其次尸體標本生物力學研究只能獲得標本模型的整體信息，標本模型內(nèi)部的區(qū)域性力學特性卻很難獲得；再次，尸體標本模型失去隨意肌的控制，與足踝部的自然狀態(tài)也存在差別；此外模型的制備及測量也存在一定誤差。隨著計算機技術(shù)和有限元理論的不斷發(fā)展，人們開始大量使用數(shù)字模型和有限元法分析復雜的結(jié)構(gòu)。與以往的生物力學實驗相比，有限元方法可以建立高度幾何相似及物理相似的有限元模型，既可以反映區(qū)域性力學特性，又可以獲得整體信息；既可以進行精確的數(shù)字分析，又可以進行形象的、直觀的定性研究。本研究基于CT數(shù)據(jù)建立具有高度幾何相似性的足踝部三維有限元模型，并對其有效性進行驗證，為臨床及生物力學研究提供新的研究手段。

1 資料與方法

1. 1 數(shù)據(jù)收集

男性志愿者1名，28歲，身高172 cm，體重60 kg，足部無外傷史，X線檢查排除足部腫瘤、畸形等病變。采用中國醫(yī)科大學附屬第一醫(yī)院影像中心的Brilliance 64排iCT，掃描參數(shù)：電壓120 kV，電流強度372.5 mA，螺距0.297，層厚0.625 mm，矩陣512×512，曝光時間0.4 s/r，重建層厚0.9 mm，重建方法B+YC，重建間隔0.45 mm。對志愿者右足自踝關(guān)節(jié)上10 cm脛腓骨遠端向下掃描至足底，掃描時右足保持中立位，共獲得足部橫斷面CT圖片516張，以DICOM格式存貯，刻錄成光盤。

1. 2 有限元模型建立

CT掃描的數(shù)據(jù)導入三維重建軟件Mimics10.0，通過閾值分開骨組織及軟組織，建立足踝部骨骼幾何模型；根據(jù)解剖圖譜[1]及文獻[2]在骨骼間的相應(yīng)附著點建立韌帶模型，骨骼以STL格式輸出，韌帶以IGES格式輸出。然后將骨骼導入計算機輔助設(shè)計軟件Geomagic Studio10.0中，對模型進行除噪點、平滑，根據(jù)各關(guān)節(jié)面的幾何形狀，利用軟件加厚功能在各骨面生成關(guān)節(jié)軟骨，最后擬合曲面，骨骼及關(guān)節(jié)軟骨以IGES格式輸出。將骨骼、關(guān)節(jié)軟骨及韌帶導入Ansys 13.0的Workbench模塊進行材料賦值、網(wǎng)格劃分、設(shè)置邊界條件等處理，最后進行計算分析（圖1）。

圖1 足踝有限元模型

1. 3 材料參數(shù)

骨性結(jié)構(gòu)和軟骨均模擬為各向同性的線彈性材料，韌帶使用只承受拉應(yīng)力單軸梁單元來模擬。骨、軟骨、韌帶的材料屬性由參考文獻[3,4]確定，見表1。

表1 材料屬性

1. 4 劃分網(wǎng)格

1. 5 設(shè)置邊界條件

本研究所建立的軟骨較厚，與骨骼有交疊現(xiàn)象，通過布爾運算減去交疊部分再和骨骼邊界設(shè)置為“綁定（bonded）”，二者就緊密地結(jié)合為一體了。韌帶與骨骼之間的接觸也設(shè)置為綁定，但如果韌帶與骨骼之間距離過大則無法綁定，調(diào)整韌帶止點，使之與骨骼綁定可靠。關(guān)節(jié)軟骨之間采用有摩擦的面面接觸，軟骨面的間隙＜0.1 mm，關(guān)節(jié)之間由于關(guān)節(jié)滑液的存在，使得其間摩擦力非常小，依據(jù)參考文獻[5]設(shè)置摩擦系數(shù)為0.01，用以模擬關(guān)節(jié)軟骨之間的滑動特性。

由于Workbench劃分網(wǎng)格功能十分強大，故對骨骼、軟骨網(wǎng)格的劃分采用“自動生成網(wǎng)格（automatic）”即可，網(wǎng)格尺寸選擇為“粗糙（coarse）”即可，這樣劃分的網(wǎng)格基本滿意，且節(jié)點和單元數(shù)量不多，在保證計算精度的前提下，提高了計算的收斂速度。對線體網(wǎng)格的劃分仍可采用自動生成網(wǎng)格功能，但劃分的單元數(shù)量選擇為“1”，這樣每條韌帶只用一個單元模擬，使得計算效果更好。

1. 6 有限元模型的驗證

參照Anderson等[6]實驗所設(shè)置的的邊界條件，模擬人體站立狀態(tài)下的力學傳遞，對模型脛、腓骨下端的上截面施加600 N垂直壓縮載荷，固定跟骨，約束距骨，測量踝關(guān)節(jié)下關(guān)節(jié)面的接觸應(yīng)力、接觸面積，并將結(jié)果與前者進行對照。同時測量踝關(guān)節(jié)周圍主要韌帶的位移（拉伸長度），對踝關(guān)節(jié)的生物力學反應(yīng)進行自身驗證。

2 結(jié)果

本研究建立的正常人足踝部三維有限元模型形態(tài)還原性好，重建效果較理想，能任意旋轉(zhuǎn)，可獲得詳細、滿意的三維信息。模型計算脛骨下關(guān)節(jié)面最大應(yīng)力（3.74 MPa）及應(yīng)力分布均與以往研究結(jié)果具有相似性[6,7]（圖2、3，表2），認為模型有效。模型測得踝關(guān)節(jié)周圍主要韌帶的位移（拉伸長度）見圖4，其中脛跟韌帶位移最大（5.24 mm），脛舟韌帶（5.04 mm）、跟腓韌帶（2.05 mm）次之，下脛腓聯(lián)合前韌帶位移最?。?.92 mm）。

圖2 Anderson等對兩個踝關(guān)節(jié)脛骨下關(guān)節(jié)面在標本和有限元模型中分別測量的應(yīng)力分布結(jié)果示意圖

圖3 本研究建立的有限元模型脛骨下關(guān)節(jié)面應(yīng)力分布結(jié)果示意圖

表2 脛骨下關(guān)節(jié)面應(yīng)力分布對照

3 討論

自從Courant于1943[8]年提出有限元概念以來，有限元理論及其應(yīng)用得到了迅速的發(fā)展，人們開發(fā)了多種逆向工程軟件和計算機輔助設(shè)計軟件，并且不斷更新軟件功能，以彌補軟件本身的局限性。同時，人們也嘗試不同方法建模，研究不同材料的屬性以及劃分網(wǎng)格方法，探索不同的邊界設(shè)置條件及加載方法，以達到更好的仿真效果。隨著技術(shù)的不斷更新，生物數(shù)字仿真學也在不斷地進步，但是由于生物組織幾何上的不規(guī)則以及非均勻性、各向異性、非線性的材料特性，并且與應(yīng)變速率、加載歷史等因素具有相關(guān)性，目前很難精確模擬。足踝部骨骼形態(tài)極不規(guī)則，數(shù)量多、體積小，關(guān)節(jié)多且曲面復雜，韌帶、肌腱繁多，且具有復雜的運動學與動力學特性，給有限元建模及模擬帶來更多的麻煩。

為了克服足踝部有限元分析的各種困難，本研究做了不同的嘗試。首先，采用多條線體模擬韌帶，在Mimics中利用軟件輔助設(shè)計（MedCAD）功能，根據(jù)解剖圖譜[1]及文獻[2]在骨骼間生成多條線單元模擬韌帶，避免了單一線體或彈簧模擬韌帶產(chǎn)生的應(yīng)力集中，也避免彈簧所需要的預加載。韌帶以IGES格式輸出，可反復導入有限元分析軟件，避免調(diào)整模型時反復制備線體，減少了工作量。韌帶使用只受拉不受壓的單軸梁單元來模擬其特有的力學屬性，使用一個單元對線體進行網(wǎng)格劃分，在不影響計算結(jié)果的同時優(yōu)化了收斂效果。

其次，在計算機輔助設(shè)計軟件Geomagic中根據(jù)各關(guān)節(jié)面的幾何形狀，利用軟件加厚功能在各骨面生成實體模擬關(guān)節(jié)軟骨。生成的各軟骨厚度均勻，構(gòu)造曲面片時采用單層結(jié)構(gòu)，這樣在與對應(yīng)的骨骼做布爾運算后不會產(chǎn)生畸形結(jié)構(gòu)，網(wǎng)格劃分時也不會產(chǎn)生不良單元。在Geomagic中可清晰查看軟骨間的接觸，調(diào)整軟骨，使相鄰軟骨之間只接觸而幾乎無交疊現(xiàn)象（交疊＜0.1 mm），提高計算的收斂率和有效性。參考文獻[5]設(shè)置關(guān)節(jié)間摩擦系數(shù)為0.01，用以模擬關(guān)節(jié)軟骨之間的滑動特性，保證了關(guān)節(jié)的有效滑動，使其更接近自然關(guān)節(jié)的狀態(tài)。

圖4 踝關(guān)節(jié)周圍主要韌帶的位移

由于受逆向工程軟件和計算機輔助軟件自身局限性，加之建模中也會受原始影像資料精確程度、模型重建方法、網(wǎng)格劃分水平、材料屬性賦值、研究者經(jīng)驗和熟練程度等因素影響，有限元計算的結(jié)果可能會與力學模型結(jié)果之間存在差別。因此，建立的有限元模型需要進行有效性驗證，一般是將其與力學數(shù)據(jù)直接比較。本模型計算中立位踝關(guān)節(jié)垂直加載應(yīng)力主要分布于脛骨下關(guān)節(jié)面中部及前外側(cè)，最大應(yīng)力2.49～3.74 MPa。Anderson等[6]研究顯示兩具尸體標本模型脛骨下關(guān)節(jié)面最大應(yīng)力為3.69 MPa和2.92 MPa，兩個有限元模型脛骨下關(guān)節(jié)面最大應(yīng)力為3.74 MPa和2.74 MPa。劉清華等[7]研究顯示有限元模型脛骨下關(guān)節(jié)面最大應(yīng)力為3.97 MPa。以上研究脛骨下關(guān)節(jié)面應(yīng)力分布結(jié)果與本研究相似，說明本模型有效。

模型同時也測量了踝關(guān)節(jié)周圍四條主要韌帶的位移（拉伸長度）。從文獻[9]可知，踝關(guān)節(jié)周圍韌帶位移與踝關(guān)節(jié)及距下關(guān)節(jié)有關(guān)。踝關(guān)節(jié)為屈戌關(guān)節(jié)，距骨體及滑車前寬后窄，負重時下脛腓關(guān)節(jié)有分離趨勢，下脛腓聯(lián)合前后韌帶十分堅韌，能夠限制此趨勢，維持踝穴的穩(wěn)定、進而維持踝關(guān)節(jié)的穩(wěn)定，因此下脛腓聯(lián)合前韌帶位移很小，本文測得的下脛腓聯(lián)合前韌帶位移最小，亦證實如此。由于外踝較內(nèi)踝較長并靠后，故踝軸向后下方成角，當人體中立位負重時，小腿有外翻趨勢。距下關(guān)節(jié)的前、中、后三個關(guān)節(jié)繞同一單軸運動，類似一斜鉸鏈，由后外下方斜向前內(nèi)上方，其運動軸無論在正位或側(cè)位均與矢狀位相交為一定角度。由于這一斜軸關(guān)系，站立位負重時進一步加重小腿的外翻，同時小腿亦出現(xiàn)了輕微的內(nèi)旋，Liacouras和Wayne[10]等研究及Marqueen[11]等研究模擬中立位負重均顯示脛骨有0°～0.5°的內(nèi)旋。因此外側(cè)的跟腓韌帶位移小于內(nèi)側(cè)韌帶，脛跟韌帶位移大于脛舟韌帶。本文計算的踝關(guān)節(jié)周圍四條主要韌帶的位移符合踝關(guān)節(jié)的生物力學反應(yīng)，進一步說明了本模型的有效性，可以用于臨床及相應(yīng)的生物力學研究。

目前對足踝部生物力學的研究雖然取得了很大的研究成果和進展，但仍然只是做近似的模仿。本研究基于CT序列圖像，對建模方法進行了嘗試和創(chuàng)新，建立正常人體足踝部三維有限元模型，較真實、準確地模擬了足踝部的解剖形態(tài)，同時具有良好的生物力學特性。不足之處在于對模型進行了簡化，僅建立了后足模型，且未考慮到肌肉等動力性結(jié)構(gòu)，以后的工作需要進一步研究，使其成為足踝部生物力學研究的更好平臺。

[1]高士濂.實用解剖圖譜.下肢分冊.第2版.上海:上?？茖W技術(shù)出版社,2004:275-281.

[2]Inhauser CW,Siegler S,Udupa JK,et al.Subject-specific models of the hindfoot reveal a relationship between morphology and passive mechanical properties.J Biomech,2008,41(6):1341-1349.

[3]Cheung JT,Zhang M,Leung AK,et al.Three-dimensional finite element analysis of the foot during standing:a material sensitivity study.J Biomech,2005,38(5):1045-1054.

[4]Bandak FA,Tannous RE,Toridis T.On the development of an osseo-ligamentous finite element model of the human ankle joint.Int J Solids Struct,2001,38:1681-1697.

[5]Anderson DD,Goldsworthy JK,Shivanna K,et al.Intra-articular contact stress distributions at the ankle throughout stance phase-patient-specific finite element analysis as a metric of degenerationpropensity.Biomech Model Mechanobiol,2006,5(2-3):82-89.

[6]Anderson DD,Goldsworthy JK,Li W,et a1.Physical validation of a patient-specific contact finite element model of the ankle.J Biomech,2007,40(8):1662-1669.

[7]劉清華,余斌,金丹,等.正常人足踝部有限元模型的構(gòu)建研究.中華創(chuàng)傷骨科雜志,2010,12(2):174-177.

[8]Courant R.Variational method for solution of problems of equilibrium and vibrations.Bull Am Math Soc,1943,49:1-23.

[9]郭世紱.骨科臨床解剖學.泰安:山東科技出版社,2001:989.

[10]Liacouras PC,Wayne JS.Computational modeling to predict mechanical function of joints:application to the lower leg with simulation of two cadaver studies.J Biomech Eng,2007,129(6):811-817.

[11]Marqueen T,Owen J,Nicandri G,et a1.Comparison of the syndesmotic staple to the transsyndesmotic screw:a biomechanical study.FootAnkle Int,2005,26(3):224-230.