有限元分析法在腰椎生物力学中的研究进展

冯其金，赵玲娟，郑昆仑，谷福顺

1972年Belytschko等[1]首次将有限元分析应用于脊柱力学领域，历经几十年发展，有限元模型已由二维线性发展到三维非线性，模型中不仅包括了骨性结构，也加入了肌肉组织特性。同时，有限元分析由静态响应向动态响应过渡。有限元方法用来分析生理和病理过程中的力学变化，为设计和改进手术器械提供理论依据，可以术前模拟手术过程，预测手术效果，优化手术方案。本文以有限元方法用于腰椎各解剖部位生物力学为例，对目前的有限元研究近况进行统一分析。

1 有限元法的概念及原理

有限元方法的基本原理是把一个由无限个质点构成并且有无限个自由度的连续体划分为有限个小单元体组成的集合体。此过程称为有限元模型的离散化，并用这些离散化的单元模型代替原有的物体，根据原有物体的几何材料特性及受力条件采用不同的单元种类，单元之间通过节点相连接，节点之间相互传递的力称为节点力。每个单元的物质特性及节点载荷、边界条件明确后，通过节点、位移与节点力之间的关系式计算出每个单元的刚度矩阵。若干个单元的刚度矩阵集合成构件的总刚度矩阵，并通过数学形式表达出来。适用于复杂的几何形状和边界条件，而且能处理各种复杂的材料。它通过有限单元、节点、自由度等实验条件控制，用计算机模拟人体体内情况。

2 有限元方法在腰椎各解剖结构的应用

1975年Liu等[2]首先报道了腰椎的三维有限元模型；1979年Hakim等[3]则在一节腰椎模型加入了后部附件，分析了整个运动节段静力学和动力学的问题之后，一些研究相继用有限元方法分析了腰椎活动节段在不同类型载荷作用时的应力分布。

2.1 有限元法在腰椎椎体方面的研究 Eswaran等[4]认为皮质骨最大的承载比例是54%，出现在椎体横切面最小的部位，松质骨则是89%，出现在靠近终板的部位。皮质骨承载功能还与其距离终板的位置相关，靠近终板约承担34%，而在上下终板之间约为63%。颜文涛等[5]根据CT图像建立L4～5有限元模型，并对模型施予不同载荷，使模型产生不同的轴向位移，通过分析，符合腰椎生物力学特性的应力分布规律，并与体外实验结果相吻合，从而证明其模型真实有效。陈浩等[6]通过建立L4～5运动节段的流固耦合有限元模型，监测24 h椎间盘的高度变化，将所建模型椎间盘高度变化与体外实验相比较，指出腰椎流固耦合有限元模型较体外实验更加符合人体的生物力学特点，证明其有效。

2.2 有限元法在腰椎间盘方面的研究 闫家智等[7]研究表明，在施加轴向压缩力时，腰椎纤维环后部纤维的应力较高，最大应力集中于髓核和终板中央，应力随轴向压缩力的增加而增大。El-Rich等[8]研究表明，俯屈时应力集中于L2椎弓根区域的上面，其次是L2下方的终板。伸展时最大应力分布于L2椎弓根区域的下面。椎间盘应力在俯屈时增大，而关节接触应力在伸展时增加。Guo等[9]研究认为如果在去除载荷时终板的渗透性增加，则在轴向和径向应变增大，丢失的液体将以更快的速度恢复。Fan等[10]通过研究在不同工作模式下的不同节段的退变椎间盘生物力学效应时指出随着休息频率的增加，腰椎间盘的轴向应力及液体流失下降，各生物力学参数的变化率亦不同，并认为在总的静息时间不变的情况下，高频的休息是可取的。如果健康人得不到足够的休息则他们会出现类似于轻中度腰椎间盘退变患者由于缺少休息的影响。然而，在椎间盘严重退化的情况下，一次休息是没什么意义的，所以合理的工作休息模式，可能会帮助不同程度的椎间盘退变患者实现更好的恢复。

