饱和 非饱和土体非稳定渗流数值分析

　1999年12月

水　　利　　学　　报

SHU IL I 　　XU EBAO

第12期

饱和2非饱和土体非稳定渗流数值分析

Ξ

吴梦喜　高莲士

(清华大学水电系,北京,100084)

摘　要　本文对一般的非饱和渗流有限元计算方法加以改进,有效地消除了非饱和渗流数值计算中存在的数值弥散现象.同时还提出了一种简便有效的逸出面处理新方法,并给出了非饱和非稳定渗流计算的实例.关键词　非饱和,非稳定,渗流,逸出面,数值方法.中图号　TV139114

地下水运动、堤坝渗流一般是饱和2非饱和问题,建筑物地基和道路路基的固结、堤坝的变形都包含这一问题.因而,研究饱和2非饱和渗流问题是很有工程意义的.

1　非饱和非稳定渗流有限单元法

孔隙水非稳定渗流的微分方程为:

99x (k x (n ,s )9H 9x )+99y (k y (n ,s )9H 9y )+99z (k z (n ,s )5H 5z )=9(sn )

9t ,H (x ,y ,z ,t )=H 1(x ,y ,z ,t ),在s 1上,

k x

9H 9x cos ( n ,x )+k y 9H 9y cos ( n ,y )+k z 9H 9z

cos ( n ,z )=q ,在s 2上,H (x ,y ,z ,t )=z (x ,y ,z ,t ),在s 3上,H (x ,y ,z ,t 0)=H 0(x ,y ,z ,t 0).

(1)

式中H 为水头,H =u/r w +z ,u 为孔隙水压力,r w 为水的容重,z 为位置水头,s 为饱和度,直

角坐标轴x 、y 、z 为渗透主方向,k x 、k y 、k z 分别为沿主方向的渗透系数,s 1为已知水头边界,s 2为已知流量边界,s 3为逸出段边界,H 0为已知水头,q 为边界流量,n 为土体孔隙度,cos ( n ,x )等为边界面外法线方向的方向余弦,t 为时间.

饱和度是吸力(孔隙气压力与水压力之差)的函数,假设孔隙气压力瞬时消散,孔隙气压力可以看成一常量,此时有:

9(s ・n )9t =s 9n 9t +n 9s 9t =s 9n 9t +ns ′(u )9u

9t

,

(2)其中:

s ′(u )=9s

9u

.

(3)　　不考虑土骨架变形对渗流的影响,式(2)可表示为

9(s ・n )9t =ns ′(u )9u

9t ,

(4)则不考虑骨架变形的孔隙水渗流微分方程可表示为:

[99x (k x (s )9u 9x )+99y (k y (s )9u 9y )+99z (k z (s )9u 9z )]/r w +99z (k z (s ))=s ′(u )n 9u

9t

.

(5)

—

83—Ξ本文于1998年8月18日收到.

应用变分法,并对时间取隐式差分,可得式(5)的有限元方程:

[k s w ]e (1r w {u}e +{z })+µ

e

s ′(u )n Δt

[N ]T

[N ]d x d y d z{u}e -µ

e

s ′

(u )n Δt

[N ]T

[N ]d x d y d z{u 0}e =-µ

s 2

∩e q w [N ]T

d s.

(6)

其中:

[k s

w ]e =µ

e

[B ′]T

[S ][B ′]d x d y d z ,

(7)

[S ]=

k x (s )

00k y (s )

k z (s

)

,[B ′]=

N 1,x …N n ,

x N 1,y …N n ,y N 1,z

…N n ,z

,[N ]=[N 1　N 2　…　N n ],N i 为单元形函数.

在非饱和渗流数值分析中,当入渗处土体饱和度比较低时,常出现该处附近孔隙负压计算结果分

布规律较乱,与实际物理现象不符的现象,即所谓的数值弥散现象[6,7].为消除这一现象,文献[6]采用变坐标的特征有限元法,文献[7]采用在G alerkin 有限元法的权函数中加上一个扰动的SU P G 有限元法.这些方法十分复杂,不便于推广应用.本文认为产生这一现象的原因是有限元迭代计算过程中s ′(u )取值不当所造成,本文对式(6)作如下的处理,较好地解决了计算中的数值弥散问题.

由于土体孔隙水压力大于0时,孔压变化不引起饱和度变化,可设如下的函数:

β(u )=1,　u <0时

0.　u ≥0时.

(8)

不妨取式(3)为如下的形式:

s ′(u )=s (u )-s (u 0)

