方法文章

一种用于定位上颌牙齿阻抗中心的有限元方法

9.4K 次观看

⸱

DOI：

10.3791/60746

⸱

2020年4月8日

 ,  ,  ,  ,  , 
1Division of Orthodontics, University of Connecticut Health, 2Private Practice, Miami, FL, 3Department of Biomedical Engineering, University of Connecticut, 4Department of Mechanical Engineering, University of Connecticut

通讯作者： Madhur Upadhyay <maupadhyay@uchc.edu>

本文内容

摘要

本研究概述了利用上颌及上颌牙齿的低剂量三维锥形束患者影像建立有限元模型所需的工具。随后，利用这些患者模型精确定位所有上颌牙齿的CRES。

摘要

阻力中心（CRES) 被视为可预测牙齿移动的基本参考点。用于估算 C 的方法RES 牙齿的研究方法从传统的放射学和物理测量，到对模型或尸体标本进行的体外分析不等。基于模型和单颗牙齿的高剂量微计算机断层扫描（micro-CT）图像进行有限元分析的技术已展现出广阔的应用前景，但针对新型低剂量、低分辨率锥形束计算机断层扫描（CBCT）图像的研究仍较为有限。此外，CRES 仅有少数特定牙齿（如上颌中切牙、尖牙和第一磨牙）的特征被描述过，其余牙齿则在很大程度上被忽视。此外，还有必要描述确定C点的方法学RES 详细说明，以便于复制和进一步拓展。

本研究利用常规的锥形束CT（CBCT）患者影像，开发用于获取上颌牙齿阻抗中心（CRES）有限元模型的工具与工作流程。通过对CBCT容积图像进行分割操作，提取与确定上颌牙齿CRES相关的三维（3D）生物结构。使用3matic软件对分割出的结构进行清理，并将其转换为由四面体（tet4）三角面构成的虚拟网格，最大边长为1 mm。随后，将这些模型进一步转化为最大边长为1 mm的实体四面体体积网格，用于有限元分析。采用工程软件Abaqus对模型进行预处理，建立装配体，并设定材料属性、相互作用条件、边界条件及载荷施加方式。在分析过程中，施加的载荷可模拟系统中的应力与应变，从而辅助定位CRES。本研究是实现牙齿移动精确预测的第一步。

引言

牙齿或牙段的阻抗中心（CRES）类似于自由物体的质心，这一术语借用于刚体力学领域。当单一力作用于阻抗中心（CRES）时，牙齿将沿力的作用线方向发生平移运动1,2。阻抗中心（CRES）的位置不仅取决于牙齿的解剖结构和物理特性，还受其周围环境（例如牙周膜、邻近骨组织、邻牙）的影响。牙齿属于受约束的物体，因此其阻抗中心（CRES）与自由物体的质心类似。在矫治器的操作过程中，大多数正畸医生会考虑力矢量与单个牙齿或牙组阻抗中心（CRES）之间的关系。实际上，当物体受到单一力作用时，其是否发生倾斜移动或整体移动，主要取决于该物体阻抗中心（CRES）的位置以及力矢量与阻抗中心（CRES）之间的距离。若能准确预测这一关系，治疗效果将显著提升。因此，对阻抗中心（CRES）进行精确估算可极大提高正畸牙齿移动的效率。

几十年来，正畸领域一直在重新审视关于特定牙齿、牙段或牙弓的CRES位置的研究1,2,3,4,5,6,7,8,9,10,11,12。然而，这些研究在方法学上存在诸多局限性。大多数研究仅针对少数牙齿确定了CRES位置，而忽略了大多数牙齿。例如，上颌中切牙和上颌切牙区段已被广泛研究。相比之下，关于上颌尖牙和第一磨牙的研究较少，而对其余牙齿则尚无相关研究。此外，许多研究基于牙齿的通用解剖数据、二维（2D）X线片测量结果以及2D绘图的计算来确定CRES的位置8。另外，部分现有文献使用通用模型或牙模的三维（3D）扫描数据，而非真实的人体数据4,8。随着正畸学逐步转向采用3D技术进行牙齿移动规划，重新审视这一概念至关重要，以建立对牙齿移动的三维、科学理解。

随着技术进步带来计算能力和建模能力的提升，创建和研究更复杂模型的能力也随之增强。计算机断层扫描和锥形束计算机断层扫描（CBCT）的引入，使模型和计算从二维世界迈向三维。计算能力与软件复杂性的同步提升，使研究人员能够利用三维放射影像提取精确的解剖模型，并在高级软件中用于分割牙齿、骨骼、牙周膜（PDL）以及各种其他结构7,8,9,10,13,14,15。这些分割出的结构可转化为虚拟网格，供工程软件使用，以计算在施加特定力或位移时系统的响应。

本研究提出了一种特定且可重复的方法学，可用于分析作用于从活体患者CBCT图像提取的模型上的假设性正畸力系统。采用该方法学，研究人员可估算不同牙齿的阻力中心（CRES），并综合考虑牙体结构的生物学形态，例如牙齿解剖结构、根数及其在三维空间中的方向、质量分布以及牙周附着结构。该过程的总体流程如图1所示，旨在帮助读者理解生成三维牙体模型以定位阻力中心（CRES）所涉及的逻辑步骤。

方案

已获得机构审查委员会豁免，用于评估口腔颌面放射学部存档的锥形束CT（CBCT）数据集（IRB编号：17-071S-2）。

