人體C67頸椎段的建模和驗(yàn)證及不同工況下的韌帶損傷分析*

宋光輝，王廣亮，李大鵬，韓曉強(qiáng)，李志剛

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

前言

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

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

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

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

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

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

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

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

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

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

當(dāng)模型的幾何形狀較復(fù)雜時(shí)，常規(guī)的網(wǎng)格劃分軟件不易劃分出高質(zhì)量規(guī)則六面體單元。本文中使用結(jié)構(gòu)化網(wǎng)格劃分軟件Truegrid（XYZ）進(jìn)行網(wǎng)格劃分，其基本思想是將復(fù)雜的三維幾何模型分割成若干相對(duì)規(guī)則的子區(qū)域，逐塊進(jìn)行網(wǎng)格劃分，然后進(jìn)行網(wǎng)格對(duì)接。對(duì)于本文的C67頸椎段，由于橫突孔周圍的骨骼較細(xì)且有更大的曲率變化，為避免相對(duì)的投影幾何面夾角過大、截面積突變，在分段時(shí)必須對(duì)其進(jìn)行細(xì)節(jié)分割處理。根據(jù)各個(gè)部位的結(jié)構(gòu)特點(diǎn)，將成人C6分割為24個(gè)區(qū)域，C7分割為16個(gè)區(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)格，并通過移動(dòng)和旋轉(zhuǎn)操作將建立的蝶形網(wǎng)格定位到幾何模型中對(duì)于網(wǎng)格投影最有利的位置；其次，分別將頂點(diǎn)、邊界和面投影到對(duì)應(yīng)幾何的點(diǎn)線面上；最后，對(duì)網(wǎng)格進(jìn)行調(diào)整、順滑。

圖3 Trigurid塊投影流程圖

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

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

重復(fù)上述步驟依次完成24個(gè)C6區(qū)域和16個(gè)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)纖維薄膜組成。參照?qǐng)D譜并按照椎體上下幾何面，建立貼合椎體上下曲面的終板軟骨曲面，按照約50%的截面積對(duì)曲面進(jìn)行分割得出髓核區(qū)域［15］。終板軟骨和纖維環(huán)纖維采用殼單元，髓核和纖維基質(zhì)采用六面體單元。將纖維環(huán)基質(zhì)的實(shí)體網(wǎng)格提取出殼單元，獲得各基質(zhì)層側(cè)邊的殼網(wǎng)格。纖維環(huán)纖維薄膜共有4層，根據(jù)殼網(wǎng)格從內(nèi)向外的順序定義為第1至第4層。椎間盤有限元模型如圖6所示。椎間盤與椎骨之間通過Contact_tied的方式進(jìn)行連接。

圖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é)面的接觸沖擊。本文中對(duì)關(guān)節(jié)軟骨建立實(shí)體單元，使用Contact_tied方式連接到關(guān)節(jié)突上，定義關(guān)節(jié)軟骨的面面接觸，關(guān)節(jié)囊韌帶簡(jiǎn)化為梁?jiǎn)卧?。關(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)。韌帶的建立一般有兩種方式：一種是按其較為真實(shí)的幾何結(jié)構(gòu)建立韌帶；另一種是使用等效的離散梁?jiǎn)卧㈨g帶。前者能較為準(zhǔn)確地反映韌帶幾何形狀，但由于韌帶的不規(guī)則結(jié)構(gòu)，故其對(duì)材料屬性要求較高，需要對(duì)應(yīng)部分的應(yīng)力-應(yīng)變?cè)囼?yàn)數(shù)據(jù)。而目前已有的數(shù)據(jù)多數(shù)是韌帶宏觀的力-位移曲線。考慮到上述限制，本模型中韌帶使用一維離散梁?jiǎn)卧M(jìn)行建模。韌帶模型如圖8所示。

圖8 韌帶模型

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

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

1.2 材料屬性

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

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

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

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

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

1.2.2 韌帶材料屬性

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

韌帶使用離散梁?jiǎn)卧M(jìn)行模擬，該模型可通過直接輸入力-位移曲線用于表征韌帶在拉伸過程中拉伸力的非線性變化特點(diǎn)。

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

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

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

1.2.3 其他部位材料屬性

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

2 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

圖15 韌帶響應(yīng)

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

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

4 結(jié)論

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