基于蒙特卡羅法的航空發(fā)動機穩(wěn)態(tài)空氣系統(tǒng)網(wǎng)絡(luò)計算

胡肖肖，李育隆，吳 宏

（1.中國艦船研究設(shè)計中心，武漢430064；2.清華大學(xué)航天航空學(xué)院，北京100084；3.北京航空航天大學(xué)能源與動力工程學(xué)院，北京100191）

0 引言

在航空發(fā)動機空氣系統(tǒng)設(shè)計中，穩(wěn)態(tài)空氣系統(tǒng)的參數(shù)值作為空氣冷卻系統(tǒng)設(shè)計中直接參考的關(guān)鍵參數(shù)，關(guān)系到空氣系統(tǒng)設(shè)計是否有效和發(fā)動機能否安全可靠地工作，有著重要的研究價值。復(fù)雜的空氣系統(tǒng)的影響因素較多，幾乎不可能對其進行嚴格的幾何求解。但可將穩(wěn)態(tài)下空氣在二次空氣系統(tǒng)內(nèi)的流動簡單看作1維不可壓縮流動，并離散為由不同元件和節(jié)點組成的網(wǎng)絡(luò)系統(tǒng)[1-2]，進而開展計算。目前，在穩(wěn)態(tài)空氣系統(tǒng)的計算研究中，常使用流體網(wǎng)絡(luò)法[3-6]，即數(shù)值求解網(wǎng)絡(luò)系統(tǒng)的連續(xù)性方程組和能量方程組。流體網(wǎng)絡(luò)算法相對成熟，但求解過程是否收斂、收斂速度及計算精度在很大程度上依賴于研究人員的求解經(jīng)驗，并且在求解復(fù)雜網(wǎng)絡(luò)和復(fù)雜部件時難以獲得收斂解。目前，空氣系統(tǒng)穩(wěn)態(tài)算法主要用于航空發(fā)動機典型設(shè)計點狀態(tài)的計算。

經(jīng)過文獻調(diào)研發(fā)現(xiàn)，流體網(wǎng)絡(luò)的研究不僅存在于航空發(fā)動機空氣系統(tǒng)領(lǐng)域，在市政管網(wǎng)領(lǐng)域也有對供水及供燃氣管網(wǎng)中1維不可壓縮穩(wěn)態(tài)流動的研究，其中蒙特卡羅法的應(yīng)用引起了研究人員關(guān)注。例如，吉慶豐等[7]將蒙特卡羅法應(yīng)用于供水管網(wǎng)分析，建立了相應(yīng)的隨機游動模型，結(jié)果表明：蒙特卡羅法計算簡便、靈活，所需計算機內(nèi)存少，便于對多水源的復(fù)雜管網(wǎng)進行水力計算。白建輝等[8]為求解天然氣管網(wǎng)節(jié)點的數(shù)學(xué)模型方程組，將蒙特卡羅方法的計算結(jié)果作為牛頓法節(jié)點壓力的初始值，高效率地實現(xiàn)了天然氣管網(wǎng)穩(wěn)態(tài)分析。結(jié)果表明：對于同一算例，用改進算法的計算效率比使用單一算法的明顯提高。

為了解決目前穩(wěn)態(tài)空氣系統(tǒng)算法在求解復(fù)雜網(wǎng)絡(luò)時不容易獲得收斂解的問題，本文嘗試引入蒙特卡羅方法[9-10]對航空發(fā)動機穩(wěn)態(tài)空氣系統(tǒng)網(wǎng)絡(luò)進行計算，并用商業(yè)軟件flowmaster驗證了計算結(jié)果的可靠性。為航空發(fā)動機穩(wěn)態(tài)空氣系統(tǒng)網(wǎng)絡(luò)計算提供了1種新方法，并為進一步研究及應(yīng)用提供了參考。

1 計算方程及方法

1.1 方法原理

蒙特卡羅方法亦稱為隨機模擬[10]，基本思想是建立1個概率模型或隨機過程，使其參數(shù)等于問題的解。通過對模型/過程的觀察或抽樣試驗來計算所求參數(shù)的統(tǒng)計特征，最后給出解的近似值。本文以概率的思想表示流量的分配關(guān)系，通過模擬大量粒子的不斷遷移來實現(xiàn)系統(tǒng)內(nèi)的壓力分配，從而求出空氣系統(tǒng)內(nèi)節(jié)點的壓力，再計算流經(jīng)各元件的流量。

