颈椎及其附属结构具有独特的解剖结构和生理功能，有多根重要神经和血管穿行于躯干和颅脑之间，是整个脊柱最容易发生损伤的区域。一旦受到损害，造成的危害极大，甚至威胁生命。

颈椎区域的生物力学表现较为复杂，因此，选择合适的生物力学模型和研究方法是探求颈椎损伤的诊断评估机制和治疗预后判断的重要途径^[1,2]。近年来，科研人员采用了多种方法进行人体颈椎在体和离体标本力学实验，取得较大进展，且各种生物力学模型也日趋完善。在不同类型的模型中，计算机模型不仅能同时反映外部运动行为和内部力学响应，获得体内、外实验研究无法实现或难以精确的结果，而且也可以为实验研究提供指导性建议^[3,4]。

本课题组基于人体高分辨率螺旋CT 图片数据，采用逆向工程原理和解剖知识构建高精度、非线性三维有限元模型(finite element model)的方法，并通过与相关权威文献比较，以验证该模型的可靠性。现向相关的研究者介绍这种方法。

选择健康女性志愿者1例，年龄46岁，体质量63 kg，身高165 cm，问诊无不适主诉，动、静态触诊及X线检查排除颈椎明显退变、病变和外创史等情况。采用GE Light Speed VCT 64层螺旋CT扫描及断层图像灰度处理。血管造影剂：碘海醇；扫描部位：枕骨底(C₀)上2 mm至第一胸椎(T₁)下2 mm获得

体层图像；CT扫描条件：140 kV、200 mA、层厚0.625 mm。 首先将CT数据导入比利时Materialise公司的交互式医学影像控制系统(Mimics 13.1软件)进行边界拾取，阈值界定为226Hu以上。然后依次经过区域增长、编辑蒙罩、腔隙填补、边界划分等操作进行分割、修补，计算生成颈椎3D模型；然后导入逆向工程软件Geomagic 8.0，历经点阶段、多边形阶段、成形阶段三个处理阶段，转化为NURBS曲面模型；最后将模型导入有限元前处理软件Hypermesh 11.0，进行拓扑分区及网格划分。

根据研究需要，忽略骨髓和胶原纤维影响。骨骼(皮质骨、松质骨)、终板和横韧带采用黏弹性材料，其中皮质骨采用平均厚度为1mm的C3D6单元；松质骨采用C3D4单元；终板采用0.1mm的C3D6单元；横韧带采用M3D4R单元，厚度为0.5mm。椎间盘(髓核、纤维环)和小关节面软骨采用不可压缩的超弹材料，其中椎间盘(含髓核和纤维环)采用增强沙漏控制的C3D8R；小关节软骨采用0.1mm的C3D6单元。各组织结构的材料属性参数详见表1。各材料基于应变能理论Mooney-Rivilin超弹性材料公式运算：W=C₁₀(I₁-3)+C₀₁(I₂-3)+1/d(J-2)²，其中C₁₀、C₀₁为材料常数，I₁、I₂为应力张量的第1、第2不变量。

本研究不涉及破坏，忽略韧带材料的塑性区和破坏区，中性区采用抛物线方程拟合，直线段采用线性方程拟合(图1)，中性区终点 (d_n、f_n) 和失效点定义弹性区的斜率 (d_f、f_f)^[4]。韧带(除横韧带)采用SPRINGA弹簧单元。预设中性区为失效变形的1/3，公式：d_n/d_f=1/3；中性区高度为失效载荷的1/10，公式：f_n/f_f=1/10。本模型设置的韧带参数参考已发表文献^{[5, 6]}。见表2。

虚线代表机械力对韧带作用的载荷—位移曲线；实线代表拟合机械力对韧带的材料属性的载荷—位移曲线，此图非精确比例[4] 图1 机械力对韧带作用的载荷—位移曲线 Fig.1 Force-deformation curves for ligaments from mechanical testing

图选项

所有关节(寰枕关节、寰枢关节、寰椎前弓和齿状突关节、齿状突与横韧带关节以及小关节突关节)的关节面均定义为滑动接触关系，摩擦系数为0.1。

