基于相互积分的应力强度因子数值分析

乐京霞 周 恒

（武汉理工大学交通学院 武汉 430063）

0 引 言

80%以上的结构严重破坏是由于疲劳和断裂破坏引起的［1］．对于断裂问题，常采用应力强度因子（stress intensity factor，SIF）这一物理量来反映裂纹尖端弹性应力场的强弱．传统求解SIF方法很多，如应力函数法、积分变换法、有限单元法、边界元法和边界配置法等．以上方法理论成熟，但需要繁琐的推导和计算．文献［2］给出了部分简单裂纹模型的SIF计算公式．然而真实构件中的裂纹形态分布各异，手册不可能完全覆盖．因此，利用商用有限元分析软件对裂纹尖端应力场进行数值模拟进而计算SIF的方法，由于其简便易行，逐渐成为工程中解决疲劳断裂问题常用方法．

ANSYS求解断裂力学问题的步骤为：首先进行线弹性分析或弹塑性静力分析，然后用特殊的后处理命令、或宏命令计算所需的断裂参数．其自带的KCALC命令可直接由裂尖应力／位移场插值求解出SIF．在新推出的高版本ANSYS中增加了一种基于相互积分理论的数值计算命令CINT．本文利用ANSYS有限元分析软件，基于上述两种不同的断裂力学理论模型，分别求解了张开型二维裂纹及三维裂纹的SIF值，并将数值计算结果与理论解析解和裂纹扩展实验结果进行比较，得出了更为准确的张开型裂纹应力强度因子数值计算方法．

1 基于相互积分法理论的应力强度因子计算

相互作用积分法假定含裂纹弹性介质受到2个载荷的共同作用：有单独真实载荷引起的裂纹尖端力学场，称为真实应力σki，σkj；uk，i为真实场中的位移；σauxkj称 为 辅 助 场 的应力；uauxk，i为 辅 助场中的位移；qi，j为裂纹扩展矢量；qn为裂纹扩展法向；δi，j为克罗次符号［3］．将 真实场 和辅助 场叠加代入J积分，再分离出由真实场和辅助场分别引起的J积分，剩余的真实场和辅助场相互作用项即为相互作用积分I，其表达式为

相互作用积分I和应力强度因子相关联的定义式为

式中：Kaux分别为辅助场中的应力强度因子；E＊为平面应力和平面应变下的弹性模量；μ为剪切模量．

2 二维裂纹应力强度因子数值分析

2．1 基于ANSYS的数值模拟

数值模拟实验的对象为对边受拉应力σ的中心裂纹板，板的高度、宽度和厚度分别为2 H，2W和B，中心裂纹的总长为2a．假设板的材料各向同性、均匀、线弹性，其弹性模量和泊松比分别是E＝2×105MPa和ν＝0．25［4］．数值模拟根据受力不同和试样几何尺寸不同总共进行了6次，表1给出的是不同的实验模型的具体参数．利用模型的对称性，取模型的1／4进行分析，见图1．在裂纹尖端第一排的单元要有奇异性，采用平面单元（PLANE183）进行网格划分，各个模型的网格尺寸相同，方便计算结果的对比．采用Plane strs w／thk，不考虑厚度方向的影响，简化为平面应力模型．模型左边所有的节点X方向以及下面裂纹尖端以右的所有节点的Y方向进行约束，对模型的所有的上边缘进行加载．

表1 数值模拟所选取的模型

图1 中心裂纹板有限元模型

针对上述中心裂纹板模型，利用ANSYS软件进行有限元分析，图2给出了裂尖区域应力云图．然后分别采用位移插值法和相互积分法，计算得出各个模型的应力强度因子．高版本ANSYS提供的CINT命令可基于相互积分法求出SIF，本文考虑裂尖奇异性和应力集中等影响，取第二圈和第三圈单元积分值的算术平均作为最终SIF．

图2 裂尖区域von－Mises应力云图

