离散元课件(三)

1、二二 基本原理基本原理-块体元离散元法块体元离散元法 安全与防灾工程研究所安全与防灾工程研究所刘刘 军军 根据离散化模型中所采用的单元种类分别介绍离根据离散化模型中所采用的单元种类分别介绍离散元法的基本原理：散元法的基本原理：颗粒元 二维圆盘单元 三维圆球单元块体元 多边形单元 多面体单元二二 基本原理基本原理基本思路基本思路将所研究的区域或散体体系剖分或看作多边形块体单元的集合；单元间的接触方式包括：角-角接触 、角-边接触、边-边接触三种形式；随着单元的平移和转动，各个单元间的接触关系由运动方程动态调整；最终，块体单元达到平衡状态或一直运动下去；单元性质可以假设为刚性（对于应力水平低的问题

2、是合理的)。 二二 基本原理基本原理-多边形单元离散元法多边形单元离散元法基本思路基本思路在解决连续介质力学问题时，除了边界条件外，还有3个方程必须满足：平衡方程、变形协调方程（弹性力学中的几何方程）、本构方程（弹性力学中的物理方程）；对于离散单元法而言，由于介质被假定为离散体（颗粒或块体）的集合，所以离散体单元间不存在变形协调的约束，但是其运动必须遵循经典力学定律； 二二 基本原理基本原理-多边形单元离散元法多边形单元离散元法运动描述运动描述 处于一个理想散体中的任意一个多边形单元,具有3个自由度,2个平动自由度与一个转动自由度,可通过NewtonNewton第二定律第二定律分别描述。二二

3、基本原理基本原理-多边形单元离散元法多边形单元离散元法运动描述运动描述运动方程（平动和转动）运动方程（平动和转动）：基本原理基本原理-多边形单元离散元法多边形单元离散元法式中， 为质量矩阵；iM为第j个单元给第i个单元的总作用力。ijpF接触模型接触模型CundallCundall模型（弹簧阻尼器模型）模型（弹簧阻尼器模型）：基本原理基本原理-多边形单元离散元法多边形单元离散元法与颗粒元相同，不再重复介绍。 接触模型接触模型简化的解析模型简化的解析模型：基本原理基本原理-多边形单元离散元法多边形单元离散元法式中， 与 分别为法向与切向阻尼系数， 为静摩擦系数， 为材料的弹性模量， 为重叠部分的

4、面积， 为特征长度。 应用举例应用举例 基本原理基本原理-多边形单元离散元法多边形单元离散元法应用举例应用举例 基本原理基本原理-多边形单元离散元法多边形单元离散元法应用举例应用举例UDEC 基本原理基本原理-多边形单元离散元法多边形单元离散元法综述综述问题的复杂性问题的复杂性： 计算量大； 接触的力学模型 块体元的计算机生成； 模拟结果的可视化； 但更为复杂的是块体间的接触发现算法。 基本原理基本原理-多面体单元离散元法多面体单元离散元法块体接触的直接判别块体接触的直接判别 接触方式接触方式：判别块体间是否存在接触的最简单方法是检查接触发生的所有可能性。对于三维块体的接触有很多方式：点-点，

5、点-边，点-面，边-边，边-面，面-面。 基本原理基本原理-多面体单元离散元法多面体单元离散元法块体接触的直接判别块体接触的直接判别 直接判别的计算效率举例直接判别的计算效率举例： 基本原理基本原理-多面体单元离散元法多面体单元离散元法如果块体 有 个顶点、 条边、 个面，第二个块体 有 个顶点、 条边、 个面，则用直接法判别接触存在的计算次数为 用直接法，两个立方体间存在676种接触的能性。 AAvAeAfBBvBfBBBAAAfevfevnBe块体接触的直接判别块体接触的直接判别 直接判别计算效率的改善直接判别计算效率的改善： 基本原理基本原理-多面体单元离散元法多面体单元离散元法事实上，

6、并不需要如此多次判别，因为有些接触类型可以并入到其他类型中，这样只需要检测点-面、边-边接触两种类型即可，而其他类型可用下面的方法通过点-面、边-边接触两种类型来体现，即（1）当在同样位置存在三个或更多个点-面接触时，说明在该位置存在点-点接触；（2）当两个点-面接触在同一位置同时存在时，表明在该位置存在点-边接触；（3）当两个块体间存在两个边-边接触时，表明存在边-面接触；（4）当三个或更多个边-边接触存在或三个或更多个点-面接触存在时，表明两个块体间存在面-面接触。块体接触的直接判别块体接触的直接判别 直接判别计算效率的改善直接判别计算效率的改善： 基本原理基本原理-多面体单元离散元法多面