将椎动脉模型导入Hypermesh 11.0软件，划分为内部流体网格与血管壁固体网格；在CFD模块下划分为 ：Inflow(血管入口)、Outflow(血管出口)、Fluid(流体)和Wall(血管壁)，其中Inflow、Outflow和Wall采用S3单元，Fluid采用C3D4单元。血管壁设定为均匀的、各项同性的弹性材料，弹性模量为3 MPa，泊松比为0.45。将椎动脉的血管壁模型导入Ansys 13.0的Fluent模块。流体(血液)属性设定为无梯度变化的层流方式，具有各向同性、不可压缩、恒定密度和黏度的牛顿流体^[13]。流体采用动态网格，离散系数为1。根据生理情况，设颈椎活动时间为0.5 s，时间步数为101。本实验不考虑流场对血管壁的作用力。

首先，设置约束T₁下面的6个自由度。其次，在颅底旋转轴上方选择一参考点，C₀上表面所有单元节点与此参考点的分布耦合(Distribution Coupling)。活动方向参考三维全局坐标(X-Y平面为水平面、X-Z平面为冠状面、Y-Z平面为矢状面)，前屈后伸方向与冠状面平行，旋转方向参考颈曲切线方向，侧弯时垂直于颈曲切线方向并与矢状面平行。设纯扭矩为1.5 Nm，测得颈椎相对活动度(range of motion，ROM)后转化为弧度，换算公式：弧度=角度/180×π，速度值(velocity，v)计算公式：v=弧度/时间。

完成上述全颈椎(固体)模型和椎动脉血流(流体)模型设置后，运行多物理场耦合软件MpCCI 4.2，选择Abaqus和Fluent作为耦合求解器。然后分别将固体模型文件和流体模型文件同步提交，进行计算。

模型导入Abaqus 6.9.1软件，进行运算和验证。在1.0Nm和1.5Nm扭矩下，颈椎模型分别模拟前屈、后伸、左侧屈、右侧屈和左右轴向旋转6种工况^{[9, 14]}。并测算模型各功能节段的活动度。测算方法：前屈、后伸和左侧屈、右侧屈活动角度：测量各椎体上切迹斜率，计算各椎体活动角度，获得相邻椎体活动角度差；左、右旋转活动角度：将棘突与椎体后缘中点连线，计算椎体旋转活动角度。本次建立的颈椎三维有限元模型分别在相同条件设置下与前期文献的研究结果进行比较，即1.0 Nm纯扭矩下C₂-C₇与Panjabi等^[14]对比；1.5 Nm纯扭矩下C₀-C₂与Zhang等^[9]对比。

模型共涉及179577个节点、604399个单元、85个组件、6种单元类型、10种材料和29种材料属性。模型清晰完整地呈现了椎体、椎动脉、椎间盘、韧带、小关节软骨和终板等关键解剖结构，逼真地反映出颈椎及附属结构的空间位置关系^[15]，见图2。

A：全颈椎(C0-T1)三维有限元模型；B：双侧椎动脉及流体模型． 图2 全颈椎和双侧椎动脉及流体三维有限元模型 Fig.2 Finite element model of the cervical spine and the bilateral vertebral artery

图选项

该模型的Mises应力云图直观呈现Von Mises等效应力变化数值，蓝色表示低压力值，红色表示高压力值(图3)，生理活动下的最大应力集中的椎体区域及椎动脉受压情况见表3。颈椎和椎动脉模型在前屈、后伸、左侧屈、右侧屈和左右轴向旋转情况下表面应力的变化见图3。进一步将该模型与前期文献的研究结果进行比较分析，结果显示模型与文献报道的活动度变化趋势一致，差异在可接受范围内。见图4。

各自由度左边为颈动脉模型，右边为椎动脉模型． 图3 颈椎及椎动脉模型在6个自由度活动下表面的应力变化图 Fig.3 Stress changes on the surface of the FEM with corresponding velocity loads under the six degrees of freedom

图选项

图4 各自由度下颈椎有限元模型各节段的活动度比较 Fig.4 Comparison of the ROMs of the cervical spine FEMs' each segment under the six degrees of freedom

