王博，王俊，杜冬冬

（浙江大學生物系統(tǒng)工程與食品科學學院，杭州310058）

玉米籽粒收獲，即指在田間一次完成玉米摘穗、脫粒、清選等作業(yè)，直接獲得玉米籽粒的過程，具有收獲周期短、生產成本低等優(yōu)點，是今后我國玉米機械化收獲的發(fā)展方向[1-3]。然而，在籽粒收獲過程中，玉米籽粒與機械相互碰撞，造成大量破損，嚴重影響玉米的經濟價值[4-5]。因此，研究玉米籽粒碰撞損傷的動態(tài)過程，揭示碰撞損傷的機制，對設計低損傷玉米籽粒收獲裝置、提高收獲效率及推進玉米生產全程機械化發(fā)展都具有重要的現實意義。

谷物的碰撞損傷受很多因素的影響，目前國內外關于谷物碰撞損傷的研究主要通過理論[6-9]及試驗的方法[10-11]進行，這些研究證明了含水率及碰撞速度對谷物的碰撞損傷具有顯著影響，并初步得到了影響規(guī)律。然而，受研究方法的限制，此類研究都缺少對谷物碰撞過程的動態(tài)分析，無法揭示碰撞損傷過程中應力、應變及能量等關鍵參數的變化。

有限單元法是一種計算機數值仿真分析方法，近年來開始被應用于預測農產品內部機械損傷[12-13]，包括谷物的受載損傷。李心平等[14]將玉米籽粒簡化為平面模型，利用二維有限元法預測了靜載下的籽粒內部應力分布；XU等[5,15]將水稻籽粒簡化為橢球模型，研究了籽粒與釘齒碰撞時籽粒內的應力分布。然而，以往針對谷物的有限元研究都將模型做了大量簡化，但谷物的幾何特征是有限元分析的重要參數，過度簡化模型會影響到解算結果的精度；不僅如此，含水率作為影響谷物碰撞損傷的最重要的因素之一，目前還未見借助有限元方法進行深入研究的報道。因此，在現有研究的基礎上，本文進一步利用三維掃描技術獲取玉米籽粒的三維點云數據，逆向建立玉米籽粒的三維真實實體模型，并建立其與脫粒元件碰撞的有限元模型，以此分析不同含水率的玉米籽粒在不同速度碰撞過程中內部應力、接觸力及能量等參數的變化，計算其碰撞損傷臨界速度，進而從損傷最小的角度優(yōu)化碰撞參數，為揭示玉米碰撞損傷機制及收獲裝置的優(yōu)化設計等提供依據。

1 材料與方法

1.1 材料與儀器

本文以形狀規(guī)則、無外表損傷的玉米品種鄭單958（半馬齒型，千粒質量307 g）的扁平籽粒為原型進行逆向建模（扁平籽粒約占該品種玉米籽粒的90%[16]）。該籽粒長9.22 mm，寬6.70 mm，厚4.48 mm。采用加拿大Creaform公司生產的手持式三維CCD（charge coupled device，電荷耦合器件）掃描儀采集玉米籽粒的三維點云數據，該掃描儀測量速度為5.5×105次/s，掃描區(qū)域為143 mm×108 mm，景深為100 mm，測量精度為0.1 mm。

1.2 三維實體模型逆向重建

三維掃描過程主要根據玉米籽粒的外形來創(chuàng)建其幾何表面的點云，用于進一步插補成籽粒的表面形狀。為了能夠將不同部位三維點云數據在程序中自動配準，本文在掃描過程中采用直徑為6 mm的反光角點作為配準基準[17]；考慮到玉米籽粒的實際尺寸，掃描時將角點布置在玉米四周的載物臺面上，并保證每次掃描有4個以上的角點。

逆向重建過程如圖1A所示：利用VXelement三維掃描軟件和反光角點輔助配準的方法配準不同部位的點云，提取玉米籽粒的三維點云數據；將該點云數據保存為obj格式后導入逆向工程軟件Geomagic Studio 2012進行降噪、分離和封裝等處理，得到玉米籽粒的曲面網格線條組模型；再將該線條組模型轉換為iges格式后導入三維建模軟件Siemens NX 8.5中，通過曲線組依次創(chuàng)建片體并進行縫合，最終得到可用于建立有限元模型的玉米籽粒的三維實體模型。本文重建的玉米籽粒三維實 體模型及尺寸如圖1B所示。

