工業熱爐膛內熱輻射場分析及反演重建

禹 健， 劉 鑫， 安永泉

(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圖像復原模型的約束條件， 對火焰探頭所獲紅外圖像進行迭代復原， 得到爐膛熱輻射場的精確數值解. 可據此分析爐膛內部三維溫度場輻射特性與調節燃煤量， 送風量， 引風量一、 二次風配比的關系， 為最優工況設計提供可靠的理論依據和接近實際物理空間結構的約束條件.