黑土區玉米秸稈-土壤混料離散元模型參數標定

田辛亮 叢 旭 齊江濤 郭 慧 李 茂 范旭輝

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

0 引言

保護性耕作技術是以機械化為手段，利用秸稈或殘茬覆蓋于地表，對于提高土壤有機質含量及降低土壤風蝕水蝕具有重要作用[1-2]。玉米秸稈是一種重要的資源，對提高土壤肥力具有重要作用[3]。秸稈還田是一項有效的保護性耕作措施，可提高土壤的蓄水保墑能力，改善土壤結構，降低秸稈焚燒對環境的危害[4-5]。采用松耙聯合整地機將玉米秸稈切斷，部分玉米秸稈混入耕層，可達到加速秸稈腐爛、培肥地力的目的[6]。玉米秸稈-土壤混料特性復雜，采用常規方法難以準確獲得相關接觸參數，故利用離散元仿真模擬的方法研究土壤及玉米秸稈的相互作用規律，對秸稈還田及黑土區土壤保護具有重要意義[7]。

隨著信息技術的發展，離散元法在農業裝備研究中應用越來越廣泛[8-9]。構建離散元仿真模型，本征參數可通過試驗測出，由于顆粒間接觸力學特性復雜，部分接觸物理參數難以通過常規試驗獲取，因此近幾年國內外學者針對部分接觸物理參數難以準確獲取的問題，提出了基于離散元法對顆粒物料進行虛擬標定的方法[10-11]。目前，研究人員針對土壤[12-14]、種子[15-16]、秸稈[17]、肥料[18-19]等利用離散元模型進行了仿真參數標定。文獻[20]提出了一種通過試驗與模擬相結合系統地標定沙土顆粒相互作用參數的方法；文獻[21]利用試驗測定及仿真模擬相結合的方法對微型薯顆粒離散元參數進行了標定和校準；文獻[16]結合臺架試驗和仿真試驗，對三七種子進行仿真參數標定，并選取三七精密排種器進行了驗證試驗；文獻[17]利用離散元軟件對玉米秸稈的接觸物理參數進行仿真標定；文獻[18]對4種不同測試方法進行篩選試驗，提出了一種基于顆粒物料整體特性的摩擦因數標定方法；文獻[22]將非線性彈性和線性滯回彈簧(塑性)接觸模型的模擬結果與掃描式耕作工具的試驗結果進行比較，得到了土壤接觸物理參數的變化規律。綜上所述，國內外研究人員的離散元仿真參數研究對象主要是單獨的土壤、種子、秸稈、肥料等，但對玉米秸稈-土壤混料的研究卻鮮有報道。黑土區保護性耕作模式中玉米秸稈-土壤混料沒有明顯規律可尋，且有限元玉米秸稈-土壤混料模型不能模擬土壤及秸稈的運動過程，因此有必要對其展開研究。

本文以黑土區保護性耕作模式中玉米秸稈-土壤混料為研究對象，結合物理試驗與仿真試驗，應用EDEM軟件對玉米秸稈-土壤混料的接觸參數進行標定。以堆積角為響應值，采用圖像處理方法測量混料堆積角，依次通過Plackett-Burman篩選試驗、爬坡試驗和Box-Behnken試驗對玉米秸稈-土壤混料堆積角進行離散元仿真和校準，以確定較優的玉米秸稈-土壤混料離散元仿真參數值。

1 材料與方法

1.1 試驗材料

本文研究對象為玉米秸稈-土壤混料，來源于長春市農安縣開安鎮試驗基地(44°10′N，120°10′E)，該基地種植作物為玉米，土壤類型為黑鈣土，土壤pH值為7.01，年降水量535 mm。試驗樣品取自2020年秋季松耙聯合整地機作業后，經田間取樣檢測，玉米秸稈在土壤中的混入深度主要為0～20 cm，故本文取樣深度為0～20 cm。取樣測量利用S型法按照“隨機、等量和多點”的原則進行[23]，選取500 g玉米秸稈-土壤混料為樣品，將樣品中的土壤與秸稈進行分離，分別測量土壤與秸稈的密度；依據機具作業后秸稈實際形狀，秸稈密度按照文獻[17]的方法進行測量。采用干燥法測量土壤與秸稈的含水率；將采集土樣進行干燥，使用不同孔徑的細篩進行篩分，結果表明土壤粒徑分布范圍為2～4 mm。樣品的基本性質如表1所示。

