基于軸向功能梯度模型的人體脛骨運動力學特性分析

于 程， 李雅峰,2*， 孫 杰， 張 璽， 楊 威， 紀錚釗

(1.天津工業大學機械工程學院， 天津 300387； 2.天津工業大學現代機電裝備技術天津市重點實驗室， 天津 300387；3.天津市天津醫院創傷骨科， 天津 300211)

人體下肢屈曲過程中股骨內收，膝關節外展內扣的異常運動模式稱為動態膝外翻，其作為一種不正確的運動姿勢被認為是導致下肢損傷的危險因素之一，對日常生活中的下蹲、落地、跑步和減速等活動有直接的影響[1]。研究動態膝外翻時人體下肢的生物力學特性對預防運動損傷、輔助康復治療等有重要的理論和實際價值。

中外很多學者對動態膝外翻的發生原因和損傷機制做了廣泛的研究。Stickler等[2]對50例健康女性志愿者進行髖關節肌力與動態膝外翻程度的相關性分析，發現髖外展肌無力會導致股骨過度內收，進而造成動態膝外翻；Numata等[3]利用冠狀平面的二維運動學分析來評估291名運動員動態膝外翻角與非接觸性前交叉韌帶損傷之間的關系，經過三年隨訪發現存在動態膝外翻的運動員更容易出現前交叉韌帶損傷；Ford等[4]對81名運動員進行動態膝外翻運動測試，發現女性運動員動態膝外翻角度普遍大于男性運動員；Graci等[5]對20名髕-股疼痛患者進行動態膝外翻矯正訓練，對訓練前后的痛感進行對比評估發現疼痛的減輕與動態膝外翻程度變小有關；Bersini等[6]建立了相對完整的下肢多體動力學模型，通過仿真發現出現膝外翻時受力最大的韌帶為前交叉韌帶，同時得到了進行深蹲訓練時的脛-股接觸力；Tamura等[7]測量了單腿著地時下肢的垂直地面反作用力和角脈沖，發現動態膝外翻可能會降低著陸減速階段下肢的緩沖能力，從而增加下肢關節和骨骼所承受的沖擊。

人體膝關節的組成結構復雜，對脛骨作用力也受到多種因素的影響。本文中正常位姿和運動狀態下對于脛骨平臺施加的肌骨力和韌帶力數據，由相關文獻中的人體下肢運動實驗和膝關節骨骼力學實驗得到，并結合骨科醫院統計數據作為邊界載荷分布條件輸入模型進行分析。目前關于動態膝外翻的研究大多集中于前交叉韌帶損傷和性別差異性分析，對于運動時的人體主要承載骨骼——脛骨的受力狀況卻鮮有研究報道。本文旨在建立人體脛骨三維有限元模型，根據受力特點施加外載，以模擬正常和動態膝外翻位姿下的屈曲運動，分析脛骨的應力分布情況。

1 材料和方法

1.1 構建脛骨模型

在天津醫院骨科完成人體脛骨的CT(computed tomography)掃描采集。志愿者身體信息：男性，年齡29歲，體重65 kg，身高175 cm(被采集人對實驗知情同意)。實驗設備：GE Optima CT660 64排螺旋CT機。掃描條件設置：切片厚度0.625 mm，切片增量0.625 mm，圖像大小為512×512像素。共計787張圖像，圖像輸出格式為DICOM格式。

首先將DICOM格式的CT數據導入MIMICS中生成原始脛骨模型。將模型以二進制STL格式導入Geomagic里進行后續處理：在多邊形階段執行消除釘狀物、填充孔洞、光滑等操作；如圖1所示，在曲面處理階段繪制輪廓線、構建NURBS曲面片、構造柵格最終擬合得到脛骨精細模型。最后將脛骨精細模型以二進制STL格式導入3-matic中進行網格優化、體網格劃分等操作。采用十節點四面體單元生成脛骨有限元模型，有限元模型包含181 274個單元和326 216個節點。將模型導入ANSYS Workbench進行材料賦值以及設定邊界條件。

圖1 曲面處理階段流程

1.2 材料賦值和邊界條件

目前有多種關于骨骼材料的賦值方法。有研究將骨骼視為剛體，采用均一賦值，這種賦值方法材料屬性過于單一，不符合骨骼的實際情況[8-9]；大部分學者將骨骼分為密質骨和松質骨分別進行材料賦值，但是材料參數直接引自其他文獻，通用性有待考察[10-11]；還有學者采用灰度賦值，這種方法會在骨髓腔中生成大量的低彈性模量物質，對骨骼的應力分布有一定的影響[12]。基于上述不足，本文提出結構灰度賦值法對脛骨進行材料賦值。

脛骨在軸向上可分為一體兩端，如圖2所示，分別為兩端的小梁區和中間的骨干皮質區。盡管在骨小梁和皮質區之間沒有確切的界線，但是一些研究人員認為兩端的小梁區占脛骨總長度的20%[13-14]。本文中所構建的脛骨模型全長約為340 mm，因此三個區域的分界線設定為z1=68 mm和z2=282 mm。在兩端小梁區各均等提取10層斷層圖像，骨干區均等提取20層斷層圖像。利用MIMICS軟件，在每層斷層圖像上均等設置十條路徑進行CT值探測如圖3所示；圖4為圖3所示斷層圖像的CT值分布情況。

