基于SPH/FEM耦合算法的土壤切削仿真与研究

刘亚超，王秀山，冯 敏，张合虎

(河南农业大学 a.机电工程学院；b.园艺学院，郑州 450000)



基于SPH/FEM耦合算法的土壤切削仿真与研究

刘亚超a，王秀山a，冯 敏b，张合虎a

(河南农业大学 a.机电工程学院；b.园艺学院，郑州 450000)

ANSYS/LS-DYNA仿真软件，采用SPH/FEM耦合算法模拟旋转刀具切削土壤作业过程。在SolidWorks中建立刀具—土壤模型，用ANSYS/LS-DYNA进行前处理，修改K文件建立耦合模型，在LS-PREPOST中查看LS-DYNA971求解结果，对数据进行二次处理，分析作业过程能耗变化及应力分布，为后期优化刀具提供依据。文中土壤本构模型采用MAT_147(MAT_FHWA_SOIL)材料，对材料屈服准则进行详细阐述，相比标准Mohr-Coulomb屈服准则,修改后的Mohr-Coulomb屈服准则引入了修正系数Ahyp，分析了Ahyp取值对屈服面的影响，得出Ahyp与土壤内聚力和内摩擦角关系；阐述SPH算法和SPH/FEM耦合算法,分别用FEM算法、SPH算法和SPH/FEM耦合算法模拟土壤切削过程。结果表明：在土壤切削前期，网格没有发生畸变时，3种方法模拟结果相近，随着有限元网格发生畸变，FEM算法产生了误差，验证了SPH/FEM耦合算法在土壤切削仿真过程中的可行性与准确性。

ANSYS/LS-DYNA；SPH/FEM耦合算法；旋转刀具；土壤切削

0 引言

目前，我国农业生产的除草方式主要为化学除草，长期使用化学除草危害性非常大，导致杂草抗药性增加、土壤生态环境恶化和各种灾害频发等，这与我国人民日益追求的绿色健康生活形成了突出的矛盾。因此，研究除草机器人代替传统的人工物理除草是今后的趋势，国内外的研究均围绕此方面展开[1-7]。本课题组目前正在研究仿生除草机器人，如图1所示。仿生除草机通过视频采集系统采集根部及株间杂草分布图像并传输至图像数据处理系统，经处理后，将执行指令反馈给除草系统；除草系统根据指令驱动滑块带动刀具切入土壤，将杂草根部切断，并对土壤进行松土。本文研究刀具切削土壤过程中土壤应力变化、挤压变形、失效破坏及刀具应力分布对预测土壤结构变化、减少切削阻力、节约能耗和优化刀具设计具有重要意义。

随着计算机技术的不断发展，基于ANSYS/LS-DYNA的动力学仿真技术也越来越成熟。国内外利用FEM有限元法在模拟土壤切削方面已取得较大进展[8-13]。夏俊芳等[14]利用ANSYS/LS-DYNA对螺旋刀辊进行动力学仿真，所得结果与实际理论值相符。刘晓红等[15]对振动深松铲切削土壤进行有限元模拟分析，得出了切削过程功耗变化和等效应力情况。朱留宪[16]利用SPH算法对微耕机旋耕刀进行有限元仿真与优化。丁峻宏等[17]对自适应网格划分、ALE及SPH等3种方法进行仿真模拟对比，结果表明：3种方法在解决土壤切削时各有特点。

图1 除草机实验平台

近年来，土壤切削仿真模拟大多采用FEM有限元法，或SPH光滑粒子法，此两种方法都存在一些缺陷：FEM有限元法在模拟大变形时容易发生网格畸变(见图2)，需对大变形区域进行网格细化，计算时间会大大增加，效率随之降低；SPH光滑质子法解决了FEM有限元法网格畸变的问题，但相比FEM算法计算效率大打折扣，且SPH算法边界处理不如FEM算法方便。因此，为了发挥两者优势，研究人员[18-21]开始对SPH/FEM耦合算法进行研究，并进行实例仿真模拟，仿真结果与实验结果相符。

图2 有限元网格发生畸变

