基于并行有限质点法的界面断裂-接触行为分析

唐敬哲，汪 伟，郑延丰,2，罗尧治,2

(1. 浙江大学空间结构研究中心，浙江，杭州 310058；2. 浙江省空间结构重点实验室，浙江，杭州 310058)

基于向量式力学所提出的有限质点法(Finite Particle Method，FPM)是一种面向工程应用、面向结构复杂行为分析的新型数值计算方法[1]。该方法的本质是对质点的动力平衡状态的描述，控制方程质点间相互独立，并以离散运动路径的方式处理几何非线性，在结构复杂非线性问题分析中具有优势和灵活性。目前该方法已被应用于机构运动[2−3]、结构倒塌破坏[4−5]、固体弹塑性大变形[6]、薄膜结构形态分析[7−8]、结构多尺度精细化[9]等结构复杂行为的分析。目前，有限质点法基于GPU并行加速的相关研究也已经初见成效，这种并行策略也已经应用在有限质点法通用计算平台的研发中[10]。借此，有限质点法的研究也逐渐向大型工程问题的精细化数值模拟的方向所发展。

界面断裂、接触等强非线性复杂行为广泛存在于结构破坏、连续倒塌、上下部结构协同等实际工程的精细化分析中。有限质点法对于强非线性行为模拟的适应性以及并行计算能力的优势使其以较低的成本完成界面行为的细观模拟。界面力学行为的模拟主要面临两个问题：一是界面作用对侦察；二是作用力计算。其中界面作用对侦查一向是最消耗时间和计算资源的步骤。目前的大多数的侦查算法都依照两个步骤进行，即先通过某种空间排序的手段将搜索区域局部化，再在局部范围进行精细的作用关系判断[11]。较有代表性的易于并行实现的侦查算法就是盒围法[12−13]。这种方法将空间点按某种排序算法进行空间排序，置于空间网格中，将侦察范围逐渐局部化。特别在将结构进行质点化离散的诸多数值方法中，如物质点法、离散元法和光滑粒子流体动力学等，大多使用这种方法的变体。如Chen等[14]使用桶排序的盒围法对有限元法与物质点法耦合的点-面接触进行了研究。Gan等[15]使用基数排序实现了离散元的颗粒系统的并行化接触侦察。Xia和Liang[16]基于四叉树邻域搜索对基于光滑例子流体动力学的盒围法进行了优化。Hu等[17]拓展出了条带盒围法进一步提高了搜索效率，成功应用在有限元与光滑例子流体动力学耦合的模拟中。尽管这些侦察算法具有易于并行实现的优点，也有不少效果显著的并行化研究[15− 16,18]，但具体的并行实现仍然是界面作用侦察算法计算效率问题的研究热点与难点。

在界面作用力计算方面，界面的接触与断裂行为采用不同的作用力模型。界面接触作为常见的界面行为，其常用的力学模型主要有两种，即适用于隐式求解的拉格朗日乘子法和通用的罚函数法，这二者都是通过对接触区域施加位移限制条件来达到接触处理的作用。显然，罚函数法适用于有限质点法的显式时间积分框架，因此有过一些成功的应用，例如喻莹和罗尧治[19]在离散化的梁杆系结构中，提出了针对空间点-线、线-线接触的侦测算法和弹性接触力模型。除了罚函数法，张鹏飞等[20]引入了显式的防御节点法，确保点-面接触过程中不会发生穿透，应用于薄壳结构的接触分析中。郑延丰[21]同样使用了防御节点法，并使用图形处理器(Graphics Processing Unit,GPU)对精细化接触行为的模拟进行了加速。对于界面断裂，学者们通常使用基于内聚力模型的各种变体[22− 24]来描述断裂过程中的裂尖应力。特别地，张鹏飞[25]基于此模型开发了可动态插入的有限质点法断裂单元模拟裂纹展开。这些实现充分利用了有限质点法将结构进行质点化离散后在强非连续问题分析中的优势，但也都缺乏普遍的适用性。

