基于CFD的棱柱型模块式高温气冷堆典型事故研究

黄 政，周 喆，王贺南，刘国明，陈巧艳

(中国核电工程有限公司，北京 100840)

作为第四代先进反应堆技术之一，高温气冷堆(HTGR)设计紧凑简化，且具有良好的固有安全性，是国际公认的安全性好、用途广泛的先进堆型[1]。基于该技术开发的微型核电装置能量密度高、体积重量小，从而能为孤立偏远地区提供机动灵活、清洁可靠的电源供应，对于经济建设具有重要的应用前景。

根据燃料组件形状的差异，高温气冷堆的堆芯设计分为球床型和棱柱型两大类。其中，球床型高温气冷堆的研究以清华大学为代表，如10 MWt的高温气冷实验堆HTR-10[2]及石岛湾200 MWe高温气冷堆示范工程(HTR-PM)[3-4]。与此同时，中核集团也在积极开展棱柱型高温气冷堆的研发设计。其堆芯设计采用卧式放置方式，功率密度约为3 MW/m3。位于中心的活性区由若干个六棱柱形的燃料组件排列成六边形形状。燃料芯块由TRISO包覆颗粒和圆柱形SiC基体组成。慢化剂为耐高温的棱柱石墨块，每个组件的棱柱石墨块内包含独立的6个燃料芯块通道和12个氦气冷却剂通道。活性区在径向方向外侧包裹有反射层结构，反应性控制转鼓机构位于反射层内；反射层外则为反应堆压力容器。压力容器外侧的反应堆腔室所形成的流道构成非能动余热排出系统(RCCS)，通过对流和辐射换热作用将反应堆的热量传递给进入RCCS流道的自然循环空气以及RCCS腔室墙体，并最终导入外部环境。该系统用于维持正常运行条件下反应堆腔室的温度，以及在事故条件下将堆芯热量导入最终热阱(外部大气环境)，从而保证堆芯结构的完整性。

高温气冷堆的事故发展进程涉及多个物理场的相互耦合作用，包括流体(氦气冷却剂)流动、流体以及固体(燃料、石墨慢化剂、控制棒和压力容器壁)之间的传热、中子动力学、放射性裂变产物迁移、石墨氧化化学腐蚀等。其中，流体域和固体域的温度是判断燃料和设备完整性的主要验收准则。为满足安全审评、建造许可和运行监管的要求，需针对各类包络假设事故工况开展全面系统的分析评价，从而验证堆型设计的固有安全特性，并为相关安全系统(如余热排出系统)和设备(如压力容器和控制转鼓)的设计提供必要的输入条件。但现阶段针对棱柱型设计的HTGR的计算分析工具和针对该堆型的事故分析研究进展相对滞后。首先，当前用于HTGR分析的计算程序相对较少，主要分为两大类：集总参数程序(如TINTE[5]、MGT-3D[6]、THERMIX[7]、GAMMA[8]、MELCOR[9]、RELAP-3D[10])和基于计算流体力学(CFD)的程序(如Pronghorn[11])。但上述程序大多以二维集总参数程序为主，无法给出局部精细的三维温度场分布；且多数程序只适用于球形燃料，无法应用于棱柱型堆芯设计(如TINTE、THERMIX、RELAP-3D和Pronghorn)。另一方面，一些学者对HTGR的各类典型事故工况开展了数值分析，但主要针对的是球床式堆型设计。如针对清华大学高温气冷实验堆HTR-10，王岩等[12]采用RETRAN-02热工程序计算了蒸汽发生器传热管发生断裂事故下的进水，Gou等[13]和Chen等[14]分别采用THERMIX程序对由于控制棒抽出引入正反应性且未能停堆(ATWS)的试验工况开展了模拟和模型验证；针对石岛湾高温气冷堆示范工程HTR-PM，徐伟等[15]采用TINTE-TIIXUW程序对热气导管双端断裂导致的进气事故的氧化后果进行了评估，Zheng等[16]采用TINTE和SPECTRA两种程序对丧失强迫冷却事故工况(包括失压DLOFC和未失压PLOFC两种情形)进行了模拟，论证了设计的安全性。但目前针对棱柱型高温气冷堆事故的分析仍相对较少，这主要是由于棱柱型设计起步相对较晚，适用的分析程序和模型也正在开发中。相关的工作包括Beeny等[17]使用MELCOR程序针对MHTGR概念设计方案，分析了DLOFC和PLOFC两种事故工况下的堆芯瞬态响应进程。

