穆童，孟鴿，謝里陽,*，張建波，石朝成

1. 東北大學(xué) 機械工程與自動化學(xué)院，沈陽 110819 2. 航空動力裝備振動及控制教育部重點實驗室，沈陽 110819 3. 上海航空材料結(jié)構(gòu)檢測股份有限公司，上海 201210

早在19世紀(jì)，人們就已經(jīng)認(rèn)識到了材料的疲勞問題，W?hler在1860年就繪制出了最早的S-N(應(yīng)力-壽命)曲線[1]。此后，對疲勞問題的研究逐步深入。目前，在常規(guī)疲勞問題的研究方面已經(jīng)形成了較為完善的理論體系。但是，振動疲勞的研究則起步較晚，直到1958年，美國學(xué)者Crandall將振動與疲勞問題結(jié)合起來[2]，從此人們才真正開始關(guān)注到振動疲勞問題。1963年，Crandall和Mark將振動疲勞定義為結(jié)構(gòu)在振動激勵下產(chǎn)生的響應(yīng)所造成的可以累加損傷的強度破壞[3]。

在疲勞壽命預(yù)測問題上，工程中常用的求解方式主要有2種，即時域分析和頻域分析。日常對各種物理量的認(rèn)知普遍是從時域角度出發(fā)的，表述時以時間作為橫坐標(biāo)，體現(xiàn)出了物理量在一定時間范圍內(nèi)的變化情況；而頻域方法表達的則是物理量與頻率之間的關(guān)系。相比之下，簡單的、規(guī)律性的變化用時域求解方式較為直觀、準(zhǔn)確，而當(dāng)面對一些復(fù)雜的或隨機的載荷時，頻域求解方式則能更加清晰明了地表達載荷特征，且求解速度更快。

在振動疲勞壽命預(yù)測理論模型方面，Bendat和Piersol于1966年率先提出了使用功率譜密度直接計算疲勞損傷的方法，其計算方法基于窄帶隨機過程的特點，使用功率譜密度函數(shù)估計隨機應(yīng)力的幅值分布，進而計算出結(jié)構(gòu)的疲勞損傷[4]；但是Bendat的研究是在窄帶隨機過程的基礎(chǔ)上進行的，而實際工程中的隨機振動普遍是以寬帶分布的形式存在的，若直接使用則會出現(xiàn)較大誤差。Wirsching等在Bendat的研究基礎(chǔ)上進行了修正，提出了新的計算模型[5-6]。此外，Chaudhury和Dover也基于Bendat所提出的計算模型進行了修正[7]。1985年，Dirlik使用了計算機模擬的方式對寬帶隨機過程的幅值分布進行了擬合，總結(jié)出了一個適用于寬帶分布的經(jīng)驗公式[8]，該公式對于帶寬的敏感性低，在實際使用時具有比其他模型更高的精度，因此應(yīng)用較為廣泛。Benasciutti和Tovo同樣在其研究中提出了一個適用于寬帶隨機過程的幅值分布模型，并且也具有較高的擬合精度[9]。

1 隨機振動疲勞理論

1.1 隨機過程

隨機過程可分為平穩(wěn)過程與非平穩(wěn)過程，生活中常見的一般為平穩(wěn)隨機過程，該條件下反映隨機過程x(t)在t1時刻自身前后相關(guān)性的自相關(guān)函數(shù)Rx不隨時間變化[10]，即

Rx(t1,t1+τ)=Rx(τ)=E[x(t)x(t+τ)]

(1)

式中：E為期望；t為時間；τ為時間間隔。

自相關(guān)函數(shù)在時域內(nèi)描述了隨機過程的統(tǒng)計特征，而在頻域內(nèi)描述隨機過程用的則是功率譜密度函數(shù)(Power Spectral Density，PSD)，二者組成一個傅里葉變換對[11-12]：

(2)

(3)

式中：Sx(ω)為功率譜密度函數(shù)；ω為圓頻率。

值得注意的是，功率譜的定義區(qū)間為正負(fù)2個實數(shù)域，稱作雙邊功率譜，而現(xiàn)實中不存在負(fù)頻率，因此定義單邊功率譜密度Gx(ω)為

(4)

此外，除了可以通過自相關(guān)函數(shù)定義隨機信號的功率譜密度外，還可以通過濾波法與快速傅里葉變換法來定義功率譜密度，公式分別為

(5)

(6)