为了提高界面力学模型的广泛适用性，学者们在界面区域构造了特殊的单元，如接触单元、界面单元等，在局部范围内同时完成作用对侦查和作用力计算的过程，早期有Goodman等[26]的零厚度接触面单元、Desai等[27]的薄层单元等。在此基础上，诸多商用软件也出现了通用的单元种类，较有代表性的有LS-DYNA中的接触单元[28]、FLAC3D中的界面单元[29]以及ABAQUS中的内聚单元[30]等，但也存在若干局限性。例如，多以有限元和有限差分法为基础的商业软件，对于存在强非线性的局部界面行为的处理较为复杂，经常出现难以收敛的情况，且计算效率较低。另外，也缺乏统一的单元形式来模拟界面粘结、接触、断裂的耦合情况。

为了充分发挥有限质点法在处理强非线行为中的优势以及并行优势，本文提出了一种适用于有限质点法并行计算体系下的三角形界面单元，能够通用化地模拟界面之间的粘结、断裂、接触等界面作用行为以及之间的耦合。此界面单元基于GPU并行模式进行高效的界面作用侦察，使用内聚力模型和罚函数法作为界面断裂和接触力计算的力学模型，具有广泛的适用性。本文首先详述了其并行化的作用对侦察，并介绍了界面在粘结、断裂、以及接触这三种作用下的状态判据和作用力计算方法。最后，依托于有限质点法通用计算平台FPM[10]，通过数值算例的形式对此界面单元的计算效果进行了正确性和有效性验证。

1 考虑断裂-接触的界面单元

1.1 三角形界面单元

本文对于界面相互作用的模拟融合了粘结、断裂和接触这三种行为，为了区别于仅考虑接触的单元类型，本文采用了界面单元这一描述方式。本文将发生相互作用的边界面离散化为一系列的三质点三角形界面单元，根据可能发生的作用情况设置为目标界面单元与从界面单元。每个界面单元将其面积加权均分给其质点，即成为每个界面单元质点的代表面积。图1展示了目标界面单元、从界面单元以及界面质点的几何关系。

图1 有限质点法界面单元Fig. 1 FPM interface element

区别于其他有限质点法单元，界面单元质点并不具有独立的质量，而是通过与壳、四面体实体或六面体实体等其他有限质点法单元共用质点而为离散模型的边界面赋予界面作用模拟的能力。界面作用限定于从界面单元的质点(后文中将简称为从界面质点)与目标界面单元之间。每个时刻针对每一个从界面质点搜索与其实际发生作用的目标界面单元，生成一个个点-面形式的界面作用对，作用对的点-单元之间产生法向或切向的作用力，实现粘结、开裂、接触等行为的模拟。

出于罚函数法对显式数值算法的普遍通用性，本文选用罚函数法作为有限质点法界面模拟力学描述的基本形式，并在此基础上添加基于点-面自由度耦合的界面粘结模拟与基于内聚力模型的断裂模拟。界面作用对中的从界面质点p与目标界面单元et之间的界面力学模型如图2所示。如果考虑界面粘结的影响，质点p与粘结点pc保持共同运动，以加速度相同来计算粘结应力。当由于粘结所产生的界面等效应力超出了给定的临界值σcr时，则判定该界面作用对的粘结失效，转为内聚力开裂模型。裂缝完全展开后，继而转为接触模型。该接触模型等效于在作用发生点pc添加了沿法向界面弹簧Kn和切向非线性弹簧Ks，如图2所示。

图2 有限质点法界面单元力学模型Fig. 2 Mechanical model of the FPM interface element

本节将详细介绍界面作用力的并行化计算过程，包括作用对侦察、界面状态判断以及作用力计算这几个步骤。

1.2 界面作用对侦察

本文假定每一个参与作用对侦察的从界面质点在某一时刻仅与一个目标界面单元发生作用，因此作用对侦查的目标即为针对每一个从界面质点确定唯一的点-面作用对。在所有可能发生相互作用的点-面对中找到真正发生作用的点-面对是任何侦察算法中最重要且最耗时的一步。为了进一步加速侦察的效率，本文在经由GPU并行的有限质点法架构[10]中，进行并行化的作用对侦察。

