[科普中国]-数值岩土力学

简介

岩石和土都是经历过变形、遭受过破坏的地质体, 受其成因、组成、结构、演化过程、生成年代及其所处的大地构造环境等诸多因素的影响, 岩土材料具有高度的非连续性、非均匀性和各向异性等鲜明的地质特征, 在力学性质上表现出强烈的非线性和非弹性。岩土工程是一门综合应用岩石力学、土力学、工程地质学与基础工程学等基本知识和手段解决实际工程建设中有关岩体与土体变形及稳定问题的学科, 它的主要任务是在极其复杂的地质条件、自然环境和人类活动中确保岩体和土体不会因强度不足或变形过大而使岩体和土体本身发生局部或整体的失稳破坏, 使与岩体和土体紧密相依的建筑物或构筑物不会失去正常稳定的条件或丧失工程功能。

近30 年来, 随着大型水利水电工程、土木工程、采矿工程、公路工程、铁路遂道工程等的迅速发展,作为岩土工程基础学科的岩土力学得到了长足的发展。世纪性和特大型的工程实践活动为岩土力学的发展赋予了巨大的动力, 同时也提出了许多复杂的岩土力学新课题。

解决岩土力学问题的方法主要有实验方法、理论分析方法和数值模拟方法三大类。这三类方法相辅相成, 互为补充。其中数值模拟是解决岩土工程问题的较有效的手段, 已被学术界和工程界广泛接受, 作为一种力学状态的分析工具, 它越来越多的应用于岩土体的稳定分析、岩土工程设计和岩土工程基本问题分析中。

近30 年来, 岩土力学数值分析方法得到了迅速发展, 出现了有限差分法(FDM)、有限单元法(FEM )、边界元法(BEM ) 、无限元法(IEM) 、拉格朗日元法(LEM) 、刚体弹簧模型或刚性有限元法(RBSM 或RFEM)、离散元法(DEM) 、非连续变形分析方法(DDA)、无单元法(EFM) 、流形元法(MM)及其耦合的数值计算方法和以数值模拟为主的渐进破坏模型等各种数值分析技术。这些数值分析方法可大致分为连续变形分析方法和非连续变形分析方法两大类, 它们分别将岩土体介质抽象为连续介质模型和离散体系模型两大类进行数值计算1。

连续变形数值分析方法这类方法主要包括有限差分法、有限元法、边界元法、无限元法等, 其中以有限元法应用最为广泛,这类方法主要分析岩土介质的连续小变形和小位移特性。有限元法在连续性分析方面取得了很大的成功, 但在解决前处理问题、应力与应变解答不连续问题和进行任意路径开裂计算等方面还存在着一些局限。为了充分考虑岩土介质的非连续性、非均匀性和多相性等物理特性, 必须对这些连续变形分析方法, 特别是有限单元法进行深入的改进和发展。

以连续介质大变形分析为目的的拉格朗日元法在实际工程中也得到了较好的应用。拉格朗日元法运用流体力学中跟踪质点运动的物质描述方法, 即拉格朗日拖带坐标系方法, 利用差分格式, 按显示时步积分方法进行迭代求解, 根据构形的变化不断更新坐标系, 以此模拟岩土介质的有限变形和大位移行为。基于拉格朗日元理论, 美国的ITASCA 咨询集团于20 世纪80 年代编写的专用程序FLAC 现已广泛应用于边坡、基础、坝体、隧道、地下采场和硐室等岩土工程分析中。拉格朗日元法可以同时考虑岩土体的材料非线性和几何非线性, 并能跟踪物体变形的全过程, 适于分析岩土力学中的大变形问题。这种方法避免了有限元法进行大型矩阵的复杂计算, 但时间步长的选择成了一个新的突出问题, 时步过大会导致解答的不稳定, 时步太小则会使计算时间过长2。

