图1为挡土墙主动土压力的计算模型。设挡土墙高度为H，地面均布超载为q，墙背水平倾角为α，地面水平倾角为β(负值表示坡面向下倾斜)，填土重度为γ，黏聚力为c，内摩擦角为φ，填土与墙背之间的黏结力为c_w，外摩擦角为δ。基本假定如下：

图  1  本文计算模型

Figure 1.  Calculation model used in this paper

1)按平面应变问题考虑。

2)墙后填土为服从Mohr-Coulomb屈服准则的理想弹塑性材料，而且是各向同性的、均匀的以及不可压缩(膨胀)的理想连续介质，由3个常参数γ、c、φ表征。

3)墙体发生充分的变位使墙后土楔ABC处处处于塑性平衡状态。

4)墙土接触面AB满足库仑摩擦定律，由2个常参数c_w、δ表征。

5)弹塑性区界面AC为平面且平行于地面DE。AC面上作用均布压力：

其中：

式中：h₀为塑性临界深度；z₀为弹性覆盖层的垂直深度，即通常所谓的地面张拉裂缝深度。

图1中塑性区内任意一点P的摩尔应力圆如图2所示，通过P点有两个极限滑裂面即α滑移线和β滑移线，它们与大主应力方向的夹角均为μ＝π/4-φ/2。设大主应力方向与y轴的夹角为θ，规定θ以y轴正向逆时针旋转到指定方向为正，则两条滑移线与y轴的夹角分别为θ-μ和θ＋μ，故滑移线方程为：

图  2  P点摩尔应力圆

Figure 2.  Mohr′s stress circle at point P

根据摩尔应力圆，可得正应力和剪应力

式中：p＝(σ_x＋σ_y)/2为平均应力；R＝p sin φ＋c cos φ为应力圆半径。

将σ_x、σ_y及τ_xy代入以下静力平衡方程：

式中：X、Y分别为水平和垂直方向的单位体积力。

经整理化简后可得到沿α滑移线和β滑移线的极限平衡方程：

式中：

由于沿α滑移线有dS_α＝dy/ cos (θ-μ)，沿β滑移线有dS_β＝dy/ cos (θ＋μ)，故方程式(7)可进一步改写为：

沿α滑移线和β滑移线的土反力分布可按下式计算：

土反力分布r_α和r_β的作用线与＋y方向的夹角分别为θ＋μ和θ-μ。

滑移线曲率半径可按下式计算：

R_α以α滑移线曲率中心处于＋y方向为正，R_β以β滑移线曲率中心处于＋x方向为正。

如图3所示，M为墙土接触面AB上的任意一点，设该点主动土压力强度为p_a，则有：

图  3  M点摩尔应力圆

Figure 3.  Mohr′s stress circle at point M

式中：σ_m、τ_m分别为边界AB上的正应力和剪应力。又由摩尔应力圆可得：

联立式(13)和式(14)并消去p_a，可得边界AB的应力边界条件：

式中：R_M＝p_M sin φ＋c cos φ，R_w＝p_M sin δ＋c_w cos δ。

当c_wtanφ＝ctan δ时，该边界条件与平均应力p_M无关，可作进一步简化。

又联立式(13)和式(14)并消去θ_M，可得到关于p_a的一元二次方程，当R_M≥R_w时有两个实根，在主动极限平衡状态下，应取较小根，即：

可见，一旦求解了边界AB上的平均应力p_M，则主动土压力可由上式计算。特别地，当c_w＝c＝0时，式(16)可简化为：

式中： 。

应当指出，Sokolovskii(1965)给出了一个类似的公式，但其右端项的分母中少了一个cos δ。

如图4所示，N为弹塑性区交界面AC上的任一点，需求解应力边界值p_N和θ_N。根据摩尔应力圆方程：

图  4  N点摩尔应力圆

Figure 4.  Mohr′s stress circle at point N

式中：R_N＝p_N sin φ＋c cos φ，当 cos ²β(tan|β|- tan φ)≤c时，以上关于p_N的一元二次方程有两个实根，在主动极限平衡状态下，应取较小根，即：

式中： ，又由应力圆可得τ_n＝R_N sin 2(β-θ_N)。

现分别考虑以下两种情况：

1)当 cos ²β(tan|β|-sin φ)≤c cos φ时：

2)当 cos ²β(tan|β|- sin φ)>c cos φ时：

特别地，当c_w＝c＝0时，式(20)或式(21)可进一步简化为：

