核反应堆系统多维度多物理场耦合有限元分析研究

巫英伟,贺亚男,章 静,田文喜,苏光辉,秋穗正

(西安交通大学 核科学与技术学院,陕西省先进核能技术重点实验室,陕西 西安 710049)

核能作为一种清洁、高效的能源形式,对满足人类不断增长的能源需求发挥着重要作用。核能技术的发展与安全运行密不可分,这需要对核反应堆进行深入的研究和精确的分析。核反应堆系统庞杂且运行环境严苛,存在多物理场耦合的复杂现象,其涉及的物理过程涵盖了中子输运、热传导、燃料性能演变、冷却剂流动换热等诸多方面,这导致反应堆安全运行面临诸多挑战性难题。采用多物理场耦合分析方法解决此类问题可提高核反应堆运行的安全性和经济性。

多物理场耦合分析方法通常采用高保真数值计算模型及先进求解技术,结合不同物理学模型,将各种物理过程和空间尺度进行有效耦合,并依托超算平台及大规模并行技术,构建多物理场耦合分析平台,以全面、系统地模拟核反应堆的行为。在过去十几年里,多物理场耦合分析方法在核工程领域取得了显著进展,并被广泛用于预测和模拟复杂系统的行为。美国开展了CASL[1]、NEAMS[2]和NRC-CRAB计划[3],搭建了VERA[4]和CRAB等堆芯多物理场耦合计算的软件集成平台;欧洲以SALOME为基础建立了堆芯多物理场耦合计算体系[5];韩国则建立了DeCART-MATRA软件包以实现多物理场耦合的计算[6]。国内方面,各科研院所与高校展开合作,开展了多物理场耦合高保真数值反应堆研发工作。中国原子能科学研究院、中国科学院与北京科技大学建立了数值反应堆原型系统CVR(China Virtual Reactor)[7]。北京应用物理与计算数学研究所同上海核工程研究设计院、西安交通大学针对国和一号(CAP1400)建立了数值反应堆[8]。中国核动力研究设计院同哈尔滨工程大学同样也开展了相关工作[9]。

有限元方法作为广泛应用于各个工程技术领域的数值分析方法,在固体导热、弹塑性力学、失稳变形、热工水力、裂纹扩展等复杂现象的数值模拟分析具有一定的优势。因此,不少机构基于有限元方法开发了燃料性能分析程序,如美国的BISON[10]和FALCON[11]、阿根廷的DIONISIO[12]、法国的ALCYONE[13]以及韩国的PRIME[14]等。除此之外,也有不少学者在商用有限元程序COMSOL和ABAQUS等的基础上进行二次开发,实现对燃料性能的模拟[15-17]。有限元分析凭借其灵活的建模、准确的数值求解等优势,将会是未来燃料性能分析的研究趋势,因此其在各多物理场耦合框架中也作为主流的分析工具。

传统系统分析程序与子通道程序框架固定、缺乏灵活的程序架构以及先进的数值算法和物理模型,而先进的多物理场耦合平台可以实现各物理模型的灵活植入,充分发挥多物理场耦合的优势,可弥补传统分析程序在软件结构及算法上的不足。因此,搭建多物理场耦合平台,针对耦合问题中的关键技术开展研究,对加快我国自主化多物理场耦合平台的开发进程具有重要意义。

西安交通大学核反应堆热工水力研究室(Nuclear THermal-hydraulic Laboratory, XJTU-NuTHeL)结合多物理场耦合平台与有限元方法的优势,开展了高保真模拟工具的开发,包括系统分析程序NUSAC(Nuclear System Analysis Code)和子通道分析程序FLARE(Fully-coupled transient code for Liquid-metal-cooled Advanced REactor)[18-19]、针对多种燃料的燃料性能分析程序BEEs[20-22],并构建了基于有限元方法的多物理场耦合模拟体系[23],实现了堆内关键现象的精细分析。

