工業(yè)熱爐膛內(nèi)熱輻射場分析及反演重建

禹 健， 劉 鑫， 安永泉

(1. 山西大學 自動化與軟件學院， 山西 太原 030006； 2. 中北大學 信息與通信工程學院， 山西 太原 030051)

0 引 言

我國目前電力生產(chǎn)仍以燃煤發(fā)電為主. 提高燃煤發(fā)電機組的發(fā)電效率以及對爐膛燃燒火焰場有效監(jiān)視、 診斷和控制， 是降低供電煤耗的重要技術手段之一. 大型電站鍋爐爐膛內(nèi)的火焰溫度分布是燃料在經(jīng)過高溫化學反應、 流動以及傳熱傳質等過程后的綜合體現(xiàn). 對三維爐膛溫度場輻射分布機理的研究， 可探索燃燒現(xiàn)象的本質和規(guī)律， 幫助運行操作人員直觀觀察、 調(diào)整和優(yōu)化， 同時， 為燃燒設備的機械設計提供依據(jù).

工業(yè)熱爐膛數(shù)學模型的建立， 最先使用的是零維模型(Zero dimensional model)， 即“充分攪拌的爐膛”模型[1]. Salter和Costik引入“假想斷面”的概念[2]， 簡化了多段加熱爐熱輻射場問題的求解. 在工業(yè)爐膛溫度場重建與監(jiān)測的研究中， 聲學CT法是一個重要的分支， 多采用聲波在氣相介質中沿路徑傳播的原理[3-4]. 考慮到傳統(tǒng) CT 算法中忽略了火焰吸收的影響， 目前越來越多的學者從熱輻射傳遞過程的分析出發(fā)， 通過對爐內(nèi)流體介質處理， 從邊界輻射的角度來實現(xiàn)爐膛三維溫度場的重建[5-6]， 還有部分學者采用CFD數(shù)值模擬進行富氧燃燒的傳熱計算[7]， 基于火焰圖像處理的燃燒診斷開展研究[8-10].

本文研究借鑒“假想斷面”的概念， 引入網(wǎng)格定位模型， 按照吸收發(fā)射系數(shù)計算方法得到網(wǎng)格單元的輻射強度， 沿虛擬射線路徑積分光譜輻射分布得到爐膛熱輻射場的分布機理和解算方法， 并基于此模型的結果解析紅外圖像， 且從多幀圖像反卷積復原， 迭代得到熱輻射場的分布與重建.

1 物理模型

圖 1 為單爐膛、 Ⅱ型四角切圓類工業(yè)熱爐本體及火焰探頭布置示意圖. 關注光學探頭視界內(nèi)爐膛主柱體內(nèi)介質. 火焰探頭光學孔徑為 60 mm×60 mm 的正方形[10]. 燃燒爐膛的三維溫度分布區(qū)域為冷灰斗以上、 折焰角以下. 以熱爐幾何中心為所分析三維柱體介質的中心點， 尺寸為DX×DY×DZ(標高DZ∈[12 m,27.6 m]， 整個爐膛橫截面DX×DY=17 m×8.475 m)[11]， 其內(nèi)折射率、 溫度分布分別為n(x,y,z)，T(x,y,z).對于在三維空間任意連續(xù)分布的溫度T(x,y,z)， 將溫度場進行離散， 可直接求解輻射的彎曲路徑. 令介質與環(huán)境的界面是半透明鏡反射界面.

圖1 工業(yè)熱爐本體及火焰探頭布置示意圖Fig.1 Industrial furnace body and arrangement of flame probes

紅外CCD火焰探頭共16支， 分4層， 每支火焰探頭均垂直伸入爐膛， 主要監(jiān)測上、 下兩層燃燒器區(qū)域、 爐膛火焰中心區(qū)域和折焰角以下區(qū)域[12].

將三維爐膛內(nèi)熱流體介質離散， 沿爐寬X， 爐深Y， 爐高Z3個坐標方向， 將溫度場計算區(qū)域劃分為 10×10×18 的網(wǎng)格. 網(wǎng)格單元的立體尺寸為 1×0.80×0.80 m3. 爐膛空間共3 600個網(wǎng)格單元， 如圖2(a)所示. 沿爐體底法向直徑取截面， 如圖2(b).

圖2 爐膛介質Fig.2 Medium in the furnace

2 基于吸收發(fā)射系數(shù)的輻射計算模型

虛擬射線路徑是指火焰探頭探測口徑內(nèi)的輻射積分路徑， 為沿CFD網(wǎng)格計算時的光線蹤跡.

首先建立CFD網(wǎng)格光線蹤跡模型， 定位虛擬射線路徑途經(jīng)的網(wǎng)格(i,j,k)的空間位置.

2.1 CFD網(wǎng)格光線蹤跡模型

對于一維介質， 溫度分布符合已知函數(shù)規(guī)律時， 采用蒙特卡洛彎曲光線蹤跡法， 可得到光線軌跡的解析解.

2.1.1 蒙特卡洛彎曲光線蹤跡法

