一种基于距离势函数三维任意凸多边形块体离散单元法的制作方法

技术特征：

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

步骤一，将研究对象离散为多面体块体单元体系，并根据研究对象的大小，确定计算区域；

步骤二，确定计算时间步长；

步骤三，在当前时间步，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测，并依据块体单元形心之间的距离、块体单元各顶点与其对应的形心之间的最大距离确定每个块体单元及与之接触的块体单元；

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

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

步骤六，根据刚体运动方程，计算出每个块体单元当前时间步的加速度，再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移；

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

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

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

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

其中，Δt是计算时间步长；ξ是系统的阻尼比，m是块体单元质量，c是阻尼系数，k是刚度系数。

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

步骤3.1，对每个块体单元，计算其最大外接球半径；

步骤3.2，所有块体单元根据其最大外接球半径，按从大到小的顺序分成n组，具体为：

1)第1组中块体单元的最大外接球半径为d₁，d₁的取值范围为：D≤d₁＜D/α；

2)第2组中块体单元的最大外接球半径为d₂，d₂的取值范围为：D/α≤d₂＜D/α²；

3)第3组中块体单元的最大外接球半径为d₃，d₃的取值范围为：D/α²≤d₃＜D/α³；

4)以此类推，第n-1组中块体单元的最大外接球半径为d_n-1，d_n-1的取值范围为：D/α^n-2≤d_n-1＜D/α^n-1；

5)第n组中块体单元的最大外接球半径为d_n，则d_n的取值范围为：d_n≥D/α^n-1；

其中，n＝int[log(D/d)/logα+0.416]+1,D为所有块体单元中最大块体单元外接球半径，d为所有块体单元中最小块体单元外接球半径，α＝2；

步骤3.3，对第1组的块体单元与所有块体单元进行接触检测，具体步骤包括：

1)取l_max＝2d₁，将计算区域划分为以l_max为边长的正方体网格；

2)根据公式：

$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)$

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

将所有块体单元映射到网格上，其中，(x_K,y_K,z_K)为第K个块体单元的形心坐标，K＝1,2,…,N，N为块体单元个数，x_min、y_min、z_min分别为计算区域在x、y、z方向的最小值；

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

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

3-2)比较dist与l_avmax，若dist≤2l_avmax，则判断此两个块体单元接触，否则判断此两个块体单元不接触；

步骤3.4，将第一组的块体单元从所有块体单元之中去除，对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测，具体步骤包括：

1)取l_max＝2d₂，将计算区域划分为以l_max为边长的正方体网格；

2)根据步骤3.3中2)至3)的方法，将第2组到第n组的块体单元映射到网格上，并对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测；

步骤3.5，以此类推，将已进行接触检测的块体单元组从所有块体单元中去除，并进行下一组块体单元与剩余所有块体单元组的接触检测，直至完成n组块体的接触检测。

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

步骤4.1，确定目标块体单元的最大内切球直径；

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

步骤4.3，计算所有与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩。

5.根据权利要求4所述的一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，步骤4.2中，计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩的具体步骤为：

1)设块体单元β_t与块体单元β_c为两个相接触的块体单元，其中，块体单元β_c为接触单元，β_t为目标单元；

2)定义距离势函数：其中，表示β_t内部点p的势函数值，h_p-1,h_p-2,…,h_p-M分别表示点p至β_t各底面的距离，M为β_t的底面数，h为β_t的最大内切球半径；

3)根据势函数法计算法向接触力，计算法向接触力，计算公式为：

其中，f_n是法向接触力；分别为β_c、β_t重叠区域内的点分别在β_t与β_c的距离势函数梯度；V为β_t与β_c的重叠区域；

应用高斯公式，将法向接触力体积分计算公式转化为重叠区域边界面的积分：

其中，n_τ为重叠区域边界面的外法向矢量，即为边界面上法向接触力方向矢量；S为重叠区域的边界面；

4)应用3)中的法向接触力，计算公式计算当前时间步β_c嵌入β_t引起的法向接触力及力矩，具体计算过程为：

4-1)以β_t各顶点的角平分线为界、其内切圆半径为控制距离，将β_t分割为若干子多面体；分别以β_c各底面切割β_t的各子多面体，确定β_t与β_c重叠区域的边界面；由此，β_t与β_c重叠区域法向接触力的计算转化为β_c各底面与β_t的各子多面体之间的相互作用；

4-2)根据距离势函数定义，距离势函数在子多面体内线性分布，在β_c的底面上建立局部坐标系，定义该底面上的距离势函数计算公式为：

其中，(x_c,y_c)为该底面上点的局部坐标；a、b、c为距离势函数方程系数；

4-3)将β_c中一底面在β_t各子多面体内的重叠面的任意一顶点作为固定点，其余各顶点与固定点的连线为分界线，将重叠面分割为若干三角形，由2)距离势函数计算公式，分别计算各三角形顶点在β_t内的距离势函数值；并根据各三角形顶点在β_t内的距离势函数值，联立方程组求得a、b、c；

