土壤—开沟部件接触的动态仿真方法的研究

曹中华，王　攀，吴先兵，李　果，陈　建

(西南大学 工程技术学院，重庆　400715)



土壤—开沟部件接触的动态仿真方法的研究

曹中华，王攀，吴先兵，李果，陈建

(西南大学 工程技术学院，重庆400715)

摘要：早期土壤—开沟部件接触的研究主要采用试验研究及传统的分析方法。近些年，计算机技术的发展突飞猛进，仿真技术也得到了很大的进步。为此，主要分析了FEM、DEM、无网格法、CFD、ANN等5种数值模拟方法及其中两种方法的耦合在土壤—开沟部件接触的动态仿真中的应用以及各种方法在应用中存在的差异，并对存在的一些问题进行了分析。最后，提出了未来发展的重点和发展趋势。

关键词：仿真技术；开沟部件—土壤接触系统；差异；发展趋势

0引言

研究土壤—开沟部件接触的目的是为了获得开沟部件最优的工作参数且探索土壤—开沟部件之间存在的关系，最终达到减小牵引力以减小能耗的目的。假设在耕作过程中，自变量—土壤初始状态为Si、开沟部件的形状为Ts、部件移动方式为Tm、因变量—土壤最终状态为Se、耕作阻力为F，Gill和Vanden Berg[1]把这5者之间的关系用两个函数f和g表示为

F=f(Si,Ts,Tm)

Se=g(Si,Ts,Tm)

根据上式，选取Si、Ts、Tm使得F最小，同时得到满足农艺要求的Se，即设计耕作部件的问题等价于求解方程的最优解。

将上述模型具体化，是一个非常困难的问题，也一直是研究者努力的目标。虽然函数定性说明了影响土壤最终状态的相关因素，但不能描述土壤—开沟部件的接触过程。为此，本文对目前各种数值模拟方法在土壤—开沟部件接触的动态仿真方面应用的特点进行了系统介绍，并分析了各种数值模拟方法在仿真中存在的差异及适用性，找出一种理论上更适合对土壤—开沟部件接触的过程进行动态仿真的方法，以期望对开沟部件的优化提供理论依据，为土壤动力学的研究产生积极的影响。

1土壤—开沟部件接触的研究方法

采用试验、数值模拟两种方法来研究上述两个函数关系：试验包括土槽试验及田间试验；数值模拟方法主要包括传统分析方法、有限单元法、离散单元法、无网格法、计算流体力学及人工神经网络。

1.1试验、传统分析方法的应用

早期的研究主要采用了试验及传统分析方法，依据下列两条技术路线进行：①通过对田间试验的数据进行理论修正,找到能对土壤状态进行描述的特定参数及它们之间存在的特定关系;②以试验分析为基础，假设土壤—开沟部件接触符合一定关系，进一步计算出能够描述开沟过程的理论模型及推求它们之间的相互作用[2]。

传统分析方法的基本思想是假设土体为刚塑性体，通过静力平衡原理求解以自由面、开沟部件与土体的相互作用边界及土壤切削试验中找出的土体破裂线(面)组成的力隔离体,求得土体对开沟部件的阻力[3]。

1918年，E.A. White 在康奈尔大学完成的博士论文“关于犁体及其对土垡作用的研究”被认为是一个里程碑，标志着土壤耕作工作的研究开始了较理论性的探讨[4]。20世纪20年代，前苏联郭辽契金理论上对土壤—开沟部件的接触进行了分析，从而为犁体曲面的优化设计、开沟过程的动态分析奠定了基础。自1925年后， Nichols及他的同事做了大量的研究工作，并且发表了一系列的文章；英国、德国的科学家也进行了相关研究[1]。我国陈秉聪、邵耀坚、曾德超和罗锡文等人在此领域也做了大量卓有成效的研究工作[5-6]。

然而，传统分析方法有其明显的局限性：它不是基于随刀具入土角度、土壤含水量、土壤密度等因素而变化的土壤实际失效模式，而土壤的失效取决于刀具的形状、耕作速度及土壤的物理特性[7]。因而，基于假设基础上被动土压力理论及土壤失效模式的传统分析方法在耕作机具的优化方面作用有限[8]。