1.2 基本方程

計算航空發(fā)動機空氣系統(tǒng)時，依據(jù)氣體動力學(xué)和傳熱學(xué)基礎(chǔ)理論，一般將冷卻空氣系統(tǒng)抽象為由節(jié)流單元與腔室組成的流體網(wǎng)絡(luò)。對于計算中的流熱耦合[11-13]（大多是弱耦合[14]），可將溫度作為已知條件，分別求解流體網(wǎng)絡(luò)和熱網(wǎng)絡(luò)。

現(xiàn)采用蒙特卡羅法求解溫度條件已知的流體網(wǎng)絡(luò)。定義空氣系統(tǒng)內(nèi)任一節(jié)點m 處流進節(jié)點的流量為正，流出節(jié)點的流量為負，則節(jié)點處質(zhì)量流量總和為0。連續(xù)性方程為

式中：mk為第k 段流動分支的流量。能量方程為

式中：Pin和Pout分別為第k 段分支兩端節(jié)點處進、出口壓力；ρ 為管內(nèi)流動介質(zhì)的密度；L、Dh、A 分別為管段的長度、直徑、橫截面積；ε 為管段的阻力系數(shù)（ε 的大小與雷諾數(shù)和粗糙度有關(guān)）；ω 為旋轉(zhuǎn)角速度；rin、rout分別為進、出口半徑。

簡單流體網(wǎng)絡(luò)如圖1所示，設(shè)節(jié)點Q1是流體網(wǎng)絡(luò)中任意1個節(jié)點，與之相連的分支管路有n個，相鄰的n個節(jié)點分別表示為Q11、Q12、…Q1n。

圖1 簡單流體網(wǎng)絡(luò)

聯(lián)立式（1）和式（3），可以化簡得節(jié)點Q1處壓力P（Q1）的方程為

1.3 計算方法

求解空氣系統(tǒng)問題即根據(jù)系統(tǒng)的流阻特性和幾何特性，對空氣系統(tǒng)進行流量分配和壓力分配。把流體看作大量粒子組成的集合體，粒子從進口依次經(jīng)過不同節(jié)點流至出口，遇到分叉處，粒子以不同的比例分別流入下游各節(jié)點。這個過程可以用概率的思想描述為：粒子以一定的概率遷移到相鄰節(jié)點，不斷遷移直至移動到壓力已知點，從而實現(xiàn)全場的壓力分配。其中，概率的大小由節(jié)點間元件的屬性決定，也就是式（5）中的轉(zhuǎn)移概率。

引入蒙特卡羅方法進行求解，求解過程為：構(gòu)造隨機游動模型，進行大量重復(fù)試驗，可以得到不同的試驗值，取其平均值作為待求節(jié)點的壓力值，具體用Fortran語言編程實現(xiàn)。例如求解空氣系統(tǒng)網(wǎng)絡(luò)的某一節(jié)點Q1的壓力值時，構(gòu)造以下隨機游動模型：某個粒子自Q1節(jié)點開始游動，按照轉(zhuǎn)移概率x1（i）向與Q1節(jié)點相鄰的節(jié)點隨機移動一步，若粒子第一步到達的節(jié)點是Q1m，再按轉(zhuǎn)移概率x1m（i）向與Q1m節(jié)點相鄰的節(jié)點隨機移動一步，如此重復(fù)下去，直至粒子遷移到壓力已知節(jié)點結(jié)束。假設(shè)共有N個粒子進行N 次隨機移動，就會獲得N個隨機變量值ε1，ε2，…εN，其平均值便可作為P（Q1）的近似值[15]。

2 計算算例與分析

2.1 計算算例

以如圖2所示的某空氣系統(tǒng)網(wǎng)絡(luò)為例，元件1與進口相連，進口處壓力為P1，元件6、10、17分別與3個出口相連，出口壓力分別為P7、P11、P18。圖中將不同的元件進行編號，已知流體進、出口的壓力值，求旋轉(zhuǎn)狀態(tài)下空氣系統(tǒng)中各節(jié)點的壓力值。

圖2 某空氣系統(tǒng)網(wǎng)絡(luò)

空氣系統(tǒng)內(nèi)包含不同的元件，如不同形式的管、封嚴壁齒、小孔、腔等。為驗證蒙特卡羅方法的可行性，把所有元件都當(dāng)作阻力元件，只要給定阻力系數(shù)的值，就可以確定壓力和流量的關(guān)系。

