木薯Manihot esculenta Crantz，也称树薯，呈灌木状，茎秆直立、木质、高2~5 m，茎秆下有块状根。木薯用茎秆繁殖，开沟后将种茎以平放、斜插或直插的方式种植^[1]。近年来，国家各类项目资金支持开展木薯收获机械研究，取得较大进展，但木薯种植机械研究处于起步阶段，进展缓慢，国内还没有技术成熟、可推广的木薯种植机。目前采用切段的茎秆作为种子的种植机械主要有实时切种式和预先切种式。实时切种式机型需要配备专门人工将茎秆喂入切段刀内，切成段后自由落入种沟中，劳动强度大且喂入不连续，易造成伤种、漏种。现有木薯种植机均采用此种形式，难以满足种植要求。预先切种式机型将预先切成段的种茎放入种箱，通过排种装置实现自动排种，可极大降低人工劳动强度，但存在一些待解决的关键技术问题，需进行机理分析和结构创新，且要进一步降低种茎损伤率和漏种率^[1-2]。木薯种茎为一定长度的圆柱状茎秆，表面具有种芽凸起^[2-3]，采用现有的预切种式播种器对种茎进行播种存在充种效果差、合格指数低等技术缺点，作业性能还需提高_。分析种茎运动规律是供种机构、调姿机构、控种机构设计和排种机理研究的重要依据，当前除材料力学特性外，还需对运动和碰撞相关力学特性进行深入分析，但未见木薯种茎相关特性的研究报道。因此，研究木薯种茎在播种过程中的碰撞损伤过程，对设计出低损伤、高效播种装置及推进木薯播种全程机械化发展具有重要意义。

农业物料碰撞损伤受很多因素影响，国内外关于农业物料碰撞损伤的研究主要通过理论及试验方法进行，证明含水率及碰撞速度对农业物料碰撞损伤具有显著影响。然而，受研究方法的限制，此类研究都缺少对农业物料碰撞过程分析，无法揭示碰撞损伤过程中应力、应变及能量等关键参数的变化^[4-5]；特别是针对木薯种茎有限元碰撞过程的研究更少。有限元法是一种计算机数值仿真方法，近年来广泛应用于预测农产品内部机械损伤^[6-7]。徐立章等^[8]建立了稻谷3层椭球体结构CAD模型，在前处理软件Hyper Mesh中划分网格，设置单元、材料属性和接触碰撞参数，生成了稻谷的有限元模型，利用LS -DYNA软件仿真分析稻谷与钉齿碰撞过程。陈燕等^[9]为减小荔枝在收获、储运过程中的机械损伤，对荔枝的宏观和微观力学特性进行研究，运用有限元法建立荔枝压缩力学模型，比较荔枝垂直和水平受压的试验值和仿真值。张荣荣等^[10]根据板栗几何尺寸和物理参数，用ANSYS建立板栗的有限元模型。通过文献梳理可知，有限元可仿真分析农业物料碰撞过程，能快速、直观地反映碰撞过程中内部应力及接触力的变化情况^[11-12]，但是目前针对木薯种茎碰撞方面的研究报道较少。

因此，在广西特色农作物自动化装备研究团队前期预切种式木薯排种机构研究^{[2-3, 13]}的基础上，本文利用三维扫描技术逆向建立木薯种茎三维模型，通过基于Hyper Mesh和LS-DYNA的种茎碰撞有限元分析，明确在不同试验因素水平下的种茎种芽−播种部件斜面接触等效应力、种芽−茎秆交接处等效应力、种芽−茎秆交接处应变及种芽−斜面接触应变的变化规律，通过二次旋转正交组合试验，寻求较小损伤精密播种方式下的最优参数。以期为揭示木薯种茎碰撞机制及播种装置优化等提供依据。

1.   材料与方法

1.1   材料与仪器

选择广西地区种植面积最大的‘桂热4号’(广西亚热带作物研究所提供)为试验对象，挑选茎秆直挺、茎皮和种芽完好的种茎(平均含水率68.46%，w)为原型进行逆向建模，如图1所示。为精确建立种茎模型，采用手持式三维CCD扫描仪采集种茎三维点云数据^{[1, 14]}，该扫描仪测量速度为5.5×10⁵次/s，测量精度为0.1 mm。

