人體脊椎腰段有限元建模及其力學(xué)分析

陳陽國 劉 彤,2 王汝恒 陳 科

(1. 西南科技大學(xué)土木工程與建筑學(xué)院 四川綿陽 621010；2. 中國物理研究院培訓(xùn)中心 四川綿陽 621900)

脊椎作為人體的中軸支柱，具有保持人體穩(wěn)定、維持各種運(yùn)動狀態(tài)、承受荷載的功能。腰椎作為脊椎的重要組成部分，至上而下由L1-L5 5個椎體4個椎間盤組成，其位置位于脊椎的下部，在運(yùn)動、負(fù)荷和保護(hù)人體等方面起著重要的作用。國內(nèi)外學(xué)者長期以來為研究人體脊椎的力學(xué)特性作了大量工作。由于脊椎結(jié)構(gòu)的復(fù)雜性，通常的力學(xué)方法無法直接對脊椎進(jìn)行力學(xué)研究，隨著科學(xué)的發(fā)展，有限元分析法成為研究脊椎力學(xué)行為的重要手段之一。Belytschko[1]在1972年第一次提出運(yùn)用有限元分析法作為脊椎力學(xué)的研究方法。經(jīng)過30多年的發(fā)展，有限元分析法在脊椎生物力學(xué)方面的研究已日益成熟，在國外已被應(yīng)用于脊椎模具的開發(fā)以及輔助臨床。我國在這一方面的研究起步較晚，直至2008年汪正宇等[2]才建立了脊椎T1至尾椎的有限元模型。近年來，覃春鈺[3]建立脊椎腰段 L4-L5的有限元模型，分析脊椎腰段在生理載荷下的力學(xué)行為，郭立新等[4]建立了詳細(xì)的人體腰骶關(guān)節(jié) L5-S1的三維非線性有限元模型，為人體脊椎腰骶段關(guān)節(jié)的生物力學(xué)研究和器械植入提供了更為準(zhǔn)確的計(jì)算模型。由于技術(shù)等各方面的因素影響，以往建立的脊椎三維重建模型的精度和準(zhǔn)確性并不高?；谏鲜鑫墨I(xiàn)調(diào)研，本研究通過應(yīng)用 Mimics，Geomagic以及Ansys軟件建立精確的脊椎腰段 L1-L5的有限元模型，模擬脊椎腰段在壓縮、彎曲和扭轉(zhuǎn)荷載作用下的力學(xué)行為。

1 有限元建模

本文基于螺旋CT掃描技術(shù)，選取一名健康男性志愿者脊椎腰段的CT圖像，運(yùn)用Mimics,Geomagic Studio和Ansys軟件建立脊椎腰段有限元模型[5-9]。

1.1 模型建立流程

(1)CT數(shù)據(jù)導(dǎo)入Mimics進(jìn)行三維模型的逆向重建；(2)三維模型導(dǎo)入3-matic進(jìn)行模型優(yōu)化；(3)優(yōu)化模型導(dǎo)入Geomagic Studio進(jìn)行模型的進(jìn)一步優(yōu)化；(4)完整模型導(dǎo)入Ansys進(jìn)行。

1.2 脊椎腰段三維模型的建立

將CT數(shù)據(jù)以dicom文件格式導(dǎo)入Micmic17.0(圖1)，設(shè)定合適的閾值，利用閾值提取工具(Thresholding)提取脊椎腰段輪廓，將提取的輪廓刪除、修補(bǔ)及填充形成蒙皮(mask)(圖2)，將所得蒙皮進(jìn)行三維重建得到脊椎腰段的三維模型，經(jīng)過優(yōu)化操作最后建立較為完整脊椎腰段的三維模型(圖3)。

圖1 健康志愿者脊椎CT螺旋掃描斷層在Mimics17.0軟件中進(jìn)行閾值分割Fig.1 Thresholding segmentation of healthy volunteers’ spinal CT helical scan tomographic segmentation in Mimics17.0 software

圖2 經(jīng)過Edit Masks形成的蒙皮Fig.2 The Mask byEdit Masks

圖3 脊椎腰段的三維模型(不包含椎間盤)Fig.3 The three-dimensional model of the lumbar spine (without intervertebral discs)

1.3 脊椎腰段三維模型的優(yōu)化

