潘慕絢，陳強龍，周永權，周文祥,*，黃金泉

1. 南京航空航天大學 能源與動力學院，南京 210016 2. 中國航發(fā)控制系統(tǒng)研究所，無錫 214063

航空發(fā)動機是一類強非線性、時變、復雜的氣動熱力學系統(tǒng)，高精度和高置信度的非線性數(shù)學模型是其特性分析、先進控制策略、控制器設計、故障診斷、仿真與試驗驗證的重要依據(jù)。

早期發(fā)動機建模方法的研究側重于發(fā)動機穩(wěn)態(tài)性能的匹配，例如美國早期建立的單軸渦噴發(fā)動機、渦槳發(fā)動機模型，美國國家航空航天局(NASA)發(fā)布的動態(tài)通用發(fā)動機(DYNamic Generalized Engine,DYNGN)模型和數(shù)字化渦扇發(fā)動機模型(DIGital Turbofan Engine Models, DIGTEM)[1-2]，德國MTU(Motoren-und Turbinen-Union)公司推出的GasTurb模型和用于進行大涵道比渦扇發(fā)動機仿真的商用模塊化航空推進系統(tǒng)仿真(Commercial Modular Aero Propulsion System Simulation, CMAPSS)模型[3]，采用面向對象方法、組態(tài)技術等建立的渦扇、加力渦扇等不同類型發(fā)動機模型[4-6]。在給定初猜值、計算精度后，根據(jù)發(fā)動機氣動熱力學特性建立部件進出口參數(shù)之間的數(shù)學關系，采用牛頓-拉普森等迭代算法求解發(fā)動機流程中氣動參數(shù)[7]。模型中旋轉部件的參數(shù)，如發(fā)動機流量、效率通過部件特性圖獲得。這些模型中，大多假設發(fā)動機工作中氣體為定常流動，任意時刻部件吸、放熱均處于動態(tài)平衡狀態(tài)，從而僅考慮轉子動態(tài)[1,8]。因此，模型中所包含的對象動態(tài)信息十分有限，尤其是壓力等高頻動態(tài)信息。此外，雖然通過特性圖修正[9]、初猜值先驗優(yōu)化等手段提高了模型精度，但是對于初猜值的依賴一定程度上仍限制了模型精度。

實際上，隨著現(xiàn)代飛行器的機動性和靈活性需求日益提高，未來先進發(fā)動機其動態(tài)工作過程更加復雜多變。例如，傳統(tǒng)小涵道比渦扇發(fā)動機的工作速域、空域加寬，過渡過程要求更加嚴苛，因此溫度、壓力、轉速等關鍵參量的動態(tài)變化劇烈。又如，變循環(huán)發(fā)動機、組合沖壓發(fā)動機等先進發(fā)動機模態(tài)轉換時，發(fā)動機流道結構變化引起參數(shù)顯著的動態(tài)變化。發(fā)動機中除了轉子慣性之外，還存在氣流質量、能量的堆積與存儲[10]。隨著氣流速度、溫度梯度的增加，由此引起的溫度、壓力等參數(shù)動態(tài)過程愈加顯著[10-11]。這意味著發(fā)動機實際包含了“或快”“或慢”等豐富的動態(tài)信息，它們對于發(fā)動機動態(tài)性能分析、各類主動控制、直接性能量控制的有效實現(xiàn)有著直接影響。因此，針對航空發(fā)動機的動態(tài)特性研究受到越來越多的關注。

國內外研究者考慮發(fā)動機中氣流在容腔中質量和能量堆積效應，建立了容腔中氣流質量、動量和能量方程，獲得發(fā)動機動態(tài)參數(shù)模型[12-13]。這里，容腔可以是發(fā)動機某部件，如外涵道、燃燒室/加力燃燒室、壓氣機內級間容腔等，也可以是部件間容腔。根據(jù)選取的容腔大小的不同，可獲得不同截面上氣動參數(shù)。如文獻[13-14]將壓氣機每一級視為一個容器建立“級模型”，分析發(fā)動機喘振和失速。級模型中容腔多，具有所需特性數(shù)據(jù)多，計算量大，耗時長的特點。因此，面向控制系統(tǒng)設計的發(fā)動機建模中，需要權衡模型精度和計算實時性后，選擇合適的容腔及個數(shù)。例如，文獻[15] 考慮渦扇發(fā)動機主要部件間的容積效應，建立核心機和外涵道靜壓一致的假設，將模型狀態(tài)變量從9個減少為6個，減少了模型計算時間。文獻[16-17]中依據(jù)模型變量個數(shù)選取容腔，建立質量和能量方程，獲得動態(tài)模型。

