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

黎 閆

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

第四代核能系統(tǒng)國際論壇(GIF)將第四代核電廠的開發(fā)目標(biāo)設(shè)定為四個方面:核能的可持續(xù)發(fā)展、提高安全性和可靠性、提高經(jīng)濟(jì)性、防止核擴散[1],對鈉冷快堆的安全性提出了更高的要求,增設(shè)非能動停堆裝置可以有效提高鈉冷快堆的固有安全性[2]。目前國際上已經(jīng)對非能動停堆裝置進(jìn)行了大量研究[3,4],液體懸浮式非能動停堆組件是其中技術(shù)成熟度相對較高,且有實堆經(jīng)驗的一種[5,6],主要用于對鈉冷快堆發(fā)生無保護(hù)失流事故的緩解[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é)流部件。冷卻劑從非能動停堆組件管腳入口進(jìn)入到組件內(nèi)部后主要分為兩部分,一部分冷卻劑由非能動棒上的小孔流入到非能動棒內(nèi)對吸收體棒束進(jìn)行冷卻,經(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)熱模型對不同溫度液鈉下非能動停堆組件水力特性進(jìn)行計算是十分必要的。

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

2 物理模型

2.1 控制方程

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

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

(2)動量守恒方程:

其中,F=Ff+Flocal

(3)能量守恒方程:

式中:ρ——密度;

t——時間;

v——速度;

p——壓力;

g——重力加速度;

Ff——沿程摩擦阻力;

Flocal——局部阻力;

q——源項。

考慮到吸收體棒束結(jié)構(gòu)特點,采用二維非穩(wěn)態(tài)導(dǎo)熱模型進(jìn)行模擬,在柱坐標(biā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的流道進(jìn)入橫截面積為A1的流道,局部阻力系數(shù)為:

2)突擴

對于流體從橫截面積為A1A1的流道進(jìn)入橫截面積為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視為一維流動,對水力件進(jìn)行一維空間離散,如圖4所示。控制體(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 固體域空間離散

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

3.2.1 質(zhì)量方程

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

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

3.2.2 動量方程

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

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

3.2.3 能量方程

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

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

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

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

3.2.4 導(dǎo)熱方程

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

其中,

3.3 邊界條件處理

3.3.1 流動邊界條件

(1)流量入口邊界條件

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

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

(2)壓力邊界條件

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

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

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

(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 程序驗證

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

4.1 程序?qū)Ρ闰炞C

將本程序的水力計算結(jié)果與Mathematica程序進(jìn)行對比分析,計算非能動停堆組件內(nèi)關(guān)鍵部件及非能動停堆組件整組件在360℃液態(tài)鈉中的壓降,在計算中保證結(jié)構(gòu)尺寸和介質(zhì)工況的一致性。各關(guān)鍵部件及非能動停堆組件壓降計算結(jié)果如表1所示。臨界流量(即非能動棒剛好懸浮于釋放位時對應(yīng)的流量)是非能動停堆組件的關(guān)鍵指標(biāo)之一,利用本程序?qū)υ撝笜?biāo)值進(jìn)行計算,得到結(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

同時,對本程序的熱工計算模塊進(jìn)行驗證,考慮非能動停堆組件在反應(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 試驗對比驗證

為了對本程序進(jìn)行驗證,除了與現(xiàn)有相對成熟的計算程序的計算結(jié)果對比驗證,同時與非能動停堆組件的水力試驗結(jié)果進(jìn)行對比。非能動停堆組件水力試驗采用非能動停堆組件的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)熱過程進(jìn)行求解,基于此完成了非能動停堆組件熱工水力設(shè)計程序的開發(fā)。通過與現(xiàn)有成熟程序的計算結(jié)果和非能動停堆組件的水力試驗結(jié)果進(jìn)行對比,本程序計算非能動停堆組件熱工水力結(jié)果的相對偏差不大于7%,滿足工程設(shè)計計算的要求,同時大大縮短了建模時間,驗證了程序的準(zhǔn)確性和高效性,可為后續(xù)非能動停堆組件熱工水力設(shè)計工作提供有效的計算工具。