顆粒肥料離散元仿真邊界參數系統化研究

劉彩玲 魏 丹 宋建農 黎艷妮 都 鑫 張福印

(中國農業大學農業部土壤-機器-植物系統技術重點實驗室， 北京 100083)

0 引言

離散元法(Discrete element method，DEM)是一種基于分子動力學的顆粒物料分析方法，近年來，隨著計算機技術的發展，該方法在農業裝備研究上應用越來越多[1-3]。顆粒肥料是一種容易施用的優質高效肥料，便于機械化作業，在施肥機械上有著廣泛的應用，因此，全面系統研究顆粒肥料離散元仿真參數，有助于離散元法在施肥裝備研究中的應用[4]。

在離散元顆粒系統模擬過程中，顆粒物料固有的物理力學參數是影響顆粒運動特性的重要因素，正確的模型參數對模擬結果的正確性具有重要的意義[5-6]。NAKASHIMA等[7]通過模擬重力對沙堆顆粒休止角的影響指出，沙堆堆積角受重力影響很小，主要與顆粒半徑、摩擦因數等因素有關。ZHOU等[8]模擬了玻璃圓球顆粒在長方體容器中落料堆積過程，發現顆粒間及顆粒與容器壁滑動和滾動摩擦因數是影響堆積形態的關鍵因素。PERSSON等[9]在評價離散元預測顆粒系統流動性時觀測到顆粒間不同動靜摩擦因數值直接影響落料速率及堆積角變化。GONIVA等[10]采取顆粒堆積角為模擬指標發現，采用具有滾動阻力模型的離散元方法對模擬結果有很大改善。韓燕龍等[6]通過對稻谷休止角的仿真，對稻谷顆粒間滾動摩擦因數進行了研究；崔濤等[11]基于高速攝像技術研究了玉米種子滾動摩擦特性，并通過玉米種子休止角試驗對EDEM邊界參數進行標定；王云霞等[12]研究了3種不同形態玉米模型的休止角仿真試驗，確定了休止角與顆粒間摩擦特性的相關模型。以上國內外眾多學者的研究表明，離散元模型邊界參數對仿真精確性具有十分重要的影響，而且通過離散單元法研究散體顆粒堆積過程能有效模擬物料物理力學參數、顆粒接觸模型變化時對散體運動特性的影響情況。

離散元仿真模型的邊界參數包括材料的基礎物理參數(如泊松比、剪切模量、實體密度等)和材料間接觸參數(如碰撞恢復系數、靜摩擦因數以及滾動摩擦因數)。目前，仿真邊界參數主要取自試驗測定或經驗值。LIU等[13]基于彈跳試驗進行了水稻恢復系數的測量；BALEVIIUS等[14]通過滑動試驗測量了豌豆與有機玻璃間的靜摩擦因數；GONZLEZ-MONTELLANO等[15]對離散元仿真中部分參數的直接測量方法進行研究，結果表明，對于玉米及橄欖核這類外形不規則的顆粒材料，直接測量誤差很大。對于顆粒間接觸特性參數，韓燕龍等[6]將稻谷粘結在傾斜板面上制成顆粒板，測定稻谷間靜摩擦因數；崔濤等[11]將玉米種子粘貼在硬紙板上測定玉米種子間滾動摩擦因數，采用顆粒板法無法保證完全一致的顆粒層厚度且表面也有凹凸，種子運動過程中的彈跳、碰撞無法避免，因此直接測定比較困難，誤差較大。

綜上所述，離散元仿真模型邊界參數是影響顆粒運動的重要因素，邊界參數的確定對仿真的準確性具有重要的意義，但目前邊界仿真參數主要取自經驗值，大部分參數直接測定比較困難，特別是對于材料間接觸參數的測定誤差較大，邊界參數的研究還缺少系統性。顆粒堆積表征了顆粒物料流動、摩擦等特性，是研究散體力學問題的重要內容[16]，同時離散元法是對顆粒隨機堆積過程進行計算模擬，認識堆積的細觀力學機理和評估所采用模型適用性的有效方法[17]，基于此，本文以反映散體顆粒群堆積過程宏觀特征的休止角為宏觀現象，以尿素顆粒為研究對象，利用離散元仿真軟件EDEM系統研究尿素顆粒仿真邊界參數對休止角的影響，對重要邊界參數進行分析與確定，以提高離散元仿真精度，縮短設計周期。

