基于概率密度演化的斜拉索多级变幅时变疲劳可靠度分析

刘纲，唐伟，高凯

(1.重庆大学 山地城镇建设与新技术教育部重点实验室，重庆 400045；2.重庆大学 土木工程学院，重庆 400045)

斜拉桥是各国大跨桥梁常用桥型之一，斜拉索是斜拉桥的主要承力构件。在车辆、风等交变荷载作用下，斜拉索易出现疲劳损伤[1]，进而引起主梁破坏甚至垮塌事故，如2018年意大利莫兰迪高架桥，因其斜拉索疲劳损伤导致坍塌，直接经济损失多达上亿欧元，造成了多达43人死亡的重大事故。在疲劳损伤过程中，由于拉索材料、外部荷载与服役环境等因素影响，斜拉索破坏往往具有较强随机性[2-3]，为此，能较好应对这一问题的可靠度分析法已成为斜拉索疲劳损伤评估的一种常用方法[4]。疲劳可靠度分析常采用应力-强度干涉模型[2]，经典疲劳可靠度求解方法主要包括一次二阶矩法(FOSM)和蒙特卡洛法(MCS)[4-5]。RATHOD等[3]建立了考虑材料随机性概率损伤模型，并利用FOSM分析了常幅和3级变幅荷载作用下弹簧的疲劳可靠度，结果与试验值接近。XΙANG等[4-5]对比了MCS和FOSM在分析焊接构件疲劳可靠度方面的性能，结果表明在多级变幅作用下，MCS疲劳可靠度预测精度高于FOSM。以上方法能够评估材料或简单结构在常幅或荷载幅级变化不多情况下的疲劳可靠度，但对于随机参数较多的复杂构件，如斜拉索等构件，FOSM需对随机变量进行当量正态化处理，将大幅降低计算精度；MCS法虽精度较高，但其对计算性能要求较高，计算量大，较难应用于实际工程[2,6]。近些年，LΙ等[7-9]提出采用概率密度演化方法(PDEM)求解结构随机可靠度，该方法将概率守恒理论与数值分析方法相结合，将复杂的高维随机问题进行离散处理，在保证精度条件下，能大幅提高计算效率，已在结构动力可靠度、随机振动分析方面得到广泛应用。目前，已有学者将PDEM应用于材料和简单构件的疲劳可靠度分析。例如，XU[2]基于Miner准则构建材料疲劳损伤方程，利用PDEM实现了常幅和2级高-低变幅的材料疲劳可靠度求解，结果表明疲劳可靠度与实验值数据拟合较好且计算效率较高。对于3级以上的变幅加载，由于服役期内各级变幅荷载作用次数关于总荷载作用次数的导数无法确定，导致其偏微分方程无法求解。为规避此问题，FU等[6]基于PDEM静力可靠度求解思路分析了随机风荷载作用下风机塔筒的疲劳性能，得到了塔筒疲劳可靠度。但PDEM静力可靠度求解需对每个时间节点进行单独分析，计算效率将大幅降低。综上所述，PDEM已在材料和构件疲劳可靠度求解中得到初步应用，但在多级变幅荷载作用下疲劳时变可靠度求解时，还存在偏微分方程难确定和计算效率较低等问题。针对以上问题，基于线性疲劳累积损伤模型和PDEM算法，提出适用于多级变幅荷载作用下斜拉索疲劳时变可靠度分析方法，并通过LY12-CZ材料的5级变幅疲劳试验数据和某大桥斜拉索疲劳可靠度分析，验证所提方法的适用性。

1 斜拉索疲劳可靠度分析模型

实际工程中，构件往往受到多级变幅荷载作用，这也是材料或构件发生疲劳损伤的主要原因。而S-N曲线来源于常幅荷载试验，无法通过其直接获取变幅荷载作用下的疲劳可靠度分析模型。目前主要通过疲劳累积损伤准则来解决这一问题[10]，即：1)确定各级荷载作用幅值大小和作用次数，利用S-N曲线获得各作用幅值对应的疲劳寿命；2)确定变幅荷载作用下疲劳累积损伤准则，将各荷载幅值作用次数与其对应疲劳寿命按线性或非线性累积损伤原则进行求和，得到变幅荷载作用下疲劳损伤累积模型。

1.1 拉索疲劳寿命分布及累积损伤模型

斜拉索在服役期内始终处于高应力状态，拉应力均值将影响斜拉索的疲劳寿命，因此，SUH等[11]建议斜拉索S-N曲线取为：

式中：Nf为等效应力幅Δσeq作用下的疲劳寿命；m和C为材料参数；Δσ为应力幅；σb为材料强度极限；σm为平均应力。对式(1)两边取对数，可得：

研究表明，材料参数C服从对数正态分布[11-14]，其分布参数可通过常幅疲劳试验获得。参数m一般为常数，文献[10]建议取m=3.645。根据式(2)可得常幅荷载作用下斜拉索疲劳寿命与C同服从对数正态分布，其均值μ和方差σ分别为：

