CN109284537B - 一种可变形二维任意圆化凸多边形离散单元法 - Google Patents

Abstract

本发明公开了一种可变形二维任意圆化凸多边形离散单元法，将块体离散单元进行圆化处理，并且考虑将离散单元法与几何非线性有限元结合，根据圆化凸多边形离散单元法求解接触力的计算方法，得到作用于块体单元的接触力；利用形函数将接触力转化成等效的单元节点的外荷载；建立可变形有限元方程组，通过求解得出离散元内部的应力和变形。本发明解决了现有的圆化凸多边形离散单元法不可变形的问题，完善了理论体系；将圆角化的多边形离散元与有限元结合，使数值模拟更加符合实际，提高了离散单元法数值模拟的可靠性与准确性；可以精确地捕捉离散体系的运动过程，精确反应离散单元内部的真实应力和变形状态。

Description

一种可变形二维任意圆化凸多边形离散单元法

技术领域

本发明属于可变形离散元技术领域，特别涉及了一种可变形二维任意圆化凸多边形离散单元法。

背景技术

离散元可以分为两大类：颗粒离散元和块体离散元。块体离散元又可以分为两大类：基于嵌入深度的模型和基于嵌入体积的模型。传统的离散元均是基于嵌入深度的，即在接触点处布置弹簧，通过定义弹簧刚度得到接触力。优点是计算速度快，缺点是接触力为集中力，且需要区分不同的接触类型，如点点，点边，点面，边边，边面，面面等。对于尖角处接触力方向不易确定，模型的鲁棒性不好，往往需要平滑尖角。基于嵌入深度的离散元，接触力均是集中力，接触力只与嵌入深度有关而与接触面积无关，这一点并不符合客观事实。颗粒离散元计算速度快，但代表性差；块体离散元可模拟任意形状，但计算效率不高。圆化凸多边形离散单元取了两者的优点，基本思想是圆化块体的边、角，可以用类似颗粒离散元的接触力计算方法来进行圆角化的块体接触力计算。

目前，英国A.MUNJIZA教授提出了基于势函数法的可变形离散元，结合离散单元法与有限单元法解决了可变形离散元问题。Munjiza利用显式解法求解有限元，避免了求解有限元非线性方程组的迭代过程。Munjiza实现了传统离散元的可变形，但是仍然存在一些问题，只能应用大小均匀的三角形或者四边形单元，一方面模型与实际情况不符，另一方面在实际应用时，均一化的单元尺寸以及最简单的单元形式会大大增加划分块体单元的数量，降低计算效率。并且工程实际中的离散元块体并不是一直保持着有棱角的状态，会随着磨损的增加，块体棱角出现磨碎，圆角化的离散元更加符合工程实际。

发明内容

为了解决上述背景技术提出的技术问题，本发明旨在提供一种可变形二维任意圆化凸多边形离散单元法，解决现有技术中圆化凸多边形离散单元不可变形的问题，使数值模拟更加符合实际。

为了实现上述技术目的，本发明的技术方案为：

一种可变形二维任意圆化凸多边形离散单元法，包括以下步骤：

(1)首先建立离散的多边形块体体系，根据研究对象的大小确定计算区域，再将离散的多边形块体单元划分成有限元网格，划分的网格单元用于可变形有限元的计算；

(2)确定时间步长；

(3)在当前时间步，确定接触单元，采用NBS接触检测方法对所有块体离散单元外围一层单元进行接触检测，得到每个块体外围一层的接触单元以及与之可能接触的单元，并且隶属于同一个离散单元的接触单元不进行接触检测；

(4)根据步(3)的接触检测结果，对可能相互接触的接触单元进行接触力计算，得到当前时间步每个接触单元所受的接触力；

(5)将步骤(4)计算得出的接触力以及作用在圆化凸多边形离散单元系统上的外力用形函数转化成载荷的等效节点力矢量；

(6)由步骤(5)计算得出的载荷的等效节点力矢量以及系统应力场的等效节点矢量，求解动力控制方程，得到当前时间步系统的各变量值，其中包括块体单元的位移；

