有限元法在牙颌面功能矫形方面的应用

王 敏，王 博，王英男，李惠山，陈秉智，赵 伟，刘红彦

(1.大连医科大学 附属第一医院 口腔正畸科，辽宁 大连116011;2.大连市友谊医院 口腔科，辽宁 大连116001;3.大连交通大学轨道交通学院，辽宁大连116028;4.大连理工大学工业装备结构分析国家重点实验室，辽宁大连116024)

有限元法是将连续的弹性体分割成有限个单元，以其结合体来代替原弹性体，并逐个研究每个单元的性质，以获得整个弹性体的力学分析方法［1］。是当今工程分析中广泛应用于分析结构应力应变的一种数值算法［2］。1943年Courant R［3］首先提出有限元法基本思想，1956年Turner MS 等［4］在航空工业的飞机结构分析中首次成功应用有限元法。1960年Clough［5］明确提出Finite Element Method(FEM)这一概念，使人们更清楚地认识到有限元法的特性和功用。随着计算机计算速度的提高和运算软件的发展，有限元法完成了从二维到三维的飞跃，并能应用到各种复杂问题的研究中。自1973年Thresher［6］首先将有限元法应用于口腔医学，有限元法已经成为口腔生物力学研究领域中一种有效的分析工具，为口腔疾病治疗、医疗器械的优化设计等提供理论依据。

1 实验方法

传统口腔医学中有限元建模多采用磨片法、切片法、描片法以及标本的人工测读等方法［7］。这些方法的缺点在于出错较高，无法准确表达细部结构特征。近年来，有限元程序及软件的强大建模功能，可以逼真地建立三维人体骨骼、肌肉等器官组织的模型，并赋予其准确的生物力学特性。同时，基于CT 输出的DICOM 结果研发的接口工具能够快速准确的将CT 结果转化为有限元模型。近年来，MRI由于可以观测到软骨及软组织，也成为了有限元模型建立的一个来源［8］。

模型的约束方式:根据不同的试验目的，常对模型的不同部位进行约束，以更有效地进行试验。例如，在研究正畸矫治弓丝入槽后，牙齿的受力情况时，根据对称性原则，在正中联合施加对称约束;在髁突处将6 个自由度全部约束;在托槽和牙齿之间的节点处以梁单元(Beam)连接固定，即可以模拟真实情况下，口腔内托槽在牙面上的粘接固位［9-10］。在此不同部位都相应地经行约束后，即可通过软件计算出试验结果。

2 下颌功能矫形

针对下颌的功能矫形矫治器的原理就是髁突软骨在生长期或生长完成后对功能刺激有适应性改建的能力。因此，明确各种矫治器产生的生物机械力在髁突软骨处的应力分布，能够更好的指导临床上矫治器的应用。

2.1 下颌后缩的功能矫形

周学军等［11］在建立的下颌骨三维有限元模型上分别模拟不同程度的下颌前伸状态，结论提示临床咬合重建下颌前伸距离不应过大，对临床有很强的指导意义。但应力分布的规律不变:以张应力为主，主要集中于后部，应力值由前向后逐渐增大，由内向外也逐渐增大，最大应力集中于后外侧部。

根据Anurag Gupta 等［12-13］的实验，下颌前伸后可在髁突的后上方产生拉应力，前方和前上方产生压应力，关节窝的后部也有拉应力，这可能与这些区域的生长和细胞增殖有关。他们更进一步建立多个咬合打开高度的模型，发现随着咬合重建高度的增加，对于关节的应力分布更有利，提高了下颌骨对于功能矫治器的髁突反应。所以，在下颌后缩的功能矫形中，各位学者所得的关于应力分布的结论基本一致，都是在髁突的后上方产生拉应力，前上方产生压应力，下颌角部由于有肌肉附着也有应力分布，同时在颈部产生最大的应力集中。这些应力集中区域可产生功能改型，从而发挥矫治效应。在力学机制应用上，相关的有限元实验也提供了很多临床使用的建议和参考，但是下颌功能矫治器种类繁多，并没有各种矫治器对比的文献出现，矫治器的效应的差别没有文献能够体现。

