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

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

《基于奇异元计算分析裂纹尖端应力强度因子.docx》由会员分享，可在线阅读，更多相关《基于奇异元计算分析裂纹尖端应力强度因子.docx（14页珍藏版）》请在冰豆网上搜索。

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

51卷第3期（总第192期）中国造船Vol.51No.3（SerialNo.1922010年9月SHIPBUILDINGOFCHINASep.2010

文章编号：

1000-4882（201003-0056-09

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

陈景杰，黄一，刘刚

（大连理工大学船舶工程学院,辽宁大连116024

摘要

采用两种奇异单元模拟裂纹尖端应力应变场的奇异性，建立了相应的计算裂纹尖端应力强度因子的ANSYS有限元模型。

通过数值计算，分别考察了这两种有限元模型中裂纹尖端附近区域网格参数的变化对应力强度因子计算精度的影响，比较了应力强度因子对各个参数的敏感程度。

发现采用20节点奇异元的有限元模型计算的应力强度因子几乎与网格参数无关，其计算结果更稳定可靠。

该模型能够用于船舶及海洋工程结构中含裂纹构件的应力强度因子计算。

关键词：

船舶、船舶工程;应力强度因子;奇异单元;有限元模型;网格参数中图分类号：

U661.4文献标识码：

A

0引言

在复杂的海洋环境载荷作用下，船舶及海洋工程结构物在其服役期内会发生结构疲劳损伤[1]，而关键承载构件上疲劳裂纹的出现是结构疲劳损伤的必然结果。

对于那些含有疲劳裂纹的结构，应力强度因子是评估它们断裂失效的重要参量，在裂纹体分析中占据着重要地位。

自断裂力学问世以来，已经产生了众多理论和数值解法求解应力强度因子，较为典型的计算方法有有限元法、边界元法、解析法、混合法、权函数法和线弹簧法等。

利用这些方法，已经获得许多裂纹体模型的应力强度因子的大小，应力强度因子手册[2]中收编了许多种典型裂纹体模型应力强度因子的解，但由于裂纹体几何形式及所受载荷的复杂性，很多情况下应力强度因子的解难以从现有手册查到。

而有限元方法不受裂纹体几何及载荷形式的限制，因而在断裂力学中得到了非常广泛的应用。

位移法是用有限元法计算裂纹尖端应力强度因子的常用方法，它是利用奇异单元来精确描述线弹性范围内裂纹尖端的应力和应变场的奇异性。

取裂纹尖端第一层单元为奇异元，其余各层单元为相应等参元，用这样的方法建立裂纹体的有限元模型[3~11]来计算裂纹尖端的应力强度因子。

目前有两种奇异单元被广泛使用，分别为20节点奇异元和12节点奇异元，它们是由20节点等参元蜕化、转变获得，形式如图1所示。

文中以ANSYS软件为基础，分别采用20节点和12节点奇异单元建立了含表面裂纹的平板的两种有限元模型。

对有限元网格技术及裂纹体网格剖分的要求作了深入研究，评估了两种有限元模型计算的应力强度因子对网格的敏感程度，进而判断出哪种有限元模型计算的应力强度因子更稳定可靠，这为有限元法确定船舶及海洋工程中的局部复杂结构裂纹尖端应力强度因子提供了依据。

收稿日期：

2009-12-19；修改稿收稿日期：

2010-04-08

51卷第3期（总第192期陈景杰，等：

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

（a）20节点的奇异元（b）12节点的奇异元（c）20节点奇异元围绕裂尖形式

图1奇异单元类型及其围绕裂尖形式

1有限元模型及网格参数

对图2所示的含半椭圆表面裂纹的平板模型进行有限元计算和分析，在有限元模型中考虑了裂纹前缘周围区域网格的变化对应力强度因子的影响，分析不同奇异单元所建的有限元模型计算结果的稳定性。

含有表面裂纹的有限平板示于图2，图中的各个参数列于表1。

在平板长度两端受到垂直裂纹面方向的均匀拉伸载荷作用。

根据模型及载荷的对称性采用奇异单元法建立1/4有限元模型，如图3所示。

图4给出了图3中裂纹尖端附近区域网格剖分情况、裂纹尖端端面形状及影响应力强度因子精度的主要参数，参数含义见表1。

图2含有表面裂纹的有限平板几何模型图3有限平板的有限元模型

（a）沿半裂纹尖端曲线（b）裂纹深度方向横截面区域附近的有限元模型num_crack=16的有限元形状，α=45°,num_tip

=5

（c）裂纹尖端半端面的局部放大图图4裂纹尖端附近区域的有限元模型

1/4节点中间节点角节点

裂纹尖端

58中国造船学术论文表1含有表面裂纹的平板的主要参数

