（岩土工程专业论文）降雨条件下的边坡渗流数值模拟及稳定性分析.pdf

降雨条件下的边坡渗流数值模拟及稳定性分析 摘要 降雨入渗时边坡的稳定性问题不仅牵涉到饱和非饱和的渗流规律，而且涉 及含水率变化对非饱和土的强度参数的影响作用。本文针对不同的土体性质，分 析了不同降雨条件对土坡渗流及稳定的影响，运用饱和非饱和渗流理论及降雨 入渗理论，分析了渗流场的变化过程，根据计算出来的瞬态渗流场利用非饱和土 抗剪强度理论，对土坡的瞬态稳定性进行了分析。研究的主要内容有： 1 介绍了降雨条件下的饱和非饱和渗流基本理论，介绍了非饱和渗流的控 制方程，给出了降雨条件下的边界条件。 2 基于g a l e r k i n 加权余量法得出了饱和非饱和渗流计算有限元列式。针对 不同的土体性质，利用简单的土质边坡有限元模型运用s e e p w 渗流计算程序计 算在给定降雨条件下的孔隙水压力的分布，探讨降雨强度、降雨持时、雨型、土 壤渗透系数以及边坡坡度等参数的变化对渗流场的影响。 3 介绍了饱和非饱和渗流条件下边坡稳定分析理论。运用g e o s l o p e ( s l o p e w 模块) 软件，将降雨入渗作为土坡稳定分析的重要因素，将渗流场 的孔隙水压力引入基于非饱和土抗剪强度理论的抗滑稳定计算中。定量的得出在 各降雨因素及边坡坡度对非饱和土坡稳定性影响的程度；探讨各参数变化对边坡 安全系数的影响，并对边坡的稳定性进行了评价。分析结果表明：在甘肃大部分 地区边坡坡度比降雨对边坡的稳定性影响更大；在甘肃大部分地区在可能发生的 降雨条件下，坡角大于6 了的边坡就要考虑对边坡进行支护。 关键词：饱和非饱和渗流；降雨入渗；孔隙水压力；有限元；边坡稳定 硕一j ：学位论文 a b s t r a c t t h e s t a b i l i t yp r o b l e m o f s l o p e s u n d e rr a i n f a l li n v o l v e sn o t0 n l yt h e s a t u r a t e d u n s a t u r a t e ds e e p a g el a w ，b u ta l s ot h ee f f e c t0 fw a t e rc o n t e c tc h a n g e so nt h e i n t e n s i v ep a r a m e t e r0 fu n s a t u r a t e ds o i l f o rd i f f e r e n tt h es o i lc h a r a c t e r i s t i c s ， t h e s e e p a g e a n d s t a b i l i t y 0 f s l o p e u n d e rr a i n f a l l i s a n a l y z e d b a s e d 0 nt h e s a t u r a t e d - u n s a t u r a t e ds e e p a g et h e o r ya n dr a i n f a l li n f i l t r a t i o nt h e o r y ，t h ec h a n g eo f s e e p a g ef i e l d si sa n a l y z e d o nt h eb a s i so fc a i c u l a t e dt r a n s i e n ts e e p a g ef i e l da n dt h e s h e a r s t r e n g t ht h e o r y o fu n s a t u f a t e ds o i l ，t h et r a n s i e n ts t a b i l i t y0 fs l o p ei s a 肋l y z e d t h em a i nc o n c l u s i o n sa r ea sf o l l o w s ： 1t h es a t u t a t e d u n s a t u r a t e ds e e p a g et h e o r yu n d e rr a i n f a l lc o n d i t i o na n dt h e g o v e r n i n ge q u a t i o no fu n s a t u r a t e ds e e p a g ea r ei n t r o d u c e d t h eb o u n d a r yc o n d i t i o n s a r ea t t a i n e d 2t h ef e mf b r m a t0 fs a t u r a t e d - u n s a t u r a t e ds e e p a g ed i f f c r e n t i a le q u a t i o ni s d r a w e db a s e do ng a l e r k i n sm e t h o do fw e i g h t e dt e s i d u a l f o rt h ed i f f e r e n tp r o p e n y s o i l ，p o r ew a t e rp r e s s u r ef i e l du n d e rd i f f e r e n tr a i n f a l lc o n d i t i o n sa r ec a l c u l a t e db y a d o p t i n g f e mm o d e lo f s i m p l e e a r t h s l o p e u s i n gs e e p wn o wc a l c u l a t i o n p r o c e d u r e s t h e i m p a c t s o f s e e p a g e f i e l du n d e rt h ev a r i e sw i t hr a i n f a l l i n t e n s i t y ，d u r a t i o n0 fr a i n f a l l ，r a i nt y p e ，s o i lp e r m e a b i l i t yc o e f f i c i e n t sa n ds l o p ea n g l e a r ed i s c u s s e d 3t h eb a s a lt h e o r y0 fs a t u r a t e d u n s a t u r a t e ds l o p e s s t a l b i l i t y i si n t r o d u c e d r e g a r d i n gr a i n f a l ls e e p i n ga sa ni n l p o r t a n tf a c t o r0 fs o i ls 1 0 p es t a b i l i t ya n a l y s i s ， t h e d i s t r i b u t i o n0 fp o r ew a t e rp r e s s u r eo fs e e p a g ef i e l di sc a l c u l a t e da n di n t r o d u c e di n t o t h ef o r m u l ao fs a f e t yf a c t o ru s e i n gg e o s l o p e ( s l o p e wm o d u l e ) s o f e w a r e t h e i m p a c t s 0 fn o n - s a t u r a t e ds 1 0 p es t a b i l i t ya r e q u a n t i t a t i v e l y d r a w e du n d e rt h e i n n u e n c eo fr a i n f a l la n ds l o p eg r a d i e n t ；t h ev a f i e so fp a r a m e t e r sa r ed i s c u s s e do nt h e i m p a c to fs l o p es a f e t ya n dt h es t a b i l i t yo ft h es l o p ei se v a l u a t e d t h er e s u l t ss h o w t h a ti nm o s tp a r t s0 fg a n s us 1 0 p eg r a d i e n th a v eag r e a t e ri n l p a c to nt h es t a b i l i t yo f s l o p et h a nr a i n f a l l ；u n d e rt h ep o s s i b l er a i n f a l lc o n d i t i o n s ，t h es l o p e so fm o r et h a n 6 5 0s l o p eg r a d i e n ta r ec o n s i d e r e ds u p p o r t i n g k e y w o r d s ： s a t u r a t e d - u n s a t u r a t e ds e e p a g e ， r a i ni n f i l t r a t i o n ， p o r ew a t e fp r e s s u r e ， f e m ，s l o p es t a b i i i t y i i 兰州理工大学 学位论文原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所 取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任 何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡 献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的 法律后果由本人承担。 作者签名：励呜日期：的年f 月卢日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，即： 学校有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许 论文被查阅和借阅。本人授权兰州理工大学可以将本学位论文的全部或部 分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段 保存和汇编本学位论文。同时授权中国科学技术信息研究所将本学位论文 收录到中国学位论文全文数据库，并通过网络向社会公众提供信息服务。 作者签名： 导师签名： 瞄怕岛 勿匆哆 日期：岛年 月户日 日期：略年多月移日 硕十学位论文 1 1 课题背景及意义 第1 章绪论 我国是一个多山、多高原的国家，斜坡地质灾害问题比较严重。不论是自然 斜坡的滑动、崩塌，还是人类工程活动形成的斜坡失稳，都对经济建设和人民生 命财产造成巨大损失。据统计，我国的滑坡等地质灾害主要分布在2 4 个省、市、 自治区，其中以四川、湖北、广西、云南、甘肃等省滑坡发生频率较高。表1 1 中总结了我国部分近代滑坡发生的诱因【。 表1 1 我国部分近代滑坡与降雨的关系统计表 在世界其它地区，每年也都发生为数众多的非饱和残积土斜坡的破坏，如巴 降雨条件下的边坡渗流数值模拟及稳定性分析 西，南非和远东地区等，对这些破坏原因的研究已有多年。这些斜坡发生破坏的 主要原因是降雨【2 7 】。如香港地区，在1 9 7 2 年6 月1 6 1 8 日期间的一次暴雨引发 了滑坡泥石流，由于滑坡泥石流造成的伤亡总数达2 5 0 人，成为香港历史上滑坡 泥石流灾害最惨痛的一页【引。又如2 0 0 6 年2 月1 7 日在菲律宾的中东部莱特省圣 伯纳德镇以及周边地区发生特大泥石流灾害，截至到搜救工作结束时为止，有近 3 0 0 0 人失踪。引起这次泥石流灾害的重要原因之一就是连续数日的暴雨。在我 国，由于降雨引起的滑坡、泥石流灾害发生也十分频繁，灾害造成的损失也十分 巨大。可以看出，滑坡受降雨影响密切，泥石流直接受控于降雨。 从表1 1 中可以看出，滑坡、泥石流很多都出现在非饱和土区域内。在非饱 和土内，由于土体中气相的存在，导致非饱和土中产生了负的孔隙水压力，即基 质吸力【引，降雨诱发非饱和土斜坡失稳的原因是雨水入渗引起负孔隙水压力的升 高直至消失，导致非饱和土的抗剪强度降低【9 1 。随着降雨因素改变，非饱和土内 的负孔隙水压力的变化将直接影响非饱和土的抗剪强度变化，从而导致了非饱和 土斜坡安全性的变化。因此，研究降雨因素对非饱和土内孔隙水压力的变化，是 研究降雨对斜坡稳定性影响的基础。 降雨引起非饱和土斜坡失稳的原因是，当降水渗入到非饱和土体内时非饱和 土土体内的负孔隙水压力开始升高。负孔隙水压的升高导致了土体的抗剪强度的 降低。当滑裂面上的抗剪强度减小到比下滑力还小的时候，斜坡就会发生失稳破 坏1 9 1 。深入研究降雨引起斜坡失稳的规律并建立定量的分析模型对于滑坡、泥石 流灾害的预测和预防有重要的指导意义。 本文从降雨的雨型对非饱和土斜坡内的负孔隙水压力的影响入手，对在不同 的雨型、初始地质条件下，非饱和土斜坡内孑l 隙水压力的变化规律进行研究。在 此基础上得到对孔隙水压力变化有显著影响的降雨因素和初始的地质条件，为滑 坡、泥石流的预报和防治工作提供一些有参考价值的依据。 1 2 降雨条件下土坡稳定分析的研究现状 国内外学者在对降雨条件下的土坡的稳定性分析上做了很多研究。代表性学 者的主要方法如下： f l e d l u n d 等强调了负孔隙水压力对土坡稳定性的影响，并重新对土坡稳定分 析的安全系数计算公式进行了推导，以便把正负孔隙水压力都包括进去。他通过 对香港地区某一边坡的算例进行了暴雨情况下的渗流和边坡稳定分析，在分析时 采用稳态条件下的孔隙水压力分布作为暴雨期间瞬态分析的初始条件，用 g a l e r k i n 有限元法对渗流场进行了模拟，通过渗流分析得出暴雨在开始后不同时 刻的孔隙水压力分布情况，从而分析了土坡在不同时刻的稳定性。 l u m b 【1 0 】率先研究了香港地区降雨和滑坡的关系，他采用了一种简化的一维 2 硕一l ：学位论文 垂直入渗模型，根据抗剪强度和饱和度的经验关系研究了地质条件和降雨特征对 斜坡稳定性的影响。分析出香港地区斜坡的稳定性不仅是由土的入渗能力来控制 的，还取决于降雨的类型和降雨量的大小，但在研究中未考虑水平渗流分量的影 响，对入渗量的假定与实际情况有出入。 s a m m o r i 和t s n b o y a n m a 【1 1 】采用g a l e r k i n 有限元法模拟了恒定降雨强度条件 下的边坡暂态渗流过程，并对边坡稳定性进行了参数研究，分析出斜坡稳定的影 响因素有土坡的渗透状况，斜坡的长度以及土层的厚度等等。 a l o n s o 【1 2 】等在考虑土体类型、降雨历时、降雨强度和土体渗透性能的基础 上采用考虑空气压力变化涡合的分析方法计算了二维非饱和渗流场，进行了土坡 二维非饱和渗流和极限平衡法的联合分析，研究了降雨入渗对边坡稳定性影响。 s h i m a d a i ”】等用g a l e r k i n 有限元法模拟了二维非饱和渗流，并采用刚体弹簧 模型进行了斜坡稳定性数值分析。得出不同的降雨强度和基质吸力的变化对斜坡 安全系数的影响情况。 s u n 【1 4 】等发展了应力和两相流涡合理论，分析了降雨条件下土坡的浅层破 坏。在对实际边坡分析时，分别研究了在考虑和不考虑空气压力影响条件下的土 坡渗流情况，并得出考虑空气压力的分析是合理的结论。 n g 和s h i 【”】等对香港地区利用有限元方法对降雨情形和初始条件的暂态渗 流对斜坡稳定性进行了研究，分析出降雨特征、水文地质条件和坡面防渗情况对 斜坡安全系数的影响。 吴宏伟、陈守义【1 6 】等针对香港地区降雨的特点和地下水位分布的特点，用 河海大学硕士学位论文法深入系统地研究了降雨强度、降雨历时、土体渗透性及 阻水层埋藏条件等对暂态渗流场和斜坡稳定性的影响，得到了一些有益的结论， 但在研究中未考虑前期降雨的影响。 姚海林、陈守义等【r 7 l 对非饱和膨胀土的暂态渗流场进行了参数分析研究， 分析了降雨和土体渗透性对土坡稳定的影响，在分析裂隙的作用时，分析了裂隙 的开裂对土坡的孔隙水压力和体积含水率的影响。 史弘鹤、王殿春【1 8 】在研究雨水入渗和蒸发过程的基础上，提出了一种考虑 土体气水的涡合暂态渗流场模型，用此模型分析了降雨入渗引起的暂态渗流场 变化，并分析了渗流场变化对土坡稳定的影响。 陈守义等【1 9 l 提出了一种考虑入渗和蒸发影响的土坡稳定性分析方法。采用 数值法求解任意给定入渗和蒸发边界条件下斜坡土体的瞬态含水率分布和斜坡 瞬态抗剪强度参数分布，在此基础上再使用常规土坡稳定分析方法计算土坡瞬态 安全因数。 高润德等【2 0 】针对南水北调中线工程河南段某一黄土边坡，研究了雨水入渗 作用下的土体渗透性和坡顶垂直裂缝对边坡稳定性的影响。 3 降雨条件下的边坡渗流数值模拟及稳定性分析 其他一些学者也对饱和非饱和渗流计算，降雨条件下土坡的渗流稳定进行 了探讨【2 1 - 2 5 1 。 1 3 本文主要研究内容及研究方法 在对降雨入渗的分析研究上，本文针对某个特定的区域对不同的土体性质 ( 比如具有不同的饱和渗透系数) 在降雨条件下的土坡渗流及稳定分析研究，具体 工作如下： 1 ) 介绍了降雨条件下的饱和非饱和渗流基本理论，介绍了非饱和渗流的控 制方程，给出了降雨条件下的边界条件。 2 ) 本文简要的叙述了非饱和土中二维渗流控制方程的g a l e r k i n 有限元解 法，提出了利用g e o s l o p e 软件的s e e p w 进行二维瞬态渗流计算并进行了建 模分析。 3 ) 根据甘肃大部分地区的土水特性资料和降雨资料，对简单的二维均质各 向异性的土质边坡降雨入渗条件下的饱和非饱和非稳定渗流进行分析，并且分 析了降雨、土体渗透系数及边坡坡度对孔隙水压力的影响。 4 ) 在渗流场分析完成的基础上，对不同坡度的边坡在不同降雨条件下的稳 定性进行比较分析。 4 硕士学位论文 第2 章降雨条件下饱和非饱和渗流基本理论 2 1 饱和非饱和渗流方程【2 6 3 4 1 2 1 1 饱和渗流基本微分方程 由文献可知，饱和条件下，考虑土体和水压缩性的非稳定渗流微分方程式为： 击 ，芸) + 喜 ，等) + 鲁 ：警) - 墨詈 c2 - 1 ) 式中七，七，七：为三个主轴方向上的渗透系数，j l 为渗流区域中各点的总水头，s ，为 单位贮存量( 量纲为【1 l 】) ，即单位体积的饱和土体，当下降1 个单位水头时，由于土体压 缩和水的膨胀所释放出来的贮存水量，它与多孔介质和流体的弹性特征有关。 2 1 2 非饱和土的达西定律 达西定律在层流条件下不仅适用于饱和土中水的流动，对于非饱和土中水的 流动也同样适用。在非饱和土中，渗透系数一般不能假定为常数，它是变量，是 非饱和土含水率或基质吸力的函数。而在特定的含水率条件下，对于作用于非饱 和土上的不同水力梯度，渗透系数七。是常数，也就是说，通过非饱和土的水流 速与水力梯度呈线性比例关系，而渗透系数是一常数，这与饱和土中的情况相同。 这说明达西定律同样适用于非饱和土。 非饱和土的达西定律可以表示为： y 。一七。筌 ( 2 2 )y 一七w 【z 。z ) 式中：，为水的流速，七。为水相的渗透系数，譬为y 方向的水力梯度。 对于饱和渗流分析，土中每一点的渗透系数为一不随时间变化的常数。但在 非饱和土中渗透系数常常是各点互不相同的。非饱和土的渗透系数空间变化的主 要原因是孔隙流体( 即孔隙水) 的非均质体积分布。虽然土颗粒分布可能是均质 的，但孔隙流体的体积分布可能由于基质吸力的空间变化而显示为非均质性能。 非饱和土中的渗透性状况可以分为非均质各向同性，非均质各向异性和渗透性随 空间连续变化三类，这种分类以渗透性变化模式为依据。 2 1 3 非饱和渗流基本微分方程 将达西定律和质量守恒定律结合起来可以推导出二维问题非饱和土体水流 运动的基本方程： 詈一书卦讣p ) 詈】 3 ， 5 降雨条件下的边坡渗流数值模拟及稳定性分析 式中：h 为总水势，h = h + z ，h 为压力水头( 饱和区为正，非饱和区为负) ，z 为位置水头，秒为体积含水率。在饱和区渗透系数为饱和值，与h 无关，在非饱 和区，渗透系数是压力水头h 的函数；七，七，为x ，y 方向的渗透系数。 对于饱和土体来说，土体孔隙为水充满，此时含水率不再变化，即：竺；0 ， 由式( 2 3 ) 可以得出饱和土体水流动的控制方程： 坤以) 卦扑缸) 针。 仁4 ， 此时x ，y 方向的渗透系数为常数，饱和渗流可看成是非饱和渗流的特例。 由于体积含水率是压力水头h 的函数，而h 又是时间的函数，所以有： 旦卫堕 一i 一一 i 一1l at8hat 、 1 一般情况下，土水特征曲线表示的是含水率和基质吸力的函数，则有： 一生上塑 ( 2 6 ) 一一一一 z - nl a ( “。一口。)p 。ga 、7 式中：p 。为水的密度，g 为重力加速度。 定义比水容重研：为体积含水率p 对基质吸力 。一“w ) 的偏导数的负值，即 朋：w 。一了与，则册：”表示为土水特征曲线斜率的绝对值。故： d 怛口一“。j 等= p ，肿：w - c 伪) ( 2 7 ) 式( 2 - 7 ) 表示为容水度的表达式。容水度的定义是压力水头减小一个单位时， 从单位体积土壤水中释放出来的水体体积。 由于h = h + z 且丝o ，则塑；丝。 综合上述可得饱和非饱和渗流的微分方程为： 丢卜o ) 警卜号卜仍) 等 2 p 。胛z ”詈 c 2 固 对于一般各向异性土体，主渗透方向与x 轴成一定的角度，则式( 2 8 ) 变为： 昙卜仍) 警+ b ) 警卜专卜仁) 等+ b 伪) 等卜“胛z ”警 ( 2 - 9 ) 式中t 、七，为沿x 、y 方向的渗透系数，k = 七lc o s 2 口+ 七2s i n 2 口， 七 一七1s i i l 2 口+ 七2 c o s 2a ，b = 七”= 瓴一七2 ) c o s 口s i l l 口，七l 为最大主渗透系数，七2 为最 小主渗透系数，口为最大主渗透方向与x 轴的夹角。 2 1 4 定解条件 在饱和- 非饱和渗流分析中，定解条件的方程式为式( 2 1 0 ) 。式( 2 - 1 0 ) 中l 为 水头分布规律已知的边界，日。为边界水头；r 2 为流量已知的边界，吼为单位时间 边界法向流量，c o s 仿，工) ，c o s 仁，y ) 等于边界面外法线方向的方向余弦；r 3 为渗出 6 硕士学位论文 面边界，z 为渗出面节点坐标，可归为第一类边界条件；日。为初始时刻的水头值， 称为初始条件。 一 | i l ，) ，f ) 一日1 0 ，) ，t ) ，g ，y ) r 1 t 芸听力嵋考c o 昨川力如) ，) r 2( 2 - 1 。) j i i o ，) ，f ) az o ，y ，f ) ，o ，y ) r 3 o ，) ，f o ) 一日0 0 ，) ，f o ) 在边界节点处已知水头为第一类边界条件，而己知通过边界的流量称为第二 类边界条件正的节点流量表示节点处有入渗，负节点流量表示该节点处有蒸发 或蒸腾。当通过边界的流量为零时，表示不透水边界。 2 2 降雨入渗基本理论 2 2 1 概述 研究降雨入渗对土坡稳定有着特别的意义。对于岩土质边坡，降雨会降低岩 土体的抗剪强度，抬高地下水位使得孔隙水压力升高，另外长时间高强度降雨会 使得稳定地下水位以上区域出现暂态饱和区，则相应区域会出现暂态孔隙水压力 升高的情况。由于许多滑坡事故并不是出现在地下水位抬升到相应高程后才出现 的，所以对降雨入渗引起的暂态水荷载需要引起足够的重视。 早在上个世纪初期，人们就开始研究降雨入渗问题随着科学技术的发展， 入渗分析已经从一维发展到了二维和三维，但是应用的比较广泛的还是一维的降 雨入渗模型，因为一维入渗的分析结果基本上反映了入渗的规律。本节阐述了降 雨入渗的机理，降雨入渗影响土坡稳定的机理以及降雨条件下边界条件的确定。 降雨入渗的过程是一个随时间和空间变化十分强烈的动态过程。雨水入渗量 受降雨强度，降雨历时，土体物理性质，入渗面坡度，植被情况等诸多因素的影 响。当降雨强度超过土体的入渗能力时，将在表面产生径流或积水，向土体内部 形成不断扩大的饱和区，此时土体中水分的运动是饱和非饱和运动。这种模型 可称为“积水模型”当降雨强度低于土体的入渗能力时，入渗过程受供水能力的 控制，可称之为“降水模型”一般降雨过程可分两个阶段：开始时，地表的含水 率梯度很大，入渗率也很高，一般而言大于降雨强度，这一阶段称为通量控制阶 段，这一阶段是无压入渗或自由入渗。随着入渗的进行，含水率梯度不断减小， 入渗率也不断降低，当小于降雨强度时，开始形成地表径流或积水，转为剖面控 制阶段，为有压入渗。因此实际的入渗量一般取决于土体的初始含水率，降雨强 度和持时以及表面径流量。 根据c o l e m a n 和b o d m a n 的研究，当均质土体地表有积水入渗时，典型含水 率分布剖面可分为四个区：饱和区，过渡区，传导区，湿润区。湿润区的前缘称 7 降雨条件下的边坡渗流数值模拟及稳定性分析 为湿润锋。其中，各个部分的特征如下： ( 1 ) 饱和区：孔隙被水充满或处于饱和状态，该区域通常只有几毫米厚，这与 积水时间有关； ( 2 ) 过渡区：该区域含水率随深度增加迅速下降，一般向下延伸几厘米厚； ( 3 ) 传导区：该区域含水率随深度增加变化很小，通常传导区是一段较厚的高 含水率非饱和区； j ( 4 ) 湿润区：该区域含水率随深度增加从传导区较高含水率值急剧下降到接 近初始含水率值； ( 5 ) 湿润锋：湿润锋在干土和湿土之间形成一个陡水力梯度的锋面。 单位时间内通过地表单位面积入渗到土体中的水量称为入渗率。一般来说， 实际降雨入渗过程中，一开始地表含水率较低，一旦有水分入渗时，就形成很高 的吸力梯度，这时入渗率很高，然后随着地表吸力梯度不断减小，入渗率逐渐降 低至一个稳定值，即饱和水力传导度。 研究入渗问题主要是分析土体剖面中含水率分布随时间变化及湿润锋前移 的规律。 土壤水动力学p j 详细介绍了这方面的解析解和半解析解的模型及解 法 本文主要是用有限元数值解法求解非稳定的饱和非饱和渗流问题来研究降 雨入渗过程。 2 2 2 降雨入渗的定解条件的处理 1 降雨入渗下初始条件的处理 进行有地表入渗的饱和非饱和非稳定渗流场模拟时，首先要给定计算区域 的初始水头场。边坡的边界条件较难给定，特别是土坡两侧的状态很难给定。在 缺乏实测值的情况下，非饱和区的初始水头场很难准确地确定。因此需要一定的 假设，即此区域不因降雨而和相邻区域发生水交换。这意味着只由初始条件决定。 雨水的入渗对初始的渗流场比较敏感，因此本文采用的方法是利用求解稳定状态 下的渗流场来获得初始的水头值。 2 降雨入渗的边界条件的处理 图2 1 给出了土质边坡在降雨入渗条件下的边界条件。其中a b 、g h 段是水 头边界，h 已知，在此边界面上，j i i z + 盟为常数；b c 、g f 段为未知水头边界， y 一 假设流量q = o ；a h 段为不透水边界，此边界上娑，o 咖 对坡面边界的边界条件按以下方法处理： 设降雨强度为r ( t ) ，在降雨的初期，由于雨强小于土体的入渗能力，在地表 8 硕士学位论文 将不会形成径流或积水，此时，地表入渗率与降雨强度相等。 当降雨强度大于土体的入渗能力时，超出入渗能力的降雨则转化为地表积水 或径流，而此时地表的总水头维持一个不变值： h = z + d 。( 2 1 1 ) 式中：z 为位置水头，d 为积水深度，在斜坡面上可认为积水深度为0 。 c b a 图2 1 土质边坡在降雨入渗条件下边界条件的区段划分 9 f h 降雨条件下的边坡渗流数值模拟及稳定性分析 第3 章降雨入渗条件下饱和非饱和渗流分析 3 1 降雨入渗条件下饱和非饱和渗流有限元原理【姗9 】 降雨入渗过程实质上是入渗水分在非饱和区运动的过程，而降雨入渗边界又 是一个流量边界，只是这个流量并不是不变的，在计算过程中需要根据含水率调 整入渗流量，从而实现对降雨入渗问题的数值模拟。 3 1 1 饱和非饱和渗流的定解问题 、以二维问题为例，饱和非饱和渗流问题所应满足的控制方程和边界条件为： - * 蝶膀卜帆讣；詈一 a ， j lg ，) ，f ) 日。g ，y ，flg ，y ) s l ( 3 1 1 b ) 七假( 钆) 罢c o s 伍，x ) + 七眇( 钆) 票c o s 仿，y ) = 吼g ，y ，f ) o ，) ，) s ： ( 3 1 1c ) 批 卵 。 j i ig ，y ，t ) = z g ，) ，fl ( z ，y ) s 3 ( 3 1 1 d ) b ，y ，f o ) 一日ob ，y ，f oj ( 3 1 1 c ) 式中：而为水头分布规律已知的边界，甄为边界水头，式( 3 1 1 b ) 称为第 二类边界条件；j 2 为流量情况已知的边界，吼为单位时间边界法向流量，0 0 s 仿，石) ， s 伍，y ) 为边界面外法线方向余弦，式( 3 1 1 c ) 称为第二类边界条件；屯为渗出 面边界，z 为渗出面节点坐标，式( 3 1 1 d ) 可以归为第一类边界条件；日。为初始 渗流场水头式( 3 1 1 e ) 为初始条件。 3 1 2 有限元列式的推导( g a l e r 埘n 法) 在直接求解压力水头函数j l i ( x ，y ，t ) 有困难的情况下，可以通过选取试函数 j l ( x ，y t ) ( 要求它严格满足边界条件) 作为饱和非饱和渗流偏微分方程边值问题的 近似解。由于它是一个近似解，代入有关方程后的值不会为0 ，把这个值记为 r ( r 一0 ) 。如使r 趋近于零，则这个水头近似函数j l i ( x ，y ，t ) 就为可能逼近微分方程 精确解h ( x ，y t ) 的解。我们要求满足 1 i 黜幽= q ( 3 1 2 ) 式中：w 为权函数。这个通过余项与权函数的乘积在区域内的积分为零而推求原 问题的方法，就称为加权余量法。若权函数为特定的函数，例如取形函数n ，则 称g a l e r k i n 法 本节以有限元法中最简单和常用的3 结点三角形单元为例，对二维域进行离 l o 坝士学位论文 散，取其中一单元，其结点编码为i ，j ，m ，见图3 1 1 所示。 假设三角形单元i j m 内的水头模式为一次多项式： 一口l + 口2 工+ 口3 y ( 3 1 3 ) 将i ，j ，m 节点坐标代入方程( 3 1 3 ) 得到各节点的水头值： i h i 口l + 4 2 而+ 口3 ) ，i h j 一口l + 口2 z ，+ 口3 y ， ( 3 1 4 ) ih m - 口l + 口2 工一+ 口3 ) ，_ 吩功 图3 1 - 1三角形单兀离散及单个三角形i j m 从方程( 3 1 4 ) 可以解出待定常数口l ，口2 ，口3 ，再将口l ，口2 ，口3 代到方程 ( 3 1 3 ) 即可将该三角形单元内的水头函数表示为三角形三个节点水头的函数： j i i 。【j ，肼】 日 。【】( 日p ( 3 1 5 ) 【日，j 式中，札( f ，j ，历) 为单元的形函数，其值为： fm 一( 4 l + 岛x + c j y ) 2 ，- ( 口，+ 6 ，+ c ，y ) 2 ( 3 1 6 ) 卜庸i ( 口。+ 6 。工+ c 庸y ) 2 该形函数具有如下性质： m ) ，沪驴忙：二 ( 3 1 7 ) 式中：系数口j ，饥，c ，( i = 1 ，2 ，3 ) 和三角形单元的面积都是坐标x ，y 的函数。 4 j 一工y 一一工_ ) ，饥一) ，一) ，一，ci - 工一石， ( 3 ：1 8 ) 其它系数按照i ，j ，m 的次序轮换排列。三角形单元的面积为： 一丢巨主；三i 。丢g c ，一c t 6 ，) ； c3 1 9 ) 将水头函数代入控制方程中，得到残余项是： r b ，y ) _ 丢( 七孵p w ) 芸) + 号( 七町饥) 号) 一几印；詈加 ( 3 加) 1 1 寻求权函数西j 使 p b ，y 净，k y 胁j o ( j = 1 ，2 ，3 ) 令弩函数垂，等于插值函数，有： ! 陋( k 帆) 罢产+ 专卜町帆学) 一“朋；詈卜仁，y m 一。 此处，j = 1 ，2 ，3 对应单元的i ，j ，m 结点( 以下也同) ：考虑到： ，丢( k 芸) 一鼬七眦静k 誓芸 p ，昙( k 芸卜。麝( ，k 芸卜一p 精笔墨 二昙( ，七。芸卜a 正，七孵芸s 仿，x 拯 房( 坼号卜- 以k 考c o s 仿，y 迹 式中：s 为围绕区域a 的边界。 由式( 3 1 1 ) 中的流量边界条件： 妒凡罢哪g ，工) + p 凡芸c o s 仁，y k 。荫，凼 ( 3 1 1 2 ) ( 3 1 1 3 ) ( 3 1 1 4 ) ( 3 1 1 5 ) ( 3 1 1 6 ) ( 3 1 1 7 ) 式中：吼为单元外部垂直于单元边界方向的水流速率。 将式( 3 1 1 4 ) 一( 3 1 1 7 ) 代入式( 3 1 1 2 ) 中，得到： 啦警芸等爹卜一豇以嬲+ p 删，警 8 ) 再将式( 3 1 5 ) 代入式( 3 1 1 8 ) 得： 办盯警丢( ；| ；m 础) ) 盖等象砉札毗卅 。3 1 9 ) 一和，心+ p - 胂；以云( 荟m 吼o ) 卜。o 由于形函数m b ，y ) 只是坐标的函数，缮点水头也( f ) 只是时间的函数，上式 砉办孵誓警毗m 咿等等珥g 和 妒肛料，五掣，j 掣斗。 1 2 ( 3 1 2 0 ) 硕士学位论文 他】r k 】 8 k 日p + ，【】r a 【m 巫娑一， 】r 吼豳。o ( 3 1 2 1 ) oi| 。 【b 】r k 如k 式中： 臼p 嚣：i ，单元结点上的水头向量 1 日i a 。p 。g 柳； 。 吼为单元外部垂直于单元边界方向的水流速率。 对式( 3 1 2 1 ) 进行数值积分，可以得到更为简单的表达式 【d l h p + 陋l 忸j e ， f k ( 3 1 2 2 ) 式中： 【d 1 传递刚度矩阵( 一般有限元分析中称刚度矩阵) ，即陋】r k p k 乩槲槲元贮赢铘泅 位p 单元结点处水头对时间的导数矩阵，即型丝 ，) | 反映边界条件的流量矢量列阵，即r 【r 吼豳 进行整体分析，可以得到饱和非饱和渗流问题的求解方程： p 】 + 陋酗 ；仁) ( 3 1 2 3 ) 由于渗透系数是基质吸力的函数，亦即与结点处的水头值有关，故式( 3 1 2 3 ) 是非线性的，数值弥散和波动现象在流动问题的高度非线性系统中是经常遇到 的，尤其入渗处土体饱和度较低时，容易产生，这时，边界处单元尺寸缸，每 要取的很小，入渗初期的时间步长f 也要取的较小，在具体计算时，可能需要 经过数次试算、调整才能找到一个较好的方案 3 1 3 有限元方程在时间上的差分近似格式 对时间微分的一般处理方法是在某种特定时刻首先用有限元方法求得该点 的水头，然后在时间域内用有限差分方法向前推进。 降雨条件下的边坡渗流数值模拟及稳定性分析 用有限差分法来近似表达时间的微分，即： 堕。堡丝二丝! ( 3 1 2 4 ) 一。- 。一 j i 二， 砸f 考虑到在时段内水位的变化比较缓慢，h 值可近似地写成时段初和时段 末水头只和只+ 山的加权平均的形式，即： 日一仅h r + + q 一口归， ( 3 1 2 5 ) 式中：口为参数，0s 口s 1 。 将式( 3 1 2 4 ) 和式( 3 1 2 5 ) 代入求解方程( 3 1 2 3 ) 中，得到： p 舢+ ( 1 一口蚺叫声 - 仁 ( 3 m 6 ) 由于选取的a 值不同，加权平均值的含义也不同。从而就形成不同的有限元 方程在时间上的差分格式。 ( 1 ) 当口0 时，百a 日，即取时段初水头值日近似地代表时段的平均 水头值，在这种情况下，有： p 陋产萨 = 仁) ( 3 抛7 ) 式( 3 1 2 7 ) 称为显式差分格式。 ( 2 ) 当口= 1 时，厅一日& ，即取时段末水头值日& ，近似地作为时段垃的 平均水头值，在这种情况下，有： p + 陋】f ，一仁) m 8 ) 式( 3 1 2 8 ) 称为隐式差分格式。 ( 3 ) 当口= o 5 时，亓。皇尝，即取时段初和时段末水头值平均值近似地 作为时段的平均水头值，在这种情况下，有： p p 产 + 陋p 萨) 。仁) 他9 ， 式( 3 1 2 9 ) 称为中心差分格式。 3 1 4 数值实施 3 1 4 1 边界条件的处理 在计算土坡暂态渗流场时，边界条件的处理是一个重要的问题，国内外对降 雨条件下非饱和土土体稳定的研究，一般在处理计算入渗量时，通常有这样几种 处理方法： ( 1 ) 以降雨强度作为流量边界的水流计算入渗速率。 ( 2 ) 把降雨强度按一定不变比例减小，作为流量边界的水流计算入渗速率。 ( 3 ) 考虑降雨强度与土体渗透性的关系，计算入渗速率。当降雨强度高于饱 1 4 硕士学位论文 和渗透系数七。时，计算入渗速率按七，取值，否则按降雨强度取值。 ( 4 ) 当降雨强度小于土壤表面的入渗能力时，计算入渗速率取为降雨强度， 边界条件为第二类边界条件；当降雨强度大于土壤表面的入渗能力时，入渗的强 度就等于土壤的入渗能力，边界条件转换为第一类边界条件。 上述几种情况中，方法( 1 ) 没有考虑到降雨入渗随土壤入渗能力的变化而改 变的特点。方法( 2 ) 对边界条件的处理过于简单。第( 3 ) 种方法和第( 4 ) 种方法近似， 但仍有差异。为了较好的模拟降雨入渗不同阶段的变化情况，本文采取了第( 4 ) 种方法的思想，即考虑入渗率是随土壤入渗能力的变化而改变的特点，对坡面的 流量以边界条件的形式加以控制。并根据上节所述入渗理论调整实际的入渗流 量，实现降雨入渗流量的模拟。 将分布在单元面上的流量转化为结点的入渗流量。设某单元面有实际入渗流 量口。( f ) ，则入渗流量列阵可通过下式计算 r ) e 一心o ) ( r 出 ( 3 1 3 0 ) ， j 式中：s 为承受降雨入渗的单元面，忸z 表示t 时刻的单元入渗流量列阵，如 。 为函数向量 实际入渗流量留。( f ) ，是随着降雨强度和降雨持时的变化而变化，假设某时刻 降雨强度，o ) ，则降雨强度f ，( f ) 在法线上的分量为： 口，o ) = ，o _ ： ( 3 1 3 1 ) 要确定实际降雨入渗流量，还需要计算实际入渗能力。单元面上任意一点处 的压力水头为： j i 一 怕y ( 3 1 3 2 ) 则可以得到任意一点两个方向上的实际入渗能力的计算公式： 尺j o ) ，k 伽) 掣 日) | ( 3 1 3 3 ) 一 ” 然后根据下式转化为法线方向上的入渗能力尺( f ) 尺o ) 一r o h ( 3 1 3 4 ) 比较吼o ) 和尺母) ，以确定实际入渗流量吼o ) 。在实际计算中，有时坡面会出 现正压力水头。所以，应随时根据坡面土体的压力水头对流量边界进行调整。当 坡面结点压力水头为负压时，继续按照上述方法进行计算，但是当坡面压力水头 为正时，则将该坡面定为定水头边界。一部分雨水沿坡面流失，会在坡面形成一 薄层水膜，此时按给定水头边界处理由于水膜很薄，计算中取水头值等于地表 高程。 3 1 4 2 迭代计算方法 由于饱和非饱和渗流有限元计算格式的非线性，介质的渗透系数七、体积 变化系数研。需根据渗流场中的压力水头定，而压力水头又是待求量，需迭代求 降雨条件下的边坡渗流数值模拟及稳定性分析 解。为提高收敛速度，减少工作量，这里采用p i c a r d 迭代方法。对非稳定流问 题，具体迭代方法和步骤如下： ( 1 ) 设定一个迭代收敛条件：i j l 乏t 一时i 一 ，其中上标七代表第七时步，下 标，+ 1 和z 分别代表第“1 次和第，次迭代序号，下标m 缸表示求全域内的最大值， 即所有结点压力水头前后两次迭代结果之间的最大差值，p 为前后两次迭代计算 结果问的误差允许值。 ( 2 ) 由气时刻的压力水头场j i l i ，确定各单元的渗透系数七和体积变化系数m ；， 计算出有限元格式中的各个系数，然后进行下一步计算。 ( 3 ) 将晚1 与群进行比较，若满足收敛条件，则j i l 等一 0 结束迭代，进入下 一时步，否则，把j i l 等作为砧“，重复第( 2 ) 步和第( 3 ) 步直至满足收敛条件或迭代 次数超过规定的上限值( 对非稳定流问题，若在规定的迭代次数内不能满足收敛 条件，还需缩小时步长，重新累计时间和时步，重新开始迭代) 。 上述迭代步骤使总体流动方程线性化，然后用高斯消元法联立求解。收敛速 度取决于渗透系数非线性的程度和问题在空间上的离散化情况。较陡的渗透函数 需要较多的迭代次数和较大的设定收敛误差，离散化中较小的单元尺寸和较短的 时间步长有助于较快的收敛，且误差较小，但费时较多。 3 1 4 3 从属量的求解 一旦获得系统各结点的收敛水头，便可认为己求得非稳态渗流方程。某时间 步长内达到收敛后，便可由收敛的结点水头算出其他从属量，如孔隙水压力，继 而可以求出水头梯度、水流速度等。 结点孔