枕寰樞關(guān)節(jié)三維有限元模型的建立及其有效性驗(yàn)證

牛雷，賀瑞，張文志，段麗群，李旭

(安徽省立醫(yī)院，安徽 合肥 230000)

寰樞關(guān)節(jié)為一復(fù)合關(guān)節(jié)[1]，形態(tài)結(jié)構(gòu)特殊，功能復(fù)雜，由寰樞側(cè)塊關(guān)節(jié)和齒狀突前后關(guān)節(jié)組成，寰樞關(guān)節(jié)僅僅靠韌帶和關(guān)節(jié)囊連接。對這些韌帶、關(guān)節(jié)囊維持寰樞椎旋轉(zhuǎn)穩(wěn)定作用的準(zhǔn)確分析，有助于理解寰樞關(guān)節(jié)旋轉(zhuǎn)脫位的機(jī)制和寰樞關(guān)節(jié)旋轉(zhuǎn)脫位后生物力學(xué)特性變化，從而為預(yù)防、診斷和治療提供有意義的信息。新鮮尸體模型接近人體內(nèi)生理環(huán)境，是研究寰樞椎旋轉(zhuǎn)脫位的生物力學(xué)性能的理想材料。但尸體容易“疲勞”，重復(fù)性差，難以收集，價(jià)格昂貴，限制了其應(yīng)用，且尸體標(biāo)本生物力學(xué)特性與新鮮標(biāo)本存在較大差異[2]。

有限元分析可以模擬體內(nèi)特定結(jié)構(gòu)的力學(xué)特征，計(jì)算應(yīng)力分布及應(yīng)變，較物理模型或尸體模型有力學(xué)性能評估全面、重復(fù)性好、實(shí)驗(yàn)時(shí)間短、費(fèi)用少等優(yōu)點(diǎn)；已有報(bào)道構(gòu)建了枕寰樞復(fù)合體有限元模型[3-5]，為枕頸交界區(qū)疾病復(fù)雜的生物力學(xué)研究以及旋轉(zhuǎn)脫位機(jī)制的探討帶來了新的機(jī)遇。但上述的模型結(jié)構(gòu)簡單，采用單軸彈簧單元建立韌帶結(jié)構(gòu)模型，簡化了韌帶-骨骼之間的連接關(guān)系，損失了模型的幾何信息，也不能計(jì)算不同運(yùn)動狀態(tài)下韌帶的應(yīng)變及應(yīng)力分布，不利于模型的進(jìn)一步深入發(fā)展應(yīng)用。

在本研究中，我們構(gòu)建包含部分枕骨以及寰樞椎椎體的三維非線性有限元模型，并添加所有的關(guān)節(jié)軟骨、關(guān)節(jié)囊以及韌帶結(jié)構(gòu)，對韌帶及關(guān)節(jié)囊設(shè)計(jì)采用具有二維結(jié)構(gòu)的膜單元，較其他研究優(yōu)化了模型的幾何形狀真實(shí)度，并且可以進(jìn)行韌帶結(jié)構(gòu)應(yīng)力變化及應(yīng)力分布的計(jì)算，通過對模型施加靜態(tài)力矩載荷使其產(chǎn)生三維六自由度的角位移活動，并與文獻(xiàn)中尸體標(biāo)本實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，對模型有效性進(jìn)行驗(yàn)證。

1 資料與方法

1.1 掃描對象 男性，年齡24歲，身高174 cm，體重65.5 kg，既往無影響枕寰樞復(fù)合體運(yùn)動的病史；影像學(xué)篩查各種測量指標(biāo)都在正常范圍內(nèi)，寰齒前間距在過伸過屈位片變化值小于1 mm。

1.2 CT掃描 使用由安徽省立醫(yī)院影像科提供的GE公司的32排多層螺旋CT對志愿者行枕外隆凸上緣和中顱窩鞍背以上至C7-T1椎間盤的水平面連續(xù)軸位薄層掃描。共獲得枕頸部軸位掃描圖像968層，掃描圖像數(shù)據(jù)轉(zhuǎn)換為標(biāo)準(zhǔn)DICOM格式，保存后直接刻錄光盤存儲。

