一种基于3D NURBS的方法，用于对颗粒介质的DEM模拟中的任意形状颗粒进行直接数学描述

《Powder Technology》：A 3D NURBS-based approach to direct mathematical description of arbitrary-shaped particles for DEM simulation of granular media

　　亚历山德拉·L·热列兹尼亚科娃，俄罗斯科学院力学问题研究所，俄罗斯莫斯科维尔纳茨基街101/1号，邮编119526

摘要
在离散元方法中，准确描述现实颗粒复杂几何形状的数学模型对于可靠模拟颗粒材料的物理行为至关重要。本研究提出了一种利用非均匀有理B样条函数高精度重建任意形状颗粒表面的通用方法。基于NURBS的自由形态颗粒描述已完全集成到三维离散元方法中。同时，还开发了用于确定颗粒间/与壁面接触、相互穿透深度以及相应相互作用力的模型，这些都是离散元方法技术所依赖的。一种新方法通过利用NURBS曲面的控制点及构成表面块来构建递归边界体积层次结构，从而确定交叉状态并逐步细化一般非凸NURBS颗粒的接触区域。通过一些实例验证了该方法在数值重建任意形状容器中真实颗粒堆积时的计算效率与可靠性。实验结果表明，所模拟的堆积在堆积密度、孔隙比和安息角等结构特性方面与真实情况在定性和定量上均相符。此外，该方法还能为各种固体材料（如金属、岩石和干面制品）的动态及准静态行为提供足够精确的数值预测。

引言
自Cundall和Strack首次开展相关研究以来，离散元方法已成为模拟颗粒材料在不同载荷和环境条件下的物理响应的最常用方法。这一计算技术得益于创新数值方法的开发以及处理能力的飞速提升，已成功用于模拟多种类型颗粒的机械行为。传统的离散元模拟基于对每个非球形颗粒的平动和转动运动的追踪，同时整合了包含碰撞力的微观运动方程。建立合适的碰撞模型以计算颗粒间/与界面间的接触力是该方法面临的关键挑战。首先，模型需具备足够精确度、计算效率及可靠性以实现接触检测；其次，需要构建合理的碰撞力模型关系，这些关系通常涉及重叠体积和相对速度；最后，还需根据单个颗粒的机械特性（如质量、接触刚度、摩擦系数等）来确定消除交叉所需的单颗粒响应程度。

对真实颗粒几何形状的数学描述同样极为重要。颗粒结构越复杂，计算其交叉与相互作用所需的资源就越多。此外，方法的准确性、通用性及效率取决于任意形状颗粒表面的几何表示方式。因此，开发一种可高精度重建复杂颗粒的通用方法，并能直接应用于离散元等方法中，已成为当前的重要研究方向。

在粉末技术领域的诸多应用中，球形颗粒的近似模型被广泛且成功地使用[1,[3], [4], [5], [6]]。由于能够快速检测和分析接触，这种方法大大降低了计算成本。在这种情况下，滚动阻力被用作模拟颗粒形状的一种方法[3,7]。然而，要准确描述单个颗粒的机械响应，则必须考虑实际复杂的颗粒形状。

本研究关注的是在外部力作用下、处于任意形状边界内的复杂形状颗粒的堆积形态及其机械特性。文献综述表明，在离散元模拟中，已有许多精确且高效的方法可用于描述通常通过粘弹性摩擦接触（其微观接触特性包括泊松比、杨氏模量、恢复系数、摩擦系数等）相互作用的非球形颗粒。这些数值方法在描述复杂不规则颗粒的数学方式上各有不同。

在三维离散元建模中常用的形状抽象形式包括：团簇球形（多球形）[8,9]、多面体型[[10], [11], [12], [13], [14], [15]]、基于三角形网格的模型[16,17]以及基于体素的模型[[18], [19], [20], [21], [22]]。此外，还有一些较少使用且鲜为人知的方法，例如[23]。