当图1中的α滑移线MN无限逼近A点时，M点和N点的应力状态一般不相同，即平均应力p_M→A和p_N→A，方向角θ_M→A和θ_N→A一般不相等，称A点为应力奇点。现令dy＝0，并对式(9)进行积分，可得到应力奇点的应力边界条件：

式中：C_α为待定积分常数，θ_A在θ_M→A和θ_N→A之间连续变化。

由此可见，应力奇点实际上就是α滑移线退化为点的必然产物，是一种特殊的应力边界，具有应力连续、应力梯度无穷大以及α滑移线曲率无穷大的特性。

在计算主动土压力之前，务必确定塑性临界深度h₀。首先求解塑性临界压力 ，计算步骤如下：

1)假设图1中存在弹性覆盖层，即h₀>0，然后令式(16)的p_a等于0，可求得当M点无限逼近A点时的平均应力：

2)将p_M→A代入式(15)求得相应的θ_M→A。

3)将p_M→A和θ_M→A代入式(23)得到积分常数C_α。

4)联立式(23)、式(19)以及式(20)或(21)，采用迭代计算求解塑性临界压力 。

在某些特别情况下无需迭代计算，例如当β＝0时， 可直接按下式计算：

当满足以下条件时：

可按下式计算：

式中：应力圆半径R_M→A＝p_M→A sin φ＋c cos φ；计算结果与文献[31]的解析解完全相同。

一般而言，塑性临界压力 是关于α、β、φ、c及c_w的函数，而与H、q、γ、δ无关。可以证明，当c_w＝c＝0时，恒有 ；当c>0时，无量纲数 是关于α、β、φ及 的函数。

作用于AC面上的均布压力 应取q和 二者中的较大值，即 。如地面超载q< ，则 ，故欲使A点进入主动极限平衡状态尚需由土重补充不足部分的压力，即存在弹性覆盖层，其垂直深度z₀可按式(1)计算，再将z₀代入式(2)可求得塑性临界深度h₀，即

如 ，则 ，由于z₀和h₀等于0，故不存在弹性覆盖层。

当图1中的α滑移线MN无限逼近A点时，如M点和N点的应力状态相同，则A点为非应力奇点，而将p_M→A＝p_N→A和θ_M→A＝θ_N→A称为非奇异条件，它与参数α、β、φ、δ、c、c_w及 有关。应分别考虑以下两种情况：

当 cos ²β(tan|β|- sin φ)≤c cos φ时，非奇异条件可改写为：

当β cos ²β(tan|β|- sin φ)>c cos φ时，

式中：ρ₁、ρ₂按下式计算：

若 且满足式(26)，例如α＝π/2，c_w＝0或α＝π/4＋φ/2，β＝0，c_w＝c，则恒有Δ＝π，非奇异条件满足；若 且α＝π/2，c_w＝0，β＝δ或α＝π/4＋φ/2，β＝0，c_w＝c，δ＝φ，则恒有Δ＝π，非奇异条件也满足。当c_w＝c＝0时，式(29)或(30)可简化为

上式就是经典库仑主动土压力的非奇异条件。应当指出，经典库仑土压力不一定满足Δ＝π，而朗肯土压力则总是满足Δ＝π。

由以上滑移线方程、极限平衡方程以及应力边界条件构成了一个静定可解的边值问题，一般采用有限差分法求解。求解时需事先构建有限差分网格，它由两族互为相交的、假想的滑移线组成，如图5所示。一般地，设α族滑移线的总条数为n，第一条α滑移线A′A″的结点数为m(m≥2)，则β族滑移线的总条数为2(n-2)＋m，网格结点总数为n(n＋m-1)，而结点按图中顺序依次进行编号。本文计算机程序暂定n不大于99，结点总数不大于9 999个。

图  5  差分网格示意图

Figure 5.  Finite difference mesh

现研究当第一条α滑移线A′A″无限逼近应力奇点A时的情况。根据式(19)和式(20)或(21)可计算A″点的p_A″→A和θ_A″→A，将它们代入式(23)求得积分常数C_α，于是应力奇点A的应力边界条件就确定了。又联立式(15)、式(23)，采用迭代算法可求得A′点的p_A′→A和θ_A′→A。一般情况下θ_A′→A和θ_A″→A不相等，若将方向角差θ_A″→A-θ_A′→A平分为m-1等份，并假定等分点和A′A″上的结点相对应，则第一条α滑移线的m个结点可由下式确定：

