刘晓川，王存海，黄 勇，朱克勇

(1.北京航空航天大学航空科学与工程学院，北京 100191，2.北京科技大学能源与环境工程学院，北京100083)

辐射传输方程描述了辐射能量在介质中的传递，在许多科学和工程领域具有重要作用，例如大气辐射传输[1]、光学层析成像[2]、天体物理学[3]及核工程[4]等.辐射传输方程是一个高维、复杂的积分微分方程，辐射强度涉及波长、时间、空间和角度等，求其解析解十分困难.学者们提出发展了很多种数值方法来求解辐射传输方程，如蒙特卡洛法[5]，离散坐标法[6]，有限体积法[7]，有限元法[8]等.

近年来，利用格子-Boltzmann方法(LBM)来求解辐射传输方程吸引了许多学者的兴趣.LBM起源于格子气自动机，已经发展成为了一种计算流体力学的有力数值工具[9].并且，LBM已经被拓展到求解许多线性和非线性系统问题，例如声子输运[10]，波传播[11]，反应扩散[12]，对流扩散[13]等.相比于其他的求解辐射传输方程的数值方法，LBM不需要计算大量的光线轨迹，也不需要离散复杂的偏微分方程.LBM具有容易实现，高并行效率等优点.目前，对于利用LBM来求解辐射方程还不完善，发展完善的LBM用于求解辐射传输方程是必要的.

Mishra等[14]假定了可调节的虚拟光速和辐射平衡条件，将LBM推广到分析参与性介质中的辐射问题.Ma等[15]基于辐射流体力学，提出了一维辐射的格子-Boltzmann模型.Zhang等[16]通过采用全隐式后项差分格式处理辐射方程中的瞬态项，将LBM扩展到求解参与性介质中的一维瞬态辐射传输.Mink等[17]在将P1近似应用辐射传输方程的基础上提出了一种三维的格子-Boltzmann模型，然而此模型仅适用于光学厚介质.Yi等[18]通过引入虚拟的扩散项，将辐射传输方程视为一种特殊的对流扩散方程，从而提出了一种二维稳态辐射传输方程的格子-Boltzmann模型.Wang等[19]将瞬态辐射传输方程处理为双曲守恒方程，然后提出了一种求解瞬态辐射和中子输运的格子-Boltzmann模型.

目前，求解辐射方程的多松弛的格子-Boltzmann模型还未见报道.本文提出了一种多松弛格子-Boltzmann模型(multiple-relaxation-time lattice Boltzmann model).基于扩散尺度下的Maxwell迭代，辐射传输方程可以严格地从格子Boltzmann方程推导得出，并且不引入任何限制和近似.本文发展的多松弛格子-Boltzmann模型可以精确地求解参与性介质内的多维瞬态及稳态辐射传输问题.数值结果表明该模型具有二阶精度和收敛速率.并且，相比于单松弛模型，多松弛模型具有更好的稳定性.该模型可以进一步推广到求解参与性介质内的辐射传输问题.

1 辐射传输方程的多松弛格子-Boltzmann模型

1.1 辐射传输方程

考虑吸收、发射和散射介质内的辐射传输方程，其离散坐标形式可以写为[20]

(1)

公式中：cL为介质内的光速；I为辐射强度；r为位置坐标；β=ka+ks为衰减系数；Ωm=μmi+ηmj+ξmk为离散方向，源项S可以表示为

(2)

公式中：N为总的离散方向，m=1，2，…，N，m′=1，2，…，N；wm′为对应方向的权重.

考虑漫发射和反射壁面，边界条件可以写为

(3)

公式中：εw为发射率；ρw为反射率；Iext为外部入射辐射强度.

1.2 多松弛格子-Boltzmann模型

瞬态辐射常常发生于极短的时间内，在瞬态辐射的模拟中，通常引入无量纲时间来避免过小的时间步长.将无量纲时间t*=cLt/LR代入方程(1)中，得到时间无量纲形式的辐射传递方程[21]

(4)

公式中

F(r，Ωm，t*)=LRS(r，Ωm，t*)-LRβI(r，Ωm，t*)

(5)

公式中：LR为介质的参考长度.

本文提出的时间无量纲形式的辐射传输方程的格子Boltzmann方程如下

(6)

公式中：fi(r，t*)为分布函数；M为变换矩阵；S=diag(s0，s1，…，sn)为松弛参数矩阵，平衡函数的表达式为

(7)

辐射强度可以由平衡函数给出，关系如下