1977年以来，尤其是20世纪90年代以来，随着试验手段不断完善与计算机性能迅速提高，特别是由于有限单元法、离散单元法等数值模拟方法的引入，土壤动力学领域的研究取得了较大的进展。

在试验方法方面，田间试验及土槽试验继续被使用，而且一些新的方法陆续也被引入：如X射线透射、扫描电镜及各种先进的测力装置，通过这些装置来观察并记录土壤最终状态及开沟部件的响应[9-10]。

1.2有限单元法(Finite Element Method，FEM)的应用

有限单元法是用较简单的问题代替复杂问题后再求解。有限单元法的基本思想是将物体离散为有限个简单单元的组合，原物体可以通过这些单元的组合来模拟。对于土壤开沟问题，有限单元法的应用非常广泛，因为开沟部件的切向和法向的应力分布以及土壤的应力分布、应变情况、破坏的具体位置可以被较准确的预测[11]。

Mouazen等[12]认为用有限元法研究土壤—开沟部件接触问题时应当考虑以下几点：土壤—开沟部件接触的相互作用；土壤机械特性参数的选择；土壤的失效形式。郭志军[13]等按照建立土壤—开沟部件接触的有限元模型、土壤的种类的选取，及土壤本构关系的确立、室内试验和田间试验的对比分析、有限元模型参数的选取及分析内容的部分研究成果进行了归纳。根据是否考虑耕作部件对土壤的动力影响，徐中华等[14]将研究分为准静态和动态仿真两个部分：前者因耕作速度较低，采用准静态模拟；后者因耕作速度较高在分析模型中加入了动力成分。用有限元方法分析土壤—开沟部件接触相互作用的前提是要建立精确的土壤本构关系。目前常采用的土壤本构关系有Drucker—Prager非线性弹塑性模型、Duncan—chang双曲线应力—应变模型，线性弹塑性本构模型及粘弹塑性帽盖模型。

2002年，Abdul[15]等以8节点线性等参数单元建立了4种不同角度组合的深耕铲的三维非线性有限元模型，同时考虑模型几何、物理、边界非线性，以及牵引力、土壤最终状态及深松铲表面压力分布等问题，对土壤模型选用Drucker—Prager非线性弹塑性模型，并且选用砂壤土进行室内土槽试验来与有限元模型进行对比分析。2003年，Mootaz Abo—Elnor[16]提前定义了土壤的失效面，土壤本构关系使用Dancan-Chang双曲线应力—应变模型，用有限元软件Abaqus建立了土壤—开沟部件的三维有限元模型并进行了开沟过程的动态仿真，分析开沟速度及加速度对开沟阻力的影响。结果表明：开沟阻力受开沟加速度的影响显著，开沟阻力受开沟速度的影响不大；且随着加速度的增加，开沟阻力也逐渐增大。陆怀民[17]在不考虑土壤与刀面间的摩擦的基础上，采用粘弹塑性帽盖模型，建立了粘性土壤—工程机械切削部件接触的有限元模型，并进行了切削过程的仿真分析，仿真与试验均表明：在切刀水平位移相同时，切削速度越高，土壤的反作用力就越大。因此，对于耕作机械，在动力不足时，可以适当降低速度以减少切削阻力。2014年，佟金[18]等利用有限元方法，土壤模型为Drucker-Prager非线性弹塑性模型，在Abaqus软件中建立了土壤—凸齿镇压器接触的三维有限元模型，并对其运动过程进行动态仿真。本文使用任意拉格朗日一欧拉方法对模型进行自适应网格划分，以解决土体局部变形引起单元畸变而导致分析中断的问题和单元扭曲导致分析不收敛的问题。根据凸齿镇压器的两种工作模式，对模型设置不同的边界条件，探讨不同载荷、凸齿镇压器沉降量及牵引力之间存在的因果关系。为了优化凸齿镇压器的表面结构，改进其几何结构、选择合理的作业参数以提高土壤表面的微形貌从而使得机具工作简便、能耗降低。本文通过土槽试验来验证有限元仿真的准确性，搭建了基于室内土槽的凸齿镇压器牵引试验平台，通过土槽试验对有限元分析结果的有效性进行验证，以实现改进凸齿镇压器的几何结构、表面结构的优化及作业参数的选择，达到提高土壤表面微形貌加工作业质量、降低能耗的目的。2015年，Tagar等[19]以Drucker-Prager非线性弹塑性模型为土壤本构模型，对土壤在土槽和田间的破碎模式进行了有限元仿真,结果表明：有限元仿真模拟土壤破碎是一个很好的工具,土槽的仿真结果比在田间的更符合实际，无论是在土槽还是在田间，土壤的破碎和土壤的稠度有直接关系。

