离散元课件(二forcopy)

1、二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 岩土工程研究所岩土工程研究所 刘刘 军军 根据离散化模型中所采用的单元种类分别介绍离根据离散化模型中所采用的单元种类分别介绍离 散元法的基本原理：散元法的基本原理： 颗粒元 二维圆盘单元 三维圆球单元 块体元 多边形单元 多面体单元 二二 基本原理基本原理 基本假设 假定速度和加速度在每个时间步长内为常量 ； 选取的时间步长应该足够小以至于在单个时间步长 内扰动的传播不会超过当前与之相邻的粒子 。 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 运动描述运动描述 处于一个理想散体中的任意一个颗粒,具有6个自由 度,3个

2、平动自由度与三个转动自由度,可通过NewtonNewton 第二定律第二定律分别描述。 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 运动描述运动描述 平动方程平动方程： 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 gFF V i k j ijdijc i i m t m i 1 , d d 式中， 与 分别为颗粒 的质量和速度。 为时间， 为颗粒的重力， 与 分别为颗粒 与 的 接触力与粘性接触阻尼力， 为所有与颗粒接触的颗 粒总数。 i m i Vit g i m ijc, F ijd, F ij i k 运动描述运动描述 接触力的分解接触力的分解： 颗粒

3、与 间的接触力可分解为法向与切向接触力，即 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 ji ijctijcnijc, FFF 同理，粘性接触阻尼力也可分解为法向与切向分量形 式，即 ijdtijdnijd, FFF 运动描述运动描述 接触力产生的力矩接触力产生的力矩： 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 颗粒间的接触力作用在两个颗粒的接触点上，而不是作 用在颗粒的中心，所以这些接触力（除法向接触力 外）将会对颗粒产生力矩 ， ijcn, F i T ijdtijctii, FFRT 式中， 为从颗粒 的质心指向接触点的矢量，其幅值为 （颗粒的半径）。

4、i Ri i R 运动描述运动描述 转动方程转动方程： 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 转动方程可以表示为 i k j i i i t I 1 d d T 式中， 与 分别为颗粒 的转动惯量与角速度，对于 球形颗粒 为 i I i i 2 5 2 iii RmI i I 接触模型接触模型 综述综述： 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 关于接触力的计算模型已有大量的研究成果，目前仍旧是 一个活跃的研究领域，特别是对于切向力的计算方法。 对于理想散体颗粒（无粘连）：对于理想散体颗粒（无粘连）：采用Hertz理论理论描述法向作用，而采 用Min

5、dlin与Deresiewicz理论描述切向作用； 对于存在粘连的散体颗粒：对于存在粘连的散体颗粒：法向接触力法向接触力根据在Hertz理论理论基础上考虑 粘连力的JKR（Johnson-Kendall-Roberts）理论确定，切向接触力切向接触力增量 则根据把Savkoor和Briggs理论与Mindlin和Deresiewicz理论相结合形 成的理论确定。 接触模型接触模型 两个处于接触颗粒单位法向和切向向量两个处于接触颗粒单位法向和切向向量： 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 ii R/R n单位法向向量 单位切向向量 nn nn t ijij ijij VV

6、 VV 单位切向量之所以通过两个颗粒的相对速度来计算， 是因为接触力与粘性阻尼力的方向与相对速度的方 向相同。 接触模型接触模型 两个处于接触颗粒接触点的相对速度两个处于接触颗粒接触点的相对速度： 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 法向相对速度为 iijj RRVVV ijij nn ijijn, VV 切向相对速度为 nn ijij VVV ijt, nn ijijt, VV 或者写为 接触模型接触模型 法向接触力计算模型法向接触力计算模型 HertzHertz模型模型： 2 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 为颗粒i与j接触时的侵入深度 n

