基于达菲变换消除“奇点”影响的辐射直接交换面积改进算法

Improved Algorithm for Radiation Direct Exchange Area Based on Duffy Transformation to Eliminate Effects of “Singularities”

摘要

应用积分法求解段法中辐射直接交换面积（DEA）时由于“奇点”的存在计算误差大，为提高“奇点”处计算精度往往需要付出大量的计算时间且收效甚微.为此，基于达菲（Duffy）变换方法，推导了“奇点”处辐射DEA达菲变换积分式，用以消除“奇点”对辐射DEA计算精度的影响.分别采用改进的方法与直接高斯积分(DGI)法求解了二维方腔内的辐射DEA.结果表明，应用达菲变换后，在同样的计算时间内可将辐射直接交换面积气体段及表面段的完整性验证误差分别由3.73%，6.70%降低到(4.88×10^-4)%，(2.98×10^-5)%；而在满足相同计算精度情况下，二者计算时间最大相差13 693倍.可见，所提改进算法在计算速度和计算精度上都有显著的优势.

Abstract

When the integral method is used to solve the direct exchange area （DEA） of radiation in the zone method， the relative error is large due to the existence of “singularities”， and it usually requires a lot of computing time and has little effect to improve the calculation accuracy at “singularities”. Therefore， based on Duffy transformation method， the Duffy transformation integral formula is derived to eliminate the influence of the singularity on the calculation of DEA. The improved method and the direct Gaussian integral （DGI） method are used to solve the DEA of radiation in a two-dimensional square cavity. The results show that the completeness verification errors of the gas-zone and the surface-zone can be decreased from 3.73% and 6.70% to （4.88×10^-4）% and （2.98×10^-5）% respectively in the same calculation time by applying Duffy transform. In the case of the same calculation accuracy， the maximum difference in calculation time between them is 13 693 times. Therefore， the improved algorithm has significant advantages in computation speed and accuracy.

Graphical abstract

关键词

达菲变换 / “奇点”问题 / 直接交换面积 / 段法 / 辐射换热

Key words

Duffy transformation / singularity problem / direct exchange area（DEA） / zone method / radiative heat transfer

热传导、对流和辐射是热量传递的3种基本方式，随着温度的升高，辐射换热逐渐占据主导地位，准确计算辐射换热是一个非常值得研究的课题.目前计算辐射换热的方法有：段法^［1］，离散坐标法^［2-3］，蒙特卡罗法^［4-5］等.其中，段法能够有效描述段与段之间的辐射换热过程，其核心在于采用辐射全交换面积来表征不同区段间的辐射热传递.目前求解辐射全交换面积主要采用两种方法：一是基于蒙特卡罗法的直接计算^［6］；二是通过高斯数值积分法先求解辐射直接交换面积，再进一步得到全交换面积.尽管积分法在求解直接交换面积时需要消耗较多计算资源，但由于其针对特定几何构型只需单次计算即可获得辐射全交换面积，这一显著优势使得段法在工业热工装备的辐射传热分析中得到广泛应用.

在采用积分方法求解辐射直接交换面积时，在邻近的表面段与表面段之间、邻近的气体段与表面段之间及气体段自身辐射的计算中存在“奇点”情况，导致辐射直接交换面积计算误差偏大，一般需要对此进行特殊的处理.为了解决这个问题，国内外学者提出了多种改进方法.Tian等^［7］采用有限容积法与中点结合的方法计算辐射直接交换面积，在自身对自身的辐射段以及相邻的段采用有限容积法以减少辐射直接交换面积的影响.Ebrahimi等^［8］采用简化段法，将段更加地细分，使段对自身的辐射变成了多个气体段非自身的辐射以提高“奇点”处的精度.文献［9-10］使用降重积分法降低了积分的重数，减少了计算量并且略微改善了“奇点”处的计算精度，之后使用改进蒙特卡罗法大幅降低了计算时间，尤其是气体吸收系数较小时的计算时间^［11］.

