以三维非线性有限元法分析融合腰椎的生物力学变化

doi:10.3969/j.issn.2095-4344.2013.35.009

摘要/Abstract

摘要：

背景：MATLAB具有大型数值计算、数学绘图和简单有限元分析的能力，建模速度快，能够识别BMP和JPG格式的灰度图，并可将识别的数据直接转换为ANSYS三维有限元软件可识别的格式，避免了常规在Autocad软件中图像的重新定位和二次处理而产生的人为误差。
目的：寻求一种简便、快捷、准确的方法建立各种腰椎融合模型，对腰椎融合后腰椎活动节段的生物力学进行分析。
方法：利用薄层CT技术，结合Matlab(Matrix Laboratory即矩阵实验室)科学计算软件，辅助Ansys有限元分析软件建立各种腰椎融合模型。对所建立的模型施加各种载荷，分析其生物力学变化。
结果与结论：对所建立的融合模型施加轴向、屈曲及伸展载荷后，生物力学分析得出所有的融合模型中椎体间融合稳定性最好。和关节间融合相结合，椎体间，侧后方，后方融合模型的轴向位移较单独的椎体间，侧后方和后方融合模型分别减少5%、1%和4%。在伸展-屈曲载荷条件下旋转角度分别减少23%、11%和45%。应力向融合部位集中的现象表明融合块可增加载荷的转移。说明薄层CT和Matlab软件的辅助使得Ansys有限元软件建立腰椎融合模型的速度和精度提高。小关节融合和椎体间，侧后方，后方融合相结合，可使腰椎获得更好的稳定性，这种稳定性的增加在后方融合中更为显著。后方融合模型的应力分布更加合理。

关键词: 骨关节植入物, 脊柱损伤基础实验, 生物力学, 有限元方法, 腰椎融合, 小关节融合, 融合模型, 载荷, 矩阵实验室

Abstract:

BACKGROUND: MATLAB has capabilities of large numerical calculation, mathematical drawing and simple finite element analysis. It can establish models rapidly and can be able to identify the grayscale with BMP and JPG format, and it can directly transform the identified data into ANSYS finite element software-readable format, thus avoiding personal error produced by the repositioning and secondary treatment in the Autocad software.
OBJECTIVE: To find a simple, convenient and accurate method to construct the model of lumber fusion and to analyze the biomechanics of lumbar motion segment after lumbar fusion.
METHODS: Lamellar CT and Matlab (Matrix Laboratory) scientific computing software combining Ansys finite element software was used to construct the models of lumber fusion. Then the models were loaded to analyze the biomechanical change of the fusion model.
RESULTS AND CONCLUSION: The established models were loaded with axial, bending and stretching loads, and the biomechanical analysis showed that interbody fusion had the best stability among all the fusion models. Combined with joint fusion, the axial displacements of interbody, rear side and rear fusion models were decreased by 5%, 1% and 4% than that of simple interbody fusion, posterolateral fusion and posterior fusion models. Under the stretch-buckling load conditions, the rotation angles were reduced by 23%, 11% and 45%. Stress concentration to the fusion parts showed fusion block could increase the load displacement. The technology of lamellar CT, Matlab software and Ansys finite element software can accelerate the construction of lumber fusion model and make the model more accurate. Facet joint fusion combined with interbody fusion, posterolateral fusion and posterior fusion can get better lumbar stability, and this increased stability is more significant in the rear fusion. Stress distribution of posterior fusion is more reasonable.

Key words: bone and joint implants, basic experiment of spinal cord injury, bone tissue engineering, biomechanics, finite element method, lumbar fusion, facet joint fusion, fusion model, loading, matrix laboratory

中图分类号:

R318

张延辉，姜宏春，李静，张宝弟. 以三维非线性有限元法分析融合腰椎的生物力学变化[J]. 中国组织工程研究, 2013, 17(35): 6273-6280.

Zhang Yan-hui, Jiang Hong-chun, Li Jing, Zhang Bao-di. Biomechanical changes of lumbar segment after fusion analyzed with three-dimensional nonlinear finite element[J]. Chinese Journal of Tissue Engineering Research, 2013, 17(35): 6273-6280.

图/表（结果） 6

