基于EDEM 的粉壟和旋耕作業混肥效果研究

單發科， 康朔， 朱建錫， 王永維， 王俊*

（1.浙江大學生物系統工程與食品科學學院，杭州 310058； 2.浙江省農業機械研究院，浙江 金華 321051）

土壤施肥是提高作物產量的重要措施，化肥顆粒的分布范圍和作物根系的吸收范圍是否匹配及化肥分布的均勻性是影響化肥利用率的重要因素[1-2]。利用耕作機械將撒施在土壤表面的化肥顆粒與土壤混合是常用的施肥方法，不同耕作方式下的化肥分布存在差異，對作物生長產生不同影響[3-5]。石麗紅等[6]研究發現，在水稻生產中，適當的肥料分布有利于提高水稻產量；張擁兵等[7]研究發現，不同施肥方式對葡萄果實品質影響較大，而施肥方式的不同也導致肥料分布的不同；雷金銀等[8]研究發現，不同耕作方式改變了養分在土壤內的分布，提高了化肥利用率?；试谕寥乐械姆植挤秶绊懽魑锷L，而耕作方式則直接影響化肥顆粒在耕層內的分布[9-11]。因此，有必要研究不同耕作方式對化肥分布的影響。

離散元分析近年來被廣泛應用于農業機械的研究中，目前針對施肥效果的仿真研究也多采用離散元法（discrete element method，DEM）。EDEM（expanded DEM）軟件近年來被廣泛應用于機械優化設計和作業過程仿真模擬等[12-13]研究，并取得一定進展。在EDEM 中，土壤由離散的顆粒組合而成，可以通過定義顆粒的大小、尺寸、位置組合以及結合方式等模擬不同的土壤特性[14-16]，DEM 是模擬耕作過程中機具與土壤相互作用的有效方式[17-18]?；暑w粒為球狀，在EDEM 中可以很方便地建立其仿真模型。

通過田間試驗的探究不同耕作方式的混肥效果時步驟復雜且耗時耗力，為研究帶來極大困難。因此，本研究通過分析化肥顆粒分布規律建立仿真預測模型并對旋耕、粉壟2 種耕作方式的適用性進行了驗證。本研究旨在通過對仿真參數進行標定，基于EDEM 建立一種適用于不同耕作方式的離散元高精度混肥仿真模型。

1 材料與方法

1.1 刀具及整機結構

粉壟機具為浙江省農業機械研究院制造的五軸粉壟機，其主要工作部件為5 把螺旋粉壟刀具。粉壟刀具主體由單螺旋刀片和側切削刃共同組成，刀軸直徑（D）80 mm，螺旋刀片半徑（a）80 mm，總高度（Hz）400 mm，螺旋角（α）28°，螺旋線圈數（n）為3。

旋耕機為河北耕耘農業機械制造有限公司生產的170 型旋耕機，用Handyscan 手持式3D 掃描儀獲得170 型旋耕機的旋耕刀三維點云模型后進行逆向建模。

使用SolidWorks2018 建立三維模型，對實際機具進行1∶1 原尺寸還原，并將殼體、箱體、懸掛等非入土結構簡化處理，保留耕作刀具。結合農藝要求，2種機具的耕作參數如表1所示。

表1 耕作參數Table 1 Tillage parameters

1.2 主要材料參數標定

土壤樣本取自浙江省金華市婺城區徐村農業試驗基地的試驗田，取樣方法參照《農業機械試驗條件測定方法》[19]。使用與環刀同尺寸的光滑薄壁塑料管獲取自然狀態下不同土層深度（0—10、10—20、20—30、30—40 cm）的原樣土壤樣本?；蕵颖救∽詫嶋H耕作時使用的顆粒狀復合肥（硫酸鉀型復合化肥，火山能公司）。

土壤屬性參數包括彈性模量、剪切模量、泊松比、含水率、密度等。使用TMS-PRO 質構儀（美國FTC 公司，檢測精度>0.015%）于實驗室標定土壤彈性模量，得到力-位移曲線。質構儀的壓縮速度為15 mm·min-1，并設置在正壓力到達500 N 時停止壓縮。每層土壤深度下做4 組重復，共做16 組壓縮試驗，試驗中一組土壤內部摻雜石塊，所得曲線出現明顯偏差，將其剔除后共獲得15 條曲線。參照相關研究進展[13-15,23-25]選取土壤泊松比，通過計算獲得土壤剪切模量。

式中，E為土壤彈性模量；G為土壤剪切模量；F為正壓力，N；S為土樣截面積，m3；dL為壓縮距離，m；L為土樣原長，m；v為土壤泊松比。

