基于滚动摩擦作用下玉米颗粒在模拟筒仓卸料中的力链分析

李坤由，潘 樊，孙慧男

（郑州中粮科研设计院有限公司，河南 郑州 450001）

粮食属于散体颗粒，其摩擦特性要比工程材料之间或固体干摩擦复杂得多，特别是粮食颗粒的不规则性对其在筒仓卸料动态过程中所呈现的摩擦特性有很大影响，而颗粒的摩擦特性又直接或间接作用于颗粒的流动状态，从而影响着粮仓的储、卸安全。

在卸粮过程中，颗粒的流态及应力等散体动力学特征变化对筒仓安全具有重要意义。Weinhart[1]运用离散元法捕捉了 DEM 模拟中的颗粒应力；R. Kobyłka[2]通过模拟储料卸粮过程中颗粒的速度、应力分布探究侧壁压力波动和摩擦特性的演化机制。Yung[3]通过实验和模拟研究发现柔性球在刚性平面上运动时的滚动摩擦与滚动速度呈非线性关系。刘万峰[4]，崔涛[5]，Ketterhagen[6]等设计了不同测定颗粒滚动摩擦的小型实验并通过离散元法进行了验证，为数值模拟提供了基础参数。在颗粒流动现象的研究中，除了筒仓外形和导流体等外部因素，颗粒摩擦特性等内部因素的重要性也逐渐被重视。

Weizman[7]在进行数值模拟时提出将旋转速度转化为平移速度，确定了刚性球在柔性表面上滚动摩擦的变化系数。Silber[8]等研究确定颗粒堆中颗粒空间分布规律依赖颗粒摩擦系数。Goniva[9]发现在单喷口流化床仿真模拟中引入滚动摩擦后能够极大地改善模拟结果。YongZeng[10]将卸粮速度场的波动与仓内力链周而复始的消失现象结合在一起。研究[11-14]将颗粒的宏观现象与细观变化对应起来，然而，即使在单一球模型情况下，筒仓中由于颗粒滚动摩擦耦合效应导致的颗粒力链分布变化和拱效应也没有得到很好的解释或研究。

基于以上分析，对不规则玉米颗粒在筒仓卸料动态过程的摩擦特性的研究具有重要意义。但流态作为可视化的宏观变化现象，在实验中很难对其进行定量分析，因此本文借助离散元建立模型研究卸粮动态过程中的细观变量，搭建宏-细观的桥梁，探究基于滚动摩擦作用下的不规则玉米颗粒在筒仓卸料动态模拟中力链演化机理，以期为筒仓的结构设计和优化提供理论基础。

1 材料与方法

1.1 实验模型

1.1.1 构建模型

1.1.1.1 颗粒模型 本文模拟颗粒原型为我国种植面积较广的德美亚一号玉米品种。样品选取的玉米籽粒的三维几何尺寸如图1所示，上底宽I、下底宽L、高h。在模拟仿真中采用多球法近似玉米颗粒的形状，共计六个球模型相互组合。

图1 玉米颗粒外形尺寸图Fig.1 Dimensions of corn grain

1.1.1.2 筒仓模型 实验室的小型筒仓采用亚克力材料制作，见图2a，得益于其良好的透明性，可以直接观察到筒仓内颗粒在填充和排放过程中的运动。筒仓边长为L=120 mm，仓高H=250 mm；其壁面与水平的夹角为60°，出料口边长D=50 mm。按照实验筒仓尺寸，采用3Dmax绘制同等尺寸的筒仓模型，如图2b所示。

图2 筒仓尺寸图Fig.2 Diagram of silo size

1.1.2 颗粒接触模型

