连续-非连续性岩体的数值方法

1、连续连续-非连续性岩体的数值方法非连续性岩体的数值方法沈振中水利水电工程学院水利水电工程学院主要内容主要内容 概述（岩体数值方法）概述（岩体数值方法） 非连续岩体数值方法非连续岩体数值方法 模拟裂缝扩展模拟裂缝扩展 算例和应用算例和应用一、岩体的数值方法一、岩体的数值方法适用于连续介质适用于连续介质 有限差分法（有限差分法（FDM ） 有限单元法（有限单元法（FEM ） 边界单元法（边界单元法（BEM ） 半解析元法半解析元法 无限元法（无限元法（IFEM ） 数值流形方法（数值流形方法（NMM） 无单元法（无单元法（EFM） 适用于非连续介质适用于非连续介质刚体极限平衡法（刚体极限平衡法（L

2、EM）刚体界面元法（刚体界面元法（RBSM）关键块理论（关键块理论（KBT）离散元法（离散元法（DEM）离散有限单元法（离散有限单元法（DFEM）非连续变形分析（非连续变形分析（DDA）数值流形方法（数值流形方法（NMM）无单元法（无单元法（EFM）Discretization concepts -1FaultsJointsFig. 1-1 The fractured rock mass Discretization concepts -2Fig. 1-2 The continuum-based approachFDM or FEMJoint elementDiscretization con

3、cepts -3Fig. 1-3 The continuum-based approachBEMRegion1Region3Region4Region2Discretization concepts -4Fig. 1-4 The discontinuum-based approachDEMBlocksTable 1-1 Numerical methods in modeling jointed rock massInvestigatorYearDimensionsType of constraintsCundall19712DEM (explicit)Burman19712explicitCh

4、appel19722explicitCundall19762explicitHocking, Kawai et al.19772explicitLotstedt19792implicitCundall19802explicitDowding et al.19832explicitBelytschko19842explicitHocking et al.19853explicitShi and Goodman19882DDA (implicit)Cundall19882explicitGhaboussi et al.19902explicitLin19903explicitShi19922NMM

5、 (implicit)Wang & Garga19932explicitShi20013DDA (implicit)图图1-5 分析方法与分析方法与Q值的关系（值的关系（Barton，1995）Q为岩体分类的一个指标：为岩体分类的一个指标：当当Q100，基于连续介质的计算方法（如：，基于连续介质的计算方法（如：FEM）当当Q=0.1100，基于连非续介质的计算方法（如：，基于连非续介质的计算方法（如：DEM）樱井建议：樱井建议：岩体工程数值分析岩体工程数值分析方法的选用，不仅方法的选用，不仅要根据岩体内部的要根据岩体内部的非连续性特点，还非连续性特点，还要根据结构的尺寸要根据结构的尺寸效应。效

6、应。二、非连续岩体的数值方法二、非连续岩体的数值方法 刚体极限平衡法刚体极限平衡法(LEM) 刚体界面元法刚体界面元法(RBSM) 关键块理论关键块理论(KBT) 离散单元法离散单元法(DEM) 离散有限单元法离散有限单元法(DFEM) 非连续变形分析非连续变形分析(DDA) 数值流形方法数值流形方法(NMM) 无单元法无单元法(EFM)图图2-1 非连续块体位移示意图非连续块体位移示意图1、刚体极限平衡法、刚体极限平衡法(LEM)Limit Equilibrium Method是对岩体的简化系统进行极限平衡分析。理论是对岩体的简化系统进行极限平衡分析。理论简单，概念清晰。简单，概念清晰。图图

7、2-1-1 刚体极限平衡法示意图（刚体极限平衡法示意图（Hoek & Bray, 1977）(Limit Equilibrium Method)刚体极限平衡法刚体极限平衡法(LEM)的改进的改进Sarma提出了对滑坡体进行斜分条的改进极限平衡法，沿提出了对滑坡体进行斜分条的改进极限平衡法，沿滑坡体进行斜分条，以模拟断层节理等不连续面，且假定滑坡体进行斜分条，以模拟断层节理等不连续面，且假定条块侧面也达到了极限平衡，这样通过静力平衡条件即可条块侧面也达到了极限平衡，这样通过静力平衡条件即可唯一地确定边坡的安全系数与超载系数。唯一地确定边坡的安全系数与超载系数。Donald和陈祖煜和陈祖煜将将Sa

