基于离散元方法的粮食自动分级研究

徐登峰　成荣　朱煜　等

摘要：以浅圆仓中垂直落料的大豆及杂质为研究对象，以软球假设为基础建立颗粒的接触力学模型，采用离散元方法对自动分级过程进行分析，研究颗粒形状、计算时间步长对结果的影响，结果表明，最终分析结果与实际扦样试验数据一致，离散元方法应用于粮食自动分级研究有较好的有效性，具有重要的实际意义。

关键词：自动分级；离散元；接触力学；瑞利波

中图分类号： S226.5文献标志码： A文章编号：1002-1302（2015）01-0396-03

收稿日期：2014-01-14

作者简介：徐登峰（1974—），男，四川华银人，博士，助理研究员，主要从事超精密机械结构及动力学研究。Tel：（010） 62780925；E-mail：xudf@u-precision.com。

通信作者：成荣，博士，在清华大学材料学院从事博士后研究工作。Tel：（010）62781457；E-mail：chengr@tsinghua.edu.cn。 浅圆仓由于容量大、占地面积小等优点，在粮食储存中得到广泛应用。由于粮食中会掺杂很多杂质，在进仓落料及堆积过程中会产生杂质聚集现象，即自动分级现象[1-2]，这会导致粮食孔隙密度减小、带菌量增大、易发生霉变板结等，长时间储粮会对浅圆仓内粮食造成大量损失。因此，对浅圆仓粮食自动分级现象进行研究，避免分级现象的产生具有重要的实际意义[3-5]。目前，针对粮食自动分级现象的研究主要通过试验手段[6-9]，如周延智等研究浅圆仓的杂质分布规律、张峻岭等研究布粮器在浅圆仓防分级中的应用等[10-11]，这些研究方法会耗费大量的人力和物力资源，且研发周期较长。为降低成本和缩短研发周期，本试验采用离散元方法对粮食自动分级现象进行分析研究。

1离散元方法及接触力学模型

1.1离散元方法

离散元方法（discrete element method，DEM）是建立在牛顿第二运动定律基础上研究非连续性颗粒物质结构和运动规律的一种数值方法，根据颗粒间重叠量计算接触力，依次更新每个颗粒的速度和位置，进而确定性地演化整个颗粒系统[12]，其具体求解过程采用显式解法，先利用中心差分法将运动方程进行离散，然后在每个时间步长内进行一次迭代，根据前一次迭代所得到的颗粒位置，由物理方程求出接触力，作为下一次迭代的出发点，再用来求出颗粒的新位置，如此反复进行迭代，直至最后达到稳定流动过程为止[11]。在大豆入仓过程中，大豆、杂质和粮仓相互之间会发生接触碰撞，采用离散元方法进行计算，可准确有效地得到大豆及杂质在任意时刻的速度和位置，最终获得粮仓内的实时颗粒分布状况，从而研究粮食自动分级现象。

1.2颗粒接触力学模型

颗粒接触理论严密，求解过程相对繁琐，在不产生显著误差的条件下，有必要予以适当简化处理。软球模型是目前常用的简化模型，把颗粒间法向力用弹簧和阻尼器模型描述，切向力用弹簧、阻尼器和滑动器描述，并引入弹性系数和阻尼系数等参量[13-14]，依据颗粒间法向重叠位移和切向位移计算接触力，不考虑颗粒表面变形和接触力加载历史，计算强度较小，适合于工程问题的数值计算[15-16]。

以软球模型假设为基础的颗粒之间法向及切向力力学模型如图1所示[17]，当颗粒之间的切向作用力小于或等于最大静摩擦力时，颗粒之间受到法向作用力和切向作用力；当颗粒之间的切向作用力大于最大静摩擦力，颗粒之间会发生相对滑动，从而受到法向作用力及滑动摩擦力。

通过作用力模型计算可得，颗粒之间法向作用力表达式为：

Fn=Knδ3/2n+Dnδ1/4nδn。（1）

式中：Fn表示颗粒之间法向作用力；δn表示颗粒之间法向重叠位移；Dn表示颗粒之间法向阻尼系数。

根据Hertz接触理论，确定法向刚度系数Kn表达式[18]为：

Kn=431-υ2iEi+1-υ2jEj-1ai+ajaiaj-1/2。（2）

式中：υi表示颗粒i的泊松比；Ei表示颗粒i的弹性模量；υj表示颗粒j的泊松比；Ej表示颗粒j的弹性模量；αi表示颗粒i的球体半径；αj表示颗粒j的球体半径。

