基于離散元法的玉米脫粒過程分析

付 宏，呂 游，李艷雙，于建群

（1.吉林大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，長春130012；2.吉林大學(xué)生物與農(nóng)業(yè)工程學(xué)院，長春130022）

我國年產(chǎn)玉米近一億五千萬噸，其中大部分玉米都是收獲晾曬后，再進(jìn)行脫粒操作。因此，玉米脫粒機(jī)的使用非常廣泛。但由于脫粒過程的復(fù)雜性，到目前為止，國內(nèi)外對玉米脫粒過程的研究，大都采用試驗方法、統(tǒng)計分析方法或傳統(tǒng)連續(xù)介質(zhì)力學(xué)的分析方法［1－3］。試驗方法和統(tǒng)計分析方法費時費力，所得結(jié)果一般也不具有普遍意義，且不能洞察玉米脫粒過程的物理機(jī)理。傳統(tǒng)連續(xù)介質(zhì)力學(xué)的分析方法，只能分析脫粒過程中單個玉米穗、玉米芯或玉米籽粒的受力和運(yùn)動情況，或把玉米穗、玉米芯和玉米籽粒當(dāng)成一個整體來分析，這與實際脫粒過程中，玉米穗、玉米芯和玉米籽粒群體的相互作用和運(yùn)動過程差別較大。

本文采用顆粒動力學(xué)方法——離散元法研究玉米的脫粒過程，建立了一種分析玉米脫粒過程的新方法和軟件。

1 玉米穗和脫粒機(jī)的建模方法

玉米穗由玉米芯和生長在芯上的籽粒組成，且通常玉米穗成圓臺形，而玉米籽粒截面近似梯形，如圖1所示。本文采用球顆粒聚合方法，建立了玉米穗和玉米籽粒的分析模型如圖2所示。

圖1 實際玉米穗和玉米籽粒形狀Fig.1 Real shapes of a corn ear and a corn grain

圖2 玉米穗和玉米籽粒分析模型Fig.2 Analysis models of a corn ear and a corn grain

為分析玉米穗的脫粒過程，建立圖3所示的全局坐標(biāo)系x-y-z和局部坐標(biāo)系X-Y-Z。其中全局坐標(biāo)系是大地坐標(biāo)系，局部坐標(biāo)系的原點與玉米穗的質(zhì)心重合，X軸與玉米穗軸線重合，且與某一玉米籽粒中心重合，Z軸方向符合右手螺旋法則。為了模擬玉米籽粒與玉米芯之間的連接，采用一個球（稱連接球）模擬籽粒果柄，且籽粒果柄與玉米籽粒是一個整體，如圖4所示，玉米籽粒脫落后，連接球被刪除。玉米芯組成單元之間的連接也采用同樣的連接球方法處理。

圖3 全局坐標(biāo)系和局部坐標(biāo)系Fig.3 Global coordinate system and local coordinate system

在分析玉米脫粒過程時，首先在全局坐標(biāo)系的玉米穗入口隨機(jī)產(chǎn)生一點，作為待生成玉米穗的質(zhì)心，然后以該質(zhì)心為坐標(biāo)原點，產(chǎn)生局部坐標(biāo)系X-Y-Z并生成一個玉米穗，下一時步再隨機(jī)產(chǎn)生一點，并用同樣方法生成另一個玉米穗，如此反復(fù)。其中玉米籽粒分析模型的高度、厚度和上、下底寬度尺寸、玉米穗分析模型的長度和大小頭直徑，均可通過程序界面由輸入?yún)?shù)控制。

圖4 玉米芯、玉米籽粒和籽粒果柄Fig.4 Corn cob，corn grain and kernel stem

玉米脫粒機(jī)的建模方法采用本項目組提出的基于CAD模型的邊界建模方法［4］，即把與玉米穗接觸的脫粒機(jī)中零件表面離散成圖元，如平面、球面、柱面、錐面等，并添加運(yùn)動屬性和材料特性，由此建立脫粒機(jī)的離散元法分析模型。以滾筒式玉米脫粒機(jī)為例，其三維CAD模型和三維離散元法分析模型如圖5和圖6所示。

