顆粒肥料離散元仿真摩擦因數標定方法研究

溫翔宇 袁洪方 王 剛 賈洪雷

(1.吉林大學生物與農業工程學院， 長春 130022； 2.吉林大學工程仿生教育部重點實驗室， 長春 130022)

0 引言

科學施肥、提高肥料的利用率是當前急需解決的問題。為了尋求更加科學的施肥方式，國內外學者在精確施肥、變量施肥的理論上進行了大量研究，形成了相對成熟的體系，肥料不再直接從肥箱排施到土壤中，增加了控制肥料排施過程[1-3]的中間環節。通過試驗可以測定施肥作業效果，但排肥過程卻很難直接觀測。因此，通過基于離散元法的軟件對肥料顆粒建模，模擬肥料排施過程進行數字化設計，探究影響肥料排施效果的因素，可以提高研發效率，節約成本[4]。

在離散元分析軟件EDEM中進行仿真試驗前，需要設定待研究材料的物理特性參數(泊松比、剪切模量、密度)以及各個材料間的接觸參數(碰撞恢復系數、靜摩擦因數和滾動摩擦因數)，這些參數對仿真試驗結果有著十分重要的影響[5-7]。張銳等[8]通過SLB-1型應力應變測試三軸剪切試驗儀，測定沙土的彈性模量，計算出沙土的泊松比；石林榕等[9]通過HDV-1K電動雙柱拉壓力測試儀對種薯樣體進行靜態壓縮試驗，可通過公式計算出種薯的彈性模量、泊松比、剪切模量。材料密度一般可通過天平和量筒應用排液法測定，但各個材料間的接觸參數卻難以通過試驗直接測定，需要在EDEM軟件中進行仿真標定試驗，標定的接觸參數是否準確直接影響仿真試驗結果的可靠性。

GELDART等[10]分析了8種粉末顆粒堆積角的測量方法，不同的測試方法對休止角測量結果有較大影響，這是因為不同測試方法，影響休止角測量結果的顯著因素以及因素顯著程度不同。劉彩玲等[11]利用Plackett-Burman試驗方法進行尿素顆粒的休止角仿真試驗，篩選影響顯著的邊界參數依次為尿素顆粒間滾動摩擦因數、顆粒間靜摩擦因數和顆粒與ABS板間靜摩擦因數，通過自制靜摩擦因數測量儀和虛擬仿真標定方法對顯著因素進行研究，并對確定值進行了堆積過程的仿真和試驗驗證，仿真休止角與實際試驗休止角相對誤差僅為0.36%。袁全春等[12]提出一種通過仿真試驗建立回歸模型并結合物理試驗尋優的方法，對有機肥離散元模型參數進行標定，在標定的參數下進行仿真驗證試驗，仿真休止角與實際休止角的相對誤差為0.42%。

目前，對于顆粒物料的參數標定多數只選用一種方法來標定多個接觸參數，但不同的測試方法所得到的結果也反映了顆粒物料的特性。因此，本文通過Plackett-Burman因素顯著性篩選試驗對多種顆粒特性測試方法進行研究，選出不同方法所對應的顯著因素，進行綜合標定，以實現標定參數可以反映顆粒物料的整體特性，提高離散元仿真試驗結果的精確度。

1 材料與方法

1.1 試驗材料基本參數

試驗用肥料選用氮質量分數46%的大顆粒尿素，因顆粒肥料在實際的流動排施過程中，主要工作環境是在施肥管中，所以試驗用圓筒等與肥料接觸的材料，均選用與施肥管同種材質的軟質PVC材料，進行參數標定試驗，便于為后續模擬肥料在實際工況中的流動情況研究提供數據基礎。通過查閱資料確定所選材料的泊松比、剪切模量的數值范圍，用天平、量筒通過排液法測定所選尿素顆粒的密度，用精度0.01 mm的游標卡尺測定尿素顆粒的粒徑范圍，計算顆粒球形率。材料的基本參數如表1所示。

1.2 試驗方法

