基于有限元法对裂纹尖端应力强度因子的计算分析

袁浩，李菁，谢禹钧，侯汶雨

(1. 辽宁石油化工大学 机械工程学院，辽宁 抚顺 113001； 2. 黑龙江省北方工具有限公司 机械制造公司，黑龙江 牡丹江 157000； 3. 中国石油天然气股份有限公司 抚顺石化分公司，辽宁 抚顺113001)

0 引言

由于加工工艺以及材料本身缺陷等原因，结构构件中会不可避免的存在缺陷或裂纹，从而影响结构系统的安全性，甚至会引起结构失效。通过对大量破坏事故研究表明，这些裂纹的存在会导致含裂纹构件的实际强度远低于理论强度[1]。传统的强度观点通常把材料视为理想材料，即材料是连续、均匀、各向同性的，但实际工况中材料很难达到理想状态[2]。为了确保含裂纹构件长期稳定地安全运行，必须对不可避免存在的裂纹对构件的影响进行预判，从而将发生损失的风险降至最低。在断裂力学问题的分析中，应力强度因子K是预判含裂纹构件发生断裂和裂纹发生扩展速率的首要判据[3-4]。获得应力强度因子的方法大致上可分为解析法、数值法和实验法。有限元数值法以计算机为平台，利用计算机的计算能力和强大的建模能力，可以解决工程中复杂的几何条件和边界条件下的实际问题，而且有限元法不仅局限于线弹性问题，在研究弹塑性断裂力学、疲劳和蠕变裂纹扩展速率等问题方面也同样适用，已经成为获得应力强度因子的主要途径。

本文以薄板中心穿透裂纹为研究对象，利用ANSYS为平台建立裂纹模型，分别采用J积分法、位移外推法、相互作用积分法等3种方法计算裂纹尖端的应力强度因子K，并对3种方法求出的结果进行比较分析。还研究了受力、裂纹长度、含裂纹薄板的几何参数等因素对应力强度因子的影响。

1 含裂纹薄板有限元模型的建模方法

在以往的研究中，大多数裂纹的有限元模型的建立采用的是对称半裂纹的建模方法，本文采用8节点四边形实体PLANE183单元，通过创建位置重合的关键点模拟裂纹尖端的方法建立全裂纹模型。

裂纹尖端附近应力场表达式为[4]：

(1)

图1 正常单元退化为奇异单元

图2 裂纹尖端奇异单元排列

需要注意的是，在用相互作用积分法和位移外推法计算K值时，裂纹面要与局部坐标x轴方向保持一致，裂纹面法相向与局部坐标系y轴保持一致[7]。在运用J积分法计算应力强度因子时要保证裂纹面与全局坐标系的x轴保持一致。而在裂纹尖端单元的选择上，除J积分法在裂纹尖端不必采用奇异单元划分网格外，其余两种方法均需要在裂纹尖端采用奇异单元[8]。

2 计算应力强度因子的几种方法

由式(1)可知，裂纹尖端应力场的强弱主要由K值来描述，它的大小可以说明裂尖附近区域的安全程度，是预防含裂纹构件的失效和评估结构安全寿命的重要参量。对于应力强度因子的计算通常有数值法和解析法，由于数值法中的有限元法有计算方便、精度高等特点，如今已被大量地运用在断裂参量的计算中。根据有限元思想的基本原理，有相互作用积分法、位移外推法和J积分法等多种方法可被用于应力强度因子的计算中。

2.1 相互作用积分法

Wang S.S等[9]首先运用相互作用积分法计算了混合型裂纹的应力强度因子。相互作用积分法与J积分法类似，都是通过对积分区域进行面积分或体积分来获取应力强度因子的方法，其定义式为:

(2)

(3)

由式(3)得，真实场和辅助场共同作用的J积分为：

(4)

减除出真实场和辅助场的J积分项可得附加项式(2)与应力强度因子之间的关系式为：

(5)

(6)

2.2 位移外推法

另一种方法是通过裂纹尖端节点位移间接推出应力强度因子，根据Ⅰ型裂纹尖端附近的位移表达式[1]：

(7)

式中，u1=μ,u2=υ，E为弹性模量，r为节点到裂纹尖端的距离。

k=3-4v平面应变状态

从式(7)可以看出，应用有限元法求出ui后，便可求得应力强度因子KⅠ。

(8)

由于节点的位移在裂纹尖端处最为明显，所以式(8)只有在裂纹尖端(r→0)时才能获得较为准确的结果，但当r→0时，KⅠ→∞，这在数值模拟中很难完成的。但在r很小的范围内，可近似地把r与KⅠ看成线性关系。因此，需要沿裂纹面选取多个较小的r值计算节点位移，从而得到多个相应的KⅠ，然后将这些KⅠ值运用曲线拟合或插值的方法最终外推到r→0，便可得到裂纹尖端的KⅠ值。Chen[10]通过计算出图1的H、D两1/4分点处的应力强度因子进行插值而得到更为精确的应力强度因子表达式：

(9)

式中,uH、uD为H、D两节点的位移;L为裂纹面单元的边长。式(9)较于式(8)，将裂纹尖端附近1/4点处节点作为一个插值点，与裂纹尖端附近应力和应变场奇异的特点相符合，因此式(9)的计算结果更为精确。这种由裂纹尖端附近的节点位移间接推导出裂纹尖端应力强度因子的方法就是位移外推法。

