多面体颗粒的接触识别及离散元动力学建模

洪 俊 李建兴 沈 月 王 潇

(1东南大学土木工程学院， 南京 211189)(2东南大学江苏省工程力学分析重点实验室， 南京 211189)

对于多面体形状的颗粒，接触识别算法是数值模拟中最复杂的问题.为避免复杂的接触识别，目前模拟软件基本采用球体代替多面体，或者采用大量球体组合来拟合多面体.球体在模拟时的计算效率很高，但是可能导致某些不真实的情况出现.在耗散型颗粒系统中，例如沙堆，球形比多面体颗粒更容易滚动，模拟中会出现不正常的静安息角和非常低的摩擦系数[8].球体组合可以较为灵活地拟合不规则颗粒的形状，但会引起颗粒表面粗糙化和波纹化，同时，精确模拟单个颗粒需要大量的球单元，导致计算效率较低.Cleary[9]指出，多面体的尖角问题是球体组合模拟最为困难的事情.因此，如何直接使用多面体颗粒进行建模是增强模拟精度的关键.

随着计算机图形学的进步，多面体接触识别算法也得到了迅速发展，并部分应用到离散元模拟中.Cundall[10]提出的公共面法引入了多面体之间的公共平面，通过测试2个多面体与该平面的关系来判断它们的接触类型.尽管研究人员对这种方法进行了很多改进[11-12]，但在确定公共面后，仍需要对接触类型进行判别，并确定接触点.GJK算法[13]目前应用较多，其优点是在计算少量颗粒时效率较高，但是在2个多面体几何体尺寸相差较大的情况下，失效概率较大，并且只能对凸多面体的接触类型进行判断.射线穿透法[14]的原理较为简单，能够判断非凸体颗粒的相互接触，而不用将非凸体分解为凸体，适用于对几何形状较为复杂的多面体颗粒进行建模.

本文首先将多面体颗粒接触识别分成粗检和细检2个过程，阐述了连接线性目录法，并对射线穿透法进行了改进；然后，基于离散元法基本原理，建立了颗粒间的接触力模型；最后，通过自编程序对规则凸多面体和复杂不规则多面体颗粒系统自由堆积过程进行了数值模拟.

1 接触识别算法

多面体颗粒间的接触识别算法较为复杂.颗粒系统包含大量颗粒，如果采用两两颗粒判断接触的方法，计算量将非常庞大.为了加速计算，先采用邻居目录算法粗检出可能存在的接触颗粒对，然后对所有可能接触的颗粒对进行接触与否的细检验.

1.1 邻居目录算法

邻居目录算法的目的是粗检出可能存在接触的颗粒对.在邻居目录算法中，较为有效的算法有Verlet邻居目录法(VL)、连接单元法(VC)和连接线性目录法(LLL).文献[15]的研究结果表明，在计算大量多面体颗粒碰撞时，连接线性目录法的计算效率远高于Verlet邻居目录法和连接单元法.本文采用连接线性目录法建立颗粒系统的邻居目录.

连接线性目录法的本质是一种边界盒相交检测算法.以二维情况为例，其基本步骤如下：

① 围绕每个颗粒设置边界盒(见图1)，使边界盒恰好容纳颗粒，并且边界盒边缘与系统的全局坐标轴保持平行.

② 将边界盒起点和终点分别投影到x，y坐标轴上.颗粒i(i=1，2,3，4)的边界盒起点bi和终点ei在x，y轴上投影分别见图1.将bi和ei按坐标轴上的顺序存储在线性目录中，生成2个列表.

图1 边界盒示意图

每个线性目录列表都包含2n个元素，对应于系统中颗粒数量的2倍.如果某个颗粒的起点和终点间存在另一个颗粒的起点(如图1中3/4颗粒在x轴上的投影)、终点(如图1中1/4颗粒在x轴上的投影)或起点和终点(如图1中2/3颗粒在x轴上的投影)，那么这2个颗粒的边界盒在该轴的投影存在重叠.当2个边界盒在x，y坐标轴上都存在重叠时(图1中3/4颗粒对)，说明边界盒相交，2个颗粒互为邻居.使用稀疏矩阵的形式创建邻居目录链表，链表中存储了互为邻居的3/4颗粒对(见图2).

