顧秀聰， 王 軻， 李奕霖

(南京航空航天大學(機械結構力學及控制國家重點實驗室振動工程研究所)， 江蘇 南京 210016)

0 引 言

隨著現(xiàn)代工業(yè)技術的飛速發(fā)展， 飛行器所處的工作環(huán)境趨向高速、 高溫和高壓， 因此， 其所處的振動環(huán)境也日趨復雜， 由此造成的結構失效也越來越多， 因此振動疲勞壽命研究受到眾多學者和工程界的廣泛關注[1]， 工程實際中， 結構與構件失效很多都是由疲勞載荷引起的[2]. 據(jù)統(tǒng)計， 在機械零件失效中大約有70%屬于疲勞破壞， 尤其是在航空航天領域[3]. 目前， 進行疲勞壽命預估計算主要借助一些市面流行的疲勞壽命預估軟件， 例如MSC軟件中搭載的MSC.fatigue模塊， Ansys中的Ncode模塊等[4].

目前， 振動疲勞壽命預估主要分時域壽命預估方法和頻域壽命預估方法. 時域方法一般采用的方法是在結構危險位置的應力時域信號基礎上[5-7]， 采用循環(huán)計數(shù)方法將信號轉化為由許多個不同正弦波載荷所組成的加載序列， 對載荷進行雨流計數(shù)后， 對每一個循環(huán)的損傷進行計算， 并通過損傷累計準則計算對應的總損傷及壽命[8]， 使用較多的為miner線性損傷累積準則. 在整個計算流程中， 通過雨流計數(shù)法構建疲勞載荷譜是最為重要的環(huán)節(jié). 但是在處理長時間歷程、 大時域樣本數(shù)據(jù)的時域信號時， 常規(guī)雨流計數(shù)法處理速度緩慢， 處理時間長， 因此， 需要提高雨流計數(shù)處理的效率. 本文利用Matlab平臺， 以雨流并行處理技術為基礎， 編寫了時域振動疲勞壽命計算程序， 并通過航空典型結構九宮格加筋板驗證了程序計算的精確性.

1 振動疲勞壽命時域計算理論

1.1 雨流計數(shù)法

雨流計數(shù)法如圖 1 所示， 是以雙參數(shù)法為基礎的一種計數(shù)法， 考慮了動強度(幅值)和靜強度(均值)2個變量， 符合疲勞載荷固有特性， 其主要功能是把實測載荷歷程簡化為若干個載荷循環(huán)， 供疲勞壽命估算和編制疲勞試驗載荷譜使用. 1982年， 閻楚良院士在論文中通過用雨流計數(shù)法對隨機載荷時間歷程統(tǒng)計處理程序的研究， 系統(tǒng)地討論了隨機疲勞載荷的統(tǒng)計處理及其物理模型， 其中包括對隨機載荷時間歷程的壓縮處理、 峰谷值的檢測、 雨流計數(shù)法的圖解模型等， 成為雨流計數(shù)的基本理論. 現(xiàn)代對雨流計數(shù)的應用也多基于此理論[9].

圖 1 雨流計數(shù)法示意圖

對于振動疲勞問題， 工程實際中的振動信號載荷歷程普遍較長， 時域樣本規(guī)模較大， 特別是對于含有低頻分量的信號更是如此. 對于這類信號進行傳統(tǒng)的雨流計數(shù)處理， 計算效率低下， 計算時間過長. 因此， 亟需研究開發(fā)一種能夠處理長時間歷程振動信號的有效算法和處理方法.

1.2 累積損傷準則

疲勞累積損傷理論研究的是變幅疲勞載荷作用下疲勞損傷的累計規(guī)律和疲勞破壞準則. 累積損傷模型都認為， 在超出疲勞極限的應力水平下， 每一循環(huán)載荷都會給承載結構造成一定的損傷增量.

目前， 常用的疲勞累積理論主要有： 線性損傷累積理論、 修正的線性損傷累積理論、 非線性損傷累積理論以及概率累積損傷理論等.

線性疲勞損傷累積理論是指在循環(huán)載荷下， 疲勞損傷是可以線性累加的， 各個應力之間相互獨立和互不相關， 當累加的損傷到達某一臨界值時， 就認為試件發(fā)生疲勞破壞. 線性累積損傷理論中典型的是Miner理論. 該理論認為： 當?shù)趇級載荷循環(huán)了ni次時， 該級載荷對結構或試件造成的損傷

