不同材料属性赋予下椎体的静力学分析

doi:10.12307/2022.644

摘要/Abstract

摘要：

文题释义：
材料属性赋予：由于骨骼的非线性及骨质分布不均匀，学者们在进行椎体模拟分析时进行了不同方法的材料属性赋予来提高分析效率。
图像分割：根据CT数据灰度阈值进行区域生长、形态学操作、布尔操作、动态区域生长、多层编辑等进而提取出所需要的部位。

背景：人体骨骼具有非线性强、骨质分布不均匀等特点，不同的弹性模量赋值导致静力学结果差异性大。
目的：分析3种材料属性赋予方法在椎体数值模拟分析中的差异性。
方法：基于逆向工程原理，通过图像分割建立腰椎L4、L5及椎间盘模型，对椎体模型赋予3种材料属性：①灰度值赋值：通过Mimics利用灰度值进行梯度赋值；②均质及灰度值混合赋值：对椎体前端皮质骨进行单一弹性模量赋予，对前端松质骨及后端椎弓根、横突、棘突等利用灰度值进行赋予；③皮质骨及松质骨赋值：采用皮质骨及松质骨两种弹性模量。对3种赋值方法下的有限元模型进行轴向载荷100-500 N加载，以及分别施加2，4，6，8，10 N·m的弯矩模拟椎体进行前屈、后伸、左右侧弯、左右旋转。
结果与结论：不同工况加载下的椎体轴向位移及偏转角度有明显差异，均方根误差分析显示，灰度值赋值法相比于其他两种方法更接近于真实实验数据，不会在材料赋值交界处产生明显的应力阶跃现象，以及不会产生明显的应力集中现象，更符合椎体的力学性质，从而可以更准确地进行脊柱生物力学分析。

https://orcid.org/0000-0003-1663-537X (关天民)

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

关键词: 材料属性, 生物力学, 灰度值, 弹性模量, 腰椎, CT数据, 逆向工程

Abstract: BACKGROUND: Human bones are characterized by strong nonlinearity and uneven bone distribution. Different elastic modulus assignments lead to great differences in static results.
OBJECTIVE: To analyze differences of the three material attribute assigning methods in the numerical simulation analysis of vertebral body.
METHODS: Based on the principle of reverse engineering, lumbar L4, L5 and intervertebral disc models were established through image segmentation, and three material attributes were assigned to the vertebral body models: (1) Gray value assignment: Gradient value was assigned by using gray value through Mimics. (2) Homogeneous and gray value mixed assignment: A single elastic modulus was given to the anterior cortical bone of the vertebral body, and gray values were used to give to the cancellous bone of the front end and the pedicle, transverse process and spinous process of the posterior end. (3) Cortical bone and cancellous bone assignment: Two elastic moduli of cortical bone and cancellous bone were used. The finite element model under the three assignment methods was loaded with an axial load of 100-500 N. The bending moments of 2, 4, 6, 8, and 10 N·m were applied to simulate the vertebral body to perform forward flexion, backward extension, left and right lateral bending and left and right rotation.
RESULTS AND CONCLUSION: Vertebral body under different loading conditions had obvious differences in axial displacement and deflection angle. Root mean square error analysis showed that compared with the other two methods, data obtained by grey value assignment method were more close to the real experimental data. It did not produce obvious stress step phenomenon in material assignment junction, and could not produce obvious stress concentration phenomenon, was more in line with the mechanical properties of vertebral body. Thus, more accurate statics analysis of the spine can be carried out.

Key words: material properties, biomechanics, grey value, elastic modulus, lumbar vertebra, CT data, reverse engineering

中图分类号:

关天民, 陈向禹, 朱晔, 任东. 不同材料属性赋予下椎体的静力学分析[J]. 中国组织工程研究, 2022, 26(21): 3357-3362.

Guan Tianmin, Chen Xiangyu, Zhu Ye, Ren Dong. Different material attributes are assigned to the vertebral body for statics analysis[J]. Chinese Journal of Tissue Engineering Research, 2022, 26(21): 3357-3362.