非连续变形数值分析方法界面单元有限元法为了推广有限单元法, 使其能够处理简单的非连续问题, 人们提出了各种能够反映非连续变形性质的简单力学模型和特殊界面单元, 并将其应用于单一性非连续界面力学行为的模拟中。这些模型主要包括联结单元(Linkage Element)、无厚度接触单元或节理单元(Joint Element)、薄层单元(Thin-LayerElemen)和接触-摩擦单元(Contact-Friction InterfaceElement)等。

界面单元有限元法将岩土介质视为准连续介质, 仍以连续分析为主, 但可以对个别具有控制作用的宏观非连续面的变形与破坏力学行为进行重点分析, 在工程中得到了广泛应用。同时, 该方法也存在诸多的局限性, 如只能对原生的非连续界面进行计算, 对次生的非连续界面无法进行处理、界面单元数目不能设置太多、界面弹簧刚度的选取较为困难等。因此, 该方法在处理复杂的非连续变形问题时仍显得无能为力。

刚体-弹簧模型或刚性有限元法刚体-弹簧模型将离散后的块体视为刚体, 块体之间用界面上的法向弹簧和切向弹簧相联结, 以块体形心处的刚性位移为基本变量, 用分片的刚体位移模式逼近实际整体位移场, 以块体间的联结弹簧反映结构内部的弹性, 并用界面应力表征结构内部的应力。这类模型主要包括刚体-弹簧模型(RBSM)、刚性有限元、分块刚体位移-界面应力元和块体-弹簧模型等。此类模型的界面特性均服从Coulomb 摩擦定律, 对体系的静力学约束条件考虑得比较充分, 在连续状态的应力分析方面可给出较高的计算精度, 对于临界状态能够估算出极限荷载,并可有效的用于少量块体界面间的摩擦接触分析。

但这类方法过分强调岩土体结构面的作用, 对结构体的变形没有给予足够重视。虽然这种方法可以用于原生界面的破坏分析, 但不能模拟实际岩土体的破坏发展过程, 也无法模拟破坏发展导致的次生界面的非连续变形行为和块体失稳后的运动过程。

基于块体理论的非连续变形分析方法基于块体理论的非连续变形分析方法主要以离散块体系统为研究对象, 针对岩土体的非连续与非均匀的特点, 将岩土体视为完全非连续介质, 对构成离散系统的各个子块的运动和变形进行数值分析。其中研究较多且具有代表性的方法有离散单元法、非连续变形分析方法和流形元方法等。

离散单元法(DEM)是由美国学者Cundall 等提出的。它是以被软弱结构面切割而成的离散块体为基本单元进行分析的, 在块体间的接触约束下运用牛顿运动定律描述各块体的运动过程。离散单元法可分为动态松驰法和静态松驰法两种。目前离散单元法大多采用动态松驰法。动态松驰法把非线性静力学问题转化为动力学问题进行求解, 用显式中心差分法近似对运动方程进行积分计算, 并假设块体在运动时将动能转化为热能耗散掉, 在计算中引入人工粘性阻尼, 使系统达到平衡、运动趋于稳定。经过近30 年的发展, 离散单元法已成为模拟岩土体非连续大变形行为的十分有效的数值方法之一。

目前, 二维离散单元法已趋于成熟, 一些商业化软件已经进入实用阶段, 用于解决一些工程实际问题。自1988 年Cundall 发表三维离散元的具体算法以来, 三维离散单元法也得到了迅速发展。国内刘连峰和焦玉勇两位学者的工作较有代表性。但是由于三维问题算法的复杂性, 特别是块体间接触判断算法的复杂性, 使得从二维离散单元法向三维离散单元法的过渡花费了近20 年的时间, 而且目前仍处于研究阶段。

