不同含水率羊粪离散元参数通用标定方法研究

朱新华 伏胜康 李旭东 魏玉强 赵 伟

(1.西北农林科技大学机械与电子工程学院， 陕西杨凌 712100； 2.四川航天职业技术学院， 成都 610100)

0 引言

粪肥施肥体量大、工费高，而高效的施肥机械较为缺乏。目前，借助离散元法分析肥料颗粒间及肥料与装置间复杂的相互作用已成为施肥机械设计、优化的重要手段，可节约研发成本、缩短周期，提高装置性能[1-4]。离散元法的应用以准确的离散元参数为前提[5-6]，快速准确标定离散元仿真参数是其在农业工程领域应用的关键[7]。羊粪是我国北方来源广、抗生素等污染小的有机肥源，其有机质含量高达40%～65%[8]，高于牛粪、马粪、猪粪和鸡鸭鹅等粪肥[9]，在实际生产中受到重视，但腐熟羊粪的离散元参数标定尚未见报道。当前，粪肥离散元参数主要采用虚拟标定试验获得[10-11]，但这些研究所构建的标定方法存在通用性和实用性不足的问题。

国外学者以液态有机肥流体特性的研究为主，而对固态粘湿有机肥离散元参数标定方法的研究尚未见报道[7]。国内学者基于堆积角试验对有机肥离散元参数标定的研究较多[12-13]，这些研究表明运用堆积角试验标定粪肥离散元参数能得到较为准确的结果。但上述研究均针对特定含水率下的粪肥进行研究，而实际粪肥的含水率范围较大，因此，该方法的通用性不足[7]。物料的离散元参数受含水率影响较大，特定含水率下标定的离散元参数不能真实反映其它含水率的粪肥与机械之间的作用规律。

基于堆积角试验，罗帅等[14]对不同含水率的蚯蚓粪基质离散元参数进行标定，建立了堆积角与离散元参数的关系模型；在含水率40%～85%范围内，王黎明等[15]探索了通过堆积角快速标定猪粪离散元参数的方法。尽管上述研究考虑了含水率对粪肥离散元参数的影响，提高了堆积角-离散元参数模型和标定方法的通用性，但存在以下不足：粪肥堆积角测量的偶然误差大，必须多次测量，且粪肥堆积操作不便，以堆积角-离散元参数模型预测粪肥离散元参数，实用性不足；研究中仅建立单一的肥料颗粒模型，与实际肥料的粒度分布差异较大[7]，且对所建模型的准确性缺少工程性验证(即用另外一种试验方法验证)，所得离散元参数可靠性不足。

因此，针对不同含水率下畜禽粪肥的离散元参数快速标定方法通用性和实用性不足的问题，本文以含水率为31.91%～58.48%的腐熟羊粪为研究对象，通过圆筒提升法[16]建立粪肥的含水率-堆积角模型；以堆积角为响应值，通过Plackett-Burman试验、爬坡试验和Box-Behnken试验建立堆积角-离散元参数模型；进而经推导建立羊粪含水率-离散元参数模型，并运用抽板法[16]进行工程性验证，以期构建一种通用性和实用性强的羊粪离散元参数标定方法。

1 材料与方法

1.1 试验材料与粒度分布

腐熟羊粪取自陕西省杨凌区西北农林科技大学第三试验站，原始含水率46.81%，密度0.47 g/cm3。在腐熟粪肥通常含水率范围(30%～60%)[17-19]内通过自然风干法和添加纯净水的方法配制5个含水率梯度，计算方法为[14]

(1)

式中ms——配制样品的加水质量，kg

m0——样品质量，kg

ω0——原始含水率，%

ω1——配制样品的目标含水率，%

每次堆积角测量的同时，均取一定量羊粪样品通过烘干法测量其含水率，测得供堆积角试验的5个含水率依次为：31.91%、39.33%、46.81%、53.32%和58.48%，如图1所示。取300 g原始含水率下的羊粪样本，利用振动筛(ZFJ-Ⅱ型，浙江省上虞市大亨桥化验仪器厂)对其筛分，分析其粒径分布(表1)。

表1 羊粪粒径分布Tab.1 Particle size distribution of goat manure

图1 不同含水率的羊粪试样Fig.1 Samples of goat manure with different moisture contents

1.2 羊粪堆积角试验及含水率-堆积角模型构建

采用圆筒提升法测定羊粪堆积角[16]。图2为试验装置，由电动式拉力试验机(HSV型)、圆筒(内径100 mm、长300 mm)和平台组成[16,20]。圆筒和平台均为钢材，圆筒直立于平台上。试验前在圆筒内填满羊粪，试验时通过拉力试验机以20 mm/s的速度提升圆筒[16]，使羊粪自圆筒下端漏出，形成肥堆，如图3所示。用数显倾角仪从肥堆的4个不同方向测量堆积角。以相同的方法依次对5种含水率羊粪进行堆积角测量，每组试验重复5次，结果取平均值。

