基于离散元的鱼饲料仿真参数标定与试验

郭三琴，万鹏，杨俊，邓志勇，汪帅帅，朱明，谭鹤群; GUO Sanqin; WAN Peng; YANG Jun; DENG Zhiyong; WANG Shuaishuai; ZHU Ming; TAN Hequn

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

基于离散元的鱼饲料仿真参数标定与试验 PDF

- ORCID：
郭三琴 ¹
- ORCID：
万鹏 ^1,2
✉
- ORCID：
杨俊 ¹
- ORCID：
邓志勇 ¹
- ORCID：
汪帅帅 ¹
- ORCID：
朱明 ^1,2
- ORCID：
谭鹤群 ^1,2

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

中图分类号： S200.1； TH6

最近更新：2022-12-08

DOI：10.13300/j.cnki.hnlkxb.2022.06.034

摘要

为全面研究鱼饲料的离散元仿真物性参数，以渔用硬颗粒3.0中号饲料为研究对象，通过密度测量试验与单轴平板压缩试验测得硬颗粒饲料的密度、弹性模量、剪切模量、泊松比数值分别为1 096 kg/m³、336.32 MPa、126.65 MPa、0.33。同时在EDEM仿真软件中建立硬颗粒饲料离散元模型，结合台架试验与仿真试验对硬颗粒饲料-钢以及硬颗粒饲料-尼龙之间的碰撞恢复系数、静摩擦因数、滚动摩擦因数进行标定，其值分别为0.380、0.538、0.028，0.405、0.356、0.033。结合堆积角试验，通过最陡爬坡试验和Box-Behnken响应面优化试验，标定硬颗粒饲料间的碰撞恢复系数为0.364、静摩擦因数为0.236、滚动摩擦因数为0.13。运用自主研发的投饵机进行参数验证试验，以硬颗粒饲料分布变异系数C_V为评价指标，得出实际硬颗粒饲料分布变异系数与仿真硬颗粒饲料分布变异系数的相对误差为1.709%，表明本研究中建立的渔用硬颗粒饲料的离散元模型以及标定的硬颗粒饲料物性参数可用于EDEM离散元仿真试验。

关键词

鱼饲料; 颗粒饲料; 离散元; 本征参数; 接触参数; 标定; 仿真参数

投饵机的结构是影响投饲均匀性的重要因素。在对投饵机结构进行优化时通常需要研究饲料颗粒与饲料颗粒间、饲料颗粒与投饵机接触部件间的相互作用。由于饲料颗粒与饲料颗粒间、饲料颗粒与投饵机接触部件间的相互作用情况复杂，基于传统理论分析和经验公式设计的投饵机与实际生产出来的投饵机其投饲效果之间通常存在较大差异^［

1-2］。随着数值模拟软件在农业工程领域的不断应用，利用离散元法（discrete element method，DEM）对农业装备关键部件的结构进行辅助设计可以显著提高设计效果^{［参考文献 3-4}3-4］。基于离散元法对投饵机投饲过程进行仿真模拟时，为确保仿真结果准确度，需选取正确的离散元模型，设定精确的物性参数^{［参考文献 5-6}5-6］。物性参数包括本征参数和接触参数，本征参数是指物体本身的特征参数，与实际数值一致，可通过台架试验测定。由于颗粒模型与真实颗粒形状不完全一致，导致仿真接触参数与真实值存在较大差异，且有些参数直接测定难度较大^{［参考文献 7

百度学术}7］，因此需要结合台架试验与仿真试验对接触参数进行标定。

目前，国内外学者基于离散元法已对多种农业散粒体物料进行了参数标定。Coetzee等^［

8］结合剪切和压缩试验标定了玉米离散元模型的摩擦因数和刚度。焦俊等^{［参考文献 9

百度学术}9］利用单轴压缩试验标定了莲藕主藕体的本征参数。张锐等^{［参考文献 10

百度学术}10］利用试验堆积角和模拟堆积角对比的方法标定了沙土间接触参数。彭飞等^{［参考文献 11

百度学术}11］通过注入截面法直接获取大猪饲料颗粒的堆积轮廓图，并对其进行休止角测定，从而标定了大猪饲料间的接触参数。此外，部分学者利用台架试验与仿真试验相结合的休止角试验对水稻^{［参考文献 12

百度学术}12］、小麦^{［参考文献 13

百度学术}13］、微型马铃薯^{［参考文献 14

百度学术}14］、麻山药^{［参考文献 15

百度学术}15］、三七种子^{［参考文献 16

百度学术}16］等农产品进行了物性参数标定。综上所述，国内外研究人应用离散元仿真参数标定的对象主要为农产品、土壤、粉末状饲料等，但对渔用颗粒饲料的物性参数未有较为全面的研究，相关参考参数缺乏。因此，对渔用硬颗粒饲料的物性参数进行全面标定具有重要意义。

本研究以淡水养殖渔用硬颗粒3.0中号饲料为研究对象，通过台架试验测定硬颗粒饲料的本征参数，在EDEM中建立硬颗粒饲料的离散元模型，结合台架试验与仿真试验，标定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的接触参数，结合堆积角试验，通过最陡爬坡试验、Box-Behnken响应面优化试验，确定硬颗粒饲料间的最优接触参数组合，通过自主设计的投饵机进行参数验证试验，以期得到全面的硬颗粒饲料物性参数，进一步为投饵机结构设计与优化提供参考。

1 材料与方法

1.1 物料本征参数测定

物料的本征参数主要包括密度、剪切模量、弹性模量、泊松比等^［

