宋光輝，王廣亮，李大鵬，韓曉強，李志剛

（1.北京交通大學(xué)機械與電子控制工程學(xué)院，北京 100044； 2.首都醫(yī)科大學(xué)附屬北京兒童醫(yī)院，北京 100045）

前言

我國平均每年交通事故50萬起，所導(dǎo)致的死亡人數(shù)超過10萬人。統(tǒng)計數(shù)據(jù)表明，每5 min就會有人因交通事故而死亡，每年因為交通事故造成的經(jīng)濟損失達到數(shù)億元。在2011年，我國由于交通事故導(dǎo)致的死亡人數(shù)為6.24萬人，大約占到世界交通死亡人數(shù)的16%［1］。汽車碰撞中頸部尤其是韌帶具有較高的損傷風(fēng)險，統(tǒng)計數(shù)據(jù)顯示頸部的損傷概率在交通事故中呈現(xiàn)上升趨勢，且交通事故中頸部的損傷引起后遺癥的概率高達50%［2］。因此，研究典型載荷形式下頸部不同韌帶響應(yīng)形式和特點具有重要意義。

尸體試驗是生物力學(xué)領(lǐng)域研究的黃金法則。目前已有較多文獻通過尸體試驗研究了成人頸椎的傷害機理和承載情況。Camacho［3］利用 10個成人尸體頸部樣本進行了準靜態(tài)的彎曲試驗，得到C7-T1/C67/C56/C45/C34/C23/O-C2角位移響應(yīng)與不同彎矩之間的定量關(guān)系和標(biāo)準差，同時擬合了角位移與彎矩的關(guān)系曲線；Nightingale等人［4］利用16個51-74歲的未防腐處理的成人男性頸椎，共計41個O-C2/C45/C67頸椎段樣本，進行了前屈后伸的純彎曲試驗，獲得了前屈后伸工況下O-C2/C45/C67的角位移響應(yīng)和標(biāo)準差，并擬合了其角位移響應(yīng)與彎矩的非線性關(guān)系；指出在前屈后伸工況下，男性上頸椎段的運動幅度明顯大于下頸椎的運動幅度，男性上頸椎段的剛度明顯大于女性上頸椎段的剛度；Wheeldon等人［5］利用成人 C23/C34/C45/C56/C67頸椎段樣本進行了準靜態(tài)的彎曲試驗，得出成人各頸椎段在不同彎矩作用下的角位移響應(yīng)的平均值和標(biāo)準差；Luck［6］對約 5個月到18歲的樣本進行了非破壞準靜態(tài)的彎曲試驗和動態(tài)拉伸破壞試驗等，并利用試驗數(shù)據(jù)研究了頸部的彎曲響應(yīng)與彎矩的定量關(guān)系，研究了嬰幼兒到成年全年齡段不同頸椎段拉伸失效力和失效位移與年齡的關(guān)系等。這些有限的試驗數(shù)據(jù)為理解頸椎關(guān)節(jié)的損傷形式和機理提供了良好的基礎(chǔ)，同時為有限元模型的驗證提供了寶貴的基礎(chǔ)數(shù)據(jù)。

進行尸體試驗的難度較大，隨著技術(shù)的發(fā)展，有限元分析已應(yīng)用到頸椎的損傷機理研究中。Hakin和King等人［7］建立了一塊單獨椎骨的有限元模型；Yang和King［8］又將這塊椎骨模型加以完善，建立了一段包含椎間盤關(guān)節(jié)面的完整脊椎運動節(jié)段模型；Goal等人［9］第一次用梁單元模擬韌帶，研究了頸椎在前屈、后仰和軸向轉(zhuǎn)動時的響應(yīng)；Panzer等人［10］建立了成人C45頸椎段有限元模型，研究了在前屈、后伸、軸向轉(zhuǎn)動、側(cè)向彎曲4種工況下椎間盤、韌帶、椎骨關(guān)節(jié)的載荷分配情況，得到的結(jié)論是，在軸向轉(zhuǎn)動和側(cè)向彎曲工況下椎間盤為主要的承載結(jié)構(gòu)，在后伸工況下椎骨關(guān)節(jié)為主要的承載結(jié)構(gòu)，而在前屈工況下韌帶為主要承載結(jié)構(gòu)。Czyz等人［11］建立了人體頸髓有限元模型，并利用豬頸髓試驗驗證了模型的有效性。付裕等人［12］建立了椎間盤有限元模型，研究纖維環(huán)的應(yīng)力分布。髓核和終板軟骨采用四面體網(wǎng)格，髓核采用簡單的不可壓縮材料，得到在壓縮力的作用下，由于纖維環(huán)向前膨脹產(chǎn)生的拉應(yīng)力在椎體的前方最大，屈曲時應(yīng)力主要集中在后部，扭轉(zhuǎn)時最大應(yīng)力出現(xiàn)在椎間盤后側(cè)的結(jié)論；王芳等人［13］建立了中國人全頸椎有限元模型用于揮鞭傷的損傷分析，其韌帶采用線彈性材料，椎間盤不區(qū)分髓核、纖維環(huán)和纖維基質(zhì)等結(jié)構(gòu)。曹立波等人［14］建立了首個中國50百分位人體的頸椎有限元模型，并利用頸椎離體墜落試驗對模型的有效性進行驗證，該模型松質(zhì)骨和皮質(zhì)骨都采用各向同性材料，韌帶和纖維環(huán)采用線彈性材料。

