水位上升对土坡稳定性的影响

柴学锐，孙冬梅，S. Semprich

(1.天津大学 水利工程仿真与安全国家重点实验室，天津 300072；2.格拉茨技术大学 土力学及基础工程研究所，奥地利 格拉茨 A-8010)

坡外部水位升高，岸坡土体内的地下水位也随之发生改变，原来被气占据的孔隙空间逐渐被水所补充，在这种“水驱替气”的过程中，孔隙压力的变化将影响岸坡的稳定性。因此，研究水位上升作用下岸坡的渗流场变化情况对进一步研究岸坡稳定性变化规律至关重要[1-2]。

水位上升作用下的岸坡是一种常见的饱和-非饱和土坡。在以往的饱和-非饱和渗流问题中，一般仅研究水相而忽略气相流动，尽管这种仅考虑水相的方法简单易行并且已经积累了一定经验，但得出的结论与实际工况相比有很大出入。因此，为了更好地模拟非饱和区流体的运移规律，国内外有学者认为采用水-气二相流模型分析饱和-非饱和渗流问题更符合实际[3-7]。

基于水-气二相流模型，分析了水位上升过程中，饱和-非饱和土坡内渗流场的变化规律；采用Bishop简化法[8]推导出同时考虑孔隙气压力、毛细压力和坡外水压力等贡献的土坡安全系数计算公式，并计算了不同时刻边坡的安全系数；分析安全系数的变化规律及各个力对安全系数的贡献，为进一步研究水位上升作用下土坡的失稳机理提供理论依据。

1 水-气二相流模型

近年来，多孔介质的多相流理论研究和应用越来越受到学者重视[9]。对于岩土工程渗流问题而言，通常认为渗流系统绝对恒温，土体骨架不变形，处于平衡状态的水相(l)和气相(g)在多孔介质中流动，且满足Darcy定律。

1.1 基本控制方程

当非饱和土体中流体的扩散作用忽略不计时，饱和-非饱和水-气二相流模型的基本控制方程为

(1)

1.2 本构方程

1.2.1 毛细压力-饱和度关系模型

近年来，许多学者均提出了毛细压力和饱和度的关系模型，如Brooks and Corey(BC)模型[10]和Van Genuchten(VG)模型[11]等，其中BC模型的进气值不具备连续性，而Van Genuchten(VG)模型不仅其表达形式简单，且模型的关系曲线中进气值和残余水饱和度处过度平滑，因此在水-气二相流模型中Van Genuchten(VG)模型得到广泛应用

(2)

式中：pc为水相和气相交界面处的毛细压力，等于孔隙水压力和孔隙气压力的差值；P0为进气值；参数Sle为有效水饱和度，Sle=(Sl-Slr)/(Sls-Slr)，其中Sl为水饱和度，Slr为残余水饱和度，Sls为饱和水饱和度；λ、ζ为模型的拟合参数，其值与土壤自身性质有关；Pmax为最大毛细压力。

1.2.2 相对渗透率和饱和度关系

土体孔隙中同时存在气相流体和水相流体，一种流体流动将会干扰另一种流体的流动，当水饱和度变化时，水相和气相之间的交界面也将随之变化，水流通道的弯曲度将会改变，导致相对渗透率发生改变。1976年，Mualem提出了相对渗透率和饱和度的关系表达式，将反应毛细压力和饱和度关系的VG模型表达式带入相对渗透率和饱和度的关系表达式中，可得到TOUGH2中相对渗透率和饱和度关系的改进模型，即VGM模型[12]

(3)

式中：krl为水相相对渗透率；krg为气相相对渗透率；τ为迂曲因子，其取值大小与土体孔隙弯曲度有关，这里取0.5。

1.3 边界条件处理方法

边界条件实际上是指研究区域边界所处的条件。在数学模型中，通过在渗流场边界施加一薄层虚拟单元来实现各种边界条件的数值模拟。其中，第一类边界条件(Dirichlet条件)和第二类边界条件(Neumann条件)为水-气二相流模型根据液相和气相划分的边界条件。

1.3.1 第一类边界条件

1.3.2 第二类边界条件

第二类边界条件主要用来描述系统与外界的流量交换情况，且第二类条件可以是常量，也可以是随时间变化的量。将上升引起的加载作用等效为边界单元上有一个流入的流量，该值可通过式(4)求解