肥料顆粒的自然休止角能反映其流動、摩擦等特性[18]，本文采用分體圓筒法、傾斜法和抽板法進行尿素顆粒休止角試驗，應用Design-Expert 8.0.6軟件進行Plackett-Burman多因素顯著性篩選試驗設計，在EDEM離散元仿真軟件中按試驗設計安排仿真試驗，待顆粒群穩定后，測定顆粒休止角，由于每組試驗參數不同，顆粒堆形態并不總是傾斜的直線，若用屏幕量角器等手段測量顆粒堆的休止角會產生很大的誤差，所以對休止角的測量方法為[19-20]：采用Matlab數字圖像處理軟件讀取顆粒堆邊緣圖像(如圖1a所示)，對圖像進行灰度化、二值化處理輸出二值圖(如圖1b所示)，提取二值化圖像邊界輪廓(如圖1c所示)，掃描輪廓圖像上每一個像素，記錄白點的坐標及個數，通過最小二乘法對所記錄的白點進行線性擬合，并用紅色線畫出擬合直線(如圖1d所示)，在Matlab軟件工作區讀取擬合直線的斜率，可求得顆粒堆邊緣輪廓的傾角。對于分體圓筒法，選擇顆粒堆積成的圓錐體相互垂直的4個方向的母線測定其傾角，求平均值即為肥料顆粒的自然休止角，對于傾斜法和抽板法，取顆粒群左、右兩側傾斜的邊緣輪廓測定其傾角，取平均值為肥料顆粒的自然休止角，每組仿真試驗重復3次取平均值。

將測定的休止角數值填入Design-Expert 8.0.6軟件中，進行方差分析，分析3種方法影響休止角的顯著因素，并增加斜面法測定組合顆粒與PVC材料靜摩擦因數仿真試驗，同樣進行Plackett-Burman篩選試驗設計，驗證斜面法標定組合顆粒與PVC材料靜摩擦因數的可行性。通過真實試驗測定3種休止角測試方法對應的尿素顆粒休止角以及斜面法的斜面傾角，每組試驗重復3次取平均值，將其作為目標值，分別對不同方法的顯著因素進行仿真標定試驗，最后用標定的參數進行仿真驗證試驗，與真實試驗進行對比，驗證標定參數的準確性。

2 休止角仿真試驗

2.1 仿真模型建立

真實試驗采用干燥的肥料顆粒作為試驗材料，含水率0.37%，肥料表面幾乎無粘附力，因此仿真采用Hertz-Mindlin無滑動接觸模型，并設置重力加速度方向。

肥料顆粒的粒徑分布對顆粒群體積密度有直接的影響，較小的顆粒可以填充在大顆粒間的空隙中，增大顆粒間的接觸面積，仿真試驗顆粒群粒徑分布與真實情況不一致，會對堆積試驗結果產生影響。因此，通過篩分法統計肥料顆粒的粒徑分布，發現EDEM軟件中提供的固定粒徑、平均分布、正態分布3種粒徑分布方式難以單一地反映本次試驗中肥料樣本真實的粒徑分布，通過預試驗建立兩種尿素顆粒模型，分別命名為D Particle與X Particle，D Particle的顆粒模型為半徑1.5 mm的球形顆粒，在顆粒工廠中設置正態分布的顆粒平均直徑為肥料顆粒模型直徑的1倍，標準差0.28 mm，最小顆粒直徑為平均直徑的1倍，最大顆粒直徑為平均值的1.666倍，X Particle的顆粒模型為半徑0.8 mm的球形顆粒，在顆粒工廠中設置正態分布的顆粒平均直徑為肥料顆粒模型直徑的1倍，標準差0.80 mm，最小顆粒直徑為平均直徑的1倍，最大顆粒直徑為平均值的1.875倍，仿真顆粒粒徑分布與篩分法得到的肥料粒徑分布接近。

按照不同的試驗方法分別設置對應的仿真幾何模型，分體圓筒法仿真設置如圖2所示，無底圓筒與底座直徑均為100 mm，無底圓筒高180 mm，底座高40 mm。初始狀態時，無底圓筒與底座貼合，裝置內部生成尿素顆粒，待大顆粒尿素在重力作用下沉降后，無底圓筒以0.001 m/s的速度沿豎直方向提升，逐步與底座分離，失去圓筒壁面支撐的尿素顆粒滑落至底座外部，在EDEM仿真試驗中設置計算域，不對滑落至底座外部的顆粒進行計算，以減少仿真工作量，提高仿真速度。待尿素顆粒在底座穩定堆積后，測定尿素顆粒休止角。

圖2 分體圓筒法休止角仿真

