基于扩展有限元和极限学习机预测斜裂纹板疲劳寿命*

宋威震，杨庆勇，倪 伟，胡永敢，何水龙，蒋占四

(桂林电子科技大学 机电工程学院，广西 桂林 541004)

0 引言

Paris[1]第一次提出了疲劳裂纹增长的模型用于疲劳寿命，在Paris模型中，需要精确的应力强度因子去预测裂纹增长和疲劳寿命[2]。应力强度因子[3]通过裂纹的长度来计算，虽然有很多方法来计算应力强度因子，但数值方法是用的最多的。数值方法比如：有限元法[4]、无网格伽辽金法[5]、边界元法[6]。在这些方法中有限元法是用得最多和最有力的工具，但是有限元法在处理像裂纹这种强不连续问题时，必须将裂纹尖端设置为单元结点，使得裂纹面和单元的边重合，并且需要在裂纹尖端附近的高应力区划分很高密度的网格。

为了解决有限元法这种问题，Belytschko[7]教授和其课题小组提出了扩展有限元法。从提出以后，该法已经被许多研究人员应用于断裂力学问题的模拟。

前欧拉法是一种经典的预测模型，是疲劳裂纹扩展的一种常用方法，然而，前欧拉法存在一个缺点，例如它需要一个小的步长来保持精度[8]。近似模型可以解决这一问题，如多项式响应面[9]、Kriging[10]、支持向量机[11]、极限学习机[12-13]等。前欧拉法预测裂纹疲劳寿命，缺点是步长小精度低。如果只用中点法积分Paris模型预测裂纹疲劳寿命，当中点处的裂纹尺寸未知时，无法计算中点处应力强度因子，用极限学习机和中点法相结合可以解决这一问题。所以本文选用极限学习机和中点法相结合来预测裂纹板疲劳寿命。

1 疲劳寿命预测

1.1 疲劳裂纹扩展模型

Paris创建了疲劳裂纹扩展模型，并得到了广泛的使用，本文使用Paris模型。Paris模型如下：

(1)

ΔK=Kmax-Kmin

(2)

其中，Kmax和Kmin在给定的加载循环中，在最大和最小施加应力的位置处的应力强度因子。

1.2 等效应力强度因子

本文主要考虑斜裂纹，斜裂纹的应力强度因子存在Ι型应力强度因子和ΙΙ型应力强度因子，即ΔKΙ和ΔKΙΙ，这时，就必须考虑等效应力强度因子，本文选用Tanaka[14]提出的方法，等效应力强度因子定义为：

(3)

式中，ΔKΙ和ΔKΙΙ分别是Ι型和ΙΙ型应力强度因子范围。

1.3 前欧拉法

目前，裂纹扩展问题最简单的计算方法是前欧拉法。前欧拉法为：

yn+1=yn+hf(xn,yn)

(4)

其中，h是从xn到xn+1的距离，f(xn,yn)是y(x)在xn处的斜率。因此，从xn到xn + 1的值采用xn处的斜率进行线性近似，其只需要当前步的数据就可以预测下一步的数据点，然而其需要进行很多次的迭代运算，运算量大。前欧拉法还有一个缺点就是要维持一定的精度必须保持较小的步长增长。

1.4 中点法

有许多方法允许较大的步长用于积分，本文选用中点法。中点法[15]通常比前欧拉法有更好的精度。因为它是在当前点xn和下一个数据点xn+1之间的中点处取得斜率。即xn+1/2处的斜率。公式可以表示为：

(5)

该方法的精度与h2成正比，因此，精度相同的条件下中点法比前欧拉法允许更大的步长。

2 积分疲劳裂纹扩展模型

2.1 疲劳裂纹扩展模型的数值积分

本文在模拟裂纹扩展过程中，设裂纹在扩展时步长为定裂纹尺寸[8]扩展，即固定的裂纹增量Δa=1mm、Δa=2mm。用前欧拉法积分Paris模型，步长为定裂纹尺寸扩展时，每周期疲劳寿命循环数ΔNi可以表示为:

(6)

