基于CT和MRI膝關節有限元模型建立與不同載荷下生物力學特性分析

杭聲琪，張文韜，王泰淇，任建澄，李富超，趙改平，馬童

1.上海理工大學醫療器械與食品學院，上海200093；2.上海理工大學環境與建筑學院，上海200093；3.同濟大學附屬楊浦區中心醫院骨科，上海200090

前言

膝關節是人體下肢活動的樞紐關節，其解剖結構和運動功能比較復雜，脛股關節和髕股關節通過肌肉、肌腱和韌帶連接組成，保證膝關節運動的穩定性和靈活性。膝關節在運動中承受較大載荷，是易發生損傷和病變的部位，對膝關節生理病理和治療機制的研究一直備受醫生和學者們的關注。隨著計算機技術和數值計算方法的快速發展，采用有限元方法研究膝關節的生理功能、損傷預防和治療成為重要的輔助手段。

采用有限元方法研究不同載荷下骨骼、韌帶及半月板的生物力學特性一直是眾多臨床和科研工作者研究的熱點。如Bendjaballah等［1］研究正常膝關節施加5、10和15 N·m 內外翻力矩時，股骨旋轉角度和前后交叉韌帶張力的生物力學特性。Li 等［2］建立了正常膝關節的有限元模型，在脛骨上加載從0到10 N·m內旋外旋力矩，研究脛骨內旋和外旋角度的變化。Gardiner等［3］建立膝關節屈曲0°、30°、60°的有限元模型，施加10 N·m 外翻力矩，研究內側副韌帶的應變情況。Liu 等［4］在股骨上施加10 N·m 內翻外旋力矩，研究后交叉韌帶重建的膝關節內翻外旋角度。Vairi等［5］建立膝關節正常、韌帶缺失和重建的有限元模型，在股骨上施加1 150 N 壓縮力、134 N 抽屜力以及10 N·m 外翻力矩，分析前交叉韌帶斷裂與重建后膝關節的力學特性。Luczkiewicz 等［6］建立了正常膝關節有限元模型，在股骨上施加1 000 N載荷，研究擠壓變形后不同橫截面半月板的應力變化情況等。膝關節有限元模型進行網格劃分時，國內外學者多數采用四面體網格進行劃分［7-9］。本文基于人體膝關節的CT 和MRI 圖像建立幾何模型，采用六面體單元網格建立包括骨骼、軟骨、半月板和韌帶在內的三維有限元模型。施加載荷，進行膝關節有限元模型的有效性驗證與分析，得出軟骨和半月板的接觸應力、接觸面積、Von Mise 應力及位移情況，研究加載不同大小的內翻、外翻力矩時股骨內翻、外翻角度、半月板和韌帶的應力變化情況等。

1 材料和方法

1.1 研究對象

選取一名膝關節健康女性志愿者，年齡52歲，體質量60 kg，使用CT 對膝關節直立時進行掃描，得到層厚0.75 mm、大小為512×512 像素的圖像，使用MRI 對同一膝關節進行掃描，得到層厚0.6 mm、大小為256×256像素的圖像。

1.2 模型與網格劃分

將CT 和MRI 圖像導入Mimics（Materialise NV,Leuven,Belgium），在CT 圖像中選取骨骼部位的閾值，進行區域生長自動分割骨骼區域，包括股骨、脛骨、髕骨和腓骨。由于在CT 圖像中軟組織顯示界限不清楚，因而在MRI 圖像中提取。在MRI 圖像T1序列中，根據解剖學位置和MRI 影像顯示，分割股骨軟骨、半月板、脛骨軟骨、腓骨軟骨、前交叉韌帶（ACL）、后交叉韌帶（PCL）、內側副韌帶（MCL）、外側副韌帶（LCL）、髕韌帶和髕骨上方相連的股四頭肌腱等。三維重建出各部分的幾何模型后，將三維模型導入Geomagic Studio 12（Raindrop Geomagic,Inc.）中進行進一步修整模型以及曲面的生成，重建的全膝關節幾何模型如圖1a 所示。利用hypermesh13.0進行六面體網格劃分時，從源面生成二維四邊形網格，然后沿著源面與目標面生成三維六面體網格，網格從源面映射到目標面。譬如股骨的六面體網格劃分示意圖如圖1b 所示。經過網格收斂性測試，將網格單元大小設置骨骼為2.0 mm，軟骨、韌帶、半月板為1.0 mm，單元類型為C3D8R，網格經檢查無穿透或交叉。膝關節六面體網格有限元模型如圖1c所示。

圖1 膝關節幾何模型和有限元模型Fig.1 Geometric model and finite element model of the knee joint

1.3 材料屬性

