基于機(jī)器學(xué)習(xí)和數(shù)值模擬的鍋爐水冷壁熱流密度分布預(yù)測模型

董凌霄，梁永，楊家輝，靳曉靈，杜勇博，鄧?yán)?，車得?/p>

(西安交通大學(xué) 能源與動力工程學(xué)院，陜西 西安，710049)

煤炭在我國一次能源供應(yīng)體系中具有重要地位，燃煤發(fā)電也在我國的電源結(jié)構(gòu)中占據(jù)著主導(dǎo)地位[1]?！疤歼_(dá)峰、碳中和”目標(biāo)的提出加速了我國經(jīng)濟(jì)和能源結(jié)構(gòu)的轉(zhuǎn)型[2-3]。盡管可再生能源持續(xù)強(qiáng)勁增長，風(fēng)能與太陽能裝機(jī)容量持續(xù)增加，但燃煤機(jī)組在以新能源發(fā)展為主體的新型電力系統(tǒng)下的作用非常重要。為適應(yīng)風(fēng)力和太陽能發(fā)電等新型發(fā)電系統(tǒng)的波動性，火力發(fā)電在電網(wǎng)調(diào)峰中的作用將更加明顯，這對燃煤機(jī)組變負(fù)荷運(yùn)行能力提出了更高要求。爐內(nèi)工作過程由燃燒和爐內(nèi)傳熱組成，爐內(nèi)傳熱不僅直接影響爐膛溫度，還決定著離開爐膛進(jìn)入對流受熱面的能量[4]。為保障鍋爐在變負(fù)荷運(yùn)行下的安全性，確保受熱面尤其是水冷壁處在合理的耐溫極限內(nèi)是關(guān)鍵。

水冷壁熱流密度與鍋爐輸入熱量和爐內(nèi)燃燒情況密切相關(guān)。水冷壁局部熱流密度可用來計算工質(zhì)參數(shù)沿程變化，校驗(yàn)水循環(huán)動力和水冷壁壁溫[5-7]，其分布規(guī)律也是研究爐膛內(nèi)換熱的重要指標(biāo)。水冷壁熱流密度分布可通過實(shí)驗(yàn)測量[8]、經(jīng)驗(yàn)總結(jié)[9]和數(shù)值模擬[10-11]等途徑獲取。馬達(dá)夫等[12]使用熱電偶測量了某600 MW 超臨界W 火焰鍋爐水冷壁向火面和背火面的壁面溫度，并使用背火面溫度代替工質(zhì)溫度，基于傳熱學(xué)原理計算了水冷壁熱流密度分布。爐內(nèi)的惡劣環(huán)境會導(dǎo)致測點(diǎn)布置和儀表安裝存在困難，并且溫度測量的精度和可靠性也難以保證，這限制了實(shí)驗(yàn)測量方法的應(yīng)用。一些水動力計算標(biāo)準(zhǔn)給出了爐內(nèi)熱流密度的經(jīng)驗(yàn)分布，然而適用范圍有限。ZHANG等[13-15]利用數(shù)值模擬手段，分別針對600 MW超臨界拱形燃燒鍋爐、600 MW四角切圓煤粉鍋爐和超臨界二氧化碳鍋爐獲得了較全面的水冷壁熱流密度分布。與實(shí)驗(yàn)測量和經(jīng)驗(yàn)總結(jié)相比，數(shù)值模擬具有技術(shù)成熟、成本低和數(shù)據(jù)完備等優(yōu)點(diǎn)，但模擬過程復(fù)雜、耗時長，該缺點(diǎn)同樣限制了該方法的直接應(yīng)用。值得注意的是，機(jī)器學(xué)習(xí)算法在燃煤電廠建模中得到了廣泛應(yīng)用。以往的研究主要將機(jī)器學(xué)習(xí)應(yīng)用在NOx排放預(yù)測[16]、水冷壁積灰結(jié)渣在線監(jiān)測[17-18]、鍋爐效率優(yōu)化[19]等方面，預(yù)測水冷壁熱流密度分布的研究尚未見報道。將機(jī)器學(xué)習(xí)和數(shù)值模擬方法相結(jié)合，為快速、準(zhǔn)確預(yù)測水冷壁熱流密度分布提供了新思路。