用干燥法测得该品种玉米的含水率为14.23%，由于其水含量较低，可忽略颗粒间的黏附力。假设颗粒堆积运动过程中，其应力、位移、速度是通过颗粒与接触体之间产生的微小交叠量值的不同来确定的。基于此假设，在EDEM软件中建立了与相对速度相关的 Hertz-Mindlin（no slip）接触力学模型来定义颗粒的接触。该模型将颗粒在接触点的碰撞力分解为法向力和切向力，如图 3所示，切向和法向的接触力简化成弹簧和阻尼器的并联。

图3 颗粒接触模型Fig.3 Particle contact model

其中，wi是角速度，rad/s；Tt是切向扭矩，N· m；Tr是滚动摩擦力矩，N· m；Rr是有效滚动接触半径，m；wrel是两颗粒间的相对角速度，rad/s；rμ是滚动摩擦系数。

1.1.3 材料性能参数

玉米颗粒的泊松比按照Tao等[15]描述的结果作为参考值。依据标准ASAE S 368.4(2006)[16]，对单个颗粒进行剪切实验分析，确定玉米的剪切模量，见表1。

表1 材料物理特性参数Table 1 Material physical properties parameters

通过借鉴Gonzalez[11]设计的测定玉米颗粒间以及颗粒与接触材料间静摩擦系数和恢复系数的装置，测得玉米颗粒间及其与亚克力材料间的静摩擦系数，见表2。

表2 接触参数Table 2 Contact parameter

1.2 实验方法

1.2.1 实验过程

以黄色玉米籽粒为主材料，加入少量黑色玉米籽粒做分隔层，以便能够清晰的观察到卸粮过程中颗粒的流动状态。采用中心装料的方式装载颗粒，在装料过程中，由于颗粒的下落具有一定的速度，彼此之间会发生碰撞从而使得不同颜色颗粒发生轻微混合，故颜色界限并不完全水平。

装料完成后，先将满载的模型仓静置2 h，然后打开数码相机，抽走出料口处的挡板使颗粒自然下卸，记录卸粮全过程，实验流程如图4所示。

图4 实验流程图Fig.4 Test flow chart

为减少实验的偶然性误差，进行多次重复实验。

1.2.2 模拟的合理性验证

将 1.1.3中列举的材料性能参数作为基础参数输入到EDEM软件中，进行玉米颗粒的装卸模拟，进行五次平行重复实验，取其平均结果与实验平均结果做对比。

从装料方面看，筒仓内实验装料高度为245.2 mm，模拟装料高度为250 mm；从卸粮时间来看，从抽出底板到颗粒完全流出漏斗，实验的平均时长为6.3 s，模拟的平均时长为6.6 s，模拟误差为4.7%。

对比实验和模拟中颗粒的流动行为（图5），可以观察到，随着卸粮进行，料床高度下降。在0.25T时刻，上部黑色颗粒层较为平直均匀，而下部呈现出明显的V形，颗粒聚拢运动明显。随着时间的增加，筒仓上部黑色颗粒层弯曲明显，在0.75T时，模拟及实验的筒仓与料斗连接处都出现颗粒滞留现象，即出现死料区，整个卸粮过程中实验与模拟的颗粒流动行为基本吻合，验证了模拟的准确性。

图5 实验与模拟流态对比图Fig.5 Comparison of experimental and simulated flow patterns

2 结果与分析

2.1 宏观现象分析

仓壁与颗粒间的摩擦系数保持不变（为非零定值0.031 1），改变颗粒间的滚动摩擦系数μ，发现μ值不但会影响颗粒在筒仓中的流动形态，还对卸粮的时间步长有一定的影响。

增大颗粒间摩擦会在两个方面影响最终时间步长。首先，它会产生较小的体积密度值，从而间接导致T的增加；此外，Balevicius[17]等发现颗粒间的滚动摩擦会影响 Beverroos[18]方程的C常数，从而影响卸料速度，即随着摩擦的增加，卸料速率减小，卸料时间步长增加。

