基于LS-DYNA三七挖掘铲仿真分析及试验

解开婷，张兆国*，王一驰，唐金鑫，余小兰

（1.昆明理工大学农业与食品学院，昆明 650500；2.云南省高校中药材机械化工程研究中心，昆明 650500）

三七是我国名贵中药材，具有较高药用价值[1]，因其根系复杂，种植深度较深[2]，土壤质地黏稠，收获三七挖掘阻力大，现有根茎类挖掘铲难以满足挖掘需求[3]。为实现挖掘铲减阻，需深入研究三七挖掘铲与土壤作用情况，为后期优化三七挖掘铲提供基础。

近年来，有限元法逐渐应用于农业机械仿真试验，蒋建东等应用有限元法仿真分析土壤旋切振动减阻，得出较优旋切碎土方式[4]。周文琪等应用有限元仿真分析液肥穴施肥机扎穴针体与土壤相互作用，建立可靠ANSYS针体模型[5]。Mouazen等应用有限元仿真分析土壤塑性，设计得到地域性较强耕作工具[6]。侯俊铭等利用有限元法建立花生壳和花生仁模型，得出不同加载方式下花生变形量，为改进花生脱壳装备提供理论基础[7]。以上学者研究表明，有限元分析方法可获取物体受力数值，观察不同物体间相互作用与形变，适用于难以直接求取物体间相互作用的试验研究。但目前尚未见利用有限元仿真方法，探究三七挖掘铲与三七种植土壤的相互作用。鉴于此，本文以三七挖掘铲为研究载体，采用CFD-POST与LS-DAYN对挖掘铲作业过程作动态数值模拟，分析挖掘铲运动特性与土壤受扰动情况，通过台架试验验证仿真模型准确性，为后期优化三七挖掘铲奠定理论基础。

1 仿真试验材料参数设定与模型构建

1.1 土壤基本物理参数试验与设定

三七种植土壤为黏重红壤，以云南弥勒东山（见图1）采集的三七种植土壤为研究对象，通过标准土壤筛分试验得出70%土壤颗粒粒径小于1 mm，利用激光粒度仪测定粒径小于1 mm土壤样本，得出96%土壤粒径介于2.38×10-3~2.13×10-1mm，土壤平均粒径为0.01 mm。参照《NY/T1121.23-2010土粒密度测定法》《DIN ISO 11465-1996土质·依据重量测定干物质和含水量·称重法》分别测得土壤密度、试验土槽土壤含水率如表1所示。查阅相关文献获取含水率为（20±1）%时土壤弹性模量为2.3×104MPa，泊松比为0.4[8]。

土壤剪切模量可由弹性模量与泊松比近似确定[9]，如式（1）所示。

式中，G为土壤剪切模量（Pa）；E为弹性模量（Pa）；λ为泊松比；∂为内摩擦角（°）。

表1 土壤材料参数Table 1 Material physical parametersof soil

采用快剪试验测定21.3%含水率土壤分别在5×10-2、7.5×10-2、10×10-2MPa垂直压力下的土壤抗剪切强度，根据试验数据与式（2）、（3），得出土壤抗剪切强度与垂直压力间数学关系式（4）。

式中，Gp为土壤抗剪切强度（MPa）；c为土壤黏聚力（MPa）；p为外施垂直压力（MPa）。由计算得，土壤内摩擦角为18.87°。

1.2 Mohr-Coulomb屈服准则

Mohr-Coulomb屈服准则是Mothr准则的特殊形式，即材料在复杂应力状态下发生破坏，材料在外荷载作用下，沿某一斜截面上的剪应力达到某一极限值造成[10-11]。Mohr-Coulomb屈服准则能够体现材料的塑性变形特性，且参数较少，容易测量，广泛应用于土壤、岩土一类材料的仿真分析。

真实准确的土壤材料参数设定是保证仿真结果准确前提条件，为准确描述土壤在挖掘铲作用下“掘裂”特性，调用ANSYS/LS-DYNA中，遵循Mohr-Coulomb屈服准则关键字材料MAT_147界定土壤参数。

