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

连续-非连续性岩体的数值(shùzí)方法沈振中水利水电工程学院共八十七页主要(zhǔyào)内容概述（岩体数值方法）非连续岩体数值方法模拟(mónǐ)裂缝扩展算例和应用共八十七页一、岩体的数值(shùzí)方法适用于连续介质有限(yǒuxiàn)差分法（FDM）有限单元法（FEM）边界单元法（BEM）半解析元法无限元法（IFEM）数值流形方法（NMM）无单元法（EFM）适用于非连续介质刚体极限平衡法（LEM）刚体界面元法（RBSM）关键块理论（KBT）离散元法（DEM）离散有限单元法（DFEM）非连续变形分析（DDA）数值流形方法（NMM）无单元法（EFM）共八十七页Discretizationconcepts--1FaultsJointsFig.1-1Thefracturedrockmass共八十七页Discretizationconcepts--2Fig.1-2Thecontinuum-basedapproach——FDMorFEMJointelement共八十七页Discretizationconcepts--3Fig.1-3Thecontinuum-basedapproach——BEMRegion1Region3Region4Region2共八十七页Discretizationconcepts--4Fig.1-4Thediscontinuum-basedapproach——DEMBlocks共八十七页Table1-1NumericalmethodsinmodelingjointedrockmassInvestigatorYearDimensionsTypeofconstraintsCundall19712DEM(explicit)Burman19712explicitChappel19722explicitCundall19762explicitHocking,Kawaietal.19772explicitLotstedt19792implicitCundall19802explicitDowdingetal.19832explicitBelytschko19842explicitHockingetal.19853explicitShiandGoodman19882DDA(implicit)Cundall19882explicitGhaboussietal.19902explicitLin19903explicitShi19922NMM(implicit)Wang&Garga19932explicitShi20013DDA(implicit)共八十七页图1-5分析方法与Q值的关系（Barton，1995）Q为岩体分类的一个指标(zhǐbiāo)：当Q<0.1orQ>100，基于连续介质的计算方法（如：FEM）当Q=0.1～100，基于连非续介质的计算方法（如：DEM）共八十七页樱井建议：岩体工程数值分析方法的选用，不仅(bùjǐn)要根据岩体内部的非连续性特点，还要根据结构的尺寸效应。共八十七页二、非连续(liánxù)岩体的数值方法刚体极限平衡法(LEM)刚体界面元法(RBSM)关键块理论(KBT)离散单元法(DEM)离散有限单元法(DFEM)非连续变形分析(fēnxī)(DDA)数值流形方法(NMM)无单元法(EFM)共八十七页图2-1非连续块体(kuàitǐ)位移示意图共八十七页1、刚体(gāngtǐ)极限平衡法(LEM)LimitEquilibriumMethod是对岩体的简化系统进行(jìnxíng)极限平衡分析。理论简单，概念清晰。图2-1-1刚体极限平衡法示意图（Hoek&Bray,1977）(LimitEquilibriumMethod)共八十七页刚体极限平衡(pínghéng)法(LEM)的改进Sarma提出了对滑坡体进行斜分条的改进极限平衡法，沿滑坡体进行斜分条，以模拟断层节理等不连续面，且假定条块侧面也达到了极限平衡，这样通过静力平衡条件即可唯一地确定(quèdìng)边坡的安全系数与超载系数。Donald和陈祖煜将Sarma的静力平衡方程转化为微分方程，并通过求解该微分方程的闭合解得到边坡的安全系数。共八十七页2、刚体(gāngtǐ)界面元法（RBSM）

由Kawai于1976年提出，最初被用于梁板结构。Kawai认为刚架结构的抗弯刚度远小于其法向刚度，因而假定其轴向刚度为无穷大，轴向变形可以忽略不计。在两个刚性杆中用转向弹簧连接作为模拟梁的基本单元，在两个刚性三角板中间用转动弹簧连接作为模拟板的基本单元，从而较多地减少(jiǎnshǎo)了自由度。

Kawai于1977年推广到平面应变问题中，根据节理划分单元，认为每个单元是刚性的，相邻单元通过接触面中点的法向、切向和转动弹簧相连，从而用来分析线性或几何非线性问题。比较接近刚体界面元法了(RigidBodySpringMethod)。

