考虑渗流-应力耦合的均质土坝边坡稳定分析

聂相田, 庄濮瑞, 焦延涛, 王博

(1.华北水利水电大学 水利学院,河南 郑州 450046; 2.水资源高效利用与保障工程河南省协同创新中心,河南 郑州 450046; 3.河南省水环境模拟与治理重点实验室,河南 郑州 450046)

渗流问题是影响土坝安全的重要因素之一。只有深入分析土坝的渗流分布,才能确保土坝运行安全[1]。而土坝的渗流场与应力场相互作用,是一种耦合关系[2]。在分析土坝边坡的稳定性时,若只考虑单独的渗流场或应力场进行计算,得到的结果必然与实际情况存在差异。近年来,已有许多学者在土坝边坡或者其他岩质边坡的稳定性分析中,考虑渗流场与应力场的耦合作用,并取得了一些研究成果。如苗丽等[3]基于多孔介质的渗流特性和土的非线性本构关系,研究了渗流场与应力场的耦合作用,并对某非均质土坝的边坡稳定情况进行了计算分析。李宗坤等[4]探讨了渗流场与应力场相互作用的机理,提出了考虑渗流与应力直接耦合作用的边坡稳定分析方法。贾善坡等[5]将损伤力学理论引入到渗流-应力耦合计算分析中,对比利时核废料库开挖过程中围岩的损伤演化以及渗流场和应力场耦合过程进行了计算分析,得到了围岩损伤特性、孔隙压力以及渗透性的变化规律。柳厚祥等[6]考虑应力场与渗流场的耦合特性,对某尾矿坝的非稳定渗流进行了计算分析。宋传旺等[7]采用有限元软件MIDAS建立了某尾矿坝的计算模型,探究了应力场和渗流场耦合作用对该尾矿坝稳定性的影响。沈振中等[8]应用迭代解耦技术建立了基于无单元法的应力场与渗流场耦合分析计算模型,并研制开发了计算程序。

综上所述,目前关于岩土体边坡稳定性分析的研究成果中,考虑渗流-应力耦合影响的已有不少,而且多数是结合有限元强度折减法进行计算的。采用有限元强度折减法进行边坡稳定性计算时,边坡失稳判据大致有3种：判据一,数据计算是否收敛;判据二,塑性区是否贯通;判据三,特征点位移突变[9]。根据上述判据,目前大多采用Drucker-Prager模型或Mohr-Coulomb模型进行岩土体的位移和应力计算,如文献[4][9][10]等。尽管采用Drucker-Prager模型或Mohr-Coulomb模型等理想弹塑性模型能很好地描述岩土体的非线性特性,但此类模型存在一些缺陷,如屈服面存在尖角，导致计算烦琐,且收敛速度缓慢,甚至不收敛[4]。因此,当采用此类模型进行边坡稳定性分析,并结合判据一和判据二进行边坡失稳破坏判断时,容易误判,甚至结果失真。而邓肯-张(Duncan-Chang)模型是一种基于三轴试验数据得出的增量弹性模型,不仅适用于描述土体的非线性变形特征,而且可以在一定程度上反映土体的弹塑性变形,同时,该模型所需参数不多,各参数的物理意义明确,获得参数的途径亦较简单,计算也很容易收敛。因此,邓肯-张模型在岩土体的应力变形计算分析中得到了广泛的应用。基于邓肯-张模型的边坡稳定有限元强度折减法也有学者进行过研究[11-12],但考虑渗流-应力相互耦合、相互影响,并基于邓肯-张模型的边坡稳定有限元强度折减法的研究则相对较少。因此,本文基于ANSYS软件,利用其APDL命令流语言并结合Fortran语言,编制模块化的渗流-应力耦合计算程序,并基于邓肯-张模型的边坡稳定有限元强度折减法对某均质土坝的边坡稳定性进行计算分析,以期为土石坝的边坡稳定分析提供一些参考。

1 计算原理

1.1 渗流计算原理

ANSYS软件本身并没有集成渗流计算模块,但集成了温度场分析模块,鉴于温度场的计算原理与渗流场的相似,故可以借用其温度场分析模块来进行渗流场的计算分析。关于渗流场与温度场计算原理的相似性,文献[1]和[13]均从基本理论、微分方程、初始边界条件3个方面进行了详细的论证,本文不再赘述。

1.2 应力变形计算原理