表1 樣品基本性質Tab.1 Basic properties of sample

1.2 試驗方法

運用堆積角物理試驗與EDEM仿真試驗相結合的方法對玉米秸稈-土壤混料的離散元模型參數進行標定，采用圓筒提升法[24]進行玉米秸稈-土壤混料堆積角試驗，測量玉米秸稈-土壤混料堆積角的實際均值。基于EDEM軟件進行仿真試驗，應用Design-Expert軟件進行Plackett-Burman篩選試驗設計，并進行多因素顯著性分析篩選出對玉米秸稈-土壤混料有顯著影響的參數；利用最陡爬坡試驗確定顯著性影響參數最優區間，根據Box-Behnken試驗建立并優化玉米秸稈-土壤混料堆積角與顯著性參數的回歸模型，并對參數進行綜合優化得到最優參數組合。最后利用最優參數組合進行仿真試驗，對比玉米秸稈-土壤混料仿真堆積角與實際堆積角的差異，驗證標定模型參數的準確性。

1.3 玉米秸稈-土壤混料堆積角物理試驗

玉米秸稈-土壤混料堆積角物理試驗裝置如圖1所示，由WDW-20型萬能試驗機、鋼制圓筒、鋼板組成。鋼制圓筒高度與直徑比為3∶1[25-26]，圓筒高度為210 mm，直徑為70 mm，圓筒底部與鋼板(300 mm×300 mm)接觸，圓筒內裝填玉米秸稈-土壤混料。利用WDW-20型萬能試驗機以0.05 m/s的速度勻速提升圓筒[27]，混料經圓筒底部落于鋼板上，待混料堆穩定后利用相機垂直拍照，相機鏡頭與萬能試驗機正面平行，相機底面與萬能試驗機正面垂直，得到混料堆的正視圖像。應用Matlab軟件進行灰度處理，提取邊緣輪廓并進行線性擬合以獲得玉米秸稈-土壤混料的堆積角。重復試驗5次求平均值，得到玉米秸稈-土壤混料的堆積角為40.17°。

2 堆積角仿真試驗

2.1 接觸模型選取

玉米秸稈-土壤混料因其本身具有一定的濕度，故顆粒間存在著粘附現象。Hertz-Mindlin with JKR接觸模型是一種凝聚力接觸模型，該模型用表面能代表顆粒間的相互吸引力，適用于模擬濕粘性顆粒和團聚的物料[28]，模型的法向彈性力的實現基于Johnson-Kendall-Roberts(JKR)理論，基于法向重疊量、表面能及相互作用，計算公式為

(1)

(2)

(3)

(4)

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

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

δ——兩接觸顆粒之間法向重疊量，m

α——兩接觸顆粒間的接觸圓半徑，m

γ——表面能，N/m

R*——等效接觸半徑，m

R1、R2——兩顆粒的接觸半徑，m

ν1、ν2——兩接觸顆粒泊松比

E1、E2——兩接觸顆粒彈性模量，Pa

當γ為0時，JKR法向力與Hertz-Mindlin法向力FHertz相同，即

(5)

若顆粒不是直接接觸，該模型也可以提供吸引凝聚力，顆粒間的非零凝聚力最大間隙為

(6)

(7)

式中δc——顆粒間有非零凝聚力時的法向最大間隙，m

αc——顆粒間有非零凝聚力時的切向最大間隙，m

當顆粒并非實際接觸且間隙小于δc時，凝聚力達到最大值，即

(8)

當選用JKR模型模擬含水率較高的顆粒時，2個顆粒分開所需要的力取決于表面能和接觸角，計算公式為

(9)

式中θ——兩顆粒的接觸角

玉米秸稈-土壤混料本身含水率相對較高，混料顆粒間受水分子影響存在著粘附現象，而JKR模型中表面能可較好模擬顆粒之間的粘結力，故本文選用Hertz-Mindlin with JKR接觸模型，并依此進行參數標定。

2.2 仿真模型建立

