一种可变形三维任意圆化凸多面体块体离散单元法的制作方法

本发明涉及一种可变形离散单元法，具体涉及一种可变形三维任意圆化凸多面体块体离散单元法，属于可变形离散单元技术领域。

背景技术：

离散元单元法是专门用来解决不连续介质问题的数值模拟方法，该方法可以精确捕捉块体系统分离、滑移破坏、倾覆旋转等非连续变形特性。而可变形的离散元可以被压缩、分离或滑动。离散元可以较真实地模拟岩体中的大变形特征，因此离散单元法在理论研究和工程应用等诸多领域得到长足发展。离散单元可以分为两大类:颗粒离散元单元和块体离散单元。颗粒离散单元计算速度快，但代表性差；块体离散单元可模拟任意形状，但计算效率不高。圆化多面体离散单元法取了两者的优点。基本思想是将块体的边、角圆化，则可以用类似颗粒离散单元的接触力计算方法来进行圆角化的块体接触力计算。

目前，英国a.munjiza教授提出了基于势函数法的可变形离散单元fdem，结合离散单元法与有限单元法解决了可变形离散单元问题。但是没有考虑单元的圆角化，并且工程实际中的离散单元块体并不是一直保持着有棱角的状态，会随着磨损的增加，块体棱角出现磨碎，圆角化的离散单元更加符合工程实际,但是目前圆角化的离散单元法并没有考虑可变形，因此不完全符合工程实际。

技术实现要素：

本发明所要解决的技术问题是提供一种可变形三维任意圆化凸多面体块体离散单元法，解决三维任意凸多面体可变形离散单元接触力计算的问题，改善离散单元法的不足。本发明的方法解决了现有技术中圆化多面体离散单元不可变形，使数值模拟更加符合实际。

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

本发明提供了一种可变形三维任意圆化凸多面体块体离散单元法，包括以下步骤：

步骤一，选取研究对象，建立可变形离散单元系统的模型，可变形离散单元系统包括：多个离散单元以及将离散单元剖分网格后形成的有限单元；

步骤二，确定模型信息的时间步长δt；

步骤三，在当前时刻t，对离散单元外围一层的网格单元进行接触检测，将其中一个离散单元外围一层产生接触的网格单元定义为目标单元，与之接触的离散单元外围一层的网格单元则定义为接触单元；

步骤四，在当前时间步，确定接触单元与目标单元，判断圆化多面体的接触方式，采用三维圆化多面体离散单元法，计算当前时间步作用于目标单元与接触单元的总法向接触力和总切向接触力；

步骤五，将步骤四计算得出的接触力的合力用形函数转化成网格单元当前时刻荷载的等效节点力矢量

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

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

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

进一步的，所述步骤一中，计算时间步长δt须满足：

δt＝min(δtd,δts)

δts≤l/c

其中，δtd为离散单元的计算时间步长；ξ为系统的阻尼比，m为块体单元质量，c为阻尼系数，k为刚度系数，δts为有限单元的时间步长，l为所有有限单元的最小边长，c的取值为10000。

进一步的，在所述步骤三中，计算当前时间步作用于目标单元与接触单元的总法向接触力和总切向接触力具体过程为：

(3-1)确定三维情况下多边形的接触方式

(3-2)计算接触力的法向，表达式如下：

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

(3-3)计算法向接触力，表达式如下：

其中，是骨架h1、h2之间的法向接触力，kn是离散单元的法向刚度，δ(h1,h2)是两个圆化凸多面体离散单元p1、p2之间的重叠距离，r1、r2分别为p1、p2的圆化半径；

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

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

fs＝fs'+δfs

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

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

(3-6)计算当前时间步的接触力f＝fn+fs。

进一步地，所述步骤六中，建立系统增量形式的动力控制方程，并求解得到下一时刻t+δt有限单元的位移的步骤包括：

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

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

