离散统一气体动理学格式稳态辐射输运应用

张华波, 周瑞睿, 李思达, 孙亚松,3

(1. 西北工业大学动力与能源学院, 陕西西安 710129; 2. 上海理工大学能源与动力工程学院, 上海 200093; 3. 西北工业大学太仓长三角研究院, 江苏太仓 215400)

引 言

在天体物理、 惯性约束聚变、 冶金行业、 太阳能热利用和发动机高温区域分析等诸多领域的研究中热辐射具有重要地位[1-3]. 以航空发动机中的燃烧室为例, 其核心区温度逾2 400 K[4], 热辐射是主要的能量传递方式之一. 刻画热辐射传递过程的辐射传递方程是一类较为复杂的微分积分方程, 难以获得通用的解析解[1,5], 如何准确、 高效地求解辐射传递方程是近几十年来热辐射研究的重要内容.

由于热辐射高温高压的实验条件通常较为苛刻, 越来越多的国内外学者采用数值手段对热辐射问题进行了研究. 求解热辐射的方法可按照计算和实现形式的不同分为两大类: 第一类是基于Lag-range坐标系的随机模拟方法, 如Monte Carlo法(Monte Carlo method, MCM)[6-7]、 离散传递法(discrete transfer method, DTM)[8-9]、 区域法[10-11]等. 该类方法计算量相对较大, 对于结构相对简单的模型都需要很大的计算量, 特别是对于光学厚区域的情况, 将更加难以计算. 另外, 采用Monte Carlo法在光学厚区域还会面临很大的统计噪声问题. 第二类是基于Euler坐标系的方程离散求解方法, 如离散坐标法(discrete ordinates method, DOM)[12-13]、 谱元法(spectral element method, SEM)[12-14]、 有限元法(finite element method, FEM)[15-16]、 无网格法[17-19]、 有限体积法(finite volume method, FVM)[20-22]、 球形谐波法(又称PN法)[23-24]等. 这类方法直接求解辐射传递方程或其近似方程, 计算效率高. 但由于这类方法需要对空间立体角进行离散处理, 从而引入射线效应影响计算结果的精度, 其中PN法没有射线效应; 同时对空间位置的离散, 也会引起假散射， 对计算结果产生影响. 另外出于稳定性的考虑, 采用显式格式需要时间步长小于光子平均自由程的平方、 空间尺寸小于光子平均自由程. 采用隐式格式虽然没有时间步长的严格限制, 但是光学厚区域网格数量的大幅度增加仍然将导致计算量急剧上升.

上述求解方法在光学厚区域的热辐射求解面临一定困难, 可以通过一些近似处理的方法来解决[25], 例如基于扩散方程的P1近似和扩散近似. 但是这类近似方法只能在光学厚度极大时才能够获得较为准确的结果. 对于光学厚度变化的情况, 可能会出现较大的数值误差. 对于参与性介质辐射, 辐射介质的光谱性质、 温度、 压力、 射线成分等存在较大差异导致光学厚度发生很大变化, 从而使得热辐射呈现多尺度特征. 近年来, 针对热辐射的多尺度问题, 主要发展了3类统一考虑不同光学尺度下的热辐射计算方法.

第1类方法是在区域分解法的框架下, 实现不同尺度区域的单独建模和耦合求解[26], 即在光学薄区域采用辐射输运模型进行计算, 在光学厚区域采用扩散模型进行计算, 并在两个区域之间引入一个缓冲区进行耦合计算. 但是区域分解法需要清晰划分不同的区域来准确模拟热辐射过程. 而对于实际物理问题, 确定具体区域边界具有很大的难度. 同时不同区域之间的边界条件处理也是一个难点. 第2类方法是通过对辐射强度的分解, 建立相应的介观分量方程和宏观分量方程, 混合传输-扩散(hybrid transport-diffusion, HTD)[27]模型实现两个方程的解耦, 提升计算效率. 第3类方法是采用渐进保持格式[2], 通过在网格单元上进行直接离散建模, 可获得不同尺度的辐射问题解, 不必依据光子自由程进行网格处理. 该类方法, 对于不同的尺度辐射问题, 可采用统一的求解格式, 同时相比于前两类方法也没有不同模型耦合的困难. 因此, 渐进保持格式在多尺度辐射计算中获得了很大的成功.

为实现多尺度问题的求解, 近年来, 在计算流体力学领域, 研究者们提出了基于Boltzmann方程求解的各类气体动理学格式[2,28-31], 如LBM, GKUA, UGKS, DUGKS, UGKWP等. 辐射是以光速进行传输能量的, 在一般工程应用中认为热辐射是稳态问题, 对于稳态问题, 上述气体动理学格式是通过时间演化更新瞬态计算直至收敛到稳态计算结果的, 而我们希望实现一种直接针对稳态辐射传递方程进行计算, 得到稳态辐射问题结果的方法.