已有研究^［12］针对“奇点”问题采取了改进方法，包括按经验取近似、局部网格加密以及舍弃存在“奇点”的值，但这些方法难以有效消除“奇点”的影响.本文采用达菲变换这一数学方法，推导出可以避免“奇点”的积分计算式，计算了辐射直接交换面积，降低了“奇点”处的积分重数.与直接高斯积分法相比，该方法在计算精度和计算速度上均具有显著优势，为段法在在线控制中的应用奠定了基础.

1 辐射直接交换面积与“奇点”问题

在段法中，将温度场划分为表面段与气体段，采用A_{s_i{s}_j}代表表面段s_i与表面段s_j之间的辐射直接交换面积，采用A_{g_i{s}_j}代表气体段g_i与表面段s_j之间的辐射直接交换面积，采用A_{g_i{g}_j}代表气体段g_i与气体段g_j之间辐射直接交换面积，在求解出辐射直接交换面积之后，再求解辐射全交换面积，即可用于计算辐射场中段与段之间的辐射换热量.在过去的认知中，求解辐射直接交换面积是耗费计算时间最大的部分，“奇点”的存在会显著地降低计算精度，为了获取可以接受的计算精度而不得不加大计算量.

1.1 辐射直接交换面积的计算

辐射换热的计算公式为

A_{s_i{s}_j}=\underset{l_{\mathrm{1}}}{\int }\underset{l_{\mathrm{2}}}{\int }\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{\mathrm{2}r}\mathrm{d}{l}_{\mathrm{2}}\mathrm{d}{l}_{\mathrm{1}}，(1)

A_{g_i{s}_j}=\underset{F}{\int }\underset{l}{\int }{\kappa }_{\mathrm{a}}\mathrm{c}\mathrm{o}\mathrm{s}\theta \frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{\mathrm{2}r}\mathrm{d}l\mathrm{d}F，(2)

A_{g_i{g}_j}=\underset{F_{\mathrm{1}}}{\int }\underset{F_{\mathrm{2}}}{\int }{\kappa }_{\mathrm{a}}^{\mathrm{2}}\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{\mathrm{2}r}\mathrm{d}{F}_{\mathrm{2}}\mathrm{d}{F}_{\mathrm{1}}.(3)

式中：r为微元之间的距离；{\kappa }_{\mathrm{a}}为气体吸收系数；l是二维模型中表面段的长度；l₁和l₂分别是表面段s_i 和s_j 的长度；F是二维模型气体段的面积；F₁和F₂分别是气体段g_i 和g_j 的面积.

在正交网格中，按照图1所示几何位置，其具体的积分公式如下：

互相平行的表面段之间的直接交换面积为

A_{s_i{s}_j}={\int }_{x_{\mathrm{1},a}}^{x_{\mathrm{1},b}}{\int }_{y_{\mathrm{2},a}}^{y_{\mathrm{2},b}}\frac{Y^{\mathrm{2}}}{\mathrm{2}{r}^{\mathrm{3}}}{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}\mathrm{d}{x}_{\mathrm{2}}\mathrm{d}{x}_{\mathrm{1}}.(4)

式中，r=\sqrt[]{(x_{\mathrm{2}}-x_{\mathrm{1}}{)}^{\mathrm{2}}+Y^{\mathrm{2}}}.

互相垂直的表面段之间的直接交换面积为

A_{s_i{s}_v}={\int }_{x_a}^{x_b}{\int }_{y_a}^{y_b}\frac{\left|y(x-X)\right|}{\mathrm{2}{r}^{\mathrm{3}}}{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}\mathrm{d}y\mathrm{d}x.(5)

式中，r=\sqrt[]{{(x-X)}^{\mathrm{2}}+y^{\mathrm{2}}}.

气体段与表面段之间的直接交换面积为

A_{g_i{s}_j}={\kappa }_{\mathrm{a}}({\int }_{y_a}^{y_b}{\int }_{x_{\mathrm{1},a}}^{x_{\mathrm{1},b}}{\int }_{x_{\mathrm{2},a}}^{x_{\mathrm{2},b}}{\kappa }_{\mathrm{a}}y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}{x}_{\mathrm{2}}\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}y.(6)

