水 利 学 报
2002年12月
文章编号:0559-9350(2002) 12-0011-07
*SHUILI XUEBAO 第12期粘性土室内平面应变试验的颗粒流模拟
廖雄华, 周健, 徐建平, 林利敏
31浙江财经学院, 浙江杭州 310012) 1123(11同济大学地下建筑与工程系, 上海 200092; 21武汉轨道交通有限公司总工办, 湖北武汉 430017;
摘要:作为颗粒流理论实际应用的前期可行性研究, 针对具有凝聚力的内摩擦材料粘性土进行了室内平面应变试验结果的颗粒流仿真。对三维颗粒流数值仿真技术(PFC 2D ) 的4组数值仿真结果和粘土试样室内平面应变试验的实测曲线进行了对比, 研究表明利用颗流理论所建立起来的PFC 2D 数值仿真试验模型是能够通过改变计算模型中颗粒单元的性质, 以及颗粒集合体的级配特征等给出与真实材料土工试验类似的本构行为的, 这种可行性最根本的意义在于基于颗粒流理论的数值仿真试验能够突破常规土工试验在仪器设备能力、试验条件上的局限性, 对实际土样的本构行为做出理论方面的预测。此外针对颗粒流仿真试样的细观力学特征与物理试样宏观力学响应之间的关系进行了参数研究, 并揭示了一些规律。
关键词:颗粒流理论; 粘土; 平面应变试验; 本构行为; 应力-应变曲线
中图分类号:TU411文献标识码:A
目前测定土抗剪强度的土工试验方法主要有直剪试验、常规(轴对称) 三轴试验以及平面应变试验等。由于仪器和试验方法固有的缺陷, 在土的抗剪强度的测定上存在很多不足, 如直剪试验的剪切破坏面人为固定并非最薄弱面, 同时不能控制排水条件; 常规三轴试验虽然克服了直剪试验的缺点却采用轴对称的应力状态, 而实际上工程问题多为平面应变或三维空间问题。即使对于平面应变试验或真三轴试验也还是存在一些由于仪器设备、试验条件本身所产生的强度测试误差问题。
在数字信息技术高速发展的背景下, 作为对土工试验技术的一种补充, 开展土工试验的数值仿真研究十分必要。颗粒流理论及其数值仿真技术(如PFC) 作为一种特殊的离散单元法其计算单元都是刚性的且单元的几何性态都为纯粹的球形或圆盘形, 尤其适用于散粒介质的力学分析。作为颗粒流理论实际应用的可行性研究, 对岩土材料室内试验结果的仿真是一个前期基础性的研究。文献[1]曾对砂土的平面应变试验进行了颗粒流模拟, 研究结果表明颗粒流方法对于无凝聚力的内摩擦材料砂性土类的室内试验应力-应变关系模拟是可行的。论文针对具有凝聚力的内摩擦材料粘性土类进行试样的室内平面应变试验的颗粒流仿真。
1 粘性土平面应变试验的颗粒流仿真过程
111 颗粒流试样的准备 所选择的粘性土为上海地区原状土, 为第º层褐黄色粉质粘土, 其物理力学性质指标以及试验结果可参见文献[2, 3]。
11111 计算参数的选择 根据文献[2, 3]的试验条件, 试样的尺寸取为70mm @25m m 。为了更好地逼近原土样在微观上的各向异性和不均匀性, 在生成二维颗粒流仿真模型(PFC ) 试样时设定颗粒收稿日期:2001-08-20
基金项目:上海市博士后科学基金沪科2000(478号) 资助; 国家自然科学基金项目资助(50178054)
作者简介:廖雄华(1969-) , 男, 安徽芜湖人, 博士, 现为同济大学地下建筑与工程系博士后, 主要从事岩土工程数值分析和
散体介质力学方面的研究工作。2D
图1 粘土室内平面应变试验的颗粒流PFC 模拟流程
试样是由不同半径的颗粒单元所组成, 颗粒半径R 的分布采用从R min 到R max 的均匀分布。经过大量颗粒流仿真方法PFC 试样的仿真试验, 对于文献[2]的试验情况仿真试样均取R min =5mm, R max /R min =217。
在PFC 程序中, 对于有内部凝聚力的材料是无法给定孔隙率的, 程序将根据R min 和R max /R min 自动设定。但由于PFC 所采用的颗粒均为球形颗粒, 且每个颗粒都与其周围其它颗粒相接触, 因此可以通过删除部分颗粒单元的办法来增大其孔隙比, 另外如果在PFC 试样中采用不均匀颗粒的组合方式, 则其孔隙比还将有更复杂的变化。
PFC 颗粒试样是通过颗粒间的相互作用来表达整个宏观试样的应力响应的, 对于粘性土颗粒试样PFC 所定义的颗粒间相互作用有:在接触处的法向接触力F ij 、切向接触力F ij , 接触摩擦力F ij , 下标表示力由第i 个颗粒单元通过接触作用于第j 个颗粒单元上; 根据颗粒流理论的接触本构假设, 上述各个相互作用力分量可分别通过法向刚度k ij 、切向刚度k ij 和摩擦系数L ij 和法向相对位移U ij 和切向相对位移U ij 按下式计算:
F ij =k ij U ij , F ij =k ij U ij , F ij =L ij F ij =L ij k ij U ij n s s s s s f n n n s [4]n s n 2D n s f 2D [4]2D 2D (1)
2D
n 法向接触力沿两颗粒单元圆心的连线, 切向接触力与摩擦力则与之垂直。对于粘性土颗粒试样PFC 给定颗粒之间的接触模量E c 以及定义法向刚度、切向刚度之比k n /k s , 在文中的计算中k ij =k , L ij =
L 。另外对于粘性土还需在PFC 试样的颗粒间定义法向、切向连接强度R c 和S c 以及各自的标准偏差。颗粒的比重可根据实际土样的干容重取为1189kN P m 。
表1 PFC 2D 数值试样的基本输入参数
围压R 3
/kPa
25
100
150
200试样尺寸/mm(高@宽) 70@2570@2570@2570@25粒径R min /mm[1**********]5
2D 3R min R max [1**********]7摩擦系数[**************]0颗粒法向接触刚度E c /Pa 4135@1064135@1064135@1064135@106颗粒刚度比[**************]0接触强度接触强墙体速率/Pa 113@104113@104113@104113@104度偏差[**************]5/mm [1**********]10收敛控制[**************]0初始锁定应力/Pa 75@10375@10375@103100@103从表1可以看出, PFC 颗粒流数值仿真试样的基本性质数据是相同的, 试样在不同围压下的仿真数
据差别主要在于颗粒试样的摩擦系数、试样加载过程中最终的墙体移动速率、试样初始的锁定应力状
|R y -R x |态, 以及加载试样最终收敛时的控制应力比, [A 。(R d ) max
11112 PFC 试样的颗粒生成过程 首先定义墙体, 共4道, 其包围的矩形为70m @25mm, 见图2(a) 。在PFC 生成墙体的过程中需要控制:(1) 上下墙体与试样宽度之比以及两侧墙体与试样高度之比要适当; (2) 颗粒与上下墙体以及两侧墙体间的接触刚度同颗粒自身间的法向接触刚度之比要适当。否则, 可能会出现颗粒在墙体运动的加载过程中由于没有墙体的约束而溢出、作自由运动的现象, 或穿透墙体而向四周散逸, 导致平衡迭代发散、实验失败。接着PFC 在上述给定的墙体范围内生成颗粒, 把半径在指定区间[R min , R max ]的颗粒往区域内填充, 如果没有已生成颗粒与之重叠, 则生成此颗粒, 否则, 改变颗粒的位置重试, 最后通过循环来消除试样内部非均匀应力。与无粘性的砂土颗粒试样不同, 粘性土颗粒试样的孔隙比是不能人为给定的, 而是由程序自动设定一个值如n =16%, 由此可以近似计算出试样的颗粒数目
N =[4]2D 2D 2D :(2) R min +R max , R =2式中:b 为试样的宽度; h 为试样的高度; N 为颗粒单元总数。
图2给出了用颗粒流程序模拟室内双轴试验整个试样的颗粒分布情况。图2(b) 为生成小粒径颗粒后颗粒的分布情况, 通过半径的膨胀, 颗粒分布情况见图2(c) 。半径膨胀后, 颗粒间将出现非平衡力, 通过循环让试样颗粒自动达到平衡状态(为加速达到平衡状态, 可以考虑增大颗粒间的摩擦力) , 此时试样内部颗粒分布情况见图2(d) 。
此外, 对于粘性颗粒试样, 为了消除试样内部可能存在的初始应力集中的不合理现象, 需要在试
2D 样加载之前先对试样施加一个初始的锁定应力(PFC 的模型参数) , 并在该应力作用下使试样达到初
始的均匀应力状态, R , 如x =R y =R 0, 一般R 0的取值大小相对于试样材料的峰值强度应该足够小小于单轴强度的1%。
经过上述工作基本上建立起一个能够和实际粘土试样比拟的PFC 数值仿真试样。
112 颗粒流数值仿真试验的实现
2D 11211 载荷与围压的施加和保持 PFC 是通过一套数值伺服系统让顶部和底部墙体作相对运动来施
加试样荷载的, 并同时调整两侧墙体的位移, 以保持试样的围压R x 恒定, 加载后的试样颗粒分布情况如图2(e)
所示。2D [4]
(a)墙体 (b) 墙体+颗粒(未膨胀) (c) 颗粒膨胀后的试样 (d) 平衡后颗粒试样
2D (e) 加载后的颗粒试样图2 室内双轴试验的颗粒流模拟11212 应力和应变的计算 在整个加载、卸载和再加载的过程中, PFC 程序可以自动记录试样的水
平E x 和竖向应变E y , 以及R x 和R y 。R x 和R y 由墙体和颗粒的接触力的总和除以颗粒和墙体接触的面积, 即颗粒试样的侧面面积和顶底部面积来计算:
R x =f x
2 r h , R y =f y
2 r b
(3)
式中:f x 为左侧或右侧墙体与颗粒单元的x 方向接触力, f y 为顶部或底部墙体与颗粒单元的y 方向接触力, r 为颗粒球体的平均半径, 在PFC 中则 r 取单元厚度。应变按下式计算:
E x =; E y =b h (4) 2D
式中:$b 和$h 分别为试样的x 和y 方向变形, 可以由墙体的位置来计算。
以(R y -R x ) 和E y 的记录作图便可得到试样的应力-应变关系曲线, 此即(R 1-R 3) ~E a 关系曲线。
2 粘土本构行为的颗粒流仿真研究
由于颗粒流试样实际上都是基于散粒介质本身微观特性基础上建立的, 故PFC 数值模型内颗粒间所发生的相互作用反映了散粒介质骨架的细观力学行为, 即基于颗粒流试样所模拟出的粘土介质宏观应力应变响应特性其实就是土体骨架所具有的本构行为。在文中的加载方式下, 可以不考虑瞬态流-固耦合问题, 颗粒流仿真试样得出的应力-应变曲线就是有效应力与应变的关系曲线。但在更深入的颗粒流试样地震荷载的响应特性分析, 则进一步的多相颗粒流理论研究是必需的。
211 轴向应力-应变关系曲线的PFC 仿真结果 以文献[2, 3]上海粘土固结不排水平面应变实验的部分数据结果为仿真对象, 给出4种围压下的应力-应变(R 1-R 3) ~E a 曲线对比。从图3围压为25kPa 的应力-应变曲线的比较来看, PFC 的数值仿真结果与实际试验实测曲线吻合得很好。但从围压分别为100kPa 、150kPa 、250kPa 的应力-应变曲线的比较还可看出, 细观颗粒各项指标相同的试样却无法得到较为理想的仿真结果, 见图4~图6。这说明, 颗粒流试样尚不能很好地反应物理试样客观上的唯一性, 不同的颗粒试样得到相似仿真曲线的数值试验从问题的另一个角度更进一步表明了该
点。2D 2D
图3 围压为25kPa 时的应力-
应变曲线对比图4 围压为100kPa 时的应力-
应变曲线对比
图5 围压为150kPa 时的应力-
应变曲线对比图6 围压为250kPa 时的应力-应变曲线对比
从见图3~图6的粘土试样实验曲线还能看出, 试样在曲线的起始阶段表现出一定的初始强度(初始切线模量较大) 。这实际上是由于粘土试样内部的细观结构特性而导致的。而对于颗粒流仿真试
样来说, 由于仅采用了单一的圆形颗粒单元来模拟整个结构, 因此无法揭示真实粘土试样内部的细观结构特性。显然, 该方面进一步的研究是必要的, 如果采用多个圆形单元组合成一系列不同形状的组合颗粒块来模拟, 有望得到更好的结果。
212 颗粒流模型若干参数对粘土本构行为影响的仿真研究 所谓的颗粒流数值仿真试样其实就是通过不断改变颗粒单元及其集合体组构方面的细观力学性质, 从而使之逼近真实材料的宏观力学响应。因此, 颗粒流理论对岩土试样的数值仿真过程就是一个不断修改、调试各项细观力学参数(包括颗粒粒径) 的过程, 现就此做若干参数研究如下(假定围压为25kPa) 。
21211 摩擦系数对试样本构行为的影响 图7给出颗粒间摩擦系数在011~019, 9组颗粒流试样应力-应变本构行为的仿真曲线。从不同摩擦系数应力-应变关系曲线的对比可以发现, 在其它因素不变的情况下随摩擦系数的增加, 应力-应变关系曲线的峰值会提高, 并且初始弹性模量也有一定的提高。从图中还可以发现, 最大轴向允许应变随着摩擦系数的增加而提高。另外, 摩擦系数越大, 软化现象越明显, 且有小的摩擦系数对应理想弹塑性关系, 而提高摩擦系数导致软化现象。
21212 加载速率(墙体移动速度) 对试样应力-应变曲线的影响 从图8中可以看出, 颗粒流试验仿真过程中墙体的移动速率(加载速率) 对仿真试样的应力-应变曲线影响不大, 尤其是对峰值强度以前的曲线部分几乎没有任何影响。而在峰值强度以后对渐进破坏部分的曲线以及允许应变有些程度不大的影响, 并且对于一组确定的颗粒流试样, 似乎存在一个最优的墙体移动速率(加载速率
)
图7 摩擦系数对仿真试样应力-应变曲线的影响
图8 墙体移动速率对仿真试样
应力-应变仿真曲线的影响
21213 颗粒接触刚度对试样应力-应变曲线的影响 从图9可以看出, 颗粒间接触刚度E c 对仿真试样的应力-应变曲线有较大的影响。显然, 随着E c 的增加, 试样的初始弹性模量增大, 其出现峰值强度对应的应变将降低; E c 对试样的峰值强度有影响, 但并不大。另外颗粒间接触刚度的大小对仿真试样应力-应变曲线的峰后形态以及允许应变也有一定的影响, 但并不象颗粒间摩擦系数那样具有
本质性的影响。
图9 颗粒接触刚度对仿真试样应力-应变仿真曲线的影响
21214 最大与最小颗粒粒径之比对试样应力-应变曲线的影响 仿真试样的最大与最小颗粒粒径比(与试样的孔隙比程函数关系) , 从图10可以看出, 对仿真试样应力-应变曲线的峰值强度及其对应的应变, 以及允许应变都有较大的影响。比值越大峰值强度越高、软化现象也越明显; 而对于允许应变来说, 最大与最小颗粒粒径比的影响是非线性的, 在其它因素确定的情况下同样存在一个最优值。
图10 最大与最小颗粒粒径之比对颗粒流应力-应变仿真曲线的影响
3 结论和前景展望
以文献[2]的工作为依据, 利用颗粒流计算机仿真程序对文献[2]的粘土平面应变试验结果进行了部分仿真, 并针对颗粒流试样的微观参数进行了多组试样的参数研究, 关于岩土试样宏观力学响应问题的颗粒流细观力学途径方面取得了以下研究结论, 如表2所示。
表2 岩土试样宏观力学响应PFC 2D 数值仿真的参数研究结果
初始弹性模量
摩擦系数f
墙体移动速率v
R m ax /R min
颗粒接触刚度E c 增大(小幅) 无影响无影响增大(大幅) 理想塑性y 软化趋势明显无影响有趋势影响不大最大轴向允许应变降低(大幅) 有非线性影响有非线性影响降低(小幅) 峰值强度增大(大幅) 无影响增大(大幅) 影响不大峰值强度对应的应变无影响无影响有非线性影响降低(大幅)
此外, 作者的研究还揭示出, 岩土试样本构响应问题的颗粒流仿真研究中存在的部分问题, 例如颗粒流细观力学模型的非唯一性, 以及粘土试样固结问题的颗粒流仿真途径还尚待解决, 这包括流固耦合、孔压和排水边界条件等若干细节的颗粒流仿真。
最后值得指出的是, 用颗粒流理论与方法来模拟岩土类散粒介质宏观力学响应的数值途径, 正如文中所介绍的, 是基于细观层面上所有散粒之间的机械类相互作用, 通过接触力、摩擦系数和接触强度来表现颗粒间的接触应力、内摩擦特性和凝聚力。之所以说是细观而非微观(如颗粒直径达到01001mm 的尺度) 主要是因为:(1) 模拟的颗粒直径过小时, 计算量将大幅增加, 仿真受到软硬件条件的限制而不可行; (2) 目前的颗粒流理论和方法还不能充分考虑粘粒效应, 如水化膜、范德华力等物化作用。另外, 由于尚未完全研究清楚颗粒流数值仿真试样与粘土物理试样之间的颗分对应关系, 所以在模拟过程中还具有一定多耙试探性(trial and validation) , 比如本文就取粒径015~1135mm 。可以明确的是, 不同粒径范围的散粒集合体的宏观力学响应是不同的, 这正如粗粒土与细粒土的力学性质差异。综上所述, 作者的工作还只能是一个初步的可行性研究, 该研究还有待深入。
参 考 文 献:
[1] 周健, 池毓蔚, 池永, 徐建平. 砂土双轴试验的颗粒流模拟[J].岩土工程学报, 2001, (6) :701-704.
[2] 董建国, 李蓓, 袁聚云, 赵锡宏. 上海浅层褐黄色粉质粘土剪切带形成的试验研究[J]. 岩土工程学报,
2001, (1) :23-27.
[3] 李蓓. 平面应变仪的研制及上海地区粘性土剪切带形成的试验研究[D].上海:同济大学, 1999.
[4] Itasca Consulting Group, Inc.. PFC 2D User . s M anual (Version 210) [Z]. Minnesota, USA, 1999.
[5] Cundall P A, Strack O D L. A Discrete Numerical Model for Graunlar Assemblies. Geotechnique, 1979, 29(1) :47-65.
[6] Huang H, Detournay E, Bellier B. Discrete element modelling of rock cutting [J]. Rock Mechanics for Industry, 1999,
1(1) :123-130.
[7] 冯国栋. 土力学[M]. 北京:水利电力出版社, 1986.
[8] 周健, 池永, 池毓蔚, 徐建平. 颗粒流方法及PFC 2D 程序[J].岩土力学, 2000, (3) :271-274.
Simulation of plane strain test of clay by means of particle flow code
LI AO Xiong -hua , ZHOU Jian , XU Jian -ping , LIN L-i min
(11Tongji Universit y , Shanghai 200092, China ;
21Gene ral O ffice o f Rail way Traffic Company Ltd , Wuhan City , W uhan 430017, China ;
31Zhe jian finance and e conomics college , Hangzhou 310012, China ) 1123
Abstract :The simulation of plane strain test for clay using particle flow code method is carried out. The comparison of simulation with test result shows that the model based on particle flow theory can reflect the constitutive behaviors of soil by means of changing the particle . s properties and size grading characteris -tics during modeling. The model can be used to predict the constitutive behaviors of clay on the level of theoretical. The relationship between microscopic para meter of 2-D particle flow and mechanical response of physical sample are also investigated. Some basic laws are found.
Key words :particle flow theory; clay; plane strain test; constitutive behavior; stress -strain curve (上接第10页)
参 考 文 献:
[1] 李锦秀, 等. 三峡库区江段纵向离散系数研究[J]. 水利学报, 2000, (8) :84-86.
[2] EPA/600/3-87/007. The Enhanced Stream Water Quali ty Models QUAL2E and QUAL2E -UNCAS [M]. Docu menta -
ti on and user Manual, 1987.
[3] 余常昭, M 马尔柯夫斯基, 李玉梁著. 水环境中污染物扩散输移原理与水质模型[M ].中国环境科学出
版社, 1989.
[4] R M Wrigh t, A J Mcdonell. In -stream Deoxygenation rate prediction [J]. Proc. ASCE. J. Env. Div. 105, 1979, 323
-333.
[5] 傅国伟编著. 河流水质数学模型及其模拟计算[M]. 北京:中国环境科学出版社, 1987.
Numerical simu lation of water quality for the Three Gorges Project Reservoir
LI Jin -xiu , LI AO Wen -geng , HUANG Zhen -li 112
(11China Institute o f W ate r Resourc es and Hydropowe r Rese arch , Bei jing 100083, China ;
21Executive O ffic e o f State Council for Thre e Gorges Projec t Construction Committe e , Be ijing 100038, China )
Abstract :The 1-D unsteady model based on finite difference method with four -point implicit scheme simulating the water flow and water quality in the reservoir of Three Gorges Project is established. The major indices of water quality including 14items are investigated. The effects of significant change of flow condition in the reservoir due to impounding are taken into account. The coefficients of the model includ -ing the longitudinal dispersion coefficient, decay rate of B OD 5and atmospheric re -aeration rate of DO etc. are established empirically. The model is applied to predic t the variation trend of water quality in the res -ervoir.
Key words :Three Gorges Project; water quality model; coefficients of water quality; empirical formula
水 利 学 报
2002年12月
文章编号:0559-9350(2002) 12-0011-07
*SHUILI XUEBAO 第12期粘性土室内平面应变试验的颗粒流模拟
廖雄华, 周健, 徐建平, 林利敏
31浙江财经学院, 浙江杭州 310012) 1123(11同济大学地下建筑与工程系, 上海 200092; 21武汉轨道交通有限公司总工办, 湖北武汉 430017;
摘要:作为颗粒流理论实际应用的前期可行性研究, 针对具有凝聚力的内摩擦材料粘性土进行了室内平面应变试验结果的颗粒流仿真。对三维颗粒流数值仿真技术(PFC 2D ) 的4组数值仿真结果和粘土试样室内平面应变试验的实测曲线进行了对比, 研究表明利用颗流理论所建立起来的PFC 2D 数值仿真试验模型是能够通过改变计算模型中颗粒单元的性质, 以及颗粒集合体的级配特征等给出与真实材料土工试验类似的本构行为的, 这种可行性最根本的意义在于基于颗粒流理论的数值仿真试验能够突破常规土工试验在仪器设备能力、试验条件上的局限性, 对实际土样的本构行为做出理论方面的预测。此外针对颗粒流仿真试样的细观力学特征与物理试样宏观力学响应之间的关系进行了参数研究, 并揭示了一些规律。
关键词:颗粒流理论; 粘土; 平面应变试验; 本构行为; 应力-应变曲线
中图分类号:TU411文献标识码:A
目前测定土抗剪强度的土工试验方法主要有直剪试验、常规(轴对称) 三轴试验以及平面应变试验等。由于仪器和试验方法固有的缺陷, 在土的抗剪强度的测定上存在很多不足, 如直剪试验的剪切破坏面人为固定并非最薄弱面, 同时不能控制排水条件; 常规三轴试验虽然克服了直剪试验的缺点却采用轴对称的应力状态, 而实际上工程问题多为平面应变或三维空间问题。即使对于平面应变试验或真三轴试验也还是存在一些由于仪器设备、试验条件本身所产生的强度测试误差问题。
在数字信息技术高速发展的背景下, 作为对土工试验技术的一种补充, 开展土工试验的数值仿真研究十分必要。颗粒流理论及其数值仿真技术(如PFC) 作为一种特殊的离散单元法其计算单元都是刚性的且单元的几何性态都为纯粹的球形或圆盘形, 尤其适用于散粒介质的力学分析。作为颗粒流理论实际应用的可行性研究, 对岩土材料室内试验结果的仿真是一个前期基础性的研究。文献[1]曾对砂土的平面应变试验进行了颗粒流模拟, 研究结果表明颗粒流方法对于无凝聚力的内摩擦材料砂性土类的室内试验应力-应变关系模拟是可行的。论文针对具有凝聚力的内摩擦材料粘性土类进行试样的室内平面应变试验的颗粒流仿真。
1 粘性土平面应变试验的颗粒流仿真过程
111 颗粒流试样的准备 所选择的粘性土为上海地区原状土, 为第º层褐黄色粉质粘土, 其物理力学性质指标以及试验结果可参见文献[2, 3]。
11111 计算参数的选择 根据文献[2, 3]的试验条件, 试样的尺寸取为70mm @25m m 。为了更好地逼近原土样在微观上的各向异性和不均匀性, 在生成二维颗粒流仿真模型(PFC ) 试样时设定颗粒收稿日期:2001-08-20
基金项目:上海市博士后科学基金沪科2000(478号) 资助; 国家自然科学基金项目资助(50178054)
作者简介:廖雄华(1969-) , 男, 安徽芜湖人, 博士, 现为同济大学地下建筑与工程系博士后, 主要从事岩土工程数值分析和
散体介质力学方面的研究工作。2D
图1 粘土室内平面应变试验的颗粒流PFC 模拟流程
试样是由不同半径的颗粒单元所组成, 颗粒半径R 的分布采用从R min 到R max 的均匀分布。经过大量颗粒流仿真方法PFC 试样的仿真试验, 对于文献[2]的试验情况仿真试样均取R min =5mm, R max /R min =217。
在PFC 程序中, 对于有内部凝聚力的材料是无法给定孔隙率的, 程序将根据R min 和R max /R min 自动设定。但由于PFC 所采用的颗粒均为球形颗粒, 且每个颗粒都与其周围其它颗粒相接触, 因此可以通过删除部分颗粒单元的办法来增大其孔隙比, 另外如果在PFC 试样中采用不均匀颗粒的组合方式, 则其孔隙比还将有更复杂的变化。
PFC 颗粒试样是通过颗粒间的相互作用来表达整个宏观试样的应力响应的, 对于粘性土颗粒试样PFC 所定义的颗粒间相互作用有:在接触处的法向接触力F ij 、切向接触力F ij , 接触摩擦力F ij , 下标表示力由第i 个颗粒单元通过接触作用于第j 个颗粒单元上; 根据颗粒流理论的接触本构假设, 上述各个相互作用力分量可分别通过法向刚度k ij 、切向刚度k ij 和摩擦系数L ij 和法向相对位移U ij 和切向相对位移U ij 按下式计算:
F ij =k ij U ij , F ij =k ij U ij , F ij =L ij F ij =L ij k ij U ij n s s s s s f n n n s [4]n s n 2D n s f 2D [4]2D 2D (1)
2D
n 法向接触力沿两颗粒单元圆心的连线, 切向接触力与摩擦力则与之垂直。对于粘性土颗粒试样PFC 给定颗粒之间的接触模量E c 以及定义法向刚度、切向刚度之比k n /k s , 在文中的计算中k ij =k , L ij =
L 。另外对于粘性土还需在PFC 试样的颗粒间定义法向、切向连接强度R c 和S c 以及各自的标准偏差。颗粒的比重可根据实际土样的干容重取为1189kN P m 。
表1 PFC 2D 数值试样的基本输入参数
围压R 3
/kPa
25
100
150
200试样尺寸/mm(高@宽) 70@2570@2570@2570@25粒径R min /mm[1**********]5
2D 3R min R max [1**********]7摩擦系数[**************]0颗粒法向接触刚度E c /Pa 4135@1064135@1064135@1064135@106颗粒刚度比[**************]0接触强度接触强墙体速率/Pa 113@104113@104113@104113@104度偏差[**************]5/mm [1**********]10收敛控制[**************]0初始锁定应力/Pa 75@10375@10375@103100@103从表1可以看出, PFC 颗粒流数值仿真试样的基本性质数据是相同的, 试样在不同围压下的仿真数
据差别主要在于颗粒试样的摩擦系数、试样加载过程中最终的墙体移动速率、试样初始的锁定应力状
|R y -R x |态, 以及加载试样最终收敛时的控制应力比, [A 。(R d ) max
11112 PFC 试样的颗粒生成过程 首先定义墙体, 共4道, 其包围的矩形为70m @25mm, 见图2(a) 。在PFC 生成墙体的过程中需要控制:(1) 上下墙体与试样宽度之比以及两侧墙体与试样高度之比要适当; (2) 颗粒与上下墙体以及两侧墙体间的接触刚度同颗粒自身间的法向接触刚度之比要适当。否则, 可能会出现颗粒在墙体运动的加载过程中由于没有墙体的约束而溢出、作自由运动的现象, 或穿透墙体而向四周散逸, 导致平衡迭代发散、实验失败。接着PFC 在上述给定的墙体范围内生成颗粒, 把半径在指定区间[R min , R max ]的颗粒往区域内填充, 如果没有已生成颗粒与之重叠, 则生成此颗粒, 否则, 改变颗粒的位置重试, 最后通过循环来消除试样内部非均匀应力。与无粘性的砂土颗粒试样不同, 粘性土颗粒试样的孔隙比是不能人为给定的, 而是由程序自动设定一个值如n =16%, 由此可以近似计算出试样的颗粒数目
N =[4]2D 2D 2D :(2) R min +R max , R =2式中:b 为试样的宽度; h 为试样的高度; N 为颗粒单元总数。
图2给出了用颗粒流程序模拟室内双轴试验整个试样的颗粒分布情况。图2(b) 为生成小粒径颗粒后颗粒的分布情况, 通过半径的膨胀, 颗粒分布情况见图2(c) 。半径膨胀后, 颗粒间将出现非平衡力, 通过循环让试样颗粒自动达到平衡状态(为加速达到平衡状态, 可以考虑增大颗粒间的摩擦力) , 此时试样内部颗粒分布情况见图2(d) 。
此外, 对于粘性颗粒试样, 为了消除试样内部可能存在的初始应力集中的不合理现象, 需要在试
2D 样加载之前先对试样施加一个初始的锁定应力(PFC 的模型参数) , 并在该应力作用下使试样达到初
始的均匀应力状态, R , 如x =R y =R 0, 一般R 0的取值大小相对于试样材料的峰值强度应该足够小小于单轴强度的1%。
经过上述工作基本上建立起一个能够和实际粘土试样比拟的PFC 数值仿真试样。
112 颗粒流数值仿真试验的实现
2D 11211 载荷与围压的施加和保持 PFC 是通过一套数值伺服系统让顶部和底部墙体作相对运动来施
加试样荷载的, 并同时调整两侧墙体的位移, 以保持试样的围压R x 恒定, 加载后的试样颗粒分布情况如图2(e)
所示。2D [4]
(a)墙体 (b) 墙体+颗粒(未膨胀) (c) 颗粒膨胀后的试样 (d) 平衡后颗粒试样
2D (e) 加载后的颗粒试样图2 室内双轴试验的颗粒流模拟11212 应力和应变的计算 在整个加载、卸载和再加载的过程中, PFC 程序可以自动记录试样的水
平E x 和竖向应变E y , 以及R x 和R y 。R x 和R y 由墙体和颗粒的接触力的总和除以颗粒和墙体接触的面积, 即颗粒试样的侧面面积和顶底部面积来计算:
R x =f x
2 r h , R y =f y
2 r b
(3)
式中:f x 为左侧或右侧墙体与颗粒单元的x 方向接触力, f y 为顶部或底部墙体与颗粒单元的y 方向接触力, r 为颗粒球体的平均半径, 在PFC 中则 r 取单元厚度。应变按下式计算:
E x =; E y =b h (4) 2D
式中:$b 和$h 分别为试样的x 和y 方向变形, 可以由墙体的位置来计算。
以(R y -R x ) 和E y 的记录作图便可得到试样的应力-应变关系曲线, 此即(R 1-R 3) ~E a 关系曲线。
2 粘土本构行为的颗粒流仿真研究
由于颗粒流试样实际上都是基于散粒介质本身微观特性基础上建立的, 故PFC 数值模型内颗粒间所发生的相互作用反映了散粒介质骨架的细观力学行为, 即基于颗粒流试样所模拟出的粘土介质宏观应力应变响应特性其实就是土体骨架所具有的本构行为。在文中的加载方式下, 可以不考虑瞬态流-固耦合问题, 颗粒流仿真试样得出的应力-应变曲线就是有效应力与应变的关系曲线。但在更深入的颗粒流试样地震荷载的响应特性分析, 则进一步的多相颗粒流理论研究是必需的。
211 轴向应力-应变关系曲线的PFC 仿真结果 以文献[2, 3]上海粘土固结不排水平面应变实验的部分数据结果为仿真对象, 给出4种围压下的应力-应变(R 1-R 3) ~E a 曲线对比。从图3围压为25kPa 的应力-应变曲线的比较来看, PFC 的数值仿真结果与实际试验实测曲线吻合得很好。但从围压分别为100kPa 、150kPa 、250kPa 的应力-应变曲线的比较还可看出, 细观颗粒各项指标相同的试样却无法得到较为理想的仿真结果, 见图4~图6。这说明, 颗粒流试样尚不能很好地反应物理试样客观上的唯一性, 不同的颗粒试样得到相似仿真曲线的数值试验从问题的另一个角度更进一步表明了该
点。2D 2D
图3 围压为25kPa 时的应力-
应变曲线对比图4 围压为100kPa 时的应力-
应变曲线对比
图5 围压为150kPa 时的应力-
应变曲线对比图6 围压为250kPa 时的应力-应变曲线对比
从见图3~图6的粘土试样实验曲线还能看出, 试样在曲线的起始阶段表现出一定的初始强度(初始切线模量较大) 。这实际上是由于粘土试样内部的细观结构特性而导致的。而对于颗粒流仿真试
样来说, 由于仅采用了单一的圆形颗粒单元来模拟整个结构, 因此无法揭示真实粘土试样内部的细观结构特性。显然, 该方面进一步的研究是必要的, 如果采用多个圆形单元组合成一系列不同形状的组合颗粒块来模拟, 有望得到更好的结果。
212 颗粒流模型若干参数对粘土本构行为影响的仿真研究 所谓的颗粒流数值仿真试样其实就是通过不断改变颗粒单元及其集合体组构方面的细观力学性质, 从而使之逼近真实材料的宏观力学响应。因此, 颗粒流理论对岩土试样的数值仿真过程就是一个不断修改、调试各项细观力学参数(包括颗粒粒径) 的过程, 现就此做若干参数研究如下(假定围压为25kPa) 。
21211 摩擦系数对试样本构行为的影响 图7给出颗粒间摩擦系数在011~019, 9组颗粒流试样应力-应变本构行为的仿真曲线。从不同摩擦系数应力-应变关系曲线的对比可以发现, 在其它因素不变的情况下随摩擦系数的增加, 应力-应变关系曲线的峰值会提高, 并且初始弹性模量也有一定的提高。从图中还可以发现, 最大轴向允许应变随着摩擦系数的增加而提高。另外, 摩擦系数越大, 软化现象越明显, 且有小的摩擦系数对应理想弹塑性关系, 而提高摩擦系数导致软化现象。
21212 加载速率(墙体移动速度) 对试样应力-应变曲线的影响 从图8中可以看出, 颗粒流试验仿真过程中墙体的移动速率(加载速率) 对仿真试样的应力-应变曲线影响不大, 尤其是对峰值强度以前的曲线部分几乎没有任何影响。而在峰值强度以后对渐进破坏部分的曲线以及允许应变有些程度不大的影响, 并且对于一组确定的颗粒流试样, 似乎存在一个最优的墙体移动速率(加载速率
)
图7 摩擦系数对仿真试样应力-应变曲线的影响
图8 墙体移动速率对仿真试样
应力-应变仿真曲线的影响
21213 颗粒接触刚度对试样应力-应变曲线的影响 从图9可以看出, 颗粒间接触刚度E c 对仿真试样的应力-应变曲线有较大的影响。显然, 随着E c 的增加, 试样的初始弹性模量增大, 其出现峰值强度对应的应变将降低; E c 对试样的峰值强度有影响, 但并不大。另外颗粒间接触刚度的大小对仿真试样应力-应变曲线的峰后形态以及允许应变也有一定的影响, 但并不象颗粒间摩擦系数那样具有
本质性的影响。
图9 颗粒接触刚度对仿真试样应力-应变仿真曲线的影响
21214 最大与最小颗粒粒径之比对试样应力-应变曲线的影响 仿真试样的最大与最小颗粒粒径比(与试样的孔隙比程函数关系) , 从图10可以看出, 对仿真试样应力-应变曲线的峰值强度及其对应的应变, 以及允许应变都有较大的影响。比值越大峰值强度越高、软化现象也越明显; 而对于允许应变来说, 最大与最小颗粒粒径比的影响是非线性的, 在其它因素确定的情况下同样存在一个最优值。
图10 最大与最小颗粒粒径之比对颗粒流应力-应变仿真曲线的影响
3 结论和前景展望
以文献[2]的工作为依据, 利用颗粒流计算机仿真程序对文献[2]的粘土平面应变试验结果进行了部分仿真, 并针对颗粒流试样的微观参数进行了多组试样的参数研究, 关于岩土试样宏观力学响应问题的颗粒流细观力学途径方面取得了以下研究结论, 如表2所示。
表2 岩土试样宏观力学响应PFC 2D 数值仿真的参数研究结果
初始弹性模量
摩擦系数f
墙体移动速率v
R m ax /R min
颗粒接触刚度E c 增大(小幅) 无影响无影响增大(大幅) 理想塑性y 软化趋势明显无影响有趋势影响不大最大轴向允许应变降低(大幅) 有非线性影响有非线性影响降低(小幅) 峰值强度增大(大幅) 无影响增大(大幅) 影响不大峰值强度对应的应变无影响无影响有非线性影响降低(大幅)
此外, 作者的研究还揭示出, 岩土试样本构响应问题的颗粒流仿真研究中存在的部分问题, 例如颗粒流细观力学模型的非唯一性, 以及粘土试样固结问题的颗粒流仿真途径还尚待解决, 这包括流固耦合、孔压和排水边界条件等若干细节的颗粒流仿真。
最后值得指出的是, 用颗粒流理论与方法来模拟岩土类散粒介质宏观力学响应的数值途径, 正如文中所介绍的, 是基于细观层面上所有散粒之间的机械类相互作用, 通过接触力、摩擦系数和接触强度来表现颗粒间的接触应力、内摩擦特性和凝聚力。之所以说是细观而非微观(如颗粒直径达到01001mm 的尺度) 主要是因为:(1) 模拟的颗粒直径过小时, 计算量将大幅增加, 仿真受到软硬件条件的限制而不可行; (2) 目前的颗粒流理论和方法还不能充分考虑粘粒效应, 如水化膜、范德华力等物化作用。另外, 由于尚未完全研究清楚颗粒流数值仿真试样与粘土物理试样之间的颗分对应关系, 所以在模拟过程中还具有一定多耙试探性(trial and validation) , 比如本文就取粒径015~1135mm 。可以明确的是, 不同粒径范围的散粒集合体的宏观力学响应是不同的, 这正如粗粒土与细粒土的力学性质差异。综上所述, 作者的工作还只能是一个初步的可行性研究, 该研究还有待深入。
参 考 文 献:
[1] 周健, 池毓蔚, 池永, 徐建平. 砂土双轴试验的颗粒流模拟[J].岩土工程学报, 2001, (6) :701-704.
[2] 董建国, 李蓓, 袁聚云, 赵锡宏. 上海浅层褐黄色粉质粘土剪切带形成的试验研究[J]. 岩土工程学报,
2001, (1) :23-27.
[3] 李蓓. 平面应变仪的研制及上海地区粘性土剪切带形成的试验研究[D].上海:同济大学, 1999.
[4] Itasca Consulting Group, Inc.. PFC 2D User . s M anual (Version 210) [Z]. Minnesota, USA, 1999.
[5] Cundall P A, Strack O D L. A Discrete Numerical Model for Graunlar Assemblies. Geotechnique, 1979, 29(1) :47-65.
[6] Huang H, Detournay E, Bellier B. Discrete element modelling of rock cutting [J]. Rock Mechanics for Industry, 1999,
1(1) :123-130.
[7] 冯国栋. 土力学[M]. 北京:水利电力出版社, 1986.
[8] 周健, 池永, 池毓蔚, 徐建平. 颗粒流方法及PFC 2D 程序[J].岩土力学, 2000, (3) :271-274.
Simulation of plane strain test of clay by means of particle flow code
LI AO Xiong -hua , ZHOU Jian , XU Jian -ping , LIN L-i min
(11Tongji Universit y , Shanghai 200092, China ;
21Gene ral O ffice o f Rail way Traffic Company Ltd , Wuhan City , W uhan 430017, China ;
31Zhe jian finance and e conomics college , Hangzhou 310012, China ) 1123
Abstract :The simulation of plane strain test for clay using particle flow code method is carried out. The comparison of simulation with test result shows that the model based on particle flow theory can reflect the constitutive behaviors of soil by means of changing the particle . s properties and size grading characteris -tics during modeling. The model can be used to predict the constitutive behaviors of clay on the level of theoretical. The relationship between microscopic para meter of 2-D particle flow and mechanical response of physical sample are also investigated. Some basic laws are found.
Key words :particle flow theory; clay; plane strain test; constitutive behavior; stress -strain curve (上接第10页)
参 考 文 献:
[1] 李锦秀, 等. 三峡库区江段纵向离散系数研究[J]. 水利学报, 2000, (8) :84-86.
[2] EPA/600/3-87/007. The Enhanced Stream Water Quali ty Models QUAL2E and QUAL2E -UNCAS [M]. Docu menta -
ti on and user Manual, 1987.
[3] 余常昭, M 马尔柯夫斯基, 李玉梁著. 水环境中污染物扩散输移原理与水质模型[M ].中国环境科学出
版社, 1989.
[4] R M Wrigh t, A J Mcdonell. In -stream Deoxygenation rate prediction [J]. Proc. ASCE. J. Env. Div. 105, 1979, 323
-333.
[5] 傅国伟编著. 河流水质数学模型及其模拟计算[M]. 北京:中国环境科学出版社, 1987.
Numerical simu lation of water quality for the Three Gorges Project Reservoir
LI Jin -xiu , LI AO Wen -geng , HUANG Zhen -li 112
(11China Institute o f W ate r Resourc es and Hydropowe r Rese arch , Bei jing 100083, China ;
21Executive O ffic e o f State Council for Thre e Gorges Projec t Construction Committe e , Be ijing 100038, China )
Abstract :The 1-D unsteady model based on finite difference method with four -point implicit scheme simulating the water flow and water quality in the reservoir of Three Gorges Project is established. The major indices of water quality including 14items are investigated. The effects of significant change of flow condition in the reservoir due to impounding are taken into account. The coefficients of the model includ -ing the longitudinal dispersion coefficient, decay rate of B OD 5and atmospheric re -aeration rate of DO etc. are established empirically. The model is applied to predic t the variation trend of water quality in the res -ervoir.
Key words :Three Gorges Project; water quality model; coefficients of water quality; empirical formula