基于响应曲面法的土壤离散元模型的参数标定研究*

刘坤宇，苏宏杰，李飞宇，焦巍

(1. 中国农业科学院草原研究所/农业农村部草原畜牧业装备科学观测实验站，呼和浩特市，010010；2. 内蒙古农业大学机电工程学院，呼和浩特市，010010)

0 引言

离散元方法是分析和求解一些复杂的离散系统的动力学问题所提出的一种数值计算方法[1]，1971年由Peter Cundall首次提出[2]，并应用于岩土力学。土壤作为一种典型的离散系统，在农业机械耕作中与土壤接触时，土壤—机械相互作用模型的建立非常困难。而在离散元中，可以将土壤看作为由大量离散的独立运动颗粒组成的整体。离散元法能够从微观的角度出发，在研究农业机械对土壤的扰动量，土壤对农业机械的作用力等研究时能够更加精确，可以直观的反应每个颗粒所受的力的大小、位移等[3-5]。

现如今离散元法在农业机械中的应用越来越广泛，李博，王燕[6-7]采用离散元方法，对深松铲进行仿真，分析深松铲在耕作中所受的阻力以及松土效果等研究，但在仿真分析过程中土壤的力学参数对仿真试验结果影响较大。

为解决仿真土壤力学参数对仿真结果的影响[8-11]，针对内蒙古呼和浩特地区沙壤土进行土壤标定。采用Hertz-Mindlin with JKR接触模型，以真实的土壤直剪切试验，堆积试验来获取土壤的泊松比，以及堆积角。以堆积角为响应值，基于响应面优化，标定了土壤离散元的参数。

采用实测试验与离散元仿真模拟相结合的方法，对研究土壤与农机具接触部件的影响有着至关重要的意义，如深松机，移栽机，旋耕机等农业装备，在微观角度可以更直接的分析土壤与农机具相互作用的运动规律。在离散元仿真过程中，土壤的离散元模型建立的准确度，对仿真结果影响巨大，所以准确地建立土壤的离散元模型有着重要的研究意义。

1 材料与方法

1.1 土壤基本参数

试验时间：2019年10月17日。试验地点：内蒙古呼和浩特市耕地试验田，采用5点取样法，取0～50 cm 土壤，将土壤带回实验室，进行相关参数测量，为后续深松机，移栽机等农业机械离散元仿真提供支持。本次试验用烘干法测得土壤含水率为11.62%，测得容重为1.36 g/m3。土壤密度为1 165 kg/m2。使用分级筛对土壤进行颗粒分级，用电子秤称取不同粒径对应质量，计算得到颗粒粒径的对应质量分数，如表1所示。

表1 颗粒粒径分布Tab. 1 Particle size distribution

1.2 土壤直剪试验

试验采用BZJ-2型应变控制式直剪仪，如图1所示，测试时对土壤施加4种不同压力，分别施加100 kPa、200 kPa、300 kPa、400 kPa的垂直压力进行剪切，每次转动手轮一圈，记录量力环的读数，使得土样在3～5 min剪破。当土样剪切结束后，取走砝码、透水石，然后清扫剪切盒里面的土样。旋转手轮使其百分表归零。每组试验重复3次，记录试验结果。

图1 BZJ-2型应变控制式直剪仪

按式(1)、式(2)计算土壤剪应力，可得剪应力(抗剪强度)测试结果如表2。

σ=P/S

(1)

式中：σ——剪切面法向垂直应力，kPa；

S——土壤的剪切截面面积，m2。

根据库伦剪切公式

τ=C+σ×tanφ

(2)

式中：τ——土壤剪应力，kPa；

C——土体粘聚力，kPa；

σ——剪切面的法向应力，kPa；

φ——土壤内摩擦角，(°)。

绘制抗剪强度与垂直应力的关系曲线如图2所示。根据土壤的剪切强度与垂直压力关系，获取土壤的内摩擦角19°，土壤的内聚力9.06，摩擦系数0.344，根据材料力学可以求得土壤的泊松比如式(3)所示[12]。

表2 剪应力(抗剪强度)测试结果Tab. 2 Shear stress (shear strength) test results

图2 抗剪强度与垂直应力的关系曲线

(3)

其中：K0=1-sinφ。

通过试验得到内摩擦角为19°，计算得到土壤泊松比为0.40。根据相关文献[13-14]与样品土壤的特性，选用剪切模量为1.2×103kPa，土壤泊松比一般为0.25～0.45之间，因此可选用实际试验计算得到的泊松比。

