林 超，楊紅義，周志偉

(中國(guó)原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

鈉冷快堆作為核能系統(tǒng)國(guó)際論壇(GIF)提出的6種第4代核電站堆型之一，具有核能的可持續(xù)發(fā)展、更高的安全性和可靠性、更高的經(jīng)濟(jì)性、防止核擴(kuò)散等優(yōu)點(diǎn)。同壓水堆相比，鈉冷快堆能將鈾資源的利用率由1%提高到60%以上[1]?；谏鲜鰞?yōu)點(diǎn)，中國(guó)正在大力發(fā)展鈉冷快堆技術(shù)，并計(jì)劃于2022年建成中國(guó)示范快堆。

中國(guó)示范快堆和實(shí)驗(yàn)快堆均屬于池式鈉冷快堆，其堆芯組件按冷卻方式可分為兩種類型：1) 強(qiáng)迫循環(huán)冷卻組件，包括燃料組件、轉(zhuǎn)換區(qū)組件、控制棒組件等，這些組件的發(fā)熱量較大，通過(guò)一次鈉泵提供足夠的冷卻劑流量對(duì)其進(jìn)行冷卻；2) 自然循環(huán)冷卻組件，包括乏燃料組件、鋼組件反射層和碳化硼組件，這些組件的發(fā)熱量較小，所需冷卻劑的流量也較小，無(wú)需泵提供強(qiáng)迫循環(huán)流量，而是依靠其自身發(fā)熱產(chǎn)生溫差驅(qū)動(dòng)流體進(jìn)入組件內(nèi)部對(duì)其進(jìn)行冷卻。

對(duì)于自然循環(huán)冷卻組件，無(wú)法對(duì)組件流量進(jìn)行精確分配，而是依靠組件自身的發(fā)熱進(jìn)行流量的自動(dòng)分配。其流量來(lái)源于堆芯漏流，但漏流不僅可進(jìn)入自然循環(huán)組件內(nèi)部，也能進(jìn)入組件盒間形成盒間流。因此，自然循環(huán)組件熱工分析十分復(fù)雜[2-3]。

目前，對(duì)于已建成的中國(guó)實(shí)驗(yàn)快堆(CEFR)和正在建設(shè)的示范快堆的設(shè)計(jì)，自然循環(huán)冷卻組件的計(jì)算分析采用一維方法。然而，一維的計(jì)算分析方法有著較大的誤差，無(wú)法精確得出自然循環(huán)冷卻組件的流量和熱點(diǎn)溫度。對(duì)于現(xiàn)有的鈉冷快堆子通道分析程序，如SUPERENERGY、COBRA-Ⅳ[4]、MATRA-LMR[5]等，均無(wú)法準(zhǔn)確模擬低雷諾數(shù)(Re)下的交混，并且無(wú)法自動(dòng)計(jì)算每盒自然循環(huán)冷卻組件的流量，因此目前尚無(wú)專門針對(duì)鈉冷快堆自然循環(huán)冷卻組件相關(guān)的子通道分析程序。

本文基于鈉冷快堆自然循環(huán)組件的特點(diǎn)以及子通道方法，開(kāi)發(fā)鈉冷快堆自然循環(huán)組件穩(wěn)態(tài)子通道分析程序，并對(duì)其進(jìn)行驗(yàn)證。

1 程序模型

1.1 單相流穩(wěn)態(tài)子通道守恒方程

程序中子通道方法的主要假設(shè)為：流體是二維的，對(duì)于任意子通道僅考慮軸向和橫向兩個(gè)方向。

對(duì)于鈉冷快堆，冷卻劑溫度距其沸點(diǎn)有較大的裕度，因此無(wú)需考慮相變。單相穩(wěn)態(tài)質(zhì)量守恒方程為：

(1)

式中：下標(biāo)i代表與該子通道相鄰的子通道；m為子通道軸向質(zhì)量流量，kg/s；w為橫流質(zhì)量流量，kg/(m·s)。

穩(wěn)態(tài)能量守恒方程為：

(2)

穩(wěn)態(tài)軸向動(dòng)量守恒方程為：

(3)