本文旨在介绍XJTU-NuTHeL为建立多物理场耦合平台所进行的热工流体计算模型的开发、燃料性能分析技术的研究以及多物理场耦合框架的建立等工作,为核反应堆多维度多物理场耦合高保真数值模拟分析提供有益参考。

1 热工流体计算模型开发

核反应堆热工水力分析是核反应堆设计、安审、运行的重要组成部分,因此热工流体的精确计算对于核能的发展具有重要意义。单相及两相流动是核反应堆正常运行和瞬态工况下非常重要的物理现象,准确模拟流动现象是核反应堆系统分析程序的关键。针对不同的尺度,堆芯热工水力分析可分为系统安全分析和子通道分析两种。随着数值计算方法和多物理场耦合平台的发展,开发模块化、更高精度的热工水力成为国际热点,如RELAP7[24]、SAM[25]、MOLTRES[26]。相较而言,当前国内在高阶、高精度、高保真的热工水力分析程序开发方面研究较少。XJTU-NuTHeL采用高精度间断有限元方法实现了单相及两相物理模型计算,开发了核反应堆系统安全分析程序NUSAC和子通道分析程序FLARE。

1.1 核反应堆系统两相流分析模型开发

在核反应堆系统流动传热模拟中,XJTU-NuTHeL基于先进多物理场耦合框架针对压水堆和先进反应堆开展了一系列研究,开发了核反应堆系统安全分析程序NUSAC。在单相流模型中,NUSAC采用高阶连续伽辽金有限元离散和SUPG(Streamline-Upwind-Petrov-Galerkin)与PSPG(Pressure Stabilizing-Petrov Galerkin)稳定性算法处理了不可压缩流体流动数值不稳定性问题[27]。其中单相流守恒方程为:

(1)

式中:ρ、t、u、p、T和z分别为密度、时间、流速、压力、温度和z轴位置;ψ为形函数;圆括号(F,ψ)代表函数F在求解域内的积分离散;g、f和De分别为重力加速度、阻力系数和水力直径;τPSPG和τSUPG分别为与扰动相关的稳定系数;H、cp和q‴分别为流体的焓、比定压热容和体积释热率。

在两相流模型中,NUSAC采用了一维两流体六方程模型和经过广泛验证的与RELAP5相同的本构模型,其中两相流守恒方程为:

(2)

式中:下标l、g、i、w分别代表液相、汽相、汽-液交界面、壁面;α和Г分别为空泡份额和单位体积相变流量;F、e、Q分别为阻力、比内能和热源;h*和h′分别为流体与壁面传质的相焓和流体与相间传质的相焓。

同时,在两相流数值方法中,NUSAC引入间断重构伽辽金数值离散方法和Roe-type对流项数值算法实现了高阶空间离散格式,显著减小了两相流的数值扩散。在单相和两相流数值求解过程中,NUSAC通过自动微分方法形成雅可比矩阵,缩短了程序开发周期,避免了人为错误构建雅可比矩阵引入的误差[18]。此外,NUSAC采用了PJFNK(Preconditioned Jacobian Free Newton-Krylov)全耦合求解方法,在保证精度的同时提高了JFNK方法的收敛性[19]。

XJTU-NuTHeL通过一系列国际通用的两相流测试基准题对NUSAC进行了广泛而全面的验证。通过水龙头问题验证了两相流方程的正确性,如图1所示。在此基础上,NUSAC模拟了Bartolomei实验的若干组工况[28],并与实验值和RELAP5计算值进行了对比,如图2所示,验证了NUSAC处理两相过冷沸腾工况的能力。NUSAC通过模拟FRIGG实验,验证了NUSAC处理饱和沸腾工况的能力,如图3所示。过冷沸腾工况和饱和沸腾工况含有大量本构关系式,如汽-液界面阻力、汽-液界面传热等,通过这两个基准实验也验证了NUSAC本构模型植入的正确性。

图1 水龙头问题验证[19]

图3 FRIGG饱和沸腾实验验证[19]

