CN116822221B - 一种基于相互侵入势的离散岩块间接触力计算方法 - Google Patents

Abstract

本发明公开了一种基于相互侵入势的离散岩块间接触力计算方法，包括以下步骤：S1、使用软件建立离散岩块的初始几何模型，并计算指定参数；S2、根据初始几何模型和指定参数判断接触类型；S3、对不同接触类型的离散岩块，基于接触重叠域计算相互侵入势；S4、根据相互侵入势计算法向接触力和切向接触力；S5、根据法向接触力和切向接触力，得到外部荷载作用在相互接触的岩块上的接触力，完成相互侵入势的离散岩块间接触力计算，本发明比以往的分布式接触力模型具有更加严格的几何学基础，其计算相互侵入势与岩块的形状和岩块的形心无关。

Description

一种基于相互侵入势的离散岩块间接触力计算方法

技术领域

本发明涉及岩土工程和工程地质数值计算方法与模拟领域，特别涉及一种基于相互侵入势的离散岩块间接触力计算方法。

背景技术

在离散岩块稳定性分析与运动演化数值模拟中，由几何关系，二维岩块之间的接触可分为“凸凸接触”、“凹凸接触”、“凸边接触”和“边边接触”等四类。岩块间接触力的计算在数值模拟中非常重要。目前，主要采用集中式接触力模型和分布式接触力模型。当采用集中式接触力模型时，由于“凸凸接触”类似于“针尖对针尖”的接触，具有不稳定性。为保证数值模拟的顺利进行需要将“凸凸接触”转换为“凸边接触”；但是，如何合理地完成该转换学界还没有统一的标准和公认的解决方案，而分布式接触力模型可以规避这个棘手的问题。

由于岩块之间的接触总是出现在它们边界附近的局部区域，与岩块的形状、大小、姿态和形心位置等无关。但是，以往的分布式接触力模型具有岩块的形状和形心位置等依耐性，并且，当两个相互接触的岩块的尺寸相差较大时，即使接触重叠域的几何形状完全相同也可能计算出不同的接触力，这显然是不合理的。因此，需要发展跟合理的分布式接触力模型。

发明内容

针对现有技术中的上述不足，本发明提供的一种基于相互侵入势的离散岩块间接触力计算方法解决了现有技术求解岩块间接触力需要依赖岩块的形状和形心的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种基于相互侵入势的离散岩块间接触力计算方法，包括以下步骤：

S1、使用软件建立离散岩块的初始几何模型，并计算指定参数；

S2、根据初始几何模型和指定参数判断接触类型；

S3、对不同接触类型的离散岩块，基于接触重叠域计算相互侵入势；

S4、根据相互侵入势计算法向接触力和切向接触力；

S5、根据法向接触力和切向接触力，得到外部荷载作用在相互接触的岩块上的接触力，完成相互侵入势的离散岩块间接触力计算。

进一步地：当岩块1上凸角V1与岩块2上凸角V2之间的距离d＜2d₀时，凸角V1与凸角V2构成凸凸接触；

当岩块1上凹角V1与岩块2上凸角V2之间的距离d＜2d₀时，凹角V1与凸角V2构成凹凸接触；

当岩块1上凸角V1与岩块2上某条边之间的距离d＜2d₀时，凸角V1与这条边构成凸边接触；

当岩块1上某条边与岩块2上另一条边近似平行，且它们之间的距离d₁＜2d₀时，这两条边构成边边接触；

其中，d₀为距离阈值。

进一步地：所述步骤S3中，两岩块接触时，基于接触重叠域的相互侵入势分为第一侵入势和第二侵入势/>

第一侵入势为岩块1侵入到岩块2之中，岩块2对岩块1的“抵抗”，其定义在岩块1上；

第二侵入势为岩块2侵入到岩块1之中，岩块1对岩块2的“抵抗”，其定义在岩块2上。

进一步地：所述步骤S3中，当接触类型为凸凸接触时，相互侵入势计算包括以下步骤：

S301、找到构成两凸角的四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S302、测量岩块1上凸角V₁到岩块2上凸角V₂的两条边的距离d₁₁和d₁₂；