图/表（结果） 11

2.1 轴向加载验证及椎体应力应变情况如图6所示，同BROWN等[24]实验进行对比，随着轴向载荷的增加，椎体L4上端面位移随之增加，几乎呈线性关系，说明3种材料属性赋值方法的有限元模型在正常受力状态下都表现出弹性性质。

椎体表面及截面Mises 应力分布云图如图7所示，从表面应力云图可以看出，3种材料属性赋值方法随着载荷的增加，应力首先产生于皮质骨处，A计划组较B、C计划组集中应力范围较小，从截面图可以看出，B、C计划组在皮质骨-松质骨交界面产生了明显的应力阶跃现象，B计划组皮质骨-松质骨应力阶跃现象较C计划组小。

2.2 不同工况下各组椎体应力应变情况
2.2.1 前屈、后伸工况应力应变情况前屈弯矩加载偏转角度变化如图8a所示，椎体Mises应力分布云图如图9a所示。从偏转角度曲线中可以看出，随着弯矩的增加，3种方法偏转角度呈线性增加，A计划组的偏转角度较 C、B计划组明显大一些，HEUER等[25]、NACHEMSON等[26]为相同工况下人体实验数据。Mises应力云图显示随着载荷的增加，集中应力逐渐产生在椎体前端面，A计划组集中应力保持梯度变化，B、C计划组应力集中较明显。
后伸弯矩加载偏转角度如图8b所示，Mises应力分布云图如图9b所示。3种方法下集中应力均产生于椎弓根周围，C计划组集中应力最明显，B、A计划组其次，符合椎体力学特性；从椎弓根截面可以看出，C计划组在皮质骨-松质骨材料赋予交界面产生了明显的应力阶跃现象，A、B计划组由于椎弓根周围采用梯度赋予材料属性，在皮质骨-松质骨交界面没有出现明显的应力阶跃现象。

2.2.2 侧弯工况下应力应变情况在侧弯工况加载下，左侧弯偏转角度曲线如图10a所示，右侧弯偏转角度曲线如图10b所示，左、右侧弯曲线均表现出线性关系。左侧弯 Mises 应力云图如图11所示，集中应力分布在偏转的椎体一侧，这符合椎体的力学特性；3种材料属性相比，依旧表现出在皮-松质骨交界面处B、C计划组出现了明显的应力阶跃现象。

2.2.3 旋转工况应力分布在旋转工况加载下，左旋转偏转角度曲线如图12a所示，右旋转偏转角度曲线如图12b所示，左右旋转曲线均表现出线性关系，两侧旋转角度曲线相似。左旋转Mises应力云图如图13所示，载荷加载下集中应力产生于椎弓根周围，分布特点符合相应工况下特征，依旧表现出在皮质骨到松质骨的材料赋予交界面处 C计划组产生明显的应力集中现象。

2.3 不同工况下各组椎体偏转角度及均方根误差分析将3种方法所得的偏转角度曲线与同类加载方式HEUER等[25]、NACHEMSON等[26]的人体实验数据对比，均方根误差结果如表3所示。

对相同工况下同一文献数据中3种方法的均方根误差进行排名，例如在前屈中与HEUER等实验数据进行对比，均方根误差为A计划组2.336 6<C计划组2.637 8<B计划组2.763 5，误差越小说明与实验数据越接近，依次得分为A计划组得3分，C计划组得2分，B计划组得1分。对其他工况下依次进行排名，总分为A计划组为27分，B计划组为20分，C计划组为25分。
2.4 不同工况下椎体最大应力情况 6种载荷工况下椎体最大应力值如图14所示。

3种方法最大应力值对比发现，除左、右侧弯加载，最大应力绝大部分均出现在B计划组中；在左、右侧弯工况中，最大应力出现在C计划组，但是B计划组最大应力值与C计划组非常接近。可以得出，3种方法中B计划组更容易产生应力集中，C、A计划组次之。

