二維參與性介質(zhì)輻射換熱的球諧有限元法

王希影

(中航空天發(fā)動(dòng)機(jī)研究院有限公司，北京100028)

二維參與性介質(zhì)輻射換熱的球諧有限元法

王希影

(中航空天發(fā)動(dòng)機(jī)研究院有限公司，北京100028)

采用球諧函數(shù)與有限元相結(jié)合的方法(即球諧有限元法)，對(duì)二維介質(zhì)內(nèi)輻射換熱問(wèn)題進(jìn)行研究。以規(guī)則和非規(guī)則形狀幾何體輻射問(wèn)題為例，使用球諧有限元法可得到較準(zhǔn)確的計(jì)算結(jié)果，且與其它數(shù)值計(jì)算方法(如球諧有限差分法、離散坐標(biāo)有限元法等)相比具有以下優(yōu)勢(shì)：球諧有限元法方法簡(jiǎn)單，計(jì)算速度快，適用于大型工程應(yīng)用；球諧有限元法可用于計(jì)算各種形狀的幾何體，計(jì)算應(yīng)用范圍廣。

航空發(fā)動(dòng)機(jī)；燃燒室；球諧函數(shù)法；球諧有限元法；輻射換熱；參與性介質(zhì)

1 引言

航空發(fā)動(dòng)機(jī)燃燒室內(nèi)存在復(fù)雜的流動(dòng)與換熱過(guò)程，準(zhǔn)確模擬燃燒產(chǎn)物與發(fā)動(dòng)機(jī)內(nèi)部結(jié)構(gòu)的換熱過(guò)程，可為燃燒室內(nèi)壁等處于嚴(yán)酷工作條件下的結(jié)構(gòu)件的設(shè)計(jì)提供重要依據(jù)。燃燒產(chǎn)物中包含炭黑、二氧化碳?xì)怏w、水蒸氣等，這些介質(zhì)在高溫下具有非常強(qiáng)的輻射能力，研究這類(lèi)具有吸收、發(fā)射、散射能力的介質(zhì)內(nèi)的輻射換熱過(guò)程，通常需要求解輻射傳遞方程。但由于該方程本身及介質(zhì)輻射特性的復(fù)雜性，理論分析解僅存在于極少數(shù)簡(jiǎn)單或極限問(wèn)題。為此，國(guó)內(nèi)外學(xué)者發(fā)展了許多高效數(shù)值方法來(lái)求解參與性介質(zhì)內(nèi)的輻射傳遞問(wèn)題，如Monte-Carlo法、區(qū)域法(ZM)、離散坐標(biāo)法(DOM)、球諧函數(shù)法(或PN法)、有限元法(FEM)、DRESOR法及上述方法相結(jié)合使用的方法等[1～6]。其中球諧函數(shù)法采用將輻射強(qiáng)度隨空間位置和方向變化進(jìn)行變量分離的處理方法，利用球諧函數(shù)的正交性，將積分-微分形式的輻射傳遞方程轉(zhuǎn)化為相對(duì)簡(jiǎn)單的偏微分方程組，理論上提供了一種獲得任意高階(精度)的近似求解方法。P1近似相對(duì)簡(jiǎn)單，且與標(biāo)準(zhǔn)能量方程的求解方

法兼容性好，因而在實(shí)際問(wèn)題求解中應(yīng)用最廣。然而隨著階數(shù)N和幾何體維數(shù)的增加，其數(shù)學(xué)復(fù)雜性迅速增加[7]。

球諧函數(shù)法的優(yōu)點(diǎn)在于其不需要空間方向離散，且便于與其它空間離散方法(如FEM、DOM)結(jié)合使用，球諧有限元法即結(jié)合了球諧函數(shù)和有限元空間離散。由于有限元法可處理復(fù)雜幾何體[8]，因此球諧有限元法也可用來(lái)處理各種形狀的幾何體。球諧有限元法一般用于偶對(duì)稱(chēng)中子傳輸、核測(cè)井和醫(yī)療成像[9，10]，目前將其應(yīng)用在參與性介質(zhì)內(nèi)輻射傳輸計(jì)算的文獻(xiàn)較少。

本文首先對(duì)P1有限元法(P1-FEM)進(jìn)行理論推導(dǎo)，然后通過(guò)編制相應(yīng)計(jì)算程序，分別對(duì)二維規(guī)則和非規(guī)則幾何體內(nèi)參與性介質(zhì)的輻射傳輸問(wèn)題進(jìn)行數(shù)值計(jì)算。

2 P1有限元法基本理論推導(dǎo)

在發(fā)射、吸收、散射介質(zhì)內(nèi)，輻射傳遞方程為：

邊界條件為[7]：

將輻射傳輸方程簡(jiǎn)化為只與空間坐標(biāo)有關(guān)的二階微分方程，然后對(duì)空間坐標(biāo)采用有限元離散。令，對(duì)直角坐標(biāo)系下的P1方程即方程(2)進(jìn)行Galerkin變分，得：

將式(3)代入式(4)中可得：

