基于改進牛頓法的航空發動機起動建模技術

姜 波，王繼強，何小龍，胡忠志

（南京航空航天大學江蘇省航空動力系統重點實驗室，南京210016）

0 引言

航空發動機由靜止逐漸加速至慢車狀態這一過程稱為起動過程。這個過程發動機氣動熱力參數變換劇烈，如何保證這個過程成功、快速的完成對延長發動機壽命、保證飛機安全有著重要的意義。對該過程進行數學抽象，建立起動過程的數值計算模型不僅能夠為設計人員提供許多有價值的信息，而且能夠讓設計人員在發動機設計的初期就展開合理的起動過程控制規律優化研究。因而建立能夠準確反映發動機起動過程的數值仿真模型對發動機的理論研究及工程應用都有著極大地意義。

在20世紀80年代，國外已經展開了發動機慢車以下狀態的氣動熱力學建模研究，M.Yunis和R.K.A-grawal在1980年基于流量平衡的假設，通過外推方法獲得低轉速部件特性從而建立起了單轉子發動機的起動模型[1]。在90年代NASA的Chappell采用慢車以上建模的方法，建立了包含高空起動、風車起動和地面起動的AVEST-V3模型[2]。此外在2001年和2007年，英國弗吉尼亞理工大學[3]英國克蘭菲爾德大學[4]的學者在上述研究的基礎上針對部件特性外推公式及零轉速附近特性修正展開了深入的研究。國內學者對此研究展開的時間較晚，但也取得了大量的研究成果，如：基于大量發動機試車數據，使用系統辨識的方法建立的實時起動模型[5-8]；基于容積動力學和部件特性外推的方法建立雙軸混合排氣的渦扇發動機的地面及空中起動模型[9-10]，利用壓氣機及渦輪特性數據的脊柱特性建立能夠準確模擬單軸發動機風車起動的部件級模型[11]等。但這些研究大都基于大量的試驗數據并且在起動初始點需要根據經驗反復的試湊出十分準確的初始值，工作量大且過程較為復雜。

本文利用NASA公開的PW公司JT9D大涵道比雙轉子渦扇發動機慢車以上部件級模型，通過部件特性外推、起動附件建模、修正部件總壓恢復系數和改進模型求解算法建立了JT9D發動機起動至最大轉速整個過程的仿真模型，并通過仿真結果驗證了模型的準確性及改進求解算法的有效性。

1 JT9D發動機部件級模型

建立的起動過程模型是基于美國PW公司于20世紀60年代年設計制造的大涵道比雙轉子渦扇發動機JT9D的慢車以上狀態的部件級模型，具體結構如圖1所示，具體的建模過程這里不再贅述，參照文獻[12]。

圖1 JT9D發動機部件結構

2 雙軸渦扇發動機起動過程

雙軸渦扇發動機在起動時沒有空氣流經整臺發動機，因此燃燒室無法點火。只有當轉子達到一定的轉速后，燃燒室內的氣流才能建立穩定燃燒所需的氣流壓力和溫度。因而發動機地面起動過程必須依靠起動機帶轉，一般情況下，雙軸渦扇發動機的起動過程（如圖2所示）分為以下3個階段。

圖2 雙軸渦扇發動機起動過程扭矩特性

第1階段：燃燒室未供油，起動機帶動高壓轉子，起動機扭矩MST大于起動過程阻力矩MC/ηm，高壓轉子轉速逐漸增加到n1。

第2階段：燃燒室開始供油并點火，渦輪開始產生力矩，高壓轉子轉速逐漸增加到n2。

第3階段：起動機脫開，渦輪產生足夠的力矩并將高壓轉子轉速獨立帶轉至n3。

3 起動建模關鍵問題及技術

3.1 低轉速部件特性獲取

發動機的部件特性在高轉速時差別較大，但在低轉速范圍內卻有著很大的相似性[1]。因而，采用部件特性指數外推的方法，能夠在缺少試驗特性的情況下得到較為合理的低轉速部件特性。采用的指數外推公式如下

式中：p、q、r為外推指數；cor表示換算狀態；ref表示參考狀態；Wc為換算流量；N為轉速；PR為增壓比（渦輪特性為膨脹比ER）；η為效率。以風扇為例，選取換算轉速為1的特性線為參考線。利用原有的特性數據通過不斷變換外推指數的數值，使得不同換算轉速線上最高效率狀態點連成的線近似為線性（效率僅在低轉速Ncor＜0.8時成線性）進而確定外推指數。不同外推指數下風扇流量、壓比、效率、隨換算轉速變化分別如圖3～5所示，不同外推指數下的外推結果對比如圖6所示。