非连续变形分析方法(DDA)是继离散单元法之后, 于20 世纪80 年代末发展起来的一种更新的用于模拟散体系统力学响应的数值分析方法。1988年, 石根华博士发表了博士学位论文“DiscontinuousDeformation Analysis :A New Numerical Model fo rthe Static and Dy namics of Block Systems” , 这标志着DDA 方法的诞生。该方法得到了国际认可, 受到了美国岩石力学权威学者Cook , Goodman , Desai 和Zaman 的极力推崇。为了推广DDA 方法, 从1996年起, 每两年召开一次国际非连续变形分析会议, 以进行这一领域研究成果的及时交流。这一方法在其诞生后的10 余年时间里, 发展速度非常快。与有限单元法相似, 非连续变形分析方法用隐式方程表达,所不同的是它引入了运动方程, 用最小势能原理把块体之间的接触问题和块体本身的变形问题统一到矩阵的求解中, 理论严密, 精度较高, 而且把静力和动力、正分析和反分析统一起来, 不仅可以计算破坏前的小位移, 也可以计算破坏后的大变形, 对滑动、崩塌、爆炸和贯入等问题也十分有效, 是一种不同于DEM 的新的数值计算方法。

由于DDA 方法在运动约束方面做得比较充分,理论上比较严密, 且静力和动力分析采用了统一的数值计算格式。因此, 在结构、岩体和土体的非连续大变形力学过程模拟方面发挥了较大潜力。

流形方法有限单元法将岩体理想化为完全连续介质,离散单元法等则将岩体抽象为完全非连续介质即离散介质, 这其实是走向了问题的两个极端。将连续变形与非连续变形统一起来的方法将更适于进行实际岩土体变形的分析。因此, 近年来另一种更新的数值方法--流形方法引起了国内外学者的广泛兴趣。

流形方法是石根华博士应用流形的覆盖技术建立的一种把有限元法、非连续变形分析方法和解析方法包含在内的全新的统一计算方法。由于流形方法可在统一的理论框架下处理连续与非连续性变形问题而引起了许多学者的兴趣, 成为目前计算岩土力学的热门课题。

流形方法的优点主要表现在它具有相对完善的非连续变形处理功能, 可以在统一的数学理论框架下同时处理连续问题与非连续问题。流形方法较有限单元法更适于进行开裂模拟, 但由于受网格连接与单元划分的限制, 流形方法在开裂计算上仍存在一定的困难。

无单元类方法对于岩体这种介于连续介质与完全非连续介质之间的材料, 沟通连续变形与完全离散体的非连续变形的桥梁是岩体中的地质非连续面的扩展、贯通等的破坏行为, 它是联接岩体变形过程中初始阶段终了阶段的一条纽带, 在整个岩体变形过程中, 这些非连续面的变形、扩展以及失稳对整个工程岩体结构的稳定是至关重要的。而现有的基于网格划分的方法均是基于单元离散思想的, 在分析裂纹扩展时, 都面临两方面的问题:

1)网格重新划分的问题;

2)裂纹扩展结果受到网格划分形式严重影响的问题。

无单元类方法因具有无须单元网格划分、前后处理简单、较传统有限单元法更适合断裂问题的计算分析等优点而受到学术界的广泛关注。无单元类方法在进行裂纹扩展模拟时不再存在传统有限单元法的重新剖分网格的困难, 而仅仅在裂尖局部区域内布置节点, 大大简化了前处理过程, 在剖分策略上, 无单元法较传统的有限单元法更适于断裂问题的计算分析。

对于以无单元伽辽金法(EFGM)为代表的无单元类方法, 现有研究成果大多集中在裂纹扩展的模拟方面, 在将无单元方法应用于岩土力学的数值计算方面仍存在很多困难, 如非连续材料插值函数的构造、摩擦接触问题和多体相互作用问题的处理等还需作进一步的研究。

