基于離散元的稻板田旋耕功耗預測模型研究

祝英豪 夏俊芳 曾 榮 鄭 侃 杜 俊 劉政源

(1.華中農業大學工學院， 武漢 430070； 2.農業農村部長江中下游農業裝備重點實驗室， 武漢 430070)

0 引言

長江中下游多熟制稻作區主要實施稻油、稻麥水旱輪作種植模式。稻板田栽播冬油菜、冬小麥前的耕整不同于北方的以松代犁，水旱輪作制的土壤耕整不僅需要符合旱作物栽培根系生長環境要求，同時也應避免犁底層過度破壞導致蓄水困難，影響下茬水稻種植。由于北方模式不適用，為解決南方稻板田土壤粘結嚴重、茬口農時矛盾、前茬稻稈量大株高等耕作難點，以旋代犁的土壤旋耕秸稈旋埋技術及裝備為其提供了新的解決方案[1-2]。其核心部件旋埋刀輥通過在傳統旋耕刀輥內加裝切幅較寬的螺旋橫刀，起到防纏繞、秸稈旋埋以及提升碎土率的作用，但由此導致的高功耗問題尤為突出。因此研究旋埋刀輥的減阻降耗對于稻板田耕整及后續油麥栽培具有重要應用價值[3]。由于旋埋刀輥結構復雜、各因素存在交互作用、優化項目主次不明確、刀輥試制成本高，以及稻板田耕作時節限制等因素影響，僅通過田間試驗途徑研究刀輥減阻降耗比較困難，所以研究初期亟需構建稻板田旋耕功耗預測模型，以指導刀輥優化設計。

離散元法(Discrete element method，DEM)是一種用于分析與求解復雜離散系統動力學問題的數值方法，廣泛應用于散粒結構的模擬研究[4-5]。土壤是農業工程領域中最常見的顆粒狀物料之一，普遍存在的各向異性、非均勻及不連續等特征使土壤耕作系統更加復雜[6]。離散元法能有效解決耕作部件與土壤相互作用的非線性問題，獲取田間試驗無法測量的數據信息，對耕作部件的性能進行預測及參數優化，是當前農機具研究的重要輔助手段。國內外學者基于離散元耕作模型開展了大量研究，主要涉及耕后力學指標評估和土壤形態及細觀位移分析，其中深松鏟與土壤的接觸模型相對成熟，土壤無滑移摩擦接觸模型[7]、延遲彈性接觸模型[8-9]和粘結接觸模型[10-12]應用廣泛。為增加深松行為的模擬精度，根茬[13]以及土層差異情況[14-15]也被考慮在模型構建中。深松為土壤小擾動耕作行為，其耕作模型不適用于旋耕。而針對旋耕作業的離散元模型研究報道較少，根據文獻[16-17]，目前模型應用局限于單回轉面旋耕(1～3把旋耕刀)及土槽試驗，無法推演旋耕刀輥的田間真實作業情況。

本文根據稻板田土壤性質及旋埋刀輥旋耕作業特點，構建刀輥-土壤相互作用模型，旨在預估旋埋刀輥作業功耗，為后續旋埋刀輥減阻降耗的優化設計提供理論依據，同時獲取的功耗數據可為刀輥扭力分析提供參考。

1 稻板田旋耕模型構建

土壤接觸模型通過微觀層面描述土壤顆粒間的動力學行為以表達宏觀土壤在載荷下復雜的應力應變特性。農具作用于土壤，土壤不斷形變直至破壞失效，同時土壤將產生的反作用力以耕作阻力的形式反饋給農具。耕作過程由土壤和農具共同參與，因此為確保仿真正確度，構建離散元模型時在考慮土壤性質的同時也需要關注農具結構及耕作方式的特殊性。

1.1 土壤顆粒接觸模型選取

