CN105912852A - 一种基于距离势函数任意凸多边形块体离散单元法 - Google Patents

Abstract

本发明公开了一种基于距离势函数任意凸多边形块体离散单元法，方法主要包括以下步骤：先建立矿石‑边界离散模型；采用NBS接触检测方法对所有块体单元进行接触检测，得到每个块体单元及与之接触的块体单元；对相互接触的块体单元通过块体离散单元距离势函数进行接触力计算，得到作用于块体单元的接触力以及力矩；进而通过velocity verlet算法计算块体单元的速度与位移，最后逐步更新矿石块体单元的位移即是其具体运动形态。本发明采用距离势函数的定义，明确了势函数的物理意义，通过基于距离势函数的块体离散单元法实现计算任意凸多边形块体的运动形态，可以用于岩土矿石的运动形态研究，为实际矿场工程提供实际指导意义。

Description

一种基于距离势函数任意凸多边形块体离散单元法

技术领域

本发明涉及一种块体离散单元法，具体涉及一种岩土工程中通过基于距离势函数的块体离散单元法实现计算二维非连续介质的大规模任意凸多边形块体的运动形态，属于块体离散元模型技术领域。

背景技术

离散单元法的思想源于较早的分子动力学，1971年由美国的Cundall P.A.教授首先提出用以模拟非连续介质力学行为的数值分析方法。这种方法是将研究对象划分为大量的离散块体(颗粒)，每个块体拥有相应的自由度。块体之间通过接触力、力矩以及摩擦力相互作用，同时施加合适的力或位移的边界条件。在给定初始条件及外力作用下，基于牛顿第二运动定律并采用显示时步步进的动态松弛法，确定下一时刻块体的运动状态。

离散单元法原理直观简单，由于离散元单元具有更真实地表达岩体的几何力学特性的能力，便于处理以所有非线性变形和破坏都集中的节理面上为特征的岩体破坏问题，被广泛地应用于模拟边坡、滑坡、矿场等力学过程的分析和计算中。从本质上来说，岩土材料都是由离散的、尺寸不一、形状各异的颗粒或块体组成。离散元法的单元从几何形状上分类可分为颗粒元和块体元两大类，块体元中对于二维问题可以是任意多边形元。矿场中矿石的运动形态采用块体离散元法，它的基本思想是把不连续体分离为刚性块体的集合，使各个块体满足运动方程，用时步迭代的方法求解各块体的运动方程，继而求得不连续体的整体运动形态。

本申请主要运用块体离散元法研究矿场中矿石在一定形状边界内的运动形态，如矿石以一定速度落入槽框内的运动形态，或者某放矿场漏斗内充满了已崩大小不等、不规则的矿石，漏斗内嵌在槽框内，漏斗与槽框分别是固定的，研究对象为，将漏斗闸门打开，矿石从漏孔中放出后落在槽框中的矿石的整体运动形态等。研究矿石的运动形态主要是研究矿石运动过程中块体离散单元之间的接触力，块体离散单元法中接触力主要有角-角接触力、角-边接触力和边-边接触力。目前，对于块体离散单元法中接触力计算主要有两种：一是引入刚度系数，通过嵌入量计算方法，该方法在计算中可能存在能量不守恒的问题，同时在处理角-角接触时，采用圆角化进行规避，过程复杂；二是由英国A.MUNJIZA教授提出的基于势函数接触力计算方法，通过将研究对象划分为三角形块体单元，并以单元形心为基础建立势函数定义，该方法能够保证系统能量守恒，同时很好地规避了角-角接触的问题，但是存在着势函数物理意义不明，对于求解任意形状的凸多边形块体接触力的情况，更是无能为力。

综合目前的研究，本发明提供一种新的势函数法，解决平面任意凸多边形块体离散单元接触力计算的问题，并明确了势函数的物理意义，有利于研究矿石在一定形状边界内的运动过程形态。

发明内容

本发明的目的在于克服现有技术中的不足，提供了一种基于距离势函数任意凸多边形块体离散单元法，解决了现有技术中势函数计算方法中物理意义不明确的技术问题。

为解决上述技术问题，本发明提供了一种基于距离势函数任意凸多边形块体离散单元法，包括以下步骤：

步骤一，建立矿石-边界离散模型；基于离散单元法理论，对在一定形状边界内的多个矿石，建立矿石-边界的离散元模型，矿石模型为多个多边形块体单元，多边形块体单元的参数包括多边形的质量、块体间的阻尼系数及块体单元的刚度系统；

步骤二，确定模型采集矿石信息的迭代时间步长；

步骤三，在当前时间步，采用NBS接触检测方法对所有块体单元进行接触检测，得到每个块体单元及与之接触的块体单元；

步骤四，根据步骤三的检测结果，对相互接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体(Target)，与此块体单元相接触的块体单元为接触块体(Contactor)，基于所提出的块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；

步骤五，重复步骤四计算当前时间步内每个块体单元所受的接触力和力矩；

步骤六，根据牛顿第二定律，计算出每个块体单元当前时间步的加速度，再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移；

步骤七，根据步骤六中块体单元的位移，更新每个块体单元顶点、形心及内点的坐标，完成当前时间步的计算；

步骤八，重复步骤三至七计算下一时间步，直至计算完所有时间步。

进一步的，在所述步骤二中，模拟采集矿石信息离散单元计算时间步长为