2.1.1 TWIN -BOCK 矫治器:早期的有限元法的应用主要在二维的层面上开展，GD Singh 和WJ Clark［14-15］使用有限元尺度分析法分析了用TWINBOCK 矫治器治疗后的患者的头颅侧位片，研究了下颌的位置改变和软组织的改变，发现使用TWINBOCK 矫治器达到矫治下颌形态的目的主要反映在髁颈部的定位生长伴随喙突结构改变中，可能涉及发展调制髁突软骨，重构分支，并有新骨沉积在牙槽骨。软组织方面，唇肌的力量可帮助软组织的改进。

2.1.2 Frankel 矫治器:白丁等［16］应用二维有限元分析法研究了12 例采用Frankel 矫治器典型的安氏Ⅱ类1 分类错颌患者发现功能矫治器使前面部增高的效应强于对后面高的影响，使下颌角增大，同时能增加下颌骨量，促进下颌骨水平生长，从而矫正Ⅱ类骨骼关系。对比头影测量技术，二维有限元法和有限元尺度分析(FESA)可以在不需要任何边界条件和参照系的条件下，分析任何形状的有限单元的变形［17］，从而反映颅面结构内部的各向异性变形。同时，在设计实验时，可将样本随机对照分组，减少误差［18］。

2.1.3 activator 矫治器:Cagri Ulusoya 等［19］利用体外的人类下颌骨建立了三维有限元模型，并在模型上固定了Class II activator 矫治器的下半部分，对比研究了Class II activator 和Class II activator 高位牵引的应力分布，发现肌肉附着处的应力最大，表现为髁突顶端和下颌角处应力值最大，同时通过实验证明两种装置都能通过激活咀嚼肌改变下颌骨生长方向继而改变下颌骨形态。

2.1.4 Herbst 矫治器:胡林华等［20-24］使用活体CT建立模型，用柔索约束、受压间隙元等形式模拟咀嚼肌、韧带的边界约束，极大提高了模型的相似性，并通过一系列的实验研究Herbst 矫治器在不同的合重建状态下髁突软骨表面应力分布的影响以及相关口颌肌肉和韧带的约束反力变化规律。确定了合重建的生理承受范围，认识到了茎突下颌韧带、蝶下颌韧带、颞肌后束和颞肌前束在Herbst 矫治器所引起的下颌骨功能改建中的重要作用。

2.1.5 Forsus 矫治器:颜功兴等［25］研究了在Forsus矫治力作用下下颌骨的应力分布情况和位移分布的情况，认为最大应力出现在弓丝托槽，下颌骨上的应力集中主要位于髁突颈和牙槽嵴部位，下颌骨上最大的Von Mises 等效应力出现在髁突颈处，下颌骨最大的位移都基本出现在颏部，髁突颈局部发生了塑性变形。结论Forsus 矫治力的加载角度最好为0° ～25°。

2.2 下颌过度生长的功能矫形

Faruk Ayhan Basciftci 等［26］通过在包含颞下颌关节的下颌骨三维有限元模型上加载各个方向的颏兜作用力，得出的结论有:当力的方向通过髁突时下颌骨可发生向下向后的位置改变，当力的方向通过或位于冠状突上时，下颌骨会发生向下向后的位置改变，不管力的方向如何，髁突和冠状突都有最小的位置改变，最大的应力分布区域在髁突和下颌升支后缘。杨辉等［27］用三维有限元法观察颏兜力作用下颞下颌关节的受力情况和整个下颌骨的应力分布，证实了颏兜力不会造成关节损伤，周学军等［28］通过对大小不同的通过髁突中心的作用力的有限元分析得出的结论是通过髁突中心的作用力无论大小均难以抑制下颌髁突的生长，仅能改变生长方向，下颌骨发生一定的生长变形。这些三维有限元实验的结果和临床多年观察的结果一致。