相比于有限元法，离散元法计算比较简单，数据量较少，无需建立大型的刚度矩阵，允许单元发生大的平移、转动，可以求解一些含有复杂物理力学模型的非线性问题；而且离散单元法适合于界面弱连接的非连续介质问题或连续体到非连续体转化的材料损伤破坏问题。

1.3离散单元法(Distinct Element Method，DEM)的应用

离散单元法DEM是一种显式求解的数值模拟方法，在最基本的牛顿第二定律基础之上，它将土壤认定为不连续的独立质点的集合[20]，具有牢固的理论依据。离散体单元法将单元分为圆盘单元和多边体单元两类，为了研究方便且考虑土壤颗粒形状的特殊性，目前研究者将用于分析土壤的离散单元法中的单元普遍视为圆盘形[21-22]。

2001年，Horner等[23]建立了一种能够反映犁与土壤之间相互作用的3D离散单元法模型。该模型具有较好的定性结果，但不能反映出定量结果。2002年，Owen等[24]提出了离散单元法适用于土壤之间相互作用的柔性或刚性建模。土壤的非线性行为和土壤之间的交互作用可以很容易地用离散单元法描述出来。2003年，徐泳[25]等提出对土壤—耕作部件接触的相互作用采用离散元法进行仿真分析。其根据土壤的力学特性，建立了合理的土壤颗粒—耕作部件接触的模型，选择合理的算法，并利用高速计算机进行仿真。结果表明：土壤的运动特征、受力、变形及能量耗散都有特定的规律，因此耕作过程的细观力学机理被读者了解。2007年，Asaf等[26]基于一种原位测试和逆向技术提出了一种确定离散单元法模型参数的方法。2013年，Korne′l等[27]开发了土壤交互作用的三维离散单元法仿真模型，并通过模拟和试验的对比证明了该模型的准确性。2014年，Elvis等[28]用离散单元法创建了土壤不可逆耕作的预测模型。

目前，在采用有限元法、离散元法处理时遇到了问题：①高速撞击往往引起结构的几何畸变(如侵彻、穿透)，采用传统有网格的数值方法，无法避免网格畸变问题，给分析计算带来困难；②材料成型过程中，采用网格模拟材料流动变形带来的不便；③动态裂纹扩展，以及“激波”现象的数值模拟，采用网格方法必须实施网格的重新生成，给问题的处理造成了困难；④流体—结构相互作用、相变问题的分析中界面的追踪模拟；⑤求解某些特殊问题(如奇异性问题、高振荡问题等)，传统的数值方法计算精度和收敛性差，且实施不便；⑥网格类的自适应算法与技术的实现比较困难，而无网格的自适应实现相比要简单便利。

1.4无网格法(Meshless Method)的应用

无网格法克服了对网格的依赖，在涉及网格畸变、移动等问题中显示了明显的优势,包括光滑粒子动力学法(Smooth Particle Hydrodynamics, SPH)及多象限法等，目前应用较多的是SPH。进行土壤的耕作时，土壤会发生破坏以至大变形，所以可以选用SPH法，因为SPH法不适用网格，就不会有网格畸变问题的发生。同时，SPH法也适用于材料在高加载速率下的断裂等问题的研究。因此，SPH非常适于耕作过程的动力学描述。

