基于SPH算法的纵向开沟过程模拟

王 永，王泽河，马胜涛，弋景刚，袁永伟

(1.河北农业大学机电工程学院，河北保定071001;2.石家庄邮电职业技术学院，河北石家庄050000)

目前，田间试验是农机具设计和研究的主要途径［1］，虽然试验效果直观且可靠性强，但这种研究途径季节性强且费时费力，难以观测到土壤运动轨迹、土壤应力变化等微观变化状况［2－4］。近年来，通过采用数值模拟技术的有限元分析方法模拟农机具作用下土壤运动情况并获得力学数据，不仅可以对农机具进行结构设计和优化，而且可以为零部件设计加工提供可靠依据［5－6］。SPH算法，即光滑粒子流体动力学方法，是通过相互作用质点组来描述连续固体，各个质点上包含质量、速度等物理量，通过跟踪质点运动的轨迹及求解质点的组动力学方程，获得整个系统的力学方法。与自适应网格、任意拉格朗日－欧拉(arbitrary Lagrangian－Eulerian，ALE)等算法相比，SPH算法是一种纯Lagrangian算法，不需要借助网格就可以避免因大变形引起的网格变形、精度破坏等问题［7－9］。

纵向旋转开沟刀辊是新型纵向开沟装置的关键部件，借助ANSYS/LS_DYNA软件，采用SPH算法对刀辊切削土壤开沟过程进行数值模拟。通过ANSYS/LS_PREPOST对土壤应力、抛撒特性和功耗等参数进行分析，通过理论计算和试验，验证SPH算法的可行性，通过研究开沟刀辊切削土壤的作用机理，以期为我国棚室秸秆反应堆新型纵向旋转开沟装置的理论和优化设计提供可靠依据。

1 有限元模型的建立

1.1 开沟刀结构参数及简化

根据秸秆生物反应堆开沟农艺要求设计开沟刀回转半径245 mm，转速540 r/min，入土深度350 mm，刀辊前进速度2 km/h，开沟刀通过螺栓安装在刀辊刀座上，随着刀辊一起转动，其结构如图1所示。其中，A点到B点边界曲线为侧切刃，采用阿基米德螺旋线。

为降低对计算机硬件的配置要求，减少计算量和求解时间，对开沟刀的部分倒角、圆角等结构进行简化，模型见图2。

1.2 刀辊有限元模型

ANSYS/LS_DYNA主要包括直接建模和间接建模2种方式，考虑到刀辊结构及装配复杂性，选择SolidWorks 2015进行刀辊实体建模，并对刀辊模型进行装配，由于开沟较深，开沟刀对双螺旋对称分布在刀辊上且轴向分布1组螺旋板(图3)。其中，刀辊旋转周向切土及抛土，刀辊螺旋轴向升土。

将刀辊模型导出为*.x_t或者*.igs格式，然后导入到HyperMesh中进行网格划分，考虑到刀辊结构及开沟刀曲面的复杂性，定义网格尺寸为6 mm。刀辊模型仿真参数见表1。

表1 刀辊开沟过程动态模拟仿真参数

1.3 土壤材料模型

土壤是包括水、孔隙气和土壤颗粒的三相材料，在很多实例的应用中，土壤被看作连续介质。利用 ANSYS/LS_PREPOST提供的MAT147(*MAT_FHWA_SOIL)材料作为土壤模型，不仅在 Abbo等的基础上［10］修正了莫尔 －库伦(Mohr－Coulomb)准则，而且扩展了应力软化、变形率、运动学硬化和含水率的影响及单元删除，仿真结果更加接近实际土壤特性。修正后的Mohr－Coulomb准则屈服面为双曲线拟合，表述式为

式中:p为压力;φ为内摩擦角;J2为偏应力张量的第二不变量;K(θ)为偏应力面中的极角;Fn为内聚力;a、h1、y、q为定义修正后屈服面和Mohr－coulomb标准屈服面之间贴合度的参数。

结合我国河北地区蔬菜大棚秸秆生物反应堆土壤特性，设计土壤材料模型参数如表2所示。

2 边界条件设置

2.1 SPH 算法简述

在1977年，Lucy等提出SPH算法，运算固定质量的可动点，能够有效解决土壤松散、物性分散和多相的问题［11］。SPH方法可用于求解在任意时刻的一批粒子的能量和速度。要显示粒子的运动情况，首先要构造一个粒子运动的近似函数:

式中:W为插值核函数，定义辅助函数θ，则

式中:v为空间维数，h为光滑长度，W(x，h)是尖峰函数，x、y为设置的坐标值。

SPH中3次B样条是最常用的光滑核，定义为

表2 土壤材料模型参数

式中:C为常数，由空间维数确定;μ为函数θ(μ)的自变量。

将开沟刀辊有限元模型导入到前后处理软件LS_PREPOST中，定义刀辊为刚体。根据表1、表2设置参数，添加SPH土壤模型，SPH粒子均匀分布，约束土壤2个侧面和1个底面，定义刀辊为主接触面，土壤为从接触面。添加刀辊与土壤之间的接触为自动点面接触，在土壤表面单元切削失效的情况下，AUTOMATIC_NODES_TO_SURFACE接触能够根据程序自动定义下一个新接触面，最后设置滑动界面惩罚因子为0.2，动摩擦因子为0.18，静摩擦因子为0.2。在考虑边界条件及刀辊切削方式前提下，建立仿真模型如图4所示。其中，土体为尺寸500 mm×400 mm×500 mm的长方体模型。

