此文章已经发表于中华医学杂志2011年第31期

有限元分析在脊柱生物力学研究中越来越广泛，但是大部分有限元法的研究均是在静态加载下的静力分析，而对于临床上常见的胸腰椎爆裂性骨折损伤而言，主要是受到垂直冲击暴力而致。研究胸腰椎爆裂性损伤过程中的应力分布及变化，以及进一步研究损伤机制，需要在有限元模型上模拟体外的损伤过程，进行冲击试验并进行动态分析。据我们所知，目前还没有报道一个完整的可进行动态冲击实验并进行动态分析的胸腰椎三维非线性有限元模型，本研究拟建立一个胸腰椎（T11-L3）的三维非线性有限元模型，为下一步在有限元模型上进行动态冲击的生物力学实验奠定基础。上海市同济医院脊柱外科曾至立

1、材料和方法

1.1 胸腰椎有限元模型的实体建模

1.1.1样本数据采集

建模数据采集自接近国人50 百分位体型的健康成年男性志愿者(年龄25 岁,身高172cm ,体重64kg，无脊柱疾病史和明显脊柱外伤史,经摄片排除脊柱畸形及其他病变)一成年男性志愿者的CT断层扫描图像。采用GE公司Light Speed 16排螺旋CT扫描机，人体位于扫描视野的中心，扫描条件为：120kV，250 mA，层厚0.625mm，螺距1.25mm，自上而下进行螺旋扫描，扫描范围包括T11-L3椎体共得到519张CT图片，通用DICOM3.0格式保存。

1.1.2 三角网格模型的建立

通用DICOM3.0格式读入医学有限元建模软件Simpleware2.0中（Simpleware Ltd., UK，它包含3个模块建模模块ScanIP、组合装配模块ScanCAD与有限元前处理模块ScanFE）。首先将原始图像重新采样，适当调整原始图像的分辨率从而调整最终生成的实体模型的精细度。再基于CT灰度值将图像通过阈值工具分割出感兴趣的区域，并且精细区分出不同骨性结构之间的界限。最后运用对感兴趣区域进行填充，得到胸腰椎的三角网格模型。对于椎间盘，由于其CT灰度与周围软组织十分接近，无法通过阈值分割直接获得，所以按照解剖知识通过对椎体之间的空间的填充获得椎间盘的轮廓。对原始图像的降噪滤波、边界处理、阈值分割与特征提取，最终建立起胸腰椎的三角网格模型。该模型包含T11、T12、L1、L2、L3共5个椎体以及之间的4个椎间盘（图1）。把生成的三角网格模型作为Geomagic8(Raindrop Geomagic, Inc., USA)的导入文件。

1.1.3实体模型的建立

将由Simpleware2.0中生成三角网格模型导入逆向工程软件Geomagic8中进行修补与优化。首先在保证模型形态精度的前提下减少三角面片的数量，为进一步优化模型提供基础。然后填充模型表面空洞，去除模型表面非特征性的尖刺和压痕，平滑松弛表面，防止非特征性的高曲率及自相交面的产生，实现用光顺平滑的NURBS曲面来拟合离散的模型表面三角面片，最后生成具有G1连续的NURBS曲面模型即实体模型（图2）。

          

图1 胸腰椎的三角网格模型                图 2 腰胸腰椎的NURBS实体模型

1.2有限元分析前处理

以通用IGES格式将胸腰椎的实体模型导入有限元前处理软件  HyperMesh10.0中。通过对5节椎骨以及4个椎间盘网格的划分、材料属性定义和边界条件的设定，完成有限元分析前处理。

1.2.1网格划分