S303、根据d₁₁和d₁₂计算凸角V₁的第一侵入势其公式如下：

S304、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S305、测量岩块2上凸角V₂到岩块1上凸角V₁的两条边的距离d₂₁和d₂₂；

S306、根据d₂₁和d₂₂计算凸角V₂的第二侵入势其公式如下：

S307、根据和/>得到线段AV₂上的第二侵入势/>和线段BV₂上的第二侵入势/>

其中，线段AV₂上的第二侵入势为从/>到/>呈线性分布；线段BV₂上的第二侵入势/>为从/>到/>呈线性分布。

进一步地：所述步骤S3中，当接触类型为凹凸接触时，相互侵入势计算包括以下步骤：

S311、找到构成凹角和凸角的四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S312、测量岩块1上凹角V₁到岩块2上凸角V₂的两条边的距离d₁₁和d₁₂；

S313、根据d₁₁和d₁₂计算凹角V₁的第一侵入势其公式如下：

S314、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S315、测量岩块2上凸角V₂到岩块1上凹角V₁的两条边的距离d₂₁和d₂₂；

S316、根据d₂₁和d₂₂计算凸角V₂的第二侵入势其公式如下：

S317、根据和/>得到线段AV₂上的第二侵入势/>和线段BV₂上的第二侵入势/>

其中，线段AV₂上的第二侵入势为从/>到/>呈线性分布；线段BV₂上的第二侵入势/>为从/>到/>呈线性分布。

进一步地：所述步骤S3中，当接触类型为凸边接触时，相互侵入势计算包括以下步骤：

S321、找到构成岩块1上凸角的两条边与岩块2上接触边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S322、测量岩块1上凸角V₁到岩块2上接触边的距离d₁；

S323、根据d₁计算凸角V₁的第一侵入势其公式如下：

S324、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S325、测量岩块2上线段AB的中点M到凸角V₁上两条边的三条边的距离d₂₁和d₂₂；

S326、根据d₂₁和d₂₂计算线段AB的中点M的第二侵入势其公式如下：

S327、根据和/>得到线段AM上的第二侵入势/>和线段BM上的第二侵入势/>

其中，线段AM上的第二侵入势为从/>到/>呈线性分布；线段BM上的第二侵入势/>为从/>到/>呈线性分布。

进一步地：所述步骤S3中，当接触类型为边边接触时，相互侵入势计算包括以下步骤：

S331、找到构成岩块1与岩块2上相关四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S332、测量岩块1上接触边的端点C到岩块2上接触边距离d_1C；

S333、根据d_1C计算端点C的第一侵入势其公式如下：

S334、测量岩块1上接触边的端点D到岩块2上接触边距离d_1D；

S335、根据d_1D计算端点D的第一侵入势其公式如下：

S336、根据和/>得到线段AC上的第一侵入势/>线段AD上的第一侵入势/>和线段CD上的第一侵入势/>

其中，线段AC上的第一侵入势为从/>到/>呈线性分布；线段AD上的第一侵入势/>为从/>到/>呈线性分布；线段CD上的第一侵入势/>呈线性分布；

S337、测量岩块2上线段AB的中点M到岩块1上端点C和端点D的三条边的距离d₂₁、d₂₂和d₂₃；

S338、根据d₂₁、d₂₂和d₂₃计算线段AB的中点M的第二侵入势其公式如下：

S339、根据和/>得到线段AM上的第二侵入势/>和线段BM上的第二侵入势/>

其中，线段AM上的第二侵入势呈线性分布；线段BM上的第二侵入势呈线性分布。

进一步地：所述步骤S4包括以下分步骤：

S41、根据A、B两点的第一侵入势和第二侵入势，计算起点为A、终点为B的有向线段AB上任意一点Q的第一侵入势和第二侵入势，其公式如下：

其中，为Q的第一侵入势，/>为Q的第二侵入势，a＝|AQ|/|AB|，为点Q在线段AB上的位置参数；

S42、计算Q的第一侵入势和第二侵入势导致的第一法向接触力和第二法向接触力，其公式如下：

其中，df_n1为Q点处的第一法向接触力，df_n2为Q点处的第二法向接触力，p_n为法向接触弹簧值；

