任意三维裂纹扩展分析-0319

任意三维疲劳裂纹扩展分析

1.前言

在工程实际中，真实的结构总是存在众多缺陷或裂纹，对于一个含裂纹或缺陷的构件，多在其服役荷载远低于容许强度的情况下就发生了破坏。实际工程结构在经受长时间多因素综合作用下，产生变形、裂纹等缺陷，从而导致整个结构的失效。结构的失效主要由疲劳引起，其最终失效形式即为断裂，有大约80%以上的工程结构的断裂与疲劳有关，由疲劳引起的巨大经济损失及灾难性的后果不胜枚举。

我们通常不能仅仅因为某个构件出现了裂纹就简单的认为该构件不安全或不可靠，尤其是对于大型设备的重要构件，因为这将使企业耗费高昂的成本。对于出现的裂纹，以往多采用以下几种处理办法：一是对出现裂纹的构件进行更换，这对于含裂纹但仍能工作的构件是一个巨大的浪费。二是强行停止使用进行维修，这样会带来巨大的经济损失；三是冒险继续使用，但这样会带来巨大风险，甚至会造成人员伤亡。所以，人们更想知道，出现的裂纹是否会在既定载荷（包括疲劳载荷在内的任意载荷）下扩展成不安全或失效的临界尺寸，因此，出现了疲劳裂纹扩展分析。疲劳裂纹扩展分析是采用断裂力学的理论和方法对含裂纹等缺陷构件的失效过程进行分析，以评估产品的安全性和可靠性，可以进行损伤容限评估和剩余寿命预测等，已经在化工机械、飞行器、核工业等各个工程领域得到了广泛应用，并得到了世界各国政府及学术机构的重视。

2.疲劳裂纹扩展分析软件

在工程实践中，疲劳裂纹扩展分析已成为评估产品性能、改良产品设计和提高服役寿命的一个重要工具。目前，疲劳裂纹扩展分析主要有解析法和数值法这样两种方法，下面分别介绍这两种方法。

1）解析法

解析法主要依据相应的规范和经验公式，将复杂的三维问题简化为二维问题，并对复杂的裂纹形状和荷载状态进行简化，然后用经验的方法对裂纹安全性进行评估。但对于大量结构复杂的工程实际问题却无能为力，况且其简化后的分析准确度及是否真实逼近服役情况也值得探讨。

目前，工程上有几款基于解析法而开发的裂纹扩展分析软件，它们主要应用于航空标准结构的裂纹扩展分析，包括DARWIN、NASGRO、AFGROW等。这些软件内嵌了航空结构多种形式的标准裂纹库，通过修改相应的模型尺寸、边界条件、载荷、裂纹位置和尺寸等参数即可根据内含的公式或内插表快速得出断裂力学结果，用来计算或查找标准航空结构中给定裂纹尺寸、载荷和形状的应力强度因子，仅能计算裂纹库里已有的裂纹模型的应力强度因子，并且适用于相对简单的几何和载荷，往往忽略真实的条件，如温度、非平面裂纹、复杂形状的裂纹、几何形状复杂的部件、部件之间的接触、残余应力和局部应力集中等。如要获得较为准确的结果，需要利用实验数据或其它方法对计算结果进行修正，但修正系数的取值往往很难确定，要靠经验来判断，并不具备求解复杂结构中三维裂纹扩展的能力。

2）数值法

近年来，随着有限元软件的发展，基于数值法的裂纹扩展分析软件已成功应用于解决工

程实际问题。市场上已有几款用于裂纹扩展分析的商用软件，它对含裂纹等缺陷产品进行计算机仿真，模拟产品的失效过程，由于计算精度高，使用方便，在工程应用上使用较为广泛。

裂纹扩展分析软件是疲劳裂纹扩展分析的一个强有力工具，在工业设计和设备安全性评估中有着重要的地位和作用，有着巨大的市场和广阔的发展前景。其主流发展方向是基于有限元法的三维裂纹扩展分析软件，包括基于扩展有限元法（XFEM）的ABAQUS、采用自适应网格划分的新一代FRANC3D、使用Crack-Block技术的Zencrack及等。