式中：u為軸向流速，m/s；p為子通道的壓力，Pa；f為摩擦系數(shù)；CT為動(dòng)量修正常數(shù)；K為軸向流動(dòng)形阻系數(shù)；Dh為子通道水力直徑，m。

穩(wěn)態(tài)橫向動(dòng)量守恒方程為：

(4)

式中：l為間隙控制體等效長(zhǎng)度，m；KG為橫向流動(dòng)形阻系數(shù)；ρ為子通道密度，kg/m3；下標(biāo)i、j代表與該間隙相鄰的子通道編號(hào)。橫流方程針對(duì)的是所有子通道之間的間隙。

在進(jìn)行守恒方程求解時(shí)，采用有限差分法對(duì)空間項(xiàng)進(jìn)行離散。

1.2 軸向壓降模型

對(duì)于本程序，由于需根據(jù)自然循環(huán)組件壓降為組件分配流量，因此軸向壓降模型的選取尤為重要。自然循環(huán)組件Re較低，跨越了層流區(qū)、過(guò)渡區(qū)和紊流區(qū)，因此本程序的軸向壓降關(guān)系式使用詳細(xì)CT模型[6]，該模型覆蓋了全Re范圍，并對(duì)過(guò)渡區(qū)進(jìn)行了合理修正，并且該程序與實(shí)驗(yàn)結(jié)果之間的偏差較小[7-8]，其表達(dá)式如下：

fi=Cfi/Rem

(5)

式中：m=1(湍流)或0.18(層流)；Cf為各類子通道摩擦系數(shù)因子；下標(biāo)i代表不同子通道類型。Cf主要由光棒束的摩擦阻力與繞絲的橫掠阻力疊加而成，其值與當(dāng)?shù)豏e、子通道類型、繞絲螺距、棒束直徑、棒束中心距等參數(shù)有關(guān)。對(duì)于不同類型的子通道，其Cf的計(jì)算方法不同。對(duì)于過(guò)渡區(qū)，摩擦系數(shù)計(jì)算方法如下：

fTran=fLam(1-ψ)1/3(1-ψλ)+fTurψ1/3

(6)

ψ=lg(Re/ReL)/lg(ReT/ReL)

(7)

式中：下標(biāo)Tran代表過(guò)渡區(qū)，Lam代表層流區(qū)，Tur代表湍流區(qū)，L代表層流區(qū)到過(guò)渡區(qū)的轉(zhuǎn)換點(diǎn)，T代表過(guò)渡區(qū)到湍流區(qū)的轉(zhuǎn)換點(diǎn)；λ為擬合常數(shù)，一般取13。

1.3 橫向壓降模型

本文采用gunther關(guān)系式[9]研究棒束區(qū)的軸向壓降關(guān)系。

對(duì)層流Re<200，有：

(8)

對(duì)湍流Re≥200，有：

(9)

式中：Dv為橫流流道的水力直徑；ReDv為橫流雷諾數(shù)；P為相鄰棒中心距。

1.4 換熱模型

換熱模型采用Mikityuk關(guān)系式[10]。Mikityuk在2009年對(duì)國(guó)際上的液態(tài)金屬傳熱關(guān)系式進(jìn)行了分析對(duì)比，并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了重新擬合。相對(duì)于其他換熱關(guān)系式，此擬合關(guān)系式精度最好，該關(guān)系式如下：

Nu=0.047[1-e-3.8(P/D-1)](Pe0.77+250)

(10)

式中：Nu為努塞爾數(shù)；Pe為貝克來(lái)數(shù)；P/D為中心距直徑比。該公式的適用范圍:30

1.5 交混模型

子通道之間的交混采用MIT交混模型[11]。鈉冷快堆內(nèi)自然循環(huán)組件內(nèi)流體的Re較低，最低約為103。MIT交混模型覆蓋的Re范圍較廣，Re最低可達(dá)520[11]。其余交混模型，如Rogers-Tahir、Rower-Angle、Seale等，所覆蓋的Re最低范圍均在104量級(jí)[12]，因此不適用于自然循環(huán)冷卻組件。

根據(jù)MIT交混模型，在低Re自然對(duì)流情況下，交混系數(shù)表達(dá)式為：

Γm=ε1+εM+κυ

(11)

