新疆果園深施散體廄肥離散元參數(shù)標(biāo)定研究

韓樹杰 戚江濤 坎 雜 李亞萍 蒙賀偉

(1.石河子大學(xué)機(jī)械電氣工程學(xué)院， 石河子 832003； 2.農(nóng)業(yè)農(nóng)村部西北農(nóng)業(yè)裝備重點(diǎn)實(shí)驗(yàn)室， 石河子 832003)

0 引言

新疆地區(qū)廄肥資源較為豐富，施用廄肥可促進(jìn)果園土壤微生物的生命活動(dòng)，對(duì)提高土壤肥力、改善土壤結(jié)構(gòu)、補(bǔ)充土壤養(yǎng)分和創(chuàng)造良好的土壤生態(tài)環(huán)境具有重要作用[1-3]。

近年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，離散元法(Discrete element method, DEM)在農(nóng)業(yè)裝備研究上應(yīng)用越來越廣泛[4-8]。利用離散元法全面、系統(tǒng)研究散體物料和機(jī)械裝置之間的相互作用機(jī)理和物料的運(yùn)動(dòng)狀態(tài)，不僅可以優(yōu)化機(jī)械裝置的結(jié)構(gòu)參數(shù)和工作參數(shù)，而且可以提高研發(fā)效率、改善機(jī)械作業(yè)性能、節(jié)約成本。物料離散元參數(shù)標(biāo)定是研究物料與機(jī)械裝置之間相互作用的基礎(chǔ)。GRIMA等[9]利用崩塌試驗(yàn)中的顆粒堆休止角對(duì)干顆粒在離散元仿真中所需滾動(dòng)摩擦因數(shù)進(jìn)行了標(biāo)定；BOAC等[10]運(yùn)用離散方法模擬了精選油籽顆粒的材料和相互作用特性；溫翔宇等[11]對(duì)顆粒肥料進(jìn)行離散元仿真，標(biāo)定了顆粒肥的摩擦因數(shù)；袁全春等[12]對(duì)有機(jī)肥散體顆粒離散元模型進(jìn)行參數(shù)標(biāo)定，仿真休止角與實(shí)際休止角的相對(duì)誤差僅為0.42%；羅帥等[13]基于JKR粘結(jié)模型標(biāo)定了蚯蚓糞基質(zhì)的離散元參數(shù)，休止角仿真結(jié)果與實(shí)際試驗(yàn)結(jié)果較為接近；文獻(xiàn)[14-17]分別對(duì)粘性土壤、沙土顆粒土壤的模型參數(shù)進(jìn)行了標(biāo)定，結(jié)果表明，虛擬仿真結(jié)果與實(shí)測(cè)值之間差異較小。

本文在以上研究基礎(chǔ)上，以新疆果園深施散體廄肥為研究對(duì)象，采用仿真試驗(yàn)與物理試驗(yàn)相結(jié)合的方法，對(duì)散體廄肥的離散元參數(shù)進(jìn)行標(biāo)定，從而獲得散體廄肥的離散元模型，以期為全面、系統(tǒng)研究肥料和施肥裝置之間的相互作用機(jī)理以及物料的運(yùn)動(dòng)狀態(tài)提供參考。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)所用廄肥為新疆生產(chǎn)建設(shè)兵團(tuán)第八師一二一團(tuán)葡萄園用廄肥，由純牛糞發(fā)酵，發(fā)酵時(shí)間為3～4個(gè)月。前期試驗(yàn)結(jié)果表明，為保證施肥機(jī)作業(yè)性能，提高施肥機(jī)施肥均勻性和穩(wěn)定性，施肥前需對(duì)廄肥進(jìn)行粉碎處理。用環(huán)刀法測(cè)定所選廄肥的密度；選取500 g廄肥作為試驗(yàn)樣本，使用土壤標(biāo)準(zhǔn)粒徑篩(篩孔直徑0.25～5 mm)確定其粒徑分布情況；用MH45型含水率測(cè)定儀(質(zhì)量測(cè)量精度0.001 g,含水率測(cè)量精度0.01%)測(cè)定廄肥的含水率；堆積角試驗(yàn)中所用鋼板材質(zhì)為45號(hào)鋼，參照文獻(xiàn)[18-20]確定所選材料的泊松比、剪切模量，材料的基本性質(zhì)如表1所示，本征參數(shù)如表2所示。