NUSAC也可用于反应堆系统稳态和瞬态计算,根据美国核能署主蒸汽管道断裂事故中的堆芯设计参数[29],模拟了稳态及降功率瞬态的响应,计算结果如图4所示。图4中:T_l为一回路流体温度,K;T_second为蒸汽发生器二次侧流体温度,K。可以看出,系统回路稳态温度分布符合预期,验证了NUSAC系统模拟的能力。同时,在稳态工况下引入降功率瞬态变化,在堆芯通道内冷却剂进出口温度和包壳最高温度分布响应如图5所示,NUSAC计算值与RELAP5计算值符合良好,验证了NUSAC的正确性。

a——计算对象示意图;b——冷却剂温度分布云图

a——冷却剂进出口温度;b——包壳最高温度

1.2 液态金属快堆子通道分析模型开发

液态金属快堆由于具有良好的固有安全性以及核燃料增殖能力,被认为是最有发展前景的第4代反应堆[30]。为了保证快堆的安全运行,有必要对其堆芯进行热工水力安全分析。子通道分析方法是核反应堆堆芯热工水力的重要分析手段,具有比系统分析方法更高的计算精度,而且相较于CFD分析方法,其计算资源消耗更小。因此,XJTU-NuTHeL建立了适用于液态金属快堆的子通道分析模型,开发了全耦合子通道瞬态分析程序FLARE。

FLARE的控制方程如下。

质量方程:

(3)

轴向动量方程:

(4)

横向动量方程:

(5)

能量方程:

(6)

通过上述四方程模型,FLARE在计算中考虑了通道间的横向交混、湍流交混以及横向导热。除此之外,FLARE具有钠、铅以及铅铋的物性状态函数,可以对上述3种冷却剂进行计算。同时,由于液体金属快堆具有特殊的绕丝结构,FLARE本构模型中还考虑了绕丝对通道几何、摩擦阻力以及湍流交混的影响,如图6所示。程序基于有限元与间断有限元方法相结合的混合方法对控制方程进行离散,并通过PJFNK算法对离散方程组进行全耦合求解。

图6 带绕丝燃料组件典型示意图

本文采用19棒束实验[31]对FLARE进行验证,图7为通道编号划分以及测量通道示意图。图8展示了FLARE对组件出口归一化温度的计算结果以及与SUBAC程序[32]和Pronghorn-SC程序[33]的计算结果对比。由图8可见,在高流量下轴向导热占主导,并由于几何结构不同,角、边和中心通道的出口温度存在着明显的差异。随着入口流量的减小,横向导热逐渐发挥作用,各通道的出口温度逐渐趋于一致。从图8可看到,FLARE的计算结果与实验值符合良好,相对误差在10%以内,并且与SUBAC以及Pronghorn-SC计算结果接近,说明FLARE能够对液态金属快堆堆芯内的热工水力参数进行准确计算。

图7 棒束实验通道编号示意图

a——高流量工况;b——中流量工况;c——低流量工况

2 燃料性能分析技术研究

核反应堆燃料作为反应堆内裂变场所以及阻挡裂变产物释放的第1道屏障,开展其服役性能分析与评估至关重要。传统的燃料性能分析程序大多采用1.5维的有限差分方法开发,受限于方法的适用性,此类程序的研究对象仅限于棒状燃料。考虑到有限元方法在固体力学方面的广泛应用及其建模对象的灵活性,XJTU-NuTHeL基于有限元开发了燃料性能分析程序BEEs。该程序不仅具备针对传统棒状燃料在稳、瞬态下的分析功能[34-35],还能够针对事故容错燃料(Accident Tolerant Fuels, ATF)[20]以及其他几何形状(如环形[22]、板状[21,36]等)的燃料开展多物理场耦合分析。

2.1 包覆颗粒弥散燃料性能分析

