基于统一标定的势接触力计算

严成增，郑 宏，葛修润

（中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071）

1 引 言

岩体中存在节理、断层等大量不连续面，对这类问题，基于连续介质力学的数值类方法，如有限元通常难于处理。为了处理这类非连续介质力学问题，诞生了离散元类数值方法。其中最有代表性的为Cundall[1]提出离散元方法和石根华[2]提出非连续变形分析方法（DDA）。这两种方法在计算接触力方面的共同点是均引入了刚度参数，通过嵌入量计算接触力，不同的是Cundall提出的离散元允许块体间相互嵌入，而石根华[2]提出DDA则通过开闭迭代，使得块体间不嵌入。邬爱清等[3]对DDA中接触力的计算过程予以详细的介绍，指出通过目前的接触状态，按总体方程求解，当系统所有接触部位满足不嵌入无张拉准则时，此时计算出的接触力分量则为真实接触力分量，否则通过增删弹簧继续迭代求解。郑宏等[4]为避免引入刚度参数和开-闭迭代，将互补理论引入到DDA中，将DDA中各潜在接触对以互补形式来表现，通过求解非线性互补方程组来获得问题的解答，该方法与传统的接触力计算方法有很大的不同。徐栋栋等[5]将Munjiza基于势的接触力计算方法引入到流形元中，用于处理流形元的接触问题。以上计算接触力的方法均存在一些问题，例如Cundall的离散元方法，不能处理角-角接触的问题，在其开发的二维版本的离散元程序UDEC中通过圆角化来规避这一问题[6]。在三维情形下，虽提出了公共面法来求解接触力，但仍存在公共面不惟一[7]，以及无法处理某些极端接触问题[8]。石根华的DDA也不能很好地处理角-角接触的问题，同时当问题扩展到三维状态时，接触判断变得相当复杂，这也是三维DDA发展较为缓慢的一个重要原因。

与以上各种计算接触力的方法不同，Munjiza在其提出的有限元-离散元（FDEM）耦合分析方法中，采用了一种非常新颖的基于势的接触力计算方法[9－10]。有关该方法的相关发展可参见文献[11－12]。该方法计算的接触力为分布力，不会产生应力集中，对点-点接触无需作特殊处理，而常规的离散元(DEM)为了处理点-点接触及避免角点处产生应力集中需要对角点进行圆角化。基于势的接触力计算法虽然规避了传统DEM在处理点-点接触上存在的困难，但仍存在势的物理意义不够明确，相同的嵌入量计算的接触力大小不一致等问题。在已有研究的基础上，本文提出了一种基于统一标定的势接触力计算方法，该方法重新定义了势函数，使得势的物理意义更为明确，该方法中势就是对嵌入的一种表征；同时，在相同的嵌入量时，计算的接触力大小是一致的。在解决了原有势接触力计算方法物理意义不够明确，对于相同嵌入量，其接触力的大小不相同等问题的同时，保留了原有Munjiza势接触力方法的所有优点，即接触力为分布力，对点-点接触无需特殊处理，接触力所做的总功为0，系统的能量不会因为接触力的存在而增加或者减少。

2 FDEM势接触力

如图1所示，将接触块体分别称之为目标块体cβ 和靶块体βt，目标块体嵌入到靶块体的面积dA所引起的接触力为[10]

式中：Pt为重叠区域内位于靶块体 βt内的一点；φt(Pt)分别为 Pt点在靶块体 βt内的势，grad为梯度；pn为法向罚参数。

图1 接触力计算示意图Fig.1 Potential based contact force between two elements

同理，靶块体βt嵌入到目标块体βc中的面积dA 引起的接触力为

特别说明，式（2）计算的接触力为目标块体βc提供势对靶块体βt的作用力，于是靶块体βt对目标块体βc的反作用力大小与式（2）相同，方向相反，即为

于是，重叠区域dA 引起的总接触力（靶块体对目标块体）为[10]

于是，整个重叠区域所引起总的接触力可通过对上式积分得到

对式（5）应用格林公式，可将面积分转化为对重叠区域边界的线积分：

Munjiza等[9－10]指出，势函数φ 在离散的单元边界上应为常量，只有这样才能保证接触力做的总功为0，这也是基于势的接触力对势函数定义的惟一要求。据此，Munjiza[10]定义了如下势函数：

式中：P为靶单元中长度为L 的边上的任意一点。

