不同含水率羊糞離散元參數通用標定方法研究

朱新華 伏勝康 李旭東 魏玉強 趙 偉

(1.西北農林科技大學機械與電子工程學院， 陜西楊凌 712100； 2.四川航天職業技術學院， 成都 610100)

0 引言

糞肥施肥體量大、工費高，而高效的施肥機械較為缺乏。目前，借助離散元法分析肥料顆粒間及肥料與裝置間復雜的相互作用已成為施肥機械設計、優化的重要手段，可節約研發成本、縮短周期，提高裝置性能[1-4]。離散元法的應用以準確的離散元參數為前提[5-6]，快速準確標定離散元仿真參數是其在農業工程領域應用的關鍵[7]。羊糞是我國北方來源廣、抗生素等污染小的有機肥源，其有機質含量高達40%～65%[8]，高于牛糞、馬糞、豬糞和雞鴨鵝等糞肥[9]，在實際生產中受到重視，但腐熟羊糞的離散元參數標定尚未見報道。當前，糞肥離散元參數主要采用虛擬標定試驗獲得[10-11]，但這些研究所構建的標定方法存在通用性和實用性不足的問題。

國外學者以液態有機肥流體特性的研究為主，而對固態粘濕有機肥離散元參數標定方法的研究尚未見報道[7]。國內學者基于堆積角試驗對有機肥離散元參數標定的研究較多[12-13]，這些研究表明運用堆積角試驗標定糞肥離散元參數能得到較為準確的結果。但上述研究均針對特定含水率下的糞肥進行研究，而實際糞肥的含水率范圍較大，因此，該方法的通用性不足[7]。物料的離散元參數受含水率影響較大，特定含水率下標定的離散元參數不能真實反映其它含水率的糞肥與機械之間的作用規律。

基于堆積角試驗，羅帥等[14]對不同含水率的蚯蚓糞基質離散元參數進行標定，建立了堆積角與離散元參數的關系模型；在含水率40%～85%范圍內，王黎明等[15]探索了通過堆積角快速標定豬糞離散元參數的方法。盡管上述研究考慮了含水率對糞肥離散元參數的影響，提高了堆積角-離散元參數模型和標定方法的通用性，但存在以下不足：糞肥堆積角測量的偶然誤差大，必須多次測量，且糞肥堆積操作不便，以堆積角-離散元參數模型預測糞肥離散元參數，實用性不足；研究中僅建立單一的肥料顆粒模型，與實際肥料的粒度分布差異較大[7]，且對所建模型的準確性缺少工程性驗證(即用另外一種試驗方法驗證)，所得離散元參數可靠性不足。

因此，針對不同含水率下畜禽糞肥的離散元參數快速標定方法通用性和實用性不足的問題，本文以含水率為31.91%～58.48%的腐熟羊糞為研究對象，通過圓筒提升法[16]建立糞肥的含水率-堆積角模型；以堆積角為響應值，通過Plackett-Burman試驗、爬坡試驗和Box-Behnken試驗建立堆積角-離散元參數模型；進而經推導建立羊糞含水率-離散元參數模型，并運用抽板法[16]進行工程性驗證，以期構建一種通用性和實用性強的羊糞離散元參數標定方法。

1 材料與方法

1.1 試驗材料與粒度分布

腐熟羊糞取自陜西省楊凌區西北農林科技大學第三試驗站，原始含水率46.81%，密度0.47 g/cm3。在腐熟糞肥通常含水率范圍(30%～60%)[17-19]內通過自然風干法和添加純凈水的方法配制5個含水率梯度，計算方法為[14]

(1)

式中ms——配制樣品的加水質量，kg

m0——樣品質量，kg

ω0——原始含水率，%

ω1——配制樣品的目標含水率，%

每次堆積角測量的同時，均取一定量羊糞樣品通過烘干法測量其含水率，測得供堆積角試驗的5個含水率依次為：31.91%、39.33%、46.81%、53.32%和58.48%，如圖1所示。取300 g原始含水率下的羊糞樣本，利用振動篩(ZFJ-Ⅱ型，浙江省上虞市大亨橋化驗儀器廠)對其篩分，分析其粒徑分布(表1)。

表1 羊糞粒徑分布Tab.1 Particle size distribution of goat manure

圖1 不同含水率的羊糞試樣Fig.1 Samples of goat manure with different moisture contents

1.2 羊糞堆積角試驗及含水率-堆積角模型構建

采用圓筒提升法測定羊糞堆積角[16]。圖2為試驗裝置，由電動式拉力試驗機(HSV型)、圓筒(內徑100 mm、長300 mm)和平臺組成[16,20]。圓筒和平臺均為鋼材，圓筒直立于平臺上。試驗前在圓筒內填滿羊糞，試驗時通過拉力試驗機以20 mm/s的速度提升圓筒[16]，使羊糞自圓筒下端漏出，形成肥堆，如圖3所示。用數顯傾角儀從肥堆的4個不同方向測量堆積角。以相同的方法依次對5種含水率羊糞進行堆積角測量，每組試驗重復5次，結果取平均值。