本文的研究关注上述格式中的DUGKS格式, 该格式简洁易行, 主要包含通过沿着特征线的离散实现界面通量的重构和通过时间离散实现时间层的更新两个部分, 具体的实现过程在文献[1]中有详细的叙述. 最近, Zhou等[32]按照气体动理学格式的核心思想, 发展了求解稳态中子输运问题的SDUGKS格式. 本文从描述辐射输运现象的气体动理学理论模型方程出发, 进一步发展了SDUGKS格式应用于稳态多尺度辐射问题计算分析, 实现对稳态方程的直接求解. SDUGKS同DUGKS的整体步骤一致, 包含单元界面重构和单元更新两个内容. 其中单元界面重构与DUGKS一致, 采用沿特征线的积分离散并插值求解, 在单元数据更新上则采用隐式增量格式进行计算, 通过两部分的循环迭代, 达到稳态辐射的收敛解.

本文将首先介绍稳态辐射模型, 并详细论述SDUGKS格式离散和求解多尺度热辐射问题的过程, 通过多组典型多尺度热辐射算例验证模型的正确性和有效性. 在此基础上, 通过改变衰减系数分布构造不同物理尺度算例和多尺度算例进一步检验SDUGKS格式的渐进保持性质. 最后, 总结SDUGKS格式求解多尺度热辐射问题的优缺点, 并给出今后的发展方向.

1 基于DUGKS格式的稳态热辐射模型

1.1 稳态辐射传递方程及角度离散

稳态辐射传递方程为

式中,I(x,s) 表示空间位置x，s方向上的辐射强度值;βx为空间位置x处的衰减系数;S为源项, 写为

式中,ω为散射反照率;Ib=σT4为黑体辐射强度,σ为Stephen Boltzmann常数,T为介质温度;Φ为散射项函数, 对于各向同性散射Φ(s′,s)=1.

漫灰边界的辐射边界条件为

I(xw,s)=εwIb(xw)+

式中,εw为壁面发射率,ρw为漫反射率,nw为壁面的单位内法向量.

采用离散坐标法进行角度离散后的辐射传递方程为

S(x,sk)=(1-ω)Ib(x)+

其中，下标k,m{1,2,…,M}表示离散角度的序号，sk,sm分别表示离散角方向和散射入射方向.

利用Gauss公式, 上式可改写为

=-βxI(x,sk)+βxS(x,sk)

其中，Vc表示单元体积，nf表示单位外法线向量，xf表示单元界面位置.

1.2 基于DUGKS格式的界面重构

SDUGKS采用对离散辐射传递方程沿着特征线积分的方案实现空间离散来求解界面值

-βxf-L·sk(I(xf-L·sk,sk)-S(xf-L·sk,sk))

整理上式, 界面值为

I(xf,sk)=I+(xf-L·sk,sk)

其中,I+为辅助辐射强度分布函数

I+(xf-L·sk,sk)=Lfβxf-L·skS(xf-L·sk,sk)+

[1-Lfβxf-L·sk]I(xf-L·sk,sk)

即

I+(x,sk)=LfβxS(x,sk)+[1-Lfβx]I(x,sk)

对周围单元的I+值斜率的加权可得

式中

辐射强度存在明显的方向性, 可能存在相对于空间位置的大梯度变化. 因此, 需要采用Van Leer限制器格式来抑制可能出现的数值振荡

为了保证界面重构的辐射强度为正值, 积分长度L的取值不超过半网格尺度和最大衰减系数的倒数, 即

1.3 隐式增量格式的单元更新

为了对单元上的辐射强度进行更新, 可采用如下增量形式的辐射传递方程

(1)

式中, Resl(xi,sk)为残差源项

Resl(xi,sk)=βxSl(xi,sk)-βxIl(xi,sk)-

最后, 通过下式对中心处辐射强度进行更新

Il+1(xi,sk)=Il(xi,sk)+ΔIl(xi,sk)

求解方程(1), 需采用如下的迎风格式

当f=1时, 为阶梯格式,f=0.5时, 为菱形格式[13].

将上述迎风格式代入增量形式的辐射传递方程, 沿辐射传递方向依次扫描计算可以获得单元中心辐射强度的增量[12-13].

1.4 求解步骤

SDUGKS求解稳态辐射传递方程的步骤如下

(1)初始化辐射强度分布函数;