2．2 与线弹性断裂力学解析解的对比分析

对于无限远受拉中心裂纹板问题，其应力强度因子已经有理论解，参看下式．

表2 应力强度因子的数值分析结果与理论解比较／应力强度因子

由表2可见，数值分析结果和理论解的误差不超过1%，且模型载荷与裂尖应力强度因子的线性关系［5］，验证了 ANSYS求解张开型裂纹SIF的可靠性．其中，基于相互积分的数值计算对网格精度要求低，计算更为简便．

3 三维裂纹应力强度因子数值分析

3．1 应力强度因子的测量实验

本文根据国标GB／T 6398—2000，通过裂纹扩展速率测试实验推算出三维裂纹的应力强度因子，计算公式如下［6］．

式中：ΔK是应力强度因子幅值；α＝2a／W．

所用试件类型为标准的单边缺口三点弯曲试件．厚度方向B＝36mm，试件高度W＝75mm，试件跨距S≈4B＝300mm，预制初始疲劳裂纹长度为9mm．加载装置为JXG－200高频疲劳试验机，材料为Q370qE高强度钢板母材．实验载荷最大值Pmax＝35kN，最小值Pmin＝3．5kN，载荷比为R＝0．1，载荷幅值为ΔP＝31．5kN．

3．2 基于ANSYS的数值模拟

对于上述实验模型，采用ANSYS的SOLID186体单元，沿试件厚度B方向划分3层网格进行有限元计算．试件所受载荷幅值为31．5kN，假定两端分别约束X，Y，Z方向和X，Z方向的位移．采用弹性模量E＝2×105MPa和泊松比ν＝0．3进行计算．分别采用位移插值法和相互积分法，计算得出各个模型的应力强度因子．计算结果和实验结果的比较见表3．

表3 应力强度因子幅值的数值分析结果和实验值比较 MPa·m1／2

表3中列出的是试件表面Z＝75mm处的计算结果．从以上的计算中可以看出，位移差值法求取的ΔK不超过3%，而相互积分法计算误差比前者小一个数量级，更接近与实验测量值．另外，在采用位移插值法求解SIF时，沿着裂尖厚度方向每次只能计算出一个值，需要输入裂纹前缘不同的位置坐标才能计算出不同位置的SIF．但是，相互积分法只需要一次后处理，就可以输出沿着裂纹厚度方向SIF的系列值，更加适用于三维裂纹的数值分析．

4 结束语

文中对二维及三维张开型裂纹应力强度因子数值分析方法研究，并将数值计算结果同线弹性断裂力学解析解和疲劳裂纹扩展实验实测值进行了对比分析，得出如下结论：对于二维裂纹，位移差值法和相互积分法和理论值的误差均不大于1%；对于三维裂纹，采用相互积分法求取的应力强度因子的误差比位移差值法的计算结果小一个数量级．鉴于相互积分法只需要一次后处理，就可以输出沿着裂纹厚度方向应力强度因子的分布，为三维裂纹扩展的数值模拟提供了参考．本文只对静载荷和I型裂纹进行了模拟，今后，还可对动态载荷和复合型裂纹扩展进行更进一步的探索．

［1］崔维成，蔡新刚，冷建兴．船舶结构疲劳强度校核研究现状及我国的进展［J］．船舶力学，1998（4）：63－81．

［2］仇仲翼．应力强度因子手册［M］．北京：科学出版社，1993．

［3］于红军．含复杂界面非均匀材料断裂力学研究［D］．哈尔滨：哈尔滨工业大学，2010．

［4］解 德，钱 勤，李长安．断裂力学中的数值计算方法及工程应用［M］．北京：科学出版社，2009．

［5］张晓敏．断裂力学［M］．北京：清华大学出版社，2012．

［6］国家质量技术监督局．GB／T 6398—2000金属材料疲劳裂纹扩展速率试验方法［S］．北京：中国标准出版社，2000．