基于奇异元计算分析裂纹尖端应力强度因子

1、51卷 第3期（总第192期） 中 国 造 船 V ol.51 No.3 (Serial No. 192 2010年9月 SHIPBUILDING OF CHINA Sep. 2010文章编号：1000-4882 (2010 03-0056-09基于奇异元计算分析裂纹尖端应力强度因子陈景杰，黄 一，刘 刚(大连理工大学 船舶工程学院, 辽宁 大连 116024摘 要采用两种奇异单元模拟裂纹尖端应力应变场的奇异性，建立了相应的计算裂纹尖端应力强度因子的ANSYS 有限元模型。通过数值计算，分别考察了这两种有限元模型中裂纹尖端附近区域网格参数的变化对应力强度因子计算精度的影响，比较了应力强度因子对

2、各个参数的敏感程度。发现采用20节点奇异元的有限元模型计算的应力强度因子几乎与网格参数无关，其计算结果更稳定可靠。该模型能够用于船舶及海洋工程结构中含裂纹构件的应力强度因子计算。关 键 词：船舶、船舶工程; 应力强度因子; 奇异单元; 有限元模型; 网格参数 中图分类号：U661.4 文献标识码：A0 引 言在复杂的海洋环境载荷作用下，船舶及海洋工程结构物在其服役期内会发生结构疲劳损伤1，而关键承载构件上疲劳裂纹的出现是结构疲劳损伤的必然结果。对于那些含有疲劳裂纹的结构，应力强度因子是评估它们断裂失效的重要参量，在裂纹体分析中占据着重要地位。自断裂力学问世以来，已经产生了众多理论和数值解法求解

3、应力强度因子，较为典型的计算方法有有限元法、边界元法、解析法、混合法、权函数法和线弹簧法等。利用这些方法，已经获得许多裂纹体模型的应力强度因子的大小，应力强度因子手册2中收编了许多种典型裂纹体模型应力强度因子的解，但由于裂纹体几何形式及所受载荷的复杂性，很多情况下应力强度因子的解难以从现有手册查到。而有限元方法不受裂纹体几何及载荷形式的限制，因而在断裂力学中得到了非常广泛的应用。位移法是用有限元法计算裂纹尖端应力强度因子的常用方法，它是利用奇异单元来精确描述线弹性范围内裂纹尖端的应力和应变场的奇异性。取裂纹尖端第一层单元为奇异元，其余各层单元为相应等参元，用这样的方法建立裂纹体的有限元模型31

4、1来计算裂纹尖端的应力强度因子。目前有两种奇异单元被广泛使用，分别为20节点奇异元和12节点奇异元，它们是由20节点等参元蜕化、转变获得，形式如图1所示。文中以ANSYS 软件为基础，分别采用20节点和12节点奇异单元建立了含表面裂纹的平板的两种有限元模型。对有限元网格技术及裂纹体网格剖分的要求作了深入研究，评估了两种有限元模型计算的应力强度因子对网格的敏感程度，进而判断出哪种有限元模型计算的应力强度因子更稳定可靠，这为有限元法确定船舶及海洋工程中的局部复杂结构裂纹尖端应力强度因子提供了依据。收稿日期：2009-12-19；修改稿收稿日期：2010-04-0851卷 第3期 (总第192期 陈

5、景杰，等：基于奇异元计算分析裂纹尖端应力强度因子 57（a ）20节点的奇异元 （b ）12节点的奇异元 （c ）20节点奇异元围绕裂尖形式图1 奇异单元类型及其围绕裂尖形式1 有限元模型及网格参数对图2所示的含半椭圆表面裂纹的平板模型进行有限元计算和分析，在有限元模型中考虑了裂纹前缘周围区域网格的变化对应力强度因子的影响，分析不同奇异单元所建的有限元模型计算结果的稳定性。含有表面裂纹的有限平板示于图2，图中的各个参数列于表1。在平板长度两端受到垂直裂纹面方向的均匀拉伸载荷作用。根据模型及载荷的对称性采用奇异单元法建立1/4有限元模型，如图3所示。图4给出了图3中裂纹尖端附近区域网格剖分情况、

6、裂纹尖端端面形状及影响应力强度因子精度的主要参数，参数含义见表1。 图2 含有表面裂纹的有限平板几何模型 图3 有限平板的有限元模型 （a ）沿半裂纹尖端曲线 （b ）裂纹深度方向横截面 区域附近的有限元模型n um_crack=16 的有限元形状，=45°,n um_tip =5（c ）裂纹尖端半端面的局部放大图 图4 裂纹尖端附近区域的有限元模型 1/4节点 中间节点角节点 裂纹尖端58 中 国 造 船 学术论文 表1 含有表面裂纹的平板的主要参数符号 参数 2L 平板长度 2W 平板宽度 T 平板厚度 2a 表面裂纹长度 b 表面裂纹最大深度 L tip裂纹尖端附近区域大小 裂

7、纹尖端单元角度 n um_tip 围绕裂纹尖端单元层数 n um_crack 裂纹尖端曲线长度所分的份数 L in 裂纹尖端奇异单元长度 L在L tip 内最外层单元长度2 网格参数变化对计算结果影响对于采用奇异单元建立的有限元模型，分析裂纹尖端附近网格变化对计算结果的影响。在这方面，目前主要是研究了围绕裂纹前缘的单元层数n um_tip及奇异单元尺寸L in 与裂纹深度b 比这两个因素的影响。关于n um_tip，文献3-8都采用了3层以上的网格，文献4更是明确建议裂纹前缘网格层数不能少于3层。关于L in 的选择问题，文献4建议L in / b应该小于0.1，文献3建议该比值应该介于0.0

8、5和0.15之间，而有限元分析软件ANSYS 的技术文件则建议L in /b 小于1/8。本文不仅研究了上述两个因素对应力强度因子的影响，而且进一步研究裂纹尖端单元角度、裂纹尖端附近区域大小L tip 和裂纹尖端曲线的份数n um_crack对该因子的影响。 2.1的变化对应力强度因子的影响在图2所示的几何模型中，取2a =40mm，b =10mm，2W =120mm，2L =200mm，T =20mm，均匀拉应力为100MPa ，建立两种奇异元的有限元模型。在n um_tip=5、L out /L in =1、 L tip /d =1/10和n um_crack=16的两种有限元网格中，研究

9、裂纹尖端附近区域分别为9°、15°、22.5°、45°时单元角度的变化对有限元模型计算的应力强度因子的影响。图5给出了不同对应的裂纹尖端附近区域的横截面。根据上述的模型尺寸及网格参数，两种有限元模型计算的应力强度因子如图6和图7所示，它们给出了20节点和12节点奇异元计算的应力强度因子对的敏感程度。在图6和图7中SIF 表示裂纹最深位置尖端的应力强度因子，在后面各图中的SIF 具有相同含义。 根据图6中裂纹尖端应力强度因子随单元角度的变化情况，不难得出在20节点奇异元的有限元模型中，裂纹尖端角度的变化对应力强度因子几乎没有影响，为了提高计算效率，在后面对

10、于20节点奇异元的有限模型只采用=45°来计算裂纹尖端的应力强度因子。 （a ）=45° （b ）=22.5° （c ）=15° （d ）=9°图5 取不同值时裂纹尖端附近区域的横截面51卷 第3期 (总第192期 陈景杰，等：基于奇异元计算分析裂纹尖端应力强度因子 59 图6 不同值时20节点奇异元计算的 图7 不同值时12节点奇异元计算的应力强度因子 应力强度因子图7给出了12节点奇异元的有限元模型计算的应力强度因子随着的增加而不断减小的变化趋势，说明裂纹尖端角度的变化对该模型计算结果的影响不可忽略。结合计算精度和计算效率，在后面对于12节

11、点奇异元的有限元模型应取裂纹尖端单元角度=22.5°来计算裂纹尖端的应力强度因子。 2.2 n um_tip的变化对应力强度因子的影响几何模型与2.1节相同，在20节点和12节点奇异元的有限元网格中，取L out /L in =1、 L tip /d =1/10和n um_crack=16，而分别为45°和22.5°来分析围绕裂纹尖端单元层数的变化对有限元模型计算的应力强度因子的影响。文中考虑了n um_tip=3、5、8、20的情况。20节点奇异元模型对应的裂纹尖端附近区域的横截面形状如图8所示，不同层数下两种有限元模型计算的应力强度因子如图910所示。 （a

12、）n um_tip=3 （b ）n um_tip=5 （c ）n um_tip=8 （d ）n um_tip=20图8 n um_tip取不同值时裂纹尖端附近区域的横截面形状 图9 不同n um_tip值20节点奇异元计算的 图10 不同n um_tip值12节点奇异元计算的应力强度因子 应力强度因子图9给出了不同n um_tip下20节点奇异元有限元模型计算的应力强度因子的变化情况，说明围绕裂纹尖端单元层数的变化对应力强度因子的影响可以忽略，在计算裂纹尖端的应力强度因子有限元模型中n um_tip可取320之间的任意值。60 中 国 造 船 学术论文 对于12节点奇异元有限元模型计算的应力强

13、度因子（SIF ）结果见表2。可以看到当n um_tip=20时，应力强度因子值有较大的变化。为了进一步分析围绕裂纹尖端的单元层数对有限元模型计算的应力强度因子的影响，又取n um_tip=11、15、18、19的情况进行计算，计算结果见图10。表2 不同n um_tip值下12节点奇异元模型计算的应力强度因子 单位：1/2MPa mm n um_tipum_tipSIF20 576.6983从图10可以看出，当n um_tip> 18时，围绕裂纹尖端单元层数的变化对计算的应力强度因子有一定的影响，而n um_tip在318的范围内，采用12节点奇异元的有限元模型计算的应力强度因子几乎相

14、等，说明采用12节点奇异元的有限元模型计算裂纹尖端应力强度因子时，n um_tip可取318的范围内的任意值。2.3 L out /L in 的变化对应力强度因子的影响几何模型不变，基于上述研究结果，在20节点和12节点奇异元的有限元模型中取n um_tip=5，分别为45°和22.5°，同时假定L tip /d =1/10和n um_crack=16，研究不同L out /L in 下两种有限元模型计算应力强度因子的变化情况，结果如图12和图13所示。图11给出了L out /L in 分别为1、3、5和10时有限元模型中裂纹尖端附近区域的横截面形状。 （a ）L out

15、 /L in =1 （b ）L out /L in =3 （c ）L out /L in =5 （d ）L out /L in =10图11 Lout /L in 取不同值时裂纹尖端附近区域的横截面形状 图12 不同L out /Lin 值时20节点奇异元计算的 图13 不同L out /Lin 值时12节点奇异元计算的应力强度因子 应力强度因子图12和图13分别描述了20节点和12节点奇异元对应的有限元模型计算的应力强度因子随L out /Lin 的变化情况。根据图12中的计算结果，最大值与最小值的相对误差为0. 35%，说明不同L out /L in 的变化对于20节点有限元模型计算的应力

16、强度因子几乎没有影响，采用该模型计算应力强度因子时，L out /L in 可取130内的任意值。后面取L out /L in =1对应的有限元模型计算裂纹尖端的应力强度因子。在图13中，采用12节点奇异元的有限元模型计算的应力强度因子具有随着L out /L in 的增大而增大51 卷 第 3 期 (总第 192 期 陈景杰，等：基于奇异元计算分析裂纹尖端应力强度因子 61 的趋势，说明这种有限元模型计算的应力强度因子对 Lout /Lin 参数具有较大的灵敏性，不同的网格剖分 形式会产生不同的计算的结果，导致计算结果不可靠。与 Newman-Raju 法1计算的应力强度因子进行 比较，当

17、Lout /Lin=4 时，两种方法计算结果的相对误差不足 0.1%，说明对于 12 节点奇异元的有限元模 型，通过改变 Lout /Lin 值的大小能够提高模型的计算精度。但是经研究发现 Lout /Lin 不仅与裂纹形状比 b/a 有关，还与外载荷加载形式有关，从而很难确定适当的 Lout /Lin 来建立相应有限元模型。 Ltip/b 的变化对应力强度因子的影响 几何模型不变，根据前面研究结果，在 20 节点奇异元的有限元模型中，取 =45°，num_tip=5，Lout /Lin=1； 12 节点奇异元的有限元模型中， =22.5°， um_tip=5， out /

18、Lin=4， 在 取 n L 同时两种模型中令 num_crack =16，研究 Ltip/b 在 1/801/4 范围内的变化时对有限元模型计算的应力强度因子的影响，结果如图 14 和图 15 所示。 590 600 595 SIF / MPa·mm 1/ 2 2.4 SIF / MPa·mm 1/ 2 585 590 585 580 575 580 575 570 0.00 0.05 0.10 0.15 0.20 0.25 0.30 L tip/b 图 14 不同 Ltip/b 值时 20 节点奇异元计算的 应力强度因子 0.00 0.05 0.10 0.15 L ti

19、p/b 0.20 0.25 0.30 图 15 不同 Ltip/b 值 12 时节点奇异元计算的 应力强度因子 从图 14 可以看出，随着 Ltip/b 的增加，20 节点奇异元有限元模型计算的应力强度因子几乎不变， 说明 Ltip 的变化对有限元模型计算裂纹尖端应力强度因子几乎没有影响，Ltip/b 可以取 1/401/4 范围内 的任意值。 图 15 给出了在不同 Ltip/b 下，12 节点奇异元的有限元模型计算的应力强度因子的变化曲线。当 Ltip/b<0.05 时，计算的应力强度因子随之 Ltip/b 的减小迅速增加，为了使计算结果稳定，Ltip/b 应取 0.050.25 范

20、围内的值进行裂纹尖端应力强度因子的计算。 2.5 num_crack 的变化对应力强度因子的影响 几何模型不变，在 20 节点奇异元的有限元模型中裂纹尖端附近区域的网格参数分别为 =45°， L 在 n num_tip=5， out /Lin=1 和 Ltip/b=0.1； 12 节点奇异元的有限元模型中网格参数分别为 =22.5°， um_tip=5， Lout /Lin=4 和 Ltip/b=0.1， 研究裂纹尖端曲线的份数变化对应力强度因子的影响。 在此取 num_crack 为 8， 12， 16 和 20 进行分析，计算结果如图 1617 所示。 图 1617 给

21、出了不同 num_crack 值对 20 节点和 12 节点奇异元模型计算应力强度因子的影响。图 16 说明在 20 节点奇异元模型中，num_crack 值的变化对应力强度因子没有影响，在模型剖分中可取任意值。 而图 17 给出在 12 节点奇异元有限元模型中，num_crack 的变化对应力强度因子具有一定影响，当 num_crack>12 时，影响趋势减缓，结合计算效率，对于 12 节点奇异元的有限元模型应在 1220 的范围内 对 num_crack 取值来计算应力强度因子。 62 中 国 造 船 学术论文 590 585 580 575 570 5 10 15 n um_cra

22、ck 图 16 不同 num_crack 值时 20 节点奇异元计算的 应力强度因子 590 585 580 575 570 20 25 SIF / MPa·mm1/ 2 SIF / MPa·mm1/ 2 5 10 15 n um_crack 20 25 图 17 不同 num_crack 值时 12 节点奇异元计算的 应力强度因子 综上所述，12 节点奇异元具有较少的节点，并且仅需与 8 节点单元相连接，大大减少了有限元模 型中的节点个数，能提高计算效率，但经过上述分析可知，这种奇异单元的有限元模型计算应力强度 因子的稳定性与裂纹尖端附近区域网格形式有关，导致计算结果不可

23、靠。而 20 节奇异元虽然节点数较 多，但该模型计算的应力强度因子十分稳定，它不受有限元模型中网格剖分形式的影响。对于图 2 所 示的模型，采用 20 节点奇异元有限元模型计算的应力强度因子为 580.561MPa mm1/2 与 Newman-Raju 二者相对误差为 2.68%。 这说明计算裂纹尖端应力强度 法12计算的结果 596.548 MPa mm1/2 进行比较， 因子的 20 节点奇异元有限元模型不仅具有很好的稳定性而且具有较高的计算精度， 它的计算结果是正 确可靠的。 3 结 论 本文采用有限元数值计算，研究了 20 节点和 12 节点奇异元这两种有限元模型中网格参数的变化 对

24、应力强度因子的影响。通过分析比较可得到如下结论：在计算裂纹尖端应力强度因子的有限元模型 中，20 节点奇异元的有限元模型是稳定可靠的，它所计算的应力强度因子与裂纹尖端附近区域的网格 状态无关，在数值模拟分析中能够得到稳定的计算结果；12 节点奇异元的有限元模型计算结果与裂纹 尖端附近区域的网格状态密切相关，对于相同问题，在数值模拟中不同剖分形式得到的计算结果差别 较大，显示该模型在应用中不够稳定，难以处理复杂结构的裂纹问题。因此对于船舶及海洋工程结构 物等复杂结构的裂纹问题宜采用 20 节点奇异元的有限元模型来确定裂纹尖端的应力强度因子， 以便于 评估结构所处的安全状态。 参 考 文 献： 1

25、 王一飞, 衣高洁, 黄小平, 崔维成. 疲劳裂纹扩展模型中表征裂纹闭合水平参数的确定J.中国造船, 2008, 49(4: 45-51. 2 3 中国航空研究院.应力强度因子手册M. 北京: 科学出版社, 1993. INGRAFFEA A R, MANU C. Stress intensity factor computation in three dimensions with quarter-point elements J. International Journal for Numerical Methods in Engineering, 1980, 15: 1427-1445.

26、 4 CAO J J, YANG G J, PAKER J A, BURDEKIN F M. Crack modeling in FE analysis of circular tubular joints J. 51 卷 第 3 期 (总第 192 期 陈景杰，等：基于奇异元计算分析裂纹尖端应力强度因子 63 Engineering Fracture Mechanics, 1998, 61:537-553. 5 CHIEW S P, LIE S T, LEE C K, HUANG Z W. Stress intensity factors for a surface crack in a t

27、ubular T-joint J. International Journal of Pressure Vessels and Piping, 2001, 78:677-685. 6 BOWNESS D, LEE M M K. Prediction of weld toe magnification factors for semi-elliptical cracks in T-butt joints J. International Journal of Fatigue, 2000, 22:369-387. 7 LEE M M K, BOWNESS D. Fatigue analysis o

28、f tubular joints using finite elements J. Fatigue in Offshore Structures, 1996, 1: 223-257. 8 SEETHARAMAN S, SREEDHAR D S. Fatigue analysis and fracture mechanics J. Fatigue in Offshore Structures, 1996, 1:259-309. 9 WANG X, LAMBENT S B. On the calculation of stress intensity factors for surface cra

29、cks in welded pipe-plate and tubular joints J. International Journal of Fatigue, 2003, 25:89-96. 10 卢炎麟, 王自力, 杜水友. T 型管节点焊趾表面裂纹扩展行为的三维断裂力学分析J. 中国造船, 1995, 2:51-59. 11 BOWNESS D, LEE M M K. The development of an accurate model for the fatigue assessment of doubly curved cracks in tubular joints J. In

30、ternational Journal of Fatigue, 1995, 73:129-147. 12 NEWMAN J C. RAJU I S. An Empirical Stress-Intensity Factor Equation for the Surface Crack J. Engineering Fracture Mechanics, 1981, 15(1-2:185-192. Analysis of Finite Element Model for Calculating Stress Intensity Factor Based on Crack-tip Singular

31、 Element CHEN Jing-jie HUANG Yi LIU Gang (School of Naval Architecture and Ocean Engineering, Dalian University of Technology, Dalian 116024, China Abstract In this paper, two kinds of crack-tip singular element are employed to describe the singularity of stress-strain field near the crack tip, and then the corresponding finite element models for calculating stress intensity factor