伍 權(quán),朱 萌,陳 磊,卿浩然,周琳剛,許 凱, 蘇 勝,胡 松,劉 輝,向 軍
(1.國家能源集團新能源技術(shù)研究院有限公司,北京 102209; 2.華中科技大學(xué)煤燃燒國家重點實驗室,湖北 武漢 430074)
超臨界二氧化碳(S-CO2)具有穩(wěn)定、易得、無毒、壓縮性好和流動性強等特點,相比于傳統(tǒng)水循環(huán),基于S-CO2工質(zhì)的新型動力循環(huán)具有體積小、效率高、靈活性好[1]等優(yōu)勢。目前,基于S-CO2循環(huán)的研究在光熱發(fā)電[2-3]、火力發(fā)電[4-5]、核能[6-7]以及余熱利用[8-9]等方面發(fā)展迅速。S-CO2新型動力循環(huán)的工質(zhì)阻力特性會對S-CO2循環(huán)效率以及各部件設(shè)計產(chǎn)生較大影響,若無法準(zhǔn)確計算循環(huán)管路和主要部件的阻力,S-CO2系統(tǒng)將無法準(zhǔn)確設(shè)計并高效可靠運行。因此,S-CO2的流動阻力特性可靠評估是S-CO2技術(shù)商業(yè)化推進中至關(guān)重要的一環(huán)。
由于湍流的脈動性和新型動力系統(tǒng)實驗數(shù)據(jù)缺乏,現(xiàn)有研究尚未形成對S-CO2摩擦阻力特性的統(tǒng)一理論。Kirillov等人[10]認(rèn)為管截面密度分布是影響非等溫流體阻力特性的關(guān)鍵因素,提出了壁面密度(依靠內(nèi)壁溫計算的流體密度)和主流密度比修正項來對等溫流體公式進行修正。Fang等人[11]利用現(xiàn)有文獻中的S-CO2以及其他流體數(shù)據(jù),采用大量函數(shù)類型進行回歸測試,提出了壁面黏度與主流黏度比修正項以及膜密度與擬臨界點密度比的指數(shù)修正項。張海松等[12]研究了S-CO2在擬臨界點附近的阻力系數(shù)峰值現(xiàn)象,認(rèn)為傳熱惡化會對流動阻力特性造成負(fù)面影響,并采用超臨界沸騰數(shù)SBO和普朗特數(shù)比作為修正項計算S-CO2摩擦阻力。然而,上述各修正項可能在一定的參數(shù)范圍內(nèi)適用,而在其他參數(shù)范圍將不再是決定性因素,這給S-CO2流動阻力計算的工程化應(yīng)用造成困難,設(shè)計人員難以根據(jù)需要采取合理修正項。
基于此,人工神經(jīng)網(wǎng)絡(luò)算法被認(rèn)為是實現(xiàn)S-CO2阻力預(yù)測的一種有效途徑。該方法不考慮流動傳熱的基本理論,而是通過大量訓(xùn)練,學(xué)習(xí)數(shù)據(jù)之間的特征和規(guī)律,并基于該規(guī)律預(yù)測其他數(shù)據(jù)的特征。這種方式的泛化能力強,不再局限于某一特定公式,可以預(yù)測更廣參數(shù)范圍的阻力特征。同時,采用神經(jīng)網(wǎng)絡(luò)可以直接將測量值例如溫度、壓力等作為輸入,阻力作為輸出,不再求解不同狀態(tài)下的特征參數(shù)例如雷諾數(shù)、普朗特數(shù)等,有利于工程化實現(xiàn)。Ma等人[13]收集了14篇文獻中的超臨界水傳熱數(shù)據(jù),利用BP(誤差反向傳播)神經(jīng)網(wǎng)絡(luò)進行了數(shù)據(jù)訓(xùn)練以預(yù)測其對流換熱系數(shù),結(jié)果表明其均方差為4.13%,相關(guān)系數(shù)R為0.971 7。Balcilar等人[14]基于垂直光滑銅管中的R134a冷卻數(shù)據(jù),訓(xùn)練了多個神經(jīng)網(wǎng)絡(luò)模型以預(yù)測其傳熱和壓降特性,預(yù)測誤差可保證在5%以內(nèi)。Shadloo等人[15]采用多層感知器神經(jīng)網(wǎng)絡(luò)對長管道中的氣體/非牛頓流體兩相流的阻力特性進行了研究,預(yù)測值與實驗數(shù)據(jù)的偏差在大多數(shù)情況下不超過5%?,F(xiàn)有研究中針對S-CO2阻力特性進行神經(jīng)網(wǎng)絡(luò)預(yù)測的相關(guān)文獻相對較少,神經(jīng)網(wǎng)絡(luò)對于其流動阻力的預(yù)測能力有待驗證。此外,現(xiàn)有非絕熱S-CO2實驗數(shù)據(jù)[12,16-18]主要集中在CO2臨界點(7.38 MPa、31 ℃)[19-20]附近,且質(zhì)量流速較低(<1 500 kg/(m2·s))、流體溫度較低(<200 ℃),無法滿足光熱、火電以及核能等領(lǐng)域的實際需求。
聚焦于解決S-CO2發(fā)電領(lǐng)域中S-CO2流動摩擦阻力特性的預(yù)測問題,本文對更寬參數(shù)范圍(750~ 2 200 kg/(m2·s)、10~20 MPa、200~340 kW/m2、60~500 ℃)的S-CO2阻力特性進行實驗測量,分析壓力、質(zhì)量流速以及熱流密度對摩擦阻力的影響,并根據(jù)實驗數(shù)據(jù),建立貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)模型,以實現(xiàn)對S-CO2摩擦阻力特性的準(zhǔn)確預(yù)測。
S-CO2流動阻力特性研究實驗系統(tǒng)如圖1所示。該系統(tǒng)首先通過加壓柱塞泵將儲罐中的低壓CO2(純度99.99%)增壓至S-CO2循環(huán)回路。S-CO2循環(huán)回路由循環(huán)泵克服流動阻力,S-CO2依次流過質(zhì)量流量計、回?zé)崞骼涠?、預(yù)熱器、實驗段、回?zé)崞鳠岫艘约袄鋮s器,最終再次進入循環(huán)泵形成完整循環(huán)。冷卻器采用與常規(guī)火電站相同的水冷循環(huán)實現(xiàn)對S-CO2的冷卻,以保證進入循環(huán)泵的S-CO2在50 ℃以內(nèi)。
實驗時,通過加壓柱塞泵和泄壓閥配合實現(xiàn)對實驗循環(huán)壓力的調(diào)節(jié);通過控制循環(huán)泵頻率實現(xiàn)對實驗循環(huán)流量的調(diào)節(jié);通過控制預(yù)熱器和實驗段的加熱功率實現(xiàn)對實驗段進口溫度以及實驗段熱流密度的調(diào)節(jié)。系統(tǒng)設(shè)計參數(shù)為:最大工質(zhì)壓力30 MPa,最大加熱功率150 kW,最大流量1 800 L/h。
實驗段設(shè)計如圖1所示。實驗段管材為Inconel 625合金鋼,管內(nèi)徑為12 mm,外徑為20 mm,在2個電極板之間的加熱段長度為1 870 mm。
除加熱段外,兩端各留有1 000 mm的過渡段,用于穩(wěn)定加熱段的流動。實驗段采用電阻發(fā)熱原理,將低電壓大電流電源連接到實驗段兩端,即2個電極板處,鋼管本身的電阻使其產(chǎn)生均勻的加熱熱流。同時,實驗段管外壁包裹了保溫層以降低散熱損失。加熱段沿管長方向均勻布置有28個熱電偶,用于測量管外壁溫度。實驗段的入口和出口分別安裝有測量工質(zhì)溫度的熱電偶和測量工質(zhì)壓力及壓差的取壓孔。實驗段的流向為垂直向上,與重力方向相反。
實驗中,工質(zhì)溫度、壁溫測量分別采用Omega 工程公司生產(chǎn)的K型鎧裝熱電偶、K型熱電偶絲,測溫精度為0.75%。工質(zhì)壓力和壓差測量采用北京俄華通儀表技術(shù)有限公司生產(chǎn)的壓力、壓差變送器,測壓精度為0.25%F.S。工質(zhì)流量測量采用北京首科實華自動化設(shè)備有限公司生產(chǎn)的質(zhì)量流量計,測量精度0.2%。通過模擬信號將上述數(shù)據(jù)傳輸至由研華科技(中國)有限公司生產(chǎn)的數(shù)據(jù)采集卡,并轉(zhuǎn)換為數(shù)字信號傳入電腦。該數(shù)據(jù)采集卡對模擬信號的識別精度為0.1%,數(shù)據(jù)采集間隔為2 s。當(dāng)實驗循環(huán)各關(guān)鍵參數(shù)穩(wěn)定后,開始進行數(shù)據(jù)采集,選取穩(wěn)定狀態(tài)下的60 s數(shù)據(jù)取平均后作為該次測量的最終數(shù)據(jù),以避免波動誤差。
2.1.1 阻力計算
由于實驗管段豎直,因此實驗系統(tǒng)中測量的壓差包括摩擦阻力ΔPf、重力壓降ΔPg和加速壓降ΔPa。摩擦阻力ΔPf的測量值可由式(1)計算。
式中:ΔP為總壓降,kPa。
摩擦阻力ΔPf還可通過理論計算得到:
式中:f為摩擦阻力系數(shù);L為管長,m;D為管內(nèi)徑,m;G為流體的質(zhì)量流速,kg/(m2·s);ρ為流體密度,kg/m3。
對于重力壓降ΔPg,由于S-CO2的密度變化在實驗參數(shù)范圍內(nèi)是非線性的,因此直接通過簡單平均處理存在較大誤差,本文采用目前認(rèn)可度較高的Ornatskiy公式[21]來計算:
式中:h為比焓,kJ/kg;g為重力加速度,m/s2。
對于加速壓降ΔPa,采用式(4)計算:
流體比焓、密度等物性參數(shù)通過NIST REFPROP 9.1軟件查詢獲取。實驗段沿長度方向不同位置z的主流比焓為:
式中:hb,z為位置z的主流比焓,kJ/kg;z為沿長度方向距實驗段入口的距離,m。
實驗段的加熱熱流密度指內(nèi)壁面?zhèn)鬟f給流體的熱流密度:
式中:qwall為內(nèi)壁面熱流密度,kW/m2。
由于實驗只能測量管外壁溫,因此管內(nèi)壁溫需要通過間接計算得到。假定管周向和軸向的導(dǎo)熱是均勻的,可將本問題簡化為內(nèi)熱源一維穩(wěn)態(tài)導(dǎo)熱,最終內(nèi)壁溫計算式如下:
式中:t為溫度,℃;λ為金屬導(dǎo)熱系數(shù),kW/(m·℃);d為管外徑,m。
2.1.2 實驗參數(shù)與不確定度
實驗主要從質(zhì)量流速、流體壓力以及熱流密度3個方面來探索摩擦阻力的差異。此外,流體溫度和壁面溫度也是測量參數(shù)。其中質(zhì)量流速、熱流密度和摩擦阻力等屬于間接測量值,壓力和溫度等屬于直接測量值??紤]到數(shù)據(jù)采集卡存在識別精度,即使是直接測量值仍需做不確定度計算。表1給出了各參數(shù)在實驗時的數(shù)值范圍以及不確定度。
表1 實驗參數(shù)范圍與不確定度 Tab.1 Range and uncertainty of experimental parameters
各參數(shù)不確定度采用式(8)[16]計算:
BP神經(jīng)網(wǎng)絡(luò)具有信號的前向傳播和誤差的反向傳播2個過程。當(dāng)輸入信號經(jīng)過隱含層、輸出層獲得的輸出信號與期望輸出不符時,將輸出誤差反向從輸出層逐層傳回輸入層,每層每個神經(jīng)元所獲得的誤差數(shù)據(jù)將作為調(diào)整其自身權(quán)值的依據(jù),最終使得誤差逐漸下降至目標(biāo)范圍或不再下降,訓(xùn)練完成。
本文BP神經(jīng)網(wǎng)絡(luò)模型采用MATLAB編程實現(xiàn),設(shè)計了2層前饋神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)(圖2),可以實現(xiàn)對任意“一個有限空間到另一個有限空間的連續(xù)映射函數(shù)”的擬合[22]。更高的層數(shù)雖然理論上擬合能力更強,但存在過擬合以及訓(xùn)練難度增加的問題。隱含層神經(jīng)元采用tansig非線性激活函數(shù),輸出層神經(jīng)元采用purelin線性激活函數(shù):
在常規(guī)BP神經(jīng)網(wǎng)絡(luò)中,一般采用神經(jīng)網(wǎng)絡(luò)響應(yīng)的均方差作為目標(biāo)函數(shù),通過最優(yōu)化算法實現(xiàn)均方差的逐漸減小,當(dāng)均方差小于所設(shè)定的范圍或不再繼續(xù)下降時,認(rèn)為神經(jīng)網(wǎng)絡(luò)已完成訓(xùn)練。然而,該方式易出現(xiàn)過擬合現(xiàn)象,使得神經(jīng)網(wǎng)絡(luò)泛化能力較差,無法達到較好的預(yù)測效果。與常規(guī)BP神經(jīng)網(wǎng)絡(luò)不同,貝葉斯正則化算法將目標(biāo)函數(shù)修改為:
式中:α和β為正則化系數(shù);Ew為神經(jīng)網(wǎng)絡(luò)中所有權(quán)值的均方差;Ed為常規(guī)神經(jīng)網(wǎng)絡(luò)的目標(biāo)函數(shù)。
設(shè)計α和β的意義在于平衡神經(jīng)網(wǎng)絡(luò)擬合效果與泛化能力之間的矛盾。當(dāng)α過大時,訓(xùn)練具有更好的泛化能力,但擬合精度不高,可能忽略了較多的數(shù)據(jù)特征;當(dāng)β過大時,對訓(xùn)練集的響應(yīng)均方差最小,但容易導(dǎo)致過擬合而對非訓(xùn)練集無法實現(xiàn)較好的預(yù)測效果。在貝葉斯理論中,假定訓(xùn)練集與權(quán)重集的先驗概率服從高斯分布,最優(yōu)的權(quán)值應(yīng)具有最大后驗概率,而最大后驗概率在上述條件下可以等價于最小化F,求得F最小時的α和β為:
式中:γ為有效網(wǎng)絡(luò)參數(shù)的數(shù)量;n為數(shù)據(jù)個數(shù)。該算法更適合訓(xùn)練具有噪聲和數(shù)據(jù)量少的數(shù)據(jù)集[23]。由于阻力實驗各參數(shù)均存在不確定度,且數(shù)據(jù)集規(guī)模對于可訓(xùn)練億級數(shù)據(jù)的神經(jīng)網(wǎng)絡(luò)來說相對較小,因此采用貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)作為阻力預(yù)測模型,以提高模型泛化能力。
在模型中,將流體壓力、質(zhì)量流速、熱流密度和主流溫度4個參數(shù)作為輸入?yún)?shù),摩擦阻力系數(shù)作為輸出參數(shù)。為提高模型預(yù)測精度,對輸入?yún)?shù)和輸出參數(shù)進行了歸一化處理,歸一化后各參數(shù)范圍均在0~1。歸一化方法如下:
實驗總共包含194組數(shù)據(jù)集,模型訓(xùn)練集設(shè)定為165組,測試集為29組(隨機選?。?。隱含層神經(jīng)元個數(shù)的選取采用經(jīng)驗公式和重復(fù)訓(xùn)練測試聯(lián)合完成,常用公式為[24]:
式中:Nhidden為隱含層神經(jīng)元個數(shù);Nin、Nout分別為輸入、輸出參數(shù)個數(shù);X為經(jīng)驗系數(shù),一般取1~10。本模型中,輸入?yún)?shù)有4個,輸出參數(shù)為1個,因此推薦的隱含層神經(jīng)元個數(shù)為3~12。為選擇最佳隱含層神經(jīng)元個數(shù),進行重復(fù)訓(xùn)練測試,每種神經(jīng)元個數(shù)方案重復(fù)訓(xùn)練100次(每次的測試集是隨機的)并對100次訓(xùn)練的性能指標(biāo)取平均作為最終值來對比不同神經(jīng)元個數(shù)的效果。性能評價指標(biāo)為Pearson相關(guān)系數(shù)R:
式中:x、y分別為模型計算值和實驗測量值;n為數(shù)據(jù)個數(shù)。相關(guān)系數(shù)R的取值范圍為0~1,其值越高,表明模型預(yù)測效果越好。
圖3顯示了不同隱含層神經(jīng)元個數(shù)對相關(guān)系數(shù)R的影響。由圖3可以發(fā)現(xiàn),初始時隨著隱含層神經(jīng)元個數(shù)的增加,R逐漸增加,但當(dāng)隱含層神經(jīng)元個數(shù)超過11時,R幾乎不再增加。因此,選取隱含層神經(jīng)元個數(shù)為11,以確保擬合精度并避免過擬合。在該隱含層神經(jīng)元個數(shù)下的100次重復(fù)訓(xùn)練中,選取R最高的訓(xùn)練神經(jīng)網(wǎng)絡(luò)作為后續(xù)計算分析使用。
為驗證貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)在預(yù)測S-CO2流動阻力時的泛化優(yōu)勢,采用新測試數(shù)據(jù)集(不在原165組訓(xùn)練集和29組測試集之內(nèi))進行性能測試。新測試數(shù)據(jù)集的參數(shù)范圍見表2。
表2 新測試數(shù)據(jù)集的參數(shù)范圍 Tab.2 Parameter range of new test data set
新測試數(shù)據(jù)集對于模型屬于未知數(shù)據(jù)集,總共35組,其中21組為文獻數(shù)據(jù),剩余14組為本實驗系統(tǒng)提供。此外,有27組(21組文獻+6組本實驗系統(tǒng))為超過原始194組數(shù)據(jù)集參數(shù)范圍限制的數(shù)據(jù)(即超限點)。
圖4給出了貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)與常規(guī)BP神經(jīng)網(wǎng)絡(luò)對新測試數(shù)據(jù)集的預(yù)測結(jié)果對比。常規(guī)BP神經(jīng)網(wǎng)絡(luò)除未使用貝葉斯正則化算法外,其余設(shè)定以及訓(xùn)練均與貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)相同。圖4中E±10%表示與實驗值誤差在10%以內(nèi)的預(yù)測值個數(shù)占總個數(shù)的百分比,E±30%表示與實驗值誤差在30%以內(nèi)的預(yù)測值個數(shù)占總個數(shù)的百分比。
由圖4可以看出:貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)的相關(guān)系數(shù)、E±10%和E±30%等指標(biāo)均優(yōu)于常規(guī)神經(jīng)網(wǎng)絡(luò);對于超限點,貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)同樣具有較好的預(yù)測性能,而常規(guī)BP網(wǎng)絡(luò)的誤差則相對較大。綜上,貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)相比常規(guī)BP神經(jīng)網(wǎng)絡(luò)具有更好的泛化能力,可滿足工業(yè)應(yīng)用要求。
在以往的研究中,非等溫流體(管壁存在加熱或冷卻)經(jīng)驗公式主要是在等溫流體公式的基礎(chǔ)上進行修正,考慮管截面方向的物性差異。Filonenko公式[25]是目前應(yīng)用最廣的等溫流體公式,在其后的很多研究者基于該公式進一步修正,獲得非等溫條件下的經(jīng)驗公式。Kirillov等人[10]給出了基于密度比修正的非等溫流體公式,張海松等[12]給出了基于超臨界沸騰數(shù)SBO和普朗特數(shù)比的非等溫流體公式。Churchill[26]和Haaland[27]給出了另外2種常見的等溫流體公式。為驗證BP神經(jīng)網(wǎng)絡(luò)模型的實際預(yù)測效果,將其計算結(jié)果與表3中的5個典型摩擦系數(shù)公式計算結(jié)果進行對比,如圖5所示。
表3 典型摩擦系數(shù)公式 Tab.3 Typical friction coefficient formula
由圖5可知:當(dāng)摩擦阻力較小時,各公式的預(yù)測效果均較好,但當(dāng)摩擦阻力逐步增加后,預(yù)測值和實驗值差距開始變大;在本實驗數(shù)據(jù)范圍內(nèi),3個等溫流體公式的預(yù)測值基本一致,3個評價指標(biāo)值也基本相同;5個典型公式計算值與摩擦阻力實驗值的相關(guān)系數(shù)均在0.98~0.99,屬于極高相關(guān);在5個典型公式中,Kirillov等人的非等溫流體公式預(yù)測效果較好,E±10%=28.35%,E±30%=86.08%;與5個典型公式相比,BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測效果最好,相關(guān)系數(shù)可達0.998 6,E±10%=98.97%,E±30%=100%。
圖6給出了不同質(zhì)量流速G(765.56、1 445.56、2160.71 kg/(m2·s))下S-CO2摩擦阻力實驗值以及預(yù)測值。由圖6可知,隨著質(zhì)量流速G的增加,S-CO2摩擦阻力增加,且摩擦阻力-主流溫度曲線的斜率也逐漸增大。這是由于摩擦阻力與質(zhì)量流速的平方G2呈正比關(guān)系。溫度變化主要影響摩擦阻力系數(shù)和密度,將主流溫度作為自變量時,其斜率主要受G2(正相關(guān))、密度(負(fù)相關(guān))以及摩擦阻力系數(shù)(正相關(guān))的影響。因此,G增加,其斜率也逐漸增大。
圖7給出了不同熱流密度(208.56、279.14、335.10 kW/m2)下S-CO2摩擦阻力實驗值以及預(yù)測值。由圖7可知,隨熱流密度的增加,摩擦阻力下降。熱流密度主要影響了管截面上的流體物性分布,即對于同一主流溫度,熱流密度主要影響摩擦阻力系數(shù)。隨著熱流密度的增加,根據(jù)實驗現(xiàn)象,其摩擦阻力系數(shù)減小。因此,在BP神經(jīng)網(wǎng)絡(luò)模型中,考慮熱流密度的影響是必要的。
圖8給出了不同流體壓力(10.47、15.53、19.52 MPa)下S-CO2摩擦阻力實驗值以及預(yù)測值。由圖8可以看出,隨著壓力的增加,摩擦阻力逐漸減小,且該現(xiàn)象在較低壓力下更明顯。根據(jù)摩擦阻力公式,當(dāng)管長L、管徑D和質(zhì)量流速G不變時,摩擦阻力主要受摩擦阻力系數(shù)f(正相關(guān))和密度ρ(負(fù)相關(guān))影響。隨壓力的增加,S-CO2密度增加,黏度增加。密度增加會使S-CO2體積流量下降,流速下降,摩擦阻力下降;但黏度增加會導(dǎo)致流體間切應(yīng)力增加,流體間能量損耗增大,摩擦阻力系數(shù)增加。綜合來看,密度影響更關(guān)鍵(可能的原因是黏度對摩擦阻力系數(shù)的影響并非線性正比例關(guān)系且密度同樣影響摩擦阻力系數(shù)),因此摩擦阻力下降。
1)S-CO2摩擦阻力隨熱流密度的增加而減小,即其摩擦阻力系數(shù)隨熱流密度的增加而減小,這有利于S-CO2循環(huán)發(fā)電的大規(guī)模工程應(yīng)用。S-CO2摩擦阻力隨流體壓力的增加而減小,主要原因在于壓力增加導(dǎo)致流體密度增加,而黏度增加導(dǎo)致的阻力增加不占主導(dǎo)地位。S-CO2摩擦阻力隨質(zhì)量流速的變化情況符合傳統(tǒng)規(guī)律。
2)完成數(shù)據(jù)訓(xùn)練的貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)模型能夠有效預(yù)測S-CO2的摩擦阻力,該方法泛化能力強、擬合度高且可以避免求解不同狀態(tài)下的SCO2特征參數(shù),簡化了工程應(yīng)用,為S-CO2阻力預(yù)測提供了新的思路。