目前，关于SPH/FEM耦合算法模拟土壤切削仿真的参考资料很少，本文基于ANSYS/LS-DYNA采用SPH/FEM耦合算法模拟旋转刀具切削土壤作业过程，利用LS-PREPOST后处理软件读取LS-DYNA971计算结果，进行数据研究和二次处理[22-23]；对刀具切削土壤作业过程进行显示动力学仿真，分析作业过程中应力分布、功耗变化，为后期优化除草机提供参考和依据。

1 算法简介

1.1 SPH方法基本理论

传统的有限元法和有限差分法虽然已经得到广泛应用，但在某些方面具有局限性，尤其是在处理运动不连续时需要每一步重划网格以保证网格线在整个求解过程中与不连续相符。但是，这样会引起很大困难，如在求解过程中需要网格映射，会导致精度降低和计算程序复杂化，同时相应费用将会增加[24]。因此，网格划分限制了它在高速冲击侵彻、爆炸效应、磨蚀及材料等方面的应用。

为了克服和消除以上限制，光滑质点动力学法(SPH)应运而生。光滑质点动力学法的离散化不使用单元，而是使用固定质量的可动点，质量固定在质点的坐标系，所以SPH也算是拉格朗日型，其基本方程也是守恒方程和固体材料本构方程，标准单元和SPH节点的拉格朗日代码非常相似[25]。由于SPH法不使用网格，没有网格畸变问题，可以处理大变形问题，其可以简单而精确地实现复杂的本构行为，也适用于材料在高速加载速率下断裂这个困难问题。

SPH的基础是插值理论[24]，各质点的相互作用借助于插值函数来描述，利用插值函数给出量场在一点的核心估计值(Kemel Estimate),将连续介质动力学的守恒定律由微分方程形式转换为积分形式，进而转换为求和。

SPH方法中质点近似函数为

∏h∫(x)=∫f(y)W(x-y,h)dy

(1)

其中，W是核函数(插值核)；h是光滑长度，光滑长度随时间和空间变化。

核函数W使用辅助函数θ进行定义，有

(2)

其中，d是空间维数；W(x,h)是尖峰函数，如图3所示。

图3 核函数

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

(3)

其中，C是常量，由空间维数决定。

光滑长度h对计算效率和精度有重要影响，目的是为了保持质点邻域内有足够的质点，以确保质点连续变量近似有效，如果采用固定的光滑长度材料膨胀时会导致数值畸变[24]。当质点分离时，光滑长度增加；当质点汇聚时，光滑长度减少，其变化范围在一定范围，可表示为

C1·h0

(4)

其中，h0是初始光滑长度；C1、C2是缩放因子；如果C1=C2，h则为固定光滑长度，不随时间和空间变化。

另外，SPH算法中的边界处理不如有限元网格成熟，有限元网格模型可以通过对边界节点进行约束来处理边界，但在SPH方法中不适合。因此，对于SPH方法，边界需要进行特别处理，通过添加关键字在边界创建虚拟质点来定义边界对称面，如图4所示。

对靠近边界处的质点，通过自身映射自动创建虚拟质点，虚拟质点具有质点相同的特性，因此可以对质点产生近似对称作用[24]。

图4 SPH边界处理

1.2 SPH/FEM算法

图5为SPH/FEM的耦合示意图。其中，左边部分是SPH质点，右边部分为FEM网格。两个部分通过罚函数约束来实现力学参数传递；网格边界区域和SPH质点接触，通过接触类型中的“点-面”固连接触进行耦合定义[24]。

图5 SPH/FEM的耦合示意图

在耦合算法中，耦合界面附近的粒子密度相对单元密度越大，算法的计算精度越精确，通常粒子密度与单元密度相当或稍大时，即可获得比较理想的计算结果；但是，耦合界面处的区域内总存在一定的应力振荡。实验表明：随着光滑长度的增大，粒子区域内的应力振荡幅度减小[19]。为了保证耦合计算精确和稳定，文中最小光滑长度取0.2倍粒子间距，最大光滑长度取2.0倍粒子间距。

固连接触失效准则[26]为

(5)

其中，fn为法向力；fs为剪切力；fn.fail为法向失效力；fs.fail为剪切失效力；m1为法向力指数；m2为剪切力指数。

