颗粒堆底面压力凹陷的离散元模拟分析

刘媛媛，陈 可，郭一辰，戴北冰

(1.华北科技学院，北京 东燕郊 065201；2. 太原理工大学，山西 太原 030024；3.中山大学，广东 广州 510275)

0 引言

颗粒材料在自然界中普遍存在，比如砂石岩土材料、谷物颗粒等等。颗粒材料集合体可以呈现不同于液体及固体材料的特性，比如颗粒堆底面压力的最大值并非对应于颗粒堆的形状顶点，而是处于形状顶点的周围的环状区域，即“应力凹陷”现象。对“应力凹陷”现象的试验探索和理论解释具有重要的基础理论认识意义，也一直是颗粒材料相关的研究热点之一[1-6]。离散元法(Discrete Element Method)是20世纪70年代发展起来的用于计算散体介质系统力学行为的一种数值方法，可以合理科学的模拟颗粒介质的微观力学特征，近年来成为研究颗粒体系相关问题的常用方法之一。

本文利用离散元手段模拟颗粒堆的形成过程，对形成的模拟颗粒堆的底部压力分布与实际物理模型试验结果进行对比，从细观力链分布及其它细观指标的分析来进一步解释压力凹陷现象底层的原因，有助于更全面深入的理解认识颗粒体系的基本力学特性。

1 离散元分析方法介绍及模型建立

1.1 特点

离散元方法基本假定是块体单元为理想刚体，各块体的运动只是空间位置的平移和绕形心的转动，其自身的形状和大小保持不变；块体之间的接触视为角-边接触或边-边接触；块体之间的接触作用力由节理面的刚度、接触点的相对位移及有关的阻尼力确定。离散元方法基本原理是集合中每一个单元都是独立的，每个单元都具有相应的尺寸、质量、转动惯量和接触参数等属性。它以牛顿第二定律和力-位移定理为基础，对每一个单元首先确定与之接触的单元，根据单元之间的重叠量，运用力-位移定理计算单元之间的接触力，从而得到单元的合力和合力矩，之后用牛顿第二定律确定单元的运动规律，如此循环计算，直到系统中所有颗粒都计算完毕。颗粒流(Particle Flow Code)是目前应用比较广泛的离散元商用软件之一。本次研究采用PFC2D进行颗粒堆试验模拟分析。在PFC计算中不需要给材料定理宏观本构关系和对应的参数，这些传统的力学特性和参数通过程序自动获得，而定义它们的是颗粒和水泥的几何和力学参数，如颗粒级配、刚度、摩擦力、粘结介质强度等细观力学参数。颗粒流PFC2D有不同的接触模型，分别用于模拟不同的接触材料。一般常用的有线性接触模型、线性接触黏结模型、平行黏结模型等。本文研究模拟对象为砂土，砂土之间的接触接近于完全刚性变形，故选取线性模型[7]。

1.2 参数选取

在整个试验中，通过赋予细观参数来影响整个宏观的模拟合理性，模型中的各参数指标及取值见表1。所模拟的对象为砂土材料，取颗粒半径为0.6～0.75 cm范围内的随机分布[7]。参考相关研究工作[7,8]，颗粒间滑动摩擦系数取值0.977，颗粒和墙体之间的滑动摩擦取值0.7，颗粒间切向刚度和法向刚度取值均为1×107N/m，颗粒与墙之间切向刚度和法向刚度取值均为1×108N/m。通过姚玉相等[7]和周国庆等[9]在颗粒模拟中分析，如果仅使用滑动摩擦很难使系统在合理的分析步内达到平衡，故加入阻尼比为0.3来使系统耗能更快，减小迭代次数。

表1 模型细观参数

1.3 模型建立

本次模型的建立主要分为以下步骤：①以墙构件建立下料漏斗边界，漏斗下落口暂时封闭。边界多边形宽为W为50 cm，底部宽W0为15 cm，倾斜角度θ为30°，侧边长度H为40 cm，下料口距离底板高度H0为40 cm；②在漏斗边界区域上部，在0.6～0.75 cm的半径范围内随机生成1176个颗粒，施加重力加速度9.81 m/s2，颗粒落入漏斗内并达到重力场下的初始平衡；③删除漏斗下方边界墙，颗粒由漏斗下部出口自由落体；④模拟固定中心点源方式，在底板墙上堆积形成颗粒堆并达到平衡稳定状态。

图1 模型建立示意图

2 模型计算结果分析及讨论

2.1 颗粒堆形态

