基于离散元的包膜肥料Bonding模型参数标定

都 鑫 刘彩玲 姜 萌 袁 昊 戴 磊 李方林

(中国农业大学工学院， 北京 100083)

0 引言

肥料在农业生产中有着不可替代的作用，中国化肥使用量大但肥料利用率低[1]。为保证农业生产的可持续发展，养分效率高、环境友好型的控释肥料得到广泛关注[2-6]。控释肥料采用聚合物包膜，可定量控制肥料养分释放量和释放期，使养分供应与作物各生育期需肥规律吻合，节肥增效作用显著[7]。

包膜肥料的养分释放特性与包膜层材料、结构紧密相关[8]，排肥器在排施肥料颗粒过程中会对部分肥料颗粒造成不同程度的机械损伤，导致包膜层破坏进而影响肥料养分释放特性。为设计适用于包膜肥料无损排施的排肥器，考虑包膜肥料与排肥器、包膜肥料间的相互挤压力、摩擦力等因素，改进排肥方式、结构参数等。其中肥料颗粒间的相互作用力难以直接测得，拟采用离散元方法对此过程进行数值模拟[9]。

自CUNDALL等[10]提出适应于土力学的离散元方法后，该方法已广泛应用于散粒体物料的混料、搅拌、输送、充填、破碎等过程[11-14]。POTYONDY等[15]提出粘合粒子模型(BPM)，将许多小颗粒通过Bond键粘结形成被破碎的物体，在外力作用下Bond键发生形变，当力或力矩达到一定程度时，粘结键断裂，小颗粒互相分开发生破碎。目前物料破碎数值模拟被广泛应用于岩石钻探、土壤挖掘、制沙、磨粉等[16-20]研究中，但在农业领域中应用较少，尚未见到包膜肥料破碎数值模拟方面的相关报道。

为保证数值模拟可真实反映肥料颗粒挤压破碎情况，需对包膜肥料颗粒破碎模型进行参数标定。王云霞等[21]通过堆积角试验标定玉米种子的种间静摩擦因素和滚动摩擦因数；鹿芳媛等[22]基于摩擦角试验标定了不同含水率下的水稻芽种主要接触参数；刘彩玲等[23]基于休止角试验系统标定了大颗粒尿素的仿真边界参数；袁全春等[24]基于JKR模型标定了有机肥散体颗粒离散元参数；于庆旭等[25]通过碰撞弹跳试验、斜面滑移试验和斜面滚动试验测定三七种子与ABS塑料间接触参数，基于响应面优化方法得到三七种子间接触参数。

本文在国内外专家学者的研究基础上[26-30]，以包膜肥料颗粒为研究对象，采用单轴压缩试验测得包膜肥料颗粒的极限破碎位移和极限破碎载荷，通过建立与真实试验相同的仿真模型，利用Placket-Burman试验、Steepest ascent试验和Box-Behnken试验逐步减少因素数量、缩小因素范围，最终优化得到与真实试验结果最接近的Bonding模型参数组合；为进一步验证最优参数组合的正确性，利用外槽轮排肥器排施包膜肥料颗粒，并统计对比仿真和实际条件下肥料颗粒破碎率，以期为研究包膜肥料颗粒机械破碎机理、优化设计无损排施的新型排肥器提供参考。

1 包膜肥料颗粒Bonding模型建立

1.1 Bonding模型

POTYONDY等[15]设计了Bonding接触模型以模拟物料破碎。在Bonding接触模型中，被破碎的物料(颗粒、块体等)由N个较小的颗粒构建组合而成，小颗粒之间通过Bond键固定在一起。Bond键具有与有限元方法相近的力学特性，在外力作用下会发生形变(拉伸、压缩、扭转)。当Bond键形变产生的力或力矩达到一定程度时，粘结键断裂，小颗粒互相分开发生破碎。Bonding接触模型假设两个颗粒之间的粘结键是一个虚拟的平板(圆柱)，如图1所示。

图1 Bonding模型原理与颗粒粘结示意图Fig.1 Bonding model and particle bonding schematic