1. 体积选择与标准

  1. 获取头面部的锥形束CT（CBCT）图像16。
  2. 检查图像的牙齿排列情况、缺牙状况、体素大小、视野范围以及图像的整体质量。
  3. 确保体素大小不超过 350 µm（0.35 mm）。

2. 牙齿和骨骼的分割

  1. 将CBCT图像的原始DICOM文件导入Mimics软件进行分割（图2）。点击图像 > 裁剪项目，裁剪图像以仅包含上颌骨及上颌牙齿。
    注意：视野应足够大，以完整包含上颌骨和上颌牙齿。确保图像包括牙冠、硬腭至鼻底、上颌窦、上颌牙齿的面部表面，以及硬腭后缘和上颌结节。
  2. 右键单击掩膜标签，为图像创建一个新掩膜。根据目标牙齿，将掩膜重命名为左侧UL1、UL2、…、UL7，右侧UR1、UR2、…、UR7。
  3. 在掩膜后的CBCT图像上识别目标牙齿（参见各视图）。使用清除掩膜工具擦除掩膜。由于牙齿与骨组织的灰度值相近，软件可能无法区分两者。
    注意：Mimics中的阈值工具无法将牙齿和骨组织分别分割，因此需要采用其他分割方法。
  4. 点击多切片编辑工具（Ctrl + M），选择视图（轴向、冠状或矢状），根据需要手动高亮（即绘制）部分切片。
    注意：高亮更多切片可为结构增加更多细节。
  5. 点击插值工具，填充被跳过的切片之间的体积，并应用。
  6. 右键单击掩膜，选择计算3D体积的选项，生成该牙齿的3D体积。
  7. 对上颌牙弓中的每一颗牙齿重复步骤2.2–2.6。
  8. 选择所有上颌3D牙齿（UL7–UR7），右键选择平滑。将平滑因子设为0.4，迭代次数设为4。
  9. 为分割上颌骨，右键单击掩膜标签，为图像创建一个新的掩膜。
  10. 从预设阈值集的下拉菜单中选择自定义，调整阈值以包含完整的上颌骨。应用阈值前，请务必勾选填充孔洞选项。
    注意：皮质骨中≤1 mm的小孔可接受，因其可在后续阶段轻松去除。
  11. 点击动态区域生长工具，以填充掩膜中可见的大孔洞。除勾选多层选项外，还需选择上颌骨掩膜作为该工具的目标。将最小值设为50，最大值设为150。按住Control键的同时，点击掩膜中未被高亮的皮质骨区域。
  12. 右键单击上颌骨掩膜，使用平滑掩膜功能。为获得最佳效果，重复此步骤3次。
  13. 右键单击掩膜，选择计算3D体积的选项，生成上颌骨的3D体积。
  14. 选择上颌3D骨结构，右键选择平滑。将平滑因子设为约0.4，迭代次数设为4。
  15. 选择上颌3D骨结构，右键选择包裹。将最小细节设为0.2 mm，间隙闭合距离设为1 mm，并勾选保护薄壁选项。点击确定。
  16. 将上颌3D骨结构重命名为"Maxilla"。

3. 清理与网格划分

  1. 选择3D对象并复制（Ctrl + C).
  2. 打开 3matic 软件，然后粘贴（Ctrl + V选定的3D对象将作为3D结构出现在3matic的对象树和工作区中图3).
  3. 点击 固定 从工具栏中选择标签并使用 平滑 选项。在 操作步骤 框选所需的 3D 对象或实体 应用 默认参数。
  4. 点击 完成 从工具栏中选择标签并使用 局部平滑 选项。在 操作步骤 框选所需的3D对象或实体。使用光标手动平滑所需区域。
  5. 复制牙齿。在对象树中选择所有牙齿，右键单击，然后选择 重复.
  6. 选择 所有重复牙齿，分组并命名文件夹 "第1组"原始样本将作为最终用于分析的牙齿。
  7. 对于第1组中的重复牙齿，单击 曲线模块 和 创建曲线 选项。对所有重复的牙齿，手动绘制一条围绕釉牙骨质界（CEJ）的曲线。
  8. 选择 曲线, 轮廓，以及 边界 实体在以下 平滑曲线 选项
  9. 将冠部和根部表面分别选择并分离为独立的部分 按曲线分割曲面 选项并左键单击三维对象以选择。
  10. 通过在釉牙骨质界（CEJ）处将牙齿分为牙根和牙冠，从牙齿的牙根结构中获取牙周膜（PDL）。
    1. 将第1组的3D对象（在步骤3.6中生成）复制为第2组。对于第2组，在对象树框中，单击 对象从表面列表中删除牙冠表面。对第2组中的所有对象执行此步骤。
    2. 对于第2组，点击 设计模块 > 中空. 应用所需的参数（表1).
    3. 点击 固定模块 > Fix Wizard单击各个部件，更新并按照给出的指示进行操作。
    4. 对所有部件重复步骤 3.10.3。将组 2 中的所有部件重命名 "UL1_PDL" 到 "UL7_PDL" 和 "UR1_PDL" 到 "UR7_PDL".
  11. 在第1组中，从对象树框中单击 对象从表面列表中删除根面。
  12. 选择 填充孔洞（常规） 选项并选择轮廓。单击 不良轮廓 和 应用整个空间将被填满。
  13. 选择 设计模块 > 局部偏移 并选择整个牙冠表面。勾选以下选项： 方向 （选择外部） 偏移距离 （选择 0.5），以及 缩短距离 （选择 2.0） 应用.
  14. 重复步骤 3.13。
  15. 对上颌牙弓的每一颗牙齿重复步骤 3.11–3.14。
  16. 重新网格化图3)
    1. 点击 重绘模块 > 创建非流形装配体 > 主要实体 > 上颌骨 从对象树中选择 3.4 节（原始牙齿）中所有对象的相交实体并选择 应用.
    2. 点击 Remesh 模块. 拆分非流形装配体。
    3. 使用来自第1组的所有对象作为相交实体，重复步骤3.16.1-3.16.2 应用.
    4. 作为可选步骤，仅在需要时选择 完成模块 > 修剪 > 实体 > 上颌骨. 选择多余的结构（即噪声）并 应用.
    5. 点击 固定模块 > 修复向导 > 上颌骨 > 更新按照给出的说明进行操作。
    6. 使用第2组中的所有对象作为相交实体，重复步骤3.16.1 应用.
    7. 点击 Remesh 模块 > 自适应重网格划分. 从 3.16.6 中选择所有相交的实体 应用.
    8. 点击 Remesh 模块 > 非流形分裂装配.
    9. 点击 重绘模块 > 创建非流形装配体 > 主要实体 > 单个物体 （PDL）从对象树的第2组中选择 相交实体 > 选择相应对象 来自步骤 3.4（对应于该牙齿类型）以及 应用.
    10. 点击 重网格模块 > 自适应重网格划分. 从 3.16.9 中选择相交的实体 应用.
    11. 点击 重绘模块 > 非流形分裂装配.
    12. 对每颗牙齿重复步骤 3.16.9–3.16.11。
  17. 点击 重绘模块 > 保持质量的三角面片简化在对象树中选择所有实体（即牙齿、牙周膜和上颌骨）并 应用.
  18. 点击 重绘模块 > 创建体积网格 > 选择实体选择 网格参数.
  19. 对所有实体（即牙齿、牙周膜和上颌骨）重复步骤 3.18。
  20. 手动将输入文件（.inp）从 3Matic 导出至 Abaqus图4).

