CFD-BP 神經(jīng)網(wǎng)絡(luò)耦合算法對焚燒爐的一次風(fēng)管道配風(fēng)的模擬研究

周洪權(quán)，陳德珍，尹麗潔，張 智，戴小云，陸子葉

（1.同濟(jì)大學(xué)熱能與環(huán)境工程研究所，上海 200092；2.上海環(huán)境衛(wèi)生工程設(shè)計院有限公司，上海，200232）

1 引言

隨著國民經(jīng)濟(jì)高速發(fā)展，城市生活垃圾產(chǎn)量逐年增長，2021 年我國城市生活垃圾清運(yùn)量為2.49 億t，無害化處理率為99.9%，其中焚燒無害化處理量為1.80億t，占無害化處理量的72.4%［1-2］。“十四五”規(guī)劃提出，到2025年底，全國城鎮(zhèn)生活垃圾焚燒處理能力達(dá)到80萬t/d左右，遠(yuǎn)大于目前中國城鎮(zhèn)生活垃圾無害化處理能力［3］。隨著垃圾分類收運(yùn)、分類處理的推行，焚燒爐爐前垃圾的水分降低、熱值升高［4］，垃圾理化特性的改變使得現(xiàn)有焚燒爐在運(yùn)行過程中的負(fù)荷、風(fēng)量等需要進(jìn)行相應(yīng)的優(yōu)化調(diào)整［5］。

爐排式焚燒爐由于物料預(yù)處理簡單、單臺處理量大、爐渣體積小等特點，是我國生活垃圾焚燒爐中應(yīng)用最多的爐型［6］。目前對爐排式焚燒爐研究多采用爐排上產(chǎn)物預(yù)測+燃燒室內(nèi)氣相燃燒的方式，燃燒室內(nèi)的煙氣溫度、成分等取決于爐排上產(chǎn)物的預(yù)測。實際爐排上產(chǎn)物的成分與含量除了受燃燒室高溫?zé)煔獾挠绊懲?，還與一次風(fēng)分配密切相關(guān)。

爐排式焚燒爐的一次風(fēng)通常由爐排下方向上吹送，經(jīng)過燃料床層后進(jìn)入爐膛內(nèi)部，一次風(fēng)除了能為料層的固相燃燒提供氧氣之外，還具有干燥垃圾、析出揮發(fā)份、使料層穩(wěn)定著火以及冷卻爐排的作用。垃圾焚燒爐一次風(fēng)布置方式大多為母管制，從每段爐排的一次風(fēng)母管進(jìn)行分支，通過分支管路將一次風(fēng)送到爐排下部風(fēng)室，運(yùn)行過程中通過調(diào)節(jié)閥門，改變一次風(fēng)入爐量，從而改變?nèi)紵隣顩r。單個母管上通常連結(jié)3～9個支管，每個支管上裝有風(fēng)量調(diào)節(jié)閥門，由于整個一次風(fēng)系統(tǒng)是連通的，單個支管的風(fēng)量調(diào)節(jié)影響所有支管出口的風(fēng)量分配。隨著生活垃圾焚燒爐逐漸大型化，爐膛寬度逐漸變寬，燃燒區(qū)域面積變大，受母管制布風(fēng)結(jié)構(gòu)的限制，對于區(qū)域化燃燒狀況調(diào)整影響較大，調(diào)整也異常艱難［7-8］。一般管路配風(fēng)設(shè)計都是基于經(jīng)驗參數(shù)方法，查詢沿程阻力系數(shù)和局部阻力系數(shù)，來進(jìn)行復(fù)雜管路系統(tǒng)的配風(fēng)設(shè)計計算，目前尚沒有普適的經(jīng)驗公式能夠預(yù)測不同閥門開度時單個支管風(fēng)量的變化對整個一次風(fēng)管道配風(fēng)的影響。