2.3 J积分法

对于二维问题，Rice[4]提出了与路径无关的J积分回路积分(图3)。J积分与应力强度因子K一样具有反映场强的性质，其值可以表征裂纹尖端应力应变场的强弱。张巨伟等[11]也证明了J积分在用ANSYS计算应力强度因子中的可行性和准确性。它的定义表达式为：

(10)

式中：Γ为起始于裂纹下表面，终止于裂纹下表面的任意一条围绕裂纹尖端的逆时针回路；Ti为回路Γ上任意一点处的应力分量；W为回路Γ上任意一点处的应变能密度；n为回路Γ上外法线方向的单位矢量；ui为回路Γ上任意一点处的位移分量；ds为回路Γ上的弧元。经过推导证明，在线弹性的情况下，J积分的值与裂纹扩展能量释放率GⅠ、应力强度因子K之间有如下关系：

图3 J积分回路示意图

(11)

运用ANSYS软件的后处理功能，经过映射路径等一系列ANSYS指令计算出裂纹尖端的J积分，然后将所得的 J积分值代入到式(11)中便得到裂纹尖端的应力强度因子。

3 计算实例

如图4所示一含有中心裂纹长、宽分别为2H=120mm、2W=100mm的薄板，其裂纹长为2a=20mm，材料特性为：弹性模量E=2×1011MPa，泊松比υ=0.3。薄板承受均匀载荷σ=10MPa。

在ANSYS环境中，建立全裂纹模型，建模中采用PLANE183单元，该单元为8节点二次实体单元，在裂纹尖端附近则采用退化的PLANE183单元，如图5、图6所示。图中能清楚地看出围绕在裂纹尖端的是三角形的奇异单元，以及经过退化后在1/4分点处的中间节点，而在远离裂纹尖端处则采用正常的四边形PLANE183单元。通过利用ANSYS求解器以及通用后处理的功能，运用GUI和命令流的操作便可获得分别用3种方法就算出的KⅠ值。需要说明的是，当使用相互作用积分法计算KⅠ时，由于ANSYS14.0暂无相互作用积分相应的GUI操作，所以只能采用命令流的方式求解。

图4 有限宽中心裂纹模型

图5 全裂纹模型的网格划分

图6 裂纹尖端的奇异单元

为了研究载荷和裂纹模型几何尺寸对上述3种求取KⅠ方法的影响，本文采用控制变量的方法分别以不同的载荷、裂纹半长、薄板长度以及薄板宽度作为变量进行求解，对比相互作用积分法、位移外推法、J积分法的计算结果，并根据应力强度因子手册[12]中理论值表达式：

(12)

计算相对误差，以得到更为直观的结果，如表(1)-表(4)所示，其中KⅠ的单位为MPa·m1/2,表中“()”内为与理论值的相对误差。

表1 不同载荷条件KⅠ的值

表2 不同裂纹长度条件下KⅠ的值

表3 不同薄板宽度条件下KⅠ的值

表4 不同薄板长度条件下KⅠ的值

图7 不同边界条件下应力强度因子KⅠ曲线比较

4 结语

1) 通过分析J积分法、位移外推法、相互作用积分法3种计算应力强度因子的方法，发现3种方法与解析法理论值的相对误差都在合理范围内，证明了以上3种方法在计算应力强度因子方面的可行性与准确性。

2) 3种方法中，J积分法与相互作用积分法较位移外推法求得的应力强度因子结果与理论值的相对误差较小、精度较高，但由于J积分法需要定义积分路径，所以相互作用积分法在效率方面要优于J积分法。因此，在选择应用有限元法原理计算应力强度因子的这3种方法中，相互作用积分法更有优势。

[1] 程靳，赵树山. 断裂力学[M]. 北京：科学出版社，2006.

[2] 丁遂栋，孙利民. 断裂力学[M]. 北京：机械工业出版社，1997.

[3] Lin X B, Smith R A. Finite Element Modelling of Fatigue Crack Growth of Surface C racked Plates[J]. Engineering Fracture Mechanics , 1999, 63 : 503-522 .

[4] Rice J R. A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks[J]. J ApplMech , 1968, 35 : 379-386 .

[5] 周益春, 郑学军. 材料的宏微观力学性能[M]. 北京: 高等教育出版社, 2009 .

[6] Fehl B D, Truman K Z. An Evaluation of Fracture Mechanics Quarter Point Displacement Techniques Used for Computing S tress Intensity F actor s[J]. Engineering Structures, 1999, 21 : 406-415 .

[7] 宫经全，张少钦，李禾，等. 基于相互作用积分法的应力强度因子计算[J]. 南昌航空大学学报(自然科学版)，2015，29(1)：45-48.

[8] 王伟. Ⅰ- Ⅱ复合型裂纹 J 积分计算[J]. 石油化工设备，2011，40(2)：55-57.

[9] Wang S S，Yau J F，Corten H T． A mixed-mode crack analysis of rectilinear anisotropic solids using conservation laws of elasticity[J]. International Journal of Fracture，1980，16( 3) : 247-259.

[10] Chen L S, Kuang J H. A Modified Linear Extrapolation Formula for Determination of Stress Intensity Factor[J]. Int Journ of Fracture, 1992, 54 : 3-8 .

[11] 金晓娜, 张巨伟. J 积分在二维断裂算例中的计算分析[J]. 机械工程与自动化，2008,149(4):57-59.