人体C67颈椎段的建模和验证及不同工况下的韧带损伤分析*

宋光辉，王广亮，李大鹏，韩晓强，李志刚

（1.北京交通大学机械与电子控制工程学院，北京 100044； 2.首都医科大学附属北京儿童医院，北京 100045）

前言

我国平均每年交通事故50万起，所导致的死亡人数超过10万人。统计数据表明，每5 min就会有人因交通事故而死亡，每年因为交通事故造成的经济损失达到数亿元。在2011年，我国由于交通事故导致的死亡人数为6.24万人，大约占到世界交通死亡人数的16%［1］。汽车碰撞中颈部尤其是韧带具有较高的损伤风险，统计数据显示颈部的损伤概率在交通事故中呈现上升趋势，且交通事故中颈部的损伤引起后遗症的概率高达50%［2］。因此，研究典型载荷形式下颈部不同韧带响应形式和特点具有重要意义。

尸体试验是生物力学领域研究的黄金法则。目前已有较多文献通过尸体试验研究了成人颈椎的伤害机理和承载情况。Camacho［3］利用 10个成人尸体颈部样本进行了准静态的弯曲试验，得到C7-T1/C67/C56/C45/C34/C23/O-C2角位移响应与不同弯矩之间的定量关系和标准差，同时拟合了角位移与弯矩的关系曲线；Nightingale等人［4］利用16个51-74岁的未防腐处理的成人男性颈椎，共计41个O-C2/C45/C67颈椎段样本，进行了前屈后伸的纯弯曲试验，获得了前屈后伸工况下O-C2/C45/C67的角位移响应和标准差，并拟合了其角位移响应与弯矩的非线性关系；指出在前屈后伸工况下，男性上颈椎段的运动幅度明显大于下颈椎的运动幅度，男性上颈椎段的刚度明显大于女性上颈椎段的刚度；Wheeldon等人［5］利用成人 C23/C34/C45/C56/C67颈椎段样本进行了准静态的弯曲试验，得出成人各颈椎段在不同弯矩作用下的角位移响应的平均值和标准差；Luck［6］对约 5个月到18岁的样本进行了非破坏准静态的弯曲试验和动态拉伸破坏试验等，并利用试验数据研究了颈部的弯曲响应与弯矩的定量关系，研究了婴幼儿到成年全年龄段不同颈椎段拉伸失效力和失效位移与年龄的关系等。这些有限的试验数据为理解颈椎关节的损伤形式和机理提供了良好的基础，同时为有限元模型的验证提供了宝贵的基础数据。

进行尸体试验的难度较大，随着技术的发展，有限元分析已应用到颈椎的损伤机理研究中。Hakin和King等人［7］建立了一块单独椎骨的有限元模型；Yang和King［8］又将这块椎骨模型加以完善，建立了一段包含椎间盘关节面的完整脊椎运动节段模型；Goal等人［9］第一次用梁单元模拟韧带，研究了颈椎在前屈、后仰和轴向转动时的响应；Panzer等人［10］建立了成人C45颈椎段有限元模型，研究了在前屈、后伸、轴向转动、侧向弯曲4种工况下椎间盘、韧带、椎骨关节的载荷分配情况，得到的结论是，在轴向转动和侧向弯曲工况下椎间盘为主要的承载结构，在后伸工况下椎骨关节为主要的承载结构，而在前屈工况下韧带为主要承载结构。Czyz等人［11］建立了人体颈髓有限元模型，并利用猪颈髓试验验证了模型的有效性。付裕等人［12］建立了椎间盘有限元模型，研究纤维环的应力分布。髓核和终板软骨采用四面体网格，髓核采用简单的不可压缩材料，得到在压缩力的作用下，由于纤维环向前膨胀产生的拉应力在椎体的前方最大，屈曲时应力主要集中在后部，扭转时最大应力出现在椎间盘后侧的结论；王芳等人［13］建立了中国人全颈椎有限元模型用于挥鞭伤的损伤分析，其韧带采用线弹性材料，椎间盘不区分髓核、纤维环和纤维基质等结构。曹立波等人［14］建立了首个中国50百分位人体的颈椎有限元模型，并利用颈椎离体坠落试验对模型的有效性进行验证，该模型松质骨和皮质骨都采用各向同性材料，韧带和纤维环采用线弹性材料。

