流体动力学 格划分技术在骨骼有限元建模中的应用

邓达人1，孟春玲1，冯敏山2，张刚1

(1．北京工商大学材料与机械工程学院，北京100048;

2．中国中医研究院望京医院，中医正骨技术北京市重点实验室，北京100102)

摘要:

目的：研究一种高质量、高效率的骨骼分 方式。

方法：结合骨结构的特点，将流体动力学(computational fluid dynamics，CFD)的分 方法引入到骨生物力学的 格划分中，用六面体单元模拟皮质骨，四面体单元模拟松质骨。

结果：采用CFD 格划分技术能够划分出高质量的六面体单元，并且较好地模拟骨骼的结构特征，可以实现计算机自动分 ，分 花费时间约为传统分 方法的1/5，有限元模型计算结果基本符合尸体实验测量值。

结论：CFD 格划分技术可应用于骨生物力学领域中，为复杂人体骨骼的重构探索了一条有效途径。

关键词:骨生物力学;皮质骨;有限元模型;流体动力学; 格划分

在人体生物力学分析中，由于人体的解剖结构十分复杂，且构成人体的三维实体表面往往是不规则的自由曲面，故如何建立合理有效的有限元模型一直是一项非常重要而烦琐的工作［1］。 格单元质量的好坏将直接决定有限元的计算结果， 格质量可能导致应力计算结果偏差达到7%～100%［2］，低质量的 格甚至会导致有限元计算无法完成。

自1972年现代生物力学之父冯元桢首次将有限元法应用到生物力学研究以来，各国研究者先后采用了多种 格划分方法，其中主要有四面体 格自动划分法和六面体映射 格划分法［3-5］。田卫军等［6］基于逆向技术制作骨盆环四面体有限元模型;曹立波等［7］建立及验证儿童假人头部有限元六面体模型。作为最简单的三维单元，四面体单元是常应变、常应力单元，计算精度较低，而且对于复杂的接触条件常伴随计算不收敛，无法满足复杂载荷条件下的力学分析;当三维实体表面是十分复杂的自由曲面时，若全部采用六面体单元，人工分区比较困难，不得不对外观做较多简化， 格划分耗时也较长，甚至在曲率变化复杂的畸形表面难以实现。

骨结构分为松质骨和皮质骨，皮质骨硬度和密度均很高，是骨的主要受力部分。同时，骨骼之间的非线性接触问题一直是骨生物力学研究的重点，皮质骨表面应力梯度较大，故皮质骨需要采用精度比较高的 格。本研究结合骨结构的特点，将流体动力学(computer fluid dynamics，CFD)流体边界层的 格划分技术［8］，引入到骨生物力学的骨骼 格划分中。把皮质骨的外表面划分为二维四边形 格，以该 格做参考 格向模型内部映射出若干层三维六面体单元，用六面体单元模拟皮质骨;对于受力较小的松质骨采用四面体单元。二维四边形 格划分曲面时不用考虑人工分区和 格映射，受到曲面曲率变化的影响较小，无需对模型外观做太多简化［9］;CFD 格划分法吸收四面体单元分 简单快速和六面体单元计算精度高的特点。研究工作为复杂人体骨骼的重构探索了一条有效途径，为人体骨骼的生物力学分析奠定基础。

1 骨模型传统 格划分方法

1． 1 四面体 格自动划分法

四面体 格划分方法是通过在所选区域内进行布点，然后将它们联接成四面体，该方法是目前最流行的实体 格划分方法之一，几乎在所有的有限元前处理软件中都可以实现该方法的 格划分［10］。但是在 格三角化过程中，使用该方法对所生成的单元形状难以控制， 格计算精度相对较低。如果通过增加 格密度提高精度又会导致计算量大大增加，对研究骨骼之间的接触等复杂问题难以完成计算［11］。例如:田卫军等［6］基于逆向技术的骨盆环三维重构与结构分析，其有限元模型采用四面体单元，椎体与椎间盘、耻骨联合处以及小关节之间采用GLUE连接，将各个部分的模型粘连成一个整体，无法研究骨骼之间的接触问题［12］。