長江中下游稻作區土壤普遍為粘性質地，稻油、稻麥水旱輪作導致土壤常年干濕交替，對于粘土而言，土壤性質受含水率影響較大。在工程應用中，粘土因含水率表現出的軟硬稀稠程度通常由液塑限指標衡量，該指標將土壤分為固態、塑態和液態[18]。土壤的固塑液狀態關系到耕作時土壤破壞形式，固態土壤強度大，受載后易產生裂紋并擴散破碎；塑態土壤可塑性強，能承受較大塑性形變而不離散；液態土壤呈泥漿狀，幾乎不具備抵抗外力和維持原有形狀的能力。在夏秋水稻收獲季節，由于氣候和地表裸露等因素，土壤含水率處于一年中最低水平，此時粘土濕脹干縮特性導致的田面土壤裂縫是稻板田最顯著的特征。根據2017年至2019年在華中農業大學現代農業科技試驗基地中稻收獲后所得試驗檢測結果顯示，稻板田耕層土壤含水率在20.38%～28.24%之間。為進一步評估稻板田土壤破壞特征，使用五點法取田間土壤，利用LP-100D型數顯式土壤液塑限聯合測定儀(上海路達實驗儀器有限公司)，依據SL 237—1999《土工試驗規程》要求，對土壤塑限進行測量，結果如圖1所示，其中塑限均值為24.45%，液限均值為43.28%。

圖1 液塑限試驗結果Fig.1 Result of liquid plastic limit test

由圖1可知，稻板田含水率區間(20.38%～28.24%)與塑限區間(21.34%～26.57%)十分接近，與液限區間(41.67%～45.28%)差距明顯，此時土壤僅可承受較小的形變并很快破碎，因此在進行離散元建模時可將稻板田土壤作固態處理。

HertzMindlin with Bonding顆粒接觸模型[19]能夠同時反映土壤的不連續性與團聚特點，有效解決農具與土壤相互作用的非線性問題，該模型不僅繼承了HertzMindlin(no slip)的算法，而且通過在土壤顆粒間設置粘結鍵的方式，添加額外的力與力矩以約束顆粒運動，粘結鍵受載斷裂后，不能再次生成，符合稻板田土壤破碎后保持松散狀態的力學行為特征。該模型共有5個參數，其中法向、切向剛度以單位步長為間隔不斷迭代并更新粘結鍵所受的載荷(力與力矩)；法向、切向臨界應力為判斷粘結鍵是否斷裂的閾值；粘結半徑為顆粒間產生粘結鍵所需要的最大間距。在離散元軟件中，稻板田土壤顆粒由粘結狀態到失效破壞的過程模型為：

粘結鍵受載迭代原理

(1)

式中Fn、Ft——粘結鍵法向、切向受力，N

Tn、Tt——粘結鍵法向、切向扭矩，N·m

vn、vt——法向、切向速度，m/s

ωn、ωt——法向、切向角速度，rad/s

Sn、St——法向、切向剛度，N/m3

A——接觸面積，m2

J——慣性矩，m4

RB——粘結變徑，m

δt——單位時間步長，s

粘結鍵斷裂條件

(2)

式中σmax——法向臨界應力，Pa

τmax——切向臨界應力，Pa

1.2 旋耕模型參數類型及標定方案

離散元仿真參數由本征參數(土壤顆粒與旋埋刀輥材料的密度、泊松比、剪切模量)、材料接觸參數(土壤顆粒間以及土壤顆粒與旋埋刀輥間的恢復系數、靜摩擦因數、滾動摩擦因數)和接觸模型參數(HertzMindlin with Bonding顆粒接觸模型的5個參數)組成，通過借鑒文獻成果、實際試驗以及虛擬標定的方法獲取。

本征參數僅與物料的材料屬性相關，而其他仿真參數受物料尺寸、結構和接觸關系等因素影響。對于土壤而言，尺寸小而尺度范圍大，且外部輪廓多樣，粒徑配比[20]和剛性重疊球體簇[21-22]的方法雖可縮小土壤實物與模型的差距，但在離散元建模時仍不可避免忽略土壤大部分真實的微觀結構，采用放大尺寸的圓球模擬土壤逐漸成為業界普遍共識，與接觸相關的實測值對于簡化后的圓球狀土壤模型反而失去了適用性。參數虛擬標定是將仿真試驗和真實試驗相結合，通過數學方法得到仿真參數與試驗指標之間的關系，最終以仿真值逼近試驗值的形式獲取與研究物料相匹配的參數。該方法完全適用于土壤這種復雜的研究對象，并已實現沙土[23]、沙壤土[24]、粘壤土[20,25]等多種土壤模型的構建。但一種土壤模型往往側重于一種或一類土壤性質的表達，并不能包含該類土壤的所有特性，其模型參數應根據應用情況不同而有所改變，同時僅靠土壤模型無法描述土壤與應用設備間的接觸情況，因此土壤模型區別于土壤應用模型。通過文獻研究發現，現有土壤模型的構建均基于堆積角試驗，堆積角往往是由松散顆粒物料在堆積過程中受重力影響形成。測量稻板田土壤堆積角的前提需對土壤進行破碎，且破碎程度影響測量結果，獲取的接觸參數僅能表征土壤在某一破碎程度下的堆積特征。稻板田土壤被農具擾動前應被視為土粒群相互粘結的緊實整體，耕后破碎情況復雜，顯然堆積角試驗不適用于稻板田旋耕研究。

