广西滨海地区甘蔗耕种土壤离散元仿真模型的参数标定

陶利民，石艳琪，黄柳益，林圣强，刘友建，李霄楠

（北部湾大学 机械与船舶海洋工程学院，广西 钦州 535011）

广西滨海地区日照充足、雨水丰盈，宜种甘蔗，海岸带土壤资源具有很好的农业开发应用前景。甘蔗属主要农业经济作物，研发适用性强、更具推广潜力的甘蔗种植机械，可促进糖料蔗生产、加快临海大工业发展进程。开沟、覆土作为甘蔗精准播种的重要环节，面临着耕种土壤回流、黏附严重、作业阻力大等难题。为更切合实际的种植环境，针对广西滨海地区甘蔗耕种土壤的特性，应用离散元法对土壤离散元仿真模型的参数进行了标定［1］。

离散元仿真模型参数的精确度取决于颗粒建模参数的准确性，主要分为本征物理参数和材料接触参数。离散元法可直观地显现土壤的力学特性和动态行为、可用于仿真分析散粒体为主的农业物料与耕种机械间的相互作用关系［2］。国内外专家基于离散元法对不同地区土壤的特性进行了深入研究，为数值模拟仿真和耦合探究奠定了前期基础。

Ucgul等［3］整合了non-linear elastic与hysteretic spring接触模型，考虑土壤塑性形变对有、无黏结力的沙滩沙增大粒径后进行参数标定。Chen等［4］分析粒径对土壤流动性的影响过程，对沙壤土模型进行有效验证。邢洁洁等［5］针对海南地区砖红土壤离散元模型，通过仿真和土槽试验验证破土阻力的一致性，选用Hertz-Mindlin with JKR接触模型，确保参数可靠。石林榕等［6］选取6种含水率的西北旱区土壤验证延迟弹性模型Hysteretic Spring Contact Mode（HSCM）和线性内聚力模型Liner Cohesion Mode（LCM）的可行性，以入土阻力为指标进行鸭嘴插入土壤和仿真试验验证了参数有效性。向伟等［7］对南方黏壤土参数进行了标定，成穴和土槽对比试验的结果表明标定准确可行。

鉴于土壤分布地区的差异性、土壤本身结构的复杂性，上述研究并不适宜本区域甘蔗种植的开沟覆土作业所需。本文以广西滨海区域的甘蔗耕种土壤为研究对象，选取“Hertz-Mindlin with JKR”接触模型，采用实测、仿真方法对相关接触参数进行标定，以便合理地设计开沟覆土装置结构，使各项性能满足实际工况的要求；为分析土壤动态扰动行为和覆土结果提供基础数据。依据土壤性质、接触特性，对土壤颗粒间、土壤与开沟、覆土装置间进行准确的参数标定显得十分必要。

1 材料与方法

采集土壤样本，测量其密度、含水率，并标定土壤物理参数。随后依据土壤质地选择土壤颗粒的接触模型，并测量土壤堆积角，为后续参数标定提供数据参考［8-9］。

1.1 材料基本参数

试验土壤取自广西钦州市钦南区久隆镇甘蔗地（108°46′11″N，22°4′56″E），在所选取5个测试点取得0~30 cm土壤层间的土壤基本数据。用烘干法将被测土壤内置于实验室恒温烘箱105 ℃条件下烘至恒重，测得含水率16.3%；采用环刀（100 cm3）与电子天秤（精度0.01 g）测得密度为1790 kg/m3。土壤和Q235号钢的其他参数参照文献［10-11］选取（表1）。

表1 离散元模型仿真参数表

1.2 土壤颗粒间接触模型选取

广西滨海地区土壤以红壤、赤红壤和砖红壤为主［8］。该地区甘蔗耕种的土壤土质疏松、透气性好、根系易吸收水分。由于不同类型土壤之间的黏结性、黏附性和土壤塑性的不同，对触土部件产生的影响也不同［9］。故接触模型的选择应考虑具有较高含水率、黏结性较强的土壤对耕种质量的动态影响。Hertz-Mindlin with Bonding黏结接触模型主要针对黏性土，通常用来模拟破碎工况，Hysteretic Spring延迟弹性模型描述了粒子的塑性变形，应用于可压缩材料。Hertz-Mindlin with JKR Cohension模型赋予土壤颗粒间的吸引凝聚力，考虑了湿颗粒之间的黏结力，故本研究选用其作为土壤离散元仿真的接触力学模型，更贴合仿真效果。

1.3 土壤堆积角物理试验

土壤堆积角试验用于评估松散土的流动性和颗粒间摩擦。取含水率为16.3%的土壤，搭建实测试验平台如图1所示。将漏斗（进料口径170 mm、漏斗出料口径17 mm）固定于铁架台一侧，为取得良好的堆积效果，多次预试验后，取漏斗出料口径与接料台间距为300 mm。开始试验之前，塞住漏斗口，放入被测土壤后，再打开漏斗口使其自然顺畅流出，待土壤静置后，从试验台不同方位平行角度拍摄静态堆积角，取其平均值。