本文中首先建立了C67颈椎段高质量网格有限元模型；然后，对模型进行了前屈、后伸、轴向转动和侧向弯曲等多个工况的尸体试验验证；最后，分析在不同载荷形式下不同韧带的响应和损伤特点。

1 C67颈椎段有限元模型的建立

颈椎的样本从邯郸市中心医院获得，其中获取CT的完整流程得到了伦理委员会的审批，CT数据显示颈椎无损伤和畸变，划分成人C67颈椎段有限元网格模型的流程如图1所示。

图1 成人C67颈椎段有限元网格流程图

1.1 椎骨网格划分

1.1.1 椎骨几何重建和几何区域划分

将颈椎CT扫描获得的DICOM格式文件导入Mimics软件中，进行断层图像分割和三维重建，获得三维模型后，将模型输出为 STL格式并导入Geomagic Studio，对颈椎段表面进行松弛和光滑处理，最终得到椎骨的三维模型。

当模型的几何形状较复杂时，常规的网格划分软件不易划分出高质量规则六面体单元。本文中使用结构化网格划分软件Truegrid（XYZ）进行网格划分，其基本思想是将复杂的三维几何模型分割成若干相对规则的子区域，逐块进行网格划分，然后进行网格对接。对于本文的C67颈椎段，由于横突孔周围的骨骼较细且有更大的曲率变化，为避免相对的投影几何面夹角过大、截面积突变，在分段时必须对其进行细节分割处理。根据各个部位的结构特点，将成人C6分割为24个区域，C7分割为16个区域。椎骨的三维模型和区域划分结果如图2所示。

图2 椎骨的三维模型和区域划分结果

1.1.2 椎骨有限元模型的建立

以C7中一段为例说明网格划分的流程，如图3所示。首先，导入被划分区域的几何文件，并根据分块的大小设定block的网格数量；然后，将建立的block体生成蝶形网格，并通过移动和旋转操作将建立的蝶形网格定位到几何模型中对于网格投影最有利的位置；其次，分别将顶点、边界和面投影到对应几何的点线面上；最后，对网格进行调整、顺滑。

图3 Trigurid块投影流程图

相邻的两块投影完后需要通过网格对接将相邻区域同一个截面上的网格节点合并，图4为C7棘突处相邻两个区域合并节点前后的对比图。

图4 网格对接对比

重复上述步骤依次完成24个C6区域和16个C7区域的网格投影，最终得到高质量C67的椎骨网格（图5）。在其表面提取一层壳单元，得到皮质骨和终板的壳单元网格。

图5 C67椎骨有限元网格

1.1.3 椎间盘有限元模型的建立

椎间盘包括纤维环、髓核和终板软骨。其中纤维环由纤维环基质和纤维环纤维薄膜组成。参照图谱并按照椎体上下几何面，建立贴合椎体上下曲面的终板软骨曲面，按照约50%的截面积对曲面进行分割得出髓核区域［15］。终板软骨和纤维环纤维采用壳单元，髓核和纤维基质采用六面体单元。将纤维环基质的实体网格提取出壳单元，获得各基质层侧边的壳网格。纤维环纤维薄膜共有4层，根据壳网格从内向外的顺序定义为第1至第4层。椎间盘有限元模型如图6所示。椎间盘与椎骨之间通过Contact_tied的方式进行连接。

图6 椎间盘模型

1.1.4 颈椎关节有限元模型的建立

颈椎关节由关节软骨、关节囊韧带及其所包含的关节滑液组成，颈椎关节可用于缓冲两关节面的接触冲击。本文中对关节软骨建立实体单元，使用Contact_tied方式连接到关节突上，定义关节软骨的面面接触，关节囊韧带简化为梁单元。关节软骨的厚度依据Yoganandan等［16］提出的成人关节软骨厚度。关节面模型如图7所示。

图7 关节面网格模型

1.1.5 韧带有限元模型的建立