设N₁为边界AC的结点，N₂为紧邻N₁的一个结点，且N₁N₂为α滑移线。已知结点N₂的x₂、y₂、p₂和θ₂，可采用迭代算法求解结点N₁的x₁、y₁、p₁和θ₁。当 cos ²β(tan |β|- sin φ)≤c cos φ时，基本方程如下：

式中：θ₁₂＝(θ₁＋θ₂)/2。当 cos ²β(tan|β|- sin φ) >c cos φ时，只需相应改变上式中的第4个方程即可。

设N₁为边界AB的结点，N₂为紧邻N₁的一个结点，且N₁N₂为α滑移线。已知N₂的x₂、y₂、p₂和θ₂，可采用迭代算法求解N₁的x₁、y₁、p₁和θ₁，基本方程如下：

式中：θ₁₂＝(θ₁＋θ₂)/2。

设N₁、N₂、N₃为求解域内3个相邻的结点，而且N₁N₂为α滑移线，N₁N₃为β滑移线。已知N₂的x₂、y₂、p₂和θ₂，以及N₃的x₃、y₃、p₃和θ₃，可采用迭代算法求解N₁的x₁、y₁、p₁和θ₁，基本方程如下：

式中：θ₁₂＝(θ₁＋θ₂)/2，θ₁₃＝(θ₁＋θ₃)/2。

采用有限差分法求解以上极限平衡边值问题，根据AB和AC应力边界条件的具体情况，可分为3种基本解法。现以图5为例分述如下：

第一种解法：AB边界各结点的x、y已知，AC边界各结点的p、θ已知。求解过程如下：

1)求解第二条α滑移线，需采用以下迭代计算：

①给结点6的平均应力 置初值，并置k＝0。

②采用3PD差分依次求解结点7至11。

③由2PDA差分求解结点12。

④比较 和按式(19)计算的p_N，若| -p_N|≤10^-8，则迭代计算终止，否则，确定高一次 ，并置k＋1给k，转第②步继续计算。

2)按以上同样的方法，可依次求解第三条和第四条α滑移线。

第二种解法：AC边界上各结点的x、y、p和θ已知。求解过程如下：

1)求解第二条α滑移线，采用3PD差分依次求解结点11至7，由2PDB差分求解结点6。

2)求解第三条α滑移线，采用3PD差分依次求解结点20至14，由2PDB差分求解结点13。

3)求解第四条α滑移线，采用3PD差分依次求解结点31至23，由2PDB差分求解结点22。

第三种解法：AC边界上各结点的x、y、p和θ已知，AB边界的p与θ关系已知。文献[29,30]采用分区解法依次求解A″CF″区的Cauchy问题、A′A″F″F′区的Riemann问题以及A′BF′区的混合问题，求解过程中反复调用3PD和2PDB差分。

应用表明，以上3种解法的计算结果是相同的，但各有优缺点。第一种解法的最后一条α滑移线总是通过墙踵点B，但求解每条α滑移线时需反复迭代，计算量较大，而且在某些情况下，例如当 cos ²β(tan|β|- sin φ)>c cos φ时，迭代计算对初值选取较为敏感而难以收敛。第二和第三种解法本质上是一致的，计算量较小且稳定性较好，但最后一条α滑移线一般不通过B点，需采用线性内插求得。本文计算机程序提供这3种算法。

按以上方法求解的塑性区一般由ABF′(Ⅰ区)、AF′F″(Ⅱ区)以及ACF″(Ⅲ区)三部分组成，图6给出几种常见的滑移线场：

图  6  几种滑移线场

Figure 6.  Several kinds of slip-line field

1)如图6(a)所示，三个区相互之间无面重叠，这种情况的Δ>π。

2)如图6(b)所示，Ⅱ区退化为一条β滑移线，三个区相互之间无面重叠，这种情况的Δ＝π。

3)当Δ<π时，滑移线场通常如图6(c)所示，滑移线场沿AF′和AF″出现折叠现象，即三个区相互之间有面重叠；同时应力场也出现折叠现象，即重叠区中同一点对应3个不同的应力状态；另外，作用于塑性区ABC上的所有外力不能满足静力平衡条件。但是，一些研究者^[16,19,23]认为这是一种应力间断现象，还给出了所谓的“应力间断条件”和“间断线”。

4)当c_w＝c≠0，δ<φ时，即使Δ≥π，也常会出现折叠现象，如图6(d)所示：滑移线场Ⅰ区自身折叠，Ⅰ区和Ⅱ区有面重叠，重叠区中同一点对应3个不同应力状态，静力平衡条件不能满足。

