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

应用三维有限单元法计算应力强度因子 摘要　描述了两种基于有限单元计算面形裂纹应力强度因子的方法，建议了一种创造三维有限单元网格的途径。计算方法的精度通过和其它解析解或数值解的比较得到了说明。　　关键词　应力强度因子　有限元　损伤容限设计　断裂评定　　 　　无论在损伤容限设计还是在缺陷评定阶段，工程师们需要知道正在分析的构件中裂纹的应力强度因子，因为判断含裂纹构件的断裂，或者计算剩余疲劳寿命大多依赖于这一参量。因此，在断裂力学发展中，如何求取应力强度因子一直是一个重要的课题。当前已有许多方法可用来计算应力强度因子，较为典型的有解析法、边界配位法、有限单元法、边界元素法、体力法、权函数法和线弹簧模型。利用这些方法，大量的应力强度因子解已经获得，已出版的应力强度手册［1］中收编了许多典型的解。尽管如此，工程师们仍然会感到自己所需要的应力强度因子解很难找到，这是因为要解决的工程问题往往是一些受复杂载荷的构件，包含的裂纹也往往是一些不规则裂纹。　　本文简单介绍了两种基于三维有限单元法计算面形裂纹应力强度因子的方法。有限单元法已经成为工程设计分析领域中一个强有力的计算工具，它能模拟非常复杂的构件。基于有限元的应力强度因子计算方法，自然也将具有卓越的工程能力。除了计算方法的介绍以外，还将简单描述一种简化网格的生成方法。最后提供了一些所得到的典型应力强度因子解，并和大家熟知的解进行了比较，以说明本文所描述的方法的可靠性。1　计算方法简介　　图1表示了笔者［2，3］建议的裂纹尖端网格形 图1　裂纹前沿单元网格 式。网格由3个半环共12个20节点三维等参单元组成，每个半环有4个单元，其中第一个半环单元的中节点被移至1/4点位置，以模拟裂纹尖端应力和应变场的奇异性［4，5］。根据线弹性断裂理论［6～8］，裂纹尖端的位移场可表示为 ut=0　　(3) 式中，n和t分别为裂纹前沿的法线方向和切线方向；z为垂直于裂纹平面方向(见图2)；μ为材料剪切模量。 图2　曲线裂纹前沿坐标系 1.1　1/4点位移法　　根据式(2)，若裂纹表面上(θ=180°)某一点的垂直于裂纹平面的位移已知，那么应力强度因子 　　(4) 用从有限单元法求出的1/4点位移uz(1／4)(见图1)求应力强度因子的方法即为1/4点位移法，即 　　(5) 显然这一方法很简单，但要求裂纹尖端附近的应力应变场能较好地被模拟。1.2　三维J积分法　　J积分法实际上是一种能量方法，近来被广泛用来计算应力强度因子，因为两者可按下式转换 用有限元计算J积分的方法通常有2种：回路积分法与虚拟裂纹扩展法。后者最早由Parks［9］和Hellen［10］独立提出，主要是通过移动有限元模型的节点位置来模拟裂纹扩展。最近deLorenzi［11，12］已经根据连续介质力学成功地推导了面形裂纹的能量释放计算公式 式中符号请参见文献［11］。商用有限元分析软件ABAQUS已经收编了上述基于“虚拟裂纹扩展原理”的三维J积分方法，可以直接获得面形裂纹前沿各角节点和中节点的J积分值。2　网格生成　　为了方便三维网格的生成，笔者建议了图3所示的分块方法。整个网格由裂纹块和非裂纹块 图3　结构分块示意图 组成，裂纹块包含了裂纹，其网格需要进行精心划分，因为它对应力强度因子有相当大的影响。相对于裂纹块，非裂纹块显得不是很重要，一般网格可以划分得粗一些，但必须使载荷和约束有一个比较合理的分配。裂纹块和非裂纹块网格的组合可采用“多点位移约束法”，使其位移在其接触面上基本保持连续，从而使组合后的网格保持变形协调。对于复杂含裂纹构件，只要合理地设计裂纹块和非裂纹块即可大大简化网格的生成。这种网格形式也特别适合于疲劳裂纹扩展的模拟计算，疲劳裂纹扩展后的新网格重新生成可以只对裂纹进行，而保持非裂纹块网格不变［2，3］。　　图4表示了一个笔者建议的三维裂纹块网格的生成过程。首先生成一个8节点等参单元二维网格，在这个网格里定义开裂表面和裂纹前沿，然后将其扩展成一个三维网格，单元的扩展有2种，见图4。一个非裂纹前沿单元扩展成一个20节点单元，而一个裂纹前沿单元扩展成6个20节点单元。围绕裂纹前沿，网格呈放射形。最后组合成的裂纹尖端网格具有图1所示的结构。显然，这种网格扩展方法进一步简化了三维网格的建立。 图4　三维裂纹块网格的形成 3　计算程序　　笔者最近发展了如图5所示的计算应力强度因子的计算机程序。对于一些常见的工程裂纹，如平板中的半椭圆表面裂纹、压力容器中的半椭圆表面裂纹、平板孔边的角裂纹、圆棒中的表面裂纹等，只要输入一些关键的几何参数和载荷参数，程序将自动输出应力强度因子结果。该程序目前主要和商用软件ABAQUS相连接，在UNIX操作系统下运行。 图5　应力强度因子计算流程图 4