鄱阳湖区洪水全过程下堤防稳定性有限元模拟计算分析

熊 威，王 姣，王萱子，胡 强，彭圣军

(江西省水利科学院，江西 南昌 330029)

0 引 言

堤防是鄱阳湖滨湖地区主要的防洪工程，与枯水期相比，鄱阳湖在汛期的水面面积与水位差别较大[1]，且汛期降雨频繁，易造成堤防渗透破坏及边坡失稳，因此，开展鄱阳湖区堤防边坡稳定性研究具有重要意义。方景成等[2]开展了库水位下降速率和降雨强度与边坡安全系数的敏感性分析。刘俊新等[3]对边坡稳定性影响因素进行了数值计算，认为降雨强度和历时及长期强度的影响最明显。王乐等[4]基于Geo-slope 2012软件，对不同类型降雨与库水位联合作用下的土质边坡滑坡失稳过程进行了分析。倪沙沙[5]考虑了基质吸力对坝体材料抗剪强度的影响，将暂态渗流场与极限平衡法相结合，计算了坝坡的稳定性。岑威钧等[6]考虑气相与基质吸力的影响，开展了堤防渗流及稳定性分析。吕雪超[7]采用数值模拟与模型实验相结合的方法，分析了堤防瞬态渗流特性及堤坡稳定性。已有研究的对象大部分为大坝及土质边坡，部分对堤防的研究也仅单独考虑水位升降或降雨的影响。本文基于2020年鄱阳湖超历史大洪水背景，模拟了堤防在洪水全过程作用下的孔隙水压力(以下简称“孔压”)与安全系数变化，并考虑降雨强度与降雨历时等不利因素影响，分析了鄱阳湖堤防边坡的稳定性变化规律，可为堤防防汛抢险提供参考。

1 工程概况

该堤防工程位于鄱阳湖流域饶河支流昌江下游，堤线总长9.607 km，圩堤内集水面积56.0 km2，保护面积15.3 km2，保护耕地1万hm2(1.5万亩)，保护区内人口1.8万。堤防建于20世纪60年代，后经多次除险加固达到现有规模，堤身为典型的二元结构；上层为透水性较低的黏土，下层为透水性较强的砂壤土。堤顶高程24.50 m，宽8.0 m，临水坡坡比1∶3，背水坡坡比1∶3.5，警戒水位19.50 m。

2 研究方法

2.1 饱和-非饱和渗流方程

饱和-非饱和渗流微分方程的张量表示[4]如下：

(1)

式中：xi，xj为空间坐标张量；kij为渗透张量；kr为透水率；Hc为水头；q为源汇项；C(Hc)为容水度；θ为水头函数；n为孔隙率；Ss为单位贮水量；t为时间。

Fredlund-Xing[8-9]模型是进行土体非饱和渗流计算的主要方法之一，其控制方程为

(2)

式中：θw，θs分别为体积含水量、饱和体积含水量；Cφ为修正项；φ为负孔压；a，n和m为经验拟合参数。

由式(2)可知，通过土体的饱和体积含水量θs可估算出渗透系数函数：

(3)

式中：kw与ks分别为渗透系数与饱和渗透系数；j与N为最小负孔压与最大负孔压；i为单位间距；θs为饱和体积含水量；y为虚变量；ψ为第j步的负孔压；θ0为初始体积含水量。

2.2 非饱和抗剪强度

土体孔隙水压力的增大会降低土的基质吸力，导致其抗剪强度削弱。非饱和土抗剪强度τ的计算采用Fredlund[10]提出的公式：

τ=c′+(σn-ua)tanφ′+(ua-uw)tanφb

(4)

式中：c′与φ′分别为有效黏聚力与有效内摩擦角；σn为法向应力；ua与uw分别为孔隙气压力与孔隙水压力；φb为受基质吸力变化影响后的内摩擦角。

2.3 饱和-非饱和边坡抗滑安全系数

在非饱和土的情况下采用简化毕肖普法计算边坡稳定安全系数Fs的表达式[11]如下：

Fs=

(5)

3 计算模型及工况

3.1 有限元模型及边界条件

堤防典型横断面见图1。堤身材料参数见表1。采用Geo-Studio软件建立有限元模型，如图2所示，网格单元类型为三角形单元和四边形单元，全局单元尺寸为1 m，模型共有493个节点，432个单元。为监测堤身不同部位孔压的变化，在堤防上下游设置监测点A，B，C(图1)。将堤防的迎水坡边界(DEI)设置为水头，水位随着时间上升和下降，堤顶(DG)和背水坡(GFH)设置为降雨边界。