尽管上述方法已在离散元技术中成功解决了众多问题，但仍有特定缺陷尚未解决。全面评估揭示了这些不足，主要体现在计算性能与成本、接触检测精度及相互作用力定义、接触模型实现难度、表面平滑度、方法的通用性、所能处理的配置类型、同一方法下描述不同颗粒几何形态的一致性、输入数据预处理难度以及是否需要额外工具等方面。

多球形描述可能是应用最广泛的方法，因为它能够表示任意几何形状的颗粒。该方法非常适合再现玉米等形状复杂且表面较为光滑的食品颗粒的几何结构[8,9]。天然颗粒的三维结构是通过重叠的球体团块来生成的，这些球体可以部分相互渗透，从而使聚集体恰好适配所要复现的颗粒体积。该方法另一个优点在于，碰撞检测和力计算的算法实现相对简单。颗粒间的微观相互作用通过球体之间的接触来实现。由于球体碰撞的物理机制已被深入研究，因此在离散元模拟中存在大量经过验证的球体-球体接触模型。Wei等人[8]和Zhou等人[9]采用了粘弹性摩擦接触模型，这类模型的参数包括杨氏模量、泊松比、恢复系数以及滚动摩擦系数和滑动摩擦系数。接触力的法向分量和切向分量也取决于成颗粒的相互作用的球体对之间的重叠距离、法向相对速度和切向相对速度[8,9]。

准确描述颗粒表面对于精确再现颗粒堆积的诸多特性至关重要，如堆积形态、安息角、孔隙率分布[8]以及体积密度[9]。此外，接触模型的准确性也会显著影响颗粒的堆积行为，尤其是当该模型能够充分考虑静摩擦、滚动摩擦和恢复系数时。而颗粒集合体的堆积行为则受材料密度、杨氏模量和泊松比等机械特性的影响较小[8,9]。

多球形模型的通用性与简洁性也带来了相应代价。该方法难以处理尖锐的角边以及精细的几何细节，这限制了其精确度，可能导致结果与实验数据不符，或者需要增加构成颗粒的球体数量，进而大幅提高计算成本。计算效率不足是团簇球形模型应用范围受限的主要缺点。此外，多球形方法还需要额外工具来优化用重叠球体填充复杂形状颗粒内部空间，即便使用大量球体（通常每颗颗粒的球体数量不超过一百个），也无法保证达到所需的精度。

多面体型[[10], [11], [12], [13], [14], [15]]模型与基于网格的模型[16,17]虽有很多相似之处，但在模拟具有非凸角几何形状和尖锐边缘的复杂非球形颗粒方面具有很强的通用性。通过使用精细的三角形网格，还可以相当准确地近似光滑表面。不过，随着表面网格三角形面数的增加，计算效率会显著下降。此外，这类方法还需要网格生成工具来准备输入数据。

在岩土力学、地质工程和建筑行业中所关注的天然颗粒，如岩石、巨石、石英砂颗粒等，其形状极为不规则，具有尖锐的角边和粗糙的表面。对于这类颗粒，在模拟颗粒材料振动堆积密实化过程时，精确描述其形状在离散元建模中具有重要意义[10]。

尽管如此，由于接触点处碰撞法向量的定义存在模糊性，接触力的计算精度仍不及解析方法。Xie等人[13,14]提出了两种碰撞模型，用于在离散元框架内分析由不同多面体和超椭球形颗粒组成的复杂混合物的动态行为。他们假设先将超椭球体初步转换为多面体。第一种模型采用重叠体积法，第二种则采用最深点法。随着多面体面数的增加，计算效率会迅速下降[13,14]。该方法的性能还取决于将光滑颗粒（如球体、超椭球体等）转换后的多面体面数以及光滑颗粒的数量[14]。多面体模型能够描述含有光滑颗粒和角状颗粒的混合体系的颗粒几何结构。Xie等人[14]就在半自动磨机中使用球形介质和多面体矿石进行研磨过程的离散元模拟中运用了这一优势。在描述易碎颗粒方面，多面体技术也取得了良好效果，例如在模拟铁路道砟堆积的离散元模拟中[15]。

