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

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

背景技术：

离散元单元法是专门用来解决不连续介质问题的数值模拟方法，该方法可以精确捕捉块体系统分离、滑移破坏、倾覆旋转等非连续变形特性。而可变形的离散元可以被压缩、分离或滑动。离散元可以较真实地模拟岩体中的大变形特征，因此离散单元法在理论研究和工程应用等诸多领域得到长足发展。

目前由英国a.munjiza教授提出有限离散单元法，通过将研究对象划分为大小均匀的四面体块体单元，并以单元形心为基础建立势函数定义以此计算单元间的接触力。

a.munjiza教授提出了基于势函数法的可变形离散元，结合离散单元法与有限单元法解决了可变形离散元问题。munjiza利用显式解法求解有限元，避免了求解有限元非线性方程组的迭代过程。但仍存在一些问题：应用大小均匀的四面体单元，一方面模型与实际情况不符，另一方面在实际应用时，均一化的单元尺寸以及最简单的单元形式会大大增加划分块体单元的数量，降低计算效率。基于距离势函数的离散单元法解决了这些问题，但是并没有考虑离散元的可变形，因此并不完全符合工程实际。

技术实现要素：

本发明所要解决的技术问题是现有技术中基于距离势函数任意凸多边形块体离散单元法不可变形，提供了一种基于距离势函数二维可变形凸多边形块体离散单元法，改善离散单元法的不足使数值模拟更加符合实际。

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

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

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

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

步骤三：对离散单元外围一层的网格单元进行接触检测；

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

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

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

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

进一步地，在步骤二中，所述时间步长δt须满足：

δt＝min(δtd,δts)，

δts≤l/c，

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

进一步地，在步骤四中，计算当前时刻每个接触单元与目标单元的接触力时采用基于距离势函数的离散单元法的接触力算法。

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

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

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

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

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

进一步地，在步骤六中，更新每个网格单元的顶点的坐标公式为：

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

本发明所达到的技术效果：本方法实现了离散元体系的可变形，使得离散元模型更加精确得反应块体内部的应力和应变情况，可用于模拟更多的工程实际问题；实现了不同大小、形态单元的接触检测及接触力计算问题，减少了实际划分单元的数量，提高了计算效率。

附图说明

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

图2本发明方法具体实施例接触单元与目标单元的接触重叠示意图；

图3～图6实施例岩质边坡不同计算时刻的滑坡破坏的过程示意图；

其中图3为具体实施例的岩质边坡稳定时的状态；

图4～图6是具体实施例中边坡受外部条件影响下沿滑动面产生滑坡的运动过程。

具体实施方式

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

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

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

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

δt＝min(δtd,δts)

δts≤l/c

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

步骤三，离散单元外围一层的网格单元进行接触检测确定目标单元和接触单元的方法为将离散单元外围一层产生接触的网格单元定义为目标单元，与之接触的离散单元外围一层的网格单元则定义为接触单元。如图2所示根据基于距离势函数任意凸多边形块体离散单元法的定义，计算当前时间步作用于目标单元和接触单元的法向接触力、切向接触力；

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

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

步骤五，由步骤四计算得出的载荷的等效节点力矢量求解系统增量形式的动力控制方程，得到有限单元的位移；

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

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

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

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

在其它的具体实施例中，可以通过本发明提出的系统增量形式的动力控制方程求解有限单元的其它变量值如每个网格单元t+δt时刻的速度或者加速度。

步骤六，根据步骤五中有限单元的位移，更新每个有限单元顶点的坐标等几何信息，完成当前时间步的计算；

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

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

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

实施例:

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

图3为岩质边坡稳定时的状态。

图4到图6是边坡受外部条件影响下沿滑动面产生滑坡的运动过程。当滑坡体离散单元受重力开始滑动时，便与坡身发生接触，滑坡体离散单元与滑坡体离散单元之间、滑坡体离散单元与坡身离散单元之间由于接触产生接触力，再将接触力通过载荷的等效节点矢量转化到有限单元上去，使得有限单元的位移发生变化，更新了网格单元的节点坐标，从而离散单元的节点坐标得到了更新，再次进行接触检测，发生新的接触，产生新的接触力，有限单元的位移再次发生改变，以此循环直到运动结束。基于本发明所提供的距离势函数二维可变形块体离散单元法模拟岩质边坡滑坡过程，能清楚地描述岩质边坡受不利荷载影响下，沿滑动面破坏的过程，可以很好的分析岩质边坡在荷载作用下是否安全，若产生滑坡破坏过程、以及滑坡体形成的堆积体的形态、体积、规模等都能很直观的展示。

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