基于Hyper Mesh和LS-DYNA的玉米籽粒碰撞損傷動(dòng)態(tài)過程的有限元分析1

王博，王俊，杜冬冬

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

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

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

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

1 材料與方法

1.1 材料與儀器

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

1.2 三維實(shí)體模型逆向重建

三維掃描過程主要根據(jù)玉米籽粒的外形來創(chuàng)建其幾何表面的點(diǎn)云，用于進(jìn)一步插補(bǔ)成籽粒的表面形狀。為了能夠?qū)⒉煌课蝗S點(diǎn)云數(shù)據(jù)在程序中自動(dòng)配準(zhǔn)，本文在掃描過程中采用直徑為6 mm的反光角點(diǎn)作為配準(zhǔn)基準(zhǔn)[17]；考慮到玉米籽粒的實(shí)際尺寸，掃描時(shí)將角點(diǎn)布置在玉米四周的載物臺(tái)面上，并保證每次掃描有4個(gè)以上的角點(diǎn)。

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

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

1.3 有限元模型的建立

1.3.1 單元類型與材料屬性的設(shè)定

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

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

1.3.2 接觸設(shè)置

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

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

1.3.3 單元大小選取與模型校驗(yàn)

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

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

1.4 有限元仿真分析

為了分析玉米籽粒-脫粒元件碰撞損傷的動(dòng)態(tài)過程，探究不同影響因素下玉米籽粒的損傷臨界速度，優(yōu)化碰撞參數(shù)，將不同的含水率（11.78%、17.63%、23.45%、29.31%、34.73%）和碰撞速度（4、6、8、10、12 m/s）寫入HyperMesh的不同的有限元模型中。處理完成后生成顯式動(dòng)力分析軟件LS-DYNA計(jì)算所用的輸入文件（k文件），然后將該k文件提交到LS-DYNA進(jìn)行解算，解算結(jié)果綜合利用Hyperview、Hypergraph和Ls-Prepost軟件進(jìn)行查看和分析。

1.5 參數(shù)優(yōu)化仿真試驗(yàn)

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

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

2 結(jié)果與討論

2.1 網(wǎng)格密度選取

將不同網(wǎng)格密度（0.10、0.15、0.20、0.25、0.30 mm）的玉米籽粒有限元模型利用LS-DYNA進(jìn)行解算，得到的不同密度下的最大接觸力、最大Von Mises應(yīng)力及最大變形結(jié)果如圖3所示?？梢钥闯觯畲蠼佑|力、最大Von Mises應(yīng)力及最大變形均隨著網(wǎng)格的不斷細(xì)化而逐漸逼近精確解，在單元大小為0.15 mm時(shí)得到的最大接觸力、最大Von Mises應(yīng)力及最大變形分別為36.94 N、53.05 MPa和0.85 mm，與單元大小為0.10 mm的模型解算結(jié)果的最大值僅相差2%，但卻大大縮短了結(jié)算時(shí)間。因此，本文選取網(wǎng)格大小為0.15 mm的有限元模型進(jìn)行計(jì)算和分析。

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

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

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

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

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

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

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

在收獲期間碰撞速度為8 m/s時(shí)不同含水率（11.78%、17.63%、23.45%、29.31%、34.73%）下的最大Von Mises應(yīng)力解算云圖和含水率為23.45%時(shí)不同碰撞速度（2、4、6、8、10、12 m/s）下的最大 Von Mises應(yīng)力解算云圖結(jié)果如圖5所示，經(jīng)統(tǒng)計(jì)后的最大應(yīng)力變化曲線如圖6所示?？梢钥闯觯涸谙嗤氏拢S著碰撞速度增大，玉米籽粒的最大Von Mises應(yīng)力值逐漸增大；在相同碰撞速度下，隨著含水率增大，最大Von Mises應(yīng)力逐漸減小。這與相關(guān)文獻(xiàn)中的試驗(yàn)結(jié)果[7-9,11,15]一致。將不同含水率下玉米籽粒的強(qiáng)度極限帶入圖6中擬合出的方程，得到各含水率下的碰撞臨界速度結(jié)果，如表3所示。

圖6 不同含水率和碰撞速度下的最大Von Mises應(yīng)力變化曲線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īng)面分析及參數(shù)優(yōu)化

根據(jù)正交中心組合設(shè)計(jì)的兩因素五水平正交試驗(yàn)，試驗(yàn)方案包括13個(gè)試驗(yàn)點(diǎn)，其中含8個(gè)分析因子，5個(gè)零點(diǎn)估計(jì)誤差。試驗(yàn)方案和響應(yīng)值見表4。

表4 試驗(yàn)設(shè)計(jì)方案和響應(yīng)值分析結(jié)果Table 4 Analysis result of experimental design and response values

2.4.1 響應(yīng)回歸模型的建立與顯著性檢驗(yàn)

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

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

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

表5 試驗(yàn)結(jié)果統(tǒng)計(jì)分析Table 5 Statistical analysis for experimental results

圖7 交互因素及水平對(duì)響應(yīng)指標(biāo)的影響Fig.7 Effects of interactive factors and levels on response indexes

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

2.4.2.1 參數(shù)優(yōu)化

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

式中：Xi為各因素，Yjmin為響應(yīng)的最小值。

2.4.2.2 仿真試驗(yàn)驗(yàn)證

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

3 結(jié)論

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

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

圖8 仿真試驗(yàn)驗(yàn)證結(jié)果Fig.8 Simulation test results

表6 優(yōu)化條件下各評(píng)價(jià)指標(biāo)的仿真試驗(yàn)值與預(yù)測(cè)值Table 6 Simulative and predictive values of each evaluation index under optimized condition

2）最大接觸力、最大Von Mises應(yīng)力隨含水率的增大而減小，隨碰撞速度增大而增大，最大變形則相反。玉米籽粒在含水率為11.78%、17.63%、23.45%、29.31%、34.73%時(shí)的臨界損傷速度分別為5.51、6.75、8.15、9.36、10.57 m/s。

3）在單因素研究的基礎(chǔ)上，開展響應(yīng)面試驗(yàn)研究，建立了二次多項(xiàng)式回歸模型，分析了含水率和碰撞速度交互項(xiàng)對(duì)最大接觸力、最大Von Mises應(yīng)力和最大變形的影響，得到各因素對(duì)最大接觸力Y1、最大應(yīng)力Y2和最大變形Y3的影響大小排序，分別為碰撞速度＞含水率、含水率＞碰撞速度和碰撞速度＞含水率。

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

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