人体脊椎腰段有限元建模及其力学分析

陈阳国 刘 彤,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.