$\mathrm{\Δ}t\≤2\·\sqrt{\frac{m}{k}}\·(\sqrt{1+{\mathrm{\ξ}}^2}-\mathrm{\ξ})$

其中：ξ是系统的阻尼比，m是块体单元质量，c是阻尼系数，k是刚度系数。

进一步的，所述步骤三中，NBS接触检测法进行接触检测的具体步骤包括：

步骤3.1，对每个块体单元，计算其形心坐标及各顶点到形心之间的距离l；

步骤3.2，选取形心到顶点之间的最大距离为l_max，记d＝2·l_max，将离散模型计算区域划分为以d为边长的正方形网格；

步骤3.3，根据公式：

$x_k=1+\mathrm{int}(\frac{x_k-x_{min}}{l_{max}}+0.5)$

$y_k=1+\mathrm{int}(\frac{y_k-y_{\min }}{l_{max}}+0.5)$

将所有块体单元映射到网格上，其中，x_k、y_k为块体单元k的形心坐标，k＝1,2,……,N，N为块体单元个数，x_min、y_min为计算区域在x、y方向的最小值；

步骤3.4，对每个块体单元与其周边所有相邻网格内块体单元进行接触检测，当循环到某个块体单元时，即以此块体单元所在格子为中心，进行判断此块体单元是否与该格子内部以及相邻格子内部单元接触，判断此块体单元与相邻的一个块体单元是否接触的具体过程为：

1)分别计算两个块体单元的形心之间的距离dist、两块体单元各顶点与其对应的形心之间的最大距离分别为l_1max、l_2max，计算l_1max、l_2max的平均值为l_avmax；

2)比较dist与l_avmax的大小，如果dist≤2l_avmax，则判断此两个块体单元接触，如果dist>2l_avmax，则判断此两个块体单元不接触；

步骤3.5，遍历所有块体单元，记录每个块体单元及与之接触的块体单元。

进一步的，在所述步骤四中，目标块体的法向接触力、切向接触力以及力矩具体计算步骤如下：

步骤4.1，确定目标块体各顶点对应的内点，其中内点定义为在多边形顶点角平分线上至顶点所在两条边距离为h的点为多边形顶点所对应的内点，其中h为多边形最大内切圆直径：

1)根据最小二分法，取h_min＝0，h_max＝2l_max迭代计算，其中l_max为目标块体的形心到顶点间的最大距离，确定目标块体最大内切圆直径h；

2)根据点到直线的距离公式：求出各顶点对应的内点坐标；其中，A、B、C是顶点所在两条边任一直线方程参数；x₀、y₀是顶点对应的内点坐标；

步骤4.2，根据步骤4.1所求目标块体的内点，连接顶点及其对应的内点将目标块体划分为n个子块体，n为目标块体的顶点数；

步骤4.3，定义距离势函数，循环计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩；

步骤4.4，更新所有与目标块体相接触的块体单元作用于目标块体的接触力和力矩。

进一步的，在所述步骤4.3中，计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩，以相互接触的块体单元1和块体单元2说明计算过程，对此两个彼此接触的块体单元，以块体单元1为目标单元，求解由块体单元2所引起、作用于目标单元的接触力和力矩过程为：

1)计算块体单元2各底边与块体单元1各子块边线的交点，确定块体单元2与块体单元1的重叠区域；

2)定义距离势函数：其中，是块体单元内部点p的势函数值，h_p-12、h_p-23……h_p-n1是点p至块体单元各子块底边的距离，如h_p-12是点p至块体单元的12边的距离；n为块体单元1的顶点数；h为块体单元最大内切圆半径；计算出各交点在块体单元1子块中的距离势函数值，若交点在块体单元1底边上则势为0，在块体单元1内点与顶点连线上则势为1；

3)根据距离势函数值在子块体内呈线性分布的特征，分别计算各子块的重叠区域边界势函数值：其中，为边界的势函数值，i＝1，……，M-1，M为交点以及块体单元2在块体单元1内的顶点个数之和，为交点p_i的势函数值，为交点p_i+1的势函数值，l_pipi+1为相邻两交点p_i和p_i+1之间的距离；

4)计算当前时间步，块体单元1内子块的法向接触力及力矩，具体计算过程为：

4-1)基于重叠区域边界势函数法向接触力计算公式为：其中，f_n是重叠区域边界上的法向接触力；n_Γ是边界的外法向向量，即法向接触力方向矢量；p_n为罚函数，为边界上各点的势函数值；根据3)所求的重叠区域边界势函数值，法向接触力计算公式简化为计算作用于重叠区域边界的法向接触力；

4-2)力的作用点取重叠区域边界势函数分布图形心在对应边界上的投影点A；

4-3)由块体单元2侵入块体单元1，所引起的作用于块体单元2的法向接触力为f_n；根据牛顿第三定律可知，作用于块体单元1的法向接触力为-f_n；

4-4)分别计算由法向接触力引起的，作用于块体单元1、2的力矩：M_j＝(f_n)_j×(n_cent-A)_j，j＝1,2；其中，(f_n)_j为块体单元1、2所受的法向接触力；(n_cent-A)_j为由点A指向块体单元1、2的形心的矢量；

4-5)重复4-1)到4-4)，计算块体单元2在块体单元1所有子块内的法向接触力及力矩并求和，得到由单元2侵入目标单元1，所引起的作用于块体单元1、2的法向接触力及力矩；

