土壤介质的空间变异性及非线性特征，给土壤-机械相互作用模型的建立带来了很大的困难^[1]。研究者多采用基于连续介质的有限元方法研究耕作过程中的土壤动态行为，该方法将土壤视为整体，然而实际耕作中，土壤是以颗粒群的形式运动的，存在土层的分离和混合、裂缝的出现及土壤颗粒的流动，采用有限元法进行土壤运动模拟具有局限性^[2]。离散元法是由Cundall等^[3]提出的分析散体行为的数值方法，可用来模拟颗粒材料间的微观和宏观变形，也适合仿真土壤和刚性体间的相互作用。已有利用离散元法对土壤与耕作机具间的作用过程的研究^[4-6]，证实了离散元仿真能够模拟耕作过程，可对土壤颗粒在耕作过程中的受力和运动情况进行分析。

土壤颗粒间受水分和化学物质的作用存在黏附现象，简化的接触模型难以准确模拟农业生产常见的黏性土壤的力学行为。张锐^[7]将土壤颗粒之间的液桥力和黏附力引入到模型中，修正了传统的离散元模型，使模拟结果与土壤特性更贴近。Kwork等^[8]在数值分析中创新性地将沙粒作为团聚体而不是单一的独立个体，为黏性土壤的离散元仿真提供了借鉴。采用三维球元的黏连颗粒模型和考虑土壤含水的湿颗粒模型进行土壤动力学的离散元模拟将推动土壤动力学的微观机理分析，有助于土壤耕作部件设计的创新^[9]。

进行离散元仿真需要选择适合的接触模型并确定相应的接触参数，但离散元仿真需要的参数很难直接获得，需要通过参数标定来确定。本研究针对样品土壤，采用考虑颗粒间黏结力的“Hertz-Mindlin with JKR”接触模型，进行土壤堆积角仿真试验，确定样品土壤的离散元接触参数与接触模型参数，为进一步研究黏性土壤与触土部件的关系提供基础。

1 材料与方法 1.1 接触模型的选取

接触模型的选取对模拟结果的准确性有很大影响。经典的Hertz接触模型的主要缺点是仅考虑弹性变形，不涉及由于范德华力作用导致的黏结力，因此不适合模拟黏性土壤。JKR（Johnson-Kendall-Roberts）接触模型是一种黏结性颗粒接触模型，颗粒之间的相互吸引力用表面能代表^[10]。本文仿真环境为EDEM 2.6，此版本软件中引入了“Hertz-Mindlin with JKR”接触模型，该接触模型基于JKR理论建立，在Hertz理论的基础上考虑了湿颗粒间黏结力对颗粒运动规律的影响，是1个凝聚力接触模型，适用于模拟颗粒间因静电、水分等原因发生明显黏结和团聚的物料，如农作物和泥土等^[11]。为表征颗粒间凝聚力，该模型采用了JKR法向弹性接触力F_JKR来计算：

$F_\text{JKR}=\frac{4E^*}{3R^*}α^3-4α^{\frac{3}{2}}\sqrt{πγE^*}，$

式中，F_JKR为JKR法向弹性接触力（N），γ为表面张力（N·m^-1），E^*为等效弹性模量（Pa），α 为切向重叠量（m），R^*为等效接触半径（m）。等效弹性模量E^*与等效接触半径R^*定义为:

$\frac{1}{E^*}=\frac{1-v^2_1)}{E_1}+\frac{1-v^2_2}{E_2}，\\ \frac{1}{R^*}=\frac{1}{R_1}+\frac{1}{R_2}，$

式中，E₁为接触颗粒1的弹性模量（Pa），v₁为接触颗粒1的泊松比，R₁为接触颗粒1的接触半径（m），E₂为接触颗粒2的弹性模量（Pa），v₂为接触颗粒1的泊松比，R₂为接触颗粒1的接触半径（m）。

