水位变化引起的瞬态渗流是影响岩土边坡安全性的一个重要因素，也是一个典型的流固耦合问题。据统计，超半数的河流、水库边坡失稳与水位变化相关^[1-2]。国内学术界对在边坡稳定问题中如何考虑渗流的影响存在不同看法：有的学者认为在土坡稳定分析中应当将渗透力作为一种体力加以考虑^[3-4]，也有的学者认为在土坡稳定问题中考虑渗透力属于化简为繁^[5-8]、概念不清^[9]，双方为此争论了30余年^[10-11]。这一争论的核心在于边坡稳定问题的研究对象是土骨架还是土水混合体。目前普遍被接受的观点认为：若将土水混合体作为研究对象，则渗透力不应被添加到分析中；若将土骨架作为研究对象，则应当考虑渗透力对土骨架的作用。这两种研究方法在力学上是等价的，对于同一问题应当得出相同的安全系数^[12]。

关于水位下降对边坡稳定性的影响已有一些成果，但国内外学者所采用的分析方法各不相同，得到的结论也并不统一。相关研究大部分以分析浸润线^[13-16]或孔压^[17-20]对边坡稳定性的影响为主，较少考虑渗透力的影响。笔者将土骨架作为研究对象，将渗透力考虑为外荷载作用于土骨架上，采用强度折减法^[13]，对坡面水位下降过程中的边坡稳定问题进行研究，主要讨论不同渗透系数K、水位降幅Δh与水位降速v对非饱和土质边坡安全稳定的影响。

1 渗透力的定义与计算

如何定义渗透力是岩土力学界广泛讨论的一个问题，国内外不少学者对渗透力的定义、力学性质、组成部分、计算方法等问题发表过不同看法^{[9-11, 21-26]}。在很多文献和土力学教材中，渗流施加于土骨架上的拖曳力被定义为渗透力^{[22, 25-26]}。这种定义并不准确。以圆球形的土颗粒为例，水流经球形土颗粒类似于流体力学中黏性不可压缩流体绕球缓慢流动问题，如图 1所示。斯托克斯于1851年给出了蠕动流绕球流动阻力D的解析解：

$ D = 6\pi \mu aU, $

(1)

式中：μ为流体的动力黏度；a为球体半径；U为流体的流速。值得注意的是，在斯托克斯阻力公式中，蠕动流作用于圆球颗粒表面的力D是由作用于颗粒表面切向的摩擦阻力4πμaU和作用于颗粒表面法向的压强阻力2πμaU两部分组成。在岩土渗流问题中，若仅将渗流作用于颗粒表面切向的拖曳力定义为渗透力，则忽视了渗流对颗粒表面法向的压力作用。“渗流作用于土骨架上的推动力和拖曳力合称为渗透力”^{[21, 23-24]}的定义更为明确和完善。对于微观土颗粒而言，渗透力是作用于颗粒表面的面力；对于宏观土体而言，渗透力可以被当成体力处理。

对于渗透力的计算，现有文献给出的数学表达式较为一致。单位体积内的土骨架受到的渗透力为：

$ \mathit{\boldsymbol{j}} = {\gamma _{\rm{w}}}\mathit{\boldsymbol{i,}} $

(2)

式中：γ_w为水的容重，N/m³；i为水力梯度，其分量表达式为：

$ {{i_n} = \partial H/\partial {x_n}(n = 1,2,3),} $

(3)

$ {H = (z - {z_0}) + \frac{{{p_{\rm{w}}}}}{{{\gamma _{\rm{w}}}}} + \frac{{v_{\rm{w}}^2}}{{2g}},} $

(4)

式中；H为总水头；z为测压点高度；z₀为零势面的高度；p_w为孔隙水压；v_w为孔隙水流速；g为重力加速度。总水头由位置水头、压强水头和流速水头三部分组成。在岩土渗流问题中，v_w通常是个小量，流速水头v_w²/2g可忽略不计，因此总水头可近似表示为：