圖2 人體脛骨模型

圖3 CT圖像上探測路徑示意圖

圖4 CT值分布情況

提取所有數據，篩選去除髓腔中的數據點，利用均值法求得各區域的平均CT值。進而利用經驗公式[15](式1)得到各區域的表觀密度和楊氏模量。脛骨模型的材料參數如表1所示。

表1 脛骨模型各區域材料參數

(1)

式(1)中：D為表觀密度，g/cm3；Gv為灰度值(Gv=CT值+1 024)；E為楊氏模量，Pa。

由于骨骼的力學性能取決于其局部鈣化的程度，而鈣化程度在空間上是逐漸變化的[17]。因此兩個區域的材料屬性應該是逐漸變化的，而非斷崖式變化。為了更真實地模擬骨骼的特性，避免在區域分界處出現異常數據，如圖5所示，本文中在區域分界處截取長為l=10 mm的區域，將其材料屬性設置為功能梯度材料。骨干區和近端小梁區的體積分數在Z方向上按照下列關系分布[18]，即

圖5 脛骨近端功能梯度區域

(2)

式(2)中：Vc為骨干區的體積分數；Vt為近端小梁區的體積分數；h為功能梯度區中某一位置的高度；m為定義材料成分變化的常數，根據脛骨材料特性此處取值為0.5。

為了簡化計算，骨骼模型中的孔隙率忽略不計。功能梯度材料的楊氏模量和泊松比滿足關系式[19]

(3)

ν=νtVt+νcVc

(4)

式中：Eo為功能梯度區的楊氏模量；ν為功能梯度區的泊松比；Ec、νc分別為脛骨骨干區的楊氏模量和泊松比；Et、νt分別為脛骨近端小梁區的楊氏模量和泊松比。

如圖6所示，在脛骨平臺上分別施加980、1 360、1 820、2 360 N的垂直載荷以模擬膝關節屈曲30°、60°、90°、120°運動[20]。由于動態膝外翻導致下肢力線發生偏移，脛骨平臺外側受力增加，因此設定動態膝外翻屈曲模型脛骨平臺外側承受60%的外載[21-22]。脛骨近端施加4 N·m的內旋力矩Me，以模擬屈曲時脛骨的內旋運動[23]。在脛骨近端施加值為50 N、沿x軸負方向的外載F，以模擬動態膝外翻屈曲時脛骨在股骨的牽引下的內扣動作。將脛骨遠端25 mm長的區域設置為固定約束。

圖6 邊界條件示意圖

1.3 模型驗證

本研究依據Yamada[24]的實驗數據對所建脛骨模型進行三點彎曲驗證。對脛骨模型兩端施加約束，限制脛骨整體的轉動和位移。用直徑25 mm的剛性圓柱塊，以0.01 m/s的速度對脛骨模型中段施加外載。模型驗證分A-P(加載方向為前方指向后方)和L-M(加載方向為外側指向內側)兩個方向進行。輸出接觸載荷和變形量之間變化曲線與Yamada[24]實驗結果對比驗證。

2 結果

由圖7可見，非功能梯度模型的等效應力在區域分界處出現了明顯的跳躍，而功能梯度模型在區域分界處的過渡則相對平滑，沒有出現等效應力值突變的現象。

圖7 脛骨功能梯度模型與非功能梯度模型應力線圖對比

準靜態三點彎曲仿真結果如圖8所示，本研究仿真結果與Yamada[24]實驗結果在變化趨勢以及數值大小上都有很高的一致性，驗證了所建脛骨模型的有效性。

圖8 準靜態三點彎曲仿真結果與實驗結果對比

經有限元分析獲得的脛骨模型在不同狀態下的應力云圖如圖9所示。可以看出，動態膝外翻組脛骨模型的等效應力分布范圍明顯大于正常位姿組，整個脛骨內側面基本都有較大的應力分布。在4種屈曲角度下，正常位姿組的應力集中現象主要發生在脛骨中下段的正面，應力峰值大致在脛骨中下1/3交界處出現；動態膝外翻組的應力集中現象主要發生在脛骨中段內側，大致沿脛骨內側緣分布，由脛骨中段向兩端擴散并且逐漸降低。隨著屈曲角度從30°增大到120°，各組發生應力集中現象的位置并沒有明顯的變化。

圖9 兩種不同位姿屈曲30°～120°應力云圖

沿軸線方向在脛骨表面過應力峰值點進行應力值探測，繪制等效應力曲線如圖10所示。在達到應力峰值之前，正常位姿組的等效應力曲線近似為二次曲線。動態膝外翻組的等效應力變化近似為線性增長。就整體趨勢而言，高屈曲度(90°、120°)的等效應力變化速度明顯高于同組的低屈曲度(30°、60°)。

