農(nóng)業(yè)領(lǐng)域離散元法（DEM）參數(shù)標(biāo)定現(xiàn)狀與發(fā)展

吳正陽(yáng)，謝方平,2※，梅玉茹，王修善,2，劉大為,2，鄔 備,2

（1.湖南農(nóng)業(yè)大學(xué)機(jī)電工程學(xué)院，長(zhǎng)沙 410128；2.智能農(nóng)機(jī)裝備湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410128）

引言

離散元法（DEM）是一種將介質(zhì)整體視為若干顆粒單元集合的數(shù)值模擬方法，利用牛頓第二定律和接觸模型，以每個(gè)顆粒的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)行為描述整體特征。從1979年提出至今，ITASCA公司開(kāi)發(fā)的PFC-2D/3D、DEM Solution公司推出的EDEM、ESSS公司研發(fā)的Rocky-DEM、開(kāi)源LIGGGHTS、ESyS-Particle等DEM軟件為DEM的廣泛應(yīng)用奠定了基礎(chǔ)。

DEM已被廣泛應(yīng)用于葡萄和玉米脫粒[1,2]、玉米和水稻清選[3-6]、土壤與農(nóng)機(jī)具相互作用[7-15]等農(nóng)業(yè)領(lǐng)域。由于農(nóng)業(yè)物料外形各異、尺寸不一，農(nóng)田土壤粒度分布復(fù)雜、孔隙度不均勻，DEM模擬的物料外形尺寸或土壤粒度和孔隙度存在一定的誤差。該誤差會(huì)導(dǎo)致仿真對(duì)象的物理特性與實(shí)際情況存在較大差異。為在DEM中再現(xiàn)研究對(duì)象的實(shí)際物理特性，需要確定合理的DEM參數(shù)以彌補(bǔ)差異。因此，DEM參數(shù)標(biāo)定是DEM準(zhǔn)確仿真的重要前提，也是大部分DEM仿真研究的重要組成部分。

DEM相關(guān)的文獻(xiàn)出版數(shù)量呈指數(shù)式增長(zhǎng)[16]，而其中只有小部分涉及詳細(xì)的DEM參數(shù)標(biāo)定方法[17]，且有關(guān)詳細(xì)的DEM參數(shù)標(biāo)定文獻(xiàn)大都出自南非Stellenbosch大學(xué)機(jī)電工程系的 Corné Coetzee[1,19-28]。大部分與DEM仿真研究相關(guān)的文獻(xiàn)沒(méi)有提及如何獲得參數(shù)[3-5]或者直接引用其他來(lái)源的參數(shù)[8,10,11,15,29]。

當(dāng)前DEM接觸模型發(fā)展迅速[12,30,31]，但與模型對(duì)應(yīng)的接觸參數(shù)的系統(tǒng)化標(biāo)定方法發(fā)展相對(duì)緩慢，在一定程度上影響了DEM仿真研究的發(fā)展。通過(guò)總結(jié)近10年國(guó)內(nèi)外研究采用的不同標(biāo)定流程以及標(biāo)定方法，從選擇接觸模型、參數(shù)獲取方法、顆粒建模、標(biāo)定實(shí)驗(yàn)與仿真試驗(yàn)等方面，分析和討論目前應(yīng)用于農(nóng)業(yè)領(lǐng)域的DEM參數(shù)標(biāo)定方法，為研究人員未來(lái)改進(jìn)或開(kāi)發(fā)標(biāo)定流程奠定基礎(chǔ)，并為開(kāi)發(fā)更加系統(tǒng)化的且能經(jīng)過(guò)實(shí)際應(yīng)用驗(yàn)證的標(biāo)定方法提供有價(jià)值的參考。

1 選擇接觸模型

選擇合適的接觸模型是DEM仿真的重要前提。DEM接觸模型是一種基于胡克定律的計(jì)算方法，用于描述相互接觸顆粒間的力學(xué)變化規(guī)律。依據(jù)不同的研究對(duì)象或不同的研究目的，適用的接觸模型有所不同，而不同的接觸模型輸入的參數(shù)也有所不同。以EDEM_v2018軟件為例，常用的接觸模型輸入?yún)?shù)及其適用范圍如表1所示。

表1 EDEM中常用接觸模型輸入?yún)?shù)及其適用范圍

2 參數(shù)獲取方法

2.1 直接測(cè)量方法

直接測(cè)量方法是直接測(cè)量接觸水平的參數(shù)，例如靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)。Ucgul等[32]利用如圖1所示的斜面試驗(yàn)測(cè)量沙-鋼的滾動(dòng)摩擦系數(shù)。其中，鋼球是直徑為18.98mm質(zhì)量為28g經(jīng)過(guò)拋光處理的40號(hào)鋼，沙盤(pán)是裝有壓實(shí)平整沙子的托盤(pán)。該試驗(yàn)通過(guò)緩慢增大沙盤(pán)傾斜角，獲得鋼球剛開(kāi)始滾動(dòng)時(shí)的傾斜角ε，再由式（1）計(jì)算沙-鋼的滾動(dòng)摩擦系數(shù)。

圖1 斜面試驗(yàn)示意圖[32]

式中，rμ是沙-鋼滾動(dòng)摩擦系數(shù)；ε是鋼球剛開(kāi)始滾動(dòng)時(shí)沙盤(pán)的傾斜角。

Barrios等[33]通過(guò)如圖2所示的自由落體試驗(yàn)測(cè)量鐵礦石-鐵礦石、鐵礦石-鋼和鐵礦石-橡膠之間的恢復(fù)系數(shù)。使單個(gè)鐵礦石在真空中自由下落，通過(guò)高速攝影分別捕捉碰撞前后5幀的圖像估算碰撞前后速度（v0和v′），最后由v'/ v0計(jì)算得恢復(fù)系數(shù)。