圖3 不同外推指數下風扇流量隨換算轉速變化

圖4 不同外推指數下風扇壓比隨換算轉速變化

圖5 不同外推指數下風扇效率隨換算轉速變化

在確定了外推指數后，以換算轉速為0.6的試驗特性線為參考，外推換算轉速為0.5的特性線并與原有試驗特性對比，結果顯示：按照不同換算轉速線上的最高效率點連成的線近似為線性，所選取的外推指數能更好的反映低轉速的部件特性。據此確定外推指數見表1。

圖6 不同外推指數下的外推結果對比

表1 部件特性外推指數

在確定部件特性的外推指數后，以各部件試驗的最低狀態換算轉速線為外推參考線，外推各部件的風扇換算流量-壓比特性、風扇換算流量-效率特性、高壓渦輪換算流量-膨脹比、高壓渦輪換算流量效率特性特性分別如圖7～10所示。

圖7 風扇換算流量-壓比特性

3.2 起動機建模

JT9D發動機使用美國漢密爾頓公司提供的PS-700空氣渦輪起動機[13]，主要參數見表2。

空氣渦輪起動機在發動機起動過程的輸出扭矩與高壓轉子轉速近似為線性關系[14]。起動機扭矩特性如圖11所示。

圖8 風扇換算流量-效率特性

圖9 高壓渦輪換算流量-膨脹比特性

表2 PS-700空氣渦輪起動機主要參數

根據扭矩特性按照如下公式建立起動機模型

式中：P為起動機輸出功率；n2為起動機脫開時高壓轉子轉速；N為高壓轉子轉速；T為起動機輸出扭矩。

圖11 起動機扭矩特性

3.3 改進牛頓法及其應用

在求解JT9D發動機共同工作非線性方程組時NASA提供了1種傳統的Newton-Raphson迭代求解方法。此方法在求解發動機慢車以上狀態共同工作非線性方程組時能夠在選定較為準確的迭代初始值情況下達到快速收斂的效果[15-17]。但發動機起動過程部件特性的相對值變化劇烈，即使在給出相對準確的迭代初始值情況下也需要多次求解才能達到收斂精度的要求。并且每次迭代求解共同工作非線性方程組時計算機都需要對模型進行氣動熱力計算，然后再對平衡方程偏導數組成的Jacobian矩陣進行中心差分，龐大的計算量降低了模型的執行效率，使得計算時間變長，模型的實時性降低。

本文所述的改進牛頓法是參考廣泛用于國外先進數值仿真程序，如美國航空航天局NASA的NPSS（Numerical Propulsion Simulation System）[18]、荷蘭國家航空航天實驗室（NLR）的GSP（Gas turbine Simulation Program）[19]等在求解非線性方程組時使用的Broyden法[20]：即在傳統Newton-Raphson方法基礎上改變Jacobian矩陣的迭代更新方式，進而減少模型的氣動熱力計算次數，提高模型的仿真求解速度。計算方法如下：

對于發動機的共同工作非線性方程組

式中：x=（x1,x2,…,xn）T，n 為平衡方程數，xi為發動機的某參數；εi為平衡方程的誤差量；fi（x）代表不同的平衡方程，且F（x）=（f1（x1）,f2（x2）,…,fn（xn））T，Jacobian矩陣為

其中：偏導數項在模型求解計算時需根據中心差分求得。根據Newton-Raphson求解方法

式中：k為迭代次數，此時構造如下矩陣

根據Broyden法求解原理[20]，Jacobian矩陣迭代更新的修正項為

使用Sherman-Morrison公式解析的求Jacobian矩陣的逆，從而可得Jacobian矩陣在更新時

當 |F（xk）|＜δ（設定誤差允許量）時，認為迭代收斂。

這種改進的牛頓法在求解發動機共同工作非線性方程組時第1步仍采用中心差分的方法得到Jacobian矩陣，而后續的迭代更新只需通過計算ΔJk就可完成。這樣模型的動態仿真過程就避免了每次通過中心差分方法重新計算Jacobian矩陣而帶來的龐大計算量，從而大大提高了需要多次迭代求解發動機共同工作方程的起動過程模型的執行效率。

3.4 冷轉動過程模型

發動機在點火前的工作狀態即起動過程的第一階段稱為冷轉動過程。這一過程渦輪和壓氣機都是負載，建模方法如下：

首先計算起動機和壓氣機的扭矩，然后根據轉子動力學計算高壓轉子轉速

