基于改进T-B法的宽带高斯过程疲劳损伤分析

2024-02-26袁奎霖孙卓成

船舶力学 2024年2期

袁奎霖，孙卓成

（大连理工大学工业装备结构分析国家重点实验室船舶工程学院，辽宁大连 116024）

0 引言

船舶与海洋工程结构物服役于恶劣的海洋环境中，长期遭受如风载荷、波浪载荷、海流载荷等多种随机交变载荷作用，极易产生疲劳损伤。在初期设计阶段，通常采用基于功率谱的频域方法计算随机应力响应造成的疲劳损伤。当结构应力是一个零均值窄带高斯随机过程时，其应力范围服从Ray⁃leigh 分布，Bendat 等[1]根据这种特性推导了在频域内的窄带疲劳损伤解析解。然而，实际情况下船海结构物在多种外部载荷和自身弹性变形影响下，结构应力往往是一个宽带高斯随机过程，其应力范围的概率分布难以确定。工程上往往仍采用基于窄带假设的Rayleigh 分布模型计算宽带疲劳损伤，但当带宽系数变大时将会高估疲劳损伤。因此，研究宽带高斯随机过程下疲劳损伤的频域评估方法具有一定的理论意义和实际应用价值。

目前，针对宽带高斯应力造成的疲劳损伤问题，学者们提出了一系列近似方法。Wirsching 和Light[2]假设宽带应力过程的应力范围仍服从Rayleigh 分布，在此基础上引入一个可以考虑带宽影响的修正系数对窄带近似疲劳损伤解进行修正。W-L 修正系数已经在船舶与海洋工程结构物疲劳规范[3-4]中得到广泛应用。由于宽带高斯应力的雨流幅值概率密度函数难以从应力响应功率谱直接推导，也有学者试图借助时域模拟仿真建立雨流幅值的经验概率密度函数，以此建立更为可靠的宽带疲劳损伤预报模型。Dirlik[5]采用不同参数形式的功率谱密度函数，通过蒙特卡罗时域模拟仿真对雨流计数后的应力范围进行统计分析，提出了一个雨流幅值的经验概率密度函数。Dirlik模型由一个指数分布和两个Rayleigh 分布组成描述雨流幅值分布，由此推导出的疲劳损伤解析解的计算精度较W-L修正系数法有明显改善。Zhao 和Baker[6]提出了由一个Weibull 分布和一个Rayleigh 分布线性组合而成的雨流幅值概率密度的近似模型，但该模型对于带宽系数较大的功率谱而言，计算效果仍有待改善。Tovo和Benasciutti[7-8]提出了T-B 法即将雨流计数损伤表示为由一个权重系数bTB控制的范围计数损伤和窄带假设损伤的线性组合，并通过大量数值实验得到了权重系数bTB的近似公式。已有研究[9-10]表明，T-B法具有适用范围广泛且计算精度较高的特点，对于不同谱型都具有较好的鲁棒性。

然而，T-B 法中参数bTB的近似公式仅与谱矩有关，并没有考虑S-N曲线斜率参数m的影响，因而随着斜率参数m的增大，T-B 法与雨流法的差异亦逐渐增大。针对该问题，本文在原有T-B 法基础上，在权重系数中引入S-N曲线斜率参数m加以修正，并通过大量数值实验结果拟合得到一个新的权重系数bMTB的非线性函数模型。在此基础上提出一种基于改进T-B法的宽带高斯疲劳损伤分析方法，并对该方法的有效性进行验证。

1 宽带高斯随机过程疲劳损伤分析方法

1.1 谱矩和带宽系数

对于单边功率谱密度函数为SXX(ω)的平稳高斯随机过程X(t) 而言，其谱矩定义为

式中，ω为角频率，单位为rad/s。对于高斯过程，其平均跨零率v0和平均峰值率vp可由谱矩表示：

此外，谱矩可以对随机过程的带宽进行表征，即带宽系数可表示为

