颗粒肥料离散元仿真边界参数系统化研究

2018-09-17刘彩玲宋建农黎艳妮张福印

农业机械学报 2018年9期

刘彩玲魏丹宋建农黎艳妮都鑫张福印

(中国农业大学农业部土壤-机器-植物系统技术重点实验室，北京 100083)

0 引言

离散元法(Discrete element method，DEM)是一种基于分子动力学的颗粒物料分析方法，近年来，随着计算机技术的发展，该方法在农业装备研究上应用越来越多[1-3]。颗粒肥料是一种容易施用的优质高效肥料，便于机械化作业，在施肥机械上有着广泛的应用，因此，全面系统研究颗粒肥料离散元仿真参数，有助于离散元法在施肥装备研究中的应用[4]。

在离散元颗粒系统模拟过程中，颗粒物料固有的物理力学参数是影响颗粒运动特性的重要因素，正确的模型参数对模拟结果的正确性具有重要的意义[5-6]。NAKASHIMA等[7]通过模拟重力对沙堆颗粒休止角的影响指出，沙堆堆积角受重力影响很小，主要与颗粒半径、摩擦因数等因素有关。ZHOU等[8]模拟了玻璃圆球颗粒在长方体容器中落料堆积过程，发现颗粒间及颗粒与容器壁滑动和滚动摩擦因数是影响堆积形态的关键因素。PERSSON等[9]在评价离散元预测颗粒系统流动性时观测到颗粒间不同动静摩擦因数值直接影响落料速率及堆积角变化。GONIVA等[10]采取颗粒堆积角为模拟指标发现，采用具有滚动阻力模型的离散元方法对模拟结果有很大改善。韩燕龙等[6]通过对稻谷休止角的仿真，对稻谷颗粒间滚动摩擦因数进行了研究；崔涛等[11]基于高速摄像技术研究了玉米种子滚动摩擦特性，并通过玉米种子休止角试验对EDEM边界参数进行标定；王云霞等[12]研究了3种不同形态玉米模型的休止角仿真试验，确定了休止角与颗粒间摩擦特性的相关模型。以上国内外众多学者的研究表明，离散元模型边界参数对仿真精确性具有十分重要的影响，而且通过离散单元法研究散体颗粒堆积过程能有效模拟物料物理力学参数、颗粒接触模型变化时对散体运动特性的影响情况。

离散元仿真模型的边界参数包括材料的基础物理参数(如泊松比、剪切模量、实体密度等)和材料间接触参数(如碰撞恢复系数、静摩擦因数以及滚动摩擦因数)。目前，仿真边界参数主要取自试验测定或经验值。LIU等[13]基于弹跳试验进行了水稻恢复系数的测量；BALEVIIUS等[14]通过滑动试验测量了豌豆与有机玻璃间的静摩擦因数；GONZLEZ-MONTELLANO等[15]对离散元仿真中部分参数的直接测量方法进行研究，结果表明，对于玉米及橄榄核这类外形不规则的颗粒材料，直接测量误差很大。对于颗粒间接触特性参数，韩燕龙等[6]将稻谷粘结在倾斜板面上制成颗粒板，测定稻谷间静摩擦因数；崔涛等[11]将玉米种子粘贴在硬纸板上测定玉米种子间滚动摩擦因数，采用颗粒板法无法保证完全一致的颗粒层厚度且表面也有凹凸，种子运动过程中的弹跳、碰撞无法避免，因此直接测定比较困难，误差较大。

综上所述，离散元仿真模型边界参数是影响颗粒运动的重要因素，边界参数的确定对仿真的准确性具有重要的意义，但目前边界仿真参数主要取自经验值，大部分参数直接测定比较困难，特别是对于材料间接触参数的测定误差较大，边界参数的研究还缺少系统性。颗粒堆积表征了颗粒物料流动、摩擦等特性，是研究散体力学问题的重要内容[16]，同时离散元法是对颗粒随机堆积过程进行计算模拟，认识堆积的细观力学机理和评估所采用模型适用性的有效方法[17]，基于此，本文以反映散体颗粒群堆积过程宏观特征的休止角为宏观现象，以尿素颗粒为研究对象，利用离散元仿真软件EDEM系统研究尿素颗粒仿真边界参数对休止角的影响，对重要边界参数进行分析与确定，以提高离散元仿真精度，缩短设计周期。