切向重叠量(α)与法向重叠量(δ)的关系如下：

$δ=\frac{α^2}{R^*}-\sqrt{\frac{４πγα}{E^*}}，$

即使颗粒间没有直接接触，该模型也提供了颗粒间相互吸引的凝聚力计算方法。颗粒间具有非0凝聚力时的最大间隙通过下式计算：

$δ_c=\frac{α^2_c}{R^*}-\sqrt{\frac{４πγα_c}{E^*}}，\\ α_c=[\frac{９πγR^{*2}}{2E^*}(\frac{3}{4}-\frac{1}{\sqrt{2}})]，$

式中，δ_c为颗粒间具有非零凝聚力时的法向最大间隙（m），α_c为颗粒间具有非零凝聚力时的切向最大间隙（m）。当δ<δ_c，模型返回0。当颗粒并非实际接触并且间隙小于δ_c时，凝聚力达到最大值，颗粒非实际接触的最大凝聚力F_pullout计算公式为：

$F_\text{pullout}=-\frac{3}{2}γR^*。$

1.2 模型参数的类型及标定

离散元仿真中需要的模型参数大致分为3类：材料本征参数，包括泊松比、剪切模量和物料密度，这是材料自身的特性参数，通常较为固定，可通过试验或者文献获得；材料接触参数，包括碰撞恢复系数、静摩擦系数和滚动摩擦系数，材质、形状、湿度等对这类参数影响较大，通常需要通过试验测定或虚拟试验进行标定；接触模型参数。某些特殊的模型还需要额外的模型参数。这些参数由于是模型化的，很难与物料特性直接换算，通常采用虚拟试验进行标定。

利用仿真软件进行一些基本的物料参数的虚拟试验（如堆积角试验），不断调整模型参数，使模拟出来的物料性质与真实情况相一致，则认为该离散元模型参数是符合实际情况的。

1.2.1 堆积角物理试验

堆积角也称为休止角，是表征颗粒物料流动、摩擦等特性的宏观参数。颗粒状物质被倾倒于水平面上堆积为锥体，堆积物的表面与水平面所成内角即为堆积角，其与颗粒密度、颗粒的表面积、颗粒形状以及该颗粒物质的摩擦系数相关，因此堆积角试验经常被用作颗粒物料的离散元参数标定^[12-14]。本文采用文献[15]的方法进行土壤堆积角试验，试验装置如图 1所示。

土壤样品采自华南农业大学试验田，土壤质地为砂壤土，土壤堆积密度为913 kg·m^-3，土壤含水率(w)为16.2%，土壤粒径分布为 0.02~2.00 mm，砂粒质量占比为65%。试验测得土壤堆积角为40.45°。

1.2.2 堆积角虚拟试验几何模型建立及初始参数设定

根据土壤堆积角的测量方法，在EDEM 2.6软件中建立堆积角虚拟试验几何模型。模型由漏斗、接料盘及土壤颗粒组成（图 2）。模型中的漏斗顶部直径为250 mm，底部直径为60 mm，高为300 mm，接料盘直径为300 mm。漏斗底部距离接料盘400 mm。土壤颗粒几何模型采用圆球模型。实际测试中，土壤粒径与物料堆积角呈正相关关系，在虚拟仿真试验中，土壤粒径的设定也对仿真结果及仿真时间具有非常显著的影响。粒径设置过小，会大量增加仿真计算时间，粒径设置过大，又会影响仿真的真实性，为平衡计算时间与仿真结果的真实性，一般在粒径设置时对其进行比例放大，使物料堆积形状在宏观上与实测试验相符。本研究对土壤粒径与堆积角测量装置的物理尺寸均进行了比例放大，将土壤颗粒半径设置为10 mm，即满足了仿真计算的性能要求，又使得堆积角仿真结果不发生明显失真。

1.2.3 材料本征参数的获得