7、体单元离散元法即使按照上面的方法进行类型合并，接触判别的次数仍然有 对于两个立方体的接触，判别次数为240。ABBAABBAfefefvfvn块体接触的直接判别块体接触的直接判别 分析分析： 基本原理基本原理-多面体单元离散元法多面体单元离散元法接触检验的次数直接取决于所要判别块体的边、顶点与面的平均数接触检验的次数直接取决于所要判别块体的边、顶点与面的平均数量量；模拟中不仅要确定块体间是否接触，还要确定接触时的详细信息详细信息，如：侵入深度、接触面的法向向量以及接触点等；对某些类型的接触检测是很困难的对某些类型的接触检测是很困难的，例如，在点-面接触检测中，不仅要检查点位于该面的上方或下方，

8、还要检验点是否位于该面投影的边界内，而这在数值计算中并不是通过简单的判别就能实现的。 公共面法公共面法 基本思路基本思路： 基本原理基本原理-多面体单元离散元法多面体单元离散元法构造一个公共面，通过公共面把两个块体构造一个公共面，通过公共面把两个块体所占据的空间分为两部分所占据的空间分为两部分；分别检验每个块体与公共面的接触情况分别检验每个块体与公共面的接触情况。公共面的构造方法可以用一个悬在两个未接触块体间的金属盘来说明，当两个块体未接触时，金属盘与任何一个块体都不接触，但是，随着两个块体逐步靠近直至接触时，金属盘在两个块体的作用下发生扭转直至完全被两个块体夹紧。无论这两个块体的形状和位向如

9、何，金属盘被夹紧后，总会在一个特定位置达到稳定，而金属盘达到稳定的位置恰恰就是处于接触中的两个块体的接触面。 公共面法公共面法 基本思路基本思路： 基本原理基本原理-多面体单元离散元法多面体单元离散元法进一步对金属盘与两个块体间的相对位置进行分析，当两个块体逐步靠近但还没有接触时，金属盘在块体的推力作用下发生移动和扭转，这时，金属盘总是位于两个块体中间的某个位置，这样就可以很容易地通过把两个块体到金属盘的距离相加而求得两个块体间的两个块体间的空隙尺寸空隙尺寸。 公共面法公共面法 优点优点： 基本原理基本原理-多面体单元离散元法多面体单元离散元法通过金属盘例子，块体间接触的很多问题都可以得到简化

10、。只需要对点只需要对点- -面接触进行判别面接触进行判别（通过点积），因为块体（或子块体）都是凸多面体，面、边接触都可以通过点-面接触的数量来体现；检测的次数线性地依赖于顶点的数量检测的次数线性地依赖于顶点的数量（直接法的检测次数与顶点数为二次方关系）。可以分别检测块体 、 的顶点与公共面的接触情况，检测次数为 对于两个立方体的情况，公共面法的检测次数为16，而直接法为240。ABBAvvn公共面法公共面法 优点优点： 基本原理基本原理-多面体单元离散元法多面体单元离散元法对于点-面接触类型，没有必要检测该接触是否位于面的周边以内，因为公共面的位置在不断变化（见后面叙述），如果两个块体都与如果

11、两个块体都与公共面接触，那么这两个块体必然接触，如果两个块体没有接触，公共面接触，那么这两个块体必然接触，如果两个块体没有接触，则肯定与公共面不接触则肯定与公共面不接触；公共面的法向矢量就是接触的法向矢量公共面的法向矢量就是接触的法向矢量，不需要额外进行计算；既然公共面的法向是唯一的公共面的法向是唯一的，那么就可以排除接触法向矢量的不连续变化。公共面的法向矢量可能化发生迅速变化（如点-点接触），但是不会因为接触类型的改变而发生跳跃式变化。可以很容易地确定两个未接触块体间的最大空隙可以很容易地确定两个未接触块体间的最大空隙：只需要把两个块体距公共面的距离相加即可。 公共面法公共面法 确定公共面位

