三维有限元动物椎体模型的建立及应力分析

冯其金，赵玲娟，谷福顺，魏万利，郑昆仑，郭 锐，李瑞峰

有限元分析是一种运用数字技术模拟真实的物理系统，利用简单并且相互作用的元素，就可以通过利用有限数量的未知量从而得到无限未知量的真实系统。早在40多年前，Belytschko等[1]第一个利用有限元方法应用于脊柱生物力学研究，并取得了一定成果，建立了脊柱三维有限元模型。Hkami等[2]则在分析了动力学和静力学的问题后，在单节腰椎的有限元模型中加入了后部附件。自此以后三维有限元方法分析了腰椎运动节段的生物力学性质愈加普遍。本研究通过CT断层扫描、计算机辅助设计、图像数字化处理，利用有限元软件Mimics10.01建立绵羊单椎体，椎体+间盘，椎体+间盘+前、后纵韧带、棘上韧带三种不同的三维有限元模型，并观察不同载荷下有限元模型应力变化情况。

1 材料与方法

1.1 样本选择 获得新鲜的成年健康绵羊脊柱标本，范围从颈椎到尾椎。剔除脊柱标本上附着的肌肉组织，脊柱椎体序列正常，曲度自然，椎间盘没有明显退变征象。排除脊柱椎体病变及损害。

1.2 样本制作 采用16排螺旋CT机（西门子公司）对绵羊脊柱标本进行常规扫描，扫描范围为T10～Sl椎体。扫描参数：有效层厚0.75 mm，图像重建矩阵512 × 512，图像重建速度16幅秒，像素尺寸为0.313 mm×0.313 mm，旋转时间0.6 s/360°。共扫描时间43 s，螺距1.5，窗宽1500，窗位450。得到DICOM格式的连续断面图像。

1.3 步骤 生物力学有限元分析（ANSYS分析）分为三个阶段：前处理、加载并求解和后处理。前处理就是创建实体模型以及有限元模型，包括创建实体模型、定义单元属性、划分有限元网格、修正模型等几项内容。前处理后，约束边界，加载，通过计算机设置加载载荷。后处理包括通用后处理和时间历程后处理。通用后处理—用来观看整个模块在某一时刻的结果，时间历程后处理—用来观看模型在不同时间段或载荷步上的结果，常用于处理瞬态分析和动力分析结果。

1.4 图像采集 CT断层扫描所得影像保存为DICOM格式文件，将所得的DICOM格式数据通过Mimics10.01软件的输入口导入，可得到横断面图像，Mimics10.01软件可自动生成冠状面图和矢状面图。这三个视窗分别显示三个不同位置，三视图相互关联，用鼠标及定位工具快速定位。

1.5 椎体的三维有限元模型建立 根据建模情况，从DICOM格式文件中筛选合适的CT扫描图像数据。骨与软组织的灰度值在CT图像上区别明显，可以通过阙值分析来提取相应的组织，设置灰度值，设置阙值范围，建立新的Mask。软件可能自动生成骨结构轮廓，椎体骨骼之间是相连的，椎体骨骼图像中存在大量的噪点和空洞（来源于CT图像上阙值差异）。使用Mimics中的Edit mask功能键可消除噪点，填补空洞，去处伪影。

本模型为多椎体模型，每次仅保留需建模的骨Mask，将其不需要的东西全部擦除。对于各部位图像的处理较为复杂，需对横断面、冠状面、矢状面的图像逐一修改。操作前需仔细鉴别椎体间相连的Mask，以决定是否擦除。随后可通过区域增长，做成椎体的几何蒙罩Mask。

留下椎体骨的几何蒙罩Mask，再利用Mimics中的Calculate-3D功能进行计算，计算生成三维模型。此时的三维模型仅是面，不是三维实体，生成的三位椎体表面凹凸不平，选中Mimics软件中的Magic模块中的Smoothing按钮进一步对模型Smoothing处理。

图1 椎体三维模型

经过以上操作，生成一个椎体的三维模型（图1）。对于多个椎体的多联模型，需多次重复操作。

以上是Mimics处理医学软件的第一步，目前状态下生成的椎体三维模型，仅是壳模型或面模型，尚未得出所需要的三维实体模型，无法进行数值模拟分析。下一步继续通过Mimics软件的接口模块，把模型导入逆向工程软件，对这些简单的面模型进行实体化，进一步处理生成高质量的三维实体模型。在Mimics软件中，在菜单栏点Export，将椎体模型分别以Point cloud格式输出并保存。操作每次只能处理一个椎体，需多次重复操作。

