快中子通量实验堆失流事故三维数值模拟

周 磊，张大林，刘雅鹏，王式保，王心安，王成龙，田文喜，秋穗正，苏光辉

(西安交通大学 动力工程及多相流国家重点实验室，核科学与技术学院，陕西 西安 710049)

为准确评估钠冷快堆在各类事故工况下的瞬态响应以及各类参数能否满足事故工况下的验收准则，通常采用系统程序针对钠冷快堆进行建模计算以获取事故工况下钠冷快堆的各种关键参数，如西安交通大学开发的THACS[1-2]、中国原子能科学研究院开发的FR-Sdaso[3]和FRSYS[4]、华北电力大学开发的SAC-CFR[5-6]，以及美国阿贡国家实验室开发的SAS4A/SASSYS-1[7]、俄罗斯国家科学中心物理动力研究院开发的GRIF[8]、法国原子能委员会开发的CATHARE[9]等，这类系统程序通常采用集总参数法或一维分析方法，以及有限的二维和伪三维分析方法以分析各类工况下的各种现象，为保证计算结果的可靠性，通常需要大量的实验数据用以提供各种经验关系式或物理模型。该类系统程序虽能满足钠冷快堆的设计和安全分析，但开发验证成本高，适用范围相对窄，同时在一定程度上忽略了反应堆系统局部现象，牺牲了核动力系统的经济性。伴随着高性能计算机的发展，对于核动力系统中部分系统与部件采用的全三维模拟逐渐成为可能。相较于系统程序，三维CFD方法能有效捕捉反应堆复杂的三维流动换热现象，识别潜在的安全影响因素，充分挖掘反应堆的经济性[10]。但不同于系统程序针对特定反应堆长期开展的大量验证与确认工作以及采用偏保守的物理模型，三维CFD方法从物理模型上、数值算法上具有更加广泛的适用范围与适用对象，与真实的物理现象更为接近，但针对不同流态(如层流、湍流、过渡流)、不同工质(如水、液态钠、液态铅、熔盐等)、不同工况(如强迫循环、自然循环工况)、不同流动现象(如回流、射流、绕流等)、不同流型(如单相流、气液两相流、气固两相流等)，CFD方法仍缺乏普适的湍流模型和数值算法，特别是在针对反应堆全堆级别的事故工况下的瞬态模拟验证结果仍有限。目前已有的全堆数值模拟[11-12]采用商业软件针对实验快堆的失给水工况和主泵卡死等事故工况进行全三维模拟，但建模时对堆芯组件结构进行了较大的简化，无法评估钠冷快堆重要的盒间流现象[13]在瞬态过程中对反应堆安全性能的影响，此外模拟结果缺乏足够的实验数据的验证，另一种较通用的方法是采用跨维度耦合[14-15]的方式进行模拟，针对核心关注部分采用三维CFD方法进行模拟，而针对其余部分采用系统程序进行模拟，该种方法在一定程度上结合了系统程序与CFD方法的优势，保证了精度，但对于两种不同程序的耦合又会带来收敛问题，同时对于瞬态变化剧烈的工况，难以获取较好的模拟结果。

因此，为确保CFD方法在针对钠冷快堆的低Pr流体液态金属钠[16]，以及钠冷快堆自然循环工况下[17]计算的准确性，本文利用快中子通量实验堆(FFTF)进行失流事故实验，并针对其建立相应的三维模型，其中冷池和热池按照原型尺寸建模，堆芯组件简化为多孔介质，同时保留堆芯盒间几何特性，针对无保护失流事故进行900 s的瞬态模拟，并将模拟结果与实验数据进行对比，初步验证CFD方法在针对钠冷快堆自然循环工况模拟中的可行性与可靠性。

1 FFTF系统及失流事故简介

