工业热炉膛内热辐射场分析及反演重建

禹 健， 刘 鑫， 安永泉

(1. 山西大学 自动化与软件学院， 山西 太原 030006； 2. 中北大学 信息与通信工程学院， 山西 太原 030051)

0 引 言

我国目前电力生产仍以燃煤发电为主. 提高燃煤发电机组的发电效率以及对炉膛燃烧火焰场有效监视、 诊断和控制， 是降低供电煤耗的重要技术手段之一. 大型电站锅炉炉膛内的火焰温度分布是燃料在经过高温化学反应、 流动以及传热传质等过程后的综合体现. 对三维炉膛温度场辐射分布机理的研究， 可探索燃烧现象的本质和规律， 帮助运行操作人员直观观察、 调整和优化， 同时， 为燃烧设备的机械设计提供依据.

工业热炉膛数学模型的建立， 最先使用的是零维模型(Zero dimensional model)， 即“充分搅拌的炉膛”模型[1]. Salter和Costik引入“假想断面”的概念[2]， 简化了多段加热炉热辐射场问题的求解. 在工业炉膛温度场重建与监测的研究中， 声学CT法是一个重要的分支， 多采用声波在气相介质中沿路径传播的原理[3-4]. 考虑到传统 CT 算法中忽略了火焰吸收的影响， 目前越来越多的学者从热辐射传递过程的分析出发， 通过对炉内流体介质处理， 从边界辐射的角度来实现炉膛三维温度场的重建[5-6]， 还有部分学者采用CFD数值模拟进行富氧燃烧的传热计算[7]， 基于火焰图像处理的燃烧诊断开展研究[8-10].

本文研究借鉴“假想断面”的概念， 引入网格定位模型， 按照吸收发射系数计算方法得到网格单元的辐射强度， 沿虚拟射线路径积分光谱辐射分布得到炉膛热辐射场的分布机理和解算方法， 并基于此模型的结果解析红外图像， 且从多帧图像反卷积复原， 迭代得到热辐射场的分布与重建.

1 物理模型

图 1 为单炉膛、 Ⅱ型四角切圆类工业热炉本体及火焰探头布置示意图. 关注光学探头视界内炉膛主柱体内介质. 火焰探头光学孔径为 60 mm×60 mm 的正方形[10]. 燃烧炉膛的三维温度分布区域为冷灰斗以上、 折焰角以下. 以热炉几何中心为所分析三维柱体介质的中心点， 尺寸为DX×DY×DZ(标高DZ∈[12 m,27.6 m]， 整个炉膛横截面DX×DY=17 m×8.475 m)[11]， 其内折射率、 温度分布分别为n(x,y,z)，T(x,y,z).对于在三维空间任意连续分布的温度T(x,y,z)， 将温度场进行离散， 可直接求解辐射的弯曲路径. 令介质与环境的界面是半透明镜反射界面.

图1 工业热炉本体及火焰探头布置示意图Fig.1 Industrial furnace body and arrangement of flame probes

红外CCD火焰探头共16支， 分4层， 每支火焰探头均垂直伸入炉膛， 主要监测上、 下两层燃烧器区域、 炉膛火焰中心区域和折焰角以下区域[12].

将三维炉膛内热流体介质离散， 沿炉宽X， 炉深Y， 炉高Z3个坐标方向， 将温度场计算区域划分为 10×10×18 的网格. 网格单元的立体尺寸为 1×0.80×0.80 m3. 炉膛空间共3 600个网格单元， 如图2(a)所示. 沿炉体底法向直径取截面， 如图2(b).

图2 炉膛介质Fig.2 Medium in the furnace

2 基于吸收发射系数的辐射计算模型

虚拟射线路径是指火焰探头探测口径内的辐射积分路径， 为沿CFD网格计算时的光线踪迹.

首先建立CFD网格光线踪迹模型， 定位虚拟射线路径途经的网格(i,j,k)的空间位置.

2.1 CFD网格光线踪迹模型

对于一维介质， 温度分布符合已知函数规律时， 采用蒙特卡洛弯曲光线踪迹法， 可得到光线轨迹的解析解.

2.1.1 蒙特卡洛弯曲光线踪迹法

图 3 所示为厚度为d的一维半透明灰体层介质中的光线踪迹.边界的壁面的折射率分别为n0和nL.吸收系数κα和散射系数κs为常数， 介质温度T沿Z方向的变化函数为Ti(z)=(Ti+1-Ti)z+Ti.

(m0+T(z1))T(z2)cosθ1].

(1)

(m0+T(z1))T(z2)cosθ1].

(2)

(3)

由体元i和面元l发射并被其吸收的辐射能为

(4)