圖10 兩種不同位姿屈曲30°～120°等效應力曲線

對兩種位姿進行定量分析，圖11為脛骨最大等效應力對比。屈曲30°～120°的過程中，正常位姿組的最大等效應力分別為19.096、27.034、34.971、43.298 MPa，最大等效應力增幅分別為29%、22.69%、19.23%；動態膝外翻組的最大等效應力分別為31.203、44.802、53.762、59.74 MPa，最大等效應力增幅分別為30.35%、16.66%、10%。在4種屈曲角度下動態膝外翻組最大等效應力比正常位姿組分別增加了38.8%、39.65%、34.95%、27.52%。

圖11 兩種不同位姿屈曲30°～120°最大等效應力對比

提取脛骨平臺下方2 mm處橫截面的等效應力如圖12所示。可以清楚看出，正常位姿組脛骨內側的等效應力分布范圍明顯大于脛骨外側；而動態膝外翻組正好與之相反，脛骨外側的等效應力分布范圍明顯大于脛骨內側。

圖12 不同屈曲角度下脛骨平臺下方2 mm處的等效應力分布

如圖13所示，對比脛骨平臺下方2 mm處橫截面在不同屈曲角度下的最大等效應力可知，正常位姿組在四種屈曲角度下最大等效應力均為脛骨內側M大于脛骨外側L；而動態膝外翻組為脛骨外側大于脛骨內側。兩種位姿各屈曲角度下的最大等效應力在整體程度上相差不大。屈曲30°～120°的過程中，正常位姿組的最大等效應力分別為1.368 1、1.834 9、2.128 3、2.395 MPa；動態膝外翻組的最大等效應力分別為1.246 8、1.785 8、2.068 2、2.324 9 MPa。

圖13 兩種位姿屈曲脛骨平臺下方2 mm處最大等效應力對比

3 討論

準確預測等效應力在骨骼中的分布情況對預防運動損傷有重要意義。與傳統生物力學實驗手段相比，有限元仿真法能夠更為方便有效地預測人體骨骼在生理性和特殊受力狀態下的應力分布情況[25]。有限元仿真結果的準確性很大程度上取決于邊界條件的準確性，本文中較為全面地考慮了脛-股接觸力、肌肉力、韌帶力和重力等力的作用，保證了仿真的準確性。正常位姿屈曲時，脛骨內側承受相對較大的載荷，與文獻[26]中的實驗研究結果相似。Sobhani等[27]通過統計學分析發現髖內收膝外翻與脛骨內側損傷疼痛有正相關性，這與本文的研究結果相吻合。

圖7清楚地顯示了脛骨非功能梯度模型的等效應力在區域分界處有一個跳躍。由于區域分界處材料特性是不連續，應變場可能會出現跳躍現象，但是應力場應該是連續的(從靜力平衡開始)。因此，這些跳躍式錯誤的。跳躍兩邊的應力有很大的梯度，這些數據可能會誤導人們錯誤的預測等效應力的大小和位置。研究發現，劃分功能梯度區可以避免異常應力出現，使分界處的應力變化更為平滑。

本研究還發現在相同屈曲角度下，動態膝外翻屈曲所產生的最大等效應力以及應力分布范圍均大于正常位姿屈曲。且隨著屈曲角度的加大，最大等效應力也相應地升高。最大等效應力出現的位置也由正常位姿屈曲時的脛骨中下段正面轉移到脛骨中段內側。究其原因應該是髖內外展肌力矩不平衡使得脛骨承受了異常的側向力。相比于脛骨中下段，脛骨中段的橫截面積更小，因此動態膝外翻屈曲更容易對脛骨造成實質性損害。另外，由于下肢力線發生偏離，脛骨平臺內外側的受力狀況也隨之改變。正常位姿屈曲時，脛骨平臺內側的等效應力大于外側；而動態膝外翻屈曲時，脛骨外側承載更大地負荷，等效應力也更大，與腓骨對脛骨近端的支撐力相互作用，長此以往，會對脛骨外側平臺造成壓陷或剪切損傷。

4 結論

利用志愿者的下肢CT數據構建了較為精細的人體脛骨幾何模型。考慮到脛骨材料的軸向不均勻性，利用ANSYS等軟件中劃分了五個區域建立了脛骨軸向功能梯度有限元模型，分析比較了正常位姿屈曲和動態膝外翻屈曲脛骨模型上等效應力的分布情況，得到了如下結論。

(1)功能梯度模型能夠有效避免因材料屬性不連續而導致的異常應力。因此，對模型進行功能分級，使模型材料特性連續且平滑的變化，能夠得到更接近實際情況的有限元仿真結果。

(2)相較于正常位姿，動態膝外翻屈曲時，脛骨的應力分布情況發生改變，最大等效應力顯著增加，脛骨骨干中段以及脛骨平臺外側為易損區域。因此，日常生活及體育運動時應糾正位姿，避免出現動態膝外翻。同時加強對髖外展肌的鍛煉，必要時通過膝支具來避免損傷。