Bui等[29-30]在2003-2006年间进行了一系列关于土体中固-液两相和固-液-气三相相互作用的SPH数值模拟研究分析，提出了SPH模型。土体长时间搁置就会板结，2006年，钟江[31]等利用SPH方法建立了土壤高速凿切的有限元数值仿真模型，研究反旋凿切土壤的动态耕作过程，为更好地疏松土壤、降低能耗提供了理论基础。2010年，Binesh等[32]第一次在增强土壤弹塑性分析中应用了无网格法，通过实验验证了该方法得出的结果与有限元分析法的结果吻合。1.5计算流体力学(Computational Fluid Dynamics, CFD)的应用

Karmakar等[33]将土壤当作粘塑性材料，使用CFD对耕作刀具周围土壤变形模式进行了研究，确定了耕作刀具周围土壤流动具有“塑性流动”和“塞流”两种模式，同时指出了土壤失效前耕作刀具的临界速度和深度。同年，Karmakar等[34]从流体流动角度，构建了土壤—开沟部件接触的动力学模型，通过对接触系统的分析可以为开沟刀具的优化提动理论基础，可以实时了解沟垄平整度、压力分布及能量需求等的情况；基于土壤的粘-塑性参数，可以建立不同类型的土壤模型。在后续的研究中，Karmakar等[35]对耕作刀具在高速耕作过程中，刀具表面的压力分布和土应力模式进行了研究，指出虽然矩形刀具是非流线体，但它比流线体刀具更能抵抗土壤的压力和粘性阻力。2009年，Karmakar等[36]在不同的耕作速度、耕作深度、土壤粘度下进行了土槽试验，得到的实验数据与利用CFD得到的模拟数据进行了比较，结果显示：当耕作刀具耕作较浅耕作深度和耕作速度较低时，两者的数据相关系数最大。2008年，Barker[37]运用计算流体力学方法模拟土壤和耕作刀具交互运行的稳态和瞬态，可预测土壤与耕作刀具接触面的负荷。

1.6人工神经网络(Artificial Neural Network，ANN)的应用

1999年，Zhang[38]等应用ANN来模拟土壤-机具的相互作用及土壤行为。2001年，Jayasuriya[39]等认为ANN能处理土壤动力学中模糊和非均匀变量的输入问题。

1.7几种数值模拟耦合方法的应用

2012年，Mabssout等[40]将应用于土壤动力学问题分析的两种方法：Taylor-SPH与Runge-Kutta进行了比较，并指出Taylor-SPH方法稳定性、效率更高。

2主要问题和未来发展趋势

2.1存在的问题

研究者们对土壤—开沟部件接触的仿真分析进行了大量的研究和改进，但是至今没有一种本构关系能适应于所有的土壤状况。主要问题是：①没有考虑土壤内部出现裂纹的模拟研究；②模拟边界条件都是人为确定的，这样就会带来计算的误差；③对于高速切削(切削速度接近或达到土壤塑性应力波传播速度的土壤切削)，目前的研究还不成熟；④对于保证开沟部件的优化设计降低剪切破碎阻力的研究较少。这些都是当前土壤—开沟部件接触系统存在的主要问题。

2.2研究重点及发展趋势

基于以上土壤—开沟部件接触系统存在的问题，目前急需解决的重点是：①考虑精确的数学模型、选择合理的算法、精确的模拟边界条件和施加载荷以建立土壤—开沟部件接触的动态仿真模型，处理土壤的大变形问题。②采用合适的破碎模型，完整地分析土壤的破碎过程，对土体的内部应力和位移进行一定的预测。③土壤种类多样，通过试验研究结果和数值模拟软件分析结果的仿真对比来寻找一种本构关系可以适合于多种土壤状况的本构模型。④通过大量的实验和计算建立耕作部件合理的优化数学模型来寻找最合适的开沟部件形状，以达到对实践的指导。

3结论

通过研究FEM、DEM、无网格法、CFD、ANN这5种数值模拟方法以及其中两种方法的耦合在土壤—开沟部件接触系统的动态仿真中的应用和差异，说明了数值模拟结果的准确性，并且数值模拟花费小、周期短，有其特定的优势，应用前景广泛。

参考文献：

[1]Gill WR,Vanden Berg G E.Soil dynamics in tillage and traction[C]//Agriculture Handbook, United States Department of Agriculture, 1967.

[2]Araya K, Gao R. A non-linear three-dimensional finite element analysis ofsubsoiler cutting with pressurized air injection[J].Agric. Engng Res, 1995, 61: 115-128.

