基于JKR粘結模型的蚯蚓糞基質離散元法參數標定

羅　帥　袁巧霞　GOUDA Shaban　楊龍元

(1.華中農業大學工學院， 武漢 430070； 2.農業部長江中下游農業裝備重點實驗室， 武漢 430070；3.本哈大學農學院， 本哈 13736)

0　引言

蚯蚓堆肥是處理畜禽糞便的有效方式[1]，經過蚯蚓過腹處理的禽畜糞便具有良好的孔隙特性和豐富的營養物質，成為優良的作物育苗和栽培基質，即蚯蚓糞基質。蚯蚓堆肥在禽畜糞便的綜合利用中效果顯著，逐漸成為研究的熱點[2]。但目前蚯蚓堆肥過程中從布料、取料、蚓糞分離到蚯蚓糞基質的應用，基本上靠人工進行，大量人工成本已成為制約蚯蚓堆肥的主要因素，因此需研究蚯蚓堆肥處理及蚓糞利用相關設備。而研發其設備，就應該對蚯蚓糞基質相關特性進行充分了解。

蚯蚓糞基質是一種典型的散體物料，離散元法(Discrete element method, DEM)被廣泛用于這些散體物料特性的研究[3-5]。JKR顆粒粘結模型是由JOHNSON等[6]于1971年提出并以JOHNSON、KENDALL和ROBERTS等的名字命名的模型，該模型引入了顆粒間表面能(Surface energy)的概念，適用于模擬細小和潮濕材料顆粒間的黏聚作用。JKR粘結模型自被提出以來，已經被許多學者引入到不同散體顆粒間的相互作用模型中，如黃豆和槭樹殘渣的混合攪拌運動的數值模擬[7]、型砂的流動性仿真模擬[8]、土壤堆積的數值模擬[9]、煤粉和生物質混合物休止角的推導[10]等，這些研究模擬結果和試驗結果均較為吻合。

休止角(Angle of repose，AoR)作為散體物料的一種固有屬性，與散體物料顆粒本身的各種屬性息息相關，常被用于散體物料顆粒參數標定[11-14]。蚯蚓糞基質含水率變化時，其物理參數也會發生變化，這些變化可以通過休止角的變化來體現。

本文在對國內外散體物料參數標定和休止角測定方法進行總結的基礎上，結合蚯蚓糞基質堆積體實際輪廓，提出一種休止角測定方法，基于JKR粘結模型進行蚯蚓糞基質的離散元法參數標定，設計休止角測定儀并通過試驗測定不同含水率下蚯蚓糞基質的休止角，得到蚯蚓糞基質顆粒含水率與休止角之間的關系。通過測定蚯蚓糞基質含水率來預測其休止角，繼而合理推定蚯蚓糞基質顆粒間滾動摩擦因數等其他參數。

1　材料與方法

1.1　試驗材料

試驗中所用的蚯蚓糞基質原料取自武漢市東西湖區新溝鎮蚯蚓養殖基地，以牛糞為主要原料，經大平2號蚯蚓過腹處理，充分腐熟后得到可以用于植物育苗的蚯蚓糞基質原料。基質顆粒近似為球體。考慮蚯蚓糞基質實際應用時的處理方法，試驗時將蚯蚓糞基質進行粉碎處理，全部通過孔徑為2 mm的標準篩(浙江上虞市五四儀器廠生產)，以滿足育苗缽成型試驗粒徑要求。蚯蚓糞基質原料pH值為5.89，電導率為0.22 mS/cm。

在不同含水率下休止角的測定試驗中，將蚯蚓糞基質在遮蔭條件下自然風干，使蚯蚓糞基質含水率緩慢下降，不定時取樣分別測定休止角與含水率。

1.2　試驗設備

參考目前廣泛使用的休止角測定方法，設計了圖1所示的鋼質休止角測試儀。落料漏斗由高度調節螺桿上的螺母限制與堆積底座間的高度差。試驗中設置漏斗出料口與堆積底座上表面距離為150 mm。堆積底座上端為直徑150 mm的圓板。試驗時，漏斗中的蚯蚓糞基質顆粒經漏斗口落于堆積底座上，在側面對堆積體進行拍照，對照片進行處理以獲取蚯蚓糞基質顆粒的休止角。

圖1　休止角測試儀Fig.1　AoR determination instrument1.落料漏斗　2.高度調節螺桿　3.堆積底座

1.3　離散元模型

