[发明专利]一种基于距离势函数任意凸多边形块体离散单元法

摘要

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

主权项

一种基于距离势函数任意凸多边形块体离散单元法，其特征是，包括以下步骤：步骤一，建立矿石‑边界离散模型；基于离散单元法理论，对在一定形状边界内的多个矿石，建立矿石‑边界的离散元模型，矿石模型为多个多边形块体单元，多边形块体单元的参数包括多边形的质量、块体间的阻尼系数及块体单元的刚度系数；步骤二，确定模型采集矿石信息的迭代时间步长；步骤三，在当前时间步，采用NBS接触检测方法对所有块体单元进行接触检测，得到每个块体单元及与之接触的块体单元；步骤四，根据步骤三的检测结果，对相互接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体(Target)，与此块体单元相接触的块体单元为接触块体(Contactor)，基于所提出的块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；定义距离势函数：其中，是块体单元内部点p的势函数值，h_p‑12、h_p‑23……h_p‑n1是点p至块体单元各子块底边的距离，n为块体单元的顶点数；h为块体单元最大内切圆半径；在所述步骤四中，目标块体的法向接触力、切向接触力以及力矩具体计算步骤如下：步骤4.1，确定目标块体各顶点对应的内点，其中内点定义为在多边形顶点角平分线上至顶点所在两条边距离为h的点为多边形顶点所对应的内点，其中h为多边形最大内切圆半径：1)根据最小二分法，取h_min＝0，h_max＝2l_max迭代计算，其中l_max为目标块体的形心到顶点间的最大距离，确定目标块体最大内切圆半径h；2)根据点到直线的距离公式：求出各顶点对应的内点坐标；其中，d＝h，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至块体单元各子块底边的距离，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的切向接触力