本文首先通過數(shù)值模擬計算某660 MW超臨界直流鍋爐在48 種工況下的水冷壁熱流密度分布，為機(jī)器學(xué)習(xí)提供原始樣本數(shù)據(jù)；其次，通過皮爾遜相關(guān)分析和緊鄰成分分析確定特征子集，利用貝葉斯優(yōu)化結(jié)合交叉驗(yàn)證的方法優(yōu)化支持向量機(jī)、決策樹和提升樹算法模型；最后，基于優(yōu)選的算法模型，建立最終的水冷壁熱流密度分布預(yù)測模型。

1 鍋爐水冷壁熱流密度模擬計算

1.1 鍋爐結(jié)構(gòu)

圖1所示為某600 MW 超臨界直流鍋爐結(jié)構(gòu)及燃燒器布置示意圖。該爐膛寬度和深度分別為18 816 mm 和17 696 mm，水冷壁下集箱標(biāo)高為8 300 mm，爐頂管中心標(biāo)高為72 350 mm，爐底冷灰斗角度為55o。在爐膛冷灰斗進(jìn)口到爐膛標(biāo)高49 720 mm 處采用傾角為13.95°的螺旋管圈布置，螺旋管圈上部采用垂直管圈均勻布置在鍋爐四面墻上，螺旋式水冷壁和垂直上升式水冷壁之間采用混合集箱進(jìn)行連接轉(zhuǎn)換。沿著煙氣流動方向依次布置著6片分隔屏過熱器、20片后屏過熱器、末級再熱器和末級過熱器。噴嘴由煤粉噴嘴(A-F)、二次風(fēng)噴嘴(UFA、AUX 和CCOFA)和燃盡風(fēng)噴嘴(SOFA)組成，以四角切圓方式組織煤粉的燃燒。鍋爐變負(fù)荷運(yùn)行時，給煤量、燃燒器投運(yùn)層數(shù)和過量空氣系數(shù)等運(yùn)行參數(shù)都隨之變化，數(shù)值模擬中采用的運(yùn)行參數(shù)如表1所示。

表1 不同負(fù)荷下的運(yùn)行參數(shù)Table 1 Operating parameters under different loads

圖1 某600 MW鍋爐結(jié)構(gòu)及燃燒器布置示意圖Fig.1 Boiler structure and burners arrangement of 600 MW

1.2 水冷壁熱流密度分布計算

1.2.1 網(wǎng)格劃分與邊界條件

采用ICEM中的六面體結(jié)構(gòu)化網(wǎng)格對該直流鍋爐進(jìn)行網(wǎng)格劃分。計算區(qū)域由冷灰斗底部至鍋爐后煙井。在冷灰斗區(qū)域，由于爐內(nèi)參數(shù)變化不大，網(wǎng)格設(shè)置較稀疏。在燃燒器區(qū)域，由于燃燒反應(yīng)劇烈、溫度等物理量變化較快，對該區(qū)域進(jìn)行網(wǎng)格加密。采用標(biāo)準(zhǔn)κ-ε模型對湍流流動進(jìn)行計算[20]；采用雙步競爭反應(yīng)模型描述揮發(fā)分[21]，焦炭燃燒采用擴(kuò)散/動力模型[22]進(jìn)行分析；煤粉顆粒在爐內(nèi)的燃燒和流動過程采用隨機(jī)軌道模型進(jìn)行分析[23]；采用離散坐標(biāo)輻射模型計算輻射換熱過程。在本研究中，燃燒器及燃盡風(fēng)入口均為質(zhì)量流量入口，一次風(fēng)和二次風(fēng)入口質(zhì)量流量分別為128.86 kg/s 和354.38 kg/s，一次風(fēng)和二次風(fēng)的溫度分別為353 K 和612 K。在模擬過程中，過量空氣系數(shù)為1.20，爐膛壁面條件選用定壁溫邊界條件，壁面溫度取水冷壁入口工質(zhì)溫度和出口工質(zhì)溫度的平均值，即670 K。該600 MW 超臨界直流鍋爐的設(shè)計煤種為神府東勝煤，本文共模擬了3個實(shí)驗(yàn)煤種用于后續(xù)特征選擇及預(yù)測模型，詳細(xì)煤質(zhì)參數(shù)如表2所示。