為縮短數值模擬時間，需要對模型進行適當的簡化。將休止角測定儀僅保留漏斗口和堆積底座，得到其簡化模型(圖2a)。經觀察，蚯蚓糞基質顆粒大部分為近似的球體(圖2b)，因此以半徑為1 mm的球體為顆粒原型(圖2c)，同時在生成顆粒時，將原型顆粒半徑變化范圍設置為滿足平均值，標準差為0.1 mm的標準正態分布，得到蚯蚓糞基質顆粒的離散元簡化模型。

圖2　測定裝置及顆粒離散元模型Fig.2　Measuring device and DEM model of particles

綜合考慮數值模擬的效率和準確性，時間步長通常設置為雷利時間步長的5%～40%，本文在模擬中視雷利時間步長的具體值，設置時間步長為雷利時間步長的5%～15%。為保證仿真結果的準確性，網格尺寸設置為最小顆粒半徑的2倍。顆粒生成后經漏斗口落于堆積底座上自然堆積，仿真時間統一設置為3 s。

JKR顆粒粘結模型將表面能引入顆粒間的相互作用，其簡化模型如圖3所示。

圖3　JKR粘結模型Fig.3　JKR contact model

1.4　休止角測定方法

休止角測定較廣泛的方法主要有兩種。一種方法是通過顆粒體高度H與底面直徑D來計算休止角，即休止角

(1)

此種方法目前應用最為廣泛[15-17]，但由于顆粒位置的不規則性，堆積體底面邊緣較分散，不是連續的圓面，不便于度量；且堆積體頂端往往也不是規則的錐形，堆積體高度H難于確定。有研究者[18]提出將包含95%顆粒數目的圓作為堆積體底面邊界圓，以該圓直徑作為休止角底面直徑。該方法在數值模擬中可行，但實際操作中難以找到邊界圓。

另一種方法是選取堆積體輪廓中較為平直的一段曲線，以與該段曲線最接近的直線傾角作為該顆粒堆積體休止角。目前此方法主要用數字化圖像分析(Digital image analysis，DIA)來實現，王云霞等[19]在對玉米種子休止角提取時用直線擬合堆積體輪廓線。但是由于所取曲線段的位置因人而異，此種方法結果受人為影響較大。FRACZEK等[20]在研究中亦指出了該方法的這一弊端。

觀察蚯蚓基質堆積體和文獻中其他散體物料堆積體的實際形狀，發現大多數散體物料堆積體兩端輪廓線近似為水平距離的凹函數，居中輪廓線近似為凸函數，而坡體中間段近似為直線。據此，擬采用高斯分布擬合堆積體輪廓線來獲取堆積體休止角。主要思路是通過圖像處理提取堆積體的輪廓點坐標后用高斯分布進行擬合，以擬合曲線拐點處的切線與水平軸的夾角為堆積體的休止角。其具體方法是，通過相機對顆粒堆積體進行拍照，得到圖4a；通過Photoshop軟件快速選擇提取堆積體輪廓并對圖片角度進行校正，得到圖4b；通過Matlab軟件將所獲得的圖像依次進行灰度處理(圖4c)、二值化處理(圖4d)，再通過Photoshop軟件描邊工具提取輪廓曲線(圖4e)；最后將圖像導入Origin軟件通過圖片文件數字化工具Digitizer獲取輪廓點坐標，用高斯分布對其進行擬合(圖4f)。圖像處理的思路參考了FADAVI等[21]種子輪廓處理的方法。

圖4　蚯蚓糞基質休止角的計算過程Fig.4　Calculation process of AoR of vermicomposting nursery substrate

高斯分布的一般方程為

(2)

式中y0、A、w、xc——常數

本文定義擬合方程拐點處切線與x軸所夾銳角為該散體的休止角。

對式(2)進行求導，分別求其一階和二階導數

(3)

(4)

曲線的拐點出現在二階導數為零處，于是令f″(x)=0，式(4)有2個解

(5)

將式(5)代入式(3)，即可得到拐點處的斜率

(6)

則休止角的計算公式為

θ=arctan|Sl|

(7)

休止角測定試驗在華中農業大學工科基地溫室大棚進行。試驗中，基質在遮蔭下自然風干。不定時取樣測定休止角，每次測定重復3次；同時另取樣測定含水率。

1.5　參數標定試驗設計

