目的 通过有限元方法计算裂纹扩展过程中裂纹尖端应力强度因子的变化和裂纹扩展的路径。方法 用ANSYS软件建立釉牙本质界结构的有限元模型，在模型的单边缘设置长度为0.1 mm的初始裂纹，该数值试样的底端固定，顶端受有б=1 MPa的均匀拉伸载荷，加载频率为5 Hz。利用ANSYS软件计算裂纹在扩展过程中裂纹尖端的应力强度因子和裂纹在扩展过程中的偏斜角度。结果 裂纹到达釉牙本质界时应力强度因子会突然降低并且在釉牙本质界内随裂纹的延伸而降低。裂纹在釉质和牙本质内的偏斜角度小于15°，当裂纹扩展到釉牙本质界时会出现大的偏斜角度，裂纹尖端出现应力降低的现象。结论 由于裂纹尖端应力强度因子在釉牙本质界处的降低以及裂纹偏斜，釉牙本质界可能具有阻断裂纹扩展的作用。
釉牙本质界； 有限元； 应力强度因子； 裂纹扩展
R 322.4+1 A 10.7518/hxkq.2014.05.009
Mechanism of the dentino-enamel junction on the resist-crack propagation of human teeth by the finite element method Zheng Jingjing1， Hou Tiezhou1， Tao Hong2， Guo Xueyan1， Wu Cui1. （1. Dept. of Endodontics， Hospital of Stomatology， Xi’an Jiaotong University， Xi’an 710004， China； 2. Dept. of Stomatology， The First Affiliated Hospital of Xi’an Jiaotong Univer-sity， Xi’an 710061， China）
Objective This study aims to identify the crack tip stress intensity factor of the propagation process， crack propagation path， and the changes in the shape of the crack tip by the finite element method. Methods The finite element model of dentino-enamel junction was established with ANSYS software， and the length of the initial crack in the single edge was set to 0.1 mm. The lower end of the sample was fixed. The tensile load of 1 MPa with frequency of 5 Hz was applied to the upper end. The stress intensity factor， deflection angle， and changes in the shape of the crack tip in the crack propagation were calculated by ANSYS. Results The stress intensity factor suddenly and continuously decreased in dentino-enamel junction as the crack extended. A large skewed angle appeared， and the stress on crack tip was reduced. Conclusion The dentino-enamel junction on human teeth may resist crack propagation through stress reduction.
dentino-enamel junction； finite element； stress intensity factor； crack propagation
有限元法（finite element method，FEM）是建立在数学分析理论上的一种数值分析法。有限元法进行裂纹扩展的模拟分析一般可分为应力分析—断裂分析—建立新模型3个过程的循环。首先通过应力分析得出位移、应力、应变等力学数值。在应力分析的基础上，选择合适的断裂准则，判断裂纹是否扩展以及扩展的方向，并且可以确定每一步的扩展量，对疲劳裂纹扩展，还可以计算每一步设定扩展量的疲劳循环次数。在得出新的裂纹扩展量和扩展方向后，计算新裂纹尖端的位置即裂纹扩展的长度和偏斜角度，通过对该处的有限元网格重新划分，模拟新裂纹以及裂纹尖端，将上一步应力分析的结果映射到新的网格上，从而实现有限元模型的更新。
釉牙本质界对裂纹的阻断性能被认为是由釉质和牙本质微结构、性能逐渐变化所致，而不是两种不同材料突然转变所致。有关釉牙本质界对裂纹阻断机制的相关研究仍然存在很多观点。由于釉牙本质界位置的特殊性和牙齿发育结构的影响，制取符合条件的标准试样难度较大，目前对釉牙本质界的研究较少。多数研究仅限于探讨釉牙本质界区域的显微结构，而真正反映其力学性能的研究非常少。本研究利用有限元法分析了釉牙本质界处裂纹尖端应力强度因子以及裂纹扩展路径的变化，探讨了釉牙本质界对裂纹扩展的阻断作用。
1 材料和方法
1.1 建立分析模型
在x-y平面内，建立单向拉伸Ⅰ型裂纹的有限元模型，使裂纹平行于x轴。其中有限元模型的宽度c为1，有限元模型的长度d为10，初始裂纹a为0.1。
1.2 材料的参数
有限元模型中各种材料均视为各向同性的均质非线性弹性材料，釉质、牙本质的弹性模量分别为46.98、20 GPa；泊松比均为0.3。牙体组织脆性较大，可认为釉质、釉牙本质界、牙本质的拉伸强度约等于屈服极限，分别为11.5、33.9、46.9 MPa。
1.3 网格划分
由于裂纹尖端应力、应变存在奇异性，因此在裂纹尖端的有限元划分时应与模型的其他部分不同。因此，裂纹尖端网格的划分必须采取特殊处理，即裂纹尖端的奇异单元。其方法为：在裂纹的尖端采用圈奇异单元，其他部分用普通单元与之协调。
1.4 有限元法求解应力强度因子的理论分析
本研究选择节点法计算裂纹尖端的应力强度因子K值。使用ANSYS有限元分析软件计算应力强度因子K值时，只需将裂纹尖端进行节点划分并求解，再沿裂纹设定相应的路径，可自动求得该模型、边界条件下裂纹尖端的应力强度因子K值。
1.5 有限元分析模型的构建
为了对张开型裂纹应力强度因子进行有限元求解，本研究构建二维模型，PLANE183为平面8节点单元。KI有限元模型宽度为1 mm，其中釉质的宽度为0.45 mm，釉牙本质界的宽度为0.05 mm，牙本质的宽度为0.45 mm；KⅡ有限元模型宽度为2.05 mm，其中釉质的宽度为0.45 mm，釉牙本质界的宽度为0.05 mm，牙本质的宽度为1.35 mm。
1.6 加载、求解裂纹尖端的应力强度因子和计算裂
纹的偏斜角度
1.6.1 裂纹尖端应力强度因子的计算 实验中所选用的材料为离体牙釉牙本质界处的牙体硬组织，假定所选用的材料为各向同性材料，加载频率（R）为
5 Hz，该数值试样的底端固定，顶端受有б=1 MPa的均匀拉伸载荷。裂纹类型为张开型边裂纹。目前ANSYS软件可以根据材料的参数、加载频率和拉伸载荷自动
生成KI和KⅡ值。
1.6.2 偏斜角度的计算 裂纹偏斜角度可根据计算公式计算得出：θ=2arctg，λ=KⅡ/KI。
2 结果
2.1 裂纹扩展过程中裂纹尖端应力强度因子的变化
裂纹扩展长度a为0.1、0.2、0.25、0.3、0.4、0.45、0.45、0.5、0.55、0.65、0.75 mm，其应力强度因子分别为0.277、0.269、1.817、2.189、5.1、6.75、1.103、0.83、0.17、0.203、2.958 MPa·m1/2。有2个a=0.45，分别将其按照不同的材料进行计算。
改变模型的大小，c=2.05 mm，d=10 mm，a为0.1、0.2、0.3、0.4、0.45、0.5、0.55、0.65、0.75、0.85、0.95 mm，其应力强度因子分别为0.033、0.437、1.741、5.152、7.026、0.545、0.272、0.352、0.607、1.407、2.537 MPa·m1/2。
不同模型下裂纹尖端的应力强度因子大小的比较见图1。由图1可见，在釉质中应力强度因子随裂纹长度的增长而增大，当裂纹扩展到釉牙本质界时应力强度因子明显变小，并且裂纹在釉牙本质界内扩展时应力强度因子随裂纹的增长而降低，随着裂纹的扩展裂纹到达牙本质内，应力强度因子则又出现增大的趋势。可见，应力强度因子和裂纹长度的关系在复合材料中受到材料性质的影响。应力强度因子除了受到裂纹长度、材料性质的影响，还会受到模型大小的影响。在同等裂纹长度下，裂纹尖端的应力强度因子会随着模型的增大而降低。
图 1 不同模型下裂纹尖端的应力强度因子大小的比较
Fig 1 Comparison of the stress intensity factor of the crack tip in
different models
2.2 裂纹偏斜角度
裂纹在釉质和牙本质内扩展时偏斜的角度小于15°，故认为在釉质和牙本质内裂纹是呈直线扩展的。而裂纹在通过釉牙本质界时，会发生较大角度的偏斜。裂纹在釉质、釉牙本质界、牙本质内扩展的模拟路径见图2。
图 2 裂纹扩展路径的散点图
Fig 2 The scatter diagram of crack propagation path
3 讨论
本实验通过有限元法计算裂纹扩展过程中的不同断裂指标，即裂纹尖端的应力强度因子、裂纹尖端在路径上的应力分布、裂纹在扩展过程中的偏斜角度以及裂纹张开角度的变化。通过这些不同指标分析釉牙本质界对裂纹的阻断机制。
由于釉牙本质界宽度较窄，且结构脆性大，试样制作难度大，因此，对釉牙本质界的研究较局限。国内外研究者通过扫描电镜观察得出釉牙本质界是一个复杂的、连接硬而脆的釉质和软而韧的牙本质的至关重要的结构，对维持牙齿的生物完整性具有至关重要的作用，因为它连接着不同的钙化组织。釉牙本质界的结构一般被描述为凹面向釉质凸面向牙本质的贝壳状结构。这些结构被认为能够加强釉质和牙本质两种不同钙化组织的粘接。在本研究中通过有限元法建立的模型交界面不是凹凸面而是直线，因为在有限元模型中如果将交界面设为较小的凹凸面，则需要更小的网格划分，从而增加软件的计算功能，增加计算误差。另外，釉牙本质界可能是牙齿硬组织结构中最韧的组织，可以阻止裂纹从釉质穿过釉牙本质界到达牙本质，使裂纹局限在釉质内。采用扫描电镜对釉牙本质界的结构进行观察，Xu等采用微压痕技术估价了釉质、牙本质以及釉牙本质界附近的硬度、模量和韧性。Dong等研究结果显示，釉牙本质界的断裂韧性是（1.5±0.28） MPa·m1/2，同时通过扫描电镜观察试样断面发现，裂纹在通过釉牙本质界时出现一定程度的偏斜。由此，推测釉牙本质界可能具有阻止裂纹扩展的作用。本研究通过有限元法直接计算裂纹尖端的应力强度因子、裂纹扩展过程的偏斜角度以及裂纹尖端应力分布等指标，直接证明釉牙本质界对裂纹的阻断机制。裂纹尖端的应力强度因子在釉质内裂纹的不断延伸而增加，但是当裂纹扩展到釉牙本质界时会出现瞬间的降低并且呈不断降低趋势，裂纹到达牙本质时又会出现增加的趋势。而且釉牙本质界的断裂韧性又大于釉质的断裂韧性。裂纹要通过釉牙本质界继续向牙本质内扩展需要更大的外载荷，而口腔内的咀嚼载荷变化不大，因此，牙齿的裂纹多局限于釉质内很少通过釉牙本质界扩展到牙本质，除非在受到较大的冲击载荷即外伤或瞬间咬到硬物会造成牙齿表面裂纹的扩展或牙齿的劈裂。
另外，通过裂纹尖端在裂纹路径上的应力分布曲线可得知裂纹扩展到釉牙本质界时会降低，即釉牙本质界具有分散应力的作用。裂纹的扩展多是因为裂纹尖端是应力集中区，裂纹尖端附近的应力远大于结构内部其他部位的应力。在本研究中得出当裂纹扩展到釉牙本质界处时出现了应力降低的现象，所以，证明釉牙本质界具有分散应力的作用，从而在一定程度上阻断了裂纹的继续扩展。
本研究发现裂纹扩展到釉牙本质界时，裂纹的扩展路径发生了明显的偏斜即裂纹在釉质内近似呈直线扩展，而裂纹扩展到釉牙本质界处时裂纹出现明显的偏斜，扩展路径发生变化，在裂纹扩展的过程中可能起到一定的阻断作用。
王铎. 断裂力学. 哈尔滨：哈尔滨工业大学出版社， 1989：9-20.
黄克智， 余寿文. 弹塑性断裂力学. 北京： 清华大学出版社， 1985.
Xu HH， Smith DT， Jahanmir S， et al. Indentation damage and mechanical properties of human enamel and dentin. J Dent Res， 1998， 77（3）：472-480.
Giannini M， Soares CJ， de Carvalho RM. Ultimate tensile strength of tooth structures. Dent Mater， 2004， 20（4）：322-
329.
Tracey DM， Cook TS. Analysis of power type singularities using finite elements. Int J Numer Meth Eng， 1977， 11
（8）：1225-1233.
Cloitre T， Panayotov IV， Tassery H， et al. Multiphoton ima-ging of the dentine-enamel junction. J Biophotonics， 2013， 6（4）：330-337.
Bhaskar SN. Orban’s oral histology and embryology. 11th ed. Chicago： Mosby Year Book， 1990：70-71.