目前,描述岩土体应力-应变关系的数学模型很多,主要包括弹性模型和弹塑性模型两大类,但真正在工程实践中得到广泛应用的却为数不多。其中,邓肯-张模型因其具有多种优点而获得了最为广泛的应用,它包括初始的E-v模型和修正后的E-B模型。目前ANSYS软件的各版本中尚未嵌入邓肯-张模型,但由于其计算公式简单,可以通过APDL语言进行编程实现。本文即通过APDL语言创建邓肯-张E-v模型的宏文件,然后通过调用该宏文件来实现土体的应力计算。关于邓肯-张模型的详细介绍,可参考文献[1]和[13]。

1.3 渗流-应力耦合计算原理

在渗流过程中,渗流给土体介质施加渗透作用力,应力或应变引起其渗透性质变化,继而影响渗流场的分布,从而反过来影响岩土体的应力场[14]。目前,进行渗流-应力耦合计算大致有两种途径：直接耦合和间接耦合。直接耦合计算需求出包含有渗流场与应力场相关未知量的复杂数学模型的解析解,这是很难实现的。因此,本文采用间接耦合方法,通过迭代计算来求得渗流-应力耦合作用下的渗流场和应力场。

1.3.1 渗流场对应力场的影响

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

(1)

p=γw(H-y)。

(2)

式中：fw为节点所受面力;NT(x)为形函数;Γ为渗流边界;p为节点静水压力列阵,γw为水的容重;H为节点水头列阵;y为节点y向坐标列阵。

动水压力可视为与自重类似,计算时按体力计入节点荷载,即：

fs=∬ΩNT(x)fpdΩ，

(3)

(4)

式中：fs为节点所受体力;Ω为渗流计算区域;fp为动水压力。

1.3.2 应力场对渗流场的影响

应力场的变化会引起岩土体的体积及孔隙发生变化,从而导致其渗透系数发生改变。应力与渗透系数的关系用文献[14]推荐的关系式来表达,即：

K=K0exp[-β(σ-pe)]。

(5)

式中：K为土体单元的渗透系数列阵;σ为土体单元的应力列阵;pe为土体单元的静水压力列阵,本文取每个单元4个节点静水压力的平均值作为土体单元的静水压力值;K0为土体的初始渗透系数;β为经验系数,取为3e-7。

1.3.3 耦合分析求解步骤

步骤1选取一个初始渗透系数,利用有限元方法求解渗流场,得区域内的水头分布H(x,y,z)。

步骤2将渗流场分析得到的水头H(x,y,z)分别代入式(2)与式(4)中,求解各节点的静水压力p和动水压力fp。

步骤3将各节点的静水压力p和动水压力fp分别代入式(1)和式(3),根据有限元法求解应力场,获取应力分布σ(x,y,z)。

步骤4根据各单元4个节点的静水压力值求得各单元的静水压力值pe,将σ(x,y,z)和pe代入式(5),求解新的渗透系数,并替代初始的渗透系数,重新进行渗流场的求解。

步骤5重复步骤1—4，进行迭代计算,直至满足如下精度要求：

|Hn+1(x,y,z)-Hn(x,y,z)|≤εH。

(6)

式中：Hn+1(x,y,z)和Hn(x,y,z)分别为第n+1次和第n次迭代求得的渗流场水头分布;εH为水头的求解精度。

1.4 基于邓肯-张模型的强度折减法

有限元强度折减法得到的安全系数是基于强度储备概念的安全系数,黏聚力c、内摩擦角φ按下式进行折减[10]：

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

(7)

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

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

(8)

(9)

式中：Ei为初始弹性模量;Rf为破坏比;(σ1-σ3)f为破坏时的主应力差,可由Mohr-Coulomb定律推导得到。

将式(9)及Ei的表达式代入式(8),可得切线模量的最终表达式为：

(10)

式中：k为初始切线模量;n为模量指数;Pa为大气压强,取101.4 kPa。

从式(10)可以看出,只要对式中的黏聚力c和内摩擦角φ进行折减,即可实现与Drucker-Prager模型或Mohr-Coulomb模型相似的材料强度折减效果,从而保证有限元强度折减法的顺利实现,也就是说,基于邓肯-张模型的有限元强度折减法在原理上是可行的。

2 工程实例分析

2.1 工程实例

王庄水库大坝为均质连续的各向同性土坝,坝高22.50 m,坝顶高程613.95 m,坝底高程591.00 m。坝顶总长128.00 m,坝顶宽4.00 m,大坝上游坡比为1∶2.5,下游坡比为1∶2,下游设贴坡排水。校核洪水位612.50 m,设计洪水位611.40 m,正常蓄水位609.50 m,3种工况下下游均无水。限于篇幅,本文中只给出运行期正常蓄水位工况的计算结果,所谓运行期即指水库蓄水后坝体和坝基形成稳定的渗流场。因此,本次计算假定土体为饱和土体,土料的饱和渗透系数取为2.23×10-8m/s。对于稳定渗流计算而言,浸润面以上节点对于流场计算贡献不大,可称之为虚点,本文参考相关文献[1]的经验,采用ANSYS提供的死活单元技术,将处于浸润线上部的单元“杀死”,使其不参与计算,而只“激活”浸润线下部的单元网,当进行应力场分析时重新“激活”浸润线上部单元,使其参与计算。参考相关文献[3,15]及试验资料,本文采用邓肯-张E-v模型进行应力计算,所用的饱和土体材料参数见表1。

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