針對蚯蚓糞基質的研究目前尚不完善，其離散元模型的參數尤其缺乏。考慮到蚯蚓糞基質物理性質與土壤較為接近，且含水率變化時其泊松比、剪切模量等參數均會發生變化，本文主要參照文獻[22-29]中各類土壤的參數值確定試驗中各參數的取值或范圍。表1給出了各待標定參數的高低水平。

其他參數[27]為鋼泊松比0.3、鋼剪切模量7.9×1010Pa、鋼的密度7 865 kg/m3、重力加速度9.81 m/s2。

本文試驗參數較多，參照文獻[30-31]的試驗設計，先進行Plackett-Burman試驗，篩選出對結果影響顯著的試驗因素，再進行Box-Behnken試驗，得到休止角和顯著性參數之間的回歸模型。通過測定不同含水率下蚯蚓糞基質的休止角，得到休止角隨含水率變化的關系曲線。

表1　離散元法參數標定試驗標定參數Tab.1　DEM calibration parameters need to be calibrated in experiment

2　結果與分析

2.1　Plackett-Burman試驗

表2是使用Minitab設計的Plackett-Burman試驗及試驗所得休止角。

表3給出了由Plackett-Burman試驗得出的參數顯著性分析數據。根據試驗結果，顆粒間靜摩擦因數、顆粒間滾動摩擦因數、JKR表面能對休止角的影響是顯著的，其他參數影響較小。

2.2　Box-Behnken試驗

根據Plackett-Burman試驗的結果，進行Box-Behnken試驗，以得到休止角與3個顯著參數之間的關系模型。

表4是Minitab設計的Box-Behnken試驗及試驗所得休止角。重點考察顆粒-顆粒靜摩擦因數(A)、顆粒-顆粒滾動摩擦因數(B)和JKR表面能(C)這3個對休止角影響顯著的參數。使用Minitab軟件建立休止角與3個顯著參數間的二次多項式回歸方程

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

表3　Plackett-Burman試驗參數顯著性分析Tab.3　Analysis of significance of parameters in Plackett-Burman test

θ=10.39+3.66A+115.01B+30.447C-6.34A2-
56.31B2-3.477C2+2.17AB+1.632AC-22.25BC

(8)

表4　Box-Behnken試驗設計及結果Tab.4　Design and results of Box-Behnken test

為0.998 1，均接近于1。表5是該模型方差分析的結果。模型的P值小于0.000 1，極為顯著；失擬項的P值大于0.05，對結果不顯著。模型擬合較好，可靠度較高，可以用來預測休止角。

表5　Box-Behnken試驗休止角二次多項式回歸模型方差分析Tab.5　Box-Behnken test AoR corner quadratic polynomial regression model variance analysis

為驗證模型的有效性，在休止角測定試驗中隨機選取了2組試驗結果(休止角分別為43.05°和36.38°(圖5))，應用Minitab的響應優化器工具，在試驗參數的取值范圍內以休止角分別對回歸模型進行尋找最優解，使用最優解進行離散元模擬，模擬所得休止角分別為43.71°和36.46°，均接近于實際休止角，差異分別為1.53%和0.22%。模擬結果與試驗結果差別均較小，認為模型是有效的。最優解分別為顆粒間靜摩擦因數為1.00，顆粒間滾動摩擦因數為0.11，JKR表面能為0.89 J/m2；顆粒間靜摩擦因數為0.65，顆粒間滾動摩擦因數為0.12，JKR表面能為0.5 J/m2。

圖5　蚯蚓糞基質堆積體仿真結果與實際堆積體圖像Fig.5　Simulation results of vermicomposting nursery substrate accumulation and image of actual accumulation

2.3　含水率與休止角關系試驗

不同含水率下基質的休止角變化如圖6所示。蚯蚓糞基質的休止角隨含水率的降低而降低，這與其他物料[17，32]的研究結果相似。目前，研究者[33-34]一般將散體含水率和休止角之間的關系用線性來描述。但是由蚯蚓糞基質休止角隨含水率變化的散點圖可以發現，隨含水率降低，休止角減小的速率有放慢的趨勢。本文分別用直線和指數函數對散點圖進行擬合，得到擬合方程

θl=19.848 02Cm+30.808 68(R2=0.91)

(9)

(10)

指數方程式(10)擬合度更優，更符合兩者之間的變化趨勢。對于蚯蚓糞基質的休止角，可以測定其含水率，根據式(10)計算預測得到其休止角。