管道、孔、腔室滿足式（2），可將其轉(zhuǎn)化為如式（4）所示的壓力和流量的線性關(guān)系式。

對于封嚴篦齒結(jié)構(gòu)

式中：N 為篦齒齒數(shù)；F 為側(cè)面積；b 為經(jīng)驗系數(shù)；R 為理想氣體常數(shù)；T 為溫度。將封嚴篦齒公式化為線性形式

旋轉(zhuǎn)阻力系數(shù)為

式中：u 為旋轉(zhuǎn)切速度；vax為孔口的軸向流速；Cd為流量系數(shù)；Cdrot為旋轉(zhuǎn)阻力系數(shù)；k 為阻力系數(shù)。

2.2 計算結(jié)果與分析

分別設(shè)定3組不同工況，即不同的進出口壓力與旋轉(zhuǎn)速度，計算得到空氣系統(tǒng)各節(jié)點的壓力，3組工況及計算結(jié)果如下：

（1）工況1：T=419.2K，P1=52619Pa，P7=52093Pa，P11=49408Pa，P18=20332Pa，n=0r/min。節(jié)點的壓力計算結(jié)果如圖3（a）所示。

（2）工況2：T=435.0K，P1=57488Pa，P7=56913Pa，P11=54718Pa，P18=22295Pa，n=20900r/min。節(jié)點壓力計算結(jié)果如圖3（b）所示。

（3）工況3：T=450.2K，P1=62077Pa，P7=61457Pa，P11=58434Pa，P18=24348Pa，n=22000r/min。節(jié)點壓力計算結(jié)果如圖3（c）所示。

圖3 蒙特卡羅法與flowmaster壓力計算結(jié)果比較

蒙特卡羅方法通過大量的簡單重復(fù)抽樣進行計算，程序簡單，不會出現(xiàn)不收斂的情況。計算發(fā)現(xiàn)，經(jīng)過較少次游動，所得的壓力值已經(jīng)比較接近真實值。

從圖3中可見，采用蒙特卡羅方法得到的空氣系統(tǒng)的各節(jié)點壓力值與參考結(jié)果相差較小，最大誤差僅為0.628%。準(zhǔn)確的壓力計算可為空氣系統(tǒng)設(shè)計提供重要依據(jù)，運用蒙特卡羅方法求得各節(jié)點的壓力值后，再根據(jù)流量與元件兩端壓力值的關(guān)系，計算流經(jīng)各元件的流體流量。

在工況1條件下，根據(jù)如圖3（a）所示的壓力值計算出流經(jīng)各元件的流量值，如圖4所示。從圖4中可見，流經(jīng)各元件的流量與flowmaster的結(jié)果基本一致，但流經(jīng)元件12、15的流量與flowmaster計算結(jié)果有較大誤差，最大為0.00367kg/s。這是由于蒙特卡羅方法是隨機算法，其收斂是概率意義下的收斂,不能解決精確度要求較高的問題。在計算過程中，以壓力為控制變量，經(jīng)過不斷迭代，得出誤差較小的壓力值，但無法避免誤差。由于壓力變化會引起較大的流量變化，因此流量計算結(jié)果也會有一定誤差。

圖4 蒙特卡羅法與flowmaster流量計算結(jié)果比較

3 結(jié)論

通過大量的簡單重復(fù)抽樣可實現(xiàn)運用蒙特卡羅方法求解穩(wěn)態(tài)空氣系統(tǒng)，該方法思路清晰、使用簡便，對初值的設(shè)定不敏感，不會出現(xiàn)無法求解復(fù)雜空氣系統(tǒng)的情況。在求解過程中發(fā)現(xiàn)，經(jīng)過數(shù)次游動即可得到比較接近壓力值的解。用flowmaster軟件對該計算方法進行驗證，壓力計算結(jié)果最大偏差為0.628%，說明2種方法的計算結(jié)果吻合較好，也初步驗證了蒙特卡羅方法求解航空發(fā)動機穩(wěn)態(tài)空氣系統(tǒng)網(wǎng)絡(luò)的可行性。

根據(jù)壓力值計算出的流量有一定誤差，結(jié)合其他方法（如牛頓-拉夫遜法），有可能提高計算精度，這需要進一步研究優(yōu)化。

