动态冲击兔脊柱单椎体：不同冲击速率、网格数量及材料属性变化对椎体应力分布的有限元分析

doi:10.3969/j.issn.2095-4344.1417

中国组织工程研究 ›› 2019, Vol. 23 ›› Issue (30): 4828-4835.doi: 10.3969/j.issn.2095-4344.1417

• 复合支架材料 composite scaffold materials • 上一篇下一篇

动态冲击兔脊柱单椎体：不同冲击速率、网格数量及材料属性变化对椎体应力分布的有限元分析

韩世冰¹，张绪树²，郭媛²，杜一铭²

太原理工大学，¹机械与运载工程学院，²生物医学工程学院，山西省太原市 030024

收稿日期:2019-05-06 出版日期:2019-10-28 发布日期:2019-10-28
通讯作者: 张绪树，博士，副教授，太原理工大学生物医学工程学院，山西省太原市 030024
作者简介:韩世冰，男，1990年生，河北省衡水市人，汉族，太原理工大学在读硕士，主要从事生物力学方面的研究。
基金资助:
国家自然科学基金项目(11472185)，课题名称：高处坠落时脊柱损伤机制研究，项目负责人：郭媛；国家自然科学基金项目(11772214)，课题名称：人工膝关节设计参数与置换者步态和受力特征关系的生物力学研究，项目参与者：郭媛

Dynamic impact on the rabbit spine single vertebral body: finite element analysis of the stress distribution of the vertebral body with different impact rates, mesh numbers and material properties

Han Shibing¹, Zhang Xushu², Guo Yuan², Du Yiming²

¹College of Mechanical and Vehicle Engineering of Taiyuan University of Technology, Taiyuan 030024, Shanxi Province, China; ²College of Biomedical Engineering of Taiyuan University of Technology, Taiyuan 030024, Shanxi Province, China

Received:2019-05-06 Online:2019-10-28 Published:2019-10-28
Contact: Zhang Xushu, MD, Associate professor, College of Biomedical Engineering of Taiyuan University of Technology, Taiyuan 030024, Shanxi Province, China
About author:Han Shibing, Master candidate, College of Mechanical and Vehicle Engineering of Taiyuan University of Technology, Taiyuan 030024, Shanxi Province, China
Supported by:
the National Natural Science Foundation of China, No. 11472185 (to GY); the National Natural Science Foundation of China, No. 11772214 (to GY)

摘要/Abstract

摘要：

文章快速阅读：

文题释义：

单椎体：脊柱主要是通过颈椎、胸椎、腰椎、骶骨和尾骨组成，单椎体就是组成脊柱的颈椎、胸椎和腰椎。不同动物的单椎体数量不同，例如人共有7块颈椎、12块胸椎、5块腰椎，而兔子却有7块颈椎、12块胸椎、7块腰椎。

冲击：将实验样品置于冲击试验机(摆锤式和落锤式)上，对实验样品进行外力重撞，从而确定样品的力学性能。

背景：车祸和高空坠落时有发生，其中大部分造成严重的脊柱损伤，因此研究脊柱在冲击载荷作用下的损伤特性至关重要，也为脊柱的损伤防护和修复提供生物力学依据。

目的：通过有限元方法分析兔脊柱单椎体受到冲击载荷时的损伤特性及规律，分析不同冲击速率、网格数量及材料属性变化对椎体应力分布的影响。为单椎体的动态冲击实验及脊柱节段实验和模拟分析提供依据。

方法：利用Mimics、HyperMesh和Abaqus建立兔单椎体T₁₂的有限元模型，首先进行网格的敏感性分析，选择合适的网格数量，然后在不同冲击速率下对椎体进行竖直方向的冲击有限元分析和模拟，并通过调节材料属性数值可模拟骨质疏松患者的椎体在受到冲击载荷下的应力分布，最后通过兔单椎体的冲击实验对模型进行验证。

结果与结论：①通过敏感性分析发现当网格数量为74 224时，有限元分析的结果比其他网格数量下的分析结果要合理，因此实验以网格数量为74 224的前提下进行分析；椎骨材料跟大多数脆性材料一样，当受到垂直载荷作用时可发生45°角的破坏形式，并且最大应力值主要集中在椎体两端，应力传递方向与骨小梁的走向趋势一致；②结果表明，通过敏感性分析选取合适的网格数量，使有限元分析结果更加可靠；通过实验验证了该有限元模型的可靠性，有限元分析结果可在一定程度上反映椎体在不同冲击速率下和骨质疏松椎体在受到冲击下的真实情况；有限元分析方法是一种比较经济、参数易于调节和控制、适用性较强的研究方法。

