滕懷亮，李本威，韋 祥，張 赟，王永華

(海軍航空大學(xué) 航空基礎(chǔ)學(xué)院，山東 煙臺 264001)

載荷譜是航空發(fā)動機整機及零部件的應(yīng)力分析、壽命監(jiān)控、結(jié)構(gòu)設(shè)計以及可靠性分析的基礎(chǔ)條件。美軍的B-1轟炸機、F-16戰(zhàn)斗機等發(fā)動機均出現(xiàn)了耐久性問題，其原因主要是與原設(shè)計任務(wù)譜相比，發(fā)動機的實際使用載荷譜條件更為嚴(yán)苛，循環(huán)情況更復(fù)雜[1-3]。進(jìn)行發(fā)動機零部件試驗時，根據(jù)已有的飛行數(shù)據(jù)編制載荷譜是需要首先解決的具有挑戰(zhàn)性的工作[4]。國內(nèi)宋迎東[5-6]，趙福星[7]等對發(fā)動機載荷譜的相關(guān)工作進(jìn)行了深入研究，但渦軸發(fā)動機的載荷使用特點與渦噴、渦扇發(fā)動機相比有較大不同[8]，在編制載荷譜時應(yīng)予以特別關(guān)注。

編制渦軸發(fā)動機載荷譜的一個主要問題在于確定其主要工作狀態(tài)及持續(xù)時間。顧永根[9]，楊興宇[10]等在劃分發(fā)動機狀態(tài)時采取確定某參數(shù)保持平穩(wěn)的時間段，并將其作為一個狀態(tài)的方法。宋迎東等[11]采用平滑處理和多參數(shù)同時穩(wěn)定的方式識別發(fā)動機載荷狀態(tài)。上述文獻(xiàn)中確定不同發(fā)動機不同參數(shù)穩(wěn)定狀態(tài)范圍的上下限主要依賴經(jīng)驗，或直接采用設(shè)計狀態(tài)，而渦軸發(fā)動機的實際載荷與設(shè)計狀態(tài)相比有很大差異。本文依據(jù)雨流計數(shù)低幅值轉(zhuǎn)速循環(huán)的分布規(guī)律，利用聚類的方法確定了發(fā)動機實際使用過程中的主要工作狀態(tài)。

編制載荷譜的另一個主要問題是處理眾多的發(fā)動機載荷循環(huán)類型，一般采用壓縮等效的方法，將眾多的循環(huán)類型歸類為幾個典型的循環(huán)類型[10,12]。壓縮等效時需要依據(jù)根據(jù)疲勞、蠕變和熱沖擊損傷一致的原則[8]，保證載荷譜和實際飛行任務(wù)譜的一致性。史海秋等[13]研究了某型發(fā)動機承力機匣疲勞載荷譜的編制方法，但只分析了與轉(zhuǎn)速成線性關(guān)系的應(yīng)力循環(huán)的等壽命線形式，沒有給出通用的等壽命帶劃分方法。本研究推導(dǎo)了等壽命線的轉(zhuǎn)速循環(huán)峰谷值表達(dá)形式，給出了基于等壽命線和發(fā)動機實際主要工作狀態(tài)的距離判別方法，將發(fā)動機的眾多循環(huán)狀態(tài)歸類為11種主要循環(huán)。

1 渦軸發(fā)動機主要工作狀態(tài)

該型直升機飛行過程中多數(shù)不會用到發(fā)動機的設(shè)計轉(zhuǎn)速狀態(tài)ng,d，因此多數(shù)飛行剖面的主循環(huán)峰值轉(zhuǎn)速nmp小于設(shè)計轉(zhuǎn)速，且發(fā)動機實際使用的主要工作狀態(tài)普遍小于設(shè)計狀態(tài)，需要獲取發(fā)動機實際使用過程中的主要載荷狀態(tài)及循環(huán)情況。根據(jù)疲勞理論，包含相同的循環(huán)數(shù)的兩個載荷譜對零部件造成的疲勞損傷基本相同，所以要提取實際飛行任務(wù)數(shù)據(jù)的載荷循環(huán)。小幅值循環(huán)反映了轉(zhuǎn)速在某一狀態(tài)附近的小幅波動[14-15]，其密集分布的區(qū)域?qū)嶋H上反映出發(fā)動機轉(zhuǎn)速主要工作在以該類循環(huán)的均值為代表的工作狀態(tài)，因而可以利用小幅值循環(huán)的分布規(guī)律確定發(fā)動機的主要工作狀態(tài)。

