摘 要：基于MSC Marc软件，本文通过编程实现了参数化、模块化建模与自动化分析；采用模块化有限元建模技术和接触技术仿真半椭圆表面裂纹(前缘)，采用在裂纹前缘形成辐射状奇异单元网格的建模方法和接触技术来施加边界条件的方法，实现裂纹前缘奇异应力场和裂纹模型远场应力的模拟。在此基础上，运用裂纹虚拟扩展技术，计算裂纹前缘的J积分, 通过线弹性有限元理论，换算出裂纹前缘的K因子分布，对裂纹前缘离散节点实行正交扩展，通过拟合裂纹前缘，实现了裂纹连续扩展仿真。仿真结果与试验结果基本吻合，说明本文所采用的方法具有很好的实用性。
关键词：接触技术；奇异单元；裂纹前缘；半椭圆表面裂纹；MSC Marc
国内外对于裂纹的研究，多数集中在理论推导和试验上[1]，只能通过理论或经验公式指导工程结构的可靠性设计。许多研究者[2-6]介绍了表面裂纹的应力强度因子的各种理论计算方法和分布特征来表征裂纹的扩展。随着计算机技术和有限元技术的快速发展，使裂纹及其扩展仿真成为可能，为含裂纹结构可靠性研究创造了条件[7-11]。
裂纹虚拟扩展技术和J 积分能量定义是有限元法计算裂纹前缘应力强度因子分布的理论基础；而有限元法计算裂纹前缘K 因子分布规律是实现裂纹连续扩展的前提条件，从计算裂纹前缘的J 积分（能量释放率）分布着手，在线弹性情况下可以换算裂纹前缘的应力强度因子分布。本文通过软件实现裂纹前缘的不同步扩展。事先不给出裂纹前缘曲线的形状假设，而将裂纹前缘离散成系列点，计算裂纹前缘应力强度因子的分布， 即局部（各个离散点）应力强度因子，并对裂纹前缘的各个离散点单独实行正交扩展，拟合得到的代表新裂纹前缘的系列离散点，来获得裂纹前缘曲线的形状表达，实现裂纹前缘的不同步扩展(裂纹前缘各部分扩展速度不同)。这样可获得更加精确的典型裂纹的扩展特性和寿命。
1 裂纹连续扩展模型的建立
为了真实地仿真裂纹尖端应力场，本文采用MSC MARC 仿真裂纹、施加载荷与边界条件，采用参数化、模块化建模技术建立裂纹连续扩展模型。图1 是含有半椭圆表面裂纹体。
图1 含有半椭圆表面裂纹块体
1.1 裂纹的参数化、模块化建模
在MSC MARC 软件中，含有半椭圆表面裂纹的四分之一模型被划分为7 个模块，各模块分别表示为：裂纹载体tets，裂纹体包括三个子模块：crack1 模块、crack2 模块、fill 模块， z方向对称面cracksym， x 方向对称面symmx，载荷施加面moving。它们之间的位置关系如图2 所示。把裂纹体镶嵌入裂纹载体中，从而可单独对裂纹体划分特定的映射网格，裂纹载体划分为四面体网格。
图2 各个模块之间的位置关系
1.2 裂纹前缘及其节点的定义
通过接触技术定义，z 方向对称面cracksym 只对裂纹韧带部分而不对裂纹自由表面运用对称边界条件，即cracksym 仅与裂纹的tets 和crack1 模块粘连在一起，而与crack2 和fill 模块没有连接关系。在模块crack1 和crack2 交界的曲线上，两模块的节点是重合的，但是crack1的节点属于裂纹韧带，与对称面cracksym 粘连（Glue）在一起；crack2 的节点属于裂纹自由表面与cracksym 没有连接关系；这样在模块crack1 和模块crack2 交界曲线处形成了三维的裂纹前缘。相应的裂纹前缘曲线的离散节点用来求解局部J 积分。
2 半椭圆裂纹模型的前缘离散正交扩展与拟合
图3-1、图3-2、图3-3、图3-4 给出了半椭圆裂纹前缘扩展拟合过程的图形表示。图3-1表示按弧长均匀离散裂纹前缘，使用离散节点代替裂纹前缘；图3-2 表示前缘离散节点正交扩展得到新前缘的离散点，新获得的点不能直接用来表示裂纹前缘，要获得裂纹前缘的数学描述，需要进行数据拟合；图3-3 表示椭圆拟合新的离散节点，获得椭圆前缘参数；图3-4 表示按弧长均匀离散新拟合的椭圆裂纹前缘；重复图3－2、图3－3、图3－4 所示的操作就实现了裂纹的连续扩展。
2.1 正交扩展计算裂纹前缘节点坐标
裂纹前缘的离散节点数在20～30 个节点之间。为了清楚说明问题，图形中裂纹前缘的节点取为5 个。设椭圆裂纹长半轴为A ，与x 轴平行，短半轴为B ，与y 轴平行； 坐标原点、椭圆的中心为O，如图4 所示。
图4 裂纹前缘节点正交扩展
裂纹前缘系列离散点的坐标序列记为(x1,y1),...(xi,yi)...(xn,yn)。在点(xi,yi)作椭圆的切线 L1i ，则点(xi,yi)的导数就是切线 L1i 的斜率，记为 k1i ，对椭圆方程两边求x 的导数，求得：
过(xi,yi)点与切线 L1i 垂直的直线 L2i 的斜率记为k2i=-1/k1i；在点(xi,yi)，裂纹的扩展量为△ai ，需要求解在直线 L2i 上，距点(xi,yi) 长△ai 的点的坐标，方法如下：
在程序中有效点的判断方法为：将(xi1,yi1)，(xi2,yi2)代入椭圆方程，如果大于1就是点在椭圆外，即点(xi1,yi1)；如果小于1 就是点在椭圆内，即点(xi1,yi2)。裂纹是向外扩展，点(xi1,yi1)为有效点。
2.2 拟合裂纹前缘及其步长控制
对裂纹前缘系列离散点进行椭圆拟合，变量代换可将椭圆拟合转化为线性拟合：
由最小二乘法直线拟合可以算出a和b ，这样可以解出椭圆的长短轴A 和B 。
图5 椭圆拟合裂纹前缘的误差分析
如图5 所示，坐标系中心为O，椭圆中心也为O(0,0)，裂纹前缘正交扩展后节点的坐标序列记为(x1,y1),...(xi,yi),...(xn,yn)，过裂纹正交扩展获得新节点与椭圆原点作直线，新节点到椭圆原点的距离为
拟合椭圆的交点到椭圆原点的距离OP。P 点的坐标是直线与椭圆的交点。定义裂纹前缘节点的相对偏差：
通过控制（8）式表示的误差不超过1％来确定裂纹扩展的步长，同时保证了拟合裂纹前缘的精度。
3 半椭圆表面裂纹连续扩展仿真试验
3.1 仿真试验参数
仿真试验采用ZG25 材料作裂纹块体，其扩展速率为：da/dN =1.44×1010（△K）2.79；门槛值△Kth=6.3246 Mpa·m1/2。取含有半椭圆表面裂纹块体的四分之一为有限元模型，厚度t 为10mm，宽度B 为15mm，长度L 为15mm，参见图1；a 为裂纹深度（短半轴），c 为裂纹长度（长半轴），a0 和c0 为初始值。裂纹前缘划分为22 份，即离散为23 个节点；每个节点设置4 个J 积分半径，对应可以计算4 个J 积分值。按照脉动载荷计算K 因子范围。
仿真试验的载荷参数如表1 所示。载荷条件为拉伸、弯曲、拉弯三种，确定载荷的大小是使裂纹前缘有4 个离散点局部K 因子大于门槛值。采用相同初始半椭圆表面裂纹模型：a0 = 2mm，a0/c0= 0.8；进行等步长裂纹扩展仿真。表1 仿真试验载荷参数
 3.2 仿真试验结果
