曹惠玲，劉乃彬，薛 鵬，于天浩

（中國民航大學(xué)a.航空工程學(xué)院；b.工程技術(shù)訓(xùn)練中心，天津 300300）

航空發(fā)動機單元體性能與整機性能之間有很大聯(lián)系，其工作狀態(tài)的好壞直接影響到飛機的可靠運行以及機組成員和乘客的生命安全，因此，準(zhǔn)確分析單元體性能及其與整機性能的關(guān)系十分重要，可為發(fā)動機維修和故障診斷提供指導(dǎo)。由于發(fā)動機的部件特性和控制規(guī)律很難獲得，目前主要集中在通過實際QAR數(shù)據(jù)和試車臺數(shù)據(jù)進行擬合，建立各單元體效率與流量對發(fā)動機排氣溫度偏差值和燃油流量偏差值的貢獻率模型，確定單元體性能對整機性能的影響關(guān)系。唐海龍等[1]應(yīng)用蒙特卡洛概率設(shè)計方法建立渦軸發(fā)動機部件機非確定性非線性性能模型，利用發(fā)動機各部件效率和流量，研究某工況下各部件性能的非確定性對發(fā)動機性能的影響。趙琳等[2]通過辨識建模的方法先建立穩(wěn)態(tài)模型，在其基礎(chǔ)上又建立動態(tài)模型，對發(fā)動機多組試車數(shù)據(jù)進行辨識檢驗，研究發(fā)動機全包線、全工況的實時性能參數(shù)變化。航空發(fā)動機是一個非常復(fù)雜的熱力機械系統(tǒng)，各單元體工作過程不同，但彼此有大量聯(lián)系，其工作過程可用相應(yīng)的熱力方程表示，但方程存在強烈的非線性，要獲得發(fā)動機性能數(shù)據(jù)，需要大量的迭代計算求解[3]。因此需要從發(fā)動機原理出發(fā)，結(jié)合實際QAR數(shù)據(jù)，將非線性方程轉(zhuǎn)化為線性方程，去除迭代計算，建立反映單元體性能與整機性能關(guān)系的數(shù)學(xué)模型。Urban[4-5]于20世紀(jì)60年代提出了一種線性模型氣路分析法，應(yīng)用小偏差理論研究燃氣輪機故障，為發(fā)動機非線性的氣動熱力方程轉(zhuǎn)化為線性方程提供了依據(jù)。

根據(jù)發(fā)動機工作原理以及各單元體的熱力方程，應(yīng)用小偏差理論，建立了各單元體的小偏差模型。又根據(jù)各單元體間的流量平衡關(guān)系和壓比平衡，結(jié)合上述各單元體模型，建立了發(fā)動機部件性能和整機性能參數(shù)關(guān)系的數(shù)學(xué)模型。通過此模型研究部件特性與整機性能之間的關(guān)系。

1 小偏差數(shù)學(xué)理論基礎(chǔ)

小偏差法是一種近似計算方法，通過自變量發(fā)生的微小變化確定因變量的相應(yīng)變化量[6]。設(shè)發(fā)動機系統(tǒng)的熱力參數(shù)之間存在以下函數(shù)關(guān)系

其中：x為性能參數(shù)；y為可測參數(shù)（用于部件性能參數(shù)對整機參數(shù)的求解結(jié)果）。在基值（x10，x20，…，xm0）領(lǐng)域內(nèi)存在n+1階的連續(xù)偏導(dǎo)數(shù)，則函數(shù)y可在基值領(lǐng)域內(nèi)按泰勒級數(shù)[7-8]展開，即

忽略上式右邊的二次項及以后各項，簡化得

上式表示性能參數(shù)偏差和可測參數(shù)偏差之間的線性關(guān)系方程，且各參數(shù)之間保持原有的量綱單位，為簡化計算，減少不同量綱之間的影響，用式（3）等號兩端分別除以式（1）等號兩端進行無量綱化處理，得到表達式為

其中：a1，a2，…，am為方程系數(shù)為可測參數(shù)偏差量;為性能參數(shù)偏差量。

2 發(fā)動機數(shù)學(xué)模型建立

2.1 主要部件組成及工作過程

