钠冷快堆自然循环组件子通道程序初步研发

林 超，杨红义，周志伟

(中国原子能科学研究院 反应堆工程技术研究部，北京 102413)

钠冷快堆作为核能系统国际论坛(GIF)提出的6种第4代核电站堆型之一，具有核能的可持续发展、更高的安全性和可靠性、更高的经济性、防止核扩散等优点。同压水堆相比，钠冷快堆能将铀资源的利用率由1%提高到60%以上[1]。基于上述优点，中国正在大力发展钠冷快堆技术，并计划于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个流程。总体计算流程图如图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 总结

本文基于钠冷快堆自然循环组件的运行工况，开发了钠冷快堆堆芯自然循环冷却组件子通道分析程序。基于61棒单组件模型，通过将本程序的结果与COBRA程序进行比较，验证了程序对自然循环冷却组件的适用性。基于多盒组件模型，初步验证了本程序具备自然循环冷却组件的流量分配和盒间换热计算的功能，并且和一维方法相比，本程序可获得组件内部详细的温度、压力和流速分布。本程序能为池式快堆自然循环冷却组件提供有效的设计和分析工具。