对比三种情况下的颗粒流动状态（图6），可以看出：摩擦系数越小，颗粒的运动越自由，不同颜色颗粒间的混合越剧烈，即在摩擦系数为 0状态下，任意时刻，黑色颗粒“越界现象”更为明显。且摩擦系数越小，黑色颗粒带弯曲得越早，即颗粒的流动更早显现出管状流的特征，在2 s、3.5 s时刻，摩擦系数为0的第一层黑色颗粒带最早开始出现弯曲，而摩擦系数越大，其第一层黑色颗粒带弯曲越不明显，进入 5 s后，颗粒基本都已进入管状流状态。故对于有漏斗的筒仓来说，颗粒间摩擦的减少会改变整体流和管状流之间的极限，从而增加管状流的面积。

图6 三种滚动摩擦系数下颗粒流态对比Fig.6 Comparison of particle flow regimes under three rolling friction coefficients

2.2 颗粒力链的空间分布

2.2.1 切片参数

在筒仓仓体中心处建立 X-Y-Z方向上7-1-10的检测网格，用以检测筒仓内部各空间方位上的颗粒流动状态及其细观参数变化。网格具体划分情况及编号见下图7。

图7 切片点位图Fig.7 Location map of monitoring points

2.2.2 力链的空间分布及概率分析

鉴于筒仓上部颗粒会快速流失且下部靠近漏斗处颗粒速度较大，取Z=3、4、5区域的颗粒接触力进行分析比较。

图8左侧为0时刻颗粒法向力/颗粒平均法向力，比值f在Z=3处水平面的空间分布情况，右侧为f的概率分布情况。

图8 力链分布图（0时刻）Fig.8 Distribution map of force chain(T=0)

颗粒力链在空间分布上呈现出中间力链较强，四周力链弱的特点，形成强力链的颗粒构成颗粒体的主体骨架，可承受较大的荷载，形成弱力链的颗粒作为必要的填充，保证了颗粒体系的整体稳定性。

在右侧f的概率分布上，f ＞ 1的部分为强力链，f ＜ 1的部分为弱力链，f =1附近概率密度较大。这是由于颗粒体内部存在拱效应，在上部颗粒自重下颗粒体内部出现各向异性的力链网络，该力链网络作为荷载的传递路径有强弱之分。对应图中的力链变化趋势，法向力越接近平均法向力，即f越接近1，其分布越广泛，分布概率越大，在法向力与平均法向力相等时其分布概率达到峰值；而在f =1两侧的弱力链和强力链在颗粒体中的分布均较少，使得远离f =1时力链的分布概率较小。

此外所取的上下三个区域的强力链存在互补的现象，其强力链与弱力链的占比也随颗粒所在区域高度的不同发生变化，对应图中 Z=3、4、5时强力链分别占 64%、55%、37%，即颗粒越靠近筒仓底部（Z=3），其强力链占比越高，颗粒挤压力越大，形成拱的概率越高。

改变摩擦系数后，不同高度处的玉米颗粒在不同时刻的力链空间分布及概率分布情况如下图9、10所示。

图9为Z=4的力链分布情况。横向对比不同时间点的力链分布发现，从0.2 T到0.4 T不同摩擦系数下的玉米颗粒的强力连都呈增长状态，说明在Z=4高度处的玉米颗粒随时间增加应力在不断聚集，未观察到成拱塌陷现象。

图9 力链分布图（Z=4）Fig.9 Distribution map of force chain (Z=4)