$ H = (z - {z_0}) + \frac{{{p_{\rm{w}}}}}{{{\gamma _{\rm{w}}}}}。$

(5)

2 土体平衡微分方程与强度理论 2.1 平衡微分方程

以土骨架为研究对象，考虑渗流条件下的平衡方程，以有效应力表示平衡条件，根据Biot固结理论^[27-28]：

$ \sigma _{ij,j}^\prime + {\gamma _{\rm{w}}}\frac{{\partial {\kern 1pt} H}}{{\partial {\kern 1pt} {x_i}}} + b_i^\prime = 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (i,j = 1,2,3), $

(6)

式中：σ′_{ij, j}为有效应力张量；b′_i为对应于浮容重的单位体积土体体力，可显式表达为b′_i=0, 0, γ_w－γ_sat，γ_sat为土体的饱和容重。由式(5)和(6)可得：

$ \sigma _{ij,j}^\prime + \frac{{\partial {p_{\rm{w}}}}}{{\partial {x_i}}} + {b_i} = 0, $

(7)

式中：b_i为对应于饱和容重的单位体积土体体力，可显式表达为b_i=(0, 0, －γ_sat)。式(6)和(7)两种平衡方程表述方式是等价的，仅在选取水头种类和土体容重方面有所区别。

2.2 非饱和土的有效应力原理及强度准则

饱和土中只存在固、液两相，其有效应力可以表述为：

$ \sigma _{ij}^\prime = {\sigma _{ij}} - {u_{\rm{w}}}, $

(8)

式中：σ_ij为土体中的总应力；u_w为孔隙水压力。考虑到非饱和土中存在固、液、气三相，Bishop基于式(8)推广出了非饱和土的有效应力公式^[29]：

$ \sigma _{ij}^\prime = {\sigma _{ij}} - [\chi {u_{\rm{w}}} + (1 - \chi ){u_{\rm{a}}}]{\delta _{ij}} = {\bar \sigma _{ij}} + \chi s{\delta _{ij}}, $

(9)

式中：u_a为土体中的气压；δ_ij为Kronecker记号；$ \bar{\sigma}_{i j}=\sigma_{i j}-u_{\mathrm{a}} \delta_{i j}$为净应力；s=u_a－u_w为吸力。χ为一个介于0与1之间的与土体饱和度相关的变量，在实际应用中χ取值通常为有效饱和度S_e。当土体处于饱和状态，χ的值取为1，此时式(9)退化为式(8)。

Mohr-Coulomb强度准则是岩土工程中常用的强度条件。Khalili等^[30]、Wheeler等^[31]及Cui等^[32]的工作表明，通过引入Bishop有效应力，M-C强度条件可以同时描述土体在饱和与非饱和状态下的强度：

$ \tau = {c^\prime } + {\sigma ^\prime }{\kern 1pt} {\rm{tan}}{\kern 1pt} {\kern 1pt} {\varphi ^\prime } = ({c^\prime } + {S_{\rm{e}}}s{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{tan}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi ^\prime }) + \bar \sigma {\kern 1pt} {\kern 1pt} {\rm{tan}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi ^\prime }, $

(10)

式中：τ和$\bar \sigma $分别表示作用于土体剪切破坏面上的剪应力和净应力；c′和φ′为土体在饱和状态下的黏聚力和摩擦角。式(10)中的(c′+S_estan φ′)可以理解为非饱和土体的表观黏聚力，它一般随土体基质吸力增大(对应饱和度的下降)而增大。基质吸力与饱和度之间的关系可以由非饱和土的土水特征曲线来描述。

3 土水特征曲线

非饱和土的强度受饱和度或吸力变化的影响，饱和度与吸力之间关系由土水特征曲线(soil-water characteristic curve, SWCC)描述。笔者主要关注水位下降过程中的非饱和土干化过程，相关的饱和度吸力关系可采用Brooks & Corey^[33]提出的用有效饱和度S_e表示的土水特征曲线经验公式。该公式较为简单实用，并与常规的试验数据也较吻合，具体公式为

