基于主要狀態識別的渦軸發動機載荷譜研究

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

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

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

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

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

1 渦軸發動機主要工作狀態

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

1.1 雨流計數

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

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

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

圖1 雨流計數示意圖

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

表1 雨流計數結果

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

圖2 雨流矩陣示意圖

圖3 轉速循環在峰谷值等壽命線中的分布圖

1.2 主要工作狀態識別

對于轉動件，影響蠕變的主要參數為轉速。轉速保持相對穩定的階段即為發動機的穩定狀態，使用頻次大的轉速狀態則為發動機的主要工作狀態。

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

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

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

圖4 穩定狀態示意圖

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

圖6 數據分類情況及聚類中心圖

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

表2 主要狀態的時間分配情況

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

表3 主要循環類型

2 循環類型壓縮

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

2.1 等壽命線相關參數標準化

為在統一標準下進行比較，需將燃氣渦輪轉速ng換算至標準大氣條件下的換算轉速ng,c，并利用該型渦軸發動機燃氣渦輪最大轉速將其化為標準化轉速n，n∈[0,1]。

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

(1)

(2)

圖7 正弦轉速循環與應力循環曲線

根據發動機旋轉件的局部應力-應變特點，離心應力與轉速的平方成線性關系：

σ=σmv+a·n2

(3)

式中：σmv為0轉速時的殘余應力；a為轉速和應力之間的比例系數；n為轉速。

若該型發動機最大轉速ng,d對應的應力為σmp時，則標準化后的應力循環對應的均值σm、幅值σa及強度極限σb分別為

式中，σb為材料的強度極限，根據需要分析的旋轉件材料進行選取。

如圖8所示，等壽命線的表達式變為式(4)

(4)

圖8 等壽命線

2.2 等壽命線轉換

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

(5)

(6)

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

代入式(4)，可得

(7)

令

式(7)變為式(8)

Na=k(Nm-b)

(8)

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

圖9 轉速循環等壽命線

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

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

(9)

其中ki為不同等壽命線的斜率，i=1,2,…,10時分別對應轉速谷值為0、轉速峰值為0.1nmp、0.2nmp、…、0.9nmp的脈動循環所在等壽命線，如k9=kPB,k10=kPA?？梢缘玫揭苑骞戎当硎镜牡葔勖€如圖10所示，3條紅色邊界線nmax=1，nmin=0和nmax=nmin分別表示峰值轉速為最大轉速的循環、谷值轉速為最小轉速的循環和循環幅值為0的轉速循環?？梢钥闯?，推導出的等壽命循環的峰谷值關系滿足公式(9)，該公式的優勢在于在考慮疲勞損傷的同時，沒有忽略高均值低幅值循環造成的與狀態變化相關的熱沖擊及蠕變等損傷，可以在同一標準下根據主要循環狀態將造成熱沖擊損傷較大、疲勞損傷較小的高均值低幅值循環考慮進來，使得載荷譜更接近實際載荷狀態。

2.3 判別歸類

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

3 編制模擬任務試車譜

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

3.1 數據預處理

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

圖11 編制載荷譜流程框圖

3.2 生成載荷譜

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

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

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

3.3 載荷譜生成及方法驗證

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

表4 主要循環類型的計數結果

圖13 模擬任務載荷譜

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

表5 主要狀態誤差對比

表6 主要狀態誤差對比

從表6可以看出，同一飛行架次的各循環類型數量較為接近，第2和3飛行架次的2發循環數與1發循環數和3發發循環數有部分區別，主要由于2發進氣環境不同導致發動機循環狀態與1發循環數和3發發循環數不完全匹配。

4 結論

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

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

3) 該型渦軸發動機實際使用的地慢與空慢狀態與設計轉速較為接近，實際使用的巡航及最大載荷狀態遠小于設計狀態；實際使用過程中的主要循環類型為空慢-最大-空慢和巡航-最大-巡航。