因此，为评估正在研发的棱柱型HTGR事故条件下的安全特性、给关键设备系统设计提供必要参数，本文基于通用CFD程序COMSOL Multiphysics建立详细三维计算模型，采用C++编写点堆求解模块并实现中子物理-传热-流动的耦合计算。以正常运行工况下的稳定状态结果作为初始条件，分别对3个典型事故瞬态工况开展模拟，包括未失压丧失强迫流动冷却(PLOFC)事故、未失压丧失强迫流动冷却且未能停堆(PLOFC+ATWS)事故，以及反应性引入且未能停堆(RIA+ATWS)事故。在此基础上，针对压力容器壁与非能动余热排出系统的辐射发射率开展敏感性分析，研究该安全系统的换热能力对事故缓解效果的影响。

1 数值模型

本研究所采用的COMSOL为通用CFD程序，该程序并没有针对棱柱型高温气冷堆的专用模型，因此需进行必要的模型开发，包括几何实体建模、氦气和燃料等材料物性定义、中子动力学模型开发以及实现中子物理-流动-传热耦合计算等。

1.1 几何实体与网格

本文所建立的反应堆堆芯和燃料组件几何实体如图1所示。

图1 棱柱形HTGR堆芯和燃料组件示意图Fig.1 Schematic diagram of prismatic HTGR core and fuel assembly

为节省计算量和资源，利用轴对称特性，本研究只针对1/6堆芯建立三维几何模型，计算域包括燃料组件(燃料通道、冷却剂通道、石墨慢化剂)以及侧反射层。由于控制转鼓详细结构和物性参数尚未定型，暂不模拟该结构。对于组件之间以及侧反射层内的细小装配缝隙，本研究也暂不考虑；缝隙内的氦气流量将以旁流的方式从堆芯流量中扣除。由于压力容器壁面厚度相对于堆芯径向尺寸较薄，且导热性较好，因此在建模时选择采用零厚度壁面来近似，而不必画出其具体实体；其导热热阻通过COMSOL的thin layer特性来添加。

剖分网格时，首先在径向方向上采用非结构化三角形网格进行划分。为更准确模拟氦气冷却剂的近壁面流动特性，在冷却剂通道壁面还划分了6层边界层。在轴向方向上采用结构化网格拉伸，共划分了10层。网格质量总体良好，能满足计算要求。初步计算时划分的网格单元总数约为50万。

1.2 流动传热模型

对于燃料芯块，采用物性均匀化的方式来处理，而不具体模拟其内部的每个TRISO球形颗粒。燃料芯块与燃料通道之间缝隙内的传热(包括热辐射、氦气导热)采用等效热阻串/并联的方式来计算[17]。其中，燃料芯块与通道内表面之间辐射换热的等效热导hrad表达式[18]为：

(1)

根据正常运行工况下氦气流量的名义值估算，氦气流动为湍流流型，计算时采用标准k-ε湍流模型。借助COMSOL提供的多物理场耦合接口，流体和固体之间的流动-传热以耦合的方式来计算。

对于压力容器壁面与非能动余热排出系统之间的换热，采用对流边界条件来等效模拟，即给定1个恒定外部环境温度(30 ℃)以及1个保守的对流换热系数(5 W/(m2·K))。

1.3 中子点堆模型

对于未能正常停堆的ATWS类事故以及反应性引入相关的事故，需计算反应堆中子通量(功率)的瞬态变化过程。本研究采用中子点堆模型来计算，即假设反应堆内的中子通量(功率)分布形状在整个事故进程中始终不变。点堆模型无法描述中子的空间分布特征，但计算速度最快，也是主流事故分析系统程序(如MELCOR、RELAP、MAAP等)普遍采用的模型。

中子点堆模型是包含1组瞬发中子(式(2))和6组缓发中子(式(3))的核子密度动态平衡方程：

(2)

(3)

其中，n(t)和Ci(t)分别为瞬发中子和第i组缓发中子的密度；ρ(t)为总的瞬时反应性，由以下几部分构成：

ρ(t)=ρext+ρfuel+ρgraphite

(4)

式中：ρext为外部引入的反应性(如控制棒插入/抽出)；ρfuel、ρgraphite分别为由多普勒效应所引起的燃料和石墨慢化剂的温度反馈反应性。