隨著數(shù)值模擬模型的完善和計算機(jī)性能的提升，通過數(shù)值模擬的方法可以預(yù)測不同閥門開度時各個支管出口的風(fēng)量。巴海濤等［9］采用數(shù)值模擬的方法研究了船舶風(fēng)管中的空氣流動，調(diào)整了風(fēng)管走向及出口，改善了機(jī)艙風(fēng)管的風(fēng)量分配。王淑坤等［10］通過模擬汽車的空調(diào)風(fēng)道，改善了各出風(fēng)口的風(fēng)量分配情況。Kumar等［11］通過數(shù)值模擬方法研究了在不同通風(fēng)條件下紡織品蒸制過程中的風(fēng)速與溫度情況，給出了通風(fēng)優(yōu)化的建議。Zhao 等［12］通過數(shù)值模擬衛(wèi)星發(fā)射前清潔存放循環(huán)風(fēng)空調(diào)系統(tǒng)噴射方式，研究了在不同衛(wèi)星布置下循環(huán)風(fēng)量的控制策略，為減小循環(huán)風(fēng)量提供了一種可行方法。

由于爐排式焚燒爐一次風(fēng)管道系統(tǒng)的支管數(shù)量多，調(diào)節(jié)某一支管的閥門開度時，都需要重新計算，并且某一支管閥門開度的改變影響所有支管出口的流量分配，風(fēng)管閥門開度與它們各自的出口風(fēng)量是一個多變量互相耦合的關(guān)系，僅通過數(shù)值模擬很難尋到最優(yōu)工況。BP（back propagation）神經(jīng)網(wǎng)絡(luò)是一種按照誤差逆向傳播算法訓(xùn)練的多層前饋神經(jīng)網(wǎng)絡(luò)，是實現(xiàn)機(jī)器學(xué)習(xí)的主要算法之一［13］。如圖1所示，BP神經(jīng)網(wǎng)絡(luò)基于現(xiàn)場多維的離散數(shù)據(jù)，通過卡方檢驗的方法對數(shù)據(jù)降維處理，通過反向傳播不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，使網(wǎng)絡(luò)的誤差平方和最?。?4］。神經(jīng)網(wǎng)絡(luò)模型發(fā)展多變量非線性相關(guān)關(guān)系的度量來尋求與其他數(shù)據(jù)之間的函數(shù)關(guān)系，隨著數(shù)據(jù)量的增加函數(shù)關(guān)系也會變得越來越復(fù)雜，未來對于數(shù)據(jù)的預(yù)測也會越來越準(zhǔn)確。

圖1 BP神經(jīng)網(wǎng)絡(luò)示意圖Fig.1 Schematic diagram of BP neural network

神經(jīng)網(wǎng)絡(luò)模型可以克服CFD建模時多變量互相耦合、互相影響的問題，可以提高預(yù)測效率。景亮等［15］建立了出菇房溫度CFD模型，引入溫度均勻性指標(biāo)，通過CFD模擬得到不同工況下的模擬數(shù)據(jù)，并將這些數(shù)據(jù)作為神經(jīng)網(wǎng)絡(luò)的樣本數(shù)據(jù)，訓(xùn)練神經(jīng)網(wǎng)絡(luò)模型以預(yù)測出菇房的溫度均勻性指標(biāo)。楊磊磊等［16］利用數(shù)值模擬方法獲得了壓力容器中的溫度場分布，并基于數(shù)值模擬結(jié)果建立了BP神經(jīng)網(wǎng)絡(luò)模型進(jìn)行溫度預(yù)測。Mukesh等［17］利用雙流體模型對循環(huán)流化床中軸向固含率分布進(jìn)行計算，并基于數(shù)值模擬結(jié)果建立數(shù)據(jù)集訓(xùn)練神經(jīng)網(wǎng)絡(luò)模型對固含率進(jìn)行預(yù)測。神經(jīng)網(wǎng)絡(luò)的對象是數(shù)據(jù)，在諸多相互關(guān)聯(lián)的數(shù)據(jù)中，僅僅靠人工篩選和學(xué)習(xí)，很難發(fā)現(xiàn)數(shù)據(jù)對結(jié)果的影響。神經(jīng)網(wǎng)絡(luò)通過采集CFD 模擬得到的數(shù)據(jù)，積累歷史數(shù)據(jù)，使用大量的數(shù)據(jù)來訓(xùn)練模型，可在后續(xù)的計算中替代CFD模擬，降低CFD模擬的計算量。

