一种用于脆性材料和层状结构断裂分析的并行加速梁-粒子模型的开源C++实现

《Engineering Analysis with Boundary Elements》：An open-source C++ implementation of parallel accelerated beam-particle model for fracture analysis of brittle materials and laminated structures

　　

# 摘要

本文提出了一种二维梁-粒子耦合数值框架，用于精确捕捉脆性材料和层状工程结构的力学响应、渐进性断裂和灾难性破坏。该框架采用欧拉-伯努利梁单元表示粒子间的粘结作用，结合力-位移本构关系和应变-转动耦合破坏准则，模拟粘结损伤与断裂过程。排序-扫描算法优化了接触检测，相较于全局遍历降低了计算成本，可容纳多顶点多边形粒子，并以更鲁棒的方式确定不同构型下接触力的方向和大小。隐式后向差分预测-校正方案增强了数值稳定性，基于OpenMP的CPU并行计算提升了大规模系统的计算效率。通过岩石拉压、不同骨料含量混凝土、混凝土三点弯曲以及多接触层状结构的压缩试验验证了该模型。结果表明，该模型再现了骨料含量对混凝土刚度、强度和断裂模式的影响，以及层状结构中的脆性弯曲断裂和渐进性界面滑移。文中讨论了二维模型的尺寸局限性和三维建模的技术瓶颈。改进后的梁-粒子模型在所考虑的基准试验中表现出满意的精度和计算效率，为岩土工程和结构工程中的损伤机制分析与稳定性评估提供了可靠的数值支撑。

# 引言

脆性材料（如天然岩石、混凝土、陶瓷和砌体）及层状复合结构是土木、水利、交通、能源和岩土工程中的关键组成部分[1],[2]。这些材料广泛应用于大坝、隧道、桥梁、核电站基础、边坡支护系统和地下储能设施等重大工程设施中。在长期静力或突发性动力荷载（如地震、爆破、冲击和疲劳）作用下，这些材料和结构容易发生细观缺陷萌生[3]、裂纹扩展、损伤累积和贯通[4]，最终导致整体灾难性失稳（图1、图2）。其渐进性断裂过程、损伤演化和失稳机制是工程力学、计算力学和岩土工程等交叉学科领域的核心科学问题，这些问题直接影响重大工程项目的长期安全与风险防控。

由于岩石、混凝土和砌体等典型脆性材料具有固有非均质性、多相复合结构以及复杂分布的细观孔隙、界面和缺陷[5],[6],[7]，其断裂行为表现出强非线性、路径随机性和不可逆过程特征。传统连续介质力学方法（如有限元法[8],[9]和有限差分法[10],[11]）依赖连续介质假设，在模拟整个裂纹过程（从萌生、扩展到分叉、贯通及最终结构失稳）时存在固有局限性：(i) 在大变形和断裂区域容易发生后处理网格畸变，导致计算过早终止或精度急剧下降[12]；(ii) 裂纹追踪与表征依赖于网格拓扑，使得复杂多路径裂纹的真实模拟受到阻碍[13]；(iii) 无法准确反映细观骨料-基体界面、粒子接触和层间滑移等不连续效应，导致对材料真实断裂机制的描述出现偏差，难以满足脆性材料和层状结构高保真断裂模拟的需求[14]。

离散数值方法天然适用于模拟不连续介质的力学行为[16],[17]，已成为克服上述局限性的核心技术途径。Bolander等人系统回顾了应用于混凝土材料和结构的各类离散力学模型，特别关注了基于粒子的格子模型。该文献从模型分类、域离散化方法、弹性与断裂行为表征、动力加载、钢筋模拟和多场耦合分析等多个方面比较了不同离散方法的特点和适用性。还总结了它们相对于连续介质模型的优势和局限性，最后呈现了离散模型在工程应用中的现状[18]。其中，梁-粒子模型（BPM）因其独特的建模优势脱颖而出[19],[20],[21]。该模型将离散粒子表征与梁单元的粘结刚度特性相结合，采用离散粒子描述材料的细观基体、骨料或结构块体，利用梁单元传递粒子间的拉、压、弯和剪力相互作用。这种能力使得对内部粘结损伤演化、界面滑移与分离以及多体接触等关键行为的模拟更加高效，在建模脆性材料的断裂模式、损伤进展和结构失稳方面相较于基于连续介质的方法具有明显优势。

