二维理想颗粒受剪行为的细观机理*

郭 鸿，陈 茜，李 军，陈栋梁

(1.陕西理工大学土木工程与建筑学院，陕西 汉中 723000；2. 西北农林科技大学水利水电科学研究院，陕西 杨凌 7121001)

二维理想颗粒受剪行为的细观机理*

郭 鸿1,2，陈 茜1，李 军1，陈栋梁1

(1.陕西理工大学土木工程与建筑学院，陕西 汉中 723000；2. 西北农林科技大学水利水电科学研究院，陕西 杨凌 7121001)

通过对二维理想颗粒(圆盘)的直接剪切细观力学机理分析，研究了颗粒受剪时的运动规律，建立了能够描述等径颗粒情况下最密排列颗粒的最大剪应力公式。并通过对另外四种不同级配的圆盘颗粒的直接剪切离散元模拟，探讨了所建立的最大剪应力公式的应用。结果表明：在竖向荷载一定时，最密排列颗粒的最大剪应力取决于边界颗粒影响参数和非边界颗粒影响参数，且这两个参数在最密排列等径颗粒的情况下数值保持一致；对于一般颗粒级配，边界影响参数基本保持不变，非边界影响参数随着d10，d30和d60的增加而线性减小。

颗粒材料；离散元法；直接剪切；颗粒级配

颗粒材料广泛的存在于自然界中，在岩土工程、道路工程以及铁路工程中的有着极广泛的应用，滑坡、地基剪切失稳等现象都与颗粒材料的剪切特性有直接关系[1-4]。目前的对于颗粒材料剪切特性的研究大致有两种方法，即宏观室内试验以及细观数值模拟[5-9]。近年来兴起的离散元方法对研究颗粒材料起到了很大的推动作用，其基本方法是针对材料定义不同的颗粒接触模型，利用离散元模拟软件来探索其细观指标诸如力链、配位数、孔隙率、局部剪应力和位移场等对宏观相应的影响[10-13]。然而，由于目前还少有文献在颗粒材料在受剪时，将其单个颗粒的受力及运动特征与宏观层面剪切应力建立定量关系，所以针对不同的颗粒材料，仍然需要试验测定或离散元模拟的方法确定其剪切特性，诸如内摩擦角等。理想圆形颗粒是颗粒材料最简化的物理模型，分析其受力规律有望会对实体颗粒材料的研究提供有效的思路，目前已有研究针对理想颗粒的力学特性研究[14-19]，主要探讨颗粒形状、级配等因素对剪切特性的影响，但是存在的问题仍然是尚未建立颗粒细观受力与宏观相应之间的定量关系，基本都是从力链演化等角度定性描述。

基于此，本文以二维等径颗粒为研究对象，通过分析单颗粒在剪切过程中的受力平衡，推导了二维等径颗粒规则排列时的剪力模型。并且运用基于PFC2D的数值剪切试验，对所建立的模型，针对其他不同级配的二维理想颗粒，进行了应用探讨和验证。

1 颗粒剪切运动受力分析

1.1 边界颗粒

颗粒受剪时，同时存在法向力和切向剪力。分析单颗粒在受剪时的受力规律，如图1(a)、(b)和(c)所示。颗粒在边界时的初始状态，上部颗粒处于平衡稳定状态，随着剪力Fi的增大，颗粒发生缓慢逆时针旋转，同时逐渐离开右下方所接触的颗粒，当刚好离开右下方的颗粒时，剪力达到最大值，随后剪力减小，颗粒“跳跃”至前方另外两个颗粒之间，此时颗粒材料已经剪坏。竖向荷载作用在边界单个颗粒上的分力Pi仍然是竖直分向，其值为

Pi=pD/m

(1)

式中p为竖向总应力，D为试样水平方向总宽度，m为水平方向的颗粒总数。

当单个颗粒在同时受到Pi和Fi作用下，它的运动从初始稳定点逐渐运动至另一个稳定点，在该运动过程中，颗粒被剪坏。如图1(d)受力分析，这种情况时其水平方向的动态平衡方程为