3 上颌功能矫形

在上颌骨有限元模型建立中，边界约束条件对最终的计算结果有很大的影响，赵志河等［29］采用枕骨大孔处全方位约束，Cattaneo PM 等［30］采用咬肌附丽处与咬合面加载，而约束位于模型的断面;Gross MD 等［31］在研究中采用咬合面加载，咀嚼肌附着处约束的方法。白石柱等［32］对上颌骨的有限元研究中在咬合面加载，在咀嚼肌附着处采用约束的方式来模拟咀嚼肌的影响。骨缝在矫形力的传递和分布中起了重要作用，骨缝形态的复杂性决定了骨缝对矫形力的反应也是复杂的，这除了与矫形力的部位、方向和大小有关外，还与骨的几何形状、骨缝的走行方向、生物力学特性有关［33］。在评价骨缝生物力学方面，FEM 有其独特优势。Tanne 等［34-36］在其上颌前牵引和后牵引的有限元法研究中均发现骨缝线处应力轮廓出现中断现象，骨缝两侧骨界面的应力类型很不一致。张强等［37］研究证实后牵引矫形力的方向对上颌腭部位移影响明显;而对应力分布影响较小。

Wang D 等［38］建立单侧完全性唇腭裂患者上颌复合体三维有限元模型，利用三维有限元应力分析法研究扩弓力在UCLP 患者上颌复合体的传递模式、骨缝应力状态及颅颌面各骨块位移特性，表明蝶骨的翼板近颅底处是约束上颌扩弓最主要的因素，上颌前缘也受较大应力，如果医师能松解翼板及打断上颌前缘，临床将获得更大的扩弓量，上颌骨、鼻腔及颧骨扩弓后也将趋于稳定，其在保持期复发回缩的可能性将减小。Provatidis C 等［39］建立了相对精确的上颌骨模型，包括完整的牙列和牙周韧带，来评估最真实的矫治力的效应，结果显示相关的应力和应变集中于颅骨，上颌骨以及牙周韧带处。

雷勇华等［40］模拟前牵引加腭部扩弓治疗，得出上颌前牵引联合应用腭部扩弓治疗唇腭裂患者，有利于全面改善上颌骨的发育的结论。

由于上颌骨形态比较复杂，骨缝关系难以简单分类，腭裂形态种类也比较多，可以通过在模型上改变一些细节，比如连接关系，单元的形式和数值，可以做出个体化的模型，然后达到对不同类型病人的临床指导［41-42］。建立多个模型来达到模拟临床上不同患者的形态，来尽量真实的反应客观情况，这也将是未来研究发展的趋势。

综上所述，有限元法作为一种经典的实验力学的研究方法，已经在牙颌面功能矫形中取得了很多成果，随着计算机技术以及各种模拟数字人体技术的发展，预计在不久的将来，就可以实现口腔环境的全真模拟;并且，在准确使用有限元技术的基础上，结合生物力学知识，有可能通过计算机准确的估计出各种口腔治疗时，牙齿以及骨组织的受力情况及变化趋势。

［1］ 龙驭球. 有限元法概论［M］. 北京:高等教育出版社，1978:1 -128.

［2］ 罗阳，府伟灵，王旭全，等.有限元分析全环与半环型外固定架的力学稳定性差异［J］.重庆医学，2008，37(7):740.

［3］ Courant R. Variational Method for Solutions of Problems of Equilibrium and Vibrations［J］. Bull Am Math Soc，1943，49:1 -3.

［4］ Turner MS，Clough RW，Marcin HC，et al. Stiffnessand deflection analysis complex structure［J］. J Aero Sci，1956，23(9):805 -809.

［5］ Clough RW. The Finite Element Method in Plane Stress Analysis［J］. Proceeding of the 2nd ASCE Conference on Electronic Computation Pittsburgh，PA，1960，9:345.