表2 實(shí)驗(yàn)煤種煤質(zhì)參數(shù)Table 2 Fuel properties of bituminous coal

1.2.2 模型驗(yàn)證

本文選用了網(wǎng)格數(shù)為1 337 476、1 752 193、2 067 250 和2 478 099 個這4 種網(wǎng)格進(jìn)行網(wǎng)格無關(guān)性檢驗(yàn)。由4種網(wǎng)格計算得到的爐膛截面平均溫度沿爐高方向變化趨勢基本一致，當(dāng)網(wǎng)格數(shù)量從1 337 476 個增加到2 067 250 個時，計算結(jié)果存在較大幅度的變化，但當(dāng)網(wǎng)格數(shù)量從2 067 250 個增加到2 478 099 個時，計算結(jié)果基本不變，這說明當(dāng)網(wǎng)格數(shù)量為2 067 250 個時，可獲得網(wǎng)格獨(dú)立的解?？紤]到計算時間與計算精度，最終選擇2 067 250個的網(wǎng)格進(jìn)行計算。

該鍋爐的蒸汽側(cè)溫度熱力計算結(jié)果與實(shí)際監(jiān)測結(jié)果的相對誤差很小，說明計算得到的受熱面出口煙氣溫度和吸熱量與運(yùn)行結(jié)果接近，可作為驗(yàn)證數(shù)值模擬可靠性的參考。屏底煙溫和末過(熱段)進(jìn)口煙溫的數(shù)值模擬與熱力計算結(jié)果對比見表3。從表3可以看出，在100%的BMCR負(fù)荷下，計算結(jié)果偏差最大為2.92%；除此之外，水冷壁吸熱量在100%的BMCR 負(fù)荷工況下模擬值為659.33 MW，熱力計算結(jié)果為609.86 MW，相對誤差為8.11%。從表3還可見：低負(fù)荷的末過(熱段)進(jìn)口煙溫最大相對誤差在8%左右，高負(fù)荷以及低負(fù)荷下屏底煙溫相對誤差均在5%以內(nèi)。根據(jù)以往的研究，認(rèn)為該計算模型是可靠的[24-27]。

表3 屏底煙溫和末過(熱段)進(jìn)口煙溫的數(shù)值模擬與熱力計算結(jié)果對比Table 3 Comparison between numerical calculation and thermal calculation for smoke temperature at screen and final superheater inletK

1.3 樣本建模與數(shù)據(jù)選取

在“雙碳”目標(biāo)和以新能源發(fā)展為主體的新型電力系統(tǒng)下，由于新能源的波動性，燃煤機(jī)組需要頻繁參與電網(wǎng)負(fù)荷深度調(diào)峰，鍋爐變負(fù)荷運(yùn)行已成為常態(tài)。當(dāng)鍋爐負(fù)荷變化時，爐內(nèi)溫度發(fā)生變化，進(jìn)而影響水冷壁熱流密度的分布。不同煤種之間灰分、水分、揮發(fā)分等煤質(zhì)參數(shù)差別較大，同一鍋爐燃用不同煤種時燃料燃燒特性也會不同。增加煤中水分和灰分含量或者降低揮發(fā)分和熱值均會導(dǎo)致煤粉的著火時間推遲，火焰拖長且火焰中心上移，從而引起水冷壁熱流密度變化。鍋爐采用的空氣分級燃燒等低NOx技術(shù)也會影響水冷壁的熱流密度分布。

綜上所述，影響水冷壁熱流密度的因素主要包括鍋爐負(fù)荷、燃料特性及運(yùn)行方式3個方面。運(yùn)行方式參數(shù)包含一次風(fēng)率、燃燒器擺動角度、SOFA擺動角度、主燃區(qū)過量空氣系數(shù)、風(fēng)溫、風(fēng)速等。通過數(shù)值模擬，為機(jī)器學(xué)習(xí)預(yù)測模型提供原始的樣本數(shù)據(jù)集。設(shè)計數(shù)值模擬的工況，共計48組工況，部分?jǐn)?shù)值模擬工況見表4。

