一种新的高斯多模态随机疲劳损伤频域分析方法

郑向远，高 山，李 炜

(1.清华大学深圳国际研究生院 海洋科学与技术学部，广东 深圳 518055；2.中国电建集团华东勘测设计研究院有限公司，杭州 311122； 3.大连理工大学 船舶工程学院，辽宁 大连 116024)

对于服役期间一直遭受循环荷载的结构，疲劳损伤破坏是其结构失效的一种主要模式.当可供分析的应力时程较长时，通过时域的雨流计数法[1]可以得到不同大小的应力循环的计数，进而依据Miner线性累积损伤理论[2]和S-N曲线[3]得到疲劳累积损伤.然而，在实际工程应用中，应力时程的获取一般基于有限元计算，当需要校核的工况繁多时，例如海上采油平台在设计时需要校核很多海况[4]，计算代价则异常庞大，仍然坚持使用时域分析方法不切实际，因此转而使用频域分析方法，通过结构物的响应功率谱来得到疲劳损伤[5-6]是更为实际的途径.

当结构响应的随机应力是一个窄带(narrow-banded, NB)高斯随机过程时，可以认为其雨流幅值的分布服从瑞利分布，疲劳损伤在频域内存在解析解.然而，当随机应力是一个宽带(wide-banded, WB)高斯随机过程时，其雨流幅值的概率分布尚不可推导.学者们因此提出了很多不同的经验型方法来得到宽带高斯随机过程的疲劳损伤，其中较为常用的有Single-moment (SM)方法[7]，Dirlik方法[8]，Zhao-Barker方法[9]，以及近年发展的Tovo-Benasciutti方法[10]，Park从雨流矩拟合的角度提出的JB法[11]等.

此外，结构物的应力响应在很多情形下会呈现出高斯多模态的特征.所谓的多模态特征是指响应在功率谱上呈现出多个分开的显著峰值，模态的多少由结构的频率响应函数和外部激励共同决定.例如，浮式平台的系泊缆索由于同时遭受波浪和平台运动的作用，其应力谱呈现出双模态特征[12]；海洋浮式风机在风浪耦合的联合作用下，其塔柱的弯矩谱会在浮式平台纵摇运动的频率、波浪频率和塔柱一阶固有频率出现3个峰值而呈现出三模态特征[13].

目前，高斯双模态随机过程的疲劳损伤分析已经积累了不少宝贵的研究成果.Jiao等[14]通过将双模态的疲劳损伤划分为大应力循环和小应力循环2个部分，提出了最早的理论框架(JM法).基于该理论框架，很多学者[15-17]对其中的应力循环计数、幅值概率分布等细节进行了修正和改良.Low[18]更是引入了相位角参数用于描述高频(high-frequency, HF)应力与低频(low-frequency, LF)应力间的相位差问题，大大提升了JM法的精度.除此之外，Han等[19]从正弦波叠加的角度出发，提出了一种新的疲劳损伤模型.

对于高斯三模态随机疲劳损伤的分析，Gao等[20]对JM法进行了拓展(GM法)，使其可以应用到3个频率模态的情形.Low[21]同样基于自己的双模态方法提出了三模态的疲劳损伤概率模型.Park等[22]也通过大量的数值分析证明JB法可以应用到三模态的情况.但是，这些现有的方法依然存在着一定的误差与不足[20]，对结构的疲劳寿命估计产生较大不确定性.

本研究致力于在频域内解决多模态的随机疲劳损伤的精确分析问题，从分割功率谱[23]的思路出发，提出了一种新的适用于双模态和三模态高斯随机过程的疲劳分析方法.该方法先是将功率谱分割成很多个极窄的频带，每个频带都可以认为是一个理想的窄带过程，因此，每个窄带过程造成的疲劳损伤可以由基于瑞利分布的窄带方法计算得到.而后，引入耦合因子ξ表述任意2个极窄频带间由于相互耦合影响而造成的疲劳损伤.最后，通过合理的非线性组合方式，将这些疲劳损伤组合在一起，得到蔚为完整的疲劳损伤.

1 频域疲劳分析理论

1.1 高斯窄带随机疲劳分析

对于一个严格的窄带高斯过程Y(t)，其峰值和谷值对称分布，且相继出现，因此可认为其雨流幅值分布与其峰值分布相同，均为瑞利分布.对于高斯随机过程，应力范围是应力幅值的两倍，服从瑞利分布，进而得到单位时间内的平均疲劳损伤为

(1)

式中：Γ(.)表示伽马函数；λ0表示0-阶谱矩，即方差.n-阶谱矩的定义为