由Mimics建立的模型由于表面存在微小的毛刺和孔洞，不利于有限元模型的網(wǎng)格劃分和計(jì)算，因此需要將模型導(dǎo)入3-matic中進(jìn)行模型優(yōu)化。將三維模型導(dǎo)入3-matic，進(jìn)行局部光滑(smooth)、壓縮(push)、張拉(pull)等操作，去除模型表面的毛刺、填補(bǔ)孔洞以及優(yōu)化模型輪廓。經(jīng)過處理后的模型建立面網(wǎng)格。選定三維模型，檢視三角面片質(zhì)量，根據(jù)具體情況適當(dāng)選擇光順(Fix-Smooth)、縮減三角面片(Fix-Reduce)以及手動修改工具，對面網(wǎng)格進(jìn)行優(yōu)化(圖4)。最后將模型文件保存為“STL”格式。

圖4 3-matic優(yōu)化的模型Fig.4 The optimized model in the 3-matic

1.4 脊椎腰段三維模型的進(jìn)一步優(yōu)化

在3D-matic中建立的脊椎腰段模型曲面較為復(fù)雜，部分曲面的曲率較大，不利于有限元計(jì)算，因此需要將模型進(jìn)一步優(yōu)化。

將“STL”格式的脊椎腰段模型文件導(dǎo)入Geomagic Studio，使用“網(wǎng)格醫(yī)生”對模型的曲面片進(jìn)行光滑和刪除釘狀物處理，使用“優(yōu)化邊緣”、“松弛”和“砂紙”等功能進(jìn)一步處理模型。處理脊椎腰段模型的曲面，設(shè)置合理的粒度與曲率級別探測曲率，合理升降約束，構(gòu)造曲面片。采用“編輯曲面片”去除相交路徑和較小的曲面片角度，運(yùn)用“松弛曲面片”處理曲面片后構(gòu)造格柵，重復(fù)“編輯曲面片”使得生成的格柵不含相交格柵，最后擬合曲面，生成NURBS曲面(圖5)，將模型文件保存為“igs”格式。

圖5 NURBS曲面模型Fig.5 The NURBS surface model

1.5 人體脊椎腰段椎間盤重建

上述提及僅僅完成了脊椎腰段椎體的三維重建，而完整的脊椎腰段模型應(yīng)包括椎間盤。椎間盤由纖維環(huán)和髓核組成，并且纖維環(huán)和髓核在椎間盤減緩沖擊和均布荷載方面起到關(guān)鍵作用。因此，本文對椎間盤的有限元建模將在有限元軟件Ansys中完成。隨后，將“igs”格式的模型文件導(dǎo)入Ansys17.2中，通過前處理功能對該模型進(jìn)行布爾操作，分離出上下椎體的上表面和下表面并分割出上下終板結(jié)構(gòu)。以終板為基礎(chǔ)再進(jìn)行一系列的布爾操作，分割出位于終板之間的椎間盤，將椎間盤模型按照實(shí)際椎間盤中纖維環(huán)和髓核的比例進(jìn)行scale命令，分割出含有纖維環(huán)和髓核的接近實(shí)體椎間盤的椎間盤模型(圖6)。重復(fù)布爾操作，得到L1-L5包含5個椎體和4個椎間盤的人體脊椎腰段的有限元模型。至此，人體脊椎腰段的三維實(shí)體模型建立工作全部完成(圖7)。

2 脊椎腰段有限元分析

2.1 單元類型和材料參數(shù)

脊椎腰段各部分的單元類型和材料參數(shù)如表1所示。其中，為模擬髓核不可壓縮的性質(zhì)，設(shè)置了較小的彈性模量和接近0.5的泊松比[10-12]。在劃分網(wǎng)格時，將椎體網(wǎng)格邊長設(shè)置為2 mm，椎間盤和終板網(wǎng)格單元邊長設(shè)置為0.5 mm。由于通過CT數(shù)據(jù)建立的三維模型不規(guī)則，因此對模型采用自由網(wǎng)格劃分，圖8為網(wǎng)格劃分后的模型。

2016年10月15日，教育部部長陳寶生在華中師范大學(xué)召開的武漢高等學(xué)校工作座談會上首次提出，高校要進(jìn)一步轉(zhuǎn)變理念，做到“四個回歸”（分別是回歸常識、回歸本分、回歸初心、回歸夢想）［1］。之后，在新時代全國高等學(xué)校本科教育工作會議等多次會議、談話中，陳寶生部長都反復(fù)強(qiáng)調(diào)了高校（或高等教育）要做到（或推進(jìn)）“四個回歸”?！八膫€回歸”以質(zhì)樸的語言點(diǎn)醒了身處改革發(fā)展浪潮之中的中國高校，及時為我國高等教育的發(fā)展敲響了警鐘，也為貫徹落實(shí)黨的十九大精神和全國教育大會精神提供了基本遵循。學(xué)習(xí)領(lǐng)會、貫徹推進(jìn)“四個回歸”在全國高等教育領(lǐng)域蔚然成風(fēng)。