式中，r=\sqrt[]{(x_{\mathrm{2}}-x_{\mathrm{1}}{)}^{\mathrm{2}}+y^{\mathrm{2}}}.

气体段之间的直接交换面积为

A_{g_i{g}_j}={\int }_{x_{\mathrm{1},a}}^{x_{\mathrm{1},b}}{\int }_{y_{\mathrm{1},a}}^{y_{\mathrm{1},b}}{\int }_{x_{\mathrm{2},a}}^{x_{\mathrm{2},b}}{\int }_{y_{\mathrm{2},a}}^{y_{\mathrm{2},b}}{\kappa }_{\mathrm{a}}^{\mathrm{2}}\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{\mathrm{2}r}\mathrm{d}{y}_{\mathrm{2}}\mathrm{d}{x}_{\mathrm{2}}\mathrm{d}{y}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{1}}.(7)

式中，r=\sqrt[]{(x_{\mathrm{2}}-x_{\mathrm{1}}{)}^{\mathrm{2}}+(y_{\mathrm{2}}-y_{\mathrm{1}}{)}^{\mathrm{2}}}.

上述计算是二、三、四重积分，为获得更高精度，对其进行复化分段时，计算量随着积分重数呈指数增长，这就是段法中直接交换面积计算量大的原因.

1.2 “奇点”问题

如式 （5）~式 （7）所示，对于气体段自身辐射，气体段与相邻表面段以及相邻的垂直表面段，此时分母中r=\mathrm{0}，那么这个值无法计算，已有研究有不同的处理方式，包括按照经验对r赋值或在物理模型上对段做一微小位移以及舍弃奇异点（设为0），但这些处理方式都会引入一定的误差；即使不是在r=\mathrm{0}处，还有一些r的值非常小并且作为分母进行计算，在计算机的计算中，会引起精度损失，最好的方法还是采用数学方法对其变换以彻底地解决这个问题.如图1a所示，可用式（4）计算平行的段之间的辐射直接交换面积，因为平行段之间距离的存在，可以避免r=\mathrm{0}的发生，在此处无需特殊处理.

2 “达菲”变换与辐射直接交换面积计算

已有研究^［13］对“奇点”问题即奇异积分的改进方法是作变量替换，这个方法被称为达菲变换，Taylor^［14］扩大了应用范围，Feischl等^［15］尝试将其应用到段法中.对于本文研究的正交网格，仅仅需要对以下3个辐射直接交换面积进行变换即可，而且达菲变换也仅仅适合定义域在（0，C）的积分区域的一类问题，其中C代表确定的常数.

2.1 “达菲”变换原理

文献［13］提到的达菲变换方法，是针对一系列如下所示问题：

I={\int }_{\mathrm{0}}^C\frac{f(x)}{x^{\alpha }}\mathrm{d}x.(8)

式中，α是常数，取t=x^{\alpha }即可解决数值积分时的奇异性问题，这就是最基本的达菲变换，但是辐射直接交换面积的计算过程中的分母是r=\sqrt[]{(x_{\mathrm{1}}-x_{\mathrm{2}}{)}^{\mathrm{2}}+（y_{\mathrm{1}}-y_{\mathrm{2}}{)}^{\mathrm{2}}}，在达菲变换之前要作一定的处理.

第一个特殊处理是将r中的x_{\mathrm{1}}-x_{\mathrm{2}}替换为x=x_{\mathrm{1}}-x_{\mathrm{2}}，y_{\mathrm{1}}-y_{\mathrm{2}}替换为y=y_{\mathrm{1}}-y_{\mathrm{2}}，这样就将r=\sqrt[]{(x_{\mathrm{1}}-x_{\mathrm{2}}{)}^{\mathrm{2}}+（y_{\mathrm{1}}-y_{\mathrm{2}}{)}^{\mathrm{2}}}转变为r=\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}}.

第二个特殊处理是在定义域上进行划分，使划分后的定义域适合使用达菲变换，其定义域为

