基于SAC-3D系统软件的FFTF未能紧急停堆失流实验的数值模拟

陆道纲，吕思宇,*，隋丹婷

(1.华北电力大学 核科学与工程学院，北京 102206；2.非能动核能安全技术北京市重点实验室，北京 102206)

美国快通量试验堆(FFTF)是由美国西屋电气公司设计建造的一座以液态金属钠为冷却剂的快中子反应堆。该快堆在1986年进行了一系列无保护瞬态工况的实验来验证FFTF堆型的非能动安全性。其中有13组实验为未能紧急停堆的失流(LOFWOS)事故工况，目的是验证液态金属冷却快堆的安全裕度。2017年，国际原子能机构(IAEA)启动了FFTF LOFWOS基准题验证的国际合作研究项目(CRP)。该项目的主要目的是提高成员国在快堆模拟和设计领域的分析能力，促进成员国通过国际合作研发实现快堆技术发展的技术进步[1]。共有来自13个国家的25个组织参与了该项目的验证工作。华北电力大学作为项目参与单位参与了FFTF CRP中的中子与热工水力系统分析验证工作。

针对钠冷快堆，华北电力大学在自主开发的SAC-CFR[2-4]系统分析软件的基础上，开发了具有三维时空动力学中子物理计算模块的系统分析软件SAC-3D[5]。本文应用SAC-3D对FFTF堆芯进行详细建模计算，以得到反应堆堆芯稳态功率分布、瞬态事故中反应性变化、组件峰值温度、开放式测试组件PIOTA(proximity instrumented open test assembly)出口处冷却剂温度瞬态变化曲线等关键参数。

1 FFTF未能紧急停堆的失流实验

1.1 FFTF概述

FFTF是一座热功率为400 MW回路式快中子反应堆，堆芯使用MOX燃料，并以液态金属钠作冷却剂。FFTF的冷却剂系统中有3套环路，每套环路均由1条一回路与1条二回路构成，在一回路中布置有1台中间热交换器(IHX)和1台主泵，图1为FFTF冷却剂系统示意图[1]。FFTF的反应堆容器为圆柱形，底部是球形封头。冷钠通过3条一回路的冷管段经冷却剂入口流入反应堆容器底部入口腔室，向上流经堆芯支撑结构，通过堆芯组件、径向屏蔽和旁路流道被堆芯加热，随后热钠从堆芯容器出口流出，进入到3条一回路热管段，流经主泵进入到IHX，通过IHX的壳侧将热量传递给管侧的二回路冷却剂，冷却后的液态钠经一回路冷段管道再次流入反应堆容器。由于FFTF不产生电力，二回路中被IHX加热的冷却剂流出IHX后直接进入废热交换器(DHX)，堆芯产生的热量最终将通过DHX传递到最终热阱大气中。

图1 FFTF冷却剂系统示意图Fig.1 Coolant system overview of FFTF

1.2 未能紧急停堆的失流实验