(2)计算辅助辐射强度分布函数值;

(3)插值计算特征位置处的辅助辐射强度分布函数值;

(4)重构单元边界上的辅助分布函数值;

(5)计算残差源项;

(6)求解增量形式的辐射传递方程;

(7)更新单元中辐射强度值, 返回第(2)步迭代计算, 直至误差满足收敛标准.

2 数值算例

选用多组一维和二维多尺度热辐射算例, 通过与谱元法(spectral element method, SEM)、 阶梯格式(STEP)、 解析解(exact)和细网格基准解(bench-marks)对比, 验证SDUGKS格式的正确性和有效性及其渐进保持性质.

2.1 一维算例

如图1所示, 对于无限大平板之间的热辐射问题, 左侧为热边界, 右侧为冷边界.

图1 一维模型图Fig. 1 One-dimensional model

2.1.1 辐射热平衡问题

对于均匀纯吸收性介质的辐射热平衡问题, 光学厚度τ取0.01, 0.1, 1和10, 空间和角度离散网格为20×20.图2给出了SDUGKS格式获得的辐射热流与温度分布. 从图2(a)中可以看出, 随着光学厚度的增加无量纲热流减小. 且在不同光学厚度下, SDUGKS计算的辐射热流与SEM的计算结果十分吻合, 它们之间的最大相对误差为1.263 6%. 从图2(b)中可以看出, 在介质与壁面处存在温度的不连续, 随着光学厚度的增大这种不连续将减少. 此外, 两种方法获得的无量纲温度的最大相对误差为0.419 3%. 因此, 对于该类热辐射问题, SDUGKS格式是准确可靠的.

2.1.2 介质温度恒定的热辐射问题

在纯吸收介质的介质温度恒定的情况下, 可获得辐射传递方程的精确解. 空间和角度离散网格数均为20.图3给出了不同光学厚度下SDUGKS格式、 离散坐标法和解析解获得的无量纲辐射热流分布. 随着光学厚度的增大, 无量纲辐射热流的梯度增大, 当τ取1和10时, 阶梯格式结果出现了明显的误差, 而SDUGKS格式的计算结果与精确解吻合良好. 对于无量纲辐射热流分布, 阶梯格式与精确解的最大相对误差分别为6.152 9%和12.383 7%, 而SDUGKS格式和精确解的最大相对误差分别为1.434 4%和3.104 2%. 可见, SDUGKS格式展现出具有多尺度渐进保持性. 即能够在同一套网格下, 获得不同光学厚度的渐进解.

(a) Distribution of dimensionless radiative heat flux

(b) Distribution of dimensionless temperature图2 一维辐射热平衡计算结果Fig. 2 Numerical results of one-dimensional equilibrium radiation

图3 不同方法获得一维纯吸收介质的无量纲辐射热流Fig. 3 Dimensionless radiative heat flux of one-dimensional purely absorbing media by different methods

2.1.3 衰减系数随空间变化的一维热辐射问题

对于一维无散射热辐射问题, 选取其衰减系数空间分布如下

式中,α为调节参数, 用于调节控制衰减系数的变化范围. 在本算例中,α分别取1 000和100.

图4给出了空间网格数和角度离散方向数均为20条件下的无量纲温度分布. 为了验证SDUGKS格式的正确性, 将超密网格(空间网格2 000和角度离散200)下的阶梯格式结果作为基准解. 从图中可以看出, 在衰减系数较大以及衰减系数变化剧烈的区域, SDUGKS比DOM的精度更高. 不需要随着光学厚度的变化而改变网格, 就能够实现多尺度热辐射计算.

(a) Distribution of dimensionless incident radiation (α=1 000)

(b) Distribution of dimensionless incident radiation (α=100)图4 衰减系数随空间变化的一维多尺度热辐射结果对比Fig. 4 Comparison of one-dimensional multiscale radiation heat transfer results for the case of attenuation coefficient variation with spatial position

2.2 二维算例

如图5所示, 考虑二维热辐射问题. 底部边界为热边界, 其余边界壁面为冷壁面, 内部介质为冷介质. 经过网格无关性验证后, 空间网格和角度网格分别为Nx×Ny=30×30和Nφ×Nθ=40×40.

图5 二维多尺度热辐射模型Fig. 5 Two-dimensional multiscale radiation model

2.2.1 二维均匀介质内热辐射

图6给出了均匀介质的无量纲温度和无量纲热流分布图. 从图中可以看出, 无量纲温度和无量纲热流均呈现左右对称分布. 这是因为, 辐射从底边边界均匀地向空间传输, 又因为方形空间的左右对称性, 所以辐射能在空间上呈现左右对称分布特征.图6(d)给出了SDUGKS与间断谱元法的结果对比. 两者之间的最大相对误差为3.37%.