4. 有限元分析

注意：所有自定义 Python 脚本均可在补充附件中找到。这些脚本是使用 Abaqus 中的宏管理器功能生成的。

  1. 预处理设置
    1. 打开 Abaqus 并选择 Standard Model。点击 File > Set the Work Directory > Select Location for File Storage。
    2. 点击 File > Run Script，然后选择 Model_setup_Part1.py。
    3. 在模型目录中指定加载 .inp 文件至 Abaqus 的文件路径。
    4. 点击 Models > Simulation > Parts > Maxilla > Surfaces。
    5. 在对话框中将表面命名为 "UL1 _socket"。
    6. 在 Select the Region of the Surface 下选择 By Angle，并输入角度值 "15"。
    7. 确保牙槽的所有区域均已被选中，完成后点击 Done。
    8. 对各个牙槽重复步骤 4.1.4–4.1.7。
    9. 点击 Models > Simulation > Parts，然后选择 UL1 > Surfaces。将表面命名为 "UL1"。
    10. 在 Select the Region of the Surface 下选择 "Individually"，从屏幕上选择牙齿，然后点击 Done。
    11. 对所有牙齿重复步骤 4.1.9–4.1.10。
    12. 点击 Models > Simulation > Parts，然后选择 UL1_PDL > Surfaces。将表面命名为 "UL1_PDL_inner"。
    13. 在 Select the Region of the Surface 下选择 By Angle，并输入角度值 "15"。
      注意：若在最终模拟过程中发现错误，应减小角度并重新选择表面。
    14. 确保牙周膜（PDL）的整个内表面区域均被选中，完成后点击 Done。
    15. 选择 UL1_PDL > Surfaces，将表面命名为 "UL1_PDL_outer"。
    16. 在 Select the Region of the Surface 下选择 By Angle，并输入角度值 "15"。
      注意：若在最终模拟过程中发现错误，应减小角度并重新选择表面。
    17. 确保牙周膜（PDL）的整个外表面区域均被选中，完成后点击 Done。
    18. 对所有牙周膜（PDL）重复步骤 4.1.13–4.1.19。
    19. 点击 File > Run Script，然后选择 Model_setup_Part2.py。
    20. 点击 Models > Simulation > BCs。将边界条件命名为 BC_all，然后选择 Step 为 Initial。在类别下选择 "Mechanical"，在 "Types of Selected Step" 下选择 "Displacement/Rotation"。点击 Continue。
    21. 在 Select Regions for the Boundary Condition 下选择 By Angle，并输入角度值 "15"。勾选 Create Set，并为14颗牙齿分别选择各自的牙槽。完成后点击 Done。
      注意：此操作有助于模拟牙齿的瞬时移动。
    22. 点击 Models > Simulation > Assembly > Sets > Create Set。将集合命名为 "U1_y_force"。
    23. 在 Select the Nodes for the Set 下选择 Individually。
      注意：在随机选择的牙冠节点上施加1牛顿的集中力，方向为正Y方向（模拟远中移动力）或正Z方向（模拟压低力）。
    24. 在上颌中切牙（U1）唇侧牙冠中心选择一个节点，然后点击 Done。
    25. 点击 Sets > Create Set。将集合命名为 "U1_z_force"。
    26. 重复步骤 4.1.23–4.1.24。
    27. 对所有牙齿重复步骤 4.1.22–4.1.26。
      注意：在为特定牙齿生成集合（如步骤 4.1.25）之前，需先对该牙齿执行 Instance > Resume 操作。
  2. 模型设置
    1. 点击 Models > Simulation > Assembly > Instances。选择 All Instances，然后点击 Resume。
    2. 点击 Tools > Query > Point/Node。在随机选择的中切牙中心选择一个节点，然后点击 Done。
    3. 在页面底部的命令中心中，复制步骤 4.2.2 所选节点的 X、Y 和 Z 坐标。
    4. 在垂直工具栏中选择 Translate Instance，并在屏幕上选择整个装配体（即所有实例），然后点击 Done。
    5. 在 Select a Start Point for the Translation Vector 框中，粘贴步骤 4.2.3 复制的坐标，或手动输入 X、Y 和 Z 值，然后点击 Enter。
    6. 在 Select an End Point for the Translation Vector or enter X,Y,Z 中输入坐标 "0.0"、"0.0" 和 "0.0"，然后点击 Enter。
    7. 对于 Position of Instance，点击 Ok。
    8. 点击 Tools > Query > Point/Node，选择位于中切牙中线正上方的一个节点，然后点击 Done。
    9. 在页面底部的命令中心中，复制步骤 4.2.8 所选节点的 X、Y 和 Z 坐标。
    10. 在垂直工具栏中选择 Translate Instance，并在屏幕上选择整个装配体（即所有实例），然后点击 Done。
    11. 将复制的坐标粘贴至 Select a Start Point for the Translation Vector - or Enter X,Y,Z 框中，然后点击 Enter。
    12. 在 Select an End Point for the Translation Vector - or enter X,Y,Z 中输入步骤 4.2.9 复制的坐标，并将 X 坐标修改为 0.0，然后点击 Enter。
    13. 对于 Position of Instance，点击 Ok。
    14. 点击 File > Run Script，然后选择 Model_setup_Part3.py。插入或修改材料属性。
    15. 点击 Models > Simulation > Materials，然后点击 Bone/PDL/Tooth。输入组织特异性属性。
    16. 点击 File > Run Script，然后选择 Functions.py。
  3. 模型处理
    1. 点击 File > Run Script，然后选择 Job_submission.py。
      注意：作业模块用于在模型上设置一个或多个操作，作业管理器则用于启动模型分析、显示进度并提示完成状态。
    2. 在标题为 Suppress All 的对话框中，根据约束条件输入牙齿的侧别（L 或 R）（在 Models > Simulation > Constraints 下）。点击 Ok。
    3. 在标题为 Job Submission 的对话框中输入 "Y"，以运行指定牙齿的分析。点击 Ok。
    4. 在标题为 Directions for Analysis 的对话框中输入 "Y"，以指定力的施加方向。点击 Ok。
  4. 用于 CRES 估算的后处理
    1. 选择 File > Run Script > Bulk_process.py。
    2. 在标题为 Analyze Multiple Jobs 的对话框中，为指定的牙齿输入 "Y"。点击 Ok。
    3. 在标题为 Directions for Analysis 的对话框中，为指定力的施加方向输入 "Y"。点击 Ok。
    4. 在标题为 Get Input 的对话框中，输入如命名实例所示的 Specific Tooth Number（例如 UL1 或 UL5 等）。点击 Ok。
    5. 在命令框中核对 Force About Point 和 Estimated Location 的坐标。若两者不一致，则重复步骤 4.3.1–4.4.4。
      注意：在运行每一步的作业后，将在 Abaqus 界面内运行一个用 Python 编写的用户自定义算法，用于分析反作用力系统及加载后产生的力矩。该算法会自动建议一个新的节点位置以施加载荷，使得力系统中的力矩接近零。该过程为迭代进行，直至找到或估算出施加力时能产生最接近零力矩的节点位置。该算法的详细说明见讨论部分。

