基于離散元法的鋤板入土作業參數優化分析*

鄧成志，王善文，黎展鵬，陳翀，魏志強，段廷億

(西南大學工程技術學院，重慶市，400715)

0 引言

我國農業機械化進程不斷深化，在很大程度上提高了農業生產有效性，改善了農民生活水平。微耕機具有結構簡單、操作簡單及效益高等優勢，無論是大田、小田，還是平川、山地，微耕機都能夠自如操作[1]。微耕機在丘陵地區作為動力機械備受認可[2]。

莖稈類農作物收獲后根系殘留在田間不僅會破壞土壤養分平衡，降低土壤的有效養分，致使土壤中真菌數量明顯上升，同時還改變了土壤微生物的多樣性，降低了土壤酶的活性，導致作物生長不良，產量和品質下降，病蟲害發生嚴重等現象[3-5]。目前對于莖稈類農作物根系清除作業設備的研究更多的局限于配套機構設計或檢測上，李華[6]設計了一種將煙兜與煙稈連同拔起的拔稈機，只給出了該機器的設計方案和煙稈的檢測方法；陳學深等[7]設計了一種根土分離裝置，并未對該裝置與土壤相互作用進行有關分析。

近年來，離散元在農業耕作領域的應用越來越廣泛[8]。Shmulevich等[9]通過仿真和試驗，利用離散元法分析得出由于土壤在刀具前進時的堆積作用，使得刀具在運動過程中受到的水平力逐漸增大，確定了利用離散元法分析土壤耕作過程的合理性。Tamás等[10]基于離散元法建立了深松鏟前傾角和耕作速度對水土流失影響的預測模型，通過土槽試驗證明了預測模型與實際情況的相對誤差最小可以達到4%，解決了其他幾種分析土壤—刀具相互作用的方法費時費力的問題。方會敏等[11]利用EDEM軟件構建土壤—旋耕刀相互作用和秸稈—土壤—旋耕刀相互作用的三維離散元模型，得到了旋耕刀在有秸稈覆蓋的情況下扭矩和力明顯高于無秸稈覆蓋的結果。馬躍進等[12]設計了一種凸圓刃式深松減阻鏟，通過離散元仿真分析和田間試驗分析，試驗結果表明該鏟比常規深松鏟耕作阻力平均降低了10.24%。

綜上來看，針對丘陵山區的除根設備以及對根土共存土壤力學屬性的研究較少。本文應用離散元軟件EDEM建立了適用于存在根土錨固現象的土壤模型。采用JKR模型來表征土壤的應力—應變關系并設置仿真所需參數。通過Box-Behnken響應面分析法設計仿真試驗并建立了三元二次回歸數學模型，并以此模型確定了鋤板最優入土作業參數組合并試驗驗證。

1 根系掘除設備結構及工作原理

1.1 根系掘除設備結構

根據丘陵山區田間地理情況和煙草種植規模，結合杠桿原理和四桿機構特點，本文設計了一種新的適應于丘陵山區的根系掘除設備，主要包括3部分：行走部分、液壓泵站、根系掘除機構，如圖1(a)所示。行走部分以微耕機(1WG6.3-110FC-Z)作為行走動力，液壓泵(306右旋油泵動力機：額定壓力16 MPa，工作轉速800～2 000 r/min，油箱容積3.5L)提供入土動力，根系掘除機構為作業機構。除此之外，其他還包括HSG40液壓缸(行程500 mm，壓力16 MPa)、DSG-02-3C2-LW-DC24電磁換向閥(中位保壓)、MRV-02P疊加式溢流閥(作為調壓閥，改變液壓缸活塞桿的速度)、油路塊、4根高壓油管、鉛酸蓄電池(單塊電壓12 V，數量3塊)、控制系統(包括應變片信號動態檢測系統)等。鋤板結構示意圖如圖1(b)所示。圖2為煙草根系的外形和尺寸示意圖。

(a) 根系掘除設備結構圖

圖2 煙草根系的形狀和尺寸

1.2 工作原理