其中，i是当前量，Δa为裂纹增长量，Ni是当前裂纹尺寸时的循环数，Ni+1是下一裂纹尺寸时的循环次数，ΔNi为Ni和Ni+1的差值，是两个评价点之间的周期数，ΔKi为循环次数Ni处的应力强度因子幅值。

虽然前欧拉法可以用来积分Paris模型，但是使用当前点的斜率预测下一个裂纹尺寸的应力强度因子，将导致较大的误差。而且前欧拉法存在收敛速度慢的缺点。因此，需要一种积分方法允许大步长又不降低精度。

中点法用来积分Paris模型，对于裂纹扩展时固定裂纹增量Δa，每周期疲劳寿命循环数ΔNi可以表示为：

(7)

其中，ΔKi+1/2是当前裂纹尺寸和下一裂纹尺寸中点所对应的应力强度因子幅值。当裂纹尺寸i+1/2处的尺寸未知时无法计算ΔKi+1/2，后面提出一种预测模型计算应力强度因子幅值使其满足中点积分法。图1为前欧拉法和中点法示意图。

图1 前欧拉法和中点法示意图

2.2 基于扩展有限元和极限学习机的预测裂纹剩余疲劳寿命步骤

本文采用极限学习机工具包ELM预测裂纹在i+1/2处的应力强度因子。极限学习机预测模型中D为训练样本，在i+1/2处裂纹的长度为测试样本。

(8)

其中，xi为裂纹长度，yi为应力强度因子(由扩展有限元法求解得到)，D作为训练样本，L为样本数。通过极限学习机得到的测试样本输出值即为裂纹i+1/2处对应的应力强度因子。

预测裂纹剩余疲劳寿命基本步骤：首先进行裂纹尖端应力强度因子的扩展有限元迭代求解，然后根据预测模型建立裂纹长度与应力强度因子(扩展有限元法获得)的非线性关系，求取中点应力强度因子(即式(7)中Δki+1/2);最后，根据应力强度因子与疲劳寿命的函数关系(即式(7))预测疲劳寿命循环数N。

3 数值分析

3.1 张力作用下的中心斜裂纹

图2 张力作用下中心斜裂纹

该模型的解析解为：

(9)

a为裂纹尺寸的一半。

图3 应力分布σyy

图4中，计算值是由扩展有限元法计算得到，理论值是由式(9)计算得到。从图4可以看出，Ι型应力强度因子的增长速率比ΙΙ型快。由扩展有限元法计算的应力强度因子的值与理论公式计算的应力强度因子的值基本吻合，证明了该程序的正确性，也证明了网格大小划分合适，为下面其他计算提供保障和正确性。

良肢位是为了保持肢体的良好功能而将其摆放在一种体位或姿势,是从治疗护理的角度出发而设计的一种临时性体位,能有效的固定偏瘫后的关节,也能预防患者在以后的生活中出现病理性运动模式[4]。通过早期对脑卒中偏瘫患者进行良肢摆放,发现可以有效的缓解患者的肌肉紧张,使患肢肌肉处于放松的状态,减少了患肢肌肉发生痉挛的机率,促进患肢的血液循环,有效预防偏瘫后遗症的发生,提高了患者治疗的积极性和患者的生活质量,增加了患者的护理服务满意度,使护患关系变得更加和谐,也让患者恢复了自信,较少了患者焦虑、抑郁的心里状态,对未来生活充满希望。

图4 应力强度因子和裂纹长度

图5 裂纹扩展长度和应力强度因子

图6 裂纹扩展速率和裂纹扩展长度

图7显示的是裂纹扩展速率和应力强度因子，从图中可以看出，应力强度因子越大，裂纹扩展速率越大。

图7 裂纹扩展速率和应力强度因子

图8显示的是裂纹扩展按定裂纹增量Δa=1mm(2mm，3mm)增长时，不同的裂纹增量和疲劳寿命之间的关系。从图中可以看出，当步长(Δa的值)固定时，裂纹增量Δa为3mm，Δa为2mm，Δa为1mm时，前欧拉法(式(6))预测疲劳寿命循环数大约为29482，26505，24454。对比标准值即裂纹增量Δa为0.1mm时，疲劳寿命循环数大约为21927，前欧拉法预测的误差分别为34%，20%,11.5% 这就造成了当定步长裂纹尺寸增长时，步长(Δa的值)较大时，前欧拉法预测裂纹剩余疲劳寿命误差较大。

