考虑渗流—应力耦合的黏土心墙坝边坡稳定分析

2021-05-31崔志瑞郭莹莹

广东水利水电 2021年5期

黄香，崔志瑞，郭莹莹，李颖

(1.三峡金沙江云川水电开发有限公司，云南昆明 650224；2.华北水利水电大学，河南郑州 450046)

1 概述

对于黏土心墙坝而言，渗流与稳定问题关系到大坝能否安全运行[1]。而边坡稳定分析是黏土心墙坝设计的基础，通过边坡稳定分析可反映出大坝设计方案的合理性，为大坝结构设计提供科学依据[2]。目前关于黏土心墙坝边坡稳定的研究成果较多，如文献[3-4]等均采用传统的边坡稳定计算方法，对一些黏土心墙坝的边坡稳定情况进行了计算分析。许多研究成果[5-6]表明，在天然状态下，渗流与土体孔隙空间相互关联，渗流场与应力场是一种耦合关系。对土石坝(包括黏土心墙坝)进行边坡稳定分析时，只考虑单独的渗流场或应力场进行计算，得到的结果会与实际情况存在差异，采用传统的边坡稳定计算方法，会存在一些不足。另一些文献[7-8]考虑了渗流—应力双场耦合的作用，并对黏土心墙坝的应力变形分布规律和工作性态进行了计算分析，但对双场耦合作用对边坡稳定的影响并没有进行更深入的研究。综合分析文献可知，目前考虑渗流—应力双场耦合的土石坝边坡稳定分析，较多研究偏重于尾矿坝或均质土坝[9-10]，而考虑双场耦合的黏土心墙坝边坡稳定的研究则比较少。

另外，从研究方法层面，目前大多数研究采用有限元强度折减法进行土石坝边坡稳定计算分析。通过此法进行边坡稳定计算时，边坡失稳判据大致有三种：判据一，数据计算是否收敛；判据二，塑性区是否贯通；判据三，特征点位移突变[11]。根据判据进行土石坝边坡稳定计算分析时，大多数研究采用Drucker-Prager模型等弹塑性模型进行坝体的位移和应力计算[9-11]。尽管上述理想弹塑性模型能很好描述岩土体的非线性特性，但此类模型如屈服面存在尖角导致计算繁琐，且收敛速度缓慢，甚至不收敛的缺陷[9]。因此，当采用此类模型进行边坡稳定计算分析，并结合判据一和判据二进行边坡失稳破坏判断时，容易造成误判。而邓肯—张模型做为一种增量型弹性模型，可克服上述弹塑性模型计算时不易收敛的缺陷，而且可以在一定程度上反映土体的弹塑性变形，目前邓肯—张模型在岩土体的应力变形计算分析中得到了广泛的应用，但考虑渗流—应力双场相互耦合，相互影响，并基于邓肯—张模型的边坡稳定有限元强度折减法的研究则相对尚较少。

综上所述，本文利用ANSYS软件的APDL命令流语言并结合Fortran语言，编制模块化的渗流—应力耦合计算程序，并基于邓肯—张模型的边坡稳定有限元强度折减法，对某黏土心墙坝的边坡稳定进行计算分析，以期为黏土心墙坝的边坡稳定分析提供借鉴和参考。

2 计算原理

2.1 渗流—应力耦合计算原理

目前，渗流—应力双场耦合计算方法主要为直接耦合法和间接耦合法，考虑直接耦合计算时需建立复杂的数学模型且求解困难，因此，本文采用间接耦合方法进行双场的耦合计算。

1)渗流场对应力场的影响

渗流场对应力场的影响有动水压力(渗透体积力)和静水压力(渗透压力)两种形式。静水压力(渗透压力)为作用与接触面上的力，计算时可按面力作用形式计入节点荷载，即：

(1)

式中NT(x)为形函数；Γ为渗流边界；p为渗透压力，其可表示为：

p=γw(H-y)

(2)

式中γw为水的容重；H为水头；y为位置y向坐标。

动水压力(渗透体积力)可视为与自重类似，计算时按体力计入节点荷载，即：

fs=∬ΩNT(x)fpdΩ

(3)

式中NT(x)为形函数；Ω为渗流计算区域；fp为动水压力，其可表示为：

(4)

考虑渗流场影响的应力场控制方程，则可表示为：

(5)

2)应力场对渗流场的影响

在荷载作用下，应力场的变化，会导致岩土体的体积及孔隙发生变化，从而导致其渗透系数发生改变。应力与渗透系数的关系可用文献[6]推荐的关系式表达，即：

K=K0exp[-β(σ-p)]

(6)

式中σ为土体应力；K0为初始渗透系数；β为经验系数，取为3e-7；p为静水压力。

考虑应力场影响的渗流场控制方程，可表示为：

(7)

式中H为水头分布；K(σij)为渗透系数，可由式(6)表示；Ω为渗流计算区域；Γ1为已知流量边界；Γ2为自由面边界；Γ3为已知水头边界；n2和n3为对应边界的法线方向。

2.2 基于邓肯—张模型的强度折减法原理

利用有限元强度折减法时，黏聚力c、内摩擦角φ按下式进行计算：

cf=c/Fs,φf=arctan(tanφ/Fs)

(8)

式中Fs为折减系数(即安全系数)；cf、φf分别为折减后的材料内聚力和内摩擦角。

对于邓肯—张模型，其切线模量Et可由下式表示：

(9)

式中Ei为初始弹性模量；Rf为破坏比；两者的具体表达式可参见文献[5]。而上式中(σ1-σ3)f表示破坏时的主应力差，可由Mohr-Coulomb定律推导得到，具体为：

(10)

将式(10)及Ei的表达式代入式(9)，可得切线模量的表达式为：