微耕機作為行走動力，同時其傳動箱輸出軸連接皮帶盤，通過帶傳動將動力傳給液壓泵。微耕機繼續前進，行駛至合適位置，按下液壓電磁閥控制按鈕，液壓缸伸長使得鋤板入土，鋤板轉到指定角度后，液壓缸停止伸長并鎖死，此時完成鋤板入土作業。再次按下液壓電磁閥控制按鈕，液壓缸收攏，回到初始位置。鋤板入土運動圖如圖3所示，其中v表示整機前進速度，ω表示鋤板繞轉動點旋轉的平均角速度，h表示鋤板刀尖距離地面的高度。

圖3 鋤板運動圖

2 離散元仿真模型

鋤板的入土過程中始終伴隨著大量的顆粒運動，在顆粒運動過程中始終存在著重力、摩擦力等多種作用力的綜合作用，采用EDEM軟件仿真分析鋤板入土過程中的受力情況，能真實了解影響鋤板入土作業的主要影響因素及參數的合理范圍。

2.1 根土共存土壤模型

土壤單顆粒模型、微團粒土壤模型如圖4所示。

(a) 單顆粒模型

本文選用球形顆粒模擬田間土壤顆粒，離散元仿真中顆粒尺寸增大會導致仿真耗時呈幾何級數增加，故仿真中的顆粒尺寸受計算時間和存儲空間限制[13-14]，本文選定土壤顆粒半徑[14-15]為4 mm，處于可接受范圍內。根系錨固現象會改變根系周圍土壤的抗剪強度[16-17]，本文將根系及其周圍的土壤視為大的團粒結構(團粒結構可提高土壤的抗剪強度，土粒越大，形狀越不規則，表面越是粗糙，則抗剪強度越大[18])作為根土復合體[19-20]的土壤模型。在EDEM中采用3、6、9個圓球進行多球填充建立微團粒土壤模型。

2.2 鋤板入土作業仿真模型

重慶市北碚區縉云山丘陵山區農田土壤為黃壤土[21]，具有一定粘性，在EDEM軟件中Hertz-Mindlin with JKR粘結模型用來模擬土壤顆粒間的相互作用[22]。鋤板(普通碳素鋼)和土槽(1 000 mm×500 mm×800 mm)在Creo3.0中繪制后導入EDEM；設置土壤顆粒參數，定義不同作業參數的鋤板；設置生成 80 000 個土壤顆粒來模擬土壤環境。仿真模型及土槽尺寸如圖5所示。離散元仿真參數主要分為材料參數和接觸參數兩大類，其中材料參數包括土壤、鋤板的泊松比、剪切模量、密度等；接觸參數包括土壤—土壤、土壤—鋤板的恢復系數和靜、滾動摩擦因素等。模型仿真參數[23-28]如表1所示。

圖5 3D仿真模型

表1 模型仿真參數

2.3 土壤顆粒動力學接觸模型

將土壤視為充分多的球顆粒群，每一個顆粒作為一個單元體，賦予物理性質并給予一定的質量和速度。鋤板入土過程中，土壤顆粒間會產生相互作用，可以看作是多個土壤顆粒之間發生相互碰撞，機械系統中以各種形式出現的物理現象均可以使用質量、彈性和阻尼3個要素來描述[28]。圖6(a)表示土壤顆粒碰撞的振動模型，其中s1表示顆粒1的等效位移。鋤板入土過程中土壤顆粒間也會發生土壤擾動現象，圖6(b)代表轉動模型，θ1表示顆粒1的旋轉角度。

由圖6(a)得到x、y、z方向動力學方程如式(1)所示，根據圖6(b)可得轉動力學方程如式(2)所示。

(a) 振動模型

(1)

(2)

式中：Fx、Fy、Fz——顆粒1在x、y、z方向分力；

m——土壤顆粒的等效質量；

cx、cy、cz——x、y、z方向阻尼系數；

Kx、Ky、Kz——x、y、z方向彈性系數；

sx、sy、sz——顆粒1在x、y、z方向相對位移；

Mx、My、Mz——顆粒1在x、y、z方向力矩；

r——旋轉半徑。

2.4 土壤—鋤板動力學接觸模型

建立土壤顆粒與鋤板之間的接觸模型，可以將其視為由質量、彈簧及阻尼器組成的一個簡單自由度系統，如圖7所示。

