基于離散元方法的糧食自動分級研究

徐登峰　成榮　朱煜　等

摘要：以淺圓倉中垂直落料的大豆及雜質為研究對象，以軟球假設為基礎建立顆粒的接觸力學模型，采用離散元方法對自動分級過程進行分析，研究顆粒形狀、計算時間步長對結果的影響，結果表明，最終分析結果與實際扦樣試驗數據一致，離散元方法應用于糧食自動分級研究有較好的有效性，具有重要的實際意義。

關鍵詞：自動分級；離散元；接觸力學；瑞利波

中圖分類號： S226.5文獻標志碼： A文章編號：1002-1302（2015）01-0396-03

收稿日期：2014-01-14

作者簡介：徐登峰（1974—），男，四川華銀人，博士，助理研究員，主要從事超精密機械結構及動力學研究。Tel：（010） 62780925；E-mail：xudf@u-precision.com。

通信作者：成榮，博士，在清華大學材料學院從事博士后研究工作。Tel：（010）62781457；E-mail：chengr@tsinghua.edu.cn。 淺圓倉由于容量大、占地面積小等優點，在糧食儲存中得到廣泛應用。由于糧食中會摻雜很多雜質，在進倉落料及堆積過程中會產生雜質聚集現象，即自動分級現象[1-2]，這會導致糧食孔隙密度減小、帶菌量增大、易發生霉變板結等，長時間儲糧會對淺圓倉內糧食造成大量損失。因此，對淺圓倉糧食自動分級現象進行研究，避免分級現象的產生具有重要的實際意義[3-5]。目前，針對糧食自動分級現象的研究主要通過試驗手段[6-9]，如周延智等研究淺圓倉的雜質分布規律、張峻嶺等研究布糧器在淺圓倉防分級中的應用等[10-11]，這些研究方法會耗費大量的人力和物力資源，且研發周期較長。為降低成本和縮短研發周期，本試驗采用離散元方法對糧食自動分級現象進行分析研究。

1離散元方法及接觸力學模型

1.1離散元方法

離散元方法（discrete element method，DEM）是建立在牛頓第二運動定律基礎上研究非連續性顆粒物質結構和運動規律的一種數值方法，根據顆粒間重疊量計算接觸力，依次更新每個顆粒的速度和位置，進而確定性地演化整個顆粒系統[12]，其具體求解過程采用顯式解法，先利用中心差分法將運動方程進行離散，然后在每個時間步長內進行一次迭代，根據前一次迭代所得到的顆粒位置，由物理方程求出接觸力，作為下一次迭代的出發點，再用來求出顆粒的新位置，如此反復進行迭代，直至最后達到穩定流動過程為止[11]。在大豆入倉過程中，大豆、雜質和糧倉相互之間會發生接觸碰撞，采用離散元方法進行計算，可準確有效地得到大豆及雜質在任意時刻的速度和位置，最終獲得糧倉內的實時顆粒分布狀況，從而研究糧食自動分級現象。

1.2顆粒接觸力學模型

顆粒接觸理論嚴密，求解過程相對繁瑣，在不產生顯著誤差的條件下，有必要予以適當簡化處理。軟球模型是目前常用的簡化模型，把顆粒間法向力用彈簧和阻尼器模型描述，切向力用彈簧、阻尼器和滑動器描述，并引入彈性系數和阻尼系數等參量[13-14]，依據顆粒間法向重疊位移和切向位移計算接觸力，不考慮顆粒表面變形和接觸力加載歷史，計算強度較小，適合于工程問題的數值計算[15-16]。

以軟球模型假設為基礎的顆粒之間法向及切向力力學模型如圖1所示[17]，當顆粒之間的切向作用力小于或等于最大靜摩擦力時，顆粒之間受到法向作用力和切向作用力；當顆粒之間的切向作用力大于最大靜摩擦力，顆粒之間會發生相對滑動，從而受到法向作用力及滑動摩擦力。

通過作用力模型計算可得，顆粒之間法向作用力表達式為：

Fn=Knδ3/2n+Dnδ1/4nδn。（1）

式中：Fn表示顆粒之間法向作用力；δn表示顆粒之間法向重疊位移；Dn表示顆粒之間法向阻尼系數。

根據Hertz接觸理論，確定法向剛度系數Kn表達式[18]為：

Kn=431-υ2iEi+1-υ2jEj-1ai+ajaiaj-1/2。（2）

式中：υi表示顆粒i的泊松比；Ei表示顆粒i的彈性模量；υj表示顆粒j的泊松比；Ej表示顆粒j的彈性模量；αi表示顆粒i的球體半徑；αj表示顆粒j的球體半徑。