式中，等號(hào)右側(cè)分別代表湍流和繞絲交混系數(shù)、熱羽(溫差)交混系數(shù)、幾何形狀交混系數(shù)。對(duì)于熱羽導(dǎo)致的交混，MIT交混模型給出了一個(gè)狀態(tài)切換函數(shù)，該函數(shù)的數(shù)值在0～1之間，并且隨著理查森數(shù)(Ri)的變化而變化。

1.6 組件間傳熱模型

本文在處理組件間傳熱時(shí)，將組件間的鈉劃分為控制體，以計(jì)算該控制體與其周圍組件邊子通道或角子通道之間的換熱。針對(duì)所有組件的每個(gè)參與組件間傳熱的邊、角子通道，分別求出它們通過(guò)組件間傳熱導(dǎo)致的換熱量，然后更新相應(yīng)子通道的能量方程。組件間控制體劃分如圖1所示[13]。由于快堆盒間流速很低，因此本文在處理盒間傳熱時(shí)將盒間流視為滯止?fàn)顟B(tài)，簡(jiǎn)化為固體控制體，該簡(jiǎn)化方法是保守的。

圖1 組件間控制體劃分Fig.1 Sub-channel meshing for inter-assembly

組件內(nèi)部邊子通道或角子通道與盒間控制體之間的熱阻采用下式進(jìn)行計(jì)算[14]：

R=1/h+dhex/λhex+dNa/2λNa

(12)

式中：h為冷卻劑與盒壁之間的換熱系數(shù)，W/(m2·K)，可通過(guò)Mikityuk關(guān)系式得出；dhex為盒壁厚度，m；λhex為盒壁的熱導(dǎo)率，W/(m·K)；dNa為盒間隙厚度，m；λNa為鈉的熱導(dǎo)率，W/(m·K)。求得盒間子通道與相鄰組件內(nèi)部子通道之間的熱阻后，通過(guò)盒間子通道的能量守恒方程即可求出盒間子通道的溫度以及換熱量，從而得出能量方程中的Q盒間。

1.7 流量分配模型

所有自然循環(huán)組件的盒內(nèi)流和盒間流均為并聯(lián)通道。因此，自然循環(huán)組件盒內(nèi)和盒間的入口、出口均可視為等壓面，因此采用壓差平衡的方法為自然循環(huán)組件及盒間流分配流量。通過(guò)調(diào)整每盒自然循環(huán)冷卻組件的流量及盒間流量，使盒內(nèi)、盒間的出入口壓差均相等，從而完成流量分配。

在計(jì)算入口壓差時(shí)，需分別計(jì)算重力壓差、摩擦阻力及局部阻力。對(duì)于重力壓差，根據(jù)盒內(nèi)、盒間沿軸向溫度分段計(jì)算。盒內(nèi)的摩擦阻力系數(shù)采用1.2節(jié)的軸向壓降模型。對(duì)于盒間流，其Re較低，處于層流狀態(tài)，且其流道結(jié)構(gòu)較為簡(jiǎn)單，因此盒間摩擦阻力系數(shù)ξ通過(guò)下式計(jì)算[15]：

ξ=64/Re

(13)

由于盒內(nèi)、盒間流道較為復(fù)雜，其局部阻力特性無(wú)法通過(guò)查詢阻力手冊(cè)的方法得出，因此需通過(guò)CFD計(jì)算來(lái)獲得，作為本程序的輸入。

獲得所有自然循環(huán)冷卻組件盒內(nèi)流和盒間流入出口壓差后，通過(guò)質(zhì)量流量平均的方法求出平均壓差。如果有任意通道的壓差與平均壓差之間的偏差不滿足收斂準(zhǔn)則，則通過(guò)下式進(jìn)行流量的重新分配：

(14)

式中：Mi為盒間總流量或其中1盒自然循環(huán)組件盒內(nèi)流量，kg/s；M′i為根據(jù)壓差結(jié)果重新分配的流量，kg/s；Δpi為該通道的壓差，Pa；Δpave為所有通道的質(zhì)量流量平均壓差，Pa；α為流量分配加速收斂因子，該因子需根據(jù)具體問(wèn)題來(lái)進(jìn)行具體設(shè)置。

2 計(jì)算流程