图选项

结果显示，椎体在侧弯时上位椎体压力集中于同侧小关节，下位椎体压力多集中于椎体后侧及对侧小关节；屈伸时压力集中于椎体；旋转时压力集中于对侧下关节突关节。血管壁在侧弯时对血管压力最大，后伸次之，前屈和旋转较小。颈椎活动过程中，最大应力几乎均集中在椎动脉的寰椎段(V3段)，即C₂横突孔内的血管壁。

颈椎各活动工况下，血管壁因颈椎活动而发生变形，并对流体产生影响。记录出口血流速度-时间变化数据后，两侧椎动脉最小流速比值及其对应时间见表4。

本文介绍的有限元建模方法将人体颈椎及椎动脉的CT图像数据分别导入相关软件，包括Mimics 13.1、Hypermesh 11.0、Abaqus 6.9.1、Ansys 13.0、MpCCI 4.2等，依次进行固体和流体模型创建、有限元静态模型试算和动态模型的流固耦合计算，实现了生理活动时(前屈、后伸、左侧屈、右侧屈和左、右旋转)的全颈椎及椎动脉表面的力学状态和血管内血流的流场变化。本模型各自由度下活动度与文献结果的趋势一致，但文献未涉及全颈椎(均为部分功能节段)，根据颈椎节段间的叠加效应，本模型活动度数值偏大，因此可以认为其真实有效。

椎动脉是保障脑干、小脑和大脑后部血供的主要血管，约占脑血流量的11.5％。椎动脉的V3段有3～5个生理弯曲，这些弯曲的优势为在颈椎旋转活动时椎动脉具备一定缓冲能力。但这种特殊的解剖结构也会造成血流速度不均匀，易发生血管管壁硬化和堵塞等不良倾向^[16]。另外，由于该区域的椎动脉受结缔组织牵拉固定，若寰枢关节出现不稳定，椎动脉易被明显牵拉或形成折叠样弯曲，导致椎动脉发生栓塞^[17]。

大量临床和基础研究表明，生理状态的椎动脉直接压迫现象与患者的临床表现没有平行关系^[18,19,20]。本研究亦发现血流速度与血管壁应变并非为线性关系，最大应变值通常在极限位出现，而血流速度在颈椎活动时便出现大幅波动。据报道，当人体的椎动脉平均血流速度下降17%～30%，椎-基底动脉供血范围内器官就会因缺血而影响其功能，并产生眩晕发作症状^[21,22]。

虽然在理论上本方法创建的三维有限元模型克服了其他实验模型的某些局限，但实际应用时还需要谨慎地对待真实模型与虚拟模型的问题。在未来研究中，还需要解决下列问题：(1)模型的可验证性与可重复性。本模型参数设置均来源于文献，以便获得更好的重复性和有效性。而如果依据被选取志愿者本身椎动脉内的血流速和舒张压、收缩压的实际测量结果确定模拟的边界和初始条件，这种情况下生物力学结果会有更好的相关性^[23]，但又会带来可靠性验证的问题。因此，如何把握和权衡参数的设置是后续研究的重要方向。(2)模型的简化。动脉血管流固耦合是血流动力学研究的热点和难点，虽然在理论上本模型能够实现颈椎和椎动脉力学关系的仿真，但根据特定研究目的，可能还需要考虑血管几何形状、动脉壁解剖分层、非线性弹性材料性质和三维残余应变、血液的非牛顿性和脉动性、血管壁与血液的相互耦合作用等方面的问题^[24]。另外，肌肉的体外实验数据较为少见^[25]，完全模仿真实颈椎肌肉生物力学规律目前难以实现。

因此，基于充分了解三维有限元模型的优势来制定研究目的，合理调整理论参数和真实参数，得到多类型实验结果的验证，优化资源配置，扬长避短，能够进一步提高全颈椎有限元模型的理论推理能力，也能为临床诊断、治疗以及安全性预测提供参考依据，为深入了解颈椎及其附属结构的生物力学机制提供理想的理论研究平台。