顆粒之間切向作用力為靜摩擦和動摩擦值之間的最小值，表達式為：

Ft=min（Ktδt+Dtδ1/4tδ′t，μFn）。（3）

式中：Ft表示顆粒之間切向作用力；δt表示顆粒之間切向相對位移；Dt表示顆粒之間切向阻尼力。

根據Mindlin-Dereiewicz接觸理論，確定切向剛度系數Kt表達式[15]為：

Kt=8δ1/2n1-υ2iGi+1-υ2jGj-1ai+ajaiaj-1/2。（4）

式中：Gi表示顆粒i的剪切彈性模量；Gj表示顆粒j的剪切彈性模量。

2離散元分析

2.1顆粒建模方式選擇

根據大豆和雜質實物（圖2）實際取樣統計發現，大豆主要為球形，雜質以四面體和球形為主。

以多個球形模型擬合四面體模型（圖3），分別以四面體模型和球形模型進行建模，在同等條件下進行驗證性計算，并將計算結果中的雜質分布狀況繪制云圖（圖4、圖5）。由圖4、圖5可見，雜質采用四面體建模和采用球形建模，對雜質的分布狀況影響不大，采用球形模型計算可以大大提高計算速度[15]。因此，離散元分析中顆粒均采用球形建模方式。

2.2分析模型的建立

考慮糧食自動分級的實際因素，在分析過程不考慮空氣阻力影響。按照實際糧倉尺寸進行建模，顆粒入倉方式為頂部垂直落料。根據實際大豆與雜質的質量百分比設定大豆與顆粒數目，大豆、雜質與糧倉模型的物理性質參數見表1。

表1大豆、雜質和糧倉的物理性質參數endprint

顆粒密度

（kg/m3）剪切模量

（Pa）泊松比大豆1 2681.5×1080.4雜質2 0003.1×1070.3糧倉4 0005.0×1090.33

2.3時間步長計算

2.3.1瑞利時間步長顆粒發生接觸碰撞時，總能耗的70%通過瑞利波（Rayleigh wave）而消耗。因此，可根據沿固體球形顆粒表面傳播的瑞利波波速來確定臨界時間步長[18-19]。不同顆粒組成的系統，瑞利時間步長計算公式為：

Δt=πR0.163υ+0.877ρΔGmin。（5）

式中：R表示顆粒半徑；υR表示瑞利波速；υ表示泊松比；G表示剪切模量；ρ表示密度。

兩顆粒間的接觸作用僅限于發生碰撞的兩顆粒上，而不應該通過瑞利波而傳遞到其他顆粒上，因此，時間步長應小于瑞利波傳遞半球所需要的時間。在實際計算時，依據顆粒運動劇烈程度選取合適的時間步長，以保證顆粒系統演變的計算穩定性[12，16]。

對淺圓倉內實際物料分布進行扦樣試驗，得到雜質分布圖（圖6）。由圖6可見，雜質主要分布在倉中心約15.0 m高的圓臺形區域內，超過15.0 m高度雜質數量減少，堆積在一個底面角約為40°、近似圓錐的區域內。

2.3.2計算時間步長確定設定圓臺上的表面直徑D為雜質分布直徑， 以雜質分布直徑D和計算所需時間t為評價指

標來選擇合適的計算時間步長Δts。當Δts設定為瑞利時間步長Δt的40%以上時，無法保證計算穩定性，過小，又會導致模型計算時間過長[11，20]。為確定合適的計算時間步長，令：

Δts=αΔt，α∈（0.05，0.4）。（6）

在此范圍內，采用Latin HyperCube方法設計采樣點進行模型計算（表2）。采用Radial Basis Function（RBF）擬合方法[21]，擬合計算時間和計算精度（雜質分布直徑）隨計算時間步長的變化關系。由圖7、圖8可見，計算所耗時間和計算時間步長接近線性關系，當計算時間小于瑞利時間的175%，計算精度（雜質分布直徑）趨于收斂。綜合考慮模型的計算精度及計算效率，設定計算時間步長為瑞利時間步長的17.5%。

表2不同計算時間步長下的計算效率和精度

α（%）t（h）D（m）37.242410.7031.543311.5727.793512.1822.683712.2122.273812.2317.344212.2613.285112.276.455912.28

3結果與分析

利用得到的所有參數，進行顆粒堆積離散元計算，對計算結果進行數據分析：將糧倉進行分區，提取每個分區內的雜質