符号参数2L平板长度2W平板宽度T平板厚度2a表面裂纹长度b表面裂纹最大深度Ltip

裂纹尖端附近区域大小α

裂纹尖端单元角度num_tip围绕裂纹尖端单元层数num_crack裂纹尖端曲线长度所分的份数Lin裂纹尖端奇异单元长度L

在Ltip内最外层单元长度

2网格参数变化对计算结果影响

对于采用奇异单元建立的有限元模型，分析裂纹尖端附近网格变化对计算结果的影响。

在这方面，目前主要是研究了围绕裂纹前缘的单元层数num_tip及奇异单元尺寸Lin与裂纹深度b比这两个因素的影响。

关于num_tip，文献[3-8]都采用了3层以上的网格，文献[4]更是明确建议裂纹前缘网格层数不能少于3层。

关于Lin的选择问题，文献[4]建议Lin/b应该小于0.1，文献[3]建议该比值应该介于0.05和0.15之间，而有限元分析软件ANSYS的技术文件则建议Lin/b小于1/8。

本文不仅研究了上述两个因素对应力强度因子的影响，而且进一步研究裂纹尖端单元角度α、裂纹尖端附近区域大小Ltip和裂纹尖端曲线的份数num_crack对该因子的影响。

2.1

α的变化对应力强度因子的影响

在图2所示的几何模型中，取2a=40mm，b=10mm，2W=120mm，2L=200mm，T=20mm，均匀拉应力为100MPa，建立两种奇异元的有限元模型。

在num_tip=5、Lout/Lin=1、Ltip/d=1/10和num_crack=16的两种有限元网格中，研究裂纹尖端附近区域α分别为9°、15°、22.5°、45°时单元角度的变化对有限元模型计算的应力强度因子的影响。

图5给出了不同α对应的裂纹尖端附近区域的横截面。

根据上述的模型尺寸及网格参数，两种有限元模型计算的应力强度因子如图6和图7所示，它们给出了20节点和12节点奇异元计算的应力强度因子对α的敏感程度。

在图6和图7中SIF表示裂纹最深位置尖端的应力强度因子，在后面各图中的SIF具有相同含义。

根据图6中裂纹尖端应力强度因子随单元角度α的变化情况，不难得出在20节点奇异元的有限元模型中，裂纹尖端角度的变化对应力强度因子几乎没有影响，为了提高计算效率，在后面对于20节点奇异元的有限模型只采用α=45°来计算裂纹尖端的应力强度因子。

（a）α

=45°（b）α=22.5°（c）α=15°（d）α=9°

图5α取不同值时裂纹尖端附近区域的横截面

51卷第3期（总第192期陈景杰，等：

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

图6不同α值时20节点奇异元计算的图7不同α值时12节点奇异元计算的

应力强度因子应力强度因子

图7给出了12节点奇异元的有限元模型计算的应力强度因子随着α的增加而不断减小的变化趋势，说明裂纹尖端角度的变化对该模型计算结果的影响不可忽略。

结合计算精度和计算效率，在后面对于12节点奇异元的有限元模型应取裂纹尖端单元角度α=22.5°来计算裂纹尖端的应力强度因子。

2.2num_tip的变化对应力强度因子的影响

几何模型与2.1节相同，在20节点和12节点奇异元的有限元网格中，取Lout/Lin=1、Ltip/d=1/10和num_crack=16，而α分别为45°和22.5°来分析围绕裂纹尖端单元层数的变化对有限元模型计算的应力强度因子的影响。

文中考虑了num_tip=3、5、8、20的情况。

20节点奇异元模型对应的裂纹尖端附近区域的横截面形状如图8所示，不同层数下两种有限元模型计算的应力强度因子如图9~10所示。

（a）num_tip=3（b）num_tip=5（c）num_tip=8（d）num_tip=20

图8num_tip取不同值时裂纹尖端附近区域的横截面形状

图9不同num_tip值20节点奇异元计算的图10不同num_tip值12节点奇异元计算的

应力强度因子应力强度因子

图9给出了不同num_tip下20节点奇异元有限元模型计算的应力强度因子的变化情况，说明围绕裂纹尖端单元层数的变化对应力强度因子的影响可以忽略，在计算裂纹尖端的应力强度因子有限元模型中num_tip可取3~20之间的任意值。

60中国造船学术论文对于12节点奇异元有限元模型计算的应力强度因子（SIF）结果见表2。

可以看到当num_tip=20时，应力强度因子值有较大的变化。

为了进一步分析围绕裂纹尖端的单元层数对有限元模型计算的应力强度因子的影响，又取num_tip=11、15、18、19的情况进行计算，计算结果见图10。

表2不同num_tip值下12节点奇异元模型计算的应力强度因子单位：

1/2MPamm⋅

num_tip