圖6 椎間盤重建模型Fig.6 The reconstruction model of intervertebral disc

圖7 人體脊椎腰段有限元模型Fig.7 The finite element model of the lumbar spine of human body

2.2 邊界條件

脊椎各部位傳力以均布荷載的形式傳力，因此本研究在L1椎體上表面上方建立一個質(zhì)量點(diǎn)，通過質(zhì)量點(diǎn)對模型施加壓縮、彎曲、扭轉(zhuǎn)荷載。約束L5下表面的6個自由度，使L5下表面處于完全固定狀態(tài)。

表1 脊椎腰段各部分的單元類型和材料參數(shù)Table 1 Unit type and material parameters of each part of the lumbar spine

圖8 網(wǎng)格劃分后的有限元模型Fig.8 The finite element model after mesh generation

2.3 模型有效性驗(yàn)證

本研究通過對模型施加3 mm的位移荷載，得到如圖9所示的荷載位移-曲線與有關(guān)學(xué)者已做試驗(yàn)所得到的荷載位移曲線相似[13]，初步驗(yàn)證了模型的有效性。

圖9 荷載-位移曲線Fig. 9 Load-displacement curves

2.4 脊椎腰段受力分析

本研究通過對模型施加500 N的壓力、10 N·m的彎矩、10 N·m的扭矩模擬脊椎腰段在受到壓力、彎矩和扭矩作用下的力學(xué)行為。

2.4.1 施加500 N壓縮荷載時的計(jì)算結(jié)果

當(dāng)施加500 N壓縮荷載時，脊椎腰段變形情況如圖10和圖11所示，模型偏心呈后仰趨勢，至上而下位移逐漸減小，最大位移發(fā)生在L1椎體棘突上為3.65 mm。可見，由于脊椎特有的生理弧度，脊椎腰段在承受壓縮荷載時伴隨著小幅度的后仰。

圖10 脊椎腰段整體變形Fig.10 The overall deformation of the lumbar spine

圖11 脊椎腰段整體位移場Fig.11 Overall displacement field of lumbar spine

圖12(a)是椎間盤纖維環(huán)等效應(yīng)力云圖。所有纖維環(huán)在壓縮荷載作用下應(yīng)力分布較為均勻，隨著脊椎腰段的后仰，在L1-L2椎間盤的纖維環(huán)腹側(cè)區(qū)域開始出現(xiàn)壓應(yīng)力，隨著后仰幅度的增大，壓應(yīng)力向下傳遞擴(kuò)散，同時在L3-L4椎間盤纖維環(huán)背部區(qū)域出現(xiàn)拉應(yīng)力，最終在L4-L5椎間盤纖維環(huán)腹側(cè)出現(xiàn)8.57 MPa的最大拉應(yīng)力，在L4-L5椎間盤背側(cè)出現(xiàn)5.72 MPa的最大拉應(yīng)力。

圖12 纖維環(huán)應(yīng)力應(yīng)變云圖和髓核應(yīng)力云圖Fig. 12 Stress and strain nephogram of fiber ring and stress nephogram of nucleus pulposus

圖12(b)是椎間盤纖維環(huán)的應(yīng)變云圖，整體上所有椎間盤的纖維環(huán)應(yīng)變分布均勻，伴隨著應(yīng)力的產(chǎn)生，最終在L4-L5椎間盤腹側(cè)區(qū)域出現(xiàn)3.43的壓應(yīng)變，在背側(cè)出現(xiàn)1.53的拉應(yīng)變。

2.4.2 施加X方向10 N·m的彎矩時的計(jì)算結(jié)果

當(dāng)施加X方向10 N·m的彎矩時，脊椎腰段的變形情況如圖13和圖14所示。模型整體發(fā)生前傾，最大位移出現(xiàn)在L1椎體棘突末端為45.82 mm。整個模型的位移場自上而下分布逐漸密集，其中L4,L5椎體和椎間盤位移基本為零，而L1,L2椎體和椎間盤位移很大，其中L1椎體位移場最為密集，最大位移發(fā)生在L1椎體右側(cè)。可見脊椎腰段在彎矩的作用下可以進(jìn)行大幅度的活動，而這一行為的基礎(chǔ)是椎間盤的彈性性能。