参考文献

[1] 孙枫原,李宗远,何希,等.有限元分析在脊柱侧凸生物力学研究中的应用及进展[J].中国组织工程研究,2019,23(32):5221-5226.
[2] LIANG L, LIU ML, MARTIN C, et al. A deep learning approach to estimate stress distribution: a fast and accurate surrogate of finite-element analysis. J R Soc Interface. 2018;15(138):20170844.
[3] 赵鹏飞,陈玲,门玉涛.有限元分析肌肉力对腰椎内固定系统的影响[J].中国组织工程研究,2018,22(23):3654-3658.
[4] 刘慧,沈国权,张喜林,等.肌肉加载下腰椎间盘突出的有限元研究[J].医用生物力学,2019,34(5):493-499.
[5] CHOI HW, KIM YE. Effect of lumbar fasciae on the stability of the lower lumbar spine. Comput Method Biomec. 2017;20(13):1431-1437.
[6] ANDERSSON G, SCHULTZ AB. Effects of fluid injection on mechanical properties of intervertebral discs. J Biomech. 1979;12(6):453-458.
[7] TENCER AF. Some static mechanical properties of the lumbar intervertebral joint, intact and injured. J Biomech Eng. 1982;104(3):193-201.
[8] 文毅,苏峰,刘肃,等.L4-5椎体有限元模型建立及退变椎间盘力学分析[J].中国组织工程研究,2019,23(8):1222-1227.
[9] 薛志鹏,李泰贤,李龚,等.基于CT灰度值赋值的股骨头坏死有限元模型对比[J].中国组织工程研究,2020,24(3):395-400.
[10] 罗林聪,马立敏,林泽,等.基于AnyBody骨骼肌肉多体动力学分析的有限元仿真[J].医用生物力学,2019,34(3):237-242+250.
[11] 张聪,赵岩,杜小宇,等.青少年特发性脊柱侧凸腰主弯患者腰椎-骨盆的生物力学分析[J].中国组织工程研究,2020,24(8):1155-1161.
[12] 成博,赵文,罗伟,等.基于CT三维重建模型在数字化骨科的应用研究[J].武汉理工大学学报,2017,39(2):89-94.
[13] RASMUSSEN J, TRHOLM S, ZEE MD. Computational analysis of the influence of seat pan inclination and friction on muscle activity and spinal joint forces. Nt J Ind Ergonom. 2009;39(1):52-57.
[14] 李丹.基于腰椎多层螺旋CT扫描三维形态学分析的腰椎材料、形态及结构属性变化与骨折相关性的FEA研究[D].长春:吉林大学,2011.
[15] MORGAN EF, BAYRAKTAR HH, KEAVENY TM. Trabecular bone modulus-density relationships depend on anatomic site. J Biomech. 2003;36(7):897-904.
[16] RHO JY, KUHN-SPEARING L, ZIOUPOS P. Mechanical properties and the hierarchical structure of bone. Med Eng Phys. 1998;20(2):92-102.
[17] 聂文忠.脊柱胸腰部的生物力学建模与应用研究[D].上海:上海交通大学, 2009.
[18] 苏晋.腰椎有限元模型的建立与生物力学分析[D].大连:大连医科大学,2010.
[19] 贾少薇,张顺心,范顺成,等.脊柱侧凸腰骶椎结构的有限元分析及其变形趋势[J].医用生物力学,2017,32(3):235-241.
[20] 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;18(11):1531-1541.
[21] ROHLMANN A, BOUSTANI HN, BERGMANN G, et al. Effect of a pedicle-screw-based motion preservation system on lumbar spine biomechanics: A probabilistic finite element study with subsequent sensitivity analysis. J Biomech. 2010;43(15): 2963-2969.
[22] 国家技术监督局.成年人人体质心.中华人民共和国国家标准GB/T 17245-1998.1998.
[23] YAMAMOTO I, PANJABI MM, CRISCO T, et al. Three-dimensional movements of the whole lumbar spine and lumbosacral joint. Spine (Phila Pa 1976). 1989; 14(11):1256-1260.
[24] BROWN T, HANSEN RJ, YORRA AJ. Some mechanical tests on the lumbosacral spine with particular reference to the intervertebral discs;a preliminary report. J Bone Joint Surg Am. 1957;39(5):1135-1164.
[25] HEUER F, SCHMIDT H, KLEZL Z, et al. Stepwise reduction of functional spinal structures increase range of motion and change lordosis angle. J Biomech. 2007; 40(2):271-280.
[26] NACHEMSON AL, SCHULTZ AB, BERKSON MH. Mechanical properties of human lumbar spine motion segments. Influence of age, sex, disc level, and degeneration. Spine (Phila Pa 1976). 1979;4(1):1-8.
[27] CHEN YC, TU YK, ZHUANG JY, et al. Evaluation of the parameters affecting bone temperature during drilling using a three-dimensional dynamic elastoplastic finite element model. Med Biol Eng Comput. 2017;55(11):1949-1957.
[28] 朱康平,祝建雯,曲恒磊.国外生物医用钛合金的发展现状[J].稀有金属材料与工程,2012,41(11):2058-2063.