圖5 滾筒式玉米脫粒機(jī)的三維CAD模型Fig.5 3D CAD model of a drum type corn sheller

圖6 滾筒式玉米脫粒機(jī)的三維離散元法分析模型Fig.6 Analysis model of a drum type corn sheller for 3D DEM

2 脫粒過程的作用力計算方法

玉米脫粒過程中的作用力包括籽粒與籽粒之間、籽粒與玉米芯之間以及籽粒與脫粒部件之間的作用力，還包括玉米芯與玉米芯之間和玉米芯與脫粒部件之間的作用力。在采用顆粒聚合體方法建立玉米穗分析模型時，當(dāng)籽粒未脫落前，籽粒與其生長的玉米芯之間的作用力為連接力（吸引力），同一玉米芯組成單元之間的作用力也為連接力（吸引力）；而籽粒與其他籽粒之間、籽粒與其生長的玉米芯以外的玉米芯之間、籽粒與脫粒部件之間的作用力為接觸斥力；玉米芯與其他玉米芯間、玉米芯與脫粒部件之間的作用力也為接觸斥力。由于玉米籽粒和玉米芯的建模方法均采用球顆粒聚合的方法，因此計算籽粒和玉米芯上的作用力，可采用球顆粒的計算方法，這樣就避免了求解非線性方程組。下面以某一玉米籽粒P為例（見圖4），介紹接觸作用力的計算方法。

首先在全局坐標(biāo)系x-y-z中建立如圖7所示的局部坐標(biāo)系X′-Y′-Z′。其中X′軸為玉米籽粒P組成球i的球心與另一個接觸籽粒或玉米芯組成球j的球心連線，Y′軸平行于全局坐標(biāo)系x-y平面，Z′軸符合右手螺旋法則。在X′-Y′-Z′局部坐標(biāo)系中，當(dāng)i球與j球的中心距小于兩球的半徑之和時，則兩球接觸，其法向疊合量δn為

式中：Ri、Rj、（xi，yi，zi）和（xj，yj，zj）分別為球i和j的半徑及中心點在全局坐標(biāo)系x-y-z中的坐標(biāo)。當(dāng)采用線性黏彈性力學(xué)模型計算兩球的接觸斥力時，兩球間的法向接觸作用力為

圖7 兩顆粒組成球接觸時的局部坐標(biāo)系Fig.7 Local coordinate system when two balls contact

當(dāng)球i與脫粒部件某一表面接觸時，如以平面為例，局部坐標(biāo)系X′軸為球i的球心到該平面的垂線，Y′和Z′軸的確定方法及法向和切向接觸作用力的計算方法與上述相同。

采用線性黏彈性力學(xué)模型計算球i與球j的連接力的方法與上述相同，但剛度系數(shù)和阻尼系數(shù)選取不同，計算斥力時選取接觸剛度系數(shù)和接觸阻尼系數(shù)，計算連接力時選取連接剛度系數(shù)和連接阻尼系數(shù)。而且在籽粒沒脫落前，不采用庫侖－莫爾準(zhǔn)則對切向力進(jìn)行修正。

當(dāng)采用上述方法計算某籽粒與其生長的玉米芯之間的法向及切向連接力和玉米芯組成單元之間的法向及切向連接力時，還需分別與由試驗測定的連接作用力最大值作比較，如果計算的連接力大于試驗測定的連接力最大值，則玉米籽粒從玉米芯上脫離——脫粒，或玉米芯組成單元之間分離——破碎。

當(dāng)采用上述方法計算作用在玉米穗、脫落的玉米籽?；蚱扑榈挠衩仔締卧厦總€組成球的法向和切向作用力后，還需把作用力分別移到玉米穗、玉米籽粒或玉米芯單元的質(zhì)心上，并變換到全局坐標(biāo)系下，然后再求出作用在玉米穗、玉米籽?；蚱扑橛衩仔締卧|(zhì)心上的合力和合力矩。

3 顆粒運(yùn)動計算方法