(2)

当随机过程Y(t)不再是窄带时，由式(1)得到的疲劳损伤是不够准确的，因此，式(1)又称之为宽带高斯过程的窄带近似疲劳损伤.Vanmarcke带宽系数常用于描述随机过程的带宽性质[24]：

(3)

(4)

当随机过程趋向于理想窄带时，Vanmarcke带宽系数δ趋向于0.工程上，认为δ<0.1的随机过程可以近似作为窄带过程进行处理[13].

1.2 高斯双模态随机疲劳分析

工程中的大多数结构响应并不会像宽带白噪声一样，能量均匀地分布在一个较宽的频带里，而多是集中出现在分隔开的2个、3个或是多个频带之中.如果一个随机过程在功率谱上表现为2个分隔足够远的独立频带，称该随机过程为理想双模态随机过程.图1给出了一个高斯双模态矩形功率谱[10].可看到，随机过程Y(t)的功率谱中包含着2个分开的频带，分别对应于低频XLF(t)和高频XHF(t)，因此

Y(t)=XLF(t)+XHF(t).

(5)

图1 典型的矩形双模态谱

1.2.1 功率谱分割法

不同于从概率理论角度出发的传统方法，Benasciutti等[25]提出了一种从分割功率谱出发的频域分析方法.提出将功率谱分割成很多份非常窄的频带，继而每个频带都可认为是一个对应于频率为ωi的理想窄带高斯过程.这些分割得到的窄带高斯过程造成的疲劳损伤均可由式(1)得到：

(6)

式中：λ0,i即为第i个窄带的0阶谱矩；ν0,i表示第i个窄带对应的随机过程的平均过零率，由于每一份切割的窄带都可认为是一个理想的窄带高斯过程，所以ν0,i=ωi/2π.随后，Benasciutti等[25]又从多轴疲劳的损伤理论中获得启发，提出使用Projection-by-Projection (PbP)规则，对式(6)得到的窄带损伤进行组合.

(7)

式中num表示分割的频率窄带的份数，对式(7)展开得到：

(8)

注意到基于PbP规则的功率谱分割法即为基于λ2/k的单矩Single-moment (SM)方法，更多细节参考文献[25].

1.2.2 基于功率谱分割的高斯双模态耦合法

图2 双模态高斯随机疲劳分析中的模态耦合

事实上，基于PbP规则的功率谱分割法忽略了双模态中高频与低频间的相互作用[26]，这导致其在很多时候都无法准确描述双模态响应中“高频骑在低频”这一典型特征.基于此，本文提出一种新的基于功率谱分割的高斯双模态疲劳分析方法，旨在能准确考虑高频模态与低频模态之间的相互耦合作用.式(9)和图2分别是该新方法的表达式和图解.

(9)

(10)

(11)

(12)

就目前而言，ξLF&HF的具体表达式尚无法理论推导出.因此，本研究借助蒙特卡洛模拟得到了表征任意2个窄带模态间耦合程度ξ的经验表达式，它是2个频率模态中高频模态与低频模态的频率比值γ、高频模态与低频模态能量比值β以及S-N曲线中材料系数k的函数.本研究通过快速傅里叶变换(fast Fourier transform, FFT)，由蒙特卡洛模拟生成了大量的双模态高斯随机过程(低频与高频均是窄带)，其中，三参量的取值范围是γ=2,3,…,15、β=0.05,0.1,0.2,…,2.0、k=3,4,…,9.为确保生成的时间历程用于疲劳分析的可靠性，每段用于分析的双模态高斯时程包含至少107个低频应力循环，且每个高频应力循环中包含至少32个数据点.此外，时间历程的偏度控制在[-0.03, 0.03]，峰度控制在[2.9, 3.1]以确保高斯性[27].在由基于雨流计数法的时域方法得到这些随机过程的单位时间疲劳损伤后，代入式(9)～(12)便反推出时域下的ξ值.其中，拟合得到ξ的近似表达式为

ξ=[P1+P2ln(γ)+P3ln(β)+P4[ln(γ)]2+

P5[ln(β)]2+P6ln(γ)ln(β)]/

[1+P7ln(γ)+P8ln(β)+P9[ln(γ)]2+

P10[ln(β)]2+P11ln(γ)ln(β)].

(13)

式(13)形式的确定是先通过固定参数k，仅对γ和β进行拟合.在大量的不同的两参数函数中寻找出最佳的两参数形式后，再引入系数Pu=a0,u+a1,uk+a2,uk2(u=1,2,…,11)这一关于k的二次函数来计入k对于ξ的影响，最后由Levenberg-Marquardt算法对式(13)进行非线性拟合，得到系数a0,u、a1,u和a2,u，拟合结果见表1、2.

