第1章 离散单元法及其应用概述
1.1 离散单元法简述
自然界的宏观物质均由一系列细微观粒子构成，当不考虑单一物质颗粒在外力作用下的运动和变形特性对材料宏观力学行为的影响时，一般将研究对象抽象为连续体并采用连续介质力学方法进行研究。连续介质力学是近代固体力学的理论基础，通过建立各种物质的力学模型，把各种物质的本构关系用数学形式确定下来，在给定的初始条件和边界条件下求出问题的解。
常用的连续介质分析方法包括有限元法、有限差分法、边界元法等。在开展岩土工程数值模拟时，连续介质分析方法具有计算效率高、可构建复杂模型等优点，但同时也存在诸多缺陷，如不能反映岩土材料细微观结构之间的复杂相互作用，无法再现岩土材料的破裂孕育演化过程。在这一背景下，离散单元法应运而生。
离散单元法(distinct element method，DEM)是由Cundall［1］在1971年基于分子动力学原理提出的一种离散体物料的分析方法。离散单元法的基本思想是将求解空间离散划分成若干个块体单元或者颗粒单元，并定义单元之间存在接触作用，根据力-位移法则和牛顿第二定律建立各单元之间的运动方程，采用时步迭代的方法进行求解，从而求取“非连续体”的运动形态。该方法是继连续介质力学方法后，用于分析岩土力学问题的又一种强有力的数值计算方法。
离散单元法最早应用于具有裂隙、节理的岩体问题研究，将岩体视为被裂隙、节理切割的若干块体的组合体，基于岩体的变形主要依赖于软弱结构面(裂隙、节理等)的客观事实，提出了将岩块假定为刚体，以刚性元及其边界的几何方程、运动方程和本构方程为基础，采用动态松弛迭代格式，建立节理岩体非连续介质大变形的差分方程并进行求解。根据所采用的求解算法，离散单元法分为隐式离散单元法和显式离散单元法。
根据离散体自身的几何特征，可分为块体和颗粒两大分支。Cundall等［2，3］改进了最初的刚体离散元模型，融合了岩块自身变形，开发了可变形的块体模型通用离散单元法程序(universal distinct element code，UDEC)，并将其推广应用至模拟爆炸运动以及岩块破碎的过程等。在此同时，Cundall等［4，5］开发了二维圆形颗粒(ball)软件，用于研究颗粒介质的物理力学行为，其结果与其他学者采用光弹技术获得的试验结果较为吻合，为研究颗粒散体介质材料的力学行为开辟了新的途径。
20世纪90年代以来，基于离散元理论开发的商业软件和开源软件发展迅速，其中以美国依泰斯卡(Itasca)公司和英国德颐姆方案(DEM-Solutions)公司开发的系列软件最具特色且应用最为广泛。美国依泰斯卡公司以解决岩土工程问题为目标，旗下离散元软件包括基于不规则形状块体单元的通用离散单元法程序和三维离散单元法程序(3 dimension distinct element code，3DEC)软件，以及基于圆盘颗粒单元的二维颗粒流(particle flow code 2 dimension，PFC2D)和基于球形颗粒单元的三维颗粒流(particle flow code 3 dimension，PFC3D)软件。英国德颐姆方案公司以颗粒处理和生产操作为目标，开发了颗粒流软件EDEM，通过模拟散体物料加工处理过程中颗粒体系的行为特征，协助设计人员对各类散料处理设备进行设计、测试和优化。同时，中国科学院基于连续介质力学的离散单元法(continuum-based discrete element method，CDEM)开发的力学分析系列软件GDEM、英国洛克菲尔德(Rockfield)公司开发的有限元/离散元耦合软件ELFEN、加拿大多伦多大学基于有限元/离散元耦合方法开发的地质力学软件Y-Geo、石根华建立的非连续变形分析(discontinous deformation analysis，DDA)方法、南京大学开发的矩阵离散元软件(fast GPU matrix computing of discrete element method，MatDEM)、Olivier和Janek采用C++和Python语言编写的开源离散元软件YADE等也得到了较为广泛的应用。
目前，基于离散元理论开发的软件为解决众多涉及颗粒、结构、流体与电磁及其耦合等综合问题提供了有效的平台，已成为过程分析、设计优化和产品研发的有力工具。其中，UDEC、3DEC和PFC软件为岩土及类岩石材料的力学行为基础理论研究(破裂机制与演化规律、颗粒类材料动力响应等)和工程应用研究(地下灾变机制、堆石料特性、矿山崩落开采、边坡岩土解体、爆破冲击等)提供了有效手段，应用示例如图1.1.1所示。
图1.1.1 离散单元法的工程应用
1.2 离散元颗粒流理论及PFC软件简述
颗粒流理论是离散单元法的一个重要分支。在颗粒流理论中，物体的宏观本构行为通过单元间细观接触模型实现。在具有颗粒结构特性岩土介质中的应用，就是从其细观力学特征出发，将材料的力学响应问题从物理域映射到数学域内进行数值求解。例如，物理域内的复杂实物颗粒被简化为数学域内的颗粒单元，并通过颗粒单元来构建所需几何形状的试样，颗粒间的相互作用通过接触本构关系定义，通过选择合适的本构模型及调试恰当的参数匹配材料的力学特性。
PFC软件是基于颗粒流理论的基本原理和显式差分法开发的细观力学分析软件［6］，其将介质整体离散为圆盘形(disk)或球形(sphere)颗粒单元进行分析，从细观角度探索研究对象的受力、变形、运动等力学响应。美国依泰斯卡公司于1994年首次推出颗粒流模拟软件PFC(2D/3D)1.0版本，截至目前已更新至6.0版本。其建立的计算模型由颗粒、接触及墙体构成。在二维分析时，离散颗粒为单位厚度的圆盘，在三维分析中为实心圆球。每个离散单元均为具备有限质量的刚性体，颗粒单元的直径及排列分布可根据需求设定，通过调整颗粒尺寸及粒径分布可以控制模型的孔隙率和非均匀性。墙体是面(facet)的集合，面可以组成任意复杂多变的空间多边形，在PFC2D模型中面以线段的形式表示，在PFC3D模型中则为三角形。
颗粒间的接触模型是PFC模型的核心要素，分为非黏结模型与黏结模型两类，其中非黏结模型主要用于模拟散体材料，描述其变形和运动，黏结模型在此基础上加入了强度的限制，主要用于模拟岩石及类岩石材料。对于黏结模型，当颗粒之间接触承受的应力大于其黏结强度时，黏结断裂，形成微破裂［7］。当微破裂逐渐增多时，颗粒相互运动，模型发生变形和位移，实现岩土体损伤破坏机制模拟。
截至PFC 6.0版本，内嵌的非黏结模型主要包括线性模型(linear model)、线性滚动阻滑模型(rolling resistance linear model)、赫兹模型(Hertz model)、滞回模型(hysteretic model)以及伯格斯模型(Burgers model)。非黏结模型从早期的单纯线性关系发展到考虑黏滞、滚动阻力等因素，再到考虑颗粒间的范德瓦耳斯力，模型种类越来越丰富，应用领域也更为广泛。
内嵌的黏结模型主要包括线性接触黏结模型(linear contact bond model)、线性平行黏结模型(linear parallel bond model)、光滑节理模型(smooth joint model)、平节理模型(flat joint model)、黏性线性滚动阻滑模型(adhesive rolling resistance linear model)以及软化黏结模型(soft bond model)，可统称为黏结颗粒体模型(bonded-particle model，BPM)。早期开发的BPM(线性接触黏结模型和线性平行黏结模型)主要用于模拟具有黏结特性的材料，随着对脆性岩石本质特征模拟需求的增加，依泰斯卡公司对原有线性平行黏结模型的本构关系进行改进，包括：引入平行黏结因子使平行黏结和接触黏结承载应力有先后之分；引入力矩贡献因子减小力矩对应力的贡献；引入黏结安装间距以提高颗粒配位数；使用含张拉截断的莫尔-库仑强度准则将黏结剪切强度与接触正应力相关联等。
除了对黏结本构关系的改进外，在改变单元形状方面也开展了许多工作，无论是“丛(clump)”、“簇(cluster)”还是等效晶质模型(grain-based model，GBM)，其本质都是将多个颗粒集合为一个几何形状更加多样化的单元，从而解决使用球形颗粒模拟时存在的固有不足。上述改进可增强单元间的自锁效应，但通常需要更细小的颗粒直径，从而导致计算效率相对较低。
上述黏结模型中，有两种模型被命名为节理模型，分别是光滑节理模型和平节理模型。这两种节理模型既能表征非黏结状态，也可表征黏结状态。平节理模型可用于模拟含微裂纹的硬脆性岩石，光滑节理模型通常用以表征岩体中的节理、层理以及预制裂纹等，在岩石力学及岩体工程的结构特征模拟中发挥了重要作用。
1.3 PFC软件在岩土工程领域的应用
工程应用是岩土力学学科发展的根本目的，通过岩土体的基本力学性质与基础理论研究，掌握岩土体力学性状的基本规律，进而开展工程应用研究，达到指导、优化工程实践的目的。
岩土体经受长期的地质构造作用，受结构面和节理面等弱面切割，常呈现明显的非连续性特点。PFC软件采用颗粒构建计算模型，考虑到岩土体结构的非均质、非连续等复杂特性，颗粒间的黏结会受外力作用产生微裂纹并产生不同类型的破坏，从而实现对模型内部破裂孕育和演化过程的模拟，适用于岩土体破裂机制、裂纹孕育演化规律和工程稳定性的研究。近年来，随着计算机运算性能的大幅提升，基于颗粒流理论的PFC软件已在各类岩土力学基础理论及工程方面得到广泛应用。
1.3.1 岩石力学基础理论研究
岩体可视为由岩块和结构面网络构成的复杂系统，通常表现出较强的非连续性，国内外学者针对岩体变形、强度特性及破坏机理开展了广泛研究。
岩体变形及其裂纹扩展是岩石力学的重要研究方向之一。吴顺川等［8］开展了卸载岩爆数值试验研究，得出不同应力状态下的岩样细观损伤特征；Zhang等［9，10］使用平行黏结模型研究了含预制裂隙类岩石材料的裂纹扩展现象；Yang等［11］研究了含两个不平行裂隙红砂岩单轴压缩过程中的破裂行为，揭示了其损伤破坏机理；丛宇等［12］以大理岩为例，定量研究了岩石类材料宏细观参数间的关系；Duan等［13］开展了不同受压状态下花岗岩颗粒流数值模拟试验，研究了裂隙发展过程及岩石破坏机理；Cao等［14］模拟了含多裂隙脆性岩石材料的峰值强度和破坏模式。
一般来说，脆性岩体抗拉强度远小于抗压强度，工程岩体受拉应力破坏现象显著，因此开展岩体抗拉强度研究具有重要意义。Cai等［15］通过有限元/离散元耦合的方法对巴西试验的破裂进程进行模拟，研究了岩石各向异性、预裂纹长度和方向对张拉裂纹、抗拉强度的影响；Xu等［16］研究了细观结构和细观参数对岩石巴西抗拉强度的影响；Wu等［17］研究了平台巴西圆盘的破坏过程，讨论了荷载条件和圆盘几何参数对试验结果的影响；Zhang等［18］提出采用柔性颗粒边界直接生成巴西圆盘，能有效降低巴西圆盘试样的各向异性；Ma等［19］通过数值试验研究了完整的巴西试验岩石破坏行为以及材料性质对抗拉强度的影响。
同时，节理作为一种重要的地质结构面，随机分布于岩体中，影响岩体的力学性质，诸多研究者使用PFC软件研究了节理的剪切性质。Morgan等［20，21］研究了局部细观物理力学参数、颗粒粒径分布特征、粒间摩擦强度等参量对颗粒材料剪切带形成、发展过程的影响；周喻等［22］从宏细观角度探讨了节理在直剪试验过程中的力学演化特征和破坏机制；夏才初等［23］在PFC2D模型中生成粗糙节理面，并通过模拟直剪试验研究其剪切性质；Bahaaddini等［24，25
展开