包覆颗粒弥散燃料是众多ATF燃料选型方案之一,结构包含三向同性TRISO(TRI-structural ISOtropic fuel)球形包覆颗粒、NITE-SiC(Nano-Infiltration and Transient Eutectic-phase Silicon Carbide)基体,如图9所示。TRISO颗粒包含UO2芯核和4层包覆层,各包覆层都有特殊的功能。第1层为疏松热解碳层(Buffer),能够为气态裂变产物提供储存空间并缓冲芯核的热膨胀与辐照肿胀;第2层为内致密热解碳层(Inner Pyrolytic Carbon, IPyC),主要是防止裂变产物对SiC层的侵蚀,并承受部分内压;第3层为化学气相沉积法制成的碳化硅层(Chemical Vapor Deposited Silicon Carbide, CVD-SiC),其阻挡固态裂变产物的能力较强,能够大幅度改善燃料滞留裂变产物的能力。此外,SiC的高导热性和比热容也能显著降低燃料的峰值温度。第4层为外致密热解碳层(Outer Pyrolytic Carbon, OPyC),与IPyC材质相同,主要作用是保护承压的SiC层,并且在SiC失效时阻挡裂变产物的释放。包覆颗粒弥散燃料具有良好的导热性能和耐辐照性能,并能有效阻止裂变产物的扩散。

图9 TRISO包覆颗粒弥散燃料形式

包覆颗粒弥散燃料结构极其复杂,高燃料体积份额的燃料芯块可包含上千个TRISO颗粒,而每个TRISO颗粒具有多层结构,给开展其燃料模拟带来了较大挑战性。此外,对包覆颗粒弥散燃料进行行为模拟需要处理复杂的非线性模型,同时考虑到辐照、热、力多场耦合作用,增加了分析的复杂性。XJTU-NuTHeL基于有限元方法开发了包覆颗粒弥散燃料性能分析程序,建立了燃料芯块-代表性体积元-TRISO颗粒的多尺度模拟方法(图10)。考虑到TRISO颗粒对NITE-SiC基体的变形影响可忽略,模拟过程中首先对移除颗粒的芯块总体进行计算并基于温度分布等结果划分子区域,再从子区域中选取代表性体积单元进行计算分析,对其中TRISO颗粒性能进行精细评估。

图10 燃料性能多尺度模拟方法示意图

XJTU-NuTHeL选用国际原子能机构(IAEA)针对高温气冷堆公开的基准题(IAEA Coordinated Research Program, CRP-6)[37]对程序进行了系列验证。其中基准题案例CASE-8是辐照过程中温度和内压周期变化条件下的TRISO燃料颗粒,因涉及辐照-热-力耦合而被国际上众多程序选为验证算例。图11展示了IPyC和CVD-SiC层内侧的切向应力结果,由图11可知,本程序的计算结果与国际上其他程序的计算结果[38-40]符合良好,验证了程序进行燃料模拟的能力。

a——IPyC层内侧的切向应力;b——CVD-SiC层内侧的切向应力

图12示出周向非均匀换热条件下包覆颗粒弥散燃料性能分析的典型结果。程序实现了燃料芯块-代表性体积元-TRISO颗粒的多尺度热力学三维精细模拟(图12a～c),获取了燃料各部分的温度和应力分布情况。燃料芯块总体温度受外表面周向非均匀换热影响而呈现“偏心”分布。由于热解碳材料随中子辐照而产生收缩,CVD-SiC层在IPyC和OPyC的辐照变形下受压,同时随着燃耗的加深在Buffer层和IPyC层之间会出现间隙。此外,程序实现了针对包覆颗粒弥散燃料的失效评估,得到燃料失效概率的变化情况(图12d)。

a——芯块总体温度;b——代表性体积元应力;c——TRISO应力;d——燃料失效概率

2.2 板状燃料元件性能分析

板状燃料元件具有结构紧凑、功率密度高、换热能力强等特点,在研究堆、船用动力堆等堆型具有广泛应用。同时,板状燃料元件包壳与芯体直接接触,对温度、应力应变演变及裂变气体释放等均有影响。根据燃料元件芯体的结构,可将板状燃料分类为单片式燃料元件和弥散型燃料元件。其中,单片式燃料元件采用整片的陶瓷燃料板作为燃料芯体,其芯体厚度一般较薄;弥散型燃料元件则采用陶瓷燃料颗粒弥散在金属基体中的复合材料作为燃料芯体。

