考虑渗流-应力耦合的不同降雨条件下非饱和土边坡稳定性分析

周亦良,贺 其,李 明,孙志远,何丰瀚

(南华大学 土木工程学院,湖南 衡阳 421001)

0 引 言

我国地形地貌多样,地质构造复杂,伴随着地震和极端天气的影响,地质灾害频发,给人民生命财产和生活环境造成了严重的威胁和破坏。根据我国自然资源部公布的全国地质灾害灾情,2021年全国共计发生地质灾害4 772起,造成80人死亡、11人失踪,直接经济损失32亿元。其中斜坡地质灾害(崩塌、滑坡、泥石流)占93.4%。根据多年的地质灾害发生规律,降雨是斜坡地质灾害中的重要致灾因子,尤其是强度大、历时长的暴雨或特大暴雨更是导致边坡失稳的主要诱导因素[1-2]。

自然环境中边坡土体大多数处于非饱和状态,降雨条件下边坡稳定性问题本质上是雨水入渗导致孔隙水压力增大,非饱和土基质吸力降低,土体抗滑力下降,诱发边坡失稳破坏[3]。李宁等[4]考虑坡面倾斜的影响及非饱和土特性,进行了低强度长时间降雨和高强度短时间降雨两种情况下浅层边坡稳定性计算,得到稳定性系数与降雨强度、历时的变化关系;O.DAVIES等[5]基于现场监测数据建立水力耦合模型,验证分析边坡稳定性,确定降雨、孔隙水压力变化与边坡失稳的相关性;L.J.WANG等[6]模拟了长时间强降雨条件下边坡渗流和变形的耦合过程,表明边坡失稳破坏是由于基质吸力减小和地下水位显著升高引起的;汪丁建等[7]基于Green-Ampt模型,分析了湿润锋深度随降雨时间的变化过程,在此基础上利用整体极限平衡法分析了滑坡稳定性;王一兆等[8]通过分析降雨期和停雨期的渗透系数和孔隙水压力,得出了边坡浅层稳定性的变化规律;吕雨桦等[9]探索了大暴雨(150 mm/d)工况下非饱和土边坡渗流场、应力场、变形场的变化过程及边坡稳定性的演化规律;L.WANG等[10]强调了土水特征曲线和渗透系数对非饱和土渗流的影响,提出一种贝叶斯方法表征土水特征曲线和渗透系数,能更准确预测边坡稳定性;Z.SU等[11]考虑渗流场和应力场的耦合效应,研究了不同降雨强度、类型对下游坝坡稳定性的影响,认为无论水位升降,降雨都会降低下游坝坡的安全系数。

已有研究表明,降雨条件下非饱和土边坡稳定性主要受到边坡雨水入渗的影响。在降雨入渗过程中,坡体内的渗流场会随之发生改变,渗流会产生渗透力,对边坡原有的应力场产生影响,反过来应力场的改变又会引起孔隙率等变化,对渗流场产生影响,渗流场与应力场是互相影响且可逆的[12-13]。因此,渗流-应力耦合作用是降雨边坡必须要考虑的因素,但目前却较少涉及考虑渗流应力耦合作用的同时进行不同降雨工况下边坡稳定性的深入研究。笔者基于非饱和土流固耦合理论,通过Geo-Studio有限元数值分析软件建立非饱和土边坡模型,对非饱和土边坡渗流场、应力场进行直接耦合分析,并在此基础上探讨了不同降雨强度、不同降雨类型、不同降雨持时对边坡稳定性的影响,以期为降雨条件下非饱和土边坡稳定性研究提供参考。

1 渗流-应力耦合原理

基于地下水连续运动方程,根据达西定律,得到二维渗流公式:

(1)

式中:kx、ky分别为x、y方向的渗透系数;uw为孔隙水压力;θw为体积含水量;γw为水的重度;t为时间。

