田辛亮 叢 旭 齊江濤 郭 慧 李 茂 范旭輝
(1.吉林大學(xué)工程仿生教育部重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)春 130022; 2.吉林大學(xué)生物與農(nóng)業(yè)工程學(xué)院, 長(zhǎng)春 130022;3.吉林省農(nóng)業(yè)機(jī)械研究院, 長(zhǎng)春 130022)
保護(hù)性耕作技術(shù)是以機(jī)械化為手段,利用秸稈或殘茬覆蓋于地表,對(duì)于提高土壤有機(jī)質(zhì)含量及降低土壤風(fēng)蝕水蝕具有重要作用[1-2]。玉米秸稈是一種重要的資源,對(duì)提高土壤肥力具有重要作用[3]。秸稈還田是一項(xiàng)有效的保護(hù)性耕作措施,可提高土壤的蓄水保墑能力,改善土壤結(jié)構(gòu),降低秸稈焚燒對(duì)環(huán)境的危害[4-5]。采用松耙聯(lián)合整地機(jī)將玉米秸稈切斷,部分玉米秸稈混入耕層,可達(dá)到加速秸稈腐爛、培肥地力的目的[6]。玉米秸稈-土壤混料特性復(fù)雜,采用常規(guī)方法難以準(zhǔn)確獲得相關(guān)接觸參數(shù),故利用離散元仿真模擬的方法研究土壤及玉米秸稈的相互作用規(guī)律,對(duì)秸稈還田及黑土區(qū)土壤保護(hù)具有重要意義[7]。
隨著信息技術(shù)的發(fā)展,離散元法在農(nóng)業(yè)裝備研究中應(yīng)用越來(lái)越廣泛[8-9]。構(gòu)建離散元仿真模型,本征參數(shù)可通過(guò)試驗(yàn)測(cè)出,由于顆粒間接觸力學(xué)特性復(fù)雜,部分接觸物理參數(shù)難以通過(guò)常規(guī)試驗(yàn)獲取,因此近幾年國(guó)內(nèi)外學(xué)者針對(duì)部分接觸物理參數(shù)難以準(zhǔn)確獲取的問(wèn)題,提出了基于離散元法對(duì)顆粒物料進(jìn)行虛擬標(biāo)定的方法[10-11]。目前,研究人員針對(duì)土壤[12-14]、種子[15-16]、秸稈[17]、肥料[18-19]等利用離散元模型進(jìn)行了仿真參數(shù)標(biāo)定。文獻(xiàn)[20]提出了一種通過(guò)試驗(yàn)與模擬相結(jié)合系統(tǒng)地標(biāo)定沙土顆粒相互作用參數(shù)的方法;文獻(xiàn)[21]利用試驗(yàn)測(cè)定及仿真模擬相結(jié)合的方法對(duì)微型薯顆粒離散元參數(shù)進(jìn)行了標(biāo)定和校準(zhǔn);文獻(xiàn)[16]結(jié)合臺(tái)架試驗(yàn)和仿真試驗(yàn),對(duì)三七種子進(jìn)行仿真參數(shù)標(biāo)定,并選取三七精密排種器進(jìn)行了驗(yàn)證試驗(yàn);文獻(xiàn)[17]利用離散元軟件對(duì)玉米秸稈的接觸物理參數(shù)進(jìn)行仿真標(biāo)定;文獻(xiàn)[18]對(duì)4種不同測(cè)試方法進(jìn)行篩選試驗(yàn),提出了一種基于顆粒物料整體特性的摩擦因數(shù)標(biāo)定方法;文獻(xiàn)[22]將非線性彈性和線性滯回彈簧(塑性)接觸模型的模擬結(jié)果與掃描式耕作工具的試驗(yàn)結(jié)果進(jìn)行比較,得到了土壤接觸物理參數(shù)的變化規(guī)律。綜上所述,國(guó)內(nèi)外研究人員的離散元仿真參數(shù)研究對(duì)象主要是單獨(dú)的土壤、種子、秸稈、肥料等,但對(duì)玉米秸稈-土壤混料的研究卻鮮有報(bào)道。黑土區(qū)保護(hù)性耕作模式中玉米秸稈-土壤混料沒(méi)有明顯規(guī)律可尋,且有限元玉米秸稈-土壤混料模型不能模擬土壤及秸稈的運(yùn)動(dòng)過(guò)程,因此有必要對(duì)其展開(kāi)研究。
本文以黑土區(qū)保護(hù)性耕作模式中玉米秸稈-土壤混料為研究對(duì)象,結(jié)合物理試驗(yàn)與仿真試驗(yàn),應(yīng)用EDEM軟件對(duì)玉米秸稈-土壤混料的接觸參數(shù)進(jìn)行標(biāo)定。以堆積角為響應(yīng)值,采用圖像處理方法測(cè)量混料堆積角,依次通過(guò)Plackett-Burman篩選試驗(yàn)、爬坡試驗(yàn)和Box-Behnken試驗(yàn)對(duì)玉米秸稈-土壤混料堆積角進(jìn)行離散元仿真和校準(zhǔn),以確定較優(yōu)的玉米秸稈-土壤混料離散元仿真參數(shù)值。
本文研究對(duì)象為玉米秸稈-土壤混料,來(lái)源于長(zhǎng)春市農(nóng)安縣開(kāi)安鎮(zhèn)試驗(yàn)基地(44°10′N,120°10′E),該基地種植作物為玉米,土壤類型為黑鈣土,土壤pH值為7.01,年降水量535 mm。試驗(yàn)樣品取自2020年秋季松耙聯(lián)合整地機(jī)作業(yè)后,經(jīng)田間取樣檢測(cè),玉米秸稈在土壤中的混入深度主要為0~20 cm,故本文取樣深度為0~20 cm。取樣測(cè)量利用S型法按照“隨機(jī)、等量和多點(diǎn)”的原則進(jìn)行[23],選取500 g玉米秸稈-土壤混料為樣品,將樣品中的土壤與秸稈進(jìn)行分離,分別測(cè)量土壤與秸稈的密度;依據(jù)機(jī)具作業(yè)后秸稈實(shí)際形狀,秸稈密度按照文獻(xiàn)[17]的方法進(jìn)行測(cè)量。采用干燥法測(cè)量土壤與秸稈的含水率;將采集土樣進(jìn)行干燥,使用不同孔徑的細(xì)篩進(jìn)行篩分,結(jié)果表明土壤粒徑分布范圍為2~4 mm。樣品的基本性質(zhì)如表1所示。
表1 樣品基本性質(zhì)Tab.1 Basic properties of sample
運(yùn)用堆積角物理試驗(yàn)與EDEM仿真試驗(yàn)相結(jié)合的方法對(duì)玉米秸稈-土壤混料的離散元模型參數(shù)進(jìn)行標(biāo)定,采用圓筒提升法[24]進(jìn)行玉米秸稈-土壤混料堆積角試驗(yàn),測(cè)量玉米秸稈-土壤混料堆積角的實(shí)際均值?;贓DEM軟件進(jìn)行仿真試驗(yàn),應(yīng)用Design-Expert軟件進(jìn)行Plackett-Burman篩選試驗(yàn)設(shè)計(jì),并進(jìn)行多因素顯著性分析篩選出對(duì)玉米秸稈-土壤混料有顯著影響的參數(shù);利用最陡爬坡試驗(yàn)確定顯著性影響參數(shù)最優(yōu)區(qū)間,根據(jù)Box-Behnken試驗(yàn)建立并優(yōu)化玉米秸稈-土壤混料堆積角與顯著性參數(shù)的回歸模型,并對(duì)參數(shù)進(jìn)行綜合優(yōu)化得到最優(yōu)參數(shù)組合。最后利用最優(yōu)參數(shù)組合進(jìn)行仿真試驗(yàn),對(duì)比玉米秸稈-土壤混料仿真堆積角與實(shí)際堆積角的差異,驗(yàn)證標(biāo)定模型參數(shù)的準(zhǔn)確性。
玉米秸稈-土壤混料堆積角物理試驗(yàn)裝置如圖1所示,由WDW-20型萬(wàn)能試驗(yàn)機(jī)、鋼制圓筒、鋼板組成。鋼制圓筒高度與直徑比為3∶1[25-26],圓筒高度為210 mm,直徑為70 mm,圓筒底部與鋼板(300 mm×300 mm)接觸,圓筒內(nèi)裝填玉米秸稈-土壤混料。利用WDW-20型萬(wàn)能試驗(yàn)機(jī)以0.05 m/s的速度勻速提升圓筒[27],混料經(jīng)圓筒底部落于鋼板上,待混料堆穩(wěn)定后利用相機(jī)垂直拍照,相機(jī)鏡頭與萬(wàn)能試驗(yàn)機(jī)正面平行,相機(jī)底面與萬(wàn)能試驗(yàn)機(jī)正面垂直,得到混料堆的正視圖像。應(yīng)用Matlab軟件進(jìn)行灰度處理,提取邊緣輪廓并進(jìn)行線性擬合以獲得玉米秸稈-土壤混料的堆積角。重復(fù)試驗(yàn)5次求平均值,得到玉米秸稈-土壤混料的堆積角為40.17°。
玉米秸稈-土壤混料因其本身具有一定的濕度,故顆粒間存在著粘附現(xiàn)象。Hertz-Mindlin with JKR接觸模型是一種凝聚力接觸模型,該模型用表面能代表顆粒間的相互吸引力,適用于模擬濕粘性顆粒和團(tuán)聚的物料[28],模型的法向彈性力的實(shí)現(xiàn)基于Johnson-Kendall-Roberts(JKR)理論,基于法向重疊量、表面能及相互作用,計(jì)算公式為
(1)
(2)
(3)
(4)
式中FJKR——JKR法向彈性力,N
E*——等效彈性模量,Pa
δ——兩接觸顆粒之間法向重疊量,m
α——兩接觸顆粒間的接觸圓半徑,m
γ——表面能,N/m
R*——等效接觸半徑,m
R1、R2——兩顆粒的接觸半徑,m
ν1、ν2——兩接觸顆粒泊松比
E1、E2——兩接觸顆粒彈性模量,Pa
當(dāng)γ為0時(shí),JKR法向力與Hertz-Mindlin法向力FHertz相同,即
(5)
若顆粒不是直接接觸,該模型也可以提供吸引凝聚力,顆粒間的非零凝聚力最大間隙為
(6)
(7)
式中δc——顆粒間有非零凝聚力時(shí)的法向最大間隙,m
αc——顆粒間有非零凝聚力時(shí)的切向最大間隙,m
當(dāng)顆粒并非實(shí)際接觸且間隙小于δc時(shí),凝聚力達(dá)到最大值,即
(8)
當(dāng)選用JKR模型模擬含水率較高的顆粒時(shí),2個(gè)顆粒分開(kāi)所需要的力取決于表面能和接觸角,計(jì)算公式為
(9)
式中θ——兩顆粒的接觸角
玉米秸稈-土壤混料本身含水率相對(duì)較高,混料顆粒間受水分子影響存在著粘附現(xiàn)象,而JKR模型中表面能可較好模擬顆粒之間的粘結(jié)力,故本文選用Hertz-Mindlin with JKR接觸模型,并依此進(jìn)行參數(shù)標(biāo)定。
應(yīng)用離散元仿真軟件及EDEM建立仿真模型,土壤顆粒采用隨機(jī)分布的方法生成[12],以球體作為仿真顆粒模型,使用不同孔徑的細(xì)篩進(jìn)行篩分后確定土壤粒徑為2~4 mm。玉米秸稈形狀近似于長(zhǎng)條狀,松耙聯(lián)合整地機(jī)對(duì)全量秸稈覆蓋的耕地進(jìn)行作業(yè),秸稈在作業(yè)過(guò)程中經(jīng)過(guò)2次切斷,依據(jù)隨機(jī)多點(diǎn)的原則對(duì)作業(yè)后的秸稈長(zhǎng)度進(jìn)行測(cè)量,確定秸稈長(zhǎng)度為40 mm[4]。目前玉米秸稈在離散元軟件中建模時(shí)被認(rèn)為剛性體,通過(guò)球面堆積的方法組合形成玉米秸稈,故離散元模型設(shè)置秸稈長(zhǎng)度為40 mm,由半徑為3 mm的12個(gè)基本球體相連接組合而成。在SolidWorks軟件中建立三維鋼制圓筒模型,把模型另存為IGS格式,導(dǎo)入EDEM軟件中;鋼板模型由EDEM軟件生成,材料屬性設(shè)置為鋼,仿真試驗(yàn)中物料的本征參數(shù)如表2所示。
表2 物料本征參數(shù)Tab.2 Material intrinsic parameters
仿真過(guò)程中,在鋼制圓筒上設(shè)置顆粒工廠,顆粒的生成方式設(shè)置為static,數(shù)據(jù)保存間隔為0.01 s,重力加速度為9.81 m/s2。根據(jù)田間篩分實(shí)測(cè)結(jié)果,土壤與秸稈質(zhì)量比為99∶1。土壤和秸稈按照質(zhì)量比采用分層交錯(cuò)生成的方式,各生成7層。為保證土壤與秸稈混合的均勻性,參照文獻(xiàn)[25]的試驗(yàn)方法,采用旋轉(zhuǎn)圓筒的方法實(shí)現(xiàn)土壤與秸稈的均勻混合。經(jīng)多次重復(fù)試驗(yàn)結(jié)果對(duì)比,顆粒生成后設(shè)置圓筒轉(zhuǎn)速為60 r/min,圓筒內(nèi)物料均勻混合時(shí)間為10 s。待顆粒達(dá)到穩(wěn)定狀態(tài)后,圓筒以0.05 m/s的速度勻速豎直向上運(yùn)動(dòng),顆粒在重力作用下滑落在鋼板上,直至所有顆粒停止運(yùn)動(dòng)并形成穩(wěn)定顆粒堆并測(cè)量堆積角,如圖2、3所示。
在前期反復(fù)試驗(yàn)的基礎(chǔ)上,參照文獻(xiàn)[12,17,21,23,25]的參數(shù)選取范圍,確定玉米秸稈-土壤混料接觸物理參數(shù)范圍如表3所示。
表3 離散元仿真參數(shù)與水平Tab.3 Parameters and level required in DEM simulation
利用Design-Expert軟件進(jìn)行Plackett-Burman篩選試驗(yàn)設(shè)計(jì),以玉米秸稈-土壤混料堆積角為試驗(yàn)指標(biāo),篩選出對(duì)指標(biāo)存在顯著性影響的參數(shù)。篩選試驗(yàn)共計(jì)12組,a~j分別為秸稈-土壤恢復(fù)系數(shù)、秸稈-土壤靜摩擦因數(shù)、秸稈-土壤滾動(dòng)摩擦因數(shù)、土壤-土壤恢復(fù)系數(shù)、土壤-土壤靜摩擦因數(shù)、土壤-土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、土壤-鋼滾動(dòng)摩擦因數(shù)、土壤JKR表面能、秸稈-土壤JKR表面能的參數(shù)水平值,k為虛擬空白列。
Plackett-Burman篩選試驗(yàn)設(shè)計(jì)如表4和圖4所示,利用Design-Expert軟件對(duì)表4中數(shù)據(jù)進(jìn)行處理,得到各參數(shù)對(duì)堆積角影響顯著性順序如表5所示。由表5可知,顯著性由大到小排序?yàn)椋和寥?土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、秸稈-土壤滾動(dòng)摩擦因數(shù)、土壤JKR表面能、土壤-土壤恢復(fù)系數(shù)、秸稈-土壤JKR表面能、土壤-鋼滾動(dòng)摩擦因數(shù)、秸稈-土壤靜摩擦因數(shù)、秸稈-土壤恢復(fù)系數(shù)、土壤-土壤靜摩擦因數(shù)。其中對(duì)堆積角貢獻(xiàn)度較大的4個(gè)因素分別為:土壤-土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、秸稈-土壤滾動(dòng)摩擦因數(shù)、土壤JKR表面能,合計(jì)貢獻(xiàn)度93.45%,其他6個(gè)因素對(duì)堆積角貢獻(xiàn)度較小,均小于5%。因此,結(jié)合貢獻(xiàn)度和顯著性順序,選取土壤-土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、秸稈-土壤滾動(dòng)摩擦因數(shù)、土壤JKR表面能這4個(gè)因素進(jìn)行后續(xù)最陡爬坡試驗(yàn)。
表4 篩選試驗(yàn)設(shè)計(jì)及結(jié)果 Tab.4 Results of Plackett-Burman design (P-BD) test
在分析Plackett-Burman試驗(yàn)結(jié)果的基礎(chǔ)上,對(duì)篩選出的4個(gè)貢獻(xiàn)度較大的因素(土壤-土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、秸稈-土壤滾動(dòng)摩擦因數(shù)、土壤JKR表面能)進(jìn)行最陡爬坡試驗(yàn),計(jì)算實(shí)際堆積角與仿真堆積角的相對(duì)誤差,以確定仿真參數(shù)較優(yōu)范圍區(qū)間。仿真過(guò)程中,對(duì)堆積角貢獻(xiàn)度較小的參數(shù)取中間水平值,即秸稈-土壤恢復(fù)系數(shù)0.4、秸稈-土壤靜摩擦因數(shù)0.5、土壤-土壤恢復(fù)系數(shù)0.4、土壤-土壤靜摩擦因數(shù)0.6、土壤-鋼滾動(dòng)摩擦因數(shù)0.15、秸稈-土壤JKR表面能0.5 J/m2。最陡爬坡試驗(yàn)設(shè)計(jì)及結(jié)果如表6所示,結(jié)果顯示,在第4組試驗(yàn)中,堆積角相對(duì)誤差最小,可以確定最優(yōu)區(qū)間在第4組參數(shù)值附近,因此將第4組參數(shù)值作為中間水平,第3組和第5組參數(shù)值分別作為低水平和高水平進(jìn)行后續(xù)Box-Behnken試驗(yàn)。
表5 Plackett-Burman篩選試驗(yàn)結(jié)果分析Tab.5 Analysis of Plackett-Burman design test results
表6 最陡爬坡試驗(yàn)設(shè)計(jì)及結(jié)果Tab.6 Design and results of climbing test
為尋求仿真試驗(yàn)中土壤-土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、秸稈-土壤滾動(dòng)摩擦因數(shù)、土壤JKR表面能的最優(yōu)參數(shù)組合,基于爬坡試驗(yàn)結(jié)果,以堆積角為試驗(yàn)指標(biāo),根據(jù)Box-Behnken試驗(yàn)原理[29-30]進(jìn)行四因素三水平試驗(yàn)設(shè)計(jì),共進(jìn)行29次試驗(yàn),試驗(yàn)方案如表7所示。
利用Design-Expert軟件對(duì)試驗(yàn)結(jié)果進(jìn)行多元擬合分析,建立堆積角與4個(gè)自變量間的二階回歸模型,其二次回歸方程為
Y=41.48-0.21c+1.80f+3.24g-1.56i+
0.90cf-0.72cg-0.50ci-1.85fg-
2.20fi-1.19gi-2.46c2-0.22f2-
1.51g2-1.68i2
(10)
表7 試驗(yàn)方案及結(jié)果Tab.7 Test scheme and results
在保證模型顯著的情況下,對(duì)模型進(jìn)行優(yōu)化,去除對(duì)堆積角影響不顯著的項(xiàng)(c、cf、cg、ci、f2),優(yōu)化后方差分析結(jié)果如表9所示,試驗(yàn)精密度為20.095,提高了試驗(yàn)精密度,說(shuō)明模型可較好預(yù)測(cè)堆積角與4個(gè)因素的關(guān)系。優(yōu)化后二次回歸方程為
表8 回歸模型方差分析 Tab.8 ANOVA of quadratic model
Y=41.34+1.80f+3.24g-1.56i-1.85fg-
2.20fi-1.20gi-2.42c2-1.47g2-1.63i2
(11)
表9 回歸模型優(yōu)化方差分析Tab.9 ANOVA of optimization quadratic model
本試驗(yàn)以堆積角為響應(yīng)值,利用Design-Expert軟件對(duì)數(shù)據(jù)進(jìn)行多元回歸擬合,生成響應(yīng)面如圖5所示,進(jìn)一步分析各影響因素對(duì)響應(yīng)值的交互作用。由圖5a可知,土壤-土壤滾動(dòng)摩擦因數(shù)f和土壤-鋼靜摩擦因數(shù)g響應(yīng)面曲線走勢(shì)相似,表明兩個(gè)因素對(duì)堆積角影響基本相同;由圖5b可知,相對(duì)于土壤JKR表面能i,土壤-土壤滾動(dòng)摩擦因數(shù)f的響應(yīng)面曲線較陡,表明其對(duì)堆積角影響更為顯著;由圖5c可知,相對(duì)于土壤JKR表面能i,土壤-鋼靜摩擦因數(shù)g的響應(yīng)面坡度較大,表明其對(duì)堆積角影響較大。
在因素水平范圍內(nèi),利用Design-Expert軟件尋優(yōu)功能對(duì)參數(shù)進(jìn)行優(yōu)化,以目標(biāo)值40.17°對(duì)模型進(jìn)行求解,對(duì)若干組解進(jìn)行仿真驗(yàn)證,選取與物理試驗(yàn)形狀最為接近的一組解,即秸稈-土壤滾動(dòng)摩擦因數(shù)0.16、土壤-土壤滾動(dòng)摩擦因數(shù)0.24、土壤-鋼靜摩擦因數(shù)0.75、土壤JKR表面能0.67 J/m2,其他非顯著性參數(shù)取中間水平值。利用優(yōu)化后的參數(shù)進(jìn)行3次重復(fù)仿真試驗(yàn),得到仿真堆積角平均值為40.83°,與實(shí)際物理堆積角的相對(duì)誤差為1.64%,試驗(yàn)對(duì)比見(jiàn)圖6,進(jìn)一步證明了仿真試驗(yàn)的有效性。
(1)針對(duì)黑土區(qū)保護(hù)性耕作中的玉米秸稈-土壤混料,基于離散元EDEM仿真軟件,選用Hertz-Mindlin with JKR接觸模型對(duì)玉米秸稈-土壤混料進(jìn)行離散元仿真并標(biāo)定相關(guān)參數(shù)。
(2)采用物理試驗(yàn)與仿真試驗(yàn)相結(jié)合的方法,利用Design-Expert軟件,由Plackett-Burman試驗(yàn)篩選出對(duì)秸稈-土壤混料堆積角有顯著影響的因素為土壤-土壤滾動(dòng)摩擦因數(shù)、土壤-鋼靜摩擦因數(shù)、土壤-秸稈滾動(dòng)摩擦因數(shù)、土壤JKR表面能。通過(guò)Box-Behnken試驗(yàn),建立堆積角與4個(gè)因素間的二階回歸模型并進(jìn)行方差及回歸模型交互效應(yīng)分析。
(3)以玉米秸稈-土壤混料物理堆積角為目標(biāo),對(duì)影響參數(shù)進(jìn)行優(yōu)化,得到最優(yōu)參數(shù)組合為秸稈-土壤滾動(dòng)摩擦因數(shù)0.16、土壤-土壤滾動(dòng)摩擦因數(shù)0.24、土壤-鋼靜摩擦因數(shù)0.75、土壤JKR表面能0.67 J/m2。通過(guò)仿真試驗(yàn)對(duì)最優(yōu)參數(shù)組合進(jìn)行驗(yàn)證,最優(yōu)參數(shù)組合堆積角與實(shí)際物理堆積角相對(duì)誤差為1.64%,驗(yàn)證了仿真模型參數(shù)的可靠性。