8、rma的静力平衡方程转化为微分方程，并通过求解该的静力平衡方程转化为微分方程，并通过求解该微分方程的闭合解得到边坡的安全系数。微分方程的闭合解得到边坡的安全系数。 2、刚体界面元法（、刚体界面元法（RBSM） 由由Kawai于于1976年提出，最初被用于梁板结构。年提出，最初被用于梁板结构。Kawai认为刚架结构的抗弯刚度远小于其法向刚度，因而假定其轴认为刚架结构的抗弯刚度远小于其法向刚度，因而假定其轴向刚度为无穷大，轴向变形可以忽略不计。在两个刚性杆中向刚度为无穷大，轴向变形可以忽略不计。在两个刚性杆中用转向弹簧连接作为模拟梁的基本单元，在两个刚性三角板用转向弹簧连接作为模拟梁的基本单元，在

9、两个刚性三角板中间用转动弹簧连接作为模拟板的基本单元，从而较多地减中间用转动弹簧连接作为模拟板的基本单元，从而较多地减少了自由度。少了自由度。 Kawai于于1977年推广到平面应变问题中，根据节理划分年推广到平面应变问题中，根据节理划分单元，认为每个单元是刚性的，相邻单元通过接触面中点的单元，认为每个单元是刚性的，相邻单元通过接触面中点的法向、切向和转动弹簧相连，从而用来分析线性或几何非线法向、切向和转动弹簧相连，从而用来分析线性或几何非线性问题。比较接近刚体界面元法了性问题。比较接近刚体界面元法了(Rigid Body Spring Method)。 Belytschko 于于1986年提

10、出了分析节理岩体的刚体界面年提出了分析节理岩体的刚体界面元法，该法要求块体之间接触关系保持不变，且为边边接触元法，该法要求块体之间接触关系保持不变，且为边边接触模式，这种模型只局限于研究节理岩体的小变形问题。模式，这种模型只局限于研究节理岩体的小变形问题。图图2-2-1 具有具有2个三角形单元的刚体界面元示意图个三角形单元的刚体界面元示意图(Kawai,1977,1991)（Rigid Body Spring Method）图图2-2-2 变形后两个单元的相对位置（变形后两个单元的相对位置（RBSM）3、关键块理论、关键块理论(KBT)在洞室、露天边坡及基础等岩体工程中，常发生某些被节在洞室、

11、露天边坡及基础等岩体工程中，常发生某些被节理裂隙完全分割的块体滑落，这些本身在几何形状上具备理裂隙完全分割的块体滑落，这些本身在几何形状上具备滑动可能的块体称为滑动可能的块体称为关键块关键块。关键块理论（关键块理论（Key Block Theory）由石根华和）由石根华和Goodman提出并发展，它适用于切割体为凸体的情况，只能分析滑提出并发展，它适用于切割体为凸体的情况，只能分析滑动而不能考虑倾倒转动，分析力的作用时采用静力平衡公动而不能考虑倾倒转动，分析力的作用时采用静力平衡公式。关键块理论中关键块的寻找是其关键，整体系统的稳式。关键块理论中关键块的寻找是其关键，整体系统的稳定以分析关键块

12、的稳定为判别标准。定以分析关键块的稳定为判别标准。4、离散单元法（、离散单元法（DEM）离散元法离散元法(Distinct Element Method)是是Cundall于于1971年年提出来的。它也像有限元那样，将区域划分成若干单元，提出来的。它也像有限元那样，将区域划分成若干单元，单元之间可以看成是角角接触、角边接触或边边接单元之间可以看成是角角接触、角边接触或边边接触，由于单元受节理等不连续面控制，在以后的运动过程触，由于单元受节理等不连续面控制，在以后的运动过程中，单元节点可以分离，且随着单元的平移和转动，允许中，单元节点可以分离，且随着单元的平移和转动，允许调整各个单元之间的接触关

13、系。单元之间相互作用的力由调整各个单元之间的接触关系。单元之间相互作用的力由力和位移的关系求出，单元的运动则完全根据单元所受的力和位移的关系求出，单元的运动则完全根据单元所受的不平衡力和不平衡力矩的大小按牛顿定律确定。块体单元不平衡力和不平衡力矩的大小按牛顿定律确定。块体单元可能达到平衡状态，也可能一直运动下去。可能达到平衡状态，也可能一直运动下去。图图2-4-1 离散元法未知量示意图（离散元法未知量示意图（Cundall, 1971）(Distinct Element Method)图2-4-2 DEM计算过程简图力力位移法则：位移法则：假定所有变形发生在角假定所有变形发生在角边接触上，力是