1.3 骨性結(jié)構(gòu)的建立 將CT圖像導(dǎo)入Mimics軟件中，采用軟件默認(rèn)的CT骨質(zhì)閾值范圍226～3 071，自動定位形成矢狀位、冠狀位和軸位三個(gè)方向的各層骨組織輪廓的曲線，針對每個(gè)圖像進(jìn)行編輯，刪除不需要的區(qū)域，去除離散噪點(diǎn)。根據(jù)矢、冠、軸三個(gè)方向的二維圖像，刪除樞椎C2以下、枕骨髁外側(cè)、頸靜脈孔前方、枕外隆凸以上的圖片，再分別建立不同mask，并將分離的每塊骨骼分別放入其中，獨(dú)立進(jìn)行建模，可以得到各自的渲染模型(見圖1)。

1.4 劃分網(wǎng)格，添加韌帶 將骨骼的模型導(dǎo)入Hypermesh軟件，生成包括部分顱骨、完整寰椎及樞椎的三維實(shí)體模型。并進(jìn)一步劃分網(wǎng)格。根據(jù)頸椎韌帶的解剖學(xué)結(jié)構(gòu)，并參考文獻(xiàn)中報(bào)道的各韌帶的起始點(diǎn)、附著位置[5-6](見表1)，在枕骨髁(C0～C2)骨性模型基礎(chǔ)上根據(jù)各韌帶起止點(diǎn)及橫截面積，補(bǔ)充建立相關(guān)韌帶結(jié)構(gòu)。各韌帶均以膜單元模擬，厚度1 mm，材料屬性：各向異性粘彈性，根據(jù)參考文獻(xiàn)賦予各結(jié)構(gòu)相應(yīng)的材料參數(shù)[5～7](見表1)。

在Hyper Mesh軟件中完成了13條頭頸部韌帶的建模，對模型中關(guān)節(jié)接觸：寰枕關(guān)節(jié)、寰樞關(guān)節(jié)、齒狀突前后的關(guān)節(jié)軟骨及關(guān)節(jié)之間的接觸，接觸屬性定義為無摩擦的面面接觸[8-9]，13條韌帶在附著點(diǎn)與椎體或枕骨緊密連接，故本文綁定約束對其進(jìn)行設(shè)置。最終模型見圖2。

體網(wǎng)格生成：將Hyper Mesh軟件中生成的文件導(dǎo)入“ABAQUS”中，生成模型面網(wǎng)格模型，各椎體骨骼又被分為皮質(zhì)骨及松質(zhì)骨。皮質(zhì)骨被賦以殼單元，松質(zhì)骨采用四面體單元進(jìn)行網(wǎng)格劃分，為達(dá)到更好的解剖真實(shí)性，對韌帶及關(guān)節(jié)囊結(jié)構(gòu)采用膜單元進(jìn)行劃分。根據(jù)文獻(xiàn)報(bào)道設(shè)定皮質(zhì)骨、松質(zhì)骨及韌帶材料參數(shù)[9-12]。

1.5 模型有效性的驗(yàn)證 對所建立的C0～C2有限元模型進(jìn)行前屈、后伸、左右側(cè)屈、左右旋轉(zhuǎn)6個(gè)方向運(yùn)動狀態(tài)進(jìn)行驗(yàn)證，加載1.53 N·m的扭矩，運(yùn)用有限元軟件Abaqus進(jìn)行計(jì)算，測量并比較相同載荷下正常上頸椎模型在屈、伸、側(cè)屈、旋轉(zhuǎn)等情況下寰樞椎的活動度，與Panjabi 9實(shí)驗(yàn)結(jié)果對比驗(yàn)證。

2 結(jié) 果

本研究建立的寰樞椎模型，除了包括骨性結(jié)構(gòu)外，還包括關(guān)節(jié)囊、韌帶等復(fù)雜的軟組織結(jié)構(gòu)，共有513 601個(gè)網(wǎng)格單元，節(jié)點(diǎn)共計(jì)1 003 183個(gè)(見表2，見圖3)。

對模型施加靜態(tài)載荷以進(jìn)行三維六自由度的運(yùn)動。當(dāng)施加載荷達(dá)到1.53 Nm時(shí)，將枕寰樞復(fù)合體活動度和Panjabi參考值對比顯示，本模型各個(gè)方向運(yùn)動范圍均位于尸體實(shí)驗(yàn)結(jié)果估計(jì)參考值95%范圍內(nèi)(見表3)。

