腕關節有限元骨性建模及力學分析

張浩，朱建民，馬南，解品亮，鐘務學，徐靈軍，趙海洋

(1.江蘇大學臨床醫學院，江蘇 鎮江212001;2.上海市徐匯區中心醫院骨科，上海200031)

手是由骨、關節、韌帶、肌肉和肌腱組成的，發揮手的任何功能都必須依賴于這些結構共同而復雜的配合來完成。只有深入地探索和研究其生物力學原理及各種疾病的發生、發展原理，手外科的整體理論深度和治療水平才有可能得到進一步發展［1]，從而更好地應用于臨床。本研究探索在個人計算機平臺上依據活體CT 圖像結合解剖資料，構建腕關節三維有限元模型的可行性及其應用價值。

1 材料與方法

1.1 影像設備與圖像處理軟件

儀器:荷蘭Philips64 排螺旋CT;計算機硬件:Dell latitude D630 工作站;CPU:Intel(R)Core(TM)i3504@ 3. 07 GHz;硬盤:320G;內存:12G;顯卡:NVIDIA GeForce310;相關軟件:Mimics10. 01，ANSYS10.0，GEOMAGIC。

1.2 方法

1.2.1 模型的建立 CT 圖像采集于1 名29 歲健康男性志愿者(圖1)，無腕部疾病史及外傷史，使用64層螺旋CT 掃描機對被試者上至尺橈骨中段，下至近端指骨進行薄層掃描。掃描條件:掃描電壓120 kV，層厚0.67 mm，層間隔0.67 mm，掃描圖像以DICOM格式保存并刻盤，并導入三維醫學圖像建模軟件Mimics 進行腕部骨骼三維立體模型重建。

圖1 腕部64 排128 層螺旋CT 掃描圖像Fig 1 Scanning images by 64 row 128 layer spiral CT

本研究主要包括骨和軟骨組織三維有限元的研究。從CT 掃描圖像中可以很清楚得到骨與軟骨的輪廓分界，人工自由選擇比較合適的單位擬定為軟骨數據，其余部分擬定為骨骼數據(圖2)。測量CT 掃描圖像中軟骨輪廓，并參考文獻［2]中軟骨厚度的數據，設置軟骨層厚度為1.2 mm。將CT 掃描得到的圖像以DICOM 3.0 格式刻盤，導入到三維醫學圖像建模軟件Mimics 10.01，根據骨與軟骨輪廓分界將骨與軟骨分割，然后建立實體三維模型，并以IGES 格式導出。并再次導入到ANSYS 軟件中，對模型進行網格化，并對網格化后的模型進行網格優化(圖3)獲得不包括軟骨及其附屬韌帶等的完整腕部骨骼有限元模型(圖4a)，并在此模型中根據軟骨層厚度進行手繪各軟骨面(圖4b)，同時記錄各骨模型的單元位置及單元編號，用于其后的材料屬性的賦值，以及接觸分析中接觸對的定義。最終得到完整的腕關節整體三維有限元模型，包括近端尺橈骨、腕骨、掌骨、近節指骨及各軟骨面。

圖2 參照CT 掃描圖像選取節點Fig 2 Selecting node according to the CT scan image

圖3 舟骨模型的網格劃分及網格優化Fig 3 Griddivision and grid optimization of scaphoid model

圖4 腕部整體有限元模型Fig 4 The whole finite element model

1.2.2 模型的材料屬性 生物組織屬于各向異性的非線性體，而且腕部結構復雜，影響因素眾多。本研究著重討論靜態載荷下的腕舟骨的受力情況(主要為舟、月骨與橈骨接觸面的接觸應力分布情況)，同時參考文獻方法，模型中僅考慮了骨與軟骨材料［3-5]，以減少計算量。