表4 部分?jǐn)?shù)值模擬工況Table 4 Part of numerical simulation conditions

2 機(jī)器學(xué)習(xí)預(yù)測模型

2.1 特征選擇

機(jī)器學(xué)習(xí)模型預(yù)測性能的上限由數(shù)據(jù)質(zhì)量決定。在數(shù)據(jù)采集完成后，詳細(xì)分析數(shù)據(jù)，進(jìn)一步選擇與提取特征，從而找到對響應(yīng)變量影響最大的特征。在采集完成由48 組工況模擬結(jié)果組成的數(shù)據(jù)集后，需對數(shù)據(jù)進(jìn)行相關(guān)性及近鄰成分分析，最終確定預(yù)測模型的輸入變量。

皮爾遜相關(guān)系數(shù)用于衡量2個變量間線性相關(guān)性，其值介于-1～1 之間。當(dāng)2 個變量之間的皮爾遜相關(guān)系數(shù)大于0.4時，認(rèn)為變量之間具有中度或強(qiáng)相關(guān)性[28]。皮爾遜相關(guān)系數(shù)r的計算公式如下：

式中：n為樣本數(shù)量；Xi為變量X的第i個觀測值；-X為變量X的樣本平均數(shù)；Yi為變量Y的第i個觀測值；-Y為變量Y的樣本平均數(shù)。

皮爾遜相關(guān)分析結(jié)果如表5所示。從表5 可見：鍋爐負(fù)荷、給煤量、總風(fēng)量、一次風(fēng)量、二次風(fēng)量、二次風(fēng)溫、二次風(fēng)速以及前墻水冷壁網(wǎng)格單元幾何中心縱坐標(biāo)y等參數(shù)與該鍋爐前墻的壁面熱流密度分布具有很強(qiáng)的相關(guān)性，但煤質(zhì)參數(shù)、SOFA風(fēng)量、燃燒器擺動角度、SOFA風(fēng)擺動角度、一次風(fēng)速、SOFA風(fēng)速和前墻水冷壁網(wǎng)格單元幾何中心橫坐標(biāo)z等參數(shù)與鍋爐前墻壁面的熱流密度的相關(guān)性不大，這與前人的研究結(jié)果[29-30]存在一定的矛盾。由于鍋爐燃燒系統(tǒng)具有多變量、非線性以及強(qiáng)耦合等特點(diǎn)，需要進(jìn)一步選擇特征。

表5 各個特征的皮爾遜分析結(jié)果Table 5 Results of Pearson's analysis for each sign

在皮爾遜相關(guān)分析的基礎(chǔ)上，利用近鄰成分分析算法進(jìn)一步選擇特征。近鄰成分分析方法可以最大程度地提高回歸和分類算法的預(yù)測準(zhǔn)確性。近鄰成分分析特征選擇的結(jié)果如圖2所示(其中，特征編號見表5)，編號為1(C 質(zhì)量分?jǐn)?shù))、12(給煤量)、13(總風(fēng)量)、15(二次風(fēng)量)、17(燃燒器擺動角度)、18(SOFA 風(fēng)擺動角度)、23(水冷壁網(wǎng)格幾何中心坐標(biāo)y)、24(水冷壁網(wǎng)格幾何中心坐標(biāo)z)可作為特征選擇的結(jié)果。

圖2 近鄰成分分析的結(jié)果Fig.2 Results of neighborhood component analysis

全面考慮煤質(zhì)參數(shù)、機(jī)組負(fù)荷、配風(fēng)情況以及網(wǎng)格單元幾何中心坐標(biāo)等關(guān)鍵因素，結(jié)合模型復(fù)雜度和精準(zhǔn)度，取2種方法得到的特征子集的并集為預(yù)測模型的輸入變量。由于鍋爐負(fù)荷和給煤量這2個特征之間并不相互獨(dú)立，刪除鍋爐負(fù)荷這個冗雜的特征，最終確定共11 個特征作為水冷壁熱流密度預(yù)測模型的輸入變量，包括C、給煤量、總風(fēng)量、一次風(fēng)量、二次風(fēng)量、二次風(fēng)溫、燃燒器擺動角度、SOFA風(fēng)擺動角度、二次風(fēng)速以及水冷壁網(wǎng)格單元幾何中心坐標(biāo)y和z。