在每一个时间步，首先对侦察范围进行区域分割，进行一次全局搜索。区域分割的目的是为了将庞大的侦察范围拆解为大量的可同时并独立地进行侦察的较小空间，方便并行执行以提高效率。参考在散粒系统中常用的“桶排序(Bucket Sort)”算法的思想，整个接触侦察区域根据坐标范围被划分为若干立方体网格，即若干个“桶”容器。若给定单元格边长l，可以计算在全局坐标系三个方向下的单元格个数Ni,j,k。以x方向上的网格数Ni为例：

式中：xmax和xmin分别为整个界面单元分布区域的包围盒x坐标的最大值和最小值，此包围盒尺寸在边界位置宜适当放大； f loor()为向下取整函数。这样，整个侦察区域即被分割为Ni×Nj×Nk个包围盒单元格，被分割区域内的任何一点均隶属于唯一的包围盒单元格中。

界面作用对侦察的过程分为两步。首先，将参与接触搜索的所有界面质点按照当前时刻的位置坐标装入划分好的包围盒单元格中，如图3所示。在某时刻的搜索开始阶段，对于任意界面质点p，根据其在上一时刻结束后的位置坐标可以确定其当前时刻的包围盒单元格归属Bp(i,j,k)，(i,j,k)表示该单元格在所有单元格中的位置索引。以x方向上的索引i为例：

图3 界面作用对侦察：全局搜索Fig. 3 Global search of interaction pairs

式中，xp为质点p在上一时刻结束后坐标。此过程以参与侦察的质点数为并行度，每一个GPU线程负责确定一个质点的包围盒单元格归属。

下一步，对于任意目标界面单元e，包络其三角形几何形状的包围盒单元格集合Be={B(i,j,k)}即为该单元的包围盒单元格集，被放入这个集合中的所有从界面质点即为潜在作用质点。这一过程以参与侦察的目标界面单元数为并行度，即每一个GPU线程负责确定一个目标界面单元的包围盒单元格集并将这些网格中的质点索引放入各个单元独有的数据结构中进行存储。至此，界面作用对的侦察从全局范围缩小至由各个目标界面单元的包围盒单元格集所控制的局部范围内，如图4所示。

图4 界面作用对侦察：局部搜索Fig. 4 Local search of interaction pairs

经过了全局和局部搜索的筛查，对于每个目标界面单元，需对其包围盒单元格集中的从界面质点进行逐一判断以最终确定实际发生界面作用的点-面作用对。如图5所示，为了保证界面作用方向的统一性，对任意一个界面单元e，其外法向ne定义三个质点逆时针围绕为正向，而任意一个界面质点p的外法向np取为与之相连的所有界面单元外法向的加权平均值：

图5 Inside-Outside算法示意图Fig. 5 Diagram of the Inside-Outside algorithm

式中：v12和v23代表单元e三个质点按照逆时针方向相连的边向量；Ae表示单元e的面积。利用式(3)可以唯一确定质点的法线方向，界面作用只发生在外法线方向相向的从界面质点与目标界面单元之间。

常规的点-面投影算法可以判断点的投影是否落在三角形内部以进一步判断点面是否发生相互作用，但在目标界面外凸或内凹时可能会发生重判或漏判。本文参考“Inside-Outside”判断法[31]来唯一确定每个从界面质点的目标界面单元归属。

如图5所示的质点p，对于单元e的三边分别判断其内外状态，计算：

式中，np为质点p的外法向。若dij≤0则认为质点p对于边ij处于In状态，反之为Out状态。如果质点p对于单元e的三边同时处于In或Out状态，则判定质点p归属于单元e。

在确定了质点的单元归属后，采用点-面投影算法计算质点在单元内的投影位置。考虑到三角形单元的形函数即为面积坐标，投影点的位置使用面积坐标进行表示：