圖1 玉米籽粒的三維重建過程（A）及實體模型（B）Fig.1 Three dimensional reconstruction process(A)and real model(B)of maize kernel

1.3 有限元模型的建立

1.3.1 單元類型與材料屬性的設定

前處理工作是有限元分析的重要環(huán)節(jié)，針對玉米籽粒復雜的幾何特征，本文使用HyperMesh軟件進行前處理。將上文建立的籽粒三維實體模型保存為iges格式，導入到前處理軟件HyperMesh中建立玉米籽粒的碰撞有限元模型。

設置正確合理的單元類型和材料屬性是進行有限元仿真分析的前提。在HyperMesh中，針對實體“Volume”的六面體和四面體單元都具有很高的精度，為了避免六面體網格在復雜曲面仿真時出現網格退化，影響計算精度等問題，本文選擇適用于復雜幾何體的四面體網格對玉米籽粒進行網格劃分。由于不同含水率籽粒的力學特性參數差異顯著，因此需要設置不同的材料參數。本文將玉米籽粒視為各向同性的線性彈性，在HyperMesh中通過“Collectors/Creat/Mats”來創(chuàng)建、存儲和管理材料的彈性模量、泊松比、密度等參數，使用的單位系統(tǒng)為：長度，mm；時間，s；質量，t；應力，N；彈性模量，MPa。玉米及元件的相關材料屬性如表1所示[18-21]。

1.3.2 接觸設置

為了分析方便，本文將脫粒元件釘齒設為固定，玉米籽粒以一定的速度與脫粒元件發(fā)生碰撞。在設置接觸時以脫粒元件表面為從面，籽粒表面為主面，接觸類型為AUTOMATIC_SURFACE_TO_SURFACE，取接觸懲罰因子為0.1，脫粒元件與籽粒的摩擦系數為0.35[16]。根據不同的仿真試驗條件設定時間步長、終止時間等控制卡片參數。

表1 有限元模型的材料屬性Table1 Material propertiesof finiteelement models

1.3.3 單元大小選取與模型校驗

單元大小直接影響到有限元分析的最終結果，理論上單元越小越能逼近真實值，但過小的單元會使得有限元解算時間成倍增加。同時，為了驗證有限元分析結果的可靠性，對模型進行能量校驗十分必要。因此，本文在開展系統(tǒng)仿真前利用所建立的籽粒模型進行了預分析。如圖2所示，本文將玉米籽粒模型用不同尺寸（0.30、0.25、0.20、0.15、0.10 mm）的單元進行劃分，在有限元分析后提取最大接觸力、最大Von Mises應力及最大變形[22]，對比響應值隨網格大小的變化情況，用于選取合適的網格大小，同時解算出碰撞過程中的能量變化，對其進行評估以校驗模型[23]。

圖2 玉米籽粒的網格劃分密度Fig.2 Mesh partitioning density of maizekernel

1.4 有限元仿真分析

為了分析玉米籽粒-脫粒元件碰撞損傷的動態(tài)過程，探究不同影響因素下玉米籽粒的損傷臨界速度，優(yōu)化碰撞參數，將不同的含水率（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 參數優(yōu)化仿真試驗

為了優(yōu)化籽粒與脫粒元件的碰撞參數，尋求較小損傷的碰撞方式，仿真試驗選取含水率X1、碰撞速度X2為影響因素，設計兩因素五水平的正交旋轉中心組合優(yōu)化試驗[24-26]，因素水平編碼如表2所示，仿真試驗共進行13組。試驗方案設計及結果分析應用Design-Expert 10.0軟件完成。

表2 因素水平編碼表Table2 Factorsand codelevelsof tests

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的有限元模型進行計算和分析。

圖3 不同網格大小的玉米籽粒有限元模型仿真結果Fig.3 Simulation results of finite element size on the stress progression of maizekernel

2.2 碰撞動態(tài)行為分析及模型可靠性檢驗

