《离散元 课件.ppt》由会员分享，可在线阅读，更多相关《离散元 课件.ppt（103页珍藏版）》请在三一办公上搜索。
1、第6章 离散单元法,The Theory of DEM(Discrete Element Method)and Its Applications,俞 缙,华侨大学岩土工程研究所,离散单元法的基本理论及应用,历史由来及研究现状 基本原理 程序设计及商业软件介绍 应用,离散单元法的基本理论及应用,历史由来及研究现状 产生背景 发展现状存在的问题,离散单元法的基本理论及应用,基本原理Cundall二维圆盘单元离散元法三维球体单元(颗粒元)离散元法多边形单元离散元法多面体单元离散元法接触发现算法Cundall公共面法其他方法,离散单元法的基本理论及应用,程序设计及商业软件介绍 程序设计方法 UDEC
2、3DEC PFEC2D与PFC3D,离散单元法的基本理论及应用,应用 理想散体岩土介质中的应用 堆积问题磨矿铁路道砟,离散单元法的基本理论及应用,应用 在连续岩土介质中的应用岩石力学性质分析颗粒破碎模拟边坡稳定性分析洞室围岩稳定性分析桩与土介质相互作用模拟分层岩土介质中侵彻与爆炸模拟,一 历史由来及研究现状,产生背景 散粒岩土材料在自然界中普遍存在 从本质上讲，岩土材料都是由离散的、尺寸不一、形状各异的颗粒或块体组成的，例如，土就是松散颗粒的堆积物，同样，天然岩体也是由被结构面切割而成的大小不一、形态各异的岩石块体所组成。散粒岩土材料的力学特性有着重要的工程应用，如泥砂的沉淀，土堤、土（岩）坡
3、、铁路道渣等的稳定性研究，散粒岩土材料的力学特性研究是岩土力学中最基本的、也是最重要的问题之一。,一 历史由来及研究现状,产生背景 用连续介质力学研究散粒岩土材料力学特性的不足 连续介质力学把散粒体作为一个整体来考虑，研究的重点放在建立粒子集合的本构关系，从粒子集合整体的角度研究散粒体介质的力学行为。不足：不能体现颗粒间的复杂相互作用及高度非线性行为；不能真实刻画散体材料的流动变形特征。,用多体动力学描述散粒体的力学行为的困难（1）对于被研究的多粒子系统而言，已经存在的接触不断地分开，而新的接触频繁的形成，在多体动力学中，接触的分开与形成都需要改变控制方程；（2）即使接触网络保持相同，在每一个
4、接触中，也可能发生在依附与滑动间的过渡，而这种过渡也会导致系统运动方程的改变。因而：多体动力学方法只能描述少数散体体系的力学行为，对于大量散体组成的岩土材料则相当困难。,一 历史由来及研究现状,分子动力学方法的引入 分子动力学模拟是一种用来计算一个经典多体体系的平衡河传递性质的方法。所谓的经典意味着颗粒体系的运动遵守经典力学定律。该方法最初是用来描述分子运动的（当处理一些较轻的原子或分子时，才需要考虑量子效应）。分子动力学方法模拟分子的运动时，邻近分子间存在吸引或排斥力。该方法可以模拟大量分子的运动。去除分子间作用力，把分子动力学中的小尺度粒子作为散体岩土材料中的颗粒，并入颗粒间及颗粒与边界间
5、的相互作用描述，即是Cundall离散元法的最初思路。,一 历史由来及研究现状,早期的离散单元法 离散元法的思想源于较早的分子动力学（Molecular Dynamics）。1971年Cundall提出适于岩石力学的离散元法；1979年Cundall和Strack又提出适于土力学的离散元法，并推出二维圆盘（Disc）程序BALL和三维圆球程序TRUBAL（后来发展为商业软件PFC-2D/3D），形成较系统的模型与方法，被称为软颗粒模型。1988年Cundall所在的ITASCA咨询公司推出针对三维块体元的3DEC程序。至此，离散元的理论体系基本形成。,一 历史由来及研究现状,早期的离散单元法
6、Cundall称之为“Distinct Element Method”，随着该方法的推广，有的学者称其为“Discrete Element Method”，缩写形式均为DEM。最初，离散元的研究对象主要是岩石等非连续介质的力学行为，它的基本思想是把不连续体分离为刚性元素的集合，使各个刚性元素满足运动方程，用时步迭代的方法求解各刚性元素的运动方程，继而求得不连续体的整体运动形态。离散元法允许单元间的相对运动，不一定要满足位移连续和变形谐调条件，计算速度快，所需存储空间小，尤其适合求解大位移和非线性问题。,一 历史由来及研究现状,早期的离散单元法 主要思路：用Newton定律描述颗粒运动，通过颗粒
7、间及颗粒与边界间的相互作用传递载荷，求解方法是解藕的。理论难点：接触力模型（Contact Force Model)与接触发现算法(Contact Detection Algorithm)。,一 历史由来及研究现状,离散单元法 的研究现状离散元法自问世以来，在岩土工程和粉体(颗粒散体)工程这两大传统的应用领域中发挥了其它数值算法不可替代的作用。,一 历史由来及研究现状,离散单元法 的研究现状岩土工程中的应用 由于离散元单元具有更真实地表达节理岩体的几何特点能力，便于处理以所有非线性变形和破坏都集中在节理面上为特征的岩体破坏问题，被广泛地应用于模拟边坡、滑坡和节理岩体下地下水渗流等力学过程的分析
8、和计算中；离散元法还可以在颗粒体模型基础上通过随机生成算法建立具有复杂几何结构模型，通过单元间多种连接方式来体现土壤等多相介质间的不同物理关系，从而更有效地模拟土壤的开裂、分离等非连续现象，成为分析和处理岩土工程问题的不可缺少的方法。,一 历史由来及研究现状,离散单元法 的研究现状粉体工程中的应用 其次，在粉体工程(过程)方面，颗粒离散元被广泛地应用于粉体在复杂物理场作用下的复杂动力学行为的研究和多相混合材料介质或具有复杂结构的材料其力学特性的研究中它涉及到粉末加工、研磨技术、混合搅拌等工业加工和粮食等颗粒离散体的仓储和运输等生产实践领域中。,一 历史由来及研究现状,离散单元法 的研究现状离散
9、元理论研究的发展 近30年来，离散元法的应用领域在不断扩大，它自身的内涵也发生了变化，以致于目前很难对离散元法给出一个严格的定义。下面，我们从离散元法的离散模型特点及便于甄别与其它数值计算方法的关系的角度给予离散元法一个比较宽松的定义。,一 历史由来及研究现状,离散单元法 的研究现状离散元理论研究的发展数值方法通常将实际具有无限自由度的介质近似为具有有限自由度的离散体(或网络)的计算模型(有限离散模型)进行计算。有限离散模型具有三个要素：单元(或网络)、节点和节点间的关联。,一 历史由来及研究现状,离散单元法 的研究现状离散元理论研究的发展离散元单元的形状有形形色色，但它只有一个基本节点(取单
10、元的形心点)，是一种物理元(physicalelement)这种单元与有限元法、边界元法等数值方法采用的由一组基本节点联成的单元(一般称为网络元，mesh element)相比有明显的不同。另外，离散元法的节点间的关联又具有明确的物理意义，同差分法等数值方法从数学上建立节点间的关联又有明显的差异因此，我们可以将离散法简单地定义为：通过物理元的单元离散方式并构成具有明确物理意义的节点关系来建立有限离散模型的数值计算方法。,一 历史由来及研究现状,离散化模型在物体的离散化方面，离散元法的离散思想同有限元法有着相似之处：将所研究的区域划分成各种单元，并通过节点建立单元间的联系离散元法的单元从几何形状
11、上分类可分为颗粒元和块体元两大类，如图1所示块体元中最常用的有4面体元、6面体元；对于二维问题可以是任意多边形元，但应用范围不广每个离散单元只有一个基本节点(取形心点)颗粒元主要是采用球体元；对于二维问题采用圆盘形单元还有人采用椭球体单元和椭圆形，但不常用,二 基本原理,离散化模型,二 基本原理,图1 颗粒元与块体元示意图,二 基本原理-球形颗粒元离散元法,俞 缙,华侨大学岩土工程研究所,根据离散化模型中所采用的单元种类分别介绍离散元法的基本原理：颗粒元二维圆盘单元三维圆球单元块体元多边形单元多面体单元,二 基本原理,基本假设假定速度和加速度在每个时间步长内为常量；选取的时间步长应该足够小以至
12、于在单个时间步长内扰动的传播不会超过当前与之相邻的粒子。,二 基本原理-球形颗粒元离散元法,运动描述 处于一个理想散体中的任意一个颗粒,具有6个自由度,3个平动自由度与三个转动自由度,可通过Newton第二定律分别描述。,二 基本原理-球形颗粒元离散元法,运动描述平动方程：,基本原理-球形颗粒元离散元法,式中，与 分别为颗粒 的质量和速度。为时间，为颗粒的重力，与 分别为颗粒 与 的接触力与粘性接触阻尼力，为所有与颗粒接触的颗粒总数。,运动描述接触力的分解：颗粒 与 间的接触力可分解为法向与切向接触力，即,基本原理-球形颗粒元离散元法,同理，粘性接触阻尼力也可分解为法向与切向分量形式，即,运动
13、描述接触力产生的力矩：,基本原理-球形颗粒元离散元法,颗粒间的接触力作用在两个颗粒的接触点上，而不是作用在颗粒的中心，所以这些接触力（除法向接触力 外）将会对颗粒产生力矩，,式中，为从颗粒 的质心指向接触点的矢量，其幅值为（颗粒的半径）。,运动描述转动方程：,基本原理-球形颗粒元离散元法,转动方程可以表示为,式中，与 分别为颗粒 的转动惯量与角速度，对于球形颗粒 为,接触模型综述：,二 基本原理-球形颗粒元离散元法,关于接触力的计算模型已有大量的研究成果，目前仍旧是一个活跃的研究领域，特别是对于切向力的计算方法。对于理想散体颗粒（无粘连）：采用Hertz理论描述法向作用，而采用Mindlin与
14、Deresiewicz理论描述切向作用；对于存在粘连的散体颗粒：法向接触力根据在Hertz理论基础上考虑粘连力的JKR（Johnson-Kendall-Roberts）理论确定，切向接触力增量则根据把Savkoor和Briggs理论与Mindlin和Deresiewicz理论相结合形成的理论确定。,接触模型两个处于接触颗粒单位法向和切向向量：,基本原理-球形颗粒元离散元法,单位法向向量,单位切向向量,单位切向量之所以通过两个颗粒的相对速度来计算，是因为接触力与粘性阻尼力的方向与相对速度的方向相同。,接触模型两个处于接触颗粒接触点的相对速度：,基本原理-球形颗粒元离散元法,法向相对速度为,切向相
15、对速度为,或者写为,接触模型法向接触力计算模型 Hertz模型：,基本原理-球形颗粒元离散元法,为颗粒i与j接触时的侵入深度,式中,接触模型法向接触力计算模型 Cundall模型：,基本原理-球形颗粒元离散元法,式中，为法向弹簧刚度。,接触模型法向接触力计算模型 法向粘性接触阻尼力：,二 基本原理-球形颗粒元离散元法,式中，为法向粘性接触阻尼系数。,接触模型切向接触力计算模型 综述：,二 基本原理-球形颗粒元离散元法,处于接触中的两个颗粒的切向作用，从本质上讲，是一种摩擦行为，按照摩擦机理，摩擦力包括：滑动摩擦、滚动摩擦与静摩擦，其中滑动摩擦与静摩擦属于切向摩擦力；滚动摩擦是由于法向接触应力的
16、不均匀分布产生的。介绍两个切向接触力模型：Coulomb准则Mindlin与Deresiewicz切向接触力模型,接触模型切向接触力计算模型 Coulomb准则：,二 基本原理-球形颗粒元离散元法,在离散元模拟中，一般用Coulomb准则这种简单的形式描述，静摩擦的详细刻画需要涉及切向位移甚至可能要考虑时间依赖效应。,式中，为静摩擦系数，切向摩擦力的方向为与相对滑动的趋势相反。,接触模型切向接触力模型 Mindlin与Deresiewicz模型：,二 基本原理-球形颗粒元离散元法,式中 为颗粒 与 间的累积切向位移矢量,接触模型切向接触力模型 阻尼力：,二 基本原理-球形颗粒元离散元法,式中，
17、为切向粘性接触阻尼系数。,二 基本原理-球形颗粒元离散元法,计算模型总结运动方程接触力的计算法向接触力,切向接触力,二 基本原理-球形颗粒元离散元法,运动方程的求解,一般采用两种方法求解运动方程：中心差分法Verlet积分法,二 基本原理-球形颗粒元离散元法,运动方程的求解中心差分法,运动方程可由Verlet显式积分求解。通过积分可获得粒子的新位置，积分时需要粒子的当前及上一步长的位置数据，而不需要粒子的速度数据。,二 基本原理-球形颗粒元离散元法,运动方程的求解Verlet积分法,二 基本原理-球形颗粒元离散元法,运动方程的求解Verlet积分法,定义：,从而导出Verlet方程为,接触发现
18、算法 在一个由众多颗粒组成的体系中，直接判别颗粒是否接触需要耗费大量的计算时间，因而，为了节约计算时间，提高计算效率，一般不直接判别任意两个颗粒间是否存在接触，而是分两个步骤判别颗粒间的接触是否存在：首先，对一个颗粒，判别其潜在的邻居个数，然后，准确确定该颗粒与每个邻居是否接触。虽然在确定邻居数目时也要耗费一定的计算时间，但是仍旧比逐个准确判别颗粒间接触是否存在要节约时间。因而，接触发现算法的效率在多颗粒体系力学行为模拟中至关重要。,二 基本原理-球形颗粒元离散元法,接触发现算法介绍三种针对球形颗粒的接触发现算法：Verlet邻居目录法连接单元法 边界盒法,二 基本原理-球形颗粒元离散元法,接
19、触发现算法Verlet邻居目录法,二 基本原理-球形颗粒元离散元法,当需要判别体系中某个颗粒的邻居数量时，在该粒子周围构建一个球（称之为参考球，称该颗粒为核心颗粒），参考球半径为体系中最大粒子半径的若干倍，那么参考球所包围的所有粒子为该球中心粒子的邻居。参考球半径的选取取决于粒子的运动速度及体系中粒子的密度。对于每个粒子，都可生成一个邻居粒子的目录。为了得到邻居目录，对每一个粒子而言，所有标号大于该粒子的粒子都必须被检验，判断是否位于该粒子的参考球中，而对于标号小于该粒子标号的粒子则没有必要被检验，因为邻居是互相的，没有必要对一个邻居对检验两次。,Verlet邻居目录法及粒子存储目录,二 基本
20、原理-球形颗粒元离散元法,接触发现算法Verlet邻居目录法,二 基本原理-球形颗粒元离散元法,接触发现算法Verlet邻居目录法,对于n个颗粒组成的体系，用Verlet邻居目录法需要 次计算，也就是说计算次数仍旧为 量级。然而，并不需要在每个时间步长上都对邻居目录进行更新。更新的频率取决于体系中粒子的密度、粒子的运动速度以及参考球的尺寸。参考球的半径也可以根据颗粒体系的稠密程度及运动速度进行动态调整，并且参考球半径与邻居目录的更新频率呈反比关系，参考球半径越小，邻居目录的更新频率越高；但是参考球半径越大，则有更多的粒子位于球体内，所以判别是否为邻居就需要较长的时间。,接触发现算法连接单元法,
21、二 基本原理-球形颗粒元离散元法,将颗粒体系所占据的空间划分成规则的网格，对于三维问题，可以划分为 个立方体单元，对于二维体系，则划分为 个正方形单元，对于颗粒体系所占据空间形状不规则时，也可采用其他形状单元划分。但是所有单元的尺寸必须大于粒子的尺寸。与Verlet邻居目录法的主要区别在于：相邻单元法中的单元不依附于粒子，单元不随粒子的运动而运动。如果粒子就当前的位置被分配到某个单元，显然，只有在同一个单元或直接相邻单元内的粒子间才可能发生相互作用，也就是说，只有相邻单元内的粒子才能成为邻居。,接触发现算法连接单元法,二 基本原理-球形颗粒元离散元法,例如，对一个二维体系而言，只有在9个不同的
22、单元内可能包含邻居粒子，对于三维体系，则只有在27个不同单元内可能包含邻居粒子。与前面介绍的Verlet邻居目录法相同，每个粒子对只需要检验一次，这样，没必要对9个单元都进行检验，只需检验中心单元及邻居单元的一半即可，即在二维情况下，只需检验5个单元，在三维情况下，只需检验14个单元。,接触发现算法连接单元法,二 基本原理-球形颗粒元离散元法,2D体系中的连接单元法,接触发现算法边界盒法,二 基本原理-球形颗粒元离散元法,这个方法与前述两种方法不同。首先，在每一个粒子周围构建一个边界盒，边界盒的尺寸按这样的方式选取：使每个粒子刚好放进它的边界盒内。边界盒的边为直线，并且与体系的坐标轴平行。在判
23、别颗粒的邻居时，把边界盒投影到体系的坐标轴上。通过边界盒在坐标轴上投影的起点和终点来判别是否为邻居。,接触发现算法边界盒法,二 基本原理-球形颗粒元离散元法,a 围绕在每个粒子周围的边界盒,b两个不同时刻粒子的边界盒在轴上的投影,在判别颗粒的邻居时，把边界盒投影到体系的坐标轴上。例如，图a表明了粒子及边界盒的位置，图b为图a中的边界盒在体系x轴的投影。通过边界盒在坐标轴上投影的起点和终点来判别是否为邻居，出于这个原因，投影的起点和终点序列被存储在目录中。,接触发现算法边界盒法,二 基本原理-球形颗粒元离散元法,对于3D体系，必须把边界盒子在三个坐标轴上投影，所以生成三个目录。如果一个粒子的边界
24、盒在某个坐标轴上投影的起点和终点间包含另外一个粒子边界盒投影的起点、终点或起点和终点，就说明这两个粒子的边界盒在该坐标轴上的投影发生了重叠。如果两个边界盒的投影在每一个坐标轴上都发生重叠，那么就说明这两个边界盒发生了接触。,接触发现算法边界盒法,二 基本原理-球形颗粒元离散元法,存储目录的更新检查一个边界盒在每一个坐标轴的投影是否在另外一个边界盒投影的起点和终点之间仍旧需要花费很多计算时间。但是，尽管这些目录在每个计算时间步长上都必须被更新，所用的计算时间仍旧可以被减少到与体系中粒子数量成比例的量级，因为在每一个新步长上所需要做的只是对旧目录的更新，并且对旧目录的更新只是对旧目录的重新分类，因
25、而，更新的过程非常简单，只需要对旧目录进行顺序地检查，判断是否在顺序上有新的变化即可。这种变化仅仅是位置的变化。,基本原理-球形颗粒元离散元法,算法设计,基本原理-球形颗粒元离散元法,计算所需参数,只有当选取的时间步长小于临界时间步长时，数值计算才是稳定的，这是模型采用显式积分的结果。临界时间步长可以通过计算估计，然后，根据临界时间步长确定计算所采用的时间步长。可以采用质量为m、用刚度k为的弹簧与地面相连接的单自由度体系来估计，这种情况下，临界步长为:,二 基本原理-球形颗粒元离散元法,时间步长的确定,三 基本原理-块体元离散元法,根据离散化模型中所采用的单元种类分别介绍离散元法的基本原理：颗
26、粒元二维圆盘单元三维圆球单元块体元多边形单元多面体单元,三 基本原理,基本思路将所研究的区域或散体体系剖分或看作多边形块体单元的集合；单元间的接触方式包括：角-角接触、角-边接触、边-边接触三种形式；随着单元的平移和转动，各个单元间的接触关系由运动方程动态调整；最终，块体单元达到平衡状态或一直运动下去；单元性质可以假设为刚性（对于应力水平低的问题是合理的。,三 基本原理-多边形单元离散元法,基本思路在解决连续介质力学问题时，除了边界条件外，还有3个方程必须满足：平衡方程、变形协调方程（弹性力学中的几何方程）、本构方程（弹性力学中的物理方程）；对于离散单元法而言，由于介质被假定为离散体（颗粒或块
27、体）的集合，所以离散体单元间不存在变形协调的约束，但是其运动必须遵循经典力学定律；,三 基本原理-多边形单元离散元法,运动描述 处于一个理想散体中的任意一个多边形单元,具有3个自由度,2个平动自由度与一个转动自由度,可通过Newton第二定律分别描述。,三 基本原理-多边形单元离散元法,运动描述运动方程（平动和转动）：,三 基本原理-多边形单元离散元法,式中，为质量矩阵；,为第j个单元给第i个单元的总作用力。,接触模型Cundall模型（弹簧阻尼器模型）：,三 基本原理-多边形单元离散元法,与颗粒元相同，不再重复介绍。,接触模型简化的解析模型：,基本原理-多边形单元离散元法,式中，与 分别为法
28、向与切向阻尼系数，为静摩擦系数，为材料的弹性模量，为重叠部分的面积，为特征长度。,应用举例,基本原理-多边形单元离散元法,应用举例,基本原理-多边形单元离散元法,应用举例UDEC,基本原理-多边形单元离散元法,综述问题的复杂性：计算量大；接触的力学模型块体元的计算机生成；模拟结果的可视化；但更为复杂的是块体间的接触发现算法。,基本原理-多面体单元离散元法,块体接触的直接判别 接触方式：判别块体间是否存在接触的最简单方法是检查接触发生的所有可能性。对于三维块体的接触有很多方式：点-点，点-边，点-面，边-边，边-面，面-面。,基本原理-多面体单元离散元法,块体接触的直接判别 直接判别的计算效率举
29、例：,基本原理-多面体单元离散元法,如果块体 有 个顶点、条边、个面，第二个块体 有 个顶点、条边、个面，则用直接法判别接触存在的计算次数为 用直接法，两个立方体间存在676种接触的能性。,块体接触的直接判别 直接判别计算效率的改善：,基本原理-多面体单元离散元法,事实上，并不需要如此多次判别，因为有些接触类型可以并入到其他类型中，这样只需要检测点-面、边-边接触两种类型即可，而其他类型可用下面的方法通过点-面、边-边接触两种类型来体现，即（1）当在同样位置存在三个或更多个点-面接触时，说明在该位置存在点-点接触；（2）当两个点-面接触在同一位置同时存在时，表明在该位置存在点-边接触；（3）当
30、两个块体间存在两个边-边接触时，表明存在边-面接触；（4）当三个或更多个边-边接触存在或三个或更多个点-面接触存在时，表明两个块体间存在面-面接触。,块体接触的直接判别 直接判别计算效率的改善：,基本原理-多面体单元离散元法,即使按照上面的方法进行类型合并，接触判别的次数仍然有 对于两个立方体的接触，判别次数为240。,块体接触的直接判别 分析：,基本原理-多面体单元离散元法,接触检验的次数直接取决于所要判别块体的边、顶点与面的平均数量；模拟中不仅要确定块体间是否接触，还要确定接触时的详细信息，如：侵入深度、接触面的法向向量以及接触点等；对某些类型的接触检测是很困难的，例如，在点-面接触检测中
31、，不仅要检查点位于该面的上方或下方，还要检验点是否位于该面投影的边界内，而这在数值计算中并不是通过简单的判别就能实现的。,公共面法 基本思路：,基本原理-多面体单元离散元法,构造一个公共面，通过公共面把两个块体所占据的空间分为两部分；分别检验每个块体与公共面的接触情况。公共面的构造方法可以用一个悬在两个未接触块体间的金属盘来说明，当两个块体未接触时，金属盘与任何一个块体都不接触，但是，随着两个块体逐步靠近直至接触时，金属盘在两个块体的作用下发生扭转直至完全被两个块体夹紧。无论这两个块体的形状和位向如何，金属盘被夹紧后，总会在一个特定位置达到稳定，而金属盘达到稳定的位置恰恰就是处于接触中的两个块
32、体的接触面。,公共面法 基本思路：,基本原理-多面体单元离散元法,进一步对金属盘与两个块体间的相对位置进行分析，当两个块体逐步靠近但还没有接触时，金属盘在块体的推力作用下发生移动和扭转，这时，金属盘总是位于两个块体中间的某个位置，这样就可以很容易地通过把两个块体到金属盘的距离相加而求得两个块体间的空隙尺寸。,公共面法 优点：,基本原理-多面体单元离散元法,通过金属盘例子，块体间接触的很多问题都可以得到简化。只需要对点-面接触进行判别（通过点积），因为块体（或子块体）都是凸多面体，面、边接触都可以通过点-面接触的数量来体现；检测的次数线性地依赖于顶点的数量（直接法的检测次数与顶点数为二次方关系）
33、。可以分别检测块体、的顶点与公共面的接触情况，检测次数为 对于两个立方体的情况，公共面法的检测次数为16，而直接法为240。,公共面法 优点：,基本原理-多面体单元离散元法,对于点-面接触类型，没有必要检测该接触是否位于面的周边以内，因为公共面的位置在不断变化（见后面叙述），如果两个块体都与公共面接触，那么这两个块体必然接触，如果两个块体没有接触，则肯定与公共面不接触；公共面的法向矢量就是接触的法向矢量，不需要额外进行计算；既然公共面的法向是唯一的，那么就可以排除接触法向矢量的不连续变化。公共面的法向矢量可能化发生迅速变化（如点-点接触），但是不会因为接触类型的改变而发生跳跃式变化。可以很容易
34、地确定两个未接触块体间的最小空隙：只需要把两个块体距公共面的距离相加即可。,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,公共面的定位和移动只是一个几何问题，但是与力学计算一样，需要在每个时间步长上对其进行更新。该算法可以用一句话来描述：“把公共面与最接近公共面的顶点间的间隙最小化”。公共面的任何平移和旋转都会减小与最接近顶点间的距离或保持不变。对于已经重叠的块体，仍旧可以采用该方法，但这时的空隙和最近距离用负值表示，此时该算法可描述为：“把公共面与最接近顶点间的重叠值最大化”。,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,1、并把公共面放置在两个块体质心
35、连线的中点，其法向向量为从一个块体的质心指向另一个块体的质心，即,，,式中，；为公共面的单位法向量；为公共面上的参考点；与 分别为块体 与 质心的位矢。,该算法通过平移或转动公共面以使公共面与块体间的间隙最小或最大化。参考点 的作用包含两个方面：（a）公共面绕该点转动；（b）两个块体接触时的法向和切向力的作用点。,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,块体间的公共面,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,2、公共面的平动,公共面的平动可以分解为沿公共面法向和切向的平动。沿公共面法向和切向平动的步骤为：（1）首先确定块体顶点与公共面的最小距离,
36、公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,2、公共面的平动,（2）判别块体间的间隙 的值，如果，说明块体间发生接触，则把参考点 移至否则，把参考点移动到两个最邻近顶点的中间位置,式中，与 分别为块体、上最接近公共面的顶点。,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,3、公共面的转动,公共面经历过平动后，如果两个块体存在接触，那么该面上的参考点 就是接触力的作用点。但是接触力的法向，即公共面的法向矢量，不一定就是平动后的公共面法向矢量，因为平动后的公共面不一定满足裂隙最大化的条件。所以，必须旋转公共面才能满足该要求。公共面的平动只通过一个步骤即可完成，
37、但是其旋转必须要经过迭代才能确定。因为块体上最接近公共面的顶点将随公共面的旋转而改变。,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,3、公共面的转动,过公共面上的参考点 在当前公共面内构造两个方向任意但互相垂直的矢量作为转轴，然后，令公共面的法向量绕两轴发生微小转动，转动角正负各一次，迭代的每一步令公共面的法向矢量转动四次。如果 与 为所构造的两个互相正交的向量，则四次转动为式中，是控制摄动大小的参数。,公共面法 确定公共面位置的算法：,基本原理-多面体单元离散元法,3、公共面的转动,确定公共面转动的迭代算法（其中 一般取为50，一般取为0.010,计算循环,基本原理-多面
38、体单元离散元法,仍旧采用显式差分策略求解块体单元运动的动态过程。在每个时间步长上，用运动方程来描述块体的运动。在把块体作为刚体时，本构方程用接触力和侵入深度间的关系来表示。通过对运动方程的积分来更新块体的位置、速度及侵入深度增量；然后，通过接触力与侵入深度的关系来确定新的接触力，在新接触力的基础上开始下一步长的循环。,计算循环,接触力的计算（本构关系）相对速度,基本原理-多面体单元离散元法,把在接触点处块体B相对于块体A的速度定义为接触速度，可由下式计算,式中：和 分别为块 和块 质心的位置矢量；和 分别为块 和块 质心的速度矢量；和 则分别为块 和块 的角速度矢量；为置换张量；下标 取值范围
39、为13，并用其表明矢量或张量在全局坐标系中的分量（下标符合爱因斯坦求和约定）。,接触力的计算（本构关系）侵入深度,基本原理-多面体单元离散元法,侵入深度增量矢量可以表示为,可以分解为垂直于公共面的法向和沿着公共面的切向分量。法向侵入深度增量分量为,切向侵入深度增量矢量分量为,接触力的计算（本构关系）法向和切向力,基本原理-多面体单元离散元法,法向力增量以压力为正，表示为,切向力增量矢量为,式中，和 分别为法向和切向接触刚度，其量纲为；为接触面积。这样，可以通过接触力和接触面积直接计算接触力。,接触力的计算（本构关系）法向和切向力,基本原理-多面体单元离散元法,对于面-面接触类型接触面积可以直接
40、计算，只需计算公共接触面面积。对于非面-面接触类型的接触面积计算则非常困难，在当前的研究中，为了减少对接触刚度临界值的限制，在上述接触力计算式中假设接触面积值很小，一般可取为处于接触中的两个块体表面积的平均值的百分之一，或由用户自己定义。,接触力的计算（本构关系）法向和切向力,基本原理-多面体单元离散元法,得到法向和切向接触力增量后，就可以计算总法向力和切向力，即,接触力的计算（本构关系）法向和切向力,基本原理-多面体单元离散元法,同时，根据Coulomb准则对法向和切向力进行调整，如果法向力超过了节理的抗拉强度，即，那么法向力和切向力为0，否则，计算最大切向力,然后计算切向力的绝对值,如果切
41、向力的绝对值大于，即，那么切向力被减小为一个限制值，即,块体所受合力（包括力矩）块体质心受力,基本原理-多面体单元离散元法,块体A对块体B作用的接触力矢量可以表示为,这样，就可以更新块体A所受的合力和力矩，即,同理，对于块体B有,块体所受合力（包括力矩）块体的位置更新,基本原理-多面体单元离散元法,单个块体的平动方程可以表示为 式中，与 分别为块体质心的加速度和速度，为与块体质量成比例的粘性阻尼常数，为施加在块体上的合力，包括接触力和外力，为块体质量，为重力加速度矢量。,块体所受合力（包括力矩）块体的位置更新,基本原理-多面体单元离散元法,用Euler方程描述块体的无阻尼刚体转动，这里的转动是
42、指针对块体主惯性轴的转动，即,式中，、与 分别为块体的主惯性矩，、与 分别为关于主轴的角速度，、与 分别为相对于块体主轴的力矩。,运动方程的简化,基本原理-多面体单元离散元法,刚性块体模型更适合于描述准静态问题，在准静态问题中，运动方程中的转动方程可以被简化。由于在准静态问题中块体运动速度很小，所以上述转动方程中的非线性项可以被忽略，从而使方程解耦。同时，由于惯性力与块体所受合力比较而言非常小，所以惯性力张量对转动的影响并不是必须考虑的，所以，只考虑了一个近似的惯性矩，在块体的质心到各个顶点的平均距离基础上计算。这样，增加一个粘性阻尼项转动方程简化为,运动方程的求解,基本原理-多面体单元离散元
43、法,采用中心差分求解运动方程，下面的表达式为 t时刻的平动和转动速度,加速度为,把上述四式代入到平动和简化后的转动方程中，得到,式中,运动方程的求解,基本原理-多面体单元离散元法,平动和转动位移增量为,块体质心坐标被更新为,块体的顶点坐标被更新为,时间步长的确定,基本原理-多面体单元离散元法,同样采用与单自由度体系类比的方法确定临界时间步长的估计值。可以采用块体体系中的最小块体质量（）与最大法向或切向刚度（）来计算临界时间步长的估计值,式中，是用户自行定义的参数，它的物理意义是：上述时间步长估计式是由单自由度体系得来的，事实上，一个块体可能同时与几个块体接触，所以可通过参数 对其进行修正。一般情况下，取0.1可以确保数值计算的稳定性。,