渐进破坏模型近年来, 一些学者从细观力学的角度, 利用扫描电镜对岩石的细观时效损伤特性和损伤力学行为(包括岩石的细观组构、初始细观损伤特性、单调拉伸和压缩加载条件下的岩石强度、变形和细观破坏特性, 以及细观条件下岩石的时效损伤发展、蠕变特性和破坏性态等)进行了大量的细观试验分析, 这些研究成果得到了国内外学者的广泛关注。

岩石细观损伤演化实验表明, 岩石断裂破坏的实质是岩石在受力过程中微裂纹的萌生、扩展直至贯通的结果, 是岩石微观结构变形破坏累积的宏观反映。现有的各种数值计算方法, 仍然只能处理连续介质力学问题, 即使一些非连续介质力学方法(如DEM 和DDA)也难以考虑岩石从细观破裂到宏观破坏的过程。对于岩石介质而言, 模拟基于细观分析的岩石宏观变形、损伤直至破坏、失稳的过程以及研究地基、滑波、地震和岩爆等非稳定现象仍然是一个难题。

渐近破坏模型是通过对材料微结构的演化和局部变形行为的研究, 对材料在变形、损伤与破坏时所产生的本质现象和因果关系进行模拟的方法。与各种已有的数值计算方法不同, 它注重对材料的细观断裂机理与断裂规律的研究, 侧重对事物内在机理的模拟。渐进破坏模型是在对复合材料和岩石等脆性材料的破坏研究中提出的。这类方法一般包括两个方面, 即应力分析和破坏分析。应力分析可以采用解析方法和有限元法, 破坏分析是根据一定的破坏准则检查材料结构中是否有单元破坏。对于岩石等脆性材料而言, 东北大学岩石破裂与失稳研究中心研制开发出了相应的岩石渐进破坏分析软件RFPA。RFPA 软件通过考虑材料的非均匀性来模拟材料非线性, 通过单元的弱化来模拟材料的变形和破坏等非连续力学行为。该法在对裂纹扩展、岩层移动、震源孕育模式、微震活动和复合材料破坏等的数值模拟方面取得了良好效果。

耦合方法除了上述各种数值方法与计算模型外, 耦合计算方法也得到了发展与应用, 如有限元与边界元的耦合、有限元与离散元的耦合、离散元与边界元的耦合等, 这些耦合方法可分别发挥各种方法的优点并进行耦合提高。

这些数值方法中, 以有限元法为代表的连续变形分析方法在实际岩土工程中的应用较为广泛, 特别是对于岩体比较完整、节理裂隙密集且符合统计规律、裂隙不贯通的岩体更是如此。从学科发展来看, 在普遍意义上更符合实际岩体变形特点的非连续变形分析方法将具有广阔的发展空间。上述岩土力学数值方法的发展过程也表明了这种趋势, 但连续性分析、开裂以及非连续性分析的统一问题长期没有得到很好解决, 有待于进行更深入的理论研究3。

总结随着岩土工程领域的不断扩展与延伸, 数值分析方法得到了长足的发展, 分析模式不断改进, 分析精度不断提高。但岩土工程材料是一种复杂的地质材料, 具有高度非连续性、非均匀性以及各向异性等地质特点, 在力学性质上表现出强烈的非线性、非弹性等力学行为。因此, 只有不断改进数值分析的思路与方法, 使其能更好的模拟实际岩土材料的工程力学性质, 才能更好的指导岩土工程实践。

无单元类方法的无网格技术和流形方法的有限覆盖技术是近年来新兴数值分析方法中独具特色的新技术。前者无需有限单元网格划分, 对求解域的剖分仅用一系列点进行, 此法简单易行, 特别适于进行断裂扩展分析。后者经大量计算表明是一种切实可行的统一处理连续与非连续变形问题的有效途径。两种方法各有优点, 互补性很强。连续与非连续变形统一分析技术、裂纹扩展分析的无单元技术和非连续变形分析方法中的散体接触、运动处理算法等将为发展新兴的更为综合的数值算法提供极好的工具和思路, 将会有力的推动计算力学与计算岩土力学的发展2。