［6］ Thresher RW. The Stress Analysis of Human Teeth［J］. J Biomech，1973，6:443.

［7］ Nakano H，Watahiki J，Kubota M，et al. Orthod Cranofac Res，2003，6(Suppl 1):168 -172.

［8］ Nishio C，Tanimoto K，Hirose M，et al. Stress analysis in the mandibular condyle during prolonged clenching:a theoretical approach with the finite element method［J］. Proc Inst Mech Eng H，2009，223(6):739 -748.

［9］ Ludwig B，Baumgaertel S，Zorkun B，et al. Application of a new viscoelastic finite element method model and analysis of miniscrew-supported hybrid hyrax treatment［J］. Am J Orthod Dentofacial Orthop，2013，143(3):426 -435.

［10］ Canales C，Larson M，Grauer D，et al. A novel biomechanical model assessing continuous orthodontic archwire activation［J］. Am J Orthod Dentofacial Orthop，2013，143(2):281 -290.

［11］ 周学军，赵志河，赵美英，等. 不同程度下颌前伸的三维有限元分析［J］. 安徽医科大学学报，2004，39(6):419 -422.

［12］ Anurag Gupta，Virender S Kohli，Pushpa V. Hazarey，Om P. Kharbanda，Amit Gunjal. Stress distribution in the temporomandibular joint after mandibular protraction:A3 -dimensional finite element method study. Part 1［J］.Am J Orthod Dentofacial Orthop，2009，135:737 -748.

［13］ Anurag Gupta，Virender S Kohli，Pushpa V. Hazarey，Om P. Kharbanda，Amit Gunjal. Stress distribution in the temporomandibular joint after mandibular protraction:A3 -dimensional finite element method study. Part 2［J］.Am J Orthod Dentofacial Orthop，2009，135:749 -756.

［14］ GD Singh，WJ Clark. Localization of mandibular changes in patients with Class II Division 1 malocclusions treated with Twin-block appliances:Finite element scaling analysis［J］. Am J Orthod Dentofacial Orthop，2001;119:419 -425.

［15］ GD Singh，WJ Clark. Soft tissue changes in patients with Class II division 1 malocclusions treated using Twin Block appliances:finite - element scaling analysis［J］. Eur J Orthod，2003，25:225 -230.

［16］ 白丁，赵美英. Frankel 矫治器治疗安氏Ⅱ类1 分类错合的有限元分析［J］. 华西口腔医学杂志，1995，14(2):97 -100.

［17］ Singh D，Medina LE，Hang WM. Soft tissue facial changes using Biobloc appliances:geometric morphometrics［J］. Int J Orthod Milwaukee，2009，20(2):29 -34.

［18］ NganP，scheick J，Florman M.A tensor analysis toevaluate the effect of high—pull headgear on Class Ⅱmlaocclusion［J］.Am J Orthod，1993，103:267.

［19］ Cagri Ulusoya，Nilüfer Darendeliler. Effects of Class II activator and Class II activator high -pull headgear combination on the mandible:A 3 - dimensional finite element stress analysis study［J］.Am J Orthod Dentofacial Orthop 2008，133:490.e9 -490.e15.

［20］ 胡林华，赵志河，宋锦磷，等. 颞下颌关节-下颌骨-Herbst 矫治器系统三维有限元模型的建立［J］.广东牙病防治，2004，12(3):224 -225.

［21］ 周学军，赵志河，赵美英，等. 包括下颌骨的颞下颌关节三维有限元模型的建立［J］. 实用口腔医学杂志，2000，16(1):17 -19.

［22］ 宋锦磷，赵志河，胡林华，等. Herbst 矫治器在不同合重建时对口颌肌肉和韧带约束反力的影响［J］.华西口腔医学杂志，2001，19(1):43 -45.

［23］ 周学军，赵志河，赵美英，等.下颌骨三维有限元模型的边界约束设计［J］.华西口腔医学杂志，1999，19(1):29-32.