本文中首先建立了C67頸椎段高質(zhì)量網(wǎng)格有限元模型；然后，對模型進行了前屈、后伸、軸向轉(zhuǎn)動和側(cè)向彎曲等多個工況的尸體試驗驗證；最后，分析在不同載荷形式下不同韌帶的響應(yīng)和損傷特點。

1 C67頸椎段有限元模型的建立

頸椎的樣本從邯鄲市中心醫(yī)院獲得，其中獲取CT的完整流程得到了倫理委員會的審批，CT數(shù)據(jù)顯示頸椎無損傷和畸變，劃分成人C67頸椎段有限元網(wǎng)格模型的流程如圖1所示。

圖1 成人C67頸椎段有限元網(wǎng)格流程圖

1.1 椎骨網(wǎng)格劃分

1.1.1 椎骨幾何重建和幾何區(qū)域劃分

將頸椎CT掃描獲得的DICOM格式文件導(dǎo)入Mimics軟件中，進行斷層圖像分割和三維重建，獲得三維模型后，將模型輸出為 STL格式并導(dǎo)入Geomagic Studio，對頸椎段表面進行松弛和光滑處理，最終得到椎骨的三維模型。

當(dāng)模型的幾何形狀較復(fù)雜時，常規(guī)的網(wǎng)格劃分軟件不易劃分出高質(zhì)量規(guī)則六面體單元。本文中使用結(jié)構(gòu)化網(wǎng)格劃分軟件Truegrid（XYZ）進行網(wǎng)格劃分，其基本思想是將復(fù)雜的三維幾何模型分割成若干相對規(guī)則的子區(qū)域，逐塊進行網(wǎng)格劃分，然后進行網(wǎng)格對接。對于本文的C67頸椎段，由于橫突孔周圍的骨骼較細且有更大的曲率變化，為避免相對的投影幾何面夾角過大、截面積突變，在分段時必須對其進行細節(jié)分割處理。根據(jù)各個部位的結(jié)構(gòu)特點，將成人C6分割為24個區(qū)域，C7分割為16個區(qū)域。椎骨的三維模型和區(qū)域劃分結(jié)果如圖2所示。

圖2 椎骨的三維模型和區(qū)域劃分結(jié)果

1.1.2 椎骨有限元模型的建立

以C7中一段為例說明網(wǎng)格劃分的流程，如圖3所示。首先，導(dǎo)入被劃分區(qū)域的幾何文件，并根據(jù)分塊的大小設(shè)定block的網(wǎng)格數(shù)量；然后，將建立的block體生成蝶形網(wǎng)格，并通過移動和旋轉(zhuǎn)操作將建立的蝶形網(wǎng)格定位到幾何模型中對于網(wǎng)格投影最有利的位置；其次，分別將頂點、邊界和面投影到對應(yīng)幾何的點線面上；最后，對網(wǎng)格進行調(diào)整、順滑。

圖3 Trigurid塊投影流程圖

相鄰的兩塊投影完后需要通過網(wǎng)格對接將相鄰區(qū)域同一個截面上的網(wǎng)格節(jié)點合并，圖4為C7棘突處相鄰兩個區(qū)域合并節(jié)點前后的對比圖。

圖4 網(wǎng)格對接對比