结果

为了验证操作步骤部分（步骤2）中所述的分割和手动勾画方法，从一个干燥颅骨中取出上颌第一磨牙，并拍摄了锥形束CT（CBCT）图像。按照步骤2所述，使用图像处理与编辑软件Mimics对该牙齿进行手动勾画。随后进行网格划分，利用3matic软件对分割后的模型进行清理，并将其导入Abaqus进行分析。我们未发现牙齿的有限元模型与实验室中实际测量的牙齿在线性及体积测量值之间存在显著差异（补充文件4）。

为验证用户自定义算法在确定物体CRES时的有效性，在脚本创建初期使用了一个梁被包裹于鞘内的简化模型（图5A）。钢制包层被约束为三个方向的位移自由度，梁与鞘界面处的节点被绑定在一起。随机选择施力节点，并以迭代方式应用子程序，直至解收敛。在该简化模型中，长度为30单位、宽度为10单位的梁被包裹于鞘内。通过遵循既定算法及其计算步骤，预测了该模型梁的CRES（图5B）。结果与理论计算一致（见补充文件3）。因此，该用户自定义算法在此简化模型中得以建立并验证，随后被用于上颌牙齿CRES的测定。

表2列出了为各结构指定的材料属性。牙周膜（PDL）和骨组织材料属性建模方式的差异可能会影响牙齿阻抗中心（CRES）的最终位置。与纤维取向相关的PDL各向异性、泊松比的差异、加载模式及载荷大小也可能产生影响。PDL的材料属性根据Ogden模型设定为非线性超弹性（µ1 = 0.07277，α1 = 16.95703，D1 = 3 x 10-7）22,23。同时指定了各组织的密度：骨组织为1.85 g/cm3；牙齿为2.02 g/cm3；PDL为1 g/cm3（即水的密度，因为PDL主要由水组成）24,25。

