基于三维有限元的凤亭河水库坝坡稳定性分析*

(1.广西水利科学研究院,广西 530023;2.中国水利水电科学研究院, 北京 100038;3.流域水循环模拟与调控国家重点实验室，北京 100038;4.水利部水工程建设与安全重点实验室，北京 100038;5.北京鼎创力维科技有限公司，北京 100085)

工程上常用的土石坝坝坡稳定分析方法有：传统刚体极限平衡法、基于有限元分析的刚体极限平衡法、塑性极限分析和模糊极值理论。刚体极限平衡法发展较为成熟，因此目前在坝坡稳定分析的工程实践中基本上都采用这种方法。《碾压式土石坝设计规范》(DL/T 5395—2007)中明确指出：坝坡抗滑稳定计算应采用刚体极限平衡法，但按照刚体极限平衡的思想进行力学分析，没有考虑土体本身的应力—应变关系。塑性极限分析和模糊极值理论考虑了土体本身的应力—应变关系，可以假定土体为刚塑性体，但计算较为复杂。基于有限元分析的刚体极限平衡法结合圆弧滑动法进行分析，假定滑动面为圆弧面，仍采用条分法进行计算，所不同的是滑裂面处的内力系由该处的应力转化而成，这和土体的实际工作状态相符。最危险滑弧仍采用试算法，但如果试算的滑弧数目过多，则会出现工作量大、效率低的现象。为解决这个问题，在假定滑面时，首先根据数值求解的应力分布成果分析坝体可能的滑裂面位置，限定最危险滑裂面搜索的范围，可减小工作量。计算得到的安全系数与传统极限平衡法相比略小[1]，这是因为在计算内力时充分考虑了土体内部变形对滑裂面的影响。该种分析方法可以为类似工程的建设提供理论参考。

1 基本理论与计算方法

1.1 邓肯-张E-B模型

增量型的应力—应变关系式一般为：

{Δσ}=[D]{Δε}

(1)

邓肯-张采用下面的双曲线方程表示由三轴试验得到的土体应力—应变曲线：

(2)

(3)

式中 (σ1-σ3)f——土体破坏时的主应力差,Pa，

(σ1-σ3)ult——最大主应力差,Pa，取双曲线的渐近线所对应的主应力差。

Rf——破坏比;

S——应力水平。

S的物理意义为土体中实际的主应力差(σ1-σ3)与极限状态下的主应力差(σ1-σ3)f之比，S>1表示土体已被剪坏，S=1表示土体处于极限状态，S<1表示土体未被剪坏。因此S是坡体内是否产生极限平衡区的一项重要指标。

1.2 基于有限元法边坡稳定性分析

通过有限元分析可以获得任意滑动面上单元高斯点上正应力σn、剪应力τn。沿滑动面上某一微段(单元)积分(见图1)，可以得到该微段(单元)的阻滑力Fz与滑动力Fh，具体计算公式如下：

(4)

(5)

式中φ——任意滑动面上的内摩擦角,°;

c——任意滑动面上的黏聚力,kN/m2。

图1 滑面示意图

1.3 刚体极限平衡理论

由于下文实例中边坡滑坡体表层为松散介质材料，刚体极限平衡稳定性分析时原则上采用适合于土体介质的稳定性计算方法[2]。对堆积体边坡主要考虑采用Bishop法、Jan-bu法和Morgenstern-Price法，本文选用Bishop法，安全系数公式为：

(6)

其中：

(7)

式中K——安全系数，由迭代得到；

c——土体黏聚力,kN/m2；

l——滑块的底边长度，m；

W——土条的重量,kg；

α——土条底面的倾斜角,°；

φ——土体的内摩擦角,°。

2 工程实例概况

凤亭河水库位于广西壮族自治区南宁市良庆区大塘镇，工程坐落于珠江流域西江水系郁江支流八尺江上游，总库容5.07亿m3，属多年调节水库。水库正常高水位175.00m，设计洪水位177.61m(P=1%)，校核洪水位178.88m(P=0.05%)，凤亭河水库主坝一座，为碾压式均质土坝，坝顶高程180.30m，最大坝高53.62m，坝顶长192.3m，坝顶宽8.0m。坝体结构布置如图2所示。

凤亭河水库坝区地处低山丘陵的地貌分区，一般河谷呈U形，峡谷呈V形，河床高程123～126m，山顶呈浑圆状居多，脊状罕见，山坡平缓，冲沟发育。

