岳明楷，金 浩，刘桐蔚，李文杰，黄永忠，恽 迪，周进雄,*

(1.西安交通大学 航天航空学院 机械振动与强度国家重点实验室，陕西 西安 710049；2.中国核动力研究设计院，四川 成都 610065；3.西安交通大学 能源与动力学院，陕西 西安 710049)

核能以其持久、经济、安全和清洁等优势被人们认为是当前最具开发价值和发展潜力的新型能源，其持续发展不仅能保证国民经济的平稳发展，还对国家安全有重要的战略意义[1]。

热管冷却反应堆一回路系统不采用传统压水堆、液态金属冷却反应堆等使用冷却剂回路的布置方式，而采用热管组件将堆芯产生的热量传递到堆外的能量转换系统。20世纪60 年代，LANL( Los Alamos National Laboratory)提出了采用高效的热管作为导热元件的新型空间核反应堆设计概念[2]。LANL的研究者们进一步研究设计，逐渐形成燃料棒、热管相互间隔排布的热管冷却反应堆结构设计理念[3-4]。目前，已开展的热管冷却反应堆方案与关键技术研究已经证实了热管冷却反应堆具有固有安全特性、成熟技术基础的可借鉴性以及较小的比质量等特性[5-6]。2018年3月，Kilopower反应堆概念的地面原型堆项目KRUSTY成功完成了有核实验，标志着世界上第1个热管堆的成功运行[7]。热管堆便于小型化，核能具有可离线性、高功率密度、高可靠性、运行时间长，未来在偏远地区供能供热、科考供能、空间核动力等领域具有应用潜力[8-9]。

通过对热管堆的燃料性能进行分析，可促进对热管冷却反应堆特性的了解，进一步可帮助我们对固态堆芯设计参数进行优化和改进，从而提高热管堆的优势同时尽可能地减弱其热电转换效率低的缺点。

在热管堆工作寿期内固态堆芯长时间处于高温、辐照环境工况中，涉及间隙传热、材料蠕变等复杂的热力耦合行为。本文介绍Ross和Stoute经典间隙传热模型并对固态堆芯所用材料在高温、辐照环境中的蠕变、肿胀等力学行为进行简要说明。综合考虑间隙传热模型及极端工况下材料的力学行为，基于ABAQUS结合GAPCON、CREEP子程序的编制进行二次开发，建立固态堆芯等效二维截面热力耦合计算模型，并对5 a寿期内的工作性能进行数值模拟。基于参数化建模有限元方法对固态堆芯尺寸进行参数优化，通过深度神经网络建立堆芯最高温度、最大应力两个优化目标与设计参数的代理模型。基于代理模型，使用NSGA-Ⅱ算法获得pareto前沿解集，以获取两个优化目标权重相等情况下的最优堆芯设计。

1 固态堆芯燃料组件性能分析

热管反应堆固态堆芯燃料组件由堆芯基体和按一定规律间隔排布在其中的燃料棒、热管构成。其中1种典型的固态堆芯排布方式如图1所示，固态堆芯由外围中子反射层、反应性控制鼓、堆芯基体、热管、燃料棒等部分组成。燃料棒由柱状裂变燃料芯块和合金包壳封装组成。热管是1963 年LANL开发的传热元件，通过管内工质的相变过程的高潜热特性和毛细作用实现热量从蒸发段到冷凝段的传递，是一种高效的被动传热方式。对于如图1所示的热管冷却反应堆，燃料棒、热管以蜂窝状六边形排布方式间隔排布在堆芯基体中，这样的排布方式能最大程度上使燃料棒均匀地临近热管，使燃料棒释放热量高效地传递给热管。热管蒸发段在堆芯内与核燃料棒一起排布在固态堆芯基体中，冷凝段布置在热电转换系统中，反应堆工作时，核燃料内的裂变能量传递给金属热管的蒸发段，热管将热量从堆芯传导至能量转换系统转换为电能或将残余的热量排出至最终热阱[10-11]。

图1 典型热管冷却反应堆结构Fig.1 Typical form of heat-pipe-cooled reactor structure

热管冷却反应堆固态堆芯工作状态涉及传热、力学、核裂变反应物理等诸多领域的问题，各因素耦合相互影响，对这些问题的仿真和分析将有助于更好地设计固态堆芯。

1.1 含气隙结构的传热

