禹 健， 劉 鑫， 安永泉

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

0 引 言

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

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

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

1 物理模型

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

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

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

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

圖2 爐膛介質(zhì)Fig.2 Medium in the furnace

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

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

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

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

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

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

圖 3 所示為厚度為d的一維半透明灰體層介質(zhì)中的光線蹤跡.邊界的壁面的折射率分別為n0和nL.吸收系數(shù)κα和散射系數(shù)κs為常數(shù)， 介質(zhì)溫度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分別為標(biāo)準(zhǔn)大氣參考折射率和溫度， 1 200 K～1 400 K,δ0∈[10,50]. 對(duì)于工業(yè)熱爐流場(chǎng)高溫氣體介質(zhì)， 如果考慮某個(gè)獨(dú)立網(wǎng)格， 溫度線性變化， 其他物性均勻， 網(wǎng)格與網(wǎng)格之間邊界的吸收忽略， 取Dij=0.

2.1.2 網(wǎng)格定位

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

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

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

圖4 光線追跡Fig.4 Ray tracing

xi(z)為每個(gè)網(wǎng)格單元中的光線方程， 為z的函數(shù)， 當(dāng)規(guī)定網(wǎng)格內(nèi)溫度函數(shù)Ti(z)=(Ti+1-Ti)z+Ti時(shí)， 得到光線偏移量xi(z)的解析表達(dá)式為

(5)

(6)

其邊界條件為x(0)=0，x′(0)=tanθ0.

T(0)為界面1(z=0)處介質(zhì)的溫度.

若∑ΔXk=d1×d2(d1∈N)， 所追擊光線超出給定的CFD數(shù)據(jù)網(wǎng)格范圍， 計(jì)算無效.

2.2 基于分子原子輻射的高溫氣體輻射計(jì)算

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

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

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

3) 電子能級(jí)、 振動(dòng)能級(jí)和轉(zhuǎn)動(dòng)能級(jí)密度分布符合Boltzmann分布， 以電子溫度Te， 振動(dòng)溫度Tv和轉(zhuǎn)動(dòng)平均溫度TR表示.

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

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

對(duì)于束縛-束縛躍遷，

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

(7)

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

由愛因斯坦系數(shù)關(guān)系A(chǔ)ij/Bij=8πhv3，gijBij=gjBji, 可得

(8)

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

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

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

(9)

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

(10)

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

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

(11)

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

(12)

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

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

kv=N·σff(v)=

(13)

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

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

(14)

通過以上模型， 得到爐膛溫度最高點(diǎn)1 400 K時(shí)， N2， O2和CO2的全光譜吸收發(fā)射系數(shù)分布和低能級(jí)態(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)格定位和分子原子輻射的高溫氣體輻射計(jì)算模型， 可以得到 (1)各組分的高低態(tài)能級(jí)分布(圖 5～圖 7)； (2) 原子輻射吸收系數(shù)(式(6)～式(12))和分子輻射吸收系數(shù)(式(13)).

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

(15)

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

(16)

可得出虛擬射線蹤跡的流場(chǎng)氣動(dòng)熱輻射強(qiáng)度.

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

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

3 序列圖像復(fù)原

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

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

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

(17)

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

MRF模型實(shí)質(zhì)為能量函數(shù)最小的后驗(yàn)估計(jì)[18]

(18)

式中：vc(f)由歸一化函數(shù)Z、 能量函數(shù)U(f)和先驗(yàn)?zāi)Ｐ蚉(f)確定， 三者關(guān)系如式(19)所示.

(19)

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

圖9 工業(yè)熱爐溫度場(chǎng)輻射分布分析模型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為水平線場(chǎng)和垂直線場(chǎng)的邏輯非.

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

4 結(jié)果與討論

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

圖10 爐膛熱輻射場(chǎng)反衍重建Fig.10 Reconstruction of furnace thermal radiation field

將熱輻射計(jì)算模型計(jì)算結(jié)果P1(f)和同一位置火焰得到的紅外圖像g1作為復(fù)原模型WSMRF的初始約束和輸入， 得到估計(jì)幀g2，P1(f) 和g2共同輸入奇異值特征計(jì)算SVD(Singular Value Decomposition)模型， 若矢量平均余弦相似度大于閾值， 則進(jìn)行下一次迭代， 其初始約束為P2(f)=P1(f)*g2， 輸入為g2， 循環(huán)迭代， 直至ΔP滿足要求.

圖 11 為火焰探頭與爐膛中軸夾角為45°時(shí)， 中心截面的二維熱輻射場(chǎng)強(qiáng)度分布. 網(wǎng)格尺寸為1 m×0.85 m×0.847 m， 追蹤區(qū)域內(nèi)的網(wǎng)格個(gè)數(shù)為10×16.

圖11 熱輻射場(chǎng)二維輻射強(qiáng)度(探頭與爐膛中軸夾角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ù)原幀F(xiàn)ig.12 Original frame and restore frame at 250 MW 45°

圖13 250 MW 30°火焰探頭原始幀與復(fù)原幀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)格個(gè)數(shù)為10×10×16的初始熱輻射場(chǎng)條件下， 矢量平均余弦相似度閾值為0.01時(shí)， 迭代復(fù)原的圖像.

5 結(jié) 論

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