1.1 雨流計數(shù)

循環(huán)計數(shù)的方法主要分為單參數(shù)計數(shù)法和雙參數(shù)計數(shù)法，雙參數(shù)計數(shù)法中的雨流計數(shù)因其與應(yīng)力變化過程一致而在工程中得到了廣泛應(yīng)用[14]。如圖1所示，引用文獻(xiàn)[14]中的例子對雨流計數(shù)法進(jìn)行檢驗和說明，其中1、3、5、7、9為谷值點，2、4、6、8為峰值點，雨流計數(shù)規(guī)則可以總結(jié)為雨流的起點和終點兩個部分：

1) 起點：雨流從各峰谷值1-9點內(nèi)側(cè)開始向下流；

2) 終點：從峰值開始的雨流到峰值更大的點對應(yīng)的時間即停止；從谷值開始的雨流到谷值更小的點對應(yīng)的時間即停止；遇到上面落下的雨流即停止。

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

保留起點時間作為循環(huán)開始時間，終點時間與起點時間之差作為循環(huán)持續(xù)時間。相鄰兩個雨流若峰谷值相同則構(gòu)成一個全循環(huán)，如2-3-2’，若不相同，則自成一個半循環(huán)，如1-2-4。經(jīng)驗證，程序統(tǒng)計結(jié)果與手工統(tǒng)計結(jié)果一致，得到全循環(huán)：2-3-2’，5-6-5’，7-8-7’和半循環(huán)：1-2-4，4-7-9。選取某次飛行的轉(zhuǎn)速-時間剖面，對其進(jìn)行雨流計數(shù)后，共得到1861個循環(huán)，并獲得了每個循環(huán)的幅值、均值、循環(huán)開始時間和循環(huán)持續(xù)時間、循環(huán)類型(半循環(huán)：用0.5表示或全循環(huán)：用1表示)，部分循環(huán)信息如表1所示。

表1 雨流計數(shù)結(jié)果

統(tǒng)計所有循環(huán)均值、幅值的不同區(qū)間的循環(huán)數(shù)，如圖2所示，高均值、低幅值的循環(huán)數(shù)最多，即主要工作轉(zhuǎn)速區(qū)間的小幅值循環(huán)較多；高幅值循環(huán)數(shù)較少，即對零部件造成疲勞損傷較大的循環(huán)較少。將各個循環(huán)以“谷峰值對”的形式作為橫縱坐標(biāo)在圖中表示出來，如圖3所示。從圖3可以發(fā)現(xiàn)，實際載荷循環(huán)包含大量接近零幅值的小幅值循環(huán)，且分布呈現(xiàn)圍繞幾個主要工作狀態(tài)聚集的簇狀分布。

圖2 雨流矩陣示意圖

圖3 轉(zhuǎn)速循環(huán)在峰谷值等壽命線中的分布圖

1.2 主要工作狀態(tài)識別

對于轉(zhuǎn)動件，影響蠕變的主要參數(shù)為轉(zhuǎn)速。轉(zhuǎn)速保持相對穩(wěn)定的階段即為發(fā)動機的穩(wěn)定狀態(tài)，使用頻次大的轉(zhuǎn)速狀態(tài)則為發(fā)動機的主要工作狀態(tài)。