β(u )・u -β(u 0)・u 0,当

β(u )u -β(u 0)u 0≠0时,

s (u -0101)-s (u 0)

(β(u )・u -0101)-β(u 0)・u 0

,当β(u )u -β(u 0)u 0=0时.

(9)

则式(6)可改写为:

[k s

w ]e (

1

r w

{u}e +{z })+

µ

e

s ′

(u )n Δt

[N ]T

[N ]d x d y d z [β(u )]{u}e -µ

e

s ′

(u )n Δt

[N ]T

[N ]d x d y d z [β(u 0)e ]{u 0}e =-µ

s 2

∩e q w {N }T

d s

(10)

其中:

[β(u )]=

β(u 1)

β(u 2)

ω

β(u n )

,[β(u 0)]=

β(u 01)

β(u 02)

ω

β(u 0n )

.

式(10)即为不考虑土骨架变形对渗流影响的饱和-非饱和非稳定渗流有限元方程组.

2　有限元计算中几个问题的处理

渗流数值分析的有限元法中,逸出面边界的确定和非线性迭代的处理是计算的关键.下面介绍本文的处理方法.211　逸出面边界条件处理的新方法　渗流数值分析的有限单元法中,已知水头和已知流量边界,将边界条件代入有限元方程中即可,很容易处理,而对逸出面边界,由于其范围要在计算中确定,因而较难处理.过去常常按以下两种方法处理:(1)在计算数据中给出可能逸出带内的边界结点,并给出

其转换为逸出点的先后趋势,并在计算过程中判断后决定是否转换为逸出点,那些水头值等于和超出结点高程的边界点按逸出点对待,按第一类边界点处理[1,2,3].这需要根据地形地质条件和渗流知识给出这些数据.由于计算时不能实行逆转换,即不能将已定为逸出点的结点转换为非逸出点,从第一类边界点中消去,而三维问题的复杂性又增加了试算过程中预测的难度[3];(2)取域内近边界面浸润线上一点,取域内浸润线外延线与过该点平行于边界的向下射线(边界坡角小于等于90度)或铅垂线(坡角大于90度)的夹角的平分线与边界的交点定为出逸点,出逸点下的结点为逸出面结点[4],该法在三维问题中难于应用.

显然,逸出面边界的水是由域内流出边界的.根据这一原理,本文提出了如下的逸出面边界处理方法.

(1)第一次计算,对于初始时步,将可能的逸出面结点全部作为未知点,对于其余时步,将时步

初的逸出面结点作为已知点,经过有限元计算求出渗流场分布.

(2)将水面以上的正压边界点全部按逸出点处理,求出渗流场分布后,再对水面以上零压边界点进行判断,若流通量大于等于零,即:

-k x 9H 9x cos (n ,x )-k y 9H 9y cos (n ,y )-k z 9H 9z

cos (

n ,z )≥0,

(11)则该点为逸出面结点.对于等参单元:

9H 9x 9H 9y 9H

9z

=[B ′]e {H}e .

(12)若式(11)成立,说明该处水流向域外,则结点为逸出面结点,否则,该点不是逸出面结点.

(3)将前一次确定的逸出面结点作为已知点处理.重复第3步,直到计算收敛,进行下一时间步计算.

212　时步选择与非线性迭代　非稳定渗流,是随时间推移而发展的,因而计算需分时段进行,每一

时段分若干时步计算.采用合适的时步,既能保证迭代收敛,又可避免步长太小而引起过大的累积误差,并节省机时.每一时段的初始时步由用户输入,若迭代收敛,则取此值为该时段的迭代步长,反之,则步长减半,直到获取收敛步长为止.

非饱和土的渗流问题是一个非线性问题,饱和度、渗透系数、逸出面范围等只能在计算过程中迭代确定.迭代可分2层进行,内层迭代确定渗流参数,该层迭代中边界条件不变,渗流参数可采用中点饱和度值所对应的参数,直到前后两次计算各结点孔压变化均小于某一该定小值;外层迭代确定渗流逸出面边界,该层迭代根据内层的计算结果,调整逸出面边界范围,直至迭代收敛.由于非饱和土的饱和度和渗透系数是孔压的函数,其关系曲线可由试验测定[10],并可按离散点输入程序,计算中按线性插值取值.

3　计算实例

表1　砂土的土水特性关系

饱和度吸力/kPa 相对渗透系数k r

10101820132501970101650194015001910120014511001130118

1010

0101

　　为检验上述方法的正确性,本文用文献[8]中介绍的砂槽模型的试验资料、假想均质粘土坝进行了渗流分析.算例中土料各向同性,且不计变形对渗流的影响.