ABAQUS XFEM、新一代FRANC3D和Zencrack分别代表了三种裂纹网格建模技术，其中，Zencrack软件由于没有自己的用户界面，使得易用性受到很大限制，其采用的Crack-Block 技术并不能保证网格的质量，在复杂三维结构中的裂纹网格更新往往不能实现，因此，本文主要基于前两种软件来做探讨。

3.三维疲劳裂纹扩展分析的关键要素

在工程应用上，几乎所有的裂纹扩展都是三维的，使用数值方法计算三维裂纹扩展有许多困难和不确定性因素，其中有两个要素必须具备：

1）精确计算三维裂纹前缘的应力强度因子

2）根据裂纹在三维空间扩展的判据计算裂纹扩展寿命

任何裂纹扩展分析都必须基于断裂力学参数的计算和使用，而疲劳裂纹扩展寿命对应力强度因子等断裂力学参数较敏感，有一个公认的法则，即一个应力强度因子(SIF)25%的偏差，将会成倍增加(或减少)裂纹扩展的寿命。因此，如何获得精确的应力强度因子的结果，成为考虑损伤容限和剩余强度等现代设计的重要输入条件。

对于疲劳裂纹扩展寿命预测，通常需要两步来进行。首先，计算裂纹尺寸和应力强度因子之间的关系；其次，利用裂纹扩展速率公式来计算裂纹尺寸和施加载荷循环次数之间的关系。而如何获得裂纹尺寸和应力强度因子的关系曲线至关重要，这是准确计算疲劳裂纹扩展寿命的基础，而实际工程结构的裂纹扩展是空间三维的，如何合理地确定三维裂纹扩展的K vs. a关系曲线是准确计算三维裂纹扩展寿命的关键。

4.三维疲劳裂纹扩展分析的建模方法

进行三维裂纹扩展分析的前提是获得精确的应力强度因子的结果（我们将精确定义为与基准/解析解的误差在1-1.5%或更小的范围内），而利用有限元法计算精确的应力强度因子，是比较困难的，这就对其建模和计算方法提出了更高的要求。总体来说，计算精确的SIF的关键因素包括：在裂纹尖端植入奇异单元；在裂纹区域划分足够密的网格；使用保守的积分计算应力强度因子以及采用对称网格来减少局部离散误差等。

4.1在裂纹尖端植入奇异单元

对于传统的有限元来说，这意味着在裂纹尖端布置1/4节点的奇异单元，可以是退化的20节点六面体单元（a）或15节点的楔形单元，这两类单元均可给出精确的结果。

对于扩展有限元（XFEM ），这意味着奇异附加函数： []???

???=ψ2cos sin ,2sin sin ,2sin ,2cos θθθ

θθ

θr

4.2 裂纹区域的网格密度

对于裂纹区域的网格密度，有一个通用的“经验法则”，即，如果要采用有限元法来获得精确的结果，单元的尺寸应该至少10倍小于要划分网格的特征尺寸。对于局部裂纹区域网格划分来说这仍然是一个一般原则。

自适应网格划分方法会自动在裂纹尖端划分细化的网格，而在远离裂纹的位置使用普通结构分析中常用的网格单元尺寸，细化网格单元会过渡或“tie ”到远离裂纹的较大单元上。

在裂纹扩展过程中，自适应网格划分方法会继续沿着预测的裂纹路径来细化网格，保证

裂纹区域具有足够密的网格。

传统有限元分析中的单元尺寸原则同样适用于XFEM ，要想获得精确的应力强度因子结果，理想的局部裂纹区域的单元尺寸应该是至少10倍小于特征裂纹尺寸。

4.3 计算应力强度因子的保守积分