14、由变形产生的。边接触上，力是由变形产生的。力和位移采用增量，即一个位移变化将产生一个力，且作力和位移采用增量，即一个位移变化将产生一个力，且作用在接触已存在的力上。在一个时步内，对一个给定的接用在接触已存在的力上。在一个时步内，对一个给定的接触，切向和法向位移增量由（触，切向和法向位移增量由（u，v，r）增量计算，而新的）增量计算，而新的切向和法向力根据原作用力通过力切向和法向力根据原作用力通过力位移关系求得。然后位移关系求得。然后这些接触力被重新分解为等效的这些接触力被重新分解为等效的X、Y分力和转动分力和转动R，并与，并与其他其他X、Y分力和转动叠加作用到每个块体上。分力和转动叠加作用到每

15、个块体上。运动法则：运动法则：某一块体上，力和运动的总和被用于计算在某一块体上，力和运动的总和被用于计算在X、Y和和R 方向的加速度，方向的加速度，通过积分求得速度和位移，有了这个新的位移，就可进行下一步计算。通过积分求得速度和位移，有了这个新的位移，就可进行下一步计算。在时间域内，在时间域内，DEM采用显式数值方法，块体间的相互作用由接触和块采用显式数值方法，块体间的相互作用由接触和块体的运动检查算出。因此，体的运动检查算出。因此，DEM可以跟踪每个块体每一步的行为。可以跟踪每个块体每一步的行为。图图2-4-3 DEM块体接触时的力学模型块体接触时的力学模型DEM采用显式方法求解方程，优点是

16、在解运动方程时，由采用显式方法求解方程，优点是在解运动方程时，由于假设作用力右端项是已知的，不需要解联立方程。但是于假设作用力右端项是已知的，不需要解联立方程。但是法向接触弹簧的长度是块体位移的函数，而在接触计算中法向接触弹簧的长度是块体位移的函数，而在接触计算中位移是未知的，因此，在接触发生时，假设是不满足的。位移是未知的，因此，在接触发生时，假设是不满足的。第第n-1步法向接触弹簧的长度产生的法向接触力是第步法向接触弹簧的长度产生的法向接触力是第n步假步假设的力。这种逼近不能产生正确的接触力，有时会破坏能设的力。这种逼近不能产生正确的接触力，有时会破坏能量守恒而使得解答发散。于是，需要引入

17、一个阻尼，使得量守恒而使得解答发散。于是，需要引入一个阻尼，使得系统稳定。这样，系统稳定。这样，DEM的接触模型可以表示如下。的接触模型可以表示如下。离散单元法的离散单元法的动态松弛法动态松弛法和和静态松弛法静态松弛法动态松弛法动态松弛法是把非线性静力学问题化为动力学问题来求解，是把非线性静力学问题化为动力学问题来求解，用显式中心差分法来近似积分运动方程，并用适当的阻尼来用显式中心差分法来近似积分运动方程，并用适当的阻尼来吸收系统的动能，使系统的振动尽可能快地消失，同时场函吸收系统的动能，使系统的振动尽可能快地消失，同时场函数收敛于静态值。它按时步在计算机上迭代求解，整个计算数收敛于静态值。它

18、按时步在计算机上迭代求解，整个计算过程只需要直接代换，即利用前一步迭代的函数值近似求解过程只需要直接代换，即利用前一步迭代的函数值近似求解新函数值，对非线性问题也可以考虑。新函数值，对非线性问题也可以考虑。静态松弛法静态松弛法是直接寻找块体失去平衡后再次达到平衡时的是直接寻找块体失去平衡后再次达到平衡时的力力位移关系位移关系，采用隐式法联立平衡方程组，并以完全消除块，采用隐式法联立平衡方程组，并以完全消除块体的残余力和力矩为目标进行迭代求解。静态松弛法避免了体的残余力和力矩为目标进行迭代求解。静态松弛法避免了动态松弛法中的难点，即粘性阻尼的确定及计算时步的选取，动态松弛法中的难点，即粘性阻尼的

19、确定及计算时步的选取，但它在求解联立平衡方程组时，有时会碰到数值奇异或病态但它在求解联立平衡方程组时，有时会碰到数值奇异或病态问题，有待进一步改进。目前，工程中广泛使用的离散单元问题，有待进一步改进。目前，工程中广泛使用的离散单元法多采用动态松弛法。法多采用动态松弛法。5、离散有限单元法、离散有限单元法(DFEM) 有限元法在模拟岩体的非连续特性时存在一些缺限，如块体的接触、有限元法在模拟岩体的非连续特性时存在一些缺限，如块体的接触、相对位移、分离等。近相对位移、分离等。近10多年来，一些学者对多年来，一些学者对FEM作了一些修正，形作了一些修正，形成离散优先单元法（成离散优先单元法（Disc