图2 圆筒提升法堆积角测试装置Fig.2 Cylinder lifting resting angle tester1.拉力试验机 2.圆筒 3.平台

图3 堆积角测量Fig.3 Measurement of resting angle

基于测量结果，建立羊粪含水率-堆积角模型，并分析其相关性。

1.3 羊粪离散元模型构建

(1)接触模型。羊粪流动性差，颗粒间存在粘附力，故采用Hertz-Mindlin with JKR模型[21-22](简称JKR模型)。该接触模型内置于EDEM(2018，EDEM Solutions)软件中，是一种能体现颗粒粘结性的接触模型，模型基于Johnson-Kendall-Roberts理论建立，在Hertz-Mindlin接触模型的基础上考虑了接触区域内范德华力对颗粒运动的影响，适用于模拟干粉状物料或含湿颗粒，并通过数值surface energy表征颗粒间的粘附力，可较好地模拟羊粪。

(2)颗粒粒径分布及材料本征参数。根据测定的羊粪粒径分布，同时考虑合理简化模型，在EDEM仿真中以球形颗粒模拟羊粪，颗粒粒径设置4种，分布如图4所示。

图4 羊粪颗粒模型Fig.4 Particle models of goat manure

羊粪和钢材的本征参数如表2所示。

表2 羊粪和钢材的本征参数Tab.2 Intrinsic parameters of goat manure and steel material

(3)接触参数和JKR模型参数。结合EDEM通用颗粒材料数据库(Generic EDEM material model database，GEMM数据库)，根据羊粪、钢材的本征参数及文献[10-11,13-14,16]的推荐值，确定接触参数及JKR模型参数取值范围如表3所示。

表3 羊粪与钢材接触参数和JKR模型参数范围Tab.3 Contact parameters and JKR model parameter range of goat manure and steel material

1.4 羊粪离散元参数标定及堆积角-离散元参数模型构建

1.4.1羊粪堆积角仿真方法

仿真过程设置与物理试验一致。设置圆筒提升速度20 mm/s，生成羊粪颗粒总质量3 kg，生成速度0.2 kg/s，仿真固定时间步长为Rayleigh时间步长的20%，为7.4×10-6s，数据保存时间间隔为0.05 s，网格大小为1.5 mm，仿真总时长为15 s。仿真几何模型及堆积角测量如图5所示，运用EDEM后处理Clipping模块在坐标为(1，0，0)和(0，1，0)羊粪堆中垂面处设置Clipping Planes Group截面，运用Tools模块中的量角器Protractor，依次选定Clipping Planes Group截面上羊粪堆单侧的斜边、夹角和水平底面线上的3个羊粪颗粒作为堆积角R1的3个点，获得堆积角数值。

图5 圆筒提升法仿真Fig.5 Simulation of cylinder lifting method

1.4.2羊粪堆积角-离散元参数模型构建

应用Design-Expert 8.0(STAT-EASE)软件设计Plackett-Burman(PBD)试验，以堆积角为响应值，从羊粪的密度、泊松比、剪切模量以及表3所示的接触参数和JKR模型参数共10个参数中筛选对堆积角影响显著的离散元参数。各参数取高、中、低3个水平，分别以编码1、0和-1表示，如表4所示，试验方案及结果见2.2节。

表4 Plackett-Burman试验因素编码Tab.4 Factors and levels of Plackett-Burman test

为缩小参数范围并确定最优区间，针对显著性参数设计爬坡试验。试验时，非显著性参数取中间水平值，显著性参数按照设计的步长逐步增大，将仿真堆积角与中间水平含水率(46.81%)羊粪物理试验堆积角(41.00°)相对误差最小时的参数作为显著性参数最优区间的中间水平值，以确定最优取值区间。爬坡试验方案及结果见2.3节。

基于爬坡试验确定的离散元参数最优区间设计Box-Behnken试验，并对试验结果进行模型优化，以建立堆积角与显著性参数的关系模型。采用3个中心点进行误差估计，共15组试验，每组试验重复3次，取堆积角平均值作为试验结果。试验中各显著性参数取爬坡试验优化后的高、中、低3个水平，分别以编码1、0和-1表示，非显著性参数取值同爬坡试验一致，试验方案及结果见2.4节。