采用環刀法測定土壤容重，所用環刀為100 cm3定體積環刀（南京土壤儀器廠）。土壤含水率通過烘箱干燥法測定，將環刀取得的土樣放入烘箱（GZX-9140ME 數顯吹風干燥機，上海博迅實業有限公司）105 ℃烘干至恒重，土樣做2 份平行測定?；谝韵鹿接嬎阃寥廊葜睾屯寥篮蔥21]。

式中，ρs為土壤容重，g·cm-3；m1為環刀的重量，g；m2為土壤和新鮮土壤的重量，g；V為環刀體積，cm3；Ws為樣品的土壤含水量g?kg-1；Ws為土壤含水率；Wf和Wd為土壤的鮮重和干重，g。

使用排油法測定化肥顆粒真實密度。供試材料的其余屬性參數參考相關文獻[22-24]。

不同材料間接觸參數包括土-土間接觸參數、土-肥間接觸參數、土-鋼間接觸參數和肥-鋼間接觸參數；材料間接觸參數包括碰撞恢復系數、滑動摩擦系數和滾動摩擦系數。采用下落法標定碰撞恢復系數，以土-肥間接觸參數標定為例，將實地取得的原樣土壤平鋪在平板上，并進行壓實以模擬真實土況；利用高速相機，記錄從固定高度自由下落的化肥顆粒在與土壤平面發生碰撞后反彈的最高高度，并通過式（5）進行計算。

式中，ε為碰撞恢復系數；?H為反彈高度，m；H為下落高度，mm。

每種材料間碰撞恢復系數做20 組重復，取平均值。采用自行搭建的摩擦系數測定實驗平臺（圖1）測定摩擦系數，通過改變表面覆蓋材料測定不同材料間的摩擦系數[25]。測定滑動摩擦系數時，將一種待測材料覆蓋固定于測試平臺之上，另一種待測材料利用環刀制成大小均勻的圓餅狀置于其上，以保證材料間不發生滾動；測定滾動摩擦系數時，將一種待測材料覆蓋固定于測試平臺之上，另一種待測材料制成大小均勻的圓球狀置于其上，以保證材料間不發生滑動。測定過程中記錄所有材料質心的位置，并通過高速相機驗證材料間的運動情況，通過數顯式旋鈕調速器調節轉速，測量物體被拋離的距離和拋離高度，計算摩擦系數。每組實驗通過多次重復試驗以減小誤差。

圖1 摩擦系數測定平臺及測定原理Fig. 1 Coefficient measuring platform and measuring principle

式中，μ為摩擦系數；r為旋轉半徑，m；h為墜落高度，m；l為拋離距離，m。

1.3 EDEM 仿真模型建立

使用離散元分析軟件EDEM2018建立仿真模型，包括刀具模型和導入以及土槽模型。

依據刀具尺寸和耕深，設置土槽尺寸為400 mm×500 mm×1 200 mm。試驗土槽每層土壤厚度為10 cm，每層生成15 萬個土壤顆粒，共計60 萬個[26]。結合試驗田土壤特性選取合適的模型： 0—10 cm 土壤松散，呈粉壤土性質，選取Hertz-Mindlin 無滑動模型；10—40 cm 土層，土壤間存在明顯黏結現象，呈黏性土性質，選取Hertz-Mindlin with Bonding 黏結模型對土壤顆粒進行約束，根據BPM 黏結模型，法向剛度和切向剛度由式（7）和式（8）計算[27-29]，黏結半徑計算公式（9）可由土壤含水率的計算公式推導[18]。

式中，kn為法向剛度，N·m-1；ks為切向剛度，N·m-1；a 為常數，0—10 cm 土層取2/3，10—40 cm土層取1/2；ν1為顆粒1 的泊松比；ν2為顆粒2 的泊松比；E1為顆粒1 的彈性模量，Pa；E2為顆粒2 的彈性模量，Pa；r1為顆粒1 的半徑，mm；r2為顆粒2的半徑，mm；ω為土壤含水率；m1為土壤顆粒質量，kg；m2為水分質量，kg；ρ1為壤顆粒干密度，kg·m-3；ρ2為水密度，kg·m-3；R為土壤顆粒半徑，m；RB為土壤顆粒黏結半徑，m。

