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