（地质工程专业论文）雨水入渗对土坡稳定的影响分析与模拟.pdf

桂椿工学院硕土，：芋位论文 摘要 雨水入渗诱发的滑坡在世界上分布最广，发生频率最高，危害相当大，引起 了各国政府和机构的广泛重视。大量统计资料表明，绝大多数的滑坡是发生在降 阿期间或降阿之后。雨水入渗使土体饱和度增加，吸力锐减并引起抗剪强度大幅 下降。为了分析降雨入渗对坡体稳定的影响，需要对雨水入渗过程进行正确的分 析和模拟。 本文以土质边坡为研究对象，首先阐述了饱和一非饱和渗流分析基本理论和 微分方程：并对影响边坡稳定的各类因素进行了论述，主要从水对土体抗剪强度 的影响、雨水的冲刷作用、雨水对土体的化学作用和降雨的强度、持时、雨型等 几个方面深入研究了雨水对边坡稳定的影响。针对目前边坡非饱和渗流一般采用 不考虑体变及气帽流动的雨水入渗计算方法，详细介绍了f l a c 2 0 中渗流计算的水 气两相流基本控制方程和数值计算公式，根据重庆某一土质边坡的特点，采用 f l a c 2 0 二维快速拉格朗日分析程序建立了降雨入渗边坡的两相流模型，进行了水 气两相流一应力的耦合计算及稳定性计算，讨论了边坡稳定状态分析和降雨持时、 降雨强度对边坡稳定性的影响。得到了一些有益的结论： m 水入渗浸润造成水分填充土体孔隙，空气随之排出。气相渗透系数是与 水相渗透系数成反比的，在降雨初期，边坡土体透气性较好，随着降雨入渗，空 气向坡面外流出，而随着表层土体的饱和度增加，气相的渗透系数降低，空气一 部分流出坡面以外，一部分则压入边坡以内，随着表层土体的饱和度的增加，气 体通过表层渗透出去的越来越少，水分浸润产生的气体绝大部分压入边坡内非饱 和区域。 降f ：阿一部分被表层的非饱和土体吸收，增加了含水量，提高了饱和度，另 一部分则流向坡脚，使得坡脚水位迅速抬高，极易引起坡脚发生破坏。 随着降冉强度和持时的提高，降雨入渗深度越深，表层土体的含水量提高， 基质吸力减小，边坡稳定系数越小。当降阿持时增加时，边坡的滑动破坏面有向 浅层移动的趋势。 关键词：| 二坡；稳定性；m 水入渗； 两千h 流 桂林工学院硕士学位论文 a b s t r a c t s l o p es t a b i l i t yi s s u eh a sb e e nam a j o rs t u d yo fr o c k m e c h a n i c sa n ds o i lm e c h a n i c s s o i ls l o p es t a b i l i t yp r o b l e m sr e l a t et oa r t i f i c i a la n dn a t u r a lf i e l d ss u c ha s m i n i n g e n 百n e e r i n g ，b r i d g ee n g n e e r i n g ，h y d r o l o g i c a le n g i n e e r i n g ，m u n i c i p a lw o r k s ，c a n a l a n dp o r ta n dw a s t er e s i d u ea n dg a r b a g em a n a g e m e n te n g i n e e r i n g ，h i l l s i d ea n d r i v e r s i d e ，e t c d e v e l o p m e n to ft h ew e s t e r nr e g i o nw i l lb ef a c e dw i t hal o to fs l o p e d i s a s t e r s ，s l o p es t a b i l i t yw i l lb e c o m em o r ei m p o r t a n t t h el a n d s l i d ew h i c hi si n d u c e db yr a i n w a t e ri n f i l t r a t i o nh a sb e e nf o u n di nt h e w o r l d ，a n dh a st h ef r e q u e n c yh i g ha n dt h eh a r mi sq u i t eb i g ，a n dh a sc a u s e dt h ev a r i o u s c o u n t r i e s g o v e r n m e n ta n dt h e o r g a n i z a t i o nw i d e l yt a k e s s t a t i s t i c ss h o wt h a tt h e o v e r w h e l m i n gm a j o r i t yl a n d s l i d ei so c c u r si nr a i n f a l ld u r i n go ra f t e rr a i n f a l l s a t u r a t e d s o i li n c r e a s e ss ot h a tr a i n w a t e ri n f i l t r a t i o n ，s h e a rs t r e n g t hw h i c hi sc a u s e db yt h ed r o po f s u c t i o ns u b s t a n t i a l l yr e d u c e d i no r d e rt oa n a l y z et h ei m p a c to fr a i n f a l lo nt h es l o p e s t a b i l i t y ，w en e e dt oa n a l y s i sa n ds i m u l a t i o nt h er a i n w a t e ri n f i l t r a t i o np r o c e s sc o r r e c t l y t h i sa r t i c l et a k e st h es o i ls l o p ea st h er e s e a r c ho b j e c t ，f i r s te x p o u n d e do nt h e s a t u r a t e d - u n s a t u r a t e df l o wa n a l y s i so fb a s i ct h e o r ya n dd i f f e r e n t i a le q u a t i o n s s o i l w a t e r c h a r a c t e r i s t i cc u r v ea n d p e r m e a b i l i t y m a t h e m a t i c a l e x p r e s s i o n a n d s a t u r a t e d u n s a t u r a t e df l o we q u a t i o n sa n db o u n d a r yc o n d i t i o n s ；s l o p es t a b i l i t yo ft h e v a r i o u si n f l u e n c i n gf a c t o r s ( i n t e r n a la n de x t e r n a l ) i sd i s c u s s e d ，s u c ha st h ei m p a c to ft h e w a t e ro nt h es h e a rs t r e n g t ho fs o i l ，r a i n f a l le r o s i o n ，t h ec h e m i s t r y o fr a i no f ft h es o i l ，a n d r a i n f a l li n t e n s i t y , d u r a t i o ne t cm a k e st h ed i f f e r e n c et os l o p es t a b i l i t y p r e s e n t l y , f o rt h e c a l c u l a t i o no nt h eu n s a t u r a t e ds e e p a g eo fs l o p e sd u r i n gt h er a i n f a l l ，t h et h ec a l c u l a t i o n o ft h ea i r - w a t e rt w o p h a s ef l o wc o n t r o le q u a t i o n sa n dn u m e r i c a lf o r m u l ao ff i a 乎oi s i n t r o d u c e d a c c o r d i n gt o as o i ls l o p ec h a r a c t e r i s t i c so fc h o n g q i n gd i s t r i c t ，f a s t l a g r a n g i a na n a l y s i so fc o n t i n u ai sa p p l i e dt os e t u pt o f i n i s ht h et w o - p h a s ef l o w c a l c u l a t i o n s r a i n f a l li n t e n s i t ya n dd u r a t i o nt i m em a k e st h ed i f f e r e n c et os l o p es t a b i l i t y s o m eb e n e f i c i a lc o n c l u s i o n si so b t a i n e d ：f i r s t ，r a i nw a t e ri n f i l t r a t i o nf i l l ss o i lp o r e s ，a n d t h ea i ri so u t g a sp e r m e a b i l i t yc o e f f i c i e n ti si n v e r s e l yp r o p o r t i o n a lt ow a t e rv a p o r p e r m e a b i l i t yc o e f f i c i e n t a tt h eb e g i n n i n go ft h er a i n f a l l ，g a sm o v e s t ot h es o i lo ft h e s l o p es u r f a c ea st h es e e p a g eo ft h er a i n ，a n dw i t ht h ei n c r e m e n to ft h es a t u r a t i o n ，s o m e “ 桂林工学院硕士学位论文 o ft h eg a sm o v e so u to ft h es l o p e ，w h i l es o m ei sp r e s s e di n t oi t w i t ht h es a t u r a t i o no f t h es o i ls u r f a c e ，g a si n f i l t r a t i o nf r o mt h es u r f a c ei st h el e s sa n dl e s s ，t h em o s tg a so f w a t e ri n f i l t r a t i o ni n t ot h es l o p ei si n t ot h en o n s a t u r a t e dr e g i o n s e c o n d ，o n ep a r to ft h e r a i n w a t e ri sa b s o r b e db yt h eu n s a t u r a t e ds u r f a c es o i l ，w h i c hi n c r e a s e st h em o i s t u r e c o n t e n ta n ds a t u r a t i o n a n dt h eo t h e rp a r tf l o w st ot h es l o p et o e ，w h i c hh e a d su pt h e s l o p et o ea n di n d u c e ss l o p et o ef a i l u r e t h i r d ，w i t ht h ei n c r e m e n to ft h er a i n f a l li n t e n s i t y a n dt h ep a s s a g eo ft i m e ，t h ed e p t h so fi n f i l t r a t i n gi n t ot h ee m b a n k m e n tb e c o m e sl a r g e r ， t h es u c t i o na n dt h es l o p e s t a b i l i t ys a f e t yf a c t o rb e c o m e sl e s s k e yw o r d s ：s o i ls l o p e ；s t a b i l i t y ； r a i n w a t e ri n f i l t r a t i o n ；t w o p h a s ef l o w 桂林工学院硕士学位论文 研究生学位论文独创性声明和版权使用授权说明 独创性声明 本人声明：所呈交的论文是我个人在陈先华教授指导下进行的研究工作及取得的 研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含他人已 经发表或撰写过的研究成果，也不包含为获得桂林工学院或其它教育机构的学位或证 书而使用过的材料。对论文的完成提供过帮助的有关人员已在论文中作了明确的说明 并致以了谢意。 学位论文作者( 签字) ：盘鱼盟生 签字日期：丝堕：墨 版权使用授权说明 本人完全了解桂林工学院关于收集、保存、使用学位论文的规定，即：按照学校 要求提交学位论文的印刷本和电子版本；学校有权保存学位论文的印刷本和电子版， 并提供目录检索与阅览服务；学校可以采用影印、缩印、数字化或其它复制手段保存 论文；在不以赢利为目的前提下，学校可以公布论文的部分或全部内容。( 保密论文在 解密后遵守此规定) 学位论文作者( 签字) ： 指导教师签字： 签字日期： 哳嗍佛 罩复基毖： 桂林工学院硕士学位论文 1 1 研究背景与意义 第1 章绪论 地球表面的土大多是非饱和土，典型的非饱和土包括失陷性黄土、膨胀土、 热带残积土和人工填土。从力学角度来说，非饱和土与饱和土不同的最主要特征 之一是土中存在负的孔隙水压力。降雨入渗引起的负孔隙水压力的变化是造成许 多土坡失稳产生滑坡的原因。 雨水入渗诱发的滑坡在世界上分布最广，发生频率最高，危害相当大，引起 了各国政府和机构的广泛重视。大量统计资料表明，绝大多数的滑坡是发生在降 雨期间或降雨之后。“中国地质灾害数据库”记录了1 9 4 9 年至1 9 9 5 年期问发生在 我国的滑坡灾害，其中6 8 5 的滑坡是由降雨引起的。其它诱发因素( 如库水位升 降、地震、洪水、人工边坡、采矿等) 导致滑坡的总数仅占同一期间发生在我国 滑坡灾害的3 1 5 。中国典型滑坡一书中列举了9 0 多个滑坡实例。1 ，其中有 9 5 以上的滑坡都与降雨有着密切关系。中国重大地质灾害实例分析一文中所 例举的2 7 顶中国重大地质灾害中有1 5 项是由于暴雨引发的“。 大量土质边坡滑动的事例表明，滑坡大都发生在雨季，尤其是暴雨、连续雨 的参与。例如1 9 8 2 年四川的万县云阳等地普降大雨，降雨量6 0 0 7 0 0 m m ，结果诱 发了全地区大小数万处滑坡。1 9 7 9 年甘肃洒勒山地区降雨比往年高3 4 1 ，其中 7 - 9 月三个月的降雨量达5 0 9 m m ，较历年同期高6 2 ，在这一年的9 1 0 月份在洒 勒山山顶出现了裂缝。1 9 8 2 年1 0 1 2 月降雪量是历年降雪量的7 倍，到来年春季 雪融化，融水下渗引起坡体内的一系列变化，最终使得滑坡发生“1 。2 0 0 4 年，全 国共发生中等规模以上地质灾害8 7 5 起，其中滑坡5 7 2 起。同年9 月5 日，重庆 市和四川省部分地区遭遇百年一遇暴雨袭击，形成群发滑坡崩塌灾害，造成8 5 人 死亡，经济损失巨大这些滑坡都说明大气降雨是大多数滑坡的主要促发因素， 降雨通过一系列复杂的影响，导致滑坡产生。降雨入渗，使土体湿润，容重增加， 抗剪强度降低，对于边坡的稳定产生极不利的影响。所以很有必要深入研究降雨 对土坡稳定性的影响。部分近代滑坡与降雨的关系统计见表1 1 。 桂林工学院硕士学位论文 表1 1 部分近代滑坡与降雨的关系统计表“1 滑坡位置及类型体积( 1 0 m 3 )发生时间滑前降雨特征 恩施屯堡杨家山土滑6 3 0 01 9 8 0 8 4滑动当天降雨1 8 0 思施屯堡桃树坪土滑 1 9 81 9 9 8 8 0 滑前3 天连续降雨1 8 l m 长阳县城王家湾土滑3 6 1 9 9 6 7 5 7 月4 日降雨1 1 7m 株归文化乡百日场土滑 2 0 01 9 9 6 7 7 7 月2 4 日连续降雨2 7 7m 柿归万塘村孔家橹土滑2 51 9 9 6 7 ，47 月2 4 日连续降雨2 7 7m 巴东沿渡河镇罗链溪土滑 51 9 8 1 8 78 月6 日降雨1 5 2 5m m 云阳安阳镇土滑2 5 1 9 9 3 6 4 6 月2 - - 3 日连续降雨1 7 7m 降雨入渗诱发的滑坡不仅与降雨强度有关，还与降雨历时、地表径流排泄条 件、初始地下水位、滑坡岩土体的渗透特性等因素关系密切”。1 。在有膨胀土、残 积士等特殊土的山区，雨季时常会发生各类大小不一的滑坡。降雨入渗对土坡的 影响主要是考虑由于降雨导致渗流场的变化而引起的作用在土体上动水荷载和静 水荷载的增大以及由于土体中基质吸力的变化而引起的土体抗剪强度参数的降低 捌。水的渗入使土体质量增大。由于土体的毛细压力，非饱和区土体不是完全干 燥的，非饱和土体的抗剪强度与土体的饱和度密切相关，而且随着土体饱和度的 变化，土体抗剪强度参数变化十分敏感，在降雨入渗过程中，随着雨强和时间变 化，土体饱和区以及非饱和区的含水量也在不断发生变化，即：坡体内土体的抗 剪强度也在不断发生变化。 我国南方以及中南，西南部分地区广泛分布着膨胀土，残积土等粘性土。那 里的山区或人工边坡在雨季常发生滑坡，而且以浅层滑坡最为多见。这类滑坡产 生的机理是由于近地表浅层土多为非饱和粘性土，天气变化对其力学性状有很大 影响。雨水入渗使土体饱和度增加，吸力锐减并引起抗剪强度大幅下降。为了分 析降雨入渗对坡体稳定的影响，需要对雨水入渗过程进行正确的分析和模拟“”。 1 2 国内外研究现状 研究降雨入渗与滑坡关系的重点是研究降雨对边坡土体含水量和强度的影 响，在此基础上研究边坡的稳定性，即从内因上研究降雨对边坡稳定性的影响。 在降雨入渗前，边坡土体尤其是上部土体处于非饱和状态，在降雨过程中， 雨水从两个方面诱发滑坡：对边坡土体，降雨入渗起加载作用，即雨水使土体 2 桂林工学院硕士学位论文 的含水量增大，重度变大，从而使滑移面的剪应力增大；降雨改变边坡土体的 力学性能，导致其内聚力下降，基质吸力减小，抗剪强度降低。边坡土体的自重 增加和强度降低这两个不利因素在降雨过程中同时影响边坡的稳定性，达到一定 程度就会引发滑坡。如果把降雨条件下边坡土体的重量和强度看作是一种动态的 分布场，那么研究的重点是边坡的重量场和强度场在降雨过程中的变化情况，只 要得到了边坡实际的重力场和强度场分布，就可以利用较为成热的方法来分析边 坡的稳定性，进而得出降雨入渗对边坡稳定性的影响规律。近年来，人们在这一 方面进行了大量的研究工作。 1 2 1 饱和一非饱和渗流的研究现状 在早期的渗流分析中，人们一般考虑的是饱和区内的渗流场情况，岩土工程 中的渗流问题往往是通过绘制流网来求解的，这种方法只能求解边界条件、土体 状况简单的饱和稳定渗流问题“”。随着计算机技术的发展，人们对渗流问题的数 、值模拟进行了探讨，变网格法和固定网格法开始用于对渗流场的解析中来。但是 这两种解法由于计算过程繁琐，计算量大，而且没有考虑到非饱和区的水的流动， 不能全面真实地反映地下水的渗流状态。 一 国外从7 0 年代开始考虑非饱和区域的流动，即把饱和区和非饱和区耦合在一 起进行整体分析。压力水头在饱和区为正，而在非饱和区为负值，零压力面就是 自由水面，即饱和区与非饱和区的分界面。这样计算域内不再有自由水面边界， 使得计算简化，程序处理也比较简便。 1 9 7 3 年，n e u m a n 呻1 对土坝进行了饱和一非饱和渗流数值模拟，首先将有限单 元法用到饱和一非饱和渗流中，他提出了用不变网格法分析有自由面渗流的 g a l e r k i n 方法； l a m 和f r e d l u n d “”应用饱和一非饱和有限元渗流分析程序对一些堤坝渗流的 经典问题进行了求解； r a h a r d j o 等“”把非饱和土力学理论应用在新加坡残积土的研究，在对现场吸力 检测和饱和一非饱和有限元数值模拟时发现降雨入渗对现场吸力的影响取决于地 表的边界条件。 国内近十几年来在这方面也开展了很多工作： 张家发“9 1 建立了考虑降雨入渗补给条件的三维饱和一非饱和非稳定渗流数学 模型，并利用编写的程序对三峡船闸高边坡的饱和一非饱和渗流场进行了模拟，但 3 桂林工学院硕士学位论文 对复杂边界条件的处理还需进一步研究； 吴梦喜和高莲士“采用了改进的有限元法有效地消除了非饱和渗流中的数 值弥散现象，并对逸出面边界提出了新的处理方法； 彭华、陈尚法和陈胜宏。”对饱和一非饱和渗流的有限元法进行了改进，提出 了修正容水度的技术和加速迭代收敛技术的新方法，在此理论的基础上对水布娅 滑坡的渗流场进行了数值模拟，提出了排水优化方案： 朱文彬、刘宝探。”利用有限元法和有限差分法模拟了公路边坡在降雨条件下 的饱和一非饱和渗流过程； 陈虹、陈彤汹1 对飞来峡水利枢纽工程纵向导流围堰典型剖面处的渗流问题进 行了有限元数值模拟，分析了降雨条件和上游水位骤降对渗流场的影响，并提出 了一个简易的非饱和参数关系曲线的拟合方法： 朱岳明、龚道勇等o “矧对三维饱和一非饱和渗流的降雨入渗边界条件以及非 饱和逸出面边界条件进行了分析，提出并定义了极限含水率和蒸逸系数两个新概 念： 彭刚、罗先启、田斌嘶1 对黄腊石滑坡群石榴树包滑坡区进行了非饱和非稳定 渗流分析，研究了地下水位因降雨而变化的过程，给出了地下水位的变化规律。 1 2 2 雨水入渗条件下土坡稳定分析的研究现状 国内外学者在对降雨条件下的土坡的稳定性分析上做了很多研究： l u m b ”7 率先研究了香港地区降雨和滑坡的关系。他提出一种简化的一维垂直 入渗模型并用于计算湿润峰推进的速度及入渗停止后含水量再分布过程。然后根 据抗剪强度与饱和度的经验关系，并假定破坏面平行坡面，研究了地质条件和降 雨特性对斜坡稳定性的影响。考虑的具体影响因素包括土的入渗能力、土的初始 饱和度、降雨强度与降雨持时曲线、暴雨前降雨量、间歇性降雨以及坡面防渗情 况等。l u m b 研究的结论为：香港斜坡的稳定性由土的入渗能力来控制，某次特定 暴雨对斜坡稳定性的影响取决于该次暴雨过程属于多少年一遇的水乎，还取决于 该次暴雨前的降雨量。但l u m b 在研究中未考虑水平渗流分量的影响，且假定导水 率k 及扩散度d 均为常量，这些假定与实际情况差别较大。 f , e d l u n d 等“1 反复强调负孔隙水压力对土坡稳定性的影响。并重新研究了边 坡稳定分析中的安全因数计算公式，以便把j 下的和负的孔隙水压力都包括进去。 他通过香港陡坡的算例进行了暴雨情形下的瞬态渗流分析和边坡稳定分析。分析 4 桂林工学院硕士学位论文 中采用模拟降雨前稳态渗流条件下的孔隙水压力分布作为暴雨期间瞬态分析的初 始条件。渗流分析中假定孔隙空气压力不变。用g a l e r k i n 有限元法模拟稳态和瞬 态渗流场，通过渗流分析可以得出暴雨开始后不同时刻的孔隙水压力分布，从而 分析土坡在不同时刻的稳定性。分析结果表明：由于临界滑动面比较浅，负孔隙 水压力在抗剪强度发挥中起主要作用，所以暴雨期间安全因数下降很显著，并有 可能导致斜坡的破坏。 s a m m o r i 和t s n b o y a n m a 1 用g a l e r k i n 有限元法模拟暂念渗流过程，并对边坡 稳定性进行了参数研究。渗流分析采用的控制方程为r i c h a r d 方程，稳定性分析采 用简化b i s h o p 法。选择的变化参数包括斜坡长度、土层深度、横截面形状和土的 性质，但降雨强度保持不变。分析结果表明：对斜坡稳定性不利的参数取值情况 是较低的导水率，较长的斜坡，较浅的土层深度和凹形的斜坡表面。 a l o n s o 等。”进行了土坡二维非饱和渗流和极限平衡法的联合分析，渗流分析 中采用了考虑空气压力变化的耦合型控制方程。考虑的影响因素包括土的类型、 降雨持时、降雨强度、水分保持曲线的形状和土的渗透性。主要研究结论为：安 全因数变化所需的时间与土的渗透系数成反比，且饱和渗透系数k 。 o ，在非饱和区巩，) o 5 0 m m 地区f 王兰生等，1 9 9 1 ) 。 由表3 - 1 中数据可见，临界值是随地区地质、自然地理等背景情况而有所不同 的，应根据地区的实际情况确定。据柳源嘞】暴雨滑坡类型不同，临界降雨强度就 不同：i 晦界暴雨强度与地层岩性和结构有关；临界暴雨强度与植被有关；临界暴 雨强度与一个地区常年降水量有关。 特大暴雨的袭击，引发某一地区大量边坡失稳，但也是一次对区域边坡中潜 在不稳定的变形体的大清理，宏观上改善了地区边坡的稳定状况。因而，可以认 为只有边坡又经历了一段演变过程之后，再度遭受类似等级暴雨的袭击，才有可 能造成区域性大量边坡失稳现象的重现。据此，可以根据气象部门暴雨概率分析 数据，参照地区暴雨历史调查资料，对地区暴雨群体滑坡事件作出预测和趋势分 析。当然，也不排除个别遭受暴雨袭击尚未失稳的边坡，在今后的一般性暴雨中 失稳；或者局部小范围集中降雨引发的滑坡，但毕竟是一种局部现象，不会影响 影响对暴雨滑坡事件发生概率或周期的分析预测。 表3 - 1国内外暴雨滑坡起动降雨量建议判据 一次降雨过程雨量降雨强度 国家或地区 总量m m 占多年平均的 日降雨量m m时降雨量m m 巴西2 5 0 3 0 0 1 2 美国 2 5 0 日本 1 5 0 - - 3 0 0 加拿大 2 5 0 中国香港 1 0 0 7 0 四川盆地 5 0 桂林工学院硕士学位论文 第4 章f l a c 2 d 软件及两相渗流计算原理 4 1f l a c 2 d 软件简介 n a 乎o ( f a s tl a g , r a n # a na n a l y s i so fc o n t i n u a ，连续介质快速拉格朗日分析) 是由美国明尼苏达大学和美国i t a s c ac o n s u l t i n gg r o u pi n c 。开发的有限差分计算程 序，主要适用地质和岩土工程的力学分析。f l a c 求解使用了混和离散化技术、动 态松弛方法和有限差分方法3 种技术，从而将连续介质的动态演化过程转化为离 散节点的运动方程和离散单元的本构方程求解。更为精确和有效地模拟计算材料 的塑性破坏和塑性流动。它全部使用动力运动方程，较好地模拟系统的力学不平 衡到平衡的全过程，即首先由节点的应力和外力( 或速度) 变化和时间步长利用虚功 原理求节点不平衡力和速度，再根据单元的本构方程，由节点速度求单元的应变 增量应力( 或位移) 增量和总应力，进而进入新的循环。虽然其基本原理类同于离散 元法，但却可与有限元一样适用于多种材料模式与边界条件非规则区域的连续问 题求解，而且计算中利用的“混合离散化”技术，比常规有限元“降低完整性” 的方法在力学上更合理；同时，f l a c 求解中使用显式差分方法，不形成刚度矩阵， 可节约计算机存储空间，减少运算时问，提高解题速度，便于在微机上实现非线 性大变形问题的求解。主要可以模拟边坡的稳定状况，给出边坡滑移迹线，动态 的描述边坡的运动规律，给出预警信息：模拟洞室的开挖过程，对支护参数进行 优化设计等等。 f l a c 程序建立在拉格朗日算法基础上，特别适合模拟大变形和扭曲。f l a c 采用显式算法来获得模型全部运动方程( 包括内变量) 的时间步长解，从而可以追踪 材料的渐进破坏，这对研究岩土工程问题非常重要。f l a c 程序具有强大的后处理 功能，用户可以直接在屏幕上绘制或以文件形式创建和输出打印多种形式的图形。 使用者还可根据需要，将若干个变量合并在同- - n 图形中进行研究分析。 f l a c t 6 。1 可以模拟流体通过可渗透固体的流动，例如模拟地下水在土体中流 动。与f l a c 中一般的力学计算无关，流动建模既可以由其本身完成，也可以同 力学建模并行完成，以便获得流一固耦合作用的效果。固结就是一种流一固耦合 作用，在这一过程中，孔压的逐渐消散导致土体中发生沉降，这种性状的类型包 含了两种力学效应。首先，孔隙水压的变化导致有效应力的变化，有效应力又影 桂林y - 学院硕士学位论文 响固体的力学响应，如有效应力的减小可能产生塑性屈服。其次，区域内的流体 也会通过改变孔隙压力来响应力学体积的变化。 f l a c 程序设有多种本构模型： 空单元模型： 各向同性弹性材料模型； 横观各向同性弹性材料模型； 德鲁克一布拉格( d r u c k e r - p r a g e r ) 弹塑性模型： 莫尔一库伦( m o h r - c o u l o m b1 弹塑性模型： 遍布节理材料模型； 应变软化硬化塑性材料模型； 双线性应变硬化软化一致节理模型； 双屈服塑性材料模型； 改进的剑桥( c a m c l a y ) 模型。 f l a c 2 d 的求解使用了如下3 种数值计算方法： ( 1 ) 离散模型方法：连续介质被离散为若干互相连接的四节点单元，作用力均 被集中在节点上。 ( 2 ) 有限差分方法：变量关于空间和时间的一阶导数均采用有限差分来近似。 ( 3 ) 动态松弛方法：用质点运动方程求解，通过阻尼使系统衰减至平衡状态。 4 2f u c 2 d 基本原理 f l a c 基于显式差分法来求解运动方程和动力方程。首先将计算区域离散化， 分成若干二维单元，单元之间由节点联结，节点受荷载作用后，其运动方程可以 写成时间步长为t 的有限差分形式。在某一微小的时段内，作用于该节点的荷载 只对周围若干节点有影响。根据单元节点的速度变化和时段a t ，可求出单元之间 的相对位移，进而求出单元应变；利用单元材料的本构关系即可求出单元应力a 随着时段的增长，这一过程将扩展到整个计算区域。在此基础上，求出单元之间 的失衡力( 不平衡力) 。将失衡力重新作用到节点上，再进行下一步的迭代过程，直 到失衡力足够小或节点位移趋于平衡为止。f l a c 程序采用最大不平衡力来描述计 算的收敛过程。如果单元的最大不平衡力随着时步增加而逐渐趋于极小值，则计 算是稳定的，否则，计算就是不稳定的。 桂林工学院硕士学位论文 圈4 1 物理同格 ( | 四边形睢元 ：一y 嘲4 2 散学网格 固圈艮 c h ) 兰角彤舡元c ) 单元节点向量 嘲4 一： 计努荦元 拉格朗日法用差分方法分析，因此首先要生成网格，将物理网格( 图4 1 ) 映射 在数学网格上( 图4 2 ) ，这样数学网格上的某个编号为i ，j 的结点就与物理网格上 相应的结点的坐标x ， y 相对应，也可以想象为数学网格是一张橡皮做的网，经 拉扯以后可以变成物理网格的形状。 f _ a c 程序的基本原理与离散元相似，但它像有限元法那样适用于那些包含多 种材料模型以及复杂边界条件的连续介质问题的求解：拉格朗日元法用差分方法 求解，因此首先要将求解的区域划分成四边形的网格，计算单元如图4 - 3 所示，在 边界等不规则的地方也可以用三角形网格来拟合。 在内部，f l a c 将每一个单元再分为两组覆盖( 共4 个) 的常应变三角形单元， 如图4 3 b 所示。由高斯散度定理的一般形式推导得出三角形差分方程，其形式如 下： l i 啦一 羔鲥 式中，f 在封闭曲面边界周围的积分 吩曲面的单位法向量 ，标量、向量或张量 ( 4 1 ) 桂林工学院硕士学位论文 毛坐标向量 d s 增量弧长 对表面积a 积分。 定义在区域面积a 上梯度的均值f 如下： 2 知善么 山 将式( 4 2 ) 代入式( 4 1 ) 得： ( 卦知善凼 c 对于一个三角形子单元，式( 4 3 ) 的有限差分形式可变为： 5 三c 眦 。 4 ， 血是三角形一边的长度，对三边进行求和，( 厂) 取三边的平均值。 用每一边速度矢量的均值代替了，利用式( 4 4 ) ，将一个区域的应变速度毛用 节点速度表示，公式如下： 芸一击科”州2 饥缸 ( 4 - 5 ) 岛= 牾剖 6 ， 其中，1 和2 是三角形边界的两个连续节点( 见图4 - 3 c ) ，正，m 是节点1 的速 度。值得注意的是，如果节点之间的速率呈线性变化，式( 4 - 5 ) 与精确积分导出的 表达式是相同的。 式( 4 - 5 ) 和式( 4 6 ) 可以用来推导就变张量的所有分量。为了打印和作图的需要， 最大剪应变就用叫平面内6 9 第- - 个应变不变量分量来表示。 最大剪应变t 的表达式如下： r 一言【一勺) 2 + 镌】1 ，2 ( 4 7 ) f l a c 提供了联系单元应变和应变速率的单元变量，在当前时步，变量s s r 利 用基于式( 4 - 5 ) 的节点速率计算出某一时刻的最大剪应变率，而变量s s i 则利用基于 桂林工学院硕士学位论文 式( 4 7 ) 的节点位移计算出某一时刻的最大剪应变。用同样的方法计算出体积应变 率v s r 和体积应变v s i ，f i s h 函数f s r 和f s i 计算了全应变速率张量和应变增量张量 的分量。 4 3 两相流 f l a c 同时提供了可选的二相流模型，可在多孔介质中存在两种不能混合的流 体进行流体流动或完全耦合模拟。进行与应力耦合计算时，遵循以下过程，即有 效应力改变产生体变，体变产生流体压力的变化，从而反过来影响有效应力。同 时可以通过f i s h 语言实现对孔隙比的更新，从而带入新的孔隙比而对土体的渗透 性与基质吸力的瞳线作自动的调整。以下对f l a c 2 d 中两相流的基本理论作详细 介绍。 4 3 1 基本控制方程 f l a c 中的二相流控制方程包括以下公式。 ( 1 ) 达西定律 水相和气相流体输运公式由达西定律描述如下： # 一k t x 7 丢x ( 4 8 ) j ( p - p w g k x ) 彰一喵告群击( o 一& g t 以) “- 9 ) j s瓠| 式中，是饱和渗透系数，表示为张量形式；是流体相对渗透系数，是饱 和度的函数；是动粘滞度：p 是孔隙压力；p 是流体密度；g 是重力。在f l a c 中，渗透系数定义为固有渗透系数与动粘滞度的比率。 ( 2 ) 相对渗透系数公式 相对渗透率在v a i lg e n u c h t e n 经验公式里，是与饱和度有关的。见图4 4 ： ? 一【1 一( 1 一抽) 4 】2 ( 4 - 1 0 ) 衫一o - & ) 1 - ( 1 一蔓“) 4 】2 口 ( 4 - 1 1 ) 在这些公式中，4 ，b 和c 是常量，是有效饱和度，定义如下： 足一等 ( 4 1 2 ) 桂林工学院硕士学位论文 式中，是水相残余饱和度，即在较高毛细压力下依然存在的残余饱和度， 被认为是只要有水存在的情况下就存在。 图4 - 4 相对渗透系数。衫与有效饱和度s e 的关系曲线 ( 3 ) 毛细压力定律 毛细压力定律表示为两相流孔隙压力随饱和度的变化情况，以下式表示： 只一只一只( s 。) ( 4 - 1 3 ) 在f l a c 中，这个经验公式用v a ng e n u c h t e n 公式表达，见图4 5 。 只( & ) 一e o i s i ”。- 1 1 “ ( 4 1 4 ) 参数最对于细颗粒材料来说更大，其与材料性质的相关性可以用l x v e r e t t 比 例公式估计。这个公式由量纲分析导出，有以下形式： 昂2 赤 ( 4 1 5 ) 式中，仃是表面张力，是基质属性；茁是固有渗透率，n 是孔隙率。参数a 有 时在文献中代替岛。这两个比例系数的关系是： e o 。旦星 口 式1 j 凡是水的密度，g 是重力。 ( 4 1 6 ) 桂林工学院硕士学位论文 图4 - 5 毛细压力和有效饱和度s 。的关系 ( 4 ) 饱和度 当孔隙由这两种流体完全填满时，有： t s ，+ sg 一1 ( 5 ) 流体平衡公式 对于微可压缩流体，平衡关系是： 亟；一掣+ q 7 o ta t 鲁一一誓+ 衫出a x ： 其中，考是流体含量单位体积的变化量，吼是体积流体源强度。 ( 6 ) 流体本构方程 流体本构方程如下： ( 4 1 7 ) ( 4 1 8 ) ( 4 1 9 ) ( 4 2 0 ) ( 4 2 1 ) 式中，0rk 是流体体积弹性模量，是体积应变。将等式( 4 。2 0 ) 代x , ( 4 - 1 8 ) ， 等式( 4 2 1 ) 代入( 4 1 9 ) ，整理得： 拈一西 拈一甜 瓯 墨 一 一 堡以晦百 叫 叫 亟甜生甜k一雕巧一胛 暑 i 监以嵋百 & 墨 桂林工学院硕士学位论文 仃陪一等+ 堡a tj 1t 一匿+ 毛割 c 4 彩， 在只考虑渗流的计算中，譬项可以省略，等式( 4 2 2 ) ，( 4 2 3 ) ，( 4 1 3 ) 和( 4 1 7 ) 组成了一个以己，0 ，& ，长为未知量的非线性方程组。 二相流分析采用了流体与力学耦合计算的附加方程，形式如下： ( 7 ) 动量平衡 平衡方程如下： 孥+ p g i = p - 矾万t p g(424dx,a t ) - - 羔 ) 式中，p 是容重，应是流速。对于二相流计算，有： p p t t n 咚。p 。七sg pe 、 4 2 5 ) 式中，“，n 是流体密度，以是基质干密度。 ( 8 ) 力学本构公式 多孔固体的本构反应增量有如下形式： a h = h 世口，8 q ，心 “ 。vv 式中，a t r 叫是有效应力增量，日是本构公式的泛函形式， 对于二相流，有效应力的增量定义如下： a t r j 。= 厶o q ，世6 式中， ( 4 2 6 ) r 是过程参数。 ( 4 2 7 ) p s 。蛾+ 蝇 ( 4 - 2 8 ) 由此公式可以看出，假设总应变保持不变，变形将随着孔隙压力改变而改变。 因为本构模型包含塑性，b i s h o p 有效应力可以用来检验失效。 ( 9 ) 协调方程 流速梯度和应变率的关系如下： 岛一糟+ 割 z 9 ， 3 1 桂林工学院硕士学位论文 4 3 2 流体平衡定律的节点方程 在f l a c 中，流体平衡方程( 4 2 2 ) 和( 4 2 3 ) 是通过离散的有限差分法来求解的。 在两楣流计算中，每个节点处都定义有两种孔隙压力( 只，只) 以及饱和度( & ) 三 种参量，而且假定这三个参量满足一定线性三角关系。数值方案只限于流体连续 方程下的节点计算。节点的水头分担根据逆流加权技术计算。微分方程通过显式 法及时求解。 对于一个有限的体积n 平衡定律如下所示： 乏哦+ 哦一万1 + 隗缸心，】 ( 钙o ) 足s s b p 。+ 峨= 一万1 + 眩出+ s g 吲 ( 4 3 1 ) 式中，v 是节点体积，o 是节点流速。 式( 4 1 7 ) ：f 1 j ( 4 - 1 3 ) 的增量形式如下： 0 一a p w ；只岱。 ( 4 3 2 ) 。 a s 。+ b s g 一0 ( 4 - 3 3 ) 式中，芝表示只( & ) 对应s w 的导数。 上述的这两个关系式可以用于消去式( 4 3 0 ) 和( 4 3 0 e e n 弓和& 。通过一些 变换后，可以得到如下的式子： 。 志越+ 战= 一专+ 蛾出+ 毛a v j ( 4 - 3 4 ) 毒蛾+ ( 等。卜。_ - 专+ 嗡出+ & ， “_ 3 5 ) 最后，可以得到哦和b s 。的表达式： 哦一一嵩 瓯( - 一等卜卜篇一等卜】( 4 - 3 6 ， 峨一一嵩p 乏一q 羔卜篇5 ，& 毒一瓦1 c 4 拼， 桂林y - 学院硕士学位论文 d 。笠+ 旦一立监 ( 4 3 8 ) k ，k | k 。k g 不排水系数口是为了数值计算的方便丽引入的，在力学耦合中取恒值等于1 ， 在稳态流计算中为0 。 对公式( 4 1 4 ) 取s ，的导数，采用式( 4 1 2 ) 计算： 足乏- 只半( 1 一) “41 - “4 r ( 4 - 3 9 ) 从上面的表达式可以看出，虽然毛细压力曲线的导数在s w = 1 处有一垂直切 线，式( 4 3 6 ) 和( 缸3 8 ) 表示，区域墨e 是有界的，而且在饱和的情况下其值趋近于0 。 在f l a c 的显式数值计算中，通过时间t 处的己知值与由式( 4 - 3 6 ) 和( 4 3 7 ) 求的 增量相加