2.1 成功建立模型正常L₄_-₅活动节段完整有限元模型共包括47 397个实体单元，598个缆单元，585个壳单元，384个目标单元，232个接触单元，总计共67 259个结点。各部位的单元类型，单元数和材料属性见表1。

2.2 模型的验证 轴向压缩载荷条件下，正常L₄-L₅节段有限元模型的预测结果与其他模型的研究结果相似^[11]。正常L₄-L₅有限元模型的压缩-轴向位移曲线位于 Brown^[9]、Markolf^[10-11]、Virgin等^[12]的实验研究曲线之间，见图5，压缩载荷增加，L₄-L₅节段椎间盘刚度亦轻度增加；压缩-椎间盘内压力曲线位于Nachesmon^[13]与 Rolander^[14]的试验研究曲线之间，见图6，且呈线性走行。

2.3 融合模型施加载荷后的变化

轴向压缩的变化：每个模型由轴向压缩引起的L₄椎体上表面的轴向位移和屈曲旋转见图7。

轴向位移：轴向压缩载荷下椎体间融合合并小关节融合的轴向位移最小，而后按照轴向位移从小到大的顺序模型依次是椎体间融合无小关节融合、侧后方融合合并小关节融合、侧后方融合无小关节融合、后方融合合并小关节融合、小关节融合、后方融合无小关节融合和正常模型。合并小关节融合的模型较无小关节融合的同样模型轴向位移小，稳定性好。椎体间融合合并小关节融合、侧后方融合合并小关节融合和后方融合合并小关节融合的稳定性较椎体间融合无小关节融合、侧后方融合无小关节融合和后方融合无小关节融合分别增加5%、1%和4%。这些数据表明小关节融合可以增加稳定性。

伴随的屈曲旋转：轴向压缩载荷下，正常模型和椎体间融合模型都很少有屈曲旋转，但在其他模型中，可以看到屈曲旋转。轴向压缩载荷下小关节融合模型的屈曲旋转角度最大，其他模型从大到小依次是后方融合合并小关节融合、侧后方融合合并小关节融合、侧后方融合无小关节融合、和后方融合无小关节融合。

屈曲和伸展载荷下的运动：椎体间融合合并小关节融合屈曲和伸展方向的稳定性最好，其他模型稳定性由大到小依次是椎体间融合无小关节融合、侧后方融合合并小关节融合、后方融合合并小关节融合、小关节融合、侧后方融合无小关节融合和正常模型，见图8。后方融合合并小关节融合、小关节融合和侧后方融合无小关节融合的旋转角度几乎相同。与椎体间融合无小关节融合、侧后方融合无小关节融合、后方融合无小关节融合相比，椎体间融合合并小关节融合、侧后方融合合并小关节融合、后方融合合并小关节融合屈曲和伸展方面的稳定性分别增加了23%、11%和45%。以上数据表明小关节融合可增加在屈曲和伸展方面的稳定性，这种稳定性的增加在后方融合模型中最为明显。

Von Mises应力分布：见图9-11。轴向压缩载荷条件下Von Mises应力分布，椎体间融合无小关节融合模型的应力集中在椎体，见图9。

侧后方融合合并小关节融合模型的应力集中于椎体的前部和侧后方融合的骨棒，见图10。后方融合合并小关节融合模型的应力在椎体的前部和终板比较集中。轴向载荷下，后方融合的骨棒所承受的应力较小。

屈曲和伸展载荷条件下的Von Mises应力分布几乎相同。屈曲和伸展载荷条件下，椎体间融合无小关节融合模型的应力集中于椎体的前部和后部，见图9。侧后方融合合并小关节融合模型的应力集中在上方椎体，侧后方融合的骨棒和下方椎体，见图10。后方融合合并小关节融合模型中上方椎体，终板和棘突是应力集中区。见图11。

参考文献

[1] Goel VK, Kong W, Han JS, et al.A combined finite element and optimization investigation of lumbar spine mechanics with and without muscles. Spine (Phila Pa 1976). 1993; 1;18(11):1531-1541.

[2] Rybicki EF, Simonen FA, Weis EB Jr. On the mathematical analysis of stress in the human femur.J Biomech.1972; 5: 203.

[3] Huiskes R，Chao EY. Asurvey of Finite element analysis in orchopedic biomechanicsi the first fecade. J Biomech. 1983;16(6):385-409.