以5个含水率下羊粪物理试验堆积角为目标，对Box-Behnken试验得到的堆积角-离散元参数模型求解寻优得到离散元参数的最优组合。基于离散元参数最优组合(非显著性参数取值同爬坡试验一致)，运用EDEM进行仿真并测量其堆积角，对比仿真堆积角与物理试验堆积角以验证堆积角-离散元参数模型的准确性。

1.5 羊粪含水率-离散元参数模型建立及工程性验证试验

通过堆积角-离散元参数模型能够推测离散元参数，但粪肥堆积角测量误差大、过程繁琐，而含水率测量精度高、简便。因此，基于羊粪含水率-堆积角模型和羊粪堆积角-离散元参数模型推导，建立羊粪含水率-离散元参数模型，以构建通用性和实用性好、准确度高的离散元参数标定方法。

为评价基于圆筒提升法建立的含水率-离散元参数模型的可靠性，通过另外一种堆积角试验方法——抽板法[17, 20]对模型进行工程性验证。验证过程中，仿真与物理试验方法一致，图6为基于抽板法的堆积角仿真与物理验证试验。

图6 基于抽板法的堆积角仿真与物理验证试验Fig.6 Resting angle simulation and physical verification test based on side plate lifting method

在含水率31.91%～58.48%范围内随机配制3个不同含水率的羊粪，并分别测量其物理试验堆积角；根据推导出的含水率-离散元参数模型求解，得到3个含水率下的羊粪离散元参数最优组合；以离散元参数最优组合进行仿真，分别得到仿真堆积角。对比物理试验堆积角与仿真堆积角以验证羊粪含水率-离散元参数模型的准确性。每组试验重复5次，结果取平均值。

2 结果与分析

2.1 含水率-堆积角模型构建物理试验测定的不同含水率羊粪样本的堆积角如表5所示。

表5 不同含水率羊粪的堆积角Tab.5 Resting angle of goat manure with different moisture contents

根据表5，以羊粪样本的含水率与堆积角拟合多项式，得到含水率-堆积角模型

R1=0.005 5x2-0.202 2x+38.329

(2)

式中x——羊粪含水率，%

模型相关系数达0.999 9，拟合曲线如图7所示。在含水率31.91%～58.48%范围内，羊粪堆积角随含水率升高而增大。羊粪属于散粒体物料，随着含水率的增大，羊粪颗粒间的粘度增大，流动性变差，故其堆积角变大。

图7 羊粪含水率与堆积角关系曲线Fig.7 Relationship curve between moisture content and resting angle of goat manure

2.2 显著性参数筛选结果与分析

基于Design-Expert 8.0软件设计PBD试验，试验方案与结果如表6所示，A～J为因素编码值。

表6 Plackett-Burman试验方案与结果Tab.6 Scheme and results of Plackett-Burman test

表7为PBD试验结果显著性分析。

表7 Plackett-Burman试验结果显著性分析Tab.7 Significance analysis of Plackett-Burman test results

根据表7，JKR表面能J(P<0.01)对羊粪堆积角影响极显著，颗粒间滚动摩擦因数F和静摩擦因数E(P<0.05)对羊粪堆积角影响显著，且3个参数对堆积角的影响均为正效应。分析其原因为：①随着JKR表面能增大，羊粪颗粒间的粘附性能增强，从而影响羊粪颗粒的流动性，因此，堆积角增大。②由于颗粒流动表现为颗粒与颗粒间的滚动和滑动，因此，影响颗粒流动性的主要参数为颗粒间滚动摩擦因数和静摩擦因数。当滚动摩擦因数和静摩擦因数增大时，颗粒间的滚动和滑动受阻，流动性变差，堆积角增大。③模型中的羊粪由4种不同粒径的球型颗粒组成，而球体模型的运动方式以滚动为主，因此，滚动摩擦因数比静摩擦因数对堆积角的影响更显著。

为方便后续试验，在爬坡试验、Box-Behnken试验和验证试验中只考虑3个显著性参数，不显著因素取中间水平值[13]。

2.3 离散元参数最优区间确定

爬坡试验方案与结果如表8所示，可知，当JKR表面能为0.05 J/m2、羊粪颗粒间滚动摩擦因数为0.15、静摩擦因数为0.30时，仿真堆积角与物理试验堆积角相对误差最小，为2%，即该组参数为显著性参数最优区间的中间水平值，对应最优区间分别为[0，0.10]、[0，0.30]、[0.10，0.50]。

表8 爬坡试验方案与结果Tab.8 Scheme and results of climbing test

2.4 堆积角-离散元参数模型建立

表9 Box-Behnken试验方案与结果Tab.9 Scheme and results of Box-Behnken test

表10 Box-Behnken试验模型方差分析Tab.10 ANOVA of Box-Behnken model