式中：ΔT為起動機扭矩與負載扭矩差；J為高壓轉子轉動慣量；N為高壓轉子轉速；t為時間。

之后，確定1個合適的點火狀態，根據此時的轉速 nR，WR和壓比 PR利用式（13）、（14）計算對應時刻下的參數

最后，低壓轉子轉速及其他截面狀態參數按高壓轉子轉速進行線性插值即可得到。

3.5 部件總壓恢復系數修正

NASA公開的JT9D發動機模型各部件的總壓恢復系數都取了定值，這在慢車以上狀態建模時是適用的。但在發動機的起動過程中由于各部件的流量、壓比、轉速都很小，真實的總壓恢復系數更接近于1，并且隨著轉速的增大而逐漸減小[9]。因而在建立起動模型時，各部件的總壓恢復系數在點火前取近似略小于1的常數，而在點火狀態和慢車狀態之間時按照轉速進行線性插值

式中：a為點火狀態；b為慢車狀態；σ為總壓恢復系數。

4 仿真結果分析

NASA公開的JT9D發動機模型的仿真范圍為：70%NHmax～100%NHmax。按照上述建模方法，補充建立包含起動過程的全狀態性能模型，并在地面標準狀況下給定供油規律（如圖12所示），仿真結果分別如圖13～15所示。

圖12 供油曲線

圖13 轉速變化曲線

圖14 推力變化曲線

圖15 燃燒室出口溫度變化曲線

由仿真結果可以看出，起動過程的第1階段燃燒室沒有供油。高壓轉子由起動機帶動，轉速緩慢上升并帶動低壓轉子轉速上升，此時推力也緩慢上升，燃燒室出口溫度基本不變。在第16 s時開始供油，轉子轉速和燃燒室出口溫度急劇上升至慢車狀態。在第40 s以后增加供油，發動機由慢車狀態逐漸運行到最大轉速狀態，轉子轉速、推力及燃燒室出口溫度都增加至100%。整個模型能很好的反映起動過程的各階段并且與原有慢車以上狀態模型有效地連接，驗證了建模方法的有效性。

為了進一步驗證建模方法的準確性，將模型的仿真結果與僅有的ICAO（國際民航組織）公布的JT9D-7J發動機數據進行對比，結果見表3。

表3 7%推力狀態下（52%NHmax）模型仿真結果驗證

各部件總壓恢復系數修正前（修正前采用固定值σb）與修正后高壓轉子轉速的仿真結果對比如圖16所示，燃燒室點火至慢車狀態過程的高壓轉子轉速最大相對誤差為1%。

圖16 總壓恢復系數修正前后仿真結果對比

同時，模型在點火后由各部件特性迭代求解的初始點計算時間最長，最難收斂。因而在此處用改進牛頓法進行驗證效果最為明顯。

單個步長Jacobian矩陣計算時間見表4。表4計算的硬件環境為：AMD A6-3400M（CPU），主頻1.4 G Hz，內存4 GB。由結果可以看出改進牛頓法的求解速度要明顯高于傳統的Newton-Raphson方法。并且，基于改進牛頓法的原理，當發動機共同工作方程和獨立變量個數增加以及迭代初始值偏離真值程度越大時，這種效果將會變得更加明顯。

表4 單個步長的Jacobian矩陣計算時間

5 結論

通過對JT9D發動機部件級起動過程建模、使用改進牛頓法替換傳統的Newton-Raphson方法仿真求解，得到以下結論：

（1）采用部件特性指數外推法并使用不同換算轉速線上最高效率狀態點連成的線，近似為線性來確定外推指數，可以得到合理的各部件低轉速部件特性；

（2）起動附件建模及轉子冷轉動過程建模能夠有效地幫助建立完整的起動過程模型。

（3）發動機起動過程各部件總壓恢復系數按照轉速變化修正能夠有效的提高模型的精度。

（4）采用改進牛頓法替代傳統的Newton-Raphson方法可以改善模型的執行效率，大幅提高模型的仿真求解速度。

[1]Agrawal R K,Yunis M.A generalized mathematical model to estimate gas turbine starting characteristics[J].Physical Review B,1981（1）:5668-5687.

[2]Chappell M A,Mclaughlin P W.Approach of modeling continuous turbine engine operation from startup to shut down[J].Journal of Propulsion&Power,2015,9（3）:466-471.

[3]Sexton W R.A new method to control turbofan engine starting by varying compressor surge valve bleed[D].Virginia:Virginia Polytechnic Institute and State University,2001.

