棱柱状高温堆堆芯高效热工计算模型开发及验证

郭娟娟，童 朴，孔德泰，何 昉，黄善仿，王 侃

（1.中核战略规划研究总院有限公司，北京 100048；2.清华大学，北京 100084）

高温气冷堆作为GIF 公布的六种第四代反应堆候选堆型之一[1]，具有固有安全性高、发电效率高、可用作核能制氢等特点。国际上广泛研究的高温堆有球床式和棱柱式两种形状的堆芯，均使用TRISO 燃料球，采用氦气作为冷却剂，石墨作为慢化剂[2]。其中，由清华大学研发的球床模块式高温气冷堆研究成熟，已在山东荣成建造并成功进行一次冷试，其堆芯采用球床式，燃料球由TRISO 颗粒组成，冷却剂氦气与燃料直接接触；棱柱状高温堆的堆芯则采用石墨棱柱作为基体，如图1 所示，基于TRISO 颗粒的圆柱状燃料芯块置于石墨基体的燃料通道中，氦气不和燃料直接接触，通过冷却剂通道将热量带走。

高温堆堆芯的热工计算方法根据堆芯形状有所不同，研究较为成熟的球床堆堆芯的冷却剂和燃料球直接接触，热工计算通常采用多孔介质法[3]，建模时不区分燃料和冷却剂；但是棱柱状高温堆堆芯的冷却剂和燃料棒不直接接触，横向温差大，多孔介质法不适用，需要对固体及气体分别建立方程，一般采用CFD 方法计算。但是CFD 方法计算量大，计算效率低，全堆尺度的模型计算时间过长，并且流体的三维NS 方程计算时间占比较大。由于棱柱状高温堆堆芯的冷却剂通道狭长，湍流强度大，单个冷却剂通道的径向温差小，因此可以将冷却剂通道的径向温度作平均，等效成为“一维流体”模型，而燃料和石墨基体的固体求解仍然采用三维导热模型计算。

本文针对棱柱状高温堆堆芯热工计算现有的问题，提出“一维流体、三维固体”热工模型，基于国际上通用的MOOSE 平台[4]，开发计算程序，计算MHTGR-350 基准题的单通道及组件级别算例，并与CFD 软件OpenFOAM[5]的计算结果进行对比，验证该模型的准确性和高效性。

1 基于MOOSE 的热工计算模型开发

本小节详细描述基于MOOSE 平台开发“一维流体、三维固体”简化热工模型的过程，包含流体传热、固体导热以及流固共轭传热模型。“一维流体”是指将圆形冷却剂通道的径向温度作平均，得到轴向温度随通道高度变化的分布曲线；“三维固体”是指燃料及石墨基体采用经典的三维导热模型求解，得到三维的固体温度分布；而流体和固体之间的共轭传热用边界条件简化处理，利用经验关系式求出流固接触面的湍流换热系数，并与流体的轴向温度分布作为三维固体导热求解的边界条件。

本文的热工模型用于棱柱状高温堆耦合计算，因此将程序命名为PHEASANT（Prismatic High Temperature Reactor Coupling System of Neutronics and Thermal-Fluid），包含中子物理等计算功能。

1.1 MOOSE 平台简介

Multiphysics Object Oriented Simulation Environment（MOOSE）是美国爱达荷国家实验室开发的面向对象的多物理耦合平台，提供功能齐全的的底层求解器，可自动实现并行计算、网格I/O 及多维度计算，大大简化了开发者的工作。MOOSE 框架如图2 所示。

求解器及libmesh 库位于最底层，提供底层算法、网格I/O、有限元离散库、并行计算等功能；MOOSE 框架位于中间层，提供方程的残差余量、对角元计算等功能；位于最顶层的计算模块或程序则需要用户开发，用户需要将求解问题对应的方程整理成为kernels、边界条件、初始条件等元素，完成模型开发，实现最终的计算功能。

1.2 流体传热模型

流体的轴向温度分布采用能量守恒方程（1）求解，

式中：cp——氦气比热容；

m˙——氦气质量流量；

Q˙(z)——z高度的功率密度，这三个变量均为已知量，待求的ΔT(z)为z高度处的流体温差，根据已知的进口流体温度，即可得到冷却剂通道轴向任意高度处的流体温度分布T(z) 。湍流换热系数采用经典的Dittus-Boelter 关系式求解，关系式的推导涉及到的参数有流体密度ρ(z) 、动力粘度μ(z) 和热导率λ(z) ，本文将以上三个参数拟合成为T(z) 的多项式。以下为经验关系式的推导过程：

其中：

v(z)——氦气在高度z处的流速；

m˙——氦气质量流量；

s——冷却剂通道的径向面积。

