应用三维有限单元法计算应力强度因子

1、应用三维有限单元法计算应力强度因子摘要描述了两种基于有限单元计算面形裂纹应力强度因子的方法，建议了一种创造三维有限单元网格的途径。计算方法的精度通过和其它解析解或数值解的比较得到了说明。关键词应力强度因子有限元损伤容限设计断裂评定无论在损伤容限设计还是在缺陷评定阶段，工程师们需要知道正在分析的构件中裂纹的应力强度因子，因为判断含裂纹构件的断裂，或者计算剩余疲劳寿命大多依赖于这一参量。因此，在断裂力学发展中，如何求取应力强度因子一直是一个重要的课题。当前已有许多方法可用来计算应力强度因子，较为典型的有解析法、边界配位法、有限单元法、边界元素法、体力法、权函数法和线弹簧模型。利用这些方法，大量的应

2、力强度因子解已经获得，已出版的应力强度手册1中收编了许多典型的解。尽管如此，工程师们仍然会感到自己所需要的应力强度因子解很难找到，这是因为要解决的工程问题往往是一些受复杂载荷的构件，包含的裂纹也往往是一些不规则裂纹。本文简单介绍了两种基于三维有限单元法计算面形裂纹应力强度因子的方法。有限单元法已经成为工程设计分析领域中一个强有力的计算工具，它能模拟非常复杂的构件。基于有限元的应力强度因子计算方法，自然也将具有卓越的工程能力。除了计算方法的介绍以外，还将简单描述一种简化网格的生成方法。最后提供了一些所得到的典型应力强度因子解，并和大家熟知的解进行了比较，以说明本文所描述的方法的可靠性。1计算方法

3、简介图1表示了笔者2，3建议的裂纹尖端网格形图1裂纹前沿单元网格式。网格由3个半环共12个20节点三维等参单元组成，每个半环有4个单元，其中第一个半环单元的中节点被移至1/4点位置，以模拟裂纹尖端应力和应变场的奇异性4，5。根据线弹性断裂理论68，裂纹尖端的位移场可表示为ut=0(3)式中，n和t分别为裂纹前沿的法线方向和切线方向；z为垂直于裂纹平面方向(见图2)；为材料剪切模量。图2曲线裂纹前沿坐标系1.11/4点位移法根据式(2)，若裂纹表面上(=180°)某一点的垂直于裂纹平面的位移已知，那么应力强度因子(4)用从有限单元法求出的1/4点位移uz(14)(见图1)求应力强度因子

4、的方法即为1/4点位移法，即(5)显然这一方法很简单，但要求裂纹尖端附近的应力应变场能较好地被模拟。1.2三维J积分法J积分法实际上是一种能量方法，近来被广泛用来计算应力强度因子，因为两者可按下式转换用有限元计算J积分的方法通常有2种：回路积分法与虚拟裂纹扩展法。后者最早由Parks9和Hellen10独立提出，主要是通过移动有限元模型的节点位置来模拟裂纹扩展。最近deLorenzi11，12已经根据连续介质力学成功地推导了面形裂纹的能量释放计算公式式中符号请参见文献11。商用有限元分析软件ABAQUS已经收编了上述基于“虚拟裂纹扩展原理”的三维J积分方法，可以直接获得面形裂纹前沿各角节点和中

5、节点的J积分值。2网格生成为了方便三维网格的生成，笔者建议了图3所示的分块方法。整个网格由裂纹块和非裂纹块图3结构分块示意图组成，裂纹块包含了裂纹，其网格需要进行精心划分，因为它对应力强度因子有相当大的影响。相对于裂纹块，非裂纹块显得不是很重要，一般网格可以划分得粗一些，但必须使载荷和约束有一个比较合理的分配。裂纹块和非裂纹块网格的组合可采用“多点位移约束法”，使其位移在其接触面上基本保持连续，从而使组合后的网格保持变形协调。对于复杂含裂纹构件，只要合理地设计裂纹块和非裂纹块即可大大简化网格的生成。这种网格形式也特别适合于疲劳裂纹扩展的模拟计算，疲劳裂纹扩展后的新网格重新生成可以只对裂纹进行，

6、而保持非裂纹块网格不变2，3。图4表示了一个笔者建议的三维裂纹块网格的生成过程。首先生成一个8节点等参单元二维网格，在这个网格里定义开裂表面和裂纹前沿，然后将其扩展成一个三维网格，单元的扩展有2种，见图4。一个非裂纹前沿单元扩展成一个20节点单元，而一个裂纹前沿单元扩展成6个20节点单元。围绕裂纹前沿，网格呈放射形。最后组合成的裂纹尖端网格具有图1所示的结构。显然，这种网格扩展方法进一步简化了三维网格的建立。图4三维裂纹块网格的形成3计算程序笔者最近发展了如图5所示的计算应力强度因子的计算机程序。对于一些常见的工程裂纹，如平板中的半椭圆表面裂纹、压力容器中的半椭圆表面裂纹、平板孔边的角裂纹、圆

7、棒中的表面裂纹等，只要输入一些关键的几何参数和载荷参数，程序将自动输出应力强度因子结果。该程序目前主要和商用软件ABAQUS相连接，在UNIX操作系统下运行。图5应力强度因子计算流程图4结果比较4.1无限大体中的椭圆形埋藏裂纹图6表示了沿椭圆形埋藏裂纹前沿的应力强度因子结果。椭圆裂纹短长轴之比为a/c=0.5。结果表明，根据位移法和J积分法获得的结果非常相近。两者的误差约为1%。图中的虚线为Irwin6的解析解，本文结果略小于Irwin的解，最大误差小于3%。图6埋藏椭圆裂纹前沿应力强度因子变化4.2有限厚度平板中的半椭圆表面裂纹近20年来，这一裂纹得到了深入的研究，大量的应力强度因子结果已经

8、发表13，其中，Newman_Raju14根据他们自己的有限单元分析结果提出的应力强度因子方程，被认为具有比较可靠的精度，并且易编程。这一方程已经在工程中得到了广泛的应用。图7比较了笔者和Newman_Raju的结果，笔者的结果由1/4点位移法获得。裂纹形状为半圆形a/c=1。结果表明，无论在裂纹深度点还是在表面点，笔者的结果和NewmanRaju的解误差约在5%以内，全面的比较见文献2。图7平板中半椭圆表面裂纹前沿应力强度因子比较(a)拉伸载荷(b)弯曲载荷4.3压力容器中的半椭圆形内表面裂纹图8所示为压力容器筒体的局部网格及应力强度因子结果。筒体的壁厚和内径之比为t/Ri=0.1，载荷为均

9、匀内压力。RajuNewman15的结果也示于图8中。结果表明，笔者的结果和RajuNewman的结果基本一致，只是当裂纹深度达到4/5的壁厚时，RajuNewman的结果相对低一些。从图8a中可以看出，筒体网格由裂纹块和非裂纹块组成。进一步的比较见文献16。图8受内压容器内表面裂纹(a)网格(b)应力强度因子比较4.4孔边的半椭圆表面裂纹图9比较了笔者的结果17和ShivarkumarNewman18及Zhao等人19，20的结果。ShivarkumarNewman的结果基于有限元分析，而Zhao等人用的是权函数方法。可见，结果同样符合得较好。图9孔边表面裂纹应力强度因子比较5结束语本文描述了两种基于有限单元法计算面形应力强