小型渦扇發動機特性計算收斂性方法研究

朱閃閃，曹永華，謝中敏

（江蘇航空職業技術學院航空工程學院，江蘇 鎮江 212134）

0 引言

航空發動機是一種復雜的裝置，它的研制水平代表了一個國家的工業、科技以及國防實力的水平，目前，只有美國、英國、中國等少有的幾個國家具有這樣的自主研發能力[1]。其余各個國家的航空界正在努力地尋找一種效率高、成本低、風險低的解決措施。在這一解決措施中，計算機仿真具有很大的優勢，通過計算機仿真可以有效地縮短研制周期，減少研制成本[2]。陳玉春等[3]提出了變步長的牛頓-拉夫遜法對航空發動機模型進行求解，但是這種方法僅在單軸渦噴發動機上得到了應用，不具備普遍性。施洋等[4]提出罰函數法等優化方法來代替發動機非線性方程組的迭代求解過程，從而完成航空發動機性能的仿真模擬計算。這些方法對計算初值沒有要求，在一定程度上緩解了計算過程中的發散問題，但是其計算速度非常緩慢，在實際應用中適用性有限。所以本文首先利用Matlab/Simulink 軟件對DGEN380 發動機進行建模仿真，然后利用牛頓-拉夫遜方法對穩態模型進行求解，針對該算法中可能存在結果不收斂的原因進行分析總結，最后對穩態求解過程中的收斂性進行仿真分析研究。

1 發動機模型的建立

目前發動機的建模大體可以分為機理建模和實驗建模[5，6]，機理建模是以發動機運行過程中所遵循的熱力學過程所建立的數學模型，而試驗建模（系統辨識）是根據發動機的實驗數據進行分析而獲得的數學模型，對實驗數據有著較高的要求。而機理建模過程簡便，所得到的模型精度也較高，并且可以在發動機運行全包線內對發動機穩態、動態特性進行分析。這種方法廣泛地被國內外研究人員所使用。本文將結合DGEN380 發動機設計點實驗數據建立該發動機的穩態模型。

1.1 發動機結構與建模假設

本文的建模對象是DGEN380 的齒輪傳動風扇發動機（GTF），按照其工作原理主要劃分為進氣道、風扇、齒輪箱[6]、壓氣機等部件。其結構劃分簡圖如圖1 所示。

圖1 DGEN380 發動機結構簡圖

建立航空發動機模型時，作以下假設：

（1）將發動機內的氣體流動視為一維流動；

（2）忽略發動機放氣和引起的影響；

（3）不考慮容積慣性和燃燒延遲的影響；

（4）各個部件中的比熱值取定值。

1.2 DGEN380 發動機部件級模型的建立

（1）初始變量模塊

該模塊主要是由DGEN380 發動機的高、低壓轉速以及四個旋轉部件的壓比組成的初值變量。

（2）減速齒輪箱模塊

根據理論公式，如果當前低壓軸的物理轉速為nl，那么風扇的物理轉速為

（4）流量和效率插值模塊

該模塊包括風扇的特性數據，因此要根據當前風扇的換算轉速和壓比來求得風扇當前的效率和換算流量[7]，即：

從而根據式（4）得到風扇實際的效率和流量：

式中k為空氣的絕熱指數，在此取1.4，R為氣體常數，在此取287.04 J/kg·K。

在simulink 環境下建立好的DGEN380 發動機的部件級模型如圖2 所示。

圖2 DGEN380 發動機部件級模型

2 發動機穩態模型的建立與求解

2.1 穩態模型的建立

航空發動機處于穩態時，需要滿足一定的穩態平衡條件，選取高、低壓轉速、風扇、壓氣機、高、低壓渦輪的壓比作為猜值變量，即[nL，πF，nH，πC，πHT，πLT]，故建立以下平衡方程[8]。

式中，ηhpm，ηlpm為高，低壓轉子的機械效率。分別取0.93和0.81。

根據前面所建立的發動機整機模型可知，但飛行條件和燃油流量已知并保持一定時，只需要知道6 個狀態參數[nL，πF，nH，πC，πHT，πLT]即可得到發動機當前所有的熱力學參數，所發動機的穩態模型可以用這六個參數來表示，即以下非線性方程：

求解上式（10）即是對發動機穩態模型求解，從而得到各工作參數。

2.2 穩態模型的求解

航空發動機穩態模型求解一般采用數值法。如N-R法[9]、改進的牛頓法[10]等。其中N-R 法最為常用，本文選取N-R 法對發動機的穩態模型進行求解。

根據上式（10），令F（x）= [f1，f2，f3，f4，f5，f6]，x= [nL，πF，nH，πC，πHT，πLT]T。N-R 的核心公式為x(k+1)=xk-J-1F（x(k)），其中xk為上一步初猜值的解，為xk+1迭代一步之后修正過后的解。J為雅可比矩陣，其表達式為：

