張 晶，郭 媛，張緒樹，張玉文

(太原理工大學(xué) 應(yīng)用力學(xué)與生物醫(yī)學(xué)工程研究所，太原 030024)

肩關(guān)節(jié)是人體結(jié)構(gòu)中最復(fù)雜、活動(dòng)度最大的關(guān)節(jié)之一[1]。肩部損傷常發(fā)生在各種碰撞中，無(wú)論是汽車行駛安全方面還是骨科領(lǐng)域均將其視為一種嚴(yán)重的關(guān)節(jié)損傷。由于肩部經(jīng)常是乘客和車廂或氣囊等約束系統(tǒng)作用時(shí)的一個(gè)主要載荷路徑，肩部的響應(yīng)會(huì)對(duì)身體其他部分造成影響，傳向胸部的載荷直接決定于肩部的響應(yīng)[2-4]。研究肩關(guān)節(jié)在側(cè)向碰撞工況下的響應(yīng)對(duì)預(yù)測(cè)碰撞對(duì)乘客的影響有重要意義。

近年來(lái)，采用有限元等數(shù)值方法來(lái)模擬肩關(guān)節(jié)的力學(xué)行為的研究受到普遍關(guān)注[3]。HAUG et al[5]構(gòu)建了完整的運(yùn)動(dòng)人體肩關(guān)節(jié)剛體有限元模型，該模型采用非線性彈簧模擬肩部關(guān)節(jié)的連接，可用來(lái)研究上肢受力和關(guān)節(jié)扭矩。INGRAM D[6]用有限元方法對(duì)肩部肌骨系統(tǒng)進(jìn)行了研究，分析了肩關(guān)節(jié)形態(tài)學(xué)結(jié)構(gòu)的功能，所建模型可再現(xiàn)肩部的動(dòng)態(tài)行為。DUPREY et al[7]根據(jù)COMPIGNE et al[8]尸體滑車實(shí)驗(yàn)建立模型，研究了肩部在不同沖擊力作用下的力學(xué)響應(yīng)并與實(shí)驗(yàn)進(jìn)行了對(duì)比。雖然這些模型和尸體實(shí)驗(yàn)進(jìn)行了驗(yàn)證，但存在大量簡(jiǎn)化，如均未考慮軟骨和肌肉的影響，同時(shí)模型網(wǎng)格較為粗糙，未能精確反應(yīng)肩關(guān)節(jié)幾何特征，影響仿真結(jié)果的準(zhǔn)確性。國(guó)內(nèi)楊曉霞[9]構(gòu)建了肩胛骨、鎖骨的三維有限元模型，研究了盂肱關(guān)節(jié)在4種外展角度工況下肩鎖關(guān)節(jié)的應(yīng)力分布情況。楊志勇等[10]建立了鎖骨的三維有限元模型并進(jìn)行軸向加載，發(fā)現(xiàn)鎖骨應(yīng)力主要集中在鎖骨干中約1/3區(qū)域，且下表面應(yīng)力集中較上表面明顯。禹銘楊[11]建立了不包括肌肉的肩部模型，研究了肩部形態(tài)學(xué)參數(shù)對(duì)肩峰撞擊癥的影響。國(guó)內(nèi)研究主要集中于用靜態(tài)加載來(lái)模擬肩部各種受力工況，且沒(méi)有關(guān)于肩部側(cè)向撞擊的報(bào)道。

本文旨在建立肩部的有限元模型，模擬肩部在側(cè)向沖擊作用下的力學(xué)響應(yīng)，為預(yù)測(cè)汽車碰撞或其他肩部側(cè)向沖擊情形下的肩部損傷提供參考。

1 材料和方法

1.1 幾何模型的建立

對(duì)一位肩部無(wú)畸形、病變成年男性志愿者的左側(cè)肩部進(jìn)行CT掃描，掃描層厚0.5 mm，得到CT斷層掃描圖像，以DICOM格式保存。