引言

随着计算机技术的快速发展，有限元法已经成为研究生物力学最有效的方法之一[1]，能够针对复杂结构及系统进行应力应变分析[2]，如椎体有限元模型的建立[3]、椎弓根内固定进行分析[4]、椎间盘突出问题等脊柱生物力学分析[5-8]。脊柱稳定由椎体、椎间盘、韧带及肌肉等组织共同维持，腰骶椎-骨盆承受着身体上部质量，椎体L4、L5处于脊柱的最下端，是传递载荷重要的枢纽，因此对于腰椎L4、L5的生物力学分析显得至关重要。椎体由外到内由皮质骨和松质骨构成，由于骨质材料属性的各向异性及骨质差异性较大，为提高计算效率，学者们对于椎体材料属性赋予方式也略有不同。2018年，文毅等[8]运用CT薄层扫描技术成功建立了可用于腰椎体及退变椎间盘力学分析的有限元模型，采用均分法将椎体体网格灰度值均分为10种材料属性，通过增加模拟体质量进行载荷的加载，分析椎间盘不同部位von-Mises应力变化情况。2019年，薛志鹏等[9]基于CT灰度值-表观密度-弹性模量的转化公式进行材料属性赋予，集中对比了20例正常股骨头与股骨头坏死的有限元模型，多样本的正常股骨头与股骨头坏死比较显示灰度值赋值具有较好的准确性和稳定性。2019年，罗林聪等[10]对肱骨进行皮质骨、松质骨两种材料属性赋值方式，并进行有限元模型与AnyBody骨骼肌肉的多体动力学结合分析，实现了肘关节屈曲运动过程中肱骨应力、位移的一个连续动态的分析，更符合人体生理解剖要求。2020年，张聪等[11]基于CT图像数据建立了青少年特发性脊柱侧凸患者腰椎-骨盆三维有限元模型，对椎体采用皮质骨、松质骨两种均质材料属性进行赋值，通过不同工况的加载分析表明，不同载荷会引起椎体或椎间盘不同位置受力及位移的相应变化。
不同的椎体材料属性赋予方式对有限元计算结果有很大影响，因此通过有限元系统对比椎体不同材料属性赋值方法显得十分重要。该文将基于逆向工程原理建立椎体模型，对椎体材料属性赋予分别采用灰度值赋值方法、均质及灰度值混合赋值方法、皮质骨-松质骨赋值方法，进行腰椎静力学对比分析，观察腰椎在3种赋值方法下不同工况的应力应变情况。
中国组织工程研究杂志出版内容重点：生物材料；骨生物材料; 口腔生物材料; 纳米材料; 缓释材料; 材料相容性；组织工程