應用離散元仿真軟件及EDEM建立仿真模型，土壤顆粒采用隨機分布的方法生成[12]，以球體作為仿真顆粒模型，使用不同孔徑的細篩進行篩分后確定土壤粒徑為2～4 mm。玉米秸稈形狀近似于長條狀，松耙聯合整地機對全量秸稈覆蓋的耕地進行作業，秸稈在作業過程中經過2次切斷，依據隨機多點的原則對作業后的秸稈長度進行測量，確定秸稈長度為40 mm[4]。目前玉米秸稈在離散元軟件中建模時被認為剛性體，通過球面堆積的方法組合形成玉米秸稈，故離散元模型設置秸稈長度為40 mm，由半徑為3 mm的12個基本球體相連接組合而成。在SolidWorks軟件中建立三維鋼制圓筒模型，把模型另存為IGS格式，導入EDEM軟件中；鋼板模型由EDEM軟件生成，材料屬性設置為鋼，仿真試驗中物料的本征參數如表2所示。

表2 物料本征參數Tab.2 Material intrinsic parameters

2.3 EDEM軟件仿真設置

仿真過程中，在鋼制圓筒上設置顆粒工廠，顆粒的生成方式設置為static，數據保存間隔為0.01 s，重力加速度為9.81 m/s2。根據田間篩分實測結果，土壤與秸稈質量比為99∶1。土壤和秸稈按照質量比采用分層交錯生成的方式，各生成7層。為保證土壤與秸稈混合的均勻性，參照文獻[25]的試驗方法，采用旋轉圓筒的方法實現土壤與秸稈的均勻混合。經多次重復試驗結果對比，顆粒生成后設置圓筒轉速為60 r/min，圓筒內物料均勻混合時間為10 s。待顆粒達到穩定狀態后，圓筒以0.05 m/s的速度勻速豎直向上運動，顆粒在重力作用下滑落在鋼板上，直至所有顆粒停止運動并形成穩定顆粒堆并測量堆積角，如圖2、3所示。

在前期反復試驗的基礎上，參照文獻[12,17,21,23,25]的參數選取范圍，確定玉米秸稈-土壤混料接觸物理參數范圍如表3所示。

表3 離散元仿真參數與水平Tab.3 Parameters and level required in DEM simulation

3 結果與分析

3.1 Plackett-Burman篩選試驗及顯著性分析

利用Design-Expert軟件進行Plackett-Burman篩選試驗設計，以玉米秸稈-土壤混料堆積角為試驗指標，篩選出對指標存在顯著性影響的參數。篩選試驗共計12組，a～j分別為秸稈-土壤恢復系數、秸稈-土壤靜摩擦因數、秸稈-土壤滾動摩擦因數、土壤-土壤恢復系數、土壤-土壤靜摩擦因數、土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、土壤-鋼滾動摩擦因數、土壤JKR表面能、秸稈-土壤JKR表面能的參數水平值，k為虛擬空白列。

Plackett-Burman篩選試驗設計如表4和圖4所示，利用Design-Expert軟件對表4中數據進行處理，得到各參數對堆積角影響顯著性順序如表5所示。由表5可知，顯著性由大到小排序為：土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、秸稈-土壤滾動摩擦因數、土壤JKR表面能、土壤-土壤恢復系數、秸稈-土壤JKR表面能、土壤-鋼滾動摩擦因數、秸稈-土壤靜摩擦因數、秸稈-土壤恢復系數、土壤-土壤靜摩擦因數。其中對堆積角貢獻度較大的4個因素分別為：土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、秸稈-土壤滾動摩擦因數、土壤JKR表面能，合計貢獻度93.45%，其他6個因素對堆積角貢獻度較小，均小于5%。因此，結合貢獻度和顯著性順序，選取土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、秸稈-土壤滾動摩擦因數、土壤JKR表面能這4個因素進行后續最陡爬坡試驗。

表4 篩選試驗設計及結果 Tab.4 Results of Plackett-Burman design (P-BD) test

3.2 最陡爬坡試驗分析

在分析Plackett-Burman試驗結果的基礎上，對篩選出的4個貢獻度較大的因素(土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、秸稈-土壤滾動摩擦因數、土壤JKR表面能)進行最陡爬坡試驗，計算實際堆積角與仿真堆積角的相對誤差，以確定仿真參數較優范圍區間。仿真過程中，對堆積角貢獻度較小的參數取中間水平值，即秸稈-土壤恢復系數0.4、秸稈-土壤靜摩擦因數0.5、土壤-土壤恢復系數0.4、土壤-土壤靜摩擦因數0.6、土壤-鋼滾動摩擦因數0.15、秸稈-土壤JKR表面能0.5 J/m2。最陡爬坡試驗設計及結果如表6所示，結果顯示，在第4組試驗中，堆積角相對誤差最小，可以確定最優區間在第4組參數值附近，因此將第4組參數值作為中間水平，第3組和第5組參數值分別作為低水平和高水平進行后續Box-Behnken試驗。