图  1  木薯种茎结构图

1：种芽，2：韧皮层，3：表皮，4：木质层，5：髓体

Figure  1.  Structure diagram of cassava seed stem

1: Bud, 2: Bast, 3: Epidermis, 4: Xylem, 5: Medulla

下载: 全尺寸图片 幻灯片

1.2   木薯种茎三维实体模型逆向重建

考虑木薯种茎实际尺寸，扫描时将角点布置在种茎四周的载物台面上，并保证每次扫描4个以上角点，如图2所示。

图  2  木薯种茎三维实体模型逆向建模过程

Figure  2.  Reverse modeling process of cassava seed stem 3D solid model

下载: 全尺寸图片 幻灯片

木薯种茎包括“主茎秆+种芽”，主茎秆木质部平均外径27 mm、韧皮部平均厚度1 mm、髓体外径9 mm。种芽为圆台体，高度(h)为5 mm，顶圆与底圆半径(r、R)分别为2.5和3.0 mm，如图3所示。种芽呈螺旋状分布，如图4所示，平均螺距为60 mm。由于种茎茎秆芯部髓体柔软，不能承受载荷，建模过程中将茎秆简化为仅由木质部和韧皮部黏连形成的空心柱体。

图  3  种茎主茎秆(A)和种芽(B)几何结构参数

Figure  3.  Parameters of the geometrical structure of the main culm (A) and the bud (B) of the seed stem

下载: 全尺寸图片 幻灯片

图  4  木薯种茎实物凸起种芽分布(A)及其建模(B)

Figure  4.  Distribution of raised buds in physical cassava seed stem (A) and modeling (B)

下载: 全尺寸图片 幻灯片

1.3   有限元模型建立

1.3.1   单元类型与材料属性设定

将生成模型保存为IGES格式，导入到Hyper Mesh后建立木薯种茎有限元模型，主茎秆的木质部取其单元类型为SHELL163。种芽和韧皮层在种茎物理结构上融为一体，取其单元为SOLIDl64。为避免产生沙漏现象^[10]，本文在*SECTION_SOLID面板的ELFORM中选用代号为2的多点高斯分算法为求解积分算法，同时将播种部件单元类型设为SOLIDl64。木薯种茎各部分相关材料属性如表1所示。

表  1  木薯种茎各部位材料参数

Table  1.  Material parameters of cassava seed stem

部位 Part	密度($ \rho $)/ (×10−10 t·mm−3) Density	径向弹性 模量($ {E_a} $)/MPa Radial elasticity modulus	轴向弹性 模量($ {E_c} $)/MPa Axial elasticity modulus	同性平面 泊松比($ {\mu _{ab}} $) Homogeneous plane Poisson’s ratio	异性平面 泊松比($ {\mu _{ac}} $) Anisotropic plane Poisson’s ratio	轴向剪切模 量($ {G_{ab}} $)/MPa Axial shear modulus	径向剪切模量 ($ {G_{bc}} $)/MPa Radial shear modulus
木质部 Xylem	8.3	50.01	25.04	0.42	0.35	17.60	92.75
韧皮部 Phloem	5.1	1.78	12.24	0.38	0.31	0.68	4.67
种芽 Seed bud	7.5	45.00	45.00	0.35	0.35	2.50	2.50

下载: 导出CSV 

| 显示表格

1.3.2   网格划分

考虑木薯种茎的生物学形态特征，采用3D_solid map面板的网格划分工具将二维网格分别沿茎秆的径向和种芽的轴向进行映射，形成三维的六面体网格，如图5A所示。木薯种茎网格划分优化处理过程如图5B所示，经测算统计，有限元模型共由11 652个单元、13 608个节点组成。

图  5  木薯种茎网格划分过程

Figure  5.  Process of meshing cassava seed stems

下载: 全尺寸图片 幻灯片

1.3.3   接触设置