12、置的算法确定公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法公共面的定位和移动只是一个几何问题，但是与力学计算一样，需要在每个时间步长上对其进行更新。该算法可以用一句话来描述：“把公共面与最接近公共面的顶点间的间隙最大化”。公共面的任何平移和旋转都会减小与最接近顶点间的距离或保持不变。对于已经重叠的块体，仍旧可以采用该方法，但这时的空隙和最近距离用负值表示，此时该算法可描述为：“把公共面与最接近顶点间的重叠值最大化”。 公共面法公共面法 确定公共面位置的算法确定公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法1、并把公共面放置在两个块体质心连

13、线的中点，其法向向量为从一个块体的质心指向另一个块体的质心，即 ， 2/iiiBACzii/Zn 式中， ； ； 为公共面的单位法向量； 为公共面上的参考点； 与 分别为块体 与 质心的位矢。 iiiABZiizZZ2iniCiAiBAB该算法通过平移或转动公共面以使公共面与块体间的间隙最小或最该算法通过平移或转动公共面以使公共面与块体间的间隙最小或最大化大化。参考点 的作用包含两个方面：（a）公共面绕该点转动；（b）两个块体接触时的法向和切向力的作用点。iC公共面法公共面法 确定公共面位置的算法确定公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法块体间的公共面 公共

14、面法公共面法 确定公共面位置的算法确定公共面位置的算法：基本原理基本原理-多面体单元离散元法多面体单元离散元法2、公共面的平动公共面的平动 公共面的平动可以分解为沿公共面法向和切向的平动。沿公共面法向和切向平动的步骤为：（1）首先确定块体顶点与公共面的最小距离 BiiBdVn min AiiAdVn max公共面法公共面法 确定公共面位置的算法确定公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法2、公共面的平动公共面的平动 （2）判别块体间的间隙 的值，如果 ，说明块体间发生接触，则把参考点 移至否则，把参考点移动到两个最邻近顶点的中间位置ABdd 0ABddiC2/

15、:iBAiiddnCC2/minmaxBiAiivvC式中， 与 分别为块体 、 上最接近公共面的顶点。 maxAivminBivAB公共面法公共面法 确定公共面位置的算法确定公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法3、公共面的转动公共面的转动 公共面经历过平动后，如果两个块体存在接触，那么该面上的参考点 就是接触力的作用点。但是接触力的法向，即公共面的法向矢量，不一定就是平动后的公共面法向矢量，因为平动后的公共面不一定满足裂隙最大化的条件。所以，必须旋转公共面才能满足该要求。公共面的平动只通过一个步骤即可完成，但是其旋转必须要经过迭代才能确定。因为块体上最接

16、近公共面的顶点将随公共面的旋转而改变。 iC公共面法公共面法 确定公共面位置的算法确定公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法3、公共面的转动公共面的转动 过公共面上的参考点 在当前公共面内构造两个方向任意但互相垂直的矢量作为转轴，然后，令公共面的法向量绕两轴发生微小转动，转动角正负各一次，迭代的每一步令公共面的法向矢量转动四次。如果 与 为所构造的两个互相正交的向量，则四次转动为式中， ， 是控制摄动大小的参数。 iCipiqzkiii/:pnnzkiii/:pnnzkiii/:qnnzkiii/:qnn221kzk公共面法公共面法 确定公共面位置的算法确定

17、公共面位置的算法： 基本原理基本原理-多面体单元离散元法多面体单元离散元法3、公共面的转动公共面的转动 确 定 公 共 面转 动 的 迭 代算法（其中 一 般 取为50， 一般取为0.010 maxkmink计算循环计算循环 基本原理基本原理-多面体单元离散元法多面体单元离散元法仍旧采用显式差分策略求解块体单元运动的动态过程。在每个时间步长上，用运动方程来描述块体的运动。在把块体作为刚体时，本构方程用接触力和侵入深度间的关系来表示。通过对运动方程的积分来更新块体的位置、速度及侵入深度增量；然后，通过接触力与侵入深度的关系来确定新的接触力，在新接触力的基础上开始下一步长的循环。计算循环计算循环块

18、体质心块体质心力更新力更新接触力接触力更新更新相对接触相对接触速度更新速度更新块体位置块体位置更新更新接触力的计算（本构关系）接触力的计算（本构关系）相对速度相对速度 基本原理基本原理-多面体单元离散元法多面体单元离散元法把在接触点处块体B相对于块体A的速度定义为接触速度，可由下式计算 式中： 和 分别为块 和块 质心的位置矢量； 和 分别为块 和块 质心的速度矢量； 和 则分别为块 和块 的角速度矢量； 为置换张量；下标 取值范围为13，并用其表明矢量或张量在全局坐标系中的分量（下标符合爱因斯坦求和约定）。 iAiBABAix Bix ABAiBiABijkekji,kkAjijkAikkB