I={\int }_{\mathrm{0}}^C{\int }_{\mathrm{0}}^x\frac{f(x,y)}{r}\mathrm{d}y\mathrm{d}x.(9)

式中，r=\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}}.

在本文中存在“奇点”的积分式中，经过这两步变换之后，就可以将原有积分式降重，这个简化是因为积分式的特殊性而出现的，不是所有的达菲变换都可以降低积分重数，在这里可以有效地降低计算量.

经过这两个特殊处理后，即可使用达菲变换以消除奇异性.接下来针对辐射直接交换面积求解中存在“奇点”的3种情况进行推导.

2.2 互相垂直的表面段之间的辐射直接交换面积的积分式推导

如图1b所示的垂直表面段之间辐射直接交换面积的计算公式为式（5），当两个互相垂直的表面段相邻时就会存在“奇点”问题，此时辐射交换面积计算的积分区域为x\in (\mathrm{0},A)，y\in (\mathrm{0},A).则可得

A_{s_{\mathrm{1}}{s}_v}={\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A\frac{\left|xy\right|}{\mathrm{2}{r}^{\mathrm{3}}}{\mathrm{e}}^{-{\kappa }_{a}r}\mathrm{d}y\mathrm{d}x.(10)

式中，r=\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}}.

首先对其定义域进行分割，划分为适合达菲变换的形式，即

A_{s_{\mathrm{1}}{s}_v}={\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^x\frac{xy}{r^{\mathrm{3}}}{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}\mathrm{d}y\mathrm{d}x.(11)

式中，r=\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}}.

然后对其采用达菲变换y=ux，可得

A_{s_{\mathrm{1}}{s}_v}={\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{\mathrm{1}}\frac{x^{\mathrm{3}}u}{r^{\mathrm{3}}}{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}\mathrm{d}u\mathrm{d}x.(12)

式中，r=x\sqrt[]{\mathrm{1}+u^{\mathrm{2}}}.

为了统一r的计算形式，作变量替换可得

A_{s_{\mathrm{1}}{s}_v}={\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{\mathrm{1}}\frac{x{y}^{\mathrm{3}}}{r^{\mathrm{3}}}{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}\mathrm{d}x\mathrm{d}y.(13)

式中，r=y\sqrt[]{\mathrm{1}+x^{\mathrm{2}}}.

2.3 相邻的气体段与表面段之间的辐射直接交换面积的积分式推导

如图1c所示，相邻的气体段与表面段之间的辐射直接交换面积计算公式为式（6），当气体段与表面段相邻时就会存在“奇点”问题，此时辐射交换面积计算的积分区域为x_{\mathrm{1}}\in (\mathrm{0},A)，x_{\mathrm{2}}\in (\mathrm{0},A)，y\in (\mathrm{0},A).则可得

A_{g_{\mathrm{1}}{s}_{\mathrm{1}}}={\kappa }_{\mathrm{a}}({\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\kappa }_{\mathrm{a}}y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}{x}_{\mathrm{2}}\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}y.(14)

式中，r=\sqrt[]{(x_{\mathrm{2}}-x_{\mathrm{1}}{)}^{\mathrm{2}}+y^{\mathrm{2}}}.

此时需要对方程进行处理，以满足达菲变换的要求，令x=x_{\mathrm{1}},u=x_{\mathrm{2}}-x_{\mathrm{1}}\to \mathrm{d}u=\mathrm{d}{x}_{\mathrm{2}}

A_{g_{\mathrm{1}}{s}_{\mathrm{1}}}=\frac{{\kappa }_{\mathrm{a}}}{\mathrm{2}}{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\int }_{-x}^{A-x}y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}u\mathrm{d}x\mathrm{d}y.(15)

式中，r=\sqrt[]{u^{\mathrm{2}}+y^{\mathrm{2}}}.

对积分区间进行适当划分后可得