um_tip

SIF

20576.6983

从图10可以看出，当num_tip>18时，围绕裂纹尖端单元层数的变化对计算的应力强度因子有一定的影响，而num_tip在3~18的范围内，采用12节点奇异元的有限元模型计算的应力强度因子几乎相等，说明采用12节点奇异元的有限元模型计算裂纹尖端应力强度因子时，num_tip可取3~18的范围内的任意值。

2.3Lout/Lin的变化对应力强度因子的影响

几何模型不变，基于上述研究结果，在20节点和12节点奇异元的有限元模型中取num_tip=5，α分别为45°和22.5°，同时假定Ltip/d=1/10和num_crack=16，研究不同Lout/

Lin下两种有限元模型计算应力强

度因子的变化情况，结果如图

12和图

13所示。

图11给出了Lout/Lin分别为1、3、5和10时有限元模型中裂纹尖端附近区域的横截面形状。

（a）Lout/Lin=1（b）Lout/Lin=3（c）Lout/Lin=5（d）Lout/Lin=10

图11L

out/Lin取不同值时裂纹尖端附近区域的横截面形状

图12不同Lout/Lin值时20节点奇异元计算的图13不同Lout/Lin值时12节点奇异元计算的

应力强度因子应力强度因子

图12和图13分别描述了20节点和12节点奇异元对应的有限元模型计算的应力强度因子随Lout/Lin的变化情况。

根据图12中的计算结果，最大值与最小值的相对误差为0.35%，说明不同Lout/Lin的变化对于20节点有限元模型计算的应力强度因子几乎没有影响，采用该模型计算应力强度因子时，Lout/Lin可取1~30内的任意值。

后面取Lout/Lin=1对应的有限元模型计算裂纹尖端的应力强度因子。

在图13中，采用12节点奇异元的有限元模型计算的应力强度因子具有随着Lout/Lin的增大而增大

51卷第3期（总第192期陈景杰，等：

基于奇异元计算分析裂纹尖端应力强度因子61的趋势，说明这种有限元模型计算的应力强度因子对Lout/Lin参数具有较大的灵敏性，不同的网格剖分形式会产生不同的计算的结果，导致计算结果不可靠。

与Newman-Raju法[1]计算的应力强度因子进行比较，当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/Lin=4，在取nL同时两种模型中令num_crack=16，研究Ltip/b在1/80~1/4范围内的变化时对有限元模型计算的应力强度因子的影响，结果如图14和图15所示。

590600595SIF/MPa·mm1/22.4SIF/MPa·mm1/25855905855805755805755700.000.050.100.150.200.250.30Ltip/b图14不同Ltip/b值时20节点奇异元计算的应力强度因子0.000.050.100.15Ltip/b0.200.250.30图15不同Ltip/b值12时节点奇异元计算的应力强度因子从图14可以看出，随着Ltip/b的增加，20节点奇异元有限元模型计算的应力强度因子几乎不变，说明Ltip的变化对有限元模型计算裂纹尖端应力强度因子几乎没有影响，Ltip/b可以取1/40~1/4范围内的任意值。

图15给出了在不同Ltip/b下，12节点奇异元的有限元模型计算的应力强度因子的变化曲线。

当Ltip/b<0.05时，计算的应力强度因子随之Ltip/b的减小迅速增加，为了使计算结果稳定，Ltip/b应取0.05~0.25范围内的值进行裂纹尖端应力强度因子的计算。

2.5num_crack的变化对应力强度因子的影响几何模型不变，在20节点奇异元的有限元模型中裂纹尖端附近区域的网格参数分别为α=45°，L在nnum_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进行分析，计算结果如图16~17所示。

图16~17给出了不同num_crack值对20节点和12节点奇异元模型计算应力强度因子的影响。

图16说明在20节点奇异元模型中，num_crack值的变化对应力强度因子没有影响，在模型剖分中可取任意值。

而图17给出在12节点奇异元有限元模型中，num_crack的变化对应力强度因子具有一定影响，当num_crack>12时，影响趋势减缓，结合计算效率，对于12节点奇异元的有限元模型应在12~20的范围内对num_crack取值来计算应力强度因子。

62中国造船学术论文59058558057557051015num_crack图16不同num_crack值时20节点奇异元计算的应力强度因子5905855805755702025SIF/MPa·mm1/2SIF/MPa·mm1/251015num_crack2025图17不同num_crack值时12节点奇异元计算的应力强度因子综上所述，12节点奇异元具有较少的节点，并且仅需与8节点单元相连接，大大减少了有限元模型中的节点个数，能提高计算效率，但经过上述分析可知，这种奇异单元的有限元模型计算应力强度因子的稳定性与裂纹尖端附近区域网格形式有关，导致计算结果不可靠。

