劉文政 何 進 李洪文 李學強 鄭 侃 魏忠彩

(1.中國農(nóng)業(yè)大學工學院， 北京 100083； 2.山東省馬鈴薯生產(chǎn)裝備智能化工程技術(shù)研究中心， 德州 253600；3.山東希成農(nóng)業(yè)機械科技有限公司， 德州 253600)

0 引言

馬鈴薯作為全球四大主糧作物之一，對于保障糧食安全具有重要意義[1]。種植是馬鈴薯生產(chǎn)的重要一環(huán)，其產(chǎn)量和質(zhì)量與種薯密切相關(guān)，脫毒微型馬鈴薯(簡稱“微型薯”)是由馬鈴薯脫毒組織經(jīng)煉苗、移栽、保護地隔離種植等培育而成，其種出的后代具有產(chǎn)量高、品質(zhì)好、商品率高等優(yōu)點[2]。為提升國內(nèi)馬鈴薯的生產(chǎn)水平和經(jīng)濟效益，應大幅提升微型薯的種植面積和水平，其種植機械化勢在必行[3]。隨著計算機技術(shù)的快速發(fā)展，以離散單元法為基礎(chǔ)的數(shù)值模擬仿真軟件EDEM在農(nóng)業(yè)裝備研究上已得到運用[4-8]。本文以微型薯為研究對象，對其離散元仿真參數(shù)進行全面系統(tǒng)的研究，以促進離散元法在微型馬鈴薯播種裝備研發(fā)中的應用。

離散元仿真時，需定義仿真模型的物性參數(shù)，主要包括本征參數(shù)(如密度、剪切模量、泊松比等)以及接觸參數(shù)(如顆粒與顆粒、顆粒與材料間的彈性碰撞恢復系數(shù)、靜摩擦因數(shù)、滾動摩擦因數(shù)等)，其中物料本征參數(shù)與真實值基本一致，而由于顆粒仿真模型與真實顆粒在幾何形態(tài)上存在差異，使得仿真接觸參數(shù)與真實值存在誤差，需對接觸參數(shù)進行重新標定。目前這些參數(shù)主要通過直接測量和虛擬標定獲取[9-10]。對離散元仿真參數(shù)的直接測量，國內(nèi)外學者已做了大量研究[11-12]。然而，由于真實顆粒的各向異性，使得顆粒仿真模型與真實顆粒在外形上存在差異，直接測量數(shù)據(jù)變化差異較大，同時有些參數(shù)很難通過試驗直接測量，因此，一些學者對所需的微觀參數(shù)進行虛擬標定[13-15]。

目前，國內(nèi)外種子顆粒標定主要集中在玉米[16-18]、水稻[19]、小麥[10]、大豆[20]等方面，而對微型薯種子顆粒進行離散元仿真參數(shù)標定方面的研究鮮有報道。微型薯相比于玉米、稻麥等在質(zhì)量和形態(tài)方面均具有較大差異，其種子顆粒尺寸較大，需進行大量相關(guān)試驗來探索較為復雜的運動特征。本文在現(xiàn)有相關(guān)研究的基礎(chǔ)上，提出一種微型薯等大顆粒種子參數(shù)標定方法。

1 微型薯模型建立

1.1 微型薯物理模型創(chuàng)建

進行仿真試驗時，應建立微型薯種子顆粒仿真物理模型。微型薯作為農(nóng)業(yè)散體物料，其籽粒間差異性大、外形輪廓復雜，為提高離散元顆粒建模準確性，采用在EDEM軟件中利用球形顆粒組合的方法建立微型薯種子仿真模型[21-22]。

微型薯選用國家馬鈴薯工程技術(shù)研究中心培育的質(zhì)量分布在3～10 g的希森3號脫毒微型馬鈴薯，該品種質(zhì)量好，產(chǎn)量高，具有微型馬鈴薯的典型特征。試驗測定其含水率為66.97%～78.23%，平均值為73.90%；密度為1 049.60～1 085.90 kg/m3，平均值為1 077.03 kg/m3；微型薯泊松比為0.48[23]，剪切模量為34.30 MPa[24]。因微型薯種子物料形態(tài)差異很大，為確定微型薯種子物理模型，隨機選取500粒微型薯種子，通過游標卡尺(精度0.02 mm，量程0～150 mm)對顆粒三維尺寸(長l、寬w、厚t)進行測量，依據(jù)測量結(jié)果發(fā)現(xiàn)微型薯種子整體呈橢球形，根據(jù)其形態(tài)特征，引入球度Sp對微型薯顆粒進行種子形狀分類，分為橢球形(Sp<0.85)和類球形(Sp>0.85)2級。其中，球度Sp為[11]

(1)

圖1 微型薯分類及離散元模型Fig.1 Classification and simulation model of potato minituber

為進一步確定微型薯種子物理模型，依據(jù)測量結(jié)果對微型薯橢球形種子按照其長軸尺寸是否大于35 mm分為2級，即微型薯種子顆粒共分為橢球小粒、橢球大粒和類球形(圖1a)，且數(shù)量比為20∶4∶1。為便于在EDEM中建立微型薯種子顆粒模型，對于橢球形種子而言，設(shè)定其長軸長度a=l，短軸長度b=(w+t)/2；對于類球形種子而言，其直徑d=(l+w+t)/3。由此得橢球小粒種子平均長軸長度為29.20 mm，平均短軸長度為18.24 mm；橢球大粒種子平均長軸長度為39.32 mm，平均短軸長度為20.54 mm；類球形種子平均直徑為20.46 mm。根據(jù)上述測量結(jié)果，在EDEM的原型顆粒模型創(chuàng)建中利用多球面組合功能建立橢球大粒、橢球小粒及類球形微型薯種子顆粒離散元模型，如圖1b所示。

1.2 微型薯接觸模型選取

試驗過程中，微型薯除了顆粒間的接觸，還會與其他材料接觸產(chǎn)生力的作用。本文試驗接觸材料選用農(nóng)業(yè)裝備中常用的Q235鋼制材料，其密度為7 850 kg/m3，泊松比為0.28，剪切模量為8.20×104MPa[25]。

在仿真過程中，微型薯顆粒模型作為離散體具有真實微型薯種子顆粒的幾何和物理兩類基本特征，采用離散元方法模擬微型薯仿真顆粒的運動過程。在開展相關(guān)仿真試驗時，依據(jù)微型薯物料特性，其表面粘附力較小，且假設(shè)顆粒在運動過程中位移、力、速度等參數(shù)的變化是通過顆粒間或顆粒與接觸材料之間產(chǎn)生的微小交疊量的不同來確定的，根據(jù)牛頓第二定律，每個微型薯顆粒模型在力和扭矩的作用下發(fā)生運動和旋轉(zhuǎn)?；谝陨霞僭O(shè)，微型薯顆粒的相關(guān)運動過程采用Hertz- Mindlin無滑動接觸力學模型[26-29]。