17］。本研究以通威126型渔用3.0中号柱形硬颗粒饲料为试验对象，开展水产饲料的本征参数测定研究。通过测定，试验所用硬颗粒饲料样本的含水率为9%~10%，直径为3.1~3.3 mm，长度为6.6~8.6 mm。

1）密度测定。随机选取完整硬颗粒3.0中号饲料100颗，利用电子分析天平（日本岛津AUY220型，精度0.1 mg）测量每颗饲料质量m_i，利用数显游标卡尺（精度0.01 mm）测量对应饲料颗粒的直径d_i和长度l_i，计算每颗饲料的体积v_i，进而计算得出每颗饲料的密度ρ_i。计算所选100颗硬颗粒饲料样本的密度及平均值，可得硬颗粒饲料的平均密度为ρ =（1 096±27） kg/m³。

2）弹性模量、剪切模量与泊松比的测定。随机选取完整硬颗粒3.0中号饲料20颗，利用质构仪（TMS-Pro型，美国FTC公司）开展单轴平板压缩试验（图1）。设置平板压缩速度为1 mm/min，起始力为0.2 N，在确保饲料不发生破碎断裂的情况下最大压缩距离设定为1 mm，测量每颗饲料压缩前后的直径与长度，计算得到饲料压缩前后的直径及长度变化量。

图1 单轴平板压缩试验

Fig.1 Uniaxial plate compression test

对选取的20颗硬颗粒饲料样本进行试验，参考文献［

18］中物理公式计算每颗硬颗粒饲料的弹性模量、剪切模量、泊松比的值并计算各参数的平均值，可得硬颗粒饲料的弹性模量、剪切模量、泊松比分别为（336.32±9.80） MPa、（126.65±4.13） MPa、0.33±0.02。

1.2 物料与接触面间接触参数测定

1）饲料颗粒仿真模型。渔用硬颗粒饲料外观呈较为规则的柱形，随机选取完整的硬颗粒3.0中号饲料100颗，量取每颗饲料直径d_i与长度l_i，计算平均值，可得硬颗粒饲料平均直径d为3.25 mm，平均长度l为7 mm。根据饲料的实际形状及参数，采用球颗粒组合的方式^［

19］，建立硬颗粒饲料模型，如图2所示。

图2 硬颗粒饲料离散元几何模型

Fig.2 DEM model of hard pellet feed

硬颗粒饲料球颗粒直径取3.25 mm，球颗粒接触模型选择Hertz-Mindlin无滑移模型。结合上文对硬颗粒饲料本征参数的测定结果，并参照文献［

20-21］中接触材料钢材、尼龙的本征参数，确定硬颗粒饲料与接触材料的本征参数，如表1所示。

表1 饲料与材料的本征参数

Table 1 Intrinsic parameter of feed and materials

材料Materials

泊松比

Poisson’s ratio

剪切模量/MPa

Shear modulus

密度/(kg/m³)

Density

硬颗粒饲料

Hard pellet feed

0.33

126.65

1 096

钢材Steel

0.30

7×10⁴

7 800

尼龙Nylon

0.40

1×10⁸

1 500

2）碰撞恢复系数。利用碰撞弹跳试验测定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的碰撞恢复系数。试验操作过程如图3A所示。以硬颗粒饲料-钢为例，将钢板水平放置于试验台上，将渔用硬颗粒3.0中号饲料从高度H₁=200 mm处释放，饲料自由落体接触到钢板时反向弹起，使用高速摄像机采集记录饲料反弹时的最高位置h₁；同时，为保证台架试验与仿真试验的一致性，颗粒下落时，保持柱形颗粒的下底面接触钢板，进行10次台架重复试验，计算平均值，可得硬颗粒饲料-钢、硬颗粒饲料-尼龙实际最高反弹高度分别为h_a=（31.70±3.65） mm、h_b=（35.60±2.88） mm。

图3 碰撞恢复系数标定试验

Fig.3 Calibration experiment of restitution coefficient

A：台架试验 Bench test； B：仿真试验 Simulation test； 1.硬颗粒饲料Hard pellet feed； 2.坐标纸 Coordinate paper； 3.钢板Steel plate.

EDEM仿真试验如图3B所示，选取已建立的硬颗粒饲料模型，输入已标定硬颗粒饲料和钢的本征参数，设置饲料颗粒进行初速度为0的自由落体运动，由于硬颗粒饲料与钢之间的静摩擦因数、滚动摩擦因数，硬颗粒饲料间的碰撞恢复系数、静摩擦因数、滚动摩擦因数对反弹高度没有影响，为避免干扰将这5个参数均设置为0。

为了确定硬颗粒饲料在不同材料表面的碰撞恢复系数，以实际最高反弹高度值h_a、h_b为目标值，在EDEM中设置不同的参数进行大量预仿真试验，确定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的碰撞恢复系数x₁_a、x₁_b均在0.35~0.45范围内，取步长为0.05，每个组合进行6组仿真试验，y₁_a表示对应硬颗粒饲料-钢之间的碰撞恢复系数x₁_a的仿真最高反弹高度值，y₁_b表示对应硬颗粒饲料-尼龙之间的碰撞恢复系数x₁_b的仿真最高反弹高度值，试验设计方案与结果如表2所示。

表2 碰撞恢复系数仿真试验方案

Table 2 Simulation test scheme and results of restitution coefficient

序号 No.	x₁_a(x₁_b)	y₁_a/mm	y₁_b/mm
1	0.35	27.43	27.56
2	0.37	30.17	30.22
3	0.39	33.22	33.39
4	0.41	36.37	36.46
5	0.43	39.63	39.68
6	0.45	43.20	43.19

对表2中的数据进行曲线拟合，构建硬颗粒饲料和接触材料间的碰撞恢复系数与仿真最高反弹高度值之间的拟合方程如下：