1 休止角测定方法

采用ABS无底圆柱桶法测定肥料颗粒休止角，休止角仿真设置如图1所示。参考文献[17]确定圆筒直径44 mm、高135 mm，提升速度为0.05 m/s。研究中对尿素颗粒进行6组试验，每组测定均从垂直交叉的4个方向拍摄图像并求其平均值，最后将6组平均值作为最终测定的休止角。

图1 休止角仿真设置Fig.1 Rest angle simulation settings1.ABS板 2.ABS无底圆筒 3.肥料堆

为减小测量休止角误差，利用Matlab读取颗粒堆积图像(图2)，依此对图像进行灰度处理、二值化处理，最后提取边界点，边界点上的连线即为颗粒堆积的边界曲线，利用最小二乘法对边界点进行直线拟合，拟合的直线斜率即为所要测得的休止角的正切值。

图2 尿素颗粒堆积角测量Fig.2 Measurement of stacking angle of urea particles

2 颗粒及接触力学模型

2.1 试验对象及颗粒仿真模型的建立

试验选取干燥的大颗粒尿素为研究对象，外观形态接近球体，采用精度为0.01 mm的千分尺测得其算术平均粒径为3.35 mm，大颗粒尿素直径分布范围为1.7～4.75 mm，千粒质量23.762 g，体积密度为1 330 kg/m3。

以3.35 mm作为颗粒模型直径，在EDEM中构建如图3a所示的颗粒离散元模型。根据大颗粒尿素直径分布范围，设置颗粒工厂颗粒尺寸按正态分布方式生成，将参数Mean设置为1，即正态分布的颗粒平均直径为尿素颗粒模型直径，标准差 Std Dev 设置为0.1，Capped用来设置生成颗粒的尺寸极限，最小为平均值的0.5倍，最大为平均值的1.4倍，尺寸极限与筛分法确定的尿素颗粒尺寸分布一致。具体设置如图3b所示。

Scale By指颗粒尺寸缩放标准为颗粒半径。

图3 尿素颗粒模型以及生成方式Fig.3 Urea particle model and its generation method

时间步长与颗粒粒径、运动速度、密度、剪切模量有关，通常设定为Rayleigh步长的5%～30%，本模型中颗粒排列紧密，确定时间步长为Rayleigh步长的25%。根据圆筒上升速度并考虑颗粒达到稳定状态的需要，设定仿真总时长为15 s。

2.2 颗粒堆积的接触力学模型

假设肥料表面没有粘附力，颗粒运动特性取决于其物体或相邻颗粒间接触点位置会发生的微小弹性形变，因此颗粒堆积过程选择Hertz-Mindlin无滑动接触力学模型(图4)，接触点法向和切向碰撞力均简化为弹簧和阻尼器的并联。

图4 接触力学模型Fig.4 Model of contact force

基于上述模型，颗粒在堆积过程中主要受自身重力、颗粒间法向力Fn、切向力Ft、法向阻尼力Fnd和切向阻尼力Ftd，其中切向力与库伦摩擦力μsFn(μs为静摩擦因数)有关，各力大小为

(1)

Ft=-Stδ

(2)

(3)

(4)

(5)

式中E*——等效弹性模量

R*——等效半径

α——法向重叠量

ε——恢复系数δ——切向重叠量

Sn——法向刚度

m*——等效质量

vnrel——法向相对速度

St——切向刚度

vtrel——切向相对速度

另外，颗粒还受到切向力力矩和滚动摩擦力矩，其计算公式为

Tt=Ri(Ft+Ftd)

(6)

Tn=-μrFnRiωi

(7)