当两个颗粒发生相对运动时，平板(圆柱)承受拉伸、弯曲和剪切，键产生力的计算公式为

(1)

其中

(2)

式中 δFn——法向力，N

δFt——切向力，N

vn——法向速度，m/s

vt——切向速度，m/s

Sn——法向刚度，N/m

St——切向刚度，N/m

δt——仿真时间步长，s

A——颗粒间的接触面积，m2

δMn——法向力矩，N·m

δMt——切向力矩，N·m

ωn——法向角速度，rad/s

ωt——切向角速度，rad/s

J——颗粒惯性矩，mm4

RB——粘结半径，mm

当力和力矩达到极限值或Fraction颗粒间的距离大于设定的接触半径时，粘结键发生断裂，法向及切向剪切力计算公式为

(3)

式中σmax——法向剪切力，N

τmax——切向剪切力，N

1.2 肥料颗粒离散元建模

选用山东农泱生物科技有限公司生产的贝可施包膜控释肥料(图2a)，肥料含水率0.88%，真实密度1.46 g/cm3，平均三轴尺寸4.08 mm×3.97 mm×3.89 mm，等效直径为3.98 mm，球形率为0.975，选取长、宽、高与等效直径相近的肥料颗粒建立轮廓模型。为研究肥料颗粒破碎特性，采用粘接颗粒模型(Bonded particle method, BPM)建立肥料颗粒离散元模型(图2b)，通过Bond键将多个直径相等的小球形填充颗粒粘结成一个大颗粒。

图2 肥料颗粒及离散元模型Fig.2 Fertilizer particles and discrete element model

填充颗粒直径越小，粘结成一个大颗粒所需的颗粒数就越多。如图3a所示，填充的球形颗粒半径为0.4 mm，球形颗粒数量为68，粘结键个数为534；如图3b所示，填充的球形颗粒半径为0.2 mm，球形颗粒数量为582，粘结键个数为3 611；如图3c所示，填充的球形颗粒半径为0.1 mm，球形颗粒数量为5 581，粘结键个数为39 119。

图3 不同直径球形颗粒填充图Fig.3 Diagrams of particles with different diameters

由于小球形颗粒间彼此相互独立，在仿真模拟中需要对每个小球形颗粒的受力、运动情况进行计算，因此，颗粒数量越多，仿真计算时间越长，综合考虑后选取球形颗粒半径为0.2 mm的粘结模型进行Bonding模型参数标定。

1.3 单轴压缩试验

1.3.1真实试验

选取直径为3.98 mm的肥料颗粒，采用万能试验仪对肥料颗粒进行单轴压缩试验，以速度0.1 mm/s进行加载，直至试样出现破坏后停止，单轴压缩试验中试件加载过程如图4所示。试验重复20次取平均值，得到肥料颗粒的位移载荷曲线如图5所示。

图4 万能试验仪Fig.4 Universal tester1.压板 2.底板 3.控制面板

图5 力-位移曲线Fig.5 Force-displacement curve

由图5可以看出，肥料颗粒力-位移曲线分为3部分，第1部分为试验开始阶段，试样局部发生棱角破坏，但由于肥料颗粒球形度较高，力先快速增大后变缓；第2部分为主体发生破坏的过程，当所施加的载荷超过试样所能承受的极限时，试样主体发生破坏，力达到最大值，此时为破坏点；第3部分随位移增加，力持续降低，试样破碎成块状。肥料颗粒加载位移为0.45 mm时加载应力最大值为58.668 N，即20粒直径为3.98 mm的包膜肥料颗粒平均极限破碎载荷为58.668 N，此时对应的极限破碎位移为0.45 mm。

1.3.2仿真试验

采用离散元软件EDEM 2018建立与实际试验一致的单轴压缩仿真模拟试验，如图6所示。仿真中用到的包膜肥料颗粒本构参数、包膜肥料间和包膜肥料与ABS间的接触参数由实验室标定得到，ABS的本构参数如表1所示[23]。EDEM中设置仿真总时间为8 s，时间步长为7.5×10-7s，每隔0.001 s保存1次数据，网格尺寸取3.5倍最小球形单元尺寸。

