对于服役期间一直遭受循环荷载的结构，疲劳损伤破坏是其结构失效的一种主要模式.当可供分析的应力时程较长时，通过时域的雨流计数法^[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)，其峰值和谷值对称分布，且相继出现，因此可认为其雨流幅值分布与其峰值分布相同，均为瑞利分布.对于高斯随机过程，应力范围是应力幅值的两倍，服从瑞利分布，进而得到单位时间内的平均疲劳损伤为

$ {\bar D_{{\rm{NB}}}} = \frac{{{\nu _0}}}{C}{(2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sqrt {2{\lambda _0}} )^k}\Gamma \left( {1 + \frac{k}{2}} \right), $

(1)

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

$ {\lambda _n} = \int_0^\infty {{\omega ^n}} G(\omega ){\rm{d}}\omega = \int_0^\infty {{{(2\pi f)}^n}} G(f){\rm{d}}f, $

(2)

其中：G(ω)为功率谱密度；ω为角频率，rad/s；f为频率，Hz；${\nu _0} = {(2\pi )^{ - 1}}\sqrt {{\lambda _2} \cdot \lambda _0^{ - 1}} $表示平均过零率；C和k是S-N曲线^[3]中的材料参数，N=C·S^-k，表示在应力范围S这一水平下，材料将在遭受N次应力循环后发生疲劳破坏.

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

$ {\delta = \sqrt {1 - \alpha _1^2} ,} $

(3)

$ {{\alpha _1} = \frac{{{\lambda _1}}}{{\sqrt {{\lambda _0}{\lambda _2}} }}.} $

(4)

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

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

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

$ Y(t) = {X_{{\rm{LF}}}}(t) + {X_{{\rm{HF}}}}(t). $

(5)

图 1中，下标ub和lb分别表示一个频带的上边界和下边界，下标chr表示该频带的特征频率，对于矩形谱而言，特征频率即为该频带的中间值.此外，定义高频与低频间的频率比：γ^HL=ω_chr^HF/ω_chr^LF；高频模态与低频模态间的能量比(0阶谱矩之比)：β^HL=λ₀^HF/λ₀^LF.

1.2.1 功率谱分割法

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

$ {\bar d_i} = \frac{{{\nu _{0,i}}}}{C}{(2\sqrt {2{\lambda _{0,i}}} )^k}\Gamma (1 + \frac{k}{2}), $

(6)

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

$ {\bar D_{{\rm{Total}}}} = {\bar D_{{\rm{PbP}}}} = {(\sum\limits_{i = 1}^{{\rm{num}}} {\bar d_i^{2/k}} )^{k/2}}, $

(7)

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

$ \begin{array}{l} {{\bar D}_{{\rm{ Total }}}} = {\left\{ {\mathop \sum \limits_{i = 1}^{{\rm{ num }}} {{\left[ {\frac{{{\omega _i}}}{{2\pi C}}{{(2\sqrt {2G({\omega _i})\Delta \omega } )}^k}\Gamma \left( {1 + \frac{k}{2}} \right)} \right]}^{2/k}}} \right\}^{k/2}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\frac{{{{(2{\kern 1pt} {\kern 1pt} \sqrt 2 )}^k}}}{{2\pi C}}\Gamma \left( {1 + \frac{k}{2}} \right){{[\sum\limits_{i = 1}^{{\rm{ num }}} {\omega _i^{2/k}} G({\omega _i})\Delta \omega ]}^{k/2}} = }\\ {\frac{{{{(2{\kern 1pt} {\kern 1pt} \sqrt 2 )}^k}}}{{2\pi C}}\Gamma \left( {1 + \frac{k}{2}} \right){{({\lambda _{2/k}})}^{k/2}} \approx {{\bar D}_{{\rm{SM}}}},} \end{array} \end{array} $

(8)

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

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

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

$ {\bar D_{{\rm{ Total }}}} = \frac{{{{(2\sqrt 2 )}^k}}}{{2\pi C}}\Gamma \left( {1 + \frac{k}{2}} \right){(\lambda _{2/k}^{{\rm{LF}}} + \lambda _{2/k}^{{\rm{HF}}} + \lambda _{2/k}^{{\rm{LF\& HF}}}{\rm{ }})^{k/2}}. $

(9)

图 2中，低频模态和高频模态都被分割成了num份极窄的频带，每个频带对总疲劳损伤的贡献在式(9)中以λ_2/k^LF和λ_2/k^HF的形式呈现.

