二维参与性介质辐射换热的球谐有限元法

王希影

(中航空天发动机研究院有限公司，北京100028)

二维参与性介质辐射换热的球谐有限元法

王希影

(中航空天发动机研究院有限公司，北京100028)

采用球谐函数与有限元相结合的方法(即球谐有限元法)，对二维介质内辐射换热问题进行研究。以规则和非规则形状几何体辐射问题为例，使用球谐有限元法可得到较准确的计算结果，且与其它数值计算方法(如球谐有限差分法、离散坐标有限元法等)相比具有以下优势：球谐有限元法方法简单，计算速度快，适用于大型工程应用；球谐有限元法可用于计算各种形状的几何体，计算应用范围广。

航空发动机；燃烧室；球谐函数法；球谐有限元法；辐射换热；参与性介质

1 引言

航空发动机燃烧室内存在复杂的流动与换热过程，准确模拟燃烧产物与发动机内部结构的换热过程，可为燃烧室内壁等处于严酷工作条件下的结构件的设计提供重要依据。燃烧产物中包含炭黑、二氧化碳气体、水蒸气等，这些介质在高温下具有非常强的辐射能力，研究这类具有吸收、发射、散射能力的介质内的辐射换热过程，通常需要求解辐射传递方程。但由于该方程本身及介质辐射特性的复杂性，理论分析解仅存在于极少数简单或极限问题。为此，国内外学者发展了许多高效数值方法来求解参与性介质内的辐射传递问题，如Monte-Carlo法、区域法(ZM)、离散坐标法(DOM)、球谐函数法(或PN法)、有限元法(FEM)、DRESOR法及上述方法相结合使用的方法等[1～6]。其中球谐函数法采用将辐射强度随空间位置和方向变化进行变量分离的处理方法，利用球谐函数的正交性，将积分-微分形式的辐射传递方程转化为相对简单的偏微分方程组，理论上提供了一种获得任意高阶(精度)的近似求解方法。P1近似相对简单，且与标准能量方程的求解方

法兼容性好，因而在实际问题求解中应用最广。然而随着阶数N和几何体维数的增加，其数学复杂性迅速增加[7]。

球谐函数法的优点在于其不需要空间方向离散，且便于与其它空间离散方法(如FEM、DOM)结合使用，球谐有限元法即结合了球谐函数和有限元空间离散。由于有限元法可处理复杂几何体[8]，因此球谐有限元法也可用来处理各种形状的几何体。球谐有限元法一般用于偶对称中子传输、核测井和医疗成像[9，10]，目前将其应用在参与性介质内辐射传输计算的文献较少。

本文首先对P1有限元法(P1-FEM)进行理论推导，然后通过编制相应计算程序，分别对二维规则和非规则几何体内参与性介质的辐射传输问题进行数值计算。

2 P1有限元法基本理论推导

在发射、吸收、散射介质内，辐射传递方程为：

边界条件为[7]：

将辐射传输方程简化为只与空间坐标有关的二阶微分方程，然后对空间坐标采用有限元离散。令，对直角坐标系下的P1方程即方程(2)进行Galerkin变分，得：

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

对于二维问题的有限元求解域的离散，有限元法可选择不同的离散单元形状[8]，如三角形、四边形、矩形等，但考虑对复杂几何形状的适应性，本文选择了三角形单元，如图1所示。

对于有限元形函数的选择，可以有高次单元和等参单元。本文中选择了较简单的自然坐标系下的等参单元。对于二维情况，形函数可由下面公式给出[8]：

式中：Ni+Nj+Nk=1；Δ为单元面积；aα、bα、cα(α=i，j，m)的定义为[8]

形函数的导数有[8]：

求解域内任意单元e的控制方程变分计算基本公式可写成：

对于等参单元，运用式(6)～式(8)，进而得到以下有限元单元矩阵：

式中：

边界单元控制方程变分计算基本公式可写成：

这里多了一项边界jm的线性积分项计算，对于边界上的线性插值函数：

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

从而得到边界单元系数矩阵：

将式(11)和式(15)带入式(10)进行求解即可。