式中μr——滚动摩擦因数

Ri——质心到接触点间距离

ωi——接触点处物体单位角速度

颗粒在其所受力和力矩作用下发生移动和滚动。当切向力大于库伦摩擦力时颗粒在接触表面发生滑动，切向力矩与滚动摩擦力矩综合作用决定颗粒的滚动运动。

3 边界参数的系统化研究

为研究各边界参数对肥料休止角的影响是否显著，设计了休止角影响因素的Plackett-Burman筛选试验。

3.1 休止角影响因素的筛选与分析

3.1.1休止角影响因素的Plackett-Burman筛选试验

应用Design-Expert 8.0软件进行Plackett-Burman试验设计，参考文献[4,18-21]确定ABS仿真参数为：泊松比0.394，剪切模量8.9×108Pa，密度为1 060 kg/m3；尿素密度为1 330 kg/m3，其余待定8个仿真边界参数作为筛选对象，并给出各因素的范围如表1所示。仿真中时间步长设置为Rayleigh步长的25%，仿真计算网格尺寸为最小颗粒直径的4倍。

表1 Plackett-Burman试验因素Tab.1 Parameters of Plackett-Burman test

3.1.2筛选试验结果与分析

Plackett-Burman试验设计与仿真休止角结果如表2所示，试验中具有12组-1与+1水平点，1组中心点仿真试验。

表2 Plackett-Burman试验设计与结果Tab.2 Design and results of Plackett-Burman test

利用Design-Expert 8.0软件对试验结果进行方差分析，结果(表3)表明，对休止角影响重要性次序依次为尿素颗粒间滚动摩擦因数、尿素颗粒间静摩擦因数、尿素颗粒与ABS间静摩擦因数，其中，前两者高度显著。

表3 Plackett-Burman试验显著性分析Tab.3 Significance analysis of Plackett-Burman test

注：*表示显著，** 表示高度显著。

为探讨3种显著性边界参数对休止角的影响，模拟了不同参数下尿素颗粒堆积状态的变化，确定仿真参数为[4，19-21]：尿素剪切模量3.48×107Pa，泊松比0.25，尿素颗粒与颗粒间、与ABS间碰撞恢复系数分别为0.35和0.60，尿素颗粒与ABS间滚动摩擦因数为0.41，作为不变参数时颗粒间静摩擦因数为0.27，颗粒与ABS板间静摩擦因数为0.33，颗粒间滚动摩擦因数为0.22，其余参数见3.1.1节，仿真结果如图5～7所示。

图5 不同颗粒间滚动摩擦因数颗粒堆积图Fig.5 Particle packing diagrams of different rolling friction coefficients between particles

图6 不同颗粒间静摩擦因数颗粒堆积图Fig.6 Particle packing diagrams of different static friction coefficients between particles

图7 不同颗粒与ABS板间静摩擦因数颗粒堆积图Fig.7 Particle packing diagrams of different static friction coefficients between particles and ABS plates

图5a、5b、5c分别是尿素颗粒间滚动摩擦因数为0.01、0.26、0.50时尿素颗粒堆积状态图，休止角随颗粒表面滚动摩擦因数增大而增大，当滚动摩擦因数较小时颗粒堆较分散，主要由于颗粒下落堆积过程中，颗粒间接触作用剧烈，中心颗粒易于形成稳定的堆积状态，对边界颗粒产生排挤作用，边界颗粒向外部扩散呈松散状态，外层颗粒具有较大的旋转动能，对颗粒的转动惯量抑制作用较小，因此较小的颗粒间滚动摩擦因数时形成的休止角较小，颗粒堆较为分散。相反，较大的颗粒间滚动摩擦因数不利于边界颗粒的扩散，降低了外层颗粒旋转动能，对颗粒的转动惯量抑制作用较大，从而使颗粒沿中心轴即颗粒堆的高度方向堆积，利于形成较大的休止角。

