董 楨,周文祥,潘慕絢 ,黃開明
(1.南京航空航天大學江蘇省航空動力系統重點實驗室,南京210016;2.中國航發湖南動力機械研究所,湖南株洲412002)
非線性部件級模型是航空發動機及其控制系統研制過程中不可或缺的仿真工具。模型的精度是建模過程中的1項關鍵技術指標,直接決定該模型的應用價值和適用范圍。借助高精度的航空發動機性能仿真模型,通過計算機仿真及優化,可以大大縮短發動機數控系統的設計及測試周期,降低研發成本,同時還可避免不必要的試車風險[1-2]。
影響發動機部件級模型精度的1個重要因素是發動機部件特性。對發動機而言,零部件制造和裝配過程必然存在著公差[3],并且在發動機壽命期內自然磨損、疲勞、腐蝕或積垢等原因也會引起相關部件的性能退化[4],從而導致設計出的發動機期望部件特性與實際裝機的部件特性之間產生差異。若以期望部件特性為依據建立的數學模型表現出的整機性能,與實際發動機整機性能之間也必然存在一定差異。因此,模型修正技術對于提高發動機部件級模型精度非常重要。
國外學者早在20世紀90年代初就開始進行修正發動機部件特性和模型參數的方法研究。Stamatis A等[5-6]提出1種基于發動機實際測量數據的方法,將發動機模型的仿真輸出和試車數據誤差的加權平方和作為目標函數,采用非線性廣義最小殘差法求解共同工作方程;同樣在試車數據的基礎上,Joachim Kurzk[7-8]著眼于大氣條件對發動機各截面總溫、總壓及流量等氣路參數的影響,給出詳細的修正方法。之后,優化算法逐漸引入部件特性修正中。Y.G.Li[9]采用牛頓-拉夫森算法求解不可測參數的修正因子,得到準確的不可測參數的真實值;Changduk Kong[10-12]從部件特性圖著手,發現壓氣機流量可以近似表示為壓比的3次多項式,效率可以近似表示為流量的3次多項式,其中多項式系數通過遺傳算法優化得到。國內學者在發動機部件級模型修正方面也進行了深入研究。段守付等[13]在模型中引入加權函數和一系列修正因子,改變修正因子的取值使加權函數最小,實現對發動機部件特性和模型的最優修正;陳策等[14]建立BP神經網絡對部件特性數據進行識別學習,實現特性數據的精確插值和對未知特性數據的推測;吳虎等[15]提出1種預測發動機部件特性的自適應模型方法,以通用特性為基礎,運用單純形優化方法預測出不同飛行條件下渦扇發動機的部件特性。上述修正方法在原理上具有相似性,但在求解非線性方程組或優化算法方面存在一定差異。
本文利用具有高效率和全局搜索能力的粒子群優化算法(Particle Swarm Optimization,PSO)來修正渦軸發動機的部件特性。在基本不改變部件級模型核心計算代碼的情況下,使用PSO算法計算出部件特性圖和其他參數的修正因子。
擬建模的渦軸發動機由進氣道(包括粒子分離器)、組合壓氣機(包括軸流壓氣機和離心壓氣機)、燃燒室、燃氣渦輪、動力渦輪和尾噴管組成,如圖1所示,各截面定義見表1。

圖1 渦軸發動機結構

表1 渦軸發動機各截面定義
渦軸發動機部件級模型的輸入參數是大氣壓力、環境溫度、高度、馬赫數、燃氣渦輪轉速和動力渦輪轉速。輸出參數是燃油流量、各截面氣動參數(總溫、總壓等)和性能參數(動力渦輪功率、耗油率等)。在建立所有部件的氣動熱力學模型后,選定4個仿真初猜值(壓氣機、燃氣渦輪和動力渦輪壓比及燃氣渦輪轉速),其對應的4個共同工作方程為
燃燒室出口和燃氣渦輪進口燃氣物理流量平衡方程

燃氣渦輪出口和動力渦輪進口燃氣物理流量平衡方程

動力渦輪出口和尾噴管進口燃氣物理流量平衡方程

軸功率平衡方程

