玉米籽粒收获，即指在田间一次完成玉米摘穗、脱粒、清选等作业，直接获得玉米籽粒的过程，具有收获周期短、生产成本低等优点，是今后我国玉米机械化收获的发展方向^[1-3]。然而，在籽粒收获过程中，玉米籽粒与机械相互碰撞，造成大量破损，严重影响玉米的经济价值^[4-5]。因此，研究玉米籽粒碰撞损伤的动态过程，揭示碰撞损伤的机制，对设计低损伤玉米籽粒收获装置、提高收获效率及推进玉米生产全程机械化发展都具有重要的现实意义。

谷物的碰撞损伤受很多因素的影响，目前国内外关于谷物碰撞损伤的研究主要通过理论^[6-9]及试验的方法^[10-11]进行，这些研究证明了含水率及碰撞速度对谷物的碰撞损伤具有显著影响，并初步得到了影响规律。然而，受研究方法的限制，此类研究都缺少对谷物碰撞过程的动态分析，无法揭示碰撞损伤过程中应力、应变及能量等关键参数的变化。

有限单元法是一种计算机数值仿真分析方法，近年来开始被应用于预测农产品内部机械损伤^[12-13]，包括谷物的受载损伤。李心平等^[14]将玉米籽粒简化为平面模型，利用二维有限元法预测了静载下的籽粒内部应力分布；XU等^{[5, 15]}将水稻籽粒简化为椭球模型，研究了籽粒与钉齿碰撞时籽粒内的应力分布。然而，以往针对谷物的有限元研究都将模型做了大量简化，但谷物的几何特征是有限元分析的重要参数，过度简化模型会影响到解算结果的精度；不仅如此，含水率作为影响谷物碰撞损伤的最重要的因素之一，目前还未见借助有限元方法进行深入研究的报道。因此，在现有研究的基础上，本文进一步利用三维扫描技术获取玉米籽粒的三维点云数据，逆向建立玉米籽粒的三维真实实体模型，并建立其与脱粒元件碰撞的有限元模型，以此分析不同含水率的玉米籽粒在不同速度碰撞过程中内部应力、接触力及能量等参数的变化，计算其碰撞损伤临界速度，进而从损伤最小的角度优化碰撞参数，为揭示玉米碰撞损伤机制及收获装置的优化设计等提供依据。

1 材料与方法 1.1 材料与仪器

本文以形状规则、无外表损伤的玉米品种郑单958（半马齿型，千粒质量307 g）的扁平籽粒为原型进行逆向建模（扁平籽粒约占该品种玉米籽粒的90% ^[16]）。该籽粒长9.22 mm，宽6.70 mm，厚4.48 mm。采用加拿大Creaform公司生产的手持式三维CCD（charge coupled device，电荷耦合器件）扫描仪采集玉米籽粒的三维点云数据，该扫描仪测量速度为5.5×10⁵次/s，扫描区域为143 mm×108 mm，景深为100 mm，测量精度为0.1 mm。

1.2 三维实体模型逆向重建

三维扫描过程主要根据玉米籽粒的外形来创建其几何表面的点云，用于进一步插补成籽粒的表面形状。为了能够将不同部位三维点云数据在程序中自动配准，本文在扫描过程中采用直径为6 mm的反光角点作为配准基准^[17]；考虑到玉米籽粒的实际尺寸，扫描时将角点布置在玉米四周的载物台面上，并保证每次扫描有4个以上的角点。

逆向重建过程如图 1A所示：利用VXelement三维扫描软件和反光角点辅助配准的方法配准不同部位的点云，提取玉米籽粒的三维点云数据；将该点云数据保存为obj格式后导入逆向工程软件Geomagic Studio 2012进行降噪、分离和封装等处理，得到玉米籽粒的曲面网格线条组模型；再将该线条组模型转换为iges格式后导入三维建模软件Siemens NX 8.5中，通过曲线组依次创建片体并进行缝合，最终得到可用于建立有限元模型的玉米籽粒的三维实体模型。本文重建的玉米籽粒三维实体模型及尺寸如图 1B所示。

1.3 有限元模型的建立 1.3.1 单元类型与材料属性的设定

前处理工作是有限元分析的重要环节，针对玉米籽粒复杂的几何特征，本文使用HyperMesh软件进行前处理。将上文建立的籽粒三维实体模型保存为iges格式，导入到前处理软件HyperMesh中建立玉米籽粒的碰撞有限元模型。