結合實際施肥，通過設置面顆粒工廠將化肥顆粒隨機撒施到土壤表面。為提高仿真精度，同時控制計算時長，結合總計生成的顆粒數量及相關參考文獻對顆粒半徑進行放尺，放尺后土壤顆粒半徑為4 mm，顆粒半徑對仿真精度影響不大[30-31]。土、肥顆粒模型生成后，讓其自然沉降3 s，生成黏結鍵[30]。

1.4 仿真模型準確性驗證

通過耕作過程牽引力值驗證仿真模型參數標定結果的準確性。設置仿真和田間試驗時機具工作參數相同，粉壟耕作前進速度為0.3 m·s-1，刀具轉速為300 r·min-1，刀具傾角0°；旋耕耕作前進速度0.2 m·s-1，刀具轉速400 r·min-1。由EDEM軟件輸出刀具前進方向合力，即為仿真刀具水平方向牽引力。粉壟機械實際安裝過程中，根據不同旋向分別按照120°和180°間隔安裝刀具，因此仿真過程按5 把刀具排布角度進行，并將所得前進方向合力求和，得機具整體所受牽引力。

田間試驗在浙江省金華市婺城區徐村農業試驗基地進行，耕作過程中的牽引力信號通過板環式拉壓力傳感器（南京天光電氣科技有限公司，量程2 t，電壓輸出范圍0～20 mV）測得，通過數據采集系統（DH5902N 堅固型動態信號測試分析系統，江蘇東華測試技術股份有限公司）采集電信號并按比例轉化為待測力，所測得耕作力被用作驗證仿真模型計算結果的準確性。試驗所采用的板環式傳感器可以測量沿傳感器兩銷孔連線方向的力，當機具入土后保持勻速水平前進時，認為傳感器測量處于平衡態，此時可根據機具與傳感器形態位置計算牽引力。

1.5 混肥效果優化試驗

結合仿真實驗對2 種耕作方式下的混肥效果進行分析，以混合的均勻性為混肥效果評價指標，均勻性通過計算各土層化肥顆粒數變異系數（C）評價，C越小混肥更均勻。

式中，xi為第i層顆粒數，p為每層顆粒數平均值，n為層數。

結合實際農藝經驗和耕作參數（表1），選取各因素取值范圍。通過所建立的仿真模型進行試驗，獲取2 種耕作方式的分布規律預測模型，通過該預測模型計算獲得最優工作參數。試驗過程如圖2 所示，在仿真模型中設置粉壟刀具運動為旋轉下降和旋轉前進2 部分，旋耕刀具運動為旋轉前進，刀具工作參數根據正交試驗表設置，因素的編碼值和變異系數如表2 所示，其中，設置粉壟刀具傾斜角度向前進方向傾斜為負，反之則為正。

1.6 混肥效果田間驗證

田間驗證試驗在浙江省金華市婺城區徐村農業試驗基地進行。試驗田長期進行農業活動，前茬作物為水稻，水稻收割完成后使用鏵式犁對試驗田進行翻耕，并施用有機肥對試驗田理化性質進行調整改造，晾曬1 月后進行田間試驗。試驗用拖拉機型號為久保田M954KQ。由于實驗所施用的化肥易溶于水，難以統計各層的數量分布，因此選用形狀、大小、密度與化肥顆粒相似的聚丙烯顆粒以1∶1 的施用比例代替化肥。耕作前采用與撒肥相同的方法將聚丙烯顆粒均勻撒施到試驗區域土壤表面，分別使用旋耕機和粉壟機對田地進行翻耕，將顆?；烊胪寥?。取樣時，使用10 cm×10 cm×10 cm 合金鋼方框以10 cm 為梯度在0—40 cm 土層的耕作擾動區域內取樣，不同層的土壤分別裝入不同的取樣袋中。分別通過5 和2 mm 孔徑篩子篩選，手動挑揀出聚丙烯顆粒，計數并統計顆粒分布情況。

2 結果與分析

2.1 參數標定結果

結合圖3 的結果發現，本研究所取土壤樣本在壓縮過程中的力學特性表現類似于塑性材料，存在彈性階段（OA）和屈服階段，且進入屈服階段以后，試樣越壓越扁，抗壓能力繼續增強，因此土壤彈性模量類比線彈性材料，通過A 處獲得的軸向載荷計算得到。

圖3 土壤力學分析Fig. 3 Soil mechanics analysis

最終得到的材料屬性參數和接觸參數如表3和表4 所示?？梢钥闯?，0—10、10—30 和30—40 cm 土層土壤物性存在明顯差異，隨著土層深度的增加，土壤被壓實，土壤含水率也增加。最終確定黏結模型參數如表5 所示?；暑w粒的半徑通過實際取樣，測量得平均值為1.94 mm，因此在EDEM 中設定模型半徑為2 mm，根據實際施肥量計算得到試驗區域面積內施用的顆粒數目為2 000個。