在Geomagic Studio中导入要处理的椎体文本文件，Geomagic Studio软件自动把文件打开，以点云格式，选择单位为毫米，出现模型图像。对模型图像进行点云降噪处理，去除多余点云数据，删除多余噪点，模型外形较整齐。

进行统一采样和封装，控制采样次数，在不影响椎体形状的前提下提高计算速度。进入多边形编辑阶段，对模型行松弛、简化、砂纸处理，可以使模型表面更加光滑，棱角过渡更加圆润，同时减少三角片数量。注意砂纸作用，把握其中的“度”。建立曲面片，在曲面片下拉菜单中选择构造曲面片，如果构建的曲面中有干涉或者交叉情况，软件会提醒修改。随后采用手动修改方式，点击曲面片下拉菜单，选择编辑曲面片，曲面片就处于编辑状态，通过用鼠标拖拽曲面片节点位置，调整曲面片形状，尽量保证四边形的规则，构建好曲面片的模型。

构建格栅。在格栅下拉菜单中，选择构建格栅，通过编辑，松弛和简化格栅，使得格栅形状规则化。如果检查到格栅之间有干涉，软件也会自动提醒，此时可通过编辑曲面片，调整曲面片形状来消除干涉，生成格栅后的模型。

在NURBS下拉菜单中选择拟合生成曲面。仅有表面张力值一项需要设置，表面张力值越大，则生成的实体模型表面越光滑。在Geomagic Studio操作中，软件默认操作顺序，即通过树形图来管理构建模型过程。以上步骤先后序不可以打乱，最后形成三维实体模型，分别保存为多个通用iges格式文件。

用ANSYS软件处理。把之前保存好的多个椎体模型iges通用格式文件导入到有限元分析软件ANSYS中，注意一定要把文件导入在ANSYS工作目录下，所用文件名只能是英文和数字，否则ANSYS无法找到并打开文件。操作步骤：打开ANSYS Product Launcher，设置好工作路径，点击运行软件，点击File下拉菜单，选择Import，导入iges格式文件，弹出对话框，导入多个通用iges格式文件，出现导入后的模型，多个椎体保持了原有的相对位置关系，见图2。

图2 椎体的实体模型

图3 间盘的实体模型

1.6 有间盘的椎体三维有限元模型的构建 间盘结构包括终板，纤维环和髓核。间盘的阈值较低，无法在MIMICS中识别，因而创建间盘模型较为困难。文中采用了两种方法：（1）利用椎体的间隙在MIMICS中创建间盘。优点是间盘的创建方法与椎体的重建方法相同，缺点是较难保证与椎体的粘接条件。同时此种方法重建的间盘结构复杂单元数量较多。（2）直接在有限元软件中建立间盘。优点是能够保证与椎体的粘接，结构规则，单元数量少，缺点是对操作人员的要求相对较高，文中采用了在有限元软件中创建椎间盘的方法，间盘的实体模型见图3。

1.7 有间盘和有前、后纵韧带、棘上韧带的椎体三维有限元模型的构建 以带有椎间盘的多椎体三维有限元模型为基础，在应用ANSYS软件支持下，将每个椎体前缘两端及间盘顶点相连接模拟前纵韧带，韧带的直径为5 mm。重复同样操作，将每个椎体后缘两端及间盘顶点相连接模拟后纵韧带。同样操作，完整建立带有椎间盘的和带有前、后纵韧带，棘上韧带的多椎体三维有限元模型，见图4。

文中的韧带采用了空间梁单元，韧带与椎体的连接采用了耦合节点位移的方法。

图4 椎体结构中韧带的有限元模型

2 结果

2.1 间盘对脊柱承载能力的影响 在压缩载荷作用下带有间盘的三维椎体模型的仿真结果见图5和图6。从图5可见，受到压缩载荷时椎体的最大应力位于第五椎体，数值为26.6 MPa，同时可以看出，横突和棘突的应力数值较小。从图6可见椎间盘的最大应力位于第四和第五椎体的椎间盘，数值为29.9 MPa。为了进行对比，作者去掉了椎间盘，用耦合节点位移的方法连接椎体。计算结果显示，锥体的应力变化较小，但是脊柱轴向位移的数值显著下降，说明椎间盘具有增大脊柱压缩变形的作用。

图5 椎体应力分布图

图6 椎间盘应力分布图

