闸基渗透破坏的机理及数值模拟

∂h ∂h ∂ s [kij kr (hc ) c + kis3 kr (hc )] + S = [C (hc ) + β S S ] c ∂xi ∂x j ∂t
s
（2）
式中： kij 为饱和渗透系数张量； k r 为相对透水率，为非饱和土的渗透系数与同一种土饱和 时的渗透系数的比值，在非饱和区 0 < kr < 1 ，在饱和区 kr = 1 ； C 为容水度， C =

闸基渗透破坏机理及数值模拟
苏敏
河海大学水利水电工程学院，南京 (210098)
E-mail：loveguo0725@
摘 要： 本文主要阐明了闸基渗透破坏的机理， 以及编制了判断闸基是否发生渗透变形破坏 的计算程序。 关键词：闸基，非稳定渗流，渗透坡降， 渗透变形，有限元法
i=1，2，…8
（9）
通过解方程组（8） ，得出闸基非稳定渗流场的水头分布值。我们知道 J = − 在三个方向的渗透坡降为
-3-
dh ，那么单元 ds

Jx = −
则
dh dh dh , J y = − , Jz = − dx dy dz
2 2 2
(10) (11)
⎧ ∂h ⎫ ⎪ ∂x ⎪ ⎧Jx ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ∂h ⎪ e ⎨ J y ⎬ = − ⎨ ⎬ = −[ B]{h} ∂ y ⎪ ⎪ ⎪ ⎪ ⎩Jz ⎭ ⎪ ∂h ⎪ ⎪ ⎪ ⎩ ∂z ⎭
⎧ ∂N i ⎫ ⎧ ∂N i ⎫ ⎪ ∂ξ ⎪ ⎪ ∂x ⎪ ⎪ ⎪ ⎪ ⎪ −1 ⎪ ∂N i ⎪ ⎪ ∂N i ⎪ 其中 [ B ] = ⎨ ⎬ = [J ] ⎨ ⎬ ⎪ ∂y ⎪ ⎪ ∂η ⎪ ⎪ ∂N i ⎪ ⎪ ∂N i ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂z ⎭ ⎩ ∂ζ ⎭
−4 −4
/ m 。初始时刻上游水位 3.0
米，下游无水，然后上游水位骤升到 5.0 米。 （单元个数 2178，结点个数 3708，网格图见图
图2
试验样砂颗分曲线
由图 2 知， d85 = 0.33mm ， d15 = 0.15mm ， d f = 1.3 d85 d15 = 0.29mm ，则 p f =78 ﹪＞35﹪，是流土型的土，会发生流土破坏。通过查规范 SL265-2001[8],表 6.0.4 水平段和出 口段允许渗流坡降值，当渗流出口没有设置滤层时， J 允许 = 0.35 。
5.2 算例二：某闸基，除了上游水位骤升至 6.0 米，其他情况与算例一均相同。
-6-