5)以块体单元2为目标块体，以块体单元1为接触块体，重复1)至4)，分别求出由块体单元1侵入块体单元2，所引起的作用于块体单元1以及块体单元2的法向接触力及力矩；并与4)所求的法向接触力和力矩求和，得到当前时间步，两个彼此接触的块体单元1、2的法向接触力(f_n)₁、(f_n)₂，以及其引起的力矩(M_n)₁、(M_n)₂；

6)计算当前时间步，块体单元1、块体单元2的切向接触力及其引起的力矩，具体过程为：

6-1)考虑到重叠区域很小，定义重叠区域任一边界中点B为此边界切向力作用点；

6-2)计算切向位移增量：Δδ＝(Δv·n_s)n_s*Δt；其中，n_s是切向单位向量，与边界法向量垂直；Δν是块体单元1相对于块体单元2的速度，计算公式为：

Δv＝v_c1-v_c2+ω₁×r₁-ω₂×r₂，

其中ν_c1是块体单元1形心的线速度，ν_c2是块体单元2形心的线速度，ω₁是块体单元1的角速度，ω₂是块体单元2的角速度，r₁、r₂分别是点B到块体单元1、2形心的矢量；

6-3)计算切向接触力：f_s＝f_s'+K_sΔδ，其中，K_s为切向刚度系数；f_s'为上一时间步目标块体所受切向接触力；Δδ为目标块体切向位移增量；

同时，切向接触力需满足库伦摩擦定律：若切向接触力f_s>(F_s)_max时，令f_s＝(F_s)_max；其中，是最大静摩擦角；c为凝聚力，F_n为单元的法向接触力；

则作用于块体单元1的切向接触力为f_s，根据牛顿第三定律，作用于块体单元2的切向接触力为-f_s。

6-4)由切向接触力引起的力矩：(M_s)_j＝(f_s)_j×(n_cent-B)_j，j＝1,2；

其中，(f_s)_j为块体单元1、2的切向接触力；(n_cent-B)_j为由点B指向块体单元1、2的形心的矢量；

7)根据5)和6)的计算结果，分别计算当前时间步块体单元1、2所受的力矩：(M_j)＝(M_n)_j+(M_s)_j，j＝1,2。

进一步的，在所述步骤4.4中，更新所有与目标块体相接触的块体单元作用于目标块体的接触力和力矩，具体过程如下：

1)重复步骤4.3，循环与目标块体相接触的每个块体单元作用于该块体单元的切向接触力、法向接触力和力矩；

2)当前时间步，目标块体所受的切向接触力、法相接触力和力矩分别为：

${\left(F_n\right)}_t=\underset{m}{\Σ}\left(f_n\right)$

${\left(F_s\right)}_t=\underset{m}{\Σ}\left(f_s\right)$

${\left(M\right)}_t=\underset{m}{\Σ}\left(M_i\right)$

其中t为当前时间，m为与目标块体相接触的块体单元个数，(F_n)_t为当前时间目标块体所受的切向接触力，(F_s)_t为当前时间目标块体所受的法相接触力，(M)_t为当前时间目标块体所受的力矩。

进一步的，在所述步骤六中，当前时间步的目标块体加速度计算公式为：

m_k(a_k)_t＝F_k

I_k(α_k)_t＝M_k

其中，I_k是块体单元k的主惯性矩，(a_k)_t是块体单元k当前时间步的加速度，(α_k)_t是块体单元k当前时间步的角加速度，F_k是块体单元k当前时间步所受的合力，M_k是块体单元k当前时间步所受的力矩，k＝1,2,……,N，N为块体单元个数；

根据求解得到的加速度，采用velocity verlet算法计算块体单元下一时间步的速度和位移：

$v(t+\frac{1}{2}\mathrm{\Δ}t)=v\left(t\right)+\frac{1}{2}a\left(t\right)\·\mathrm{\Δ}t$

$v(t+\mathrm{\Δ}t)=v(t+\frac{1}{2}\mathrm{\Δ}t)+\frac{1}{2}a(t+\mathrm{\Δ}t)\·\mathrm{\Δ}t$

$r(t+\mathrm{\Δ}t)=r\left(t\right)+v\left(t\right)\·\mathrm{\Δ}t+\frac{1}{2}a\left(t\right)\·{\mathrm{\Δt}}^2$

其中，Δt是时间步长，v是块体单元的速度，r是块体单元的位移。

进一步的，所述步骤七中，更新每个块体单元的顶点、形心以及内点的坐标，公式为：

x_k(t+Δt)＝x_k(t)+r_k(t+Δt)

y_k(t+Δt)＝y_k(t)+r_k(t+Δt)

其中，x_k、y_k是块体单元k的顶点、形心以及内点的坐标，r_k为块体单元k的位移,k＝1,2,……,N，N为块体单元个数。

与现有技术相比，本发明所达到的有益效果是：

1)本方法采用距离势函数的定义，明确了势函数的物理意义，可以实现二维非连续介质大规模任意凸多边形块体离散单元接触力计算问题，计算过程满足能量守恒；

2)本方法采用块体离散元建模，适用范围广，尤其适用于岩土矿石的运动形态研究，在实际矿场开矿工程中有实际指导意义。

附图说明

图1是NBS接触检测方法中块体单元映射至网格的示意图；

图2是块体单元分为子块体的示意图；

图3是两块体单元接触的重叠区域示意图；

