钢结构中疲劳裂纹增长概率模型和Bayesian动态预测

陈梦成，方 苇，杨 超，谢 力

(1.华东交通大学 省部共建轨道交通基础设施性能监测与保障国家重点实验室，江西 南昌 330013；2.华东交通大学 土木建筑学院，江西 南昌 330013)

近50年来，随着我国改革开放和经济实力提升，基础设施建设得到蓬勃发展。大量既有桥梁结构在使用环境和使用荷载作用下结构性能逐渐下降，产生疲劳裂纹并不断扩大，严重影响了桥梁的耐久性和适用性，甚至有部分桥梁尚未达到设计使用年限就已经出现重大安全隐患。准确地预测桥梁工程中钢结构性能退化和疲劳寿命评估是进行桥梁结构维修与决策、延长桥梁使用寿命的关键[1-2]。工程结构可靠度理论的研究是个长久的课题，早在20世纪70年代发达国家就已开始进行服役工程构件及结构系统的可靠性分析及预测的研究[3-4]，从20世纪90年代开始，结构时变可靠性评估已成为国内外学者的研究热点[5-6]，但结构时变可靠性的计算方法还远不成熟。

结构可靠性的研究通常是通过考虑影响结构性能退化的各不确定因素，推算结构抗力的时变概率特性[7-10]，从而建立结构抗力随机过程模型[7-8,11-12]。但在具体的实践中，因桥梁结构耐久性问题的时间周期长、影响因素复杂、不确定性大、不同结构从材料到构件再到结构设计的差异性等因素，很难用某一种模型来概括。已有的工程结构失效预测和可靠性评估模型也无法做到充分考虑实际工程中各种复杂因素对各参数的影响。因此研究结构可靠性就要求建立合理的结构性能随机过程模型[13]。现代无损检测技术常用于桥梁结构系统的健康监测中，受此启发，众多研究者提出利用健康监测数据进行桥梁结构的疲劳寿命预测[14-15]，而贝叶斯更新理论可以实现耦合服役工程结构性能退化过程(模型的不确定性)和监测数据，动态更新结构性能退化随机过程模型。Bayesian更新理论的数值计算方法中，马尔卡夫链蒙特卡罗(Markov Chain Monte Carlo，简称MCMC)抽样方法最为流行，因为它利用计算机的强大计算能力和仿真能力[16-17]。

疲劳是导致桥梁结构失效的原因之一，而疲劳裂纹在既有桥梁结构中非常普遍，其增长过程是一个随机过程，疲劳裂纹的增长会引起结构性能退化。本文拟针对桥梁结构中因疲劳裂纹损伤而引起的结构性能退化以及时变可靠性预测的科学问题，采用贝叶斯更新理论，开展相关研究。文中介绍了Bayesian更新理论在疲劳裂纹增长模型的参数更新和构件寿命预测的应用中如何实施，并进行了具体的推导论证。最后结合具体算例，采用Matlab编程，对疲劳裂纹损伤构件在未来一段时间内的裂纹增长和时变可靠性进行预测，以检验本文方法的可行性和有效性。

1 Bayesian更新定理

(1)

依据上面说明，Bayesian更新公式可用文字表述为：更新概率(后验概率)=新观测信息的条件概率(似然度×先验概率)/新观测信息的全概率(标准化常量)，也就是说，后验概率与先验概率和似然度的乘积成正比。另外，比例P(B/Ai)/P(B)有时也被称作标准似然度，因此，Bayesian更新公式可进一步用文字表述为：后验概率= 标准似然度×先验概率。

2 构件疲劳裂纹损伤性能随机退化建模

2.1 疲劳裂纹扩展确定性模型

对于一个含有表面裂纹的钢构件，依据线性弹性断裂力学(LEFM)理论，裂纹端部应力强度因子K在某一应力水平σ作用下为

(2)

式中：Y为应力强度修正系数，反映了构件和裂纹的几何形状；a为半裂纹长度。

在构件承受静荷载(一般使用环境)时，只有其应力水平达到临界应力σc时，即裂纹端部的应力强度因子达到临界值Kc时，才会立即发生失稳断裂，此时相对应的临界裂纹长度为ac。当静应力水平降低到σ0(σ0<σc)，则构件不会发生失稳破坏。但如果构件承受一个与静应力σ0大小相等的往复循环应力，则初始裂纹a0在循环应力σ0作用下发生缓慢扩展。根据著名的Paris半经验公式[18]，疲劳裂纹扩展速率为

(3)

假设疲劳裂纹扩展过程中，Δσ和Y保持不变，则可取ΔK(a)=a，积分式(3)可得

(4)

式中：a(0)为N=0时构件的初始疲劳裂纹长度；a(N)为荷载循环次数为N时的裂纹长度。

2.2 疲劳裂纹损伤性能退化概率模型