S43、根据Q点处的第一法向接触力和第二法向接触力，计算整条有向线段AB上的第一、二侵入势导致的第一、二法向接触力，其公式如下：

f_n1＝∫_ABdf_n1

f_n2＝∫_ABdf_n2

其中，f_n1为有向线段AB上的第一侵入势导致的第一法向接触力，f_n2为有向线段AB上的第二侵入势导致的第二法向接触力；

S44、根据f_n1和f_n2，得到有向线段AB上的第一切向接触力和第二切向接触力，其公式如下：

f_s1＝μf_n1

f_s2＝μf_n2

其中，f_s1为有向线段AB上的第一切向接触力，f_s2为有向线段AB上的第二切向接触力，μ为摩擦系数；

S45、根据有向线段AB的中点M的瞬时速度，得到作用在岩块2上的切向接触力的方向sign₂，其表示式如下：

其中，+表示从A点指向B点；-表示从B点指向A点；和/>分别表示岩块1和岩块2上的M点的瞬时速度在有向线段AB上的投影。

进一步地：所述步骤S5外部荷载作用在相互接触的岩块上的接触力包括：

作用在岩块1上的因有向线段AB的第一侵入势导致的第一法向和切向接触力、作用在岩块2上的因有向线段AB的第一侵入势导致的第一法向和切向接触力的反作用力、作用在岩块2上的因有向线段AB的第二侵入势导致的第二法向和切向接触力以及作用在岩块1上的因有向线段AB的第二侵入势导致的第二法向和切向接触力的反作用力。

本发明的有益效果为：

1.比现有的分布式接触力模型具有更加严格的几何学基础，其相互侵入势与岩块的形状和岩块的形心无关；

2.相互侵入势具有彻底的局部性质，当且仅当接触重叠域几何全等时，由相互侵入势计算的接触力才相同。

附图说明

图1为基于相互侵入势的离散岩块间接触力计算方法的流程图。

图2为二维岩块间的四种接触类型。

图3为凸凸接触的第一、二侵入势示意图。

图4为凹凸接触的第一、二侵入势示意图。

图5为凸边接触的第一、二侵入势示意图。

图6为边边接触的第一、二侵入势示意图。

图7为本发明中由侵入势导致的法向接触力的示意图。

图8为本发明中切向接触力方向示意图。

图9为本发明具体应用实例1提供的岩块斜面滑动验证算例的几何模型。

图10为本发明具体应用实例1提供的滑动岩块形心距离随时间变化的数值解和解析解对比图。

图11为本发明具体应用实例2提供的岩块倾倒失稳运动演化算例的几何模型。

图12为本发明具体应用实例2提供的岩块倾倒失稳运动演化算例的计算结果。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，在本发明的一个实施例中，提供一种基于相互侵入势的离散岩块间接触力计算方法，包括以下步骤：

S1、使用软件建立离散岩块的初始几何模型，并计算指定参数；

S2、根据初始几何模型和指定参数判断接触类型；

S3、对不同接触类型的离散岩块，基于接触重叠域计算相互侵入势；

两岩块接触时，基于接触重叠域的相互侵入势分为第一侵入势和第二侵入势

第一侵入势为岩块1侵入到岩块2之中，岩块2对岩块1的“抵抗”，其定义在岩块1上；

第二侵入势为岩块2侵入到岩块1之中，岩块1对岩块2的“抵抗”，其定义在岩块2上；

S4、根据相互侵入势计算法向接触力和切向接触力；

S5、根据法向接触力和切向接触力，得到外部荷载作用在相互接触的岩块上的接触力，完成相互侵入势的离散岩块间接触力计算。

如图2所示，在本实施例中，当岩块1上凸角V1与岩块2上凸角V2之间的距离d＜2d₀时，凸角V1与凸角V2构成凸凸接触；

当岩块1上凹角V1与岩块2上凸角V2之间的距离d＜2d₀时，凹角V1与凸角V2构成凹凸接触；

当岩块1上凸角V1与岩块2上某条边之间的距离d＜2d₀时，凸角V1与这条边构成凸边接触；

当岩块1上某条边与岩块2上另一条边近似平行，且它们之间的距离d₁＜2d₀时，这两条边构成边边接触；