圖6　蚯蚓糞基質含水率和休止角的關系曲線Fig.6　Relationship curves between vermicomposting nursery substrate moisture content and AoR

3　結論

(1)提出了一種散體顆粒休止角獲取方法，即提取堆積體輪廓，用高斯分布對輪廓曲線進行擬合，以擬合曲線拐點的傾角作為休止角。試驗表明，該方法擬合效果良好，為堆積體輪廓與蚯蚓糞基質相似的散體的休止角的測定提供了一種思路。

(2)將JKR粘結模型用于蚯蚓糞基質顆粒，采用離散元法對蚯蚓糞基質顆粒參數進行標定試驗，篩選出對休止角影響顯著的參數(即顆粒間靜摩擦因數、顆粒間滾動摩擦因數、JKR表面能)，在此基礎上對休止角試驗數據進行二次多項式回歸分析，建立了休止角與3個顯著參數間的回歸模型。經試驗驗證，模型結果與試驗結果較為吻合。

(3)測定不同含水率下蚯蚓糞基質顆粒的休止角，建立了含水率與休止角之間的指數函數，對比發現較傳統的直線擬合效果更優。由此函數及休止角與3個顯著參數之間的回歸模型，可以通過測定蚯蚓糞基質含水率、預測休止角，繼而合理推定顆粒間靜摩擦因數、顆粒間滾動摩擦因數、JKR表面能等其他參數。

1SANGWAN P, KAUSHIK C P, GARG V K. Vermicomposting of sugar industry waste (press mud) mixed with cow dung employing an epigeic earthwormEiseniafetida[J]. Waste Management & Research, 2010, 28(1): 71-75.

2HUANG K, XIA H. Role of earthworms’ mucus in vermicomposting system: biodegradation tests based on humification and microbial activity[J]. Science of the Total Environment, 2018, 610: 703-708.

3ANAND A, CURTIS J S, WASSGREN C R, et al. Experimental study of wet cohesive particles discharging from a rectangular hopper[J]. Industrial & Engineering Chemistry Research, 2015, 54(16): 4545-4551.

4BRIEND R, RADZISZEWSKI P, PASINI D. Virtual soil calibration for wheel-soil interaction simulations using the discrete-element method[J]. Canadian Aeronautics and Space Journal, 2011, 57(1): 59-64.

5劉凡一,張艦,李博, 等. 基于堆積試驗的小麥離散元參數分析及標定[J]. 農業工程學報,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)

6JOHNSON K L, KENDALL K, ROBERTS A D. Surface energy and the contact of elastic solids[J]. Proceedings of the Royal Society of London Series A, Mathematical and Physical Sciences, 1971, 324(1558):301-313.

7CUNHA 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.

8張帥, 單忠德, 張杰. 基于離散元方法的型砂流動性仿真研究[J]. 鑄造技術, 2016, 37(2):288-291.

ZHANG Shuai, SHAN Zhongde, ZHANG Jie. Simulation of self-hardening resin sand mobility based on DEM[J]. Foundry Technology, 2016, 37(2):288-291. (in Chinese)

9胡紅. 玉米行間定點扎穴深施追肥機設計與研究[D]. 北京：中國農業大學, 2017.

HU Hong. Design and research of targeted hole-pricking and deep-application fertilizer between maize rows [D]. Beijing: China Agricultural University, 2017. (in Chinese)

10GUO Z, CHEN X, LIU H, et al. Theoretical and experimental investigation on angle of repose of biomass-coal blends[J]. Fuel, 2014, 116(1): 131-139.

11GHAZAVI M, HOSSEINI M, MOLLANOURI M. A comparison between angle of repose and friction angle of sand[C]∥The 12th International Conference for International Association for Computer Methods and Advances in Geomechanics (IACMAG), 2008.

12FIELKE J M, UCGUL M, SAUNDERS C. Discrete element modeling of soil-implement interaction considering soil plasticity, cohesion and adhesion[C]∥2013 ASABE International Meeting, 2013.

13FRIEDMAN S P, ROBINSON D A. Particle shape characterization using angle of repose measurements for predicting the effective permittivity and electrical conductivity of saturated granular media[J]. Water Resources Research, 2002, 38(11): 11-18.

14XIAO X, TAN Y, ZHANG H, et al. Experimental and DEM studies on the particle mixing performance in rotating drums: effect of area ratio[J]. Powder Technology, 2017, 314: 182-194.