2.2 韧带对脊柱承载能力的影响 带有前、后纵韧带的多椎体三维有限元模型在拉伸载荷状况下，脊柱整体，椎体，椎间盘和韧带的应力分布如图7～10所示。从图8可以看出，椎体处局部的应力较高外，大部分区域的应力水平较低，最高应力的数值为205 MPa，位于第五椎体，与压缩载荷出现最大应力的位置接近，作者认为应该属于局部的应力集中。从图9可以看出，椎间盘的应力最大值出现在L3椎体和L4椎体之间的间盘，数值为71.2 MPa，高应力位置与脊柱受到压缩时的位置相同。从图10可以看出，韧带的应力数值较小，前纵韧带高应力位于L3椎体和L4椎体之间的间盘位置，数值为15.1 MPa，后纵韧带的高应力位置靠前，位于L2椎体和L3椎体之间的间盘位置，数值为14.2 MPa。

图7 拉伸载荷下整体应力分布图

图8 拉伸载荷下椎体应力分布图

图9 拉伸载荷下椎间盘应力分布图

图10 拉伸载荷下前、后纵韧带应力分布图

3 讨论

椎体三维有限元模型的建立为脊柱生物力学研究提供了有效的手段，研究者多利用CT扫描图像建立有限元模型，从多个角度、不同方法验证三维有限元分析的准确性、实用性。蔡芳芳等[3]通过机械建模的方式建立了脊柱L1- L5节段的有限元模型，但忽略肌肉和韧带对脊柱的拉力与支撑作用。陈浩等[4]通过扫描健康成人胸腰椎体，建立三维有限元模型，然后对模型施加前屈、后伸、左右侧屈和左右旋转载荷，T12-L2节段的力矩-旋转角度曲线与文献体外实验数据基本吻合。苏晋等[5]利用自编软件和HyperMesh软件，建立了包括椎体、间盘、韧带等结构的腰椎有限元模型，引入接触理论处理上下关节突之间的面-面接触问题，生理载荷下的模拟结果与体外生物力学实验结果更加吻合。秦计生等[6]模拟腰椎受轴向、前弯、侧弯和后伸4种载荷下的生物力学反应，结果表明L4/5间盘的纤维环应变较大，纤维环在各种载荷下均出现较明显的应力集中，纤维环局部应力集中和应变较大是导致其易破裂的原因之一。吴小辉等[7]研究表明三维有限元力学仿真能够有效评估人体胸腰段受力情况，同时为胸腰椎内固定系统稳定性提供理论依据。王建平、顾菊平[8]研究有限元分析结果显示前屈时椎体应力主要集中在椎体前端，而椎间盘的应力分布集中在前外侧。晏礼等[9]所建立的三维有限元模型包括T12-L2椎体，3个椎间盘及前、后纵韧带、黄韧带、棘上、棘间韧带、关节囊韧带等重要结构。在动态冲击试验下观察应力云图，发现应力主要集中于椎弓根及其周围骨皮质，并于椎弓根外下缘与椎体交界处应力最大，沿此处向椎体扩散，而松质骨应力明显小于皮质骨，这与其他研究结果相一致。

季蓉蓉、程学进[10]利用UG软件的缝合及建模相关功能完成椎体及椎间盘实体模型的三维重建，完成了在不同载荷下腰椎L4-L5节段的后伸与侧弯分析，通过与相关的实验结果进行比较，验证模型有效性。徐浩等[11]建立了L3-L5三维有限元模型，该模型根据腰椎的解剖特点模拟了椎体的皮质骨、松质骨、纤维环、髓核、终板、小关节、椎弓根、椎板、横突、棘突、前纵韧带、后终韧带、黄韧带、棘上韧带和棘间韧带等三维结构。对该模型进行前屈、后伸、侧屈和旋转工况下的L3-L5节段活动范围的测定，结果表明该模型在不同工况下的腰椎活动范围，与文献一致，因此认为该模型在一定条件下有效，可以应用于临床和实验研究。

闫家智等[12]建立退变椎问盘L3-L4节段的有限元模型，对模型施加0.3 MP轴向压缩载荷，比较正常和退变椎间盘纤维环、髓核、软骨终板、关节突关节的应力，结果显示与正常椎间盘比较，退变椎间盘纤维环周边承受较大应力，髓核的压应力明显降低，软骨终板的应力集中在外周偏后，关节突关节面的应力明显增大。徐飞等[13]建立腰椎有限元模型，在此模型的椎体表面加载500 N轴向压力以模拟正常人站立时的情况，再分别在15 NM力矩下进行屈、伸、侧弯以及扭转四个动作来验证模型的有效性。建立的腰椎有限元模型各方向的位移符合真实情况，并且椎间盘应力分布接近现实。倪伟峰等[14]研究表明动力固定通过改善一个脊柱节段的活动和载荷转导，可以延缓一个轻微退变的椎间盘的退变过程。Rohlmann等[15]建立了L1-L5三维非线性有限元模型，该模型在L3/4节段模拟装有双侧单节段后路内固定，建有椎间盘模型，在该模型上加载载荷，经研究表明动力内固定减小了固定节段间旋转度，减小了站立状态、伸展位时健康椎间盘内的压力。