基于該模型，微型薯顆粒在顆粒間或顆粒- 材料間所受法向力Fn和切向力Ft滿足函數(shù)關(guān)系式

(2)

式中E*——等效彈性模量

R*——等效接觸半徑

δn——法向重疊量

δt——切向重疊量

St——切向剛度

微型薯顆粒間或顆粒- 材料間法向阻尼力Fn，d及切向阻尼力Ft，d滿足

(3)

式中β——阻尼比m*——等效質(zhì)量

Sn——法向剛度

υn,rel——法向相對速度

υt,rel——切向相對速度

微型薯顆粒間或顆粒- 材料間的切向力Ft還受庫倫摩擦力μsFn限制，其中μs為靜摩擦因數(shù)。微型薯仿真顆粒在運動中不可避免受到滾動摩擦力的影響，而滾動摩擦可以通過接觸表面上的力矩Ti來說明，其中，力矩Ti與滾動摩擦因數(shù)μr、顆粒質(zhì)心至接觸點間距Ri、顆粒在接觸點處的角速度ωi存在關(guān)系

Ti=-μrFnRiωi

(4)

式中μr為無量綱參數(shù)，該參數(shù)值的不同會影響顆粒阻抗?jié)L動程度的大小。

2 微型薯離散元參數(shù)獲取模型

2.1 微型薯碰撞恢復系數(shù)測定模型

碰撞恢復系數(shù)是表征物體碰撞時變形恢復能力的參數(shù)，為碰撞前后的兩物體在接觸點處的法向相對分離速度與法向相對接近速度之比[30]。根據(jù)碰撞恢復系數(shù)物理學定義建立其測定模型，分別對微型薯- 鋼板、顆粒間碰撞恢復系數(shù)進行分析和標定。

2.1.1微型薯- 鋼板碰撞恢復系數(shù)

利用微型薯顆粒的碰撞彈跳試驗對微型薯- 鋼板碰撞恢復系數(shù)進行測定[26]，整個運動過程由高速攝像機(美國Vision Research公司，Phantom v 9.1)采集記錄，如圖2a所示。微型薯- 鋼板碰撞恢復系數(shù)可表示為顆粒與鋼板碰撞前后在碰撞接觸點處的法向方向瞬時分離速度v1與瞬時接觸速度v0之比，計算公式為

(5)

式中H0——微型薯顆粒下落高度

H1——微型薯顆粒與鋼板碰撞反彈最大高度

圖2 微型薯- 鋼板碰撞恢復系數(shù)測定試驗Fig.2 Measurement test of potato minituber- steel plate coefficient of restitution