采用向前差分法[11]對各個變量求偏導，對初猜解向量中的各項進行擾動求出雅可比矩陣各列，如第一次擾動nL，將擾動后的向量[nL+ △nL，πF，nH，πC，πHT，πLT]代入式（11），求出雅可比矩陣的第一列，依次類推，得到最終的雅可比矩陣。

3 穩態求解中可能存在的問題分析

3.1 迭代中斷

要保證迭代過程可以順利地進行，首先保證迭代過程不會出現中斷，而在仿真的過程中，經常會出現迭代中斷的問題，原因大致可以總結為以下兩點：

（1）插值過程無法進行

部件級模型都是基于部件特性來建立的，在整機模型中，風扇、壓氣機、高壓渦輪、低壓渦輪內的部件特性的插值能否順利地進行，決定著迭代是否完成。

以風扇為例，其部件特性的插值如上式（3）所示，如果在插值的過程中，插值的輸入量超界將會使插值計算過程無法進行，最終導致穩態求解無法完成。而產生插值中斷的原因主要是由于獨立變量設置的不合理，當設置的獨立變量相比于真實值偏差較大時，會使迭代計算求解過程中斷。

（2）迭代過程中出現奇異

迭代過程中出現計算奇異的情況主要是由于開負平方根造成的，例如風扇求解流量時：

如果上式的溫度出現負值的情況，將會導致計算無法進行，迭代中斷。出現這種奇異一般是獨立變量設置不當造成的。表1 是在發動機設計點處設置不同的獨立變量對計算結果的影響。

表1 初值對計算結果的影響

從以上可以看出，合適的獨立變量是進行迭代的第一步，所以對獨立變量限制的設置顯得至關重要。

3.2 迭代不收斂

大量的實驗數據表明，基本的N-R 法雖然完成了迭代的過程，但是大多數情況下結果都是不收斂的，即穩態共同工作方程的求解精度不滿足要求，這樣殘差會隨著迭代過程不斷地增大，導致最后的結果發散，這一結果的主要原因有：

（1）獨立變量初值設置的不合理

由于發動機模型的特殊性，在求解迭代過程中也許會使求解結果與現實狀態不符，如燃燒室的出口總溫小于進口總溫，這時的結果只具有數學意義而沒有物理意義的參數[10]。雖然完成了迭代過程，但最后結果是發散的。

（2）基本N-R 迭代算法本身的不足

N-R 方法的關鍵是牛頓步長的選取，一般情況下，牛頓步長選取是獨立變量的0.001 倍。而如果獨立變量的殘差較大時，擾動量較小將起不到什么作用，如果獨立變量的殘差較小，擾動量如果太大將會導致特性圖插值超界。所以如果采用固定步長，很容易導致最后結果的不收斂。

根據以上所述，導致發動機穩態求解最后結果不收斂的情況大致可以歸結為兩點，變量初值設置的不合適，二是迭代算法中牛頓步長選取有問題。限制獨立變量X 的初值的目的是為了保證獨立變量的初值不會導致部件插值出現外插值的情況，具體的做法是：例如風扇和壓氣機，根據特性圖中nccor的范圍，可以對獨立變量進行限制。對于高、低壓渦輪，根據特性圖中Wgcor和ntcor的范圍，可以對獨立變量進行限制：Wgcor如果超限，則對其進行限制；根據渦輪進口總溫和nc可以計算得到ntcor，如果ntcor超過范圍，則限制nc。其他部件的計算根據具體情況采取參數限制措施。獨立變量初值的限制，保證N-R 法計算的第一步能夠完成。而對于牛頓步長的選取，本文采用變牛頓步長法使得迭代過程能夠在保證精度的情況下順利地進行。以下根據這兩個問題進行仿真算例分析。

4 仿真算例及結果分析

根據以上建立的發動機穩態模型，對發動機的穩態求解過程進行仿真，本文選取飛行包線內比較有代表意義的5 個點，見表2。在基本牛頓-拉夫遜方法的基礎上，采用獨立變量初值的設置方法以及牛頓變步長的方法改善發動機的穩態求解的迭代過程。表3 為兩種方法的對比結果。

表2 油門桿位置對應的燃油流量

表3 不同方法對求解結果的影響

從以上結果不難看出，對獨立變量的限制可以加快迭代速度，改善迭代效果。迭代算法的改進（即在迭代過程中采用變步長），可以使迭代朝著更加有效的方向進行。所以以上方法不僅有效地改善了收斂性，還對收斂速度有了進一步的改善。

5 結論

本文對DGEN380 發動機部件特性計算中的不收斂問題進行了系統的研究后，以牛頓- 拉夫遜方法為基礎，對該方法進行了改進并取得了非常好的效果，改善了算法迭代收斂性的問題。認為可以不必使用其它計算量大的求解算法。同樣這些改進辦法也可以為其他類型的航空發動機的特性計算提供參考。