1.3 土壤堆积角试验

土壤的堆积角是农机具—土壤相互作用过程中影响土壤应变的重要参数，堆积角试验选用FT-104B型堆积角测定仪，如图3所示，测试时，将土壤倒入漏斗内，使土壤受重力自然降落并堆积，土壤的堆积斜面与水平面的夹角就是堆积角，重复5次试验。

图3 FT-104B型堆积角测定仪

为防止肉眼观测带来的数据误差，将堆积角利用Origin软件读取颗粒单侧图像，通过获取颗粒边界点，进行线性拟合，拟合曲线的斜率即为堆积角。具体图像拾取与拟合如图4所示。最终取得平均值为43.54°。

(a) 获取颗粒边界点

2 离散元虚拟标定

2.1 土壤颗粒模型

在Solidworks软件中建立堆积角试验装置，并将其导入EDEM中，试验装置如图5所示。

图5 仿真试验装置图

模型中，漏斗顶部直径为125 mm，底部直径为25 mm。颗粒在漏斗上方，颗粒工厂内生成颗粒，然后自由下落，经过漏斗底部，落入接料板上。颗粒的生成方式为Dynamic，总生成颗粒重量为500 g。仿真完成后，测量颗粒堆积角度。并进行响应面试验。

在离散元(EDEM 2018)模拟过程中，土壤模型的准确性是离散元仿真模拟的基础。为提高土壤模型的准确性，缩短仿真时间，采用离散元中自带的球形颗粒代替土壤颗粒形状，选用标准球形颗粒，颗粒半径设置为1 mm。设置瑞利时间步长为25%。

2.2 离散元仿真基本参数

内蒙古呼和浩特地区土壤类型为粘壤土，土壤具有一定粘性，采用EDEM中自带Hertz-Mindlin with JKR接触模型，颗粒之间相互吸引力用颗粒表面能表示[15]。JKR模型适用于颗粒间有明显粘性和团聚力的物料[16]。土壤堆积试验中，试验装置材料为钢材，设置密度为7 850 kg/m3，泊松比为0.3，剪切模量为7.0×107kPa。

2.3 Plackett-Burman试验

本次仿真使用EDEM2018版本，在GEMM数据库中输入实测土壤密度与堆积角度，获取仿真参数的范围值，土壤JKR表面能0～11.25 J/m2，土壤—土壤滚动摩擦0.1～0.2，土壤—土壤静摩擦0.32～1.16，土壤—土壤恢复系数0.15～0.75。通过查阅文献[13-14]获得以下参数范围，土壤—钢材恢复系数0.2～0.5，土壤—钢材静摩擦0.5～1.2，土壤—钢材滚动摩擦0.05～0.2，土壤泊松比0.2～0.5。

应用Design Expert软件进行Plackett-Burman试验设计，选取上述8个真实参数，与3个虚拟参数，每个参数选取高低水平，分别编码-1和+1表示，如表3所示，共进行12组试验，每组试验仿真3次，选取平均值为试验堆积角，Plackett-Burman试验设计与结果如表4所示。

表3 Plackett-Burman试验参数列表Tab. 3 List of Plackett-Burman test parameters

Plackett-Burman试验设计与结果如表4所示。运用Design Expert软件对PB试验结果进行方差分析，得到各参数的影响效果如表5所示。

表4 Plackett-Burman试验设计与结果Tab. 4 Plackett-Burman test design and results

表5 Plackett-Burman试验参数显著性分析Tab. 5 Significance analysis of Plackett-Burman test parameters

2.4 最陡爬坡试验

根据表5可知，X2土壤接触模型JKR表面能(参数A)、X3土壤—土壤恢复系数(参数B)、X4土壤—土壤静摩擦因数(参数C)对颗粒堆积角影响显著，其余参数影响不显著，因此针对参数A、参数B、参数C进行最陡爬坡试验，试验设计与结果如表6所示。

表6 最陡爬坡试验设计与结果Tab. 6 Design and results of the steepest climbing test

由表6可得，堆积角误差的趋势是由大减小然后再次增大的，在3号水平时，误差达到最小值，因此最优值区间在3号水平附近，随后选取3号水平为中心点，2号、4号水平分别选取为低、高水平进行后续的Box-Behnken试验和响应面设计。

2.5 Box-Behnken试验和响应面设计

2.5.1 Box-Behnken试验