考虑到该方程组具有较大的刚性特性，本研究选择采用具有预测-校正的Gear数值算法来求解，该算法能保证计算的精度，同时又避免时间步长过小的问题[19]。采用C++语言编写了求解代码，并编译成动态链接库(DLL)，以外部函数的形式供COMSOL程序在瞬态计算过程中调用。事故瞬态计算过程中，点堆中子动力学模块首先从COMSOL程序获取上一时刻的燃料芯块以及石墨慢化剂的平均温度，再分别通过燃料以及慢化剂的温度反馈系数计算得到反馈的反应性，然后求解点堆方程组得到新时刻的中子密度，最后转换成新的堆芯功率水平返回给COMSOL程序，以内热源的方式添加到燃料芯块计算域中，从而实现耦合计算。

2 正常运行稳态工况分析

事故工况的瞬态计算需以正常运行条件下的稳态结果(温度分布和冷却剂流场)作为初始值。模拟正常运行工况时，反应堆裂变功率是根据中子三维输运计算得到的分布形状按照归一化后的比例来进行分配的。功率分布形状大致为：径向方向由堆芯中心向四周功率逐渐降低；轴向方向则近似呈截断余弦分布，功率峰值大致位于中段位置。氦气冷却剂通道入口为给定流量以及温度边界条件，出口则为定压边界条件。压力容器外壁面设置为对流换热边界条件(环境温度为30 ℃，等效对流换热系数为5.0 W/(m2·K))，用于模拟非能动余热排出系统带热。圆柱θ方向的两个对称截面设为对称边界；其余边界则均为绝热条件。

正常运行条件下稳态工况的温度分布示于图2，包括整体温度(含侧反射层)、燃料和氦气冷却剂通道各自的温度，以及典型燃料(位于堆芯中心)及其邻近冷却剂通道沿轴向(y方向)从入口到出口的温度分布。整体上，堆芯温度最高的区域位于堆芯中心且靠近氦气出口的位置。氦气冷却剂逐渐被燃料加热，并沿轴向流动方向单调上升。而燃料芯块沿轴向方向的温度分布则是由轴向功率分布和氦气冷却作用共同决定的：相比于轴向功率截断余弦分布的峰值位置，在氦气的冷却作用下，温度最高位置由中部移至靠近氦气出口的位置(图2d)。具体数值为：出口处氦气的最高温度940.1 ℃，燃料最高温度1 004.6 ℃、平均温度765.5 ℃。

a——整体计算域温度；b——冷却剂温度；c——燃料温度；d——典型燃料和冷却剂轴向温度

计算采用的网格单元总数约为50万。为确保网格划分的无关性，本文采用不同的网格数重复开展上述计算。选取燃料最高温度作为衡量指标，该指标随网格数的变化趋势示于图3。可看出，当网格单元数超过80万后，计算结果已趋于稳定，变化已不明显。因此综合考虑计算的精度和效率，本文后续瞬态计算中将采用单元总数为90万的网格划分。

图3 网格无关性分析结果Fig.3 Mesh-independence analysis result

3 事故瞬态分析

针对棱柱型高温气冷堆的堆芯和主回路设计特性，参考同类堆型事故分析清单并结合工程经验判断，本文选取3个具有代表性的事故工况开展计算分析。更系统全面的事故工况谱分析将在后续研究中进一步开展。如前文所述，以稳态的计算结果作为瞬态计算的初始值。事故分析的验收准则为：整个事故进程中的燃料芯块最高温度不超过燃料温度的设计限值1 620 ℃[15-16]。因与堆内其他构件/设备(如控制转鼓、压力容器等)以及与辐射剂量相关的限制准则目前尚未明确，故本文暂不予考虑。

3.1 PLOFC事故

该事故是由于主回路系统设备发生故障(如氦汽轮机卡轴)但未发生破口所引起的。分析时保守假设0时刻堆芯即丧失全部强迫循环流量，氦气流速瞬时降低为0。控制棒插入后成功实现停堆。之后仅依靠非能动余热排出系统带出堆芯热量。整个瞬态计算时间为10万s。

PLOFC事故主要计算结果如图4所示。由图4a可看出，整个事故进程中，燃料最高温度从开始就一直处于下降的趋势，温度最大值即为0时刻的1 004.6 ℃，始终低于燃料温度限值(1 620 ℃)。这是由于丧失氦气强迫循环流动后，氦气的对流换热作用消失，传热以导热为主，因此沿着燃料通道轴向的导热使得轴向温度的分布逐渐展平(图4b)，轴向温度峰值逐渐由靠近氦气出口位置移至中段部位(图4c)，从而有利于燃料通道最高温度的降低。另一方面，非能动余热排出系统不断从堆芯吸收并导出热量。此外，加上衰变热功率水平的持续下降，使得堆芯燃料的整体温度不断降低。约13 111 s时，余热排出系统的带热功率便开始超过衰变热功率(图4d)，最终导致堆芯燃料长期维持在一低于限值的安全状态。