$ {S_{\rm{e}}} = \frac{{{S_{\rm{r}}} - {S_{{\rm{rr}}}}}}{{1 - {S_{{\rm{rr}}}}}} = \left\{ {\begin{array}{*{20}{c}} {1\quad }&{s \le {s_{\rm{b}}},}\\ {{{({s_{\rm{b}}}/s)}^\lambda }}&{s > {s_{\rm{b}}}}。\end{array}} \right. $

(11)

式中：s_b为进气值，代表空气开始进入土体时对应的吸力；S_r为土体饱和度；S_rr为土体的残余饱和度；λ为反映土体孔隙大小分布的参数，定义为有效饱和度和基质吸力曲线的负斜率，其值小于1，且孔隙尺寸分布越均匀则λ的值越大。图 2为用式(11)描述的上海黏土SWCC，其进气值s_b和λ的取值分别为74 kPa和0.5。

4 水位下降过程的安全系数

在边坡水位下降过程中，超孔隙水压的消散总是滞后于坡面水位的下降，浸润线呈现上凸的形态。坡面水位下降中产生的渗透力对边坡稳定具有消极影响。但水位下降也会提高土体内部的有效应力；土体中的非饱和区域也会随着浸润线的降低而增大，非饱和土基质吸力作用增强，对边坡稳定产生积极影响。两种因素的共同作用决定了非饱和土质边坡稳定性的变化。因此，以强度折减因子作为评估边坡稳定性的安全系数，通过模拟坡面水位下降的瞬态渗流过程，计算水位下降过程中边坡的瞬时安全系数F_s，进而分析边坡稳定性的变化规律。

水位下降过程中的瞬态渗流场是一个非稳定渗流场，孔压的消散、饱和度的降低都与时间相关。对于瞬态条件下的非饱和土质边坡稳定性问题，常用的有限元计算软件通常不支持直接通过强度折减法计算安全系数。笔者采用的分析方法为：根据边坡水位下降过程的不同工况，结合式(7)(9)和(11)对坡面水位下降进行流固耦合分析，获得土体在该工况下瞬时的孔隙水压、饱和度和有效应力等数据，然后结合式(2)计算出等效作用于各个结点的瞬态渗透力，进而结合土体平衡微分方程(7)和M-C强度准则对土骨架进行强度折减计算，最终得到该瞬态渗流下土坡的瞬时安全系数F_s。在一个坡面降水位的过程中，每一个时刻都对应了一个F_s，这些瞬时安全系数中的最小值F_{s, min}则对应了该降水位过程中最危险的状态，即为整个降水位过程的安全系数。

5 有限元模拟与分析

笔者使用有限元计算软件ABAQUS对坡面水位下降过程中的稳定性问题进行数值模拟，数值模型如图 3所示：右侧为与土接触端面，高度为35 m，水位高32 m且在模拟过程中保持不变；左侧为变水位坡面，初始水位与右侧相同，在模拟坡面降水位的过程中下降至20 m。采用平面四结点单元，模型共3 900个单元，4 041个结点。模型为均质边坡，采用上海黏土的参数，干密度为1 450 kg/m³，弹性模量为7.0 MPa，泊松比为0.375，孔隙比为1.0，黏聚力为90 kPa，内摩擦角为26°，SWCC如图 2所示。在初始水位，土坡的初始安全系数F_s0为1.556。

渗透系数K、水位降幅Δh与水位降速v是影响边坡稳定的重要因素。渗透系数决定了超孔隙水压力的消散速度，水位下降幅度和速度将对土体饱和度的分布以及水力梯度的变化产生重要影响，这3个因素共同决定了水位下降对边坡稳定是否有利。数值模拟将采用控制变量的方法对几组不同的工况进行模拟，分析这3个因素在土坡稳定中所起到的作用。

黏土的渗透系数K通常较小，分别选取1.0×10^-8 m/s、5.0×10^-8 m/s和1.0×10^-7 m/s 3组渗透系数进行模拟。在实际的水库边坡降水位问题中，水位下降速度v通常在1.0 m/d至2.0 m/d之间。为了更清楚地研究水位下降速度对边坡安全的影响，将下降速度的模拟范围扩大至5.0 m/d。不同的水位降速可以通过设置孔压边界条件获得。在分析中，当土体强度参数折减到塑性区贯通或塑性应变过大导致计算无法收敛时，即得到相应的安全系数。图 4为对应于渗透系数1.0×10^-7 m/s、水位降速1.0 m/d、水位降幅12 m时的土坡在强度折减到临界状态时的等效塑性应变云图，从图上可以看出，一个滑移面已经在坡内生成，坡脚附近出现了较大的塑性应变。

第1种工况，将模型的渗透系数K指定为1.0×10^-8 m/s，研究坡面水位下降过程中，非饱和土质边坡的瞬时安全系数F_s随水位降幅Δh与水位降速v的变化规律。在不同的水位降速下，F_s随时间的变化情况如图 5所示。可以看出，随着坡面降水位过程的推进，F_s呈现出先降后升的趋势。在水位下降的初期，F_s较初始时刻产生了一个明显的下降，在达到安全系数最小值F_{s, min}后随着坡面水位的继续下降逐渐上升。F_{s, min}的大小与水位降速v相关，v越大则F_{s, min}的值越小，到达该最小值点的时间也越短。此外，在相同的水位降幅下，水位降速v越小的土坡其安全系数越大，且水位下降越慢这一特点越明显。这两点都表明水位下降速度越快对土坡的安全越不利。

F_s的上升主要集中在水位降幅4~8 m的阶段，并逐渐越过初始安全系数线。在水位降速v较大的算例中，当Δh超过8 m后，F_s上升的趋势明显减缓。水位降速v更大的土坡，而不是更小的土坡，F_s升至F_s0上方用时更短。这一现象可以解释为：在相同时间内，虽然更大的水位降速v会产生更大的坡面水位降幅，在土坡内产生更高水平的渗透力，但更大的水位降幅也会导致土体内的水渗出更快，使有效应力的增加、非饱和区域的扩大、饱和度的降低都变得更显著，当有利影响居于主导地位时就会出现这一现象。这说明在降水位过程中，边坡的稳定性是由渗透力和土坡强度的提高这两个因素共同决定的。

图 6记录了在第1种工况下坡脚附近的第756结点(即图 4中的高亮结点)所受到渗透力的水平分量F_jx。在相同的时间内，水位降速大的土坡内部的渗透力水平普遍高于水位降速慢的土坡。在3个算例中都出现了756结点的F_jx在曲线右端下降的现象，这是因为在Δh接近12 m的时候，边坡左侧的水位逐渐下降至756结点之下，756结点附近的水力梯度和渗流场都发生了显著的变化，从而导致了曲线右端F_jx的下降。此外，从图 6还能看出，在相同的水位降幅下，水位降速v越大则渗透力水平越高。通过对比图 5和图 6可以发现，渗透力水平最高点出现的时间与F_{s, min}点出现的时间并不重合，这再次说明坡面降水位引起的渗透力并不是决定土坡安全的唯一因素。

第2种工况，指定坡面水位降速v为1 m/d，研究土坡的瞬时安全系数F_s随渗透系数K和水位降幅Δh的变化规律，计算结果如图 7所示。不同渗透系数的计算结果都表明，随着坡面水位的下降，F_s先降后升，F_{s, min}出现在水位下降的早期，这一特征与上一个工况相同。渗透系数高的土坡，其F_{s, min}的值也较大，出现的时间也较早，F_s从F_{s, min}升至F_s0上方所用时间也更短。从图 7中还能看出，K越小对土坡的安全越不利，K越大的土坡对坡面水位下降带来F_s提升的响应也越快。

图 8记录了该工况下756结点渗透力水平分量F_jx的变化曲线，从图中可以看出，在相同的水位降幅下，渗透系数小的土坡其内部的渗透力水平更高；随着坡面水位的下降，渗透力水平也基本保持提高的趋势，只在曲线右端出现了下降或增长速度的降低，这与前一工况相类似。

在以上2种工况的分析中，非饱和土质边坡的F_s随着坡面水位的下降都表现出了先降后升的特点，这表明坡面降水位的初期是整个降水位过程中最危险的阶段。第3种工况考虑不同的渗透系数K及水位降速v与边坡安全系数最小值F_{s, min}之间的关系，计算结果如图 9所示。可以看出，F_{s, min}随水位降速v的增大而降低，随渗透系数K的增大而升高。对于渗透性较好的土体，F_{s, min}随水位降速v变化的幅度并不十分明显；而对于渗透性较差的土体，F_{s, min}对v的变化较为显著。此外，水位降速较低时，不同渗透系数土坡的F_{s, min}差距并不很大，但随着水位降速的提高，不同渗透性土坡F_{s, min}的差距变得相当可观。在任一坡面降水位条件下，F_{s, min}始终小于初始安全系数F_s0，这表明坡面降水位的过程对土坡安全的影响是不利的，在处理实际问题中须要加以重视。

图 10记录了对应于不同条件下土坡最危险状态的756结点渗透力水平分量。与图 9对比可以发现，在土坡最危险状态时，F_{s, min}的大小与其内部的渗透力水平的高低呈负相关关系。对于渗透性较好的土坡，这种负相关关系较弱；而对于渗透性较差的土体，这种关系就显得比较明显。可见对于渗透系数较低的土坡，进行坡面降水位操作时需要着重控制水位降速以尽量降低土坡内部的渗透力水平，防止土坡安全系数下降过大；对渗透系数较高的土坡，可供其选择的水位降速区间相对更宽一些。

6 结论

笔者主要研究了考虑瞬态渗流的非饱和土质边坡在水位下降过程中的安全问题，提出了一种边坡降水位条件下考虑渗透力的非饱和土质边坡安全系数有限元计算方法。通过数值模拟发现，在坡面降水位过程中，土坡的安全系数是由土体内部的渗透力等不利因素和土体强度提高等有利因素共同作用的结果，处于主导地位的一方决定了土坡是否安全。数值模拟结果表明：在坡面水位下降的初始阶段，土坡中的孔隙水压力尚未有足够的时间消散，土坡中的有效应力变化不明显，渗透力在此阶段处于主导地位，使边坡的稳定性有所下降；在水位下降的后期，随着孔隙压力的逐渐消散，土坡中的有效应力增大成为影响边坡稳定的主要因素，边坡的稳定性提升。体现边坡稳定性的安全系数(强度折减因子)F_s呈现出先下降后上升的变化趋势。

在坡面降水位的过程中，存在一个安全系数的最小值F_{s, min}，其大小、出现时刻与土体渗透系数K和水位降速v相关。坡面降水过程中土坡的稳定性受渗透系数和水位降速的共同影响。F_{s, min}的大小随K的增大而增大，随v的增大而减小；出现的时刻随K的增大而减小，也随v的增大而减小。对于渗透系数较低的土坡，随着水位降速的增大，F_{s, min}会出现明显的下降；对于渗透系数较高的土坡，F_{s, min}随着水位降速增大的变化相对不明显。渗透系数较低土坡的安全系数从F_{s, min}恢复到初始水平需要耗费更长的时间。

在面对水库降水位等实际问题时，应着重关注水位下降初始阶段的边坡安全。降水位初期是边坡最危险的阶段，特别是对于渗透系数较小的土质边坡。合理控制水位下降速度，给超孔隙水压力的消散预留更充裕的时间，可以起到降低土体渗透力水平，增大非饱和区域面积，降低土体饱和度的效果，对于提升边坡整体强度、防止边坡破坏具有重要意义。