如圖4所示，黑色表示載荷變化過程，彩色表示雨流計數(shù)的循環(huán)情況，該過程中有一個大幅值循環(huán)，其余均為圍繞某一均值的小幅值循環(huán)。小幅值轉(zhuǎn)速循環(huán)實際為轉(zhuǎn)速較為穩(wěn)定狀態(tài)下的循環(huán)，所以可以通過對小幅值循環(huán)進(jìn)行聚類進(jìn)而確定發(fā)動機的主要工作狀態(tài)。k-means聚類算法具有速度快和可擴展性良好等特點，在數(shù)據(jù)處理等領(lǐng)域得到了廣泛的應(yīng)用[16]。進(jìn)行聚類要選定分類數(shù)據(jù)集的屬性，根據(jù)小幅值循環(huán)的峰谷值分布情況，采用峰谷值對作為分類屬性，但k均值聚類算法的劣勢在于需要預(yù)先確定分類數(shù)目，在不知道最佳分類數(shù)的情況下，可以采用肘部法確定最佳分類數(shù)[17]，首先定義一個數(shù)據(jù)對象和其所在的聚類域的目標(biāo)函數(shù)，然后通過不同分類數(shù)目對應(yīng)的目標(biāo)函數(shù)值的大小評價聚類效果。常用的目標(biāo)函數(shù)是中心誤差平方和J。

其中：xi為第j類數(shù)據(jù)中的第i個數(shù)據(jù)向量；μj表示第j類數(shù)據(jù)的中心向量；向量xi的維數(shù)為數(shù)據(jù)的屬性值個數(shù)，在本例中即為峰值和谷值兩個；Nj為第j類數(shù)據(jù)中的數(shù)據(jù)個數(shù)；k表示數(shù)據(jù)分類的數(shù)目。

計算不同分類數(shù)的目標(biāo)函數(shù)值即平均畸變程度J值，其變化規(guī)律如圖5所示。選取肘部位置的分類數(shù)5作為最佳分類數(shù)，得到目標(biāo)數(shù)據(jù)的分類情況及各類別中心如圖6所示。

圖4 穩(wěn)定狀態(tài)示意圖

圖5 不同聚類數(shù)的平均畸變程度曲線

圖6 數(shù)據(jù)分類情況及聚類中心圖

根據(jù)聚類中心峰谷值可以得到峰谷值均值對應(yīng)的主要狀態(tài)，如表2所示，其中實際使用主要狀態(tài)1(0.591 4)和狀態(tài)2(0.759 0)與發(fā)動機設(shè)計地慢(0.598 8)和設(shè)計空慢(0.748 5)狀態(tài)非常接近。但使用時間最長的巡航狀態(tài)(0.827 3、0.852 1)和最大狀態(tài)(0.901 9)均遠(yuǎn)小于發(fā)動機的設(shè)計狀態(tài)(即起飛狀態(tài))(1.000 0)。

表2 主要狀態(tài)的時間分配情況

這是因為該型直升機發(fā)動機三發(fā)總功率比直升機所需功率大很多，所以執(zhí)行訓(xùn)練任務(wù)時，發(fā)動機的實際使用狀態(tài)要小于設(shè)計使用狀態(tài)，單純使用巡航、最大連續(xù)、起飛等設(shè)計狀態(tài)進(jìn)行分析會與實際使用載荷狀態(tài)有較大出入。依據(jù)主要載荷狀態(tài)，可以將發(fā)動機主要循環(huán)類型分為表3所示的11類。

表3 主要循環(huán)類型

2 循環(huán)類型壓縮

