一种优化的颗粒破碎替代方法

《International Journal of Solids and Structures》：An optimized replacement method for particle breakage

　　吴萌萌 摘要 粒子破碎在决定颗粒材料的力学行为中起着关键作用。然而，现有的用于模拟粒子破碎的数值方法往往存在人为的质量损失、不真实的碎片形态以及不连续的力传递网络等问题，这些缺陷限制了其准确性和可靠性。为解决这些问题，本研究提出了一种基于功率图（一种广为人知的Voronoi剖分法的加权形式）的优化粒子替换方法，并将其应用于离散元方法框架中，以真实地模拟粒子破碎过程。所提出的方法采用基于Weibull最弱环节理论的粒子破碎准则，并根据实验测得的粒子尺寸分布生成碎片，从而确保质量守恒并保留原始粒子的形态特征。通过使用沙粒的一一映射进行离散元模拟，以及与高分辨率X射线计算机断层扫描技术的三轴剪切实验直接对比，验证了该方法的有效性。结果表明，该优化方法能够准确再现实验中观察到的宏观力学响应、粒子尺寸变化以及微观力学特性，同时有效克服了传统替换方法的局限性。这项研究为探究粒子破碎的微观力学机制提供了一种稳健且具有物理意义的工具，有望提升地质工程领域数值模拟的预测能力。

引言 随着土壤及其他颗粒系统力学行为研究的深入，粒子破碎的重要性日益凸显（如Bandini和Coop，2011；Lv等人，2020）。粒子破碎是指在外部应力作用下单个颗粒发生破碎或粉碎的现象，通常发生在施加的载荷超过颗粒强度时（如Einav，2007；de Bono和McDowell，2016）。这一过程通常由高围压、重复加载或强烈剪切引发，从而导致更细小颗粒的产生。粒子破碎会导致一系列问题，包括压缩性增加、抗剪强度降低以及渗透性变化，进而改变土壤的基本工程性质，可能影响地质工程的稳定性与性能（如Wang等人，2022；Li等人，2025）。在地质工程之外，粒子破碎在粉末技术领域也具有重要意义，诸如粉碎、研磨以及颗粒材料的压实等工艺在工业上有着重要应用价值（如Bu等人，2023；Li等人，2022；Petit等人，2025）。因此，结合地质工程与粉末技术的研究成果，可以更全面地理解不同颗粒系统中的破碎机制。

近年来，关于粒子破碎的理论研究与实验研究都取得了显著进展。在理论方面，人们提出了多种模型来描述不同加载条件下粒子破碎的机制与演化过程（如Wittel等人，2008；Gou等人，2022）。在实验方面，一维压缩试验和单颗粒破碎试验等实验室研究为了解影响粒子破碎的因素及其对土壤行为的影响提供了重要线索（如Yamamuro等人，1996；Meng等人，2022）。先进成像技术的发展，如X射线计算机断层扫描，进一步推动了粒子破碎研究，使得人们能够直接且无损地观察颗粒内部的结构变化（如Li等人，2019；Wang等人，2023）。例如，Zhao等人（2020）通过一维固结试验和单颗粒试验研究了粒子破碎的机制。在此基础上，Wu和Wang（2024）开展了三轴剪切条件下的粒子破碎研究，并提出了一种粒子追踪方法，从而能够更详细地分析破碎过程以及颗粒在变形过程中的形态演化。

这些先进的实验观察表明，颗粒集合体中单个颗粒的力学响应和破碎倾向深受其局部接触结构的影响，而局部接触结构通常用协调数来量化。文献中普遍认为，随着协调数的增加，颗粒会表现出缓冲效应——外部载荷在多个接触点上的分布降低了内部拉应力集中程度，从而提高了颗粒的表观抗压强度（Jaeger，1967；Tsoungui等人，1999；McDowell和Humphreys，2002）。相反，在无围压系统或协调数较低的稀疏颗粒流中，由于存在集中的点载荷，颗粒更容易发生灾难性破坏（Sukumaran等人，2006）。Cil和Alshibli（2014）通过3D X射线断层扫描提供了实验证据，证实接触的空间分布和数量对于决定单个沙粒的断裂模式和最大承载能力至关重要。

尽管实验室研究和成像技术为了解接触结构提供了重要线索，但它们往往难以捕捉内部力的动态、连续演化过程。因此，离散元方法的出现使得从颗粒尺度上研究粒子破碎成为可能，为了解颗粒材料中的破碎微观机制提供了宝贵见解（如Adesina等人，2024；Wu等人，2023）。为将上述接触结构效应纳入考虑，人们提出了多种数值方法，试图将协调数引入粒子破碎准则中。例如，Marketos和Bolton（2007）以及Ben-Nun和Einav（2010）提出了基于协调数的存活概率函数，研究表明，随着邻居数量的增加，颗粒破碎的概率显著下降。