雙軸渦扇發(fā)動機目前應(yīng)用較廣泛，其主要部件有進氣道（Inlet）、風(fēng)扇（Fan）、低壓壓氣機（LPC）、高壓壓氣機（HPC）、燃燒室（C/C）、低壓渦輪（LPT）、高壓渦輪（HPT）、尾噴管（Nozzles）。將發(fā)動機工作過程理想化，忽略漏氣、壓氣機引氣的損失；不考慮外涵道影響，將發(fā)動機簡化成只有1個涵道，氣體的流動為一元流動；尾噴管為不可調(diào)收斂性噴管且氣體完全膨脹。在此基礎(chǔ)上根據(jù)部件工作過程建立發(fā)動機故障診斷數(shù)學(xué)模型。發(fā)動機具體站位劃分如圖1所示。

圖1 發(fā)動機站位圖Fig.1 Engine station diagram

2.2 單元體模型建立

應(yīng)用小偏差理論對發(fā)動機各單元體的熱力方程進行偏差處理，得到每個單元體的小偏差方程，如式（8）～（19）所示。其中Ki為各偏差方程系數(shù)；Ti*和Pi*分別為各單元體的截面參數(shù)溫度和壓力，如圖1所示；π為壓比，不同單元體的壓比分別用不同的下角標(biāo)表示；N1為低壓轉(zhuǎn)速；N2為高壓轉(zhuǎn)速。首先，風(fēng)扇工作過程的熱力方程如下。

耗功方程為

效率方程為

其中：T2*為風(fēng)扇進口溫度；T23*為風(fēng)扇出口溫度（低壓壓氣機進口溫度）為風(fēng)扇壓比；P2*為進口壓力；P23*為風(fēng)扇出口壓力（低壓壓氣機進口壓力）；γ為空氣系數(shù)；R為熱力學(xué)常數(shù)。

對上式進行無量綱化，得到風(fēng)扇耗功小偏差方程為

同理，可得效率的小偏差方程為

同理，可得低壓壓氣機小偏差方程為

高壓壓氣機小偏差方程為

低壓渦輪小偏差方程為

燃燒室燃燒過程能量平衡方程為

其中：Wf為燃油流量；Hu為燃油低熱值；ηB為燃燒效率；Cp為空氣定壓比熱容；Gc為發(fā)動機進氣量；f為油氣比。

對應(yīng)的小偏差方程為

2.3 性能參數(shù)變化量分離

航空發(fā)動機發(fā)生故障時，其單元體性能參數(shù)的變化包含兩方面：①各單元體的自身故障使其特性發(fā)生變化，造成各單元體效率和流量的改變，如圖2～圖3中的AA1過程，A為發(fā)動機發(fā)生故障前的工作點，部件故障引起特性線移動到A1，此時流量和效率的變化量為和δηF0；②發(fā)生故障后，各單元體要重新匹配達到新的平衡工作點，也會引起效率和流量的改變[9]，如圖2～圖3中，A1A2′過程為發(fā)生故障后，部件重新匹配，滿足功率與轉(zhuǎn)速平衡的條件，移到了新平衡點A2′。通過單元體通用特性曲線工況點的變化可以確定發(fā)生故障時單元體性能參數(shù)的獨立變化量，排除其他單元體性能參數(shù)的影響。

圖2 發(fā)生故障時折合流量與效率和轉(zhuǎn)速的工況點變化Fig.2 Change of flow rate,efficiency and rotation speed in case of failure

圖3 發(fā)生故障時折合流量與壓比和轉(zhuǎn)速的工況點變化Fig.3 Change of flow rate,pressure ratio and rotation speed in case of failure

將A1A2′過程分為沿新的等轉(zhuǎn)速線上變化的A1A2和沿平衡工作線變化的A2A2′兩個過程。其中A1A2過程的流量和效率分別為GA1和ηA1，A2A2′過程的流量和效率為GA2和ηA2。所以由發(fā)生故障前的點A到新的平衡工作點A2′，部件的流量和效率變化分別由3部分組成，用小偏差方程表示為

A1A2過程在等轉(zhuǎn)速線上移動，轉(zhuǎn)速不變，流量和效率為關(guān)于壓比的方程，表示為

對應(yīng)的小偏差方程為

A2A2′過程關(guān)于壓比的小偏差方程為

又由于此過程轉(zhuǎn)速不斷變化，設(shè)πA2A2′=φ（nA2A2′），其偏差方程為

將式（28）代入式（26）～（27），有

由于AA1過程為等壓比過程，其變化量為0，所以風(fēng)扇總壓比變化量的偏差方程由兩部分組成，即

將式（24）、式（29）、式（31）～（32）代入式（20），得到風(fēng)扇流量的小偏差方程為

同理得：