本研究中首先根據(jù)實際一次風(fēng)管尺寸建立計算模型，通過CFD 的數(shù)值模擬得到不同支管上閥門不同開度時的出口流量分配結(jié)果并建立數(shù)據(jù)集，選取部分?jǐn)?shù)據(jù)作為訓(xùn)練樣本，其余數(shù)據(jù)作為測試集，然后使用BP 神經(jīng)網(wǎng)絡(luò)模型建立焚燒爐一次風(fēng)風(fēng)量分配預(yù)測模型并以結(jié)果的均方誤差評估，實現(xiàn)理想情況風(fēng)量分配下各個支管的閥門開度的預(yù)測。

2 計算模型及數(shù)據(jù)處理

2.1 CFD計算模型

研究對象為500 t/d 的生活垃圾焚燒爐一次風(fēng)分配管道，計算模型如圖2所示，一次風(fēng)由總管進(jìn)入3 個支路，其中干燥段和燃燒Ⅰ段在第1 支路上，燃燒Ⅱ-1 段在第2 支路上，燃燒Ⅱ-2 段和燃燼段在第3 支路上。爐排深度方向按照三列布置，每列的干燥段下有一個配風(fēng)倉室，每列燃燒Ⅱ-1和燃燒Ⅱ-1段下各有兩個配風(fēng)倉室，每列燃燼段下有一個配風(fēng)倉室，每個倉室由3個支路引出通過支管連接，共計18個支管。計算模型如圖2所示，采用六面體網(wǎng)格，網(wǎng)格數(shù)約為131萬。主要模型參數(shù)見表1。

表1 主要模型參數(shù)Tab.1 Main model parameters

圖2 計算模型Fig.2 Schematic diagram of geometric model

2.2 控制方程

計算過程中假定整個一次風(fēng)系統(tǒng)中氣體溫度恒定，氣相控制方程主要包括質(zhì)量守恒方程和動量守恒方程。

氣相質(zhì)量守恒方程：

式中：ρ為氣相的密度，單位為kg/m3；u為速度矢量，單位為m/s。

氣相動量守恒方程：

式中：ρ為壓力；i、j為坐標(biāo)系中的方向；g為重力的分量，單位為m/s2；Fi為阻力等其他源項；μ為流體的粘度，單位為kg·m-1·s-1。

采用k-ε 模型求解湍流的影響，其中湍動能k和耗散率ε方程如下：

式中：GK是速度梯度產(chǎn)生的湍動能；Gb為浮力作用產(chǎn)生的湍流動能；YM為可壓縮湍流中的脈動擴(kuò)張；G1ε、G2ε為系數(shù)；Sε為源項。

2.3 邊界條件

采用三維穩(wěn)態(tài)模型，選擇Coupled 模型計算壓力與速度之間的耦合。設(shè)定氣體入口為流量入口邊界條件，出口為壓力出口邊界條件，壁面為無滑移絕熱邊界條件。

閥門開度是調(diào)節(jié)風(fēng)量的重要參數(shù)，閥門開度實際影響的是管道的流通面積，在模擬過程種，通過在管道中設(shè)置擋板調(diào)節(jié)管道的流通面積，計算不同閥門開度的影響。

2.4 BP神經(jīng)網(wǎng)絡(luò)模型

2.4.1 BP神經(jīng)網(wǎng)絡(luò)模型的建立

神經(jīng)網(wǎng)絡(luò)中隱含層的神經(jīng)節(jié)點數(shù)與BP 神經(jīng)網(wǎng)絡(luò)的非線性映射能力呈正相關(guān)關(guān)系，網(wǎng)絡(luò)結(jié)構(gòu)設(shè)置影響預(yù)測結(jié)果的準(zhǔn)確程度［18］。與單個隱含層數(shù)量相比，擁有兩個及以上隱含層單元數(shù)的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)在進(jìn)行訓(xùn)練測試時更容易出現(xiàn)局部極值，本研究選用單一隱含層的三層BP 神經(jīng)網(wǎng)絡(luò)，輸入層共18 個節(jié)點，分別為18 個支管出口的風(fēng)量，輸出層同樣為18個節(jié)點，分別為18個出口對應(yīng)的閥門開度。