其中，最常用的两个带宽系数α1和α2定义如下：

此外，工程上常用的与α1和α2有关的另外两种带宽系数定义如下：

式中，δ称为Vanmarcke 带宽系数[11]，取值范围为0 ≪δ≪1；ε称为Wirsching 带宽系数[2]，取值范围为0 ≪δ≪1。当α1和α2越趋近于0，随之δ和ε越趋近于1 时，该随机过程越趋近于理想宽带随机过程，反之则越趋近于窄带随机过程。工程上一般认为δ大于0.1时，可将一个随机过程视作宽带随机过程。

1.2 高斯随机疲劳损伤分析

目前工程上的疲劳损伤分析主要是基于S-N曲线方法进行的。材料的S-N曲线表示了该材料在恒幅循环载荷作用下应力范围与疲劳失效循环次数之间的关系，其一般表达式如下：

式中，S代表应力范围，N代表在特定的应力水平作用下疲劳失效时的循环次数，m和K分别是材料参数。

当结构遭受多级应力载荷时，工程上常采用Miner线性累积损伤理论计算总疲劳损伤：

式中，ni表示在应力范围Si作用下的循环次数，Ni表示在相同应力范围Si作用下疲劳失效时的循环次数，n代表载荷作用时长T内的总循环次数为Sm的数学期望。

当结构应力是一个随时间变化的随机过程时，其雨流循环的概率密度函数为fs(S)，则在随机响应作用时间T下的疲劳损伤可表示为

当随机过程X(t)是一个严格的零均值窄带高斯过程时，其峰值和谷值在每一个循环中对称出现，由此可以认为其雨流幅值分布与其峰值分布均服从Rayleigh分布：

其应力范围S是应力幅值R的2倍，也服从Rayleigh 分布，并且平均峰值率vp等于平均跨零率v0。由式（8）可得到作用时间T范围内的窄带疲劳损伤解析解为

式中，Γ(∙)表示gamma函数。目前,工程中常用Rayleigh分布模型来近似代替宽带过程，将一个零阶谱矩和平均跨零率与原始随机过程相等的窄带过程来近似代替实际宽带过程，即窄带近似方法。

1.3 宽带高斯疲劳损伤分析方法

对于宽带高斯随机过程，学者们根据不同功率谱生成大量随机时间历程，力争频域结果与雨流计数时域结果的误差最小化，提出了多种宽带高斯疲劳损伤的计算公式。

1.3.1 Wirsching-Light方法（W-L）

Wirsching 和Light[2]基于等效近似窄带过程的概念，提出了一种估算宽带高斯疲劳损伤的经验模型：

式中，DWL为W-L法估算的宽带高斯疲劳损伤，DNB为式（10）计算的窄带高斯疲劳损伤，ρWL是W-L法的修正因子，它是关于带宽系数ε和S-N曲线斜率参数m的函数，即

式中：

由上式可知，W-L方法认为宽带高斯过程的疲劳损伤仅与λ0、λ2和λ4三个谱矩参数和单斜率S-N曲线斜率参数m有关。

1.3.2 Dirlik方法（DK）

Dirlik[5]提出了雨流幅值分布的半经验公式，该分布由一个指数型分布和两个Rayleigh 型分布组成：

式中：

进而可以得到宽带高斯过程的疲劳损伤为

由上式可知，DK 方法的雨流幅值分布仅与两个带宽参数α1和α2有关，而α1和α2又与λ0、λ1、λ2和λ4四个谱矩参数有关。

1.3.3 Zhao-Baker方法（Z-B）

Zhao和Baker[6]提出了一个由Weibull分布和Rayleigh分布线性组合的雨流幅值分布，公式如下：

式中：

由式（8）和式（17）得到疲劳损伤公式如下：

1.3.4 Tovo-Benasciutti方法（T-B）

Tovo 和Benasciutti[7-8]通过理论分析发现，对于一个稳态高斯随机过程，其雨流计数疲劳损伤DRFC总是介于基于范围计数的疲劳损伤DRC和窄带假设疲劳损伤DNB之间，即