颗粒之间切向作用力为静摩擦和动摩擦值之间的最小值，表达式为：

Ft=min（Ktδt+Dtδ1/4tδ′t，μFn）。（3）

式中：Ft表示颗粒之间切向作用力；δt表示颗粒之间切向相对位移；Dt表示颗粒之间切向阻尼力。

根据Mindlin-Dereiewicz接触理论，确定切向刚度系数Kt表达式[15]为：

Kt=8δ1/2n1-υ2iGi+1-υ2jGj-1ai+ajaiaj-1/2。（4）

式中：Gi表示颗粒i的剪切弹性模量；Gj表示颗粒j的剪切弹性模量。

2离散元分析

2.1颗粒建模方式选择

根据大豆和杂质实物（图2）实际取样统计发现，大豆主要为球形，杂质以四面体和球形为主。

以多个球形模型拟合四面体模型（图3），分别以四面体模型和球形模型进行建模，在同等条件下进行验证性计算，并将计算结果中的杂质分布状况绘制云图（图4、图5）。由图4、图5可见，杂质采用四面体建模和采用球形建模，对杂质的分布状况影响不大，采用球形模型计算可以大大提高计算速度[15]。因此，离散元分析中颗粒均采用球形建模方式。

2.2分析模型的建立

考虑粮食自动分级的实际因素，在分析过程不考虑空气阻力影响。按照实际粮仓尺寸进行建模，颗粒入仓方式为顶部垂直落料。根据实际大豆与杂质的质量百分比设定大豆与颗粒数目，大豆、杂质与粮仓模型的物理性质参数见表1。

表1大豆、杂质和粮仓的物理性质参数endprint

颗粒密度

（kg/m3）剪切模量

（Pa）泊松比大豆1 2681.5×1080.4杂质2 0003.1×1070.3粮仓4 0005.0×1090.33

2.3时间步长计算

2.3.1瑞利时间步长颗粒发生接触碰撞时，总能耗的70%通过瑞利波（Rayleigh wave）而消耗。因此，可根据沿固体球形颗粒表面传播的瑞利波波速来确定临界时间步长[18-19]。不同颗粒组成的系统，瑞利时间步长计算公式为：

Δt=πR0.163υ+0.877ρΔGmin。（5）

式中：R表示颗粒半径；υR表示瑞利波速；υ表示泊松比；G表示剪切模量；ρ表示密度。

两颗粒间的接触作用仅限于发生碰撞的两颗粒上，而不应该通过瑞利波而传递到其他颗粒上，因此，时间步长应小于瑞利波传递半球所需要的时间。在实际计算时，依据颗粒运动剧烈程度选取合适的时间步长，以保证颗粒系统演变的计算稳定性[12，16]。

对浅圆仓内实际物料分布进行扦样试验，得到杂质分布图（图6）。由图6可见，杂质主要分布在仓中心约15.0 m高的圆台形区域内，超过15.0 m高度杂质数量减少，堆积在一个底面角约为40°、近似圆锥的区域内。

2.3.2计算时间步长确定设定圆台上的表面直径D为杂质分布直径， 以杂质分布直径D和计算所需时间t为评价指

标来选择合适的计算时间步长Δts。当Δts设定为瑞利时间步长Δt的40%以上时，无法保证计算稳定性，过小，又会导致模型计算时间过长[11，20]。为确定合适的计算时间步长，令：

Δts=αΔt，α∈（0.05，0.4）。（6）

在此范围内，采用Latin HyperCube方法设计采样点进行模型计算（表2）。采用Radial Basis Function（RBF）拟合方法[21]，拟合计算时间和计算精度（杂质分布直径）随计算时间步长的变化关系。由图7、图8可见，计算所耗时间和计算时间步长接近线性关系，当计算时间小于瑞利时间的175%，计算精度（杂质分布直径）趋于收敛。综合考虑模型的计算精度及计算效率，设定计算时间步长为瑞利时间步长的17.5%。

表2不同计算时间步长下的计算效率和精度

α（%）t（h）D（m）37.242410.7031.543311.5727.793512.1822.683712.2122.273812.2317.344212.2613.285112.276.455912.28

3结果与分析

利用得到的所有参数，进行颗粒堆积离散元计算，对计算结果进行数据分析：将粮仓进行分区，提取每个分区内的杂质