1.3 FEM、SPH、SPH/FEM耦合3种算法仿真对比

3种算法模拟土壤切削时刀具为首次切削土壤，此时最大的等效应力分别为3.047E-03、3.334E-03、3.339E-03。结果相差不大，证明耦合算法的正确性，如图6～图8所示。在前期研究时发现，土壤切削前期，网格没有发生畸变时，3种方法模拟结果相近，随着有限元网格发生畸变造成FEM算法产生了误差。

图6 FEM算法仿真

图7 SPH算法仿真

图8 SPH/FEM耦合算法

2 SPH/FEM耦合模型

2.1 作业原理

除草系统中的滑块左右张合运动带动左右刀具相向运动，刀片在电机的带动下，以平行于地面的状态高速切入土壤中，将杂草的根部切断，同时对土壤进行松土。图9为刀具工作示意图；图10为刀具在ADAMS仿真中的工作路径。

图9 刀具工作示意图

图10 刀具株间除草运动轨迹

2.2 模型建立

为了简化模型、提高效率，本文只研究一侧刀具在土壤中的工作状况，为了减少刀具入土前的无用时间，建立模型时刀具应紧靠土壤模型但不接触。

通常ANSYS/LS-DYNA建立复杂模型比较困难，本文先在SolidWorks进行建立0.06×0.06×0.08的刀具-土壤模型，并转换为Parasolid文件；通过与ANSYS/LS-DYNA无缝连接，将文件导入ANSYS/LS-DYNA进行前处理。定义模型为SOLID164单元，并划分网格，保存K文件，然后用UltraEdit对K文件进行修改建立耦合模型。图11为SolidWorks模型，图12为有限元网格模型，图13为耦合模型。

图11 SolidWorks模型

图12 有限元网格模型

图13 耦合模型

2.3 模型参数设定

2.3.1 刀具模型

刀具在土壤中运动时几乎不会发生形变，定义刀具模型材料为刚性体MAT_RIGID，采用自由网格划分法划分网格。刀具材料参数如下：

密度/kg·m-3：7 890

弹性模量/Pa：2.1E+11

泊松比：0.3

2.3.2 土壤模型

2.3.2.1 土壤材料

土壤模型参数的选取直接决定了仿真结果的准确性，为了能真实地反映出刀具切削土壤实际情况，选取ANSYS/LS-DYNA中土壤专用材料MAT_147(MAT_FHWA_SOIL)[27-28]。这种材料考虑了土壤的密度、土粒相对密度、水密度、粘塑性、体积弹性模量、剪切模量、内摩擦角、土壤含水率及孔隙水效应等，可以很好地表现出土壤在失效后的特征，材料达到失效后可以及时删除失效单元。

土壤的失效准则采用修正的Mohr-Coulomb(莫尔-库伦原则)屈服准则[28]，相比标准Mohr-Coulomb屈服准则，修正的Mohr-Coulomb屈服准则引入了Drucker-Prager修正系数Ahyp，修正后的方程为

(6)

其中，σy为屈服面力；P为压力；φ为内摩擦角；J2为第二不变应力偏量；K(θ)偏平面角函数；Ahyp为修正系数；C为内聚力。

当Ahyp为0时，修正的Mohr-Coulomb屈服准则及为标准屈服准则,如图14为标准Mohr-Coulomb屈服准则与修正的Mohr-Coulomb屈服准则线性表示。

图14 标准Mohr-Coulomb屈服准则与

由图14可知：修正的 Mohr-Coulomb屈服准则在屈服面力为0时，“脊”部更为平滑，因此更接近真实失效条件。

选择适当的Ahyp系数对仿真的稳定性很重要，研究人员[28]发现，Ahyp系数选择与内摩擦角和内聚力可以用一个合理的关系表示，有

(7)

随着Ahyp增大，σy随之增大。如果Ahyp较大，则明显偏离Mohr-Coulomb屈服准则的包络线。为了保持Mohr-Coulomb屈服准则的包络线，取Ahyp值为1.0E+2Pa。图15为Ahyp影响下的屈服面。

图15 Ahyp影响下的屈服面

2.3.2.2 土壤参数