式中，x1、x2、x3、xp分别为单元e的三个质点以及质点p的位置坐标。由式(5) 即可求得投影点pc在单元e中的面积坐标值Li。

得到了单元内投影点位置后，可以计算质点p对单元e的法向嵌入深度(或间隙)gn，gn取发生嵌入为负，产生间隙为正：

式中，npc为单元e内投影点(即界面作用发生点)的外法向，可以取单元e的外法向。为了将目标界面平滑化，本文将npc按照单元e三个质点的质点外法向按形函数(即面积坐标Li)进行单元内插：

这样处理可以使得界面单元内任意一点法线在整个目标界面内保持光滑连续，有效减小相邻界面单元作用力的突变。

上述过程以并行方式进行，每个线程负责对一个目标界面单元包围盒单元格集中的所有从质点依次进行循环，计算此从质点是否归属于该目标界面单元，如果确定归属则进而计算二者之间的法向嵌入深度gn。如果gn≤0，则发生相互作用，从而使得每个从质点的点-面作用对得以唯一确定。

下面讨论特殊情况。一个从质点可能同时处于若干目标界面单元的局部侦察范围。由于并行实现时各个线程之间具有同时性与独立性，不同线程针对不同单元的局部侦察过程可能同时对同一个从质点进行归属计算。如果某从质点的外法向恰好指向目标界面单元的质点或两个相邻目标界面单元的公共边，这种情况在界面初始粘结时十分常见，此时“Inside-Outside”算法会出现重判。这种情况可以借由并行算法中的原子操作进行处理。

当某个目标界面单元判定某从质点与之发生接触后，会将嵌入深度、单元索引、投影点面积坐标、投影点外法向写入该从质点独有的数据结构中。如果发生了重判，则会出现同一时刻不同线程(即不同目标界面单元)向同一个内存地址(即某从质点记录接触信息的数据结构)进行写入的操作。本文使用原子交换操作(AtomicExchange)[32]保证同一时刻只有一组接触信息成功写入该质点的数据结构中。此过程如图6所示。

图6 并行化的界面作用对侦察Fig. 6 Parallel procedures for searching interaction pairs

经过“Inside-Outside”算法以及并行原子写操作的双重保障，每个从质点的唯一点-面作用对得以最终确定，可以进行后续的作用力计算。

1.3 界面粘结状态判断

经过了界面作用对的侦察与最终确定，对于任意从界面质点p，如果它对单元e发生了嵌入，则p−e作用对成立，此单元标记为et。由于一个时刻界面作用对对于每个从界面质点是唯一的，界面作用力的计算将以界面质点为并行度执行，即每个线程负责计算一个界面作用对的作用力。

每个从界面质点均存在一个粘结状态指示量sbond，并在时间循环开始前根据界面设置情况设置为sbond=1 (初始考虑粘结)与sbond=0(初始不考虑粘结)。对于考虑粘结的情况，只需要首个迭代步中进行界面作用对的侦察与确定即可，这是因为粘结的存在会使得点-面相对位置保持不变。

对于界面作用对中的从界面质点p与界面作用发生点pc，结合有限质点法的点公式，有：

式中：下标p和pc分别为从界面质点与界面作用发生点；力向量的上标 e xt 、 int 和 b ond分别为质点的外力、单元变形所产生的单元内力以及界面作用的粘结力。界面作用发生点pc虽然并非在网格离散时生成的质点，但可将其视为et上为了抵御质点p的作用所产生的防御点，其质量、力等物理量可由三角形界面单元的形函数得到：

式 中：Li为pc在et中 的面积 坐 标；mi和Fi分别 为et的三个质点的质量与合力。

粘结作用力的计算原理为从界面质点p与界面作用发生点pc的加速度相同，辅以粘结作用力对于两质点互为相互作用力的关系：

1.4 界面开裂作用力计算

当界面等效应力 σeq超出了界面的临界应力σcr时，界面开始发生断裂。本文使用外加内聚力模型[33]对裂缝开展过程进行模拟，该模型中变形与应力的本构关系如图7所示。