[4] Liu YK, Ray G, Hirsch C. The resistance of the lumbar spine to direct shear. Orthop Clin North Am. 1975;6(1):33-49.

[5] Chen CS,Cheng CK,Liu CL,et al.Stress analysis of the disc adjacent to interbody fusion in lumbar spine. Med Eng Phys.2001; 23:483-491.

[6] Yoganandan N, Kumaresan S, Voo L,et al. Finite element of the human lower cervical spine. Parametric analysis of the C4 -C6 unit. J Biomech Eng. 1997;119(1):87-92.

[7] Wu JS , Chen JH.Clarification of th mechanial behaviour of spinal motion segments through a three dimensional poroe2lastic mixed finite element model. Med Eng Phys. 1996 ;18(3):215-224.

[8] Gilbertson LG, Goel VK, Kong W, et al. Fininte element methods in spine biomenanics research Crilt. Crit Rev Biomed Eng.1995;23(5-6):411-473.

[9] Brown T, Hansten RJ, Yorra AJ. Some mechanical tests on the lumbosacral spine with particular reference to intervertebral discs. J Bone Joint Surg Am.1957; 39(5):1135 - 1164.

[10] Markolf KL. Deformation of the thoracolumbar intervertebral joints in response to external load : a biomechanical study using autopsy material. J Bone Joint Surg Am.1972;54(3) : 511- 533.

[11] Markolf KL, Morris JM. The structural components of the intervertebral disc. A study of their contributions to the ability of the disc to withstand compressive forces. J Bone Joint Surg Am.1974;56(4):675-687.

[12] Virgin WJ. Experimental investigations into the physical properties of the intervertebral disc. J Bone Joint Surg Br.1951;33(4):607 - 611.

[13] Nachemson A.Mechanical stresses on lumbar disks. Curr Pract Orthop Surg. 1966;3:208-224.

[14] Rolander SD. Motion of the lumbar spine with special reference to the stabilizing effect of posterior fusion. An experimental study on autopsy specimens. Acta Orthop Scand. 1966:Suppl 90:1-144.

[15] Saito T ,Tamatouro T ,Shikara J ,et al. Analysis and prevention of spinal column deformity following cervical laminectomy , pathogenic analysis of post laminectomy deformities. Spine, 1991;16 :494-502.

[16] Kleinberger M. Application of finite element techniques to the study of cervical spine mechanics. In : Proceedings of the 37th Steppe Car Crash Conference. San Antonio Tesas,1993; 7-8.

[17] Bozic KJ, Keyak JH, Skinner HB, et al. Three dimensional finite element modeling of a cervical vertebra : An inverstigation of burst fracture mechanism. J Spinal Disorders.1994; 7:102-110.

[18] Voo L,Denman J, Kumaresan S ,et al. Development of 3-D finite element model of cervical spine. Adv Bioengineer Development of finite element model of cervical spine. Avd Bioeng.1995;31:13-14.

[19] Kumaresan S,Yoganadan N,Pintar FA,et al. Finite element model of cervical laminectomy with graded facetectomy. Spinal Disorder.1997; 10 (1) :40 - 46.

[20] Kumaresan S, Yoganandan N, Pintar FA. Finite element analysis of anterior cervical spine interbody fusion. Biomed Mater Eng.1997; 7:221-230.

[21] 刘耀升,陈其昕,廖胜辉,等.椎问盘高度降低及退变对腰椎生物力学影响的有限元分析[J].中国临床解剖学杂志,2006,24(5): 566-570.

[22] Lee CK,Langrana NA.Lumbosacral spinal fusion.A biomechanical study.Spine.1984;9:574-581.

引言

有限元方法是生物力学试验研究方法的必要补充和有效的替代工具，它不仅能较好地模拟复杂的力学环境，而且可以提供其他实验方法不能直接测量的结构内部的力学反应，从而获得全域性信息。MATLAB科学计算软件辅助建模，为腰椎有限元模型的建立开辟了一条新途径。

MATLAB全称matrix laboratory，原意为矩阵实验室，具有大型数值计算、数学绘图和简单有限元分析的能力，建模速度快，能够识别BMP和JPG格式的灰度图，并可将识别的数据直接转换为ANSYS三维有限元软件可识别的格式，避免了常规在Autocad软件中图像的重新定位和二次处理而产生的人为误差。