结合本地土壤情况并实地测量土壤参数(见图16)，参考Soil Material Model 147[27-28]及《农业机械设计手册》[29]，土壤的主要参数如表1所示，采用kg-m-s单位制。

图16 土壤含水率测试

名称单位数值土壤密度kg/m32350土粒相对密度kg/m32．65体积模量Pa2．8E+07剪切模量Pa1．8E+07内摩擦角rad1．1内聚力Pa6．2E+03含水量0．064AhypPa1．0E+2

2.4 修改K文件

1)删除大变形区域的拉格朗日有限元网格，相应区域添加SPH粒子及相关Part信息。

2)添加SPH相关关键字。添加*CONTROL_SPH用于设定SPH质点计算，添加关键字*SECTION_SPH用来定义SPH质点属性，质点光滑长度CSLH取值范围在1.05～1.3；小于1是不允许的，大于1.3将会增加计算时间[24]，本文中CSLH取值1.2。

3)定义相关接触。添加*CONTACT_ERODING_NODES_TO_SURFACE关键字用于定义刀具与SPH质点之间的接触为侵蚀接触算法;添加*CONTACT_TIED_NODES_TO_SU-RFACE_OFFSET关键字用于定义SPH质点与有限元网格的固连接触。

4)处理边界约束。考虑土壤与大地接触的实际情况，添加*BOUNDARY_NON_REFLECTING关键字,定义模型中土壤有限网格前后面及底部为无反射界面。SPH无法添加无反射边界，需添加*BOUNDARY_SPH_SYMMETRY_PLANE来定义其边界。

5)定义速度载荷曲线。刀具移动速度和旋转速度在模拟仿真中是恒定不变的，添加载荷曲线*DEFINE_CURV*BOUNDARY_PRESCRIBED_MOTION_RIGID关键字定义刀具旋转角速度。最后保存K文件，提交到ANSYS/LS-DYNA971进行求解。

3 仿真结果分析

3.1 土壤等效应力分析

当刀具开始切削土壤时，土壤先发生形变，随着刀具切入土壤部分的增加，变形量也随之增加，直到达到屈服准则，土壤开始被破坏，其应力是先增加后减少最后趋于平稳。图17为刀具作业过程中土壤应力变化情况。

(a) 运行时间0.00225s (b) 运行时间0.00325s

(c) 运行时间0.00450s (d) 运行时间0.00650s

图17(a)为刀具刚刚切入土壤，图17(b)为刀具部分切入土壤，图17(c)为刀具全部切入土壤，图17(d)为刀具从土壤中切出。从图17中可以看出：刀具刚开始切入土壤时等效应力较大，随着刀具切入土壤部分增多，等效应力趋于平稳。整个切削过程中，土壤受力没有大的波动，切削过程较为平稳。由于受到土壤压力及切削力，图中最大应力均位于刀具内侧曲面处，因此在加工刀具时应注意进行局部强化处理。

图17(d)中，SPH整块粒子向上翘起，说明刀具切削土壤时对土壤起到了松土的效果，达到设计要求。土壤切削过程中，应力向运动方向扩散，有利于未切削区域土壤松散，减少切削阻力。

3.2 切削能耗分析

刀具切削土壤时，刀具的动能保持不变，切削能耗主要包括为克服土壤变形和土壤运动时需要的能耗。随着刀具切入土壤，土壤能耗逐渐增加，最后趋于平稳，如图18所示。

将能耗曲线对时间进行求导，得出功率变化曲线，如图19所示。

图18 土壤能耗变化曲线

图19 土壤功率变化曲线

当刀具切入土壤时，功率迅速增加，随着刀具切出土壤，功率随之迅速降低。刀具全部切入时，功率达到最大40W。由此可知，每个高速扭矩电机的输出功率至少为40W。

4 结论

1)采用SPH/FEM耦合算法模拟刀具在绕自身中心旋转的情况下水平运动切削土壤的作业过程，在大变形区域采用SPH算法，其余区域采用有限元网格算法。结果表明：耦合算法相比有限网格法可以很好地解决网格畸变的问题，相比SPH算法减少了计算时间，提高了工作效率。

