鉛鉍自然循環熱工水力分析程序開發

程道喜，杜開文，于 松，齊曉光，翟偉明，周 平，馬曉瑤，張學鋒

(中國原子能科學研究院，北京 102413)

鉛基快堆(LFR)作為第4代核能系統選定的最具有發展前景的反應堆堆型之一，具有良好的增殖核燃料和嬗變核廢料的潛力，并擁有突出的經濟性和安全性。鉛或鉛鉍作為反應堆冷卻劑，擁有良好的自然循環特性，對于鉛基反應堆事故工況下的余熱導出，可以完全依靠一回路系統冷卻劑的自然循環來實現；甚至小型鉛基反應堆可以設計成純自然循環反應堆。因此鉛或鉛鉍的自然循環特性對于鉛基反應堆的正常運行和事故工況余熱導出均具有非常重要的意義。

目前針對鉛鉍自然循環的分析，國內外開展了基于已有程序的適用性評價以及程序開發。如瑞典皇家理工學院Ma等[1-2]開展了基于鉛冷熱工水力回路TALL的自然循環特性分析實驗，模擬了瞬態事故，包括失熱阱事故，一、二次失流事故，超功率以及啟動、停止的瞬態過程等，并利用TRAC/AAA程序等進行了相應的數值模擬分析，驗證了程序對于鉛鉍系統的模擬能力；意大利ENEA的Agostini等[3]在CHEOPE裝置上開展了自然循環實驗并利用RELAP5 MOD3.2開展了相應的分析，取得了較好的結果。Brogohain等[4]開發了鉛鉍自然循環程序，并利用矩形回路自然循環實驗結果進行了初步驗證。Zhao等[5]開發了純自然循環鉛冷反應堆的穩態熱工水力分析程序，并與MARS程序開展了初步對比驗證；Fiori和Zhou[6]驗證了RELAP5/SCDAP對于鉛鉍熱工水力行為的預測能力。陸道綱等[7]基于MATLAB Simulink編制了含有散熱損失模型的鉛鉍自然循環回路計算程序，并利用實驗結果進行了驗證，結果表明鉛鉍回路散熱會對分析結果產生較大的影響。

目前，對于鉛或鉛鉍作為冷卻劑的自然循環特性尚未得到充分認識，上述針對鉛鉍系統的分析程序最初均是針對水等介質開發，對于鉛鉍系統的適用性仍需進行深入研究及驗證。為更好支持鉛鉍反應堆自然循環熱工水力分析工作，本文在研究國內外快堆系統分析程序模型的基礎上，建立堆芯、熱交換器等一回路關鍵部件的熱工水力模型，開發集成鉛鉍工質物性參數及流動換熱關系式的自然循環瞬態分析程序，并基于鉛鉍回路內自然循環實驗數據對程序進行初步驗證。

1 數學物理模型

本程序主要開展一回路系統的自然循環熱工水力分析，不包括中子模塊，并且將二回路視為邊界條件，給定換熱器二次側入口流量及溫度。

1.1 一回路水力模型

閉合回路內自然循環的水力計算，需在整個冷卻劑流動的閉合回路上對動量方程進行積分。每個控制單元動量方程為：

(1)

其中：Δz為控制單元長度；a為流動面積；W為質量流量；t為時間；ρ為流體密度；p為壓力；g為重力加速度；θ為方向角度；f為摩擦阻力系數；De為當量直徑；K為局部阻力系數；i為控制單元號。在整個冷卻劑閉合環路上積分可得：

(2)

其中：j為環路上設備或部件號；ji代表第j個部件的第i個控制單元；Δppump為泵壓頭。

程序中泵模型采用離心泵。泵壓頭可通過額定轉速下的泵特性曲線H0=f(W0)，結合泵轉速與揚程、流量的關系計算得到，如式(4)所示。

(3)

(4)

其中：H為泵揚程；n為轉速；下標0代表額定參數。

1.2 堆芯熱工水力模型

堆芯采用平行多通道模型，每個通道代表1類組件。燃料元件內熱量通過導熱方式進行傳遞，包殼與組件內冷卻劑之間通過對流換熱進行熱量傳遞。柱坐標下，燃料元件導熱方程為：

(5)

堆芯組件內冷卻劑能量方程為：

(6)

組件盒壁的能量方程為：

(7)

其中：下標b、f、s分別代表燃料元件、組件盒內冷卻劑以及組件盒壁；T為溫度；cb為燃料或包殼的比熱；cs為組件盒壁的比熱；r為徑向坐標；z為軸向坐標；λb為燃料或包殼導熱系數；λf為冷卻劑導熱系數；qbv為體積釋熱率；h為冷卻劑比焓；hfb、afb分別為元件棒與冷卻劑的換熱系數和換熱面積；hsf、asf分別為組件盒內壁與冷卻劑的換熱系數和換熱面積；Φs為組件盒壁與外部環境的熱交換量，即：

(8)

其中：hsa、asa分別為組件盒壁與堆內冷卻劑的換熱系數和換熱面積；R為組件盒壁外部與環境大氣的熱阻，根據外壁空氣的自然對流散熱來計算。