膝關節的股骨、脛骨和腓骨視為各向同性線彈性材料，各個關節軟骨也視為各向同性線彈性材料［10］，數值如表1所示。韌帶為各向同性的超彈性材料，采用Neo-Hooke 本構模型，其能量密度函數為Ψ=C1×(I1-3)，C1表示初始剪切模量，I1表示Cauchy-Green 應變張量的第一修正不變量［11］。半月板為橫觀各向同性彈性材料，即半月板所處平面內的力學性質與垂直面的不同［12］。半月板周向上彈性模量為Ec，徑向和軸向的彈性模量為Ea=Er，各方向泊松比和剪切模量數值如表1所示。內外側半月板前后角與脛骨平臺間通過韌帶相連，采用剛度為2 000 N/mm線性彈簧單元模擬［13］。內外側半月板前角間相連的橫韌帶，也采用線性彈簧單元模擬，彈簧單元的剛度為900 N/mm［13］。

表1 膝關節材料屬性Tab.1 Material properties of the knee joint

1.4 邊界條件和載荷

使用ABAQUS（Simulia,Providence,RI）軟件設置各結構的相互作用、邊界條件以及載荷，股骨軟骨和內外側半月板之間設置2個接觸對，股骨軟骨和內外側脛骨軟骨之間設置2個接觸對，內外側半月板和內外側脛骨軟骨之間設置2 個接觸對，共計6 個接觸對。各部位的接觸摩擦系數設為0.02，相互作用屬性設置為有限滑動，主表面與從表面之間采用默認的“硬接觸”。所有的軟骨與相應的骨骼之間設置綁定約束，韌帶的兩端與骨骼之間也設置綁定約束。股骨的6個自由度不受固定，將邊界條件設置在脛骨遠端和腓骨遠端，固定6個自由度。將股骨的網格設置耦合約束，耦合點設置在股骨兩髁的中心點。耦合點加載形式如下幾種載荷：1 000 N 壓縮力以及股骨前向134 N 抽屜力［18］，為了進行模型驗證；1 000 N 壓縮力和分別為5、10、15 N·m外翻力矩以及1 000 N壓縮力和分別為5、10、15 N·m 內翻力矩［14］，為了對比施加不同內翻和外翻力矩時的生物力學特性。

2 結果

2.1 膝關節接觸應力和接觸面積

膝關節施加1 000 N 軸向壓縮載荷和股骨后向134 N 抽屜力的作用，股骨軟骨、內外側半月板和內外側脛骨軟骨上的接觸應力（CPRESS）如圖2所示。膝關節接觸應力體現了兩個結構之間擠壓時接觸面上產生的應力。計算結果顯示，股骨軟骨上最大接觸應力值為4.47 MPa，在內側和外側的中心部位出現應力集中，最大接觸應力分布于內側。內外側半月板最大接觸應力值分別為3.25和2.83 MPa，內側半月板接觸應力分布于中后部，分布面積較大且均勻，外側半月板的接觸應力分布于后部，半月板前部的接觸應力相對較小。內外側脛骨軟骨的最大接觸應力值分別為2.70和2.52 MPa，位于中心和股骨軟骨與半月板接觸的地方，外側脛骨軟骨的接觸應力最大位置呈現在與股骨軟骨接觸的地方，其中內側脛骨軟骨接觸應力比外側脛骨軟骨的分布要廣泛。各個部位的接觸面積與文獻中結果［16-20］對比如圖3所示。在接觸面積方面，內外側脛骨軟骨的接觸面積分別為505.47和413.3 mm2。

圖2 膝關節股骨軟骨、半月板和脛骨軟骨的接觸應力Fig.2 Contact stresses of femoral cartilage,meniscus and tibial cartilage of the knee joint

圖3 膝關節接觸應力的對比Fig.3 Comparison of contact stresses of the knee joint

2.2 膝關節Von Mises應力和位移

膝關節的Von Mises 應力分布如圖4所示。股骨軟骨、內外側半月板和內外側脛骨軟骨最大Von Mises 應力分別為2.22、2.44、2.25、2.07 和1.64 MPa。分布上，股骨軟骨最大Von Mises 應力分布于內側部位，半月板主要分布于中部以及前后角位置。脛骨軟骨應力主要分布于與股骨軟骨接觸的中心位置，以及與半月板接觸位置。位移方面，膝關節股骨相對于脛骨前向上的位移為4.19 mm。

圖4 膝關節股骨軟骨、半月板和脛骨軟骨Von Mises應力分布Fig.4 Von Mises stress distribution of femoral cartilage,meniscus and tibial cartilage of the knee joint

2.3 膝關節不同內外翻力矩下應力結果

