鈉冷快堆自然循環組件子通道程序初步研發

林 超，楊紅義，周志偉

(中國原子能科學研究院 反應堆工程技術研究部，北京 102413)

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

中國示范快堆和實驗快堆均屬于池式鈉冷快堆，其堆芯組件按冷卻方式可分為兩種類型：1) 強迫循環冷卻組件，包括燃料組件、轉換區組件、控制棒組件等，這些組件的發熱量較大，通過一次鈉泵提供足夠的冷卻劑流量對其進行冷卻；2) 自然循環冷卻組件，包括乏燃料組件、鋼組件反射層和碳化硼組件，這些組件的發熱量較小，所需冷卻劑的流量也較小，無需泵提供強迫循環流量，而是依靠其自身發熱產生溫差驅動流體進入組件內部對其進行冷卻。

對于自然循環冷卻組件，無法對組件流量進行精確分配，而是依靠組件自身的發熱進行流量的自動分配。其流量來源于堆芯漏流，但漏流不僅可進入自然循環組件內部，也能進入組件盒間形成盒間流。因此，自然循環組件熱工分析十分復雜[2-3]。

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

本文基于鈉冷快堆自然循環組件的特點以及子通道方法，開發鈉冷快堆自然循環組件穩態子通道分析程序，并對其進行驗證。

1 程序模型

1.1 單相流穩態子通道守恒方程

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

對于鈉冷快堆，冷卻劑溫度距其沸點有較大的裕度，因此無需考慮相變。單相穩態質量守恒方程為：

(1)

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

穩態能量守恒方程為：

(2)

穩態軸向動量守恒方程為：

(3)

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

穩態橫向動量守恒方程為：

(4)

式中：l為間隙控制體等效長度，m；KG為橫向流動形阻系數；ρ為子通道密度，kg/m3；下標i、j代表與該間隙相鄰的子通道編號。橫流方程針對的是所有子通道之間的間隙。

在進行守恒方程求解時，采用有限差分法對空間項進行離散。

1.2 軸向壓降模型

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

fi=Cfi/Rem

(5)

式中：m=1(湍流)或0.18(層流)；Cf為各類子通道摩擦系數因子；下標i代表不同子通道類型。Cf主要由光棒束的摩擦阻力與繞絲的橫掠阻力疊加而成，其值與當地Re、子通道類型、繞絲螺距、棒束直徑、棒束中心距等參數有關。對于不同類型的子通道，其Cf的計算方法不同。對于過渡區，摩擦系數計算方法如下：

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

(6)

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

(7)

式中：下標Tran代表過渡區，Lam代表層流區，Tur代表湍流區，L代表層流區到過渡區的轉換點，T代表過渡區到湍流區的轉換點；λ為擬合常數，一般取13。

1.3 橫向壓降模型

本文采用gunther關系式[9]研究棒束區的軸向壓降關系。

對層流Re<200，有：

(8)

對湍流Re≥200，有：

(9)

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

1.4 換熱模型

換熱模型采用Mikityuk關系式[10]。Mikityuk在2009年對國際上的液態金屬傳熱關系式進行了分析對比，并對實驗結果進行了重新擬合。相對于其他換熱關系式，此擬合關系式精度最好，該關系式如下：

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

(10)

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

1.5 交混模型

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

根據MIT交混模型，在低Re自然對流情況下，交混系數表達式為：

Γm=ε1+εM+κυ

(11)

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

1.6 組件間傳熱模型

本文在處理組件間傳熱時，將組件間的鈉劃分為控制體，以計算該控制體與其周圍組件邊子通道或角子通道之間的換熱。針對所有組件的每個參與組件間傳熱的邊、角子通道，分別求出它們通過組件間傳熱導致的換熱量，然后更新相應子通道的能量方程。組件間控制體劃分如圖1所示[13]。由于快堆盒間流速很低，因此本文在處理盒間傳熱時將盒間流視為滯止狀態，簡化為固體控制體，該簡化方法是保守的。

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

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

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

(12)

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

1.7 流量分配模型

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