表1 廄肥基本性質(zhì)Tab.1 Basic properties of manure

表2 廄肥、鋼本征參數(shù)Tab.2 Intrinsic parameters of manure and steel plate

1.2 試驗(yàn)方法

肥料顆粒的自然堆積角能反映其流動(dòng)、摩擦等特性[21],本文采用物理試驗(yàn)與仿真試驗(yàn)相結(jié)合的方法[22-23]對(duì)廄肥離散元模型參數(shù)進(jìn)行標(biāo)定，采用注入法對(duì)廄肥進(jìn)行堆積角試驗(yàn)，并在EDEM離散元仿真軟件中進(jìn)行仿真，應(yīng)用Design-Expert 8.0.6軟件進(jìn)行Plackett-Burman多因素顯著性篩選試驗(yàn)與分析，得出對(duì)堆積角有顯著性影響的參數(shù)；在此基礎(chǔ)上，通過Box-Behnken響應(yīng)面分析法建立并優(yōu)化廄肥堆積角與顯著性參數(shù)的回歸模型，以實(shí)際堆積角為目標(biāo)值對(duì)回歸方程求解尋優(yōu)，得到顯著性參數(shù)最優(yōu)值。最后在最優(yōu)參數(shù)下進(jìn)行仿真試驗(yàn)，對(duì)比廄肥仿真堆積角和實(shí)際堆積角，驗(yàn)證標(biāo)定的廄肥離散元模型參數(shù)的準(zhǔn)確性。

1.3 物理模型

結(jié)合文獻(xiàn)[24-26]對(duì)堆積角的相關(guān)研究，試驗(yàn)參照GB 11986—89/ISO 4324—1977《表面活性劑粉體和顆粒休止角的測(cè)量》，采用注入法測(cè)量散體廄肥的堆積角，測(cè)量裝置如圖1所示，漏斗下口內(nèi)徑為10 mm，錐度為60°，托盤直徑為100 mm，高度為25 mm，漏斗的下端口與托盤上表面距離為75 mm。試驗(yàn)時(shí)，漏斗中的廄肥顆粒經(jīng)漏斗口落于托盤上，最終在托盤上形成穩(wěn)定的顆粒堆，在側(cè)面對(duì)堆積角拍照，采用Matlab對(duì)圖像進(jìn)行處理以獲得廄肥顆粒的堆積角。重復(fù)5次試驗(yàn)取其平均值，試驗(yàn)得到廄肥的堆積角為35.47°。

2 廄肥堆積角仿真試驗(yàn)

2.1 顆粒間接觸模型

在EDEM接觸模型理論中，Hertz-Mindlin with JKR Cohesion模型是一個(gè)凝聚力接觸模型，該接觸模型在Hertz接觸理論的基礎(chǔ)上結(jié)合JKR理論，考慮濕顆粒間粘結(jié)力對(duì)顆粒運(yùn)動(dòng)的影響，適用于模擬顆粒間因水分發(fā)生明顯粘結(jié)和團(tuán)聚的物料[14]。在該模型中，法向彈性接觸力的實(shí)現(xiàn)基于Johnson-Kendall-Roberts理論，切向彈性力、法向耗散力和切向耗散力均與Hertz-Mindlin(no slip)接觸模型中的計(jì)算方法一致，在Johnson-Kendall-Robert理論中，JKR法向彈性力的實(shí)現(xiàn)基于重疊量δ、相互作用參數(shù)和表面能。計(jì)算式為

(1)

(2)

其中

(3)

(4)

式中FJKR——JKR法向彈性力，N

α——相互接觸2個(gè)顆粒的接觸圓半徑，m

γ——表面能，N/m

E*——等效彈性模量，Pa

R*——等效半徑，m

Ei、Ej——相互接觸2個(gè)顆粒的彈性模量，Pa

vi、vj——相互接觸2個(gè)顆粒的泊松比

Ri、Rj——相互接觸2個(gè)顆粒的半徑，m

當(dāng)γ=0時(shí)，力變?yōu)镠ertz-Mindlin法向力

(5)

即使顆粒并不是直接接觸，該模型也提供吸引凝聚力，顆粒間具有非凝聚力的最大間隙為