基于网格的表示方法被用于描述任意形状颗粒的表面，并已成功集成到离散元模拟中[16]。Zhan等人[16]提出了基于网格方法中最完善的接触力模型之一。该模型的接触检测与力计算算法如下：将网格节点视为颗粒表面的接触点，并为其赋予相应的单位法向量。通过遍历所有点，找出位于相邻颗粒内部的点（即活跃接触点），从而实现接触检测[16]。两个相互作用颗粒之间的接触力由相互穿透深度决定，该深度可通过属于不同颗粒的最近活跃接触点之间的距离来估算。在线性弹性接触模型中，法向相互作用力由阻尼/恢复系数、相对法向速度和接触刚度共同决定，而切向相互作用力则取决于摩擦系数、相对切向速度和接触刚度。

基于网格方法的一个显著优点是，可以直接利用表面网格计算颗粒的质量-惯性属性，如质量、质心、转动惯量等。此外，基于网格的描述方法常用于表示岩石颗粒和沙粒的复杂形状，而这些形状是通过球谐分析重新构建的[17]。近年来，基于计算机断层扫描图像的实时颗粒形状重建技术的发展，使得能够考虑包括表面粗糙度在内的多尺度形态特征，从而再现真实不规则颗粒的剪切强度行为、膨胀行为和破碎行为[[24], [25], [26]]。通过开发包含真实沙粒形状的离散元模型，现已能够根据沙粒的形状、大小和含量，再现沙土混合物的复杂行为[[27], [28], [29]]。