針對發(fā)動機過渡態(tài)中部件吸/放熱過程，研究者們主要討論了能量儲存和釋放引起的溫度動態(tài)變化，尤其是高溫燃氣與機匣、葉片間的熱交換導致的溫度動態(tài)[18-19]，進一步由此引起的部件效率改變對部件出口參數(shù)產生影響[20-21]。這些研究中將機匣、輪盤和轉子葉片視為整體，未考慮不同部件材料和換熱系數(shù)對熱傳遞的影響，因而所得溫度精度不高。文獻[18]還針對葉片、渦輪盤以及機匣，分別構建了對應的熱傳遞方程，建立了發(fā)動機起動過程的渦輪部件的熱慣性模型。

盡管發(fā)動機動態(tài)過程受到多種動力學影響，但是目前發(fā)動機建模研究中未能綜合考慮轉子慣性、質量和能量堆積、能量存儲與釋放等引起的參數(shù)動態(tài)。如文獻[10,17]建模中僅考慮了轉子動力學和容積動力學，而忽略了熱慣性對參數(shù)動態(tài)的影響。文獻[20]中僅討論了個別部件中熱慣性對性能參數(shù)穩(wěn)態(tài)值的影響。文獻[22]中轉子動力學迭代模型(以下簡稱迭代模型)考慮精度和計算量，計算步長取為1 ms。該模型在主頻3.3 GHz CPU上單次流路計算平均耗時為0.1 ms，限制了其實時模擬發(fā)動機中高頻信號。因此，目前研究中的迭代模型和動態(tài)模型都存在發(fā)動機中關鍵參數(shù)的動態(tài)特性不夠完整，模型實時性難以提高，缺乏模型過渡態(tài)過程的試驗驗證等問題。

1 渦扇發(fā)動機多動力學模型

發(fā)動機在工作中滿足如下假設：

1) 氣體沿發(fā)動機軸向為一維流動，同一截面上氣體參數(shù)均勻。

2) 不考慮氣體動力學方程中氣體的黏性和慣性力。

3) 發(fā)動機各截面絕熱系數(shù)ka是該截面氣體總溫T*和余氣系數(shù)α的函數(shù)，即ka=ka(T*,α)。

1.1 轉子轉速動態(tài)模型

發(fā)動機的動態(tài)特性包括2個轉子的力矩平衡方程，即

(1)

(2)

式中：t為時間；nH、nL分別為發(fā)動機高、低壓轉子轉速；JH、JL分別為發(fā)動機高、低壓轉子的轉動慣量；ΔMH、ΔML分別為高、低壓轉子的剩余力矩。

由發(fā)動機原理可知，高、低壓轉子的剩余力矩分別為

ΔMH(t)=MTH(t)-MCH(t)-Mfr,H(t)

(3)

ΔML(t)=MTL(t)-MCL(t)-Mfr,L(t)

(4)

式中：MTH、MTL分別為高、低壓渦輪輸出力矩；MCH、MCL分別為壓氣機、風扇的輸入力矩；Mfr,H、Mfr,L分別為高、低壓轉子的摩擦力矩，由于其值相對較小，通常視為常數(shù)或忽略不計。

1.2 基于容積效應的部件性能動態(tài)模型

圖1 主燃燒室及其相關參數(shù)Fig.1 Main combustor and its relative parameters

燃燒室進口總焓為燃油熱值、燃油完全燃燒的焓和進口空氣焓構成，即

(5)

式中:HC,in為燃燒室進口總焓；ηC為燃燒室完全燃燒效率；Hf為燃油熱值；hC,in為燃燒室進口空氣焓；燃油完全燃燒的焓hp,C定義為

(6)

令

(7)

式中：hC,out為燃燒室出口焓；kaC,out為燃燒室出口氣體絕熱系數(shù)；αC,out為燃燒室出口余氣系數(shù)。燃燒室中燃氣容積慣性的非定常過程可描述為[12]

(8)

(WC,in(t)+Wf(t)-WC,out(t))+

(9)

式中：RC,out為燃燒室出口氣體常數(shù)；VC為燃燒室體積。

為避免DC計算過程中的微分求解，在此采用以下近似的方法求解:

(10)

進一步，令

(11)

(WC,in(k)+Wf(k)-WC,out(k-1))

(12)

由式(11)、式(12)可將式(8)、式(9)分別離散化為