5)图6(e)是一种更为畸形的滑移线应力场：溢出边界，Ⅰ区自身面重叠，三个区有面重叠。

笔者认为，真实应力场中的点(x，y)与应力(p，θ)一一对应，且作用于塑性区上的所有外力必须满足静力平衡条件。因此，除图6(a)、图6(b)两种情况外，其它只是该边值问题在数学意义上的虚解而不符合实际情况。

对于图1所示的挡土墙主动土压力的极限平衡边值问题，当X＝0，Y＝γ>0时，可采用以下无量纲分析法求解。首先引进3个新的无量纲变量x′、y′和p′：

同时定义3个无量纲参数η₁、η₂和η₃：

滑移线方程可改写为以下无量纲形式：

沿α和β滑移线的极限平衡方程可改写为：

式中：R′＝p′ sin φ＋η₁ cos φ，而A′(θ)、B′(θ)按下式计算：

边界AB的应力边界条件可改写为：

式中：R′_M＝p′_M sin φ＋η₁ cos φ；R′_w＝p′_M sin δ＋η₂ cos δ。

边界AC的应力边界条件可改写为：

当|η₃|≤η₁ cos φ时，

当η₁ cos φ<|η₃|≤η₁时，

式中：R′_N＝p′_N sin φ＋η₁ cos φ。

应力奇点A的应力边界条件可改写为：

经变换后的式(40)(式(49)与原来相比有着相似的数学形式，它们构成了一个完整的无量纲应力边值问题，所确定的(x′，y′，p′，θ)称为规格化滑移线应力场，可采用前面的差分法求解。规格化滑移线应力场取决于参数α、β、φ、δ、η₁、η₂及η₃，当β＝0时，仅取决于α、φ、δ、η₁、η₂。根据式(37)和式(38)，规格化滑移线应力场可转化为实际滑移线应力场，进而按式(16)计算p_a，或直接按下式计算：

沿α滑移线和β滑移线的土反力分布可按下式计算：

实际应力场可由下式确定：

式中：σ′_x、σ′_y、τ′_xy为规格化正应力和剪应力，有：

根据以上无量纲分析，关于主动土压力极限平衡问题的几何力学相似原理可表述为：若两个问题的α、β、φ、δ、η₁、η₂及η₃相等，则滑移线场几何相似，相似系数为两者(H-h₀)之比值；对应于相同规格化坐标的两点的剪应力相似，正应力减去 cos ²β后相似，主动土压力强度减去 后相似，相似系数为两者γ(H-h₀)之比值；总主动土压力减去 后相似，相似系数为两者γ(H-h₀)²之比值。不难证明，其逆命题也成立。

设 以及q′＝ ，还可进一步证明以下3个命题：

推论1：若两个问题的α、β、φ、δ、c′、c′_w、q′相等，则它们几何力学相似，而且正应力、主动土压力强度、总主动土压力分别相似。除β＝0，q> 的情况外，逆命题也成立。

推论2：若两个问题的 ，则它们几何力学相似的充要条件是α、β、φ、δ、c′、c′_w相等。

推论3：若两个问题的α、β、φ、δ、η₁、η₂及η₃相等，则它们的奇异性相同即Δ相等。

当体力X≠0，Y≠0时，也可采用无量纲分析法求解；对于无重土X＝Y＝0的情况，只需将以上各式中的γ取任意非零值如水容重等以及将式(42)、式(43)右端项取0即可。

下面给出几个特例的一些数值分析结果：

情况1：当γ＝0， ，Δ≥π时，有以下一些结论：

1)主动土压力强度处处相等，除了几种特殊情况有解析解外，一般需由迭代计算求得。

2)Ⅰ区或Ⅲ区为均匀应力场，即p、θ处处相等，α滑移线和β滑移线为直线。

3)Ⅱ区的β滑移线为直线且相交于A点，每条β滑移线上的p、θ值不变，α滑移线为对数螺旋线，中心为A点，极坐标方程为：

式中：r₀、θ₀分别为初始极半径和大主应力方向角。特别地，当Δ＝π时，Ⅱ区退化为一条β滑移线。

4)Hencky第一定理适用。

5)Hencky第二定理不完全适用，仅限于Ⅱ区沿β滑移线有：

沿α滑移线因R_β处处无穷大、 在数学上无意义而没有类似的关系式。

情况2：当δ＝φ＝0，β＝0，Δ≥π时，有以下一些结论：

