基于核独立快速多极子方法（KIFMM）增强的离散单元法（DEM）复杂多物理场系统建模

《Advanced Powder Technology》：Modeling complex multiphysics systems with discrete element method enriched with the kernel-independent fast multipole method

　　

该论文描述了离散单元法（Discrete Element Method, DEM）代码MercuryDPM与核独立快速多极子方法（Kernel-Independent Fast Multipole Method, KIFMM）实现的耦合。该组合模拟框架能够处理大量的多尺度问题，包括在细观尺度上的颗粒机械相互作用以及在粗观尺度上各种性质的长程相互作用。这些问题包括粉末、黏土和颗粒中的静电相互作用，铁磁颗粒中的磁相互作用，以及小行星云中的引力相互作用。刚性团块（rigid clumps）的形式主义已成功与KIFMM结合，使得能够处理涉及具有任意电荷分布的非球形颗粒之间的复杂长程相互作用问题。研究人员的方法能力通过几个应用实例得到展示，包括立方盒中相互作用的排斥电荷系统、复杂形状偶极子颗粒的聚集以及空间尘埃引力坍缩成星状体。

论文解读：基于核独立快速多极子方法（KIFMM）增强的离散单元法（DEM）复杂多物理场系统建模

研究背景与意义：

离散单元法（Discrete Element Method, DEM）自P. Cundall提出以来，已成为模拟颗粒固体力学性能的重要工具，广泛应用于宽粒径分布、复杂非球形形状、可变形颗粒、复杂接触力学以及范德华（vdW）相互作用结合的软物质预测建模等领域。现有的DEM项目（如MercuryDPM、LIGGGHTS、YADE等）展现了其架构的成熟性、高效性、跨兼容性及在不同高性能计算架构上的高效并行化能力。然而，DEM的一个严重局限性在于其无法有效处理涉及长程相互作用的问题。现有DEM代码的核心是用于检测颗粒间相互作用的流程，无论是面向密集颗粒组件的扫掠和剪枝算法（sweep-and-prune algorithm），还是面向颗粒流动模拟的链接单元算法（linked cell algorithm），均旨在高效检测短程（机械）相互作用。分析长程相互作用本质上是不同的问题，需要不同的相互作用近似、区域划分和并行化等方法，目前尚缺乏能高效结合DEM与长程相互作用的代码。当颗粒通过某种性质的长程力相互作用且不能忽略时，若需考虑所有可能的成对相互作用，会导致力计算复杂度为O(N²)。尽管通常只需以特定精度评估力，此类问题存在渐近更快的方法（如Barnes-Hut树算法、快速多极子方法（FMM）及后续的H-、H²-和HSS-矩阵等分层矩阵框架），可将复杂度降至O(N log N)甚至O(N)。此外，许多DEM应用需要模拟由颗粒产生或作用于颗粒的连续场（如黏性流场、周围介质中的弹性位移、极化产生的静电场），这类问题中颗粒相互作用由底层偏微分方程（PDEs）介导，可重新表述为颗粒表面的边界积分方程（BIEs），虽避免了体积离散，但会导致复杂度为O(N²)的稠密线性系统。快速求和法（如用于长程颗粒相互作用的方法）也可用于加速BIE公式，将计算成本降低至近线性复杂度。此类快速求和框架的关键要求包括：支持广泛核且无需核特定的解析展开、可扩展的并行实现、低场评估成本（或在BIE公式情况下的矩阵-向量乘积）。这些考量促使研究人员在本工作中使用核独立快速多极子方法（KIFMM）。广泛使用的MD模拟框架（如LAMMPS）虽支持长程相互作用（通常通过Ewald-type方法处理静电），但依赖周期性边界条件且对均匀系统最高效；相比之下，此处采用的快速多极子框架无需周期性，更适用于非均匀颗粒分布、复杂几何形状和更广泛的相互作用核。本研究聚焦于加速具有长程相互作用的DEM模拟，特别关注使用KIFMM模拟长程力时，由刚性团块（rigid clumps）建模的非球形颗粒的模拟。

主要关键技术方法：