由公式（2）即可得到雷诺数Re(z)：

其中：

d——冷却剂通道的直径。

由公式（3）即可以得到普朗特数Pr（z）：

将公式（3）和（4）代入Dittus-Boelter 湍流换热关系式（5）中，得到Nu(z)：

根据公式（6）可以求得换热系数h（z）。接下来对换热系数进行修正，利用已知的热流密度分布q（z），求出壁面固体温度随高度的变化Tw（z），如公式（7）所示：

利用上式求出的壁面温度可以得到修正系数Ct(z)：

以上的h'(z) 和T(z) 是流固共轭传热边界所需要的两个变量，作为固体导热求解的边界条件。

流体计算模型属于利用已知条件和显式公式推导变量的过程，没有实际的数值计算，因此在MOOSE 计算模型中，公式（1）～（10）属于辅助变量系统，方程中的项是Auxkernel，T(z)，h'(z)变量是Auxvariable，在数值计算开始时程序已计算出辅助变量，作为数值求解的边界条件。

1.3 固体导热模型

“三维固体”模型采用经典的导热方程，如公式（11）所示。求解的固体区域包含燃料及石墨基体，本文求解单通道及三维组件级别的小尺度模型，暂不考虑辐射对传热的影响。

Q˙——燃料区域的功率。

在MOOSE 中，将方程（11）按照Galerkin有限元离散方法处理为弱解形式的残差方程[6]，得到公式（12）：

1.4 流固边界共轭传热模型

热工计算中，流体和固体之间的共轭传热十分重要，结合1.2 节和1.3 节的推导过程，固体壁面处的热流密度和换热系数及流体温度之间的关系可以表示为：

式中：q——热流密度；

h——换热系数h'(z) 。

结合公式（13）和公式（14），公式（12）变形为：

公式（15）即为简化模型的核心方程，实现了流体和固体之间的共轭传热。

1.5 PHEASANT 程序实现

PHEASANT 程序的框架如图3 所示，最底层的MOOSE 框架支撑程序实现相应计算功能，经过有限元离散的方程中的kernels、BCs 等分别归类到PHEASANT 程序相应的块中，增加或改变问题方程只需要在相应块中增加或更改类即可。因此，不同方程若含有相同的kernels 或者边界条件类型等，可以通用。

本文中，公式（15）中的kernels 和边界条件BCs 分别在PHEASANT 相应的块中开发好，材料的密度、导热系数等参数在material 块中给定，选择Newton 求解器并行计算。

2 基准题验证

2.1 MHTGR-350 基准题简介

MHTGR-350 MW 基准题[7]的原型堆是MHTGR-350 MW 概念堆，由美国、韩国、德国等国家的高校或国家实验室共同参与开发，致力于为棱柱状高温堆堆芯模拟计算建立相同的计算条件，从而比较不同程序之间的计算结果。该基准题热功率350 MW，属于棱柱状高温气冷堆，如图4 所示，氦气作为冷却剂，石墨作为慢化剂，每个组件都是六棱柱状，由TRISO颗粒制成的棒状燃料插入燃料通道中，氦气通过冷却剂通道将热量带走。

本文计算基准题的热工稳态算例，采用给定的功率分布及流体固体物性参数。已知的功率分布如图5 所示，是基准题的物理临界算例的功率分布结果[8]，该算例为平衡循环寿期末的临界问题，仍有部分控制棒插入，因此功率曲线不完全满足余弦分布。

计算所需的边界条件列于表1，是基准题热工稳态算例设定的参数。氦气、燃料和石墨的物性参数分别列于表2 和表3 压力取6.39 MPa 后拟合成温度T（单位为K）的多项式。表3 中的石墨物性参数采用全堆平均温度所对应的值[7]。

表1 MHTGR-350 MW 基准题热工计算主要参数[7]Table 1 Thermal-hydraulics properties in the steady-state calculation of MHTGR-350 MW benchmark

表2 MHTGR-350 MW 基准题的氦气物性参数[7]Table 2 Helium properties of MHTGR-350 MW benchmark[7]

表3 MHTGR-350 MW 基准题燃料及石墨物性参数[7]Table 3 Fuel and graphite properties of MHTGR-350 MW benchmark[7]

2.2 单通道计算及验证

本小节以CFD 软件OpenFOAM 为基准，采用单通道级别算例验证本文提出的“一维流体、三维固体”热工模型的准确性。

在堆芯中，每个燃料柱均为六棱柱状，径向的每个冷却剂通道周围有六个燃料棒，以此为最小单元规律排列。为简化计算、便于结果对比，以冷却剂通道为中心将此单元按照面积等效原则简化成为圆形通道模型，如图6 所示。