(6)

(7)

式中δc——顆粒間具有非凝聚力時(shí)的法向最大間隙，m

αc——2個(gè)顆粒的接觸半徑，m

當(dāng)顆粒并非實(shí)際接觸并且間隙小于δc時(shí)，凝聚力達(dá)到最大值

(8)

摩擦力的計(jì)算和Hertz-Mindlin(no slip)接觸模型不同，其取決于JKR法向力的正向排斥部分。因此，該模型在接觸力的凝聚力分量更大時(shí)提供一個(gè)更大的摩擦力。

因此，經(jīng)過處理后的廄肥具有散粒體的物料特性，廄肥顆粒之間受水分子和化學(xué)物質(zhì)的影響也具有粘結(jié)力特性。為了更準(zhǔn)確地模擬廄肥的真實(shí)狀態(tài)，本文選擇Hertz-Mindlin with JKR Cohesion 模型進(jìn)行仿真模擬。

2.2 幾何模型建立

在仿真過程中，顆粒和裝置模型對(duì)仿真結(jié)果有很大的影響，在保證仿真模型尺寸與物理模型一致的基礎(chǔ)上，建立了簡(jiǎn)化的廄肥和漏斗裝置模型。經(jīng)過粉碎處理的散體廄肥顆粒近似球形，將廄肥顆粒模型設(shè)置為球形。在SolidWorks軟件中根據(jù)物理試驗(yàn)?zāi)Ｐ徒⒎抡婺Ｐ?，把模型轉(zhuǎn)換為STL格式，導(dǎo)入EDEM 2.7軟件中作為漏斗裝置仿真模型，本文采用EDEM軟件內(nèi)置的Hertz-Mindlin(no slip)接觸模型作為顆粒與裝置之間的接觸模型，顆粒模型及漏斗模型如圖2、3所示。

2.3 仿真參數(shù)設(shè)置

根據(jù)試驗(yàn)所測(cè)定的廄肥顆粒粒徑分布范圍，為確保仿真結(jié)果準(zhǔn)確性的同時(shí)提高仿真效率，僅生成粒徑分布占比較大區(qū)間的顆粒。在EDEM中采用隨機(jī)分布，將生成的球顆粒半徑限制在0.5～1.25倍的初始球半徑之間。仿真中，動(dòng)態(tài)生成漏斗中顆粒，在漏斗正上方建立顆粒工廠，設(shè)置為虛擬，生成總質(zhì)量為0.6 kg，生成速率為0.2 kg/s，數(shù)據(jù)保存時(shí)間間隔為0.01 s，固定時(shí)間步長(zhǎng)是瑞利時(shí)間步長(zhǎng)的20%，網(wǎng)格尺寸取2倍最小球形單元尺寸?；贓DEM內(nèi)嵌Hertz-Mindlin with JKR Cohesion模型進(jìn)行廄肥接觸參數(shù)的仿真標(biāo)定，其模型參數(shù)(JKR表面能)是表征所研究物料含水率效果的重要參數(shù)。通過大量預(yù)試驗(yàn)確定了JKR表面能的取值范圍。綜合對(duì)比文獻(xiàn)[12,18-20]中所研究肥料、土壤與本文所研究廄肥特性的差異，確定廄肥仿真接觸參數(shù)的取值范圍，如表3所示。

表3 仿真參數(shù)取值范圍Tab.3 Simulation parameters

2.4 仿真試驗(yàn)與結(jié)果分析

2.4.1Plackett-Burman篩選試驗(yàn)

Plackett-Burman篩選試驗(yàn)通過考察目標(biāo)響應(yīng)與各因子間的關(guān)系，比較各個(gè)因子2水平間的差異性來確定因子顯著性。本文Plackett-Burman試驗(yàn)設(shè)計(jì)以廄肥堆積角為響應(yīng)值，對(duì)仿真接觸參數(shù)及接觸模型參數(shù)進(jìn)行篩選。試驗(yàn)接觸參數(shù)高水平設(shè)置為低水平的2倍，根據(jù)文獻(xiàn)[13-14]中因素高低水平取值方法，確定本文因素水平如表4所示。

表4 因素水平Tab.4 Factors and levels