式中:Wg為燃氣物理流量,下標數字為相應截面編號;Ngt為燃氣渦輪輸出功率;Ncp為組合壓氣機消耗功率;Nex為附件提取功率;ηmgt為燃氣渦輪軸機械效率;εi(i=1,2,3,4)為平衡方程的殘差。
利用牛頓-拉夫森法迭代求解共同工作方程組就可以得到穩態工作點。由上述4個共同工作方程結合各部件氣動熱力計算方程,構成發動機的穩態模型。在發動機穩態工作時,若共同工作方程的殘差|εi|<ε0,i=1,2,3,4(其中 ε0為控制精度,本文取 10-3),則認為方程收斂,停止迭代。
采用粒子群優化算法修正發動機部件級模型,其原理是利用PSO優化發動機模型的修正因子,使目標函數最小化或使適應度函數最大化。
對渦軸發動機部件級模型而言,待修正的參數包括旋轉部件的部件特性(包括壓比、換算質量流量和效率)、典型流道部件的總壓恢復系數、壓氣機引、放氣比例以及燃燒室效率。
修正設計點時,對上述參數都進行修正;修正非設計點時,本文只修正旋轉部件的部件特性,其他待修正參數采用設計點的修正結果。
選擇燃油物理流量Wf、組合壓氣機增壓比πcp、發動機進口空氣物理流量Wa2、動力渦輪進口燃氣總溫T45以及輸出功率Ne作為目標參數。
發動機測量參數的誤差定義為

式中:ycal為計算參數;yact為基準參數。
目標函數定義為

式中:ai為各目標參數的權重系數,本文在優化過程中暫取1。
粒子群算法是模擬鳥群覓食的過程[16-17]。每個優化問題的解都是搜索空間的1只鳥,稱為“粒子”。每個粒子具有由優化函數確定的適應度值。每個粒子也具有1個速度,決定搜索的方向和距離。所有粒子追隨當前最優粒子在解空間中搜索。
假設D維空間中的粒子群中的第i個粒子的位置表示為

速度表示為

當前粒子群中的最優個體粒子表示為

歷代全局最優個體粒子表示為

粒子群中的所有粒子根據式(11)和(12)更新速度和位置

式中:r1和r2取(0,1)中的隨機數;c1和c2為學習因子,通常取2;ω為慣性權重,取較大的值意味著粒子群具有較強的全局搜索能力,取較小的值意味著粒子群局部搜索能力達到最優。
所有粒子的速度和位置具有上下限,由實際問題決定。由式(11)和(12)組成的粒子群算法稱為基本粒子群算法。
在此基礎上,如果慣性權重的值根據迭代次數而減少,則算法被稱為慣性權重線性遞減粒子群算法[18]。一般來說,ω的初始值設為0.9,然后根據迭代次數將其線性減小到0.4。粒子群優化算法的慣性權重線性遞減是粒子群算法中比較常見的算法模型。本文采用慣性權重線性遞減粒子群優化算法優化渦軸發動機部件特性。
慣性權重定義為

式中:i為當前迭代次數;n為總迭代次數。
本文粒子群總優化代數設置為300,所有待修正因子優化范圍統一設置為[0.98,1.02]。適應度在初始100代以內急劇增加,然后減慢,在150代左右達到最大,如圖2所示。
部件特性修正流程如圖3所示。
需要特別說明的是,粒子群優化迭代過程中更新的是部件特性圖,而不是部件流量或效率特性的修正因子,這樣做的好處是一旦優化結束,可以直接輸出修正后的部件特性圖和其他有關模型參數,如總壓損失、引放氣比例等。

圖2 粒子群適應度變化
壓氣機流量-壓比特性如圖4所示。圖中描繪出待修正的一些穩態工作點,A為設計點,B、C、D為非設計點。對于設計點A,部件特性圖將根據模型計算的結果與發動機設計點參數進行整體縮放。非設計點則是根據相對換算轉速從高到低的順序依次進行修正。例如,當修正工作點B時,保持L1不變并縮放L2,L4,…,Ln,以免影響設計點A的仿真精度。然后固定L1和L2,繼續修正L4及更低的等換算轉速線來修正工作點C。
上述修正方法假定在壓氣機2條等換算轉速線之間只存在1個待修正工作點,這一假設條件在模型修正過程中有時是不滿足的。圖4顯示在壓氣機L2和L42條等換算轉速線之間存在C和D 2個工作點,針對這種情況,可以先選擇距離設計點A較遠的工作點D作為待修正點,對L4及其以下的等換算轉速線進行整體修正,然后用工作點C來檢驗修正好的部件特性是否達到仿真精度要求。如若不滿足要求,可以在點C和點D之間補充1條等換算轉速線L3。這時問題轉化成2條等換算轉速線之間只有1個待修正點的情況,可根據前述方法對補作的特性線L3及其以下特性線進行修正,修正完畢后,繼續根據點D的仿真結果修正L4及其以下轉速特性線,如圖5所示。

圖4 壓氣機流量-壓比特性