2.3 有限元法在腰椎关节突关节方面的研究 Du等[11]研究发现在不同的载荷条件下，小关节的生物力学关系不同。载荷作用下增大了小关节接触面积和接触压力。在腰椎侧弯时同侧小关节变化与上述相同，但对侧则相反。但是在扭转载荷下小关节几乎没有影响。然而随着载荷的增加对下关节的影响越小,并指出在腰椎的生物力学研究中，特别是在负载下的生物力学研究中，应充分考虑不对称性对关节突关节的影响。Holzapfel[12]研究了腰椎小关节曲率在应力分析中的作用，认为在分析矢状面的剪切力/位移和接触应力分布时，忽视小关节的曲率会发生较大的偏差。Claeson等[13]通过研究腰椎小关节囊韧带在脊柱运动中的形变指出关节突关节囊的前后表面中部平面内和通过平面的不均衡的形变很明显。由于关节间隙较窄，故前表面的形变更明显。而后表面的形变较分散，因为其较大的关节间隙增加了形变能力,并认为这中形变在脊柱运动中起着重要作用。Woldtvedt等[14]运用有限元分析法分析关节突软骨厚度及分布特点与腰椎活动的关系，得出关节突软骨分布厚度及分布对腰椎各个方向的活动度影响较小，而对关节突间的应力改变影响较大的结论。

2.4 有限元法在腰椎旁肌肉和韧带的研究 Zander等[15]将肌肉作用力加入腰椎的三维非线性有限元模型中，设计了不同的上身屈曲角度，评估在上身前倾过程中，肌肉作用力对腰椎间盘纤维环的应力分布的影响。此研究发现屈曲角度固定，施加载荷不同，椎间盘所受应力大小和分布差异显著；当脊柱承受矢状面屈曲力矩时，纤维环腹侧承受压应力，纤维环背侧承受张应力；当脊柱同时承受体质量和背侧肌肉作用力时，应力叠加，导致椎间盘腹侧压应力增加，椎间盘背侧张应力减小。分析显示，肌肉作用力加入前后，椎间盘的应力分布差异显著。Kaminska等[16]认为腰椎等稳定性与肌肉骨骼的紊乱程度相关，运用有限元分析法对其相关性进行分析，指出腰椎间盘对应力大小与腰椎负荷成正比，腰椎前屈时应力最大。Zander等[17]建立了L3/L4有限元模型，模拟依次切断韧带，计算剩余韧带的应力。分析发现：屈曲时，后纵韧带的受力最大，后伸和侧曲时，前纵韧带受力最大，旋转时，关节囊韧带受力最大；屈曲时，棘间韧带对运动节段内的旋转影响最大，后伸和屈曲时，前纵韧带对旋转影响最大，关节囊韧带对轴向旋转影响最大，韧带的刚度明显影响节间的活动范围和韧带的应力。

3 有限元法在腰椎病理状态及手术方式中的应用