式中：y(f,t,b)為濾波后中心頻率為f、時間長度為t、帶寬為b的信號的幅值；fft為對信號進行快速傅里葉變換。

快速傅里葉變換法即為將時間長度趨于無限的隨機信號通過傅里葉變換完成從時域到頻域的轉(zhuǎn)換。

基于信號的功率譜密度，定義平穩(wěn)隨機過程的譜矩mi為

(7)

其各階譜矩在時域內(nèi)的意義為[13]

m0=E[x2(t)]

(8)

(9)

(10)

定義了譜矩之后，可用各階譜矩來表達隨機過程的一些特征，如隨機過程的均方根值σRMS:

(11)

隨機過程的均方根值是均方值的開根，均方值代表了功率譜密度曲線與橫軸包圍的面積，面積越大，說明隨機譜的能量越大，因此隨機振動的均方值和均方根值在某種意義上表明了隨機振動的強弱，二者越大說明振動越強。

引入不規(guī)則因子γ來描述隨機過程的不規(guī)則性，其公式為

(12)

不規(guī)則因子的取值范圍為[0,1]，當(dāng)其取值趨近于1時，即意味著在時域歷程中幾乎每次正向穿越零水平后都會對應(yīng)出現(xiàn)一個波峰，也就意味著此時隨機過程為窄帶隨機過程，當(dāng)γ=1時，隨機過程則變成了簡單的簡諧波；相反地，當(dāng)不規(guī)則因子趨近于0時，隨機過程也隨之趨近于寬帶過程，當(dāng)γ=0時，隨機過程為白噪聲。

另外，帶寬系數(shù)ε也可以衡量隨機過程為寬帶還是窄帶過程：

(13)

與不規(guī)則因子相反，當(dāng)帶寬系數(shù)趨近于1時，隨機過程趨近于寬帶過程，反之則趨近于窄帶過程。

反映單位時間內(nèi)峰值次數(shù)的峰值期望率Vp[14]為

(14)

反映單位時間內(nèi)正向穿越均值的正斜率穿越期望值V+為

(15)

1.2 隨機疲勞

頻域疲勞壽命預(yù)測，廣泛應(yīng)用于解決結(jié)構(gòu)承受隨機載荷時的疲勞問題。經(jīng)典疲勞理論是建立在應(yīng)力或應(yīng)變的時域分析基礎(chǔ)之上的，主要思路為對零部件工作狀態(tài)中的危險點進行應(yīng)力應(yīng)變時間歷程的求解，結(jié)合Miner準(zhǔn)則與響應(yīng)歷程中的循環(huán)次數(shù)來計算累積損傷量，從而推算出預(yù)估的疲勞壽命；而頻域疲勞分析方法的基本思想則是從載荷的功率譜密度等統(tǒng)計特征出發(fā)，通過傳遞函數(shù)求解出零部件危險點的響應(yīng)應(yīng)力應(yīng)變的統(tǒng)計特征，并根據(jù)這些統(tǒng)計特征估算出各水平的應(yīng)力在單位時間內(nèi)的期望循環(huán)次數(shù)的概率密度函數(shù)，然后使用Miner準(zhǔn)則從統(tǒng)計學(xué)角度求出疲勞損傷的期望值，最終完成疲勞壽命的預(yù)測。時域方法能夠清晰明了地反映出結(jié)構(gòu)的壽命周期中的受載過程，可將加載順序等頻域方法無法考慮到的因素考慮進來，壽命預(yù)測更為精準(zhǔn)，但疲勞分析時域方法的計算量非常大，計算起來難度較大；而頻域分析方法的切入點則是載荷的統(tǒng)計學(xué)特征，該方法無法準(zhǔn)確地描述應(yīng)力應(yīng)變歷程，只能通過統(tǒng)計學(xué)的方式計算疲勞損傷，因此準(zhǔn)確性上要低于時域方法，但其優(yōu)勢在于計算量小，求解速度快，更加適合工程應(yīng)用。

根據(jù)Miner準(zhǔn)則可知，應(yīng)力應(yīng)變每循環(huán)一次都會造成一定損傷，不同大小的應(yīng)力應(yīng)變將造成不同程度的損傷，當(dāng)累積損傷量到達一定值時，零件將會發(fā)生疲勞失效[15]。S-N曲線中某應(yīng)力幅S對應(yīng)的疲勞循環(huán)次數(shù)為N(S)，該水平應(yīng)力在某時段內(nèi)實際的循環(huán)次數(shù)為n(S)，規(guī)定零件總損傷量為1時發(fā)生疲勞失效，則這一時段內(nèi)該水平應(yīng)力造成的損傷量D為[16-19]