4-4)分别在4-3)中的若干三角形上建立局部坐标系，并通过坐标变换，将4-2)的底面坐标系变换到三角形局部坐标系上，坐标变换的形函数分别为：N(1)＝1-ζ-η，N(2)＝ζ，N(3)＝η，ζ,η∈(0,1)；(ζ,η)为三角形上的点在三角形局部坐标系内的坐标；

根据3)中法向接触力积分公式，计算作用于三角形面的法向接触力：

其中，f_n,tri为三角形面上的法向接触力；|J|为积分坐标变换的雅可比行列式；k_n为法向刚度；

将4-2)中距离势函数计算公式代入，积分化简得到作用于三角形面的法向接触力：

$f_{n,tri}=n_{\mathrm{\τ}}{k}_n(\frac{A}{6}+\frac{B}{6}+\frac{C}{2})$

其中A＝[a(x₂-x₁)+b(y₂-y₁)]|J|,B＝[a(x₃-x₁)+b(y₃-y₁)]|J|,C＝[ax₁+by₁+c]|J|；(x₁,y₁)、(x₂,y₂)、(x₃,y₃)分别为三角形三个顶点在底面局部坐标系中的坐标值；

4-4)由作用于三角形面的法向接触力引起的接触力矩为：

将4-2)中距离势函数计算公式代入简化得到作用于三角形面的法向接触力引起的力矩为：

$M_{\mathrm{\ζ}}=n_{\mathrm{\τ}}{k}_n(\frac{A}{24}+\frac{B}{12}+\frac{C}{6})$

$M_{\mathrm{\η}}=n_{\mathrm{\τ}}{k}_n(\frac{A}{12}+\frac{B}{24}+\frac{C}{6})$

其中，M_ζ、M_η分别为作用于三角形面的法向接触力对三角形局部坐标系ζ轴、η轴的力矩；

4-5)重复4-3)到4-4)，并对计算得到的法向接触力和力矩分别进行矢量求和，得到β_c的一个底面与β_t间的法向接触力及力矩；

4-6)重复4-2)到4-5)，分别计算β_c所有底面与β_t的法向接触力和力矩，并分别对计算得到的法向接触力和力矩进行矢量求和，得到由β_c嵌入β_t所引起的法向接触力及合力矩；

5)以块体单元β_c为目标单元、β_t为接触单元，求出由β_t嵌入β_c所引起的法向接触力及力矩，并与4)所求由β_c嵌入β_t所引起的法向接触力和力矩分别进行矢量求和，得到当前时间步β_c与β_t间的法向接触力F_n以及其引起的合力矩M_n；

6)计算当前时间步，作用于β_c、β_t的切向接触力及其引起的力矩，具体过程为：

6-1)计算β_t相对于β_c的切向位移增量：Δδ_s＝(Δv·n_s)n_s·Δt；其中，n_s是切向单位向量，与当前时间步法向接触力方向垂直；Δν是β_t相对于β_c的速度；

6-2)计算β_t的切向接触力增量：Δf_s＝k_s·Δδ_s，其中，k_s为切向刚度系数；

6-3)将上一时间步的切向接触力转换到当前时间步切向接触力增量方向：

f_s'＝r×f_s”

其中，f_s'为方向转换后上一时间步切向接触力；f_s”为方向转换前上一时间步切向接触力；r为旋转矩阵；

6-4)分别计算当前时间步，β_t的切向接触力，计算公式为：

F_s＝f_s'+Δf_s

同时，切向接触力需满足库伦摩擦定律：其中，为最大静摩擦角；c_凝为凝聚力；F_n为法向接触力；(F_s)_max为允许的最大切向接触力，即最大静摩擦力；若切向接触力F_s＞(F_s)_max时，令F_s＝(F_s)_max；

6-5)计算由切向接触力引起的作用于β_t的力矩：其中，F_s为β_t所受切向接触力；为由切向接触力作用点P_pτ指向β_t形心的矢量；

7)根据5)和6)的计算结果，计算当前时间步β_t受的力矩，计算公式为：

M＝M_n+M_s

8)同理，根据步骤6)和7)的计算公式，计算当前时间步β_c所受的力矩。

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

1)重复步骤4.2，循环与目标块体可能相接触的每个块体单元，若两块体没有重叠区域则判断两单元不接触，若有重叠区域则求出作用于该块体单元的切向接触力、法向接触力和力矩；

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

${\left(F_n\right)}_t=\underset{m^{\′}}{\mathrm{\Σ}}{F}_n$

${\left(F_s\right)}_t=\underset{m^{\′}}{\mathrm{\Σ}}{F}_s$

${\left(M\right)}_t=\underset{m^{\′}}{\mathrm{\Σ}}M$

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

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

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

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

z_t(t+Δt)＝z_t(t)+(r(t+Δt))_z

其中，(x_t,y_t,z_t)是块体单元的顶点、形心的坐标，(r(t+Δt))_x、(r(t+Δt))_y、(r(t+Δt))_z分别为块体单元的位移在x_t,y_t,z_t方向的分量。