圖2 自由落體試驗(yàn)[33]

鄭侃等[8]利用MXD-01型摩擦系數(shù)測(cè)量?jī)x測(cè)量不同土層顆粒與破土刃之間的靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)，賀一鳴等[34]用比重瓶法測(cè)量土壤的密度，González-Montellano等[35]同樣使用比重瓶法測(cè)量玉米和橄欖的密度。

諸如密度、摩擦系數(shù)、恢復(fù)系數(shù)等參數(shù)可以直接測(cè)量，而要測(cè)量類(lèi)似Edinburgh Elasto-Plastic Adhesion（EEPA）模型定義的范德華型恒定拉脫力可能比較困難，該模型定義的塑性變形比、表面能和切向剛度因子甚至無(wú)法直接測(cè)量。因此，直接測(cè)量法的應(yīng)用有一定局限性。

2.2 批量標(biāo)定方法

批量標(biāo)定方法是利用實(shí)地測(cè)量或?qū)嶒?yàn)室實(shí)驗(yàn)方法測(cè)定特定物料的物理特性，再遵循實(shí)驗(yàn)室或?qū)嵉販y(cè)量的設(shè)置和步驟，并在仿真中模擬實(shí)際情況，然后迭代更改DEM參數(shù)值，直到仿真的整體響應(yīng)與實(shí)際測(cè)量結(jié)果匹配。這個(gè)迭代更改的過(guò)程又可以細(xì)化為基于經(jīng)驗(yàn)法的直接調(diào)整參數(shù)、基于優(yōu)化算法的參數(shù)尋優(yōu)、基于響應(yīng)面法的試驗(yàn)設(shè)計(jì)等。

直接調(diào)整參數(shù)方法一般不涉及參數(shù)獲取過(guò)程，通過(guò)調(diào)整DEM參數(shù)直至仿真的物理特性與實(shí)測(cè)情況匹配。Ucgul等[32]通過(guò)直接調(diào)整仿真時(shí)間步長(zhǎng)、顆粒間靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)，實(shí)現(xiàn)仿真試驗(yàn)中的休止角和能量-貫入深度曲線與實(shí)驗(yàn)室測(cè)量情況基本相同。Janda和Ooi[31]直接調(diào)整EEPA模型參數(shù)，使仿真試驗(yàn)的無(wú)側(cè)限抗壓強(qiáng)度與實(shí)測(cè)的無(wú)側(cè)限抗壓強(qiáng)度呈現(xiàn)良好的一致性。

優(yōu)化算法的參數(shù)獲取方法，通過(guò)揭示研究的物理特性與DEM參數(shù)之間的映射關(guān)系，并通過(guò)這種映射關(guān)系找到最接近實(shí)際情況的參數(shù)組合。苑進(jìn)等[36]通過(guò)參數(shù)敏感性分析，得到對(duì)摻混均勻度影響顯著的參數(shù)是顆粒-傳送帶的靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和顆粒間的滾動(dòng)摩擦系數(shù)。利用相關(guān)向量機(jī)（RVM）揭示這3個(gè)參數(shù)與摻混均勻度的映射關(guān)系，并以平均相對(duì)誤差為適應(yīng)度函數(shù)，通過(guò)遺傳算法（GA）對(duì)初始種群選擇、交叉、變異直至種群平均適應(yīng)度收斂于種群最大適應(yīng)度，從而得到最大適應(yīng)度的參數(shù)組合。Simone等[37]通過(guò)遺傳編程（GP）對(duì)單軸壓縮實(shí)驗(yàn)的敏感性參數(shù)進(jìn)行標(biāo)定。

響應(yīng)面法的試驗(yàn)設(shè)計(jì)方法通過(guò)系統(tǒng)地調(diào)整參數(shù)并重復(fù)試驗(yàn)，擬合試驗(yàn)結(jié)果與參數(shù)之間的多項(xiàng)式函數(shù)，并通過(guò)多項(xiàng)式求解符合實(shí)際情況的最佳參數(shù)組合。賀一鳴等[34]以JKR表面能、顆粒-顆粒靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)為試驗(yàn)因素，休止角為響應(yīng)，按照二次正交旋轉(zhuǎn)組合設(shè)計(jì)3中心點(diǎn)的試驗(yàn)，得到上述參數(shù)與休止角的歸回模型。以實(shí)測(cè)休止角為目標(biāo)得到最終參數(shù)組合，其仿真休止角與實(shí)測(cè)休止角基本一致。李俊偉等[37]以休止角為參考，利用基于響應(yīng)面法的Box-Behnken試驗(yàn)設(shè)計(jì)得到顆粒間JKR表面能、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)對(duì)休止角響應(yīng)的二次回歸模型，獲得類(lèi)似的結(jié)果。王憲良等[29]利用類(lèi)似的方法獲得與實(shí)測(cè)休止角、黏聚力和內(nèi)摩擦角匹配的顆粒半徑、顆粒間靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)。

2.3 方法小結(jié)

對(duì)于直接測(cè)量方法，摩擦系數(shù)、恢復(fù)系數(shù)、剪切模量等參數(shù)比較容易測(cè)量，細(xì)觀參數(shù)的測(cè)量則比較困難，數(shù)值化的參數(shù)甚至無(wú)法直接測(cè)量。再者，顆粒形狀、尺寸和接觸方式的高度理想化決定了即使可以直接準(zhǔn)確地測(cè)量所有參數(shù)，也不一定意味著建立的DEM模型會(huì)顯現(xiàn)相同的整體性能。對(duì)于批量標(biāo)定方法，直接調(diào)整參數(shù)滿足一個(gè)以上條件是比較難實(shí)現(xiàn)的，響應(yīng)面法的試驗(yàn)設(shè)計(jì)次數(shù)會(huì)隨著輸入?yún)?shù)的增加而急劇增加，優(yōu)化算法的計(jì)算成本也會(huì)隨之劇增。因此，通過(guò)簡(jiǎn)單的直接測(cè)量確定一部分輸入?yún)?shù)，再對(duì)剩余輸入?yún)?shù)進(jìn)行批量標(biāo)定是一種節(jié)約研究成本、提高工作效率的有效手段。

