关于腰椎三维有限元模型的建模及有效性验证*

2015-06-01黄宇峰潘福敏赵卫东巴兆玉马敏吴德升

生物骨科材料与临床研究 2015年5期

黄宇峰潘福敏赵卫东巴兆玉马敏* 吴德升*

黄宇峰1潘福敏1赵卫东1巴兆玉1马敏2* 吴德升1*

目的利用三维有限元原理，通过三维有限元软件创建正常人腰骶椎三维有限元模型，并对模型进行有效性考证。方法选取一名健康成年男性志愿者，在平卧非承重条件下采用螺旋 CT行薄层扫描，利用数字软件Simple ware2.0、Geo magic studio9.0、Hypermesh10.0、ANS YS9.0及Abaqus6.10-l等建立正常人腰骶椎（L3-S1）模型。对健康完整模型进行四种状态（前屈、后伸、侧屈和轴向旋转）的加载分析，并将L3/4在不同工况下的角位移与以往文献结果比较，来验证模型的有效性。结果建立健康成年男性的L3-S1三维有限元模型，包括72158个节点，173035个单元；模型L3/4节段在前屈、后伸、侧屈及旋转运动时的角位移数值与其他实验数据对比相似性好，各椎间盘最大Von Mises应力与文献结果一致。结论通过有效性验证，可认为该基于正常人体L3-S1腰骶部的三维有限元模型在一定的条件下是有效的，可用于生物力学实验。

腰椎；有限元模型；生物力学

有限元方法是将一个连续的无限自由度问题简化为离散的有限自由度问题求解的一种数值分析方法。这一新兴的科学方法解决了许多用传统力学分析的数学工具无法求解的骨科学问题[1]。特别是对脊柱力学研究方面，有着独特的优越性，并能提供全域性信息。而腰椎的各个组成部分表面形态极不规则，椎间盘的结构复杂，特别是胶原纤维走向引起的几何非线性效应的改变，这给有限元模型的建立造成了一定的难度。当研究下腰椎的多节段运动时，非线性几何效应就变得非常重要，由于模型几何结构的原因导致的椎体偏移也将相当明显。因此腰骶椎的三维有限元模型包括非线性接触条件、非线性材料属性以及非线性集合效应。直到近年来，随着计算机技术、软件以及有限元理论的高速发展，针对腰骶椎模型的有限元分析才逐渐具可行性。

本研究通过数字化软件建立正常成年男性下腰椎三维有限元模型（L3-S1），并对模型有效性进行验证。为腰骶椎的生物力学研究提供有效的方法和途径。

1 实验对象及设备

1.1 实验对象

选取一名健康成年男性志愿者，30岁，身高175cm，体重75kg。否认腰椎疾病及外伤史，X线、CT及临床检查排除脊柱损伤、肿瘤转移、畸形等病理变化。志愿者自愿参与本研究，并签署知情同意书。

1.2 设备及相关软件

飞利浦Brilliance 64排多层螺旋CT，连续扫描数据采集，采集层厚0.625mm。Simpleware2.0用于三角网格模型的建立。GeomagieStudio 9.0用于三维几何实体模型的建立。ANSYS9.0用于三维网格划分与构建有限元模型。HyPer-Mesh10.0用于有限元网格的划分。Abaqus6.10-l用于有限元模型的建立和力学分析。上述数字化软件由上海大学力学实验中心提供。

2 实验方法

2.1 成年男性下腰椎三维有限元模型的建模

2.1.1 CT扫描体位及扫描范围

患者自然仰卧位，正交定位线置于身体中线，受检区下置校正躯体模块。从第三腰椎椎体终板上缘沿人体躯干长轴以层厚0.625mm向远端螺旋扫描直到第一骶椎下缘。共得到二维图像数据400幅，实际建模过程共采用其中间有效部位269幅。所有图像数据均采用DICOM文件形式存盘。

2.1.2 正常人下腰椎三维网格模型的建立

将DICOM格式的CT扫描二维图像数据导入有限元建模软件simpleware2.0中，使用ScanlP模块对原始图像进行重新采样，并将分辨率调整至合适水平，避免影响最后建立的实体模型的精细程度。再经过图像切割、降噪等处理后得到一个合适的兴趣区域。最后运用填充工具及闭值工具根据CT灰度值对该重要区域进行相关阈值设定，提取指定灰度范围内的数据从而标识、分割出腰骶椎中各块特征性的骨骼并分别建立蒙皮。计算蒙皮生成腰骶椎（L3-S1）的点云模型，并将该三角网格模型以STL格式保存。

2.1.3 正常人下腰椎三维实体模型的建立