骨性结构区分皮质骨（厚度为1mm）和松质骨、后部结构以及骨性终板（厚度为0.5mm），这些结构均采用一阶六面体单元C3D8进行网格划分。椎间盘（图3）位于骨性终板之间，根据组织学发现按比例由髓核（44%）和纤维环(56%)组成。髓核由5层8节点实体单元组成，纤维环有8层相互交叉胶原纤维层组成，纤维和椎间盘平面呈35°角[1]，相邻的二层纤维走向相互交叉，呈110°。髓核与纤维环采用一阶杂交的减缩积分六面体单元C3D8RH进行网格划分，胶原纤维采用2节点桁架TRUSS单元。按照文献描述的韧带的起止位置[2]，并赋予所有韧带1mm的厚度，通过4节点减缩积分壳单元S4R模拟了前纵韧带、后纵韧带、黄韧带、棘上韧带、棘间韧带、横突间韧带，采用3节点S3R壳单元模拟小关节囊韧带。单元数及类型汇总见表1

 

图3椎间盘由髓核和纤维环组成，纤维环由8层胶原纤维和基质组成

图4 胸腰椎有限元模型的网格，包括椎体、椎间盘和七条重要韧带。

 

表1 有限元模型各结构单元类型及单元数

结构

松质骨

皮质骨

后部结构

纤维环

髓核

胶原纤维

终板

关节囊韧带

其他韧带

单元类型

C3D8

C3D8

C3D8

C3D8RH

C3D8RH

T3D2

C3D8

S3R

S4R

单元数

24288

7268

9068

4830

2760

11036

4048

6529

9060

 

1.2.2 材料属性设置

根据文献资料[3]，骨性结构假设为各向同性材料，按照Johnson-Cook定律赋予对称性弹塑性特性；韧带按照Maxwell-Kelvin-Voigt粘弹性材料定律赋予粘弹性属性，韧带纤维组织使用非线性材料属性；椎间盘按照Mooney-Rivlin材料定律赋予超弹性属性。椎体、椎间盘及韧带的材料特性总结如表（表2、3、4）。

表2  骨性结构材料属性

材料属性

皮质骨

松质骨

终板

骨密度（kg/mm3）

1.83E-06

0.17E-06

1.06E-06

杨氏模量，E(MPa)

14000

291

10000

泊松比，v

0.3

0.25

0.3

屈服应力  a（MPa）

110

1.92

6

硬化模量  b（MPa）

100

20

100

硬化指数 n

0.1

1

1

失效塑变应变 eP

9.68E-03

14.5E-03

0.02

最大应力（MPa）

155

2.23

7.5

应变率系数 c

1

1

3

表3 椎间盘材料属性

材料属性

髓核

纤维环基质

胶原纤维

密度（Kg/mm3）

1.00E-06

1.20E-06

非线性曲线

泊松比，v

0.495

0.45

 

C10

0.12

0.18

 

C01

0.03

0.045

 

 

 

 

 

表4 韧带材料属性

材料属性

ALL

PLL

ITL

ISL

LF

SSL

JC

密度（Kg/mm3）

1.0E-06

1.0E-06

1.0E-06

1.0E-06

1.0E-06

1.0E-06

1.0E-06

杨氏模量，E(MPa)

11.4

9.12

11.4

4.56

5.7

8.55

22.8

泊松比，v

0.4

0.4

0.4

0.4

0.4

0.4

0.4

切线模量，Et（MPa）

10.0

9.0

11.0

4.0

5.0

8.0

22.0

切线泊松比，vt

0.42

0.42

0.42

0.42

0.42

0.42

0.42

粘度系数，ho

28

28

28

28

28

28

28

Navier氏常数，l

1.0E-06

1.0E-06

1.0E-06

1.0E-06

1.0E-06

1.0E-06

1.0E-06

ALL:Anterior longitudinal ligaments前纵韧带、 PLL: posterior longitudinal ligaments后纵韧带、 FL :flavum ligaments黄韧带、 SSL: supraspinous ligaments棘上韧带、 ISL: interspinous ligaments棘间韧带、  ITL: intertransverse ligaments 横突间韧带 、CL: capsular ligament关节囊韧带

 

1.2.3相互作用方式设定

在椎间盘与骨性结构之间使用共节点连接方式。韧带于骨性结构的界面使用Tied焊接方式，阻止在有限元模拟过程中出现组织之间的相互移动。小关节突关节面之间采取无摩擦通用接触。

模型分析使用ABAQUS 6.9的EXPLICIT的显式动力学有限元求解器进行。

1.2.4载荷与边界条件设置

