基于随机散粒体模型的堆石体真三轴数值试验研究.pdf

第 34 卷 第 4 期 岩 土 工 程 学 报 Vol.34 No.4 2012 年 .4 月 Chinese Journal of Geotechnical Engineering Apr. 2012 基于随机散粒体模型的堆石体真三轴数值试验研究 周 伟 1,2，刘 东1,3，马 刚1,2*，常晓林1,2 1. 武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉 430072；2. 水工岩石力学教育部重点实验室，湖北 武汉 430072； 3. 中国水电顾问集团中南勘测设计研究院，湖南 长沙 410014 摘 要基于三维变形体离散单元法，模拟堆石体的真三轴试验，研究堆石体在三向不等应力状态下的强度和变形特 性。真三轴数值试验装置采用六刚性板加荷方式，等中主应力比路径加载。数值试验结果表明真三轴数值试验能较 好反映堆石体在三向不等应力状态下的应力和变形规律，数值模拟得到宏观应力变形特性与试验规律相似；等中主应 力比参数 2313 bσσσσ−−的大小对堆石在三向应力状态下应力和强度特性有显著的影响， 堆石体在 3 个加载方 向的变形也随该参数的变化而变化， 应力比参数b从 0 到 1 变化过程中中主应变方向先压缩后膨胀， 小主应变方向一直 处于压缩状态； 堆石体的内摩擦角随着b值的增加而增大， 基本符合 Lade-Duncan 破坏准则； 在细观层面上， 围压越高， b值越大，颗粒配位数越大；加载过程中，颗粒接触法向和法向接触力各向异性程度加强，各向异性主方向角由水平向 转向大主应变方向，试样各向异性系数的演化规律和试样的宏观应力变形曲线相对应，试样的宏观力学特性与细观组 构存在内在的关联。 关键词堆石体；真三轴；散粒体模型；细观力学；数值试验 中图分类号TU43 文献标识码A 文章编号1000–4548201204–0748–08 作者简介周 伟1975– ，男，博士，教授，博士生导师，主要从事高坝结构数值仿真方面的教学与研究工作。 E-mail zw_mxx。 Numerical simulation of true triaxial tests on mechanical behaviors of rockfill based on stochastic granule model ZHOU Wei1, 2, LIU Dong1, 3, MA Gang1, 2, CHANG Xiao-lin1, 2 1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China; 2. Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering Wuhan University, Ministry of Education, Wuhan 430072, China; 3. Hydrochina Zhongnan Engineering Corporation, Changsha 410014, China Abstract The three-dimensional deable discrete element is employed to study the strength and deation characteristics of rockfill under three-dimensional unequal stress states by means of true triaxial numerical tests. In addition, the six-rigid-plate loading and the constant value of the intermediate principal stress ratio loading paths are adopted. The simulated results show that the true triaxial numerical tests can preferably reflect the stress and deation of rockfill under the three-dimensional unequal stress states, and the numerical results are the same as the test ones. The stress ratio parameters significantly affect the stress intensity of rockfill under the three-dimensional stress state and the three-dimensional strain also changes with the increase of the ratio parameters. To be more specifical, when the intermediate principal stress ratio parameter b increases from 0 to 1, the second principal strain direction changes from compression to expansion while the third principal strain direction has always been in the state of compression, and the the internal friction angle of rockfill increases with the increase of b, which basically cons to the Lade-Duncan failure criterion. In the mesoscopic level, the higher the confining pressure is, the larger the value of b and coordination number of particles are. In the loading process, the anisotropy degree of contact normal orientation and contact normal force are strengthened, the main direction of contact force turns to the loading direction, the anisotropic coefficient evolution curve is similar to the stength curve, in other words, the strength of the rockfill samlpes has some relationship with the normal contact. Key words rockfill; true triaxial test; granule model; mesomechanics; numerical simulation 0 引 言 堆石体由于取材方便，对地形、地质条件有良好 ─────── 基金项目国家自然科学基金项目（50979082） 收稿日期2011–03–31 *通讯作者 第 4 期 周 伟，等. 基于随机散粒体模型的堆石体真三轴数值试验研究 749 适应性，在堆石坝、海岸工程和路堤建设等得到广泛 应用， 这就迫切需要我们全面了解堆石体的力学特性。 堆石体一种复杂的不连续颗粒材料，它在外荷载作用 下的响应是颗粒运动和颗粒间接触作用的综合体现。 目前，对于堆石体的力学性能的研究主要是以常规三 轴试验为主，试样处于轴对称的应力状态，然而在实 际工程中，堆石体一般处于三向不等应力状态下，因 而常规三轴试验结果不能准确、全面地反映堆石力学 性能。 由于真三轴试验[1]操作比常规三轴试验复杂，设备 的设计、制作和加载都比较困难，在岩土工程领域，国 内外学者对于颗粒材料的真三轴试验研究还相对较少。 李广信[2]对中密砂进行了真三轴试验，研究了中主应力 对砂土强度曲线的影响；朱俊高等[3]、徐志伟等[4]对土 体和粉砂侧向变形性状进行了真三轴试验， 表明土体侧 向变形可能收缩，可能膨胀；Shi等[5]、施维成等[6]对粗 粒料进行了真三轴试验研究， 结果表明中主应力对其强 度和变形有较大影响； Ko等[7]采用中主应力从小主应力 开始增大的加载方式进行真三轴试验， 探索了砂试样的 变形和破坏过程， 以及内摩擦角在加载过程中先变大后 有所减小的变化规律。Lade等[8-10]也在无黏性材料的应 力应变特性上做了真三轴试验研究， 指出中主应力比的 增大会使中主应变从膨胀变为压缩， 并在真三轴试验结 果基础上提出了Lade-Duncan破坏准则；日本的NAKAI 等[11-12]对砂石的试验研究中也得到了相似的结论。 近年来，计算机科学和数值计算方法的发展使我 们能够从细观角度来研究堆石体的力学性能，特别是 离散元法（DEM）的迅展发展，使我们可以更方便的 从细观结构方面对堆石颗粒材料进行研究。与室内试 验相比，堆石体离散元数值模能够克服试验设备，时 间和成本的制约，同时又能和室内试验相互印证、补 充，所以，有必要进行堆石体真三轴数值试验。在真 三轴试验数值拟方面，殷建华等[13]采用混合边界对土 体的应力应变特性进行真三轴数值试验研究，得出混 合边界能够很好的解决真三轴试验当中的边角效应问 题； Callisto等[14]对自然和重塑状态下的比萨黏土分别 进行了常规三轴和真三轴的数值模拟试验，研究了土 体的各向异性和本构关系；Thornton[15]对散粒体材料 进行真三轴数值仿真，研究结果表明，离散元法能够 较真实的模拟散粒体材料的应力应变扩张现象，颗粒 间的运动和接触行为影响着颗粒材料的力学性能； Belheine[16]应用三维离散元法进行土体排水剪数值仿 真试验，从细观参数层面分析了土体颗粒的强度变形 特性；Tang-Tat Ng[17-18]对颗粒材料在不同应力路径下 的力学性能做了真三轴数值模拟研究，验证了颗粒材 料的各向异性，同时得出不同的应力路径对于颗粒材 料力学性能有一定的影响。 Thornton.c和Tang-Tat N采用颗粒离散元进行真三 轴数值试验研究，然而对于堆石体这种多面体颗粒不 容易真实模拟其形态。本文基于变形体离散元法，在 随机数值模拟技术基础上生成三维多面体颗粒，颗粒 的物理力学性能符合Weibull分布函数， 采用周伟等[19-21] 建立的随机颗粒不连续变形SGDD模型来模拟堆石 体的真三轴试验。本文通过真三轴数值试验研究了堆 石体的变形及不同应力路径下的应力应变关系，探讨 了堆石体细观层面的力学性能。 1 堆石体随机散粒体模型 堆石体在本质上是一种不连续介质，由块石、碎 石等多种成分混合而成，属于无黏性颗粒材料，堆石 体的外观形态主要为有棱角的凸多面体，颗粒之间能 相互咬合，并且堆颗粒之间的接触一般为点接触。 本文基于变形体离散元法，采用随机模拟技术， 结合堆石颗粒的特征参数，包括级配、尺寸、孔隙率， 生成堆石体颗粒及其在空间中的分布，由于堆石颗粒 是通过在空间中随机布点形成多面体形状，能够比较 真实地反映堆石的实际形态。堆石体颗粒之间的相互 作用是通过虚拟的法向和切向弹簧来传递，颗粒间的 接触为线性刚度接触模型，颗粒的物理力学参数符合 Weibull概率分布。 2 堆石体真三轴试验细观数值模拟 2.1 试样制备 根据图 1 所示的试样级配曲线，采用随机模拟技 术生成多面体颗粒，再压缩成指定大小的试样，试样 尺寸为 120 mm60 mm120 mm， 最大粒径为 max d 10 mm。在建模的时候，考虑到计算能力的限制，去 掉粒径小于 1 mm的小颗粒，采用较大粒径颗粒等体 积替代，堆石体随机颗粒模型如图 2 所示。由于模型 颗粒形状极不规则，采用二阶四面体单元离散，堆石 体试样的相关参数见表 1。 图 1 堆石料级配曲线 Fig. 1 Particle size distribution curves of rockfill 750 岩 土 工 程 学 报 2012 年 图 2 随机颗粒模型 Fig. 2 Stochastic granule model 表 1 试样颗粒数值模型相关参数 Table 1 Some parameters of granular samples 孔隙率 颗粒数 单元数 节点数 35 5796 90619 257609 2.2 数值模拟加载方式 真三轴试验加载方式一般有 3 种刚性板加载方 式、柔性带加载方式和复合加载方式。本文采用刚性 板加载方式，加载模具如图 3 所示。为了模拟室内真 三轴试验的试验条件，模型上下两端及两侧为 6 块刚 性板，底部刚性板全约束，顶板铅直向自由约束，侧 板分别在加压方向自由约束。 图 3 加载装置 Fig. 3 Loading device 本文采用等应力比参数b加载路径， 它反应了中主 应力与大小主应力相对关系，其表达式如下 23 13 b σσ σσ − − ， 1 由（1）式变形得 2313 bσσσσ−− ， 2 加载过程中保持小主应力 3 σ不变， 同时增加 1 σ和 2 σ， 保持固定的应力比参数，式（2）两端对时间t求导， 可得 21 ddddtbtσσ ， 3 即可知， 等b加载过程中， 第二主应力的加载速率是第 一主应力加载速率的b倍。 数值模拟开始时，先对试样三向施加围圧 3 σ进行 固结，小主应力方向（Z方向）保持围压不变，顶板 进行轴向 （Y方向） 应变加载， 加载速率为0.0005 mm/ 步，按照式（1）保持b值一定，中主应力（X方向） 和小主应力采用应力加载。 2.3 细观参数取值 堆石体细观数值模拟中涉及到的物理参数见表 2，其中 u ϕ为颗粒间的摩擦角； n k， s k分别为法向和 切向接触刚度； E u， E σ分别为细观单元弹性模量概 率分布的均值和标准差； f u， f σ分别为细观单元抗压 强度概率分布的均值和标准差；rho是细观单元弹性 模量与抗压强度概率分布之间的相关系数； v /cuσ 是概率分布的不均匀系数。 表 2 细观参数取值 Table 2 Values of mesoscopic parameters u ϕ / n k /Nm-1 s k /Nm-1 E u /GPa f u /MPa v c rho 23 111097.510930 120 0.15 0.8 3 宏观力学响应分析 真三轴数值试验时采用相同的数值试样，以及同 样的加载方式和加载速率，排除设备、人为因素的干 扰。分别进行了不同围压（0.2，0.3，0.4 MPa）下的 等b试验，每个围压下b的取值分别为0，0.25，0.5， 0.75，1.0。 3.1 堆石体试样的变形 堆石体试样在加载过程中的变形如图4（以围压 为0.3 MPa，应力比参数0.5b 为例）所示，试样在 加载过程中，随着加载的进行，颗粒不断旋转、滑移， 结构不断调整，试样在整体上表现为小主应力方向发 生膨胀，中主应力方向压缩量很小。 图 4 堆石体试样的变形 Fig. 4 Deation of numerical samples 3.2 宏观应力–应变关系 堆石试样在各围压 3 σ和不同b值下的偏应力 13 σσ−与轴向应变 1 ε关系曲线如图5所示。由图可 见不同围压下的曲线规律相似，初始加载阶段，堆 石体进一步密实，偏应力随轴向应变增长很快；随着 加载的继续进行，偏应力和轴向应变都增大，并呈现 明显的非线性关系，达到峰值后出现应变软化，偏应 第 4 期 周 伟，等. 基于随机散粒体模型的堆石体真三轴数值试验研究 751 力有所减小；在一定围压 3 σ下，偏应力 13 σσ−和轴 向应变 1 ε曲线随着b值增大逐渐变陡，峰值强度也会 有所增加；同一b值情况下，围压越大曲线越陡，峰值 强度也会增加。 图 5 不同围压下轴向应变–偏应力曲线 Fig. 5 Deviational stress-axial strain curves under different .confining pressures 本文统一规定应变为正表示压缩，负表示膨胀。 中主应变 2 ε和轴向应变 1 ε关系如图6，从图中看出， 中主应变 2 ε随着b值的增大由膨胀转为压缩，在 0.25b 以后，试样在中主应力方向处于膨胀状态。 小主应变 3 ε与轴向应变 1 ε关系曲线如图7所示， 3 ε随 着b值的增加逐渐减小， 且均为负值， 表现为膨胀特性。 图 6 不同围压下 1 ε和 2 ε关系曲线 Fig. 6 Relationship between major principal strain and intermediate principal strain under different confining pressures 图 7 不同围压下 1 ε和 2 ε关系曲线 Fig. 7 Relationship between major principal strain and minor principal strain under different confining pressures 752 岩 土 工 程 学 报 2012 年 从体积应变 v ε和轴向应变 1 ε关系曲线（如图8） 可以看出，不同围压下的曲线变化规律基本相同，随 着加载的进行，堆石体从剪缩转为剪胀，围压越大， 堆石体从剪缩转为剪胀所对应的轴向应变 1 ε越大，剪 胀变化越缓慢，围压越低，试样剪胀越强烈；同一围 压下，b从0到1变化过程中， 堆石体从剪缩转为剪胀 所对应的轴向应变 1 ε减小，剪胀特性更加明显。 图 8 不同围压下 v ε和 1 ε曲线 Fig. 8 Relationship between volumetric strain and major principal strain under different confining pressures 综合分析数值模拟结果，得出以下几点认识偏 应力与轴向应变呈明显的非线性关系，达到强度峰值 后出现应变软化；在一定围压下，b值越大，偏应力峰 值出现越早，曲线越陡，小主应变方向由压缩转为膨 胀，小主应变方向则不断压缩；相同b值下，随着围压 的增加，偏应力曲线明显变陡，峰值强度增加， 2 ε– 1 ε和 3 ε– 1 ε应变曲线变化规律相似。 试样在加载过程 中由压缩逐渐变为膨胀，围压越高，b值越大，剪缩越 明显，剪缩阶段，颗粒间咬合作用增大，试样偏应力 强度增加较快，峰值强度后试样出现软化。 3.3 内摩擦角 内摩擦角ϕ由Mohr-Coulomb强度理论表达式为 13f 13 f arcsin σσ ϕ σσ − ， 4 式中， 1f σ， 3f σ分别为试样的峰值大主应力和小主应 力。 数值模拟得到的堆石体摩擦角ϕ和应力比b的关 系曲线如图9所示，可以看出根据Lade-Duncan破坏 准则所得到的堆石体内摩擦角与数值模拟的变化规律 基本相似，Lade-Duncan破坏准则比较适合描述堆石 体在真三轴应力状态下的强度规律。从图中看出，在 一定围压下，内摩擦角ϕ随着b的增大而增大，在初 始阶段ϕ值增长较快，b为0.25以后ϕ值增长较慢； 同时，在一定b值下，围压越高，内摩擦角ϕ值越小， 在加载初期增长越慢。 图 9 内摩擦角ϕ和应力比参数 b 曲线 Fig. 9 Relationship between friction angle and parameter b under .different confining pressures 4 细观力学响应分析 4.1 各向异性演化规律 第 4 期 周 伟，等. 基于随机散粒体模型的堆石体真三轴数值试验研究 753 对于堆石颗粒，其宏观变形和强度特性与其与 堆石的细观组构的演化规律密切相关。以围压 3 0.3 MPaσ，应力比参数0.5b 加载路径下试样为 例，分析堆石在加载初始 1 0ε、峰值强度时 1 7ε以及加载完成 1 15ε后的细观组构各向 异性演化， 为了便于观察堆石体的各向异性演化规律， 将其投影到YZ平面上，如图11所示，分别给出了在 这3个阶段中颗粒的接触法向、法向接触力的各向异 性玫瑰图，以及相应的傅立叶函数拟合结果（图例为 E和fn），其中a是傅立叶系数，反映了细观组构各 向异性的程度， n θ为细观组构参数各向异性的主方 向。图中每2绘制一个区间，统计接触法向和法向接 触力，图中的接触法向表示落入该角度区间的接触点 总数，法向接触力为落入该角度区间内接触点处法向 接触力的平均值，并投影到YZ平面。 分析图10可知，加载初始阶段 1 0ε，接触 法向的主方向位于水平向， 接触法向的曲线近似圆形， 表示加载初期，在YZ平面上接触法向基本是各向同 性的；在加载过程中，接触法向主方向逐渐向竖直方 向偏转，竖直向接触数目增多，接触法向曲线和拟合 曲线呈“花生状”图形，呈现出明显的各向异性，从 拟合的参数a和θ也可以得出相同的结论；加载完成 后，接触法向曲线和峰值强度时的曲线相似，接触数 目有所减少。法向接触力在加载过程中主方向角由水 平向逐渐倾向竖直方向，堆石颗粒试样到达峰值强度 时，法向接触力明显增大，曲线呈“花生状” ，呈现明 显各向异性，加载完成后法向接触力有所减小。 图 10 颗粒接触法向、法向接触力各向异性演化 Fig. 10 Rose diagram of contact orientation anisotropy and contact .normal force anisotropy 图11给出了颗粒间法向接触力各向异性系数的 变化规律，从图中可以看出加载过程中，颗粒间法 向接触力各向异性程度先增强后又逐渐平稳，堆石试 样的宏观强度变化规律保持较好的一致性，各向异性 系数达到最大值时的轴向应变与峰值偏应力对应的轴 向应变基本相同，试样的宏观强度与法向接触存在内 在的关联。 图 11 颗粒法向接触力各向异性系数 Fig. 11 Variation of contact normal force anisotropy coefficient during loading process 4.2 颗粒配位数 图12表示在加载过程中颗粒配位数的变化曲线。 可以看出，随着加载的进行，配位数先明显增大，后 期逐渐减小，这主要是由于在加载初期，堆石试样进 一步密实，配位数增加；加载后期，堆石发生明显剪 胀，配位数随之减小；相同围压下，b值越大，颗粒接 触越紧密，颗粒配位数越大；围压越高，试样越密实， 颗粒配位数越大。相同轴向应变 1 ε，b值越大，颗粒 配位数越大，围压越高，颗粒配位数越大。 754 岩 土 工 程 学 报 2012 年 图 12 不同围压下颗粒配位数和轴向应变 1 ε曲线 Fig. 12 Relationship between coordination number and axial strain .under different confining pressures 5 结 论 本文在变形体离散元法的基础上，对堆石体进行 真三轴数值模拟试验， 采用等应力比参数的加载路径， 研究堆石体在三向不等应力条件下的应力变形特性， 得到了以下成果和结论 （1） 在围压 3 σ一定条件下， 应力比参数b值越大， 偏应力 13 σσ−和轴向应变 1 ε曲线增长越快，达到峰 值强度后出现应变软化，0b 时， 即为常规三轴试验， 强度最小；b值大小对主应变有显著影响， 中主应变方 向在b值从0增加到1过程中， 会先膨胀再压缩， 小主 应变方向会一直处于膨胀。在加载过程中，堆石试样 会先发生剪缩后转为剪胀，剪缩阶段，堆石强度不断 增加。在1b 时，试样的 2 ε和 3 ε变化显著加快，此时 试样的剪胀特性加剧，试样会在较短时间里达到峰值 强度，并出现应变软化。 （1）围压一定条件下，内摩擦角ϕ随b值的增加 而增加，ϕ值在0.25b 后时增长变缓；同时，围压 越高，ϕ值越小；试验ϕ值变化规律基本符合 Lade-Duncan破坏准则，但普遍比基于Lade-Duncan破 坏准则所得内摩擦角小。 （2）b值一定时，随着围压的增加，偏应力峰值 会明显增大， 2 ε和 3 ε变化规律相似， 最大应变值相近； 围压越高，试样剪缩越明显，出现剪胀时对应的轴向 应变 1 ε越大。 （3）在细观层面上，围压越高，b值越大，颗粒 的配位数越大，颗粒间的接触作用更显著，从而使堆 石试样的强度相应增加。加载过程中，接触法向和法 向接触力的主方向均从水平开始倾向坚直加载方向， 各向异性程度也显著增强，加载完成后，接触法向和 法向接触力相应减小，各向异性程度有所减弱。 （4） 根据各向异性参数的变化规律再结合配位数 的变化情况，表明堆石试样宏观强度的变化与颗粒内 部法向接触力及各向异性程度存在内在关联。 参考文献 [1] 朱思哲, 刘 虔, 包承纲. 三轴试验原理与应用技术[M]. 中国电力出版社, 2003. ZHU Si-zhe, LIU Qian, BAO Cheng-gang. Triaxial test principle and application technology[M]. Beijing China Electric Power Press, 2003. in Chinese [2] 李广信. 土的三维本构关系的探讨与模型难[D]. 北京 清 华大学水利水电工程系, 1985. LI Guang-xin. A study of three-dimensional constitutive relationship of soils and an examination of various models[D]. Beijing Hydraulic Engineering Department, Tsinghua University, 1985. in Chinese [3] 朱俊高, 卢海华, 殷宗泽. 土体侧向变形性状的真三轴试 验研究[J]. 河海大学学报, 1995, 236 28–33. ZHU Jun-gao, LU Hai-hua, YIN Zong-ze. Lateral deation of soil in true triaxial test[J]. Journal of Hohai University, 1995, 236 28–33. in Chinese [4] 徐志伟, 殷宗泽. 粉砂侧向变形特性的真三轴试验研究[J]. 岩石力学与工程学报, 2000, 195 626–629. XU Zhi-wei, YIN Zong-ze. Study on deation characteristic of silt by true triaxial test[J]. Chinese Journal of Rock Mechanics and Engineering, 2000, 195 626–629. in Chinese [5] SHI Wei-cheng, ZHU Jun-gao. Strength and deation beation behaviour of coarse-grained soil by true triaxial tests[J]. J Cent South Univ Technol, 2010, 17 1095–1102. [6] 施维成, 朱俊高, 刘汉龙. 粗粒土应力诱导各向异性真三 轴试验研究[J]. 岩土工程学报, 2008, 3010 810–814. SHI Wei-cheng, ZHU Jun-gao, LIU Han-long. Stress induced anisotropic behaviour of coarse-grained soil by true triaxial tests[J]. Chinese Journal of Geotechnical Engineering, 2008, 3010 810–814. in Chinese [7] KO H, SCOTT R F. Deation sand and failure[J]. ASCE, 1968, 944 883–898. [8] LADE P V, DUNCAN J M. Cubical triaxial tests on cohesionless soil[J]. Journal of the Soil Mechanics and Foundations Division, 1973, 9910 793–812. 第 4 期 周 伟，等. 基于随机散粒体模型的堆石体真三轴数值试验研究 755 [9] LADE P V. Assessment of test data for selection of 3-D failure criterion for sand[J]. International Journal for Numerical and Analytical s in Geomechanics, 2006, 30 307–333. [10] LADE P V, DUNCAN J M..Elastoplastic stress-strain theory for cohesionless soil[J]. Journal of Geotech Engng Div, ASCE, 1975, 101 1037–1053. [11] NAKAI T, MATSUOKA H. Shear behavior of sand and clay under three-Dimensional stress condition[J]. Soils and Foundations, Japanese Siciety of siol Mechanics and Foundation Engineering, 1983, 23 27–42. [12] MATSUOKA h, NAKAI T. Stress-deation and strength characteristics of soil under three different principal stresses[C]// Proc Japan Soc Civil Engrs, 1974, 232 59–70. [13] YIN Jian-hua, CHENG Chun-man, KUMRUZZAMAN M D. New mixed boundary, true triaxial loading device for testing three-dimensional stress-strain-strength behaviour of geoma- terials[J]. Can Geotech, 2010, 47 1–15. [14] CALLISTO L, GAJO A, WOOD Muir. Simulation of triaxial and true triaxial tests on natural and reconstituted Pisa clay [J]. Gotechnique, 2002, 52 649–666. [15] THORNTON C. Numerical simulations of deviatoric shear deation of granular medial[J]. Gotechnique, 2000, 50 43–53. [16] BELHEINE N, PLASSIARD J P. Numerical simulation of drained triaxial test using 3D discrete element modeling[J]. Computers and Geotechnics, 2009, 36 320–331. [17] TANG-TAT N G. Behavior of gravity deposited granular material under different stress paths[J]. Canadian Geotechnical Journal, 2005, 42 1644–1655. [18] TANG-TAT N G. Macro-and micro-behavior of granular materials under different sample preparation s and stress paths[J]. International Journal of Solids and Structures, 2004, 4121 5871–5884. [19] 周 伟, 常晓林, 周创兵, 等. 堆石体应力变形细观模拟 的随机散粒体不连续变形模型及其应用[J]. 岩石力学与工 程学报, 2009, 283 491–499. ZHOU Wei, CHANG Xiao-lin, ZHOU Chuang-bing, et al. Stochastic granule discontinuous deation model of rockfill and its application[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 283 491–4995. in Chinese [20] CHANG X L, ZHOU W. A contact model on basis of the augmented Lagrange and engineering a