粗粒料工程力学性质的细观模拟

基金项目:新疆维吾尔自治区自然科学基金面上项目（2013211A027）

作者简介:齐阳（1988），男，河南舞阳人，硕士研究生，主要从事水利水电工程数值分析研究。Email:yqi_qy@163.com

通讯作者：唐新军（1959），男，陕西蒲城人，教授，博士生导师，主要从事当地材料坝设计理论教学与科研工作。Email:tangxj59@sina.comDOI：10.13476/j.cnki.nsbdqk.2014.04.024

摘要：以三维离散元颗粒流理论为基础，PFC3D_EV为工具，改进内置程序代码，对比粗粒料室内三轴固结排水试验结果，生成颗粒接触点数均匀、各向同性应力相同的三轴试验数值模型。引入clump颗粒提高了数值试样内部颗粒形状的复杂度，对比了不同围压下数值试验应力-应变曲线、变形曲线与室内试验的差异，并探讨了剪切带发展规律及颗粒形状对粗粒料强度、变形特性的影响。结果表明：颗粒形状是影响粗粒料工程力学特性的主要因素，提高试样内部颗粒形状的复杂性可获得与物理试验拟合较好的应力-应变曲线；数值模型暂无法实现颗粒破碎状态，造成大应变时剪胀偏大；位移场发展变化过程显示应变软化加速了剪切带的形成；高围压下剪切带明显，且围压越高厚度越薄。

关键词：颗粒流；粗粒料；强度；剪切位移场；细观模拟

中图分类号：TU431文献标志码：A文章编号：16721683（2014）04011105

Mesoscopic Simulation of Mechanical Behaviors of Coarse Grained Materials

QI Yang,TANG Xinjun,LI Xiaoqing

(School of Water Conservancy and Civil Engineering,Xinjiang Agricultural University,Urumqi 830052，China)

Abstract:Based on the results of laboratory triaxial consolidateddrained tests on coarse grained materials and threedimensional discrete particle flow theory,a numerical model of triaxial test for coarse grained materials with uniform particle contact points and same isotropic stresses was developed using the tool of PFC3D_EV and improvement of program code.Clump particles were introduced to improve the complexity of internal particle shape of the sample.The differences of stressstrain curves and volumetric strain  axial strain curves of triaxial tests under different confining pressures were compared,and the impacts of particle shape and the development law of shear bands on the strength and deformation characteristics of coarse grained materials were discussed.The results showed that (1) particle shape is the main factor on the mechanical behavior of coarse grained materials.The stressstrain curve obtained from the numerical simulation matches well with the experimental results if the complexity of internal particle shape is improved;(2) the shear dilatancy of numerical simulation is a little larger because the particle deformation and breakage cannot be performed in numerical simulation;(3) the development of displacement field indicates that strain softening accelerates the formation of shear bands;and (4) the thickness of shear bands becomes thinner with the decreasing of confining pressure.

Key words:particle flow;coarse grained materials;strength;shear bands;mesoscopic simulation

土石坝工程对坝壳料等粗粒料力学特性的精度要求越来越高。颗粒形状、级配、排列方式（紧密、疏松）、接触方式、颗粒强度等众多因素决定了粗粒料细观结构的复杂性。传统土力学研究方法中基于宏观连续性假设的力学模型不能理想地解释粗粒料等离散颗粒体所表现出的力学响应。但是基于离散元颗粒流理论的数值模拟，能够将细观结构变化和宏观力学响应联系起来，从土体内部细观结构力学的角度，分析土体的力学性质和变形特性。

虽然利用颗粒流仿真模拟土体工程力学性质已成为国内外岩土研究的新趋势［16］，但是，目前对粗粒料大三轴试验的模拟还不多见，邵磊［7］模拟了某工程堆石料的三轴试验，但是没有考虑颗粒形状对粗粒料力学特性的影响；张超［8］虽然考虑了颗粒形状的影响，但是采用了黏结模型，将所有离散颗粒黏结为一个整体，不仅背离了堆石体材料离散性、无黏结性的基本条件，同时也造成了初始弹性模量与室内试验偏差较大。

本文通过改进PFC3D_EV软件内置程序代码，生成clump颗粒，提高数值试样与实际土体在颗粒形状、排列方式等方面的契合度，并在不应用黏结模型的情况下对室内大三轴试验进行了数值模拟，分析剪切带发展规律及其颗粒形状对粗粒料工程力学特性的影响。

1颗粒流的基本理论

