基于块体离散元法的多面体接触重叠算法

刘广煜， 徐文杰， 周 乾， 冯泽康

(清华大学 水沙科学与水利水电工程国家重点实验室，北京 100084)

1 引 言

自然界中的颗粒材料，如沙子和砾石，是由几何上带棱角的非球形颗粒组成的，而颗粒的几何形状对于颗粒系统的物理力学特性有非常重要的影响。近年来，非连续介质力学方法的发展为颗粒材料的研究提供了更多手段。在这些方法中，最常用的是离散元法，由Cundall等[1]首次提出，在海洋工程、岩土工程及化学工程等领域有广泛的应用[2-5]。

根据颗粒的几何形状，离散元法可分为球体离散元法和块体离散元法。球体离散元法由于接触检测和接触力计算简单，得到广泛的发展和应用。为更加真实地体现颗粒几何形状的影响，球体离散元法在颗粒碰撞时引入了弯矩和扭矩[6]。但是弯矩和扭矩的计算需要引入额外的细观参数，这使参数反演过程更加复杂[7]。球体颗粒簇模型是离散元法中另一种研究颗粒几何形状的有效方法[8]，但是该方法的计算成本较高。块体离散元法解决了球体颗粒难以表征颗粒形状的问题，能够模拟颗粒的复杂形状，且无需引入复杂的细观参数[9,10]。

在离散元法中，颗粒假定为刚体，并允许颗粒之间有微小的重叠。基于重叠部分计算嵌入深度、法线方向、接触面积和接触点等接触特性，可求得接触力。球体离散元法计算碰撞球体间接触特性相对简单，而块体离散元法计算碰撞块体间接触特性难度大且耗时。公共平面法是目前应用较为广泛的一种块体接触算法[10,11]，该方法对接触关系进行分类，如点-点、点-棱和点-面等，分别求解接触特性及接触力。公共平面法为了确保数值模拟的稳定性，会针对某些特殊接触类型(如棱-棱接触)进行特殊处理，导致接触计算耗时较长。众所周知，多面体在顶点处存在法线方向不连续的问题[9]，因而两碰撞多面体在顶点处发生小的相对运动时，公共平面法不能保证接触法线方向的平滑演化，可能导致在数值模拟中引入额外的能量。为解决这类问题，冯云田等[9,12]提出了能量守恒接触模型，证明了碰撞的多面体颗粒利用重叠多面体计算接触力是符合能量守恒定律的。重叠多面体的计算目前是一个难点。Yade作为一种开源和扩展性强的离散元法框架，基于CGAL算法库计算重叠多面体[13]，但法线方向是通过最小二乘拟合得到，因而该方法耗时且模拟的块体颗粒数有限。为提高计算效率，基于图像处理器GPU(Graphic processing unit)，Govender等[14]发展了并行离散元计算软件 Blaze -DEM，该软件显式计算得到重叠多面体，能够实现百万级多面体颗粒模拟计算。但由于接触计算精度较低且接触模型少，难以在岩土工程得到应用。

本文在耦合模拟器的块体离散元框架(CoSim-DEM)下开发了多面体颗粒的接触重叠算法，CoSim-DEM是基于GPU并行计算的面向对象的可扩展性强的数值模拟仿真平台。接触重叠算法是基于几何对偶理论[15]衍生的多面体接触理论，是一种更精确的接触算法。基于重叠多面体，可以计算接触点、接触面积、法线方向和嵌入深度等接触特性，进而用于接触力和力矩求解计算。

2 多面体接触重叠算法

在块体离散元中，颗粒的迭代计算主要包括邻居搜索、接触检测、接触力计算和插值计算四个步骤。其中，邻居搜索就是通过均匀网格法[16]或边界包围盒层次树法[17]等找到每个颗粒的邻近颗粒列表。接触检测基于邻居搜索得到的列表进行重叠多面体计算，进而求解接触力。多面体接触重叠计算算法流程如图1所示。

图1 多面体重叠接触算法流程

本文仅考虑凸多面体的接触算法，因为对于任意形状的复杂块体，可以由若干个凸多面体拼接而成，因此凸多面体的算法同样适用于复杂块体的碰撞。基于几何对偶理论衍生出来的多面体接触理论，将位于初始坐标系的两多面体转换到对偶空间中，此时两多面体的所有面转换成了对偶空间中的点，在对偶空间中计算所有点的凸包，再将该凸包转换到初始坐标系下生成一个新的多面体，此时初始坐标系下的这个多面体即为两碰撞多面体的重叠体。多面体接触重叠算法基于该理论，融合了GJK算法和快速凸包算法，基于重叠体计算接触面积、接触法线方向、嵌入深度和接触点等接触特性，可以提高计算效率，保证计算结果精确且高效。该算法主要分为GJK算法、几何对偶算法、快速凸包算法和接触特性计算四步。