[4]Howard,Jason.Sub-idle modeling of gas turbines:altitude relight and wind milling.[D].Cranfield:Cranfield University,2007.

[5]于達仁,郭鈺鋒,牛軍,等.渦噴發動機風車啟動工況的神經網絡建模[J].推進技術,2001,22（3）:183-186.YU Daren,GUO Yufeng,NIU Jun,et al.Turbojet modeling in wind milling based on radial basis function networks[J].Journal of Propulsion Technology,2001,22（3）:183-186.（in Chinese）

[6]卓剛.航空發動機智能建模與故障診斷研究[D].南京：南京航空航天大學,2004.ZHUO Gang.The research on aeroengine intelligent modeling and fault diagnosis[D].Nanjing:Nanjing University of Aeronautics and Astronautics.（ in Chinese）

[7]劉建勛,王劍影,李應紅,等.某型發動機起動模型的支持向量機辨識及應用[J].推進技術,2004,25（5）:401-404.LIU Jianxun,WANG Jianying,LI Yinghong,et al.An identification model of aero engine starting based on support vector machine and its application[J].Journal of Propulsion Technology，2004，25（5）:401-404.（in Chinese）

[8]李永進,張海波,張天宏.一種考慮非線性余項的機載發動機自適應模型建立及其在尋優控制中的應用 [J].推進技術,2016,37（1）:172-180.LI Yongjin，ZHANG Haibo，ZHANG Tianhong.Establishment and application in performance seeking control of an on-board adaptive aero-engine model considering nonlinear remainders[J].Journal of Propulsion Technology，2016，37（1）:172-180.（in Chinese）

[9]周文祥,黃金泉,竇建平.渦扇發動機部件級起動模型[J].航空動力學報,2006,21（2）:248-253.ZHOU Wenxiang,HUANG Jinquan,DOU Jianping.Development of component-level startup model for a turbofan engine[J].Journal of Aerospace Power.2006,21（2）:248-253.（in Chinese）

[10]周文祥,黃金泉.渦扇發動機高空起動模型研究[J].航空動力學報,2007,22（8）:1384-1390.ZHOU Wenxiang,HUANG Jinquan.Research on the startup model of turbofan engine at high altitude[J].Journal of Aerospace Power,2007,22（8）:1384-1390.（in Chinese）

[11]施洋,屠秋野,嚴紅明,等.一種航空發動機全狀態性能模型[J].航空動力學報,2017,32（2）:373-381.SHI Yang,TU Qiuye,YAN Hongming,et al.A full states performance model for aero engine[J].Journal of Aerospace Power,2007,22（8）:1384-1390.（in Chinese）

[12]AIAA.A process for the creation of T-MATS propulsion system models from NPSS data[C]//.AIAA Joint Propulsion Conference.Cleveland,OH.2014.

[13]胡曉煜.世界中小型航空發動機手冊[M].北京.航空工業出版社,2006:443-445.HU Xiaoyu.World minitype aeroengine manual[M]Beijing:Aviation Industry Press,2006:443-445.（ in Chinese）

[14]楊帆,樊丁,彭凱,等.基于試車數據的航空發動機起動過程建模[J].空軍工程大學學報:自然科學版,2013（6）:1-4.YANG Fan,FAN Ding,PENG Kai,et al.Modeling of an aero-engine in starting process based on the test data[J].Journal of Air Force Engineering University:Natural Science Edition,2013（6）:1-4.（in Chinese）

[15]Sellers J F,Daniele C J.Dyngen:a program for calculating steady-state and transient performance of turbojet and turbofan engines[R].NASA-TN－D-7901,1975.

[16]Daniele C J,Krosel S M,Szuch J R,et al.Digital computer program for generating dynamic turbofan engine models [R].NASA-TM-83446,1983.

[17]李慶揚,莫孜中,祁力群.非線性方程組的數值解法[M].北京：科學出版社,1999:108-113.LI Qingyang,MO Zizhong,QI Liqun.Numerical solutions of nonlinear equations[M].Beijing：Science Press,1999:108-113.（in Chinese）

[18]Claus R W,Evans A L,Lylte J K,et al.Numerical propulsion system simulation [J].Computing Systems in Engineering,1991,2（4）:357-364.

[19]Visser W P J,Broomhead M J.GSP,a generic object-oriented gas turbine simulation environment[C]//ASME Turbo Expo 2000:Power for Land,Sea,and Air.2000:V001T01A002.

[20]Broyden C G.Quasi-Newton methods and their application to function minimization[J].Mathematics of Computation,1967,21（99）:368-381.