傾斜法仿真設置如圖3所示，圓筒直徑為100 mm，高度180 mm，將圓筒上下兩端封閉并水平放置，圓筒容器生成占其容積1/2～1/3的尿素顆粒，設置圓筒繞水平軸以2(°)/s的速度緩慢旋轉，當顆粒的表面產生滑動時，測定其表面的傾斜角。

圖3 傾斜法休止角仿真

抽板法仿真設置如圖4所示，長方體盒子在XY面投影為邊長80 mm的正方形，長方體高度150 mm，初始狀態長方體內部封閉并生成尿素顆粒，待顆粒在重力作用下沉降后，長方體一側的活動板以0.002 m/s的速度沿豎直方向提升，失去壁面支撐的顆粒滑出長方體，待穩定后測定顆粒堆休止角。

圖4 抽板法休止角仿真

斜面法仿真設置如圖5所示，為了保證尿素顆粒在斜面上滑動而非滾動，采用4個半徑2 mm的圓球建立組合顆粒模型，通過調整圓球坐標位置，設置圓球球心均在同一水平面上且4個圓球間兩兩相互接觸(如圖5a所示)，設置生成顆粒數量為1，斜面初始角度為0，斜面繞其一側的水平軸以1(°)/s的速度緩慢旋轉(如圖5b所示)，待組合顆粒在傾斜的斜面上產生滑動趨勢時，記錄此時斜面的傾斜角度。

圖5 斜面法仿真

2.2 仿真試驗

考慮到現階段軟質PVC材料因加工工藝以及增塑劑的成分和添加比例的不同，會對材料的本征參數產生影響，因此將軟質PVC材料的泊松比、剪切模量也作為影響因素，共選取10個因素進行顯著性篩選試驗，篩選影響顯著的因素，試驗因素與水平見表2。

表2 Plackett-Burman試驗因素水平

應用Design-Expert 8.0.6軟件進行Plackett-Burman試驗設計，試驗方案如表3所示，a～k表示各試驗因素對應水平值，L列作為空白列，將每組試驗對應的仿真參數輸入EDEM軟件中，進行休止角仿真試驗。

表3 Plackett-Burman試驗方案

分體圓筒法、傾斜法、抽板法、斜面法4種顆粒特性測試方法的仿真結果見表4。

表4 Plackett-Burman試驗結果

對分體圓筒法、傾斜法、抽板法和斜面法測定的結果進行方差分析，得出所選的10個因素對休止角影響的顯著程度，方差分析結果如表5～8所示。

分體圓筒法方差分析結果表明，尿素顆粒間靜摩擦因數對分體圓筒法測定的休止角結果影響極顯著，尿素顆粒間滾動摩擦因數對休止角影響顯著，其他因素對休止角影響不顯著，分體圓筒法與無底圓筒法的差別在于：無底圓筒法錐形顆粒堆底面與試驗材料接觸，顆粒與試驗材料間的摩擦因數對休止角產生影響，而分體圓筒法中與錐形顆粒堆下方接觸的是顆粒本身，所以在無底圓筒法試驗中顆粒間靜摩擦因數、顆粒間滾動摩擦因數和顆粒與試驗材料間靜摩擦因數對休止角影響顯著[21]，而分體圓筒法中僅尿素顆粒間靜摩擦因數和尿素顆粒間滾動摩擦因數對休止角影響顯著。

表5 分體圓筒法試驗方差分析結果

注：*表示顯著(p<0.05)，** 表示極顯著(p<0.01)。下同。

表6 傾斜法試驗方差分析結果

表7 抽板法試驗方差分析結果

傾斜法方差分析結果表明，尿素顆粒間靜摩擦因數對傾斜法測定的休止角結果影響顯著，尿素顆粒與PVC材料間靜摩擦因數對休止角影響極顯著，其他因素對休止角影響不顯著，仿真過程中顆粒在重力作用下相互擠壓，邊緣的顆粒被擠壓在PVC圓筒內壁面上，在顆粒與PVC材料間靜摩擦力和顆粒間靜摩擦力的作用下，顆粒群隨著圓筒繞水平軸轉動，處于相對靜止狀態，直至顆粒群表面與水平面傾角達到一定值時，顆粒間以及顆粒與PVC材料間的力學平衡被打破，顆粒群表面開始滑動，使休止角始終在小范圍內波動，維持穩定值。

從抽板法方差分析數據可以看出，尿素顆粒間靜摩擦因數和尿素顆粒與PVC材料間靜摩擦因數的正負效應、均方和與其他因素有顯著差異，這與傾斜法分析的結果相似，但p仍未達到顯著要求，而且抽板法試驗模型也不顯著。