图4是图3中边界67的势函数分布示意图；

图5是本发明实施例一矿石-边界的离散元模型示意图；

图6-图10是实施例一矿石的运动过程示意图；

图11是实施例一矿石块体单元运动过程中动能变化图；

图12是实施例二中矿石-漏斗的离散单元模型示意图；

图13-图19是实施例二放矿时的矿石运动过程示意图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明提供了一种基于距离势函数任意凸多边形块体离散单元法，包括以下步骤：

步骤一，建立矿石-边界离散模型；基于离散单元法理论，对在一定形状边界内的多个矿石，建立矿石-边界的离散元模型，将矿石模型为多边形块体单元，多边形块体单元的参数包括多边形的质量、块体间的阻尼系数、块体单元的刚度系数；

步骤二，确定模型采集矿石信息离散单元计算时间步长；

所述迭代时步采用多边形块体的基本运动方程为：

$m\stackrel{\·\·}{r\left(t\right)}+cr\left(t\right)+k\stackrel{\·}{r\left(t\right)}=F\left(t\right)$

其中：m是块体单元质量，r(t)是位移，t是时间，c是阻尼系数，k是刚度系统，F(t)是块体单元收到的外力。要使求解稳定，模拟采集矿石信息离散单元计算时间步长须满足：

$\mathrm{\Δ}t\≤2\·\sqrt{\frac{m}{k}}\·(\sqrt{1+{\mathrm{\ξ}}^2}-\mathrm{\ξ})$

其中：ξ是系统的阻尼比，

通常为了保证稳定性，一般取迭代时间步长为0.1Δt。

步骤三，在当前时间步，采用NBS接触检测方法，对所有块体单元进行接触检测，得到每个块体单元及与之接触的块体单元；

Munjiza提出的NBS(No Binary Search)接触检测法，首先将单元映射到规则格子中，然后在相邻格间进行接触对判断，NBS对以上建立离散元模型进行接触检测的具体步骤包括：

步骤3.1，计算每个块体单元形心坐标和此形心到其各顶点之间的距离l；

步骤3.2，选取所有块体单元中形心到顶点间的最大距离为l_max，记d＝2l_max，将离散模型计算区域划分为以d为边长的正方形网格；

步骤3.3，根据公式：

$x_k=1+\mathrm{int}(\frac{x_k-x_{min}}{l_{max}}+0.5)$

$y_k=1+\mathrm{int}(\frac{y_k-y_{\min }}{l_{max}}+0.5)$

将所有块体单元映射到网格上，其中，x_k、y_k为块体单元k的形心坐标，k＝1,2,……,N，N为块体单元个数，int()函数为将数值向下取整数，x_min、y_min为计算区域在x、y方向的最小值；

将所有块体单元映射到网格上，如图1所示；

步骤3.4，对每个块体单元与其周边所有相邻网格内块体单元进行接触检测，当循环到某个块体单元时，即以此块体单元所在格子为中心，进行判断此块体单元是否与该格子内部以及相邻格子内部单元接触，判断此块体单元与相邻的一个块体单元是否接触的具体过程为：

1)分别计算两个块体单元的形心之间的距离dist、两块体单元各顶点与其对应的形心之间的最大距离分别为l_1max、l_2max，计算l_1max、l_2max的平均值为l_avmax；

2)比较dist与l_avmax的大小，如果dist≤2l_avmax，则判断此两个块体单元接触，如果dist>2l_avmax，则判断此两个块体单元不接触；

步骤3.5，遍历所有块体单元，记录每个块体单元及与之接触的块体单元；

步骤四，根据步骤三的检测结果，对每个块体单元及与之接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体(Target)，与此块体单元相接触的块体单元为接触块体(Contactor)，基于所提出的块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；

具体步骤如下：

步骤4.1，确定目标块体各顶点对应的内点，其中内点定义为在多边形顶点角平分线上至顶点所在两条边距离为h的点为多边形顶点所对应的内点，其中h为多边形最大内切圆直径：

1)根据最小二分法，取h_min＝0，h_max＝2l_max迭代计算，其中l_max为目标块体的形心到顶点间的最大距离，确定目标块体最大内切圆直径h；

2)根据点到直线的距离公式：求出各顶点对应的内点坐标；其中，A、B、C是顶点所在两条边中任一条直线方程参数；x₀、y₀是顶点对应的内点坐标；

步骤4.2，根据步骤4.1所求目标块体的内点，连接顶点及其对应的内点将目标块体划分为n个子块体，n为目标块体的顶点数；图2以任意四边形单元为例，通过内点与顶点的连线，将四边形块体分成4个子块体；

步骤4.3，计算由相接触的块体单元作用于目标块体的接触力和力矩，以相互接触的两个块体单元为例，以块体单元1为目标块体(Target)，块体单元2为接触块体(Contactor)，定义距离势函数，计算由接触块体(块体单元2)所引起、作用于目标块体(块体单元1)的法向接触力、切向接触力以及力矩：

1)计算块体单元2各底边与块体单元1各子块边线的交点，确定块体单元2与块体单元1的重叠区域；如图3所示，块体单元2与块体单元1的交点分别为p₁、p₂、p₃、p₄和p₅；

