（岩土工程专业论文）岩土介质静态松驰界面分析方法和程序及其应用.pdf

摘。要 随着计算机科学的发展，依赖于计算机科学的数值计算在岩土工程 研究领域愈来愈得到重视，数值计算和分析已成为解决各类岩土工程问 题的重要和有力的手段。本论文在对目前国内外岩土数值分析方法的综 合分析基础上，研究了基于静态松弛的界面分析方法、相关算法程序， 并对其在岩土工程中的应用问题进行了探索。 论文首先论述了国内外岩土工程松弛数值方法的研究与发展。之后 从力学的角度来分析和介绍了静态松弛力分配法和静态松弛界面单元法 的基本原理，给出了力分配法的算例，与有限元计算做对比，验证了力 分配法的正确性和可行性。应用界面单元法分析了悬臂梁受拉和受压时 的界面变形问题，研究了悬臂梁的变形情况。 研制了基于静态松弛界面单元法的应用程序，并针对该程序开发了 配套的前处理程序和结果显示的后处理程序。前处理程序包括模型绘制， 网格划分，参数输入和约束施加。提出了一种满足一定约束条件下任意 平面三角网格自动划分的简化方法，研究了其相应的算法。该方法可以 根据所给出的尺寸和边界进行网格划分，自动记录节点坐标和单元编号， 将结果进行输出和显示，并给出了算例；后处理程序可以观察单元和节 点的位移、应力大小、过程的动态显示。文章给出了程序研制思路，关 键步的算法和流程图。 应用所研制的程序对岩土工程中的静力压桩问题进行了深入研究， 分析了轴力和桩侧摩阻力的变化规律。从桩顶荷载，桩体下沉的关系等 方面对静压桩的荷载传递机理进行了探讨。 关键词：静态松弛；界面；力分配法；三角网格；静压桩 a b s t r a c t w i 也m ed e v e l o p m e n to fc o m p u t e rs c i e n c e ，t h en u m e r i c a lm e t h o d sd 印e n d i n go n e o m p u t e rs c i e n c ea r eb e i n gp a i dm o r e 绷l dm o r e 甜e n t i o ni nt h eg e o t e c l l l l i c a le n g i n e e r i n g f e s e a r c hf i e l d ，柚dh a 、，ea l r e a d yb e c o m ei m p o r t a i l ta i l de 仃e c t i v em e a n sf o rs o l v i n gv a r i o u s g e o t e c m l i c a le n g i n e e r i n gp r o b l e m t l l i st l l e s i s ，b a s e do nc o m p r e h e n s i v ea i l a j y s i so fm e d o m e s t i c 龇l di n t e m a t i o n a lr o c ka r l d s o i lm m l e i l i c a lm e t h o d sa tp r e s e n t ，h a ss t u d i e do n i n t e r f a c ea n a l y s i sm e m o d 讹s t a t i cb e h a v i o r ，r e l e v a n ta l g o r i t l m la n dc o m p u t e rp r o c e d u r e ， 肌di t sa p p l i c a t i o np r o b l e m si ne n g i n e e r i n g f i r s t ，饰s t l l e s i sd i s c u s s e dt h ec u r r e n tr e s e a r c ha i l d d e v e l o p m e mo fr e l a x a t i o n n u m e r i c 甜m e t h o d sf o rr o c ka i l d s o i le n g i n e e r i n gb o t hi n l 觚da 1 1 da b r o a d l a t e r ，i nt e n i l so f m e c h a l l i c s ，b a s i cp r 王n c i p l e so ff o r c ed i s t r i b u t i o nm e t h o d 锄di n t e r f a c ea n a l y s i sm e t l l o db y s t a t i cr e l a x a t i o na r ei n t r o d u c e da n da i l a l y z e d a ne x 锄p l ei sc 撕e do u tb yf o r c ed i s t r i b u t i o n m e t h o d c o m p a r e d 谢t hf i i l i t ee l e m e mm e t h o d ，t h cr e s u i t ss h o wt h ee x a c t n e s sa n df e a s i b i l i t y m e a n w h i l e ，t 1 1 ed e f o m l a t i o np r o b l e mo fi n t e m 屺ew h e nac a n t i l e v e rb e 锄i sb o m i i lt e n s i o n a n dp r e s s e di sa n a l y s i sb yp r o p o s e dm e t h o d t h ea p p l i c a t i o nc o m p u t e rp r o g 舢lf o ri n t e r f a c e sa 1 1 a i y s i sb a s e do ns t a t i cr e l a x a t i o n s d e v e l o p e d ，i n c l u d i n gi t sr e l e v a n tp o s tp r o c e s s o ra n dp r e p r o c e s s o rp r o g r a m s t h ep r e p r o c e s s o rp r o g r ；h ni n c l u d e sd r a w i n ga n db u i l d i n gt h en u m e r i c a lm o d e l ，n u i i l e r i c a lg r i dm e s h ， i 印嘣洫gp a r a i i l e t e r sa i l da p p l y i n gr e s t r a i l l s as i m p l i f i e dt r i a n g u l a rm e s hg e n e r a t i n gm e t l l o d f o ra r b i t r a r yp l a l ls h a p ei sp u tf o r w a r d ，w h i c hc a i is 缸i s 母c e r t a i nr e s t r a i n sc o n d i t i o n ，a i l di t s c o m s p o n d i i l ga l g o r i t h n li ss t u d i e d t h er n e s hg e n e r a t i n ga l g o r i 伽 i lc a i ld i v i d et h ez o n ei n t o t r i a n g l em e s ha c c o r d i n gt 0m es i z ea n db o r d e rb ep r 0 v i d e d ，卸t o m a t i c a l l yr e c 0 r dn o d a l c o o r d i n a t e sa n du n i ts e f i a lm n b e r ，e x p o r t 锄ds h o wt l l er e s u l t s o m ee x 锄p l e sa r ep r e s e n t e d i 1 1t l l e s i s t h ep o s tp r o c e s s o rc a ns h o wt h eu n “a i l dn o d 甜d i s p l a c e m e n t s ，s 缸e s s ，a n dd y n 锄i c p r o c e d u r eo fc a l c u l a t i n g t h et l l e s i sa l s 0p r e s e n t st 1 1 e c o u r s et os h o wt 1 1 ep r o 铲a m ，t i l e a l g o r i t h m sa 1 1 dn o w c h a n so ft h ek e ys t 印s t h ed e v e l o p e dc o m p u t e rp r o g r 锄i su s e df o ra i l a l y z i n gm ep r o b l e mo fj a c k e dd o w n p i l e a x i a lf i o r c ea n ds k i n 衔c t i o no f t h ep i l ea r eo b t a i n e da i l dt l l e i rr e g u l 捌t i e sa r es t u d i e d b y 咖d y 试gt h er e l a t i o nb e t w e e nt l l el o a da p p l i e dt ot h et o po fp i l ea n ds e t t l e m e n t ，t h el o a d 缸觚s f e rm e c h a i l i s mo fi a c k e dd o 、v np i l ei sd i s c u s s e d k e y w o r d s ：s t a t i cr e i a x a t i o n ；i n t e i f a c e ：f b r c ed i s t r i b u t i o nm e t h o d ；t r i a n g l em e s h ；j a c k e d d o w np i l e i i 长沙理工大学 学位论文原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的 研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或 集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集体，均 已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。 作者签名：寺慈珂王 日期：2 。产5 。月2 ，厂日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保 留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借 阅。本人授权长沙理工大学可以将本学位论文的全部或部分内容编入有关数据库 进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 本学位论文属于 1 、保密口，在年解密后适用本授权书。 2 、不保密团。 ( 请在以上相应方框内打“ ) 储签名：和陲 刷醛轹同，谚詹专 日期：7 年厂月矿日 醐谚r 町日 第一章绪论 j l - - 一 1 1 目l j吾 岩土材料具有高度的非连续性、非均匀性和各向异性等鲜明的地质特征，在 力学性质上表现出强烈的非线性和非弹性。许多国内学者f 1 ，2 】均对岩土力学研究现 状、最新进展以及岩土工程技术的发展现状作过详细综述。由于经典力学对解析 解的设定非常严格，使其用于岩土工程的求解极为有限。随着计算机科学的发展， 依赖于计算机科学的数值计算在岩土工程研究领域愈来愈得到重视，数值计算和 分析己成为解决各类岩土工程问题的重要和有力的手段，数值计算也向计算规模 大型化、问题复杂化发展，如采用更为复杂的非线性本构模型、岩土介质细观计 算分析以及并行计算等等。然而，随着计算规模大型化、问题复杂化的增长，也 使计算分析对工程适应能力表现出了其局限性，如计算速度和计算耗时的问题； 数据信息量大、难检查； 问题的复杂化( 如复杂的本构方程) 出现的计算结果收敛 问题和计算结果是否唯一，以及如何用试验验证等问题。经过大量的岩土工程数 值分析，人们开始认识到在计算分析时单纯地追求问题或本构模型的复杂，并非 能得到与实际更相符合的结果，而简单实用、物理概念一目了然的方法更容易被 工程师所接受，因而使得本文对静力松弛数值分析方法的研究工作得以进行。 1 2 国内外研究状况 解决岩土力学问题的方法主要有实验方法、理论分析方法和数值模拟方法三 大类。其中数值模拟是解决岩土工程问题较有效的手段，已被学术界和工程界广 泛接受。目前，岩体工程课题中，计算技术已遥遥领先，各种数值分析方法日异 推新且方兴未艾。这些技术与方法的发展，为研究岩体工程问题提供了强有力的 工具。我国学者【3 ，4 】都曾对岩土力学数值分析方法与最新进展进行过系统地评述。 总的来说，主要分为连续变形数值分析方法和非连续变形数值分析方法两大类， 以及各种方法的耦合，数值计算从确定性向非确定性的过渡【5 】等。 1 3 岩土工程问题的数值分析方法 1 3 1 连续变形数值分析方法 连续变形数值分析方法主要包括有限单元法、有限差分法、边界元法、无限 元法等，其中以有限元法应用最为广泛，此类方法主要针对岩土介质的连续小变 形和小位移特性进行分析。有限元法在连续性分析方面取得了很大的成功，但在 解决前处理问题、应力与应变解答不连续问题和进行任意路径开裂计算等方面还 存在着一些局限。为了充分考虑岩土介质的非连续性、非均匀性和多相性等物理 特性，必须对这些连续变形分析方法，特别是有限单元法进行深入的改进和发展。 以连续介质变形分析为目的快速拉格朗日分析法( f l a c ) ，可以同时考虑岩土体的 材料非线性和几何非线性，并能跟踪物体变形的全过程，适于分析岩土力学中的 大变形问题。这种方法避免了有限元法进行大型矩阵的复杂计算，但时间步长的 选择成了一个新的突出问题，时间步长过大会导致解答的不稳定，时间步长太小 则会使计算时间过长。 1 3 1 1 、有限单元法 有限单元法是随着电子计算机的发展而迅速发展起来的一种现代计算方法。 计算的一般思路为：首先对物体进行离散化，将结构离散为有各种单元组成的计 算模型，这一步称作单元划分。离散后单元与单元之间利用单元的节点相互连接 起来；单元节点的设置和单元性质、数目等应视问题的性质、描述变形形态的需 要和计算进度而定( 一般情况单元划分越细则描述变形情况越精确，即越接近实际 变形，但计算量越大) 。所以有限元中分析计算所获得的结果只是近似的。如果划 分的单元数目非常多而且又合理，则所获得的结果就与实际情况相符合。然后分 析单元特性，选择位移模式。在有限单元法中，选择节点位移作为基本未知量时 称为位移法；选择节点力作为基本未知量时称为力法；取一部分节点力和一部分 节点位移作为基本未知量时称为混合法。位移法易于实现计算自动化，所以，在 有限单元法中位移法应用范围最广。当采用位移法时，物体或结构物离散化之后， 就可以把单元总的一些物理量如位移，应变和应力等由节点位移来表示。这时可 以对单元中的位移的分布采用一些能逼近原函数的近似函数予以描述。通常，在 有限单元法中，我们就将位移表示为坐标变量的简单函数。这种函数称为位移模 式或位移函数。根据单元的材料性质、形状、尺寸、节点数目、位置及其含义等， 找出单元节点力和节点位移的关系式，这是单元分析中的关键一步。此时需要应 用弹性力学中的几何方程和物理方程来建立力和位移的方程式，从而导出单元刚 度矩阵，这是有限单元法的基本步骤之一。假定力是通过节点从一个单元传递到 另一个单元中去的。因而，这种作用在单元边界上的表面力、体积力和集中力都 需要等效的移到节点上去，也就是用等效的节点力来代替所有作用在单元上的力。 最后对单元组集求解未知节点位移，利用结构力的平衡条件和边界条件把各个单 元按原来的结构重新连接起来，形成整体的有限元方程，解有限元方程得出位移。 应用有限元法能较好地反映复杂的岩体介质，可以解弹性、塑性和粘性问题， 可以解动载、静载问题，可以解流体问题，也便于施工程序的模拟。解弹塑性问 2 题时采用增量理论，是交剐度法。用迭代法求勰时，可以瘸常透度法或交刚度法， 从道理讲后者较前者收敛快，由于解题中源程序处理的差异，常刚度法省机时且 稳定性好。 需要指出的是，对岩体工程至关重要的结构面模拟问题节理元( 固体或流 体) ，在理论上尚欠成熟，无厚度元的质量为零，使动力问题求鳃困难。虽然可以 采用节理单元模拟不连续面，但问题之一是节理单元的数量不能太多，之二是节 理单元本质上仍然是一种有限单元，单元开裂后的力学行为无法模拟下去，因此 节理单元有限元法所麓模拟的岩体的罩# 连续大变形是十分有限的。计算中多组结 构面相交时，剖分处理要非常谨慎，它可能使岩体相互嵌入或分离，甚至求解失 败骖，钍。 1 3 1 2 、边界元法 边界元法是和有限元法并行发展的另一种数值方法。它以定义在边界上的边 界积分方程为控制方程，通过对边界分元差值离散，化为代数方程组求解。它与 基于偏微分方程的区域解法褶毖，由予降低了赫题昭维数，焉显著降低了富由度 数，边界的离散也比区域的离散方便的多，可用较简单的单元准确地模拟边界形 状，最终得到阶数较低的线性代数方程组。又囊于它剩孀微分算子的解析的基本 解作为边界积分方程的核函数，因而具有解析与数值相结合的特点，通常具有较 高的精度。特别是对于边界变量变化梯度较大的闻题，如应力集中问题，或边界 变量出现奇异性的裂纹问题，边界元法被公认为比有限元法更加精确高效。由于 边界元法所利用的微分算子基本能自动满足无限远处的条件，因丽特别便于处理 无限域以及半无限域闯题。单元划分与有限元寄所不同的是，边赛元法只需在边 界上进行离散化。因为_ i 程中大多力学问题都可以归结为l a p l a c e 或n a v i e r 方程， 这强种编微分方程的特点在予，只要能找到满足给定方程又符合边乔条件的函数， 这个函数便是所求的唯一解。把求解区域的边界离散后，将求函数解简化为求单 元节点上的函数值，降低了方程组的维数，使得工作量比有限元大为减少，这一 点在无限或半无限空间中计算尤为明显。但是，一般边界元法得到的线性代数方 程组的关系矩阵是一个不对称满阵，不便应用在墨前有限元法中已成熟的对稀疏 对称矩阵的线性代数方程组的一系列解法，也不适宜求解比表觚积大的计算域。 对多种介质的区界相当于两条边界，因此对计算域内介质分区较多( 如含有多种岩 石成分及较多的断层、节理等结构面，不同的风化分带) 的情况不宜用边界元方法 求解。对于复杂的非线性问题以及模拟分步开挖及施工过程等，不如有限元法方 便。经过4 0 多年的发展，在罄体力学领域，边界元法已经被定为有限元法的一种 重要补充【8 ，叭。 1 3 1 3 、有限差分法 3 差分法是微分方程的一种近似数值解法( 1 0 】。它不是去寻求函数式的解答，而 是去求出函数在一些网格节点上的数值。具体地讲，差分法就是把微分用有限差 分代替，把倒数用有限差商代替，从而把基本方程和边界条件( 一般均为微分方程) 近似地改用差分方程( 代数方程) 来表示，把求解微分方程的问题改换成为求解代 数方程的问题【1 1 1 。这一方法简便直观，且在流体动力学、热交换和固体力学领域 内常被应用。传统差分法的差分网格为规则均匀，由于其在处理岩体内的节理、 复杂边界和非均质材料时存有缺陷，故不适合模拟实际的岩石力学问题。 1 3 1 4 、无限元法 无限单元法是为了克服有限单元在求解无界域问题时的不足才提出的。然而， 在有限单元尺寸和数目之“有限”与无界域之“无限”间必须进行折中和妥协，其中 最自然的办法就是保留有限元法的所有特征，丽对无界域进行近似处理，即人为 截取“足够大”的区域进行几何上的网格剖分，同时，在“人为”边界上施加相应的 近似约束边界条件。这种方法虽然不需引入新概念，但在理论上会带来较大的误 差：区域较大能减少理论误差，但数值计算规模将以数倍增加。实际上，处理无 界域问题还可以采用这样一种方式，引入一种几何上无限大的“有限”单元，再对 其在物理上进行界定，这就是无限元法思想的由来。无限元( i n n n i t ee l e m e n t ) 这一 术语1 9 7 7 年由b e t t e s s 和z i e n k i e w i c e 【1 2 】第一次使用，在概念上它是有限元的延伸， 是一种几何上趋于无限的单元，即它所占的区域是无限的；又由于无限元必须反 映近场的边界特征或与模拟近场的有限元结合，它实际上只在一个方向趋于无限， 因而又被称为半无限元。由于有限元的概念涵盖所有占非无穷小区域的单元，广 义的讲，无限元仍然属于有限元的范畴【1 3 14 1 。总之，无限元为克服有限元在解决 无界域问题时而提出，常常与常规有限元同时用来解决更复杂的无界问题，是对 有限元方法的“协调”与生俱在，比边界元等其它求解无界域问题的数值方法更有 优势【1 5 】。 1 3 2 非连续变形数值分析方法 在岩土工程中，常会遇到各种非连续变形力学问题，这是由工程岩土的特征 决定的。岩土是在漫长的地质历史发展过程中形成的地质体的一部分，它的成因 和构造复杂，岩性多样。为了综合研究岩土的工程特性，近年来人们提出了岩土 结构这一概念，认为岩土是有结构的，岩土的力学作用主要受岩土结构面控制， 岩土的变形与破坏一般都发生在结构面上。所谓结构面是在地质发展历史中，岩 土内形成具有一定方向、一定规模、一定形态和特性的面、缝隙以及带状的地质 界面，如层面、片理、节理、软弱夹层，以及断层破碎带等。岩土就是由结构面 及所包围的结构体共同组成的。 正是由于岩土的这种结构特点，在建筑物与岩土的交界面上、结构的分缝面 4 上等均具有普遍的滑移或脱开等非连续变形特性，使应力或位移不连续，而这些 交界面上的非连续力学行为极大地影响着结构的受力和变形，往往成为问题分析 的关键。随着以有限元为代表的连续变形力学离散模型的日趋完善和广泛运用， 学术界和工程界已不再满足在连续变形力学框架下得到上述问题的解答。因此， 近年来，先后出现了不少的理论和模型来解决这些非连续变形力学问题【1 6 l ，迄今 为止已有数十种方法。其中应用较广的有：块体理论；非连续变形分析方法( d d a ) ； 离散元法；快速拉格朗日分析( f l a c ) 法等。 1 3 2 1 、块体理论 块体理论是7 0 年代由石根华提出，后来他在美国与r e g o o d m a n 教授合作 完善起来的一种岩体稳定性分析方法。块体理论的标志性论著为块体理论及其 在岩体工程中的应用( b l o c kt h e o r ya n di t sa p p l i c a t i o nt or o c ke n g i n e e r i n g ，19 8 5 ) 。 该书的出版具有划时代意义，使块体理论作为岩体工程稳定性分析的一种有效方 法得到了国际岩石力学界的普遍赞誉和推崇。 块体理论是以岩体完全被交错的结构面所切割为前提的，它适用于节理及块 状岩体。石氏块体理论的基本假定为： ( 1 ) 结构面为平面并贯穿所研究的岩体区域； ( 2 ) 结构体为刚体，即不考虑块体的自身变形和结构面的压缩变形； ( 3 ) 岩体的失稳是岩体在各种荷载作用下沿着结构面产生的剪切滑移。 块体理论有两个基本定理，即有限性定理和可动性定理，这两个定理已经给 予了严格的数学证明，是块体理论的核心【1 7 】。 1 3 2 2 、非连续变形分析法( d d a ) 不连续变形分析方法d d a ( d i s c o n t i n o u sd e f o r m a t i o na n a l y s i s ) 是石根华提出 的一种基于岩土介质非连续性提出的分析块体系统运动和变形的新数值分析方 法。d d a 理论的基本内容是：以自然存在的岩土被节理面或断层面等结构面切割 形成不同的块体单元，块体的运动和变形由刚体位移、转动、正应变和剪应变组 成；非连续变形分析以各个块体的位移为未知量，通过块体的接触和几何约束形 成一个块体系统，块体单元受非连续面的控制，在块体运动的过程中单元之间可 以接触也可以分离，在块体运动过程中，满足块体之间不侵入和不承受拉伸力的 条件；总体平衡方程由系统的最小势能原理求得，求解方程组就可得到当前时步 的位移场、应力场、应变场和块体间的作用力，反复形成和求解总体平衡方程式， 即可得到多个时步后块体的变形、位移和应力情况，也可求得块体系统最终达到 平衡时应力场、位移场以及运动过程中各块体的相对位置和接触关系【1 8 五0 1 。在已 知块体系统中知道每个块体的几何条件、边界条件及力学条件后。就可以采用非 连续变形分析理论计算块体系统中每个块体的位移、应力和应变，从而确定块体 间相对移动和滑动。因此，d d a 法可以模拟出岩石块体的移动、转动、张开、闭 合等全部过程。据此，可以判断出岩土的破坏程度、破坏范围，从而对岩土整体 和局部的稳定性作出正确评价。 1 3 2 3 、离散单元法 离散单元法是由c u n d a l l 于2 0 世纪7 0 年代提出的。离散单元法是针对节理 岩体提出的一种适用于模拟岩土体大位移的数值方法。它是在块体准刚性假设的 前提下，以牛顿第二定律为理论基础建立起来的，以微小运动状态的求解模拟岩 体的大位移。其基本思想是将所研究的物体看作被节理在内部分割的离散状块体 单元，块体间通过边界接触力相互联系，边界接触力用物理方程来描述。当所研 究的块体在力系作用下或边界条件发生变化时，单元就在牛顿第二运动定律下发 生平移或转动，允许调整各个块体单元的接触关系，接触关系的调整产生新的接 触力，在该力的作用下产生新的的位移。如此循环：在不平衡力作用下，块体产 生运动。不断调整位移、调整接触力，最终块体可能达到新的平衡状态，也可能 一直运动下去【2 1 1 。 由于岩块或颗粒组合体被模拟成通过角或边的接触面相互产生作用，块体之 间边界的相互作用可以体现其不连续性和节理的特性，使用显式的时步迭代算法， 对于岩基和节理块体，允许有大的位移，转动和使用各种非线性本构关系。因此 可以求解单元块的失稳过程。如果单元的分割与结构面的强度指标选取能满足工 程精度要求，则用这一方法分析的岩体失稳过程和形态是基本符合实际情况的。 但是，岩石介质种类繁多，性质极为复杂。计算用主要参数( 如阻尼参数、刚 度系数等) 的选取还没有统一和完善的确定方法；计算时步的选择原则是在数学方 程趋于收敛的条件下确定的，与实际中的“时间”概念如何联系起来，并合理地考 虑时间，这是需要进一步研究的问题。刚性单元在模拟坚硬岩体之间的作用方面 比较适用，对于软硬悬殊的岩体的变形模拟显得过于简单，计算结果与实际差距 太大。目前迭代计算次数太多，运算时间冗长，不同的时步单位及不同的单元划 分可能得出相差较大的结果。 目前离散单元法大多采用动态松驰法【2 2 1 。动态松驰法把非线性静力学问题转 化为动力学问题进行求解，用显式中心差分法近似对运动方程进行积分计算，并 假设块体在运动时将动能转化为热能耗散掉，在计算中引入人工粘性阻尼，使系 统达到平衡、运动趋于稳定。经过近3 0 年的发展，离散单元法已成为模拟岩土体 非连续大变形行为的十分有效的数值方法之一。目前，二维离散单元法已趋于成 熟，一些商业化软件已经进入实用阶段，用于解决_ 些工程实际问题。自1 9 8 8 年c u n d a l l 发表三维离散元的具体算法以来，三维离散单元法也得到了迅速发展。 国内刘连峰【2 3 1 和焦玉勇【2 4 】两位学者的工作较有代表性。但是由于三维问题算法的 6 复杂性，特别是块体间接触判断算法的复杂性，使得从二维离散单元法向三维离 散单元法的过渡花费了近2 0 年，而且目前仍处于研究阶段。 1 3 2 4 、快速拉格朗日分析( f l a c ) 法 快速拉格朗日法是一种基于显式有限差分法的数值分析方法，它可以模拟岩 土或其他材料的力学行为。快速拉格朗日分析将计算区域划分为若干个四面体单 元，每个单元在给定的边界条件下遵循指定的线性或非线性本构关系，如果单元 应力使得材料屈服或产生塑性流动，则单元网格可以随着材料的变形而变形，这 就是所谓的拉格朗日算法，这种算法非常适合于模拟大变形问题。快速拉格朗日 分析采用了显式有限差分格式来求解场的控制微分方程，并应用了混合单元离散 模型，可以准确地模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的 弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点【2 5 ，2 6 1 。但是， 同有限元一样，计算边界、单元网格划分带有很大的随意性。 1 3 2 5 、数值流形法 数值流形法是石根华应用流形的覆盖技术建立的一种把有限元法、非连续变 形分析法和解析法包含在内的全新的统一计算方法。它包融了有限元与d d a 法， 是岩石力学中一种新的较通用的数值分析方法。由于流形法可在统一的理论框架 下处理连续与非连续性变形问题而引起了许多学者的兴趣，成为目前计算岩土力 学的热门课题。 数值流形法的优点主要表现在它具有相对完善的非连续变形处理功能，可以 在统一的数学理论框架下同时处理连续问题与非连续问题。数值流形法的重大意 义是把连续和非连续变形的力学问题统一起来，指出了一个新的方向，昭示了现 代数学手段和岩土数值模拟结合的必要性。数值流形法较有限单元法更适于进行 开裂模拟，但由于受网格连接与单元划分的限制，流形法在开裂计算上仍存在一 定的困难【2 7 1 。 1 3 2 6 无单元类方法 无单元类方法因具有无须单元网格划分、前后处理简单、较传统有限单元法 更适合断裂问题的计算分析等优点而受到学术界的广泛关注。无单元类方法在进 行裂纹扩展模拟时不再存在传统有限单元法的重新剖分网格的困难，而仅仅在裂 尖局部区域内布置节点，大大简化了前处理过程，在剖分策略上，无单元法较传 统的有限单元法更适于断裂问题的计算分析。 对于以无单元伽辽金法( e f g m ) 为代表的无单元类方法，现有研究成果大多集 中在裂纹扩展的模拟方面，在将无单元方法应用于岩土力学的数值计算方面仍存 在很多困难，如非连续材料插值函数的构造、摩擦接触问题和多体相互作用问题 7 的处理等还需作进一步研究。 。 1 3 3 各种方法的耦合 各种数值方法的耦合应用，能够扬长避短，改善精度，提高计算效率。 1 3 3 1 、有限单元法与边界单元法的耦合 利用边界元法处理无限域边界的优点以及有限元法求解非线性问题的灵活 性，可以克服有限元法的“边界效应”及边界元法处理多介质的困难。然而有限元 与边界元结合程序庞大，则是使用者熟悉程序的难题。 1 3 3 2 、离散单元法与有限单元法的耦合 离散单元法与有限单元法的耦合计算比较容易。只要使交界面上的有限单元 的节点与离散单元的角点重合，并保证它们的位移和力连续，就可通过节点力和 位移的相互传递将离散单元与有限单元耦合起来。 1 3 3 2 、离散单元与边界单元法的耦合 对于无限域问题，在开挖边界周围一定范围内用离散单元法离散，而在此范 围以外，则用边界单元法离散。具体做法是：首先使边界单元的节点在交界面上 与离散单元的角重合，保证几何相容。由于两个块体的角对应于边界单元的一个 节点，边界单元的节点位移可以取这两个相邻块体角位移的平均值，以保证运动 相容。其次，将边界单元法计算出的节点力反向施加到与其相连的两个块体上。 由于牵引力沿边界单元线性分布，因此节点力可按边界单元的长度大小比例分配 到相应的块体形心上，以保证交界面上力的平衡。 1 3 4 数值计算从确定性向非确定性过渡 1 3 4 1 、随机有限元法 对岩体工程进行有限元分析时出现的矛盾，一方面，随着有限元理论和方法 的发展，本构关系研究的深化，一批高精度单元的提出，计算精度越来越高；另 一方面，在介质的材料性能参数，边界条件及载荷方面又存在有很大的统计随机 性。显然，有限元计算的高精度被随机性的影响大大降低了。 随机有限元分析法是确定性有限元技术与概率统计方法相结合的一种数值方 法。随机有限元法多采用以摄动法为基础的一次二阶矩分析，求出响应( 应力、位 移等) 的均值、方差、协方差，并据此求出岩体结构内的局部破坏概率。然后利用 动态规划法等寻求最可能的滑移面，得到岩体结构的整体可靠度。 1 3 4 2 、模糊有限元法 据模糊数学的观点，岩体的不确定性和非精确性叫做模糊性。这种模糊性： 8 表现在岩体力学性质的不确定性和不精确性，如岩体的弹性性质，不连续特征 和整体强度等显示出很大范围的模糊性，而且也不能用试验的方法对整个岩体进 行试验并精确地确定其性质和强度；岩体工程的模糊性还表现在边界条件上， 亦即边界位移条件和约束条件本身是不确定的；荷载的不确定性，如很难确定 岩体中的原岩应力和残余应力等。 模糊有限单元法能让使用者输人岩体力学参数，边界位移条件和边界荷载条 件的模糊范围得到输出所求量的模糊范围，并能对计算结果进行模糊分析，例如 资格函数分析，敏感性分析，局部安全系数分析，贡献权系数分析以及有限元计 算结果的可靠性分析等。但本法目前主要用于平面问题计算。数值计算的非确定 性方法除了以上介绍的两种之外，还有摄动边界元，模糊边界元法等。 1 4 节理岩体的界面及其数值模拟 土木、水利等工程中所接触到的岩体总是不连续的，而含有节理的岩体是最常 见的不连续体节理的存在将减小岩体的强度，增加岩体的变形，致使岩体的强度 和变形表现为不均匀性和各向异性，还能使岩体中的岩块常常发生开裂、滑移和运 动【2 引。因此，针对不连续体的数值分析界面处理是一个很重要的问题。人们早期 用常规的有限单元法、有限差分法等模拟岩石材料。有许多程序包含了节理单元 或“滑动线 1 2 9 30 1 。后来人们引入“离散单元”的概念( c u n d a l l ，1 9 7 1 ) ，成功地模拟 了界面之间的各种复杂运动。当然，基于“离散单元”模型的模拟是适合于岩体富 含节理的情况。总的说来，岩体的不连续面的数值模拟方法可分成两类：一类是 节理单元模型方法；另一类是离散模型方法。 ( 一) 、应用于节理岩体的有限元数值仿真 2 0 世纪6 0 年代，将有限元应用于节理岩体的数值模拟中最有代表性的是r e g 0 0 d m a n 提出的节理单元，也称接触面单元，用构建零厚度的矩形单元模拟岩石 内部的节理。由于它是建立在连续性假定的基础上，故不能用于大变形、滑移及 完全脱开的情形，同时，由于节理单元的零厚度，故容易出现数值病态。 o c z i e n k i e w i c z 等在接触面单元的基础上建议采用六节点节理单元，并赋予单元 一个小的厚度【3 1 1 ，但这同样不能消除数值病态。j g h a b o u s s i 等【3 2 1 将节理两个相 对面的相对位移作为独立的未知量，基于塑性理论提出了节理单元法，即在性质 差异悬殊的两相邻单元界面上，设置法向与切向弹簧，一旦法向弹簧受拉则视单 元界面发生分离，若切向相互作用超过单元界面的极限强度则认为界面产生相对 滑移。其后，c s d e s a i 等【3 3 】又提出了薄层单元法。即假定界面是由特殊介质组成 的、厚度很薄( 约为周边单元平均宽度的1 1 0 0 1 1 0 ) 的实体单元，而这种单元的 本构性质取决于周边介质力学性态，由试验直接测得。以上这些方法均是将岩体 视为准连续介质，以连续分析为主，并对个别具有控制作用的宏观非连续面的变 9 形与破坏等力学效应给予重点分析，但由于一些自身的缺陷而不能处理复杂的非 连续变形问题，如只能计算原生的非连续节理面而对次生的非连续节理面则需重 新划分网格，处理也较为困难；节理面单元不宜设置太多，否则整体劲度矩阵容 易出现病态；选取节理面的法向刚度较为困难。为克服这些缺陷，t b e l v s c h k o 【3 4 】 等提出加密有限元用以分析节理岩体的扩展，即用任意的不连续函数及其导数直 接表示目标函数，而不必让有限元网格去适应节理，也无需通过有限元网格的再 划分来模拟节理的扩展。用仅包含节点函数值的位移函数来表达岩石节理面，并 在节理面及节理端点分别引入跳跃函数和缝端函数，以提高有限元分析节理岩体 的能力。结合有限元在处理材料的非均质、非线性上的优势，使“e n r i c h e df e m 可用于节理岩体的非线性分析。 ( 二) 、离散模型界面分析 由于岩石力学问题所遇到的断层、节理等一些特殊的结构，很多情况下连续 体力学求解方法明显不能适应，离散单元计算方法在这种问题中表现出了更大的 优势和灵活性。离散单元法的理论基础是结合本构关系的牛顿第二定律，可以采 用动态松弛法或静态松弛法求解。目前，大多采用动态松弛法将非线性静力问题 转化为动力问题来求解【35 1 。在“离散单元”模型的描述上，均有一些共同点：即单 元是离散的，理论上单元可产生任何形式的运动，不连续面就是分割各单元的分 界线，不连续面的力学行为就是离散块体单元外边界在其运动式变形域内的相互 的一些力学行为。从“离散单元”的建模意图和其所具有的块体力学机制看到，其 对岩体的不连续面模拟是很充分的，已得到普遍承认。值得一提的是，作为连续 介质分析方法f l a c ，其对不连续面的模拟就是按照离散单元的处理模式进行的。 上述第一类方法是将不连续面通过节理单元这种特殊的方式来进行不连续面 的模拟。由于节理单元掺入了一些人为的假定，且对系统方程组很敏感，很难模 拟不连续面的大的错动和张裂等变形形式，因此只能在与连续岩块变形的同量级 上反映不连续面的变形；第二类方法能模拟不连续面的这种错综复杂的变形形式， 然而在岩土工程中真实的这种离散的区域是很难遇到的，其单元模型的划分也远 比连续介质模型困难，且对岩体这类原始数据不充足的问题，很难真实地反映其 物理力学规律，这是此类方法计算结果与真实情况出入较大的主要原因，也是第 二类方法至今还没有有限元方法应用广泛的原因之一。 对不连续面( 严格地说是含不连续面的岩体) 的模拟，f l a c 方法为我们开辟了 一条新途径。f l a c 将连续的分析方法和离散单元方法的界面处理方式两者相结 合，使不连续面的模拟有了新的突破。 静力松弛拉格朗日分析方法，同样也将结合连续与离散两类方法的优点，模 拟真实的岩体断续变形。 l o 1 5 松弛数值方法综述 所谓松弛法，指的是块体位移的求解方法。松弛数值方法包括动态松弛离散 单元法【3 6 1 、静态松弛离散单元法【3 7 1 、静态同步松弛离散单元法【3 8 1 和求解连续介 质模型的静态松弛拉格朗日方法等【3 9 l 。 c u n d a l l ( 19 7 1 ) 提出的离散单元法是基于动态松弛解法的，它采用牛顿第二 定律或的差分格式求解块体的位移，并用时步积分累积求算块体的大位移，因而 动态松弛法存在着动能的耗散和时间步长的选取问题，一般认为，时步和阻尼的 确定是动态松弛解法中最为困难的问题之一。针对这一弊端，s t e wa r t ( 1 9 8 2 ) 提 出了静态松弛离散单元法。静态松弛法根据块体产生不平衡力后达到再平衡时的 力位移关系建立平衡方程组，通过求解方程组得到块体的位移。由于采用隐式格 式和静力解法，静态松弛法完全避开了时步和阻尼【4 们，并且计算效率大大提高。 静态松弛法与动态松弛法的根本区别就在于块体位移的求解方法不同。众所 周知，动态松弛法的出发点是牛顿第二定律，由加速度一速度_ 位移的求解格式 使得阻尼系数和时间因素的引入成为不可避免，继而导致系统振荡控制、加快收 敛速度、防止块体相互嵌入等一系列问题。与此相对照，静态松弛法则采取了较 为简单的处理方法：根据产生不平衡力的块体达到再平衡时的力位移关系建立平 衡方程组，由方程组隐式求解块体的位移。除块体位移的求解外，两种松弛法的 其他方面几乎是相同的。 静态同步松弛法的求解，实际上是一个块体系统在平衡和失衡两种状态之间 反复交替、最终达到平衡的过程。系统在外力( 如重力、边界力、地震力、水压力 等) 和边界条件等因素的作用下，块体单元受到不平衡力的作用，要达到平衡状 态，必须以失衡块体发生位移来消除不平衡力。而位移的产生必然导致块体间接 触条件的变化，接触力发生变化，块体又处于新的不平衡状态，其空间位置又相 应发生变化，如此循环往复。所谓同步松弛，是为了避免求解顺序不同带来的误 差，使系统块体同时运动，这也符合岩体变形的实际。静态同步松弛法的计算过 程包括数次迭代计算，每次迭代的位移被记录在单元的角点坐标上，经过多次迭 代，块体系统的最终状态就很容易描述了。因此，对岩体大变形的模拟实际上是 数次小位移累积的结果。 快速拉格朗日分析方法( f l a c 【4 1 1 ，该方法其实质是将逐块离散元动态松弛法 引伸到连续介质问题中) 。f l a c 方法是属于动态松弛方法，被描绘为进行“拉格 朗日分析”的“显式的有限差分程序 【4 2 1 。这些计算特点决定了其与逐块动态松 弛的离散单元法有很多计算上的共同点，如时步参数的选定、阻尼项的选取等。 f l a c 又可称动态松弛( d y n a m i cr e l a x a t i o n ) 的快速拉格朗日分析方法。 本文所述的连续介质的静态松弛( s t a t i cr e l a x a t i o n ) 快速拉格朗日分析方法【4 引， 或s f l a c 。是与f l a c 方法平行的一种连续介质的松弛计算分析方法，与静态 松弛的离散单元法类似，通过建立差分网格节点( 离散元中是“块体”) 的逐点平衡 方程，通过逐点解平衡方程形成一个松弛的过程，最后求得系统的解。而静态松 弛界面单元法，是基于静态松弛拉格朗日分析方法而提出，所以又称静态松弛界 面单元法，该方法避免了出现人为的界面刚度系数砌、魅，界面上节点可以错动 和分开，结合了