2)FEM在处理大变形时需要格外注意网格畸变问题，SPH算法在处理大变形时不用考虑网格畸变问题，但计算效率不及FEM,在处理边界问题上也不如其成熟，无法添加无反射边界，需要用“BOUNDARY_SPH_SYMMETRY_PLAN”来定义边界。

3)比较分析了FEM、SPH、SPH/FEM耦合3种算法模拟土壤切削过程，在土壤切削前期，网格没有发生畸变时，3种方法模拟结果相近；随着有限元网格发生畸变造成FEM算法产生了误差，从而验证了SPH/FEM耦合算法在土壤切削仿真过程中的可行性与准确性。

4)仿真得出刀具在切削过程中的土壤受力及应力分布情况，最大应力均位于刀具内侧曲面处，为加工刀具及局部强化处理提供了依据。土壤切削过程中，应力波向运动方向扩散，有利于未切削土壤松散破碎，节约切削能耗。整个切削过程中，土壤受力没有大的波动，切削过程比较平稳。

本实验基于ANASYS/LS-DYNA动力学仿真软件，采用SPH/FEM耦合算法模拟旋转刀具切削土壤，目前处于研究阶段，随后将进行一系列系统仿真并与实验结果进行对比，并对所建立模型进行修改完善，也为今后同类问题研究提供参考。

[1] 李江国，刘占良，张晋国，等.国内外田间机械除草技术研究现状[J].农机化研究,2006(10)：14-16.

[2] J K Kouwebhoven.Intra-raw Mechanica Weed Control-possibilities and Problems [J]. Soil and Tillage Research,1997,41:87-104 .

[3] 张朋举，张纹，陈树人，等.八爪式株间机械除草装置虚拟设计与运动仿真木[J].农业机械学报,2010,41(4)：56-59.

[4] Z Gobor，P Schulze Lammers，M Martinov. Development of a mechatronic intra-row weeding system with rotational hoeing tools: Theoretical approach and simulation[J].Computers and Electronics in Agriculture，2013，98：166-174.

[5] 陈勇,田磊,郑加强.基于直接施药方法的除草机器人[J].农业机械学报,2005,36(10):91-93.

[6] 梁远，汪春，张伟，等. 3ZCS-7型复式中耕除草机的设计[J].农机化研究,2010,32(6)：21-24.

[7] 韩豹，申建英，李悦梅.3ZCF-7700型多功能中耕除草机设计与试验[J].农业工程学报,2011,27(1)：124-129.

[8] 齐龙.基于ANSYS/LS-DYNA的松土刀切削有限元仿真[J].农机化研究，2015，37(7)：48-52.

[9] 马爱丽，廖庆喜，田波平，等.基于ANSYS/LS-DYNA的螺旋刀具土壤切削的数值模拟[J].华中农业大学学报，2009，28(2)：248-252.

[10] 周明，张国忠，姚兴林.有限元法在土壤切削过程模拟中的应用[J].农机化研究，2009，31(9)：190-192.

[11] 周明，张国忠，许绮川，等.土壤直角切削的有限元仿真[J]. 华中农业大学学报，2009，28(4)：491-494.

[12] 刘增汉，杨坚，周仕成，等.基于LS-DYNA的土壤—甘蔗—切割器系统的动力学仿真研究[J].农业装备与车辆工程，2010(5)：29-31.

[13] 夏哲浩，姚立红，阚江明.基于ANSYS/LS-DYNA旋转刀具切削土壤与木材的数值模拟[J].森林工程，2016,32(1)：43-47.

[14] 夏俊芳，何小伟，余水生，等.基于ANSYS/LS-DYNA的螺旋刀辊土壤切削有限元模拟[J].农业工程学报,2013(10)：34-39.

[15] 刘晓红，邱立春.振动深松铲土壤切削有限元模拟分析-基于ANSYS/LS-DYNA[J].农机化研究，2017，39(1)：19-23.

[16] 朱留宪.基于SPH算法的微耕机旋耕刀有限元仿真与优化[D].重庆：西南大学，2012.

[17] 丁峻宏，金先龙，郭毅之，等. 土壤切削大变形的三维数值仿真[J]. 农业机械学报， 2007， 38(4)：118-121.