改变摩擦系数，玉米颗粒在筒仓中的力链空间分布在不同时间点有明显差别。0.2 T时刻，摩擦系数为0.06时颗粒的强力链分布最少，通过概率统计其中 f≥1的强力链占比25%，弱力链占据大部分剖面；随时间增加到0.4 T，应力聚集，强力链占比达到50%，强力链增长速度达到100%。摩擦系数为0时，0.2 T时刻强力链占比为38%，较μ=0.06时增长了52%，说明在这一高度降低摩擦系数玉米颗粒进行应力聚集的时间提前；随着时间的增加，摩擦系数为0的玉米颗粒应力在不断聚集，在0.4 T时达到53%，但其增长速度却小于μ=0.06的情况，应力增长比较平稳，不易出现成拱现象。摩擦系数为0.102时，0.2 T时刻强力链占比为48%，较μ=0.06时增长了92%；随着时间的增加，强力链占比仍在不断增大，在0.4 T时达到78%，远高于其他两种摩擦情况，但其增长速率上仍小于μ=0.06，说明在这一高度增大滚动摩擦系数玉米颗粒进行应力聚集的现象明显，但应力增长速率相对并不高。

图10为Z=3高度的力链分布情况，该高度位于仓身靠下 1/3处，由于更接近漏斗口，在前期卸料时流动截面收缩，颗粒迅速集聚，导致颗粒间应力增加，即在0.2 T时刻Z=3高度处的强力链占比都比同条件下Z=4高度的占比高。其中以μ=0.06条件下的力链变化最为显著，强力链由25%增长至 92%，颗粒迅速集聚，应力激增，随着时间的增加，强力链占比减小，出现应力消散，形成了起拱-塌陷效应。削弱颗粒间滚动摩擦μ=0，颗粒的强力链占比由Z=4高度的38%到Z=3高度的 49%，应力有所增长，但其增长幅度远小于μ=0.06条件下的应力变化；随着卸料的进行，强力链仍在小幅度增长，未出现应力消散，说明减小滚动摩擦玉米颗粒在Z=3到Z=4高度保持较稳定的畅通卸料，未出现起拱的应力激增，以及拱塌陷的应力衰减。这一应力变化现象与该摩擦条件下卸料时长最短相吻合。增大颗粒间滚动摩擦μ=0.102，在0.2 T时刻，强力链的占比随其下落高度的变化基本保持不变，随着卸料的进行，0.4 T时刻强力链占比为31%，相比0.2 T时刻的49%减少了37%，出现应力消散现象，但与其他两种摩擦条件不同的是，该时刻Z=3高度的强力链相比Z=4高度的反而有所减小，由31%减少至28%，说明增大颗粒间滚动摩擦不但会增加供效应，且堆积成拱高度距离漏斗口更高。

图10 力链分布图（Z=3）Fig.10 Distribution map of force chain (Z=3)

3 结论

本文利用EDEM建立了筒仓模型和玉米仿真颗粒模型，模拟装卸粮实验，并以物理卸粮实验验证了模型及模拟的准确性。在此基础上对模拟数据进行分析，探究了不同滚动摩擦系数条件下，玉米颗粒卸料时的力链分布规律，分析了筒仓内起拱-塌陷现象的应力变化。模拟得出以下结论：

（1）改变颗粒间的滚动摩擦系数，会影响颗粒在筒仓中的流动形态。颗粒间滚动摩擦系数越小，颗粒的流动由整体流转变为管状流的时间越早。对于有漏斗的筒仓来说，颗粒间摩擦的减少会改变整体流和管状流之间的极限，从而增加产生管状流的面积。

（2）0 T时刻筒仓中下部三个区域的强力链存在互补的现象，其强力链与弱力链的占比也随颗粒所在区域高度的不同发生变化，即颗粒越靠近筒仓底部，其强力链占比越高，颗粒挤压力越大，形成拱的概率越高。

（3）卸料初期，摩擦系数为0.06时颗粒的强力链分布最少，但随时间增加，玉米颗粒迅速集聚，应力增幅最大，后期，随着颗粒高度的减小，出现应力消散，呈现起拱-塌陷效应；降低摩擦系数，玉米颗粒进行应力聚集的时间提前，应力增长比较平稳，未出现起拱的应力突增，以及拱塌陷的应力衰减；增大滚动摩擦系数，玉米颗粒进行应力集聚的现象明显，滚动摩擦不但会增加供效应，且出现成拱高度距离漏斗口更高。