編制發(fā)動機載荷譜的另一個主要問題是處理類型眾多的循環(huán)。單次飛行的轉(zhuǎn)速數(shù)據(jù)，忽略雨流計數(shù)的小幅值循環(huán)后，仍然有幾十個循環(huán)類型，在工程應(yīng)用中對所有載荷剖面進(jìn)行零部件和壽命分析是非常困難的，所以需要進(jìn)行循環(huán)類型壓縮。壓縮時需要遵循疲勞、蠕變及熱沖擊損傷等效的原則，蠕變和熱沖擊損傷等效時需要保證載荷狀態(tài)及循環(huán)類型較為接近。疲勞損傷的等效需要依據(jù)等壽命線，將非脈動循環(huán)轉(zhuǎn)化為脈動循環(huán)。為在歸類轉(zhuǎn)速循環(huán)時將與轉(zhuǎn)速幅值有關(guān)的疲勞損傷和與發(fā)動機狀態(tài)有關(guān)的熱沖擊損傷和蠕變損傷統(tǒng)一起來，需要推導(dǎo)等壽命線的轉(zhuǎn)速表達(dá)形式。

2.1 等壽命線相關(guān)參數(shù)標(biāo)準(zhǔn)化

為在統(tǒng)一標(biāo)準(zhǔn)下進(jìn)行比較，需將燃?xì)鉁u輪轉(zhuǎn)速ng換算至標(biāo)準(zhǔn)大氣條件下的換算轉(zhuǎn)速ng,c，并利用該型渦軸發(fā)動機燃?xì)鉁u輪最大轉(zhuǎn)速將其化為標(biāo)準(zhǔn)化轉(zhuǎn)速n，n∈[0,1]。

如圖7所示，某轉(zhuǎn)速循環(huán)的最小值(即谷值)為nmin，最大值(即峰值)為nmax。應(yīng)力幅σa，平均應(yīng)力σm如式(1)、式(2)所示

(1)

(2)

圖7 正弦轉(zhuǎn)速循環(huán)與應(yīng)力循環(huán)曲線

根據(jù)發(fā)動機旋轉(zhuǎn)件的局部應(yīng)力-應(yīng)變特點，離心應(yīng)力與轉(zhuǎn)速的平方成線性關(guān)系：

σ=σmv+a·n2

(3)

式中：σmv為0轉(zhuǎn)速時的殘余應(yīng)力；a為轉(zhuǎn)速和應(yīng)力之間的比例系數(shù)；n為轉(zhuǎn)速。

若該型發(fā)動機最大轉(zhuǎn)速ng,d對應(yīng)的應(yīng)力為σmp時，則標(biāo)準(zhǔn)化后的應(yīng)力循環(huán)對應(yīng)的均值σm、幅值σa及強度極限σb分別為

式中，σb為材料的強度極限，根據(jù)需要分析的旋轉(zhuǎn)件材料進(jìn)行選取。

如圖8所示，等壽命線的表達(dá)式變?yōu)槭?4)

(4)

圖8 等壽命線

2.2 等壽命線轉(zhuǎn)換

定義轉(zhuǎn)速的幅值函數(shù)Na和均值函數(shù)Nm如式(5)、式(6)所示：

(5)

(6)

則由式(1)—式(3)可得

代入式(4)，可得

(7)

令

式(7)變?yōu)槭?8)

Na=k(Nm-b)

(8)

一般取b=1.2，即等壽命線輻射點為P(1.2,0)，由式(8)可知圖8中的等壽命曲線123就轉(zhuǎn)化為圖9中PB、PC、PD。

圖9 轉(zhuǎn)速循環(huán)等壽命線

其中，OA線表示谷值為0的脈動循環(huán)，AF線表示峰值為最大轉(zhuǎn)速nmp的循環(huán)，所以O(shè)AF三角形內(nèi)側(cè)表示旋轉(zhuǎn)機械工作區(qū)域。ABCD分別對應(yīng)峰值轉(zhuǎn)速為nmp、0.9nmp、0.8nmp、0.7nmp，谷值為0的脈動循環(huán)。ABCD對應(yīng)的轉(zhuǎn)速幅值函數(shù)Na的值可以由式(7)計算得到，進(jìn)而求得不同壽命的等壽命線所對應(yīng)的斜率kPA,kPB,kPC,kPD。

由Na=k(Nm-b)可以得到峰谷值之間的函數(shù)關(guān)系式

(9)