(a) 振動模型

3 仿真結果分析

3.1 仿真因素與評價指標

以前進速度、轉動角速度、離地高度作為控制因素，以鋤板所受的合力作為考察指標。仿真試驗重點考察前進速度、轉動角速度、離地高度對入土作業性能影響規律。鋤板入土前，設定土壤顆粒數量為80 000個，每種顆粒生成速率為5 000個/s。

仿真試驗時，為了準確體現鋤板所受合力變化，待土壤總動能為0時，鋤板開始運動，每隔0.02 s進行一次測量，每組試驗重復3次，最終取3組試驗的時間點對應數據平均值作為鋤板入土過程所受最大合力值。

3.2 試驗設計

采用三因素三水平的Box-Behnken Design方法對計算工況進行設計，方法參見文獻[29]，并結合農藝知識[7]，確定了各變量的水平范圍，如表2所示。采用響應面法分析3個因子對響應值(鋤板所受的最大合力值)的影響。

表2 試驗因素與水平表

3.3 試驗結果分析

3.3.1 數學模型的建立與檢驗

采用Box-Behnken中心組合進行三因素三水平的試驗設計，試驗設計共17個試驗點，序號1～12為分析因子，13～17為零點估計誤差。仿真試驗設計及結果見表3。

表3 響應面試驗方案及結果

方差分析結果見表4，在此模型中，各因素及其的交互作用、平方項對數學方程影響顯著。經過顯著性檢驗，該模型的F值為318.11，P值小于0.000 1，說明模型極其顯著。該數學模型決定系數R2為0.998 1，表明該模型擬合程度良好。通過仿真試驗得到基于實際因素的最終回歸方程如式(3)所示。根據模型各因素回歸系數絕對值大小，可以得到各因素的影響順序為前進速度>轉動角速度>離地高度。

表4 回歸分析結果

實際因素方程

F=2 507.975-17 504.75v+111.325ω-11.378 5h-

452vω+33.35vh-0.371ωh+38 757.5v2+

1.498 25ω2+0.054 63h2

(3)

3.3.2 兩因素交互作用對最大合力的影響

根據模型回歸分析結果可知，前進速度與轉動角速度交互項、前進速度與離地高度交互項、轉動角速度與離地高度交互項對鋤板所受最大合力均影響顯著，P值均小于0.05。應用軟件Design-Expert 10繪制前進速度與轉動角速度交互作用的響應曲面，可以直觀地看到這兩個參數之間的交互效應。

如圖8所示，離地高度取值為100 mm，在轉動角速度處于高水平時，前進速度對最大合力的影響顯著，表現為響應曲面速度曲線比較陡，說明了在轉動角速度較大時，增大前進速度能極為明顯改變鋤板的受力情況；而在轉動角速度比較小時，增大前進速度也能夠明顯改變鋤板的最大合力，同時也證明了二者的交互作用達到顯著。

(a) 離地高度為50 mm

從圖9可以看出，當轉動角速度為50°/s時，隨著兩參數取值的增加，所受最大合力值均呈現減小趨勢，與此同時也能明顯看出前進速度對最大合力的影響程度明顯大于離地高度。從圖10可以看出，當前進速度為0.4 m/s時，隨著兩參數的增加，所受的最大合力變化不明顯；當離地高度分別處于較高和較低水平時，轉動角速度對所受最大合力的影響程度有明顯的差異。

(a) 轉動角速度為40°/s

(a) 前進速度為0.3 m/s

3.3.3 參數最優組合的確定及仿真驗證

在軟件中設置每個參數的取值范圍，測試裝置前進速度的取值范圍是0.3～0.5 m/s，轉動角速度的取值范圍是40～60°/s，離地高度的取值范圍是50～150 mm，最大合力的最小值的取值范圍是560～2 832 N。對最大合力的回歸方程進行綜合分析并求解，試驗結果如表5所示。

表5 試驗優化結果

綜合考慮v、ω、h對F的影響并尋優，應用Design-Expert軟件響應曲面分析法求解最大合力值最小時各因素最佳參數：v=0.41、ω=40、h=115，即前進速度為0.41 m/s，轉動角速度為40°/s，離地高度為115 mm，在此條件下，鋤板所受最大合力為563.36 N。