针对板状燃料的上述特点,XJTU-NuTHeL基于有限元平台建立了燃料元件的材料物性和行为模型,采用全耦合方法实现了燃料元件性能分析[41]。为了准确模拟弥散型芯体复合结构,建立了一系列等效物性模型,包括等效热导率、等效热膨胀系数、等效杨氏模量、等效泊松比、等效本征应变等。图13示出寿期初和寿期末的燃料元件中心截面处的燃耗、温度和Von Mises应力分布。在寿期初,燃料元件温度分布是力学变形的主导因素,在功率分布和冷却剂自下而上的对流换热的综合影响下,元件上部温度较高,热膨胀剧烈,导致寿期初中上部应力较大。而在燃耗较深的情况下,长期辐照造成的本征应变在力学方面逐渐占据了主导作用,元件中部燃耗最深的部位产生了很大的辐照肿胀,导致寿期末燃料元件中部的应力显著上升。

a——寿期初;b——寿期末

在燃料元件及辐照-力-热耦合分析的基础上,开发锆合金的氧化和吸氢模型,实现了核反应堆恶劣环境下燃料包壳的氧化和吸氢模拟。针对1/4燃料元件模拟获得的代表性计算结果如图14所示。氧化层厚度分布主要受温度影响,在元件中部最厚。吸氢首先在氧化部位发生,然后沿浓度梯度和温度梯度扩散。最终燃料包壳边缘处氢浓度达到饱和固溶氢,并析出形成氢化物。

图14 燃料包壳氧化和吸氢模拟计算结果

3 多物理场耦合框架网格映射技术方案及评估

随着高精度数值算法、大规模并行计算及计算机软件技术的发展,实现反应堆不同系统部件精细建模及多物理场分析,获取高保真计算结果成为可能。当前,国内外已开发多个多物理场耦合平台用于反应堆多场耦合计算,如MOOSE[42-43]、SALOME[44-45]、MpCCI[46]、COMSOL[47]、LIME[48]和CVR[7,49]等。传统的直接耦合方式存在着数据交换效率偏低、耦合程序需采用相同网格等不足,需要建立多物理场耦合框架,采用准确高效的网格映射算法,实现反应堆多物理场耦合问题在统一框架下的求解。为此,XJTU-NuTHeL基于MOOSE/LIBMESH平台搭建了堆芯多物理场耦合框架,用于实现反应堆的高保真多物理场耦合计算。

网格映射技术主要包括两方面内容,一方面是外部程序与多物理场耦合框架之间的网格映射,另一方面是框架内程序之间的网格映射。在外部程序数据映射方面,本文采用了逐一映射的方式,即在框架内构建与外部程序采用相同网格划分形式的数据存储网格,之后将外部程序的计算结果根据点/单元逐个映射到网格中。网格映射方案如图15所示。

a——求解域;b——外部求解离散方式及值分布;c——有限元网格及值分布

框架内程序之间的网格映射采用网格映射算法实现。为了获得堆芯多物理场耦合计算中高效映射算法,本文采用单块板状燃料作为对象,研究了不同算法的传递误差和在千万级-百万级网格之间的网格映射效率。网格映射算法的效率评估结果列于表1。网格映射算法效率从高到低依次为:反距离权重插值[50]、网格函数传递[51]、L2映射[52]、最近邻传递。需要说明的是,为了避免计算资源对效率评估结果的影响,本研究中对每种算法所采用的计算资源均为该方法在超算节点上测试所需的最少计算资源,超算单个节点的核数为64。

表1 网格映射算法效率评估结果