本文所建模的对象为由IAEA FFTF合作研究项目所提供的第13组未能紧急停堆的失流实验 (LOFWOS Test #13)，该实验的反应堆主要初始状态参数列于表1，堆芯功率为正常运行功率的50%，冷却剂流量维持在正常运行功率等级。实验开始后，3个一回路主泵同时关闭，所有二回路冷却机泵保持正常工作状态，反应堆保护系统不动作，控制棒、安全棒维持实验开始前的状态，模拟由一回路主泵失效引起的LOFWOS事故工况。该实验是为了验证FFTF堆型的非能动安全裕度，本文针对该工况进行建模并开展稳态、瞬态计算。

表1 FFTF LOFWOS Test #13初始状态参数Table 1 Initial condition of FFTF LOFWOS Test #13

2 计算方法

2.1 堆芯模型

基于SAC-3D对FFTF堆芯的建模分为两部分：1) 中子物理模型；2) 热工水力模型。SAC-3D的中子物理计算模块基于高阶六边形节块展开法[6]开发，该模块需要堆芯材料的宏观截面作输入项进行初始化。本文中堆芯材料均匀化宏观截面库的制作使用的是欧洲反应堆分析优化计算系统ERANOS 2.0中的栅元计算模块ECCO[7]和ENDF/B-Ⅶ[8]核数据库。FFTF基准题技术规格书中提供了全堆芯所有组件不同高度区域材料的详细组分及原子密度数据[9]，为减少运算负担，本文根据材料所属的组件类型、区域、关键元素原子密度等变量将所有材料整合为14种代表性材料进行栅元计算。对于燃料组件的燃料区域及安全棒/控制棒组件的吸收体区域的材料，本文进行真实的几何描述；对于其他栅元，本文采用均匀化材料假设。在栅元计算中，本文对可裂变材料采用ECCO中的6步计算步骤；对不可裂变的材料，其中安全棒/控制棒吸收体材料采用4步计算步骤，其余材料采用2步计算步骤[10]。由于次临界材料没有内中子源，本文在计算中将采用与其临近的材料的中子能谱作为其外中子源进行栅元计算。14种材料的栅元计算在500、750和1 000 K的温度下得到了33群组件均匀化宏观截面库，瞬态计算中其余温度下的材料截面数据将由线性插值得到。

堆芯中子物理计算建模的范围为径向包含所有组件，轴向从堆芯活性区域底部到活性区域顶部。在FFTF中子物理建模中，径向上每个组件均设置为1个计算节点；轴向如图2所示，根据模型高度及节点数量敏感性分析结果划分为38个计算节点，计算节块高度为2.5 cm，模型高度为95 cm，与堆芯活性区域高度一致。

图2 FFTF堆芯中子物理模型轴向节点划分Fig.2 Schematic of neutronics calculation model nodes division of FFTF core in axial direction

堆芯热工水力模型采用单通道模型[2-3,11]，不考虑组件内部的交叉流动，同时不考虑燃料元件的轴向导热，因此，通道内冷却剂稳态工况的动量方程可表示为：

pi,in-pi,ex=f(Li,De,i,Ai,Wi,μi,ρi)

(1)

其中：pi,in、pi,ex分别为第i个通道的进、出口压力；Li、De,i、Ai分别为通道的长度、当量直径和流通截面积；Wi、μi、ρi分别为质量流量、黏性系数和密度。式(1)右边的压降包括提升压降、加速压降、摩擦压降和局部压降。

燃料元件内能量以导热方式传递，燃料元件-包壳间的热传递采用间隙导热模型，包壳与冷却剂间的热传递方式为对流换热。柱坐标下燃料元件导热微分方程的形式为：

(2)

其中：qv为体积释热率，W/m3；κ为热导率，W/(m2·℃)；r为径向坐标。

根据组件类型、在堆芯所处位置、内部冷却剂的质量流量，将FFTF堆芯分为9个通道(表2)。其中1～6通道为燃料组件，7～9通道分别代表控制棒/安全组件通道和两种屏蔽组件通道。

表2 热工水力模型堆芯通道初始质量流量对比Table 2 Initial mass flow comparison between calculation and measurement results

2.2 GEM组件瞬态计算处理方法

