岩石力学与工程学报 22(4)：625～629
2003年4月 Chinese Journal of Rock Mechanics and Engineering April ，2003
2002年5月14日收到初稿，2002年7月6日收到修改稿。 * 国家自然科学基金重大研究计划资助项目(90102002)。
作者 戚国庆 简介：男，35岁，博士，19年毕业于河海大学工程勘测系水文地质及工程地质专业，现任高级工程师，主要从事地质工程的研究工作。
岩质边坡降雨入渗过程的数值模拟
*
戚国庆 黄润秋 速宝玉 胡云进 詹美礼
(成都理工大学环境与土木工程学院 成都 610059) (河海大学水电学院 南京 210098)
摘要 对岩质边坡降雨入渗的数学模型及模拟方法进行了探讨。在此基础上，研究了非饱和-饱和渗流模型中的非线性方程组及其解法，给出了Picard 迭代法求解该方程组的方法、步骤，并对边坡降雨入渗的规律进行了研究，阐述了岩质边坡降雨入渗过程中基质吸力的变化、暂态饱和区的形成、发展以及暂态水压力的分布、变化。用上述方法分析、模拟了某露天矿边坡降雨入渗的全过程。
关键词 数值分析，非饱和-饱和渗流，基质吸力，暂态饱和区，暂态水压力，Picard 迭代法。 分类号 O 241，TD 824.7 文献标识码 A 文章编号 1000-6915(2003)04-0625-05
NUMERIC SIMULATION ON RAINFALL INFILTRATION ON ROCK SLOPE
Qi Guoqing 1， Huang
Runqiu 1， Su Baoyu 2， Hu Yunjin 2， Zhan Meili 2 (1College of Environment and Civil Engineering ，Chengdu University ， Chengdu 610059 China )
(2College of Hydropower Engineering ，Hohai University ， Nanjing 210098 China
)
Abstract The mathematical model of rainfall infiltration on the rock slope and the numerical simulation method
for this model is discussed. Based on this ，the system of nonlinear equations in saturated-unsaturated seepage model is researched. The solution method and process of the system of nonlinear equations with Picard iteration method is given. The pattern of rainfall infiltration on the rock slope is analyzed. The variation of matrix suction ，the formation and development of transient saturated zone ，and the distribution and variation of transient water pressure in the process of rainfall infiltration on the rock slope are described. All the stages of rainfall infiltration on rock slope in a given open pit mine is simulated with the presented method.
Key words numerical analysis ，saturated-unsaturated seepage ，matrix suction ，transient saturated zone ，
transient water pressure ， Picard
iteration
1 引 言
基于“天然岩坡和铁路、公路、露天矿等开挖边坡的滑坡以及岩体中地下硐室的塌方多发生在雨季[1]”的事实，近年来，人们开始研究降雨对边坡稳定性的影响，并认为：降雨入渗会导致边坡非饱和区水压力的暂时升高，产生暂态的附加水荷载(暂态水压力)；同时，在非饱和区中出现暂态饱和
区。暂态饱和区与暂态水压力虽然是暂时的，但对
边坡稳定而言却是至关重要的，因为，暂态附加水荷载远比稳态水荷载增量大，常成为边坡失稳的控制因素，而且，裂隙发育的岩质边坡的渗流场受降雨影响往往比土坡更敏感[2
，3]
。
然而，目前对这种暂态水压力的升高尚无可应用的定量评价成果，仅是定性地认为它主要取决于裂隙网络的几何特征、裂隙的产状、风化程度、降雨强度和历时。故已有的工程设计对这种暂态饱和
区及暂态水压力均采用假定值[4]。例如，美国一些工程采用安全饱和边坡考虑。对于高边坡工程，采用这一水压力分布使加固设计过于保守。我国一些边坡工程常将静水压力乘以折减系数来进行边坡设计，如漫湾折减系数为0.4；五强溪水电站则采用类似Hoek建议图形，但取值较小；三峡工程采用折减系数为0.3；目前采用强排水体系而不考虑暂态水压力。不同方法采用的边坡暂态水压力分布彼此差别极大，且这些方法均缺少理论分析及实测数据。由此可见，如何恰当估计降雨在边坡体内的入渗及其所产生的影响，是边坡稳定性评价和加固、治理的关键问题。
事实上，降雨入渗首先使得边坡非饱和区的含水量增高、基质吸力降低，然后形成暂态饱和区，在暂态饱和区中出现暂态水压力。依据文[5]，笔者认为：降雨入渗所产生的结果对边坡稳定性影响最大的应是非饱和区基质吸力的降低，其次是暂态水压力的升高。本文旨在研究岩质边坡降雨入渗的分析、模拟方法，并对岩质边坡在降雨入渗情况下非饱和区基质吸力的变化、暂态饱和区的形成、发展以及暂态水压力分布、变化情况进行分析、模拟。以期能为降雨入渗对岩质边坡稳定性影响的定量评价提供一种新方法。
2 模拟岩质边坡降雨入渗的数学模型
降雨渗入边坡岩体到达潜水面经历了一个非饱和-饱和渗流过程[6]。目前已有求解裂隙岩体非饱和渗流场的数学模型，归纳起来可划分为以下4种[7]：(1) 等效连续介质模型，就是根据流量相等原则把岩块-裂隙系统等效成连续介质[8]，它适用于裂隙发育的大区域准恒定流问题的数值模拟；(2) 离散裂隙网络模型，认为岩块本身不透水，整个地下水运动是通过裂隙网络来进行的[4]，较适合于求解岩块致密、裂隙稀疏的小区域的非饱和渗流问题；(3) 双重介质模型，是一种双连续介质模型，认为岩块孔隙系统和裂隙系统(把裂隙网络等效为连续介质)均连续充满整个研究域[9]，它具有较好的拟真性(相对于等效连续介质模型)和较好的可操作性(相对于离散裂隙网络模型)；(4) 离散介质-连续介质耦合模型，其基本思路是，用离散裂隙网络模型描述主干裂隙中的水运动，用等效连续介质模型描述次要裂隙和孔隙中的水运动[10]。这样，既能反映裂隙的导水作用，又能体现岩块的贮水作用。但是，该模型数学处理比较困难，水交换量难以准确确定。迄今为止，在工程及生产科研中应用最为广泛的是等效连续介质模型。
2.1控制方程及其代数方程组
裂隙发育的岩质边坡的降雨入渗问题可采用等效连续介质模型进行分析，其控制方程如下：
t
h
S
h
C
R
x
h
x
K
i
j i
i j i
h
x∂
∂
+
=
+
+
∂
∂
∑∑
==
∂
∂
]
)
(
[
)]
(
)
(
[3
3
1
3
1
β
(1) 式中：h 在饱和区为压力水头，在非饱和区为基质
势；K(h)
i j
为渗透系数，在非饱和区为基质势h的函数；C(h)为容水度，它反映了基质势h与饱和度的关系；R 为源(汇)项；S为贮水系数；β在饱和区为0 ，在非饱和区为1；x1，x2，x3分别表示笛卡尔坐标系中的轴x，y，z，轴z(x3)为正向向上的铅直轴。
采用有限单元法对降雨入渗过程进行数值模拟，首先采用Galerkin有限单元法将控制方程式(1)离散为代数方程组：
−
−
−
=
Δ
++
=
∑n
n
n
k
m
nm
NP
m
nm
D
B
Q
h
t
A1
1
]
1
2
1
[
k
m
nm
nm
NP
m
h
F
t
A]
1
2
1
[
1
Δ
−
∑
=
(2)
式中：k
m
k
m
h
h，
1+分别为k+1，k时刻的水势，在非饱和区为基质势，在饱和区为压力水头；t
Δ为时间步长；]
1
2
1
[
1
nm
NP
m
nm
F
t
A
∑
=
Δ
+为系数矩阵，且
∑∑∑∫∫∫
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
∂
∂
∂
∂
=
==
e
3
1
3
1
e
e
e d
)
(
e
G
x
N
x
N
h
K
A
i j G j
m
i
n
j i
k
nm
⎪⎩
⎪
⎨
⎧
≠
=
+
=
∑∫∫∫
m
n
m
n
G
N
N
S
h
C
F
G
m
n
k
nm
d
]
)
(
[
e
e
e
e
e
e
β
S
n
h
K
h
N
x
h
K
N
Q
i
j
i
k
k
m
m
i
j i
k
i
n
n
d
)
(
)
(
)
(
3
1
3
e
e
e
3
1
e
e
2
r
⎥
⎦
⎤
⎢
⎣
⎡
+
∂
∂
=
∑
∑
∑∫∫
=
=
•
Γ
∑∑∫∫∫
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
∂
∂
=
=
e
3
1
e
3
e
e
d
)
(
i G i
n
i
k
n
G
x
N
h
K
B
G
N
R
D
G
n
n
d
e
e
e
∑∫∫∫
=
2.2非线性方程组求解方法讨论
第22卷 第4期 戚国庆等. 岩质边坡降雨入渗过程的数值模拟 • 627 •
由于非饱和-饱和渗流有限元计算格式(式(2))
中的系数矩阵依赖于未知量(k m k m h h ，1+)而变化，即在
非饱和区介质的渗透系数K (h )ij 和容水度C (h )为介质的基质势h 的函数，故式(2)为一非线性方程组，需通过迭代法来求解。收敛速度较快的Newton-
Raphson 法在迭代过程中必须求解Jacobi 矩阵，不但迭代每一步的工作量很大，而且有时无法求得
Jacobi 矩阵(当根据基质势h 和试验数据来线性内插K (h )ij 和C (h )时)或Jacobi 矩阵很难用解析表达式来表示(当K (h )ij 和C (h )与h 的关系式很复杂时)。为确保收敛速度并减少工作量，本文采用Picard 迭代 法。具体迭代方法和步骤如下：
(1) 设定一个迭代收敛条件：max
1k
l k l h h ++≤ε，其中上标k 代表第k 时步，下标l +1和l 分别代
表第l +1和l 次迭代得到的值，ε为前后两次迭代计算结果间的误差允许值。
(2) 根据t k －1和t k 时刻的渗流场h k －1和h k ，线
性外推待求时刻t k +1的初始渗流场，即=+10k h h k (=k 0时)或=+1
k h h k +(Δt k /Δt k －1)(h k －h k －1)(k >0时)，这里Δt k =t k －1－t k 。 (3) 由t k+1/2
时刻的渗流场2)(12/1+++=k l k k l h h h 确定出各单元的渗透系数K (h )i j 和容水度C (h ) 后，计算有限元计算格式中的A nm ，F nm ，Q n ，B n 和
D n ，用Crout 分解法求解式(2)，得到新的渗流 场11++k l h 。
(4) 将111+++k l k l h h 和进行比较，若max
1
1
1++++k l
k l h h ≤ε，结束迭代进入下一个时步。否则，继续重复进行第(3)和(4)步直至满足收敛条件，即可求得降雨 入渗过程中边坡非饱和-饱和渗流场、非饱和带暂态饱和区、暂态水压力以及基质吸力。
3 岩质边坡降雨入渗的规律研究
3.1 某露天矿岩质边坡降雨入渗的数值模拟
某露天采矿边坡由花岗岩(γ)，变质岩(Su)，矿层(Fe)组成(如图1)。地下水埋深60 m 左右。非饱和带位于强风化岩体内[11]。
非饱和带渗透系数：=K 3.94×10－
2×|h |
－1.047
(m/d)，饱和渗透系数=K 3.94×10－
2(m / d)，容水度C (h )=0.011 9×|h |
－0.672
；在饱和带中：花岗岩(γ)渗
透系数=K 5.59×10－
2(m/d)，变质岩(Su)渗透系数
图1 某露天矿边破状况
Fig.1 Situation of the slope in a certain open pit mine
=K 8.21×10
－2
(m/d)，矿层(Fe)渗透系数=K 4.06
×10－
2 (m/d)。依据所选参数及初始条件、边界条 件，模拟某矿边坡在7 d 最大降雨为676.1 mm 的情况下，非饱和渗流场的变化情况(如图2)。 3.2 岩质边坡降雨入渗规律的研究
模拟结果反映了边坡降雨入渗状态下非饱和-饱和渗流场的变化规律。
(1) 边坡非饱和区，降雨历时1 d 以后，首先在饱和区逸出点以上附近出现暂态饱和区。随降雨历时的增加，暂态饱和区范围沿坡面线附近区域不断扩展、延伸。至降雨历时4 d 后，在整个坡面线附近区域形成了一个由暂态饱和区构成的“饱和壳”，并随着降雨过程的延续向边坡内部扩展，影响范围
逐渐扩大(如图2)。
(2) 随着降雨不断入渗，暂态饱和区内出现正的水压力，此即为通常所说的暂态水压力。暂态水压力在坡面线附近数值较小，在暂态饱和区下缘与非饱和区接触的地带为零。最大值出现在非饱和区中部偏上部位(如图2中A -A ′断面)。在整个7 d 的降雨过程中，暂态水压力随降雨历时的增加而逐渐增大，其最大值为25.66 kPa(相当于2.566 m 压力 水头)。暂态水压力随降雨时间的变化如图3。
(3) 随着降雨入渗边坡，非饱和区上部的含水量逐渐增大，基质吸力逐渐降低。非饱和区最大基质吸力随降雨时间的变化如图4。
(4) 在7 d 连续降雨过程中，
暂态饱和区未延伸到潜水面，因而饱和区渗流场的变化不大，潜水面升高很小，这与以往运用饱和渗流模型模拟降雨入渗的结果有很大差别。
100
200
300
500
水平距离/m
高程/m
Fig.3 Variation of transient water pressure in the process of rainfall
图4 降雨过程中，基质吸力的变化Fig.4 Variation of matrix suction in the process of rainfall
降雨历时/d 01234567
降雨历时/d
• 626 • 岩石力学与工程学报 2003年
4 结论
(1) 文章对岩质边坡降雨入渗的模拟方法进行了探讨、研究，对降雨入渗过程中，边坡非饱和区基质吸力的变化，暂态饱和区的形成、发展、变化及暂态水压力的大小、分布等进行了数值模拟和分析。通过分析发现边坡降雨入渗过程的规律为：非饱和带基质吸力不断降低；暂态饱和区首先从饱和区逸出点以上附近及坡顶开始形成，再向边坡内部扩展，影响范围逐渐扩大；在暂态饱和区中出现暂态水压力，其最大值出现在非饱和区中部偏上部分，随降雨历时的增加而逐渐增大；而潜水面的升高幅度不大，这与以往运用饱和渗流模型对降雨入渗直接补给潜水的认识是完全不同的。
(2) 针对非饱和-饱和渗流数值模拟中求解非线性方程组的问题，文章采用了Picard迭代法求解的方法、步骤。通过文中的算例分析和工程应用表明，Picard迭代法用于求解非饱和-饱和方程不仅收敛快，而且工作量小。
参考文献
1 孙广忠. 中国典型滑坡[M]. 北京：科学出版社，1998
2 费里泽R A，彻里J A. 地下水[M]. 北京：地震出版社，1987
3 陆晓敏，任青文，盛芳. 裂隙岩质边坡的弹粘塑性变形及稳定性
分析[J]. 岩石力学与工程学报，2002，21(4)：493～497
4 张有天，刘中. 降雨过程裂隙网络饱和/非饱和、非恒定渗流分
析[J]. 岩石力学与工程学报，1997，16(2)：104～111
5 弗雷德隆德D G，拉哈尔佐H. 非饱和土力学[M]. 陈仲颐译. 北
京：中国建筑工业出版社，1997
6 张蔚榛. 地下水与土壤水动力学[M]. 北京：中国水利出版社，1996
7 胡云进，速宝玉，詹美礼. 裂隙岩体非饱和渗流研究综述[J]. 河海
大学学报，2000，28(1)：45～51
8 Dykhuizen R C. Transport of solutes through unsaturated fractured
media[J]. Water Res.，1987，21(12)：1531～1539
9 Salam A A，Chrysikopoulos C V. Unsaturated flow in a quasi-three-
dimensional fractured medium with spatially variable aperture[J].
Water Resour. Res.，1996，32(6)：1 531～1 540
10 王媛，速宝玉，徐志英. 三维裂隙岩体渗流耦合模型及其有限元
模拟[J]. 水文地质工程地质，1995，22 (3)：1～5
11 戚国庆. 裂隙岩体非饱和带地下水渗流研究及其在露天矿边坡稳
定性评价中的应用[硕士学位论文][D]. 南京：河海大学，2000
关于在学术论文中规范关键词选择的规定(试行)
一、发表在中国科协系统学术期刊中所有学术论文，必须在摘要后列出不少于4个关键词。从技术角度考虑，没有关键词的论文应列入非学术论文类。
二、这些关键词按以下顺序选择：
第一个关键词列出该文主要工作或内容所属二级学科名称。学科体系采用国家技术监督局发布的学科分类与代码（国际GB/T13745-92）。
第二个关键词列出该文研究得到的成果名称或文内若干个成果的总类别名称。
第三个关键词列出该文在得到上述成果或结论时采用的科学研究方法的具体名称。对于综述和评述性学术论文等，此位置分别写“综述”或“评论”等。对科学研究方法的研究论文，此处不写被研究的方法名称，而写所应用的方法名称。前者出现于第二关键词的位置。
第四个关键词列出在前三个关键词中没有出现的，但被该文作为主要研究对象的事或物质的名称，或者在题目中出现的作者认为重要的名词。
如有需要，第五、第六个关键词等列出作者认为有利于检索和文献利用的其他关键词。
中国科协学会学术部
2002年8月22日
第22卷第4期戚国庆等. 岩质边坡降雨入渗过程的数值模拟 • 629 •下载本文