仿真试验结果如图6 所示。
由图6-1、图6-2、图6-3 中裂纹前缘变化规律可见，在裂纹扩展相同长度（长轴c方向）时，拉伸载荷作用下，裂纹沿深度方向（短轴a 方向）扩展较快，而弯曲载荷作用下则较慢，拉弯载荷作用下介于两者之间；当裂纹深度超过一定尺寸后，三种载荷条件下裂纹沿深度方向的扩展均逐渐放慢，而沿长度方向逐渐加快。
图6-4 的a/t－a/c 关系也反映出在拉伸载荷作用下，裂纹沿深度方向扩展较快，在弯曲载荷作用下裂纹沿长度方向扩展较快，拉弯载荷作用下介于两者之间。
3.3 仿真结果与试样实验结果比较
本文进行了拉伸载荷下的裂纹扩展实验，将仿真结果与实验结果进行比较，得出图7 所示结果。文献[12]是对含有不同a0/c0 的初始半椭圆表面裂纹的四点弯曲试样进行裂纹扩展直至断裂的实验，裂纹前缘的形状采用AC 电位法监测，采用降载勾线法在试样断面上留下海滩状条纹。从图7 中可以看出，本文半椭圆表面裂纹连续扩展模型的仿真结果与试验结果基本吻合，这反映了本文的半椭圆表面裂纹连续扩展模型是有效可行的。
图7 本文仿真结果与文献[12]实验结果对照
4 结论
本文建立了3D裂纹连续扩展模型，通过编制接口程序，得到不同步扩展的裂纹前缘。实现了裂纹连续扩展模型的总体框架。并进行了半椭圆连续扩展仿真与实验验证，结论如下：
1.裂纹前缘的演变轨迹呈半椭圆形，这说明本文采用半椭圆曲线拟合裂纹前缘是合理的。
2.在裂纹扩展深度超过一定尺寸后，在深度方向的裂纹扩展逐渐放慢，在宽度方向的扩展加快，这与本文半椭圆裂纹扩展仿真结果一致。
3.载荷作用下，本文半椭圆表面裂纹连续扩展模型的仿真结果与裂纹扩展特性试验结果基本吻合，这说明了本文建立的半椭圆表面裂纹连续扩展模型是有效可行的。
参考文献
[1] 高庆, 工程断裂力学[M]．重庆：重庆大学出版社, l986
[2] Smith.RA, LIN Xiao-bin. The Calculation of Stress Intensity Factor using 3D Finite Element Method[J]. CHINA MECHANICAL ENGINEERING. 1998, 9(11), 39-42
[3] Newman Jr J C, Raju I S. An empirical stress intensity factor equation for the surface crack[J]. Engng Fracture Mech, 1981, 15; 185-1 92
[4] 陈道礼. 基于有限元分析的表面裂纹应力强度因子的拟合估算[J]. 机械设计与制造, 2003,3, 6-7
[5] 钟明, 张永元. 用时域边界元法分析半圆表面裂纹的动态应力强度因子[J]. 应用数学和力学, 2001, 22(11), 1211-1216
[6] 曹宗杰, 闻邦椿, 王志超. 一种计算三维裂纹应力强度因子的新方法[J]. 力学季刊,2001,22(4),401-407.
[7] 吴志学. 表面裂纹疲劳扩展的数值模拟[J]. 应用力学学报, 2006, 23(4), 563-567
[8] Hellen T K. On the method of virtual crack extension[J]. Int J Numer Meth Eng. 1975,9(1) 187—207
[9] C.G.Hwang, P.A.Wawrzynek, A.R.Ingraffea. On the virtual crack extension method for calculating the derivatives of energy release rates for a 3D planar crack of arbitrary shape under model-I loading. Engineering Fracture Mechanics[J]. 68(2001), 925-947
[10] J.Z. Yang, J.F.Chen, G..D. Holt. Efficient evaluation of stress intensity factor using virtual crack extension technique. Computers & Structures[J]. 79(2001),2705-2715
[11] A.Bacila, X.Decoopman, et al. Computer simulation of fatigue crack propagation under random loading conditions[J]. International Journal of Fatigue, 29 (2007) 1772–1780
[12] McFadyen NB, Bell R, Vosikovsky O. Fatigue crack growth of semi-elliptical surface cracks[J].Int J Fatigue, 1990, 12
|