a——燃料温度；b——典型燃料轴向温度；c——100 000 s时刻燃料温度；d——衰变热和余热排出系统功率

3.2 PLOFC+ATWS事故

该事故是在3.1节工况的基础上叠加了由于紧急停堆系统故障而导致的停堆失败。此后反应性的控制主要依靠堆芯(包括燃料和慢化剂)自身的温度反馈效应。根据核设计计算可得到堆芯总温度反馈系数。此外，非能动余热排出系统仍假设正常工作，把堆芯热量带入最终热阱。

燃料温度以及堆芯和非能动余热排出系统功率随时间的变化示于图5。可看出，0时刻氦气强迫循环流量丧失后，由于堆芯热量无法及时导出，燃料整体经历了一个快速升温过程，最高温度达到1 073.3 ℃。随后在堆芯负反应性温度反馈的作用下，反应堆功率快速下降。另一方面，随着燃料轴向导热将温度分布逐渐展平，加上余热排出系统不断从堆芯导出热量以及反应堆功率水平的降低，燃料温度也开始下降。当燃料温度下降到一定程度后，在堆芯负反应性温度反馈的作用下又会引起反应性的增加，从而导致反应堆功率和燃料温度的再次升高。之后该波动现象又重复出现数次，但波动幅度逐渐降低。最后反应堆功率长期维持在与余热排出功率相同的安全水平(约1%初始额定功率)。整个事故过程中燃料的最高温度为1 073.3 ℃，并未超过燃料的温度限值(1 620 ℃)。表明依靠反应堆自身的固有安全性，在该事故条件下仍能确保燃料结构的完整性，从而避免放射性裂变核素的大量释放。

图5 PLOFC+ATWS事故工况计算结果Fig.5 Computational result of PLOFC+ATWS accident

3.3 RIA+ATWS事故

该事故假设反应堆运行期间，由于某个控制转鼓发生故障转动而引入了正反应性。与此同时紧急停堆控制棒由于故障未能插入堆芯，导致反应堆紧急停堆失败。本文假设1个控制转鼓发生故障转动，其旋转速度为1(°)/s，共旋转了180°(对应该控制转鼓的中子吸收体完全移出)，因此整个反应性的引入时间为180 s。此后反应性的控制同样是依靠堆芯(包括燃料和慢化剂)的温度反馈效应。此外，由于停堆信号触发开启了汽轮机旁通阀，压气机将发生惰转，分析时保守假设0 s时刻的氦气强迫循环流量即降为0。非能动余热排出系统仍假设正常工作。

RIA+ATWS事故主要计算结果如图6所示。从图6可看出，由于0～180 s期间外部正反应性的不断引入，导致功率升高；但由于堆芯温度升高所引入的负反应性基本与控制转鼓匀速转动所引入的正反应性大致相等，因此0～180 s期间反应堆总的功率升高幅度较小，基本保持不变。另一方面，由于氦气强迫循环流量的丧失，导致热量无法及时导出，因此早期燃料通道整体温度有一快速上升过程，燃料通道的最高温度达到1 280.5 ℃。随后由于控制鼓故障转动导致的正反应性引入消失，在堆芯负反应性温度反馈的作用下，反应堆功率开始快速下降。当燃料温度下降到一定程度后，在堆芯负反应性温度反馈的作用下又会引起反应性的增加，从而导致反应堆功率和燃料温度的再次升高。之后该波动现象又重复出现数次，但波动幅度逐渐降低。此外，由于堆芯整体温度的升高，余热排出系统的功率也相应增大。最后长期阶段，反应堆功率大致维持在与余热排出功率相同的水平(约1.3%初始额定功率)，堆芯始终处于安全状态。整个事故进程中燃料通道的最高温度为1 280.5 ℃，未超过事故条件下的燃料温度限值(1 620 ℃)，且仍具有较大裕量。

图6 RIA+ATWS事故工况计算结果Fig.6 Computational result of RIA+ATWS accident

4 压力容器壁面发射率敏感性分析

