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

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

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

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

1 有限元建模

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

1.1 模型建立流程

(1)CT數據導入Mimics進行三維模型的逆向重建；(2)三維模型導入3-matic進行模型優化；(3)優化模型導入Geomagic Studio進行模型的進一步優化；(4)完整模型導入Ansys進行。

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

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

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

圖2 經過Edit Masks形成的蒙皮Fig.2 The Mask byEdit Masks

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

1.3 脊椎腰段三維模型的優化

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

圖4 3-matic優化的模型Fig.4 The optimized model in the 3-matic

1.4 脊椎腰段三維模型的進一步優化

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

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

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

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

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

2 脊椎腰段有限元分析

2.1 單元類型和材料參數

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

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

圖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椎體上表面上方建立一個質量點，通過質量點對模型施加壓縮、彎曲、扭轉荷載。約束L5下表面的6個自由度，使L5下表面處于完全固定狀態。

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

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

2.3 模型有效性驗證

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

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

2.4 脊椎腰段受力分析

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

2.4.1 施加500 N壓縮荷載時的計算結果

當施加500 N壓縮荷載時，脊椎腰段變形情況如圖10和圖11所示，模型偏心呈后仰趨勢，至上而下位移逐漸減小，最大位移發生在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)是椎間盤纖維環等效應力云圖。所有纖維環在壓縮荷載作用下應力分布較為均勻，隨著脊椎腰段的后仰，在L1-L2椎間盤的纖維環腹側區域開始出現壓應力，隨著后仰幅度的增大，壓應力向下傳遞擴散，同時在L3-L4椎間盤纖維環背部區域出現拉應力，最終在L4-L5椎間盤纖維環腹側出現8.57 MPa的最大拉應力，在L4-L5椎間盤背側出現5.72 MPa的最大拉應力。

圖12 纖維環應力應變云圖和髓核應力云圖Fig. 12 Stress and strain nephogram of fiber ring and stress nephogram of nucleus pulposus

圖12(b)是椎間盤纖維環的應變云圖，整體上所有椎間盤的纖維環應變分布均勻，伴隨著應力的產生，最終在L4-L5椎間盤腹側區域出現3.43的壓應變，在背側出現1.53的拉應變。

2.4.2 施加X方向10 N·m的彎矩時的計算結果

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

圖15(a)是椎間盤纖維環的等效應力云圖，由圖可得，L1-L2和L2-L3椎間盤纖維環應力較為復雜，應力分布主要集中在纖維環背側區域，L3-L4和L4-L5椎間盤纖維環應力分布在纖維環背側和腹側。其中最大應力出現在L3-L4椎間盤纖維環的腹側偏后區域為1.32 MPa，出現該現象的原因可能是椎體生理弧度的變化造成的。

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

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

圖15 纖維環應力應變云圖和髓核應力云圖Fig. 15 Stress and strain nephogram of fiber ring and stress nephogram of nucleus pulposus

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

2.4.3 施加Z方向10 N·m的扭矩時的計算結果

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

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

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

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

圖18 椎體的扭轉角度Fig.18 The torsion angle of the vertebral body

圖19(a)和圖19(b)為纖維環應力應變云圖，由圖可知，纖維環整體應力應變分布較為均勻，應力應變主要分布在纖維環外側，左右兩側的應力應變較大，內側較小。最大應力應變出現在L3-L4椎間盤纖維環右側區域為0.58 MPa和0.151 MPa。

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

圖19 纖維環的應力應變云圖和髓核應力云圖Fig. 19 Stress and strain nephogram of fiber ring and stress nephogram of nucleus pulposus

2.5 脊椎腰段受力分析小結

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

3 結論

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

[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] 汪正宇,劉祖德,王哲,等. 脊柱側凸有限元模型的建立和參數優化[J]. 北京生物醫學工程,2008,(1):28-32,60.

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

[4] 郭立新. 脊椎腰骶關節的有限元模型及其有效性驗證[J]. 中國生物醫學工程學報,2006,(4):426-429.

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

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

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

[8] 李丹. 基于腰椎多層螺旋CT掃描三維形態學分析的腰椎材料、形態及結構屬性變化與骨折相關性的FEA研究[D].吉林長春：吉林大學,2011.

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

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

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

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