(13)

(14)

式中：τ為仿真步長。

1.3 基于熱力學的部件性能動態(tài)模型

發(fā)動機過渡態(tài)工作中，高溫氣流與部件間的熱量交換過程存在熱慣性，通過熱交換動態(tài)方程能夠更真實反映發(fā)動機中相關截面的溫度動態(tài)。本文以高壓渦輪(圖2)為例，建立其能量方程，獲得高壓渦輪部件出口總溫。

考慮參與渦輪換熱量大的部件：渦輪轉子葉片、靜子葉片、渦輪盤和機匣。將轉子與靜子葉片、渦輪盤與機匣分別視為一體，則它們與高溫燃氣及它們之間的換熱量為[18]

(15)

(16)

(17)

圖2 高壓渦輪熱交換模型Fig.2 Heat exchange model for high pressure turbine

考慮熱傳遞引起的熱慣性，葉片溫度、機匣/渦輪盤溫度和渦輪出口燃氣溫度的動態(tài)特性可描述為

(18)

(19)

(20)

將式(15)～式(17)代入式(18)～式(20)中，并采用歐拉法離散化后可得

(21)

(22)

(23)

1.4 渦扇發(fā)動機部件級多動力學模型

某渦扇發(fā)動機結構及其截面定義如圖3所示。綜合考慮容腔進出口參數(shù)重要性、與其他參數(shù)的耦合程度、容腔大小和模型中待求參量數(shù)量等因素，在此選取渦扇發(fā)動機中5個容腔：燃燒室Ⅰ、高低壓渦輪間容腔Ⅱ、內涵噴口容腔Ⅲ(低壓渦輪出口和混合室進口前容腔)、加力燃燒室Ⅳ、外涵道Ⅴ。容腔選取的依據(jù)為：① 第Ⅰ、Ⅵ、Ⅴ容腔容積較大，能量、質量堆積效應較為明顯；② 第Ⅱ、Ⅲ容腔容積雖然較小，但是對高頻信號影響明顯；③ 基于第Ⅰ～Ⅴ容腔的動力學求解可避免與容腔相關的部件出口參數(shù)的迭代求解。

1) 風扇與壓氣機模型

(24)

(25)

式中：下標std表示設計點參數(shù)。

由壓氣機流量特性圖fW,Comp和效率特性圖fη,Comp可知壓氣機出口流量W3和效率ηc分別為

圖3 渦扇發(fā)動機截面示意圖Fig.3 Sectional diagram of turbofan engine

圖4 航空發(fā)動機多動力學模型流路計算框圖Fig.4 Calculation process of aeroengine multi-dynamic model flow path

(26)

(27)

式中：cW,Comp和cη,Comp分別為壓氣機流量和效率修正系數(shù)。

由πc、ηc和焓熵轉換函數(shù)fH2T，可得壓氣機出口總溫，即

(28)

式中：壓氣機出口焓h3的計算可參考文獻[23]。

2) 燃燒室(容腔Ⅰ)

由式(13)、式(14)和圖4可得

(29)

(30)

式中：

(31)

(W3(k)+Wf(k)-W4(k-1))

(32)

(33)

3) 高、低壓渦輪模型

高低壓渦輪之間和低壓渦輪之后存在容腔Ⅱ和容腔Ⅲ，且高低壓渦輪中存在的熱慣性，因此高低壓渦輪部件建模采用“串聯(lián)”思想，依據(jù)1.2節(jié)和1.3節(jié)中的建模原理，逐次考慮熱慣性和容積效應(圖5)，建模過程如下。

(34)

(35)

(36)

(37)

式中：W42(k)為高壓渦輪出口流量：C42為高壓渦輪出口截面的氣體比熱。

如圖5所示，高壓渦輪出口參數(shù)經熱慣性修正后，進一步考慮容積效應對其影響，由式(11)～式(14)，可得其出口參數(shù)為

(38)

(39)

式中：

圖5 高、低壓渦輪部件性能參數(shù)計算Fig.5 Calculation of performance parameters for high-pressure and low-pressure turbines

(40)

(41)

4) 外涵道模型

(42)

(43)

式中：

(44)

(45)

(46)

5) 其他部件模型

加力燃燒室模型與燃燒室的類似，在此不再給出。此外，渦扇發(fā)動機中進氣道、混合室以及尾噴管這些管道類部件建模中暫未考慮動力學特性，僅根據(jù)壓力損失、流量平衡等建立輸入輸出特性間的代數(shù)模型，文獻[24-25]中對此有詳細描述，本文也不再給出。