$ {\lambda _{2/k}^{{\rm{LF}}} = \sum\limits_{i = 1}^{{\rm{ num }}} {\rm{s}} {\rm{m}}_i^{{\rm{LF}}},\lambda _{2/k}^{{\rm{HF}}} = \sum\limits_{i = 1}^{{\rm{num}}} {\rm{s}} {\rm{m}}_i^{{\rm{HF}}},} $

(10)

$ {\left\{ {\begin{array}{*{20}{l}} {{\rm{sm}}_i^{{\rm{LF}}} = {{(\omega _i^{{\rm{LF}}})}^{2/k}}G(\omega _i^{{\rm{LF}}})\Delta {\omega ^{{\rm{LF}}}},}\\ {{\rm{sm}}_i^{{\rm{HF}}} = {{(\omega _i^{{\rm{HF}}})}^{2/k}}G(\omega _i^{{\rm{HF}}})\Delta {\omega ^{{\rm{HF}}}}.} \end{array}} \right.} $

(11)

此外，低频中的第i个频带(ω_i^LF)会和高频中的第i个频带(ω_i^HF)发生耦合，对这些耦合项加和即得到了式(9)中的交叉项λ_2/k^{LF & HF}.

$ \lambda _{2/k}^{{\rm{LF \& HF}}} = {\xi ^{{\rm{LF \& HF}}}}{\rm{ }}\sum\limits_{i = 1}^{{\rm{ num }}} {\sqrt {{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_i^{{\rm{HF}}}} } . $

(12)

在计算λ_2/k^{LF & HF}时，式(12)中包含了1个耦合因子ξ^{LF & HF}，它表示在双模态疲劳分析中低频模态和高频模态间的相互耦合程度，是高低频模态之间频率比γ^HL，能量比β^HL以及材料系数k的三参量函数.

就目前而言，ξ^{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.为确保生成的时间历程用于疲劳分析的可靠性，每段用于分析的双模态高斯时程包含至少10⁷个低频应力循环，且每个高频应力循环中包含至少32个数据点.此外，时间历程的偏度控制在[-0.03, 0.03]，峰度控制在[2.9, 3.1]以确保高斯性^[27].在由基于雨流计数法的时域方法得到这些随机过程的单位时间疲劳损伤后，代入式(9)~(12)便反推出时域下的ξ值.其中，拟合得到ξ的近似表达式为

$ \begin{array}{*{20}{l}} {\xi = [{P_1} + {P_2}{\rm{ln}}(\gamma ) + {P_3}{\rm{ln}}(\beta ) + {P_4}{{[{\rm{ln}}(\gamma )]}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {P_5}{{[{\rm{ln}}(\beta )]}^2} + {P_6}{\rm{ln}}(\gamma ){\rm{ln}}(\beta )]/}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [1 + {P_7}{\rm{ln}}(\gamma ) + {P_8}{\rm{ln}}(\beta ) + {P_9}{{[{\rm{ln}}(\gamma )]}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {P_{10}}{{[{\rm{ln}}(\beta )]}^2} + {P_{11}}{\rm{ln}}(\gamma ){\rm{ln}}(\beta )].} \end{array} $

(13)

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

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

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

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

$ \lambda _{2/k}^{{\rm{LF \& HF}}} = \sum\limits_{j = 1{\rm{ }}}^M {\lambda _{2/k,j}^{{\rm{LF \& HF}}}} , $

(14)

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

$ \lambda _{2/k,j}^{{\rm{LF \& HF}}} = \xi _j^{{\rm{LF \& HF}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt {{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_{ji}^{{\rm{HF}}}/M} } , $

(15)

$ {\rm{sm}}_{ji}^{{\rm{HF}}} = {(\omega _{ji}^{{\rm{HF}}})^{2/k}}G(\omega _{ji}^{{\rm{HF}}})\Delta {\omega ^{{\rm{HF}}}}. $

(16)

式(15)、(16)中，M表示高频模态中包含的子模态个数，ω_ji^HF表示其中第j个子模态所属的第i个极窄频带对应的角频率.式(15)左端λ_{2/k, j}^{LF & HF}表示第j个高频子模态与低频模态之间的疲劳损伤耦合影响，式(15)右端ξ_j^{LF & HF}的表示第j个高频子模态与低频模态的耦合系数，它是关于高频子模态与低频模态频率比γ_j^HL=ω_{j, chr}^HF/ω_chr^LF和能量比β_j^HL=(λ_{0, j}^HF)/λ₀^LF的函数，由式(13)计算得到.而λ₀^LF和λ_{0, j}^HF则分别是低频和第j个高频子模态的零阶谱矩.注意到，由于耦合中来自于低频模态的贡献不应该被多次计算，因此，式(15)右端根号里的矩积需要除以M.

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

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

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

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

$ Y(t) = {X_{{\rm{LF}}}}(t) + {X_{{\rm{MF}}}}(t) + {X_{{\rm{HF}}}}(t). $

(17)

图 4可知，低频应力作为载波出现，中频和高频叠加其上，所有的应力循环都同时受到低频、中频、高频的影响.因此，对于三模态高斯随机过程，其疲劳损伤来自于一阶项：低频λ_2/k^LF、中频λ_2/k^MF、高频λ_2/k^HF；二阶耦合项：低频与中频的耦合λ_2/k^{LF & MF}、低频与高频的耦合λ_2/k^{LF & HF}、中频与高频的耦合λ_2/k^{MF & HF}；三阶耦合项：低频、中频、高频之间的耦合λ_2/k^{LF & MF & HF}.

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

$ \begin{array}{*{20}{c}} {{\gamma ^{{\rm{HL}}}} = \omega _{{\rm{chr}}}^{{\rm{HF}}}/\omega _{{\rm{chr}}}^{{\rm{LF}}};{\gamma ^{{\rm{HM}}}} = \omega _{{\rm{chr}}}^{{\rm{HF}}}/\omega _{{\rm{chr}}}^{{\rm{MF}}};}\\ {{\gamma ^{{\rm{ML}}}} = \omega _{{\rm{chr}}}^{{\rm{MF}}}/\omega _{{\rm{chr}}}^{{\rm{LF}}},} \end{array} $

(18)

$ \begin{array}{*{20}{c}} {{\beta ^{{\rm{HL}}}} = \lambda _0^{{\rm{HF}}}/\lambda _0^{{\rm{LF}}};{\beta ^{{\rm{HM}}}} = \lambda _0^{{\rm{HF}}}/\lambda _0^{{\rm{MF}}};}\\ {{\beta ^{{\rm{ML}}}} = \lambda _0^{{\rm{MF}}}/\lambda _0^{{\rm{LF}}}.} \end{array} $

(19)

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

$ {\bar D_{{\rm{Total}}}} = \frac{{{{(2{\kern 1pt} {\kern 1pt} \sqrt 2 )}^k}}}{{2\pi C}}\Gamma \left( {1 + \frac{k}{2}} \right) \cdot {(\lambda _{2/k}^{{\rm{Total}}})^{k/2}}, $

(20)

$ \begin{array}{l} \lambda _{2/k}^{{\rm{ Total }}} = (\lambda _{2/k}^{{\rm{LF}}} + \lambda _{2/k}^{{\rm{MF}}} + \lambda _{2/k}^{{\rm{HF}}}) + \\ (\lambda _{2/k}^{{\rm{LF\& MF}}} + \lambda _{2/k}^{{\rm{LF\& HF}}} + \lambda _{2/k}^{{\rm{MF\& HF}}}) + \lambda _{2/k}^{{\rm{LF\& MF\& HF}}}. \end{array} $

(21)

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

$ \left\{ \begin{array}{l} \lambda _{2/k}^{{\rm{LF\& MF}}} = {\xi ^{{\rm{LF\& MF}}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt {{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_i^{{\rm{MF}}}} } ,\\ \lambda _{2/k}^{{\rm{LF\& MF}}} = {\xi ^{{\rm{LF\& MF}}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt {{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_i^{{\rm{MF}}}} } ,\\ \lambda _{2/k}^{{\rm{LF\& HF}}} = {\xi ^{{\rm{LF\& HF}}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt {{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_i^{{\rm{HF}}}} } . \end{array} \right. $

(22)

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{sm}}_i^{{\rm{LF}}} = {{(\omega _i^{{\rm{LF}}})}^{2/K}}G(\omega _i^{{\rm{LF}}})\Delta {\omega ^{{\rm{LF}}}},}\\ {{\rm{sm}}_i^{{\rm{MF}}} = {{(\omega _i^{{\rm{MF}}})}^{2/K}}G(\omega _i^{{\rm{MF}}})\Delta {\omega ^{{\rm{MF}}}},}\\ {{\rm{sm}}_i^{{\rm{HF}}} = {{(\omega _i^{{\rm{HF}}})}^{2/K}}G(\omega _i^{{\rm{HF}}})\Delta {\omega ^{{\rm{HF}}}}.} \end{array}} \right. $

(23)

式中耦合因子ξ^{LF & MF}、ξ^{MF & HF}以及ξ^{LF & HF}均可由式(13)得到.对于低频、中频以及高频的三阶耦合项λ_2/k^{LF & MF & HF}，受高阶波浪摄动理论的启发，并经过大量的尝试，得到表达式：

$ \lambda _{2/k}^{{\rm{LF\& MF\& HF}}} = {\xi ^{{\rm{LF\& MF\& HF}}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt[3]{{{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_i^{{\rm{MF}}}{\rm{sm}}_i^{{\rm{HF}}}}}} , $

(24)

$ {\xi ^{{\rm{LF\& MF\& HF}}}} = \sqrt[3]{{\left| {{\xi ^{{\rm{LF\& MF}}}}{\xi ^{{\rm{LF\& HF}}}}{\xi ^{{\rm{MF\& HF}}}}} \right|}}. $

(25)

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

$ \lambda _{2/k}^{{\rm{LF\& HF}}} = \sum\limits_{j = 1}^M {\lambda _{2/k,j}^{{\rm{LF\& HF}}}} , $

(26)

$ \lambda _{2/k}^{{\rm{MF\& HF}}} = \sum\limits_{j = 1}^M {\lambda _{2/k,j}^{{\rm{MF\& HF}}}} . $

(27)

其中λ_{2/k, j}^{LF & HF}和λ_{2/k, j}^{MF & HF}分别表示低频模态和中频模态与第j个高频子模态的耦合：

$ \lambda _{2/k,j}^{{\rm{LF\& HF}}} = \xi _j^{{\rm{LF\& HF}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt {{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_{ji}^{{\rm{HF}}}/M} } , $

(28)

$ \lambda _{2/k,j}^{{\rm{MF\& HF}}} = \xi _j^{{\rm{MF\& HF}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt {{\rm{sm}}_i^{{\rm{MF}}}{\rm{sm}}_{ji}^{{\rm{HF}}}/M} } , $

(29)

$ {\rm{sm}}_{ji}^{{\rm{HF}}} = {(\omega _{ji}^{{\rm{HF}}})^{2/k}}G(\omega _{ji}^{{\rm{HF}}})\Delta {\omega ^{{\rm{HF}}}}. $

(30)

式(28)中，ξ_j^{LF & HF}和ξ_j^{MF & HF}由式(13)计算得到，前者是γ_j^HL=ω_{j, chr}^HF/ω_chr^LF和β_j^HL=(λ_{0, j}^HF)/λ₀^LF的函数，后者是γ₂^HM=ω_{j, chr}^HF/ω_chr^MF和β_j^HM=(λ_{0, j}^HF)/λ₀^MF的函数.λ_{0, j}^HF表示第j个高频子模态的0阶谱矩.类似的，三阶耦合项λ_2/k^{LF & MF & HF}改写为

$ \left\{ \begin{array}{l} \lambda _{2/k}^{{\rm{LF\& MF\& HF}}} = \sum\limits_{j = 1}^M {\lambda _{2/k,j}^{{\rm{LF\& MF\& HF}}}} ,\\ \lambda _{2/k,j}^{{\rm{LF\& MF\& HF}}} = \xi _j^{{\rm{LF\& MF\& HF}}}\sum\limits_{i = 1}^{{\rm{num}}} {\sqrt[3]{{{\rm{sm}}_i^{{\rm{LF}}}{\rm{sm}}_i^{{\rm{MF}}}{\rm{sm}}_{ji}^{{\rm{HF}}}/M}}} . \end{array} \right. $

(31)

其中λ_{2/k, j}^{LF & MF & HF}表示低频模态、中频模态与第j个高频子模态的耦合，其耦合因子为

$ \xi _j^{{\rm{LF\& MF\& HF}}} = \sqrt[3]{{\left| {\xi _{}^{{\rm{LF\& MF}}} + \xi _j^{{\rm{MF\& HF}}} + \xi _j^{{\rm{LF\& HF}}}} \right|}}. $

(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.

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个算例.

Case1：频率比ω_chr^LF:ω_chr^MF:ω_chr^HF=1:3:9, 3个频率模态的能量相等，且带宽系数δ均为0.05；

Case2：频率比ω_chr^LF:ω_chr^MF:ω_chr^HF=1:3:6，3个频率模态的能量相等，且带宽系数δ均为0.05.

很多研究表明，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法和新提出的模态耦合法均能给出非常准确的疲劳计算结果，然而模态耦合法的计算量要少许多，且更加简单，易于编程.

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

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.

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

3 结语

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

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

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