非能動停堆組件熱工水力設(shè)計程序開發(fā)與驗證

黎 閆

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

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

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

1 非能動停堆組件特點

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

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

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

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

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

2 物理模型

2.1 控制方程

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

(1)質(zhì)量守恒方程:

(2)動量守恒方程:

其中,F=Ff+Flocal

(3)能量守恒方程:

式中:ρ——密度;

t——時間;

v——速度;

p——壓力;

g——重力加速度;

Ff——沿程摩擦阻力;

Flocal——局部阻力;

q——源項。

考慮到吸收體棒束結(jié)構(gòu)特點,采用二維非穩(wěn)態(tài)導(dǎo)熱模型進行模擬,在柱坐標系中導(dǎo)熱方程如下:

(4)導(dǎo)熱方程

式中:c——比熱容;

T——溫度;

λ——導(dǎo)熱系數(shù);

S——源項。

2.2 經(jīng)驗關(guān)系式

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

(1)沿程阻力系數(shù)

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

(2)局部突縮/突擴

1)突縮

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

2)突擴

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

(3)直角轉(zhuǎn)彎

直角阻力系數(shù)

(4)對流換熱關(guān)系式

本程序中采用了液態(tài)鈉的Aoki換熱公式:

式中:Nu——努塞爾數(shù);

Pe——佩克萊特數(shù)。的關(guān)系式如下:

3 數(shù)值模型

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

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

3.1 空間離散

3.1.1 流體域空間離散

本程序?qū)⑵胀ü艿滥Ｐ椭械牧黧w視為一維流動,對水力件進行一維空間離散,如圖4所示?？刂企w(j)代表了計算空間劃分的最小幾何單位,同一個水力件的網(wǎng)格為等距劃分,控制體j長度為Zj,控制面(j+1/2)規(guī)定了與各節(jié)點相對應(yīng)的控制體的分界面位置,節(jié)點位于控制體的質(zhì)心。

圖4 水力件網(wǎng)格劃分示意圖Fig.4 Meshing of the hydraulic structure

3.1.2 固體域空間離散

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

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

圖5 熱構(gòu)件節(jié)點示意圖Fig.5 The calculation process

3.2 方程離散

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

3.2.1 質(zhì)量方程

對質(zhì)量方程,對控制體j進行差分,可得到控制體密度與速度的關(guān)系如下:

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

3.2.2 動量方程

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

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

3.2.3 能量方程

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

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

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

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

3.2.4 導(dǎo)熱方程

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

其中,

3.3 邊界條件處理

3.3.1 流動邊界條件

(1)流量入口邊界條件

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

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

(2)壓力邊界條件

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

3.3.2 導(dǎo)熱邊界條件

本程序采用附件源項法對邊界節(jié)點進行離散化,即將邊界條件中所規(guī)定的進入或?qū)С鲇嬎銋^(qū)域的熱量作為與邊界相鄰的控制體的當量源項,以節(jié)點內(nèi)側(cè)r方向內(nèi)側(cè)作為邊界為例,對各邊界條件下邊界節(jié)點(i,j)進行離散。

(1)第二類邊界條件

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

(2)共軛換熱

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

對于熱構(gòu)件2則有:

(3)流-固換熱

假設(shè)熱構(gòu)件內(nèi)側(cè)為流體,換熱系數(shù)hf,流體溫度Tf,相對于第二類邊界條件來說,離散公式中只有系數(shù)ai,j和Si,j發(fā)生變化,具體如下:

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

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

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

4 程序驗證

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

4.1 程序?qū)Ρ闰炞C

將本程序的水力計算結(jié)果與Mathematica程序進行對比分析,計算非能動停堆組件內(nèi)關(guān)鍵部件及非能動停堆組件整組件在360℃液態(tài)鈉中的壓降,在計算中保證結(jié)構(gòu)尺寸和介質(zhì)工況的一致性。各關(guān)鍵部件及非能動停堆組件壓降計算結(jié)果如表1所示。臨界流量(即非能動棒剛好懸浮于釋放位時對應(yīng)的流量)是非能動停堆組件的關(guān)鍵指標之一,利用本程序?qū)υ撝笜酥颠M行計算,得到結(jié)果如表2所示。由計算結(jié)果可知,本程序與現(xiàn)有成熟程序的水力計算結(jié)果相對偏差很小,壓降最大相對偏差為6.77%,臨界流量相對偏差為3.37%。

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

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

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

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

4.2 試驗對比驗證

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

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

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

5 結(jié)論

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