材料本征参数通常较为固定，可通过文献[6]获得，将漏斗与底板的材质设置为玻璃，土壤材料本征参数：泊松比为0.38；剪切模量为1×10⁶ Pa；密度为1 850 kg·m^-3。玻璃材料本征参数：泊松比为0.25；剪切模量为1×10⁸ Pa；密度为2 500 kg·m^-3。

1.2.4 材料接触参数与接触模型参数的获得

EDEM通用颗粒材料数据库（Generic EDEM material model database，GEMM）由英国DEM Solution公司于2015年推出，是世界上首个离散元颗粒模型数据库，包含了数千种具有代表性的颗粒物料，如岩石、土壤、矿石等。将仿真规模、材料堆积密度和堆积角输入GEMM数据库，便可得到恢复系数、静摩擦系数、滚动摩擦系数与JKR表面能的参考值范围；将土壤堆积角试验中测得的土壤堆积密度与土壤堆积角输入GEMM数据库中，便可获得材料接触参数与接触模型参数的参考值范围。土壤-土壤接触参数与接触模型参数参考值范围如下：恢复系数为0.15~0.75；静摩擦系数为0.44~1.16；滚动摩擦系数为0.05~0.25；JKR表面能为3.5~10.5 J·m^-2。

1.3 参数标定试验设计

得到离散元接触参数与接触模型参数参考值范围后，采用Box-Behnken方法进行参数标定试验设计。试验因素为JKR表面能、恢复系数、静摩擦系数和滚动摩擦系数。试验指标为堆积角，利用软件的量角器工具，分别从锥面的x方向和y方向测量堆积角，取平均值。利用design-expert软件设计了4因素3水平共29个试验点进行响应面分析，29个试验点分为2类：一类是析因点，共24个；一类是零点为区域的中心点，零点重复5次，用于估计试验的误差。因素及水平设置如表 1所示。

2 结果与分析

试验因素组合及仿真试验结果见表 2。

2.1 堆积角回归模型

对表 2 中的试验数据进行多元回归拟合分析，得到堆积角的回归模型：

y=37.6+2.53x₁-1.81x₂+0.47x₃+0.20x₄+1.31x₁x₂-0.46x₁x₃-3.22x₁x₄+1.88x₂x₃+3.80x₂x₄-0.16x₃x₄-0.34x₁²-2.13x₂²-0.90x₃²+0.57x₄²+3.26x₁²x₂-1.10x₁²x₃+1.19x₁²x₄-1.08x₁x₂²-2.14x₁x₃²-1.30x₂²x₃+1.38x₂²x₄+2.58x₂x₃²-0.79x₁²x₂²+1.19x₁²x₃²，

回归方程的决定系数R²=0.998 4，校正决定系数adj-R²=0.989 1，说明回归方程的拟合度很好，可以用该回归方程代替真实试验点结果进行分析。模型的方差分析结果见表 3，P <0.05表示该因素对试验指标具有显著影响。由表 3可知，在给定的试验因素参考值范围内，x₁、x₂、x₃对堆积角的影响极显著，影响顺序依次为x₁、x₂、x₃；x₄对堆积角的影响不显著。试验因素的一次交互作用项中，x₃x₄对堆积角的影响不显著，其余均显著。试验因素的二次方项中，x₁²对堆积角的影响不显著，其余均显著。试验因素的二次方与一次方项的交互作用均极显著。试验因素的二次方项交互作用中，x₁²x₂²的影响不显著，x₁²x₃²的影响显著。模型的P=0.000 2<0.05，说明该模型极显著，可以根据模型对堆积角进行预测。

以上分析表明，JKR表面能、恢复系数、静摩擦系数和动摩擦系数对堆积角的影响为二次多项式，且存在复杂的一次与二次交互作用，这也解释了黏性土壤离散元参数标定试验的复杂性，对于农业耕作中常见的黏性土壤来讲，在进行离散元仿真时，很难通过少量的虚拟试验来获得理想的标定结果。借助GEMM这样的数据库工具可缩小参数的标定范围，并结合适合的试验设计方法可以减少参数标定工作量，避免标定过程中参数调整的盲目性。

