基于离散元的种植体力学性能研究

姜 杨，俞经虎，钱善华

(1. 江南大学机械工程学院，江苏 无锡 214122； 2. 江南大学江苏省食品加工技术装备重点实验室，江苏 无锡 214122)

1 引言

随着人们的生活条件得到改善，口腔问题也越来越受到重视。种植牙因在临床应用中具有良好的修复效果，咀嚼效率高，美观度好的特点，取得了较高的认可度。很多种植产品已经投入市场，种植成功率也到90%以上。但根据临床资料显示，仍然有种植失败的情况发生，这在长期临床实践来说仍旧是无法避免的问题。

种植牙植入成功的标准在于植入物和骨组织是否产生良好的生物相容性。咀嚼过程中，种植体和骨组织间的接触应力过大会导致骨水平下降，造成种植失败。因此研究咀嚼过程中种植牙和骨组织的应力分布，改善应力分布情况，对提高种植成功率有重要意义。由于数字化设计的发展，越来越多的人使用有限元法(FEM)来研究种植牙的生物力学。Dantas等通过有限元法比较常规植入物和牙根模拟植入物对骨组织应力分布的影响，提出牙根模拟植入物周围应力分布更均匀。赵吉奎等通过有限元法对种植体的结果进行了分析，提出等腰梯形螺纹能更好的传递力。而在大多数研究中都将牙齿的咬合力简化成垂直和斜向的加载，而在实际咀嚼过程中，种植体受到的咬合力是随时间变化的，各部分之间受力也存在区别。Sun等提出了用二维离散元软件PFC模拟咀嚼过程，研究了第一次咬合情况下食物形态的变化，探讨了食物和牙齿间的作用关系。吉林大学学者利用离散元理论建立了食物咀嚼破碎模型，分析了咬合过程中食物的变化情况，提出食物破碎模型可以用于模拟食物的破碎过程。蔡柳等通过离散元和有限元耦合模拟了刮板输送机运输煤料的过程，为输送机的结构优化提供了依据。

以往多数学者在研究种植义齿应力分布时选择静载荷作为加载方式，但咀嚼行为是一个动态过程，与实际情况不符。基于以上研究，本文提出一种新的载荷加载方式，先利用 CT 扫描重建磨牙齿冠模型，应用 Solidworks 软件建立种植体、皮质骨和松质骨模型，用 Hypermesh 软件对三维模型进行网格前处理。利用三维离散元软件 EDEM 建立了食物离散元模型，并模拟口腔咀嚼过程，通过 Add-in 耦合接口将牙冠上的压力输入 ANSYS Workbench 进行力学分析。

2 DEM-FEM联合分析模型

2.1 离散元模型

离散元主要是用于建立食物模型。食物被咀嚼后发生破碎，因此本文将食物建立为球体粘结模型，采用 Hertz-Mindlin 接触模型，微颗粒的接触部分重叠产生粘结键，颗粒之间发生粘结，颗粒间有粘结力和力矩，随着时间步的增加力和力矩如式(1)所示

(1)

式中：，-分别为法向力，切向力；

，-分别为法向力矩，切向力矩；

，-分别为法向刚度，切向刚度；

，-分别为颗粒的法向和切向速度；

，-分别为法向和切向角速度；

-时间步；

在受到外力作用后，颗粒被挤压，当颗粒之间的应力超过临界值时，粘结键被破坏，颗粒之间发生相对运动。法向和切向应力的最大值如下式(2)

(2)

本文模型选用大豆，以子叶所在面为正视面，用游标卡尺测量大豆的三维尺寸，测10组，将大豆近似为球形，计算出等效直径D=6.95mm。将食物材料作为均质的且各向同性。对食物颗粒建立粘结模型，大豆的密度ρ由质量与体积的比值计算得出1.15g/cm，泊松比v取经验值0.4，剪切模量G由式(3)计算得出

(3)

为弹性模量，测量方法参考文献[10]通过质构仪实验测得，经过多次实验剪切模量取均值71Mpa。

2.2 有限元模型

种植牙模型由牙冠、基台和种植体三部分组成，用 Solidworks 软件进行建模，种植体模型参考瑞士ITI种植体，见图1。选用标准种植体，体部直径D=4.10mm，植体长度L=10.00mm，颈部选用常规颈，直径D=4.80mm，穿龈高度L=2.80mm，颌龈高度L=3.00mm。螺纹部分螺距p=1.25mm，螺高L=0.40mm，螺深D=0.35mm。骨组织模型仿照下颌横截面，表面覆盖一层厚为2mm的皮质骨，内部是松质骨。

