一種新的高斯多模態隨機疲勞損傷頻域分析方法

鄭向遠，高 山，李 煒

(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)此外，如何將模態耦合法拓展到非高斯的隨機疲勞分析，是未來的重要研究工作.