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

Abstract

本发明公开了一种基于距离势函数三维任意凸多边形块体离散单元法，将研究对象离散为块体单元体系；根据本发明提供的非均匀块体离散单元接触检测方法，对所有块体单元进行接触检测；定义距离势函数，对相互接触的块体单元进行法向接触力计算，得到作用于块体单元的法向接触力及其引起的力矩；基于本发明提供的切向接触力方向矢量计算方法，根据切向相对位移的增量，计算得到作用于块体单元的切向接触力及其引起的力矩；由刚体运动方程求解块体单元的加速度与角加速度；采用velocity verlet算法计算块体单元的速度与位移，根据各块体单元的运动和相对位置，描述整个介质的变形和运动状态。

Description

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

技术领域

本发明涉及一种块体离散单元法，具体涉及一种基于距离势函数三维任意凸多边形块体离散单元法，属于块体离散元技术领域。

背景技术

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

离散单元法原理直观简单，以每个块体的运动为基础，构建整个系统的运动情况。能够模拟研究物体在动态条件下的变形和破坏过程，特别是大变形以及断裂等。目前离散单元法已经在散体力学、岩土工程和地质工程等领域得到广泛应用并取得重要成果。

目前，离散单元法主要分为两种。一是由Cundall教授提出传统离散单元法。该方法在计算中可能存在能量不守恒的问题；接触检测判断需考虑角-角接触，角-面接触，面-面接触等不同的块体接触形式，过程繁琐；无法处理角-角接触的情况，需采用圆角化进行规避，过程复杂；接触检测方面，目前普遍采用的方法是公共面法，但是公共面法仍然存在很多问题，例如公共面的正确选取、唯一性以及旋转时的迭代误差等。二是由英国A.MUNJIZA教授提出有限离散单元法，通过将研究对象划分为大小均匀的四面体块体单元，并以单元形心为基础建立势函数定义以此计算单元间的接触力；接触检测方面提出了针对大小均匀单元接触检测的线性检测方法—NBS检测方法。

与Cundall提出的传统离散单元法相比，有限离散单元法很好地规避了不同接触形式的接触检测问题，解决了传统方法不能处理角-角接触的问题，且整个过程满足能量守恒的要求。但仍存在一些问题：应用大小均匀的四面体单元，一方面模型与实际情况不符，另一方面在实际应用时，均一化的单元尺寸以及最简单的单元形式会大大增加划分块体单元的数量，降低计算效率；单元的刚度以及法向接触力的计算受单元形式影响；接触力计算中没有考虑切向接触力；由于任意形状的凸多面体单元的形心无法将单元分割为体积相等的子单元，基于单元形心的势函数定义无法解决任意形状的凸多面体块体接触力的情况。

发明内容

本发明所要解决的技术问题是提供一种基于距离势函数三维任意凸多边形块体离散单元法，解决三维任意凸多面体块体离散单元接触力计算的问题，改善有限离散单元法的不足。本发明的方法解决了现有技术中单元尺寸均一、形式单一、单元刚度及接触力的计算受单元形式影响、没有考虑切向接触力以及无法应用于任意多面体单元等技术问题。

本发明为解决上述技术问题采用以下技术方案：

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

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

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

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

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

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

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

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

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

步骤九，根据每个时间步内每个块体单元的位移以及顶点、形心的坐标，模拟研究对象的滑坡过程。

作为本发明的进一步优化方案，步骤二中的计算时间步长须满足为：

其中，Δt是计算时间步长；ξ是系统的阻尼比，

m是块体单元质量，c是阻尼系数，k是刚度系数。

作为本发明的进一步优化方案，步骤三中，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测的具体步骤包括：

步骤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,y_K,z_K)为第K个块体单元的形心坐标，K＝1,2,…,N，N为块体单元个数，x_min、y_min、z_min分别为计算区域在x、y、z方向的最小值；

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

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

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

步骤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.2，定义距离势函数，循环计算当前时间步每个与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩；

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

作为本发明的进一步优化方案，步骤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)中距离势函数计算公式代入，积分化简得到作用于三角形面的法向接触力：

其中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_ζ、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所受的力矩。

作为本发明的进一步优化方案，在所述步骤4.3中，计算有与目标块体相接触的块体单元作用于目标块体的接触力和力矩，具体过程如下：

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

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

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

作为本发明的进一步优化方案，所述步骤七中，更新每个块体单元的顶点、形心的坐标，公式为：

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方向的分量。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：本方法采用非均匀块体离散单元接触检测法和距离势函数的定义，实现了不同大小、形态单元的接触检测及接触力计算问题，减少了实际划分单元的数量，提高了计算效率；单元刚度、法向接触力的计算不随单元形式的变化而产生差异，且考虑了切向接触力的影响，计算更符合实际，提高了离散单元数值模拟的准确性与可靠性；可以实现三维非连续介质大规模任意凸多面体块体离散单元接触力计算问题，计算过程满足能量守恒。

附图说明

图1是两接触块体示意图。

其中，1-块体单元1，2-块体单元2。