根据试验结果及方差分析，E、F、J、E2、F2均对羊粪堆积角影响极显著，EJ对羊粪堆积角影响显著。建立的羊粪堆积角-离散元参数模型为

R1=37.04+4.96E+7.50F+11.91J-
2.88EJ+5.51E2+6.82F2

(3)

2.5 堆积角-离散元参数模型验证

以5个含水率羊粪样本物理试验堆积角为目标求解寻优，分别获得5组离散元参数(JKR表面能、颗粒间滚动摩擦因数和静摩擦因数)的最优组合，如表11所示。

表11 不同含水率下羊粪离散元参数最优值及验证试验结果Tab.11 Optimal value of goat manure discrete element parameters and verification test results under different moisture contents

当含水率为58.48%时，仿真和物理试验获得的羊粪堆积角如图8所示。由表11可知，由堆积角-离散元参数模型确定3个显著性离散元参数的最优组合后，各含水率羊粪仿真堆积角与物理试验堆积角的相对误差小于等于2.42%，表明建立的模型可用于离散元参数确定。

图8 仿真与物理试验获得的堆积角Fig.8 Resting angle obtained by simulation and physical test

2.6 含水率-离散元参数模型构建及验证

为准确、快速标定离散元参数，基于含水率-堆积角模型和堆积角-离散元参数模型，由式(2)、(3)推导获得含水率与离散元参数(JKR表面能、颗粒间滚动摩擦因数和静摩擦因数)的模型，即含水率-离散元参数模型

4.96E+7.50F+11.91J-2.88EJ+5.51E2+
6.82F2=0.005 5x2-0.202 2x+1.289

(4)

根据式(4)，可通过含水率直接得到离散元参数组合目标值。即借助Design-Expert 8.0软件中的模型目标值寻优求解，得到离散元参数组合。

基于羊粪含水率-离散元参数模型，获得3组随机含水率下的羊粪离散元参数最优组合，通过抽板法物理试验与仿真分别获得羊粪物理试验堆积角与仿真堆积角，如图9所示。每组试验重复5次，结果取平均值，测量结果如表12所示。

图9 抽板法试验与仿真对比Fig.9 Comparisons of test and simulation with side plate lifting method

表12 抽板法验证试验结果Tab.12 Verification test results of side plate lifting method

由表12可知，3组试验中仿真堆积角与物理试验堆积角的相对误差小于等于5.37%。表明建立的含水率-离散元参数模型可靠。

羊粪含水率-离散元参数模型的相对误差略大于堆积角-离散元参数模型，这是由于模型验证试验中采用了不同于圆筒提升法的抽板法。抽板法在此作为工程验证，能进一步反映模型的可靠性。两个模型的相对误差均较小，说明建立的两种模型均可用于离散元参数的预测。

建立的羊粪含水率-堆积角模型与罗帅等[14]研究的蚯蚓粪基质含水率与堆积角关系变化趋势一致，且筛选得到的显著性参数相同，分别为JKR表面能、颗粒间滚动摩擦因数和静摩擦因数。建立的堆积角-离散元参数模型相对误差小于等于2.42%，均小于含水率25%～65%的蚯蚓粪基质[23]以及含水率40%～85%的猪粪[15]堆积角预测模型误差。建立的含水率-离散元参数模型通过工程性验证，相对误差小于等于5.37%，表明模型可靠。由于含水率测量比堆积角测量准确性高、易操作，考虑到离散元参数预测模型的实际应用，建立的含水率-离散元参数模型比堆积角-离散元参数模型更具实用性。两种模型均能适应有机肥料不同含水率下离散元参数标定，具有通用性。

3 结论

(1)基于圆筒提升法物理试验建立了31.91%～58.48%范围内羊粪含水率-堆积角模型，模型相关系数为0.999 9。

(2)通过Plackett-Burman试验和爬坡试验确定了含水率31.91%～58.48%的羊粪离散元显著性参数及其最优区间：JKR表面能为0～0.10 J/m2、羊粪颗粒间滚动摩擦因数为0～0.30、静摩擦因数为0.10～0.50；通过Box-Behnken试验建立了羊粪堆积角-离散元参数模型，模型P<0.000 1，相对误差小于等于2.42%。

(3)基于含水率-堆积角模型和堆积角-离散元参数模型推导，建立了含水率-离散元参数模型；采用抽板法试验验证了羊粪含水率-离散元参数模型的可靠性，相对误差小于等于5.37%。

(4)基于建立的含水率-离散元参数模型和堆积角-离散元参数模型，均可求解出特定含水率下羊粪离散元显著性参数的最优组合，对不同含水率羊粪参数标定具有通用性和实用性。