式中，Fm为模型屈服表面力（MPa）；J2为第2个应力偏张不变量（MPa）；k（θ）为应力罗德角函数；A为修正相似因数，2 360。

因ANSYS/LS-DYNA软件中无材料MAT_147，需利用后处理软件LS-PREPOST修改ANSYS/LSDYNA中生成的K文件，修改MAT_147材料主要参数为：土壤密度、剪切模量、内摩擦角、含水率和黏聚力，其余参数参照MAT_FHWA_SOIL中默认值设定。

1.3 挖掘铲结构参数设计

三七种植土壤较为黏重，要求挖掘铲有足够强度，采用凸面挖掘铲结构简单，强度高，且便于后续简化仿真模型，提高仿真运算效率。挖掘铲铲片数量n由三七种植垄宽与挖掘铲单铲铲宽共同决定。

式中，F为整块挖掘铲所受土壤阻力（N）；f为单片挖掘铲所受土壤阻力（N）。

分析单片凸面铲受力情况，凸面铲侧视图与俯视图如图2所示。

图2 凸面铲Fig.2 Convex digger

当挖掘铲入土作业时，凸面铲单片铲尖受到来自三七与土壤阻力f，将f沿铲尖分解出两个力f1和f2，得到：

式中，μ为摩擦因素，一般取0.58~0.71；

根据摩擦公式：

式中，β0为铲尖夹角（°）；φ为对凸面铲最大摩擦角（°）。

计算得β0≤120°满足条件，所以取铲尖夹角为90° 。

根据图3土壤在挖掘铲上受力分析可得：

图3 受力分析图Fig.3 Stressanalysis diagram

解式（9）~（11）得：

式中，F1为铲面受到的水平力（N）；F2为铲面摩擦力（N）；G为土壤受到的重力（N）；FN为支反力（N）；γ为铲面倾角（°）。

为满足条件，凸面三七挖掘铲入土角需满足10°≤γ2<γ1≤20°。查阅相关文献，三七挖掘铲单片铲宽选取60 mm、挖掘铲垂直总高度取68 mm可满足挖掘条件[12]。

为保证挖掘过程中将三七挖掘出土壤并提升一定高度，根据公式：

式中，b为铲长（mm）；b1为一阶凸面铲长度（mm）；b2为二阶凸面铲长度（mm）；h为挖掘铲垂直总高度（mm）；h1为一阶凸面铲高度（mm）；h2为二阶凸面铲高度（mm）；γ1为一阶凸面铲与水平面夹角（°）；γ2为二阶凸面铲与水平面夹角（°）。

计算得挖掘铲总长度为302.3 mm，近似取值302 mm。

表2 挖掘铲材料参数Table 2 Material physical parameters of excavation

1.4 挖掘铲-土壤模型

本文利用三维制图软件CATIA建立挖掘铲及土壤模型，保存为igs格式，设置其初始位置如图4所示。仿真整块挖掘铲，对电脑配置要求较高，仿真时间过长，因此仅模拟三片铲片入土作业，土壤模型尺寸设置为400 mm×400 mm×800 mm。将模型导入ANSYS中，采用Mapped方式划分几何体，通过Sizing控制网格划分质量，将土壤与挖掘铲模型划分为26 375个实体单元，得出划分后网格如图5所示。

图4 仿真三维模型Fig.4 3Dsimulation model

设置挖掘铲与土壤间接触为面对面接触中Eroding接触，确保在模型外元件失效后，其余元件仍可接触。设置挖掘铲沿Z方向自由度为0，定义挖掘铲为主动面，设置X向平动约束，土壤为从动面设置固定约束。团队前期深入研究挖掘铲运动参数[13]，得出挖掘铲最优运动参数为：前进速度0.6 m·s-1，挖掘铲与土壤夹角20°，挖掘深度为30 cm，因此定义最优运动参数为仿真参数，方向为X正方向，求解时间为2 s，结果输出步数为400步。

图5 有限元网格划分结果Fig.5 Result of finite element meshing

图6 挖掘铲运动方向Fig.6 Dig the direction of shovel movement