图7 外加内聚力模型Fig. 7 Extrinsic cohesive model

建立接触点pc处的局部坐标系 (ξ,η,npc)，则p与pc的相对剪切变形gs为：

式 中：xp和xpc分别 为p与pc的 位移；xpc同样 在et内根据形函数插值求得。

类似于界面等效应力的定义，这里使用嵌入深度gn与相对剪切变形gs定义界面等效变形geq：

式中， max()函数是为了保证只有界面间隙计入等效变形而界面嵌入不计入。根据线性外加内聚力模型，开裂所产生的内聚应力 σcrack与等效变形之间geq的关系满足线性关系：

式中，gc为等效变形的临界值，可由断裂释放能Ec计算，gc=2Ec/σcr。

式中，Ap为质点p的代表区域面积。当geq＜gm时则视为卸载，此时用gm替换式(16)中的geq即可。当geq＞gc时，界面完全断裂，记scrack=0，界面内聚力为0，彻底转为接触模型。

1.5 界面接触作用力计算

式中：gn由式 (6) 计算得来；Kn为界面法向接触刚度；Ap为质点p的代表区域面积。

对于切向接触力，其计算通常基于摩尔库仑准则的基本形式，等效于弹塑性状态下的切向弹簧，屈服函数为：

式中，fsmax为界面可以提供的最大切向接触力，可由界面粘聚力c和摩擦角φ计算而出：

由于常用的双线性塑性形式的摩尔库仑准则在时间步长不合理或刚度设置不合理时可能会在相对剪切变形接近转向时难以收敛，本文使用反正切函数型的平滑摩擦力模型计算切向接触力，在每一个途径单元内部按照如下两个步骤进行：

1) 计算t时刻 (ξ,η) 平面内p与pc的相对剪切变形增量 ∆gs：

式中：vs=∆gs/∆t为相对剪切变形速度；C为无量纲的平滑参数，一般取0.1～0.001。可以看出，式 (21)实质为变刚度的摩尔库仑准则。

将以上所描述的界面作用对的侦察、确定与作用力计算的计算过程嵌入并行有限质点法的计算框架[10]中，即可到如图8所示的界面单元计算流程。

2 数值算例

本课题组在先前的研究中设计研发了采用GPU加速的有限质点法并行求解系统，并形成了有限质点法通用计算平台FPM。本节的算例均使用此平台进行加速计算。关于有限质点法的并行化以及FPM平台的相关信息，可参考文献[10]。

2.1 基础-土体接触

为了验证本文的界面单元在处理接触问题上的效果，本例考虑Qian等[34]使用的一个简化的基础-土体的接触模型。如图9所示，基础为6 m×6 m×0.6 m的立方体；上部结构的作用被简化为作用在基础上的法向均布荷载1 MPa；下部的土层简化为10 m×10 m×4 m的立方体，四周侧向约束，底部竖向约束。相互接触的两层三角形界面单元分别被设置在基础下表面和土层上表面，土层上表面设置为目标界面单元。

基础与土层采用相同的弹性材料，其弹性模量为100 MPa，泊松比为0.2。接触界面摩擦角为5.71°，即等效的摩擦系数为0.1。本例仅考虑接触作用，不考虑界面粘结与断裂。作为目标界面的土层上表面的接触单元网格尺寸为0.4 m，作为从界面的基础下表面的网格为0.2 m，对应的临界步长约为0.2 ms。使用有限质点法对本例进行计算时，计算步长取为0.1 ms，计算总时长1 s，均布荷载缓慢斜坡加载至基础上表面。为了加速收敛速度，取虚拟的全局质量阻尼系数100。作为对比，本例使用通用有限元软件的罚函数法接触模型，接触刚度取558 MPa，其余计算参数均保持一致。为了验证界面单元在不同网格形状上的适应性，在FPM平台中分别使用四面体实体单元和六面体实体单元进行建模分析。

