玉米秸秆饲料除尘筛出物离散元数值模拟参数标定

肖子卿 田海清* 张 涛 王 迪 盛 越 李大鹏 刘 飞 周建成

(1.内蒙古农业大学 机电工程学院，呼和浩特 010018；2.内蒙古领异农科技发展有限公司，呼和浩特 010200)

我国秸秆资源丰富，将秸秆资源转化为饲草料，不仅能解决畜牧养殖饲草料缺乏的问题，而且还能改善生态环境。现有秸秆饲料收获机可对田间秸秆收获并破碎为直接可饲喂的段状或丝状，但收获机捡拾的秸秆中掺杂大量尘土，从而导致收获的秸秆饲料清洁程度较低，适口性较差。因此对收获后的秸秆饲料进行除尘提高清洁度，对提高饲料利用率意义重大。

目前带有除尘装置的饲料收获机市场所占比例逐年增加。如张开飞等设计的小麦收获机复合除尘装置，实现了割台静电除尘和喷淋除尘相结合的复合除尘；杜洋研制的秸秆除尘设备，实现了对秸秆的净化，进而保障秸秆二次加工产物的质量；李嘉诚等设计的新型无尘脱粒机，采用空气幕防尘技术从源头对秸秆除尘；类似的还有秸秆收获机螺旋输送除尘装置、具有除尘功能的滚筒筛和带有除尘结构的饲料粉碎机等。近几年结合收获、除尘、打捆为一体的饲料收获机也开始兴起，市场上出现的有顺邦4JL-2.2型秸秆饲料收获机和绿萌机械9YFG-2.2P型秸秆揉碎除尘方草捆打捆机等。带有除尘装置的饲料收获机虽逐渐增加，但对除尘机理的研究还较为少见。研究除尘过程中秸秆饲料与尘土的运动分离机理，对优化设计除尘装置，提升除尘效率有积极意义。

本研究前期对收获机收获后的饲料采用标准筛筛分时发现，所筛出的杂质为尘土与秸秆碎末的混合物。通过试验方法分析秸秆饲料与尘土运动分离机理较为困难，而离散元法是模拟混合物内部运动状态的有效方法，国外相关研究在运用此方法对固体颗粒、粉末等混合物的内部运动及参数标定中已取得较多成果。国内也有研究将此方法应用在混合物中存在秸秆时的运动状态分析及参数标定中，如张涛等为研究玉米秸秆与揉碎机的相互作用及运动规律，标定了玉米秸秆与揉碎机锤片、秸秆间的接触参数。王韦韦等标定了玉米秸秆粉料与不锈钢板、秸秆粉料间的接触参数；马彦华等标定了苜蓿秸秆与45钢、苜蓿秸秆间的接触参数。综上，秸秆已有较为完整的物理参数与数值模拟参数，但对于尘土与秸秆碎末之间的物理参数与数值模拟参数的研究仍较为少见。

本研究拟采用物理试验与数值模拟相结合的方法对玉米秸秆饲料的筛出物(尘土与秸秆碎末)进行离散元数值模拟参数标定。利用EDEM软件对尘土与秸秆碎末运动情况进行数值模拟，以休止角为响应值，分别通过Plackett-Burman试验、最陡爬坡试验和Box-Behnken试验确定较优的尘土与秸秆碎末的数值模拟参数值，以期为尘土与秸秆运动分离机理分析及除尘装置的优化提供数据支持。

1 材料及方法

1.1 试验材料

2020年11月进行秸秆收获作业试验，试验地点位于内蒙古自治区呼和浩特市土默特左旗大圪贲村(121°58′ E，36°15′ N)试验田，试验所用秸秆饲料和土壤样品均取自收获作业后的试验田。根据GB/T 5917.1—2008 《饲料粉碎粒度测定两层筛筛分法》选用筛孔直径为3和7 mm的标准筛对秸秆饲料进行筛分，筛出物包括尘土和秸秆碎末。利用游标卡尺测得秸秆碎末平均直径为2 mm，平均长度5 mm。利用标准土壤筛对尘土进行筛分，尘土平均粒径为0.5 mm。

1.2 试验方法

