卸荷岩体有限元分析的基本方法

卸荷岩体有限元分析的基本方法李建林三峡大学土木水电学院2004年11月1引言自然界中的岩体是经过了各种各样的成岩运动及地质构造作用后的产物，岩体在漫长的地质年代构造作用过程中，岩体不断地受到加荷与卸荷的反复作用，从而形成了大小不等、方向各异、性质不同的各类结构面。即现今岩体是经过了各种各样作用后的损伤岩体，其存在岩体中的应力是残余应力（如图1所示）。图中c点是经过构造运动作用后，岩体的现今应力00,00状态，其值为（），为残留在岩体中的地应力，属岩体多次作用后的残余应力，为历次构造作用后达到的总变形。这种变形在历次作用后已损失掉，无法考虑。如果岩体开挖卸荷，则是残余应力进一步释放的过程。此时，曲线从c点沿cd向下移现今岩体初始应力状态C点Oσ'b(如边坡开挖卸荷)按卸荷实验曲线c()eO1df(如基础岩体受荷)先卸荷后加荷εhε,0σg0'c历史作用b按常规加荷实验曲线图1构造运动作用示意图应力应变图1构造运动示意图动，如果岩体卸荷为零，则曲线达到d点；如果继续卸荷使之产生移动，则曲线继续向下延伸至e点，岩体拉断，即达到残余抗拉强度。如果岩体受力曲线沿c点向上移动，则表现为加荷，如果加荷至b点时，达到了岩体的残余抗压强度，曲线延伸至c点破坏。如果岩体先开挖卸荷，然后再加载，则其受力过程为cfgh，可见，卸荷岩体力学与常规岩体力学相比，具有如下特点①应力路径不同；②力学参数不同；③屈服条件不同；④分析方法不同。因此，卸荷作用后岩体的本构关系、力学参数、破坏机理等与常规的分析方法是不同的。对于岩体的卸荷分析，应依据卸荷岩体力学方法来进行，才能使之与岩体的受力特性一致。2岩体卸荷分析的本构关系由于岩体开挖卸荷展现出明显的各向异性特性，并且岩体的卸荷应力应变关系一般为非线性，为了有效地描述各向异性非线性岩体的力学行为，就必须研究各向异性非线性本构关系。对于正交异性非线性岩体，其本构关系（平面应变）可写成如下形式3岩体卸荷应力应变关系岩体开挖卸荷是从σ0（初始地应力）开始的，如果岩体卸荷后产生拉应力，则拉应力最大为Rt（岩体抗拉强度），即卸荷应力△σ从0开始到σ0+Rt，根据岩体卸荷试验，典型的卸荷曲线如图2所示。图2典型的卸荷曲线4卸荷岩体力学参数分析与选择4.1基本方法研究以现场地质调查、材料试验为基础，与数值计算相结合，根据岩体结构面的分布状况、级别及岩体的不同性质，对岩体力学参数采用分区分级计算，最终对整个工程岩体力学参数进行评价。数值模拟的分级方法如图3所示。由图3，分级计算步骤如下：①首先根据工程地质勘探成果对给定的整个工程岩体各个区域分区分块。分区分块主要考虑岩石的类型及强度、岩石的风化程度、岩体的地应力分布情况、结构面的走向和倾向、结构面的间距及连通率、结构面的力学参数、地下水情况等有效采用数值模拟的计算试块，该计算试块应该能较够全面、真实地反映其对应岩应岩体的岩性组合和岩体结构特征。②根据每个区块的边界条件和外荷条件，对各区块内的分级计算试块逐级进行大三轴试验数值模拟，对每一区块的岩体参数进行评定。计算中应模拟区块内每一条裂隙、每一个层面。对同一区块不同级别计算试块计算顺序是先小尺寸后大尺寸，先小裂隙后大裂隙。第一级试块内为岩石和小节理裂隙材料，有关的计算参数和本构关系直接采用室内试验成果，通过第一级计算可得到第一级试块材料的力学参数。第二级计算试块内含第一级计算试块材料和中等裂隙材料。第一级试块材料的力学材料由前面计算所得；根据以上方法，依次计算各区块内不同级别的试块，后一级别的计算试块的材料力学性能采用前一级别试块的计算值，直至计算到该区块最大级别试块，最大区块的计算所的的参数就为该区块内的力学参数值。③把各区块和大断层及软弱厚层组合在一起，建立一整体计算试块，对整体计算试块的计算，可得到该区块整体的力学参数。图3参数数值模拟框图通过以上计算分析，计算顺序是先小裂隙、后大裂隙、先小尺寸、后大尺寸、先计算出来的参数作为下一级计算参数。基本思路如图4所示。由此计算，最终可得到该区块岩体的力学参数。裂隙-4，宽度8cm,间距10m裂隙-3，宽度4cm,间距4m裂隙-2，宽度1cm,间距1m裂隙-1，宽度1cm,间距1m裂隙-5，宽度10cm,间距15m123abcdfe图4分级计算示意图4.2岩体的等效变形参数计算试件在给定的应力增量条件下，通过弹塑性有限元分析计算，可计算各方向上的应变增量，由此计算其相当某单一应力状态下的材料的等效变形参数。在主应力空间，弹性应力状态下，应力应变关系服从广义虎克定律：平面应变状态下，计算试件在应力增量、的作用下，通过有限元计算得到该试件的，，未知量E，和通过三个方程可求解泊松比和变形模量，如下式所示。0312123由所求的变形参数，可得到单向应力状态下的应变增量和应力增量之间的关系式：E5卸荷岩体分析的基本方法岩体开挖卸荷，不同区域的岩体其卸荷程度不同，并随着开挖的进行，卸荷区域和卸荷程度不断变化。因此，卸荷分析的模拟过程也是不断调整的。应根据开挖进程和不同卸荷区域及程度取与之对应的卸荷岩体力学参数，这种模拟方法可通过图5进行模拟。图5计算分析图卸荷区扩展否将岩体视为连续介质或各向异性体计算初始应力场选择初始岩体参数计算开挖边界释放荷载计算开挖释放应力场组合、比较、划分卸荷区分区选择岩体参数计算开挖卸荷分区应力应变场开挖卸荷应力应变场停止是否6工程实例现以三峡工程永久船闸高边坡为例进行卸荷岩体的稳定分析。6.1计算条件6.1.1地应力船闸区内宏观构造应力场，最大主应力方向为NWW向而近EW向。船闸所处山体地区，上、下游方向被河湾切割，NWW向地质构造应力在地表已被释放，因此，在高程77m以上，该地区最大主应力方向以调整至NEE向（Y轴）。由于花岗岩的风化作用，地表约30~40m深度内的岩层，经风化后卸荷，在此地层内，仅为自重应力场。因此，根据工程地质研究资料，各方向地应力取值如图6所示。图6地应力取值图6.1.2力学参数1)初始弹性模量及初始变形模量。初始弹性模量及初始变形模量如表1所示。表1各向异性弹性模量和变形模量2）卸荷阶段的变形模量。岩体卸荷各阶段的变形模量根据前面分区块计算所得的结果并利用表2选用。表2边坡岩体卸荷后变形模量减低百分比表3岩体强度参数3）岩体强度参数。卸荷岩体强度参数如表3所示。6.1.3计算剖面及方案（1）计算剖面选取较为典型的三闸首16号剖面为非线性分析剖面，计算域取南北长1200m，底部高程–100m，如图7所示。图7三峡工程永久船闸高边坡16号剖面断面图（据设计报告）2)地应力施加方法首先将底部边界的垂直及右侧边界的水平方向固定，有限元计算时地应力的施加方法如下：将图6中地应力分布值等效移置为左边界节点，在重力及节点力作用下计算初始位移场及应力场；将上述计算结果得到的左边界固定（即为地应力作用下完成的初始位移），在此条件下模拟边坡开挖。3)计算方案计算中考虑三种方案如下：①常规弹塑性方案——方案一；②正交异性弹性方案——方案二；③非线性正交异性方案——方案三，抗拉强度取Rt=1.5MPa；④非线性正交异性卸荷方案—方案四，此时，抗拉强度取Rt=1.5MPa。6.2计算结果及分析1)计算位移场分析根据上述各方案的计算，最大位移出现在闸室侧墙的顶部，各方案的最大位移值如表4所示。方案四的位移等值线图如图8所示。表4各方案最大位移量2（a）Y方向水平位移（b）垂直位移图816号剖面方案4相对位移等值线图（cm）计算结果表明，常规弹塑计算方法的最大位移为3cm左右，当按正交异性弹性计算时，最大位移达到了6.8cm，再考虑非线性时，最大位移值达到了13cm，进一步考虑卸荷时，最大位移值达到了16cm。对于卸荷岩体，若考虑各向异性，计算位移值比常规大一倍，在考虑非线性时，大二倍，进一步考虑卸荷时，大四倍。因此，对于卸荷岩体的计算，当然不能按常规计算方法进行。2)计算应力场及分析永久船闸区岩体初始应力状态，大主应力在近地表处接近铅直，深部接近水平，小主应力方向近地表处接近水平，深部接近垂直，船闸底板（+112m高程）处最大主应力值约为8MPa左右，小主应力值约为3MPa。随着边坡的开挖下降，岩体的应力状态随之不断调整，局部区域由压应力状态变为拉应力状态，中墩上部由双向应力状态趋向双向拉应力状态、单向拉应力状态或零应力状态，在左右闸室底板与直立墙拐角处为压应力集中区，已处于压剪破坏状态。正交异性计算结果表明，在直立墙以上的斜坡岩体中均有拉应力区存在，但只是局部的，在微新岩体中一般在15m以内，而在微风化及岩基中则表现为拉剪破坏状态，且范围较大。在直立墙中，其上部约1／4范围为拉应力集中区，量值在0.6MPa左右，在离墙20余米的墙体中，有呈三角形的拉应力区，一直至墙体下方，该区内外均为压应力。中墩上部1／4范围均为拉应力区，拉应力值一般在1MPa以内，局部达3.5MPa左右，计算16剖面最大及最小主应力等值线如图9所示，由图可见拉应力的分布规律，此处不再叙述。（a）大主应力（b）小主应力图916号剖面主应力等值线图（MPa）考虑岩体拉坏时，在各种抗拉强度下，超过强度的拉应力点拉裂，拉应力区向邻近扩展，且表现成为岩体位移的增大。由上分析可见，船闸边坡的开挖卸荷发育过程，是应力状态调整及演变的根本原因，流变过程的发育使边坡岩体抗拉强度发生变化，岩体强度不断劣化，进而使边坡的变形量增加。因此，边坡的变形量的控制是岩体稳定研究和控制的关键。6.3与实测结果比较随着永久船闸开挖地进行，地质条件不断揭露，为准确而及时地校核对岩体的正确认识程度，检验边坡工程设计和施工工艺的合理性，监测岩体的变形是必要的，这也是岩体力学系统工程方法论很重要的一个工作环节，所测得变形量是在边坡稳定性研究过程中岩体主要的控制参变量。为此，根据永久船闸边坡有限的变形观测成果进行分析，并与计算成果进行比较。边坡观测点是随着边坡的不断开挖而逐渐设立的，如200m高程的TP42GP02观测点，于1995年11月份设立并开始观测的，此时的开挖高程由250m下降至170m高程，因此，1995年11以前的变形损失掉了，没有获得观测成果，该观测点观测到的变形是开挖到170m以下的影响变形，以及岩体的蠕变变形，也就是说，各观测点测得位移没有考虑的该点瞬时变形和部分蠕变变形。这说明对每一观测点的观测结果，我们必须分析其设立时间以及边坡的开挖的开挖高程等条件对应地进行分析比较，而不能笼统地拿来就比较。为了便于对比分析，我们以观测点TP40GP02、TP41GP02和TP42GP02所在剖面进行分析比较（如图10所示）图10测点布置及开挖观测情况TP41GP02、TP42GP02于1995年11月开始观测，此时的开挖高程为170m，TO40GP02于1996年4月开始观测，此时高程155m，至1996年11月（开挖高程+140m）的各观测点观测水平位移如下：TP40GP0222.88mm；TP41GP0224.53mm；TP42GP0217.38mm计算时，采用卸荷非线性岩体力学分析方法，其中考虑高程为170m、155m和140m这3个开挖步。考虑岩体抗拉强度Rt=1.5MPa。为了与观测结果比较，边坡计算开挖至+140m高程，计算成果如表5所示。表5计算成果与观测成果比较（挖至+140m高程）采用计算理论：卸荷非线性岩体力学理论；有限元有365个节点，572个单元，分8步开挖。计算表明，由卸荷岩体力学理论计算结果与实测比较值吻合。TP40GP02观测点，当抗拉强度Rt=1.5MPa时，按对应计算位移为18.2mm，与实测值很吻合。而总位移为73.3mm和76.6mm，由此可见，观测点设立前的变形发展也是比较大的。7结论根据上述分析可知，岩体卸荷力学特性与加载力学特性有重要区别，其区别在于岩体卸荷本构关系、卸荷力学参数、岩体的破坏机理及准则和分析方法等不同。因此，卸荷岩体中岩体本构关系、力学参数、岩体的破坏机理及准则和计算方法等是岩体力学研究的新内容，它必须有一