改进颗粒组构力学模型模拟筒仓卸粮成拱细观机理.PDF

第34卷 第20期 农 业 工 程 学 报 Vol.34 No.20 286 2018年 10月 Transactions of the Chinese Society of Agricultural Engineering Oct. 2018 改进颗粒组构力学模型模拟筒仓卸粮成拱细观机理 冯 永，李 萌 （河南工业大学土木建筑学院，郑州 450001） 摘 要：现有研究表明筒仓卸粮成拱和粮食的内外摩擦密切相关，但现行相关离散元模拟采用单一圆形颗粒，模拟粮食的真实接触面积要小很多，不能客观反映卸粮过程粮食的摩擦情况，也较难还原卸粮成拱现象的细观动态过程。该文针对离散元模拟中圆形颗粒的内摩擦力小于真实粮食内摩擦力的缺陷，在已有PFC离散元程序基础上，添加了黏度系数较大的微型颗粒模拟粉尘，建立了一种改进颗粒组构力学模型，采用几何方法判断圆形颗粒间的接触情况，推导出基本单元间力-位移关系。基于典型事故案例和室内试验成果，采用建立的颗粒组构力学模型模拟了卸粮成拱动态过程中圆形大颗粒间以及圆形小颗粒与仓壁之间的力-位移关系。研究表明在卸粮过程中，切应力在剪切位移达到0.3 mm的过程中，迅速提高，达到最大值切应力的60%，所得的切应力位移图与其应力特征曲线与试验成果基本吻合。以试验结果曲线各点值为标准值，改进后模拟结果曲线值的标准差相比改进前减小37%，说明曲线相似度更高，利用该模型可更加客观反映筒仓卸粮成拱的动态细观机理。该文提出的改进颗粒组构力学模型，不仅可用于模拟卸粮成拱机理模拟，而且对于模拟散颗粒流动特性、散体-仓壁相互作用机理都具有一定借鉴意义。 关键词：模型；离散元；力学特性；卸粮成拱；细观机理 doi：10.11975/j.issn.1002-6819.2018.20.036 中图分类号：TU312+.3 文献标志码：A 文章编号：1002-6819(2018)-20-0286-08 冯 永，李 萌. 改进颗粒组构力学模型模拟筒仓卸粮成拱细观机理J. 农业工程学报，2018，34(20)：286293. doi：10.11975/j.issn.1002-6819.2018.20.036 http:/www.tcsae.org Feng Yong, Li Meng. Simulation of meso-mechanism of silo unloading grain aching based on improved particle composition mechanical model J. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(20): 286293. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2018.20.036 http:/www.tcsae.org 0 引 言 近年来筒仓卸粮成拱相关事故频发，主要起因为大多数筒仓仍采用依靠重力的卸料方法，粮食颗粒在多因素相互作用下的散体物料力学问题下形成拱起，导致应该由仓底承受的压力传递到仓壁，形成破坏。现有研究表明卸粮成拱和粮食的内外摩擦力密切相关，任杰等1-4主要从内摩擦力对贮料固结条件的影响分析了成拱的主要原因。Ayuga等5-10研究了粮食剪胀角及摩擦系数对仓壁侧压力的影响。目前研究多在于宏观机理，随着散体力学的发展，徐泳等11-16均认识到：从细观力学角度看，粮食颗粒之间、粮食与仓壁的摩擦作用（粮食内、外摩擦）是粮食散粒卸料成拱的直接因素。而离散元模拟可以更加客观的反映粮食散体在仓内的细观力学机制，因此针对性开展相关研究重要现实意义。 颗粒流程序（particle plow code，PFC）基于牛顿第二定律和力与位移的关系，李坤蒙等17-24主要用于研究散粒体或可简化为散粒体的系统的分析。目前PFC离散收稿日期：2018-04-10 修订日期：2018-07-01 基金项目：国家自然科学基金：筒仓卸粮成拱及其对仓壁超压作用的动态演进机制研究（51708182）；河南省教育厅基础研究项目：多因素作用下的仓内粮食结拱及塌陷机理研究（16B560002） 作者简介：冯 永，副教授，博士，主要从事粮食力学性质及仓储相关研究工作。Email：51433719qq.com 元模拟主要以单一圆形颗粒单元为基础，而筒仓卸粮成拱主要涉及粮食物理力学特性、筒仓结构特征、散料与仓壁间的相互作用等诸多因素25-27。因此导致离散元模拟粮食的真实接触面积要小很多，不能客观反映粮食的摩擦情况，导致模拟的颗粒摩擦力偏小，模拟卸粮过程无法成拱。 鉴于此，本文基于PFC改进了一种颗粒组构力学模型（将黏度小直径大的圆形颗粒与黏度大直径小的圆形颗粒），编写相应的接触判断以及力-位移关系，对筒仓卸粮过程的细观机理进行动态模拟，深层次揭示其形成机制。 1 改进颗粒组构力学模型 原PFC程序通过生成均匀的二维圆形颗粒试样，用来分析小米大豆的力学特性等功能。改进后的PFC离散元程序的接触模型在基本构架上与Cundall等28-29所提出的模型类似。 改进后程序有筒仓仓壁、圆形大颗粒与圆形小颗粒3种基本单元，其中筒仓仓壁单元用弧面刚性边界表示，粮食颗粒则被模拟化简成为一个具有质量、接触特性和运动惯性的弹性球体，简化接触模型分为法向接触与切向接触，简化接触图如图1所示。 PFC离散元模型数值计算主要参数包括：颗粒正切向刚度比Kn/Ks、颗粒之间接触模量Ec、正向强度Tn与正向强度偏差Tdn、切向强度Ts与切向强度偏差Tds、颗粒第20期 冯 永等：改进颗粒组构力学模型模拟筒仓卸粮成拱细观机理 287 之间与颗粒仓壁之间的摩擦系数u。如图1所示：在颗粒间的本构模型接触中，A代表大型颗粒单元体，B代表小型颗粒单元体，其中Kn，Ks，Tn，Ts，u决定了颗粒之间接触的刚度与强度；弹簧单元连接A、B两颗粒，表示接触点运动前的线弹性行为，运动产生的能量消耗与准静态变形由阻尼器单元模拟产生，通过分离器表示颗粒接触之前作用力为零；颗粒切向接触间的滑片单元遵循Mohr-Coulomb摩擦定律，由此计算确定A，B颗粒之间切应力的数值，切应力应呈现出线弹性增长直到峰值，达到峰值之后颗粒间发生相对滑移，切应力恒定，切应力的数值应由两颗粒之间的黏度系数与摩擦系数共同决定。 注：Kn为颗粒正向刚度，N·m-1；Ks为颗粒切向刚度，N·m-1；Tn为正向强度，N·m-2；Ts为切向强度，N·m-2；u 为摩擦系数；A代表小型颗粒单元体，B代表大型颗粒单元体。 Note: Kn is positive stiffness of the particles, N·m-1; Ks is tangential stiffness of particles, N·m-1; Tn is positive strength, N·m-2; Ts is tangential strength, N·m-2; u is coefficient of friction; A represents the small particle unit; B represents the large particle unit. 图1 改进的颗粒组构力学模型示意图 Fig.1 Schematic diagram of improved mechanical model of particle structure 2 基本单元间的接触判断 2.1 大小圆颗粒间的接触判断 目前离散元模拟软件有2种生成圆形颗粒的方法，不同方法生成的基本颗粒单元之间的接触判断与求解方法不同，分述如下。 一种是由四段圆弧首尾连接形成一个组合体单元，组合体单元可以根据组成圆弧的半径和圆心坐标，与相邻组合体单元利用圆与圆之间的几何接触关系直接进行判断接触情况并求解,计算量大并对硬件条件要求高。 另外一种是软件直接生成圆形颗粒，但是颗粒单元之间的接触情况比较复杂，目前针对其接触情况有两种计算方法。方法1利用两圆的标准方程，直接求解出交点坐标，并将2个交点连线的中点认定为两圆的接触点，这种方法虽然编写程序较为简单但计算量巨大而且由于颗粒较多，精确解很难求出，如若两圆接触面很小还会导致无解。方法2假设相交重叠的两圆为圆A、圆B，把圆A进入圆B区域内最远的点，与圆B进入圆A区域内最远的点连线，将连线的中点作为两圆接触点的几何平均法，但是如果2个最远点的法向不平行，就不能确定两接触圆形颗粒的共用切平面。 基于以上分析，综合考虑编写程序的难易程度、计算机的运算能力和接触点判断精度等因素，本文提出的改进颗粒组构力学模型采用由Lin等30所提出的几何平均法来求解判断圆形颗粒之间接触情况与接触点坐标。 2.2 圆颗粒与仓壁之间的接触判断 改进后离散元模型中圆形颗粒与仓壁之间的接触情况与接触点坐标计算，可以通过坐标变换公式将圆形颗粒所在局部坐标系的方程转化为整体坐标系的方程，然后与代表仓壁的平面直线方程联立求解得到接触点坐标。推导过程如下所示：假设整体坐标系中为半径为Ri的圆形颗粒i的圆心坐标为（xci，yci），横轴与x轴的正向（逆时针方向为正）夹角为i；仓壁单元的起始点坐标为（Xcj，Ycj），倾角为w。当w=90°，仓壁的直线表达式为 x=Xcj （1） 否则，仓壁的直线表达式为 y=tanw(xXcj)+Ycj （2） 局部坐标系中圆形颗粒i的方程为 12222 ii RYRX （3） 通过坐标转换公式，得到X、Y在整体坐标系的方程 ( )cos ( )sin( )sin ( )cosci i ci ici i ci iX x x y yY x x y y （4） 将式（4）整体代入式（3）中，与仓壁对应的方程式（1）联立并求解，若无解，则圆形颗粒与仓壁之间没有接触；若只有一组解，则此解为接触点坐标；若有两组解，则取两点平均值为接触点坐标（xc，yc）。 3 基本单元间的力与位移关系 3.1 大小圆间接触点的力与位移关系 在求解得到大小圆形颗粒之间的接触点坐标后，在一个时步之内，根据接触点位移的变化情况，可以得出此时步内接触力的增量，从而求解出颗粒间接触点的应力数值，然后计算得到这个颗粒所有接触情况的接触力大小，得到此颗粒单元上不平衡的合力与合力距，再使用牛顿第二定律求出该颗粒的线加速度和角加速度，积分后就得到该颗粒在此时步内的位移与速度，推导过程如下所示。 已知圆形大颗粒i与圆形小颗粒j 的圆心坐标分别为（xci，yci），（xcj，ycj），半径分别为RA、RB，横轴与x轴正向（逆时针方向为正）的夹角分别为 i、 j ；对应的x、y 向线速度及角速度分别为vxi、vyi、 i、vxj、vyj、j ；大小颗粒间的接触点坐标为（xc，yc），相交最远点（2个）向量为e1、e2，大小颗粒接触情况如图2所示。 图2中大小颗粒接触点向量 1 22 1 e ei e jt ei e j （5） 式中e1=cos；e2=sin，为相交最远点利用几何平均法求得的接触点与X轴正向的夹角。圆形大颗粒i在接触点处的半径rci由式（6）求得。 2 2( ) ( )ci c ci c cir x x y y （6） 农业工程学报（http:/www.tcsae.org） 2018年 288 注：xci，yci，xcj，ycj为圆心坐标，mm；ir，jr 为接触半径的向量，mm；i ，j 为角速度，rad·s-1；vxi、vyi；vxj、vyj为线速度，m·s-1；e，e1，e2为交点向量，i 、j 为与X轴夹角，rad。 Note: xci, yci, xcj, ycj are center coordinates, m; irand jr are contact radius vector, m; i and j are angular velocity, rad·s-1; vxi, vyi, vxj, vyj are linear velocity, m·s-1; e, e1 and e2 are vector of intersection points; i and j are angle with the X axis, rad. 图2 大小颗粒接触情况 Fig.2 Contact between large round particle and small round particle 接触半径的对应的向量为 ( ) ( ) cos sini c ci c cici i ci ir x x i y y jr i r j （7） 接触点处圆形大颗粒i相对于圆形小颗粒j的速度为 d c ci cjv v v （8） 圆形大颗粒i接触点处的速度利用式（8）计算 ( ) ( sin ) ( cos )ci xi yi i ixi i ci i yi i ci iv v i v j k rv r i v r j （9） 因此，大小颗粒i、j之间接触点处的相对速度为 d d dsin sincos cosc cx cyxi i ci i xj j cj jyi i ci i yj j cj jv v i v jv r v r iv r v r j （10） 将式（10）中相对速度分别向接触点的切向、法向方向投影可得 1 22 1 d d d dd d d dcn n cx cycs s cx cyv v e v e v ev vt v e v e （11） 对相对切向速度分量和相对法向速度分量在一个时步t进行积分，得到相应的位移增量 ddnsn v ts v t （12） 通过力-位移定律计算，得到（12）式中位移增量分别产生力的切、法向力增量 n ns sF k nF k s （13） 式中kn为接触点的法向刚度系数、ks为接触点的切向刚度系数。 在前一个时步上分别将第t 时步计算得到的切向力增量(Fn)t 与法向力增量(Fs)t 加上，便可得到该时刻此触点处的切向力、法向力 n n nt t t ts s st t t tF F FF F F （14） 根据摩尔-库仑定律判断可能最大剪力值 max tans n m mtF F c （15） 式中m为2种圆形颗粒内摩擦角值的较小者，(º)；cm为2种圆形颗粒黏聚力值的较小者，kPa。 若(Fs)t的绝对值大于(Fs)max，则颗粒之间产生相对滑动，取(Fs)t=(Fs)max。将所取的(Fn)t、(Fs)t由局部坐标系转化到整体坐标系中，得出圆形大颗粒i在此接触点所受的作用力分别为 1 22 1sin cosxi n syi n sxi ci i yi ci iF Fe FeF Fe FeM F r F r （16） 将圆形大颗粒i中所有接触力与力矩求和，得到合力Fxi，Fyi 与合力矩M。据牛顿第二定律，计算出此时刻圆形大颗粒i的线加速度和角加速度分别为 1 22 1( )( )( sin cos )i i xi n si i yi n si i xi ci i yi ci imx F Fe Femy F Fe FeI M F r F r （17） 式中mi为圆形大颗粒i的质量，g；Ii为圆形大颗粒i的惯性矩，mm4；ix为圆形大颗粒i在t 时刻x向的加速度，m/s2；iy为圆形大颗粒i在t时刻y向的加速度，m/s2；i为圆形大颗粒i在t时刻角加速度，rad/s2。 因此，圆形大颗粒i在下一个时步中间时刻的速度为 1/2 1/21/2 1/21/2 1/2t t t t tt t t t tt t t t tx x x ty y y tt （18） 利用式（19）可计算出圆形大颗粒i的位移 1/21/21/2t t t t t tt t t t t tt t t t t tx x x ty y y tt （19） 此外，改进后的PFC程序中还可以通过分别或同时更改速度、接触、局部阻尼等物理参数，以及重力等因素对颗粒的力-位移的影响，由于篇幅有限，不再深入探讨。 3.2 圆与仓壁接触点间力与位移关系 圆形颗粒与仓壁的力-位移关系与大小圆形颗粒之间的力-位移关系原理相似，但圆与线之间接触点计算与程序编写相较简单，且2种颗粒与仓壁接触情况仅参数不同，推导过程完全相同。图3为圆形颗粒i与仓壁j的接触情况示意图。 仓壁j的向量如下 1 21 2 2 1 w wi w je ei e j wi wj （20） 圆颗粒i接触点的速度与式（7）相同，仓壁j上接触点的速度可表示为 ( ) ( ) ( )cj xj yj j wxj j c cj yj j c cjv v i v j k rv y y i v x x j （21） 第20期 冯 永等：改进颗粒组构力学模型模拟筒仓卸粮成拱细观机理 289 则圆i相对于仓壁j的接触点处相对速度 d d dsin sincos cosc cx cyxi i ci i xj j cj jyi i ci i yj j cj jv v i v jv r v r iv r v r j （22） 将接触点的相对速度转化为切向与法向方向，接触点法、切向力以及其增量的求解，还有圆形颗粒的运动情况均与3.1小节中大小圆形颗粒之间的求解过程相同。对于仓壁j，在接触点处x、y方向所受的合力与合力矩求解方程分别为1 22 12 1x n sy n sn c cj c cjF Fe FeF Fe FeM F x x e y y e （23） 式中力矩M为接触点相对仓壁的起始坐标(xcj,ycj)求解，接触点的切向力对仓壁的力矩没有影响。 注：wr 为接触半径的向量，mm；w，e 为交点向量。 Note: wr is contact radius vector, mm; w, e are vector of intersection points. 图3 圆形颗粒与仓壁的接触情况 Fig.3 Contact between round particles and silo wall 4 筒仓卸粮成拱细观机理模拟分析 为验证改进后离散元程序的合理性，现取2014年2月黑龙江省兰西县宏鼎粮食购销有限公司生筒仓卸粮成拱事故中筒仓（混凝土筒仓，高12 m，直径5 m）为模拟与试验样本，缩小10倍后，经过几何相似性计算，通过PFC软件分别生成与物理试验所用颗粒大小、数目及力学参数基本相同的颗粒进行模拟卸粮成拱，并沿高度设2列12组共24个监测点，监测卸粮瞬间颗粒间切应力数值。 由于试验与模拟成拱位置不尽相同，模拟成拱后，在拱脚附件由下至上取3个监测点（B1,B2,B3）的图像，和改进前模拟对应位置3个监测点（A1,A2,A3）图像与试验所得拱脚附近3个监测点（C1,C2,C3）图像结果进行趋势对比分析。 4.1 改进前模拟 模型建立：生成筒仓模型尺寸（高×直径）为1 200 mm×500 mm，加入圆形颗粒及模拟过程中的所用参量见表1、2。 表1 改进前颗粒参量 Table 1 Unimproved particle parameters 颗粒种类 Particle types 颗粒A Particle A 颗粒尺寸 Particle size/mm 5 切向接触刚度 Tangential contact stiffness/(N·m-1) 1.5×108 法向接触刚度 Normal contact stiffness/(N·m-1) 1.0×108 密度Density/(kg·m-3) 1.2×103 摩擦系数/（颗粒生产时） Friction coefficient/(in particle production) 1.0×10-10 摩擦系数/（试验过程中） Friction coefficient/(in test process) 0.5 颗粒间的黏聚力 Particle cohesion/kPa 0.5 表2 计算时步参量 Table 2 Calculation step parameters 类目 Category 数值 Numerical value 时步Time step/s 2×10-4 加速度Acceleration/(m·s-2) 0 质量阻尼系数 Mass damping coefficient 刚度阻尼系数 Stiffness and damping coefficient 0 局部阻尼系数 Local damping coefficient 0 图4a为装粮后试样的接触力链图，图中深颜色的粗线条为颗粒间力的传递方向。可以看出，在筒仓装满粮食后，颗粒力传递主要发生在纵向竖直向下。在静止期，大小颗粒各个部位都是受压状态，在上部荷载作用下孔隙率达到最小。 a. 装粮后模型 a. Model after loading grain b. 卸粮后模型 b. Model after unloading grain 图4 改进前筒仓模型装粮与卸粮后的接触力链图 Fig.4 Contact force chain diagram of unimproved silo model after loading grain and after unloading grain 图4b为普通模型卸粮过程接触力链图，图中深色力链基本竖直向下，在颗粒下落过程中力链变形分叉，向仓壁发展，存在局部变形交替发展的现象。即下落开始时出现多条裂缝，且裂缝的发展此消彼长。从微观角度上，由于颗粒之间只存在正应力，缺少黏结力，导致这种局部变形的交替发展，颗粒内部结构应力的不断调整，粮食颗粒间相对位置的变化，导致微裂隙的数量、方位、大小及发展程度等多种因素的不断改变，最终无法形成拱，直至粮食颗粒完全流出筒仓。 农业工程学报（http:/www.tcsae.org） 2018年 290 由图5可知，改进前模型中，在卸粮瞬间，切应力变化幅度较小，证明颗粒之间的黏结力较弱，卸粮过程顺畅，与现实粮仓卸粮过程相差较大，难以出现事故中的成拱现象。 图5 改进前模型中监测点切应力与位移关系曲线图 Fig.5 Relationship between displacement and shear stress of monitoring points in unimproved model 4.2 改进后模拟 模型建立：生成筒仓模型尺寸（高×直径）为1 200 mm×500 mm，圆形颗粒及模拟过程中的所用参量见表1、2，新加入的小型颗粒B参量见表3，大小颗粒A/B数量比为101。 图6a为装粮后试样的接触力链图，与图4a基本相似，此时颗粒静止，在重力作用下，接触力方向主要为纵向向下。颗粒处于三向受压状态，由于存在小颗粒，孔隙率小于图4a，颗粒更为密实。 表3 改进后颗粒参量 Table 3 Improved particle parameters 参数 Parameters 颗粒B Particle B 颗粒尺寸Particle size/mm 1 切向接触刚度 Tangential contact stiffness/(N·m-1) 1.5×108 法向接触刚度 Normal contact stiffness/(N·m-1) 1.0×108 密度Density/(kg·m-3) 1.2×103 摩擦系数/（颗粒生产时） Friction coefficient/(in particle production) 1.0×10-10 摩擦系数/（试验过程中） Friction coefficient/(in test process) 0.5 颗粒间的黏聚力Particle cohesion/kPa 5 a. 装粮后模型 a. Model after loading grain b.卸粮成拱后模型 b. Model after unloading grain arching 注：图6b中黑色粗体线段的方向代表每个颗粒力的传递。 Note: As is shown in Fig.6b, the direction of the black bold line segment represents the transmission of each particle force. 图6 改进后装粮与卸粮后模型的接触力链图 Fig.6 Contact force chain diagram of improved silo model after loading grain and after unloading grain arching 图6b为卸粮成拱后试样的接触力链图，从颗粒的位移，拱的发生发展位置来看，模拟的结果和试验的结果基本吻合。成拱瞬间在自重荷载作用下首先产生压缩,随着荷载的增大，在拱角处的裂隙被压实，在剪切应力作用下，剪切带自外至里，从2个拱脚处开始到完全贯通，粮食整体停止下滑，力链由纵向竖直向下变为横向，应力集中在拱角处，施加在仓壁上。图6b右侧为局部放大的颗粒接触示意图，图6b中黑色粗体线段的方向代表每个颗粒力的传递。从图6b右侧可以更清楚地看到，除了在底部成拱处力横向传递，其余各处颗粒力的传递方向基本没变，保持在纵向竖直方向。 由图7可知，改进后模型中，在卸粮瞬间，切应力迅速提高，很快达到峰值，证明加入黏结力大的小型颗粒之后，颗粒之间的黏结力较强，随着卸粮过程进行，出现了此次事故中的成拱现象。 图7 改进后模型中监测点切应力与位移关系曲线图 Fig.7 Relationship between displacement and shear stress of monitoring points in improved model 4.3 试验验证 由于样本过大，试验采用几何相似原理按比例（101）制作试验模型，筒仓模型高度1 200 mm，直径500 m。试验通过在仓壁布置压力盒来检测筒仓卸粮结拱过程中，颗粒力链的传递与切应力的变化。 试验颗粒为小麦，由小麦与粉尘2种颗粒组成。测量装置为压力盒，在仓壁内侧对称布置2列，底部间隔为5 cm共6个，上部间隔为10 cm，共6个，合计24个。试验过程如图8a所示。在粮食储存运输过程中，粮食颗粒之间摩擦，形成粉尘，经过压实固结，卸粮过程中便出现了图8b中的成拱现象。 a. 试验过程 a. Experimental procedure b. 卸粮成拱试验 b. Unloading arching test 图8 试验过程及装置 Fig.8 Experimental procedure and device 第20期 冯 永等：改进颗粒组构力学模型模拟筒仓卸粮成拱细观机理 291 4.4 数据对比及误差分析 以3号位置探头（A3、B3、C3）为例，对比5、图7、图9：当剪切位移为0.2 mm时，改进前模型中的切应力几乎没有增长，仅为8 kPa，远小于改进后模型的40 kPa，与试验的50 kPa；当剪切位移为1 mm时，改进前模型中的P3探头切应力增长缓慢只达到30 kPa，此时改进后模型的剪切应力已迅速增长达到90 kPa，而试验中切应力也达到了110 kPa。 图9 试验过程中监测点切应力与位移关系曲线图 Fig.9 Relationship between displacement and shear stress of monitoring points in experimental procedure 以3号位置探头试验结果曲线各点值C3i为标准值，利用式（24），求得改进前曲线各点值A3i 的标准差为13.86 kPa，改进后曲线各点值B3i的标准差为8.73 kPa，改进后模型标准差值相比改进前减小37%，拟合效果提高显著。 23 3 3123 3 3111Ni iiNi iiA C ANB C BN （24） 同样对比1、2号位置探头数据发现，改进前模型中的切应力在剪切位移较小时，增长较为缓慢，与试验结果不符。改进后模型中，切应力在剪切位移前期增长迅速，与试验结果基本吻合，可以反映出卸粮成拱过程的切应力与位移关系。 对比图7、图9发现，模拟曲线与试验曲线基本吻合，数值上整体小于试验结果，但差值在10 kPa以内。分析影响粮食内外摩擦力的因素后发现，模拟过程只考虑了粉尘对颗粒的影响，忽略了含水率对粮食内外摩擦力的影响。因此，模拟过程中切应力数值整体略小于试验过程，对比图5、图7说明，相比于含水率，粉尘的缺少对粮食内外摩擦力的影响效果更大。 5 结 论 1）由于PFC的单一圆形颗粒单元不能客观反映粮食的摩擦情况，本文加入小黏度大颗粒与大黏度小颗粒，模拟真实粮食内摩擦角，通过几何方法推导判断颗粒之间、颗粒与仓壁之间的接触情况，建立了一种改进颗粒组构的离散元模型，揭示了卸粮成拱动态过程中圆形大颗粒间以及圆形小颗粒与仓壁之间的力-位移关系，明晰了筒仓卸粮成拱机理。 2）卸粮过程中，当剪切位移为0.2 mm时，改进前模型切应力不足10 kPa，而改进后模型则达到40 kPa；当剪切位移为1 mm时，改进前模型中切应力缓慢增长到30 kPa，此时改进后模型的剪切应力已迅速增长达到90 kPa。试验结果曲线为标准值，改进后模型中的切应力与位移关系曲线图相比改进前标准差减小37%，整体拟合度相比改进前提高显著。 本文提出的改进颗粒组构力学模型，不仅可用于模拟卸粮成拱机理模拟，而且对于模拟散颗粒流动特性、散体-仓壁相互作用机理都具有一定借鉴意义。 参 考 文 献 1 任杰，韩阳. 粮食结拱机理及防范技术研究现状J. 河南科技，2015(13)：102104. Ren Jie, Han Yang. Research status of grain arching mechanism and prevention technologyJ. Journal of Henan Science and Technology, 2015(13): 102104. (in Chinese with English abstract) 2 安蓉蓉. 粮食的内摩擦角、弹性模量及体变模量的试验研究D. 南京：南京财经大学，2010. An Rongrong. Experiment Study on Angle of Internal Friction, Elastic Modulus and Bulk Strain Modulus of GrainD. Nanjing: Nanjing University of Finance and Economics, 2010. (in Chinese with English abstract) 3 彭政，王璐珠，蒋亦民. 颗粒物质与固体交界面静摩擦系数的测量与分析J. 山东大学学报：理学版，2011，46(1)：4245. Peng Zheng, Wang Luzhu, Jiang Yimin. Measurement and analysis of static friction coefficient on a granular-solid interfaceJ. Journal of Shandong University: Natural Science, 2011, 46(1): 4245. (in Chinese with English abstract) 4 丁盛威. 基于PFC3D的筒仓贮料侧压力离散元分析D. 武汉：武汉理工大学，2014. Ding Shengwei. Distinct Element Analysis on Lateral Pressure of Silo Based on PFC3DD. Wuhan: Wuhan University of Technology, 2014.(in Chinese with English abstract) 5 Ayuga F, Guaita M, Aguado P J, et al. Discharge and the eccentricity of the hopper influence on the silo wall pressureJ. Journal of Engineering Mechanics, 2001, 127(10): 1067-1074. 6 杨鸿，杨代恒，赵阳. 钢筒仓散料静态压力的三维有限元模拟J. 浙江大学学报：工学版，2011，45(8)：14231429. Yang Hong, Yang Daiheng, Zhao Yang. Three-dimensional finite element simulation of static granular material pressure for steel silosJ. Journal of Zhejiang University: Engineering Science, 2011, 45(8): 14231429. (in Chinese with English abstract) 7 周长东，郭坤鹏，孟令凯，等. 钢筋混凝土筒仓散料的静力相互作用分析J. 同济大学学报：自然科学版，2015，43(11)：16561661，1669. Zhou Cha