3 顆粒建模

在進(jìn)行標(biāo)定實(shí)驗(yàn)和仿真試驗(yàn)之前，除了選擇合適的接觸模型，確定適宜的參數(shù)獲取方法，還應(yīng)謹(jǐn)慎考慮研究對(duì)象的形狀和粒度分布。在大多數(shù)DEM軟件中，設(shè)定離散顆粒的基本單元為球體，因?yàn)榍蛐晤w粒有助于高效檢測(cè)和判斷接觸，且不用考慮方向。研究表明[39]，非球形單元的計(jì)算成本是球形單元的2至3倍，甚至是10倍。本文中討論的是以3D球體為基本單元的顆粒建模，不包括2D圓盤(pán)。

3.1 顆粒形狀

農(nóng)業(yè)領(lǐng)域的DEM研究對(duì)象形狀各異。由多個(gè)球形顆粒組合成團(tuán)塊，是目前DEM研究人員表征各種形狀的常用方法。該團(tuán)塊由兩個(gè)或多個(gè)球形顆粒剛性連接而成，球形顆?？梢栽谌魏纬潭壬现丿B且相互之間不會(huì)產(chǎn)生接觸力，仿真過(guò)程中團(tuán)塊不管受力如何都不會(huì)破裂[40]。大多數(shù)DEM仿真研究不會(huì)提及具體的顆粒形狀建模過(guò)程，少數(shù)研究人員對(duì)比了不同顆粒形狀產(chǎn)生的影響。

Markauskas等[41]分別使用如圖3所示的4、6、8和12個(gè)球形顆粒組成的軸對(duì)稱(chēng)團(tuán)塊模擬玉米料倉(cāng)卸料的過(guò)程，通過(guò)調(diào)整團(tuán)塊顆粒的密度以達(dá)到相同的質(zhì)量，并與實(shí)際測(cè)量結(jié)果進(jìn)行比較。結(jié)果表明，如果這4種團(tuán)塊都使用相同的接觸參數(shù)，那么隨著顆粒數(shù)的增加，卸料速率也會(huì)隨之增加。通過(guò)調(diào)整4種顆粒間的摩擦系數(shù)可以使得每種團(tuán)塊模擬的卸料速率都相同。

圖3 使用4（a），8（b）和12（c）個(gè)球形顆粒組合的玉米籽粒[41]

當(dāng)研究對(duì)象形狀近似球形或橢球時(shí)，為減少仿真時(shí)間和建模難度，將近似球形或橢球的物料簡(jiǎn)化為球形顆粒或球形單元組成的橢球團(tuán)塊是一種常見(jiàn)的辦法。Markauskas和 Kaˇcianauskas[42]使用如圖4所示的不同數(shù)量球形單元組合且長(zhǎng)短軸比分別為1.5、2.35、5和10的橢球團(tuán)塊，用于模擬稻米休止角試驗(yàn)。在進(jìn)行休止角對(duì)顆粒間靜摩擦系數(shù)（0.1-1.0）和滾動(dòng)摩擦系數(shù)（0-0.3）的敏感度分析過(guò)程中發(fā)現(xiàn)：在滾動(dòng)摩擦系數(shù)小于0.3的情況下，即使采用1.0的靜摩擦系數(shù)，不論是何種形狀顆粒的試驗(yàn)休止角都無(wú)法達(dá)到實(shí)測(cè)的35°；當(dāng)滾動(dòng)摩擦系數(shù)設(shè)為0.3時(shí)，采用0.4左右的靜摩擦系數(shù)，每個(gè)形狀顆粒的試驗(yàn)休止角都可以達(dá)到實(shí)測(cè)休止角。

圖4 不同長(zhǎng)短軸比和不同球體數(shù)量的橢球稻米顆粒[42]

研究對(duì)象近似于球體或橢球，其與DEM中理想化的球形橢球形顆粒在內(nèi)部流動(dòng)方面仍有顯著差異。因此，在將顆粒模型簡(jiǎn)化為軸對(duì)稱(chēng)的橢球或球形時(shí)，摩擦系數(shù)的物理意義損失較大。李俊偉等[37]將粘土顆粒簡(jiǎn)化為球形顆粒進(jìn)行休止角的標(biāo)定實(shí)驗(yàn)，在最終的標(biāo)定結(jié)果中靜摩擦系數(shù)達(dá)到0.78和0.80，在賀一鳴等[33]的標(biāo)定結(jié)果中靜摩擦系數(shù)達(dá)到1.06，更加驗(yàn)證了文獻(xiàn)[42]的結(jié)果。

農(nóng)業(yè)物料（例如果實(shí)和種子）的形狀建模最常用的方法是將其模擬成球形顆粒或球形顆粒組成的團(tuán)塊[43]。谷物收獲后流程（例如清選和卸料）的物料建模取決于特定的谷物和實(shí)際工況。針對(duì)諸如稻谷、小麥、玉米等非球形物料，一般采用團(tuán)塊模擬[44]。使用團(tuán)塊模擬物料時(shí)，顆粒形狀建模至關(guān)重要。目前對(duì)顆粒形狀進(jìn)行詳細(xì)標(biāo)定的研究比較少見(jiàn)。大多數(shù)研究人員在使用團(tuán)塊顆粒時(shí)沒(méi)有交待詳細(xì)的建模過(guò)程或者只簡(jiǎn)單對(duì)比了不同形狀物料的特性[27]。像大豆和油菜這種近似于球形的物料，用單個(gè)球形顆粒模擬可能更加準(zhǔn)確。如果使用單個(gè)球形顆粒模擬非球形物料，那么（與用橢球模擬類(lèi)橢球物料類(lèi)似）有必要引入更大的摩擦系數(shù)限制其流動(dòng)性。