［24］ 胡林华，赵志河，宋锦磷，等. Herbst 矫治器在不同殆重建时对髁突软骨表面应力分布的影响［J］. 2001，19(1):46 -48.

［25］ 颜功兴，刘占芳. Forsus 矫治力前导下颌有限元仿真及最佳施力角度探讨［J］. 第三军医大学学报，2009，31(19):1890 -1893.

［26］ Faruk Ayhan Basciftci，Hasan H sn Korkmaz，Oguz Eraslan. Biomechanical evaluation of chincup treatment with various force vectors［J］. Am J Orthod Dentofacial Orthop，2008，134:773 -781.

［27］ 杨辉，刘洪臣. 颏兜力作用下颞下颌关节及下颌骨受力的三维有限元分析［J］.口腔颌面修复学杂志，2000，1(3):139 -142.

［28］ 周学军，赵志河，赵美英，等.不同大小颏兜牵引力的下颌骨三维有限元分析［J］. 北京口腔医学，2004，12(3):125 -129.

［29］ 赵志河，房兵，赵美英.颅面骨三维有限元模型的建立［J］.华西口腔医学杂志，1994，1(24):298 -300.

［30］ Cattaneo PM，Dalstra M，Frich LH. A three -dimensional finite element model from computed tomography data:a semi-automated method［J］.Proc Inst Mech Eng，2001，21(52):203 -213.

［31］ Gross MD，Arbel G，Hershkovitz I.Three-dimensional finite element analysis of the facial skeleton on simulated occlusal loading［J］. J Oral Rehabi，2001，2(87):684 -946.

［32］ 白石柱，李涤尘，赵铱民，等. 上颌骨有限元分析中边界约束条件的研究［J］. 临床口腔医学杂志，2006，22(12):720 -723.

［33］ Mao JJ，Wang X，Kopher RA.Biomechanics of craniofacial sutures:orthopedic implications［J］. Angle Orthod，2003，73(2):128.

［34］ Tanne K，Miyasaka J，Yamagata Y，et al. Three—dimensional model of the human craniofacial skeleton:method and preliminary results using finite element analysis［J］.J Biomed Eng，1988，10(3):240.

［35］ Tanne K，Hiraga J，Kakiuchi K，et al.Biomechanical effect of anteriorly directed extraoral forces on the craniofacial complex:A study using the finite element method［J］.Am J Orthod Dentofac Orthop，1989，95(3):200.

［36］ Tanne K，Hiraga J，Sakuda M.Effect of directions of maxillary protraction forces on biomechanical changes in craniofacial complex［J］.Eur J 0rthod，1989，11(4):328.

［37］ 张强，赵志河，王军，等. 后牵引矫形力的方向对上颌腭部位移及应力分布的影响［J］. 四川大学学报(医学版)，2004，35(5):680 -682.

［38］ Wang D，Cheng L，Wang C，et al. Biomechanical analysis of rapid maxillary expansion in the UCLP patient［J］.Medical Engineering ＆ Physics，2009，31:409 -417.

［39］ Provatidis C，Georgiopoulos B，Kotinas A，et al. On the FEM modeling of craniofacial changes during rapid maxillary expansion［J］.2007，29:566 -579.

［40］ 雷勇华，翦新春，任毕乔. 前方牵引下颅面复合体骨缝应力特征的三维有限元分析［J］. 中南大学学报(医学版)，2009，3:19 -23.

［41］ Lee H，Ting K，Nelson M，et al. Maxillary expansion in customized finite element method models［J］. Am J Orthod Dentofacial Orthop，2009，136(3):367 -374.

［42］ Pan X，Qian Y，Yu J，et al. Biomechanical effects of rapid palatal expansion on the craniofacial skeleton with cleft palate:a three -dimensional finite element analysis［J］. Cleft Palate Craniofac J，2007，44(2):149 -154.