需要强调的是，正如前文中所述，用于确定式(13)中系数的数值试验均为低频与高频是窄带的情形，但是由于本方法的出发点是功率谱分割法，所以，式(9)并非仅仅适用于低频和高频都是窄带的双模态情形，还适用“宽带低频+窄带高频”的情形.

表1 参数Pu (2≤γ≤ 4)

表2 参数Pu(4≤γ≤15)

对于高频模态是宽带的情形，式(9)的应用需要一定的变化.当高频模态是宽带时，低频与高频的耦合现象将不再是“一对一”而是“一对多”，即在一个低频应力循环上会同时叠加若干个不同频率的高频应力循环.因此，需要先将宽带高频划分成M个子模态.随后，将低频模态和每个高频子模态都分割成num个极窄频带，再计算低频模态和每个高频子模态间的耦合.图3给出了高频是宽带时的模态耦合法的具体流程图.

图3 窄带低频+宽带高频的双模态高斯疲劳分析方法

此时，式(9)中的低频与高频的耦合项变成：

(14)

其中M即为划分的高频子模态个数.

(15)

(16)

在将高频模态划分成M个子模态时，本研究采用了等频率间隔(等频宽法)的划分方式.事实上，对一个很宽的频带进行划分的方式有很多种，例如本研究使用的等频率间隔划分，还有等能量划分，等带宽系数划分等[20- 21].通过大量的试算，得到使用等频率间隔的划分方式最适合于本研究提出的模态耦合法，使用等能量划分或是等带宽系数划分在精度上并没有提升.其中的原因，一是通过等频率间隔划分得到的子模态能更好地表征高频模态的频率跨度，二是等频率间隔划分相较于等能量划分和等带宽系数划分也更为容易理解和方便编程.

1.2.3 基于功率谱分割的高斯三模态耦合方法

上一节提出的耦合因子ξ反映了2个不同频率模态在疲劳分析中的耦合作用，因此，它也可以应用于高斯三模态的随机疲劳分析之中.

图4给出了一段高斯三模态随机过程Y(t)的时间历程，表示成低频、中频和高频3个高斯随机过程的叠加，即：

Y(t)=XLF(t)+XMF(t)+XHF(t).

(17)

图4 典型的三模态高斯随机过程

Y(t)的功率谱见图5，定义频率比和能量比：

(18)

(19)

图5 典型的矩形三模态谱

图6给出了三模态高斯随机疲劳分析中的模态耦合情况，相应的，总的疲劳损伤表达式为：

(20)

(21)

式(21)中二阶耦合项类似于双模态中的低频与高频的耦合[式(12)]，表达式为：

图6 三模态高斯随机疲劳分析中的模态耦合

(22)

(23)

(24)

(25)

以上的分析均是针对高频是窄带的情形，类似于双模态方法，对于高频模态是宽带的三模态高斯随机过程，应用式(20)、(21)进行疲劳分析时，也需要先将高频划分成M个子模态，此时式(21)中，

(26)

(27)

(28)

(29)

(30)

(31)

(32)

2 算例分析

2.1 高斯双模态随机应力疲劳

本小节先是对低频模态和高频模态都是窄带的情形进行了讨论，随后对低频模态和高频模态都是宽带的情形进行了讨论.事实上，正如前文中所述，本文提出的方法还适用于低频模态是宽带，高频模态是窄带的情形，预测精度亦非常准确.由于其在方法的使用上和低频模态和高频模态均是窄带的情形一致，因此，为节约篇幅，在本文中没有赘述.

2.1.1 窄带低频+窄带高频的高斯双模态疲劳

本小节以图1所示的双模态矩形功率谱为例，讨论在低频模态和高频模态都是窄带的情形(带宽δLF=δHF=0.057 6).以时域的雨流结果作为参考值，对比了LOW法、SM法以及本文提出的新方法，其中，LOW法是在近些年研究中最为准确的双模态方法[18, 21].其中，用于时域分析的时间历程生成方法与小节1.2.2中所叙相同.

表3给出了k=3，C=1，频率比γHL=2、6、12和能量比βHL=0.05、0.4、1.2、2的结果.k=3这一材料参数的取值在土木和海洋结构物中广泛使用.可以看到，SM法由于忽略了高频与低频间的相互作用，其在βHL=0.05时给出的疲劳损伤结果远较于其他两种方法误差要大得多.对于LOW法而言，在大多数情形下都能给出相当准确的计算结果，但是在γHL=2时，其误差明显.这是因为在γHL=2时，低频模态和高频模态距离很近，LOW法中采用的一些关于相位差的假设不尽合理.同时可以看到，本文提出的模态耦合法在所有情形下都能给出非常准确的疲劳估计，相对于雨流结果的误差都小于1%.