GEM组件是FFTF堆型中特有的非能动安全组件，在LOFWOS实验中，9个GEM组件被布置在燃料组件与屏蔽组件之间。GEM组件顶部密封底部开放，内部有约0.028 3 m3的空腔，空腔内部充有氩气。在正常工况中，由GEM组件外部液态钠与一回路主泵提供的压头将GEM内部冷却剂液面维持在超出堆芯活性区域顶部30.48～40.64 cm的高度。当主泵停转发生失流事故时，冷却剂压头急剧下降，GEM组件内氩气膨胀，钠液面下降至接近堆芯活性区底部，堆芯中子泄漏率急剧上升，引入负反应性。SAC-3D三维时空动力学计算模块的粗网格节块法在每个计算节点使用的是单一材料截面，而在失流事故瞬态过程中，GEM组件内的液面高度会随着堆芯质量流量下降而连续变化，这导致在瞬态计算中会出现在液面所在的节块中同时存在两种材料的情况。因此，在对FFTF进行失流事故瞬态计算的过程中，需对GEM组件内材料的截面进行在线均匀化处理。截面均匀化方法如图3所示，本文假设组件径向边界条件为全反射并保持轴向边界条件与堆芯中子计算边界条件相同。在轴向进行一维SN计算得到轴向的中子通量分布，并按照式(3)进行截面均匀化得到该节块新的平均截面数据[12-13]。

图3 GEM组件截面均匀化方法Fig.3 Cross-section homogenization process of GEM assembly

(3)

其中：Σi包括裂变截面、散射截面、吸收截面、移出截面及总截面等；Σi1、Σi2分别为节块内自身材料与相邻节块内原材料的各种截面；h为节块高度；h′为节块底部到两种材料交界处的高度；φ为中子通量密度。

2.3 IHX模型

IHX模型[2]种将所有传热管简化为1根传热管进行传热计算，假设进口和出口区域内流体完全混合，换热管区域内的换热为充分发展段的对流换热。径向节点分为4类：二次侧流体、换热管、一次侧流体、外壳，其中换热管和外壳温度节点定义在控制体的中心，而流体内温度节点定义在控制体界面处。

一次侧流体i与i+1间的控制体能量方程为：

HptApt(Tp,i,i+1-Tt,i)-HpshApsh(Tp,i,i+1-Tsh,i)

(4)

二次侧流体i与i+1间的控制体能量方程为：

HstAst(Tt,i-Ts,i,i+1)

(5)

换热管管壁第i个控制体的能量方程为：

HpstApst(Tt,i-Ts,i,i+1)

(6)

外壳第i个控制体的能量方程为：

(7)

其中：下标p表示一次侧流体，s表示二次侧流体，t表示换热管，sh表示外壳；e为焓；W为质量流量；T为温度；M为换热管质量；c为比热容；V为控制体体积；H为流体与结构件的换热系数；A为流体与结构间的换热面积；Ti,i+1为流体控制体内的平均温度，定义式为：

(8)

IHX二次侧入口处为本文建模的系统边界，根据美国阿贡实验室提供的IHX二回路冷却剂入口处液态钠参数的瞬态数据设定边界条件。在进行瞬态求解时，对进出口区域、换热区域内的节点方程采用一阶全隐式积分，而换热量则采用显式方法确定，沿着流动方向进行求解。

3 计算结果与讨论

图4为SAC-3D对FFTF LOFWOS实验基准题进行仿真计算后得到的初始堆芯功率分布结果。通过对比本文所使用的SAC-3D、ANL使用的SAS4A/SASSYS-1和日本原子能机构JAEA使用的Super-COPD得到的功率分布结果[14-15]，得到了堆芯功率分布计算结果的变异系数分布(图5)。计算偏差最大的组件为DF4.2型燃料组件，变异系数为3.27%。在该组件周边的其他组件也存在着相对较大的变异系数。通过对比所使用的计算方法与模型，本文认为该区域的计算偏差有一部分是来自于中子物理计算方法、核数据库与热工水力简化方法的区别。由于缺少真实测量数据，无法客观评估3套软件对堆芯功率分布预测结果的准确性。本文将在IAEA公布组件功率的真实测量数据后对该部分的偏差进行进一步分析。

图4 FFTF LOFWOS初始堆芯功率分布Fig.4 Initial power distribution of FFTF LOFWOS

图5 堆芯功率计算变异系数分布Fig.5 Distribution of calculation variation coefficient of core power