表1 材料物理特性参数Tab.1 Physical characteristic parameters of material

图6 单轴压缩试验仿真图Fig.6 Simulation diagram of uniaxial compression test

2 Bonding模型参数标定

2.1 Plackett-Burman试验

Plackett-Burman试验通过考察目标响应与各因子间的关系，比较各个因子2水平间的差异性来确定因子显著性。以极限破碎位移Y1和极限破碎

载荷Y2为目标，以单位面积法向刚度X1、单位面积切向刚度X2、临界法向应力X3、临界切向应力X4和粘结半径X5为Plackett-Burman试验的5个因子，每个因子选取2个水平。试验选用N=11的Plackett-Burman设计，预留6个虚拟项作误差分析，试验因素水平如表2所示。

表2 Plackett-Burman试验因素水平Tab.2 Plackett-Burman test factors and levels

Plackett-Burman试验方案与结果如表3所示，利用Design-Expert 8.0.6软件对各因素效应分别进行方差分析和t检验，结果如表4和图7所示。

表3 Plackett-Burman试验结果Tab.3 Plackett-Burman test results

由表4和图7可知，各因素对极限破碎位移Y1影响由大到小为X5、X2、X1、X4、X3，其中X5、X2、X1的P值小于0.05，表明因素对目标影响显著，由帕累托图的t-value检验也可得出同样结论，且知X5、X1对目标的效应值为正，X2对目标的效应值为负。同理，可知各因素对极限破碎载荷Y2影响由大到小为X5、X1、X4、X2、X3，其中X5、X1、X4的P值小于0.05，表明因素对目标影响显著，由Pareto Chart的t-value检验也可得出同样结论，且知X5、X1、X4对目标的效应值均为正。

图7 帕累托图Fig.7 Pareto Charts

表4 Plackett-Burman试验结果显著性分析Tab.4 Analysis of significance of parameters in Plackett-Burman test

各因素与极限破碎位移Y1和极限破碎载荷Y2的模型为

(4)

2.2 Steepest.ascent试验

根据Plackett-Burman试验结果，X3对结果影响不显著，因此将X3取为5×105Pa。选取因素X1、X2、X4、X5初值分别为6.40×109N/m3、1.92×1010N/m3、8.00×105Pa、0.16 mm，步长分别为1.60×109N/m3、-3.20×109N/m3、2.00×105Pa、0.02 mm，进行Steepest ascent试验，以进一步找到逼近真实值的参数组合，Steepest ascent试验方案和结果如表5所示。

表5 最速上升试验方案及结果Tab.5 Steepest ascent test program and results

由1.3.1节可知，肥料颗粒实际单轴压缩极限破碎位移和极限破碎位移载荷分别为0.45 mm、58.668 N，由表5可知，仿真试验中极限破碎位移和载荷与真实值间的误差先减小后增大，第3组试验的极限破碎位移Y1和极限破碎载荷Y2与真实值最接近。

2.3 Box-Behnken试验

为得到最接近真实值的极限破碎位移Y1和载荷Y2，根据Steepest ascent试验结果，选取第2和第4组试验的因素组合为BBD试验的高低水平，试验方案及结果如表6所示，方差分析结果如表7所示。

由表7可知，Y1和Y2的二次回归模型均高度显著，且失拟项不显著，表明模型精度高。X1、X2、X4、X5均显著且对Y1的影响由大到小为X5、X1、X2、X4。X1、X2、X4、X5均显著且对Y2的影响由大到小为X5、X2、X1、X4。

表7 方差分析Tab.7 Analysis of variance

为准确求得各因素最优参数组合，需尽量缩小约束条件范围，已知Y1和Y2目标值分别为0.45 mm、58.668 N，由表6可知，试验方案2、13、15和19的Y1和Y2最接近实际值，且X5均为0.20 mm，根据方差结果可知X4对Y1和Y2的影响程度最小，所以设定约束条件范围为

(5)

