CN109408977B - 一种基于距离势函数可变形三维凸多面体块体离散单元法 - Google Patents

Abstract

本发明公开了一种基于距离势函数可变形三维凸多面体块体离散单元法，步骤一，选取研究对象，建立可变形离散单元系统；步骤二，确定可变形离散单元系统的时间步长Δt；步骤三，计算作用于接触单元与目标单元之间总的接触力；步骤四，将步骤三计算得出作用在接触单元以及目标单元上的总的接触力用形函数转化成载荷的等效节点力矢量R；步骤五，由步骤四计算得出的载荷的等效节点力矢量R，计算得出下一时刻t+Δt每个网格单元的位移，速度，以及加速度；步骤六，根据步骤五中网格单元的位移，更新下一时刻t+Δt每个网格单元节点坐标。本发明使数值模拟更加符合工程实际，可以精确地捕捉离散体系的运动过程，精确反应离散单元内部的真实应力和变形状态。

Description

一种基于距离势函数可变形三维凸多面体块体离散单元法

技术领域

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

背景技术

离散元单元法是专门用来解决不连续介质问题的数值模拟方法，该方法可以精确捕捉块体系统分离、滑移破坏、倾覆旋转等非连续变形特性。而可变形的离散元可以被压缩、分离或滑动。目前由英国A.MUNJIZA教授提出有限离散单元法，通过将研究对象划分为大小均匀的四面体块体单元，并以单元形心为基础建立势函数定义以此计算单元间的接触力。

A.MUNJIZA教授提出了基于势函数法的可变形离散元，结合离散单元法与有限单元法解决了可变形离散元问题。Munjiza利用显式解法求解有限元，避免了求解有限元非线性方程组的迭代过程。但仍存在一些问题：应用大小均匀的四面体单元，一方面模型与实际情况不符，单元的刚度以及法向接触力的计算受单元形式影响；基于距离势函数的离散单元法解决了这些问题，但是并没有考虑离散元的可变形，因此不太符合工程实际。

发明内容

目的：为了克服现有技术中存在的不足，本发明提供一种基于距离势函数可变形三维凸多面体块体离散单元法。

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

一种基于距离势函数三维可变形凸多面体块体离散单元法，包括以下步骤：

步骤一，选取研究对象，建立可变形离散单元系统；

步骤二，确定可变形离散单元系统的时间步长Δt；

步骤三，在当前时刻t，采用No Binary Method接触检测方法对离散单元外围一层的网格单元进行接触检测，将其中一个离散单元外围一层产生接触的网格单元定义为目标单元，与之接触的离散单元外围一层的网格单元则定义为接触单元，并且根据距离势函数的定义，计算作用于接触单元与目标单元之间总的接触力；

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

步骤五，由步骤四计算得出的载荷的等效节点力矢量

计算得出下一时刻t+Δt每个网格单元的位移，速度，以及加速度；

步骤六，根据步骤五中网格单元的位移，更新下一时刻t+Δt每个网格单元节点坐标。

作为优选方案，步骤1中所述可变形离散单元系统包括：多个离散单元，以及将离散单元剖分网格后形成的有限单元；离散单元剖分网格后的每个网格定义为网格单元，所述网格单元的节点坐标与有限单元的节点坐标相一致，所述离散单元的参数包括：离散单元的节点坐标、质量、阻尼比、刚度，所述有限单元的参数包括：有限单元的节点坐标、质量矩阵、阻尼矩阵、刚度矩阵。

作为优选方案，所述步骤2中计算时间步长Δt须满足：

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

其中，Δt_D为离散单元的时间步长；ξ为离散单元的阻尼比，

m为离散单元的质量，c为离散单元的阻尼系数，k为离散单元的刚度系数，Δt_s为有限单元的时间步长，L为所有有限单元的最小边长，C为常数；C的取值范围为9000-12000。

作为优选方案，所述C的取值为10000。

作为优选方案，所述

采用下式计算：

其中，

是网格单元当前时刻t载荷的等效节点力矢量，

和

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

作为优选方案，所述步骤5包括：根据当前时刻t可变形离散单元系统的动力控制方程

求解得到当前时刻t到下一时刻t+Δt的加速度增量

其中，M是有限单元的质量矩阵，D是有限单元的阻尼矩阵，K是有限单元的刚度矩阵，

是有限单元的加速度增量，

是有限单元的速度增量，Δu是有限单元的位移增量，再由广义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

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

其中，x(t)、y(t)、z(t)分别为网格单元在当前时刻t，x,y,z方向的节点坐标，(r(t+Δt))_x、(r(t+Δt))_y、(r(t+Δt))_z分别为网格单元的位移在下一时刻t+Δt，x,y,z方向的分量。

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

附图说明

图1为接触单元与目标单元的接触重叠示意图；

图2为接触单元与目标单元的接触重叠剖面示意图；

图3实施例岩质边坡0s滑坡破坏的过程示意图；

图4实施例岩质边坡1.6s滑坡破坏的过程示意图；

图5实施例岩质边坡2.5s滑坡破坏的过程示意图；

图6实施例岩质边坡3.5s滑坡破坏的过程示意图；

图7实施例岩质边坡5.3s滑坡破坏的过程示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

一种基于距离势函数可变形三维凸多面体块体离散单元法，包括以下步骤：

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

步骤二，确定可变形离散单元系统的时间步长Δt，计算时间步长Δt须满足：

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

其中，Δt_D为离散单元的时间步长；ξ为离散单元的阻尼比，