在Mimics中新建工程，導(dǎo)入包含左側(cè)肩部的CT掃描數(shù)據(jù)，設(shè)定分割閾值生成蒙板。對(duì)蒙板進(jìn)行裁剪、編輯、縮小范圍等操作，使鎖骨、肩胛骨、肱骨近端、肩部軟組織能明顯區(qū)分，在完成光滑、開閉等形態(tài)學(xué)操作之后進(jìn)行三維重建。對(duì)三維模型進(jìn)行光滑、三角形縮減、包裹等操作，以提高其質(zhì)量[12-14]。

根據(jù)韌帶解剖位置起止點(diǎn)，在骨模型表面選取相應(yīng)位置重建韌帶，考慮到韌帶對(duì)骨施加面載荷，使用一組彈簧單元模擬韌帶[14]。最終在Abaqus中重建出肩部主要韌帶：喙鎖韌帶、肩鎖韌帶、喙肩韌帶、盂肱韌帶[15]。

基于相鄰骨之間的空隙均由軟骨填充，且軟骨厚度以及空隙最小間距的一半作為假設(shè)，建立軟骨幾何模型[16]。將上述三維模型導(dǎo)入3-Matic中，根據(jù)軟骨解剖位置在相應(yīng)骨表面標(biāo)記相應(yīng)面，將標(biāo)記面從骨面上分離后，沿外法線方向偏移相應(yīng)距離生成軟骨，對(duì)實(shí)體進(jìn)行倒圓角、光滑等操作以提高其幾何質(zhì)量。

側(cè)向撞擊鋼板在Abaqus中創(chuàng)建為邊長(zhǎng)200 mm的正方形殼。

1.2 三維有限元模型的建立

將所建立的實(shí)體模型載入3-Matic中，生成體網(wǎng)格。因網(wǎng)格質(zhì)量對(duì)計(jì)算結(jié)果有重要影響，所以在生成體網(wǎng)格之前進(jìn)行了過(guò)濾鋒銳三角形、減少三角面片數(shù)量、網(wǎng)格質(zhì)量檢測(cè)、自動(dòng)面網(wǎng)格生成等操作，使網(wǎng)格質(zhì)量滿足計(jì)算要求[12]。肩部有限元模型參數(shù)見表1，完整的模型如圖1所示。

表1 肩部有限元模型參數(shù)Table 1 Parameters of the shoulder finite element model

圖1 肩部側(cè)向碰撞的有限元模型Fig.1 Finite element model of shoulder lateral impact

1.3 材料賦值

由于軟件和實(shí)驗(yàn)數(shù)據(jù)的局限性，本文涉及的生物材料屬性均假設(shè)為連續(xù)、均質(zhì)、線彈性。將劃分好網(wǎng)格的有限元模型以INP文件格式導(dǎo)入到Mimics中，基于CT圖像的灰度值，使用經(jīng)驗(yàn)公式以10種材料屬性對(duì)骨性部分進(jìn)行賦值[17]。根據(jù)參考文獻(xiàn)[18-19]，對(duì)鎖骨、肩胛依據(jù)經(jīng)驗(yàn)式(1)和(2)賦值：

D=1.067×VG-961.608 .

(1)

E=0.01×D1.86.

(2)

對(duì)肱骨依據(jù)經(jīng)驗(yàn)式(3)和(4)賦值：

D=0.624×VG-465.976 .

(3)

E=0.32×D1.41.

(4)

式中：D為材料密度，kg/m3；VG為CT圖像的灰度值；E為材料彈性模量，MPa.

肌肉、軟骨、鋼板的材料屬性見表2.韌帶設(shè)置為只能承受拉力的非線性彈簧單元，肩部主要韌帶參數(shù)見表3.

表2 材料屬性[20]Table 2 Mechanical properties of materials

1.4 接觸、載荷和邊界條件設(shè)置

