基于JKR粘结模型的蚯蚓粪基质离散元法参数标定

罗　帅　袁巧霞　GOUDA Shaban　杨龙元

(1.华中农业大学工学院， 武汉 430070； 2.农业部长江中下游农业装备重点实验室， 武汉 430070；3.本哈大学农学院， 本哈 13736)

0　引言

蚯蚓堆肥是处理畜禽粪便的有效方式[1]，经过蚯蚓过腹处理的禽畜粪便具有良好的孔隙特性和丰富的营养物质，成为优良的作物育苗和栽培基质，即蚯蚓粪基质。蚯蚓堆肥在禽畜粪便的综合利用中效果显著，逐渐成为研究的热点[2]。但目前蚯蚓堆肥过程中从布料、取料、蚓粪分离到蚯蚓粪基质的应用，基本上靠人工进行，大量人工成本已成为制约蚯蚓堆肥的主要因素，因此需研究蚯蚓堆肥处理及蚓粪利用相关设备。而研发其设备，就应该对蚯蚓粪基质相关特性进行充分了解。

蚯蚓粪基质是一种典型的散体物料，离散元法(Discrete element method, DEM)被广泛用于这些散体物料特性的研究[3-5]。JKR颗粒粘结模型是由JOHNSON等[6]于1971年提出并以JOHNSON、KENDALL和ROBERTS等的名字命名的模型，该模型引入了颗粒间表面能(Surface energy)的概念，适用于模拟细小和潮湿材料颗粒间的黏聚作用。JKR粘结模型自被提出以来，已经被许多学者引入到不同散体颗粒间的相互作用模型中，如黄豆和槭树残渣的混合搅拌运动的数值模拟[7]、型砂的流动性仿真模拟[8]、土壤堆积的数值模拟[9]、煤粉和生物质混合物休止角的推导[10]等，这些研究模拟结果和试验结果均较为吻合。

休止角(Angle of repose，AoR)作为散体物料的一种固有属性，与散体物料颗粒本身的各种属性息息相关，常被用于散体物料颗粒参数标定[11-14]。蚯蚓粪基质含水率变化时，其物理参数也会发生变化，这些变化可以通过休止角的变化来体现。

本文在对国内外散体物料参数标定和休止角测定方法进行总结的基础上，结合蚯蚓粪基质堆积体实际轮廓，提出一种休止角测定方法，基于JKR粘结模型进行蚯蚓粪基质的离散元法参数标定，设计休止角测定仪并通过试验测定不同含水率下蚯蚓粪基质的休止角，得到蚯蚓粪基质颗粒含水率与休止角之间的关系。通过测定蚯蚓粪基质含水率来预测其休止角，继而合理推定蚯蚓粪基质颗粒间滚动摩擦因数等其他参数。

1　材料与方法

1.1　试验材料

试验中所用的蚯蚓粪基质原料取自武汉市东西湖区新沟镇蚯蚓养殖基地，以牛粪为主要原料，经大平2号蚯蚓过腹处理，充分腐熟后得到可以用于植物育苗的蚯蚓粪基质原料。基质颗粒近似为球体。考虑蚯蚓粪基质实际应用时的处理方法，试验时将蚯蚓粪基质进行粉碎处理，全部通过孔径为2 mm的标准筛(浙江上虞市五四仪器厂生产)，以满足育苗钵成型试验粒径要求。蚯蚓粪基质原料pH值为5.89，电导率为0.22 mS/cm。

在不同含水率下休止角的测定试验中，将蚯蚓粪基质在遮荫条件下自然风干，使蚯蚓粪基质含水率缓慢下降，不定时取样分别测定休止角与含水率。

1.2　试验设备

参考目前广泛使用的休止角测定方法，设计了图1所示的钢质休止角测试仪。落料漏斗由高度调节螺杆上的螺母限制与堆积底座间的高度差。试验中设置漏斗出料口与堆积底座上表面距离为150 mm。堆积底座上端为直径150 mm的圆板。试验时，漏斗中的蚯蚓粪基质颗粒经漏斗口落于堆积底座上，在侧面对堆积体进行拍照，对照片进行处理以获取蚯蚓粪基质颗粒的休止角。

图1　休止角测试仪Fig.1　AoR determination instrument1.落料漏斗　2.高度调节螺杆　3.堆积底座

1.3　离散元模型