[3]姚禹肃.应用于耕具阻力预测中的力学方法[M].北京:科学出版社,1994:116-122.

[4]华云龙,张伟,董务民.力学可以为农业现代化作贡献[J].力学进展,1998, 28(3): 289-298.

[5]曾德超.机械土壤动力学[M].北京:科学技术出版社,1995.

[6]罗锡文.地面机器系统力学[M].北京:中国农业出版社, 2002.

[7]Karmakar S, Kushwaha R L.Dynamic modeling of soil-tool interaction: An overview from a fluid flow perspective[J].Terramechanics, 2006,43: 411-425.

[8]Kushwaha R L,Chi L, Shen J.Analytical and numerical models for predicting soil forces on narrow tillage tools-a review[J].Can Agric Eng，1993, 35(3):183-190.

[9]Salokhe V M, Pathak B K. Effect of aspect ratio on soil reactions on flat tines in dry sand[J].Agric. Engng Res,1993,56:179-188.

[10]Fiekle J M. Finite element modelling of the interaction of the cutting edge of tillage implements with soil[J].Agric. Engng Res, 1999, 74: 91-101.

[11]刘辉.基于LS-DYNA的旱地旋耕刀工作机理研究[D].重庆：西南大学，2012.

[12]Mouazen A M, Neményi M. A review of the finite element modeling techniques of soil tillage [J].Mathematics Computers in Simulation, 1998, 48(1): 23-32.

[13]郭志军,佟金,任露泉,等.耕作部件—土壤接触问题研究方法分析[J].农业机械学报, 2001, 32(4): 102-104.

[14]徐中华,王建华.有限元分析土壤切削问题的研究进展[J].农业机械学报,2005, 36(1): 134-137.

[15]Mouazen A M，Ramon H. A numerical statically hybrid scheme for evaluation of draught requirements of a subsoiler cutting a sandy loam soil，as affected by moisture content，bulk density and depth[J].Soil and Tillage Research，2002，63(1)：155-165.

[16]Abo-Elnor M，Hamilton R，Boyleb J T.3D Dynamic analysis of soil-tool interaction using the finite element method[J].Journal of Terramechanics，2003,40(1)：51-62.

[17]陆怀民.切土部件与土壤相互作用的粘弹塑性有限元分析[J].土木工程学报，2002,35(6)：1-4.

[18]佟金,张智泓,陈东辉,等.凸齿镇压器与土壤相互作用的三维动态有限元分析[J].农业工程学报,2014,30(10):48-59.

[19]Tagar A A, Ji C Y, Jan A, et al. Finite element simulation of soil failure patterns under soil bin and field testing conditions[J].Soil & Tillage Research，2015，145:157-170.

[20]Mouazen A M，Nemenyi M. Tillage tool design by the finite element method: Part1. Finite element modeling of soil plastic behavior[J].J.Agric Engng.Res，1999,72(1)：37-51.

[21]Hiroaki Tanaka, Koji Inooku, Osamu Sumikawa, et al.Simulation of soil behavior at subsoiling by the distinct element method[C]//Proceeding s of the 6th Asia-Pacific Conference of the International Society of Terrain-Vehicle Systems,2001: 194-200.

[22]Masato shi Momozu, Akira Oida, Hito shi Nakashima. Simulation of shear box text by the distinct element method[C]//Proceedings of the 6th Asia-Pacific Conference of the International Society of Terrain-Vehicle Systems,2001:181-188.

[23]Horner D A, Peters J F, Carrillo A. Large scale discrete element modeling of vehicle-soil interaction[J].Eng. Mech，2001, 127(10):1027-1032.

[24]Owen D, Feng R J, Souza Y, et al. The modeling of multi-fracturing solids and particulate media[C]//In: Proceedings of the 5th World Congress on Computational Mechanics,Vienna,Austria，2002.

[25]徐泳,李红艳,黄文彬.耕作土壤动力学的三维离散元建模和仿真方案策划[J].农业工程学报,2003,19(2):34-38.

[26]Asaf Z, Rubinstein D, Shmulevich I. Determination of discrete element model parameters required for soil tillage[J].Soil Tillage Res，2007, 92: 227-242.