运用休止角物理试验与数值模拟相结合的方法对玉米秸秆饲料的筛出物(尘土与秸秆碎末)进行离散元数值模拟参数标定，利用休止角测定仪进行尘土-秸秆碎末休止角测定试验，测量尘土与秸秆碎末的锥体高度，计算得出休止角。利用EDEM软件进行尘土-秸秆碎末休止角数值模拟试验，分别通过Plackett-Burman试验、最陡爬坡试验和Box-Behnken试验确定最优参数组合。应用最优参数组合进行休止角和标准筛除尘率数值模拟，与休止角和标准筛除尘率物理试验结果相对比，验证标定参数组合的准确性。

1.3 休止角测定试验

休止角测定采用FT-104B型休止角测定仪，漏斗上方进料口直径为135 mm，下方出料口直径20 mm，高度140 mm。透明环形容器位于漏斗下方50 mm处，直径为100 mm，高度25 mm。将仪器平置于操作平台后，将基板调整水平，将透明环形容器置于基板中心位置，使容器的中心与漏斗轴线一致。参考文献中对休止角测定方法，将筛选出的尘土和秸秆碎末匀速倒入仪器中的漏斗，待其流出停止后120 s，测量锥体高度。试验重复5次取平均值，锥体平均高度为69.5 mm，计算得到秸秆碎末和尘土颗粒堆物理试验休止角

α

为41.669°。图1为休止角测定装置及试验材料。

1.漏斗；2.秸秆碎末与土；3.透明环形容器；4.基板1.Funnel; 2.Straw debris and soil;3.Transparent annular container; 4.Substrate图1 休止角测定装置及试验材料Fig.1 Angle of repose measuring device and test material

物理试验休止角(

α

)计算公式为：

(1)

式中：

h

为锥尖到基板的高度，mm；

L

为锥体底端到基板的高度，mm；

R

为透明环形容器半径，mm。

1.4 接触参数测定

秸秆饲料的筛出物中有尘土与秸秆碎末，秸秆碎末中包含秸秆外皮和秸秆內瓤，所以需要测定的接触参数有尘土-外皮滚动摩擦因数、尘土-外皮静摩擦因数、尘土-內瓤滚动摩擦因数和尘土-內瓤静摩擦因数，测定结果作为数值模拟输入参数选择依据。

物理试验选用的秸秆均为饲料收获机收获后残留于田间的漏收秸秆，本研究前期已测得其外皮密度为508 kg/m、內瓤密度40 kg/m、总密度110.9 kg/m、泊松比0.3、剪切模量64 MPa。与秸秆接触的材料为钢，密度为7 865 kg/m、泊松比0.3、剪切模量79 000 MPa。

1

.

4

.

1

静摩擦因数测定试验使用CNY-1型斜面仪对尘土-外皮和尘土-內瓤静摩擦因数进行测定。斜面上的物料受3个力作用，分别为重力

G

、支持力

F

和摩擦力

F

。根据力的正交分解得：

(2)

得出：

μ

=tan

θ

(3)

式中：

G

为物料所受重力，N；

m

为物料质量，kg；

g

为重力加速度，m/s；

F

为物料所受支持力，N；

F

为物料所受摩擦力，N；

θ

为斜面仪的倾斜角度,(°)；

μ

为摩擦因数。将取自试验田的土壤样品装于直径90 mm、高10 mm的透明圆柱容器内，平整其表面并固定于斜面仪上。将秸秆切割为30 mm的规则圆柱体，为防止在测定过程中秸秆发生滚动，将2个规则秸秆段粘在一起并置于土样之上。调整斜面仪使其指针位于零刻度处，将斜面仪缓慢逆时针转动直至秸秆在土样上匀速滑动，记录斜面仪倾斜角度

θ

，根据式(3)求得静摩擦因数。重复10次试验，得出秸秆在土样上发生滑动的角度为21.801°～34.992°，尘土-外皮的静摩擦因数为0.4～0.7。去除秸秆外皮，将內瓤切割为长、宽、高均为15 mm的规则正方体，将內瓤置于土样之上，缓慢逆时针转动斜面仪直至內瓤方块在土样上缓慢滑动，记录斜面仪倾斜角度

θ

