非连续数值方法综述

非连续数值方法综述
杨凡
（河海大学水利水电学院，江苏南京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)。

Gebra(1994) [18]主要研究了块体在瞬间碰撞条件下的块体的破裂算法。

截止目前，二维DDA已发展形成了比较完善的理论体系，且在实际工程中得到了很广泛地使用。

但是，二维DDA法的研究和应用大多局限于平面应变或平面应力问题，并不能完全满足实际工程应用的要求，不能完全反映实际的变形和特性，有很多的局限性[19]。

所以近年来，国内外一些学者正致力于三维DDA的研究，并取得了很大的进步。

在接触判断方面，姜清辉（2000,2002）[20,21]分别运用罚函数法和矢量法推导了摩擦接触问题和角-面接触模型的计算方法，并给出其法向子矩阵公式和固定接触子矩阵以及摩擦接触子矩阵公式。

张旭（1999）[22]将柔性弹簧引入DDA，讨论了柔性弹簧和刚性弹簧的共同作用，并验证了其合理性。

在工程应用方面，王书法（2000）[23]用DDA法，将温度变化对岩体的力学影响转化为等效初应力的作用来求解节理岩体热应力问题，是对DDA理论体系的补充，增强了DDA方法解决实际问题的能力，拓广了DDA方法的使用范围。

靖洪文（2003）[24]将DDA法应用到地下工程中，定量研究了非连续围岩体位移影响因素的变化规律。

沈振中（2007）[25]提出了基于三维非稳定饱和-非饱和渗流分析的岩体边坡稳定非连续变形分析方法,对于研究降雨入渗条件下岩体边坡的稳定性分析和安全评价具有重要的理论和实践意义。

4 数值流形方法（NMM：Numerical Manifold Method）
流形方法是石根华(l992)应用流形的覆盖技术建立的一种把有限元法、非连续变形分析方法和解析方法包含在内的全新的统一计算方法。

在继承了DDA的接触机制和块体定义方法(即物理网格)的基础上，流形元法设置了专门用于定义块体内部插值函数的数学网格。

这样，流形元法不仅可以像DDA那样有效地模拟不连续面的接触和块体的运动，而且可以像有限单元法那样精确地求解块体内部的应力分布。

有关数值流形方法的研究主要集中在最近的几年，Terada（2004）[26]在数值流形方法理论方面进行了系统的研究，结合平面弹性变形问题对一系列不规则形状流形单元的性能进行了评估，并与有限元进行对比。

Terada（2005）[27]还研究了一般连续材料和复合材料弹塑性变形分析的数值流形方法，他将非标准几何形状的边界流形单元分解为多个三角形的组合，并采用3个高斯积分进行单元矩阵计算。

对于三维数值流行方法的研究也取得了很大的进展，骆少明（2005）[28]研究了三维数值流行方法理论，采用标准六面体数学网格构建有限覆盖系统，推导了单元刚度矩阵，详细研究了三维数值流行的Hammer积分及剖分规则，建立了线弹性变形三维数值流形分析的理论体系。

郑榕明（2004）[29]研究了基于八节点六面体等参单元形式的三维数值流行方法，推导了相应的单元矩阵，引入向量理论和迭代方法解决三维接触问题。

另外，姜清辉、彭自强等人均作了很多的研究。

但总体而言，目前三维和高阶流形方法的研究主要还是针对线弹性问题，且基本应用于岩石变形和静力结构上面。

5 无单元法（Element—Free Method）
无单元法是20世纪90年代发展起来的一种新的数值方法。

这一方法的主要思想是近似的场变量采用形函数的线性组合，在建立形函数时不借助于区域的网格，代之为在区域内设置结点并以其
一定的权函数使得子域内结点之间的联系，对某一结点的形函数的建立考虑这一结点权函数的影响区域，而没有必要像有限单元法中通过单元建立结点之间的联系[30,31]。

由于无单元法仅要求计算域的几何边界及一些计算点，而不需要也没有单元信息，从而摆脱了有限单元法中的单元限制，可大大简化有限单元法中的前处理工作。

无单元法中的几种主要算法有：光滑微湍水动力算法、有限点法、单位分解法、无单元Galerkin法、弥散单元法、局部P无网格法、再生核微团法、局部边界积分方程法，其中应用较多且较成熟的是无单元Galerkin法[32]。

在无单元法的应用研究方面，凸显了更深更广的方向趋势。

丁丽宏等[33]在广义塑性力学中应用无单元法计算了土体在均布荷载和偏荷载共同作用下的弹塑性变形；龙述尧等[34]用无单元法解决弹塑性问题；介玉新等[35]将无单元法与非线性有限元耦合应用于面板堆石坝的应力应变计算，在面板和防渗墙等结构上布置无单元节点，避免了有限元计算网格划分长宽比例过大的缺陷；陈虹等[36]通过引入滑动最小二乘法和有限差分法，得到水动力学无单元计算法并应用于复杂边界的河道水流运动方程；胡德安[37]采用了增量和修正的Newton-Raphson迭代分析的法，并在整个分析过程中所有变量的表达格式都采用全拉格朗日格式，从而实现了用无单元法求解几何非线性问题；胡云进、周维垣等[38-39]编制了三维无单元程序，并能够用于模拟拱坝三维开裂问题。