3.2 顆粒尺寸

實(shí)驗(yàn)室規(guī)模的顆粒建?？赡軐?duì)粒度分布進(jìn)行較精確的建模。對(duì)于某些大量顆粒建模（例如土槽建模，尤其是粒度分布測(cè)量困難、土壤結(jié)構(gòu)復(fù)雜的農(nóng)田土壤），按照實(shí)際粒度分布進(jìn)行準(zhǔn)確建模可能會(huì)有數(shù)以百萬(wàn)計(jì)乃至千萬(wàn)計(jì)的顆粒。由于一般難以接受如此巨大的計(jì)算成本，因此為降低計(jì)算成本，在仿真中可以增大粒徑。目前國(guó)內(nèi)外流行的可行方案有：保持模型域即幾何體尺寸不變的同時(shí)增加粒徑[32,35,37,43]；忽略尺寸小于某個(gè)特定值的顆粒[46-48]；按比例放大所有顆粒尺寸[47,49,50]。

Ucgul等[32]對(duì)無(wú)粘性、含水率為0.5%的沙灘沙進(jìn)行參數(shù)標(biāo)定。由于沙子顆粒較小，因此在仿真中使用半徑在9.5～10.5mm范圍內(nèi)的顆粒。結(jié)果表明這種放大的顆?？梢灶A(yù)測(cè)耕作部件的受力情況。李俊偉等[37]利用休止角試驗(yàn)對(duì)不同含水率的黑粘土進(jìn)行參數(shù)標(biāo)定時(shí)，使用的顆粒半徑在2～4mm的范圍內(nèi)。類(lèi)似地，這種放大的顆粒再現(xiàn)了實(shí)測(cè)休止角。

Roessler和Katterfeld[51]利用提升法的休止角試驗(yàn)對(duì)加水的濕沙進(jìn)行標(biāo)定實(shí)驗(yàn)。傳統(tǒng)的提升法休止角試驗(yàn)是典型無(wú)粘性物料的休止角測(cè)定方法。然而在對(duì)粘性物料進(jìn)行標(biāo)定實(shí)驗(yàn)時(shí)，發(fā)現(xiàn)提升空心圓柱后的物料并沒(méi)有呈現(xiàn)出良好的錐形，測(cè)量的休止角在40.4°～84.3°之間。因此，在勻速提升圓柱的過(guò)程中，每當(dāng)圓柱底面與落料平面間隙分別為80mm、100mm、120mm和140mm時(shí)，記錄下20mm、40mm、80mm和120mm高度下物料的截面直徑，并用定義的相對(duì)彎曲代替休止角，如圖5所示。進(jìn)一步對(duì)比縮小2.5倍尺寸的空心圓柱和8～100mm/s的提升速率，得到的結(jié)論為：相對(duì)彎曲與提升速度無(wú)關(guān)，在保持圓柱相同高徑比的前提下與其尺寸也無(wú)關(guān)。因而在仿真標(biāo)定試驗(yàn)中，使用相同高徑比但所有尺寸都縮小2.5倍的空心圓柱。此外，由于沙粒太小，半徑在5.5～12.7mm范圍的球形顆粒被用于模擬沙粒。沒(méi)有完全復(fù)制實(shí)際情況的粒度分布，標(biāo)定后的結(jié)果仍然與實(shí)測(cè)值呈現(xiàn)良好的一致性。這與文獻(xiàn)[52]的結(jié)論相似。

圖5 提升過(guò)程中測(cè)量指定高度的截面直徑[51]

Feng等[47]開(kāi)發(fā)了基于三個(gè)相似性原理（幾何體相似性、力學(xué)相似性和動(dòng)力學(xué)相似性）用于離散元顆粒系統(tǒng)的縮放定律。其研究結(jié)果表明：如果顆粒間作用力-重疊量具有以下式（2）所示一般形式，且α與β滿足以下式（3）所示條件，則保持力學(xué)和動(dòng)力學(xué)條件不變并以相同的比例因子縮放幾何體尺寸，仿真中的結(jié)果具有尺度不變性（不論縮放比例如何，結(jié)果都會(huì)保持不變）。

其中：k為接觸剛度，r為顆粒平均半徑，δ為重疊量，α與β均為常數(shù)。

其中：nd是常數(shù)，二維空間中為2，三維空間中為3。

由式（2）和式（3）可見(jiàn)，線性接觸定律就是當(dāng)α和β同時(shí)為1的情況下，僅在二維空間中仿真結(jié)果具有尺度不變性。三維空間中線性接觸定律就是當(dāng)α為2和β為1的情況下，如果要滿足尺度不變性，那么應(yīng)使用與粒度縮放因子成反比的乘數(shù)將接觸剛度k進(jìn)行縮放。此外，根據(jù)該文獻(xiàn)的描述，如果幾何體保持未縮放比例，也就是未保持幾何相似性，將導(dǎo)致取決于縮放因子的建模誤差。

Obermayr等[50]證明了Feng等[47]的理論，應(yīng)用α為1/2和β為3/2的三維非線性接觸定律，將幾何體按照相同的比例因子縮放，進(jìn)行顆粒半徑分別為3mm，30mm和300mm的三軸壓縮試驗(yàn)。試驗(yàn)結(jié)果表現(xiàn)為3組試驗(yàn)的應(yīng)力-應(yīng)變曲線具有良好的一致性，如圖6所示。表明三維非線性接觸定律的模型具有尺度不變性。