。重复10次试验，得出秸秆内瓤在土样上发生滑动的角度为30.964°～41.987°，尘土-内瓤的静摩擦因数为0.6～0.9。

1

.

4

.

2

滚动摩擦因数测定试验测定滚动摩擦因数的仪器也为CNY-1型斜面仪，测定方法与静摩擦因数测定方法类似。将切割好的单个秸秆段置于土样之上，调整斜面仪使其指针位于零刻度处，缓慢逆时针转动斜面仪，当秸秆开始滚动时记录斜面仪倾斜角度

θ

。重复10次试验，得出秸秆在土样上发生滚动的角度范围为11.310°～26.565°，尘土-外皮的滚动摩擦因数为0.2～0.5。去除秸秆外皮，将圆柱形的秸秆內瓤置于土样之上，将指针转动到零刻度处时开始缓慢逆时针转动斜面仪，当秸秆內瓤开始滚动时停止转动，记录斜面仪倾斜角度

θ

。重复10次试验，得出秸秆內瓤在土样上发生滚动的角度为16.699°～30.964°，尘土-內瓤的滚动摩擦因数为0.3～0.6。

为了保证数值模拟的准确性，根据物理试验测得的数据，得出测定参数的取值范围，为后续的数值模拟输入参数取值提供依据，最后通过数值模拟确定参数的最优取值。

2 秸秆与尘土离散元模型的建立及参数标定

2.1 离散元模型的建立

秸秆饲料经筛分后留下的杂质包括秸秆碎末和尘土，秸秆碎末包含秸秆外皮和秸秆內瓤，外皮和內瓤的长度与直径差异性不大。因此在EDEM软件中建立数值模拟模型时将秸秆外皮和秸秆內瓤的数值模拟模型都设置为长5 mm、直径2 mm的圆柱体，尘土的数值模拟模型设置为直径0.5 mm的球体(图2)。

图2 秸秆外皮、內瓤(a)和尘土(b)数值模拟标定模型Fig.2 Calibration model for numerical simulation of straw skin and pulp (a) and soil particle (b)

导入休止角测定仪模型，模型漏斗进料口直径为135 mm，出料口直径20 mm，下方圆盘直径50 mm、高度25 mm(图3)。设置颗粒模型参数，结合GEMM数据库与相关文献，得出部分参数取值见表1。

1.漏斗；2.秸秆碎末与土；3.圆盘1.Funnel; 2.Straw debris and soil; 3.Disk图3 休止角数值模拟模型Fig.3 Numerical simulation model of angle of repose

表1 尘土与秸秆碎末休止角数值模拟试验部分参数取值
Table 1 Some parameter values of dust and straw debris angle of repose numerical simulation test

参数Parameter数值Numericalvalue参数Parameter数值Numericalvalue秸秆外皮-外皮静摩擦因数Straw skin-straw skin static friction factor0.142秸秆外皮-外皮碰撞恢复系数Straw skin-straw skin restitution coefficient0.485秸秆外皮-外皮滚动摩擦因数Straw skin-straw skin rolling friction factor0.078秸秆外皮-內瓤碰撞恢复系数Straw skin-straw pulp restitution coefficient0.263秸秆外皮-內瓤静摩擦因数Straw skin-straw pulp static friction factor0.495秸秆内瓤-內瓤碰撞恢复系数Straw pulp-straw pulp restitution coefficient0.348秸秆外皮-內瓤滚动摩擦因数Straw skin-straw pulp rolling friction factor0.166秸秆外皮-钢碰撞恢复系数Straw skin-steel restitution coefficient0.663秸秆内瓤-內瓤静摩擦因数Straw pulp-straw pulp static friction factor0.427秸秆内瓤-钢碰撞恢复系数Straw pulp-steel restitution coefficient0.325秸秆内瓤-內瓤滚动摩擦因数Straw pulp-straw pulp rolling friction factor0.144秸秆外皮-尘土碰撞恢复系数Straw skin-soil restitution coefficient0.213秸秆外皮-钢静摩擦因数Straw skin-steel static friction factor0.226秸秆内瓤-尘土碰撞恢复系数Straw pulp-soil restitution coefficient0.135秸秆外皮-钢滚动摩擦因数Straw skin-steel rolling friction factor0.119尘土泊松比Soil poisson’s ratio0.4 秸秆内瓤-钢静摩擦因数Straw pulp-steel static friction factor0.375尘土密度ρ/(kg/m3)Soil density2 600秸秆内瓤-钢滚动摩擦因数Straw pulp-steel rolling friction factor0.150

