非连续数值方法综述

非连续数值方法综述杨凡（河海大学水利水电学院，江苏南京210098）摘要：非连续问题是岩土及水利工程中不可避免的一类难题，由于其对工程的影响巨大，近几百年来特别近一个世纪以来一直是工程界研究的一个热门话题。

从最早的非连续问题解析解法—刚体极限平衡法出发，引申出近几十年来有关非连续问题研究的热点—非连续问题的数值解法，然后对这些非连续的数值方法的基本原理和实际应用发展情况进行一一综述。

关键词：非连续；数值方法；岩石和土都是经历过变形的地质体，受其成因、组成、结构、年代等诸多因素的影响，岩土材料具有高度的非连续性、非均匀性和各向异性的特征，在力学性质上表现出强烈的非线性。

岩土工程是一门综合应用岩石力学、土力学、工程地质学等基本知识解决实际工程中有关岩体与土体变形及稳定问题的学科[1]。

岩土工程中的非连续变形问题主要是由岩石及土体中不连续面的存在引起的，岩土工程问题中的不连续面大致可分为两类，一类是指存在于岩体中的节理、软弱夹层以及土体中的剪切破坏面，另一类则是岩土结构如各类基础、挡土结构、地下结构等与岩土体之间的接触面。

显然，不连续面对岩土体或结构的受力、变形有着重要的影响，因此为使计算结果真实地反映出岩土体及结构的受力和变形情况，在计算时不能忽视不连续面的存在[2]。

对于具有不连续面的结构，在承受荷载的过程中，不连续面的状态是在不断变化的，这将影响到两侧岩土体的应力和变形，从而影响到整个体系的应力场，而应力场的改变又影响到不连续面的状态。

因此，解决岩土力学问题的关键在于对非连续变形的模拟，分析研究结构中各种不连续面的构造特点和力学性能，研究其受力状态的变化规律及其对结构整体性能的影响是工程设计中的关键研究课题之一，具有很大的学术意义和实用价值[3]。

几百年来，人们对非连续变形问题作了大量的研究工作。

最早有关非连续问题的研究主要集中在寻求解析解的层面上。

1773年，法国科学家库伦在大量实验基础上总结了著名的库伦土压理论，刚性楔体和静力平衡的应用也为后续研究奠定了一个基调。

在此基础上，瑞典圆弧法的提出表明刚体极限平衡方法的正式诞生，在接下来的几十年间该方法得到了很好的发展并出现了Bishop、Janbu、Sarma等改进的方法。

在发展过程中，极限平衡法在水利工程上得到了很好的运用。

但其缺点亦非常明显，即不考虑土体的变形，得到的结果偏安全。

在许多的工程中，其运用受到局限。

近年来，伴随着计算机技术的高速发展，有关非连续问题的数值分析逐渐发展成为热门。

非连续分析的数值方法主要包括：界面单元有限元法(FEM with interface elements)、刚性有限元法（RFEM）、离散单元法(DEM：Distinct Element Method)、非连续变形分析(DDA：Discontinuous Deformation Analysis)、数值流形方法（NMM：Numerical Manifold Method）、无单元法（(Element—Free Method）、耦合方法（Couple Methods）以及渐进破坏模型(Progressive Fracture Models)等[4]。

以下对各个方法研究进展情况进行综述。

1 界面单元有限元法(FEM with interface elements)有限单元法在连续介质上面的成功应用，使其广泛地被人们接受，但在解决前处理问题、应力与应变解答不连续问题计算等方面还存在着一些局限。

为使其能够处理简单的非连续性问题，人们提出了各种能够反映非连续性质的简单力学元件或特殊界面单元用于模拟单一非连续界面的力学行为[4]。

这类方法主要包括有：联结单元(Linkage Element)（Ngo和Scordelis，1967)，无厚度接触单元法或节理单元法(Joint Element)(Goodman，1968)、薄层单元法(Thin-Layer Element)(Desai，1982)、接触—摩擦单元(Contact—Friction Interface Element)(雷晓燕和swoboda，1994)，其中联结单元法与节理单元法是在性质差异悬殊的两相邻单元间界面上设置法向与切向弹黄，一旦法向弹黄受拉则视界面发生分离，若切向相互作用力超过界面的极限强度则认为界面产生相对滑移。

而薄层单元法假定界面是由特殊介质所组成的厚度很薄的实体单元，这种单元的本构性质取决于周边介质力学性态，由试验直接测取。

通常认为界面的法向性质与周围岩土材料相应特性相同，Desai与Nagaraj(1988)对此曾提出采用非线性方程来考虑法向应力—应变关系。

关于切向响应特性，Desai等研制了一套循环直剪试验仪来测定界面剪力与相对滑移距离之间的依赖关系，并建议用多项式函数或五参数Ramberg —Osgood模型加以描述.所谓接触—摩擦单元与接触单元的区别在于它采用节点接触应力代替节点接触力，对于接触面的几何方面与静力学方面约束条件用附加条件纳入到系统方程中，从而能够通过计算自动满足接触界面的滑移、分离以及粘结真实状态。

这类方法将岩土介质视为准连续介质，仍以连续分析为主，但可以对个别具有控制作用的宏观非连续面的变形与破坏等力学效应给予重点分析，在工程中得到了一定应用，但由于下述自身的缺陷而受到了限制，在处理复杂的非连续变形问题时显得无能为力[4]。

例如，(l)只能对原生的非连续界面进行计算，对次生的非连续界面无法处理；(2)界面单元数目不能设置得太多；(3)界面弹簧刚度的选取较为困难。

