离散元法中所用土壤参数测量及标定方法研究
庞靖１，２，林晓君１，陈松涛３，耿令新１，周浩１，金鑫１，２㊀
（１．河南科技大学农业装备工程学院，河南洛阳４７１００３；２．河南省机械
装备先进创造协同中心，河南洛阳４７１００３；３．潍柴动力股份有限公司，山东潍坊２６１０００）
摘要㊀［目的］针对离散元仿真软件涉及的土壤参数进行测量及标定㊂［方法］基于Ｈｅｒｔｚ－ＭｉｎｄｌｉｎｗｉｔｈＪＫＲ黏结模型，通过直接测量法测量土壤的固体密度㊁弹性模量和泊松比，并用堆积角和滑动摩擦角来标定土壤接触参数㊂通过中心组合试验，采用Ｄｅｓｉｇｎ－Ｅｘｐｅｒｔ８．０．６软件，以土壤休止角㊁土壤与６５Ｍｎ钢滑动摩擦角的仿真值与实测值的相对误差为优化目标进行回归分析㊂［结果］通过分析获得最优的离散元接触参数组合为土壤间恢复系数０．２８㊁静摩擦系数０．４９㊁滚动摩擦系数０．２４㊁土壤表面能０．０４Ｊ／ｍ２，土壤与６５Ｍｎ钢间恢复系数０．５９㊁静摩擦系数０．６７㊁滚动摩擦系数０．１３㊂在所测土壤参数及最优标定参数下，采用离散元仿真模拟探针入土行为，获得探针在８ｍｍ／ｓ的贯入速度下，贯入２０㊁４０㊁６０㊁８０和１００ｍｍ处仿真试验和土槽试验探针阻力相对误差分别为８．５９％㊁９．８８％㊁９．７２％㊁０．１５％㊁６．９８％，误差在可接受范围内㊂［结论］参数测量和标定方法准确可靠性，可为松软土壤
的离散元仿真提供参考㊂关键词㊀离散元法；松软土壤；参数测量；参数标定；探针贯入中图分类号㊀Ｓ１５２．９㊀㊀文献标识码㊀Ａ㊀㊀文章编号㊀０５１７－６６１１（２０２３）１８－０００６－０６
ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．０５１７－６６１１．２０２３．１８．００２㊀㊀㊀㊀㊀开放科学（资源服务）标识码（ＯＳＩＤ
）：ＳｔｕｄｙｏｎＭｅａｓｕｒｅｍｅｎｔａｎｄＣａｌｉｂｒａｔｉｏｎＭｅｔｈｏｄｓｏｆＳｏｉｌＰａｒａｍｅｔｅｒｓＵｓｅｄｉｎＤｉｓｃｒｅｔｅＥｌｅｍｅｎｔＭｅｔｈｏｄ
ＰＡＮＧＪｉｎｇ１，２，ＬＩＮＸｉａｏ⁃ｊｕｎ１，ＣＨＥＮＳｏｎｇ⁃ｔａｏ３ｅｔａｌ㊀（１．ＣｏｌｌｅｇｅｏｆＡｇｒｉｃｕｌｔｕｒａｌＥｎｇｉｎｅｅｒｉｎｇ，ＨｅｎａｎＵｎｉｖｅｒｓｉｔｙｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏ⁃ｇｙ，Ｌｕｏｙａｎｇ，Ｈｅｎａｎ４７１００３；
２．ＨｅｎａｎＭａｃｈｉｎｅｒｙａｎｄＥｑｕｉｐｍｅｎｔＡｄｖａｎｃｅｄＣｒｅａｔｉｏｎＣｏｏｒｄｉｎａｔｉｏｎＣｅｎｔｅｒ，Ｌｕｏｙａｎｇ，Ｈｅｎａｎ４７１００３；３．ＷｅｉｃｈａｉＰｏｗｅｒＣｏ．，Ｌｔｄ．，Ｗｅｉｆａｎｇ，Ｓｈａｎｄｏｎｇ２６１０００）
Ａｂｓｔｒａｃｔ㊀［Ｏｂｊｅｃｔｉｖｅ］Ｔｏｍｅａｓｕｒｅａｎｄｃａｌｉｂｒａｔｅｓｏｉｌｐａｒａｍｅｔｅｒｓｉｎｖｏｌｖｅｄｉｎｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｓｉｍｕｌａｔｉｏｎｓｏｆｔｗａｒｅ．［Ｍｅｔｈｏｄ］ＢａｓｅｄｏｎＨｅｒｔｚＭｉｎｄｌｉｎｗｉｔｈＪＫＲｂｏｎｄｍｏｄｅｌ，ｓｏｉｌｓｏｌｉｄｄｅｎｓｉｔｙ，ｅｌａｓｔｉｃｍｏｄｕｌｕｓａｎｄＰｏｉｓｓｏｎ ｓｒａｔｉｏｗｅｒｅｍｅａｓｕｒｅｄｂｙｄｉｒｅｃｔｍｅａｓｕｒｅｍｅｎｔｍｅｔｈｏｄ，ａｎｄｓｏｉｌｃｏｎｔａｃｔｐａｒａｍｅｔｅｒｓｗｅｒｅｃａｌｉｂｒａｔｅｄｂｙｓｔａｃｋｉｎｇａｎｇｌｅａｎｄｓｌｉｄｉｎｇｆｒｉｃｔｉｏｎａｎｇｌｅ．［Ｒｅｓｕｌｔ］Ｔｈｒｏｕｇｈａｎａｌｙｓｉｓ，ｔｈｅｏｐｔｉｍａｌｃｏｍｂｉｎａｔｉｏｎｏｆｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｃｏｎｔａｃｔｐａｒａｍｅｔｅｒｓｗａｓｏｂｔａｉｎｅｄａｓｆｏｌｌｏｗｓ：ｔｈｅｅｑｕａｔｉｏｎｗａｓｏｐｔｉｍｉｚｅｄｔｏｏｂｔａｉｎｔｈｅｏｐｔｉｍａｌｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｃｏｎｔａｃｔｐａｒａｍｅｔｅｒｃｏｍｂｉ⁃ｎａｔｉｏｎａｓｆｏｌｌｏｗｓ：ｔｈｅｃｏｅｆｆｉｃｉｅｎｔｏｆｒｅｓｔｉｔｕｔｉｏｎｂｅｔｗｅｅｎｓｏｉｌｐａｒｔｉｃｌｅｓｗａｓ０．２８，ｔｈｅｃｏｅｆｆｉｃｉｅｎｔｏｆｓｔａｔｉｃｆｒｉｃｔｉｏｎｗａ
ｓ０．４９，ｔｈｅｃｏｅｆｆｉｃｉｅｎｔｏｆｒｏｌｌｉｎｇｆｒｉｃｔｉｏｎｗａｓ０．２４，ｔｈｅｓｏｉｌｓｕｒｆａｃｅｅｎｅｒｇｙｗａｓ０．０４Ｊ／ｍ２，ａｎｄｔｈｅｃｏｅｆｆｉｃｉｅｎｔｏｆｒｅｓｔｉｔｕｔｉｏｎｂｅｔｗｅｅｎｓｏｉｌａｎｄ６５Ｍｎｓｔｅｅｌｉ１３１１０．５９，ｔｈｅｃｏｅｆｆｉ⁃ｃｉｅｎｔｏｆｓｔａｔｉｃｆｒｉｃｔｉｏｎｗａｓ０．６７，ｔｈｅｃｏｅｆｆｉｃｉｅｎｔｏｆｒｏｌｌｉｎｇｆｒｉｃｔｉｏｎｗａｓ０．１３．Ｕｎｄｅｒｔｈｅｃｏｍｂｉｎａｔｉｏｎｏｆｔｈｅｍｅａｓｕｒｅｄｓｏｉｌｐａｒａｍｅｔｅｒｓａｎｄｔｈｅｏｐｔｉ⁃ｍａｌｃａｌｉｂｒａｔｉｏｎｐａｒａｍｅｔｅｒｓ，ｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｓｉｍｕｌａｔｉｏｎｗａｓｕｓｅｄｔｏｓｉｍｕｌａｔｅｔｈｅｐｒｏｂｅ ｓｐｅｎｅｔｒａｔｉｏｎｂｅｈａｖｉｏｒｉｎｔｏｔｈｅｓｏｉｌ．Ｔｈｅｒｅｌａｔｉｖｅｅｒｒｏｒｓｏｆｔｈｅｐｒｏｂｅｒｅｓｉｓｔａｎｃｅｉｎｔｈｅｓｉｍｕｌａｔｉｏｎｔｅｓｔａｎｄｓｏｉｌｔａｎｋｔｅｓｔａｔ２０，４０，６０，８０ａｎｄ１００ｍｍｐｅｎｅｔｒａｔｉｏｎｓｐｅｅｄｓｏｆ８ｍｍ／ｓｗｅｒｅ８．５９％，９．８８％，９．７２％，０．１５％ａｎｄ６．９８％，ｒｅｓｐｅｃｔｉｖｅｌｙ，ｗｉｔｈａｃｃｅｐｔａｂｌｅｅｒｒｏｒｓ．［Ｃｏｎｃｌｕｓｉｏｎ］Ｔｈｅｐａｒａｍｅｔｅｒｍｅａｓｕｒｅｍｅｎｔａｎｄｃａｌｉｂｒａｔｉｏｎｍｅｔｈｏｄｄｅｓｃｒｉｂｅｄｗａｓａｃｃｕｒａｔｅａｎｄｒｅｌｉａｂｌｅ，ｗｈｉｃｈｃａｎｐｒｏｖｉｄｅａｒｅｆｅｒｅｎｃｅｆｏｒｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｓｉｍｕｌａｔｉｏｎｏｆｓｏｆｔｓｏｉｌ．
Ｋｅｙｗｏｒｄｓ㊀Ｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｍｅｔｈｏｄ；Ｓｏｆｔｓｏｉｌ；Ｐａｒａｍｅｔｅｒｍｅａｓｕｒｅｍｅｎｔ；Ｐａｒａｍｅｔｅｒｃａｌｉｂｒａｔｉｏｎ；Ｐｒｏｂｅｐｅｎｅｔｒａｔｉｏｎ
基金项目㊀国家重点研发计划项目（２０１７ＹＦＤ０７００３００）㊂
作者简介㊀庞靖（１９７７ ），男，河南洛阳人，副教授，博士，从事农业装
备作业环境与作业质量检测技术研究㊂
收稿日期㊀２０２２－０９－０６；修回日期㊀２０２３－０２－２８
㊀㊀我国农业土壤种类多样，在离散元仿真软件里要涉及很多土壤参数，而明确土壤参数是进行正确仿真的前提㊂部分科研工作者或学者没有合适的测量工具测量土壤参数，因此选择合适的土壤参数测量和标定方法可以有效解决这一问题㊂离散元是一种可以将介质整体视为若干颗粒单元集合的数值模拟方法［１］，在散落物料流动性㊁固体破碎及机器－土壤相互作用方面具有广泛应用［２
－３］
㊂由于土壤特性复杂，有
限元土壤模型准确性不高，且只能模拟土壤破坏行为，无法模拟土壤运动过程［４］，而离散元可以解决散粒之间及边界间的接触作用，利用黏连颗粒模拟生成土壤团聚体，极大提高了土壤模型的准确性［５］㊂
目前，土壤的本征参数如固体密度㊁剪切模量和泊松比
等可以通过仪器进行测量，而接触参数难以通过常规测量方法获得，因此很多学者通过离散元对物料参数进行标定的方法获得，主要包括以堆积角和滑动摩擦角为指标的参数标定法㊂冯俊小等［６］㊁刘文政等［７］㊁郝建军等［８］㊁王黎明等［９］采用离散元法分别对秸秆㊁马铃薯㊁油葵籽㊁猪粪进行了分析和标
定㊂孙景彬等［１０］以坡地黏壤土为研究对象，针对Ｈｅｒｔｚ－
ＭｉｎｄｌｉｎｗｉｔｈＪＫＲＣｏｈｅｓｉｏｎ接触模型，以土壤颗粒的仿真堆积角为响应值标定了土壤颗粒间接触参数，通过静摩擦㊁斜板及碰撞等试验得到了土壤与６５Ｍｎ钢之间静摩擦因数㊁滚动摩擦因数和恢复系数的范围；宋少龙等［１１］用Ｈｅｒｔｚ－Ｍｉｎｄｌｉｎ
（ｎｏｓｌｉｐ）作为接触模型，通过土壤堆积和滑落试验标定土壤间和土壤与６５Ｍｎ钢间的系数；张晋［１２］用烘干法测土壤含水率，用筛分法测土壤质地，针对Ｈｅｒｔｚ－Ｍｉｎｄｌｉｎｗｉｔ
ｈｂｏｎｄｉｎｇ模型，测土壤堆积角和坚实度来标定土壤之间及土壤与其他部件间的参数㊂不同的物料以及不同的土壤接触参数都存在一定的差异，标定所用的接触模型不同得到的结果也会有差距，对于坚实度特别小的松软土壤如耕后土壤，上述研究结果已经不适用㊂
为了提高离散元法针对松软土壤研究的适应性和准确性，该研究选用Ｈｅｒｔｚ－ＭｉｎｄｌｉｎｗｉｔｈＪＫＲ黏结模型作为土壤接触模型，对于该模型所涉及的土壤泊松比㊁固体密度㊁剪切模量进行参数测量；采用堆积角和滑动摩擦角试验方法，对土壤间及土壤与触土部件间的恢复系数㊁动静摩擦系数和土壤表面能进行参数标定及优化，在优化后的参数下建立土槽仿真模型，并进行贯入阻力的室内试验与仿真试验对比，验证所测参数和标定参数的准确性，为松软土壤的仿真参数设置
㊀㊀㊀
安徽农业科学，Ｊ．ＡｎｈｕｉＡｇｒｉｃ．Ｓｃｉ．２０２３，５１（１８）：６－１１
提供参考
㊂１㊀土壤本征参数测量
该研究以河南孟津林沟村（３４ʎ３９ᶄ４７ᵡＮ㊁１１２ʎ２６ᶄ０４ᵡＥ）的土壤为试验对象，通过比重计法测量，得到土壤颗粒直径小于０．０１ｍｍ的占１１％，土壤颗粒直径在０．０１ ＜０．０５ｍｍ
的为
２８％，土壤颗粒直径在０．０５ ＜１．００ｍｍ的为４９％，ȡ１．００ｍｍ的为１２％，根据土壤质地分类表得到土壤质地为砂壤土㊂由于土壤泊松比㊁剪切模量等参数受土壤含水率的影响，因此该研究所测土壤均在含水率为（１５ʃ１）％的条件下测量，其他含水率下的测量方法与此相同㊂
１．１㊀固体密度的测量㊀土壤密度分为固体密度和体积密度，体积密度含有孔隙和空气，与压实程度有关，而固体密度是土壤本身的密度，与压实和空气无关㊂仿真软件ＥＤＥＭ所涉及的土壤密度为固体密度，该研究通过体积置换法［１７］
测量
土壤的固体密度，如图１所示
㊂
图１㊀体积置换法
Ｆｉｇ．１㊀Ｖｏｌｕｍｅｄｉｓｐｌａｃｅｍｅｎｔｍｅｔｈｏｄ
㊀㊀取少量土壤并将土壤充分烘干，得到土的质量（Ｍｓ）为
１２３．２ｇ，所用容器的体积（Ｖｒ）为４００ｃｍ３㊂先将水注满容器，得到水的质量（Ｍｗ）为４１０ｇ，由此得到水的密度：
ρｗ＝
ＭｗＶｒ
（１）
式中：ρｗ为水的密度（ｇ／ｃｍ３
）；Ｍｗ为注满水的质量（ｇ）；Ｖｒ为容器体积（ｃｍ３
）㊂
将烘干的土壤放入容器中，再将容器注满水，用水置换
出烘干土壤中的间隙和空气㊂称量注满水后土壤颗粒与水的总质量（Ｍᵡ）为４８４ｇ，因此得到补充水的质量（总空隙充水质量）和体积：
Ｍᶄｗ
＝Ｍᵡ－Ｍｓ（２）Ｖᶄｗ
＝Ｍᶄｗρｗ
（３）
式中：Ｍᶄｗ为补充水的质量（ｇ）；Ｍᵡ为土壤与水的总质量（ｇ）；
Ｍｓ为土壤的质量（ｇ）；Ｖᶄｗ为补充水的体积（ｃｍ３）㊂
恒容容器的总体积与补充水的体积之差得到土壤的体积和固体密度：
Ｖｓ＝Ｖｒ－Ｖᶄｗ（４）ρｓ＝
ＭｓＶｓ
（５）
式中：Ｖｓ为土壤体积（ｃｍ３）；ρｓ为土壤固体密度（ｇ／ｃｍ３）㊂
土壤密度一般在２．６ ２．８ｇ／ｃｍ３［１３
－１４］
，将该试验所用土
壤通过体积置换法得到土壤的固体密度为２．５６６ｇ／ｃｍ３，符合基本的土壤密度㊂
１．２㊀泊松比㊀泊松比是反映材料横向正应变与轴向正应变的绝对值的比值，通过直接剪切试验（图２）可以测量土壤泊松比
㊂
图２㊀直剪试验Ｆｉｇ．２㊀Ｄｉｒｅｃｔｓｈｅａｒｔｅｓｔ
通过绘制抗剪强度与垂直压力的关系曲线图，得到土壤的内摩擦角（ϕ）为２６ʎ，通过公式（６）和（７）得到土壤泊松比为０．３６㊂
Ｋ＝１－ｓｉｎϕ（６）ｖ＝Ｋ
１＋Ｋ
（７）
式中：Ｋ为土壤测压力系数；ϕ为土壤的内摩擦角（ʎ）；ｖ为土壤泊松比㊂
１．３㊀弹性模量与剪切模量㊀弹性模量是衡量材料抵抗弹性变形能力大小的标尺㊂试验时先将土壤做成标准大小的圆
柱形土样，土样的直径（Ｄ）和高（Ｌ）均为５０ｍｍ，用万能试验机（ＤＮＳ０２－１ＫＷ）以１ｍｍ／ｓ的速度对土样施加载荷，并读取力（Ｆ）和变形（ΔＬ）的数据，直至土壤应力呈下降趋势达到
压溃效果，共做５组试验，由公式（８） （１０）计算出土壤弹性模量的平均值为８ＭＰａ，并由公式（１１）得到土的剪切模量为
２．９９ＭＰａ，试验过程如图３所示㊂
ε＝ΔＬＬ（８）Ａ＝πＤ２æèç
ö
ø
÷
２
（９）Ｅ＝Ｆ
Ａˑε
（１０）Ｇ＝
Ｅ２ˑ（１＋ｖ）
（１１）
式中：ΔＬ为土样受压后的变形量（ｍ）；Ｌ为土样原来高度（ｍ）；Ｄ为土样直径（ｍ）；Ａ为土样接触面积（ｍ２）；Ｆ为土样受到的轴向载荷（Ｎ）；Ｅ为弹性模量（Ｐａ）；Ｇ为剪切模量（Ｐａ）㊂
图３㊀土样压溃试验
Ｆｉｇ．３㊀Ｃｒｕｓｈｉｎｇｔｅｓｔｏｆｓｏｉｌｓａｍｐｌｅ
７
５１卷１８期㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀庞靖等㊀离散元法中所用土壤参数测量及标定方法研究
２㊀土壤参数仿真标定
２．１㊀土壤堆积试验㊀通过测量堆积角，可以用来标定土壤颗粒间的碰撞恢复系数㊁静摩擦系数㊁动摩擦系数和土壤表面能㊂采用漏斗测定土壤颗粒的堆积角，如图４所示，漏斗的出口直径为２７ｍｍ，出口距水平面的高度为７５ｍｍ，试验中对土壤颗粒进行５组试验，每组试验都从４个方向测量角度并
求平均值，最后得到５组试验的平均值作为休止角，通过试验测得含水率为（１５ʃ１）％时土壤休止角为
３６．４３ʎ㊂
图４㊀堆积试验Ｆｉｇ．４㊀Ｓｔａｃｋｉｎｇｔｅｓｔ
土壤颗粒直径大小不一，土壤粒径是非常小的，在仿真计算时土壤模型的尺寸一般会比真实的土壤颗粒大，根据工况可以将土壤粒径放大１０ ５０倍，在此设定土壤颗粒粒径在
０．５ １．５ｍｍ㊂
通过文献［１５－１７］以及ＥＤＥＭ里的ＧＥＭＭＷｉｚａｒｄ材料
库，设定土壤颗粒间的恢复系数Ａ（０．２ ０．６）㊁静摩擦系数Ｂ（０．３ ０．７）㊁动摩擦系数Ｃ（０．１ ０．４）㊁土壤表面能Ｄ（０．０２
０．１０Ｊ／ｍ２
）㊂最陡爬坡试验可以较快地确定因素最优值所在
区间，由表１可知，随着参数的增大，堆积角呈增大趋势，２号水平最接近试验结果，因此选择１㊁２㊁３组试验所选的水平进行中心组合试验（ＣＣＤ）㊂试验因素和水平编码如表２所示，根据试验设计，确定编码系数γ为２，中心组合试验结果如表
３所示，其中相对误差（δ）的计算方法由下式确定：
δ＝ａ－ａ１ａˑ１００％
vimeo 0006解决方法（１２）
式中：ａ为休止角的试验值；ａ１为休止角的仿真值；δ为相对误差㊂
表１㊀最陡爬坡试验方案及结果
Ｔａｂｌｅ１㊀Ｔｅｓｔｓｃｈｅｍｅａｎｄｒｅｓｕｌｔｓｏｆｓｔｅｅｐｅｓｔｃｌｉｍｂｉｎｇ
试验序号ＴｅｓｔＮｏ．ＡＢＣＤＪ／ｍ２堆积角Ｓｔａｃｋｉｎｇａｎｇｌｅʊʎ１０．２０．３０．１０．０２１６．３９２０．３０．４０．２０．０４３４．４７３０．４０．５０．３０．０６４３．０３４０．５０．６０．４０．０８４９．１２５
０．６
０．７
０．５
０．１０
５７．３１
㊀㊀应用软件ＤｅｓｉｇｎＥｘｐｅｒｔ８．０．６对试验结果进行分析，得到
二次回归的模型㊂该二次回归模型的方差分析如表４所示，
该回归模型Ｐ＜０．０００１，说明休止角相对误差与所得回归方程关系是十分显著，失拟项Ｐ＝０．４５２３＞０．０５，说明所得回归方程的非正常误差所占比例很小㊂该试验的相关系数（ｒ）为
０．９９５７，因此所得回归方程可靠度较高㊂
表２㊀中心组合试验因素水平编码
Ｔａｂｌｅ２㊀Ｆａｃｔｏｒｌｅｖｅｌｃｏｄｉｎｇｏｆｃｅｎｔｒａｌｃｏｍｂｉｎａｔｉｏｎｔｅｓｔ
编码ＣｏｄｅＡＢＣＤʊＪ／ｍ２
２０．４００．５００．３００．０６１０．３５０．４５０．２５０．０５００．３００．４００．
２００．０４－１０．２５０．３５０．１５０．０３－２
０．２０
０．３０
０．１０
０．０２
表３㊀中心组合试验及结果
Ｔａｂｌｅ３㊀Ｃｅｎｔｅｒｃｏｍｂｉｎａｔｉｏｎｔｅｓｔａｎｄｒｅｓｕｌｔｓ
试验序号
ＴｅｓｔＮｏ．
ＡＢＣＤＪ／ｍ２相对误差（δ）ʊ％１０．２５０．３５０．１５０．０３２８．７４２０．３５０．３５０．１５０．０３２７．４８３０．２５０．４５０．１５０．０３２６．６５
４０．３５０．４５０．１５０．０３２９．１８５０．２５０．３５０．２５０．０３－２．３６６０．３５０．３５０．２５０．０３－８．９９７０．２５０．４５０．２５０．０３－１．９９８０．３５０．４５０．２５０．０３－７．８４９
０．２５０．３５０．１５０．０５
２１．０７１００．３５０．３５０．１５０．０５２１．８５１１０．２５０．４５０．１５０．０５２４．１１１２０．３５０．４５０．１５０．０５２４．９２１３０．２５０．３５０．２５０．０５－４．３２１４０．３５０．３５０．２５０．０５－７．９１１５０．２５０．４５０．２５０．０５－５．８２１６０．３５０．４５０．２５０．０５－３．１２１７０．２００．４００．２００．０４７．０３１８０．４００．４００．２００．０４１０．５１１９０．３００．３００．２００．０４８．８５２００．３００．５００．２００．０４１１．６０２１０．３００．４００．１００．０４４７．９８２２０．３００．４００．３００．０４－１１．８３２３０．３００．４００．２００．０２１０．７９２４０．３００．４００．２００．０６３．１２２５０．３００．４００．２００．０４６．６８２６０．３００．４００．２００．０４７．７７２７０．３００．４００．２００．０４３．７１２８０．３００．４００．２００．０４６．１２２９０．３００．４００．２００．０４６．６４３０
０．３０
０．４０
０．２０
０．０４８．８７㊀㊀通过回归分析得到的误差回归方程为
δ＝２３４．０３－１８２．４３Ａ－３４４．４２Ｂ－７０８．３５Ｃ－１４５４．７５Ｄ＋２７２．２５ＡＢ－４０５．７５ＡＣ＋１４８８．７５ＡＤ－２２．７５ＢＣ＋１０３３．７５ＢＤ＋２５１３．７５ＣＤ＋１５３．６３Ａ２＋２９９．１３Ｂ２＋１０８４．１２Ｃ２－６９６．８８Ｄ２
８
㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀安徽农业科学㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀２０２３年
表４㊀中心组合试验回归模型方差分析
Ｔａｂｌｅ４㊀Ｖａｒｉａｎｃｅａｎａｌｙｓｉｓｏｆｒｅｇｒｅｓｓｉｏｎｍｏｄｅｌｏｆｃｅｎｔｒａｌｃｏｍｂｉｎａｔｉｏｎ
ｔｅｓｔ
方差来源Ｓｏｕｒｃｅｏｆｖａｒｉａｎｃｅ
平方和ＳＳ自由度ＤＦ均方ＭＳＦ值ＦｖａｌｕｅＰ值Ｐｖａｌｕｅ模型Ｍｏｄｅｌ５９１８．７６１４４２２．７７１２４．７３＜０．０００１Ａ７．３７１７．３７２．１７０．１６１０Ｂ２１．７７１２１．７７６．４２０．０２２９Ｃ１３０．３３１１３０．３３３８．４５＜０．０００１Ｄ２１．９９１２１．９９６．４９０．０２２３ＡＢ７．４１１７．４１２．１９０．１５９９ＡＣ１６．４６１１６．４６４．８６０．０４３６ＡＤ８．８７１８．８７２．６２０．１２６７ＢＣ０．０５１０．０５０．０２０．９０３３ＢＤ４．２７１４．２７１．２６０．２７９１ＣＤ２５．２８１２５．２８７．４６０．０１５５Ａ２４．０５１４．０５１．１９０．２９１８Ｂ２１５．３４１１５．３４４．５３０．０５０４Ｃ２２０１．４８１２０１．４８５９．４４＜０．０００１Ｄ２
０．１３１０．１３０．０４０．８４５
５
残差Ｒｅｓｉｄｕａｌ
５０．８４１５３．３９失拟项Ｍｉｓｓｉｎｇｉｔｅｍ
３５．７４
１０
３．５７
１．１８
０．４５２３
㊀㊀应用ＤｅｓｉｇｎＥｘｐｅｒｔ８．０．６软件以休止角的相对误差为目标对回归方程求解寻优，得到土壤间相关系数（ｒ）的最优值，土壤－土壤恢复系数Ａ为０．２８，土壤－土壤静摩擦系数Ｂ为０．４９，土壤－土壤动摩擦系数Ｃ为０．２４，土壤表面能Ｄ为０．０４Ｊ／ｍ２㊂
２．２㊀土壤滑落试验㊀通过土壤滑落试验（图５），测量滑动摩擦角可以标定土壤与触土部件的碰撞恢复系数㊁静摩擦系数和动摩擦系数㊂采用的钢板长２５０ｍｍ㊁宽１９０ｍｍ，在钢板一侧放置少量的土壤颗粒，斜面沿转轴缓慢旋转，当土壤由一侧滑落到另一侧时，测定土壤的滑动摩擦角㊂试验中进行５组试验，求平均值得到土壤滑动摩擦角为３２．６３ʎ㊂
㊀㊀在仿真中为了平衡仿真时间仍将土壤颗粒半径设置为０．５ １．５ｍｍ，根据文献［１８－２１］设定土壤与触土部件的恢复
系数Ｅ（０．２ ０．６）㊁静摩擦系数Ｆ（０．３ ０．７）㊁动摩擦系数Ｇ（０．０１ ０．２０）㊂由表５可知，４号水平最接近试验结果，因此选择３㊁４㊁５组试验所选的水平进行中心组合试验（ＣＣＤ）㊂
试验因素水平编码如表６所示，编码系数γ为
１．６８２㊂
图５㊀土壤滑落试验Ｆｉｇ．５㊀Ｓｏｉｌｓｌｉｄｉｎｇｔｅｓｔ
表５㊀最陡爬坡试验方案及结果
Ｔａｂｌｅ５㊀Ｔｅｓｔｓｃｈｅｍｅａｎｄｒｅｓｕｌｔｓｏｆｓｔｅｅｐｅｓｔｃｌｉｍｂｉｎｇ
试验序号ＴｅｓｔＮｏ．
ＥＦＧ滑动摩擦角Ｓｌｉｄｉｎｇｆｒｉｃｔｉｏｎａｎｇｌｅʊʎ
１０．２０．３０．０１１９．４８２０．３０．４０．０６２３．７８３０．４０．５０．１１２８．３６４０．５０．６０．１６３２．０９５
０．６
０．７
０．２１
３７．５３
㊀㊀中心组合试验结果如表７所示，其中相对误差（δ）的计算方法与休止角的计算方法一致㊂应用软件ＤｅｓｉｇｎＥｘｐｅｒｔ
８．０．６对试验结果进行分析，得到中心组合试验的回归模型㊂该回归模型的方差分析结果如表８所
示，可以看出回归模型的Ｐ＜０．０００１，说明滑动摩擦角的相对误差与所得回归方程关系是显著的；失拟项Ｐ＝０．０６０２＞０．０５，说明所得回归方程的非正常误差所占比例很小㊂该试验的相关系数（ｒ）为０．９７９０，因此所得回归方程可靠度较高㊂
㊀㊀通过回归分析得到的误差回归方程为
δ＝５１．７８＋３９０．４８Ｅ－４２３．９７Ｆ＋１６９．９１Ｇ－２４１．６７ＥＦ－
２４０．２８ＥＧ－３１２．５０ＦＧ－２０６．４４Ｅ２＋４２５．０６Ｆ２＋１３８．２４Ｇ２
表６㊀中心组合试验因素水平编码
Ｔａｂｌｅ６㊀Ｆａｃｔｏｒｌｅｖｅｌｃｏｄｉｎｇｏｆｃｅｎｔｒａｌｃｏｍｂｉｎａｔｉｏｎｔｅｓｔ
编码ＣｏｄｅＥＦＧ１．６８２０．６００．７００．２１１０．５６０．６６０．１９０
０．５０
０．６０
０．１６－１
０．４４０．５４０．１３－１．６８２
０．４０
０．５０
０．１１
㊀㊀应用ＤｅｓｉｇｎＥｘｐｅｒｔ８．０．６软件以滑动摩擦角的相对误差为目标对回归方程求解寻优，得到土壤与钢的相关系数（ｒ）的最优值，土壤－６５Ｍｎ钢恢复系数Ｅ为０．５９，土壤－６５Ｍｎ钢
静摩擦系数Ｆ为０．６７，土壤－６５Ｍｎ钢动摩擦系数Ｇ为０．１３㊂３㊀验证试验
３．１㊀堆积角和滑动摩擦角验证试验㊀通过ＤｅｓｉｇｎＥｘｐｅｒｔ８．０．６软件对土壤的接触参数及土壤与６５Ｍｎ钢的接触参数寻优，用最优结果进行仿真试验，通过５组重复的仿真试验，
得到仿真值与实际值的对比，如表９所示㊂
９
５１卷１８期㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀庞靖等㊀离散元法中所用土壤参数测量及标定方法研究
表７㊀中心组合试验及结果
Ｔａｂｌｅ７㊀Ｃｅｎｔｅｒｃｏｍｂｉｎａｔｉｏｎｔｅｓｔａｎｄｒｅｓｕｌｔｓ
试验序号ＴｅｓｔＮｏ．
ＥＦＧ相对误差（δ）ʊ％
１０．５６０．６６０．１９－５．８２２０．５６０．６６０．１３３．８６３０．５６０．５４０．１９６．９３４０．５６０．５４０．１３１２．１４５０．４４０．６６０．１９－２．７６６０．４４０．６６０．１３２．９７７０．４４０．５４０．１９４．２９８０．４４０．５４０．１３９．９９９０．４００．６００．１６１．６５１００．６００．６０．１６０．３７１１０．５００．５００．１６１７．４７１２０．５００．７００．１６－２．８２１３０．５００．６００．１１６．９３１４０．５００．６００．２１－０．０９１５０．５００．６００．１６４．７２１６０．５００．６００．１６２．５４１７０．５００．６００．１６２．５４１８０．５００．６００
．１６３．４３１９０．５００．６００．１６２．５４２０
０．５０
０．６０
０．１６
２．５４
表８㊀中心组合试验回归模型方差分析
Ｔａｂｌｅ８㊀Ｖａｒｉａｎｃｅａｎａｌｙｓｉｓｏｆｒｅｇｒｅｓｓｉｏｎｍｏｄｅｌｏｆｃｅｎｔｒａｌｃｏｍｂｉｎａｔｉｏｎ
ｔｅｓｔ
方差来源Ｓｏｕｒｃｅｏｆｖａｒｉａｎｃｅ平方和ＳＳ自由度ＤＦ均方ＭＳＦ值ＦｖａｌｕｅＰ值Ｐｖａｌｕｅ模型Ｍｏｄｅｌ５１１．０２９５６．７８２５．６１＜０．０００１Ｅ１５．２８１１５．２８６．８９０．０２５４Ｆ
１５．８５１１５．８５７．１５０．０２３３Ｇ
０．８８１０．８８０．４００．５４１７ＥＦ６．０６１６．０６２．７３０．１２９４ＥＧ１．５０１１．５００．６８０．４３０５ＦＧ
２．５３１２．５３１．１４０．３１０４Ｅ２７．７５
１７．７５
３．５０
０．０９１１Ｆ
２
３２．８５１３２．８５１４．８２０．００３２Ｇ
２
０．２２１
０．２２０．１００．７６０７残差Ｒｅｓｉｄｕａｌ２２．１７１０２．２２失拟项Ｍｉｓｓｉｎｇｉｔｅｍ
１８．１９５３．６４４．５８０．０６０２
表９㊀试验结果对比
Ｔａｂｌｅ９㊀Ｃｏｍｐａｒｉｓｏｎｏｆｔｅｓｔｒｅｓｕｌｔｓ
试验Ｔｅｓｔ
堆积角Ｓｔａｃｋｉｎｇａｎｇｌｅʊʎ滑动摩擦角Ｓｌｉｄｉｎｇｆｒｉｃｔｉｏｎａｎｇｌｅʊʎ
仿真试验Ｓｉｍｕｌａｔｉｏｎｔｅｓｔ３５．８８３３．５３室内试验Ｉｎｄｏｏｒｔｅｓｔ
３６．４３３２．６３相对误差Ｒｅｌａｔｉｖｅｅｒｒｏｒʊ％
１．５１
２．７３
㊀㊀将仿真结果与实际堆积结果进行对比，如图６所示，结果显示优化后土壤的仿真堆积角与实际堆积角的角度差距较小，表明该组仿真设置有一定的准确性㊂
３．２㊀探针贯入验证试验㊀为了进一步验证土壤参数的准确性，采用ＥＤＥＭ软件对探针贯入土壤过程进行仿真模拟，土壤颗粒半径仍设置为０．５ １．５ｍｍ，土壤的接触模型为Ｈｅｒｔｚ－ＭｉｎｄｌｉｎｗｉｔｈＪＫＲ黏结模型，仿真参数设定以测量和标定结果
为准，探针仿真模型采用ＳｏｌｉｄＷｏｒｋｓ软件创建的．ｘ＿ｔ文件直接导入，土壤为直径８０ｍｍ㊁高１５０ｍｍ的圆柱形，设置探针入土深度为１００ｍｍ，入土速度为８ｍｍ／ｓ，仿真时间为１２．５ｓ，瑞丽时间步长为１０％㊂在所测土壤参数和最优标定参数组合下的探针贯入１㊁５和１０ｃｍ的试验过程如图７所示㊂从图
７可以看出，随着探针入土深度的增加，土壤扰动范围增大，与孙文峰等［２２－２４］
的研究结果相似㊂因为探针尖头部分直径
大于探杆直径，因此土壤扰动区域基本分布在探针尖头部分㊂
㊀㊀室内试验如图８所示，所用到的设备有万能试验机
（ＤＮＳ０２－１ＫＷ）㊁探针（尖头最粗直径１４ｍｍ）㊁桶（直径
８０ｍｍ，高度４００ｍｍ）㊂将含水率为（１５ʃ１）％的土壤装入桶内，装入深度为１５０ｍｍ，不经过压实作用，体积密度约为
１．１５ｇ／ｃｍ３㊂利用万能试验机设置探针入土速度与仿真速度一致，在８ｍｍ／ｓ速度下贯入，贯入深度为１００ｍｍ，将探针的受力情况传到电脑存储，试验重复３次，分别记下探针入土深度为２０㊁４０㊁６０㊁８０和１００ｍｍ处的阻力值，最后取平均值为阻力值，试验结果见表１０所示㊂由表１０可知，随着探针入土深度的增加，探针阻力在逐渐增加，且仿真结果和室内试验结果的阻力增加幅度大致相同
㊂
图６㊀实际试验（ａ）和仿真试验（ｂ）对比
Ｆｉｇ．６㊀Ｃｏｍｐａｒｉｓｏｎｂｅｔｗｅｅｎａｃｔｕａｌｔｅｓｔ（ａ）ａｎｄｓｉｍｕｌａｔｉｏｎｔｅｓｔ（ｂ）
４㊀结论
（１）针对土壤的本征参数（固体密度㊁弹性模量㊁剪切模
量和泊松比），用体积置换法㊁直剪试验和土样压溃试验直接测量得到㊂
（２）针对土壤间的接触参数㊁表面能以及土壤与６５Ｍｎ钢的接触参数，采用ＥＤＥＭ软件进行土壤堆积仿真模拟和土壤在６５Ｍｎ钢板的滑落试验模拟㊂以堆积角和滑动摩擦角为指标，通过中心组合试验标定接触参数和表面能，获得土壤间的恢复系数０．２８㊁静摩擦系数０．４９㊁动摩擦系数０．２４和表面能０．０４Ｊ／ｍ２，土壤与６５Ｍｎ钢的恢复系数０．５９㊁静摩擦系数
０
１㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀安徽农业科学㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀２０２３年

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。