如图2所示，Ai(i=1，2，3)为P 点与三角形3条边所构成的子三角形面积，A为三角形的面积。这样在三角形中心点的势等于1，而在边界上的势等于0。为了便于理解，以三角形的中心点C 将三角形分为3个子三角形，要据式（7）求解势，只需先判断点P 在那个子三角形内，例如点P 在C12子三角形内，则该点的势据式（7）可进一步改写为

式中：hP-12表示点P 到12边的距离；h0-12表示0点到12边的距离，亦即三角形12边的高。

图2 三角形单元内一点P 的势Fig.2 Potential at point P in a triangle element

这种势函数的定义，虽然满足了在边界上势为常量的要求，但存在物理意义不够明确，相同的嵌入量计算出的接触力大小不等的问题。

对该势函数稍作分析，便可发现其中的问题，如图3所示，左右两组接触单元，除了靶单元的大小不同外，两组接触单元的嵌入量完全相同，其中目标单元标号为L 的边与靶单元的12、1′ 2 ′边平行，且该边到12、1′ 2 ′边的距离相等，亦即嵌入量相同，L 边的外法线向量为 nL。以图3靶单元中长度为L的边为例，据式（6）、（8）对该边积分得到的势接触力，其中图3中右边的长度为L 的边积分得到的势接触力为

式中：pn为法向罚参数。

同理，图3中左边的长度为L 的边积分得到的

势接触力为

图3 嵌入量相同的两组接触Fig.3 Embedding the same amount of two group contacts

此处仅仅是对重叠区域边界的一条边作了分析，总的接触力需要对整个重叠区域的边界进行积分，与上述分析类似，对图3所示的情况，原有程序计算出目标单元受到的总接触力如表1所示（均已分配到单元的节点上），采用的罚参数pn=67.3 GPa。表1表明，在嵌入量完全相同的情况下，目标单元所受的总接触力并不相同，而这是与事实不相符合的。也就是说，此处的势函数的定义物理意义是不明确的，无法对嵌入量这一指标基于统一的标准进行表征。为此，在接下来的部分，提出了基于统一标定的势接触力计算方法，采用了新的势函数定义，使得势可以基于统一标准对嵌入量进行表征，在相同嵌入量的情形下，计算出的接触力大小相等。

表1 基于原有势函数计算的接触力Table 1 Results of contact force with original potential function

3 统一标定的势接触力

3.1 基本思路

为了解决原有势接触力计算方法存在的问题，对势函数进行了重新定义，为了使得势能够对嵌入量基于统一标准进行表征，三角形单元内一点的势与该点到3条边的最短距离成正比，于是可以定义如下势函数：

式中：k为比例常数，hP-01、hP-12、hP-20分别表示点P 到三角形01边、12边、20边的垂直距离，如图4所示。

图4 三角形单元内一点P 的新势函数Fig.4 A new potential function at point P in a triangle element

为了使得整个离散系统的单元的势基于统一标准对嵌入量进行表征，同时保证势φ为无量纲量，于是定义一个标尺长度H，认为当一点的嵌入量（一点到三角形3条边的最大距离即为该点的嵌入量）为H 时，该点的势为1。于是据式（11）有k=1/H，这样势函数的表达式可进一步写为

此处需要特别说明，H为统一标定嵌入量，可以根据需要设定H 的值。为了使得法向罚参数 pn具有与传统DEM中法向刚度类似的意义，取H=1 m，这样法向罚参数就表示一点嵌入为1 m时，提供与法向罚参数大小相同的压力，H 也可设置为所有三角形单元中最短的高，来防止小尺寸的单元过大地嵌入。

为了便于理解，对式（12）作进一步的分析，依据三角形3条角平分线的交点（因为角平分线上的点到两边的距离相等），该交点即为三角形的内心I，将三角形分为3个子三角形，如图4所示。这样要求一点的势，只需先判断该点位于那个子三角型内，比如点P 位于I12子三角形内，则P 点势据式（12），可写为

据以上分析可知，三角形单元位于内心的点的势最大，但该点处的势不一定为1，同时不同三角形单元内心点处的势的大小不同，这与Munjiza定义的势函数在三角形内的中心处势最大，且不同三角形在中心处的势均为1是不同的。比较式（8）和式（13）可以发现，两者的主要不同在于，一个分母是 h0-12，一个分母是H，由于 h0-12的大小会随三角形的不同而不同，而采用H，则因为所有三角形单元都基于统一的标准，这样的势函数定义基于统一的标准表征了嵌入量，物理意义明确。同时，新定义的势函数在三角形的边界上势为0，完全满足势函数在边界上为常量这一要求。

