非能動停堆組件熱工水力設計程序開發與驗證

黎 閆

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

第四代核能系統國際論壇(GIF)將第四代核電廠的開發目標設定為四個方面:核能的可持續發展、提高安全性和可靠性、提高經濟性、防止核擴散[1],對鈉冷快堆的安全性提出了更高的要求,增設非能動停堆裝置可以有效提高鈉冷快堆的固有安全性[2]。目前國際上已經對非能動停堆裝置進行了大量研究[3,4],液體懸浮式非能動停堆組件是其中技術成熟度相對較高,且有實堆經驗的一種[5,6],主要用于對鈉冷快堆發生無保護失流事故的緩解[7,8]。中國正在設計及建造的中國示范快堆也將采用液體懸浮式非能動停堆組件(簡稱非能動停堆組件),非能動停堆組件通過組件自身水力特性實現其懸浮、液力自緊和落棒等功能,非能動停堆組件的熱工水力設計研究是研發設計中的核心。

非能動停堆組件結構復雜,在非能動停堆組件的方案設計過程中需要開展大量的熱工水力計算,因此開發具備快速的瞬態和穩態計算能力,且滿足工程設計精度要求的非能動停堆組件熱工水力設計程序是很有必要的。本文在研究非能動停堆組件結構特點的基礎上,建立了非能動停堆組件的物理模型和數值模型,基于此開發了一維的非能動停堆組件熱工水力設計程序,并完成程序驗證。

1 非能動停堆組件特點

非能動停堆組件主要包括兩部分:可移動的非能動棒(移動體)和固定的組件套筒,非能動棒可以在組件套筒的導向管內上下移動。以中國示范快堆中非能動停堆組件為例,其內部流道如圖1所示。

為滿足非能動停堆組件的各項功能要求,需要在非能動停堆組件內部設置小孔(如②、③、○17等)、節流件(如⑦等)和狹縫(如○13、○16等)等特殊的節流部件。冷卻劑從非能動停堆組件管腳入口進入到組件內部后主要分為兩部分,一部分冷卻劑由非能動棒上的小孔流入到非能動棒內對吸收體棒束進行冷卻,經由節流件從出口流出,另一部分冷卻劑由非能動棒與組件套筒之間的環形縫隙流出,兩者在組件出口處匯合后流出組件,非能動停堆組件內部流道形成了包含分流、匯流、沿程、突縮和突擴等結構的水力管網。

圖1 非能動停堆組件內部流道示意圖Fig.1 Flow channel of the passive shutdown assembly

同時,在反應堆正常運行時,非能動停堆組件內吸收體棒束會產生釋熱,吸收體棒束處截面如圖2所示,由內而外分別為吸收體芯塊、包殼、套管、導向管和外套管,釋熱功率在徑向和軸向上非均勻分布,不同固體結構之間及固體結構與流體之間都存在換熱。非能動停堆組件內冷卻劑經加熱后,可使組件出、入口冷卻劑溫度產生約50℃的溫差,液鈉的密度和運動粘度是影響非能動停堆組件水力特性的關鍵因素,因此通過建立導熱模型對不同溫度液鈉下非能動停堆組件水力特性進行計算是十分必要的。

圖2 吸收體棒束截面圖Fig.2 The cross section of the absorber rod

2 物理模型

2.1 控制方程

將非能動停堆組件內部各流道結構等效為普通管道,采用單相流模型進行模擬,將普通管道模型中的流體視為一維流動,可得到三大基本方程,即質量方程、動量方程和能量方程,如下:

(1)質量守恒方程:

(2)動量守恒方程:

其中,F=Ff+Flocal

(3)能量守恒方程:

式中:ρ——密度;

t——時間;

v——速度;

p——壓力;

g——重力加速度;

Ff——沿程摩擦阻力;

Flocal——局部阻力;

q——源項。

考慮到吸收體棒束結構特點,采用二維非穩態導熱模型進行模擬,在柱坐標系中導熱方程如下:

(4)導熱方程

式中:c——比熱容;

T——溫度;

λ——導熱系數;

S——源項。

2.2 經驗關系式

在動量和能量方程中,需要考慮局部阻力和沿程阻力,其中局部阻力主要包括由管道突縮、突擴引起的局阻及管道發生直角轉彎產生的局阻。本文通過查閱水力手冊[9,10],選取適用于非能動停堆組件內結構的經驗公式,同時需要選取適用于非能動停堆組件的對流換熱關系式,具體如下:

(1)沿程阻力系數

通過查水力手冊,圓形流道紊流狀態(4 000＜Re＜106)的摩擦阻力系數如下:

(2)局部突縮/突擴

1)突縮

對于流體從橫截面積為A2的流道進入橫截面積為A1的流道,局部阻力系數為:

2)突擴

對于流體從橫截面積為A1A1的流道進入橫截面積為A2的流道,局部阻力系數:

(3)直角轉彎

直角阻力系數

(4)對流換熱關系式

本程序中采用了液態鈉的Aoki換熱公式:

式中:Nu——努塞爾數;

Pe——佩克萊特數。的關系式如下:

3 數值模型

根據非能動停堆組件物理模型特點,將計算域分為流體域和固體域,相應的流體域結構為水力件,固體域結構為熱構件,如圖3所示。流體域入口為流量入口邊界條件,出口為壓力邊界條件,熱構件與熱構件之間、熱構件與水力件之間存在換熱,基于此建立非能動停堆組件的數值模型。

圖3 非能動停堆組件計算模型示意圖Fig.3 The calculation model of the passive shutdown assembly

3.1 空間離散

3.1.1 流體域空間離散

本程序將普通管道模型中的流體視為一維流動,對水力件進行一維空間離散,如圖4所示。控制體(j)代表了計算空間劃分的最小幾何單位,同一個水力件的網格為等距劃分,控制體j長度為Zj,控制面(j+1/2)規定了與各節點相對應的控制體的分界面位置,節點位于控制體的質心。

圖4 水力件網格劃分示意圖Fig.4 Meshing of the hydraulic structure

3.1.2 固體域空間離散

在柱坐標中對非能動停堆組件熱構件采用二維離散化處理,如圖5所示,同一熱構件在r和z方向的步長均為等步長,r方向步長為Δr,z方向步長為Δz,體積單元為ΔV。

表3還表明，在不同分位點上，β的估計值存有比較明顯的差異：第一，不同分位點上，對于不同期貨合約下，β存在明顯差異；第二，隨著分位數的增大，不同期貨合約下，β有變小的趨勢，由此說明，隨著原油期貨價格的上漲，期貨價格對現貨價格的引導作用不斷下降。

圖5 熱構件節點示意圖Fig.5 The calculation process

3.2 方程離散

本程序采用交錯網格思想,對三大基本方程用半隱解法進行離散[11],即在控制體處求解質量和能量方程,在接管(即控制面)處求解動量方程,動量守恒方程是以接管為中心的兩個半控制體為單元進行差分的。在交錯網格中,壓力、溫度、密度和能量等參數都定義在控制體的中心,速度則定義在接管上。以普通管道中的控制體j為例對三大基本方程進行離散,流通面積為Aj,當前時刻上標為1,上一時刻上標為0。三大基本方程離散結果如下:

3.2.1 質量方程

對質量方程,對控制體j進行差分,可得到控制體密度與速度的關系如下:

施主元法同樣適用于動量方程和能量方程的離散中。

3.2.2 動量方程

對動量守恒方程,對各接管進行離散處理,可得到接管速度與壓力的關系式。對上游接管j-1/2進行差分,結果如下:

對下游接管j+1/2進行差分,結果如下:

3.2.3 能量方程

對能量方程,對控制體j進行差分如下:

將離散后的質量、動量方程代入能量方程,轉化為只與控制體壓力相關的非線性方程組,可形成關于控制體壓力的非線性矩陣:

其中等式左邊即為Ej,用牛頓迭代法[12]對該矩陣進行求解,則可將關于壓力的非線性矩陣處理成牛頓迭代法的Jacobi方程組:

經計算可得到相應控制體壓力分布值,將得出的壓力值代入質量、動量和狀態方程,即可得到控制體密度、接管處速度、控制體焓值及流體溫度等其他量。

3.2.4 導熱方程

根據節3.1.2中固體域的空間離散,在柱坐標中對導熱方程進行離散,假定時間步長為Δt,體積源項為s,溫度為T,對節點(i,j)做積分,可得到如下的離散方程形式:

其中,

3.3 邊界條件處理

3.3.1 流動邊界條件

(1)流量入口邊界條件

與普通接管的區別在于,入口邊界控制體的流量和溫度恒定,假設接管j+1/2為入口邊界相關接管,則動量方程差分格式變化如下:

接管j-1/2的動量表達形式不變。

(2)壓力邊界條件

對于壓力邊界的處理,在內迭代形成壓力矩陣時,壓力邊界控制體的系數為1,其余系數為0,殘差項為0,以控制體j為壓力邊界控制體為例,關于壓力Jacobi方程組中,有

3.3.2 導熱邊界條件