2.2 JKR表面能对土壤颗粒堆积形状的影响

Hertz-Mindlin with JKR模型中用JKR表面能表征颗粒接触的黏性，因此该参数对于黏性土壤的离散元仿真具有重要作用。当恢复系数、静摩擦系数与滚动摩擦系数都为中心水平时，不同JKR表面能形成的土壤颗粒堆积状态如图 3所示。由图 3a可见，当JKR表面能为3.5 J·m^-2时，土壤颗粒间无显著的黏结现象，土壤颗粒的流动性较好，因此形成的土壤颗粒堆积锥面较为平滑。由图 3b可见，当JKR表面能为7.0 J·m^-2时，部分土壤颗粒间出现黏结现象，若干颗粒黏结在一起，形成小块的颗粒聚集体，影响了土壤颗粒的流动性，因此土壤颗粒堆积锥面出现凸起与凹陷现象。由图 3c可见，当JKR表面能为10.5 J·m^-2时，土壤颗粒间的黏结现象更加明显，小块的颗粒聚集体黏附在一起形成了大块的颗粒聚集体，土壤颗粒的流动性很差，因此土壤颗粒堆积锥面产生更大的凹陷。

由以上分析可知，JKR表面能反映了颗粒间黏结力的大小，对颗粒的流动性具有显著影响，从而决定了土壤颗粒的堆积形态。随着JKR表面能的增加,颗粒间黏结力增大，颗粒流动性下降，土壤颗粒的堆积高度和堆积角均升高。

2.3 参数优化与验证

利用design-expert的优化功能，在试验因素取值范围内以堆积角40.45°为目标对回归模型进行寻优，所得到优化解并非唯一解，而是若干组解。对这些优化解进行堆积角仿真试验验证，选取与土壤堆积角试验堆积形状近似的一组优化解，其中JKR表面能为7.91 J·m^-2，恢复系数为0.66，静摩擦系数为0.83，动摩擦系数为0.25。此优化解下的堆积角仿真结果为39.73°。优化解下的堆积角仿真试验结果与物理试验结果对比如图 4所示。图 4表明，通过对堆积角仿真试验获得的堆积角回归模型进行寻优，可得到一组优化参数，利用该优化参数获得的堆积角仿真试验结果与堆积角物理试验结果在堆积角度和堆积形状上都具有较高的相似性。这表明，对于本试验中所使用的样品土壤，利用此优化参数进行堆积角仿真试验所模拟出来的物料性质与真实情况相一致，该参数可为样品土壤与触土部件间相互作用的进一步的离散元仿真提供依据。

3 结论

Hertz-Mindlin with JKR 模型在Hertz理论的基础上考虑了湿颗粒间黏结力对颗粒运动规律的影响，以该模型为接触模型对黏性土壤进行堆积角仿真试验，实现了黏性土壤的离散元参数标定。JKR表面能参数反映了颗粒间黏结力的大小，对颗粒的流动性具有显著影响，从而决定了土壤颗粒的堆积形态。随着JKR表面能的增加,颗粒间黏结力增大，颗粒流动性下降，土壤颗粒的堆积高度和堆积角均升高。

堆积角回归模型的方差分析表明：该模型极显著，可以根据模型对堆积角进行预测，试验因素对堆积角的影响为二次多项式，且存在复杂的一次与二次交互作用，因此很难通过少量虚拟试验获得理想的标定结果。借助GEMM数据库并结合Box-Behnken试验设计方法能够缩小参数标定范围，避免了标定试验的盲目性，提高了标定效率。堆积角仿真试验与物理试验在堆积角度和堆积形状上都具有较高的相似性，本研究参数可为样品土壤与触土部件间相互作用的进一步离散元仿真提供依据。