高等土力学高等土力学 (13).pdf

第 49 卷第 8 期2 0 1 6 年 8 月土木工程学报CHINA CIVIL ENGINEEING JOUNALVol 49AugNo 82016基金项目:国家自然科学基金(51479099，51179092，51379103)和国家“973”计划(2013CB036404)作者简介:张丙印，博士，教授收稿日期:2016-02-18土工构筑物和边坡工程发展综述 作用机理与数值模拟方法摘要:结合我国土工构筑物和边坡工程的建设实践，主要针对近年来有关作用机理与数值模拟方法方面的研究成果进行了讨论和总结，主要包括土工构筑物工作性态和变形时空发展规律、破坏机制和数值模拟方法、土-结构接触特性和多体接触分析、土质边坡(支挡结构)稳定分析和安全评价等。论文重点介绍了一些数值计算新方法在土工计算分析中的应用。关键词:土工构筑物;土质边坡;应力变形分析;破坏机制;接触特性中图分类号:TU4文献标识码:A文章编号:1000-131X(2016)08-0001-16Action mechanism and numerical simulation methods of soil structures andslope engineering:the state-of-the-artZhang Bingyin1Wen Yanfeng2Zhu Benzhen3Yu Yuzhen1(1 State Key Laboratory of Hydroscience and Engineering，Tsinghua University，Beijing 100084，China;2 China Institute of Water esources and Hydropower esearch，Beijing 100044，China;3 Northwest esearch Institute，Co，Ltd of CEC，Lanzhou 730000，China)Abstract:The research results of action mechanism and numerical simulation methods of soil structures and slopes inChina were reviewed and discussed The main contents include characteristics of temporal and spatial progress of workingstate and deformation of soil structures，failure mechanism and numerical simulation methods，contact properties of soil-structure interface and contact analysis of multi-body，and stability analysis and safety assessment of soil slopes andretaining structuresSome new numerical methods and the application to geotechnical analyses were introducedemphaticallyKeywords:soil structures;soil slopes;stress deformation analysis;failure mechanism;contact propertyE-mail:byzhang tsinghua edu cn引言土工构筑物与边坡专题涉及堤坝、尾矿库、路基、机场高填方、边坡和挡墙等。我国西部水能资源丰富，但均位于高海拔地区，地质条件复杂，地震烈度高，自然环境恶劣，交通不便。土石坝是当地材料坝，对基础适应性强，抗震能力高，造价低，是各种坝型中竞争力很强的坝型。特高土石坝主要包括土质心墙堆石坝和混凝土面板堆石坝两大坝型。我国 2000 年以前建设的土石坝多为100m 级的高坝。其中，最高的是小浪底心墙堆石坝(坝高 160m)和天生桥一级面板堆石坝(坝高 178m)。进入 21 世纪，已建成瀑布沟(坝高 186m)和糯扎渡(坝高 261 5m)等高心墙堆石坝，以及洪家渡(坝高179 5m)、三板溪(坝高185 5m)、水布垭(坝高233m)等高面板堆石坝。正在规划设计或建设长河坝(坝高240m)、双江口(坝高 314m)、两河口(坝高 295m)、如美(坝高 315m)、古水(坝高 244m)、茨哈峡(坝高253m)和马吉(坝高 270m)等高土石坝。我国土石坝建设正在由 200m 级逐渐向 300m 级的超高坝发展。我国已建尾矿库 12600 多座，每年新建尾矿坝300 多座，坝高、危险病库量大，300 米级高尾矿坝不断兴建。我国尾矿坝灾害频发，重特大事故高发，伤亡人数和财产损失呈增长态势，造成了极其严重的社会影响。此外，大量的路基和机场等高填方工程也正在DOI:10.15951/j.tmgcxb.2016.08.001 2 土木工程学报2016 年不断兴建。高速铁路和机场等对工后变形的要求严苛，变形预测和控制是这些高填方工程的重点难点课题。采用各种试验技术、分析手段与计算方法研究土工构筑物的工作性态，分析评价其安全性，是进行工程建设的前提。本论文重点综述近年来我国在土工构筑物和边坡工程机理分析与数值模拟分析方法方面所取得的一些进展。1工作性态和时空发展规律土工构筑物的工作性态主要是指其应力-变形-孔压的分布规律以及构筑物的安全性状态。土工构筑物工作性态的准确预测需要准确把握填筑土石料的变形特性，并建立合理的计算模型和方法。1 1填筑土石料变形机理与多尺度力学性能填筑土石坝等土工构筑物的主要材料通常是就地开采的土石料。其中，堆石料属典型的大粒径颗粒材料。组成颗粒形状各异且具无序、随机等诸多不确定性，其物理力学性质尚无法通过物理或数值模型准确模拟。堆石料的力学行为依赖于应力路径，还存在剪胀等特性，加上遇水软化、长期流变、颗粒破碎等复杂因素影响，其变形机理和本构关系仍然是未来进一步研究的重点课题。土的本构模型是现代土力学的核心内容之一1。自 20 世纪 60 年代以来，国内外学者们发展和建立了数以万计的基于各种理论的土的本构模型。这些研究工作无疑极大地促进了对土体基本特性认识的深度。但是，其中的大多数模型主要被用来解释和模拟试验中土体特性和现象，较少被用于土工构筑物的变形数值计算。也有不少在理论上看起来完美、具有复杂表达形式的模型，应用于土工计算时的精度却令人失望。目前国内在土工变形计算中，常用的本构模型主要包括邓肯-张模型、修正剑桥模型、基于 DP 准则或莫尔库仑准则的理想弹塑性模型、沈珠江双屈服面模型2 等。这些模型主要是建立在室内常规三轴试验的基础上。受制于试件尺寸的限制，室内堆石料试验存在较大的缩尺效应。此外，真三轴条件下的研究工作也不多见。在特高土石坝中，高应力水平可导致堆石体颗粒发生较为严重的破碎3-4，坝体填筑和水库蓄水过程可使得堆石料的应力路径非常复杂，这使得已有的这些模型具有较大的局限性。堆石坝坝体是典型的颗粒堆积体，其宏观变形是细观颗粒位置调整和破碎的结果。基于细观组构特征和细观变形机理研究堆石料的本构特性是近年来的一个新趋势。学者们采用多种方法研究了土石料的结构特征，发现土石料的粒度分布特征在宏观统计层次上表现出良好的自组织性5-6。随着数值计算技术的飞速发展，基于数值分析技术的岩土介质细观力学试验作为一门新兴的试验方法深受国内外研究者的青睐，为堆石料细观力学行为及变形破坏机理的研究提供了有力的手段7-9。不断发展的数字图像技术也已应用于岩土材料细观变形特性的研究。土石料的力学特性具有明显的多尺度效应。随着计算力学和数值模拟技术的飞速发展，多尺度模拟和计算逐渐成为一个正在迅速发展的热点和前沿研究领域，而且也逐渐被应用到岩土材料的研究中10。1 2坝体变形时空演变规律对于特高土石坝、高填方路基和机场等，变形控制无疑是工程建设中的核心关键问题。经过近年来200m 级高堆石坝的建设实践，我国在高土石坝变形控制方面取得了不少宝贵经验11-13。高土石坝的变形控制可概括为如下三个方面:变形总量控制，这是变形控制的基础;变形协调控制，使坝体各部分变形平顺，防止发生不利的应力或应变状态;有害变形控制，通过调控施工填筑过程和蓄水过程，让坝体变形尽量发生在对大坝安全影响相对较小的阶段。高土石坝工程规模宏大，结构复杂，其变形的演化规律受到坝料、分区、水位变动、环境、气候等众多复杂因素的影响。总体看我国 200m 级高土石坝工程的建设和运行历史尚短，相关的监测资料还比较缺乏，有关坝体变形时空演变规律的研究成果还较少。我国近期建设的几座高堆石坝在取得成功及宝贵经验的同时，部分工程出现了坝体变形偏大、发生较大规模坝体裂缝、面板挤压破坏、防渗体系破损、渗漏量较大等问题。一些高堆石坝在经过十多年的挡水运行后，其变形的量值变化也未稳定，运行期坝体或面板结构破损的情况时有发生。表 1 给出了国内外典型高面板堆石坝变形的沉降监测特性。为对大坝变形进行有效控制，关键是对不同条件下大坝的工作性态进行相对准确的预测14-16。在土石坝变形计算方面，现有的分析方法和计算程序许多是 20 世纪 80 年代针对百米级土石坝研制的，虽已成功地应用于许多高土石坝的建设，但也存在较大的问题。经验表明，目前土石坝变形预测的精度总体较低，且存在所谓的“低坝算大，高坝算小”的现象。特高堆石坝的建设要求更高的计算精度和更完善的分析方法，相关研究工作十分必要。第 49 卷第 8 期张丙印等土工构筑物和边坡工程发展综述 作用机理与数值模拟方法 3 表 1国内外典型高面板堆石坝变形特性Table 1Settlements of some typical CFDs in the world坝名坝高(m)最大沉降(cm)沉降坝高比(%)阿里亚160 0358 02 24阿瓜米尔巴186 0170 00 91天生桥一级178 0347 0200洪家渡179 5132 20 74三板溪185 5156 30 84水布垭233 0220 01 00工程经验表明，除了施工期的应力加卸载变形之外，非荷载变化(湿化和流变变形等)所导致的工后变形也是土工构筑物变形的重要组成部分。在这方面，相关的研究工作多是结合土石坝工程进行的。国外学者在 20 世纪 70 年代初提出了土石坝湿化变形的计算方法。国内湿化变形的研究起步于“七五”攻关期间，目前已建立起多种可进行土石坝湿化变形计算的模型和方法，如基于双线法的初应变求解方法、殷宗泽的增量初应力法、沈珠江湿化模型等。在土石坝流变变形方面，沈珠江17 最早进行了堆石料单级加载流变试验。根据试验成果，提出了流变变形速率随时间衰减的关系式以及最终流变变形与应力状态间的关系式等重要规律，并据此建立了在我国得到广泛应用至今的土石坝流变变形的计算框架。米占宽18、孙振远和李国英19、李海芳以及徐泽平20 等学者进行了类似的试验工作并提出了改进的模型计算公式等。分级加载流变试验更能反映实际土工建筑物分级填筑过程的特性。程展林和丁红顺21 针对多种堆石料进行了一系列的分级加载流变试验研究。陈晓斌和张家生22-23 针对红砂岩路堤填料进行了分级加载的流变试验。近年来，不少学者开始研究各种环境因素变化对堆石料变形的影响。王海俊和殷宗泽24 以饱水-自由排水的方式研究试样在饱和-非饱和循环下的变形特性。孙国亮、张丙印和张建红等25-27 研制堆石料风化仪探讨了堆石料的风化变形特性。高土石坝的实际施工过程和运行方式决定了其流变变形和加载变形在时间上是交替发生的。现有堆石料流变试验和计算模型，没有考虑流变变形与瞬时加载变形之间的相互影响。为此，钱晓翔28 等进行了分级加载三轴流变试验，研究了复杂加载条件下土体的流变特性及计算模型。1 3超静孔压演变规律和多场耦合分析在心墙流固耦合的计算中，一般认为心墙土料的渗透系数为常量，不随心墙的受力变形而发生变化。但是，在我国多座高土石坝的施工和运行过程中，实测到的超静孔隙水压力及其变化过程与计算值相差较大。例如，坝高 160m 的小浪底斜心墙堆石坝已运行 12 年，其心墙内某些部位的孔隙水压力至今仍然较高，水头甚至可超出坝顶高程29。坝高 261 5m 的糯扎渡心墙堆石坝已于 2012 年底完成主体工程施工，但施工期和运行期心墙内实测孔压值很高，比计算值总体偏高约一个数量级30。这些现象根据已有的理论和方法难以解释和模拟，其原因可能是在施工和蓄水过程中，坝体内的物态场、渗流场、变形场和应力场以及周围环境之间存在着复杂的耦合关系。1 3 1黏性土渗透性研究进展针对黏性土的渗透性，过去几十年学者们进行了大量的研究工作。关于孔隙比对黏性土渗透性的影响，Marshall31-32 等的试验研究结果表明，孔隙比大小对黏性土渗透性具有重要影响，渗透性随孔隙比的减小而迅速减小。刘杰33 针对黄河小浪底工程的黏性土进行了试验研究，提出了一个土体渗透系数与孔隙率的经验关系式，并认为影响压实黏性土渗透性的重要因素是压实含水量和压实功能以及由此引起的土的不同结构形态。谢康和等34 对萧山饱和软弱黏性土进行的固结试验表明，其渗透系数随着固结过程的进行而减小，且孔隙比与渗透系数基本满足对数关系。由于实际土体大多处于受力发生变形的过程中，因此，有学者从力学变形的角度对黏性土的渗透性进行了研究。Chan 和 Kenney35 等利用三轴仪对某黏性土进行了渗透试验，采用立方体试样在不改变受力条件的前提下分别测定平行于和垂直于层面方向的渗透系数。Carpenter 和 Stephenson36 采用三轴仪研究了某黏性土试样尺寸、水力坡降、渗透历时以及渗流出口处的有效应力等因素对渗透系数的影响，但同样没有考虑应力状态的影响。李平等37 采用三轴仪对试样等向固结后进行渗透试验，结果表明，固结后的渗透系数与对应的孔隙比在半对数坐标系下呈近似直线关系;土体渗透系数随固结度和周围压力的增大而减小。总体看来，关于黏性土渗透性的研究，主要集中于土体物理力学特性对其渗透性的影响。但这些试验研究中，土样大多没有受力变形，不能反映土体受力变形后的渗透特性。关于力学状态对渗透性影响的试验大多是在土体等向固结的条件下进行的，没有考虑土体实际上承受不等向应力而产生的剪切变形的影响。因此，关于土体发生剪切变形特别是大剪切变形条件下的渗透特性的研究还有待深入。1 3 2应力-变形-渗流多场耦合作用规律研究 4 土木工程学报2016 年雷红军等38 研制了三轴渗透试验装置，研究了应力和变形状态发生变化时压实黏性土的渗透特性。试验结果表明:同一围压下，在试样受力发生剪切变形的起始阶段，渗透系数迅速减小，随着轴向应变的增加，其变化速率越来越慢，最后基本上趋于稳定;低围压下试样的渗透系数变化幅度较大，最大变化幅度可达数十倍;不同围压下的试样，对应于相同的轴向应变，围压小的试样渗透系数较大。土体在受力和发生变形后渗透性变化的机理在于:试样的体积压缩引起土体孔隙比的变化;剪切引起土颗粒的移动错位及重新排列，即土体结构发生变化，可用剪应力水平的变化描述。建议用下述指数函数的数学模型描述黏性土大剪切变形后的渗透系数:k=k0exp(ae+bSl)(1)式中:k 为试样的渗透系数;e 为孔隙比;Sl为试样的剪应力水平;k0、a 和 b 为试验系数。可将拟合公式计算值与试验实测结果进行比较。图 1 给出了土体渗透系数实测值与模型计算值的比较，二者符合较好，表明所提出的渗透系数数学模型具有良好的适用性。1 3 3考虑多场耦合作用的超静孔压计算方法在传统基于 Biot 固结方程的有限元计算中，大多假设渗透系数恒定不变。可将上述渗透系数模型引入 Biot 渗流固结理论，以反映黏性土应力变形场与渗流场之间的相互影响。图 1不同围压三轴试验中渗透系数随孔隙比的变化Fig 1elationship between permeability and void ratio采用上述计算方法对两河口水电站心墙堆石坝进行了计算分析。结果表明，相比初始渗透系数，心墙底部渗透系数相对可降低近 100 倍。吴永康等39 结合糯扎渡心墙堆石坝的坝体监测资料，考虑应力变形场与渗流场耦合作用，对渗透模型参数进行了反演分析。图 2 是心墙典型渗压计 2012 年 11 月 10 日监测值与根据反演参数得到的计算值的空间分布对比图。总体来看，考虑多场耦合作用的计算值与监测值能较好的符合，而不考虑耦合作用的计算结果则与监测值具有明显的差别。图 2典型渗压计监测值与计算值分布对比Fig 2Distribution of calculated and observed water head2破坏机制和数值模拟方法土工构筑物的破坏状态对应了土体变形或承载的极限状态。在该状态下，土体通常伴随剪切带和裂缝扩展等应变局部化现象，以及大变形和非稳定等问题。在传统土力学研究和土工数值模拟计算中，这些问题通常是难点研究课题。近年来，各种数值模拟计算技术的迅猛发展，为土体破坏现象的模拟提供了有力的手段。本节主要介绍近年来在土体张拉裂缝和剪切带模拟等方面的一些主要进展。2 1土体的张拉断裂特性与裂缝模拟计算裂缝是土工构筑物常见的隐患和引起发生破坏的主要原因之一。土体在绝大多数情况下处于受压状态，其在受拉条件下的特性长期未受到重视。土体裂缝通常是土工构筑物发生破坏的预兆，因此，研究其发生扩展机理，发展其模拟计算方法具有非常重要的意义。2 1 1土体张拉特性和裂缝计算研究进展土是碎散颗粒的集合，颗粒之间的相互联系薄弱，其强度主要由颗粒间的相互作用力所决定。土的强度受到其外部和内部、微观和宏观众多因素的影响，是十分复杂的研究课题。目前，一般采用单轴拉伸试验、三轴伸长试验、土梁弯曲试验、径向和轴向压裂试验以及空心圆柱内压劈裂试验等研究土体的张拉强度特性。一些典型的研究工作包括:丁金粟40、张小江41、余湘娟等42、骆亚生43、Tang 和 Graham44、党进谦45、张少宏等46、Nahlawi 等47、郭飞等48、朱安龙49、王俊杰等50、Tamrakar 等51 等。综合文献中的相关研究工作，有关土体拉伸特性的研究主要集中在通过试验研究其抗拉强度和峰值应变的特性。由于试验方法的限制，一些试验采用了第 49 卷第 8 期张丙印等土工构筑物和边坡工程发展综述 作用机理与数值模拟方法 5 间接拉伸方法。早期的单轴拉伸试验设备多采用应力控制式加载方式，无法测定土样在达到峰值抗拉强度后拉断过程中的应力应变关系曲线，故通常会得到土体受拉时发生脆性破坏的结论。在基于线弹性断裂力学研究方面，多集中于针对断裂韧度 KIC及其影响因素的研究层面开展工作，对裂缝发生和扩展过程模拟计算方法的研究工作开展相对较少。裂缝产生的主要原因是土工构筑物的不均匀变形。目前用于判断是否会发生裂缝的方法主要有变形倾度法和有限元方法等。前者是一种经验的方法。彭翀和张宗亮等52 通过在有限元计算程序中嵌入变形倾度计算模块，发展了基于有限元变形计算的变形倾度有限元法。长期以来，在土工数值计算分析中很少关注黏性土的张拉破坏及裂缝问题，目前尚没有成熟的数值模拟计算方法。2 1 2土体张拉断裂特性和裂缝扩展机理研究在总结国内外黏性土抗拉特性试验设备经验的基础上，清华大学岩土所53-56 研制了能用于测定黏性土拉伸全过程曲线的位移加载式刚性试验设备，包括应变控制式单轴拉伸仪、位移控制式土梁弯曲仪和位移控制式土工三轴拉压仪。使用研制的上述试验设备，对糯扎渡、双江口和两河口高心墙堆石坝心墙黏土料进行了大量的单轴拉伸试验、断裂韧度试验、三轴拉伸试验和土梁弯曲试验。根据试验成果研究了黏性土的张拉断裂特性。将黏性土材料的裂缝扩展过程划分为微裂、亚临界扩展和失稳扩展三个阶段。图 3 给出了拉伸试验断裂发生后试样的照片和特性。图 3黏性土张拉断裂区及特性Fig 3Characteristics of extension area in clayey soil压实黏土张拉断裂区顶端可分为近似完好区、裂缝损伤扩展区和宏观裂缝区三个区域。土体张拉断裂过程符合准脆性材料的基本特征57。在土样发生张拉断裂的过程中，断裂区土体和非断裂区土体表现出完全不同的应力变形特性。在非断裂区，土体承受的轴向拉力减小，发生卸载回弹。在断裂区，土体发生破损，裂缝逐步发展，是土体张拉变形的集中发生区，只有张拉断裂区土体的应力应变关系才真正是土体抗拉应力应变的全过程曲线。根据上述分析，建议了受拉条件下应力应变全过程曲线的整理方法，发展了拉伸条件下土体的非线性本构模型。根据多种土料三轴拉伸试验结果，发现随室压大小的不同，土样会发生如下三种不同形式的破坏58:室压较小时的张拉破坏;室压较大时的剪切破坏;中间围压范围的张拉-剪切复合破坏。为了能够统一描述土体张拉-剪切复合破坏的现象，借鉴闭合应力点概念，提出了一个新的联合强度准则。土石结构局部土体的受拉状态总是由受压状态演化而来，进行分析时需采用能够综合考虑土体拉压变形特性的本构模型。张琰等59 根据压实黏土的常规三轴试验和三轴拉伸试验结果，发展了能够统一考虑黏性土拉压特性的邓肯-张扩展本构模型。2 1 3基于有限元-无网格法耦合的裂缝模拟方法裂缝是一种位移不连续现象，在数值计算中需要引入一定的方法对其进行处理。无网格法易于处理裂缝等移动边界问题，在土体裂缝模拟计算方面具有优势。但是，无网格法的计算效率较有限元方法为低。为了发挥无单元法和有限元法各自的优势，张琰55、彭翀等58 针对土体三维张拉裂缝问题，发展了基于无单元法的弥散裂缝模型，提出了径向基点插值无单元法与有限元直接耦合的计算方法。在实际土工构筑物计算中，采用的计算网格一般尺寸较大，难以准确模拟裂缝可能的起始位置和发展方向。为此，提出了在裂缝计算过程中自动判别可能开裂区域、进行有限元网格至无网格节点的转化以及节点自动加密的方法。计算开始时使用应力变形分析常用的有限元模型，在模拟过程中通过一定的标准来判断裂缝可能开裂区，并将其从有限元网格转化为无网格节点，然后进行局部自动加密。网格依赖性问题是裂缝计算中的难题。对此发展了基于无网格法和开裂势函数的三维复杂裂缝追踪算法60。该算法可通过一次计算获得计算域内所有可能的开裂方向。厦门翔安隧道连接厦门岛东部的五通和对岸翔安大陆的西滨，是中国大陆地区第一座大断面海底隧道。在隧道开挖过程中，地层发生了较大的变形，造成地表发生了三维复杂形态的张拉裂缝。采用发展的基于有限元-无网格耦合的土体裂缝三维模拟方法进行了计算分析。图 4 给出了典型结果，可见裂缝的发生形态与实测裂缝非常相似。小浪底水利枢纽大坝为壤土斜心墙堆石坝，最大坝高160m。2001年7月，在该坝坝顶下游侧发现一 6 土木工程学报2016 年图 4地表裂缝开裂形态与模拟计算结果Fig 4Observed and simulated cracks near ground surface条长约 100m，最大宽度约 10mm 的非连续纵向裂缝。2004 年 8 月测得裂缝长度 627m，裂缝缝口宽度 2 100mm。采用建立的三维张拉裂缝计算方法对该裂缝的开裂过程进行了计算分析。结果表明，由模拟计算所得裂缝的开展过程和发生形态与实际观测结果基本一致。2 2土体剪切带形成机理与模拟计算2 2 1土体剪切带形成机理和模拟计算研究进展人们很早就认识到土体渐进破坏特性与剪切带扩展过程的密切关联。20 世纪 60 年代，Skepmton61 和 Bjerrum62 等对超固结黏土边坡渐进破坏问题进行研究时，认识到剪切带上的强度并不是同时发挥的，并指出土体应变软化特性以及由此引发的剪切带渐进扩展过程对于土坡渐进破坏问题研究的重要意义。此后几十年，剪切带萌生和扩展问题一直是岩土工程领域的研究热点之一。Palmer 和 ice63 等提出的基于断裂力学方法预先假定剪切带的扩展方向，仅仅考虑已有剪切带尖端的剪应力集中，不考虑应力重分布导致的主应力轴偏转现象，总体上为纯概念模型，难以实用。为解决这些问题，学者们求助于描述和处理不连续位移场的应变局部化理论，如 Cosserat 理论64 和梯度塑性理论65 等。但是，与常规有限元结合时，这些理论常需要改变裂缝区材料的本构关系，以反映场函数的非连续性。这类方法通常应用于弱不连续位移场的模拟，其剪切带路径受网格形状的影响大，且存在应力锁死现象。嵌入非连续方法66 在单元内嵌入不连续的位移场，用于分担单元总体应变，从而允许单元产生较大畸变，在网格依赖性和应力锁死方面有改进。然而，这类方法并不协调一致，剪切带尖端难以确定，也就无法描述尖端应力集中和重分布对剪切带扩展的影响。扩展有限元法(XFEM)67 可将物体的材料非线性和土体结构的边界非线性区别对待，既可较好地描述内部运动的不连续位移场，又可相对独立地选择合适的本构模型，具备模拟剪切带等运动的不连续位移场的能力。与嵌入非连续方法相比，该方法网格依赖性较小，对剪切带形状和界面接触状态的描述也更准确。综合分析，XFEM 整体上较优。但目前 XFEM 等多用于模拟金属、混凝土和岩石等材料中的张拉裂缝，较少用于土质材料中的剪切带问题，其自身还需要进一步改进和完善。此外，便于嵌入 XFEM 的剪切带扩展准则也有待研究。2 2 2土体剪切带形成机理试验研究喻葭临等68 采用压力室模型箱研究了人工结构性土体中剪切带扩展机理。试验考虑了无预设剪切带、两种预设剪切带倾角和两种加载方式，共进行了 5个模型试验。由试验现象可知，预设剪切带的存在使剪切带并没有突然贯穿整个研究区域，而是存在一个扩展的过程。剪切带前端区域内的位移场经历了从均匀变形到弱不连续性最后到强不连续的过程。综合比较 5 个试验结果可知，土体内部材料性质的不均匀性和变异性，比如材料软化形成的弱面等，会导致尖端的应力集中和重分布，从而推动已有剪切带向前扩展，且只需较小的外部荷载作为驱动力。预设剪切带试样的扩展区域存在明显的共性。剪切带无论以剪切型裂缝或是张拉型裂缝向前扩展，均发生在同侧，如图 5 中的半圆区域所示，而在尖端另一侧的土体始终保持均匀变形状态。显然，这种状况仍然是尖端应力集中和重分布的结果。(a)SB-BA-In(b)SB-BA-De图 5剪切带扩展区域的共性Fig 5Evolution area of the shear band2 2 3尖端应力集中和扩展控制区域为掌握更一般的扩展规律，需要对剪切带尖端的应力集中和重分布状态做更全面的考察。在如图 6 所第 49 卷第 8 期张丙印等土工构筑物和边坡工程发展综述 作用机理与数值模拟方法 7 示的土坡中，假定由于某种原因已经存在剪切带。若应用 Bjerrum 的渐进扩展理论，只需考察 O 点或者已有剪切带延长线上的 A 点是否达到破坏状态，只用到剪切面上的正应力和剪应力，而忽略了其它应力的影响。但当尖端的应力集中比较严重时，尖端两侧的应力状态有较大差异。将 A、B、C 三点的应力状态用 Mohr 圆绘出，如图7 所示，从中可知 B 点较 A、C 点更有可能先达到破坏状态，这使得 B 点的应力状态在控制剪切带扩展时刻和方向上较 A、C 点占优。而这很好地解释了模型试验预设剪切带扩展区域的共性特点。图 6土坡中的初始缺陷示意图Fig 6Schematic of the initial defect in the slope图 7尖端附近各点对应的 Mohr 圆Fig 7Mohr s circles of typical points near the crack tip结合试验结果和以上分析可以推论，剪切带的扩展具有区域选择性，最终控制剪切带扩展的是尖端附近的局部区域而非全域。以往的分析往往将剪切带尖端附近的全域作为控制区域，实际上这是忽略这种差异性和非连续性的平均化手段。但当这种差异性足够大时，平均可能得不到合理的结果。为方便叙述，将剪切带尖端扩展的控制区域称为扩展控制域。该方法中提出的扩展控制域具有区域平均的意义，且具有一定程度的可选择性。在数值计算过程中，不再将尖端附近如图 6 所示的 O 附近所有区域都用来计算，而是可在已有剪切带尖端附近选择若干区域，如计算以 A、B、C 为顶点的三个半圆所占据的区域的加权平均应力状态，然后结合强度理论选择其中最先达到强度的区域平均应力并以此计算剪切带的扩展方向。可见，该分析方法的概念较简单且容易与数值模型相对应，可较好地反映该方法的区域控制特点。2 2 4土中不连续位移场的模拟方法扩展有限元法在描述连续-非连续位移场构造上具备显著的优势，并在张拉型裂缝扩展问题的模拟研究中取得了成功。然而，应用于描述土中剪切带时需要改进扩展有限元法中非连续区域的积分方案以便保障积分点上非线性材料内变量的一致性;同时构建基于扩展有限元法的界面接触算法以便合理地模拟剪切带形成后的界面行为。Song 等69 在模拟平面裂缝问题时，使用双线性四边形单元和 1 点减缩积分方案，提出一种无需分片积分的简易积分方法。但由于每个单元仅有一个积分点，为保证矩阵非奇异，还需一定的控制措施以避免出现零能模式，这增加了建模的复杂性。为此发展了一种更具普适性的简易积分方案。该积分方案采用了双线性四边形单元和常规 4 点积分方法，不用分片积分，仍采用常规高斯积分方法即可获得矩阵各系数最终结果。此时，常规自由度对应积分项的有效积分点仍是 4 个，可以有效避免零能模式的出现。该积分方案在不连续区域仍然使用和常规有限元方法一致的积分点，无需增加新的积分点，积分点上反映材料非线性的内变量的一致性可得到保证，常规有限元中适用的各种应力积分方法可以直接应用于扩展有限元方法中非连续区域的积分，这就较好地解决了扩展有限元方法对于非线性材料的描述和兼容问题。通过引入恰当的非连续区域本构关系构建与扩展有限元相适应的界面接触算法。与传统的接触力学方法相比，该方法无需增加特殊的界面单元，也不涉及复杂的接触条件变化和接触状态搜索，原理简单，建模容易。张剪是指裂缝面上接触点对在法向远离的同时切向距离也有变化的状态，此时的裂缝亦称为复合型张拉裂缝，可采用黏聚裂纹模型70 来描述。压剪对应摩擦接触状态，是更为复杂的现象。Willner71 提出了基于罚函数法的一种连续的摩擦定律。在静摩擦和滑动摩擦的过渡过程中，摩擦力被认为保持连续。2 2 5土中不连续位移场的模拟算例Carsington 土堤坝是位于英国中部的 Derbyshire的一座黏土心墙坝，坝轴线长 1200m，最大坝高 37m。1984 年 6 月 1 日，坝顶已经填筑至距离设计高程不到 8 土木工程学报2016 年1m。因监测水平位移大和下雨而停工，6 月 4 日早晨在坝顶观察到沟壑状的张拉型裂缝，滑坡随后迅速向上游方向滑动，到 6 月 6 日才最终停止。该滑坡体的初始冠顶长约 190m，最终冠顶长约 500m。对于 Carsington 堤坝的失稳过程，Skempton 等人72 经过调查后认为，滑裂面首先在心墙中出现，并产生向下和向上游的滑动，促使滑坡体沿着黄黏土夹层滑动，并刺穿盖重边墙形成巨大位移。Potts 等73 基于 Skempton 的分析，利用包含软化机制的弹塑性模型和常规有限元研究分析了堤坝的失稳机理。喻葭临等68 模拟了 Carsington 堤坝的失稳过程。计算中考虑了分层填筑坝体以及黄黏土夹层的软化和错动。数值模拟较好地追踪了 Carsington 堤坝心墙中的剪切带路径和形成过程，计算结果与调查结果吻合良好。图 8 给出了最终的破坏形态，与 Skempton 和Potts 等的分析结论一致。图 8Carsington 堤坝最终破坏状态Fig 8Failure state of the Carsington Dam3土-结构接触特性和多体接触分析在土工构筑物中经常会遇到不同材料的交界面，也即不同材料的接触界面。当接触界面两侧材料的特性存在较大差异时，在界面两侧常会出现较大的剪应力并发生位移不连续现象，从而导致十分复杂的应力和变形性态。土与结构接触面是土与结构之间应力和变形相互传递的重要媒介，在土木、水利、建筑、交通、市政等工程中普遍存在。如何合理地描述这些接触界面的力学特性是多年来得到广泛关注的热点研究课题之一。3 1接触面力学试验和接触面本构模型当前针对土-结构接触特性的研究，多集中于接触面特性试验研究和接触面本构模型研究等方面。在这些方面，清华大学岩土工程研究所的张建民研究团队从不同的角度完成了一系列具有代表性的工作74-77。该研究团队先后研制了大中小型接触面力学试验系统，系统探讨了各向异性、剪切体变、渐进破坏、静动力加载、复杂加载路径等众多复杂因素对接触特性的影响，发展了系列化的三维接触面弹塑性数值模拟技术。在接触面单元的基本框架内，殷宗泽等78、雷晓燕等79、李守德等80 也对接触特性和接触面本构模型进行了研究。3 2传统接触问题模拟方法在岩土工程数值计算中，目前常用的数值模拟计算方法主要是以 Goodman 单元81 为代表的无厚度类型的接触面单元和以 Desai 薄层单元82 为代表的有厚度类型的接触面单元。Goodman 接触面单元法能够模拟接触面的滑移和张开，概念明确，应用方便。该法自提出至今已在岩土工程计算领域得到广泛应用近50 年，积累了丰富的经验。但是，该法在法向劲度系数的处理等方面存在一定的缺陷，易造成迭代过程的震荡等。薄层单元可以模拟土与结构物接触面的黏结、滑动、张开和闭合等各种接触状态。但是该法也存在单元厚度的选取不易合理确定等问题。前述两种接触面单元法是目前在岩土工程应力变形计算分析中普遍采用的方法。该类方法将变形特性相差较大的不同材料分别当作同一个连续体的两个不同材料分区，并在两者之间设置接触面单元以模拟材料性质的突变。在这类方法中，对位移和应变的描述仍然是以连续介质力学的基本模式为基础，侧重于描述不同材料分界面上应力的传递，而对于材料界面上位移不连续现象的描述则较为粗糙，尤其是对发生大规模滑移和脱开问题的描述常导致计算收敛性不好和计算结果发生震荡等问题。近年来，基于非线性接触力学的模拟方法得到了迅速的发展，并被逐步应用于岩土工程领域83-84。3 3现代接触计算力学方法及实现接触问题是非线性数值分析中最具挑战性的课题之一。计算接触力学方法自提出以来就一直受到广泛关注，目前仍是数值计算领域内的研究热点，不断有最新的研究成果涌现。直观地讲，接触问题反映的是不同物体之间的相互作用。计算接触力学将发生相互接触作用的不同物体处理为各自独立的接触体。3 3 1接触问题的描述和虚功原理两个物体 A 和 B 接触的一般情形如图 9 所示。接触边界条件可概括为:不可贯入条件，即运动过程中独立的固体间不会发生相互贯穿;法向压力条件，是指在不考虑界面黏结作用的情况下，接触面之间的法向应力不能为拉应力;摩擦力条件，是指接触点对间的切向相互作用力应由法向压力和相对运动模式决定。在有限元范畴内，现代接触分析算法的核心在于构建相应的虚功原理。考虑 A-B 两物体间的接触问题，忽略一般的阻尼因素，将接触面边界视为面力边第 49 卷第 8 期张丙印等土工构筑物和边坡工程发展综述 作用机理与数值模拟方法 9 图 9物体间的相互接触Fig 9Contact of bodies界，则该问题在 t+t 时刻位形下的虚功原理可以表示为:A，Br=t+tWrintt+tWrLt+tWrIt+tWrC=0(2)式中:Wint、WL、WI和 WC分别表示内能、外荷载虚功、惯性力虚功和接触力虚功的变分。特别地，将接触力虚功的变分 WC展开如下:t+tWC=A，Br=t+tSrct+tFriurit+tdS(3)式中:t+tF 为接触力分布，u 为质点位移的变分。3 3 2接触非线性方程的求解方法接触分析区别于一般有限元分析的根本问题是需要求解接触力虚功，其主要困难来自于接触非线性，即平衡时刻的接触边界条件在计算前是未知的。目前常用的非