图2 邻居目录链表

在每个时间步长内，通过检查边界盒投影的线性目录列表来新建邻居目录链表需要耗费大量时间，可以采用分类更新算法来提高效率.原因在于，时间步长小，线性目录列表的顺序变化不大，故只需要对旧的邻居目录链表进行更新，而不必新建链表.更新邻居目录链表时，只需考虑以下4种情况：

1) 2个起点发生交换，表明2个边界盒在该坐标轴的投影已经重叠并且继续重叠.

2) 2个终点发生交换，表明2个边界盒在该坐标轴的投影已经重叠并且继续重叠.

3) 起点和与之最接近的终点发生交换，表明2个边界盒在该坐标轴的原本重叠的投影发生分离，须删除该重叠.

4) 终点与另一个边界盒的起点发生交换，表明出现了一个新的重叠.

对于前2种情况，只需要在线性目录列表中交换顺序，不需更新邻居目录链表，因为颗粒的邻居状态没有改变.对于第3种情况，如果2个边界盒在任意一个坐标轴上的重叠被删除，则边界盒不相交，解除这2个颗粒的邻居关系.第4种情况则需要添加新的重叠，并检查两边界盒在其他轴上的重叠情况.

邻居目录建立后，细检验只需测试目录中的颗粒对是否发生碰撞，而不必遍历所有颗粒对进行测试.因此，接触识别的计算效率将显著提升.

1.2 射线穿透法

使用射线穿透法对邻居目录中可能存在接触的颗粒对进行接触与否的细检验.射线穿透法的基本原理是：对于空间内的一点，生成任意方向的射线.当射线与多面体的边界或面的非相切交点个数为奇数时，该点位于多面体内；交点个数为偶数时，该点位于多面体外.如图3所示，由点P生成的任意方向的射线与多面体的交点个数为奇数，则点P位于多面体内；由点Q生成的任意方向的射线与多面体的交点个数为偶数，则点Q位于多面体外.由于有效检测的需要，通常会在点P和点Q多次生成射线.

图3 射线穿透法基本原理

接触识别步骤如下：

① 排除被测试多面体上不在另一个多面体边界盒内的所有点.如图4所示，被穿透多面体1相对待检测多面体2，只需检测边界盒内的7个点.本文在射线穿透法充分利用连接线性目录法中的边界盒，精简了编程代码，提高了计算效率.

② 对于边界盒内的检测点，需对多面体1的所有面生成射线并进行测试.具体的射线穿透算法程序编制过程参考文献[14].

图4 2个多面体相交的平面示意图

射线穿透法本质是点-面接触识别算法，下面对多面体接触特殊情况(点-点、点-边、边-边、边-面和面-面接触)分别进行讨论.

对于点-点和点-边接触，考虑变量存在精度以及接触后外法线方向无法确定问题，采用延迟识别法侦测下一时间步中的微小穿透现象，将被检测的点或边所在多面体上距检测点最近的表面作为接触面，则点-点和点-边接触转化为点-面接触问题.

射线穿透法中边-边接触问题无法处理.如图5所示，采用在多面体边上附加检测点的方法，对原算法进行改进.通过判断附加检测点侦测是否和被检测多面体发生接触，将边-边接触转化为点-面接触问题.改进算法解决了射线穿透法无法处理边-边接触的问题，但增加了计算量，同时附加检测点数目和位置也会影响计算结果的精度.

图5 附加检测点示意图

对于边-面和面-面接触问题，利用边上原有顶点和附加检测点将边-面和面-面接触问题转化为点-面接触问题.

2 离散元法

