基于比例型Paris 公式和逆高斯过程的金属疲劳裂纹扩展随机模型

陈 龙，黄天立，周 浩

(中南大学土木工程学院，湖南，长沙 410075)

金属构件的疲劳可靠性及剩余寿命预测一直是抗疲劳设计的主要任务之一，而疲劳裂纹扩展分析作为金属构件疲劳寿命预测的主要内容，其相关问题一直是该领域的研究热点。由于外部的不确定因素和材料内部的非均质等影响，疲劳裂纹扩展过程中的波动范围通常较大，如著名的Virkler 金属疲劳裂纹扩展试验数据[1]等。探讨这种随机性的原因及规律，对于桥梁、飞机、船舶等安全性和可靠性要求较高的构件设计和寿命预测是十分重要的，也是必需的。

大量研究表明[2-5]，金属疲劳裂纹扩展过程中的随机性主要包括两种：1)材料内部阻力的不恒定、裂纹面的不规则导致裂纹扩展路径的不规则波动，处理此类随机性的通常做法是，在确定性疲劳裂纹扩展速率模型的基础上引入随机因子，借助数学理论和方法对概率寿命进行估计，形成了基于随机过程的疲劳裂纹扩展模型，蒋鸣晓等[2]、吴圣川等[3]对这类模型进行了综述，其中包括B-model、连续的马尔科夫过程模型、相关过程模型等，杨冰等[4]在Elber 方程的基础上，假设疲劳裂纹扩展速率服从一个平稳的对数正态分布，建立了相应的裂纹扩展随机模型，刘勇等[5]基于Paris公式，利用概率密度演化方法，解释了疲劳裂纹扩展过程中的波动行为，并对构件的疲劳寿命进行预测；2)材料微观结构、缺陷分布的不均匀性导致了样本间裂纹扩展速率的差异，处理这种随机性的通常方法是，在第一类模型基础上，随机化代表材料固有属性的参数(如Paris 公式中的材料参数C和m)，从而引入随机效应来表征样本之间的异质性。

单个样本的疲劳裂纹扩展路径在宏观上表现为循环荷载下的间接性增长，这类由一系列微小增量累积形成的损伤很适合采用独立增量过程模拟[6-7]。近年来，维纳(WP)、伽马(GP)和逆高斯(IGP)等具有独立增量性质的随机过程，由于其在退化建模方面表现出的优良特性，被广泛用于疲劳裂纹扩展的随机建模。Wang 等[8]采用幂函数作为退化轨道函数，联合维纳过程，建立了疲劳裂纹扩展随机效应模型，并推导了剩余寿命概率密度函数的近似解析解。Ye 等[9]利用指数函数描述疲劳裂纹的平均扩展路径，结合伽马过程，建立了相应的简单模型和随机效应模型，并采用了半参数方法估计模型参数。Peng 等[10]利用逆高斯过程描述疲劳裂纹的随机扩展，同时比较了幂函数和指数函数两种退化轨道下的模型拟合效果，结果表明，指数型退化轨道下的随机模型能更好地解释裂纹扩展的随机行为。Hao 等[11]提出了一个具有偏正态随机效应的扩展逆高斯过程模型，并将其应用到疲劳裂纹的剩余寿命预测，结果表明，相对于传统的逆高斯随机效应模型，其所提出的扩展模型拟合效果更好。

虽然采用独立增量过程描述裂纹扩展行为的研究较多，但仍存在以下两个问题：1)描述裂纹平均扩展路径的轨道函数仅采用两种形式，即幂函数和指数函数，由此建立的扩展模型实质上是一个统计模型，无法从退化机理层面描述疲劳裂纹的扩展过程；2)上述研究中都采用“时间间隔[t,t+δt] 内裂纹长度增量 δa与当前状态a(t)无关，仅取决于当前t和时间间隔 δt”的基本假设。实际上，疲劳裂纹扩展过程中不规则波动的本质是，由于靠近裂纹尖端的材料性质的局部变化。众多关于疲劳裂纹扩展的研究表明，裂纹扩展速率da/dt是一个与当前状态a有关的函数，裂纹长度增量 δa本质上是依赖当前状态a的。因此，为了更好地解释金属疲劳裂纹扩展中的随机性，有必要建立一个更贴近裂纹扩展机理的退化模型。

本文针对金属疲劳裂纹扩展特性，以“首次达到给定裂纹长度a的时间t(a)”作为随机过程的描述量，采用比例型Paris 公式作为退化轨道函数描述裂纹的平均扩展轨迹，建立基于逆高斯过程的单体样本疲劳裂纹扩展随机模型和考虑样本异质性的裂纹扩展随机效应模型，分别采用最大似然估计(MLE)和期望最大算法(EM)推导了单样本模型和随机效应模型的参数估计公式。最后，利用所提出的模型来拟合Virkler 金属疲劳裂纹扩展试验数据，并进行拟合优度分析。

