应力强度因子的求解方法的综述

摘要：应力强度因子是结构断裂分析中的重要物理量，计算应力强度因子的方法主要有数学分析法、有限元法、边界配置法以及光弹性法。本文分别介绍了上述几种方法求解的原理和过程，并概述了近几年来求解应力强度因子的新方法，广义参数有限元法，利用G*积分理论求解，单元初始应力法，区间分析方法，扩展有限元法，蒙特卡罗方法，样条虚边界元法，无网格—直接位移法，半解析有限元法等。

关键词：断裂力学；应力强度因子；断裂损伤；

Solution Methods for Stress Intensity Factor of Fracture Mechanics

Shuanglin LU

(HUANGSHI Power Survey&Design Ltd.)

Abstract: The solution methods for stress intensity factor of fracture mechanics was reviewed, which include mathematical analysis method, finite element method, boundary collocation method and photo elastic method. The principles and processes of those methods were introduced, and the characteristics of each method were also simply analyzed in this paper.

Key words: fracture mechanics; stress intensity factors

中图分类号：TU318 文献标识码：A

0 引言

断裂力学的基础理论最初起源于1920年Griffith的研究工作[1]。Griffith在研究玻璃、陶瓷等脆性材料的断裂现象时，认为裂纹的存在及传播是造成断裂的原因。裂纹的扩展过程,从能量的观点来看,存在着两种完全对抗的因素：一种是阻止裂纹扩展的因素，另一种是推动裂纹扩展的因素。Griffith由此建立了材料的脆性断裂判据[1]：

(1)

在(1)式中： —断裂应力；E—材料的弹性模量； —材料的表面能；a—裂纹长度的一半。

Griffith判据并不能完全成功地应用于金属断裂问题。1949年, Orowan考虑到裂纹释放的应变能不仅转化成表面能，也同时转化成使裂纹顶附近材料发生塑性变形所需要的功。因此，Orowan对Griffith判据进行修正并得到了具有塑性变形的金属材料的断裂判据[1]：

(2)

在(2)式中： —断裂应力；E—材料的弹性模量； 为塑性功；a—裂纹长度的一半。

1975年，Irwin认为裂纹是脆性断裂破坏的要害，而裂纹顶端区域的应力场又是其中的核心。从(1)、(2)可以看出： 是一个常数，也就是说与载荷条件、式样尺寸、裂纹大小毫不相干，是只由材料的固有性质决定的不变值。当 大于这个值时裂纹就快速扩展，因而，这个常数才真正代表了材料对断裂的抵抗能力。于是，Irwin对应提出了一个崭新的物理量—应力强度因子。

由裂纹尖端的应力应变的表达式[2]可以看出：裂纹尖端附近各点的应力、应变和位移均由应力强度因子K唯一确定，因此，如何计算K值是断裂力学中的一个重要内容。目前，对于无限体中的简单裂纹和有限边界的贯穿裂纹，确定K值的主要方法有：数学分析、数值计算、试验标定以及光弹性法等。

1数学分析法

1.1复变函数法

对于平面弹性问题，利用复变函数能够很方便的求得裂纹尖端应力应变场。在文献[2]中详细给出了针对 型裂纹，利用威斯特葛尔德(Westergaard)应力函数求解应力分量的过程，最后得到各应力分量的表达式为：

(3)

根据(3)式可以由胡克定律得到应变分量，然后再根据应变与位移之间的关系式可以得到位移分量的表达式。由上所述可以看出，只要知道了ZI函数的表达式，应力分量、应变分量和位移分量都可以求出来了。

因此，用复变函数法求解应力强度因子的思想就是，针对不用的裂纹情况构造出满足相应边界条件的复变解析函数，并由此复变函数求得裂纹尖端的应力应变场，最后由应力强度因子的表达式求得K值。复变函数法在弹性平面问题的应用中比较方便，但对于弹塑性或三维空间问题，该方法就不再实用，其主要原因是构造满足边界条件的复变函数很困难。文献[3]和文献[4]中给出了利用复变函数法求解正交各向异性含内部裂纹板、带单裂纹无限平板中作用有集中力和力矩以及带单裂纹无限弹性体作用有纵向集中力等情况下应力强度因子的计算方法。