軟骨、軟組織與所附著骨面之間的接觸均設(shè)為固定約束，以保證它們之間沒(méi)有相對(duì)位移。使用Abaqus中的面-面接觸算法模擬關(guān)節(jié)與軟骨之間的接觸，由于人體中關(guān)節(jié)面之間存在潤(rùn)滑，因此將關(guān)節(jié)面之間設(shè)為無(wú)摩擦的接觸關(guān)系。同樣地使用面-面接觸算法模擬肩部軟組織與鋼板的接觸，其中法向采用硬接觸，切向接觸設(shè)置為摩擦系數(shù)為0.6的庫(kù)倫摩擦[21]。

將鋼板置于yoz平面，并將其6個(gè)自由度完全約束。調(diào)整肩與鋼板的相對(duì)位置，使肩部軟組織外側(cè)最大程度接近鋼板，以縮短計(jì)算時(shí)間。肩部施加沿x軸方向的速度場(chǎng)，大小為3.5 m/s，分析時(shí)間設(shè)置為0.02 s.

根據(jù)運(yùn)動(dòng)生物力學(xué)知識(shí)可知，肩胛骨活動(dòng)受到肩胛胸壁關(guān)節(jié)的限制，使肩胛骨在運(yùn)動(dòng)中和胸壁緊密切合，防止了翼狀飄起，近似可認(rèn)為其沿肩胛胸壁方向無(wú)線位移，故本文限制了肩胛骨前骨面上部分節(jié)點(diǎn)沿y軸方向的線位移。鎖骨由胸鎖關(guān)節(jié)固定于第1肋骨上，鎖骨近端受到胸鎖關(guān)節(jié)的限制使其旋轉(zhuǎn)運(yùn)動(dòng)較多而線位移較少，模擬中限制鎖骨近端z軸方向的線位移來(lái)代表胸鎖關(guān)節(jié)對(duì)其的約束[9,15,22]。整個(gè)動(dòng)態(tài)碰撞求解過(guò)程在Abaqus顯示求解器中完成，耗時(shí)29 h.求解完成之后進(jìn)入Abaqus可視化結(jié)果處理頁(yè)面，對(duì)結(jié)果進(jìn)行分析。

模型驗(yàn)證參照楊曉霞[9]肩鎖關(guān)節(jié)力學(xué)測(cè)試實(shí)驗(yàn)的固定和加載方式(如圖2)，將肩胛骨與胸壁接觸一側(cè)的部分節(jié)點(diǎn)進(jìn)行固定約束，分別用70,80,90,100 N的拉力沿x軸負(fù)向作用于鎖骨近端，使用Abaqus求解器中的隱式方法完成模擬。選取如下測(cè)量點(diǎn)：點(diǎn)1為鎖骨外端中點(diǎn)，點(diǎn)2為鎖骨外端中點(diǎn)在肩峰上對(duì)應(yīng)點(diǎn)(內(nèi)外方向)，點(diǎn)3為喙突尖頂點(diǎn)。測(cè)量標(biāo)志點(diǎn)的位移，將標(biāo)志點(diǎn)位移曲線與實(shí)驗(yàn)結(jié)果對(duì)比。

圖2 沿x軸橫向拉力加載的實(shí)驗(yàn)和模擬Fig.2 Experiment and simulation of x axial tension

2 實(shí)驗(yàn)結(jié)果

2.1 靜態(tài)加載位移情況

施加橫向拉力下目標(biāo)點(diǎn)的位移結(jié)果見表4.從模擬和實(shí)驗(yàn)的對(duì)比圖(圖3)可見，隨著拉力的增加，目標(biāo)點(diǎn)位移也在增加，喙突尖頂點(diǎn)位移最大，而鎖骨外端中點(diǎn)在肩峰上對(duì)應(yīng)點(diǎn)位移較小，位移整體趨勢(shì)與實(shí)驗(yàn)結(jié)果一致。因本模型材料屬性均設(shè)定為線彈性，故位移結(jié)果表現(xiàn)出一定的線性特征。