非饱和土中体积含水量表示土体中储水体积与总体积的比值,体积含水量的改变依赖于应力状态的改变和土体的性质。在渗流-应力耦合计算中,假定孔隙气压保持不变,土体材料特性不变,那体积含水量的改变可用体积应变和孔隙水压力变量来表示:

Δθw=βΔεv-ωΔuw

(2)

(3)

(4)

式中:εv为体积应变;v为泊松比;KB为体积模量;H表示为非饱和土与基质吸力相关的模量;E为土结构的弹性模量;R为与体积含水量随基质吸力变化而变化的模量。

(5)

将式(2)代入(5)中,分部积分得到:

(6)

使用有限元近似法,方程表示为:

(7)

式中:B为梯度矩阵;Kw为渗透系数矩阵;N为形函数行矢量;m为各向同性单元张量;δ为节点位移;Vn为边界流量。

在Δt时间内渗流方程利用差分方法表示为:

(8)

式中:Lf为渗流耦合矩阵;Kf单元刚度矩阵;MN为质量矩阵;Q为边界节点的渗流。

利用位移增量和孔隙水压力增量作为场变量,土体应力平衡方程可表示为:

K[Δδ]+Ld[Δuw]=F

(9)

式中:K为刚度矩阵;Ld为耦合矩阵;F为外部荷载。

联合渗流方程(8)和平衡方程(9)可得有限元分析的耦合方程,通过相应的边界条件即可对该耦合方程进行求解。

2 案例分析

边坡稳定性分析模型选自澳大利亚计算机应用协会ACADS考核题,为了减小边界约束条件对边坡应力-变形场的影响,将边坡分析区域范围进行了扩大,几何断面如图 1。土体天然重度取20 kN/m3,内摩擦角16.5°,黏聚力5 kPa,饱和含水量25%,泊松比0.3,饱和渗透系数2×10-6m/s。模型边界条件两侧采用水平约束,底部边界水平、竖向约束,初始地下水位如图1。根据我国气象主管部门降水强度等级划分,采用的降雨强度分别为小雨(5 mm/d)、中雨(20 mm/d)、大雨(40 mm/d)、暴雨(80 mm/d)。模型计算中降雨强度均会小于土体入渗能力,斜坡水力边界直接设置为单位流量,入渗大小等于降雨强度。同时选取了坡脚A点(30,20)、坡面中间B点(40,25)、坡顶C点(50,30)以及边坡表层内部D点(40,24)、E点(40,23),进行降雨过程中孔隙水压力和应力变化的分析。

图1 边坡计算模型

以均质粉土边坡为例,利用软件自带的样本函数,设定饱和含水量为25%,拟合得到土体的土水特征曲线如图2,然后基于体积含水量函数以及土体饱和渗透系数2×10-6m/s,利用Van Genuchten 模型估算渗透系数曲线如图3。土体材料采用线弹性模型,其中弹性模量函数如图4。

图2 土-水特征曲线

图3 渗透系数曲线

图4 弹性模量函数曲线

非耦合作用实现过程是首先在 Seep/W中先做瞬态渗流分析,得到不同步长下渗流场的水压分布;然后将渗流场中的孔隙水压力值输入到Slope/W模块中完成边坡稳定性分析。耦合作用实现过程是首先在Sigma/W原位分析模块中通过给定的地下水位位置以及最大允许孔隙水压力的初始条件,计算初始渗流场和应力场;然后考虑不同降雨条件(强度、持时、类型),在Sigma/W 应力-孔隙水压力耦合分析模块中对整个降雨过程进行瞬态渗流及应力分析,得到每个时间步长(2 h)的孔隙水压力和应力分布;最后通过将节点水头、应力数据导入Slope/W模块进行边坡稳定性分析,采用有限元应力单元法得到每个时间步长的边坡稳定性安全系数。

2.1 渗流和应力对比分析

图5给出了耦合作用下暴雨(80 mm/d)工况下边坡孔隙水压力0、12、24 h时刻的分布图。降雨开始后,坡面位置处孔隙水压力迅速增大,等值线分布密集,从坡面向坡内呈现出高-低-高的特点。表明雨水入渗,边坡表层土体最先做出响应,土体孔隙水压力急剧上升。随着降雨的持续,降雨入渗深度增加,相同位置的孔隙水压力随时间增加而增大。

