有密集排水孔的三维饱和-非饱和渗流场分析

第23卷 第12期
岩石力学与工程学报 23(12)：2027～2031
2004年6月 Chinese Journal of Rock Mechanics and Engineering June ，2004
2002年5月21日收到初稿，2002年7月3日收到修改稿。

* 教育部年轻骨干教师项目“降雨入渗对边坡稳定性影响研究”和南京市留学回国人员基金资助项目。

作者 陈建余 简介：男，1974年生，博士，主要从事水工结构工程及岩土工程的渗控、稳控等方面的研究工作。

E-mail ：2000chenjianyu@ 。

有密集排水孔的三维饱和-非饱和渗流场分析
*
陈建余
(广州市水利局 广州 510640)
摘要 针对工程实际存在的饱和-非饱和边坡岩体等计算域内有密集排水孔的问题，引入排水子结构技术，推导了相应的计算公式，并开发了相关的程序，同时也简介了降雨入渗和饱和-非饱和渗流逸出面的正确模拟问题。

最后对长江三峡永久船闸边坡降雨入渗饱和-非饱和渗流场问题进行了求解，分析和论证了工程中的有关渗流控制情况。

关键词 渗流力学，有限单元法，饱和-非饱和，边界条件，排水子结构
分类号 O 357.3，P 641.2 文献标识码 A 文章编号 1000-6915(2004)12-2027-05
3D FEM ANALYSIS OF SATURATED-UNSATURATED SEEPAGE FIELD
WITH DENSELY-DISTRIBUTED DRAINAGE HOLES
Chen Jianyu
(Water Conservancy Bureau of Guangzhou City ， Guangzhou 510640 China )
Abstract To the unsaturated-saturated seepage field with densely-distributed drainage holes existing in the computational domain ，the drainage substructure is introduced. The involved formulae are deduced and the com-putation program is developed. Meanwhile ，rainfall infiltration and precise simulation of the exit face of the satu-rated-unsaturated seepage are also introduced. Finally ，the problem of the slope seepage field of the permanent shiplock of the Three Gorges Project is calculated. The effectiveness and propriety of the seepage control measures in the slope are analyzed and verified.
Key words seepage mechanics ，FEM ，saturated-unsaturated ，boundary condition ，drainage substructure
1 引 言
在堤坝、岩土质边坡等工程中经常涉及到用有
限元法求解三维饱和-非饱和非恒定流渗流场问题，Miller 于1962年提出了非饱和介质渗透系数是含水量或压力水头的函数，文[1]据此于1973年先后提出了一维、二维的饱和-非饱和渗流计算有限单元法，Akai 于1979年提出了三维饱和-非饱和渗流计
算的有限单元法，另外，文[2，3]都对此进行了相当深入的研究，取得了一定的成果。

目前，这一问题的基本理论及其求解方法已取得了比较一致的认同，但对渗流场的边界条件，包括蒸发与入渗边界以及饱和-非饱和逸出面，还有计算域内部存在众多排水孔时排水孔内部边界的精细模拟也是求得渗流场高精度解的保证之一，边界条件的全真模拟成为解决三维饱和-非饱和渗流场问题的难点和关键。

本文在文[4]提出的改进排水子结构法的基础上，首次
• 2028 • 岩石力学与工程学报 2004年
将排水子结构技术应用到饱和-非饱和非恒定流渗流场问题中，推导了相应的公式，并开发了相应的程序，同时，结合长江三峡永久船闸边坡渗流场问题，验证了相关理论、算法和程序是科学和正确的。