及颗粒数目，计算得到杂质数目所占比例；采用拉格朗日二次插值方法对数据进行处理，得到杂质在粮仓内的分布云图（图9）。由图9可见，杂质分布主要集中在粮仓的中心高约14.4 m的圆台形区域内，超过14.4 m高度杂质数量减少，堆积在一个近似圆锥的区域内，圆锥底面直径约12.26 m，底面角约41°，这与粮仓实际物料分布扦样试验的结果（图6）基本一致。采用离散元方法计算，能够准确有效地进行粮食自动分级现象研究，该方法可以为解决粮食自动分级现象提供有效的分析手段。

4结论

以浅圆仓中垂直落料的大豆及杂质为研究对象，建立以软球模型假设为基础的颗粒间法向及切向力切换力学模型，推导出作用力表达公式。在此基础上，基于离散元方法对自动分级过程进行分析计算，详细论述颗粒建模方法对计算速度及分析结果的影响，对计算精度影响较大的计算时间步长，引入Latin HyperCube分析方法，得到计算时间及计算精度随步长的变化关系，确定计算时间步长值，且最终计算结果与扦样试验数据一致。离散元方法在粮食自动分级研究工作中具有较好的有效性。

参考文献：

[1]杨文生，张中，张成. 浅圆仓储粮自动分级研究现状[J]. 粮食加工，2011，36（5）：74-76.

[2]吕新. 浅圆仓储粮中心杂质集中问题的解决方案[J]. 粮食流通技术，2003（1）：19-20.

[3]王永昌. 浅圆仓粮食自动分级问题的讨论[J]. 粮食流通及仓储，2002（12）：11-13.

[4]李洪昌，李耀明，唐忠，等. 基于EDEM的振动筛分数值模拟与分析[J]. 农业工程学报，2011，27（5）：117-121.

[5]张来林，张爱强，朱彦，等. 粮食自动分级的类型与预防措施[J]. 粮食储藏，2008（6）：28-30.

[6]姜永顺. 几种浅圆仓进仓设备的比较[J]. 粮食与饲料工业，2002（9）：16-17.

[7]张峻岭，陈艺，张卓青. 布粮器在浅圆仓入粮防分级中的研究与应用[J]. 粮食流通技术，2012（6）：16-22.

[8]王永昌，姚文冠，李军五. 浅圆仓粮食入仓自动分级原因分析及解决措施[J]. 粮食流通技术，2010（3）：28-29.

[9]程绪铎，黄之斌，石翠霞，等. 筒仓中入粮自动分级的动力学分析及缓解措施研究[J]. 粮食储藏，2012（1）：16-19.

[10]周延智，李松伟，曾卓. 浅圆仓杂质分布规律的研究[J]. 粮油仓储科技通讯，2009（4）：14-16.

[11]张峻岭，陈艺，张卓青. 布粮器在浅圆仓入粮防分级中的研究与应用[J]. 粮食流通技术，2012（6）：16-22.

[12]Richarda P，Taberlet N. Recent advances in DEM simula-dons of grains in a rotating drum[J]. Soft Matter，2008（4）：345-1348.

[13]焦玉勇，葛修润，谷先荣. 三维离散元法中的数据结构[J]. 岩土力学，1998，19（2）：74-79.

[14]焦玉勇，葛修润. 基于静态松弛法求解的三维离散单元法[J]. 岩石力学与工程学报，2000，19（4）：453-458.

[15]孙其诚，王光谦. 颗粒物质力学导论[M]. 北京：科学出版社，2009：41-42.

[16]Teufelsbauer H，Wang Y，Chiou M C，et al. Flow-obstacle interaction in rapid granular avalanches：DEM simulation and comparison with experiment[J]. Granul Matter，2009，11（4）：209-220.

[17]Grima A P，Wypych P W. Investigation into calibration of discrete element model parameters of scale-up and validation of particle-structure interactions under impact conditions[J]. Power Technology，2011，212（1）：198-209.

[18]吴清松，胡茂彬. 颗粒流的动力学模型和实验研究进展[J]. 力学进展，2002，32（2）：251-260.

[19]李婉宜，曾攀，雷丽萍，等. 离散颗粒流动堆积行为离散元模拟及实验研究[J]. 力学与实践，2012，34（1）：20-25.

[20]徐泳，黄文彬. 颗粒离散元法建模和仿真的若干进展[J]. 过程工程学报，2002，2（增刊）：530-536.

[21]杨洋，唐寿高. 颗粒流的离散元法模拟及其进展[J]. 中国粉体技术，2006，12（5）：38-43. 丁克奎，钟凯文. 基于“3S”的精准农业管理系统设计与实现[J]. 江苏农业科学，2015，43（1）：399-401.endprint