将制作三角网格模型中采集生成的点云文件导入逆向工程软件GeomagicStudio9.0中进行修补与优化。针对模型采用了点、多边形和成形三个处理阶段。先根据曲率连续性原则处理模型表面空洞，使模型表面更为松弛平滑，防止非特征性的高曲率以及自相交面的产生。随后按照模型表面曲率变化情况描绘轮廓线，通过对每块封闭轮廓线（嵌片）进行参数设置来调整模型表面面片的精细程度，最后生成G1连续的非均匀有理B样条曲面（NURBS），保存为IGES格式。

2.1.4 定义模型单元及网格划分

将之前建模过程中获得的正常人腰骶椎椎体、椎间盘的三维实体模型数据导入HyperMesh10.0软件，使用其前处理工具调整修改各部分网格的分割密度，以此来划分并修改网格。我们采用6节点五面体（C3D6）结构定义骨性终板与椎体后结构，椎体松质骨模拟为10节点四面体单元（C3D10），椎间盘纤维环及髓核定义为三维8结点杂交单元（C3D8H），韧带、关节囊模拟为仅受拉杆单元（Truss）。

2.1.5 定义材料属性完成建模

数据在完成网格划分后被导入有限元软件ANSYS9.0中进行赋值。腰骶椎椎体的骨性后结构包括椎板、棘突、横突以及小关节等，每种骨性结构均包含皮质骨与松质骨，但其体积较椎体而言较小不会影响计算结果，因此在本研究中，我们对于后结构的两种骨质未予以区分，而是采用同一种材料系数，提高了计算机响应速度。建模单元、材料参数和韧带的横截面积如表1所示，所有参数均来源于文献报道[2-5]。

表1 各结构建模单元、材料属性参数与韧带的横截面积

2.2 正常人下腰椎三维有限元模型的有效性验证

有限元分析的核心是模拟，所建立的有限元模型要尽可能地与真实结构相近。因此，必须对其真实性、精确性进行检验。在本研究中，我们选取正常腰椎模型中L3-L4-L5三个节段，边界条件设定为在模型底部固定，最近端的 L3可自由活动。在模型L3上终板施加垂直向下的500N预载荷（约为该志愿者体重的2/3），附加扭矩为10N.m，模拟腰椎前屈、后伸、左右侧屈和左右旋转四种运动。使用Abaqus6.10-1的job模块，对腰椎中立、前屈、后伸、旋转、侧弯五种工况进行求解，求解结束，在visualization模块中选中各椎间盘，查询并记录各椎间盘最大Von Mises应力。

同时，记录模型在各个运动方向上的最大角位移作为下腰椎活动范围（Range of motion,ROM）。具体方法为在远端椎体相对固定的前提下，在上一相对活动的椎体上表面选择2个节点，将节点空间坐标的连线在加载负荷前后形成的夹角记录为，即最大角位移。根据相关文献[6]，我们将局部坐标系原点定义为远端椎体的终板中心，平行于终板在矢状位上为X轴，垂直于终板方向为Z轴方向，根据右手原则，将平行于终板向右为Y轴方向，如图1（彩图见插页）。

图1 局部坐标系示意图。原点为远端椎体的终板中心，X轴为平行于终板在矢状位方向上，Z轴为垂直于终板方向，根据右手原则，Y轴为平行于终板向右方向

2.3 实验结果

2.3.1 建立了正常人下腰椎三维有限元模型

本研究根据真实健康人体脊柱腰骶椎 CT资料，使用Simpleware2.0、GeomagicStudio9.0、Ansys9.0、HyperMesh 10.0及Abaqus6.10-1等有限元分析软件建模并优化L3-S1节段，整个模型共包括72158个节点，173035个单元，得到几何外观与实际标本高度相似的模型，如图2-5（彩图见插页）。

2.3.2不同工况椎间盘最大Von Mises应力的测量结果

对健康完整模型进行四种状态（前屈、后伸、侧屈和轴向旋转）的加载分析，测量本模型L3/4节段椎间盘在四种不同工况下的最大VonMises应力，并与其他以往的研究结果比较来验证模型的有效性。结果显示，垂直加载状态下，椎间盘应力集中于后半部分，靠近椎弓根处压应力最大，由后向前逐渐减小，椎间盘前缘处的压应力为最小值；前屈运动时，椎间盘前缘区应力较大，且为压应力。后伸位时与之相反，椎间盘后侧为压应力，前侧为拉力；侧屈运动时椎间盘压缩区应力较大且为压应力，至对侧逐渐减小且变为拉力。本实验数据与以往其他报道[7-9]的数据相比略小，尤其在后伸运动中，但总体而言与以往其他试验数据相似性好，表明本模型在一定条件下有效性高，如表2所示。