2 渗流场求解基本理论
2.1 微分方程
多孔介质不可压缩时饱和-非饱和非恒定流渗
流场求解的微分方程[1，
5]为
t h S h C S h k k x h
h k k x i j
ij i ∂∂+=++∂∂∂∂c s c c r s 3c c r s ))(()]()([β
(1)
式中：s
ij
k 为饱和渗透系数张量；r k 为相对透水率，
为非饱和土的渗透系数与同一种土饱和时的渗透系数的比值，在非饱和区0＜r k ＜1，在饱和区=r k 1；
C 为容水度，c
h C ∂∂=θ
，在正压力区0=C ；c h 为压 力水头；θ为含水率；β为饱和-非饱和常数，非饱
和区等于0，饱和区等于1；s S 为弹性贮水率，饱和土体s S 为一个常数，非饱和土体=s S 0，当忽略
骨架及水的压缩性时可设=s S 0；S 为源汇项；s
ij k ，
r k ，C ，s S 等参数可通过试验获得。

2.2 定解条件
(1) 初始条件：
)()0 (0c c t x h x h i i ，，= (2)
(2) 边界条件：
)(|)0 (01c c 1t x h x h i Γi ，，
= (3) n q n h k k x h
h k k Γi i j
ij
=+∂∂−2|)]()([c r s 3c c r s (4)
3|)]()([c r s 3c c r
s Γi i j
ij n h k k x h
h k k +∂∂−≥0
且
0|3c =Γh (5)
)(|)]()([4c r s 3c c r s t q n h k k x h
h k k Γi i j
ij
θ=+∂∂−
且
4|c Γh ＜0 (6)
式中：n q 为已知边界流量；)()(e θθλθ−=t q ，λ为渗流场中非饱和逸出面上的水量逸出(入)渗系数，
e θ为介质在长期的环境湿度条件a θ下所拥有的极
限含水率，)(a e e θθθ=，需由试验测定，本文仅考虑
降雨入渗边界，)()(t i t q r =θ，)(t i r 为入渗率；i n 为外法线方向余弦；0t 为初始时刻；1Γ为已知水头边界；2Γ为流量边界；3Γ为饱和逸出面边界；4Γ为非饱和逸出面边界。

3 渗流场有限元法求解
对于式(1)～(6)组成渗流场问题，由Galerkin 加权余量法可得求解的积分方程为 ∫−⎪⎭
⎪⎫⎪⎩⎪⎨⎧∂∂−⎥
⎦⎤⎢⎣⎡+∂∂∂∂
e n i j ij i
v N t h h C h k k x h h k k x c c c r s
3c c r s d )()()( −⎪⎭⎪⎬⎫+⎪⎩⎪⎨⎧⎥⎦⎤⎢⎣⎡+∂∂∫s N q n h k k x h h k k n n i Γi j ij e d )()(2 c r s
3c c r s 0d )()()(4 c r s
3c c r s =⎪⎭
⎪⎬⎫+⎪⎩⎪⎨⎧⎥⎦⎤⎢⎣⎡+∂∂∫s N t q n h k k x h h k k n i Γi j ij e θ (7) 对式(7)左第1项用格林第1公式，同时取=c h m m h N c ，此处m N 为形函数，m h c 为结点压力水头，
并令v x N x N h k k A i n j m NE
e e ij d )(][c r 1 s
∂∂∂∂=∑∫=，
=
][S ∑
=NE
e 1
∫Δe t h C c )
(v N N n m d ，
i n NE
e e i x N
h k k F ∂∂⋅−=∑∫=)(][c r
1
s 3−
v d ⋅∑∫
=NE
e Γe n q 1
2
s N t q s N n NE
e Γn e d )(d 1
4∑∫=−θ
则由式(7)可以得到式(1)时间差分型式的有限元求解格式为
}{}]{[}]){[]([c c F h S h S A t t t +=+Δ− (8)
4 排水子结构技术[4
～6]
4.1 排水孔渗流行为分析
计算域内大量的排水孔对渗流场的渗流特性影响极大，对排水孔的准确模拟直接关系到渗流场的精确求解，文[3]于1992年率先提出了排水子结构的概念，文[4]在此基础上提出了在饱和渗流场求解中理论和算法上都严密的改进排水子结构技术，而本文则试图将文[4]提出的改进排水子结构技术引入到饱和-非饱和非恒定流渗流场问题的求解中，解决饱和-非饱和非恒定流渗流场中密集型排水孔的
第23卷 第12期 陈建余. 有密集排水孔的三维饱和-非饱和渗流场分析究 • 2029·
精细模拟问题。

排水孔的渗流行为可分为3种情况：(1) 排水孔位于渗流自由面上方，排水孔内部边界全部按非饱和逸出面处理；(2) 排水孔与渗流自由面相交，排水孔内部边界出现非饱和逸出面、饱和逸出面和已知水头边界3种情况；(3) 排水孔位于渗流自由面下方，排水孔内部边界全部按已知水头边界处理。

第2种情况的排水孔边界条件包含了第1，3种 情况，因此，只需分析第2情况(图1)。

Fig.1
(1) ⎩⎨⎧−)(c r h k [(2) 满足下式的要求： ⎭
⎬⎫⎩⎨⎧+⎦
⎤⎢⎣⎡∂∂∂∂∂∂−][][)(s
c c c
c r zz yz
xz k k k k z h y
h x
h
h k · T ][z y x n n n ≥0 (11)
0c =h (12) (3) f i 或ik 段为已知水头边界，其渗流行为满
足下式的要求：
z z h −=0c (13)
4.2 排水子结构及其出口传导矩阵和右端项
子结构的划分型式一般没有特别的限制，但需要将计算域排水孔实际分布划分为合理的母单元，图2为常用的包含3个排水孔和5个排水孔的排水子结构划分型式。

对于任意某一排水子结构，假定初始条件)()0 (0c c t x h x h i i ，，
=，有限元方程格式为 }{}]{[t c q h K =
式中：
NES
NES
[K ]为子结构有限元方程总刚度矩阵，
}{q 为右端项。

设子结构内部结点未知水头列阵和右端项为h 1和q 1，出口结点未知水头列阵和右端项为h 2和q 2(含邻
域单元的结点流量贡献)，将[K ]及}{q 作相应分块，则有
⎭
⎬⎫
⎩⎨⎧=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡212122211211q q h h k k k k (14) )(21211
11
1h k q k h −=− (15)
• 2030 • 岩石力学与工程学报 2004年
2222q h k ′=′ (16)
k k k k k 1
11212222−−=′ (17)
q k k q q 111
2122−−=′
(18)
式中：11k ，12k ，21k 和22k 为按子结构内结点划分及出口结点而对子结构传导矩阵进行分块的子矩阵；22
k ′，q ′为子结构静凝聚后出口传导矩阵及出口右端项。

5 工程算例
5.1 工程概述
长江三峡永久船闸系开挖在大坝左岸山体中，分北线和南线的5级双线船闸，闸室壁面开挖成直立陡高坡，闸面筑有厚约1.8 m 的贴壁式混凝土衬砌边墙，墙后山体开挖成多台阶的岩体边坡。

边坡地下水主要来自山体中的地下水、地表降雨入渗及闸室高水位运行时内水外渗。

闸室中高水位骤降 时，闸室衬墙有可能承受较大的外水荷载，为此，在墙体背后及边坡岩体内布置了多种排水渗控措施。

主要渗控措施有直墙后贴面网格式排水管排水系统，沿闸室水流向及高程向布置成4 m ×7.5 m 的网格型式，排水管直径按=φ30 cm 考虑。

图3为 第2闸室边坡体计算模型及单元网格剖分图。

图3 计算模型及单元网格分剖分图 Fig.3 Calculation model and its mesh drawing
5.2 计算结果分析
图4为无降雨时第2闸室北线边坡稳定渗流
场等压力线图，这说明边坡岩体中的排水幕布置至关重要，排水渗控措施很好地隔断了山体地下水 对边墙的作用。

图5为雨强为100 mm/d 降雨历时 4 d 时边坡渗流场等压力线图，这说明此时渗流饱和
区从顶层排水廊道以上区域越过排水幕逐渐向边 坡体表面扩展，同时渗流自由面也逐渐上抬，依然显现了排水廊道及排水孔等排水渗控措施的重要 性。

图4 无降雨时边坡稳定渗流场等压力线图 Fig.4 Pressure isolines of steady seepage field of the slope
without rainfall
图5 雨强为100 mm/d 降雨历时4 d 时边坡渗流场等压力
线图
Fig.5 Pressure isolines of slope seepage field after 4 days of
rainfall at 100 mm/d
6 主要结论
(1) 饱和-非饱和非恒定流渗流场精确计算的关
单位：m
单位：m
第23卷第12期陈建余. 有密集排水孔的三维饱和-非饱和渗流场分析究• 2031·
键是计算域边界条件的正确严密模拟，引入排水子结构技术可以很好地解决密集排水孔的精细模拟问题。

(2) 用排水子结构技术可以准确模拟排水孔的几何形状、分布特征、是否失效等所影响的排水孔不同渗流行为，理论上严密，技术上可靠，操作上可行。

(3) 饱和渗流逸出面是一个可变流量边界，目前一般只考虑降雨入渗边界，对于其他蒸发或入渗边界由于缺乏准确的计算参数通常只能作近似处理，这在一定的程度上影响了饱和-非饱和渗流场的计算精度。

(4) 通过对长江三峡永久船闸边坡饱和-非饱和非恒定流渗流场的计算分析，证明了研发的程序计算结果符合工程实际，在理论和算法上是正确的。

参考文献
1 Neuman S P，Narasimhan T N，Witherspoon P A. Application of mixed
explicit-implicit finite element method to nonlinear diffusion-type problems[A]. In：Finite Elements in Water Resources：Proc. of the First International Conference on Finite Elements in Water Resources[C].
London：Pentech Press，1977，1 153～1 186
2 毛昶熙，段祥宝，杨祖贻. 渗流数值计算与程序应用[M]. 南京：河
海大学出版社，1999
3 吴良骥. 饱和-非饱和区中的渗流问题的数值模型[J]. 水利水运科
学研究，1985，(2)：1～12
4 朱岳明，张燎军. 渗流场求解的改进排水子结构法[J]. 岩土工程学
报，1997，19(2)：69～76
5 王镭，刘中，张有天. 有排水孔幕的渗流场分析[J]. 水利学报，
1992，(4)：15～20
6 姜晋庆，张铎. 结构弹塑性有限元分析法[M]. 北京：宇航出版
社，1990
2020年我国全面实现小康综合国力跻身世界“三强”
中国科学院可持续发展战略研究组3月2日发布的《2004中国可持续发展战略报告》提出了2020年我国全面实现小康社会总体目标，以供有关部门作决策参考：用20年左右的时间，GDP总量在2000年的基础上翻两番，基本完成工业化，我国综合国力进入世界前3名的行列，同时基本达到可持续发展良性循环的要求。

课题组首席科学家牛文元说，到那时，我国的城市化率将从现在的36%提高到55%以上，具有容纳7亿至7.5亿人口的城市容量，形成结构合理、功能互补、整体效益最大化的大、中、小“城市体系”。

“国家发展将主要依靠3大组团式城市群、7大沿江沿海沿路的城市带、若干个具有竞争能力的中心城市圈”。

牛文元说，到2020年，我国东部沿海地区要走过3大“零增长”台阶，即依次在2010年左右实现人口自然增长率的“零增长”；2015年左右实现资源和能源消耗速率的“零增长”；2020年左右实现城市生态环境退化速率的“零增长”。

中部地区再用15年走过3大“零增长”的台阶，西部地区再用30年走过3大“零增长”的台阶。

到2020年，4大系数也将在我国发生变化：恩格尔系数不超过0.35，基尼系数保持在0.35～0.40的水平，人文发展指数不低于0.80，城乡二元结构指数限制在1.80以下。

到2020年，中国人平均预期寿命将达到76岁，人均受教育年限超过10年，科技创新能力指数平均达25，科技在整个国民经济中的贡献率达到65%以上。

(摘自2004年3月3日《长江日报》第3版)。