同時(shí)根據(jù)動態(tài)數(shù)據(jù)，描繪不同運(yùn)動方向的力矩-角位移曲線；橫坐標(biāo)為加載力矩的大小，縱坐標(biāo)為力矩對應(yīng)的角位移(見圖4)。由圖4可見，四個(gè)方向運(yùn)動的轉(zhuǎn)矩-旋轉(zhuǎn)角曲線都表現(xiàn)非線性關(guān)系，此非線性關(guān)系與尸體實(shí)驗(yàn)結(jié)果[1，8，13-16]基本一致。根據(jù)文獻(xiàn)中寰齒間距的測量方法[17]，測得本模型中立位寰齒間距是1.2 mm，在不斷增大向前剪切力的作用下，動態(tài)記錄寰齒間距大小的變化，本模型在橫韌帶臨界失效狀態(tài)下，寰齒間距至2.8 mm，也與尸體標(biāo)本實(shí)驗(yàn)一致。

3 討 論

3.1 有限元軟件的特點(diǎn) 隨著有限元方法的迅速發(fā)展[18]，有限元輔助分析軟件的功能也越來越全面，但他們的側(cè)重點(diǎn)不同，如ANSYS應(yīng)用涵蓋了電磁、流體和等多種物理耦合場。涵蓋領(lǐng)域較廣；ABAQUS對解決結(jié)構(gòu)力學(xué)尤其是在處理非線性領(lǐng)域問題時(shí)具有很強(qiáng)的優(yōu)勢，但覆蓋物理場范圍較；%90%90狹窄，基于尸體研究證明枕寰樞結(jié)構(gòu)具有明顯的非線性特點(diǎn)，我們選取了擁有突出非線性求解能力的ABAQUS建立和運(yùn)算模型。

圖1 椎體渲染模型 圖2 完整枕寰樞復(fù)合體模型示意圖

表1 各韌帶的起止點(diǎn)及材料參數(shù)

表2 C0、C1、C2體網(wǎng)格信息

表3 模擬加載1.53 N·m時(shí)得到的枕寰樞復(fù)合體活動度和Panjabi參考值(°)

a 骨性結(jié)構(gòu)網(wǎng)格 (后面觀) b 骨性結(jié)構(gòu)網(wǎng)格 (側(cè)面觀) c 韌帶結(jié)構(gòu)網(wǎng)格 (前面觀) d 韌帶結(jié)構(gòu)網(wǎng)格 (后面觀) e 韌帶結(jié)構(gòu)后面觀 (剖視) f 韌帶結(jié)構(gòu)后面剖視(去除覆膜)

圖3 劃分網(wǎng)格后的模型示意圖

a 前驅(qū)運(yùn)動 b 后伸運(yùn)動 c 旋轉(zhuǎn)運(yùn)動 d 側(cè)屈運(yùn)動

圖4 模型前屈、后伸、旋轉(zhuǎn)和側(cè)屈運(yùn)動的力矩-角位移曲線

3.2 枕寰樞復(fù)合關(guān)節(jié)研究有限元建模進(jìn)展 枕寰樞復(fù)合關(guān)節(jié)結(jié)構(gòu)復(fù)雜，有限元建模也經(jīng)歷了由單一椎體到完整有韌帶結(jié)構(gòu)的模型的過程[14]，既往研究由于側(cè)重點(diǎn)不同，對復(fù)合體解剖結(jié)構(gòu)的建模方法不盡相同，對材料的定義也不一致，限制其廣泛應(yīng)用。周圍韌帶在枕寰樞復(fù)合體中非常重要，由于韌帶結(jié)構(gòu)在CT圖像上均不能清晰顯示，文獻(xiàn)中其構(gòu)建方法大都是參考解剖學(xué)文獻(xiàn)，在骨性模型上逐一添加，橫韌帶前方與齒狀突相關(guān)節(jié)，限制其后移，在韌帶建模過程中具有特殊性[3，7]，大多數(shù)的文獻(xiàn)都是以固體單元來模擬橫韌帶，和齒狀突的關(guān)系以關(guān)節(jié)之間的“接觸對”來建立。就目前的有限元研究而言，很少有有限元模型去研究寰樞椎的復(fù)合性損傷。因此，構(gòu)建一個(gè)合理的寰樞椎復(fù)合體三維有限元模型并在此基礎(chǔ)上研究其損傷機(jī)制是非常有必要的。

3.3 本模型的特點(diǎn) 此研究為使骨骼結(jié)構(gòu)更貼近實(shí)際，在重建過程中，對各層圖像均僅作平滑處理，使其幾何形狀高度逼真，是后續(xù)成功驗(yàn)證有效性的基礎(chǔ)。對于韌帶模型的建立，采用了具有二維結(jié)構(gòu)的膜單元，和既往研究彈簧單元相比，不僅優(yōu)化了模型的幾何形態(tài)真實(shí)度，還可分析各韌帶在不同運(yùn)動狀態(tài)下的應(yīng)變及應(yīng)力分布，以分析其損傷機(jī)制，這是以前的枕寰樞復(fù)合體三維有限元模型所不能達(dá)到的。

