黑水虻生物转化餐厨垃圾有机肥离散元模型参数标定

宋世圣，孙松林，方芹，彭才望，周婷，朱海英

（湖南农业大学机电工程学院，长沙 410128）

随着人口数量的增加和人们饮食习惯的变化，餐厨垃圾的量日益增多。餐厨垃圾中含有大量的有机质，如处理不当易造成环境破坏，从而影响人们的健康，餐厨垃圾如何有效且洁净地处理就显得尤为重要[1-4]。

黑水虻近年来被广泛关注，其通过虫体的取食与消化，可将餐厨垃圾转化为含有高质量蛋白质的生物质[5-7]。现有关于黑水虻生物转化餐厨垃圾有机肥的研究主要集中在黑水虻幼虫的生长发育以及生物转化餐厨垃圾的转化率等方面[8-9]，在黑水虻生物转化餐厨垃圾有机肥的离散元仿真模型方面研究极少，影响了后续收集、输送、筛分机械的研究，降低了生产效益。

离散元素法(discrete element method，DEM)是新型数值计算方法，其认为整体是由大量离散个体组成，同时各个个体之间存在靠拢与离散2种形式，可用于分析与求解复杂系统的动力学问题[10-11]。物料离散元参数标定是研究物料与机械相互作用的基础，韩树杰等[12]利用漏斗注入法对新疆果园深施散体厩肥进行了离散元参数标定，为机械化深施厩肥离散元仿真提供了理论依据；
王韦韦等[13]对玉米秸秆粉料致密成型的离散元模型进行了参数标定，通过模孔压缩对比试验验证了其标定参数的准确性，可为秸秆粉料致密成型的仿真过程提供基础参考；
石辰风等[14]基于“Hertz-Mindlin with JKR Cohesion”模型对4种中药浸膏粉末进行了参数标定，同时对比了不同硬件配置对仿真计算结果的影响，可以对中药制药中的药粉混合和输送过程进行精确离散元模拟；
彭才望等[15]利用圆筒提升法对黑水虻处理的猪粪有机肥进行了参数标定，为该物料收集机械的设计提供一定理论参考；
邢洁洁等[16]对海南热区砖红壤颗粒进行了离散元模型参数标定，通过破土阻力试验，验证了其标定参数的可靠性，在离散元软件中，可以较好地模拟触土部件与海南热区砖红壤的接触作用。

目前，关于黑水虻生物转化餐厨垃圾的研究大多围绕在转化效率[17]及生态效益[18]等方面，关于该物料的转运、输送研究较少。而通过对餐厨垃圾有机肥离散元模型参数的标定，可以对机械化收集、转运物料的过程进行离散元模拟仿真，极大地减少收集机械的设计研究时间，降低劳动力成本，提高生产效益。因此，本研究利用仿真与物理相结合的方法，构造餐厨垃圾有机肥模型，采用图像处理方法精确测量堆积角，依次通过Plackett-Burman试验、最陡爬坡试验、Box-Behnken试验确定物料的离散元参数数值选择，可为黑水虻生物转化餐厨垃圾有机肥在收集、转运过程中的离散元仿真提供一定理论依据。

1.1 试验材料

黑水虻生物转化餐厨垃圾有机肥（以下简称餐厨垃圾有机肥）取自湖南农业大学黑水虻科研基地（28.18°N，113.06°E），由含水率70%左右的餐厨垃圾经黑水虻幼虫处理10～12 d后形成，该餐厨垃圾主要成分包括肉类、青菜、辣椒、土豆、莲藕、主食等。餐厨垃圾中因含有骨头等硬物，会影响后续处理设备的正常运行[17-18]，因此在餐厨垃圾处理前，使用小型粉碎机进行粉碎，黑水虻生物转化处理后的餐厨垃圾有机肥呈颗粒状，经烘干法测得其含水率为37.5%。由于在餐厨垃圾中存在辣椒籽等黑水虻无法转化的杂质，通过预测量发现，有机肥颗粒粒径一般在2 mm以下，杂质颗粒粒径一般在2.6 mm以上，因此使用10目标准筛（筛孔2 mm）对有机肥进行筛分，保留筛下物作为堆积角试验的材料。将颗粒分作10组，每组200个颗粒，使用数显游标卡尺（精度0.02 mm）对粒径进行测量。粒径分布如图1所示，走势同正态分布，大多数颗粒粒径在1.1～1.3 mm之间，粒径均值1.2 mm，标准差为0.25 mm。