2.2 3種算法模型的建立和評估

利用MATLAB 自帶的randperm 函數(shù)，生成偽隨機(jī)數(shù)，并劃分?jǐn)?shù)據(jù)集，rng 設(shè)置為42。在9 600個樣本中，選擇85%的樣本作為訓(xùn)練集，剩余的15%樣本作為測試集。在使用算法構(gòu)建模型之前，需要對數(shù)據(jù)進(jìn)行歸一化預(yù)處理，使樣本每個特征值都映射到[-1，1]之間，不僅可加快計算時的預(yù)處理速度，還可以防止屬性值數(shù)量級較小的特征被淹沒，造成原始數(shù)據(jù)信息丟失。

式中：xscale、xmax、xmin和分別為屬性x歸一化后的數(shù)值、最大值、最小值和均值。

采用貝葉斯優(yōu)化算法并結(jié)合五折交叉驗(yàn)證，對支持向量機(jī)、決策樹、提升樹進(jìn)行超參數(shù)優(yōu)化。在優(yōu)化過程中，三者超參數(shù)選擇不同，迭代次數(shù)均設(shè)置為30 次。支持向量機(jī)優(yōu)化時，超參數(shù)為“框約束”和“核尺度”，這2 個參數(shù)分別以一定的函數(shù)關(guān)系對應(yīng)支持向量機(jī)算法的懲罰系數(shù)C和核系數(shù)g，搜索范圍均為0.001～1 000，優(yōu)化結(jié)果如下：框約束為19.53，核尺度為0.80。決策樹模型優(yōu)化時超參數(shù)為“最小葉”，超參數(shù)搜索范圍為1～4 080，最終最小葉優(yōu)化結(jié)果為1。提升樹優(yōu)化時超參數(shù)選擇最小葉、學(xué)習(xí)器數(shù)量和學(xué)習(xí)率，搜索范圍分別為1～4 080，10～500和0.001～1，優(yōu)化結(jié)果如下：最小葉為113，學(xué)習(xí)器數(shù)量為448個，學(xué)習(xí)率為0.31。

圖3 和圖4所示分別為3 種算法模型在訓(xùn)練集和測試集上的預(yù)測效果。從圖3 和4 可以看出：3種方法模型在訓(xùn)練集上的預(yù)測值與模擬值相對誤差較小，測試集上的預(yù)測值與目標(biāo)值也非常接近，與數(shù)值模擬的結(jié)果吻合程度較高，且沒有出現(xiàn)過擬合現(xiàn)象。3 種模型的預(yù)測性能(均方根誤差、均方誤差和決定系數(shù))如表6所示，支持向量機(jī)模型、決策樹模型和提升樹模型的R2依次為0.983 6、0.970 7和0.985 0，均比較接近1。因此，利用支持向量機(jī)、決策樹和提升樹建立的3種預(yù)測模型均具有較高的預(yù)測精度和較強(qiáng)的泛化能力，3種算法建立的水冷壁熱流密度預(yù)測模型是可行和可靠的。

表6 不同學(xué)習(xí)算法的預(yù)測性能Table 6 Predictive performance of different learning algorithms

圖3 不同算法預(yù)測模型在訓(xùn)練集上的預(yù)測效果Fig.3 Prediction of different algorithm prediction models on training set

圖4 不同算法預(yù)測模型在測試集上的預(yù)測效果Fig.4 Prediction of different algorithm prediction models on test set

從預(yù)測性能來看，提升樹算法建立的模型具有更高的預(yù)測準(zhǔn)確性，但與支持向量機(jī)算法構(gòu)建的模型相比并沒有明顯差距。支持向量機(jī)是基于小樣本統(tǒng)計學(xué)習(xí)理論建立的，具有不受限于樣本總量的特點(diǎn)，被認(rèn)為是針對小樣本條件下分類和回歸的最佳學(xué)習(xí)方法。鍋爐水冷壁熱流密度的模擬計算存在計算量大、耗時多等問題，難以獲得較大數(shù)量的樣本數(shù)據(jù)。綜合考慮樣本容量和預(yù)測性能，可采用支持向量機(jī)預(yù)測鍋爐水冷壁熱流密度。

