一种基于瓦赫普雷斯形状函数用于组合有限元–离散元方法的能量守恒光滑接触势

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：An energy-conserving smooth contact potential based on Wachspress shape functions for the combined finite–discrete element method

　　

摘要
本文提出了一种能量守恒的光滑接触势模型，用于在联合有限元-离散元方法中模拟任意凸多边形与凸多面体之间的接触相互作用。该方法基于Wachspress形函数构建接触势场，使接触力能够作为全局连续标量泛函的精确梯度被一致地导出。与传统依赖于分片子区域划分的势方法不同，新方法确保了势场的几何与力学连续性，消除了单元内部子域边界上的梯度不连续性。这有效地抑制了接触力大小与方向中的非物理振荡，降低了网格和局部几何特征的影响。该模型严格满足能量守恒，因为接触力源于与路径无关的势能量泛函。势场的固有光滑性还使得无需进行子单元细分即可直接进行边界积分，从而提供了流程简洁且计算高效的方案。通过一系列基准问题的验证表明，该方法在数值稳定性、精度和鲁棒性方面相较传统方法均有显著提升。

引言
接触相互作用是多体系统中最基本且最具挑战性的问题之一，广泛出现在固体力学[1,2]、颗粒系统[3,4]、断裂传播[5,6]和冲击动力学[7,8]等领域。接触直接介导力的传递与能量的演化，使系统具有固有的非线性和非光滑特性。在大变形且拓扑结构持续演化的情况下，这一问题尤为突出，频繁的接触状态切换和持续变化的接触区域及其几何关系对数值方法的稳定性和一致性提出了严格要求。因此，开发既能保持变分一致性又能保证计算效率的稳健接触算法，始终是计算接触力学中的关键问题。

现有的接触方法主要分为非光滑动力学方法[9,10]和罚函数方法[5,11]。非光滑方法通过拉格朗日乘子或变分不等式施加非穿透约束，在原则上可实现精确的约束满足，但会导致鞍点问题，并在接触状态频繁转换时面临严重的收敛困难[11,12,13]。罚函数方法通过引入接触势将不等式约束重构为基于能量的框架，避免了复杂的约束求解过程，在显式动力学格式中具有更高的计算效率[7]。在该框架下，接触力由接触势对广义坐标的梯度导出[14]，这意味着势函数的数学性质对数值解的质量与稳定性起着关键作用。传统接触理论要求接触势至少满足C1连续性，以确保力的连续演化和动力学稳定性；否则，状态转换期间的梯度不连续性将导致力的突变和虚假振荡[15,1]。

尽管存在上述要求，现有的罚函数方法仍然依赖于分片构造，即在不同的子域上定义不同的势函数，并仅在界面上弱地施加连续性[16,17]。尽管易于实现，但这种方法牺牲了子域边界处的光滑性。这些界面处的梯度不连续性会引发虚假的力振荡[18,19]，在显式动力学中尤其有害，因为力的跳变会损害时间积分的稳定性[20]。已有若干补救措施被提出，包括基于归一化距离度量的连续势[21,22]、纳入局部几何特征的网格无关方法[17]以及基于采样积分的全局光滑势场[18]。然而，这些方法仍存在若干固有局限性。例如：(i) 在多体问题中，力可能随接触面积的增加而减小，违反物理一致性；(ii) 在整个域上严格的C1连续性仍未实现，界面和几何转换处仍存在不连续性；(iii) 密集的采样使得该方法在许多体系统上计算代价过高；(iv) 对于具有尖锐特征或近似平行面的多边形或多面体[23,24]，法向量不连续性使力的稳定评估变得复杂[19,25,26]，而传统势在拓扑快速演化期间仍易受域不连续性和奇点的影响[27]。