为缩短数值模拟时间，需要对模型进行适当的简化。将休止角测定仪仅保留漏斗口和堆积底座，得到其简化模型(图2a)。经观察，蚯蚓粪基质颗粒大部分为近似的球体(图2b)，因此以半径为1 mm的球体为颗粒原型(图2c)，同时在生成颗粒时，将原型颗粒半径变化范围设置为满足平均值，标准差为0.1 mm的标准正态分布，得到蚯蚓粪基质颗粒的离散元简化模型。

图2　测定装置及颗粒离散元模型Fig.2　Measuring device and DEM model of particles

综合考虑数值模拟的效率和准确性，时间步长通常设置为雷利时间步长的5%～40%，本文在模拟中视雷利时间步长的具体值，设置时间步长为雷利时间步长的5%～15%。为保证仿真结果的准确性，网格尺寸设置为最小颗粒半径的2倍。颗粒生成后经漏斗口落于堆积底座上自然堆积，仿真时间统一设置为3 s。

JKR颗粒粘结模型将表面能引入颗粒间的相互作用，其简化模型如图3所示。

图3　JKR粘结模型Fig.3　JKR contact model

1.4　休止角测定方法

休止角测定较广泛的方法主要有两种。一种方法是通过颗粒体高度H与底面直径D来计算休止角，即休止角

(1)

此种方法目前应用最为广泛[15-17]，但由于颗粒位置的不规则性，堆积体底面边缘较分散，不是连续的圆面，不便于度量；且堆积体顶端往往也不是规则的锥形，堆积体高度H难于确定。有研究者[18]提出将包含95%颗粒数目的圆作为堆积体底面边界圆，以该圆直径作为休止角底面直径。该方法在数值模拟中可行，但实际操作中难以找到边界圆。

另一种方法是选取堆积体轮廓中较为平直的一段曲线，以与该段曲线最接近的直线倾角作为该颗粒堆积体休止角。目前此方法主要用数字化图像分析(Digital image analysis，DIA)来实现，王云霞等[19]在对玉米种子休止角提取时用直线拟合堆积体轮廓线。但是由于所取曲线段的位置因人而异，此种方法结果受人为影响较大。FRACZEK等[20]在研究中亦指出了该方法的这一弊端。

观察蚯蚓基质堆积体和文献中其他散体物料堆积体的实际形状，发现大多数散体物料堆积体两端轮廓线近似为水平距离的凹函数，居中轮廓线近似为凸函数，而坡体中间段近似为直线。据此，拟采用高斯分布拟合堆积体轮廓线来获取堆积体休止角。主要思路是通过图像处理提取堆积体的轮廓点坐标后用高斯分布进行拟合，以拟合曲线拐点处的切线与水平轴的夹角为堆积体的休止角。其具体方法是，通过相机对颗粒堆积体进行拍照，得到图4a；通过Photoshop软件快速选择提取堆积体轮廓并对图片角度进行校正，得到图4b；通过Matlab软件将所获得的图像依次进行灰度处理(图4c)、二值化处理(图4d)，再通过Photoshop软件描边工具提取轮廓曲线(图4e)；最后将图像导入Origin软件通过图片文件数字化工具Digitizer获取轮廓点坐标，用高斯分布对其进行拟合(图4f)。图像处理的思路参考了FADAVI等[21]种子轮廓处理的方法。

图4　蚯蚓粪基质休止角的计算过程Fig.4　Calculation process of AoR of vermicomposting nursery substrate

高斯分布的一般方程为

(2)

式中y0、A、w、xc——常数

本文定义拟合方程拐点处切线与x轴所夹锐角为该散体的休止角。

对式(2)进行求导，分别求其一阶和二阶导数

(3)

(4)

曲线的拐点出现在二阶导数为零处，于是令f″(x)=0，式(4)有2个解

(5)

将式(5)代入式(3)，即可得到拐点处的斜率

(6)

则休止角的计算公式为

θ=arctan|Sl|

(7)

休止角测定试验在华中农业大学工科基地温室大棚进行。试验中，基质在遮荫下自然风干。不定时取样测定休止角，每次测定重复3次；同时另取样测定含水率。

1.5　参数标定试验设计

针对蚯蚓粪基质的研究目前尚不完善，其离散元模型的参数尤其缺乏。考虑到蚯蚓粪基质物理性质与土壤较为接近，且含水率变化时其泊松比、剪切模量等参数均会发生变化，本文主要参照文献[22-29]中各类土壤的参数值确定试验中各参数的取值或范围。表1给出了各待标定参数的高低水平。

其他参数[27]为钢泊松比0.3、钢剪切模量7.9×1010Pa、钢的密度7 865 kg/m3、重力加速度9.81 m/s2。