表4 沿x軸橫向拉力加載方式下目標(biāo)點(diǎn)的位移Table 4 Displacement of target point under x axial tension load

圖3 沿x軸拉力加載下目標(biāo)點(diǎn)位移Fig.3 Displacement of target point under x axial tension force

2.2 側(cè)向沖擊后肩部整體受力情況

本文研究的碰撞過(guò)程歷時(shí)20 ms：0～12 ms時(shí)肩部與墻面發(fā)生接觸；13～20 ms時(shí)肩部進(jìn)入無(wú)外力作用，即變形恢復(fù)階段。模擬展現(xiàn)了碰撞的整個(gè)過(guò)程。

碰撞中von Mises應(yīng)力最大值發(fā)生在11 ms時(shí)，大小為70.12 MPa，位置在肩胛頸內(nèi)側(cè)，這與臨床觀察肩胛頸易發(fā)生骨折位置相符[23]。整個(gè)碰撞過(guò)程中最大von Mises應(yīng)力變化如圖4所示。結(jié)合整個(gè)碰撞的動(dòng)態(tài)過(guò)程可知，最大von Mises應(yīng)力起初出現(xiàn)于肩部肌肉，肩部肌肉的變形引起肱骨的受力變化，最大應(yīng)力由肱骨向肩胛骨轉(zhuǎn)移；3～11 ms最大應(yīng)力主要發(fā)生在肩胛骨上，這反映出在碰撞中肩胛骨起著主要承擔(dān)載荷的作用。碰撞作用的后期，肩部反彈與剛性面脫離接觸，肩部?jī)?chǔ)存的能量釋放，肩部各組成部分之間相互作用，最大應(yīng)力發(fā)生在肱骨和鎖骨遠(yuǎn)端軟骨上。

圖4 最大von Mises應(yīng)力-時(shí)間曲線Fig.4 Max von Mises stress versus time curve

從圖5最大法向接觸力-時(shí)間曲線可看出，肩部與剛性墻接觸后，法向接觸力開始上升，此時(shí)肩部速度在下降，肩部肌肉變形壓縮，在3 ms時(shí)法向接觸力達(dá)到峰值12.77 kN，這表示軟組織被最大程度壓縮，肩部各組織大部分區(qū)域速度降為0.隨后肩部反彈，法向接觸力開始下降；在7 ms時(shí)法向接觸力下降到低谷后又開始上升，從碰撞運(yùn)動(dòng)過(guò)程推測(cè)這可能是因?yàn)殡殴呛图∪庋豿軸正向恢復(fù)變形速度大于肩部沿x軸負(fù)向反彈速度，肩部和剛性墻再次發(fā)生了接觸壓縮，從而使接觸力上升。在13 ms接觸力下降為0，意味著肩部肌肉最下層和鋼板脫離，沖擊結(jié)束。

圖5 最大法向接觸力-時(shí)間曲線Fig.5 Max normal contact force versus time curve

2.3 側(cè)向沖擊后肩部各骨骼受力情況

2.3.1 肱骨骨應(yīng)力分布

整個(gè)碰撞中肱骨所受最大等效應(yīng)力發(fā)生在15 ms，大小為47.15 MPa.從圖6肱骨截面圖可見最大應(yīng)力位于內(nèi)表面，并且在整個(gè)碰撞中肱骨應(yīng)力集中主要發(fā)生在其內(nèi)部。碰撞中應(yīng)力逐漸由肱骨干移向肱骨近端，隨著力的傳遞，應(yīng)力集中在頸干交界部位和肱骨干上部，在其內(nèi)側(cè)產(chǎn)生聚積。同時(shí)，與關(guān)節(jié)盂相接觸的肱骨關(guān)節(jié)面部分，應(yīng)力也在逐漸增加，但肱骨頭上的應(yīng)力明顯大于肱骨干，據(jù)此推測(cè)肱骨頭和肱骨頸干交界部位更易發(fā)生骨折。