(16)

式中：N(S)可以由材料疲勞曲線的表達式獲得，工程上常用的疲勞曲線冪函數(shù)公式為

SmN=C

(17)

其中：m與C均為材料常數(shù)；N為疲勞壽命對應(yīng)循環(huán)數(shù)。

若某工況下應(yīng)力為連續(xù)分布，則時間t內(nèi)幅值在(S,S+ΔS)的應(yīng)力循環(huán)次數(shù)可以記為

n(S)=vtp(S)ΔS

(18)

式中：v為單位時間內(nèi)所有應(yīng)力雨流幅值循環(huán)次數(shù)；p(S)代表幅值為S的應(yīng)力分布密度；ΔS為應(yīng)力區(qū)間。因此，連續(xù)分布的應(yīng)力幅造成的疲勞損傷為

(19)

則用時間表示的零件的疲勞壽命tL為

(20)

由此可見，要想在頻域范圍內(nèi)合理預(yù)測零部件的疲勞壽命，需要了解響應(yīng)應(yīng)力的幅值分布情況，即式(20)中的p(S)，而這正是頻域疲勞方法的重難點所在，目前的研究中最為常用的是Dirlik應(yīng)力幅值分布模型。

Dirlik使用蒙特卡洛方法在計算機中建立了時域模擬信號，對大量不同不規(guī)則因子的PSD曲線進行了模擬仿真，總結(jié)出了一套經(jīng)驗公式。Dirlik將幅值范圍概率密度函數(shù)假設(shè)為了一個指數(shù)分布與2個瑞利分布的和，公式為

(21)

式中：

Dirlik法是頻域疲勞壽命預(yù)測研究中使用最廣泛的方法，對于大部分寬帶隨機過程都有很好的適應(yīng)性，因此工程中一般采用該模型進行隨機振動疲勞壽命預(yù)測。雖然其只是經(jīng)驗公式，沒有理論支撐，但在絕大多數(shù)情況下，都有比其他模型更好的計算精度[12-14,16-20]。

以目前的研究與應(yīng)用來看，多數(shù)研究人員將頻域幅值分布模型直接應(yīng)用在了單一S-N曲線表達式所描述的應(yīng)力分布區(qū)間內(nèi)，但使用基于Dirlik等模型的頻域方法時，計算出的應(yīng)力幅分布區(qū)間寬泛，而且會有大量的小幅值應(yīng)力出現(xiàn)。對于一些常用的碳鋼、合金鋼、高強鋼等材料，其S-N曲線在雙對數(shù)坐標(biāo)系下會出現(xiàn)明顯的拐點，如果用單一的S-N曲線表達式來描述其疲勞性能將導(dǎo)致一部分應(yīng)力幅值對應(yīng)的損傷出現(xiàn)一定的誤差，降低了壽命預(yù)測模型的準(zhǔn)確度。因此在使用頻域幅值分布模型時，可以采用分段計算疲勞損傷的方法，這樣將得到更加精確的計算結(jié)果。

在考慮超高周疲勞的影響時，大部分材料的S-N曲線將呈現(xiàn)出三段甚至更多段的形式[21]，此時如果使用單一冪函數(shù)擬合的S-N曲線將會存在明顯的精度問題，顯然也應(yīng)該分開計算疲勞損傷。本文僅以兩段式S-N曲線進行說明。

使用雙對數(shù)坐標(biāo)系下的2條直線來擬合S-N曲線，其表達式為

(22)

式中：下標(biāo)a、b用以區(qū)分兩段表達式參數(shù);S0為兩段表達式的分隔點。

此時零件的疲勞損傷表達式為

(23)

對應(yīng)的疲勞壽命為

(24)

2 加速試驗

在機械裝備產(chǎn)品的研發(fā)過程中，可靠性試驗是不可或缺的一環(huán)。在產(chǎn)品前期的設(shè)計研發(fā)階段中，一般需要對產(chǎn)品的可靠性及壽命進行預(yù)估，并在之后通過可靠性試驗的方式對其進行驗證。通常來講，機械產(chǎn)品的疲勞壽命普遍較高，如果將其按照正常的工作應(yīng)力水平進行壽命試驗，需要很長的時間，因此可以在試驗時適當(dāng)提高激勵水平，使結(jié)構(gòu)工作在比正常工況更加嚴(yán)酷的環(huán)境中，達到更早地出現(xiàn)失效的目的，以此節(jié)省大量的試驗時間。在試驗結(jié)束后通過等效關(guān)系便可推算出正常應(yīng)力水平下的工作壽命，這一試驗方法即加速壽命試驗。