其中，d₀为距离阈值。

如图3所示，当接触类型为凸凸接触时，相互侵入势计算包括以下步骤：

S301、找到构成两凸角的四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S302、测量岩块1上凸角V₁到岩块2上凸角V₂的两条边的距离d₁₁和d₁₂；

S303、根据d₁₁和d₁₂计算凸角V₁的第一侵入势其公式如下：

S304、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S305、测量岩块2上凸角V₂到岩块1上凸角V₁的两条边的距离d₂₁和d₂₂；

S306、根据d₂₁和d₂₂计算凸角V₂的第二侵入势其公式如下：

S307、根据和/>得到线段AV₂上的第二侵入势/>和线段BV₂上的第二侵入势/>

其中，线段AV₂上的第二侵入势为从/>到/>呈线性分布；线段BV₂上的第二侵入势/>为从/>到/>呈线性分布。

如图4所示，当接触类型为凹凸接触时，相互侵入势计算包括以下步骤：

S311、找到构成凹角和凸角的四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S312、测量岩块1上凹角V₁到岩块2上凸角V₂的两条边的距离d₁₁和d₁₂；

S313、根据d₁₁和d₁₂计算凹角V₁的第一侵入势其公式如下：

S314、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S315、测量岩块2上凸角V₂到岩块1上凹角V₁的两条边的距离d₂₁和d₂₂；

S316、根据d₂₁和d₂₂计算凸角V₂的第二侵入势其公式如下：

S317、根据和/>得到线段AV₂上的第二侵入势/>和线段BV₂上的第二侵入势/>

其中，线段AV₂上的第二侵入势为从/>到/>呈线性分布；线段BV₂上的第二侵入势/>为从/>到/>呈线性分布。

如图5所示，当接触类型为凸边接触时，相互侵入势计算包括以下步骤：

S321、找到构成岩块1上凸角的两条边与岩块2上接触边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S322、测量岩块1上凸角V₁到岩块2上接触边的距离d₁；

S323、根据d₁计算凸角V₁的第一侵入势其公式如下：

S324、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S325、测量岩块2上线段AB的中点M到凸角V₁上两条边的三条边的距离d₂₁和d₂₂；

S326、根据d₂₁和d₂₂计算线段AB的中点M的第二侵入势其公式如下：

S327、根据和/>得到线段AM上的第二侵入势/>和线段BM上的第二侵入势/>

其中，线段AM上的第二侵入势为从/>到/>呈线性分布；线段BM上的第二侵入势/>为从/>到/>呈线性分布。

如图6所示，当接触类型为边边接触时，相互侵入势计算包括以下步骤：

S331、找到构成岩块1与岩块2上相关四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S332、测量岩块1上接触边的端点C到岩块2上接触边距离d_1C；

S333、根据d_1C计算端点C的第一侵入势其公式如下：

S334、测量岩块1上接触边的端点D到岩块2上接触边距离d_1D；

S335、根据d_1D计算端点D的第一侵入势其公式如下：

S336、根据和/>得到线段AC上的第一侵入势/>线段AD上的第一侵入势/>和线段CD上的第一侵入势/>

其中，线段AC上的第一侵入势为从/>到/>呈线性分布；线段AD上的第一侵入势/>为从/>到/>呈线性分布；线段CD上的第一侵入势/>呈线性分布；

S337、测量岩块2上线段AB的中点M到岩块1上端点C和端点D的三条边的距离d₂₁、d₂₂和d₂₃；

S338、根据d₂₁、d₂₂和d₂₃计算线段AB的中点M的第二侵入势其公式如下：

S339、根据和/>得到线段AM上的第二侵入势/>和线段BM上的第二侵入势/>

其中，线段AM上的第二侵入势呈线性分布；线段BM上的第二侵入势呈线性分布。

所述步骤S4包括以下分步骤：

S41、根据A、B两点的第一侵入势和第二侵入势，计算起点为A、终点为B的有向线段AB上任意一点Q的第一侵入势和第二侵入势，其公式如下：

其中，为Q的第一侵入势，/>为Q的第二侵入势，a＝|AQ|/|AB|，为点Q在线段AB上的位置参数；