\{\begin{matrix} y_{1 a} = 224.88 x_{1 a}^{2} - 22.23 x_{1 a} + 7.65 \\ y_{1 b} = 200.58 x_{1 b}^{2} - 3.89 x_{1 b} + 4.31 \end{matrix}

（1）

式（1）中，2个二次拟合方程的决定系数R²值依次为0.999 9、0.999 8，均接近1，表明这2个拟合方程均具有较高的可靠性。将通过台架试验测得实际最高反弹高度h_a（31.70 mm）、h_b（35.60 mm）分别带入式（1）中2个二次拟合方程，得出饲料与接触材料之间的碰撞恢复系数x₁_a、x₁_b分别为0.38、0.405；将所得结果带入EDEM软件进行验证，每组试验重复5次，计算平均值，可得仿真最高反弹高度值分别为 $h_{a}^{'}$ =31.74 mm、 $h_{b}^{'}$ =35.75 mm，与目标值的相对误差分别为0.13%、0.42%，表明标定后的仿真结果与实测值基本一致，从而确定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的碰撞恢复系数分别为0.38、0.405。

3）静摩擦因数。利用斜面滑移试验测定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的静摩擦因数。试验操作过程如图4所示。

图4 静摩擦因数标定试验

Fig.4 Calibration experiment of static friction coefficient

A：台架试验 Bench test； B：仿真试验 Simulation test；1.钢板 Steel plate；2.硬颗粒饲料 Hard pellet feed.

以硬颗粒饲料-钢为例，将渔用硬颗粒3.0中号饲料放置于水平放置的钢板上，再将钢板一端绕另一端的底边缓慢匀速转动，直到硬颗粒饲料开始发生滑动，停止转动并记录钢板最大倾斜角度；为防止试验过程中硬颗粒饲料滚动，放置饲料时，保证柱形饲料中心轴与钢板旋转轴垂直。每个组合进行10次台架重复试验，计算平均值，可得硬颗粒饲料-钢、硬颗粒饲料-尼龙实际最大倾斜角分别为α_a=（27.91±1.84）°、α_b=（19.44±1.2）°。

EDEM仿真试验时，选取已建立的硬颗粒饲料模型，输入已标定硬颗粒饲料与钢的本征参数以及硬颗粒饲料与钢之间的碰撞恢复系数，由于硬颗粒饲料与钢之间的滚动摩擦因数，硬颗粒饲料间的碰撞恢复系数、静摩擦因数、滚动摩擦因数对钢板最大倾斜角没有影响，为避免干扰，将这4个参数均设置为0。

为了确定硬颗粒饲料在不同材料表面的静摩擦因数，设置饲料初速度为0，以实际最大倾斜角α_a、α_b作为目标值，在EDEM中设置不同的参数进行大量预仿真试验，确定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的静摩擦因数x₂_a、x₂_b范围分别为0.40~0.65、0.25~0.50，取步长为0.05，每个组合进行6组仿真试验；y₂_a表示对应硬颗粒饲料-钢之间的静摩擦因数x₂_a的仿真最大倾斜角，y₂_b表示对应硬颗粒饲料-尼龙之间的静摩擦因数x₂_b的仿真最大倾斜角，试验设计方案与结果如表3所示。

表3 静摩擦因数仿真试验方案

Table 3 Simulation test scheme and results of static friction coefficient

序号 No.	硬颗粒饲料-钢 Hard pellet feed and steel		硬颗粒饲料-尼龙 Hard pellet feed and nylon
序号 No.	x₂_a	y₂_a/(°)	x₂_b	y₂_b/(°)
1	0.40	21.2	0.25	14.0
2	0.45	23.8	0.30	16.6
3	0.50	26.3	0.35	19.1
4	0.55	28.3	0.40	21.6
5	0.60	30.5	0.45	23.8
6	0.65	32.6	0.50	25.9

对表3中的数据进行曲线拟合，构建硬颗粒饲料与接触材料间的静摩擦因数与仿真最大倾斜角之间的关系拟合方程如下：

\{\begin{matrix} y_{2 a} = - 26.43 x_{2 a}^{2} + 72.95 x_{2 a} - 3.71 \\ y_{2 b} = - 26.43 x_{2 b}^{2} + 67.59 x_{2 b} - 1.27 \end{matrix}

（2）

式（2）中，2个二次拟合方程的决定系数R²值依次为0.999 5、0.999 9，均接近1，表明这2个拟合方程均可靠有效。将通过台架试验测得实际最大倾斜角α_a（27.91°）、α_b（19.44°）依次带入式（2）中对应的2个拟合方程，得出饲料与接触材料之间的静摩擦因数分别为x₂_a=0.538、x₂_b=0.356，将2个静摩擦因数分别带入EDEM软件进行验证，每组试验重复5次，计算平均值，得出钢板与尼龙板最大倾斜角分别为 $α_{a}^{'}$ =28.2°、 $α_{b}^{'}$ =19.0°，与目标值的相对误差依次为1.04%、2.26%，表明标定后的仿真结果与实测值基本一致，从而确定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的静摩擦因数分别为0.538、0.356。

4）滚动摩擦因数。利用斜面滚动试验测定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的滚动摩擦因数，试验操作过程如图5所示。

图5 滚动摩擦因数标定试验

Fig.5 Calibration experiment of rolling friction coefficient

A：台架试验 Bench test； B：仿真试验 Simulation test； 1.倾斜钢板 Inclined steel plate； 2.平置钢板 Flat steel plate； 3.硬颗粒饲料 Hard pellet feed.

以硬颗粒饲料-钢为例，将折弯过的钢板水平放置在试验台上，倾斜钢板倾斜角β=20°，将渔用硬颗粒3.0中号饲料从高度H₂=20 mm处释放，饲料沿斜面做初速度为0的滚动运动，在摩擦阻力的作用下，饲料滚动一段距离后停止。使用刻度尺测量饲料水平滚动的距离s，为保证实际试验环境与仿真试验的一致性，放置饲料时，保证柱形饲料中心轴与钢板折弯轴平行。进行10次台架重复试验计算平均值，可得硬颗粒饲料-钢、硬颗粒饲料-尼龙的实际水平滚动距离分别为s_a=（278.52±5.25） mm、s_b=（234.97±6.08） mm。

EDEM仿真试验时，选取已建立的硬颗粒饲料模型，输入已标定硬颗粒饲料与钢的本征参数、硬颗粒饲料与钢之间的碰撞恢复系数和静摩擦因数，由于硬颗粒饲料间的碰撞恢复系数、静摩擦因数、滚动摩擦因数对饲料在钢板上的水平滚动距离无影响，为避免干扰，将这3个参数均设置为0。为了确定硬颗粒饲料在不同材料表面的滚动摩擦因数，设置饲料初速度为0，将实际水平滚动距离s_a、s_b作为目标值，在EDEM中设置不同的参数进行大量预仿真试验，确定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的滚动摩擦因数x₃_a、x₃_b范围分别为0.026~0.036、0.03~0.055，取步长为0.05，每个组合进行6组仿真试验；y₃_a表示对应硬颗粒饲料-钢之间的滚动摩擦因数x₃_a的仿真水平滚动距离、y₃_b表示对应硬颗粒饲料-尼龙之间的滚动摩擦因数x₃_b的仿真水平滚动距离，试验设计方案与结果如表4所示。

表4 滚动摩擦因数仿真试验方案

Table 4 Simulation test scheme and results of rolling friction coefficient

序号 No.	硬颗粒饲料-钢 Hard pellet feed and steel		硬颗粒饲料-尼龙Hard pellet feed and nylon
序号 No.	x₃_a	y₃_a/mm	x₃_b	y₃_b/mm
1	0.026	297.74	0.030	262.14
2	0.028	281.66	0.035	218.48
3	0.030	259.50	0.040	188.47
4	0.032	230.98	0.045	166.06
5	0.034	224.31	0.050	140.20
6	0.036	210.07	0.055	124.28

对表4中的数据进行曲线拟合，构建硬颗粒饲料与接触材料间的滚动摩擦因数与仿真水平滚动距离之间的关系拟合方程如下：

\{\begin{array}{l} y_{3 a} = 319933.04 x_{3 a}^{2} - 28956.91 x_{3 a} + 837.2 \\ y_{3 b} = 110960 x_{3 b}^{2} - 14840.61 x_{3 b} + 605.48 \end{array}

（3）

式（3）中，2个二次拟合方程的决定系数R²值依次为0.984 4、0.997 6，均接近1，表明这2个拟合方程均可靠有效。将通过台架试验测得实际水平滚动距离s_a（278.52 mm）、s_b（234.97 mm）依次带入式（3）中2个二次拟合方程，得出饲料与接触材料之间的滚动摩擦因数分别为x₃_a=0.028、x₃_b=0.033，将2个滚动摩擦因数分别带入EDEM软件进行验证，每组试验重复5次，计算平均值，可得仿真水平滚动距离分别为 $s_{a}^{'}$ =281.66 mm、 $s_{b}^{'}$ =237.13 mm，与目标值的相对误差分别为1.13%、0.92%，表明标定后的仿真结果与实测值基本一致，从而确定硬颗粒饲料-钢、硬颗粒饲料-尼龙之间的滚动摩擦因数分别为0.028、0.033。

2 结果与分析

2.1 饲料颗粒间接触参数标定

本研究采用圆筒提升法进行堆积角试验对硬颗粒饲料间的接触参数进行标定。由于前文已完成硬颗粒饲料与接触材料间的接触参数标定，进行硬颗粒饲料间接触参数标定时接触材料选择1种即可。本研究选用内径60 mm、长度70 mm的钢管作为接触材料进行堆积角试验，对硬颗粒饲料间的接触参数进行标定。

1）饲料堆积角试验。以渔用硬颗粒3.0中号饲料为研究对象进行堆积角试验，具体试验过程如图6所示。将钢管竖直固定在万能材料试验机上，钢板水平放置。试验开始时，保持钢管下端面与钢板表面接触，在钢管内填充饲料颗粒与钢管上端面平齐；开启电源，设定万能材料试验机提升速度为300 mm/min，钢管随万能材料试验机匀速缓慢上升，饲料颗粒自由扩散在钢板上形成稳定堆积角。

图6 堆积角试验

Fig.6 Stacking angle test

A：台架试验 Bench test； B：仿真试验 Simulation test； 1.万能材料试验机Universal material testing machine； 2.钢管Steel pipe； 3.硬颗粒饲料 Hard pellet feed； 4.钢板Steel plate.

采集饲料颗粒堆积图像，利用灰度化、二值化等图像处理算法提取堆积边缘线并进行线性拟合，则拟合直线与水平面间的夹角即为实际堆积角，如图7所示。

图7 堆积角图像处理

Fig.7 Image processing of stacking angle

A：原始图像Original image；B：二值化处理 Binarization treatment； C：轮廓拟合直线 Contour fitting line.

称量可得钢板上硬颗粒饲料的质量为136.5 g，进行10次重复试验，统计每次试验的饲料堆积角并计算其平均值，可得硬颗粒饲料实际堆积角γ=（31.1±1.70） °。

采用EDEM软件进行仿真试验（图6B），选取建立的硬颗粒饲料模型，输入已标定的硬颗粒饲料与钢的本征参数以及硬颗粒饲料与钢之间的接触参数，设置饲料颗粒质量为136.5 g，生成时间1.5 s，待颗粒生成且稳定后，钢管以300 mm/min匀速上升，钢管中饲料颗粒缓慢扩散堆积，形成稳定堆积角。获取仿真试验得到的饲料颗粒堆积图像进行灰度化、二值化等处理，提取堆积边缘线进行直线拟合，拟合直线与水平面间的夹角即为仿真堆积角γ´。

为了确定硬颗粒饲料间的碰撞恢复系数、静摩擦因数、滚动摩擦因数的取值，以台架试验测得实际堆积角γ作为目标堆积角，以仿真堆积角与目标堆积角的相对误差δ作为评价指标，设置不同的参数进行大量预仿真试验，确定硬颗粒饲料间的碰撞恢复系数、静摩擦因数、滚动摩擦因数范围分别为0.20~0.50、0.10~0.40、0.05~0.20，在此范围内进行最陡爬坡试验。

2）最陡爬坡试验。以台架试验获得的硬颗粒饲料实际堆积角γ作为目标堆积角，硬颗粒饲料间碰撞恢复系数X₁、静摩擦因数X₂、滚动摩擦因数X₃为试验因素，仿真堆积角γ´与目标堆积角γ的相对误差δ作为评价指标进行最陡爬坡试验，通过寻找最为接近目标值的仿真堆积角，确定硬颗粒饲料间3个接触参数的最优组合。仿真堆积角与目标堆积角相对误差δ计算式为：

δ = \frac{|γ^{'} - γ|}{γ} \times 100 %

（4）

式（4）中，γ´为仿真堆积角，（°）；γ为目标堆积角，（°）。

最陡爬坡试验设计方案与结果如表5所示。由表5可看出，随着试验因素值增大，硬颗粒饲料仿真堆积角与目标堆积角相对误差δ呈现出先减小后增大趋势，且在4号试验得到最小相对误差值。因此，将4号组作为0水平，3号和5号分别作为低水平和高水平进行BBD（Box-Behnken design）响应面试验设计。

表5 最陡爬坡试验方案

Table 5 The scheme and results of steepest climbing test

序号 No.	X₁	X₂	X₃	γ΄/(°)	δ/%
1	0.20	0.10	0.050	21.4	31.19
2	0.25	0.15	0.075	24.5	21.22
3	0.30	0.20	0.100	29.0	6.75
4	0.35	0.25	0.125	32.0	2.81
5	0.40	0.30	0.150	36.5	14.79
6	0.45	0.35	0.175	40.0	22.25
7	0.50	0.40	0.200	45.0	30.89

3）方差分析与回归模型建立。基于Design-Expert 11进行BBD响应面试验设计，以硬颗粒饲料仿真堆积角γ´与目标堆积角γ的相对误差 δ 作为试验指标，X₁、X₂、X₃分别表示硬颗粒饲料间碰撞恢复系数、静摩擦因数、滚动摩擦因数，试验因素编码如表6所示。试验中心设置5组重复，试验结果如表7所示。

表6 试验因素编码

Table 6 Coding of factors

编码 Code	X₁	X₂	X₃
-1	0.30	0.20	0.100
0	0.35	0.25	0.125
1	0.40	0.30	0.150

表7 响应面试验设计与结果

Table 7 Design and results of Box-Behnken test

序号 No.	X₁	X₂	X₃	γ΄/(°)	δ/%
1	0.30	0.20	0.125	34.8	11.897
2	0.40	0.20	0.125	32.5	4.502
3	0.30	0.30	0.125	33.8	8.682
4	0.40	0.30	0.125	33.5	7.717
5	0.30	0.25	0.100	28.0	9.968
6	0.40	0.25	0.100	33.4	7.395
7	0.30	0.25	0.150	34.5	10.932
8	0.40	0.25	0.150	29.2	6.109
9	0.35	0.20	0.100	28.6	8.039
10	0.35	0.30	0.100	32.1	3.125
11	0.35	0.20	0.150	30.0	3.537
12	0.35	0.30	0.150	33.6	8.039
13	0.35	0.25	0.125	31.7	1.929
14	0.35	0.25	0.125	30.5	1.929
15	0.35	0.25	0.125	31.9	2.572
16	0.35	0.25	0.125	31.5	1.286
17	0.35	0.25	0.125	30.4	2.251

应用Design-Expert 11软件对试验结果进行多元回归拟合分析，可得硬颗粒饲料仿真堆积角与目标堆积角相对误差δ的回归方程，对回归方程进行显著性检验，结果如表8所示。由表8可知，二次回归模型（P<0.000 1）极显著，失拟项（P=0.824 9>0.05）不显著，回归方程不失拟，拟合效果好，具有预测意义。X₁、X₁X₂、X₂X₃、X₁²、X₂²、X₃²对δ 影响极显著，X₁X₃影响显著，X₂、X₃影响不显著，各因素对试验指标δ 的影响由大到小分别为：碰撞恢复系数、静摩擦因数、滚动摩擦因数。经逐步回归法剔除影响不显著项，得各因素与试验指标δ 回归方程为：

表8 回归方程方差分析

Table 8 Variance analysis of regression equation

方差来源 Source of variation	平方和 Quadratic sum	自由度 Freedom	F值 F value	P值 P value
模型 Model	191.637	9	133.73	<0.000 1^**
X₁	31.031	1	194.89	<0.000 1^**
X₂	0.021	1	0.13	0.725 9
X₃	0.001	1	0.01	0.938 7
X₁X₂	10.336	1	64.92	<0.000 1^**
X₁ X₃	1.266	1	7.95	0.025 8^*
X₂ X₃	22.165	1	139.21	<0.000 1^**
X₁²	87.592	1	550.12	<0.000 1^**
X₂²	11.394	1	71.56	<0.000 1^**
X₃²	17.635	1	110.76	<0.000 1^**
残差 Residual	1.115	7
失拟项 Lack of fit	0.205	3	0.30	0.824 9
纯误差 Pure error	0.910	4
总和 Sum	192.751	16

注Note： “*”表示显著，“**”表示极显著。 “*”means significant，“**”means very significant.

\begin{array}{l} δ = 1.99 - 1.97 X_{1} + 1.61 X_{1} X_{2} - 0.56 X_{1} X_{3} + \\ 2.35 X_{2} X_{3} + 4.56 X_{1}^{2} + 1.65 X_{2}^{2} + 2.05 X_{3}^{2} \end{array}

（5）

由表8可知，虽然X₂、X₃对试验指标δ 无显著影响，但2个因素的平方项X₂²、X₃²及因素间交互项X₁X₂、X₂X₃对δ 影响极显著，交互项X₁X₃对δ 影响显著，因此不可忽视这2个因素对试验指标的影响。

4）响应面分析和参数优化。采用Design-Expert 11软件进行响应面分析，以直观分析因素间交互作用对试验指标的影响情况。设定硬颗粒饲料间碰撞恢复系数X₁为0.35、静摩擦因数X₂为0.25、滚动摩擦因数X₃为0.125，因素间交互作用对硬颗粒饲料堆积角相对误差δ 的影响如图8所示。

图8 试验因素对硬颗粒饲料堆积角相对误差的影响Fig.8 Effects of experimental factors on relative error of hard pellet feed stacking angle

A：X₃=0.125；B：X₂=0.25；C：X₁=0.35.

由图8可知，随着碰撞恢复系数X₁、静摩擦因数X₂、滚动摩擦因数X₃的增加，硬颗粒饲料堆积角相对误差δ均呈现出先减小后增大的趋势，且3个因素对δ的影响显著性由大到小依次为X₁、X₂、X₃，与方差分析所得结果一致。

应用Design-Expert 11软件的参数优化模块，以表6中各因素范围为约束条件，以试验指标δ 回归方程式（5）为目标函数，求解其最小值，得出硬颗粒饲料间接触参数最优参数组合为碰撞恢复系数0.364，静摩擦因数0.236，滚动摩擦因数0.13，此时，堆积角相对误差为1.721%。将所得的最优参数组合进行仿真堆积角试验验证，得硬颗粒饲料仿真堆积角γ´=31.6°，与硬颗粒饲料实际堆积角γ（31.1°）之间的相对误差为1.608%，与优化结果基本吻合，表明所得饲料间最佳参数组合可用于EDEM仿真试验。

2.2 优化参数验证

为了验证渔用硬颗粒饲料离散元模型以及仿真参数的准确性，采用自主研发设计的投饵机样机进行硬颗粒饲料抛撒试验。该投饵机主要由料箱、下料装置、抛撒装置组成，工作原理为料箱中的饲料靠自身重力作用落入下料装置中下料辊的凹槽内，下料装置旋转下料，凹槽内饲料随之旋转至下料装置出料口时，饲料自由落入抛撒装置中匀速转动的离心转盘上，饲料在离心力的作用下沿离心转盘边缘的切线方向抛出，抛出的饲料经圆台挡板改变其运动方向，最终落到收集区域内，设置转盘转速为300 r/min。其中，下料辊、离心转盘为尼龙材料，其他部件均为钢材。硬颗粒饲料抛撒试验如图9所示。

图9 硬颗粒饲料抛撒试验

Fig.9 Centrifugal dispersion test of the hard pellet feed

A：台架试验Bench test； B：仿真试验Simulation test； 1.料箱Hopper； 2.下料辊Blanking roller；3.离心转盘Centrifugal turntable； 4.圆台挡板Frustum baffle； 5.收集区域 Collection area； 6. 饲料质量检测区域 Feed quality inspection area； 7. 落料圆板 Circular plate for collecting feed.

验证试验以硬颗粒饲料分布变异系数C_V为评价指标，饲料收集区域为直径4 m的圆面，其网格划分方式如图10所示。

图10 区域划分网格布置图

Fig.10 Area division grid layout

1.离心转盘Centrifugal turntable；2.计算网格Computational grid； 3.承载面 Bearing surface.

将圆面沿径向划分为面积相等的5个圆环面，沿周向划分为8个面积相等的圆锥面，即将整个饲料收集区域划分为40个面积相等的网格区域，每个网格区域的面积为0.314 m²。试验时，使收集区圆面圆心位于离心转盘中心线上，圆面距离离心转盘底部的距离为1.0 m。由于饲料颗粒下落接触地面时会产生弹跳、滑动、滚动等运动，影响试验结果的准确性，且与饲料颗粒落入水面时的实际状态不符。因此，本研究选用加绒地毯对下落的饲料颗粒进行收集。

抛撒试验完成后，统计硬颗粒饲料落在每个网格区域内的质量m_i，并计算饲料分布变异系数C_V。

\{\begin{matrix} C_{V} = \frac{S D}{\bar{m}} \times 100 % \\ S D = \sqrt[]{\frac{\sum_{i = 1}^{n} (m_{i} {- \bar{m})}^{2}}{n - 1}} \end{matrix}

（6）

式（6）中， C_V 为饲料分布变异系数，%；SD 为标准差，g；n 为网格区域数量，n=40； $\bar{m}$ 为n个网格区域内收集的饲料质量平均值，g；m_i 为第i个网格区域收集的饲料质量，g。

选用通威126型渔用硬颗粒3.0中号饲料1 000 g进行饲料抛撒试验，将饲料颗粒装入料箱，开启电源进行抛撒试验，待饲料颗粒全部抛撒完后关闭电源，统计落在每个网格区域内的硬颗粒饲料质量m_i，并计算分布变异系数C_V，试验重复5次计算平均值，得实际硬颗粒饲料分布变异系数C_V=1.092%。

根据台架试验参数进行仿真饲料抛撒试验，如图9B所示。将投饵机模型以及硬颗粒饲料模型导入EDEM离散元软件，同时输入已标定的硬颗粒饲料物性参数；在投饵机下方放置直径为4 m的落料圆板，设置硬颗粒饲料与落料圆板间的碰撞恢复系数为0、静摩擦因数为1、滚动摩擦因数为1，确保饲料落在落料圆板上时不发生运动。设置落料圆板圆心位于离心转盘中心线上，圆板上表面位于离心转盘底部的距离为1 m。在落料圆板上方建立饲料质量检测区域，绘制如图10所示的饲料质量检测区域模型，并将其导入EDEM软件后处理模块Setup Selections。设置饲料质量检测区域模型的位置参数，使饲料质量检测区域模型下表面与落料圆板上表面完全重合。最后设置饲料质量检测区域中饲料颗粒的检测方式为各个网格区域内的饲料累计质量。仿真试验结束后，统计饲料检测区域中经过各网格区域内的硬颗粒饲料模型的质量m_i'，并计算其分布变异系数C_V'，试验重复5次计算平均值，得仿真硬颗粒饲料分布变异系数C_V'=1.073%。

根据台架试验和仿真试验结果可知，硬颗粒饲料在进行台架试验和仿真试验时的分布变异系数相对误差为1.74%，且饲料颗粒在收集区域中的分布情况基本相同，这表明硬颗粒饲料的离散元模型以及标定后的硬颗粒饲料的物性参数可用于离散元仿真试验。

3 讨论

本研究通过试验测定了硬颗粒饲料的本征参数。根据密度公式求出硬颗粒饲料的密度为1.096 kg/m³，通过单轴平板压缩试验结合物理公式得出硬颗粒饲料的弹性模量、剪切模量、泊松比分别为336.32 MPa、126.65 MPa、0.33。建立了硬颗粒饲料的离散元模型，结合台架试验与仿真试验标定了硬颗粒饲料与接触材料间的接触参数，即硬颗粒饲料-钢板、硬颗粒饲料-尼龙板的碰撞恢复系数、静摩擦因数、滚动摩擦因数分别为0.380、0.538、0.028，0.405、0.356、0.033。采用圆筒提升法进行堆积角试验，基于Design-Expert 11软件对试验数据进行方差分析、响应面分析和参数优化，得出硬颗粒饲料间最优接触参数组合，并对最优参数组合进行仿真堆积角试验验证，从而确定硬颗粒饲料间的碰撞恢复系数为0.364、静摩擦因数为0.236、滚动摩擦因数为0.13。利用自主研发的投饵机进行硬颗粒饲料抛撒试验，以硬颗粒饲料分布变异系数C_V为评价指标，设置相同的试验参数进行台架试验与仿真试验，试验结果显示，实际硬颗粒饲料分布变异系数为1.092%，仿真硬颗粒饲料分布变异系数为1.073%，即台架试验与仿真试验的硬颗粒饲料分布变异系数相对误差为1.709%，表明本研究所建立的渔用硬颗粒饲料离散元模型以及所标定的硬颗粒饲料物性参数可用于离散元仿真试验。

参考文献References

于建群，付宏，李红，等.离散元法及其在农业机械工作部件研究与设计中的应用［J］.农业工程学报，2005，21（5）：1-6.YU J Q，FU H，LI H，et al.Application of discrete element method to research and design of working parts of agricultural machines［J］.Transactions of the CSAE，2005，21（5）：1-6（in Chinese with English abstract）. [百度学术]

曾智伟，马旭，曹秀龙，等.离散元法在农业工程研究中的应用现状和展望［J］.农业机械学报，2021，52（4）：1-20.ZENG Z W，MA X，CAO X L，et al.Critical review of applications of discrete element method in agricultural engineering［J］.Transactions of the CSAM，2021，52（4）：1-20（in Chinese with English abstract）. [百度学术]

顿国强，范芯蕊，杨永振，等.基于离散元的插装式大豆排种器改进设计及试验［J］.华中农业大学学报，2020，39（3）：135-144.DUN G Q，FAN X R，YANG Y Z，et al.Improved design and experiment of a plug-in soybean seed metering device based on discrete element［J］.Journal of Huazhong Agricultural University，2020，39（3）：135-144 （in Chinese with English abstract）. [百度学术]

MCCARTHY J J，KHAKHAR D V，OTTINO J M.Computational studies of granular mixing［J］.Powder technology，2000，109（1/2/3）：72-82. [百度学术]

王韦韦，蔡丹艳，谢进杰，等.玉米秸秆粉料致密成型离散元模型参数标定［J］.农业机械学报，2021，52（3）：127-134.WANG W W，CAI D Y，XIE J J，et al.Parameters calibration of discrete element model for corn stalk powder compression simulation［J］.Transactions of the CSAM，2021，52（3）：127-134（in Chinese with English abstract）. [百度学术]

MAK J，CHEN Y，SADEK M A.Determining parameters of a discrete element model for soil-tool interaction［J］.Soil and tillage research，2012，118：117-122. [百度学术]

CAMPOS A V P，OOLIVEIR A O S，FERREIR A L V，et al. Dem simulations of the dynamic angle of repose of acerola residue： a parametric study using a response surface technique［J］. Blucher chemical engineering proceedings，2015，1（2）：11326-11333. [百度学术]

COETZEE C J，ELS D N J.Calibration of discrete element parameters and the modelling of silo discharge and bucket filling［J］.Computers and electronics in agriculture，2009，65（2）：198-212. [百度学术]

焦俊，张国忠，杜俊，等.莲藕主藕体弯曲破坏离散元仿真分析［J］.华中农业大学学报，2021，40（5）：217-225.JIAO J，ZHANG G Z，DU J，et al.Discrete element simulation of bending failure of main lotus root［J］.Journal of Huazhong Agricultural University，2021，40（5）：217-225（in Chinese with English abstract）. [百度学术]

张锐，韩佃雷，吉巧丽，等.离散元模拟中沙土参数标定方法研究［J］.农业机械学报，2017，48（3）：49-56.ZHANG R，HAN D L，JI Q L，et al.Calibration methods of sandy soil parameters in simulation of discrete element method［J］.Transactions of the CSAM，2017，48（3）：49-56（in Chinese with English abstract）. [百度学术]

彭飞，王红英，方芳，等.基于注入截面法的颗粒饲料离散元模型参数标定［J］.农业机械学报，2018，49（4）：140-147.PENG F，WANG H Y，FANG F，et al.Calibration of discrete element model parameters for pellet feed based on injected section method［J］.Transactions of the CSAM，2018，49（4）：140-147（in Chinese with English abstract）. [百度学术]

MIYAMOTO T，NOGUCHI S，MATSUTOMO S.Numerical simulation and experimental verification for rice using distinct element method［J］.Communications in computer and information science，2012， 323： 472-479. [百度学术]

刘凡一，张舰，李博，等.基于堆积试验的小麦离散元参数分析及标定［J］.农业工程学报，2016，32（12）：247-253.LIU F Y，ZHANG J，LI B，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 with English abstract）. [百度学术]

余参参，段宏兵，蔡兴奎，等.基于离散元仿真参数的微型薯物料测定［J］.华中农业大学学报，2021，40（1）：210-217.YU C C，DUAN H B，CAI X K，et al.Discrete element simulation parameters-based measurement of materials for potato minituber［J］.Journal of Huazhong Agricultural University，2021，40（1）：210-217 （in Chinese with English abstract）. [百度学术]

郝建军，龙思放，李浩，等.机收麻山药离散元模型构建及其仿真参数标定［J］.农业工程学报，2019，35（20）：34-42.HAO J J，LONG S F，LI H，et al.Development of discrete element model and calibration of simulation parameters for mechanically-harvested yam［J］.Transactions of the CSAE，2019，35（20）：34-42（in Chinese with English abstract）. [百度学术]

于庆旭，刘燕，陈小兵，等.基于离散元的三七种子仿真参数标定与试验［J］.农业机械学报，2020，51（2）：123-132.YU Q X，LIU Y，CHEN X B，et al.Calibration and experiment of simulation parameters for Panax notoginseng seeds based on DEM［J］.Transactions of the CSAM，2020，51（2）：123-132（in Chinese with English abstract）. [百度学术]

吴孟宸，丛锦玲，闫琴，等.花生种子颗粒离散元仿真参数标定与试验［J］.农业工程学报，2020，36（23）：30-38.WU M C，CONG J L，YAN Q，et al.Calibration and experiments for discrete element simulation parameters of peanut seed particles［J］.Transactions of the CSAE，2020，36（23）：30-38（in Chinese with English abstract）. [百度学术]

廖宜涛，王在腾，廖庆喜，等.果荚初期饲料油菜茎秆离散元接触模型参数标定［J］.农业机械学报，2020，51（S1）：236-243.LIAO Y T，WANG Z T，LIAO Q X，et al.Calibration of discrete element model parameters of forage rape stalk at early pod stage［J］.Transactions of the CSAM，2020，51（S1）：236-243（in Chinese with English abstract）. [百度学术]

石林榕，吴建民，孙伟，等.基于离散单元法的水平圆盘式精量排种器排种仿真试验［J］.农业工程学报，2014，30（8）：40-48.SHI L R，WU J M，SUN W，et al.Simulation test for metering process of horizontal disc precision metering device based on discrete element method［J］.Transactions of the CSAE，2014，30（8）：40-48（in Chinese with English abstract）. [百度学术]

赵思琪，丁为民，赵三琴，等.基于EDEM-Fluent的气动式鱼塘投饲机性能优化［J］.农业机械学报，2019，50（9）：130-139.ZHAO S Q，DING W M，ZHAO S Q，et al.Performance analysis and optimization of pneumatic fishpond feeder based on EDEM-fluent coupling［J］.Transactions of the CSAM，2019，50（9）：130-139（in Chinese with English abstract）. [百度学术]

霍星辰.双孢菇菌种离散元分析参数与排种器优化试验研究［D］.北京：中国农业科学院，2019.HUO X C.The experimental research on discrete element analysis parameters of the Agaricus bisporus seeds and the seed-metering device optimization［D］.Beijing：Chinese Academy of Agricultural Sciences，2019（in Chinese with English abstract）. [百度学术]

联系我们

电话： 027-87287437，027-87287364
邮政编码： 430070
E-mail：hnlkxb@mail.hzau.edu.cn
地址：武汉市洪山区狮子山街1号

友情链接

期刊订阅 Email Alert 国际标准刊号：ISSN 1000-2421 国外邮发代号：BM 3816 国内刊号：CN 42-1181/S 国内发行代号：38-120