传热问题在核工程中普遍存在，是固态堆芯重要研究问题。在热管冷却反应堆固态堆芯中，核燃料棒中的裂变燃料芯块发生核反应，产生的热量依次经过燃料包壳、堆芯基体传递到热管，经由热管内传热介质排出堆芯。整个传热路径经过各零件，零件装配所致的间隙对传热过程有明显的影响。传热路径示意图如图2所示。

图2 固态堆芯内传热示意图Fig.2 Schematic diagram of heat transfer in solid core

发生在固体内部的传热问题由傅里叶热传导方程控制：

(1)

式中：k为热导率；T为温度；Q为热源功率密度；ρ为密度；c为比热容；t为时间。

间隙是传热路径的重要组成部分，Ross-Stoute提出的间隙传热机制认为，透过间隙表面传递的热量由各表面的温度和等效间隙传热系数决定，间隙传热模型[12]如下：

qg=hgap(T1-T2)

(2)

式中：qg为间隙热通量；hgap为间隙传热系数；T1、T2为各间隙表面温度。其中，hgap由间隙内的传热机制决定，主要由三部分构成，包括间隙内气体的导热、间隙表面之间的辐射换热以及可能发生的间隙闭合部分的接触导热，可由下式表示：

hgap=hg+hs+hr

(3)

式中：hg为间隙气体传热系数；hs为间隙表面接触传热系数；hr为间隙表面辐射换热系数。

1) 间隙气体传热模型

hg与气隙的大小以及两表面的粗糙度有关，间隙表面附近气体温度并不等于固体表面温度，这种温度不连续的现象被称为温度跳跃。间隙表面的粗糙度也是影响间隙气体导热的重要因素。

气体传热系数计算模型[12]如下：

(4)

式中：kg为气体热导率；d为气隙宽度；Cr为粗糙度系数；R1、R2分别为间隙表面粗糙度；g1、g2为间隙表面跳跃距离。

跳跃距离的计算由Kennard模型[13]给出：

(5)

式中：amix为混合气体组分调节系数；Tg为混合气体平均温度；p为混合气体压力；n为混合气体中所包含的气体分子种类；Xi为第i种气体的分子份额；Mi为第i种气体的相对分子质量。

在间隙气体导热模型中，将气体假设为不会流动的固体进行分析，忽略在加热情况下可能发生流动而产生对流换热的作用。气体的热导率随着温度、气体成分和反应堆操作条件的变化而变化。例如对于由He、Xe、Kr等组成的混合气体，使用下式计算热导率：

kg=(kHe)x(kXe)y(kKr)z

(6)

式中，x、y、z为各气体组分的摩尔分数。

2) 间隙接触导热模型

当间隙表面局部可能发生接触时，接触区域对间隙传热起到了重要作用。由于接触表面并非理想光滑的，有一定的表面粗糙度，造成实际上只有接触区域的部分点实际挤压接触在一起，而且还存在一定的气隙。因此现有的接触导热模型均认为导热系数与接触压有关[12]：

(7)

3) 辐射换热模型

在气隙换热过程中，温度较高的情况下，辐射换热也是间隙换热的重要组成部分。依据波尔兹曼辐射换热定律有[12]：

(8)

(9)

式中：σ为Stephan-Boltzman常数；Fe为与两表面发射率相关的辐射函数；e1、e2为间隙表面辐射率。

1.2 堆芯材料的蠕变

对于反应堆部件材料，由于长期处于高温、高辐照的工况中，其形变成因包含了弹性变形、热膨胀、蠕变和可能存在裂变产物引起肿胀等多种因素[14-19]，材料总应变为：

εtotlp=εep+εtp+εcp+εsw

(10)

式中：εtotlp为总应变；εep为弹性应变；εtp为热膨胀应变；εcp为蠕变应变；εsw为肿胀应变。

对于UO2燃料芯块，εcp和εsw由Allison给出了经验模型[14]。芯块蠕变模型为：

(11)

燃料芯块的肿胀行为主要由固体裂变产物所导致的体积生长和裂变气体所在燃料芯块内部产生的气泡导致的体积膨大构成，芯块肿胀模型为：

(12)

(13)

e-0.016 2(2 800-T)-8.0×10-27F

(14)

对于并不发生裂变反应的包壳、基体等不锈钢材料，不考虑肿胀，仅考虑蠕变[15]，包壳蠕变模型为：

(15)

1.3 热力耦合问题