19、jijkBiiACexBCexV接触力的计算（本构关系）接触力的计算（本构关系）侵入深度侵入深度 基本原理基本原理-多面体单元离散元法多面体单元离散元法侵入深度增量矢量可以表示为可以分解为垂直于公共面的法向和沿着公共面的切向分量。法向侵入深度增量分量为 切向侵入深度增量矢量分量为tVUiiiinnUUjijisinnUUU接触力的计算（本构关系）接触力的计算（本构关系）法向和切向力法向和切向力 基本原理基本原理-多面体单元离散元法多面体单元离散元法法向力增量以压力为正，表示为 切向力增量矢量为cnnnAUKFcsissiAUKF式中， 和 分别为法向和切向接触刚度，其量纲为 ； 为接触面积。这

20、样，可以通过接触力和接触面积直接计算接触力。 nKsKmPa/cA接触力的计算（本构关系）接触力的计算（本构关系）法向和切向力法向和切向力 基本原理基本原理-多面体单元离散元法多面体单元离散元法对于面-面接触类型接触面积可以直接计算，只需计算公共接触面面积。对于非面-面接触类型的接触面积计算则非常困难，在当前的研究中，为了减少对接触刚度临界值的限制，在上述接触力计算式中假设接触面积值很小，一般可取为处于接触中的两个块体表面积的平均值的百分之一，或由用户自己定义。 cnnnAUKFcsissiAUKF接触力的计算（本构关系）接触力的计算（本构关系）法向和切向力法向和切向力 基本原理基本原理-多面

21、体单元离散元法多面体单元离散元法得到法向和切向接触力增量后，就可以计算总法向力和切向力，即 cnnnAUKFcsissiAUKFnnnFFF :sisisiFFF :接触力的计算（本构关系）接触力的计算（本构关系）法向和切向力法向和切向力 基本原理基本原理-多面体单元离散元法多面体单元离散元法同时，根据Coulomb准则对法向和切向力进行调整，如果法向力超过了节理的抗拉强度，即 ，那么法向力和切向力为0，否则，计算最大切向力cnTAFtanmaxnsFcAF然后计算切向力的绝对值 2/1sisisFFF如果切向力的绝对值大于 ，即 ，那么切向力被减小为一个限制值，即 sFmaxssFFmaxs

22、ssisiFFFF/:max块体所受合力（包括力矩）块体所受合力（包括力矩）块体质心受力块体质心受力 基本原理基本原理-多面体单元离散元法多面体单元离散元法块体A对块体B作用的接触力矢量可以表示为这样，就可以更新块体A所受的合力和力矩，即 同理，对于块体B有 siiniFnFFiAiAiFFF :kjjijkAiAiFACeMM :iBiBiFFF :kjjijkBiBiFBCeMM :块体所受合力（包括力矩）块体所受合力（包括力矩）块体的位置更新块体的位置更新 基本原理基本原理-多面体单元离散元法多面体单元离散元法单个块体的平动方程可以表示为 式中， 与 分别为块体质心的加速度和速度， 为与

23、块体质量成比例的粘性阻尼常数， 为施加在块体上的合力，包括接触力和外力， 为块体质量， 为重力加速度矢量。 iiiigmFxx/ ix ix iFmig块体所受合力（包括力矩）块体所受合力（包括力矩）块体的位置更新块体的位置更新 基本原理基本原理-多面体单元离散元法多面体单元离散元法用Euler方程描述块体的无阻尼刚体转动，这里的转动是指针对块体主惯性轴的转动，即 1232311MIII2313122MIII3121233MIII式中， 、 与 分别为块体的主惯性矩， 、 与 分别为关于主轴的角速度， 、 与 分别为相对于块体主轴的力矩。 1I2I3I1231M2M3M运动方程的简化运动方程的简化 基本原理基本原理-多面体单元离散元法多面体单元离散元法刚性块体模型更适合于描述准静态问题，在准静态问题中，运动方程中的转动方程可以被简化。由于在准静态问题中块体运动速度很小，