在計算入口壓差時，需分別計算重力壓差、摩擦阻力及局部阻力。對于重力壓差，根據盒內、盒間沿軸向溫度分段計算。盒內的摩擦阻力系數采用1.2節的軸向壓降模型。對于盒間流，其Re較低，處于層流狀態，且其流道結構較為簡單，因此盒間摩擦阻力系數ξ通過下式計算[15]：

ξ=64/Re

(13)

由于盒內、盒間流道較為復雜，其局部阻力特性無法通過查詢阻力手冊的方法得出，因此需通過CFD計算來獲得，作為本程序的輸入。

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

(14)

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

2 計算流程

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

圖2 總體計算流程圖Fig.2 Flow chart of overall calculation

2.1 單組件穩態迭代計算

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

2.2 組件間換熱迭代計算

在進行組件間換熱迭代計算時，根據上一次單組件穩態迭代計算得出的組件邊角子通道的溫度和相應熱阻，得出組件間的換熱量，并通過多次迭代，使相鄰組件間的溫差和傳熱量相匹配。組件間換熱迭代計算流程圖如圖4所示。

2.3 自然循環冷卻組件流量分配迭代計算

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

圖3 單組件穩態迭代計算流程圖Fig.3 Flow chart of single assembly calculation

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

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

3 算例測試與初步驗證

3.1 單組件算例

以CEFR乏燃料組件典型工況為基礎，采用鈉冷快堆自然循環組件穩態子通道分析程序及COBRA程序分別進行計算，并將計算結果進行對比。CEFR燃料組件由61根棒束構成，棒束外徑為6 mm，下軸轉、活性區、上軸轉高度分別為250、450、100 mm；棒束采用繞絲進行徑向定位，繞絲外徑為0.95 mm；組件盒壁內對邊距為56.6 mm。組件的子通道劃分如圖6所示。

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

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

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

圖7 本程序和COBRA程序計算結果對比Fig.7 Comparison between results of this code and COBRA code

3.2 多組件算例

1) 計算輸入

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

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

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

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

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

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

2) 計算結果

對于CEFR多組件模型，本程序計算耗時約為1 400 s。表2列出了溫度、流量和組件間傳熱量(正為傳入，負為傳出)計算結果。程序根據自然循環組件壓降，自動為每盒自然循環組件分配了流量。其中，Ⅲ型鋼組件(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，相較于自然循環組件自身功率，其影響不可忽略。

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

對于不同類型的組件，由于其阻力特性的差異，流量分配不完全與功率相關。如乏燃料組件功率是碳化硼組件的4～8倍，然而二者流量相差不大，這是因為乏燃料組件內部有61根棒，其棒束區阻力系數大于僅有7根棒的碳化硼組件。

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

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

對于子通道溫度計算結果，以相鄰的19號和25號組件為例進行分析。19號組件的9號子通道和25號組件的118、119號子通道相鄰，這3個子通道的溫度結果如圖10a所示，其中nc表示多組件自然循環模型計算結果，single表示單組件計算結果。結果表明，由于盒間傳熱的影響，19-9子通道出口溫度升高了約21 ℃，25-118和25-119出口溫度降低了約13 ℃。25號組件子通道出口溫度結果如圖10b所示。結果表明，由于25號組件邊角子通道的溫度遠高于與其相鄰的19號和20號組件的邊角子通道溫度，因此25號組件部分邊角子通道溫度較單組件計算結果明顯下降；同時，由于受子通道徑向導熱的影響，與這些邊角子通道相連的內子通道溫度也有所下降。

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

表3 自然循環組件與盒間流壓差計算結果Table 3 Pressure drop of natural convection assembly and inter-assembly flow

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

4 總結

本文基于鈉冷快堆自然循環組件的運行工況，開發了鈉冷快堆堆芯自然循環冷卻組件子通道分析程序?；?1棒單組件模型，通過將本程序的結果與COBRA程序進行比較，驗證了程序對自然循環冷卻組件的適用性?；诙嗪薪M件模型，初步驗證了本程序具備自然循環冷卻組件的流量分配和盒間換熱計算的功能，并且和一維方法相比，本程序可獲得組件內部詳細的溫度、壓力和流速分布。本程序能為池式快堆自然循環冷卻組件提供有效的設計和分析工具。