大型土石坝三维饱和开挖及渗流作用下围堰稳定性数值模拟研究

大型土石坝三维饱和开挖及渗流作用下围堰稳定性数值模拟研究

在节水工程中，洪水广泛存在于水库中。堤坝渗漏会导致管涌,引起滑坡甚至溃坝。长期以来,堤坝的渗流研究只局限于浸润面以下的饱和区,而忽略了浸润面以上水体运动状态的研究。事实上,上游水位的不断变化,导致饱和与非饱和区之间相互变化。如果忽略非饱和区而至分析饱和区水流的运动规律,经典的自由面逼近方法便不能计算非饱和区孔隙内的负压。在浸润面,即饱和区与非饱和区的交接面,通常作用着一个穿越它的局部加速度,故浸润面可解释为一阶加速度。渗透系数是含水量的函数,非饱和区的含水量与渗透系数和饱和区的渗透系数与储水系数同样重要,上述两个区域的水体是相互联系的。因此,只将饱和区作为研究对象是不全面的。有关该问题的计算方法大致可以分为两类:一类是解析解;一类是数值解。数值解通常有有限单元法或有限差分法。本文以有限差分法为基础,考虑流体在孔隙介质中的流动状态对孔隙介质的影响。有限差分法中在处理流体方面比较成熟的商业软件由美国Itasca公司研制,现对其饱和流—固耦合的基本方程做以简单介绍。1多孔介质流固耦合控制方程FLAC3D模拟多孔介质(如土体)中流体流动时,流体的模拟独立于力学计算。其主要通过孔隙水压力的消散引起土体中位移的变化,这一过程包含两种力学效果。第一,孔隙水压力的变化引起结构体中有效应力的变化;第二,孔隙水压力的变化又引起流体区域的变化。流体在孔隙介质中的流动依据Darcy定律,流固耦合过程满足Biot方程。流体在多孔介质中流动时,主要引起以下几个变量的改变:孔隙压力,饱和状态和渗透流量。这些变量通过流体质点平衡方程、Darcy定律间的相互关系来描述流体的流动。本构方程主要表现孔隙压力,饱和状态,体积应变化关系,近而实现流体——固体之间的耦合。下面将介绍一下FLAC3D中流体——固体之间相互耦合的差分方程:1.1外法线矢量的计算一般来讲,四面体的四个节点一般是指从1到4的四个节点号,面n相对于节点n,脚标f是面f相关的变量值。在四面体内假定孔隙水压力是按线性变化和流体密度为常量。根据孔隙水压力的节点值,利用高斯定律,水头梯度可写为(p-ρfxigi),j=-13V4∑l=1(pl-ρfxligi)n(l)jS(l)(1)(p−ρfxigi),j=−13V∑l=14(pl−ρfxligi)n(l)jS(l)(1)式中n(l)是面l的外法线矢量,S是四面体的面积,V是四面体的体积。为了提高精度值xi-x1i1i,式中x1i1i是四面体的一个角点坐标,将xi代入式(1),得到如下方程(p-ρfxigi),j=-13V4∑l=1p*ln(l)jS(l)(2)(p−ρfxigi),j=−13V∑l=14p*ln(l)jS(l)(2)式中p*l被定义为p*l=pl-ρf(xli-x1i)gi(3)p*l=pl−ρf(xli−x1i)gi(3)1.2局部节点压力对于全耦合流体(s是常数且等于1),流体的连续方程为:qi,i+b*=0(4)qi,i+b∗=0(4)式中b*=SΜ∂p∂t-q*v(5)b*=SM∂p∂t−q∗v(5)是体积力,也被用于力学分析的节点公式中,q*v=qv-α∂e∂t+β∂Τ∂t(6)q∗v=qv−α∂e∂t+β∂T∂t(6)首先,考虑一个四面体。利用这种算法,节点流量,Qnene[m3/s],n=1,4,等于四面体流量和被测体积的流体源,b*,则有如下的表达式:Qne=Qnt-q*vV4+mndpndt(7)Qne=Qnt−q∗vV4+mndpndt(7)式中Qnt=qin(n)iS(n)3(8)Qnt=qin(n)iS(n)3(8)mn=V4Μn(9)mn=V4Mn(9)四面体流量(式8)的各分量通过流动法则(式1)和水头压力梯度联系在一起。反之,水头压力梯度的各分量也可以通过四面体节点孔隙压力求得。为节约计算时间,局部矩阵被组装。单个区域的流量式通过流量求和得到Qnznz,式(8)给出了该方程的具体表述。局部区域矩阵[M]通过8个节点值{Qz}与8个节点压力水头{p*}得到。因为这些短阵是对称的,36个分量被计算;这些值在计算开始时被存储,且在大应变模式中每10步更新一次。根据区域矩阵的定义,有:Qnz=Μnjp*j(10)Qnz=Mnjp*j(10)式中{p*}是区域局部节点压力水头。总节点值QnΤnT通过重叠区域对节点的贡献计算得到,对于每个总节点n:QnΤ=Cnjp*j(11)式中[C]是总刚矩阵,{p*}节点压力水头的总矢量。根据流量平衡原理,有:-∑Qne+∑Qnw=0(12)式中∑代表所有区域在节点n的贡献之和。将式(7)和式(6)代入式(12)得到:dpndt=-1∑mn[QnΤ+∑Qnapp](13)式中QnΤ是节点孔隙压力的函数由式(11)和式(3)求出,∑Qnapp由下式给出:∑Qnapp=-∑[qvV4+Qw]n(14)在FLAC3D中,Biot模量是节点的属性,则式(9)可以写成下式:1∑mn=ΜnVn(15)式中Vn=∑(V4)n(16)将式(14)代入式(13),则有:ddt[pn-pnv]=-ΜnVn[QnΤ+∑napp](17)式中pnv=-ΜnVn[∑(α∈V4)n](18)下面就是求解差分方程式(17),此方程有两种解法:显式有限差分法和隐式有限差分解法,与普通有限差分解法相同。1.3节点内部流体的充放电测量对于饱和流体,在单个区域内的节点流量{Qz}和节点孔隙水压力{p}通过式(10)联系在一起,其矩阵形式可表示为:{Qz}=[Μ]{p-ρfxigi}(19)式中短阵[M]是区域几何形状和饱和度相关的已知矩阵,这一方程可以推广到非饱和流体,按下面方式进行改变:(a)考虑部分区域饱和,重力ρfxigi项乘以区域的平均饱和度ˉs。(b)节点流量要乘以相对渗透系数ˉk。(c)节点入口流速根据局部饱和状态确定。对于饱和或非饱和流体,流体连续方程的显式有限差分格式为:ΔpnΜ+Δss=-1snV[(QΤ+∑Qapp)Δt+s(∑αΔeVΤ4)](20)式中角标t和角标n被省略,VT是四面体的体积。在饱和时,s=1式(20)写成Δp=-ΜV[(QΤ+∑Qapp)Δt+s(∑αΔeVΤ4)](21)在非饱和状态时,p=0式(20)被写成如下:Δs=-1nV[(QΤ+∑Qapp)Δt+s(∑αΔeVΤ4)](22)这种技术利用了变饱和度的方法,在非饱和节点处(孔隙压力等于0处)其改变根据式(22),在饱和节点处(孔隙压力等于0处)孔隙压力的增加根据式(13)求得。2泰山污水储水站周围耳桥的饱和非饱和流固耦合值的模拟2.1围堰标准断面布设泰山抽水蓄能电站进、出水隧洞从京沪铁路、京福公路(104国道)底下通过。期间还将穿过断层破碎带,围岩受多条节理裂隙切割。其中围堰坐落在F4断层影响带,其施工导流采用全断面全年断流的土石围堰加高压射灌浆防渗心墙形式,围堰最大堰高10m,最大挡水水头约为9m,整个围堰要求在一个枯水期内建成。围堰轴线全长690.5m,围堰水深在标号S4处最深,约9m。围堰为临时构筑物,围堰类型的选择应以施工方便,拆除简单,就地取材,造价低为原则。根据施工方的围堰选型,采用高压喷射灌浆防渗心墙结合土工膜柔性防渗心墙围堰。围堰填料采用上游水库开挖的砂粘土、砂石土的混合料,个别块石粒径达到50cm。围堰填筑采用自然堆坡,迎水坡坡角为1ue5fe1.5,背水坡坡角为1ue5fe1.7。堰顶土体利用行走的重车压实。为了便于工程车的行走,堰顶顶部宽度为6m。堰体填料约为9.45万m3,高压喷射灌浆防渗心墙约为13362延米。土工膜施工约2054m2,混凝土底座约为102m3,图1给出了围堰标准断面图。以下计算确定围堰浸润线的位置,渗流场的渗透坡降、渗流量及其渗流对围堰和基础边坡的稳定影响;由于围堰及基础具有明显的三维特性,故对其进行三维流固耦合分析,同时也考虑围堰和基础共同作用下开挖对二者边坡的影响。2.2主上层基岩分布围堰堰体防渗结构上部采用防渗土工膜,围堰堰体下部和堰基采用高压喷射灌浆防渗心墙方案。围堰堰体由强透水的堆石料和石渣料堆筑而成。围堰处的河床覆盖层由沙砾石层组成,一般厚度2～5m,透水性强,但已被防渗墙完全隔断。其下的全风化或强风化岩层厚度一般在7～15m之间,其中防渗心墙深入其中4.2m左右,没有完全隔断下部的强风化层。根据地质勘察得到的吕荣值可以换算出基岩的渗透系数,并按照渗透系数值划分为弱透水层和强透水层。堰下基岩上部一定范围内属弱透水层,是整个围堰中渗透性最强的连通区域,因此渗流将主要在该部分基岩内发生。基岩分布具有明显的三维特征。主要材料性质见表1。设计单位在初设时,围堰下岩基不设置帷幕。本文渗流分析按岩基不设帷幕条件计算,取防渗墙嵌入基岩1.0m。3大坝三维界面分析3.1渗流分析的有限差分有限差分法计算深度取至基岩弱透水层一定深度,由于背水坡侧抽水后开挖约8m左右,故堰基取为12m,此时计算深度已至弱透水层一定深度,可视为不透水层;另外,堰基两侧也各取一定深度(基岩深度取12m),以考虑在基岩内出现绕岩体的渗流情形,更远处的绕流将十分微弱,从而可视为不透水层。渗流分析采用六面体单元,图2给出了含防渗心墙围堰的三维渗流分析的有限差分网格。围堰的三维弹塑性有限差分分析选用Mohr-Coulomb屈服准则,渗流分析选用Darcy定律。耦合分析根据采用前述的基本原理进行分析。计算范围以围堰轴线z轴,垂直地表方向为y轴,坐标原点选在S4处。考虑围堰主要用于挡水,故在S2～S4方向水深小于2.95m处围堰没有计算。S2～S4方向长取为245.07m至水深等于2.95m处;同理,由于S5～S6段水位比较浅,最深处2.95m,不会对围堰的稳定性造成危害,故S5～S6段没有被列入计算范围。S5～S4段长取为197.17m。围堰高度随低部基础的变化而变化。覆盖层也随地表高度的变化而变化,取其平均厚度为4m。以下均为风化层,取为8m。共18270个区域,16048个节点。浸润线的确定是通过透水地基上的均质土堤,由于受地基透水的影响,堤身浸润线降低。3.2基础开挖对围堰水压力的影响图3和图4分别给出了围堰及堰基在不开挖和开挖时最大水深处孔隙水压力和塑性区图。图中孔隙水压力等于零处即是浸润线的位置。图3给出了围堰在不开挖时的塑性破坏区和孔隙水压力的变化关系,从图中可知围堰由于受渗流的影响,而导致堰体在背水坡坡脚处发生塑性破坏,加之渗流影响,堰基处也发生了贯通破坏,和围堰的塑性破坏连在一起,导致围堰和堰基边坡共同形成正反向滑弧破坏。但在开挖侧预留平台没有形成贯通破坏。从孔隙水堰压力看,可以明显看出浸润线的位置,在浸润线以上,属于非饱和区,在浸润线以下属于饱和区。浸润线是饱和区与非饱和区的分界线,具有明显的空间分布特性。将其与基础开挖相比如图4,围堰塑性区没有太大的变化,说明围堰在不开挖时,结构变形破坏已基本完成,开挖对围堰及基础结构破坏的影响是很小的,主要是因为开挖结构位于浅地表的缘故。但预留平台的塑性破坏范围有所增加,导致平台形成贯通破坏。在开挖后,孔隙水压力产生了明显的变化,压力水头有明显下降,尤其是在靠近开挖边坡一侧,下降更为明显,表明这一区域流量明显加大。非饱和区有所增加,饱和区也随之减少。开挖后,渗流发生了变化,塑性破坏区在背水坡坡脚处和围堰本身出现,且形成贯通型破坏,说明背水坡侧形成滑坡破坏。在迎水坡侧,由于渗流和开挖的共同作用,导致塑性破坏区连在一起,因此形成围堰边坡的滑移破坏和基础边坡在预留平台出的小面积滑坡,这样围堰就很有可能产生劈裂破坏。所以,在围堰抽水开挖时,应加强背水坡脚和预留平台的加固;在迎水坡侧,应对坝体和基础的采取适当的加固措施,以防事故的发生。从计算的渗流量来看,在不开挖条件下,围堰最大剖面单位流量为6.27e-4m2/s,出口坡降为0.19。在开挖条件下,围堰最大剖面单位流量为1.06e-2m2/s,出口坡降为0.39。二者相差比较大,说明开挖对围堰及基础的渗流产生很大的影响。4“三维”流—结论基于非均质各向异性的饱和非饱和土的三维稳定渗流问题的基本控制方程,采用有限差分进行流——固耦合分析,对泰安抽水蓄能电站进出水口围堰和基础共同作用进行了三维流——固耦合分析。计算结果具有良好的规律性,从而表明这一方法特别适用于分析形状复杂、结构分区多且各分区渗透性差异较大的大型土石坝流固耦合分析。对泰安抽水蓄能电站进出水口围堰和基础共同作用进行了开挖和不开挖两种形式的三维流——固耦合分析。结果表明,基础开挖后,出口流量和坡降都有较大增加,根据规范如果在开挖后不设置防渗结构,围堰可能发生流图型破坏。因此有必要对围堰和基础设置防渗心墙。由于围堰是由沙砾石堆积而成。没有进行特殊处理,渗透性比较强,结构强度低,故应在背水坡坡脚处加设滤层或设置加固棱体。在迎水坡侧堰体和堰基出现反向划弧破坏,应在堰体和堰基中采用砂浆或混凝土封堵,以确保工程的安全。开挖对围堰的安全影响不大,主要是因为堰基位于地表,地应力对基