材料方法

1.1 设计生物力学数值计算分析。
1.2 时间及地点实验于2021年在大连交通大学人机工程及医疗器械研究实验室完成。
1.3 材料选取22岁健康男性志愿者的一组脊柱CT数据，层厚1.0 mm，像素512×512，通过DICOM 3.0格式保存。受试者对实验知情同意。
1.4 实验方法

1.4.1 椎体提取及椎间盘的建立 Mimics是一款基于阈值的快速分割软件，通过查阅文献[12-15]，该文采用灰度值范围为226-1 753 Hu进行椎体提取。由于CT精度不同及提取过程中的人为不当操作，导致椎体表面不光滑，严重的会出现较大孔洞，影响有限元网格划分及静力学分析，因此使用Geomagic wrap软件进行表面处理及优化，效果见图1所示。通过3-matic软件进行椎间盘的构建，建立上下软骨终板、纤维环、髓核，此方法建立的椎间盘上下表面与上下椎体端面贴合性较好。

1.4.2 三种椎体材料属性赋予方法对优化后的椎体及椎间盘进行曲面重构生成实体，通过Hypermesh软件进行网格划分，单元网格质量检查标准如下：翘曲度<15°，单元长度比≤5°，三角形网格单位最小角≥30°，三角形网格单位最大角≤120°，雅克比> 0.6，网格锥度<0.6。为使静力学分析对比更加准确，3种方法赋值采用同一个网格划分的腰椎文件。
A计划——灰度值赋值方法：研究表明骨密度与弹性模量之间存在一定的线性关系，这也为区分椎体不同位置的材料特性提供了依据[16]。水的灰度值为0，脂肪灰度值为-90至-70，实验将使用Mimics软件通过灰度值Hu的计算进行梯度弹性模量的赋值，CT数据灰度值Hu与椎骨弹性模量之间的关系如公式(1)(2)所示[17]。

其中ρ为CT图像某一点的表观密度，单位为g/cm3；Hu为椎体的灰度值；E为弹性模量，单位为MPa。
通过Mimics软件的材料赋值模块可以得到CT数据上每一点的灰度值，通过公式计算从而对每个网格单元进行材料属性的赋予，由椎体外部向内密度及弹性模量值逐渐递减，弹性模量赋值结果如图2所示。

B计划——均质及灰度值混合法：皮质骨分布在骨的表面，致密、具有较高的硬度，椎体前端皮质骨分布较均匀，松质骨多孔结构表现出的各向异性比皮质骨更加明显。椎体的后部结构包括椎弓根、横突和棘突部分，材料属性很复杂，但其皮质骨的厚度更薄，因此后部结构可视为不区分皮质骨和松质骨[18]。松质骨密度与弹性模量存在着一定的线性关系，所以建立了图3中的模型，将椎体分割为3个部分，前端皮质骨、前端松质骨、椎体后部，各个部位材料属性赋予如图3所示，椎体前端皮质骨采取单一弹性模量赋值为 12 000 MPa，椎体前端松质骨及椎体后部利用Mimics软件通过灰度值进行梯度赋值。

C计划——皮质骨、松质骨赋值法：骨是一种复合材料，骨质分为皮质骨及松质骨，皮质骨分布在椎体表层，抗压缩、抗弯曲能力较强；松质骨结构复杂，由骨小梁交织组成。学者们对椎体材料属性赋予时采用皮-松质骨两种材料[19-20]，该文第3种方法进行1 mm厚皮质骨的建模，其余内部为松质骨，皮质骨弹性模量为12 000 MPa，松质骨弹性模量为400 MPa，弹性模量赋予效果如图4所示。

1.4.3 椎间盘材料属性赋予椎间盘是椎体间的重要连接，椎间盘富有弹性，抗压、减缓冲击。椎间盘纤维环包含纤维环纤维和纤维环基质，纤维环纤维弹性模量较纤维环基质大很多，研究显示纤维环可以用单一的材料属性来进行赋予。对于该文中椎间盘的材料通过查阅文献[19-21]，材料属性赋予如表1所示。