图1 实测试验平台

在堆积角的测量试验中，土壤堆积边缘参差不一，很难使用常规角度软件实现精准测量。本次测量应用Matlab软件进行图像边缘检测处理，后利用高斯分布拟合来测量堆积角［12］。具体步骤为：先将原始图像转换成灰度图；并通过Otsu大津算法，将图像的兴趣区域进行有效分割，以此作为后续特征提取的目标对象；随后进行边缘检测，选取了比Canny算子效果更好、更符合试验要求的Zerocross算子；最后，在Origin软件中对导入的加深轮廓图像进行相对应轮廓点的提取，采用Gauss分布对图像轮廓进行拟合并取拐点处的切线与X轴的夹角为土壤堆积角。堆积角的测量均值为35.8°，过程如图2所示。

图2 土壤堆积角测量过程

2 离散元参数标定

2.1 土壤接触参数标定

将上述所测的密度、堆积角和仿真规模输入GEMM（颗粒材料模型数据库）中，以堆积角为试验指标，进行EDEM仿真标定试验。并以实测堆积角度为标准进行优化，从而确定土壤间的最优参数组合。

2.1.1 参数标定试验 在UG中建立堆积角仿真试验平台后导入EDEM中，因颗粒工厂通常被颗粒粒径限制，为合理、高效地生成所需土壤颗粒，本文在仿真中设置半径10 mm的默认球形颗粒为供试土壤进行仿真［13］。仿真试验平台如图3所示。将仿真规模、颗粒密度和堆积角等输入颗粒材料模型数据库GEMM（颗粒材料数据库）中。当JKR表面能太高、恢复系数太小、静摩擦系数太大时，土壤易因黏重而塌陷偏移［14］，因此，参数范围取值：恢复系数为0.35~0.75；静摩擦系数为0.32~0.80；滚动摩擦系数为0.05~0.15；JKR表面能为4~15 J/m2。依据上述参考值，应用Design-Expert软件的Box-Behnken优化方法设计4因素3水平共29组参数组合，选取堆积角为评价指标进行标定试验，利用CAD测量堆积角的角度。试验因素及水平如表2所示，试验结果见表3。

表2 土壤堆积角仿真试验因素及水平

图3 仿真试验平台

对表3中数据进行多元回归拟合分析，建立堆积角和变量x1、x2、x3、x4的二阶回归模型：

表3 土壤堆积角仿真试验设计及结果

回归方程（1）的决定系数R2=0.9979，校正决定系数Adj R2=0.9957，均接近于1，说明回归方程的拟合程度极佳，根据模型的方差分析结果可知（表4）：模型P＜0.0001，表明该回归模型极显著；失拟项P＞0.05，失拟不显著显示该模型拟合程度高且可靠。在4个因素给定范围内：JKR表面能x1、土壤恢复系数x2、土壤滚动摩擦系数x4对堆积角的影响极其显著，土壤静摩擦系数x3较为显著；在一次交互作用项中，x1x2、x2x4的P值＜0.05，表明此交互项对堆积角的影响尤为显著，x1x3、x1x4略微显著；试验因素的二次方项中，x22对堆积角的影响很显著，x32与x42均显著，x12不显著。

表4 堆积角回归模型的方差分析

2.1.2 参数优化及验证 应用Design-Expert软件中的优化模块Numerical，以含水率为16.3%的土壤堆积角值35.8°为目标对堆积角回归模型进行寻优，得到最优解：土壤颗粒之间的JKR表面能8.89 J/m2，恢复系数0.51、静摩擦因数0.72、滚动摩擦因数0.08。该最优解下的堆积角仿真值为36.06°，与实测值（35.8°）相对误差为0.73%。仿真试验与实测试验的形态对比如图4所示。

图4 土壤堆积角物理与仿真试验堆型

2.2 土壤与Q235号钢接触参数测定

对于土壤与触土部件材料Q235号钢间的接触参数测定：基于斜面滑动法原理进行不同含水率下静摩擦系数的测定与分析，并采用斜板滚动试验对滚动摩擦系数进行测定。

2.2.1 静摩擦系数测定验证 土壤与Q235号钢之间静摩擦系数的测定选用斜面滑动法［15］，因广西滨海土壤质地黏重、易吸湿，相比旱土、沙壤土等的颗粒特性更繁杂，其静摩擦特性也会随含水量发生变动。为探究其关系，遂配置不同含水率（10.15%、10.62%、12.03%、13.01%、14.60%、16.30%）的1 cm3的土块进行试验。含水率为16.3%的土壤与Q235钢静摩擦系数μ1为0.453。

基于不同含水率下的土壤静摩擦系数变化图可知（图5）：含水率在10.0%~16.5%范围变化时，静摩擦系数在范围在0.40~0.46之间。随着含水率的增大，土壤随黏着力的出现，摩擦系数也随倾角的增大而增大。当超过一定程度后，随着含水率的增加，静摩擦便会减少产生黏附现象，试验中未达到临界值［16］。结合毗邻北部湾的广西滨海地区甘蔗耕种土壤环境现状，为降低湿润多雨气候对田间耕种环境的影响，应避免耕种装置受土壤摩擦阻力过大而导致作业质量变差。