圖6 不同粒徑的應(yīng)力-應(yīng)變曲線[50]

4 標(biāo)定實(shí)驗(yàn)與仿真試驗(yàn)

選擇合適的接觸模型，確定適宜的參數(shù)獲取方法，針對(duì)研究對(duì)象建立顆粒模型，還應(yīng)根據(jù)工況進(jìn)行標(biāo)定實(shí)驗(yàn)和仿真試驗(yàn)。國(guó)內(nèi)外研究人員依據(jù)研究工況，利用休止角和土壤-觸土部件相互作用[28,32,34,37,51-53,62-72]、直剪切試驗(yàn)[28,29,53-60]、壓縮試驗(yàn)[61,62]等對(duì)DEM材料參數(shù)（泊松比、顆粒半徑、剪切模量等）和模型接觸參數(shù)的標(biāo)定提出了諸多方法。

4.1 休止角和土壤-觸土部件相互作用

休止角（Angle of repose）是反應(yīng)物料流動(dòng)性的重要指標(biāo)，國(guó)內(nèi)外學(xué)者常用的測(cè)量方法如圖7所示。其中：注入法是指將物料從一定高度的漏斗均勻注入到圓盤(pán)中，直到圓盤(pán)上形成穩(wěn)定的錐形體，測(cè)量錐形體底角作為休止角；排出法是指將物料填充在圓筒中，然后勻速提升圓筒，測(cè)量形成錐形體的底角作為休止角，或者將物料填充至方盒中，通過(guò)撤去方盒的一個(gè)側(cè)面排出物料，測(cè)量物料排出后形成的斜坡角作為休止角；傾斜角法是指將物料填充至托盤(pán)中，緩速傾斜托盤(pán)，把物料處于剛好不下滑狀態(tài)下的傾斜角視為休止角；轉(zhuǎn)動(dòng)圓筒法（Rotating drum）是指將物料置于密封圓筒中，通過(guò)勻速轉(zhuǎn)動(dòng)圓筒，測(cè)量圓筒中形成的傾斜角作為休止角。針對(duì)同種顆粒材料，采用不同測(cè)量方法獲得的休止角也會(huì)存在差異。

圖7 常用休止角測(cè)量方法

賀一鳴等[34]以注入法測(cè)量的休止角為參考，在仿真中遵循實(shí)際實(shí)驗(yàn)設(shè)置和步驟，利用基于響應(yīng)面法二次正交旋轉(zhuǎn)組合試驗(yàn)設(shè)計(jì)得到顆粒間JKR表面能、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)對(duì)休止角響應(yīng)的二次回歸模型，并以實(shí)測(cè)休止角為目標(biāo)對(duì)上述4個(gè)參數(shù)進(jìn)行標(biāo)定。李俊偉等[37]使用類(lèi)似的方法進(jìn)行標(biāo)定實(shí)驗(yàn)。Combarros等[53]利用轉(zhuǎn)動(dòng)圓筒法測(cè)量休止角，在仿真試驗(yàn)中發(fā)現(xiàn)顆粒間滾動(dòng)摩擦系數(shù)和靜摩擦系數(shù)是休止角的敏感性參數(shù)，而單個(gè)實(shí)驗(yàn)參考不足以確定這兩個(gè)參數(shù)。因此將排出法的休止角加入其中，以轉(zhuǎn)動(dòng)滾筒和排出法測(cè)量的休止角為目標(biāo)確定了上述兩參數(shù)。劉凡一等[63]采用排出法測(cè)量小麥的休止角，并對(duì)顆粒間靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和小麥-有機(jī)玻璃靜摩擦系數(shù)進(jìn)行標(biāo)定。

另外，不少學(xué)者在進(jìn)行基于休止角實(shí)驗(yàn)的DEM參數(shù)標(biāo)定中發(fā)現(xiàn)：顆粒間靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)是影響休止角的主要因素[28,51-53,69-72]。這可能為休止角實(shí)驗(yàn)的參數(shù)標(biāo)定節(jié)約參數(shù)敏感性分析的環(huán)節(jié)。

Mak等[66]對(duì)土壤-觸土部件相互作用進(jìn)行研究，確定建模粘性土壤所需的顆粒間法向剛度、切向剛度、摩擦系數(shù)、粘結(jié)法向剛度、粘結(jié)切向剛度、粘結(jié)法向強(qiáng)度、粘結(jié)剪切強(qiáng)度、粘結(jié)鍵半徑系數(shù)、局部阻尼系數(shù)和粘性阻尼系數(shù)。只對(duì)顆粒間的法向剛度進(jìn)行標(biāo)定，其他參數(shù)均取自文獻(xiàn)。在標(biāo)定法向剛度過(guò)程中，沒(méi)有進(jìn)行獨(dú)立的標(biāo)定實(shí)驗(yàn)，而是將與實(shí)際實(shí)驗(yàn)相同外形尺寸和相同行進(jìn)速度的仿真葉片分別作用于5個(gè)由不同法向剛度顆粒組成的土壤模型，通過(guò)對(duì)比葉片水平、垂直方向的阻力和土壤擾動(dòng)寬度與McKye公式[67]預(yù)測(cè)值的誤差，將這3個(gè)誤差平均值最小情況下的顆粒法向剛度作為標(biāo)定值。Li等[15]采用完全相同的方法獲取顆粒間法向剛度，如圖8所示，不同的是將葉片水平、垂直方向的阻力和土壤擾動(dòng)寬度與UEE公式[68]預(yù)測(cè)值進(jìn)行比較。文獻(xiàn)[15,66]均未與實(shí)際實(shí)驗(yàn)進(jìn)行比較。

圖8 葉片入土試驗(yàn)[15]