在离散元和基于粒子的数值模型发展历程中，早期离散元方法（DEM）采用圆形圆盘或球形粒子作为基本单元，通过粒子间接触相互作用模拟不连续介质行为，但仅能表征粒子间的接触压缩[22],[23]，无法有效反映材料内部实际存在的抗拉、抗弯和抗剪粘结能力，因此与脆性材料的力学性质和断裂模式存在较大偏差。为弥补这一核心缺陷，全球研究者相继引入梁单元[4],[24]、弹簧单元[25]和粘结线段[26],[27]来描述粒子间粘结效应，逐步发展出多类耦合粒子模型。作为代表性模型，梁-粒子模型通过粒子与粘结梁的结合，能更真实地反映脆性材料的抗拉强度、弯曲效应、断裂韧性和复杂断裂模式，显著增强了基于粒子的模型对非均质脆性材料的表征能力。

尽管传统梁-粒子模型已取得重大进展，但在应用于大规模工程仿真和高精度断裂分析时仍面临四个关键技术瓶颈：(i) 接触检测算法效率低下——传统全局遍历搜索的时间复杂度为O(N2)，导致大规模粒子系统或多边形粒子场景下计算时间急剧增加，难以满足高效模拟需求；(ii) 缺乏多交点接触几何判定与力学计算的统一且精确的求解框架，导致接触力计算精度和稳定性不足；(iii) 大多数主流模型采用显式积分方案，数值稳定性差且临界时间步长严重受限，导致迭代次数多、大规模模型求解效率低；(iv) 大多数模型未实现高效并行加速，单线程计算难以支撑工程尺度大规模模拟，限制了适用性扩展和算力适配能力。

近期在多边形粒子建模[28],[29]（图3）、高效近邻搜索[30]、隐式积分和并行计算[31],[32]方面的进展为提升基于粒子的模拟精度和效率开辟了新的途径。特别是通过质心Voronoi镶嵌（CVT）生成的多边形粒子，相比传统圆形粒子能更好地表征材料的真实细观形态。排序-扫描算法将接触检测复杂度从O(N2)降低到接近线性水平，隐式后向差分预测-校正方案改善了数值稳定性并允许显著更大的时间步长。此外，基于OpenMP的多核并行化具有部署简单、硬件兼容性强和计算加速效果显著的优点，无需复杂集群配置即可实现高效的CPU级并行，为大规模梁-粒子模拟提供了低成本、高可行性的加速方案。

目前已有多种开源数值工具可用于模拟颗粒和准脆性材料，包括GranOO和Yade-DEM（显式DEM框架）以及OOFEM（有限元求解器）。与这些工具相比，本文提出的模型具有以下特点：(i) 混合粒子-梁公式，可捕捉超越标准接触力的拉、弯和剪相互作用；(ii) 用于精细粘结断裂的应变-转动耦合破坏准则；(iii) 适用于多交点的多边形粒子统一接触算法；(iv) 隐式后向差分积分方案，提升数值稳定性并允许更大的时间步长；(v) 基于OpenMP的并行加速。据我们所知，现有开源DEM实现中没有任何一个在单一框架中整合了所有这些方面。OOFEM基于连续介质，属于互补而非直接可比较的对象。

当前研究在接触算法、积分方案、粒子离散化和并行策略等各个方面均取得了渐进性突破。然而，将这些五项关键技术完整整合到统一的梁-粒子模型框架中——并针对岩石、混凝土和砌体结构进行系统验证——的开源实现仍然稀缺。现有模型大多存在功能碎片化、精度与效率难以平衡以及开源可扩展性不足等问题，因此无法为脆性材料和层状结构断裂机制研究提供通用、可靠且高效的数值工具。