(4)

式中：m为源汇项，流入为正；b为拟合系数，拟合系数的取值与土壤的性质有关；V为边界单元的体积；C为边界单元的压缩性；ΔP为水压力的变化量；Δt为水压力变化所需的时间。以往的研究通常将水位变化过程等效为一系列特征水位，这种处理方法不能很好地模拟边界上相态的转变，而将水位上升过程等效为边界单元上有一个流入的流量，可以实现水位的连续变化和边界上相态的转变。

1.4 初始条件

初始条件即为t=0时刻土体中各个状态变量的分布情况。求解土体稳定渗流问题时，只需给出相应的边界条件即可，不需给定初始条件，但是求解土体非稳定渗流问题时，除了给出边界条件之外，还需给定初始条件。对于坡外水位上升作用下的非稳定渗流问题，迭代求解的时间取决于给定的初始条件的合理性，当初始条件接近真实值时，迭代收敛速度会加快，所需要的模拟时间会减少，因此初始条件的给出应该以接近实际物理意义为准则。

2 非饱和土坡稳定性分析方法

2.1 非饱和土坡抗剪强度公式

Fredlund等[13]认为非饱和土的抗剪强度由净法向应力、基质吸力和有效黏聚力组成，由此提出了非饱和土的抗剪强度公式

τ=c′+(σn-pg)tanΦ′+(pg-pw)tanΦb

(5)

图1 Bishop简化法下的土条作用力Fig.1 Forces acting on a slice in simplified Bishop′s Method

式中：τ为非饱和土的抗剪强度；c′为有效黏聚力；σn为总法向应力；pg=pg-patm为孔隙气压力；pw=pw-patm为孔隙水压力；Φ′表示有效内摩擦角；Φb表示抗剪强度随基质吸力增加的曲线倾角；(σn-pg)为净法向应力；(pg-pw)为基质吸力，其值等于负的毛细压力[14]。

2.2 非饱和土坡安全系数计算公式

Bishop简化法常用来计算某一危险圆弧滑动面上的安全系数。图1为一均质非饱和土坡，土坡中任一土条i上的作用力如图所示：bi为土条i的宽度；αi为土条i的底面与水平方向的夹角；βi为土坡侧面与水平方向的夹角；Wi为土条i的重力；Ni为土条i底部的正法向应力；Ti为土条i底部的切向力；hi为土条i的高度；pi为作用在土条i上的坡外水压力。

根据土条i在垂直方向上受力平衡

(6)

根据圆心力矩平衡

(7)

式中：di为坡外水压力pi到圆心的垂直距离，R为半径。

根据极限平衡条件，土条i底面的抗剪力Ti可表示为

(8)

则非饱和土坡圆弧滑动面上的安全系数F可表示为

(9)

3 算例分析

3.1 算例说明

图2 边坡几何形状剖面图(单位：m)Fig.2 Cross section of soil slope

图2为一个各向同性的匀质土坡，该土坡长58.0 m，宽30.0 m，左侧坡高25.0 m，右侧坡高10.0 m，坡面倾角40°，坡内坡外的水位水平，且水位高程均为10.0 m。为了保证模拟结果的精确程度，对土坡进行网格剖分，若网格尺寸剖分过小，则计算时间较长且计算结果精度并未有明显提升，若网格尺寸剖分过大，将导致模拟结果精确程度降低，因此，根据笔者已有经验将网格剖分成水平方向上从0.5 m到1.0 m不等，垂直方向的网格尺寸均为0.5 m。在土体处于完全饱和的状态下，采用Bishop简化方法计算得到最危险圆弧滑动面的圆心坐标为(43.0，31.1)，半径为21.1 m，且其中一条半径与水平方向垂直。在该土坡上设置一个追踪面A-A，该追踪面处于土坡的中间偏左位置，其水平坐标是25.0 m。

表1 土体特性参数Tab.1 Parameters of different soils

整个模拟过程中，设定研究区处于恒温系统中，表1给出土体的特性参数取值。

3.2 渗流场分析

