基于貝葉斯正則化BP神經網絡的 超臨界二氧化碳流動阻力特性預測

伍 權，朱 萌，陳 磊，卿浩然，周琳剛，許 凱， 蘇 勝，胡 松，劉 輝，向 軍

（1.國家能源集團新能源技術研究院有限公司，北京 102209； 2.華中科技大學煤燃燒國家重點實驗室，湖北 武漢 430074）

超臨界二氧化碳（S-CO2）具有穩定、易得、無毒、壓縮性好和流動性強等特點，相比于傳統水循環，基于S-CO2工質的新型動力循環具有體積小、效率高、靈活性好[1]等優勢。目前，基于S-CO2循環的研究在光熱發電[2-3]、火力發電[4-5]、核能[6-7]以及余熱利用[8-9]等方面發展迅速。S-CO2新型動力循環的工質阻力特性會對S-CO2循環效率以及各部件設計產生較大影響，若無法準確計算循環管路和主要部件的阻力，S-CO2系統將無法準確設計并高效可靠運行。因此，S-CO2的流動阻力特性可靠評估是S-CO2技術商業化推進中至關重要的一環。

由于湍流的脈動性和新型動力系統實驗數據缺乏，現有研究尚未形成對S-CO2摩擦阻力特性的統一理論。Kirillov等人[10]認為管截面密度分布是影響非等溫流體阻力特性的關鍵因素，提出了壁面密度（依靠內壁溫計算的流體密度）和主流密度比修正項來對等溫流體公式進行修正。Fang等人[11]利用現有文獻中的S-CO2以及其他流體數據，采用大量函數類型進行回歸測試，提出了壁面黏度與主流黏度比修正項以及膜密度與擬臨界點密度比的指數修正項。張海松等[12]研究了S-CO2在擬臨界點附近的阻力系數峰值現象，認為傳熱惡化會對流動阻力特性造成負面影響，并采用超臨界沸騰數SBO和普朗特數比作為修正項計算S-CO2摩擦阻力。然而，上述各修正項可能在一定的參數范圍內適用，而在其他參數范圍將不再是決定性因素，這給S-CO2流動阻力計算的工程化應用造成困難，設計人員難以根據需要采取合理修正項。

基于此，人工神經網絡算法被認為是實現S-CO2阻力預測的一種有效途徑。該方法不考慮流動傳熱的基本理論，而是通過大量訓練，學習數據之間的特征和規律，并基于該規律預測其他數據的特征。這種方式的泛化能力強，不再局限于某一特定公式，可以預測更廣參數范圍的阻力特征。同時，采用神經網絡可以直接將測量值例如溫度、壓力等作為輸入，阻力作為輸出，不再求解不同狀態下的特征參數例如雷諾數、普朗特數等，有利于工程化實現。Ma等人[13]收集了14篇文獻中的超臨界水傳熱數據，利用BP（誤差反向傳播）神經網絡進行了數據訓練以預測其對流換熱系數，結果表明其均方差為4.13%，相關系數R為0.971 7。Balcilar等人[14]基于垂直光滑銅管中的R134a冷卻數據，訓練了多個神經網絡模型以預測其傳熱和壓降特性，預測誤差可保證在5%以內。Shadloo等人[15]采用多層感知器神經網絡對長管道中的氣體/非牛頓流體兩相流的阻力特性進行了研究，預測值與實驗數據的偏差在大多數情況下不超過5%。現有研究中針對S-CO2阻力特性進行神經網絡預測的相關文獻相對較少，神經網絡對于其流動阻力的預測能力有待驗證。此外，現有非絕熱S-CO2實驗數據[12,16-18]主要集中在CO2臨界點（7.38 MPa、31 ℃）[19-20]附近，且質量流速較低（<1 500 kg/(m2·s)）、流體溫度較低（<200 ℃），無法滿足光熱、火電以及核能等領域的實際需求。

聚焦于解決S-CO2發電領域中S-CO2流動摩擦阻力特性的預測問題，本文對更寬參數范圍（750~ 2 200 kg/(m2·s)、10~20 MPa、200~340 kW/m2、60~500 ℃）的S-CO2阻力特性進行實驗測量，分析壓力、質量流速以及熱流密度對摩擦阻力的影響，并根據實驗數據，建立貝葉斯正則化BP神經網絡模型，以實現對S-CO2摩擦阻力特性的準確預測。

1 實驗系統

1.1 S-CO2循環系統

S-CO2流動阻力特性研究實驗系統如圖1所示。該系統首先通過加壓柱塞泵將儲罐中的低壓CO2（純度99.99%）增壓至S-CO2循環回路。S-CO2循環回路由循環泵克服流動阻力，S-CO2依次流過質量流量計、回熱器冷端、預熱器、實驗段、回熱器熱端以及冷卻器，最終再次進入循環泵形成完整循環。冷卻器采用與常規火電站相同的水冷循環實現對S-CO2的冷卻，以保證進入循環泵的S-CO2在50 ℃以內。