颈椎韧带较多且在CT扫描图像中难以观察重建。因此，基于解剖学图谱中韧带的位置关系建立韧带结构。韧带的建立一般有两种方式：一种是按其较为真实的几何结构建立韧带；另一种是使用等效的离散梁单元建立韧带。前者能较为准确地反映韧带几何形状，但由于韧带的不规则结构，故其对材料属性要求较高，需要对应部分的应力-应变试验数据。而目前已有的数据多数是韧带宏观的力-位移曲线。考虑到上述限制，本模型中韧带使用一维离散梁单元进行建模。韧带模型如图8所示。

图8 韧带模型

为保证模型与已有试验对象的一致性，成人C67颈椎段不建立肌肉模型。模型建立完成后，进行网格质量检查，发现成人C67模型最小雅各比为0.31，具体的网格数量和质量信息如表1所示。从表1看出，整个模型具有较高的网格质量，能保证后期模型的计算精度与效率。

表1 网格质量信息表

1.2 材料属性

1.2.1 松质骨材料属性

目前，颈部有限元模型中关于松质骨材料的表征尚未有统一的模式，部分文献认为松质骨的材料属性使用各向同性材料［17-19］，另一部分文献认为其属于各向异性材料［20-21］。表2给出了文献中的两种松质骨材料模型。

为获取颈椎松质骨相关力学特点，使用猪的松质骨作为等效样本进行两个方向下的压缩试验，利用圆钻将不同部位下的椎骨按照轴向、径向取出，试验样本的获取过程如图9所示。

表2 松质骨材料模型

对轴向和径向的试验样件进行了准静态压缩试验，试验照片如图10所示。结果表明，松质骨在轴向和径向方向材料属性并不相同，呈现出明显的各向异性的特征，因此模型采用各向异性材料（*MAT_OPTION TROPIC_ELASTIC）进行建模。试验测试结果发现各向异性程度与文献［17］所介绍的成人松质骨两个方向较为接近，因此选取该文献的松质骨材料数据作为最终使用的材料属性。

1.2.2 韧带材料属性

图9 椎骨原始材料获取过程

韧带使用离散梁单元进行模拟，该模型可通过直接输入力-位移曲线用于表征韧带在拉伸过程中拉伸力的非线性变化特点。

图10 样件准静态压缩试验

试验结果表明，韧带的力随位移变化的曲线可近似视为“S”形，包括3个区域。将韧带的拉伸力随应变的变化曲线进行归一化，如图11所示。需要3个控制点，第1个阶段是0点到A点的非线性区域，第2个阶段从A到B的线性区域，第3个区域是B到C的非线性区域。通过A，B，C 3点的数据，进行近似回归，即可得到颈椎所有韧带的力-位移曲线［22］。获得的所有成人韧带的A，B和C点的参数如表3所示。

图11 韧带力-应变归一化后的典型曲线

1.2.3 其他部位材料属性

皮质骨、终板、终板软骨、纤维基质、纤维环纤维、髓核和小关节软骨的材料本构与单元类型如表4所示。

2 模型验证

2.1 成人C67准静态前屈-后伸验证

按照Camacho［3］提出的试验工况进行仿真来验证模型的逼真度。在前屈后伸仿真中，约束C7的所有自由度，然后对 C6分别施加±0.5，±1.0，±1.5，±2.0 N·m的弯矩，得到C6的前屈、后伸角位移响应，如图12和图13所示。

表3 成人相关韧带的A，B和C的数值

可以看出：对于前屈工况，C6前屈角位移响应落在 Camacho和 Nightingale［4］的试验区间内，相对于Wheeldon［5］的试验，有限元模型的前屈角位移响应普遍偏小，存在这种误差的原因可能是试验加载位置测量不准确，仿真加载位置与试验加载的真实位置不一致［10］，也可能是椎骨的中性区与模型的中性区不一致［27］，样本之间的个体差异也可能是误差的来源；对于后伸工况，成人C6角位移响应均落在Camacho，Nightingale和Wheeldon的试验区间内；总体来说，成人C67颈椎段有限元模型能很好地重现成人C67颈椎段的前屈后伸工况。

2.2 成人C67准静态轴向转动和侧向弯曲验证