图1 粒径分布Fig.1 Particle size distribution

1.2 接触参数测定

1.2.1 静摩擦系数的测定 静摩擦系数是物料受到最大静摩擦力与法向正压力的比值，可以较好地反映物料与接触表面的摩擦特性[19]。首先将倾角仪调0，为防止有机肥颗粒较小而产生的滑动摩擦及有机肥颗粒间的黏结性，因此，使用上下无底的方形模具插入餐厨垃圾有机肥堆体，使其形成长4 mm、宽4 mm、高5 mm的有机肥立方体[27]，放置在光滑的钢板（图2）上，钢板下放置1个滑块，通过推动滑块调节钢板倾斜角度，在有机肥立方体与钢板产生相对滑动的瞬间，停止推动，立即记录倾角仪读数，试验重复10次，通过式（4）得到餐厨垃圾有机肥与钢板的静摩擦系数。同时在测定餐厨垃圾有机肥之间的静摩擦系数时，只需在钢板上利用双面胶粘上0.5 mm厚的有机肥颗粒排即可，粘结时尽可能让颗粒排布紧密，以免黏结剂对摩擦系数产生影响[25,28]。

图2 静摩擦系数测定原理Fig.2 Principle of static friction coefficient measurement

式中，G为有机肥立方体重力（N）；
N为斜面对有机肥立方体的支撑力（N）；
f为斜面对有机肥立方体的最大静摩擦力（N）；
μ1为静摩擦系数；
θ1为有机肥立方体滑动瞬间的倾角（°）。

1.2.2 滚动摩擦系数的测定 餐厨垃圾有机肥相较于钢板弹性模量较小，因此在钢板滚动时，可忽略微观变形和弹性滞后，只有少许滑动影响试验结果，且倾斜角度越大，滑动摩擦影响越大[19]。因此在进行滚动摩擦系数测定试验前，使用有机肥圆柱体进行预滚动，测得有机肥滚动的最小角度为5.7°。

利用直径为10 mm的圆形环刀将餐厨垃圾有机肥制成半径5 mm、高20 mm的圆柱体，由钢板顶端滚动至钢板底端，利用相机拍摄有机肥圆柱体在钢板滚动的全过程，将视频导入Premiere 2018软件中，选取不同帧数记录1、2、3、4点的位置（图3），1、2距离为L1，3、4距离为L2，由于L1、L2距离较小，可以求得中点a、b的瞬时速度，通过式（5）得到餐厨垃圾有机肥与钢板之间的滚动摩擦系数[20]，试验重复10次。有机肥与有机肥之间的滚动摩擦系数测定同1.2.1。

图3 滚动摩擦系数测定原理Fig.3 Measuring principle of rolling friction coefficient

式中，μ2为滚动摩擦系数；
m为有机肥的质量（kg）；
h为有机肥下落高度（m）；
vb为b点瞬时速度（m·s-1）；
va为a点瞬时速度（m·s-1）；
L3为a点到b点位移（m）；
g为重力加速度9.81（m·s-2）。

1.3 堆积角物理测定

有机肥堆积角试验中漏斗下端管口直径18 mm，漏斗下端管口与钢板距离100 mm。试验时，利用挡板将漏斗口堵住，将有机肥颗粒缓慢倒入漏斗，待稳定后移开挡板，颗粒经漏斗口落在钢板上形成颗粒堆，待该颗粒堆稳定后，利用相机拍摄颗粒堆视图，将照片导入Photoshop 2018中，在+x及+y方向分别标注角度取平均值即为该颗粒堆的堆积角，重复10次取平均值，测得堆积角为41.34°。

1.4 EDEM仿真模型构建及仿真参数设置