本程序在進(jìn)行自然循環(huán)組件計(jì)算時(shí)，分為單組件穩(wěn)態(tài)迭代計(jì)算、組件間換熱迭代計(jì)算、自然循環(huán)冷卻組件流量-盒間流量分配迭代計(jì)算3個(gè)流程?？傮w計(jì)算流程圖如圖2所示。其中對(duì)于第1輪組件間換熱迭代計(jì)算，將組件間換熱量置0。對(duì)于第1輪自然循環(huán)冷卻組件流量-盒間流量分配迭代計(jì)算，需手動(dòng)輸入自然循環(huán)組件盒內(nèi)總流量與盒間流量的比值，同時(shí)按照自然循環(huán)組件功率比例和自然循環(huán)組件總流量給自然循環(huán)組件分配初始計(jì)算流量。

圖2 總體計(jì)算流程圖Fig.2 Flow chart of overall calculation

2.1 單組件穩(wěn)態(tài)迭代計(jì)算

燃料組件的子通道分析屬于三維數(shù)值模擬范疇，但如果采用全場(chǎng)求解的方法將會(huì)導(dǎo)致求解矩陣龐大，所以本文采用逐層計(jì)算的方法，每層迭代收斂后再沿軸向計(jì)算下一層，最后得到整個(gè)組件的熱工水力特性。單組件穩(wěn)態(tài)迭代計(jì)算流程圖如圖3所示。

2.2 組件間換熱迭代計(jì)算

在進(jìn)行組件間換熱迭代計(jì)算時(shí)，根據(jù)上一次單組件穩(wěn)態(tài)迭代計(jì)算得出的組件邊角子通道的溫度和相應(yīng)熱阻，得出組件間的換熱量，并通過(guò)多次迭代，使相鄰組件間的溫差和傳熱量相匹配。組件間換熱迭代計(jì)算流程圖如圖4所示。

2.3 自然循環(huán)冷卻組件流量分配迭代計(jì)算

自然循環(huán)冷卻組件流量分配迭代計(jì)算根據(jù)上一輪迭代的盒內(nèi)、盒間流量和壓降修正各通道的流量，直至所有通道的壓差與質(zhì)量流量平均壓差之間的偏差小于1×10-4。自然循環(huán)冷卻組件流量分配迭代計(jì)算流程圖如圖5所示。

圖3 單組件穩(wěn)態(tài)迭代計(jì)算流程圖Fig.3 Flow chart of single assembly calculation

圖4 組件間換熱迭代計(jì)算流程圖Fig.4 Flow chart of inter-assembly heat transfer calculation

圖5 自然循環(huán)冷卻組件流量分配迭代流程圖Fig.5 Flow chart of flow distribution for natural convection assembly

3 算例測(cè)試與初步驗(yàn)證

3.1 單組件算例

以CEFR乏燃料組件典型工況為基礎(chǔ)，采用鈉冷快堆自然循環(huán)組件穩(wěn)態(tài)子通道分析程序及COBRA程序分別進(jìn)行計(jì)算，并將計(jì)算結(jié)果進(jìn)行對(duì)比。CEFR燃料組件由61根棒束構(gòu)成，棒束外徑為6 mm，下軸轉(zhuǎn)、活性區(qū)、上軸轉(zhuǎn)高度分別為250、450、100 mm；棒束采用繞絲進(jìn)行徑向定位，繞絲外徑為0.95 mm；組件盒壁內(nèi)對(duì)邊距為56.6 mm。組件的子通道劃分如圖6所示。

額定工況下，乏燃料組件功率約為5～10 kW。本文在計(jì)算時(shí)，乏燃料組件入口溫度取358 ℃，功率取6.4 kW，流量取0.06 kg/s。組件的軸向功率分布和徑向功率分布根據(jù)乏燃料組件在額定工況下的典型分布給出。在使用COBRA程序進(jìn)行子通道對(duì)比計(jì)算時(shí)，軸向壓降模型也選用詳細(xì)的CT模型，交混因子分別選取0.005和0.01。

圖6 CEFR燃料組件子通道示意圖Fig.6 Sub-channel meshing for CEFR fuel assembly