同理可得其余單元體流量和效率的分離偏差方程。將分離的性能參數(shù)偏差方程[10]代入2.2節(jié)建立的各單元體小偏差方程，得到各單元體偏差方程，如表1所示。

2.4 矩陣模型建立

發(fā)動機在工作過程中，各單元體之間存在內(nèi)在聯(lián)系。運用雙軸燃氣輪機中轉(zhuǎn)子功率相等、各單元體間流量平衡、壓比平衡等條件，對式（19）、式（35）～（50）進行聯(lián)立求解。求解流程如圖4所示。

表1 單元體偏差方程Tab.1 Engine unit deviation equation

圖4 矩陣模型建立過程Fig.4 Building process of matrix model

通過以上流程圖的求解，得出關(guān)于低壓轉(zhuǎn)速（N1）、高壓轉(zhuǎn)速（N2）、排氣溫度EGT（T5*）和燃油流量（Wf）的表達式，即

將上述結(jié)果以矩陣表示，求得方程形式為y=Ax

為求解上述數(shù)學(xué)模型的系數(shù)矩陣A，需要各單元體進出口的總溫和總壓，雖然可以通過理論計算獲得，但在部件特性未知的情況下，很多因素需要假設(shè)，導(dǎo)致計算結(jié)果誤差很大。為了使計算結(jié)果與發(fā)動機實際工作情況相近，這些參數(shù)可通過發(fā)動機實際運行QAR數(shù)據(jù)獲得。每次飛行時所記錄的QAR數(shù)據(jù)都可以認(rèn)為是發(fā)動機的一次正常飛行試驗數(shù)據(jù)。選用PW4000系列發(fā)動機QAR記錄中某一巡航狀態(tài)的無故障數(shù)據(jù)，并根據(jù)參數(shù)修正公式對所選參數(shù)進行修正，結(jié)果如表2所示。

表2 修正后的QAR數(shù)據(jù)Tab.2 Corrected QAR data

各單元體的實際效率根據(jù)發(fā)動機非設(shè)計點模型計算得出[11-12]。將所選數(shù)據(jù)和計算獲得的效率代入偏差方程系數(shù)Ki中，得到該巡航狀態(tài)的系數(shù)矩陣A，如表3所示。

表3 數(shù)學(xué)模型矩陣Tab.3 Mathematical model matrix

3 模型驗證與應(yīng)用

為驗證模型的有效性，采集PW4077D指印圖[13]上有關(guān)單元體的性能衰退征兆量，如圖5所示。采集數(shù)據(jù)如表4所示。

圖5 PW4077D指印圖Fig.5 PW4077D fingerprints

表4 單元體性能參數(shù)偏差Tab.4 MDL performance parameter change

表5 整機性能參數(shù)變化量及誤差比較結(jié)果Tab.5 Parameter change and comparison results of whole engine

由模型計算的參數(shù)變化與采集的性能參數(shù)變化可知，模型計算的結(jié)果與采集的數(shù)據(jù)較接近，大部分相對誤差小于10%，個別誤差不超過15%。誤差產(chǎn)生的原因可能為沒有考慮壓氣機的引氣損失、各單元體的機械效率以及假設(shè)了氣體在噴管內(nèi)完全膨脹，噴管為單一涵道，忽略了內(nèi)外涵道之分。雖然這些會對模型的計算結(jié)果造成一定的誤差，但通過此方法建立的模型具有可行性。根據(jù)此模型結(jié)合單元體發(fā)生不同故障或同一故障嚴(yán)重程度不同時性能的衰退量可以計算整機性能參數(shù)的變化量。建立了單元體性能與整機性能之間的聯(lián)系。

4 結(jié)語

1）在熱力計算基礎(chǔ)上結(jié)合實際QAR數(shù)據(jù)建立了數(shù)學(xué)模型，并通過模型得到了轉(zhuǎn)速、EGT（T5*）、燃油流量等性能參數(shù)的偏差值，偏差值誤差大部分小于10%，表明計算結(jié)果比較準(zhǔn)確，從而通過模型確定了單元體性能與整機性能參數(shù)之間的變化規(guī)律。通過模型得到EGT（T5*）、FF的變化值，可以確定單元體對整機性能影響的大小，決定拆發(fā)時間，為發(fā)動機以單元體為單位的重點維修提供依據(jù)。

2）應(yīng)用小偏差理論，對非線性方程進行線性化處理，相較于傳統(tǒng)部件級建模方法利用非線性方程組對性能參數(shù)變化的求解，減少了大量的迭代計算，提高了計算速度，節(jié)省了計算時間。