断裂参数的数值计算方法分析

刘丹丹，李晓川

(沈阳工业大学建筑工程学院，辽宁 沈阳 110870)

从20世纪50年代开始，固体力学的一个分支断裂力学快速发展了起来，并为现代化的生产提供了理论依据。半个多世纪过去了，断裂力学[1]这一领域的知识理论已经发展的相当成熟。断裂力学有助于人们认识和理解材料中裂纹的产生、扩展和最终破坏的过程，并为结构损伤容限设计提供理论基础。

断裂力学中的3个基本断裂参数是应力强度因子、路径无关积分（J积分）和应变能释放率。应力强度因子的计算主要依赖于裂纹前端的局部应力场，它描述了弹性裂纹尖端应力场的强弱。J积分描述的是由于裂纹的存在所吸收的能量，而能量释放率是描述产生新的裂纹面所需要的能量。这两个基于能量的参数是等价[2]的。随着计算机硬件和软件的迅速发展，用数值方法计算断裂参数就变的切实可行。很多数值方法被尝试应用于断裂参数的计算，如有限差分法、边界元法和发展起来的无网格法，但是由于缺少软件的支持，这些方法应用实例相对缺乏。而有限元法成功的应用于很多工业部门，并发展了针对三个基本断裂参数的数值计算技术。用应力或应变外推法计算应力强度因子，等效积分区域技术计算J积分，全局或者局部的虚拟裂纹扩展技术和虚拟裂纹闭合法计算应变能释放率。

1 外推法

1.1 基于应力的外推法

对于线弹性材料，应力在裂纹尖端处是趋于无穷大的，也就是应力奇异性。所以网格尺寸越小，裂纹前端单元内的积分点的应力值就会越大。从这可以知道，应力值对于网格的尺寸是不收敛的。于是引入了应力强度因子的概念来克服由于应力奇异带来的数学上的困扰，用来描述裂纹尖端附近应力的强弱。外推法是借助于有限元的方法来计算应力强度因子，在有限元分析中，裂尖前面的单元里积分点上的应力值以及其对应的积分点坐标值是很容易就可以直接读出的。虽然无法用数值方法直接计算裂尖处的应力强度因子，可是裂尖前端那些非奇异的应力值却是可知的。对于每一个ri（距离裂纹尖端的极半径）大于零都有一个非奇的应力值σyi和对应的KIi：KIi=σyi，然后可以构造数据对（ri,KI）i，用最小二乘法拟合数据点。在这里假定了两者之间线性关系为：

上式中的N 是外推法所采用的数据对象，其截距的物理意义就是所要计算的应力强度因子。

可是这个方法存在几个问题。第一，没有明确规定采取多少对数据点进行外推计算得到最后的精度高。一般情况下是采用许多点直至KI对N 收敛于某个值才行。第二，靠近裂尖处的数据反而是引起误差的主要来源，而相反的距离裂尖越远越满足线性的假设，这样使得外推法所采用的点往往不包括贴近裂尖的点，说明线性外推对裂尖来说并不很精确，而在断裂分析中又恰好最应该多关注裂纹尖端，因为应力强度因子的作用就是描述裂尖应力信息，计算一个用来描述裂尖应力场的量所采用的方法却不能依靠裂尖信息，而必须用远场数据来计算，这在逻辑上存在悖论。

1.2 基于应变的外推法

在很多的商业软件中，位移是求解的基本变量，而应力是通过应变和位移联系起来是次要的变量，所以位移的外推法的计算精度要高。对于确定的距离裂纹尖端r 处，裂纹后端垂直位移v 的数据可以在有限元分析中直接读取。和基于应力的外推法一样，也可以构造数据对（ri,KIi）：

然后利用最小二乘法来拟合数据点，得到相应的应力强度因子数值。同样的，基于应力的外推法存在的问题也存在于此方法中。如果想解决问题，使用多项式外推应该是最直接的方法，可是并不是所有的情况都存在确定的多项式。因此，引入奇异单元和叠合单元[3]理论上在裂尖处采用这些单元，在其他区域采用常规的单元，可以极大地减少网格的数量同时保持计算精度。然而采用这些单元在使用时非常的不方便，尤其是对于裂纹扩展的问题，而且适用范围也很有限，无法满足实际情况多样性和复杂性的需求。

2 等效积分区域法

J积分是由Rice[4]提出来的一个处理非线性断裂问题的断裂参数，这个参数的引入是基于能量守恒的概念，因此对裂纹尖端应力奇异性的依赖程度比较弱，不需要对裂纹尖端处的单元进行特殊的处理，而是从全局的角度分析处理断裂问题。J积分的数学表达式为：

上式中，ui是位移矢量的分量；ds是沿着积分路径的微小增量；w 是应变能密度因子，其定义是：

上式中σij和εij分别是应力张量和应变张量。

由于Rice证明了J积分的数值不随着积分回路的改变而改变，因此它也叫做路径的无关积分。在这里可以看到J积分的表达式中，只有应力、应变和位移三场没有材料属性。但是在计算应力的时候是需要用到材料各向同性或者异性的本构关系。所以，这个表达式只是使用于一般的热弹塑性和各向异性材料的问题。实际上这个表达式并不适用于数值计算，为此Shih和Raju[5]等人提出了等效积分区域法来达到对J积分进行数值计算。基于用裂尖附近的一个有限区域来代替积分回路进行J积分计算的思想，表达式可以转化为：

