1.2 非连续介质力学数值方法综述
1.极限分析法
极限平衡法是岩土工程稳定分析中最为广泛使用的一种方法。该法以摩尔—库仑抗剪强度理论为基础，对滑动岩体进行力的平衡分析，结合结构面的强度参数得到抗滑稳定安全系数。该法不涉及荷载施加的过程，直接分析结构的破坏极限状态，概念清晰，理论简单。
1996年，Duncan［14］列表总结了二十多种常用的极限分析方法，给出了每种方法采用的不同假定，并比较了各种方法的优缺点；1980年潘家铮［15］在详细分析建筑物和地基抗滑稳定分析的各种方法后，在其著作中提出了“潘家铮最大最小原理”，即：①滑坡如能沿许多滑面滑动，则失稳时，它将沿抵抗力最小的一个滑面破坏（最小值原理）；②滑坡体的滑面确定时，滑面上的反力（以及滑坡体内的内力）能自行调整，以发挥最大的抗滑能力（最大值原理）。后来该原理被陈祖煜［16］所证明，从而使极限平衡分析方法和结构分析方法一样成为严格的理论分析体系。
极限平衡方法也存在缺点［17-18］：①平衡的条件未被全部满足；不能考虑全部的力平衡和力矩平衡条件，不是严格意义上的极限平衡法。②即使只满足部分平衡条件，不同的方法所采取的假定也互不相同；这些假定纯粹是为了使问题静定而引入的，没有任何物理基础。③假定滑动方向；由于实际滑坡空间形态较为复杂，准确地确定滑动方向是很困难的。
2.刚体弹簧元法
刚体弹簧元法（Rigid Body-Spring Model，RBSM）是1976 年日本东京大学教授Kawai［19］首次提出来的一种非连续介质力学数值方法，1978年该方法由周维垣、杨若琼［20］首次引入我国。此法首先将结构体离散化为一系列块体，每个块体包含六个自由度，块体与块体之间用弹簧连接，每个块体本身均是一个刚性体。该方法本质上是由有限元法演绎而来的，与有限元不同的是，刚体弹簧元在块体形心处插值，用块体形心的位移作为基本未知量，用分片的刚体位移去逼近实际整体位移场。结构内部变形通过弹簧的变形来体现，结构内部应力通过单元交界面上的面力来体现。
该法能反映结构面的不连续特征，可以模拟开裂、错动和滑移等问题，但是其缺点是：①获得的运动方程随着块体间接触点位置的变化而变化，给最终线性方程组的求解带来困难；②块体为刚体，变形均集中到交界面的弹簧中，这种变形等效存在一定误差；③与离散元等方法刚度的取值一样，刚体弹簧元特征长度的取值存在较大随机性，对应力和位移结果影响较大；④用于求解空间块体大范围运动等强几何非线性问题时存在一定的难度。当然，随着刚体弹簧元研究的进一步深入，这些缺点也正逐步被克服。
自从该法引进我国之后，国内众多学者对其进行了深入研究。张建海等［21］研究了二维和三维刚体弹簧元方法，并将之应用到锦屏高边坡和溪洛渡的拱坝坝肩稳定分析中；钱令希、张雄［22］等给出了刚体弹簧元的数学解释，并指出刚体弹簧元给出的应力精度不会低于甚至会高于位移精度；赵宁、卓家寿［23-24］基于刚体弹簧元提出了改进的刚体界面元方法，并与有限元耦合，研究了Koyna重力坝的自振特性和动力响应；陆晓敏［25］考虑了块体本身的变形，并应用到裂隙岩体的数值模拟中；陈欣［26］研究了界面元法的开裂和损伤梯度模型，并应用到小湾拱坝和锦屏拱坝的三维破坏分析中。
3.离散单元法
离散单元法是1971年由Cundall［27］提出的一种分析节理岩体的数值计算方法，最初是为了模拟岩质边坡的破坏过程。离散元法的基本原理是：①将整个研究区域根据其自身的节理裂隙或者人为分割而分成若干个相互独立的块体，如果需要考虑变形，则在块体内部进一步细分三角形（二维）或者四面体（三维）网格，使得刚性体的运动满足牛顿第二定律；②块体与块体之间用法向弹簧和切向弹簧连接，界面破坏准则满足摩尔—库仑定律；③基于牛顿第二定律采用显式的时步步进的动态松弛算法求解块体系统的运动和变形；④块体与块体之间的接触面可以打开、滑移与翻转，应力—位移关系控制块体之间的接触变化情况。
后来Cundall［28］又提出了可变形的二维离散元模型，并开发了商用程序UDEC，在这一程序中，Cundall进一步在块体内部差分三角形网格，将针对块体的平动和转动方程统一为针对节点的平动方程，变形体的计算采用有限差分法，接触检索依然在块体的表面三角形中进行，这一程序的问世，进一步推动了离散元的发展，使得其处理范围从刚体进一步发展到变形体，因此得到了广泛的应用。
在离散元方法中，计算时步、阻尼、结构面刚度的选取往往关系着离散元计算的成败，很多学者对此进行了大量的研究。魏群［29］借助于计算机和CAD技术，直接将不连续岩体的破坏过程，应力应变场，位移场及力场的分布迅速直观地显示出来，并用离散元进行岩体抗滑稳定的初步研究，用高速摄影机和激光离散方法做了验证试验，计算结果和试验有着相当好的吻合；1991年，魏群在《散体单元法的基本原理数值方法及程序》［30］一书中对离散元的计算方法和数据结构做了很大的改进：采用“静力抗衡”和“失稳破坏”两个阶段区别处理的方法研究散体结构的渐进破坏过程，在试验验证的基础上，提出了模拟渐进破坏过程的散体单元法；从随机理论出发，提出了对任意形状边界内任意形状离散颗粒组合体的数值分析方法，其中包括一种新的“椭圆散体单元法”等，进一步完善和改进了非连续颗粒体的应力应变数值计算模型；提出了岩体锚杆支护系统的散体单元分析方法，并完成了散体单元变形及破坏过程的计算机模拟技术以及CAD 动态显示等研究，编制了大型电算程序，对实际工程的高边坡稳定进行了计算分析，另外，该程序可用于研究节理岩体或颗粒组合体变形和稳定的宏观预测，也可用于微观非连续散粒体的内部力学机理研究。
1985年Cundall［31］又提出了应用于岩石力学分析的三维刚体离散元法并于1988年在公开刊物上发表文献对其进行了详细的阐述，包括数据结构、公共面算法、运动学和动力学方程等方面。1991年王泳嘉教授［32］将三维离散元引入国内后，鲁军［33］研究了三维块体的自动剖分，并根据Cundall的思想建立了三维刚体离散元模型；崔玉柱［34］对其进行了进一步的发展，提出了三个等效原则；王泳嘉［35］基于Cundall的思想自行研制了三维刚体离散元软件TRUDEC并将其应用于采矿研究中；焦玉勇等［36］基于三维刚体离散单元法提出了地下水、锚杆的三维分析方法，给出了相应的算法和公式。张冲［37］等研究了三维变形体离散元，并进行了拱坝破坏模拟和拱坝—坝肩动静力稳定分析研究。张楚汉、金锋等［38］在岩石和混凝土离散—接触—断裂分析及高坝地基、边坡稳定分析方面做了大量的研究工作，并且在三维离散元数值仿真应用方面取得了卓有成效的工作。
在三维接触关系搜索上，Cundall［31，39］提出了著名的公共面接触算法，该算法很容易扩展到三维，并被用于3DEC软件中，但该算法的效率依赖于公共面初始位置与最终位置的交角，此外，公共面位置的确定需要经过多次迭代。2004年，Erfan［40-41］等提出了一套不需要迭代的快速公共面接触算法，指出块体间最终的公共面就是少数几个确定的参与面之一，再按一定的方向对几个平面一一检验就可以确定公共面，极大地提高了公共面算法的效率。Cundall［28］提出了二维的修圆算法，鲁军［33］等将其扩展到三维，该算法将多边形的棱角修理成了圆角，从而保证了边与边之间的连续过渡；张冲［37］在修圆算法基础上，进一步优化，建立了一套基于可见性和后修圆的块体接触算法，提高了计算效率。
4.不连续变形分析
不连续变形分析方法是继离散元法之后，石根华在20世纪80年代［42-46］后期提出的一种新型分析裂隙岩体的非连续介质力学数值方法，该方法平行于有限元，把结构面切割而成的块体单元作为基本单元，假定描述块体运动的位移形函数，选择位移作为联立方程式的未知数，利用最小势能原理建立控制方程。
DDA 理论的基本特点有：①满足三个收敛准则，分别是位移收敛准则，惯性收敛准则，开—闭迭代收敛准则；②有一套完备的块体系统接触搜索和接触开—闭模拟理论；③单纯形积分技术，是一种解析法积分，确保DDA算法收敛性和精确性。
DDA自问世以来已经获得广泛的应用。石根华［42-46］将DDA应用在岩土工程中桥梁、隧洞、边坡等变形破坏问题上；Amadei［47］研究岩质边坡块体破坏问题；Chen等［48］将DDA方法应用到岩石边坡工程的稳定性分析中；Yeung［49］利用DDA方法计算了斜面上滑块的滑动、阶梯形斜面上的方形滑块群的滑动及倾倒、三支点梁以及地下矿井顶板的破坏，并将模拟结果与理论解进行了比较；邬爱清，丁秀丽等在DDA的应用方面做了很多工作，如地下洞室围岩稳定分析［50］、岩质边坡稳定分析［51］、千将坪滑坡稳定性分析［52］、唐家山堰塞湖稳定性分析［53］等；杨军、陈鹏万［54］等将DDA应用到爆破的数值分析上；张国新等［55］对DDA进行扩展，推导了水压、扬压力等荷载矩阵以及抗滑稳定局部、整体安全系数的求解方法，在此基础上研究了重力坝抗滑稳定问题；张国新等［56］在DDA方法基础上，考虑裂隙网络中地下水的流动，以及渗流压力与岩石变形的耦合作用，研究了岩石边坡稳定问题；黄涛、张国新［57］采用三阶多项式作为块体的位移模式，能够精确地描述块体内部的位移和应力状态，实现块体初应力传递和曲边接触等问题模拟；为研究土的复杂变形机理问题，张国新、李广信等［58］采用DDA方法对土的平面应变试验进行了模拟。
在三维DDA方法的研究方面，国内外学者做了一些研究工作，主要集中在基本理论研究和接触算法模拟方面。Gen-hua Shi［59-60］发表了三维不连续变形分析基本理论，以及三维DDA方法中的节理产生、块体切割和块体搜索问题；Yeung［61］利用楔形模型、块体理论、三维DDA三种方法研究楔形体稳定问题，验证三维DDA方法有效性和正确性；Zhao［62］发展了球形颗粒单元的DDA方法并对球形颗粒材料体进行了计算；姜清辉［63］对三维DDA方法中的基本理论公式进行了深入地推导和研究；刘君［64］研究了三维DDA方法与有限单元法的耦合算法；Moosavi［65］利用理论分析结果对动态三维DDA方法进行了验证；张扬、邬爱清［66］等研究了三维高阶DDA静力问题；黄涛、张国新［67］提出了面—面接触模型，在此基础开发了三维DDA软件，并对某水电站边坡进行了稳定性分析。
5.关键块体理论
关键块体理论最早由Warburton［68-69］和石根华、Goodman［70］各自独立提出并发展的一种分析裂隙岩体稳定性的方法。它不需要进行任何的应力与应变分析，而是通过拓扑学知识来判断那些由若干裂隙面和临空面所形成的、不受几何约束的且沿某个方向具有平动和转动倾向的岩石块体，这些岩石块体被称为“关键块体”。关键块体破坏后常常要引发周围更多的块体的破坏，因此，寻找关键块体，及时对其进行支护、加固，可防止滑坡、坍塌等事故的发生。它是岩体工程中边坡及地下洞室围岩稳定分析和支护设计中非常有效的方法。这种方法非常适合硬岩工程，并不适合软岩工程。这主要是因为后者中岩体介质的应力、应变及破坏占据主要位置。
关键块体理论在国内外工程实例中得到应用，并取得了一定的成果。Stone 和 Young［71］、Mauldon［72］、Hatzor［73］、Jakubowski 和 Tajdus［74］、Kusznaul 和 Goodman［75］等研究随机性块体理论问题；Karaca 和 Goodman［76］研究了水荷载作用下关键块体的稳定问题；Yarahlnadi、Bafghi、Verdel［77］提出了“关键块体群方法”（Key Group Method），该方法将相邻的块体组合在一起，判断组合后的“大块体”是否稳定；Wibowo［78］研究了块体的二次破坏问题；裴觉民［79］及张奇华［80］研究了水电工程地下厂房关键块体稳定性问题。
6.数值流形方法
为了统一解决连续和非连续变形的力学问题，石根华于20世纪90年代，在其提出的非连续变形分析方法基础上，利用流形分析中的有限覆盖技术创立了一种新的数值计算方法——数值流形法，并率先应用在块体与节理岩体的变形模拟中［81-83］。此方法克服了以往方法的不足，既吸收了有限元等连续力学方法的优点，又吸收了DDA等不连续方法的长处。其像有限元法那样在连续介质中划分单元，使该方法与有限元同等的精度计算连续域内的应力分布，物理网格和数学网格的分离及一套完善的接触理论的引入，使该法能很好地模拟接触、物体运动等不连续力学问题，覆盖的使用又使该方法较容易和客观地模拟裂缝的产生和扩展。
流形方法以拓扑流形和微分流形为基础，它有分开的且独立的两套网格：数学覆盖和物理网格。数学覆盖定义近似解的精度，由用户选择且独立于物理网格划分，可以是规则或不规则的格子，物理网格则包括块体的物理边界、裂缝、块体和不同材料区域的交界面等，它不能人为选择。物理网格对数学覆盖再剖分就形成了覆盖材料全域的流形方法求解的物理覆盖系统。对物理覆盖系统上的每一个物理覆盖可以定义各自独立的覆盖位移函数，然后在几个覆盖的公共区域（即流形单元）内将其所有覆盖上的独立覆盖位移函数加权求和就能形成适应于该域的总体位移函数。借助于建立的物理覆盖系统，流形方法能够很方便地分析有裂缝的结构块体和明显可见的大变形和大位移；而且流形方法覆盖位移函数的采用也能够为偏微分方程的求解创造出更高的效率，提供更好的算法。
近年来，有关数值流形方法及其应用的研究已涉及很多方面，如裂纹扩展跟踪［84-86］、P型自适应分析［87］、渗流问题的求解［88-89］、非线性大变形行为的模拟［90］等，展现出较为广阔的应用前景。Shyu和Salami［91］借鉴有限元法的研究经验，在数值流形方法中引入了四节点等参单元；1997年，王水林等［92］研究了流形方法在模拟裂纹扩展中的实现过程；Chen和Ohnishi等［93］用完备的高阶覆盖函数来提高流形方法的求解精度，并给出了其具体实现过程；张国新对石根华的数值流形法做了较大扩充，采用了二阶流形法［94-95］，引入Newmark 方法分析动力学问题，并与子域奇异边界元法相结合，模拟裂缝沿任意方向扩展及结构的地震破坏问题［96］；对流形元考虑裂纹渗流和变形耦合，用断裂力学方法求解应力强度因子，采用杆单元模拟钢筋作用并成功模拟钢筋混凝土结构破坏［97］；增加了带强度节理的模拟［98］及裂纹在块体内的扩展和完整块体破坏模拟的功能［99］等；李海枫、张国新等系统地研究了有限元网格覆盖下的三维数值流形单元的生成问题［100］；2009年11月25—27日，在新加坡召开第9届非连续变形分析国际研讨会（ICADD9），新加坡南洋理工大学马国伟研究团队公布关于三维数值流形方法的最新研究成果。
7.颗粒体离散元
自1971 年Cundall 提出适用于岩石力学的块体离散元法后，1979年，Cundall 和Strack［101-103］又联合提出了适用于土力学的颗粒体离散元方法，并开发了二维圆盘程序Ball和三维圆球程序Trubal，这两个程序在Itasca公司的进一步开发下，形成了最后的商业球体离散元软件，即PFC2D［104］和PFC3D［105］。其基本思想和刚体离散元是一致的，以三维为例，两者的区别在于研究的对象不同，刚体离散元研究的是块体，而颗粒体离散元研究的是球。正是由于颗粒体离散元研究的是球，所以接触关系检索的难度大大降低，从而在同样的计算条件下，能够计算的球体数量飞速提高，因此，和块体离散元相比，球体直径可以很小，数量可以很多。也正是由于这个特点，使得颗粒体离散元方法被广泛应用于细观力学层次的研究中。需要指出的是，Cundall提出的颗粒体离散元方法中，球一般是指刚性球，但球与球之间是可以发生相互嵌入的，嵌入量的多少也决定了球与球之间的相互作用力，一般学术界将这种允许嵌入的接触称为软接触，因此Cundall提出的颗粒体离散元称为软球模型离散元。与软球模型离散元相对，1980年，Campbell［106］进一步提出了硬球模型离散元，并应用于剪切流的分析中。
自从颗粒体离散元提出后，国际上众多学者纷纷对其进行了研究。英国Aston 大学引入Cundall的Trubal程序，并改造成Granule程序，用其模拟干-湿、弹性-塑性和颗粒两相流问题［107］；Venugopal等［108］用自编的颗粒体离散元程序研究了料仓的进料和拌料过程；Qiming等［109］用PFC2D程序研究了波在散粒体中的传播过程；Deluzarche 等［110］用PFC2D程序进行了堆石坝堆石、密实和失稳研究。
我国颗粒体离散元研究始于20世纪80 年代，王泳嘉［111］引入Cundall的颗粒体离散元方法，并将之应用于煤矿出煤、边坡稳定分析中。李世海等［112］用PFC方法对边坡开挖进行了研究。张洪武［113］考虑了接触面粘连的求解技术，并将之用于细观非线性问题的分析中。周健［114］推导了颗粒流显式求解的稳定条件，并对土中渗流问题进行了模拟。
8.离散裂隙网络方法
离散裂隙网络方法是利用裂隙网络来考虑裂隙岩体中液体流动和传导的一种离散模型方法。这种方法早在1980年左右就被建立［115-125］，经过这么多年的发展，已应用到土木、环境、水库以及其他岩土工程领域。这方面研究可参见Bear［126］、Sahimi［127］、Adler和Thovert［128］等的著作，以及美国国家科学研究委员会1996年的研究报告［129］。最近的关于DFN的文献综述可参见Berkowitz的文献［130］。
目前，有很多不同的DFN生成和计算分析代码。其中，最著名的是FRACMAN/MAFIC［131］和NAPSAC［132-135］。