在传递精度方面,同样选取板状燃料作为研究对象,假设传递的变量分别为余弦形式的函数,源网格为六面体网格(HEX8),目标网格为四面体网格(TET4),源网格与目标网格节点数分别为6 080和65 567。不同映射算法的传递误差列于表2。其中,L2误差为误差平方之和,H1误差为变量梯度误差平方之和。由表2可见,网格函数传递算法在该类网格上具有最小的传递误差,L2映射次之,最近邻传递的误差最大。综合考虑,网格函数传递算法能够兼顾精度和效率,在结构化网格中映射变量时更具优势。

表2 不同映射算法的传递误差

4 堆芯多物理场耦合框架搭建及典型应用

基于上述映射算法,XJTU-NuTHeL开发了堆芯多物理场耦合框架,并实现了板状燃料性能分析程序BEEs-Plates[21,36]、蒙特卡罗中子物理程序OpenMC[53]和一维流体域分析程序NUSAC[18,27]在框架内的集成,初步实现了开源CFD软件OpenFOAM[54]的耦合。除此之外,在该框架中还实现了由中国核动力研究设计院开发的热工水力程序CORTH[55]和中子学程序CORCA[56]的耦合计算。框架的结构如图16所示。

基于BEEs-Plates/OpenMC/NUSAC耦合程序[23]开展了JRR-3研究堆组件级核-热-力耦合计算。为了考虑燃耗对燃料性能的影响,在每一计算时刻将预先完成的燃耗计算中的结果映射到网格中,之后传递给BEEs-Plates进行燃料性能计算。程序采用Picard迭代方式进行求解,程序结构如图17所示,计算总时长为2.0×107s。对25%富集度的U3Si2-Al弥散燃料组件的中子物理计算结果如图18所示。可以看到,由于计算中在组件长度(z轴方向)和宽度方向(x轴方向)设置了反射边界条件,组件边缘区域的功率密度显著高于中心区域。由于235U消耗速度较快,功率越高的区域铀原子密度越小。

图17 BEEs-Plates/OpenMC/NUSAC耦合程序结构及数据传递的方向

图18 组件功率密度(a)和235U原子密度分布(b)

部分燃料性能参数的计算结果如图19所示。由图19a可看出,在组件包壳与连接部件接触区域容易发生应力集中现象。当组件的平均燃耗达到64.77 GW·d/tU时,组件最大应力达到329.96 MPa。从图19b中体积应变的分布结果可知,组件的最大体积应变为10%,约为最小体积应变的10倍。此外,与截面上功率的分布趋势一样,组件的体积应变分布也呈现出中心低、边缘高的趋势,说明温度对燃料的应变影响较大。

图19 组件应力(a)和体积应变(b)分布

5 结论

针对核反应堆高精度数值模拟及多场耦合分析需求,基于开源有限元计算框架MOOSE,XJTU-NuTHeL团队完成了相关热工流体计算模型开发,构建了系统分析程序NUSAC,其中一维单相有限元离散计算依靠相关稳定性算法,一维两相计算采用高精度间断有限元离散,实现了核反应堆系统相关稳态及瞬态计算;构建了子通道分析程序FLARE,可高精度准确预测液态金属快堆堆芯内的热工水力参数。开展了相关燃料性能分析技术研究,构建了燃料性能分析程序BEEs,实现了涉及包覆颗粒弥散、板状等多种类型、多类工况的燃料高维度精细化模拟。搭建了多物理场耦合框架,测试评估框架内不同物理场求解程序的网格映射与数据传递算法,实现了板状燃料组件级核-热-力耦合计算。

综上所述,针对核反应堆多维度多物理场耦合分析的有限元技术可以重构传统核反应堆热工水力分析程序,实现在数值算法、离散精度等方面的优化,为核反应堆多维度耦合问题奠定了基础;可以优化传统燃料性能分析技术,突破传统的棒状燃料元件的研究范围,充分应用到各类先进燃料概念的开发验证中,高效解决了核反应堆工程中的各类精细化模拟问题;可以为多物理场耦合问题提供合理统一的计算框架,通过在不同物理场程序间提供的网格映射与数据传递接口,充分考虑不同物理场之间的影响,有助于更好地理解和研究核反应堆复杂多场环境下的相关现象。