[18] 张志春，强洪夫，高巍然.光滑粒子流体动力学—有限元接触算法研究[J].高压物理学报，2011，25(2)：97-102.

[19] 肖毅华.有限元法与光滑粒子法的耦合算法研究[D].长沙：湖南大学，2012.

[20] 吴文峰.基于SPH-FEM耦合方法的微细铣削仿真与实验研究[D].南京：南京航空航天大学，2013.

[21] 刘春含.有限元和光滑粒子耦合方法在爆炸问题中的应用研究[D].长沙：湖南大学，2013.

[22] 尚晓江，苏建宇.ANSYS/LS-DYNA 动力分析方法与工程实例[M].北京：中国水利水电出版社，2005：102-116.

[23] 白金泽.LS-DYNA3D理论基础与实例分析[M].北京：科学出版社，2005：204-231.

[24] 李裕春，时党用，赵远. ANSYS+11.0LS-DYNA基础理论与工程实践[M].北京：中国水利水电出版社，2005：434-431.

[25] 宿崇，唐亮，侯俊铭，等.基于FEM与SPH耦合算法的金属切削仿真研究[J].系统仿真学报，2009,21(16)：5002-5005.

[26] 段念，王文珊，于怡青，等.基于FEM与SPH耦合算法的单颗磨粒切削玻璃动态过程仿真[J].中国机械过程，2013,24(20)：2716-2721.

[27] Lewis B A. Manual for LS-DYNA soil material model 147(FHWA-HRT-04-095)[R]. Department of Transportation:Federal,Highway Administration,U.S.A.,2004.

[28] Reid J D,Coon B A,Lewis B A,etc. Evaluation of LS-DYNA soil material model 147(FHWA-HRT-04-095)[R].Department of Transportation: Federal Highway Administration,U.S.A.,2004.

[29] 中国农业机械化科学研究院.农业机械设计手册(上)[K].北京：中国农业科学技术出版社，2007：237-238.Abstract ID:1003-188X(2017)07-0021-EA

The Soil Cutting Dynamics Simulation and Research Based on SPH/FEM Coupling Algorithm

Liu Yachaoa, Wang Xiushana, Fen Minb, Zhang Hehua

(a.Mechanical and Electrical Engineering College; b.College of Horticulture, Hennan Agricultural University, Zhengzhou 450000,China)

In this paper, based on ANSYS/ LS-DYNA ,The SPH/FEM coupling algorithm are adopted to simulate the rotary cutter cutting soil process,tools-soil model is established in the SolidWorks, Before processing in the ANSYS/LS-DYNA and the mod-el is set up by mod-ifeied the K file,check the LS-DYNA971 results and the datais secondary processing by th-e LS-PREPOST, analysis of process energy consumption ch-ange and stress distribution, provide the basis for later optimization tool;The soil constitutive model adopts MAT_147 (MAT_FHWA_SOIL) materials, the material yieldc-riterion are expounded in details,co-mpare with the standard mohr-coulomb, Ahyp correction coefficient is introduced of The modified mohr-coulomb，analyzed the Ahyp values influence on the yield surface, Ahyp relationship with cohesion and internal friction angle are presented also;Explain SPH algorithm and SPH/FEM coupling algori-thm, using FEM algorithm respectively, SPH algorithm and SPH/FEM coupling algorithm to simulate the soil cutting process,the results show that in the early stage of the soil cutting, gridwithout distortion occurs, three ways of the simulation results, the finite element mesh distort on caused by FEM algorithm error, thus in SPH/FEM coupling algorithm is verified the feasibility and accuracy in the process of soil cutting simulation.

ANSYS/LS-DYNA；SPH/FEM coupling algorithm；rotational cutting tool; soil cutting

2016-08-30

国家自然科学基金项目(U1204524)；河南省高等学校骨干教师资助项目(2014GGJS-034)

刘亚超(1990-),男，河南安阳人，硕士研究生，( E-mail)liuyachao0808@163.com。

王秀山(1975-)，男，郑州人，副教授，博士，( E-mail)towxs@163.com。

S224.1

A

1003-188X(2017)07-0021-07