式中： ΔV为体元的体积； ΔA为面元的面积； 辐射分配系数Dij定义为从面元或体元发出的总辐射中被面元或体元吸收的份额；δ0=(n0-1)T0，n0和T0分别为标准大气参考折射率和温度， 1 200 K～1 400 K,δ0∈[10,50]. 对于工业热炉流场高温气体介质， 如果考虑某个独立网格， 温度线性变化， 其他物性均匀， 网格与网格之间边界的吸收忽略， 取Dij=0.

2.1.2 网格定位

特征射线在控制体P中的行程δs即辐射积分曲线轨迹S(s)， 通过前述建立的介质离散物理模型进行分段追迹.

对图2(b)中的截面进行射线追迹， 设网格单元(正方网格)大小为d， 节点1-1处对应的温度为T11， 节点i-j处温度为Tij，Z为炉高方向， 如图 4 所示.

设火焰探测器与炉膛中心轴夹角为θ0， 通过节点1-1后折射角为θ1， 第k次(k定义为光线Z向折变次数)折射时， 所对应折射角为θk； 第k个折变点处对应的坐标为(X,Z).令ΔXk表示光线平移量(偏差量)， ∑ΔXK表示光线的总实际偏移量.

图4 光线追迹Fig.4 Ray tracing

xi(z)为每个网格单元中的光线方程， 为z的函数， 当规定网格内温度函数Ti(z)=(Ti+1-Ti)z+Ti时， 得到光线偏移量xi(z)的解析表达式为

(5)

(6)

其边界条件为x(0)=0，x′(0)=tanθ0.

T(0)为界面1(z=0)处介质的温度.

若∑ΔXk=d1×d2(d1∈N)， 所追击光线超出给定的CFD数据网格范围， 计算无效.

2.2 基于分子原子辐射的高温气体辐射计算

工业热炉内辐射场发射吸收光谱计算的初始条件为：

1) 考虑N2， O2和CO23种主要成分；

2) 考虑束缚-束缚跃迁、 束缚-自由跃迁和自由-自由跃迁3类能级跃迁形式；

3) 电子能级、 振动能级和转动能级密度分布符合Boltzmann分布， 以电子温度Te， 振动温度Tv和转动平均温度TR表示.

高温下空气组分将发生多种形式的量子辐射跃迁过程， 在三温度模型[14]基础上， 计算非平衡辐射吸收系数首先需确定各种组分辐射的波长.考虑每种辐射情况的高态和低态能级， 对于原子谱线辐射， 可根据其确定谱线直接得出； 对于分子波段辐射， 可通过转动常数Be，αe和振动常数ωe，ωexe，ωeye，ωeze间接得出[15].所涉及客观参数通过高温气体辐射数据HITRAN2008和HITEMP2010得到.其次进行强度因子的计算.对于原子谱线辐射， 即为爱因斯坦系数A； 对于分子谱线辐射， 考虑电子跃迁过程， 为∑|Re|2以及Franck-Condon因子； 对于自由-束缚和自由-自由连续光谱， 则由吸收截面给出[16]. 在前两项工作基础上， 得到随波长变化的气体总组分发射系数和吸收系数.

2.2.1 高温气体原子辐射吸收系数

对于束缚-束缚跃迁，

k(v)=(NiBij-NjBij)hcvijF(v)，

(7)

式中：F(v)是谱线线形因子；vij为谱线中心波数(cm-1)；Bij,Bji分别为Einstein诱导发射和吸收系数；Ni,Nj为原子处在上、 下能态时的数密度(cm-3)；i为所求原子从束缚态到自由态可能发生的跃迁数目.

由爱因斯坦系数关系Aij/Bij=8πhv3，gijBij=gjBji, 可得

(8)

式中：gi,gj分别为上下态简并度；Aij为Einstein自发发射系数.

类氢近似时， 束缚-自由跃迁， 多电子原子的光电吸收截面可写成如下形式

σbf(v)=σk(v)·gbf(v)，

(9)

式中：gbf(v)为光电吸收的gaunt因子；σk(v)为经典Kramer光电吸收横截面.

(10)

式中：β为束缚能级的层数；v为频率；me为电子质量；Z为原子电荷数目；e为电子电荷量.

求得光电吸收截面后， 易得吸收系数k(v).

(11)

类氢近似时， 自由-自由跃迁， 自由-自由吸收截面可写成如下形式[16]

(12)

式中：Ne为电子数密度(cm-3)；gff(v)为韧致过程的gaunt因子.

逆韧幅吸收系数为

kv=N·σff(v)=

(13)

2.2.2 高温气体分子辐射吸收系数

分子振-转谱线吸收系数为

(14)

通过以上模型， 得到炉膛温度最高点1 400 K时， N2， O2和CO2的全光谱吸收发射系数分布和低能级态分布， 如图 5～图 7 所示.