餐厨垃圾有机肥颗粒外观近似球形，因此将有机肥颗粒仿真模型设置为球形，有机肥颗粒平均粒径为1.2 mm，设置动态颗粒工厂，颗粒尺寸以正态分布生成，粒径标准差0.25 mm，选择Cap，Lower设置为0.5，Upper设置为1.5。在SolidWorks中根据物理试验模型1∶1进行构建，以IGS文件格式导入EDEM 2020，模型材料为钢，泊松比为0.30，密度为7 860 kg·m-3，剪切模量为7.9×1010Pa[15]，如图4所示。

图4 堆积角仿真模型Fig.4 Simulation model of accumulation angle

设置factory在漏斗上方以virtual状态生成颗粒，物理有机肥颗粒百粒重在0.511～0.585 g之间，均值为0.534 g。在物理试验时取100 g有机肥颗粒进行堆积角试验，为与物理堆积角试验条件一致，在构建仿真模型时，有机肥颗粒共18 000个，生成速率为9 000个·s-1，生成颗粒时，挡板堵住漏斗口，2.5 s时待颗粒在漏斗中稳定后，移开挡板，使颗粒由漏斗口落在钢板上。固定时间步长为Rayleigh时间步长的22%，Cell Size为3 Rmin，数据保存的时间间隔为0.01 s。由于餐厨垃圾有机肥较易粘结且具有一定湿度，故选取“Hertz-Mindlin with JKR”模型作为此次离散元仿真试验的接触模型[21]。

1.5 参数标定试验设计

根据斜面试验及已有研究进展[12-16,22-25]设置参数水平（表1），利用Design Expert 10.0软件设计Plackett-Burman试验进行优化，以相对误差为评价指标，当相对误差最小时说明该组参数即为各组参数中的最优选择，从而得到参数的最优区间。根据上述试验确定参数范围，以显著性参数为试验因素，以堆积角为响应值，设计Box-Behnken试验。在Design Expert软件的响应优化模块中，以餐厨垃圾有机肥物理堆积角41.34°为目标值对回归模型进行寻优求解，得到参数的最优组合，在该参数组合下进行仿真试验，重复试验3次，与物理试验结果进行对比。

表1 Plackett-Burman试验参数Table 1 Factors of Plackett-Burman Design

1.6 数据处理

采用Design Expert 10.0软件进行试验设计及数据分析，采用Origin 2017进行绘图，采用Matlab 2018b进行图像处理。

2.1 Plackett-Burman筛选试验

Plackett-Burman试验筛选参数结果如表2所示，可以看出，仿真试验堆积角数值跨度较大，在30.2°～62.85°之间，与物理试验堆积角的相对误差较大，因此需进行方差分析筛选显著性参数。

表2 Plackett-Burman试验设计及结果Table 2 Design and results of Plackett-Burman Design

方差分析结果如表3所示，可以看出，模型P＜0.05，表明模型显著。对响应值有正向效应的因素有A、E、G、H，说明在其他因素不变时，随着该因素的增大，响应值也随之增大。对响应值有负向效应的因素有B、C、D、F、J、K，说明在其他因素不变时，随着上述因素的增大，响应值会逐渐减小。D的P值小于0.01，说明D因素对于堆积角的影响极显著，G、J的P值均小于0.05，说明其对堆积角的影响显著，其余因素对堆积角的变化均无显著影响，比较各参数的P值大小，各因素显著性排序为D、J、G、B、A、H、K、C、E、F。故选取对堆积角影响较大的D、G、J进行后续试验。

表3 Plackett-Burman试验结果显著性分析Table 3 Significance analysis of Plackett-Burman test results

2.2 最陡爬坡试验结果分析

由Plackett-Burman试验筛选得到3个显著性参数，根据标准化效应在一定步长范围内增大或减小，可以使仿真的结果快速逼近物理试验值[25]。在仿真试验中，改变显著参数D、G、J的取值，而其他非显著参数取中间值，即A为0.2、B为7 MPa、C为2 000 kg·m-3、E为0.8、F为0.15、H为0.65、K为0.1 J·m-2。通过表4可以看出，随着参数G的增大和参数D、J的减小，仿真堆积角与物理堆积角相对误差的走势呈现先减小后增大的趋势，其中最小误差为4.74%，最大误差则达到48.16%。而第2组参数组合下得到的堆积角相对误差最小，因此分别选择D为0.58、G为0.42、J为0.09作为后续Box-Behnken试验的中心参数组合水平，第1与第3组的参数组合分别为低水平和高水平。