實驗時，通過加壓柱塞泵和泄壓閥配合實現對實驗循環壓力的調節；通過控制循環泵頻率實現對實驗循環流量的調節；通過控制預熱器和實驗段的加熱功率實現對實驗段進口溫度以及實驗段熱流密度的調節。系統設計參數為：最大工質壓力30 MPa，最大加熱功率150 kW，最大流量1 800 L/h。

1.2 實驗段

實驗段設計如圖1所示。實驗段管材為Inconel 625合金鋼，管內徑為12 mm，外徑為20 mm，在2個電極板之間的加熱段長度為1 870 mm。

除加熱段外，兩端各留有1 000 mm的過渡段，用于穩定加熱段的流動。實驗段采用電阻發熱原理，將低電壓大電流電源連接到實驗段兩端，即2個電極板處，鋼管本身的電阻使其產生均勻的加熱熱流。同時，實驗段管外壁包裹了保溫層以降低散熱損失。加熱段沿管長方向均勻布置有28個熱電偶，用于測量管外壁溫度。實驗段的入口和出口分別安裝有測量工質溫度的熱電偶和測量工質壓力及壓差的取壓孔。實驗段的流向為垂直向上，與重力方向相反。

1.3 數據采集

實驗中，工質溫度、壁溫測量分別采用Omega 工程公司生產的K型鎧裝熱電偶、K型熱電偶絲，測溫精度為0.75%。工質壓力和壓差測量采用北京俄華通儀表技術有限公司生產的壓力、壓差變送器，測壓精度為0.25%F.S。工質流量測量采用北京首科實華自動化設備有限公司生產的質量流量計，測量精度0.2%。通過模擬信號將上述數據傳輸至由研華科技（中國）有限公司生產的數據采集卡，并轉換為數字信號傳入電腦。該數據采集卡對模擬信號的識別精度為0.1%，數據采集間隔為2 s。當實驗循環各關鍵參數穩定后，開始進行數據采集，選取穩定狀態下的60 s數據取平均后作為該次測量的最終數據，以避免波動誤差。

2 數據處理與神經網絡模型

2.1 數據處理

2.1.1 阻力計算

由于實驗管段豎直，因此實驗系統中測量的壓差包括摩擦阻力ΔPf、重力壓降ΔPg和加速壓降ΔPa。摩擦阻力ΔPf的測量值可由式(1)計算。

式中：ΔP為總壓降，kPa。

摩擦阻力ΔPf還可通過理論計算得到：

式中：f為摩擦阻力系數；L為管長，m；D為管內徑，m；G為流體的質量流速，kg/(m2·s)；ρ為流體密度，kg/m3。

對于重力壓降ΔPg，由于S-CO2的密度變化在實驗參數范圍內是非線性的，因此直接通過簡單平均處理存在較大誤差，本文采用目前認可度較高的Ornatskiy公式[21]來計算：

式中：h為比焓，kJ/kg；g為重力加速度，m/s2。

對于加速壓降ΔPa，采用式(4)計算：

流體比焓、密度等物性參數通過NIST REFPROP 9.1軟件查詢獲取。實驗段沿長度方向不同位置z的主流比焓為：

式中：hb,z為位置z的主流比焓，kJ/kg；z為沿長度方向距實驗段入口的距離，m。

實驗段的加熱熱流密度指內壁面傳遞給流體的熱流密度：

式中：qwall為內壁面熱流密度，kW/m2。

由于實驗只能測量管外壁溫，因此管內壁溫需要通過間接計算得到。假定管周向和軸向的導熱是均勻的，可將本問題簡化為內熱源一維穩態導熱，最終內壁溫計算式如下：

式中：t為溫度，℃；λ為金屬導熱系數，kW/(m·℃)；d為管外徑，m。

2.1.2 實驗參數與不確定度

實驗主要從質量流速、流體壓力以及熱流密度3個方面來探索摩擦阻力的差異。此外，流體溫度和壁面溫度也是測量參數。其中質量流速、熱流密度和摩擦阻力等屬于間接測量值，壓力和溫度等屬于直接測量值。考慮到數據采集卡存在識別精度，即使是直接測量值仍需做不確定度計算。表1給出了各參數在實驗時的數值范圍以及不確定度。

表1 實驗參數范圍與不確定度 Tab.1 Range and uncertainty of experimental parameters

各參數不確定度采用式(8)[16]計算：

2.2 神經網絡模型