1． 2 六面体映射 格划分法

六面体 格划分一般采用映射 格划分方法。该方法通过适当的映射函数将待剖分物理域映射到参数空间中形成规则参数区域，然后对规则参数区域进行 格剖分，最后将参数域的 格反向映射回物理空间，从而得到物理域的有限元 格［9］。目前常用的 格划分方法中，六面体映射 格划分法的算法简单、计算速度快、单元质量好、单元密度可控，因此在各个领域得到了广泛的应用。

骨骼之间的接触问题一直是骨生物力学研究的重点，由于该接触问题为非线性问题，皮质骨表面应力梯度较大，故皮质骨需要采用精度比较高的 格;目前而言，六面体单元精度比四面体单元高，因此，皮质骨需要划分为六面体单元，但是六面体分 无法实现类似于四面体的自动分 ，当三维实体的表面曲率变化十分复杂时，六面体映射 格划分法需要创建大量的辅助线，剖分模型 格映射的区域十分困难，故不得不对模型外观做较多的简化，划分高质量的六面体单元需耗费大量的时间［13］，甚至在曲率变化复杂的畸形表面难以实现。特别在研究畸形骨骼模型时，需要对 格模型进行反复的更改，导致分 的工作量非常大。

2 CFD 格划分法

2． 1 CFD 流体边界层分

骨是一种具有独特结构的高密度结缔组织，大体结构分为松质骨(骨松质)和皮质骨(骨密质)，前者分布于骨的内部占全身骨量20%，密度低于皮质骨且富有弹性;后者位于骨外侧，占全身骨量80%，其硬度和密度均很高，主要用于发挥骨的支撑作用，是骨生物力学研究的重点之一［14］。

本研究根据骨结构组成及特点，将CFD中流体边界层的 格划分技术引入到骨生物力学中。CFD 格划分法是利用大型有限元前处理软件Hyper Mesh的CFD-Tetramesh功能，将流体管道模型的外表面划分为二维四边形 格，以该二维 格做参考 格，向模型内部映射出若干层三维六面体单元作为流体边界，再将模型内部划分为四面体单元作为流体域［9］。因此，在对骨骼进行有限元分 时，可以使用CFD 格划分法，将皮质骨的部分采用精度高、不易分 的六面体单元，将松质骨的部分采用精度较低、容易分 的四面体单元。

2． 2 CFD 格划分法的分 过程

以某医院提供的枕寰枢复合体模型为研究对象，采集CT断层图片，应用Mimics医学影像处理软件完成枕寰枢复合体的三维重建并获取点云数据［见图1(a)］。

将枕寰枢复合体的点云数据导入三维建模软件SolidWorks，进行逆向建模，得出骨骼三维实体模型［见图1(b)］。

使用大型有限元前处理软件HyperMesh对骨骼的三维模型进行CFD 格划分，将骨骼模型的几何面划分为二维四边形 格见图1(c)］;

以二维四边形 格为参考，向模型内部映射出3层高质量的六面体单元作为皮质骨［见图1(d)］;

再将模型其余部分划为四面体单元作为松质骨［见图1(e)］，

完成骨骼模型的 格划分［见图1(f)］。

CFD 格划分法极大减小了分区的难度，保证了模型内部的 格质量;该方法根据模型的几何面生成二维四边形 格，保留了原模型的几何特征，降低简化对模型精度的影响，由于二维四边形 格划分对象是曲面，划分曲面时无需考虑 格映射，曲面的复杂程度对二维四边形 格划分造成的影响较小，故采用二维四边形 格划分几何面比采用六面体单元划分几何实体的难度小;该方法能够在保持模型几何面形状特征不变的前提下，调整任意位置六面体单元层的厚度和层数，实现模拟皮质骨厚度不均的目的。图2所示为枢椎 格模型的上关节面剖视图，该关节面皮质骨的最大厚度被调整为正常厚度的1．5倍;此外，通过控制外部六面体层厚度、层数和 格密度以及内部四面体的 格生长率，能够再次调整骨骼模型的 格质量和 格密度。

