随着计算机水平的进步，有限元极限分析法［2］发展迅速，通过把有限元思想、数学线性规划技术与极限分析理论相结合，可以直接计算出二维复杂模型下的下限或上限解，这在岩土力学领域具有重要意义。此后，大量学者开展了相关研究。Sloan［3-4］在最速下降活动集法的基础上，提出并建立了求解效率更高的线性规划理论，而且可以在台式计算机上算出中等规模的二维极限分析问题；王均星等［5］成功在有限元下限分析中引入孔隙水压力；李春光等［6］通过建立四边形网格积分意义上平衡方程的弱形式，建立了基于四边形网格的下限有限元法计算模型；赵明华等［7］运用有限元极限分析法探讨了双孔隧道在支护反力及地面超载共同作用下的相互规律；孙锐等［8］基于复合锥优化技术，提出一种基于广义Hoek-Brown屈服准则的极限分析自适应下限有限元计算方法；许晓亮等［9］结合随机场理论，提出了网格自适应的边坡可靠性随机有限元极限分析方法。
在线性规划理论的基础上，有限元极限分析理论在中小规模的二维问题中能够取得良好结果。但是，同样的方法并不能很好地推广运用到三维问题中，因为三维问题的计算规模过于庞大，而且三维问题中6个线性莫尔库仑平面相交会导致屈服面梯度不连续，这给应用常规的线性规划方法求解带来了较大的困难［10］。
徐干成等［11］按照π平面上与摩尔库伦准则面积相等原理提出等面积圆屈服准则［12］来逼近，计算结果与真实的摩尔库伦理论相当，取得了良好的逼近效果，但没有学者从三维下限分析进行研究。本文采用等面积德鲁克-普拉格（Drucker-Prager，DP）屈服准则来代替摩尔库伦屈服准则，并采用非线性规划的二阶锥规划成功构建了三维问题的极限分析下限法计算模型，可以求解三维边坡的应力场和安全系数，2个算例证明了该方法的可行性。
1 DP屈服准则
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]的表达式为：
[α=2sinφ33-sinφ,k=3cαcotφ] （4）
若与摩尔库伦屈服面内角点外接，则[α]和[k]的表达式为：
[α=2sinφ33+sinφ, k=3cαcotφ] （5）
若与摩尔库伦屈服面内切，则[α]和[k]的表达式为：
[α=sinφ9+3sin2φ, k=3cαcotφ] （6）
若与摩尔库伦屈服面面积相等，则[α]和[k]的表达式为：
[α=23sinφ2π39-sin2φ,k=3cαcotφ] （7）
本文采用等面积DP屈服准则来开展研究。
<G:\武汉工程大学\2023\第3期\孙 聪-1.tif>
图1 三维应力空间中的莫尔-库仑屈服面和等效面积
D-P屈服面
Fig. 1 View of Mohr-Coulomb yield surface and equivalent area D-P yield surface in three-dimensional space of
principal stresses
2 DP准则的二阶锥形式表达
令：
[σ=σxσyσzτyzτzxτxyT] （8）
则第一应力不变量可写成：
[I1=eσ] （9）
其中：
[e=111000] （10）
第二应力不变量可表示为：
[J2=σTMσ] （11）
[M=163I-1I] （12）
其中，I是3×3的单位矩阵，1是3×3的元素全为1的矩阵。
式（11）同样可以表达成：
[J2=σTATAσ] （13）
[A=263I-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,iσy,iσz,iσyz,iσzx,iσxy,iT] （17）
那么单元中任意一点p（x，y，z）处的应力表示为：
[σ=i=14Niσi] （18）
式（18）中[Ni]为四面体单元的形函数，是线性的，表达式为：
[Ni=ViV] （19）
[V=161111x1x2x3x4y1y2y3y4z1z2z3z4] （20）
V是四面体单元的体积，[Vi]是单元内一点p和除节点i外的其他节点组成的四面体的体积，[Vi]和V满足：
[i=14Vi=V] （21）
<G:\武汉工程大学\2023\第3期\孙 聪-2.tif>
图2 下限分析的四节点线性应力四面体
Fig. 2 Four-node linear stress tetrahedron for
lower bound limit analysis
3.2 应力平衡
四面体单元内的应力平衡方程可表达为：
[?σx?x+?τxy?y+?τxz?z+bx=0?τxy?x+?σy?y+?τxz?z+by=0?τxz?x+?τxz?y+?σz?z+bz=0] （22）
其中，b=［bx，by，bz］T是单元内的体积力。将式（18）代入式（22）可得：
[Bσe=b] （23）
[σe=σ1σ2σ3σ4T] （24）
[B=B1B2B3B4] （25）
[Bi=?Ni?x?Ni?z?Ni?y?Ni?y?Ni?z?Ni?x?Ni?z?Ni?y?Ni?x] （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）
<G:\武汉工程大学\2023\第3期\孙 聪-3.tif>[3][6][2][5][1][4][a][b]
图3 相邻四面体公共面上的应力连续条件
Fig. 3 Stress continuous condition between adjacent
tetrahedrons
3.4 应力边界条件
由于四面体单元是线性单元，可以将每个节点i的应力边界条件转换为：
[σT, i=ti] （29）
其中表面力[ti=ti,1ti,2ti,3T]。
3.5 三维下限法的二阶锥规划模型
通过引入平衡方程、应力边界条件和应力连续条件，三维问题的下限极限分析可表达为：
[Maximize λ] （30）
[Subject to Bσ+λb=0]（平衡方程） （31）
[σT=λt]（边界条件） （32）
[σT, a=σT, b]（应力连续条件） （33）
[Aσ]+[αbσ-k≤0 ]（屈服条件） （34）
其中[λ]为超载系数。
可以看到，这是一个典型的二阶锥规划，可以通过现成的凸规划软件（比如斯坦福大学Stephen P. 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为不同数量网格离散后，采用本文二阶锥优化方法计算得到的对应安全系数，可以清楚地看到，有限元下限法得到的结果和单元剖分数量密切相关。
表2 算例1不同单元尺寸时的安全系数
Tab. 2 Effects of mesh refinement on safety factor for slope 1
[编号 单元数量 安全系数 1 42 1.362 6 2 3 879 1.714 6 3 9 326 1.730 2 4 16 262 1.731 0 5 20 380 1.734 0 6 27 559 1.734 6 ]
从表2可以看出，随着单元数量的增加，安全系数从低到高逐渐收敛到1.734 6（从下方收敛），并且与使用莫尔库仑屈服准则（安全系数<1.734 4）的弹塑性有限元法获得的解［14］非常一致。
4.2 带有拐角的三维边坡
算例2为一转角为90°的代表性三维边坡，如图6所示，材料参数见表3。Nian等［15］使用满足摩尔库伦屈服准则的弹塑性有限元法（非关联流动法则）对本例进行了研究，该边坡的边界条件为：底面和图6（a）中阴影部分全固定，图6（b）中阴影部分仅固定垂直于平面方向。
<G:\武汉工程大学\2023\第3期\孙 聪-6-1.tif><G:\武汉工程大学\2023\第3期\孙 聪-6-2.tif>[ a ][ b ][10][10][60][60][40][40][20][20]
图6 算例2边坡的尺寸（单位：m）
Fig. 6 Geometry of slope 2（Unit：m）
表3 算例2边坡的材料参数
Tab. 3 Material property of slope 2
[黏聚力c / kPa 摩擦角φ / （°） 重度γ / （kN/m3） 40 20 20 ]
为了研究对安全系数的影响，同样采用6种不同尺寸的网格进行划分，图7为代表性的3种网格数量下的剖分示意图，表4为不同网格数量对应的安全系数。可以清楚地看出，当单元数足够多时，安全系数将接近稳定值。
表4 算例2不同单元尺寸时的安全系数
Tab. 4 Effects of mesh refinement on safety factor for slope 2
[编号 单元数量 安全系数 1 1 309 2.570 0 2 2 911 2.712 9 3 8 450 2.808 4 4 18 162 2.845 8 5 27 688 2.859 2 6 32 102 2.861 5 ]
从表4可以看出，根据本文下限分析方法计算得到的安全系数为2.861 5，与Nian等［15］使用ABAQUS莫尔-库仑屈服准则有限元法得出的结果（2.860）几乎相同。
5 结 论
（1）对于三维下限分析，DP准则可以写成二阶锥形式，这样就可以将三维下限问题表达成一个二阶锥规划模型，从而可以采用开源的免费代码CVX来求解。
（2）在求解三维边坡安全系数时，采用等面积圆DP准则极限分析下限法得到的解与采用摩尔库伦弹塑性有限元得到的解非常接近。