(1)

則所有各級載荷累計造成的損傷

(2)

那么疲勞壽命

(3)

式中：Tm為m級載荷所經(jīng)歷的總時間或循環(huán)數(shù)；DCR為損傷臨界值， 一般為1[10-11].

2 振動疲勞壽命時域預估流程

利用軟件進行時域疲勞壽命預估， 需要獲得模型的功率譜密度曲線并通過時域模擬獲得時域信號， 因此需要建立模型的有限元模型， 通過計算有限元的頻響， 結合隨機分析方法確定模型危險點并獲得危險點的PSD曲線. 利用三角級數(shù)法對PSD曲線進行時域模擬， 獲得時域信號. 對時域信號進行雨流計數(shù)處理， 獲得疲勞載荷譜， 再結合材料的S-N曲線和Miner線性損傷累計理論獲得結構的疲勞壽命. 計算流程如圖 2 所示.

圖 2 時域疲勞壽命預估流程圖

2.1 傳統(tǒng)雨流計數(shù)法的程序編制

傳統(tǒng)雨流計數(shù)法對首輪循環(huán)提取后剩余的發(fā)散收斂波是以半循環(huán)配對形式給出， 帶來的計數(shù)誤差較大. 本文采用波型封閉式計數(shù)模型[12]， 將載荷歷程從最高峰值處截開， 使左段的起點與右段的末點首尾相接， 構成收斂發(fā)散波， 再對重構后的收斂發(fā)散波進行雨流計數(shù)處理.

對相鄰的峰谷點XI,XI+1,XI+2,XI+3, 若XI

SA=|XI+2-XI+1|,

SM=(XI+2+XI+1)/2,N=1.(4)

若XI>XI+1,XI≥XI+2且XI+1≥XI+3， 則

SA=|XI+2-XI+1|,

SM=(XI+2+XI+1)/2,N=1,(5)

式中：SA為變程循環(huán)；SM為相對變程循環(huán)的平均值；N為變程循環(huán)次數(shù).

計完一次循環(huán)后， 從計算機內(nèi)存中抹去XI+1和XI+2兩點， 繼續(xù)判斷XI，XI+3，XI+4，XI+5.

若不符合上述循環(huán)提取條件， 則后移判斷XI+1，XI+2，XI+3，XI+4， 并XI將存放于Xs(t)數(shù)組中.

對剩余的峰谷值數(shù)列Xs(t)再進行循環(huán)提取， 直到所有循環(huán)提取結束.若循環(huán)提取結束時剩余的峰谷值點為Xs1,Xs2,Xs3， 根據(jù)波形封閉處理，Xs1與Xs3一定相等， 則直接提取此循環(huán)， 即SAs1=|Xs1-Xs2|.

綜上所述， 便得到了循環(huán)次數(shù)N、 循環(huán)幅值SA、 循環(huán)均值SM.

2.2 雨流計數(shù)法并行計算的程序編制

并行計算的編程實現(xiàn)可通過可進行并行編程預計算的軟件平臺實現(xiàn)， 在本文中初始編程采用的是matlab平臺， matlab平臺的并行語句為“parfor”語句. ”parfor”語句為并行計算的“for”語句， 但在循環(huán)中處理的數(shù)據(jù)不會出現(xiàn)數(shù)據(jù)的依賴， 例如循環(huán)內(nèi)部不能出現(xiàn)數(shù)據(jù)的傳遞等.

在進行雨流并行前， 需要先對原時域信號進行分段. 根據(jù)上文所述， “parfor”語句循環(huán)之間不能出現(xiàn)數(shù)據(jù)的依賴， 所以在分段點處會出現(xiàn)載荷循環(huán)的丟失， 當時域信號較長時， 丟失載荷循環(huán)的影響可以忽略不計. 同時， 為控制丟失載荷循環(huán)造成的影響， 需要對分段點進行選擇. 本文采用的準則是各分段點需要在等分點左右5%的范圍內(nèi)， 且應力絕對值最小.

分段后， 對各分段信號同時進行雨流計數(shù)處理， 得到各分段信號的載荷循環(huán)幅值SA1,SA2,SA3, …,SAn， 載荷循環(huán)均值SM1，SM2，SM3， …，SMn， 循環(huán)次數(shù)N1，N2，N3， …，Nn，n為分段數(shù). 再將上述分段信號的疲勞譜整合至一個疲勞譜. 利用Goodman修正法修正循環(huán)幅值SA， 再確定雨流分級數(shù)， 根據(jù)修正后疲勞譜中幅值SA的最大值和最小值進行分級， 得到分級疲勞譜.