FFTF[18]是由美国西屋公司设计建造的一座400 MW热功率，采用氧化物燃料、液态金属钠冷却的实验堆。FFTF位于华盛顿汉福德研究中心，于1980年首次临界，一直运行到1992年，并于1993年关闭，期间进行了大量实验与测试，为钠冷快堆的燃料、材料、部件积累了大量经验。1986年，为测试FFTF非能动安全特性，在该堆上进行了大量无保护瞬态实验，其中包含了13次无紧急停堆的失流事故实验(LOFWOS)，为液态金属堆的程序验证、钠冷快堆的固有安全性提供了重要支撑。2018年，在IAEA的组织下，选取LOFWOS#13作为国际基准题，多个国家和组织参与了此次基准题计算，西安交通大学作为参与单位之一，开展了系统程序和三维瞬态计算。

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

FFTF是一座环路式钠冷快堆，其典型环路及主要部件示于图1。FFTF包含2个回路、3个环路，液态钠从反应堆容器离开，分别进入3个热腿，在经过每条环路上的中间热交换器时将热量传递给二回路中无放射性的钠，最后通过二回路的空气热交换器将热量导入环境中，因此该堆并不直接用来发电。FFTF还采用了独特的安全设计以保证反应堆的安全，如较为特殊的气体膨胀组件，能减小无保护失流事故的事故后果，该组件插入到堆芯径向反射区的靠近内堆芯区，如图2所示，组件内部顶端通过不锈钢插件限制了约0.028 3 m3的氩气，造成顶端氩气无法泄漏，底端与栅板联箱连通，液钠可进入组件。因此，组件的液钠高度与入口钠的压力及氩气温度有关，额定工况下，由于主泵提供较大的扬程，使得组件内液钠的高度高于堆芯活性区，而当堆芯温度升高或失去强迫循环后，液钠的高度下降，堆芯径向中子泄漏增加，为反应堆提供负反馈，从而保障反应堆的安全性。堆芯组件分布示于图3，图中DF为燃料组件，CR、SR为控制棒组件，GEM为气体膨胀组件。图3中具有特殊标记的组件PIOTA为装有温度测点的实验组件，LOFWOS#13实验开始前保持50%额定功率、100%额定流量、约7.5 h，使得整个系统达到热平衡状态。为保证实验完整运行，反应堆安全系统在实验前进行了调整，避免控制棒过早插入，对堆内流量和通量反应堆关闭信号也进行了关闭或旁通。实验开始的初始时刻，主回路的3个环路主泵同时关闭，转入惰转，约90 s后主泵惰转结束，反应堆转入自然循环，整个实验过程中二回路泵仍保持额定运行状态。在实验进行过程中未进行其他动作与干预。实验开始前反应堆稳态参数列于表1。

图2 气体膨胀组件Fig.2 Gas expansion module

图3 堆芯组件类型及分布Fig.3 Assembly type and distribution in core

表1 反应堆稳态参数Table 1 Steady state parameter of reactor

1.1 网格模型和计算边界

FFTF堆本体如图4所示，堆内构件众多，结构复杂，尺寸跨度大，完全按照原始几何建模划分网格，将导致计算量巨大，现阶段难以完成瞬态计算。因此，针对真实反应堆结构进行了一定程度的简化。本文建模中忽略了一系列与主要流道不相关的几何腔室与结构，如图5所示，仅保留了与流动密切相关的几何结构，如冷池、热池、堆芯等关键结构，此外将各分隔板简化为无厚度的面，对于位于栅板联箱中的堆芯支持结构以及入口管脚，该结构复杂且数目众多，主要起流量分配的作用，因此在几何建模过程中忽略了该部分结构，并将该部分结构的压降及流量分配的作用归到燃料组件的入口处。堆芯的各类组件采用绕丝棒束组件，流道狭长，针对每盒组件进行了多孔介质简化，最终的网格模型如图6所示，其中堆芯组件采用六面体网格，其余部分采用四面体网格，整体网格呈混合网格，各类型网格面连续一致。

图4 FFTF堆本体Fig.4 Overview of reactor vessel

图5 FFTF简化几何模型Fig.5 Simplified geometry model of FFTF

1.2 数学模型与边界条件

1) 堆芯组件模型