图6a、6b、6c分别是尿素颗粒间静摩擦因数为0.20、0.45、0.70时颗粒堆积状态图，相比滚动摩擦因数，静摩擦因数较小时颗粒堆分散现象更为明显，较小的颗粒间静摩擦因数使颗粒堆发生变异或不易形成堆状，主要由于颗粒间静摩擦因数较小时，库仑摩擦力较小，切向力较小，导致切向力不足以支持颗粒自身重力，颗粒在接触表面易发生滑动，不易形成稳定堆积体；颗粒间静摩擦因数较大时，库伦摩擦力较大，颗粒在接触表面相对静止，形成稳定堆积体。

图7a、7b、7c分别是尿素颗粒与ABS板间的静摩擦因数为0.10、0.35、0.60时颗粒的堆积状态图，当颗粒与接触面间的静摩擦因数较小时，颗粒在接触面间不能形成稳定的颗粒堆，当摩擦因数为0.35时，所形成的尿素颗粒堆较为集中，边缘无分散现象。在散粒体堆积的过程中，重力作用使中底部散体颗粒间相互挤压变形，形成强力链，且相互链接形成架拱保证颗粒堆的稳定性。当颗粒与接触面间的静摩擦因数较小时，底部颗粒与接触面间的作用力不足以支撑上部颗粒的挤压变形而使底部的强力链破坏，进而颗粒堆分散；随着摩擦因数增大，底部颗粒与接触面间作用力变大，底部颗粒易于形成强力链，使得颗粒向高度方向堆积，形成较大休止角，颗粒受到的接触力作用线处于摩擦角的范围之内，力链中的颗粒处于自锁状态，提高了颗粒与接触面间接触稳定性，形成稳定的堆积体。

3.2 重要仿真边界参数的确定

上述分析结果表明，尿素颗粒间滚动摩擦因数以及颗粒间、颗粒与ABS间静摩擦因数是影响颗粒堆积过程的重要因素，因此三者取值是提高仿真精度的关键。

3.2.1尿素颗粒间静摩擦因数的标定

颗粒间细观参数很难通过试验方法直接准确测定，目前常用的方法是通过标定来确定[22]。采用自然休止角作为宏观特征，依据第1节中介绍的休止角测定方法，在EDEM中进行不同尿素颗粒间静摩擦因数下休止角仿真测定试验，尿素颗粒间滚动摩擦因数和尿素颗粒与ABS间静摩擦因数分别为0.30和0.33，其余仿真参数同3.1.2节，每组参数下取6组重复试验平均值作为最终测定的休止角，得休止角ψ与静摩擦因数μspp的关系如图8，其拟合方程为

ψ=35.883μspp+18.772 (R2=0.984 1)

(8)

图8 颗粒间静摩擦因数与休止角关系曲线Fig.8 Relationship curves between urea-urea static friction coefficient and repose angle

二者呈线性变化趋势，方程拟合精度较高(R2=0.984 1)，因此可用于标定尿素颗粒间静摩擦因数。

在与仿真相同试验条件下测定尿素颗粒实际休止角，试验对象同2.1节，休止角测定方法如第1节所述，得到尿素颗粒实际休止角均值30.43°。将试验测定值代入拟合方程(8)，计算得到标定的尿素颗粒间静摩擦因数μspp=0.32。

3.2.2尿素颗粒与ABS间静摩擦因数的测定

利用实验室自制静摩擦因数测量仪测定尿素颗粒与ABS板间静摩擦因数(图9)。将ABS板固定在支撑板上，为避免肥料颗粒滚落而非滑落，尽量选择扁平肥料(图9b)，调节丝杠螺母机构，当颗粒刚好下滑时记录升降板角度，其正切值即为尿素颗粒与ABS板间静摩擦因数，测量10次求平均值，为0.33。

图9 静摩擦因数测量Fig.9 Measurement of static friction coefficient

3.2.3尿素颗粒间滚动摩擦因数的标定