20、rete Finite Element Method）主要有两类：）主要有两类： 采用特殊单元模拟非连续特性。采用特殊单元模拟非连续特性。 采用接触判断模拟非连续特性而不增加附加单元。采用接触判断模拟非连续特性而不增加附加单元。(1) Goodman(1976)最早提出最早提出 “joint element”，用于模拟非连续缝面，用于模拟非连续缝面，就是众所周知的就是众所周知的“Goodman element”。Nakazawa(1979) 利用利用FEM讨论不同摩擦条件下的接触分析，用于研究非线性接触问题。讨论不同摩擦条件下的接触分析，用于研究非线性接触问题。Lei（2001）介绍了一种简单

21、的界面单元用于接触分析，他将节点位移作）介绍了一种简单的界面单元用于接触分析，他将节点位移作为未知量，模拟两个具有初始接触或间隙的块体的接触滑动、解耦和再为未知量，模拟两个具有初始接触或间隙的块体的接触滑动、解耦和再结合。结合。 缺点：根本点都是基于小位移和小变形，难以很好地解决大位移问缺点：根本点都是基于小位移和小变形，难以很好地解决大位移问题。题。(2)(2)为模拟非连续体的大位移特性，出现了利用法向和切向弹簧来进行接触为模拟非连续体的大位移特性，出现了利用法向和切向弹簧来进行接触计算，但是块体之间不设特殊单元，有两类计算，但是块体之间不设特殊单元，有两类：(a)接触预确定接触预确定FEM

22、和和(b)接接触搜索算法触搜索算法FEM。(a) Chaudhary（1986）、）、Laursen and Simo（1993）给出了一种）给出了一种FEM 方式来处理多块体、大变形摩擦接触问题，它采用连续理论，应用通常的方式来处理多块体、大变形摩擦接触问题，它采用连续理论，应用通常的FEM离散方法，继承了连续离散方法，继承了连续FEM 的优点。的优点。Solberg等人（等人（1998）利用固体）利用固体力学中无摩擦接触问题的数值解答，将微分力学中无摩擦接触问题的数值解答，将微分代数方程的分析思想应用于代数方程的分析思想应用于FEM。这种方法采用。这种方法采用预确定接触模式（预确定接触模式

23、（pre-determined contact）来求解来求解问题，可以模拟块体间的接触、相对滑动、分离等。问题，可以模拟块体间的接触、相对滑动、分离等。(b) 为模拟真正的非连续问题（接触是每步都在变化，且预先未知的）为模拟真正的非连续问题（接触是每步都在变化，且预先未知的） ， Belytshko and Neal（1991）提出了一种）提出了一种弹球算法弹球算法“pinball algorithm” 。该方法将块体视为弹球，并置于块体表面，认为只有在弹球交迭的地方发该方法将块体视为弹球，并置于块体表面，认为只有在弹球交迭的地方发生接触，从而把复杂形状的接触判断转化为圆形的接触判断。由于过于

24、粗生接触，从而把复杂形状的接触判断转化为圆形的接触判断。由于过于粗略，略，Zhong（1993）介绍了一种基于几何形状的接触算法来解决接触）介绍了一种基于几何形状的接触算法来解决接触/嵌入嵌入的大变形问题。的大变形问题。图2-5-1 二维“弹球”模型示意图(b) Belytshko and Neal（1991）提出了一种）提出了一种弹球算法弹球算法“pinball algorithm” 。该方法将块体视为弹球，并置于块体表面，认为只有在弹球。该方法将块体视为弹球，并置于块体表面，认为只有在弹球交迭的地方发生接触，从而把复杂形状的接触判断转化为圆形的接触判断。交迭的地方发生接触，从而把复杂形状的

25、接触判断转化为圆形的接触判断。由于过于粗略，由于过于粗略，Zhong（1993）介绍了一种基于几何形状的接触算法来解）介绍了一种基于几何形状的接触算法来解决接触决接触/嵌入的大变形问题。嵌入的大变形问题。6、非连续变形分析、非连续变形分析(DDA) 石根华和石根华和Goodman于于1989年提出了一种称为非连续变形分析年提出了一种称为非连续变形分析 (Discontinuous Deformation Analysis DDA)的数值方法，它是在的数值方法，它是在非连续体位移反分析法的基础上推广而来的一种正分析方法，它可以从非连续体位移反分析法的基础上推广而来的一种正分析方法，它可以从块体结