本程序采用附件源項法對邊界節點進行離散化,即將邊界條件中所規定的進入或導出計算區域的熱量作為與邊界相鄰的控制體的當量源項,以節點內側r方向內側作為邊界為例,對各邊界條件下邊界節點(i,j)進行離散。

(1)第二類邊界條件

根據第2節可知,非能動停堆組件中芯塊結構中心為絕熱邊界條件,即第二類邊界條件,假定熱流密度為q,則離散方程與內部節點具有相同的方程形式如公式(18),系數作如下變化:

(2)共軛換熱

非能動停堆組件中存在熱構件之間的共軛換熱,如包殼與芯塊之間的換熱,假定節點(i,j)內側連接其他熱構件,記為熱構件2,設其節點溫度為T2,導熱系數為λ2,步長為Δr2,則相對于第二類邊界條件來說,節點(i,j)的離散公式中只有系數ai,j和Si,j發生變化,具體如下:

對于熱構件2則有:

(3)流-固換熱

假設熱構件內側為流體,換熱系數hf,流體溫度Tf,相對于第二類邊界條件來說,離散公式中只有系數ai,j和Si,j發生變化,具體如下:

對于流體域的水力件則有:

本程序計算流程圖如圖6所示。

圖6 程序計算流程圖Fig.6 Calculation process

4 程序驗證

為了驗證程序的準確性,主要通過與現有成熟程序的計算結果及非能動停堆組件的水力試驗結果進行對比,完成程序的驗證工作。

4.1 程序對比驗證

將本程序的水力計算結果與Mathematica程序進行對比分析,計算非能動停堆組件內關鍵部件及非能動停堆組件整組件在360℃液態鈉中的壓降,在計算中保證結構尺寸和介質工況的一致性。各關鍵部件及非能動停堆組件壓降計算結果如表1所示。臨界流量(即非能動棒剛好懸浮于釋放位時對應的流量)是非能動停堆組件的關鍵指標之一,利用本程序對該指標值進行計算,得到結果如表2所示。由計算結果可知,本程序與現有成熟程序的水力計算結果相對偏差很小,壓降最大相對偏差為6.77%,臨界流量相對偏差為3.37%。

表1 關鍵部件及整組件壓降計算結果對比Table 1 The pressure drop of key components and the passive shutdown assembly

表2 臨界流量計算結果對比Table 2 The critical flow rate of the passive shutdown assembly

同時,對本程序的熱工計算模塊進行驗證,考慮非能動停堆組件在反應堆正常運行時的釋熱功率約為160.594 k W,其中吸收體棒束段功率約為42.052 k W,非能動停堆組件冷卻劑入口溫度為358℃,流量為4.3 kg/s。非能動停堆組件出口溫度計算結果與CFD模擬計算結果對比如表3所示。

表3 非能動停堆組件出口溫度計算結果對比Table 3 The temperature of the passive shutdown assembly outlet

4.2 試驗對比驗證

為了對本程序進行驗證,除了與現有相對成熟的計算程序的計算結果對比驗證,同時與非能動停堆組件的水力試驗結果進行對比。非能動停堆組件水力試驗采用非能動停堆組件的1∶1模擬件在93℃水回路試驗段中開展,測量不同流量下的組件壓降,同時測量了非能動棒剛好懸浮于釋放位時對應的臨界流量。計算中保持幾何尺寸及介質工況的一致性。非能動停堆組件在不同流量下壓降的試驗結果與計算結果如表4所示,非能動停堆組件臨界流量的計算結果如表5所示。由計算結果可知,本程序的水力計算結果與試驗結果相對偏差很小,壓降的最大相對偏差為6.42%,臨界流量的相對偏差為1.18%。

表4 非能動停堆組件壓降試驗與計算結果對比Table 4 The pressure drop of the passive shutdown assembly

表5 非能動停堆組件臨界流量試驗與計算結果對比Table 5 The temperature of the passive shutdown assembly outlet

5 結論

本文通過分析非能動停堆組件工作原理及結構特點,將非能動停堆組件內部流道等效為一維流動,利用交錯網格思想,采用半隱差分格式和施主元法,聯立三大基本方程求解壓力矩陣,并建立導熱方程對導熱過程進行求解,基于此完成了非能動停堆組件熱工水力設計程序的開發。通過與現有成熟程序的計算結果和非能動停堆組件的水力試驗結果進行對比,本程序計算非能動停堆組件熱工水力結果的相對偏差不大于7%,滿足工程設計計算的要求,同時大大縮短了建模時間,驗證了程序的準確性和高效性,可為后續非能動停堆組件熱工水力設計工作提供有效的計算工具。