\begin{array}{l}{A}_{g_{\mathrm{1}}{s}_{\mathrm{1}}}=\frac{{\kappa }_{\mathrm{a}}}{\mathrm{2}}{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{A-u}y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}x\mathrm{d}u\mathrm{d}y+\\ \frac{{\kappa }_{\mathrm{a}}}{\mathrm{2}}{\int }_{\mathrm{0}}^A{\int }_{-A}^{\mathrm{0}}{\int }_{-u}^{A}y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}x\mathrm{d}u\mathrm{d}y.\end{array}(16)

对变量x进行积分后可得

\begin{array}{l}{A}_{g_{\mathrm{1}}{s}_{\mathrm{1}}}=\frac{{\kappa }_{\mathrm{a}}}{\mathrm{2}}{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A(A-u)y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}u\mathrm{d}y+\\ \frac{{\kappa }_{\mathrm{a}}}{\mathrm{2}}{\int }_{\mathrm{0}}^A{\int }_{-A}^{\mathrm{0}}(A+u)y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}u\mathrm{d}y.\end{array}(17)

因为此时y>\mathrm{0}，继续化简可得

\begin{array}{l}{A}_{g_{\mathrm{1}}{s}_{\mathrm{1}}}={\kappa }_{\mathrm{a}}A{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{A}y\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}u\mathrm{d}y-\\ {\kappa }_{\mathrm{a}}{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{A}uy\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}u\mathrm{d}y.\end{array}(18)

式中，r=\sqrt[]{u^{\mathrm{2}}+y^{\mathrm{2}}}.

对式（18）中的两个部分按照式（10）~式（13）的过程作达菲变换后结果为

A_{g_{\mathrm{1}}{s}_{\mathrm{1}}}={\kappa }_{\mathrm{a}}({\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{\mathrm{1}}(A+Ax-\mathrm{2}xy)y^{\mathrm{2}}\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r^{\mathrm{2}}}\mathrm{d}x\mathrm{d}y).(19)

式中，r=y\sqrt[]{\mathrm{1}+x^{\mathrm{2}}}.

2.4 气体段对自身的辐射直接交换面积的积分式推导

如图1d所示的两个气体段之间的辐射直接交换面积计算公式为式（7），对于一个段既是发射段又是吸收段的情况，即气体段自身辐射的计算，此时也存在大量“奇点”，积分区域为x_{\mathrm{1}}\in (\mathrm{0},A)，x_{\mathrm{2}}\in (\mathrm{0},A)，y_{\mathrm{1}}\in (\mathrm{0},A)，y_{\mathrm{2}}\in (\mathrm{0},A).可得

A_{g_{\mathrm{1}}{g}_{\mathrm{1}}}={\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\kappa }_{\mathrm{a}}^{\mathrm{2}}\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{\mathrm{2}r}\mathrm{d}{y}_{\mathrm{2}}\mathrm{d}{x}_{\mathrm{2}}\mathrm{d}{y}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{1}}.(20)

式中，r=\sqrt[]{(x_{\mathrm{2}}-x_{\mathrm{1}}{)}^{\mathrm{2}}+(y_{\mathrm{2}}-y_{\mathrm{1}}{)}^{\mathrm{2}}}.

首先作如下积分变换u=X-x,v=Y-y\to \mathrm{d}u=\mathrm{d}X,\mathrm{d}v=\mathrm{d}Y，则可得