图2是块体单元2侵入块体单元1所形成的重叠区域示意图。

图3是目标单元1分割成子块体示意图。

图4是块体单元2重叠区域的底面被目标单元1的子块体分成四部分的示意图。

图5是构成重叠区域单元2的底面ABCD上是势函数分布示意图。

图6-图9是实施例岩质边坡不同计算时刻的滑坡破坏的过程示意图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

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

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

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

要使求解稳定，计算允许的最大时间步长须满足：

其中，Δt是计算时间步长；ξ是系统的阻尼比，

m是块体单元质量，c是阻尼系数，k是刚度系数。

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

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

步骤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,y_K,z_K)为第K个块体单元的形心坐标，K＝1,2,…,N，N为块体单元个数，x_min、y_min、z_min分别为计算区域在x、y、z方向的最小值；

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

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

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

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

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

2)重复步骤3.3，将第2组到第n组的块体映射到网格上，并对第2组块体与剩余所有块体进行接触检测；

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

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

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

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

下面，以可能相互接触的块体单元1和块体单元2具体说明计算过程，如图1所示，对此两个彼此接触的块体单元。

以块体单元1为目标单元，求解由块体单元2所引起、作用于目标单元1的接触力和力矩过程为：

1)确定块体单元2和块体单元1在彼此内部的顶点以及计算块体单元2各底面与块体单元1各底面的交点，确定块体单元2与块体单元1的重叠区域，如图2所示；如果没有重叠区域则判断俩单元不接触。

2)定义距离势函数：

其中，

表示目标单元1内部点p的势函数值，h_p-1,h_p-2,…,h_p-M分别表示点p至目标单元1各底面的距离，M为目标单元1的底面数，h为目标单元1的最大内切球半径。

计算出重叠区域各端点在目标单元1中的距离势函数值，若端点在目标单元1的底面上则势函数值为0，若端点到单元底边上距离等于最大内切球半径，则势函数值为1。

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

其中，f_n是法向接触力；

分别为β_c、β_t重叠区域内的点分别在β_t与β_c的距离势函数梯度；V为β_t与β_c的重叠区域。

应用高斯公式，分别在重叠区域各面上建立局部坐标系，将重叠区域内法向接触力积分计算转化为重叠区域边界面接触力积分计算：

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

4)计算当前时间步块体单元2嵌入目标单元1引起的法向接触力及力矩，具体计算过程为：

4-1)根据以目标单元1各顶点的角平分线为界、其内切圆半径为控制距离，将目标单元1分割为若干子多面体。以图3所示的八面体为例，通过八面体的各顶点角平分线及内切圆半径，将八面体分成8个子块体。

4-2)分别以块体单元2的各底面切割目标单元1的各子多面体，确定目标单元1与块体单元2重叠区域的边界面，由此，β_t与β_c重叠区域法向接触力的计算转化为β_c各底面与β_t的各子多面体之间的相互作用。

如图4所示，底面ABCD被目标单元1的八个子多面体切割为4个区域。由距离势函数计算公式可以看出，距离势函数在重叠区域内部线性分布，在重叠区域的边界面上表现为空间平面。如图5所示的底面ABCD四个子区域OAB、OBC、OCD、OAD的势函数分布示意图。

因此，在块体单元2的底面上建立局部坐标系，定义该底面上的距离势函数计算公式为：

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

4-3)对块体单元2中任一底面在目标单元1的各子多面体内的重叠面的任意一顶点作为固定点、其余各顶点与固定点的连线为分界线，将重叠面分割为若干三角形。

根据距离势函数计算公式，分别计算各三角形顶点在目标单元1内的距离势函数值；并根据各三角形顶点在目标单元1内的距离势函数值，联立方程组即可求得距离势函数方程系数a、b、c。

4-4)将重叠区域边界面局部坐标系下法向接触力计算公式的二重积分进一步进行坐标变换。将积分变换到4-3)所分割的三角形局部坐标系内，坐标变换的形函数分别为：N(1)＝1-ζ-η，N(2)＝ζ，N(3)＝η，ζ,η∈(0,1)；(ζ,η)为三角形上的点在三角形局部坐标系内的坐标。

三角形面上任意一点的坐标(x_p,y_p)可以表示为：

x_p＝N(1)x₁+N(2)x₂+N(3)x₃

y_p＝N(1)y₁+N(2)y₂+N(3)y₃

其中，(x_i,y_i)是三角形面上三个顶点的坐标，i＝1,2,3。

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

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

将4-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₃)分别为三角形三个顶点在块体单元2中任一底面局部坐标系中的坐标值；

4-5)求解由作用于三角形面的法向接触力引起的接触力矩，计算公式为：

将4-2)中距离势函数计算公式代入简化可得法向接触力引起的力矩为：

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

4-6)重复4-3)到4-5)，并对计算得到的法向接触力和力矩分别进行矢量求和，得到块体单元2的一个底面与目标单元1间的法向接触力及力矩；