1.3 熱交換器模型

圖1 熱交換器模型節點劃分Fig.1 Nodal diagram for heat exchanger

熱交換器為管殼式換熱器，可采用單管模型來進行模擬。熱交換器徑向上可劃分為4個區域：一次側流體、換熱管、二次側流體以及換熱器殼體。熱交換器模型節點劃分如圖1所示。考慮對流換熱以及流體軸向導熱，忽略換熱管以及殼體軸向導熱，分別列出4個區域的控制方程：

一次側流體：

(9)

二次側流體：

(10)

(11)

殼體：

(12)

其中：下標p、se、t、sh分別代表一次側流體、二次側流體、換熱管以及殼體；ct為換熱管壁固體比熱；csh為殼體固體比熱；mp為一次側流體質量流量；mse為二次側流體質量流量；a對于流體為流通面積，對于換熱管或殼體為橫截面積；hpt、apt分別為一次側流體與換熱管換熱系數及面積；hst、ast分別為二次側流體與換熱管之間的換熱系數及面積；hss、ass分別為二次側流體與殼體之間的換熱系數及面積；λ為導熱系數；Φ為殼體與外部的熱交換量，對于“浸泡”在冷卻劑內部的熱交換器，如池式堆中間熱交換器(IHX)等，Φ為殼體與堆池內冷卻劑的熱交換量，對于外部熱交換器，Φ為殼體(添加保溫層)與環境的熱交換器量，表達式與式(8)相同。對于圓形的熱交換器殼體，殼體外部與環境大氣的熱阻R的計算方式為：

(13)

其中：dsh1、dsh2、λsh分別為殼體內徑、外徑和導熱系數；dh1、dh2、λsh分別為保溫層(若有)內徑、外徑和導熱系數；ha為保溫層外表面與大氣的自然對流換熱系數。

1.4 管道模型

管道模型中考慮流體軸向導熱以及管道外壁與周圍介質的熱交換，忽略管壁軸向導熱。管道外壁與周圍介質的熱交換同樣采用式(8)中的方法。

1.5 流動及換熱模型

本程序中采用的鉛鉍換熱關系式列于表1，圓管及環隙流動阻力關系式以及典型結構局部阻力采用Idelchik阻力手冊關系式[8]，棒束阻力采用Cheng-Todreas關系式[9]。

表1 鉛鉍換熱關系式Table 1 Lead-bismuth heat transfer correlation

2 求解方法及程序結構

2.1 求解方法

為保證穩定性，動量方程及各部件的能量方程采用半隱格式進行離散。考慮到動量方程的求解與整個系統溫度分布有關，能量方程的求解與質量流量有關，需采用合適的方法來對動量方程和能量方程解耦。在自然循環建立到穩定的過程中，回路內冷卻劑的熱工水力參數變化較為緩慢，同時在研究的溫度范圍內，各物性參數變化率不大，因此對于此類緩慢變化的物理過程，可通過先求解動量方程，再利用求解動量方程的結果來求解能量方程，并循環迭代直至滿足收斂準則的方法來獲得最終解。本程序中求解步驟如下：

1) 利用前一時間步的溫度分布、流量求解動量方程，得到環路內冷卻劑的流量；

2) 根據求解動量方程獲得的質量流量，從堆芯開始，依次求解各部件的能量方程，獲得溫度分布；

3) 根據新的溫度分布重新求解動量方程，得到冷卻劑流量，并再次依次求解能量方程，循環迭代多次直至滿足收斂標準(溫度、流量相對變化小于限值)，最終獲得的當前時間步的各部件溫度分布及流量等；

目前地方高校科研成果專利保護和轉化不足的主要原因之一就是資金不足。由于專利轉化的成本高，特別是我國高等院校資金主要源于國家撥款，國家對科研的支持是當前高校科研的主要動力。缺乏良好的資本市場和金融環境，嚴重阻礙了科技成果轉化的速度與水平。相比之下，日本不僅在全國各類大學中建立知識產權管理本部，并給予資金上的大力支持，同時成立了區域知識產權戰略總部，在產學研合作開發項目給予資金上的支持。 這些資金扶持的舉措積極推進了大學的創新能力。缺乏風險投資、足量投入，科技成果轉化市場的有序良性發展難以得到保障和維系。

4) 對于各時刻重復1～3過程直至結束時刻。瞬態計算求解流程如圖2所示。

圖2 瞬態計算求解流程Fig.2 Flowchart of transient solving

2.2 程序構架

本程序采用模塊化結構，包括物性模塊、輸入輸出模塊、部件模塊、關系式模塊、求解模塊等，程序結構如圖3所示。程序采用Fortran95語言編制。

圖3 程序結構Fig.3 Structure of code

3 程序初步驗證

圖4 鉛鉍自然循環實驗裝置Fig.4 Experimental device for natural circulation of lead-bismuth