A_{g_{\mathrm{1}}{g}_{\mathrm{1}}}={\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^A{\int }_{-x}^{A-x}{\int }_{-y}^{A-y}\frac{{\kappa }_{\mathrm{a}}^{\mathrm{2}}}{\mathrm{2}r}{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}\mathrm{d}v\mathrm{d}u\mathrm{d}y\mathrm{d}x.(21)

式中，r=\sqrt[]{{(X-x)}^{\mathrm{2}}+{(Y-y)}^{\mathrm{2}}}.

按照2.3节类似的推导过程可得

A_{g_{\mathrm{1}}{g}_{\mathrm{1}}}=\mathrm{4}{\kappa }_{\mathrm{a}}^{\mathrm{2}}({\int }_{\mathrm{0}}^A{\int }_{\mathrm{0}}^{\mathrm{1}}y(A-y)(A-xy)\frac{{\mathrm{e}}^{-{\kappa }_{\mathrm{a}}r}}{r}\mathrm{d}x\mathrm{d}y.(22)

式中，r=y\sqrt[]{\mathrm{1}+x^{\mathrm{2}}}.

图2是使用直接高斯积分的传统算法和基于“达菲变换”的改进算法计算辐射直接交换面积的流程图.

3 结果与讨论

3.1 计算精度比较方法

由于多重数值积分的误差难以分析，已有研究利用物理意义对计算精度进行分析^［13］，辐射直接交换面积一定满足一些条件：

\left.\begin{array}{c}{A}_{s_i{s}_j}=A_{s_j{s}_i},\\ A_{g_i{s}_j}=A_{s_j{g}_i},\\ A_{g_i{g}_j}=A_{g_j{g}_i}.\end{array}\right\}(23)

式（23）说明了辐射直接交换面积具有“交换性”.

\left.\begin{array}{c}{\displaystyle \sum _{k=\mathrm{1}}^M}{A}_{s_i{s}_k}+{\displaystyle \sum _{k=\mathrm{1}}^N}{A}_{s_i{g}_k}=l,\\ {\displaystyle \sum _{k=\mathrm{1}}^N}{A}_{g_i{g}_k}+{\displaystyle \sum _{k=\mathrm{1}}^M}{A}_{g_i{s}_k}=\mathrm{\pi }{\kappa }_{\mathrm{a}}F.\end{array}\right\}(24)

式中，M，N分别为表面段和气体段的数量.式（24）说明辐射直接交换面积具有“完整性”.

蒙特卡罗法将能量赋予到能束中并对能束进行追踪，完整性验证结果较好，而其交换性则因为随机性的影响较差^［16］，所以使用交换性或者使用多次计算的方差对蒙特卡罗法的精度进行评价更加合理；对于积分法求解辐射直接交换面积来说，其一定可以满足交换性，可以借此对几何位置相同的直接交换面积进行简化计算，但较难满足完整性验证，所以使用完整性对段法的计算精度进行评价更合理.即辐射直接交换面积的误差定义为

\left.\begin{array}{c}{e}_{\mathrm{s}}=\frac{\left|l-({\displaystyle \sum _{k=\mathrm{1}}^M}{A}_{s_i{s}_k}+{\displaystyle \sum _{k=\mathrm{1}}^N}{A}_{s_i{g}_k})\right|}{l},\\ e_{\mathrm{g}}=\frac{\left|\mathrm{\pi }{\kappa }_{\mathrm{a}}F-({\displaystyle \sum _{k=\mathrm{1}}^M}{A}_{g_i{s}_k}+{\displaystyle \sum _{k=\mathrm{1}}^N}{A}_{g_i{g}_k})\right|}{\mathrm{\pi }{\kappa }_{\mathrm{a}}F}.\end{array}\right\}(25)

式中，e_s和e_g分别为表面段和气体段的完整性误差.

3.2 计算精度分析

选择一些算例对达菲变换的计算精度进行分析.在接下来的计算中，对于气体段假设其为各向同性的灰色气体，{\kappa }_{\mathrm{a}}=\mathrm{0.5}{\mathrm{m}}^{-\mathrm{1}}.

3.2.1 试验算例1

为了分析达菲变换的效果，需要设定一个以达菲变换后的积分式为主的算例，选择一个填充了吸收性介质的1 m×1 m的二维方腔，这样仅仅有一个相互平行的表面段之间的直接交换面积A_{s_{\mathrm{1}}{s}_{\mathrm{2}}}使用直接高斯积分计算，其他的直接交换面积都使用达菲变换后的积分式进行计算.

如图3a所示，按照这种段的划分方法，气体段为1 m×1 m，表面段长1 m.