(7)根据步骤(6)中得到的块体单元位移更新块体单元的几何信息，几何信息包含每个块体单元顶点和形心的坐标完成当前时间步的计算；

(8)重复步骤(3)-(7)计算下一时间步，直至计算完所有时间步。

进一步地，在步骤(2)中，时间步长Δt需要满足以下条件：

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

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

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

进一步地，在步骤(4)中，分别计算每个接触单元所受的法向接触力和切向接触力，则法向接触力和切向接触力的合力即为接触单元所受的接触力。

进一步地，步骤(4)的具体过程如下：

(41)首先定义圆化凸多边形离散单元：设H是多边形，B是半径为R的圆，则将H和B的闵可夫斯基和记为P，P即为圆化凸多边形离散单元，H是P的骨架；圆化半径R＝hc，h为多边形H的最大内切圆半径，c为圆化半径的系数；(42)确定圆化多边形离散单元的骨架间的最小距离：

(42-1)确定二维情况下多边形有两种接触方式：点-点接触和点-线接触；

(42-2)计算两种接触方式下骨架间的最小距离；

(43)基于嵌入深度的离散元计算公式，计算当前时间步接触力的法向以及法向接触力：

(43-1)计算接触力的法向：

其中，

是骨架H₁、H₂之间的接触力法向，

是骨架H₁、H₂之间最小距离对应的点坐标，

为骨架H₁、H₂之间的最小距离；

(43-2)计算法向接触力：

其中，

是骨架H₁、H₂之间的法向接触力，K_n是离散单元的法向刚度，δ(H₁,H₂)是两个圆化凸多边形离散单元P₁、P₂之间的重叠距离，

R₁、R₂分别为P₁、P₂的圆化半径；

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

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

F_s＝f′_s+Δf_s

其中，f′_s为上一时间步的切向接触力，Δf_s为切向接触力增量，Δf_s＝k_s·Δδ_s，k_s为切向刚度系数，Δδ_s为切向位移增量，Δδ_s＝(Δv·n_s)n_s·Δt，n_s是切向单位向量，与当前时间步接触力法向垂直，Δv是离散单元间的相对速度；

同时，当切向接触力F_s大于最大静摩擦力(F_s)_max时，令F_s＝(F_s)_max，

为最大静摩擦角，c_凝为凝聚力；

(45)计算当前时间步的接触力F＝F_n+F_s。

进一步地，步骤(6)的具体过程如下：

(61)建立参考坐标系，选择变形前的构形为参考构形；

(62)采用广义Newmark法进行时间域离散，预测上一时间步的力学量；

(63)由动力控制方程

其中，m是单元的质量，u是位移，c是阻尼系数，k是刚度系数，f是单元所受的外荷载，计算出每个块体单元的加速度增量

(64)再由广义Newmark法进行时间域离散计算出每个块体单元当前时刻的位移、速度加速度。

进一步地，在步骤(7)中，按下式更新每个块体单元的顶点或形心的坐标：

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是骨架的接触方式示意图；

图4是圆化凸多边形离散单元嵌入示意图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

本发明提出了一种可变形二维任意圆化凸多边形离散单元法，如图1所示，具体步骤如下。

步骤一，首先建立离散的多边形块体体系，根据研究对象的大小，合理确定计算区域，再将离散的多边形块体单元进一步划分成有限元网格，划分的网格单元用于可变现有限元的计算。

步骤二，确定时间步长，优选时间步长Δt满足：

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

其中，Δt_D为离散元的计算时间步长；ξ为系统的阻尼比，

m为块体单元质量，c为阻尼系数，k为刚度系数，Δt_s为有限元的时间步长，L为所有有限单元的最小边长，C的取值为10000。

步骤三，假设t时刻已经计算完毕，在当前时间步，确定接触单元，采用NBS接触检测方法对所有块体离散单元外围一层单元进行接触检测，得到每个块体外围一层的接触单元及与之可能接触的单元，并且隶属于同一个离散元的接触单元不进行接触计算和接触检测；