2 挖掘铲-土壤互作分析

2.1 挖掘铲入土过程分析

以挖掘铲前进速度0.6 m·s-1，入土角20°，挖掘深度30 cm为参数开展挖掘铲与土壤互作仿真分析。在作业过程中挖掘铲承受3个方向阻力，分别是挖掘铲前进方向（X轴）所受的阻力，即挖掘铲土壤阻力，Y轴方向土壤对挖掘铲的压力，即挖掘铲铲起物重力，Z轴方向土壤对挖掘铲的摩擦力。

图7表示不同时刻挖掘铲不同部位所受阻力情况。t=0 s时，挖掘铲铲尖刚接触土壤，此时在X轴方向，挖掘铲铲尖所受土壤阻力为750 N；t=0~0.4 s时，挖掘铲沿X方向移动，此时铲尖未全部入土，铲尖所受阻力相对集中，挖掘铲横梁所受阻力较大，中间与横梁衔接部分铲片所受阻力最大，同时与铲尖接触的土壤单元失效，继而破碎；t=0.4~0.8 s时，铲尖已全部入土，此时挖掘铲整体所受阻力比较均匀；t=0.8~1.2 s时，随挖掘铲移动距离增加，土壤单元不同程度破碎，t=1.2~1.6 s时，此时土体不再均匀，挖掘铲所受阻力集中点发生移动；t=2.0 s时，挖掘铲完成一次挖掘作业。

图7 各时刻挖掘铲阻力变化情况Fig.7 Excavate the changeof shovel stress at each moment

图8表示不同时刻土壤破碎情况。t=0 s时，土壤产生以挖掘铲铲尖为中心，向四周阶梯扩散的楔形破碎痕迹，此时土壤主要发生X方向形变，Y轴与Z轴方向土壤未发生明显形变；t=0~0.4 s时，随挖掘铲位移变化，土壤受扰动范围变大，加剧土壤形变程度，铲尖接触到的土壤单元变形最大，受到的应力也最大。被挖掘铲抬起的土壤向铲面两侧滑落，因此Z方向土壤形变范围也随之变大。此时铲尖未完全进入土体，Y轴方向仍未产生明显形变；t=0.4~0.8 s时，此时铲尖已完全入土，土壤最大形变位置随挖掘铲移动而向周围扩散，相邻土壤单元也随之发生形变。由于挖掘铲前进过程中受到X方向阻力，部分土壤向挖掘铲前进方向的反方向移动，从铲片间隙间落下并堆积，堆积土壤对挖掘铲后方土壤产生作用力，使得后方土壤发生Y方向形变；t=0.8~1.2 s时，土壤受力范围继续扩大，承受阻力值也随之增加，Y方向土壤形变加剧，表明挖掘铲作业过程中，漏土效果显著；t=1.2~1.8 s时，土壤在X、Y、Z轴方向形变均随挖掘铲运动继续扩张；t=2.0 s时，挖掘铲完成一次仿真作业，未与挖掘铲接触部分土壤未产生明显形变。

图8 各时刻土壤破碎情况Fig.8 Soil fragmentation at each moment

2.2 挖掘铲所受三向阻力分析

为节省仿真运算时间，仅模拟3片铲片挖掘过程。如图9所示，在挖掘铲初始速度为0.6 m·s-1，挖掘铲挖掘深度为20 cm，入土角为20°情况下，获得挖掘铲在一个铲土过程中三向阻力变化关系。

图10表示土壤在挖掘铲作用下，土壤移动轨迹云图。t=0~0.25 s时，由于土体模型在挖掘铲未运动时是完整几何体，当挖掘铲铲尖刚开始运动时，铲尖接触部分由完整土壤单元直接破碎，形变最明显，对应图10中红色部分。因此0~0.25 s内，Z轴方向阻力波动最大。0.25 s以后，挖掘铲所受三向阻力变化规律明显，因此选取0.25 s后数据作分析。t=0.25~2.0 s时，X轴方向阻力随挖掘铲位移增加而呈微波波动上升趋势，由于挖掘铲匀速运动，阻力变化幅度不明显，最大阻力值为831.25 N；随挖掘铲前进，铲面所承受土壤增多，因此Y方向阻力随作业距离增加而增加；Z轴方向应力因土块向挖掘铲四周堆散，入土瞬间波动较大，后趋于0，对挖掘铲作用最小。