颗粒离散元法将材料理想化为相互独立、相互接触、相互作用的颗粒集合，从细观角度对材料的力学特性进行分析。它在计算过程中交替应用颗粒体的牛顿第二定律和接触的力—位移关系，其中牛顿第二定律用来决定每一个颗粒的运动和旋转行为，而力—位移关系则用来更新由每一对接触产生的接触力。无论何时，颗粒体内部力都需处于一种平衡状态。

1.1接触本构模型

坝壳料等粗粒料没有黏结性，属摩擦材料，因此只能应用接触刚度模型和滑动模型来模拟其本构行为。

接触刚度模型用于揭示接触力-位移关系，主要由以下四个向量方程表征：

Fn=Fni+Fsi（1）

Fni=KnUnni（2）

Fsi←Fsi+ΔFsi（3）

ΔFsi=-KsΔUsi（4）

式中：Fn为接触力；Fni、Fsi为接触力法向、切向分量；Kn、Ks为法向、切向接触刚度值；Un为接触实体法向相对位移；ΔUsi为切向位移增量；ni为接触平面单位法向向量。

滑动模型，即接触实体可在剪切强度范围内发生滑动，其滑动条件为

|Fsi|＞Fsmax滑动（5）

Fsmax=μ|Fni|（6）

式中：μ为接触摩擦系数，其值为两接触实体摩擦系数中的较小值，其它符号同前。

1.2接触本构模型的适用性

用球体单元代表颗粒介质是PFC模型亟需克服的缺陷［911］。周健［12］、罗勇［13］等人将砂粒简化为球体颗粒，应用上述本构模型进行模拟试验，证明是可行的。但是堆石体具有颗粒形状不规则，排列紧密，咬合力大，抗剪强度高，易破碎等特点，仅应用上述模型很难达到实际应有的抗剪强度。为解决此问题，邵磊［7］、张超［8］等人采用了引入黏结模型的方法，虽然获得了与室内试验相接近的加载曲线，但该模型并不能反映粗粒料实际的细观结构变化。本文尝试生成clump颗粒，从问题的根源出发，通过提高试样内部颗粒形状和排列方式的复杂性来解决数值模拟试验颗粒间咬合力弱的问题。

2某工程坝壳料三轴试验的数值模拟

2.1室内试验数据

对某工程坝壳料进行了大三轴固结排水试验，试样尺寸为：高600 mm、直径300 mm，制样干密度通过相对密度试验测得为2.16 g/cm3，试样级配曲线见图1。

图1室内三轴试验及数值模拟试验试样颗粒级配曲线

Fig.1Gradation composition curves for the laboratory triaxial

experiment and numerical simulation

2.2数值模拟试验过程

2.2.1试样生成过程

为了探讨颗粒形状的影响，数值模型生成两种细观参数相同的试样，纯球体颗粒试样和clump试样（球体颗粒、clump颗粒混掺）。两种试样的生成步骤如下：

（1）生成压力室。压力室由上下加载墙和控制围压的圆筒墙构成。加载墙为刚体墙，其刚度稍大于颗粒单元的刚度；圆筒为柔性墙，其刚度值一般设定为颗粒刚度值的1/10。试样生成阶段控制压力室高度为600 mm，直径为300 mm。

（2）按级配曲线生成特定孔隙率的试样。为提高计算效率，控制数值试样最小颗粒直径为1 mm，具体级配曲线见图1。在有限空间直接生成标准试样往往会造成颗粒溢出，先缩减再膨胀可避免此问题。纯球体试样由以下步骤生成：确定半径缩减系数并指定孔隙率；计算颗粒集合体总体积并按级配曲线确定各粒径区间体积；用区间中间半径估计该区间颗粒个数，参考该个数反复调试，直到按半径均匀分布生成的该组颗粒满足体积百分比要求；按缩减系数将各颗粒半径膨胀至标准尺寸。clump试样则是在纯球体试样的基础上将部分球体颗粒等体积替换为clump颗粒获得。

（3）设定初始各向同性应力值。该步主要用来削弱模拟过程中（卸荷、细观参数赋值等阶段）内嵌应力的影响。各向同性应力值由三向轴向应力的均值确定，该值宜设为无侧限抗压强度的百分之一。

（4）减少“浮点”颗粒个数。“浮点”颗粒定义为接触点数小于3的颗粒。本步主要用来增加颗粒平均接触点数，提升试样内部结构的均匀性与稳定性。

（5）细观参数赋值后，运行至颗粒平衡状态即得到最终数值试样。

压力室与clump试样见图2，其中深色部分为clump颗粒，浅色部分为球体颗粒。