根据Miner准则，斜拉索的多级变幅疲劳损伤模型为：

式中：Da为累积疲劳损伤；ni和Nfi分别为第i级应力幅Δσi作用的次数和疲劳寿命；s为应力幅级数；n为各级应力幅的总加载次数。

1.2 斜拉索疲劳损伤极限状态函数

根据式(4)可知，斜拉索疲劳累积损伤为荷载作用次数n的单调函数。当累积损伤超过临界损伤时，则认为斜拉索失效，故其疲劳损伤极限状态函数可定义为：

式中：Y(n)为疲劳累积损伤极限状态；Dc为疲劳临界损伤，服从均值为1，标准差为0.3的对数正态分布[10,13]；F为描述随机参数与疲劳损伤状态的关系；Θ为随机参数，Θ=[(Nf1,Nf2,…,Nfs),Dc]。对式(5)关于n求导，可得斜拉索疲劳损伤状态函数的变化率Y˙(n)为：

2 基于PDEM的疲劳可靠度分析

2.1 多级变幅疲劳损伤的PDEM方程

XU等[2,6]已证明疲劳损伤过程为随机保守系统，即满足概率密度守恒原理。因此将式(5)和式(6)代入广义概率密度偏微分演化方程可得：

式中：p YΘ(y,θ,n)为[Y(n),Θ]的联合概率密度函数。

要求解偏微分方程(7)，须确定Y˙(n)的具体表达式，但各应力幅作用次数ni随服役时间不断变化，即∂ni/∂n的取值无法确定。因此，难以直接利用PDEM方法计算多级变幅斜拉索疲劳可靠度。

为解决该问题，笔者提出采用应力幅出现频率gi=ni/n来量化各级应力幅作用次数。考虑到斜拉索上作用的荷载具有一定周期性，可将作用次数等效为日均作用次数nd或年均作用次数ny。以nd为例，通过统计可得各应力幅日均出现频率gdi=ndi/nd，则在任意服役年限t内荷载作用次数n(t)=nd×365t，各应力幅作用次数为ni(t)=ndi×365t=gi(t)×n(t)，因此，gi(t)=ni(t)/n(t)，代入式(5)可得：

将式(8)代入式(6)可得：

进一步推导可知gi(t)关于n(t)的导数为：

利用链式求导法则可得：

将式(11)代入式(10)，由gdi×n(t)=ni(t)可得式(10)等于0，将结论代入式(9)可简化为：

由于极限状态函数是关于自变量时间t的函数，根据链式求导法则可将式(12)转换为关于时间t的变化率，即：

因此，通过gdi可确定Y˙(t)的具体表达式，即将其代入式(7)可得：

p YΘ(y,θ,t)初始条件为

式中：δ[·]为Dirac函数；pΘ(θ)为Θ的联合概率密度函数。将p YΘ(y,θ,t)在随机参数概率空间ΩΘ上积分可得Y(t)的时变概率密度函数：

对式(16)积分则可得疲劳时变可靠度：

式中：y≥0表示拉索累积疲劳损伤未超过临界损伤值；Pr{·}表示斜拉索可靠概率。

2.2 数值求解方法

基于PDEM的疲劳偏微分方程求解步骤如下：

1)通过雨流计数法获得各等效应力幅及作用次数；

2)根据各级等效应力幅日均作用次数ndi与日均荷载作用次数nd获得gi。再将各级等效应力幅代入式(3)得到Nfi的分布参数，即可得到疲劳损伤过程中随机参数向量Θ=[(Nf1,Nf2,…,Nfs),Dc]；

3)根据参数向量及其概率密度函数，在随机空间ΩΘ上获取均匀样本集Θq=[(Nf1,q,Nf2,q,…,Nfs,q),Dcq]，其中q=1,2,…,Nsel为样本数量，将样本代入式(8)和式(13)可得Y(t)和Y˙(t)；

4)将Y(t)和Y˙(t)分别代入式(13)和式(15)，利用总变差减小法(TVD)求解式(14)得到p YΘ(y,θ,t)；

5)通过式(16)对p YΘ(y,θ,t)在随机空间ΩΘ内积分，再利用式(17)即可获得多级变幅疲劳时变可靠度。

等效应力幅计算及偏微分方程求解原理及具体流程可参考文献[2,15,16]。

3 算例分析

3.1 材料LY12-CZ可靠度分析

为验证PDEM方法在材料多级变幅疲劳荷载作用下的可行性与计算效率，利用LY12-CZ的5级变幅循环加载试验数据进行验证，并将计算结果与MCS进行对比，参数见表1所示。具体试验数据和加载细节详见文献[17]。

表1 LY12-CZ材料应力幅分布参数Table 1 LY12-CZ material stress amplitude distribution parameters

利用2.2节计算步骤分析LY12-CZ时变疲劳可靠度，PDEM法抽取样数为159个，计算耗时26 s；MCS抽样数为105个，计算耗时156 s，约为PDEM法的6倍。2种方法得到的时变可靠度如图1所示。