为了解决少数非标准问题，人们还开发了特殊方法。Wang等人[18]提出了一种用于重建构成颗粒介质的COVID-19病毒颗粒复杂形态的方法。该方法采用基于所谓体素（三维像素）的表示方式。通过对颗粒构型和堆积空间的数字化处理[19]，显著降低了以下方面的复杂性：任意形状颗粒的数学描述难度、碰撞和重叠检测算法的复杂性，以及基于体素的模拟的计算成本。该方法允许直接使用由扫描电子显微镜生成的3D生物结构图像中的数字体素3D数据。通过使用三维离散矩阵（立方晶格）来描述颗粒系统，可以实现対带有刺突蛋白的复杂形状病毒颗粒的建模。各个被离散化的病毒颗粒通过将零值元素改为1值而位于该三维矩阵中。动态自堆积过程则是通过不断调整每个任意形状颗粒的位置来实现的。颗粒会沿着随机选定的、未被相邻颗粒阻挡的方向移动一个单元格的距离。为了考虑重力作用，会提高颗粒向下移动的概率。这种数字堆积算法能够很好地预测任意几何形状容器中复杂形状颗粒的堆积密度[[18], [19], [20], [21]]。该方法还能帮助确定复杂形状不规则颗粒的堆积比例与其球形度、圆度及尺寸之间的关联[22]。不过，基于体素的方法[[18], [19], [20], [21], [22]]无论颗粒形状多么简单或复杂，都需要大量体素（每颗颗粒需数万个体素）才能呈现出所需的颗粒结构。尽管这种方法具有很高的灵活性，但仅适用于极其复杂的非规则物体建模。从计算机断层扫描得到的数字化结构可直接作为离散元模拟的输入数据。此外，该方法的碰撞模型选择较为有限。接触力预测算法通常是基于重叠体积的计算，这类方法在模拟堆积行为时精度一般，需要经过仔细的拟合和校准。上文提到的传统方法（如多面体法、三角形网格法、体素化技术等）都采用直观简单的模型来表示颗粒的几何形态。而若运用更复杂的数学理论，则能进一步提升离散元方法的功能。通过对颗粒表面的解析描述，可以精确确定接触点处的曲率及表面法线，这对于计算接触力而言非常重要。隐式表面模型[30]则是利用隐函数对颗粒几何形态进行解析表达，颗粒表面由函数值为零的空间中的点集来描述。虽然该方法能高精度地再现各种曲线表面，但一旦出现细微的不规则性，方程的复杂性就会急剧上升，从而导致离散元模拟的计算成本大幅增加。隐式表面模型[30]在模拟轴对称、非球形且具有封闭旋转面的颗粒系统时应用效果良好。然而，用隐函数描述颗粒形状存在诸多弊端：很难找到能定义出所需形状表面的隐函数；当需要对几何形态进行微调时也会面临很大困难；此外，对用隐函数定义的几何对象进行矩阵变换（平移、旋转、缩放等）也颇为棘手。对于椭球形[[31], [32], [33]]、超椭球形[34,35]以及多超椭球形[36,37]等复杂颗粒表面的解析表达同样基于隐函数。这类技术非常适合用于模拟由人工制造的颗粒所构成的散装材料的动力学行为。与天然颗粒不同，人工颗粒的几何形状是预先设计好的，旨在为颗粒材料赋予特定的宏观特性。基于椭球形的方法能够模拟颗粒的无序堆积行为，因为这类颗粒的形状、大小以及扁平度/拉长度和长宽比都很容易被控制或调整[[31], [32], [33]]。基于椭球形的方法不仅在描述复杂颗粒的几何形态方面精度很高，在确定接触力方面也同样出色。因此，使用椭球形模型能够很好地再现颗粒系统在各种边界条件下的宏观特性（如孔隙率、堆积密度和取向有序度[[31], [32], [33]]）以及微观特性（例如接触力、呈长之字链状的力网络[31]）。在离散元技术框架下，以超椭球形作为颗粒的数学表示形式，大大拓展了椭球形方法的适用范围，可描述更多种类的几何形状[34,35]。超椭球形颗粒的表面可以通过类比椭球形表面的通用方程进行解析表达，这些方程能够涵盖多种非球形颗粒的形态。而且，超椭球形颗粒还可以像多球形技术那样进行组合。与多球形颗粒模型相比，该方法在精度和计算时间方面都具有显著优势，因此很适用于制药领域[36,37]。刘和赵[36]以及刘等人[37]在利用离散元模拟片剂和胶囊时，采用了该模型来描述多种颗粒几何形态，这些复杂结构可以被看作是由相互重叠的超椭球形颗粒构成的团块。对颗粒表面的解析表达使得能够精确确定碰撞点附近的法向力和切向力分量。这些方法还提供了多种经过验证的相互作用力模型，可用于精确模拟颗粒介质的堆积和流动行为。使用椭球形/超椭球形颗粒的方法的计算效率高于基于多球形颗粒的技术[34,37]。不过，超椭球形模型仅适用于特定类别的颗粒形态，这就限制了该方法的应用范围——比如对于具有平面、尖边和尖角的颗粒几何形态，即便超椭球形描述能够覆盖相当多的颗粒形状，该方法仍难以适用。基于参数的曲面方法是描述复杂表面的有效工具。该方法利用依赖于两个变量的多项式函数来确定物理空间中参数曲面的点集。贝塞尔[38]提出了相应的方程，其系数代表了空间中的控制点，这些控制点能够精确界定复杂曲线或曲面的形状。贝塞尔曲面的方程是通过伯恩斯坦多项式来确定的，而这些多项式被视为控制点之间的插值函数。贝塞尔曲面上的任意一点都是定义它的各控制点的加权之和。然而，基于贝塞尔曲面的方法也存在诸多缺陷：很难精确控制曲面的几何形态，任何一个控制点的位移都会对整个曲面产生全局影响；为满足一定的几何约束条件，往往需要使用高阶贝塞尔曲面，这会降低基于贝塞尔曲面运行的算法的效率和稳定性。此外，该方法在形状表达方面也存在明显局限，无法准确描述尖边和尖角。贝塞尔曲面的阶数不仅会随着控制点数量的增加而上升，还会随着所描述几何形态的复杂程度提升而增加。另外，该方法也无法覆盖工业应用中常见的许多几何形状。为了避免使用高阶贝塞尔曲面，可以采用由平滑连接的低阶片段构成的曲面，这些分段低阶片段是基于B样条构建的。B样条是贝塞尔曲线的扩展形式，它同样可以通过控制点来定义曲线的形状，但其阶数并不取决于控制点的数量。B样条方程的形式与贝塞尔曲线类似，只不过其插值函数并非伯恩斯坦多项式，而是根据参数值递归定义的。利用B样条基函数得到的参数曲线和曲面被称为非均匀有理B样条（NURBS）。虽然参数曲线和曲面最初是为计算机辅助设计领域开发的，但安德拉德等人[39]以及刘等人[40]将基于B样条的颗粒形状表示方法应用于离散元方法中，用于模拟颗粒力学行为。安德拉德等人[39]将离散元方法与基于B样条的二维颗粒形态描述方法相结合，从而再现了天然沙子的复杂颗粒形状。刘等人[40]进一步改进了安德拉德等人提出的方法，在颗粒形状描述和接触处理方面取得了更好效果。他们的接触算法是基于对两条任意B样条曲线之间交点的显式计算，该方法还包含了离散元建模的所有标准流程。尽管基于B样条的离散元方法[39,40]能够精确再现几乎任何二维几何形状，但其计算成本非常高。将该方法扩展到三维几何形态后，上述问题会更加严重[41]。曲面与曲面之间的接触检测是一项极具挑战性的任务，相应的三维算法相比曲线与曲线之间的接触检测要复杂得多，计算成本也更高[41]。因此，检测颗粒间接触的难度必然会导致模拟过程的巨大计算开销。此外，该方法也无法完全消除颗粒形状方面的限制。现有研究表明，要构建一种高精度、计算效率高且应用范围广、可直接用于模拟任意形状颗粒堆积行为的模型，目前还远未解决这一问题。因此，开发一种能够避免以往方法缺陷的新方法显得极为重要。NURBS技术能够精确地表示现实世界中颗粒的复杂几何形态，有望解决这一难题。本研究的主要目标就是克服现有障碍，开发一种基于NURBS的通用方法，用以精确表示现实中复杂的颗粒形态，进而用于在离散元方法中模拟颗粒在任意形状容器中的堆积行为。本研究的创新之处在于，能够在保持合理计算成本的前提下，将基于NURBS的离散元建模能力扩展到三维任意形状的非凸颗粒。这项研究弥补了以往仅针对二维形态的研究向三维复杂几何形态扩展方面的不足。与之前的二维研究相比，本研究采用了三维NURBS曲面来直接对任意形状颗粒进行数学描述，并优化了三维非凸几何形态的接触检测功能。本文的结构如下：第2部分阐述了基于NURBS的任意三维颗粒表面重建的数学基础，同时介绍了碰撞检测算法的实现细节以及相互作用力的计算方法。第3部分则用于对该方法进行实验验证，通过实际应用该模型，展示了其在模拟置于任意形状容器中的复杂颗粒的堆积行为方面的应用潜力。第4部分会对该模型的功能、属性和特点进行讨论。最后第5部分会对所提出方法的计算效率进行评估，具体是从颗粒数量和碰撞次数两个方面来进行分析。在颗粒配置的复杂性和多样性方面，所提出的方法可与以往基于参数曲面的离散元方法以及作为多边形表示特例的基于网格的离散元模型相媲美。在颗粒离散元建模过程中，碰撞检测是计算成本最高的步骤。

结论：本文提出了一种基于NURBS的通用离散元方法，可用于处理具有任意几何形状、结构极为复杂的3D颗粒的随机堆积问题。该方法采用分段多项式数学表达式来精确描述颗粒的自由形态。同时，文中还提出并测试了一种用于计算颗粒质量-惯性特性（如体积、质量、质心以及各向惯性主矩）的技术。此外，还介绍了一种处理碰撞的方法。

CRediT作者贡献说明：Alexandra L. Zheleznyakova：写作——审稿与编辑，写作——初稿撰写，可视化，验证，软件，方法论，研究，概念设计。

利益冲突声明：作者声明不存在任何可能影响本文研究成果的已知财务利益或个人关系。

致谢：本研究是在俄罗斯科学院基础研究政府计划的支持下完成的（合同编号：124012500442-3）。