對模型有效性驗(yàn)證，和尸體標(biāo)本實(shí)驗(yàn)采用一樣的方法，計(jì)算其角位移隨力矩大小動態(tài)變化，驗(yàn)證其仿真程度。由本模型的力矩-角位移曲線可以看出，在各種運(yùn)動的運(yùn)動初期，較小的載荷即可引起較大范圍的活動，曲線較陡直，當(dāng)運(yùn)動范圍增大到一定程度，曲線較為平坦，從而證明模型的非線性特征明顯，說明模型符合人體上頸椎的真實(shí)運(yùn)動規(guī)律；通過和Panjabi的結(jié)果進(jìn)行對比，我們模型在前屈后伸、左右旋轉(zhuǎn)等方面均和Panjabi的結(jié)果類似，說明模型的運(yùn)動學(xué)特征較為可靠。

我們建立的有限元模型，韌帶及關(guān)節(jié)囊具有二維特征，更加逼真；通過運(yùn)動學(xué)特征進(jìn)行驗(yàn)證，具有上頸椎典型的非線性特征，符合人體上頸椎的真實(shí)運(yùn)動規(guī)律，在一定程度上可用于寰樞椎脫位模型建立等后續(xù)研究。

參考文獻(xiàn)：

[1]Izzo RG，Ambrosanio A.Biomechanics of the spine Ⅲ：the cranio-cervical junction[J].Neuroradiol J，2007，20(2)：209-217.

[2]韓應(yīng)超.上頸椎韌帶對寰樞椎穩(wěn)定性影響的生物力學(xué)研究[J].中華骨科雜志，2013，33(6)：628-634.

[3]Brolin K，Halldin P.Development of a finite element model of the upper cervical spine and a parameter study of ligament characteristics[J].Spine，2004，29(4)：376-385.

[4]Graham RS，Oberlander EK.Validation and use of a finite element model of C2for determination of stress and fracture patterns of anterior odontoid load[J].J Neurosurg，2000，93(1 Suppl)：117-125.

[5]任中武，倪斌.上頸椎三維非線性有限元模型的建立及其有效性驗(yàn)證[J].脊柱外科雜志，2007，5(3)：159-162.

[6]陳金水，倪斌.寰樞椎脫位三維非線性有限元模型的建立和分析[J].中國脊柱脊髓雜志，2010，20(9)：749-753.

[7]Wang H，Chen B.Application of finite element analysis in Chinese cervical manipulation biomechanics[J].Shengwu Yixue Gongchengxue Zazhi，2013，30(5)：1123-1126.

[8]Rocha R，Sawa AG.Atlantoaxial rotatory subluxation with ligamentous disruption：A biomechanical comparison of current fusion methods[J].Neurosurgery，2009，64(3 Suppl)：137-143.

[9]Zhang H，Bai J.Development and validation of a finite element model of the occipito-atlantoaxial complex under physiologic loads[J].Spine，2007，32(9)：968-974.

[10]Ha SK.Finite element modeling of multi-level cervical spinal segments(C3～C6)and biomechanical analysis of an elastomer-type prosthetic disc[J].Med Eng Phys，2006，28(6)：534-541.

[11]Ng HW，Teo EC.Nonlinear finite-element analysis of the lower cervical spine(C4～C6)under axial loading[J].J Spinal Disord，2001，14(3)：201-210.

[12]Panzer MB，Cronin DS.C4～5segment finite element model development，validation，and load-sharing investigation[J]，J Biomech，2009，42(4)：480-490.

[13]Panjabi M，Dvorak J.Three-dimensional movements of the upper cervical spine[J].Spine，1988，13(7)：726-730.

[14]Huggare J，HoughtonP.Associations between atlantoaxial and craniomandibular anatomy[J].Growth Dev Aging，1996，60(1)：21-30.

[15]Pang D，Li V.Atlantoaxial rotatory fixation：Part 1：biomechanics of normal rotation at the atlantoaxial joint in children[J].Neurosurgery，2004，55(3)：614-625.

[16]JainVK.Atlantoaxial dislocation[J].Neurol India，2012，60(1)：9-17.

[17]Tubbs RS，Salter EG.The Accessory atlantoaxial ligament[J].Neurosurgery，2004，55(2)：400-402.

[18]Trivedi S.Finiteelement analysis：A boon to dentistry[J].J Oral Biol Craniofac Res，2014，4(3)：200-203.