1.2.3載荷設置與邊界條件 將模型中的橈骨和尺骨近端在X/Y/Z 3 個方向完全約束［3]，對近節指骨各施加軸向共100 N 的壓力載荷(圖5)［6]。為了更好地模擬腕管載荷傳遞情況，本模型將相鄰腕骨、橈骨與尺骨、橈骨與腕骨、尺骨與腕骨、腕骨與掌骨各關節處的軟骨單元設置為接觸單元。根據解剖實驗數據得到載荷作用下各骨周圍可能與之接觸的骨骼分布情況，設置各軟骨接觸對?？紤]軟骨之間的摩擦系數較小，且關節中的組織滑液有減小摩擦的作用，因此本文沒有考慮摩擦因素的影響。

圖5 有限元模型約束Fig 5 The finite element model been restricted of the wrist joint with prosthesis in ANYSYS

2 結果

建立了包括橈、尺骨遠端，腕骨，掌骨，近端指骨，各關節軟骨面，三角纖維軟骨復合體等結構的三維有限元模型，成功模擬了掌骨受軸向壓力情況，并得到各腕骨應力分布圖。見圖6。尤其真實有效模擬軸向壓力加載下舟狀骨應力分布最為集中地分布于腰部(圖7)，這明確解釋了臨床工作中舟狀骨骨折多見于腰部，并且有效驗證了本模型的真實性及有效性。

圖6 模型整體和局部應力云圖Fig 6 the whole and local stress nephogram of model

圖7 舟骨應力云圖Fig 7 The stress nephoqram of scaphoid

3 討論

有限元素法是一種有效的離散化數值計算方法，Courant 于1943年首創了有限元單元法，隨著醫學與工程力學的結合，有限元分析作為一種新的生物力學研究手段，逐漸被應用到醫學領域，成為目前生物力學研究中先進有效的方法。Brekelmans等［7]和Rybicki 等［8]在1972年首次將有限元法應用于骨科生物力學的研究。

腕關節各腕骨形態極不規則，不能用簡單的圖形或公式計算出來，目前的建模方法通常是通過CT、X 線、磁共振、超聲等其他醫學影像資料得到研究對象的連續切面的二維圖像，然后用計算機三維重建技術來得到三維模型。本研究利用64 排螺旋CT 進行薄層掃描，基于Mimics 軟件首先繪制出三維模型，然后放入GEOMAGIC 軟件中進行網格劃分，再導入Mimics 中，然后使用布爾運算將其與相接觸的骨骼相減，然后再一次導入GEOMAGIC 中進行光滑處理和重新的網格劃分(此時的模型由于進行了布爾運算，與之前模型有了很大的差別，故需要重畫網格)，然后再導入Mimics 中，進行整體放大1.02 倍的操作，這樣就可以解決之前模型各組件之間完全沒有重合的問題。通過以上操作后模型不僅有重合部分，而且重合部分很小［9-10]。

本研究首次將尺骨和橈骨遠端、腕骨以及掌骨、近端指骨、各關節軟骨建立成三維有限元模型，模擬軸向荷載從掌骨經腕骨傳遞到遠端尺骨和橈骨的受力情況。本模型對掌骨受到100 N 軸向壓力荷載時各腕骨接觸應力的分布情況進行了三維有限元分析，并獲得應力分布云圖。結果表明，模型中應力集中分布區域，正是我們臨床工作中骨折常見部位，如腕舟骨，應力主要分布于腰部，臨床中舟骨骨折多見于腰部.這對我們進一步探討外荷載作用下整個腕部結構的力學行為提供了一個可操作的平臺和可靠的理論依據。

然而本研究在約束的施加方面發現新問題。首先，軟骨最多只能在3 個節點上施加約束，不然會導致計算不收斂，具體原因暫不清楚。其次，腕骨、掌骨、指骨，初次試驗于每塊骨骼的下端施加約束，在計算后發現，無論施加多大的力(100 ～1000 N)，模型的最大變形和最大應力，以及應力的分布都是相同的。我們推測其原因是，模型滿足本構關系，即應力應變的乘積為常數。在約束較多的情況下，模型只能有非常小的位移。本研究中位移為0.98 mm，在這樣的情況下，應力值相應確定，而施加軸向壓力超出此范圍，則無法計算，是否可以說明在此情況下已導致骨折，還有待于進一步試驗，也是我們下一步研究的一部分。