26、构的几何参数、力学参数、外荷载及约束情况计算出块体的位移、块体结构的几何参数、力学参数、外荷载及约束情况计算出块体的位移、变形、应力、应变以及块体间离合的情况。变形、应力、应变以及块体间离合的情况。 DDA将系统视为由常应变并且可作刚性位移的块体所组成，建立将系统视为由常应变并且可作刚性位移的块体所组成，建立起类似于有限元的平衡方程（其矩阵含义与有限元不同）起类似于有限元的平衡方程（其矩阵含义与有限元不同） ，DDA在建在建立平衡方程时考虑如下两个假定（约束）：立平衡方程时考虑如下两个假定（约束）： 块体间无拉应力作用。块体间无拉应力作用。 块体间不能有叠合。块体间不能有叠合。一旦有拉力或叠合

27、产生，即在相应点上作用一反向弹簧以使拉力或叠合一旦有拉力或叠合产生，即在相应点上作用一反向弹簧以使拉力或叠合消失。从这两个约束出发，利用最小位能原理，消失。从这两个约束出发，利用最小位能原理，DDA建立了平衡方程。建立了平衡方程。其刚度矩阵由材料特性、约束情况及变形情况所决定。位移列阵包含块其刚度矩阵由材料特性、约束情况及变形情况所决定。位移列阵包含块体的平动位移、转角及应变，力列阵由各块体所受外力集成。体的平动位移、转角及应变，力列阵由各块体所受外力集成。图图2-6-1 DDA中块体的变形中块体的变形图图2-6-3 DDA采用的块体未知量采用的块体未知量平衡方程平衡方程图图2-6-4 DDA

28、块体接触的不同类型块体接触的不同类型图图2-6-5 DDA块体接触时的力学模型块体接触时的力学模型DDA与与FEM 、DEM的比较的比较DDA与有限元法（与有限元法（FEM）有相似之处，但它属于离散计算的）有相似之处，但它属于离散计算的范畴，它的单元可以任意凹凸，也可以分离，接触面上采用范畴，它的单元可以任意凹凸，也可以分离，接触面上采用摩尔摩尔库仑准则，对每级荷载以时步求解平衡方程。它的未库仑准则，对每级荷载以时步求解平衡方程。它的未知量是块体位移，而有限元的未知量是结点位移。知量是块体位移，而有限元的未知量是结点位移。DDA与离散元法（与离散元法（DEM）不同，）不同，DDA是位移解法，而

29、离散是位移解法，而离散元法通过调整接触力来使系统达到平衡，属力法范畴。此外，元法通过调整接触力来使系统达到平衡，属力法范畴。此外，DDA是一种隐式解法，而离散元法是显式解法。是一种隐式解法，而离散元法是显式解法。DDA既可用于静力问题，也可用于动力问题，已在一些工程既可用于静力问题，也可用于动力问题，已在一些工程中得到应用。目前二维程序比较成熟，但模型比较简单，应中得到应用。目前二维程序比较成熟，但模型比较简单，应用还不十分广泛。三维用还不十分广泛。三维DDA还仅仅处于开发阶段。还仅仅处于开发阶段。DDA只是只是近几年才传入我国，有待于进一步发展，是很有潜力的一种近几年才传入我国，有待于进一步

30、发展，是很有潜力的一种算法。算法。7、数值流形方法（、数值流形方法（NMM）1995年，石根华提出了年，石根华提出了DDA与与FEM的统一形式的统一形式：数值流形数值流形方法方法(NMM) 。NMM以数学流形为基础，计算结构体的位移和变形。以数学流形为基础，计算结构体的位移和变形。NMM中的网格就是中的网格就是数学覆盖数学覆盖，这些数学覆盖互相重叠并且覆盖，这些数学覆盖互相重叠并且覆盖整个计算区域，在每个数学覆盖上定义互相独立的位移近整个计算区域，在每个数学覆盖上定义互相独立的位移近似函数。这些数学覆盖被物理边界切割而形成似函数。这些数学覆盖被物理边界切割而形成物理覆盖物理覆盖，物理覆盖的重叠

31、区域形成单元。将这些覆盖上的位移函数物理覆盖的重叠区域形成单元。将这些覆盖上的位移函数结合起来形成计算域上的全域位移近似函数，在每个单元结合起来形成计算域上的全域位移近似函数，在每个单元上的近似函数就是形成此单元的若干个互相重叠的覆盖上上的近似函数就是形成此单元的若干个互相重叠的覆盖上的近似函数的加权平均，并利用最小位能原理形成整体平的近似函数的加权平均，并利用最小位能原理形成整体平衡方程。在衡方程。在NMM中，积分方法采用了单纯形上的解析积分中，积分方法采用了单纯形上的解析积分形式，这是与传统数值方法的一个很大的不同。形式，这是与传统数值方法的一个很大的不同。图2-7-1 NMM 的求解过程