LS_PREPOST中设置完成后，生成K文件，最后提交ANSYS/LS_DYNA求解器进行运算。模型求解时间一般比较长，本次仿真设置时长为600 ms，求解时间约为63 h。

3 结果与分析

3.1 切削过程分析

图5是纵向刀辊开沟过程。开沟刀辊初始位于土壤右侧，与土壤无接触(图5－a)，当t=40 ms时，开沟刀辊与土壤接触，土壤受到刀辊上端开沟刀的挤压和剪切作用，被破坏(图5－b);随着刀辊不断前进，刀辊与土壤接触面积增大，在螺旋排布的刀片作用下土壤受到轴向推力，同时在刀辊旋转的作用下，土壤一部分被甩到刀辊一侧，另一部分被推送到前方(图5－c、图5－d);当t=360 ms时，刀辊完全接触土壤，进入稳定状态，刀辊整体旋转完成切土和抛土作业，深层土壤在刀辊螺旋板作用下轴向输送到土壤上层(图5－e、图5－f)，整个开沟过程，由于刀片循环渐进式开沟且螺旋板缠绕在刀辊上，刀辊轴向升土和周向切土、抛土效果好，效率高。

3.2 土壤应力变化

图6为刀辊切削过程中不同时刻土壤等效应力变化云图，随着刀辊不断切削土壤，土壤所受等效应力缓慢增大，涉及的土壤面积不断变大，当t=20 ms时最大等效应力为12.86 MPa;t=40 ms 时最大等效应力为 15.78 MPa;t=120 ms时最大等效应力为30.16 MPa;t=200 ms时最大等效应力为39.74 MPa;t=280 ms时，刀辊完全切入土壤，最大等效应力为43.21 MPa，且趋于稳定。由整个应力云图可以看出，开始时，刀辊未完全进入土壤且循环渐进式切削土壤，土壤等效应变变化很大，最后趋于稳定。整个开沟过程比较平稳，最大等效应力由刀片对土壤的强烈挤压集中在与开沟刀辊接触面上。

3.3 切削功耗分析

开沟过程中，总的能量分为刀辊旋转、土壤运动所需的动能和开沟过程中刀辊与土壤相互作用所损耗的内能。刀辊转速一定，相应动能为定值，整个过程的功耗来源主要是土壤运动。总功耗变化规律如图7所示，在开沟初期，刀辊循环渐进式切土，引起土壤功耗随着时间延长近似呈线性增大趋势，随着刀辊完全切入土壤，土粒相互作用结合力减小，总能量稳定在8.23左右，整个过程符合纵向开沟刀辊设计要求。

3.4 理论验证

功耗计算经验公式:

其中，Kλ=KgK1K2K3K4。

式中:N为切削土壤和抛送土壤所需的功耗，kW;d为开沟刀等效入土深度，cm;vm为机组前进速度;B为耕幅，m;Kλ为旋耕比阻，N/cm2;Kg为含水率16%、开沟深度为350 mm时的土壤开沟比阻，N/cm2;K1为开沟深度修正系数;K2为土壤物质成分修正系数;K3为土壤含水量修正系数;K4为作业方式修正系数。

结合农业机械手册及田间土壤条件，选择设计参数:Kg=15 N/cm2，K1=0.95，K2=0.96，K3=0.9，K4=0.71，d=20 cm，vm=2 km/h，B=0.45 m，由此可得功耗 N=7.91 kW。与图7稳定时总能量进行对比，相对误差为(8.23－7.91)/8.23=0.04≤0.05，在允许范围内，说明采用 SPH 算法进行数值模拟是可行的。

3.5 试验验证

田间试验在河北省保定市蔬菜大棚试验地进行，土层土壤含水率16%，土壤紧实度2.1 MPa。刀辊安装在具有3点悬挂的纵向开沟装置上(图8)。机具前进速度2.0 km/h，切土深度350 mm，刀轴转速540 r/min。通过位于开沟装置前方的田间综合测试系统进行试验，重复4次。试验结果表明，刀辊通过性良好，试验过程中由六分力显示的牵引力如表3所示。由公式N=Fv可得试验过程中平均功耗(N)为8.25 kW，与图7结论一致。说明采用SPH算法进行数值模拟是可行的。

4 结论

采用SPH算法研究纵向旋转开沟刀辊切削土壤开沟过程，建立刀辊切削数值仿真模型，直观显示了刀辊作用下的土壤运动状况。

表3 作业过程中产生的牵引力

对开沟刀辊作用下的土壤开沟仿真模拟显示，整个开沟过程中最大等效应力为43.21 MPa，且切削过程中土壤所受等效应力比较平稳。由纵向刀辊循环渐进式切土，总功耗随时间呈线性增加且增大到8.23 kW时趋于稳定。通过理论和试验验证，误差值在允许范围内，验证了通过采用SPH算法模拟纵向刀辊切削土壤开沟过程的可行性。

SPH算法研究结果表明，纵向旋转开沟刀辊切削过程中，一部分土壤因开沟刀螺旋分布及螺旋板分布而轴向上升，同时，一部分土壤经由刀辊的旋转抛洒到沟的一侧，虽然仍有少量土壤回落到沟渠里，整体开沟效果符合设计要求，为新型纵向旋转开沟装置的设计提供可靠依据。