及顆粒數目，計算得到雜質數目所占比例；采用拉格朗日二次插值方法對數據進行處理，得到雜質在糧倉內的分布云圖（圖9）。由圖9可見，雜質分布主要集中在糧倉的中心高約14.4 m的圓臺形區域內，超過14.4 m高度雜質數量減少，堆積在一個近似圓錐的區域內，圓錐底面直徑約12.26 m，底面角約41°，這與糧倉實際物料分布扦樣試驗的結果（圖6）基本一致。采用離散元方法計算，能夠準確有效地進行糧食自動分級現象研究，該方法可以為解決糧食自動分級現象提供有效的分析手段。

4結論

以淺圓倉中垂直落料的大豆及雜質為研究對象，建立以軟球模型假設為基礎的顆粒間法向及切向力切換力學模型，推導出作用力表達公式。在此基礎上，基于離散元方法對自動分級過程進行分析計算，詳細論述顆粒建模方法對計算速度及分析結果的影響，對計算精度影響較大的計算時間步長，引入Latin HyperCube分析方法，得到計算時間及計算精度隨步長的變化關系，確定計算時間步長值，且最終計算結果與扦樣試驗數據一致。離散元方法在糧食自動分級研究工作中具有較好的有效性。

參考文獻：

[1]楊文生，張中，張成. 淺圓倉儲糧自動分級研究現狀[J]. 糧食加工，2011，36（5）：74-76.

[2]呂新. 淺圓倉儲糧中心雜質集中問題的解決方案[J]. 糧食流通技術，2003（1）：19-20.

[3]王永昌. 淺圓倉糧食自動分級問題的討論[J]. 糧食流通及倉儲，2002（12）：11-13.

[4]李洪昌，李耀明，唐忠，等. 基于EDEM的振動篩分數值模擬與分析[J]. 農業工程學報，2011，27（5）：117-121.

[5]張來林，張愛強，朱彥，等. 糧食自動分級的類型與預防措施[J]. 糧食儲藏，2008（6）：28-30.

[6]姜永順. 幾種淺圓倉進倉設備的比較[J]. 糧食與飼料工業，2002（9）：16-17.

[7]張峻嶺，陳藝，張卓青. 布糧器在淺圓倉入糧防分級中的研究與應用[J]. 糧食流通技術，2012（6）：16-22.

[8]王永昌，姚文冠，李軍五. 淺圓倉糧食入倉自動分級原因分析及解決措施[J]. 糧食流通技術，2010（3）：28-29.

[9]程緒鐸，黃之斌，石翠霞，等. 筒倉中入糧自動分級的動力學分析及緩解措施研究[J]. 糧食儲藏，2012（1）：16-19.

[10]周延智，李松偉，曾卓. 淺圓倉雜質分布規律的研究[J]. 糧油倉儲科技通訊，2009（4）：14-16.

[11]張峻嶺，陳藝，張卓青. 布糧器在淺圓倉入糧防分級中的研究與應用[J]. 糧食流通技術，2012（6）：16-22.

[12]Richarda P，Taberlet N. Recent advances in DEM simula-dons of grains in a rotating drum[J]. Soft Matter，2008（4）：345-1348.

[13]焦玉勇，葛修潤，谷先榮. 三維離散元法中的數據結構[J]. 巖土力學，1998，19（2）：74-79.

[14]焦玉勇，葛修潤. 基于靜態松弛法求解的三維離散單元法[J]. 巖石力學與工程學報，2000，19（4）：453-458.

[15]孫其誠，王光謙. 顆粒物質力學導論[M]. 北京：科學出版社，2009：41-42.

[16]Teufelsbauer H，Wang Y，Chiou M C，et al. Flow-obstacle interaction in rapid granular avalanches：DEM simulation and comparison with experiment[J]. Granul Matter，2009，11（4）：209-220.

[17]Grima A P，Wypych P W. Investigation into calibration of discrete element model parameters of scale-up and validation of particle-structure interactions under impact conditions[J]. Power Technology，2011，212（1）：198-209.

[18]吳清松，胡茂彬. 顆粒流的動力學模型和實驗研究進展[J]. 力學進展，2002，32（2）：251-260.

[19]李婉宜，曾攀，雷麗萍，等. 離散顆粒流動堆積行為離散元模擬及實驗研究[J]. 力學與實踐，2012，34（1）：20-25.

[20]徐泳，黃文彬. 顆粒離散元法建模和仿真的若干進展[J]. 過程工程學報，2002，2（增刊）：530-536.

[21]楊洋，唐壽高. 顆粒流的離散元法模擬及其進展[J]. 中國粉體技術，2006，12（5）：38-43. 丁克奎，鐘凱文. 基于“3S”的精準農業管理系統設計與實現[J]. 江蘇農業科學，2015，43（1）：399-401.endprint