2.1 GJK算法

对于两个块体A和B，其接触关系在欧氏空间中可以分为碰撞、触碰和分离三种类型(图2)。

图2 两多面体间接触状态

GJK算法计算成本低且通用性高，通过迭代，能够检查两个块体是否发生碰撞。该算法利用多面体A和B各自的顶点点集a和b，通过寻找Minkowski差分A-B(A-B={a-b,a∈A,b∈B})内部或外部边界上的点V，找到距离原点O最小的点，如图2所示。在确定两块体的接触状态为碰撞或触碰时，执行仿射变换将点V转换为v′(v′∈A，v′∈B)，该点将用于后续对偶空间变换算法。

2.2 对偶空间变换算法

GJK算法筛除了分离状态的多面体对，而对于存在接触的多面体对，需要对其进行几何对偶化，将其转化到对偶空间中，用于后续凸包的计算。几何对偶化具体操作是将初始坐标系的多面体的面变换到对偶空间的点。对偶空间的点vdual,j可表示为

d=(ci-v)·ni，vdual,j=ni/d

(1，2)

式中ci和ni分别为多面体每个面的面心和法线方向，v为重叠体内部的任意一点。GJK算法计算出的点v′不能直接作为重叠多面体的内部点，因为v′有可能在重叠体的面上(图3)。因此，Koziara[18]提出了向重叠多面体内部延伸点v′的迭代算法，如图3所示，最终延伸方向nt可以是若干个多面体的面法线反方向(图3的n1，n2和n3)的组合，进而得到重叠体内部点v。延伸的距离由自定的阈值决定，可根据精度要求进行设置。根据式(1，2)，基于重叠体内部的该点v，多面体A和B的面将转换为对偶空间中的散点vdual,j。

图3 重叠体内部点v计算算法

2.3 快速凸包算法

在对偶空间中，通过快速凸包算法，搭建包含所有散点vdual,j的最小凸多面体。图4显示了基于给定点集构建最终凸包的整个迭代过程。如图4(a)所示，创建一个初始四面体作为初始凸包，给定散点集由此分成两组，第一组是初始四面体四个不共面的顶点(图4(a)的VT k)；另一组是位于四面体面外(图4(a)的fk)的点(图4(a)的VF i,j)。初始四面体搭建完成后，需要迭代计算建立整个凸包。

为说明迭代算法，迭代计算凸包的其中一个迭代步如图4(b)所示。VF 1,1是f1外部且距离f1最远的点。算法中，f1将标记且删除，因为该面不是凸包的边界。遍历f1的所有相邻面(图4(b)的f2，f3和f4)，满足式(3)的面也会标记删除。然后基于点VF 1,1和已删除面(f1)的边，创建新的面fN 1，fN 2和fN 3。此时更新后的凸包由原始的面(图4(b)的f2，f3和f4)和新创建的面(图4(b)的fN 1，fN 2和fN 3)构成。最终，当凸包外部没有点时，迭代终止，此时整体凸包构造完成(图4(d))。

图4 快速凸包算法

nf·(VF i,j-cf)>0

(3)

式中cf和nf分别为凸包表面的面心和法线方向。

2.4 接触特性计算

图5 两碰撞多面体的重叠体

(4)

(5)

kn=(EA+EB)·a/(RA+RB)

(6)

(7)

ks=(EA·vA+EB·vB)·a/(RA+RB)

(8)

3 多面体接触重叠算法验证

为验证接触重叠算法的适用性和准确性，模拟分析了两个算例，分别是颗粒碰撞的基准测试和砌体结构破坏过程，对比了模拟结果与试验结果。

3.1 块体碰撞测试

构造了三种块体碰撞类型的基准测试，分别是点-点、点-面和棱-棱碰撞，如图6所示。在这些基准测试中，具有相同质量的颗粒在重力作用下从相同高度(1 m)落下，并与下方固定的颗粒相撞，此时回弹系数设置为0.1。此外，颗粒的几何形状均采用轴对称多面体，包括八面体、六面体和三棱柱。为研究回弹系数的影响，针对点-面碰撞类型采用了三个回弹系数进行模拟，并对比了不同回弹系数下碰撞力的演化过程。表1是模拟中使用的计算参数。

图6 颗粒碰撞基准测试

表1 颗粒碰撞的数值计算参数