由于非能动余热排出系统的详细设计细节尚未确定，因此前文计算中仅采用一个简化的对流换热边界条件来初步模拟空气的自然循环带热作用。但真实事故条件下，由于反应堆压力容器壁面的温度通常较高，因此辐射换热也是余热排出系统的一个重要带热机理。故本文进一步考虑辐射效应，即在前文自然对流换热的基础上，叠加压力容器壁面与非能动余热排出系统构成的反应堆腔室之间的辐射换热作用。计算时，采用漫灰表面辐射模型，并假设非能动余热排出系统构成的反应堆腔室与外部环境良好换热，其表面温度恒定为30 ℃，且该反应堆腔室的表面积远大于压力容器壁面面积。

以PLOFC事故为例，计算得到的不同表面发射率下燃料最高温度以及非能动余热排出系统的总带热功率(对流+辐射)示于图7。其中，图7c、d分别为事故早期(0 s)和长期(100 000 s)阶段非能动余热排出系统的对流和辐射各自带走的功率比例。图7表明：考虑辐射换热机理后，非能动余热排出系统的带热能力明显增强。发射率越大，增强效果越显著，对应的燃料最高温度越低，事故缓解效果越好。特别是事故早期阶段，反应堆压力容器壁面温度相对较高，且单位面积辐射换热功率与温度的4次方呈正比，因此此时辐射换热的带热效果不可忽略，对于早期阶段的反应堆热量导出具有重要意义；且发射率越大，辐射换热在非能动余热排出系统换热总功率中的占比也越大。事故后期阶段，随着反应堆压力容器壁面温度的降低，不同算例计算得到的非能动余热排出系统的总带热功率差别不大。但因为存在辐射换热，相较于不考虑辐射换热的情形，压力容器壁面的温度能降至更低的水平。由此可看出，非能动余热排出系统带热能力的准确模拟对于事故分析具有重要影响，因此在后续设计定型后还需对其进一步开展详细建模、模型校核与事故分析研究。

a——燃料最高温度；b——不同发射率的RCCS功率；c——事故初始阶段(0 s)；d——事故长期阶段(100 000 s)

反应性引入ATWS事故下的计算结果示于图8。与PLOFC事故工况下的结果相比，该类ATWS事故下的非能动余热排出系统换热能力对事故后果的影响趋势正好相反。即考虑辐射换热后，虽然非能动余热排出系统的换热能力增强，但长期稳定后的反应堆功率水平以及对应的燃料最高温度反而都会增大。这是因为考虑辐射效应后，堆芯温度下降幅度过大，由于温度负反馈效应而引入的正反应性也会相应增大。由于最终达到长期稳定状态后，反应堆功率基本与非能动余热排出系统的换热功率水平相匹配和平衡，因此当非能动余热排出系统换热能力增大时，长期阶段的堆芯平衡功率增大，所对应的燃料最高温度升高，但仍远低于其事故条件下的限值(1 620 ℃)。所以在工程设计和确定非能动余热排出系统的排热能力时，需同时兼顾和分析PLOFC以及ATWS类事故条件下的堆芯响应，根据不同事故工况下的燃料温度变化特性，确定和优化实际的排热功率需求。

a——燃料最高温度；b——反应堆功率

5 结论

本文基于通用CFD程序COMSOL，针对正在研发的棱柱型模块式高温气冷堆建立了详细三维模型，开发了点堆中子动力学模块，分别针对正常运行稳态工况和选取的3个典型事故瞬态工况开展了计算分析，并探讨了考虑辐射换热后不同发射率下非能动余热排出系统的换热效果，得到如下结论。

1) 在本文所选取的事故工况下，燃料最高温度均低于其安全限值(1 620 ℃)，且具有较大的裕量，因此均能保证堆芯燃料结构的完整性。其中，反应性引入且未能停堆工况的事故后果相对最为恶劣，早期燃料最高温度能升高到1 280.5 ℃，但仍满足事故验收准则；通过堆芯的负温度反应性反馈以及非能动余热排出系统的带热作用，最终仍可实现堆芯长期处于安全可控状态，表明该堆型设计具有较高的固有安全特性。

2) 辐射带热是非能动余热排出系统的一个重要换热机理，特别是在事故早期阶段，对于快速带走堆芯热量具有重要作用。其中，反应堆压力容器表面发射率是影响辐射换热的一个重要参数。后续还需对其进一步开展详细建模与校核。

3) PLOFC事故与ATWS事故工况下，非能动余热排出系统换热能力对燃料最高温度的影响趋势相反，在进行该系统设计时需同时对这两类事故开展分析。

本文开发的计算模型可作为后续事故分析模型开发的基础，从而应用于后续开展该堆型的详尽安全分析和评审。但该模型仍需进一步进行验证与校核，特别是针对非能动余热排出系统。