为应对上述研究空白和工程需求，本文开发了一个完整的、开源的、二维并行加速的梁-粒子模型框架，采用C++编程语言实现，专门面向混凝土和岩石等脆性材料的渐进性断裂与损伤演化分析，以及砌体类层状结构的研究。该框架的核心技术创新包括：采用欧拉-伯努利梁单元精确表征粒子间的力学粘结作用，建立应变-转动耦合破坏准则来描述粘结梁的断裂和损伤行为；引入排序-扫描算法实现多边形粒子的高效接触检测，提出多交点接触条件下的统一接触力计算方法以保证不规则粒子接触计算的准确性；采用隐式后向差分预测-校正积分方案改善数值稳定性并增大时间步长，从而克服显式积分的效率瓶颈；以及基于OpenMP实现接触力和梁单元内力计算的多核并行加速，有助于减少当前模拟中中等规模模型的计算时间。总之，本工作首次在一个开源C++框架中完整整合并协同优化了此前分散的五项关键技术——高质量多边形粒子离散化（CVT）、高效排序-扫描接触检测、多交点接触力统一计算方法、隐式后向差分预测-校正积分以及基于OpenMP的多核并行加速（见表1）。与现有梁-粒子模型或DEM代码相比，本文提出的框架为工程尺度模拟提供了一个功能完备、可扩展且易于应用的开源工具。

为系统验证模型的精度和适用性，本文开展了一系列经典数值实验，包括斜长角闪岩单轴压缩和巴西劈裂、不同骨料含量混凝土单轴压缩、混凝土三点弯曲断裂以及多接触层状结构压缩，全面验证了模型在强度特性、变形行为、断裂模式和渐进性损伤演化方面的可靠性。本研究的目标是为所考察的工况提供具有合理预测能力且竞争效率的通用开源可扩展数值模拟工具，以揭示脆性材料和层状复合结构的断裂机制、损伤防控和工程安全评估，同时为梁-粒子模型未来向三维扩展、多场耦合分析和工程尺度灾难模拟奠定坚实的理论和技术基础。

需要强调的是，本文建模策略并非旨在精确复制真实的细观结构，因为在工程实践中这些结构往往未知且测量代价过高。相反，我们的目标是通过统计上代表材料本质几何特征的人工微结构来捕捉主导的宏观断裂现象学。这一理念与离散元范式一致，即使用简化的细观模型来再现控制结构破坏的关键物理机制。

## 章节片段

### 梁-粒子模型：物理公式与数值方法

本节依次介绍梁-粒子模型的理论公式、接触检测算法的实现、粒子间相互作用的计算、系统动力学的求解、数值积分方案的应用，以及总体并行计算流程和支持性计算环境的介绍。最后介绍……

### 数值模拟实验

本节基于梁-粒子离散数值模型，对岩石、混凝土和层状接触结构开展了一系列多条件力学试验和对比数值模拟。首先，进行斜长角闪岩实验室单轴压缩试验以确定试样的基本物理力学参数并分析其脆性剪切破坏特性。基于……构建梁-粒子模型。

### 局限性

本研究主要基于二维梁-粒子模型，系统分析了岩石和混凝土等脆性材料的力学响应和渐进性断裂机制。

需要指出的是，本二维框架采用合成微观结构（CVT）是一种有意的建模选择，旨在面向机理的参数化研究，其合理性在于此类合成组装体能够在统计意义上再现主导性的断裂现象。①**结论**在本文中，我们系统构建了一种面向脆性材料及层状结构力学分析的二维梁-颗粒计算模型，并完善了一个涵盖理论构建、数值算法与模拟验证的完整研究框架。混凝土及类砌体结构的力学响应与渐进断裂机理得到了有效揭示。核心结论如下：- 全面耦合的梁-颗粒模型 **CRediT作者贡献声明** 赵宁宁：撰写——审阅与编辑，撰写——初稿，可视化，验证，软件，资源，方法论，调查研究，形式分析，数据整理，概念设计。李晶：形式分析。**利益冲突声明** 作者声明不存在已知的可能影响本文所述工作的竞争性经济利益或个人关系。具体而言：- 本研究未获得任何资金资助。- 通讯作者赵宁宁与合著者李晶与任何具有经济利益关联的组织或实体不存在隶属关系或参与关系（如咨询费、教育资助、演讲人团队参与、会员资格、雇佣关系等）。**致谢** 本研究未获得公共、商业或非营利部门任何资助机构的专项拨款。赵宁宁 | 李晶