约束L3下终板6个方向自由度作为边界条件，采用Distribution Coupling的方式约束T11上终板作为从节点，建立其与主节点（旋转轴上某一点）的分布式载荷，将载荷加载于主节点上。

2、模型的验证

2.1 椎间盘的载荷-位移实验

参照有关文献[3]在有限元模型上模拟体外实验过程，在准静态应力下分析椎间盘的载荷-位移曲线（图5），将模型中的每个运动单元（二个椎体和其之间的椎间盘，但不包括后部结构）逐步进行准静态面加载，加载方式1.267mm/s，加载时间1s，得出每个椎间盘的载荷-位移曲线(图6)和文献中的结果[3]比较。

图5 有限元模型上椎间盘载荷-位移实验

图6 椎间盘的载荷-位移实验验证（*数据来源于参考文献[3]）

2.3 冲击载荷下大体应力分布

将T11上终板单元采用offset偏置的方法，在已建立的有限元模型T11椎上方覆盖厚3mm的刚体板，覆盖范围包括椎体及后部结构。依据文献[4]建立半径11mm的密度均匀的球形刚体，刚板的密度设置为7.8E-6kg/mm3 ，钢球的密度为9.97E-4Kg/mm3（质量为5.56Kg）。模拟体外的冲击锤，钢球以6.0m/s的速度冲击钢板中央，按照能量计算公式E=(1/2)×m×v2冲击能量为100J垂直轴向冲击有限元模型，动态观察应力分布云图。

 

图7 100J轴向冲击下应力云图（左为整体右侧面；右为L1皮质骨）

3 结果

建立了胸腰椎的三维有限元模型，包括T11、T12、L1、L2和L3五个脊椎及四个椎间盘和前纵韧带、后纵韧带、黄韧带、棘上韧带、棘间韧带、囊韧带及横突间韧带七条重要韧带。共包含71939个节点和78887个单元。模型各个部分分别设定了材料参数，赋予骨性结构弹塑性、椎间盘超弹性、韧带粘弹性材料属性。

在和体外实验在相同条件下进行了有限元模型的准静态椎间盘的载荷-位移实验，每个椎间盘共获得20个数据，得到的椎间盘载荷-位移曲线图（图 6）和文献报道体外实验结果[3]非常相近；在动态冲击实验下观察应力云图，发现应力集中位于椎弓根周围、小关节突关节及棘突和上关节突之间的椎板中间，应力最大值出现在椎弓根外下缘和椎体交界的皮质骨区，动态观察发现应力最早集中于该部位，且应力沿皮质骨向椎体逐步扩散；松质骨的应力明显较皮质骨小，在中央偏后方应力较大，动态观察发现应力最早集中出现中央偏后部位并向四周扩散（图7），这和Dai[5]的有限元研究结果一致。

4 讨论

4.1三维有限元模型的建模方法

有限元分析法是计算力学中的一种重要方法, 它是20世纪40年代初兴起的应用数学、现代力学及计算机科学等学科相互渗透、综合利用的一门边缘科学。1972年Brekelmans[6]首次将有限元法用于骨的应力分析。近20年来, 有限元分析法在脊柱生物力学研究中的应用日益广泛和深入。由于脊柱解剖结构复杂性,应用实体建模法较直接生成法更合适 [7]。

为了精确的描述胸腰椎的几何形态，对建模方法进行改良，将模型的建立分为三角网格模型与NURBS曲面模型两个主要步骤来实现。三角网格模型是用大量的三角网格对原物件表面进行离散的近似的三维模拟，对复杂几何拓扑结构有较强描述力，可以很大程度上还原标本的解剖形状。但是由于三角网格模型被处理为离散的几何体，表面的曲率不连续[8]，直接导入有限元软件中无法生成体，导致无法对模型作进一步的分析计算。采用NURBS曲面拟合离散的三角网格,则能构成具有G1连续的、封闭的曲面，由封闭的曲面构成了实体，进而进行有限元分析。这样结合了三角网格对复杂标本解剖形状的描述能力和NURBS曲面的连续性，可以保证了后期有限元分析的精度和可行性。