1 疲劳裂纹扩展随机模型

1.1 比例型Paris 公式

迄今为止，在断裂力学基础上发展出来的疲劳裂纹扩展模型较多，其中Paris-Erdogan 公式[12]因其形式简单和实用，得到了工程界的普遍应用。其表达式如式(1)所示：

式中： da(t)/dt为裂纹扩展的速率；a为疲劳裂纹的长度；t为试件经受疲劳荷载的加载时间，通常以循环加载的次数表示；参数C′和m′代表材料参数，其值取决试件的材料和几何尺寸； ΔK(a)为应力强度因子，在实验条件和试件尺寸确定的情况下； ΔK(a) 为一个关于裂纹长度a的函数。

图1 给出了经典的Virkler[1]疲劳裂纹扩展数据。该试验样本采用68 块应用于飞机结构的2024-T3 铝合金板件，试件形状呈矩形，长度、宽度和厚度分别为558.8 mm、152.4 mm、2.54 mm。金属板在试验时受到频率20 Hz，应力幅为24.14 MPa，应力比为0.2 的正弦循环荷载。试验之前在样本中心形成长度为 2a0=18 mm的初始裂纹长度，裂纹在疲劳荷载作用下向两边扩展延伸，同时只记录一边的裂纹长度数据，当疲劳裂纹长度达到49.8 mm时，试验结束。

图1 Virkler 疲劳试验数据Fig. 1 Experimental dataset from Virkler fatigue test

图2 l og10C′ 和 m′的相关性Fig. 2 Correlation between l og10C′ andm′

1.2 单样本疲劳裂纹扩展随机模型

由于疲劳裂纹长度在宏观上随时间的推移不断增长，因此，以往的研究多将裂纹长度a(t)视为随时间t变化的一种随机过程，记做 {A(t),t＞0}，并假设其过程是一个独立增量过程，由此来模拟疲劳裂纹的间接性增长。然而实际上，疲劳裂纹长度的增量 δa之间是相互依赖的，这种建模方式与裂纹扩展的本质是相违背的。Guida 等[14]认为，研究疲劳裂纹随机扩展问题的最佳方法是，将“达到给定裂纹长度a所需要的时间t(a)”视为一个随机过程，记为 {T(a),a≥a0}。若随机过程{A(t),t＞0}是一个状态依赖过程，则相应的“逆过程” {T(a),a≥a0} 的时间增量 δt是相互独立的。反之，若 {T(a),a≥a0}被建模成独立增量过程，则描述状态随时间变化的随机过程 {A(t),t＞0}应该是一个纯状态依赖的过程[15]。

本文假设 {T(a),a≥a0}服从一个非齐次的逆高斯过程(IGP)，由于传统的IGP[16]描述的是状态随时间的随机变化的过程，因此，按照文献[15]对伽马过程的一般定义修改方法，本文对IGP 的一般定义进行修改，

1.3 考虑样本差异的疲劳裂纹扩展随机效应模型

在实际工程中，由于工艺、设计以及试验环境等因素的差异，同类样本不同个体之间的裂纹扩展过程会产生较大的差异，这种差异性在退化数据上体现为裂纹扩展路径的分散性，因此有必要对其进行描述。假设单路径随机扩展模型中的参数C为一个随机变量，将个体特有的随机效应集成到退化过程中以表征个体差异性，使退化模型能够同时描述裂纹扩展过程中的两种随机性。

为了保证参数C的非负性以及数学推导上可行性，Wang 等[17]建议随机参数应服从一个截断正态分布TN(ω,κ-2)κ＞0，概率密度函数如下所示：

首达时间增量 ΔT(δa)的边缘概率密度函数以及边缘累积分布函数可分别通过下列积分计算：

式中，BvN(·)表示二元标准正态累积分布函数，其余参数如下所示：

2 随机模型参数估计

2.1 单样本随机模型参数估计

2.2 考虑样本差异的随机效应模型参数估计

3 模型拟合优度检验

利用Kolmogorov-Smirnov 统计量D 衡量经验分布和理论分布之间的差异，同时获取假设检验的P 值，对于简单假设(25)：

1) 若P 值小于等于α ，则在显著水平α 下拒绝H0；

2) 若P 值大于 α，则在显著水平 α 下接受H0。

4 疲劳裂纹随机模型拟合优度分析