本文采用三维重建的方法重建了动物脊柱段的有限元模型，考虑了椎间盘和韧带对脊柱变形和应力分布的影响。研究结果显示，受到载荷作用时，椎间盘具有吸收变形的作用，具有降低椎体的应力和保护椎体的作用。韧带在脊柱拉伸时也会承担一定的载荷，对于降低椎体的椎间盘的应力有一定的作用。脊柱有限元方法的优点是建模方便和进行参数研究方便，解决了进行人在体实验研究违背伦理和进行尸体实验时尸体标本来源困难的问题，缺点是有限元模型存在对于数据的各种简化处理以及假设，不同的研究者处理方式不同，则会对数据分析产生影响。同时，有限元模型属于数据处理方式，即随着数据计算的扩大则存在的偏差与误差就会增大，也会使得有限元分析方法的可靠性受到影响。其次，有限元方法的模型资料多数来自健康人，尸体及解剖数据较少，从而也减弱了其临床可靠性。

有限元分析要想对研究更有意义，就必须与临床资料进行比较相互补充。有限元分析方法虽然还存在着不足，但是随着人们对生物力学的研究逐步深入，有限元技术的发展会更加精准，从而在医学领域发挥其不可替代的作用。

[1] Belytschko TB, Andriacchi TP, Schuhz AB, et al.Analog studies of forces in human spine: computational techniques [J].J Biomech, 1973, 6(4): 361-371.

[2] Hakim NS, King AI.A three dimensional finite element dynamic response analysis of a vertebra with experimental verification [J].J Biomech, 1979, 12(4): 277-292.

[3] 蔡芳芳, 孙东明, 解京明, 等.L1-L5段腰椎简化模型的建立及其有效性验证 [J].第二军医人学学报, 2013, 34(12): 1384-1386.

[4] 陈浩, 张锦洪, 贺增良, 等.脊柱胸腰段三维有限元模型的建立与验证 [J].江苏大学学报(医学版), 2014, 24(1): 72-75.

[5] 苏晋, 赵文志, 陈秉智, 等.建立全腰椎有限元接触模型 [J].医用生物力学, 2010, 25(3): 200-205.

[6] 秦计生, 王昱, 彭雄奇, 等.全腰椎三维有限元模型的建立及其有效性验证 [J].医用生物力学, 2013, 28(3): 321-325.

[7] 吴小辉, 刘小聪, 等.人体胸腰段三维有限元力学的仿真研究[J].微创医学, 2017, 12(3): 315-319.

[8] 王建平, 顾菊平.下腰椎三维有限元模型的建立及应力分析[J].控制工程, 2013, 20(1): 76-83.

[9] 晏礼, 宋文慧, 王春强, 等.胸腰段脊柱有限元模型的建立及其生物力学分析 [J].中国临床研究, 2013, 26(7): 629-631.

[10] 季蓉蓉, 程学进.腰椎L4-L5节段有限元分析 [J].机械设计与制造, 2014, 1(1): 229-231.

[11] 徐浩, 张秋林, 唐昊, 等.腰椎L3-L5三维有限元模型的有效性验证 [J].中国组织工程研究, 2013, 17(35): 6261-6266.

[12] 闫家智, 吴志宏, 汪学松, 等.腰椎间盘退变后应力变化的有限元分析 [J].中国医学科学院学报, 2009, 31(4): 464-467.

[13] 徐飞, 陈安民, 董永辉, 等.腰椎有限元模型建立方法的研究[J].生物骨科材料与临床研究, 2012, 9(6): 5-7.

[14] 倪伟峰, 徐建广.有限元法在脊柱结构和腰椎融合生物力学评价中的应用 [J].中国组织上程研究与临床康复, 2008, 12(30):5949-5952.

[15] Rohhnann A, Burra NK, Zander T, et al.Comparison of the effects of bilateral posterior dynamic and rigid fixation devices on the loads in the umbar spine: a finite element analysis [J].Eur Spine J,2007, 16(8): 1223-1231.