一种用于三维离散元（3D DEM）模拟的高效密集球体堆积方法，适用于复杂域环境

《Powder Technology》：An efficient dense spheres packing method in complex domains for 3D DEM simulations

　　李乃平|刘和|夏明|曲同明|冯宇涛
中南大学肝脏疾病研究中心，中国长沙410083

**摘要**
高效生成具有理想粒径分布的逼真3D粒子组装体是一个重大挑战，这主要是由于复杂的内部和外部几何约束。本研究提出了一种将球体密集填充到任意3D形状中的方法。通过明确考虑边界域，该方法显著减少了最终组装体表面的间隙和空洞。所提出的方法可以通过三步程序处理从STL文件导入的任意复杂几何形状：（1）生成初始粒子填充；（2）压缩初始填充；（3）填充剩余空间。通过几个示例展示了该方法的高效率、低硬件要求以及对复杂几何形状的强适应性。值得注意的是，该算法仅使用2GHz笔记本电脑在1218.2秒内就在圆柱形域中生成了1,084,791个球体，填充密度达到58.6%。

**引言**
离散元方法（DEM）是模拟不连续材料行为的成熟工具[3]、[7]、[8]、[9]、[16]、[28]、[29]。使用DEM首先需要通过粒子填充来生成几何模型。由于复杂的内部和外部域约束，高效创建具有特定粒径分布的大规模逼真3D组装体是一个重大挑战。
作为最简单的3D形状，球体是最广泛使用的离散元素。在过去的二十年里，已经开发了三种主要的球体填充算法[33]：动态方法[15]、[27]、构造方法[5]、[22]、[24]以及耦合的构造-动态方法[30]。虽然动态方法可以使用商业或开源软件在复杂域中生成3D组装体，但对于涉及数百万个球体的系统来说效率较低，因为需要大量的数值步骤，生成逼真组装体会非常耗时。
相比之下，构造方法效率很高，因为它仅依赖于几何条件而不需要进行动态模拟。这一类别包括两种主要方法：基于网格的方法和前进前沿方法（AFM）。基于网格的方法首先将3D域离散为体积网格，然后将球体放置在网格顶点或四面体中心[2]、[14]。这种方法易于实现，并且可以使用现有的FEM网格生成软件处理任意几何形状。然而，得到的粒子半径分布可能不符合目标，而且由于球体不一定接触，可能会出现较大的空洞。为了解决这些问题，冯等人[6]引入了用于填充圆盘的AFM。此后，该方法已被扩展到复杂域中的椭圆和凸多边形[4]、[33]、[36]。然而，将AFM适应到3D环境中并不简单，因为构建和更新前进前沿很复杂。为了避免这种复杂性，罗和王[20]提出了一种用于无界3D域的球体填充算法，无需显式构建前进前沿，尽管没有考虑边界条件。刘等人[19]将种子扩展方法与AFM结合，使用局部Delaunay镶嵌和距离函数填充有界域。他们的方法重新计算边界附近球体的位置和半径以确保精确接触。李和纪[18]通过空间网格提高了效率以简化前沿检测。同样，托等人[32]引入了一种基于解析求解三角测量方程的顺序填充算法。然而，复杂几何形状中的填充问题尚未完全解决。后来，洛萨诺等人[21]提出了一种高效的算法，用于将具有预定尺寸分布的数百万个球体填充到任意几何形状中。它使用有符号距离场来确定球体是否位于域内。但是，与边界重叠的球体会被移除，这可能会在表面附近留下空洞。除了这两种构造方法外，韩等人[10]还提出了一种基于几何的压缩算法，以其简单性和处理简单形状的高效率而著称。
对于简单几何形状，构造方法和动态方法都表现良好。然而，对于复杂的3D几何形状，构造方法比动态方法更具挑战性。据作者所知，目前仍缺乏一种在复杂3D几何形状中具有预定粒径分布的高效简单球体填充算法。本文通过扩展韩等人[10]最初提出的基于几何的压缩算法来解决这一空白。立体光刻（STL）文件格式是一种表示复杂几何形状的多功能标准，广泛应用于CAD、3D打印和计算机图形学。本研究采用STL文件格式来定义这些几何形状。与原始方法相比，所提出的方法具有三个关键优势：（1）使用改进的方案大大提高了简单形状的填充效率；（2）明确的边界考虑最小化了最终组装体表面的间隙和空洞；（3）该方法可以直接处理从STL文件导入的任意复杂几何形状。

**本文的其余部分组织如下：**
第2节描述了从STL文件导入的复杂几何形状的初始填充生成过程。第3节介绍了压缩和摇动初始配置以增加复杂几何形状中的填充密度。第4节提出了一种新的填充方案，用于在压缩后填充剩余的空洞。然后，第5节通过涉及简单和复杂几何形状的几个数值示例展示了所提出方法的性能和有效性。第6节讨论了所提出方法的各个方面。最后，第7节总结了本文的主要发现。

**示例和性能**
通过在不同域中的球体填充示例展示了所提出方法的多功能性和能力。该算法在配备2.0GHz AMD R5-7530U CPU的笔记本电脑上使用自制的Fortran95代码实现。

**势能容忍度的敏感性分析**
在算法2和算法3中，“势能”的变化作为收敛标准，其中球体的“高度”hc和“势能”pc定义为：
$$
hc = x_c - x_0 \cdot v_g, \quad pc = h_c \cdot r_c^3
$$
其中$x_c$和$r_c$分别是球体的中心坐标和半径；$x_0$是参考点的坐标，在这里取为边界长方体的最低点（即图1中的[xmin, ymin, zmin]）。
系统的总“势能”Pc则为：
$$
P_c = \sum_p c_p
$$

**结论**
本研究提出了一种基于几何的压缩算法，用于在复杂3D几何形状中密集填充球体。所提出的方法通过三步程序处理从STL文件导入的任意对象：生成初始填充、压缩初始填充和填充剩余空间。在初始填充生成过程中，引入了一种新方案来确定新放置的球体是否位于复杂几何形状的封闭表面内。与其他方法相比，...

**作者贡献声明**
李乃平：撰写——原始草稿、资源、方法论、数据管理。
刘和：可视化、软件、形式分析。
夏明：撰写——审阅与编辑、项目管理、资金获取、概念化。
曲同明：撰写——审阅与编辑、验证、方法论、研究。
冯宇涛：撰写——审阅与编辑、监督、资金获取。

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

**致谢**
本研究得到了国家自然科学基金（编号52574167和12072217）、湖南省自然科学基金（编号2025JJ50215）和湖南省教育厅科学研究基金（编号25A0141）的财政支持。作者感谢匿名审稿人的建设性评论和建议。