本研究需标定的试验参数取值分别为尘土-钢滚动摩擦因数、尘土-钢静摩擦因数、尘土-外皮滚动摩擦因数、尘土-外皮静摩擦因数、尘土-内瓤滚动摩擦因数、尘土-内瓤静摩擦因数、尘土-尘土滚动摩擦因数、尘土-尘土静摩擦因数、尘土接触模型JKR表面能系数、尘土-尘土碰撞恢复系数及尘土-钢碰撞恢复系数。

将筛分后的尘土和秸秆碎末称重，将秸秆碎末进行人工挑选，剩余的尘土进行称重，得出尘土所占比率。重复3次试验，求其平均值，得出尘土占比率为9.34%。尘土具有一定的粘性，在数值模拟过程中选用Hertz-Mindlin with JKR接触模型，设置生成秸秆外皮和秸秆內瓤颗粒的质量均为15.69 g，生成尘土质量为3.27 g。时间步长5.026×10s，总时间5 s。在漏斗上方建立颗粒工厂，每1 s生成颗粒质量为100 g，5 s后所有颗粒从漏斗下端流出，在下方圆盘上形成稳定的颗粒堆，通过后处理测定休止角。

2.2 数值模拟参数的标定

2

.

2

.

1

显著性因素的筛选

结合物理试验所获得的接触特性参数以及查阅文献，确定数值模拟中各参数的取值范围见表2，利用Design-Expert软件中的Plackett-Burman进行试验设计，以数值模拟颗粒堆休止角为响应值，筛选出对数值模拟休止角影响最显著的因素。

表2 休止角数值模拟影响因素Plackett-Burman试验参数取值
Table 2 Influencing factors of numerical simulation of angle of repose and values of Plackett Burman test parameters

试验参数Test parameter低水平Low-level高水平High-level尘土-钢滚动摩擦因数X1Soil-steel rolling friction factor X10.050.20尘土-钢静摩擦因数X2Soil-steel static friction factor X20.51.2尘土-外皮滚动摩擦因数X3Soil-straw skin rolling friction factor X30.20.5尘土-外皮静摩擦因数X4Soil-straw skin static friction factor X40.40.7尘土-內瓤滚动摩擦因数X5Soil-straw pulp rolling friction factor X50.30.6尘土-內瓤静摩擦因数X6Soil-straw pulp static friction factor X60.60.9尘土-尘土滚动摩擦因数X7Soil-soil rolling friction factor X70.10.2尘土-尘土静摩擦因数X8Soil-soil static friction factor X80.321.16

表2(续)

试验参数Test parameter低水平Low-level高水平High-level尘土接触模型JKR表面能系数X9JKR surface energy coefficient of soil contact model X9011.25尘土-尘土碰撞恢复系数X10Soil-soil restitution coefficient X100.150.75尘土-钢碰撞恢复系数X11Soil-steel restitution coefficient X110.20.5

在Plackett-Burman试验设计中，将表2内11个试验参数取值的最高值和最低值分别标为+1、-1 水平，设置一个中心点即选取每个试验参数最高值和最低值的平均值作为0水平，共进行13次数值模拟。数值模拟得出的尘土与秸秆碎末数值模拟模型堆见图4，利用EDEM软件中后处理的量尺工具测量颗粒堆的高度，计算得出数值模拟休止角

α

。

H为数值模拟颗粒堆高度；R为圆盘半径。H is numerical simulation of particle height; R is the radius of disk.图4 尘土与秸秆碎末数值模拟模型堆Fig.4 Numerical simulation model pile of soil and straw debris

漏斗下方的圆盘半径

R

为25 mm，已知颗粒堆的高度，利用式(4)求得颗粒堆的休止角(

α

)：

(4)

式中：

H

为数值模拟颗粒堆高度，mm；