其中ki為不同等壽命線的斜率，i=1,2,…,10時分別對應(yīng)轉(zhuǎn)速谷值為0、轉(zhuǎn)速峰值為0.1nmp、0.2nmp、…、0.9nmp的脈動循環(huán)所在等壽命線，如k9=kPB,k10=kPA。可以得到以峰谷值表示的等壽命線如圖10所示，3條紅色邊界線nmax=1，nmin=0和nmax=nmin分別表示峰值轉(zhuǎn)速為最大轉(zhuǎn)速的循環(huán)、谷值轉(zhuǎn)速為最小轉(zhuǎn)速的循環(huán)和循環(huán)幅值為0的轉(zhuǎn)速循環(huán)?？梢钥闯觯茖?dǎo)出的等壽命循環(huán)的峰谷值關(guān)系滿足公式(9)，該公式的優(yōu)勢在于在考慮疲勞損傷的同時，沒有忽略高均值低幅值循環(huán)造成的與狀態(tài)變化相關(guān)的熱沖擊及蠕變等損傷，可以在同一標(biāo)準(zhǔn)下根據(jù)主要循環(huán)狀態(tài)將造成熱沖擊損傷較大、疲勞損傷較小的高均值低幅值循環(huán)考慮進(jìn)來，使得載荷譜更接近實際載荷狀態(tài)。

2.3 判別歸類

圖10 峰谷值表示的等壽命曲線

3 編制模擬任務(wù)試車譜

根據(jù)飛行數(shù)據(jù)編制模擬任務(wù)試車譜的流程如圖11所示。雨流計數(shù)得到載荷循環(huán)信息，對小幅值循環(huán)進(jìn)行k-means聚類，確定主要狀態(tài)和主要循環(huán)類型，依據(jù)等壽命線和主要循環(huán)類型，通過距離判別壓縮大幅值循環(huán)類型，得到模擬任務(wù)載荷譜。

3.1 數(shù)據(jù)預(yù)處理

飛參數(shù)據(jù)中存在很多等值點，將轉(zhuǎn)速大小相同的相鄰各點合并為1個點，該轉(zhuǎn)速點對應(yīng)的時間為各相鄰等轉(zhuǎn)速點的時間均值。同時，需要將超出轉(zhuǎn)速范圍及傳感器誤采集等原因造成的不符合實際情況的數(shù)據(jù)剔除，或?qū)ζ溥M(jìn)行平滑處理。雨流計數(shù)前需要先提取峰谷值，判斷第i個轉(zhuǎn)速點ng(i)是否為峰谷值，即判斷相鄰兩點與該點的差值ng(i)-ng(i-1)和ng(i)-ng(i+1)是否同號[18]：若兩轉(zhuǎn)速差異號，則ng(i)非峰谷值；若兩轉(zhuǎn)速差均為正，則ng(i)為峰值，若均為負(fù)，則ng(i)為谷值。

圖11 編制載荷譜流程框圖

3.2 生成載荷譜

根據(jù)疲勞理論，小于疲勞極限的載荷循環(huán)造成的損傷可以忽略，所以只需保留造成疲勞損傷較大的大幅值循環(huán)。疲勞極限可以作為舍棄小幅值循環(huán)，保留大幅值循環(huán)的臨界點。在具體分析時應(yīng)根據(jù)分析對象的不同，利用零部件材料S-N曲線確定應(yīng)力疲勞極限，計算其對應(yīng)的轉(zhuǎn)速。本文取疲勞極限對應(yīng)轉(zhuǎn)速為0.2npm，經(jīng)雨流計數(shù)后舍棄小幅值循環(huán)，得到大幅值速循環(huán)和主要循環(huán)狀態(tài)在等壽命線中的分布情況如圖12所示。

圖12 大幅值循環(huán)和主要循環(huán)類型在等壽命圖中的分布

