（岩土工程专业论文）运用离散单元法研究土体管涌机理.pdf

摘要 管涌是渗透变形的一种重要的破坏形式，给江河两岸的人民带来了巨大的经 济损失和人员伤亡。目前对于管涌的研究从不同的方向对管涌的原理进行了揭 示，采用的具体方法也不一样，但是由于管涌是由于颗粒的运移引起的，因此如 果能够从颗粒受力机理上进行探索，进而揭示管涌的机理可以说是比较有意义的 探索。 基于这样的想法，联系现有的数值计算方法，考虑离散单元法作为研究散粒 体比较有效的工具，因此设想利用离散元的方法来进行管涌机理的研究。 本文的工作主要包括： ( 1 ) 研究了离散元方法目前在工程中的运用，以及目前管涌研究所采取的 方法，得出离散元方法研究管涌的优势所在。 ( 2 ) 针对目前离散元方法中考虑流体作用较少的现状，仔细分析了离散元 计算的基本原理，分析将流体的计算方法，分析流体计算时应该采用的计算方法， 并在分析颗粒与流体间作用的基础上，将颗粒对流体间的作用作为流体方程中的 一项体积力，在此前提下计算流体动力学方程，同时颗粒受到流体施加的力，将 其作为颗粒受到的外力施加在颗粒上，通过不断的计算迭代，得出最终的结果。 ( 3 ) 选用颗粒流软件p f c 对管涌的发展过程进行数值模拟，采用粗细两种 颗粒进行模拟，尝试得出几何条件对管涌发生的影响，验证了s s t o m l i n s o n 和 yp v a l i d 通过反滤层试验得出的管涌发生的几何条件。 ( 4 ) 考虑管涌发生的几何条件，及水力条件，进行管涌计算模拟，得到了 和实际比较接近的结论，说明了离散元方法运用于管涌研究中的可行性，并借助 于p f c 软件的优势，得出颗粒起动力，与颗粒受到的阻力进行比较，得出管涌的 大致判别方式。 ( 5 ) 针对影响管涌的条件进行了定性的分析，讨论了计算区域的几何尺寸， 也即横截面积以及渗径对管涌的影响作用，以及渗透水头的施加速度对管涌发展 的影响。 关键词：管涌离散元p f c 颗粒孔隙率 ；基金项目：高等学校阵二i ：学科点专项科研基金资助( 2 0 0 4 0 2 9 4 0 0 3 ) a b s t r a c t b e c a u s ep i p i n g ，a so n eo ft h em o s ti n l p o r t a n td e s t i 。u c t i o nf o n no fs e e p a g e m e t 锄。叩h o s e ，h a sb r o u 曲ta b o u th u g el o s si ne c o n o m ya n ds a f e t y ，al o to fr e s e a r c h w o r ka b o u tt h ep i p i n gp r i n c i p l eh a sb e e nd o n ea td i a 、e r e n tp o i n t s ，i nw h i c ht h e d i s s i m i l a rm e t h o d sa r ea d o p t e d b e c a u s et h ep i p i n gi sc a u s e db yt h em i g r a t i o no ft h e n n ep a n i c l e st h r o u 曲t h ev o i d sf o 肿e db yt h ec o a r s ep a n i c l e s ，i tw i l lb em o s t s i g n m c a n tt os t u d yt h ep i p i n gm e c h a n i s mb ym i c r o c o s m i cm e t h o d b a s e do nt h ea b o v e ，r e l a t e dt ot h ea n a l y s i so ft h ec u n n td i h 色r e n tn u m e r i c a l m e t h o d s ，t h ed i s t i n c te l e m e n tm e t h o di sr e g a r d e da sam o r ee f r e c t i v et o o lt o 铲a n u l a r m a t e r i a l s t h e r e f o r et h ed e mi sf i r s t l yu s e df o rt h er e s e a u r c ho fp i p i n gm e c h a n i s m t h em a j o rw o r ki nt h et h e s i si sd e s c r i b e da sf o l l o w s ： f i r s t l y t h ec u r r e n tu t i l i z a t i o no fd i s t i n c te l e m e n tm e t h o di n t 1 1 e p r a c t i c ei s s t u d i e d ，a sw e l la st h ep r e s e n tr e s e a r c hm e t h o df o rp i p i n g ，t h es u p e r i o r i t ) ro fd i s t i n c t e l e m e n tm e t h o du s e di np i p i n gr e s e a r c hi sc o n c l u d e d s e c o n d l y ，b e c a u s et h en u i de f r e c ti ss e l d o mc o n s i d e r e di nt h ec u r r e n tr e s e a r c h 、v o r ka b o u td e m ，t h eb a s i cc o m p u t a t i o n a lp r i n c i p l eo ft h ed i s t i n c te l e m e n tm e t h o di s c a r e 如l l ya n a l y z e d t h ec o m p u t a t i o n a lm e t l l o df o rn u i di sp r e s e n t e dw h i c hm a yb e m o r es u i t a b l ef o rt h eu s eo fd e m o nt h eb a s i so ft 1 1 ea i l a l y s i so fm ei n t e r a c t i o nf o r c e b e t w e e ng m u l a rp a r t i c l e sa n dn u i d s ， t h eh y d r o d y n a m i c se q u a t i o n si sc o n d u c t e dt o d e t e m l i n et 1 1 ef o r c eo fn u i d sa c t i n go nt l l ep a n i c l e s ，a t 也es 锄et i m et h ef o r c ei sa l s o a c t e do nt 1 1 eg r a n u l a rp a n i c l e sa se x t e m a lf o r c e s t h ef i n a lr e s u l t sc a nb eo b t a i n e db y t h ea b o v ei t e r a t i o n s t h ec o u r s eo fp i p i n gi ss i m u l a t e db ys e l e c t i n gp a n i c l en o wc o d e ，c o a r s ea n d f i n ep a r t i c l e sa r es e l e c t e dt os i m u l a t e ，t h ee f r e c t so fg e o m e t r i cc o n d i t i o nt op i p i n ga r e a t t e n l p t e dt og e to u t s s t o m l i n s o na n dy p v a l i d se x p e r i m e n tt h r o u 曲f i l t e rl a y e r i n v e s t i g a t i o n t b o u tt 1 1 eg e o m e t r i cc o n d i t i o no fp i p i n gi sv a l i d a t e d f i n a l l y ，w i t hc o n s i d e r a t i o no ft h eg e o m e t r yc o n d i t i o na n dt h ew a t e rp o w e r c o n d i t i o n ，t h es i m u l a t i o no fp i p i n gi sp u to n ，a n dt h ec o n c l u s i o n sa p p r o a c h e dt ot h e p a n i c ei so b t a i n e d t h ef e a s i b i l i t yo fu s i n gd i s t i n c te l e m e n tm e m o di nm ep i p i n g r e s e a r c hi sp r o v e d w i t ht 1 1 es u p e r i o r i t yo ft h ep f cs o r w a r e ，t 1 1 ed r a gf o r c ei s d e t e m l i n e d ，a n db yt h ec o m p a r i s i o nw i t ht h er e s i s t a n c ef o r c e ，t h ea p p r o x i m a t e f o n n u l af o rd e t e n n i n gt h ep o s s i b i l i t yo ft h ep i p i n gi sp r e s e n t e d k e yw o r d s ：p i p i n g ； d e m ；p f c ； g r a n u l a r ；p o r o s i t y f u n dp r o j e c t s ：s u p p o n e db ye d u c a t i o np r o g r 锄f o rd o c t o r a t e so fs p e c i a lr e s e a r c h 如n d so fc o l l e g e sa n du n i v e r s i t y ( 2 0 0 4 0 2 9 4 0 0 3 ) 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论 文中不包含其他人已经发表或撰写过的研究成果。与我一同工作的同事 对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。 如不实，本人负全部责任。 论文作者( 签名) ： 2 0 0 6 钤殷归 学位论文使用授权说明 河海大学、中国科学技术信息研究所( 含万方数据库) 、国家图书 馆、中国学术期刊( 光盘版) 电子杂志社有权保留本人所送交学位论文 的复印件或电子文档，可以采用影印、缩印或其他复制手段保存论文。 本人电子文档的内容和纸质论文的内容相一致。除在保密期内的保密论 文外，允许论文被查阅和借阅。论文全部或部分内容的公布( 包括刊登) 授权河海大学研究生院办理。 论文作者( 签名) ： 2 0 0 6 年岁月州汨 第一章绪论 第一章绪论 第一节研究背景与意义 近年来，长江淮河及洞庭湖一带连续遭受历史罕见的洪水灾害，干堤溃决给 两岸人民的生命财产安全带来巨大的损失：特别是长江中下游及沿岸平原地区是 中国人口最密集，经济最发达的地区。这些地区的洪涝灾害最严重，严重影响着 我国国民经济和社会的持续发展。因此防洪成了汛期当地甚至全国的重要工作。 堤防主要是第四系河湖相冲积层的基础上逐步加高而形成的，在筑堤过程中一般 是就地取材，地基大都未作处理，因而常形成隐患。 堤防溃堤主要有两种形式：一种是河流静水压力和动水压力的破坏作用， 另一种是地下水渗流破坏j 前者可以看成为外力作用，如崩岸中的窝崩即为一典 型方式。后者在渗流力作用下，由于汛期江河水位上涨，江水渗透穿过堤身、地 基，在堤内坡，或坡脚逸出，形成潮湿、湿润、松散或浸润水流，称散浸。其进 一步发展则会形成漏洞，从清水漏洞发展到浑水漏洞，再发展则形成管涌。管涌 继续发展则可能现成管漏。管漏是指土中洞穴管道中的水流一般为紊流作用，而 管涌在考虑渗透水流的作用时，一般认为属于层流的范畴。 渗流的几种破坏形式中破坏作用比较常见的是流土和管涌，据统计，1 9 9 8 年长江中游在地基问题引起的险情中9 0 是由流土和管涌所致，就其宏观上的统 计得出的原因有三：一是地基条件差，堤下多为沙基或人工填土地基，且多数堤 防未作地基处理；二是堤身主体结构质量差，与其历史成因有关，大多为原有堤防 加高而成；三是堤后地形地貌所致，如出现坑塘，覆盖层局部比较薄弱等。当然 流土与管涌并不是只是在堤防中发生，在基坑开挖过程中，随着越来越多的超高 层建筑的兴起，在基坑开挖的时候需要相应地加大开挖深度，这必然会导致坑内 外水头差的增大，极易导致流土以及管涌破坏。 管涌和流土是工程地质、水利工程领域是经常提到的渗透变形现象。按国标 g b t 5 0 2 9 7 9 8 和g b t 5 0 2 8 6 9 8 的定义：土体中渗透水流的动水压力在某处变大， 将微小土体带出的现象叫流土，可发生在无粘性土和粘性土中。当动水压力不足 以带走整个微小土块，而是将其中较小颗粒带走的现象称为管涌。它主要发生在 无粘性土和分散性粘土中。管涌又可以分为发展型及不发展型两种。不发展型管 河海人学硕：匕学位论文 涌危害较小，但是在很多情况下，管涌发生后随着细颗粒的不断析出和流失，孔 洞的直径不断变大，它的深度也逐渐向堤身或地基土内部不断延伸，形成集中渗 透通道，严重时引发堤防不均匀沉降及整体失稳。 在理论上，对于流土型渗流破坏的研究，自从1 9 4 3 年太沙基根据单位土体 受力平衡提出著名的临界水力梯度后，经其本人和其合作者及其它学者的努力， 已经形成比较完善的理论，得到发生流土的临界水力梯度，但是对于管涌的理论 研究至今尚不成熟，有待进一步加以研究，本文拟针对实际生活中危害较大的管 涌现象进行力学机理方面的探讨，以得出一定的成果。 目前用临界状态的相关公式来研究管涌大都面临一个问题：当在确定流土的 破坏时是对单位土体进行受力分析，考虑的是土体在水中的浮容重和作用于该体 积的渗透力相平衡的原理，从而得出流土的抗渗比降。而在进行管涌的破坏机理 研究时采用的大体都是针对于单个颗粒的力，包括渗流力以区别与流土的破坏。 但是管涌是细颗粒在孔隙中移动并被带出的现象，而不是土块整个的启动。所以 必须考虑颗粒在移动过程中所受到周围颗粒及水体对其的作用力，即土水之间的 相互影响是不可忽视的因素。另外，流土的发生是一个突发的事件，某处的水头 梯度一旦超过临界水头梯度，就会立即在粘性土中发生局部隆起，而管涌的发展 过程是细颗粒在粗颗粒的空隙中被不断被带出，具有很强的时间效应，是一渐变 的过程( 汪时敏1 9 9 3 ) ，随着水头梯度的升高，土体中的颗粒从细到粗依次起动 并被水流带走最后形成管涌通道，造成管道漏水。有时出于简化模型考虑，管涌 产生初期，将涌水点附近局部区域的透水层看作是各向均质，实际上在管涌尚未 发生及其发展过程中，各项指标如渗透系数、孔隙率等都是随时间和空间变化的 量。在实际工程中我们发现管涌的发生与地层中土的组成成分、结构、土的级配、 水力梯度、管涌发生的距离、深度、表面覆盖粘土层的内摩擦角、覆盖层厚度、 粘滞系数、土的饱和度、固结系数、浸泡时间等因素有关，是一个多元的复杂问 题。 如果从单个颗粒受力平衡的角度着手，就会发现在进行受力分析时，单个颗 粒在力的作用下，产生加速度，然后产生速度，发生位移，必然和周围颗粒发生 作用。当然可以看到，对于较少的颗粒集合，处理起来不太困难，但是对于数目 巨大的颗粒集合体，即散粒体，需要有一比较好的计算方法对在渗流力的作用下， 颗粒受力r 运动状态进行动态模拟。夺l 北情况下，以往在进行宏观分析时的一些 第一章绪论 影响因数就不需要作为主要的条件来进行考虑，以往在进行宏观分析时的一些宏 观影响因素就不需要作为主要的条件进行考虑，必须从细观的角度来考虑结构的 稳定性。本文从微观或者细观上来进行分析，从而离散元的概念得以引入。 第二节管涌的研究现状 管涌临界水力梯度的理论研究至今尚不成熟，其原因是发生管涌的渗流机理 还没有在理论上取得突破，并且由于颗粒大小、级配、形状各不相同，采用实验 手段在探索其规律方面也有很大的困难。因而对管涌的预测一直是个比较有意义 的课题。从半个世纪前，地质界和土工界开始综合分析管涌的破坏机理，以及模 拟管涌发展过程。在此期间通过大量实验和分析土体的颗粒级配曲线，发展了常 见的几种方法。针对管涌大多发生在无粘性土中的实际状况，目前大多数的研究 都集中在无粘性土中，如果管涌在粘性土中发生，其产生的机理与在无粘性土中 的完全不同，目前一般认为起因可归结为土体裂隙等损伤的存在，并进一步引起 粘性力的丧失从而引起的侵蚀破坏。当然如果从微观上来看，不管管涌发生在何 种类型的土中，其基本的原理都是土颗粒在孔隙孔道中的运动，不同之处在于运 动颗粒的直径和所受到的周围介质对其的作用力的变化有许多学者从事孔隙中 颗粒的运动的研究如：r a a t s ( 1 9 8 1 ) ，m c d o w e l l - b o y ee ta l ( 1 9 8 6 ) 等，当然通过模拟 孔隙中颗粒的运动来讨论反滤层的设计中的基土颗粒的启动( h o n j oa n d v e n e z i a n o1 9 8 9 ；i i l d r a r a t n aa n dv a f a i1 9 9 7 ) ，得到了许多启发性的结果，这些成 果表明，从孔隙中细颗粒的受力运动入手来研究渗透变形是一个行之有效的途 径。 传统研究判定管涌的方法有： ( 1 ) 依斯托美娜的不均匀系数法 ( 2 ) a b e r g ( 1 9 9 3 ) 的颗粒- 孑l 隙链索模型 ( 3 ) 巴特拉雷夫等的土体空隙直径与填料粒径对比法 ( 4 ) 刘杰【1 j 和沙金煊土体细粒含量判别法 ( 5 ) 水利水电科学研究院方法 并由此得出了相关的一些临界水力梯度的实用公式。如：依斯托美娜公式、 刘杰公式、沙金煊公式等。 ”- i 卜：，。 河海大学硕二l ：学位论文 由于土体中颗粒间的相对位置是随机的，孔隙的形状是不规则的，以及颗粒 之间的相互作用力带有不确定性，从而不同的颗粒具有不同的运动状态，因此有 学者将管涌是看作是一随机的过程，认为必须采用概率概念来分析管涌问题。如 曹敦履【2 】【3 】8 0 年代中期就提出了渗流管涌的随机模型，根据己知参数如：渗流坡 降及临界渗流坡降的分布函数，用以计算管涌破坏的概率并估计土的渗流稳定 性，从概率的范畴上取得了很好的效果。 目前在工程中普遍使用的数值计算方法为有限单元法，当然，目前研究管涌 绝大部分也都是运用有限元的方法对工程进行模拟计算，而且这些都是从宏观上 对管涌进行分析。 殷建华【4 】利用有限元方法，模拟了管涌区长度和渗透性对堤身渗流的影响， 发现管涌区越长和渗透系数越大，流向或流入管涌区的水量就越多，从而导致渗 流速度的增大，作用在土颗粒上的作用力就更大，所以推测会有更大的颗粒启动， 从而导致管涌区范围的增大以及渗透系数的提高。从中也可以看出了解管涌发展 过程中的渗流场的变化情况十分必要。 张家剔5 】针对长江中下游堤防工程现状概化出典型条件，采用有限元数值模 拟方法对地基渗流场进行了模拟计算。谢红、张家发通过数值模拟分析了心墙堆 石坝、坝基及两岸山体渗流分布规律和防渗帷幕的作用效果。 朱伟【6 】结合日本阿武限川的地基渗漏防治工程中遇到的问题，应用有限元饱 和一非饱和渗流解析，对地基渗透破坏发生机制及其影响因素作了分析和讨论。 其在假设表层为弱透水层下卧较厚强透水层，也即所谓的工程中研究管涌常用的 双地基层。在此基础上在考虑管涌的临界水头的影响因素时，考虑表土层渗透系 数和下卧层的渗透系数比值，在此比值下，在考虑降雨入渗及洪水波峰形状的影 响情况下就表土层的厚度对管涌发生的影响做了定性的分析。 其他一些学者则分别从不同的出发点对管涌发生的机理进行了探讨。 陈建生、李兴文【_ 7 j 等对堤防渗流管涌发生后产生集中渗漏通道的机理进行了 深入的分析探讨，认为：管涌发生之前，覆盖层下部的静水压力基本上就是堤防 上游江水位，因为粘土的渗透系数比透水层中砂土要小几个数量级，这时江水与 j 地下水的水头差基本上都作用在覆盖层上。一旦覆盖层被水压力顶穿，四面八方 的水都集中涌向管涌处。在管涌发生初期采用井流理论模拟江水向管涌口补给时 地下水的流场情况，确定管涌盯始发生时饷范围以及逐步的发展过程，根据管涌 4 第一章绪论 的临界水力梯度，通过模型可以求出管涌初期的临界面，管涌初期涌砂区的范围 较大，由于管涌离堤坝最近的地层中的水力梯度最大，被带走的砂也最多，涌水 量增加而水力梯度减小，造成接近管涌初期临界面附近的水力梯度达不到地层 颗粒移动的临界水力梯度，造成临界面向里缩小，最终形成了集中渗漏带。 何永红【8 】则认为凡是渗流就有渗流力，当动水压力达到一定的值后，土中某 些颗粒就会克服阻力，在渗透水头作用下被带走，发生潜蚀，使结构变松，强度 降低，进而可以产生渗透变形，但是其认为，要使颗粒在渗透力的作用下松动以 至悬浮，必须克服颗粒间的内摩擦力与以及颗粒的水下浮容重。土体中所有单个 颗粒，在某种程度上所产生的渗流运动阻力会引起局部压力差，局部压力差之和 形成渗流中的总的压力差。但是，颗粒云以事实上在水头达到一定的值以前在研 究区域内就存在颗粒的运移，而且运移的范围也不仅仅是在管涌口附近的区域， 再整个研究范围内都可能发生颗粒的相对运动。在土体内部，存在小颗粒在渗透 力的作用下试图克服阻力起动，从而导致颗粒互相碰撞，这一颗粒间的接触作用 力在研究管涌发生的过程中是不能不予以考虑的。 国家自然科学基金委员会与长江科学院联合资助的重大项目“洪水特性与减 灾方法研究”的子课题“堤防破坏机理和安全评价方法研究”中同济大学的白冰 与周健【9 】针对流砂与管涌等破坏机理，建立了变化形成和发展的理论模型，应用 “颗粒流理论”考虑地下水管状渗流与渗流介质相互作用及水流的贴壁挤土效 应，建立管涌形成模型，并分析破坏的发生和发展过程。在模拟管涌形成与发展 过程中，引进“颗粒流理论”计算方法将渗流介质看作由颗粒组成的流体，分析 土颗粒渗流场的变化特征。将已有的管状渗流模型与地一卜水渗流模型进行耦合， 考虑管状渗流与渗透介质的共同作用。 李守刹1 0 】，张晓海，刘志祥对于基坑开挖工程中管涌的发生过程进行了研究， 提出使用有限元迭代搜索方法模拟基坑工程管涌发生过程。 兰州大学的刘忠玉在其博士论文无粘性土中管涌的机理研究【1 1 j 及无粘 性土管涌的毛管模型及其应用【1 2 j 文中对管涌发生阶段中土水相互作用进行了相 当的研究并取得了一定的成果。其在分析单个颗粒引起的压力降和考虑可动颗粒 的随机性，建立了管涌的随机模型，可以模拟管涌过程中士体颗粒级配变化。通 过计算，发现拖曳力也即渗流力的标准差对可动颗粒的启动、管涌临界水力梯度 i ? -、 有很大影响。 ： 河海大学硕士学位论文 在水头值达到一定值时才能发生渗透，且在发生渗流后，由于颗粒及水体之 间的各种作用力的影响也不是呈线性关系。在粗颗粒土中，只有在小的水力梯度 下，渗流速度与水力梯度才成线性关系。而在较大的水力梯度作用下，水在土体 中的流动进入紊流状态，达西定律此时不再适用，从层流到紊流，这之间必然存 在一动态转换过程，在此过程中，砂土中各处的颗粒组成和渗透系数必然是随时 间的变化量。因此在研究管涌和渗流时最关切的参数是土的渗透性即渗透张量。 对于砂土在层流状态时一般认为符合达西定律，即渗流速度和水头梯度之间呈线 性关系： y = 七f ( 1 2 - 1 ) 式中，七为渗透系数，f 为水力梯度。 对于密实粘土，由于吸着水具有较大的粘滞阻力，只有达到一定的起动水头 才能发生流动。用有限元来模拟管涌时其原理是将计算区域划分为有限单元来进 行求解及模拟，必然面临当单元内颗粒组成发生变化时如何有效地模拟渗透系数 变化后的真实情况的问题。 由于有限元方法的对象是假定的连续介质，但是，由于研究管涌等渗透变形 得具体形式的破坏过程时，介质作为连续介质是不合理的，因此有限单元法在这 一领域并不能很好地发挥作用。也即运用有限元方法来进行管涌现象的研究在原 理上是值得探讨的。因而，理想的方法是将有限单元取无限小，单个颗粒作为一 个分析单元最能模拟实际的情况。总的来说管涌的临界水力梯度的理论研究至今 尚不成熟，目前只能在宏观上对其进行研究探索，管涌的微观渗透机理还没有在 理论上突破，并且实验手段也不能很好的模拟这种随机性比较强的状况。因此对 管涌机理的研究具有一定的挑战性。但是可以看到，越来越多的研究者都已经将 注意力放在从微观上揭示管涌的机理。 第三节本文的研究内容 本文在广泛调研国内外文献资料的基础上，拟对无粘性土中的管涌机理进行一定 的研究和探讨。主要内容如下： 首先，从管涌的特点出发，在分析比较现有的数值计算方法的基础上，选 用离散单元法对管涌的机理进行尝试性的研究，在学习了离散元的相关原理后， 。针对将离敞单兀法运用卡管涌研究中会遇到的问题进行了定的探讨。 第一章绪论 第二，根据流体动力学的基本原理，并且在考虑了多孔介质特性的基础上， 将流体动力学的基本方程进行扩展，在流体动力学方程中，考虑将多孔介质对流 体的阻力项作为一项原力考虑，进而进行流体与多孔介质的耦合计算。 第三，在分析了离散元相关软件的基础上，选用颗粒流软件p f c 3 d ( p a r t i c l e n o wc o d e ) 对管涌的发展过程进行数值模拟，首先尝试进行颗粒间级配组合即几 何条件对管涌的影响，采用粗细两种颗粒对管涌现象进行模拟，然后分析颗粒组 成对流场中速度场压力场的影响，进而得出计算域内水力梯度分布，计算出流域 内每个颗粒实际受到的拖曳力，研究细颗粒受到的阻力，得出需要克服的起动阻 力，和受到的拖曳力进行比较。在分析颗粒起动的临界拖曳力的时候，将颗粒的 临界拖曳力分布和流场中实际施加在颗粒上的拖曳力进行比较，宏观上按照两者 的重合量作为判断流场中颗粒流失多少的依据，也可以直接将计算单元内颗粒的 流失量作为判断管涌发生与否的标准。 第四，在颗粒组成一定的基础上，采用p f c 研究影响管涌的条件，分别探 讨了渗径长度，截面尺寸以及渗流坡降升高速度对管涌的影响，得出了相关的曲 线，和前人得出的结果比较接近，从而也说明了将p f c 程序运用于管涌的研究 是值得尝试的。 第二章离散元法及散粒体特性概述 第二章离散元法及散粒体特性概述 第一节d e m 的基本原理 离散元法一d e m ( t h ed i s t i n c te l e m e n tm e t h o d ) 是7 0 年代才发展起来的适用与 散体力学数值模拟的计算方法，主要研究离散体的运动变形，接触颗粒间的相互 作用等，适用于模拟离散颗粒在拟静态或动态条件下的变化过程，考虑到散体的 离散特性，研究方法与传统的连续介质模型分析或模拟完全不同，而是以固体颗 粒的运动学和接触力学理论建立数值模型。 离散单元法的理论基础【1 3 】1 4 】是结合不同本构关系的牛顿第二定律；有限单元 法的理论基础是基于最小势能原理；边界单元法的基础是b e t l i 互等定理。连续 介质问题求解时，除满足边界条件外，还必须满足三大方程即：平衡方程、变形 协调方程、本构方程。离散单元法在处理离散体的时候认为离散体是( 块、盘、 球) 的集合，故块与块之间没有变形协调的约束，继而采用显式差分法进行求解。 2 1 1 离散单元法的发展及研究现状 离散单元方法一般首先由c u n d a u 【1 5 】于1 9 7 1 年首先提出，最初是用来分析 岩石边坡的运动。1 9 7 4 年二维的离散单元法程序趋于成熟，在以后的阶段c u n d a u 和s t r a c k 开发了二维圆形块体的b 创一l 程序，用于研究颗粒介质的力学行为， 令人欣慰的是计算结果和d 诧s c h e r 等人用光弹性技术得出的实验结果极为吻合， 这也就说明了离散元方法由理论到实践的可行性。1 9 8 0 年c u n d a l l 开始研究块体 在受力后变形以及根据破坏准则允许断裂的离散单元法，后来发展了u d e c 程 序，三维离散单元法的开发在u d e c 之后，基本原理和u d e c 相同，称为3 d e c ， 在前面的基础上，结合岩土工程中土体的特点，玎a s c a 公司开发出了针对土体 颗粒的p f c 程序，用来模拟外力作用下圆颗粒的流动特性。 离散单元法最早由王咏嘉和剑万禧与1 9 8 6 年介绍给国内的学者，并在这方 面做了相当程度的研究，特别是将离散元方法运用与放矿工程【1 6 】取得了相当多的 成果；i 中国农业大学的徐泳、黄文彬1 17 j 将其运用于模拟颗粒流动，并对颗粒中充 填流体的特性进行了模拟，基于液桥力建立了一些湿颗粒模型；中科院的李世海、 第二章离散元法及散粒体特性概述 汪远年对于离散元方法中的参数选掣1 8 j 作了深入的探讨，并将离散单元方法运用 来进行断续节理岩体随机模型三维离散元数值模拟【1 9 j ；俞良群等对料仓堆积中的 现象进行了实验乜踟及数值模拟心嵋的对比：同济大学的周健教授采用离散元中的 p f c 软件进行了针对拔桩实验【2 2 】以及对沙土的细观力学性质【2 3 】进行了模拟，取 得了相当的成果；同时其学生也运用离散单元法对渗流的颗粒流细观模拟1 2 4 j 中国 农业大学的李红艳瞄】徐春晖【2 6 】对存在填隙流体时颗粒间的相互作用，杨全文【2 7 】 针对离散元法干颗粒接触模型及微机可视化程序设计进行了相关的研究。同时国 外的学者运用离散元方法对各种情况进行了大量的模拟计算，取得了较好的效 果，如j l i 木，c w 色b 【2 8 】等对筛料的过程运用离散元的方法进行数值模拟，但是对 于考虑流体作用的研究的相对于干颗粒模型还是较少，t n o u c h i 考虑浮力作用 在颗粒流的实验及数值模拟中的影响【2 9 1 。 2 1 2 离散单元法的接触模型 接触模型是颗粒离散元法的核心，在颗粒流中有硬球、软球两种基本模型。 硬球模型是假定颗粒流中颗粒表面承受的应力较低时，颗粒之间不发生显著的塑 性变形，认为颗粒之间的碰撞是即时的，在碰撞过程中颗粒本身不会变形，而且 其缺点是只考虑两个颗粒之间的同时碰撞，因此只能用于稀疏的快速颗粒流中； 软球模型和硬球模型相对，假设颗粒碰撞能够持续一定的时间，同时可以考虑多 个颗粒的碰撞，因此具有较广的适用范围。 软球模型中假定颗粒发生碰撞保持彤状1 i 变。而且互相叠加。在现阶段很 多学者都采用弹簧一阻尼器模型【1 7 】。c u n d a u 首先提出了盘元接触模型，法向的 作用简化为一个弹簧一阻尼器，切相作用则简化为弹簧一阻尼器加一个滑动摩擦 器。允许切向滑移和法向的分离，但是实际上的圆盘颗粒极少，若法向和切向刚 度是根据三维颗粒实测得到的，则模型就需要采用三维球元模型。三维球元模型 除采用弹簧一阻尼器( 见图2 2 ) 给定参数外，也有采用接触力学结果的模型， 法向作用根据h e r t z 理论，切向作用根据m i n d l i n 和d e r e c i e w i c z 理论确定，当然 后者法向和切向的力位移关系都是非线性的。 一 。： 一 图2 1 中b 为接触区直径，“。为法向叠合量，9 为接触中心点，图2 2 中尼。、 ： - ：? ： j 、；： ：虹、d d ，、d ，。分别为法向接触刚度、切向接触剐度、接触法向阻届、切向阻 河海大学硕：i ：学位论文 两个相互接触的圆颗粒 图2 1 球元模型图2 2d e m 模型 尼，和质量阻尼，各项阻尼作用将在后面详细介绍。 在模拟颗粒实际接触时，需要做如下假定： ( 1 ) 颗粒单元为刚性体； ( 2 ) 接触发生在很小的范围内，即可以近视认为是点与点接触； ( 3 ) 接触特性为柔性接触，接触处允许有一定的叠合量； ( 4 ) 叠合量的大小与接触力有关，与颗粒大小相比，重叠量很小； ( 5 ) 根据模拟情况不同，接触处连接强度不同； ( 6 ) 颗粒单元采用圆球 2 1 3 离散单元法中的基本方程 ( 1 ) 物理方程一一力和位移的关系 f 。= k n 。 4 f s = ks 4 us 式中七。、七，分别表示法向和切向刚度 ( 2 ) 运动方程一一牛顿第二运动定律 ，：墨 优 f ，。生 肌 u i ( tl ) = di ( ta ) + ot 4 t 。 u i ( t1 ) = ui i tq ) + u1 4 t “；臼f ( + j 掣心 式中u ，矽j 矽；分别为位移j 速度、加速度，口分别为角位移， ( 2 - 1 1 ) ( 2 - 1 - 2 ) f j 为角加速度，聊 第二章离散元法及散粒体特性概述 为单元质量，r o 为起始时间，彳f 为时步；f 1 = f o + 彳f ，f = x ，y ，z 。 当然在实际颗粒计算中针对不同的情况，颗粒受到的力和力矩会不一样j 但 是其基本求解就是在力与位移之家进行不断的转换计算，受到别的力的时候将其 作为单元受到的力代入外力，中即可。 2 1 4 离散单元求解方法 离散单元法中有两种求解方法：静态松弛法和动态松弛法。 在离散元方法中，颗粒之间的作用是拟静态的，但是整个系统的运动是动态 的，在每一时步对所有的颗粒单元进行检查和迭代运算，分析颗粒可能在下一时 步和其它颗粒接触的情况，确定每个单元的接触状态，进而得出整个系统的运动 状态。动态松弛法是把非线性静力问题化为动力学问题求解的一种数值方法，其 实质是对临界阻尼振动方程进行积分。为了保证求得精确解，一般采用质量阻尼 和刚度阻尼来吸收系统的动能。将带有阻尼项的动态平衡方程，利用有限差分法 按时步在计算机上进行迭代求解的方法就是动态松弛法。其最大的优势就是可以 考虑非静态问题。用动态松弛法求解的基本运动方程为： m 【厂o ) + c u 0 ) + 足u0 ) = ，0 )( 2 1 - 4 ) 式中m 为单元的质量；u 是位移；t 是时间；c 是粘性阻尼系数；k 是刚度系 数；厂是单元素所受的外荷载。利用中心差分法，在假定f + 彳f 时刻以前的变量 只砂、“ 、矽0 一小) 、d o 一4 f ) 以及【厂o 一4 f ) 等已知的前提下得到单元在f 时刻 的速度和加速度： d ( t ) 二 u ( t + t ) 一u ( t 一t ) ( 2 t ) o ( t ) ； u ( t + t ) 2 u ( t ) + u ( t t ) ( t ) 2 八 ( 位移u ) u ( 2 1 - 5 ) ( 2 - 1 6 ) 图2 3 离散元计算中的力位移循环过程 _ ， ：一 ， l 得到f 时刻j 的速度，角速度后，可以得到颗粒在下一时刻的新的位置，由各球的 1+ 。 ：、：tt ： 、 ： ：t 一1 1 1 河海大学硕士学位论文 新的位置坐标可以确定相邻颗粒是否接触或原接触点是否脱离，相互接触的球会 产生图2 1 中的假性重叠量，再由接触模型求出接触力和力矩。总之通过计算循 环( 如图2 3 ) 可以求解任意时刻的力、位移、加速度。 2 1 5 计算参数的选取 ( 1 ) 阻尼的选取 在颗粒的碰撞过程中，运动的颗粒最终必然要趋于稳定，也即动能转化为 热能，因而系统中单元的运动是一不可逆的过程为了耗散动能，离散元计算中 引入阻尼来吸收系统的动能。阻尼系数的考虑有质量阻尼和刚度阻尼，两者的作 用分别为限制颗粒的绝对运动和相对运动，前者等价于将整个系统浸泡在粘性液 体中，在物理上可以理解为将块体单元和不动点连接起来，后者则可以理解为用 粘性活塞将接触单元连接起来，从而使单元的相对运动受到阻碍。质量阻尼产生 与物体的绝对速度成正比而方向相反的阻尼力，质量阻尼系数为： d 。；口m( 2 1 - 7 ) 刚度阻尼系数为： d 。= 肚。( 2 1 - 8 ) d ，= 户k ， ( 2 1 - 9 ) 式中，a ，卢为r a y l e i 曲阻尼常数，口；豺r 1 1 i 。a 。j n ，户；兄i i l i 。砑f i i i 。，式中厂m i 。二 系统基频，为系统的弹模、密度、泊松比及基于系统基频的波长的函数。a 。；。一 系统临界阻尼比，m 为单元质量。由系统宏观参数则大体可以确定阻尼各项参数。 阻尼系数产生与刚度成正比且与接触点处相对运动速度成正比的阻尼力。 而与对于做低频运动的系统，质量阻尼比较有效，而对于快速高频振动的系统， 刚度阻尼的作用则显得及其明显。阻尼的选取在时步确定的情况下会影响系统的 收敛程度，对于连续弹性介质来说，振型频率和阻尼比能够经过计算与实验获得， 相对容易得到，但是对于非连续介质，由于振型不确定，只能采取尝试的方法来 获得接近容易达到收敛状态时的临界阻尼。因为欠阻尼将使求得的解为非静态， 而过阻尼则使得达到平衡：状态所需迭代次数过多。临：界阻尼系数为： 第二章离散元法及散粒体特性概述 c 。：2 ；忑；2 小。式中m 为单元质量，七为弹性刚度系数，考= c c 。为阻尼比， 脯和志2 志艄哟黝椭弧孰尼o ( 2 ) 时步的选取 在离散元计算过程中，除了阻尼外，对系统的计算同样重要的是计算时步的 选取。由于系统采用差分法对微分运动方程求解，因而每一迭代时步的取值对求 解的稳定性有很大的影响，为了确保数值计算的稳定，一般都要球时步取非常小 的值。时步取值最主要的依据是颗粒在个时步内只和周围颗粒发生作用，对其 它颗粒的运动不产生影响，从而避免由于时步过大，系统内的力学性质重新调整， 丧失了离散元模拟颗粒任意时刻的位移、受力的优势。由于离散单元特别是本文 所采用的颗粒流方法来模拟管涌的过程中，用颗粒球来模拟实际的砂土颗粒，颗 粒本身的体积非常小，因而在计算中为了实际模拟颗粒间的作用，我们按照下面 的标准来选取时步：理论上单自由度的具有质量m 刚度七的弹性振动系统，在 运用中心差分进行计算时，得出出 m a x 址，因而我们可 以将川i n 出作为选取计算时步的标准。当然由于颗粒大小不一，由于有的部分颗 粒质量很小，按照上式计算则会使得计算时步取得非常小，从而影响计算收敛的 速度，可以利用虚拟密度的方法来进行调整。 虚拟密度的原理就是将我们数值上相等的质量m ；与惯性质量聊：加以区分， m ；豇= f + m ( 2 1 1 0 ) 式中f 是单元受到的除重力外的合外力。如惯性质量取系统内所有单元质量 的平均值j 则上式可以改写成： 历：詈+ 譬 ( 2 - 1 ) m ，m ， 、7 式中所，表示系统内所有单元的平均质量，从式2 2 2 可以得到i 各单元在相 - 。： ?i，；- - 等的外力作用下会有相同的加速摩，旦各尊元的重力仍签具有原来的作用效果。 由于时步由具有最小质量的单完疾走：卤此，在冈f j 焘一定的宿况下，遗过增 ：，j 。 一 ，j 河海大学硕士学位论文 加最小单元的惯性质量，增大时步，提高计算的收敛速度。 ( 3 ) 刚度及摩擦系数的选取 在用离散元模拟颗粒运移的过程中，需遵循一定的规则，在进行数值模拟时， 对于采用简单的模型模拟所得出的结果必需和由实验室所得的实际结果相对应， 这样的参数才是合理的，在采用离散元法进行计算时，当各种参数不是很齐备的 时候，可以通过采用简化模型进行数值模拟，通过数值模拟的方法来得出宏观效 应和实际实验相同的参数，对于本文的情况，选取不具有粘性的砂土颗粒进行模 拟，则模拟结果必需和三轴实验进行对比，在孔隙率一定的情况下，可以通过调 整接触杨氏弹模e 。( 不同于杨氏模量e ) 、法向刚度和切向刚度的比值尼。肛，以及 摩擦系数来达到和三轴相一致的结果，在进行模拟的同时发现：在微观层面泊 松比没有意义物理意义。 值得注意的是在离散元模拟颗粒的过程中会用到摩擦系数的概念，而且以摩 擦力作为判断颗粒之间是否在切向发生滑移的依据，形式上和c o u l o m b 破坏准 则相同，但是太沙基对讨论传统土压力的假定进行分析时，就提到c o u l o m b 导 入时存在着基本的错误，忽视了砂土由个别颗粒组成的事实，而将其看作具有一 定力学性质的均匀物体，也即用表征体积进行考虑，从而在宏观上研究土体的时 候选用的就是代表土体宏观特性的参数，而没有考虑颗粒间的接触方式。从图 2 4 中我们可以看到，在散粒体中，颗粒之间在剪切面上存在三中可能情况，一 是产生接触，二是保持接触，再就是接触消失，同时可以看到，颗粒之间的接触 面和结构剪切面的方向并不相同，所以实际颗粒之间的摩擦系数并不能用简单地 认为就是宏观表现出的值，而应该是根据颗粒实际的表面情况，以及其它影响因 素确定，在实际研究散粒体受力的过程中，颗粒间的内摩擦系数除了和颗粒表面 有关，还和颗粒受到的中间应力、有效应力等因素相关。 一么：” r 与强x 蝣! 弓姆镂丢罗磷沁9 ：曼事季p 。接触 i j ；觐 f ； 扮接触 。一铲- 、- 参骞1 傣登登