32、NMM基本概念四边形覆盖图2-7-2 NMM 例子示意图三角形数学覆盖三角形数学覆盖流形单元含有裂隙的边坡覆盖系统含有裂隙的边坡覆盖系统（1）含裂隙的边坡）含裂隙的边坡（2）覆盖系统）覆盖系统（3）流形单元）流形单元 NMM 中全局函数的形成NMM与与DDA、FEM的关系的关系NMM使用数学覆盖系统，使得连续体、非连续体的整体平使用数学覆盖系统，使得连续体、非连续体的整体平衡方程都可以用统一的形式来表达。有限元法衡方程都可以用统一的形式来表达。有限元法FEM和非连和非连续变形分析续变形分析DDA都是都是NMM的特殊形式，可以说的特殊形式，可以说NMM使使DDA和有限元法和有限元法FEM在理论基

33、础上和表达形式上得到了统在理论基础上和表达形式上得到了统一，这是一个很大的贡献。一，这是一个很大的贡献。目前，目前，NMM中的覆盖只能采用一些特殊的形状，如正三角中的覆盖只能采用一些特殊的形状，如正三角形、正六边形等，而实现任意形状的覆盖则有一定的困难。形、正六边形等，而实现任意形状的覆盖则有一定的困难。因为虽然在理论上可以采用任意形状的覆盖，但真正操作因为虽然在理论上可以采用任意形状的覆盖，但真正操作起来，对任意形状的覆盖，权函数的选取及解析积分的求起来，对任意形状的覆盖，权函数的选取及解析积分的求解都有一定的困难。在这些方面，还有很多工作要做。所解都有一定的困难。在这些方面，还有很多工作要

34、做。所以，目前工程界使用较多的仍然是以，目前工程界使用较多的仍然是DDA。NMMDDANMMFEMNMMFEM，三角形单元NMMFEM，四边形单元，四边形单元8、无单元法（、无单元法（EFM）无单元法属于连续介质的范畴，但是，它可以比较方便地无单元法属于连续介质的范畴，但是，它可以比较方便地处理裂缝扩展跟踪问题，而且，也是近来的研究热点。因处理裂缝扩展跟踪问题，而且，也是近来的研究热点。因此，这里作一简单介绍。此，这里作一简单介绍。无单元法的思想最早由无单元法的思想最早由Nayroles等人于等人于1992年提出，他们年提出，他们称之为虚拟单元法（称之为虚拟单元法（Diffuse Elemen

35、t Method）。）。Nayroles的近似比较简单，使用了低阶积分，边界条件的的近似比较简单，使用了低阶积分，边界条件的引入也不准确，而且在拟合函数的求导过程中忽略了一项，引入也不准确，而且在拟合函数的求导过程中忽略了一项，计算比较粗糙。但是，他们第一次在偏微分方程的数值解计算比较粗糙。但是，他们第一次在偏微分方程的数值解法中引入了法中引入了滑动最小二乘法滑动最小二乘法的思想。事实上，滑动最小二的思想。事实上，滑动最小二乘法很早就有人提出来了，如乘法很早就有人提出来了，如Lancaster等，只是在等，只是在Nayroles之前，它只被应用于曲线、曲面拟合领域。之前，它只被应用于曲线、曲面

36、拟合领域。滑动最小二乘法原理滑动最小二乘法原理正规方程（法方程）无单元伽辽金法（无单元伽辽金法（EFGM ）Belytschko等人于等人于1994年对年对Nayroles的方法进行了改进，提出了无的方法进行了改进，提出了无单元伽辽金法（单元伽辽金法（Element-free Galerkin Methods），主要作了如），主要作了如下改进：下改进： 将计算域上的积分分离为很多子域上的积分，在每个子域上采用高将计算域上的积分分离为很多子域上的积分，在每个子域上采用高阶高斯积分。阶高斯积分。 在拟合函数的导数中加入被在拟合函数的导数中加入被Nayroles忽略了的项。忽略了的项。 用拉格朗日乘

