農業領域離散元法（DEM）參數標定現狀與發展

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

（1.湖南農業大學機電工程學院，長沙 410128；2.智能農機裝備湖南省重點實驗室，長沙 410128）

引言

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

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

DEM相關的文獻出版數量呈指數式增長[16]，而其中只有小部分涉及詳細的DEM參數標定方法[17]，且有關詳細的DEM參數標定文獻大都出自南非Stellenbosch大學機電工程系的 Corné Coetzee[1,19-28]。大部分與DEM仿真研究相關的文獻沒有提及如何獲得參數[3-5]或者直接引用其他來源的參數[8,10,11,15,29]。

當前DEM接觸模型發展迅速[12,30,31]，但與模型對應的接觸參數的系統化標定方法發展相對緩慢，在一定程度上影響了DEM仿真研究的發展。通過總結近10年國內外研究采用的不同標定流程以及標定方法，從選擇接觸模型、參數獲取方法、顆粒建模、標定實驗與仿真試驗等方面，分析和討論目前應用于農業領域的DEM參數標定方法，為研究人員未來改進或開發標定流程奠定基礎，并為開發更加系統化的且能經過實際應用驗證的標定方法提供有價值的參考。

1 選擇接觸模型

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

表1 EDEM中常用接觸模型輸入參數及其適用范圍

2 參數獲取方法

2.1 直接測量方法

直接測量方法是直接測量接觸水平的參數，例如靜摩擦系數和滾動摩擦系數。Ucgul等[32]利用如圖1所示的斜面試驗測量沙-鋼的滾動摩擦系數。其中，鋼球是直徑為18.98mm質量為28g經過拋光處理的40號鋼，沙盤是裝有壓實平整沙子的托盤。該試驗通過緩慢增大沙盤傾斜角，獲得鋼球剛開始滾動時的傾斜角ε，再由式（1）計算沙-鋼的滾動摩擦系數。

圖1 斜面試驗示意圖[32]

式中，rμ是沙-鋼滾動摩擦系數；ε是鋼球剛開始滾動時沙盤的傾斜角。

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

圖2 自由落體試驗[33]

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

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

2.2 批量標定方法

批量標定方法是利用實地測量或實驗室實驗方法測定特定物料的物理特性，再遵循實驗室或實地測量的設置和步驟，并在仿真中模擬實際情況，然后迭代更改DEM參數值，直到仿真的整體響應與實際測量結果匹配。這個迭代更改的過程又可以細化為基于經驗法的直接調整參數、基于優化算法的參數尋優、基于響應面法的試驗設計等。

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

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

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

2.3 方法小結

對于直接測量方法，摩擦系數、恢復系數、剪切模量等參數比較容易測量，細觀參數的測量則比較困難，數值化的參數甚至無法直接測量。再者，顆粒形狀、尺寸和接觸方式的高度理想化決定了即使可以直接準確地測量所有參數，也不一定意味著建立的DEM模型會顯現相同的整體性能。對于批量標定方法，直接調整參數滿足一個以上條件是比較難實現的，響應面法的試驗設計次數會隨著輸入參數的增加而急劇增加，優化算法的計算成本也會隨之劇增。因此，通過簡單的直接測量確定一部分輸入參數，再對剩余輸入參數進行批量標定是一種節約研究成本、提高工作效率的有效手段。

3 顆粒建模

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

3.1 顆粒形狀

農業領域的DEM研究對象形狀各異。由多個球形顆粒組合成團塊，是目前DEM研究人員表征各種形狀的常用方法。該團塊由兩個或多個球形顆粒剛性連接而成，球形顆粒可以在任何程度上重疊且相互之間不會產生接觸力，仿真過程中團塊不管受力如何都不會破裂[40]。大多數DEM仿真研究不會提及具體的顆粒形狀建模過程，少數研究人員對比了不同顆粒形狀產生的影響。

Markauskas等[41]分別使用如圖3所示的4、6、8和12個球形顆粒組成的軸對稱團塊模擬玉米料倉卸料的過程，通過調整團塊顆粒的密度以達到相同的質量，并與實際測量結果進行比較。結果表明，如果這4種團塊都使用相同的接觸參數，那么隨著顆粒數的增加，卸料速率也會隨之增加。通過調整4種顆粒間的摩擦系數可以使得每種團塊模擬的卸料速率都相同。

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

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

圖4 不同長短軸比和不同球體數量的橢球稻米顆粒[42]