(8)

LBM方法中采用DmQn格子模型，对于一维和二维问题，本文分别采用D1Q3和D2Q9模型.对于D1Q3模型，其格子信息为

[c0，c1，c2]=eic=[0 1 -1]c

(9)

(10)

(11)

对于D2Q9模型，其格子信息为

(12)

(13)

(14)

1.3 从格子Boltzmann方程到辐射传输方程

本节基于扩散尺度Δt*=γ(Δx)2下的Maxwell迭代，不引入任何限制和假设，从多松弛格子-Boltzmann模型严格推导得出辐射传输方程.这种扩散尺度是针对模型中的无量纲时间步长和空间步长的尺度.

首先，令f(r，t*)=(f0(r，t*)，f1(r，t*)，…，f8(r，t*))T，ω=(ω0，ω1，…，ω8)T，时间无量纲形式的辐射传递方程(6)可以写成矢量形式

f(r+ciΔt*，t*+Δt*)-f(r，t*)=-M-1SM[f(r，t*)-feq(r，t*)]+Δt*ωF(r，Ωm，t*)

(15)

方程(15)左边应用Taylor展开，

(16)

其中微分算子Ds

(17)

矩阵

Ex=diag(e0，x，e1，x，…，e8，x)Ey=diag(e0，y，e1，y，…，e8，y)

(18)

公式中：p和q均为非负整数.

令m=M·f，meq=M·feq，将Taylor展开形式代入方程(15)并整理得到

(19)

其中

(20)

(21)

m=meq-S-1Lm+γ(Δx)2FS-1Mω

(22)

基于扩散尺度下的Maxwell[22]迭代，从m0=meq开始，方程(19)经过三次迭代得到：

(23)

根据矢量方程(23)的第零项及各算子作用结果，可以得到辐射传递方程

(24)

至此，我们从多松弛格子-Boltzmann模型出发，基于扩散尺度下的Maxwell迭代，严格推导得出了辐射传输方程，并且可以从方程(24)理论上得出该模型具有二阶的精度.一般而言，对于对流扩散问题，计算流体力学等问题的LB模型，其中的松弛系数与宏观方程中的扩散系数，流体黏性系数等有定量关系.需要指出的是，根据从多松弛格子-Boltzmann模型严格推导得出辐射传输方程可知，本文提出的多松弛格子-Boltzmann模型中的松弛参数均是自由的，与其他参数无关.对于一维和二维LB模型，我们取如下的松弛参数矩阵

S=diag(1，sr，1)

(25)

S=diag(1，1，1，sr，1，sr，1，1，1)

(26)

对流扩散方程的多松弛LB模型也采用了同样的处理方法，其中一维模型中的松弛参数s1，二维模型中的松弛参数s3和s5与扩散系数有关，而其他的松弛参数均取1.由于松弛矩阵中的松弛参数有无限种组合方式，因此出于通用性考虑，我们选择了这种处理方法.同时需要指出的是当松弛参数矩阵中的松弛系数相同时，多松弛模型退化到单松弛模型，即松弛矩阵中的松弛参数均为sr.

2 结果及分析

2.1 具有高斯型发射场的一维无限大平板

考虑一充满吸收发射性介质的一维无限大平板内的辐射传递问题，平板内具有一高斯型发射场，该问题由如下方程控制

(27)

考虑如下边界条件

I(0，ξ)=β-1e-b2/α2，ξ>0

(28)

该问题存在解析解形式，其表达式如下

(29)

考虑方向ξ=1.0，a=0.02，b=0.5，采用LBM来模拟衰减系数为β=1，10和50 m-1时介质内辐射强度的分布，取100个格子，无量纲时间步长取Δt*=0.000 1，单松弛模型得到的结果和解析解对比，如图1所示，LBM得到的辐射强度分布和解析解得到的辐射强度分布吻合地很好.

图1 衰减系数为β=1，10和50m-1时LBM得到的辐射强度分布和解析解对比

接下来，我们进一步研究一维多松弛模型的稳定性和精度.为了研究稳定性，我们考虑衰减系数为10 m-1的情况，取100个格子，研究不同松弛参数下所允许的最大时间步长.数值解和解析解的相对误差定义为

(30)