1.2积分变换法

弹性理论已经证明，常体力下弹性平面问题存在应力函数，称为Airy应力函数，为双调和函数[5]。对于平面问题，可用Laplace Transform和Fourier Transform来解答应力场强度因子。鉴于求解方程为 4Ψ=0（Ψ为Airy应力函数）很困难，故可考虑Fourier Transform来解断裂力学问题。首先对Ψ取Fourier变换，记为 ，即：

(4)

于是， 应满足方程：

(5)

用降阶法可以求出方程(5)的通解为：

(6)

由(6)式结果来求解应力分量如下：

(7)

其相应的位移场为：

(8)

经过反演分析即可得出Ψ以及σ，μ等全部场量。

如用Fourier变换仍求解椭圆形裂缝问题得KI，则由：

(9)

一旦两个材料参数m、s确定，则KⅠ、KⅡ的数值可以根据下列公式十分容易地求得：

(10)

在式(10)中：σ为材料的抗压强度；l为裂纹长度。

2边界配置法

由弹性力学可知，二维弹性力学问题的应力函数 为双调和函数，即满足微分方程式： 。当裂纹表面满足边界条件 ， ， 时，有Williams无穷级数的应力函数[6,7]：

(11)

其中：

(12)

在(12)式中： 为偶函数部分，相当于Ⅰ型裂纹里对称加载； 为奇函数部分，相当于Ⅱ型裂纹里反对称加载。

应用复应力强度因子公式：

(13)

注意到(12)式中的Cj=-Cj/2=-Cn以及Dj=-Dj/2=-Dn，因此有，Cj/2 -C1和Dj/2 -D1

故有：

(14)

即：

(15)

因此，要计算应力强度因子KⅠ、KⅡ，则先要求解(12)式。为此，需要由边界条件建立含有Ci、Di的线性方程组，求解此方程组以确定系数C1、D1。

由弹性力学可知，弹性力学问题的解必须满足平衡条件和边界条件。这里，在边界上取2m个配置点，对于每一个配置点i可以提出两个边界条件：

(16)

在(16)式中： ， 分别为含裂纹体的应力函数及其法向偏导数； ， 分别为非裂纹体的应力函数及其法向偏导数。

因此，对于2m个配置点便可以建立4m个类似的边界条件，由4m个方程式组成线性方程组。解此线性方程组即可求得4m个未知量 的值。采用边界配置法就是将(11)式或(12)式截断，然后由边界上的2m个配置点处4m个边界条件去确定其中的4m个待定常数Cj、Dj，把问题归结于求解4m个线性方程组，用计算机及程序计算很方便。

3有限元法

随着有限元法的发展，有限元在断裂力学中的应用越来越普及。近些年，计算机技术得到迅猛发展，很多大型通用软件，如ANSYS、ADINA以及MSC/Nastran等都具有计算各算各种断裂参数的功能，因此利用有限元计算断裂力学中的应力强度因子也得到广泛的应用。

构件中的裂纹可以抽象为二维或三维模型，如图1所示。求解断裂力学问题的步骤包括先进行弹性或弹塑性静力分析，然后用特殊的后处理命令或宏命令计算所需的断裂参数。在有限元中主要采用1/4法计算应力强度因子。根据县弹性断裂理论，裂纹尖端的位移场可以表示为[7]：

(17)

在(17)式中：u、v和w为对应于裂纹尖端局部坐标的位移；r和θ是计算点在局部柱坐标的坐标值；G是剪切模量,v是泊松比；对于平面应力 ，而对于平面应变 。 、 和 分别为Ⅰ型、Ⅱ型和Ⅲ型裂纹的应力强度因子；0(r)是高阶无穷小量。根据公式(17),如果裂纹表面(θ=±180°)某一点垂直于裂纹平面的位移已知，可以导出对称裂纹的应力强度因子计算公式：