2.3 渗流控制的原则
渗流控制包括控制理论和技术两个方面[4]。 控制理论是指基本原理和方法， 包括控制技 术的合理使用原则和优化设计。 控制技术是指基本原理的实施措施， 如灌浆、 防渗墙等技术。 渗流控制主要包括渗流量及安全控制两方面，安全多指渗透稳定。 众所周知，1910 年布莱总结了修建在不同地基上的 200 多座水工建筑物的运行经验， 首次发表了闸坝地基渗流控制法则，并称为爬行比理论。 布莱法则问世后， 1922 年 H. H 巴甫洛夫斯基通过渗流场的理论分析， 阐明了渗流场的 水力比降以进口、 出口为最大， 因此提出渗流控制应以控制渗流出口的水力比降不超过土壤 允许水力比降值为准则，并建议以此确定水工建筑物不透水轮廓的长度，即 J出口 ≤ J 允许 ， 式中 J出口 为地下轮廓渗流出口的水力比降； J 允许 为地基实际的允许水力比降。
1. 前言
水闸建成挡水后， 在上下游水位差的作用下将在闸基及两岸土体内产生渗流。 闸基渗流 属于有压渗流，而绕过两岸连接结构的渗流具有自由水面(浸润面)，属于无压渗流。闸基渗 流对闸室底板产生向上的渗透压力，减小了闸室的有效重量，对闸室的抗滑稳定不利。两岸 的绕渗，不仅对水闸岸、翼墙底面产生向上的渗透压力，而且会对墙背产生侧向压力，影响 岸、翼墙的稳定和强度安全。当渗透坡降和渗透流速超过某一限度时，会引起闸基及两岸土 的渗透变形，若渗透变形不能终止而继续发展，将导致水闸塌陷破坏[1]。本文基于闸基渗透 破坏的基本原理，研发了闸基是否会发生管涌或流土特性的渗透变形破坏的数值计算程序， 通过几个算例，也验证了程序的正确性。
2. 闸基渗透破坏的基本理论
2.1 渗透破坏的机理
通过工程实践经验和从室内模型试验观察到[2]，土体和水工建筑物的渗透破坏总是开始 于渗流出口，然而向深部发展。其原因不仅是渗流出口的水力比降最大，而且是表面处于无 应力状态。另外，如果土层内部有可移动的细颗粒而无出口，最终将自行停止，不会给建筑 物带来危害。工程破坏实例还表明，渗流出口的概念不仅是指下游河床表面，当渗流由弱透 水层进入无压的强透水层时， 强透水层将是弱透水层的渗流出口。 因为在渗流进入强透水层 后，渗透压力可能全部消失，在弱透水层中出现较大的水力比降。若强透水层的孔隙过大， 渗流就会将细颗粒带进强透水层， 使弱透水层遭破坏。 工程实际中有些防渗铺盖的渗透破坏 就是由于铺盖的土颗粒进入强透水的砂砾石层所造成的， 此时强透水层成为防渗铺盖的渗流 出口。
3.3 有限元法求解
对整个渗流计算域进行单元离散，单元内压力水头插值公式为（采用八节点） ：
h = ∑ N i (ξ ,η , ζ )hi
i =1
8
（5）
由 Galerkin 加权余量法可得求解的积分方程：
∫ { ∂x
e
∂
s k r ( hc ) [ k ij i
∂ hc ∂h + k is3 k r ( h c )] − C ( h c ) c } N n d v − ∂x j ∂t
开始 输入原始数据及信息 信息处理 形成带宽 形成总渗透矩阵及自由项 计 算 50 步
解线性代数方程组
计算每个单元形心 处的渗透坡降 < 否 去掉此单元，形成新的 单元信息及边界条件
输出新单元各点水头值
停止
图1
闸基非稳定渗流破坏计算程序框图
单元的渗透坡降取其形心处的渗透坡降，表达式：
-4-