在不同仿真設置下，碰撞參數值差異較大，但總體變化趨勢相似。為了分析和研究玉米-脫粒元件碰撞的動態(tài)過程，以含水率23.45%、碰撞速度8 m/s的有限元模型解算結果為例，玉米-脫粒元件碰撞過程中籽粒截面的Von Mises應力云圖及Von Mises應力、接觸力、能量變化如圖4所示?？梢钥闯觯蚜Ｅc脫粒元件的碰撞過程大致可分為碰撞前、碰撞中和碰撞后3個階段。解算時間0～58.59μs為碰撞前階段，如圖4A（a）、圖4B和圖4C所示，該階段籽粒與脫粒元件之間沒有發(fā)生接觸，籽粒的應力和接觸力均為0，此時能量全部以動能的形式儲存在玉米籽粒中。解算時間58.59～137.54μs是碰撞中階段，其中：在解算時間58.59～113.17μs過程中，籽粒從開始接觸脫粒元件[圖4A（a）]到最大變形處[圖4A（c）]，該過程中籽粒內部的Von Mises應力、接觸力均逐漸增大；在任一時刻[以圖4A（b）中所示85.88μs為例]，接觸區(qū)域中心處的Von Mises應力值最大，沿四周擴散并逐漸減小，該過程中玉米籽粒的動能逐漸轉化為內能（主要是彈性勢能），同時，由于籽粒和脫粒元件間的摩擦，一部分籽粒動能轉換為界面滑移能。因此，該階段即為玉米籽粒損傷產生的階段，損傷區(qū)域為接觸的應力集中區(qū)域。如圖4A（e）、圖4B和圖4C所示，解算時間137.54～161.17μs為碰撞后階段，在該過程中，玉米籽粒的應力、接觸力均逐漸減小，玉米籽粒在完全離開脫粒元件后內部依舊殘存少量應力，籽粒的內能逐漸轉化為動能，但籽粒與脫粒元件停止接觸，故界面滑移基本不變。

沙漏能是檢測有限元模型和結果可靠性的重要依據，在進行有限元分析時，一般要求沙漏能不超過5%～10%[23,27]。由圖4C可以看出，在本文的有限元分析過程中，沙漏能最大時僅占總能量的4.67%，說明本文的有限元模型是可靠的。

2.3 單因素分析及臨界速度判定

圖4 仿真動態(tài)過程的可視化結果Fig.4 Visual result of simulation dynamics process

圖5 不同含水率和碰撞速度下的最大Von Mises應力分布圖Fig.5 Contour plot of maximum Von Misesstressunder different moisturecontentsand impact velocities

在收獲期間碰撞速度為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所示，經統(tǒng)計后的最大應力變化曲線如圖6所示?？梢钥闯觯涸谙嗤氏拢S著碰撞速度增大，玉米籽粒的最大Von Mises應力值逐漸增大；在相同碰撞速度下，隨著含水率增大，最大Von Mises應力逐漸減小。這與相關文獻中的試驗結果[7-9,11,15]一致。將不同含水率下玉米籽粒的強度極限帶入圖6中擬合出的方程，得到各含水率下的碰撞臨界速度結果，如表3所示。

圖6 不同含水率和碰撞速度下的最大Von Mises應力變化曲線Fig.6 Change curve of the maximum Von Mises stress under different moisturecontentsand impact velocities

表3 不同含水率的玉米籽粒的損傷臨界速度判定Table 3 Critical velocities of maize kernel with different moisturecontents

2.4 響應面分析及參數優(yōu)化

根據正交中心組合設計的兩因素五水平正交試驗，試驗方案包括13個試驗點，其中含8個分析因子，5個零點估計誤差。試驗方案和響應值見表4。

表4 試驗設計方案和響應值分析結果Table 4 Analysis result of experimental design and response values

2.4.1 響應回歸模型的建立與顯著性檢驗

根據表4中的試驗設計及結果，運用Design-Expert 10.0軟件進行回歸擬合，建立最大接觸力Y1、最大應力Y2和最大變形Y3對籽粒含水率X1、碰撞速度X2的二次多項式回歸模型，回歸模型的顯著性檢驗如表5所示，剔除不顯著因素后的回歸模型如式（1）～（3）所示。

由表5可知：最大接觸力Y1、最大應力Y2和最大變形Y3的響應面模型均在α=0.01水平上顯著；Y1、Y2和Y3模型的決定系數分別為0.98、0.98和0.96，回歸模型對樣本點的擬合程度較好。因此，可用此模型對最大接觸力Y1、最大應力Y2和最大變形Y3進行分析和預測。由表5及式(1)～(3)可知，各因素對最大接觸力Y1、最大應力Y2和最大變形Y3的貢獻率大小排序分別為：碰撞速度＞含水率、含水率＞碰撞速度和碰撞速度＞含水率。

