边坡降雨入渗的饱和-非饱和有限元数值分析

2010-07-09方燕

湖南工程学院学报（自然科学版） 2010年3期

方燕

(金华广播电视大学理工学院,金华321000)

0 引言

降雨入渗时边坡的稳定问题是一个牵涉到饱和-非饱和状态水的渗流和含水量变化的复杂问题.对于这种非典型饱和流固耦合现象,雨水入渗的瞬态渗流场、土坡变形以及各种边界条件的影响是很重要的影响因素.前人做了很多研究工作,Fredlund[1]利用非饱和土壤水运动理论和土固结理论,推算出饱和-非饱和渗流控制方程:文献[2]对高程不同的各类边坡,对其基质吸力进行了分析;文献[3]结合强度折减法计算边坡稳定性中的渗流情况;文献[4]在程序子上进行二次开发,编写了SLOPE/W和SEEP/W的二次程序,用以计算渗流问题.本文利用运动饱和-非饱和渗流-应力场耦合的方法,着重模拟降雨条件下土质边坡不同阶段的渗流场分布情况,对工程算例中的位移,孔隙压力等参数进行分析,为研究降雨引起的边坡失稳提供了有价值的参考依据.

1 降雨条件下边坡的渗流控制方程

1.1 饱和-非饱和渗流数学模型

结合Fredlund[1]的理论,多孔介质中各向异性饱和-非饱和渗流控制方程[5]为:其中,kij——介质饱和渗透张量;kr(hp)——介质相对渗透率,在非饱和区的值介于0和1之间,饱和区的值为 1;h——总水头;Q——源汇项;C(hp)——容水度,在非饱和区的为体积含水率),饱和区中的值是0;β——选项系数,在非饱和区的值是0,在饱和区的值是1;Ss——饱和岩土体单位贮水率,在非饱和区的值是0,在饱和区为某一常数;xi——坐标;t——时间变量.方程的定解条件为:

初始条件:

1.2 降雨入渗边界的处理方法

降雨入渗时雨水的入渗量与很多因素有关,其中包括降雨持续时间、降雨强度、土壤的渗透性、吸力状况等[6].通过达西定律可知地表上不同方向的

当降雨完全入渗时,不会形成地表径流.这时,边界条件属于第二类边界;

式中:Γ5——降雨入渗面边界;qr——降雨强度.

如果只有部分降雨渗入到土体中,证明土壤的入渗能力小于降雨强度,此时会形成径流,当作是第一类边界计算:

1.3 渗流一应力场耦合

为了使力学平衡理论与渗流理论结合,可利用孔隙水压力和节点位移的自由度作空间离散,得出如下矩阵[7]:

平衡方程为:入渗透能力为:

利用一个时间步长内形成的迭代求解基础,放入到具体的求解过程中,结合已知的各种条件,当单元产生塑性变形的时候,用弹塑性矩阵将弹性矩阵用取代,即能进行迭代计算.

2 数值算例分析

2.1 有限元计算模型

本文以一边坡作为研究对象,尺寸如下:边坡的底部长为160 m,左侧为总坡高,高度为40 m,右侧坡高为20 m,其中坡度比为1∶2,左侧边界距离坡顶为80 m,右侧边界距离破趾为40 m,边坡土体的渗透系数为2.4 m/day,泊松比v=0.3,变形模量E=70 MPa,饱和容重 19.7 k N/m3,Gamma容重15.9 kN/m3,粘聚力c=12 k N/m2,摩擦角φ=200.计算过程采用莫尔-库伦本构模型,兼顾考虑非饱和特性,所建有限元模型如图1所示,共生产5352个单元,5545个节点.

图1 有限元模型网格图

2.2 降雨入渗时的边界条件

在有限元计算过程中,设置边坡底部为零流量边界,左侧边界水头为32 m,右侧边界水头为16 m,边坡的上部由于考虑有降雨的情况出现,降雨强度为 15 mm/h,历时60 h,故设置为定流量边界,因此在整个降雨入渗分析过程中,边坡顶及坡面都会受到降雨作用.由于有定水头边界的存在,首先要进行非降雨情况的稳定渗流计算,再将降雨的因素加入到模型中,并设置时步,从0～60 h,统计出孔隙水压力、边坡横向位移、竖向位移及整体稳定性的情况.

2.3 结果分析

图2 降雨后不同时间内孔隙压力图(单位:kPa)

在计算过程中,本文对某些感兴趣的单元点进行了实时监测.如图2所示,通过所布置的监测点数据以及孔隙压力等值线图,显示出降雨入渗后会引起非饱和土体中孔隙水压力上升,随着土体中的含水量增加,浸润面也随之上升.以坡顶的监测点为例,孔隙压力由-66 kPa增至-62 k Pa再至-47 k Pa,规律的趋势显示非饱和去基质吸力会随着降雨不断入渗而下降,最终会导致边坡土体的总粘聚力也相应减少.

图4 降雨后不同时间内竖向位移图(单位:cm)

如图3和图4两个位移等值线图所示,水平位移呈现出中部增加明显,坡趾处等值线最为密集,位移梯度较大,由12～60 h,位移由5.75 mm增加到15.3 mm.竖向位移的分布遍及很大的区域,降雨60 h之后,坡顶竖向位移沉降了大约9.66 mm,其中坡趾位移最小.通过对比水平与竖向位移可见一个规律:相对于深层区域,边坡的浅层区域更容易发生破坏.