心墙堆石坝渗流-应力耦合作用下的渗流分析

心墙堆石坝渗流-应力耦合作用下的渗流分析

1应力场对渗流场的影响为了确保水库的安全运行，特别是水库的稳定运行，应特别考虑渗流问题，并进行科学、客观的渗流分析。目前,坝体渗流分析大多是独立于应力分析直接计算的,即计算时不考虑坝体应力改变和孔隙变形对渗流场的影响,这使计算过程简化,但不能反映在应力场作用下的渗流场。事实上,坝体内应力场和渗流场是相互影响的,渗流通过施加于某作用面上的渗透压力和在渗流区域内分布的渗流体积力而影响土体的应力分布,应力使土体产生形变、改变孔隙率或孔隙分布而影响土体的渗透性,从而影响渗流场,最终使坝体处于双场耦合平衡状态。因此,有必要将渗流-应力耦合作用纳入到坝体渗流场的分析中,研究应力场对渗流场的影响。当前岩土体渗流-应力耦合分析中,关于渗流场对应力场影响的研究较多,关于应力场对渗流场影响的研究较少,且基本都集中在裂隙岩体的应力变化对裂隙渗流的影响。在土体应力场对渗流场影响的研究方面,柴军瑞等提出均质土坝流固耦合分析的连续介质数学模型,推求渗透系数随应力场变化的规律;杨林德等推导饱和土体各向异性渗流直接耦合的有限元法计算公式,模拟应力作用下各向异性弹性多孔介质中的孔压变化;罗晓辉等指出,深基坑大变形耦合作用使得墙背土体的透水性能降低,孔压聚集;陈晓平等得出非均质土坝在渗流-应力耦合作用下的孔压分布和心墙附近的水力坡降;周创兵等指出,面板堆石坝在耦合作用下的流量变化及垫层、过渡层和堆石区最大水力坡降的变化。上述研究都指出了应力对渗流场的影响,但有的只考察了应力对渗透系数或孔隙水压力等渗流要素的影响,没有全面地反映应力场影响下渗流场的变化;有的考察了应力场影响下的多个渗流要素,但主要针对特定水位下的稳定渗流场,没有探讨蓄水位上升时,不断变化的坝体应力场对非稳定渗流场的影响。本文利用有限元方法计算分析心墙堆石坝分层填筑后蓄水过程中坝体渗流场分布,较为系统地研究和比较了各向同性和各向异性两种情况下不同蓄水时期坝体应力场对坝体孔隙水压力、渗透流速和坝体不同区域水力坡降的影响规律。2计算原理2.1增加一气多聚碳质权的本构方程Frelund和Rahardjo在Terzaghi土体有效应力原理基础上增加了气相,得到非饱和土体的本构方程为式中:{Δσ}为应力增量矩阵;[D]为土体本构矩阵;{Δε}为应变增量矩阵;其中H为与基质吸力相关的非饱和土结构模量;ua为孔隙气压力;uw为孔隙水压力;{Δua}为孔隙气压力增量。2.2基质吸力随时间的变化Frelund和Rahardjo研究指出,与饱和土体中渗流一样,非饱和土中渗流也是符合达西定律和连续性方程。只是在非饱和土中渗透系数不再是常数,而是随基质吸力变化的函数,土中体积含水率随时间变化而改变。考虑土中渗流的各向异性,将达西定律导入渗流连续性方程,可得到非饱和土的渗流控制方程为式中:kx、ky分别为x和y方向的渗透系数;γw为水的重度;θw为体积含水率;t为时间。2.3基本变量选取在流-固耦合计算中,土体平衡方程和渗流控制方程是同时求解的。应力场求解的是节点处的位移增量,渗流场中求解的是节点处的水头值,即耦合计算时联立式(3)和式(4)求解得到位移增量{Δδ}和孔隙水压力增量{Δuw}分别作为应力场和渗流场的基本变量。式中:[K]=[B]T[D][B]为刚度矩阵;[B]为应变矩阵;[Ld]=[B]T{m}N,为耦合矩阵;{m}T=1,1,1,0;N为形函数行向量;[ΔF]为节点外荷载增量。β为体变系数;[Lf]=∫NT{m}T[B]dV,为流体耦合矩阵;[Kf]=∫[B]T[Kw][B]dV,为单元刚度矩阵;[Kw]为渗透系数矩阵;ω为孔隙水压力系数;[MN]=NTN为质量矩阵。3计算模型、方案和参数3.1坝体材料设计本文的研究对象为心墙堆石坝,由于掺砾土心墙变形小、强度高,在国内外高坝中广泛使用,因此,本次研究中的心墙土料采用掺砾黏土料,坝壳为堆石料。为研究方便,将堆石坝体的反滤层和过渡层忽略,假设其与堆石的特性相近,设计一个坝高100m、坝体最大剖面如图1所示、仅由砾石土心墙和堆石坝壳组成的坝体。计算中将上述模型均匀划分为三角形和四边形单元,共2896个单元。3.2考虑应力场-固耦合-非耦合为了研究和比较各向同性和各向异性2种条件下坝体应力场对渗流场的影响规律,本文设计的计算方案分为考虑应力场的流-固耦合和不考虑应力场的非耦合渗流计算;在耦合与非耦合计算中,又分别考虑渗流的各向同性和各向异性,详见表1。3.3流-固耦合计算方法各种计算方案的初始应力场都为大坝分层填筑完成时的应力场和位移场。初始渗流场都为蓄水位为0m时对应的孔隙水压力场,非饱和区的负孔隙水压力场的绝对值随坝高线性增大,计算时未限定最大负孔压值。流-固耦合计算时,对于应力、应变分析,坝基设为x、y两向约束,坝顶和坝坡面都是自由边界。对于渗流分析,迎水面为变水头边界,蓄水位随时间从0m线性增加至90m(历时270d),水位以下部分的坝面作用有相应的静水压力。背水面为可变水头边界,即当总水头值低于高程时,为不透水边界;高于高程时,为自由出流边界。坝顶和坝基为不透水边界。在非耦合计算中,渗流初始条件和边界条件与耦合计算中相同。3.4土-水特征曲线本次研究中坝体土料采用摩尔-库仑弹塑性模型,坝体土料特性参数如表2所示。坝料的土-水特征曲线(SWCC)是参考其他学者对类似土料的室内试验成果取得;渗透函数根据vanGenuchten提出的由土-水特征曲线推算方法求得。计算所用坝料的土-水特征曲线、渗透函数如图2所示。应力场对渗流场的影响可通过渗透系数随应力的变化来反映。考虑应力影响的渗透系数k可由下式计算:式中:κ(σ′)为渗透系数修正函数;σ′为土体有效应力;k0为未考虑应力影响的渗透函数。应力的变化会引起土体体积的改变,进而改变土体的孔隙比和渗透性。由室内单向固结试验结果,可用下式求得土体的渗透系数k:式中:vc为固结系数;av为压缩系数;0e为初始孔隙比。由压缩定律可知:式中:e为土体孔隙比。不考虑应力影响时土体的渗透系数k0为初始孔隙比e0对应的渗透系数,即由式(7)、(8)可得渗透系数修正函数κ(σ′)为式中:av0为初始压缩系数。计算时采用的渗透系数修正函数如图3所示。竖直和水平方向的渗透系数均按上述函数进行修正。4结果和讨论4.1耦合与非耦合分析图4为渗流各向同性(kx=ky)、耦合与非耦合情况下蓄水至不同时刻时坝体孔隙水压力分布图。图5为各向同性坝体不同蓄水时刻的有效应力场分布。由图4可看出,不同时刻时耦合计算所得孔压分布与非耦合不同。蓄水初期耦合计算得出的孔隙水压力明显低于非耦合结果(见图4(a)),随着水位上升两者的差值逐渐减小,当蓄水至90m时,耦合与非耦合计算得出的饱和区孔隙水压力场基本相同(见图4(c))。这是由于蓄水初期,上游坝壳浸水少,坝体中的孔隙水压力小,坝体竖向有效应力较大,土体压密变形大,因此受应力影响的渗透系数显著低于不考虑应力的渗透系数值,致使耦合分析所得的浸润线较低,孔隙水压力也较低;随着水位上升,上游坝壳和心墙内的竖向有效应力大幅减小(见图5),受应力影响的渗透系数与不考虑应力的渗透系数值相差较小,因此,耦合与非耦合结果相近。图6为渗流各向异性(kx=2ky)条件下蓄水至不同时刻时坝体孔隙水压力分布图。由图可看出,应力场对坝体孔压的影响规律与各向同性时类似。各向异性条件下,水头前锋行进较各向同性条件快,在蓄水后期考虑坝料的各向异性时耦合作用对孔压影响没有各向同性时显著。4.2心墙内各渗流速度的变化规律图7是同一时刻坝体A-A截面(见图1)上不同点的流速随高程变化曲线。从图可看出,各向同性和异性(kx=2ky)条件下应力场对流速影响规律相似,考虑应力场影响均使得渗透流速减小,低于非耦合结果。在自由面处即饱和区和非饱和区的交界面上,流速都出现突变现象,当x方向的渗透系数大于y方向时,渗透流速较大,自由面处渗透流速的突变量较小。图8是坝体不同点的流速随蓄水时间的变化曲线。由图8(a)可看出,在各向同性(kx=ky)和异性(kx=2ky)两种情况下,心墙底部D点(见图1)流速变化规律类似,水流前锋到达D点后该点流速逐渐增长。由图8(b)可看出,在各方案下心墙墙趾E点流速变化规律亦类似。水流前锋到达时流速从0陡增至峰值,随着时间增长流速下降且逐渐趋于稳定。耦合计算所得D、E点流速变化相对非耦合结果有一定的滞后,且渗流稳定后耦合作用下的流速小于非耦合流速。这是由于在应力作用下坝体压密使渗透流速减小,致使渗流量减小,延缓了水流前锋行进。4.3各向异性条件下心墙渗流出逸点的水力坡降特征水力坡降是影响坝体渗流稳定和渗透破坏的重要因素,因此,有必要研究坝体应力场对水力坡降的影响规律,为科学合理的渗流稳定评价提供参考。图9是各向同性和异性条件下坝体不同部位点水力坡降随时间变化图。由图可看出,在各向同性(kx=ky)和异性(kx=2ky)两种情况下,坝体中同一点的水力坡降随时间的变化规律相差较小。上游坝壳底部中心A点(见图1)的水力坡降从0陡增至峰值,然后迅速降低。当水流前锋达到时,心墙底部中心B点的水力坡降先陡增,然后逐渐增大,与此点邻近的D点的渗透流速的变化规律类似。心墙下部水流出逸点C的坡降先急剧升高到峰值,然后略降后趋于稳定。在水流前锋突至时,A点耦合计算所得水力坡降峰值明显高于非耦合解。水流前锋达到时,水力坡降会突然升高,这是许多大坝在蓄水初期易发生水力劈裂的原因之一。随着蓄水时间的增长,蓄水后期上游坝壳和心墙中耦合和非耦合所得的水力坡降趋近一致。而对于心墙和下游坝壳交界面水流出逸点C,耦合计算得出的水力坡降仍比非耦合结果高(见图9(e)、(f)),因此,不考虑应力影响的渗流分析得出的出逸点的水力坡降值偏小,对渗透破坏的判别来说偏于危险。心墙下游出逸区域是影响心墙坝渗透稳定的关键区域,因此,对该区域的水力坡降峰值及其变化规律应进行更深入的分析研究。图10为不同蓄水位下心墙出逸区最大水力坡降值imax。图中各方案下imax随蓄水位大致呈线性上升趋势,耦合作用下imax明显高于非耦合结果。同一蓄水位下,各向同性条件下imax高于各向异性的结果,从危险区域的水力坡降峰值看,渗流各向同性坝体比各向异性更易发生破坏。总体来说,各向同性和异性条件下非耦合作用下坝体不同区域点水力坡降,尤其是心墙下游渗流出逸区域的水力坡降峰值均明显低于耦合解,这说明在渗流分析时考虑应力的影响是必要的,否则会得出偏于危险对工程不利的结果。需要指出,本文对坝体渗流场-应力场进行二维数值分析有一定的局限性,渗流各向异性研究时只能考虑平面上的各向异性,不能完全反映实际情况。5应力场作用与非耦合作用渗流场的影响(1)坝体在蓄水初期耦合计算得孔隙水压力均明显低于非耦合结果,随着水位上升两者的差值逐渐减小,当蓄水位接近坝顶时,耦合与非耦合作用下饱和孔隙水压力场基本相同。渗流各向异性条件下,坝体应力场对渗流场的影响规律与各向同性相似。(2)渗流各向同性和异性坝体在