表2 不同研究间L3/4椎间盘最大应力比较（MPa）

图2 椎间盘有限元模型图3椎间盘实体模型剖视图

图4 正常腰骶椎模型正面及背面观

图5 正常腰骶椎模型

表3 不同工况下模型的活动角度

2.3.3 不同工况下模型活动角度的测量结果

本研究对正常人腰椎模型进行四种运动状态（前屈、后伸、侧屈和轴向旋转）的加载分析，边界条件设定为在模型底部固定，最近端的L3可自由活动。在模型L3上终板施加垂直向下的500N预载荷（约为该志愿者体重的2/3），附加扭矩10N.m，模拟腰椎前屈、后伸、左右侧屈和左右旋转四种运动，使用有限元分析软件 Ansys9.0测得其活动角度与旋转刚度如表3所示，其中侧屈活动角度取左侧屈与右侧屈均值，旋转运动角位移采用旋左与旋右状态下的均值。经过加载运算后得到四种不同工况下的位移分布云图，如图 6（彩图见插页）。将结果与其他以往的研究结果[10-13]比较后，本模型在各种工况下的旋转刚度介于不同研究结果之间，可认为模型有效，结果如表4所示。

图6 正常腰骶椎模型位移云图。左上图为在前屈状态下位移云图；右上图为在后伸状态下位移云图；左下图为在侧屈状态下位移云图；右下图为在旋转状态下位移云图

表4 不同研究间腰椎有限元模型平均旋转刚度的比较（N.m/°）

3 讨论

利用三维有限元法建立生物力学模型具有其独特的优点。其操作简便，建模灵活多变，分析研究的可重复性高，在反复利用同一个有限元模型过程中不会造成模型生物力学属性的改变。同时又避免了在体外尸体实验中难以获取实验标本、费用昂贵等缺点。尤其是能较好地观察外力作用下脊柱内部结构的应力变化情况，这是尸体实验、动物实验等其他传统实验方法所不能比拟的[14]。

现如今三维有限元分析目前虽然已广泛应用于腰椎的生物力学之中，但由于腰椎几何形态高度不规则，不同结构之间材料系数差异较大，往往需要对腰椎的有限元模型进行一定程度的简化与假设才能使之应用于实验研究中。因此，当任何一个三维有限元模型建立之后，对其进行应力分析之前都应当验证模型的有效性，只有在确认有限元模型准确、有效的前提下，其后的分析才会有意义[15]。在目前的研究中应用较多的有以下几种有效性验证方法：自身验证，包括模型几何形态的精确程度与材料参数的准确性；通过与其他三维有限元分析实验结果比对验证，这是目前应用最多的一类验证方法；与尸体模型的生物力学实验结果或者体内实验结果对照，对三维有限元模型施加与尸体模型或实验相同的载荷或尽可能模拟体内实验载荷，如果二者的活动范围相符，则该有限元模型可被认为具有一定的有效性。

本研究通过灵活运用多种三维有限元软件根据正常人体腰骶椎CT扫描资料完成了L3-S1三维有限元模型的造模。由于腰骶椎表面形态高度不规则，并且骨性后结构包括椎板、棘突、横突以及小关节等，每种骨性结构均包含皮质骨与松质骨两种不同密度的材料，但其体积较椎体而言较小不会影响计算结果，为了提高计算机响应速度及简化计算过程，我们并未区分后结构的两种骨质，而是设置为同一种材料系数。

在本研究中，我们根据真实的健康人体脊柱腰骶椎CT资料，使用多个有限元分析软件建模并优化L3-S1节段，整个模型共包括72158个节点，173035个单元，丰富的节点与单元数目保证了模型几何形态的精确性。同时，所涉及到的材料参数与韧带的横截面积均来源于被较多研究者引用的文献，准确性高。最终得到了几何外观与实际标本高度相似的模型。由于不同研究可能对模型施加的载荷并不相同，如果单纯地比较绕轴转角就不够精确，因此本研究在与其他实验对比验证活动度时采用了平均旋转刚度来验证模型的有效性，平均刚度是力矩与载荷形成的最大角位移之比，单位为N.m/°。我们将结果与其他以往的实验结果比较后发现本模型在前屈、后伸、侧屈以及旋转四种工况下FSU的平均旋转刚度介于不同的研究结果之间。此外，我们还测量了本模型L3/4节段椎间盘在四种不同工况下的最大VonMises应力，并与其他以往的研究结果进行比较，发现本实验数据与以往其他报道的数据相比略小，尤其在后伸运动中，但总体而言相似性好。