图5 不同时刻边坡孔隙水压力分布

图6(a)对比了暴雨(80 mm/d)工况下考虑/不考虑渗流-应力耦合作用的不同位置A,B,C,D,E点的孔隙水压力时程变化特征。坡面位置A,B,C点的孔隙水压力变化趋势基本相同,随着降雨持续而快速增长,然后逐渐趋于平稳,呈现暂态饱和。坡内位置D点孔隙水压力降雨在开始后0～4 h基本保持不变,4 h后孔隙水压力随着降雨持续逐渐增长,到降雨结束仍保持上升趋势。坡内位置E点孔隙水压力从开始到降雨结束基本保持平稳,降雨入渗尚未影响到该深度。表明随着降雨持续,坡面最先做出响应后,雨水入渗水分向坡内迁移,入渗影响深度逐渐增加,引起体积含水量和孔隙水压力的增加。

图6 不同特征点孔隙水压力、有效应力和降雨时间的变化关系

与非耦合分析相比,耦合分析的孔隙水压力增长速度要大于非耦合分析,相同时间同一位置耦合计算的孔隙水压力比非耦合的稍大,这种差异在开始快速阶段最明显。表明应力场作用加快了土体渗流速度,使得土体孔隙水压力增长速度以及基质吸力消散速度比非耦合计算的更快。

图6(b)中分别给出了耦合作用下A,B,C,D,E点在Y方向有效应力时程变化特征。对比图6(a)可以看出,相应位置的有效应力与孔隙水压力呈负相关关系。随着降雨的进行,边坡表层非饱和土体的体积含水量增加,基质吸力降低,孔隙水压力逐渐增大,破坏了其原始应力平衡状态,对应有效应力(u-uw)下降,对边坡抗剪强度产生不利影响。

2.2 不同降雨条件下边坡稳定性分析

为研究不同降雨强度对边坡稳定性的影响,采用单一变量法,降雨类型和持时相同,降雨强度不同,取降雨持时为24 h,降雨强度取5、20、40、80 mm/d。图7(a)～图7(c)给出了不同降雨强度下特征点B,D,E的孔隙水压力与降雨持时的关系曲线。当降雨强度为80 mm/d时,坡面位置B点孔隙水压力增长最迅速,并随降雨的持续而趋于暂态饱和。随着降雨强度下降(40～5 mm/d),降雨结束时刻孔隙水压力大小分别为-22.9,-35.9,-62.2 kPa,仍保持着增长趋势。不同降雨强度下D点孔隙水压力变化规律更加明显,在降雨强度为80 mm/d时,孔隙水压力还保持着增长趋势,但降雨为5 mm/d,可以看到孔隙水压力基本保持平稳。坡内E点孔隙水压力基本不变,与边坡表层土体反应迅速形成鲜明对比。说明在高强度降雨作用下表层孔隙水压力受降雨入渗影响明显而迅速;随降雨强度下降,雨水入渗能力和影响深度减小,土体孔隙水压力增长速度和大小也会显著减小。

图7 不同降雨强度下孔隙水压力、稳定性安全系数和降雨时间关系

在耦合与非耦合对比分析中,不同降雨强度条件下边坡的稳定性安全系数都会随降雨的持续而不断减小,随着降雨强度的增加而减小,耦合分析得到的安全系数比非耦合分析的安全系数偏低〔图7(d)〕,与文献[13]结论一致。正如耦合分析雨水下渗引起孔隙水压力上升速度更快,相同时间内基质吸力消散更快,使得边坡稳定性更低,同时雨水入渗形成的渗透力作用于土体骨架,使边坡产生更大的变形,不利于边坡稳定。

