一种改进的渗流系数方法

一种改进的渗流系数方法

1变单元渗透系数法由于自由面界面的需要，需要对自由面界面进行迭代求解。求解该问题的有限元法以往采用传统的变网格法,该算法最大的优点是,渗流自由面和逸出点可以随着求解渗流场的迭代过程逐步稳定自行形成,迭代过程是收敛的。但该算法对复杂夹层和复杂排水系统的水工结构处理起来太困难了,几乎不可能实现;另外初始渗流自由面位置的假定要求也很高,如果初始位置与最终自由面位置相距甚远,则极易造成单元严重畸变,影响计算的精度;由于网格的变形,渗流与应力耦合分析也无法在同一个网格下进行,应用受到一定的限制。为了克服变网格法存在的问题,国内外学者致力于固定网格法的研究。自1973年Neuman提出用不变网格分析有自由面渗流的Galerkin法以来,出现了多种固定网格法。较有代表性的有剩余流量法、单元渗透矩阵调整法和初流量法。剩余流量法通过计算自由面单元内过自由表面的流量,来修改各节点的势,直至过自由面的流量小于某一允许值时为止。整个迭代过程依赖于第一次有限元计算的结果,计算精度较差。单元渗透矩阵调整法跨自由面单元按复合材料单元处理。单元通过节点与外部联系,其内部各点的参数可表示为节点处参数的多项式插值,是坐标的连续可导函数。复合材料单元渗透系数在复合缝面突变,其单元渗透矩阵不能反过来代表这一特性,不能真实反映这一区域的透水特性,且矩阵主元常不占优,影响总渗透矩阵的占优特性。采用复合材料单元,计算精度和计算稳定性都受到影响。初流量法在计算跨自由面单元的节点初流量时,自由面以下高斯点不予计算,实际上是把跨自由面单元当成复合材料单元处理,计算精度受到影响。本文在上述研究的基础上,提出了变单元渗透系数法。变单元渗透法系数是一种真正的固定网格法,迭代计算过程中对网格不作任何改动,只是在迭代过程中改变在自由面之上的单元的渗透系数,即可由程序自动求出自由面的位置。2固定网格法改变单元的渗透系数2.1me—数学模型的建立考虑土和水的压缩性,符合达西定律的二向非均质各向异性土体渗流问题的微分方程:∂∂x(Κx∂Η∂x)+∂∂y(Κy∂Η∂y)=Ss∂Η∂t(1)∂∂x(Kx∂H∂x)+∂∂y(Ky∂H∂y)=Ss∂H∂t(1)将整个渗流域离散为有限个单元,各单元内的水头函数可近似表示为:Ηe(x‚y)=me∑i=1ΝiΗiHe(x‚y)=∑i=1meNiHi式中:me——单元结点数;Hi——单元各结点水头;Ni——单元形函数。把单元e作为整个渗流域内的一个子域,并对该子域上的泛函式进行求导,则可得到方程组:Κ{Η}+S{∂Η∂t}+Ρ{∂Η∂t}={F}(2)K{H}+S{∂H∂t}+P{∂H∂t}={F}(2)对时间项取隐式有限差分,则式(2)变为:(Κ+1△tS){Η}t+△t+1△tΡ{Η}t+△t-1△tS{Η}t-1△tΡ{Η}t={F}(3)不计时间项且S和P矩阵均为零时,得到稳定渗流有限单元法计算公式:K{H}={F}(4)这样,把问题化为n介线性代数方程组的求解问题,通过求解上式可得到结点水头H。对于有压流,由于不存在自由面,边界是确定的,求解方程(3)可得到问题的解。由于土和水的压缩系数很小,在计算非稳定渗流时一般不考虑土和水对渗流自由面的影响,而是考虑由于自由面变化而产生的流量补给关系对于自由面的影响,即考虑给水度的影响。渗流域边界位置不能预先确定,因此确定渗流域边界是解题的关键,本次计算程序采用二维四结点等参单元和三角形单元进行计算。2.2单元渗流计算精度的确定变单元渗透系数法不同于别的固定网格法,它克服了别的固定网格法有时还需重新划分网格,迭代时需要调整自由面节点位置,增加计算工作量,容易造成单元出现畸异的缺陷。具体做法如下:1)由于自由面的位置不确定,在上游水位与下游水位之间先估计出初始自由面的区域。2)对全区域进行网格划分,在预估的初始自由面域内将网格划分得很密。第一次渗流计算时,全区域的渗透系数为给定的渗透系数。3)由于自由面上的节点的水头等于其位置势,自由面以上的节点水头小于其位置势,自由面以下的节点水头大于其位置势,将自由面之上的所有节点进行标识,然后判断处于自由面之上的所有单元。4)把位于自由面之上的单元的渗透系数改为一个很小的数,由于自由面之上的单元在第二次以后的迭代计算过程中,不需要参与计算,为了让自由面之上的单元对以后的计算不产生影响,方便编程序,所以把自由面以上的单元在迭代过程中都乘以一个很小的数,比如为1.0×10-18m/s,这么小的数,足以在计算过程中把自由面以上的单元渗流结果忽略,也就是说在计算过程中不再考虑自由面以上单元渗流的影响。5)将本次求出的节点势与上一次迭代求出的节点势比较,亦即判断:|Ηi+1j-Ηij|≤ε1(5)式中:i——迭代计算次数;j——节点号;ε1——同一节点两次计算的水头值的误差,在计算精度要求较高时,可将ε1取的小一些,如本文例子计算取10-4。若式(5)不满足,则重复步骤3)—5);若式(5)满足,则结束迭代。当渗流域内节点均满足式(5)时,渗流场中所有满足|Η-Ζ|≤ε2(6)的节点的连线即为自由面。在利用本方法进行渗流计算中,要求根据在已知边界水头连线方向估计渗流面的变化区域,对该区域的单元进行加密,单元愈密则计算精度愈高,由于在计算非稳定渗流时,网格加密的范围将很大,这增加了计算迭代的时间。本文提出的方法与计算饱和—不饱和渗流模型很相似,本方法在计算迭代过程中不考虑非饱和区的影响,但得出的结果还是令人满意的。3变单元渗透系数法的验证3.1计算结果以下游水位为例考察一各向同性且均质介质的河间地块,长100m,高80m,无源汇,无垂向补给,上游河流稳定水位为h1=60m,下游水位为h2=20m。整个计算区域的网格划分以及最终自由面的位置见图1和图2。计算结果与文献用改进初流量法计算的结果符合得很好,但是改进初流量法需要迭代10次时两次最大水头差才达到0.085m,而用改变渗透系数法计算,迭代7次时两次最大水头差已达到0.05m,由此可见,改变单元渗透系数法的效果是十分明显的。3.2全期水位上升过程中非稳定渗流的变化文献①介绍了一个室内非稳定模型试验,长×宽×高为315cm×23cm×33cm的均质砂槽物理模型。试验介质是平均粒径约为15mm的均质砂。取两种情况的试验结果与计算结果相比较。一是初期水位在离槽底10cm处,瞬时水位上升而产生的非稳定渗流变化;二是由于瞬时水位下降而产生的非稳定渗流的变化。计算域划分为180个三角形单元,计算中采用:μs=0.44,kx=kz=0.33cm/s,Ss=0。水位骤降时的浸润线试验与计算结果的比较见图3,图中浸润线的时间点从上到下依次为60、120、240、300和2400s。水位从10cm骤升至30cm时的浸润线试验与计算结果的比较见图4,图中浸润线的时间点从下到上依次为30、60、120、600和2400s。从中可以看出计算结果与实验结果非常一致。可见本文提出的方法对非稳定渗流计算也是非常有效的。3.3大坝初始条件以新疆伊犁特克斯河恰甫其海水利枢纽工程的粘土心墙土石坝为例,该工程位于特克斯河与小吉尔乃朗河交汇口东侧约20km处,是国家特大型水利枢纽工程。坝高96m,坝顶宽度为12m,上游边坡为1∶2.5,下游边坡为1∶1.7,心墙的坡率为1∶0.3。砂砾石坝壳的渗透系数为2.6×10-4m/s,心墙的渗透系数为1.0×10-7m/s,正常蓄水位为89.0m,正常坝后尾水位为10.05m。大坝模型见图5,共划分出4365个三角形单元、2089个结点。现以两种不同工况来模拟大坝的非稳定渗流,一种是上游水位下降引起的非稳定渗流,另一种是上游水位上升引起的非稳定渗流。上游水位从89.0m以速度3.0m/d下降到10.8m水位,下游水位保持在10.05m的水位。假定大坝底部按不透水边界处理,初始条件为大坝上游水位保持在89.0m时的稳定渗流。浸润线变化见图6,图中的数字是自由面位置的时刻(单位:d)。上游水位从0m以速度3m/d上升到89.0m水位。下游水位以1m/d速度上升到10.05m的水位。浸润线变化见图7,图中的数字是自由面位置的时刻(单位:d)。从图6和图7可以看出,上游水位下降或上升时,自由面在初期变化很快。由于坝壳与心墙内的渗透系数相差很大,上游水位下降时,水主要滞留在心墙内,心墙内自由面的下降变化很陡,自由面变化很慢,心墙内的孔隙水压力一直很高。水位上升时,上游坝坡的自由面变化比较平缓,但心墙内变化很陡,水基本上是从心墙内垂直下去的。通过非稳定渗流计算,动态了解坝壳、心墙等部位的渗透规律和浸润面变化规律,为大坝的设计以及动力计算分析提供定量化的依据。4单元渗透矩阵计算的一般规则通过对固定网格算法中的传导矩阵调整法的研究,本文提出了迭代收敛快、只需划分一次单元、结果稳定的改变单元渗透系数法。改