上式中δij是Kronic函数，q 是辅助函数，它只是一个数学处理，目的是为了让积分表达式更利于采用数值计算方法。J积分是远场回路积分，它具有较好的精度，也不需要对裂尖处的单元进行特殊处理，不过由于没有分离断裂模式，对于混合断裂型问题而言，还需要额外处理才能进行完整的分析。就算对于平面问题，它的计算过程都显得相当的复杂。所以尽管它能够用于裂纹扩展问题但是实际上操作很不方便，而且在有限元的分析中实用性很低。

3 虚拟裂纹扩展法

3.1 全局的虚拟裂纹扩展法

Irwin[6]提出了应变能释放率的概念。考虑有一个厚度为B的二维裂纹体，其裂纹长度是a，那么应变能释放率就是产生面积为ΔA 的新裂纹面所需要的能量，于是：

上式中∏=U-W 是裂纹体的势能，U 是裂纹体的应变能，W是外力功，Δa 是裂纹的微小扩展量。

从上式可以看出，应变能释放率的计算要求Δa 趋近于零，很明显对于有限元分析这种数值方法来说，无法达到这个极限。所以，采用两步分析过程的虚拟裂纹扩展法。在第一步的分析中，通过有限元分析获得裂纹长度是a 的裂纹体的势能（∏1=U1-W1）；第二步分析中，通过有限元分析获得裂纹长度是a+Δa 的裂纹体的势能（∏2=U2-W2）。若和网格尺寸相关的Δa足够小，那么就可以很好的将应变能释放率近似为：

全域虚拟裂纹扩展法有一个局限性就是只能够得到总的应变能释放率，不能分离断裂模式。

3.2 局部的虚拟裂纹扩展法

基于将裂纹闭合一个微小的扩展增量所需要的功和势能的改变是等效的这一思想，可以用裂纹的闭合积分来计算裂尖的能量释放率，则局部的形式为：

其中应变能分量为：

假设裂纹的方向是沿着X 轴方向，σyy是张开断裂模式中沿着闭合裂纹面上的法向应力，txy是滑移模式中沿着闭合裂纹面上的切向应力。Δu 和Δv 是闭合裂纹张开后裂尖后面闭合裂纹面上的位移分量。上两式可以用两个步骤的分析过程计算出来。如果网格尺寸相关的Δa 很小，就可以将应变能释放分量近似为：

如果直接用上两式计算应变能释放率分量则需要节点上的应力值，因为公式中是沿着闭合裂纹线对应力进行数值积分而闭合裂纹线一般位于有限元单元的边上。在有限元分析中，当把应力外推到单元的节点上，或者当单元接近裂纹尖端时，所得到的应力时非常不精确的。为了避免使用不精确的应力值，采用节点力来代替应力的积分。用节点力和节点位移计算出的应变能释放率为：

用节点力和节点位移计算应变能释放率有很多的优点是：①公式中只有节点力和节点位移这两个基本变量，任何有限元计算的商业软件都可以直接输出；②公式不需要对应力进行数值积分，这样减少了计算量；③有限元网格尺寸的影响小，为了保证近似精度，需要采用合理的网格密度，但是用比较粗糙的网格也能获得令人满意的结果；④无论是线性的还是非线性的材料都可以采用这个方法。

而此方法的缺点是分析过程必须是2步，且裂纹长度在这两步中是不同的，这样增加了网格准备工作量，尤其对于三维问题时，同时，对于裂纹扩展问题，由于需要用上一步的计算结果来不断的准备新的网格，使得研究起来相当的不方便。

4 虚拟裂纹闭合法

为了避免虚拟裂纹扩展法的不便之处。Rybicki和Kanninen提出[7]适用于二维平面问题的一步分析法即修正的裂纹闭合积分，也就是后来的虚拟裂纹闭合法。

虚拟裂纹闭合法的基本假设是虚拟裂尖和实际裂尖后面张开位移近似相同。于是，公式转化为：

其中，Fy1和Fx1为裂纹尖端节点Y 和X 方向的节点力；Δu3,4和Δv3,4为紧挨裂纹尖端后面节点在X 和Y 方向的相对位移；B 为裂纹体的厚度。Δa 为裂纹尖端前面的虚拟裂纹扩展量。

虚拟裂纹闭合法在保留了虚拟裂纹扩展法的所有优点的同时，克服了它的所有缺点，只要求一步的有限元分析。这使得此方法有很大的吸引力，尤其是那些用有限元计算商业软件对断裂问题进行分析研究的工程师们。

5 结 论

由上述分析可知，在结合有限元分析断裂问题的数值计算方法中，虚拟裂纹闭合法非常简便高效。

[1]丁遂栋.断裂力学[M].北京:机械工业出版社,1991.

[2]沈成康.断裂力学[M].上海：同济大学出版社,1996.

[3]Blandford,G.E..Two-dimensional stress intensity factor computations using the boundary element method.International Journal for Numerical Methodsin Engineering,1981.

[4]Rice,J.R..APath Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks,Journal of Applied Mechanics,1968.

[5]Moura,B.,and Shih,C.F..A Treatment of Crack Tip Contour Integrals.International Journal of Fracture,1987.

[6]Irwin GR.Crack-extension forcefor a part-through crack in aplate.J Appl Mech,1962.

[7]Rybicki EF,Kanninen MF.A finite element calculation of stress intensity factors by a modified crack closure integral.Engineering Fracture Mechanics,1977(9).