3 计算模型及结果

例1：规则形状——正方形封腔

如图2所示，一个二维正方形半透明灰体介质被发射率为0.8的不透明边界包围，底面温度1 000 K，其它壁面温度500 K，光学厚度τL=1.0。吸收和散射系数均为0.5 m-1，求正方形封腔内介质的温度分布。

二维正方形封腔的网格分布如图3所示，共有648个三角形网格单元，总节点数为359，边界单元数为68。采用P1-FEM计算得到的封腔内的各向同性介质温度分布如图4所示，并将沿直线y=0.782 5处的温度分布，与采用P1-有限差分(P1-FDM)及文献[6]中DOM-FEM的计算结果进行比较(图5)。结果显示，P1-FEM与P1-FDM、DOM-FEM的计算结果规律一致，说明了P1-FEM程序的正确性；P1-FEM与DOM-FEM(文献[6]已证明DOM-FEM具有高精度，在此假定其计算结果精确)温度场计算的最大误差为6.53%，且主要来自于球谐函数展开部分的计算。

不同网格数量下P1-FEM、P1-FDM与DOM-FEM的计算时间比较见表1(Intel(R)Core(TM)i7-2600：8G内存)。其中P1-FEM和DOM-FEM采用图3所示的三角形网格，网格单元数分别为100、648和2 916；而P1-FDM为计算方便多采用正交四边形网格，网格单元数分别为100、625和2 916。可见，P1-FDM的计算速度最快，但其最大缺点是不能计算非规则几何体的辐射换热，限制了其使用范围；P1-FEM的计算速度快于DOM-FEM，且在网格数目较大时体现得更为明显，说明P1-FEM可以很好地应用于实际工程计算。

例2：非规则形状——非规则四边形

图6所示的非规则四边形，左壁面温度较高为1 000 K，其它壁面为300 K，介质温度也为300 K，壁面为黑体壁面。介质为反照率0.5的吸收、各向同性散射介质。衰减系数为1.0 m-1，求底面无量纲热流密度分布。

计算结果如图7所示，可见P1-FEM在计算无量纲热流分布时与文献[11]的结果有一定差距，这主要是由于球谐函数截断误差所致。球谐函数产生的误差分为假散射和射线效应两类，且都随着PN法中N值的不断增大而减小[12]。但本文并不推荐通过增大N值来减小误差，原因为当N增大时会使数学复杂性增大。如当N=3时，控制方程将由方程(2)一个二阶偏微分方程变为由四个二阶偏微分方程组成的方程组，边界条件方程也由方程(3)一个一阶偏微分方程变为由四个一阶偏微分方程组成的方程组[13]，编程计算过程复杂，计算速度很慢。

4 结束语

本文对用于求解辐射传输的球谐有限元法进行了理论推导，并采用P1-FEM对规则形状和不规则形状几何体内的辐射传输进行了计算。结果表明，球谐有限元法原理简单，计算速度快、精度较好，适用于解决工程问题；与P1方法相比，球谐有限元法可用于计算各种形状的几何体，拓宽了计算应用范围。因此，工程上在计算精度要求不太高的情况下，推荐使用P1-FEM求解复杂形状区域内的辐射换热问题。

[1]谈和平，夏新林，刘林华，等.红外辐射特性与传输的数值计算——计算热辐射学[M].哈尔滨：哈尔滨工业大学出版社，2006：139—190.

[2]张建强，朱谷君.燃烧室中辐射热流分布的蒙特卡罗计算[J].航空动力学报，1999，14(3)：251—254.

[3]毛文懿，林宇震，许全宏，等.运用区域法模型计算燃烧室内一维辐射换热[J].航空动力学报，2010(3)：515—520.

[4]齐宏，阮立明，谭建宇.矩形介质内辐射换热的有限元法[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]谢仲生，邓力.中子输运理论数值计算方法[M]．西安：西北工业大学出版社，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]张昊春，易洪亮，谈和平.球谐函数求解辐射传输方程的假散射和射线效应[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-)，黑龙江大庆人，工程师，博士，主要从事航空发动机传热和燃烧研究。