图9 挖掘铲三向阻力变化曲线Fig.9 Three-direction resistance curve of excavating shovel

图10 挖掘铲作业过程中土壤轨迹云图Fig.10 Soil trajectory cloud map during excavating and shoveling operation

3 试验验证

3.1 台架试验

使用温室三七收获机前端机架作为挖掘铲载体，连接挖掘铲，该三七收获机采用液压控制，前端机箱中安装压力调节表，通过压力调节可控制收获机行驶速度，具有较好稳定性，因此选择该三七收获机有利于提高试验准确性，搭建土槽试验台并开展挖掘铲性能试验。试验前测得土槽土壤含水率为8.1%，通过洒水与翻整，调制土槽土壤水分至21%，与采样地土壤含水率基本相同，以处理部分为土槽试验区开展试验。三七收获机挖掘铲性能试验台如图11所示。测试试验台由温室三七收获机前端机架、挖掘铲连接机架、凸面挖掘铲、数据采集仪、应变片、加速度传感器、振动传感器及笔记本电脑组成。因挖掘铲作业过程中与土壤密切接触，应变片直接贴于挖掘铲上方，导致试验数据偏差较大，甚至损坏应变片，因此试验前，将应变片粘贴于挖掘铲与机架连接板内侧。为保证试验结果准确性，左右两侧均粘贴4片应变片，利用4个数据采集通道，对挖掘铲机架左右两侧受力应变实时监测。利用三轴加速度传感器、振动传感器分别采集挖掘过程中三向加速度与挖掘铲振动情况。通过计算机接收数据采集仪采集的试验数据。

3.2 压力标定试验

以挖掘铲挖掘深度、入土角度为影响因素，挖掘铲所受最大土壤阻力为测量指标，采用半桥测量法获取具体阻力值，连接高精度“S”型拉力传感器与挖掘铲，对其施加不同拉力，得出不同拉力对应的应变片应变。采用数显推拉力计实时显示并记录外施拉力数值，应变片接线方式如图12所示，计算机采集应变片对应应变。为直观获得应变片应变与受力间规律，采用一阶线性拟合法，利用Origin软件求取二者间函数关系，得到应力拟合曲线如图13所示，得出压力标定公式：

式中，P为应变片所受拉力（N）；εr为应变片应变。

图11 挖掘铲性能试验台架Fig.11 Excavation shovel performancetest bench

3.3 挖掘铲性能台架试验与结果分析

3.3.1 三向阻力对比分析

试验得出挖掘铲作业过程中所受三向阻力如图14所示。0~1 s内，三七收获机处于启动准备状态，此时三七挖掘铲未受到三向阻力；1~5 s内，挖掘铲开始作业，铲尖接触土壤并向前运动，此时土槽土壤开始发生形变，因此该时间段内三向阻力变化幅度不稳定，规律不明显，因此选取5~20 s数据作对比分析；5~20s内，三七收获机匀速前进，X轴方向阻力平缓上升，最大阻力值为2433.38N，整体变化幅度小，但实际土槽土壤土质不均匀、机架发动机及风扇工作产生振动等因素，导致机架载体也随之振动，因此变化幅度没有仿真结果中X轴方向阻力变化均匀；Y轴方向阻力随铲面承受土壤的增多而增大，由于铲片间隙存在漏土，部分被挖掘铲抬起的土壤掉落回土槽中，因此Y方向阻力并非绝对线性增长，呈微波抖动上升趋势；Z轴方向阻力主要来源于两侧土壤与铲片间摩擦力，铲片侧面积较小，因此Z轴方向阻力波动较为平稳；t=20 s时，三七收获机暂停前进，因惯性作用，收获机发动机与风扇等部件继续工作至25 s，此时，受机架本身振动影响，还存在三向阻力。X轴、Z轴阻力变化不明显。Y轴方向因未继续前进，铲面土壤部分落下且未掘其新土，其阻力缓慢下降。实际采集的挖掘铲阻力整体波动较大，但对比仿真结果，X轴、Y轴与Z轴阻力变化规律一致。