本文试验参数较多，参照文献[30-31]的试验设计，先进行Plackett-Burman试验，筛选出对结果影响显著的试验因素，再进行Box-Behnken试验，得到休止角和显著性参数之间的回归模型。通过测定不同含水率下蚯蚓粪基质的休止角，得到休止角随含水率变化的关系曲线。

表1　离散元法参数标定试验标定参数Tab.1　DEM calibration parameters need to be calibrated in experiment

2　结果与分析

2.1　Plackett-Burman试验

表2是使用Minitab设计的Plackett-Burman试验及试验所得休止角。

表3给出了由Plackett-Burman试验得出的参数显著性分析数据。根据试验结果，颗粒间静摩擦因数、颗粒间滚动摩擦因数、JKR表面能对休止角的影响是显著的，其他参数影响较小。

2.2　Box-Behnken试验

根据Plackett-Burman试验的结果，进行Box-Behnken试验，以得到休止角与3个显著参数之间的关系模型。

表4是Minitab设计的Box-Behnken试验及试验所得休止角。重点考察颗粒-颗粒静摩擦因数(A)、颗粒-颗粒滚动摩擦因数(B)和JKR表面能(C)这3个对休止角影响显著的参数。使用Minitab软件建立休止角与3个显著参数间的二次多项式回归方程

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

表3　Plackett-Burman试验参数显著性分析Tab.3　Analysis of significance of parameters in Plackett-Burman test

θ=10.39+3.66A+115.01B+30.447C-6.34A2-
56.31B2-3.477C2+2.17AB+1.632AC-22.25BC

(8)

表4　Box-Behnken试验设计及结果Tab.4　Design and results of Box-Behnken test

为0.998 1，均接近于1。表5是该模型方差分析的结果。模型的P值小于0.000 1，极为显著；失拟项的P值大于0.05，对结果不显著。模型拟合较好，可靠度较高，可以用来预测休止角。

表5　Box-Behnken试验休止角二次多项式回归模型方差分析Tab.5　Box-Behnken test AoR corner quadratic polynomial regression model variance analysis

为验证模型的有效性，在休止角测定试验中随机选取了2组试验结果(休止角分别为43.05°和36.38°(图5))，应用Minitab的响应优化器工具，在试验参数的取值范围内以休止角分别对回归模型进行寻找最优解，使用最优解进行离散元模拟，模拟所得休止角分别为43.71°和36.46°，均接近于实际休止角，差异分别为1.53%和0.22%。模拟结果与试验结果差别均较小，认为模型是有效的。最优解分别为颗粒间静摩擦因数为1.00，颗粒间滚动摩擦因数为0.11，JKR表面能为0.89 J/m2；颗粒间静摩擦因数为0.65，颗粒间滚动摩擦因数为0.12，JKR表面能为0.5 J/m2。

图5　蚯蚓粪基质堆积体仿真结果与实际堆积体图像Fig.5　Simulation results of vermicomposting nursery substrate accumulation and image of actual accumulation

2.3　含水率与休止角关系试验

不同含水率下基质的休止角变化如图6所示。蚯蚓粪基质的休止角随含水率的降低而降低，这与其他物料[17，32]的研究结果相似。目前，研究者[33-34]一般将散体含水率和休止角之间的关系用线性来描述。但是由蚯蚓粪基质休止角随含水率变化的散点图可以发现，随含水率降低，休止角减小的速率有放慢的趋势。本文分别用直线和指数函数对散点图进行拟合，得到拟合方程

θl=19.848 02Cm+30.808 68(R2=0.91)

(9)

(10)

指数方程式(10)拟合度更优，更符合两者之间的变化趋势。对于蚯蚓粪基质的休止角，可以测定其含水率，根据式(10)计算预测得到其休止角。

图6　蚯蚓粪基质含水率和休止角的关系曲线Fig.6　Relationship curves between vermicomposting nursery substrate moisture content and AoR

3　结论

(1)提出了一种散体颗粒休止角获取方法，即提取堆积体轮廓，用高斯分布对轮廓曲线进行拟合，以拟合曲线拐点的倾角作为休止角。试验表明，该方法拟合效果良好，为堆积体轮廓与蚯蚓粪基质相似的散体的休止角的测定提供了一种思路。

(2)将JKR粘结模型用于蚯蚓粪基质颗粒，采用离散元法对蚯蚓粪基质颗粒参数进行标定试验，筛选出对休止角影响显著的参数(即颗粒间静摩擦因数、颗粒间滚动摩擦因数、JKR表面能)，在此基础上对休止角试验数据进行二次多项式回归分析，建立了休止角与3个显著参数间的回归模型。经试验验证，模型结果与试验结果较为吻合。