ANSYS是目前世界顶端的有限元商业应用程序，是融结构、流体、电场、磁场、声场分析于一体的大型通用有限元分析软件。该应用程序的主要特点是它能与多数CAD软件接口，实现数据的共享和交换，具有用户界面好，前、后处理和图形功能完备，并且使用高效的有限元系统。能够高效地求解各类结构的静力、动力、振动、线性和非线性问题；能分析系统受到外力负载所出现的反应，根据该反应可知道系统受到外力负载后的状态，进而进行应力分析。

腰椎融合可为腰段脊柱提供良好的稳定性。融合方式根据融合的位置不同分为椎体间融合，后方融合，侧后方融合及小关节融合。人们已经利用尸体的脊柱进行了腰椎融合生物力学分析。Goel等^[1]于1988年描述了利用椎体间骨移植物和钢板螺钉进行脊柱操作的机械力学。尔后，有人报道用有限元方法进行脊柱操作^[2-4]，但很少有人对腰椎融合进行报道。Chen等^[5]报道了相邻节段椎体间融合的应力分析。接着他们报道了侧后方融合和后方融合的对比，但是没有涉及到小关节融合。然而基于三维有限元分析的椎体间融合，侧后方融合和后方融合的比较研究是缺乏的。实验对不同的融合方式(包括椎体间融合，侧后方融合，后方融合及小关节融合)进行研究，将小关节融合和其他3种融合方式相结合，对其进行验证和各种不同载荷条件下生物力学的分析。

材料方法

设计：单一样本观察。

时间及地点：实验于2011年9月至2012年10月在牡丹江医学院完成。

对象：选取1名健康成年男性志愿者作为测试对象，对其腰椎进行X射线扫描以排除腰椎的异常。

方法：

扫描数据：测试对象取仰卧位，使用美国GE Lightspeed16层螺旋CT机。扫描条件：选择L₄-L₅骨组织窗扫描，层厚1 mm，床进速度1 mm/s。

图像预处理：将获得的图像储存在与CT机相连的电脑中，见图1，转移到移动硬盘中带回。以jpg格式的文件储存在工作机中，L₄-L₅节段CT图像共48张。将其转移到photoshop图像处理软件中进行图像的处理，把jpg格式文件转换为bmp格式的灰度图，尺寸为512×512，见图2。由于从CT获得的图像直接从电脑上拷贝，不涉及到扫描的步骤，所以可省略在photoshop软件中图像精确对位的步骤。

重建数据的采集和数据转换：利用Matlab软件中的imread函数读入现有的48幅L₄-L₅节段CT图像。用 mat2gray函数将图像的类型由uint8转换为double.而后用cat函数将转化后的每一幅图像集中到一个三维体数据集中，并将此数据集保存为mat文件。mat 文件是一个二进制的、与平台无关的数据文件，Matlab使用其存放工作空间中的各种变量和矩阵、数组。把三维数据集保存为mat文件的好处是：由于所用到的图片多达48 幅，数据量很大，如果每次运行程序时都重新读入这些图片将很费时。若将其保存为mat文件只需在下次运行时，用load 命令将其载入即可，以达到节省时间之目的。

三维模型的构建：将三维数据集存储成ANSYS可以直接识别的命令格式，存储形式为*.dat或*.txt，其中关键点的存储形式为k、M、x、y和z，其中M为关键点的参考序号。在ANSYS有限元软件中读取MATLAB转换的文本文件，采用自底向上的方法，即从最低级图元向上构造模型，首先定义关键点，然后依次是相关的线、面、体。得到L₄-L₅活动节段完整有限元模型，利用有限元软件的自由网格划分功能进行网格划分，见图3。

有限元模型的单元划分：L₄_-₅活动节段完整有限元模型共包括47 397个实体单元，598个缆单元，585个壳单元，384个目标单元，232个接触单元，总计共67 259个结点。模型包含了几何非线性、材料非线性、接触非线性3种结构非线性类型。其中椎体皮质骨用6结点三角面薄壳单元模拟，单元厚度设定为1.0 mm，椎体松质骨、终板、椎间盘纤维环、基质以及后部结构，包括椎弓根、峡部、终板、棘突、横突、小关节突以及关节软骨均由10结点四面体单元模拟。