Sairyo等[18]经有限元分析认为，单侧的椎体滑脱使对侧峡部应力增加，导致对侧峡部应力性骨折和硬化，提示在临床上单侧椎体滑脱的患者，特别是那些活动水平不高的患者，如果有持续的下腰痛主诉，要考虑对侧是否发生了应力性骨折。潘志强等[19]通过相关研究发现不管是单节段L4-5，还是双节段L4-S1固定后，邻近L3-4椎间盘有效应力在前屈、后伸、垂直压缩、侧屈、旋转时均大于未融合固定，差异显著；各加载方式两两比较，加载平均应力由大到小依次为旋转、侧屈、后伸、前屈和垂直压缩。并指出L4-5或L5-S1节段的旋转力矩会明显提高L3-4椎间盘的应力，随后依次是侧屈、后伸、前屈和垂直压缩的应力，充分证实旋转是引起椎间盘慢性损伤及退变的最大风险因素，所以双节段腰椎融合手术较单节段腰椎融合手术更容易发生邻近节段病。Wang等[20]经过研究指出在退行性变引起的脊柱侧弯情况下当Cobb角较小时，关节突关节的凸侧受力较大。当Cobb角大于20°，关节突关节凹侧受力较大。在轴向旋转的情况下，同一水平的关节突关节面受压不同。并认为用非对称加载负荷时，关节面压缩变形出现在凹侧，并降低对关节的影响来限制椎体的旋转和移位。小关节接触力的不对称性加载加速腰椎的不对称性。费琦等[21]研究了对于术后引起的腰椎不稳的单节段和多节段的后路腰椎融合术的有限元分析。指出由于行椎板及小关节部分切除引起的腰椎不稳可以通过腰椎间融合来治疗。但是腰椎融合会增加邻近节段退变的可能性，因为相邻椎间盘的应力将增大。

Ambati等[22]通过研究经椎间孔椎体间融合术的不同术式及内固定物数量和方式的有限元分析指出，在所有重建方式下L4/L5的活动度均减小，在左侧屈曲时，双侧后路椎弓根螺钉固定的稳定性高于单侧固定。单侧固定时将保留原活动度的50%，而双侧固定时则仅有10%的原活动度。单侧固定时后植入物应力在屈曲和侧曲时分别是正常的6倍和4倍。并认为经椎间孔双侧后路内固定增强了融合结构的稳定性，并减少了后部内固定物的应力。在双侧椎弓根螺钉固定时椎间内植物的形状和数量不会影响节段的稳定性。Newcomb等[23]通过研究螺钉位置对腰椎椎弓根螺钉内载荷传递的影响的有限元分析指出椎弓根螺钉的位置在影响螺钉松动及体内断钉的风险方面起着重要作用，包括矢状面和轴面的角度。放置一颗在轴面的横向螺钉可以通过减少松质骨应力及螺钉应力来降低断钉与螺钉松动的风险。Tang[24]通过比较后路椎间融合与经椎间孔椎间融合对腰椎退变性变化的有限元分析指出PLIF和TLIF两种模型下的椎间盘内压力与邻近节段间的旋转度均高于正常。但TLIF模型在各个方向载荷下时其椎间盘内压力与椎间旋转度均低于PLIF组。并认为椎间孔椎间融合术优于后路椎间融合术。

4 有限元法在腰椎生物力学研究的问题与展望

有限元法是研究人体腰椎生物力学的有力工具，是一种行之有效的方法。经过几十年的发展，有限元法在脊柱力学分析方面的研究更趋深入。但是有限元分析也受到各种因素的限制，比如模型的构建，不同模型的几何形状及材质的定义不同等。其与体外实验相比，有限元法有其优缺点。体外实验对一些周围软组织的处理有一定困难，且很多应力估算只能通过骨表面的应变来实现，而无法对骨内部的应力进行估算。有限元法则能弥补这一点，所以有限元法应与体外实验相互补充与验证才能发挥其最大价值，并继续推动脊柱生物力学更深入的发展。

参考文献：

[1] Belytschko TB,Andriacchi TP,Schultz AB,et al.Analog studies of forces in the human spine: computational techniques[J].J Biomech,1973,6(4):361-371.

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

[3] Hakim S,King KL.Finite element Methods in spine research[J].J Biomech,1979,12(5):277.

[4] Eswaran SK,Gupta A,Adams MF,et al.Cortical and trabecular load sharing in the human vertebral body[J].J Bone Miner Res,2006,21(2):307-314.

[5] 颜文涛,赵改平,方新果,等.人体腰L4-5节段有限元建模及分析[J].生物医学工程学杂志,2014,31(3):612-618.

[6] 陈浩,于晓华,华国军,等.腰椎运动节段流固耦合有限元模型的建立与验证[J].中国组织工程研究与临床康复,2010,14(52):9706-9709.