(2)

式中N为下部颗粒对其产生的法向支持力，μ为颗粒之间的摩擦系数，θ为该颗粒和下部颗粒几何中心连线与水平方向的夹角。

竖直方向的平衡方程为

(3)

式中m为单个颗粒的质量，g为重力加速度。

由式(2)和式(3)得出边界颗粒情况下，单个颗粒移动所需要的水平剪力为

(4)

图1 边界颗粒受剪运动过程及受力分析Fig.1 Movement and force analysis of singe particle in the boundary under shear force

1.2 非边界颗粒颗粒

另一种情况就是颗粒在非边界上，其剪切运动过程如图2(a)、(b)和(c)所示。竖向荷载作用在非边界单个颗粒上的分力P′不在竖直方向，假设左上角相邻颗粒与该单个颗粒同步剪切运动，则其值为

(5)

图2 非边界颗粒受剪运动过程及受力分析Fig.2 Movement and force analysis of singe particle within the boundary under shear force

对于此种情况来说，当单个颗粒在同时受到P′和Fi作用时，它的运动从初始稳定点逐渐运动至另一个稳定点，受力分析如图2(d)示，此时，水平方向的动态平衡方程为

(6)

竖直方向的平衡方程为

(7)

由式(6)和式(7)得出非边界颗粒情况下，单个颗粒移动所需要的水平剪力为

(8)

1.3 一般情况

推广到一般的等径颗粒最密排列情景，如图3所示，若每层m个颗粒，考虑到颗粒数目的影响效应，则要使该层的所有颗粒都同步运动，则对于每层来说需要的的剪力为

F=fjFi

(9)

式中fj(j=1，2)是与边界条件有关的系数，f1表示边界层颗粒影响系数，f2表示非边界层颗粒影响系数。

假设有一共有2n层颗粒，所有层总剪力为

(10)

得到的剪应力τ为

(11)

理想情况下，规则排列的刚性颗粒在受剪时其运动是同步的，然而考虑到颗粒之间非刚性接触，其运动并非同步进行，且每层不同颗粒运动相应存在一定的相位差。由于非同步现象涉及到颗粒之间的接触本构，且影响因素诸多，故为了简单起见，本文只考虑边界层颗粒影响系数f1和非边界层颗粒影响系数f2。

2 二维等径颗粒DEM剪切试验

建立二维等径DEM(离散元)数值剪切试样，如图4所示，所有试样都是等径颗粒最密排列。保持理想试样的尺寸为6cm×2cm，设置不同粒径的二维颗粒，即5，2.5，1.667，1.25，1，0.833，0.714和0.625mm。设定颗粒间的摩擦系数为1，颗粒接触为最简线弹性模型，墙体刚度和颗粒刚度分别为500MN/m和50MN/m。

图3 m×2n理想剪切试样Fig.3 m×2nideal shear sample

图4 不同粒径的二维离散元剪切试样Fig.4 2D DEM shear sample of different particle sizes

对不同粒径的试样做离散元数值剪切试验，编写自动加压伺服程序段，使得在试样剪切前颗粒体上部的竖向压应力为100kPa，然后固定数值剪切盒的下半部分，向左移动数值剪切盒的上半部分。需要说明的是，为了提高模拟计算速度，虚拟剪切速度采用2cm/s，在下文中所有的数值试验均采用该剪切速度。

离散元数值剪切试验得到的内摩擦角和半径的关系如图5所示。由此图可以看出，粒径对摩擦角的影响十分显著，且根据拟合公式，二者基本呈线性关系。当竖直荷载为100kPa时，数值剪切试验得到的最大剪应力和颗粒半径的关系如图6所示，菱形点为DEM试验值，方形点数据为本文式(11)的计算值。在本文所设定的计算参数和DEM参数的情况下，f1=2.4,f2=3.3。说明非边界颗粒影响系数对最大剪切力的影响较大，这一方面是因为非边界颗粒的平均配位数较大，另一方面是由于内部颗粒在剪切的过程中颗粒相互咬合更多(而边界颗粒的相对滑动更多)。换言之就是非边界颗粒相对边界颗粒的运动阻力更大，而各个颗粒的定向运动阻力之和和宏观上的剪切力是平衡的，因此f2的数值大于f1。