ORCID: 0000-0002-0972-043X(韩世冰)
中国组织工程研究杂志出版内容重点：生物材料；骨生物材料; 口腔生物材料; 纳米材料; 缓释材料; 材料相容性；组织工程

关键词: 单椎体, 冲击, 有限元分析, 损伤, 速率, 落锤, 冲击实验, 生物力学

Abstract:

BACKGROUND: Car accidents and falling from high altitudes occur from time to time, most of which cause severe spinal injuries. Therefore, it is important to study the damage characteristics of the spine under impact load. It also provides a biomechanical basis for injury protection and repair of the spine.

OBJECTIVE: The finite element analysis method was used to analyze the damage characteristics and regularity of the single vertebral body of rabbit under impact load, and the effects of different impact velocities, mesh numbers and material properties on the stress distribution of the vertebral body were analyzed. These provided the basis for the dynamic impact test of single vertebral body and the segmental experiment and simulation analysis of the spine.

METHODS: The finite element model of rabbit single vertebral T₁₂was established by Mimics, HyperMesh and Abaqus. First, the sensitivity analysis of the mesh was carried out to select the appropriate number of meshes. Finite element analysis and simulation of vertical impact on vertebral bodies were conducted at different impact velocities. Stress distribution of vertebral body in osteoporotic patients under impact loading was simulated by adjusting the numerical value of material properties. The model was validated by a single vertebral impact test in rabbits.

RESULTS AND CONCLUSION: (1) Through sensitivity analysis, it was found that when the number of grids was 74 224, the results of finite element analysis were much more reasonable than those of other grids. Therefore, the analysis was based on the premise of the number of grids 74 224. Vertebral material like most brittle materials, a 45° angle of failure could occur when subjected to vertical loads. The maximum stress value was mainly concentrated at the two ends of the vertebral body. The direction of stress transmission was consistent with the trend of the trabecular bone. (2) It is concluded that the appropriate number of grids was selected by sensitivity analysis, which makes the results of finite element analysis more reliable. The reliability of the finite element model was verified by experimental verification. The results of finite element analysis can reflect the true situation of the vertebral body under different impact velocities and osteoporotic vertebral body under impact. The finite element analysis method is economical and has strong control and applicability; and the parameters are easy to adjust.

Key words: single vertebra, impact, finite element analysis, damage, speed, drop hammer, impact test, biomechanics

中图分类号:

韩世冰，张绪树，郭媛，杜一铭. 动态冲击兔脊柱单椎体：不同冲击速率、网格数量及材料属性变化对椎体应力分布的有限元分析[J]. 中国组织工程研究, 2019, 23(30): 4828-4835.

Han Shibing, Zhang Xushu, Guo Yuan, Du Yiming. Dynamic impact on the rabbit spine single vertebral body: finite element analysis of the stress distribution of the vertebral body with different impact rates, mesh numbers and material properties [J]. Chinese Journal of Tissue Engineering Research, 2019, 23(30): 4828-4835.

图/表（结果） 6

2.1 不同网格数量对椎体的影响—敏感性分析

2.1.1 网格数量对特征点O应力的影响以冲击速率为4.4 m/s为例，椎体体网格数量25 243，32 910，47 949，65 297和74 224在图1中特征点O的Von Mises应力最大值分别为14.28，16.68，15.62，13.4，19.53 MPa，最大应力值的变化较大，可见网格数量对特征点处O处的影响较大，通过图3不同数量网格划分下特征点O的应力-时间曲线所示，网格数量多时应力随时间的变化曲线比较稳定。

2.1.2 网格数量对等效塑性应变PEEQ的影响以分析时间为0.02 s时为例，效塑性应变PEEQ云图如图4所示，椎体体网格数量74 224，65 297，47 974，32 910，25 243下的等效塑性应变最大值分别为0.604 9，0.789 8，0.611 0，0.144 6，0.821 5，数值范围较大，但前3种网格的等效塑性应变变化比较小，且5种网格下的等效塑性变最大值位置一致，都集中在椎体上端。

2.1.3 网格无关性分析^[18] 表3为不同网格数量时椎体的真应变、主应力和整体位移的比较。图5为不同网格数量下的应变值和应力值。通过表3可发现随着网格数量的增加，波动比逐渐增大，当网格数量到达65 297时波动比开始减小。从图5不同网格数量下的应变值和应力值可看出，当网格数量到达47 949，65 297，74 224时，曲线明显变缓。因此，当网格数量为74 224或更高时有限元模型的网格数量较为合理。

2.2 不同冲击速率对椎体应力和位移的影响