由图4可得，达菲变换的计算结果与直接高斯积分的结果非常接近.应用完整性验证对计算精度进行评价，由表1的完整性验证结果可得，其气体段与表面段的相对精度分别由直接高斯积分的0.58%与0.36%提高到了达菲变换后的（2.5×10^-7）%与（2.1×10^-6）%，由此可得，单独使用达菲变换改善了存在“奇点”的辐射直接交换面积的完整性.

3.2.2 试验算例2

如图3b所示，选择一个5 m×5 m二维方腔，但是划分了更多的气体段与表面段，有25个气体段以及20个表面段.按照这种段的划分方法，气体段为1 m×1 m，表面段长1 m.将达菲变换与直接高斯积分结合使用以计算全部的辐射直接交换面积.

由图5可以看出，达菲变换计算出的A_{g_{\mathrm{1}}{g}_{\mathrm{1}}}，A_{g_{\mathrm{1}}{s}_{\mathrm{1}}}以及A_{s_{\mathrm{1}}{s}_{\mathrm{11}}}与直接高斯积分计算的结果非常吻合，验证了达菲变换的有效性.表2展示了使用达菲变换计算前后完整性验证，可以看到，其可以有效地提高计算的精度.将气体段与表面段的辐射直接交换面积计算的相对误差由0.093 5%，0.169%分别降低到了（1.73×10^-5）%，（3.92×10^-6）%.由此表明，达菲变换与直接高斯积分结合使用以计算辐射直接交换面积是非常有效的.

3.3 计算速度分析

表3与表4给出了达菲变换与使用不同积分步长的直接高斯积分计算上述两个算例时所需要的计算时间以及相对误差.在算例1中，同样花费0.005 s时，达菲变换的气体段与表面段的相对误差分别只有（4.81×10^-6）%，0%；而直接高斯积分的气体段与表面段的相对误差为11.6%，7.23%；在算例2中，相同计算时间的情况下，达菲变换的气体段与表面段的相对误差分别为（4.88×10^-4）%，（2.98×10^-5）%，而直接高斯积分的相对误差则分别为3.73%，6.70%.由此可见，在同样的计算时间内，达菲变换可以获取更高的精度.

算例1中，在气体段相对误差都达到0.364%时，直接高斯积分计算用时68.465 s，而同样能够达到此计算精度的达菲变换计算用时仅0.005 s，计算时间相差13 693倍；算例2中，在气体段相对误差都达到0.117%时，直接高斯积分计算用时957.069 s，而达菲变换计算用时仅0.071 s，计算时间相差13 479.8倍.由此可见，满足相对误差相同的情况下，二者计算时间最大相差13 693倍.

应用达菲变换求解辐射直接交换面积时，在相同计算时间时计算精度更高，在达到同样的相对误差时，计算时间更短，因此本文提出的基于达菲变换的辐射直接交换面积算法在计算速度上有显著优势.

4 结 论

1） 在二维模型辐射换热的计算中使用达菲变换，将相邻且垂直的表面段之间的二重积分、相邻的表面段与气体段之间的三重积分、气体段对气体段自身的四重积分都降为二重积分.

2） 使用达菲变换可以消除段法中辐射换热直接交换面积计算过程中“奇点”问题的影响，从而显著提高计算精度.在算例1中，使用达菲变换后气体段与表面段的完整性评价的相对精度由11.6%，7.23%分别提高到（4.81×10^-6）%，0%；在算例2中，使用达菲变换后气体段与表面段的完整性评价的相对精度由3.73%，6.70%分别提高到（4.88×10^-4）%，（2.98×10^-5）%.

3） 达菲变换可以显著地提高辐射直接交换面积的计算速度.算例1中气体段达到同样0.364%的相对误差时，达菲变换的计算时间为直接高斯积分的1/13 693；算例2中同样达到0.117%的相对误差时，达菲变换的计算时间为直接高斯积分的1/13 479.8.这为段法处理更复杂的模型以及更广泛的应用提供了基础.

参考文献

[1] Larsen M E， Howell J R. The exchange factor method： an alternative basis for zonal analysis of radiating enclosures［J］. Journal of Heat Transfer， 1985， 107（4）： 936-942.