對(duì)于本單組件模型，本程序和COBRA程序計(jì)算耗時(shí)均為5 s以內(nèi)，占用內(nèi)存約100 M。二者得出的結(jié)果對(duì)比示于圖7。結(jié)果表明：本程序得出的平均出口溫度與COBRA程序相比基本無(wú)差別，平均入出口壓差相差約2.0 Pa，子通道最高溫度相差約1.1 ℃。然而，本程序得出的子通道出口的溫度分布更加平坦，子通道出口最高溫度和最低溫度相差12.3 ℃；對(duì)于COBRA-0.01和COBRA-0.005，其值分別為22.7 ℃和23.8 ℃。原因?yàn)橄啾扔贑OBRA程序，本程序所采用的MIT交混模型不僅能計(jì)算Re對(duì)交混的影響，還能計(jì)算軸向溫差、繞絲對(duì)交混的影響，因此本程序得出的交混系數(shù)更大，從而導(dǎo)致子通道之間的溫差更小。由于本程序選取的交混模型適用于低Re，因此本程序的子通道出口溫度分布結(jié)果較COBRA程序更為可信。

圖7 本程序和COBRA程序計(jì)算結(jié)果對(duì)比Fig.7 Comparison between results of this code and COBRA code

3.2 多組件算例

1) 計(jì)算輸入

多組件算例選取CEFR典型的自然循環(huán)組件為研究對(duì)象。根據(jù)CEFR的堆芯布置，計(jì)算選取了29盒組件，組件分為4類：Ⅱ型鋼反射層組件、Ⅲ型鋼反射層組件、碳化硼組件和乏燃料組件。組件盒壁厚為1.2 mm，盒間距為2 mm。組件布置如圖8所示。

Ⅱ型鋼反射層組件為強(qiáng)迫循環(huán)組件，它們是自然循環(huán)組件的邊界；其余3類組件均為自然循環(huán)組件。其中鋼反射層組件、碳化硼組件均為7棒組件，棒外徑為20 mm，無(wú)繞絲，子通道劃分如圖9所示。

圖8 組件排列示意圖Fig.8 Assembly configuration diagram

圖9 鋼反射層組件、碳化硼組件子通道劃分Fig.9 Sub-channel meshing for stainless steel and BC4 assemblies

組件的功率、強(qiáng)迫冷卻組件的流量和每盒組件內(nèi)部的軸向功率分布、徑向功率分布以CEFR堆芯的實(shí)際情況作為輸入。組件入口溫度為360 ℃。表1列出了每盒組件的功率和流量，其中流量為“-”代表自然循環(huán)組件，其流量由程序計(jì)算得出。組件編號(hào)方法為從上至下，從左至右，即左上角組件為1號(hào)組件，右下角組件為29號(hào)組件。根據(jù)CEFR堆芯漏流量和本計(jì)算模型自然循環(huán)組件功率占堆芯的比例，本模型中自然循環(huán)組件和盒間總流量約為2.11 kg/s。盒內(nèi)、盒間的局部阻力系數(shù)通過(guò)CFD得出，作為本算例的輸入。

表1 組件功率和流量Table 1 Power and flow rate of each assembly

2) 計(jì)算結(jié)果

對(duì)于CEFR多組件模型，本程序計(jì)算耗時(shí)約為1 400 s。表2列出了溫度、流量和組件間傳熱量(正為傳入，負(fù)為傳出)計(jì)算結(jié)果。程序根據(jù)自然循環(huán)組件壓降，自動(dòng)為每盒自然循環(huán)組件分配了流量。其中，Ⅲ型鋼組件(5～13)的流量為0.075～0.100 kg/s，碳化硼組件(14～23)的流量為0.042～0.053 kg/s，乏燃料組件(24～29)的流量為0.050～0.054 kg/s，盒間流量為0.52 kg/s。組件間的傳熱量為0.1～0.5 kW，相較于自然循環(huán)組件自身功率，其影響不可忽略。

對(duì)于同類組件，流量分配基本與功率相關(guān)，功率越大，組件分配的流量越大，并且其平均出口溫度和最高子通道溫度也越大。例如對(duì)于乏燃料組件，最大功率和最小功率分別為8.3 kW和6.4 kW，程序得出的組件流量分別為0.056 kg/s和0.050 kg/s，組件平均出口溫度分別為471.3 ℃和458.5 ℃，子通道最高溫度分別為479.0 ℃和463.8 ℃。