在大变形和断裂问题中，这些挑战尤为突出，因为材料从连续介质向离散碎片过渡，而在失效后的响应阶段，碎片之间的接触相互作用主导了系统行为。联合有限元-离散元方法（FDEM）有效地捕捉了这一过渡过程，并已应用于岩石断裂[5,28]、爆破[29,30]以及多物理场耦合分析[31]。在此背景下，对碰撞、滑动和能量耗散的精确建模对于预测可靠性至关重要，这使得接触运动学和构成本质的统处理成为核心关注点。从变分角度看，接触被视为系统能量上的不等式约束。Wriggers [1]和Laursen [15]指出，稳定解要求接触势至少满足二阶可微性，以确保力的平滑演化。然而，在大多数FDEM方法中，势通过子单元积分计算，当接触点跨越单元边界时会导致不连续性。正如Konyukhov和Irslinger [32]所指出，梯度连续性的缺乏会引发力的跳变和数值振荡，在动力碎裂中还可能违反能量守恒，导致人工能量损失或增加[20]。能量一致的接触理论已逐步发展以解决这些缺陷，通过几何量之间的变分关系直接从系统总能量构造接触力，从而消除了对传统穿透度量的显式依赖[33,34,35]。尽管这些方法强制实现了严格的离散层面能量守恒并显著减少了网格依赖性，但仍在势的光滑性与几何普适性之间面临权衡。

基于上述考虑，本研究提出了一种在严格变分框架下建立的光滑接触势模型。该方法基于Wachspress形函数[36]构建接触势场，能够适配任意凸多边形和凸多面体几何形状，且不对几何表示方法产生依赖。通过引入高阶连续性，所提出的模型在整个域上实现了势及其梯度的平滑过渡，有效消除了传统分片势在内部边界上因梯度不连续性而引发的力振荡，并减轻了跨越单元边界穿越所引起的力波动。至关重要的是，该模型严格满足能量和动量守恒条件，提供了稳定且物理一致的接触响应，同时显著降低了对网格离散化的敏感性。所提出的模型不仅适用于一般的接触力学问题，还特别为联合有限元-离散元方法量身定制，以应对连续向不连续过渡过程中快速演化拓扑下鲁棒的接触处理需求。

片段摘要
控制方程弱形式
FDEM在统一框架内整合了连续介质动力学与离散体动力学，能够模拟从连续变形、裂纹萌生与扩展到碎片化及块体分离的全过程。控制方程源于动量守恒：

ρü(x) ? ?·σ + b = 0，在Ωs中

其中Ωs ? ?^d表示计算域，ρ为材料密度，u为位移向量，x为坐标，σ为Cauchy应力张量，b为单位体积体力。

传统接触势
Munjiza引入了接触势场[5]，接触力由定义在相互作用体重叠区域内的势导出。二维域被离散化为三角形或多边形单元，三维域采用四面体或多面体单元。在二维中，势值取决于评价点所在的单元。理想情况下，相同的穿透深度应产生相同的势值。

能量守恒
传统势方法[14]可能在多边形单元内产生奇异区域，导致势分布不连续或不完整（图5(a)），从而在整个域上形成全局不一致的势场。相比之下，基于Wachspress形函数的所提出方法在多边形和多面体单元内产生连续且完整的势场（图5(b)），消除了内部不连续性并确保了几何……

动量守恒测试
本算例模拟两个相同四边形块的碰撞，以验证数值方法的动量全局守恒。如图15(a)所示，建立了两个边长为1m的正方形块，其质心位于同一水平轴上。块1位于左侧，以0.2 m/s的初始水平速度运动，而块2初始静止，两者之间留有1mm间隙。模型网格尺寸如……

三角形压头对对称分段多边形系统的影响
本算例研究块体与对称多块系统之间的碰撞[19]。图23(a)给出了详细的几何构型，其中系统的对称轴穿过顶部块的质心。模型网格尺寸如图23(b)所示。模拟中使用的材料参数列于表1。

在模拟中，顶部三角块以1 m/s的初始垂直速度运动，而下部块系统……

结论
本研究提出了一种针对任意凸多边形和凸多面体之间能量守恒接触相互作用的光滑接触势框架，特别强调其在联合有限元-离散元方法（FDEM）中的整合。该模型基于Wachspress形函数构建，天然保证了势场及其梯度在单元界面处的全局连续性。

与传统子单元势模型相比，……

CRediT作者贡献声明
Kunru Xie：写作—审阅与编辑、写作—原稿撰写、验证、形式化分析、数据整理、概念设计。
Xiaofeng Li：写作—审阅与编辑、监督、软件、资源、项目管理、方法论、调查、资金获取、形式化分析、数据整理、概念设计。
Haibo Li：监督、资源、项目管理、资金获取。
Dongbao Yang：验证、软件、方法论。
Zhonghao Li：验证、监督、形式化……

Kunru Xie | Xiaofeng Li | Haibo Li | Dongbao Yang | Zhonghao Li | Qi Zhao | Giovanni Grasselli