设置正确合理的单元类型和材料属性是进行有限元仿真分析的前提。在HyperMesh中，针对实体“Volume”的六面体和四面体单元都具有很高的精度，为了避免六面体网格在复杂曲面仿真时出现网格退化，影响计算精度等问题，本文选择适用于复杂几何体的四面体网格对玉米籽粒进行网格划分。由于不同含水率籽粒的力学特性参数差异显著，因此需要设置不同的材料参数。本文将玉米籽粒视为各向同性的线性弹性，在HyperMesh中通过“Collectors/Creat/Mats”来创建、存储和管理材料的弹性模量、泊松比、密度等参数，使用的单位系统为：长度，mm；时间，s；质量，t；应力，N；弹性模量，MPa。玉米及元件的相关材料属性如表 1所示^[18-21]。

1.3.2 接触设置

为了分析方便，本文将脱粒元件钉齿设为固定，玉米籽粒以一定的速度与脱粒元件发生碰撞。在设置接触时以脱粒元件表面为从面，籽粒表面为主面，接触类型为AUTOMATIC_SURFACE_TO_SURFACE，取接触惩罚因子为0.1，脱粒元件与籽粒的摩擦系数为0.35^[16]。根据不同的仿真试验条件设定时间步长、终止时间等控制卡片参数。

1.3.3 单元大小选取与模型校验

单元大小直接影响到有限元分析的最终结果，理论上单元越小越能逼近真实值，但过小的单元会使得有限元解算时间成倍增加。同时，为了验证有限元分析结果的可靠性，对模型进行能量校验十分必要。因此，本文在开展系统仿真前利用所建立的籽粒模型进行了预分析。如图 2所示，本文将玉米籽粒模型用不同尺寸（0.30、0.25、0.20、0.15、0.10 mm）的单元进行划分，在有限元分析后提取最大接触力、最大Von Mises应力及最大变形^[22]，对比响应值随网格大小的变化情况，用于选取合适的网格大小，同时解算出碰撞过程中的能量变化，对其进行评估以校验模型^[23]。

1.4 有限元仿真分析

为了分析玉米籽粒-脱粒元件碰撞损伤的动态过程，探究不同影响因素下玉米籽粒的损伤临界速度，优化碰撞参数，将不同的含水率（11.78%、17.63%、23.45%、29.31%、34.73%）和碰撞速度（4、6、8、10、12 m/s）写入HyperMesh的不同的有限元模型中。处理完成后生成显式动力分析软件LS-DYNA计算所用的输入文件（k文件），然后将该k文件提交到LS- DYNA进行解算，解算结果综合利用Hyperview、Hypergraph和Ls-Prepost软件进行查看和分析。

1.5 参数优化仿真试验

为了优化籽粒与脱粒元件的碰撞参数，寻求较小损伤的碰撞方式，仿真试验选取含水率X₁、碰撞速度X₂为影响因素，设计两因素五水平的正交旋转中心组合优化试验^[24-26]，因素水平编码如表 2所示，仿真试验共进行13组。试验方案设计及结果分析应用Design-Expert 10.0软件完成。

2 结果与讨论 2.1 网格密度选取

将不同网格密度（0.10、0.15、0.20、0.25、0.30 mm）的玉米籽粒有限元模型利用LS-DYNA进行解算，得到的不同密度下的最大接触力、最大Von Mises应力及最大变形结果如图 3所示。可以看出，最大接触力、最大Von Mises应力及最大变形均随着网格的不断细化而逐渐逼近精确解，在单元大小为0.15 mm时得到的最大接触力、最大Von Mises应力及最大变形分别为36.94 N、53.05 MPa和0.85 mm，与单元大小为0.10 mm的模型解算结果的最大值仅相差2%，但却大大缩短了结算时间。因此，本文选取网格大小为0.15 mm的有限元模型进行计算和分析。

2.2 碰撞动态行为分析及模型可靠性检验