選用何種方案標定接觸參數，關系到模型的適用性。旋耕作業是由直線牽引運動與旋轉運動復合而成，刀輥刀具由入土至出土過程土壤相互接觸形式和作用形式十分復雜，大致分為入土、出土和回程3個過程，如圖2所示，刀具在入土區間與實土接觸，對土壤有剪切撕裂、擠壓和滑移摩擦作用；出土區間與松土接觸，對土壤提升拋撒并二次切削破碎土塊；回程區間刀具與拋起的土壤發生沖擊與碰撞。旋耕過程土壤力學行為的特殊性，也說明了接觸情況的復雜性，這也是深松作業離散元模型無法直接應用于旋耕作業的原因，因此除旋耕試驗本身，很難找出能直接表征旋耕過程中土壤破碎情況及土壤與刀輥接觸情況的替代方案。

圖2 旋耕作業過程示意圖Fig.2 Schematic of rotary tillage process

根據圖3可知，旋埋刀輥由6段組成，每段刀輥結構一致，沿幅寬方向縮小刀輥尺度，以1∶6縮比試制刀輥，并安裝于旋耕測試平臺內，將平臺放置在稻板田土壤試驗土槽上，其中土壤含水率23.37%，容重1 521 kg/m3。以前進速度0.6 m/s、刀輥轉速280 r/min、耕深15 cm為工況，旋耕功耗為指標進行模型接觸參數標定參照試驗，試驗過程如圖4，每次試驗后由測試平臺的鎮壓裝置將土面恢復至原有高度以確保土壤狀態的一致性，試驗共重復3次，結果見表1。

圖3 旋埋刀輥結構圖Fig.3 Structure drawing of rotary burying blade roller

圖4 模型參數標定參照試驗Fig.4 Model parameter calibration reference test

1.3 模型參數標定

1.3.1本征參數

本征參數為材料的固有特性，有相應參考資料和較為成熟的測量方法，其中刀輥材料參考了文獻[26]；土壤泊松比和剪切模量由三軸壓縮試驗結果(彈性模量5×106Pa，內摩擦角23.8°)根據文獻[23，25]的方法近似換算測得；顆粒密度通過顆粒填充試驗由土壤容重校正獲取，校正公式為

表1 旋耕功耗試驗結果Tab.1 Test result of rotary tillage power consumption kW

(3)

式中ρ——顆粒密度，kg/m3

V——容器體積，m3

k——填充顆粒數量

r——顆粒半徑，取4×10-3m

ρ1——土壤容重，kg/m3

最終確定的稻板田土壤與刀輥材料的本征參數如表2所示。

表2 模型本征參數Tab.2 Basic properties of models

1.3.2接觸參數范圍

接觸參數中包含6個材料接觸參數和5個模型接觸參數。HertzMindlin with Bonding顆粒接觸模型在應用時，常認為粘結行為源于土粒間的液橋，所以可用粘結半徑衡量土壤含水率。所以假定土壤中的水分均勻分布，且包裹在土壤顆粒外圍形成均勻水膜，水膜厚度與顆粒半徑之和即為粘結半徑，計算公式為

(4)

式中ω——土壤含水率，%

ρ2——水密度，kg/m3

根據HertzMindlin with Bonding顆粒接觸模型應用于土壤耕作方面的研究成果[14,26-31]，綜合分析后，確定其他接觸參數的取值范圍，如表3所示。

表3 待標定接觸參數取值范圍Tab.3 Range of contact parameters to be calibrated

1.3.3接觸參數標定方法及結果

通過試驗設計與數學方法令仿真值逼近測試值最終確定接觸參數是離散元參數標定的主要途徑。最陡爬坡標定方法[20,32]準確高效，在離散元參數標定的研究中得到廣泛應用。該方法首先通過Plackett-Burman試驗提取顯著因素，然后以參照試驗結果為目標值，將各顯著因素以等步長增加的方式創建各步階參數水平并開展仿真試驗，通過相對誤差變化趨勢，進一步縮小參數范圍，最后建立多因素回歸模型并求解，完成參數標定。但最陡爬坡標定方法仍存在以下問題：Plackett-Burman試驗對于因素顯著性的判別存在爭議[33]，現有文獻也僅在定性的尺度上確定因素是否顯著；當顯著因素較多時，試驗次數會成倍增加，導致回歸模型的獲取變得異常困難。