1 休止角測定方法

采用ABS無底圓柱桶法測定肥料顆粒休止角，休止角仿真設置如圖1所示。參考文獻[17]確定圓筒直徑44 mm、高135 mm，提升速度為0.05 m/s。研究中對尿素顆粒進行6組試驗，每組測定均從垂直交叉的4個方向拍攝圖像并求其平均值，最后將6組平均值作為最終測定的休止角。

圖1 休止角仿真設置Fig.1 Rest angle simulation settings1.ABS板 2.ABS無底圓筒 3.肥料堆

為減小測量休止角誤差，利用Matlab讀取顆粒堆積圖像(圖2)，依此對圖像進行灰度處理、二值化處理，最后提取邊界點，邊界點上的連線即為顆粒堆積的邊界曲線，利用最小二乘法對邊界點進行直線擬合，擬合的直線斜率即為所要測得的休止角的正切值。

圖2 尿素顆粒堆積角測量Fig.2 Measurement of stacking angle of urea particles

2 顆粒及接觸力學模型

2.1 試驗對象及顆粒仿真模型的建立

試驗選取干燥的大顆粒尿素為研究對象，外觀形態接近球體，采用精度為0.01 mm的千分尺測得其算術平均粒徑為3.35 mm，大顆粒尿素直徑分布范圍為1.7～4.75 mm，千粒質量23.762 g，體積密度為1 330 kg/m3。

以3.35 mm作為顆粒模型直徑，在EDEM中構建如圖3a所示的顆粒離散元模型。根據大顆粒尿素直徑分布范圍，設置顆粒工廠顆粒尺寸按正態分布方式生成，將參數Mean設置為1，即正態分布的顆粒平均直徑為尿素顆粒模型直徑，標準差 Std Dev 設置為0.1，Capped用來設置生成顆粒的尺寸極限，最小為平均值的0.5倍，最大為平均值的1.4倍，尺寸極限與篩分法確定的尿素顆粒尺寸分布一致。具體設置如圖3b所示。

Scale By指顆粒尺寸縮放標準為顆粒半徑。

圖3 尿素顆粒模型以及生成方式Fig.3 Urea particle model and its generation method

時間步長與顆粒粒徑、運動速度、密度、剪切模量有關，通常設定為Rayleigh步長的5%～30%，本模型中顆粒排列緊密，確定時間步長為Rayleigh步長的25%。根據圓筒上升速度并考慮顆粒達到穩定狀態的需要，設定仿真總時長為15 s。

2.2 顆粒堆積的接觸力學模型

假設肥料表面沒有粘附力，顆粒運動特性取決于其物體或相鄰顆粒間接觸點位置會發生的微小彈性形變，因此顆粒堆積過程選擇Hertz-Mindlin無滑動接觸力學模型(圖4)，接觸點法向和切向碰撞力均簡化為彈簧和阻尼器的并聯。

圖4 接觸力學模型Fig.4 Model of contact force

基于上述模型，顆粒在堆積過程中主要受自身重力、顆粒間法向力Fn、切向力Ft、法向阻尼力Fnd和切向阻尼力Ftd，其中切向力與庫倫摩擦力μsFn(μs為靜摩擦因數)有關，各力大小為

(1)

Ft=-Stδ

(2)

(3)

(4)

(5)

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

R*——等效半徑

α——法向重疊量

ε——恢復系數δ——切向重疊量

Sn——法向剛度

m*——等效質量

vnrel——法向相對速度

St——切向剛度

vtrel——切向相對速度

另外，顆粒還受到切向力力矩和滾動摩擦力矩，其計算公式為

Tt=Ri(Ft+Ftd)

(6)

Tn=-μrFnRiωi

(7)

式中μr——滾動摩擦因數

Ri——質心到接觸點間距離

ωi——接觸點處物體單位角速度

顆粒在其所受力和力矩作用下發生移動和滾動。當切向力大于庫倫摩擦力時顆粒在接觸表面發生滑動，切向力矩與滾動摩擦力矩綜合作用決定顆粒的滾動運動。

3 邊界參數的系統化研究

為研究各邊界參數對肥料休止角的影響是否顯著，設計了休止角影響因素的Plackett-Burman篩選試驗。

3.1 休止角影響因素的篩選與分析