4-7)重复4-2)到4-6)，分别计算块体单元2所有底面与目标单元1的法向接触力和力矩，并分别对计算得到的法向接触力和力矩进行矢量求和，得到由块体单元2嵌入目标单元1所引起的法向接触力及合力矩；

5)以块体单元2为目标单元，以块体单元1为接触单元，重复1)至4)，分别求出由块体单元1侵入块体单元2，所引起的作用于块体单元2法向接触力及力矩；并与4)所求的法向接触力和力矩矢量求和，得到当前时间步块体单元1与块体单元2间的法向接触力F_n以及其引起的合力矩M_n。

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

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

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

6-3)为使切向接触力始终保持与法向接触力垂直，保持大小一致的情况下，将上一时间步的切向接触力方向转换到当前时间步：

f_s′＝r×f_s”

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

(n₁,n₂)＝n_τ1，(n₁′,n₂′)＝n_τ2，n_τ1为当前时间步块体单元1所受法向接触力的单位方向矢量，n_τ2为上一时间步块体单元1所受法向接触力的单位方向矢量；

6-4)分别计算当前时间步，块体单元1的切向接触力，计算公式为：

F_s＝f_s′+Δf_s

同时，切向接触力需满足库伦摩擦定律：

其中，

为最大静摩擦角；c_凝为凝聚力；F_n为法向接触力；(F_s)_max为允许的最大切向接触力，即最大静摩擦力；若切向接触力F_s＞(F_s)_max时，令F_s＝(F_s)_max；

6-5)计算由切向接触力引起的作用于块体单元1的力矩：

其中，F_s为单元1所受切向接触力；

为由切向接触力作用点P_pτ指向单元1形心的矢量；

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

M＝M_n+M_s

8)同理，根据步骤6)和7)的计算公式，计算当前时间步块体单元2所受的力矩。

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

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

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

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

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

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

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

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

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,yt,z_t)是块体单元的顶点、形心的坐标，(r(t+Δt))_x、(r(t+Δt))_y、(r(t+Δt))_z分别为块体单元的位移在x_t,y_t,z_t方向的分量。

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

步骤九，根据每个时间内每个步块体单元的位移以及顶点、形心的坐标，模拟研究对象的滑坡过程。

本发明的优点在于，提出非均匀块体离散单元接触检测方法，改进了NBS接触检测方法只适用于大小均匀的块体单元体系的不足，解决了原方法存在的划分单元数量过大等问题，大大提高了检测效率；提出距离势函数计算法向接触力的方法，解决了有限离散单元法中存在的单元的刚度、法向接触力的计算受单元形式影响等问题，同时将原方法推广至空间任意凸多面体单元，使数值模拟更符合实际；提出切向接触力计算方法，改进了有限离散单元法中没有考虑切向接触力的问题，完善了理论体系，提高了离散单元法数值模拟的可靠性与准确性。

实施例：

某岩质边坡由于其内部存在软弱夹层，受地震、降雨等外部条件影响下，可能会引发滑坡的地质灾害。采用本发明提供的方法，为岩质边坡建立离散单元模型，如图6所示，共划分离散单元2653个。

图6为岩质边坡稳定时的状态，图7到图9是边坡受外部条件影响下沿滑动面产生滑坡的运动过程。基于本发明所提供的距离势函数三维块体离散单元法模拟岩质边坡滑坡过程，能清楚地描述岩质边坡受不利荷载影响下，沿滑动面破坏的过程，可以很好的分析岩质边坡在荷载作用下是否安全，若产生滑坡破坏过程、以及滑坡体形成的堆积体的形态、体积、规模等都能很直观的展示。

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

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (3)

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

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

步骤二，确定计算时间步长；计算时间步长须满足为：

其中，Δt是计算时间步长；ξ是系统的阻尼比，

m是块体单元质量，c是阻尼系数，k是刚度系数；

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

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

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

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

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

4.2.2)定义距离势函数：

其中，

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

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

其中，f_n是法向接触力；

分别为β_c、β_t重叠区域内的点分别在β_t与β_c的距离势函数梯度；V为β_t与β_c的重叠区域；

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

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

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

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

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

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

4.2.4-3)将β_c中一底面

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

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

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

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

将4.2.4-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.2.4-4)由作用于三角形面的法向接触力引起的接触力矩为：

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

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

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

与β_t间的法向接触力及力矩；

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

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

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

4.2.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为切向刚度系数；

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

f_s′＝r×f_s”

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

4.2.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；

4.2.6-5)计算由切向接触力引起的作用于β_t的力矩：

其中，F_s为β_t所受切向接触力；

为由切向接触力作用点P_pτ指向β_t形心的矢量；

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

M＝M_n+M_s

4.2.8)同理，根据步骤4.2.6)和4.2.7)的计算公式，计算当前时间步β_c所受的力矩；

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

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

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

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

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

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

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

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

步骤九，根据每个时间步内每个块体单元的位移以及顶点、形心的坐标，模拟研究对象的滑坡过程。

2.根据权利要求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,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组块体的接触检测。

3.根据权利要求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方向的分量。

Images