通过中心点源法堆积并达到最终体系内部静力平衡后的颗粒堆形态如图2(a)，外部形态与本文作者开展的室内物理模型试验[10]图2(b)较为一致。其中，如图2(a)所示，模拟颗粒堆下部的底板用Facet墙结构模拟；在装料漏斗中由低到高的六个不同高度范围内，生成相应颗粒，分别命名为LevelOne到LevelSix，用不同的颜色表示。可以看出，底层区域生成的颗粒，也是最先自由落体至底面墙的颗粒，多数处在颗粒堆底部的边缘部位。生成位置较高，较迟落至底面的颗粒则多数位于颗粒堆的中部，甚至占据了底面中心部位。可见，由中心固定点源法堆积形成的颗粒堆，后续降落的颗粒会将前期已落定的颗粒推向两侧，势必会导致中心颗粒会受到两侧颗粒的指向中心方向的水平推力，而两侧颗粒会受到底面施加的指向中心方向的水平摩檫力。

2.2 内部力链

颗粒堆内部力链分布如图3所示，颗粒之间及颗粒与墙间的接触力的大小以不同的颜色表示。颗粒堆在自重荷载条件下，总体来说，浅表层区域颗粒间接触力普遍较小，而内部接触力较大。但颗粒堆内部接触力分布不均匀，其中最大值出现在堆体中心的两侧区域。力链呈网状分布，并明显存在由中心向两侧倾向的主力链，构架出类似于拱形构造力链形态。

图2 颗粒堆形态

图3 颗粒堆内部接触力链分布图

2.3 底面压力分布

提取颗粒与底部墙所有接触力的数值及相应位置，得到颗粒堆底部法向压力在径向上的分布，如图4。为了得到接触压力的相对大小及其分布的相对位置，并和室内物理模型试验进行对比，将接触力及位置分别除以γh和堆体半径进行标准化处理，得到标准化后的颗粒堆底面压力沿径向的分布如图5。由于散粒体系内接触力离散性较大，为了了解更有价值的数据分布趋势，对模拟结果数据进行相邻数据的移动平均得到趋势线，如图5中红色虚线。模拟数据趋势线与室内物理模型试验数据进行对比，吻合性较好，呈现出颗粒堆中心压力凹陷特征。

2.4 细观结构特征

细观结构的特征描述，可采用组构张量来定量分析颗粒间接触方向、颗粒朝向等细观指标的分布情况[11,12]。由于本次模拟是采用一定尺寸范围内随机分布的二维圆盘来模拟颗粒材料，并且研究焦点为颗粒堆内部粒间力的传递方向，考虑到颗粒堆为左右对称的形态，故仅对颗粒堆左侧一半的区域内颗粒接触力法向方向分布进行统计分析。根据Rothenburg等[13]提出的用来近似描述以上细观指标离散分布及各向异性的连续函数：

图4 颗粒堆底部接触压力径向分布图

图5 颗粒堆底部压力标准化值径向分布图

Ec(φ)=E0[1+ancos2(φ-φn)]

(1)

式中，Ec(φ)为颗粒接触法向的分布函数。其中an反映分布的各向异性程度，而φn对应表示各向异性的主方向，E0是各向同性状态时的分布函数，取值为常数1/2π。对模拟试验的颗粒堆体左侧一半的区域内粒间接触法向数据进行统计分析，并与上文中的连续函数对比，可反推得到粒间接触法向这一指标的各向异性特征：an为0.008；φn为30°。在模拟颗粒堆内部粒间法向力方向的分布具有一定程度的主方向的定向性，主方向大约偏离竖向30°。这种分布定向性说明颗粒堆内部对重力荷载的传递会偏离通常认为的沿重力的竖向传递，从而可一定程度上解释底部压力的分布并非与颗粒堆剖面形状一致，而产生堆顶投影底面部位的压力凹陷现象。

2.5 讨论

从颗粒堆体内部力链形态和细观结构特征分析，可反映出由中心点源方式形成的颗粒堆，在重力作用下的力的传递具有偏离竖向的主力链定向特征。这种特征与表示颗粒堆内力传递方向的Edwards拱图[8]假定不谋而合。Edwards拱图是较早的理论假定模型，如图6所示。在此模型中，进行了定性描述：当以中心点源方式堆积颗粒时，颗粒所受重力假定为以一系列的层叠的拱式结构沿斜向直线传递，会形成底面压力最大值出现在堆体顶点两侧的位置，而中心点压力为零。很显然，这个Edwards拱图假定夸大了压力凹陷的程度，但可以从一定程度上解释了颗粒堆体底面中心压力凹陷现象的可能原因，得到了众多后续研究的借鉴。

图6 颗粒堆内力传递的Edwards拱图及对应的底部压力分布图[8]

3 结论

(1) 利用离散元手段模拟颗粒堆的中心点源方式堆积过程，分析堆体底部压力和堆体内部力链的分布特征，得到与室内物理模型试验较吻合的堆体底部中心部位压力凹陷的分布形态。

(2) 基于堆体内部主力链分布形态和粒间接触力的法向角度各向异性特征分析，结合Edwards拱图模型，解释了颗粒堆体底部中心压力凹陷的原因来自颗粒集合体内部力的斜向传递特征，是颗粒堆体中土拱效应的体现。