(a) Contour plot of dimensionless temperature

(b) Contour plot of dimensionless radiative heat flux in x direction

(c) Contour plot of dimensionless radiative heat flux in y direction

(d) Comparison of radiative heat flux in y direction at upper boundary图6 二维辐射计算结果Fig. 6 Numerical results of two-dimensional radiation

2.2.2 二维多尺度热辐射

如图5所示, 小矩形内外衰减系数不等导致二维多尺度热辐射问题.图7, 8分别给出了衰减系数为β1=10,β2=0.01和β1=100,β2=0.01时无量纲温度和无量纲辐射热流分布.

(a) Contour plot of dimensionless temperature

(b) Contour plot of dimensionless radiative heat flux in x direction

(c) Contour plot of dimensionless radiative heat flux in y direction

(d) Distribution of incident radiation at x/L=0.5图7 β1=10和β2=0.01时二维辐射计算结果Fig.7 Numerical results of two-dimensional multiscale radiation for β1=10 and β2=0.01

从图7(a), (b)和(c)可以看出, 无量纲温度和无量纲辐射热流分布云图呈左右对称分布.图7(d)给出了x/L=0.5处y方向上的无量纲投入辐射.图7(d)中, 中间区域的投入辐射是大梯度变化的, 而其上下两侧投入辐射变化平缓, 这反映了光学厚区的强吸收作用. 最上段的投入辐射分布相对下段很小, 说明在光学厚区域, 热辐射遮蔽效应明显. 与均匀介质内热辐射对比发现, 多尺度介质的投入辐射在光学厚区域的前端形成了一个热区, 其值甚至高于热边界附近的投入辐射值. 这是因为在光学厚区域, 介质将辐射能吸收之后再发射出来, 有明显的辐射增强效果. 此外, 将SDUGKS格式的计算结果、 阶梯格式的结果以及基准解对比发现, SDUGKS格式的计算结果更为准确.

进一步增加中间区域的衰减系数, 考虑β1=100和β2=0.01情况. 如图8(a), (b)和(c), 在空间上存在大衰减系数差的情况, SDUGKS仍然能够获得较好的计算结果. 从图中可以看出, 中间衰减系数较大区域具有明显的辐射遮蔽效应, 经过中间区域后的温度和热流显著低于其他区域. 从图8(d)可以看出,x/L=0.5处y方向上的无量纲投入辐射变化更加剧烈. 光学厚区域前端的热区更加明显, 多尺度的辐射增强效果更加显著. 此外, 相同的网格下, SDUGKS格式的计算结果比阶梯格式结果更加接近基准解.

(a) Contour plot of dimensionless temperature

(b) Contour plot of dimensionless radiative heat flux in x direction

(c) Contour plot of dimensionless radiative heat flux in y direction

(d) Distribution of incident radiation at x/L=0.5图8 β1=100和β2=0.01时二维辐射计算结果Fig.8 Numerical results of two-dimensional multiscale radiation for β1=100 and β2=0.01

3 结论与展望

本文针对稳态多尺度辐射传热问题提出了稳态离散统一气体动理学格式(SDUGKS), 详细阐述了SDUGKS离散和求解多尺度热辐射问题的过程, 通过沿特征线积分进行辐射传递方程的离散实现单元界面重构, 通过辐射传递方程的守恒性进行单元数据的隐式更新, 实现了对辐射传递方程的多尺度求解. 通过多组一维和二维多尺度稳态热辐射算例检验了SDUGKS格式的正确性和有效性, 论证了SDUGKS的渐进保持性质. 并且, 通过与离散坐标法计算结果对比可知, 较单一尺度方法， SDUGKS求解多尺度问题时具有明显的优势.

高温热辐射在高温传热中所占比重愈加重多, 且因为辐射过程与光谱波长、 温度、 压力等密切相关, 存在明显的多尺度热辐射. 而SDUGKS在稀疏网格下, 仍然能够实现对各种辐射问题的准确计算, 有较好的计算精度和渐近保持性质, 这对于减少计算资源具有重要作用. 此外, SDUGKS在多尺度热辐射问题的研究成果将为研究多尺度条件下辐射与导热、 对流耦合传热问题提供便利.

后续的工作将从SDUGKS格式中发展高精度无振荡的插值方法, 与插值精度相适应的边界条件处理方法, 实现高精度单元数据更新方式, 并在波长相关的热辐射问题应用等方面展开.