应变软化岩体分析原理及其应用.pdf

第 35 卷第 3 期 岩 土 力 学 Vol.35 No.3 2014 年 3 月 Rock and Soil Mechanics Mar.2014 收稿日期：2013-02-20 基金项目：国家自然科学基金资助项目（No.51179185，No.41130742）。第一作者简介：王水林，男，1965 年生，博士，研究员，主要从事岩土工程数值方法及其应用研究。E-mail： 文章编号文章编号：10007598(2014)0360915 应变软化岩体分析原理及其应用应变软化岩体分析原理及其应用 王水林，郑 宏，刘泉声，郭明伟，葛修润 （中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，武汉 430071）摘摘 要要：应变软化是指应力-应变曲线中轴向应力随应变的增加而减小的现象，许多种类的岩土介质在工程扰动的作用下呈现应变软化的行为。在分析应变软化问题时，其应力-应变关系式中的切线刚度矩阵是非正定的，由此导致计算求解的困难。将岩体应变软化过程简化为一系列脆塑性过程，于是应变软化问题的求解归结为一系列脆塑性过程的分析。基于经典弹塑性力学理论，提出了应变软化过程模拟方法及其相应的有限元求解过程，编制了计算程序，研究了应变软化本构模型中不同强度弱化速率对圆形洞室围岩塑形区分布的影响，进一步分析了应变软化模型对应的隧道径向变形沿洞轴方向的分布特征，并与已有监测数据得到的分布规律进行了对比。初步的研究结果表明，应变软化模型得到的计算结果是比较合理的。关关 键键 词词：应变软化；脆塑性模型；圆形隧道；数值解；摩尔-库仑屈服准则 中图分类号中图分类号：TU452 文献标识码文献标识码：A Principle of analysis of strain-softening rock mass and its application WANG Shui-lin，ZHENG Hong，LIU Quan-sheng，GUO Ming-wei，GE Xiu-run（State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan 430071,China）Abstract:In stress-strain curve,the process of decline of uniaxial stress at increasing strain is defined as strain-softening.Many types of geomaterials behave in strain-softening way in the case of disturbance of engineering activities.When the stress-strain relationship is described in mathematical formulation,the matrix of tangential elastic moduli is no longer positive-definite.Thereby difficulties arise in finding the solution to strain-softening problem.Strain-softening process is simplified into a series of brittle-plastic behavior;and solving brittle-plastic problem comes down to obtaining a series of brittle-plastic solutions.On the basis of classic plasticity theory,the method for analyzing strain-softening behaviour is proposed;and the corresponding solution process is implemented in finite element code.Furthermore,the influences of different strength-weakening modes on the distribution of plastic zone in the surrounding rock mass of a circular tunnel are studied.The longitudinal deformation profile(LDP)obtained by strain-softening constitutive model is also analyzed;and it is compared with the measurements in a tunnel.The preliminary numerical analyses indicate that the presented results look reasonable.Key words:strain-softening;brittle-plastic model;circular tunnel;numerical solution;Mohr-Coulomb yield criterion 1 引 言 在岩土工程中，工程岩土体大多呈现应变软化现象，室内单轴（或三轴）力学试验和现场大尺寸压缩（或剪切）试验很好地佐证了岩土介质的这种力学特性。为了较为准确地分析工程岩土体的变形与受力情况，国内外众多学者采用不同的方法、从不同的角度对应变软化行为的发生机制、分析理论及方法进行了系统的研究，取得了卓有成效的成果。综合而言，应变软化研究方法可以划分为细观和宏观分析方法。针对岩土体软化与脆性破坏问题，从细观角度出发提出的模型与方法有损伤力学模型1、细观力学模型2、强度统计分布方法34和物理细胞自动机方法5等，这类方法主要考虑材料颗粒在外载荷作用下引起的变形、错位、损伤和破坏，由此研究其宏观力学响应。这些成果能很好地描述变形局部化的发生与发展机制，解释变形局部化现象，再现应变软化过程中力与变形的关系。由于该类方法中有的引入了一些难以准确确定的参数，理论与方法还有待于完善。而且，其中有些方法离工程应用还有一段距离。610 岩 土 力 学 2014 年 从研究岩土体宏观变形与破坏现象出发而提出的模型与方法多基于塑性理论。经典弹塑性力学理论在分析岩土体非线性行为和工程稳定性时得到了广泛应用，但当材料变形过程中出现应变软化时，在峰值后区，问题的求解面临两大困难：其一为边值问题出现病态。介质强度进入峰值后区后，研究对象的本构方程特性发生变化，应力增量-应变增量关系矩阵为非正定矩阵，描述问题的偏微分方程失去椭圆型方程的性质。其二为数值分析过程中计算结果对网格划分的依赖。当采用数值方法求解该类问题时，问题的病态性以解的不稳定和解的网格依赖性出现，得到的计算结果有时变得没有物理意义。例如在有限元求解过程中，当网格划分越来越精细时，塑性区（或破损区）有可能会集中在一个非常狭小的范围，由此使得计算结果失去客观性。关于应变软化本构关系与相关问题求解，文献6对峰后本构模型的研究进展进行了探讨，文献7对建立在经典弹塑性理论基础上的应变软化问题的理论及其求解方法进行了深入分析，阐述了各种不同计算方法的优缺点，指出了需要深化研究的内容。为克服经典的塑性力学理论求解应变软化问题的缺陷，国内外学者考虑材料尺寸效应，提出了各种不同的理论与方法，如梯度理论817（包括应变梯度理论和内变量梯度理论）、偶应力理论（Cosserat理论）1819、非局部应变理论2024、扰动态理论25、复合体理论2628、内嵌不连续方法2934、复合单元方法35、自适应有限单元法36等。在这些理论与方法中，大多数都引入了一个与材料颗粒尺度有关的几何量，目的是消除问题的病态性和数值解的网格依赖性。这些模型、理论和方法的提出及应用极大地推进了应变软化问题的研究，对上述工作的特点与发展，文献37进行了详细的分析和评述。基于经典弹塑性力学理论，本文将应变软化过程简化为一系列的脆塑性过程，由此，应变软化问题的求解就归结为一系列脆塑性问题的求解。提出了应变软化本构模型分析方法及其相应的有限元求解过程，编制了计算程序，并尝试在工程中进行应用。从理论上讲，本文方法的核心是求解脆塑性问题；从数值计算结果来看，数值解虽然无法消除网格依赖性，但宏观上得到的变形和塑性区与半解析解的结果是吻合的，与其他方法的解是基本一致的。2 经典弹塑性理论本构矩阵特性介绍 对弹塑性介质，假设其屈服准则表示为 (,)0F （1）式中：为应力张量；为硬化或者软化参数。对岩土类软化材料，通常取塑性剪应变。与屈服函数对应的塑性势函数为(,)Q。材料进入塑性状态后，在当前载荷增量下，如果加载条件成立，则 T()d0FD （2）式中：D为弹性应力-应变关系矩阵，d为应变增量。那么，在增量载荷作用下，介质的塑性应变增量为 pddQ （3）式中：p 为塑性应变张量；d为待求的塑性乘子。总应变增量（d）可以分解为弹性应变增量（ed）与塑性应变增量（pd）之和，即 epddd （4）进一步可以导出塑性阶段应力增量（d）与应变增量之间的关系式为 epdd D （5）其中，Teppp1QF,MDDDDDD （6）TFQMAD （7）1ddFA （8）对硬化材料 A0，对理想弹塑性材料 A=0，对应变软化材料 A0。文献38阐述了epD矩阵的特性，对关联流动法 则，从数学上证明了当材料在塑性阶段表现为应变软 化时，epD为不定矩阵。这直接说明了采用有限元方 法形成的应变软化介质的切线刚度矩阵是不定的。文献39-40也阐述了无论对于遵守关联流动法则还是遵守非关联塑性流动法则的应变软化材 料，当应变软化速率超过某一个特定值时，epD矩阵变得没有意义。应变软化材料epD矩阵的上述特 性导致了应变软化问题求解的困难。3 脆塑性分析方法简述 基于经典塑性力学理论，文献39,41从宏观唯象的角度提出了脆塑性介质的分析原理与方法，为研究脆塑性岩体的力学行为与相应工程的稳定性提供了一条非常有效的途径。为了本文的完整性，下面对脆塑性问题的求解方法39进行简要阐述。3.1 增量边值问题及其变分提法增量边值问题及其变分提法 岩土体的工程行为是一个与时间和过程密切第 3 期 王水林等：应变软化岩体分析原理及其应用 611 相关的力学问题（本文暂不考虑时间因素的影响），其边值问题按增量方式给出较为简便39。下面的表达式中字母前的符号“”表示与增量载荷对应的有限增量以区别于上一部分中微分增量的表述方法。设研究区域内岩体在某阶段处于平衡状态后的位移、应变与应力张量形式为u、，其分量分别为ui、ij、ij。下一个施工作业导致岩体内位移、应变与应力的增量为ui、ij、ij。应力增量应满足平衡条件,0ij jib （9）总应力应满足屈服条件,0F （10）位移增量与应变增量应满足几何方程,1()2iji jj iuu （11）应力增量与应变增量应满足本构方程 dep0d D （12）此外，位移增量与应力增量还必须满足边界条件 ijjiiinTuu （13）式中：nj为边界面法向方向余弦；iT为边界给定面力分量；iu为边界面给定位移分量。在材料交界面或结构面上 21()0ijijjn （14）为便于数值求解，与边值问题对应的变分方程为 TTTdddVb VT S uu（15）式中：表示相应变量的变分。对上式进行有限元离散后就可以得到增量形式的有限元方程如下：TTepT()ddd0eeeeeeSSubTS B D BNN（16）式中：为研究区域；e为划分的单元区域；S为力边界；eS为力边界单元；B为形函数导数矩阵；N为形函数矩阵；b为体力增量；T为面力增量。位移增量u为方程的基本未知量。一般采用迭代方法求解该非线性方程组。3.2 满足满足 Mohr-Coulomb 准则准则的脆塑性过程应力的脆塑性过程应力增量计算增量计算 在脆塑性过程中，应力增量（）是应力跌落增量（d）与塑性流动产生的应力增量（f）之和，即 df （17）关于应力跌落增量（d）计算，文献39认为，应力跌落过程中最大主应变保持不变，即1 0，提出了一种计算应力跌落增量的方法。下面假设在应力跌落过程中最小主应力保持不变（即3 0）的前提下，给出应变软化过程中服从Mohr-Coulomb屈服准则时应力跌落增量的计算方法。记应力张量的分量ij（i=13,j=13）在主应力空间中的主应力为(1,3)ii，Mohr-Coulomb屈服准则在主应力（规定压应力为正）空间中表示为 1313(,)0FY （18）式 中：(1sin)/(1sin)；2 cos/Yc(1sin)，c和分别为黏聚力和内摩擦角；1和 3分别为第一和第三主应力。若 峰 值 强 度 和 残 余 强 度 主 应 力 分 别 为p(1,3)ii和r(1,3)ii，则有 p1p3p1pp3ppr1r3r1rr3rr(,)0(,)0FYFY （19）式中：ppp(1sin)/(1sin)，ppp2cos/Yc p(1sin)，rrr(1sin)/(1sin)，rr2Yc rrcos/(1sin)。假设应力跌落过程中最小主应力不变，有 d33r3p0 （20）由峰值强度面和残余强度面应力关系式，可以得到 d11r1prp3rrp()YY （21）图1为主应力空间中对应于峰值强度与残余强度的应力摩尔圆，进一步认为应力从峰值强度跌 落至残余强度时Lode参数不变，即2r3r1r3r 2p3p1p3p，则有 2p3pd22r2p3r1r3r2p1p3p()（22）由此得到了主应力空间中的应力跌落增量为 drp(1,3)iiii （23）如果在应力跌落前后主应力方向保持不变，应 力跌落增量dij就可以方便地由di转换得到。塑性流动产生的应力增量fij通过对应于残 余强度的本构积分计算得到 fep0dD （24）612 岩 土 力 学 2014 年 式（23）中的应力跌落项di为0时，与脆塑性问题相关的方程就蜕化为理想弹塑性问题的对应表示式。(a)对应于峰值强度的应力摩尔圆 (b)对应于残余强度的应力摩尔圆 图图 1 应力摩尔圆应力摩尔圆 Fig.1 Mohr circles 3.3 脆塑性过程应力增量求解脆塑性过程应力增量求解 脆塑性求解过程与通常的弹塑性分析过程基本相同，其差异之处仅在于应力增量计算时增加了应力跌落增量的计算，这里给出迭代计算过程中包含本构积分计算的应力增量求解过程。对某次载荷增量的第k次迭代，迭代计算过程按照下述步骤进行：（1）计算弹性应力增量kkeD和试应力1kkkee。（2）检验试应力ke是否在峰值强度面（对应于上次载荷增量结束，且计算收敛时得到的各高斯点强度参数值，包括黏聚力、内摩擦角和剪胀角等）之外。如果在峰值强度面之内，第k次迭代对应的应力为kke；如果超出峰值强度面，则进行脆塑性分析。（3）对脆塑性过程，计算应力跌落增量d和 塑性流动引起的应力增量fep0dD，k df；对已经达到残余强度值的高斯点，只需要计算塑性流动应力增量f，即fk。（4）更新节点位移以及高斯点应力（k 1kk）。（5）计算残余力增量，判断迭代计算是否收敛；如果没有收敛，重新进入步骤（1），进行下一次迭代计算。（6）计算塑性应变，更新强度参数，检查两次更新的强度参数之差是否小于给定误差（即是否收敛），如果不小于给定误差，由本次迭代收敛计算得到的各高斯点强度参数作为本次载荷增量计算过程中的残余强度参数，重新进入步骤（1），开始新一轮迭代计算。如果已经收敛，则开始下一个载荷增量。在步骤（1）（6）的计算流程中，包含两重收敛判断，其一为非线性迭代过程中残余力的收敛，其二为软化过程中强度参数的收敛。注意：在步骤（3）中，如果高斯点在上次载荷增量中没有进入塑性状态，而本次加载过程中开始屈服时，对应于本次载荷增量的应力增量k(1)dep0d DD，其中由1p(kF)0D确定。对于理想脆塑性分析，所有的高斯点只有一个峰值强度参数和一个残余强度参数，步骤（6）中不需要更新强度参数，直接进入下一次载荷增量计算。有关脆塑性理论及其分析方法的细节可参阅文献39。4 应变软化过程模拟 岩土材料的软化过程分析及其有限元方法的实现最早见之于文献42，Lo42在考虑土体应变软化对边坡稳定性影响时，首先采用有限元方法对坡体进行了分析，然后将得到的应力分布作为极限平衡方法计算潜在滑面安全系数的已知条件。在进行应力分析过程中，为使问题可解，强度弱化段采用让材料先产生弹性变形增量，再将超出屈服面外的应力消除的方法。计算分析结果显示，由于应变软化土体（如超固结土）出现峰后强度弱化，导致边坡稳定性安全系数发生很大的变化（出现较大的下降）。受文献39,42的启发，本文基于脆塑性分析方法，将应变软化过程看作为一系列的脆塑性过程，并由此提出了一套求解应变软化问题的新方法。4.1 应变软化过程简化应变软化过程简化43 图 2 所示曲线 OPBC 为一个典型的包含应变软化过程的应力-应变曲线。岩（石）体从峰值强度 P点经历应变软化过程的 A2、A4和 A6点直至到达残余强度点 B。PA2段与 A2A4段等表示具有负切线模3r 2r 1r 0 0 3p 2p 1p 第 3 期 王水林等：应变软化岩体分析原理及其应用 613 量的线性软化段，下面图示说明将线性软化段简化为脆性跌落与塑性流动的过程。以材料从点 P 到点 A2的软化过程为例，将强度随应变的软化简化为图中的应力脆性跌落 PA1段与塑性流动 A1A2段，那么强度随塑性应变的增加而降低的 PA2段的应变软化过程就转化为脆塑性过程中的脆性跌落 PA1段与塑性流动 A1A2段。依此类推，应力从峰值点 P 到残余强度点 B 的软化过程就转化为一系列的脆性与塑性交互发生的过程。而弹脆塑性过程的数值求解已经有了较为成熟的方法39。由此简化过程，应变软化问题的求解转化为求解一系列脆塑性问题，避免了直接采用经典的弹塑性理论分析应变软化过程出现的峰值后区切线刚度矩阵的不定性问题。图图 2 应变软化段应力应变软化段应力-应变曲线的简化应变曲线的简化 Fig.2 Simplfication of stress-strain curve in strain-softening process 4.2 峰值后区强度参数弱化峰值后区强度参数弱化 在应变软化过程中，随着软化参数值的变化，岩体介质的总体强度逐渐下降，强度参数在峰值后区弱化。本文定义软化参数为塑性剪切应变，即pp13。为简单起见，假设强度参数是塑性剪切应变的双线性函数，*ppr*r()/0()（25）式中：代表屈服准则中强度参数（c,）或塑性势函数中的剪胀角（），下标“p”与“r”分别表示其为峰值和残余值；*为极限塑性剪应变，当塑性剪切应变大于其极限值后，材料进入残余强度状态，相关参数取残余值。塑性剪切应变极限值*应该由试验确定。c、和对应的塑性剪切应变极限值可以是不同的。4.3 岩土工程中应变软化问题求解过程岩土工程中应变软化问题求解过程 将应变软化过程简化为脆性跌落与塑性流动后，应变软化问题的求解就归结为一系列的脆塑性问题求解。不同之处在于脆塑性分析方法中高斯点应力进入屈服状态后直接从峰值强度屈服面跌落至残余强度屈服面，每个高斯点只有一次应力跌落计算。而在应变软化问题的求解过程中，在每个载荷增量步都需要计算应力跌落增量，直至高斯点的应力进入残余强度屈服面。对每一个处于屈服状态的高斯点而言，每次载荷增量结束后材料强度参数值将作为下一次载荷增量对应的脆塑性计算的峰值强度值。总体而言，简化为一系列脆性跌落与塑性流动过程的应变软化问题，求解思路与脆塑性问题求解思路基本相同，但过程较为复杂。图 3 为求解应变软化问题的计算流程图，程序中考虑了岩土工程开挖作业过程的模拟。图图 3 应变软化求解流程应变软化求解流程 Fig.3 Numerical procedures for dealing with strain-softening behavior 依据上述应变软化计算分析流程图，笔者编制了相应的适合求解软化岩土工程问题的有限元程序FEC_SSA（Finite Element Code for Strain Softening Analysis），该程序可以模拟分步开挖和强度参数不同演化行为的力学过程（如黏聚力弱化与摩擦系数 强化的所谓 CWFS 模型对应的应变软化过程）。目前嵌入的强度准则为 Mohr-Coulomb 准则，对应的塑性流动法则既可以是关联流动法则，也可以是非关联流动法则。程序运行过程中，载荷增量过程与强度参数更O P A2 A4 A6 B C A1 A3 A5 A7 开始 输入几何尺寸、边界条件、材料属性和外荷载数据 进行开挖计算并计算开挖力 根据指定的荷载因子施加增量载荷 求解控制方程 本构积分计算 计算残余力 检查迭代过程是否收敛 更新残余强度参数 检查两个强度更新步之间的 求解过程是否收敛 输出本次荷载增量的结果 结束 开挖循环 载荷增量循环 强度参数更新循环 迭代循环 是 是 否 否 614 岩 土 力 学 2014 年 新是同步进行的，对应的载荷增量步计算收敛后，记录每一个高斯点更新的强度参数并作为下一步载荷增量的初始强度。与常规弹塑性有限元分析相同，载荷增量步越小（对应的计算量越大），模拟应变软化过程的精度越高。5 数值算例分析 5.1 数值算例数值算例1均匀地应力场中应变软化岩体均匀地应力场中应变软化岩体内圆形洞室开挖问题内圆形洞室开挖问题 在文献40中，笔者介绍了应变软化分析方法的有限元实现过程，并通过可以获取半解析解的经典考题，检验了程序的正确性。本文在此基础上，对文献40中圆形隧道问题进行了进一步的分析，给出了不同强度弱化速率对围岩塑形区分布的影响，希望了解应变软化导致的变形局部化行为及其塑性区分布。图 4 为圆形洞室开挖模型（考虑到对称性，只取 1/4 作为研究对象）。表 1 给出了求解圆形洞室问题采用的几何与力学参数，岩体材料遵守Morh-Coulomb 屈服准则。表中塑性剪切应变控制参数*越小，材料在塑性变形过程中强度由峰值下降至残余值的“速度”越快，即应变软化速率越大；*越大，材料应变软化速率越小。理论上而言，脆塑性材料的应变软化速率为无穷大，理想弹塑性材料的应变软化速率为 0。事实上，当*取值比材料单轴压缩时峰值强度对应的轴向应变小约 3 个数量级时，就可以认为研究对象为脆塑性材料，当*取值比材料单轴压缩时峰值强度对应的轴向应变大约3 个数量级时，其本构模型计算得到的结果就与理想弹塑性模型的相同。图 5 给出了不同软化参数控制值（*=0.810-5）对应的洞室围岩塑性区形态及等效塑性半径的大小，从图中可以看出，当*逐渐变小，强度弱化速 率变大，应变软化变得明显，同时伴随变形局部化出现。虽然外载荷是轴对称的，但围岩弹塑形区的交界面并不是一个近似的圆弧（弹塑性解得到的则是一个近似的圆弧），这里用等效塑性半径大小描述塑形区的大小（等效塑性半径与开挖洞室的半径构成的圆环面积等于图中实际的塑形区面积）。(a)圆形洞室开挖模型几何尺寸 (b)圆形洞室附近局部网格划分图(11 172 个三角形单元,5 708 个节点)图图 4 圆形洞室开挖模型尺寸及有限元网格划分圆形洞室开挖模型尺寸及有限元网格划分 Fig.4 Geometry and finite element meshes of the circular tunnel 表表 1 几何与力学参数几何与力学参数 Table 1 Geometric and mechanical parameters r0/m 0/MPa E/GPa *cp/MPa cr/MPa p/()r/()p/()r/()3 5 10 0.25 0.8110-5 1.0 0.7 30 22 3.75 3.75 注：r0为圆形洞室半径；0为初始地应力；E 为弹性模量；为泊松比。(a)*=0.8 (b)*=0.008 (c)*=0.004 (d)*=0.002 (e)*=0.001 (f)*=10-5 图图 5 围岩塑形区形态与等效塑性半径围岩塑形区形态与等效塑性半径 Fig.5 Shapes of plastic zone in surrounding rock mass and equivalent plastic radii rp=4.18 m 塑性区 塑性区 rp=4.25 m rp=4.47 m rp=4.73 m rp=5.14 m rp=4.92 m 塑性区 塑性区 塑性区 塑性区 r0=3 m 50 m 50 m 第 3 期 王水林等：应变软化岩体分析原理及其应用 615 当*=0.8 时，计算得到的塑性区半径与表 1 中给出的峰值强度对应的理想弹塑性解相同。可以认为此时强度弱化速率很小，软化问题蜕化为理想弹塑性问题。当*=10-5，计算得到的塑性区半径即为表 1 中对应参数的脆塑性解，可以认为此时强度弱化速率非常快，单元内任一高斯点的应力进入屈服状态时，只要变形进一步增加，材料的强度即从峰值跌落至残余值。图6为表1中力学参数在*取值不同时对应的应力-应变曲线，图中只给出了围压为 0 MPa 与围 压为 2 MPa 时的曲线，*越小，应力-应变曲线的峰值后区越陡，对应于图 5 中变形局部化现象越明显。从图 6 中可以看出，不同的强度参数演化规律对应于不同的峰值后区应力-应变关系。反过来说，不同的全应力-应变曲线形态反映了不同岩土材料的强度演化特性。譬如，对某一种岩石（体）材料，如果通过试验获取的应力-应变关系曲线与图 6(c)中曲线形态比较一致，那么，在确定材料强度参数时，所考虑参数取值的演化规律应与对应的应力-应变关系一致。(a)*=0.8 (b)*=0.008 (c)*=0.004 (d)*=0.002 (e)*=0.001 (f)*=10-5 图图 6 不同强度弱化速率对应的应力不同强度弱化速率对应的应力-应变曲线应变曲线 Fig.6 Stress-strain curves corresponding to different strength-weakening parameters 表 2 给出了不同软化速率控制值（参数*）对应的等效塑性半径的有限元解，有限元解与半解析解的比较表明，两者基本吻合。在该算例中，图 5 与表 2 给出的塑性区形态与大小说明，因应变软化会导致变形局部化，但在复杂应力问题（包含平面或者三维问题）中，变形出现局部化后，应力自动调整进入新的平衡状态，最后得到的塑性区面积总体上与简化得到的半解析 解4344应该是相同的。5.2 数值算例数值算例2URL试验洞洞壁试验洞洞壁 V形缺口模拟形缺口模拟 加拿大原子能公司地下研究实验室（Underground Research Laboratory，简写为 URL）的试验洞是较 早开展岩石力学现场试验的实验室之一。许多学者以试验洞为对象，开展了岩石（体）变形与破坏机制、强度参数与演化规律、地应力分布、试验洞洞壁 V 形缺口模拟等研究，积累了比较详实的岩石力学试验与观察资料4547。URL 试验洞模型可以简化为一个拟平面应变问题，如图 7 所示。地应力大小在洞室横截面上的主应力为1和3，洞轴方向上的主应力为2。应力符号规定压应力为正。模型几何参数与初始地应力大小见表 3，文献 46采用了表中的 1#参数，而文献47采用的是表中 的 2#参数。0 2 4 6 8 10 00.001 0.002 0.003 0.004 0.005 1 1-3/MPa 0 MPa 2 MPa 0 2 4 6 8 10 0 0.001 0.002 0.003 0.004 0.0051 1-3/MPa 0 MPa 2 MPa 0 2 4 6 8 10 0 0.001 0.002 0.003 0.004 0.005 1 1-3/MPa 0 MPa 2 MPa 0246810 00.001 0.002 0.003 0.004 0.0051 1-3/MPa 0 MPa 2 MPa 0 0.001 0.002 0.003 0.004 0.005 10 8 6 4 2 0 1 1-3/MPa 0 MPa 2 MPa 0 0.001 0.002 0.003 0.004 0.005 1 0 MPa 2 MPa 024681-3/MPa 10 616 岩 土 力 学 2014 年 表表 2 不同软化参数控制值对应的塑性区半径对比不同软化参数控制值对应的塑性区半径对比 Table 2 Comparison of plastic radii for different softening parameters*pr/m pr*/m 0.800 4.18 4.17 0.008 4.25 4.27 0.004 4.46 4.40 0.002 4.73 4.74 0.001 4.92 5.04 10-5 5.14 5.22*=对应理想弹塑性解，pr=4.18 m；*=0 对应理想脆塑解，pr=5.14 m；pr为等效塑性半径（有限元解）；pr*为塑性半径（半解析解4344）。图图7 URL试验洞模型试验洞模型 Fig.7 A circular tunnel mode in URL 表表3 URL试验洞模型几何尺寸与力学参数试验洞模型几何尺寸与力学参数 Table 3 Geometric and mechanical parameters for URL 编号 r0/m 1/MPa 2/MPa 3/MPa E/GPa cp/MPa p/()p/()cr/MPa r/()r/()c*t/MPa 1#1.75 60 45 11 60 0.2 50 0 30 15 48 30 0.002 0.005 0.002 10.0 2#1.75 60 43 11 60 0.2 35 22 0 0.000 1 50 0 0.003 0.003 0.003 7.5 注：r0为圆形洞室半径；E 为弹性模量；为泊松比；t为拉应力极限值。Diederichs47和 Hajiabdolmajid 等46分析了洞室开挖完成后围岩表面形成的 V 形凹口，模拟结果如图 8、9 所示，他们采用的方法均为 FLAC，在强度峰值后区，均采用了黏聚力弱化与内摩擦角强化模型，但两者选择的峰值强度与残余强度参数的大小有较大差异。他们的模拟结果（塑性区大小）与现场的实际情况（见图 10）基本上是吻合的。笔者采用相同的力学参数，通过本文提出的求解应变软化问题的有限元程序，对试验洞开挖后形成的塑性区进行了模拟，图 11 为有限元网格，研究区域大小取 18 m18 m，有限元网格在洞室表面附近的的尺寸为0/40r，图 11(a)中模型共有四节点单元 13 505 个，节点数 13 616 个，图 11(b)为圆形洞室局部网格放大图。按照拟平面应变问题求解，计算过程中首先施加表 3 中给出的初始地应力（包括垂直于洞室截面方向的2），洞室一次开挖完成，开挖力分 5 次卸载直至开挖面应力为 0。得到的围岩剪破坏与拉破坏区如图 12 所示。采用 Diederichs47的参数（表 1 中 2#参数）计算得到的塑性区范围比现场观测到的塑性区小，采用 Hajiabdolmajid46的参数（表 1 中 1#参数）计算得到的塑性区范围与现场观测到的塑性区比较一致，而且与文献45给出的结果基本吻合。图图 8 Diederichs47采用采用 FLAC3D模拟模拟 URL 试验洞得到的围岩塑性区试验洞得到的围岩塑性区 Fig.8 Yielding zone around test tunnel at URL modeled by FLAC3D(Diederichs47)图图 9 Hajiabdolmajid46采用采用 FLAC2D模拟模拟 URL 的的 试验洞得到的围岩塑性区试验洞得到的围岩塑性区 Fig.9 Prediction of yielding zone around test tunnel at URL modeled by FLAC2D(Hajiabdolmajid46)r0 1 3 1 3 r=1.35a 破坏区轮廓 剪切破坏单元 拉破坏单元 第 3 期 王水林等：应变软化岩体分析原理及其应用 617 (a)照片47 (b)塑性区在断面上的分布46 图图10 URL试验洞围岩破坏区试验洞围岩破坏区 Fig.10 Failed zone around test tunnel at URL (a)有限元网格 (b)洞室区域局部网格放大图 图图11 有限元网格有限元网格 Fig.11 Finite element meshes (a)采用 Diederichs47的参数 (b)采用Hajiabdolmajid等46的参数 图图 12 有限元模拟得到的有限元模拟得到的 URL 试验洞围岩塑性区分布试验洞围岩塑性区分布 Fig.12 Predictions of plastic zone obtained by finite element simulation 为了了解图 12 中塑性区大小的差异，笔者根据 Diederichs47和 Hajiabdolmajid 等46采用的岩体力学参数值，采用变形控制方法，模拟了材料单元受三轴压缩的过程，给出了两组参数对应的应力-应 变 关 系 曲 线，如 图 13 所 示。图 13(a)为 Diederichs47采用的强度参数对应的应力-应变关系曲线，从图中可以看出，随围压的增加，对应于峰值强度的偏应力（13）增大，材料进入初始屈服时第一主应力将随围压的增大而明显增大，此外，随围压的增加，材料的力学行为从应变软化逐渐向强度硬化转化；图 13(b)为 Hajiabdolmajid 等46采用的强度参数对应的应力-应变曲线，由于材料初始摩擦角为 0，峰值强度对应的偏应力（13）与围压大小无关，材料进入初始屈服的应力不受围压影响，在进入非线性变形阶段后，材料首先呈现应变软化行为，随后材料的综合强度逐渐增加，随着围压变大，强度增加幅度增大。从强度参数演化对应的应力-应变曲线来看，图 13(a)显示的应力-应变曲线形态与一般岩石力学试验记录的应力-应变曲线形态比较相似，图 13(b)的应力-应变曲线形态在试验中几乎很少能观察到。通过比较图 13 中的两组参数对应的应力-应变曲3 1 1.3R 1.75 m 微震点 破坏区 1 3 张拉区 rf=1.2 a a=1.75 m 声发射点 0.26 剪切破坏区 r0=1.75 m 拉破坏区 剪切破坏区 r0=1.75 m 0.53 拉破坏区 618 岩 土 力 学 2014 年 线，可以看出，Diederichs47给出的岩体强度值和应变软化参数更加合理一些。但由于强度参数值略微偏高，在相同的围压下，采用 Diederichs47给出的岩体参数时，需要施加更大的压力才能使岩石进入屈服状态。(a)Diederichs47采用的强度参数 (b)Hajiabdolmajid 等46采用的强度参数 图图 13 Diederichs47和和 Hajiabdolmajid 等等46采用的采用的 强强度参数对应的应力度参数对应的应力-应变曲线应变曲线 Fig.13 Stress-strain curves corresponding to strength parameters by Diederichs47 and Hajiabdolmajid et al46 进一步检查两组计算参数得到的围岩内应力分布情况，发现离洞室开挖面很小的深度范围内单元高斯点最小主应力（第三主应力）随开挖面的距离增加而迅速增大，最大主应力的分布大体相同。虽然单元进入塑性状态（即应变软化状态）后塑性区