稳定性标准为数值解和解析解的相对误差小于10-2.表1给出了不同松弛参数下所允许的最大时间步长，不同参数的最大时间步长得到是根据我们定义的稳定性标准，然后通过数值实验得到的，可以发现多松弛模型允许的最大时间步长可以随松弛参数调整，尤其当松弛参数小于1时，所允许的时间步长大于单松弛模型，结果表明相比单松弛模型，多松弛模型可以在更大的时间步长内保持稳定，具有更好的稳定性.多松弛模型的碰撞过程发生在矩空间，与多个速度分布函数相关联，相比单松弛模型发生在速度空间的碰撞，多松弛模型本身在稳定性方面展现了很大的优势，数值结果证明了多松弛模型在稳定性上的优势.此外，表2给出了不同格子数下单松弛和多松弛模型的相对误差，可以看出多松弛模型相比单松弛模型具有更高的精度.

表1 衰减系数β=10 m-1，100个格子下，单松弛(BGK)和多松弛(MRT)模型允许的最大时间步长

表2 衰减系数β=10 m-1，不同格子数下，单松弛(BGK)和多松弛(MRT)模型的相对误差

2.2 受高斯型脉冲照射的一维纯散射介质

考虑厚度为L=1 m的一维半透明平板介质内的瞬态辐射传输问题.介质为各向同性散射，壁面和介质温度均为0 K，无发射.介质边界为透明边界，环境为真空.平板介质的衰减系数为1 m-1，右侧边界无照射，左侧边界受到如下法向平行光入射辐射的照射：

(31)

公式中：I0为脉冲辐射强度；H(t)为Heaviside阶跃函数.

图2(a)和图2(b)中的松弛参数分别取sr=0.8，LBM得到的计算结果和文献[24]中Liu和Hsu采用间断有限元方法得到的结果进行对比.由图2可知本文的计算结果和文献结果对比吻合很好，证明了本文提出的多松弛格子-Boltzmann模型可以稳定精确地求解参与性介质内一维瞬态辐射传递问题.本文的LB模型执行基于碰撞和迁徙的显式瞬态演化过程.其它基于离散化偏微分方程的数值方法，如FVM、DOM、无网格法和FEM等，一般都需要在每个时间步长内进行全局迭代收敛，在效率上并不优于LB模型.因此，本文提出的LB模型非常适合于求解瞬态辐射传输问题.在处理稳态辐射传输问题时，与基于离散化偏微分方程的数值方法相比，LB模型的缺点是计算效率较低.这是因为LBM通过依赖时间的演化过程来求解稳态问题.对于稳态问题，这些基于离散化偏微分方程的数值方法只需对稳态辐射传输方程进行全局迭代收敛即可获得收敛解，因而对稳态辐射传输问题具有较高的计算效率.

图2 高斯型脉冲照射下平板界面处时域反射率和透射率信号

2.3 二维方腔内辐射传递问题

本算例考虑二维方腔内的辐射传递问题，方腔的边长为L=1 m，衰减系数为1 m-1，方腔内介质为热介质，所有壁面保持0 K，四个壁面均为黑壁面.方腔内充满吸收散射性介质，散射为各向异性散射，散射相函数F2的不对称因子为0.669 72(见文献[25]).将计算域划分为60×60的格子，无量纲时间步长取Δt*=0.000 5，空间离散采用S8方案.松弛参数取0.8时的多松弛模型得到的方腔底部无量纲热流，如图3所示.同时文献[25]中采用离散坐标法得到的计算结果也显示在图3中作为对比，结果显示LBM得到的计算结果和离散坐标法得到的计算结果吻合得很好，进一步验证了本文发展得多松弛格子-Boltzmann模型可以准确地求解参与性介质内多维辐射传输问题.

图3 不同散射反照率下方腔底部的无量纲辐射热流

3 结 论

对于辐射方程，本文提出建立了一种多松弛的格子-Boltzmann模型.对多松弛格子-Boltzmann模型进行Maxwell迭代，可以严格地得出辐射方程.相比于已有的格子-Boltzmann模型，本文提出的模型没有任何近似和限制.

数值实验表明本文发展的多松弛格子-Boltzmann模型可以稳定精确地求解参与性介质内一维、二维瞬态和稳态辐射传输问题.同时，理论结果表明该模型具有二阶精度.本文构建的多松弛的格子-Boltzmann模型可以退化到单松弛格子-Boltzmann模型，且多松弛模型相比于单松弛模型具有更好的稳定性.采用三维的格子模型，多松弛模型也可以用来求解三维辐射传输问题，与此同时，本文提出的多松弛的格子-Boltzmann模型可以进一步推广到求解参与性介质内复杂的辐射传输问题.