圖2 圓筒提升法堆積角測試裝置Fig.2 Cylinder lifting resting angle tester1.拉力試驗機 2.圓筒 3.平臺

圖3 堆積角測量Fig.3 Measurement of resting angle

基于測量結果，建立羊糞含水率-堆積角模型，并分析其相關性。

1.3 羊糞離散元模型構建

(1)接觸模型。羊糞流動性差，顆粒間存在粘附力，故采用Hertz-Mindlin with JKR模型[21-22](簡稱JKR模型)。該接觸模型內置于EDEM(2018，EDEM Solutions)軟件中，是一種能體現顆粒粘結性的接觸模型，模型基于Johnson-Kendall-Roberts理論建立，在Hertz-Mindlin接觸模型的基礎上考慮了接觸區域內范德華力對顆粒運動的影響，適用于模擬干粉狀物料或含濕顆粒，并通過數值surface energy表征顆粒間的粘附力，可較好地模擬羊糞。

(2)顆粒粒徑分布及材料本征參數。根據測定的羊糞粒徑分布，同時考慮合理簡化模型，在EDEM仿真中以球形顆粒模擬羊糞，顆粒粒徑設置4種，分布如圖4所示。

圖4 羊糞顆粒模型Fig.4 Particle models of goat manure

羊糞和鋼材的本征參數如表2所示。

表2 羊糞和鋼材的本征參數Tab.2 Intrinsic parameters of goat manure and steel material

(3)接觸參數和JKR模型參數。結合EDEM通用顆粒材料數據庫(Generic EDEM material model database，GEMM數據庫)，根據羊糞、鋼材的本征參數及文獻[10-11,13-14,16]的推薦值，確定接觸參數及JKR模型參數取值范圍如表3所示。

表3 羊糞與鋼材接觸參數和JKR模型參數范圍Tab.3 Contact parameters and JKR model parameter range of goat manure and steel material

1.4 羊糞離散元參數標定及堆積角-離散元參數模型構建

1.4.1羊糞堆積角仿真方法

仿真過程設置與物理試驗一致。設置圓筒提升速度20 mm/s，生成羊糞顆?？傎|量3 kg，生成速度0.2 kg/s，仿真固定時間步長為Rayleigh時間步長的20%，為7.4×10-6s，數據保存時間間隔為0.05 s，網格大小為1.5 mm，仿真總時長為15 s。仿真幾何模型及堆積角測量如圖5所示，運用EDEM后處理Clipping模塊在坐標為(1，0，0)和(0，1，0)羊糞堆中垂面處設置Clipping Planes Group截面，運用Tools模塊中的量角器Protractor，依次選定Clipping Planes Group截面上羊糞堆單側的斜邊、夾角和水平底面線上的3個羊糞顆粒作為堆積角R1的3個點，獲得堆積角數值。

圖5 圓筒提升法仿真Fig.5 Simulation of cylinder lifting method

1.4.2羊糞堆積角-離散元參數模型構建

應用Design-Expert 8.0(STAT-EASE)軟件設計Plackett-Burman(PBD)試驗，以堆積角為響應值，從羊糞的密度、泊松比、剪切模量以及表3所示的接觸參數和JKR模型參數共10個參數中篩選對堆積角影響顯著的離散元參數。各參數取高、中、低3個水平，分別以編碼1、0和-1表示，如表4所示，試驗方案及結果見2.2節。

表4 Plackett-Burman試驗因素編碼Tab.4 Factors and levels of Plackett-Burman test

為縮小參數范圍并確定最優區間，針對顯著性參數設計爬坡試驗。試驗時，非顯著性參數取中間水平值，顯著性參數按照設計的步長逐步增大，將仿真堆積角與中間水平含水率(46.81%)羊糞物理試驗堆積角(41.00°)相對誤差最小時的參數作為顯著性參數最優區間的中間水平值，以確定最優取值區間。爬坡試驗方案及結果見2.3節。

基于爬坡試驗確定的離散元參數最優區間設計Box-Behnken試驗，并對試驗結果進行模型優化，以建立堆積角與顯著性參數的關系模型。采用3個中心點進行誤差估計，共15組試驗，每組試驗重復3次，取堆積角平均值作為試驗結果。試驗中各顯著性參數取爬坡試驗優化后的高、中、低3個水平，分別以編碼1、0和-1表示，非顯著性參數取值同爬坡試驗一致，試驗方案及結果見2.4節。