BP神經網絡具有信號的前向傳播和誤差的反向傳播2個過程。當輸入信號經過隱含層、輸出層獲得的輸出信號與期望輸出不符時，將輸出誤差反向從輸出層逐層傳回輸入層，每層每個神經元所獲得的誤差數據將作為調整其自身權值的依據，最終使得誤差逐漸下降至目標范圍或不再下降，訓練完成。

本文BP神經網絡模型采用MATLAB編程實現，設計了2層前饋神經網絡結構（圖2），可以實現對任意“一個有限空間到另一個有限空間的連續映射函數”的擬合[22]。更高的層數雖然理論上擬合能力更強，但存在過擬合以及訓練難度增加的問題。隱含層神經元采用tansig非線性激活函數，輸出層神經元采用purelin線性激活函數：

在常規BP神經網絡中，一般采用神經網絡響應的均方差作為目標函數，通過最優化算法實現均方差的逐漸減小，當均方差小于所設定的范圍或不再繼續下降時，認為神經網絡已完成訓練。然而，該方式易出現過擬合現象，使得神經網絡泛化能力較差，無法達到較好的預測效果。與常規BP神經網絡不同，貝葉斯正則化算法將目標函數修改為：

式中：α和β為正則化系數；Ew為神經網絡中所有權值的均方差；Ed為常規神經網絡的目標函數。

設計α和β的意義在于平衡神經網絡擬合效果與泛化能力之間的矛盾。當α過大時，訓練具有更好的泛化能力，但擬合精度不高，可能忽略了較多的數據特征；當β過大時，對訓練集的響應均方差最小，但容易導致過擬合而對非訓練集無法實現較好的預測效果。在貝葉斯理論中，假定訓練集與權重集的先驗概率服從高斯分布，最優的權值應具有最大后驗概率，而最大后驗概率在上述條件下可以等價于最小化F，求得F最小時的α和β為：

式中：γ為有效網絡參數的數量；n為數據個數。該算法更適合訓練具有噪聲和數據量少的數據集[23]。由于阻力實驗各參數均存在不確定度，且數據集規模對于可訓練億級數據的神經網絡來說相對較小，因此采用貝葉斯正則化BP神經網絡作為阻力預測模型，以提高模型泛化能力。

在模型中，將流體壓力、質量流速、熱流密度和主流溫度4個參數作為輸入參數，摩擦阻力系數作為輸出參數。為提高模型預測精度，對輸入參數和輸出參數進行了歸一化處理，歸一化后各參數范圍均在0~1。歸一化方法如下：

實驗總共包含194組數據集，模型訓練集設定為165組，測試集為29組（隨機選取）。隱含層神經元個數的選取采用經驗公式和重復訓練測試聯合完成，常用公式為[24]：

式中：Nhidden為隱含層神經元個數；Nin、Nout分別為輸入、輸出參數個數；X為經驗系數，一般取1~10。本模型中，輸入參數有4個，輸出參數為1個，因此推薦的隱含層神經元個數為3~12。為選擇最佳隱含層神經元個數，進行重復訓練測試，每種神經元個數方案重復訓練100次（每次的測試集是隨機的）并對100次訓練的性能指標取平均作為最終值來對比不同神經元個數的效果。性能評價指標為Pearson相關系數R：

式中：x、y分別為模型計算值和實驗測量值；n為數據個數。相關系數R的取值范圍為0~1，其值越高，表明模型預測效果越好。

圖3顯示了不同隱含層神經元個數對相關系數R的影響。由圖3可以發現，初始時隨著隱含層神經元個數的增加，R逐漸增加，但當隱含層神經元個數超過11時，R幾乎不再增加。因此，選取隱含層神經元個數為11，以確保擬合精度并避免過擬合。在該隱含層神經元個數下的100次重復訓練中，選取R最高的訓練神經網絡作為后續計算分析使用。

3 結果與討論

3.1 貝葉斯正則化BP神經網絡與常規BP神經網絡泛化能力對比

為驗證貝葉斯正則化BP神經網絡在預測S-CO2流動阻力時的泛化優勢，采用新測試數據集（不在原165組訓練集和29組測試集之內）進行性能測試。新測試數據集的參數范圍見表2。

表2 新測試數據集的參數范圍 Tab.2 Parameter range of new test data set

新測試數據集對于模型屬于未知數據集，總共35組，其中21組為文獻數據，剩余14組為本實驗系統提供。此外，有27組（21組文獻+6組本實驗系統）為超過原始194組數據集參數范圍限制的數據（即超限點）。

圖4給出了貝葉斯正則化BP神經網絡與常規BP神經網絡對新測試數據集的預測結果對比。常規BP神經網絡除未使用貝葉斯正則化算法外，其余設定以及訓練均與貝葉斯正則化BP神經網絡相同。圖4中E±10%表示與實驗值誤差在10%以內的預測值個數占總個數的百分比，E±30%表示與實驗值誤差在30%以內的預測值個數占總個數的百分比。