等效之后的单通道模型以冷却剂通道为中心，直接与中间石墨环相接触，最外环是燃料，等效之后的轴向高度不变，仍是活性区高度793 cm，径向尺寸列于表4 中。

表4 MHTGR-350 MW 基准题等效单通道模型尺寸Table 4 Single channel model size of MHTGR-350 MW benchmark

分别采用CFD 程序OpenFOAM 和本文的PHEASANT 程序计算此单通道算例。其中OpenFOAM 不采用简化模型，计算完整的三维流体传热、三维固体导热及流固共轭传热，直到满足收敛准则。本文的PHEASANT 程序则利用开源优势采用本文开发的“一维流体、三维固体”简化热工模型。

两个程序计算的单通道模型均采用ANSYS软件划分网格，由于本文的简化热工模型不对流体区域做数值计算，因此PHEASANT 计算的网格模型没有流体区域，如图 7 所示，而OpenFOAM 的计算模型有流体区域。

两个程序计算时，除流固接触面的边界条件不同外，其余边界条件、物性参数等条件完全相同。冷却剂和燃料温度的计算结果如图8 所示。

图8 中的冷却剂及燃料轴向温度均经过径向平均，两个程序的计算曲线基本吻合，冷却剂温度曲线的平均偏差为1.13%，最大偏差为2.30%；燃料温度的平均偏差为0.44%，最大偏差为1.01%，拟合曲线的最大偏差均小于工程应用的5%不确定度。因此，单通道的计算结果验证了本文提出的“一维流体、三维固体”热工模型的准确性。

由于单通道模型的网格量较少，两程序计算花费的时间差距较小，因此采用组件模型进一步验证该简化模型的高效性。

2.3 组件计算及验证

图9 中，红色圆点代表燃料棒通道，蓝色圆点代表氦气冷却剂通道，六边形的六个绿色顶点代表可燃毒物通道，详细的组件尺寸等参数列于表5 中。

表5 MHTGR-350 MW 基准题组件模型尺寸[7]Table 5 Size information of the MHTGR-350 MW benchmark assembly model[7]

对1/12 组件模型建模画网格，网格模型如图10 所示，其中PHEASANT 程序采用的网格不包含冷却剂通道。

在PHEASANT 和OpenFOAM 计算中，模型的三个侧面和上下底面均设置为对称边界条件，即绝热边界。冷却剂通道与石墨接触面的条件设置同两程序单通道模型的计算条件，计算结果如图11 和图12 所示。

图11 定性展示了两程序计算的堆芯出口处的径向温度分布情况，除PHEASANT 没有冷却剂通道外，其余位置的温度分布基本相同。图12 中的轴向温度经过径向平均，两曲线基本吻合，平均偏差为0.96%，最大偏差为1.50%。

反应堆设计重点关注平均温度及最值温度参数，两程序计算得出的冷却剂出口平均温度、三维组件模型中燃料及石墨平均温度、燃料最高温度如表6 所示。

表6 MHTGR-350 MW 基准题组件模型温度结果Table 6 Temperature results of MHTGR-350 MW benchmark assembly model calculation

根据能量守恒，两程序计算得出的出口处冷却剂平均温度相等，均为687.1 ℃；两程序计算出的燃料平均温度、最高温度及石墨平均温度偏差分别为0.72%、0.28%和0.57%，最大不超过1%。对两程序计算组件级别算例花费的时间进行比较，结果如表7 所示。

表7 MHTGR-350 MW 基准题组件算例计算时间对比Table 7 Calculation time of the MHTGR-350 MW benchmark assembly case

PHEASANT 及OpenFOAM 均使用20 核并行计算，分别花费3.5 min 和43 min。

由于PHEASANT 计算的组件模型不包含流体区域，因此网格量比OpenFOAM 少将近1/3，计算中也不需要对流体区域做数值计算；此外，PHEASANT 程序是基于MOOSE 平台开发，计算采用Newton 算法求解方程组，相比OpenFOAM 采用的SIMPLE 算法计算效率相对更高。因此采用“一维流体、三维固体”简化模型的PHEASANT 的计算时间比OpenFOAM节省了90%以上。

3 结论

本文提出了“一维流体、三维固体”简化热工模型，用于棱柱状高温堆堆芯热工计算。基于多物理耦合平台MOOSE 开发了PHEASANT 程序，实现了简化热工模型的计算功能。

采用PHEASANT 程序计算单通道及组件级别算例，并与CFD 程序OpenFOAM 进行比较。两者计算结果的平均偏差1%，最大偏差不超过3%，且由于PHEASANT 程序所用网格量减少1/3，计算时间减少了90%以上。对比结果充分验证了简化热工模型的准确性和高效性。