Ucgul等[32]在標(biāo)定實(shí)驗(yàn)中進(jìn)行如圖9所示的圓盤(pán)和圓錐貫入阻力試驗(yàn)，通過(guò)直接調(diào)整法使得仿真的能量-貫入深度曲線與實(shí)際情況匹配。

圖9 貫入阻力試驗(yàn)與標(biāo)定結(jié)果[32]

4.2 直剪切試驗(yàn)

Karkala等[54]使用來(lái)自文獻(xiàn)[55]的DEM參數(shù)在EDEM軟件中建立一個(gè)無(wú)側(cè)限抗壓強(qiáng)度（動(dòng)態(tài)屈服強(qiáng)度）試驗(yàn)和直剪切試驗(yàn)的對(duì)照試驗(yàn)。之后，又在這個(gè)對(duì)照試驗(yàn)的基礎(chǔ)上使用10倍的剪切速率，并將該試驗(yàn)結(jié)果與對(duì)照試驗(yàn)結(jié)果對(duì)比發(fā)現(xiàn)：加速剪切狀態(tài)下，剪切應(yīng)力確實(shí)更早地達(dá)到臨界值，而臨界剪切應(yīng)力并未有顯著變化。該結(jié)果表明提高剪切速率是縮短仿真剪切試驗(yàn)時(shí)間的有效辦法。在仿真的無(wú)側(cè)限抗壓強(qiáng)度試驗(yàn)中，壓縮后的試樣（直徑25mm、高25mm的圓柱）被去除限制壁，頂板以100mm/s的恒定速率（應(yīng)變率=4s-1）下壓直至試樣失效。期間記錄頂板的壓力和試樣端面的變形，并將壓力換算為壓應(yīng)力。分析發(fā)現(xiàn)JKR表面可能是無(wú)側(cè)限抗壓強(qiáng)度的極敏感參數(shù)。通過(guò)“試錯(cuò)法”（在參數(shù)的可行域內(nèi)等分若干組并逐一嘗試），選取一組與實(shí)測(cè)值誤差最小的JKR表面能參數(shù)。在這項(xiàng)研究中并未與實(shí)際實(shí)驗(yàn)進(jìn)行比較，4s-1的應(yīng)變率是否與實(shí)際情況差異較大，也沒(méi)有像使用10倍剪切速率那樣驗(yàn)證應(yīng)變率是否會(huì)對(duì)最終的峰值壓應(yīng)力有顯著影響。

王憲良等[29]結(jié)合直剪切試驗(yàn)和排出法的休止角試驗(yàn)，分析休止角、黏聚力和內(nèi)摩擦角對(duì)顆粒密度、半徑、泊松比以及顆粒間表面能、法向剛度、切向剛度、剪切模量、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)的敏感度。結(jié)果表明休止角、黏聚力和內(nèi)摩擦角對(duì)顆粒半徑、顆粒間靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)的敏感度都比其他參數(shù)高出一個(gè)數(shù)量級(jí)。因此，利用響應(yīng)面法可得到這3項(xiàng)參考與3個(gè)主要參數(shù)的回歸模型，優(yōu)化求解可得到唯一的參數(shù)組合，即顆粒半徑、顆粒間的靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)分別為5.7mm、0.45和0.21。在試驗(yàn)中，并未對(duì)該特定模型定義的接觸參數(shù)進(jìn)行標(biāo)定，所有接觸模型參數(shù)均取自其他文獻(xiàn)；利用輪胎-土壤接觸面應(yīng)力驗(yàn)證了參數(shù)的準(zhǔn)確性。在驗(yàn)證過(guò)程中，沒(méi)有將該參數(shù)組合的仿真值與實(shí)測(cè)值進(jìn)行對(duì)比，而是將模型預(yù)測(cè)值與實(shí)測(cè)值進(jìn)行比較。根據(jù)文獻(xiàn)[27]，離散元參數(shù)標(biāo)定的最終目的是使得標(biāo)定的參數(shù)組合能應(yīng)用于實(shí)際，而不是為了與設(shè)計(jì)的標(biāo)定試驗(yàn)匹配。

4.3 壓縮試驗(yàn)

Li等[61]利用一系列的二維雙軸DEM模擬，通過(guò)改變接觸法向剛度、接觸切向剛度和接觸滑動(dòng)摩擦系數(shù)來(lái)找到響應(yīng)面系數(shù)。對(duì)于給定的圍壓和載荷步，求解7個(gè)方程，可以找到7個(gè)響應(yīng)面系數(shù)。在給定的圍壓和10個(gè)加載步驟中的每一個(gè)加載下，法向剛度、切向剛度和摩擦系數(shù)與總體偏應(yīng)力相關(guān)。使用優(yōu)化算法，可以最小化被測(cè)應(yīng)力和預(yù)測(cè)偏應(yīng)力的均方根值，以找到最適合給定圍壓的測(cè)量數(shù)據(jù)。結(jié)果表明，法向剛度和切向剛度大致相同，且在所有情況下的差異均小于5%；摩擦系數(shù)在0.92至1.04之間變化。如果要考慮顆粒形狀，則必須增加滑動(dòng)摩擦系數(shù)，并且該摩擦系數(shù)很可能會(huì)高于兩個(gè)物理顆粒之間的測(cè)量值[62]。結(jié)果還表明接觸剛度和摩擦系數(shù)隨著圍壓的增加而增加，證明參數(shù)和試驗(yàn)參考值取決于應(yīng)力，應(yīng)仔細(xì)選擇在標(biāo)定實(shí)驗(yàn)中使用的應(yīng)力水平。

4.4 其他試驗(yàn)