式中，DNB可根据公式（10）计算，DRC的近似解如下：

基于这种思想，Tovo 和Benasciutti 提出雨流计数损伤可表示为由一个权重系数bTB控制的范围计数损伤和窄带假设损伤的线性组合，即T-B模型：

为了正确评估DRFC需要确定权重系数bTB，Tovo 和Benasciutti 经过大量数值实验得到bTB的近似公式如下：

由此可知，T-B 法的权重参数bTB仅与带宽系数α1和α2即λ0、λ1、λ2和λ4四个谱矩参数有关，与S-N曲线斜率参数m无关。需要说明的是，Tovo和Benasciutti[7-8]通过数值实验确定bTB时将S-N曲线斜率参数固定为m=3。然而，已有研究[13-16]表明，随着S-N曲线斜率参数m的增大，T-B法计算的疲劳损伤与雨流法结果的误差明显增大。如图1所示，本文在前期研究中分析了m=3和m=5时，不同带宽（α1和α2）条件下由式（23）反推的bsim值与bTB近似公式即式（24）的差别。由图1 可知，当带宽系数α2=0.1时，bsim值受m值影响较小且与bTB近似公式吻合较好；随着α2的增大，m=3和m=5对应的bsim值之间的差别变得更加明显，表明T-B法权重参数bTB的近似公式应该考虑m的影响。

图1 权重系数的模拟结果bsim与拟合公式bTB的对比Fig.1 Comparison of the weighting factor bsim from numeri⁃cal simulation and approximate formula of bTB

2 基于改进T-B法的宽带高斯疲劳损伤分析方法

2.1 参数化功率谱

针对T-B法的权重系数提出新的拟合公式，本文采用逆傅里叶变换技术，根据不同参数形式的功率谱模拟生成大量的时间序列。采用时域的雨流计数法计算疲劳损伤，考虑不同的S-N曲线斜率参数m的影响，代入式（23）中反推出新的权重系数。为了能够得到更大范围的带宽参数，同时考虑到不同谱形对疲劳损伤的影响，本文采用了如图2所示的六种具有不同形状的功率谱。

图2 不同形状的功率谱Fig.2 Illustration of spectral shapes considered in this study

图2为六种参数化功率谱的示意图。以图2（a）中的对称二次型功率谱为例，其表达式如下：

图2（b）～（f）分别为反对称二次型、线性以及常数型功率谱，表达式可参照式（25）获得。图2（a）～（e）五种谱型的ω1和ω3为固定值，分别取2π/1000 rad/s和2π rad/s，而ω2介于ω1和ω3之间，通过改变ω2、h1和h2的值可以得到不同形状的功率谱。为了使得所有谱型的零阶谱矩λ0为常数，各参数有如下关系：

其中，图2（a）～（c）的M=3，图2（d）的M=2，图2（e）的M=1。

此外，图2（f）为分离式矩形双模态谱，常被用于高斯双模态疲劳损伤分析[12-13]的数值模拟中，其形状参数定义如下:

式中，B为高频模态与低频模态功率谱下的面积（能量）之比，R为两个模态的中心频率之比。为了保证两个矩形谱是分离的，R还需要满足如下条件：

本文中ωa取5 rad/s，B=0.01～10，R=3～20，A1+A2=1。c1和c2决定了低频和高频模态的带宽特性，本文中c1和c2都取为1.1。

2.2 改进的权重系数bMTB

如式（24）所示，原T-B 法的权重系数bTB的近似公式仅是与带宽系数α1和α2相关的函数。本文通过数值模拟方法重新建立权重参数b与带宽系数α1、α2以及S-N曲线斜率参数m的函数关系。船舶与海洋工程结构为钢质焊接结构，其S-N曲线一般为以N=107为拐点的双斜率曲线。但是，为了便于研究，本文将采用单斜率的S-N曲线进行疲劳损伤计算。具体方法为选取不同的α1、α2和m值，利用2.1节中的六种谱型生成大量随机时间序列进行数值模拟试验，利用雨流计数法计算疲劳损伤DsRimFC。由式（23）可得到通过数值模拟确定权重系数bsim的表达式为

