有限元强度折减法

有限元强度折减法1背景1974年，Smith&Hobbs[1]使用有限元方法分析了φu=0条件下的边坡稳定性并与Taylar[2]的结果进行对比，得到了很好的一致性；1975年，Zienkiewicz等[3]考虑c’、φ’进行有限元边坡稳定性分析，其结果与圆弧滑面解有较好吻合；1980年Griffiths[4]验证了一系列具有不同材料特性和形状的边坡稳定性并通过与Bishop&Morgenstern[5]的结果进行了对比确定了数据的可靠性；此后也有研究证实了利用有限元方法进行边坡稳定性分析的可靠性[6,7,8,9]；在文献[9]中，引入一些案例证明了有限元强度折减法的准确性，并证明了有限元强度折减法在分析非均质边坡时相对于传统方法的优越性。2001年，郑颖人等[10]把有限元强度折减法引入国内，并对此进行了后续研究[11,12,13,14]。相较于一些传统的边坡稳定型分析方法，有限元强度折减法有以下几个优点[9]：(1)不必假设滑面的位置和形状，当土体自身强度不足以抵抗剪应力时土体失稳会自然发生。(2)由于有限元强度折减法中没有条分的概念，因此也不必假设条间力，在整体失稳之前土体都处于整体稳定状态。(3)使用有限元方法能够查看破坏过程。2有限元强度系数折减法1.模型参数边坡模型主要包括六个参数，分别是：膨胀角ψ、内摩擦角φ’、黏聚力c’、弹性模量E’、泊松比υ’、重度γ。膨胀角影响土体屈服后的体积变形，若ψ0，则土体屈服后体积减小，若ψ0则体积增大，ψ=0则体积不变。ψ=φ的情况被称之为关联流动法则，但是此时ψ值通常高于实验观测值，特别是在侧限条件下会提高土的承载力预测值。边坡稳定型问题通常是处于无侧限条件下，此时膨胀角的选取不再重要[9]，因此文献[9]选取ψ=0条件下的非关联流动法则，并且通过案例分析可以得出此膨胀角的选取可以得出准确的安全系数以及滑动面。c’和φ’指Mohr-Coulomb准则中边坡土体的有效黏聚力和内摩擦角；E’和υ’是土体材料的弹性参数，这两个参数对土体稳定性分析的影响较小；γ是土体的重度。应用有限元方法进行边坡稳定性分析中最重要的三个参数是c’、φ’、和γ。2.屈服条件(1)Mohr-Coulomb准则Mohr-Coulomb准则用大小主应力表示如式(1)所示：𝜎1′−𝜎3′2=𝜎1′+𝜎3′2𝑠𝑖𝑛𝜑′−𝑐𝑐𝑜𝑠𝜑′(1)其中，𝜎1′、𝜎3′分别指土中一点的大小主应力。在主应力空间中，如果不考虑𝜎1、𝜎2、𝜎3之间的大小关系，屈服面是一个不等角六棱锥，在π平面上是一个等边不等角六边形。(2)D-P准则(4)(6)D-P准则可以写成式(2)形式：−βI1+√J2=𝑘𝑓(2)其中I1为第一应力不变量、J2为第二偏应力不变量，β和kf为试验常数。在主应力空间中其屈服面为一个圆锥，在π平面上是一个圆形。3)D-P准则转换为Mohr-Coulomb准则首先引入参数b，如式(3)所示：b=𝜎2−𝜎3𝜎1−𝜎3(3)则，𝐼1和√𝐽2分别可转化为式(4)：𝐼1=3(𝜎1+𝜎3)2+(𝑏−12)(𝜎1−𝜎3)√𝐽2=𝜆√3(𝜎1−𝜎3)其中𝜆=√1−𝑏+𝑏2将其带入(2)，得式(5)：𝜎1−𝜎3=1.5𝛽2−𝑏𝛽+𝜆√3(𝜎1+𝜎3)+𝑘𝑓𝛽2−𝑏𝛽+𝜆√3(5)与式(1)对比可知两个准则之间的转换关系如式(6)所示：sinφ𝑏=1.5𝛽2−𝑏𝛽+𝜆√32c𝑏cosφ𝑏=𝑘𝑓𝛽2−𝑏𝛽+𝜆√3因此，当b=0时，即外角点外接DP圆的两个试验常数分别如式(7)所示，当b=1时，即内角点外接DP圆的两个试验常数分别如式(8)所示。β=2𝑠𝑖𝑛𝜑√3(3−𝑠𝑖𝑛𝜑)，𝑘𝑓=6𝑐𝑐𝑜𝑠𝜑√3(3−𝑠𝑖𝑛𝜑)(7)β=2𝑠𝑖𝑛𝜑√3(3+𝑠𝑖𝑛𝜑)，𝑘𝑓=6𝑐𝑐𝑜𝑠𝜑√3(3+𝑠𝑖𝑛𝜑)(8)σ1σ3σ2≥σ1σ3σ2≥≥σ3σ1σ2≥≥σ321≥≥231≥≥σ2σ1σ3≥≥σ1σ2σ3≥σ1σ2==b123==b0外角点外接DP圆内角点外接DP圆3.安全系数的定义(1)Mohr-Coulomb准则中的安全系数1955年，Bishop[15]首先在边坡稳定性分析中提出了抗剪强度折减的概念，在有限元强度折减法中通过将坡体的强度参数：黏聚力c和内摩擦角φ同时除一个折减系数Ft，得到一组新的c’和φ’值，作为一个新的强度参数输入进行试算，(10)当计算不收敛时，对应的Ft即为所求的安全系数，此时坡体达到极限状态，发生剪切破坏。c’=c/Ftφ’=arctan(tanφ/Ft)(2)D-P(Drucker-Prager)准则中的安全系数取Ft为D-P准则中的强度折减系数，则D-P准则可以表示为式(9)，−𝛽𝐹𝑡𝐼1+√𝐽2=𝑘𝑓𝐹𝑡(9)(3)不同屈服条件下安全系数转换[13]首先引入Mohr-Coulomb等面积圆屈服准则，在π平面上，其屈服面是一个圆，并且面积与Mohr-Coulomb准则的不等角六边形相等，Mohr-Coulomb等面积圆屈服准则中的试验参数如式(10)所示：−β=sinφ√3(√3𝑐𝑜𝑠𝜃𝛿−𝑠𝑖𝑛𝜃𝛿𝑠𝑖𝑛𝜑)𝑘𝑓=√3ccosφ√3𝑐𝑜𝑠𝜃𝛿−𝑠𝑖𝑛𝜃𝛿𝑠𝑖𝑛𝜑式中𝜃𝛿=𝑎𝑟𝑐𝑠𝑖𝑛−23𝐴𝑠𝑖𝑛𝜑+[49𝐴2𝑠𝑖𝑛2𝜑−4(𝑠𝑖𝑛2𝜑3+1)(𝐴23−1)]122(𝑠𝑖𝑛2𝜑3−1)，𝐴=√𝜋(9−𝑠𝑖𝑛2𝜑)6√3简称外接圆屈服准则为DP1准则，其试验常数分别为β1，kf1；Mohr-Coulomb等面积圆屈服准则为DP2准则，其试验常数分别为β2，kf2。把DP1准则表示为f1=√𝐽2=𝛽1𝐼1+𝑘𝑓1，DP2准则可表示为f2=√𝐽2=𝛽2𝐼1+𝑘𝑓2。令η=β1\β2=kf1\kf2=f(φ)，f1=𝛽1𝐼1+𝑘𝑓1=𝜂𝛽2𝐼1+𝜂𝑘𝑓2，所以f1f2=𝜂𝛽2𝐼1+𝜂𝑘𝑓2𝛽2𝐼1+𝑘𝑓2=𝜂=𝑓(𝜑)。由此可知，η是φ的函数，当φ取不同值时可以得到不同的η值如表1所列：表1不同内摩擦角时的η值4.失稳判据目前两个比较主流的失稳判据分别是有限元计算中力不平衡和位移的不收敛以及广义塑性应变或者等效塑性应变从坡脚到坡顶贯通。Griffiths[9]和郑颖人[11,12,13,14]都使用计算不收敛作为失稳判据。Griffiths[9]提出，当在用户定义的最大迭代数目下计算仍不收敛时，则没有任何一种应力分布方式可以同时满足Mohr-Coulomb准则以及整体稳定，这种情况可看做边坡失稳判据。边坡失稳与数值计算不收敛同时发生，并伴随着极大的节点位移，并以1000作为最大的迭代步数。郑颖人[14]提出，有限元的计算迭代过程就是寻找外力和内力达到平衡状态的过程，整个迭代过程直到一个合适的收敛标准得到满足才停止。可见，如果边坡失稳破坏，滑面上将产生没有限制的塑性变形，有限元程序无法从有限元方程组中找到一个既能满足静力平衡又能满足应力-应变关系和强度准则的解，此时不管是从力的收敛标准，还是从位移的收敛标准来判断有限元计算都不收敛。3案例分析例一，不含地基的均质边坡[9]该边坡如图1所示，有限元程序采用Mohr-Coulomb失效准则，建立平面应变条件下八节点四边形单元减缩积分计算模型，其强度参数为φ’=20°，c’/γh=0.05。边坡坡度为26.57°(2:1)，坡底水平，其边界条件为坡底约束竖直方向位移与水平方向位移，左侧约束水平方向位移，其余面为自由面。施加重力荷载后使安全系数从0.8到1.4逐步变化直至计算不收敛图1不含地基的均质边坡每一个安全系数对应的迭代次数如表2所列，当真正的安全系数接近时需要更多的迭代次数。表2例一计算结果当安全系数为1.4时，无量纲位移E’δmax/γH2突变，并且此时计算无法收敛，在此情况下有限元计算结果与Bishop&Morgenstern[5]给出的结果吻合良好，如图2所示。图2安全系数与无量纲位移边坡失稳时(FOS=1.4)节点位移矢量和网格变形如图3(a)和图3(b)所示，由此可得到边坡的潜在滑动面。图3安全系数为1.4计算不收敛时边坡变形(a)节点位移矢量(b)网格变形例二，有软弱层的不排水黏性土边坡在本案例中，使用Tresca准则(φu=0)进行总应力分析。边坡几何形状如图4所示，地基厚度与边坡高度相同，该边坡有一个软弱层，在有限元计算中，令其抗剪强度(Cu2)在一定范围内变化但其周围土体抗剪强度保持Cu1/γH=0.25不变。利用有限元方法计算该边坡的安全系数结果如图5所示，对于均质边坡情况，Cu2/Cu1=1，有限元计算结果与Taylor[2]的结论很接近，随着软弱层的强度逐渐减小，在Cu2/Cu1≈0.6时，结果发生了明显的变化。分别假定圆弧滑面和穿过软弱面的三段线滑面并利用Janbu法计算安全系数，可见在Cu2/Cu1≈0.6处也发生了滑动机制的转换，当Cu2/Cu10.6时，潜在滑面形状为圆弧，当Cu2/Cu10.6时，潜在滑面为结构软弱面。图6更加清晰的展示了这一现象，图6(a)为均质边坡(Cu2/Cu1=1)时的潜在滑面，可见此时的滑面形状为圆弧滑面，与Taylor[2]的预测相同；图6(c)为软弱层强度只有其周围土体20%(Cu2/Cu1=0.2)时的潜在滑面，此时潜在滑面沿软弱层发展；图6(b)为软弱层强度只有其周围土体60%(Cu2/Cu1=0.6)时的潜在滑面，此时圆弧滑面和沿软弱层的三段线式滑面都有可能发展，至少存在两种明显的滑动机制。图4有软弱层的不排水黏性边坡图5不同软弱层强度时的安全系数图6不同软弱层强度下的网格变形(a)Cu2/Cu1=1.0(b)Cu2/Cu1=0.6(c)Cu2/Cu1=0.2例三，不同坡度边坡安全系数计算[13]，验证Mohr-Coulomb等面积圆屈服准则均质边坡,坡高H=20m,土容重γ=25kN/m3,黏聚力c=42kPa,内摩擦角φ=17°,求坡角β分别为30°,35°,40°,45°,50°时边坡的安全系数。计算结果如表3所列。表3安全系数计算结果从表中计算结果可以看出,采用外接圆屈服准则计算的安全系数比传统的方法大许多,采用莫尔-库仑等面积圆屈服准则计算的结果与传统极限平衡方法(Spencer法)计算的结果十分接近,说明采用莫尔-库仑等面积圆屈服准则来代替莫尔-库仑不等角六边形屈服准则是可行的,这样使计算大为方便。而采用外接圆屈服准则计算的安全系数要比莫尔-库仑等面积圆屈服准则计算的结果大η(1.21)倍。例四，存在两组节理面的岩质边坡稳定性分析[12]如图7所示，岩体中存在两组方向不同的软弱结构面，贯通率100%，第一组软弱结构面倾角为30°，平均间距10m；第二组软弱结构面倾角75°，平均间距10m。岩体重度为25kN/m3，弹性模量1×1010Pa，泊松比0.2，黏聚力1MPa，内摩擦角38°，两组节理参数相同，重度为17kN/m3，弹性模量1×107Pa，泊松比0.3，黏聚力0.12MPa，内摩擦角24°。按照二维平面应变问题建立有限元模型，按照连续介质处理。通过有限元强度折减，求得坡体破坏时的运动矢量如图8所示，滑动面如图9(a)所示，它是最先贯通的塑性区，塑性区贯通并不等于破坏，当塑性区贯通后继续发展到一定程度，岩体发生整体破坏，同时出现第二条贯通的塑性面，如图9(b)所示。求得的稳定安全系数如表4所列，其中，极限平衡方法计算结果是根据最先贯通的那一条滑动面求得的。图7岩质边坡节理图9坡体破坏时的运动矢量图图8极限状态时的塑性区表4例四计算结果例五，存在接触问题的边坡稳定性分析[12]当边坡中存在如图10所示的硬性结构面时，不能按照例四中软弱结构面的方法进