步骤6.3：再由广义newmark法进行时间域离散，计算出每个有限单元下一时刻t+δt的位移；

进一步地，所述步骤七中，更新每个有限单元的顶点的坐标，公式为：

x(t+δt)＝x(t)+(r(t+δt))x

y(t+δt)＝y(t)+(r(t+δt))y

z(t+δt)＝z(t)+(r(t+δt))z

其中，(x,y,z)x、y、z是块体单元的顶点的坐标，(r(t+δt))x、(r(t+δt))y、(r(t+δt))z分别为块体单元的位移在x,y,z方向的分量。

本发明所达到的有益技术效果：本方法采用圆化多面体离散单元的定义，实现了不同大小、形态单元的接触力的计算，并且提高了计算效率；可以实现三维非连续介质大规模任意圆化凸多面体可变形离散单元接触力计算问题，更加符合工程实际。

附图说明

图1本发明具体实施例中个离散单元进一步剖分成有限单元的示意图；

图2本发明具体实施例中两个网格单元发生接触重叠的示意图；

图3本发明具体实施例中圆化多面体离散单元生成的步骤示意图；

图4是本发明具体实施例中接触方式为点-点接触方式的示意图；

图5(a)是本发明具体实施例中接触方式为点在线的端点处示意图；

图5(b)是本发明具体实施例中接触方式为点在线段上不包括端点时的示意图；图6(a)是本发明具体实施例中接触方式为点在面的顶点处示意图；

图6(b)是本发明具体实施例中接触方式为点在面上不包括面的顶点示意图；

图7(a)是本发明具体实施例中两线段对齐接触示意图；

图7(b)是本发明具体实施例中两线段交错接触示意图；

图7(c)是本发明具体实施例中线段与另一根线段的端点相交示意图；

图7(d)是本发明具体实施例中线段与另一根线段相交单不包括端点示意图；

图8(a)是本发明具体实施例中线段的两端点落在面的边上示意图；

图8(b)是本发明具体实施例中线段与面的一条边相交示意图；

图8(c)是本发明具体实施例中线段与面的两条边相交单不包括端点示意图；

图8(d)是本发明具体实施例中线段与面接触但不包括面的边示意图；

图8(e)是本发明具体实施例中线段与面的顶点接触示意图；

图9(a)是本发明具体实施例中面的四个顶点与另一个面的顶点相接触示意图；

图9(b)是本发明具体实施例中面面交错接触示意图；

图9(c)是本发明具体实施例中面完全包含另一个接触面示意图；

图10是本发明具体实施例中三维圆化多面体离散单元的点-点、点-线、点-面和线-线四种接触方式示意图；

图11是本发明具体实施例中矿石在漏斗中的初始状态示意图；

图12-图16是本发明具体实施例中漏斗口闸门开启后放矿时的几何形态示意图；其中图12是漏斗口闸门开启后放矿时(7s)的几何形态示意图；

图13是漏斗口闸门开启后放矿时(15s)的几何形态示意图；

图14是漏斗口闸门开启后放矿时(22s)的几何形态示意图；

图15是漏斗口闸门开启后放矿时(43s)的几何形态示意图；

图16是漏斗口闸门开启后放矿时(112s)的几何形态示意图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本具体实施例提供了一种可变形三维任意圆化凸多面体块体离散单元法，包括以下步骤：

步骤一，选取研究对象，建立可变形离散单元系统的放矿模型，所述可变形离散单元系统包括：多个离散单元以及将离散单元剖分网格后形成的有限单元；如图一所示，离散单元剖分网格后的每个网格定义为网格单元，这个网格单元的节点坐标与有限单元的节点坐标相一致，其中离散单元的参数包括离散单元的节点坐标、质量、阻尼比、刚度，有限单元的参数包括有限单元的节点坐标、质量矩阵、阻尼矩阵、刚度矩阵；

步骤二，确定时间步长，计算时间步长δt须满足：

δt＝min(δtd,δts)