4 田間試驗及結果

4.1 試驗目的及條件

為了驗證團粒結構的根土復合體力學屬性及離散元仿真結果的準確性，同時實際檢驗鋤板的工作性能，如圖11所示。試驗地為重慶市北碚區縉云山西南大學實踐基地煙草種植地，試驗地區土壤質地為黃壤土，土壤緊實度為0.6 MPa、土壤容重為0.99 g/cm3和土壤含水率為32%。測量設備包括：皮尺(量程30 m)、120-3AA引線應變片(量程120 Ω，敏感柵尺寸5 mm×3 mm)、E6B2-CWZ6C光電編碼器、UNI-T萬用表、PC機。利用皮尺和秒表測量整機平均前進速度，采用光電編碼器測量鋤板的平均轉速，采用電阻應變片測力方法[14]測量鋤板入土作業所受合力。

圖11 田間試驗

4.2 應變片測力標定

鋤板與轉動桿焊接，沿軸向距焊接處40 mm的轉動桿表面用砂紙打磨光亮，并利用無水乙醇擦洗打磨處，干燥后涂上一層膠水(膠水層很薄)，應變片貼于膠水處并保證粘貼處無氣泡。結合惠斯通電橋和差分放大電路，利用秤砣對重力G(重力加速度10 m/s2)與壓差ΔU的關系進行標定[11]。檢查無誤后，固定轉動桿，使鋤板中心水平，在鋤板中心處添加砝碼，做6組試驗，每組試驗重復3次，取平均值并記錄。試驗結果如表6所示，由表6得到鋤板受力F1與電壓差值ΔU之間的關系。其標定結果為

表6 鋤板標定結果

F1=158.65ΔU-2.548 1

(4)

式(4)的擬合系數R2為0.997。

4.3 田間試驗方法及驗證結果

選取煙葉收獲后的長20 m、寬15 m的試驗地塊。將液壓泵站及根系掘除機構安裝在微耕機上。通過調節微耕機發動機油門大小將微耕機前進速度控制在0.41 m/s，通過調節溢流閥將鋤板轉動角速度控制在40°/s，手動調節根系掘除機構安裝位置將離地高度控制在115 mm。應變片信號動態檢測系統采樣頻率100 Hz，田間試驗結果如圖12所示。

圖12 田間試驗結果

田間試驗實測鋤板在前進速度0.41 m/s，轉動角速度40°/s，離地高度為115 mm時所受最大合力的平均值為592.08 N，田間試驗結果同回歸模型結果相比，相對誤差為5.1%。

分析認為誤差存在的原因：(1)田間工況復雜，如存在秸稈、植物根系、碎石等，仿真中未能考慮到這些因素的存在；(2)同仿真相比，田間地表平整度不佳以及機器作業的振動，易造成數據波動；(3)整機前進速度、鋤板轉動角速度測量方法并不太精確。

總體而言，田間試驗結果與仿真模型結果基本一致，表明仿真模型模擬的根土復合體能基本模擬該試驗地的土壤力學特性。

5 結論

1) 本文創新性設計了一種適應于丘陵山區的根系掘除設備和耕作部件(鋤板)，并基于離散元法利用EDEM軟件分析根系掘除設備鋤板入土作業受力情況。

2) 基于重慶市北碚區縉云山丘陵地區土壤特性，應用離散元軟件EDEM建立了適用于存在根土錨固現象的土壤模型。采用JKR模型來表征土壤的應力—應變關系并設置仿真所需參數。通過Box-Behnken響應面分析法設計仿真試驗并建立了三元二次回歸數學模型，并以此模型確定了鋤板最優入土作業參數組合并試驗驗證。

3) 依托重慶市北碚區縉云山西南大學煙草種植地，在前進速度為0.41 m/s，轉動角速度為40°/s，離地高度為115 mm時，仿真最優值為563.36 N，試驗平均值為592.08 N，通過對比分析仿真數據數值與田間試驗鋤板所受最大合力值，兩者之間的相對誤差為5.1%。表明模擬的根土復合體的土壤力學特征基本符合重慶市北碚區縉云山西南大學煙草種植地的土壤特性。