如圖7所示，運用Design-Expert 10.0軟件繪制響應曲面圖來分析各交互因素對響應指標的影響效應。由表5中的方差分析可知，對于每一個響應，X1X2均在α=0.01水平上顯著。如圖7A、7B所示：當碰撞速度X2過高而含水率X1過低時，最大接觸力Y1和最大Von Mises應力Y2出現最大值；當碰撞速度最小而含水率最高時，上述2個響應出現最小值。如圖7C所示：當含水率和碰撞速度均最大時，最大變形Y3出現最大值；當含水率和碰撞速度最小時，最大變形Y3出現最小值。

表5 試驗結果統(tǒng)計分析Table 5 Statistical analysis for experimental results

圖7 交互因素及水平對響應指標的影響Fig.7 Effects of interactive factors and levels on response indexes

2.4.2 參數優(yōu)化與仿真試驗驗證

2.4.2.1 參數優(yōu)化

通過上述分析可知，各因素及其交互作用對預測指標的影響趨勢不同。為了尋求最佳參數組合，需要綜合考慮各因素對預測指標的影響，進行多目標參數優(yōu)化。如式（4）所示，以最小Von Mises應力、最小接觸力及最小變形為優(yōu)化目標，利用Design-Expert 10.0軟件對建立的3個回歸方程模型作最優(yōu)化求解，得到最佳參數組合為 X1=0.603，X2=－1.414，對應的實際含水率為26.99%，碰撞速度為5.17 m/s，此時最大接觸力Y1=19.30 N，最大應力Y2=40.64 MPa，最大變形Y3=0.73 mm。

式中：Xi為各因素，Yjmin為響應的最小值。

2.4.2.2 仿真試驗驗證

為了對優(yōu)化結果進行驗證，在HyperMesh軟件中設置相關參數，將2.4.2.1節(jié)中的最優(yōu)參數組合（含水率26.99%，碰撞速度5.17 m/s）按上文所述方法將k文件導入LS-DYNA中進行解算，得到的結果如圖8所示，各仿真試驗指標與理論優(yōu)化值對比如表6所示。從中可知，最大接觸力Y1的試驗值與模型優(yōu)化值的相對誤差為4.77%，最大應力Y2的試驗值與模型優(yōu)化值的相對誤差為7.17%，最大變形Y3的試驗值與模型優(yōu)化值的相對誤差為1.37%，即試驗值與模型優(yōu)化值的相對誤差均小于8%。因此，上述參數優(yōu)化模型是準確的。采用上述最優(yōu)工作參數組合能有效降低最大接觸力、最大應力和最大變形，即碰撞速度為5.17 m/s、含水率為26.99%，此時最大接觸力為19.30 N，最大應力為40.64 MPa，最大變形為0.73 mm。

3 結論

本文建立了玉米籽粒的三維實體模型，對籽粒碰撞的動態(tài)過程進行了有限元分析，研究了不同含水率、不同碰撞速度下的籽粒損傷情況，確定了不同含水率籽粒的碰撞損傷臨界速度，優(yōu)化了碰撞參數。得到的主要結論如下：

1）玉米籽粒在碰撞時，應力集中在較小的接觸區(qū)域，接觸區(qū)域受壓應力作用，應力沿四周擴散并逐漸減小，區(qū)域中心處應力值最大，損傷在接觸區(qū)域產生，并沿應力擴散方向發(fā)展。

圖8 仿真試驗驗證結果Fig.8 Simulation test results

表6 優(yōu)化條件下各評價指標的仿真試驗值與預測值Table 6 Simulative and predictive values of each evaluation index under optimized condition

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應力和最大變形的影響，得到各因素對最大接觸力Y1、最大應力Y2和最大變形Y3的影響大小排序，分別為碰撞速度＞含水率、含水率＞碰撞速度和碰撞速度＞含水率。

4）多目標優(yōu)化得到最佳碰撞參數組合為含水率26.99%，碰撞速度5.17 m/s，對應的最大接觸力Y1=19.30 N，最大應力Y2=40.64 MPa，最大變形Y3=0.73 mm，與實際仿真試驗結果誤差小于8%，證明了回歸模型的可靠性。

總之，本文對揭示玉米脫粒損傷機制、優(yōu)化設計脫粒裝置、選擇較優(yōu)的脫粒裝置工作參數具有一定意義。