坝区出露为白墨系那贞统地层，大致可以分为两部分：上部为紫红色泥岩、泥质页岩、砂岩、砂泥岩互层，厚约12.4m；下部为砾状砂岩、细砂岩、粗砂岩、夹泥岩，砾岩，厚度约26m，岩性一般为轻度变质。

3 计算模型与条件

3.1 计算模型

为消除边界影响，模型范围根据给出的资料，上下游各延伸出坝体边界100m，左右岸延伸100m，基础向深部延伸115m。模型坐标采用笛卡尔坐标系，坐标原点选取在坝顶轴线中心点处，X方向沿坝轴线方向，Y方向与水流方向一致，Z方向垂直向上。模型分别模拟了坝体主体、防渗墙、排水体、坝基不同地层等结构构造。采用八节点六面体等参单元，模型单元数113852个，节点122748个。计算网格模型如图3所示。边界条件为底部固定、周边为链杆约束。

3.2 计算参数与计算工况

a.计算参数。计算分析采用的材料参数取自广西风亭河水库管理局提供的设计报告，按照工程经验类比的方式加以确定，主要参数取值如下表1。

图3 凤亭河水库大坝稳定分析模型网格图

b.计算工况。本次计算主要考虑正常蓄水位下的坝体应力、应变及稳定性。将渗透水压力作为外荷载进行考虑。即：按照大坝正常蓄水位175.00m分析计算渗流场，将渗透力作为体力施加在节点上进行静力计算分析。

4 计算结果分析

4.1 大坝变形分析

大坝最大变形发生在河床坝段。施工完建期，坝体上游坡面向上游最大水平位移3.35cm，下游坡面向下游最大水平位移3.90cm，且均出现在坝体内部。渗流对坝体水平位移有一定影响，水平位移量从上游到下游逐渐增大，坝高1/2以上部位的水平位移增加明显，坝顶部的水平位移增加约0.274cm。上游坝面约1/3坝高、坝轴线上游侧仍然有向上游水平位移，分别比竣工期减少0.278cm和0.244cm，且均出现在坝体内部，范围与量级比施工完建期小，位移分布云图如图4所示。坝体最大垂直位移10.29cm，比竣工期减小约0.024%，减少幅度合理，产生最大沉降的平面位置偏向下游，竖向位置变化不大，仍然在坝高1/2～2/3处。

图4 河床坝段位移分布云图 (单位：m)

4.2 大坝应力分析

坝体沿坝轴线下游部分土体的最小主应力受浮托力和渗透力作用明显，在下游坝面约2/3坝高处有上抬趋势。坝体内最小主应力值普遍有所减少，断面上最小主应力最大值为0.711MPa，比竣工期最小主应力减少约0.224MPa，大主应力最大值为0.505MPa，应力分布云图如图5所示，随坝体断面高程增加而减小，渗透力和浮托力的影响对最大主应力的影响不大。计算结果符合工程规律。

图5 河床坝段应力分布云图 (单位:MPa)

4.3 大坝坝坡稳定性复核

仍然选取河床坝段、左岸岸坡坝段、右岸岸坡坝段进行计算分析。计算工况仍然按照前述的工况进行。假定若干滑面，将每个滑面所经过单元斜面上的剪应力按照有限元计算的节点应力插值积分求得，采用Bishop法进行滑裂面稳定性分析计算，带入相应的抗滑稳定计算公式求解稳定安全系数，经计算分析各断面稳定安全系数见表2。

表2 坝坡稳定安全系数统计

表2计算结果表明，凤亭河水库主坝在正常水位下工作，上、下游坝坡抗滑稳定安全系数K均大于规范规定的最小安全系数[K]，满足规范要求。

5 结 论

a.基于有限元分析的刚体极限平衡法结合圆弧滑动法进行分析，既考虑土体内部的变形也可以计算滑裂面上的稳定安全性，这和土体的实际工作状态相符。

b.风亭河水库在加固后，上游坝坡稳定安全系数普遍大于下游坝坡稳定安全系数，但均处于安全状态。

[1] 刘士乙,基于有限元极限平衡法的土工结构稳定分析研究[D].大连,大连理工大学,2015.

[2] 陈祖煜.土质边坡稳定分析—原理—方法—程序[M].北京:中国水利水电出版社,2001.