由圖4可以看出：貝葉斯正則化BP神經網絡的相關系數、E±10%和E±30%等指標均優于常規神經網絡；對于超限點，貝葉斯正則化BP神經網絡同樣具有較好的預測性能，而常規BP網絡的誤差則相對較大。綜上，貝葉斯正則化BP神經網絡相比常規BP神經網絡具有更好的泛化能力，可滿足工業應用要求。

3.2 貝葉斯正則化BP神經網絡模型與典型摩擦 系數公式對比

在以往的研究中，非等溫流體（管壁存在加熱或冷卻）經驗公式主要是在等溫流體公式的基礎上進行修正，考慮管截面方向的物性差異。Filonenko公式[25]是目前應用最廣的等溫流體公式，在其后的很多研究者基于該公式進一步修正，獲得非等溫條件下的經驗公式。Kirillov等人[10]給出了基于密度比修正的非等溫流體公式，張海松等[12]給出了基于超臨界沸騰數SBO和普朗特數比的非等溫流體公式。Churchill[26]和Haaland[27]給出了另外2種常見的等溫流體公式。為驗證BP神經網絡模型的實際預測效果，將其計算結果與表3中的5個典型摩擦系數公式計算結果進行對比，如圖5所示。

表3 典型摩擦系數公式 Tab.3 Typical friction coefficient formula

由圖5可知：當摩擦阻力較小時，各公式的預測效果均較好，但當摩擦阻力逐步增加后，預測值和實驗值差距開始變大；在本實驗數據范圍內，3個等溫流體公式的預測值基本一致，3個評價指標值也基本相同；5個典型公式計算值與摩擦阻力實驗值的相關系數均在0.98~0.99，屬于極高相關；在5個典型公式中，Kirillov等人的非等溫流體公式預測效果較好，E±10%=28.35%，E±30%=86.08%；與5個典型公式相比，BP神經網絡模型的預測效果最好，相關系數可達0.998 6，E±10%=98.97%，E±30%=100%。

3.3 不同邊界條件對S-CO2摩擦阻力的影響

圖6給出了不同質量流速G（765.56、1 445.56、2160.71 kg/(m2·s)）下S-CO2摩擦阻力實驗值以及預測值。由圖6可知，隨著質量流速G的增加，S-CO2摩擦阻力增加，且摩擦阻力-主流溫度曲線的斜率也逐漸增大。這是由于摩擦阻力與質量流速的平方G2呈正比關系。溫度變化主要影響摩擦阻力系數和密度，將主流溫度作為自變量時，其斜率主要受G2（正相關）、密度（負相關）以及摩擦阻力系數（正相關）的影響。因此，G增加，其斜率也逐漸增大。

圖7給出了不同熱流密度（208.56、279.14、335.10 kW/m2）下S-CO2摩擦阻力實驗值以及預測值。由圖7可知，隨熱流密度的增加，摩擦阻力下降。熱流密度主要影響了管截面上的流體物性分布，即對于同一主流溫度，熱流密度主要影響摩擦阻力系數。隨著熱流密度的增加，根據實驗現象，其摩擦阻力系數減小。因此，在BP神經網絡模型中，考慮熱流密度的影響是必要的。

圖8給出了不同流體壓力（10.47、15.53、19.52 MPa）下S-CO2摩擦阻力實驗值以及預測值。由圖8可以看出，隨著壓力的增加，摩擦阻力逐漸減小，且該現象在較低壓力下更明顯。根據摩擦阻力公式，當管長L、管徑D和質量流速G不變時，摩擦阻力主要受摩擦阻力系數f（正相關）和密度ρ（負相關）影響。隨壓力的增加，S-CO2密度增加，黏度增加。密度增加會使S-CO2體積流量下降，流速下降，摩擦阻力下降；但黏度增加會導致流體間切應力增加，流體間能量損耗增大，摩擦阻力系數增加。綜合來看，密度影響更關鍵（可能的原因是黏度對摩擦阻力系數的影響并非線性正比例關系且密度同樣影響摩擦阻力系數），因此摩擦阻力下降。

4 結 論

1）S-CO2摩擦阻力隨熱流密度的增加而減小，即其摩擦阻力系數隨熱流密度的增加而減小，這有利于S-CO2循環發電的大規模工程應用。S-CO2摩擦阻力隨流體壓力的增加而減小，主要原因在于壓力增加導致流體密度增加，而黏度增加導致的阻力增加不占主導地位。S-CO2摩擦阻力隨質量流速的變化情況符合傳統規律。

2）完成數據訓練的貝葉斯正則化BP神經網絡模型能夠有效預測S-CO2的摩擦阻力，該方法泛化能力強、擬合度高且可以避免求解不同狀態下的SCO2特征參數，簡化了工程應用，為S-CO2阻力預測提供了新的思路。