令式(4)中材料参数C=θ1、m-1=θ2，并假设这些参数均为不确定性参数，则对式(4)两边取自然对数，并令y=ln(a(N)/a(0))，则有

(5)

式中：yij为第i(i=1,2,…,n)个疲劳裂纹在荷载循环到tj(j=1,2,…,m)=N时对应的相对疲劳裂纹扩展长度，工程结构构件相应性能退化；tj为时间变量，表示作用在构件上的疲劳荷载循环次数，为性能变量的观测时间；η(·)为描述性能变量变化的轨迹函数，它是t的单调非线性函数，在许多情况下它也是参数θ(θ1i,θ2i)的非线性函数；θ(θ1i,θ2i)∈Θ⊂Rq为性能退化函数η(·)中的未知随机参数向量；q为参数向量位数；Θ为参数向量的取值空间；εij为反映性能退化过程中第i个疲劳裂纹在tj时刻相对扩展长度的不确定性的误差项。

参数向量的“随机性”应该理解为，对于同一轨迹的性能退化，模型参数的取值是确定的；对于不同轨迹的性能退化，模型参数的取值是随机的，这种随机性可以由某个概率分布刻画。为此，进一步假设：①对应特定的tj，随机误差项序列εij，i=1,2,…,n独立同分布，也就是说，εij在空间上是独立的；②对于特定的疲劳裂纹i，随机误差项序列εij，j=1,2,…,m是相关的；③随机误差项序列εij服从均值为0，未知协方差矩阵为Σε的多维随机变量正态分布，即N(0,Σε)。

定义yi=(yi1yi2…yim)T和ηi=(ηi1ηi2…ηim)T，其中T为矩阵转置；用fMVN(x/μ,Σ)表示均值为μ，协方差为Σ的多维随机变量正态分布概率密度函数，因此，在给定ηi的情况下，由式(4)可知

(6)

式中：Σε为随机误差项m×m阶的协方差方阵，其中矩阵元素等于ρklσkσl，(k=1，2，…，m；l=1，2，…，m)；ρkl为荷载循环次数k次和l次时随机误差项之间的相关系数；σk和σl分别为荷载循环次数k次和l次时随机误差项的标准方差。

3 疲劳裂纹损伤性能退化模型随机参数的Bayesian更新

3.1 先验分布

3.2 性能退化过程似然函数和后验分布

对于特定的疲劳裂纹i，假定在tj(j=1,2,…,m)时刻有损伤性能退化量测值yi1,yi2,…,yim。由于随机误差项序列εij，j=1,2,…,m是独立同分布的随机变量，在给定θ和σε的条件下，性能退化样本yi1:m的条件联合分布函数或似然函数为

(7)

在已知联合先验分布π(θ,σε)的基础上，得到性能退化量测值D(=yi1,yi2,…,yim)后就可以利用Bayesian更新式(1)来推断后验分布p(θ,σε|D)，即有

(8)

式中：Θ、Ω分别为θ、σε的取值空间。

θ1(1)从全条件分布π(θ1(0))π(θ2(0))π(θ3(0))L(y1:m|θ1(0),θ2(0),θ3(0))产生；

θ2(1)从全条件分布π(θ1(1))π(θ2(0))π(θ3(0))L(y1:m|θ1(1),θ2(0),θ3(0))产生；

θ3(1)从全条件分布π(θ1(1))π(θ2(1))π(θ3(0))L(y1:m|θ1(1),θ2(1),θ3(0))产生；

⋮

以上Gibbs抽样技术的关键在于仅考虑了单变量的条件分布。抽样样本可以通过选择一个建议分布来生成。这里选择随机游走链作为建议分布为

(9)

(10)

重复上述迭代过程直到各分量Θk+1的遍历均值稳定后，可生成一系列所需要采集的样本{Θ1,Θ2,Θ3,…}。利用这些后验(更新)样本，即可估计后验分布和随机参数的置信区间。

4 疲劳裂纹损伤性能退化与时变可靠性评估

4.1 单根疲劳裂纹极限状态方程

建立构件单根疲劳裂纹长度安全余量方程，作为进行可靠性分析的极限状态方程为

Z=acr-y(tk+Δt)≥0

(11)

式中：Z为疲劳裂纹长度的安全余量；acr为给定的疲劳裂纹临界值；y(tk+Δt)为疲劳裂纹在tk+Δt时刻的预测值，tk为当前时刻；y(tk+Δt)=η(tk,θ1,θ2)+εk+1，其中εk+1为疲劳裂纹损伤性能退化轨迹模型在tk+Δt的随机误差项，假设它与ε1,ε2,…,εk独立同分布，即：εk+1～N(0,σε)；Δt为未来时间区间长度。

4.2 疲劳裂纹损伤构件时变可靠性Bayesian预测