初始稳态的边界条件包括地表(如图2中ab边所示)和坡外水位以上的右侧边界(如图2中bc边所示)为大气边界，坡外水位以下(如图2中cd和de边所示)为水相边界，区域左侧(如图2中af边所示)和底部(如图2中ef边所示)为不透水边界，在此边界条件下，运行模型直至毛细压力与重力平衡。水位上升过程中，坡外水位在52 h内从10.0 m快速上升至21.2 m，上升速度为0.598×10-4m/s，该过程中，地表面和坡外水位以上的右侧边界仍为大气边界，区域左侧和底部仍为不透水边界，但原来在水位以上的右侧气相边界随着水位上升逐渐被水淹没，转变为水相边界。

模拟过程中，在边界单元上施加源汇项，以实现由于水位上升引起的加载作用导致的边界单元上的水压力变化和相态的转变。初始的源汇项值可以通过式(4)计算得到，经过试算，得出水位上升过程的源汇项值为2.70×10-44kg/s。在边界单元上施加源汇项处理水位上升引起的加载作用可以实现边界单元上水压力的连续变化，且比以往将水位上升过程简化为一系列不连续的特征水位更符合实际。

在稳定渗流状态下，各个物理量的分布状态如图3所示，非饱区的孔隙气压力水头为0，孔隙水压力水头小于0，因此毛细压力水头小于0，在饱和区，孔隙气压力水头等于孔隙水压力水头，毛细压力水头为0。观察图4，随着坡外水位上升，坡外水逐渐流入土坡内部，原来被气体填充的土壤孔隙逐渐被水填充，水饱和度逐渐增加，饱和区范围逐渐上移，且毛细压力水头，孔隙水压力水头和孔隙气压力水头均逐渐增加，其中接近土坡表面处的各物理量的变化更为显著。由图5表示的水位稳定后渗流场各个物理量分布情况可知，饱和区范围继续上移，非饱和区的水饱和度继续增加，毛细压力水头、孔隙水压力水头和孔隙气压力水头也进一步增加。

3-a 水饱和度3-b 毛细压力3-c 孔隙水压力3-d 孔隙气压力图3 初始稳态渗流场状态分布Fig.3 Distributions of water saturation, capillary pressure and pore pressures in initial steady state

4-a 水饱和度4-b 毛细压力4-c 孔隙水压力4-d 孔隙气压力图4 52 h渗流场状态分布Fig.4 Distributions of water saturation, capillary pressure and pore pressures in 52 h

图6为追踪面A-A在初始稳态，第52小时和第200小时的孔隙气压力、孔隙水压力、毛细压力和水饱和度分布情况。观察图6，水位上升以前，非饱和土体内的孔隙气压力为0，但是由于毛细压力的存在，在饱和土和非饱和土的接触面上孔隙气压力有一个微小的波动。水位上升过程中，随着坡外水体流入土坡内部，孔隙气压力的大小逐渐增加，水位稳定后，坡外水体继续流入土坡内部，土体内的孔隙气压力继续增加。初始状态下，土体内的孔隙水压力水头随着高程的降低呈现线性增加的趋势，在水位上升过程中和水位稳定后，随着坡外水体流入土坡内部，孔隙水压力水头均在逐渐增加。毛细压力为孔隙水压力与孔隙气压力的差值，且在水位上升作用下，毛细压力值逐渐增加。水位上升前，土坡的水饱和度随高程的降低逐渐升高，在水位上升过程中和水位稳定后，由于坡外水流入土坡内部，土坡的水饱和度继续增加。

图6 初始稳态,第52小时和第200小时追踪面A-A上的孔隙气压力、孔隙水压力、毛细压力和水饱和度分布情况Fig.6 Simulation results of pore-air pressure, pore-water pressure, capillary pressure and water saturation at section A-A at initial steady state, 52 h and 200 h

3.3 饱和-非饱和土坡稳定性分析

图7为土坡安全系数随时间变化过程，同时考虑孔隙气压力、毛细压力、坡外水压力和土体重力贡献的安全系数的初始值为1.54，随着水位的上升，安全系数逐渐增加，当水位升至21.2 m时，安全系增加至最大值2.67，水位稳定后，安全系数开始逐渐下降。不考虑孔隙气压力贡献的安全系数明显大于考虑孔隙气压力和毛细压力贡献的安全系数，其值从初始时的1.54上升到2.86后开始逐渐下降。不考虑毛细压力贡献的安全系数在三者中最小，其值从初始的1.06上升至2.21后开始逐渐下降。水位上升过程中，对土坡稳定性起主要作用的且有利于土坡稳定的坡外水压力逐渐增加，使得安全系数逐渐增加；当水位稳定后，坡外水压力保持不变，而不利于土坡稳定的孔隙气压力值逐渐增加使得安全系数逐渐降低。