δts≤l/c

其中，δtd为离散单元的计算时间步长；ξ为离散单元的阻尼比，m为块体单元质量，c为阻尼系数，k为刚度系数，δts为有限单元的时间步长，l为所有有限单元的最小边长，在本实施例中优选地，c的取值为10000。

步骤三，在当前时间步t时刻，确定目标单元和接触单元，对所有离散单元外围一层单元进行接触检测，得到每个块体外围一层的接触单元及目标单元，并且隶属于同一个离散单元的接触单元不进行接触计算和接触检测；如图2所示,离散单元a中标注的网格单元为接触单元，离散单元b标注的网格单元为目标单元；

步骤四，在当前时间步，确定接触单元与目标单元，判断圆化多面体的接触方式，采用三维圆化多面体离散单元法，计算当前时间步作用于目标单元与接触单元的总法向接触力和总切向接触力；

计算每个与目标块体相接触的块体单元作用于目标块体的接触力的具体步骤为：

1)首先定义圆化多面体离散单元：h是多面体，b是半径为r的球，则圆化多面体就是s和b的闵可夫斯基，记为p，p即为圆化多面体离散单元，h是p的骨架。一个圆化多面体离散单元生成的步骤如图1所示，圆化多面体离散单元由多面体单元先缩小至骨架h，然后在骨架h的基础上扩大半径为r的球体，圆化步骤如图3所示。

1-1)确定圆化半径的系数，设多面体的最大内切球半径为h，圆化半径的系数为c，则圆化半r径定义为r＝hc；

2)判断圆化多面体的接触方式，确定圆化多面体离散单元的骨架间的最小距离：

2-1)首先判断三维情况下多面体的接触方式；

图4为点-点接触；图5(a)为点在线的端点处，可归结为点点接触，图5(b)为点在线段上(不包括端点)则可归结为点线接触；图6(a)为点在面的顶点出，则可归结为点点接触，图6(b)为点在面上(不包括面的顶点)，则可归结为点面接触；图7(a)为两线段对齐接触，可归结为点点接触，图7(b)为两线段交错接触，可归结为点线接触，图7(c)为线段与另一根线段的端点相交，可归结为点线接触，图7(d)为线段与另一根线段相交(但不包括端点)，则可归结为线线接触；图8(a)为线段的两端点落在面的边上，则可归结为点线接触，图8(b)为线段与面的一条边相交，可归结为线线接触和点面接触，图8(c)为线段与面的两条边相交(不包括端点)，可归结为线线接触，图8(d)为线段与面接触(不包括面的边)，可归结为点面接触，图8(e)为线段与面的顶点接触，可归结为点线接触。图9(a)面的四个顶点与另一个面的顶点相接触，可归结为点点接触，图9(b)为面面交错接触，可归结为线线接触和点面接触，图9(c)为面包含另一个接触面，可归结为点面接触，因此，三维圆化多面体离散单元的接触方式有包括点-点，点-线，点-面，线-线；四种，如图10所示。

2-2)计算各接触情况下骨架间的最小距离，点-点最小距离可以直接求；对于点-线之间的最小距离，首先是垂足必需落在端点之内(不包括端点)，其次再求点到垂足的距离；点-面最小距离是垂足必需落在面内(不包括面的边界)，其次再求点到垂足的距离；线-线最小距离是求两条线段公垂线的长度。

3)基于传统的嵌入深度的离散单元计算公式，计算接触力的法向以及法向接触力。

3-1)接触力的法向计算公式为：

其中多面体p1的骨架为h1，圆化半径为r1。多面体p2的骨架为h2，圆化半径为r2。

3-2)接触力的计算公式为：

其中kn是法向刚度的参数；代表两个几何粒子之间的重叠距离。

3-3)多面体p1和p2之间的接触力为骨架几何元素之间的作用力之和，包括点-点，点-线，点-面，线-线，