m为离散单元的质量，c为离散单元的阻尼系数，k为离散单元的刚度系数，Δt_s为有限单元的时间步长，L为所有有限单元的最小边长，C为常数，C的取值范围为9000-12000，这里取10000。

步骤三，在当前时刻t，采用No Binary Method(NBS)接触检测方法对离散单元外围一层的网格单元进行接触检测，将其中一个离散单元外围一层产生接触的网格单元定义为目标单元，与之接触的离散单元外围一层的网格单元则定义为接触单元，如图2所示,离散单元A中阴影部分的网格单元为接触单元，离散单元B阴影部分的网格单元为目标单元，并且根据距离势函数的定义，计算作用于接触单元与目标单元之间总的接触力；

步骤四，将步骤三计算得出作用在接触单元以及目标单元上的总的接触力用形函数转化成载荷的等效节点力矢量，采用下式计算：

其中，

是网格单元当前时刻t载荷的等效节点力矢量，

和

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

步骤五，由步骤四计算得出的载荷的等效节点力矢量

根据当前时刻t可变形离散单元系统的动力控制方程

求解得到当前时刻t到下一时刻t+Δt的加速度增量

其中，M是有限单元的质量矩阵，D是有限单元的阻尼矩阵，K是有限单元的刚度矩阵，

是有限单元的加速度增量，

是有限单元的速度增量，Δu是有限单元的位移增量，再由广义Newmark法进行时间域离散，计算出每个有限单元下一时刻t+Δt的位移，速度，以及加速度，由于网格单元与有限单元的节点坐标相一致，因此计算得出下一时刻t+Δt每个网格单元的位移，速度，以及加速度；

步骤六，根据步骤五中网格单元的位移，更新下一时刻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(t)、y(t)、z(t)分别为网格单元在当前时刻t，x,y,z方向的节点坐标，(r(t+Δt))_x、(r(t+Δt))_y、(r(t+Δt))_z分别为网格单元的位移在下一时刻t+Δt，x,y,z方向的分量。

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

实施例:

某岩质边坡由于其内部存在软弱夹层，受地震、降雨等外部条件影响下，可能会引发滑坡的地质灾害。采用本发明提供的方法，为岩质边坡建立可变形离散单元系统，模拟滑坡体在重力的作用下以一定速度产生滑坡破坏的过程。如图3所示，定义滑坡体和坡身为离散单元，将滑坡体和坡身剖分网格形成有限单元，滑坡体网格单元为197个，坡身网格单元为471个。

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

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

Claims (6)

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

步骤一，选取研究对象，建立可变形离散单元系统；

步骤二，确定可变形离散单元系统的时间步长Δt；

步骤三，在当前时刻t，采用No Binary Method接触检测方法对离散单元外围一层的网格单元进行接触检测，将其中一个离散单元外围一层产生接触的网格单元定义为目标单元，与之接触的离散单元外围一层的网格单元则定义为接触单元，并且根据距离势函数的定义，计算作用于接触单元与目标单元之间总的接触力；

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

步骤五，由步骤四计算得出的载荷的等效节点力矢量

计算得出下一时刻t+Δt每个网格单元的位移，速度，以及加速度；

步骤六，根据步骤五中网格单元的位移，更新下一时刻t+Δt每个网格单元节点坐标；

所述

采用下式计算：

其中，

是网格单元当前时刻t载荷的等效节点力矢量，

和

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

2.根据权利要求1所述的一种基于距离势函数三维可变形凸多面体块体离散单元法，其特征在于：步骤一中所述可变形离散单元系统包括：多个离散单元，以及将离散单元剖分网格后形成的有限单元；离散单元剖分网格后的每个网格定义为网格单元，所述网格单元的节点坐标与有限单元的节点坐标相一致，所述离散单元的参数包括：离散单元的节点坐标、质量、阻尼比、刚度，所述有限单元的参数包括：有限单元的节点坐标、质量矩阵、阻尼矩阵、刚度矩阵。

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

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

其中，Δt_D为离散单元的时间步长；ξ为离散单元的阻尼比，

m为离散单元的质量，c为离散单元的阻尼系数，k为离散单元的刚度系数，Δt_s为有限单元的时间步长，L为所有有限单元的最小边长，C为常数；C的取值范围为9000-12000。

4.根据权利要求3所述的一种基于距离势函数三维可变形凸多面体块体离散单元法，其特征在于：所述C的取值为10000。

5.根据权利要求1所述的一种基于距离势函数三维可变形凸多面体块体离散单元法，其特征在于：所述步骤五包括：根据当前时刻t可变形离散单元系统的动力控制方程

求解得到当前时刻t到下一时刻t+Δt的加速度增量Δü；其中，M是有限单元的质量矩阵，D是有限单元的阻尼矩阵，K是有限单元的刚度矩阵，Δü是有限单元的加速度增量，

是有限单元的速度增量，Δu是有限单元的位移增量，再由广义Newmark法进行时间域离散，计算出每个有限单元下一时刻t+Δt的位移，速度，以及加速度，由于网格单元与有限单元的节点坐标相一致，因此计算得出下一时刻t+Δt每个网格单元的位移，速度，以及加速度。

6.根据权利要求1所述的一种基于距离势函数三维可变形凸多面体块体离散单元法，其特征在于：所述更新下一时刻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(t)、y(t)、z(t)分别为网格单元在当前时刻t，x,y,z方向的节点坐标，(r(t+Δt))_x、(r(t+Δt))_y、(r(t+Δt))_z分别为网格单元的位移在下一时刻t+Δt，x,y,z方向的分量。

Images