图2clump试样及压力室

Fig.2Clump sample and pressure chamber

2.2.2试样剪切过程

PFC3D_EV通过伺服系统控制上下加载墙及柔性圆筒墙，实现试样的等压固结。围压达到指定值后，加载墙脱离伺服系统的控制开始以应变控制式的方式对试样进行剪切，剪切速度通过多阶段逐渐递加的方式获得。剪切过程中，圆筒墙在伺服系统的作用下自动控制围压的误差不大于1%。为保证该过程为静力加载过程，剪切速度设置得非常小。

2.2.3细观参数与宏观力学响应之间的联系及细观参数标定（1）接触颗粒间弹性模量Ec 。该参数不仅克服了传统方法中颗粒刚度值不能随半径协同变化的缺点，而且建立了细观弹性模量与宏观弹性模量的直接联系。初始弹性模量随着Ec的增加而增大，标定出的Ec值一般大于宏观初始弹性模量。颗粒刚度值通过式（7）决定。

Kn=4REc（7）

式中：R为两接触颗粒半径的均值。

（2）法向刚度与切向刚度的比值Kn/Ks。该比值一方面用来确定切向刚度值，另一方面影响试样泊松比的大小，泊松比随着该比值的增加而增大。

（3）接触摩擦系数μ。μ主要影响峰值强度的大小，峰值强度随着摩擦系数的增加而增大。罗勇［13］等人亦通过试验发现，摩擦系数增大会加剧剪胀和应变软化，但对初始弹性模量的影响甚微。

（4）孔隙率n。孔隙率是影响试样宏观力学响应的主要因素之一。罗勇［13］等人数值模拟结果显示，低围压下，孔隙率越小，峰值强度越高，应变软化和剪胀越明显。对于数值模拟，颗粒密度为定值，孔隙率越小试样密度越大。上述结论与姜景山［14］由室内试验得出的“密实粗粒土在低围压下应力一应变曲线一般呈软化型，且常常表现出较大的体胀变形”的结论是一致的。

（5）球体、peanut、dyad颗粒在试样中所占的体积百分比v1、v2、v3。引入clump颗粒会成倍增加颗粒个数，降低计算效率，这里只生成如图3所示的两种clump颗粒。clump颗粒一般通过等体积替代球体颗粒的方式生成，即保持peanut、dyad颗粒的体积与被替代球体颗粒的体积相等。clump颗粒的体积为各组成颗粒的体积和减去重叠部分体积。clump颗粒不会出现破碎和开裂的情况。

图3两种clump颗粒

Fig.3Two kinds of clump particles

通过对以上细观参数与宏观力学响应之间联系的分析，经反复标定，得到了一组比较理想的细观参数,见表1。

2.3数值模拟结果分析

2.3.1应力-应变特性分析

由图4可知，室内试验与数值模拟试验的应力-应变曲线具有基本一致的发展趋势，都为双曲线型，且峰值主应力差（σ1-σ3）都随围压的增加而增大。但是clump试样与室内试验加载曲线拟合的更好，因为clump试样颗粒形状更加多样，排列方式更加复杂，颗粒之间咬合力更接近实际情况。由此可知，不应用黏结模型，仅通过提高数值试样与实际土体在颗粒形状、排列方面的相似度，即可获得与室内试验拟合较好的加载曲线。

表1PFC3D_EV数值模拟细观参数

Table 1Mesoscopic parameters in PFC3D_EV simulation

细观

参数Ec

/MPakn/ksμnρ

/（kg·m3）v1v2v3试样194.282.50.90.32 163100试样294.282.50.90.32 1630.350.20.45但大应变时，数值模拟和室内试验加载曲线差异较大，主要表现在数模试验的应力-应变曲线为应变软化型。从细观角度分析，主应力差达到峰值之前，为应变硬化阶段，体积变形处于剪缩状态。该阶段孔隙率不断减小，颗粒平均接触点数增多，颗粒所受限制和约束力持续增加，颗粒间咬合力不断加强。峰值出现之后为应变软化阶段，该阶段处于剪胀加剧状态，孔隙率增大，颗粒平均接触点数减少，颗粒所受限制和约束力变小，颗粒间咬合力变弱，抗剪能力下降。

2.3.2变形特性分析

粗粒料体积应变-轴向应变关系曲线见图5。对于室内试验，粗粒料体积变形有低压剪胀高压剪缩的特点；对于数值模拟试验，两种试样剪胀（缩）特性相差不大，都有随围压增加剪胀变弱的特点，与室内试验一致。但是，各围压下，两数值试样剪胀均偏大。