Roessler等[17]針對(duì)同一種干粗礫石進(jìn)行落料箱試驗(yàn)，如圖10所示。目的是為了對(duì)顆粒間的靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)進(jìn)行標(biāo)定。每種試驗(yàn)重復(fù)三次并將測(cè)量誤差考慮在內(nèi)，這樣對(duì)應(yīng)每個(gè)參考的參數(shù)可行解在等高線圖上就不是單條曲線而是誤差范圍內(nèi)的可行域。將排出法和注入法休止角的可行域疊加，發(fā)現(xiàn)重疊部分比任何單獨(dú)參考值的可行域都要小很多，得出使用單個(gè)參考標(biāo)準(zhǔn)的標(biāo)定試驗(yàn)不足以確定可行解的結(jié)論，與文獻(xiàn)[53]的結(jié)論相似。此外，還將能匹配落料箱試驗(yàn)中下部箱體排出質(zhì)量、落料孔平均質(zhì)量流率和排出法的休止角參數(shù)可行域疊加，得到比之前更小的可行域。因此，單獨(dú)的落料箱試驗(yàn)可獲取多個(gè)參考，且在考慮測(cè)量誤差的前提下能有效減少可行解的數(shù)量。

圖10 落料箱試驗(yàn)

李俊偉等[37]利用斜面試驗(yàn)，即將3D打印的半徑為5mm的土粒分別靜置在相同傾斜角的65Mn和PTEF材料表面，使其從同一高度自然滾下，測(cè)量土粒的水平滾動(dòng)距離，對(duì)觸土部件與土壤間的JKR表面能、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)進(jìn)行標(biāo)定。在該研究中，即使是3D打印的球形土粒也存在一定的孔隙，也就是說(shuō)這個(gè)土粒的堆積密度固然小于土壤的真實(shí)密度，而用于斜面試驗(yàn)仿真的顆粒密度是如何定義的無(wú)從知曉。土壤與觸土部件的相互作用，不應(yīng)是單個(gè)土粒與部件的相互作用，而應(yīng)是觸土部件作用于土壤顆粒群時(shí)，部件對(duì)顆粒群整體的宏觀響應(yīng)或某些力學(xué)方面的變化，或者是土壤顆粒群作用于觸土部件時(shí)，顆粒群對(duì)部件的作用效果。采用離散元法，本質(zhì)上是利用各個(gè)離散元素表征研究對(duì)象的整體運(yùn)動(dòng)規(guī)律。這種將物料顆粒模型簡(jiǎn)化的方法一般是可取的[74,75]，而為簡(jiǎn)化模型將實(shí)際土壤制作成球形顆粒略顯得不合時(shí)宜。3D打印球體的球面度也可能成為影響滾動(dòng)距離的主要因素，不同打印精度的土球滾動(dòng)距離差異可能較大。斜面試驗(yàn)適用于直接測(cè)量顆粒與材料的靜摩擦系數(shù)或滾動(dòng)摩擦系數(shù)[28,62-65]。

5 發(fā)展與建議

獲取DEM參數(shù)是DEM仿真的重要內(nèi)容。雖然有關(guān)參數(shù)標(biāo)定的研究發(fā)展迅速，但仍存在一些亟待解決的問(wèn)題。為推動(dòng)參數(shù)標(biāo)定方法的改進(jìn)與創(chuàng)新，提出以下建議：

（1）DEM待輸入?yún)?shù)眾多，為簡(jiǎn)化參數(shù)獲取過(guò)程，一般應(yīng)將批量標(biāo)定方法與直接測(cè)量法配合使用，例如標(biāo)定之前直接測(cè)量材料的密度、泊松比、剪切模量、摩擦系數(shù)等。農(nóng)業(yè)物料外形各異、尺寸不一，農(nóng)田土壤粒度分布復(fù)雜、孔隙度不均勻。DEM模擬的物料外形尺寸或土壤粒度和孔隙度存在一定的誤差。由于該誤差會(huì)導(dǎo)致仿真對(duì)象的物理特性與實(shí)際情況存在較大差異，因此不建議完全使用直接測(cè)量法。

（2）農(nóng)業(yè)物料形狀各異，將多個(gè)球形顆粒組合成團(tuán)塊是目前DEM研究人員表征物料形狀的常用方法。使用這種團(tuán)塊的優(yōu)點(diǎn)在于其表面仍是多個(gè)球面的組合，接觸檢測(cè)和判定仍然與球形顆粒相同。也正是因?yàn)檫@種球面的組合決定了要精確建模物料的尖銳邊緣和光滑表面就需要大量的球形顆粒，這就大幅度增加工作量和仿真時(shí)間。因此，有必要將顆粒形狀建模作為一種DEM參數(shù)。大多數(shù)研究中不會(huì)提及詳細(xì)的形狀建模過(guò)程，系統(tǒng)化的形狀建模方法仍然匱乏。在顆粒碰撞、內(nèi)部流動(dòng)等與顆粒形狀密切相關(guān)的研究中，顆粒形狀建模十分重要。一旦對(duì)某種顆粒形狀的DEM參數(shù)進(jìn)行標(biāo)定，該接觸參數(shù)可能只適合該形狀的顆粒。因此，在調(diào)整建模形狀時(shí)，也應(yīng)對(duì)DEM參數(shù)進(jìn)行調(diào)整或重新標(biāo)定。

（3）在不考慮顆粒形狀影響的情況下，簡(jiǎn)化顆粒模型是一種常用手段。例如使用球形顆粒組成的橢球模型模擬稻米和小麥，使用球形顆粒模擬土壤顆粒、油菜籽粒、大豆籽粒等。研究對(duì)象即使非常近似橢球或球形，也與DEM中理想的橢球或球形在內(nèi)部流動(dòng)性方面與實(shí)際情況有較大差異。由于使用較大的摩擦系數(shù)可限制這種流動(dòng)性，因此在參數(shù)標(biāo)定之前，建議選擇較大的摩擦系數(shù)范圍。