类似的，表4给出了k=6.5，C=1的结果，其中，k=6.5在汽车行业较为常用.当k=6.5时，应力循环与疲劳损伤间的联系存在很强的非线性，这时，SM法的误差明显增大，这说明SM法所丢失的耦合作用所占权重增大.同样的，LOW法的精度也有略微的下降，其在γHL=2，βHL=4时的误差接近20%.可喜的是，新方法依然拥有非常好的精度，误差不超过3%.

此外，本研究提出的模态耦合法在计算过程中仅涉及一维积分，因此其计算速度与JM法、SM法和TB法相仿，在一台普通的计算机下，其对于一种功率谱组合的计算时间不超过0.1 s.相对的，LOW法虽然也能给出较为准确的疲劳估计，但是其计算涉及变上限的三维积分(文献[18]中提及的通过使用级数展开对三重积分的简化仅在材料系数k为整数时可用)，对于一种功率谱组合的计算大约需要20 s.

表3 窄带低频与窄带高频矩形谱下不同方法相对于雨流结果的相对误差(k=3)

表4 窄带低频与窄带高频矩形谱下不同方法相对于雨流结果的相对误差(k=6.5)

2.1.2 宽带低频+宽带高频的高斯双模态疲劳

实际工程中，结构的双模态响应中，低频和高频经常可能出现是宽带的情形，因此本小节着重讨论讨论低频模态和高频模态都是宽带的情形.依然以图1所示的双模态矩形功率谱为例，但是带宽δLF=0.142 9，δHF=0.277 4.

表5给出了k=3，C=1，频率比γHL=3、6、12和能量比βHL=0.05、0.4、1.2、2的结果.可以看到，SM法的误差依然是3种方法中最大的，再次说明了在高斯双模态的随机疲劳分析中，高频与低频耦合影响不可忽略.还注意到，尽管LOW法是基于低频与高频都是窄带的假设推导的，但是此处对于宽带低频与宽带高频在k=3时，它仍然能给出较为准确的疲劳估计，其最大误差出现在γHL=12，βHL=0.4时，达到-8.49%.对于本文提出的模态耦合法，在对宽带高频划分成4个子模态后，可以看到，能得到非常准确的疲劳损伤结果，绝大多数误差都在3%以内，最大的误差出现在γHL= 12，βHL=1.2时，达到-5.24%.

表6给出了k=6.5，C=1的结果.由于k值的增大，疲劳损伤与应力循环关系的非线性增强，SM法和LOW法相对于雨流结果都出现了较大的偏离，最大误差分别为-36.07%和-29.18%.相对的，本文提出的模态耦合法依然能给出非常准确的疲劳损伤结果，绝大部分误差依然控制在5%以内，最大误差也不超过10%.其中，在βHL=2，γHL=12时，模态耦合法的误差达到了-9.85%，其可能的原因是这一情形下，将高频宽带划分成4个子模态(M=4)并非最优.对于M的具体取值，是本方法目前仍然需要深度研究的一个重要方向，也是未来的一大重要工作.

2.2 高斯三模态随机应力疲劳

本小节讨论新方法在应用于高斯三模态随机疲劳分析时的表现，并与工程中常用的SM法、Dirlik法等进行比较.讨论分为2个部分，第一部分以低频、中频以及高频均是窄带的理想的高斯三模态随机过程(矩形三模态谱)作为对象进行了讨论.随后，以一个漂浮式风力发电机在遭受到风浪耦合载荷时的桩柱弯矩作为对象，讨论了高频模态是宽带的三模态情形.

2.2.1 理想高斯三模态随机疲劳

对于低频、中频以及高频均是窄带的理想高斯三模态随机过程，本小节讨论了2个算例.

很多研究表明，TB法和Dirlik法精度相差无几[10,11]，且在三模态上，Dirlik法要更精确一些[22]，因此，为了简洁，本小节仅列出了新方法与SM法，Dirlik法的对比结果.此外，对于Case1这一理想三模态过程还给出了与JM法和LOW法的对比.

表7给出了Case1，k=3～6时，各方法相对于雨流结果的相对误差.可以看到GM法在k=3时的误差尚可，但是随着k值的增大，其误差迅速增大，k=6时，误差甚至超过了40%.SM法和Dirlik法的精度差不多，k=6时，其误差接近20%.LOW法和新提出的模态耦合法均能给出非常准确的疲劳计算结果，然而模态耦合法的计算量要少许多，且更加简单，易于编程.