對(duì)于二維問(wèn)題的有限元求解域的離散，有限元法可選擇不同的離散單元形狀[8]，如三角形、四邊形、矩形等，但考慮對(duì)復(fù)雜幾何形狀的適應(yīng)性，本文選擇了三角形單元，如圖1所示。

對(duì)于有限元形函數(shù)的選擇，可以有高次單元和等參單元。本文中選擇了較簡(jiǎn)單的自然坐標(biāo)系下的等參單元。對(duì)于二維情況，形函數(shù)可由下面公式給出[8]：

式中：Ni+Nj+Nk=1；Δ為單元面積；aα、bα、cα(α=i，j，m)的定義為[8]

形函數(shù)的導(dǎo)數(shù)有[8]：

求解域內(nèi)任意單元e的控制方程變分計(jì)算基本公式可寫(xiě)成：

對(duì)于等參單元，運(yùn)用式(6)～式(8)，進(jìn)而得到以下有限元單元矩陣：

式中：

邊界單元控制方程變分計(jì)算基本公式可寫(xiě)成：

這里多了一項(xiàng)邊界jm的線(xiàn)性積分項(xiàng)計(jì)算，對(duì)于邊界上的線(xiàn)性插值函數(shù)：

其中0≤g≤1且ds=sidg，可得：

從而得到邊界單元系數(shù)矩陣：

將式(11)和式(15)帶入式(10)進(jìn)行求解即可。

3 計(jì)算模型及結(jié)果

例1：規(guī)則形狀——正方形封腔

如圖2所示，一個(gè)二維正方形半透明灰體介質(zhì)被發(fā)射率為0.8的不透明邊界包圍，底面溫度1 000 K，其它壁面溫度500 K，光學(xué)厚度τL=1.0。吸收和散射系數(shù)均為0.5 m-1，求正方形封腔內(nèi)介質(zhì)的溫度分布。

二維正方形封腔的網(wǎng)格分布如圖3所示，共有648個(gè)三角形網(wǎng)格單元，總節(jié)點(diǎn)數(shù)為359，邊界單元數(shù)為68。采用P1-FEM計(jì)算得到的封腔內(nèi)的各向同性介質(zhì)溫度分布如圖4所示，并將沿直線(xiàn)y=0.782 5處的溫度分布，與采用P1-有限差分(P1-FDM)及文獻(xiàn)[6]中DOM-FEM的計(jì)算結(jié)果進(jìn)行比較(圖5)。結(jié)果顯示，P1-FEM與P1-FDM、DOM-FEM的計(jì)算結(jié)果規(guī)律一致，說(shuō)明了P1-FEM程序的正確性；P1-FEM與DOM-FEM(文獻(xiàn)[6]已證明DOM-FEM具有高精度，在此假定其計(jì)算結(jié)果精確)溫度場(chǎng)計(jì)算的最大誤差為6.53%，且主要來(lái)自于球諧函數(shù)展開(kāi)部分的計(jì)算。

不同網(wǎng)格數(shù)量下P1-FEM、P1-FDM與DOM-FEM的計(jì)算時(shí)間比較見(jiàn)表1(Intel(R)Core(TM)i7-2600：8G內(nèi)存)。其中P1-FEM和DOM-FEM采用圖3所示的三角形網(wǎng)格，網(wǎng)格單元數(shù)分別為100、648和2 916；而P1-FDM為計(jì)算方便多采用正交四邊形網(wǎng)格，網(wǎng)格單元數(shù)分別為100、625和2 916。可見(jiàn)，P1-FDM的計(jì)算速度最快，但其最大缺點(diǎn)是不能計(jì)算非規(guī)則幾何體的輻射換熱，限制了其使用范圍；P1-FEM的計(jì)算速度快于DOM-FEM，且在網(wǎng)格數(shù)目較大時(shí)體現(xiàn)得更為明顯，說(shuō)明P1-FEM可以很好地應(yīng)用于實(shí)際工程計(jì)算。

例2：非規(guī)則形狀——非規(guī)則四邊形

圖6所示的非規(guī)則四邊形，左壁面溫度較高為1 000 K，其它壁面為300 K，介質(zhì)溫度也為300 K，壁面為黑體壁面。介質(zhì)為反照率0.5的吸收、各向同性散射介質(zhì)。衰減系數(shù)為1.0 m-1，求底面無(wú)量綱熱流密度分布。