3.1.1休止角影響因素的Plackett-Burman篩選試驗

應用Design-Expert 8.0軟件進行Plackett-Burman試驗設計，參考文獻[4,18-21]確定ABS仿真參數為：泊松比0.394，剪切模量8.9×108Pa，密度為1 060 kg/m3；尿素密度為1 330 kg/m3，其余待定8個仿真邊界參數作為篩選對象，并給出各因素的范圍如表1所示。仿真中時間步長設置為Rayleigh步長的25%，仿真計算網格尺寸為最小顆粒直徑的4倍。

表1 Plackett-Burman試驗因素Tab.1 Parameters of Plackett-Burman test

3.1.2篩選試驗結果與分析

Plackett-Burman試驗設計與仿真休止角結果如表2所示，試驗中具有12組-1與+1水平點，1組中心點仿真試驗。

表2 Plackett-Burman試驗設計與結果Tab.2 Design and results of Plackett-Burman test

利用Design-Expert 8.0軟件對試驗結果進行方差分析，結果(表3)表明，對休止角影響重要性次序依次為尿素顆粒間滾動摩擦因數、尿素顆粒間靜摩擦因數、尿素顆粒與ABS間靜摩擦因數，其中，前兩者高度顯著。

表3 Plackett-Burman試驗顯著性分析Tab.3 Significance analysis of Plackett-Burman test

注：*表示顯著，** 表示高度顯著。

為探討3種顯著性邊界參數對休止角的影響，模擬了不同參數下尿素顆粒堆積狀態的變化，確定仿真參數為[4，19-21]：尿素剪切模量3.48×107Pa，泊松比0.25，尿素顆粒與顆粒間、與ABS間碰撞恢復系數分別為0.35和0.60，尿素顆粒與ABS間滾動摩擦因數為0.41，作為不變參數時顆粒間靜摩擦因數為0.27，顆粒與ABS板間靜摩擦因數為0.33，顆粒間滾動摩擦因數為0.22，其余參數見3.1.1節，仿真結果如圖5～7所示。

圖5 不同顆粒間滾動摩擦因數顆粒堆積圖Fig.5 Particle packing diagrams of different rolling friction coefficients between particles

圖6 不同顆粒間靜摩擦因數顆粒堆積圖Fig.6 Particle packing diagrams of different static friction coefficients between particles

圖7 不同顆粒與ABS板間靜摩擦因數顆粒堆積圖Fig.7 Particle packing diagrams of different static friction coefficients between particles and ABS plates

圖5a、5b、5c分別是尿素顆粒間滾動摩擦因數為0.01、0.26、0.50時尿素顆粒堆積狀態圖，休止角隨顆粒表面滾動摩擦因數增大而增大，當滾動摩擦因數較小時顆粒堆較分散，主要由于顆粒下落堆積過程中，顆粒間接觸作用劇烈，中心顆粒易于形成穩定的堆積狀態，對邊界顆粒產生排擠作用，邊界顆粒向外部擴散呈松散狀態，外層顆粒具有較大的旋轉動能，對顆粒的轉動慣量抑制作用較小，因此較小的顆粒間滾動摩擦因數時形成的休止角較小，顆粒堆較為分散。相反，較大的顆粒間滾動摩擦因數不利于邊界顆粒的擴散，降低了外層顆粒旋轉動能，對顆粒的轉動慣量抑制作用較大，從而使顆粒沿中心軸即顆粒堆的高度方向堆積，利于形成較大的休止角。

圖6a、6b、6c分別是尿素顆粒間靜摩擦因數為0.20、0.45、0.70時顆粒堆積狀態圖，相比滾動摩擦因數，靜摩擦因數較小時顆粒堆分散現象更為明顯，較小的顆粒間靜摩擦因數使顆粒堆發生變異或不易形成堆狀，主要由于顆粒間靜摩擦因數較小時，庫侖摩擦力較小，切向力較小，導致切向力不足以支持顆粒自身重力，顆粒在接觸表面易發生滑動，不易形成穩定堆積體；顆粒間靜摩擦因數較大時，庫倫摩擦力較大，顆粒在接觸表面相對靜止，形成穩定堆積體。

