中国地质大学WH弹塑性力学作业

王宇等：有限元强度折减法在边坡可靠性研究中的应用 第 卷基于弧长控制法的边坡稳定性弹塑性有限元分析王宇 120090456（中国地质大学 武汉 430074）摘 要：通过有限元强度折减，使边坡达到破坏状态时，滑动面上的位移将产生突变，产生很大的且无限制的塑性流动，有限元程序无法从有限元方程组中找到一个既能满足静力平衡又能满足应力-应变关系和强度准则的解，此时不管是从力的收敛标准，还是从位移的收敛标准来判断有限元计算都不收敛，因此采用力和位移的收敛标准作为边坡破坏的判据是合理的。在强度折减过程中讨论了基本计算方法和弧长控制法的应用，关于有限元计算过程计算过程中怎样利用弧长控制法来调整强度参数的降低量，以便精确得到土体结构的极限状态及相应的强度参数值进行了研究。以杭兰高速公路巫山到奉节段大水田边坡为例，借助PLAXIS有限元分析程序，利用基于弧长控制法的弹塑性有限元法分析了该边坡的稳定性，并与传统的Bishop法、Janbu法等方法计算所得边坡稳定系数进行了对比分析。结果表明，有限元强度折减法能更加真实地反映边坡的实际情况，求得的边坡稳定系数更接近边坡的实际稳状态，显示出其在边坡稳定性分析中的一定优势。关 键 词：边坡稳定性分析；强度折减法；弧长控制法；稳定系数Slope stability analysis based on the arc length control elastic-plastic finite element methordWANG Yu 120090456 (School of Engineering, China University of Geosciences, Wuhan 430074 China )Abstract: With the c-tan reduction，the FEM model of slope reaches instability，and the value of the nodal displacement just after slope failure has a big jump compared with that before failure.This actually means that no stress distribution can be achieved to satisfy both the yield criterion and global equilibrium.Slope failure and numerical non-convergence take place at the same time.So non-convergence in finite element program can be taken as a suitable evaluation criterion of slope failure. In strength reduction process discussed the basic method and application of the arc length control, the finite element calculation process about the process of how to make use of the arc of the strength parameters to adjust the control method, in order to reduce quantity of soil structure accurately and the corresponding intensity limit state parameters were studied. In HangLan highway wushan to fengjie of large slope, for example, with the paddy PLAXIS finite element analysis program, based on the arc length of weak elasto-plastic finite element method (fem) control are analyzed, and the stability of the slope and the Bishop method, Janbu traditional methods of calculating coefficient of income of slope stability are analyzed. Results show that the finite element strength subtraction can reflect the actual situation of slope, the slope stability coefficient of slope stability is more close to actual condition, showing the slope stability analysis in the certain advantages.Key words: slope stability analysis; strenghen reduction method; arc length control method; stability factor 弹 塑 性 力 学 第 卷0 引言边坡稳定分析是边坡设计的前提，它决定着边坡是否失稳以及边坡失稳时存在多大推力，以便为支护结构设计提供科学依据。然而这个问题至今仍未得到妥善解决，因为解决这一问题必须先要查清坡体的地质状况及其强度参数，同时又要有科学合理的分析方法1。对于均质土坡，传统方法主要有：极限平衡法，极限分析法，滑移线场法等，就目前工程应用而言，主要还是极限平衡法，但需要事先知道滑动面位置和形状。极限平衡方法的一个局限是对较复杂的土层及土工结构，计算较困难。边坡稳定分析涉及复杂的地质地形边界条件、材料的应力-应变的非线性行为、初始地应力、水压力、地震荷载的耦合分析等等，多数情况下不能获得解析解。传统的数值分析方法一般只是得出边坡应力、位移、塑性区等，不能直接与边坡稳定建立定量关系。然而有限元强度折减法却能很好解决这一问题，本文用有限元强度折减法对边坡的稳定性进行了分析，给出用有限元计算此类问题的具体方法，特别要着重讨论在此计算过程中怎样利用弧长控制法来调整强度参数的降低量，以便精确得到结构的极限状态及相应的强度参数值。本文以国家重点高速公路杭兰线巫山奉节段YK39+300边坡为例，应用弧长控制技术，结合有限元强度强度折减对其进行稳定性分析，该方法简便易行，计算速度快,稳定性好，可以推广到其它工程中去。1 弧长控制法1.1 土的应力应变模式土工结构计算中广泛应用的理想弹塑性莫尔库仑模式，因为就土工结构的破坏分析来说该模式是相当简便而又可靠的。该模式的屈服面方程可以写为 （1）当主应力不按大小排列时，屈服面需用6 个与上式类似的方程来描述。在主应力空间中莫尔-库仑屈服面是一个六棱锥面。决定塑性流动方向的塑性势面方程和屈服面方程有相同的形式，但需用剪胀角来代替内摩擦角，以给出合理的塑性体积应变。在用初应力法进行弹塑性分析时，每次迭代首先得到试探应力为 （2）当 超出屈服面时,则用下式修正该应力使其回到屈服面: （3） （4）式中： 为塑性势面方程。这里采用Koiter一般流动法则, 由一致条件决定。将一致条件在点泰勒展开并略去高阶项,有 (5) 上式是 的线性代数方程组，由此可求解。具体计算时，为更为可靠有效,可在进行上述修正之前事先判定修正后的应力将位于屈服面的哪一区域。详细介绍见文献2。这种应力应变模式的积分方法为隐式方法，它与求解总体平衡方程的初应力法相结合，给出无条件稳定的算法。如引言所述，土工结构的安全系数应定义为其极限承载力与所需承载力之比。当采用莫尔库仑模式时稳定系数Fs可表示为 （6）这里下标“a”和“r”分别表示“所能提供的”和“所需要的”。为简化计算，这里假定c和按同一比例变化，从而Fs 也等于ca与cr之比，或tana 与tanr 之比。1.2 弧长控制法的应用为计算边坡的稳定系数，首先按所给土的力学性质参数( ca ,a 等) 进行计算，以施加土体的全部荷载。然后逐步减小土的强度参数直至土体结构破坏，从而求得cr 与r 。每次减小土强度参数后的迭代计算方法与一般弹塑性有限元类似，该迭代过程中典型的应力变化大致如图1中的折线ABC 所示。图1 基本迭代法中的应力变化Fig.1 stress variation in the basic iteration method图中L1，L2 分别对应于强度参数降低前后的屈服面。在每次减小强度参数后，对于超出新屈服面的应力(如图1中A 点应力)，将按上节所述方法修正到屈服面(如图1中B点应力) 。这样修正后的应力将不再符合平衡等条件，所以需进行迭代。在迭代过程中应力在新的屈服面上调整，以使整个区域的内力和变形符合连续体的平衡、几何和物理三方面的条件。这种迭代方法可避免不平衡误差的积累。具体的推导如下：有限元平衡方程： (7)式中是外部荷载增量。在迭代过程中有，则可得迭代方程： (8)在开始迭代时i=1，取为经过修正而符合当前屈服条件的应力。随着迭代次数的增加，应力的子增量将趋近于0。当用初应力法时，记，则上式可改写成： (9)当结构趋于破坏时，在强度参数不变的情况下变形持续发展。如分别以强度参数比和位移为纵横坐标画曲线,则该曲线趋于水平。如采用上述基本算法,当土的强度参数在某一步中减小过多，则结构不可能达到平衡，从而使迭代失败，不能求出准确的安全系数值。为解决此问题，这里采用在Riks方法基础上发展的一种弧长控制法34 。弧长法是计算结构破坏荷载及跟踪结构变形曲线软化段的一种方法。其基本原理是在迭代过程中调整荷载增量的大小，以得到收敛的解。荷载调整的原则是使位移向量增量的内积在迭代过程中保持为常量。但是，在安全系数的计算过程中如直接采用这种方法，则意味着在迭代过程中不断改变强度参数的减小量，这样屈服面随迭代过程不断变化。如此进行计算，收敛性将很差。为使计算更为有效，这里进行一些简化处理。如图2 所示,设L1、L2和L3分别代表初始屈服面、试探屈服面和迭代结束后的屈服面，我们将应力点沿AD的移动看成沿折线ABCD的移动。在迭代过程中屈服面用L2，但用弧长法所求出的荷载乘子m修正所求出的应力，使其近似对应于最终屈服面L3上的应力点D。图2 使用弧长控制法中的应用变化Fig.1 stress change baseed on arc length control这样式(8) 右端的应力此时可写为 (10)这里下标A，B，C 给出应力所对应的点，m是荷载乘子，在这里为材料强度参数降低量中所实现的部分。上式右端第二项是使用弧长法时对试探屈服面上应力的修正,此种修正可在迭代结束之后进行。由于在新屈服面上的应力点并非总是对应于原来屈服面上的某一点，式(9) 的修正也并非精确，故这里的处理方法是近似的。但如果我们注意到在结构接近破坏时，土的强度参数几乎不再变化，则不难理解上述处理方法的最终误差将是很小的。将式(9)代入式(8)，有(11)在上式的推导中利用了的条件。记，则式（10）可写为 (12)其中 (13) (14)这里我们采用一种简化的弧长法，它要求相邻两次迭代所求得的位移向量增量的内积为常数，即 (15)注意到，可有将式(11) 代入上式可以导出迭代过程中依据不平衡力的大小调整荷载乘子的公式为对于材料非线性问题，上述方法的一个重要修正是用部分结点的位移来代替其中的整体结点位移向量，以便更准确地探知结构破坏或软化过程中的应变集中5 。计算是分步进行的。在每一步中将强度参数降低之后，我们保持式(13) 中的l不变，而据不平衡力的大小按式(15)对强度参数降低量进行调整,当结构破坏时我们即可得到接近水平的强度参数比- 位移曲线。由于采用理想弹塑性本构模式,再利用上述的弧长法，一般总可得到该曲线的水平段，从而定出结构的安全系数(边坡工程中通常称其为稳定系数)。 (16)2 强度折减法及判据2.1 强度折减法大量试验表明：一个纯摩擦土坡在土自重增加的实验(离心机实验)中不会发生破坏。因此，对稳定系数更恰当定义是6：稳定系数这里，S代表抗剪强度。实际抗剪强度和计算获得的保证土体平衡状态所需要的最小剪应力之比是土力学中传统上使用的安全系数。通过引入标准库伦条件，安全系数（边坡工程中常称为稳定系数）可以表达为：稳定系数这里，c和是输入强度参数，而n 是实际正应力分量。cr 和r 是不断减小到恰好足够大而能保持土平衡的抗剪参数。上面描述的原理正是PLAXIS程序中为计算整体稳定系数而使用的折减phi/c方法的基础。应用这种方法，内聚力和内摩擦角的正切将成比例的减小： 稳定系数强度参数的减小将由总乘子Msf 来控制。程序在开始计算时默认Msf =1.0 ,然后Msf按设置的数值递增至计算模型发生破坏,此时的Msf 值即为计算模型的稳定系数值。有限元强度折减法不需要对滑动面形状和位置做假定，也无需进行条分，通过强度折减使边坡达到不稳定状态,非线性有限元静力计算将不收敛，此时的折减系数就是稳定安全系数。2.2 有限元中边坡破坏的判据有限元强度折减法分析边坡稳定性的一个关键问题是如何根据有限元计算结果来判别边坡是否处于破坏状态。目前的失稳判据主要有两类：(1)在有限元计算过程中采用力和位移的不收敛作为边坡失稳的标志710。(2)以广义塑性应变或者等效塑性应变从坡脚到坡顶贯通作为边坡破坏的标志1113。以上两种判据得到的稳定系数相差不大。3 实例分析3.1 边坡概况YK39+300边坡位于国家重点公路杭州至兰州线重庆巫山奉节段重庆市巫山县龙井乡六水村境内，起止桩号为YK39+105YK39+355，长250m，路堑边坡高度约10.831.5m。坡体地层从上到下依次为：残坡积含碎石亚粘Qel+dl土；强风化钙质泥岩；强风化泥质粉砂岩（T2b2）紫红色，层状构造,裂隙较发育，局部间夹灰绿色薄层粉砂岩，岩芯破随稳定性差，岩层产状33036；弱风化泥质粉砂岩（T2b2），紫红色，层状构造，裂隙较发育，局部间夹灰绿色薄层粉砂岩，岩芯呈柱状，岩石较完整。稳定性较好，岩层产状33036。钻孔揭露无地下水。据岩土工程勘察资料，该边坡后缘出现拉张裂缝，采用极限平衡法计算所得边坡稳定系数为0.9921.089，边坡处于稳定极限状态，提出了削坡工程处治方案初步分析认为，该边坡实施削坡工程后仍存在2条潜在滑移面，现拟对削坡后边坡稳定性进行评价(图1) 。 图1 YK39+300边坡工程地质剖面图Fig.1 Engineering geological cross - section No.3.2 物理力学参数选取为了准确、可靠地对边坡进行稳定性分析，从而为边坡治理提供科学的依据,科学合理的选取边坡物理力学参数至关重要。因此，在本次计算中，土体强度参数采用原状土样试验成果的平均值，并根据原地形条件下的边坡稳定情况，反求了土体的强度参数作对比分析，各层岩土体计算参数见表1。表1 岩土体物理力学参数Table1 Physical and mechanical parameters of soil地层岩性重度/kNm-3内聚力/Mpa内摩擦角/ ()弹性模量/Kpa泊松比/碎石土17.50.141930000.27强风化钙质泥岩190.212035000.31弱风化泥质粉砂岩2222370000.28弱风化泥质粉砂岩242.23490000.304 边坡模型建立4.1 力学模型的建立考虑到岩土体材料变形的非线性特征，本文采用理想弹塑性MohrCoulomb模型，模型的破坏包络线和MohrCoulomb强度准则(剪切屈服函数)以及拉破坏准则(拉屈服函数)相对应。4.2 边界条件设定根据该边坡实际状况，模型边界约束条件如下：(1)左右边界约束为法向约束边界，取u = 0，v 0 ( u为X 方向位移，v为Y方向位移)；(2)下部边界为全约束边界，取u = v = 0；(3)上部边界为自由边界。4.3 初始应力场由于该边坡上部为松散土体,下部为基岩，边坡变形和破坏主要发生于二者的交界面处，其构造应力在过去的地质历史时期已经释放。因此，模型建立时不考虑构造应力的作用，只考虑重力场的影响。初始应力场模拟如图2。图2 初始应力矢量图Fig.2 The initial stress vector diagram4.4 结果分析通过对土体抗剪强度指标的逐步折减，直到有限元计算不收敛，此时滑面上的位移产生突变，此时得到的Msf值即为边坡的稳定系数。图3为边坡最大剪应变图，图4为边坡的总位移增量图，从图中出现一个明显的最大剪应变率集中带，即为该边坡的潜在滑动面。图4给出随强度参数的降低(即计算稳定系数的增加)边坡坡脚处A点水平位移增加的Fs-d关系曲线。看到该曲线的末段几乎绝对水平，表明土的强度参数降低到该值时，位移持续增大，亦即土体趋于破坏。土原有强度参数值与此时的强度参数值之比Fs 即为稳定系数，削坡后边坡的稳定系数为1.4224。A潜在滑面 图3 边坡最大剪应变图 图4 边坡总位移增量图Fig.3 Contours of maximum shear strain of slope Fig.4 Total displacement incremental of slope图4 坡脚Fs-d（mm）关系曲线Fig.3 Fs - d (mm) curve of slope toe为便于将有限元强度折减法计算结果与极限平衡法计算结果进行对比，建立与PLAIXS相同的模型，采用Ordinary、Bishop和Janbu 3种方法对边坡的稳定系数进行计算，计算得到边坡稳定系数1.3731.412(表2)，比有限元强度折减法求得的边坡稳定系数1.4224偏小。造成二者差值的原因是由于传统边坡稳定性分析方法在计算过程中，没有考虑土体内部的应力、应变关系，无法分析边坡破坏的发生和发展过程，在求稳定系数时通常需要假定滑动面的位置和形状。而有限元法克服了这一不足，不仅满足平衡条件，而且还考虑土体应力、应变关系，能够得到边坡在荷载作用下的应力、应变分布，模拟出边坡破坏的实际滑动面，考虑了土体的变形协调关系，所以求得的稳定系数稍大一些，但符合边坡实际稳定状态。表2基于极限平衡理论方法计算所得边坡稳定系数Table 2Slope stability factors calculated with limit equilibrium method分析方法稳定系数Ordinary1.373Bishop1.389janbu1.4125 结论（1）为准确得到土体结构的极限状态及相应的强度参数，本文旨在研究了在此问题的计算中如何应用弧长控制法来解决这一问题，从而得到一个可靠而有效的计算方法。文中借助大型有限元分析软件PLAXIS，结合强度折减法对YK39+300段削坡后的边坡进行了稳定性分析，从而实现了弧长控制这一技术，应用于边坡稳定性分析，结果可靠合理。（2）与传统的边坡稳定性分析方法相比，有限元强度折减法不仅满足平衡条件，而且还考虑了土体的本构关系以及变形对应力的影响，求稳定系数时，不需要假定滑移面的形状，也无需进行条分，能够模拟出边坡破坏的实际滑动面，所求得的稳定系数更符合边坡的实际稳定状态，显示出它在边坡稳定性分析中的一定优势。参 考 文 献1 陈祖煜.土质边坡稳定性分析- 原理方法程序M .北京:中国水利水电出版社, 2003.2 Song E X , Vermeer P A , Bonnier P G. Implementation and Application of Mohr - Coulomb Model with Ten2sion Cut - off for Limit Load Analysis. In : Proc of Inter Conf on Comp Meth in Struct and Geotech Eng. HongKong , 1994, 830-835.3 Riks E. The Application of Newtons