表8 斜面法試驗方差分析結果

斜面法方差分析表明，僅顆粒與PVC材料間靜摩擦因數對試驗結果影響極顯著，其余因素對試驗結果不顯著，驗證了斜面法對組合顆粒與PVC間靜摩擦因數標定的可行性。

3 摩擦因數標定

根據上述仿真試驗分析結果，提出一種基于顆粒物料整體特性的摩擦因數標定方法，將仿真試驗與真實試驗相結合，采用斜面法、傾斜法、分體圓筒法對尿素顆粒間靜摩擦因數、尿素顆粒間滾動摩擦因數、尿素顆粒與PVC材料間靜摩擦因數這3個因素進行標定，其余7個因素參數不變：尿素泊松比0.25，尿素剪切模量3.4×107Pa，軟質PVC泊松比0.47，軟質PVC剪切模量2×106Pa，尿素顆粒間碰撞恢復系數0.28，尿素顆粒與PVC材料間碰撞恢復系數0.35，尿素顆粒與PVC材料間滾動摩擦因數0.04，而尿素顆粒間靜摩擦因數作為不變參數時，設定數值為0.50，尿素顆粒間滾動摩擦因數作為不變參數時，設定數值為0.08。

3.1 尿素顆粒與PVC間靜摩擦因數標定

在EDEM軟件中對尿素顆粒與PVC間靜摩擦因數進行單因素試驗，因素范圍設置在0.1～0.6之間，間隔為0.1，將試驗結果繪制成散點圖，發現尿素顆粒與PVC間靜摩擦因數在0.1～0.6之間，與實際情況接近，并且線性相關性極好，在此范圍內對斜面傾角α和尿素顆粒與PVC間靜摩擦因數μp進行線性擬合(如圖6所示)，擬合方程為

α=52.5μp+1.4(0.1≤μp≤0.6) (R2=1)

(1)

圖6 尿素顆粒與PVC間靜摩擦因數擬合曲線

采用精度0.1 mg、滿量程110 g的精密電子天平，稱取4粒等效直徑4 mm、質量近似相等的尿素顆粒，在同一平面上將其粘結在一起，將組合顆粒置于粘貼有PVC材料的斜面上(如圖7所示)，斜面初始角度0，緩慢提升斜面角度，待組合顆粒在PVC斜面上開始下滑時，記錄此時斜面的傾角為22.9°，將其代入擬合方程(1)中，可求出尿素顆粒與PVC間靜摩擦因數μp為0.41。

圖7 斜面法試驗

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

在傾斜法因素顯著性篩選試驗中僅顆粒間靜摩擦因數和顆粒與PVC材料間靜摩擦因數影響顯著，因此將斜面法標定的顆粒與PVC材料間靜摩擦因數0.41作為固定值，對尿素顆粒間靜摩擦因數進行單因素仿真試驗，因素范圍設置在0.05～0.40之間，間隔為0.05，發現尿素顆粒間靜摩擦因數在0.05～0.40范圍內休止角趨近于二次項擬合(如圖8所示)，顆粒群表面傾角β和尿素顆粒間靜摩擦因數μn二次項擬合方程為

(2)

圖8 尿素顆粒間靜摩擦因數擬合曲線

真實試驗用圓筒采用軟質PVC材料圍制而成，與仿真試驗中圓筒模型尺寸一致，考慮到軟質PVC材質的圓筒易變形，在PVC圓筒外套裝一個尺寸相匹配的樹脂圓筒，使PVC圓筒在填充物料以及試驗過程中不會發生形變(如圖9所示)，在圓筒內填充其容積1/3的尿素顆粒后，將圓筒用PVC材料封閉，使圓筒水平放置并以穩定的角速度繞中心軸緩慢轉動，待顆粒堆表面開始滑動時，記錄此時傾角為35.97°，將其代入擬合方程(2)中，求解得μn為0.36和0.90，因一元二次方程的解0.90不在本試驗范圍，所以得出尿素顆粒間靜摩擦因數μn為0.36。

圖9 傾斜法休止角測定試驗

3.3 顆粒間滾動摩擦因數標定