图1裂纹的二维和三维模型

(18)

对于非对称裂纹体,其应力强度因子的计算公式为：

(19)

在(19)式中：Δu、Δv和Δw分别为两个裂纹面之间的相对位移。

由于裂纹尖端的应力和应变是奇异的,因此在进行有限单元建模或单元网格划分时,必须先在裂纹尖端位置定义应变奇异点,而且围绕裂纹定点的有限单元是二项式的奇异单元,它是把单元边上的中间点放到1/4边处。图2所示为ANSYS的2-D和3-D模型中所采用的奇异单元。

图2裂纹尖端的奇异单元

应用有限元方法计算裂纹体的应力强度因子,关键是要建立一个能够反映裂纹体特征的共线(共面)的裂纹几何模型,并确定裂纹尖端的局部坐标。在划分裂纹尖端附近的几何体时,必须选用具有奇异特征的单元。在完成静力学计算后,才能计算裂纹尖端的应力强度因子。文献[8]和文献[9]中的计算结果表明，应用有限元分析软件计算出的应力强度因子与断裂力学求得的应力强度因子非常相近，因此，利用有限元计算材料的断裂强度因子是可行的。

4光弹性法

由于光弹性法可以确定光弹性模型在裂纹尖端附近的应力变化规律，因此提供了用实验方法确定裂纹尖端应力强度因子K的基础[11]。利用透光材料制成含裂纹的试件，用激光光源照射，由于实时全息干涉原理，在照片上可以看到一组以裂纹尖端为中心的明暗交替的条纹。可以证明：条纹中光的强度和试件的主应力 、 间的关系如下：

(20)

在(20)式中： 是材料的应力—光学常数； 是光的波长； 是光波振幅； 是光的强度。

因为出现暗条纹的条件是 ，即：

也就是：

(21)

引进常量m，它与条纹序数N的关系为： ，因此：

(22)

对于张开型裂纹，在裂纹延长线上( )，由裂纹尖端应力分量的表达式可看出 xy=0，因为在裂纹上的剪应力为0，所以σx和σy就是主应力σ1和σ2。因此，由裂纹尖端应力分量的表达式可得：

(23)

由于(23)式是在双向应力σ作用下导出来的，为了得到单项拉伸下的应力场公式，可在x方向叠加一套应力 ， ， ，但这并不改变裂纹尖端的奇异性和KⅠ值，这套应力在裂纹内产生一个均匀的应力场，故x方向的合力为：

则单向拉伸时x轴上的应力为：

(24)

(25)

将(25)式代入到(22)之中得：

(26)

在远离裂纹处，只有在y方向的均匀拉应力 ，这时σ1+σ2=σy=σ,该处的m用 表示，代入到(22)式得：

(27)

联立(26)式和(27)式得：

(28)

由于一般KⅠ的表达式为：

(29)

将(29)式代入(28)式得：

(30)

由 可得：

(31)

联立(28)式和(31)式，得：

(32)

其中，N为裂纹线上距裂纹顶端为r的干涉条纹序数，N*为远离裂纹其应力等于均匀拉应力 处的条纹序数。按(32)式可以 为纵坐标， 为横坐标的直角坐标系中将实验结果画出，它是一条直线，其斜率就是Y。将Y值代入到式即可得到KⅠ值。

实时全息条纹法只能得到二维问题的裂纹尖端数值解，对于三维裂纹问题则不可行。

5 近几年求解应力强度因子的新进展

近几年来求解应力强度因子的新方法主要有广义参数有限元法[12]、利用G*积分理论求解[13]、单元初始应力法[14]，区间分析方法[15]、扩展有限元法[16]、蒙特卡罗方法[17]、样条虚边界元法[18]、无网格—直接位移法[19]、半解析有限元法[20]等。