将NURBS曲面保存为通用IGS格式后，将其导入到有限元前处理软件HyperMesh10.0中。HyperMesh10.0是一个高效的有限元前后处理器，能够建立各种复杂模型的有限元和有限差分模型，与多种cad和cae软件有良好的接口并具有高效的网格划分功能。在HyperMesh10.0中，通过高质量的六面体网格对实体模型进行划分，有利于计算机资源的节省，并为后面冲击载荷下骨折模型的建立提供模型基础。

4.2三维有限元模型的验证

有限元模型的验证问题一直是一个很困难的问题[9]。验证三维有限元模型的方法一般有[10]：1、取实际标本模型进行同样条件下的实验，观察同样的指标，对比二者的结果进行验证；2、与现存文献中类似的实验结果进行比较；3、将实验得到的最终结果与临床观察到的现象进行验证。

由于实体模型的相似性程度直接关系到后续有限元分析的可靠性， Robin发现几何形态对脊柱的力学行为有重要的影响[11]，Ha等[12]指出，在进行有限元分析时，需要尽可能精确的模拟标本的几何形状，以避免简化所带来的精度的降低。本实验模型的几何外形来源于健康成年人的无脊柱疾病史的三维超薄CT扫描数据，经过专业的医学有限元建模软件Simpleware2.0、逆向工程软件Geomagic8及有限元划分软件HyperMesh10.0，建立了具有六面体网格的、包含71939个节点和78887个单元的有限元模型，充分保证了几何相似度；模型进行了准静态的椎间盘的载荷实验，观察其载荷和位移的关系，结果和文献[3]的体外实验高度相似；模型动态冲击实验下发现应力集中位于椎弓根周围、小关节突关节及棘突和上关节突之间的椎板中间，应力最大值出现在椎弓根外下缘和椎体交界的皮质骨区，和其他胸腰椎三维有限元模型[5]获得的应力分布规律一致，并且和脊椎的骨小梁分布规律一致[13、14]，即在骨小梁密度密集及骨密度高的部位应力高。本研究通过以上三个方面保证了有限元模型的几何相似性和力学相似性和胸腰椎高度一致，为进行下一步的进行冲击实验并分析应力、得到真实的结果提供了可靠基础。

4.3 本模型的特点及意义

近些年来已经有不少学者建立了胸腰椎的有限元模型以研究临床相关问题。 Dai[5]通过一个运动节段的共293个单元及388节点简化的线性三维有限元模型对胸腰椎爆裂性骨折的机制进行了分析。Liu等[15]通过两个运动节段的线弹性有限元模型分析了在创伤载荷下退变椎间盘对胸腰段椎体应力分布的影响。Wilcox[16]等通过三个节段的有限元模型研究的胸腰段爆裂性骨折机制及椎管占位，但是在此模型中，损伤椎体的上下相邻椎体仅仅有部分被模拟，并且在受伤机制中起重要作用的韧带仅包含了后纵韧带。为了研究冲击载荷下胸腰段爆裂性骨折的机制， Qiu等[4]通过有限元方法一个运动节段的有限元模型上分析了垂直冲击载荷下应力分布和终板位移，其模型设定均为线弹性，对结果有一定影响。为了弥补既往模型的不足（相关节段较短，忽略了部分解剖结构，多为四面体网格，网格质量一般，材料线弹性等），本研究建立了包括T11-L3的全胸腰椎三维有限元模型，模型精确的区分了皮质骨、松质骨，终板，后部单元，椎间盘（区分了髓核和纤维环，纤维环由其中的纤维加强），并且模拟了力学行为重要的七条韧带，采用六面体网格，优化了网格质量，赋予模型非线性材料属性。本模型为T11-L3节段，在进一步进行胸腰椎损伤机制的研究中可以不分析受冲击脊椎（T11）和约束脊椎（L3）应力变化，而研究不受到载荷和边界条件影响的T12、L1、L2脊椎的应力分布及变化，更加具有科学性。本研究使用的HyperMesh10.0进行了高质量的网格划分，除韧带及胶原纤维外均采用更高质量的六面体网格。模型采用了采用一阶的四边形壳网格与六面体体网格，这是因为在相同阶数下，它们相对于三角形壳网格与四面体体网格有更高的精度与更小的计算代价；采用杂交单元模拟间盘结构材料的超弹属性；采用减缩积分单元是因为该单元类型在大变形工况下能有效减少单元“沙漏”现象的产生，防止由于单元的剪切自锁而导致的计算结果不收敛。此方法网格质量高，可以准确的反映应力、应变、位移等力学参数，可以大幅度降低计算资源的消耗，为冲击载荷的瞬态加载和骨折机制的模拟提供必要的基础。在材料属性方面，最大程度的模拟正常人体组织结构的非线性属性以及压缩载荷时显示的非线性行为，除了胸腰椎的骨性结构赋予弹塑性材料属性外，并赋予椎间盘超弹性的材料属性，使椎间盘接近生理状态。韧带方面，通过精确的解剖位置和非线性粘弹性的材料属性，最大限度的恢复其生理作用。本模型经过准静态的椎间盘载荷-位移实验及动态冲击分析椎体的应力分布验证了模型的力学相似性。因此，本模型是一个精细、有效、完整胸腰段的三维非线性有限元模型，可进行动态有限元分析。为进一步在有限元模型上模拟胸腰椎损伤过程，进行冲击加载试验，分析椎体、椎间盘及韧带内部的应力分布及变化规律，研究胸腰椎损伤机制奠定基础。

 