2)定义距离势函数：其中，是块体单元内部点p的势函数值，h_p-12、h_p-23……h_p-n1是点p至块体单元各子块底边的距离，如h_p-12是点p至块体单元的12边的距离；n为块体单元的顶点数；min()函数为取最小值，h为块体单元最大内切圆半径；计算出各交点在块体单元1子块中的距离势函数值，若交点在块体单元1底边上则势为0，在块体单元1内点与顶点连线上则势为1；如图3所示，p₁、p₂和p₄各交点的势由势函数定义计算，p₃、p₅点因在块体单元1的边上则其势为0；

3)根据距离势函数值在子块体内呈线性分布的特征，分别计算各子块的重叠区域边界势函数值：其中，为边界的势函数值，i＝1，……，M-1，M为交点以及块体单元2在块体单元1内的顶点个数之和，为交点p_i的势函数值，为交点p_i+1的势函数值，l_pipi+1为相邻两交点p_i和p_i+1之间的距离；如图4所示，p₁、p₂和p₃交点的边界67势函数为相应的梯形的面积；

4)计算当前时间步，块体单元1内子块的法向接触力及力矩，具体计算过程为：

4-1)基于重叠区域边界势函数法向接触力计算公式为：其中，f_n是重叠区域边界上的法向接触力；n_Γ是边界的外法向向量，即法向接触力方向矢量；p_n为罚函数，为边界上各点的势函数值；根据3)所求的重叠区域边界势函数值，法向接触力计算公式简化为计算作用于重叠区域边界的法向接触力；

4-2)力的作用点取重叠区域边界势函数分布图形心在对应边界上的投影点A；

4-3)由块体单元2侵入块体单元1，所引起的作用于块体单元2的法向接触力为f_n；根据牛顿第三定律可知，作用于块体单元1的法向接触力为-f_n；

4-4)分别计算由法向接触力引起的，作用于块体单元1、2的力矩：M_j＝(f_n)_j×(n_cent-A)_j，j＝1,2；其中，(f_n)_j为块体单元1、2所受的法向接触力；(n_cent-A)_j为由点A指向块体单元1、2的形心的矢量；

4-5)重复4-1)到4-4)，计算块体单元2在块体单元1所有子块内的法向接触力及力矩并求和，得到由单元2侵入目标单元1，所引起的作用于块体单元1、2的法向接触力及力矩；

5)以块体单元2为目标块体，以块体单元1为接触块体，重复1)至4)，分别求出由块体单元1侵入块体单元2，所引起的作用于块体单元1以及块体单元2的法向接触力及力矩；并与4)所求的法向接触力和力矩求和，得到当前时间步，两个彼此接触的块体单元1、2的法向接触力(f_n)₁、(f_n)₂，以及其引起的力矩(M_n)₁、(M_n)₂；

6)计算当前时间步，块体单元1、块体单元2的切向接触力及其引起的力矩，具体过程为：

6-1)考虑到重叠区域很小，定义重叠区域任一边界中点B为此边界切向力作用点；

6-2)计算切向位移增量：Δδ＝(Δv·n_s)n_s*Δt；其中，n_s是切向单位向量，与边界法向量垂直；Δν是块体单元1相对于块体单元2的速度，计算公式为：

Δv＝v_c1-v_c2+ω₁×r₁-ω₂×r₂，

其中ν_c1是块体单元1形心的线速度，ν_c2是块体单元2形心的线速度，ω₁是块体单元1的角速度，ω₂是块体单元2的角速度，r₁、r₂分别是点B到块体单元1、2形心的矢量；

6-3)计算切向接触力：f_s＝f_s'+K_sΔδ，其中，K_s为切向刚度系数；f_s'为上一时间步目标块体所受切向接触力；Δδ为目标块体切向位移增量；

同时，切向接触力需满足库伦摩擦定律：若切向接触力f_s>(F_s)_max时，令f_s＝(F_s)_max；其中，是最大静摩擦角；c为凝聚力，F_n为单元的法向接触力；

则作用于块体单元1的切向接触力为f_s，根据牛顿第三定律，作用于块体单元2的切向接触力为-f_s。

6-4)由切向接触力引起的力矩：(M_s)_j＝(f_s)_j×(n_cent-B)_j，j＝1,2；

其中，(f_s)_j为块体单元1、2的切向接触力；(n_cent-B)_j为由点B指向块体单元1、2的形心的矢量；

7)根据5)和6)的计算结果，分别计算当前时间步块体单元1、2所受的力矩：(M_j)＝(M_n)_j+(M_s)_j，j＝1,2。

步骤4.4，更新所有与目标块体相接触的块体单元作用于目标块体的接触力和力矩，具体过程如下：

1)重复步骤4.3，循环与目标块体相接触的每个块体单元作用于该块体单元的接触力和力矩；

2)当前时间步，目标块体所受的切向接触力、法相接触力和力矩分别为：

${\left(F_n\right)}_t=\underset{m}{\Σ}\left(f_n\right)$

${\left(F_s\right)}_t=\underset{m}{\Σ}\left(f_s\right)$

${\left(M\right)}_t=\underset{m}{\Σ}\left(M_i\right)$

其中t为当前时间，m为与目标块体相接触的块体单元个数，(F_n)_t为当前时间目标块体所受的切向接触力，(F_s)_t为当前时间目标块体所受的法相接触力，(M)_t为当前时间目标块体所受的力矩。

步骤五，更新当前时间步所有块体单元所受的力和力矩；

步骤六，根据牛顿第二定律，计算出块体单元当前时间步的加速度，在由velocity verlet算法计算出块体单元下一时间步的速度和位移；

