三维离散颗粒单元模拟无黏性土的工程力学性质.pdf

第30卷 第2期 岩 土 工 程 学 报 Vol.30 No.2 2008 年 2 月 Chinese Journal of Geotechnical Engineering Feb., 2008 三维离散颗粒单元模拟无黏性土的工程力学性质 罗 勇，龚晓南，连 峰 浙江大学岩土工程研究所，浙江 杭州 310027 摘 要基于三维离散单元颗粒流理论，引入了在极限剪力状态下颗粒间可以发生滑动的接触本构模型，建立颗粒体 试样，并赋予颗粒相应的细观结构参数进行无黏结摩擦材料的三维应力应变数值模拟。通过大量的颗粒流数值试验， 从细观力学角度对砂土的工程力学特性进行了初步的模拟研究。对砂土的室内常规三轴试验及其剪切带形成和发展进 行了数值模拟，分别对比了不同围压下三维颗粒体试样与室内三轴的应力应变关系曲线，基本再现了试样的加载曲线。 通过进一步分析颗粒细观结构参数变化对组成材料的宏观力学的影响及其剪切位移场的形成和发展规律，结果表明颗 粒介质之间的摩擦系数、颗粒组成材料的孔隙率对材料的宏观力学行为有比较明显的影响。研究的意义在于理想的颗 粒流数值模拟试验能够突破常规的室内土工试验能力及其局限性，对于砂土的细观结构的影响研究成为可能，并且揭 示了微观的一些规律，对今后的岩土工程的颗粒模拟提供了有价值的参考。 关键词颗粒流理论；滑动模型；砂土；细观参数 中图分类号TU411 文献标识码A 文章编号1000–4548200802–0292–06 作者简介罗 勇1977– ，男，贵州毕节人，博士，主要从事软淤处理、基坑工程、边坡工程等岩土工程数值分析 和岩土计算力学研究。E-mail luoyong7966。 Simulation of mechanical behaviors of granular materials by three-dimensional discrete element based on particle flow code LUO Yong，GONG Xiao-nan，LIAN Feng Institute of Geotechnical Engineering, Zhejiang University, Hangzhou 310027, China Abstract The triaxial tests of sand by PFC3D particle flow code in three dimensions were simulated. The slip model was used to constitute samples of sand. The slip model provided no normal strength in tension and allowed slip to occur by limiting the shear force. Preliminary study on mechanical properties of soil was done by PFC3D numerical simulation through endowing appropriate micro-properties. Good agreement between the results of numerical simulations and laboratory tests was shown. Macro properties of sand samples with various parameters of particle elements were presented, and the results were valuable for developing soil constitutive models. The emergence and development of shear bands were simulated numerically by PFC3D under different confining pressures. It was shown that the microscopic parameters of granular materials affected markedly the macroscopic response of the noncohesive frictional material. The meaning of this research lies in the fact that PFC mode could break through the limits of real lab tests and make the microcosmic research of soil possible. In addition, some laws were discovered, and it would be helpful to the study of the application of particle to the geotechnical problems in future. Key words particle flow code; slip model; sand; microscopic parameter 0 引 言 由于土体由破碎的固体颗粒组成，土的宏观变形 主要不是由于颗粒本身变形，而是由于颗粒间位置的 变化。这样在不同的应力水平下由相同应力增量而引 起的应变增量就不会相同，亦即表现出非线性。土体 的结构可视为一个由单粒、集粒或凝块等骨架单元共 同形成的空间结构体系，它的单元形态确定力的传递 性能和土的变形性质，它的连接方式确定土的结构强 度，它的排列方式确定土的稳定性。土体这种细观结 构的复杂性决定了土体具有复杂的工程力学特性,所 以，对于土体的这种复杂工程特性的研究可以从细观 形态学途径入手。但是这种研究非常复杂，很难建立 宏观与微观的定量关系。而理想离散单元颗粒流方法 克服了传统连续介质力学模型的宏观连续性假设，可 ─────── 收稿日期2006–12–22 第 2 期 罗 勇，等. 三维离散颗粒单元模拟无黏性土的工程力学性质 293 以从细观角度对土的工程特性进行数值模拟，并通过 颗粒结构的细观参数的研究来分析材料的宏观力学行 为。 国外学者Bardet和Proubet[1-2]应用理想的二维颗 粒集合体模拟了粒状材料中剪切带的结构,对剪切带 的厚度、带内位移、孔隙比、体应变以及颗粒旋转等 进行了研究。 周健等[3-4]应用颗粒流理论模拟了砂土和 黏性土平面应变试验的应力–应变关系曲线。张洪武 等[5]建立了非线性的接触本构模性模拟了二维平面应 变的颗粒圆筒的数值模拟。对于三维的情况的颗粒应 力–应变数值模拟的和细观参数对其影响的工作还不 多见。 本文的研究重点是基于三维颗粒流理论和开发颗 粒流数值模拟技术，对无黏性土体应力应变关系和剪 切带形成机理进行细观数值模拟，并从颗粒体的复杂 的细观结构入手，通过不同的细观参数分析宏观的力 学行为， 将土体微细观结构与宏观力学反应联系起来， 对土体工程力学特性和剪切带形成与发展等渐进破坏 过程有更深入的了解和发现。 1 离散颗粒流的基本原理 颗粒离散元法是离散元法中最理想的一种，它把 离散体看作有限个离散颗粒单元的集合体，每个颗粒 为一个单元，将材料理想化为相互独立、相互接触和 相互作用的颗粒群体。颗粒流理论在整个计算循环过 程中，交替应用力–位移定律和牛顿运动定律。通过 力–位移定律更新接触部分的接触力。 通过运动定律， 更新颗粒–颗粒与颗粒–边界的位置， 达到新的平衡。 1.1 作用力与位移的关系 颗粒流理论中，颗粒与颗粒之间的接触会产生力 的作用。将此力分为法向力 n i F和切向力 s i F ns nii FFF ， 1 法向作用力与位移的代数关系可写为 nnn ii FK U n ， 2 式中， n K为法向刚度，属于割线模量，与总位移和 力对应， n U 为法向位移量（颗粒–颗粒或颗粒–墙 体的变形重叠量）， i n为接触面上的单位法向向量， 如图1所示。 切向作用力与位移的代数关系可写为 sss ii FkU∆ −∆ ， 3 式中， s k为切向刚度，属于切线模量，与位移和力的 增量对应， s i U∆为切向位移增量，可写为 ss ii UVt∆∆ ， 4 式中， s i V为切向接触速度，t∆为的计算时步。切 向作用力写为 sss iii FFF← ∆ 。 5 通过方程（2）和（5），在每一计算时步结束， 有 [][] [][] [][][][] 3 333 [][][][] 3 333 AA i ii BB i ii AACA k jkii BBCB k jkii FFF FFF MMexxFM MMexxFM ⎫ ←− ⎪ ←⎪ ⎪ ⎬ ←−−− ⎪ ⎪ ←− ⎪⎭ ， ， ， ， 6 式中， []A i F， []B i F， [] 3 A M和 [] 3 B M分别为颗粒A和颗粒 B的力和弯矩矢量。 图 1 颗粒接触模型 Fig. 1 Contact model of ball-ball 1.2 运动法则 每一个颗粒的运动是由作用于其上的剩余力和剩 余力矩决定,可以用颗粒内一点的线速度与颗粒的角 速度来描述。 运动方程由两组向量方程表示,一组是剩 余力与平动的关系，另一组是表示剩余力矩与旋转运 动的关系 iii Fm xg 平动方程 ， 7 ii MH 旋转运动方程 。8 式中， i F为剩余力，m为颗粒质量， i g为体积力加速 度（如重力加速度）， i M为剩余力矩， i H为角动量。 2 颗粒接触本构模型 本文的研究重点是无黏结材料的力学性能，属于 摩擦材料，故引入接触摩擦滑动本构模型来模拟颗粒 之间的滑动行为。这种模型的特点是，滑动模型在相 互接触颗粒之间没有法向和切向抗拉强度，允许颗粒 在其抗剪强度范围内发生滑动，该模型适用于模拟颗 粒间不存在黏结力的散离体材料（如砂土、堆石等）。 其本构行为可以描述为 sn maxi FF ， 9 式中，为颗粒的摩擦系数，如果当 ss maxi FF，则颗 粒发生滑动，在下一时步计算开始时，通过下面公式 294 岩 土 工 程 学 报 2008 年 表 1 三维颗粒试验样本生成控制参数 Table 1 Sand properties in PFC3D sample 参数 比重 初始孔 隙率 n 颗粒半径/mm 粒径放大率 摩擦系数 颗粒接触法向刚度 /Nm -1 颗粒接触刚 度比 砂 2.63 0.36 0.075～0.1 1.66 1.2 2.5310 4 1.0 设置 s i F等于 s max F， ssss max / iii FFFF← 。 10 3 无黏性土的工程力学的细观模拟 3.1 应力–应变关系曲线的数值模拟 应力–应变的数值试验的对比依据以李广信[6]的 室内三轴常规试验进行的承德中密砂为基础。这种天 然砂土的物理力学指标为平均粒径d500.18 mm，不 均匀系数Cu2.8，颗粒比重Gs2.63，最大孔隙比e 0.80，最小孔隙比e 0.40，试样相对密度Dr 0.63， 干密ρd1.70.01 g/cm3。先建立颗粒流（PFC3D）模 型模拟其应力–应变关系曲线，并以此PFC3D模型为 基本模型进行细观参数的调整，来分析颗粒细观参数 的变化对材料宏观特性的影响。 三维颗粒试验样本的生成第一步是要建立约束 散离体的结构空间，容量空间包括一个施加围压的完 全没有摩擦的柔性圆筒和两个施加轴向载荷的无摩擦 平面刚性加载压盘；第二步是在给定参数算法的前提 下生成颗粒体试样， 其中颗粒控制参数如表1所示。 三 维颗粒试验样本的生成图示如图2所示。 图 2 三轴数值试验示意图 Fig. 2 Sketch map of triaxial test of PFC3D 图 3 三轴试验宏观应力应变关系曲线比较 Fig. 3 Comparison between PFC3D model and lab test results 建立颗粒流PFC3D数值模型模拟砂土的应力–应 变关系曲线，并以此PFC3D模型为基本模型进行反复 的细观参数的调整，来分析颗粒细观参数的变化对材 料宏观力学特性的影响。为了节约计算时间和提高效 率，颗粒流理想砂土的数值模型的尺寸大小为高度 410 -3 m，直径210-3 m的圆柱体，由2686个颗粒实体 构成。颗粒流的砂土试样分别在围压为100，300，500 kPa的固结围压下， 加载主应力， 通过编制伺服式加载 程序严格控制围压的误差在加载过程中小于1， 使试 样在每一步的平衡状态下加载，直到加载结束。经过 大量的参数试算和分析，最终确定颗粒的参数如表1 所示的情况下，颗粒流的三维数值试验应力–应变关 系基本与室内承德中密砂的三轴试验的应力–应变关 系曲线变化规律一致。 应力–应变图示如图3所示， 摩 尔强度包线如图4所示。 图 4 内摩擦角度比较 Fig. 4 Comparison of the friction angle 3.2 影响应力–应变关系曲线诸多因素 为了清楚的展现细观参数的变化对砂土的力学性 能的影响，本文选择了下面的几个因素分别对三轴应 力应变曲线进行模拟和分析。 （1）孔隙率的影响 三轴试验样本选择孔隙率分别为0.35，0.40，0.45 和0.50四种情况下进行加载试验，应力–应变曲线变 化规律如图5所示， 从图中可以看出， 和室内三轴试验 规律一样，随着孔隙率的减少，峰值强度增加，曲线 初始模量也相应提高， 应变软化加剧。 从图6的围压大 小为0.5 MPa的情况， 体积应变–轴应变关系曲线看出 随着孔隙率的增加体胀明显减少。 孔隙率为0.5的情况 下，在数值试验过程中没有体积膨胀发生，呈现出松 砂在中围压下无剪胀的特性。围压为0.1 MPa的时候， 从图7也能得出上述同样的结论， 在较低的围压下， 出 现明显的两极分化，即是颗粒越是密实（孔隙率越小 的情况） ，体胀越是明显；颗粒越是松散（孔隙率大的 情况） ，体胀就越小，甚至只表现出体缩的特征，如孔 第 2 期 罗 勇，等. 三维离散颗粒单元模拟无黏性土的工程力学性质 295 隙率为0.5的松散颗粒在围压0.5 MPa和0.1 MPa的情 况下没有体胀发生，可以这样理解，密实的砂土颗粒 间要产生相对错动就不得不向上翻越，约束压力越小 就越容易膨胀。对于松散的砂土而言，要产生相对错 动就只需挤密周围的孔隙就可以，不用上下翻越，特 别在约束压力小的情况，体胀现象会完全消失。 图 5 孔隙率对应力–应变的影响0.5 MPa Fig. 5 Stress-strain curves with various porosities 0.5 MPa 图 6 孔隙率对体积应变的影响0.5 MPa Fig. 6 Variation of volumetric strain vs. axial strain with various .porosities 0.5 MPa 图 7 孔隙率对应力–应变的影响0.1 MPa Fig. 7 Variation of volumetric strain vs. axial strain with various .porosities 0.1 MPa 从孔隙率的变化对应力–应变影响的数值试验得 出，如图5和图8所示，无论颗粒是松散的还是密实的 情况下，随着约束应力的 3 σ的增加，剪应力 13 σσ− 与其呈一定比例增加， 这也是散粒体的一个共同特性， 这也是土区别其他材料的重要特性之一[6-10]。 图 8 孔隙率对应力–应变的影响0.1 MPa Fig. 8 Variation of volumetric strain vs. axial strain with various .porosities 0.1 MPa （2）颗粒间摩擦系数的影响 由于理想的离散单元颗粒流用到的模型是规整的 球形颗粒，与真实颗粒间的咬合特征不同，要想模拟 以摩擦特性为主的材料，比如砂土、碎石等，就必须 通过研究颗粒接触之间的摩擦系数的影响， 也就是说， 摩擦系数对应力应变的关系是十分明显的。 如图9和图 11所示， 接触摩擦系数的提高， 峰值强度明显的增加， 应变软化加剧，体积膨胀，如图11所示。随着围压的 增加，摩擦系数对峰值强度的影响就越明显，应变软 化现象就越加剧，体积膨胀减少，如图12所示。 图 9 摩擦系数对应力–应变的影响0.1 MPa Fig. 9 Stress-strain curves with various friction coefficients 0.1 MPa 图 10 摩擦系数对应力–应变的影响0.5 MPa Fig. 10 Stress-strain curves with various friction coefficients 0.5 MPa 296 岩 土 工 程 学 报 2008 年 图 11 摩擦系数对体积应变的影响0.1 MPa Fig. 11 Variation of volumetric strain vs. axial strain with various .friction coefficients 0.1 MPa 图 12 摩擦系数对体积应变的影响0.5 MPa Fig. 12 Variation of volumetric strain vs. axial strain with various .friction coefficients 0.5 MPa 3.3 剪切位移场的形成与发展规律 利用颗粒流的三维数值模拟应力–应变的曲线效 果比较理想，在试验过程中跟踪PFC3D模型中颗粒的 位移和其它变量的变化情况,进而对试样中剪切带的 形成及其性状进行分析。可以进一步利用此模型来分 析在施加轴向剪力的过程中，试样的剪切位移场发展 趋势和变化规律。 如图13所示为围压为0.1 MPa时候， 剪切位移场的 形成和变化规律，其中图13（a）和图13（b）表示在 峰值强度未达到之前的剪切位移场，试样处在初始体 胀阶段，从图中的颗粒剪切位移场看出处于比较杂乱 无章的状态。由于在剪胀过程中土颗粒一般是从低势 能状态向高势能状态变化，图13（c）和图13（d）表 示随着体胀的继续发展和轴向应变的增加的剪切位移 场， 说明此时试样内部出现迅速的颗粒结构排列变化， 这种颗粒排列规律的变化将引起最终剪切带的形成。 颗粒剪切位移场的定向趋势越来越呈现出规律性，剪 切带的形状趋于明显且集中，厚度逐渐变小，处于不 稳定状态。 分析不同围压下砂土试样位移场发现，如图14随 着围压的增大，剪切带的宽度逐渐变小，这与砂土试 验结果类似。从细观机理上分析，围压越大，对土样 内部颗粒的约束越大，其细观结构发生巨大的重新排 列的机率越小，更不容易产生颗粒之间的相对滑动破 坏，发生局部破坏的范围也越小。故在很多室内砂土 三轴试验证实，随着围压的提高，砂土的内摩擦角度 将逐渐减少，破坏包线呈现明显的抛物线[11-14]。 图 13 剪切位移场的形成和变化规律 Fig. 13 Emergence and development of shear bands 图 14 不同围压下剪切位移场 Fig. 14 Shear bands of various confining pressure 4 结论与建议 本文基于三维颗粒流离散单元理论的数值试验突 破常规上土工试验的局限性，对砂土的宏观应力应变 和剪切位移场形成和发展规律的模拟，得到了与室内 试验相似的结果，基本再现了砂土室内试验的应力– 应变关系曲线的特征。在此基础上研究了颗粒结构的 细观参数对应力–应变曲线的影响，得出了与室内试 验基本一致具有普适性的结论。 （1）随着孔隙率的减少，试样初始模量增加，应 变软化加剧，体胀增加。 （2） 随着接触摩擦系数的提高， 试样峰值强度增 加，初始弹性模量也有小幅的提高，软化有加剧的趋 势，体胀轻微变缓。 （3）由于在剪胀过程中土颗粒一般是从低势能 状态向高势能状态变化，随着体胀的继续发展和轴向 应变的增加， 试样内部出现迅速的颗粒结构排列变化, 这种颗粒排列规律的变化将引起最终剪切位移场的形 成。 （4） 对于散粒体的一个共同特性， 就是随着约束 围压的提高，峰值强度和初始变形模量都增加，一般 情况下伴随有明显的应变软化现象，这也是土区别其 他材料的重要特性之一。 而对于颗粒流试样来说，由于仅采用了单一的球 形颗粒单元来模拟整个结构，因此无法全面揭示真实 第 2 期 罗 勇，等. 三维离散颗粒单元模拟无黏性土的工程力学性质 297 砂土试样内部的细观结构特性。比如说，砂土颗粒的 不规则外观对其工程力学性能的影响等，为此，在今 后的进一步模拟过程中可以考虑开发利用PFC3D的 扩展功能中的簇单元块体，建立不规则的颗粒集合刚 体,有望得到更好的结果。 类似砂土的宏观工程力学研 究，可以引入颗粒之间的接触黏结本构模型和平行黏 结本构模型分别来模拟具有黏结特性的材料和水泥固 化后的复合材料，比如黏土和水泥土等，研究其宏观 工程性质和力学损伤特性，这将是下一步的工作。 参考文献 [1] BARDET J P, PROUBET J. A numeeerical investigation of stnlcture of persistent shear bands in granular media[J]. Geotechnique, 1991, 414 599–613. [2] BARDET J P, PROUBET J. Shear-band analysis in idealized granular material[J]. J Engrg Mech, 1992, 1188 397–415. [3] 周 建, 池毓蔚, 池 永, 徐建平. 砂土双轴试验的颗粒流 模拟[J]. 岩土工程学报, 2000, 226 701–704. ZHOU Jian, CHI Yu-wei, CHI Yong, XU Jian-ping. Simulation of biaxial test on sand by partical flow code[J]. Chinese Journal of Geotechnical Engineering, 2000, 226 701–704. in Chinese [4] 周 建, 廖雄华, 池 永, 徐建平. 土的室内平面应变试验 的颗粒流模拟[J]. 同济大学学报, 2002, 309 1044–1050. ZHOU Jian, LIAO Xiong-hua, CHI Yong, XU Jian-ping. Simulating plane strain test of soils by particle flow code[J]. 2002, 309 1044–1050. in Chinese [5] 张洪武, 秦建敏. 基于非线性接触本构的颗粒材料离散元 数值模拟[J]. 岩土工程学报, 2006, 2811 1964–1969. ZHANG Hong-wu, QIN Jian-min. Simulation of mechanical behaviors of granular materials by discrete element based on mesoscale nonlinear contact law[J]. Chinese Journal of Geotechnical Engineering, 2006, 2811 1964–1969. in Chinese [6] 李广信. 高等土力学[M]. 北京 清华大学出版社, 2004 32 –113. LI Guang-xin. Advanced soil mechanics[M]. Beijing Tsinghua University Press, 2004 32–113. in Chinese [7] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 291 47– 65. [8] PULS M, PULSFORT M, WALZ B. Application of the PFC3D for determination of soil properties and the simulation of the excavation process in front of sheet pile wall constructions [C]//Numerical Modeling in Micromechnics via Particle s, London, 2004 35–44. [9] PIERCE M E, CUNDALL P A, VAN HOUT G J, LORIG L. PFC3D modeling of caved rock under draw[J]. Numerical Modeling in Micromechanics via Particle s, 2003 211–217. [10] CUNDALL P A. Numerical modeling of jointed and faulted rock[J]. Mechanics of Jointed and Faulted Rock, 1990 11– 18. [11] POTYOND D O, CUNDALL P A. A bonded-particle model for rock[J]. Int J Rock Mech Mtn Sci, 2004, 418 1329-– 1364. [12] CUNDALL P A. Numerical modeling of jointed and faulted rock[J]. Mechanics of Jointed and Faulted Rock, 1990 11– 18. [13] CATHERINE O’Sullivan, et al. Examination of the response of regularly packed specimens of spherical particles using physical tests and discrete element simulations[J]. J Eng Mech, 2004, 1404 1140–1150. [14] PFC3D-Particle flow code in 3 dimensions[M]. Minneapolis Itasca, Inc, 2003.