2.3 鍋爐水冷壁熱流密度分布預(yù)測模型的建立和評估

基于以上研究，利用LIBSVM 工具箱編寫鍋爐水冷壁熱流密度預(yù)測軟件，以方便用戶利用支持向量機(jī)建立鍋爐水冷壁預(yù)測模型。預(yù)測模型采用MATLAB R2020b 中的App Designer 作為開發(fā)工具，運(yùn)用面向?qū)ο蟮某绦蛟O(shè)計技術(shù)編制。在設(shè)計技術(shù)編制過程中，還參考并引用了Faruto 版本的LIBSVM加強(qiáng)工具箱的部分輔助函數(shù)。歸一化范圍選擇[-1，1]，超參數(shù)優(yōu)化方法選擇貝葉斯優(yōu)化，懲罰系數(shù)C和核系數(shù)g的尋優(yōu)范圍均設(shè)置為[0，100]，交叉驗(yàn)證的折數(shù)為20，得到的最佳參數(shù)為C=99.21，g=0.000 49。采用最佳參數(shù)和訓(xùn)練集對模型進(jìn)行訓(xùn)練，得到鍋爐前墻壁面熱流密度峰值預(yù)測模型。模型在訓(xùn)練集上的擬合效果如圖5所示，其中，熱流密度正值表示煙氣將熱量傳遞給水冷壁。

圖5 預(yù)測模型在訓(xùn)練集上的鍋爐水冷壁熱流密度訓(xùn)練結(jié)果Fig.5 Results of predictive modeling on training set of heat flow density distribution in boiler water-cooled walls

圖6所示為預(yù)測模型在訓(xùn)練集上的鍋爐水冷壁熱流密度訓(xùn)練結(jié)果。模型在訓(xùn)練集上的均方誤差為0.039 2，R2為0.849 6，在測試集上的均方誤差為0.007 7，R2為0.979 7。預(yù)測軟件得到的結(jié)果具有較高的預(yù)測精度和較強(qiáng)的泛化能力，可以準(zhǔn)確地預(yù)測鍋爐的水冷壁熱流密度。

圖6 預(yù)測模型在測試集上的鍋爐水冷壁熱流密度訓(xùn)練結(jié)果Fig.6 Results of predictive modeling on test set of heat flow density distribution in boiler water-cooled walls

3 結(jié)論

1)通過理論分析得到水冷壁熱流密度的各個影響因素，設(shè)定樣本建模數(shù)據(jù)的工況。在此基礎(chǔ)上，通過皮爾遜相關(guān)性和近鄰成分分析，選擇C、給煤量、總風(fēng)量、一次風(fēng)量、二次風(fēng)量、二次風(fēng)溫、燃燒器擺動角度、SOFA風(fēng)擺動角度、二次風(fēng)速以及水冷壁網(wǎng)格單元幾何中心坐標(biāo)y和z作為水冷壁熱流密度分布預(yù)測模型的輸入變量。

2)利用貝葉斯優(yōu)化結(jié)合交叉驗(yàn)證的方法對支持向量機(jī)、決策樹和提升樹3 種算法進(jìn)行超參數(shù)優(yōu)化，三者所構(gòu)建的水冷壁熱流密度分布預(yù)測模型在測試集上的均方根誤差分別為1.711 2、2.276 6 和1.622 2，R2分別為0.983 6、0.970 7 和0.985 0，均具有較高的預(yù)測準(zhǔn)確性。經(jīng)綜合考慮，選擇支持向量機(jī)算法構(gòu)建水冷壁熱流密度預(yù)測模型。

3)利用優(yōu)化后的支持向量機(jī)算法建立的水冷壁熱流密度分布預(yù)測模型，模型具有較高的預(yù)測精度和較強(qiáng)的泛化能力，可以對鍋爐的水冷壁熱流密度進(jìn)行準(zhǔn)確預(yù)測。