目前，在离散元方法框架内模拟粒子破碎主要有两种方法：粘结颗粒法和粒子替换法。在粘结颗粒法中，每个颗粒被表示为由虚拟接触或粘聚力结合在一起的多个较小子颗粒的集合体（如Zhang等人，2025；Liakas等人，2017）。当子颗粒之间的粘结在载荷作用下失效时，颗粒就会破碎成更小的碎片。虽然这种方法能够真实地反映裂纹的逐步发展过程以及碎片的复杂形态，但它通常需要大量的计算资源，且需要对粘结属性进行精确校准，这限制了其在大规模模拟中的应用。另一种方法是粒子替换法，即一旦满足预定的破碎准则，就用多个较小的碎片替代原来的颗粒来模拟破碎过程（如Tavares和Anderson，2021；Zhou等人，2020；de Bono和McDowell，2016）。这种方法在计算上更为高效，也更容易实现，因此适合用于模拟大规模的颗粒集合体。例如，在粉末技术中，离散元建模已被广泛用于开发随机粒子替换策略，以模拟动态载荷下的破碎过程（如Zhou等人，2020）。然而，粒子替换法也存在一些局限性。首先，它无法捕捉详细的裂纹扩展过程，因为破碎是通过瞬间用碎片替代颗粒来表示的。其次，如果碎片的总体积或质量与原始颗粒不一致，这种方法可能导致质量不守恒。此外，颗粒的突然替换可能会在模拟中造成人为的应力下降，而这无法准确反映实际材料中破碎的渐进过程。另外，碎片的大小和分布选择也会显著影响模拟结果，需要仔细考虑才能真实反映破碎行为。

这些局限性的主要根源在于以往模型中采用的简化几何模型；历史上，大多数粒子破碎模拟都是针对球形颗粒进行的，采用的各种替换方案都是用一组预定义的较小球体来替代原始球体（如Tavares和Anderson，2021；Zhou等人，2020）。虽然这些方法在大规模模拟中计算效率较高，但它们不太适用于不规则颗粒，因为使用内切球体会导致明显的质量损失，也无法真实反映原始颗粒的角形形态。为了解决形状不规则的问题，基于Voronoi的方法已被引入一些商业离散元软件中（如Ghazvinian等人，2014）。然而，将标准的Voronoi破碎方法应用于高度不规则的离散颗粒集合体时，往往会生成极小的碎片，这大大降低了稳定的时间步长。为了保持计算效率，这些软件往往不得不人为删除小于最小尺寸阈值的碎片，这就又导致了质量损失。此外，它们内置的能量基准则主要是为动态冲击而设计的，并非针对准静态的地质工程载荷。因此，如何在不人为删除碎片的情况下实现严格的质量守恒，同时准确捕捉由应力引起的准静态破碎机制，仍然是一个重大挑战。

为在保证计算效率的同时克服上述局限性，我们提出了一种优化粒子替换方法。在该方法中，通过最大接触力准则来确定可能发生破碎的颗粒。与传统替换方法不同，我们的方法会将符合破碎准则的颗粒分割成多个细小颗粒，而这些颗粒的大小是根据实验中测得的粒子尺寸分布来确定的。这一策略不仅保留了原始颗粒形态的关键特征，还确保了质量守恒，有效解决了模拟中的应力下降问题。此外，由于避免了人为空洞的产生，新生成的细小颗粒的运动受到限制，从而能够模拟出实际材料中观察到的逐步破碎过程。为验证所提方法的有效性，我们建立了基于沙粒一一映射的离散元模型，并将其模拟结果与先前基于X射线CT的三轴剪切实验结果进行比较，证明了OPRM在捕捉粒子破碎行为方面的能力和准确性。

材料与实验 我们使用了Cheng等人（2020）设计的微型三轴装置，在上海同步辐射设施的BL13W光束线上进行了原位三轴试验（见图1）。该实验装置与传统的三轴系统非常相似，包括一个能够提供高达5 kN轴向力的加载系统、一个最大压力为2,000 kPa的Global Digital Systems控制的围压系统，以及一个数据采集系统。

粒子破碎准则 准确模拟粒子破碎的一个关键因素是选择合适的破碎准则。文献中提出了多种准则，包括最大拉应力准则（Lobo-Guerrero和Vallejo，2005）、平均主应力与最大主应力准则（De Bono和McDowell，2016）、最大接触力准则（Hanley等人，2015；Marketos和Bolton，2007）等。其中，最大接触力准则被证明特别有效。

数值模拟验证 为评估所提出的OPRM的有效性，我们对试样在三轴剪切作用下的宏观力学响应进行了对比分析。图5展示了偏应力与体积应变的演变情况，将实验结果与在2 MPa围压下采用传统替换方法和OPRM方法得到的结果进行了比较。如图5(a)所示，所有曲线的趋势相似，即在加载过程中偏应力迅速上升至峰值后...

结论 在本研究中，我们基于功率图开发了一种优化粒子替换方法，并将其应用于离散元方法框架中，用于模拟三轴压缩条件下颗粒材料的粒子破碎过程。所提出的方法采用了基于物理原理的破碎准则，并根据实验测得的粒子尺寸分布生成碎片，从而确保质量守恒并保留原始颗粒的形态特征。通过使用沙粒的一一映射离散元模拟，并与高分辨率X射线计算机断层扫描实验结果进行对比，验证了该方法的有效性。

CRediT作者贡献声明 吴萌萌：验证、监督、项目管理、方法论、研究、资金获取、正式分析、概念构建。

利益冲突声明 作者声明自己没有已知的财务利益或个人关系可能影响到本文的研究工作。

致谢 本研究得到了国家自然科学基金（项目编号：52408394）、香港特别行政区研究资助局的一般研究基金项目编号PolyU 15227624，以及地下工程智能建造与健康运维国家重点实验室（项目编号：SDGZ2504）的支持。