图7 危险圆弧滑动面上安全系数 随时间变化曲线Fig.7 Time evolution of the safety factors on the slip surface

图8为水位上升过程中的各个力对安全系数的贡献及作用在危险圆弧滑动面上的各个力随时间的变化情况。由于土体自身重力在水位上升初期为1 640.40 kN/m，到模拟结束后为1 665.00 kN/m，上升了24.60 kN/m，较其自身而言，增加的重力可以忽略不计，因此土体自身重力对安全系数的贡献几乎不变，将不做详细讨论。

图8-a为水位上升作用下毛细压力对安全系数的贡献随时间变化的曲线及作用在危险圆弧滑动面上的总毛细压力随时间的变化规律。观察图8-a，作用在危险圆弧滑动面上总的毛细压力值小于0，随着水位上升毛细压力值逐渐下降，当水位稳定后，毛细压力值缓慢上升，因此毛细压力对安全系数的贡献先逐渐增加后逐渐减小，所以，安全系数图像中不考虑毛细压力贡献的安全系数与考虑孔隙气压力和毛细压力贡献的安全系数的差值先逐渐增加后逐渐减小。

图8-b为水位上升作用下孔隙气压力对安全系数的贡献随时间变化的曲线及作用在危险圆弧滑动面上总的孔隙气压力随时间的变化规律。如图8-b所示，水位上升过程中作用在危险圆弧滑动面上总的孔隙气压力先快速增加，即从0 kN/m增加至805 kN/m，由于孔隙气压力不利于边坡稳定，总的孔隙气压力增加导致孔隙气压力对安全系数的不利贡献从0上升至-0.16，当水位稳定后，总的孔隙气压力继续由805 kN/m逐渐增加至1 068 kN/m，因此安全系数图像中，不考虑孔隙气压力贡献的安全系数与考虑孔隙气压力和毛细压力贡献的安全系数的差值先快速增加后缓慢增加。

图8-c为水位上升作用下坡外水压力对安全系数的贡献随时间变化的曲线及作用在危险圆弧滑动面上总的坡外水压力随时间的变化规律。水位上升过程中坡外水压力从0 kN/m快速上升至960 kN/m，使得坡外水压力对安全系数的贡献从0快速增加至1.56，当水位稳定以后，由于坡外水位保持不变，因此作用在危险圆弧滑动面上总的坡外水压力值始终为960 kN/m，坡外水压力对安全系数的贡献也始终稳定在1.56。

8-a 毛细压力8-b 孔隙气压力8-c 坡外水压力图8 滑动面上各个力对安全系数的贡献及各个力随时间的变化过程Fig.8 Contributions to the safety factors from different kinds of pressures and different kinds of pressures acting on the slip surface

4 结论

通过建立水-气二相流模型，研究了水位上升作用下土坡渗流场变化规律，并进一步分析了土坡稳定性变化规律，得到以下结论：

(1)将水位上升引起的饱和-非饱和渗流问题作为水-气二相流问题研究更符合实际物理意义。从物理现象上来看，水位上升是土壤中“水驱替气”的过程，因此，随着水位上升，孔隙气压力、孔隙水压力、毛细压力和水饱和度均逐渐增加。

(2)岸坡的安全系数随着水位的上升逐渐增加，水位稳定后，安全系数逐渐下降，且考虑孔隙气压力和毛细压力贡献的安全系数大于不考虑毛细压力贡献的安全系数，但小于不考虑孔隙气压力贡献的安全系数。因此，毛细压力的存在有利于土坡稳定，而孔隙气压力的存在不利于土坡稳定。

(3)水位上升作用下，毛细压力逐渐减小，水位稳定后，毛细压力逐渐增加，因此毛细压力对安全系数的贡献先逐渐增加后逐渐减小；水位上升过程中，孔隙气压力先快速增加，当水位稳定后，孔隙气压力缓慢增加，因此孔隙气压力对安全系数的不利贡献先快速增加后缓慢增加；坡外水压力逐渐增加后稳定不变，因此，坡外水压力对安全系数的贡献先逐渐增加后保持不变。