表5 Plackett-Burman篩選試驗結果分析Tab.5 Analysis of Plackett-Burman design test results

表6 最陡爬坡試驗設計及結果Tab.6 Design and results of climbing test

3.3 Box-Behnken 試驗設計

為尋求仿真試驗中土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、秸稈-土壤滾動摩擦因數、土壤JKR表面能的最優參數組合，基于爬坡試驗結果，以堆積角為試驗指標，根據Box-Behnken試驗原理[29-30]進行四因素三水平試驗設計，共進行29次試驗，試驗方案如表7所示。

利用Design-Expert軟件對試驗結果進行多元擬合分析，建立堆積角與4個自變量間的二階回歸模型，其二次回歸方程為

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 試驗方案及結果Tab.7 Test scheme and results

在保證模型顯著的情況下，對模型進行優化，去除對堆積角影響不顯著的項(c、cf、cg、ci、f2)，優化后方差分析結果如表9所示，試驗精密度為20.095，提高了試驗精密度，說明模型可較好預測堆積角與4個因素的關系。優化后二次回歸方程為

表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 回歸模型優化方差分析Tab.9 ANOVA of optimization quadratic model

3.4 回歸模型交互效應分析

本試驗以堆積角為響應值，利用Design-Expert軟件對數據進行多元回歸擬合，生成響應面如圖5所示，進一步分析各影響因素對響應值的交互作用。由圖5a可知，土壤-土壤滾動摩擦因數f和土壤-鋼靜摩擦因數g響應面曲線走勢相似，表明兩個因素對堆積角影響基本相同；由圖5b可知，相對于土壤JKR表面能i，土壤-土壤滾動摩擦因數f的響應面曲線較陡，表明其對堆積角影響更為顯著；由圖5c可知，相對于土壤JKR表面能i，土壤-鋼靜摩擦因數g的響應面坡度較大，表明其對堆積角影響較大。

3.5 參數優化及仿真驗證

在因素水平范圍內，利用Design-Expert軟件尋優功能對參數進行優化，以目標值40.17°對模型進行求解，對若干組解進行仿真驗證，選取與物理試驗形狀最為接近的一組解，即秸稈-土壤滾動摩擦因數0.16、土壤-土壤滾動摩擦因數0.24、土壤-鋼靜摩擦因數0.75、土壤JKR表面能0.67 J/m2，其他非顯著性參數取中間水平值。利用優化后的參數進行3次重復仿真試驗，得到仿真堆積角平均值為40.83°，與實際物理堆積角的相對誤差為1.64%，試驗對比見圖6，進一步證明了仿真試驗的有效性。

4 結論

(1)針對黑土區保護性耕作中的玉米秸稈-土壤混料，基于離散元EDEM仿真軟件，選用Hertz-Mindlin with JKR接觸模型對玉米秸稈-土壤混料進行離散元仿真并標定相關參數。

(2)采用物理試驗與仿真試驗相結合的方法，利用Design-Expert軟件，由Plackett-Burman試驗篩選出對秸稈-土壤混料堆積角有顯著影響的因素為土壤-土壤滾動摩擦因數、土壤-鋼靜摩擦因數、土壤-秸稈滾動摩擦因數、土壤JKR表面能。通過Box-Behnken試驗，建立堆積角與4個因素間的二階回歸模型并進行方差及回歸模型交互效應分析。

(3)以玉米秸稈-土壤混料物理堆積角為目標，對影響參數進行優化，得到最優參數組合為秸稈-土壤滾動摩擦因數0.16、土壤-土壤滾動摩擦因數0.24、土壤-鋼靜摩擦因數0.75、土壤JKR表面能0.67 J/m2。通過仿真試驗對最優參數組合進行驗證，最優參數組合堆積角與實際物理堆積角相對誤差為1.64%，驗證了仿真模型參數的可靠性。