考虑两物体A和B的接触问题，它们当前构形记为$ {V_{\text{A}}} $和$ {V_{\text{B}}} $，边界面为$ {\varOmega _{\text{A}}} $和$ {\varOmega _{\text{B}}} $，接触面记为$ {\varOmega _{\text{C}}} = {\varOmega _{\text{A}}} \cap {\varOmega _{\text{B}}} $，如图6所示。为分析方便，将播种振动板设为固定，种茎以一定高度跌落到振动板，与之碰撞后分离。接触类型为AUTOMATIC_SURFACE_TO_SURFACE，考虑摩擦)，接触罚因子为0.1，静摩擦系数F_S和动摩擦系数F_D设为0.15。接触界面能由*DATABASE_SLEOUT输出，同时在GLSTAT中也包含总体接触界面能。

图  6  物体接触面定义

A和B为两物体，$ {V_{\text{A}}} $和$ {V_{\text{B}}} $为构形，$ {\varOmega _{\text{A}}} $和$ {\varOmega _{\text{B}}} $为边界面，$ {\varOmega _{\text{C}}} $为接触面

Figure  6.  Definition of object contact surface

A and B are two objects, with conformations of $ {V_{\text{A}}} $ and $ {V_{\text{B}}} $, boundary surfaces of $ {\varOmega _{\text{A}}} $ and $ {\varOmega _{\text{B}}} $, contact surface of $ {\varOmega _{\text{C}}} $

下载: 全尺寸图片 幻灯片

A与B接触时的非嵌入条件可表示为：

式(1)表明，物体A与B不能互相重叠，由于事先无法确定两物体在哪一点接触，因此大变形问题中无法将非嵌入条件表示成位移的代数或微分方程，只能在每一时步对比$ {{\mathit{\Omega}} _{\text{C}}} $面上物体A和B对应节点的坐标，或对比速率来实现位移协调条件

式中，下标N表示接触法线方向，$ U_N^{\text{A}} $为接触法线方向的A点位移量，$ U_N^{\text{B}} $为接触法线方向的B点位移量，$ {n^{\text{A}}} $为A点的位移代数增量，$ {\mu ^{\text{A}}} $为A点接触能耗系数，$ {\mu ^{\text{B}}} $为B点接触能耗系数。

两物体接触面力必须服从动量平衡，由于界面上没有质量，这要求两物体上面力的合力（T）为0。

在法线方向上，不考虑接触表面之间的黏性，由作用力和反作用力定律得

式中，$T_N^{\rm{A}}$和$T_N^{\rm{B}}$分别为物体A和B的法向接触力。

2.   有限元数值模拟与结果分析

有限元仿真数值模拟研究分为2步。第1步，明确不同试验因素水平下种茎种芽−斜面接触等效应力、种芽−茎秆交接处等效应力、种芽−茎秆交接处应变及种芽−斜面接触应变的变化规律；第2步，通过二次旋转正交组合试验，寻求最优参数组合。

2.1   试验指标与试验因素

2.1.1   试验指标

采用种茎种芽−播种部件斜面接触应力Y₁、种芽−茎秆连接处应力Y₂、种芽−茎秆连接处应变Y₃及种芽−斜面接触应变Y₄为试验指标进行计算，如图7所示。

图  7  试验指标示意图

Figure  7.  Schematic diagram of test indexes

下载: 全尺寸图片 幻灯片

2.1.2   试验因素

根据前期研究基础^{[2-3, 13]}，预切种振动供种式木薯播种器如图8所示。工作时，将一定量木薯种茎投至喂种箱，通过调控种量控制提升板高度，保证一定供种量；通过阶梯式振动散种机构下方振动系统，使种茎群在振动散种板和阶梯式调姿板上往复运动，通过阶梯式调姿板卡槽限位作用，实现种茎调姿定向；接着后方输送链在电动机驱动下带动捞种勺转动，种茎在种茎群间挤压力及自身重力作用下进入捞种勺，完成充种；种茎进入清种区，当捞种勺存有多根种茎，刮种部件将多余的种茎清除；种茎随捞种勺继续运动，在护种部件作用下运动至投种口，依靠自身重力掉落，完成精密播种。