圖6 最大應(yīng)力發(fā)生時(shí)肱骨應(yīng)力云圖Fig.6 Von Mises stress of humerus at the time of maximum stress

2.3.2 鎖骨應(yīng)力分布

鎖骨上的最大等效應(yīng)力發(fā)生在14 ms，大小為13.96 MPa，整個(gè)碰撞過(guò)程中鎖骨上應(yīng)力分布較為均勻。從14 ms時(shí)的鎖骨應(yīng)力云圖(圖7)可見鎖骨下表面應(yīng)力集中更明顯，最大應(yīng)力發(fā)生在鎖骨干中部約1/3處，且近端較遠(yuǎn)端更易發(fā)生應(yīng)力集中。同時(shí)觀察應(yīng)力云圖發(fā)現(xiàn)，下表面比上表面更易發(fā)生應(yīng)力集中，說(shuō)明撞擊中鎖骨下表面斷裂風(fēng)險(xiǎn)更高。

圖7 最大應(yīng)力發(fā)生時(shí)鎖骨應(yīng)力云圖Fig.7 Von Mises stress of clavicle at the time of maximum stress

2.3.3 肩胛骨的應(yīng)力分布

整個(gè)碰撞過(guò)程中最大應(yīng)力發(fā)生在肩胛骨上，在最大應(yīng)力發(fā)生時(shí)刻，肩胛骨上的應(yīng)力分布如圖8所示。肩胛頸處有明顯的應(yīng)力集中，且整個(gè)過(guò)程中肩胛頸相對(duì)于周圍的岡下窩、肩胛岡均有較大的應(yīng)力值和更明顯的應(yīng)力集中現(xiàn)象；表明力經(jīng)肩鎖和盂肱關(guān)節(jié)傳導(dǎo)時(shí)，肩胛頸易發(fā)生骨折，與WILLIAMS et al[24]統(tǒng)計(jì)的大多數(shù)肩胛骨骨折為外科頸骨骨折的結(jié)論一致。因此在肩部骨折固定時(shí)，應(yīng)首先考慮肩胛頸的完整性。

圖8 最大應(yīng)力發(fā)生時(shí)肩胛骨應(yīng)力云圖Fig.8 Von Mises stress of scapula at the time of maximum stress

3 討論

近年來(lái)，有限元分析法在生物力學(xué)研究中的應(yīng)用日益廣泛和深入，并且在處理復(fù)雜生物力學(xué)問(wèn)題時(shí)它的優(yōu)勢(shì)得到體現(xiàn)，已成為分析和預(yù)測(cè)骨應(yīng)力、應(yīng)變的重要工具。采用有限元分析法，肩關(guān)節(jié)的結(jié)構(gòu)形狀、組織性能、載荷邊界條件等均可用數(shù)學(xué)形式概括，可通過(guò)改變其中任一參數(shù)以觀察其對(duì)整體結(jié)構(gòu)的影響，從而解釋損傷中的力學(xué)變化[1,6,20]。現(xiàn)在人體生物力學(xué)有限元模型的建立大多依靠CT等影像學(xué)圖像和Abaqus等有限元軟件[1]，這樣一方面可以增加模型模擬的準(zhǔn)確性，使模型更趨于真實(shí)；另一方面簡(jiǎn)化了建模過(guò)程，節(jié)省人力物力。