图1 LY12-CZ的疲劳可靠度Fig.1 Fatigue reliability of LY12-CZ

从图1可知，PDEM法预测疲劳可靠度曲线与试验数据吻合良好；PDEM，MCS与试验数据的均方差分别为7.6×10-3和1.28×10-2，前者均方差小一个量级，这主要是因为PEDM选取样本在随机空间ΩΘ内均匀分布，即在失效区域与可靠区域均有样本，能保证前期小概率失效时，具有足够精度，如图1所示。因此，所提方法能够准确分析材料LY12-CZ在5级变幅循环加载作用下的疲劳可靠度。

3.2 斜拉索疲劳可靠度分析

3.2.1 斜拉索应力幅及随机参数确定

某双塔双索面斜拉桥跨径布置为40 m+110 m+320 m+110 m+40 m，塔高74 m，双主塔采用H型桥塔，主梁为C55预应力混凝土梁，如图2所示。全桥共布置4×25对斜拉索，每根拉索由公称直径为15.2 mm，抗拉强度为1 860 MPa的55根平行钢绞线组成。

图2 斜拉索大桥模型Fig.2 Cable-stayed bridge model

选取该桥具有典型代表的最短索B1，中长索B12和最长索B25进行疲劳可靠度分析。采用Midas建立该桥有限元模型，计算得到3根斜拉索对应的影响线如图3所示。由于缺乏本桥长期的车辆监测信息，为保证数据具有可靠性，将本桥短期监测数据与本桥处于同一线路上且相距最近的斜拉桥动态称重(WΙM)数据[18]进行对比，发现车辆信息基本一致。故将文献[18]中的动态称重(WΙM)数据作为参考。

图3 3根拉索的影响线Fig.3 Ιnfluence lines of the three cables

基于成熟的MATLAB车流模拟程序[10,13]生成各拉索的等效日均随机车流，结合斜拉索力影响线计算拉索应力时程，从而获得日均等效应力幅频率直方图，如图4所示。根据直方图可知索B1，B12，B25应力幅级数分别为10,18和16个。

图4 3根索的等效应力幅分布Fig.4 Equivalent stress amplitude distribution of three cables

参考文献[10,17]，选取斜拉索的随机参数如表2所示。

表2 斜拉索随机参数分布类型Table 2 Types of random parameter distribution of stay cables

3.2.2 结果分析

根据提出的疲劳可靠度计算方法和求解步骤，对拉索B1，B12和B25的疲劳可靠度进行分析，其中拉索B25的疲劳损伤演化过程如图5所示。

图5 B25斜拉索疲劳损伤概率密度演化曲面Fig.5 Evolution surface of fatigue damage probability density of B25 cable

将所提方法得到的疲劳可靠度与MCS进行对比，从图6可知，PDEM法与MCS法得到的时变疲劳可靠度曲线吻合较好。选取目标可靠指标β=3.5时[11,17]，最短索B1因所受应力幅较小，其疲劳寿命最长，约为204 a，与MCS计算值误差为2 a；中长索B12所受应力幅值最大，其疲劳寿命最短约为39 a，与MCS计算值误差为1 a；长索B25寿命约为120 a，与MCS计算值误差为1 a。对所选3根斜拉索，PDEM和MCS方法计算的疲劳可靠度预测误差均在2 a内，表明所提方法具有较高的计算精度。

图6 3根斜拉索的疲劳可靠度Fig.6 Fatigue reliability of three stay cables

在同一台电脑上，PDEM法与MCS法的计算时间对比如表3所示。PEDM方法计算时长最长约为85 s，而MCS法则高达978 s，是PDEM法的11倍。其原因在于MCS抽样数量是PDEM的1 000倍左右，故计算效率很低。此外，短索B1幅级数为10，耗时最短，而中长索B12幅级数为18，计算耗时最长。因此，随着幅级数不断增加，MCS计算耗时以及对电脑性能要求越高。

表3 MCS与PDEM计算耗时Table 3 MCSand PDEM calculation time-consuming

4 结论

1)通过提出的应力幅出现频率指标，解决了多级变幅累积损伤偏微分方程难确定问题，并提出了求解多级变幅疲劳可靠度的PDEM方法。

2)对于5幅级LY12-CZ材料疲劳试验，所提PDEM方法预测的疲劳时变可靠度与试验数据均方差为7.6×10-3，吻合度良好，且小于MCS所得均方差1.28×10-2，因为PDEM方法抽样样本均匀散布在状态空间内，在小概率失效时精度高于MCS。

3)对于斜拉索疲劳时变可靠度分析，本文所提方法与MCS预测结果基本一致。在给定可靠度指标β=3.5时，PDEM和MCS方法所得3根斜拉索的疲劳寿命误差均小于2 a。

4)所提PDEM方法样本抽样数量远小于MCS方法，故其计算效率约为MCS的8～11倍。