由最陡爬坡试验结果设计Box-Behnken试验，表7为Box-Behnken试验设计与结果表。根据试验结果采用Design Expert软件建立堆积角与变量A、B、C的二阶回归方程

α=45.79+4.36A+1.37B+1.15C-0.672 5AB-

1.50AC+3.74BC-1.70A2-1.44B2+0.640 5C2

式中：α——为土壤堆积角。

表7 Box-Behnken试验设计与结果Tab. 7 Box-Behnken test design and results

对表7进行方差分析，如表8所示。

表8 Box-Behnken试验方差分析表Tab. 8 Box-Behnken test analysis of variance

由表8可知，土壤接触模型JKR表面能(参数A)、土壤—土壤恢复系数(参数B)、土壤—土壤静摩擦因数(参数C)对颗粒堆积角影响显著，交互作用项BC对堆积角影响显著。从单因素角度分析，各因素对堆积角的影响顺序：土壤接触模型JKR表面能(参数A)>土壤—土壤恢复系数(参数B) >土壤—土壤静摩擦因数(参数C)。

2.5.2 回归模型交互效应分析

根据Box-Behnken试验方差分析结果可得，土壤接触模型JKR表面能(参数A)—土壤—土壤恢复系数(参数B)的交互项、土壤接触模型JKR表面能(参数A)—土壤—土壤静摩擦因数(参数C)的交互项、土壤—土壤恢复系数(参数B)—土壤—土壤静摩擦因数(参数C)的交互项对颗粒堆积角影响显著。采用Design Expert软件绘制3个交互作用的响应曲面，如图6所示。

由图6(a)可以看出，土壤接触模型JKR表面能(参数A)与土壤—土壤恢复系数(参数B)变化对应的曲面坡度较大，(参数A)与(参数B)引起的堆积角变化较大，等高线曲率平缓，表明(参数A)与(参数C)交互作用不显著。

由图6(b)表示，土壤接触模型JKR表面能(参数A)与土壤—土壤静摩擦因数(参数C)对应的曲面坡度较大表明(参数A)与(参数C)对堆积角的影响较大，图6(b)的等高线曲率平缓，表明土壤接触模型JKR表面能(参数A)—土壤—土壤恢复系数(参数B)交互影响不显著。

(a) A和B交互作用

由图6(c)可知土壤—土壤静摩擦因数(参数C)与土壤—土壤恢复系数(参数B)对应的曲面坡度较小表明(参数C)与(参数B)变化引起的堆积角变化较小，图中的等高线显示较大曲率表明(参数C)与(参数B)交互作用影响显著。

3 最优参数组合与仿真验证

应用Design Expert软件的优化功能，以实际测得堆积角43.54°为目标，使得仿真结果最接近实际堆积角43.54°，所求得到若干组解，最终选取与实测堆积角最接近一组最优解：土壤接触模型JKR表面能为3.927 J/m2、X3土壤—土壤恢复系数为0.332、土壤—土壤静摩擦因数为0.719。如图7所示为实测堆积角与采用最优参数组合的仿真对比结果。本次仿真使用最优解获得的参数，仿真中泊松比、剪切模量为采用文中1.2直剪试验中得到的参数，进行3次重复模拟试验。得到3次仿真堆积角分别为41.15°、44.20°、41.75°，取得平均值为42.36°，相对误差2.7%。通过图7可知，本次试验获得的堆积角仿真结果与真实的试验结果在堆积角角度与堆积的形态上有很高的相似性。

(a) 堆积角仿真结果

4 结论

1) 针对内蒙古呼和浩特地区沙壤土，通过土壤堆积角实测试验，获取实际堆积角度，采用实测试验与仿真相结合的方法，使用Design Expert软件依次设计Plackett-Burman试验、最陡爬坡试验与Box-Behnken试验，筛选对堆积角影响显著的物理参数，分析影响堆积角参数的交互作用，最终运用优化功能，得到土壤的最优参数组合。将最后参数输入EDEM中仿真，得到仿真结果与实际堆积角试验对比发现，堆积角角度与堆积的形态上有很高的相似性。

2) 由Plackett-Burman试验结果可知，土壤接触模型JKR表面能(参数A)、土壤—土壤恢复系数(参数B)、土壤—土壤静摩擦因数(参数C)对颗粒堆积角影响显著，其余参数影响均不显著。由Box-Behnken试验可知，土壤—土壤静摩擦因数(参数C)与土壤—土壤恢复系数(参数B)对堆积角的交互作用影响为显著。其余交互作用均布显著。