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

权利要求书

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

步骤一：建立可变形离散单元系统，所述可变形离散单元系统包含若干个离散单元以及将离散单元剖分网格后形成的有限单元；

步骤二：确定系统的时间步长Δt；

步骤三：对所有离散单元外围一层的网格单元进行接触检测确定目标单元和接触单元，将离散单元外围一层产生接触的网格单元定义为目标单元，与之接触的离散单元外围一层的网格单元则定义为接触单元，计算当前时间步作用于目标单元和接触单元的法向接触力、切向接触力：

步骤四：将步骤三计算得出作用在接触单元以及目标单元上的接触力的合力用形函数转化成网格单元当前时刻荷载的等效节点力矢量表达式如下：

其中，是网格单元当前刻荷载载荷的等效节点力矢量，和分别是当前时刻的网格单元体力和面力的载荷矢量，N是网格单元节点的形函数，V₀是网格单元的体积，A_0t是当前t时刻网格单元的表面积，A₀是网格单元的表面积；

步骤五：由步骤四计算得出的载荷的等效节点力矢量，建立系统增量形式的动力控制方程并求解得到下一时刻t+Δt有限单元的位移；

步骤六：根据步骤五有限单元的位移，更新每个网格单元顶点的坐标的信息，完成当前时刻的计算；

步骤七：重复步骤二至步骤六计算，直至计算完所有时间步。

2.根据权利要求1所述的一种基于距离势函数二维可变形凸多边形块体离散单元法，其特征在于，所述时间步长Δt须满足：

Δt＝min(Δt_D,Δt_s)，

Δt_s≤L/C，

其中，Δt_D为离散单元的计算时间步长；ξ为离散单元的阻尼比且m为离散单元块体质量，c为阻尼系数，k为刚度系数，Δt_s为有限单元的时间步长，L为所有有限单元的最小边长，C的取值为10000。

3.根据权利要求1所述的一种基于距离势函数二维可变形凸多边形块体离散单元法，其特征在于步骤三中，计算当前时刻每个接触单元与目标单元的接触力时采用基于距离势函数的离散单元法的接触力算法。

4.根据权利要求1所述的一种基于距离势函数二维可变形凸多边形块体离散单元法，其特征在于，在所述步骤五中，建立系统增量形式的运动控制方程，并求解得到下一时刻t+Δt有限单元的位移，具体包括：

步骤5.1：建立参考坐标系，选择变形前的构形为参考构形；

步骤5.2：采用广义Newmark法进行时间域离散，预测当前时刻力学量；

步骤5.3：建立系统增量形式的动力控制方程求解得到当前时刻t到下一时刻t+Δt的加速度增量其中，M是有限单元的质量矩阵，D是有限单元的阻尼矩阵，K是有限单元的刚度矩阵，是有限单元的加速度增量，是有限单元的速度增量，Δu是有限单元的位移增量；

步骤5.4：再由广义Newmark法进行时间域离散，计算出每个有限单元下一时刻t+Δt的位移，由于网格单元与有限单元的节点坐标相一致，因此计算得出下一时刻t+Δt每个网格单元的位移。

5.根据权利要求1所述的一种基于距离势函数二维可变形凸多边形块体离散单元法，其特征在于，所述步骤五中，更新每个有限单元的顶点的坐标公式为：

x(t+Δt)＝x(t)+(r(t+Δt))_x

y(t+Δt)＝y(t)+(r(t+Δt))_y

其中，x(t+Δt)、y(t+Δt)是当前时间步网格单元的顶点的坐标；x(t)、y(t)是上一时间步块体单元的顶点的坐标；(r(t+Δt))_x、(r(t+Δt))_y分别为网格单元的位移在x,y方向的分量。

6.根据权利要求1所述的一种基于距离势函数二维可变形凸多边形块体离散单元法，其特征在于，离散单元的参数包括离散单元的节点坐标、质量、阻尼比、刚度；有限单元的参数包括有限单元的节点坐标、质量矩阵、阻尼矩阵、刚度矩阵。