(3)测定不同含水率下蚯蚓粪基质颗粒的休止角，建立了含水率与休止角之间的指数函数，对比发现较传统的直线拟合效果更优。由此函数及休止角与3个显著参数之间的回归模型，可以通过测定蚯蚓粪基质含水率、预测休止角，继而合理推定颗粒间静摩擦因数、颗粒间滚动摩擦因数、JKR表面能等其他参数。

1SANGWAN P, KAUSHIK C P, GARG V K. Vermicomposting of sugar industry waste (press mud) mixed with cow dung employing an epigeic earthwormEiseniafetida[J]. Waste Management & Research, 2010, 28(1): 71-75.

2HUANG K, XIA H. Role of earthworms’ mucus in vermicomposting system: biodegradation tests based on humification and microbial activity[J]. Science of the Total Environment, 2018, 610: 703-708.

3ANAND A, CURTIS J S, WASSGREN C R, et al. Experimental study of wet cohesive particles discharging from a rectangular hopper[J]. Industrial & Engineering Chemistry Research, 2015, 54(16): 4545-4551.

4BRIEND R, RADZISZEWSKI P, PASINI D. Virtual soil calibration for wheel-soil interaction simulations using the discrete-element method[J]. Canadian Aeronautics and Space Journal, 2011, 57(1): 59-64.

5刘凡一,张舰,李博, 等. 基于堆积试验的小麦离散元参数分析及标定[J]. 农业工程学报,2016,32(12):247-253.

LIU Fanyi, ZHANG Jian, LI Bo, et al. Calibration of parameters of wheat required in discrete element method simulation based on repose angle of particle heap[J]. Transactions of the CSAE, 2016, 32(12): 247-253. (in Chinese)

6JOHNSON K L, KENDALL K, ROBERTS A D. Surface energy and the contact of elastic solids[J]. Proceedings of the Royal Society of London Series A, Mathematical and Physical Sciences, 1971, 324(1558):301-313.

7CUNHA R N, SANTOS K G, LIMA R N, et al. Repose angle of monoparticles and binary mixture: an experimental and simulation study[J]. Powder Technology, 2016, 303: 203-211.

8张帅, 单忠德, 张杰. 基于离散元方法的型砂流动性仿真研究[J]. 铸造技术, 2016, 37(2):288-291.

ZHANG Shuai, SHAN Zhongde, ZHANG Jie. Simulation of self-hardening resin sand mobility based on DEM[J]. Foundry Technology, 2016, 37(2):288-291. (in Chinese)

9胡红. 玉米行间定点扎穴深施追肥机设计与研究[D]. 北京：中国农业大学, 2017.

HU Hong. Design and research of targeted hole-pricking and deep-application fertilizer between maize rows [D]. Beijing: China Agricultural University, 2017. (in Chinese)

10GUO Z, CHEN X, LIU H, et al. Theoretical and experimental investigation on angle of repose of biomass-coal blends[J]. Fuel, 2014, 116(1): 131-139.

11GHAZAVI M, HOSSEINI M, MOLLANOURI M. A comparison between angle of repose and friction angle of sand[C]∥The 12th International Conference for International Association for Computer Methods and Advances in Geomechanics (IACMAG), 2008.

12FIELKE J M, UCGUL M, SAUNDERS C. Discrete element modeling of soil-implement interaction considering soil plasticity, cohesion and adhesion[C]∥2013 ASABE International Meeting, 2013.

13FRIEDMAN S P, ROBINSON D A. Particle shape characterization using angle of repose measurements for predicting the effective permittivity and electrical conductivity of saturated granular media[J]. Water Resources Research, 2002, 38(11): 11-18.

14XIAO X, TAN Y, ZHANG H, et al. Experimental and DEM studies on the particle mixing performance in rotating drums: effect of area ratio[J]. Powder Technology, 2017, 314: 182-194.

15PROBST K, AMBROSE K, ILELEJI K. The effect of moisture content on the grinding performance of corn and corncobs by hammer milling[J]. Transactions of the ASABE, 2013, 56(3):1025-1033.

16UNAL H, ISIK E, IZLI N, et al. Geometric and mechanical properties of mung bean (VignaradiataL.) grain: effect of moisture[J]. International Journal of Food Properties, 2008, 11(3): 585-599.

17ZAALOUK A K, ZABADY F I. Effect of moisture content on angle of repose and friction coefficient of wheat grain[J]. Misr Journal of Agricultural Engineering, 2009, 26(1): 418-427.