（4）只對(duì)敏感參數(shù)進(jìn)行標(biāo)定，是國(guó)內(nèi)外學(xué)者常用的方法。在這種情況下，對(duì)非敏感參數(shù)的確定存在一定的不確定性。目前，除JKR等定義參數(shù)較少的接觸模型以外，模型全部接觸參數(shù)的標(biāo)定方法鮮有報(bào)道。例如，EDEM軟件內(nèi)置的EEPA接觸模型定義的參數(shù)加上滾動(dòng)摩擦系數(shù)、靜摩擦系數(shù)和恢復(fù)系數(shù)，共有9個(gè)建立模型所需的接觸參數(shù)。由于目前可接受的待標(biāo)定參數(shù)至多4個(gè)，因此對(duì)多參數(shù)模型的全參數(shù)進(jìn)行標(biāo)定具有一定挑戰(zhàn)性。

（5）不論參數(shù)標(biāo)定流程如何，在標(biāo)定結(jié)束時(shí)都需要提供一組參數(shù)組合。仿真試驗(yàn)的整體響應(yīng)（如休止角）往往受一個(gè)以上參數(shù)的影響，這就意味著回歸模型至少是二元的。在大多數(shù)標(biāo)定過(guò)程中，僅考慮單個(gè)響應(yīng)獲取的參數(shù)組合一般不是唯一的。在選擇最終的一組參數(shù)組合時(shí)，這種模糊性將不利于DEM參數(shù)標(biāo)定程序的研發(fā)。為減少可行參數(shù)組合甚至獲得唯一的參數(shù)組合，建議在單個(gè)實(shí)驗(yàn)中獲取多個(gè)與研究相關(guān)的參考值，或進(jìn)行多個(gè)標(biāo)定實(shí)驗(yàn)。

（6）標(biāo)定實(shí)驗(yàn)和仿真試驗(yàn)的力學(xué)和動(dòng)力學(xué)條件應(yīng)該保持一致，在不得不改變?cè)摋l件之前，應(yīng)分析變化的力學(xué)和動(dòng)力學(xué)條件對(duì)試驗(yàn)結(jié)果的影響。

（7）參數(shù)標(biāo)定的最終目的不僅是為了使標(biāo)定的參數(shù)能匹配實(shí)驗(yàn)室結(jié)果，還要與實(shí)驗(yàn)相關(guān)的實(shí)際應(yīng)用結(jié)果相對(duì)應(yīng)。例如，匹配實(shí)驗(yàn)室測(cè)量休止角的DEM參數(shù)，還應(yīng)匹配開(kāi)溝實(shí)驗(yàn)的溝型。因此，為驗(yàn)證標(biāo)定方法的實(shí)用性，建議除了進(jìn)行實(shí)驗(yàn)室規(guī)模的驗(yàn)證外，還應(yīng)進(jìn)行實(shí)際應(yīng)用的驗(yàn)證。

6 結(jié)論

農(nóng)業(yè)物料外形各異、尺寸不一，農(nóng)田土壤粒度分布復(fù)雜、孔隙度不均勻，DEM模擬的物料外形尺寸或土壤粒度和孔隙度存在一定的誤差。因此，為在DEM中再現(xiàn)研究對(duì)象的實(shí)際物理特性，需要確定合理的DEM參數(shù)以彌補(bǔ)誤差。分析、總結(jié)了目前常用的直接測(cè)量方法和批量標(biāo)定方法。直接測(cè)量法可用于直接測(cè)定顆粒密度、剪切模量、接觸剛度、滑動(dòng)摩擦系數(shù)、恢復(fù)系數(shù)等，其優(yōu)點(diǎn)為測(cè)定的參數(shù)在不同的軟件下或不同的接觸模型中均可受用，其缺點(diǎn)為在顆粒尺寸較小或形狀復(fù)雜的情況下，測(cè)量工作將很難展開(kāi)。在顆粒形狀、尺寸和接觸方式高度理想化的情況下，精確測(cè)量參數(shù)值也很可能導(dǎo)致整體性能與實(shí)際情況存在較大差異。批量標(biāo)定方法的優(yōu)點(diǎn)在于不需要精密的儀器測(cè)量顆粒級(jí)或接觸水平的參數(shù)，只需測(cè)量一項(xiàng)或多項(xiàng)整體性能指標(biāo)。缺點(diǎn)就是重復(fù)地進(jìn)行仿真試驗(yàn)非常耗時(shí)，有時(shí)一次試驗(yàn)獲得參考數(shù)量不足以獲取更少甚至唯一的參數(shù)組合。為簡(jiǎn)化標(biāo)定流程，一般應(yīng)將批量標(biāo)定方法與直接測(cè)量方法結(jié)合使用。顆粒形狀和粒度分布也應(yīng)被視為DEM參數(shù)。在研究顆粒碰撞、流動(dòng)性等情況下，不應(yīng)忽略形狀和粒度的建模過(guò)程。標(biāo)定實(shí)驗(yàn)結(jié)果應(yīng)與研究目的對(duì)應(yīng)。此外，仿真試驗(yàn)應(yīng)遵循標(biāo)定實(shí)驗(yàn)的設(shè)置和步驟。

從接觸模型選擇、參數(shù)獲取方法、顆粒建模、標(biāo)定實(shí)驗(yàn)和仿真試驗(yàn)方面，陳述和討論了農(nóng)業(yè)領(lǐng)域DEM參數(shù)標(biāo)定現(xiàn)狀，并提出了幾點(diǎn)建議，有助于未來(lái)研究人員改進(jìn)或開(kāi)發(fā)標(biāo)定流程，為更加系統(tǒng)化且經(jīng)過(guò)實(shí)際應(yīng)用驗(yàn)證的標(biāo)定方法開(kāi)發(fā)提供有價(jià)值的參考。