研究人员开展了MercuryDPM（离散单元法代码）与PVFMM（KIFMM的公共可用库）的耦合工作。两者均使用CMake构建系统，集成较为直接。抽象颗粒类增加了“电荷”（单实值参数）属性，力计算流程扩展了累积长程相互作用力的阶段。对于每个驱动问题类，从标准Mercury3D类派生，在setupInitialConditions()方法中初始化PVFMM组件（包括构建初始树结构和预计算所选核的平移算子），在actionsAfterTimeStep()方法中执行长程力（及对应势能贡献）评估。对于团块颗粒，电荷分布在卵石（pebble）成分上，力首先施加在卵石层面，随后聚合到每个团块的总力和扭矩。由于机械碰撞的最大刚度k_DEM与长程相互作用刚度k_LR通常差几个数量级，研究人员设置了长程力计算步长（stride）S≈√(k_DEM/k_LR)（默认S=10）。PVFMM计算（时间步中最耗时的部分）以MPI并行模式执行，MercuryDPM驱动通过OpenMP采用共享内存并行。该框架支持周期性边界条件（需满足特定正则性/中性条件，如Laplace/Coulomb核要求单周期单元总电荷为零）。研究人员通过三个关键技术方案实现耦合：1) MercuryDPM粒子动力学代码的使用（支持球形、超二次曲面颗粒及刚性团块，采用Velocity Verlet算法更新位置，前向欧拉方案更新旋转，提供广泛短程接触模型）；2) 核独立快速多极子方法（KIFMM，通过PVFMM库实现，基于八叉树划分域，通过等效密度和辅助包络面表示集群场，避免核特定的解析展开，支持标量、向量或张量值核，本文使用Laplace核及其梯度评估标势和相关力场）；3) 刚性团块（rigid clumps）的应用（由N个给定密度的刚性球形颗粒（卵石，pebble）刚性连接而成，具有集体质量和惯性张量（TOI），作为刚体处理6自由度，可通过卵石分配电荷实现复杂电荷分布）。

研究结果：

1. 1.
   引言（Introduction）
   研究人员指出DEM的发展与应用扩展，以及其在长程相互作用问题上的局限，引出快速多极子方法等O(N)或O(N log N)复杂度方法，以及BIE公式中快速求和法的应用，阐述使用KIFMM的动机，对比LAMMPS等MD框架的局限性，明确本研究聚焦于DEM与长程相互作用的加速及刚性团块结合KIFMM的模拟。
2. 2.
   方法（Methods）
   2.1 MercuryDPM粒子动力学代码：介绍MercuryDPM的开源属性、开发背景（2009年特温特大学起始）、模块化面向对象C++架构、关键特征（灵活几何与边界条件、集成分析工具包、高级接触检测算法适用于强多分散颗粒）、默认颗粒表示（球形，支持超二次曲面颗粒和刚性团块），以及典型模拟流程（Velocity Verlet更新位置，前向欧拉更新旋转，丰富短程接触模型库）。
   2.2 核独立快速多极子方法（KIFMM）：介绍KIFMM用于加速DEM中长程相互作用评估，采用PVFMM库（并行实现Ying等人提出的方法）。说明直接成对求和的O(N²)成本问题，FMM核心思想（近场直接准确评估，远场集体效应压缩表示，分层组织复用降低至O(N)），八叉树划分域，流程包括后序遍历（叶到根）的S2M（源到多极）和M2M（多极到多极）平移，交互阶段M2L（多极到局部）平移，前序遍历（根到叶）的L2L（局部到局部）传播，近场直接评估。KIFMM特点：避免核特定解析展开，通过等效密度和辅助包络面表示集群场，转移操作为等效密度映射矩阵，支持多核类型，本文使用Laplace核及其梯度。
   2.3 刚性团块（rigid clumps）：介绍刚性团块定义（N个刚性球形卵石刚性连接，6自由度，集体质量和惯性张量（TOI），接触检测视为集合，时间积分视为单个颗粒），以及通过卵石分配电荷实现复杂电荷分布的能力。
   2.4 MercuryDPM驱动与PVFMM库的耦合：说明CMake集成方式，抽象颗粒类增加电荷属性，力流程扩展长程力阶段。问题类派生自Mercury3D，setupInitialConditions()初始化PVFMM（树结构、平移算子预计算），actionsAfterTimeStep()执行长程力评估（团块电荷分布到卵石，力聚合到团块）。长程力计算步长S的设置依据刚度差异，默认S=10。周期性边界条件的应用及中性条件要求。混合并行架构（PVFMM用MPI，MercuryDPM用OpenMP），支持大规模模拟（达10⁶颗粒和10⁶时间步）。未来扩展方向（BIE求解器、全并发MPI并行耦合）。