18韩燕龙, 贾富国, 唐玉荣, 等. 颗粒滚动摩擦系数对堆积特性的影响[J]. 物理学报, 2014, 63(17):173-179.

HAN Yanlong, JIA Fuguo, TANG Yurong, et al. Influence of granular coefficient of rolling friction on accumulation characteristics[J]. Acta Physica Sinica, 2014, 63(17):173-179. (in Chinese)

19王云霞, 梁志杰, 张东兴, 等. 基于离散元的玉米种子颗粒模型种间接触参数标定[J]. 农业工程学报,2016,32(22):36-42.

WANG Yunxia, LIANG Zhijie, ZHANG Dongxing, et al. Calibration method of contact characteristic parameters for corn seeds based on EDEM[J]. Transactions of the CSAE, 2016, 32(22): 36-42. (in Chinese)

21FADAVI A, MIRZABE A H, MANSOURI A. Moisture-dependent physical properties of Plantain (PlantagomajorL.) seeds by image processing analysis[J]. Agricultural Engineering International: CIGR Journal, 2015, 17(3): 353-363.

22张锐, 韩佃雷, 吉巧丽, 等. 离散元模拟中沙土参数标定方法研究[J/OL]. 农业机械学报, 2017, 48(3):49-56.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170306&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2017.03.006.

ZHANG Rui, HAN Dianlei, JI Qiaoli, et al. Calibration methods of sandy soil parameters in simulation of discrete element method[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3):49-56. (in Chinese)

23武涛, 黄伟凤, 陈学深,等. 考虑颗粒间黏结力的黏性土壤离散元模型参数标定[J]. 华南农业大学学报, 2017, 38(3):93-98.

WU Tao, HUANG Weifeng, CHEN Xueshen, et al. Calibration of discrete element model parameters for cohesive soil considering the cohesion between particles[J]. Journal of South China Agricultural University, 2017, 38(3):93-98. (in Chinese)

24潘世强. 基于离散元法的芯铧式开沟器优化设计与试验研究[D]. 长春：吉林大学, 2015.

PAN Shiqiang. Research on the optimization design and the experiment of the core ploughshare furrow opener based on the discrete element method [D]. Changchun: Jilin University, 2015. (in Chinese)

25方会敏, 姬长英, AHMED A T, 等. 秸秆-土壤-旋耕刀系统中秸秆位移仿真分析[J/OL]. 农业机械学报, 2016, 47(1):60-67. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160109&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2016.01.009.

FANG Huimin, JI Changying, AHMED A T, et al. Simulation analysis of straw movement in straw-soil-rotary blade system[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(1):60-67. (in Chinese)

26UCGUL M, FIELKE J M, SAUNDERS C. 3D DEM tillage simulation: validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil[J]. Soil and Tillage Research, 2014, 144(4): 220-227.

27UCGUL M, FIELKE J M, SAUNDERS C. Three-dimensional discrete element modelling (DEM) of tillage: accounting for soil cohesion and adhesion[J]. Biosystems Engineering, 2015, 129: 298-306.

28邵明安，王全九，黄明斌. 土壤物理学[M]. 北京：高等教育出版社， 2006.

29依艳丽. 土壤物理研究法[M]. 北京：北京大学出版社, 2009.

30ABDENACER M, KAHINA B I, AÏCHA N, et al. Sequential optimization approach for enhanced production of glutamic acid fromCorynebacteriumglutamicum2262 using date juice[J]. Biotechnology and Bioprocess Engineering, 2012, 17(4): 795-803.

31JOB J, SUKUMARAN R K, JAYACHANDRAN K. Production of a highly glucose tolerant β-glucosidase byPaecilomycesvariotiiMG3: optimization of fermentation conditions using Plackett-Burman and Box-Behnken experimental designs[J]. World Journal of Microbiology and Biotechnology, 2010, 26(8): 1385-1391.

32KALKAN F, KARA M. Handling, frictional and technological properties of wheat as affected by moisture content and cultivar[J]. Powder Technology, 2011, 213(1): 116-122.

33UNAL H, ISIK E, IZLI N, et al. Geometric and mechanical properties of mung bean (VignaradiataL.) grain: effect of moisture[J]. International Journal of Food Properties, 2008, 11(3): 585-599.

34SOLOGUBIK C A, CAMPAONE L A, PAGANO A M, et al. Effect of moisture content on some physical properties of barley[J]. Industrial Crops and Products, 2013, 43(5): 762-767.