圖 3 所示為厚度為d的一維半透明灰體層介質中的光線蹤跡.邊界的壁面的折射率分別為n0和nL.吸收系數(shù)κα和散射系數(shù)κs為常數(shù)， 介質溫度T沿Z方向的變化函數(shù)為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發(fā)射并被其吸收的輻射能為

(4)

式中： ΔV為體元的體積； ΔA為面元的面積； 輻射分配系數(shù)Dij定義為從面元或體元發(fā)出的總輻射中被面元或體元吸收的份額；δ0=(n0-1)T0，n0和T0分別為標準大氣參考折射率和溫度， 1 200 K～1 400 K,δ0∈[10,50]. 對于工業(yè)熱爐流場高溫氣體介質， 如果考慮某個獨立網(wǎng)格， 溫度線性變化， 其他物性均勻， 網(wǎng)格與網(wǎng)格之間邊界的吸收忽略， 取Dij=0.

2.1.2 網(wǎng)格定位

特征射線在控制體P中的行程δs即輻射積分曲線軌跡S(s)， 通過前述建立的介質離散物理模型進行分段追跡.

對圖2(b)中的截面進行射線追跡， 設網(wǎng)格單元(正方網(wǎng)格)大小為d， 節(jié)點1-1處對應的溫度為T11， 節(jié)點i-j處溫度為Tij，Z為爐高方向， 如圖 4 所示.

設火焰探測器與爐膛中心軸夾角為θ0， 通過節(jié)點1-1后折射角為θ1， 第k次(k定義為光線Z向折變次數(shù))折射時， 所對應折射角為θk； 第k個折變點處對應的坐標為(X,Z).令ΔXk表示光線平移量(偏差量)， ∑ΔXK表示光線的總實際偏移量.

圖4 光線追跡Fig.4 Ray tracing

xi(z)為每個網(wǎng)格單元中的光線方程， 為z的函數(shù)， 當規(guī)定網(wǎng)格內(nèi)溫度函數(shù)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數(shù)據(jù)網(wǎng)格范圍， 計算無效.

2.2 基于分子原子輻射的高溫氣體輻射計算

工業(yè)熱爐內(nèi)輻射場發(fā)射吸收光譜計算的初始條件為：

1) 考慮N2， O2和CO23種主要成分；

2) 考慮束縛-束縛躍遷、 束縛-自由躍遷和自由-自由躍遷3類能級躍遷形式；

3) 電子能級、 振動能級和轉動能級密度分布符合Boltzmann分布， 以電子溫度Te， 振動溫度Tv和轉動平均溫度TR表示.

高溫下空氣組分將發(fā)生多種形式的量子輻射躍遷過程， 在三溫度模型[14]基礎上， 計算非平衡輻射吸收系數(shù)首先需確定各種組分輻射的波長.考慮每種輻射情況的高態(tài)和低態(tài)能級， 對于原子譜線輻射， 可根據(jù)其確定譜線直接得出； 對于分子波段輻射， 可通過轉動常數(shù)Be，αe和振動常數(shù)ωe，ωexe，ωeye，ωeze間接得出[15].所涉及客觀參數(shù)通過高溫氣體輻射數(shù)據(jù)HITRAN2008和HITEMP2010得到.其次進行強度因子的計算.對于原子譜線輻射， 即為愛因斯坦系數(shù)A； 對于分子譜線輻射， 考慮電子躍遷過程， 為∑|Re|2以及Franck-Condon因子； 對于自由-束縛和自由-自由連續(xù)光譜， 則由吸收截面給出[16]. 在前兩項工作基礎上， 得到隨波長變化的氣體總組分發(fā)射系數(shù)和吸收系數(shù).

2.2.1 高溫氣體原子輻射吸收系數(shù)

對于束縛-束縛躍遷，

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

(7)

式中：F(v)是譜線線形因子；vij為譜線中心波數(shù)(cm-1)；Bij,Bji分別為Einstein誘導發(fā)射和吸收系數(shù)；Ni,Nj為原子處在上、 下能態(tài)時的數(shù)密度(cm-3)；i為所求原子從束縛態(tài)到自由態(tài)可能發(fā)生的躍遷數(shù)目.

由愛因斯坦系數(shù)關系Aij/Bij=8πhv3，gijBij=gjBji, 可得

(8)

式中：gi,gj分別為上下態(tài)簡并度；Aij為Einstein自發(fā)發(fā)射系數(shù).

類氫近似時， 束縛-自由躍遷， 多電子原子的光電吸收截面可寫成如下形式

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

(9)

式中：gbf(v)為光電吸收的gaunt因子；σk(v)為經(jīng)典Kramer光電吸收橫截面.

(10)

式中：β為束縛能級的層數(shù)；v為頻率；me為電子質量；Z為原子電荷數(shù)目；e為電子電荷量.

求得光電吸收截面后， 易得吸收系數(shù)k(v).

(11)

類氫近似時， 自由-自由躍遷， 自由-自由吸收截面可寫成如下形式[16]

(12)

式中：Ne為電子數(shù)密度(cm-3)；gff(v)為韌致過程的gaunt因子.

逆韌幅吸收系數(shù)為