图8 前欧拉法和疲劳寿命循环数

从图9可以看出,当步长尺寸固定时，Δa=2mm，用极限学习机结合中点法预测裂纹剩余疲劳寿命大约为22681，对比标准值Δa=0.1mm时的疲劳寿命循环数21927，误差为3.5%，对比相同步长尺寸Δa=2mm时，前欧拉法预测裂纹剩余疲劳寿命为26505，误差为20%。通过对比相同的步长尺寸即裂纹增量Δa，可以看出极限学习机结合中点法比前欧拉法有更高的精度，在误差方面，极限学习机比前欧拉法误差更小，同时在大步长尺寸(Δa的值较大时)增长时，极限学习机结合中点法精度仍然很高。

图10显示的是中心斜裂纹扩展后的大致路径，大致呈现Z形，这与文献[16-18]基本一致，证明本方法程序的正确性。

图9 两种方法和疲劳寿命循环数

图10 斜裂纹扩展路径

3.2 张力作用下的边缘斜裂纹

(10)

图11 张力作用下板边缘斜裂纹

图12显示的是固定裂纹长度a，改变倾斜角度β，Ι型应力强度因子和ΙΙ型应力强度因子的大小变化，从图中可以看出在19°左右，KΙ和KΙΙ在数值上相等。

图12 应力强度因子和裂纹角度

图13 裂纹扩展长度和应力强度因子

图14 裂纹扩展速率和裂纹扩展长度

图15显示的是裂纹扩展速率和应力强度因子，从图中可以看出，应力强度因子越大，裂纹扩展速率越大。

图15 裂纹扩展速率和应力强度因子

图16显示的是裂纹扩展按定裂纹增量Δa=1mm(2mm，3mm)增长时，不同的裂纹增量和疲劳寿命之间的关系。从图中可以看出，当步长(Δa的值)固定时，裂纹增量Δa为3mm，Δa为2mm，Δa为1mm时，前欧拉法(式(6))预测裂纹剩余疲劳寿命循环数大约为2494，2374，2293。对比标准值即裂纹增量Δa为0.1mm时，预测裂纹剩余疲劳寿命循环数大约为2212，前欧拉法计算的误差分别为12.5%，7.3%,3.6% 造成了当定步长裂纹尺寸增长时，步长(Δa的值)较大时，前欧拉法预测裂纹剩余疲劳寿命误差较大。

图16 前欧拉法和疲劳寿命循环数

本文同样用了极限学习机和中点法相结合的方法，来预测裂纹剩余疲劳寿命。结果如图17所示。

图17 两种方法和疲劳寿命循环数

图18 斜裂纹扩 展路径

从图17可以看出,当步长尺寸固定时，Δa=2mm ，用极限学习机结合中点法预测裂纹剩余疲劳寿命大约为2154，对比标准值的寿命循环数2212，误差为2.7%，对比相同步长尺寸Δa=2mm时，前欧拉法预测裂纹剩余疲劳寿命为2374，误差为7.3%。通过对比相同的步长尺寸即裂纹增量Δa，可以看出极限学习机结合中点法比前欧拉法有更高的精度，在误差方面，极限学习机比前欧拉法误差更小，同时在大步长尺寸(Δa的值较大时)增长时，极限学习机结合中点法精度仍然很高。

图18显示的是边缘斜裂纹扩展后的大致路径，与文献[16-18]基本一致，证明本方法程序的正确性。

4 结论

针对前欧拉法在预测裂纹疲劳寿命具有小步长、精度低的问题，提出极限学习机和中点法相结合的方法，用于预测含裂纹板结构的疲劳寿命。进行了两个数值分析，数值分析结果表明，本文提出的方法进行剩余疲劳寿命预测，精度高并且允许大步长，提高精度的同时也降低了有限元计算量。