燃烧是化石能源利用的主要方式，燃烧诊断技术的发展对于更好地认识燃烧和控制燃烧以实现碳中和目标具有重要的意义．基于辐射图像分析的燃烧检测技术是一种被动式的光学诊断技术，具有对被测对象无干扰、响应快、系统简单的特点而备受关注[1-5]．在该技术中，通过求解燃烧介质中的辐射传递过程建立燃烧空间介质发射和辐射图像间对应关系，是进行辐射图像分析获得温度场、浓度场等燃烧信息的基础．现有辐射传递过程求解方法大多以计算辐射传热量即辐射热流等辐射积分量为主，如球形谐波法、区域法、有限体积法、离散坐标法和蒙特卡罗法等．而在基于辐射图像分析的燃烧检测技术中，辐射图像每一个像素图像信息对应该像素方向上的辐射强度信息．为了分析高像素辐射图像，须要发展求解高方向分辨率辐射强度的辐射传递过程求解方法．文献[6-7]基于离散坐标法、文献[8]基于谱方法，发展了可以求解任意方向辐射强度的计算方法．文献[4]利用逆蒙特卡罗法求解方向辐射强度并用于实际的燃烧温度场测量．文献[5，9-10]基于蒙特卡罗法提出了高效计算方向辐射强度的被介质散射或壁面反射的能量分布份额(DRESOR)法．现有的高方向分辨率辐射强度计算方法都是针对灰性介质，适用于燃煤等燃烧空间存在大量固体颗粒的燃烧检测技术[1-5]．而对于燃气等燃烧空间主要为气体产物的应用场合，由于气体介质辐射的强烈光谱选择性，因此须要将气体辐射模型和辐射传递求解方法相结合，发展针对气体非灰介质的高方向分辨率辐射强度求解方法．现有气体辐射模型主要包括[11]：具有最高光谱分辨率的逐线法(LBL)模型；窄带模型，如统计窄谱带模型(SNB)等；宽带模型，如灰气体加权和模型(WSGG)等；全波段模型，如全光谱关联k分布模型(FSCK)[12]等．其中，逐线法具有最好的求解精度，但是同时需要大量的计算资源，因此仅应用于简单辐射问题求解．相较于LBL模型，SNB模型可以在大大减少计算资源同时保持较高精度，当LBL模型无法求解时，其计算结果常作为基准解．然而，SNB模型中给出的是窄带透射率而不是吸收系数，使得其只能应用于无散射性介质和边界无反射情况下辐射传递过程求解，大大限制了其应用范围[13]．将SNB模型和correlated-k相结合的SNBCK模型可用于求解散射性介质和反射性边界的辐射对象，文献[14]将DRESOR法与SNBCK模型相结合，成功求解含反射性壁面的辐射问题．SNBCK模型本质上是k-distribution模型，相较SNB模型更加复杂，应用广泛性远不及SNB模型．本研究尝试直接将SNB模型和DRESOR法相结合，以求解散射性非灰介质和反射性边界条件下的辐射传递过程．通过修改DRESOR法的求解过程，将吸收系数相关变量利用透射率进行替代，使其与SNB模型结合计算高分辨率方向辐射强度．首先介绍DRESOR法和SNB模型，以及DRESOR法和SNB模型的结合方法；然后以LBL模型为基准解，评判DRESOR法结合SNB模型的计算精度；最后得出相关结论．1 计算原理1.1　DRESOR法原理DRESOR法是一种求解高方向分辨率辐射强度的有效方法，该方法基于正向蒙特卡罗方法能束追踪的思想．不同于蒙特卡罗方法的是：蒙特卡罗法只能求解各点的辐射热流，而DRESOR法一旦求解得到DREOSR数，可以方便地计算任意点任意方向上的辐射强度．积分形式的辐射传递方程为Iη(r,ŝ)=Iwη(r,ŝ)exp(-∫0sβηds')+∫0sSη(r',ŝ)exp(-∫s'sβηds″)βηds',式中：下标η为谱带波数；Iη(r,ŝ)为在r位置处ŝ方向上的光谱辐射强度；s为积分路径总长；s'为0～s积分路径上的一点；s″为s'～s积分路径上的一点；βη为光谱消光系数，等于吸收系数kaη和散射系数ση之和；源函数Sη(r',ŝ)和边界条件Iwη(r,ŝ)可分别写为Sη(r',ŝ)=(1-ωη)Ibη(r')+ωη4π∫4πIη(r',ŝi)Φη(ŝi,ŝ)dΩi, (1)Iwη(rw,ŝ)=ε(rw)Ibη(rw)+∫n̂×ŝ0ρ(rw,ŝi,ŝ)Iη(rw,ŝ)n̂×ŝdΩ', (2)其中，ωη为散射率，rw为壁面位置，ε(rw)为壁面发射率，Ibη为黑体辐射强度，Φη(ŝi,ŝ)为散射相函数，ρ(rw,ŝi,ŝ)为壁面反射相函数．当采用DRESOR法求解辐射传递方程时，将原辐射强度以一种分能量贡献份额的方法进行了重新表示．对于计算点r在ŝ方向上的辐射强度Iη(r,ŝ)，等于各个其他单元(包括面元和体元)到达该点的直接贡献与间接贡献之和．DRESOR法中辐射能量传递过程示意图如图1所示，A为能量传递过程中直接贡献点网格，B为间接贡献点网格，C为待计算辐射强度位置，对C点的辐射能量贡献积累包含直接贡献单元(如A网格)发射经过衰减达到C的直接贡献，以及所有网格发射经过间接贡献单元(如B网格)被散射或反射再到达C点的间接贡献．所有网格发出到达B点被散射或反射的间接贡献份额采用能束统计的方式得到，用DRESOR数进行表示．10.13245/j.hust.220501.F001图1DRESOR法辐射能量传递过程示意图为了更加细致地描述DRESOR法分能量贡献计算辐射强度的思想，本研究将某位置某方向辐射强度分成6项进行讨论，对于各项同性散射和漫反射问题，常见的各项DRESOR法离散格式通常可分别表示[15]为I1(r,ŝ)=(1/π)[πε(rw)Ibη(rw)]τη(A→C)；I2(r,ŝ)=1π(∑iwπε(rw)Ibη(rw)⋅Rds(rw,rw',ŝ)Aw)τη(B→C);I3(r,ŝ)=1π(∑is4πkaηIbη(r)⋅Rds(r,rw',ŝ)Vs)τη(B→C); (3)I4(r,ŝ)=14π∑is4πkaηβηIbη(r)⋅(τη(A-→C)-τη(A+→C)); (4)I5(r,ŝ)=14π∑is∑iwπβηε(rw)Ibη(rw)⋅Rds(rw,r',ŝ)Aw(τη(B+→C)-τη(B-→C)); (5)I6(r,ŝ)=14π∑is∑is4πkaηβηIbη(r)Rds(r,r',ŝ)⋅Vs[τη(B+→C)-τη(B-→C)], (6)式中：I1为面元的直接贡献；I2为面元经面元反射的间接贡献；I3为体元经面元反射的间接贡献；I4为体元的直接贡献；I5为面元经体元散射的间接贡献；I6为体元经体元散射的间接贡献；iw为壁面单元；is为空间体积单元；Aw和Vs分别为壁面单元面积和空间单元体积；Rds为DRESOR数，用于计算式(1)和(2)中的积分部分．DRESOR数表示某个单元经另外一个单元散射/反射的能量份额，如Rds(r1,r2,ŝ)表示r1单元发出被r2单元在ŝ方向上散射/反射的份额．DRESOR数基于路径长度法计算得到，具体可以参考文献[9]；一旦计算得到DRESOR数，可以很方便计算出任意点在任意方向上的辐射强度．总辐射强度Iη(r,ŝ)可以通过各个项能量贡献相加得到，即Iη(r,ŝ)=I1(r,ŝ)+I2(r,ŝ)+I3(r,ŝ)+I4(r,ŝ)+I5(r,ŝ)+I6(r,ŝ).2.2　SNB气体辐射模型SNB模型是一种提供气体窄带平均透射率的谱带模型，将光谱划分为一定波数范围的窄带．对于Malkmus线强分布函数结合Lorentz谱线线形函数，SNB模型计算窄带平均透射率的解析表达式[16]为τ¯=exp-(β¯/π)1+2πfplk¯/β¯-1，式中：f为气体的摩尔分数；p为气体总压；l为几何路径长度；k¯为平均线强与线间距的比值或平均吸收系数；β¯为谱线重叠参数，通常可表示为β¯=2πγ¯/δ¯，其中γ¯和δ¯分别为平均谱线半宽和平均谱线间距．SNB模型的计算精度主要取决于模型参数：k¯和β¯，或k¯,γ¯和δ¯．生成SNB模型参数的方法主要有参数定义法、非线性最小二乘拟合方法和K-moment方法，本研究所采用的SNB模型参考文献[17]提供的最小二乘拟合方法生成．在最小二乘拟合方法中，将SNB模型计算的窄带透射率随路径长度变化的曲线与LBL模型计算的精确窄带平均透射率随路径变化曲线进行拟合，不断调节SNB模型参数，使两条曲线整体误差的平方和最小．SNB模型是基于均温均质气体路径下的LBL模型拟合得到的，不能直接用于非均匀介质中，须要引入近似方法进行处理．目前针对SNB模型发展了Curtis-Godson近似和Lindquist-Simmons近似[17]，本研究将结合Curtis-Godson近似处理SNB模型非均匀气体介质问题．引入Curtis-Godson近似后非均匀路径透射率可以表示为τ¯CG(s1,s2)=exp-β¯CGπ1+2πk¯CGvβ¯CG-1，式中：s1和s2分别为非均匀路径的起点和终点；v=∫s1s2f(s)p(s)ds；k¯CG=1v∫s1s2f(s)p(s)k¯(s)ds；β¯CG=1vk¯CG∫s1s2f(s)p(s)k¯(s)β¯(s)ds，其中f(s)，p(s)，k¯(s)和β¯(s)分别为路径上各空间位置处的气体介质摩尔分数、总压、窄带平均吸收系数和谱线重叠参数．2.3　DRESOR法和SNB模型结合当体元网格划分足够小、光学厚度足够薄时，有kaηΔl≈1-exp(-kaηΔl)；(7)1-exp(-βηΔl)≈βηΔl，(8)式中：kaη为网格吸收系数；βη为网格衰减系数；Δl为网格划分长度．将式(7)和(8)代入式(3)~(6)中，则I3~I6计算表达式可以分别表示为：I3(r,ŝ)=1π∑is4π1-exp(-kaηΔl)]Ibη(r)⋅Rds(r,rw',ŝ)Vsτη(B→C)/Δl; (9)I4(r,ŝ)=14π∑is4π1-exp(-kaηΔl)]⋅Ibη(r)τη(A→C); (10)I5(r,ŝ)=14π∑is∑iwπε(rw)Ibη(rw)⋅Rds(rw,r',ŝ)Awτη(B→C)Δl;I6(r,ŝ)=14π∑is∑is4π[1-exp(-kaηΔl)]⋅Ibη(r)Rds(r,r',ŝ)Vsτη(B→C). (11)在DRESOR法与SNB模型结合过程中，为了避免两段能束传递过程中透射率计算引入额外误差问题，对DRESOR法能束追踪过程进行了修改．在能束追踪过程中，将式(9)~(11)中(1-e-kaηΔl)与能束追踪过程相结合，DRESOR法中所有追踪能束起始点能量份额从1改变为(1-e-kaηΔl)，经过一段距离L后的能束剩余能量为[e-kaηL-e-kaηΔl+L]．另外，在原来的DRESOR数计算过程中，只统计每条能束从A点发出到B点被散射或被反射的能量份额，当计算辐射强度时再将所有能量汇集计算B点到C点的能量衰减．而在本方法中，能束行进到B点时直接计算经散射或反射到达C点的能量份额，将能量份额储存到DRESOR数中，再用该DRESOR数进行辐射强度计算．综上，当使用DRESOR法结合SNB模型时，各项最终离散格式表达式分别为：I1(r,ŝ)=1π[πε(rw)Ibη(rw)]τ(A→C)；I2(r,ŝ)=1π∑iwπε(rw)Ibη(rw)Rds(rw,rw',ŝ)Aw；I3(r,ŝ)=1π∑is4πIbη(r)Rds(r,rw',ŝ)Vs/Δl；I4(r,ŝ)=14π∑is4π1-exp(-kaηΔl)]⋅Ibη(r)τη(A→C);I5(r,ŝ)=14π∑is∑iwπε(rw)Ibη(rw)⋅Rds(rw,r',ŝ)AwΔl;I6(r,ŝ)=14π∑is∑is4πIbη(r)Rds(r,r',ŝ)Vs．可以看出：一旦通过能束跟踪过程计算得到DRESOR数，再结合SNB模型中给出的透射率结果，各位置任意方向光谱辐射强度均可计算得到．3 结果与分析本研究中所有计算工况均为充满介质的一维平行平板，两平板间间距为1.0 m，考虑的辐射气体主要为CO2和H2O．计算谱带范围为250~6 000 cm-1，涵盖了CO2和H2O大部分强发射谱带．在本节中，将采用DRESOR法结合SNB模型求解均匀、非均匀介质以及不同散射颗粒与壁面条件工况下的辐射传递问题，并以LBL模型计算结果作为基准解，验证所提出方法的准确性．3.1　等温均匀介质本节中的平行平板充满等温均匀CO2和H2O气体介质．在工况1中，介质温度为1 500 K，摩尔分数分别为fH2O=0.1，fCO2=0.1，壁面边界条件为冷黑壁面．工况2中气体介质温度为1 250 K，摩尔分数分别为fH2O=0.25，fCO2=0.15；边界为漫反射壁面，发射率/反射率为0.5/0.5，壁面温度为750 K；工况2在介质中加入散射颗粒，忽略颗粒的吸收和发射，颗粒的散射系数采用文献[18]使用的简单方法，假设颗粒的散射特性是各向同性，颗粒的散射系数计算公式为ση=πa2Nt2-4msin m-1-cosmm-(1-e-0.1u)，式中：a为颗粒粒径；Nt为颗粒数密度；m=2u(n-1)，其中u=2πa/λ为颗粒的尺寸参数，λ为波长．在本小节中球形颗粒粒径为50 μm，复折射率n=1.3，颗粒数密度为2.0×108 m-3．计算工况如表1所示．10.13245/j.hust.220501.T001表1等温均匀介质工况参数工况fCO2fH2O介质温度/K壁面温度/ K壁面发射率散射颗粒10.100.101 50011.0无20.250.151 2507500.5有图2展示了两个工况在x=1.0 m位置处垂直于平行平板方向光谱辐射强度及相对误差计算结果，图中：Iv为窄带辐射强度积分；Re为与LBL基准解的相对误差．2 100 cm-1和3 500 cm-1附近为两个辐射强度波峰光谱区域，在两个工况中，DRESOR法结合SNB模型在该高发射谱带内计算精度较高，误差约为1%．工况1中具有明显发射的谱带最大误差约为8%，工况2中考虑加入散射颗粒以及壁面边界条件影响，相较于工况1在大部分谱带区域误差水平有所下降，大部分谱带区域计算误差在3%以内，这是由于使用DRESOR法结合SNB计算各能量间接贡献较为准确．10.13245/j.hust.220501.F002图2工况1和工况2右边界(x=1.0 m)垂直平板方向辐射强度及相对误差图3为工况1和工况2在平行平板各位置各方向光谱积分辐射强度分布图，图中θ为方向．计算值为250~6 000 cm-1光谱辐射强度积分．由于计算对象的对称性，因此方向辐射强度结果也呈现对称性，0°方向上随着检测点向右移动，辐射强度有逐渐增大的趋势；在180°方向上呈现相反的特性，随着检测点向右移动，辐射强度逐渐减小．工况1由于工况中光学厚度较小，辐射强度随方向变化梯度较大，高辐射强度区域集中在各位置90°方向附近．相较于工况1，由于加入了介质散射的影响，因此工况2整体上各个方向处的辐射强度变化更为平缓，高辐射区域位于靠近平板中心的位置，且其范围更大．10.13245/j.hust.220501.F003图3工况1和工况2不同位置和方向光谱积分辐射强度分布3.2　非等温非均匀介质本节中1.0 m平行平板充满非均温非均匀CO2和H2O气体介质，介质的温度和摩尔分数分布情况分别采用文献[17]的函数，即T(x)=1 500-1 000(2x/L-1)2；f(x)=0.8(1-x/L)(x/L),式中：L为平行平板间距；T(x)和f(x)分别为不同位置处温度和介质摩尔分数．工况3中设置壁面条件为冷黑壁面，介质中无散射颗粒．工况4在介质中加入散射颗粒，散射颗粒参数设置同3.1节，同时综合考虑壁面的影响，壁面温度为750 K，壁面发射率/反射率为0.5/0.5．非等温非均匀介质工况参数如表2所示．10.13245/j.hust.220501.T002表2非等温非均匀介质工况参数工况fCO2fH2O介质温度壁面温度/K壁面发射率散射颗粒3非均匀非均匀非均匀11.0无4非均匀非均匀非均匀7500.5有图4展示了工况3和工况4的计算结果，在本节工况中，介质温度从500 K变化到1 500 K，介质摩尔分数从0.0变化到0.2．工况3在各个中高发射谱带计算精度良好，大部分谱带辐射强度误差在5%以内．工况4相较于工况3，加入介质散射和壁面条件因素后边界处光谱辐射强度更大，其误差曲线波动趋于平稳，计算精度表现略优于工况3．工况3和工况4的最大误差出现在2 400 cm-1谱带附近，该谱带区域辐射强度较弱，对整体计算精度影响较小．10.13245/j.hust.220501.F004图4工况3和工况4右边界(x=1.0 m)垂直平板方向辐射强度及相对误差图5展示工况3和工况4各位置各方向光谱积分辐射强度变化趋势．工况3各位置处在90°方向附近有较大的变化梯度，且从中心向两侧梯度变化有所增大．在平板中心区域辐射强度较大，主要是因为平板中心处温度和浓度相较于两侧更高．工况4在加入散射因素后各位置梯度变化明显更加趋于平缓，不同位置、方向上的辐射强度差异减小，高辐射强度区域集中于平板中心位置接近90°方向，相较于工况3其高辐射强度区域范围增大，最大辐射强度数值降低．工况3和工况4在平板两侧处均存在低辐射强度区域，且低辐射强度区域面积较均匀介质工况更大，其原因在于两侧为低温低浓度介质区域．10.13245/j.hust.220501.F005图5工况3和工况4不同位置和方向光谱积分辐射强度分布3.3　不同散射颗粒辐射系统本节探究在充满非均温非均匀CO2和H2O气体介质的一维平行平板中，添加不同粒径以及不同颗粒数密度的散射颗粒，计算其光谱辐射强度变化．参考文献[19-20]中对于炉膛内颗粒参数模拟计算以及实验测量数据，工况5中固定颗粒数密度为2.0×108 m-3，颗粒粒径依次改变为10，50，100 μm；工况6中固定颗粒粒径为100 μm，颗粒数密度依次改变为1.0×107，5.0×107，1.0×108 m-3．图6展示了工况5中不同颗粒粒径下边界光谱辐射强度及相对误差计算结果．随着粒径从10 μm增大到100 μm，散射系数从约0.1 m-1增大到约6.3 m-1，在大部分高发射谱带，谱带辐射强度随着粒径增大而减小．DRESOR法结合SNB方法在不同颗粒粒径工况下的计算误差除少数辐射强度较小光谱位置外，误差均小于5%．10.13245/j.hust.220501.F006图6工况5不同颗粒粒径下辐射强度及相对误差图7为工况6的计算结果，随着颗粒数密度从1.0×107 m-3增大到1.0×108 m-3，谱带辐射强度随着颗粒数密度增大而略微有所减小．相较于不同的颗粒粒径，颗粒数密度对辐射强度结果的影响偏小．DRESOR法结合SNB模型在处理不同颗粒数密度散射颗粒时有着较好的计算精度，大部分谱带误差小于5%．10.13245/j.hust.220501.F007图7工况6不同颗粒数密度下辐射强度及相对误差4 结语为了高效计算含散射性颗粒和反射性边界的气体非灰介质中方向辐射强度，本研究将DRESOR法和SNB模型相结合，以一维平行平板系统为研究对象，以LBL模型计算结果作为基准解，分析了在等温均匀介质与非等温非均匀介质情况下计算结果精度．结果表明：DRESOR法结合SNB模型在绝大部分具有明显发射的光谱范围内误差小于8%，且强发射波段误差进一步减小；对加入散射性颗粒或反射性壁面的辐射系统，光谱辐射强度计算误差较黑体边界的纯气体介质辐射系统误差小，绝大部分具有明显发射谱带的误差小于5%．另外，讨论了不同散射颗粒粒径与颗粒数密度对边界方向辐射强度影响．颗粒散射系数随着颗粒粒径、颗粒数密度增大而增大，在高发射谱带，光谱辐射强度随之减弱．与颗粒数密度相比，颗粒粒径变化对光谱辐射强度结果的影响较大．