利用Design-Expert 8.0.6软件进行优化求解，得到X1、X2、X4、X5分别为9.978 25×109N/m3、1.033 75×1010N/m3、1.20×106Pa、0.2 mm时，Y1、Y2分别为0.450 mm、58.668 N，对优化求解得到的参数组合进行仿真验证，得到Y1、Y2仿真值分别为0.451 mm、58.993 N，与真实值的相对误差分别为0.222%、0.554%。

3 肥料颗粒破碎率验证试验

3.1 试验条件

为进一步验证优化得到的肥料颗粒Bonding模型参数组合的可靠性，选取目前播种机上常用的外槽轮排肥器进行排肥性能试验，试验肥料为前文包膜肥料颗粒并过筛(孔径3 mm)后备用，如图8a所示。每次试验称量1 500 g粒径大于3 mm的完整肥料颗粒，经外槽轮排肥器完全排出后再次过筛(孔径3 mm)，统计粒径小于3 mm的肥料质量，在不同排肥轮转速下进行试验，每组试验重复3次取平均值。采用三维软件SolidWorks 2018对排肥器进行建模并保存为.IGS文件，将模型导入离散元软件EDEM 2018中，如图8b所示。在与真实试验相同的排肥轮转速、槽轮有效工作长度条件下进行虚拟仿真试验，各条件下肥料颗粒破碎率η计算公式为

图8 肥料颗粒破碎试验Fig.8 Fertilizer particles crushing test1.肥料 2.接肥盒 3.步进电机 4.仿真肥料颗粒 5.外壳 6.排肥轮 7.排肥舌

(6)

式中m1——粒径小于3 mm的肥料质量，g

m0——每次试验肥料总质量，取1 500 g

3.2 试验结果

在排肥轮转速分别为10、15、20、25、30 r/min时得到肥料颗粒真实破碎率和仿真破碎率如图9所示。由图9可以看出，肥料颗粒仿真破碎率与真实情况相比误差为9.42%～11.40%，说明本文标定得到的Bonding模型参数组合准确性高。

图9 不同转速下肥料颗粒破碎率Fig.9 Fertilizer particles crushing rate at different rotational speeds

肥料颗粒在排肥轮和排肥舌的挤压、摩擦作用下，颗粒间粘结键断裂，整个颗粒从中间破碎为两大块，颗粒破碎仿真过程如图10所示。由图10d～10f可以看出，在排肥轮和肥料刚接触时，肥料除受挤压力外，还伴随滚动，说明实际情况下肥料颗粒包膜材料还会发生研磨现象，此过程不会破碎肥料颗粒，但在一定程度上破坏包膜结构，进而影响肥料养分释放过程。

图10 排肥过程中肥料颗粒破碎仿真过程Fig.10 Fertilizer granules crushing simulation

4 结论

(1)通过单轴压缩试验得到包膜肥料颗粒的极限破碎位移和极限破碎载荷分别为0.45 mm、58.668 N。

(2)Placket-Burman试验结果表明，粘结半径、单位面积切向刚度和单位面积法向刚度对极限破碎位移影响显著，粘结半径、单位面积法向刚度和临界切向应力对极限破碎载荷影响显著。

(3)Box-Behnken试验结果表明，粘结半径、单位面积法向刚度、单位面积切向刚度、临界法向应力、临界切向应力和粘结半径均对极限破碎位移和极限破碎载荷影响显著。优化求解得到当单位法向刚度、单位面积切向刚度、临界切向应力和粘结半径分别为9.978 25×109N/m3、1.033 75×1010N/m3、1.20×106Pa、0.2 mm时，极限破碎位移和极限破碎载荷分别为0.450 mm、58.668 N。最优参数条件下得到极限破碎位移和极限破碎载荷与真实值的相对误差分别为0.222%、0.554%。

(4)外槽轮排肥器排肥结果显示，在排肥轮转速分别为10、15、20、25、30 r/min时肥料颗粒真实破碎率和仿真破碎率相对误差不大于11.40%，说明本文标定得到的Bonding模型参数组合准确性高。