Panjibi［28］对 C67颈椎段进行了 1.0 N·m力矩作用下的轴向转动和侧向弯曲试验，运动示意图如图14（a）所示。按照Panjibi的试验方式，对有限元模型进行仿真，仿真后得到的C6颈椎段的角位移响应如图14（b）所示。模型轴向转动和侧向弯曲的角位移落在了Panjibi的试验区间内，模型侧向弯曲刚度比Panjibi试验的抗侧弯刚度稍大。总的来说，该模型能很好地模拟成人颈部C67颈椎段的侧向弯曲和轴向转动工况。

表4 颈椎其他部位有限元模型的材料参数

图12 C6前屈角位移响应

图13 C6后伸角位移响应

本文中对成人C67颈椎段进行了前屈、后伸、轴向转动和侧向弯曲等4种工况下的准静态验证。C6角位移响应都落在试验区间内，仅侧向弯曲接近Panjabi试验的区间下界，说明在侧弯工况下，模型的刚度略高，总体来说成人C67颈椎段有限元模型能很好地反映成人C67颈椎段的运动响应，具有较高的生物力学逼真度。

图14 C67颈椎段准静态工况与仿真结果

3 成人C67颈椎段韧带对载荷形式的响应特点

为研究成人C67颈椎段韧带组织的伸长率响应的特点，将C7固定，对C6施加大小为1.0 N·m的前屈、后伸、轴向转动和侧向弯曲等4种载荷。5种韧带在4种载荷下的伸长率如图15所示。

图15 韧带响应

从图15可知：前屈工况时，前纵韧带（ALL）、后纵韧带（PLL）、黄韧带（LF）、关节囊韧带（CL）和棘间韧带（ISL）伸长率分别为 0，1.7%，6.4%，13.05%和21.4%，ISL伸长率最大，CL次之；后伸工况时，ALL和CL伸长率分别为5.16%和0.85%，其他韧带的伸长率几乎为零，ALL的伸长率最大，且明显高于其他韧带的伸长率；轴向转动时，ALL，PLL，LF，CL和ISL伸长率分别为0.59%，0，0.45%，13.47%和1.09%，CL的伸长率最大，为韧带中主要的承载对象；侧向弯曲时，CL的伸长率最大为12.87%，ALL，PLL，LF和ISL 4种韧带的伸长率很小甚至接近于0。

从上述结果还可看出：在1.0 N·m的前屈、后伸、轴向转动和侧弯等4种载荷作用下，ALL在后伸工况下的伸长率最大为5.16%，明显高于其他3种工况下伸长率的响应，说明ALL在后伸工况受力最大，最容易受到损伤；PLL在前屈工况下的伸长率最大且明显高于其他3种工况下伸长率的响应，说明PLL在前屈工况受力最大，最容易受到损伤；CL在后伸工况下的伸长率最小，在其他3种工况下伸长率较大且水平相当接近13%，说明CL在前屈、轴向转动、侧弯3种工况下都容易受到损伤；ISL在前屈工况下的伸长率最大为21.4%，其他3种工况下的伸长率响应均不明显，说明ISL在前屈工况下更容易损伤。

4 结论

建立了成人C67颈椎段有限元模型。在几何信息获取方面，利用CT扫描、图像处理和几何三维模型重构的方法，充分保证了颈椎段的几何精确度；在椎骨网格划分方面，使用分块与投影相结合的方法获得高质六面体椎骨网格，既保证了网格的质量又提高了运算的精度和速度；在材料选取方面，利用试验与文献相结合的方法选择最优的材料模型，充分保证模型材料的可靠性。同时，对模型进行了前屈、后伸、轴向转动和侧向弯曲等工况的试验验证，验证结果表明建立的成人C67颈椎段模型具有较高的生物逼真度。在模型验证的基础之上，利用该模型研究了成人C67颈椎段韧带对不同载荷形式的响应特点，结果表明，在1.0 N·m的前屈、后伸、轴向转动和侧向弯曲等4种载荷作用下：ALL在后伸工况下最容易受到损伤；PLL在前屈工况下最容易受到损伤；CL在前屈、轴向转动、侧弯3种工况下都容易受到损伤；ISL在前屈工况下更容易损伤。C67高质量颈椎段模型的建立为后期建立全颈椎模型提供了方法，另外，对C67颈椎段中韧带在不同载荷形式下的响应特点的研究有助于对颈部开展针对性的防护。