6 耦合方法（Couple Methods）
耦合方法是指对两种或者两种以上的方法进行结合运用得到的方法。

如有限元与离散元的耦合、离散元与边界元的耦合，这些耦合方式可以发挥各种方法各自的优点，以达到对具体工程的适应和良好应用。

7 渐进破坏模型(Progressive Fracture Models)
渐近破坏模型是通过对材料微结构的演化和局部变形行为的研究，对材料在变形、损伤与破坏时所产生的本质现象和因果关系进行模拟的方法[40-41]。

与各种已有的数值计算方法不同，它注重对材料的细观断裂机理与断裂规律的研究，侧重对事物内在机理的模拟。

渐进破坏模型是在对复合材料和岩石等脆性材料的破坏研究中提出的。

这类方法一般包括两个方面，即应力分析和破坏分析。

应力分析可以采用解析方法和有限元法，破坏分析是根据一定的破坏准则检查材料结构中是否有单元破坏。

对于岩石等脆性材料而言，东北大学岩石破裂与失稳研究中心研制开发出了相应的岩石渐进破坏分析软件RF2AP。

该软件通过考虑材料的非均匀性来模拟材料非线性，通过单元的弱化来模拟材料的变形和破坏等非连续力学行为。

8 结语
当前我国水电资源的开发主要集中在西南，地质复杂的高边坡及地下工程屡见不鲜。

怎样最合理地对岩土问题中的非连续问题进行模拟是近几年研究的一个热门。

伴随着锦屏一级、二级，白鹤滩等工程的进行，对非连续问题的有效分析研究又显得刻不容缓且意义重大。

以上所述的各种非连续数值方法目前还处于发展阶段，很多好的模拟方法理论体系还不够完善，例如三维DDA和数值流行方法；很多理论体系已较为完善，但是在工程的应用中还是存在很大的发展空间，例如离散元法和无单元法。

但是主要集中在理论方面，实际工程中的应用更待进一步的探索和研究。

伴随着近年
来对工程考虑的因素越来越多，多场耦合已然成为工程研究的必然趋势；在这种环境下，发展多种方法的耦合方法将是今后研究的重点。