采用CFD 格划分法对枕寰枢复合体模型的枕骨、寰椎、枢椎等复杂的三维实体进行 格划分，建立的最终模型如图3所示。

图3枕寰枢复合体的 格模型

该模型保留了骨骼的几何特征，骨骼接触面 格尺寸为0.5mm，其余部分 格尺寸为1.5mm，皮质骨厚度为1mm，由3层尺寸为0.33mm的六面体单元组成;模型包含103497个六面体单元和617784个四面体单元，单元类型包括C3D20、C3D8、C3D6、C3D4等。模型的单元质量见表1。

3 格模型验证

3． 1 有限元模型的建立

将HyperMesh制作的枕寰枢复合体 格模型，通过大型有限元软件ABAQUS进行韧带添加、材料参数赋值、接触定义、设置载荷以及约束条件等，从而完成有限元模型的建立。

3． 2 模型验证

通过计算枕寰枢复合体有限元模型在4种生理载荷下(前屈、后伸、侧屈、轴向旋转)枕寰、寰枢之间的相对运动角度，并与Panjabi等［16］1．5N·m载荷下尸体试验数据进行比较，验证有限元模型的有效性。模型在4种生理载荷下计算的枕寰(C0～1)、寰枢(C1～2)间相对运动角度与尸体试验结果见表2。

在1．5N·m载荷作用下，侧屈运动中寰枕间偏转角度比尸体实验数据大10．3°;轴向旋转中，寰枕间比尸体实验数据大0．42°，本模型其他运动结果均位于尸体试验数据范围内。当1．5N·m外载荷全部施加于模型时，模型整体的最大旋转角度发生在轴向旋转运动。

通过比较有限元模型的计算结果与尸体试验数据，证明采用CFD 格划分法制作的骨骼有限元模型，具有较高的计算精度，能够模拟复杂生理运动下的骨骼运动，对上位颈椎的临床研究具有指导作用。

3． 3 格收敛性分析

在有限元模拟运算中，不同的 格尺寸会对有限元计算精度产生较大的影响， 格尺寸越小，计算精度越高;但是过小的 格会大大增加模型的单元数目和模型解算时间步，故在建模中需要综合考虑模拟精度和计算时间成本等要素［10］。

通过比较4组不同 格尺寸的枕寰枢复合体有限元模型在相同生理运动条件下(前屈、后伸)的计算结果，验证不同 格尺寸对计算结果的影响。有限元模型的 格尺寸见表 3。

将4组有限元模型的计算结果与Panjabi等［16］在1．5N·m载荷条件下的尸体试验数据及某医院提供的1．5N·m载荷下尸体试验数据进行比较，具体数据见表4。1～4 模型的收敛时间分别为27、6、8、3h。

通过表 4 可知:

(1)模型的 格尺寸越小，计算结果越接近实验值。通过1、2 模型可以发现，随着 格尺寸的减小，模型计算精度逐渐提高，但是随着 格的减小其尺寸变化对计算精度的影响逐渐降低。同时，随着 格尺寸的减小，模型计算的时间成倍增加，大大提高了计算的时间成本。

(2)由1、3、4 模型可以发现，模型计算结果受到接触面 格尺寸的影响较大。较小的接触面 格尺寸，能够显著提升模型计算精度［12］。

4 结论

根据骨结构的特点，把CFD流体边界层的 格划分技术引入到骨生物力学的骨骼 格划分中。 格模型外观逼真，较好模拟了骨骼的生理特征;该分 方法综合了四面体 格自动划分法和六面体映射 格划分方法的优点，降低了分 的难度，提高了模型的 格质量，保证了模型的计算精度;同时受到模型几何面曲率变化的影响较小，保留了原始模型的几何特征，降低了简化对模型精度的影响;此外，可以调节任意位置边界层厚度，满足皮质骨在骨骼不同部位厚度不一致的特点;分 后的模型通过控制外部六面体层厚度、层数和 格密度以及内部四面体的 格生长率，能够再次调整骨骼模型的 格质量和 格密度。研究工作为复杂人体骨骼的重构探索了一条有效途径，为人体骨骼的生物力学分析奠定基础。