为了实现高效率、低污染物燃烧, 采用基于自发辐射的燃烧诊断技术对于认识燃烧现象、研究燃烧机理具有重要的意义[1]。基于自发辐射的燃烧诊断技术是一种被动式光学诊断技术, 它具有对诊断对象无干扰、响应速度快、系统简单等优点[2], 近年来倍受国内外学者关注。但在该技术中, 关键点是通过介质辐射传递方程建立燃烧介质发射与观测辐射图像之间的关系, 获取介质在角度和空间方向上高分辨率的辐射强度分布。

辐射传递方程用来描述辐射强度在介质中吸收、发射和散射过程。它是一种具有强对流数值特性的微分积分方程。在大多数实际工况下, 它都无法获得解析解, 只能借助数值方法获得数值解[3-4]。经过几十年发展, 主要数值方法有基于微分形式辐射传递方程全局离散的离散坐标法(DOM)[5]、有限体积法(FVM)[6]和有限元法(FEM)[7], 以及基于射线跟踪的蒙特卡罗法[8]等。但是, 这些方法多以计算辐射热流的辐射积分项为主。为获得辐射强度在角度和空间上的高分辨率分布, 需要发展一种既能够实现辐射强度在空间和角度上高阶离散, 又能够方便、快速计算的数值方法。

微分单元法[9-11]借鉴了有限元思想, 利用等参元生成高阶导数，并依据边界元思想, 对节点方程组进行配齐。它可实现物理量的高阶表达, 并且容易获得控制方程和边界条件的配置点方程。在实施过程中, 较为容易产生稳定的结果, 在节点配置方面不需要任何积分。目前, 已经广泛应用于传热学[9-10]、固体力学[11]等领域。

本文将利用微分单元法的上述优点, 分析圆柱内各向异性散射介质高温热辐射的数值特性, 构建相应热辐射模型, 实现辐射强度在角度和空间上的三维高阶离散。通过带有解析解的算例, 检验模型在角度和空间上辐射强度的高精度刻画; 通过与蒙特卡罗法对比, 证明微分单元法对各向异性散射介质中辐射强度计算的准确性; 进一步对比分析采用迎风方案前后, 圆柱内辐射强度在空间和角度上的三维分布, 验证迎风方案对数值振荡抑制的有效性。

轴对称圆柱坐标系下辐射强度为r, φ, θ的函数[12-15]。各向异性散射介质的热辐射传递方程可简化为

式中: μ=sinθcosφ, η=sinθsinφ, ξ=cosθ为方向余弦; β=κa+κs为衰减系数; Ib为黑体辐射强度; Φ(Ω, Ω')为表征散射介质特性的散射相函数, 可采用Legendre多项式表示为[16-19]

式中: pl为Legendre多项式; dl为相应系数, 当dl为0时, 即为各向同性散射介质, 本文采用5阶的Legendre多项式作为各向异性散射介质的散射相函数。

漫灰壁面的热辐射边界条件为

式中, ε为壁面发射率。

热辐射源项为

投入辐射为

微分单元法是通过在等参元内插入节点的方法对方程进行离散。采用拉格朗日插值基函数作为等参元内形函数

式中: m为等参元内节点个数; χ为节点坐标。

等参元内辐射强度可离散为

辐射强度一阶偏导数可离散为

式中, 形函数关于坐标的偏导数为

式中, J为雅可比矩阵。

微分单元法的离散节点可根据在等参元中的位置分为内部点、连接点、边界点3类。对于内部点的信息, 可采用单元内控制方程直接求解; 对于连接点的信息, 需通过连接点周围单元共同确定; 对于边界点的信息, 通过边界条件直接确定。对于热辐射虚拟边界上节点的信息, 不能通过边界点的处理方法获取, 而需要通过连接点的处理手段获得。

辐射传递方程仅存在辐射强度对各个方向的一阶偏导, 是一种没有扩散项的对流方程, 其对流项容易引起数值不稳定, 产生非物理振荡。因此, 需要在r和φ方向上采用迎风方案, 来抑制非物理振荡。如图 1所示, 当辐射强度由R1传向R2时(φ∈[0, π/2]), 靠近R1的单元为上游单元; 反之, 靠近R2的单元为上游单元。此外, 无论辐射强度是由R1传向R2还是R2传向R1, φ较大位置总是上游单元。

对于内部点, 辐射传递方程的离散形式为

如图 2所示, 对于连接点, 采用迎风方案后, 辐射传递方程采用上游单元离散为

式中: a, b, c分别表示r方向上游单元数、φ方向上游单元数、r和φ方向上共有上游单元数。

对于真实边界上的节点, 辐射强度可直接根据热辐射边界条件确定。对于虚拟边界上的节点, 辐射强度需要采用内部点和边界点相同的离散方法。在φ=π/2处存在一簇强间断奇异点。采用双层节点方案后, 奇异点上的辐射强度分为I+和I-。对于r=R1, I+为真实边界上的辐射强度, I-为虚拟边界上的辐射强度; 对于r=R2, 则相反。双层节点方案保证了间断点分别在φ∈[0, π/2]和φ∈[π/2, π]2个区间内进行了2次计算, 并且在各自区间内依据物理模型采取合适的迎风方案, 保证了辐射强度在间断点处的高精度分辨率。

基于微分单元法的圆柱内各向异性介质中热辐射模型的求解步骤为:

1) 确定物性参数及网格数、离散空间变量和角度变量, 获得对应节点。