由于核反应堆工作环境处于高温、高辐照工况下，材料属性与温度、辐照等因素相关。堆芯因温度变化引起结构形变也使间隙特征改变，间隙热导率对间隙属性敏感，从而影响堆芯传热行为。因此对堆芯性能进行分析时，必须考虑多物理场相互耦合下温度、变形、辐照等多种因素的相互作用，在ABAQUS中热力耦合分析步进行逐步迭代计算。

在固态堆芯设计过程中，非线性是建模过程的难点。材料的肿胀、蠕变等机理导致材料非线性。间隙尺寸相对较小，随着寿期的变化，部分间隙处可能闭合，导致接触非线性。

2 固态堆芯热力耦合分析有限元模型设置

常规的奥氏体不锈钢和镍基合金都是较为合适的包壳结构材料，同时考虑反应堆整体的经济性，采用商业上技术比较成熟的316不锈钢作为包壳材料，同时在堆芯材料选择采用二氧化铀作为堆芯材料。固态堆芯由燃料棒、热管和基体组成，基体整体形状为六棱柱，其中有排布燃料棒和热管的通道，燃料棒和热管均为细长的圆柱体，固态堆芯结构示意图如图3所示，其截面尺寸远小于轴向方向尺寸，其几何特点为长细比(轴向方向尺寸与截面尺寸的比值)大、具有对称性。间隙尺寸至关重要，导致全尺寸的三维模型计算将耗费大量的时间成本，不适用于参数化生成大量数据样本。

图3 固态堆芯三维模型及横截面示意图Fig.3 3D model and cross-section diagram of solid core

出于计算成本的考虑，结合上述特点，对堆芯的横截面进行广义平面应变假设，并利用堆芯截面的对称性，选取代表性栅元代表堆芯，如图4所示。基于ABAQUS在堆芯设计尺寸允许范围内选取1组参数采用热力耦合分析步进行仿真计算。尺寸参数列于表1。计算堆芯在正常工况下5 a寿期内性能变化。在间隙等细节处布置了尺寸为0.4 mm的局部细化的网格。堆芯外表面设为绝热边界条件，采用简化的恒温热管模型，燃料芯块的热功率为107W/m3。通过编写子程序GAPCON加入了上述间隙传热模型，实现间隙传热行为的数值仿真，通过子程序CREEP加入了上述蠕变、肿胀模型，实现对材料蠕变和肿胀特性的模拟。计算参数列于表2[20-26]。

a——代表性栅元选取示意图；b——设计尺寸变量示意图图4 计算模型示意图Fig.4 Schematic diagram of model

表1 固态堆芯尺寸参数Table 1 Design size parameter of solid core

3 有限元结果分析

对于上述尺寸的固态堆芯结构，进行热力耦合分析。随着堆芯工作时间的增长，燃料芯块的蠕变肿胀逐渐加剧，各部件体积逐渐增大，因此，堆芯内部状态随着寿期的增长而逐渐加剧。通过有限元数值模拟获取其在正常工况下5 a寿期时的温度场分布和面内应力分布，如图5所示。

a——温度场仿真结果；b——面内应力仿真结果图5 寿期5 a时堆芯截面热力耦合仿真结果Fig.5 Simulation results of temperature and stress at operation period of 5 years

由图5a可知，堆芯整体温度分布规律显著，燃料芯块中心温度较高，呈现沿着燃料芯块径向环状下降的特点，除了堆芯最外围的燃料芯块，其余部分的分布规律相似，幅值相近，这是由于热管和燃料棒的分布具有周期性。堆芯最外围相邻的燃料芯块的温度显著高于其他位置的燃料芯块，这是由于该位置处附近的热管平均数量较少。

由图5b可看到堆芯面内应力的分布情况，据此能对芯块开裂进行部分解释。燃料芯块中心受到向内挤压，出现明显的压应力区域，而燃料芯块外围受拉，出现环向拉应力区域，这种应力特点是造成柱状的燃料芯块径向开裂的重要原因。形成这样应力分布的主要原因是，燃料芯块产生的热量要从外表面通过气隙向外传递，从而导致燃料芯块内部温度高，膨胀较多，内外相互挤压牵拉。除了燃料芯块的热膨胀，燃料芯块的肿胀、蠕变特性导致燃料体积随着时间逐渐增长，产生的蠕变应变、肿胀应变，导致燃料半径逐渐增大，燃料芯块与包壳之间的间隙逐步减小，逐渐缩小间隙，间隙传热系数增加，有利于间隙传热。不过此设计参数下的堆芯在设计的5 a寿期内，芯块和包壳的间隙并未闭合，因而未产生接触压力，燃料包壳、热管包壳和堆芯基体的截面内应力相较于燃料芯块并不显著。