在分體圓筒法因素顯著性篩選試驗中顆粒間靜摩擦因數和顆粒間滾動摩擦因數影響顯著，將尿素顆粒間靜摩擦因數0.36和顆粒與PVC材料間靜摩擦因數0.41作為固定值，對尿素顆粒間滾動摩擦因數進行單因素仿真試驗，因素范圍設置在0.01～0.30之間，間隔為0.05，尿素顆粒休止角γ和尿素顆粒間滾動摩擦因數μs擬合曲線如圖10所示，擬合方程為

(3)

圖10 尿素顆粒間滾動摩擦因數擬合曲線

圖11 分體圓筒法休止角測定試驗

真實試驗裝置如圖11所示，將尿素顆粒填充在圓筒內部，緩慢提升PVC無底圓筒，使其與PVC底座逐漸分離，部分尿素顆粒滑落至PVC底座外，待底座上方的尿素顆粒堆穩定后，選擇圓錐體相互垂直的4個方向的母線與水平面的夾角，求平均值即為對應肥料顆粒的自然休止角，試驗重復3次取平均值，為37.26°，將其代入擬合方程(3)中，求出尿素顆粒間滾動摩擦因數μs為0.15。

3.4 驗證試驗

圖12 仿真結果與真實試驗結果對比

將所標定的顆粒與PVC材料間靜摩擦因數0.41、顆粒間靜摩擦因數0.36和顆粒間滾動摩擦因數0.15采用無底圓筒法進行仿真試驗，PVC無底圓筒直徑為80 mm，高180 mm，圓筒內部生成尿素顆粒，待大顆粒尿素在重力作用下沉降后，PVC無底圓筒以0.05 m/s的速度沿豎直方向提升，待顆粒堆穩定后測定其休止角，試驗重復3次取平均值，休止角仿真結果為30.57°(如圖12a所示)，在相同試驗條件下測定大顆粒尿素的實際休止角為31.74°(如圖12b所示)，誤差為3.69%，無顯著差異，試驗結果證明了所標定因數的有效性，標定方法滿足離散元仿真要求。

考慮到肥料在存放或使用的過程中易受潮，導致肥料含水率會發生變化，影響肥料顆粒的物理特性，本文采用無底圓筒法測定了含水率在0.37%～1.13%之間的尿素顆粒休止角，與本文所標定參數條件下的仿真試驗結果進行對比，相對誤差如表9所示。

表9 不同含水率的尿素顆粒休止角

試驗結果表明：仿真結果與真實試驗的相對誤差不大于4.59%，所標定摩擦因數滿足不同含水率下尿素顆粒的仿真研究。

4 結論

(1)對分體圓筒法、傾斜法、抽板法和斜面法4種顆粒特性測試方法進行研究，采用Plackett-Burman多因素顯著性篩選試驗設計，對離散元仿真軟件EDEM中需要設定的物理特性參數(泊松比、剪切模量、密度)以及各個材料間的接觸參數(碰撞恢復系數、靜摩擦因數和滾動摩擦因數)共10個因素進行顯著性方差分析。分析結果表明：尿素顆粒間靜摩擦因數對分體圓筒法測定的休止角結果影響極顯著，尿素顆粒間滾動摩擦因數影響顯著；對于傾斜法測定的休止角結果，尿素顆粒間靜摩擦因數影響顯著，尿素顆粒與PVC材料間靜摩擦因數影響極顯著；斜面法方差分析表明，僅組合顆粒與PVC材料間靜摩擦因數對試驗結果影響極顯著。

(2)基于上述分析結果，提出一種基于顆粒物料整體特性的摩擦因數標定方法，通過斜面法對顆粒與PVC材料間靜摩擦因數進行單因素仿真試驗，得出試驗結果的線性擬合方程，將真實試驗結果斜面傾角22.9°代入擬合方程，標定出顆粒與PVC材料間靜摩擦因數為0.41，將其作為固定值，采用傾斜法對顆粒間靜摩擦因數進行單因素試驗，與真實試驗結合，標定顆粒間靜摩擦因數為0.36，用分體圓筒法對顆粒間滾動摩擦因數進行擬合，標定為0.15。

(3)將標定摩擦因數采用無底圓筒法進行驗證試驗，仿真試驗結果為30.57°，真實試驗結果為31.74°，誤差為3.69%，無顯著差異。在含水率0.37%～1.13%條件下測定尿素顆粒的休止角，與本文所標定參數條件下的仿真試驗結果進行對比，相對誤差不大于4.59%，試驗結果驗證了所標定摩擦因數的有效性，標定方法滿足離散元仿真要求。