將本文建立多動力學模型中采用微分方程求解的模型動態(tài)參數(shù)與容積動力學模型、迭代模型中動態(tài)參數(shù)對比列于表1[8,12]。由表1可知，相較于后2種建模方法，本文方法可以獲得更多截面上的動態(tài)參數(shù)。

表1 模型參數(shù)求解方法對比Table 1 Comparison of model parameter solving methods

注：D為動態(tài)求解方法，-為非動態(tài)求解法/迭代求解法。

2 模型驗證及結果分析

針對某型渦扇發(fā)動機，采用本文建模方法建立其多動力學模型，其中取計算步長τ=1 ms，σⅠ= 0.97，σⅡ=0.98，σⅢ=0.98，σⅣ=0.95，σⅤ= 0.93，分別在設計點和非設計點進行數(shù)字仿真，多動力學模型輸出與發(fā)動機試驗數(shù)據(jù)間的相對誤差絕對值如圖6所示。由圖6可知，在設計點處，模型計算輸出與試驗數(shù)據(jù)間的相對誤差小于0.7%； 非設計點處，兩者間的相對誤差小于1.6%。 采用本文方法建立的渦扇發(fā)動機部件級模型具有良好的穩(wěn)態(tài)精度。

為了驗證本文建模方法對發(fā)動機動態(tài)過程描述的有效性和模型計算實時性，將多動力學模型、某容積動力學模型和某迭代模型的仿真結果進行對比。在標準大氣條件，高度H=0 km和馬赫數(shù)Ma=0時，令燃油量做20%的階躍，3個模型中典型參數(shù)的響應曲線如圖7所示。圖中括號內的數(shù)值表示迭代模型的計算收斂精度。3個模型在CPU主頻為3.2 GHz計算機中運行，單次流路計算耗時和平均耗時如表2所示。

圖6 設計點及非設計點參數(shù)相對誤差絕對值Fig.6 Absolute values of relative error of design point and non-design point parameters

圖7 多動力學模型、容積動力學模型與迭代模型響應曲線Fig.7 Response curves of multi-dynamic model, volume-dynamic model and iterative model

表2 單次流路計算耗時Table 2 Time consumption of single flow path calculation

ms

由表2可知，由于多動力學模型和容積動力學模型參數(shù)計算上一次通過，無迭代，其單次流路的平均計算耗時分別為9×10-3ms和9.4×10-3ms，是傳統(tǒng)迭代模型平均計算耗時的8%和8.5%。

定義動態(tài)過程最大誤差emax,Dy、誤差積分Ie,Dy和平均穩(wěn)態(tài)誤差ess,avg分別為

(47)

(48)

(49)

式中：t1為過渡態(tài)開始時刻；t2為過渡態(tài)終止時刻；t3為仿真終止時刻；y(t)為模型的輸出量；ys(t)為試車數(shù)據(jù)。

圖8 相同燃油量輸入下的模型數(shù)據(jù)與試車數(shù)據(jù)Fig.8 Model data and test data under the same fuel input

表3 動態(tài)過程誤差積分和穩(wěn)態(tài)誤差積分Table 3 Dynamic process error integral and steady state error integral

參數(shù)emax,Dy/%ess,avg/%Ie,Dy/%多動力學模型容積動力學模型迭代模型多動力學模型容積動力學模型迭代模型多動力學模型容積動力學模型迭代模型nL2.482.482.380.720.721.5114.2114.2174.3nH1.501.501.530.310.310.24107.8107.8102.7T?60.190.250.290.020.170.246.9428.2627.46p?64.925.805.472.612.702.78336.1343.1347.5

3 結 論

本文提出了一種航空發(fā)動機多動力學建模方法，建立了某混合排氣加力渦扇發(fā)動機部件級動態(tài)模型，獲得以下結論:

1) 與試車數(shù)據(jù)相比，在設計點處多動力學模型主要參數(shù)的穩(wěn)態(tài)誤差小于0.7%，非設計點處模型計算值的穩(wěn)態(tài)誤差小于1.6%。

2) 燃油階躍動態(tài)過程中，多動力學模型最大動態(tài)誤差小于5%。

3) 多動力學模型單次流路平均計算時間為0.009 ms。

4) 與傳統(tǒng)迭代模型相比，多動力學模型具有更好參數(shù)動態(tài)模擬能力和更高的實時性；與容積動力學模型相比，多動力學模型具有更精確的溫度響應和壓力響應過程和略快的計算速度。