参考文献

1.   Schmidt H, Kettler A, Heuer F, et al. Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading. Spine (Phila Pa 1976). 2007;32(7):748-55.

2.   Pintar FA, Yoganandan N, Myers T, et al. Biomechanical properties of human lumbar spine ligaments. J Biomech. 1992;25(11):1351-6.

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

4.   Qiu TX, Tan KW, Lee VS, et al. Investigation of thoracolumbar T12-L1 burst fracture mechanism using finite element method. Med Eng Phys. 2006;28(7):656-64.

5.   Dai L. Mechanism of thoracolumbar burst fractures: a biomechanical study. Chin Med J (Engl). 2002;115(3):336-8.

6.   Brekelmans WA, Poort HW, Slooff TJ. A new method to analyse the mechanical behaviour of skeletal parts. Acta Orthop Scand. 1972;43(5):301-17.

7.   朱睿, 程黎明. 骶椎切除重建的三维有限元模型建立及其研究进展. 医用生物力学. 2008;23:327-31.

8.   Sun W, Starly B, Nam J, et al. Bio-CAD modeling and its applications in computer-aided tissue engineering. Computer-Aided Design. 2005;37(11):1097-114.

9.   Pope MH, Novotny JE. Spinal biomechanics. J Biomech Eng. 1993;115(4B):569-74.

10.  李全, 蔡郑东, 都承斐, 等. 骶骨部分切除术后骨盆重建的三维有限元分析. 医用生物力学. 2008;23:47-51.

11.  Robin S, Skalli W, Lavaste F. Influence of geometrical factors on the behavior of lumbar spine segments: a finite element analysis. Eur Spine J. 1994;3(2):84-90.

12.  Ha SK. Finite element modeling of multi-level cervical spinal segments (C3-C6) and

biomechanical analysis of an elastomer-type prosthetic disc. Medical Engineering & Physics. 2006

;28(6):534-41.

13.  Heggeness MH, Doherty BJ. The trabecular anatomy of thoracolumbar vertebrae: implications for burst fractures. J Anat. 1997;191(2):309-12.

14.  Dai LY, Wang XY, Wang CG, et al. Bone mineral density of the thoracolumbar spine in relation to burst fractures: a quantitative computed tomography study. Eur Spine J. 2006;15(12):1817-22.

15.  Liu L, Pei F, Song Y, et al. The influence of the intervertebral disc on stress distribution of the thoracolumbar vertebrae under destructive load. Chin J Traumatol. 2002;5(5):279-83.

16.  Wilcox RK, Allen DJ, Hall RM, et al. A dynamic investigation of the burst fracture process using a combined experimental and finite element approach. Eur Spine J. 2004;13(6):481-8.