图5 颗粒半径和内摩擦角的关系Fig.5 Rrelation of particle radius vs internal friction angle

图6 颗粒半径和最大剪应力的关系Fig.6 Relation of particle radius vs maximum shear stress

3 非规则颗粒数值试验

实际工程中存在的材料大都属于非等径且具有复杂级配的颗粒介质。据此，本节以四种不同级配的颗粒材料为对象，进行数值剪切试验，进一步探索颗粒大小级配对最大剪应力的影响规律，同时讨论二维非等径理想颗粒情况下式(11)的应用。如图7所示，为四种级配的二维颗粒(圆盘)，采用Jiang氏方法[20]分四层生成。颗粒级配曲线如图8所示，Cc和Cu值如表1所示。

表1 四种级配砂土的Cc和Cu值

在处理非等径颗粒情况时，先计算出每种颗粒级配情况下的加权平均粒径

图7 四种颗粒级配的离散元模型Fig.7 Discrete model of the 4 types sands

图8 四种级配试样Fig.8 Four DEM samples of different particle size distributions

(12)

式中Ri代表某种级配情况下第i级的颗粒半径，Ni代表该i级颗粒的数目，n为不同粒径的总级数。

然后将处理后的等径颗粒按照矩形方式排列为2n行×m列，再用式(10)进行计算最大剪切力。在以上四种颗粒级配情况下，f1=1，f2与d10，d30和d60有关，如图9所示。可以看出，f2基本都随d10，d30和d60的增大而减小，且其线性关系比较明显。需要进一步说明的是，f1是与边界颗粒有关的系数，f2是与非边界颗粒有关的系数，所以当数值试样内的颗粒数量增大时，f1的影响相对会减弱，而f2的影响相对会加强。

图9 d10，d30和d60和f2的关系曲线Fig.9 Relation of d10，d30 and d60 vs f2

这是因为在数值试样内部，随着颗粒数量的增加，水平方向的力链数量更多，形状更细，更弯曲，如图10所示，而对于边界上的颗粒来说，由于上部竖向均布压应力，颗粒的排列相对规则，形态变化也相对不明显。f2和最大剪应力的关系如图10所示，可以看出最大剪应力随着f2的增大而线性增大，表明了力链越密集，剪应力越大(见图11)。

图10 f2和τ的关系曲线Fig.10 Relation of f2 vs τ

图11 四种级配离散元剪切过程力链图Fig.11 Force chains of four samples with different particle size distributions in DEM shear test

4 结 论

本文基于二维等径颗粒体，分析了颗粒材料在受剪时的力学机理，提出了颗粒规则最密规则排列时的最大剪应力公式。通过对四种不同级配颗粒体的数值剪切试验，发现本文提出的最大剪应力公式仍然适用，所不同的只是f1和f2，其中，两者在非等径颗粒情况下的数值都较等径颗粒情况下的数值小，而且后者随着d10，d30和d60的增大而减小。f2和最大剪应力的关系表明，力链形态越密集，剪应力越大。

[1] 郭鸿,骆亚生,程大伟,等. 倾斜推移体激发颗粒流的连续介质模型[J]. 中山大学学报(自然科学版),2015,54(6):94-98.

[2] 李洋溢. 条形基础加筋砂土地基室内模型试验的分析研究[D].上海：同济大学,2006.

[3] 曹延波,彭芳乐,小竹望,等. 加筋砂土边坡渐进性变形破坏的数值分析[J]. 岩石力学与工程学报,2010,29(S2):3905-3915.

[4] 吴剑. 滑带剪切过程的离散元模拟研究[D].北京：中国科学院研究生院,2007.

[5] 严颖,季顺迎. 颗粒形态对离散介质剪切强度的影响[J]. 岩土力学,2009,30(S1):225-230，249.