2) 对内部点、连接点和边界点上的辐射传递方程进行离散, 获取相应系数矩阵。

3) 对内部点、连接点和边界点上的源项进行计算, 并引入热辐射边界条件。

4) 求解各点上的辐射强度, 并更新投入辐射。

计算前后2次投入辐射的最大相对误差, 如果最大相对误差小于收敛标准, 则终止迭代并输出结果; 否则, 返回步骤3)进行迭代计算。

为了验证基于微分单元法的各向异性散射介质热辐射模型的正确性, 需要分别对带有解析解的各向同性介质的特殊算例和各向异性介质的一般算例进行对比。

在Φ(Ω, Ω')=1的各向同性散射介质中, 取圆柱内外半径比为1∶3, 壁面边界条件为解析解条件下的定值。将辐射源项采用解析解化简, 其辐射传递方程描述为

(12) 式的解析解为

表 1给出了不同单元和节点数下, 微分单元法的数值解和解析解之间的最大相对误差。从表中可以看出, 随着单元和节点数的增加, 最大相对误差逐渐减少。但当单元和节点数超过4×4×4时, 最大相对误差变化不明显。因此, 本文选用单元和节点数均为4×4×4进行计算。另外, 此时辐射强度的最大相对误差为0.019 2%, 这表明基于微分单元法的热辐射模型有高阶精度。

本文进一步选用文献[20]中采用蒙特卡罗法计算的圆柱内各向异性散射介质热辐射算例, 验证微分单元法辐射模型的正确性。在该算例中, 内壁面R1为冷壁面, 其辐射强度为0;外壁面作为辐射热源I0=I(R2)=1;壁面发射率均为1, 内外半径比为1∶2;衰减系数、吸收系数、散射系数分别为β=0.1 m-1, κa=0.05 m-1, κs=0.05 m-1; 各向异性散射为瑞利散射, 相应散射相函数为

微分单元法所采用的网格数为31×31×31, 蒙特卡罗法[20]采取的网格数为201×201×201, 对比θ=π/2时辐射强度在r-φ平面上的分布。从图 3可以看出, 微分单元法在较少的单元数下获得了较好的计算结果, 并对强间断点处的无量纲辐射强度I/I0进行了有效捕捉。其中, I0为外壁面的初始辐射强度值。

针对辐射热平衡问题, 内外径比为R1/R2=0.5;衰减系数为β=1 m-1; 散射反照率为ω=0.8。边界条件为壁面发射率为εw=0.5, 内外径温度比为T1/T2=0.5的漫反射边界。图 4给出了采用迎风方案前后无量纲辐射强度I/I0对比, 其中I0为T2处辐射强度。从图中可以看出, 采用迎风方案后, 无量纲辐射强度分布的数值振荡得到了有效抑制。这是由于无量纲辐射强度的传递与光的传播方式相同, 具有很强的方向性。如果不采用迎风方案, 上下游等参元都参与计算, 这会导致连接点上无量纲辐射强度出现峰值, 引起数值振荡, 产生较大的误差。因此, 迎风方案可以有效处理辐射传递方程存在的强对流数值特性问题。

图 5进一步给出了外壁面上无量纲辐射强度在角度上的分布。虽然无量纲辐射强度在φ=π/2方向上发生了强间断, 但是微分单元法依然可以实现辐射强度在角度上的高分辨率刻画, 并且对间断进行有效捕捉。

针对非辐射热平衡问题, 考虑内外径比R1/R2=0.5, 内外径温度比T1/T2=1, 介质初始温度与内径温度比Tg/T1=10的情况。分析衰减系数、散射反照率和壁面发射率等参数对介质无量纲温度分布的影响规律。

如图 6所示, 随着衰减系数的增加, 介质温度快速上升, 而介质中间区域的温度梯度变化不明显。这是因为当散射反照率固定时, 衰减系数越大意味着参与性介质发射和吸收辐射的能量越强, 辐射传热更加剧烈, 介质的温度上升。而在介质内, 其初始温度相同, 介质之间的温度梯度较小, 因此辐射传热程度较低, 介质整体温度基本保持一致。

如图 7所示, 随着散射反照率的增加, 介质温度不断降低。这是因为散射反照率表征散射辐射能在能量衰减中所占比重, 散射反照率越大, 更多能量被散射, 吸收的辐射能随之减少。而在辐射源项中, 散射的能量是由空间点之间角向辐射强度传递的, 其总量不会改变; 吸收的能量是由介质的黑体辐射力产生的, 散射反照率的减少意味着介质由黑体辐射力产生的能量增多, 因此总的辐射源项增大, 介质温度升高。

图 8给出了壁面发射率对介质温度分布的影响规律。从图中可以看出, 介质温度会随着壁面发射率的增加而减小。壁面发射率为漫反射灰壁面边界条件中表征壁面反射辐射能强弱的参数, 决定了壁面发射出黑体辐射强度的大小, 其变化区间为[0, 1], 当壁面发射率为0时, 表示全反射壁面; 当壁面发射率为1时, 表示黑体壁面。壁面发射率越大, 在边界上更多的辐射能被吸收, 反射到介质中的辐射能减少。此算例中, 介质的温度高于壁面, 壁面反射的能量越大, 则介质温度越高; 因此介质温度会随着壁面发射的增大而减小。

针对不同的散射介质, 采用内外径比R1/R2=0.5, T1/T2=0.5的黑壁面; 衰减系数为β=1 m-1, 散射反照率ω=1的介质作为算例进行了分析。如图 9所示, 本文计算了5种散射介质, 辐射强度在空间点r=(R1+R2)/2, θ=0平面上无量纲辐射强度沿φ方向的分布。

在此算例中, 外壁面为高温壁面, 因此无论何种散射介质, 辐射强度均在φ=π方向上最大, 在φ=0方向上最小。从图 9中可以看出, 相较于各向同性散射介质, 在勒让德离散的各向异性散射介质影响下辐射强度较大, 而瑞利散射介质影响下的辐射强度较小。这是因为勒让德离散使得散射相函数在各个方向上均大于1, 增强了空间点上各个方向上散射的辐射能, 使得辐射强度增大; 反之瑞利散射使得散射相函数在各个方向上均小于1, 使得辐射强度降低。

Φ(Ω, Ω')=1+0.1cosφ的各向异性散射介质强化了向外壁面的辐射能散射, 而Φ(Ω, Ω')=1-0.1cosφ强化了向内壁面的辐射能散射。由于外壁面为高温, 向内壁面散射方向上的辐射能更大, 因此在前者的影响下, 辐射强度更高; 反之, 后者的影响使辐射强度降低。

本文采用微分单元法对圆柱内各向异性介质的热辐射进行分析, 得到以下结论:

1) 通过与解析解对比发现, 微分单元法求解圆柱内热辐射问题具有较高的计算精度, 可实现辐射强度在空间和角度上的高精度刻画。

2) 与文献中蒙特卡罗法结果相比, 微分单元法对于各向异性散射介质热辐射问题具有很好的计算精度, 并且可以对辐射强度在角度上进行精细刻画, 实现对间断点的有效捕捉。

3) 对于辐射传递方程的强对流特性, 采用迎风方案, 可有效抑制数值解的非物理性振荡, 获得稳定的计算结果。