采用新的势函数的定义，可以使得在嵌入量相同的情况下，计算得到的接触力大小相同。仍然以图3为例，据式（6）和式（13）对目标单元中长为L 的边积分得到的接触力为

由式（14）可知，对图3所示的左右两种接触情况，因为嵌入量相等，即hP-12相同，这样两种情况计算得到的接触力相同，此处仅对重叠区域的一条边界作了分析，其余边界也可类似分析，对图3左右两种接触情况，基于新的势函数定义计算得到的总的接触力相同。这一点可在第4节中算例1中得到进一步证实。

3.2 实施方案

由3.1节论述可知，本文仅仅是定义了新的势函数，再无其他更改，以图5所示，要求靶单元对目标单元AB 边的作用力，由于势函数的值在AB边上为分段线性函数，基于原有势函数定义，I 点为三角形的中心点，I2为三角形01边的中线，即势在I2边与AB 边的交点处发生转折，采用新的势函数的定义，只需将三角形的中心点I 的坐标替换为三角形的内心点的坐标，来求I2与AB 边的交点，并采用式（13）计算势即可。整个过程非常容易实施，只需对原有程序作非常小的改动即可。

图5 靶单元对目标单元AB 边的接触力计算Fig.5 Contact force of target element to edge AB of contactor element

如图5所示，在求得AB 边与角平分线I2的交点P1以及AB与三角形01边的交点P2后，依据式（5）计算靶单元对目标单元AB 边的接触力为

图5中，φt在AB 边上是分段线性分布的，故为图5中所示区域的面积，于是可进一步写为

φt(A)、φt(P1)、φt(P2)由式（14）可得

4 算 例

4.1 算例1

对图3所示的算例，采用基于统一标定的势接触力计算方法进行计算，采用的罚参数仍为pn=67.3 GPa 。对图3左右两组接触对计算得到对目标单元的总接触力（均已分配到节点上）如表2所示。由表可知，采用本文的方法计算，对左右两种接触情形，计算的接触力完全一致。

表2 基于新势函数计算的接触力Table 2 Results of contact force with new potential function

4.2 算例2

图3所示的算例给出的是一种较为特殊情况，该算例主要是便于进行理论分析，不能完全体现出基于统一标定的势接触力计算方法的全部优势，同时也无法完全暴露原有势函数定义所存在的不合理之处。这里给出一个更为通用的算例，如图6所示，采用Munjiza原有的势函数定义方法，在目标单元和靶单元重叠区域保持不变的条件下，靶单元的AB 边变为A′ B 边，或者CD 边变为CD′边，计算出的总的接触力都会发生变化，因为AB 边变为A′ B边，或者CD 边变为CD′边，使得三角形的高发生了变化，必然导致重叠区域的势发生变化，从而使得在重叠区域完全一致的情形下，计算出的总接触力并不相同，而这违背了基本的常识，是让人无法理解和不可接受的。正确的情况是，只要接触的局部区域相同，那么计算出的接触力大小就应该一致。如图7所示，只要重叠区域不变，不管相互接触三角形单元的另外一条边如何，或者说三角形单元的大小、整体形状如何，计算出的接触力大小都应该相等，本文提出的基于统一标定的势接触力计算方法，对图7所示的情况，计算的总接触力是完全相同的，因为重叠区域的势只与该点到三角形单元的两条边的距离有关，因此，重叠区域的势是固定的，因而计算出的总接触力完全一致。这说明新定义的势函数，计算出的接触力具有很好的局部特性，势基于统一标准对嵌入量进行了表征，解决了原有势函数定义存在的问题，新的势函数定义使得基于势的接触力计算有了坚实的物理基础，意义重大。

图6 嵌入量相同的4组接触Fig.6 Four group contacts with same embed amount

图7 相同嵌入量的一般化Fig.7 Generalization of same embed amount

为了进一步验证上述论述，对图6所示的4种接触情况，即ABE 与CDF 接触对、ABE 与CD F′ 接触对，即A′ BE 与CDF 接触对、A′ BE 与CD F′ 接触对，对这4组接触对的总接触力采用原有势接触力计算方法和本文的方法进行计算，由于三角形单元在变化，分配到各个节点上的力也在变化，为此仅仅比较作用在目标单元上总接触力，计算结果如表3所示。由表可知，原有势函数定义计算的接触力大小不等，而新定义的势函数计算的接触力完全一致。