2.1 基于Dirlik模型的加速等效關(guān)系

除使用隨機響應(yīng)分析直接求解結(jié)構(gòu)的響應(yīng)功率譜密度外，工程上也會使用結(jié)構(gòu)的傳遞函數(shù)來得出危險點的應(yīng)力功率譜密度。在合理范圍內(nèi)，當(dāng)外載荷的施加形式、位置、方向沒有變化時，結(jié)構(gòu)的傳遞函數(shù)只取決于結(jié)構(gòu)本身，而與載荷大小無關(guān)。某點的響應(yīng)應(yīng)力PSD與在某特定位置處施加的激勵PSD之間的關(guān)系為[20]

(25)

式中：Gb(f)為結(jié)構(gòu)上b點的響應(yīng)應(yīng)力功率譜密度函數(shù)；Wa(f)為施加在結(jié)構(gòu)a處的激勵功率譜密度函數(shù)，其形式不限，可以為集中力、基礎(chǔ)運動等；Hba(f)則為結(jié)構(gòu)在a點處承受激勵時b點的頻響函數(shù)，即傳遞函數(shù)。

由式(25)可知，在一定合理范圍內(nèi)，當(dāng)結(jié)構(gòu)所受的激勵增大時，危險點處的響應(yīng)也隨之線性增大。當(dāng)激勵點處的激勵PSD的幅值W′a(f)變?yōu)樵蟹档膋倍時，危險點b處的響應(yīng)應(yīng)力PSD即G′(f)將同樣變?yōu)樵械膋倍，而根據(jù)式(7)中的譜矩定義，各階譜矩的變化情況為

(26)

隨機過程其他參數(shù)的變化情況為

(27)

(28)

(29)

(30)

(31)

式中：上標(biāo)′表示變化后的變量。對于Dirlik模型來說，當(dāng)響應(yīng)應(yīng)力PSD的幅值G′b(f)為原有幅值的k倍時，模型中各參數(shù)的變化情況為

D′3=1-D1-D2=D3

則此時應(yīng)力幅值的概率密度函數(shù)為

(32)

單位時間內(nèi)的疲勞損傷為

(33)

可推導(dǎo)出使用不完全gamma函數(shù)表達的基于Dirlik模型的損傷比例j損傷表達式為

(34)

因此壽命之間的比例關(guān)系j壽命為

(35)

其中:

式中：Γ(·)為上不完全伽馬函數(shù)；γ為下不完全伽馬函數(shù)。

2.2 算 例

某吊掛結(jié)構(gòu)局部如圖1所示，查詢機械工程材料數(shù)據(jù)手冊得知，該結(jié)構(gòu)所用材料的S-N曲線散點數(shù)據(jù)如表1所示[22]，可以看出該材料的S-N曲線在疲勞極限附近出現(xiàn)了明顯的拐點，如果用一條直線對其進行擬合則效果較差。忽略各修正系數(shù)的影響，使用雙對數(shù)坐標(biāo)系下的兩條直線對材料的疲勞特性擬合情況如圖2所示。

圖2中所對應(yīng)的兩擬合線的表達式為

對該結(jié)構(gòu)給定如圖3所示的y方向?qū)拵щS機加速度載荷，該載荷與GJB150A中所給出的噴氣式飛機外掛設(shè)備振動載荷相一致[23]。使用有限元分析得到危險點的位置如圖1所示，其響應(yīng)Mises應(yīng)力功率譜密度如圖4所示。

由圖4曲線可以看出，該結(jié)構(gòu)的響應(yīng)應(yīng)力PSD有多個尖峰，屬于寬帶隨機過程，使用Dirlik模型可以較好地擬合出應(yīng)力幅值的分布。基于該模型，使用式(24)求解出該結(jié)構(gòu)危險點的疲勞壽命為5.468×107s。

表1 材料S-N曲線散點[22]

圖2 S-N曲線擬合結(jié)果

圖3 輸入加速度功率譜密度

圖4 危險點響應(yīng)應(yīng)力功率譜密度

