基于改进的最小距离接触力势方法的三维不连续变形分析方法

《Computers and Geotechnics》：Three-dimensional discontinuous deformation analysis method based on an improved minimum distance contact force potential method

　　张克勤 | 吴伟 | 王曦 | 朱俊峰 | 朱和华
同济大学土木工程学院，上海200092，中国

摘要
接触处理是三维非连续变形分析（3D DDA）中的核心问题，而接触力势方法在DDA中得到了广泛的应用。本文提出了一种改进的最小距离接触力势方法。首先，基于原始的最小距离接触力势，通过投影判断算法解决了凹多面体凹角处接触力势的奇异性问题。其次，构建了一种基于半空间理论的快速剔除的接触检测方法，以快速筛选无效的接触对，提高接触区域识别的效率。最后，通过六组数值示例验证了该方法的性能，并明确了惩罚参数与时间步长之间的正相关关系。研究结果表明，该方法能够有效适应任何凸多边体和凹多边体，并在接触计算效率上取得了显著提升。

引言
非连续变形分析（DDA）方法作为一种离散数值方法，自Shi（1988年）提出以来，已成为求解复杂块体系统变形、碰撞和滑动问题等重要工具。与连续介质力学方法（如有限元方法）不同（Wang等人，2025年），DDA无需假设介质的连续性，而是通过描述块体之间的接触约束和位移协调来准确地捕捉非连续介质的整体运动和局部相互作用（Zhang等人，2020年）。正是这种对非连续特性的固有适应性，使得DDA能够在众多工程领域得到广泛应用，包括边坡稳定性分析（Song等人，2025年；Li等人，2025年；Wu等人，2025年；Wu等人，2025年；Wu等人，2025年）、隧道开挖与支护设计（Wu等人，2018年；Wu等人，2020年；Zhu等人，2016年；Nie等人，2025年）以及地质灾害模拟（Gao等人，2025年；Li等人，2025年）。

接触计算是DDA的关键难点和性能瓶颈，尤其是在三维场景中这一问题更为突出（Wu等人，2025年；Wu等人，2025年）。经典DDA的接触计算基于“点-线-面”接触分类系统，需要逐一判断块体之间的各种可能接触类型，如顶点-面、顶点-边和边-面。对于不同接触类型，需要构建相应的穿透距离计算模型，这容易导致接触判断错误和遗漏（Zhang等人，2025年）。为了解决这一问题，许多学者进行了广泛探索：Bao和Zhao（2012年）以及Xia等人（2021年）分别提出了改进方法来处理二维DDA中顶点对顶点接触的不确定性。Wu等人（2014年）提出了多壳覆盖方法以加速接触检测，此后Wang等人（2019年；Wang等人，2020a；Wang等人，2020b；Wang等人，2020c）在此基础上提出了全局直接搜索方法，用于高精度检测不同凸块体之间的接触。

在离散元方法（DEM）中也进行了大量关于接触计算的研究，DEM与DDA类似。DEM引入了基于脉冲的方法，以计算速度换取了一定的精度（Je等人，2025年），同时赫兹接触模型也被应用于DEM（Shane等人，2019年）。一些研究将水平集函数（Level Set function）引入DEM，形成了LS-DEM（Duriez和Bonelli，2021年），水平集函数可以方便地表示任意几何形状、识别接触并计算接触力。基于相同原理，还开发了基于水平集的四维格子弹簧模型（4D-LSM）（Wei和Zhao，2025年）。值得注意的是，4D-LSM与DDA之间也进行了大量的耦合研究（Zhao等人，2016年；Zhao，2017年）。在接触处理方面，耦合的4D-LSM&DDA进一步引入了NURBS曲线以实现精确的DDA接触检测（Wei等人，2025年）。