在玉米脫粒過程中，顆粒運(yùn)動包括玉米穗的運(yùn)動、脫落籽粒的運(yùn)動、破碎玉米芯單元的運(yùn)動和未脫落籽粒相對玉米芯的運(yùn)動。

通過上面的分析計算后，可求得作用在玉米穗、脫落玉米籽粒和破碎玉米芯單元質(zhì)心上的在t時刻下沿全局坐標(biāo)系3個坐標(biāo)軸方向所受的合力，繼而可得t＋Δt時刻顆粒質(zhì)心在全局坐標(biāo)系3個坐標(biāo)軸的新位置 （x（t＋Δt），y（t＋Δt），z（t＋Δt））為

玉米穗、脫落玉米籽粒和破碎玉米芯單元在t時刻的轉(zhuǎn)動角速度滿足如下歐拉動力學(xué)方程［5］

求解上述方程即可得到顆粒在局部坐標(biāo)系X-Y-Z中的角位移及新的局部坐標(biāo)系，限于篇幅求解過程從略。

圖8 計算玉米籽粒轉(zhuǎn)動的局部坐標(biāo)系Fig.8 Local coordinate system when calculated corn grain rotating

計算未脫落籽粒相對于玉米芯的運(yùn)動，還需建立局部坐標(biāo)系X″-Y″-Z″，參見圖3和圖4。其中X″軸與X軸重合，Y″軸為被分析玉米籽粒質(zhì)心與玉米芯中心連線，當(dāng)被分析籽粒為P時，Y″軸與Y軸重合，當(dāng)分析其他籽粒時，Y″軸需變化；Z″軸符合右手螺旋法則。在X″-Y″-Z″局部坐標(biāo)系中，求出籽粒上作用的合外力，并求出籽粒與玉米芯的連接力，根據(jù)合外力與連接力即可求出未脫落籽粒相對于玉米芯的新位置并對其更新。

通過上述計算方法按時步迭代，即可求出玉米穗、脫落玉米籽粒、破碎玉米芯單元和未脫落籽粒相對于玉米芯在每一時刻所受的力、力矩、運(yùn)動速度和新位置，由此即可分析玉米的脫粒過程。

4 軟件設(shè)計和實例驗證

在上述工作的基礎(chǔ)上，本文設(shè)計了玉米脫粒過程分析軟件，并實現(xiàn)了與CAD軟件的集成，從而開發(fā)出一種集設(shè)計和性能分析評價為一體的玉米脫粒過程分析和脫粒機(jī)優(yōu)化設(shè)計軟件［6－7］，其結(jié)構(gòu)如圖9所示。其流程為由CAD軟件進(jìn)行玉米脫粒機(jī)的設(shè)計，然后由邊界建模模塊自動建立脫粒機(jī)的離散元法分析模型，接著即可進(jìn)行脫粒過程的離散元法分析計算，計算結(jié)束后可根據(jù)計算結(jié)果文件進(jìn)行脫粒過程的性能分析，如脫凈率、籽粒損失率和籽粒沿脫粒機(jī)軸向分布等，還可以動態(tài)顯示玉米的脫粒過程。

圖9 玉米脫粒過程和脫粒機(jī)優(yōu)化設(shè)計新方法及其軟件的結(jié)構(gòu)Fig.9 Proposed method and software structure for analysis of corn threshing process and optimization of thresher

圖10和圖11為采用上述方法及其軟件，由CAD模型實現(xiàn)的兩種結(jié)構(gòu)滾筒式玉米脫粒機(jī)工作過程的二維離散元法仿真分析。圖12為采用上述方法及其軟件，由圖5所示玉米滾筒式脫粒機(jī)的三維CAD模型實現(xiàn)的玉米脫粒過程的三維離散元法仿真分析。從圖12可以看出，隨著仿真計算時間的延續(xù)，玉米籽粒逐漸從玉米穗上脫落并從脫粒機(jī)凹板上的孔漏出。圖13為軟件給出的破碎率、脫凈率和脫落籽粒沿脫粒機(jī)軸線分布的分析結(jié)果。由此初步證明了本文方法和軟件的可行性。有關(guān)仿真和試驗的詳細(xì)對比分析將另文撰述。