[6] 边学成,李伟,李公羽,等. 基于颗粒真实几何形状的铁路道砟剪切过程三维离散元分析[J]. 工程力学,2015,32(5):64-75，83.

[7] 秦建敏,张洪武. 颗粒材料的微观临界状态理论模型[J]. 岩土力学,2010,31(12):3697-3703.

[8] 秦建敏,迟璐璐. 颗粒材料剪胀性的微观力学分析[J]. 岩土力学,2013,34(5):1508-1514.

[9] 孙其诚,程晓辉,季顺迎,等. 岩土类颗粒物质宏-细观力学研究进展[J]. 力学进展,2011,41(3):351-371.

[10] 张宇,余飞,陈善雄,等.CAS-1模拟月壤动剪切模量与阻尼比的试验研究[J]. 岩土力学,2014,35(1):74-82.

[11] 蒋明镜,彭镝,申志福,等. 深海能源土剪切带形成机理离散元分析[J]. 岩土工程学报,2014,36(9):1624-1630.

[12] 蒋明镜,张望城,王剑锋. 密实散粒体剪切破坏能量演化的离散元模拟[J]. 岩土力学,2013,34(2):551-558.

[13] 刘动,陈晓平. 滑带土环剪剪切面的微观观测与分析[J]. 岩石力学与工程学报,2013,32(9):1827-1834.

[14] 李罡,刘映晶,尹振宇,等. 粒状材料临界状态的颗粒级配效应[J]. 岩土工程学报,2014,36(3):452-457.

[15] 孔亮,季亮亮,曹杰峰. 应力路径和颗粒级配对砂土变形影响的细观机制[J]. 岩石力学与工程学报,2013,32(11):2334-2341.

[16] 陆勇,周国庆,赖泽金. 砂土剪切特性的荷载与粒径效应研究[J]. 中国矿业大学学报,2014,43(2):195-202.

[17]DABEETA,WIJEWICKREMED,BYRNEP.Applicationofdiscreteelementmodelingforsimulationofcyclicdirectsimpleshearresponseofgranularmaterials[C]//NationalConferenceonEarthquakeEngineering， 2014: 21-25.

[18]CAVARRETTAI,O’SULLIVANC,IBRAIME,etal.CharacterizationofartificialsphericalparticlesforDEMvalidationstudies[J].Particuology, 2012, 10(2): 209-220.

[19]O’SULLIVANC.Particulatediscreteelementmodelling[M].Taylor&Francis, 2011.

[20]JIANGMingjing,KONRADJM,LEROUEILS.AnefficienttechniqueforgeneratinghomogeneousspecimensforDEMstudies[J].ComputersandGeotechnics, 2003, 30 (5): 579-597.

The microscopic mechanism of the shear properties of two dimensional ideal particle

GUOHong1,2,CHENXi1，LIJun1,CHENDongliang1

(1. School of Civil Engineering and Architecture, Shaanxi University of Technology,Hanzhong 723000, China；2. Institute of Water Resources and Hydropower Research, Northwest A&F University,Yangling 712100, China)

By analyzing the micromechanics of the two-dimensional ideal particles (disc) under direct shear test, the movement during the shearing is investigated. The formula described the maximum shear stress in the regular dense packing condition is established. In addition, following the direct shear simulations of the four samples with different particle size distributions, application of the formula which could predict the maximum shear stress was discussed. The results show that the maximum shear stress is depending on two parameters namely boundary effect parameter and the non-boundary effect parameters, and these two parameters keep the same under the regular dense packing when the vertical loading is constant. For the general packing, the boundary effect parameters remained unchanged, and the non-boundary effect parameter is linearly decreases with thed10,d30andd60.

particulate material; discrete element method; direct shear; particle size distribution

10.13471/j.cnki.acta.snus.2016.04.008

2015-10-26

国家自然科学基金资助项目(51178392) ；陕西理工大学科研资助项目(SLGKY15-24)

郭鸿(1984年生)，男;研究方向：离散元数值模拟；通讯作者： 骆亚生;E-mail:lyas1967@nwsuaf.edu.cn

TU442

A

0529-6579(2016)04-0047-06