图  8  预切种振动供种式木薯播种器示意图

1：喂种箱，2：种量控制提升板，3：阶梯式振动散种机构，4：振动散种板，5：振动系统，6：阶梯式调姿板，7：充种箱，8：充种板，9：输送链，10：捞种勺，11：护种部件，12：刮种部件，13：机架，14：木薯种茎群，15：电动机，16：种层挡板

Figure  8.  Diagram of the whole machine structure of the pre-cut seed vibrating cassava seeder

1: Seed feeding box, 2: Seed volume control lifting plate, 3: Stepped vibration seed dispersal mechanism, 4: Vibration seed dispersal plate, 5: Vibration system, 6: Stepped posturing plate, 7: Seed charging box, 8: Seed charging plate, 9: Conveyor chain, 10: Seed scooping spoon, 11: Seed guarding part, 12: Seed scraping part, 13: Rack, 14: Cassava seed stem cluster, 15: Electric motor, 16: Seed layer baffle plate

下载: 全尺寸图片 幻灯片

从预切种振动供种式木薯播种器结构来说，入料高度、振动频率、振幅、种茎与播种部件碰撞接触角度是影响碰撞试验指标的主要因素。但实际现场中，播种器结构形式、尺寸和种茎物料特性是确定和不可更改的，要想减轻播种过程中碰撞损伤，需从可更改的工艺参数入手。因此，以播种部件(振动板)振动频率f (40~85 Hz)、振动板安装倾斜角度$ \alpha $ (0~27°)及种茎跌落高度$ h $ (60~240 mm)作为试验因素和各因素参数选择范围。

2.2   木薯种茎与振动板碰撞全过程分析

以振动板振动频率为60 Hz、振动板安装倾斜角度为12°、种茎跌落高度为140 mm为例，进行种茎撞击全过程及速度影响分析。如图9A，1号种芽在0.029 954 ms时与振动板产生接触，而后种茎由于惯性及重力作用继续运动，振动频率为60 Hz振动板在0~20.8 ms内沿Z轴正向做速度为0.72 m/s的匀速直线运动，由于振动板为刚性体，种茎在与振动板接触后产生压溃变形。如图9B，在1.089 9 ms时编号为61631的单元处产生最大应力。接着受到振动板撞击力作用种茎反弹，种芽处应力逐渐减小。振动板继续向上运动，种茎又受到重力作用，整个茎秆产生一种沿板面往下滚动的运动趋势。如图9C，在1.869 9 ms时2号种芽与振动板发生碰撞，伴随变形和压溃吸能。但由于第1次碰撞已经把部分的动能通过压溃吸能转化为种茎内部能量，因此第2次碰撞的受损程度远远小于第1次的碰撞，因此本处仅讨论1号种芽在第2次碰撞之前(0~1.869 9 ms)的情况。

图  9  种茎撞击全过程应力变化分析

Figure  9.  Analysis of stress changes during the whole process of seed stem impaction

下载: 全尺寸图片 幻灯片

如图10，不同视角下展现木薯茎秆与振动板接触冲击的过程中种芽处的应力云图，碰撞发生开始，种芽上表皮尖端与振动板发生接触，应力集中在接触点，随时间变化应力逐渐向种芽内部扩散，种芽处也开始出现压溃变形且接触面积逐步增大。

图  10  木薯茎秆与振动板接触冲击的过程中种芽处的应力云图

Figure  10.  Stress cloud at seed sprouts during impact of cassava stalks in contact with vibration plate

下载: 全尺寸图片 幻灯片

木薯茎秆与振动板接触冲击过程中整体单元结果位移云图如图11所示，位移越小说明碰撞中产生的压溃变形越大，分析发现随着时间增加，压溃变形逐渐由种芽处过渡到种芽和茎秆连接处再逐步扩展到茎秆中段，使其产生弯曲，变形结果如图12所示。

图  11  木薯茎秆与振动板接触冲击过程中整体单元结果位移云图

Figure  11.  Overall unit resultant displacement cloud during impact of cassava stalks in contact with vibration plate