图5 不同含水率下的静摩擦系数与倾角的变化

2.2.2 滚动摩擦因素测定 土壤接触Q235钢板表面并相对其做纯滚动运动或产生纯滚动趋势时，即产生滚动摩擦力。参照斜板滚动试验进行滚动摩擦系数的测定［17］，搭建的斜板实测试验台由斜板和数显液晶倾角仪、测量直尺3部分组成，如图6a。将含水率16.3%的5 mm小球放置第一块斜板中央起始位置，以初速度为0 m/s沿斜面自由平稳滑落至第二块斜板直至静止，测量其水平滚动距离S，重复10次操作并记录，采用式（2）计算，得出S滚动距离的平均值21.3 mm，滚动摩擦系数μ2为0.158。

式（2）中：H2为小球下落前与水平面距离（mm）；L为小球在斜面的滚动距离（mm）；θ2为斜面钢板与水平面钢板的夹角（°）；S为小球在水平面钢板滚动的距离（mm）。

仿真试验台如图6b所示，为消除土壤之间未产生接触对试验结果造成的误差影响，将土壤间接触系数均设置为0。在预仿真试验结果的基础上，保证仿真与实测试验台的关键参数相统一，取弹性恢复系数为0.5，测定此条件下的仿真滚动距离，共进行10组试验，每组重复3次后取均值，利用EDEM自带的Ruler工具测定滚动距离，绘成折线图，由图7可知：仿真距离与实测距离的变化趋势基本相同，根据土壤滚动距离y和滚动摩擦因数x的关系，对仿真值进行曲线拟合得方程：

图6 滚动摩擦实测和仿真试验台

图7 不同滚动距离下滚动摩擦因数的实测值与仿真值

式（3）的决定系数R2=0.9604，接近于1，证明其拟合的可靠性和准确度高。实测滚动距离均值y=21.3 mm代入二次多项式拟合曲线方程（3）得x为20.68 mm，仿真中验证与实值的相对误差为2.91%，验证了仿真标定的准确性，遂取滚动摩擦因数μ2为0.158。

3 结果与分析

由于不同土层含水率存在差异性，土壤的平均含水率最大不超过30%［18］，使用上述所测得的土壤参数，构建不同含水率（1.0%、10.0%、16.3%、30.0%）的土层。在0.3 m/s的前进速度下，通过EDEM软件进行开沟器仿真试验，探究不同含水率下工作阻力的变化。由表5可知，开沟器在4种不同含水率的土壤下工作阻力范围为653.63~ 1387.12 N，当土壤含水率由1.0%增大至16.3%时，开沟器破土力随土壤的弹性模量与剪切模量的下降而减小［19-20］；而由16.3%增大至30.0%时，土壤颗粒之间的黏附力增大，结块的大颗粒对开沟器造成冲击，开沟工作阻力也随之增大。

表5 不同含水率的开沟工作阻力

从微观层面上分析表明：应用EDEM中的Clipping功能添加剪裁平面，后勾选show legend显示法向力的梯度表。由土壤颗粒受力云图可知（图8）：图8a中开沟器前端红色颗粒分布密集，表明其受到的阻力最大；由图8b和图8c得知，随着含水率由1.0%增大至16.3%的过程中，红色颗粒明显减少，表明随着含水率的增大工作阻力减小至适宜值。但伴随含水率从16.3%增大至30.0%，增多的红色颗粒反映出工作阻力因土壤黏性增强又持续增大（图8d）。因此，在种植甘蔗时，应尽量避免干旱或雨天时节，选择含水率适中的土壤进行耕种，使开沟覆土等作业满足工况需求。

图8 不同含水率下土壤颗粒受力情况

4 结论

确定广西滨海地区甘蔗耕种土壤的本证参数及接触模型，借助Matlab图像处理技术，选择Zerocross算子对土壤轮廓边缘进行处理，得土壤堆积角测量值为35.8°

通过4因素3水平正交组合试验，对土壤间的接触参数进行标定，以实测堆积角的角度35.8°为目标寻优得土壤间接触参数最佳组合：土壤—土壤恢复系数、土壤—土壤静摩擦系数、土壤—土壤滚动摩擦系数、表面能参数分别为：0.51、0.72、0.08、8.89 J/m2。通过实际试验和模拟仿真相结合的方法得出土壤与Q235号钢的接触参数：静摩擦系数0.453、滚动摩擦系数0.158。

通过开沟器仿真试验，验证了适宜含水率土壤对开沟作业工作阻力的影响。为减少土壤摩擦阻力在设计耕种部件时应注重土壤黏附问题，土壤含水量较大时进行松土以改善土壤透气性。本文所测得的土壤模型参数在后续基于离散元法的甘蔗种植机开沟覆土作业过程仿真和装置的设计优化发挥可观后效。