4.1 单样本疲劳裂纹扩展模型拟合优度分析

图3 材料参数的常规ML 估计值和去相关ML 估计值Fig. 3 conventional and decorrelated ML estimates of material parameters

图4 68 条疲劳路径的K-S 检验P 值Fig. 4 P-value of 68 fatigue crack growth paths by using Kolmogorov-Smirnov test

4.2 考虑样本差异的疲劳裂纹随机效应模型拟合优度分析

在实际试验过程中，同类产品不同个体之间在裂纹扩展路径上也会表现出较大的差异，通过随机化单体扩展模型中的参数C，令其服从一个截断正态分布，可以建立体现个体差异的随机效应模型。表1 给出了采用EM 算法对随机效应模型的参数结果。

表1 随机效应模型参数估计值Table 1 Parameter estimates of random effect model

在表1 模型估计的基础上，本节对随机效应模型一些关键量进行了统计拟合分析。图5 给出了均值函数和方差函数ML 估计和样本估计的对比图。当裂纹长度较小时(a≤20 mm)，均值和方差的ML 估计值和样本估计值相差较大，随机效应模型无法较好地预测样本值；当疲劳裂纹长度较大(a≥20 mm )时，均值E[T(a)]的ML 估计值和样本估计值之间的误差为±1.4%，模型和样本的拟合程度较好，而方差 Var[T(a)]的拟合程度相对于均值E[T(a)]较低。

图5 首达时间 T(a)的均值和方差的ML 估计以及样本估计Fig. 5 Empirical and ML estimates of the mean and variance of first timeT(a)

图7(a)给出了首达时间T(a)在10%、40%、70%、90%时的百分位数ML 估计和样本估计的比较图；图7(b)给出了ML 估计和样本估计的相对误差图，可以看到，裂纹扩展前期，模型估计值和样本值相差较大，拟合程度并不高；当裂纹长度a≥20 mm，相对误差基本上稳定在±2.5%左右，模型的拟合效果较好。

图6(a)给出了首次到达给定裂纹长度a=20 mm,25 mm, 30 mm, 35 mm, 40 mm 的时间T(a)的累积分布函数的ML 估计和样本估计的对比图；图6(b)给出了时间t=0.8 ×104cycle,1.2 ×104cycle,1.6×104cycle,2.0 ×104cycle,2.2 ×104cycle时，疲劳裂纹的扩展长度A(t)的累积分布函数的ML 估计和经验估计的对比图。从图6 可以看出，随机变量T(a)和A(t)的CDF 的ML 估计曲线和样本估计值的拟合程度较好。

图6 首达时间 T(a) 和 裂纹长度 A(t)的CDF 的ML 估计以及样本估计Fig. 6 Sample and ML estimates of the CDF of first passage time T (a) and crack lengthA(t)

图7 首达时间 T (a)的百分位数ML 估计与样本估计对比图以及相对误差图Fig. 7 Sample and ML estimates of percentiles of first passage time T(a) and relative percent differences

为对随机参数C的截断正态分布假设进行检验，图8 给出了随机参数C的Q-Q 图。从图8 可以看出，样本点基本分布在对角线附近。通过K-S检验得到的假设检验P 值为0.777，大于0.05 的显著性水平 α。因此，可以认为随机参数C服从一个截断正态分布。

图8 随机变量的C 的截断正态分布Q-Q 图Fig. 8 The normal Q-Q plot under normal distribution hypothesis for C

5 结论

针对金属疲劳裂纹扩展过程中的不确定性，采用比例型Pairs 疲劳公式，建立了基于逆高斯过程的单样本疲劳裂纹扩展随机模型和考虑样本差异的裂纹扩展随机效应模型，分别采用最大似然估计和EM 算法推导了单样本模型和随机效应模型的参数估计公式。最后将两种模型分别应用于Virkler 疲劳试验数据的拟合，并对拟合效果进行检验，可得出如下结论：

1)采用比例型Paris 公式的随机模型能够有效消除参数之间的强相关性。

2)通过概率积分变换判断疲劳裂纹增量是否符合逆高斯分布族，衡量了单样本疲劳裂纹扩展随机模型的拟合效果。结果表明，没有明显的证据显示单个疲劳裂纹的扩展过程不服从一个逆高斯过程。

3)利用图形法对一些关键量进行拟合分析，判断了考虑样本差异的疲劳裂纹扩展随机效应模型整体的拟合效果。结果表明，所选择关键量的ML 估计值都能够较好地拟合样本估计值。

由此表明，基于比例型Pairs 公式和逆高斯过程的随机模型能够有效地分析和解释金属疲劳裂纹扩展过程中的不确定性。