下载: 全尺寸图片 幻灯片

图  12  木薯种茎变形结果

Figure  12.  Results of cassava seed stem deformation

下载: 全尺寸图片 幻灯片

2.3   不同因素水平对试验指标的影响

为探明不同试验因素水平对试验指标的影响规律，以振动板振动频率60 Hz、振动板安装倾斜角度12°、种茎跌落高度140 mm为基准，每次模拟试验改变单个参数，找出单因素参数对试验指标的影响，每组试验重复3次，结果取平均值。跌落高度(X₁)取60、80、100、120、140、160、180、200、220、240 mm，振动板安装倾斜角度(X₂)取0、3、6、9、12、15、18、21、24、27°，振动板振动频率(X₃)取40、45、50、55、60、65、70、75、80、85 Hz。

经探明不同跌落高度、振动板安装倾斜角度及振动板振动频率下的种茎碰撞损伤规律，确定各工作参数较优水平区间为跌落高度120~180 mm、振动板安装倾斜角度15~24°及振动板振动频率为55~70 Hz。

2.4   低损伤播种方式参数组合试验与分析

为确定较小损伤播种方式下的最优参数组合，采用二次旋转正交组合试验进行。

2.4.1   试验方案与结果

试验因素编码水平如表2所示，试验方案与结果如表3所示。

表  2  试验因素与水平

Table  2.  Test factor and level

水平 Level	种茎跌落高度(X1)/mm Seed stem drop height	振动板安装倾斜角度(X2)/(°) Vibration plate mounting tilt angle	振动板振动频率(X3)/Hz Vibration plate vibration frequency
−1.682	120.00	15.00	55.00
−1	132.16	16.82	58.04
0	150.00	19.50	62.50
1	167.84	22.18	66.96
1.682	180.00	24.00	70.00

下载: 导出CSV 

| 显示表格

表  3  试验方案与结果

Table  3.  Test plan and result

序号	X1/mm	X2/(°)	X3/Hz	Y1/MPa	Y2/MPa	Y3	Y4
1	132.16	16.82	58.04	24.31	8.21	0.045	0.776
2	167.84	16.82	58.04	24.43	8.31	0.046	0.787
3	132.16	22.18	58.04	28.21	11.42	0.063	0.912
4	167.84	22.18	58.04	29.71	12.61	0.069	0.957
5	132.16	16.82	66.96	26.44	10.12	0.056	0.891
6	167.84	16.82	66.96	29.22	16.22	0.091	0.943
7	132.16	22.18	66.96	24.91	9.42	0.052	0.813
8	167.84	22.18	66.96	32.21	17.23	0.095	1.044
9	120.00	19.50	62.50	25.72	8.74	0.048	0.831
10	180.00	19.50	62.50	32.43	16.12	0.087	1.052
11	150.00	15.00	62.50	27.63	10.32	0.057	0.891
12	150.00	24.00	62.50	30.14	11.41	0.063	0.973
13	150.00	19.50	55.00	26.62	8.24	0.045	0.861
14	150.00	19.50	70.00	31.16	14.42	0.081	1.011
15	150.00	19.50	62.50	28.63	10.23	0.056	0.923
16	150.00	19.50	62.50	29.42	11.24	0.062	0.947
17	150.00	19.50	62.50	29.14	10.04	0.055	0.883
18	150.00	19.50	62.50	29.14	10.82	0.061	0.94
19	150.00	19.50	62.50	28.74	10.42	0.057	0.934
20	150.00	19.50	62.50	28.82	10.72	0.059	0.933
21	150.00	19.50	62.50	29.42	11.11	0.061	0.951
22	150.00	19.50	62.50	29.26	10.61	0.058	0.937
23	150.00	19.50	62.50	29.42	11.24	0.062	0.947

下载: 导出CSV 

| 显示表格

2.4.2   回归模型建立与显著性检验

采用Design-Expert进行模型和系数显著性检验，结果如下。

1) Y₁回归方程为：