本文对实验前900 s进行了计算。实验开始后，由于主泵的跳闸，反应堆一回路流量急剧下降，引入较大的负反应性，导致反应堆功率迅速下降(图6)。从0～25 MW局部功率区间的实验与计算结果对比(图7)可看出，SAC-3D的预测结果较实际功率高，通过分析反应性反馈计算结果发现，功率预测的偏差是因为净反应性的计算结果较实验测量值高。净反应性计算的偏差可能来自以下几方面：1) GEM组件内液面随流量变化采用的是由ANL提供的流量-液位曲线，与实验中GEM内液位实际变化情况会有出入；2) 在反应性反馈模型中本文未考虑由于组件轴向温度梯度造成的燃料弯曲现象引入的反应性反馈，该反应性反馈在中国实验快堆(CEFR)的物理设计计算中引入的是负反应性反馈[16]。图8中，在实验开始的前9 s内随着主泵停转，流量迅速下降，组件温度迅速上升，引入了较大的正的Doppler反应性，同时在GEM组件中钠液位急剧下降，堆芯中子泄漏率陡升引入了巨大的负反应性，导致净反应性迅速下降；随后由于总功率的迅速下降，组件温度出现短暂回落，引入了正的净反应性反馈，出现第1个净反应性低点；随着冷却剂流量继续下降，组件温度再次回升，GEM组件由于液面已经接近堆芯活性区底部，对净反应性的负反馈贡献减弱，此时组件的径向膨胀及控制棒驱动机构的轴向膨胀引入的负反应性占据上风，净反应性再次下降；冷却剂建立起自然循环，流量逐渐稳定，而裂变功率持续下降导致组件温度再次下降，引入正的反应性反馈，净反应性出现第2个低点，此后组件温度进入长期冷却阶段，净反应性稳定在较小的负值处。图9为GEM组件内液面高度与引入的反应性关系曲线，图10为稳态计算中两个PIOTA组件内部的温度分布。

图6 FFTF LOFWOS瞬态功率曲线Fig.6 Transient power profile of FFTF LOFWOS

图7 FFTF LOFWOS瞬态功率曲线(低功率区间)Fig.7 Transient power profile of FFTF LOFWOS (low power range)

图8 反应性瞬态曲线Fig.8 Transient reactivity profile

图9 GEM反应性-液面高度关系曲线Fig.9 GEM reactivity vs. liquid level

图10 PIOTA组件轴向温度分布Fig.10 Temperature distribution of PIOTA assembly in axial direction

4 结论与展望

1) 实验测量结果与仿真计算结果均表明FFTF堆型具有良好的非能动安全性，在50%功率未能紧急停堆的失流事故工况下可依靠自身的负反应性反馈将反应堆功率迅速降下来，并建立起冷却剂自然循环长期冷却反应堆堆芯。

2) GEM组件在FFTF堆型非能动安全设计中起着重要的作用。准确模拟该组件在事故瞬态中的行为对反应堆安全设计及分析十分重要，通过初步建模模拟，发现在一回路流量急剧下降的过程中，GEM组件中的液钠液面在迅速下降过程中会出现短暂回升的现象，该现象在反应堆事故瞬态中将引入正反应性反馈。关于GEM组件在事故瞬态中对堆芯反应性的影响本文将在之后工作中继续深入研究。

3) 通过SAC-3D对FFTF LOFWOS的建模计算结果与实验测量数据对比、软件与软件建模验证，证明了SAC-3D对回路式液态金属冷却快堆进行仿真计算与安全分析的能力。但该系统分析程序的计算方法仍需进一步完善，如堆芯采用的单通道模型忽略了组件与组件之间的质量、能量传递，在一定程度上引入了误差；三维时空动力学应用节块展开法求解中子扩散方程，该方法在强吸收体及散射各向异性较强的组件区域有着较大误差。在之后工作中，本文将对SAC-3D的计算模块进一步优化，修正近似假设带来的误差。