[27]Korne′l Tama′s, Istva′n J.Jo′ri, Abdul M. Modelling soil-sweep interaction with discrete element method[J].Soil & Tillage Research, 2013, 134: 223-231.

[28]Elvis López Bravo, Engelbert Tijskens, Miguel Herrera Suárez,et al. Prediction model for non-inversion soil tillage implemented on discrete element method [J]. Computers and Electronics in Agriculture, 2014, 106: 120-127.

[29]Bui H H, Sako K, Fukagawa R. Numerical Simulation of Soil-Water Interaction Using Smoothed Paticle Hydrodynamics(SPH)Method [C]//Proceedings of the 15th International Conference of the ISTVS,Hayama, Japan，2005(9):1-15.

[30]Bui H H, Fukagawa R. Smoothed paticle hydrodynamics for soil mechanics[C]//The proceedings of 6th European Conference on numerical methods in geotechnical engineering,London.Taylor&Francis Group, 2006 :275-281.

[31]钟江，蒋建东.基于光滑粒子流体动力学仿真的板结土壤深旋耕技术[J].机械工程学报，2010,46(19):63-69.

[32]Binesh S N, N Hataf, A Ghahramani. Elasto-plastic analysis of reinforced soils using mesh-free method[J].Applied Mathematics and Computation, 2010, 215: 4406-4421.

[33]Karmakar S, Kushwaha R L. Simulation of soil deformation around a tillage tool using computational fluid dynamics[J].Transactions of ASABE, 2005, 48(3):923-932.

[34]Karmakar S, Kushwaha R L. Dynamic modeling of soil-tool interaction:An overview from a fluid flow perspective [J].Journal of Terramechanics, 2006, 43:411-425.

[35]Karmakar S, Kushwaha R L, Lague C.Numerical modeling of soil stress and pressure distribution on a flat tillage tool using computational fluid dynamics[J].Biosystems Engineering, 2007, 97: 407-414.

[36]Karmakara S R, Ashrafizadeh R L, Kushwaha. Experimental validation of computational fluid dynamics modeling for narrow tillage tool draft[J].Journal of Terramechanics,2009,46: 277-283.

[37]Barker M E. Predicting Loads on Ground Engaging Tillage Tools Using Computational Fluid Dynamics[D].Iowa:Iowa State University,2008.

[38]Zhang Z X, Kushawaha R L. Application of neural networks to simulate soil-tool interaction and soil behavior[J].Can Agric Eng, 1999, 41(2):119-205.

[39]Jayasuriya H P W, Salokhe V M. A review of soil-tine models for a range of soil conditions [J].Agric Eng Res，2001, 79(1):1-13.

[40]Mabssout M, Herreros M I， Runge-Kutta，et al.Two time integration schemes for SPH with application to Soil Dynamics[J].Applied Mathematical Modelling, 2013, 37: 3541-3563.

Research of Dynamic Simulation Method on Soil-ditching Component Interaction

Cao Zhonghua , Wang Pan, Wu Xianbing, Li Guo, Chen Jian

(College of Engineering and Technology，Southwest University，Chongqing 400715,China)

Abstract：Early research on soil-ditching component interaction mainly uses experimental and traditional methods. Recent years, with the rapid development of computer technology, the simulation technology make great progress, too. The text mainly analyzes FEM，DEM，meshless method,CFD,ANN and the coupling of the two methods in which these five numerical simulation method, which applied to dynamic simulation of the soil-ditching component interaction; In addition, different method existing difference are analyzed, lastly, future emphasis and trend of development are proposed.

Key words：simulation technology; soil-ditching component; difference; trend of the development

文章编号：1003-188X(2016)03-0236-05

中图分类号：S222.5+2

文献标识码：A

作者简介：曹中华(1988-)，男，山西长治人，硕士研究生，(E-mail) 1534254366@qq.com。通讯作者：陈建(1957-)，男，重庆人，教授，博士，博士生导师，(E-mail)jianchen@swu.edu.cn。

基金项目：重庆市应用开发计划重点项目(cstc2013yykfB70002)

收稿日期：2015-04-13