在對稻板田土壤與旋埋刀輥接觸模型參數標定時，為避免上述問題，需對最陡爬坡試驗方法進行優化。首先保留該方法中最核心的等步長爬坡思想，由于各因素以等步長增加方式進行仿真試驗，步階次序一旦確定，就間接確定該步階次序下各因素取值，因此建立步階次序與仿真指標間的函數關系，可達成指標值與各因素參數值的關聯，起到類似回歸模型的作用。

爬坡試驗的起點參數組為A1(待標定接觸參數的下限值)，終點參數組為Am(待標定接觸參數的上限值)，則各步階參數組可表示為

(5)

式(5)中m為爬坡的步階總數(即標定試驗總次數)，n為待標定參數總數量(n=10)，步長A為

(6)

步階次序x的取值范圍為1≤x≤m，對應各接觸參數組的取值為

Ax=(x-1)A+A1

(7)

所以功耗P=f(X1，X2，X3，X4，X5，X6，X7，X8，X9，X10)可由式(7)轉換為P=f(x)。

當m越大，各組標定試驗的接觸參數取值越接近，旋耕功耗隨參數組的變化更加準確，但需要仿真的次數會增加，綜合標定準確度與工作量，令m=9，建立旋埋刀輥與稻板田土壤顆粒在各步階接觸參數下的旋耕仿真模型，如圖5所示，模型中土槽尺寸(長×寬×高)為1 200 mm×500 mm×300 mm，共計容納441 185個土壤顆粒，刀輥工況與稻板田參照試驗一致，刀輥轉速280 r/min(0～1.1 s)，垂直入土速度1.50 m/s(0～0.1 s，耕深15 cm)，前進速度0.60 m/s(0.1～1.1 s)，數據以0.001 25 s為間隔保存，即刀輥旋轉一周存儲171個數據，刀輥平穩作業區間(0.1～1.1 s)旋轉4.67圈，共計儲存800個數據，從中選取刀輥運動的4個周期求功耗均值，功耗仿真結果如表4所示。

表4 仿真標定試驗結果Tab.4 Result of simulation calibration test

圖5 離散元標定試驗仿真模型Fig.5 Discrete element calibration test simulation model

隨著步階次序增加，刀輥功耗單調遞增，對仿真結果進行擬合，如圖6所示。

圖6 仿真結果擬合曲線Fig.6 Fitting curve of simulation result

擬合方程為

P=0.754 2x+2.156 7

(8)

其中決定系數R2為0.998 7，說明擬合方程對仿真試驗值的擬合程度較好。

將稻板田參照試驗結果(表1)代入式(8)，求得步階次序x=4.26，根據式(7)，最終獲得稻板田土壤與旋埋刀輥相互作用對應的接觸參數組A4.26，如表5所示。

表5 接觸參數標定結果Tab.5 Calibration results of contact parameters

將標定后的接觸參數組A4.26輸入EDEM軟件中，進行仿真模型驗證，得到旋耕功耗為5.269 kW，仿真預測值與實測值的相對誤差為1.84%，差異較小，表明在等步長爬坡試驗中利用步階次序建立模型的方法可行，因此接觸參數組A4.26能夠反映旋埋刀輥與稻板田土壤在功耗指標下的離散元旋耕接觸情況。

2 功耗預測模型驗證

2.1 模型適用性分析

為檢驗該離散元模型的適用性，增加通用刀輥(圖7)，利用旋耕平臺分別以3種工況開展稻板田土槽試驗和仿真試驗，并計算仿真值相對于實測值的相對誤差并取絕對值，結果如表6所示。

圖7 刀輥實物圖Fig.7 Physical map of blade roller

表6 誤差對比試驗設計與結果Tab.6 Error comparision test design and result

