.. .
。
2.4 接触模型
接触模型是颗粒离散元法的核心,干颗粒模 型是接触的两圆球间在法一切向相对运动时接
触力和局部变形的拟静态关系有许多研究者仍 采用弹簧-阻尼器模型,根据经验或实验给定参数, 但不少三维研究已采用球体接触力学解析结果。 湿颗粒的接触模型是一种近似模型,依据是两刚 性圆球间存在液桥或浸渍态时,当两球对心相对 运动时的液桥表面张力和流体黏性产生的挤压 力,以及切向相对运动时的阻力。无论是干颗粒 还是湿颗粒模型,接触力和变形关系都是非线性 的 , 法 - 切向作用很难分开 , 可是单一的法向或切 向作用求解已很复杂,故近似采用叠加原理
图 5 链式结构示意图
数值计算部分所采用的结构化程序设计的 主流程图如图 6 所示。
图 4 区域法示意图
下，模拟的结果有可能偏离实际很大。加 强离散元法基础理论、基础算法及误差分 析方面的研究,并汲取有限元法等数值方 法的优点,使之既能保持在描述散体的整 体力学行为和力学演化全过程方面的优势, 又能有效描述介质局部连续处应力状态和 变形状态,使离散元法的模型建立真正满 足几何仿真,物理本构仿真,受力仿真和过 程仿真的原则,是离散元法研究领域的首 要工作。 另外,通过同实验结果、 理论解及 其它数值方法的计算结果进行比较,把握 离散元法的计算精度和计算效率,进而对 离散元法的建模和算法进行改进也是必不 可少的 [ 6 ] 。
力 学 （ molecular dynamics) [1] .1971 年 Cundall 提出了适于岩石力学的离散元 法， 而后他与 Stack 在 1979 年联合提出了 适于土力学的离散元法，并推出了二维圆 盘 （disc）程序 BALL 和三维圆球程序 TRUBAL（后发展成商业软件 PFC-2D/3D)， 形成较系统的模型与方法 [ 2 ] 。离散元分为 两个大的分支：以块体为基本单元的块体 离散元法和以圆盘为基本单元的颗粒离散 元法。其中颗粒离散元法是基于最初的圆 盘和圆球颗粒模型发展起来的，适用于颗 粒数目多且单个几何形状可用圆球近似而 不产生显著差异的情况，圆盘属于固体力 学平面问题。 1989 年英国 Aston 大学 Thornton 引 入 Cundall 的 TRUBAL 程序， 从发展颗粒接 触模型入手对程序进行了全面改造形成了 TRUBAL-Aston 版，后定名 GRANULE。它完 全符合弹塑性圆球接触力学原理，能模拟 干-湿、弹性-塑性和颗粒两相流问题。 Leeds 大学等校也利用它用于模拟。在英 国 DEM 研究较深入的还有 Surrey 大学的 Tuzun 研究组（以 DEM 模拟和实验研究见 长） ， Leeds 大 学 的 Ghadiri 研 究 组 ， Swansea 大学 Owen 的研究中心 （以有限元 -离散元法结合见长） 等。 在英国多次举办 相关主题的学术会议，促进了颗粒离散元
[8]
。
2.5 颗粒离散元程序实现
颗粒离散元数值方法需要将研究区域划分 成若干个单元,而且在计算过程中,每个单元除了 有自己的编号!质量及初始形心坐标等固定属性, 还有力 ! 速度 !加速度 ! 形心坐标及位移等多个随 时间变化的物理量,要保存的数据非常多,如何有 效的组织数据,使之既能够调用方便快捷,提高运 行速度,又能够节省内存空间,是利用该数值方法 解决问题的关键. 另外,颗粒离散元法是分时步进行计算的,虽 然力学模型简单,但计算所需机时和内存都很大, 所以目前主要应用在模拟二维问题上。 在三维模 型中,寻找相邻元素花费计算时间较多。为此,提 出了区域法和链式结构法。 区域法的基本思想是将材料所在的坐标系 , 用最大邻居距离 d max , 划分成网格 , 任意元素必 定位于某一网格中。 查找某一元素的邻居元素只 需逐一判断与其相邻的 9 个网格中的每个元素 即可,如图 4
mu （ t ) c u ( t ) ku ( t ) f ( t )
.. .
(2-3)
由于上式中右边的量都是已知的 , 因此可 以求出左边的量 u(t 十 t )。 再将 u(t+ t ) 代入下面两式中 ,就可以得到单元在 t 时 刻的速度 u (t ) 和加速度 u (t ) :
图1
计算循环
滑带土是典型的离散物质，并且目前 对于其在高速剪切条件下的动力学特性， 缺乏试验手段来进行研究。离散单元法作 为离散物质的宏观特性的微观手段，可以 为研究解释滑带剪切过程的动力现象提供 一种研究方法和途径。
2 颗粒离散元法的基本理论
2.1 离散元法的基本原理
离散元法是专门用来解决不连续介质 问题的数值模拟方法。该方法把节理岩体 视为由离散的岩块和岩块间的节理面所组 成, 允许岩块平移、 转动和变形, 而节理 面可被压缩、分离或滑动。因此, 岩体被 看作一种不连续的离散介质。其内部可存 在大位移、旋转和滑动乃至块体的分离 , 从而可以较真实地模拟节理岩体中的非线
3.2 离散元法与其他数值模拟方法 的结合
正因为离散元法有其自身无法克服的缺 陷、问题，所以有时候就必须与其他的数 值方法结合。传统的有限元法、边界元法 等数值模拟方法适合解决连续介质问题， 而离散元法适合于界面弱连接的非连续介 质问题或连续体到非连续踢转化的材料损 伤破坏问题。故将离散元法与有限元法和 边界元法结合起来便能充分发挥各自的优 势， 也可以极大地扩大该数值方法的范围。 胥建龙，唐志平 [ 9 ] 提出并建立了离散元与 离散元法的研究和应用已 40 多年的 历史了，国内外学者发表了大量的学术论 文和研究报告。但是，总体看来，大多数 论文都紧限于利用离散元法计算工程问 题，而对离散元法的理论和算法的研究文 章却很少。然而，离散元法自诞生起就存 在缺乏理础的欠 缺在块体元模型中尤为明显，离散元方法 人为假定太多，法向、切向刚度都是人为 假设的，节理的确定也是经过统计分析处 理的，这也是不真实的，在这些假定前提 有限元结合多尺度方法，将这一方法应用 于激光辐照下预应力铝板的破坏响应，, 从而得到了与实验结果较吻合的计算结 果。 金峰 [10 ] 等提出二维变形体离散元与时 域边界元的耦合模型，从而可以将非连续 体的模拟与无限域的模拟统一在一个模型 中。 雷晓燕 [11] 采用将边界元区域的全部未 知量凝聚到耦合面上的方法来实现有限元 与边界元的耦合，从而达到了节省内存的 目的。
[ 4]
。
性大变形特征。 离散元法的一般求解过程为: 将求解 空间离散为离散元单元阵, 并根据实际问 题用合理的连接元件将相邻两单元连接起 来; 单元间相对位移是基本变量,由力与 相对位移的关系可得到两单元间法向和切 向的作用力; 对单元在各个方向上与其它 单元间的作用力以及其它物理场对单元作 用所引起的外力求合力和合力矩, 根据牛 顿运动第二定律 F=ma 可以求得单元的加 速度; 对其进行时间积分, 进而得到单元 的速度和位移。从而得到所有单元在任意 时刻的速度、 加速度、 角速度、 线位移 和转角等物理量 [5] 。
离散元法研究综述
1
引言
离散元法的思想来源于较早的分子动 法的发展 [ 2 ] 。 20 世纪 90 年代以来，离散元法在国 外得到足够的重视， 并且得到了迅猛发展。 这一时期，各种离散单元法商用软件相继 出现， 美国 ITASCA 公司开发和完善了基于 圆 盘 形 和 球 形 离 散 单 元 的 PFC2.0 和 PFC3.0，这些软件在很多领域得到了广泛 应用，从而使离散元法在工程中的应用向 前迈进了一大步 [1] 。 离散元法在我国的研究和应用起步比 较晚，但是发展迅速。1986 年东北大学的 王泳嘉教授和淮南矿业学院的万禧教授在 第一届全国岩石力学数值计算及模型试验 讨论会上，首次向我国岩石力学和工程力 学界介绍了离散元法的基本原理及几个应 用例子 [3] 。自此，关于离散元理论和应用 方面的研究论文在国内刊物上不断出现。 离散元法是一种专门用于解决非连续 介质问题的有效的方法，其最初的研究对 象主要是岩石等非连续介质的力学行为。 它的基本原理是牛顿第二定律，其基本思 想是将岩体看成是由断层、节理、裂隙等 结构面切割而成的一个个刚性或者可变性 块体，块体与块体之间通过角、面或者边 进行接触， 块体可以平移、 转动或者变形， 节理面可以被压缩、分离、滑动，所有块 体镶嵌排列，在某一时刻当给定块体一个 外力或者边界位移约束，各个块体在外界
[7]
（2-1）
式中 m—单元质量； U—位移； K—刚度系数； C—阻尼系数； f—单元的外荷载。 式（2-1）的动态松弛解法就是假定 t+ t 时 刻 以 前 的 变 量 f（t) , u (t ) ， 利 u (t t ) , u (t t ) 以及 u (t t ) 等已知，
链式结构法思路如下：第一步计算时,采用 全局搜索, 并将第 i 个元素的一级相邻元素用链 表记录下来 ,在接下来的时间步长内 ,当查找第 i 个元素的相邻关系时,只判断上一步所记录的该 元素的链式结构体内的元素即可,元素的链式结 构体的大小取决于所研究的问题和所选时间的 步长。为了防止元素结构体内丢失相邻元素,应 每隔若干步再进行全局搜索一次。同理,全局搜 索的间隔长度,也与所研究的问题和所选时间步 长有关,一般通过实践可以确定。图 5 是这种方 法的示意图,图中 1,2,3 分别表示 i 元素的第一、 二、三级相邻元素链。
图 2 离散元的单元分类
2.3 离散元的计算方法-动态松弛法
离散元法的计算原理虽然很简单，但是 也必须由计算机来实现，所以会遇到很多 问题。离散元法具体的求解过程分为显式 解法和隐式解法 , 显示解法用于动力问题 的求解或动态松弛法困的静力求解 ,而隐 式解法用于求解静力问题的静态松弛法。 其中，动态松弛法是把非线性静力学问题 转化为动力学问题求解的一种数值方法， 其实质是对临界阻尼振动方程进行逐步积 分，通过质量阻尼和刚度阻尼来吸收系统 的动能，收敛于静态值。这种带有阻尼项 的动态平衡方程 , 利用有限差分法按时步 在计算机上迭代求解就是所谓的动态松弛 法。为保证准静界，需要加入质量阻尼和 刚度阻尼来吸收系统的动能，考虑阻尼后 离散单元法的基本方程为：