块体单元k当前时间步的目标块体加速度计算公式为：

m_k(a_k)_t＝F_k

I_k(α_k)_t＝M_k

其中，I_k是块体单元k的主惯性矩，(a_k)_t是块体单元k当前时间步的加速度，(α_k)_t是块体单元k当前时间步的角加速度，F_k是块体单元k当前时间步所受的合力，M_k是块体单元k当前时间步所受的力矩，k＝1,2,……,N，N为块体单元个数；

根据求解得到的加速度，采用velocity verlet算法求解单元下一时间步的的速度和位移，velocity verlet算法由块体单元在t时刻的速度和加速度，计算出t+Δt时刻的位置，公式为：

$v(t+\frac{1}{2}\mathrm{\Δ}t)=v\left(t\right)+\frac{1}{2}a\left(t\right)\·\mathrm{\Δ}t$

$v(t+\mathrm{\Δ}t)=v(t+\frac{1}{2}\mathrm{\Δ}t)+\frac{1}{2}a(t+\mathrm{\Δ}t)\·\mathrm{\Δ}t$

$r(t+\mathrm{\Δ}t)=r\left(t\right)+v\left(t\right)\·\mathrm{\Δ}t+\frac{1}{2}a\left(t\right)\·{\mathrm{\Δt}}^2$

其中，Δt是时间步长，v是块体的速度，r是块体的位移；

步骤七，根据步骤六中块体单元的位移，更新所有块体单元顶点、形心及内点的坐标，完成当前时间步的计算；

更新每个块体单元的顶点、形心以及内点的坐标，公式为：

x_k(t+Δt)＝x_k(t)+r_k(t+Δt)

y_k(t+Δt)＝y_k(t)+r_k(t+Δt)

其中，x_k、y_k是块体单元k的顶点、形心以及内点的坐标，r_k为块体单元k的位移,k＝1,2,……,N，N为块体单元个数。

步骤八，重复步骤三至七计算下一时间步，直至计算完所有时间步。

实施例一

矿石以一定速度落入槽框内的运动形态，槽框为一个矩形框，即多边形矿石块体单元在一定形状边界内的运动形态，如图5所示，四个矿石块体以0.08m/s的速度向下运动，采用本发明提供的方法，为矿石-边界建立离散元模型，现将四个矿石块体分为35个块体离散单元，即计算分析35个块体单元的运动情况。

图6至10是矿石块体开始接触槽框后，矿石块体运动过程图。在离散元模型中，所有的块体单元以0.08m/s的速度向下运动，当最下部的块体单元与边界接触时，块体单元与块体单元之间、块体单元与边界之间由于接触产生接触力和力矩，块体单元所受的合力及力矩发生改变，使块体单元的运动状态发生变化。块体飞离彼此，以恒定的速度继续运动，直到与别的块体或者边框发生接触，产生新的接触力和力矩，块体的运动状态再次发生改变，以此循环直到运动结束。

由于接触，系统的部分动能转化为势能，随着块体的分离，势能再次转化为动能，但整个过程中系统的总能量是不变的。图11是块体运动过程中动能变化图。

通过观察图6-10可以看到，本发明提出的基于距离势函数的二维块体离散单元计算方法，可以实现二维非连续介质的任意凸多边形离散单元数值接触力计算，图11的系统运动过程中动能变化图可以说明本发明所提出的离散单元计算方法的正确性。

实施例二

某放矿场漏斗内充满了已崩大小不等、不规则的矿石，将漏斗闸门打开，受自重作用，矿石从漏孔中放出后的运动形态，即多边形矿石块体单元在一定形状边界内的运动形态，采用本发明提供的方法，为矿石-边界建立离散元模型，如图12所示，共划分矿体离散单元309个，漏斗及槽框离散单元162个，漏斗及槽框固定。

图12为矿石在漏斗中的初始状态，图13到图19为漏斗口闸门开启后，放矿时的几何形态。矿体之间受重力作用，彼此接触，产生接触力和力矩，矿体的速度与加速度因此发生改变，矿体以变加速度向槽框底运动，直至与槽框底边接触，考虑摩擦和阻尼的影响，矿体的速度逐渐减小为0，矿体逐渐堆积。当所有的矿体都堆积到槽框底部时，运动结束。基于本发明所提出的距离势函数的二维块体离散单元法模拟放矿过程，能清楚地描述崩落矿体的运动规律。

综上所述，采用基于距离势函数的二维离散单元法，能够解决二维任意形状的凸多边形离散单元计算问题。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (8)

1.一种基于距离势函数任意凸多边形块体离散单元法，其特征是，包括以下步骤：

步骤一，建立矿石-边界离散模型；基于离散单元法理论，对在一定形状边界内的多个矿石，建立矿石-边界的离散元模型，矿石模型为多个多边形块体单元，多边形块体单元的参数包括多边形的质量、块体间的阻尼系数及块体单元的刚度系统；

步骤二，确定模型采集矿石信息的迭代时间步长；

步骤三，在当前时间步，采用NBS接触检测方法对所有块体单元进行接触检测，得到每个块体单元及与之接触的块体单元；

步骤四，根据步骤三的检测结果，对相互接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体(Target)，与此块体单元相接触的块体单元为接触块体(Contactor)，基于所提出的块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；

步骤五，重复步骤四计算当前时间步内每个块体单元所受的接触力和力矩；