椎间盘和髓核被模拟为不可压缩的黏弹性液体，其体积模量为1 666.70 MPa，髓核的体积占椎间盘体积的48.87%，椎间盘纤维环由基质与埋在基质中的胶原纤维所构成的混合体表示，前者由3层连续性环状体构成，纤维环纤维由8层只承受张应力的索单元构建，纤维在环状体中呈剪刀状方式走行，并与椎间盘平面成平均±30°的夹角^[39]。纤维体积占纤维环总体积的19%，各层纤维的横截面积自外向内逐渐减小，纤维环纤维的弹塑性应力应变关系曲线参考文献而设定^[6-7]。

椎间盘的上下表面由1.0 mm厚的软骨终板构成，L₄下终板与L₅上终板凹陷角分别设为156.9°和165.6°。关节突软骨层厚0.5 mm，关节软骨层的表面接触选用面-面接触单元模拟(无摩擦的滑动表面接触单元)，接触单元的初始间距为0.6 mm，摩擦因数为0。关节软骨的材料属性被处理为线弹性，其弹性模量设定10 000 MPa^[8-10]。有限元模型包含的前纵，后纵，棘上，棘间，横突间韧带以及黄韧带和关节囊韧带均采用超弹性索单元模拟，单元的位置、横截面积及长度决定于各组韧带的实际解剖结构。见表1。

椎体间融合模型：将纤维环的最外层材料变为皮质骨并把内层的纤维环和髓核变为松质骨，通过这样的转换，标准模型成为椎体间融合的模型。

侧后方融合和后方融合的模型：通过向正常模型添加骨棒建立L₄-L₅侧后方融合和后方融合的模型，见图4。骨棒由表面的皮质骨和内部的松质骨构成。

侧后方融合模型骨棒的形状是长方体，大小为 10 mm×10 mm×51 mm，由372个皮质骨单元和228个松质骨单元构成，共计600个单元，见图4A。后方融合模型的骨棒大小和形状与侧后方融合模型的骨棒相同，由278个皮质骨单元，147个松质骨单元构成，共计425个单元组成，见图4B。

侧后方融合模型的固定位置为从L₄椎体的横突到L₅椎体的横突和小关节的侧方。这和临床上的应用是一样的。后方融合模型的固定位置是从L₄椎体的棘突到L₅椎体的棘突。

小关节融合的模型：为了模拟小关节的融合，必须将正常模型小关节空白单元区用皮质骨单元充填。

边界和载荷条件：L₅椎体的下表面和下方的小关节固定不动。实验中的载荷条件是轴向压缩和屈伸活动。在L₄椎体表面的中央分十步逐渐增加轴向载荷至1 000 N。在L₄椎体表面分15步逐渐增加屈曲-伸展方向载荷至15 N·m。每个模型都施加轴向压缩载荷，屈曲和伸展载荷，并对结果进行比较和应力分布的分析。

主要观察指标：①模型验证结果。②融合模型施加载荷后的变化(包括轴向压缩的变化、屈曲和伸展载荷下的运动、Von Mises应力分布)。

讨论

利用三维有限元方法进行腰椎模型构建的方法很多^[15-20]，Chen等^[5]首先用有限元方法建立了多节段腰椎活动模型(L₁-L₅)，但建立的过程繁琐，步骤较多，而且对后方结构和椎间盘的材料特性的设定过于简单。因此所建立的模型精度有限。目前比较先进的构建方法是刘耀升等^[21]采用改良的“非种子区域分割方法”提取腰椎CT图像数据中目标区域得到二值图像，用Marching Cubes方法由二值数据生成初始表面模型。采用反映腰椎生理曲度的“最佳切割平面”从初始表面模型获得非平行的切割轮廓线并建立“分段线性子空间”，后者经仿射变换到“规则子空间”重构腰椎曲面，最后逆变换恢复腰椎原三维空间形状特征。此方法所建立的腰椎模型精确高，能较好的反映人体脊柱的真实情况，但所建立的节段尚局限于L₃-L₄水平，且操作过程所需计算机知识较高，需从事专门研究的工科人员协助方能完成。