表3 基于新、旧势函数计算的接触力比较Table 3 Comparison of contact forces between new and old potential functions

4.3 算例3

算例1和算例2均是对两种势接触力计算方法的验证和对比，此处给出一个应用性的算例，来充分展现基于统一标定的势接触力计算法的优点，由于采用统一标定，所有单元的接触力计算都基于同一标准，这样便不会出现局部某个接触对单元的接触力偏大或者偏小，而导致系统算出的接触力分布不均匀。由于采用统一标定的势函数定义，不会因为单元形状和尺寸的差异而导致接触力计算相差悬殊，对求解域的网格划分不再有严格的要求，而对原有方法则一般要求采用尺寸几乎相同的网格，否则因为尺寸的差异，导致局部即便嵌入相等，但接触力大小不一致。同时，单元的尺寸相差较大时，对小尺寸单元，相当于变相增加其刚度值，使得刚度参数失去了原有的物理意义。实际上，只要材料确定了，其刚度是不会改变的，当然不应随单元的尺寸和形状而改变。而新的势函数则不存在这一问题，刚度参数为材料的本质属性，不随单元尺寸和形状而改变。采用一个巴西圆盘劈裂的算例来加以说明，巴西圆盘的网格如图8所示，其网格尺寸和形状并不均一，用两种势接触力计算方法对该算例进行计算，由于两种方法采用不同的势函数定义，为了使得程序有相同参数输入，两者需要采用不同的法向罚参数，亦即采用原有程序计算时，刚度参数需要按照系统内所有三角形单元的高的平均值进行折算，已知所划分的三角形单元高平均值为0.24 m，原有程序输入的法向罚参数值，使得两者采用的罚参数不一致，是由于原有势函数定义方法存在变相增加刚度的问题，这样处理使得两种方法输入的实际等效刚度接近，其余的力学参数完全一致，只有这样才能对两者的计算结果进行比较。这种方法的模拟结果如图8～16所示。

由图8～16的计算结果可知，从应力分布的均匀程度上看，新方法获得应力分布更为均匀，应力云图中不同颜色的交界面更为光滑，同时在同一种颜色区域内部较少出现其他颜色的空洞，亦即代表应力值的颜色分布较均匀，如图9、10中的时步为10 000、15 000。从裂纹的扩展形态上看，新方法获得裂纹形态基本沿着圆盘的竖向直径扩展，与理论分析较为接近；而原有势接触力计算方法获得的裂纹出现了分支，圆盘下端出现了偏离竖向直径的裂纹，如图14、15中的时步为35 000、40 000所示。从两种方法计算得到荷载-位移曲线上看（如图16所示），两种方法获得的峰值强度基本一致，但新方法获得的荷载-位移曲线更为平滑，震荡幅度较小，这一点在曲线的峰后部分表现得尤为明显，这说明基于统一标定的势接触力计算方法要比原有势接触力方法计算更稳定。综上所述可知，基于新的势函数定义获得的应力分布也更均匀，较少出现应力集中，裂纹扩展形态与理论解更为接近，同时计算过程更为稳定，对网格划分的鲁棒性更高。

图8 两种方法计算结果的比较（时步为5 000）Fig.8 Comparison of two methods of calculation (step is 5 000)

图9 两种方法计算结果的比较（时步为10 000）Fig.9 Comparison of two methods of calculation (step is 10 000)

图10 两种方法计算结果的比较（时步为15 000）Fig.10 Comparison of two methods of calculation (step is 15 000)

图11 两种方法计算结果的比较（时步为20 000）Fig.11 Comparison of two methods of calculation (step is 20 000)

图12 两种方法计算结果的比较（时步为25 000）Fig.12 Comparison of two methods of calculation (step is 25 000)

图13 两种方法计算结果的比较（时步为30 000）Fig.13 Comparison of two methods of calculation (step is 30 000)

图14 两种方法计算结果的比较（时步为35 000）Fig.14 Comparison of two methods of calculation (step is 35 000)

图15 两种方法计算结果的比较（时步为40 000）Fig.15 Comparison of two methods of calculation (step is 40 000)

图16 巴西圆盘劈裂试验的荷载-位移曲线Fig.16 Load-displacement curves of the Brazilian disc tests

5 结 论