腕部功能的完整描述還應該同時包括肌腱和韌帶的作用，而本模型最大的不足就是因數據提取困難，未進行肌腱及韌帶建模。隨著實驗手段和測試技術的發展，通過建立并利用腕部結構三維有限元模型，計算和分析腕部整體結構或各組成部分如腕骨、韌帶、肌肉及肌腱等組織結構在各種力荷載作用下的應力應變及位移情況，將會成為腕部生物力學的一個重要研究方向。

另外本研究在接觸算法中應用MPC 算法，使得計算更具精度。bonded 接觸算法使用了懲罰方法，由于接觸剛度和穿透，可能會影響結果的精度，即使是線性問題，通常也需要秩代，這是因為使用了接觸剛度。MPC 算法不需要操作者手工定義MPC 方程，提高了求解效率，不再需要通過多次嘗試來提高和保證求解精度，求解方程不需要秩代，不僅可以約束平移自由度，而且可以約束旋轉自由度，使模型后續的有限元分析更具優勢。本研究建立的腕關節三維有限元模型的應用價值在于后續導入有限元軟件可以進行腕部相關的生物力學分析、計算，并作為腕關節病理有限元模型的構建基礎，為腕關節外科提供一定依據。

有限元的研究隨著計算機技術的進步日益凸顯成熟，這也將推動腕關節有限元模型的研究，本文做了這方面的初步嘗試，難免存在不足，例如因數據提取困難，未進行肌腱及韌帶建模，另外由于為單個樣本研究，普遍的統計學意義有待于驗證。

［1]顧玉東，王澍寰，侍德，等.手外科學［M].上海:上?？茖W技術出版社，2004:85.

［2]Streekstra GJ，Braskamp P，vander Leij C，et a1.Cartilage thickness measurement in the sub-millimeter range［C].∥Barillot C，Haynor DR，Hellier P，eds. Medical Image Computing and Computer-assisted Intervention. Berlin:Springer-Verlag，2004:950-958.

［3]Carrigan SD，Whiteside RA，Pichora DR，et al. Development of a three-dimensional finite element model for carpal load transmission in a static neutral posture［J]. Ann Biomed Eng，2003，31(6):718-725.

［4]Anderson DD，Deshpande BR，Daniel TE，et al. A threedimensional finite element model of the radiocarpal joint:distal radius fracture step-off and stress transfer［J]. Iowa Orthop J，2005，25:108-117.

［5]Dacho AK，Baumeister S，Germann G，et al. Comparison of proximal row carpectomy and midcarpal arthrodesis for the treatment of scaphoid nonunion advanced collapse(SNACwrist)and scapholunate advanced collapse(SLAC-wrist)in stage II［J]. J Plast Reconstr Aesthet Surg，2008，61(10):1210-1218.

［6]Ural A. Prediction of Colles＇fracture load in human radius using cohesive finite element modeling［J]. J Biomech，2009，42(1):22-28.

［7]Brekelmans WA，Poort HW，Slooff TJ，et al.A new method to analyse the mechanical behaviour of skeletal parts［J].Acta Orthop Scand，1972，43(5):301-317.

［8]Rybicki EF，Simonen FA，Weis EB Jr.On the mathematical analysis of stress in the human femur［J]. J Biomech，1972，5(2):203-215.

［9]方馳華，周五一，黃立偉，等. 虛擬中國人女性一號肝臟圖像三維重建和虛擬手術的切割［J]. 中華外科雜志，2005，43 (11):748-752.

［10]劉文軍，鐘世鎮. 虛擬中國人女性一號松質骨圖像數據的配準與三維重建［J]. 中國臨床解剖學雜志，2004，22 (4):380-383.