并行雨流計數(shù)處理流程圖如圖 3 所示， 軟件操作頁面如圖 4 所示.

圖 3 并行雨流處理計算流程圖

圖 4 軟件操作頁面

3 疲勞壽命預估程序的驗證

本文采用的模型為九宮格加筋板模型. 該模型長 1 000 mm， 寬800 mm， 整體采用殼單元進行離散， 共10 020個有限元單元. 在有限元模型中， 板和加強筋均采用殼單元， 二者之間的連接方式采用RBE2的方式進行模擬. 為模擬試驗的約束和加載方式， 將模型四周與試驗中夾具接觸部分的所有節(jié)點通過RBE2的方式連接到同一節(jié)點上， 該節(jié)點除Z軸方向之外的其他自由度全部約束. 九宮格加筋板有限元模型見圖 5.

圖 5 九宮格加筋板有限元模型

該模型所采用的的材料為2024-O鋁合金， 其材料屬性見表 1.

表 1 2024-O鋁合金材料屬性表

本文采用的S-N曲線經(jīng)驗公式為冪函數(shù)形式

SmN=C,(6)

式中：m=7.423 9，C=3.873 6×1021.

根據(jù)疲勞試驗結果， 先通過MSC.Patran對該模型進行頻率范圍65 Hz～350 Hz下的頻率響應分析， 再輸入頻率范圍為65 Hz～350 Hz， 幅值大小為0.990 3 g2/Hz， RMS值為16.8 g的加速度平直譜， 對模型進行隨機響應分析， 由結果得出， 危險位置為九宮格板中心節(jié)點處， 危險位置的局部應力云圖如圖 6 所示.

圖 6 危險位置的局部應力云圖

由高應力點的應力RMS值的各分量可以看出， 該點x方向和y方向上的應力水平均較高， 因此采用Lemaitre多軸應力等效方法. 導出該點的x方向自功率譜密度曲線Gσxxσxx(f)、y方向自功率譜密度曲線Gσyyσyy(f)、xy方向自功率譜密度曲線Gτxyτxy(f)、x方向和y方向互功率譜密度曲線Gσxxσyy(f).

根據(jù)頻域上的Lemaitre等效公式[8]

2vGσxxσyy(f)+2(1+v)Gτxyτxy(f).(7)

計算得到響應的PSD曲線， 曲線如圖 7 所示.

圖 7 Lemaitre多軸等效PSD曲線圖

利用程序的時域模擬功能， 對該PSD曲線進行時域信號模擬， 模擬得到的時域信號時長為2 s， 時間步長為0.001 s， 時域信號如圖 8 所示.

圖 8 模擬后的時域信號圖

根據(jù)該時域信號進行疲勞壽命預估， 雨流分級數(shù)設置為100， 并行線程數(shù)設置為8. 經(jīng)過計算， 所得雨流計數(shù)結果如圖 9 所示.

圖 9 雨流計數(shù)結果圖

結合S-N曲線和miner線性損傷累計理論， 該結構的疲勞壽命為5 484.9 s.

將程序計算結果與試驗結果對比， 對比結構如表 2 所示.

表 2 疲勞壽命試驗與編制程序的時域壽命計算結果對比

4 結 論

本文研究了加入并行計算的振動疲勞壽命時域預估方法及其程序編制. 該程序通過導入時域振動信號， 利用雨流計數(shù)法編制了疲勞載荷譜， 再借助S-N曲線和損傷累積理論實現(xiàn)了時域疲勞壽命的計算. 同時針對工程實際中時域信號振動時間歷程較長、 數(shù)據(jù)量較大的問題， 程序中雨流計數(shù)過程加入了并行計算的功能， 通過選擇分段數(shù)和分段點控制誤差， 提高了計算效率. 程序先通過Matlab平臺進行初始編程， 然后利用執(zhí)行文件轉換Labview編程， 并編制操作頁面， 簡化操作， 并實現(xiàn)結果的可視化. 本文通過九宮格加筋板構件的實例， 分別計算了分段并行和未并行情況下的時域壽命， 證明了所編程序的精確性， 該疲勞壽命計算軟件可以運用到工程實際中.