1)主动土压力呈线性分布。

2)Ⅰ区或Ⅲ区为线性应力场，即任意两点之间p呈线性变化而θ处处相等，α滑移线和β滑移线为直线。

3)Ⅱ区的β滑移线为直线且相交于A点，每条β滑移线上的p呈线性分布而θ不变，α滑移线为圆弧，圆心为A点。特别地，当Δ＝π时，Ⅱ区退化为一条β滑移线。

4)Hencky第一定理适用。

5)Hencky第二定理不完全适用，仅限于Ⅱ区沿β滑移线满足式(55)的关系式。

应当指出，若β≠0，则除第4点结论仍成立外，其它几点均不成立。

情况3：当α＝π/2，β＝0，c_w＝0，δ＝0时，或当c_w＝c＝0且Δ＝π时，本文滑移线解与朗肯解或经典库仑解完全一致，有以下结论：

1)主动土压力呈线性分布。

2)Ⅰ区或Ⅲ区为线性应力场，α滑移线和β滑移线为直线。

3)Ⅱ区退化为一条β滑移线。

4)Hencky第一定理适用。

5)Hencky第二定理不适用。

情况4：当α＝π/4＋φ/2，β＝0，c_w＝c，δ＝φ时，其滑移线解与文献[8]的库仑解完全一致，有以下结论：

1)主动土压力呈线性分布，与滑移线BC的土反力分布相互对称。

2)Ⅰ区和Ⅱ区退化为边界AB，此时AB为一条β滑移线。

3)Ⅲ区为线性应力场，α滑移线和β滑移线为直线。

4)Hencky第一定理适用。

5)Hencky第二定理不适用。

某挡土墙H＝10 m，c_w＝c＝0，φ＝30°，其他参数γ、q、α、β、δ以及计算结果详见表1。按本文方法计算的总主动土压力与Sokolovskii(1965)的相比，当δ＝0时两者完全一致，而当δ>0时有些偏差，这是因为Sokolovskii(1965)主动土压力计算公式的右端项分母被人为地忽略了一个cos δ。本文滑移线解大于或等于经典库仑解，当Δ＝π时两者完全一致。

某挡土墙H＝10m，q＝20 kPa，α＝95°，β＝10°，γ＝19 kN/m³，c＝15 kPa，φ＝25°，c_w＝10 kPa，δ＝10°。采用第一种解法，计算结果如下：

1)差分网格划分取n＝49，m＝17，结点总数为3 185个，AB边界上结点等间距布置。计算得到 ＝64.816 kPa，h₀＝2.396 m，Δ＝211.384°，P_a＝211.975 kN，主动土压力分布以及滑裂面BC上的土反力分布如图7(a)所示。为了图形清晰起见，图中给出的滑移线场网格为n＝25，m＝7。

图  7  算例2计算结果：(a)本文滑移线解；(b)库仑解

Figure 7.  Results for Example 2: (a)slip-line solution in this paper; (b)Coulomb′s solution

2)文献[8]假定一簇α滑移线为直线即平面滑裂面，采用极限平衡法推导了黏性土情况的库仑解。其计算结果为 ＝68.391 kPa，h₀＝2.587 m，P_a＝203.700 kN，主动土压力分布以及滑裂面BC的土反力分布如图7(b)所示，图中给出与图7(a)相同深度处的土压力强度值，总库仑主动土压力比滑移线解小3.90%。

3)作用于图7(a)塑性区上所有外力在水平、垂直方向的投影之和F_x＝-1 N，F_y＝-6 N，对B点的力矩之和M_B＝-8 N·m，可见，这些外力满足了平面力系的3个平衡条件。

4)进一步数值分析表明，Hencky第一和第二定理在这里均不成立。

1)通过引进应力奇点及其应力边界条件，极限平衡问题静定可解而不必考虑应力-应变关系。

2)塑性区内处处满足静力平衡方程和屈服条件，应力处处连续，作用于塑性区上的所有外力满足静力平衡条件。

3)主动土压力的滑移线解一般总是大于或等于库仑解，朗肯解或满足非奇异条件的经典库仑解与滑移线解一致。

4)两个问题几何力学相似的充要条件是参数α、β、φ、δ、η₁、η₂及η₃相等。

5)Hencky第一定理和第二定理不具有普遍适用性。

6)极限平衡问题是一个理想化问题，而滑移线解是该问题的理论解，它与真实状态下的极限荷载有多大相符或差异，尚有待实验和实践去检验。