基于二阶锥规划的三维边坡极限分析下限法研究.pdf

武汉工程大学学报第45卷边坡的稳定性分析一直是岩土工程领域的经典问题，其中，基于莫尔库仑屈服准则的极限平衡法应用最为广泛，它通常被用来计算边坡的安全系数，但由于其理论基础不够严谨，引入了大量条间力假设，有时可能会导致较大的计算误差1。随着计算机水平的进步，有限元极限分析法2发展迅速，通过把有限元思想、数学线性规划技术与极限分析理论相结合，可以直接计算出二维复杂模型下的下限或上限解，这在岩土力学领域具有重要意义。此后，大量学者开展了相关研收稿日期：2022-10-31基金项目：中国科学院国防科技创新基金（CXJJ-16M272）作者简介：孙聪，博士，高级工程师。E-mail：引文格式：孙聪，熊威，李春光.基于二阶锥规划的三维边坡极限分析下限法研究 J.武汉工程大学学报，2023，45（4）：468-472.基于二阶锥规划的三维边坡极限分析下限法研究孙聪1，熊威2，李春光31.武汉生态环境设计研究院有限公司，湖北 武汉 430050；2.长江勘测规划设计研究有限责任公司，湖北 武汉 430010；3.中国科学院武汉岩土力学研究所，湖北 武汉 430071摘要：由于屈服面梯度不连续，导致基于摩尔库伦屈服准则的三维边坡极限分析问题无法直接使用线性规划求解，从而计算效率低。采用等面积德鲁克-普拉格屈服准则来代替摩尔库伦屈服准则，并在考虑应力平衡、应力连续性条件和应力边界条件的基础上，用非线性规划的二阶锥规划成功构建了三维问题的极限分析下限法计算模型，可以方便地通过凸规划软件直接求解出三维边坡的应力场和安全系数，从真实解的下方收敛，两个算例证明了其可行性。关键词：下限分析；德鲁克-普拉格准则；二阶锥规划；凸规划；边坡稳定性中图分类号：TD712文献标识码：ADOI：10.19843/ki.CN42-1779/TQ.202210034Lower Bound Limit Analysis of Three-Dimensional Slopes UsingSecond-Order Cone ProgrammingSUN Cong1，XIONG Wei2，LI Chunguang31.Wuhan Ecological Environment Design and Research Institute Co.，Ltd，Wuhan 430050，China；2.Changjiang Survey，Planning and Design Research Co.，Ltd，Wuhan 430010，China；3.Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，Wuhan 430071，ChinaAbstract：Due to the discontinuous gradient of the yield surface，the limit analysis of three-dimensional slopebased on Mohr-Coulomb yield criterion can not be solved directly by linear programming，which causes aninefficient calculation.In this paper，the equal-area Drucker-Pluger yield criterion was used to replace theMohr-Coulomb yield criterion.Considering the stress balance，stress continuity conditions and stress boundaryconditions，we successfully constructed the limit analysis lower bound method calculation model of three-dimensional problems by using the second-order cone programming of nonlinear programming.The stressfield and safety factor of three-dimensional slope can be directly solved by convex programming software，from the lower convergence of the real solution.Two examples are given to prove its feasibility.Keywords：lowerlimitanalysis；Drucker-Pragercriterion；second-orderconeprogramming；convexprogramming；slope stability第45卷第4期2023年8月文章编号：1674-2869（2023）04-0468-05武汉工程大学学报Journal of Wuhan Institute of TechnologyVol.45 No.4Aug.2023第4期究。Sloan3-4在最速下降活动集法的基础上，提出并建立了求解效率更高的线性规划理论，而且可以在台式计算机上算出中等规模的二维极限分析问题；王均星等5成功在有限元下限分析中引入孔隙水压力；李春光等6通过建立四边形网格积分意义上平衡方程的弱形式，建立了基于四边形网格的下限有限元法计算模型；赵明华等7运用有限元极限分析法探讨了双孔隧道在支护反力及地面超载共同作用下的相互规律；孙锐等8基于复合锥优化技术，提出一种基于广义 Hoek-Brown屈服准则的极限分析自适应下限有限元计算方法；许晓亮等9结合随机场理论，提出了网格自适应的边坡可靠性随机有限元极限分析方法。在线性规划理论的基础上，有限元极限分析理论在中小规模的二维问题中能够取得良好结果。但是，同样的方法并不能很好地推广运用到三维问题中，因为三维问题的计算规模过于庞大，而且三维问题中 6个线性莫尔库仑平面相交会导致屈服面梯度不连续，这给应用常规的线性规划方法求解带来了较大的困难10。徐干成等11按照 平面上与摩尔库伦准则面积相等原理提出等面积圆屈服准则12来逼近，计算结果与真实的摩尔库伦理论相当，取得了良好的逼近效果，但没有学者从三维下限分析进行研究。本文采用等面积德鲁克-普拉格（Drucker-Prager，DP）屈服准则来代替摩尔库伦屈服准则，并采用非线性规划的二阶锥规划成功构建了三维问题的极限分析下限法计算模型，可以求解三维边坡的应力场和安全系数，2个算例证明了该方法的可行性。1DP屈服准则DP 屈服准则常被用来处理岩土工程中的塑性变形问题，其公式如下：J2+I1-k=0（1）其中，I1=x+y+z（2）J2=16()x-y2+()y-z2+()z-x2+2xy+2yz+2zx（3）I1和 J2分别是应力偏张量的第一不变量和第二不变量，和k是与黏聚力 c 和摩擦角 有关的常数，为正应力，为剪应力。如图1所示，摩尔库伦屈服面在应力空间中是一个不规则的六边形截面角锥体，而DP屈服面在应力空间中是一个圆锥面。假设在 平面上，DP屈服面与摩尔库伦屈服面外角点外接13，则和k的表达式为：=2sin3()3-sin,k=3ccot（4）若与摩尔库伦屈服面内角点外接，则和k的表达式为：=2sin3()3+sin,k=3ccot（5）若与摩尔库伦屈服面内切，则和k的表达式为：=sin9+3sin2,k=3ccot（6）若与摩尔库伦屈服面面积相等，则和k的表达式为：=23 sin23()9-sin2,k=3ccot（7）本文采用等面积DP屈服准则来开展研究。-1-2-3图1三维应力空间中的莫尔-库仑屈服面和等效面积D-P屈服面Fig.1View of Mohr-Coulomb yield surface and equivalentarea D-P yield surface in three-dimensional space ofprincipal stresses2DP准则的二阶锥形式表达令：=xyzyzzxxyT（8）则第一应力不变量可写成：I1=e（9）其中：e=111000（10）第二应力不变量可表示为：J2=TM（11）M=16()3I-1I（12）其中，I是 33的单位矩阵，1是 33的元素全孙聪，等：基于二阶锥规划的三维边坡极限分析下限法研究469武汉工程大学学报第45卷为1的矩阵。式（11）同样可以表达成：J2=TATA（13）A=26()3I-1I（14）J2可表达为如下结构形式：J2=A（15）表示矢量的大小。因此，D-P屈服准则可以被重新表述为以下二阶锥形式：A+b-k=0（16）3极限分析下限法的二阶锥规划模型3.1四节点四面体单元离散对于三维问题，可以用四节点四面体来构建有限单元，如图 2 所示，单元节点 i（i=1，2，3，4）处的应力为：i=x,iy,iz,iyz,izx,ixy,iT（17）那么单元中任意一点 p（x，y，z）处的应力表示为：=i=14Nii（18）式（18）中Ni为四面体单元的形函数，是线性的，表达式为：Ni=ViV（19）V=16|1111x1x2x3x4y1y2y3y4z1z2z3z4（20）V是四面体单元的体积，Vi是单元内一点 p和除节点 i外的其他节点组成的四面体的体积，Vi和V满足：i=14Vi=V（21）1423图2下限分析的四节点线性应力四面体Fig.2Four-node linear stress tetrahedron forlower bound limit analysis3.2应力平衡四面体单元内的应力平衡方程可表达为：xx+xyy+xzz+bx=0 xyx+yy+xzz+by=0 xzx+xzy+zz+bz=0（22）其中，b=bx，by，bzT是单元内的体积力。将式（18）代入式（22）可得：Be=b（23）e=1234T（24）B=B1B2B3B4（25）Bi=NixNizNiyNiyNizNixNizNiyNix（26）因为Ni是线性的，所以B是一个常数矩阵。3.3间断面上应力连续性条件对于三维问题下限法，必须保证作用在相邻四面体单元公共面上的正应力和剪应力分量相等。如图 3所示，对于四面体 a和 b，公共面 T上的任一点的应力须满足：T,a=T,b（27）由于公共面 T上的应力线性变化，这相当于对 3个公共节点施加约束：T,1=T,4；T,2=T,5；T,3=T,6（28）362514ab图3相邻四面体公共面上的应力连续条件Fig.3Stress continuous condition between adjacenttetrahedrons3.4应力边界条件由于四面体单元是线性单元，可以将每个节点i的应力边界条件转换为：T,i=ti（29）其中表面力ti=ti,1ti,2ti,3T。3.5三维下限法的二阶锥规划模型通过引入平衡方程、应力边界条件和应力连470第4期续条件，三维问题的下限极限分析可表达为：Maximize（30）Subject to B+b=0（平衡方程）（31）T=t（边界条件）（32）T,a=T,b（应力连续条件）（33）A+b-k 0（屈服条件）（34）其中为超载系数。可以看到，这是一个典型的二阶锥规划，可以通过现成的凸规划软件（比如斯坦福大学 StephenP.Boyd教授基于Matlab开发的CVX工具箱）求解。4算例4.1常规三维边坡如图 4所示，该均质边坡坡比为 2 1，Griffiths等14研究了不同的长度（L）与高度（H）比（L/H）的影响，在本文中，仅采用L/H=2的情况来研究本文方法的可行性。边坡的材料参数为：=2 000 kg/m3，c=40 kPa，=0，尺寸见表1。用 6种不同尺寸的四面体网格单元来离散该边坡，其中有代表性的 3种如图 5所示。表 2为不同数量网格离散后，采用本文二阶锥优化方法计算得到的对应安全系数，可以清楚地看到，有限元下限法得到的结果和单元剖分数量密切相关。图5算例1三种网格密度剖分Fig.5Three grid divisions of slope 1表2算例1不同单元尺寸时的安全系数Tab.2Effects of mesh refinement on safety factor for slope 1编号123456单元数量423 8799 32616 26220 38027 559安全系数1.362 61.714 61.730 21.731 01.734 01.734 6从表 2可以看出，随着单元数量的增加，安全系数从低到高逐渐收敛到 1.734 6（从下方收敛），并且与使用莫尔库仑屈服准则（安全系数1.734 4）的弹塑性有限元法获得的解14非常一致。4.2带有拐角的三维边坡算例 2为一转角为 90的代表性三维边坡，如图6所示，材料参数见表3。Nian等15使用满足摩尔库伦屈服准则的弹塑性有限元法（非关联流动法则）对本例进行了研究，该边坡的边界条件为：底面和图 6（a）中阴影部分全固定，图 6（b）中阴影部分仅固定垂直于平面方向。ab1010606040402020图6算例2边坡的尺寸（单位：m）Fig.6Geometry of slope 2（Unit：m）表3算例2边坡的材料参数Tab.3Material property of slope 2黏聚力c/kPa40摩擦角/（）20重度/（kN/m3）20网格数量3 879网格数量16 262网格数量27 559W1W2SxzHDyL/2图4算例1边坡示意图Fig.4Geometry of slope 1表1算例1边坡尺寸（单位：m）Tab.1Dimensions of slope 1（unit：m）D=1.5H15S=2H20L/2=H10W1=H10W2=H10孙聪，等：基于二阶锥规划的三维边坡极限分析下限法研究471武汉工程大学学报第45卷为了研究对安全系数的影响，同样采用6种不同尺寸的网格进行划分，图7为代表性的3种网格数量下的剖分示意图，表4为不同网格数量对应的安全系数。可以清楚地看出，当单元数足够多时，安全系数将接近稳定值。表4算例2不同单元尺寸时的安全系数Tab.4Effects of mesh refinement on safety factor for slope 2编号123456单元数量1 3092 9118 45018 16227 68832 102安全系数2.570 02.712 92.808 42.845 82.859 22.861 5从表4可以看出，根据本文下限分析方法计算得 到 的 安 全 系 数 为 2.861 5，与 Nian 等15使 用ABAQUS 莫尔-库仑屈服准则有限元法得出的结果（2.860）几乎相同。5结论（1）对于三维下限分析，DP 准则可以写成二阶锥形式，这样就可以将三维下限问题表达成一个二阶锥规划模型，从而可以采用开源的免费代码CVX来求解。（2）在求解三维边坡安全系数时，采用等面积圆 DP 准则极限分析下限法得到的解与采用摩尔库伦弹塑性有限元得到的解非常接近。参考文献1 TSCHUCHNIGG F，SCHWEIGER H F，SLOAN S W，et al.Comparison of finite-element limit analysis andstrength reduction techniquesJ.Gotechnique，2015，65（4）：249-257.2LYSMER J.Limit analysis of plane problems in soilmechanics J.Journal of the Soil Mechanics&Foundations Division，1970，96（4）：1311-1334.3SLOAN S W.Lower bound limit analysis using finiteelements and linear programming J.InternationalJournal for Numerical and Analytical Methods inGeomechanics，1988，12（1）：61-77.4SLOAN S W.Upper bound limit analysis using finiteelements and linear programming J.InternationalJournal for Numerical and Analytical Methods inGeomechanics，1989，13（3）：263-282.5 王均星，李泽，陈炜.考虑孔隙水压力的土坡稳定性的有限元下限分析 J.岩土力学，2005，26（8）：1258-1262，1268.6 李春光，朱宇飞，刘丰，等.基于四边形网格的下限原理有限元法 J.岩石力学与工程学报，2012，31（3）：461-468.7 赵明华，侯继超，彭文哲，等.基于有限元极限分析的双孔隧道稳定性分析J.地下空间与工程学报，2020，16（1）：234-242.8 孙锐，阳军生，李雨哲，等.基于广义 Hoek-Brown 屈服准则的极限分析下限有限元法J.岩土力学，2021，42（6）：1733-1742.9 许晓亮，张家富，曾林风，等.考虑网格自适应的边坡可靠性随机有限元极限分析研究 J.三峡大学学报（自然科学版），2022，44（6）：48-57.10 KRABBENHFT K，LYAMIN A V，SLOAN S W.Formulation and solution of some plasticity problemsas conic programsJ.International Journal of Solidsand Structures，2007，44（5）：1533-1549.11 徐干成，郑颖人.岩石工程中屈服准则应用的研究J.岩土工程学报，1990，12（2）：93-99.12 梁俊伟，郭昆，韩丽娟，等.基于等面积圆D-P屈服准则的波纹板挡土墙应用研究 J.青海交通科技，2019（2）：99-102.13 时卫民，郑颖人.摩尔-库仑屈服准则的等效变换及其在边坡分析中的应用 J.岩土工程技术，2003（3）：155-159.14 GRIFFITHSDV，MARQUEZRM.Three-dimensional slope stability analysis by elasto-plasticfinite elementsJ.Gotechnique，2007，57（6）：537-546.15NIAN T K，HUANGRQ，WAN SS，etal.Three-dimensional strength-reduction finite elementanalysisofslopes J.CanadianGeotechnicalJournal，2012，49（5）：574-588.本文编辑：苗变图7算例2三种网格密度剖分Fig.7Three grid divisions of slope 2网格数量1 658网格数量11 438网格数量24 787472