表3 材料屬性參數Table 3 Material property parameter

表4 材料接觸參數Table 4 Material contact parameter

表5 黏結模型參數Table 5 Bonding model parameter

2.2 模型準確性分析

粉壟耕作驗證結果如圖4A 所示，將實際情況下機具受力信號分別放縮10%，作為仿真值的10%誤差區間。粉壟耕作過程中，仿真所得刀具受力在入土初期波動較大，且呈現一定的周期性，波動范圍為-100～300 N，總體呈上升趨勢，這可能是由于機具入土初期，刀具入土深度較淺，與土壤相互作用部分較小，因此整體受力較小，刀具周期性與土壤的相互作用對受力的影響大于接觸面積的影響。在入土中期，受力大幅提升，曲線波動不明顯，這可能是由于隨著耕作深度的增加，刀具與土壤接觸面積對受力的影響大于周期性接觸的影響。在入土后期，牽引力變化趨勢和變化大小相似，但由于入土階段總體受力較小，實際值的10%誤差區間較小，仿真值僅有部分在誤差區間內。在機具平穩前進階段，仿真值基本全部位于10%的誤差區間內，牽引力大小在1 100 N 波動，這是因為機具平穩前進時，相同深度土壤性質較為相似，刀具受土壤性質變化影響小，且仿真模型土壤性質參數與實際土壤性質參數的差異不大，因此預測值與真實值擬合程度較好。

圖4 耕作牽引力曲線Fig. 4 Tillage traction force curve

旋耕機的牽引力驗證結果如圖4B所示，在同一時間長度內，仿真結果與實際信號同樣呈現較大的波動，波動原因是工作過程中旋耕刀入土數量和入土深度呈周期性波動且旋耕機耕作過程中自身抖動較大。對信號數據進行平滑處理后，發現仿真結果和試驗結果呈現相似的變化規律，波動周期大小相似但存在一定差異，差異存在原因可能是實際耕作時刀具轉速難以精確控制。由平滑曲線發現，旋耕機實際受力略大于仿真受力，這可能是由于實際工作中旋耕機存在摩擦力產生的機械損耗[19]。

綜合上述2 種模型的驗證結果，認為該仿真模型具有一定的準確性，可以用于耕作過程的模擬分析。

2.3 混肥效果優化分析

模型方差結果如表6所示，響應曲面如圖5所示。粉壟耕作模型顯著性檢驗P=0.000 3，決定系數R2=0.91，失擬項P=0.092 1，回歸模型極顯著，失擬項不顯著，擬合程度較高。由表6 可知，對均勻性的影響，A2、B2、C2極顯著，顯著性順序為A2>C2>B2。試驗因素與均勻性之間存在二次非線性關系，得到粉壟耕作（CF）的回歸模型如下。

圖5 標準差響應曲面Fig. 5 Response surface on standard deviation

表6 模型方差分析Table 6 Model analysis of variance

旋耕耕作模型顯著性檢驗P=0.000 5，決定系數R2=0.93，失擬項P=0.063 9，回歸模型顯著，失擬項不顯著，擬合程度較高。由表6 數據可知，對均勻性的影響，B、A2、B2極顯著，AB顯著，顯著性順序為A2>B2>B>AB。試驗因素與均勻性之間存在二次非線性關系，得到旋耕耕作（CX）的回歸模型如下。

結合所得回歸模型，以獲得最小變異系數為目標，優化方程求解得到最優混肥效果下的工作參數為，粉壟轉速308 r·min-1，前進速度0.28 m·s-1，傾斜角度1.7°；旋耕轉速420 r·min-1，前進速度0.25 m·s-1。根據該組參數，設計田間驗證試驗。

2.4 混肥效果田間驗證結果分析

由EDEM 后處理模塊分析化肥顆粒分布情況。分層建立土壤化肥顆粒取樣的取樣溝模型，統計每層統計模型中化肥顆粒的數量，最后按數量計算化肥顆粒在每層土層中的數量占比。仿真模型如圖6所示，最終統計結果如表7所示。

圖6 仿真結果Fig. 6 Simulation results

表7 化肥顆粒分布仿真和驗證結果Tab. 7 Simulation and verification results of fertilizer particle distribution

