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

朱閃閃，曹永華，謝中敏

（江蘇航空職業(yè)技術(shù)學(xué)院航空工程學(xué)院，江蘇 鎮(zhèn)江 212134）

0 引言

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

1 發(fā)動機模型的建立

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

1.1 發(fā)動機結(jié)構(gòu)與建模假設(shè)

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

圖1 DGEN380 發(fā)動機結(jié)構(gòu)簡圖

建立航空發(fā)動機模型時，作以下假設(shè)：

（1）將發(fā)動機內(nèi)的氣體流動視為一維流動；

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

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

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

1.2 DGEN380 發(fā)動機部件級模型的建立

（1）初始變量模塊

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

（2）減速齒輪箱模塊

根據(jù)理論公式，如果當(dāng)前低壓軸的物理轉(zhuǎn)速為nl，那么風(fēng)扇的物理轉(zhuǎn)速為

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

該模塊包括風(fēng)扇的特性數(shù)據(jù)，因此要根據(jù)當(dāng)前風(fēng)扇的換算轉(zhuǎn)速和壓比來求得風(fēng)扇當(dāng)前的效率和換算流量[7]，即：

從而根據(jù)式（4）得到風(fēng)扇實際的效率和流量：

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

在simulink 環(huán)境下建立好的DGEN380 發(fā)動機的部件級模型如圖2 所示。

圖2 DGEN380 發(fā)動機部件級模型

2 發(fā)動機穩(wěn)態(tài)模型的建立與求解

2.1 穩(wěn)態(tài)模型的建立

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

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

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

求解上式（10）即是對發(fā)動機穩(wěn)態(tài)模型求解，從而得到各工作參數(shù)。

2.2 穩(wěn)態(tài)模型的求解

航空發(fā)動機穩(wěn)態(tài)模型求解一般采用數(shù)值法。如N-R法[9]、改進的牛頓法[10]等。其中N-R 法最為常用，本文選取N-R 法對發(fā)動機的穩(wěn)態(tài)模型進行求解。

根據(jù)上式（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為雅可比矩陣，其表達(dá)式為：

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

3 穩(wěn)態(tài)求解中可能存在的問題分析

3.1 迭代中斷

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

（1）插值過程無法進行

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

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

（2）迭代過程中出現(xiàn)奇異

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

如果上式的溫度出現(xiàn)負(fù)值的情況，將會導(dǎo)致計算無法進行，迭代中斷。出現(xiàn)這種奇異一般是獨立變量設(shè)置不當(dāng)造成的。表1 是在發(fā)動機設(shè)計點處設(shè)置不同的獨立變量對計算結(jié)果的影響。

表1 初值對計算結(jié)果的影響

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

3.2 迭代不收斂

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

（1）獨立變量初值設(shè)置的不合理

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

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

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

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

4 仿真算例及結(jié)果分析

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

表2 油門桿位置對應(yīng)的燃油流量

表3 不同方法對求解結(jié)果的影響

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

5 結(jié)論

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