計(jì)算結(jié)果如圖7所示，可見(jiàn)P1-FEM在計(jì)算無(wú)量綱熱流分布時(shí)與文獻(xiàn)[11]的結(jié)果有一定差距，這主要是由于球諧函數(shù)截?cái)嗾`差所致。球諧函數(shù)產(chǎn)生的誤差分為假散射和射線(xiàn)效應(yīng)兩類(lèi)，且都隨著PN法中N值的不斷增大而減小[12]。但本文并不推薦通過(guò)增大N值來(lái)減小誤差，原因?yàn)楫?dāng)N增大時(shí)會(huì)使數(shù)學(xué)復(fù)雜性增大。如當(dāng)N=3時(shí)，控制方程將由方程(2)一個(gè)二階偏微分方程變?yōu)橛伤膫€(gè)二階偏微分方程組成的方程組，邊界條件方程也由方程(3)一個(gè)一階偏微分方程變?yōu)橛伤膫€(gè)一階偏微分方程組成的方程組[13]，編程計(jì)算過(guò)程復(fù)雜，計(jì)算速度很慢。

4 結(jié)束語(yǔ)

本文對(duì)用于求解輻射傳輸?shù)那蛑C有限元法進(jìn)行了理論推導(dǎo)，并采用P1-FEM對(duì)規(guī)則形狀和不規(guī)則形狀幾何體內(nèi)的輻射傳輸進(jìn)行了計(jì)算。結(jié)果表明，球諧有限元法原理簡(jiǎn)單，計(jì)算速度快、精度較好，適用于解決工程問(wèn)題；與P1方法相比，球諧有限元法可用于計(jì)算各種形狀的幾何體，拓寬了計(jì)算應(yīng)用范圍。因此，工程上在計(jì)算精度要求不太高的情況下，推薦使用P1-FEM求解復(fù)雜形狀區(qū)域內(nèi)的輻射換熱問(wèn)題。

[1]談和平，夏新林，劉林華，等.紅外輻射特性與傳輸?shù)臄?shù)值計(jì)算——計(jì)算熱輻射學(xué)[M].哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2006：139—190.

[2]張建強(qiáng)，朱谷君.燃燒室中輻射熱流分布的蒙特卡羅計(jì)算[J].航空動(dòng)力學(xué)報(bào)，1999，14(3)：251—254.

[3]毛文懿，林宇震，許全宏，等.運(yùn)用區(qū)域法模型計(jì)算燃燒室內(nèi)一維輻射換熱[J].航空動(dòng)力學(xué)報(bào)，2010(3)：515—520.

[4]齊宏，阮立明，譚建宇.矩形介質(zhì)內(nèi)輻射換熱的有限元法[J].計(jì)算物理，2004，24(6)：547—549.

[5]程強(qiáng)，周懷春，黃志峰，等.DRESOR法對(duì)二維輻射傳輸傳遞問(wèn)題研究[J].工程熱物理學(xué)報(bào)，2008，29(10)：1735—1738.

[6]An W，Ruan L M，Qi H，et al.Finite Element Method for Radiative Heat Transfer in Absorbing and Anisotropic Scattering Media[J].Journal of Quantitative Spectroscopy and Radiative Transfer，2005，96：409—422.

[7]Modest M F.Radiative Heat Transfer[M].2nd ed.New York：Academic Press，2003：465—483.

[8]孔祥謙.有限單元法在傳熱學(xué)中的應(yīng)用[M].北京：科學(xué)出版社，1998.

[9]De Oliveira C R E.An Arbitrary Geometry Finite Element Method for Multigroup Neutron Transport with Anisotropic Scattering[J].Prog.Nucl.Energy，1996，18：227—236.

[10]謝仲生，鄧力.中子輸運(yùn)理論數(shù)值計(jì)算方法[M]．西安：西北工業(yè)大學(xué)出版社，2005：64—128.

[11]Baek S W，Byun D Y，Kang S J.The Combined Mon?te-Carlo and Finite Volume Method for Radiation in a Two-Dimensional Irregular Geometry[J].J Heat Mass Transfer，2000，43：2337—2344.

[12]張昊春，易洪亮，談和平.球諧函數(shù)求解輻射傳輸方程的假散射和射線(xiàn)效應(yīng)[J]．計(jì)算物理，2006，23(2)：237—242.

[13]Yang J，Modest M F.Elliptic PDE Formulation of General Three-DimensionalHigher-OrderPN-Approximations for Radiative Transfer[J].Journal of Quantitative Spectros?copy and Radiative Transfer，2006，104：217—227.

Spherical Harmonics Finite Element Method for Radiative Heat Transfer in Two Dimensional Participating Media

WANG Xi-ying
(China Aviation Engine Establishment，Beijing 100028，China)

The spherical harmonics finite element method,which combines the spherical harmonics method and the finite element method together,has been applied to deal with two dimensional radiative transfer prob?lems.Taking the regular and irregular geometry cases as examples,better results are obtained with the spheri?cal harmonics finite element method.Compared with other numerical methods,such as the spherical harmon?ics finite difference method and the discrete-ordinate finite-element method,the spherical harmonics finite el?ement method has some advantages.Firstly,the method is simple and the calculation speed is fast,and it is ap?propriate for engineering application.Secondly,the spherical harmonics finite element method can be applied to all kinds of geometry shapes,the range of application is wide.

aero-engine；combustor；pherical harmonics method；spherical harmonics finite element method；radiative heat transfer；participating media

V231.1；TK124

：A

：1672-2620(2014)02-0010-04

2013-09-24；

：2014-01-20

王希影(1982-)，黑龍江大慶人，工程師，博士，主要從事航空發(fā)動(dòng)機(jī)傳熱和燃燒研究。