(11)

从式(11)可以看出，只要对式中的黏聚力c和内摩擦角φ进行折减，即可实现与Drucker-Prager模型或Mohr-Coulomb模型相似的材料强度折减效果，另外本文在进行强度折减时，还考虑将初始切线模量k进行折减，以充分满足土体的变形性质，即：

kf=k/Fs

(12)

将式(8)和式(12)代入式(11)，即可得到考虑有限元强度折减法的邓肯—张模型切线模量计算表达式，如式(13)所示：

(13)

3 工程实例及结果分析

3.1 工程实例

某黏土心墙坝，工程总库容为110.87万m3，工程规模为小(1)型，工程等别为Ⅳ等，永久性水工建筑物级别为4级，临时性水工建筑物级别为5级。跨河布置总长为185.00 m，坝顶高程为1 804.30 m，坝顶设高1.2m混凝土防浪墙，墙顶高程为1 805.50 m，坝顶宽度为5.00 m。大坝位于河谷处，最大坝高断面及材料分区如图1所示。限于篇幅，本文中只给出正常蓄水工况的计算结果，大坝正常蓄水位为1 800.08 m。根据实验资料[12]，本次计算分析时，坝体材料采用邓肯—张E-v模型，具体材料参数见表1，坝基及帷幕各材料假定为线弹性材料(见表2)。

图1 黏土心墙坝计算断面及材料分区示意(单位：m)

表1 坝体邓肯—张E-v模型材料参数

表2 坝基弹性模型材料参数

3.2 计算结果分析

进行渗流—应力耦合计算时，计算断面的有限元网格如图2所示。有限元模型的顶边界截至坝顶高程为1 840.30 m，底边界截取至坝基帷幕最深处以下，至高程为1 680.00 m处。本次有限元模拟计算，渗流计算边界条件取为：坝基底部为不透水边界，左侧为上游水头边界，右侧无水，因此不加水头，但下游侧斜边视为可能溢出边界；应力计算部分边界条件为：坝基底部为x、y向位移约束，坝基上下游侧x向位移约束。模型的坐标系选用笛卡尔直角坐标系，x轴指向下游为正，y轴向上为正。

图2 计算断面有限元网格示意

1)渗流计算结果分析

从图3可以看出:考虑和不考虑耦合作用时，坝体的水头等值线分布规律大致相同。但从中也可看出考虑耦合作用时，水头为90.737 m和92.884 m的等值线更靠近下游，这是由于考虑耦合作用后，受应力场的影响使材料的渗透系数有所增大，心墙削减水头的能力减弱。

a 考虑耦合作用

b 不考虑耦合作用

2)应力变形计算结果分析

从图4可以看出，考虑耦合后坝体的y向位移与不考虑耦合时相比变大，这是由于耦合作用渗透水压力的竖向分力使坝体的土体有下沉趋势，导致考虑耦合效应后，坝体的y向位移增大。

从图5可以看出，考虑耦合和不考虑耦合情况下，坝体的y向应力均符合一般规律。但不考虑耦合情况下，心墙y向应力小于两侧坝壳料的应力，具有拱效应，这是因为心墙的变形模量较两侧坝壳料小，心墙位移大，两侧坝壳料位移小，心墙部分自重传递到坝壳，反映了心墙内部应力小于两侧坝壳料的应力。考虑耦合情况下，心墙部位也产生了明显的拱效应，但考虑耦合情况下，心墙上下游两侧应力分布不成对称分布，且心墙上下游两侧同一高程处，上游侧y向应力大于下游侧，这是由于耦合情况下，考虑了渗透水压力的影响，上游水头较高，渗透水压力较大，导致上游土体压力大于下游土体的压力。

a 考虑耦合作用

b 不考虑耦合作用

a 考虑耦合作用

b 不考虑耦合作用

3)边坡稳定计算结果分析

考虑边坡失稳判据1和判据2的特点及局限性以及失稳判据与邓肯—张模型的相容性，本文采用判据3对边坡的失稳破坏进行判定。采用判据3时，需确定特征点后将水平或竖向位移作为位移突变的判据。因为黏土心墙变形比两侧坝壳料大，更易导致边坡的失稳破坏，故本文选取心墙顶的9 672号节点作为特征点(图2中A点)。在开始计算时，折减系数以0.1为间隔进行变化，以找出安全系数的大致范围，再以0.025或更小的间隔进行折减，直到找出更精确的安全系数。根据上述方法得到的特征点竖向位移与安全系数的关系曲线如图6所示。

图6 特征点竖向位移与安全系数的关系曲线示意

从图6可以看出，当考虑耦合时，折减系数由1.2增大为1.3时，特征点的位移有明显的突变。对折减系数进行加密计算，当安全系数在1.250～1.275之间时，有明显的的位移突变，因此，考虑耦合计算时，本文所研究黏土心墙坝的边坡稳定安全系数为1.250。同样可以看出，当不考虑耦合计算时，当折减系数由1.3增大为1.4时，特征点的位移有明显的突变。同样，对折减系数进行加密，可以看出安全系数在1.300～1.400之间时，有明显的的位移突变，因此，不考虑耦合计算时，可确定本文所研究黏土心墙坝的边坡稳定安全系数为1.375。从图6还可以看出，考虑耦合后计算所得坝体边坡的稳定安全系数比不考虑渗流—应力耦合时小，究其原因是，一方面考虑耦合作用时，渗透水压力竖向分力使坝体的土体有下沉趋势，渗流力加大了坝坡的滑移力，从而增大了坝体边坡土体的下滑趋势，另一方面考虑耦合作用时，受渗流场影响，土体的变形模量变小，土体软化，更易产生变形，对坝坡抗滑稳定不利。