在不同仿真设置下，碰撞参数值差异较大，但总体变化趋势相似。为了分析和研究玉米-脱粒元件碰撞的动态过程，以含水率23.45%、碰撞速度8 m/s的有限元模型解算结果为例，玉米-脱粒元件碰撞过程中籽粒截面的Von Mises应力云图及Von Mises应力、接触力、能量变化如图 4所示。可以看出，籽粒与脱粒元件的碰撞过程大致可分为碰撞前、碰撞中和碰撞后3个阶段。解算时间0~58.59 μs为碰撞前阶段，如图 4A（a）、图 4B和图 4C所示，该阶段籽粒与脱粒元件之间没有发生接触，籽粒的应力和接触力均为0，此时能量全部以动能的形式储存在玉米籽粒中。解算时间58.59~137.54 μs是碰撞中阶段，其中：在解算时间58.59~113.17 μs过程中，籽粒从开始接触脱粒元件[图 4A（a）]到最大变形处[图 4A（c）]，该过程中籽粒内部的Von Mises应力、接触力均逐渐增大；在任一时刻[以图 4A（b）中所示85.88 μs为例]，接触区域中心处的Von Mises应力值最大，沿四周扩散并逐渐减小，该过程中玉米籽粒的动能逐渐转化为内能（主要是弹性势能），同时，由于籽粒和脱粒元件间的摩擦，一部分籽粒动能转换为界面滑移能。因此，该阶段即为玉米籽粒损伤产生的阶段，损伤区域为接触的应力集中区域。如图 4A（e）、图 4B和图 4C所示，解算时间137.54~ 161.17 μs为碰撞后阶段，在该过程中，玉米籽粒的应力、接触力均逐渐减小，玉米籽粒在完全离开脱粒元件后内部依旧残存少量应力，籽粒的内能逐渐转化为动能，但籽粒与脱粒元件停止接触，故界面滑移基本不变。

沙漏能是检测有限元模型和结果可靠性的重要依据，在进行有限元分析时，一般要求沙漏能不超过5%~10%^{[23, 27]}。由图 4C可以看出，在本文的有限元分析过程中，沙漏能最大时仅占总能量的4.67%，说明本文的有限元模型是可靠的。

2.3 单因素分析及临界速度判定

在收获期间碰撞速度为8 m/s时不同含水率（11.78%、17.63%、23.45%、29.31%、34.73%）下的最大Von Mises应力解算云图和含水率为23.45%时不同碰撞速度（2、4、6、8、10、12 m/s）下的最大Von Mises应力解算云图结果如图 5所示，经统计后的最大应力变化曲线如图 6所示。可以看出：在相同含水率下，随着碰撞速度增大，玉米籽粒的最大Von Mises应力值逐渐增大；在相同碰撞速度下，随着含水率增大，最大Von Mises应力逐渐减小。这与相关文献中的试验结果^{[7-9, 11, 15]}一致。将不同含水率下玉米籽粒的强度极限带入图 6中拟合出的方程，得到各含水率下的碰撞临界速度结果，如表 3所示。

2.4 响应面分析及参数优化

根据正交中心组合设计的两因素五水平正交试验，试验方案包括13个试验点，其中含8个分析因子，5个零点估计误差。试验方案和响应值见表 4。

2.4.1 响应回归模型的建立与显著性检验

根据表 4中的试验设计及结果，运用DesignExpert 10.0软件进行回归拟合，建立最大接触力Y₁、最大应力Y₂和最大变形Y₃对籽粒含水率X₁、碰撞速度X₂的二次多项式回归模型，回归模型的显著性检验如表 5所示，剔除不显著因素后的回归模型如式（1）~（3）所示。

$ \begin{array}{l} {Y_1} = 36.56 - 6.78{X_1} + 11.59{X_2} - 3.56{X_1}{X_2} + 1.68{X_1}^2 - \\ 0.21{X_2}^2. \end{array} $

(1)

$ \begin{array}{l} {Y_2} = 52.42 - 19.52{X_1} + 4.96{X_2} - 6.93{X_1}{X_2} + 5.30{X_1}^2\\ - 0.41{X_2}^2. \end{array} $

(2)

$ {Y_3} = 0.87 + 0.069{X_1} + 0.082{X_2} - 0.033{X_1}{X_2}. $

(3)

由表 5可知：最大接触力Y₁、最大应力Y₂和最大变形Y₃的响应面模型均在α=0.01水平上显著；Y₁、Y₂和Y₃模型的决定系数分别为0.98、0.98和0.96，回归模型对样本点的拟合程度较好。因此，可用此模型对最大接触力Y₁、最大应力Y₂和最大变形Y₃进行分析和预测。由表 5及式(1)~(3)可知，各因素对最大接触力Y₁、最大应力Y₂和最大变形Y₃的贡献率大小排序分别为：碰撞速度＞含水率、含水率＞碰撞速度和碰撞速度＞含水率。