为了标准化力矢量并确定 CRES 的位置，建立了一个笛卡尔坐标系（X-Y-Z），其方向定义如下：Y 轴（前后方向或唇腭向轴）沿腭中缝方向，后部为正方向；Z 轴为垂直方向（上下方向或𬌗龈向轴），模型的上部或龈向部分为正方向；X 轴为横向方向（颊舌向轴），颊侧部分为正方向（图6）。

该坐标系以两种方式应用：1）建立一个全局坐标系，其原点（O）位于中切牙唇侧表面之间，切牙乳头下方，且处于X-Y平面内平分切牙间宽度和磨牙间宽度的直线上；2）为每颗牙齿构建局部坐标系，每颗牙的原点标记为“R”。每颗牙齿的“R”点定义为牙冠颊面表面的几何中心。选择该位置是为了近似模拟正畸操作中施加矫治力时托槽可能放置的最近位置。代表性结果见图7。

相对于全局和局部坐标系的CRES位置如表3和表4所示。当沿Y和Z坐标施加力系统时，沿X坐标获得的CRES位置彼此不同（表5）。然而，其平均差异较小（0.88 ± 0.54 mm）。

CT扫描数据分析工作流程；分割模型处理、体网格创建、有限元设置
图1：设计流程图。定位CRES的三步工作流程。 请点击此处查看此图的放大版本。

在Mimics软件中进行三维下颌模型分析；显示CT扫描、分割和重建过程。
图2：Mimics软件界面布局，显示上颌牙齿在三个正交视图（X-Y-Z）及体积模型中的呈现。 请点击此处查看该图的放大版本。

显示上颌扫描、网格操作和组装选项的3D建模软件界面。
图3：利用3matic软件的非流形组装生成牙周膜（PDL）的步骤。 网格重划模块（A）创建非流形组装，（B）将上颌设为主要实体，（C）将PDL设为相交实体，（D）自适应网格重划，（E）分离上颌与PDL，（F）以PDL为主要实体、选定牙齿为相交实体，重复步骤B-F，（G）创建体网格。请点击此处查看该图的放大版本。

牙科模型的有限元分析；网格图；模拟结构应力分布。
图4：Abaqus 软件界面布局。 请点击此处查看该图的放大版本。

三维建模示意图、结构分析、有限元网格、工程仿真研究。
图5：钢梁的简化模型。（A）包裹在钢套管中的梁，用于测试所定义算法的准确性。（B）由所定义算法预测的包裹梁的CRES位置。请点击此处查看该图的放大版本。

用于三维平面分析中CR和R位置的牙科空间坐标图及方程。
图6：相对于全局原点（O）和每颗牙齿的局部原点（R）的CRES估算坐标系。 本图以右侧上颌第二前磨牙为例进行说明。该方法已应用于牙弓中每一颗牙齿的分析。请点击此处查看本图的放大版本。

3D 表面模型；7 种形状标记为 A-G；结构分析；生物学研究；示意图
图 7：上颌牙齿 CRES 的三维示意图。 (A) 中切牙。 (B) 侧切牙。 (C) 尖牙。 (D) 第一前磨牙。 (E) 第二前磨牙。 (F) 第一磨牙。 (G) 第二磨牙。 请点击此处查看此图的放大版本。

中空类型：内外两侧
距离0.2
最小细节：0.05
减小：已勾选
边界清理：已勾选
清理因子：1.1

表1：中空工具参数。

结构弹性模量 (MPa)泊松比密度 (g/cm3)
牙齿170000.32.02
骨170000.31.85
PDL0.05见正文1

表2：有限元模型的材料属性。

牙齿编号牙冠长度牙根长度xyz
UL125.215.13.411.012.9
UL226.016.88.813.214.3
UL329.119.515.118.015.6
UL423.815.718.421.510.6
UL524.818.220.928.210.1
UL622.016.425.838.711.6
UL721.415.027.443.211.4
UR124.914.6-4.610.813.2
UR226.316.7-9.913.013.6
UR330.921.1-15.617.714.2
UR422.916.7-19.021.99.2
UR523.416.7-21.129.48.8
UR622.216.3-23.939.69.8
UR720.815.9-21.747.010.4

表3：上颌牙齿的CRES相对于全局原点O的三维（X-Y-Z）位置。

牙齿编号牙冠长度牙根长度xyz
UL125.215.1-1.110.99.4
UL226.016.8-5.59.410.4
UL329.119.5-5.79.313.2
UL423.815.7-6.45.79.0
UL524.818.2-6.77.09.5
UL622.016.4-6.98.310.4
UL721.415.0-8.63.37.3
UR124.914.60.510.811.1
UR226.316.75.010.39.3
UR330.921.15.78.512.0
UR422.916.75.35.39.3
UR523.416.75.36.59.1
UR622.216.35.67.810.1
UR720.815.99.54.38.6

表4：C的三维（X-Y-Z）位置RES 上颌牙齿相对于每个牙齿的局部点R的CRES 正在评估中。 此处，R 为牙冠颊面的几何中心。

牙齿编号FyFZ差异
UL1-1.36-0.800.56
UL2-5.73-5.230.5
UL3-6.00-5.450.55
UL4-6.11-6.650.54
UL5-5.95-7.401.46
UL6-6.18-7.671.49
UR10.360.670.31
UR25.234.770.46
UR35.935.380.55
UR44.576.011.44
UR55.884.691.91
UR65.195.980.79