重復(fù)上述步驟依次完成24個C6區(qū)域和16個C7區(qū)域的網(wǎng)格投影，最終得到高質(zhì)量C67的椎骨網(wǎng)格（圖5）。在其表面提取一層殼單元，得到皮質(zhì)骨和終板的殼單元網(wǎng)格。

圖5 C67椎骨有限元網(wǎng)格

1.1.3 椎間盤有限元模型的建立

椎間盤包括纖維環(huán)、髓核和終板軟骨。其中纖維環(huán)由纖維環(huán)基質(zhì)和纖維環(huán)纖維薄膜組成。參照圖譜并按照椎體上下幾何面，建立貼合椎體上下曲面的終板軟骨曲面，按照約50%的截面積對曲面進行分割得出髓核區(qū)域［15］。終板軟骨和纖維環(huán)纖維采用殼單元，髓核和纖維基質(zhì)采用六面體單元。將纖維環(huán)基質(zhì)的實體網(wǎng)格提取出殼單元，獲得各基質(zhì)層側(cè)邊的殼網(wǎng)格。纖維環(huán)纖維薄膜共有4層，根據(jù)殼網(wǎng)格從內(nèi)向外的順序定義為第1至第4層。椎間盤有限元模型如圖6所示。椎間盤與椎骨之間通過Contact_tied的方式進行連接。

圖6 椎間盤模型

1.1.4 頸椎關(guān)節(jié)有限元模型的建立

頸椎關(guān)節(jié)由關(guān)節(jié)軟骨、關(guān)節(jié)囊韌帶及其所包含的關(guān)節(jié)滑液組成，頸椎關(guān)節(jié)可用于緩沖兩關(guān)節(jié)面的接觸沖擊。本文中對關(guān)節(jié)軟骨建立實體單元，使用Contact_tied方式連接到關(guān)節(jié)突上，定義關(guān)節(jié)軟骨的面面接觸，關(guān)節(jié)囊韌帶簡化為梁單元。關(guān)節(jié)軟骨的厚度依據(jù)Yoganandan等［16］提出的成人關(guān)節(jié)軟骨厚度。關(guān)節(jié)面模型如圖7所示。

圖7 關(guān)節(jié)面網(wǎng)格模型

1.1.5 韌帶有限元模型的建立

頸椎韌帶較多且在CT掃描圖像中難以觀察重建。因此，基于解剖學(xué)圖譜中韌帶的位置關(guān)系建立韌帶結(jié)構(gòu)。韌帶的建立一般有兩種方式：一種是按其較為真實的幾何結(jié)構(gòu)建立韌帶；另一種是使用等效的離散梁單元建立韌帶。前者能較為準確地反映韌帶幾何形狀，但由于韌帶的不規(guī)則結(jié)構(gòu)，故其對材料屬性要求較高，需要對應(yīng)部分的應(yīng)力-應(yīng)變試驗數(shù)據(jù)。而目前已有的數(shù)據(jù)多數(shù)是韌帶宏觀的力-位移曲線。考慮到上述限制，本模型中韌帶使用一維離散梁單元進行建模。韌帶模型如圖8所示。

圖8 韌帶模型

為保證模型與已有試驗對象的一致性，成人C67頸椎段不建立肌肉模型。模型建立完成后，進行網(wǎng)格質(zhì)量檢查，發(fā)現(xiàn)成人C67模型最小雅各比為0.31，具體的網(wǎng)格數(shù)量和質(zhì)量信息如表1所示。從表1看出，整個模型具有較高的網(wǎng)格質(zhì)量，能保證后期模型的計算精度與效率。

表1 網(wǎng)格質(zhì)量信息表

1.2 材料屬性

1.2.1 松質(zhì)骨材料屬性

目前，頸部有限元模型中關(guān)于松質(zhì)骨材料的表征尚未有統(tǒng)一的模式，部分文獻認為松質(zhì)骨的材料屬性使用各向同性材料［17-19］，另一部分文獻認為其屬于各向異性材料［20-21］。表2給出了文獻中的兩種松質(zhì)骨材料模型。

為獲取頸椎松質(zhì)骨相關(guān)力學(xué)特點，使用豬的松質(zhì)骨作為等效樣本進行兩個方向下的壓縮試驗，利用圓鉆將不同部位下的椎骨按照軸向、徑向取出，試驗樣本的獲取過程如圖9所示。

表2 松質(zhì)骨材料模型