以5個含水率下羊糞物理試驗堆積角為目標，對Box-Behnken試驗得到的堆積角-離散元參數模型求解尋優得到離散元參數的最優組合?；陔x散元參數最優組合(非顯著性參數取值同爬坡試驗一致)，運用EDEM進行仿真并測量其堆積角，對比仿真堆積角與物理試驗堆積角以驗證堆積角-離散元參數模型的準確性。

1.5 羊糞含水率-離散元參數模型建立及工程性驗證試驗

通過堆積角-離散元參數模型能夠推測離散元參數，但糞肥堆積角測量誤差大、過程繁瑣，而含水率測量精度高、簡便。因此，基于羊糞含水率-堆積角模型和羊糞堆積角-離散元參數模型推導，建立羊糞含水率-離散元參數模型，以構建通用性和實用性好、準確度高的離散元參數標定方法。

為評價基于圓筒提升法建立的含水率-離散元參數模型的可靠性，通過另外一種堆積角試驗方法——抽板法[17, 20]對模型進行工程性驗證。驗證過程中，仿真與物理試驗方法一致，圖6為基于抽板法的堆積角仿真與物理驗證試驗。

圖6 基于抽板法的堆積角仿真與物理驗證試驗Fig.6 Resting angle simulation and physical verification test based on side plate lifting method

在含水率31.91%～58.48%范圍內隨機配制3個不同含水率的羊糞，并分別測量其物理試驗堆積角；根據推導出的含水率-離散元參數模型求解，得到3個含水率下的羊糞離散元參數最優組合；以離散元參數最優組合進行仿真，分別得到仿真堆積角。對比物理試驗堆積角與仿真堆積角以驗證羊糞含水率-離散元參數模型的準確性。每組試驗重復5次，結果取平均值。

2 結果與分析

2.1 含水率-堆積角模型構建物理試驗測定的不同含水率羊糞樣本的堆積角如表5所示。

表5 不同含水率羊糞的堆積角Tab.5 Resting angle of goat manure with different moisture contents

根據表5，以羊糞樣本的含水率與堆積角擬合多項式，得到含水率-堆積角模型

R1=0.005 5x2-0.202 2x+38.329

(2)

式中x——羊糞含水率，%

模型相關系數達0.999 9，擬合曲線如圖7所示。在含水率31.91%～58.48%范圍內，羊糞堆積角隨含水率升高而增大。羊糞屬于散粒體物料，隨著含水率的增大，羊糞顆粒間的粘度增大，流動性變差，故其堆積角變大。

圖7 羊糞含水率與堆積角關系曲線Fig.7 Relationship curve between moisture content and resting angle of goat manure

2.2 顯著性參數篩選結果與分析

基于Design-Expert 8.0軟件設計PBD試驗，試驗方案與結果如表6所示，A～J為因素編碼值。

表6 Plackett-Burman試驗方案與結果Tab.6 Scheme and results of Plackett-Burman test

表7為PBD試驗結果顯著性分析。

表7 Plackett-Burman試驗結果顯著性分析Tab.7 Significance analysis of Plackett-Burman test results

根據表7，JKR表面能J(P<0.01)對羊糞堆積角影響極顯著，顆粒間滾動摩擦因數F和靜摩擦因數E(P<0.05)對羊糞堆積角影響顯著，且3個參數對堆積角的影響均為正效應。分析其原因為：①隨著JKR表面能增大，羊糞顆粒間的粘附性能增強，從而影響羊糞顆粒的流動性，因此，堆積角增大。②由于顆粒流動表現為顆粒與顆粒間的滾動和滑動，因此，影響顆粒流動性的主要參數為顆粒間滾動摩擦因數和靜摩擦因數。當滾動摩擦因數和靜摩擦因數增大時，顆粒間的滾動和滑動受阻，流動性變差，堆積角增大。③模型中的羊糞由4種不同粒徑的球型顆粒組成，而球體模型的運動方式以滾動為主，因此，滾動摩擦因數比靜摩擦因數對堆積角的影響更顯著。

為方便后續試驗，在爬坡試驗、Box-Behnken試驗和驗證試驗中只考慮3個顯著性參數，不顯著因素取中間水平值[13]。

2.3 離散元參數最優區間確定

爬坡試驗方案與結果如表8所示，可知，當JKR表面能為0.05 J/m2、羊糞顆粒間滾動摩擦因數為0.15、靜摩擦因數為0.30時，仿真堆積角與物理試驗堆積角相對誤差最小，為2%，即該組參數為顯著性參數最優區間的中間水平值，對應最優區間分別為[0，0.10]、[0，0.30]、[0.10，0.50]。

表8 爬坡試驗方案與結果Tab.8 Scheme and results of climbing test

2.4 堆積角-離散元參數模型建立

表9 Box-Behnken試驗方案與結果Tab.9 Scheme and results of Box-Behnken test

表10 Box-Behnken試驗模型方差分析Tab.10 ANOVA of Box-Behnken model