图12 半桥接线法Fig.12 Half bridgeconnection method

图13 应力拟合曲线Fig.13 Stressfitting curve

图14 挖掘铲实际三向阻力Fig.14 Actual three-way resistance of excavating shovel

因仅模拟三片铲片运动情况，实际挖掘铲共八片铲片，因X轴方向阻力对挖掘铲作用最显著，根据式（6）分析仿真与试验结果误差，得出其相对误差为8.91%。仿真模型较准确，可为后续挖掘铲优化改进提供理论依据。

3.3.2 土壤破碎形态对比分析

如图15所示，实际挖掘铲开始运动时，土壤形成缺口，挖掘作业结束后，土壤自然堆积。

图15 土壤姿态Fig.15 Soil attitude

如图16所示，为直观对比仿真与实际挖掘中土壤形变差异，利用OpenCV灰度化经高反差保留处理的图像，提取挖掘轮廓坐标。如图17所示，利用Origin拟合仿真与实际挖掘轮廓坐标。

图16 轮廓处理流程Fig.16 Contour processing flow

图17 轮廓坐标对比Fig.17 Contour coordinate contrast

根据式（17）~（18）计算得出重合部分仿真土壤挖掘轮廓x、y轴坐标相对误差分别为17.83%、4.50%。

式中，δ%为仿真土壤挖掘轮廓相对误差，下标x、y分别表示x轴、y轴；x1、y1为仿真挖掘土壤轮廓横纵坐标值；x2、y2为实际挖掘土壤横纵坐标值。

由误差分析可知，x轴较y轴误差偏大，y轴仿真土壤形变与实际土壤形变姿态一致，相对误差较小，故本文构建的仿真模型可靠。

4 讨论与结论

a.通过实地调研采样，以云南弥勒三七种植土壤为试验样本，通过土壤基础物理机械试验，得出三七种植土壤基本物理参数。结合三七种植农艺与理论计算，设计挖掘铲主要结构参数。利用CATIA三维制图软件绘制简化的三七收获机挖掘铲与土壤三维模型，选用Mohr-Coulomb屈服准则界定土壤关键字材料MAT_147并修改关键字，得出准确的仿真参数。运用ANSYS/LS-DNYA有限元仿真软件，搭建仿真模型，微观分析挖掘铲在作业过程中动态行为特性以及三向受力情况，为后期三七收获机挖掘铲改进和优化提供方向。

b.通过仿真得出挖掘铲在前进速度为0.6 m·s-1、入土角度为20°、挖掘深度为30 cm时，挖掘铲受到的三向阻力变化情况与土壤形变程度。结果表明，挖掘铲在X轴方向上（即挖掘铲作业方向）所受最大阻力值为831.25 N，对挖掘铲影响最显著；Y轴方向阻力上升趋势最为明显，呈微波抖动上升趋势，所受最大阻力值为565.60 N；Z轴方向受力平稳并趋于0，所受最大阻力值为33.42 N，对挖掘铲影响不显著。

c.搭建三七收获机挖掘铲性能试验台，利用应变片、8通道数据采集仪和多种传感器组成动力学测试系统，得出挖掘铲作业过程中所受实际三向阻力。结果表明，X轴方向所受阻力最大值为2 433.38 N，对挖掘铲影响最显著，Y轴方向阻力呈微波抖动上升趋势，所受最大阻力值为1 923.33 N；Z轴方向受力平稳，所受最大阻力值为360.32 N，对挖掘铲影响不显著。与仿真结果对比得出，仿真试验结果变化符合实际规律，仿真挖掘铲所受阻力与台架试验结果变化趋势一致，相对误差为8.91%；x轴、y轴土壤形变相对误差分别为17.83%、4.50%。误差均低于20%，验证该模型准确性与可行性。