由于 EDEM 目前只能读取三角形的面网格，而离散元软件 EDEM 和有限元软件 Work-bench 耦合成功的条件是从 EDEM 中导出的力可以正确加载到 Workbench 对应的单元上。因此使用 Hypermesh 软件对种植体模型进行网格前处理，先划分四面体网格，如图2所示。由于咀嚼过程中与食物发生接触的只有牙冠部分，因此只需要在牙冠体网格表面划分一层三角形面单元，在 EDEM 中只需导入牙冠部分的面网格。本次研究中各部分材料均为均质、连续、各向同性的线性材料，各种材料参数如下表1所示。

图1 种植体三维模型

图2 网格模型

表1 材料的力学性能参数

2.3 有限元和离散元联合分析

DEM 和 FEM 耦合分析过程可概括为：①应用 Hypermesh 对种植体进行前处理；②使用DEM 进行咀嚼过程模拟，输出牙冠表面应力分布情况；③通过 EDEM 与 Workbench 的耦合，将力加载到种植体的网格模型上；④使用 FEM进行进行种植体的应力分析。基于该方法，本文应用 EDEM 与 Workbench 两种软件耦合对种植体进行强度分析，计算咀嚼过程中种植体的应力分布和位移变形情况情况。

3 建立咀嚼粘结模型

3.1 颗粒接触参数

大豆之间的碰撞恢复系数通过跌落试验测得，如下图3，将大豆胶粘在测试平面上作为接触表面，大豆从高处落下，与测试平面接触碰撞弹起，通过高速摄像机记录大豆弹起高度，重复实验取平均值0.356。大豆之间的摩擦系数通过堆积角仿真获得，研究表明大豆间的摩擦系数对仿真结果影响不大，因此大豆间的动摩擦系数与静摩擦系数取等值0.4。关于牙齿与大豆间的接触参数查阅文献获得如表2所示。

图3 测试碰撞恢复系数

表2 颗粒接触参数

3.2 粘结参数选取

基于离散元理论，大豆可以分为很多粘连的离散小颗粒组成的，微颗粒粒径一致且相互粘结的，微颗粒间设置bonding接触模型。粘结接触模型需要设置5个参数：法向刚度，切向刚度，临界法向应力，临界切向应力和粘结半径。微颗粒由粘结键作用粘结，当食物与牙齿面之间发生接触碰撞时，微颗粒受到法向和切向的力作用，当相邻颗粒间的法向和切向应力超过临界点时，粘结键发生断裂，食物被破碎。

大豆的法向刚度通过食品物性检测仪获得，将大豆颗粒平放，检测探头沿垂直于大豆子叶平面的方向向下压缩，大豆发生破碎变形，重复实验10次取均值，大豆的法向刚度系数可用式(4)计算

=ΔΔ

(4)

式中：Δ-垂直于大豆子叶平面的载荷增量；

Δ-大豆产生的形变增量。

切向刚度系数一般为法向刚度系数的06-1倍，本文取08，临界法向应力和临界切向应力按以下式(5)计算：

(5)

-临界法向应力，单位；

-临界切向应力，单位；

-粘结半径，单位。

通过计算得出粘结参数取值如下：粘结半径r为0.6mm，单位面积法向刚度为67693N·m，单位面积切向刚度为54154N·m，临界法向应力为8.62×10Pa，临界切向应力为6.90×10Pa。

4 咀嚼过程模拟

4.1 离散元仿真

在 EDEM 软件中设置两种尺寸颗粒(大豆完整颗粒和破碎之后的微颗粒)，设置在0.3s时将大豆颗粒替换成微颗粒粘结成的颗粒团。设置磨牙运动方式时，将咬合运动简化为上下的单向匀速运动，即下牙保持静止，上牙向下移动。从上牙开始向下运动至回到原位为一个运动周期。上牙的运动速度设置为2mm/s，连续咀嚼3个周期。将颗粒工厂(用于生成颗粒)设置在上下磨牙之间，由于食物在口腔内的咀嚼过程是随机的，因此将颗粒的生成位置设置为随机。

大豆颗粒由整体变化成颗粒团后，颗粒间产生粘结键颗粒粘结在一起，变化过程如图4(a)-(b)所示。此时上牙向下运动咀嚼大豆，粘结键被破坏断裂，颗粒团发生分散，随着上牙位移增加，断裂键增多，如图4(b)-(d)所示。

模拟过程结束后，通过 EDEM 后处理模块输出下牙在咀嚼过程中的载荷数据文件。图5为咀嚼3个周期的压力随时间变化的曲线。

图4 咀嚼过程模拟

图5 仿真结果

4.2 仿真结果验证

利用质构仪对大豆进行破碎，将压缩探头作为上磨牙，按照仿真案例中的运动方式进行压缩，如图6所示，每隔0.1s输出一次数据，得到大豆受到的压力变化曲线如图7所示。从图7可以看出，质构仪连续压缩3次的压力变化曲线趋势和离散元仿真结果变化趋势基本一致。应用 SP-SS 软件对质构仪实验和仿真得到的压力值进行一元线性回归分析，分析结果得判定系数R为0.941，拟合优度为0.939，显著性为0.000，实验值y和仿真值x之间的关系如下式，实验值与仿真值对比如图8所示。

