离散元法研究综述

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