圖15(a)是椎間盤纖維環(huán)的等效應(yīng)力云圖，由圖可得，L1-L2和L2-L3椎間盤纖維環(huán)應(yīng)力較為復(fù)雜，應(yīng)力分布主要集中在纖維環(huán)背側(cè)區(qū)域，L3-L4和L4-L5椎間盤纖維環(huán)應(yīng)力分布在纖維環(huán)背側(cè)和腹側(cè)。其中最大應(yīng)力出現(xiàn)在L3-L4椎間盤纖維環(huán)的腹側(cè)偏后區(qū)域?yàn)?.32 MPa，出現(xiàn)該現(xiàn)象的原因可能是椎體生理弧度的變化造成的。

圖13 脊椎腰段整體變形Fig.13 The overall deformation of the lumbar spine

圖14 脊椎腰段整體位移場Fig.14 Overall displacement field of lumbar spine

圖15 纖維環(huán)應(yīng)力應(yīng)變云圖和髓核應(yīng)力云圖Fig. 15 Stress and strain nephogram of fiber ring and stress nephogram of nucleus pulposus

圖15(b)是椎間盤纖維環(huán)等效應(yīng)變云圖，與應(yīng)力分布類似，在L1-L2椎間盤纖維環(huán)背側(cè)出現(xiàn)0.01的壓應(yīng)變，在L3-L4椎間盤纖維環(huán)腹側(cè)偏后區(qū)域出現(xiàn)0.33的拉應(yīng)變。圖15(c)是椎間盤髓核的等效應(yīng)力云圖，髓核應(yīng)力整體分布均勻，L1-L2和L3-L4椎間盤髓核應(yīng)力較大，其中最大應(yīng)力出現(xiàn)在L1-L2椎間盤髓核為0.24 MPa，最小應(yīng)力出現(xiàn)在L3-L4髓核為0.005 MPa。

2.4.3 施加Z方向10 N·m的扭矩時的計(jì)算結(jié)果

當(dāng)在Z方向施加10 N·m的扭矩時，脊椎腰段的變形情況如圖16和圖17所示，模型整體發(fā)生大弧度的扭轉(zhuǎn)，扭轉(zhuǎn)弧度自下而上增大，最大旋轉(zhuǎn)位移矢量發(fā)生在L1椎體棘突處為12.1 mm。位移場分布自下而上依次變密，在L1椎體處最為密集。

如圖18所示，各椎體扭轉(zhuǎn)角度自下而上依次增大，L1椎體的扭轉(zhuǎn)角度最大為46°，L5椎體的扭轉(zhuǎn)角度最小僅為0.02°。可見，椎間盤為椎體能大幅度旋轉(zhuǎn)提供了基礎(chǔ)。

圖16 脊椎腰段的整體變形圖Fig.16 The overall deformation of the lumbar spine

圖17 脊椎腰段整體位移場Fig.17 Overall displacement field of lumbar spine

圖18 椎體的扭轉(zhuǎn)角度Fig.18 The torsion angle of the vertebral body

圖19(a)和圖19(b)為纖維環(huán)應(yīng)力應(yīng)變云圖，由圖可知，纖維環(huán)整體應(yīng)力應(yīng)變分布較為均勻，應(yīng)力應(yīng)變主要分布在纖維環(huán)外側(cè)，左右兩側(cè)的應(yīng)力應(yīng)變較大，內(nèi)側(cè)較小。最大應(yīng)力應(yīng)變出現(xiàn)在L3-L4椎間盤纖維環(huán)右側(cè)區(qū)域?yàn)?.58 MPa和0.151 MPa。

圖19(c)為髓核應(yīng)力云圖，由圖可知，髓核應(yīng)力分布主要分布在髓核腹側(cè)和髓核表面。最大應(yīng)力在L3-L4椎間盤髓核表面為0.054 MPa。

圖19 纖維環(huán)的應(yīng)力應(yīng)變云圖和髓核應(yīng)力云圖Fig. 19 Stress and strain nephogram of fiber ring and stress nephogram of nucleus pulposus

2.5 脊椎腰段受力分析小結(jié)