表4 最陡爬坡试验设计与结果Table 4 Design and results of steepest ascent test

2.3 Box-Behnken试验结果分析

以3个显著性参数D、G、J为试验因素，堆积角作为响应值，设计Box-Behnken试验，其中，5组为中心水平重复试验，以探究D、G、J对堆积角的影响效果，试验设计方案及结果如表5所示。对试验结果进行多元回归分析，得到堆积角回归方程（式3）。

表5 Box-Behnken试验设计及结果Table 5 Design and results of Box-Behnken Design

方差分析结果如表6所示，可以看出，该模型P＜0.01，可以利用该回归模型对有机肥颗粒堆的堆积角进行准确的预测。失拟项P＞0.05，不显著，说明此回归方程中不存在失拟项。各因素显著性由高到低排序为J＞D＞G，其中D、J对堆积角影响极显著，G对堆积角影响显著，各因素交互作用对响应值均无显著性影响效果，D2、G2对堆积角有极显著影响，其余因素二次方项对响应值并无显著影响效果。变异系数为2.92%，说明此试验具有一定的可靠性。试验精确度Ap为14.391，说明试验具有较高的精确度。决定系数R2为0.967 8，校正决定系数为0.926 4，均接近于1，表明该回归方程拟合程度较好，可以用于堆积角预测分析。

2.4 参数最优组合确定及验证

以餐厨垃圾有机肥物理堆积角41.34°为目标值对回归模型进行寻优求解，得到最优组合为：D为0.46、G为0.51、J为0.09；
其他参数数值选取依次为：A为0.2、B为7 MPa、C为2 000 kg·m-3、E为0.8、F为0.15、H为0.65、K为0.1 J·m-2，进行仿真试验，利用EDEM后处理模块的测量工具测得3次仿真试验的堆积角分别为41.85°、43.4°、43.28°。3次仿真试验的平均值为42.84°，与物理试验堆积角（41.34°）的相对误差为3.63%，利用Matlab对图像进行处理，提取轮廓曲线（图5），可以看出，通过标定后得到的参数进行的仿真试验堆积边界与物理试验较为接近。t检验P值为0.094（＞0.05），证明仿真试验结果与物理试验结果无显著性差异。

图5 物理试验及仿真模拟对比Fig.5 Comparison of physical test and simulation test

目前围绕如何提高黑水虻转化餐厨垃圾效率的研究较多，关于机械收集、转运餐厨垃圾有机肥的研究较少。张杰等[29]对黑水虻资源化处理有机废弃物进行了阐述，认为黑水虻处理城市餐厨垃圾领域具有较大的应用前景，但在黑水虻对餐厨垃圾生物转化完成后如何回收地槽或料盘中的有机肥颗粒是目前的难题。本课题组前期设计了斗式取料机[30]及双向螺旋装置[31]，一定程度上可以收集、转运餐厨垃圾有机肥，但在设计过程中由于缺乏离散元模型参数，产生了设计进度慢、试验劳动强度大等一系列问题。因此本文以堆积角为评价指标，通过对照比较仿真试验与物理试验的堆积角数值，得到餐厨垃圾有机肥颗粒的离散元模型参数，仿真试验与物理试验的相对误差仅为3.63%，可以为有机肥收集机械的离散元仿真提供一定数据参考。

离散元法可以精确地模拟散体物料的运动形态，但农业物料的含水量往往对其物理特性具有较大影响，而本文仅对含水率为37.5%的餐厨垃圾有机肥进行离散元模型参数标定，后续将进一步探究不同含水率餐厨垃圾有机肥的离散元参数变化规律，完善本领域的离散元仿真数据库。另外需注意，在收集、转运有机肥颗粒时，物料内部存在一定数量的五龄期黑水虻幼虫，尚为活体，生物特性较为复杂，蔡静[32]对活体蝇蛆采用喷洒化学试剂及降温法使其静止，从而建立了两种姿态的蝇蛆离散元仿真模型，但由于蝇蛆与黑水虻之间的体态特征和生物特性存在较大差异，因此黑水虻幼虫与有机肥颗粒之间的接触参数有待进一步研究。