步骤六，根据牛顿第二定律，计算出每个块体单元当前时间步的加速度，再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移；

步骤七，根据步骤六中块体单元的位移，更新每个块体单元顶点、形心及内点的坐标，完成当前时间步的计算；

步骤八，重复步骤三至七计算下一时间步，直至计算完所有时间步。

2.根据权利要求1所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，在所述步骤二中，采集矿石信息离散单元计算时间步长为：

$\mathrm{\Δ}t\≤2\·\sqrt{\frac{m}{k}}\·(\sqrt{1+{\mathrm{\ξ}}^2}-\mathrm{\ξ})$

其中：ξ是系统的阻尼比，m是块体单元质量，c是阻尼系数，k是刚度系数。

3.根据权利要求1所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，所述步骤三中，NBS接触检测法进行接触检测的具体步骤包括：

步骤3.1，对每个块体单元，计算其形心坐标及各顶点到形心之间的距离l；

步骤3.2，选取形心到顶点之间的最大距离为l_max，记d＝2·l_max，将离散模型计算区域划分为以d为边长的正方形网格；

步骤3.3，根据公式：

$x_k=1+\mathrm{int}(\frac{x_k-x_{min}}{l_{max}}+0.5)$

$y_k=1+int(\frac{y_k-y_{\min }}{l_{max}}+0.5)$

将所有块体单元映射到网格上，其中，x_k、y_k为块体单元k的形心坐标，k＝1,2,……,N，N为块体单元个数，x_min、y_min为计算区域在x、y方向的最小值；

步骤3.4，对每个块体单元与其周边所有相邻网格内块体单元进行接触检测，当循环到某个块体单元时，即以此块体单元所在格子为中心，进行判断此块体单元是否与该格子内部以及相邻格子内部单元接触，判断此块体单元与相邻的一个块体单元是否接触的具体过程为：

1)分别计算两个块体单元的形心之间的距离dist、两块体单元各顶点与其对应的形心之间的最大距离分别为l_{1 max}、l_{2 max}，计算l_{1 max}、l_{2 max}的平均值为l_{av max}；

2)比较dist与l_{av max}的大小，如果dist≤2l_{av max}，则判断此两个块体单元接触，如果dist>2l_{av max}，则判断此两个块体单元不接触；

步骤3.5，遍历所有块体单元，记录每个块体单元及与之接触的块体单元。

4.根据权利要求1所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，在所述步骤四中，目标块体的法向接触力、切向接触力以及力矩具体计算步骤如下：

步骤4.1，确定目标块体各顶点对应的内点，其中内点定义为在多边形顶点角平分线上至顶点所在两条边距离为h的点为多边形顶点所对应的内点，其中h为多边形最大内切圆直径：

1)根据最小二分法，取h_min＝0，h_max＝2l_max迭代计算，其中l_max为目标块体的形心到顶点间的最大距离，确定目标块体最大内切圆直径h；

2)根据点到直线的距离公式：求出各顶点对应的内点坐标；其中，A、B、C是顶点所在两条边任一直线方程参数；x₀、y₀是顶点对应的内点坐标；

步骤4.2，根据步骤4.1所求目标块体的内点，连接顶点及其对应的内点将目标块体划分为n个子块体，n为目标块体的顶点数；

步骤4.3，定义距离势函数，循环计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩；

步骤4.4，更新所有与目标块体相接触的块体单元作用于目标块体的接触力和力矩。

5.根据权利要求4所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，在所述步骤4.3中，计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩，以相互接触的块体单元1和块体单元2说明计算过程，对此两个彼此接触的块体单元，以块体单元1为目标单元，求解由块体单元2所引起、作用于目标单元的接触力和力矩过程为：

1)计算块体单元2各底边与块体单元1各子块边线的交点，确定块体单元2与块体单元1的重叠区域；

2)定义距离势函数：其中，是块体单元内部点p的势函数值，h_p-12、h_p-23……h_p-n1是点p至块体单元各子块底边的距离，n为块体单元的顶点数；h为块体单元最大内切圆半径；计算出各交点在块体单元1子块中的距离势函数值，若交点在块体单元1底边上则势为0，在块体单元1内点与顶点连线上则势为1；

3)根据距离势函数值在子块体内呈线性分布的特征，分别计算各子块的重叠区域边界势函数值：其中，为边界的势函数值，i＝1，……，M-1，M为交点以及块体单元2在块体单元1内的顶点个数之和，为交点p_i的势函数值，为交点p_i+1的势函数值，为相邻两交点p_i和p_i+1之间的距离；

4)计算当前时间步，块体单元1内子块的法向接触力及力矩，具体计算过程为：

4-1)基于重叠区域边界势函数法向接触力计算公式为：其中，f_n是重叠区域边界上的法向接触力；n_Γ是边界的外法向向量，即法向接触力方向矢量；p_n为罚函数，为边界上各点的势函数值；根据3)所求的重叠区域边界势函数值，法向接触力计算公式简化为计算作用于重叠区域边界的法向接触力；

4-2)力的作用点取重叠区域边界势函数分布图形心在对应边界上的投影点A；

4-3)由块体单元2侵入块体单元1，所引起的作用于块体单元2的法向接触力为f_n；根据牛顿第三定律可知，作用于块体单元1的法向接触力为-f_n；

