多相流涉及两种或更多不相溶的相，在自然环境和工业应用中都很常见[1]。它们在各种应用中发挥着关键作用，包括石油和天然气开采[2]、海洋和环境建模[[3], [4], [5]]、化学工程[6]等。然而，由于复杂的界面动力学、强烈的相耦合以及涉及的空间和时间尺度范围广泛，多相流的准确数值模拟仍然极具挑战性。
在过去的几十年中，基于网格的方法被广泛用于多相流模拟，包括体积流体（VOF）方法[7]、水平集（LS）方法[8,9]、前沿追踪方法（FTM）[10]和格子玻尔兹曼方法（LBM）[11]。它们的鲁棒性和数值精度使它们成为解决各种多相流问题的有效工具。然而，它们在捕捉复杂拓扑变化（如聚合、破裂和碎裂）方面常常面临限制[7]。
与传统的基于网格的方法不同，平滑粒子流体动力学（SPH）方法是一种完全拉格朗日和无网格的方法。SPH方法将流体建模为离散的移动粒子，使其非常适合处理涉及大变形、自由表面以及界面破裂、飞溅和聚合的问题[12]。最近，SPH方法在各种流体动力学问题中得到了广泛应用，如自由表面流动[13,14]、多相流[5,15]和流体-结构相互作用[16,17]。
尽管SPH方法具有优势，但在多相流模拟过程中仍面临一些限制。特别是，准确解析界面动力学通常需要在相界面附近显著提高空间分辨率[18]。然而，大多数现有的SPH实现在整个计算域内采用均匀的空间分辨率，这导致在高分辨率应用时计算开销过大，无论是时间还是内存使用。为了解决这个问题，提出了多分辨率SPH方法，以在需要的地方提供局部细化的分辨率，而在其他地方保持较粗的分辨率。Feldman和Bonet[19]在流体流动问题中提出了一种局部细化技术，其中一个粒子被分裂成四个小粒子。Vacondio等人[20,21]通过引入粒子聚合过程开发了一种完全多分辨率SPH方法。在他们的方案中，一个粒子在二维中被细化为七个更小的粒子，在三维中被细化为十四个粒子，通过合并两个粒子来实现粗化。Barcarolo等人[22]采用了一种细化策略，其中二维中的每个母粒子被分成四个粒子，而母粒子本身不被删除，而是保持在不活跃状态。这些不活跃粒子的运动通过来自邻近粒子的插值来控制。基于这种方法，Chiron等人[23]开发了一种自适应粒子细化（APR）方法，提高了鲁棒性并放宽了约束。上述方法通常需要预先识别需要细化的区域，这适用于感兴趣区域固定的问题。然而，在许多多相流场景中，需要高分辨率的区域（如界面、剪切层或涡旋结构）会随时间动态变化。最近，引入了几种完全自适应方法来解决这一限制。Yang和Kong[24]提出了一种用于多相流的自适应分辨率SPH方案，其中细化是根据与界面的距离来确定的应用到相应区域的。Yang等人[18]进一步将这种方法扩展到自由表面流动问题。在之前的研究[18,24]基础上，Hong等人[25]将这种方法应用于工业规模行星混合器中的非牛顿浆料混合。Muta和Ramachandran[26]引入了背景粒子来指定二维SPH模拟中的空间分辨率。
虽然多分辨率或自适应分辨率SPH方法可以显著降低计算成本，但大尺度多相流模拟仍然非常具有挑战性。三维多相问题可能涉及数千万甚至数亿个粒子，即使使用多分辨率或自适应分辨率方法，对于串行程序来说仍然是难以处理的。鉴于高性能计算资源的快速发展，采用并行计算框架为实现大规模三维多相SPH模拟提供了有效的解决方案，特别是与自适应分辨率方法结合使用时。许多研究人员开发了并行SPH求解器以提高计算效率和可扩展性。现有框架采用了图形处理单元（GPU）加速[15,27,28]、带有消息传递接口（MPI）[[[29], [30], [31]]的分布式内存并行化，以及共享内存方法，如OpenMP[[32], [33], [34]]和Intel Threading Building Blocks（TBB）[35]。此外，SPH模拟中还采用了混合并行化策略，如MPI–OpenMP[36]和MPI–TBB[37]。除了并行化策略外，由于基于粒子的计算引入的同步开销，有效的负载平衡对于大规模并行SPH模拟也至关重要。Oger等人[31]通过正交递归二分（ORB）方法实现了负载平衡，使得在32,768个核心上可以模拟多达十亿个粒子的SPH。Guo等人[38]使用希尔伯特空间填充曲线实现了不可压缩SPH的负载平衡，同时利用空间局部性来提高缓存命中率。Zhang等人[39]通过使用METIS（一种基于图的划分工具）实现了负载平衡。这些方法通常通过分配给每个进程的粒子数量来量化计算工作负载。Zhu等人[30]提出了一种基于MPI的自由表面流动并行SPH框架，其中使用粒子数量和计算时间作为负载指标来实现动态负载平衡。然而，在并行框架中实现多分辨率或自适应分辨率存在一些挑战。首先，频繁的粒子分裂和合并需要复杂的进程间通信，可能会引入竞争条件。其次，如果处理不当，不同的相互作用范围可能会导致过多的通信。此外，多分辨率或自适应分辨率引起的空间非均匀性会加剧负载不平衡，降低并行效率。为了解决这些挑战，一些研究提出了支持多分辨率的并行SPH框架。在[30]的工作基础上，Zhu等人[40]将多分辨率纳入并行框架，为不同进程分配不同的分辨率，并采用了一种局部时间步进方法，其中每个进程的时间步长根据其分辨率来确定。Ji等人[37]结合了质心Voronoi粒子（CVP）域划分策略和局部嵌套层次数据结构，以实现具有自适应平滑长度的高效并行SPH模拟。然而，很少有研究关注自适应分辨率与大规模并行SPH框架的集成及其性能优化。
在这项工作中，我们提出了一种用于大规模多相流模拟的大规模MPI并行自适应分辨率SPH方法。为了实现高性能和可扩展性的大规模并行执行，我们提出了一种并行自适应分辨率方法，该方法动态识别细化区域并通过分裂和合并调整粒子分辨率，允许每个子域内存在不同的分辨率。为了防止数据竞争，我们采用了一种对称的成对合并算法，确保进程间的操作无冲突。此外，我们设计了一种基于历史的自适应幽灵粒子管理策略，根据实时局部相互作用需求调整交换区域，从而最小化不必要的通信。此外，为了解决自适应细化引入的严重粒子负载不平衡问题，我们引入了一种逐轴分层负载平衡方法，根据从背景网格采样的相互作用粒子对的数量更新子域边界。
本文的组织结构如下。第2节介绍了多相SPH模型。第3节描述了并行自适应分辨率方法。第4节介绍了高性能并行化策略。第5节提供了几个数值示例，以验证所提出方法的准确性、鲁棒性和可扩展性。第6节以一些结论性评论结束。