R

为圆盘半径，mm。计算所得的数值模拟休止角(

α

)见表3，利用Design-Expert软件对结果进行方差分析。影响休止角的各个因素显著性排序见表4：影响休止角的最显著性因素有3个，其中尘土接触模型JKR表面能系数(

X

)对休止角的影响最大，其次为尘土-尘土碰撞恢复系数(

X

)和尘土-尘土滚动摩擦因数(

X

)。在进行后续数值模拟时，其余不具显著性影响的因素，参数取值均为平均值，分别为尘土-钢滚动摩擦因数0.125、尘土-钢静摩擦因数0.85、尘土-外皮滚动摩擦因数0.35、尘土-外皮静摩擦因数0.55、尘土-內瓤滚动摩擦因数0.45、尘土-內瓤静摩擦因数0.75、尘土-尘土静摩擦因数0.74、尘土-钢碰撞恢复系数0.35。

表3 休止角数值模拟影响因素Plackett-Burman试验设计及结果
Table 3 Influencing factors of numerical simulation of angle of repose and Plackett-Burman experimental design and results

序号Serialnumber各因素水平编码 Coding of each factor levelx1x2x3x4x5x6x7x8x9x10x11休止角α2/(°)Repose angle11-111-1111-1-1-141.3482-1-11-111-1111-142.6143-1111-1-1-11-11130.9644-11-111-1111-1-147.2035-1-1-11-111-111141.348

表3(续)

序号Serialnumber各因素水平编码 Coding of each factor levelx1x2x3x4x5x6x7x8x9x10x11休止角α2/(°)Repose angle6-1-1-1-1-1-1-1-1-1-1-138.6607-111-1111-1-1-1138.66080000000000040.03091-1-1-11-111-11141.3481011-1-1-11-11-11145.00011111-1-1-11-111-144.421121-1111-1-1-11-1143.8311311-1111-1-1-11-137.954

注：,,…,分别为,,…,的水平编码。,,…,各因素的物理意义见表2，下同。
Note： , , …, are the horizontal codes of , , …, respectively. The physical meanings of the factors , , …, are shown in Table 2. The same below.

表4 休止角数值模拟影响因素显著性排序
Table 4 Significance ranking of influencing factors of numerical simulation of angle of repose

因素Factor效应Effect均方和Mean square sum影响率/%Influence ratio排序SortX12.4117.438.764X2-0.832.051.039X3-1.617.783.917X4-1.345.392.718X51.658.134.086X60.080.020.5510X72.5519.499.793X80.601.080.5411X95.92104.9852.731X10-2.6821.5010.802X11-1.8410.155.105

2

.

2

.

2

最陡爬坡试验通过Plackett-Burman试验确定了3个对休止角影响最显著的因素，以数值模拟所测休止角(

α

)和物理试验休止角(

α

)的相对误差作为评价指标，利用最陡爬坡试验可迅速确定各因素的最佳参数取值范围。最陡爬坡试验设计及结果见表5：随着尘土接触模型JKR表面能系数(

X

)、尘土-尘土碰撞恢复系数(

X

)和尘土-尘土滚动摩擦因数(

X

)的不断增大，试验所测休止角先增大后减小，整体呈现一个抛物线形状。观察各试验相对误差发现，5号试验相对误差最小，所以最佳参数取值范围在试验5附近。继续取试验4、试验5和试验6参数分别为低、中和高水平进行后续数值模拟。

表5 休止角显著因素最陡爬坡试验设计及结果
Table 5 Design and results of the steepest climbing test for the apparent factors of angle of repose

序号Serial number因素 FactorX7X9X10休止角α2/(°)Repose angle相对误差/%Relative error10.1000.1537.9548.9120.12 2.250.2740.0303.9330.14 4.500.3943.8315.2040.16 6.750.5145.0008.0050.18 9.000.6341.3480.7860.2011.250.7540.0303.93

2

.

2

.

3

Box

-

