渦扇發動機多動力學建模方法

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

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

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

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

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

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

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

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

1 渦扇發動機多動力學模型

發動機在工作中滿足如下假設：

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

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

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

1.1 轉子轉速動態模型

發動機的動態特性包括2個轉子的力矩平衡方程，即

(1)

(2)

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

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

Δ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分別為高、低壓轉子的摩擦力矩，由于其值相對較小，通常視為常數或忽略不計。

1.2 基于容積效應的部件性能動態模型

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

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

(5)

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

(6)

令

(7)

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

(8)

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

(9)

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

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

(10)

進一步，令

(11)

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

(12)

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

(13)

(14)

式中：τ為仿真步長。

1.3 基于熱力學的部件性能動態模型

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

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

(15)

(16)

(17)

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

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

(18)

(19)

(20)

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

(21)

(22)

(23)

1.4 渦扇發動機部件級多動力學模型

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

1) 風扇與壓氣機模型

(24)

(25)

式中：下標std表示設計點參數。

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

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

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

(26)

(27)

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

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

(28)

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

2) 燃燒室(容腔Ⅰ)

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

(29)

(30)

式中：

(31)

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

(32)

(33)

3) 高、低壓渦輪模型

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

(34)

(35)

(36)

(37)

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

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

(38)

(39)

式中：

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

(40)

(41)

4) 外涵道模型

(42)

(43)

式中：

(44)

(45)

(46)

5) 其他部件模型

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

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

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

注：D為動態求解方法，-為非動態求解法/迭代求解法。

2 模型驗證及結果分析

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

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

圖6 設計點及非設計點參數相對誤差絕對值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可知，由于多動力學模型和容積動力學模型參數計算上一次通過，無迭代，其單次流路的平均計算耗時分別為9×10-3ms和9.4×10-3ms，是傳統迭代模型平均計算耗時的8%和8.5%。

定義動態過程最大誤差emax,Dy、誤差積分Ie,Dy和平均穩態誤差ess,avg分別為

(47)

(48)

(49)

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

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

表3 動態過程誤差積分和穩態誤差積分Table 3 Dynamic process error integral and steady state error integral

參數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 結 論

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

1) 與試車數據相比，在設計點處多動力學模型主要參數的穩態誤差小于0.7%，非設計點處模型計算值的穩態誤差小于1.6%。

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

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

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