广义参数有限元法建立了裂尖处应力强度因子计算的W单元。利用修正的Williams级数建立裂尖附近奇异域的整移场，使得计算模型中含有与应力强度因子直接相关的参数，便于直接计算应力强度因子，从而避免了奇异单元需要外推计算且人为选择直线和计算点带来的计算误差和种种不便，不仅便于应用，而且计算精度较高。

直接计算应力强度因子的扩展有限元法以常规有限元法为基础,利用单位分解法思想,通过在近似位移表达式中增加能够反映裂纹面的不连续函数及反映裂尖局部特性的裂尖渐进位移场函数,间接体现裂纹面的存在,从而无需使裂纹面与有限元网格一致,无需在裂尖布置高密度网格,也不需要后处理就可以直接计算出应力强度因子,并且大大简化了前后处理工作。

应用区间分析方法对具有不确定参数的应力强度因子进行估计。该方法以区间数学为基础,将不确定参数描述为区间变量;再利用Taylor级数展开通过区间运算得到应力强度因子的区间范围,从而为工程设计提供可信的数据。区间分析方法优于传统的概率分析方法的是:它不需要预先知道关于不确定参数大量的统计数据信息,并且具有计算方法简便、实用和精度高的特点。

运用蒙特卡罗方法求解线性代数方程组时，仅仅是求得解得一个分量Xi，而与其它分量无关。这条性质跟我们边界配置法所需要的相一致。因为，边界配置法最后就是得到一个线性方程组，而且为了求得应力强度因子，就只是需要求出其中一个解的分量。运用蒙特卡罗方法求解由边界配置法得出的线性方程组，即求出的应力强度因子。采用差分方程求解的蒙特卡罗方法采用的游动网格为规则网格，对于不规则的几何边界问题，由于边界的近似处理，计算精度受到了很大的影响。采用不规则游动网格的蒙特卡罗方法，改善了蒙特卡罗方法求解复杂边界问题的精度，大大的扩宽了蒙特卡罗方法的应用。

采用基于Kelvin基本解的样条虚边界元法，结合位移外推法，给出了断裂问题应力强度因子的求解方法通过对两个典型断裂问题的分析，对边界子段与虚边界元的划分、小单元的采用以及拟合点位置的确定等关键问题展开了讨论，获得了相关计算参数的选取规律，为该法在断裂问题的进一步应用打下良好的基础。

计算裂纹尖端应力强度因子的无网格法一般均采用J积分方法,但由于该方法为间接求解,降低了求解精度与求解效率。文中采用无网格—伽辽金方法,选取带有扩展基的奇异基函数,以精确计算裂纹尖端位移场,并借鉴有限元法中计算应力强度因子的直接位移法,提出一种计算含裂结构裂纹尖端应力强度因子的新方法,即无网格—直接位移法。数值计算结果表明,该方法具有简捷、高效的特点,可以准确计算裂纹尖端应力强度因子。

半解析有限元法，从弹性力学哈密顿理论出发，在裂纹尖端构造了一个解析的超级单元，该单元能够准确描述平面裂纹尖端场，将该超级单元与普通有限单元相结合可求解任意几何形状和载荷的平面裂纹应力强度因子计算问题。该方法将解析法与有限元法相结合，各取所长，发挥各自的优势。

6结语

(1) 数学分析法中的复变函数法一般只能解决弹性平面问题，且比较简单方便，但对于较为复杂的三维问题则无法用复变函数去求解。此外，对于积分变换法，同样也只能较为方便的去求解弹性平面中的裂纹应力强度因子，故数学分析法有着其本质的局限性。

(2) 边界配置法其本质是根据边界上的配置点建立线性方程组，以满足平衡条件和边界条件，最后由计算机程序求解线性方程组，但其结果的精度与边界配置点的选取相关，且需要有编程的基础才能较为方便快捷的计算出应力强度因子。

(3) 有限元法是目前用的较为普遍的一种计算应力强度因子的方法，该方法能够解决平面和三维的各种裂纹情况下的应力强度因子，并且还能复合裂纹的动态应力强度因子，且计算结果的精度也比较高。因此，该方法在今后将会得到越来越广泛地应用，但利用有限元计算强度因子的建模过程很繁琐，它需要将裂纹尖端附近的网格划分的十分密集，工作量很大。