37、子法引入边界条件。用拉格朗日乘子法引入边界条件。在此基础上，在此基础上，Belytschko从变分原理出发，建立了类似于有限元的无从变分原理出发，建立了类似于有限元的无单元法基本方程。单元法基本方程。Belytschko方法的一些不足：方法的一些不足： 首先，在求解滑动最小二乘拟合函数时，需要求解线性方程组，因首先，在求解滑动最小二乘拟合函数时，需要求解线性方程组，因此，在集成刚度矩阵时，在每个高斯点上都要求解一个线性方程组，此，在集成刚度矩阵时，在每个高斯点上都要求解一个线性方程组，在求位移、应力、应变时也是如此，这是非常繁琐的，而且容易造在求位移、应力、应变时也是如此，这是非常繁琐的，而且

38、容易造成大的误差。成大的误差。 其次，用拉格朗日乘子引入边界条件后，增加了未知量，也破坏了其次，用拉格朗日乘子引入边界条件后，增加了未知量，也破坏了刚度矩阵的正定性及带形分布。增加了方程求解的时间和难度。刚度矩阵的正定性及带形分布。增加了方程求解的时间和难度。数值方法的精度比较数值方法的精度比较xy50 x1x2x3x4y3c.y1y2y4xx50 x1x2x3x4y3d.xx50 x1x2x3x4b.a. 具有一个非连续点的一维函数，细线是函数真解，粗线是数值逼近(a) FDM逼近(b) FEM逼近，一维单元：x0 x1, x1x2, x2x3, x3x4, x4x5 (c) DEM/DDA

39、逼近，一维块体：y0 x1, y1x2, y2x3, y3x4, y4x5 (d) NMM逼近，一维物理覆盖：U1 = x0 x1, U2 = x0 x2, U3 = x1x3,U4 = x2x3, U5 = y3x4, U6 = y3x5, U7 = x4x5 由于x3y3附近的突变，数学覆盖x2x4被分成2个物理覆盖（U4和U5）三、模拟裂缝扩展三、模拟裂缝扩展 DEM、DDA与模拟裂缝扩展与模拟裂缝扩展 NMM与模拟裂缝扩展与模拟裂缝扩展 EFM与模拟裂缝扩展与模拟裂缝扩展DEM、DDA与模拟裂缝扩展与模拟裂缝扩展非连续介质的数值方法，如离散元法非连续介质的数值方法，如离散元法DEM、非

40、连续变形分析、非连续变形分析DDA，都，都可以用来模拟岩体的开裂，甚至可以模拟结构大位移、大转动。但是，可以用来模拟岩体的开裂，甚至可以模拟结构大位移、大转动。但是，由于受模型本身的限制，它们都有其局限性：由于受模型本身的限制，它们都有其局限性：1. 块体是由岩体节理裂隙充分切割而形成的，块体间只简单地满足摩尔块体是由岩体节理裂隙充分切割而形成的，块体间只简单地满足摩尔库仑准则，对于结构面已将岩体完全切割成块体状的情况，采用这一库仑准则，对于结构面已将岩体完全切割成块体状的情况，采用这一理论是合适的，如用来处理节理裂隙发育地区的地下结构工程的塌方与理论是合适的，如用来处理节理裂隙发育地区的地下

41、结构工程的塌方与支护问题，是相当实用和有效的。但当岩体并未被裂隙切割成块体集合支护问题，是相当实用和有效的。但当岩体并未被裂隙切割成块体集合时或者当裂隙不发育时，这些方法就无优势可言了。时或者当裂隙不发育时，这些方法就无优势可言了。2. 模拟开裂时，受块体形状所限，开裂面只能沿块体边界，而由于裂隙模拟开裂时，受块体形状所限，开裂面只能沿块体边界，而由于裂隙扩展方向事先是无法知道的，因此，用这种方法计算岩体开裂，将随着扩展方向事先是无法知道的，因此，用这种方法计算岩体开裂，将随着不同的块体结构得到不同的破坏形式，很难得到符合实际的结果。不同的块体结构得到不同的破坏形式，很难得到符合实际的结果。可

42、见，可见，DEM、DDA在模拟切割完整的结构面的破坏时是非常有效的，在模拟切割完整的结构面的破坏时是非常有效的，但若要裂缝扩展则有一定的困难。但若要裂缝扩展则有一定的困难。 NMM与模拟裂缝扩展与模拟裂缝扩展NMM方法使连续介质与非连续介质的计算从理论上得到了方法使连续介质与非连续介质的计算从理论上得到了统一，这是一个很大的贡献。统一，这是一个很大的贡献。NMM在计算连续体与非连续在计算连续体与非连续体的大变形或进行动力分析时非常有效。在开裂跟踪方面，体的大变形或进行动力分析时非常有效。在开裂跟踪方面，NMM也得到了应用。但是，由于也得到了应用。但是，由于NMM的近似函数在数学覆的近似函数在数