根據(jù)2.3節(jié)的距離判別公式，將造成疲勞損傷較大的51個轉(zhuǎn)速循環(huán)歸類為表4所示的11類循環(huán)。由表4所列的判別歸類結(jié)果可以看出，該型渦軸發(fā)動機主要循環(huán)類型為空慢-最大-空慢和巡航-最大-巡航。

3.3 載荷譜生成及方法驗證

將判別歸類后的轉(zhuǎn)速按循環(huán)開始時間進(jìn)行排序并連接。根據(jù)蠕變理論，各級載荷大小和保持時間相同的兩個載荷譜對零部件造成的蠕變損傷相同。主要狀態(tài)的保持時間可以由屬于該類別的低幅值循環(huán)的循環(huán)時間確定，統(tǒng)計該次飛行中5種類別的低幅值循環(huán)的持續(xù)時間[19]，得到5種主要狀態(tài)的持續(xù)時間比例，如表2所示。根據(jù)總飛行時間生成平均分配各狀態(tài)的持續(xù)時間，保證每個狀態(tài)的總持續(xù)時間與實測載荷譜相同，生成如圖13所示載荷譜。

表4 主要循環(huán)類型的計數(shù)結(jié)果

圖13 模擬任務(wù)載荷譜

渦軸發(fā)動機的使用可分為訓(xùn)練飛行、任務(wù)飛行和地面維護開車三部分，訓(xùn)練飛行是使用的主要部分。選取訓(xùn)練飛行某一科目3個架次的3臺發(fā)動機共9組飛行數(shù)據(jù)進(jìn)行對比分析，比較同一發(fā)動機不同飛行架次，以及同一飛行架次三臺不同渦軸發(fā)動機的主要狀態(tài)誤差情況、循環(huán)類別分散度情況，對比情況見表5及表6所示，序號1-1表示直升機第1飛行架次的1號發(fā)動機。由于同一飛行架次的三臺發(fā)動機使用載荷狀態(tài)較為接近，從表5可以看出，該方法識別的主要狀態(tài)誤差較小，保持在5%以內(nèi)。

表5 主要狀態(tài)誤差對比

表6 主要狀態(tài)誤差對比

從表6可以看出，同一飛行架次的各循環(huán)類型數(shù)量較為接近，第2和3飛行架次的2發(fā)循環(huán)數(shù)與1發(fā)循環(huán)數(shù)和3發(fā)發(fā)循環(huán)數(shù)有部分區(qū)別，主要由于2發(fā)進(jìn)氣環(huán)境不同導(dǎo)致發(fā)動機循環(huán)狀態(tài)與1發(fā)循環(huán)數(shù)和3發(fā)發(fā)循環(huán)數(shù)不完全匹配。

4 結(jié)論

1) 實際載荷歷程通過雨流計數(shù)后得到大量接近零幅值的小幅值循環(huán)，且其峰谷值呈現(xiàn)圍繞幾個核心狀態(tài)的簇狀分布，通過對其峰谷值的聚類確定發(fā)動機的實際主要工作狀態(tài)，同一飛行架次的三臺發(fā)動機主要工作狀態(tài)的識別誤差在5%以內(nèi)；

2) 實際使用轉(zhuǎn)速載荷循環(huán)以高均值小幅值循環(huán)為主，造成疲勞損傷較大的大幅值循環(huán)較少，但高均值低幅值循環(huán)造成熱沖擊及蠕變損傷不可忽略。推導(dǎo)出以轉(zhuǎn)速循環(huán)的峰谷值表示的等壽命線，利用該曲線進(jìn)行循環(huán)類型壓縮時可同時考慮疲勞、熱沖擊和蠕變損傷；

3) 該型渦軸發(fā)動機實際使用的地慢與空慢狀態(tài)與設(shè)計轉(zhuǎn)速較為接近，實際使用的巡航及最大載荷狀態(tài)遠(yuǎn)小于設(shè)計狀態(tài)；實際使用過程中的主要循環(huán)類型為空慢-最大-空慢和巡航-最大-巡航。