2.2.1 Von Mises应力、主应力和轴向应力通过不同冲击速率下对椎体T₁₂的冲击模拟，得到了椎体的von Mises应力、主应力及轴向应力。以0.03 s时的受力情况为例，如表4所示，von Mises应力的最大值均约等于20 MPa，而椎体材料参数的屈服应力均设置为20 MPa，这说明在7种冲击速率下椎体均达到屈服极限，而且最大von Mises应力的节点位置集中在椎体两端，0.03 s时的位移云图如图6所示。随着冲击速率的增大，最大主应力和最大轴向应力逐渐增大，并且可发现同一时刻的主应力与轴向应力的最大值相差不大。

2.2.2 von Mises应力分布以10.0 m/s的von Mises应力云图为例，获得T₁₂在冲击载荷下的应力分布，横截面如图7所示。截面内von Mises应力极值出现在上下两端即A、B处，C处为椎体内部应力最大值，椎体内部其他区域的应力值较小。因此，椎体T₁₂在受到冲击载荷时其应力传递方式是上下两端受到冲击载荷，应力达到极值，然后应力值沿AB轴线向椎体内部逐渐增大，在椎体内部C处达到最大值，最后以C点为中心，再沿DE方向，水平向椎体外侧逐渐增大。

2.2.3 位移云图分布特点以4.4，6.0，7.5，9.0 m/s为例，分析了椎体在0.03 s时椎体的整体位移。通过图8可发现4种冲击速率下整体位移的最大值分别为2.934，4.109，5.972，7.750 mm，且最大值的位置都位于椎体上端。通过位移云图可发现椎体有明显的分界线并与轴线呈45°角，而分界线的位置就是椎体破坏的位置，在前期静态压缩实验中也有类似的结果。

2.2.4 特征点O在竖直方向上的Stress components应力-时间曲线分析特征点O在竖直方向上Stress components应力随时间变化的曲线，如图9所示。通过曲线可以发现特征点O处的单元在椎体受到冲击的瞬间应力达到最大值，在速率较小时最大值与冲击速率成正比，除 4.4 m/s外，其他冲击速率下的极值均在22.5 MPa左右，然后应力随后开始减小，并趋于一个稳定的值。

2.3 不同材料属性对椎体应力的影响

2.3.1 不同材料属性对轴向应力与位移的影响不同材料属性对轴向应力与位移的影响见图10所示，从图中可发现轴向应力与密度、弹性模量呈正相关关系，随着密度、弹性模量的增加其所承受的轴向应力也逐渐增加。而材料属性变化对轴向位移的影响则是随着材料密的减小，轴向位移则是先减小，后逐渐增大，约80%时为最小值。而且当密度、弹性模量较小时的轴向位移更大，当材料密度为原来的60%，70%时的轴向位移远大于100%时的轴向位移。

2.3.2 不同材料属性对特征点O处的影响以冲击速率7.5 m/s为例，不同材料属性下特征点O的应力-时间曲线如图11所示。通过分析曲线可发现不同的材料属性对特征点处的应力影响不明显。

2.3.3 不同材料属性对整体位移云图的影响以t=0.01 s时的受力为例，不同材料属性下的整体位移云图如图12所示，通过位移云图可以发现，不同材料属性下的位移云图变化较为明显，从100%到60%的5种材料属性下的位移最大值分别为4.648，4.927，4.766，4.625，4.859 mm，可见材料属性的不同对位移云图的影响较大。

2.4 实验验证根据有限元不同冲击速率对椎体应力和位移的影响得出：在7 J的冲击能量下椎体已发生完全破坏，因此在Instron CEAST 9350落地式落锤试验机上对单椎体T₁₂进行落锤冲击实验，冲击能量为7 J，实验冲头重400.3 g，冲击设备除冲头外还包括传感器，总质量 5.925 3 kg，根据牛顿第二定律可得出冲击速度为1.5 m/s。冲击效果如图13所示。有限元模拟与实验的载荷-时间曲线的对比图如图14所示，通过对比图可发现从0-0.005 s内实验的载荷曲线和有限元的载荷曲线趋势相同，主要经历了2个阶段，下关节突和椎弓的破坏为第1个阶段，第2个阶段为椎体的破坏。2个曲线图虽然趋势相同但在细节方面存在一点差异，例如载荷的大小和时间节点的对应。造成差异的主要原因是：首先，有限元模拟和分析尽可能还原实验状态但仍存在一些差异，一个是冲击载荷，有限元冲击的冲头为5 kg，冲击速度为1 m/s，相对于实验环境冲击载荷较小；其次，有限元分析设置材料属性与实验时椎体有一定的差异，由于目前国内外现有的研究资料中缺乏新西兰大白兔椎体的材料属性经验公式，而实验有限元模型的材料属性参考人股骨的经验公式，导致材料的杨氏模量变大，因此有限元分析的结果略大于实验的结果；最后，实验中冲击能量过大，导致椎体变形过大，变形程度达到90%以上，因此载荷作用时间较长，而有限元分析为简化模型，作用时间较短，椎体变形程度不足50%，因此载荷作用时间较短。在实际生活中，脊柱的椎骨受到肌肉、韧带、椎间盘等软组织的保护，椎骨变形不会达到实验中的情况，因此分析前2个峰值的变形即可。

