利用三维绘图软件构建人体腰椎L4-L5节段几何模型, 采用有限元分析方法并调整计算精度, 模拟研究椎间盘发生退变等四种工况, 如轻度退变、中度退变、重度退变以及椎间盘完全摘除时关节峡部的受力情况。通过施加在L4椎板上378.93 N的垂直载荷, 得到四种工况下关节峡部的应力云图。结果表明, 椎间盘退变程度影响小关节处的接触面积, 从而使关节峡部应力发生相应的变化。说明关节峡部是腰椎后部抗压的重要结构, 而关节峡部应力分布与大小以及腰椎后部结构刚度大小和小关节处的接触面积有关。该结论可为脊柱生物力学分析及人工椎间盘置换提供一定的理论参考依据。
引用本文： 胡迎春, 欧亚龙, 胡裔志, 余兵浩. 椎间盘退变对L4-L5在垂直载荷下生物力学行为特性的影响. 生物医学工程学杂志, 2015, 32(1): 55-58, 66. doi: 10.7507/1001-5515.20150010 复制
引言
如今，椎间盘退变性疾病严重影响着人们的正常生活。目前对椎间盘退变性疾病的研究越来越深入，生物力学研究表明，脊柱的力学特性受椎间盘高度、横截面积等形态学参数改变的影响[1-3]，椎间盘高度降低是椎间盘最常发生的形态改变。本研究采用计算机辅助设计(computer aided design, CAD)与计算机辅助工程(computer aided engineering, CAE)相结合的方法，构建人体腰椎L4-L5节段三维模型，对人体椎间盘轻度退变(lightly degenerated disc, LDD)[4]、椎间盘中度退变(moderately degenerated disc, MDD)和椎间盘重度退变(severely degenerated disc，SDD)，以及完全切除椎间盘(completely excision disc，CED)时的四种情况进行模型分析，在轴向压缩载荷条件下，对腰椎后部结构关节峡部等效应力作生物力学分析。
1 研究条件
1.1 研究对象
一例健康男性志愿者，身高168 cm，体重58 kg，既往无腰部疼痛疾病史，无腰部外伤史，X线检查腰椎未见畸形及退变。
1.2 数据采集
采用螺旋CT对自愿者L4椎体上边缘至L5椎体下边缘沿横断面扫描。要求自愿者仰卧位，腰椎位于扫描视野中心。
1.3 扫描条件
层距1 mm，共166层，扫描数据以Dicom 3.0标准直接存储。
2 有限元模型的建立
将腰椎L4-L5节段三角面片信息导入至三维软件中，采用自底向上的方法建立腰椎L4-L5节段几何模型[5-6]。将几何模型导入至Ansys Workbench中，利用有限元软件的自由网格划分功能进行网格划分，根据有限元理论可知, 将模型划分为三角形单元、单元之间以节点连接、内力通过节点在单元之间相互传递。得到腰椎L4-L5节段有限元模型，如图 1所示。
图1
正常腰椎L4-L5节段有限元模型
Figure1.
Finite element models of normal lumbar L4-L5 segment
腰椎间盘是腰椎椎体之间的韧性连接组织[7]，是脊柱运动节段的组成部分，它使脊柱获得一定程度的活动范围并起到传递和分布载荷的作用。椎间盘由中央的髓核和外周的纤维环组成[7]，纤维环包括基质和具有强化作用的胶原纤维，一系列平行的纤维构成的薄层板包埋在基质中，如图 2所示。通过分别改变椎间盘材料特性、降低椎间盘高度，模拟建立椎间盘轻、中、重度退变及完全摘除时的有限元模型。轻度退变椎间盘有限元模型各组分的材料属性是正常高度椎间盘有限元模型髓核弹性模量的2倍[8]。中度退变椎间盘的髓核及纤维环基质弹性模量均设定为正常椎间盘纤维环基质弹性模量的2倍，同时纤维环纤维体积较正常纤维环纤维体积减小25%以模拟纤维环破裂[8]。重度退变椎间盘为中度椎间盘退变的基础上，伴有椎间盘高度降低1/2[8]。各材料属性值如表 1所示[9-10]。椎间盘四种工况下的有限元模型如图 3所示。
图2
椎间盘模型
Figure2.
Intervertebral disc model
图3
椎间盘四种工况下有限元模型
Figure3.
Finite element models of intervertebral disc on four conditions
表1
L4-L5腰椎节段材料特性
Table1.
Material properties of L4-L5 lumbar spine
各接触面设置有摩擦，但关于椎体和椎间盘以及小关节接触的摩擦系数研究尚少[5]。本研究通过选定0~1间选取不同的摩擦系数，在垂直载荷下对模型进行加载分析，得出不同摩擦系数下的椎间盘最大应力及最大位移变化不是很大。最大应力绝对值误差为0.008 MPa，最大位移绝对值误差为0.012 mm。可见摩擦系数对结果影响并不大，可以在0~1之间任意取值。本研究中取摩擦系数为0.2。即在第4腰椎下终板与椎间盘、第5腰椎上终板与椎间盘、小关节接触面处的摩擦系数取为0.2进行后续计算。小关节是摩擦接触，L4上端所受外力，通过摩擦，将外力传递给L5下椎体，通过有限元单元传递内力。
物理上，接触体间不相互渗透，因此程序必须建立两表面间的相互关系以阻止分析中的相互渗透。本研究中两小关节处设置为friction(有摩擦)形式，在Ansys Workbench中，法向刚度系数是影响精度和收敛行为的重要参数，适当调整法向刚度系数，如以体积为主的问题：“Program Controlled”或手动输入“Normal Stiffness Factor”为1；以弯曲为主的问题：手动输入“Normal Stiffness Factor”为“0.001”到“0.1”之间的数值。本研究模型为弯曲为主，因此此时设置的法向刚度系数为0.005。
3 负载和边界条件
将L5下终板面完全约束，即L5下终板所有方向位移及转动完全固定，将后部的棘突x、y方向的位移及x、y方向的转动完全约束，来模拟肌肉和韧带对椎体的约束，如图 4所示。人在直立状态下，第4腰椎受到来自人体2/3体重的载荷。由研究条件中得知人体重58 kg。因此，该体重的2/3为38.67 kg，亦即378.93 N。结合人体实际运动，将378.93 N均匀施加于整个L4椎板上表面。
图4
约束和加载
Figure4.
Constraint and loading
4 结果与分析
4.1 模型验证
经前人研究成果可知，猪与人的椎体大小相似, 标本来源广，可满足各种类型实验的要求[11]。为满足模型外形特征，剔除猪骨韧带，椎旁肌肉，保留椎间盘及关节囊，设计满足约束条件及负载的实验, 在L4-L5关节峡部粘贴应变片，测量得到L5下关节峡部应变值分别为783、792、814、839、844、898、951、1 017，其中静态应变仪上显示测量数据为με级。
von Mises等效应力是基于畸变能密度理论(第四强度理论)适用于塑性材料受压，由屈服准则知：
| $ {{\sigma }_{s}}=\sqrt{\frac{1}{2}\left[{{\left({{\sigma }_{1}}-{{\sigma }_{2}} \right)}^{2}}+{{\left({{\sigma }_{2}}-{{\sigma }_{3}} \right)}^{2}}+{{\left({{\sigma }_{3}}-{{\sigma }_{1}} \right)}^{2}} \right]} $ |
将实验标本近似看成一个柱体，在Ansys Workbench中可以直接读取脊椎L4上终板处横截面积为1 767.6×10-6 m2。将实验标本近似看成一个柱体，由材料力学知识可知，在柱体横截面内，这是一个二向均匀应力状态，在这种情况下，柱体任一点的径向和周向应力值均相等，及σ1=σ2，由脊椎横截面积及垂直载荷可知：
| $ {{\sigma }_{3}}=-\frac{p}{A}=-\frac{378.93N}{1767.6\text{ }\!\!\times\!\!\text{ }{{10}^{-6}}{{\text{m}}^{2}}}=-0.2143\text{MPa} $ |
由应力云图可知理论σ3=15.242 MPa，代入式(1)中得到：σ1=σ2=15.027 7 MPa。其中，皮质骨E=12 000 MPa，μ=0.3通过广义胡可定律可知：
| $ {{\varepsilon }_{1}}=\frac{1}{E}\left[{{\sigma }_{1}}-\mu \left({{\sigma }_{2}}+{{\sigma }_{3}} \right)\right] $ |
代入数值可知ε1=8.819 7×10-4。为了使实验测量数据与理论计算数据在数量级保持一致，由实验结果可知，L5下关节峡部测量均值为, 所以误差。满足误差允许范围，数据合理。
4.2 关节峡部及纤维环基质最大应力值
在大小为378.93 N的垂直载荷作用下，椎间盘四种工况下纤维环基质等效应力分布及大小如图 5所示。由图可知，在同一载荷下，纤维环基质等效应力值随椎间盘退变趋势呈负相关，随着椎间盘退变程度增加而增加。其中重度退变模型在轴向垂直载荷下，纤维环基质等效应力变化尤为明显。
图5
椎间盘四种工况下纤维环基质应力云图
Figure5.
Fiber loop matrix stress nephograms under four conditions
在大小为378.93 N的垂直载荷作用下，椎间盘五种工况下关节峡部应力分布及其大小如图 6所示。研究结果显示，正常椎间盘腰椎关节峡部应力值为15.242 MPa，随着纤维环脱水，纤维环硬度增加，导致椎间盘轻度退变，此时关节峡部应力值为19.103 MPa，腰椎后部刚度减小。随着纤维环进一步脱水，纤维环进一步硬化，椎间盘中度退变时，椎间盘的高度降低，导致椎间孔间隙减小，后部结构中小关节接触处发生错动，使得接触面积增大，从而使得关节峡应力值减小，其值为13.669 MPa，刚度相应增大。随着椎间盘退变进一步加深，达到重度退变时，其高度减小为中度退变时的一半，椎间孔间隙进一步减小，后部结构中小关节接触面积进一步增大，从而使关节峡部所受应力值为10.413 MPa，后部结构刚度进一步增大。但是在完全摘除椎间盘后，椎间孔间隙达到最小，后部结构中小关节处错动严重，导致接触面积由最大降至最小，从而使关节峡部受到很大的应力值，此时为21.441 MPa，此时后部结构刚度降至最小。
图6
椎间盘五种工况下关节峡部应力云图
Figure6.
Joint isthmus stress nephograms under five conditions
5 结论
本研究从理论分析的角度指出关节峡部是腰椎后部承载的主要结构，椎间盘退变导致后部结构刚度变化及峡部应力值大小变化的原因是小关节发生错动，从而导致小关节接触面积的变化。椎间盘退变导致关节峡部处于高负荷工作状态下，长期以往，容易照成小关节滑脱，随着时间推移严重者将会出现关节炎症。本研究结论为脊柱生物力学分析及人工椎间盘置换提供一定的理论参考依据。
引言
如今，椎间盘退变性疾病严重影响着人们的正常生活。目前对椎间盘退变性疾病的研究越来越深入，生物力学研究表明，脊柱的力学特性受椎间盘高度、横截面积等形态学参数改变的影响[1-3]，椎间盘高度降低是椎间盘最常发生的形态改变。本研究采用计算机辅助设计(computer aided design, CAD)与计算机辅助工程(computer aided engineering, CAE)相结合的方法，构建人体腰椎L4-L5节段三维模型，对人体椎间盘轻度退变(lightly degenerated disc, LDD)[4]、椎间盘中度退变(moderately degenerated disc, MDD)和椎间盘重度退变(severely degenerated disc，SDD)，以及完全切除椎间盘(completely excision disc，CED)时的四种情况进行模型分析，在轴向压缩载荷条件下，对腰椎后部结构关节峡部等效应力作生物力学分析。
1 研究条件
1.1 研究对象
一例健康男性志愿者，身高168 cm，体重58 kg，既往无腰部疼痛疾病史，无腰部外伤史，X线检查腰椎未见畸形及退变。
1.2 数据采集
采用螺旋CT对自愿者L4椎体上边缘至L5椎体下边缘沿横断面扫描。要求自愿者仰卧位，腰椎位于扫描视野中心。
1.3 扫描条件
层距1 mm，共166层，扫描数据以Dicom 3.0标准直接存储。
2 有限元模型的建立
将腰椎L4-L5节段三角面片信息导入至三维软件中，采用自底向上的方法建立腰椎L4-L5节段几何模型[5-6]。将几何模型导入至Ansys Workbench中，利用有限元软件的自由网格划分功能进行网格划分，根据有限元理论可知, 将模型划分为三角形单元、单元之间以节点连接、内力通过节点在单元之间相互传递。得到腰椎L4-L5节段有限元模型，如图 1所示。
图1
正常腰椎L4-L5节段有限元模型
Figure1.
Finite element models of normal lumbar L4-L5 segment
腰椎间盘是腰椎椎体之间的韧性连接组织[7]，是脊柱运动节段的组成部分，它使脊柱获得一定程度的活动范围并起到传递和分布载荷的作用。椎间盘由中央的髓核和外周的纤维环组成[7]，纤维环包括基质和具有强化作用的胶原纤维，一系列平行的纤维构成的薄层板包埋在基质中，如图 2所示。通过分别改变椎间盘材料特性、降低椎间盘高度，模拟建立椎间盘轻、中、重度退变及完全摘除时的有限元模型。轻度退变椎间盘有限元模型各组分的材料属性是正常高度椎间盘有限元模型髓核弹性模量的2倍[8]。中度退变椎间盘的髓核及纤维环基质弹性模量均设定为正常椎间盘纤维环基质弹性模量的2倍，同时纤维环纤维体积较正常纤维环纤维体积减小25%以模拟纤维环破裂[8]。重度退变椎间盘为中度椎间盘退变的基础上，伴有椎间盘高度降低1/2[8]。各材料属性值如表 1所示[9-10]。椎间盘四种工况下的有限元模型如图 3所示。
图2
椎间盘模型
Figure2.
Intervertebral disc model
图3
椎间盘四种工况下有限元模型
Figure3.
Finite element models of intervertebral disc on four conditions
表1
L4-L5腰椎节段材料特性
Table1.
Material properties of L4-L5 lumbar spine
各接触面设置有摩擦，但关于椎体和椎间盘以及小关节接触的摩擦系数研究尚少[5]。本研究通过选定0~1间选取不同的摩擦系数，在垂直载荷下对模型进行加载分析，得出不同摩擦系数下的椎间盘最大应力及最大位移变化不是很大。最大应力绝对值误差为0.008 MPa，最大位移绝对值误差为0.012 mm。可见摩擦系数对结果影响并不大，可以在0~1之间任意取值。本研究中取摩擦系数为0.2。即在第4腰椎下终板与椎间盘、第5腰椎上终板与椎间盘、小关节接触面处的摩擦系数取为0.2进行后续计算。小关节是摩擦接触，L4上端所受外力，通过摩擦，将外力传递给L5下椎体，通过有限元单元传递内力。
物理上，接触体间不相互渗透，因此程序必须建立两表面间的相互关系以阻止分析中的相互渗透。本研究中两小关节处设置为friction(有摩擦)形式，在Ansys Workbench中，法向刚度系数是影响精度和收敛行为的重要参数，适当调整法向刚度系数，如以体积为主的问题：“Program Controlled”或手动输入“Normal Stiffness Factor”为1；以弯曲为主的问题：手动输入“Normal Stiffness Factor”为“0.001”到“0.1”之间的数值。本研究模型为弯曲为主，因此此时设置的法向刚度系数为0.005。
3 负载和边界条件
将L5下终板面完全约束，即L5下终板所有方向位移及转动完全固定，将后部的棘突x、y方向的位移及x、y方向的转动完全约束，来模拟肌肉和韧带对椎体的约束，如图 4所示。人在直立状态下，第4腰椎受到来自人体2/3体重的载荷。由研究条件中得知人体重58 kg。因此，该体重的2/3为38.67 kg，亦即378.93 N。结合人体实际运动，将378.93 N均匀施加于整个L4椎板上表面。
图4
约束和加载
Figure4.
Constraint and loading
4 结果与分析
4.1 模型验证
经前人研究成果可知，猪与人的椎体大小相似, 标本来源广，可满足各种类型实验的要求[11]。为满足模型外形特征，剔除猪骨韧带，椎旁肌肉，保留椎间盘及关节囊，设计满足约束条件及负载的实验, 在L4-L5关节峡部粘贴应变片，测量得到L5下关节峡部应变值分别为783、792、814、839、844、898、951、1 017，其中静态应变仪上显示测量数据为με级。
von Mises等效应力是基于畸变能密度理论(第四强度理论)适用于塑性材料受压，由屈服准则知：
| $ {{\sigma }_{s}}=\sqrt{\frac{1}{2}\left[{{\left({{\sigma }_{1}}-{{\sigma }_{2}} \right)}^{2}}+{{\left({{\sigma }_{2}}-{{\sigma }_{3}} \right)}^{2}}+{{\left({{\sigma }_{3}}-{{\sigma }_{1}} \right)}^{2}} \right]} $ |
将实验标本近似看成一个柱体，在Ansys Workbench中可以直接读取脊椎L4上终板处横截面积为1 767.6×10-6 m2。将实验标本近似看成一个柱体，由材料力学知识可知，在柱体横截面内，这是一个二向均匀应力状态，在这种情况下，柱体任一点的径向和周向应力值均相等，及σ1=σ2，由脊椎横截面积及垂直载荷可知：
| $ {{\sigma }_{3}}=-\frac{p}{A}=-\frac{378.93N}{1767.6\text{ }\!\!\times\!\!\text{ }{{10}^{-6}}{{\text{m}}^{2}}}=-0.2143\text{MPa} $ |
由应力云图可知理论σ3=15.242 MPa，代入式(1)中得到：σ1=σ2=15.027 7 MPa。其中，皮质骨E=12 000 MPa，μ=0.3通过广义胡可定律可知：
| $ {{\varepsilon }_{1}}=\frac{1}{E}\left[{{\sigma }_{1}}-\mu \left({{\sigma }_{2}}+{{\sigma }_{3}} \right)\right] $ |
代入数值可知ε1=8.819 7×10-4。为了使实验测量数据与理论计算数据在数量级保持一致，由实验结果可知，L5下关节峡部测量均值为, 所以误差。满足误差允许范围，数据合理。
4.2 关节峡部及纤维环基质最大应力值
在大小为378.93 N的垂直载荷作用下，椎间盘四种工况下纤维环基质等效应力分布及大小如图 5所示。由图可知，在同一载荷下，纤维环基质等效应力值随椎间盘退变趋势呈负相关，随着椎间盘退变程度增加而增加。其中重度退变模型在轴向垂直载荷下，纤维环基质等效应力变化尤为明显。
图5
椎间盘四种工况下纤维环基质应力云图
Figure5.
Fiber loop matrix stress nephograms under four conditions
在大小为378.93 N的垂直载荷作用下，椎间盘五种工况下关节峡部应力分布及其大小如图 6所示。研究结果显示，正常椎间盘腰椎关节峡部应力值为15.242 MPa，随着纤维环脱水，纤维环硬度增加，导致椎间盘轻度退变，此时关节峡部应力值为19.103 MPa，腰椎后部刚度减小。随着纤维环进一步脱水，纤维环进一步硬化，椎间盘中度退变时，椎间盘的高度降低，导致椎间孔间隙减小，后部结构中小关节接触处发生错动，使得接触面积增大，从而使得关节峡应力值减小，其值为13.669 MPa，刚度相应增大。随着椎间盘退变进一步加深，达到重度退变时，其高度减小为中度退变时的一半，椎间孔间隙进一步减小，后部结构中小关节接触面积进一步增大，从而使关节峡部所受应力值为10.413 MPa，后部结构刚度进一步增大。但是在完全摘除椎间盘后，椎间孔间隙达到最小，后部结构中小关节处错动严重，导致接触面积由最大降至最小，从而使关节峡部受到很大的应力值，此时为21.441 MPa，此时后部结构刚度降至最小。
图6
椎间盘五种工况下关节峡部应力云图
Figure6.
Joint isthmus stress nephograms under five conditions
5 结论
本研究从理论分析的角度指出关节峡部是腰椎后部承载的主要结构，椎间盘退变导致后部结构刚度变化及峡部应力值大小变化的原因是小关节发生错动，从而导致小关节接触面积的变化。椎间盘退变导致关节峡部处于高负荷工作状态下，长期以往，容易照成小关节滑脱，随着时间推移严重者将会出现关节炎症。本研究结论为脊柱生物力学分析及人工椎间盘置换提供一定的理论参考依据。