对比不同降雨强度可以发现,从5～80 mm/d,稳定性安全系数会随之降低,而且下降速率会逐渐增大。因为降雨时,雨水从边坡表层开始往下渗,随着降雨强度的增加,雨水入渗的速度和深度都逐渐增大,坡内的基质吸力消散区域更大,土体有效应力大范围降低,土体抗剪强度减小,边坡稳定性安全系数越低。

为了研究降雨持时影响,控制总降雨量相同(80 mm),降雨类型统一采用均匀型,即:降雨强度为80 mm/d时,降雨持时取24 h(1 d),降雨强度为40 mm/d时,降雨持时取48 h(2 d),以此类推。如图8所示,可以发现当降雨强度为80 mm/d时,降雨结束时刻(1 d)边坡安全系数为1.626;当降雨强度为5 mm/d时,降雨结束时刻(16 d)边坡安全系数为1.612,随着降雨持时越长,降雨结束时刻边坡稳定性安全系数越低。总降雨量相同的条件下,随着降雨持时增长,降雨入渗影响深度的增加,土体体积含水率增加,边坡土体的强度下降,进而导致土体抗滑力下降而滑动力增大,边坡稳定性安全系数下降。对比短时间强降雨(80 mm/d),当长期小雨作用时(5 mm/d),雨水入渗能力虽小,整个降雨过程中边坡表层也并未出现暂态饱和区,安全系数下降速度也相对较慢;但长时间雨水入渗影响范围会逐渐增大,引起地下水位高度逐渐升高,导致降雨结束时刻边坡稳定性安全系数反而更小一些。表明边坡在当地出现长期小雨情况下,发生失稳的可能性比短时间强降雨可能更大,仍然不容忽视。

图8 总降雨量相同条件下稳定性安全系数和降雨时间的关系

为了研究降雨类型影响,降雨量和降雨持时相同,降雨类型取后锋型、均匀型和前锋型[15],结果如图9。

图9 不同降雨类型下稳定性安全系数和降雨时间的关系

由图9可知:随着降雨的持续,安全系数值不断减小,但降雨类型不同,边坡稳定性安全系数下降速度不同。这是由于不同降雨类型其降雨强度随时间的变化不同,如均匀型降雨强度随时间不变,因为其稳定性安全系数基本呈线性下降;前锋型降雨强度随着降雨时间的增加而减小,因此其稳定性安全系数下降速度先快后慢,呈凹形抛物线;后锋型降雨强度随着降雨时间的增加而增大,因此其稳定性安全系数下降速度先慢后快,呈凸形抛物线。降雨结束后,边坡的稳定安全系数仍保持继续下降,在60 h左右,稳定性安全系数趋于一致。由于降雨结束后,暂态饱和区的水分会继续向边坡下部渗透,在降雨结束一段时间内对边坡稳定性造成持续性危害。

3 结 论

基于非饱和土流固耦合理论,通过对非饱和土边坡有效应力和孔隙水压力对比分析,探讨了不同降雨强度、不同降雨持时、不同降雨类型对边坡稳定性的影响,得出以下结论:

1)渗流-应力耦合作用使坡内雨水下渗速度变快,受影响土体的孔隙水压力增长速度及基质吸力消散速度比非耦合分析的更快,导致耦合分析得到的安全系数更小。因此,在进行边坡稳定分析时,应考虑到渗流-应力耦合作用的影响。

2)强降雨入渗导致非饱和土边坡表层的孔隙水压力迅速增大,出现暂态饱和区,体积含水量的增加造成基质吸力、有效应力下降,弱化土体抗剪强度。在进行边坡治理时,要特别注重坡面的防排水。

3)降雨强度越大、降雨持时越长,边坡稳定安全系数越小。降雨强度影响安全系数的下降速度,降雨时长影响安全系数的最小值。降雨类型不同稳定性安全系数下降速率不同,当边坡稳定性接近临界状态时,前锋型降雨是触发边坡失稳发生时刻最早的降雨类型。因此除了考虑临界降雨量,长时间的连绵小雨以及降雨类型都是边坡稳定性不容忽视的影响因素。