圖7a、7b、7c分別是尿素顆粒與ABS板間的靜摩擦因數為0.10、0.35、0.60時顆粒的堆積狀態圖，當顆粒與接觸面間的靜摩擦因數較小時，顆粒在接觸面間不能形成穩定的顆粒堆，當摩擦因數為0.35時，所形成的尿素顆粒堆較為集中，邊緣無分散現象。在散粒體堆積的過程中，重力作用使中底部散體顆粒間相互擠壓變形，形成強力鏈，且相互鏈接形成架拱保證顆粒堆的穩定性。當顆粒與接觸面間的靜摩擦因數較小時，底部顆粒與接觸面間的作用力不足以支撐上部顆粒的擠壓變形而使底部的強力鏈破壞，進而顆粒堆分散；隨著摩擦因數增大，底部顆粒與接觸面間作用力變大，底部顆粒易于形成強力鏈，使得顆粒向高度方向堆積，形成較大休止角，顆粒受到的接觸力作用線處于摩擦角的范圍之內，力鏈中的顆粒處于自鎖狀態，提高了顆粒與接觸面間接觸穩定性，形成穩定的堆積體。

3.2 重要仿真邊界參數的確定

上述分析結果表明，尿素顆粒間滾動摩擦因數以及顆粒間、顆粒與ABS間靜摩擦因數是影響顆粒堆積過程的重要因素，因此三者取值是提高仿真精度的關鍵。

3.2.1尿素顆粒間靜摩擦因數的標定

顆粒間細觀參數很難通過試驗方法直接準確測定，目前常用的方法是通過標定來確定[22]。采用自然休止角作為宏觀特征，依據第1節中介紹的休止角測定方法，在EDEM中進行不同尿素顆粒間靜摩擦因數下休止角仿真測定試驗，尿素顆粒間滾動摩擦因數和尿素顆粒與ABS間靜摩擦因數分別為0.30和0.33，其余仿真參數同3.1.2節，每組參數下取6組重復試驗平均值作為最終測定的休止角，得休止角ψ與靜摩擦因數μspp的關系如圖8，其擬合方程為

ψ=35.883μspp+18.772 (R2=0.984 1)

(8)

圖8 顆粒間靜摩擦因數與休止角關系曲線Fig.8 Relationship curves between urea-urea static friction coefficient and repose angle

二者呈線性變化趨勢，方程擬合精度較高(R2=0.984 1)，因此可用于標定尿素顆粒間靜摩擦因數。

在與仿真相同試驗條件下測定尿素顆粒實際休止角，試驗對象同2.1節，休止角測定方法如第1節所述，得到尿素顆粒實際休止角均值30.43°。將試驗測定值代入擬合方程(8)，計算得到標定的尿素顆粒間靜摩擦因數μspp=0.32。

3.2.2尿素顆粒與ABS間靜摩擦因數的測定

利用實驗室自制靜摩擦因數測量儀測定尿素顆粒與ABS板間靜摩擦因數(圖9)。將ABS板固定在支撐板上，為避免肥料顆粒滾落而非滑落，盡量選擇扁平肥料(圖9b)，調節絲杠螺母機構，當顆粒剛好下滑時記錄升降板角度，其正切值即為尿素顆粒與ABS板間靜摩擦因數，測量10次求平均值，為0.33。

圖9 靜摩擦因數測量Fig.9 Measurement of static friction coefficient

3.2.3尿素顆粒間滾動摩擦因數的標定

由于摩擦理論、測量設備的不完備，顆粒滾動摩擦因數還沒有準確、成熟的測量方法。目前顆粒間接觸參數(靜摩擦因數和滾動摩擦因數)采用粘結在傾斜板上的顆粒板法直接測定比較困難，因此本文采用模擬標定方法近似獲得顆粒滾動摩擦因數。

為了降低接觸面對尿素顆粒堆積過程的影響，重點考察顆粒間摩擦特性對休止角的影響，提高顆粒間滾動摩擦因數的測定精度，參考文獻[6]設計了如圖10的顆粒堆積模擬和試驗測定平臺，無底細圓管高250 mm、直徑18 mm，圓筒高40 mm、直徑190 mm，事先將400 g尿素在底部圓筒內鋪均勻，并將表面刮平，再將細圓管中填滿尿素顆粒，豎直置于底部圓盤肥料表面，以0.05 m/s的速度將細圓管緩慢提起，使尿素顆粒在肥料水平面上滾落堆積，從而獲得肥料堆如圖10所示。休止角的測定均由Matlab讀取顆粒堆積圖像，參見上述休止角的測定方法。

