王希影

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

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

王希影

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

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

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

1 引言

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

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

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

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

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

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

邊界條件為[7]：

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

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

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

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

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

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

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

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

式中：

邊界單元控制方程變分計算基本公式可寫成：

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

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

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

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

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

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

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

不同網(wǎng)格數(shù)量下P1-FEM、P1-FDM與DOM-FEM的計算時間比較見表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為計算方便多采用正交四邊形網(wǎng)格，網(wǎng)格單元數(shù)分別為100、625和2 916。可見，P1-FDM的計算速度最快，但其最大缺點是不能計算非規(guī)則幾何體的輻射換熱，限制了其使用范圍；P1-FEM的計算速度快于DOM-FEM，且在網(wǎng)格數(shù)目較大時體現(xiàn)得更為明顯，說明P1-FEM可以很好地應(yīng)用于實際工程計算。

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

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

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

4 結(jié)束語

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

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

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

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

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

[5]程強(qiáng)，周懷春，黃志峰，等.DRESOR法對二維輻射傳輸傳遞問題研究[J].工程熱物理學(xué)報，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ù)值計算方法[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ù)求解輻射傳輸方程的假散射和射線效應(yīng)[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ā)動機(jī)傳熱和燃燒研究。