3. 3.
   结果与讨论（Results and discussion）
   研究人员通过四个示例展示MercuryDPM-PVFMM模拟能力：
   3.1 盒中相同带电颗粒：模拟10⁴个相同均匀带电（半径0.0025）颗粒，排斥静电力（F_ij=r(1/(4π))(q_iq_j/r³)）加短程线性粘弹性摩擦接触模型，单位立方盒弹性壁。初始均匀分布，短暂瞬态后所有颗粒迁移到盒壁（边缘和角落浓度更高），势能演化显示快速沉积到壁后沿壁缓慢重排，最终平衡配置非结晶（强边缘效应）。模拟耗时2小时（CPU节点），参数见Drivers/LongRange/Charges/。
   3.2 相互作用偶极子的聚集：模拟2×10⁴个刚性团块（每个3个半径0.005的卵石，间距0.01，电荷-1、0、1），短程接触模型同前，单位盒弹性壁，外加均匀引力场和背景拖拽（模拟周围流体介质）。偶极子聚集成多个近刚性树枝状聚集体，缓慢沉积到盒底；能量演化显示粘附能单调降低至平台（稳定聚集体网络），初始弹性峰值后消失，动能快速淬灭，引力能随聚集体沉降降低，偶极-偶极相互作用能U_LR单调降至平衡水平。模拟耗时14小时，参数见Drivers/LongRange/DipolesAggregation/。
   3.3 偶极血小板堆叠：模拟50个血小板（每个1094个卵石，六边形密排双层，层分别为正负电荷，N=14，卵石间距0.01，半径0.008），全轴周期性边界条件，初始共取向无重叠单位盒，小时间步（确保倾斜几何团块时间积分稳定），长程计算步长S=20。结果显示偶极子自组装成近理想堆叠（穿过周期边界），平衡堆叠长度略小于周期边界距离（形成间隙），经历高度动态大尺度重排；能量演化显示长程势能U_LR、引力势能、弹性势能、动能（含较大比例旋转动能）的变化。模拟耗时20小时，参数见Drivers/LongRange/Clay/。
   3.4 弹性颗粒的引力坍缩：模拟10⁵个相同球形颗粒，长程引力（F_ij=r_ijG*m_im_j/r_ij³）加短程耗散弹性摩擦碰撞，初始位置均匀随机在(0.2,0.8)³，初速度零，立方(0,1)³弹性壁（无法用周期性边界条件，否则引力势发散）。数周期快速衰减振荡后，颗粒坍缩成近球形聚集体（“星状体”），平衡由弹性（电磁）力和引力平衡；零初始角动量（冷启动）且成对相互作用不产生角动量，星状体无净旋转和离心变形。该示例可作为天体力学模拟起点，10⁵颗粒、10⁵时间步模拟耗时12小时（普通桌面机），10⁶颗粒仍在单桌面机能力范围内（性能与无PVFMM的DEM代码相当）。参数见Drivers/LongRange/Galaxy/。
4. 4.
   结论（Conclusion）
   研究人员总结：该工作展示了结合处理短程机械相互作用的DEM（MercuryDPM）与处理长程相互作用的KIFMM的建模框架能力。通过将带电球形颗粒集合表示为刚性团块，能够准确处理复杂的、传递力矩的长程力，同时保持简单的相互作用核。当前实现成功处理天体力学、静电相互作用和偶极-偶极耦合问题；带电尘埃动力学可结合先进充电模型解决复杂多物理场现象（如光刻和等离子体加工中表面尘埃扬起、无空气行星体尘埃动力学），偶极聚集示例与黏土类材料及其相变建模相关，引力坍缩示例与尘埃聚集体天体力学（如土星环）相关。框架可进一步扩展：关键方向是KIFMM加速的边界积分方程（BIE）求解器，纳入BIE方法可实现外部场中粉末颗粒介电极化的自洽建模（此类模拟的核心问题）、磁静力学类比公式，以及刚性颗粒与可变形连续体（复合材料弹性BIE、颗粒流Stokes BIE）的耦合相互作用，从而构建广泛科学工程应用的通用多物理场模拟环境；此外，框架在提升并行性能和可扩展性方面有大量机会。所有代码见指定仓库，示例位于/Drivers/LongRange/。

（注：以上内容严格依据论文原文浓缩总结，未添加额外推测信息，专业术语均备注英文缩写，保留原文上下标格式，去除引用和图示标识，作者名保留原文。）