从细观上分析，当颗粒在荷载的作用下翻滚、跨越时容易造成剪胀，而颗粒在荷载作用下滑动挤压空隙、破碎填充空隙时往往呈现出剪缩特性［15］。低围压下数值试验和室内试验都表现出先剪缩后剪胀的特性，主要是因为剪切初期颗粒在荷载的作用下不断挤压空隙，而剪切后期低围压又不能限制颗粒翻滚、跨越的缘故。室内试验颗粒在高围压下表现出剪缩特性，除了因为高围压能够限制颗粒翻滚、跨越外，最重要的是荷载加大所引起的大量颗粒破碎。而数值试样皆为刚体颗粒，不会出现颗粒破碎现象，高围压虽然能够在一定程度上限制颗粒翻滚但却无法抑制剪胀加剧。本室内试验在600 kPa较高围压下仍表现出剪胀的趋势，与试样中含有较多质地坚硬的圆砾有关。

两数值试样剪胀（缩）特性总体趋势一致，但又略有不同。首先，clump试样剪缩持续时间稍长，主要是因为颗粒形状的复杂性增加了颗粒之间咬合力，延缓了颗粒出现翻滚、跨越现象的时间；其次，剪切后期clump试样相对球形试样剪胀偏大，这是因为clump颗粒形状不规则，颗粒接触点数较多，与周围颗粒相互作用强，一旦开始翻滚、跨越就会带动周围颗粒滚动。颗粒形状、圆度、表面粗糙度、是否易破碎等因素决定了数值试样与室内试样体变曲线存在较大的差异。两数值试样体变曲线总体趋势相同说明了仅靠更改颗粒形状无法完全改变数值试样剪胀偏大的事实。

2.4剪切位移场的变化和发展过程

PFC3D_EV不仅可以实现对加载曲线和变形曲线的理想模拟，还可以记录各组成颗粒的位移变量并以不同的颜色表现出来，可以清楚地观察到剪切过程中试样位移场的发展变化过程，从而探究剪切带的形成规律。图4室内试验与数值试验应力-应变曲线

Fig.4Stressstrain curves of laboratory experiment and numerical simulation

图5室内试验与数值试验εvεa关系曲线

Fig.5εvεa curves of laboratory experiment and numerical simulation图6为某围压下，试样剪切位移场不同阶段的发展变化情况。其中，图6(a)为固结阶段的位移场，该阶段颗粒位移方向、大小毫无规律，处于一种比较混乱的状态，固结结束后位移场归零；图6(b)为峰值强度出现之前的位移场，该阶段剪切位移场关于中心截面对称，并沿着主应力方向均匀变化；图6(c)-图6(e)为峰值强度出现之后的位移场，亦即应变软化阶段的位移场，该阶段细观结构迅速排列变化，颗粒剪切位移场的定向趋势越来越明显，剪切带的形状趋于明显且集中。剪切带形成后，试样沿剪切带分成两个部分，两侧颗粒位移方向基本上与剪切带方向平行。

图6剪切位移场的形成和变化过程

Fig.6Formation and development of shear displacement field

图7为剪切停止后不同围压下的位移场，其中显示，围压较小时，剪切带并不明显，剪切带处位移场比较混乱。围压比较大时候，剪切带明显且有比较规整的边界。比较图7后两幅图可知，剪切带随着围压的增大厚度减小，这主要是因为，围压越大对内部颗粒的限制和约束力越大，剪切破坏影响的范围相对较小。综上可知，高围压下剪切带明显，且围压越高厚度越薄。

图7不同围压下试样最终剪切位移场

Fig.7The final shear displacement fields under

various confining pressures

4结论

（1）对于粗粒料大三轴数值模拟试验，不需采用与实际背离的黏结模型，仅提高数值试样与实际土体在颗粒形状、排列等方面的相似度，即可显著提高抗剪强度，获得与室内试验拟合比较理想的应力-应变曲线。

（2）对比数值试样应力-应变曲线发现，颗粒形状主要影响峰值强度的大小，而对初始弹性模量影响很小。

（3）数值模拟与室内试验体积变形特性差异较大，主要表现在数值模拟试验剪切后期剪胀偏大。这是因为数值试样皆由刚体颗粒组成，不会出现颗粒破碎，且颗粒之间的咬合力弱于实际土体颗粒，轴向压力不断增大促使颗粒翻滚、跨越，造成大应变时剪胀偏大。