参考文献

[1]贾涛,赵晶.1053例脊柱损伤的临床流行病学研究[J].创伤外科杂志,2003,5(4):256-258.

[2]曾至立,程黎明.胸腰椎爆裂性骨折损伤机制的研究进展[J].中华创伤骨科杂志,2011,13(8):773-775.

[3]Giuseppe S,Alessandra B,Hugo G,et al.Biomechanical effects of metastasis in the osteoporotic lumbar spine: A Finite Element Analysis.BMC Musculoskelet Disord.2018;19(1):38.

[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-664.

[5]John JD, Saravana Kumar G, Yoganandan N.Cervical spine morphology and ligament property variations: A finite element study of their influence on sagittal bending characteristics.J Biomech. 2019;85:18-26.

[6]Wang Md K, Jiang PhD C, Wang PhD L,et al.The biomechanical influence of anterior vertebral body osteophytes on the lumbar spine: A finite element study.Spine J. 2018;18(12):2288-2296.

[7]蔡康健,王丽珍,姚杰,等.腰椎椎体有限元建模的最优单元尺寸和材料属性分布及建模方法[J].医用生物力学, 2016,31(2): 135-141.

[8]于研,曾至立,朱睿,等.垂直冲击下L1应力分布的有限元研究[J].同济大学学报(医学版),2012,33(3): 5-9.

[9]Su JC, Li ZD, Cao LH,et al.Three-dimensional finite element analysis of lumbar vertebra loaded by static stress and its biomechanical significance.Chin J Traumatol. 2009;12(3): 153-156.

[10]Park WM, Park YS, Kim K,et al.Biomechanical comparison of instrumentation techniques in treatment of thoracolumbar burst fractures: a finite element analysis.J Orthop Sci. 2009; 14(4):443-449.

[11]秦大平,张晓刚,聂文忠,等.不同运动状态下模拟人体腰椎结构特征变化的有限元分析[J].医用生物力学, 2017,32(4):355-362.

[12]张国栋,廖维靖,陶圣祥,等.股骨有限元分析赋材料属性的方法[J].中国组织工程研究与临床康复,2009,13(43):8436-8441.

[13]张晶,郭媛,张绪树,等.侧向撞击时肩部有限元分析[J].太原理工大学学报,2018,49(6):886-892.

[14]张晶.肩部三维有限元模型的建立及侧向撞击中的力学分析[D].太原:太原理工大学,2018.

[15]李洲,韩世冰,张玉文,等.兔单椎体压缩性骨折的应力-应变关系[J].中国组织工程研究,2019,23(4): 538-544.

[16]李洲.兔脊柱单椎体压缩性骨折的实验研究[D].太原:太原理工大学,2018.

[17]张健,赵桂平,卢天健.泡沫金属在冲击载荷下的能量吸收特性[J].西安交通大学学报,2013,47(11): 105-112.

[18]王宝中,于玉真.电动车车架数值分析下的网格特性分析与研究[J].机械设计与制造,2018,56(3):98-101.

[19]荀福兴.不同材料属性分配方法对椎体有限元模型力学性能的影响[D].广州:南方医科大学,2013.

[20]Kopperdahl DL, Morgan EF, Keaveny TM.Quantitative computed tomography estimates of the mechanical properties of human vertebral trabecular bone.J Orthop Res. 2002;20(4): 801-805.

[21]崔壮,余斌,许长鹏,等.骨密度对距下关节融合初稳性影响的三维有限元分析[J].中华创伤骨科杂志,2013, 15(3):220-224.

[22]Alonso Vazquez A,Lauge-Pedersen Hlidgren L,Taylor M.The effect of bone quality on the stability of ankle arthrodesis. A finite element study.Foot Ankle Int.2004;25(11):840-850.

[23]沈宇辉,袁高翔,郁健,等,骨质疏松性股骨粗隆间骨折类型与骨密度关系及临床意义[J].国际骨科学杂志, 2011,32(6):389-391.

引言

进入21世纪以来，随着中国经济的飞速发展，人们的出行速度越来越快，城市建筑越来越高，但由此给人们带来的人体损伤也逐年增加。据不完全统计，交通伤和坠落伤已成为人们受伤的主要形式，并且在这些伤者中大部分都会有不同程度的脊柱损伤^[1]。因此，研究交通事故、坠落等导致的脊柱损伤机制至关重要。随着有限元的发展，通过建立椎体有限模型研究脊柱损伤机制的方法越来越成熟^[2]。国内外许多学者在脊柱损伤方面进行了深入研究。Giuseppe等^[3]通过建立L₃-L₅有限元模型分析腰椎骨转移的数量与骨密度对椎体稳定性的影响。Qiu等^[4]通过对T₁₂-L₁节段有限元模型进行动态垂直冲击模拟，研究了胸腰椎交界处椎体爆裂性骨折的发生机制。John等^[5]基于网格变形的参数化和基于模型的灵敏度分析研究脊柱形态变化对屈伸效应的影响。Wang等^[6]通过有限元分析探讨前路椎体骨质增生对整个腰椎的生物力学影响。宫赫教授的团队通过最优单元尺寸和材料属性的分布方法在建模方面进行了详细介绍^[7]。于研等^[8]通过有限元研究椎体L₁在冲击载荷下的应力分布。Su等^[9-10]通过有限元分析计算了不同应力加载部位的腰椎机械承受能力。秦大平等^[11]建立正常与病变的腰椎有限元模型，分析了两者在不同运动状态下的应力特点和量效关系。由于脊柱损伤多发生在胸腰椎结合处，故而相关的研究也集中在此。

因此，实验采用有限元分析方法建立胸椎与腰椎连接处的椎体，也是在日常生活中最易发生损伤的椎体——T₁₂的三维有限元分析模型，研究单个T₁₂椎体的动态冲击特性，分析不同冲击速率、网格数量及材料属性变化对椎体应力分布的影响。

中国组织工程研究杂志出版内容重点：生物材料；骨生物材料; 口腔生物材料; 纳米材料; 缓释材料; 材料相容性；组织工程

材料方法

1.1 设计生物力学有限元分析和实验。

1.2 时间及地点兔单椎体的冲击有限元分析于2018年6至10月在太原理工大学生物医学工程学院完成；兔单椎体冲击实验于2019年4月在太原理工大学工程训练中心完成。

1.3 材料冲击实验仪器为Instron CEAST 9350落地式落锤试验机；三维重构软件Mimics 19.0(Materialise，比利时)；逆向工程软件Geomagic Design X(Geomagic，美国)；HyperMesh 14.0(Altair，美国)；Abaqus软件(Simulia，美国)。

实验动物：兔龄191 d新西兰大白兔整个脊柱，由太原理工大学生物医学工程学院提供，在山西医科大学第二医院西院进行CT扫描。

1.4 实验方法

1.4.1 几何模型的建立建模数据来自1只成年新西兰大白兔脊柱的CT断层扫描图像，经CT扫描并未发现脊柱畸形等其他病变，对C₁-L₇共26个椎体进行了扫描，共获得1 136张DICOM格式的连续断层图像(分辨率512×512，像素尺寸0.601 563 mm，层厚0.5 mm)。以T₁₂为研究对象，将CT扫描数据导入三维重构软件Mimics 19.0中，通过对骨质部分进行阈值设置并分割图像，建立了整个脊椎骨的蒙板；最后对蒙板的裁剪和修复，建立胸椎T₁₂的几何模型。

将生成的几何模型以STL的格式导入逆向工程软件Geomagic Design X中生成T₁₂的几何面片，然后将T₁₂几何面片以STL的格式导入HyperMesh 14.0中划分体网格，生成四面体单元，并在椎体上下两端各建立一个刚性圆盘，使用Abaqus分析和模拟动态落锤冲击过程，并与前期的兔单椎体压缩实验相对应。

1.4.2 有限元模型的建立 HyperMesh软件是有限元分析常用的前处理软件，网格划分功能十分强大。通过HyperMesh对几何模型进行网格划分，具体划分的网格类型为四面体，单元数量为74 224，节点数量为15 357，最小单元边长0.2 mm，雅各比系数0.6。图1为建立的整个T₁₂三维有限元分析模型，其中O点为选取椎体中央位置的特征点，单元的编号为11 728。

1.4.3 材料属性的确定将划分的网格重新导入Mimics 中赋予材料属性，椎体设置成弹塑性材料，在Mimics中赋予材料时，将椎体看成非均匀材料，材料密度的分类越多越能接近真实的材料属性，但这显然会增加有限元分析的计算量和计算时间。张国栋等^[12]在股骨有限元分析赋材料属性的方法中，通过10，20，50，100，200和400种材料属性探讨合理的材料种类数量，并指出骨骼材料赋予10种材料属性就基本能满足分析需求^[13-14]。因此实验基于材料灰度值的不同，利用经验公式将整个椎体分成10种不同的材料属性^[12]。

Density=-13.4+1 017×Grayvalue(g/m³)

E-Modulus=-388.3+5 925×Density(Pa)

公式中的Density为密度，Grayvalue为灰度值，E-Modulus为弹性模量。

在前期兔单椎体静态压缩实验获得的兔单椎体T₁₂应力-应变曲线基础上，通过分析曲线获得T₁₂的屈服应力和塑性应变分别为20 MPa和0.05^[15-16]。由于缺乏可直接参考的兔椎体密度、材料强度等材料参数的经验公式，故兔椎体材料属性采用了椎体整体压缩的平均参数。椎体的有限元模型的材料如表1所示。

为了分析有限元模型网格数量对椎体应力的影响，在保持椎体面网格数量为11 176不变的情况下，将椎体体网格数量分别划分为25 243，32 910，47 949，65 297和74 224五种，其对应的节点数量分别为7 356，8 586，11 114，14 600和15 357，然后将模型以4.4 m/s的速率进行冲击模拟，分析不同网格数量对椎体应力的影响。

在不同冲击速率下研究了椎体的最大von Mises应力、主应力和轴向应力，椎体横截面内的应力分布，以及椎体的位移分布和特征点O处应力的竖直分量随时间的变化曲线。

骨质疏松患者椎体的承载能力明显下降，此部分通过将椎体的材料属性，例如密度、杨氏模量等分别调整为原来的90%，80%，70%和60%后，来模拟不同程度骨质疏松疾病的椎体力学特性，椎体有限元模型的材料参数如表2所示。

1.4.4 接触、约束和载荷设置将椎体与上、下两刚性圆盘之间设置为通用接触，并使用Abaqus中的面—面接触算法。将椎体与刚性圆盘之间的切向方向摩擦系数象征性取为0.3，以模拟椎体与刚性圆盘之间有摩擦，目前的数据库中并未有椎骨与刚性圆盘之间的摩擦数据^[17]。法向方向设置成硬接触，考虑接触非线性。

将下面的刚性圆盘固定，并对其6个运动自由度全部约束。把上刚性圆盘作为模拟落锤竖直下落冲击实验的落锤，同时为了简化计算和增加冲击接触的稳定性，将Z轴平动自由度以外的其他5个自由度全部约束。对整个模型施加重力载荷(刚性圆盘的质量为0.44 kg，直径为30 mm)，并对上刚性圆盘施加初速度，以此来模拟落锤下落后接触椎体前的速度，将椎体置于圆盘的正中央。整个动态冲击模拟过程在Abaqus的动态显示求解器中完成，同时考虑几个非线性，求解完成后在Abaqus的可视化页面对结果进行后处理和数据分析。

1.4.5 有限元和实验静态压缩的对比将有限元分析模型的静态压缩结果与前期兔单椎体静态压缩实验结果进行对比^[15-16]，并且有限元分析采用类似于文献[15-16]中静态压缩实验的固定和加载方式。将上面刚性圆盘施加位移载荷-最大压缩位移为8 mm，模拟Instron5544材料性能实验机的实验工况进行分析。考虑到静态压缩与动态冲击的不同，将上面刚性圆盘与椎体的接触方式设为方程约束，下面刚性圆盘与椎体的接触方式设为多点约束，对比结果如图2所示。通过有限元和实验静态压缩对比图可发现，实验和有限元分析的压缩载荷随时间变化曲线的趋势比较吻合，这可在一定程度上验证有限元模型的可靠性，但是曲线数值上存在一些差异，造成差异的主要原因是有限元分析设置材料属性与实验时椎体有一定的差异，由于目前国内外现有的研究资料中缺乏新西兰大白兔椎体的材料属性经验公式，而此次实验有限元模型的材料属性参考人股骨的经验公式，导致材料的杨氏模量变大，因此有限元分析结果略大于实验结果。

1.5 主要观察指标实验和有限元分析的压缩载荷随时间变化曲线的趋势是否吻合。

中国组织工程研究杂志出版内容重点：生物材料；骨生物材料; 口腔生物材料; 纳米材料; 缓释材料; 材料相容性；组织工程

讨论

实验通过三维重构软件、逆向工程软件及有限元分析软件的有效结合和使用，建立了兔脊柱椎体T₁₂的三维有限元分析模型，在不同冲击速率、单元数量和材料属性的情况下对椎体进行了动态竖直冲击的有限元分析和模拟。多种手段的、多种方法的结合，对生物力学的研究和发展起到了巨大的推动作用。三维重构软件Mimics可以对CT、MRI等医学断层扫描数据进行分割和处理，可有效地建立组织和器官比较真实、可靠的几何模型。对于生物材料，例如脊柱椎骨大都不是一种均匀的材料，但在分析和计算中，有些学者为简化模型，在研究时将其当做一种均匀材料，也有些学者为了更接近实际材料而将其作为一种非均匀材料^[19]。宫赫教授的团队研究发现，将骨赋予均匀的材料属性和分配成非均匀的材料属性，得到的轴向刚度存在明显差异^[7]。Abaqus作为一种有限元分析软件，并不能方便、快速地赋予单元非均匀的材料属性，而Mimics可跟据材料灰度值的不同将材料赋予非均匀的材料属性^[7，20]。

实验采用了非均匀的材料属性，模拟体外环境下椎体在不同冲击速率、网格数量和材料属性情况下受到的冲击载荷作用，分析椎体应力等参量。在7种冲击速率下，Mises应力均达到屈服应力极限——20 MPa。考虑到使用的材料本构关系为弹塑性材料模型，没有考虑单元的失效，故而椎体未有发生断裂破坏现象。通过有限元分析结果也可发现，应力极值出现的位置主要集中在椎体上下两端，而于研等^[8]的研究也表明在垂直载荷作用下，椎弓根基底部是应力集中的主要区域。其他学者的实验研究也证实了这一观点，他们通过测量尸骨表面应力证实了应力最高的区域发生在椎体两端。通过分析椎体整体位移-时间曲线中的界限也可得出椎骨材料跟大多数脆性材料一样，在受到垂直载荷作用后也会发生45°角破坏的形式。

Qiu等^[4]通过对椎体微观的结构研究发现，椎体两端的骨密度、骨小梁等含量要明显高于其他部位，依据生物力学中的相关理论，骨小梁含量多的部位应力值越大，而且骨小梁的走向趋势与应力传递方向一致。通过截面的应力云图也可发现椎体在受到冲击载荷后的应力传递方式与此结论相一致。

在有限元分析中，网格划分的数量、密度、质量及单元类型等会对有限元分析结果产生重要影响。因此，实验分析并比较了网格数量对有限元分析结果的影响，发现不同数量的网格对椎体应力等参量存在一定的影响，尤其对位移的影响较为明显，而单元数量多的模型，其应力等参量随时间变化的曲线比较稳定。王宝中等^[18]在网格特性分析与研究中也指出网格数量越多其结果越精确，但会增加计算成本。因此，需要找到一个临界值来解决两者之间的矛盾，这就需要做网格无关性分析。实验通过分析数据和比较发现当网格数量为74 224时，这个网格数量下的有限元分析相对于结果与其他网格数量下的结果要合适的多。

南方医科大学的崔状等^[21]根据AlonsoVazquez等^[22]提供的方法，逐渐降低正常骨的弹性模量来建立骨质疏松模型。此次实验也通过类似方法降低骨密度和弹性模量，分别达到正常水平的90%，80%，70%和60%来模拟不同骨质疏松情况下椎体的力学性能，分析不同材料属性对有限元分析结果的影响。有限元模拟结果发现随着骨密度的降低，材料力学性能也在逐渐降低，其所能承受的应力也越来越低，而在同等冲击能量下发生的位移压缩程度则不断增大。在断裂力学中，比较疏松的材料在同等工况下其破坏程度更大。沈宇辉等^[23]也发现较低骨密度人群发生骨折的概率更大，这与此次实验结果非常相似。

实验研究了不同冲击速率、网格数量和材料属性情况下兔脊柱T₁₂椎体的应力等力学参量分布特点和规律，仍有不足与需改进之处：①模型为考虑合理的失效准则，故破坏现象未能显现；②基于灰度赋予材料属性时，未建立兔脊柱的材料密度、强度与灰度关系的经验公式，故而对有限元的分析结果会有一定的影响；③文中进行了简化研究，也未开展椎体小关节对冲击性能影响的相关研究；④对模型验证方面，仅进行了静态的压缩实验和有限元对比和分析，在下一步开展动态冲击实验后进行详细的动态模型验证。

实验得出如下结论：①通过敏感性分析选取合适的网格数量，使有限元分析的结果更加可靠；②通过实验验证了该文有限元模型的可靠性，因此该文有限元分析结果可在一定程度上反映椎体在不同冲击速率下和骨质疏松椎体在受到冲击下的真实情况；③有限元分析方法是一种比较经济、参数易于调节和控制、适用性较强的研究方法。

中国组织工程研究杂志出版内容重点：生物材料；骨生物材料; 口腔生物材料; 纳米材料; 缓释材料; 材料相容性；组织工程

动态冲击兔脊柱单椎体：不同冲击速率、网格数量及材料属性变化对椎体应力分布的有限元分析

Dynamic impact on the rabbit spine single vertebral body: finite element analysis of the stress distribution of the vertebral body with different impact rates, mesh numbers and material properties

PDF

可视化

摘要/Abstract

引用本文

使用本文

图/表（结果） 6

参考文献

相关文章 15

引言

材料方法

讨论

文章快阅

延伸阅读

编辑推荐

Metrics

本文评价

[1]	闵友江, 姚海华, 孙洁, 周璇, 余航, 孙前谱, 洪恩四. 磁共振评价“三通针法”对脊髓损伤患者脑功能的影响[J]. 中国组织工程研究, 2021, 25(在线): 1-8.
[2]	徐峰, 康辉, 魏坦军, 席金涛. 椎弓根螺钉不同固定方法治疗胸腰椎骨折的生物力学分析[J]. 中国组织工程研究, 2021, 25(9): 1313-1317.
[3]	陈心敏, 李文标, 熊凯凯, 熊晓燕, 郑利钦, 李木生, 郑永泽, 林梓凌. 钉道强化股骨近端防旋髓内钉治疗老年A3.3型股骨转子间骨折：最佳骨水泥量有限元分析[J]. 中国组织工程研究, 2021, 25(9): 1404-1409.
[4]	周继辉, 李新志, 周游, 黄卫, 陈文瑶. 髌骨骨折修复内植物选择的多重问题[J]. 中国组织工程研究, 2021, 25(9): 1440-1445.
[5]	黄登承, 王志科, 曹学伟. 体外冲击波疗法治疗中老年膝骨关节炎短期疗效对比的荟萃分析[J]. 中国组织工程研究, 2021, 25(9): 1471-1476.
[6]	蒋红英, 朱亮, 余曦, 黄靖, 向小娜, 兰正燕, 何红晨. 富血小板血浆干预脊髓损伤患者压力性损伤的作用[J]. 中国组织工程研究, 2021, 25(8): 1149-1153.
[7]	曾祯, 胡经纬, 李璇, 唐琳梅, 黄志强, 李明星. 超声造影定量分析重度失血性休克再灌注模型大鼠复苏期的肾血流灌注[J]. 中国组织工程研究, 2021, 25(8): 1201-1206.
[8]	柴乐, 吕建兰, 胡劲涛, 胡华辉, 许庆军, 余进伟, 全仁夫. 诱导急性脊髓损伤模型大鼠炎症反应信号通路的变化[J]. 中国组织工程研究, 2021, 25(8): 1218-1223.
[9]	陈继铭, 吴晓静, 刘田丰, 陈海聪, 黄成硕. 水飞蓟素对四氯化碳致小鼠肝损伤和骨代谢的影响[J]. 中国组织工程研究, 2021, 25(8): 1224-1228.
[10]	万然, 史旭, 刘京松, 王岩松. 间充质干细胞分泌组治疗脊髓损伤的研究进展[J]. 中国组织工程研究, 2021, 25(7): 1088-1095.
[11]	谢文佳, 夏天娇, 周卿云, 刘羽佳, 顾小萍. 小胶质细胞介导神经元损伤在神经退行性疾病中的作用[J]. 中国组织工程研究, 2021, 25(7): 1109-1115.
[12]	曾燕华, 郝延磊. 许旺细胞体外培养及纯化的系统性综述[J]. 中国组织工程研究, 2021, 25(7): 1135-1141.
[13]	孔德胜, 何晶晶, 冯宝峰, 郭瑞云, Asiamah Ernest Amponsah, 吕飞, 张舒涵, 张晓琳, 马隽, 崔慧先. 间充质干细胞修复大动物模型脊髓损伤疗效评价的Meta分析[J]. 中国组织工程研究, 2021, 25(7): 1142-1148.
[14]	宋成杰, 常恒瑞, 石明鑫, 孟宪中. 侧方入路腰椎融合治疗后的生物力学稳定性的研究与进展[J]. 中国组织工程研究, 2021, 25(6): 923-928.
[15]	刘钊, 徐西林, 申意伟, 张晓峰, 吕航, 赵军, 王政春, 刘旭卓, 王海涛. 塌陷预测方法联合分期分型对股骨头坏死治疗的指导作用与前景[J]. 中国组织工程研究, 2021, 25(6): 929-934.