本文基于逆向工程原理建立了人體左肩部復(fù)合體，在Mimics中將CT掃描所得連續(xù)斷層圖像進(jìn)行重建、分割、優(yōu)化處理，建模過(guò)程中始終保證軟件對(duì)數(shù)據(jù)的自動(dòng)識(shí)別優(yōu)化。模型采用了骨廣泛的力學(xué)參數(shù)，不再將骨簡(jiǎn)單區(qū)分為松質(zhì)骨和密質(zhì)骨。根據(jù)張國(guó)棟等[17]的研究，將骨有限元模型賦予10種材料屬性后可準(zhǔn)確反映其真實(shí)應(yīng)力狀態(tài)，因此本文骨的材質(zhì)同樣設(shè)為10種。本文所建模型包括了從CT圖像中提取的軟組織，更符合實(shí)際的側(cè)向撞擊情況，體現(xiàn)了軟組織對(duì)沖擊能量的吸收作用。肩部實(shí)體模型均劃分為4節(jié)點(diǎn)四面體單元，韌帶用非線性彈簧單元代替，整個(gè)模型共135 677個(gè)節(jié)點(diǎn)，516 392個(gè)單元，高度精細(xì)劃分。模型的建立較好地體現(xiàn)了肩部的復(fù)雜幾何形狀和骨密度的不均勻性，保證了模型從數(shù)據(jù)源上的準(zhǔn)確性和有效性。

靜態(tài)加載模擬所得結(jié)果與楊曉霞[9]所測(cè)結(jié)果吻合，證明所建模型合理有效。動(dòng)態(tài)碰撞所得結(jié)果表明碰撞中肩胛頸骨折風(fēng)險(xiǎn)最大，肱骨頭、肱骨頸、鎖骨下表面均為損傷易發(fā)區(qū)，由此為肩部碰撞救治損傷位置的確定提供了指導(dǎo)。同時(shí)最大等效應(yīng)力-時(shí)間曲線和接觸力-時(shí)間曲線也可為肩關(guān)節(jié)護(hù)具的開發(fā)提供參考，護(hù)具設(shè)計(jì)可圍繞降低Von Mises峰值應(yīng)力和接觸力來(lái)進(jìn)行設(shè)計(jì)。

需要指出的是，因缺少肩部側(cè)向撞擊的實(shí)驗(yàn)數(shù)據(jù)，所以所建模型只能通過(guò)對(duì)比部分結(jié)構(gòu)靜態(tài)加載位移進(jìn)行驗(yàn)證，這顯然未能完整顯示出模型與真實(shí)肩部的接近程度，今后需要進(jìn)行力、位移以及其他參數(shù)的綜合比對(duì)。側(cè)向撞擊速度選擇COMPIGNE et al[8]統(tǒng)計(jì)的肩部側(cè)向撞擊平均速度3.5 m/s，設(shè)置為純側(cè)向碰撞，實(shí)際碰撞中可能有更高的速度且方向多變，今后需要進(jìn)行更高速度不同方向碰撞結(jié)果的對(duì)比。本文的結(jié)果對(duì)于了解肩部側(cè)向撞擊過(guò)程中各組織的生物力學(xué)響應(yīng)具有一定的學(xué)術(shù)價(jià)值，對(duì)于肩部撞擊的預(yù)防保護(hù)以及撞擊后的臨床治療具有一定的參考價(jià)值。

4 結(jié)束語(yǔ)

本文采用組織逆向重建方法，利用Mimics、3-Matic、Abaqus等軟件構(gòu)建了人體左肩部各組織三維有限元模型，并驗(yàn)證了模型的有效性。純側(cè)向撞擊結(jié)果表明，肩部骨折固定時(shí)，應(yīng)首先考慮肩胛頸的完整性，肱骨頭和肱骨頸干交界為肱骨骨折高風(fēng)險(xiǎn)區(qū)，鎖骨干中部約1/3處應(yīng)力集中解釋了該處易發(fā)骨折的原因。今后還需完善模型，嘗試以實(shí)體單元來(lái)構(gòu)建韌帶和使用超彈性材料屬性，進(jìn)一步研究不同速度側(cè)向碰撞對(duì)肩部以及上軀干的影響。