式中，DNB和DRC可利用式（10）和式（22）直接计算。由此，生成了一系列与不同α1、α2和m对应的bsim值，通过非线性拟合技术确定新的权重系数函数模型bMTB：

式中，系数ai为随不同S-N曲线斜率参数m变化的待定系数。对于船舶与海洋工程领域常用的S-N曲线[4]，斜率参数m的范围为3～6，相应的拟合系数ai的具体值见表1。图3 以m=3 和m=5 为例，显示了数值模拟得到的权重系数bsim和两种拟合函数模型bTB和bMTB随α1和α2的变化规律。可以看出，新的权重系数模型bMTB与原T-B 法的权重系数模型bTB相比，与数值模拟结果的拟合效果更好，且能反映不同斜率参数m的影响。

表1 拟合系数ai(i = 1,2,3,4)Tab.1 Fitted coefficients ai(i = 1,2,3,4)

图3 权重系数b的数值模拟结果与拟合公式结果对比Fig.3 Comparison of the weighting factors of b from numerical simulation and empirical formulae

2.3 不同宽带高斯疲劳损伤评估方法的比较

为了将改进T-B 法与现有的宽带高斯疲劳损伤评估方法进行对比，本文分别采用决定系数R2和Tovo-Benasciutti[8]定义的误差指标EI作为衡量指标来评估各方法的准确性和鲁棒性。R2和EI的定义如下：

针对2.1 节中的六种谱型，S-N曲线斜率参数分别取m=3 和m=5 进行数值模拟试验，将改进T-B法与W-L 法[2]、DK 法[5]、Z-B 法[6]和T-B 法[7-8]的频域疲劳损伤结果与相应的时域雨流结果进行对比，各方法的决定系数R2和误差指标EI如表2和表3所示。

表2 不同频域方法的决定系数R2Tab.2 Determination coefficients R2 for different frequency-domain methods

表3 不同频域方法的误差指标EITab.3 Error indexes EI for different frequency-domain methods

由表2 可知，对于六种参数化功率谱的数值模拟结果而言，当m=3 时，改进T-B 法的决定系数R2的平均值为0.9996，与其他四种方法相比更接近于1；当m=5时，五种方法的决定系数R2均有不同程度的降低，而改进T-B 法的决定系数R2为0.9963，表明当斜率参数m增大时该方法仍能保持良好的计算精度。

由表3可知，当采用误差指标EI作为五种方法计算精度衡量指标时，改进T-B法的计算结果最为准确，其次为T-B法与DK法，且两者精度相当，而W-L法和Z-B法的误差最大。

为了更加直观地将改进T-B法与DK法和T-B法进行对比，采用式（36）定义的指标γ，分析三种方法的计算精度随Vanmarcke带宽系数（0.1 ≪δ≪0.95）的变化规律，如图4所示。

图4 DK法、T-B法与改进T-B法疲劳损伤评估结果对比Fig.4 Comparison of fatigue damage evaluated by DK,T-B and modified T-B methods

由图4 可知，对于m=3，DK 法、T-B 法和改进T-B 法三种方法得到的疲劳损伤与雨流损伤的比值均在1.0 附近，T-B 法的相对误差范围为-5%到5%，改进T-B 法的相对误差范围为-4%到4%；DK 法的损伤结果通常小于真实雨流损伤，并且当Vanmarcke 带宽系数大于0.7 时，最大误差超过-10%。由于疲劳损伤与应力循环的非线性关系会随着m的增大而增强，当m=5时，三种方法的误差都随着Van⁃marcke 带宽系数的递增而明显增大，改进T-B 法的最大误差为-21.19%，但仍然小于T-B 法的最大误差-25.25%和DK 法的最大误差-33.6%。结果表明，在不同参数化谱型、带宽范围以及S-N曲线斜率参数的条件下，本文所提出的新方法的计算精度较已有的宽带高斯疲劳损伤方法有了明显改善。