要了解tk+Δt时刻的疲劳裂纹损伤性能退化后的构件可靠性，就必须依据当前时刻以及过去的疲劳裂纹长度量测数据，对其可靠度进行评价。根据式(11)，构件在未来某一时刻(tk+Δt)的可靠概率为

(12)

式中：g(θ1,θ2,σε,tk,Δt,acr)=Z=acr-y(tk+Δt)，由4.1节可知服从N(0,σε)分布。

一般说来，要解析求解式(12)积分通常很困难，因此，我们用Bayesian更新式(8)、式(9)、式(10)所建立的抽样方法，以解决预测和可靠度分析问题。在获得未知随机参数和随机误差项的后验分布后，可将式(12)的求解转化为

(13)

每次达到疲劳试验的设计试验节点，就会收集到一个新的量测数据。利用收集到的历史量测数据，相应地更新似然函数，即可利用式(13)计算未来一段时间Δt(可以取不同值)内疲劳裂纹损伤构件的可靠性分析和预测，并绘制一条时变可靠度变化曲线。

5 算例

为了说明以上建立的Bayesian更新方法和MCMC抽样技术的应用，以文献[21]中的疲劳裂纹扩展量测数据为例，对疲劳裂纹损伤构件的时变可靠性进行预测。疲劳试验试件共21个，试验时长(疲劳荷载循环次数)为0.12百万次，每隔0.01百万次对每组疲劳裂纹扩展长度进行量测。所有试件在疲劳试验开始前的初始裂纹长度均为22.86 mm，设定试件疲劳失效时裂纹长度的临界值acr为40.64 mm。21个试件的疲劳裂纹扩展量测数据见表1[21]，与试件相应的疲劳失效时间见表2[22]。由表1可知，疲劳试件13～21在疲劳试验终止时，裂纹长度尚未达到临界值。利用表1量测数据、最小二乘法和统计学知识，对式(5)中的模型参数进行近似估计，可得到模型参数的估计值θ1和θ2，参见表2，表2中还列出了文献[22]给出的疲劳失效时间预测。

表1 疲劳裂纹扩展长度量测数据[21] mm

表2 疲劳裂纹损伤性能退化轨迹模型参数与疲劳失效时间预测

5.1 算例1

选择第9条试件的疲劳裂纹损伤性能退化轨迹作为分析疲劳损伤构件时变可靠性分析及预测的研究对象，其余20条试件的退化轨迹用来提取退化轨迹模型随机参数θ1和θ2的先验信息。

依据第2、3、4节的理论，采用Matlab语言编制了Bayesian更新MCMC法源程序，该程序可读性强，而且力求忠实MCMC理论框架。在数值模拟过程中，模型参数θ1、θ2和σε的Markov链的计算初始值分别取为8、5和0.5。θ1和θ2的建议分布采用正态分布，标准方差均为σ=0.1；σε的建议分布，按照其先验分布形式，假定它服从对数正态分布，取其标准方差σ′=1。Markov链总长为10 000次。获得的三个模型随机参数的后验样本的时间序列，见图1。由图1可知，在抽样过程中三个参数的后验样本很快进入稳定期。去掉初始抽样(燃烧期)1 000次后的后验分布的直方图见图2。由图2 可知，模型参数θ1、θ2的后验分布均表现出正态分布形式，σε的后验分布表现出对数正态分布形式。需要指出的是：对于马尔卡夫链燃烧期的选择，在学术上目前还没有一个统一的标准，本文算例燃烧期的选取根据时间序列图中所示的马尔卡夫链的收敛情况为初始抽样的10%；在M-H算法中，由于采用了随机游走抽样，所以对参数的建议分布的选取要求不是很严格，并不要求建议分布与其后验分布接近，而且计算结果收敛快。

图1 参数的后验样本时间序列

图2 模型参数θ1、θ2、σe的后验分布直方图

t∈(0.10,0.11]时间内疲劳裂纹损伤性能退化构件时变可靠度预测曲线见图3(a)。由图3(a)可知，t∈(0.10,0.11]时间内构件时变可靠度预测曲线随时间增长，下降非常缓慢，可靠性变化范围在1.0～0.999 9之间，这说明构件在t∈(0.10,0.11]内失效的可能性很小。图3(b)给出了t∈(0.11,0.12]时间内构件时变可靠度预测曲线的变化趋势见图3(b)。由图3(b)可知，该时变可靠度预测曲线与t∈(0.10,0.11]时间内的时变可靠度预测曲线的变化特征明显不同，下降速度非常快，当t=0.115左右时变可靠度降为0。如果取时变可靠度下降10%作为构件疲劳失效的阈值，那么从图中曲线可以看出，此时t=0.113 4，对应的可靠度为0.900 8，这与表2给出的第9条试件退化轨迹的疲劳失效时间t=0.113 0完全吻合，证实了本文方法的有效性。