实验用Mtlab科学计算软件辅助Ansys有限元软件和薄层ct技术，所建立的模型经验证符合人体标本体外生物力学实验数据，可以用来模拟正常人体L₄-L₅运动节段。并且本方法计算机操作相对简单，步骤简便，适宜医学院校学生操作，省去许多花在模型设计上的时间，可以将更多的精力放在模型的生物力学分析上。

实验腰椎融合模型的建立参考了其他文献记载的方法，已被证实是有效的。Rolander^[14]报道了在坚固的后方融合实验模型中，由于椎体弹性所造成的位移。Lee 和Langrana^[22]利用黏固剂和金属丝将人尸体的腰椎相邻的未融合节段制成融合节段。他们发现椎体间融合稳定性最好，其次是侧后方融合和后方融合。

和以前的实验报道一样，实验结果也表明椎体间融合的稳定性最好，其次是侧后方融合和后方融合。Chen等^[5]利用L₁-L₅运动节段有限元模型研究后方及侧后方腰椎融合的生物力学改变，进行了椎体间融合相邻节段的压力分析。本实验侧后方融合和后方融合的固定位置和Chen等的研究不同。Chen等的研究中侧后方融合骨条的位置在L₃-L₄椎体的横突间，不包括小关节的侧方。后方融合的骨条直接位于L₃-L₄运动节段之间薄层的后方，并用统一的材料性质模拟椎体后方的单元和骨棒，没有考虑到皮质骨和松质骨的差异。这次实验考虑到其中的不同，研究结果更为准确。Chen等的研究是用多节段椎体的模型分析相邻椎体的运动，本实验是用单节段模型进行分析，所以Chen等研究的结果和作者的实验没有可比性。但是Chen等的研究结果确实显示侧后方融合模型中骨棒所受的应力比后方模型中骨棒所受应力大。在这点上他们的实验结果和作者的结果一致。其他方面，他们的研究结果显示在屈伸运动时，侧后方融合的运动范围比后方融合的运动范围大。在这点上，作者的结论和他们的不同。Chen等的研究结果和本实验的研究结果的不同之处原因可能是两个实验所建立的融合的几何模型不同。小关节融合和侧后方融合与后方融合相结合可获得更坚强的固定，提高融合度，控制小关节疼痛。在轴向压缩载荷下，小关节融合屈曲旋转的角度最大。这表明小关节融合可限制椎体后方边缘的轴向位移，同时由于椎体骨质的弹性前方边缘可获得更大的运动幅度。

小关节融合与椎体间融合，侧后方融合和后方融合相结合，压缩载荷下轴向位移减少，减少的程度比无小关节融合分别为5%、1%和4%。同时在屈曲-伸展载荷条件下的旋转角度也相应的减小，减小的程度分别为23%、11%和45%。椎体间融合，侧后方融合和后方融合和小关节融合相结合后，活动度明显降低—尤其对于屈曲和伸展运动。作者的研究表明小关节融合和上述3种类型的腰椎融合相结合，可使腰椎获得更好的稳定性，这种稳定性的增加在后方融合中更为显著。

对于von Mises应力分布，椎体间融合无小关节融合模型的应力集中于椎体。侧后方融合合并小关节融合模型的应力集中于融合的骨棒，而后方融合合并小关节融合的应力主要位于后方的棘突。这些结果表明施加在椎体上的载荷主要由融合传导，应力主要集中于融合部位周围。后方融合合并小关节融合将施加在椎体上的载荷更有效的传导到整个椎体。

侧后方融合模型骨棒所受的应力比后方融合模型的大。由于侧后方融合模型和后方融合模型所用的骨棒大小和形状相同，这种所受应力不同的原因可能是侧后方融合模型中的骨棒比后方融合模型的骨棒更接近旋转中心。

结论：①用Matlab科学计算软件和薄层CT技术辅助Ansys有限元软件，建立正常人体L₄-L₅运动节段有限元模型的速度及精度提高。②用Ansys软件行力学分析可知：小关节融合和椎体间，侧后方，后方融合相结合，可使腰椎获得更好的稳定性，这种稳定性的增加在后方融合中更为显著。后方融合模型的应力分布更加合理。