步骤四，根据步骤三的检测结果，对可能相互接触的接触单元进行接触力计算，并基于二维圆化凸多边形离散单元法的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力；

目标块体单元的法向接触力、切向接触力具体计算步骤如下：

1)首先确定圆化凸多边形离散单元：如图2所示，H是多边形，B是半径为R的圆，则SpheroPolyhedra就是H和B的闵可夫斯基和，记为P，P即为圆化凸多边形离散单元，H是P的骨架。确定圆化半径的系数，设多边形的最大内切圆半径为h，圆化半径的系数为c，则圆化半R径定义为R＝hc；

2)确定SpheroPolyhedra离散元的骨架间的最小距离：

2-1)首先确定二维情况下多边形有两种接触方式包括点-点接触和点-线接触；

图3所示为各骨架的接触方式，a为对齐接触，实际为点-点接触；b为错开接触，变成点-线接触；c为包含接触，变成点-线接触；d为点线接触；e为顶点接触，即为点-点接触。所以，二维情况下有两种接触方式：点-点接触和点-线接触。

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

3)基于嵌入深度的离散元计算公式，计算当前时间步接触力的法向以及法向接触力，如图4所示为圆化凸多边形离散单元嵌入示意图，其中多边形P₁的骨架为H₁，圆化半径为R₁，多边形P₂的骨架为H₂，圆化半径为R₂：

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

其中，

是骨架H₁、H₂之间的接触力法向，

是骨架H₁、H₂之间最小距离对应的点坐标，

为骨架H₁、H₂之间的最小距离；

3-2)计算法向接触力：

其中，

是骨架H₁、H₂之间的法向接触力，K_n是离散单元的法向刚度，δ(H₁,H₂)是两个圆化凸多边形离散单元P₁、P₂之间的重叠距离，

R₁、R₂分别为P₁、P₂的圆化半径；

3-3)对于P₁、P₂，先以P₁为目标单元、P₂为接触单元，按照步骤3-1)-3-2)求出由P₂嵌入P₁所引起的法向接触力，再以P₂为目标单元、P₁为接触单元，求出由P₁嵌入P₂所引起的法向接触力，两次求出的法向接触力的矢量和即为当前时间步P₁与P₂间的法向接触力F_n；

4)计算当前时间步离散单元间的切向接触力：

F_s＝f′_s+Δf_s

其中，f′_s为上一时间步的切向接触力，Δf_s为切向接触力增量，Δf_s＝k_s·Δδ_s，k_s为切向刚度系数，Δδ_s为切向位移增量，Δδ_s＝(Δv·n_s)n_s·Δt，n_s是切向单位向量，与当前时间步接触力法向垂直，Δv是离散单元间的相对速度；

同时，当切向接触力F_s大于最大静摩擦力(F_s)_max时，令F_s＝(F_s)_max，

为最大静摩擦角，c_凝为凝聚力；

5)计算当前时间步的接触力F＝F_n+F_s。

步骤五，将步骤四计算得出的接触力以及作用在圆化凸多边形离散单元系统上的外力用形函数转化成载荷的等效节点力矢量，可以采用下式计算：

其中

和

分别是当前时刻的体力和面力的载荷矢量，N是单元节点的形函数，V₀是单元的体积，A₀是单元的表面积，求得单元当前时刻荷载的等效节点力矢量

步骤六、由步骤五计算得出的载荷的等效节点力矢量以及系统应力场的等效节点矢量，求解动力控制方程，得到当前时间步系统的各变量值，具体方法如下：

1、建立参考坐标系，选择变形前的构形为参考构形；

2、采用广义Newmark法进行时间域离散，预测上一时间步的力学量；

3、由动力控制方程

其中，m是单元的质量，u是位移，c是阻尼系数，k是刚度系数，f是单元所受的外荷载，计算出每个块体单元的加速度增量

4、再由广义Newmark法进行时间域离散计算出每个块体单元当前时刻的位移、速度加速度。