2.2 渗流-应力耦合计算结果分析

进行渗流-应力耦合计算时,为节省计算时间,本文只建立坝体的有限元模型。渗流计算边界条件取为：坝体底部视为不透水边界；左侧为上游水头边界；右侧无水,因此不加水头,但右侧斜边视为可能溢出面。应力计算部分边界条件为：坝体底部为x、y向位移约束。模型的坐标系选用笛卡尔直角坐标系,x轴指向下游为正,y轴向上为正。

在正常蓄水位工况下,渗流-应力耦合计算时坝体y向位移云图与非耦合计算时坝体y向位移云图的对比如图1所示。从图1中可以看出,考虑渗流-应力耦合作用后,坝体的y向位移明显增大。这是因为考虑渗流-应力耦合作用时,一方面,渗透水压力的竖向分力使坝体的土体有下沉趋势;另一方面,地下水动水压力作用于土体骨架,改变了原有土体颗粒结构的排列,从而导致坝体的y向位移明显增大。

图1 正常蓄水位工况下坝体y向位移

2.3 边坡稳定计算结果分析

根据前文所述的基于邓肯-张模型的有限元强度折减法原理,进行坝体边坡的稳定计算。考虑到边坡失稳判据一和判据二的特点及局限性,并考虑到失稳判据与邓肯-张模型的相容性,本文采用判据三(即特征点位移突变判据)来进行边坡的失稳破坏判定。采用判据三时,需确定特征点,并将其水平或竖向位移作为位移突变的判据。为此,本文选取坝顶的72号节点作为特征点,用以监控坝体位移的改变。72号节点具体位置如图2中A点所示。

采用强度折减法计算时,首先可使折减系数以0.1为间隔进行变化,大概确定安全系数范围;然后再以0.01或更小的间隔折减土体的黏聚力和内摩擦角,直到找出更精确的安全系数。

根据上述方法得到的特征点竖向位移与安全系数的关系曲线如图3所示。从图3可以看出：

1)考虑渗流-应力耦合计算时,若以0.1的间隔改变折减系数,当折减系数介于1.2和1.3之间时，特征点的竖向位移有明显的突变。为了找出更准确的安全系数,在1.2与1.3之间以0.025的间隔加密折减系数,安全系数在1.225与1.250之间时有明显的位移突变。因此,当考虑渗流-应力耦合作用时,可确定本文所研究的均质土坝的边坡稳定安全系数为1.225。

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

2)不考虑渗流-应力耦合计算时,当折减系数介于1.3和1.4之间时特征点的竖向位移有明显的突变。同样,对折减系数以0.025的间隔进行加密,安全系数在1.325与1.350之间时特征点有明显的竖向位移突变。因此,当不考虑渗流-应力耦合作用时,可确定本文所研究的均质土坝的边坡稳定安全系数为1.325。

3)考虑渗流-应力耦合作用后,坝体边坡的稳定安全系数较不考虑渗流-应力耦合作用时有所减小。原因应为考虑耦合作用时,一方面,渗流场对土体产生的渗流力在计算时直接加入到坝体原有的应力场里面进行计算,渗流力加大了坝坡的滑移力,从而降低了大坝的抗滑稳定性;另一方面,受渗流场的影响,土体出现软化,更容易发生变形,对坝坡抗滑稳定不利。因此,这也说明了考虑渗流-应力耦合作用时,渗流场对坝体边坡的稳定是不利的。

3 结论

基于渗流场与应力场的相互作用原理,本文利用ANSYS软件,提出了考虑渗流-应力耦合作用下基于邓肯-张模型的有限元强度折减法的计算方法。计算结果表明：

1)渗透水压力的竖向分力使坝体的土体有明显的下沉趋势,因此对土石坝进行应力-变形计算分析时,应充分考虑渗透水压力的影响。

2)进行边坡稳定计算分析时,考虑渗流-应力耦合作用计算得到的边坡稳定安全系数明显比不考虑渗流-应力耦合作用的小,渗流场对坝体边坡的稳定是不利的。因此,在进行土石坝边坡稳定计算时,应尽量考虑渗流-应力的耦合作用,使计算结果更为准确。