表5：当沿Y轴（Fy）和Z轴（Fz）施加力时，阻抗中心位置在X轴方向上的变化。

补充文件1：用于有限元分析（FEA）的算法Python脚本。 请点击此处查看该文件（右键单击可下载）。

补充文件2：力系分析概述。 请点击此处查看该文件（右键点击可下载）。

补充文件3：包覆在鞘管内的简单梁质心位置的理论估算。 请点击此处查看该文件（右键点击可下载）。

补充文件4：一颗拔除的上颌第一磨牙的有限元模型。 请点击此处查看该文件（右键单击可下载）。

讨论

本研究展示了一套工具，用于建立基于患者锥形束CT（CBCT）图像构建的上颌牙齿模型进行有限元分析（FEA）以确定其牙根表面中心（CRES）的一致性工作流程。对于临床医生而言，清晰且直观的上颌牙齿CRES分布图将是一种极为宝贵的临床工具，可用于规划牙齿移动并预测副作用。有限元法（FEM）于1973年引入牙科生物力学研究17，此后已被广泛应用于分析牙槽支持结构中的应力与应变场6,7,8,9,10,11,12。从工作流程中列出的步骤数量可以看出（图1），构建有限元模型是一项复杂的任务，因此必须对方法学中的某些方面进行简化。

首先，假设牙槽骨的吸收和沉积不发生，仅考虑牙齿在牙槽窝内的移动。这种类型的位移称为原发性4或瞬时牙移动18。研究发现，牙周膜（PDL）在瞬时牙位移过程中起着关键作用。为了定义牙齿移动过程中牙周膜的应力，可合理地将骨组织和牙齿视为刚性体15。因此，在本研究中，应力分布被限制在牙槽窝内。使用“创建边界条件”工具可让用户为模型设置边界条件或施加约束。选定的点被赋予零自由度，以确保模型在该区域保持刚性。由此，以往研究中用于计算骨组织变形及对变形后的牙槽骨实体单元重新划分网格的分析时间得以消除19,20。

其次，尝试将图像分辨率保持在中等水平。CBCT图像的体素大小为0.27 mm。这不仅使辐射剂量保持在最低水平，还降低了为四面体单元组装全局刚度矩阵时的计算负担。然而，其缺点在于CBCT分辨率不足以在扫描中准确且清晰地捕捉牙周膜（PDL）。这主要是因为PDL的平均厚度约为0.15 mm–0.38 mm（平均值：0.2 mm）21，而图像体素大小为0.27 mm。CBCT扫描的这一不足导致了两个问题：1）无法单独对PDL进行分割；2）由于骨组织与牙齿之间缺乏明显的灰度值变化，无法通过阈值分割法对骨组织和牙齿进行分割。因此，软件无法区分牙齿与骨组织，因为两者的灰度值相近。换句话说，Mimics无法将牙齿与骨组织分别分割。为此，开发了一种不同的分割方法。在尝试了多种工具（如Mimics中的区域生长或分割工具）后，确定最佳的牙齿分割方法是在CBCT的每一层切片上手动勾勒出牙齿结构。在此过程中，多层切片编辑工具提供了效率优势。用户无需手动标记每一层切片，而只需标记部分切片即可。因此，该方法成为牙齿分割的最佳选择，因其能够以一致的方式获得牙齿解剖结构的高质量图像，具有最高的准确性。

由于Mimics无法根据CBCT图像的低分辨率对牙周膜（PDL）进行分割，因此需要从牙齿的牙根结构生长出牙周膜。这要求在釉牙骨质界（CEJ）处将牙齿分为牙根和牙冠两部分。生成后，构建的牙周膜本质上是两个相互平行、相距0.2 mm的表面，其中一个表面与骨组织紧密接触，另一个则与牙根接触。在有限元分析中，必须将这两个表面连接在一起，以确保施加在牙齿上的载荷能够通过牙周膜传递至骨组织。工程软件会拒绝表面间距过大或过度交叉的模型，因为这会导致表面无法连接，从而使有限元分析（FEA）模型失效。

第三，所有模型表面均保持相对光滑，避免出现对整体模型分析无显著意义的微小表面形貌，例如颊侧皮质骨表面额外突出的骨性结构。解剖结构突起上的细微部分会使最终模型的网格在精细解剖区域的复杂部位减小单元尺寸，从而增加模型中的单元数量，导致有限元分析的计算量增大。

C 的位置RES 当力在Y方向和Z方向施加时，其差异体现在X方向上的位置不同。然而，该差异较小（表5）在临床上和统计学上均无显著意义。因此，C的位置RES 单向计算的结果可用于另一方向。先前的研究还表明，在三维评估中，C 的单个点RES 未观察到10,26,27因此，有研究建议，与其认为存在一个明确的CRES 更合适的术语可能是 "阻力半径"这种差异可归因于多种因素，例如根系形态、边界条件、材料特性以及载荷施加位置。

使用自定义算法分析力系统
用于确定牙齿阻抗中心（CRES）的数学概念、理论推导及计算机模拟此前已有详细描述27,28,29,30。为了分析施加不同载荷后产生的力系统，并预测牙齿的CRES位置，编写了一个自定义算法并在Abaqus中运行（参见补充代码文件）。该算法使用Python编写，以有限元分析软件的输出数据库（.odb文件）中的数据作为输入，对数据进行处理，并输出由施加载荷在系统中产生的力矩值。此外，该算法还能估算出使系统内产生较小力矩的节点位置。这使得用户可以以迭代方式运行模拟，直至估算结果收敛于单一位置。

