三维非稳定饱和-非饱和渗流有限元法及其改进

非饱和渗流场的变化过程。由图可见，等势线均光滑平顺，各时刻位势场连续变化，没有出
现振荡等现象。此外，由于孔口压力水头较低、孔口面积很小，因此孔口处的饱和区很小，
且试验开始时饱和区迅速增大，过一定时间后维持稳定。随着时间的增长，负压区的负压力
逐渐减少，在底部出现饱和区，饱和区逐渐升高，直至趋向稳定，如图2（h）所示。孔口部
3. 计算方法的改进
文献[ 7 ]中指出计算 ∂θ 时，可能出现三种情况：（1）两相邻时步间均为非饱和，且在 ∂h
非饱和吸力与含水率的关系曲线上不为同一点；（2）前一时步为饱和，后一时步为非饱和；
或前一时步为非饱和，后一时步为饱和；（3）在时间步之间均为非饱和，且在非饱和吸力与
含水率的关系曲线上为同一点。对于（1）、（2）两种情况，采用 ∂θ 为 θ (h) −θ (h0 ) ； ∂h ζ (h) ⋅ h − ζ (h0 ) ⋅ h0
∂h
[ B ] 、 [ P ] 。当计算满足收敛准则时，停止本时间步的计算。
2.3 时间步选择
非稳定渗流是随时间推移而发展的，因而计算需要分时段进行。时间步长的选取至关重 要，它直接影响到计算精度和计算工作量。如果时间步长取值太大，会使解的误差过大，计 算结果偏离真解；时间步长取值过小，计算工作量将成倍增大，不仅浪费机时，而且会导致 数值计算的误差增大。实际工程中不同研究对象的渗透性相差很大，可根据模型大小、渗透 性强弱来确定可行时间步的数量级(秒、时、天)，然后选定一初始时间步进行计算，如收敛， 则取此值为该时段的步长，否则，时间步减半，直到获取收敛的时间步为止。
ζ
(h0
)
=
⎧1 ⎨⎩0
h0 h0
< ≥
0 ；θ (h)
0
=
⎧⎪1
⎨ ⎪⎩
f
h
(h)
≥
0 h
<
0
；θ (h0 )
=
⎧⎪1 ⎨ ⎪⎩ f
h0
(h0 )
≥
0 h0
<
0
其中， f (h) 、 f (h0 ) 表示由θ − h 离散点线性插值。在计算中，根据实际工程情况设
置介质非饱和区域的负压水头的最大值 hL ，计算中如果有负压水头超过该值，则取 h = hL 。
续性，在计算中产生振荡，造成非饱和区孔隙负压计算结果分布规律较乱、与实际物理现象 不符的现象，即所谓的数值弥散现象[ 3] [ 4 ]。为消除这一现象，黄康乐[ 3]采用变坐标的特征 有限元法，朱学愚等[ 4 ]采用在Galerkin有限元法的权函数中加上一个扰动的SUPG有限元法， 这些方法较好地改善了数值弥散现象，但是均较复杂。吴梦喜[ 5 ]提出了类似式（10）函数
饱和－非饱和渗流在计算方法上把饱和－非饱和区作为一个整个流场研究，不需要进行 自由面的调整，自由面仅是负压区正压区的分界面，不再是像饱和渗流场计算时流场的一个 边界，避免了饱和渗流计算中试找自由面和选取给水度补给自由面边界的麻烦。但是，通常
k − h 和θ − h 的函数关系是由离散的点线性插值得来，θ − h 的线性插值易造成 ∂θ 的非连 ∂h
ζ m 为单元结点的局部坐标；ξ ，η ，ζ 为单元内任一点的局部作标。
由Galerkin加权余量法及格林第一公式有
[
A]{hc
}
+
[
B]
⎧ ⎨ ⎩
∂hc ∂t
⎫ ⎬ ⎭
=
[
P
]
（8）
对时间采用隐式有限差分格式，即 ∂hc ∂t
=
1 ∆t
( ) ⎡⎣ hc t+∆t
−
(
hc
) t
⎤⎦
，代入式（8）可得
kr = 1
； C = ∂θ 为容水度，θ 为介质含水率，在正压力区 C
∂h
=0
； β 为饱和－非饱和选
择常数，在非饱和区为0，在饱和区为1； Ss 为弹性贮水率，饱和土体的 Ss 为常数，在非饱
和土体中 Ss = 0，当忽略土体骨架及水的压缩性时对于饱和区也有 Ss = 0 ； Q 为源汇项； i, j = 3 为垂直坐标； t 为时间变量。
ks = 0.00077 cm/s；体积含水率θ 与负压 h 的关系，以及体积含水率与相对透水率 kr 的关
系如表1所示。考虑到试验装置的轴对称性，取以圆筒中心（孔口）为中心轴、周向10o夹角 的范围建立有限元模型，模型网格图如图1所示。径向剖面各时段的等势线图如2所示。
表 1 中细砂的土水特性关系
5
图 1 模型网格图
（a）t=600s
（b）t=1800s
-4-
http://www.paper.edu.cn
（c）t=3000s
（d）t=4200s
（e）t=4800s
（f）t=5400s
（g）t=6000s
（h）稳定期
图 2 各时段的等势线图（单位水头：m）
上述径向剖面各时刻的等势线图清晰地表现了复合土工膜破损孔口渗漏时模型内饱和
h /c -200.0 -100.0 -70.0 -50.0 -40.0 -35.0 -30.0 -25.0 -15.0 -13.0 0
m
0
0
0
0
0
0
0
0
0
0
kr
0 0.0001 0.001 0.01 0.03 0.082 0.225 0.55 0.886 0.963 1
θ
0
0.062 0.100 0.150 0.175 0.201 0.225 0.255 0.288 0.300 0.3
从而可求出单元各高斯点的压力水头，对于压力水头大于或等于零的高斯点，处于饱和区， 渗透系数取饱和值，对于压力水头小于零的高斯点，处于非饱和区，渗透系数取值可由渗透
介质 k − h 曲线得出， ∂θ 也可由θ − h 曲线得出。在计算单元矩阵[ A]e 、[B]e 、[P]e 时，
∂h
由上述方法确定各高斯点上的 k 及 ∂θ ，求出[ A]e 、[B]e 、[P]e ，再将各单元迭加，得[ A]、
⎡⎣C
(
hc
)
+
β
Ss
⎤⎦
d
Ω⎤⎦
；
e=1
∑ ∫∫∫ ∑ ∑ ∫∫∫ w∫∫ {P}
=
−
NE
⎡ ⎢
e=1 ⎣
( ) 3
k Ωe i=1 r
hc
kis3
∂N ∂xi
⎤ dΩ⎥
⎦
−
NE e=1
⎡ ⎣
Ωe NnQdΩ +
2
qn
Nn
ds
⎤ ⎦
； NE 为单元
总数；Nn
，N
m
为单元结点
n
，m
的形函数；t
为时间变量，∆t
的渗透系数 k 及 ∂θ ，进而修正矩阵[ A]、[B]、[P]，这里以高斯点为研究对象，在单元
∂h 中不同的高斯点有不同的压力水头，进而可得到不同的渗透系数 k 及 ∂θ ，这样修正得的矩
∂h
阵[ A]、[B]、[P]比较合理。在时间 t + ∆t 中 i 次迭代计算后，可得各节点的压力水头 h ，
算例中详细模拟了复合土工膜破损渗漏时的试验模型内饱和非饱和渗流场的变化过程， 计算时数值稳定，计算成果与试验结果基本吻合。该非稳定饱和―非饱和有限元程序是正确 的。
-5-
http://www.paper.edu.cn
2. 有限元法
2.1 基本方程
采用固定网格法对整个渗流计算区域进行有限元离散，单元内压力水头在空间和时间域
上的变化近似地表达为
hc = Nmhcm
Nm
=
1 8
(1+ ξξm )(1+ηηm )(1+ ζζ m )
m = 1, 2",8
（7）
式中：m 单元结点数，hcm 为单元结点 m 的压力水头，Nm 为单元结点 m 的形函数；ξm ，ηm ，
同时，在饱和与非饱和接触区域，网格大小不亦超过 hL 的一半。
4. 算例
为了检验上述方法的正确性和有效性，这里考虑某复合土工膜缺陷渗漏量室内物理模型 试验，建立数学模型对其进行计算，模拟土工膜破损后的渗流发展过程，计算渗流稳定后的 渗漏量。
此试验是在直径1m、高1m的圆桶内进行，圆桶内装满干的中细砂，圆桶材料不透水， 桶壁底部钻有小孔供排水，顶部覆盖复合土工膜，复合土工膜上中心钻小孔(2mm)，在整个 试验过程中小孔上维持1m的压力水头。中细砂饱和状态的各向同性渗透系数
位坡降较大，超过试验土料的允许渗透坡降，可能出现渗透破坏。这与试验结束时检查所见
垫层材料被渗漏水冲刷出一个“冲刷坑”的现象相吻合。稳定时渗漏量为1.29 ml s ，与试验值
1.08 ml s 基本一致。可见，本计算方法和程序是正确的。
5. 结语
对非稳定饱和―非饱和有限元分析方法加以改进，提出了容水度的改进算法，有效地消 除了计算中出现的数值弥散现象。
∂ ∂xi
⎡ ⎢⎢⎣kisj kr
(hc )
∂hc ∂x j
+ kis3kr
(
hc
)⎤⎥
⎥⎦
−Q
=
⎡⎣C (hc ) +
β Ss ⎤⎦
∂hc ∂t
（1）
式中： h 为总水头； hc 为压力水头； kisj 为饱和渗透系数张量； kr 为相对透水率，是非
饱和土的渗透系数与同一种土饱和时的渗透系数的比值，在非饱和区0 < kr < 1，在饱和区
来计算容水度 C ，能较好的解决这一问题，但表达式不严密，没有分析表达式中分母等于 零时的情况。这里在文献[ 5 ]中对计算容水度 C 的处理方法加以改进。
1. 非稳定饱和－非饱和土渗流定解问题
对于服从达西定律的非均质各向异性饱和－非饱和土，在忽略水的压缩性和孔隙气体对 水流运动的影响情况下渗流基本微分方程为[ 1]：
http://www.paper.edu.cn
三维非稳定饱和－非饱和渗流有限元法及其改进
江沆，沈振中Βιβλιοθήκη Baidu邱乾勇
河海大学水利水电工程学院，江苏南京（210098）
E-mail：jw01021514@hhu.edu.cn。
摘 要：阐述了非稳定饱和－非饱和渗流的基本方程、定解条件和有限元迭代计算过程。针 对有限元计算过程中经常出现的数值弥散现象，研究了其出现的规律和处理方法，提出了容 水度的改进算法，并编制了有限元程序。最后详细模拟复合土工膜缺陷渗漏量室内物理模型 试验，给出了试验过程中饱和－非饱和渗流场的变化过程。计算表明改进后的计算程序能有 效地消除数值弥散现象，计算成果与试验结果基本一致。 关键词：非稳定饱和－非饱和渗流；容水度；有限元法
为时间增量；(
hc
) t
，(
hc
) t
+ ∆t
分别为时间 t ， t + ∆t 时的节点压力水头； Ω 为计算空间域； s 为计算域边界。
式（9）即为非稳定饱和－非饱和渗流有限元计算的基本方程。
-2-
http://www.paper.edu.cn
2.2 迭代解法
在非稳定渗流计算的每个时间步中须进行迭代，迭代的过程就是根据压力水头确定单元
式(2) ～ (6) 中： ni 为边界面外法线方向余弦；t0 为初始时刻； hc1 为已知水头； qn 为
已知流量； qr (t ) 为降雨入渗流量； hc (t0 ) 为初始 t0 时刻渗流场水头； Γ1 为已知结点水头
边界； Γ2 为已知流量边界； Γ3 为降雨入渗边界； Γ4 为饱和逸出面边界。
hc
⎤
⎥ ⎥⎦
ni
Γ2
= qn
（4）
( ) ( ) ⎡
− ⎢kisj kr ⎢⎣
hc
∂hc ∂x j
+ kis3kr
hc
⎤ ⎥ ni Γ3 ≥ 0且hc Γ3 = 0 ⎥⎦
（5）
−
⎡ ⎢kisj ⎢⎣
kr
(
hc
)
∂hc ∂x j
( )⎤
+ kis3kr hc ⎥ ni ⎥⎦
Γ4
= qr (t )
（6）
⎛⎜⎝ [
A]
+
1 ∆t
[
B
]
⎞ ⎟⎠
{hc
} t
+
∆t
=
{P}
+
1 ∆t
[
B
]{hc
} t
（9）
其
中
：
[ ] ∑ ∫∫∫ ∑ ∑ NE ⎡
A= ⎢ e=1 ⎢⎣
( ) 3 3
k Ωe i=1 j=1 r
hc
kisj
∂Nn ∂xi
∂Nm ∂x j
⎤ dΩ⎥
⎥⎦
；
NE
[
B]
=
∑
⎡ ⎣
∫∫∫Ωe
Nn
Nm
初始条件：
hc ( xi , 0) = hc ( xi ,t0 ) ， i = 1, 2,3
（2）
边界条件：
hc ( xi , t ) ( ) Γ1 = hc1 xi , t
（3）
-1-
http://www.paper.edu.cn
( ) ( ) ⎡
−
⎢ ⎢⎣
kisj
kr
hc
∂hc ∂x j
+ kis3kr
⎪⎩ζ (h) ⋅ (h − 0.01) − ζ (h0 ) ⋅ h0
ζ (h) ⋅ h − ζ (h0 ) ⋅ h0 ≠ 0 ζ (h) ⋅ h − ζ (h0 ) ⋅ h0 = 0
（10）
式中：h
为本次时间步的压力水头，h0
为上一时间步的压力水头；ζ
(h)
=
⎧1 ⎨⎩0
h<0
；
h≥0
-3-
http://www.paper.edu.cn
当出现（3）情况时，认为 ∂θ 为零。但从非饱和吸力与含水率的关系曲线来看，但前后两 ∂h
时步为同一点时，这点的斜率不一定为零。因此，这里采用式（11）来解决这一问题。
⎧ θ (h) −θ (h0 )
∂θ ∂h
=
⎪⎪ζ ⎨ ⎪
(h) ⋅ h − ζ (h0 ) ⋅ h0 θ (h − 0.01) −θ (h0 )