基于最小距离势能的非连续变形分析方法及其并行实现

《Engineering Analysis with Boundary Elements》：A discontinuous deformation analysis method based on minimum distance potential and its parallel implementation

　　

摘要 不连续变形分析（DDA）是模拟离散块体系统运动与变形的重要数值工具。然而，其传统的隐式罚函数方法存在接触处理复杂、计算成本高的问题。为解决这些问题，本文提出了一种基于最小距离势的接触力计算方法，并实现了一个利用并行计算的DDA框架。该方法通过定义块体内任意一点到其边界的最小距离来构建接触力势场，从而实现对凸多边形和凹多边形接触问题的统一处理。此外，本文提出了一种适用于接触力势计算的空间半检测法，用于快速接触检测和排除不相交的边对。数值算例表明，所提方法能够准确模拟凹多边形之间复杂的咬合行为，并行计算对大规模块体系统表现出显著的加速效果，半空间检测法有效加速了接触计算过程。数值结果表明，对于显式接触处理的隐式时间积分，建议将接触罚刚度设置为低于材料的弹性模量，并随时间步长进行调节，以维持数值稳定性。

引言 由Shi提出的不连续变形分析（DDA）已成为模拟离散块体系统运动与变形的一种关键数值工具。该方法在岩石裂纹扩展、地下工程稳定性、边坡失稳机理以及地震响应分析等方面展现出显著优势。其核心方面——包括接触检测算法、接触力学模型——的持续改进，进一步拓展了DDA的工程适用性。

然而，与同样属于离散数值方法的离散元法（DEM）不同，DDA采用隐式分析方法，能够以相对较大的时间步长实现无条件稳定计算。为处理接触非线性，Shi提出了一种"开-闭迭代"技术。若迭代无法收敛，通常采用经验策略减小时间步长并重新计算该步。尽管如此，隐式形式需要组装并求解全局刚度矩阵，且频繁的开-闭迭代导致计算成本高昂。Zheng等人利用互补理论成功避免了开-闭迭代问题；Fan等人进一步提出了DDA的混合线性互补形式。然而，这些方法需要求解非线性方程，计算仍然昂贵。其他研究者从不同角度出发试图缓解开-闭迭代的影响：Zheng等人发展了一种基于力的DDA形式，而Yu和Yin则完全放弃了开-闭迭代，仅依据初始接触状态计算每个时间步，并使用较小的弹簧刚度以最小化误差。

随着大规模离散块体系统模拟需求的增加，并行计算已成为提高DDA和数值流形法（NMM）计算效率的关键研究方向。DDA中使用的隐式罚函数法需要在全局刚度矩阵中组装非对角子矩阵以考虑接触力，这要求形成全局刚度矩阵，从而阻碍了DDA并行计算技术的发展。现有的大多数DDA并行策略都集中在刚度矩阵求解上，而对整个DDA流程全面并行的研究有限，这主要是由于全局矩阵组装过程中存在数据竞争问题。Wang采用显式罚函数法，以避免将接触力纳入全局刚度矩阵，允许独立的块级求解，从而显著降低了并行的难度。

显式罚函数法仍将接触力视为集中荷载，需要对不同的接触情形进行分别处理，并且必须妥善处理凸块体顶点-顶点接触可能引起的奇异性问题。尽管已有大量研究提出了改进方案，但仍未形成普适性解决方案。为克服这些挑战，最初由Munjiza提出并广泛应用于流形离散元法（FDEM）的接触力势方法被引入DDA。该方法不依赖穿透深度计算接触力，而是基于在每个块体内定义的势场计算力，无需区分复杂的接触类型，提供统一处理方式，极大地简化了接触检测过程。

Munjiza的接触力势方法最初仅限于三角形单元。为将其适配于DDA，Xu等人将块体细分为三角形单元进行力计算，Yang等人也采用了该方法，但计算成本较高。Zheng等人提出了针对整个凸多边形计算接触势的改进公式。Fan等人进一步提出了一种结合接触检测的广义接触势，将方法推广至凹多边形。然而，与Munjiza的原始方法不同，他们的方法仍需要在检测阶段检测和分类不同的接触情形。除DDA外，基于势的接触力方法也被推广至物质点法（MPM）和数值流形法（NMM）。特别是针对MPM，Zhou等人通过将接触力表示为保守势场作出了杰出贡献：通过由势函数的梯度导出接触力，其稳健框架优雅地解决了法向不连续问题，并保证碰撞过程中严格的能量守恒。通过在流形单元上定义势函数，接触力可直接由单元重叠区域计算得到。

在本研究中，我们定义了基于最小距离的接触力势，使其可应用于凹多边形。此外，由于接触力向量可直接由势导出而无需组装全局刚度矩阵，我们实现了DDA的并行化版本。Python因其开发迅速和可扩展性而成为最受欢迎的编程语言之一，尽管由于其计算速度受限而很少用于高性能数值计算。尽管如此，我们使用Python开发了一个并行DDA程序，并在本文给出了相应的伪代码。

本文结构如下：第2章简要概述DDA的基本理论；第3章介绍最小距离势及其在凹多边形上的计算；第4章描述基于Python的并行实现；第5章通过数值算例验证所提方法；第6章给出结论与未来展望。

DDA简要回顾：在DDA中，块体作为最小计算单元，任意块体的位移由六个位移分量表示：d=(u0,v0,r0,εx,εy,εxy)T，其中u0、v0为块体质心的刚体位移，r0为块体绕质心的旋转角，εx、εy、εxy为块体的应变分量。块体内任意点(x,y)的位移可由变形分量d表示。

基于最小距离势的接触力计算方法：本节首先简要回顾Munjiza提出的接触力势方法，随后介绍本文采用的最小距离势及其在DDA中的计算实现。

逐块计算方案：公式(3)展示了经典DDA中线性系统的组装与求解过程。在该公式中，每个块体的刚度矩阵Kii位于主对角线上，而非对角项Kij则使用隐式罚函数法构建。本文采用的接触力势方法属于显式罚函数法：它直接计算接触力而不将其组装到全局刚度矩阵中。

数值算例：本节首先通过数值算例验证所提最小距离势方法对凹多边形的适用性。随后，采用滑块模型分析在隐式时间积分中使用显式接触力所引起的数值稳定性问题。最后，通过随机颗粒堆积模型评估不同块体数量下的计算效率。

结论：为解决传统DDA方法在接触处理和计算效率方面的局限，本研究提出了一种基于最小距离势的接触力计算模型，并结合并行计算技术，实现了块体系统的高效稳定模拟。主要结论如下：
（1）通过定义点到多边形边界的最小距离势函数，所提方法避免了接触模式识别问题……

作者贡献声明：
Keqin Zhang：审阅与编辑、原稿撰写、软件、方法、构思。
Wei Wu：审阅与编辑、监督、调查、资金获取、构思。
He Huang：监督、资金获取、数据管理。
Guangmao Mao：审阅与编辑、监督、方法、构思。
Yongtao Hao：软件、数据管理、构思。
Hehua Zhu：监督、资金获取、形式分析、构思。

Keqin Zhang | Wei Wu | He Huang | Guangmao Mao | Yongtao Hao | Hehua Zhu