图7展示了三种不同接触类型在颗粒碰撞时接触力的演化。如图7(a,b)所示，在点-点接触的情况下，尽管嵌入深度已经大于八面体颗粒尺寸的15%，接触力仅作用在垂直方向(Z方向)。在三种碰撞类型下(点-点、点-面和棱-棱接触)，X和Y方向的接触力都相同且都等于0，且由于黏滞力消耗能量，多次碰撞后，上方颗粒均可以静止于下方颗粒上方。由于颗粒形状均为轴对称图形，碰撞力在X和Y方向分量为0是合理的。进而对比了三种碰撞类型在垂直方向的接触力的演化，如图7(c)所示，由于颗粒碰撞中黏壶的作用，部分能量得到释放，最终接触力可以收敛到上部颗粒的重力。不同碰撞类型接触力收敛速度不同，由式(6)可发现，接触面积越大会导致接触刚度越大，点-面接触由于接触面积最大，因此收敛速度最快。回弹系数对接触力的收敛速度有更为显著的影响，较小的回弹系数可以明显加快接触力收敛速度，如图7(d)所示。在回弹系数为1.0时，发生多次碰撞后，上部颗粒都能够回到初始位置，每次碰撞力的演化均一致，由于颗粒系统中无摩擦和黏滞力导致的能量损失，因此碰撞过程符合能量守恒，验证了碰撞过程中没有额外的能量增加。该算例充分验证了多面体接触重叠算法的适用性，可广泛用于多体接触模拟计算。

图7 碰撞测试下不同接触类型的接触力和颗粒位置的演化

3.2 砌体结构破坏过程

砌体结构在结构工程中十分常见。为进一步验证算法在计算颗粒系统的旋转、位移和接触力等方面的准确性，采用了文献[19]的砌体结构破坏试验进行模拟对比。

Portioli等[19]使用的模型是干接缝凝块岩块制成的砌体面板(图8)。面板由两种砌块组成，尺寸分别为100 mm×200 mm×50 mm和100 mm×100 mm×50 mm。面板的尺寸是高400 mm、宽1100 mm和厚100 mm。底部左侧设置长度为 600 mm 的固定支撑台，支撑台上方的砖块仅可以沿水平方向移动。底部右侧设置长度为500 mm的升降台，可在水平和垂直方向自由移动。为监测升降台上的承载力和位移，设置机械千分表和台秤分别测量升降台的垂直位移和承载力。试验过程是升降台每次缓慢降低1 mm，待面板运动稳定后，记录升降台的位移变化，同时记录升降台的承载力变化。根据模型试验，相应的数值模型如图8(c)所示，块体之间没有设置黏聚力，其中固定支撑结构和升降台分别由固定块体和可移动块体表示。在模拟过程中，可移动块体以10-4m/s的速度向下移动。可移动块体与上方块体之间的摩擦角设置为零，以确保在水平方向上有自由位移，块体颗粒之间的摩擦力设置为20°，计算过程中监测可移动块体受到的外力变化，进而与试验结果对比。其他参数列入表2。

图8 砌体结构模型试验及离散元模型[19]

表2 砌体结构的数值计算参数

数值结果和模型试验结果在面板破坏时的比较如图9所示。可以看出，数值结果与模型试验结果吻合较好。根据变形特征，砌体结构中部发生滑动和转动破坏(图9(b)的第2部分)，表明该算法可以准确模拟块体的旋转和平移。左侧部分(图9(b)的第1部分)表现得像不可变形体，位移非常小；右侧部分(图9(b)的第3部分)也类似不可变形体，但通过阶梯状裂缝与中间部分分开。图10显示了作用在升降台的承载力随垂直位移的演化。从数值结果与试验结果的对比来看，两者吻合较好。而且，砌体结构发生破坏时数值模拟计算得到的承载力约为104 N，这也与试验结果相一致。该算例充分验证了算法的精确性。

图9 可移动块体竖直位移20 mm的模型试验与模拟结果对比[19]

图10 作用在升降台的承载力随位移的演化

4 结 论

离散元法已广泛用于颗粒材料的物理和力学研究。在块体离散元法的框架CoSim-DEM下，开发了多面体接触重叠算法，该算法能够显式计算得到精确的重叠多面体。根据计算流程依次分为四个子步骤，分别是GJK算法、对偶空间变换算法、快速凸包算法和接触特性计算。基于精确的重叠多面体，可以直接计算法线方向、嵌入深度、接触点和接触面积等接触特征，更合理地计算接触力。

采用了颗粒碰撞的基准测试验证该算法的适用性和准确性。通过计算结果可以发现，该算法能够准确计算各种接触类型的接触力，如点-点接触、点-面接触和棱-棱接触，并发现回弹系数越小，接触力收敛越快。进而对砌体结构的破坏过程进行了模拟，并与模型试验进行了对比，表明数值模拟对破坏机理的揭示和再现承载力的演化非常有效。本文开发的多面体接触重叠算法避免了多面体之间接触搜索的繁琐几何判断，其计算结果具有很好的可靠性，可应用于非规则颗粒材料的离散元模拟。