目前，计算应力强度因子的方法包括J-积分和M 积分（也称为交互积分），这种两种方法的等效域方程是计算应力强度因子最精确的方法。其中，J-积分适用于纯I 型裂纹问题，只能计算单调加载的情况，而M-积分可分别给出各项同性和一般各向异性材料中KI, KII,和KIII 的值，比J 积分更具通用性。

J 积分和M 积分的表达式分别为：

??????? ??-??=ds x q W x u J j

j

ij 11δσ

??????? ??-??+??=ds x q W x u x u M j j ij ij 1)2,1(1)1()2(1)2()1(δσσ

新一代FRANC3D 采用M-积分来计算应力强度因子，它对围绕裂纹前缘的两个单元环执行保守积分计算，积分域包括一个15节点奇异楔形单元的内环和一个20节点六面体单元的外环。FRANC3D 采用的自适应网格划分，还会在裂纹前缘周围布置第三个环，由六面体单元组成，但不参与积分计算。

保守积分的等效域方程需要使用一个物理上可以被解释为虚拟裂纹扩展的加权函数。在FRANC3D 中，通过在裂纹前缘生成没有径向畸变的奇异单元来实现该加权函数的定义，而在奇异单元周围的单元环均为标准的20节点六面体单元。

4.4 局部网格对称

由于裂纹前缘的积分只在少量的单元中进行，因此，这些单元对离散误差是高度敏感的。因为裂纹前缘区域是关于裂纹面对称和反对称的，将网格关于裂纹面和其垂直面进行对称将消除一些离散误差。不幸的是，XFEM 难以利用此对称来抵消离散误差，因此其计算精度要低于FRANC3D 的精度。

5. 任意三维裂纹扩展分析的步骤

新一代FRANC3D 分三个步骤来预测三维裂纹扩展：计算裂纹前缘上每个节点的局部裂纹扩展方向，或称扭转角度；计算每个节点的局部裂纹扩展距离；对扩展之后的新裂纹前缘进行光顺化处理，以减少不必要的数值“噪音”，并将裂纹前缘延伸到结构自由表面外。

5.1计算扭转角度（Kink Angle）

在大多数情况下，预测裂纹的扭转角度和扩展轨迹会相对来说会非常容易，特别是对疲劳裂纹的扩展计算。计算扭转角度有很多准则可以选择，包最大张应力准则、最大能量释放率准则、局部对称准则、最小应变能密度准则等。这些准则，包括任何倾向于最小化KII的准则，都能给出相似的并且精确的预测。下图为FRANC3D盲算获得的裂纹扩展趋势的结果和实际观察到的裂纹扩展的比较。

盲算结果观察到的裂纹扩展

然而，有一些情况会使得预测裂纹的扩展轨迹变得更加复杂，如裂纹扩展各向异性的阻力、非比例载荷（包括残余应力）、非常高的II型载荷（对有些材料）等。如下图为20世纪90年代初，波音公司做过的一个全尺寸窄体板测试，展示了各向异性材料对裂纹扩展趋势的影响。

5.2计算局部裂纹扩展距离

一般情况下，裂纹前缘上每个节点的扩展距离是不同的，使用疲劳裂纹扩展速率模型来计算局部裂纹扩展距离，有两个选项可供选择：

1) 指定位于应力强度因子中值（median）上的节点的扩展增量，所有其它节点的扩展

距离通过适当缩放获得，计算公式可表达为

()

()???

?

?

?

?

?

?

?

?

=

?

...

,

,

,

,

...

,

,

,

,

max

max

th

m

th

i

m

i K

K

R

K

f

K

K

R

K

f

a

a

；

2) 指定载荷的循环次数，从裂纹扩展速率公式直接计算所有节点的扩展距离，可表达

为

()...

,

,

,

,

max th

i

i

K

K

R

K

Nf

a?

?

?

=

?

。

指定扩展增量指定载荷循环次数

对于这两种方法来说，应注意局部裂纹扩展的增量不能设置的太小，这样会导致效率的

低下；也不能设置的过大，这样会导致裂纹前缘的不稳定。因此，应根据分析目的不同而合理设置裂纹局部扩展增量，如果要进行进精确寿命计算，可采用较小的增量，如果只需要查看裂纹扩展的趋势，则可设置较大的增量，以提高效率。

