岩质边坡降雨入渗过程的数值模拟.pdf

第 2 2卷第 4期 2 0 0 3年 4月 岩石力学与工程学报 C h i n e s e J o u r n a l o f R o c k Me c h a n i c s a n dE n g i n e e r i n g 2 2 4 6 2 5 6 2 9 A pr i l ，2 0 0 3 岩质边坡降雨入渗过程的数值模拟 戚 国庆黄润秋 成都理工大学环境与土木工程学院成都6 1 0 0 5 9 速宝玉胡云进詹美礼 河海大学水电学院南京2 1 0 0 9 8 摘要对岩质边坡降雨入渗的数学模型及模拟方法进行了探讨。在此基础上，研究了非饱和． 饱和渗流模型中的非 线性方程组及其解法，给出了 P i c a r d迭代法求解该方程组的方法、步骤，并对边坡降雨入渗的规律进行了研究， 阐述 了岩质边坡降雨入渗过程中基质吸力的变化、暂态饱和区的形成、发展以及暂态水压力的分布、变化。用上 述方法分析、模拟了某露天矿边坡降雨入渗的全过程。 关键词数值分析，非饱和． 饱和渗流，基质吸力，暂态饱和区，暂态水压力，P i c a r d迭代法。 分类号0 2 4 1 ，T D 8 2 4 ． 7 ． 文献标识码A 文章编号1 0 0 0 - 6 9 1 5 2 0 0 3 0 4 ． 0 6 2 5 - 0 5 ． RI C S D 【 ]【 Ⅱ A TI oN oN RAD I D唧TR棚oN oN RoCK S LoP E Qi G u o q i n g ， H u a n g R u n q iu ， S u B a o y u ， H u Y u n j in ， Z h a n Me i li I C o l l e g e o f E n v i r o n m e n t a n dC i v i l E n g i n e e r i n g ，C h e n g d u U n i v e r s i ty ， C h e n g d u 6 1 0 0 5 9 C h i na 2 C o l l e g e o fH y d r o p o w e r E n g i n e e r i n g ，Ho h a i U n i v e r s i ty， N a n j i n g 2 1 0 0 9 8 C h i n a Ab s t r a c t Th e ma t h e ma t i c a l mo d e l o f r a i n f a l l i n fi l t r a t i o n o n t h e r oc k s l o p e a n d t h e n u me r i c a l s i mu l a t i o n me t h o d f o r t h i s mo d e li S d i s c u s s e d ．Ba s e d o n t h i s ， t h e s y s t e m o f n o n l i n e a r e q u a t i o n s i n s a t u r a t e d ． u n s a t u r a t e d s e e p a g e mo d e li S r e s e arc h e d ． Th e s o l u t i o n me t h od a n d p r o c e s s o f t h e s y s t e m o f n o n l i n e ar e q u a t i o n s wi t h P i c ard i t e r a t i o n me t h o d i S g i v e n． Th e p a t t e r n 0 f r a i n f all i n fi l t r a t i o n o n t h e r oc k s l o p e i S an a l y z e d． Th e v a r i a t i o n o f ma t r i x s u c t i o n， t h e f o r ma t i o n an d d e v e l o p me n t o f t r an s i e n t s a t u r a t e d z o n e ， an d t h e d i s t il b u t i o n an d v ari a t i o n o f tran s i e n t wa t e r p r e s s u r e i n t h e p r oc e s s o f r a i n f a l l i n fil t r a t i o n o n t h e r o c k s l o p e a r e d e s c r i b e d ． Al l t h e s t a g e s 0 f r a i n f a l l i n fil t r a t i o n o n r ock s l o p e i n a g i v e n O p e n p i t mi n e i S s i mu l a t e d wi t h the p r e s e n t e d me t h o d ． Ke y wo r d s n u me r i c a l an a l y s i s ， s a tur a t e d - u n s a tur a t e d s e e p a g e ， ma t r i x s u c t i o n ， t r an s i e n t wa t e r p r e s s u r e ，P i c ard i t e r a t i o n 1 引 言 基 于 “ 天然岩坡和铁路、公路、露天矿等开挖 边坡的滑坡以及岩体中地下硐室的塌方多发生在雨 季【 l I ’’ 的事实，近年来，人们开始研究降雨对边坡 稳定性的影响，并认为降雨入渗会导致边坡非饱 和区水压力的暂时升高， 产生暂态的附加水荷载 暂 态水压力 ；同时，在非饱和区中出现暂态饱和 t r a n s i e n t s a t u r a t e d z o n e ， 区。暂态饱和区与暂态水压力虽然是暂时的，但对 边坡稳定而言却是至关重要的 ，因为 ，暂态 附加水 荷载远 比稳态水荷载增量大 ，常成为边坡 失稳的控 制因素 ，而且 ，裂隙发育的岩质边坡 的渗流场受降 雨影响往往 比土坡更敏感【 2 _J 。 然而， 目前对这种暂态水压力的升高尚无可应 用的定量评价成果，仅是定性地认 为它主要取决于 裂隙网络的几何特征、裂隙的产状 、风化程度 、降 雨强度和历时。故 已有 的工程设计对这种暂态饱和 2 0 0 2年 5月 l 4日收到初稿。2 0 0 2年 7月 6日收到修改稿． 国家 自然科学基金重大研究计划资助项 I I 9 0 1 0 2 0 0 2 ． 作者 戚国庆 简介男，3 5岁 ，博士 ．1 9 8 9年毕业于河海大学工程勘测系水文地质及工程地质专业 ．现任高级工程师．主要从事地质工程的研究工 作． 维普资讯 6 2 6 岩石力学与工程学报 2 0 0 3 年 区及暂态水压力均采用假定值【 4 J 。例如，美 国一些 工程采用安全饱和边坡考虑。对于高边坡工程 ，采 用这一水压力分布使加固设计过于保守。我国一些 边坡工程常将静水压力乘 以折减系数来进行边坡设 计 ，如漫湾 折减 系数 为 0 ． 4 五强溪 水 电站则采 用类似 H o c k建议图形，但取值较小；三峡工程采 用折减系数 为 0 ． 3 ； 目前采用强排水体系而不考虑 暂态水压力 。不同方法采用 的边坡暂态水压力分布 彼此差别极大 ，且这些方法均缺少理论分析及 实测 数据。由此可见，如何恰当估计降雨在边坡体内的 入渗及其所产生的影响， 是边坡稳定性评价和加固、 治理的关键 问题。 事实上 ，降雨入渗首先使得边坡非饱和区的含 水量增高、基质吸力降低，然后形成暂态饱和区， 在暂态饱和 区中出现暂态水压力。依据文【 5 】 ，笔者 认为 降雨入渗所产生的结果对边坡稳定性影响最 大 的应是非饱和区基质吸力的降低，其次是暂态水 压力 的升高。本文 旨在研究岩质边坡降雨入渗的分 析 、模拟方法 ，并对岩质边坡在降雨入渗情况下非 饱和区基质吸力 的变化 、暂态饱和区的形成 、发展 以及暂态水压力分布、变化情况进行分析、模拟。 以期能为降雨入渗对岩质边坡稳定性影响的定量评 价提供一种新方法 。 2 模拟岩质边坡降雨入渗的数学模型 降雨渗入边坡岩体到达潜水面经历 了一个非饱 和． 饱和渗流过程【 o J 。目前 已有求解裂隙岩体非饱和 渗流场的数学模型， 归纳起来可划分为以下 4 种l 1 等效连续介质模 型，就是根据流量相等原则把 岩块． 裂隙系统等效成连续介质【 8 J ，它适用于裂隙发 育 的大区域准恒 定流问题 的数值模拟 2 离散裂 隙网络模型，认为岩块本身不透水 整个地下水运 动是通过裂隙网络来进行的l引 ，较适合于求解岩块 致密 、 裂 隙稀疏 的小区域的非饱和渗流 问题 3 双 重介质模 型，是一种双连续介质模型 ，认为岩块孔 隙系统和裂隙系统 把裂隙网络等效为连续介质 均 连续充满整个研究域l引 ，它具有较好的拟真性 相对 于等效连续介质模型 和较好的可操作性 相对于离 散裂隙 网络模型 4 离散介质． 连续介质耦合模 型，其基本思路是，用离散裂隙网络模型描述主干 裂隙中的水运动，用等效连续介质模型描述次要裂隙 和孔隙中的水运动【 l 叭 。这样 ，既能反映裂隙的导水 作用 ，又能体现岩块 的贮水作用。但是，该模型数 学处理比较困难 ，水交换量难 以准确确定。迄今为 止 ，在工程及生产科研 中应用最为广泛的是等效连 续介质模型。 2 ． 1 控制方程及其代数方程组 裂隙发育的岩质边坡 的降雨入渗 问题可采用等 效连续介质模 型进行分析 ，其控制方程如下 喜 喜 亭 砉 圳 1 式中h在饱和区为压力水头，在非饱和区为基质 势； 为渗透系数 ，在非饱和区为基质势 h的函 数；C 为容水度 ，它反映了基质势 h与饱和度 的 关系；R 为源 汇 项；S 为贮水系数； 在饱和 区 为 0 ，在非饱和区为 1 ， ，X 3分别表示笛卡 尔坐标系中的轴 ， Y ， z ，轴 z x 3 为正 向向上 的铅直 轴。 采 用有 限单元 法对 降雨 入渗过 程进行 数值模 拟 ，首先采用 Ga l c r k i n有限单元法将控制方程式 1 离散为代数方程组 】 椿 一 一 一 吉 4 w 一 l 2 式 中， ， 分别为 1 ，k时刻的水势 ，在非饱 和区为基质势 ，在饱和区为压力水头；△f 为时间步 长 ； 畦 4 m F -- ] 为 系 数 矩 阵 ， 且 i t ， 等 l e l f’ 6t f l FJ ∑ 肌c 。 ‘ p S 。 】 Ⅳ Ⅳ t d G 刀 J 一 、 0 G 0 【 0 刀 ≠ m ∑ Ⅳ ∑． 匡 砉 c I f『J 筹a c l 。 L 6． J ∑『 『 』 R d G 2 ． 2 非线性方程组求解方法讨论 维普资讯 第2 2 卷第4期 戚国庆等．岩质边坡降雨入渗过程的数值模拟 6 2 7 由于非饱和． 饱和渗流有限元计算格式 式 2 中的系数矩阵依赖于未知量 砧“ ， 而变化， 即在 非饱和区介质的渗透系数 和容水度 C III 为介 质的基质势 h的函数，故式 2 为一非线性方程组， 需通 过迭代法来 求解 。收敛速度 较快的 Ne wt o n ． R a p h s o n法在迭代过程中必须求解 J a c o b i 矩阵，不 但迭 代每一 步的工作 量很大 ，而且有时无法 求得 J a c o b i 矩阵 当根据基质势 h 和试验数据来线性内插 和 C 时 或 J a c o b i 矩阵很难用解析表达式来 表示 当K h o 和 JII 与 h的关系式很复杂时 。 为确 保收敛速度并减少工作量，本文采用 P i c a r d迭代 法。具体迭代方法和步骤如下 1 设定一个迭代收敛条件l JII 。 JII 一≤ ，其 中上标 k代表第 k时步 ，下标 ， l 和 ， 分别代 表第 ， l 和 ， 次迭代得到 的值 ， 为前后两次迭代计 算结果间的误差允许值。 2 根据 和 ， ‘时刻的渗流场 h H 和 h ‘ ，线 性外推待求时刻 的初始渗流场，即h I I h C k 0时 或 h h k A t k ／ A t 卜’ J Il ‘ 一JI l 卜 O时 ，这里 At ‘ t k 一-- t k 。 3 由t 时刻的渗流场h f h ‘ h ，k 1 ／ 2 确 定 出各 单元 的渗透 系数 ， 和容 水度 C 后，计算有限元计算格式中的 ， ， ， 和 D ，用 Cr o u t分解 法求解 式 2 ，得 到新 的渗流 场 h 甜 。 ． 4 将h kf l l和h k， l 进行比 较， 若l h k， l1 ．j _ h l ≤ ，结束迭代进入下一个时步。否则，继续重复 进行第 3 和 4 步直至满足收敛条件， 即可求得降雨 入渗过程中边坡非饱和． 饱和渗流场、非饱和带暂态 饱和区 、暂态水压力以及基质吸力。 3 岩质边坡降雨入渗的规律研究 3 ． 1 某露天矿岩质边坡降雨入渗的数值模拟 某露天采矿边坡由花岗岩 ，变质岩 S u ，矿 层 F e 组成 如图 1 。地下水埋深 6 0 m 左右 。非饱 和带位于强风化岩体内【 l l 】 。 非饱和带渗透 系数 K3 ． 9 41 0 - 2 I h l 叫 m ， d ， 饱和渗透系数 3 ．9 4 X1 0 m， d ， 容水度 C J II 0 ． 0 1 l 9 Ih l - o 。6 ’ ；在饱和带中花 岗岩 r 渗 透系数 5 ． 5 9 1 0 _ 。 m ／ d ，变质岩 S u 渗透系数 { 醛 懈 水平距离， m 图 l 某露天矿边破状况 F i g ． 1 S i t u a t i o n o f t h e s l o p e i n a c e r t a i n o p e n p i t mi n e K8 ． 2 l l 0 - 2 m／ d ，矿层 F e 渗透系数 K4 ． 0 6 1 0 _ ‘ m／ d 。依据所选参数及初始条件 、边界条 件， 模拟某矿边坡在 7 d 最大降雨为 6 7 6 ． 1 m m的情 况下，非饱和渗流场的变化情况 如图 2 。 3 ． 2 岩质边坡降雨入渗规律的研究 模拟结果 反映了边坡 降雨入渗状态下非饱和． 饱和渗流场的变化规律。 1 边坡非饱和区，降雨历时 l d以后 ，首先在 饱和区逸出点以上 附近出现暂态饱和区。随降雨历 时的增加，暂态饱和区范围沿坡面线附近区域不断 扩展 、延伸 。至降雨历时 4 d后 ，在整个坡面线 附 近区域形成 了一个 由暂态饱和 区构成 的 “ 饱和壳 ” ， 并随着降雨过程的延续 向边坡 内部扩展，影响范围 逐渐扩大 如图 2 。 2 随着降雨不断入渗，暂态饱和区内出现正 的水压力，此 即为通常所说的暂态水压力 。暂态水 压力在坡面线附近数值较小，在暂态饱和 区下缘与 非饱和区接触的地带为零。最大值 出现在非饱和区 中部偏上部位 如图 2中 断面 。在整个 7d的 降雨过程中，暂态水压力随降雨历时的增加而逐渐 增大，其最大值为 2 5 ． 6 6 l P a 相当于 2 ． 5 6 6 m 压力 水头 。暂态水压力随降雨时间的变化如图 3 。 3 随着降雨入渗边坡，非饱和区上部的含水 量逐渐增大，基质吸力逐渐降低 。非饱和区最大基 质吸力随降雨时间的变化如图4 。 4 在 7d连续降雨过程 中， 暂态饱和区未延伸 到潜水面，因而饱和区渗流场的变化不大，潜水面 升高很小，这与 以往运用饱和渗流模型模拟降雨入 渗的结果有很大差别 。 维普资讯 6 2 8 岩石 力学与工程学 报 2 0 0 3盔 3 2 o ．o o Fi g ． 2 3 0 ． 0 芒 妄 o o 出 磁l 0 ． 0 缸 O ． O 图2 降雨过程中，某矿岩体边坡饱和． 非饱和渗流场 过程中 断面基质吸力的麦化 0一为基质吸力 0一为水头压力 同 I 降雨 l d后 降雨 2 d后 降雨 3 d后 降雨 4d后 降雨 5 d后 降雨 6d后 降雨 7d后 s e e p a g e fi e l d i n t h e p r o c e s s o f r a i n f a l l i n a r o c k s l o p e i n a c e r t a i n o p e n p i t mi n e ／ ． 一 一。 ， 一， ．／。 J 一 ／ ．／ 。 O l 2 3 4 S 6 7 降雨历 时， d 图 3 降雨过程中，暂态水压力的变化 Fi g ． 3 Va r i a t i o n o f t r a n s i e n t wa t e r p r e ssu r ei nt h e p r oc e s s o fr a i n f a l l 4 0．0 3 0O 2 0．0 蜒 糊 1 0 0 0 ．0 ’ ＼ ‘＼ 0 l 2 3 4 5 6 7 降雨历 时， d 图 4 降雨过程中，基质吸力的变化 Fi g ． 4 Va r i ati o no f ma t r i xs u c t i o ni nt h ep r o c e ss o f r a i n f a l l 、 一 ￡ 维普资讯 第 2 2 卷第 4 期 戚国庆等．岩质边坡降雨入渗过程的数值模拟 6 2 9 关于在学术论文 中规范关键词选择 的规定 试行 一 、发表在 中国科协系统学术期刊 中所有学术论文，必须在摘要后列出不少于 4个关键词。从技术角 度考虑 ，没有关键词 的论文应列入非学术论文类。 二、这些关键词按以下顺序选择 第一个关键词列出该文主要工作或内容所属二级学科名称。学科体系采用国家技术监督局发布的学科 分类与代码 国际 GB ／ T 1 3 7 4 5 ． 9 2 。 第二个关键词列出该文研究得到的成果名称或文内若干个成果的总类别名称 。 第三个关键词列出该文在得到上述成果或结论时采用的科学研究方法 的具体名称。对于综述和评述性 学术论文等 ，此位置分别写 “ 综述 ”或 “ 评论 ”等 。对科学研究方法的研究论文 ，此 处不写被研 究的方法 名称 ，而 写所应用的方法名称 。前者出现于第二关键词的位置 。 第 四个关键词列出在前三个关键词中没有出现的，但 被该文作为主要研究对象的事或物质的名称 ，或 者在题 目中出现的作者认为重要 的名词 。 如有需要，第五、第六个关键词等列出作者认为有利于检索和文献利用的其他关键词。 中国科协学会学术部 2 0 0 2年 8 月 2 2日 维普资讯