步骤七，根据步骤六中得到的块体单元位移更新块体单元的几何信息，每个块体单元顶点和形心的坐标优选按照下式进行更新：

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

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

实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (4)

1.一种可变形二维任意圆化凸多边形离散单元法，其特征在于，包括以下步骤：

(1)首先建立离散的多边形块体体系，根据研究对象的大小确定计算区域，再将离散的多边形块体单元划分成有限元网格，划分的网格单元用于可变形有限元的计算；

(2)确定时间步长；

(3)在当前时间步，确定接触单元，采用NBS接触检测方法对所有块体离散单元外围一层单元进行接触检测，得到每个块体外围一层的接触单元以及与之可能接触的单元，并且隶属于同一个离散单元的接触单元不进行接触检测；

(4)根据步(3)的接触检测结果，对可能相互接触的接触单元进行接触力计算，得到当前时间步每个接触单元所受的接触力；该步骤的具体过程如下：

(41)首先定义圆化凸多边形离散单元：设H是多边形，B是半径为R的圆，则将H和B的闵可夫斯基和记为P，P即为圆化凸多边形离散单元，H是P的骨架；圆化半径R＝hc，h为多边形H的最大内切圆半径，c为圆化半径的系数；

(42)确定圆化多边形离散单元的骨架间的最小距离：

(42-1)确定二维情况下多边形有两种接触方式：点-点接触和点-线接触；

(42-2)计算两种接触方式下骨架间的最小距离；

(43)基于嵌入深度的离散元计算公式，计算当前时间步接触力的法向以及法向接触力：

(43-1)计算接触力的法向：

其中，

是骨架H₁、H₂之间的接触力法向，

是骨架H₁、H₂之间最小距离对应的点坐标，

为骨架H₁、H₂之间的最小距离；

(43-2)计算法向接触力：

其中，

是骨架H₁、H₂之间的法向接触力，K_n是离散单元的法向刚度，δ(H₁,H₂)是两个圆化凸多边形离散单元P₁、P₂之间的重叠距离，

R₁、R₂分别为P₁、P₂的圆化半径；

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

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

F_s＝f_s'+Δf_s

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

同时，当切向接触力F_s大于最大静摩擦力(F_s)_max时，令F_s＝(F_s)_max，

为最大静摩擦角，c_凝为凝聚力；

(45)计算当前时间步的接触力F＝F_n+F_s；

(5)将步骤(4)计算得出的接触力以及作用在圆化凸多边形离散单元系统上的外力用形函数转化成载荷的等效节点力矢量；

(6)由步骤(5)计算得出的载荷的等效节点力矢量以及系统应力场的等效节点矢量，求解动力控制方程，得到当前时间步系统的各变量值，其中包括块体单元的位移；

(7)根据步骤(6)中得到的块体单元位移更新块体单元的几何信息，几何信息包含每个块体单元顶点和形心的坐标完成当前时间步的计算；

(8)重复步骤(3)-(7)计算下一时间步，直至计算完所有时间步。

2.根据权利要求1所述可变形二维任意圆化凸多边形离散单元法，其特征在于，在步骤(2)中，时间步长Δt需要满足以下条件：

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

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

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

3.根据权利要求1所述可变形二维任意圆化凸多边形离散单元法，其特征在于，步骤(6)的具体过程如下：

(61)建立参考坐标系，选择变形前的构形为参考构形；

(62)采用广义Newmark法进行时间域离散，预测上一时间步的力学量；

(63)由动力控制方程

其中，m是离散单元块体质量，u是位移，c是阻尼系数，k是刚度系数，f是单元所受的外荷载，计算出每个块体单元的加速度增量

(64)再由广义Newmark法进行时间域离散计算出每个块体单元当前时刻的位移、速度加速度。

4.根据权利要求1所述可变形二维任意圆化凸多边形离散单元法，其特征在于，在步骤(7)中，按下式更新每个块体单元的顶点或形心的坐标：

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

Images