图3 t∈(0.10,0.11]和t∈(0.11,0.12]时间内时变可靠度预测曲线

5.2 算例2

选择在疲劳试验终止时，裂纹长度尚未达到临界值的第18条试件的疲劳裂纹损伤性能退化轨迹作为预测疲劳损伤构件的疲劳裂纹增长的研究对象，构建疲劳裂纹增长概率模型随机参数的更新方法，动态预测疲劳裂纹损伤性能退化轨迹。疲劳裂纹损伤性能退化轨迹与时变可靠度预测曲线见图4。由图4(a)可知，当t=0.11时，疲劳裂纹长度为32.58 mm，实际量测得的疲劳裂纹长度为32.51 mm，绝对误差为0.07；当t=0.12时，疲劳裂纹长度为34.08 mm，实际量测得的疲劳裂纹长度为34.29 mm，绝对误差为0.21。证实了本文方法的准确性。由图4(b)可知，第18条疲劳裂纹在t=0.151 9时，疲劳裂纹长度自40.63 mm扩展至40.66 mm，达到临界值，此时试件疲劳失效。由图4(c)可知t∈(0.12，0.16]时间内的时变可靠度预测曲线与疲劳裂纹损伤性能退化轨迹预测结果吻合。

图4 疲劳裂纹损伤性能退化轨迹与时变可靠度预测曲线

图5 t∈(0.05,0.12]时间内疲劳裂纹损伤性能退化轨迹预测曲线

图5为t∈(0.05，0.12]时间内疲劳裂纹损伤性能退化轨迹预测曲线，其中第1次更新是以t=0.04时刻的已收集量测数据，相应的更新似然函数后得到的疲劳裂纹性能损伤退化轨迹曲线，第2次更新是以t=0.05时刻收集到的历史量测数据进行更新得到的退化轨迹曲线。

由图5可知，利用数据更新疲劳裂纹增长概率模型后得到的退化轨迹与实际量测得到的退化轨迹曲线相近，第1次更新后得到的退化轨迹与实际量测得到的退化轨迹的平均绝对误差为0.40 mm，均方根误差为0.44，平均绝对百分比误差为1.32%；第2次更新的退化轨迹更接近实际量测的退化轨迹，其平均绝对误差仅为0.30 mm，均方根误差为0.33，平均绝对百分比误差为0.97%。对比上述三个评价指标可知，第2次更新后得到的预测模型的平均绝对误差、均方根误差以及平均绝对百分比误差均小于第1次更新后得到的预测模型，可靠性更高。同时两次更新后得到的退化轨迹预测图相近，更进一步证实了贝叶斯更新理论在结构疲劳失效预测中的有效性和精确性，由此建立的钢构件性能退化损伤的贝叶斯动态预测模型更符合实际工程。

6 结论

本文给出了一个基于线弹性断裂力学和Paris半经验疲劳裂纹增长速率公式的概率模型，该模型可用于疲劳荷载作用下的桥梁结构中。在此基础上，进一步采用Bayesian更新理论和MCMC模拟技术，耦合实际结构新量测数据，建立了疲劳裂纹增长模型中随机参数的更新评价方法，进而动态预测疲劳裂纹损伤性能退化轨迹和结构时变可靠性。疲劳裂纹增长概率模型考虑了实际结构量测数据的系统误差和随机误差，它们在Bayesian更新过程中被处理为超级参数。通过实际应用算例研究，可以得到以下结论。

(1)本文提出的基于MCMC算法的Bayesian更新方法在进行模型参数修正的过程中充分考虑了模型参数的先验信息，并不断地利用新量测到的疲劳裂纹增长数据，有效地降低了抽样过程中不确定因素对抽样结果的影响，而且与仅考虑一次样本数据的传统统计学的计算结果相比，Bayesian更新统计学的计算结果更贴近于工程实际。

(2)用本文方法计算得到的可靠性预测曲线、疲劳裂纹损伤性能退化轨迹预测曲线与疲劳试验结果高度吻合，说明本文建议的方法能够有效预测反映结构从疲劳裂纹扩展开始至疲劳失效的全过程，适用于长期的、依时性动态过程预测且结果可靠，从而可为结构维修与养护方案、决策等提供科学依据。

(3)本文方法的准确性建立在合理的先验信息的基础上，因此对量测数据的精度和数量要求较高，此外，当考虑影响因素较多、样本数增加时，MCMC算法的计算效率较低，需要发展更高效率的抽样算法。

(4)文中的量测数据为实验条件所得，与实际工程中的桥梁结构疲劳裂纹增长有一定的区别，在未来的研究工作中，将更注重工程实际。通过分析实际监测数据的统计特征，建立考虑塑性影响的疲劳裂纹统计模型，结合数值方法建立有限元模型，进行钢结构疲劳寿命预测。