离散元法的基本思想是将非连续体分离为刚性单元的集合，使各个刚性单元满足运动方程，继而求得非连续体的整体运动形态.离散元法允许单元间的相对运动，不一定满足位移连续和变形协调条件，计算速度快，所需存储空间小，尤其适合求解大位移和非线性问题[16].离散元法的运动方程及积分方法参考文献[17]，本节仅介绍接触力计算过程.

在使用多面体离散元模拟颗粒系统时，多面体颗粒有效接触识别过程非常复杂，颗粒间接触力的计算也较为困难.目前，多面体颗粒间的接触力定义较多，结合射线穿透法的点-面接触特性，接触力可由以下公式确定：

(1)

(2)

式中，rA为穿透多面体颗粒A上检测点的位置；rB为A点在被穿透多面体表面的最近投影位置.

(3)

式中，vA(P)，vB(P)分别为多面体颗粒A，B在接触点P处的速度矢量.

(a) 凸体 (b) 非凸体

3 程序实现与算例

基于离散元法，使用C语言进行程序编制.程序包括主程序和后处理两大模块.主程序包括数据输入模块、循环计算模块和数据输出模块.数据输入模块读入颗粒系统的物理参数和几何参数；循环计算模块采用连接线性目录法和射线穿透法，对多面体颗粒的接触进行检测，并采用Verlet积分法[17]模拟动力学过程；数据输出模块记录每个时间步长内颗粒的物理信息，并输出到文件.后处理程序使用OpenGL软件处理图像，包括数据输入模块和图形演示模块.

接触识别算法和模型的有效性可以采用下面2个算例进行验证.

算例1为规则凸多面体颗粒系统自由堆积过程模拟.在一个底部半径为0.4 m、高为1.2 m的圆柱筒内，48个规则凸多面体颗粒在重力作用下做自由落体运动.颗粒尺寸是底部外切圆半径为0.057 m、高为0.114 m的六棱柱.随机选取颗粒的初始方位角，圆柱筒设置为透明.各计算参数如下：密度ρ=1 500 kg/m3；阻尼系数dAB=0.07；接触刚度kAB=10 MN/m；时间步长dt=50 ns.图7给出了t=0,0.06,0.12 s时的规则凸多面体颗粒系统整体位置.

(a) t=0 s

算例2为复杂不规则多面体颗粒系统自由堆积过程模拟(同时含有凸和凹多面体).在一个0.2 m×0.2 m×1.5 m的长方体中，90个复杂不规则多面体颗粒在重力作用下做自由落体运动.为便于观测，长方体用轮廓线代替.各计算参数如下：密度ρ=1 000 kg/m3；阻尼系数dAB=0.05；接触刚度kAB=5 MN/m；时间步长dt=50 ns.图8给出了t=0,0.1,0.3 s时的复杂不规则多面体颗粒系统整体位置.

(a) t=0 s

2个多面体颗粒系统自由堆积算例在模拟过程中均没有出现颗粒间以及颗粒和边界间的过量穿透，也未出现颗粒间碰撞时颗粒异常运动，说明连接线性目录法能够快速建立颗粒系统的邻居目录，改进的射线穿透法能有效地识别颗粒间各种复杂的接触情况，动力学模型能够真实地反映颗粒系统运动状态.

4 结论

1) 采用连接线性目录法建立颗粒系统的邻居目录，使用改进的射线穿透法判断邻居目录中颗粒对是否真正接触.结合2种不同阶段的接触算法，共用边界盒，精简了编程代码并提高了计算效率.

2) 对多面体接触类型进行了讨论，通过延迟识别法和附加检测点的方式，解决了包括边-边接触、边-面接触等几种极端接触类型的接触识别问题.

3) 通过规则凸多面体和复杂不规则多面体的颗粒系统自由堆积算例，证明连接线性目录法能够快速建立颗粒系统的邻居目录，改进的射线穿透法能有效地识别颗粒间各种复杂的接触情况，建立的动力学模型能够真实地反映颗粒系统运动状态.