S42、计算Q的第一侵入势和第二侵入势导致的第一法向接触力和第二法向接触力，其公式如下：

其中，df_n1为Q点处的第一法向接触力，df_n2为Q点处的第二法向接触力，p_n为法向接触弹簧值；

S43、如图7所示，根据Q点处的第一法向接触力和第二法向接触力，计算整条有向线段AB上的第一、二侵入势导致的第一、二法向接触力，其公式如下：

f_n1＝∫_ABdf_n1

f_n2＝∫_ABdf_n2

其中，f_n1为有向线段AB上的第一侵入势导致的第一法向接触力，f_n2为有向线段AB上的第二侵入势导致的第二法向接触力；

S44、根据f_n1和f_n2，得到有向线段AB上的第一切向接触力和第二切向接触力，其公式如下：

f_s1＝μf_n1

f_s2＝μf_n2

其中，f_s1为有向线段AB上的第一切向接触力，f_s2为有向线段AB上的第二切向接触力，μ为摩擦系数；

S45、如图8所示，根据有向线段AB的中点M的瞬时速度，得到作用在岩块2上的切向接触力的方向sign₂，其表示式如下：

其中，+表示从A点指向B点；-表示从B点指向A点；和/>分别表示岩块1和岩块2上的M点的瞬时速度在有向线段AB上的投影。

所述步骤S5外部荷载作用在相互接触的岩块上的接触力包括：

作用在岩块1上的因有向线段AB的第一侵入势导致的第一法向和切向接触力、作用在岩块2上的因有向线段AB的第一侵入势导致的第一法向和切向接触力的反作用力、作用在岩块2上的因有向线段AB的第二侵入势导致的第二法向和切向接触力以及作用在岩块1上的因有向线段AB的第二侵入势导致的第二法向和切向接触力的反作用力。

具体应用示例1：如图9所示，建立该应用实例的离散岩块的初始几何模型，可滑动岩块尺寸为0.5m×1.0m，固定岩块斜面倾角α＝35°。并设置计算参数如下：密度为1,000kg/m3，杨氏模量为10GPa，泊松比为0.35。重力加速度为-10m/s2，接触弹簧刚度为300GPa。时间步长为0.00001s，总时间为1.5s。滑面摩擦角分别采用φ＝0°、5°、15°、20°、25°和35°等不同的值进行算法验证。

该算例中，滑动岩块形心沿斜面的滑移距离d的解析表达式为：

滑动岩块形心距离d随时间t变化的数值解对比如图10所示，从该图中可以看出：对于滑面摩擦角分别为φ＝0°、5°、15°、20°、25°和35°时，相应的离散数据点均落在连续曲线上，表明本发明所提方法获得的数值解具有很好的精确度。

具体应用示例2：如图11所示，建立该应用实例的离散岩块的初始几何模型，该算例由10个岩块组成，分别标记为B1、B2、…、B10。其中，岩块B1与水平面的夹角为30°且固定不动，它们的几何尺寸参见图10。设置计算参数如下：密度为1,000kg/m3，杨氏模量为10GPa，泊松比为0.30。重力加速度为-10m/s2，接触弹簧刚度为100GPa。时间步长为0.00001s，总计算时间为0.5s，不考虑接触面的摩擦。

具体应用示例2没有解析解，由本发明所提方法得到数值模拟结果如图12所示。从该图可以看出，当总计算时间为0.5s时，岩块B2、B3、…、B8出现了倾倒失稳，岩块B9随着岩块B4和岩块B5的倾斜向右下方做相应的运动。另一方面，由于岩块B6和岩块B7向右倾斜，岩块B10逐渐脱离它们，呈自由落体之势。随着岩块B2、B3、…、B10整体沿岩块B1向右下方滑动，最右边的岩块B8与岩块B1之间出现了“凸边接触”。由于岩块B1固定不动，这对接触中的“凸角”就成为这些岩块产生旋转运动的支点。该算例表明本发明提出的方法能模拟离散岩块系统的失稳启动和运动演化全过程。