=000007361-0001

(6)

图6 实验过程

图7 实验结果

图8 结果对比

经过重复实验，结果未出现显著性差异，表明通过 EDEM 仿真得出的结果在回归模型修正后可用于分析咀嚼食物过程中种植牙的受力情况。

5 有限元分析

利用 Add-in 耦合模块将Ansys Workbench中的 EDEM 模块和 Static Structural 模块进行耦合，通过 Add-in 模块可以把颗粒作用在几何设备上的压力数据传输到有限元分析软件。如图9所示，A模块读取 EDEM 导出的数据文件，B模块用于导入 Hypermesh 前处理后的种植体模型，C模块用于分析种植牙整体受力情况。

图9 耦合流程

由于种植体是一个装配体，各组件之间存在相互的作用关系，需要设置接触关系来传递载荷。牙齿和基台用粘结剂固定，种植体和骨组织已经实现骨整合，所以将上述接触类型都设置为 Bonded 接触。EDEM 中提取的压力分布在牙冠表面上，因此导入载荷时加载位置应选择整个牙冠，如图10所示。

图10 压力输入

从压力载荷分布来看，模拟咀嚼过程中牙冠最大受力处在牙尖部位，实际咬合过程中牙尖首先接触到食物。由于在模拟咀嚼过程中，食物颗粒产生位置是不固定的，因此设置3组模拟实验进行对照。

对于种植体的有限元分析结果，将 Von-Mises 应力作为比较依据。从3组有限元分析结果可以看出：有限元模型的应力分布模式是一致的，等效应力的峰值存在较小的差别。最大应力集中出现在牙冠的牙尖部位，也就是压力载荷施加的位置，种植体上应力分布在颈部区域，皮质骨上应力分布主要分布在与牙冠接触面上，松质骨上整体应力分布程度都比较低。

如图11所示，3组模型中最大应力都集中在牙尖处，最大应力值为155Mpa。如图12所示，3组模型中种植体的最大应力值出现在中段颈部部分，与骨组织相接处，最大应力值为60Mpa，出现在近皮质骨处，靠近第一道螺纹。如图13所示，3组模型中骨组织的应力集中在皮质骨周围，与种植体接触区域，并呈环状分布，最大应力值为19Mpa。由此可以看出，单次作用于种植牙的压力载荷产生的应力未超过各部件材料的破坏强度，也小于材料的屈服强度，所以正常咀嚼情况下，静载荷不会引发种植体的断裂失效。由此可以推测种植牙的失效可能是经历了长时间不间断的作用，导致材料超出本身的疲劳极限，出现突然断裂的情况。

图11 整体应力分布

图12 种植体应力分布

图13 骨组织应力分布

种植体能够长期使用的前提条件在于骨结合成功，骨结合取决于骨应力分布及大小。骨应力过大会导致骨组织融合失败，出现种植牙松动的情况。之前有研究证明，在使用过程中，种植体上的应力集中分布在种植体与皮质骨接触的区域，应力最大值出现在种植体上部螺纹处。本文的有限元结果也证明了以上结论：种植体的应力集中多数情况下都集中在上部区域，与骨组织接触部分。这也表明压应力可能就是造成骨结合失败并造成种植体疲劳失效的原因。当种植体承受的应力过大时，会破坏种植体和骨组织之间的融合，造成种植失败，因此考虑到种植牙的长期使用要求，要优化种植体的结构，改善应力集中情况，延长种植体使用寿命。

6 结论

本文以大豆为例，利用离散元方法模拟了口腔咀嚼过程，分析了咀嚼过程中下磨牙的受力情况，用质构仪进行大豆破碎实验，对仿真值进行了修正。通过ANSYS有限元软件分析得到种植牙模型的受力情况：

1)牙尖处的载荷大于周围的载荷，牙窝处压力载荷较小，应力集中一般出现在较高侧的牙尖部位，长期使用过程中会出现磨损情况。

2)种植体的应力分布在颈部区域，应力最大值出现与皮质骨结合处，靠近第一道螺纹处，在咀嚼较硬食物时有折断风险。

3)骨组织应力集中在皮质骨和松质骨交界处，呈环状分布，皮质骨承受较大的应力，长期会导致骨组织吸收。

本文在研究过程中对有限元模型中的种植体和下颌骨结构进行了简化，存在不足之处。咀嚼不同食物时，牙齿承受的咬合力差异较大，本文仅针对大豆进行了模拟分析，后期应考虑对不同特性的食物进行模拟，更深入的进行种植牙力学性能研究。