kv=N·σff(v)=

(13)

2.2.2 高溫氣體分子輻射吸收系數(shù)

分子振-轉譜線吸收系數(shù)為

(14)

通過以上模型， 得到爐膛溫度最高點1 400 K時， N2， O2和CO2的全光譜吸收發(fā)射系數(shù)分布和低能級態(tài)分布， 如圖 5～圖 7 所示.

圖5 二氮全光譜與低能態(tài)Fig.5 Full spectrum and low-energy state of N2

圖6 氧全光譜與低能態(tài)Fig.6 Full spectrum and low-energy state of O2

圖7 二氧化碳全光譜與低能態(tài)Fig.7 Full spectrum and low-energy state of CO2

基于網(wǎng)格定位和分子原子輻射的高溫氣體輻射計算模型， 可以得到 (1)各組分的高低態(tài)能級分布(圖 5～圖 7)； (2) 原子輻射吸收系數(shù)(式(6)～式(12))和分子輻射吸收系數(shù)(式(13)).

每一標注網(wǎng)格單元的輻射強度Ii,j,k(T(i,j,k),ρ(i,j,k))為

(15)

在火焰探頭探測口徑內(nèi)沿著每個截面的虛擬射線蹤跡積分累加，

(16)

可得出虛擬射線蹤跡的流場氣動熱輻射強度.

圖8 3種組分的吸收發(fā)射合成光譜Fig.8 Absorption emission synthesis spectra of three components

火焰探測器所得的近紅外圖像因熱場非均勻性產(chǎn)生了畸變效應. 通過復原算法， 反推引起畸變效應的熱場離散模型參數(shù).

3 序列圖像復原

將每個火焰探頭的CCD靶面離散為30×30=900 個象素單元， 刷新率不超過5 s.

對于空不變線性系統(tǒng)， 圖像退化過程為

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

(17)

式中:g為退化圖像;n為加性噪聲， 假設其服從方差為σ的零均值高斯分布.

MRF模型實質為能量函數(shù)最小的后驗估計[18]

(18)

式中：vc(f)由歸一化函數(shù)Z、 能量函數(shù)U(f)和先驗模型P(f)確定， 三者關系如式(19)所示.

(19)

式中：Ω為隨機場所有可能的結構組成的集合； 簇c表現(xiàn)了位置之間的關系或紋理的基本構成.前述熱輻射場提供了輻射分布分析模型P(f)， 如圖 9 所示.

圖9 工業(yè)熱爐溫度場輻射分布分析模型Fig.9 Analysis model of radiation distribution in temperaturefield of industrial furnace

正則化參數(shù)λ取式(19)為零的最小二乘解.

(20)

式中：C1和C2分別為[-1,1]T和[-1,1]的 Toeplitz 陣；H和V為水平線場和垂直線場的邏輯非.

參數(shù)r關系到構造不連續(xù)性的閾值的大小， 用信噪比信息約束， 根據(jù)所加降晰模板人工設定.

4 結果與討論

根據(jù)多視角序列圖像， 迭代反衍得到爐膛熱輻射場， 過程如圖 10 所示， 其中ΔP為兩幀圖像的矢量平均余弦相似度參數(shù)值.

圖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， 循環(huán)迭代， 直至ΔP滿足要求.

圖 11 為火焰探頭與爐膛中軸夾角為45°時， 中心截面的二維熱輻射場強度分布. 網(wǎng)格尺寸為1 m×0.85 m×0.847 m， 追蹤區(qū)域內(nèi)的網(wǎng)格個數(shù)為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單爐膛、 Ⅱ型布置四角切圓類工業(yè)熱爐， 爐膛中軸夾角45°火焰探頭和30°火焰探頭得到的原始圖像.

圖12 250 MW 45°火焰探頭原始幀與復原幀F(xiàn)ig.12 Original frame and restore frame at 250 MW 45°

圖13 250 MW 30°火焰探頭原始幀與復原幀F(xiàn)ig.13 Original frame and restore frame at 250 MW 30°

圖12(b)和圖13(b)分別為在網(wǎng)格尺寸為1 m×0.85 m×0.847 m， 取多截面， 追蹤三維區(qū)域內(nèi)的網(wǎng)格個數(shù)為10×10×16的初始熱輻射場條件下， 矢量平均余弦相似度閾值為0.01時， 迭代復原的圖像.

5 結 論

為研究大型電站鍋爐爐膛內(nèi)的火焰溫度場分布， 發(fā)展了基于網(wǎng)格定位和原子分子輻射計算的熱輻射強度的精確遞推計算模型， 其中采用虛擬射線追跡完成輻射積分路徑的選取， 模型結果作為MRF圖像復原模型的約束條件， 對火焰探頭所獲紅外圖像進行迭代復原， 得到爐膛熱輻射場的精確數(shù)值解. 可據(jù)此分析爐膛內(nèi)部三維溫度場輻射特性與調(diào)節(jié)燃煤量， 送風量， 引風量一、 二次風配比的關系， 為最優(yōu)工況設計提供可靠的理論依據(jù)和接近實際物理空間結構的約束條件.