由于对反应堆堆内结构进行了一定程度的简化，特别是堆芯中忽略了组件插入栅板联箱中的组件分流结构，可能直接导致堆芯各类组件的流量分配与设计值产生较大差距；另外堆芯各类几何结构复杂，形状阻力难以确定，为保证简化后流量分配与设计值保持一致，在稳态下采用如下方法进行组件进口区域的多孔介质参数调整。

图6 堆本体网格模型Fig.6 Mesh model for reactor vessel

将堆芯内各组件按照各组件的实际结构沿轴向标记为进口区、出口区、绕丝棒束区3个分区类型，在初始计算前进行如下参数设置：进口区黏性阻力系数初始指定为0，惯性阻力系数初始指定为0；出口区黏性阻力系数指定为0，惯性阻力系数指定为0，并在计算过程中保持不变；绕丝棒束区黏性阻力系数和惯性阻力系数通过经验关系式获得，并在计算过程中保持不变。

在稳态调试中，每隔n个迭代步获取堆芯每盒组件的出口流量，并根据式(1)计算堆芯内各组件出口流量W的设计值与计算值之间的相对误差εi。

(1)

式中，下标i为堆芯内组件序号，cal为计算值，des为设计值。

根据式(1)计算堆芯内各组件出口流量的设计值与计算值之间的相对误差εi，通过式(2)调整堆芯各组件进口区多孔介质的黏性阻力系数和惯性阻力系数，直至各组件的流量与设计流量的误差在允许范围，认为稳态调试完成，可开展瞬态计算。

(2)

式中：1/α为黏性阻力系数；C2为惯性阻力系数；下标n为当前迭代步，m为上一迭代步，u为用户指定修正的参数。

2) 物理模型

本文湍流模型采用Realizablek-ε模型、标准壁面函数、SIMPLE算法，梯度离散格式为Green-Gauss Node Based格式，压力离散格式为二阶格式，动量、湍动能、湍流耗散率和能量方程采用二阶迎风格式，时间离散格式采用一阶显式格式，由于密度变化不大，浮升力作用采用Boussinesq假设。

3) 边界条件和初始条件

计算时，进口采用流量进口边界和进口温度边界，出口采用压力边界。其中进口温度边界采用图7所示的瞬态温度变化，该瞬态曲线根据实际反应堆实验数据获得。进口流量变化(归一化值)如图8所示，该曲线由实验所得，该值乘以稳态流量即为瞬态流量入口边界。堆芯各类组件的稳态功率分布如图9所示，而瞬态下堆芯功率变化曲线(归一化值)如图8所示，该值乘以稳态功率即为瞬态功率变化值。需要指出的是，在瞬态计算过程假设堆芯各组件功率份额与稳态下保持一致。

图7 冷管道进口瞬态温度变化Fig.7 Transient temperature change at inlet of cold leg

图8 归一化堆芯功率和入口流量Fig.8 Normalized of core power and cold leg inlet mass flow rate

图9 稳态堆芯功率分布Fig.9 Core power distribution under steady condition

2 计算结果与分析

2.1 稳态计算结果

每盒组件经过调整多孔介质参数后流量的计算值与设计值的相对误差示于图10，每盒组件的流量的相对误差均小于±2%，且绝大多数组件流量的相对误差小于±1%，因此可认为经过调整后的多孔介质较为合理，能满足堆芯的流量分配需求，堆芯稳态计算结果是合理的。稳态下堆芯各类组件的出口温度统计示于图11，可看出，堆芯组件的最高出口温度约为699 K，最低出口温度约为602 K。

2.2 失流事故计算结果

1) 热池流动换热特性分析

反应堆堆芯出口温度即热管段温度的模拟值与实验值的对比示于图12，可看出，钠池巨大的钠装量使得在整个瞬态变化过程中温度实验值变化很小，前200 s内基本保持不变，而后700 s缓慢下降，且在整个瞬态过程仅有4 K的温度下降。CFD模拟值较好地捕捉了反应堆出口温度的变化趋势。

图10 组件流量计算值与设计值的相对误差Fig.10 Relative error between calculated value and design value of assembly mass flow rate

图11 稳态下堆芯组件出口温度Fig.11 Outlet temperature of assembly under steady condition

图12 热管段温度对比Fig.12 Comparison for hot leg temperatures