對(duì)于不同類型的組件，由于其阻力特性的差異，流量分配不完全與功率相關(guān)。如乏燃料組件功率是碳化硼組件的4～8倍，然而二者流量相差不大，這是因?yàn)榉θ剂辖M件內(nèi)部有61根棒，其棒束區(qū)阻力系數(shù)大于僅有7根棒的碳化硼組件。

根據(jù)CEFR設(shè)計(jì)報(bào)告，采用一維分析得出的3類組件的平均流量分別為0.068、0.020和0.046 kg/s[3]。本程序得出的Ⅲ型鋼組件和乏燃料的流量與設(shè)計(jì)報(bào)告基本一致，而碳化硼組件流量大于設(shè)計(jì)報(bào)告中的流量。對(duì)于碳化硼組件，造成二者差異的主要原因是本程序考慮了組件間的傳熱，且相鄰組件向碳化硼組件的傳熱量占碳化硼組件本身發(fā)熱的比例較大(50%左右)。因此，本程序得出的流量分配結(jié)果更為準(zhǔn)確。

表3列出了自然循環(huán)組件、盒間流壓差計(jì)算結(jié)果。結(jié)果表明，通過(guò)程序的流量分配計(jì)算，所有自然循環(huán)組件和盒間通道的入出口壓差均為約6 811 Pa。其中，重力壓差比重最大，占入出口總壓差的98%以上；其次是摩擦阻力，最后是局部阻力。碳化硼組件和乏燃料組件的流量相當(dāng)，但前者的摩擦阻力僅為后者的1/3；然而前者溫度較低，重力壓差更大，彌補(bǔ)了二者摩擦阻力之間的差值。

對(duì)于子通道溫度計(jì)算結(jié)果，以相鄰的19號(hào)和25號(hào)組件為例進(jìn)行分析。19號(hào)組件的9號(hào)子通道和25號(hào)組件的118、119號(hào)子通道相鄰，這3個(gè)子通道的溫度結(jié)果如圖10a所示，其中nc表示多組件自然循環(huán)模型計(jì)算結(jié)果，single表示單組件計(jì)算結(jié)果。結(jié)果表明，由于盒間傳熱的影響，19-9子通道出口溫度升高了約21 ℃，25-118和25-119出口溫度降低了約13 ℃。25號(hào)組件子通道出口溫度結(jié)果如圖10b所示。結(jié)果表明，由于25號(hào)組件邊角子通道的溫度遠(yuǎn)高于與其相鄰的19號(hào)和20號(hào)組件的邊角子通道溫度，因此25號(hào)組件部分邊角子通道溫度較單組件計(jì)算結(jié)果明顯下降；同時(shí)，由于受子通道徑向?qū)岬挠绊?，與這些邊角子通道相連的內(nèi)子通道溫度也有所下降。

表2 溫度、流量和組件間傳熱量結(jié)果Table 2 Results for temperature, flow rate and inter-assembly heat transfer capacity

表3 自然循環(huán)組件與盒間流壓差計(jì)算結(jié)果Table 3 Pressure drop of natural convection assembly and inter-assembly flow

圖10 自然循環(huán)模型和單組件模型子通道溫度對(duì)比Fig.10 Sub-channel temperature comparison between natural convection model and single assembly model

4 總結(jié)

本文基于鈉冷快堆自然循環(huán)組件的運(yùn)行工況，開(kāi)發(fā)了鈉冷快堆堆芯自然循環(huán)冷卻組件子通道分析程序。基于61棒單組件模型，通過(guò)將本程序的結(jié)果與COBRA程序進(jìn)行比較，驗(yàn)證了程序?qū)ψ匀谎h(huán)冷卻組件的適用性?；诙嗪薪M件模型，初步驗(yàn)證了本程序具備自然循環(huán)冷卻組件的流量分配和盒間換熱計(jì)算的功能，并且和一維方法相比，本程序可獲得組件內(nèi)部詳細(xì)的溫度、壓力和流速分布。本程序能為池式快堆自然循環(huán)冷卻組件提供有效的設(shè)計(jì)和分析工具。