7、2/3* , 3 4 nijcn RE F 式中 )1 (2 2 * v E E ji RR R 11 * n ijjin RRRR 接触模型接触模型 法向接触力计算模型法向接触力计算模型 CundallCundall模型模型： 2基本原理基本原理- 球形颗粒元离散元法球形颗粒元离散元法 n nnijcn k , F 式中， 为法向弹簧刚度。 n k 接触模型接触模型 法向接触力计算模型法向接触力计算模型 法向粘性接触阻尼力法向粘性接触阻尼力： 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 式中， 为法向粘性接触阻 尼系数。 nn ij VF nijdn c , n c 接触模

8、型接触模型 切向接触力计算模型切向接触力计算模型 综述综述： 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 处于接触中的两个颗粒的切向作用，从本质上讲，是 一种摩擦行为，按照摩擦机理，摩擦力包括：滑动摩 擦、滚动摩擦与静摩擦，其中滑动摩擦与静摩擦属于 切向摩擦力；滚动摩擦是由于法向接触应力的不均匀 分布产生的。介绍两个切向接触力模型： Coulomb准则准则 Mindlin与与Deresiewicz切向接触力模型切向接触力模型 接触模型接触模型 切向接触力计算模型切向接触力计算模型 Coulomb准则准则： 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 在离散元

9、模拟中，一般用Coulomb准则这种简单的形 式描述，静摩擦的详细刻画需要涉及切向位移甚至可能 要考虑时间依赖效应。 式中， 为静摩擦系数，切向摩擦力的方向为与相对滑 动的趋势相反。 s ijcnsijctijcns ijcnsijctijct ijct , , , , , FFF FFF F 接触模型接触模型 切向接触力模型切向接触力模型 Mindlin与与Deresiewicz模型模型 ： 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 t maxt, maxt,t t ijcns ijct F F 2 3 , , ,min 11 式中 为颗粒 与 间的累积切向位移矢量 t

10、i j nsmaxt, v v 12 2 t N t N t nn ij V 1 接触模型接触模型 切向接触力模型切向接触力模型 阻尼力阻尼力 ： 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 式中， 为切向粘性接触阻尼系数。 nn ij VF tijdt c , t c 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 计算模型总结计算模型总结 运动方程运动方程 接触力的计算接触力的计算 法向接触力法向接触力 t maxt, maxt,t t ijcns ijct F F 2 3 , , ,min 11 gFF V i k j ijdijc i i m t m i

11、 1 , d d i k j i i i t I 1 d d T n 2/3* , 3 4 nijcn RE F n nnijcn k , F ijcnsijctijcns ijcnsijctijct ijct , , , , , FFF FFF F nn ij VF nijdn c , nn ij VF tijdt c , 切向接触力切向接触力 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 运动方程的求解运动方程的求解 gFF V i k j ijdijc i i m t m i 1 , d d i k j i i i t I 1 d d T 一般采用两种方法求解运动方程：一

12、般采用两种方法求解运动方程： 中心差分法中心差分法 Verlet积分法积分法 hXXX n nn 2 1 2 1 hXXX n nn 2 11 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 运动方程的求解运动方程的求解 中心差分法中心差分法 运动方程可由Verlet显式积分求解。通过积分可获得粒 子的新位置，积分时需要粒子的当前及上一步长的位 置数据，而不需要粒子的速度数据。 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 运动方程的求解运动方程的求解 Verlet积分法积分法 h XX X nn n 2 11 h XX X nn n 2 11 2 22 4 2

13、h XXX nnn 2 11 2 h XXX nnn 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 运动方程的求解运动方程的求解 Verlet积分法积分法 定义定义： 从而导出Verlet方程为 2 11 2hXXXX nnnn 接触发现算法接触发现算法 在一个由众多颗粒组成的体系中，直接判别颗粒是否接触需 要耗费大量的计算时间，因而，为了节约计算时间，提高计算效 率，一般不直接判别任意两个颗粒间是否存在接触，而是分两个 步骤判别颗粒间的接触是否存在：首先，对一个颗粒，判别其潜 在的邻居个数，然后，准确确定该颗粒与每个邻居是否接触。虽 然在确定邻居数目时也要耗费一定的计算时间，

14、但是仍旧比逐个 准确判别颗粒间接触是否存在要节约时间。因而，接触发现算法 的效率在多颗粒体系力学行为模拟中至关重要。 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 接触发现算法接触发现算法 介绍三种针对球形颗粒的接触发现算法： Verlet邻居目录邻居目录法法 连接单元法连接单元法 边界盒法边界盒法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 接触发现算法接触发现算法 Verlet邻居目录邻居目录法法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 当需要判别体系中某个颗粒的邻居数量时，在该粒子周围构 建一个球（称之为参考球，称该颗粒为核心颗粒）

15、，参考球半径 为体系中最大粒子半径的若干倍，那么参考球所包围的所有粒子 为该球中心粒子的邻居。参考球半径的选取取决于粒子的运动速 度及体系中粒子的密度。对于每个粒子，都可生成一个邻居粒子 的目录。为了得到邻居目录，对每一个粒子而言，所有标号大于 该粒子的粒子都必须被检验，判断是否位于该粒子的参考球中， 而对于标号小于该粒子标号的粒子则没有必要被检验，因为邻居 是互相的，没有必要对一个邻居对检验两次。 Verlet邻居目录邻居目录法法及粒子存储目录及粒子存储目录 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 接触发现算法接触发现算法 Verlet邻居目录邻居目录法法 二二 基本原

16、理基本原理-球形颗粒元离散元法球形颗粒元离散元法 接触发现算法接触发现算法 Verlet邻居目录邻居目录法法 对于n个颗粒组成的体系，用Verlet邻居目录法需要 次 计算，也就是说计算次数仍旧为 量级。然而，并不需要在每 个时间步长上都对邻居目录进行更新。更新的频率取决于体系中 粒子的密度、粒子的运动速度以及参考球的尺寸。参考球的半径 也可以根据颗粒体系的稠密程度及运动速度进行动态调整，并且 参考球半径与邻居目录的更新频率呈反比关系，参考球半径越小， 邻居目录的更新频率越高；但是参考球半径越大，则有更多的粒 子位于球体内，所以判别是否为邻居就需要较长的时间。 2/1nn 2 no 接触发现算

17、法接触发现算法 连接单元法连接单元法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 将颗粒体系所占据的空间划分成规则的网格，对于三维问题，可以划 分为 个立方体单元，对于二维体系，则划分为 个正方 形单元，对于颗粒体系所占据空间形状不规则时，也可采用其他形状 单元划分。但是所有单元的尺寸必须大于粒子的尺寸。与Verlet邻居目 录法的主要区别在于：相邻单元法中的单元不依附于粒子，单元不随 粒子的运动而运动。如果粒子就当前的位置被分配到某个单元，显然， 只有在同一个单元或直接相邻单元内的粒子间才可能发生相互作用， 也就是说，只有相邻单元内的粒子才能成为邻居。 mmmmm 接触发现

18、算法接触发现算法 连接单元法连接单元法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 例如，对一个二维体系而言，只有在9个不同的单元内可能包含邻居 粒子，对于三维体系，则只有在27个不同单元内可能包含邻居粒子。 与前面介绍的Verlet邻居目录法相同，每个粒子对只需要检验一次， 这样，没必要对9个单元都进行检验，只需检验中心单元及邻居单元 的一半即可，即在二维情况下，只需检验5个单元，在三维情况下， 只需检验14个单元。 接触发现算法接触发现算法 连接单元法连接单元法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 2D体系中的连接单元法体系中的连接单元法 接触发

19、现算法接触发现算法 边界盒法边界盒法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 这个方法与前述两种方法不同。首先，在每一个粒子周围构建一个 边界盒，边界盒的尺寸按这样的方式选取：使每个粒子刚好放进它使每个粒子刚好放进它 的边界盒内的边界盒内。边界盒的边为直线，并且与体系的坐标轴平行。在判 别颗粒的邻居时，把边界盒投影到体系的坐标轴上。通过边界盒在 坐标轴上投影的起点和终点来判别是否为邻居。 接触发现算法接触发现算法 边界盒法边界盒法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 a a 围绕围绕 在每在每 个粒个粒 子周子周 围的围的 边界边界 盒盒 b b

20、 两个不两个不 同时刻同时刻 粒子的粒子的 边界盒边界盒 在轴上在轴上 的投影的投影 在判别颗粒的邻居时，把边界 盒投影到体系的坐标轴上。例 如，图a表明了粒子及边界盒的 位置，图b为图a中的边界盒在体 系x轴的投影。通过边界盒在坐 标轴上投影的起点和终点来判 别是否为邻居，出于这个原因， 投影的起点和终点序列被存储 在目录中。 接触发现算法接触发现算法 边界盒法边界盒法 二二 基本原理基本原理-球形颗粒元离散元法球形颗粒元离散元法 对于3D体系，必须把边界盒子在三个坐标轴上投影，所以生成三 个目录。如果一个粒子的边界盒在某个坐标轴上投影的起点和终 点间包含另外一个粒子边界盒投影的起点、终点或起点和终点， 就说明这两个粒子的边界盒在该坐标轴上的投影发生了重叠。如 果两个边界盒的投影在每一个坐标轴上都发生重叠，那么就说明