根據表6可知，預測誤差均值為6.65%，范圍在3.63%～9.48%之間；旋埋刀輥與通用刀輥的平均誤差分別為5.28%和8.02%；工況Ⅰ、Ⅱ、Ⅲ的誤差均值分別為7.43%、5.10%、7.42%。旋埋刀輥的功耗預測誤差小于通用刀輥，說明用旋埋刀輥標定的離散元參數用于其他刀輥仿真研究時，可能會造成額外的誤差，可能與模型中不同接觸參數對刀輥結構變化引起功耗指標變化的敏感程度不同有關，原結構的標定結果可能忽略了對現結構較敏感的參數，造成了誤差。以這種觀點研究不同工況下的旋耕功耗時，發現工況Ⅱ的誤差較小，相對于工況Ⅰ、Ⅲ，工況Ⅱ與標定工況差距較小，由此推斷這種額外誤差也有可能存在作業工況中。

對試驗結果進行方差分析，如表7所示，在α=0.05水平下刀輥類型和作業工況對試驗結果的影響不顯著，說明從統計學觀點可認為，該離散元模型在應用過程中誤差并非來自刀輥類型和作業工況。

表7 方差分析Tab.7 Variance analysis

綜上，功耗預測的相對誤差均在合理范圍，說明稻板田旋耕功耗預測模型基本滿足應用需求，但對于一些結構或作業方式與旋耕作業差距較大的耕作部件(如犁、深松鏟等)，該模型的適用性難以保證。

2.2 原尺度旋埋刀輥功耗預測

于2019年10月在華中農業大學現代農業科技試驗基地開展試驗。所選試驗田內留有中稻秸稈，含水率為22.16%，容重為1.52×103kg/m3。刀輥安裝在幅寬2.3 m的側邊傳動旋耕機機架內，由東方紅LX954型輪式拖拉機驅動，通過調節油門與擋位改變作業工況，功耗由動態轉矩轉速傳感器(北京中航科儀測控技術有限公司，轉速測量范圍0～4 000 r/min，扭矩測量范圍0～3 000 N·m) 測量，每種工況功耗重復3次，同時根據實測工況和表2、5的仿真參數建立稻板田旋耕功耗預測模型，進行仿真。田間試驗與仿真試驗過程如圖8所示，功耗對比結果見表8。

圖8 旋埋刀輥功耗預測模型檢驗Fig.8 Power consumption prediction test of rotary buried blade roller

表8 功耗預測結果Tab.8 Power consumption prediction results

由表8可知，田間試驗功耗預測誤差范圍在2.50%～12.81%之間，均值為7.28%，刀輥結構在縮放過程誤差變化較小，說明離散元接觸參數標定試驗方案可行。根據每種工況的3次試驗可知，極差分別為2.632、4.386、6.176 kW，與對應工況功耗均值的百分比依次為8.22%、11.59%、14.21%。田塊區域土壤不均勻、機具操控不穩定等不可控因素的干擾使實測重復誤差較大，當刀輥優化前后功耗變化小于重復誤差時，單純通過少量重復的田間試驗很難驗證優化結果，該稻板田旋耕功耗預測模型精度符合應用要求，能在一定誤差范圍內獲取刀輥功耗，當模型輸入參數變化，輸出功耗隨之響應，可為后續旋耕刀具減阻降耗的優化研究提供支持。

3 結論

(1)栽播冬油菜、冬小麥前稻板田土壤的含水率通常與其塑限接近，此時土壤不能承受較大的塑性形變，符合HertzMindlin with Bonding顆粒接觸模型的模擬要求。

(2)土壤耕作由農具和土壤共同參與，參數標定時需同時考慮土壤與農具的接觸特征。鑒于旋耕作業結構與運動的特殊性，結合旋埋刀輥的結構特點，沿幅寬方向以1∶6 縮比試制刀輥，進行標定參照試驗；稻板田土壤與旋埋刀輥旋耕接觸參數的虛擬標定采用等步長爬坡試驗，通過步階次序建立了接觸參數與功耗指標之間的函數關系；結合標定參考試驗功耗值，最終確定了稻板田旋耕功耗預測模型的接觸參數取值，完成了模型的構建。

(3)為驗證該模型的適用性，在不同工況下對通用刀輥和旋埋刀輥的旋耕功耗進行預測，預測誤差范圍為3.63%～9.48%，均值為6.65%，方差分析顯示，稻板田旋耕功耗預測模型適用于不同旋耕刀輥及工況下的功耗預測。原尺度刀輥田間試驗功耗預測誤差范圍為2.50%～12.81%，均值為7.28%，刀輥結構在縮放過程誤差變化較小，說明模型能準確反映旋埋刀輥在稻板田作業的功耗情況。