圖10 由CAD模型實現(xiàn)玉米脫粒過程的二維仿真分析Fig.10 2D DEM simulation and analysis by CAD model

圖11 改變CAD模型實現(xiàn)的玉米脫粒過程二維仿真分析Fig.11 2D DEM simulation and analysis by changing CAD model

本文方法和軟件的優(yōu)點是：①在設(shè)計階段通過修改脫粒機(jī)的CAD模型、玉米籽粒、玉米芯和玉米穗的分析模型、連接和接觸作用的力學(xué)模型、離散元法計算參數(shù)，能分析不同品種玉米、不同工況、不同原理、不同結(jié)構(gòu)和尺寸的玉米脫粒機(jī)的性能，由此實現(xiàn)玉米脫粒機(jī)結(jié)構(gòu)方案和尺寸參數(shù)的優(yōu)化；②通過脫粒機(jī)的CAD模型，能進(jìn)行脫粒過程的動態(tài)仿真，由此分析玉米脫粒機(jī)的工作機(jī)理或工作過程，還可以發(fā)明新原理和新結(jié)構(gòu)的玉米脫粒機(jī)，這是現(xiàn)有玉米脫粒機(jī)的研究和設(shè)計方法不能做到的。

圖12 玉米脫粒過程的三維離散元法仿真分析Fig.12 3D DEM simulation and analysis of corn threshing process

圖13 本文軟件給出的脫粒過程分析結(jié)果Fig.13 Analysis results of threshing process by software

5 結(jié)束語

根據(jù)玉米脫粒機(jī)研究和設(shè)計中存在的問題，本文提出采用離散元法研究玉米的脫粒過程，由脫粒機(jī)的CAD模型建立了脫粒機(jī)的離散元法分析模型，采用顆粒聚合體方法建立了玉米穗的分析模型，采用接觸和連接力學(xué)模型計算脫粒過程中的作用力，研究了玉米脫粒過程的離散元法計算方法，在此基礎(chǔ)上建立了一種分析玉米脫粒過程的新方法和軟件。通過對該軟件的實例驗證，初步證明了新方法的可行性，為玉米脫粒過程的分析和脫粒機(jī)的優(yōu)化提供了一種新思路。

［1］Petkevichius S，Shpokas L，Kutzbach H D.Investigation of the maize ear threshing process［J］.Biosystems Engineering，2008，99：532-539.

［2］Miu P I，Kutzbach H D.Modeling and simulation of grain threshing and separation in threshing units-Part I［J］.Computers and Electronics in Agriculture，2008，60：96-104.

［3］毛欣，衣淑娟.軸流脫粒、分離機(jī)理及仿真研究的發(fā)展探討［J］.農(nóng)機(jī)化研究，2007，10：202-205.

Mao Xin，Yi Shu-juan.The developmental discussion of mechanism research and simulation of the axial flow threshing mechanism［J］.Journal of Agricultural Mechanization Research，2007，10：202-205.

［4］于建群，付宏，劉振宇，等.基于CAD模型的離散元法邊界建模方法：中國，200510016835.7［P］. 2006-07-28.

［5］Lin X，Ng T T.A three dimensional discrete element model using arrays of ellipsoids［J］.Geotechnique，1997，47（2）：319-329.

［6］李艷雙.顆粒聚合體的三維離散元法算法研究與軟件設(shè)計［D］.長春：吉林大學(xué)，2011.

Li Yan-shuang.Research of algorithm for particles agglomerate in 3D DEM and design of the software［D］.Changchun：Jilin University，2011.

［7］吳玄辰.基于二維離散元法的玉米脫粒過程仿真分析［D］.長春：吉林大學(xué)生物與農(nóng)業(yè)工程學(xué)院，2011.

Wu Xuan-chen.Simulation and analysis of maize threshing based on 2D DEM［D］.Changchun：Jilin University，College of Biological and Agricultural Engineering，2011.