3 算例验证

本章选取船舶与海洋工程结构常遭受的一些典型宽带响应谱，对本文所提出的改进T-B 法的有效性进行进一步验证。

3.1 双模态功率谱

对于船海结构物而言，其结构应力响应在功率谱上常常呈现出两个明显的峰值[14]，即高斯双模态特征，因此高斯双模态随机过程是高斯宽带过程的一种特例。在此，讨论一个贴近于实际海洋结构物响应的双模态功率谱[15]，其表达式如下：

式中，A是用于调整谱方差为1 的比例因子，TW=2π/ωw表示海浪周期，ωN=2π/TN是结构物的一阶固有频率，ξ为阻尼比。图5 为双模态功率谱算例，其中TW取9 s，结构振动的一阶固有周期TN为2 s，阻尼比取0.04，该功率谱的带宽系数δ= 0.5496。根据图5 中的功率谱，模拟生成高斯随机过程时域信号，由雨流法得到的疲劳损伤作为参考值，对五种频域方法的计算精度进行对比，如表4所示。表中的相对误差定义为

表4 各方法计算结果相对于雨流法的相对误差（%，双模态功率谱时）Tab.4 Relative errors of different methods compared with RFC method(%)

图5 双模态功率谱数值模拟算例（TW=9 s，TN=2 s，x=0.04）Fig.5 Bimodal spectrum used in numerical simulation in case of TW=9 s,TN=2 s,x=0.04

从表4 中可以看出，在m=3～6 时，W-L 法、DK 法、Z-B 法和T-B 法的误差随着m的增大而变大；当m=6 时，W-L 法的误差为-16.61%，DK 法为-7.01%，Z-B 法为-7.99%，T-B 法为-13.63%。本文提出的改进T-B法在m=3～6时都较为准确，所有的误差都未超过4%。

3.2 三模态功率谱

结构响应功率谱上存在三个明显峰值的三模态随机过程在海洋工程领域也是一种典型的宽带高斯随机过程[16]。图6 为一个Spar 型漂浮式风机塔柱的弯矩响应谱[17]，三个模态的特征频率分别为0.18 rad/s、0.48 rad/s 和2.4 rad/s，分别对应于纵摇固有频率、波浪频率和塔柱振动的一阶固有频率，该功率谱的带宽系数δ= 0.6994。根据图6 中的功率谱，模拟生成高斯随机过程时域信号，由雨流法得到的疲劳损伤作为参考值，对五种频域方法的计算精度进行对比，如表5 所示。可以看出在该工况下，相对于雨流法结果，W-L法和Z-B 法的结果更加保守，且相对误差随着S-N曲线斜率参数m的增大而降低，最大误差分别为18.96%和14.40%；相比之下，DK 法、T-B 法和改进T-B法低估了疲劳损伤，且误差随S-N曲线斜率参数m的增大而增大，最大误差分别为-19.33%、-21.71%和-12.16%。

表5 各方法计算结果相对于雨流法的相对误差（%，三模态功率谱时）Tab.5 Relative errors of different methods compared with RFC method(%)

图6 漂浮式Spar型风机三模态弯矩响应谱[17]Fig.6 Tri-modal bending moment spectrum of a floating Spar wind turbine[17]

以上的验证算例表明，相对于已有的宽带高斯疲劳损伤评估方法，本文所提出的改进T-B法对于船海结构物领域常见的真实响应功率谱具有更好的计算精度。需要说明的是，采用W-L 方法和Z-B方法计算本节中的双模态和三模态功率谱算例时，随着Vanmarcke 带宽系数的不同，会存在高估或低估疲劳损伤的情况，在以往研究[9]中也有类似结论。