43、学覆盖上定义，由物理边界切割数学覆盖而成的物理覆盖互相盖上定义，由物理边界切割数学覆盖而成的物理覆盖互相重叠形成单元，单元上的近似函数为各物理覆盖近似函数重叠形成单元，单元上的近似函数为各物理覆盖近似函数的加权平均，因此，虽然的加权平均，因此，虽然NMM由于采用了可以不变的数学由于采用了可以不变的数学覆盖而比传统的有限元法更为灵活，但它实际上仍未摆脱覆盖而比传统的有限元法更为灵活，但它实际上仍未摆脱单元的限制，在物理边界发生变化后，单元也应作相应的单元的限制，在物理边界发生变化后，单元也应作相应的调整。在开裂跟踪时，调整。在开裂跟踪时，NMM有以下局限。有以下局限。 NMM在模拟开裂时的局限在

44、模拟开裂时的局限1. 为了提高裂尖的计算精度，通常需要在裂尖加密网格。为了提高裂尖的计算精度，通常需要在裂尖加密网格。NMM在这种情况下就需要加密裂尖的数学覆盖，如果数学在这种情况下就需要加密裂尖的数学覆盖，如果数学覆盖不变，则达不到要求的计算精度，但加密数学覆盖与覆盖不变，则达不到要求的计算精度，但加密数学覆盖与有限元的重新划分网格是同样繁琐的。有限元的重新划分网格是同样繁琐的。2. 在数学覆盖不变的情况下，在数学覆盖不变的情况下，NMM的开裂长度应贯穿某个的开裂长度应贯穿某个数学覆盖或物理覆盖。若非如此，开裂后数学覆盖或物理覆盖。若非如此，开裂后NMM的物理覆盖的物理覆盖仍不会改变，则单元

45、也不会改变，其结果就是新开裂的裂仍不会改变，则单元也不会改变，其结果就是新开裂的裂缝得不到模拟。因此，目前用缝得不到模拟。因此，目前用NMM计算开裂时，多让裂缝计算开裂时，多让裂缝沿单元边界贯穿单元，以简化单元的重新调整过程。这在沿单元边界贯穿单元，以简化单元的重新调整过程。这在模拟受力简单的裂纹时是可以的，但对于岩土工程中复杂模拟受力简单的裂纹时是可以的，但对于岩土工程中复杂受力条件下的裂缝，由于事先无法确定开裂方向，就会产受力条件下的裂缝，由于事先无法确定开裂方向，就会产生较大的误差，难以得到符合实际的结果。生较大的误差，难以得到符合实际的结果。EFM与模拟裂缝扩展与模拟裂缝扩展无单元法（

46、无单元法（EFM）采用滑动最小二乘法拟合场函数，因此，）采用滑动最小二乘法拟合场函数，因此，在计算中只需计算域边界条件和结点，不需要任何单元信在计算中只需计算域边界条件和结点，不需要任何单元信息；另外，由于滑动最小二乘法拟合函数具有高阶连续的息；另外，由于滑动最小二乘法拟合函数具有高阶连续的特点，因此，无单元法具有前后处理简单、计算精度高的特点，因此，无单元法具有前后处理简单、计算精度高的优点。这使得无单元法特别适用于岩土工程中的稳定、开优点。这使得无单元法特别适用于岩土工程中的稳定、开裂分析，特别是在开裂问题中，无单元法可以很好地模拟裂分析，特别是在开裂问题中，无单元法可以很好地模拟裂缝尖端

47、奇异场，算出高精度的应力强度因子，较为方便裂缝尖端奇异场，算出高精度的应力强度因子，较为方便地跟踪裂缝扩展，而不存在重新划分网格的问题。地跟踪裂缝扩展，而不存在重新划分网格的问题。 四、算例和应用四、算例和应用 非连续变形分析（非连续变形分析（DDA） 数值流形方法（数值流形方法（NMM） 无单元法（无单元法（EFM）1、DDA：边坡稳定分析：边坡稳定分析El.702.80El.676.80F水库水位El.650.00FF123175 mEl.630.00图图1-1(a) DDA概化模型概化模型（日本）下小鳥湛水池不連続岩盤斜面安定性評価（日本）下小鳥湛水池不連続岩盤斜面安定性評価 图1-1(b) DDA计算的边坡不稳定的情况1、DDA：边坡落石分析：边坡落石分析图1-2(a) 9号国道A017断面边坡坍塌安全分析评价的DDA模型（Step=0 ）（日本）国道号美方郡温泉歌長、養父郡関宮町尾崎斜面落石、崩壊検討图1-2(b) Step=1421、DDA：边坡稳定