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

胡肖肖，李育隆，吳 宏

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

0 引言

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

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

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

1 計(jì)算方程及方法

1.1 方法原理

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

1.2 基本方程

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

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

式中：mk為第k 段流動(dòng)分支的流量。能量方程為

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

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

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

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

1.3 計(jì)算方法

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

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

2 計(jì)算算例與分析

2.1 計(jì)算算例

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

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

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

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

對(duì)于封嚴(yán)篦齒結(jié)構(gòu)

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

[1]陸海鷹，楊燕生，王鳴.航空發(fā)動(dòng)機(jī)空氣系統(tǒng)特性的數(shù)值模擬[J].航空發(fā)動(dòng)機(jī)，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]左承基，傅秋陽(yáng)，歐陽(yáng)明高，等.壓縮空氣發(fā)動(dòng)機(jī)系統(tǒng)的數(shù)值模擬[J].合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(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ā)動(dòng)機(jī)空氣冷卻系統(tǒng)設(shè)計(jì)中的應(yīng)用[J].航空動(dòng)力學(xué)報(bào)，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]吉慶豐，周濟(jì)人.應(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].油氣儲(chǔ)運(yùn)，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]徐鐘濟(jì).蒙特卡羅方法[M].上海：上海科學(xué)技術(shù)出版社，1985：2-18.XU Zhongji.Monte Carlo method[M].Shanghai:Shanghai scientific&Technical Publishers,1985：2-18.(in Chinese)

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

[11]朱惠人，白洛林，張麗，等.多孔隔熱壁溫度場(chǎng)氣動(dòng)傳熱耦合計(jì)算方法[J].航空動(dòng)力學(xué)報(bào)，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ì)算方法研究[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.