15PROBST K, AMBROSE K, ILELEJI K. The effect of moisture content on the grinding performance of corn and corncobs by hammer milling[J]. Transactions of the ASABE, 2013, 56(3):1025-1033.

16UNAL H, ISIK E, IZLI N, et al. Geometric and mechanical properties of mung bean (VignaradiataL.) grain: effect of moisture[J]. International Journal of Food Properties, 2008, 11(3): 585-599.

17ZAALOUK A K, ZABADY F I. Effect of moisture content on angle of repose and friction coefficient of wheat grain[J]. Misr Journal of Agricultural Engineering, 2009, 26(1): 418-427.

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

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):173-179. (in Chinese)

19王云霞, 梁志杰, 張東興, 等. 基于離散元的玉米種子顆粒模型種間接觸參數標定[J]. 農業工程學報,2016,32(22):36-42.

WANG Yunxia, LIANG Zhijie, ZHANG Dongxing, et al. Calibration method of contact characteristic parameters for corn seeds based on EDEM[J]. Transactions of the CSAE, 2016, 32(22): 36-42. (in Chinese)

21FADAVI A, MIRZABE A H, MANSOURI A. Moisture-dependent physical properties of Plantain (PlantagomajorL.) seeds by image processing analysis[J]. Agricultural Engineering International: CIGR Journal, 2015, 17(3): 353-363.

22張銳, 韓佃雷, 吉巧麗, 等. 離散元模擬中沙土參數標定方法研究[J/OL]. 農業機械學報, 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)

23武濤, 黃偉鳳, 陳學深,等. 考慮顆粒間黏結力的黏性土壤離散元模型參數標定[J]. 華南農業大學學報, 2017, 38(3):93-98.

WU Tao, HUANG Weifeng, CHEN Xueshen, et al. Calibration of discrete element model parameters for cohesive soil considering the cohesion between particles[J]. Journal of South China Agricultural University, 2017, 38(3):93-98. (in Chinese)

24潘世強. 基于離散元法的芯鏵式開溝器優化設計與試驗研究[D]. 長春：吉林大學, 2015.

PAN Shiqiang. Research on the optimization design and the experiment of the core ploughshare furrow opener based on the discrete element method [D]. Changchun: Jilin University, 2015. (in Chinese)

25方會敏, 姬長英, AHMED A T, 等. 秸稈-土壤-旋耕刀系統中秸稈位移仿真分析[J/OL]. 農業機械學報, 2016, 47(1):60-67. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160109&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2016.01.009.

FANG Huimin, JI Changying, AHMED A T, et al. Simulation analysis of straw movement in straw-soil-rotary blade system[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(1):60-67. (in Chinese)

26UCGUL M, FIELKE J M, SAUNDERS C. 3D DEM tillage simulation: validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil[J]. Soil and Tillage Research, 2014, 144(4): 220-227.

27UCGUL M, FIELKE J M, SAUNDERS C. Three-dimensional discrete element modelling (DEM) of tillage: accounting for soil cohesion and adhesion[J]. Biosystems Engineering, 2015, 129: 298-306.

28邵明安，王全九，黃明斌. 土壤物理學[M]. 北京：高等教育出版社， 2006.

29依艷麗. 土壤物理研究法[M]. 北京：北京大學出版社, 2009.

30ABDENACER M, KAHINA B I, A?CHA N, et al. Sequential optimization approach for enhanced production of glutamic acid fromCorynebacteriumglutamicum2262 using date juice[J]. Biotechnology and Bioprocess Engineering, 2012, 17(4): 795-803.

31JOB J, SUKUMARAN R K, JAYACHANDRAN K. Production of a highly glucose tolerant β-glucosidase byPaecilomycesvariotiiMG3: optimization of fermentation conditions using Plackett-Burman and Box-Behnken experimental designs[J]. World Journal of Microbiology and Biotechnology, 2010, 26(8): 1385-1391.

32KALKAN F, KARA M. Handling, frictional and technological properties of wheat as affected by moisture content and cultivar[J]. Powder Technology, 2011, 213(1): 116-122.

33UNAL H, ISIK E, IZLI N, et al. Geometric and mechanical properties of mung bean (VignaradiataL.) grain: effect of moisture[J]. International Journal of Food Properties, 2008, 11(3): 585-599.

34SOLOGUBIK C A, CAMPAONE L A, PAGANO A M, et al. Effect of moisture content on some physical properties of barley[J]. Industrial Crops and Products, 2013, 43(5): 762-767.