其中多面体p1的骨架为h1，圆化半径为r1，骨架h1(并非多面体)的几何元素集合为分别代表顶点(vertex)、棱边(edge)和表面(surface)，nv1,ne1,ns1分别代表h1的几何元素点、线、面的个数；多面体p2的骨架为h2，圆化半径为r2，骨架h2(并非多面体)的几何元素的集合为nv2,ne2,ns2分别代表h2的几何元素点、线、面的个数。

3-5)重复3-3)到3-4)，并对计算得到的法向接触力进行矢量求和，可以得到由p1嵌入p2所引起的法向接触力；

3-6)以块体单元p1为目标单元、p2为接触单元，求出由p2嵌入p1所引起的法向接触力并与3)所求由p1嵌入p2所引起的法向接触力矢量求和，得到当前时间步p1与p2间的法向接触力fn；

4)计算t+δt时刻，p1与p2间的切向接触力，具体过程为：

4-1)计算t+δt时刻，p1与p2间的切向接触力：fs＝fs'+δfs；同时，切向接触力需满足库伦摩擦定律：

其中，为最大静摩擦角；c凝为凝聚力；fn为法向接触力；(fs)max为最大静摩擦力。若切向接触力fs＞(fs)max时，令fs＝(fs)max；

5)根据3)和4)的计算结果，分别计算当前时间步块体单元p1、p2受的合力：

f＝fn+fs

步骤五，将步骤四计算得出的接触力的合力用形函数转化成网格单元当前时刻荷载的等效节点力矢量

其中，是网格单元当前刻荷载载荷的等效节点力矢量，和分别是当前时刻的网格单元体力和面力的载荷矢量，n是网格单元节点的形函数，v0是网格单元的体积，a0t是当前t时刻网格单元的表面积，a0是网格单元的表面积；

步骤六，根据步骤五计算得出的载荷的等效节点力矢量，建立可变形圆化多面体离散单元法系统的增量形式的动力控制方程，并进行求解，得到下一时刻t+δt有限单元的特定变量值包括位移；

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

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

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

步骤七，更新每个网格单元的顶点的坐标，公式为：

x(t+δt)＝x(t)+(r(t+δt))x

y(t+δt)＝y(t)+(r(t+δt))y

z(t+δt)＝z(t)+(r(t+δt))z

其中，(x,y,z)x、y、z是网格单元的顶点的坐标，(r(t+δt))x、(r(t+δt))y、(r(t+δt))z分别为块体单元的位移在x,y,z方向的分量。

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

实施例：

某放矿场漏斗内充满了已崩大小不等、不规则的矿石，将漏斗闸门打开，受自重作用，矿石从漏孔中放出后的运动形态，即多边形矿石块体单元在一定形状边界内的运动形态，采用本发明提供的方法，为矿石-边界建立离散元模型，如图11所示，将离散单元矿体共划分有限单元500个，漏斗及槽框有限单元1504个，漏斗及槽框固定，从而矿体的网格单元500个，漏斗及槽框网格单元1504个。

图11为矿石在漏斗中的初始状态，图12到图16为漏斗口闸门开启后，放矿时的几何形态。离散的矿体单元之间受重力作用，彼此接触，产生接触力，将产生的接触力转化成载荷的等效节点力，使得有限单元的位移发生变化，从而网格单元的顶点坐标得到更新，从而离散的矿石单元的节点坐标得到了更新，再次进行接触检测，离散单元产生新的接触和接触力，有限单元的位移再次发生改变，以此循环直到运动结束。矿体以变加速度向槽框底运动，直至与槽框底边接触，考虑摩擦和阻尼的影响，矿体的速度逐渐减小为0，矿体逐渐堆积。当所有的矿体都堆积到槽框底部时，运动结束。基于本发明所提出的一种可变形三维任意圆化凸多面体块体离散单元法模拟放矿过程，能清楚地描述崩落矿体的运动规律。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。