潘慕絢，陳強(qiáng)龍，周永權(quán)，周文祥,*，黃金泉

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

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

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

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

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

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

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

1 渦扇發(fā)動(dòng)機(jī)多動(dòng)力學(xué)模型

發(fā)動(dòng)機(jī)在工作中滿(mǎn)足如下假設(shè)：

1) 氣體沿發(fā)動(dòng)機(jī)軸向?yàn)橐痪S流動(dòng)，同一截面上氣體參數(shù)均勻。

2) 不考慮氣體動(dòng)力學(xué)方程中氣體的黏性和慣性力。

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

1.1 轉(zhuǎn)子轉(zhuǎn)速動(dòng)態(tài)模型

發(fā)動(dòng)機(jī)的動(dòng)態(tài)特性包括2個(gè)轉(zhuǎn)子的力矩平衡方程，即

(1)

(2)

式中：t為時(shí)間；nH、nL分別為發(fā)動(dòng)機(jī)高、低壓轉(zhuǎn)子轉(zhuǎn)速；JH、JL分別為發(fā)動(dòng)機(jī)高、低壓轉(zhuǎn)子的轉(zhuǎn)動(dòng)慣量；ΔMH、ΔML分別為高、低壓轉(zhuǎn)子的剩余力矩。

由發(fā)動(dòng)機(jī)原理可知，高、低壓轉(zhuǎn)子的剩余力矩分別為

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

(3)

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

(4)

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

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

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

燃燒室進(jìn)口總焓為燃油熱值、燃油完全燃燒的焓和進(jìn)口空氣焓構(gòu)成，即

(5)

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

(6)

令

(7)

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

(8)

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

(9)

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

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

(10)

進(jìn)一步，令

(11)

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

(12)

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

(13)

(14)

式中：τ為仿真步長(zhǎng)。

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

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

考慮參與渦輪換熱量大的部件：渦輪轉(zhuǎn)子葉片、靜子葉片、渦輪盤(pán)和機(jī)匣。將轉(zhuǎn)子與靜子葉片、渦輪盤(pán)與機(jī)匣分別視為一體，則它們與高溫燃?xì)饧八鼈冎g的換熱量為[18]

(15)

(16)

(17)

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

考慮熱傳遞引起的熱慣性，葉片溫度、機(jī)匣/渦輪盤(pán)溫度和渦輪出口燃?xì)鉁囟鹊膭?dòng)態(tài)特性可描述為

(18)

(19)

(20)

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

(21)

(22)

(23)

1.4 渦扇發(fā)動(dòng)機(jī)部件級(jí)多動(dòng)力學(xué)模型

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

1) 風(fēng)扇與壓氣機(jī)模型

(24)

(25)

式中：下標(biāo)std表示設(shè)計(jì)點(diǎn)參數(shù)。

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

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

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

(26)

(27)

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

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

(28)

式中：壓氣機(jī)出口焓h3的計(jì)算可參考文獻(xiàn)[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é)中的建模原理，逐次考慮熱慣性和容積效應(yīng)(圖5)，建模過(guò)程如下。

(34)

(35)

(36)

(37)

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

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

(38)

(39)

式中：

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

(40)

(41)

4) 外涵道模型

(42)

(43)

式中：

(44)

(45)

(46)

5) 其他部件模型

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

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

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

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

2 模型驗(yàn)證及結(jié)果分析

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

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

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

圖7 多動(dòng)力學(xué)模型、容積動(dòng)力學(xué)模型與迭代模型響應(yīng)曲線(xiàn)Fig.7 Response curves of multi-dynamic model, volume-dynamic model and iterative model

表2 單次流路計(jì)算耗時(shí)Table 2 Time consumption of single flow path calculation

ms

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

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

(47)

(48)

(49)

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

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

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

參數(shù)emax,Dy/%ess,avg/%Ie,Dy/%多動(dòng)力學(xué)模型容積動(dòng)力學(xué)模型迭代模型多動(dòng)力學(xué)模型容積動(dòng)力學(xué)模型迭代模型多動(dòng)力學(xué)模型容積動(dòng)力學(xué)模型迭代模型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 結(jié) 論

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

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

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

3) 多動(dòng)力學(xué)模型單次流路平均計(jì)算時(shí)間為0.009 ms。

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