隱含層神經(jīng)元是完成從輸入層神經(jīng)元到輸出層神經(jīng)元非線性映射的關(guān)鍵。若隱層節(jié)點數(shù)太少，網(wǎng)絡(luò)可能根本不能訓(xùn)練或網(wǎng)絡(luò)性能很差；若隱層節(jié)點數(shù)太多，增加隱藏層中神經(jīng)元的數(shù)量也并不意味著模型將具有更好的預(yù)測精度［19］，一方面使網(wǎng)絡(luò)訓(xùn)練時間延長，另一方面，訓(xùn)練容易陷入局部極值而得不到最優(yōu)點，也是在神經(jīng)網(wǎng)絡(luò)訓(xùn)練時出現(xiàn)“過擬合”的內(nèi)在原因。公式（5）［20］給出了隱含層節(jié)點數(shù)的計算方法：

式中：m和n分別為輸入層和輸出層的節(jié)點個數(shù)；a的取值范圍為［1，10］。根據(jù)輸入層神經(jīng)元個數(shù)是18，輸出層神經(jīng)元的個數(shù)是18，可得神經(jīng)網(wǎng)絡(luò)模型隱含層節(jié)點數(shù)為［7，16］。本研究從7 個節(jié)點開始，再分別增加至16 個來測量不同數(shù)量的節(jié)點數(shù)對模型誤差的影響，最后選擇均方誤差最小的數(shù)量作為隱含層節(jié)點數(shù)，均方誤差公式為：

式中：n為訓(xùn)練樣本數(shù)量；f(xi)為原始數(shù)據(jù)集中輸出數(shù)據(jù)，yi為訓(xùn)練后輸出數(shù)據(jù)。

訓(xùn)練函數(shù)選擇Levenberg-Marquardt（trainlm 訓(xùn)練函數(shù)），學(xué)習(xí)率為0.01，迭代次數(shù)設(shè)為1 000，誤差閾值設(shè)為0.001。使用平均絕對誤差（MAE）和決定系數(shù)（R2）來評估所建立的神經(jīng)網(wǎng)絡(luò)性能，公式如下：

式中：Yi為原始數(shù)據(jù)集中輸出數(shù)據(jù)，為訓(xùn)練后輸出數(shù)據(jù)，為原始數(shù)據(jù)集中輸出數(shù)據(jù)平均值，n為數(shù)據(jù)數(shù)量。

2.4.2 數(shù)據(jù)預(yù)處理

在數(shù)據(jù)集構(gòu)建過程中，設(shè)置了不同支路上閥門的不同開度，通過統(tǒng)計不同閥門開度下各個出口處的風(fēng)量，計算爐排不同段的風(fēng)量分配比例，共獲得了65 條數(shù)據(jù)，其中70%作為訓(xùn)練集，30%作為測試集。由于原始樣本中各向量的數(shù)量級差別很大，為方便計算并防止部分神經(jīng)元達(dá)到過飽和狀態(tài)，對訓(xùn)練集數(shù)據(jù)進(jìn)行了歸一化處理。本文采用MATLAB R2021a的內(nèi)置函數(shù)mapminman 函數(shù)對數(shù)據(jù)進(jìn)行歸一化處理，將數(shù)據(jù)歸一化到［0，1］，歸一化函數(shù)公式為：

式中：y為輸出值，x為所選數(shù)據(jù)，xmax、xmin為數(shù)據(jù)集中最大值、最小值，ymax、ymin為1、0。

3 結(jié)果與分析

3.1 CFD模擬結(jié)果