研究對象近似于球體或橢球，其與DEM中理想化的球形橢球形顆粒在內部流動方面仍有顯著差異。因此，在將顆粒模型簡化為軸對稱的橢球或球形時，摩擦系數的物理意義損失較大。李俊偉等[37]將粘土顆粒簡化為球形顆粒進行休止角的標定實驗，在最終的標定結果中靜摩擦系數達到0.78和0.80，在賀一鳴等[33]的標定結果中靜摩擦系數達到1.06，更加驗證了文獻[42]的結果。

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

3.2 顆粒尺寸

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

Ucgul等[32]對無粘性、含水率為0.5%的沙灘沙進行參數標定。由于沙子顆粒較小，因此在仿真中使用半徑在9.5～10.5mm范圍內的顆粒。結果表明這種放大的顆粒可以預測耕作部件的受力情況。李俊偉等[37]利用休止角試驗對不同含水率的黑粘土進行參數標定時，使用的顆粒半徑在2～4mm的范圍內。類似地，這種放大的顆粒再現了實測休止角。

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

圖5 提升過程中測量指定高度的截面直徑[51]

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

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

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

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

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

圖6 不同粒徑的應力-應變曲線[50]

4 標定實驗與仿真試驗

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

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

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

圖7 常用休止角測量方法

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

另外，不少學者在進行基于休止角實驗的DEM參數標定中發現：顆粒間靜摩擦系數和滾動摩擦系數是影響休止角的主要因素[28,51-53,69-72]。這可能為休止角實驗的參數標定節約參數敏感性分析的環節。

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

圖8 葉片入土試驗[15]

Ucgul等[32]在標定實驗中進行如圖9所示的圓盤和圓錐貫入阻力試驗，通過直接調整法使得仿真的能量-貫入深度曲線與實際情況匹配。

圖9 貫入阻力試驗與標定結果[32]

4.2 直剪切試驗

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

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

4.3 壓縮試驗

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

4.4 其他試驗

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

圖10 落料箱試驗

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

5 發展與建議

獲取DEM參數是DEM仿真的重要內容。雖然有關參數標定的研究發展迅速，但仍存在一些亟待解決的問題。為推動參數標定方法的改進與創新，提出以下建議：

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

（2）農業物料形狀各異，將多個球形顆粒組合成團塊是目前DEM研究人員表征物料形狀的常用方法。使用這種團塊的優點在于其表面仍是多個球面的組合，接觸檢測和判定仍然與球形顆粒相同。也正是因為這種球面的組合決定了要精確建模物料的尖銳邊緣和光滑表面就需要大量的球形顆粒，這就大幅度增加工作量和仿真時間。因此，有必要將顆粒形狀建模作為一種DEM參數。大多數研究中不會提及詳細的形狀建模過程，系統化的形狀建模方法仍然匱乏。在顆粒碰撞、內部流動等與顆粒形狀密切相關的研究中，顆粒形狀建模十分重要。一旦對某種顆粒形狀的DEM參數進行標定，該接觸參數可能只適合該形狀的顆粒。因此，在調整建模形狀時，也應對DEM參數進行調整或重新標定。

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

（4）只對敏感參數進行標定，是國內外學者常用的方法。在這種情況下，對非敏感參數的確定存在一定的不確定性。目前，除JKR等定義參數較少的接觸模型以外，模型全部接觸參數的標定方法鮮有報道。例如，EDEM軟件內置的EEPA接觸模型定義的參數加上滾動摩擦系數、靜摩擦系數和恢復系數，共有9個建立模型所需的接觸參數。由于目前可接受的待標定參數至多4個，因此對多參數模型的全參數進行標定具有一定挑戰性。

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

（6）標定實驗和仿真試驗的力學和動力學條件應該保持一致，在不得不改變該條件之前，應分析變化的力學和動力學條件對試驗結果的影響。

（7）參數標定的最終目的不僅是為了使標定的參數能匹配實驗室結果，還要與實驗相關的實際應用結果相對應。例如，匹配實驗室測量休止角的DEM參數，還應匹配開溝實驗的溝型。因此，為驗證標定方法的實用性，建議除了進行實驗室規模的驗證外，還應進行實際應用的驗證。

6 結論

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

從接觸模型選擇、參數獲取方法、顆粒建模、標定實驗和仿真試驗方面，陳述和討論了農業領域DEM參數標定現狀，并提出了幾點建議，有助于未來研究人員改進或開發標定流程，為更加系統化且經過實際應用驗證的標定方法開發提供有價值的參考。