饱和非饱和土体非稳定渗流数值分析_吴梦喜

程中 s′(u)取值不当所造成 , 本文对式 (6)作如下的处理 , 较好地解决了计算中的数值弥散问题 .
由于土体孔隙水压力大于 0 时 , 孔压变化不引起饱和度变化 , 可设如下的函数 :
β(u )=
1 , u 0 . u
<0 时 ≥0 时 .
不妨取式 (3)为如下的形式 :
(8)
s′(u )=
饱和度是吸力 (孔隙气压力与水压力之差)的函数 , 假设孔隙气压力瞬时消散 , 孔隙气压力可以
看成一常量 , 此时有 :
(s
· t
n) =
s
n t
+n
s t
=s
n t
+ n s′(u )
u t
,
(2)
其中 :
s′(u ) =
s u
.
不考虑土骨架变形对渗流的影响 , 式 (2)可表示为
(s
· t
n)=
图 4 第 1500d 孔压等值线 (单位 :kPa)
图 5 第 5000d 孔压等值线 (单位 :kPa)
图 6 5000d 后浸润线变化过程
ns′(u)
u t
,
则不考虑骨架变形的孔隙水渗流微分方程可表示为 :
(3) (4)
[
x(k x(s)
xu )+
y(ky (s)
yu )+
z (kz (s)
u z
)] /
r
w
+
z (kz (s))= s′(u)n
u t
.
(5)
本文于 1998 年 8 月 18 日收到 .
— 38 —
应用变分法 , 并对时间取隐式差分 , 可得式 (5)的有限元方程 :
布规律较乱 , 与实际物理现象不符的现象 , 即所谓的数值弥散现象[ 6, 7] .为消除这一现象 , 文献 [ 6]
采用变坐标的特征有限元法 , 文献 [ 7] 采用在 Galerkin 有限元法的权函数中加上一个扰动的 SUPG
有限元法 .这些方法十分复杂 , 不便于推广应用 .本文认为产生这一现象的原因是有限元迭代计算过
相对渗透系数 kr
1 9.80 ×10 -1 9.10 ×10 -1 4.80 ×10 -1 1.00 ×10 -1 6.00 ×10 -2 1.00 ×10 -2 2.00 ×10 -3 2.00 ×10 -4 2.00 ×10 -5 3.00 ×10 -6
图 2 土坝网格剖分
图 3 5000d 前浸润线变化过程
0.01
平衡态 , 在时间 0s , 上游水位骤升至 30cm .计算域划分为 300 个
10.5cm ×3.3cm 的矩形单元 .计算条件与模型试验相同 , 计算成果与试验结果比较于图 1 示 .可以
看出 , 计算结果与试验成果比较接近 .
— 40 —
图 1 自由水面数值解与砂槽模型试验结果比较
DO I :10.13243/j .cnki .slxb.1999.12.007
1999 年 12 月
水 利 学 报
SH UIL I XUEBAO
第 12 期
饱和-非饱和土体非稳定渗流数值分析
吴梦喜 高莲士
(清华大学水电系 , 北京 , 100084)
摘 要 本文对一般的非饱和渗流有限元计算方法 加以改 进 , 有效 地消除了 非饱和 渗流数 值计算 中存在 的数值 弥散现象 .同时还提出了一种简便有效的逸出面处 理新方法 , 并给出了非饱和非稳定渗流计算的实例 . 关键词 非饱和 , 非稳定 , 渗流 , 逸出面 , 数值方法 . 中图号 T V 139.14
— 39 —
其转换为逸出点的先后趋势 , 并在计算过程中判断后决定是否转换为逸出点 , 那些水头值等于和超出 结点高程的边界点按逸出点对待 , 按第一类边界点处理[ 1 , 2 , 3] .这需要根据地形地质条件和渗流知识
给出这些数据 .由于计算时不能实行逆转换 , 即不能将已定为逸出点的结点转换为非逸出点 , 从第一 类边界点中消去 , 而三维问题的复杂性又增加了试算过程中预测的难度[ 3] ;(2)取域内近边界面浸 润线上一点 , 取域内浸润线外延线与过该点平行于边界的向下射线 (边界坡角小于等于 90 度)或铅 垂线 (坡角大于 90 度)的夹角的平分线与边界的交点定为出逸点 , 出逸点下的结点为逸出面结点[ 4] , 该法在三维问题中难于应用 .
H(x , y , z , t 0)= H 0(x , y , z , t 0).
式中 H 为水头 , H =u/ rw +z , u 为孔隙水压力 , rw 为水的容重 , z 为位置水头 , s 为饱和度 , 直 角坐标轴 x 、 y 、 z 为渗透主方向 , kx 、 ky 、 kz 分别为沿主方向的渗透系数 , s1 为已知水头边界 , s2 为已知流量边界 , s3 为逸出段边界 , H0 为已知水头 , q 为边界流量 , n 为土体孔隙度 , cos (n , x ) 等为边界面外法线方向的方向余弦 , t 为时间 .
2 有限元计算中几个问题的处理
渗流数值分析的有限元法中 , 逸出面边界的确定和非线性迭代的处理是计算的关键 .下面介绍本 文的处理方法 . 2.1 逸出面边界条件处理的新方法 渗流数值分析的有限单元法中 , 已知水头和已知流量边界 , 将 边界条件代入有限元方程中即可 , 很容易处理 , 而对逸出面边界 , 由于其范围要在计算中确定 , 因而 较难处理 .过去常常按以下两种方法处理 :(1)在计算数据中给出可能逸出带内的边界结点 , 并给出
3.2 均质土坝 某假想均质粘土坝 , 坝底相对高程 0.0m , 坝高 50m , 顶宽 10m , 上游坝坡 1∶2 , 下 游坝坡 1∶1.5 , 坝体下游设水平褥垫式排水 , 伸入坝内 67.5m , 饱和渗透 系数为 3.4 ×10-8m/ s , 孔 隙率为 0.464 , 土体水分特性 参数如表 2[ 9] .初始条件假定为坝体竣工时土体的饱和度 均为 80 %, 上下 游 无水 .上游 水 位 30d 内 均匀 上 升至 45.0m , 第 5000d 上游水位 开始下降 , 30 d 内均 匀降至 5.0m , 下游水位均为 0.0m , 计算库水上升和下降后坝体 的浸润线变化过程和孔压分布情况 , 坝体计算网格剖分如 图 2 示 .图 3 为 0 至 1500d 坝体浸润线变化过程 , 图 4 为 第 1500d 时的孔压分布 , 图 5 为第 5000d 的孔压分布 , 图 6 为第 5000 至第 10000d 的浸润线变化过程 .
β(u
s(u )-s(u )· u -β(u
0) 0)·
u
0
,
当
β(u )u
-β(u 0)u0
≠0 时,
s(u -0.01)-s(u0) (β(u)· u -0.01)-β(u0)· u0
,当
β(u )u
- β(u 0 )u 0
=0时 .
(9)
则式 (6)可改写为 :
[ ksw] e(r1w {u}e +{z})+ e s′(u)Δnt [ N ] T[ N ] d x dy dz[ β(u)] {u}e -
表 2 粘土的土水特性关系
饱和度
1 0.9892 0.9763 0.9396 0.6315 0.5280 0.3793 0.2866 0.2125 0.1776 0.1547
吸力/ kPa
0 21.19 31.0 39.44 47.87 51.31 61.12 72.89 92.02 115.07 138.12
[ B′] T [ S ] [ B′] d x dy dz ,
e
0
N 1, x …
Nn, x
(7)
[ S ] = 0 ky(s) 0 ,[ B′] = N 1 , y … N n , y ,
0
0 k z(s)
N1, z … Nn, z
[ N ] = [ N 1 N 2 … N n] , N i 为单元形函数 . 在非饱和渗流数值分析中 , 当入渗处土体饱和度比较低时 , 常出现该处附近孔隙负压计算结果分
-k x
H x
cos(n
,
x)-ky
H y
cos(n
ห้องสมุดไป่ตู้
,
y)-kz
H z
cos(n
,
z)≥0
,
则该点为逸出面结点 .对于等参单元 :
(11)
H x
H y
H z
=[ B′] e{H}e .
(12)
若式 (11)成立 , 说明该处水流向域外 , 则结点为逸出面结点 , 否则 , 该点不是逸出面结点 .
(3) 将前一次确定的逸出面结点作为已知点处理 .重复第 3 步 , 直到计算收敛 , 进行下一时间步
计算 .
2.2 时步选择与非线性迭代 非稳定渗流 , 是随时间推移而发展的 , 因而计算需分时段进行 , 每一
时段分若干时步计算 .采用合适的时步 , 既能保证迭代收敛 , 又可避免步长太小而引起过大的累积误
差 , 并节省机时 .每一时段的初始时步由用户输入 , 若迭代收敛 , 则取此值为该时段的迭代步长 , 反
显然 , 逸出面边界的水是由域内流出边界的 .根据这一原理 , 本文提出了如下的逸出面边界处理 方法 .
(1) 第一次计算 , 对于初始时步 , 将可能的逸出面结点全部作为未知点 , 对于其余时步 , 将时步 初的逸出面结点作为已知点 , 经过有限元计算求出渗流场分布 .
(2) 将水面以上的正压边界点全部按逸出点处理 , 求出渗流场分布后 , 再对水面以上零压边界点 进行判断 , 若流通量大于等于零 , 即 :
之 , 则步长减半 , 直到获取收敛步长为止 .
非饱和土的渗流问题是一个非线性问题 , 饱和度 、 渗透系数 、 逸出面范围等只能在计算过程中迭
代确定 .迭代可分 2 层进行 , 内层迭代确定渗流参数 , 该层迭代中边界条件不变 , 渗流参数可采用中
点饱和度值所对应的参数 , 直到前后两次计算各结点孔压变化均小于某一该定小值 ;外层迭代确定渗
表 1 砂土的土水特性关系
型的试验资料 、 假想均质粘土坝进行了渗流分析 .算例中土料各 饱和度 吸力/ kPa 相对渗透系数 kr
向同性 , 且不计变形对渗流的影响 .
1
0
1
3.1 砂槽模型 文献 [ 8] 介绍的赤井浩一等作者的砂槽模型长
0.82 0.64
0.325 0.65
0.97 0.94
其中 :
e s′(u)Δnt [ N] T[ N ] dx dy dz [ β(u0)e] {u 0}e =- s2 ∩e qw{N }T ds
(10)
β(u 1)
β(u 0 1)
[ β(u)] =
β(u 2)
, [ β(u 0)] =
β(u 0 2)
.
β(u n )
β(u0 n)
式 (10) 即为不考虑土骨架变形对渗流影响的饱和 -非饱和非稳定渗流有限元方程组 .
地下水运动 、 堤坝渗流一般是饱和-非饱和问题 , 建筑物地基和道路路基的固结 、 堤坝的变形都 包含这一问题 .因而 , 研究饱和-非饱和渗流问题是很有工程意义的 .
1 非饱和非稳定渗流有限单元法
孔隙水非稳定渗流的微分方程为 :
x (kx(n , s)
Hx )+
y (k y(n
, s)
H y
)+
z(kz (n , s)
H z
)=
(sn t
),
H(x , y , z , t )= H 1(x , y , z , t), 在 s1 上 ,
kx
H x
cos(n
,
x)+ky
H y
cos(n
,
y)+kz
H z
cos(n
,
z)=q ,
在
s2
上
,
(1)
H(x , y , z , t )= z(x , y , z , t), 在 s3 上 ,
流逸出面边界 , 该层迭代根据内层的计算结果 , 调整逸出面边界范围 , 直至迭代收敛 .由于非饱和土 的饱和度和渗透系数是孔压的函数 , 其关系曲线可由试验测定[ 10] , 并可按离散点输入程序 , 计算中
按线性插值取值 .
3 计算实例
为检验上述方法的正确性 , 本文用文献 [ 8] 中介绍的砂槽模
[ ksw] e(r1w {u}e +{z })+ e s′(u)Δnt [ N ] T[ N ] d xdy dz {u}e -
其中 :
e s′(u)Δnt [ N ] T [ N ] d x dy dz{u 0}e =- s ∩ e qw[ N ] T ds . 2
(6)
[ ksw] e = k x(s) 0
315cm , 宽 23cm , 高 33cm .模型材料为均匀砂 , 孔隙率为 0.44 ,
0.50
0.91
0.20
饱和渗透系数为 3.3 ×10-3 m/ s , 非饱和砂土的饱和度 、 吸力 、 相
0.45
1.0
0.13
对渗透系数关系如表 1 示 .初始条件采用上下游水位均为 10cm 的
0.18 10.0