小的裂纹扩展步（稳定的裂纹前缘形状）

大的裂纹扩展步（不稳定的裂纹前缘形状）

5.3裂纹前缘拟合

预测的裂纹前缘点是一系列的数值计算的结果，将这些点连线作为新裂纹前缘有时会出现振荡，甚至呈锯齿形，利用多项式曲线来拟合裂纹前缘可以消除这些振荡。我们已经尝试了各种多项式和样条拟合曲线，遗憾的是，没有任何一种曲线能很好地拟合所有的情况，需要根据特定情况进行具体分析后选择合适的拟合方法。

6. 任意三维疲劳裂纹扩展寿命的计算和挑战

疲劳裂纹扩展寿命的计算通常需要两步来进行：

1) 计算裂纹尺寸和应力强度因子之间的关系；

2) 整合裂纹扩展速率模型来计算裂纹尺寸和施加载荷循环次数之间的关系，公式为()???= ,,,),(max th K K R a K f da

N 。

其中，如何合理地确定三维裂纹扩展的K vs. a 关系曲线是任意三维裂纹扩展寿命计算的关键。因为K vs. a 曲线是一个单自由度的曲线，而每条裂纹前缘上的应力强度因子分布是不均匀的，选取哪个点的结果作为本扩展步的K 值是不容易确定的；另外，对于裂纹扩展的距离a ，如果是两维问题，可以很容易获得，但对于三维裂纹扩展来说，问题变的很复杂，因为裂纹扩展路径是三维的，裂纹前缘上的每个节点的扩展距离和扭转方向都不一样，很难准确定义裂纹扩展的尺寸。

通常的做法是采用某种规则定义穿过裂纹前缘的一个“路径”，提取裂纹前缘与该路径交点处的K 值作为应力强度因子的历程数据。将该路径的长度定义成裂纹扩展的距离a ，从而生成K vs. a 曲线。这些规则包括：

1) 每条裂纹前缘恒定位置处（如每条裂纹前缘的中间位置）的直线连线形成路径

2) 找到距离前一条裂纹前缘指定位置最近的点，以此类推，形成路径

3) 定义一个平面，与裂纹前缘交点的连线形成路径

但是，利用每种规则都可以产生多个路径，应该如何选择单自由度路径，它如何影响三维疲劳裂纹扩展的寿命，还有待研究，因为，不同的路径预测出来的寿命并不同。下图为在拉弯载荷下的角裂纹，选择不同路径预测出来的寿命曲线对比情况。

另外，这些准则在一些复杂的情况下并不适用，或并不精确，如下图所示的情况，定义裂纹扩展的距离就很难。

为了解决以上问题，康奈尔大学断裂工作组提出了一个改进的新方法，比传统的方法更具有通用性，步骤如下：

1.对于裂纹前缘上的每个节点，找出其垂直投影到下一个裂纹前缘上的交点；

2.通过插值计算出这些点的K；

3.假设DK在当前和下一个裂纹前缘的点之间为线性变化，并通过积分获得预测的载

荷循环次数；

4.对每个裂纹扩展步的每个裂纹前缘单元预测的N进行平均，得到一个“N vs. 裂纹

扩展步”曲线。

但是，这种方法仍然是通过减少全三维的K数据来预测疲劳寿命，知识在现有的方法基础上进行得一个改进，最好方法仍然是一个开放性问题，还有待进一步研究解决。

7.结论

1.选择一个数值方法（自适应网格划分、XFEM等）只是执行疲劳寿命预测的第一步；

2.采用适当的单元、网格划分和积分技术，有限元法是可以获得精确的SIF结果的；

3.预测裂纹扭转方向和扩展轨迹是相对容易的；

4.预测的裂纹前缘形状对选择的扩展步长比较敏感；

5.通过减少全三维的K数据来预测疲劳寿命的最好方法仍然是一个值得解决的问题。