椰糠是椰子加工过程中的副产物,是从椰子外壳纤维脱落下来的一种天然有机质媒介,椰糠因其保温、保湿、疏松、透气的优良理化性质,在设施农业中正逐步取代草炭成为应用最为广泛的栽培基质[1,2,3,4]。研究表明,单一依靠椰糠作为栽培基质是存在缺陷的,通常需要与其他物料如有机肥等按照一定比例搭配形成混合基质,才能保证栽培基质的碳氮比、养分含量、保水通气性等性状更适应作物生长[5,6]。
混合机械设备是椰糠栽培基质机械化生产的关键设备之一,其生产过程属于固体混合。因为固体混合区别于液体搅拌,不具有扩散特性,须施加外力才能强制流动,固体颗粒与混合设备（或其相关接触部件）的接触关系、颗粒运动特性以及微观作用机理等直接关系到混合设备的作业性能和工作效率,所以不同的固体混合设备没有统一的制造标准,往往根据不同的功能和使用要求单独设计,而设计专用设备,就要充分了解物料特性[7,8]。
目前,实际应用的椰糠混合基质生产设备大多为直接套用市面普通的饲料混合设备或粉粒体混合设备,生产过程中存在混合不均匀或功率不符等问题。现阶段仍缺乏对椰糠基质混合过程及混合性能的系统研究,对椰糠混合基质的研究也主要集中在基质不同配比对农作物生长特性的影响规律等方面,椰糠作为一种典型散体物料,科研院所对其开展物料离散元特性方面的研究较少,而近年来,针对土壤、颗粒肥料、农业物料等离散性物料特性方面的研究相对较多,并广泛运用了离散元法（discrete element method,简称DEM）进行分析,可对椰糠物料离散元特性的分析提供重要的借鉴与参考。在土壤方面,邢洁洁等[9]通过试验并利用离散元建模仿真确定砖红壤颗粒特性最佳组合,为研发适用于海南热区砖红壤地高性能耕作装备关键触土部件的设计优化提供支撑;Ucgul等[10]利用离散元软件中Hertze- Mindlin及Hysteretic Spring接触模型,研究分析了有黏结力和无黏结力的情况下土壤颗粒的塑性形变等问题。在颗粒肥料方面,罗帅等[11]使用离散元参数标定方法筛选蚯蚓粪基质颗粒参数,为研究蚯蚓堆肥处理及蚓粪利用相关设备奠定基础;Bangura等[12]采用离散元法对肥料进行建模通过分析排肥过程比较了螺旋槽轮和直槽轮的排肥速度与排肥均匀性。在农业物料方面,郑智旗[13]运用离散元分析粉碎后的玉米秸秆在输送装置内的运动规律,通过改变抛送板叶片倾角、旋转轴转速和喂入量等模拟研究参数改变对输送性能的影响;Ghodki等[14]建立了大豆离散元模型,通过与大豆形状、静止角等特征进行比较,获得了该模型的接触参数。因此,本文在对国内外散体物料离散元参数标定方法运用进行总结的基础上,针对椰糠物料具体实际,采用物理试验和建模仿真相结合的手段,基于离散元法分析确定椰糠的颗粒物料特性参数,为后续设计高效的椰糠栽培基质混合专用设备、优化混合结构关键部件和混合工艺参数、提高椰糠栽培基质生产效率与产品品质提供理论依据。
1 材料与方法
1.1 仿真原理及模型选取
椰糠离散元模拟中需要的物料特性参数主要分为三类：一是材料本征参数,包括颗粒大小、堆积密度、剪切模量和泊松比等,这是椰糠自身的特性参数,和外界无关;二是材料接触参数,包括碰撞恢复系数、静摩擦系数和滚动摩擦系数,这是椰糠颗粒相互接触才会发生作用的参数,与接触双方都有关系;三是接触模型参数,跟接触模型选取有关,它描述椰糠颗粒之间接触时颗粒的行为,一些特殊的接触模型需要特殊的模型参数[15,16,17]。本文选取“Hertz-Mindlin (no slip)”颗粒接触模型,该模型在力的计算方面高效准确[18,19],其基本原理如图1所示。颗粒之间作用力包括法向弹性力${{F}_{n}}$、法向阻尼力$F_{n}^{d}$、切向弹性力${{F}_{t}}$、切向阻尼力$F_{t}^{d}$、滚动摩擦力${{\tau }_{i}}$。
其中
$~~\frac{1}{{{E}^{\text{*}}}}=\frac{1-\mu _{i}^{2}}{{{E}_{i}}}+\frac{1-\mu _{j}^{2}}{{{E}_{j}}}$
$\frac{1}{{{R}^{\text{*}}}}=\frac{1}{{{r}_{i}}}+\frac{1}{{{r}_{j}}}$
式中：${{E}^{\text{*}}}$为等效弹性模量,${{\delta }_{n}}$为法向重叠量,${{R}^{\text{*}}}$为等效接触半径,${{\mu }_{i}}$、${{\mu }_{j}}$为接触颗粒的泊松比,${{E}_{i}}$、${{E}_{j}}$为接触颗粒的弹性模量,${{r}_{i}}$、${{r}_{j}}$为接触颗粒的接触半径。
$F_{n}^{d}=-2\sqrt{\frac{5}{6}}\beta \sqrt{{{S}_{n}}{{m}^{*}}}v_{n}^{\overrightarrow{rel}}$
其中
$\text{ }\!\!~\!\!\text{ }\beta =\frac{\ln e}{\sqrt{\text{l}{{\text{n}}^{2}}e+{{\pi }^{2}}}}~$
${{S}_{n}}=2{{E}^{\text{*}}}\sqrt{{{R}^{\text{*}}}{{\delta }_{n}}}$
$\frac{1}{{{m}^{\text{*}}}}=\frac{1}{{{m}_{i}}}+\frac{1}{{{m}_{j}}}$
式中：$\beta $为阻尼比,${{S}_{n}}$为法向刚度,${{m}^{\text{*}}}$为等效质量,$v_{n}^{\overrightarrow{rel}}$为相对速度的法向分量,$e$为恢复系数,${{m}_{i}}{{m}_{j}}$为接触颗粒的质量。
${{F}_{t}}=-{{S}_{t}}{{\delta }_{t}}$
其中
${{S}_{t}}=8{{G}^{*}}\sqrt{{{R}^{*}}{{\delta }_{n}}}$
式中：${{S}_{t}}$为法向刚度,${{\delta }_{t}}$为法向重叠量,${{G}^{*}}$为等效剪切模量。
$F_{t}^{d}=-2\sqrt{\frac{5}{6}}\beta \sqrt{{{S}_{t}}{{m}^{*}}}v_{t}^{\overrightarrow{rel}}$
式中：$v_{t}^{\overrightarrow{rel}}$为相对速度的切向分量。
${{\tau }_{i}}=-{{u}_{r}}{{F}_{n}}{{R}_{i}}{{\omega }_{i}}$
式中：${{u}_{r}}$为滚动摩擦系数,${{R}_{i}}$为接触点到颗粒质心的距离,${{\omega }_{i}}$为颗粒在接触点处单位角速度矢量。
1.2 试验材料与物理试验
1.2.1 试验材料选取 试验中所用的材料采用椰糠砖充分泡发后的椰糠物料,椰糠砖尺寸为200 mm×100 mm×50 mm、质量为(550±10)g、泡发量为12 L/kg,相关术语取自由中国农业机械化协会设施农业分会发布的团体标准T/CAMA 12—2019《无土栽培椰糠》[20]。
1.2.2 含水率测量 在待测的椰糠物料堆中采用五点法取样,每个取样点称取3~5 g样品,分别放入MA45水分分析测试仪的样品盘中,并确保样品在样品盘中均匀分布。根据5次样品测量的平均值,计算出待测椰糠的含水率为86.57%。
1.2.3 堆积密度测量 称取一定质量的椰糠将其沿喂料漏斗自然下落到量筒中,通过用质量除以量筒内椰糠的体积读数计算出待测椰糠的堆积密度,重复试验5次取平均值为370 kg/m3。
1.2.4 真实密度测量 称取一定质量的椰糠装入盛有一定体积水的100 mL量筒中,量筒增加的体积读数就是椰糠的真实体积,用椰糠质量除以增加的体积,重复5次试验,计算取平均值得到椰糠的真实密度为1 074 kg/m3。
1.2.5 颗粒尺寸测定 称取100 g椰糠样品,使用分级筛进行颗粒分级,用电子天平称取不同粒径范围对应的质量,重复5次试验,获得椰糠粒径质量分数分布,为仿真时椰糠颗粒建模和颗粒工厂生成提供依据,结果如表1所示。
表1 椰糠颗粒分级及对应的质量分数Tab. 1 Size classification and mass fraction of coconut coir |
| 粒径Size/mm | >5 | 2~5 | 1.25~2 | 0~1.25 |
|---|---|---|---|---|
| 质量分数/% | 3.65 | 32.80 | 60.38 | 3.17 |
1.3 颗粒建模与参数选取
1.3.2 标定参数范围确定 椰糠物料特性的研究目前尚不完善,其离散元仿真的参数尤其匮乏。EDEM软件内置颗粒材料数据库（Generic EDEM Material Model Database）中包含了数千种具有代表性的颗粒物料,如岩石、土壤和矿石等。通过选择仿真规模和堆积密度范围、输入物料堆积角,获得椰糠颗粒之间的接触参数参考值范围：恢复系数为0.15~0.75,静摩擦系数为0.20~1.16,滚动摩擦系数为0.05~0.15。考虑到椰糠作为农业物料及其作为栽培基质与土壤接近,主要参照文献[25,26,27,28,29,30]中各类物料参数取值或范围,综合考虑设定椰糠泊松比为0.1~0.5,剪切模量为1.0~10.0 MPa,椰糠与钢材恢复系数为0.1~0.7、静摩擦系数为0.20~0.90、滚动摩擦系数为0.05~0.30。通过查阅文献取钢材泊松比为0.3,剪切模量为7.9×1010 Pa,密度为7 850 kg/m3,重力加速度为9.81 m/s2。各参数标定试验因素及水平如表2所示。
表2 参数标定试验因素及水平Tab. 2 Factors and levels of parameter calibration test |
| 因素 Factors | 水平Levels | ||
|---|---|---|---|
| -1 | 0 | 1 | |
| A 椰糠泊松比 | 0.1 | 0.3 | 0.5 |
| B 椰糠剪切模量/MPa | 1.0 | 5.5 | 10.0 |
| C 椰糠-椰糠恢复系数 | 0.15 | 0.45 | 0.75 |
| D 椰糠-椰糠静摩擦系数 | 0.20 | 0.68 | 1.16 |
| E 椰糠-椰糠滚动摩擦系数 | 0.05 | 0.10 | 0.15 |
| F 椰糠-钢材恢复系数 | 0.1 | 0.4 | 0.7 |
| G 椰糠-钢材静摩擦系数 | 0.20 | 0.55 | 0.90 |
| H 椰糠-钢材滚动摩擦系数 | 0.05 | 0.18 | 0.30 |
注：L为椰糠;M为钢材。 | |
Note: L is coconut coir; M is steel. |
1.4 数据处理
采用Design Expert 10软件设计试验方案并对试验数据进行方差分析和显著性分析,因试验参数较多,故先设计Plackett-Burman试验筛选出对堆积角影响显著的试验因素,然后根据试验结果对显著性影响因素高（1）、中（0）、低（-1）3个水平设计Box-Behnken试验进行响应面分析,得到堆积角和显著性影响因素之间的回归模型,最后通过参数优化寻找最优解和试验验证完成参数标定。
2 结果与分析
2.1 Plackett-Burman 试验设计与结果
表3 Plackett-Burman试验设计及结果Tab. 3 Design and results of Plackett-Burman test |
| 序号 No. | 因素Factors | 堆积角 Repose angle/(°) | |||||||
|---|---|---|---|---|---|---|---|---|---|
| A | B | C | D | E | F | G | H | ||
| 1 | 1 | -1 | -1 | 1 | -1 | 1 | 1 | 1 | 38.21 |
| 2 | 1 | 1 | -1 | -1 | -1 | 1 | -1 | 1 | 30.82 |
| 3 | 1 | 1 | 1 | -1 | -1 | -1 | 1 | -1 | 31.94 |
| 4 | -1 | -1 | -1 | 1 | -1 | 1 | 1 | -1 | 36.98 |
| 5 | 1 | -1 | 1 | 1 | 1 | -1 | -1 | -1 | 44.34 |
| 6 | 1 | -1 | -1 | -1 | 1 | -1 | 1 | 1 | 36.12 |
| 7 | -1 | -1 | 1 | -1 | 1 | 1 | -1 | 1 | 32.02 |
| 8 | -1 | 1 | 1 | -1 | 1 | 1 | 1 | -1 | 34.67 |
| 9 | 1 | 1 | -1 | 1 | 1 | 1 | -1 | -1 | 45.31 |
| 10 | -1 | 1 | 1 | 1 | -1 | -1 | -1 | 1 | 37.18 |
| 11 | -1 | 1 | -1 | 1 | 1 | -1 | 1 | 1 | 46.94 |
| 12 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | 28.53 |
通过方差分析得到8个因素对堆积角的影响效果和显著性,如表4所示。分析可得：A、B、D、E、G对椰糠颗粒的堆积角有正效应,即堆积角随该5个因素的增大而增大;C、F、H则有负效应,即堆积角随该3个因素的增大而减小;其中对堆积角影响极显著（P<0.01）的因素是椰糠-椰糠静摩擦系数（D）和椰糠-椰糠滚动摩擦系数（E）,影响显著（P<0.05）的因素是椰糠剪切模量（B）,而其余因素影响不显著;模型的P<0.01,说明回归模型极显著,预测值${{R}^{2}}$=0.994 0,说明模型能很好地反映各因素的变化情况。
表4 Plackett-Burman 试验结果方差分析Tab. 4 Variance analysis of Plackett-Burman test results |
| 因素Factors | 平方和Sum of squares | 标准化效应Stdized effects | P值 P value | 显著性排序Significance ranking |
|---|---|---|---|---|
| 模型model | 387.36 | - | 0.003 0** | - |
| A | 6.74 | 1.53 | 0.060 6 | 4 |
| B | 11.45 | 1.99 | 0.031 4* | 3 |
| C | 3.75 | -1.14 | 0.116 2 | 6 |
| D | 231.10 | 8.93 | 0.000 4** | 1 |
| E | 110.19 | 6.17 | 0.001 3** | 2 |
| F | 5.55 | -1.38 | 0.075 9 | 5 |
| G | 2.34 | 0.90 | 0.181 8 | 7 |
| H | 0.25 | -0.29 | 0.614 0 | 8 |
注：**表示该项极显著（P<0.01）,*表示该项显著（P<0.05）。 | |
Note: ** shows extremely significant (P<0.01), *shows significant (P<0.05). |
2.2 Box-Behnken试验设计与结果
表5 Box-Behnken 试验设计及结果Tab. 5 Design and results of Box-Behnken test |
| 序号 No. | 因素Factors | 堆积角 Repose angle/(°) | ||
|---|---|---|---|---|
| B | D | E | ||
| 1 | 0 | -1 | 1 | 35.34 |
| 2 | -1 | 1 | 0 | 44.05 |
| 3 | 0 | 0 | 0 | 39.38 |
| 4 | -1 | -1 | 0 | 35.96 |
| 5 | 0 | 0 | 0 | 42.93 |
| 6 | 1 | -1 | 0 | 35.65 |
| 7 | 0 | 1 | -1 | 36.63 |
| 8 | 0 | 1 | 1 | 46.04 |
| 9 | 0 | 0 | 0 | 40.74 |
| 10 | -1 | 0 | 1 | 41.44 |
| 11 | 0 | 0 | 0 | 40.70 |
| 12 | 1 | 0 | -1 | 36.24 |
| 13 | 0 | 0 | 0 | 40.74 |
| 14 | 1 | 0 | 1 | 43.67 |
| 15 | -1 | 0 | -1 | 36.91 |
| 16 | 1 | 1 | 0 | 40.43 |
| 17 | 0 | -1 | -1 | 30.65 |
2.2.2 堆积角回归模型分析 应用Design expert软件建立堆积角与3个显著性影响因素的二阶模型：
$\begin{align} & \theta =40.90-0.3B+3.69D+3.26E-0.83BD+ \\ & \ \ \ \ \ 0.73BE+1.18DE+0.26{{B}^{2}}-2.14{{D}^{2}}-1.60{{E}^{2}} \\ \end{align}$
该回归方程的P=0.001,决定系数${{R}^{2}}$= 0.9492,校正决定系数$Adj\ {{R}^{2}}$=0.8838,信噪比为13.47,说明该堆积角回归方程极显著,相关性好,可以对目标堆积角进行预测;变异系数CV=3.43%,失拟项P=0.4007>0.05,表明该方程拟合程度好,对该模型试验结果进行回归方差分析,如表6所示。结果表明,椰糠-椰糠静摩擦系数（D）、椰糠-椰糠滚动摩擦系数（E）对堆积角影响极显著;椰糠-椰糠静摩擦系数的二次方项（${{D}^{2}}$）、椰糠-椰糠滚动摩擦系数的二次方项（${{E}^{2}}$）对堆积角影响显著。
表6 Box-Behnken试验二次回归模型方差分析Tab. 6 ANOVA of Box-Behnken design quadratic model |
| 方差源 Source of variance | 平方和 Sum of squares | 自由度df | 均方 Mean square | P值 P value |
|---|---|---|---|---|
| 模型Model | 236.72 | 9 | 26.30 | 0.001 0** |
| B | 0.7 | 1 | 0.70 | 0.553 2 |
| D | 109.15 | 1 | 109.15 | 0.000 1** |
| E | 84.89 | 1 | 84.89 | 0.000 2** |
| BD | 2.74 | 1 | 2.74 | 0.258 5 |
| BE | 2.10 | 1 | 2.10 | 0.317 0 |
| DE | 5.57 | 1 | 5.57 | 0.122 9 |
| ${{B}^{2}}$ | 0.29 | 1 | 0.29 | 0.701 2 |
| ${{D}^{2}}$ | 19.24 | 1 | 19.24 | 0.013 9* |
| ${{E}^{2}}$ | 10.72 | 1 | 10.72 | 0.045 3* |
| 残差Residual | 12.68 | 7 | 1.81 | |
| 失拟项Lack of fit | 6.15 | 3 | 2.05 | 0.400 7 |
| 纯误差Pure error | 6.52 | 4 | 1.63 | |
| 总和Sum | 249.40 | 16 |
注：**表示该项极显著（P<0.01）,*表示该项显著（P<0.05）。 | |
Note: ** shows extremely significant (P<0.01), * shows significant (P<0.05). |
2.2.3 参数优化与试验验证 运用Design expert软件优化功能,以椰糠堆积角试验值45.69°为目标值,在设置的试验因素水平范围内对堆积角回归模型进行寻优,所得优化解若干组。以优化解作为因素水平值进行堆积角仿真试验,并与物理试验堆积角进行对比验证,得到试验形状相似的一组最优解,即椰糠剪切模量（B）为1.44 MPa,椰糠-椰糠静摩擦系数为（D）为1.11,椰糠-椰糠滚动摩擦系数（E）为0.15,其余非显著性因素取中间水平。此优化解下的堆积角仿真试验结果45.25°与物理试验结果45.69°对比如图3所示,仿真结果与真实情况一致,说明该优化参数可为样品椰糠与混合设备的进一步离散元仿真提供依据。
3 讨论
（1）利用椰糠+有机肥制成混合基质替代草炭,不仅可以“变废为宝”,避免生态资源的枯竭,还可以减少环境的二次污染。采用离散元法研究机械设备与物料相互作用的研究虽然不少,在混合设备领域也有相关报道,但主要集中在饲料、颗粒肥料或者混凝土的研究方面,缺乏针对椰糠基质综合利用的混合搅拌研究。现有文献中针对椰糠物料特性的研究主要分布在结构及比例、分解度、灰分含量、孔隙度、pH、矿质营养等方面,如前所述椰糠离散元模拟中所需的本征参数、接触参数等必备特性参数缺乏数据支撑,而这些参数以试验方法测量往往难以获得。本文选择利用EDEM 2018离散元软件建立椰糠颗粒模型,预设相关参数进行椰糠堆积角的模拟仿真,通过与物理试验相结合进行参数标定,确定椰糠离散元模拟所需的准确参数,研究结果在国内外椰糠加工设备基础研究方面尚未见报道。
（2）本文通过分析文献以及结合前期试验测量的堆积角,根据EDEM 2018材料库确定椰糠泊松比、剪切模量,椰糠-椰糠恢复系数、静摩擦系数、滚动摩擦系数,椰糠-钢材恢复系数、静摩擦系数、滚动摩擦系数这8个因素的取值范围。然后进行了Plackett-Burman试验设计筛选出对椰糠堆积角影响显著的接触参数和本征参数包括：椰糠-椰糠静摩擦系数、椰糠-椰糠滚动摩擦系数和椰糠剪切模量。通过Box-Behnken试验设计建立了堆积角与显著性影响因素的回归模型,方差分析表明了椰糠-椰糠静摩擦系数和椰糠-椰糠滚动摩擦系数对堆积角影响极显著,椰糠-椰糠静摩擦系数的二次方项、椰糠-椰糠滚动摩擦系数的二次方项对堆积角影响显著。以椰糠物理试验堆积角为目标,对回归方程寻优求解,得到显著性影响因素最佳组合为椰糠剪切模量1.44 MPa,椰糠-椰糠静摩擦系数1.11,椰糠-椰糠滚动摩擦系数0.15。最后通过最优组合参数试验验证,证实了标定的参数准确、可靠,建立的模型正确、可行。
（3）离散元法可以准确高效地分析颗粒的运动情况,并具有可视化特点,本文研究结果对椰糠输送、混合等机械装备的研发提供了基础数据。但农业物料的各种特性往往与含水率有很大关系,本文取样以椰糠砖泡发量为依据,虽尽可能缩小了与实际应用中的差距,但受不同栽培基质配方特性的要求,需进一步摸清不同含水率对椰糠堆积角的影响规律,从而确定不同含水率下椰糠物料的特性参数,实现工业化生产。另外,上述试验得出：对椰糠堆积角影响显著的是椰糠-椰糠静摩擦系数、椰糠-椰糠滚动摩擦系数和椰糠剪切模量,它们都是颗粒物料特性,与钢材特性关联不大。而椰糠混合基质生产肯定至少包含2种物料,如椰糠和有机肥等,那么有必要进一步研究椰糠-有机肥的颗粒接触特性,这样才能使整个离散元仿真模拟系统更加完善。