由于摩擦理论、测量设备的不完备，颗粒滚动摩擦因数还没有准确、成熟的测量方法。目前颗粒间接触参数(静摩擦因数和滚动摩擦因数)采用粘结在倾斜板上的颗粒板法直接测定比较困难，因此本文采用模拟标定方法近似获得颗粒滚动摩擦因数。

为了降低接触面对尿素颗粒堆积过程的影响，重点考察颗粒间摩擦特性对休止角的影响，提高颗粒间滚动摩擦因数的测定精度，参考文献[6]设计了如图10的颗粒堆积模拟和试验测定平台，无底细圆管高250 mm、直径18 mm，圆筒高40 mm、直径190 mm，事先将400 g尿素在底部圆筒内铺均匀，并将表面刮平，再将细圆管中填满尿素颗粒，竖直置于底部圆盘肥料表面，以0.05 m/s的速度将细圆管缓慢提起，使尿素颗粒在肥料水平面上滚落堆积，从而获得肥料堆如图10所示。休止角的测定均由Matlab读取颗粒堆积图像，参见上述休止角的测定方法。

图10 尿素颗粒堆积模拟和试验平台Fig.10 Simulation and experimental verification of urea particles’ accumulation

得到尿素颗粒间静摩擦因数和尿素颗粒与ABS间静摩擦因数分别为0.32和0.33，其余仿真参数同3.1.2节。在EDEM中进行不同肥料颗粒滚动摩擦因数下的休止角的测定，取6组重复试验平均值作为最终测定的休止角，得到肥料颗粒间休止角ψ′与滚动摩擦因数μrpp的函数关系如图11，其拟合方程为

ψ′=62.317μrpp+17.386

(9)

图11 颗粒间滚动摩擦因数与休止角关系曲线Fig.11 Relationship curves between urea-urea rolling friction coefficient and repose angle

二者呈线性变化趋势，方程拟合精度较高(R2=0.989 4)，因此可用于标定尿素颗粒间滚动摩擦因数。

利用休止角试验测定装置测定尿素颗粒休止角，试验选取2.1节中经筛分的大颗粒尿素，用烘干法测得含水率为0.11%。实际测定试验时与EDEM仿真计算保证试验条件一致，测得真实颗粒的休止角为30.30°。将试验测定值代入拟合方程(9)，计算得到标定的尿素颗粒间滚动摩擦因数μrpp=0.21。

3.3 休止角试验与仿真模型的验证

为验证尿素颗粒间滚动摩擦因数以及颗粒间、颗粒与ABS间静摩擦因数3个重要边界参数确定的可靠性，在三者分别为0.21、0.32和0.33时对尿素颗粒进行休止角的模拟仿真(其余边界参数同3.1.2节)，休止角仿真结果如图12a所示，测得休止角为30.32°，在相同试验条件下测定尿素颗粒实际休止角结果如图12b所示，测得休止角均值为30.43°，二者相差0.11°，相对误差为0.36%，模拟结果接近试验值，表明试验测定和标定后的EDEM边界参数和仿真模型的有效性，也证明对于难试验测定的边界参数通过离散元模拟进行参数标定的可行性。

图12 休止角仿真结果与试验结果对比Fig.12 Comparison of simulation and experiment results

为研究标定的3个重要边界参数对不同含水率大颗粒尿素的适应性，对不同含水率下尿素颗粒休止角进行了测量[23]。根据GB 2440—2001中农用尿素的含水率标准(小于等于1.0%)，试验中发现当含水率为1.1%时尿素颗粒开始出现粘结现象，因此文中研究的尿素颗粒最大含水率为1.2%。不同含水率下的休止角与本节中标定参数下仿真测定休止角间的相对误差如表4所示，可以看出在含水率小于1.2%时相对误差均低于3.25%，表明通过仿真得到的3个重要参数标定值适用于不同含水率下尿素颗粒的标定，因此仿真模型和颗粒间滚动摩擦因数、静摩擦因数、颗粒与ABS间静摩擦因数的标定值0.21、0.32和0.33可用于不同含水率下尿素颗粒的仿真研究。