接触力势方法最初由Munjiza（2004年）在有限离散元方法（FDEM）中提出，该方法将接触力的产生与重叠区域的势能存储相关联。由于FDEM最初仅适用于四面体元素（Munjiza和Andrews，2000年），其接触力势函数的适用范围有限。后续学者对其进行了扩展，提出了最小距离接触力势方法，将适用范围从四面体元素扩展到任意凸多边体（Liu等人，2020年；Zhao等人，2018年；Zhao等人，2020年）。该方法通过计算元素中任意点到每个面的最小距离来定义接触力势，简化了接触力的计算逻辑，无需区分特定接触类型，仅关注块体之间的重叠状态。一些研究尝试将这种最小距离接触力势方法引入DDA：Xu等人（2020年）将块体内的三角形元素分割开来计算接触力，Yang和Yang（2024年）也实现了这一方法；然而，这增加了计算成本。此外，Zheng等人（2020年）提出了改进的接触力势定义，实现了凸多边体整体接触力势场的计算。Fan等人（2018年）提出了一种结合接触力势方法和接触检测的通用接触力势，成功应用于凹多边体，但该方法仍需检测和区分不同的接触场景。

现有方法仍存在显著局限性：最小距离接触力势在凹多面体的凹角处表现出接触力奇异性现象，即势值为0，无法满足任意凸多边体和凹多边体接触分析的需求。为解决上述问题，本文提出了一种改进的最小距离接触力势方法，旨在实现任意凸多边体和凹多边体高效准确的接触计算。同时，优化了三维接触检测过程，并提出了一种基于半空间理论的快速剔除接触检测方法。

本文结构如下：第2节详细介绍了改进的最小距离接触力势的定义和实现逻辑；第3节阐述了该方法在3D DDA中的应用，包括3D DDA的基本理论、接触检测过程和接触力计算方法；第4节设计了六组数值示例以全面验证该方法的性能；第5节总结了全文的研究结论。

**改进的最小距离接触力势**
本节首先简要回顾了Munjiza提出的接触力势方法，然后介绍了本文使用改进的最小距离势。

**3D DDA的基本理论**
在DDA中，块体作为最小计算单元，任何块体的位移由12个位移分量表示：
\[d = (u_0, v_0, w_0, rx, ry, rz, \epsilon_x, \epsilon_y, \epsilon_z, \epsilon_yz, \epsilonzx, \epsilon_xy)^\top\]
其中 \(u_0, v_0, w_0\) 是块体重心的刚体位移，\(rx, ry, rz\) 是块体绕重心的旋转角度，\(\epsilon_x, \epsilon_y, \epsilon_z, \epsilon_yz, \epsilonzx, \epsilon_xy\) 是块体的应变分量。块体内任意点 \((x, y, z)\) 的位移可以用变形分量表示。

**接触力势的验证**
选择两个立方体块来验证接触力势方法的正确性。计算模型如图7所示：一个块的底部固定，忽略块体之间的摩擦，块体的弹性模量 \(E = 3 \times 10^6\) MPa，泊松比为0.2，密度为2000 kg/m3，整个计算过程中考虑了重力载荷，惩罚刚度设为 \(E \times 10^{-3}\)，时间步长为0.00001 s。

**结论**
本文重点关注3D DDA中接触计算的核心挑战，提出了一种改进的最小距离接触力势方法，实现了任意凸多边体和凹多边体接触行为的高效准确分析。主要研究结论如下：
(1) 成功构建了一种适用于任意凸多边体和凹多边体的接触力势模型。这种改进方案有效解决了接触力势的奇异性问题。

**作者贡献声明**
张克勤：撰写 – 审稿与编辑、撰写 – 初稿、软件开发、数据整理、概念化。
吴伟：撰写 – 审稿与编辑、监督、资源协调、项目管理、资金获取。
王曦：可视化、监督、资金获取。
朱俊峰：可视化、监督、资金获取。
朱和华：撰写 – 审稿与编辑、项目管理、方法论、正式分析。

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

**致谢**
本研究得到了国家自然科学基金[42272338, 41902275]和四川省高速公路建设发展集团有限公司[2024-cg-ky-17]科研项目的支持。