Behnken试验设计及结果分析对最陡爬坡试验结果分析，以试验相对误差为评价指标，确定了试验最优参数取值范围。以试验4、试验5和试验6参数分别为低(-1)、中(0)和高水平(1)进行Box-Behnken试验设计(表6)。试验中取3个中心点对试验误差进行评估，得出Box-Behnken试验设计表(表7)。Box-Behnken试验共进行15次数值模拟，当数值模拟完成后，同样利用EDEM中后处理的量尺工具测量颗粒堆高度，求得数值模拟休止角(

α

)结果见表7。

表6 休止角显著因素最优参数取值
Table 6 Optimal parameter value of apparent factor of angle of repose

因素Factors低水平Low-level中水平Medium level高水平High-levelX70.160.180.20X96.759.0011.25 X100.510.630.75

表7 休止角显著因素Box-Behnken试验设计及结果
Table 7 Box-Behnken test design and results of significant factors of repose angle

序号Serial number因素水平编码 Coding of each factor levelx7x9x10休止角α2/(°)Repose angle100043.232-1-1038.88310-139.35410142.61501-142.33601142.617-10138.668-11041.32900042.96100-1-139.0011-10-140.031200043.56131-1039.351411045.32150-1140.03

利用Design-Expert软件对其结果进行多元回归分析，得出显著性参数与数值模拟休止角(

α

)的二次多项式方程为：

α

=43

.

28+0

.

97

X

+1

.

79

X

+0

.

40

X

+
0

.

88

X

X

+1

.

16

X

X

-0

.

19

X

X

-
1

.

45

X

-0

.

62

X

-1

.

67

X

(5)

对Box-Behnken试验进行方差分析，统计结果见表8。尘土-尘土滚动摩擦因数(

X

)、尘土接触模型JKR表面能系数(

X

)、尘土-尘土滚动摩擦因数与尘土接触模型JKR表面能系数(

X

X

)、尘土-尘土滚动摩擦因数与尘土-尘土碰撞恢复系数(

X

X

)、尘土-尘土滚动摩擦因数的二次项(

X

)和尘土-尘土碰撞恢复系数的二次项(

X

)对休止角的影响显著。模型的

P

<0.05，说明休止角与回归方程关系显著；失拟项的

P

=0.33>0.05，表明方程拟合较好；变异系数CV=1.18%，表明试验可靠性较高；决定系数

R

=0.980 7，校正决定系数两者的数值均接近1，模型可以较准确的反映实际状况。

表8 休止角显著因素Box-Behnken试验设计回归模型方差分析
Table 8 Box-Behnken experimental design regression model analysis of variance for significant factors of angle of repose

方差来源Source of variation均方和Mean square sum自由度Degree of freedom均方Mean squareP值P value模型 Model60.4096.710.000 9 X77.4917.490.002 5 X925.63125.630.000 1 X101.2811.280.067 8 X7X93.1213.120.015 2 X7X105.3615.360.005 1 X9X100.1410.140.476 2 X727.7417.740.002 3 X921.4111.410.058 8 X10210.33110.330.001 2残差 Residual1.1950.24失拟项 Lack of fit1.0030.330.222 5纯误差 Pure error0.1820.09总和 Sum61.5814

2

.

2

.

4

回归模型交互效应分析根据二次多项式模型方差分析结果，尘土-尘土滚动摩擦因数与尘土接触模型JKR表面能系数(

X

X

)、尘土-尘土滚动摩擦因数与尘土-尘土碰撞恢复系数(

X

X

)的

P

<0.05，这2个交互项对颗粒堆休止角影响显著。利用Design-Expert软件对尘土-尘土滚动摩擦因数与尘土接触模型JKR表面能系数(

X

X

)以及尘土-尘土滚动摩擦因数与尘土-尘土碰撞恢复系数(

X

X

)交互作用的响应曲面进行绘制，交互项对休止角的影响见图5。通过对比响应面坡度趋势得出，尘土接触模型JKR表面能系数(

X

)所对应的响应面坡度趋势最陡，尘土-尘土滚动摩擦因数(

X

)和尘土-尘土碰撞恢复系数(

X

)所对应的响应面坡度较为平缓。这说明尘土接触模型JKR表面能系数(

X

)对颗粒堆休止角的影响显著。

图5 部分因素交互作用对休止角(α2)的影响Fig.5 Effect of interaction of some factors on angle of repose

3 参数优化及试验验证