Plackett-Burman試驗(yàn)方案及結(jié)果如表5所示，A～G為因素水平值，H～L為空白列，利用Design-Expert 8.0.6軟件[27-28]對(duì)該結(jié)果進(jìn)行方差分析，得到接觸參數(shù)和接觸模型參數(shù)顯著性如表6所示。由表6可知，JKR表面能、廄肥-廄肥恢復(fù)系數(shù)、廄肥-鋼恢復(fù)系數(shù)的P<0.05，對(duì)廄肥堆積角的影響顯著；而其他參數(shù)的P>0.05，對(duì)廄肥堆積角影響不顯著。

表5 Plackett-Burman試驗(yàn)方案及結(jié)果Tab.5 Scheme and results of Plackett-Burman test

表6 Plackett-Burman試驗(yàn)結(jié)果顯著性分析Tab.6 Significance analysis of Plackett-Burman test results

為方便后續(xù)試驗(yàn)，在Box-Behnken試驗(yàn)[29-30]中只考慮3個(gè)影響顯著的參數(shù)，不顯著因素取值分別為廄肥-廄肥靜摩擦因數(shù)0.65、廄肥-廄肥滾動(dòng)摩擦因數(shù)0.3、廄肥-鋼靜摩擦因數(shù)0.53、廄肥-鋼滾動(dòng)摩擦因數(shù)0.3，進(jìn)行響應(yīng)面試驗(yàn)設(shè)計(jì)。

2.4.2Box-Behnken響應(yīng)面試驗(yàn)及回歸模型

根據(jù)響應(yīng)面設(shè)計(jì)原理，選取顯著性參數(shù)的低、中、高3個(gè)水平進(jìn)行試驗(yàn)設(shè)計(jì)，試驗(yàn)選取5個(gè)中心點(diǎn)對(duì)誤差進(jìn)行評(píng)估。Box-Behnken試驗(yàn)參數(shù)取值如表7所示。

表7 Box-Behnken試驗(yàn)參數(shù)取值Tab.7 Parameter value of Box-Behnken test

Box-Behnken試驗(yàn)方案及結(jié)果如表8所示，重點(diǎn)考察3個(gè)對(duì)廄肥堆積角影響顯著參數(shù)，即JKR表面能、廄肥-廄肥恢復(fù)系數(shù)、廄肥-鋼恢復(fù)系數(shù)。應(yīng)用Design-Expert軟件建立堆積角θ與3個(gè)顯著性參數(shù)的二階回歸方程為

θ=30.92+1.75A+0.88D+3.23G-0.21AD- 1.07AG+0.46DG-0.63A2+0.041D2+0.69G2

(9)

表8 Box-Behnken試驗(yàn)方案及結(jié)果Tab.8 Design and results of Box-Behnken test

Box-Behnken試驗(yàn)?zāi)Ｐ头讲罘治鼋Y(jié)果如表9所示，由表9可知，該擬合模型P=0.000 3，擬合度較好；廄肥-廄肥恢復(fù)系數(shù)(A)、JKR表面能(G)P值均小于0.01；廄肥-鋼恢復(fù)系數(shù)(D)、廄肥-廄肥恢復(fù)系數(shù)和JKR表面能交互項(xiàng)(AG)P值均小于0.05，說明各個(gè)參數(shù)對(duì)堆積角的影響顯著，表明回歸模型的有效性。失擬項(xiàng)P=0.114 7>0.05，表明所得回歸方程與實(shí)際擬合中非正常誤差所占比例小，沒有彎曲失擬現(xiàn)象發(fā)生，擬合性較好。試驗(yàn)中變異系數(shù)為2.65%，說明試驗(yàn)有較高的可靠性。決定系數(shù)為

表9 Box-Behnken試驗(yàn)?zāi)Ｐ头讲罘治鯰ab.9 ANOVA of quadratic polynomial model of Box-Behnken test

θ=30.97+1.75A+0.88D+3.23G-1.07AG

(10)

表10 Box-Behnken試驗(yàn)優(yōu)化回歸模型方差分析Tab.10 ANOVA of modified model of Box-Behnken test

2.4.3回歸模型交互效應(yīng)分析