根據試驗結果及方差分析，E、F、J、E2、F2均對羊糞堆積角影響極顯著，EJ對羊糞堆積角影響顯著。建立的羊糞堆積角-離散元參數模型為

R1=37.04+4.96E+7.50F+11.91J-
2.88EJ+5.51E2+6.82F2

(3)

2.5 堆積角-離散元參數模型驗證

以5個含水率羊糞樣本物理試驗堆積角為目標求解尋優，分別獲得5組離散元參數(JKR表面能、顆粒間滾動摩擦因數和靜摩擦因數)的最優組合，如表11所示。

表11 不同含水率下羊糞離散元參數最優值及驗證試驗結果Tab.11 Optimal value of goat manure discrete element parameters and verification test results under different moisture contents

當含水率為58.48%時，仿真和物理試驗獲得的羊糞堆積角如圖8所示。由表11可知，由堆積角-離散元參數模型確定3個顯著性離散元參數的最優組合后，各含水率羊糞仿真堆積角與物理試驗堆積角的相對誤差小于等于2.42%，表明建立的模型可用于離散元參數確定。

圖8 仿真與物理試驗獲得的堆積角Fig.8 Resting angle obtained by simulation and physical test

2.6 含水率-離散元參數模型構建及驗證

為準確、快速標定離散元參數，基于含水率-堆積角模型和堆積角-離散元參數模型，由式(2)、(3)推導獲得含水率與離散元參數(JKR表面能、顆粒間滾動摩擦因數和靜摩擦因數)的模型，即含水率-離散元參數模型

4.96E+7.50F+11.91J-2.88EJ+5.51E2+
6.82F2=0.005 5x2-0.202 2x+1.289

(4)

根據式(4)，可通過含水率直接得到離散元參數組合目標值。即借助Design-Expert 8.0軟件中的模型目標值尋優求解，得到離散元參數組合。

基于羊糞含水率-離散元參數模型，獲得3組隨機含水率下的羊糞離散元參數最優組合，通過抽板法物理試驗與仿真分別獲得羊糞物理試驗堆積角與仿真堆積角，如圖9所示。每組試驗重復5次，結果取平均值，測量結果如表12所示。

圖9 抽板法試驗與仿真對比Fig.9 Comparisons of test and simulation with side plate lifting method

表12 抽板法驗證試驗結果Tab.12 Verification test results of side plate lifting method

由表12可知，3組試驗中仿真堆積角與物理試驗堆積角的相對誤差小于等于5.37%。表明建立的含水率-離散元參數模型可靠。

羊糞含水率-離散元參數模型的相對誤差略大于堆積角-離散元參數模型，這是由于模型驗證試驗中采用了不同于圓筒提升法的抽板法。抽板法在此作為工程驗證，能進一步反映模型的可靠性。兩個模型的相對誤差均較小，說明建立的兩種模型均可用于離散元參數的預測。

建立的羊糞含水率-堆積角模型與羅帥等[14]研究的蚯蚓糞基質含水率與堆積角關系變化趨勢一致，且篩選得到的顯著性參數相同，分別為JKR表面能、顆粒間滾動摩擦因數和靜摩擦因數。建立的堆積角-離散元參數模型相對誤差小于等于2.42%，均小于含水率25%～65%的蚯蚓糞基質[23]以及含水率40%～85%的豬糞[15]堆積角預測模型誤差。建立的含水率-離散元參數模型通過工程性驗證，相對誤差小于等于5.37%，表明模型可靠。由于含水率測量比堆積角測量準確性高、易操作，考慮到離散元參數預測模型的實際應用，建立的含水率-離散元參數模型比堆積角-離散元參數模型更具實用性。兩種模型均能適應有機肥料不同含水率下離散元參數標定，具有通用性。

3 結論

(1)基于圓筒提升法物理試驗建立了31.91%～58.48%范圍內羊糞含水率-堆積角模型，模型相關系數為0.999 9。

(2)通過Plackett-Burman試驗和爬坡試驗確定了含水率31.91%～58.48%的羊糞離散元顯著性參數及其最優區間：JKR表面能為0～0.10 J/m2、羊糞顆粒間滾動摩擦因數為0～0.30、靜摩擦因數為0.10～0.50；通過Box-Behnken試驗建立了羊糞堆積角-離散元參數模型，模型P<0.000 1，相對誤差小于等于2.42%。

(3)基于含水率-堆積角模型和堆積角-離散元參數模型推導，建立了含水率-離散元參數模型；采用抽板法試驗驗證了羊糞含水率-離散元參數模型的可靠性，相對誤差小于等于5.37%。

(4)基于建立的含水率-離散元參數模型和堆積角-離散元參數模型，均可求解出特定含水率下羊糞離散元顯著性參數的最優組合，對不同含水率羊糞參數標定具有通用性和實用性。