摘要
目的 建立T2～L5胸腰椎有限元模型并验证其有效性,为探究脊柱冲击载荷下的动态响应特性及损伤机制提供数值模型支撑。方法 基于CT扫描图像数据建立T2～L5胸腰椎三维有限元模型;仿真分析施加不同力矩下(屈伸、旋转和侧弯工况)T12～L1段载荷-转角曲线,并与文献报道的数据进行对比;对T2～6、T7～11和T12～L5 3段脊柱有限元模型施加不同高度下的自由落体载荷并进行仿真分析,获得轴向力峰值和弯矩,并与文献报道的数据进行对比分析。结果 T12～L1段脊柱有限元模型受不同方向力矩发生最大转角在-2.24°～1.55°,与文献数据吻合良好。在不同跌落高度下,T2～6、T7～11和T12～L5 3段脊柱有限元模型的轴向峰值力分别为1.7～5.3、1.3～5.5、1.3～7.5 kN,均处于文献数据误差范围内;脊柱与椎间盘应力云图显示,椎体由外缘最先受力,椎间盘由髓核承受主要载荷,符合实际脊柱损伤发生机制。结论 所建立的T2～L5脊柱模型能够正确模拟不同工况下脊柱的生物力学行为特性,分析结果具有有效性。
关键词:
胸腰椎脊柱
冲击载荷
脊柱损伤
动力学分析
静力学分析
脊柱生物力学
脊柱是人体支柱,具有复杂的结构和功能 。随着现代科技的发展,在交通事故、载人飞行、高空坠落等工况下产生脊柱创伤的患者越来越多。因此,研究脊柱生物力学对于脊柱损伤防护、医疗救治、器械研发等具有重要意义 。有限元法已被公认为是研究生物力学最有效的方法之一,并被广泛用于脊柱生物力学研究 。有限元模型效率高、成本低、响应特征参数提取方便,并且避免了实验相关的伦理问题 。已有不少学者采用有限元法研究脊柱在外力作用下的力学性能及响应规律 [5-11] 。有限元法的脊柱生物力学研究可分为静力分析与动力分析。静力分析主要探究力和力矩静态加载对脊柱的影响 ,动力分析多为探究各类冲击载荷对脊柱造成的损伤 [9-11] 。临床上常见的胸腰椎骨折损伤主要受垂直冲击导致。研究胸腰椎骨折损伤过程中的应力分布及损伤机制,需在有限元模型上模拟动态过程,进行冲击实验和动态响应分析。
目前,胸腰椎骨折损伤的研究集中在腰椎、胸腰椎交界处 [12-14] 。而基于T2～L5节段完整的可进行动态冲击实验分析的胸腰椎三维非线性有限元模型鲜有报道,目前缺少动态载荷下人体整段脊柱响应特性和损伤机制的有限元研究,特别是对于人体在保持直立状态发生跌落工况下,脊柱受到冲击后的动态响应特性与损伤机制尚不清楚。本文通过CT扫描获取正常成年男性志愿者的脊柱影像数据,创建T2～L5段脊柱三维有限元模型,通过静态和动态加载仿真计算及与已发表结果的比较分析,验证模型的有效性。
选取1位中国50百分位体型健康男性志愿者(身高172 cm,体质量56 kg,年龄46岁),通过CT扫描获得其胸腹部影像数据,扫描层厚为0.5 mm。椎间盘结构与韧带、小关节软骨等组织几何数据无法通过CT扫描获得,需要在用于网格划分的HyperMesh 13.0中生成。椎间盘结构由髓核和纤维环构成,应用椎体上、下表面创建椎间盘,进一步划分出髓核与纤维环,髓核体积约占椎间盘体积的44% ;基于上椎骨下关节突外侧与下椎骨上关节突内侧表面创建相邻椎骨之间关节软骨。T2～L5椎体、椎间盘、关节软骨均采用长度为1～4 mm的C3D10单元进行划分;考虑脊柱6条周围韧带,采用轴向只拉不压的非线性弹簧单元进行仿真,根据人体解剖图谱确定韧带位置。不同组织之间的边界使用公共节点进行连接,不需要其他运动约束。最终建立T2～L5脊柱模型如 图1 所示。
将划分好网格的脊柱骨骼模型导入Mimics 19.0中,应用灰度赋值法根据经验公式定义材料属性 [13,15] ：
式中: ρ 为CT图像某一点密度; Hu 为椎体CT值; E 为弹性模量。假设不同灰度值下的椎体材料均为各向同性,遵循应变率相关的Johnson-Cook本构模型 ,模型恒温下等效应力计算公式为:
式中: σ 为等效应力; A 为屈服应力; B 为硬化模量; C 为应变速率常数; n 为硬化指数; 为当前应变率; 为参考应变率。
椎间盘材料遵循一阶Mooney-Rivlin超弹模型,应变能 W 定义如下:
W = C 10 ( I 1 -3)+ C 01 ( I 2 -3)+
式中: C 10 、 C 01 为材料常数; I 1 、 I 2 为第一、第二Green应变不变量; D 1 为材料不可压缩常数。有限元模型材料参数见 表1
考虑应变速率对椎间盘力学性能有重要影响,本文椎间盘采用两组材料参数赋值。低应变率用于静态或准静态分析,高应变率用于高应变率动态分析,如碰撞、冲击等 [17-18] 。
韧带根据Chazal等 通过实验得出的非线性应力-应变曲线赋值。实验中获得的前纵、后纵、横向、棘间、棘上和黄韧带横截面积分别为35、19、2、29、29、33 mm 2 。有限元模型中韧带使用的应力-应变曲线如 图2 所示。
有限元模型对真实脊柱进行了一定简化与假设。因此,需要对模型静力学行为仿真结果的可靠性进行验证。将2个椎体和1个椎间盘的组合视为1个脊柱功能单元。1个脊柱功能单元的生物力学特征反映了整个脊柱的特征 。胸腰段被认为是脊柱冲击损伤的常见部位 [21-22] 。取胸腰椎交界处T12～L1单元与文献[23-25]的研究进行对比验证。采用右手笛卡尔坐标系,选择上椎体几何中心为坐标系原点,约束椎体L1下表面6个方向上的自由度,T12上端面耦合至上表面任一点,分别在耦合点上沿 X 、 Y 、Z轴正负向以3个增量步施加7.5 N· m力矩模拟胸腰椎前屈后伸、轴向旋转、左右侧弯6种工况[见 图3 (a)]。
为验证模型动力学分析结果的可靠性,参照文献[26]对模型施加动态载荷并进行仿真分析(见 表2 )。将T2～L5脊柱模型分为上胸椎(T2～6)、下胸椎(T7～11)和腰椎(T12～L5)3段,每段脊柱上表面耦合14.7 kg的质量点,模拟人体上半身质量的影响。进行不同高度落体实验仿真。根据掉落高度将实验分为低能实验(A)与高能实验(B)两组。脊柱保持解剖学直立状态,每段脊柱下方设置1块骨性材料碰撞面,碰撞面具有一定转角,以消除棘突较长对碰撞的影响,模拟垂直冲击载荷。为减少计算成本,将脊柱与碰撞面的距离设置为10 mm,重力加速度为9.8 m/s 2 ,赋予脊柱一定初始速度,代替实验中从不同高度自由落体的速度[见 图3 (b)]。
由T12～L1模型在6种力矩载荷下的转角值与文献[23-25]的比较结果可见,在不同力矩下的载荷-转角曲线中,T12～L1均显示出转角随力矩增加而增大的趋势。左右侧弯时,模型最大转角值小于文献实验值,负号表示方向。轴向旋转、前后屈伸时,模型产生的最大转角值落于文献实验值之间(见 图4 )。在人体关节机制理论中,关节软骨之间通过润滑液进行液固相互作用来传递载荷,本文模型关节软骨连接处采用共节点形式进行约束绑定,一定程度上限制了脊柱的运动,这可能是分析数值结果小于文献数据的原因。本文认为,模型预测的6种力矩下载荷-转角曲线与文献[23-25]的研究结果吻合良好。
图4 T12～L1有限元模型在不同力矩下载荷-转角曲线与文献[23-25]实验数据对比
比较上、下胸椎与腰椎模型在3种掉落高度下轴向力峰值与文献[26]中测得平均峰值与标准差,结果显示,上胸椎在0.40、0.60 m高度下掉落受到的最大轴向力分别为1.7、1.8 kN,下胸椎在0.40、0.60 m高度下掉落受到的最大轴向力分别为1.3、1.8 kN,腰椎在0.46、0.82 m下落受到的最大轴向力分别为1.2、3.2 kN,与文献[26]实验获得的数据较为接近。在高能测试中,上、下胸椎、腰椎轴向力峰值分别为5.3、5.5、7.5 kN,大于文献[26]实验结果。3段脊椎模型仿真结果均处于文献数据误差范围内[见 图5 (a)]。
分析可知,3段脊椎模型轴向力峰值均随掉落高度的增加而增加,与文献数据体现出良好的吻合性。高能测试中,仿真获得的轴向力大于实验数据,而低能测试中未明显体现这一现象。上、下胸椎模型在低能测试中,最大轴向力不超过2 kN,而腰椎模型在0.82 m高度下落时轴向力大于2 kN。高能测试下,胸椎模型最大轴向力在5～6 kN范围内;腰椎模型则受到更高的轴向力,位于7～8 kN范围内。仿真数据在一定程度上说明腰椎部分在同等条件下倾向于承受更大的轴向力,在活动过程中更容易造成损伤。
对上、下胸椎与腰椎模型在3种掉落高度下矢状面弯矩与文献[26]实验测得平均弯矩和标准差进行比较,结果显示,上胸椎在0.40、0.60 m高度掉落下弯矩分别为21、23 N·m,未体现出与高度的相关性,在1.50 m高度弯矩为62 N·m,明显低于96 N·m的实验值,但仍处于实验误差内。下胸椎在0.40、0.60 m高度弯矩分别为9、15 N·m,在1.50 m高度下受到弯矩为37 N·m。低能测试中下胸椎几乎不承受弯矩;高能测试中弯矩在10～20 N·m范围内,但具有高达50 N·m的误差 。仿真获得的弯矩与掉落高度呈正相关趋势,结合上胸椎与腰椎模型分析,体现出建模的一致性。然而,文献[26]中下胸椎实验值与仿真结果具有一定偏差,考虑数值模型与实际人体具有一定差异,且在材料取值、各部位接触设置间存在一定简化,此偏差有待进一步探究。腰椎在0.46、0.82、1.50 m高度下落弯矩分别为18、22、21 N·m,均处于实验误差范围内[见 图5 (b)]。
同时,比较高能测试下模型上、下胸椎和腰椎承受的轴向力峰值、弯矩和文献[26]实验数据,结果显示,模型和文献[26]中腰椎在相同条件下承受的峰值力最高,弯矩最低。有限元模型的峰值力从上胸椎、下胸椎至腰椎逐渐增高;相反的,弯矩从上胸椎、下胸椎至腰椎依次递减。由文献[26]实验中峰值力呈现趋势与模型一致,弯矩在下胸椎处达到最小值,其次在腰椎,最后在上胸椎。
由高能实验中脊柱模型在冲击载荷下的应力分布可见,上胸椎最大应力发生在T4胸椎椎弓根处,高应力自T4椎弓根底部开始沿椎弓根延伸至T4上关节突,故该部位最易发生骨折与损伤。T5胸椎椎弓根发生应力仅次于T4胸椎,且椎体部分应力分散均匀。同时,高应力从T5椎弓根底部开始沿椎弓根扩散至T5椎弓根中段,有向上关节突扩散的趋势。T6胸椎椎弓根底部与顶部存在少量范围产生高应力,但未延展至椎弓根中部,椎体后部应力略大于前部。各胸椎椎弓根处应力均为该节胸椎的最大值。本文认为,高能冲击下上胸椎受到部分损伤,T4、T5、T6胸椎发生了一定程度骨折。下胸椎最大应力发生在T11椎弓根上缘与椎体连接处,各节胸椎最大应力同样集中在每节椎弓根处。本文结果与文献[26]实验结果一致。腰椎应力集中在椎骨、椎骨后部及椎弓根处。椎间盘内应力主要通过髓核进行传递,进一步扩散至纤维环,这与文献[25]实验结论一致[见 图6 (a)]。
由高能冲击下腰椎截面应力分布可见,应力首先发生在L5底部及各段椎体外沿,随着碰撞时间增加扩散至椎体内部。本文认为,皮质骨首先承担主要应力,侧面说明了脊柱材料采用灰度赋值的有效性[见 图6 (b)]。
比较有限元模型与文献[26]采用腰椎标本测试后CT扫描的轴向图像,结果显示,L1节腰椎主要应力集中在椎体后沿、椎弓根及椎板处,推测这些地方发生骨折与损伤,与文献[26]得出 “粉碎性L1椎体、双椎弓根和椎板骨折明显”的结论较为一致。
本文基于搭建的有限元模型,探究应用于跌落事故的人类脊柱这一区域的生理损伤机制。在动力学验证中,本文有限元模型CT数据来源与文献[26]采用的测试样本较为相近,模型数据源于志愿者年龄44～56岁、经过评估被认定为正常健康的脊柱。因此,两实验之间与年龄相关的影响认为可以忽略不计。本文重点考察胸腰椎模型有效性验证及其在动力学载荷中的应用。
在碰撞过程中,有限元脊柱模型前、后部均发生塑性应变,双侧椎弓根与椎体均存在骨折现象,表明载荷由整个椎骨分担 。由于前、后部均受到损伤,临床上认为此类发生机制存在不稳定性 [27-28] 。然而,高能测试下脊柱所受弯矩大于静力学验证中施加的弯矩,即正常屈伸所受弯矩,推测除主要轴向压缩载荷的传递之外,弯矩是脊柱发生损伤的次要机制。
腰椎在高应变率下的失效载荷强度比胸椎约高26% [29-30] 。因此,Yoganandan等 研究降低了从腰椎到胸椎的跌落高度,即从0.82 m下降至0.60 m。对于能量最低的跌落实验,将高度增加至0.40 m来充分体现胸椎碰撞性能。同时,选择较低附加质量以确保样本在低能量测试中不会受到损伤,从而可以评估亚失效和失效测试中同一样本的生物力学数据。
本文模型在尽量接近真实情况进行构建的基础上,最大程度还原Yoganandan等 实验环境与边界条件,由此获得的仿真数据更具可靠性与说服力。在对比实验中,模型与实验提取的均为脊柱在掉落过程中受到的轴向力峰值,不代表骨折时受到的轴向力。另外,有限元模型具有数字化、可重复利用、效率高的优点,避免了Yoganandan等 在实验中重复加载可能对测试样本造成损伤进而影响实验结果的问题。在同样材料参数与边界条件下,有限元模型可获得确定的结果,一定程度上规避了实验中因操作、材料等造成的误差。研究进一步表明,灰度赋值法在脊柱有限元模型材料赋值中能够发挥可靠应用。未来工作将利用已验证模型探究更复杂类型载荷,探讨人体胸腰椎在载人飞行、载人航天等载荷工况下损伤情况,并以此为基础进一步扩大模型,建立人体肌肉、内脏等软组织,为胸腰椎损伤的发生发展和治疗等提供更加全面准确的参考与理论依据。
本文建立并验证了正常成年男性T2～L5段脊柱三维有限元模型,研究脊柱以直立姿势自由掉落受到垂直载荷的动力学特性。对T12～L1段胸腰椎施加不同方向力矩,获取的转角-扭矩曲线与已发表文献数据相吻合。将脊柱分为上胸椎(T2～6)、下胸椎(T7～11)和腰椎(T12～L5)三部分在低能与高能冲击载荷下进行有限元仿真。结果显示,所有脊柱的轴向峰值力都随着掉落高度的增加而增加。在高能测试下,上、下胸椎和腰椎的矢状面弯矩随着掉落高度的降低而降低。本文认为,所建有限元模型具有良好的有效性。本文结果为人体在意外工况下受到生理损伤及其发生机制提供了研究方法和理论依据。
作者贡献声明： 黄莞沨负责研究设计和实施、资料收集及论文撰写;曲爱丽、李立、汪方负责资料收集及研究指导;王冬梅负责研究指导、论文指导及修改。
,
王冬梅 1 
(3)
微信
QQ空间
微博
QQ