而20节奇异元虽然节点数较多，但该模型计算的应力强度因子十分稳定，它不受有限元模型中网格剖分形式的影响。

对于图2所示的模型，采用20节点奇异元有限元模型计算的应力强度因子为580.561MPa⋅mm1/2与Newman-Raju二者相对误差为2.68%。

这说明计算裂纹尖端应力强度法[12]计算的结果596.548MPa⋅mm1/2进行比较，因子的20节点奇异元有限元模型不仅具有很好的稳定性而且具有较高的计算精度，它的计算结果是正确可靠的。

3结论本文采用有限元数值计算，研究了20节点和12节点奇异元这两种有限元模型中网格参数的变化对应力强度因子的影响。

通过分析比较可得到如下结论：

在计算裂纹尖端应力强度因子的有限元模型中，20节点奇异元的有限元模型是稳定可靠的，它所计算的应力强度因子与裂纹尖端附近区域的网格状态无关，在数值模拟分析中能够得到稳定的计算结果；12节点奇异元的有限元模型计算结果与裂纹尖端附近区域的网格状态密切相关，对于相同问题，在数值模拟中不同剖分形式得到的计算结果差别较大，显示该模型在应用中不够稳定，难以处理复杂结构的裂纹问题。

因此对于船舶及海洋工程结构物等复杂结构的裂纹问题宜采用20节点奇异元的有限元模型来确定裂纹尖端的应力强度因子，以便于评估结构所处的安全状态。

参考文献：

[1]王一飞,衣高洁,黄小平,崔维成.疲劳裂纹扩展模型中表征裂纹闭合水平参数的确定[J].中国造船,2008,49（4:

45-51.[2][3]中国航空研究院.应力强度因子手册[M].北京:

科学出版社,1993.INGRAFFEAAR,MANUC.Stressintensityfactorcomputationinthreedimensionswithquarter-pointelements[J].InternationalJournalforNumericalMethodsinEngineering,1980,15:

1427-1445.[4]CAOJJ,YANGGJ,PAKERJA,BURDEKINFM.CrackmodelinginFEanalysisofcirculartubularjoints[J].

51卷第3期（总第192期陈景杰，等：

基于奇异元计算分析裂纹尖端应力强度因子63EngineeringFractureMechanics,1998,61:

537-553.[5]CHIEWSP,LIEST,LEECK,HUANGZW.StressintensityfactorsforasurfacecrackinatubularT-joint[J].InternationalJournalofPressureVesselsandPiping,2001,78:

677-685.[6]BOWNESSD,LEEMMK.Predictionofweldtoemagnificationfactorsforsemi-ellipticalcracksinT-buttjoints[J].InternationalJournalofFatigue,2000,22:

369-387.[7]LEEMMK,BOWNESSD.Fatigueanalysisoftubularjointsusingfiniteelements[J].FatigueinOffshoreStructures,1996,1:

223-257.[8]SEETHARAMANS,SREEDHARDS.Fatigueanalysisandfracturemechanics[J].FatigueinOffshoreStructures,1996,1:

259-309.[9]WANGX,LAMBENTSB.Onthecalculationofstressintensityfactorsforsurfacecracksinweldedpipe-plateandtubularjoints[J].InternationalJournalofFatigue,2003,25:

89-96.[10]卢炎麟,王自力,杜水友.T型管节点焊趾表面裂纹扩展行为的三维断裂力学分析[J].中国造船,1995,2:

51-59.[11]BOWNESSD,LEEMMK.Thedevelopmentofanaccuratemodelforthefatigueassessmentofdoublycurvedcracksintubularjoints[J].InternationalJournalofFatigue,1995,73:

129-147.[12]NEWMANJC.RAJUIS.AnEmpiricalStress-IntensityFactorEquationfortheSurfaceCrack[J].EngineeringFractureMechanics,1981,15（1-2:

185-192.AnalysisofFiniteElementModelforCalculatingStressIntensityFactorBasedonCrack-tipSingularElementCHENJing-jieHUANGYiLIUGang（SchoolofNavalArchitectureandOceanEngineering,DalianUniversityofTechnology,Dalian116024,ChinaAbstractInthispaper,twokindsofcrack-tipsingularelementareemployedtodescribethesingularityofstress-strainfieldnearthecracktip,andthenthecorrespondingfiniteelementmodelsforcalculatingstressintensityfactorsatthecracktipareestablishedbyusingANSYS.Inthisnumericalsimulation,theinfluencesofvariousgridparametersintheregionnearcracktipsonthecalculationaccuracyofthestressintensityfactorsareinvestigatedrespectivelyforthesetwokindsoffiniteelementmodels.Thentheparameters’sensitivitydegreestothestressintensityfactorsareinvestigatedindetail.Thestressi