根據(jù)優(yōu)化回歸模型方差分析結(jié)果，可知廄肥-廄肥恢復(fù)系數(shù)和JKR表面能的交互項(xiàng)(AG)對(duì)廄肥的堆積角影響顯著(P<0.05)。當(dāng)廄肥-廄肥恢復(fù)系數(shù)為0.35時(shí)，應(yīng)用Design-Expert軟件繪制廄肥-廄肥恢復(fù)系數(shù)和JKR表面能交互作用(AG)的響應(yīng)曲面(圖4)，可以直觀地看出兩個(gè)參數(shù)之間的交互效應(yīng)。由AG曲面可知，隨著兩個(gè)參數(shù)取值的增加，廄肥堆積角均呈現(xiàn)上升趨勢(shì)，相對(duì)于廄肥-廄肥恢復(fù)系數(shù)(A)，JKR表面能(G)的效應(yīng)面曲線比較陡，表明其對(duì)堆積角影響較為顯著。

3 最佳參數(shù)組合確定與仿真驗(yàn)證

應(yīng)用Design-Expert軟件對(duì)優(yōu)化后的回歸方程進(jìn)行尋優(yōu)求解，當(dāng)JKR表面能為0.02 J/m2，廄肥-鋼恢復(fù)系數(shù)為0.49，廄肥-廄肥恢復(fù)系數(shù)為0.34，其余非顯著性參數(shù)選取中間水平(廄肥-廄肥靜摩擦因數(shù)0.65，廄肥-廄肥滾動(dòng)摩擦因數(shù)0.3，廄肥-鋼靜摩擦因數(shù)0.53，廄肥-鋼滾動(dòng)摩擦因數(shù)0.3)時(shí)，仿真結(jié)果與實(shí)際堆積角相對(duì)誤差最小。為驗(yàn)證最優(yōu)參數(shù)組合的準(zhǔn)確性，采用上述參數(shù)值，其他參數(shù)設(shè)置不變，應(yīng)用EDEM 2.7軟件進(jìn)行堆積角仿真試驗(yàn)，3次重復(fù)仿真所得廄肥堆積角分別為34.8°、35.0°、33.7°。廄肥堆積角3次平均值為34.5°。與廄肥實(shí)際堆積角35.47°的相對(duì)誤差為2.73%，并且，從圖5a可以直觀看出，利用標(biāo)定后的廄肥參數(shù)得到的堆積邊界與廄肥物理試驗(yàn)堆積結(jié)果比較接近，說明所得廄肥參數(shù)的最優(yōu)值準(zhǔn)確可靠，仿真試驗(yàn)與物理試驗(yàn)的對(duì)照如圖5b、5c所示。

4 結(jié)論

(1)通過Plackett-Burman篩選試驗(yàn)，得到對(duì)廄肥堆積角具有顯著影響的接觸參數(shù)和接觸模型參數(shù)為廄肥-廄肥恢復(fù)系數(shù)、廄肥-鋼恢復(fù)系數(shù)、JKR表面能，而廄肥-廄肥的靜摩擦因數(shù)和滾動(dòng)摩擦因數(shù)、廄肥-鋼的靜摩擦因數(shù)和滾動(dòng)摩擦因數(shù)對(duì)堆積角無顯著性影響。

(2)通過Box-Behnken響應(yīng)曲面試驗(yàn)，得出對(duì)廄肥堆積角影響顯著的參數(shù)，建立顯著性參數(shù)與堆積角之間的二次回歸模型，并對(duì)其進(jìn)行優(yōu)化，根據(jù)其方差分析得出，3個(gè)顯著性參數(shù)的一次項(xiàng)(廄肥-廄肥恢復(fù)系數(shù)、廄肥-鋼恢復(fù)系數(shù)、JKR表面能、廄肥-廄肥恢復(fù)系數(shù)和廄肥表面能的交互項(xiàng)對(duì)廄肥堆積角影響顯著。

(3)通過對(duì)優(yōu)化后的回歸模型求解可知，當(dāng)JKR表面能為0.02 J/m2、廄肥-鋼恢復(fù)系數(shù)為0.49、廄肥-廄肥恢復(fù)系數(shù)為0.34，非顯著性參數(shù)選取中間水平時(shí)，仿真試驗(yàn)和物理試驗(yàn)得到的廄肥堆積角無顯著性差異(P>0.05)，說明采用響應(yīng)面方法分析標(biāo)定廄肥接觸參數(shù)和接觸模型參數(shù)可行。