土坡稳定的有限元塑性极限分析上限法研究

第23卷第11期岩石力学与 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 学报23(11)：1867～18732004年6月ChineseJournalofRockMechanicsandEngineeringJune，20042002年8月29日收到初稿，2002年11月8日收到修改稿。作者王均星简介：男，39岁，博士，1988年毕业于武汉水利电力大学水电系水工结构专业，现任副教授，主要从事水工结构方面的教学与科研工作。E-mail：junxing.wang@163.com。土坡稳定的有限元塑性极限分析上限法研究王均星1王汉辉2吴雅峰1(1武汉大学水利水电学院水工教研室武汉430072)(2长江水利委员会武汉430010)摘要系统地介绍了结构塑性极限分析的原理和方法。借助有限单元法和线性规划，运用塑性极限分析的上限法可求解土坡的极限承载力和安全系数，同时分析了该法的困难所在。以一经典土坡的稳定性作为算例，讨论了屈服条件近似对精度的影响和边坡土体容重对极限荷载的影响，论证了本方法的正确性。算例 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 明该方法具有较高的精度。关键词土力学，线性规划，塑性极限分析，上限法，安全系数分类号TU413.62文献标识码A文章编号1000-6915(2004)11-1867-07STABILITYANALYSISOFSOILSLOPEBYFINITEELEMENTMETHODWITHPLASTICLIMITUPPERBOUNDWangJunxing1，WangHanhui2，WuYafeng1(1SchoolofWaterResourcesandHydropower，WuhanUniversity，Wuhan430072China)(2ChangjiangWaterResourcesCommission，Wuhan430010China)AbstractThetheoryandmethodoftheplasticlimitanalysisareintroducedsystematically.Byusingthefiniteelementmethodandlinearprogram，theupperlimitmethodcanevaluatethebearingcapacityandthesafetyfactorofsoilslope.Theinfluenceofyieldsurfacelinearizationonthedegreeofaccuracyandtheinfluenceofslopesoildensityonthelimitloadarediscussed.Aclassicalslopecalculationcaseshowthatthepresentedmethodisofcorrectnessandhighprecision.Keywordssoilmechanics，linearprogram，plasticlimitanalysis，upperlimitmethod，safetyfactor1概述目前常用的边坡稳定分析方法是刚体极限平衡法和有限单元法。前者以其概念清晰、计算简便、工程实践资料丰富而被广泛采用，但未考虑材料的应力应变关系，不能给出相应的应力场和位移场；后者虽然能给出边坡的应力场和位移场，却无法给出明确的安全系数，在工程上多作为边坡稳定参考分析方法。事实上，在实际工程当中，人们关心的不是边坡失稳的整个过程，而是边坡开始滑动那一瞬间边坡的极限承载能力以及边坡的应力状态或滑动速度场，或者是边坡在正常工作下的安全度。而塑性极限分析法的最大优点就是回避了工程中最不容易弄清楚的本构关系，直接研究边坡的极限状态，因而是一种合理而且可行的方法。极限分析法从极大极小原理出发，运用上限解法和下限解法，放松极限荷载的某些约束条件，寻求问题的上限解和下限解。而问题的真正解答在上限解与下限解之间。岩土材料没有唯一的极限荷•1868•岩石力学与工程学报2004年载，借助于理想塑性材料的上下限定理，并考虑岩土材料的摩擦屈服特性及流动特性，可以推求极限荷载的近似值，由此得出的应力场为其下限解，所得的速度场为其上限解。本文在Sloan[1，2]提出的有限元塑性极限分析上限法的基础上，对速度不连续面上的约束条件加以改进，基于虚功率原理和有限元等效结点荷载的思想，建立以超载系数(极限荷载)和强度储备系数(安全系数)为目标函数的线性规划模型，借助线性规划求解边坡稳定问题。2塑性极限分析的基本原理2.1极限状态结构物内部构成的破坏结构的状态通常称为极限状态[3]，它是结构物在破坏与未破坏之间的临界状态，应该满足下述条件：(1)平衡条件；(2)边界条件；(3)屈服条件；(4)破坏机构条件。完全满足上述4个条件的解答称为极限状态的完全解或真解，只满足部分条件的解答称为极限状态的近似解，比如上限解。满足条件(2)，(3)，(4)的结构速度场为机动容许速度场，它是极限状态下速度场的上限解。2.2极限分析上限原理结构物有多种破坏形式，每个机动容许速度场对应着一个外荷载。因为结构已经破坏，因而这个荷载一定比极限荷载大。上限定理可以表述为：在所有的与机动容许速度场对应的荷载中，最小的荷载为极限荷载，与最小荷载对应的机动容许速度场为破坏速度场。2.3基本假定为了简化计算，采用如下假定[3]：(1)材料是理想刚塑性的，因此，内功功率不含弹性变形能；(2)结构物即使临近破坏状态时，变形依然很小，因而虚功率原理成立；(3)参与超载的外荷载按同一比例增加，以确定超载系数；(4)各种材料强度参数按同一比例降低，以确定材料强度储备系数。3上限法的基本方程3.1单元的离散本文仅讨论二维问题，土体用三角形单元离散，不同的单元可赋予不同的力学参数以反映其不均匀性，但同一单元的材料均一。3.2三角形单元的塑性变形速度表示假定三角形单元的塑性变形速度(以下简称速度)在单元内部呈线性分布，则单元内任一点的速度分量可表示为节点速度分量的线性函数：⎪⎪⎭⎪⎪⎬⎫==∑∑==iiiiiivNvuNu3131(1)式中：iN为节点坐标的线性函数。由于塑性极限分析不考虑本构关系，采用的是刚塑性模型，构造的是机动容许速度场，允许速度间断，因此，每个节点只属于某一个特定的单元，即不同单元的节点可具有相同的坐标，见图1。图1上限法三角形速度有限元模式Fig.1Triangleelementforvelocityfiniteelementmodelofupperlimitmethod3.3三角形单元的塑性流动约束条件对于理想刚塑性体，相关联流动法则可写成⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫∂∂=∂∂+∂∂=∂∂=∂∂=∂∂=∂∂=xyxyyyxxFxvyuFyvFxuτλγσλεσλε&&&&&&(2)式中：塑性应变率λ&≥0，同时假定拉应变为正。为了以后求解线性规划问题的方便，需将摩尔-库仑准则用若干平面线性化。对于平面应变问题，摩尔-库仑屈服准则可表示为−−+−=ϕτσσcos2[)2()(22cFxyyx0]sin)(2=+ϕσσyx(3)式中：正应力以拉为正；c，ϕ分别为材料的粘聚第23卷第11期王均星等.土坡稳定的有限元塑性极限分析上限法研究•1869·力和内摩擦角。令yxXσσ−=，xyYτ2=，cR2=·cosϕσσϕsin)(yx+−，在XY坐标系中，上述屈服条件为一个圆(见图2)，圆的半径与应力水平有关。用外接正多边形对屈服圆进行近似屈服[4]，屈服条件可写成)21(0pkFk，，，L==(4)式中：p为外接正多边形边数，即0cos2=−++=ϕτσσcCBAFxykykxkk(5)式中：⎪⎪⎭⎪⎪⎬⎫=+−=+=)/π2sin(2sin)/π2cos(sin)/π2cos(pkCpkBpkAkkkϕϕ(6)此时塑性应变率可写成⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫=∂∂=⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂=′=∂∂=∂∂==∂∂=∂∂=∑∑∑∑∑∑============pkkkkxykpkkkxykpkkkpkkyxkykpkkkpkkxxkxCFxvyuBFyvAFxu111111λτλγλσλελσλε&&&&&&&&(7)将式(1)，(5)，(6)代入式(7)，整理后可得三角形单元满足流动规则约束条件的矩阵表示形式：0212111=−xaxa(8)式中：⎥⎦⎤⎢⎣⎡=kkjjiikkjjiibcbcbccbcbcbA00000021e11a(9){}332211T1vuvuvu=x(10)⎪⎭⎪⎬⎫=+−=−=)(kjiixxcyybkjikji，，(11)式(9)中eA为三角形单元的面积，式(11)中)(iiyx，为节点坐标，而且有图2屈服条件线性化Fig.2Linearizationofyieldcondition⎥⎥⎦⎤⎢⎢⎣⎡=pkpkpkCCCCCBBBBBAAAAA..................32132132112a(12)}......{21T2pkλλλλ&&&&=x≥0(13)3.4速度不连续的约束条件假定速度不连续位于两个相邻单元的公共边上(见图3)，对应的结点编号为(1，2)和(3，4)，厚度为零。为了满足机动许可的条件，在该边上不连续的法向和切向速度间断值必须符合流动准则。摩尔-库仑屈服函数在速度不连续边所确定的局部坐标系n-s中可写为[5]0tan)(1=−+=cfnssnϕσττσ，(14)0tan)(2=−+−=cfnssnϕσττσ，(15)图3单元之间机动许可速度不连续面Fig.3Kinematicallyadmissiblevelocitydiscontinuitybetweenadjacentelements从而法向应变率与切向应变率可写为nnnffσλσλε∂∂+∂∂=2211&&&(16)sssffτλτλγ∂∂+∂∂=2211&&&(17)则法向速度和切向速度可写为=∂∂+∂∂==Δ∫∫→→tfftvnznzznztd)(limdlim22011000σλσλε&&&tzzd)(limtan0210∫+→λλϕ&&(18)=∂∂+∂∂==Δ∫∫→→tfftuszszzszsd)(limdlim22011000τλτλγ&&&tzzd)(lim0210∫−→λλ&&(19)令tkzzdlim0101∫→=λ&&，tkzzdlim0202∫→=λ&&，于是x=(σx－σy)R=2ccosϕ－(σx+σy)sinϕ摩尔-库仑屈服面x2+y2=R2线性化摩尔-库仑屈服面(p=6)y=2τxyθ21AtBs34zx，v(u1，v1)(u2，v2)(u3，v3)(u4，v4)y，u•1870•岩石力学与工程学报2004年ϕσσtan)(212211kkfkfkvnn&&&&+=∂∂+∂∂=Δ(20))(212211kkfkfkusss&&&&−=∂∂+∂∂=Δττ(21)如果速度不连续面与x轴的夹角为θ，利用坐标转换矩阵⎟⎠⎞⎜⎝⎛⎥⎦⎤⎢⎣⎡−=⎟⎠⎞⎜⎝⎛yxtsθθθθcossinsincos有⎟⎠⎞⎜⎝⎛⎥⎦⎤⎢⎣⎡−=⎟⎠⎞⎜⎝⎛vuvutsθθθθcossinsincos(22)因此，速度不连续约束可写成ϕθθtan)(cossin21kkvu&&+=Δ+Δ−(23))(sincos21kkvu&&−=Δ+Δθθ(24)式中：uΔ，vΔ为x，y方向的速度间断值。因此，速度不连续方程的矩阵表示为0323121=−xaxa(25)=21a⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡ABBAABBAθθθθθθθθsincos0000sincos00000000sincos0000sincos(26)式中：θsin−=A，θcos−=B⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−=1100tantan00001100tantan23ϕϕϕϕa(27){}{}⎪⎭⎪⎬⎫==0T)43(2)43(1)21(2)21(13T443322111kkkkkkkkvuvuvuvu，，，，&&&&xx(28)式中：k表示公共边的序号。3.5速度边界条件为了满足机动许可条件，计算的速度场必须满足已知的边界条件。设结点i位于与x轴夹角为θ的边界上，该边界上已知的切向速度和法向速度为，uv，此时结点i的速度分量必须满足以下方程：⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡−vuvuθθθθcossinsincos(29)写成矩阵形式为3131bxa=(30){}{}iivuvu，，==⎥⎦⎤⎢⎣⎡−=T1T331cossinsincosxbaθθθθ⎪⎪⎭⎪⎪⎬⎫(31)当边界上仅存在一种速度时，上述边界条件退化为一个方程。3.6内功功率与外功功率相等条件3.6.1内功功率内功功率是指应力分量在相应的应变增量上所做的功，又称内功耗散率。(1)连续体的内功功率对于每个三角形单元，内功功率可按下式计算：=++=∫∫AWxyxyyyxAxd)(eeτγσεσε&&&内=⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂+∂∂∫∫∑=AFFFApkxyxykyykxxkkde1ττσσσσλ&∫∫∑∑===e11ecos2d)cos2(ApkpkkkAcAcλϕϕλ&&(32)写成矩阵形式为2T2excW=内(33)⎪⎭⎪⎬⎫==0}{}11{cos221T2eT2pcAλλλϕ&L&&Lxc(34)(2)速度不连续面上的内功功率=+=∫∫→stslsnnzzldd)(lim000γτεσ&&内W+⎢⎣⎡⎟⎠⎞⎜⎝⎛∂∂+∂∂∫∫=lnnnzzff0221100limσλσλσ&&=+=⎥⎦⎤⎟⎠⎞⎜⎝⎛∂∂+∂∂∫∫→stcstffzlzsssdd)(limdd201002211λλτλτλτ&&&&)(d)(21201kkcllkkcl&&&&+=+∫(35)式中：l为速度不连续边的长度，c为速度不连续面的粘聚力。写成矩阵形式为3T3xcW=l内(36)式中：{}⎪⎭⎪⎬⎫==T)43(2)43(1)21((2)21(13T3}1111{21kkkkkkkkcl，，，，&&&&xc(37)3.6.2外功功率外功功率是指结点荷载(所有外荷载都通过等效移植而成，与常规有限元不同的是节点荷载不用累加)在结点速度(位移增量)上所做的功。写成矩阵≥≥第23卷第11期王均星等.土坡稳定的有限元塑性极限分析上限法研究•1871·形式为1T1xcW=外(38)式中：T22111][npnpvuvuvuLL=x为结点速度的行向量，xnpyxyxppppp{2211T1LL=c}ynpp为结点荷载的列向量；np为结点总数。3.6.3内功功率与外功功率相等的条件外内WW=(39)3.6.4屈服条件线性化对极限荷载的影响按照上限定理的推论，在其他条件相同的情况下，如果新的屈服面包含原屈服面，则新的极限荷载不小于原极限荷载。因此，屈服面用外接多边形线性化后，所求极限荷载不小于原屈服面的极限荷载。当外接多边形的边数无穷大时，新的屈服面趋向于同原屈服面重合。按照这一理论，屈服面按照本文方式线性化后所求得的极限荷载不会小于真实的极限荷载，符合上限定理的定义。4塑性极限分析线性规划模型在满足所列流动准则方程、速度边界条件和速度不连续约束条件的速度许可场中，通过虚功率原理，找出相应的破坏荷载，在众多破坏荷载中，最小的荷载即为所求极限荷载的上限值。上述所列方程或不等式都是单元节点位移分量和塑性乘子表示的线性关系式，可用线性规划方法求出极限状态的上限值。4.1求解超载系数的线性规划模型在边坡工程中，一般没有明确的外荷载，但间接等效荷载作用是存在的，比如开挖的影响，裂隙水压力的陡然增加等，可以等效当作面荷载处理。边坡的自重是边坡的主要荷载。在结构安全度的讨论中，曾定义超载安全系数和结构强度储备系数，求结构的极限荷载所对应的是结构的超载系数，而结构的强度储备系数可降低结构材料强度参数至破坏而得到。基于上述讨论，对于边坡工程的上限定理，通过虚功率原理可建立相应的求解超载系数和强度储备系数的目标函数。4.1.1求超载系数的目标函数将结点荷载分为超载部分和不超载部分，外功功率可改写为1T01T11}{}{xxkWpp+=外(40)式中：1}{p为超载的外荷载移置到结点上的等效结点荷载；0}{p为不超载的外荷载移置到结点上的等效结点荷载。由式(39)可得1T01T11}{}{xxkWpp+=内(41)由于我们求解的是破坏时结构的破坏形式，仅与}{1x的相对大小有关，令0.1}{1T1=xp(42)则1T01}{xWkp－内=(43)式中：T22111][npnpvuvuvuLL=x(44)对于边坡稳定问题，当求外荷载的最大值时，外荷载为可变荷载，容重为不变荷载；当求容重的最大值时，容重为可变荷载，外荷载为不变荷载。通过内功功率对单元和不连续面的累加及虚功率原理，可得到整体目标函数，见式(45)。4.1.2约束条件综合前面讨论的塑性流动条件、速度间断条件、速度边界条件，通过对式(8)，(25)，(29)进行集合，可得到总的约束方程，见式(46)。4.1.3线性规划1T13T32T21:minXCXCXCk−+=(45)subjectto：⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫==−=−000032313323121212111XXBXAXAXAXAXA(46)且附加约束条件0.1}{1T1=Xp(47)式中：T22111][npnpvuvuvuLL=X，=2XT)()(2)(1)1()2(2)1(1}{nepnenepλλλλλλ&L&&L&L&&，nlnlkkkkkk)21(2)21(1)43(2)43(1)21(2)21(13{1111，，，，，，&&L&&&&=XT)4,3(2)43(1}nlnlkk&&，，np为总结点数，ne为总单元数，nl为总速度不连续面数。文[1，2，6，7]总结了上述线性规划问题的特点，即它是一个约束条件个数远小于变量个数、异≥≥•1872•岩石力学与工程学报2004年常稀疏的超大规模线性规划问题。事实上，有限元塑性极限分析法的困难就在于大规模线性规划问题的求解。文[8]研究了改进单纯形法求解此类线性规划问题的时间复杂性和空间复杂性，认为对于列远远多于行的线性规划问题，改进单纯形法是比较经济的。另外，人们对改进的单纯形法积累了丰富的经验，如果在存储、误差减小方面采取有效的措施，可以取得较好的效果[9]。本文的算例采用改进单纯形法求解。根据以上思想，编制了上限程序UBK1.FOR，并调试成功。4.2求强度储备k2的最大值(安全系数)的线性规划模型对于边坡稳定分析问题，由于自重在具有下滑作用的同时也具有抗滑作用，所以用容重的超载系数定义安全系数不够准确。采用强度储备系数来定义安全系数，即002/tan/tancck==ϕϕ(48)但由流动约束条件和目标函数可知，这将导致一个非线性规划问题。文[10]中采用迭代方法求解，通过求解超载系数，使它充分接近1来求解强度储备系数，可以避免非线性规划问题。本文采用类似的迭代求解思想。设置2k的初值为0k，则以)/(tantan01kϕ−代替ϕ，以0/kc代替c，由式(45)，(46)，(47)求出超载系数1k，若0.11−k＜ε，则02kk=，结束运算；若1k＞ε+1，则令kkkΔ+=02；若1k＜ε+1，则令kkkΔ−=02，重复计算，当0.11→k时得到强度储备系数2k。上述ε为一任意小正数，且kΔ可以在计算过程中视1k与1.0的接近程度而调整。5算例图4所示均质土坡，坡比为1.75，土体的参数=γ18.0kN/m3，==ucc30kN/m2，0=ϕ，求解：图4边坡几何形状及计算简图Fig.4Geometryofslopeandsimplegraphforcalculation(1)当坡顶作用均布荷载时，均布荷载的最大值；(2)土坡在自然状态下的稳定安全系数。5.1求解极限荷载(相当于求超载系数)文[11]中给出了当0=γ时的解析解：上限解为)1(2α+=cP上，下限解为)sin1(2α+=cP下，其中α为斜坡面与y轴的夹角。此处=上P123.10kPa，=下P112.09kPa。由于结构完全左右对称，计算仅取一半，网格及边界条件如图5。表1中给出了屈服函数线性化为不同多边形时，边数与上限极限荷载的关系；表2给出了容重对极限荷载的影响。从表1可以看出，当屈服函数线性化边数L越图5边界条件及土坡计算网格Fig.5Boundaryconditionandfiniteelementmeshofsoilslope表1不同内接多边形的极限荷载Table1LimitloadatdifferentinscribedpolygonLP上/kPa误差/%4141.9615.326129.475.1712129.275.0124128.294.2148128.194.1360128.174.11表2不同边坡容重的极限荷载Table2Limitloadatdifferentweightofslopeγ/kN·m－3P上/kPa0.01129.261128.364125.096122.7210117.2112114.1916107.342097.672481.583049.2032m4m8m1∶1.751∶1.75=u0u=0,0=vu=0u=0，v=0第23卷第11期王均星等.土坡稳定的有限元塑性极限分析上限法研究•1873·大时，解越接近精确解，但L较大时，变化不太敏感。对于一般问题，12=L时可满足精度 要求 对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗 。理论分析表明，当ϕ比较大时，L应取大些。从表2可以看出，随容重增大，极限承载能力开始对容重的变化并不敏感，但当γ比较大时极限承载力随γ的增大而显著降低。图5中对应的速度场为破坏时的破坏机构场，大体上反映了滑动的趋势。5.2求安全系数根据式(42)的算法，求解上述土坡得到的安全系数11.2=k，而采用传统的条分法得到的安全系数为=k2.03。由于本文对安全系数的定义与传统的条分法是一致的，所计算的结果相近，证明了该方法的正确性。事实上，从塑性力学的观点来看，一方面，条分法属于刚体极限平衡法，它仅考虑了力的平衡条件，而不能考虑结构的屈服，其解答可归属下限解的范畴，但它只是无穷多个可能下限解中的一个，不一定最大；另一方面，由于滑动面是人为假定的，这意味着已经假定结构在极限状态下满足破坏机构条件，其解答又可归属上限解的范畴。总之，刚体极限平衡法的物理意义是比较模糊的，而上限法从上限逼近问题的真实解，物理意义明确。6结论边坡稳定的有限元塑性极限分析法有严格的理论基础和明确的物理意义，算例表明该方法具有较高的计算精度。但是，这种方法特别依赖于好的线性规划求解器，实际问题的规模往往是超巨型的，值得进一步探讨和研究。另外，本文所举算例虽然网格比较稀疏，但仍然得到了较为精确的结果。研究[1，2，5，6]表明，上限法对网格的疏密要求不高。当然，网格越细，精度越高，所花费的机时也越多。参考文献1SloanSW.Upperboundlimitanalysisusingfiniteelementandlinearprogramming[J].Int.J.AnalyticalMethodsinGeomechanics，1989，13(3)：263～282.2SloanSW.Upperboundlimitanalysisusingdiscontinuousvelocityfields[J].J.ComputerMethodsinAppliedMechanicsandEngineering，1995，127(4)：293～3143王开治，王均星.重力坝塑性极限分析[J].大坝与安全，2000，(4)：12～184龚晓南.土工计算机分析[M].北京：中国建筑工业出版社，19995ZhengX.Limitanalysisofthebearingoffissuredmaterials[J].Int.J.SolidsandStructures，2000，37(8)：1211～12436KimJ.Limitanalysisofsoilslopessubjectedtopore-waterpressures[J].J.GeotechnicalandGeoenvironmentalEngineering，1999，125(1)：49～587SloanSW.Asteepestedgeactivesetalgorithmforsolvingsparselinearprogrammingproblems[J].Int.J.forNumericalMethodsinEngineering，1988，26(12)：2671～26858卢开澄.单目标、多目标与整数规划[M].北京：清华大学出版社，19999杨冰.实用最优化方法及计算机程序[M].哈尔滨：哈尔滨船舶工程学院出版社，199410李国英，沈珠江.下限原理有限单元法及其在土工问题中的运用[J].岩土工程学报，1997，19(5)：84～8911张学言.岩土塑性力学[M].北京：人民交通出版社，1993