表2 固态堆芯材料及计算参数Table 2 Material and calculation parameters of solid core

4 固态堆芯多目标参数优化

4.1 固态堆芯优化问题

固态堆芯关键性能指标主要为堆芯最高温度和堆芯内最大应力，其受到燃料棒、热管之间排布间距以及各间隙大小等设计参数的影响。对固态堆芯尺寸进行参数优化，以期保持燃料总量不变，使优化后堆芯能获得更低的堆内温度和Mises应力。

基于ABAQUS、isight、机器学习等搭建联合仿真智能优化平台进行优化问题求解，优化技术路线如图6所示。首先通过拉丁超立方抽样在设计空间内构建样本集，利用有限元求解堆芯温度场和应力场，构成样本数据库，共获得5 800组样本数据。然后基于样本数据库，将5 000组用于机器学习训练，800组用于测试，使用DNN建立固态堆芯设计参数与堆内最高温度、最大Mises应力两个关键性能指标的代理模型。基于代理模型，使用NSGA-Ⅱ遗传算法获得连续的pareto解集，以获取指定权重下的最优解。

图6 优化流程技术路线图Fig.6 Schematic diagram of optimization process

4.2 基于DNN构建代理模型

图7为堆内最大Mises应力和最高温度两个优化目标的DNN预测值和有限元计算结果的对比。代理模型的拟合优度R2均大于0.998，且预测平均相对误差均在0.8%以下，具有良好的预测准确性。

a——最大应力预测与有限元结果对比；b——最高温度预测与有限元结果对比图7 DNN代理模型预测与有限元计算结果对比Fig.7 Comparison of DNN surrogate model prediction and finite element results

4.3 多目标优化

利用NSGA-Ⅱ算法[27]对目标函数进行优化，图8为堆内最大Mises应力和最高温度的迭代曲线。经过5 000次迭代，搜集整个设计域中的非劣解集构成较为平滑的pareto前沿。图9为最终得到的pareto解集。利用理想点法求解给定权重下的最优解，获得了两个优化目标权重相等时的最优解(240，967)。

图8 应力、温度迭代过程Fig.8 Iterative processes of stress and temperature

图9 pareto解集与1∶1权重下的最优解Fig.9 pareto front solution and optimal solution with weight ratio of 1∶1

对等权重多目标优化的固态堆芯结构进行了有限元计算验证，热力耦合计算结果如图10所示，其分布规律和初始设计十分相似。在表3中将优化后的堆芯结构和初始设计的各项参数进行对比，优化后堆内最高温度为965 K，下降了8.44%，最大Mises应力为240 MPa，下降了34.43%。

图10 初始设计与优化设计结果对比Fig.10 Result comparison between initial design and optimal design

表3 初始设计与优化设计参数对比Table 3 Parameter comparison between initial design and optimal design

5 结论

基于热管冷却反应堆固态堆芯工作特点，结合商用软件ABAQUS和其子程序的二次开发，实现了对固态堆芯内含间隙的复杂传热行为的数值模拟以及肿胀、蠕变、热应变等不同材料力学行为共同作用下堆芯力学行为数值模拟，得出结论如下。

1) 堆芯整体温度分布具有显著规律，燃料芯块中心温度较高，呈现沿着燃料芯块径向环状下降的特点，此温度场分布规律导致燃料芯块中心受到向内的挤压，而燃料芯块外围受拉，出现明显的环状拉应力区域，这些拉应力的产生是造成圆柱状的燃料芯块径向开裂的重要原因。

2) 对固态堆芯尺寸设计参数进行优化，根据有限元大量计算获得的数据样本，使用DNN机器学习建立了尺寸参数和最高温度、最大Mises应力的代理模型，使用NSGA-Ⅱ算法获得了pareto前沿解集，以及在降低最高温度和减小最大Mises应力两个优化目标权重为1∶1的情况下最优解，相较于初始设计参数，最高温度下降了8.44%，最大Mises应力下降了34.43%，改善了固态堆芯的性能。