如图 7所示，运用Design-Expert 10.0软件绘制响应曲面图来分析各交互因素对响应指标的影响效应。由表 5中的方差分析可知，对于每一个响应，X₁X₂均在α=0.01水平上显著。如图 7A、7B所示：当碰撞速度X₂过高而含水率X₁过低时，最大接触力Y₁和最大Von Mises应力Y₂出现最大值；当碰撞速度最小而含水率最高时，上述2个响应出现最小值。如图 7C所示：当含水率和碰撞速度均最大时，最大变形Y₃出现最大值；当含水率和碰撞速度最小时，最大变形Y₃出现最小值。

2.4.2 参数优化与仿真试验验证 2.4.2.1 参数优化

通过上述分析可知，各因素及其交互作用对预测指标的影响趋势不同。为了寻求最佳参数组合，需要综合考虑各因素对预测指标的影响，进行多目标参数优化。如式（4）所示，以最小Von Mises应力、最小接触力及最小变形为优化目标，利用Design-Expert 10.0软件对建立的3个回归方程模型作最优化求解，得到最佳参数组合为X₁=0.603，X₂=－1.414，对应的实际含水率为26.99%，碰撞速度为5.17 m/s，此时最大接触力Y₁=19.30 N，最大应力Y₂=40.64 MPa，最大变形Y₃= 0.73 mm。

$ \left\{ \begin{array}{l} - 1.414 < {X_i} < 1.414{\rm{ }}(i = 1, 2)\\ {Y_{j{\rm{min}}}} > 0{\rm{ }}(j = 1, 2, 3) \end{array} \right. $

(4)

式中：X_i为各因素，Y_{j min}为响应的最小值。

2.4.2.2 仿真试验验证

为了对优化结果进行验证，在HyperMesh软件中设置相关参数，将2.4.2.1节中的最优参数组合（含水率26.99%，碰撞速度5.17 m/s）按上文所述方法将k文件导入LS-DYNA中进行解算，得到的结果如图 8所示，各仿真试验指标与理论优化值对比如表 6所示。从中可知，最大接触力Y₁的试验值与模型优化值的相对误差为4.77%，最大应力Y₂的试验值与模型优化值的相对误差为7.17%，最大变形Y₃的试验值与模型优化值的相对误差为1.37%，即试验值与模型优化值的相对误差均小于8%。因此，上述参数优化模型是准确的。采用上述最优工作参数组合能有效降低最大接触力、最大应力和最大变形，即碰撞速度为5.17 m/s、含水率为26.99%，此时最大接触力为19.30 N，最大应力为40.64 MPa，最大变形为0.73 mm。

3 结论

本文建立了玉米籽粒的三维实体模型，对籽粒碰撞的动态过程进行了有限元分析，研究了不同含水率、不同碰撞速度下的籽粒损伤情况，确定了不同含水率籽粒的碰撞损伤临界速度，优化了碰撞参数。得到的主要结论如下：

1）玉米籽粒在碰撞时，应力集中在较小的接触区域，接触区域受压应力作用，应力沿四周扩散并逐渐减小，区域中心处应力值最大，损伤在接触区域产生，并沿应力扩散方向发展。

2）最大接触力、最大Von Mises应力随含水率的增大而减小，随碰撞速度增大而增大，最大变形则相反。玉米籽粒在含水率为11.78%、17.63%、23.45%、29.31%、34.73%时的临界损伤速度分别为5.51、6.75、8.15、9.36、10.57 m/s。

3）在单因素研究的基础上，开展响应面试验研究，建立了二次多项式回归模型，分析了含水率和碰撞速度交互项对最大接触力、最大Von Mises应力和最大变形的影响，得到各因素对最大接触力Y₁、最大应力Y₂和最大变形Y₃的影响大小排序，分别为碰撞速度＞含水率、含水率＞碰撞速度和碰撞速度＞含水率。

4）多目标优化得到最佳碰撞参数组合为含水率26.99%，碰撞速度5.17 m/s，对应的最大接触力Y₁=19.30 N，最大应力Y₂=40.64 MPa，最大变形Y₃= 0.73 mm，与实际仿真试验结果误差小于8%，证明了回归模型的可靠性。

总之，本文对揭示玉米脱粒损伤机制、优化设计脱粒装置、选择较优的脱粒装置工作参数具有一定意义。