311　砂槽模型　文献[8]介绍的赤井浩一等作者的砂槽模型长315cm ,宽23cm ,高33cm.模型材料为均匀砂,孔隙率为0144,

饱和渗透系数为313×10-3m/s ,非饱和砂土的饱和度、吸力、相对渗透系数关系如表1示.初始条件采用上下游水位均为10cm 的平衡态,在时间0s ,上游水位骤升至30cm.计算域划分为300个1015cm ×313cm 的矩形单元.计算条件与模型试验相同,计算成果与试验结果比较于图1示.可以

看出,计算结果与试验成果比较接近.—

4—

图1　自由水面数值解与砂槽模型试验结果比较

表2　粘土的土水特性关系

饱和度

吸力/kPa

相对渗透系数k r

1010192211199180×10-101976331109110×10-1019396391444180×10-1016315471871100×10-1015280511316100×10-2013793611121100×10-20128667212100×10-3012125921022100×10-40117761151072100×10-5011547

138112

3100×10-6

312均质土坝　某假想均质粘土坝,坝底相对高程010m ,

坝高50m ,顶宽10m ,上游坝坡1∶2,下游坝坡1∶115,

坝体下游设水平褥垫式排水,伸入坝内6715m ,饱和渗透系数为314×10-8m/s ,孔隙率为014,土体水分特性参数如表2[9].初始条件假定为坝体竣工时土体的饱和度均为80%,上下游无水.上游水位30d 内均匀上升至4510m ,第5000d 上游水位开始下降,30d 内均匀降至510m ,下游水位均为010m ,计算库水上升和下降后坝体

的浸润线变化过程和孔压分布情况,坝体计算网格剖分如

图2示.图3为0至1500d 坝体浸润线变化过程,图4为第1500d 时的孔压分布,图5为第5000d 的孔压分布,图6为第5000至第10000d 的浸润线变化过程

.

图2　

土坝网格剖分图3　5000d

前浸润线变化过程

图4　第1500d 孔压等值线(单位:kPa

)图5　第5000d 孔压等值线(单位:kPa

)

图6　5000d

后浸润线变化过程图7　第10000d 孔压等值线(单位:kPa )

　　从以上算例的浸润线变化过程图可以看出,水位变化的初期,浸润线位置变化很快,随着时间的推移,浸润线位置变化逐步变慢,直至趋于稳定.以上算例的饱和与非饱和区的孔压分布和逸出点位

置都是合理的,逸出面位置一般迭代2至3次就收敛,说明本文提出的处理方法可行.

本文提出的非饱和土渗流有限元计算中的处理方法,很好地消除了数值弥散现象,边界处理方法简单有效,具有推广价值.

参考文献

1　陈平,李祖饴.堤坝非稳定渗流计算程序的编制及其应用.水利水运科学研究.1990,(4).

2　吴梦喜,张学勤.有自由面渗流数值分析的虚单元法.水利学报,1994,(8).

3　李春华.固定网络法用于三向稳定渗流有限元的实践.北京:水利电力出版社,第四界全国水利水电工程渗流学术讨论会论文集,1993.

4　河野伊一郎.用有限元法解坝体的渗流问题.水利水运科技情报,1978年,增刊3.

5　吴良骥,G.L.Bloomaburg.饱和非饱和区中渗流问题的数值模拟.水利水运科学研究,1985,(6).

6　黄康乐.求解非饱和土壤水流问题的一种数值方法.水利学报,1987,(9).

7　朱学愚,谢春红,钱孝星.非饱和流动问题的SU PG有限元法.水利学报,1994,(6).

8　高骥,雷光耀,张锁春.堤坝饱和2非饱和渗流的数值分析.岩土工程学报,1988,(6).

9　杨代泉,沈珠江.非饱和土一维广义固结的数值计算.水利水运科学研究,1991,(4).

10　Fredlund D G,Rahardio H.非饱和土土力学.陈仲颐等译,北京:中国建筑工业出版社,1997.

S aturated2unsaturated unsteady seepage numerical analysis

Wu Mengxi G ao Lianshi

(Tsi nghua U niversity)

Abstract　This paper presents an adapted finite element method for unsaturated analysis.Numer2 ical confusion phenomenon in unsaturated seepage numerical analysis is avoided effectively.A new method for seepage spill boundary handling is also presented.Some analysis examples are given.

K ey w ords　unsaturated,unsteady,seepage,spill boundary,numerical method.

—

—

2

4