4-4)分别计算由法向接触力引起的，作用于块体单元1、2的力矩：M_j＝(f_n)_j×(n_cent-A)_j，j＝1,2；其中，(f_n)_j为块体单元1、2所受的法向接触力；(n_cent-A)_j为由点A指向块体单元1、2的形心的矢量；

4-5)重复4-1)到4-4)，计算块体单元2在块体单元1所有子块内的法向接触力及力矩并求和，得到由单元2侵入目标单元1，所引起的作用于块体单元1、2的法向接触力及力矩；

5)以块体单元2为目标块体，以块体单元1为接触块体，重复1)至4)，分别求出由块体单元1侵入块体单元2，所引起的作用于块体单元1以及块体单元2的法向接触力及力矩；并与4)所求的法向接触力和力矩求和，得到当前时间步，两个彼此接触的块体单元1、2的法向接触力(f_n)₁、(f_n)₂，以及其引起的力矩(M_n)₁、(M_n)₂；

6)计算当前时间步，块体单元1、块体单元2的切向接触力及其引起的力矩，具体过程为：

6-1)考虑到重叠区域很小，定义重叠区域任一边界中点B为此边界切向力作用点；

6-2)计算切向位移增量：Δδ＝(Δv·n_s)n_s*Δt；其中，n_s是切向单位向量，与边界法向量垂直；Δν是块体单元1相对于块体单元2的速度，计算公式为：

Δv＝v_c1-v_c2+ω₁×r₁-ω₂×r₂，

其中ν_c1是块体单元1形心的线速度，ν_c2是块体单元2形心的线速度，ω₁是块体单元1的角速度，ω₂是块体单元2的角速度，r₁、r₂分别是点B到块体单元1、2形心的矢量；

6-3)计算切向接触力：f_s＝f_s'+K_sΔδ，其中，K_s为切向刚度系数；f_s'为上一时间步目标块体所受切向接触力；Δδ为目标块体切向位移增量；

同时，切向接触力需满足库伦摩擦定律：若切向接触力f_s>(F_s)_max时，令f_s＝(F_s)_max；其中，是最大静摩擦角；c为凝聚力，F_n为单元的法向接触力；

则作用于块体单元1的切向接触力为f_s，根据牛顿第三定律，作用于块体单元2的切向接触力为-f_s。

6-4)由切向接触力引起的力矩：(M_s)_j＝(f_s)_j×(n_cent-B)_j，j＝1,2；

其中，(f_s)_j为块体单元1、2的切向接触力；(n_cent-B)_j为由点B指向块体单元1、2的形心的矢量；

7)根据5)和6)的计算结果，分别计算当前时间步块体单元1、2所受的力矩：(M_j)＝(M_n)_j+(M_s)_j，j＝1,2。

6.根据权利要求5所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，在所述步骤4.4中，更新所有与目标块体相接触的块体单元作用于目标块体的接触力和力矩，具体过程如下：

1)重复步骤4.3，循环与目标块体相接触的每个块体单元作用于该块体单元的切向接触力、法向接触力和力矩；

2)当前时间步，目标块体所受的切向接触力、法相接触力和力矩分别为：

${\left(F_n\right)}_t=\underset{m}{\Σ}\left(f_n\right)$

${\left(F_s\right)}_t=\underset{m}{\Σ}\left(f_s\right)$

${\left(M\right)}_t=\underset{m}{\Σ}\left(M_i\right)$

其中t为当前时间，m为与目标块体相接触的块体单元个数，(F_n)_t为当前时间目标块体所受的切向接触力，(F_s)_t为当前时间目标块体所受的法相接触力，(M)_t为当前时间目标块体所受的力矩。

7.根据权利要求1所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，在所述步骤六中，当前时间步的目标块体加速度计算公式为：

m_k(a_k)_t＝F_k

I_k(α_k)_t＝M_k

其中，I_k是块体单元k的主惯性矩，(a_k)_t是块体单元k当前时间步的加速度，(α_k)_t是块体单元k当前时间步的角加速度，F_k是块体单元k当前时间步所受的合力，M_k是块体单元k当前时间步所受的力矩，k＝1,2,……,N，N为块体单元个数；

根据求解得到的加速度，采用velocity verlet算法计算块体单元下一时间步的速度和位移：

$v(t+\frac{1}{2}\mathrm{\Δ}t)=v\left(t\right)+\frac{1}{2}a\left(t\right)\·\mathrm{\Δ}t$

$v(t+\mathrm{\Δ}t)=v(t+\frac{1}{2}\mathrm{\Δ}t)+\frac{1}{2}a(t+\mathrm{\Δ}t)\·\mathrm{\Δ}t$

$r(t+\mathrm{\Δ}t)=r\left(t\right)+v\left(t\right)\·\mathrm{\Δ}t+\frac{1}{2}a\left(t\right)\·{\mathrm{\Δt}}^2$

其中，Δt是时间步长，v是块体单元的速度，r是块体单元的位移。

8.根据权利要求1所述的一种基于距离势函数任意凸多边形块体离散单元法，其特征是，所述步骤七中，更新每个块体单元的顶点、形心以及内点的坐标，公式为：

x_k(t+Δt)＝x_k(t)+r_k(t+Δt)

y_k(t+Δt)＝y_k(t)+r_k(t+Δt)

其中，x_k、y_k是块体单元k的顶点、形心以及内点的坐标，r_k为块体单元k的位移,k＝1,2,……,N，N为块体单元个数。