（1）指出了原有势函数定义，存在物理意义不明确，计算出的接触力与物理直观不符，同时势无法对嵌入量进行统一表征，接触力的计算没有建立在一个统一的标准之上等问题。

（2）提出了一种新的势函数定义方法，三角形单元内一点的势与该点到三角形单元3条边的最短距离成正比，三角单元内位于内心点处的势最大，但并不一定为1，且每个三角形内位于内心点处的势的值不一定相同。

（3）新提出的势函数，对嵌入量进行了统一的表征，物理意义明确、直观，使得基于势的接触力计算建立在坚实的物理基础之上，意义重大。

（4）新定义的势函数，解决原有势函数所存在的重大问题，只要重叠区域不变，计算的总的接触力将是不变的，接触力计算具有局部特性，而与划分的三角形单元的整体形状和大小无关，即新的接触力计算方法，具有更好的网格鲁棒性。

（5）基于统一标定的势接触力计算方法，计算出的接触力大小分布更为均匀，不会在局部产生过大的接触力，使得所有单元的接触力计算都基于一个统一的标准，不容易产生应力集中，计算更为稳定。

[1]CUNDALL P A.A computer model for simulating progressive large scale movement in block rock system[C]//Nancy:Symposium ISRM，1971，2:129－136.

[2]石根华.数值流形方法与非连续变形分析[M].裴觉民译.北京:清华大学出版社，1997:192－194.

[3]邬爱清，丁秀丽，卢波，等.DDA方法块体稳定性验证及其在岩质边坡稳定性分析中的应用[J].岩石力学与工程学报，2008，27(4):664－664.WU Ai-qing，DING Xiu-li，LU Bo，et al.Validation for rock block stability and its application to rock slope stability evaluation using DDA method[J].Chinese Journal of Rock Mechanics and Engineering，2008，27(4):664－664.

[4]郑宏，江巍.基于互补理论的非连续变形分析方法[J].中国科学(E辑:技术科学)，2009，39(10):1702－1708.ZHENG Hong，JIANG Wei.Discontinuous deformation analysis based on complementary theory[J].Science China (Series E:Technological Sciences)，2009，39(10):1702－1708.

[5]徐栋栋，孙冠华，郑宏.接触处理Munjiza方法在数值流形法中的应用[J].岩土力学，2013，34(2):526－534.XU Dong-dong，SUN Guan-hua，ZHENG Hong.Application of Munjiza’s algorithm to contact treatment in numerical manifold method[J].Rock and Soil Mechanics，2013，34(2):526－534.

[6]Itasca Consulting Group Ltd.User manual of UDEC code[S].Minneapolis:Itasca Consulting Group Ltd.，2011.

[7]陈文胜，王桂尧，刘辉，等.岩石力学离散单元计算方法中的若干问题探讨[J].岩石力学与工程学报，2005，24(10):1639－1644.CHEN Wen-sheng，WANG Gui-yao，LIU Hui，et al.Insight into some aspects of discrete element numerical methods for rock mass[J].Chinese Journal of Rock Mechanics and Engineering，2005，24(10):1639－1644.

[8]CUNDALL P A.Formulation of a three-dimensional distinct element model-part I:A scheme to detect and represent contact in composed of many polyhedral blocks[J].International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts，1988，25(3):107－116.

[9]MUNJIZA A，ANDREWS K R F.Penalty function method for combined finite-discrete element systems comprising large number of separate bodies[J].International Journal for Numerical Methods in Engineering，2000，49(11):1377－1369.

[10]MUNJIZA A.The combined finite-discrete element method[M].London:John Wiley &Sons，Ltd.，2004:29－32.

[11]严成增，孙冠华，郑宏，等.基于局部单元劈裂的FEM/DEM自适应分析方法[J].岩土力学，2014，35(7):2064－2070.YAN Cheng-zeng，SUN Guan-hua，ZHENG Hong，et al.Adaptive FEM/DEM analysis method based on the local splitting elements[J].Rock and Soil Mechanics，2014，35(7):2064－2070.

[12]严成增，孙冠华，郑宏，等.三维FEM/DEM中摩擦力的实施及验证[J].岩石力学与工程学报，2014，33(6):1248－1256.YAN Cheng-zeng，SUN Guan-hua，ZHENG Hong，et al.Implementation of friction to 3D FEM/DEM and verification[J].Chinese Journal of Rock Mechanics and Engineering，2014，33(6):1248－1256.