图5 二氮全光谱与低能态Fig.5 Full spectrum and low-energy state of N2

图6 氧全光谱与低能态Fig.6 Full spectrum and low-energy state of O2

图7 二氧化碳全光谱与低能态Fig.7 Full spectrum and low-energy state of CO2

基于网格定位和分子原子辐射的高温气体辐射计算模型， 可以得到 (1)各组分的高低态能级分布(图 5～图 7)； (2) 原子辐射吸收系数(式(6)～式(12))和分子辐射吸收系数(式(13)).

每一标注网格单元的辐射强度Ii,j,k(T(i,j,k),ρ(i,j,k))为

(15)

在火焰探头探测口径内沿着每个截面的虚拟射线踪迹积分累加，

(16)

可得出虚拟射线踪迹的流场气动热辐射强度.

图8 3种组分的吸收发射合成光谱Fig.8 Absorption emission synthesis spectra of three components

火焰探测器所得的近红外图像因热场非均匀性产生了畸变效应. 通过复原算法， 反推引起畸变效应的热场离散模型参数.

3 序列图像复原

将每个火焰探头的CCD靶面离散为30×30=900 个象素单元， 刷新率不超过5 s.

对于空不变线性系统， 图像退化过程为

g(x,y)=f(x,y)·h(x,y)+n(x,y),

(17)

式中:g为退化图像;n为加性噪声， 假设其服从方差为σ的零均值高斯分布.

MRF模型实质为能量函数最小的后验估计[18]

(18)

式中：vc(f)由归一化函数Z、 能量函数U(f)和先验模型P(f)确定， 三者关系如式(19)所示.

(19)

式中：Ω为随机场所有可能的结构组成的集合； 簇c表现了位置之间的关系或纹理的基本构成.前述热辐射场提供了辐射分布分析模型P(f)， 如图 9 所示.

图9 工业热炉温度场辐射分布分析模型Fig.9 Analysis model of radiation distribution in temperaturefield of industrial furnace

正则化参数λ取式(19)为零的最小二乘解.

(20)

式中：C1和C2分别为[-1,1]T和[-1,1]的 Toeplitz 阵；H和V为水平线场和垂直线场的逻辑非.

参数r关系到构造不连续性的阈值的大小， 用信噪比信息约束， 根据所加降晰模板人工设定.

4 结果与讨论

根据多视角序列图像， 迭代反衍得到炉膛热辐射场， 过程如图 10 所示， 其中ΔP为两帧图像的矢量平均余弦相似度参数值.

图10 炉膛热辐射场反衍重建Fig.10 Reconstruction of furnace thermal radiation field

将热辐射计算模型计算结果P1(f)和同一位置火焰得到的红外图像g1作为复原模型WSMRF的初始约束和输入， 得到估计帧g2，P1(f) 和g2共同输入奇异值特征计算SVD(Singular Value Decomposition)模型， 若矢量平均余弦相似度大于阈值， 则进行下一次迭代， 其初始约束为P2(f)=P1(f)*g2， 输入为g2， 循环迭代， 直至ΔP满足要求.

图 11 为火焰探头与炉膛中轴夹角为45°时， 中心截面的二维热辐射场强度分布. 网格尺寸为1 m×0.85 m×0.847 m， 追踪区域内的网格个数为10×16.

图11 热辐射场二维辐射强度(探头与炉膛中轴夹角45°)Fig.11 Two dimensional radiation intensity of the thermalradiation field(The angle between the probe and center axis is 45°)

图12(a)和图13(a)分别为250 MW单炉膛、 Ⅱ型布置四角切圆类工业热炉， 炉膛中轴夹角45°火焰探头和30°火焰探头得到的原始图像.

图12 250 MW 45°火焰探头原始帧与复原帧Fig.12 Original frame and restore frame at 250 MW 45°

图13 250 MW 30°火焰探头原始帧与复原帧Fig.13 Original frame and restore frame at 250 MW 30°

图12(b)和图13(b)分别为在网格尺寸为1 m×0.85 m×0.847 m， 取多截面， 追踪三维区域内的网格个数为10×10×16的初始热辐射场条件下， 矢量平均余弦相似度阈值为0.01时， 迭代复原的图像.

5 结 论

为研究大型电站锅炉炉膛内的火焰温度场分布， 发展了基于网格定位和原子分子辐射计算的热辐射强度的精确递推计算模型， 其中采用虚拟射线追迹完成辐射积分路径的选取， 模型结果作为MRF图像复原模型的约束条件， 对火焰探头所获红外图像进行迭代复原， 得到炉膛热辐射场的精确数值解. 可据此分析炉膛内部三维温度场辐射特性与调节燃煤量， 送风量， 引风量一、 二次风配比的关系， 为最优工况设计提供可靠的理论依据和接近实际物理空间结构的约束条件.