3.1 休止角验证试验

以秸秆碎末和尘土颗粒堆物理试验休止角

α

=41.668°为目标，利用Design-Expert软件对显著性参数与数值模拟休止角(

α

)的二次多项式方程(式(5))进行寻优求解，在软件提供的众多参数组合中寻找若干组组合，进行数值模拟。最后选择与物理试验休止角相对误差最小的一组组合为数值模拟最优取值，即尘土-尘土滚动摩擦因数0.182、尘土接触模型JKR表面能系数9.052和尘土-尘土碰撞恢复系数0.52。为验证上述参数的准确性，以上述参数再次进行3次重复数值模拟，所得数值模拟休止角

α

分别为40.356°、40.741°和42.332°，平均值为41.143°，与物理试验颗粒堆休止角

α

的相对误差为1.26%。

3.2 标准筛除尘率验证试验

为了进一步验证数值模拟参数取值的准确性，进行标准筛除尘率对比验证试验。标准筛直径为450 mm，高度150 mm，筛孔直径分别为3 mm和7 mm。利用标准筛进行除尘物理试验，将筛出物称重得出除尘率，重复6次试验并取平均值，得出物理试验除尘率为13.24%。利用三维建模软件Solidworks绘制标准筛三维模型，模型尺寸与实物尺寸相同，具体结构见图6。

图6 标准筛三维模型Fig.6 3D model of standard sieve

将三维模型导入EDEM软件中，在标准筛上添加正弦平移方向的力，模拟手动摇晃标准筛的运动状态，频率为10 Hz，振幅8 mm。在标准筛上方设置颗粒工厂，颗粒分别为段状的秸秆外皮、段状的秸秆內瓤、尘土、碎末状的秸秆外皮和碎末状的秸秆內瓤。生成段状的秸秆外皮和内瓤模型，模型直径为3 mm，长度20 mm，质量88 g，尘土颗粒和秸秆碎末颗粒模型不变，质量分别为2 和21 g。设置时间步长为5.026×10s，总时间8 s。数值模拟结束后，在后处理面板的Setup selections选项中选择Grid bin group，网格位于标准筛下方，设置网格数为1，尺寸为长220 mm、宽170 mm、高43 mm，在数据显示窗口中选择显示统计颗粒质量为尘土颗粒和秸秆碎末颗粒，最后选择显示颗粒属性即显示从标准筛中筛出的尘土颗粒和秸秆碎末颗粒质量。数值模拟所筛出的颗粒质量为28.44 g，除尘率为14.29%。标准筛数值模拟除尘率与物理试验除尘率相对误差为7.93%，再次验证了数值模拟参数取值的准确性。

4 结 论

本研究采用物理试验与数值模拟相结合的方法对玉米秸秆饲料的筛出物(尘土与秸秆碎末)进行离散元数值模拟参数标定，分别采用Plackett-Burman试验、最陡爬坡试验和Box-Behnken试验确定对数值模拟休止角影响最显著的因素以及参数取值，主要结论如下：

1)利用斜面仪测定尘土-外皮滚动摩擦因数为0.2～0.5，尘土-外皮静摩擦因数0.4～0.7，尘土-內瓤滚动摩擦因数0.3～0.6，尘土-內瓤静摩擦因数0.6～0.9。

2)通过Design-Expert软件进行Plackett-Burman试验和最陡爬坡试验得到了3个对颗粒堆休止角影响最显著因素，分别为尘土接触模型JKR表面能、尘土-尘土滚动摩擦因数和尘土-尘土恢复系数。

3)以物理试验颗粒堆休止角为目标值，通过Box-Behnken试验以及方差分析得到了显著性参数最优解，分别为尘土接触模型JKR表面能系数9.052 J/m、尘土-尘土滚动摩擦因数0.182、尘土-尘土恢复系数0.52。对优化后的参数进行数值模拟得到休止角为41.143°。

4)对所得出的试验数据分别进行休止角验证试验和标准筛除尘率验证试验，数值模拟休止角与物理试验休止角相对误差为1.26%，标准筛数值模拟除尘率与物理试验除尘率相对误差为7.93%，进一步验证了数值模拟参数的可靠性。