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

徐登峰　成榮　朱煜　等

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

關(guān)鍵詞：自動分級；離散元；接觸力學(xué)；瑞利波

中圖分類號： S226.5文獻(xiàn)標(biāo)志碼： A文章編號：1002-1302（2015）01-0396-03

收稿日期：2014-01-14

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

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

1離散元方法及接觸力學(xué)模型

1.1離散元方法

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

1.2顆粒接觸力學(xué)模型

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

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

通過作用力模型計(jì)算可得，顆粒之間法向作用力表達(dá)式為：

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

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

根據(jù)Hertz接觸理論，確定法向剛度系數(shù)Kn表達(dá)式[18]為：

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

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

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

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

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

根據(jù)Mindlin-Dereiewicz接觸理論，確定切向剛度系數(shù)Kt表達(dá)式[15]為：

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

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

2離散元分析

2.1顆粒建模方式選擇

根據(jù)大豆和雜質(zhì)實(shí)物（圖2）實(shí)際取樣統(tǒng)計(jì)發(fā)現(xiàn)，大豆主要為球形，雜質(zhì)以四面體和球形為主。

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

2.2分析模型的建立

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

表1大豆、雜質(zhì)和糧倉的物理性質(zhì)參數(shù)endprint

顆粒密度

（kg/m3）剪切模量

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

2.3時間步長計(jì)算

2.3.1瑞利時間步長顆粒發(fā)生接觸碰撞時，總能耗的70%通過瑞利波（Rayleigh wave）而消耗。因此，可根據(jù)沿固體球形顆粒表面?zhèn)鞑サ娜鹄úㄋ賮泶_定臨界時間步長[18-19]。不同顆粒組成的系統(tǒng)，瑞利時間步長計(jì)算公式為：

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

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

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

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

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

標(biāo)來選擇合適的計(jì)算時間步長Δts。當(dāng)Δts設(shè)定為瑞利時間步長Δt的40%以上時，無法保證計(jì)算穩(wěn)定性，過小，又會導(dǎo)致模型計(jì)算時間過長[11，20]。為確定合適的計(jì)算時間步長，令：

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

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

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

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

3結(jié)果與分析

利用得到的所有參數(shù)，進(jìn)行顆粒堆積離散元計(jì)算，對計(jì)算結(jié)果進(jìn)行數(shù)據(jù)分析：將糧倉進(jìn)行分區(qū)，提取每個分區(qū)內(nèi)的雜質(zhì)

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

4結(jié)論

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

參考文獻(xiàn)：

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

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

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

[4]李洪昌，李耀明，唐忠，等. 基于EDEM的振動篩分?jǐn)?shù)值模擬與分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2011，27（5）：117-121.

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

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

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

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

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

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

[11]張峻嶺，陳藝，張卓青. 布糧器在淺圓倉入糧防分級中的研究與應(yīng)用[J]. 糧食流通技術(shù)，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]焦玉勇，葛修潤，谷先榮. 三維離散元法中的數(shù)據(jù)結(jié)構(gòu)[J]. 巖土力學(xué)，1998，19（2）：74-79.

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

[15]孫其誠，王光謙. 顆粒物質(zhì)力學(xué)導(dǎo)論[M]. 北京：科學(xué)出版社，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]吳清松，胡茂彬. 顆粒流的動力學(xué)模型和實(shí)驗(yàn)研究進(jìn)展[J]. 力學(xué)進(jìn)展，2002，32（2）：251-260.

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

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

[21]楊洋，唐壽高. 顆粒流的離散元法模擬及其進(jìn)展[J]. 中國粉體技術(shù)，2006，12（5）：38-43. 丁克奎，鐘凱文. 基于“3S”的精準(zhǔn)農(nóng)業(yè)管理系統(tǒng)設(shè)計(jì)與實(shí)現(xiàn)[J]. 江蘇農(nóng)業(yè)科學(xué)，2015，43（1）：399-401.endprint