综上所述，本研究成功建立了正常人下腰椎三维有限元模型，该模型在一定条件下具有较高的有效性，能为脊柱的生物力学研究提供一种简便、科学、独特的研究方法。相信在不久的将来，采用计算机技术辅助的三维有限元分析进行脊柱生物力学研究，将成为解剖学、生物力学和临床医学工作者的一个重要的研究领域[16]。

[1] 李伟，陆皓，孙康，等.复合材料与金属材料髋关节假体应力分布的三维有限元分析[J].生物骨科材料与临床研究，2005，2（3）：1-4.

[2] Kim K,Lee SK,Kim YH.The biomechanical effects of variation in the maximum forces exerted by trunk muscles on the joint forces and moments in the lumbar spine:a finite element analysis[J]. Proc Inst Mech Eng H,2010,224(10):1165-1174.

[3] Xu HG,Yang XM,WuTL,et al.Finiteelementanalysisof the screw in percutaneous axial lumbosacral interbody fusion[J].Orthop Surg,2010,8,2(3):207-210.

[4] Goto K,Tajima N,Chosa E,et al.Effect of lumbar spinal fusion on theother lumbar intervertebral levels（three-dimensional finite element analysis）[J].J Orthop Sci,2005,8(4):577-584.

[5] Chen CS,Chen CK,Liu CL.A biomechanical comparison ofposterolateralfusion and posteriorfusion in the lumbarspine[J].J Spinal Disorders Tech,2002,15(1):53-63.

[6] Rousseau MA,Bradford DS,Bertagnoli R,et al.Disc arthroplasty design influences intervertebral kinematics and facet forces[J]. Spine J,2006,6(3):258-266.

[7] 陈肇辉.腰椎棘突间撑开装置的三维有限元分析及临床应用研究[PhD].第二军医大学，2010.

[8] Wang X,Dumas GA.Evaluation of effects of selected factors On intervertebral fusion：a simulation study[J].Med Eng Phys,2005, 27(3):197-207.

[9] El RiehM,Arnoux P,Wagnac E,etal.Finiteelementinvestigation of the loading rate effect on the spinal load-sharing changes under impact conditions[J].J Biomech,2009,42（9）:1252-1262.

[10]HethJA,HitchonPW,Goel VK,etal.A biomechanical comparison between anterior and transverse interbody fusion cages[J].Spine, 2001,15,26（12）:261-267.

[11]闫家智，林欣，王磊，等.计算机辅助的腰椎单节段融合后邻近节段的力学变化研究[J].徐州医学院学报，2011，31（9）:573-576.

[12]张德胜，宋岳明.L3-L5三维非线性有限元模型的建立及临床意义[J].生物医学工程学杂志，2006，23（6）:1250-1252.

[13]Vadapallis,SairyoK,GoelVK,et al.Biomechanical rationale for using Polyetheretherketone spacers for lumbar inter body fusion-A finite element study[J].Spine,2006,31（26）:992-998.

[14]Daniel G.T.Strange,Sandie T.Fisher,Philip C.Boughton,et al. Restoration of compressive loading properties of lumbar discs with a nucleus implant.A finite element analysis study[J].Spine, 2010,10(7):602-609.

[15]龚志强，姜晓幸，冯振洲，等.腰椎单双侧椎弓根螺钉固定后的有限元分析比较[J].复旦学报（医学版），2013，40（2）:187-192.

[16]苏春涛，眭承志.腰椎三维有限元生物力学分析的技术应用进展[J].中国康复医学杂志，2013，28（5）:479-482.

Establishment of normal lumbosacral vertebral three-dimensional finite element

Objective Establish normal lumbosacral vertebral three-dimensional finite element model by using three-dimensionalfinite element software,then ver ify thevalidityofthemodel.Methods Choose aheal thy maleadult volunteer, giving a spiral CT scanning under the condition of nonbearing supine posture,digital software including Simpleware2.0, Geomagic studio9.0,Hypermesh10.0,ANSYS9.0 and Abaqus6.10-l is used to establish normal lumbosacral vertebral（L3-S1）model.And to verify the validity of the model.Results Normal L3-S1three-dimensional finite element model is established,including 72158 nodes,173035 units;the angular displacement value of Model L3/4segment in forward bends,rear protraction,lateral flexion and rotation motion is similar when compared with other experimental data,the maximum Von Mises stress of intervertebral disc is consistent with the literature results.Conclusion With the normal L3-S1lumbosacral vertebral three-dimensional finite element model validation,under certain conditions,this model can be considered efficacious to apply in biomechanical experiments.

Lumbar spine;Finite element model;Biomechanics

R318

10.3969/j.issn.1672-5972.2015.05.005

swgk2015-05-00105