使用FPM平台计算所得的结构位移云图如图10所示，接触面所在平面与Y=0平面交界线位置的界面单元的法向接触应力分布如图11(a)和图11(b)所示，并与Qian等[34]的结果进行了对比。可以看到，本文的有限质点法界面单元的模拟结果呈现出了弹性力学解中的中心应力均匀两端应力集中的状态，且相较于有限元结果，本文中经由法向量插值的曲面连续化处理使得法向接触应力的分布更加均匀。另外，Qian等[34]在模拟中使用了接触面的曲面光滑算法RPIM，目的是让接触面上的接触力分布更为平滑连续。而本文采用了基于形函数插值的法向量连续化处理，也同样产生了一定的应力平滑效果，虽然平滑效果有限，但计算代价极低。本例说明了本文实现的有限质点法界面单元在处理接触问题上的精度和有效性。

图10 基础-土体接触：四面体网格模型变形云图Fig. 10 Interface contact analysis between foundation and soil: displacement contour of the model with tetrahedron mesh

在四面体模型中，实体单元表面为三角形，用一个界面单元进行覆盖；六面体单元表面为四边形，使用两个三角形界面单元进行覆盖。两种网格划分的模型使用FPM平台计算所得的中心面接触应力分布如图11(c)所示。除了四面体网格由于网格走向不对称所造成的应力分布不完全对称的情况之外，两种界面单元覆盖方式下的接触应力分布具有极高的吻合度，证明了本文的三角形界面单元具有普遍的几何适应性，能够在建模时以覆盖在其他三维单元表面上的形式赋予其表面界面作用分析的能力。

图11 基础-土体接触：Y=0平面界面单元接触应力分布Fig. 11 Interface contact analysis between foundation and soil: contact stress distribution on the Y=0 plane

2.2 砖块粘合界面剪切

本例对采用文献[22]所使用的Beyer等[35]进行的砖块粘合界面剪切试验进行数值模拟，用于测试本文界面单元在处理界面粘结后剪切开裂并耦合接触时的模拟效果。原始试验装置如图12(a)所示，三块砖块之间用砂浆粘结并列放置，两侧使用加载装置进行水平向的施压。中间砖块顶端靠近界面的位置进行了竖向的位移约束，两侧砖块在靠近接触面的位置作用了竖向的移动支座，用于对界面施加剪切作用。

根据此试验装置，本文建立了简化的有限质点法数值模型，相关尺寸和边界条件如图12(b)所示。数值模型采用对称的一半建立，砖块使用四面体实体单元进行建模，中央半块砖块的顶部和左侧施加对应法向的固定约束，右侧砖块底部施加竖向的剪切位移约束，速度固定为25 mm/s，右侧作用法向压力P，分别取值为0.2 MPa、0.4 MPa和0.65 MPa。界面左右位置的两个面覆盖了界面单元，界面单元尺寸在3 mm～4 mm，砖块的实体网格在界面位置与界面单元网格重合，边界端网格尺寸适当放大，如图13所示，模型的临界步长约为0.38 µs。

图12 砖块粘合界面剪切：试验装置图与简化的数值模型Fig. 12 Shear test of masonry wallettes: experiment setup and simplified numerical model

图13 砖块粘合界面剪切：FPM平台中网格划分Fig. 13 Shear test of masonry wallettes: the FPM mesh

根据材性试验，砖块弹性模量为14 GPa，泊松比为0.15，密度为940 kg/m3。界面断裂释放能为750 J/m2，界面粘结力的临界值取为0.25 MPa，界面摩擦系数为0.77。根据文献[24]所建议的接触刚度的计算公式，可以得到近似的接触刚度2590 GPa。考虑到接触位置过刚会导致临界时间步长显著减小，实际计算中减小两个量级。计算中步长保持恒定为0.05 µs，通过保持位移加载速度的25 mm/s，模拟支座位移加大到10 mm的过程。在FPM平台中计算时，需要设置两个分析过程，首先仅对模型右侧作用法向压力P；待稳定后再施加底部的位移荷载。由于本例为轴压状态下的剪切断裂，求解过程中设置接触面法向始终保持粘结，允许切向剪切断裂。图14展示了不同条件下FPM平台所计算得出的界面剪切应力-剪切变形曲线，其中应力和变形均提取自作用位移荷载的第二个分析过程，应力值取界面中心0.15 m高度范围内界面单元质点应力的平均值，此过程界面法向压力保持恒定。

