[发明专利]一种可变形二维任意圆化凸多边形离散单元法

权利要求书

1.一种可变形二维任意圆化凸多边形离散单元法，其特征在于，包括以下步骤：

(1)首先建立离散的多边形块体体系，根据研究对象的大小确定计算区域，再将离散的多边形块体单元划分成有限元网格，划分的网格单元用于可变形有限元的计算；

(2)确定时间步长；

(3)在当前时间步，确定接触单元，采用NBS接触检测方法对所有块体离散单元外围一层单元进行接触检测，得到每个块体外围一层的接触单元以及与之可能接触的单元，并且隶属于同一个离散单元的接触单元不进行接触检测；

(4)根据步(3)的接触检测结果，对可能相互接触的接触单元进行接触力计算，得到当前时间步每个接触单元所受的接触力；该步骤的具体过程如下：

(41)首先定义圆化凸多边形离散单元：设H是多边形，B是半径为R的圆，则将H和B的闵可夫斯基和记为P，P即为圆化凸多边形离散单元，H是P的骨架；圆化半径R＝hc，h为多边形H的最大内切圆半径，c为圆化半径的系数；

(42)确定圆化多边形离散单元的骨架间的最小距离：

(42-1)确定二维情况下多边形有两种接触方式：点-点接触和点-线接触；

(42-2)计算两种接触方式下骨架间的最小距离；

(43)基于嵌入深度的离散元计算公式，计算当前时间步接触力的法向以及法向接触力：

(43-1)计算接触力的法向：

其中，是骨架H₁、H₂之间的接触力法向，是骨架H₁、H₂之间最小距离对应的点坐标，为骨架H₁、H₂之间的最小距离；

(43-2)计算法向接触力：

其中，是骨架H₁、H₂之间的法向接触力，K_n是离散单元的法向刚度，δ(H₁,H₂)是两个圆化凸多边形离散单元P₁、P₂之间的重叠距离，R₁、R₂分别为P₁、P₂的圆化半径；

(43-3)对于两个相互接触的圆化凸多边形离散单元P₁、P₂，先以P₁为目标单元、P₂为接触单元，按照步骤(43-1)-(43-2)求出由P₂嵌入P₁所引起的法向接触力，再以P₂为目标单元、P₁为接触单元，求出由P₁嵌入P₂所引起的法向接触力，两次求出的法向接触力的矢量和即为当前时间步P₁与P₂间的法向接触力F_n；

(44)计算当前时间步离散单元间的切向接触力：

F_s＝f_s'+Δf_s

其中，f_s'为上一时间步的切向接触力，Δf_s为切向接触力增量，Δf_s＝k_s·Δδ_s，k_s为切向刚度系数，Δδ_s为切向位移增量，Δδ_s＝(Δv·n_s)n_s·Δt，n_s是切向单位向量，与当前时间步接触力法向垂直，Δv是离散单元间的相对速度；

同时，当切向接触力F_s大于最大静摩擦力(F_s)_max时，令F_s＝(F_s)_max，为最大静摩擦角，c_凝为凝聚力；

(45)计算当前时间步的接触力F＝F_n+F_s；

(5)将步骤(4)计算得出的接触力以及作用在圆化凸多边形离散单元系统上的外力用形函数转化成载荷的等效节点力矢量；

(6)由步骤(5)计算得出的载荷的等效节点力矢量以及系统应力场的等效节点矢量，求解动力控制方程，得到当前时间步系统的各变量值，其中包括块体单元的位移；

(7)根据步骤(6)中得到的块体单元位移更新块体单元的几何信息，几何信息包含每个块体单元顶点和形心的坐标完成当前时间步的计算；

(8)重复步骤(3)-(7)计算下一时间步，直至计算完所有时间步。