(4) 光弹性法属于实验分析的一种方法，因此，用光弹性法首先需要有实验仪器上的支持。此外，光弹性法中的实时全息条纹法只能得到二维问题的裂纹尖端应力强度因子的数值解，无法解决三维裂纹问题。现在通过三维光弹性冻结技术通过切片后在偏振光场中观察到裂纹尖端附近应力场的最大剪应力条纹，可以估算出KⅠ。

(5) 对于应力强度因子的计算还有其它的一些计算方法，如柔度标定法、权函数法以及无网格直接位移法等，针对不同的问题应该选择合适、简便的方法去求解，因为每种计算方法都有着其各自的局限性或缺点。总之，随着断裂理论的不断发展，应力强度因子的求解方法将越来越成熟，计算三维空间复杂裂纹应力强度因也得到一定的发展，也是今后发展的重要方向。

参考文献：

[1] 雷振德. 解读应力强度因子[J]. 武钢职工大学学报,2000,11(1):68~76.

[2] 尹双增. 断裂、损伤理论及应用[M]. 清华大学出版社,1992.

[3] 崔德渝, 张 行. 求解正交各向异性含内部裂纹板应力强度因子的复变函数—分区广义变分解法[J]. 航空学报,1991,12(1):87~94.

[4] 陈宜周. 利用特定复变函数积分计算裂纹端的应力强度因子[J]. 江苏工学院学报,1986,7(3):1~6.

[5] 李炳锋, 董世方. 应力强度因子的积分变换与数值求解法的比较[J]. 湖南理工学院学报(自然科学版),2005,18(4):26~28.

[6] 熊先仁, 袁去惑, 杨菊梅等. 应力强度因子的边界配置法程序及工程应用[J]. 江西科学,1996,14(2): 106~113.

[7] 徐德福, 罗景文. 用边界配置法计算正交各向异性材料单边裂纹试件的应力强度因子[J]. 固体力学学报,1982,(2): 283~291.

[8] 陈家权, 沈炜良, 徐家园等. 应力强度因子的有限元计算[J]. 装备制造技术,2003,(4): 6~9.

[9] 赵海涛, 战玉宝, 杨永腾. 基于ANSYS的应力强度因子计算[J]. 煤矿机械,2007,28(2): 22~23.

[10] 徐 慧, 伍晓赞, 程仕平等. 复合裂纹的应力强度因子有限元分析[J]. 中南大学学报(自然科学版),2007,38(1): 79-83.

[11] 汤庆辉, 叶 彬. 光弹性法测定多裂纹应力强度因子研究[J]. 洪都科技,2006,(1):12~15.

[12] 杨绿峰, 徐华,李冉,等. 广义参数有限元法计算应力强度因子[J]. 工程力学, 2009,(03):48-54

[13] 孔令超. G*积分理论求解应力强度因子及有限元法验证[J].机械强度，2008，30(1)157-161

[14] 樊鸿, 张盛. 求解应力强度因子的单元初始应力法[J]. 武汉理工大学学报,2008,30(11):90-93

[15] 王军,邱志平,王晓军. 应力强度因子的区间分析方法[J]. 航空学报,2008,29(3):611-615

[16] 董玉文,余天堂. 任青文直接计算应力强度因子的扩展有限元法[J]. 计算力学学报. 2008,25(1):72-77

[17] 吕国志. 关于蒙特卡罗方法求解应力强度因子问题的探讨[D]. 西安：西北工业大学. 2007.03

[18] 苏成, 郑淳. 应力强度因子计算的样条虚边界元法[J]. 工程力学，2007，24（8）：49-53

[19] 梁尚清, 黄其青, 殷之平等. 计算应力强度因子的无网格—直接位移法[J]. 机械强度，2006，28(3):383-386

[20] 王承强, 郑长良. 平面裂纹应力强度因子的半解析有限元法[J]. 工程力学，2005，22（1）33-37