该算法获取每个加载步骤中节点的坐标、各节点的总位移以及各节点处因施加载荷而产生的反作用力。将系统中每个节点处与原始载荷施加方向相同的反作用力以及方向相反的反作用力分别进行求和，以确定在模拟过程中作用于牙齿上的合力矢量。针对每个节点处的反作用力，计算其相对于载荷施加点的力矩，并以与反作用力相同的方式对这些力矩进行求和。因此，可计算出与原始载荷方向相同的合力矢量及其关于载荷施加点所产生的力矩，同时也可计算出反方向的合力矢量及其相应的力矩。由于系统处于静力平衡状态，所有力与力矩的总和为零。然而，通过此种方式对反作用力和力矩进行分解，可进一步计算出这些合力在系统中作为支点作用的有效作用位置，而这两个支点之间的中心点即为更接近CRES的载荷施加点的近似位置。

为了进行这些计算，需将所得力矩的大小除以其对应作用力的大小，从而得到从支点到力作用点的距离（R矢量）的大小。R矢量的方向通过力矩矢量与力矢量的叉积确定，其中所有矢量必须相互正交，单位矢量则通过叉积的大小进行归一化得到。将此前计算出的R矢量大小乘以R单位矢量，即可获得每个支点相对于初始力作用点在三维空间中的坐标总体估计值。这两个矢量之间的中点即为下一次迭代中下一个力作用点位置的估计值。更多详细信息见补充文档2。

当系统中产生的力矩总和接近零时，即可确定CRES的估计值。在本研究中，该值通过取计算所得力矩中最小的正值和最小的负值对应的X分量，并对二者取平均来确定。由于节点位置是随机生成的，且任意两个节点之间存在固有的距离（0.5 mm），因此很难找到恰好产生零力矩的位置（表5）。

局限性
尽管我们已尽最大努力，本研究仍存在一些局限性。首先，由于牙周膜（PDL）无法在锥形束CT（CBCT）上直接显示，因此无法单独进行分割，而是通过牙齿根面以0.2 mm的均匀厚度生成。有限元分析研究表明，均匀建模与非均匀建模相比会影响有限元分析（FEA）的结果，且非均匀建模更为优越30,31。其次，构建精确模型所需的步骤较多，耗时较长。这限制了模型的生成速度，从而限制了这些工具在临床中为每位患者制定个体化治疗方案的应用前景。此外，生成此类模型所需的软件价格昂贵，通常仅限于教育机构或大型企业才具备相应资源。再者，模型构建完成后，运行有限元分析还需要非常强大的计算能力。因此，在相关技术尚未广泛普及之前，该方法尚无法成为可行的治疗规划工具。

未来的研究应着重利用这些模型对上颌牙齿进行有限元分析，以确定牙弓及各牙组的阻抗中心（CRES），尤其是正畸治疗中常被移动的牙组，例如拔牙病例中的前牙段，或开合患者中用于压入的后牙段。一旦确定了这些模型的CRES，应基于更多的锥形束CT（CBCT）图像构建更多模型，以扩充现有数据。当积累足够数量的CRES位置数据后，可生成热图以显示CRES的总体位置分布，从而为临床医生提供极具价值的参考依据。

披露

作者无任何利益冲突需要披露。

致谢

作者谨此感谢 Charles Burstone 基金会奖项对本项目的支持。

材料

本文使用的材料清单
姓名公司目录编号评论
3-matic 软件Materialise, Leuven, Belgium.清洁与网格化处理
Abaqus/CAE 软件，版本 2017Dassault Systèmes Simulia Corp., Johnston, RI, USA.有限元分析
Mimics 软件，版本 17.0Materialise, Leuven, Belgium.牙齿与骨骼的分割