2 刚性有限元法（RFEM）刚性有限元法最早由日本Kwaai教授于1976年提出，该方法将离散后的块体视为刚体，块体之间用界面上的法向弹簧和切向弹簧相联结，以块体形心处的刚性位移为基本变量，用分片的刚体位移模式逼近实际整体位移场，以块体间的联结弹簧反映结构内部的弹性，并用界面应力表征结构内部的应力[5]。

刚性有限元中的模型主要包括刚体—弹簧模型(RBSM)、刚性有限元、分块刚体位移—界面应力元和块体—弹簧模型等。

此类模型的界面特性均服从Coulomb摩擦定律，对体系的静力学约束条件考虑得比较充分，在连续状态的应力分析方面可给出较高的计算精度，对于临界状态能够估算出极限荷载，并可有效的用于少量块体界面间的摩擦接触分析。

但这类方法过分强调岩土体结构面的作用，对结构体的变形没有给予足够重视。

虽然这种方法可以用于原生界面的破坏分析，但不能模拟实际岩土体的破坏发展过程，也无法模拟破坏发展导致的次生界面的非连续变形行为和块体失稳后的运动过程。

3 离散单元法(DEM：Distinct Element Method)离散单元法(DEM：Distinct Element Method)是由美国Cundall等(1971)提出，它特别适用于含有结构面的结构的应力分析，最初用来分析岩石边坡的渐近破坏[6][7]。

它是以软弱结构面切割而成的离散块体为基本单元，其几何形状取决于岩土结构中非连续面的空间位置及其产状，在块体间的接触约束下运用牛顿运动定律描述各块体的运动过程。

离散元法有动态松驰法和静态松驰法两种。

目前离散单元法大都采用动态松驰法。

动态松驰法是把非线性静力学问题转化为动力学问题来求解，用显式中心差分法来近似地对运动方程进行积分计算，并假设块体在运动时动能将转化为热能而耗散掉，在计算中引入人工粘性阻尼以使系统达到平衡、运动趋于稳定[8]。

经过近40年的发展，离散单元法已成为模拟岩土体非连续大变形的有效的数值方法之一。

目前，二维、三维离散单元法均趋于成熟，主要表现在商业化软件己经进行入实用化阶段。

开发离散元商用程序最有名的公司要属ITASCA公司，该公司分别开发了二维的UDEC和三维的3DEC块体离散元程序，在模拟节理岩石在准静或动载条件下力学过程及采矿过程的工程问题上有很好的应用。

另外，颗粒流程序PFC2D和PFC3D对破损累计导致的破裂、动态破坏和地震响应等问题有很好的模拟。

国内有关离散元的研究起步较晚，王泳嘉[9]首次向我国岩石力学与工程界介绍了离散元法的基本原理及几个应用例子。

但是发展很快，东北大学开发了用于土木工程设计的离散元法软件系统2D- Block、三维离散单元软件TRUDEC[10,11]。

现在很多大学和研究所均有人从事离散元法的研究和应用。

工作,成果显著。

4 非连续变形分析(DDA：Discontinuous Deformation Analysis)非连续变形分析(DDA：Discontinuous Deformation Analysis)方法是继离散单元法之后，从20世纪80年代末以来发展起来的一种更新的模拟散体系统力学响应的数值分析方法。

非连续变形分析方法用隐式方程，它类似于有限元，所不同的是它引进了运动方程，用最小势能原理把块体之间的接触问题和块体本身的变形问题统一到矩阵的求解，具有完备的运动学理论、严格的平衡假定、正确的能量消耗[12,13]；理论严密，精度较高，而且把静力和动力、正分析和反分析统一起来，不仅可以计算破坏前的小位移，也可以计算破坏后的大位移，对滑动、崩塌、爆炸及贯入等问题十分有效。

在DDA提出后的几年内，其立刻得到了很快的发展，成果集中在二维DDA上，主要体现在以下两个方面：(l)块体接触判断计算方法的改进；(2)加入块断裂力学，考虑块体破裂；石根华(1988)[12]最早用罚函数法来处理非连续面上的接触问题,该方法有其优缺点。

优点是采用罚函数总体平衡方程的方程数量不增加而且容易得到解答，但缺点也很明显，即接触问题的解答精度与所选取的罚数关系密切，罚函数方法仅能近似地满足接触限制并且接触力必须先计算出来。

为克服此缺点，Lin(1994)[14]最早应用扩展Lagrange乘子法来代替罚函数法进行接触处理。

后来，为了更准确地求解块体系统的接触力，我国学者蔡永恩(1998)[15]采用Lagrange乘子法和分区解法，在石根华(1988) 的DDA的基础上提出了Lagrange非连续变形分析，栾茂田、黎勇(1999)[16]提出基于接触力元的非连续变形分析方法。

这些方法的特点是将接触力或接触界面的流动距离作为未知参数，纳入总体方程进行求解，从而得到接触界面上的真实接触力。

但由于接触力或流动距离等参量增加了总体方程个数，增加了总体方程求解的计算量。

当考虑的块体较多时，这种情况会更为严重。

这些方法的总体方程一般不对称[4]。

在石根华(1988)的DDA中，没有考虑块体本身的破裂问题，而在实际岩体工程中，岩体中裂纹的扩展以及块体在受力过程中的破坏必须加以考虑。

Ke(1993)[17]为了计算块体内的各种应力状态以及模拟完整块体中裂纹的扩展，发展了基于人工节理的非连续变形分析方法(Artificial Joint-Based DDA)。