经方差分析，回归模型通过显著性检验，P<0.01，模型有效。由式(5)可知，Y₁与X₁、X₂和X₃都分别呈开口向下的二次函数关系。各试验因素对Y₁的贡献率排序为X₁>X₂>X₃。

2) Y₂回归方程为：

经方差分析，回归模型通过显著性检验，P<0.01，模型有效。各试验因素对Y₂的贡献率排序为X₁>X₃>X₂。

3) Y₃回归方程为：

经方差分析，回归模型通过显著性检验，P<0.01，模型有效。各试验因素对Y₃的贡献率排序为X₁>X₃>X₂。

4) Y₄回归方程为：

经方差分析，回归模型通过显著性检验，P<0.01，模型有效。各试验因素对Y₄的贡献率排序为X₁>X₃>X₂。

响应面如图13、14所示，如图13A，固定X₃为62.50 Hz，当X₁为158.92~167.84 mm和X₂为20.84~22.18°时，Y₁存在最优值。固定X₁，随着X₂的增加，Y₁增加；固定X₂，随着X₁增加，Y₁增加。如图13B所示，固定X₁为150.00 mm，当X₂为20.84~22.18°和X₃为64.73~66.96 Hz时，Y₁存在最优值。固定X₂，随着X₃增加，Y₁增加；固定X₃，随着X₂增加，Y₁增加，同理分析图14。

图  13  各因素对种茎种芽−播种部件斜面接触的等效应力的影响

Figure  13.  Effects of various factors on the equivalent stress of the contact between the seed stem bud and sowing part on the slope

下载: 全尺寸图片 幻灯片

图  14  各因素对种芽−茎秆交接处的等效应力的影响

Figure  14.  Effects of various factors on the equivalent stress at the junction of seed bud and stalk

下载: 全尺寸图片 幻灯片

2.4.3   参数优化与仿真试验验证

1)参数优化。利用Design-Expert软件中Optimization模块对建立的回归方程模型做最优化求解，得到最优工作参数，各试验因素的边界条件：

结果表明：当X₁为167.83 mm、X₂为22.18°、X₃为66.96 Hz，Y₁为32.64 MPa、Y₂为17.08 MPa、Y₃为0.094及Y₄为1.049。

2)仿真试验验证。为对优化结果进行验证，在Hyper Mesh软件中设置相关参数，将参数优化中的最优参数组合导入LS -DYNA中进行解算，比较各仿真试验指标与理论优化值。Y₁的试验值(32.64 MPa)与模型优化值(33.56 MPa)的相对误差为2.82%、Y₂的试验值(17.08 MPa)与模型优化值(16.16 MPa)的相对误差为5.39%、Y₃的试验值(0.094)与模型优化值(0.087)的相对误差为7.29%以及Y₄的试验值(1.049)与模型优化值(1.025)的相对误差为2.29%，试验值与模型优化值的相对误差均小于7.50%，因此，上述参数优化模型是准确的。

本研究得出，在对播种器进一步优化设计中，要进行木薯种茎在机械结构作用下的运动和碰撞动力学分析，优化种茎运动轨迹；研究结果为减少碰撞并降低碰撞力提供了重要的参考，同时为深入分析木薯种茎的碰撞破坏模型和破坏机理、得出造成种茎破坏的主要因素提供借鉴。

3.   结论

1)通过开展单因素试验，分析工作参数对Y₁、Y₂、Y₃及Y₄的影响规律，确定各参数较优水平区间跌落高度为120~180 mm、振动板安装倾斜角度为15~24°及振动板振动频率为55~70 Hz。

2)在单因素研究基础上，开展二次旋转正交组合试验。得到各因素对Y₁的贡献率为X₁>X₂>X₃；各因素对Y₂的贡献率为X₁>X₃>X₂；各因素对Y₃的贡献率为X₁>X₃>X₂；各因素对Y₄的贡献率为X₁>X₃>X₂。

3)当X₁为167.83 mm、X₂为22.18°、X₃为66.96 Hz，Y₁为32.64 MPa、Y₂为17.08 MPa、Y₃为0.094及Y₄为1.049，模型预测结果与实际仿真结果相近，证明了回归模型的可靠性。