Belytschko于1986年提出了分析节理岩体的刚体界面元法，该法要求块体之间接触关系保持不变，且为边边接触模式，这种模型只局限于研究节理岩体的小变形问题。共八十七页图2-2-1具有2个三角形单元(dānyuán)的刚体界面元示意图(Kawai,1977,1991)（RigidBodySpringMethod）共八十七页图2-2-2变形(biànxíng)后两个单元的相对位置（RBSM）共八十七页3、关键(guānjiàn)块理论(KBT)在洞室、露天边坡及基础等岩体工程中，常发生某些被节理裂隙完全分割的块体滑落，这些本身在几何形状上具备滑动可能的块体称为关键块。关键块理论（KeyBlockTheory）由石根华和Goodman提出并发展，它适用于切割体为凸体的情况，只能分析滑动而不能考虑倾倒(qīngdǎo)转动，分析力的作用时采用静力平衡公式。关键块理论中关键块的寻找是其关键，整体系统的稳定以分析关键块的稳定为判别标准。共八十七页4、离散(lísàn)单元法（DEM）离散元法(DistinctElementMethod)是Cundall于1971年提出来的。它也像有限元那样，将区域划分成若干单元，单元之间可以看成是角－角接触、角－边接触或边－边接触，由于单元受节理等不连续面控制，在以后的运动过程中，单元节点可以分离，且随着单元的平移和转动，允许调整各个单元之间的接触关系(guānxì)。单元之间相互作用的力由力和位移的关系(guānxì)求出，单元的运动则完全根据单元所受的不平衡力和不平衡力矩的大小按牛顿定律确定。块体单元可能达到平衡状态，也可能一直运动下去。共八十七页图2-4-1离散(lísàn)元法未知量示意图（Cundall,1971）(DistinctElementMethod)共八十七页图2-4-2DEM计算(jìsuàn)过程简图共八十七页力—位移法则：假定所有变形发生在角—边接触上，力是由变形产生的。力和位移采用增量，即一个位移变化将产生一个力，且作用(zuòyòng)在接触已存在的力上。在一个时步内，对一个给定的接触，切向和法向位移增量由（u，v，r）增量计算，而新的切向和法向力根据原作用(zuòyòng)力通过力—位移关系求得。然后这些接触力被重新分解为等效的X、Y分力和转动R，并与其他X、Y分力和转动叠加作用(zuòyòng)到每个块体上。共八十七页运动法则：某一块体上，力和运动的总和被用于计算在X、Y和R方向的加速度，通过积分求得速度和位移，有了这个新的位移，就可进行(jìnxíng)下一步计算。在时间域内，DEM采用显式数值方法，块体间的相互作用由接触(jiēchù)和块体的运动检查算出。因此，DEM可以跟踪每个块体每一步的行为。共八十七页图2-4-3DEM块体接触时的力学(lìxué)模型共八十七页DEM采用显式方法求解方程，优点是在解运动方程时，由于假设作用力右端项是已知的，不需要解联立方程。但是法向接触弹簧(tánhuáng)的长度是块体位移的函数，而在接触计算中位移是未知的，因此，在接触发生时，假设是不满足的。第n-1步法向接触弹簧(tánhuáng)的长度产生的法向接触力是第n步假设的力。这种逼近不能产生正确的接触力，有时会破坏能量守恒而使得解答发散。于是，需要引入一个阻尼，使得系统稳定。这样，DEM的接触模型可以表示如下。共八十七页离散单元法的动态(dòngtài)松弛法和静态松弛法动态松弛法是把非线性静力学问题化为动力学问题来求解，用显式中心差分法来近似积分运动方程，并用适当的阻尼来吸收系统的动能，使系统的振动尽可能快地消失，同时场函数收敛于静态值。它按时步在计算机上迭代求解，整个计算过程只需要直接代换，即利用前一步迭代的函数值近似求解新函数值，对非线性问题也可以考虑。静态松弛法是直接寻找块体失去平衡后再次达到平衡时的力－位移关系，采用隐式法联立平衡方程组，并以完全消除块体的残余力和力矩为目标进行(jìnxíng)迭代求解。静态松弛法避免了动态松弛法中的难点，即粘性阻尼的确定及计算时步的选取，但它在求解联立平衡方程组时，有时会碰到数值奇异或病态问题，有待进一步改进。目前，工程中广泛使用的离散单元法多采用动态松弛法。共八十七页5、离散有限(yǒuxiàn)单元法(DFEM)

有限元法在模拟岩体的非连续特性时存在一些缺限，如块体的接触、相对位移、分离等。近10多年来，一些学者对FEM作了一些修正，形成离散优先单元法（DiscreteFiniteElementMethod）主要有两类：①采用特殊单元模拟非连续特性。②采用接触判断模拟非连续特性而不增加附加单元。(1)Goodman(1976)最早提出“jointelement”，用于模拟非连续缝面，就是众所周知的“Goodmanelement”。Nakazawa(1979)利用FEM讨论(tǎolùn)不同摩擦条件下的接触分析，用于研究非线性接触问题。Lei（2001）介绍了一种简单的界面单元用于接触分析，他将节点位移作为未知量，模拟两个具有初始接触或间隙的块体的接触滑动、解耦和再结合。

缺点：根本点都是基于小位移和小变形，难以很好地解决大位移问题。共八十七页(2)为模拟非连续体的大位移特性，出现了利用法向和切向弹簧来进行接触计算，但是块体之间不设特殊单元，有两类：(a)接触预确定FEM和(b)接触搜索算法FEM。(a)

Chaudhary（1986）、LaursenandSimo（1993）给出了一种FEM方式来处理多块体、大变形摩擦接触问题，它采用连续理论，应用通常的FEM离散方法，继承了连续FEM的优点(yōudiǎn)。Solberg等人（1998）利用固体力学中无摩擦接触问题的数值解答，将微分—代数方程的分析思想应用于FEM。这种方法采用预确定接触模式（pre-determinedcontact）来求解问题，可以模拟块体间的接触、相对滑动、分离等。(b)为模拟真正的非连续问题（接触是每步都在变化，且预先未知的），BelytshkoandNeal（1991）提出了一种弹球算法“pinballalgorithm”。该方法将块体视为弹球，并置于块体表面，认为只有在弹球交迭的地方发生接触，从而把复杂形状的接触判断转化为圆形的接触判断。由于过于粗略，Zhong（1993）介绍了一种基于几何形状的接触算法来解决接触/嵌入的大变形问题。共八十七页图2-5-1二维“弹球”模型(móxíng)示意图(b)BelytshkoandNeal（1991）提出了一种弹球算法“pinballalgorithm”。该方法将块体视为弹球，并置于块体表面，认为只有在弹球交迭的地方发生接触，从而把复杂(fùzá)形状的接触判断转化为圆形的接触判断。由于过于粗略，Zhong（1993）介绍了一种基于几何形状的接触算法来解决接触/嵌入的大变形问题。共八十七页6、非连续(liánxù)变形分析(DDA)

石根华和Goodman于1989年提出了一种称为非连续变形分析(DiscontinuousDeformationAnalysis—DDA)的数值方法，它是在非连续体位移反分析法的基础上推广而来的一种正分析方法，它可以从块体结构的几何参数、力学参数、外荷载及约束情况计算出块体的位移、变形、应力、应变以及块体间离合的情况。

DDA将系统视为由常应变并且可作刚性位移的块体所组成，建立起类似于有限元的平衡方程（其矩阵含义与有限元不同），DDA在建立平衡方程时考虑(kǎolǜ)如下两个假定（约束）：①块体间无拉应力作用。②块体间不能有叠合。一旦有拉力或叠合产生，即在相应点上作用一反向弹簧以使拉力或叠合消失。从这两个约束出发，利用最小位能原理，DDA建立了平衡方程。其刚度矩阵由材料特性、约束情况及变形情况所决定。位移列阵包含块体的平动位移、转角及应变，力列阵由各块体所受外力集成。共八十七页图2-6-1DDA中块体(kuàitǐ)的变形共八十七页图2-6-3DDA采用(cǎiyòng)的块体未知量共八十七页平衡(pínghéng)方程共八十七页图2-6-4DDA块体接触(jiēchù)的不同类型共八十七页图2-6-5DDA块体接触(jiēchù)时的力学模型共八十七页DDA与FEM、DEM的比较(bǐjiào)DDA与有限元法（FEM）有相似之处，但它属于离散计算的范畴，它的单元可以任意凹凸，也可以分离，接触面上采用摩尔—库仑准则，对每级荷载以时步求解平衡方程。它的未知量是块体位移，而有限元的未知量是结点(jiédiǎn)位移。DDA与离散元法（DEM）不同，DDA是位移解法，而离散元法通过调整接触力来使系统达到平衡，属力法范畴。此外，DDA是一种隐式解法，而离散元法是显式解法。DDA既可用于静力问题，也可用于动力问题，已在一些工程中得到应用。目前二维程序比较成熟，但模型比较简单，应用还不十分广泛。三维DDA还仅仅处于开发阶段。DDA只是近几年才传入我国，有待于进一步发展，是很有潜力的一种算法。共八十七页7、数值(shùzí)流形方法（NMM）1995年，石根华提出了DDA与FEM的统一形式：数值流形(liúxínɡ)方法(NMM)。NMM以数学流形为基础，计算结构体的位移和变形。NMM中的网格就是数学覆盖，这些数学覆盖互相重叠并且覆盖整个计算区域，在每个数学覆盖上定义互相独立的位移近似函数。这些数学覆盖被物理边界切割而形成物理覆盖，物理覆盖的重叠区域形成单元。将这些覆盖上的位移函数结合起来形成计算域上的全域位移近似函数，在每个单元上的近似函数就是形成此单元的若干个互相重叠的覆盖上的近似函数的加权平均，并利用最小位能原理形成整体平衡方程。在NMM中，积分方法采用了单纯形上的解析积分形式，这是与传统数值方法的一个很大的不同。共八十七页图2-7-1NMM的求解(qiújiě)过程共八十七页NMM基本概念——四边形覆盖(fùgài)图2-7-2NMM例子(lìzi)示意图共八十七页三角形数学(shùxué)覆盖含有裂隙的边坡覆盖系统(xìtǒng)（1）含裂隙的边坡（2）覆盖系统（3）流形单元

共八十七页共八十七页NMM中全局函数(hánshù)的形成共八十七页NMM与DDA、FEM的关系(guānxì)NMM使用数学覆盖系统，使得连续体、非连续体的整体平衡(pínghéng)方程都可以用统一的形式来表达。有限元法FEM和非连续变形分析DDA都是NMM的特殊形式，可以说NMM使DDA和有限元法FEM在理论基础上和表达形式上得到了统一，这是一个很大的贡献。目前，NMM中的覆盖只能采用一些特殊的形状，如正三角形、正六边形等，而实现任意形状的覆盖则有一定的困难。因为虽然在理论上可以采用任意形状的覆盖，但真正操作起来，对任意形状的覆盖，权函数的选取及解析积分的求解都有一定的困难。在这些方面，还有很多工作要做。所以，目前工程界使用较多的仍然是DDA。共八十七页NMM

DDA共八十七页NMM

FEM共八十七页NMM

FEM，三角形单元(dānyuán)共八十七页NMM

FEM，四边形单元(dānyuán)共八十七页8、无单元(dānyuán)法（EFM）无单元法属于连续介质的范畴，但是，它可以比较方便地处理裂缝扩展跟踪问题，而且，也是近来的研究热点。因此，这里作一简单介绍。无单元法的思想(sīxiǎng)最早由Nayroles等人于1992年提出，他们称之为虚拟单元法（DiffuseElementMethod）。Nayroles的近似比较简单，使用了低阶积分，边界条件的引入也不准确，而且在拟合函数的求导过程中忽略了一项，计算比较粗糙。但是，他们第一次在偏微分方程的数值解法中引入了滑动最小二乘法的思想。事实上，滑动最小二乘法很早就有人提出来了，如Lancaster等，只是在Nayroles之前，它只被应用于曲线、曲面拟合领域。共八十七页滑动(huádòng)最小二乘法原理共八十七页共八十七页正规(zhèngguī)方程（法方程）共八十七页共八十七页无单元(dānyuán)伽辽金法（EFGM）Belytschko等人于1994年对Nayroles的方法进行了改进，提出了无单元伽辽金法（Element-freeGalerkinMethods），主要作了如下改进：将计算(jìsuàn)域上的积分分离为很多子域上的积分，在每个子域上采用高阶高斯积分。在拟合函数的导数中加入被Nayroles忽略了的项。用拉格朗日乘子法引入边界条件。在此基础上，Belytschko从变分原理出发，建立了类似于有限元的无单元法基本方程。Belytschko方法的一些不足：首先，在求解滑动最小二乘拟合函数时，需要求解线性方程组，因此，在集成刚度矩阵时，在每个高斯点上都要求解一个线性方程组，在求位移、应力、应变时也是如此，这是非常繁琐的，而且容易造成大的误差。其次，用拉格朗日乘子引入边界条件后，增加了未知量，也破坏了刚度矩阵的正定性及带形分布。增加了方程求解的时间和难度。共八十七页数值(shùzí)方法的精度比较xy50x1x2x3x4y3c.y1y2y4xx50x1x2x3x4y3d.xx50x1x2x3x4b.a.具有一个非连续(liánxù)点的一维函数，细线是函数真解，粗线是数值逼近(a)FDM逼近(b)FEM逼近，一维单元：x0x1,x1x2,x2x3,x3x4,x4x5(c)DEM/DDA逼近，一维块体：y0x1,y1x2,y2x3,y3x4,y4x5(d)NMM逼近，一维物理覆盖：U1=x0x1,U2=x0x2,U3=x1x3,U4=x2x3,U5=y3x4,U6=y3x5,U7=x4x5由于x3y3附近的突变，数学覆盖x2x4被分成2个物理覆盖（U4和U5）共八十七页三、模拟裂缝(lièfèng)扩展DEM、DDA与模拟裂缝(lièfèng)扩展NMM与模拟裂缝扩展EFM与模拟裂缝扩展共八十七页DEM、DDA与模拟(mónǐ)裂缝扩展非连续介质的数值方法，如离散元法DEM、非连续变形分析DDA，都可以用来模拟岩体的开裂，甚至可以模拟结构大位移、大转动。但是，由于受模型本身的限制，它们都有其局限性：1.块体是由岩体节理裂隙充分切割而形成的，块体间只简单地满足摩尔－库仑准则，对于结构面已将岩体完全切割成块体状的情况，采用这一理论是合适的，如用来处理节理裂隙发育地区的地下结构工程的塌方与支护问题，是相当实用和有效的。但当岩体并未被裂隙切割成块体集合时或者当裂隙不发育时，这些方法就无优势可言了。2.模拟开裂时，受块体形状所限，开裂面只能沿块体边界，而由于裂隙扩展方向事先是无法知道的，因此，用这种方法计算岩体开裂，将随着不同的块体结构得到不同的破坏形式，很难得到符合实际的结果。可见(kějiàn)，DEM、DDA在模拟切割完整的结构面的破坏时是非常有效的，但若要裂缝扩展则有一定的困难。共八十七页NMM与模拟裂缝(lièfèng)扩展NMM方法使连续介质与非连续介质的计算从理论上得到了统一，这是一个很大的贡献。NMM在计算连续体与非连续体的大变形或进行动力分析时非常有效。在开裂跟踪方面，NMM也得到了应用(yìngyòng)。但是，由于NMM的近似函数在数学覆盖上定义，由物理边界切割数学覆盖而成的物理覆盖互相重叠形成单元，单元上的近似函数为各物理覆盖近似函数的加权平均，因此，虽然NMM由于采用了可以不变的数学覆盖而比传统的有限元法更为灵活，但它实际上仍未摆脱单元的限制，在物理边界发生变化后，单元也应作相应的调整。在开裂跟踪时，NMM有以下局限。共八十七页NMM在模拟开裂(kāiliè)时的局限1.为了提高裂尖的计算精度，通常需要在裂尖加密网格。NMM在这种情况下就需要加密裂尖的数学覆盖，如果数学覆盖不变，则达不到要求的计算精度，但加密数学覆盖与有限元的重新划分网格是同样繁琐的。2.在数学覆盖不变的情况下，NMM的开裂长度应贯穿(guànchuān)某个数学覆盖或物理覆盖。若非如此，开裂后NMM的物理覆盖仍不会改变，则单元也不会改变，其结果就是新开裂的裂缝得不到模拟。因此，目前用NMM计算开裂时，多让裂缝沿单元边界贯穿(guànchuān)单元，以简化单元的重新调整过程。这在模拟受力简单的裂纹时是可以的，但对于岩土工程中复杂受力条件下的裂缝，由于事先无法确定开裂方向，就会产生较大的误差，难以得到符合实际的结果。共八十七页EFM与模拟裂缝(lièfèng)扩展无单元法（EFM）采用滑动(huádòng)最小二乘法拟合场函数，因此，在计算中只需计算域边界条件和结点，不需要任何单元信息；另外，由于滑动(huádòng)最小二乘法拟合函数具有高阶连续的特点，因此，无单元法具有前后处理简单、计算精度高的优点。这使得无单元法特别适用于岩土工程中的稳定、开裂分析，特别是在开裂问题中，无单元法可以很好地模拟裂缝尖端奇异场，算出高精度的应力强度因子，较为方便地跟踪裂缝扩展，而不存在重新划分网格的问题。共八十七页四、算例和应用(yìngyòng)非连续变形(biànxíng)分析（DDA）数值流形方法（NMM）无单元法（EFM）共八十七页1、DDA：边坡稳定(wěndìng)分析图1-1(a)DDA概化模型(móxíng)（日本）下小鳥ダム湛水池不連続岩盤斜面の安定性評価

共八十七页图1-1(b)DDA计算(jìsuàn)的边坡不稳定的情况共八十七页1、DDA：边坡落石分析(fēnxī)图1-2(a)9号国道A017断面(duànmiàn)边坡坍塌安全分析评价的DDA模型（Step=0）（日本）国道９号美方郡温泉歌長、養父郡関宮町尾崎斜面落石、崩壊検討共八十七页图1-2(b)Step=142共八十七页1、DDA：边坡稳定(wěndìng)分析图1-3(a)金居原水电站水库边坡潜在滑动面稳定分析的DDA模型(móxíng)（Step=0）（日本）金居原発電所潜在的岩盤地すべり安定率研究

共八十七页图1-3(b)Step=800共八十七页1、DDA：边坡稳定(wěndìng)分析图1-4(a)殿山水电站水库岩体边坡滑动分析的DDA模型（SlidingFace:B）（日本）殿山発電所下露地(lùdì)すべりDDA解析、安定性評価

共八十七页图1-4(b)滑动体的位移随滑动面上接触参数(cānshù)变化的情况共八十七页1、DDA：洞室围岩(wéiyán)稳定分析图1-5(a)龙滩水电站地下厂房(chǎngfáng)开挖的DDA模型（开挖前）仅考虑断层的情况共八十七页图1-5(b)开挖(kāiwā)后共八十七页图1-6(a)开挖(kāiwā)前考虑断层和节理的情况：节理平均密度为1条/2.5m和1条/5.8m共八十七页图1-6(b)开挖(kāiwā)后共八十七页2、NMM：边坡稳定(wěndìng)分析图2-1NMM概化模型（其中的三角形是数学覆盖(fùgài)）（日本）下小鳥ダム湛水池不連続岩盤斜面の安定性評価Step=100，稳定情况共八十七页图2-2边坡滑动(huádòng)共八十七页图2-3边坡失稳共八十七页图2-4拱座变形(biànxíng)失稳共八十七页图2-5挤压(jǐyā)变形破坏共八十七页图2-6裂缝(lièfèng)扩展共八十七页3、EFM：半无限平面(píngmiàn)受集中力作用考虑平面应变情况。分别用有限元法及无单元法进行计算。有限元采用96个节点、25个8结点矩形(jǔxíng)二次单元，无单元法采用同一组节点，积分域为10×10个1m×1m的矩形子域，每个子域采用2×2高斯积分，基函数选用一次基。图3-1(a)半无限平面受集中力作用共八十七页半无限平面(píngmiàn)