CN106529146A - 一种基于距离势函数三维任意凸多边形块体离散单元法 - Google Patents
一种基于距离势函数三维任意凸多边形块体离散单元法 Download PDFInfo
- Publication number
- CN106529146A CN106529146A CN201610952844.5A CN201610952844A CN106529146A CN 106529146 A CN106529146 A CN 106529146A CN 201610952844 A CN201610952844 A CN 201610952844A CN 106529146 A CN106529146 A CN 106529146A
- Authority
- CN
- China
- Prior art keywords
- block
- contact force
- unit
- block unit
- moment
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000011439 discrete element method Methods 0.000 title abstract 2
- 238000004364 calculation method Methods 0.000 claims abstract description 88
- 238000000034 method Methods 0.000 claims abstract description 76
- 238000001514 detection method Methods 0.000 claims abstract description 41
- 238000006073 displacement reaction Methods 0.000 claims abstract description 15
- 230000001133 acceleration Effects 0.000 claims abstract description 5
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 4
- 230000008569 process Effects 0.000 claims description 22
- 230000009466 transformation Effects 0.000 claims description 9
- 238000013016 damping Methods 0.000 claims description 6
- 238000013507 mapping Methods 0.000 claims description 6
- 230000003068 static effect Effects 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 5
- 238000005345 coagulation Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 3
- 230000003993 interaction Effects 0.000 claims description 2
- 239000011435 rock Substances 0.000 description 8
- 238000005457 optimization Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 230000006378 damage Effects 0.000 description 4
- PCTMTFRHKVHKIS-BMFZQQSSSA-N (1s,3r,4e,6e,8e,10e,12e,14e,16e,18s,19r,20r,21s,25r,27r,30r,31r,33s,35r,37s,38r)-3-[(2r,3s,4s,5s,6r)-4-amino-3,5-dihydroxy-6-methyloxan-2-yl]oxy-19,25,27,30,31,33,35,37-octahydroxy-18,20,21-trimethyl-23-oxo-22,39-dioxabicyclo[33.3.1]nonatriaconta-4,6,8,10 Chemical compound C1C=C2C[C@@H](OS(O)(=O)=O)CC[C@]2(C)[C@@H]2[C@@H]1[C@@H]1CC[C@H]([C@H](C)CCCC(C)C)[C@@]1(C)CC2.O[C@H]1[C@@H](N)[C@H](O)[C@@H](C)O[C@H]1O[C@H]1/C=C/C=C/C=C/C=C/C=C/C=C/C=C/[C@H](C)[C@@H](O)[C@@H](C)[C@H](C)OC(=O)C[C@H](O)C[C@H](O)CC[C@@H](O)[C@H](O)C[C@H](O)C[C@](O)(C[C@H](O)[C@H]2C(O)=O)O[C@H]2C1 PCTMTFRHKVHKIS-BMFZQQSSSA-N 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000004134 energy conservation Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000011229 interlayer Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000000329 molecular dynamics simulation Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
- G16Z99/00—Subject matter not provided for in other main groups of this subclass
Landscapes
- Force Measurement Appropriate To Specific Purposes (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
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组中块体单元的最大外接球半径为d1,d1的取值范围为:D≤d1<D/α;
2)第2组中块体单元的最大外接球半径为d2,d2的取值范围为:D/α≤d2<D/α2;
3)第3组中块体单元的最大外接球半径为d3,d3的取值范围为:D/α2≤d3<D/α3;
4)以此类推,第n-1组中块体单元的最大外接球半径为dn-1,dn-1的取值范围为:D/αn-2≤dn-1<D/αn-1;
5)第n组中块体单元的最大外接球半径为dn,则dn的取值范围为:dn≥D/αn-1;
其中,n=int[log(D/d)/logα+0.416]+1,D为所有块体单元中最大块体单元外接球半径,d为所有块体单元中最小块体单元外接球半径,α=2;
步骤3.3,对第1组的块体单元与所有块体单元进行接触检测,具体步骤包括:
1)取lmax=2d1,将计算区域划分为以lmax为边长的正方体网格;
2)根据公式:
将所有块体单元映射到网格上,其中,(xK,yK,zK)为第K个块体单元的形心坐标,K=1,2,…,N,N为块体单元个数,xmin、ymin、zmin分别为计算区域在x、y、z方向的最小值;
3)对每个块体单元与其周边所有相邻网格内块体单元进行接触检测,当循环到某个块体单元时,即以此块体单元所在网格为中心,判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触,其中,判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触的具体过程为:
3-1)分别计算两个块体单元的形心之间的距离dist、两个块体单元各顶点与其对应的形心之间的最大距离,并计算两个块体单元各顶点与其对应的形心之间的最大距离的平均值lavmax;
3-2)比较dist与lavmax,若dist≤2lavmax,则判断此两个块体单元接触,否则判断此两个块体单元不接触;
步骤3.4,将第一组的块体单元从所有块体单元之中去除,对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测,具体步骤包括:
1)取lmax=2d2,将计算区域划分为以lmax为边长的正方体网格;
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的势函数值,hp-1,hp-2,…,hp-M分别表示点p至βt各底面的距离,M为βt的底面数,h为βt的最大内切球半径;
3)根据势函数法计算法向接触力,计算法向接触力,计算公式为:
其中,fn是法向接触力;分别为β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的底面上建立局部坐标系,定义该底面上的距离势函数计算公式为:
其中,(xc,yc)为该底面上点的局部坐标;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)中法向接触力积分公式,计算作用于三角形面的法向接触力:
其中,fn,tri为三角形面上的法向接触力;|J|为积分坐标变换的雅可比行列式;kn为法向刚度;
将4-2)中距离势函数计算公式代入,积分化简得到作用于三角形面的法向接触力:
其中A=[a(x2-x1)+b(y2-y1)]|J|,B=[a(x3-x1)+b(y3-y1)]|J|,C=[ax1+by1+c]|J|;(x1,y1)、(x2,y2)、(x3,y3)分别为三角形三个顶点在底面局部坐标系中的坐标值;
4-4)由作用于三角形面的法向接触力引起的接触力矩为:
将4-2)中距离势函数计算公式代入简化得到作用于三角形面的法向接触力引起的力矩为:
其中,Mζ、Mη分别为作用于三角形面的法向接触力对三角形局部坐标系ζ轴、η轴的力矩;
4-5)重复4-3)到4-4),并对计算得到的法向接触力和力矩分别进行矢量求和,得到βc的一个底面aβc1与βt间的法向接触力及力矩;
4-6)重复4-2)到4-5),分别计算βc所有底面与βt的法向接触力和力矩,并分别对计算得到的法向接触力和力矩进行矢量求和,得到由βc嵌入βt所引起的法向接触力及合力矩;
5)以块体单元βc为目标单元、βt为接触单元,求出由βt嵌入βc所引起的法向接触力及力矩,并与4)所求由βc嵌入βt所引起的法向接触力和力矩分别进行矢量求和,得到当前时间步βc与βt间的法向接触力Fn以及其引起的合力矩Mn;
6)计算当前时间步,作用于βc、βt的切向接触力及其引起的力矩,具体过程为:
6-1)计算βt相对于βc的切向位移增量:Δδs=(Δv·ns)ns·Δt;其中,ns是切向单位向量,与当前时间步法向接触力方向垂直;Δν是βt相对于βc的速度;
6-2)计算βt的切向接触力增量:Δfs=ks·Δδs,其中,ks为切向刚度系数;
6-3)将上一时间步的切向接触力转换到当前时间步切向接触力增量方向:
fs'=r×fs”
其中,fs'为方向转换后上一时间步切向接触力;fs”为方向转换前上一时间步切向接触力;r为旋转矩阵;
6-4)计算当前时间步,βt的切向接触力,计算公式为:
Fs=fs'+Δfs
同时,切向接触力需满足库伦摩擦定律:其中,为最大静摩擦角;c凝为凝聚力;Fn为法向接触力;(Fs)max为允许的最大切向接触力,即最大静摩擦力;若切向接触力Fs>(Fs)max时,令Fs=(Fs)max;
6-5)计算由切向接触力引起的作用于βt的力矩:其中,Fs为βt所受切向接触力;为由切向接触力作用点Ppτ指向βt形心的矢量;
7)根据5)和6)的计算结果,计算当前时间步βt受的力矩,计算公式为:
M=Mn+Ms
8)同理,根据6)和7)的计算公式,计算当前时间步βc所受的力矩。
作为本发明的进一步优化方案,在所述步骤4.3中,计算有与目标块体相接触的块体单元作用于目标块体的接触力和力矩,具体过程如下:
1)重复步骤4.2,循环与目标块体可能相接触的每个块体单元,若两块体没有重叠区域则判断两单元不接触,若有重叠区域则求出作用于该块体单元的切向接触力、法向接触力和力矩;
2)当前时间步,目标块体所受的切向接触力、法相接触力和力矩分别为:
其中,m'与目标块体相接触的块体单元个数,(Fn)t为当前时间步目标块体所受的法向接触力,(Fs)t为当前时间目标块体所受的切向接触力,(M)t为当前时间步目标块体所受的力矩。
作为本发明的进一步优化方案,所述步骤七中,更新每个块体单元的顶点、形心的坐标,公式为:
xt(t+Δt)=xt(t)+(r(t+Δt))x
yt(t+Δt)=yt(t)+(r(t+Δt))y
zt(t+Δt)=zt(t)+(r(t+Δt))z
其中,(xt,yt,zt)是块体单元的顶点、形心的坐标,(r(t+Δt))x、(r(t+Δt))y、(r(t+Δt))z分别为块体单元的位移在xt,yt,zt方向的分量。
本发明采用以上技术方案与现有技术相比,具有以下技术效果:本方法采用非均匀块体离散单元接触检测法和距离势函数的定义,实现了不同大小、形态单元的接触检测及接触力计算问题,减少了实际划分单元的数量,提高了计算效率;单元刚度、法向接触力的计算不随单元形式的变化而产生差异,且考虑了切向接触力的影响,计算更符合实际,提高了离散单元数值模拟的准确性与可靠性;可以实现三维非连续介质大规模任意凸多面体块体离散单元接触力计算问题,计算过程满足能量守恒。
附图说明
图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组中块体单元的最大外接球半径为d1,d1的取值范围为:D≤d1<D/α;
2)第2组中块体单元的最大外接球半径为d2,d2的取值范围为:D/α≤d2<D/α2;
3)第3组中块体单元的最大外接球半径为d3,d3的取值范围为:D/α2≤d3<D/α3;
4)以此类推,第n-1组中块体单元的最大外接球半径为dn-1,dn-1的取值范围为:D/αn-2≤dn-1<D/αn-1;
5)第n组中块体单元的最大外接球半径为dn,则dn的取值范围为:dn≥D/αn-1;
其中,n=int[log(D/d)/logα+0.416]+1,D为所有块体单元中最大块体单元外接球半径,d为所有块体单元中最小块体单元外接球半径,α=2;
步骤3.3,将第1组的块体单元与所有块体进行接触检测,具体步骤包括:
1)取lmax=2d1,将计算区域划分为以lmax为边长的正方体网格;
2)根据公式:
将所有块体单元映射到网格上,其中,(xK,yK,zK)为第K个块体单元的形心坐标,K=1,2,…,N,N为块体单元个数,xmin、ymin、zmin分别为计算区域在x、y、z方向的最小值;
3)对每个块体单元与其周边所有相邻网格内块体单元进行接触检测,当循环到某个块体单元时,即以此块体单元所在网格为中心,进行判断此块体单元是否与该网格内部以及相邻网格内部单元接触,初步判断此块体单元与相邻的一个块体单元是否接触的具体过程为:
3-1)分别计算两个块体单元的形心之间的距离dist、两块体单元各顶点与其对应的形心之间的最大距离l1vmax、l2vmax,计算l1vmax、l2vmax的平均值为lavmax;
3-2)比较dist与lavmax的大小,如果dist≤2lavmax,则初步判断此两个块体单元可能接触,如果dist>2lavmax,则判断此两个块体单元不接触;
步骤3.4,将第一组的块体单元从总体之中去除,对第2组块体单元与剩余第2到第n组的块体单元进行接触检测,具体步骤包括:
1)取lmax=2d2,将计算区域划分为以lmax为边长的正方体网格;
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的势函数值,hp-1,hp-2,…,hp-M分别表示点p至目标单元1各底面的距离,M为目标单元1的底面数,h为目标单元1的最大内切球半径。
计算出重叠区域各端点在目标单元1中的距离势函数值,若端点在目标单元1的底面上则势函数值为0,若端点到单元底边上距离等于最大内切球半径,则势函数值为1。
3)根据势函数法计算法向接触力,计算法向接触力,计算公式为:
其中,fn是法向接触力;分别为β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的底面上建立局部坐标系,定义该底面上的距离势函数计算公式为:
其中,(xc,yc)为该底面上点的局部坐标;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);(ζ,η)为三角形上的点在三角形局部坐标系内的坐标。
三角形面上任意一点的坐标(xp,yp)可以表示为:
xp=N(1)x1+N(2)x2+N(3)x3
yp=N(1)y1+N(2)y2+N(3)y3
其中,(xi,yi)是三角形面上三个顶点的坐标,i=1,2,3。
根据3)中法向接触力积分公式,并计算作用于三角形面的法向接触力:
其中,fn,tri为三角形面上的法向接触力;|J|为积分坐标变换的雅可比行列式;kn为法向刚度。
将4-2)中距离势函数计算公式代入,积分化简可得作用于三角形面的法向接触力计算公式为:
其中A=[a(x2-x1)+b(y2-y1)]|J|,B=[a(x3-x1)+b(y3-y1)]|J|,C=[ax1+by1+c]|J|;(x1,y1)、(x2,y2)、(x3,y3)分别为三角形三个顶点在块体单元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间的法向接触力Fn以及其引起的合力矩Mn。
6)计算当前时间步,块体单元1与块体单元2之间的切向接触力及其引起的力矩,具体过程为:
6-1)计算切向位移增量:Δδs=(Δv·ns)ns·Δt;其中,ns是切向单位向量,与当前时间步法向接触力方向垂直;Δν是块体单元1相对于块体单元2的速度。
6-2)计算单元1的切向接触力增量:Δfs=ks·Δδs,其中,ks为切向刚度系数。
6-3)为使切向接触力始终保持与法向接触力垂直,保持大小一致的情况下,将上一时间步的切向接触力方向转换到当前时间步:
fs'=r×fs”
其中,fs'为方向转换后上一时间步切向接触力;fs”为方向转换前上一时间步切向接触力;r为旋转矩阵,(n1,n2)=nτ1,(n1',n2')=nτ2,nτ1为当前时间步块体单元1所受法向接触力的单位方向矢量,nτ2为上一时间步块体单元1所受法向接触力的单位方向矢量;
6-4)分别计算当前时间步,块体单元1的切向接触力,计算公式为:
Fs=fs'+Δfs
同时,切向接触力需满足库伦摩擦定律:其中,为最大静摩擦角;c凝为凝聚力;Fn为法向接触力;(Fs)max为允许的最大切向接触力,即最大静摩擦力;若切向接触力Fs>(Fs)max时,令Fs=(Fs)max;
6-5)计算由切向接触力引起的作用于块体单元1的力矩:其中,Fs为单元1所受切向接触力;为由切向接触力作用点Ppτ指向单元1形心的矢量;
7)根据5)和6)的计算结果,计算当前时间步1受的力矩,计算公式为:
M=Mn+Ms
8)同理,根据步骤6)和7)的计算公式,计算当前时间步块体单元2所受的力矩。
步骤4.3,计算所有与目标块体相接触的块体单元作用于目标块体的接触力和力矩,具体过程如下:
1)重复步骤4.2,循环与目标块体可能相接触的每个块体单元,若两块体没有重叠区域则判断两单元不接触,若有重叠区域则求出作用于该块体单元的切向接触力、法向接触力和力矩;
2)当前时间步,目标块体所受的切向接触力、法相接触力和力矩分别为:
其中,m'与目标块体相接触的块体单元个数,(Fn)t为当前时间步目标块体所受的法向接触力,(Fs)t为当前时间目标块体所受的切向接触力,(M)t为当前时间步目标块体所受的力矩。
步骤五,重复步骤四计算当前时间步内每个块体单元所受的接触力和力矩。
步骤六,由刚体运动方程,计算出每个块体单元当前时间步的加速度,再由velocity verlet算法计算出每个块体单元下一时间步的速度和位移。
步骤七,根据步骤六中块体单元的位移,更新每个块体单元顶点、形心的坐标,完成当前时间步的计算。
更新每个块体单元的顶点、形心的坐标,公式为:
xt(t+Δt)=xt(t)+(r(t+Δt))x
yt(t+Δt)=yt(t)+(r(t+Δt))y
zt(t+Δt)=zt(t)+(r(t+Δt))z
其中,(xt,yt,zt)是块体单元的顶点、形心的坐标,(r(t+Δt))x、(r(t+Δt))y、(r(t+Δt))z分别为块体单元的位移在xt,yt,zt方向的分量。
步骤八,重复步骤三至七计算下一时间步,直至计算完所有时间步。
本发明的优点在于,提出非均匀块体离散单元接触检测方法,改进了NBS接触检测方法只适用于大小均匀的块体单元体系的不足,解决了原方法存在的划分单元数量过大等问题,大大提高了检测效率;提出距离势函数计算法向接触力的方法,解决了有限离散单元法中存在的单元的刚度、法向接触力的计算受单元形式影响等问题,同时将原方法推广至空间任意凸多面体单元,使数值模拟更符合实际;提出切向接触力计算方法,改进了有限离散单元法中没有考虑切向接触力的问题,完善了理论体系,提高了离散单元法数值模拟的可靠性与准确性。
实施例:
某岩质边坡由于其内部存在软弱夹层,受地震、降雨等外部条件影响下,可能会引发滑坡的地质灾害。采用本发明提供的方法,为岩质边坡建立离散单元模型,如图6所示,共划分离散单元2653个。
图6为岩质边坡稳定时的状态,图7到图9是边坡受外部条件影响下沿滑动面产生滑坡的运动过程。基于本发明所提供的距离势函数三维块体离散单元法模拟岩质边坡滑坡过程,能清楚地描述岩质边坡受不利荷载影响下,沿滑动面破坏的过程,可以很好的分析岩质边坡在荷载作用下是否安全,若产生滑坡破坏过程、以及滑坡体形成的堆积体的形态、体积、规模等都能很直观的展示。
综上所述,采用基于距离势函数的三维离散单元法,能够解决三维任意形状的凸多面体离散单元计算问题。
以上所述,仅为本发明中的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉该技术的人在本发明所揭露的技术范围内,可理解想到的变换或替换,都应涵盖在本发明的包含范围之内,因此,本发明的保护范围应该以权利要求书的保护范围为准。
Claims (7)
1.一种基于距离势函数三维任意凸多边形块体离散单元法,其特征在于,包括以下步骤:
步骤一,将研究对象离散为多面体块体单元体系,并根据研究对象的大小,确定计算区域;
步骤二,确定计算时间步长;
步骤三,在当前时间步,基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测,并依据块体单元形心之间的距离、块体单元各顶点与其对应的形心之间的最大距离确定每个块体单元及与之接触的块体单元;
步骤四,根据步骤三的检测结果,对相互接触的块体单元进行接触力计算,当循环到某个块体单元时,以此块体单元为目标块体单元,与此块体单元相接触的块体单元为接触块体单元,进一步进行接触判断,并基于块体离散单元距离势函数的定义,计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩;
步骤五,重复步骤四计算当前时间步内每个块体单元所受的法向接触力、切向接触力和力矩;
步骤六,根据刚体运动方程,计算出每个块体单元当前时间步的加速度,再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移;
步骤七,根据步骤六中块体单元的位移,更新每个块体单元顶点、形心的坐标,完成当前时间步的计算;
步骤八,重复步骤三至七计算下一时间步,直至计算完所有时间步。
2.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法,其特征在于,步骤二中的计算时间步长须满足为:
其中,Δt是计算时间步长;ξ是系统的阻尼比,m是块体单元质量,c是阻尼系数,k是刚度系数。
3.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法,其特征在于,步骤三中,基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测的具体步骤包括:
步骤3.1,对每个块体单元,计算其最大外接球半径;
步骤3.2,所有块体单元根据其最大外接球半径,按从大到小的顺序分成n组,具体为:
1)第1组中块体单元的最大外接球半径为d1,d1的取值范围为:D≤d1<D/α;
2)第2组中块体单元的最大外接球半径为d2,d2的取值范围为:D/α≤d2<D/α2;
3)第3组中块体单元的最大外接球半径为d3,d3的取值范围为:D/α2≤d3<D/α3;
4)以此类推,第n-1组中块体单元的最大外接球半径为dn-1,dn-1的取值范围为:D/αn-2≤dn-1<D/αn-1;
5)第n组中块体单元的最大外接球半径为dn,则dn的取值范围为:dn≥D/αn-1;
其中,n=int[log(D/d)/logα+0.416]+1,D为所有块体单元中最大块体单元外接球半径,d为所有块体单元中最小块体单元外接球半径,α=2;
步骤3.3,对第1组的块体单元与所有块体单元进行接触检测,具体步骤包括:
1)取lmax=2d1,将计算区域划分为以lmax为边长的正方体网格;
2)根据公式:
将所有块体单元映射到网格上,其中,(xK,yK,zK)为第K个块体单元的形心坐标,K=1,2,…,N,N为块体单元个数,xmin、ymin、zmin分别为计算区域在x、y、z方向的最小值;
3)对每个块体单元与其周边所有相邻网格内块体单元进行接触检测,当循环到某个块体单元时,即以此块体单元所在网格为中心,判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触,其中,判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触的具体过程为:
3-1)分别计算两个块体单元的形心之间的距离dist、两个块体单元各顶点与其对应的形心之间的最大距离,并计算两个块体单元各顶点与其对应的形心之间的最大距离的平均值lavmax;
3-2)比较dist与lavmax,若dist≤2lavmax,则判断此两个块体单元接触,否则判断此两个块体单元不接触;
步骤3.4,将第一组的块体单元从所有块体单元之中去除,对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测,具体步骤包括:
1)取lmax=2d2,将计算区域划分为以lmax为边长的正方体网格;
2)根据步骤3.3中2)至3)的方法,将第2组到第n组的块体单元映射到网格上,并对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测;
步骤3.5,以此类推,将已进行接触检测的块体单元组从所有块体单元中去除,并进行下一组块体单元与剩余所有块体单元组的接触检测,直至完成n组块体的接触检测。
4.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法,其特征在于,步骤四中,目标块体单元的法向接触力、切向接触力以及力矩具体计算步骤如下:
步骤4.1,确定目标块体单元的最大内切球直径;
步骤4.2,定义距离势函数,循环计算当前时间步每个与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩;
步骤4.3,计算所有与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩。
5.根据权利要求4所述的一种基于距离势函数三维任意凸多边形块体离散单元法,其特征在于,步骤4.2中,计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩的具体步骤为:
1)设块体单元βt与块体单元βc为两个相接触的块体单元,其中,块体单元βc为接触单元,βt为目标单元;
2)定义距离势函数:其中,表示βt内部点p的势函数值,hp-1,hp-2,…,hp-M分别表示点p至βt各底面的距离,M为βt的底面数,h为βt的最大内切球半径;
3)根据势函数法计算法向接触力,计算法向接触力,计算公式为:
其中,fn是法向接触力;分别为β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的底面上建立局部坐标系,定义该底面上的距离势函数计算公式为:
其中,(xc,yc)为该底面上点的局部坐标;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)中法向接触力积分公式,计算作用于三角形面的法向接触力:
其中,fn,tri为三角形面上的法向接触力;|J|为积分坐标变换的雅可比行列式;kn为法向刚度;
将4-2)中距离势函数计算公式代入,积分化简得到作用于三角形面的法向接触力:
其中A=[a(x2-x1)+b(y2-y1)]|J|,B=[a(x3-x1)+b(y3-y1)]|J|,C=[ax1+by1+c]|J|;(x1,y1)、(x2,y2)、(x3,y3)分别为三角形三个顶点在底面局部坐标系中的坐标值;
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间的法向接触力Fn以及其引起的合力矩Mn;
6)计算当前时间步,作用于βc、βt的切向接触力及其引起的力矩,具体过程为:
6-1)计算βt相对于βc的切向位移增量:Δδs=(Δv·ns)ns·Δt;其中,ns是切向单位向量,与当前时间步法向接触力方向垂直;Δν是βt相对于βc的速度;
6-2)计算βt的切向接触力增量:Δfs=ks·Δδs,其中,ks为切向刚度系数;
6-3)将上一时间步的切向接触力转换到当前时间步切向接触力增量方向:
fs'=r×fs”
其中,fs'为方向转换后上一时间步切向接触力;fs”为方向转换前上一时间步切向接触力;r为旋转矩阵;
6-4)分别计算当前时间步,βt的切向接触力,计算公式为:
Fs=fs'+Δfs
同时,切向接触力需满足库伦摩擦定律:其中,为最大静摩擦角;c凝为凝聚力;Fn为法向接触力;(Fs)max为允许的最大切向接触力,即最大静摩擦力;若切向接触力Fs>(Fs)max时,令Fs=(Fs)max;
6-5)计算由切向接触力引起的作用于βt的力矩:其中,Fs为βt所受切向接触力;为由切向接触力作用点Ppτ指向βt形心的矢量;
7)根据5)和6)的计算结果,计算当前时间步βt受的力矩,计算公式为:
M=Mn+Ms
8)同理,根据步骤6)和7)的计算公式,计算当前时间步βc所受的力矩。
6.根据权利要求5所述的一种基于距离势函数三维任意凸多边形块体离散单元法,其特征在于,在所述步骤4.3中,计算有与目标块体相接触的块体单元作用于目标块体的接触力和力矩,具体过程如下:
1)重复步骤4.2,循环与目标块体可能相接触的每个块体单元,若两块体没有重叠区域则判断两单元不接触,若有重叠区域则求出作用于该块体单元的切向接触力、法向接触力和力矩;
2)当前时间步,目标块体所受的切向接触力、法相接触力和力矩分别为:
其中,m'与目标块体相接触的块体单元个数,(Fn)t为当前时间步目标块体所受的法向接触力,(Fs)t为当前时间目标块体所受的切向接触力,(M)t为当前时间步目标块体所受的力矩。
7.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法。其特征在于,所述步骤七中,更新每个块体单元的顶点、形心的坐标,公式为:
xt(t+Δt)=xt(t)+(r(t+Δt))x
yt(t+Δt)=yt(t)+(r(t+Δt))y
zt(t+Δt)=zt(t)+(r(t+Δt))z
其中,(xt,yt,zt)是块体单元的顶点、形心的坐标,(r(t+Δt))x、(r(t+Δt))y、(r(t+Δt))z分别为块体单元的位移在xt,yt,zt方向的分量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610952844.5A CN106529146B (zh) | 2016-11-03 | 2016-11-03 | 一种基于距离势函数三维任意凸多边形块体离散单元法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610952844.5A CN106529146B (zh) | 2016-11-03 | 2016-11-03 | 一种基于距离势函数三维任意凸多边形块体离散单元法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106529146A true CN106529146A (zh) | 2017-03-22 |
CN106529146B CN106529146B (zh) | 2020-10-09 |
Family
ID=58326750
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610952844.5A Active CN106529146B (zh) | 2016-11-03 | 2016-11-03 | 一种基于距离势函数三维任意凸多边形块体离散单元法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106529146B (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108446498A (zh) * | 2018-03-21 | 2018-08-24 | 河海大学 | 一种替换式逐步细化边坡块体离散元数值模型构建方法 |
CN109284537A (zh) * | 2018-08-24 | 2019-01-29 | 河海大学 | 一种可变形二维任意圆化凸多边形离散单元法 |
CN109408973A (zh) * | 2018-10-30 | 2019-03-01 | 河海大学 | 一种基于距离势函数二维可变形凸多边形块体离散单元法 |
CN109408977A (zh) * | 2018-10-30 | 2019-03-01 | 河海大学 | 一种基于距离势函数可变形三维凸多面体块体离散单元法 |
CN109492285A (zh) * | 2018-10-30 | 2019-03-19 | 河海大学 | 一种可变形三维任意圆化凸多面体块体离散单元法 |
CN110211235A (zh) * | 2019-05-14 | 2019-09-06 | 河海大学 | 基于并行rcb三维势函数离散元的放矿仿真方法 |
CN116822221A (zh) * | 2023-06-30 | 2023-09-29 | 中国科学院、水利部成都山地灾害与环境研究所 | 一种基于相互侵入势的离散岩块间接触力计算方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2004215822A (ja) * | 2003-01-14 | 2004-08-05 | Yaskawa Electric Corp | 訓練装置における抵抗力発生方法 |
CN103377307A (zh) * | 2012-04-16 | 2013-10-30 | 利弗莫尔软件技术公司 | 用于形成包含堆积在任意形状的体积中的多分散球形颗粒的计算机化模型的方法和系统 |
CN105701349A (zh) * | 2016-01-13 | 2016-06-22 | 河海大学 | 非均匀颗粒离散单元快速线性接触检测方法 |
CN105808858A (zh) * | 2016-03-10 | 2016-07-27 | 大连理工大学 | 一种重力堆积结构和随机泡沫结构的构建方法 |
CN105912852A (zh) * | 2016-04-08 | 2016-08-31 | 河海大学 | 一种基于距离势函数任意凸多边形块体离散单元法 |
JP2016181257A (ja) * | 2015-03-24 | 2016-10-13 | 東レ株式会社 | ファンデルワールス相互作用から生じる圧力の計算方法および、高分子−低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法 |
-
2016
- 2016-11-03 CN CN201610952844.5A patent/CN106529146B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2004215822A (ja) * | 2003-01-14 | 2004-08-05 | Yaskawa Electric Corp | 訓練装置における抵抗力発生方法 |
CN103377307A (zh) * | 2012-04-16 | 2013-10-30 | 利弗莫尔软件技术公司 | 用于形成包含堆积在任意形状的体积中的多分散球形颗粒的计算机化模型的方法和系统 |
JP2016181257A (ja) * | 2015-03-24 | 2016-10-13 | 東レ株式会社 | ファンデルワールス相互作用から生じる圧力の計算方法および、高分子−低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法 |
CN105701349A (zh) * | 2016-01-13 | 2016-06-22 | 河海大学 | 非均匀颗粒离散单元快速线性接触检测方法 |
CN105808858A (zh) * | 2016-03-10 | 2016-07-27 | 大连理工大学 | 一种重力堆积结构和随机泡沫结构的构建方法 |
CN105912852A (zh) * | 2016-04-08 | 2016-08-31 | 河海大学 | 一种基于距离势函数任意凸多边形块体离散单元法 |
Non-Patent Citations (1)
Title |
---|
张冲 等: "基于可见性和后修圆的三维块体域接触算法", 《岩石力学与工程学报》 * |
Cited By (13)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108446498B (zh) * | 2018-03-21 | 2021-09-14 | 河海大学 | 一种替换式逐步细化边坡块体离散元数值模型构建方法 |
CN108446498A (zh) * | 2018-03-21 | 2018-08-24 | 河海大学 | 一种替换式逐步细化边坡块体离散元数值模型构建方法 |
CN109284537B (zh) * | 2018-08-24 | 2020-12-29 | 河海大学 | 一种可变形二维任意圆化凸多边形离散单元法 |
CN109284537A (zh) * | 2018-08-24 | 2019-01-29 | 河海大学 | 一种可变形二维任意圆化凸多边形离散单元法 |
CN109408977A (zh) * | 2018-10-30 | 2019-03-01 | 河海大学 | 一种基于距离势函数可变形三维凸多面体块体离散单元法 |
CN109492285A (zh) * | 2018-10-30 | 2019-03-19 | 河海大学 | 一种可变形三维任意圆化凸多面体块体离散单元法 |
CN109408973A (zh) * | 2018-10-30 | 2019-03-01 | 河海大学 | 一种基于距离势函数二维可变形凸多边形块体离散单元法 |
CN109408977B (zh) * | 2018-10-30 | 2022-07-22 | 河海大学 | 一种基于距离势函数可变形三维凸多面体块体离散单元法 |
CN109492285B (zh) * | 2018-10-30 | 2022-08-26 | 河海大学 | 一种可变形三维任意圆化凸多面体块体离散单元法 |
CN110211235A (zh) * | 2019-05-14 | 2019-09-06 | 河海大学 | 基于并行rcb三维势函数离散元的放矿仿真方法 |
CN110211235B (zh) * | 2019-05-14 | 2022-08-19 | 河海大学 | 基于并行rcb三维势函数离散元的放矿仿真方法 |
CN116822221A (zh) * | 2023-06-30 | 2023-09-29 | 中国科学院、水利部成都山地灾害与环境研究所 | 一种基于相互侵入势的离散岩块间接触力计算方法 |
CN116822221B (zh) * | 2023-06-30 | 2024-02-23 | 中国科学院、水利部成都山地灾害与环境研究所 | 一种基于相互侵入势的离散岩块间接触力计算方法 |
Also Published As
Publication number | Publication date |
---|---|
CN106529146B (zh) | 2020-10-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106529146B (zh) | 一种基于距离势函数三维任意凸多边形块体离散单元法 | |
Favier et al. | Shape representation of axi‐symmetrical, non‐spherical particles in discrete element simulation using multi‐element model particles | |
Lu et al. | Critical assessment of two approaches for evaluating contacts between super-quadric shaped particles in DEM simulations | |
Zhao et al. | A novel Quaternion integration approach for describing the behaviour of non-spherical particles | |
CN105912852B (zh) | 一种基于距离势函数任意凸多边形块体离散单元法 | |
Latham et al. | New modelling and analysis methods for concrete armour unit systems using FEMDEM | |
Lemos | Contact representation in rigid block models of masonry | |
Fan et al. | Investigations of a new scheme for wave propagation | |
CN107766287B (zh) | 一种应用于爆炸冲击工程中的基于物质点法的随机动力学分析方法 | |
CN105260499A (zh) | 一种三维柔性梁系统力学的运动仿真方法及系统 | |
Owen et al. | Numerical modeling of industrial applications with multi-fracturing and particulate phenomena | |
Chen | Discrete element method for 3D simulations of mechanical systems of non-spherical granular materials | |
Azimi et al. | Efficient dynamics modeling for rover simulation on soft terrain | |
CN114154384A (zh) | 一种三维立方体空间的球形颗粒随机填充算法 | |
Lisitsa et al. | GPU-based implementation of discrete element method for simulation of the geological fault geometry and position | |
Ivanov | On the impulsive dynamics of M-blocks | |
CN109408977B (zh) | 一种基于距离势函数可变形三维凸多面体块体离散单元法 | |
Giorgi et al. | NLFK4ALL: An open-source demonstration toolbox for computationally efficient nonlinear Froude-Krylov force calculations | |
CN109492285B (zh) | 一种可变形三维任意圆化凸多面体块体离散单元法 | |
CN109408973B (zh) | 一种基于距离势函数二维可变形凸多边形块体离散单元法 | |
RU2611892C1 (ru) | Способ трехмерного моделирования заданного гидрогеологического объекта, реализуемый в вычислительной системе | |
Soloviev | Transformation of frequency-magnitude relation prior to large events in the model of block structure dynamics | |
CN109284537B (zh) | 一种可变形二维任意圆化凸多边形离散单元法 | |
Sypniewska-Kamińska et al. | Double pendulum colliding with a rough obstacle | |
Tijskens et al. | Strategies for contact resolution of level surfaces |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |