朱閃閃，曹永華，謝中敏

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

0 引言

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

1 發(fā)動(dòng)機(jī)模型的建立

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

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

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

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

建立航空發(fā)動(dòng)機(jī)模型時(shí)，作以下假設(shè)：

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

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

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

（4）各個(gè)部件中的比熱值取定值。

1.2 DGEN380 發(fā)動(dòng)機(jī)部件級(jí)模型的建立

（1）初始變量模塊

該模塊主要是由DGEN380 發(fā)動(dòng)機(jī)的高、低壓轉(zhuǎn)速以及四個(gè)旋轉(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)扇實(shí)際的效率和流量：

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

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

圖2 DGEN380 發(fā)動(dòng)機(jī)部件級(jí)模型

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

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

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

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

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

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

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

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

根據(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]對各個(gè)變量求偏導(dǎo)，對初猜解向量中的各項(xiàng)進(jìn)行擾動(dòng)求出雅可比矩陣各列，如第一次擾動(dòng)nL，將擾動(dòng)后的向量[nL+ △nL，πF，nH，πC，πHT，πLT]代入式（11），求出雅可比矩陣的第一列，依次類推，得到最終的雅可比矩陣。

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

3.1 迭代中斷

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

（1）插值過程無法進(jìn)行

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

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

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

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

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

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

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

3.2 迭代不收斂

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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