從表7 可以看出，在預測準確性方面，旋耕后的化肥分布預測在0—10 和20—30 cm 土層較精準，預測誤差分別為4.47%和4.18%，而10—20 cm土層的預測誤差為8.65%。旋耕耕作方式下的化肥分布總體呈由上而下逐漸減少趨勢，20—30 cm土層化肥顆粒數量占比較少。對其結果進行分析，原因可能是仿真過程中旋耕耕作不能到達20 cm 以下土層深度，但存在少量化肥顆粒沿耕作產生的空隙滑落至耕作區底部；實際取樣過程中，耕作機具耕深不穩定，時常大于設定耕深，因此實際取樣的顆粒數占比比仿真更高，從而導致顆粒數量波動較大。粉壟耕作過程的預測結果在10—20和20—30 cm土層較為精準，分別為2.32%和0.59%，但0—10 和30—40 cm 土層誤差較大，分別為8.03%和6.31%?；暑w粒分布在數量上大體呈現由上到下逐層減少的規律。對出現的較大誤差進行分析，0—10 cm 土層范圍內實際測得數量小于仿真數量，原因可能是該土層內化肥顆粒數量較多，田間取樣時，少量顆?；渲?0—20 cm 土層，而10—20 cm 土層實際取樣所得顆粒數大于仿真顆粒數也印證了這一點。30—40 cm土層誤差較大原因可能為，一方面，深層土層內，化肥顆粒數量較少，取樣時數量很少的幾?；暑w粒數量波動就會導致較大誤差的產生；另一方面，同樣存在耕深不穩定現象導致實際取樣中30—40 cm土層化肥顆粒數量占比更高。

綜合表7 數據，混肥效果仿真預測值和測量值最大誤差為8.65%，最小誤差為0.59%，且仿真與實際測得的化肥顆粒分布呈現相似的規律，理論分析與實際結果非常接近，說明該預測模型具有可靠性。

3 討論

相較于粉壟耕作方式，旋耕在耕層內能夠近似實現化肥的均勻分布，適用于根系淺、根數量多的須根系作物如水稻、小麥、蔥等。粉壟耕作的影響區域較大，在每個鉆削部件的工作直徑范圍內都產生影響，其土壤運動范圍在深度上較大，因此便于實現地表撒肥深施。粉壟耕作下耕層內的化肥數量近似倒三角分布，相較于其他耕作方式，能夠實現不同土壤深度下分布不同化肥顆粒量，更加貼合部分根系較深作物的生長過程，如茄子、青花菜、蘿卜等。因此，結合不同作物的生長習性和種植農藝，利用預測模型對特定耕作條件下的化肥分布進行模擬預測，根據預測結果選取最優的耕作方式，可以提高農業施肥效率，改善施肥效果。

目前已有研究利用離散元方法對耕作機具進行仿真優化，陳新予等[32]應用EDEM 軟件對旋耕機拋土效果進行了驗證；鮑洋清[33]研制了新型仿生深松鏟，基于EDEM 軟件探究了深松鏟的減阻作用和土壤顆粒的運動規律；高亮[34]建立了耕刀切割土壤的離散元模型，并對耕刀形狀進行優化，進一步探究了耕刀形狀與功耗間的關系?，F有研究多從減小耕作阻力和降低能耗角度出發，模擬耕作機具與土壤間相互作用，優化工作部件設計參數，從而實現對機具結構的改進，鮮有從化肥混合效果角度出發，對機具結構進行優化的研究。

本文使用離散元法模擬了粉壟機具和旋耕機具對地表撒施化肥顆粒的混合過程和耕作過程，通過優化計算獲得最優混肥效果下的工作參數。對比田間試驗結果和仿真試驗結果，發現每層化肥顆粒占比的仿真預測值和測量值最大誤差為8.65%，最小誤差為0.59%，證明了預測模型的準確性。研究結果表明，不同耕作機具作用下的化肥顆?？臻g分布大不相同。

綜上所述，一方面，本研究結果在旋耕機具設計時可以以化肥顆粒在耕層內混合更均勻為目標對刀具結構等進行優化，在粉壟機具設計時可以以自下而上，由少到多、等差分布為目標進行優化，或結合作物根系分布范圍，以化肥分布與根系分布相匹配為目標，對刀具長度、螺旋葉片寬度等進行優化；另一方面，在機具應用時，可以通過重新標定、設置不同的仿真參數，實現不同土壤條件下的工況模擬，通過虛擬試驗方法，以達到所需的化肥混合效果為目的，對機具工作時的工作參數進行優化，提高作業質量，從改善化肥顆?？臻g分布角度出發，對旋耕、粉壟等耕作機械設計和應用提出新的要求。