s −[kij kr (hc )
∂hc + kis3 kr (hc )]ni |Γ2 = qn ∂x j ∂hc + kis3 kr (hc )]ni |Γ3 ≥ 0且hc |Γ3 = 0 ∂x j
（4）
s −[kij kr (hc )
式中： qn 为已知边界流量； ni 为外法线方向余弦；t0 为初始时刻；Γ1 为已知水头边界；Γ 2 为流量边界；为饱和逸出面边界， Γ 3 这里将非饱和逸出面边界简单地视为不透水面边界。
2.2 渗透破坏的形式
根据多年的研究,从渗透破坏形成机理的角度可以将渗透破坏形式分为流土、管涌两种 形式[3]。 流土,在渗透力作用下,土体中的某一颗粒群同时起动而流失的现象称为流土。发生流土 破坏的土又叫流土型土,这种破坏形式在粘性土和无粘性土中均可以发生。发生流土破坏的 土体,其中的土颗粒之间都是相互紧密结合的,相互之间具有较强的约束力,可以承受的水头 较大。但是,流土破坏的危害性却是最大的。因为,流土破坏一旦发生,它是土体的整体破坏, 流土通道会迅速向上游或横向延伸,一旦抢险不及时或措施不得当,就有造成土体结构破坏, 引发溃堤灾难发生的危险。 管涌,在渗透力的作用下,土体中的细颗粒(填料颗粒) 沿着土体骨架颗粒间的孔道移动 或被带出土体,这种现象叫管涌。发生管涌破坏的土又叫管涌型土,管涌的发生表明土体内有
3. 闸基渗流场有限单元法求解
3.1 数学模型
饱和-非饱和渗流连续性方程为[5]：
−
∂ ∂ ( ρ vi ) + S = ( ρ nS w ) ∂xi ∂t
（1）
式中： ρ ——水的密度； vi ——达西流速； n ——饱和含水率； S ——源汇项； S w ——土 体饱和度； xi ——直角坐标，i=1,2,3。 由达西定律 ， 以及考虑到体积含水率与压力水头之间存在函数关系 ， 可得多孔介质不 可压缩时饱和-非饱和、稳定-非稳定渗流场求解的微分方程为：
∑∫
ne
C (hc ) N m N n dv ∆t e =1 e
e =1 e ne ∂N n dv − ∑ ∫ qn N n ds ∂xi e =1 Γe 2
[ F ] = −∑ ∫ kis3 kr (hc )
得渗流场求解的有限单元法迭代方程为：
[ A]{hc } + [ S ]{
∂hc } = [F ] ∂t hc
∂ hc + k is3 k r ( h c )] n i + q n } N n d s = 0 {[ k k r ( h c ) ∫ ∂x j Γ2
s ij
（6）
s kr (hc ) [ A] = ∑ ∫ kij e =1 e ne
ne
∂N m ∂N n dv ∂x j ∂xi
式中： [ S ] =
k +1
（7）
采用隐式差分格式，有
[ A]{hc } + [ S ]{
k +1
− hc } = [F ] ∆t
k
（8）
4. 程序设计及部分说明
在程序中，单元采用六面体八节点等参单元[6]，其形函数为
1 N i (ξ ,η , ζ ) = (1 + ξiξ )(1 + ηiη )(1 + ζ iζ ) 8
-5-

图 3 计算模型网格图
4.5
3.5
2.5

0.5
图4
2h 渗流场等水头分布图（单位：米）
4.5
3.5
2.5
1.5
0.5
图5
100h 渗流场等水头分布图（单位：米）
图 4 是闸基非稳定渗流 2 小时后的等水头分布图， 图 5 是闸基非稳定渗流 100 小时后的 等水头分布图，由这两图可以看出，在上游水位从 3.0 米骤升到 5.0 米时，闸基土体的渗透 坡降没有超过允许值， 闸基没有发生渗透破坏。
J = Jx + J y + Jz
根据巴甫洛夫斯基的渗流控制理论，只要 J出口 ≤ J 允许 则闸基渗流是安全的，反之，若 每个单元的渗透坡降 J ≥ J 允许 ， 则可以把此单元看作是危险单元， 这种单元也是最容易发 生渗透破坏的。本程序根据这个原理，从非稳定渗流场初始时刻算起，计算时间步长为 2h 的渗流场水头值，并判断单元的渗透坡降是否超过允许值，若超过，则剔除此单元，反之保 留。通过这样的计算，可以看到闸基渗透破坏的发生与发展全过程。单元破坏后会形成新的 边界条件，而渗流破坏总是开始于渗流出口，所以在程序中，把下游水位作为破坏单元周围 的水头边界条件。闸基非稳定渗流破坏计算程序框图见图 1。
s
C 和 S S 等材料特性参数均需通过试验获得。
3.2 定解条件
(1) 初始条件：
-2-

hc ( xi , 0) = hc ( xi , t0 ), i = 1, 2,3
(2) 边界条件：
（3）
hc ( xi , t ) |Γ1 = hc1 ( xi , t )
-1-

一部分细颗粒没有紧密接触,甚至是处于自由状态,粗颗粒无法制约细颗粒,所以抗渗强度低。 它通常发生在砂砾石地层中。由于管涌破坏并不会直接造成土体结构破坏,所以其危险性不 大。但是管涌型土和流土型土的区别并不绝对的。当管涌型土发生管涌破坏,管涌通道贯通 上游后,若渗流坡降继续增大,管涌型土也会发生流土破坏。所以可以说,渗透破坏的最终形式 为流土破坏。
∂θ ， ∂t
在正压区 C = 0 ； hc 为压力水头； θ 为体积含水率； β 为饱和-非饱和选择常数，在非饱和 区等于 0，在饱和区等于 1； S S 为弹性贮水率，饱和土体的 S S 为一个常数，在非饱和土体 中 S S = 0 ，当忽略土体骨架及水的压缩性时则有 S S = 0 ； S 为源汇项。一般地 kij ， k r ，
坡降。
（12）
（13）
在形心处 ξ = η = ζ = 0 ，当已知八个节点的水头值时，就可以求出单元形心处的渗透
5. 算例
5.1 算例一：某闸基，底板长 20 米，厚 1 米，无上下游铺盖，闸基向上下游各延伸 30 米，
向地基延伸 10 米，取宽度 1 米的计算段。闸基土的颗分曲线如图 2（采用毛昶熙在文[7]中 试验砂样） ，渗透系数为 1.6 ×10 m / s ，单位贮存量为 2 × 10 3） 。