對軸向和徑向的試驗樣件進行了準靜態(tài)壓縮試驗，試驗照片如圖10所示。結(jié)果表明，松質(zhì)骨在軸向和徑向方向材料屬性并不相同，呈現(xiàn)出明顯的各向異性的特征，因此模型采用各向異性材料（*MAT_OPTION TROPIC_ELASTIC）進行建模。試驗測試結(jié)果發(fā)現(xiàn)各向異性程度與文獻［17］所介紹的成人松質(zhì)骨兩個方向較為接近，因此選取該文獻的松質(zhì)骨材料數(shù)據(jù)作為最終使用的材料屬性。

1.2.2 韌帶材料屬性

圖9 椎骨原始材料獲取過程

韌帶使用離散梁單元進行模擬，該模型可通過直接輸入力-位移曲線用于表征韌帶在拉伸過程中拉伸力的非線性變化特點。

圖10 樣件準靜態(tài)壓縮試驗

試驗結(jié)果表明，韌帶的力隨位移變化的曲線可近似視為“S”形，包括3個區(qū)域。將韌帶的拉伸力隨應(yīng)變的變化曲線進行歸一化，如圖11所示。需要3個控制點，第1個階段是0點到A點的非線性區(qū)域，第2個階段從A到B的線性區(qū)域，第3個區(qū)域是B到C的非線性區(qū)域。通過A，B，C 3點的數(shù)據(jù)，進行近似回歸，即可得到頸椎所有韌帶的力-位移曲線［22］。獲得的所有成人韌帶的A，B和C點的參數(shù)如表3所示。

圖11 韌帶力-應(yīng)變歸一化后的典型曲線

1.2.3 其他部位材料屬性

皮質(zhì)骨、終板、終板軟骨、纖維基質(zhì)、纖維環(huán)纖維、髓核和小關(guān)節(jié)軟骨的材料本構(gòu)與單元類型如表4所示。

2 模型驗證

2.1 成人C67準靜態(tài)前屈-后伸驗證

按照Camacho［3］提出的試驗工況進行仿真來驗證模型的逼真度。在前屈后伸仿真中，約束C7的所有自由度，然后對 C6分別施加±0.5，±1.0，±1.5，±2.0 N·m的彎矩，得到C6的前屈、后伸角位移響應(yīng)，如圖12和圖13所示。

表3 成人相關(guān)韌帶的A，B和C的數(shù)值

可以看出：對于前屈工況，C6前屈角位移響應(yīng)落在 Camacho和 Nightingale［4］的試驗區(qū)間內(nèi)，相對于Wheeldon［5］的試驗，有限元模型的前屈角位移響應(yīng)普遍偏小，存在這種誤差的原因可能是試驗加載位置測量不準確，仿真加載位置與試驗加載的真實位置不一致［10］，也可能是椎骨的中性區(qū)與模型的中性區(qū)不一致［27］，樣本之間的個體差異也可能是誤差的來源；對于后伸工況，成人C6角位移響應(yīng)均落在Camacho，Nightingale和Wheeldon的試驗區(qū)間內(nèi)；總體來說，成人C67頸椎段有限元模型能很好地重現(xiàn)成人C67頸椎段的前屈后伸工況。

2.2 成人C67準靜態(tài)軸向轉(zhuǎn)動和側(cè)向彎曲驗證

Panjibi［28］對 C67頸椎段進行了 1.0 N·m力矩作用下的軸向轉(zhuǎn)動和側(cè)向彎曲試驗，運動示意圖如圖14（a）所示。按照Panjibi的試驗方式，對有限元模型進行仿真，仿真后得到的C6頸椎段的角位移響應(yīng)如圖14（b）所示。模型軸向轉(zhuǎn)動和側(cè)向彎曲的角位移落在了Panjibi的試驗區(qū)間內(nèi)，模型側(cè)向彎曲剛度比Panjibi試驗的抗側(cè)彎剛度稍大?？偟膩碚f，該模型能很好地模擬成人頸部C67頸椎段的側(cè)向彎曲和軸向轉(zhuǎn)動工況。

表4 頸椎其他部位有限元模型的材料參數(shù)

圖12 C6前屈角位移響應(yīng)

圖13 C6后伸角位移響應(yīng)