[7] 闫家智,吴志宏,汪学松,等.腰椎间盘退变后应力变化的有限元分析[J].中国医学科学院学报,2009,31(4):464-467.

[8] El-Rich M,Arnoux PJ,Wagnac E,et al.Finite element investigation of the loading rate effect on the spinal load-sharing chanes under impact conditions[J].J Biomech,2009,42(9):1252-1262.

[9] Guo LX,Li R,Zhang M.Biomechanical and fluid flowing characteristics of intervertebral disc of lumbar spine predicted by poroelastic finite element method[J].Acta Bioeng Biomech,2016,18(2):19-29.

[10] Fan R,Gong H,Qiu S,et al.Effects of resting modes on human lumbar spines with different levels of degenerated intervertebr al discs: a finite element investigation[J].BMC Musculoskelet Disord,2015,16(1):221.

[11] Du CF,Yang N,Guo JC,et al.Biomechanical response of lumbar facet joints under follower preload:a fi nite element study[J].BMC Musculoskelet Disord,2016,17(1):126.

[12] Holzapfel GA,Stadler M.Role of facet curvature for accurate vertebral facet load analysis[J].Eur Spine J，2006,15(6):849-856.

[13] Claeson AA,Barocas VH.Computer simulation of lumbar fl exion shows shear of the facet capsular ligament[J].Spine J，2017,17(1):109-119.

[14] Woldtvedt DJ,Womack W,Gadomski BC.et al.Finite element lumbar spine facet contact parameter predictions are affected by the cartilage thickness distribution and initial joint cap size[J].J Biomech Eng,2011,133(6):325-331.

[15] Zander T,Rohlmann A,Calisse J,et al.Estimation of muscle forces i n the lumbar spine during upper-body inclination[J].Clin Biomech(Bristol,Avon),2001,16(1):S73-80.

[16] Kam iń ska J, Roman-Liu D, Zagrajek T,et al.Differences in lumbar spine load due to posture and upper limb external load[J].Int J Occup Saf Erqon,2010,16(4):421-430.

[17] Zander T,Rohlmann A,Bergmann G.Analysis of simulated single ligament transection on the mechanical behaviour of a lumbar functional spinal unit[J].Biomed Tech (Berl),2004,49(1-2):27-32.

[18] Sairyo K,Katoh S,Sasa T,et al.Athletes with unilateral spondylolys is are at risk of stress fracture at the contralateral pedicle and pars interarticularis:a clinical and biomechanical study[J].Am J Sports Med,2005,33(4):583-590.

[19] 潘志强,叶立民,王求永,等.腰椎融合术后临近节段椎间盘应力的有限元分析[J].中外医学研究,2014,12(17):20-22.

[20] WangL,ZhangB,ChenS,et al.A Validated Finite Element Analysis of Facet Joint Stress in Degenerative Lumbar Scoliosis[J].World Neurosurg,2016,95(11):126-133.

[21] 费琦,赵凡,杨雍,等.腰椎后路融合手术对失稳模型节段稳定性及相邻节段力学的影响[J].中华医学杂志,2015,95(45):3681-3686.

[22] Ambati DV,Wright EK Jr,Lehman RA Jr,et al.Bilateral pedicle screw fixation provides superior biomechanical stability in transforaminal lumbar interbody fusion:a fi nite element study[J].Spine J,2015,15(8):1812-1822.

[23] Newcomb AG,Baek S,Kelly BP,et al.Effect of screw position on load transfer in lumbar pedicle screws:a non-idealized finite element analysis[J]. ComputMethodsBiomech Biomed Engin,2017,20(2):182-192.

[24] Tang S.Comparison of posterior versus transforaminal lumbar interbody fusion using finite element analysis.Influence on adjacent segmental degeneration[J].Saudi Med J,2015,36(8):993-996.