試驗選取類球形微型薯種子顆粒從一定高度H0釋放，使其跌落至水平放置的鋼板上，用高速攝像機采集到種子顆粒彈起的最高點，通過坐標紙讀取該點的高度H1(H1

試驗中選取已建立的類球形微型薯顆粒仿真模型，顆粒以初速度零且距離水平放置的鋼板450 mm處做自由落體運動。在微型薯- 鋼板仿真接觸參數(shù)中，靜摩擦因數(shù)、滾動摩擦因數(shù)對顆粒反彈高度無影響，為避免干擾將這兩個參數(shù)設(shè)為零；參照真實試驗測定的微型薯- 鋼板碰撞恢復系數(shù)取值范圍，設(shè)定微型薯- 鋼板仿真碰撞恢復系數(shù)在0.485～0.570范圍內(nèi)取值，并通過EDEM離散元分析模塊對相應顆粒反彈高度進行測定(圖2b)，經(jīng)過預仿真試驗確定仿真試驗設(shè)計方案與結(jié)果如表1所示。以恢復系數(shù)x1為試驗因素，反彈高度y1為評價指標，對結(jié)果進行曲線擬合，建立擬合方程

(6)

表1 微型薯- 鋼板碰撞恢復系數(shù)仿真試驗方案與結(jié)果Tab.1 Simulation test plan and results of potatominituber- steel plate coefficient of restitution

決定系數(shù)R2=0.999 9，該值接近1，表明擬合方程可靠度高。將y1=124.80 mm代入方程(6)得x1=0.523，設(shè)定微型薯- 鋼板恢復系數(shù)為0.523并分別在H0為400、450、500 mm仿真試驗條件各做3次重復仿真試驗，得顆粒平均仿真反彈高度H′1分別為113.233、124.778、139.987 mm，取平均值124.810 mm，與真實試驗條件得到的反彈高度相對誤差分別為1.53%、0.008%、1.93%。由此表明標定后的仿真結(jié)果與試驗結(jié)果基本吻合，故選取顆粒- 鋼板仿真恢復系數(shù)為e′w1=0.523。

通過標定后的微型薯- 鋼板仿真碰撞恢復系數(shù)e′w1=0.523，真實試驗測量中測定ew1=0.526，兩者相對誤差僅為0.57%，發(fā)現(xiàn)標定后的微型薯- 鋼板仿真碰撞恢復系數(shù)與真實試驗條件下測定的結(jié)果十分接近，其原因可能在于微型薯顆粒在與鋼板碰撞過程中因接觸時間極短、瞬間碰撞力極大，利用仿真進行模擬的物理環(huán)境條件與真實試驗條件趨于一致，故通過仿真標定的顆粒- 鋼板碰撞恢復系數(shù)與真實試驗條件所測值出現(xiàn)了近乎相等的情況；同時，在真實試驗中通過對不同高度微型薯開展碰撞恢復系數(shù)測定試驗發(fā)現(xiàn)，隨著釋放高度的增加，所得恢復系數(shù)雖變化不大但相應地會略有減小，分析其原因可能在于微型薯釋放高度越高其與鋼板碰撞時變形量越大，顆粒所需恢復時間就越長，而此時顆粒與鋼板接觸時間極短，微型薯顆粒因變形較大消耗較多能量，進而導致恢復系數(shù)值的降低。

2.1.2微型薯顆粒間恢復系數(shù)

為測定微型薯顆粒間碰撞恢復系數(shù)，引入兩顆類球形微型薯進行碰撞試驗[26]，如圖3a所示，其中顆粒a與顆粒b分別用尼龍繩套接，尼龍繩另一端固結(jié)于角鋼上，角鋼掛接在墻體上，并保證在顆粒自然懸垂狀態(tài)下，顆粒a、b處于同一高度且均在尼龍繩的束縛下做半徑r=640 mm的擺動。試驗過程中，以顆粒a、b 自然懸垂點為基點，提升顆粒a(保持尼龍繩處于拉伸狀態(tài))至基點垂直方向高度H0，本試驗分別設(shè)定H0為80、100、120 mm，顆粒b自然懸垂并保持靜止。以初速度為零釋放顆粒a，運動至最底端與顆粒b發(fā)生碰撞，顆粒a、b均繞著尼龍繩擺動至最高點，此時在垂直方向，顆粒a、b距基點高度分別為Ha、Hb，整個過程由高速攝像機記錄。根據(jù)碰撞恢復系數(shù)物理學定義得微型薯顆粒間碰撞恢復系數(shù)計算公式為

(7)

式中v1——碰撞前瞬間顆粒a瞬時速度

v′——碰撞后瞬間顆粒a瞬時速度

v0——碰撞后瞬間顆粒b瞬時速度

顆粒a在H0為80、100、120 mm等條件下各進行30次重復碰撞試驗得Ha為(5.01±1.48) mm、(5.93±2.08) mm、(6.56±1.70) mm，Hb為(43.01±6.75) mm、(52.67±11.55) mm、(60.46±8.94) mm，仿真過程中分別對應選取Ha為5.01、5.93、6.56 mm，Hb為43.01、52.67、60.46 mm作為驗證值，得顆粒a在3種釋放高度下微型薯顆粒間恢復系數(shù)ew2依次為0.483±0.009、0.480±0.011、0.476±0.007。由此表明，改變顆粒a初始釋放高度對于獲取微型薯顆粒間碰撞恢復系數(shù)值無顯著影響，故選取H0=100 mm作為初始條件開展標定仿真試驗，此條件下測得ew2取值范圍為0.131～0.589。

圖3 顆粒間恢復系數(shù)測定試驗Fig.3 Measurement test of particle- particle coefficient of restitution

考慮到在EDEM軟件中模擬微型薯顆粒在繩的束縛下擺動難度較大，仿真試驗采用等效顆粒碰撞模型，即仿真顆粒通過半徑r1=640 mm的半圓形管道進行約束，且穿過管道中心線的平面與水平面垂直，管道內(nèi)徑D2=10.50 mm，如圖3b所示。仿真過程中選取類球形微型薯顆粒仿真模型，分別在距軌道最低點垂直高度為100 mm的管道位置處及軌道最低點各生成1粒微型薯顆粒模型a、b且均保持靜止。經(jīng)預仿真試驗發(fā)現(xiàn)，管道本征參數(shù)對顆粒碰撞過程影響不顯著，將管道本征參數(shù)參照鋼板設(shè)定，并將微型薯- 管道接觸參數(shù)中的碰撞恢復系數(shù)、靜摩擦因數(shù)、滾動摩擦因數(shù)均設(shè)定為0。以初速度為0釋放微型薯顆粒a，模擬顆粒間碰撞過程，并利用EDEM軟件分析模塊求出兩顆粒碰撞后上升的最大高度Ha、Hb，經(jīng)過預仿真試驗并結(jié)合真實顆粒碰撞試驗條件下測得的顆粒間恢復系數(shù)ew2取值范圍確定試驗方案并得出結(jié)果(表2)。以微型薯顆粒間碰撞恢復系數(shù)x為試驗因素，顆粒a、顆粒b碰撞后上升最大高度ya、yb為評價指標，對結(jié)果進行曲線擬合，得擬合方程

(8)

表2 顆粒間碰撞恢復系數(shù)仿真試驗方案與結(jié)果Tab.2 Simulation test plan and results of particle- particlecoefficient of restitution

通過標定后的顆粒間仿真碰撞恢復系數(shù)e′w2=0.478，真實試驗測量中得ew2=0.480，兩者相對誤差為4.17%，發(fā)現(xiàn)標定后的顆粒間仿真碰撞恢復系數(shù)與真實試驗條件下測定的結(jié)果差異較小，其原因可能在于微型薯顆粒與顆粒間的碰撞過程因接觸時間極短、瞬間碰撞力極大，同時微型薯顆粒內(nèi)部具有相對較均勻的密度及質(zhì)量分布，利用仿真進行模擬的物理環(huán)境條件與真實試驗條件趨于一致，故通過仿真標定的顆粒間碰撞恢復系數(shù)與真實試驗條件所測值出現(xiàn)了差異較小的情況。

2.2 微型薯- 鋼板摩擦因數(shù)

建立微型薯- 鋼板摩擦因數(shù)測定模型，依據(jù)能量守恒定律，利用斜面滑動法、滾動法分別對微型薯- 鋼板靜摩擦因數(shù)和滾動摩擦因數(shù)進行分析和標定。

2.2.1微型薯- 鋼板靜摩擦因數(shù)

斜面滑動法是測量靜摩擦因數(shù)常用方法，質(zhì)量為m的物體重力分解為2個力：平行于斜面的力F(N)和垂直于斜面的力T(N)。當斜面傾角α(°)小于滑動臨界角時，F(xiàn)小于物體與斜面間的靜摩擦力f(N)，物體保持靜止，隨著斜面角度α的增加，F(xiàn)值越來越大，當α大于物體滑動臨界角時，F(xiàn)>f，物體將開始沿著斜面下滑。其中，靜摩擦因數(shù)μs與斜面傾角α的關(guān)系為

μs=f/T=mgsinα/(mgcosα)=tanα

(9)

試驗中取一張Q235鋼制平板作為試驗斜面，為防止單粒微型薯在斜面上滾動，將4顆類球形微型薯顆粒粘結(jié)在一起(半徑r=10 mm)放置在平板上(圖4a)，平板一側(cè)與水平試驗臺始終貼合并保持不動，緩慢勻速地將平板另一側(cè)抬起，當微型薯在平板上開始滑動時，利用數(shù)顯角度尺(精度0.05°，量程0°～360°)測量平板與試驗臺的夾角θ，重復試驗30次，取平均值得θ=(32.80±0.05)°，仿真過程中選取θ=32.80°，通過試驗測定微型薯- 鋼板靜摩擦因數(shù)μs1=0.646±0.058。

圖4 顆粒- 鋼板靜摩擦因數(shù)測定試驗Fig.4 Measurement test of potato minituber- steel plate coefficient of static friction1.組合顆粒 2.鋼制平板

在EDEM軟件中添加1個方形平板(長400 mm、寬500 mm)，設(shè)置其本征參數(shù)與Q235鋼板相同，并采用多球面組合方法生成4球組合微型薯仿真顆粒以模擬真實試驗條件下的粘結(jié)微型薯顆粒(圖4b)。在離散元仿真過程中，微型薯- 鋼板間仿真接觸參數(shù)設(shè)定如下：碰撞恢復系數(shù)選取已標定值0.523，滾動摩擦因數(shù)取值為0以消除其對斜面滑動仿真試驗的影響。通過改變方形平板的傾斜角度，在對應每個傾斜角度下通過二分法確定顆粒- 鋼板靜摩擦因數(shù)，試驗方案與結(jié)果如表3所示。以顆粒- 鋼板靜摩擦因數(shù)x2為試驗因素，傾斜角度y為評價指標，對試驗結(jié)果進行曲線擬合，得回歸方程

(10)

方程的決定系數(shù)R2=1，表明回歸方程準確有效。將y=32.80°代入方程(10)得x2=0.644，并以顆粒- 鋼板靜摩擦因數(shù)為0.644進行3次重復斜面滑動仿真試驗，得鋼制平板傾斜角度分別為32.75°、32.78°、32.78°，取平均值32.77°，與試驗所得鋼制平板傾斜角度相對誤差為0.09%。由此表明標定后的仿真結(jié)果與試驗結(jié)果基本一致，故選取微型薯- 鋼板靜摩擦因數(shù)μ′s1=0.644。

耳石癥是一種當頭部快速移動時出現(xiàn)眩暈和位置性眼震的眩暈癥,發(fā)病率隨年齡增長,且女性發(fā)病率高于男性,嚴重影響患者生活質(zhì)量。關(guān)于耳石癥的發(fā)病機理,目前主流觀點認為是由于耳石變形脫落(可能由重大外力撞擊或衰老導致)移動到其它平衡器官造成眩暈。臨床上常采用手法復位來使耳石復位進行治療,Epley手法復位和Barbecue翻滾手法復位是當前臨床上較為常見的用于治療耳石癥的復位手法,Epley手法復位通過移動患者頭部,在重力的作用下使耳石復位,Barbecue翻滾手法復位是根據(jù)耳石假說和前庭結(jié)構(gòu)解剖學關(guān)系來移動患者頭部使耳石復位。

表3 微型薯- 鋼板靜摩擦因數(shù)仿真試驗方案與結(jié)果Tab.3 Simulation test plan and results of potatominituber- steel plate coefficient of static friction

通過標定得微型薯- 鋼板靜摩擦因數(shù)μ′s1=0.644，真實試驗條件下測得μs1=0.646，兩者相對誤差僅為0.31%，發(fā)現(xiàn)標定后的微型薯- 鋼板靜摩擦因數(shù)與真實試驗條件下測定的結(jié)果差異較小，經(jīng)分析其原因可能在于通過斜面法測定微型薯- 鋼板靜摩擦因數(shù)，試驗過程中微型薯顆粒受到重力和靜摩擦力的共同作用始終與鋼板表面保持貼合，當鋼板斜面傾角達到滑動臨界角時，微型薯顆粒開始沿著斜面下滑，整個過程并沒有滾動摩擦力作為混合摩擦力的影響因素，顆粒受力相對單一，利用仿真模擬所創(chuàng)造的物理環(huán)境條件與真實試驗條件趨于一致，故通過仿真標定顆粒- 鋼板靜摩擦因數(shù)與真實試驗條件所測值基本一致。

2.2.2微型薯- 鋼板滾動摩擦因數(shù)

滾動摩擦是指當一個物體在另一物體表面作無滑動的滾動或有滾動的趨勢時，由于物體在接觸部分受壓發(fā)生形變而產(chǎn)生阻礙滾動的作用，其一般用阻力矩來度量，力的大小與物體的性質(zhì)、表面形狀以及滾動物體的質(zhì)量有關(guān)[11]。微型薯顆粒在鋼制平板上的滾動將產(chǎn)生滾動摩擦，試驗分別選取橢球形微型薯顆粒(長軸長度29 mm、短軸長度18 mm)和類球形微型薯顆粒(半徑10 mm)為試驗對象，從一定傾斜角度β的鋼制平板上在某一位置以初速度為0沿著斜面向下滾動(顆粒在斜面上的滾動距離為S)，并最終滾落至水平放置的鋼制平板上，由于受到滾動摩擦，微型薯顆粒在水平鋼板上滾動一段距離之后最終靜止，采用卷尺測定顆粒在水平鋼板上的滾動距離L，以L為評價指標，對微型薯- 鋼板滾動摩擦仿真系數(shù)進行標定(圖5a)。試驗過程中，假設(shè)所選取的微型薯種子為理想的橢球和圓球顆粒，微型薯顆粒做的是純滾動故認為其所受阻力僅為滾動摩擦力而不考慮滾動過程中靜摩擦力的影響，由能量守恒定律得[31]

mgSsinβ=μrmg(Scosβ+L)

(11)

通過預滾動試驗發(fā)現(xiàn)，當斜置鋼板傾斜角度β過小時，因選取的微型薯并不是理想的橢球或圓球，顆粒在斜面上無滾動或滾動一段距離S之后將靜置停留在斜面上；當斜置鋼板傾斜角度β過大時，顆粒滾落至平行放置的鋼板上時將產(chǎn)生彈跳，影響試驗結(jié)果，綜合考慮，本試驗選取斜置鋼板傾斜角度β=20°。同時，微型薯在斜置鋼板上運行距離S不能過小也不能過大，當斜向運行距離S過小時，顆粒滾落至水平鋼板上水平速度太小導致水平滾動位移L過短，不利于試驗測量；當斜向運行距離S過大時，顆粒滾落至水平鋼板上速度較大導致水平滾動位移L較長，由于微型薯顆粒不是理想的橢球體或圓球體，顆粒水平運行軌跡將不能保證在一條直線上，因此為提升試驗的可靠性及準確性，經(jīng)過預試驗，設(shè)定顆粒在斜面上滾動距離S=50 mm。通過大量預試驗發(fā)現(xiàn)兩類微型薯顆粒滾動距離基本一致，但由于橢球形微型薯顆粒因形狀及質(zhì)量分布不均，在滾動過程中極易出現(xiàn)運動軌跡不在同一條直線且軌跡線較復雜的現(xiàn)象，綜合考慮，本文以類球形微型薯為試驗對象開展微型薯- 鋼板滾動摩擦因數(shù)標定仿真試驗。根據(jù)上述試驗要求，進行30次重復滾動測量試驗，取平均值得顆粒在水平鋼板上的滾動距離為L=639 mm。將上述測定值代入式(11)中，得真實試驗條件下微型薯- 鋼板滾動摩擦因數(shù)μr1=0.010 8。

圖5 顆粒- 鋼板滾動摩擦因數(shù)測定試驗Fig.5 Measurement test of particle- steel plate coefficient of rolling friction1.傾斜鋼板 2.平置鋼板 3.微型薯顆粒 4.卷尺

在離散元仿真中，通過EDEM軟件分別添加傾斜放置(角度20°)和水平放置的矩形平板，且傾斜放置的平板底端與水平放置的平板相接觸，設(shè)置兩平板的本征參數(shù)與Q235鋼相同，沿著傾斜平板面向上距離底端S=50 mm處生成一半徑r=10 mm的球形微型薯仿真顆粒并以初速度為0沿著斜面向下滾動(圖5b)。顆粒在實際滾動過程中同時受到靜摩擦的作用，故在仿真模擬中微型薯- 鋼板接觸參數(shù)如下：碰撞恢復系數(shù)、靜摩擦因數(shù)均選取上文標定值e′w1=0.523、μ′s1=0.644。經(jīng)預仿真試驗得出微型薯- 鋼板滾動摩擦因數(shù)取值范圍為0.020 0～0.027 0，仿真試驗方案與結(jié)果如表4所示。以微型薯- 鋼板滾動摩擦因數(shù)x3為試驗因素，以微型薯仿真顆粒在水平平板上的滾動距離L為評價指標，對結(jié)果進行曲線擬合，建立擬合方程

(12)

方程決定系數(shù)R2=0.999 9，接近1，表明擬合方程準確可靠。將L=639 mm代入方程(12)得x3=0.022 1，設(shè)定微型薯- 鋼板滾動摩擦因數(shù)為0.022 1并在試驗對象為橢球形和類球形微型薯仿真顆粒兩種條件下各進行3次重復斜面滾動仿真試驗，得微型薯仿真顆粒在水平放置的平板上平均滾動距離分別為627 mm和638 mm，與試驗所得微型薯在水平鋼板上的滾動距離相對誤差各為1.88%、0.16%。由此表明標定后的仿真結(jié)果與試驗結(jié)果基本一致，故選取微型薯- 鋼板仿真滾動摩擦因數(shù)為μ′r1=0.022 1。

表4 微型薯- 鋼板滾動摩擦仿真試驗方案與結(jié)果Tab.4 Simulation test plan and results of potatominituber-steel plate coefficient of rolling friction

通過標定得微型薯- 鋼板間仿真滾動摩擦因數(shù)μ′r1=0.022 1，真實試驗測定得μr1=0.010 8，兩者相對誤差高達104.63%，發(fā)現(xiàn)標定后微型薯- 鋼板間滾動摩擦因數(shù)與真實試驗條件測定值差異較大，分析其原因可能是：①真實微型薯顆粒并不是理想的橢球或圓球體，在真實試驗條件下并不能保證顆粒在滾動過程中其運行軌跡為一條直線，用卷尺測定其滾動距離存在誤差，同時所選取的微型薯顆粒表皮以及鋼板表面特性與仿真模型有所差異。②在真實試驗條件下微型薯顆粒重力勢能的損耗并不只是來自滾動摩擦，還應包括顆粒- 鋼板間所產(chǎn)生的靜摩擦以及顆粒在滾動過程中與鋼板碰撞所損耗的能量，而仿真條件下所設(shè)定的則是理想的物理環(huán)境條件，即仿真模擬的環(huán)境條件與真實試驗條件差異較大，故出現(xiàn)標定值與真實試驗值相對誤差較高的情況。

2.3 微型薯顆粒間摩擦因數(shù)

堆積角是表征顆粒物料流動、摩擦等特性的宏觀參數(shù)，其中靜摩擦因數(shù)、滾動摩擦因數(shù)等接觸參數(shù)對結(jié)果影響顯著[22,31-33]，上文已對微型薯- 鋼板靜摩擦因數(shù)、滾動摩擦因數(shù)等仿真參數(shù)進行了標定，本試驗模型則用來對微型薯顆粒間靜摩擦因數(shù)、滾動摩擦因數(shù)等仿真參數(shù)進行分析和標定。

2.3.1堆積角測量試驗

采用箱體抽板法進行試驗[31,34](圖6a)，測量裝置(Q235鋼制材料)由箱體(長400 mm、寬300 mm、高500 mm)、擋板(長400 m、高500 mm)和底板(長1 000 mm、寬1 000 mm)組成，為使顆粒種子在堆積過程中充分擴散并在水平底板上形成無約束顆粒堆，應以較緩慢的速度向上抽提擋板[35]，設(shè)定向上抽提速度為0.05 m/s，待種群穩(wěn)定后，顆粒種群所形成的斜面與水平底板平面的夾角即為種群堆積角，試驗重復10次。

圖6 箱體測量裝置Fig.6 Measurement device with box1.擋板 2.箱體 3.種子 4.底板

圖7 堆積角圖像處理Fig.7 Image analysis by Matlab for stacking angle

堆積角試驗所采集的圖像，利用Matlab進行處理，因微型薯種子顆粒較大，應依次對圖像進行灰度處理、二值化處理、孔洞填充、提取邊界曲線，最后利用最小二乘法對邊界曲線進行直線擬合，如圖7所示，其中橫縱坐標分別為圖像水平像素點和垂直像素點，不具有實際量綱，擬合的直線斜率即為堆積角正切值，將正切值轉(zhuǎn)化為堆積角，由此得堆積角均值為(26.32±1.20)°，試驗中堆積角目標值選取26.32°。

2.3.2仿真參數(shù)分析和標定

進行仿真試驗時，以建立的微型薯仿真模型為基準，為避免生成過小的顆粒，根據(jù)實際測得的微型薯顆粒最小、最大尺寸，將單元球半徑限制在0.8～1.2倍的初始半徑之間。此外，由于應力波在微型薯顆粒中傳播受仿真參數(shù)的影響，各仿真中瑞利時步可能不同，因此在所有仿真中統(tǒng)一取20%瑞利時步。仿真中網(wǎng)格尺寸取3倍最小球形單元尺寸[23,36-37]。堆積角仿真試驗通過改變顆粒間靜摩擦因數(shù)、滾動摩擦因數(shù)，對堆積角進行標定，其它離散元仿真參數(shù)均采用上述所確定的經(jīng)驗值和標定參數(shù)，如表5所示。

表5 離散元仿真參數(shù)Tab.5 Parameters of discrete element simulation

將箱體測量裝置的Solidworks三維幾何模型轉(zhuǎn)換成.stp格式文件并導入EDEM中，以確保仿真中試驗裝置與實際試驗中所用試驗裝置保持一致；創(chuàng)建3個顆粒工廠，分別生成橢球小粒、橢球大粒、類球形3類微型薯顆粒，該顆粒工廠為虛擬表面(長400 mm，寬300 mm)，生成后的顆粒以0.5 m/s的速度下落并充滿上述試驗箱體內(nèi)，直至達到穩(wěn)定狀態(tài)，設(shè)置顆粒生成時間為3 s、生成方式為Dynamic，對于橢球小粒、橢球大粒及類球形3類微型薯顆粒，其生成速率分別為4 000、800、200顆/s。顆粒生成后，在EDEM中以0.05 m/s的速度垂直向上提升擋種板，種子顆粒將從箱體底端緩慢流出，最終在水平放置的鋼制底板上形成穩(wěn)定的顆粒堆，為確保顆粒堆達到穩(wěn)定狀態(tài)，將仿真運行時間設(shè)定為10 s(圖6b)。穩(wěn)定的顆粒堆形成之后，運用Matlab對堆積角圖像進行處理并讀取角度值。

經(jīng)過大量預仿真試驗，確定微型薯顆粒間靜摩擦因數(shù)仿真值范圍0.300～0.325，顆粒間滾動摩擦因數(shù)仿真值范圍0.030 0～0.032 6，應用Design-Expert軟件進行通用旋轉(zhuǎn)中心組合試驗設(shè)計，選取微型薯顆粒間靜摩擦因數(shù)和滾動摩擦因數(shù)為試驗因素，堆積角和仿真時間為評價指標，設(shè)置的仿真試驗因素編碼如表6所示。

表6 仿真試驗因素編碼Tab.6 Coding of factors

試驗方案及結(jié)果如表7所示，A、B為因素編碼值。其中仿真時間是指堆積角仿真試驗時從運行開始至顆粒堆達到穩(wěn)定狀態(tài)的時間，選用仿真時間作為評價指標是考慮到不同的試驗組合運行時間有所差異，仿真運行時間越長表明仿真耗散的時間越多，同時生成的仿真文件所占存儲空間越大，本仿真試驗設(shè)定全部顆粒的平均速度小于1.0×10-4m/s時顆粒堆達到穩(wěn)定狀態(tài)，找出這一時間點即為模擬仿真時間。

表7 仿真試驗方案設(shè)計與結(jié)果Tab.7 Simulation experiment plan and result

利用Design-Expert 8.0.6軟件對仿真試驗結(jié)果進行方差分析，發(fā)現(xiàn)兩試驗模型均為極顯著(P<0.01)，說明此試驗合理有效。對指標有顯著影響的所有因素均已考慮到，得出擬合較好且具有實際分析意義的因素編碼值的回歸方程[38]

(13)

利用Design-Expert 8.0.6軟件的優(yōu)化模塊，使仿真結(jié)果最接近試驗所得的微型薯顆粒堆積角，對回歸模型進行有約束目標的優(yōu)化求解[39]。

在優(yōu)化模塊的條件設(shè)置中對兩評價指標的重要程度進行設(shè)定，其中堆積角設(shè)定為“+++++”，仿真時間設(shè)定為“+”。由此得到優(yōu)化結(jié)果為：微型薯顆粒間靜摩擦因數(shù)為μ′s2=0.325，微型薯顆粒間滾動摩擦因數(shù)為μ′r2=0.030 0，此時堆積角仿真結(jié)果為θ=26.32°，仿真時間為t=5.466 36 s。

將標定后的顆粒間摩擦因數(shù)代入EDEM軟件中進行3次重復堆積角仿真試驗，得堆積角分別為26.16°、26.39°、26.45°，取平均值26.33°，與試驗測量值相對誤差為0.04%；仿真時間分別為5.250 01、5.250 01、5.250 01 s，取平均值5.250 01 s，小于優(yōu)化結(jié)果所得的仿真時間5.466 36 s。由此表明標定的仿真結(jié)果與試驗結(jié)果基本吻合，故微型薯顆粒間靜摩擦因數(shù)、微型薯顆粒間滾動摩擦因數(shù)的仿真參數(shù)標定值分別為μ′s2=0.325、μ′r2=0.030 0。

采用堆積角試驗對微型薯顆粒間靜摩擦因數(shù)和滾動摩擦因數(shù)進行標定而不同于上述顆粒- 鋼板摩擦因數(shù)模型中利用斜面滑動與滾動法對其參數(shù)進行試驗測量與仿真模擬相結(jié)合的方法來標定，是由于考慮到微型薯種子顆粒形狀不規(guī)則，種子表面有凹凸，在進行相關(guān)測定時會產(chǎn)生太多干擾因素，直接測定微型薯顆粒間相關(guān)摩擦因數(shù)誤差太大，故利用顆粒種群堆積角等宏觀現(xiàn)象來表征其相關(guān)物理特性。

3 驗證試驗

3.1 圓筒提升法堆積角試驗

采用無底圓筒提升法進行堆積角驗證試驗[13,17,34](圖8a)，測量裝置(Q235鋼制材料)由無底圓筒(內(nèi)徑250 mm、高度750 mm)和底板(長1 000 mm、寬1 000 mm)組成，以0.05 m/s速度向上提升圓筒，待微型薯顆粒種群穩(wěn)定后，測定顆粒種群所形成的斜面與水平底板平面的夾角即為種群堆積角，試驗所采集的圖像利用Matlab進行處理，試驗重復10次，求得堆積角平均值為(23.96±1.15)°，其目標值取23.96°。將上述標定獲取的離散元本征參數(shù)和接觸參數(shù)輸入EDEM中進行無底圓筒提升法堆積角仿真試驗(圖8b)，其中參數(shù)設(shè)定方法同2.3.2節(jié)，重復仿真試驗3次，測定結(jié)果為：24.10°、25.88°、22.39°，求均值得(24.12±1.75)°，取值24.12°與試驗得到的堆積角相對誤差為0.67%，表明標定后的仿真結(jié)果與試驗結(jié)果基本吻合，由此說明通過標定的方法可以找到仿真模型與試驗微型薯顆粒間在物理特性上的對應關(guān)系，根據(jù)等效原則建立了仿真與實際試驗之間的聯(lián)系。

圖8 圓筒測量裝置Fig.8 Measurement device with steel cylinder1.無底圓筒 2.底板 3.種子

圖9 落種試驗Fig.9 Falling seeds test1.種箱 2.振動排序單元 3.落種單元 4.開溝單元

3.2 落種試驗

為進一步驗證參數(shù)的準確性及適用性，選用一款微型馬鈴薯振動排序播種裝置(圖9a)進行落種試驗。該播種裝置主要由種箱、振動排序單元、落種單元以及開溝單元構(gòu)成，其工作原理為微型薯種子從種箱排種口排出并落至振動排序單元，微型薯在以周期性振動的振動排序單元作用下實現(xiàn)薯種的單列排序并輸送至落種單元。落種單元將單列排序的薯種播落至開溝單元所開的種溝內(nèi)，完成整個播種過程。其中種箱和振動排序單元均由Q235鋼制材料制成。本文在真實試驗條件下將1500粒微型薯倒入種箱內(nèi)，薯種顆粒從種箱排種口排出至振動排序播種單元內(nèi)，并設(shè)置振動排序單元進行3個周期的振動，待種群穩(wěn)定后觀察微型薯顆粒分布情況(圖9b)。依據(jù)真實試驗相關(guān)參數(shù)進行微型薯仿真落種試驗(圖9c)，兩種條件下均做3次重復試驗。選取落入播種單元內(nèi)的顆粒數(shù)量m、顆粒種群在落種口形成的堆積角α作為種子分布情況的2個關(guān)鍵特征進行對比，測量值如表8所示。

表8 真實試驗與仿真種群關(guān)鍵特征參數(shù)對比Tab.8 Comparison of key feature parameters foractual and simulation results

結(jié)果表明測得播種單元內(nèi)的顆粒數(shù)量在仿真與真實試驗條件下的種量相對誤差為4.4%；微型薯種群堆積角在仿真與真實試驗條件下的堆積角相對誤差為4.8%。仿真條件下微型薯顆粒分布情況與真實試驗條件基本相同，說明通過標定后的相關(guān)微型薯顆粒種子物性參數(shù)可以應用于離散元仿真中，為后續(xù)仿真模擬提供基礎(chǔ)。

4 結(jié)論

(1)以微型薯種子作為研究對象，創(chuàng)建微型薯模型，以此為基礎(chǔ)建立微型薯離散元參數(shù)獲取模型，利用試驗與仿真模擬相結(jié)合的方法，根據(jù)各仿真接觸模型的相關(guān)試驗數(shù)據(jù)進行回歸分析，建立回歸方程，依次求取微型薯- 鋼板碰撞恢復系數(shù)為0.523，微型薯顆粒間碰撞恢復系數(shù)為0.478，微型薯- 鋼板靜摩擦因數(shù)為0.644，微型薯- 鋼板滾動摩擦因數(shù)為0.022 1，微型薯顆粒間靜摩擦因數(shù)為0.325，微型薯顆粒間滾動摩擦因數(shù)為0.030 0。

(2)采用無底圓筒提升法進行堆積角試驗、采用落種法觀察微型薯顆粒種群分布情況，對離散元仿真參數(shù)進行驗證，結(jié)果為采用無底圓筒提升法獲取堆積角為24.12°，相對誤差為0.67%；落種法觀察顆粒種群分布情況，對比2個關(guān)鍵特征，其數(shù)值相對誤差均在4.80%以內(nèi)。由此表明該方法可全面、系統(tǒng)地找到微型薯顆粒離散元仿真參數(shù)，從而為微型馬鈴薯相關(guān)播種機具設(shè)計和優(yōu)化提供了理論依據(jù)。

1 NIEDERHAUSER J S. International cooperation and the role of the potato in feeding the world[J]. American Potato Journal，1993，70(5): 385-403.

2 張勇，李麗霞，韓宏宇，等. 馬鈴薯微型原種播種機的研究現(xiàn)狀與展望[J]. 安徽農(nóng)業(yè)科學，2015，43(4): 372-375.

3 呂金慶，楊穎，李紫輝，等. 舀勺式馬鈴薯播種機排種器的設(shè)計與試驗[J]. 農(nóng)業(yè)工程學報，2016，32(16): 17-25.

Lü Jinqing, YANG Ying, LI Zihui, et al. Design and experiment of cup-belt type potato seed-metering device[J]. Transactions of the CSAE, 2016,32(16): 17-25. (in Chinese)

4 王云霞，梁志杰，崔濤,等. 玉米分層施肥器結(jié)構(gòu)設(shè)計與試驗[J/OL]. 農(nóng)業(yè)機械學報，2016，47(增刊): 163-169. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=2016s025&journal_id=jcsam.DOI：10.6041/j.issn.1000-1298.2016.S0.025.

WANG Yunxia, LIANG Zhijie, CUI Tao, et al. Design and experiment of layered fertilization device for corn[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(Supp.): 163-169. (in Chinese)

5 閆銀發(fā)，孟德興，宋占華，等. 槽輪式補飼機顆粒動力學數(shù)值模擬與試驗[J/OL]. 農(nóng)業(yè)機械學報，2016，47(增刊): 249-253. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=2016s038&flag=1&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2016.S0.038.

YAN Yinfa, MENG Dexing, SONG Zhanhua, et al. Particle kinetic simulation and experiment for flute-wheel feeding machine[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(Supp.): 249-253. (in Chinese)

6 劉月琴，趙滿全，劉飛,等. 基于離散元的氣吸式排種器工作參數(shù)仿真優(yōu)化[J/OL]. 農(nóng)業(yè)機械學報，2016，47(7): 65-72. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160710&journal_id=jcsam.DOI：10.6041/j.issn.1000-1298.2016.07.010.

LIU Yueqin, ZHAO Manquan, LIU Fei, et al. Simulation and optimization of working parameters of air suction metering device based on discrete element[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(7): 65-72. (in Chinese)

7 黃玉祥，杭程光，苑夢嬋，等. 深松土壤擾動行為的離散元仿真與試驗[J/OL]. 農(nóng)業(yè)機械學報，2016，47(7): 80-88. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160712&journal_id=jcsam.DOI：10.6041/j.issn.1000-1298.2016.07.012.

HUANG Yuxiang, HANG Chengguang, YUAN Mengchan, et al. Discrete element simulation and experiment on disturbance behavior of subsoiling[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(7):80-88. (in Chinese)

8 LENAERTS B, AERTSEN T, TIJSKENS E, et al. Simulation of grain-straw separation by discrete element modeling with bendable straw particles [J]. Computers & Electronics in Agriculture, 2014, 101(2): 24-33.

9 MARTINA C L, BOUVARDA D, SHIMAB S. Study of particle rearrangement during powder compaction by the discrete element method[J]. Journal of the Mechanics and Physics of Solids, 2003, 51(4): 667-693.

10 劉凡一，張艦，李博，等. 基于堆積試驗的小麥離散元參數(shù)分析及標定[J]. 農(nóng)業(yè)工程學報，2016，32 (12): 247-253.

LIU Fanyi, ZHANG Jian, LI Bo,et al. Calibration of parameters of wheat required in discrete element method simulation based on repose angle of particle heap[J]. Transactions of the CSAE, 2016, 32(12): 247-253. (in Chinese)

11 崔濤，劉佳，楊麗，等. 基于高速攝像的玉米種子滾動摩擦特性試驗與仿真[J]. 農(nóng)業(yè)工程學報，2013，29(15): 34-41.

13 王憲良，胡紅，王慶杰，等. 基于離散元的土壤模型參數(shù)標定方法[J/OL]. 農(nóng)業(yè)機械學報，2017，48(12): 78-85. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20171209&journal_id=jcsam.DOI：10.6041/j.issn.1000-1298.2017.12.009.

WANG Xianliang, HU Hong, WANG Qingjie, et al. Calibration method of soil contact characteristic parameters based on DEM theory[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(12): 78-85. (in Chinese)

14 李守巨，李德，于申基. 基于宏觀實驗數(shù)據(jù)的堆石料細觀本構(gòu)模型參數(shù)反演[J]. 山東科技大學學報: 自然科學版，2015，34(5): 20-26.

15 CAMPOS A V P, OLIVEIRA O S, FERREIRA L V, et al. Dem simulations of dynamic angle of repose of acerola residue: a parametric study using a response surface technique[J]. Blucher Chemical Engineering Proceedings, 2015, 1(2): 11326-11333.

16 WANG X, YU J, LV F, et al. A multi-sphere based modelling method for maize grain assemblies[J]. Advanced Powder Technology, 2017, 28(2):584-595.

17 WANG L, LI R, WU B, et al. Determination of the coefficient of rolling friction of an irregularly shaped maize particle group using physical experiment and simulations[J]. Particuology, 2018,38:185-195.

18 COETZEE C J, DNJ E. Calibration of discrete element parameters and the modelling of silo discharge and bucket filling[J]. Computers & Electronics in Agriculture, 2009, 65(2): 198-212.

19 邱白晶，姜國微，楊寧，等. 水稻籽粒流對承載板沖擊過程離散元分析[J]. 農(nóng)業(yè)工程學報，2012，28(3): 44-49.

20 CUNHA R N, SANTOS K G, LIMA R N, et al. Repose angle of monoparticles and binary mixture: an experimental and simulation study[J]. Powder Technology, 2016, 303: 203-211.

21 杜欣，曾亞武，高睿，等. 基于CT掃描的不規(guī)則外形顆粒三維離散元建模[J]. 上海交通大學學報，2011，45(5): 711-715.

22 石林榕，吳建民，孫偉，等. 基于離散單元法的水平圓盤式精量排種器排種仿真試驗[J]. 農(nóng)業(yè)工程學報，2014，30(8): 40-48.

SHI Linrong, WU Jianmin, SUN Wei, et al. Simulation test for metering process of horizontal disc precision metering device based on discrete element method[J]. Transactions of the CSAE, 2014, 30(8): 40-48. (in Chinese)

23 ASAE S368.4- 2000. Compression test of food materials of convex shape [S].2000.

24 YAN Z, WILKINSON S K, STITT E H, et al. Discrete element modelling (DEM) input parameters: understanding their impact on model predictions using statistical analysis[J]. Computational Particle Mechanics, 2015, 2(3): 283-299.

25 ?？担穱纻?，羅敏，等. 雙層種箱式馬鈴薯排種裝置設(shè)計與試驗[J]. 農(nóng)業(yè)工程學報，2016，32(20): 32-39.

NIU Kang, YUAN Yanwei, LUO Min, et al. Design and experiment of potato metering device with double-deck seed tank[J]. Transactions of the CSAE, 2016, 32(20): 32-39. (in Chinese)

27 王國強，郝萬軍，王繼新. 離散單元法及其在EDEM上的實踐[M]. 西安： 西北工業(yè)大學出版社，2010.

28 文愿運，劉馬林，劉榮正，等. 顆粒離散單元法數(shù)值模擬與典型實驗對比研究[J]. 中國粉體技術(shù)，2015(3): 1-5.

29 韓燕龍，賈富國，唐玉榮，等. 顆粒滾動摩擦系數(shù)對堆積特性的影響[J]. 物理學報，2014，63(17): 165-171.

HAN Yanlong, JIA Fuguo, TANG Yurong, et al. Influence of granular coefficient of rolling friction on accumulation characteristics[J]. Acta Physica Sinica, 2014, 63(17): 165-171. (in Chinese)

30 王俊，許乃章，胥芳. 桃子沖擊力學特性及其與桃子硬度的數(shù)學模型[J]. 農(nóng)業(yè)機械學報，1994，25(4): 58-62.

31 ZHOU Y C, XU B H, YU A B, et al. An experimental and numerical study of the angle of repose of coarse spheres [J]. Powder Technology, 2002, 125(1): 45-54.

32 YANG J, XUAN K, BIN L I, et al. Study on release mechanism of inhibitory components from cinnamon and clove powders [J]. Journal of Food Safety, 2012, 32(2): 189-197.

33 GHODKI B M, GOSWAMI T K. DEM simulation of flow of black pepper seeds in cryogenic grinding system [J]. Journal of Food Engineering, 2017, 196: 36-51.

34 張銳，韓佃雷，吉巧麗，等. 離散元模擬中沙土參數(shù)標定方法研究[J/OL]. 農(nóng)業(yè)機械學報，2017，48(3): 49-56. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170306&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2017.03.006.

ZHANG Rui, HAN Dianlei, JI Qiaoli, et al. Calibration methods of sandy soil parameters in simulation of discrete element method[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 49-56. (in Chinese)

35 ZHOU Y C, WRIGHT B D, YANG R Y, et al. Rolling friction in the dynamic simulation of sandpile formation [J]. Physica A: Statistical Mechanics & Its Applications, 1999, 269(2-4): 536-553.

36 賴慶輝，馬文鵬，劉素，等. 氣吸圓盤式微型薯排種器充種性能模擬與試驗[J/OL]. 農(nóng)業(yè)機械學報，2017，48(5): 44-53. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170505&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2017.05.005.

LAI Qinghui, MA Wenpeng, LIU Su, et al. Simulation and experiment on seed-filling performance of pneumatic disc seed-metering device for mini-tuber[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(5): 44-53. (in Chinese)

37 MARIGO M, STITT E H. Discrete element method (DEM) for industrial applications: comments on calibration and validation for the modelling of cylindrical pellets [J]. Kona, 2015, 32:236-252.

38 袁志發(fā)，周靜芋. 試驗設(shè)計與分析[M]. 北京:高等教育出版社，2000: 292-303.

39 徐中儒. 農(nóng)業(yè)試驗最優(yōu)回歸設(shè)計[M]. 哈爾濱:黑龍江科技出版社，1998: 181-187.