由于不同爐排段上的出口3 離入口較遠(yuǎn)，受網(wǎng)格影響大，所以采用出口3 處的風(fēng)量比例驗證網(wǎng)格質(zhì)量。圖3 給出了網(wǎng)格數(shù)分別為約7 萬、34 萬、131萬和413 萬時的不同爐排段出口3 處的風(fēng)量比例，可以看出網(wǎng)格數(shù)為7 萬和34 萬時，出口處的風(fēng)量比例與網(wǎng)格數(shù)為413 萬相差較大，網(wǎng)格數(shù)為131 萬時，出口處的風(fēng)量與網(wǎng)格數(shù)為413 萬時相近，為了提高計算效率，模擬過程中采用131的網(wǎng)格進(jìn)行計算。

圖3 不同網(wǎng)格數(shù)的出口風(fēng)量分配Fig.3 Distribution of air flow at the outlet with different grids

圖4給出了閥門全開時一次風(fēng)系統(tǒng)的速度矢量圖。從圖中可以看出，3個支路管內(nèi)的速度分布并不均勻，由于氣流的慣性，3個支路氣流都會吹到管路的右側(cè)壁面，再沿壁面反射回來，形成類似旋轉(zhuǎn)的氣流形式。18個支管內(nèi)速度分布也不均勻，與3個支路管道類似，氣流會吹到沿流動方向?qū)γ娴墓鼙谏?，其中干燥段和燃燒?1段的支管內(nèi)速度不均勻性較大。

圖4 氣體速度矢量圖Fig.4 Gas velocity distribution

圖5 分別給出了閥門全開和干燥段閥門開度93.3%時不同爐排段的風(fēng)量分配。閥門全開時6 個爐排段的風(fēng)量分配分別約為8.4%、17.7%、23.2%、19.0%、21.7%、10.0%。其中第4支路的管徑較細(xì)，阻力大，造成燃燒Ⅱ-1 段的風(fēng)量較低。圖中同時給出了現(xiàn)場運(yùn)行時的風(fēng)量分配，現(xiàn)場運(yùn)行時各支管的風(fēng)量與閥門開度見表2。

表2 現(xiàn)場運(yùn)行時各支管風(fēng)量及閥門開度Tab.2 Actual values of air volume and valve opening of each branch pipe

圖5 不同段的爐排風(fēng)量分配Fig.5 Air volume distribution of different grate sections

模擬計算得到的風(fēng)量分配與現(xiàn)場總體趨勢一致，風(fēng)量集中在燃燒段。由于現(xiàn)場運(yùn)行時爐排上堆積了料層，不同燃燒段料層的厚度和料層的孔隙度不同，同時現(xiàn)場運(yùn)行時各閥門的開度根據(jù)爐排上垃圾料層厚度和推料速度調(diào)節(jié)有所不同，因此各個支管的風(fēng)量分配比例不盡相同，表2中的數(shù)據(jù)為現(xiàn)場運(yùn)行時DCS 上實際數(shù)據(jù)。當(dāng)干燥段開度設(shè)為93.3%時，干燥段風(fēng)量相對減少，其余支路的風(fēng)量有不同程度的變化，這也說明焚燒爐一次風(fēng)管道配風(fēng)模型t各支路之間存在著耦合關(guān)系，是一個多變量互相影響的過程，當(dāng)改變一個支管的閥門開度時，其余支管的出口風(fēng)量都將產(chǎn)生變化，若想得到理想配風(fēng)比的情況，僅僅靠CFD模擬需要消耗大量的計算時間。

3.2 神經(jīng)網(wǎng)絡(luò)預(yù)測結(jié)果

不同隱含層神經(jīng)元數(shù)量得到的訓(xùn)練結(jié)果均方誤差如圖6所示，從圖中可以看出，隨著隱含層節(jié)點數(shù)的增多，訓(xùn)練結(jié)果均方誤差呈先下降后升高的趨勢，隱含層神經(jīng)元個數(shù)為13 時BP 神經(jīng)網(wǎng)絡(luò)對函數(shù)的逼近效果最好，故本研究將網(wǎng)絡(luò)的隱含層神經(jīng)元數(shù)目定為13。

圖6 不同隱含層神經(jīng)元數(shù)量訓(xùn)練結(jié)果均方誤差Fig.6 Mean square error values of different training results

