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

王希影

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

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

王希影

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

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

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

1 引言

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

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

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

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

2 P1有限元法基本理論推導

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

邊界條件為[7]：

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

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

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

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

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

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

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

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

式中：

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

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

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

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

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

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

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

如圖2所示，一個二維正方形半透明灰體介質(zhì)被發(fā)射率為0.8的不透明邊界包圍，底面溫度1 000 K，其它壁面溫度500 K，光學厚度τ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)及文獻[6]中DOM-FEM的計算結(jié)果進行比較(圖5)。結(jié)果顯示，P1-FEM與P1-FDM、DOM-FEM的計算結(jié)果規(guī)律一致，說明了P1-FEM程序的正確性；P1-FEM與DOM-FEM(文獻[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可以很好地應用于實際工程計算。

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

4 結(jié)束語

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

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

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

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

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

[5]程強，周懷春，黃志峰，等.DRESOR法對二維輻射傳輸傳遞問題研究[J].工程熱物理學報，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]孔祥謙.有限單元法在傳熱學中的應用[M].北京：科學出版社，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]謝仲生，鄧力.中子輸運理論數(shù)值計算方法[M]．西安：西北工業(yè)大學出版社，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ù)求解輻射傳輸方程的假散射和射線效應[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ā)動機傳熱和燃燒研究。