图14 砖块粘合界面剪切：界面剪切应力-剪切位移曲线Fig. 14 Shear test of masonry wallettes: shear stress plotted against shear displacement

图14(a)～图14(c)展示了当法向压应力分别保持为0.2 MPa、0.4 MPa和0.65 MPa时，随着剪切位移增大，界面剪切应力的变化曲线，并与Beyer等[35]的原始试验数据进行了比较。有限质点法模拟结果与试验结果呈现了相同的变化趋势。从模拟结果可以看出，在法向压应力保持不变的情况下，界面剪切变形从0开始增长，界面间始终存在摩擦并很快达到最大摩擦力，保持滑动，切向内聚力开始发挥作用；随着剪切变形继续增长，切向的内聚力逐渐减小，直至内聚力减小为零，界面间仅存在接触造成的滑动摩擦的作用，剪切应力大小逐渐稳定。因此，本文所实现的界面单元在受压剪切的情况下的本构曲线形式实质上是外加内聚力曲线和切向摩擦力曲线这二者的叠加。对于模拟中限定的断裂释放能以及临界应力值，可以计算得到界面内聚力在剪切位移为6 mm左右的位置失效，也与模拟结果相吻合。试验结果曲线更为平滑，这主要是由于实际情况中摩擦力和界面内聚力的耦合更加连续；而且本文采用了线性外加内聚力模型，较真实情况有一定的简化。

诸多文献同样使用此试验进行了数值模拟，这里选择了Snoozi和Molinari[36]、Spring和Paulino[24]以及Baek和Park[22]三者的模拟结果与本文的结果进行比较，如图14(d)所示。Snoozi和Molinari[36]采用了与本文相同的线性外加内聚力模型，但Snoozi和Molinari[36]采用了指数形式的摩擦力规则化函数以引入界面开裂损伤对摩擦力的影响。Spring和Paulino[24]以及Baek和Park[22]使用了基于势能的Park-Paulino-Roesler内聚力模型的变化形式，更符合真实物理过程，但计算相对复杂。需要说明的是，为了匹配试验结果，Snoozi和Molinari[36]、Spring和Paulino[24]并未采取试验得出的界面粘聚力临界值0.25 MPa，Snoozi和Molinari[36]取0.3 MPa，而Spring和Paulino[24]取0.45 MPa。

图14(e)～图14(f)进一步展示了剪切应力-剪切位移曲线随着界面临界应力 σcr和断裂能Ec取值不同而变化的情况。当界面临界应力保持不变时，断裂能越大，界面开裂过程的内聚力作用距离越长，界面完全断裂所需要变形也就越大；当断裂能保持不变时，界面临界应力越大，剪切应力峰值越大，界面完全开裂所需要的剪切变形就越小，越快进入完全接触摩擦的平台段。

总地来说，本文基于有限质点法所开发的界面单元在界面粘结-断裂-接触耦合的复杂行为时具有足够的精度以及有效性。

3 结论

基于有限质点法基本理论并依托有限质点法通用计算平台FPM，本文开发了受用于并行化的有限质点法计算体系下的三角形界面单元，用于模拟界面之间的粘结、断裂、接触等复杂的界面作用行为。该单元采用并行化的作用对侦察算法，实现了基于等加速度原理的粘结作用、基于内聚力模型的断裂作用，以及基于罚函数法的接触作用这三种作用下的状态判据和作用力计算。数值算例证明本文所开发的界面单元在处理普通接触问题以及粘结-断裂-接触耦合的复杂界面作用行为方面是正确并有效的。利用本文提出的界面单元，可对岩土和复合结构等领域的实际工程常见的复杂界面行为进行较为高效的模拟，也为实际工程中复杂界面行为的精细化大规模数值分析奠定了基础。