[2] Fiveland W A. Discrete ordinate methods for radiative heat transfer in isotropically and anisotropically scattering media［J］. Journal of Heat Transfer， 1987， 109（3）： 809-812.

[3] Kim T K， Lee H S. Radiative transfer in two-dimensional anisotropic scattering media with collimated incidence［J］. Journal of Quantitative Spectroscopy and Radiative Transfer， 1989， 42（3）： 225-238.

[4] He F， Shi J S， Zhou L， et al. Monte Carlo calculation of view factors between some complex surfaces： rectangular plane and parallel cylinder， rectangular plane and torus， especially cold-rolled strip and W-shaped radiant tube in continuous annealing furnace［J］. International Journal of Thermal Sciences， 2018， 134： 465-474.

[5] Zhou W G， Qiu T. Zone modeling of radiative heat transfer in industrial furnaces using adjusted Monte-Carlo integral method for direct exchange area calculation［J］. Applied Thermal Engineering， 2015， 81： 161-167.

[6] Li D Y， Li G J， Hong D L， et al. Improved Monte Carlo method for radiative heat transfer in semitransparent media with BRDF surface［J］. International Journal of Thermal Sciences， 2023， 187： 108152.

[7] Tian W X， Chiu W K S. Hybrid method to calculate direct exchange areas using the finite volume method and midpoint integration［J］. Journal of Heat Transfer， 2005， 127（8）： 911-917.

[8] Ebrahimi H， Zamaniyan A， Mohammadzadeh J S S， et al. Zonal modeling of radiative heat transfer in industrial furnaces using simplified model for exchange area calculation［J］. Applied Mathematical Modelling， 2013， 37（16/17）： 8004-8015.

[9] 李国军，陈海耿，伊智.积分降重法求解辐射直接交换面积［J］.东北大学学报（自然科学版）， 2010， 31（1）： 80-83.

[10] Li Guo-jun， Chen hai-geng， Yi Zhi. Calculation of radiative direct exchange areas using reduced integration scheme［J］. Journal of Northeastern University （Natural Science）， 2010， 31（1）： 80-83.

[11] Li G J， Li B W， Sun Y S. Extension of the reduced integration scheme to calculate the direct exchange areas in 3D rectangular enclosures with nonscattering media［J］. Mathematical Problems in Engineering， 2015， 2015： 703823.

[12] 李国军，钟佳琪，李定勇，等.基于改进蒙特卡洛法求解辐射传热问题的计算精度评价［J］.湖南大学学报（自然科学版）， 2022， 49（2）： 55-62.

[13] Li Guo-jun， Zhong Jia-qi， Li Ding-yong， et al. Evaluation on computational accuracy for improved Monte Carlo method of radiative heat transfer problem［J］. Journal of Hunan University （Natural Sciences）， 2022， 49（2）： 55-62.

[14] Xie G Z， Zhong Y D， Li H， et al. Near singularity cancellation in weakly singular integrals of three-dimensional boundary element method［J］. Engineering Analysis with Boundary Elements， 2020， 118： 54-59.

[15] Duffy M G. Quadrature over a pyramid or cube of integrands with a singularity at a vertex［J］. SIAM Journal on Numerical Analysis， 1982， 19（6）： 1260-1262.

[16] Taylor D J. Accurate and efficient numerical integration of weakly singular integrals in Galerkin EFIE solutions［J］. IEEE Transactions on Antennas and Propagation， 2003， 51（7）： 1630-1637.

[17] Feischl M， Führer T， Niederer M， et al. Efficient numerical computation of direct exchange areas in thermal radiation analysis［J］. Numerical Heat Transfer， Part B： Fundamentals， 2016， 69（6）： 511–533.

[18] Li D Y， Li G J， Wei L Y， et al. Reduced integration coupled with Monte Carlo ratios method for zone modeling of radiative heat transfer in reheating furnaces［J］. International Journal of Thermal Sciences， 2024， 195： 108640.

基金资助

国家自然科学基金资助项目(52106079)

中央高校基本科研业务费专项资金资助项目(N2325021)