参考文献：
[1]孔德森，栾茂田.岩土力学数值分析方法研究[J].岩土工程技术. 2005,19(5)：249-253.
[2] 毛坚强. 接触问题的一种有限元计算方法及其在岩土工程中的应用[D].西南交通大学博士学位论文. 2002.
[3]张燎军. 水工结构接触问题的力学模型及其在三峡工程中的应用[D].河海大学博士学位论文.2005.
[4]田荣. 连续与非连续变形分析的有限覆盖无单元方法及其应用研究[D].大连理工大学博士学位论文.2000.
[5]张雄. 边坡稳定性的刚性有限元评价[J].成都科技大学学报.1994，80(6)：47-52.
[6] 鲍鹏，姜忻良，崔奕.可变形体离散元模型[J] .天津大学学报.2005，38(6)：552-555.
[7] 李世海，高波，燕琳..三峡永久船闸高边坡开挖三维离散元数值模拟[J] .岩土力学.2002，23(3)：272-277.
[8 ]王濂嘉，邢纪波.离散单元法及其在岩土力学中的应用[M] .沈阳：东北工学院出版社，1991.
[9] 王泳嘉.离散单元法—一种适用于节理岩石力学分析的数值方法[A].见：陈祖煜编.第一届全国岩石力学数值计算及模型
试验讨论会论文集.江西吉安:西南交通大学出版社,1986 .
[10]王泳嘉,宋文洲,赵艳娟.离散单元法软件系统2D-Block的现代化特点[J].岩石力学与工程学
报,2000,19(supl.):1057-1060.
[11]王泳嘉,刘连峰.三维离散单元法软件系统TRUDEC的研制[J].岩石力学与工程学报,1996,15(3):201-210 .
[12]Shi G H. Block System Modeling by Discontinuous Deformation Analysis[M] .Southampton：Computational Mechanics
Publications，1993.
[13]石根华.数值流形方法与非连续变形分析[M].裴觉民译.北京：清华大学出版社，1997 .
[14]Lin C T，Amadei B and Sture S. Using an augmented Lagrangian method and block fracturing in the DDA method. Computer
Methods and Advances in Geomech.，1994，837-842.
[15]蔡永恩，梁国平，殷有泉.岩体滑波动力过程的数值模拟新方法[A].第六届全国岩土力学数值分析与解析方法讨论会议
论文集，广州，1998.
[16]栾茂田，黎勇，杨庆.非连续变形计算力学模型在岩体边坡稳定性分析中的应用[J].岩石力学与工程学报，2000，19(3)：
289-294.
[17]Ke T C. Simulated testing of two dimensional heterogeneous and discontinuous rock masses using discontinuous deformation
analysis .A dissertation submitted in Partial satisfaction of the requirements for the degree of Doctor of PhilosoPhy，Berkeley: University of Califomia，1993.
[18] Gebara J M. The finite block method: Its basis and its modification to allow the fracturing of blocks under high impact loads. A
dissertation submitted in Partial satisfaction of the requirements for the degree of Doctor of PhilosoPhy，West Lafayette : Purdue University，1994.
[19]唐瑜，苏超，姜媛媛.三维非连续变形分析方法及其研究应用进展[A] .水电站设计，2008，24(4)，89-93.
[20]姜清辉，丰定祥.三维非连续变形分析方法中角—面接触模型的研究[J].岩石力学与工程学报，2000，19(增刊)：930-935.
[21] 姜清辉，杨文柱，吴益民，孙宁.三维非连续变形分析方法中摩擦接触问题的研究[J].岩石力学与工程学报，2002，21(增
刊)：2418-2421.
[22] 张旭，于建华.柔性弹簧与刚性弹簧在DDA方法中的应用[J] .四川联合大学学报(工程科学版)，1999，3(30):122-128.
[23] 王书法，朱维申，郑少河.节理岩体热应力问题的非连续变形分析方法[J].岩石力学与工程学报，2000，19(5)：643-646.
[24] 靖洪文，许国安.地下工程破裂岩体位移规律数值分析[J] .岩石力学与工程学报，2003，22(8)：1281-1286.
[25] 沈振中，姜媛媛，李金海.降雨入渗非饱和渗流影响的岩体边坡稳定非连续变形分析[DB/OL].中国科技论文在线(http://
www.paper /downloadpaper.php? serial-number=shenzhen- zhong-2007-03-05).
[26] Terada K, Kurumatani M. Performance assessment of generalized elements in the finite cover method. Finite Elem Anal Des,
2004, 41(2): 111-132.
[27] Terada K, Kurumatani M. Finite cover method with mortal elements for elastoplasticity problems. Comput Mech,
2005,36(1):45-61.
[28] 骆少明，张湘伟，吕文阁等.三维数值流行方法的理论研究.应用数学和力学，2005,26（9）:1027-1032.
[29] 郑榕明，张勇慧.基于六面体覆盖的三维数值流行方法的理论探讨与应用.岩石力学与工程学报，2004,23
（10）:1745-1754.
[30] 王燕昌，李进.无单元方法的研究进展[M].宁夏大学出版社.2003.
[31] 胡云进，朱智兵，周维垣.无单元法对三维不连续面的模拟[J] .岩石力学与工程学报.2004，23(18):3127- 3131.
[32] 陈小虎，沈振中.无单元法的工程应用进展[J].2006,26(6):90-93
[33] 丁丽宏,刘学文,郑璐石.基于广义塑性力学的无单元法[J].2003,23(3):214-216
[34] 龙述尧,陈莘莘.弹塑性力学问题的无单元伽辽金法[J].工程力学,2003,20(4):66-70.
[35] 介玉新,高波,李广信.无单元法在混凝土面板堆石坝应力分析中的应用[J].清华大学学报:自然科学版,2003,
43(8):1092-1095.
[36] 陈虹,林建华.用无单元法求解河道水流运动方程[J]. 计算力学学报,2001,18(2):250-252.
[37] 胡德安.求解几何非线性问题的无单元迦辽金法[D].长沙:湖南大学,2003.
[38] 胡云进,周维垣,寇晓东.三维无单元法及其应用[J].岩石力学与工程学报,2004,23(7):1136-1140.
[39] 周维垣,黄岩松,林鹏.三维无单元伽辽金法及其在拱坝分析中的应用[J].水利学报,2005,36(6):644-649
[40] 刘爱华，王思敬.平面坡体渐近破坏模型及其应用[J].工程地质学报.1994，2(l):l-8.
[41] 刘爱华，颜荣贵.坡体平面渐进模型分析研究[J].湖南有色金属.1994，l0(1):5-8.
Review on Discontinuous Numerical Methods
YANG FAN
( College of Water Conservancy and Hydropower Engineering , Hohai University , Nanjing 210098, China) Abstract:The discontinuous problem in geotechnical projects and water conservancy projects is a kind of difficult and inevitable engineering problem, as its huge impact, nearly hundreds of years especially for the latest century it has been a hot topic in engineering research. Based on the earliest analytical solution,limit equilibrium method,the discontinuous numerical methods for discontinuous problems are given which is the hot point in recent decades, then the development of basic principles and applications of those discontinuous numerical methods are given.
Key words: discontinuous; numerical methods。