將該點的響應(yīng)PSD數(shù)據(jù)與材料疲勞性能數(shù)據(jù)代入式(35)中，使用MATLAB繪制出激勵與壽命放大倍數(shù)之間的加速關(guān)系如圖5所示，其中橫坐標(biāo)表示激勵PSD的放大倍數(shù)，縱坐標(biāo)表示壽命的放大倍數(shù)。從圖中可以看出，曲線大體呈兩段直線走向。

圖5 激勵與損傷的等效關(guān)系

不難理解這一規(guī)律，當(dāng)施加在結(jié)構(gòu)上的隨機載荷由小變大時，也就意味著有更多的大幅值應(yīng)力出現(xiàn)，即落入圖2所示S-N曲線左段的應(yīng)力幅占比越來越大，因此在激勵比例增加時，壽命將以更快的比例衰減，也可以說，圖5中的右半段直線對應(yīng)了所使用的S-N曲線高應(yīng)力部分，左半段則對應(yīng)了低應(yīng)力部分；而如果使用一條直線對S-N曲線進行表征，則這一等效關(guān)系也將變?yōu)橐粭l直線，即激勵與壽命呈對數(shù)線性關(guān)系，此種情況下與他人現(xiàn)有的研究結(jié)論保持一致[20]。Clevenson和Steiner曾經(jīng)對一批2024-T4鋁合金試件進行了多次隨機振動疲勞試驗，研究各種載荷參數(shù)對壽命的影響，針對2個不同形狀的PSD，在一定范圍內(nèi)改變一組試樣所受載荷大小，得到了不同均方根值下的壽命如表2[24]所示，將數(shù)據(jù)繪制在雙對數(shù)坐標(biāo)系內(nèi)如圖6所示。

圖6中，2條折線分別是根據(jù)2個PSD在各級載荷水平下的對數(shù)壽命均值點連線所得，曲線1與數(shù)據(jù)1對應(yīng)中表2中載荷均方根1列的數(shù)據(jù)，曲線2與數(shù)據(jù)2對應(yīng)表2中載荷均方根2列的數(shù)據(jù)?？梢宰⒁獾?組數(shù)據(jù)均呈現(xiàn)出較明顯的線性特征；同時，在2條折線的左上方均呈現(xiàn)出了一種向下彎折的趨向，同圖5中的曲線特征保持一致，從而證明了該方法的正確性。

表2 文獻[24]試驗數(shù)據(jù)

圖6 加速試驗結(jié)果

另外，一般加速試驗需要提高應(yīng)力水平來縮短試驗時間，即對應(yīng)圖中k>1的部分，可以看出，這一段曲線既有后續(xù)呈現(xiàn)線性的部分，也有兩段曲線的過渡部分，因此針對零部件結(jié)構(gòu)上這一點的加速關(guān)系不能用簡單的對數(shù)線性來表達。需要注意的是，假設(shè)該算例中其他條件均不變，增大原有載荷PSD的幅值，則須將曲線整體向左移，此時k>1的部分可按照對數(shù)線性近似處理；反之，需要將曲線向右平移，此時k>1的部分同時包含了曲線中兩段呈直線的部分，這樣則不能簡化處理。

由討論可以看出，在進行基于頻域疲勞模型的加速試驗時，加速關(guān)系不僅僅取決于材料性能，還與零部件所承受的載荷有關(guān)，如果按照一般的對數(shù)線性關(guān)系來計算的話，可能會出現(xiàn)較大的誤差。

3 結(jié) 論

1) 針對隨機振動中小應(yīng)力循環(huán)占比較大的情況，將Dirlik雨流幅值分布模型應(yīng)用到了雙對數(shù)坐標(biāo)系下的2段直線所表征的S-N曲線上，該求解思路使得頻域方法更加適用于隨機振動疲勞問題，且該方法還可以拓展到多段直線所表征的S-N曲線的情況中。

2) 基于所提出的多段式S-N曲線的求解思路，以2段的情況為例，推導(dǎo)了使用Dirlik模型為基礎(chǔ)的振動疲勞頻域求解方法支撐下的加速應(yīng)力等效關(guān)系式，總結(jié)出了加速壽命換算方法，為振動疲勞的加速試驗提供了理論支持。

3) 算例表明，在使用雙對數(shù)坐標(biāo)系下的2條直線對S-N曲線進行表達時，所推導(dǎo)出的應(yīng)力與損傷放大倍數(shù)的等效關(guān)系近似為2段對數(shù)直線關(guān)系，并借助相關(guān)文獻中的試驗數(shù)據(jù)對方法進行了驗證。