由于脊椎腰段特有的生理弧度導(dǎo)致其在壓縮荷載作用下偏心，出現(xiàn)小幅度的后仰，在椎間盤纖維環(huán)的背側(cè)和腹側(cè)分別出現(xiàn)壓應(yīng)力和拉應(yīng)力。在彎矩作用下，脊椎腰段出現(xiàn)大幅度的前傾，L1椎體的位移最大達(dá)到45.82 mm。整個位移場分布從下往上依次變密。椎間盤的彈性使得椎體能夠較大幅度的活動。在扭矩作用下，各椎體呈現(xiàn)大幅度扭轉(zhuǎn)的趨勢，L1-L5椎體扭轉(zhuǎn)角度分別為：46°，38°，26.7°，10.7°，0.02°，L1椎體和L5椎體扭轉(zhuǎn)角度相差較大，在本文研究內(nèi)容中椎體近似剛體，進(jìn)一步證明椎間盤是脊椎腰段能夠較大幅度活動的基礎(chǔ)。

3 結(jié)論

本文完成了對脊椎腰段模型的三維重建以及有限元模擬，得出的荷載位移曲線與相關(guān)試驗(yàn)所得出的荷載位移曲線相似，證明了該模型的有效性。本研究得出了以下結(jié)論：(1)脊椎腰段在壓縮荷載作用下因偏心作用會伴隨小幅度的后仰，椎間盤應(yīng)力應(yīng)變分布較為均勻，在纖維環(huán)和髓核的腹側(cè)和背側(cè)區(qū)域會分別出現(xiàn)拉、壓應(yīng)力和拉、壓應(yīng)變。(2)脊椎腰段在彎曲荷載作用下各椎體由下往上出現(xiàn)幅度逐漸增大的前傾。椎間盤應(yīng)力應(yīng)變分布較為復(fù)雜，最大應(yīng)力應(yīng)變在L3-L4椎間盤纖維環(huán)外側(cè)。(3)脊椎腰段在扭轉(zhuǎn)荷載作用下，各椎體呈現(xiàn)大幅度扭轉(zhuǎn)，L1椎體和L5椎體扭轉(zhuǎn)角度相差很大。椎間盤應(yīng)力應(yīng)變主要分布在纖維環(huán)和髓核的外側(cè)，左右兩側(cè)較大，內(nèi)側(cè)較小。(4)脊椎腰段在壓縮、彎曲和扭轉(zhuǎn)荷載作用下，整個模型的位移都較大，在本文研究內(nèi)容中，椎體近似剛體，因此模型的位移主要是因?yàn)樽甸g盤的位移。椎間盤為脊椎腰段在各種受力狀態(tài)下提供活動基礎(chǔ)。

[1] BELYTSCHKO,T B, SCHULTZ,A B .Analog studies of forces in the human spine: computational techniques [J].Journal of Biomechanics, 1973, 6(4):361-371.

[2] 汪正宇,劉祖德,王哲,等. 脊柱側(cè)凸有限元模型的建立和參數(shù)優(yōu)化[J]. 北京生物醫(yī)學(xué)工程,2008,(1):28-32,60.

[3] 覃春鈺. 人體脊椎腰段的三維模型構(gòu)建及有限元力學(xué)分析[D].陜西西安：西安電子科技大學(xué),2014.

[4] 郭立新. 脊椎腰骶關(guān)節(jié)的有限元模型及其有效性驗(yàn)證[J]. 中國生物醫(yī)學(xué)工程學(xué)報(bào),2006,(4):426-429.

[5] 石磊,陸聲,徐永清. 腰椎骨質(zhì)疏松三維有限元模型的建立及應(yīng)用[J]. 西南軍醫(yī),2007,(4):9-10.

[6] 王麗珍. 基于CT掃描之腰椎椎體有限元分析[D].吉林長春：吉林大學(xué),2007.

[7] 李偉. 正常腰椎及腰椎骨質(zhì)疏松三維有限元模型的建立及分析[D].河北石家莊：河北醫(yī)科大學(xué),2011.

[8] 李丹. 基于腰椎多層螺旋CT掃描三維形態(tài)學(xué)分析的腰椎材料、形態(tài)及結(jié)構(gòu)屬性變化與骨折相關(guān)性的FEA研究[D].吉林長春：吉林大學(xué),2011.

[9] 付立會. 椎間盤組織工程的加載裝置設(shè)計(jì)及其有限元建模與分析[D].天津：天津理工大學(xué),2012.

[10] 郭立新,張義民,周淑文,等. 人體腰段脊椎建模及其振動趨勢預(yù)測研究[J]. 系統(tǒng)仿真學(xué)報(bào),2009,21(20):6617-6620.

[11] 許慶慶. 基于ANSYS二次開發(fā)的脊椎模型有限元分析[D].陜西西安：西安電子科技大學(xué),2014.

[12] 張兵. 生理載荷下椎間盤力學(xué)性能的實(shí)驗(yàn)研究及有限元分析[D].天津：天津理工大學(xué),2015.