將測試集樣本的輸入數(shù)據(jù)進(jìn)行歸一化處理，并將測試集的輸出數(shù)據(jù)代入已訓(xùn)練完的網(wǎng)絡(luò)中進(jìn)行仿真測試，即利用MATLAB R2021a 中的simulink 仿真測試函數(shù)將測試集樣本數(shù)據(jù)錄入進(jìn)BP 神經(jīng)網(wǎng)絡(luò)當(dāng)中然后再經(jīng)過reverse 反歸一化函數(shù)把得出的預(yù)期結(jié)果輸出數(shù)據(jù)反歸一化來實現(xiàn)仿真測試。根據(jù)圖7所示，BP神經(jīng)網(wǎng)絡(luò)在完成了14次訓(xùn)練好已經(jīng)達(dá)到了穩(wěn)定狀態(tài)，且滿足了誤差閾值為0.001 的條件，最佳驗證性能為第8輪的訓(xùn)練。

圖7 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練結(jié)果及訓(xùn)練次數(shù)Fig.7 BP neural network training results and training times

仿真測試集樣本數(shù)據(jù)在測試完成時模型的訓(xùn)練結(jié)果對比和解釋度分別如圖8和圖9所示，可以看出基于BP 神經(jīng)網(wǎng)絡(luò)的一次風(fēng)系統(tǒng)分配預(yù)測模型的解釋度R=0.98724，說明仿真測試的預(yù)測結(jié)果與其實際值間的相關(guān)程度高，這也表明模型的訓(xùn)練是準(zhǔn)確的。

圖8 測試完成時模型的訓(xùn)練結(jié)果對比Fig.8 Comparison of training results of models

圖9 測試完成時模型的解釋度Fig.9 The interpretability of the model

仿真測試集樣本數(shù)據(jù)在測試完成時訓(xùn)練集與測試集的平均絕對誤差和決定系數(shù)分別如圖10 和圖11 所示，在進(jìn)行多次重復(fù)訓(xùn)練后，訓(xùn)練集與測試集的平均相對誤差分別約為1.24%和2.01%，訓(xùn)練集與測試集的平均相對誤差分別最低達(dá)到了0.83%和1.49%。訓(xùn)練集與測試集的決定系數(shù)分別約為0.98 和0.91，訓(xùn)練集與測試集的決定系數(shù)分別最高達(dá)到了0.988和0.981，這也表明訓(xùn)練質(zhì)量較高，模型能較好表現(xiàn)出數(shù)據(jù)之間的非線性映射關(guān)系，所建立的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)比較合適，網(wǎng)絡(luò)可以較為準(zhǔn)確地預(yù)測在理想配風(fēng)比下各個出口的閥門開度。

圖10 訓(xùn)練集與測試集的平均絕對誤差Fig.10 Average absolute error value between training and testing sets

圖11 訓(xùn)練集與測試集的決定系數(shù)Fig.11 Coefficient of determination between training and testing sets

4 結(jié)論

本文提出了一種CFD-BP 耦合算法，基于CFD計算得到一次風(fēng)系統(tǒng)中各支管的風(fēng)量分配，采用BP 算法預(yù)測不同支路閥門開度對整體管路配風(fēng)的影響，該耦合算法可為建立智能配風(fēng)設(shè)計與設(shè)備運(yùn)行智能管理提供一定的理論支撐，主要結(jié)論如下：

（1）建立了某500 t/d 的生活垃圾焚燒爐一次風(fēng)管道分配CFD 計算模型，計算了在不同閥門開度下各個支管出口的流量分配情況，并與實際結(jié)果進(jìn)行了比對，模擬計算得到的風(fēng)量分配與現(xiàn)場總體趨勢一致。

（2）以CFD 計算結(jié)果為數(shù)據(jù)庫，建立了以各個支管出口風(fēng)量為輸入量的一次風(fēng)管道分配BP 神經(jīng)網(wǎng)絡(luò)模型，模型訓(xùn)練集與測試集的平均相對誤差分別最低達(dá)到了0.83%和1.49%，訓(xùn)練集與測試集的決定系數(shù)分別最高達(dá)到了0.988 和0.981，預(yù)測效果較好。