球床式高温气冷堆堆芯持续性地卸出、装入球形燃料元件，燃料球的堆积、流动会影响反应堆堆芯中的中子通量分布及冷却剂流动状态，进而影响反应堆的功率和温度。燃料球间的传热机理和温度状态也会间接通过核反应截面、流体物性等因素影响球床的功率和温度分布。从经济性和安全性的角度出发，深入理解堆芯球床的流动和传热规律是高温气冷堆(HTGR)工程研究的重要组成内容。
球流的特殊性在于高温气冷堆中燃料球和石墨球随机堆积在堆芯区域，由于连续装料的设计，球床内燃料球按恒定的速率进行装卸，且装卸速率极慢，呈现为一种缓慢稠密循环流动形式。并且，这种流动极为特殊，与快速稀疏流不同，后者常被认为满足气体动理论。反应堆球床流动为宏观稠密复杂系统，球堆积紧密且相邻球之间的接触力持续时间相对长得多，但又会发生间断，球流运动难以用气体动理论描述、也难借鉴已有的连续介质的流动规律。
球流的复杂性在于球床堆芯堆积紧密，球间接触时间长、运动存在间歇性、球床堆积率较高，使得装卸料过程中球形元件的运动与多物理耦合较为复杂，存在丰富的非线性现象。然而目前并没有统一的方程描述宏观球流的运动，只能通过广泛的研究尽可能地发现与总结各可能因素对球流运动的影响，以深入理解球流运动规律。
一般，堆芯流动的影响因素主要包括(见图 1)：几何参数(床尺寸与球直径之比、球床壁面结构、出口大小、球床底部倾角等)、物性参数(球密度、摩擦系数、恢复系数、材料的弹性模量、Poisson比等)、运行条件(循环速率、加载比等)。实验场景复杂，较难完全复刻实验工况得到所有燃料球数据；实验周期长，限制了通过大量的实验研究得到足够数据以揭示球流运动的普适规律。这些导致球流问题的复杂性仍未得到充分的解释，需要借助更多的研究方法，计算机模拟的快速发展使之成为燃料球流动研究的主要方式之一。
1 球流研究
1.1 球流模型
从方法研究的角度，研究者尝试以上模型与DEM进行耦合。例如将计算流体动力学(CFD)与DEM进行结合以实现球床中流场的精确计算[22]，将DEM和格子Boltzmann(LBM)耦合以研究球床中的氦气流动[23]等。这2种方法均属于Euler-Lagrange混合方法，更常见于流化床、化工、采矿等领域[12]。光滑粒子流体动力学(SPH)与DEM的耦合方法则完全基于Lagrange思想[24]。为更好地与燃料球内部力学特性的研究耦合，本文团队也尝试将有限元方法(FEM)与DEM相结合，提出了DEFEM模型，可用于精细化分析燃料球的运动、变形和传热[25-26]。DEFEM方法不仅能够准确模拟燃料球的变形及燃料球内部的温度分布，也具有较好的并行性，可以借助并行算法来提高计算效率。
1.2 燃料球流动参数
开展DEM数值仿真时，燃料球物性参数的选取可能会对结果产生影响。例如，石墨的摩擦系数在不同气体氛围、不同温度条件下会发生显著变化[27-31]；较小的杨氏模量值对燃料球碰撞时间和形变量有显著影响，但却对最终的混合效果影响不大[32]；不同恢复系数对球流也存在影响[33]。Yurata等[34]基于实验结果给出了不同材料恢复系数的参数表达式，为后续DEM数值研究提供了参考数据。Lima等[35]基于实验分析了DEM模拟中恢复系数和摩擦系数的取值，提出特别要考虑材料性质及燃料球碰撞的尺度。Wang等[36]和Li等[37]均分析了恢复系数和摩擦系数分别对静态堆积球床和HTGR球流的影响，指出燃料球相对壁面的摩擦系数对流型的影响很大。对聚变反应堆中摩擦系数对球床结构的影响也有专门讨论[38-39]。
1.3 燃料球流动设计
球床的几何形状是数值研究和工程设计的重点关注内容。例如，Rycroft等[40]发现卸料口倾角越大，燃料球在球床中的总体滞留时间越短，但燃料球结晶化现象可能导致部分燃料球的滞留时间偏离均值。Khane等[41]从流型、速度分布的角度讨论了卸料口倾角的影响，指出速度曲线的径向分布呈现明显的下凹抛物线形状，且卸料口倾角越大，速度曲线下凹的现象越明显。Li等[42]发现壁面沟槽结构能够减少球床中近壁面燃料球的结晶化问题。Wu等[43]对比分析了多种卸料口漏斗形状，发现传统的锥形卸料口可能并非最优设计，外凸形的卸料口可以进一步提高流动的均匀性；还分析了球床高度、卸料口倾角和直径对球流的影响，发现卸料口倾角对燃料球流动性的影响更大，且更大的卸料口直径能够加强流动的均匀性[44]。Li等[45]对比分析了多种三维壁面结构对球流的影响，发现带沟槽的壁面结构会减缓近壁面区域的燃料球流速，但能够加强上下层燃料球的混合，有助于展平反应堆功率分布。Reger等[46]开展了HTR-PM全尺寸堆芯球流模拟，对比了不同壁面挖孔结构对流动的影响。该研究发现，壁面挖孔结构使得不同高度的近壁面孔隙率更加一致，且有助于防止燃料球的结晶化，如图 2所示。
燃料球的形状和尺寸也是数值研究的重点关注内容。Lu等[47]和Feng[48]对这个领域进行了总结性回顾和评述。Cui等[49-50]结合数值模拟和实验讨论了方形、球形燃料混合对卸料的影响，结果显示：过多的方形燃料元件会阻碍球床卸料，且增大中心流道与近壁面燃料球的速度差，导致更加严重的燃料球滞留现象。Ge等[51]讨论了椭球形、球形燃料混合对球床卸料的影响，发现当椭燃料球的长径比小于或等于2时，其体积占比增加可以促进球床卸料，而更大长径比的椭燃料球则会对球床卸料起阻碍作用。Kim等[52]讨论了不同尺寸燃料球混合对聚变反应堆球床堆积的影响，同样发现燃料球相对尺寸是一个重要因素。Wu等[53]也讨论了不同尺寸燃料球混合对球床卸料和孔隙率分布的影响，发现当装载比与其数量比一致时，球流在循环过程中的分区比例保持稳定。Pupeschi等[54]结合实验讨论了燃料球尺寸对球床的影响，通过改变床层高度与球尺寸比(H/d)进行了参数化研究，为此还采用了不同大小的商用单一粒径氧化锆球，系统评估了床层高度、球尺寸和球材料的影响。Feng等[55]讨论了燃料球尺寸以及摩擦系数对球床结构的影响，指出较低的球间摩擦系数和较宽的球尺寸分布会带来更高的球床堆积分数；随着球床宽度的增加，球床平均堆积分数逐渐增大，同时壁面效应的影响逐渐减弱。
综上，燃料球流动规律和行为的影响因素很多，而且往往与传热因素、物理因素发生深度关联。因此，燃料球流动需要根据具体工程设计中关注的重点进行针对性设计和优化。
2 燃料球传热研究
除了HTGR，堆积床中的燃料球传热问题也广泛出现在各种工程如流化床反应器[56]、聚变反应堆[57]等。有关球床内燃料球传热的数值模型可分为两大类：一类是基于Euler-Euler模型的连续论方法[58]，如双流体模型(TFM)[59]和多孔介质模型[60]；另一类则基于Euler-Lagrange框架的混合方法，如CFD-DEM、LBM-DEM。连续论方法将燃料球与流体均视为连续体，这一假设能够降低网格划分难度，同时简化复杂的流固耦合问题。因此，通常认为连续论方法具有较好的计算性能。相比连续论方法，CFD-DEM等混合方法更加契合真实的物理机制，能够有效揭示燃料球尺度的各种物理信息[61]，但该类方法依赖于准确的流固耦合模型，同时可能导致巨大的计算量。总体而言，目前已有非常多基于DEM的燃料球传热模型[62]和实验研究[62]，但是大多数模型并不适用于球床中的燃料球传热问题。一方面，HTGR球床中的燃料球最高温度可达1 600 ℃[63]，此时热辐射占主导地位，但大部分燃料球传热模型并没有充分在燃料球尺度考虑辐射换热。另一方面，HTGR中的燃料球数目较大，如HTR-PM包含约420 000个燃料球，直接套用已有的模型方法可能会带来不可接受的计算量，因此还需要对计算模型进行优化。
对于堆积球床，大部分工程研究基于多孔介质模型[60]开展全堆芯尺寸的流动传热模拟。这是因为球床流场中存在复杂的固体边界条件，很难基于燃料球解析的CFD-DEM方法开展大尺度的流动传热模拟。Sharma等[57]结合实验数据，对比分析了多孔介质模型和CFD-DEM模型在矩形堆积床热工水力特性研究中的效果，发现2种方法具有一致性。Baggemann等[64]基于多孔介质模型和MGT-3D程序分别模拟了德国SANA实验中的温度场分布，验证了程序的准确性。Novak等[65-66]同样利用SANA实验结果验证了新开发的数值程序，后续基于该程序开展了球床式高温熔盐堆中的多物理场数值研究。Zou等[67]将多孔介质模型和燃料球内部导热模型相结合，开发了球床堆中的多尺度传热模型。Nouri-Borujerdi等[68]基于改进的多孔介质模型开发了名为THPP的计算程序，采用IAEA提供的HTGR基准题进行了验证。
2.1 燃料球导热模型
导热、对流、辐射是燃料球热传导的3种方式，而基于DEM的燃料球尺度导热模型则是其中被研究最多的。
2.2.1 燃料球内导热
Luo等[69]从解析的尺度对传统的接触热阻模型进行了改进，提出了一种新型的燃料球接触面附近的导热热阻模型，该模型基于温度场和应力场的类比，使得燃料球表面的热流密度连续变化，消除了接触面边缘温度的奇异性问题，更加具有物理意义。该新型热阻模型与实验结果吻合较好。将该模型运用于HTTU(hochtemperatur-test-Urananlage) 堆的导热计算，得到了球床的温度分布和有效导热系数的分布，并可用于分析球床的动态性能。
式(1)、(2)、(3)分别给出了热-力比拟求解的温度方程、比拟边界条件、热阻解。
针对燃料球表面两接触点之间的导热热阻计算依赖于数值解，理论和解析的研究不够充分的问题，Luo等[71]继续将广义热阻模型进一步推广到三维的情况，从严格的数学推导出发，基于广义热阻理论推导出了燃料球表面任意2个接触点之间的导热热阻的解析解。该模型能够反映燃料球之间的接触面的大小和相对位置关系，能够为球床的导热问题提供参考，可依据该模型获得更精确的温度分布，并获得基本堆积方式的有效导热系数。
由于燃料球二氧化铀燃料核芯的温度不能超过限值是反应堆设计的准则之一，对核芯温度计算的建模需要结合堆芯物理热工耦合的瞬态模拟。相关研究结果可见文[79]，采用液态氟化盐作为冷却剂、内嵌包裹燃料颗粒的球形燃料元件的氟盐冷却球床堆，并进行堆芯瞬态分析。
2.2.2 燃料球间导热
HTGR球床属于准静态稠密堆积，聚变堆中的球床属于静态堆积，因此静态、直接接触模型最适合用于描述这2个领域的燃料球传热问题。典型的接触导热模型一般为
其中：qc, ji代表单位时间第j个燃料球向第i个燃料球传导的热量，Tp, i和Tp, j分别为第i、j个燃料球的温度，Hcond为接触热导率，单位为W/K。燃料球导热模型的重点在于如何计算接触热导率。在大部分模型中，接触热导率与燃料球的固体热导率、接触半径、碰撞时间等因素有关。Batchelor等[86]基于Hertz理论提出了经典的燃料球静态接触导热模型，该模型是后续诸多研究的基础。
2.2 燃料球对流模型
球与流体的对流换热是CFD-DEM研究的重要内容之一。具体而言，当流体流过燃料球表面或者壁面时，流体相和固体相在温度梯度的作用下会产生热量交换，一般采用对流换热系数来计算换热量：
其中，hconv, i为第i个燃料球的对流换热系数，Ai为第i个燃料球或壁面的表面积，Tf为当前位置处的流体温度，Ts为固体相的局部温度。
对流换热系数与诸多因素如流场流速、流体热导率、球床孔隙率、固体热导率、固体比热容等存在关联。基于实验或者如图 4的直接数值模拟(DNS)，大量研究[56-61]给出了对流换热系数与这些参数的关系式。Wakao等[89]和Gunn等[90]各自提出了随机堆积床中Nu数的拟合式，这也是最为经典的2个拟合式。Wakao等的拟合式默认孔隙率为定值(约等于0.4)，而Gunn等的拟合式则考虑了孔隙率的变化。后续Tavassoli等[91]在这2个拟合式的基础上进行了诸多修正[92-93]。Chen等[94]和Chang等[95]基于LBM和浸没边界法(IBM)的高精度模拟，拟合了特殊形状燃料元件的对流换热系数。需要注意的是，这里提到的拟合式适用范围各有不同，在球床传热计算中需要根据实际的流场参数选择合适的模型。
2.3 燃料球辐射模型
热辐射是HTGR中的燃料球传热的主要方式，重要程度甚至超过燃料球之间的接触导热。例如Wu等[96]指出，当温度大于550 ℃时，HTGR球床中的燃料球热辐射不可忽略。Forgber等[97]则发现高Pe下辐射换热量是接触导热的100倍以上。热辐射的研究包括连续性理论和离散化理论。连续性理论一般通过各种方法求解辐射传输方程(RTE)来得到辐射换热量[98]。Shi等[99]针对RTE提出了新的求解方式，通过引入UGKP格式将微观和宏观物理量联系在一起。HTGR和聚变堆中的稠密球床更适合采用离散化方法来计算热辐射。Cheng等[100]基于Voronoi单元体计算燃料球间的辐射换热，Wu等[101]在其辐射实验数据[102]的基础上提出了长程和短程辐射换热模型，以及燃料球尺度模型[101]、均匀化连续模型[103]、神经网络模型[104]、亚颗粒辐射模型(SCM)[105]及其改进型(SCRC)[106]，以及梯度提升决策树的机器学习模型[107]等多种燃料球热辐射计算方法。最近，谢宛均[108]将相关性函数与集成学习算法结合，将角系数函数分解为距离和遮挡函数，距离函数由理论推导得出，遮挡函数由集成学习模型预测，进一步提高了预测精度。以上这些工作，为球床辐射换热研究做出了重要贡献(见图 5)。
式(6)为经典的灰体辐射模型：
3 数值实现
3.1 GPU-DEM模型
当前，并行计算成为了数值研究的一个重要研究方向。传统的并行计算一般依托于开源多处理器(open multi-processing，OpenMP)或者消息传输接口(message passing interface，MPI)的架构[114]，采用CPU多核实现并行加速。随着图形处理器(GPU)硬件和人工智能技术的迅猛发展，GPU逐渐被用于科学计算领域。
3.2 GRAPE软件
基于DEM-GPU并行加速技术，本文团队提出了一套三维数值模型，并自主开发了对应的C++数值程序，取名为“面向高温气冷堆球床燃料球元件的GPU加速反应堆芯分析代码”(GPU-based reactor analysis code for pebble elements in HTGR，GRAPE)，如图 6所示。该程序针对HTR-PM设计优化，能够低成本、高效率地开展大规模燃料球流动、传热数值研究，最多支持千万量级燃料球的模拟。GRAPE可以在配有GPU的个人主机上实现核电站堆芯的数值计算。例如，模拟HTR-PM球床中42万个燃料球循环卸料100 s，时间步长取10-5 s，GRAPE在一块RTX 3090显卡上只需要运行11 h即可完成计算(作为对比，若采用12核的CPU需要大概20 d)。
4 AI技术的应用
4.1 球流研究中的机器学习应用
随着人工智能的发展，深度学习方法可以学习数据的内在规律，建立输入与输出间的广义映射关系。
近年来，颗粒物质相关的研究领域也逐步开始探索深度学习的应用。在燃料球流动领域，早期主要是利用人工神经网络对DEM参数进行校正(此时使用的神经网络层数较浅，不属于深度学习范畴)，而后逐渐使用更深的人工神经网络或者更复杂的神经网络如卷积神经网络对燃料球流动特性进行探索与预测，如图 7所示。
4.1.1 卷积神经网络在球流中的应用
目前球流研究中卷积神经网络主要用于以下几方面。
1) 模型参数校正。
2) 预测球流运动特性。
类似地，深度全连接网络也被应用于预测黏性燃料球流动在球床内的滞留时间[132]，采用数据驱动方法建立输入数据和输出数据之间的广义映射关系，其中输入数据通常是一系列变量如物理参数、燃料球流动位置坐标等。
3) 基于图像获取燃料球流动信息。
随着更先进的神经网络的出现，深度卷积神经网络在图像识别、缺陷检测、超分重建等诸多图像任务中取得了巨大成功，近年来开始有研究者对燃料球物质的图像驱动方法进行探索。例如，使用卷积神经网络(convolutional neural network, CNN) 从单张照片中自动提取燃料球的形状和尺寸分布[133]；根据燃料球的堆积状态自动提取图像特征，以实现较高准确度的有效载荷分布预测[134]；结合CNN和数字图像处理提取和评估粒子形状[135]；预测自由卸料漏斗流的卸出时间和剩余燃料球比例[136]；使用三维CNN输出不同燃料球形状以及堆积高度的填充率与质量流率[137]；基于图像获取岩石渗透率信息[138]；对混合流动燃料球体系的混合燃料球径和形状信息进行分类[139]；开发完全自动化的能够识别并跟踪燃料球的位置信息的算法[140]等。
目前，深度学习在球流研究中的应用大多使用的是全连接神经网络，输入为燃料球坐标值、物性参数、球床参数等，输出为流动指标如滞留时间、质量流率等。以图像数据驱动的CNN在球流领域的应用如在岩土、化工领域中对颗粒信息(如尺寸、类型、渗透率等)的提取方面取得了较好的效果。尽管如此，CNN在获取颗粒球信息或预测流动行为方面仍然显示出巨大潜力。
4.1.2 燃料球剩余停留时间和未来分布状态的预测
球流运动的很多信息如燃料球剩余停留时间实际上可通过图像挖掘，该特性对于高温气冷堆来说是一个重要特性，对于了解燃料燃耗、核换料循环和反应堆运行均具有重要意义。不仅如此，在料仓卸料问题中，通过图像快速预测燃料球剩余停留时间将有助于卸料进度的监测和球流状态的评估。除此之外，燃料球的分布变化规律也会反映在图像之中。由于高温气冷堆中的球流是一种稠密缓慢流动，燃料球数目多，循环速率慢，以往这方面的研究往往涉及较大的计算量或者较高的实验成本，而通过CNN方法对燃料球的未来分布状态进行图像预测，能够提高计算效率，更高效地进行参数设计，具有良好的工程应用前景。
鉴于此，本文团队利用CNN对燃料球剩余停留时间与未来分布状态进行了智能预测。
针对球流实验图像，本文团队提出了燃料球剩余停留时间预测的RT-Net[141]模型。RT-Net为直筒式架构，训练时该模型利用多分支卷积组件增强网络的表达能力，提高准确度，而在推理时将多分支卷积组件等价合并为单路径架构，在保留准确度的同时，降低计算量与显存，可以实现更快的推理速度。与VGG、AlexNet对比，RT-Net模型在准确度和计算效率上均能达到最佳的平衡，模型性能获得了提升。
本文认为：为使模型能够预测更广的未知数据，拥有更强的泛化能力，数据集的多样性扩充和应用域的拓宽将是未来有待深入的方向。Pre-Net模型没有考虑图像数据的时间序列特点，将各个图像视为独立对象，而未利用图像之间的序列关系。事实上，CNN与递归神经网络(RNN)组合可能在这些时间序列数据上表现得更好，甚至可以识别一些与时间相关的特征如球的滞留率如何随时间变化等。另外，可以使用主动学习，利用不确定度有针对性地补充数据集以提高模型的外推能力。
4.2 球床传热研究中的神经网络应用
本文团队早先开发了一套神经网络，可用于计算同尺寸燃料球间的角系数[145]。针对不同的物理情形训练了对应的神经网络，通过分情况讨论来提高神经网络的准确性。图 11给出了部分物理情形的示意图。其中，两球情形代表整个工况共有2个燃料球，球间没有阻挡；三球情形代表共有3个燃料球，2个燃料球之间存在1个阻挡燃料球；剩下的情形以此类推。共开发了一球到八球情形对应的神经网络。采用的是全连接神经网络(FCN)。FCN的迭代训练采用LM-BP算法，所有训练在个人主机上完成。所有情形对应的神经网络均包含2层隐含层。神经网络的输入为所有燃料球的相对位置关系，输出为燃料球S1相对燃料球S2的角系数，所有输入输出值都需要归一化处理。解析值和MC方法被用于神经网络的训练和检验。两球情形存在解析表达式，直接采用解析值进行FCN的训练和检验；三球情形以上没有解析表达式，采用MC方法进行FCN的训练和检验。不同情形的训练、检验数据量不同。训练数据最少为1.6万组，最多为20万组，检验数据则从1 600组到1万组不等。简而言之，情形越复杂训练数据量越大，更多细节可见文[145]。
FCN训练完成后，首先结合测试数据检验其预测效果。对于两球情形，采用解析值进行检验；对于更多球情形，采用MC程序生成的数据进行检验。所有测试数据都未参与FCN训练。图 12为其中4种情形的检验结果。数据点越接近y=x直线，表示神经网络的预测效果越好。
5 结论
本文以高温气冷堆堆芯球床流动与传热为重点，回顾总结了围绕球床燃料球高温流动与传热的耦合的数学模型、数值模拟方法及基于神经网络的深度学习技术及其应用。
后续研究主要包括：
1) 对于燃料球动力学，一方面，虽然目前的模型在大型工程问题的计算结果初步验证了其可行性，但是对冷却剂流场的耦合，特别是基于CFD-DEM的大型工程过程中两相流动的模拟方法、软件及验证，仍然需要深入研究。另一方面，对使用连续介质描述的宏观球流本构及控制方程也值得认真研究。
2) 对于燃料球传热学，耦合计算不同尺度的传热，特别是与冷却剂的耦合传热仍然是一个重要的问题。包含传热、传质及燃料球燃耗计算的多物理过程十分复杂，但是工程需求非常迫切，而现有解决方案的效率和准确程度仍然不令人满意。建立一个工业上适用、计算准确、运行高效的堆芯流动与传热的耦合计算程序是下一步研究的重要目标。