本文中對成人C67頸椎段進行了前屈、后伸、軸向轉(zhuǎn)動和側(cè)向彎曲等4種工況下的準靜態(tài)驗證。C6角位移響應(yīng)都落在試驗區(qū)間內(nèi)，僅側(cè)向彎曲接近Panjabi試驗的區(qū)間下界，說明在側(cè)彎工況下，模型的剛度略高，總體來說成人C67頸椎段有限元模型能很好地反映成人C67頸椎段的運動響應(yīng)，具有較高的生物力學(xué)逼真度。

圖14 C67頸椎段準靜態(tài)工況與仿真結(jié)果

3 成人C67頸椎段韌帶對載荷形式的響應(yīng)特點

為研究成人C67頸椎段韌帶組織的伸長率響應(yīng)的特點，將C7固定，對C6施加大小為1.0 N·m的前屈、后伸、軸向轉(zhuǎn)動和側(cè)向彎曲等4種載荷。5種韌帶在4種載荷下的伸長率如圖15所示。

圖15 韌帶響應(yīng)

從圖15可知：前屈工況時，前縱韌帶（ALL）、后縱韌帶（PLL）、黃韌帶（LF）、關(guān)節(jié)囊韌帶（CL）和棘間韌帶（ISL）伸長率分別為 0，1.7%，6.4%，13.05%和21.4%，ISL伸長率最大，CL次之；后伸工況時，ALL和CL伸長率分別為5.16%和0.85%，其他韌帶的伸長率幾乎為零，ALL的伸長率最大，且明顯高于其他韌帶的伸長率；軸向轉(zhuǎn)動時，ALL，PLL，LF，CL和ISL伸長率分別為0.59%，0，0.45%，13.47%和1.09%，CL的伸長率最大，為韌帶中主要的承載對象；側(cè)向彎曲時，CL的伸長率最大為12.87%，ALL，PLL，LF和ISL 4種韌帶的伸長率很小甚至接近于0。

從上述結(jié)果還可看出：在1.0 N·m的前屈、后伸、軸向轉(zhuǎn)動和側(cè)彎等4種載荷作用下，ALL在后伸工況下的伸長率最大為5.16%，明顯高于其他3種工況下伸長率的響應(yīng)，說明ALL在后伸工況受力最大，最容易受到損傷；PLL在前屈工況下的伸長率最大且明顯高于其他3種工況下伸長率的響應(yīng)，說明PLL在前屈工況受力最大，最容易受到損傷；CL在后伸工況下的伸長率最小，在其他3種工況下伸長率較大且水平相當(dāng)接近13%，說明CL在前屈、軸向轉(zhuǎn)動、側(cè)彎3種工況下都容易受到損傷；ISL在前屈工況下的伸長率最大為21.4%，其他3種工況下的伸長率響應(yīng)均不明顯，說明ISL在前屈工況下更容易損傷。

4 結(jié)論

建立了成人C67頸椎段有限元模型。在幾何信息獲取方面，利用CT掃描、圖像處理和幾何三維模型重構(gòu)的方法，充分保證了頸椎段的幾何精確度；在椎骨網(wǎng)格劃分方面，使用分塊與投影相結(jié)合的方法獲得高質(zhì)六面體椎骨網(wǎng)格，既保證了網(wǎng)格的質(zhì)量又提高了運算的精度和速度；在材料選取方面，利用試驗與文獻相結(jié)合的方法選擇最優(yōu)的材料模型，充分保證模型材料的可靠性。同時，對模型進行了前屈、后伸、軸向轉(zhuǎn)動和側(cè)向彎曲等工況的試驗驗證，驗證結(jié)果表明建立的成人C67頸椎段模型具有較高的生物逼真度。在模型驗證的基礎(chǔ)之上，利用該模型研究了成人C67頸椎段韌帶對不同載荷形式的響應(yīng)特點，結(jié)果表明，在1.0 N·m的前屈、后伸、軸向轉(zhuǎn)動和側(cè)向彎曲等4種載荷作用下：ALL在后伸工況下最容易受到損傷；PLL在前屈工況下最容易受到損傷；CL在前屈、軸向轉(zhuǎn)動、側(cè)彎3種工況下都容易受到損傷；ISL在前屈工況下更容易損傷。C67高質(zhì)量頸椎段模型的建立為后期建立全頸椎模型提供了方法，另外，對C67頸椎段中韌帶在不同載荷形式下的響應(yīng)特點的研究有助于對頸部開展針對性的防護。