本程序采用鉛鉍回路內自然循環實驗數據進行初步驗證。中國原子能科學研究院開展了鉛鉍自然循環實驗研究工作，獲得了不同功率啟動自然循環以及功率臺階變化對于鉛鉍自然循環影響的實驗數據。該實驗裝置主要由矩形閉合鉛鉍回路以及導熱油冷卻回路等組成，如圖4所示。模擬組件采用電加熱元件模擬堆芯燃料元件發熱，熱交換器為列管式換熱器，一次側(管側)為鉛鉍，二次側(殼側)為導熱油S750(氫化三聯苯)，最終熱阱為空冷器。

回路控制單元劃分如圖5所示。加熱器徑向控制體劃分時，電熱元件內部劃分為3層(最內側為發熱絲，其余兩層為MgO層)，包殼1層。加熱器軸向加熱區和非加熱區通過輸入參數中給定加熱絲軸向體積釋熱率分布來實現模擬；熱交換器徑向分為一次側、二次側、換熱管以及外壁。對于所有回路設備及管道外壁面散熱考慮殼體、管道與外部空氣通過自然對流方式來模擬，自然對流的傳熱關系式采用文獻[14]中的關系式。熱交換器二次側導熱油物性見生產商物性參數表[15]。

3.1 自然循環啟動瞬態分析

在鉛鉍自然循環啟動實驗中，鉛鉍回路內工質從初始靜止等溫狀態開始，瞬時啟動模擬組件電加熱和熱交換器二次側導熱油循環泵，維持電加熱功率、二次側流量以及二次側入口溫度不變。圖6為自然循環啟動過程瞬態溫度和流量的變化。

從圖6可知，在電加熱功率投入后，回路內鉛鉍自然循環流量迅速上升，100 s左右達到峰值流量，同時，電加熱組件出口溫度也同樣迅速上升達到峰值，溫度峰值時間稍早于流量峰值。

圖5 回路控制單元劃分Fig.5 Discretization of loop division of control units of loop

圖6 自然循環啟動過程瞬態溫度和流量的變化Fig.6 Transient temperature and flow rate in natural circulation startup

在相同的二次側邊界條件下，電加熱功率越大，組件峰值溫度越高，最后形成的穩定自然循環流量也越大，冷熱段溫差也越大。由于自然循環流動較弱，瞬態過程中，電加熱組件出口流體到達熱交換器入口需要一定時間，因此熱交換器入口溫度的響應較組件出口要相對滯后；同樣，熱交換器出口流體到達電加熱組件入口也需要一定時間，電加熱組件入口溫度下降較熱交換器出口滯后。在散熱較小的情況下，當穩定的自然循環建立后，冷卻劑溫度分布也逐步達到穩態，最終組件出口與換熱器進口，組件入口與換熱器出口之間的溫度差異變得很小。

從圖6可知，本程序計算結果能較為正確地反映鉛鉍自然循環回路內關鍵位置冷卻劑溫度的瞬態變化規律及相對關系，溫度計算值與實驗值偏差小于5%，流量計算值與實驗值的偏差在自然循環建立的最初階段稍大，其余時刻均小于10%。在實驗中，由于不可避免的熱損失以及熱交換器向二次側的導熱，回路內初始存在溫度差，從而使得自然循環建立速度變化。這一初始溫度差在計算中不能準確反映，導致計算結果在自然循環建立階段的偏差稍大。

3.2 功率臺階變化影響分析

功率臺階變化影響實驗中，鉛鉍工質在一定電加熱功率下建立穩定的自然循環后，不斷改變電加熱功率臺階，鉛鉍回路內工質的自然循環流動將隨著功率臺階的不同而發生變化，組件進出口以及熱交換器的溫度也將響應功率的變化。如圖7所示為功率臺階變化瞬態實驗中，組件電加熱功率隨時間的變化，功率臺階分別為5、8和12 kW。

圖7 組件功率臺階變化Fig.7 Step change of heater power

圖8為功率臺階變化時瞬態溫度和自然循環流量響應。隨著功率以臺階形式變化，溫度與流量有相應的響應。當組件電加熱功率逐步上升時，自然循環流量增加，組件出口溫度先于換熱器出口溫度上升；當組件電加熱功率逐步下降時，自然循環流量降低，組件出口溫度又先于換熱器出口溫度降低。對比結果可知，程序計算得到的結果與實驗結果符合較好，能正確預測功率臺階變化時，組件進出、換熱器進出口等關鍵位置的溫度響應情況以及回路內自然循環流量的瞬態變化規律，其中溫度相對偏差小于5%，自然循環流量相對偏差小于10%，進一步驗證了本程序能較好模擬鉛鉍自然循環瞬態過程。

圖8 功率臺階變化時瞬態溫度和流量響應Fig.8 Transient temperature and flow rate responses during step change of heater power

4 結論

1) 本文針對鉛鉍堆一回路系統主要部件，包括堆芯、熱交換器等建立了熱工水力物理模型，并自主開發了適用于鉛鉍自然循環瞬態分析計算的分析程序。

2) 利用鉛鉍回路內自然循環實驗結果開展了程序初步驗證工作。結果表明，計算得到的溫度和流量的變化規律及趨勢與實驗一致，溫度計算結果相對偏差小于5%，自然循環流量計算結果相對偏差小于10%，初步證明了程序物理模型的正確性以及程序計算鉛鉍自然循環瞬態的能力。