利用EMPS方法对核反应堆严重事故中碎片床形成特性进行的数值研究

《Progress in Nuclear Energy》：Numerical study on debris bed formation characteristics following the FCI of a nuclear reactor severe accident with the EMPS method

　　**摘要**

在燃料-冷却剂相互作用（FCI）过程中形成的碎片床是分析严重核反应堆事故的关键过程。现有的数值研究通常将熔融喷流破碎和碎片床形成过程分离开来，并依赖于理想化的碎片形状假设。本研究通过结合显式运动粒子半隐式（EMPS）方法、等高线连续体表面力（CCSF）模型、碎片识别模型和被动运动固体（PMS）模型，对FCI过程进行了连续的数值模拟。该集成方法追踪了从初始熔融注入和固化到最终碎片床积聚的动态演变过程。所建立的数值模型通过实验数据进行了验证，其宏观预测与物理测量结果一致。参数分析揭示了控制碎片床形态的机制。提高熔融铅铋共晶（LBE）温度、喷流速度和水温会增加碎片床的扩散面积。较高的LBE和水温会导致边缘壳层重新熔化并产生二次膨胀。此外，随着这三个参数的增加，碎片床的积聚高度会降低。总的孔隙率随着LBE温度和喷流速度的升高而增加，在冷却水温度较高时则会降低。所提出的建模方法为碎片床的形成提供了定量见解，并可能有助于制定严重事故管理策略。

**引言**

在核反应堆的安全分析中，特别是针对严重事故，高温核熔体与冷却剂的相互作用（FCI）是一个关键且复杂的多相现象（Berthoud, 2000; Sehgal, 2012）。如图1所示，严重的燃料棒故障会导致熔融的核燃料向下移动到剩余的水池中。在初始的热液破碎后，固化的熔融碎片不断沉降并在下层空间积聚，最终形成多孔的碎片床（Karbojian et al., 2009; Magallon, 2006）。这种碎片床的形态演变及其冷却能力具有重要意义，因为重新熔化的碎片可能会引发严重的局部热点，威胁反应堆压力容器的完整性（Yakush et al., 2015）。因此，准确解析熔融破碎的复杂动态及其后的碎片床形成过程对于缓解严重事故至关重要。

传统的严重事故系统分析代码主要依赖于半的经验性本构关系，这些关系难以解决微观尺度上的拓扑变形和非线性破碎动态。为了实现多相界面的微观追踪，基于欧拉网格的界面捕获技术，如体积流体（VOF）方法（Gopala and van Wachem, 2008; Hirt and Nichols, 1981; Rider and Kothe, 1998）、水平集方法（Osher and Sethian, 1988; Sethian and Smereka, 2003; Sussman et al., 1994）以及水平集-体积流体（CLSVOF）方法（Sussman and Puckett, 2000）已被广泛应用于FCI的预混合阶段。然而，正如Meignen等人（2014）所指出的，在喷流破碎过程中，界面剧烈撕裂、微小液滴的脱落以及拓扑结构的极端不连续性非常普遍。面对这种极端的流体变形，传统的基于网格的方法不仅容易产生严重的网格畸变和计算发散，而且平流项的离散化方案不可避免地引入了数值扩散，这导致相位界面被人为地过度平滑，从而严重影响了碎片破碎尺度的预测准确性和可靠性。

为了克服欧拉框架在极端大变形模拟中固有的拓扑重建瓶颈，纯拉格朗日无网格粒子方法展现了天然的算法优势。运动粒子半隐式（MPS）方法由Koshizuka等人（Koshizuka and Oka, 1996）首创，它将连续流体离散为具有多种物理特性的拉格朗日粒子。通过使用核函数来近似微分算子，MPS方法内在地绕过了由非线性平流项引起的数值耗散。在探讨流体动力破碎的基础研究中，Ikeda等人（2001）成功地利用MPS方法再现了冷却剂注入条件下的喷流穿透和界面不稳定性。随后，Nomura等人（2001）引入了表面张力模型来模拟液-液和气-液系统中的液滴破碎，量化了不同破碎状态下的韦伯数，并揭示了蒸汽膜对液滴破碎的显著抑制作用。Shibata等人（2004）进一步阐明了连续喷流破碎长度对韦伯数和弗劳德数的依赖性，验证了通过粒子方法获得的碎片尺寸分布与经典的Nukiyama-Tanasawa经验相关性高度一致。对于更复杂的FCI场景，Li等人（2019, 2020）基于改进的MPS框架进行了系统的数值研究。这些研究准确地揭示了惯性主导的破碎和表面张力主导的凝聚的两阶段动态机制，并确认了初始蒸汽膜的快速脱落对液滴变形的影响可以忽略不计。这些发现全面验证了MPS方法在解决FCI过程中的大规模界面变形和复杂流体动力学方面的计算能力。