图1 堤防典型横断面

表1 堤身材料参数

图2 有限元模型

3.2 土水特征曲线

土水特征曲线通过材料的饱和体积含水率和SEEP/W模块提供的黏土样条函数和砂壤土样条函数估算，渗透系数曲线采用Fredlund-Xing模型拟合，参数设置为a=12.0 kPa、n=1.0、m=1.5。堤身材料土水特征曲线拟合结果见图3，渗透系数曲线拟合结果见图4。

3.3 计算工况

为尽可能接近2020鄱阳湖超历史大洪水工况，根据实测洪水资料，将洪水过程划分为以下3个阶段：① 不同水位上升速率联合不同降雨强度作用阶段；② 持续高水位情况下降雨强度及降雨历时作用阶段；③ 不同水位下降速率与不同降雨强度作用阶段。各阶段对应工况的详细划分见表2。

图3 土水特征曲线

图4 渗透系数曲线

表2 计算工况

4 结果与分析

4.1 堤身孔压变化

图5为阶段1堤身孔压整体变化。随降雨过程的持续与水位的上升，堤身孔压变化明显。图6～8为本次洪水全过程3个监测点的孔压变化过程：孔压变化规律总体为随上游水位上升而逐渐升高，随上游水位下降而逐渐降低；不同位置孔压变化规律基本相同，上游坡孔压响应速度明显高于下游坡。未考虑降雨时，水位上升(下降)速率越大，孔压变化速率越大，但最终孔压值相同，水位上升(下降)速率不影响堤身最终孔压值。考虑降雨时的最终孔压值比未考虑降雨时大，且降雨强度越大，孔压变化速度越快，降雨历时越长、高孔压持续时间越长；降雨强度与降雨历时对监测点孔压最终状态无影响。其中，考虑降雨时，阶段3监测点的孔压先增大后减小，主要是因为降雨入渗使土体孔压增大，并随降雨的停止而逐渐减小。在本次洪水全过程中，阶段1、阶段2孔压变化幅度大于阶段3，此时若遭遇降雨，则会明显提高堤身孔压，降低土体基质吸力，进而易引起堤防出现管涌、渗漏及滑坡等险情。

图5 阶段1孔压整体变化

4.2 上下游坡安全系数变化

图6 阶段1各监测点孔压变化

图7 阶段2各监测点孔压变化

图8 阶段3各监测点孔压变化

图9～11为上下游坡安全系数的变化过程。经分析可知：在上游水位上升过程中，上游坡安全系数先增大后减小，下游坡安全系数逐渐降低；这是因为随着上游水位升高，上游坡受水压加载作用，安全系数会先增大，在上游水位上升以后，土体的强度参数及有效应力降低，从而安全系数逐渐减小。持续高水位阶段，上下游坡安全系数在降雨的作用下逐渐减小，降雨停止时会小幅上升。在考虑降雨作用的上游水位下降过程中，上游坡安全系数逐渐减小，下游坡安全系数先增大后逐渐减小并趋于稳定，主要原因是：降雨初期，堤身浸水部分土体抗滑力下降程度小于下滑力的下降程度，使下游坡稳定性在降雨初期有一定升高；随着降雨停止和雨水渗入土体内部，土体强度参数和有效应力降低，从而稳定系数减小。水位变化速率对上下游坡最终安全系数无影响；考虑降雨时的上下游坡最终安全系数比未考虑降雨时小，但降雨强度和降雨历时对上游坡最终安全系数无影响。综上所述，降雨对堤坡不同阶段安全系数有不同程度的影响，即使在水位下降阶段，上游坡遭遇降雨时安全系数仍会明显下降。

图9 阶段1安全系数变化

图10 阶段2安全系数变化

图11 阶段3安全系数变化

5 结 论

本文通过有限元模拟，分析了鄱阳湖堤防在2020年洪水全过程作用下的孔压与安全系数变化，并考虑降雨强度与降雨历时等不利因素影响，研究了堤防的边坡稳定性，得到以下结论。

(1) 水位上升(下降)速率不影响堤身最终孔压；降雨会明显提高堤身孔压、降低土体基质吸力，使堤防易出现管涌、渗漏及滑坡等险情。

(2) 水位上升(下降)速率不影响堤身最终安全系数；降雨会明显降低上下游坡安全系数。需特别注意在水位下降过程中，上游坡在遭遇降雨的情况下易发生滑坡险情，因此在洪水全过程均应加强堤防的安全巡查。