[1]陸海鷹，楊燕生，王鳴.航空發(fā)動機空氣系統(tǒng)特性的數(shù)值模擬[J].航空發(fā)動機，1997（1）：6-13.LU Haiying,YANG Yansheng,WANG Ming.The numerical simulation of the aeroengine air system [J].Aeroengine,1997（1）:6-13.(in Chinese)

[2]左承基，傅秋陽，歐陽明高，等.壓縮空氣發(fā)動機系統(tǒng)的數(shù)值模擬[J].合肥工業(yè)大學(xué)學(xué)報：自然科學(xué)版，2005，28（9）：985-988.ZUO Chengji,FU Qiuyang,OUYANG Minggao,et al.Numerical simulation of the compressed air engine system[J].Journal of Hefei University of Technology:Natural Science,2005,28（9）:985-988.(in Chinese)

[3]Majumdar A,Steadman T.Numerical modeling of pressurization of a propellant tank[J].Journal of Propulsion and Power,2001,17(2):385-390.

[4]Kritpiphat W,Tontiwachwuthikul P,Chan C.W.Pipeline network modeling and simulation for intelligent monitoring and control:a case study of a municipal water supply system[J].Industrial&Engineering Chemistry Research,1998,37(3):1033-1044.

[5]Obradovic D.Modelling of demand and losses in reallife water distribution systems[J].Urban Water,2000,2(2):131-139.

[6]陶智，侯升平，韓樹軍，等.流體網(wǎng)絡(luò)法在發(fā)動機空氣冷卻系統(tǒng)設(shè)計中的應(yīng)用[J].航空動力學(xué)報，2009，24（1）：1-6.TAO Zhi,HOU Shengping,HAN Shujun,et al.Study on application of fluid network into the design of air system in engine[J].Journal of Aerospace Power,2009，24（1）:1-6.(in Chinese)

[7]吉慶豐，周濟人.應(yīng)用蒙特卡羅法分析供水管網(wǎng)[J].灌溉排水，1998，17（4）：16-20.JI Qingfeng,ZHOU Jiren.Pipe network analysis with Monte Carlo method[J].Irrigation and Drainage,1998,17（4）:16-20.(in Chinese)

[8]白建輝，汪玉春，郜峰，等.天然氣管網(wǎng)穩(wěn)態(tài)分析綜合方法[J].油氣儲運，2009，28（2）：37-39.BAI Jianhui,WANG Yuchun,GAO Feng,et al.Combination based algorithm for steady analysis of natural gas pipeline network[J].Qil&Gas Storage and Transportation,2009,28（2）:37-39.(in Chinese)

[9]徐鐘濟.蒙特卡羅方法[M].上海：上海科學(xué)技術(shù)出版社，1985：2-18.XU Zhongji.Monte Carlo method[M].Shanghai:Shanghai scientific&Technical Publishers,1985：2-18.(in Chinese)

[10]朱本仁.蒙特卡洛方法引論[M].濟南：山東大學(xué)出版社，1987：30-42.ZHU Benren.Introduce of Monte Carlo method[M].Jinan:Shandong University Press,1987：30-42.(in Chinese)

[11]朱惠人，白洛林，張麗，等.多孔隔熱壁溫度場氣動傳熱耦合計算方法[J].航空動力學(xué)報，2003，18（2）：230-234.ZHU Huiren,BAI Luolin,ZHANG Li,et al.Coupled method of flowheat transfer to calculate temperature field of multihole wall[J].Journal of Aerospace Power,2003,18（2）:230-234.(in Chinese)

[12]Miltiadis V P.Feed back control of thermal systems modeled via the network approach[J].Journal of Dynamic Systems,Measurement,and Control,2004,126（3）:509-519.

[13]Athavale M M,Ho Y H.Analysis of coupled seals,secondary and power stream flow fields in aircraft and aerospace turbomachines[R].NASA-CR-2005-212716.

[14]吳麗君，吳宏.流熱網(wǎng)絡(luò)耦合及局部3維計算方法研究[J].科學(xué)技術(shù)與工程，2011，11（3）：521-527.WU Lijun,WU Hong.Study of thermal flow network and local three dimensional calculation method[J].Science Technology and Engineering,2011,11（3）:521-527.(in Chinese)

[15]Spiegel M R,Schiller J J,Srinivasan R A.Probability and statistics[M].New York:McGraw-Hill,2009：23-74.