股骨內外翻傾角、半月板接觸應力、前后交叉韌帶應力和內外側副韌帶峰值應力隨不同力矩的變化分布如圖5所示。隨著內翻力矩從5、10增大到15 N·m，股骨內翻傾角分別為3.49°、4.48°、4.91°，隨著外翻力矩同樣線性增大，股骨外翻傾角分別為3.22°、3.62°、4.01°。從內翻到外翻過渡過程中，外側半月板應力漸增，而內側漸減，ACL 與PCL 應力呈現出先減少再增大的趨勢，MCL應力先漸減再趨于平穩，LCL應力呈現漸增的趨勢。

圖5 膝關節內外翻力矩下股骨傾角、半月板和韌帶峰值應力變化Fig.5 Femoral inclination and the changes in peak stresses of meniscus and ligament of the knee joint under varus and valgus torques

3 討論

人體膝關節三維有限元模型采用六面體網格劃分，在膝關節直立位時，施加1 000 N 壓縮力和134 N抽屜力，研究關節軟骨和半月板接觸應力和Von Mises 應力的分布情況，脛骨軟骨接觸面積以及股骨的前后位移情況，驗證了模型的有效性。

膝關節股骨軟骨、內外側半月板和內外側脛骨軟骨接觸應力CPRESS 的計算結果與前人體外試驗和有限元分析結果［16-20］的對比如圖3所示。股骨軟骨的接觸應力與Shriram 等［16］施加1 150 N 載荷計算的最大接觸應力4.64 MPa 相一致，與Bae 等［15］計算的股骨軟骨接觸應力的分布趨勢相同。半月板接觸應力與Pe?a等［17］施加1 150 N載荷和134 N脛骨前向抽屜力的計算結果3.68 和3.82 MPa 相比減少了11.7%和25.9%，原因在于半月板設置的材料屬性以及施加的載荷不同，比Zhu 等［18］施加1 000 N 載荷的數值2.55和2.47 MPa略大，由于其未考慮134 N抽屜力的作用。Seitz 等［19］和Bao 等［20］模擬內外側脛骨軟骨接觸應力分別為2.7 和3.2 MPa，3.09 和2.97 MPa，本文與其結果相近。接觸面積比Liu［21］計算的內側脛骨軟骨和外側脛骨的平均接觸面積515.29和390.29 mm2減少了1.9%和增大了5.9%，與Geeslin 等［22］體外實驗結果538 mm2減少了22.3%，由于施加載荷的不同以及個體膝關節的差異，而與前人的有限元結果相近，具有一定的準確性。股骨軟骨Von Mises 應力與Shriram 等［16］股骨軟骨應力為2.76 Mpa 相比減少了19.6%，內外側半月板Von Mises 應力與Pe?a 等［11］在1 150 N 載荷下內外側半月板應力分別為2.75 和1.62 MPa 相比減少10.9%和增加38.9%，內外側脛骨軟骨Von Mises 應力與Pe?a［17］等在1 150 N 壓縮載荷和134 N 抽屜力作用下內外側脛骨軟骨應力分別為2.19 和1.76 Mpa 相比減少了5.5%和6.8%，存在差別的原因是施加載荷的不同所致，整體Von Mises 應力接近。臨床上，抽屜力試驗用于檢查交叉韌帶損傷以及韌帶重建后的療效測評［23］，評價膝關節的穩定性。對于抽屜力的施加，可位于脛骨或股骨，方向前向或后向。本文選取在股骨上施加后向的134 N 抽屜力，在位移的數值上，Gabriel 等［24］計算的位移為4.0 mm，位移結果與前人研究相一致。

內側副韌帶在外翻力矩增加時應力比其他韌帶應力都要大，在抵抗外翻力矩作用時起主要作用，Ren 等［25］也表示內側副韌帶傾向于抵抗膝關節的外翻運動。同時也注意到，膝關節在過度外翻時，可能會導致內側副韌帶損傷。前交叉韌帶在內翻力矩達到最大時，應力比其他韌帶應力都要大，在內翻力矩作用時起主要作用，Hinckel等［26］的研究表明小于10°的內翻是導致前交叉韌帶應力增大的主要原因。股骨傾角結果與Gardiner 等［27］體外試驗在膝關節加載10 N·m 外翻力矩時外翻傾角為3.2°的結果相近，Gardiner 的結果也表明膝關節內翻時傾角略大于外翻時傾角。為膝關節在內外翻運動時的生物力學特性和運動損傷提供理論參考。

本文采用有限元方法分析膝關節復雜結構在承受不同載荷時的生物力學特性，與膝關節體外試驗及有限元結果進行對比分析。膝關節模型包括骨骼、軟骨、半月板和韌帶等復雜結構，基于有限元方法強大的通用性和有效性，可依據臨床特性精確調整材料屬性及載荷等參數，獲取膝關節組成中骨骼、軟骨、半月板和韌帶組織的應力分布及位移等，模型可為膝關節生理病變分析、進一步的手術方案設計和植入器械的力學特性研究提供有效可行的手段。