表5 宽带低频与宽带高频矩形谱下不同方法相对于雨流结果的相对误差(k=3)

表6 宽带低频与宽带高频矩形谱下不同方法相对于雨流结果的相对误差(k=6.5)

表7 Case1中各方法的相对误差

表8给出了Case2，k=3～6时，各方法相对于雨流结果的相对误差.类似的，SM法和Dirlik法的精度差不多，k=6时，其误差在20%左右.相对的新提出的方法依然能给出非常准确的疲劳计算结果，k=6时的误差亦没有超过5%.因此，通过Case1和Case2这2个频率比不同的算例，可以证明新方法完全可以胜任这种低频、中频以及高频均是窄带的高斯三模态的随机疲劳分析问题.

表8 Case2中各方法的相对误差

2.2.2 漂浮式风机塔柱的高斯三模态疲劳

本小节以一个实际工程结构的三模态响应谱来讨论高频模态是宽带时，本文提出的模态耦合法的表现.图7是一个Spar型的NREL 5MW漂浮式风机[28]在遭受到风浪耦合荷载时其塔柱的弯矩功率谱[13]，可看到功率谱中有3个谱峰，呈现出明显的三模态特征.其中，低频模态的特征频率为0.18 rad/s，对应于该漂浮式结构的纵摇频率；中频模态的特征频率为0.48 rad/s，对应于波浪频率；高频模态的特征频率为2.4 rad/s，对应于塔柱振动的一阶固有频率.该弯矩响应的偏度是0.12，峰度是3.07，可以近似地认为这是一个高斯过程.因此，可以使用本文提出的方法对该响应进行疲劳分析.由于高频是一个宽带，故将其划分成了2个子模态，各模态的基本信息见表9.

图7 漂浮式Spar风机的三模态弯矩谱

表10给出了SM法、Dirlik法以及本文提出的模态耦合法得到的疲劳损伤，及其相对于雨流结果的误差.由表10中的结果可以看到Dirlik法精度较SM法略好，2个方法在k=3时，误差在10%左右；k=6时，Dirlik误差接近20%, SM法的误差超过25%.相对的，本文提出新方法给出的结果与雨流结果非常接近，k=3～5时，相对误差均在1%以内；k=6时，相对误差也仅有-1.03%，体现了其在工程应用中的巨大潜力.

表9 功率谱中各模态的基本信息

表10 各方法的相对误差

3 结 语

本文从功率谱分割的角度出发，提出了一种能应用于高斯双模态、三模态随机疲劳分析的频域分析方法.该方法最大的创新在于其先是通过分析典型的低频与高频均为窄带的双模态随机过程，得到了双模态响应中低频模态与高频模态间由于相互耦合作用而对疲劳损伤造成的影响.为此，提出了一个耦合因子ξ，在疲劳计算公式中引入交叉项，对这一影响定量化.事实上，这一耦合因子ξ可以用于表达高斯多模态随机疲劳中任意2个频率模态之间的耦合作用.基于此，本文进一步地将该双模态方法拓展到了高斯三模态过程的疲劳分析.

通过蒙特卡洛模拟生成众多不同的应力时程，并以时域雨流计数法结果作为参考，在双模态算例中，将本文提出的模态耦合法与LOW法、SM法进行了对比；在三模态算例中，新方法也与GM法、LOW法、SM法和Dirlik法进行了对比.高斯双模态的疲劳分析结果表明，模态耦合法在几乎所有算例下都胜于LOW法和SM法.特别的，针对于低频与高频都是宽带的算例讨论更是证明了耦合因子ξ虽然是基于窄带低频和高频的情形分析得到，其使用并没有受到必须是窄带这一限制.类似的，在三模态的疲劳分析中，模态耦合法不仅对于理想的三模态响应有着优秀的表现，在随后讨论的漂浮式风机算例中，其表现也一样令人满意.该法不仅精度高，而且编程容易，计算量少，明显优于诸多传统方法，在实际工程应用中具有巨大的潜力.

最后，不得不提的是：1)在处理高频是宽带的情形时，需要对高频模态划分成若干个子模态，通常2至4个子模态就可以获得较为精确的疲劳损伤结果，而计算量仅略微增加；2)虽然已有不少文献提出了针对高斯宽带疲劳计算的方法，但是这些方法尚不能精确处理本文所讨论的宽带多模态情形；3)此外，如何将模态耦合法拓展到非高斯的随机疲劳分析，是未来的重要研究工作.