圖10 尿素顆粒堆積模擬和試驗平臺Fig.10 Simulation and experimental verification of urea particles’ accumulation

得到尿素顆粒間靜摩擦因數和尿素顆粒與ABS間靜摩擦因數分別為0.32和0.33，其余仿真參數同3.1.2節。在EDEM中進行不同肥料顆粒滾動摩擦因數下的休止角的測定，取6組重復試驗平均值作為最終測定的休止角，得到肥料顆粒間休止角ψ′與滾動摩擦因數μrpp的函數關系如圖11，其擬合方程為

ψ′=62.317μrpp+17.386

(9)

圖11 顆粒間滾動摩擦因數與休止角關系曲線Fig.11 Relationship curves between urea-urea rolling friction coefficient and repose angle

二者呈線性變化趨勢，方程擬合精度較高(R2=0.989 4)，因此可用于標定尿素顆粒間滾動摩擦因數。

利用休止角試驗測定裝置測定尿素顆粒休止角，試驗選取2.1節中經篩分的大顆粒尿素，用烘干法測得含水率為0.11%。實際測定試驗時與EDEM仿真計算保證試驗條件一致，測得真實顆粒的休止角為30.30°。將試驗測定值代入擬合方程(9)，計算得到標定的尿素顆粒間滾動摩擦因數μrpp=0.21。

3.3 休止角試驗與仿真模型的驗證

為驗證尿素顆粒間滾動摩擦因數以及顆粒間、顆粒與ABS間靜摩擦因數3個重要邊界參數確定的可靠性，在三者分別為0.21、0.32和0.33時對尿素顆粒進行休止角的模擬仿真(其余邊界參數同3.1.2節)，休止角仿真結果如圖12a所示，測得休止角為30.32°，在相同試驗條件下測定尿素顆粒實際休止角結果如圖12b所示，測得休止角均值為30.43°，二者相差0.11°，相對誤差為0.36%，模擬結果接近試驗值，表明試驗測定和標定后的EDEM邊界參數和仿真模型的有效性，也證明對于難試驗測定的邊界參數通過離散元模擬進行參數標定的可行性。

圖12 休止角仿真結果與試驗結果對比Fig.12 Comparison of simulation and experiment results

為研究標定的3個重要邊界參數對不同含水率大顆粒尿素的適應性，對不同含水率下尿素顆粒休止角進行了測量[23]。根據GB 2440—2001中農用尿素的含水率標準(小于等于1.0%)，試驗中發現當含水率為1.1%時尿素顆粒開始出現粘結現象，因此文中研究的尿素顆粒最大含水率為1.2%。不同含水率下的休止角與本節中標定參數下仿真測定休止角間的相對誤差如表4所示，可以看出在含水率小于1.2%時相對誤差均低于3.25%，表明通過仿真得到的3個重要參數標定值適用于不同含水率下尿素顆粒的標定，因此仿真模型和顆粒間滾動摩擦因數、靜摩擦因數、顆粒與ABS間靜摩擦因數的標定值0.21、0.32和0.33可用于不同含水率下尿素顆粒的仿真研究。

表4 不同含水率尿素顆粒的休止角測定結果Tab.4 Testing results of urea particles repose angle with different moisture contents

4 結論

(1)設計了Plackett-Burman休止角仿真試驗篩選出影響顯著的邊界參數重要性順序依次為尿素顆粒間的滾動摩擦因數、顆粒間的靜摩擦因數、顆粒與ABS板間的靜摩擦因數，且單因素仿真試驗表明，3種邊界參數影響尿素顆粒堆積特性，休止角隨著3種邊界參數的增大而增大。

(2)基于離散元軟件對尿素顆粒間靜摩擦因數和滾動摩擦因數進行虛擬仿真標定，并進行堆積過程的仿真和試驗測定，仿真休止角與實際試驗休止角相對誤差為0.36%，表明標定后的邊界參數和仿真模型的有效性，為難以試驗測定的邊界參數確定提供了一種有效的解決途徑。

(3)測量了不同含水率下大顆粒尿素的休止角，與標定參數下仿真測定休止角對比相對誤差均不大于3.25%，證明了仿真模型和標定的邊界參數適合于不同含水率的尿素顆粒標定要求。