1.4.4 韧带设置及材料属性赋予根据解剖学结构对韧带进行建立，添加前纵韧带、后纵韧带、棘间韧带、棘上韧带以及横突韧带和黄韧带，韧带结构参数如表2所示[19]。椎体上下端面和软骨终板上下面的接触定义为绑定，相邻小关节面之间的接触定义为摩擦接触，摩擦系数为0.1。

1.4.5 外部载荷加载设置
有效性验证设置：①边界条件：约束椎体L5下表面6个方向上的自由度；②施加载荷：在L4椎体上端面耦合点分别施加竖直向下的轴向载荷100，200，300，400，500 N。
腰段L4-L5椎体6种工况的设置：①边界条件：约束椎体L5下表面6个方向上的自由度；②施加载荷：通过人体体质量分布计算所得[22]，在L4腰椎上端面耦合点施加300 N集中力模拟L4上端面承受的重力，同时分别在耦合点X轴正负向、Y轴正负向、Z轴正负向施加2，4，6，8，10 N·m力矩来模拟腰椎前屈、后伸、左右侧弯、左右旋转6种工况[19，23]，如图5所示。

1.5 主要观察指标不同工况下各组有限元模型的椎体应力应变情况。

讨论

有限元分析在脊柱研究中的应用较为广泛[27]，数值模拟计算不会对人体造成伤害，尤其是对于研究骨组织生物力学特性。正确的弹性模量赋值能够使有限元模型更加接近于身体的真实力学特性。对于植入体，不合适的弹性模量会产生应力遮挡等弊端，因此通过改变孔隙率控制植入材料弹性模量的多孔钛合金被研制出来，孔径适宜的孔隙也有利于骨细胞的黏附、生长和增殖[28]，所以材料属性的正确赋予显得十分重要。椎体材料属性表现为各向异性，针对椎体力学性能，该文基于逆向工程原理建立精度较高的椎体模型，还原了椎间盘结构，对L4、L5椎体进行了3种方法的材料属性赋予，并对3种材料属性赋值的椎体进行多种工况下模拟计算。
从轴向载荷加载可以看出，3种材料属性赋予的椎体轴向位移变化在BROWN等[24]实验结果范围内，方法A在500 N轴向加载时比BROWN实验的最大位移大0.135 mm。从实验结果可以看出，3种方法随着载荷加载呈现出相似的线性递增，与HEUER等[25]和NACHEMSON等[26]中的真实椎体偏转角度相似，通过与实验数据的角度偏移进行均方根误差分析，A计划方法中出现了更多次的最小均方根误差，均方根误差显示接近于实验值的程度，通过打分机制进行排名，结果总分A计划为27分，B计划为20分，C计划为25分，可以得出在多种工况下A计划方法相比于其他两种方法更接近于HEUER等[25]和NACHEMSON等[26]的实验数据。在Mises应力云图分布上，A计划方法明显表现出较B、C计划的集中应力范围较小，以及在皮质骨到松质骨材料赋予交界面不会产生应力值的大幅度递减。在椎体最大应力值对比中可以发现，相比于A、C计划，B计划中明显出现了更多次的应力集中现象。综上所述，相比于其他两种赋值方法，采用A计划方法(即梯度法)进行材料属性的赋予更接近于真实的椎骨生物力学特性。
相比于实验的昂贵成本及耗时较长，合理利用有限元计算能够满足一定的需求，可为临床医生使用工程软件提供思路。由于人体个体之间表现出个性化和较大的差异性，导致有限元结果也存在偏差，后续应进行同一椎体的实验结合有限元进行对比，此3种材料属性赋值方法作为一种参考和对比。
中国组织工程研究杂志出版内容重点：生物材料；骨生物材料; 口腔生物材料; 纳米材料; 缓释材料; 材料相容性；组织工程