圖5 補作等換算轉速線后壓氣機流量-壓比特性
為了驗證提出的部件特性修正方案是否可行,本文從發動機性能分析軟件GasTurb提取渦軸發動機通用部件特性圖,以GasTurb軟件計算結果作為模型修正基準值。在地面標況下(ISA H=0 km,Ma=0)選取燃氣發生器轉子相對物理轉速分別為1、0.975、0.941、0.907的 4個工作點 A、B、C、D 來修正通用部件特性圖。設計點A修正前、后各目標參數的穩態誤差見表2。

表2 設計點A修正前、后穩態誤差對比(ISA H=0 km,Ma=0,n g=100%,n p=100%)
由表中可知,修正設計點A后,模型精度大大提高,設計點處模型平均穩態仿真誤差由4.7%下降至0.2%。設計點A修正前、后壓氣機流量特性變化如圖6所示。

圖6 設計點A修正前、后壓氣機流量特性對比
接著修正穩態工作點B。穩態工作點B修正前、后各目標參數的穩態誤差見表3。

表3 穩態工作點B修正前、后穩態誤差對比(ISA H=0 km,Ma=0,n g=97.5%,n p=100%)
穩態工作點B修正前、后壓氣機流量特性變化如圖7所示。
由于穩態工作點C、D同時在2條等換算轉速線之間,先修正穩態工作點D。穩態工作點D修正前、后穩態誤差見表4。
修正后 A、B、D 3 個穩態工作點目標參數相對誤差都在1%以內,滿足精度要求。然后檢驗修正后的部件特性是否能夠使穩態工作點C也滿足指標要求。修正后穩態工作點C穩態誤差見表5。

圖7 穩態工作點B修正前、后壓氣機流量特性對比

表4 穩態工作點D修正前、后穩態誤差對比(ISA H=0 km,Ma=0,n g=90.7%,n p=100%)

表5 修正后穩態工作點C穩態誤差(ISA H=0 km,Ma=0,n g=94.1%,n p=100%)
從表中可見,Wf和T45的相對誤差仍高于1%。因為這2個參數是反映渦軸發動機經濟性及整機性能指標非常重要的參數,很有必要進一步提升其仿真精度,所以此處考慮在點C、D之間補作1條等換算轉速線。補作的等換算轉速線需要根據實際情況來確定,這里取相鄰等換算轉速的中值。
補作等換算轉速線后穩態工作點C、D修正前、后穩態誤差見表6、7。修正前、后壓氣機流量特性對比如圖8、9所示,其中,L3為補作的等換算轉速線。

表6 補作等換算轉速線后工作點C修正前、后誤差對比(ISA H=0 km,Ma=0,n g=94.1%,n p=100%)

表7 補作等換算轉速線后工作點D修正前、后誤差對比(ISA H=0 km,Ma=0,n g=90.7%,n p=100%)
以上仿真結果表明,修正后點A、B、C、D目標參數穩態仿真誤差均在1%以內,滿足精度要求并且結果與預期相同。

圖8 穩態工作點C修正前、后壓氣機流量特性對比

圖9 穩態工作點D修正前、后壓氣機流量特性對比
壓氣機效率特性的修正與流量特性的修正方法一致,修正前、后壓氣機效率特性對比如圖10所示。
修正渦輪部件特性時,2條等換算轉速線之間存在B、C、D 3個穩態工作點。采用與壓氣機部件特性相同的修正方法,由于等換算轉速線分布比較密集,修正穩態工作點D之后的特性能夠使這3個穩態工作點同時滿足精度要求。修正前、后渦輪流量特性對比如圖11所示,渦輪效率特性對比如圖12所示。

圖10 修正前、后壓氣機效率特性對比

圖11 修正前、后燃氣渦輪流量特性對比

圖12 修正前、后燃氣渦輪效率特性對比
本文提出基于粒子群優化算法(PSO)修正渦軸發動機部件特性的方法,通過模型計算結果與Gas-Turb仿真結果的對比,得到以下結論:
(1)用粒子群算法優化并更新已有的部件特性圖可以有效減少部件級模型的穩態仿真誤差。
(2)在模型修正過程中,不更新特性圖只更新修正因子的方法對修正因子缺乏約束,求解出的數學解可能存在對應的發動機旋轉部件效率大于1的情況,這與實際情況不相符。而本文提出的方法在每次優化修正結束后直接輸出更新的部件特性圖,可以直觀地查看更新后的部件特性,有效避免了這一問題。
本文所提出的部件特性修正方法還有不足之處,例如在修正發動機非設計點的過程中,只是對發動機部件特性中的等換算轉速線進行平移或者縮放處理,并沒有改變其形狀或趨勢,有待引入更多的高空穩態工作點,進一步開展改變等換算轉速線局部形狀或趨勢的修正方法研究。