在本发明的描述中，需要理解的是，术语“中心”、“厚度”、“上”、“下”、“水平”、“顶”、“底”、“内”、“外”、“径向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的设备或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性或隐含指明的技术特征的数量。因此，限定由“第一”、“第二”、“第三”的特征可以明示或隐含地包括一个或者更多个该特征。

Claims (3)

1.一种基于相互侵入势的离散岩块间接触力计算方法，其特征在于，包括以下步骤：

S1、使用软件建立离散岩块的初始几何模型，并计算指定参数；

S2、根据初始几何模型和指定参数判断接触类型；

所述步骤S2中，判断接触类型的方法为：

当岩块1上凸角V1与岩块2上凸角V2之间的距离d＜2d₀时，凸角V1与凸角V2构成凸凸接触；

当岩块1上凹角V1与岩块2上凸角V2之间的距离d＜2d₀时，凹角V1与凸角V2构成凹凸接触；

当岩块1上凸角V1与岩块2上某条边之间的距离d＜2d₀时，凸角V1与这条边构成凸边接触；

当岩块1上某条边与岩块2上另一条边近似平行，且它们之间的距离d₁＜2d₀时，这两条边构成边边接触；

其中，d₀为距离阈值；

S3、对不同接触类型的离散岩块，基于接触重叠域计算相互侵入势；

所述步骤S3中，两岩块接触时，基于接触重叠域的相互侵入势分为第一侵入势和第二侵入势/>

第一侵入势为岩块1侵入到岩块2之中，岩块2对岩块1的“抵抗”，其定义在岩块1上；

第二侵入势为岩块2侵入到岩块1之中，岩块1对岩块2的“抵抗”，其定义在岩块2上；

当接触类型为凸凸接触时，相互侵入势计算包括以下步骤：

S301、找到构成两凸角的四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S302、测量岩块1上凸角V₁到岩块2上凸角V₂的两条边的距离d₁₁和d₁₂；

S303、根据d₁₁和d₁₂计算凸角V₁的第一侵入势其公式如下：

S304、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S305、测量岩块2上凸角V₂到岩块1上凸角V₁的两条边的距离d₂₁和d₂₂；

S306、根据d₂₁和d₂₂计算凸角V₂的第二侵入势其公式如下：

S307、根据和/>得到线段AV₂上的第二侵入势/>和线段BV₂上的第二侵入势/>

其中，线段AV₂上的第二侵入势为从/>到/>呈线性分布；线段BV₂上的第二侵入势/>为从/>到/>呈线性分布；

所述步骤S3中，当接触类型为凹凸接触时，相互侵入势计算包括以下步骤：

S311、找到构成凹角和凸角的四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S312、测量岩块1上凹角V₁到岩块2上凸角V₂的两条边的距离d₁₁和d₁₂；

S313、根据d₁₁和d₁₂计算凹角V₁的第一侵入势其公式如下：

S314、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S315、测量岩块2上凸角V₂到岩块1上凹角V₁的两条边的距离d₂₁和d₂₂；

S316、根据d₂₁和d₂₂计算凸角V₂的第二侵入势其公式如下：

S317、根据和/>得到线段AV₂上的第二侵入势/>和线段BV₂上的第二侵入势/>

其中，线段AV₂上的第二侵入势为从/>到/>呈线性分布；线段BV₂上的第二侵入势/>为从/>到/>呈线性分布；

当接触类型为凸边接触时，相互侵入势计算包括以下步骤：

S321、找到构成岩块1上凸角的两条边与岩块2上接触边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S322、测量岩块1上凸角V₁到岩块2上接触边的距离d₁；

S323、根据d₁计算凸角V₁的第一侵入势其公式如下：

S324、根据和/>得到线段AV₁上的第一侵入势/>和线段BV₁上的第一侵入势/>

其中，线段AV₁上的第一侵入势为从/>到/>呈线性分布；线段BV₁上的第一侵入势/>为从/>到/>呈线性分布；

S325、测量岩块2上线段AB的中点M到凸角V₁上两条边的三条边的距离d₂₁和d₂₂；

S326、根据d₂₁和d₂₂计算线段AB的中点M的第二侵入势其公式如下：

S327、根据和/>得到线段AM上的第二侵入势/>和线段BM上的第二侵入势/>