参考文献

  1. Smith, R. J., Burstone, C. J. Mechanics of tooth movement. American Journal of Orthodontics. 85 (4), 294-307 (1984).
  2. Christiansen, R. L., Burstone, C. J. Centers of rotation within the periodontal space. American Journal of Orthodontics. 55 (4), 353-369 (1969).
  3. Tanne, K., Nagataki, T., Inoue, Y., Sakuda, M., Burstone, C. J. Patterns of initial tooth displacements associated with various root lengths and alveolar bone heights. American Journal of Orthodontics and Dentofacial Orthopedics. 100 (1), 66-71 (1991).
  4. Burstone, C. J., Pryputniewicz, R. J. Holographic determination of centers of rotation produced by orthodontic forces. American Journal of Orthodontics. 77 (4), 396-409 (1980).
  5. Dermaut, L. R., Kleutghen, J. P., De Clerck, H. J. Experimental determination of the Cres of the upper first molar in a macerated, dry human skull submitted to horizontal headgear traction. American Journal of Orthodontics and Dentofacial Orthopedics. 90 (1), 29-36 (1986).
  6. Tanne, K., Sakuda, M., Burstone, C. J. Three-dimensional finite element analysis for stress in the periodontal tissue by orthodontic forces. American Journal of Orthodontics and Dentofacial Orthopedics. 92 (6), 499-505 (1987).
  7. Meyer, B. N., Chen, J., Katona, T. R. Does the Cres depend on the direction of tooth movement? American Journal of Orthodontics and Dentofacial Orthopedics. 137 (3), 354-361 (2010).
  8. Kojima, Y., Fukui, H. A finite element simulation of initial movement, orthodontic movement, and the centre of resistance of the maxillary teeth connected with an archwire. European Journal of Orthodontics. 36 (3), 255-261 (2014).
  9. Reimann, S., Keilig, L., Jäger, A., Bourauel, C. Biomechanical finite-element investigation of the position of the centre of resistance of the upper incisors. European Journal of Orthodontics. 29 (3), 219-224 (2007).
  10. Viecilli, R. F., Budiman, A., Burstone, C. J. Axes of resistance for tooth movement: Does the Cres exist in 3-dimensional space? American Journal of Orthodontics and Dentofacial Orthopedics. 143 (2), 163-172 (2013).
  11. Ammar, H. H., Ngan, P., Crout, R. J., Mucino, V. H., Mukdadi, O. M. Three-dimensional modeling and finite element analysis in treatment planning for orthodontic tooth movement. American Journal of Orthodontics and Dentofacial Orthopedics. 139 (1), 59-71 (2011).
  12. Sia, S., Koga, Y., Yoshida, N. Determining the center of resistance of maxillary anterior teeth subjected to retraction forces in sliding mechanics. An in vivo study. Angle Orthodontics. 77 (6), 999-1003 (2007).
  13. Cattaneo, P. M., Dalstra, M., Melsen, B. Moment-to-force ratio, center of rotation, and force level: a finite element study predicting their interdependency for simulated orthodontic loading regimens. American Journal of Orthodontics and Dentofacial Orthopedics. 133 (5), 681-689 (2008).
  14. Tominaga, J. Y., et al. Effect of play between bracket and archwire on anterior tooth movement in sliding mechanics: A three-dimensional finite element study. Journal of Dental Biomechanics. 3, 1758736012461269(2012).
  15. Cai, Y., Yang, X., He, B., Yao, J. Finite element method analysis of the periodontal ligament in mandibular canine movement with transparent tooth correction treatment. BMC Oral Health. 15 (106), (2015).
  16. Pauwels, R., Araki, K., Siewerdsen, J. H., Thongvigitmanee, S. S. Technical aspects of dental CBCT: state of the art. Dentomaxillofacial Radiology. 44 (1), 20140224(2015).
  17. Farah, J. W., Craig, R. G., Sikarskie, D. L. Photoelastic and finite element stress analysis of a restored axisymmetric first molar. Journal of Biomechanics. 6 (5), 511-520 (1973).
  18. van Driel, W. D., van Leeuwen, E. J., Von den Hoff, J. W., Maltha, J. C., Kuijpers-Jagtman, A. M. Time-dependent mechanical behavior of the periodontal ligament. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine. 214 (5), 497-504 (2000).
  19. Bourauel, C., et al. Simulation of orthodontic tooth movements. A comparison of numerical models. Journal of Orofacial Orthopedics. 60 (2), 136-151 (1999).
  20. Schneider, J., Geiger, M., Sander, F. G. Numerical experiments on longtime orthodontic tooth movement. American Journal of Orthodontics and Dentofacial Orthopedics. 121 (3), 257-265 (2002).
  21. Ten Cate, A. R. Oral histology, development, structure and function (5th ed). , St. Louis Mosby. (1998).
  22. McCormack, S. W., Witzel, U., Watson, P. J., Fagan, M. J., Gröning, F. The Biomechanical Function of Periodontal Ligament Fibres in Orthodontic Tooth Movement. PLoS One. 9 (7), e102387(2014).
  23. Huang, H., Tang, W., Yan, B., Wu, B., Cao, D. Mechanical responses of the periodontal ligament based on an exponential hyperelastic model: a combined experimental and finite element method. Computer Methods in Biomechanics and Biomedical Engineering. 19 (2), 188-198 (2016).
  24. Yang, J. A new device for measuring density of jaw bones. Dentomaxillofacial Radiology. 31 (5), 313-316 (2002).
  25. Gradl, R., et al. Mass density measurement of mineralized tissue with grating-based X-ray phase tomography. PLoS One. 11 (12), e01677979(2016).
  26. Jiang, F., Kula, K., Chen, J. Estimating the location of the center of resistance of canines. Angle Orthodontics. 86 (3), 365-371 (2016).
  27. Nyashin, Y., et al. Center of resistance and center of rotation of a tooth: experimental determination, computer simulation and the effect of tissue nonlinearity. Computer Methods in Biomechanics and Biomedical Engineering. 19, 229-239 (2016).
  28. Toms, S. R., Eberhardt, A. W. A nonlinear finite element analysis of the periodontal ligament under orthodontic tooth loading. American Journal of Orthodontics and Dentofacial Orthopedics. 123 (6), 657-665 (2003).
  29. Osipenko, M. A., Nyashin, M. Y., Nyashin, Y. I. Centre of resistance and centre of rotation of a tooth: the definitions, conditions of existence, properties. Russian Journal of Biomechanics. 3 (1), 5-15 (1999).
  30. Dathe, H., Nägerl, H., Dietmar, K. M. A caveat concerning center of resistance. Journal of Dental Biomechanics. 4, 1758736013499770(2013).
  31. Hohmann, A., et al. Influence of different modeling strategies for the periodontal ligament on finite element simulation results. American Journal of Orthodontics and Dentofacial Orthopedics. 139 (6), 775-783 (2011).

重印与许可

标签

有限元分析锥形束计算机断层扫描三维分割四面体网格Abaqus 软件牙周膜牙齿移动预测正畸力施加
一种用于定位上颌牙齿阻抗中心的有限元方法
联系我们
向图书馆推荐

隐私政策

使用条款

政策

版权所有 © 2026 MyJoVE Corporation。保留所有权利。