在熔融喷流和液滴破碎之后，材料最终固化并沉降形成多孔的碎片床。为了解决离散颗粒之间的接触动力学问题，开发了离散元方法（DEM）（Cundall and Strack, 1979）和被动运动固体（PMS）模型（Koshizuka et al., 1998）。在上述框架的基础上，PMS模型、DEM方法及其相关的耦合技术构成了模拟碎片床积聚的基础数值方法。在流体-固体耦合模拟方面，Duan等人（2022）提出了一个多相MPS框架，将流体-固体相互作用与固体-固体碰撞结合起来，强调了大规模碎片的重新熔化可能会在较低位置引发严重的局部热点，可能导致压力容器失效。对于纯固体颗粒积聚的形态演变，许多学者基于Karbojian等人（2009）的实验获得的碎片尺寸分布进行了数值验证。Muramoto等人（2022）利用MPS-DEM耦合框架系统地验证了通过调整摩擦系数来等效模拟不规则碎片形状的可行性。为了进一步探索宏观碎片床的动态形成机制，Ding等人（2022）开发了一个MPS-DEM耦合模型，以处理复杂的流体-固体耦合和固体-固体碰撞。通过粒子坠落实验的模拟，这种耦合方法准确预测了非热传递条件下的碎片颗粒的动态沉降和积聚特性。

尽管MPS方法在捕捉多相界面方面具有固有的优势，但在实际规模的反应堆严重事故中，颗粒数量的指数级增加导致了计算成本的显著上升，这主要是由于压力泊松方程（PPE）的隐式求解。为了解决这一问题，基于状态方程（EOS）的显式求解框架被开发出来。Shakibaeinia等人（Shakibaeinia and Jin, 2010）提出了弱压缩MPS（WC-MPS）方法，通过引入人工声速将不可压缩流体转换为弱压缩流体。这种方法成功绕过了大规模的矩阵求解过程，显著提高了自由表面流动的计算效率。为了进一步抑制显式计算中的高频压力振荡，Jandaghian等人（2021）开发了一种改进的WC-MPS方法，结合了动态粒子碰撞和优化的粒子移动技术，显著提高了剧烈多相流动和大变形过程中的数值稳定性。随后，Oochi等人（2010）和Yamada等人（2011）改进了显式MPS（EMPS）方法的基本架构，并成功地将其与大规模GPU并行加速技术深度融合。基于这一基础，Zheng等人（2020）提出了一种三维耦合的EMPS-FEM（有限元方法）框架，利用虚拟网格边界模型优雅地克服了处理复杂壁面的难题，实现了大规模流体变形和弹性结构动力学的有效模拟。

现有的数值研究通常将FCI过程划分为独立的阶段，忽略了破碎过程中的动态固化，并且在积聚时受到预定义的规则形状颗粒的限制。为了克服这些关键限制，本研究采用了一种基于显式并行架构的流体-固体耦合框架（EMPS-PMS）。这种集成方法能够实现包括熔融注入、动态相变和最终床形成在内的连续模拟。通过避免理想化的碎片假设，所提出的模型能够提供关于严重事故条件下碎片床的形态演变、总孔隙率分布和整体积聚特性的定量预测。

**控制方程**

在EMPS方法中，流体域被离散为一系列拉格朗日粒子。这些粒子携带质量、速度和温度等物理属性。在本研究中，熔融的铅铋共晶（LBE）和冷却水都被视为弱压缩的粘性流体。基于拉格朗日框架，质量、动量和能量守恒方程表示如下：
$$
\mathbf{D}\mathbf{\rho}\mathbf{v} = 0, \quad \mathbf{D}\mathbf{\rho}\mathbf{u} = -\nabla\mathbf{P}, \quad \mathbf{D}\mathbf{\rho}\mathbf{T} = -\nabla\mathbf{Q}
$$
其中ρ是流体密度，u是速度。

**模拟域和边界条件**

为了再现实验中观察到的熔融固化与积聚过程，构建了一个1:1比例的粒子几何模型。如图4所示，计算域包括一个上引导释放装置和一个长度和宽度均为40厘米的下水槽，冷却水深度为10厘米。此外，初始粒子间距设置为0.005米。

**结论**

本研究利用EMPS数值方法研究了熔融LBE在水池中的沉积和碎片床形成过程。通过在不同条件下的对比模拟，阐明了LBE温度、喷流速度和冷却水温度对碎片床宏观形态和微观特性的影响机制。主要结论总结如下：
(1) 与实验数据的比较验证了所建立的模型。

**作者贡献声明**

Jinhua Que：撰写—原始草稿、软件开发、数据整理。
Yunlong Liao：撰写—审阅与编辑、可视化、验证、形式分析。
Gen Li：撰写—审阅与编辑、监督、方法论制定、资金获取、概念构思。

**利益冲突声明**

作者声明他们没有已知的可能影响本文工作的财务利益或个人关系。

**致谢**

本研究得到了中国广东省基础与应用基础研究基金（项目编号2024A1515010335）的支持。

Jinhua Que | Yunlong Liao | Gen Li
华南理工大学电力工程学院，广州，510641，中国