堆本体xz截面的温度分布云图示于图13。从图13可看出，初始阶段，由于堆芯组件流速较高，致使在热池中心区域形成1个高温区域，而热池外围温度相对较低，当堆芯流速降低后，中心高温区域逐渐萎缩，从热池顶部至热池底部逐渐形成明显的热分层现象，特别是在热池底部，由于流速低，从堆芯出来的冷流体可能积聚于此，直接导致此处的热分层最为明显。而对于冷池与其他腔室，整体温度变化不大。

热池xz截面不同时刻的速度矢量图和速度云图示于图14。在初始阶段，堆芯出口的流体向上流到液面处，在热池中心区域形成1个高速区，而在中心高速区外，形成了2个大的漩涡，而后从热管流出，同时伴随着流速的下降，中心高速区域逐渐减小，特别是在后期，流速较低时，热池上部的流体几乎保持静止，流体在几乎与出口相同高度位置向外流出反应堆。

图13 xz截面温度云图Fig.13 Temperature contour of xz cross section

图14 热池xz截面的速度矢量图和速度云图Fig.14 Vector and contour of velocity for xz cross section of hot pool

2) 堆芯流动换热特性分析

具有温度测点的组件的出口温度模拟值与实验值的对比示于图15，该组件在堆芯中的位置如图3所示。对于Row 6组件的温度预测，实验值除出现2个峰值温度外，在100～150 s左右，温度出现1个较为平缓的趋势，而后在150 s后又开始下降，模拟值没有体现该温度变化趋势，但考虑到Row 6组件在堆芯中的位置较特殊，周围为气体膨胀组件和非燃料组件，是堆内热工和物理条件最为复杂的组件，从热工角度来说，外侧非燃料组件温度较低，而内侧为燃料组件，温度较高，此处温度梯度较大，盒间的导热或对流换热作用较强烈；从物理角度来说，瞬态下，由于气体膨胀组件失去强迫循环，气体膨胀组件液位下降，堆芯中子泄漏增加，可能导致该组件的功率份额分布与稳态分布差距较大，属于物理热工强耦合作用，单独考虑热工影响，可能不能准确模拟该过程。

图15 Row 6 PIOTA组件出口温度Fig.15 Outlet temperature of Row 6 PIOTA

堆芯活性出口截面的温度分布示于图16。可看出，堆芯最热组件的位置在堆芯中不断发生变化，刚开始位于堆芯中心区域，而后靠近燃料组件的外围区域，最终呈现出中心低、中间环型区域高、外侧低的马鞍形温度分布，而从堆芯y=0截面的温度分布曲线可看出，刚开始时堆芯功率较高，y=0截面温度分布差距较大，而伴随瞬态进程，由于盒间的存在，温度分布逐渐被展平。

不同时刻下盒间出口的流动特性示于图17。在初始阶段，从外侧非燃料组件区域向下流入盒间，而从中心较高温区的燃料区域流出盒间，而在后期，由于堆芯温度分布产生了较大变化，使得盒间从外部和中心区域流入堆芯，而从外侧燃料区域流出堆芯。通过对比盒间的矢量图和堆芯xy截面的温度分布云图，可看出矢量图流动方向基本与温度保持一致，都是从低温区流入盒间，从高温区流出盒间，这是由于低温时密度较大，而高温时密度较小，在浮升力的作用下，呈现这种流动趋势。

图17 盒间出口截面流动特性Fig.17 Flow characteristic of outlet for inter wrapper

3 结论

本文采用CFD方法针对FFTF的LOSWOS#13无保护失流事故进行了900 s的瞬态模拟，并将模拟值与实验值进行了对比，得出如下结论：

1) CFD方法能较好地捕捉钠冷快堆全堆的热池内热分层等复杂的三维流动换热特性，特别是对热池巨大的热容及其在瞬态下的热惰性响应；

2) 堆芯最热组件在瞬态过程中的位置不断发生变化，堆芯盒间流的流动方向受盒内温度影响较大，且盒间在堆芯流速低、功率小的情形下，能有效展平堆芯温度。