基于軸向功能梯度模型的人體脛骨運(yùn)動力學(xué)特性分析

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

(1.天津工業(yè)大學(xué)機(jī)械工程學(xué)院， 天津 300387； 2.天津工業(yè)大學(xué)現(xiàn)代機(jī)電裝備技術(shù)天津市重點(diǎn)實(shí)驗(yàn)室， 天津 300387；3.天津市天津醫(yī)院創(chuàng)傷骨科， 天津 300211)

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

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

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

1 材料和方法

1.1 構(gòu)建脛骨模型

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

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

圖1 曲面處理階段流程

1.2 材料賦值和邊界條件

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

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

圖2 人體脛骨模型

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

圖4 CT值分布情況

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

表1 脛骨模型各區(qū)域材料參數(shù)

(1)

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

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

圖5 脛骨近端功能梯度區(qū)域

(2)

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

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

(3)

ν=νtVt+νcVc

(4)

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

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

圖6 邊界條件示意圖

1.3 模型驗(yàn)證

本研究依據(jù)Yamada[24]的實(shí)驗(yàn)數(shù)據(jù)對所建脛骨模型進(jìn)行三點(diǎn)彎曲驗(yàn)證。對脛骨模型兩端施加約束，限制脛骨整體的轉(zhuǎn)動和位移。用直徑25 mm的剛性圓柱塊，以0.01 m/s的速度對脛骨模型中段施加外載。模型驗(yàn)證分A-P(加載方向?yàn)榍胺街赶蚝蠓?和L-M(加載方向?yàn)橥鈧?cè)指向內(nèi)側(cè))兩個方向進(jìn)行。輸出接觸載荷和變形量之間變化曲線與Yamada[24]實(shí)驗(yàn)結(jié)果對比驗(yàn)證。

2 結(jié)果

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

圖7 脛骨功能梯度模型與非功能梯度模型應(yīng)力線圖對比

準(zhǔn)靜態(tài)三點(diǎn)彎曲仿真結(jié)果如圖8所示，本研究仿真結(jié)果與Yamada[24]實(shí)驗(yàn)結(jié)果在變化趨勢以及數(shù)值大小上都有很高的一致性，驗(yàn)證了所建脛骨模型的有效性。

圖8 準(zhǔn)靜態(tài)三點(diǎn)彎曲仿真結(jié)果與實(shí)驗(yàn)結(jié)果對比

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

圖9 兩種不同位姿屈曲30°～120°應(yīng)力云圖

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

圖10 兩種不同位姿屈曲30°～120°等效應(yīng)力曲線

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

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

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

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

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

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

3 討論

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

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

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

4 結(jié)論

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

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

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