（4）峰值强度出现之前，位移场沿主应力方向均匀变化；峰值强度出现后，细观结构迅速排列变化，剪切位移场的定向趋势越来越明显，剪切带的形状趋于明显且集中，即应变软化加速了剪切带的形成。

（5）高围压下剪切带明显，且围压越高剪切带厚度越薄。

参考文献（References）：

［1］Cundall P A,Strack O D L.A Discrete Numerical Model for Granular Sssemblies［J］.Geotechnique,1979,29(1):4765.

［2］Pierce M E,Cundall P A,Van Hout G J,et al.PFC3D Modeling of Caved rock under Draw［J］.Numerical Modeling inMicromechanics via Particle Methods,2003:211217.

［3］Deluzarche R,Cambou B.Discrete numerical modelling of rockfill dams［J］.International Journal for Numerical and AnalyticalMethods in Geomechanics,2006,30(11):10751096.

［4］Hosseininia E S,Mirghasemi A A.Numerical Simulation of Breakage of Twodimensional Polygonshaped Particles usingDiscrete Element Method［J］.Powder Technology,2006,（166）:100112.

［5］周健，史旦达，贾敏才,等.循环加荷条件下饱和砂土液化细观数值模拟［J］.水利学报，2007,38(6):697703.(ZHOU Jian,SHIDanda,JIA Mincai,et al.Micromechanical Simulation of Liquefaction Behavior of Saturated Sand under Cyclic loading［J］.Journalof Hydraulic Engineering,2007,38(6):697703.（in Chinese)）

［6］李凡.岩土材料破损特性的颗粒流研究［J］.土木工程学报，2007,40(9):7881.(LI Fan,A Study on the Breakage Properties ofGeological Materials using Particle Flow Simulation［J］.China Civil Engineering Journal,2007,40(9):7881.（in Chinese))

［7］邵磊，迟世春，贾宇峰.堆石料大三轴试验的细观模拟［J］.岩土力学，2009,30(增刊1）:239243.(SHAO Lei,CHI Shichun,JIA Yufeng.Mesomechanical Simulation of a Large Scale Triaxial test of Rockfill Materials［J］.Rock and Soil Mechanics,2009,30(Supp.1):239443.（in Chinese))

［8］张超，展旭财，杨春和.粗粒料强度及变形特性的细观模拟［J］.岩土力学，2013,34(7):20772083.(ZHANG Chao,ZHANXucai,YANG Chunhe.Mesoscopic Simulation of Strength and Deformation Characteristics of Coarse Grained Materials［J］.Rockand Soil Mechanics,2013,34(7):20772083.（in Chinese))

［9］Potyondy D O,Cundall P A.A Bondedparticle Model for Rock［J］.International Journal of Rock Mechanics and MiningScience,2004,41(8):13291364.

［10］Hoek E ,Brown E T.Practical Estimates of Rock Mass Strength［J］.International Journal of Rock Mechanics and MiningScience,1998,34(8):11651186.

［11］Cho N,Martin C D,Sego D C.A Clump Particle Model for Rock［J］.International Journal of Rock Mechanics and MiningScience,2007,44(7):9971010.

［12］周健，池永.砂土力学性质的细观模拟［J］.岩土力学，2003,24(6):901906.(ZHOU Jian,CHI Yong.MesomechanicalSimulation of Sand Mechanical Properties［J］.Rock and Soil Mechanics,2003,24(6):901906.（in Chinese))

［13］罗勇，龚晓南，连峰.三维离散颗粒单元模拟无黏性土的工程力学性质［J］.岩土工程学报,2008,30(2):292297.（LUOYong,GONG Xiaonan,LIAN Feng.Simulation of Mechanical Behaviors of Granular Materials by Threedimensional Discrete Element Method Based on Particle Flow code［J］.Chinese Journal of Geotechnical Engineering,2008，30(2):292297.（in Chinese))

［14］姜景山，刘汉龙，程展林等.密度和围压对粗粒土力学性质的影响［J］.长江科学院院报，2009,26(8):4650.(JIANG Jingshan,LIU Hanlong,CHENG Zhanlin,et al.Influences of Density and Confining Pressure on Mechanical Properties for CoarseGrained soils［J］.Journal of Yangtze River Scientific Research Institute,2009，26(8):4650.（in Chinese))

［15］迟明杰,赵成刚,李小军.砂土剪胀机理的研究［J］.土木工程学报,2009,42(3):99104.(CHI Mingjie,ZHAO Chenggang,LI Xiaojun.Stressdilation Mechanical of Sands［J］.China Civil Engineering Journal,2009,42(3):99104.（in Chinese))