其中，线段AM上的第二侵入势为从/>到/>呈线性分布；线段BM上的第二侵入势/>为从/>到/>呈线性分布；

当接触类型为边边接触时，相互侵入势计算包括以下步骤：

S331、找到构成岩块1与岩块2上相关四条边的交点A和交点B，并定义A、B两点的第一侵入势和A、B两点的第二侵入势；

其中， 为交点A的第一侵入势，/>为交点B的第一侵入势，/>为交点A的第二侵入势，/>为交点B的第二侵入势；

S332、测量岩块1上接触边的端点C到岩块2上接触边距离d_1C；

S333、根据d_1C计算端点C的第一侵入势其公式如下：

S334、测量岩块1上接触边的端点D到岩块2上接触边距离d_1D；

S335、根据d_1D计算端点D的第一侵入势其公式如下：

S336、根据和/>得到线段AC上的第一侵入势/>线段AD上的第一侵入势/>和线段CD上的第一侵入势/>

其中，线段AC上的第一侵入势为从/>到/>呈线性分布；线段AD上的第一侵入势/>为从/>到/>呈线性分布；线段CD上的第一侵入势/>呈线性分布；

S337、测量岩块2上线段AB的中点M到岩块1上端点C和端点D的三条边的距离d₂₁、d₂₂和d₂₃；

S338、根据d₂₁、d₂₂和d₂₃计算线段AB的中点M的第二侵入势其公式如下：

S339、根据和/>得到线段AM上的第二侵入势/>和线段BM上的第二侵入势/>

其中，线段AM上的第二侵入势呈线性分布；线段BM上的第二侵入势/>呈线性分布；

S4、根据相互侵入势计算法向接触力和切向接触力；

S5、根据法向接触力和切向接触力，得到外部荷载作用在相互接触的岩块上的接触力，完成相互侵入势的离散岩块间接触力计算。

2.根据权利要求1所述的基于相互侵入势的离散岩块间接触力计算方法，其特征在于，所述步骤S4包括以下分步骤：

S41、根据A、B两点的第一侵入势和第二侵入势，计算起点为A、终点为B的有向线段AB上任意一点Q的第一侵入势和第二侵入势，其公式如下：

其中，为Q的第一侵入势，/>为Q的第二侵入势，a＝|AQ|/|AB|，为点Q在线段AB上的位置参数；

S42、计算Q的第一侵入势和第二侵入势导致的第一法向接触力和第二法向接触力，其公式如下：

其中，df_n1为Q点处的第一法向接触力，df_n2为Q点处的第二法向接触力，p_n为法向接触弹簧值；

S43、根据Q点处的第一法向接触力和第二法向接触力，计算整条有向线段AB上的第一、二侵入势导致的第一、二法向接触力，其公式如下：

f_n1＝∫_ABdf_n1

f_n2＝∫_ABdf_n2

其中，f_n1为有向线段AB上的第一侵入势导致的第一法向接触力，f_n2为有向线段AB上的第二侵入势导致的第二法向接触力；

S44、根据f_n1和f_n2，得到有向线段AB上的第一切向接触力和第二切向接触力，其公式如下：

f_s1＝μf_n1

f_s2＝μf_n2

其中，f_s1为有向线段AB上的第一切向接触力，f_s2为有向线段AB上的第二切向接触力，μ为摩擦系数；

S45、根据有向线段AB的中点M的瞬时速度，得到作用在岩块2上的切向接触力的方向sign₂，其表示式如下：

其中，+表示从A点指向B点；-表示从B点指向A点；和/>分别表示岩块1和岩块2上的M点的瞬时速度在有向线段AB上的投影。

3.根据权利要求2所述的基于相互侵入势的离散岩块间接触力计算方法，其特征在于，所述步骤S5外部荷载作用在相互接触的岩块上的接触力包括：

作用在岩块1上的因有向线段AB的第一侵入势导致的第一法向和切向接触力、作用在岩块2上的因有向线段AB的第一侵入势导致的第一法向和切向接触力的反作用力、作用在岩块2上的因有向线段AB的第二侵入势导致的第二法向和切向接触力以及作用在岩块1上的因有向线段AB的第二侵入势导致的第二法向和切向接触力的反作用力。