旋耕后農田模型離散元參數標定

摘要：土壤因其含水率、成分、粒徑等因素不同導致接觸參數存在差異，為提高新型高密度灌木幼苗移栽機開溝覆土作業過程的準確性，對北方春季耕后農田土壤進行離散元參數標定試驗。通過物理試驗與仿真試驗相結合的方法，采用Hertz-Mindlin（no slip）接觸模型以土壤—土壤堆積角和土壤—觸土鋼滑動摩擦角為響應值，恢復系數、靜摩擦系數和滾動摩擦系數為目標值進行參數標定。首先將目標值應用Box-Behenken響應面分析法進行三因素三水平的正交旋轉試驗得到二階回歸模型；其次以實際堆積角和滑動摩擦角為響應值，通過Design-Expert尋優模塊得到最優參數組合：土壤—土壤恢復系數0.4，土壤—土壤靜摩擦系數0.8，土壤—土壤滾動摩擦系數0.39；土壤—鋼恢復系數0.27，土壤—鋼靜摩擦系數0.53，土壤—鋼滾動摩擦系數0.16；并進行物理試驗與仿真試驗的對比，結果表明：兩者誤差分別為0.68%和1.34%，在可接受范圍之內，可信度良好，具有一定工程指導意義。

關鍵詞：農田；土壤；旋耕；離散元；參數標定

中國分類號：S152.9" " " 文獻標識碼：A" " " 文章編號：2095?5553 （2024） 08?0296?06

Calibration of discrete element parameters of farmland model after rotary tillage

Guo Lingyun1， Zhao Meiqing2， Gao Youshan1， Wang Aihong1， Ning Zhiqiang1， Lü Shining1

（1. College of Mechanical Engineering， Taiyuan University of Science and Technology， Taiyuan， 030024， China;"2. Department of Mechanical and Electronic Engineering， Shanxi University of Engineering and Technology，"Yangquan， 045000， China）

Abstract： Soil contact parameters are different" due to different factors such as soil moisture content， composition， and particle size. In order to improve the accuracy of the trenching and soil covering operation of a new high?density shrub seedling transplanter， this study conducted a discrete element parameter calibration experiment on farmland soil after spring cultivation in northern China. Through a combination of physical and simulation experiments， the Hertz-Mindlin （no slip） contact model was used to calibrate the parameters with the soil?soil response angle and the soil?steel sliding friction angle as the response values， and the recovery coefficient， static friction coefficient， and rolling friction coefficient as the target values. Firstly， the Box-Behenken response surface analysis method was applied to the target value to conduct a three factor and three level orthogonal rotation test， so as to obtain a second?order regression model. Secondly， taking the actual stacking angle and sliding friction angle as the response values， the optimal parameter combinations were obtained through the Design-Expert optimization module as follows： soil?soil recovery coefficient was" 0.4， soil?soil static friction coefficient was 0.8， and soil?soil rolling friction coefficient was 0.39. Soil?steel recovery coefficient was 0.27， soil steel static friction coefficient was 0.53， and soil steel rolling friction coefficient was 0.16. A comparison between physical experiments and simulation experiments was also conducted， and the results showed that the errors of the two methods were 0.68% and 1.34%， respectively， which were within an acceptable range， and had a good reliability and certain engineering guiding significance.

Keywords： farmland; soil; rotary tillage; discrete element; parameter calibration

0 引言

我國地域遼闊，農產品種類大不相同，華北地區作為重要的經濟農業區，以種植小麥、玉米、灌木[1]為主。隨著農業機械化不斷加快，已經出現了較為成熟的移栽機械，然而種植農作物時的耕后開溝覆土是重要工序，需要明確土壤和觸土部件之間的作用機制。大田試驗只能宏觀上了解土壤的擾動情況，無法探究到二者之間內部動態行為；又因其土壤的復雜性和不確定性，目前沒有一種固定準確的參數來代表所有的土壤模型。因此基于離散元對粒子系統進行仿真，能夠獲得離散物料難以測量的接觸參數。

近年來，離散元在工業工程中的應用越來越廣泛，除應用于玉米、冶金、土石、混凝土等散體物料外，有關接觸參數標定的研究更是數不勝數，如飼料、羊糞、玉米秸稈、種子、煤等。離散元在農業工程中的應用多數在對土壤力學參數的研究，其接觸模型和標定方法也不相同。張銳等[2]利用默認Hertz-Mindlin模型通過測量無粘性沙土的實際堆積角、剪切模量等因素對標準球和非標準球進行參數標定，進一步比較兩種球型的誤差大小。李俊偉等[3]采用Hertz-Mindlin with JKR Cohesion模型，以土壤堆積角和土球在斜面的滾動距離為響應值對東北地區兩種含水率的黏重黑土之間、觸土部件和土壤顆粒之間進行離散元仿真得到接觸參數的最優解。都鑫等[4]選擇Bonding模型，以極限破碎位移和極限破碎載荷為響應值進行肥料顆粒參數標定，并且借助外槽輪排肥器驗證最優參數組合的可靠性。王憲良等[5]利用ECM彈塑性模型以土壤堆積角、黏聚力、內摩擦角為目標量，土壤顆粒半徑、靜摩擦因數、滾動摩擦因數為自變量，通過高斯—牛頓迭代法進行參數優化，用輪胎—土壤壓實模擬試驗進行驗證。以上學者在該領域取得了很好的成果，但大部分是對原始土壤的研究，有關旋耕后農田模型離散元參數標定的資料屈指可數。

因此，本文為了保證開溝覆土[6， 7]作業的準確性，在前人的研究成果基礎上對旋耕后農田模型接觸參數進行標定。首先，搭建物理試驗平臺測得實際土壤靜堆積角和實際土壤—觸土鋼的滑動摩擦角；其次，利用Boxbehken響應面分析法對待標定系數進行多因素多水平正交旋轉試驗，對二階模型尋優后獲得回歸方程，進而得到最優參數組合；最后將最優參數再次進行仿真試驗，求出與物理試驗的誤差值，證明接觸參數模型標定的準確性和可行性。

1 接觸模型的選取

本文的研究對象是經過旋耕破碎作業后的春季農用土壤，土質柔軟細膩、疏松透氣。EDEM軟件中默認的Hertz-Mindlin（no slip）模型可以用于研究耕后土壤的力學特性。

設半徑分別為[R1]、[R2]的兩球形顆粒發生彈性接觸，法向重疊量

[α=R1+R2-r1-r2] （1）

式中： [r1]、[r2]——兩顆粒球心位置矢量。

顆粒間的接觸面為圓形，接觸半徑

[a=aR*] （2）

式中： [R*]——等效粒子半徑。

[1R*=1R1+1R2] （3）

顆粒間法向力

[Fn=43E*（R*）12a32] （4）

式中： [E*]——等效彈性模量。

[1E*=1-V12E1+1-V22E2] （5）

式中： [E1]、[E2]——顆粒1和顆粒2的彈性模量；

[V1]、[V2]——顆粒1和顆粒2的泊松比。

仿真中滾動摩擦很重要，可以通過接觸表面上的力矩來說明，即

[Ti=-μrFnRiwi] （6）

式中： [μr]——滾動摩擦因素；

[Ri]——質心到接觸點的距離；

[wi]——接觸點處物體的單位角速度矢量。

2 接觸模型參數標定

2.1 試驗材料與本征參數

本文的試驗材料為春季耕后農田土壤和觸土鋼（Q235鋼），為確保試驗的嚴謹性和準確性，通過環刀法[8]采用“蛇形步點”的大面積步點方法，“隨機、等量、多點、具有代表性”的步點原則對深度為0～200 mm的土壤進行取樣，并用環刀（100 cm3）和電子天平（0.001 kg）測量土壤樣品的平均密度。為了高計算效率，縮小簡化模型，用孔隙大小為3 mm的細篩對土壤樣品進行篩分，并用電子天平稱重出200 g的樣品作為試驗材料。觸土材料是與開溝器材料屬性相同的Q235鋼，因進行開溝作業之前已經完成旋耕工作，因此強度塑性等綜合性能較好的Q235鋼完全能滿足使用要求。試驗材料的本征參數如表1所示。

2.2 物理試驗平臺的搭建和求解

散料在堆放時能夠保持自然穩定狀態的最大角度（單邊對水平面的角度）稱為安息角，又稱靜堆積角。搭建物理試驗平臺來測量靜堆積角，如圖1（a）所示，主要部件有支架、漏斗、接料盤。漏斗總高150 mm，口外徑150 mm，口徑20 mm，嘴長45 mm，接料盤為一塊方形板，漏斗口徑到接料盤的高度為100 mm。將細篩篩選過的200 g試樣土壤在漏斗頂部進行自由落體后在接料盤形成穩定的錐狀土堆，用高清數碼相機從兩個不同的方向記錄此時的堆積狀態，重復5組。將收集的堆積狀態圖經過Matlab去噪、灰度、二值處理后用Getdate軟件提取邊緣點數據，最后Origin軟件對邊緣點數據進行線性擬合求實際堆積角[Y1]的平均值，為30.68°。

滑動摩擦角的測量由圖1（b）來完成，主要由一塊平面和長度500 mm，底寬40 mm，高20 mm，厚度為1 mm，材質為Q235鋼的U型槽組成。取同樣細篩處理過的50 g試樣土壤均勻鋪灑在U型槽中，以一端為支點，在另一端用步進電機驅動絲桿的升降裝置以10°/s勻速舉升。為了更明顯地觀察土壤的動態行為，控制步進電機以0.5 s為間歇步長，同時在步進電機非工作時間觀察土壤有無滑落趨勢2～3 s，循環此工作至土壤開始整體滑動。這時用高清數碼設備記錄此時U型槽的角度，重復5組。將獲得的狀態圖導入CAD軟件中，測量其滑動摩擦角[Y2]的平均值，為30.558°。

2.3 仿真標定試驗設計

本文以無法直接測量的土壤—土壤和土壤—觸土鋼恢復系數、靜摩擦系數、滾動摩擦系數為待標定值，以實際土壤堆積角和滑動摩擦角為響應值依據Design-Expert軟件中Box-Beheken響應面優化法進行參數標定。查閱相關文獻宋少龍[9]、戴飛[10]等的土壤離散元參數標定、Sun等[11]的仿生深松器的DEM模擬研究、Ucgul等[12]的田間條件下使用全尺寸模板犁進行表土埋置的離散元建模研究、孫偉等[13]的覆土裝置覆土研究特性可知，土壤—土壤的恢復系數[x1]、靜摩擦系數[x2]、滾動摩擦系數[x3]的待標定范圍分別在0.1～0.7、0.2～0.8、0.2～0.6；土壤—觸土鋼之間的恢復系數[x4]、靜摩擦系數[x5]、滾動摩擦系數[x6]的待標定范圍分別在0.18～0.6、0.3～0.75、0.05～0.4，因素水平編碼表如表2、表3所示。

為了避免仿真試驗失真，要保證與物理試驗平臺的一致性，如圖2所示。

在測量堆積角時，首先將SolidWorks三維模型導入EDEM2018軟件，其次采用系統默認的Hertz-Mindlin模型，并通過漏斗口外徑的顆粒工廠以50 g/s生成半徑為1.5 mm、總質量為200 g的顆粒，靜置2 s后分別從z和x方向截屏，最后測每個方向雙邊仿真堆積角的平均值。在測量滑動摩擦角時，通過顆粒工廠以100 g/s生成半徑為1.5 mm、質量為50 g的顆粒，然后一端為支點以10° /s順時針轉動，當試樣土壤整體開始滑動時停止仿真并記錄測量此時U型槽的角度。

2.4 仿真試驗結果及分析

2.4.1 二階回歸模型

二次正交旋轉仿真試驗結果如表4、表5所示。

對表4、表5進行回歸模型方差分析如表6、表7所示。二者回歸模型的P值均lt;0.01，表示回歸模型差異極顯著，失擬項P值均大于gt;0.05，表示回歸模型失擬項不顯著；[x1]、[x2]、[x3]，[x4]、[x5]、[x6]中至少有兩項P值lt;0.05且[x12、x22、x32]，[x42、x52、x62]均lt;0.01；對于表6的判定系數[R2]=0.977且校正判定系數Adjusted [R2]=0.947 4；表7的判定系數[R2]=0.991 4且校正判定系數Adjusted [R2]=0.980 3，兩者的值無限逼近于1，說明回歸模型擬合程度都較好；在表6中預測判定系數Predicted [R2]=0.822 8，|[Adjusted R2-Predicted R2]|lt;0.2，說明回歸模型預測效果可靠；變異系數CV=4.2lt;10，說明回歸模型可信度高；Adeq Precision＝19.026 6gt;4，說明有良好的精確度。在表7中預測判定系數Predicted [R2]=0.894 1，|[Adjusted R2-Predicted R2]|lt;0.2，說明回歸模型預測效果可靠；變異系數CV=3.07lt;10，說明回歸模型可信度高；Adeq Precision＝30.051 8gt;4，說明有良好的精確度。

由此可得，兩個回歸模型均可預測恢復系數、靜摩擦系數、滾動摩擦系數三者與靜堆積角和滑動摩擦角的關系。

表6中回歸項[x1]、[x2]、[x12]影響極顯著，[x3]、[x1、x2]影響顯著，其余回歸項影響不顯著，由F值大小可知，影響靜堆積角的顯著順序為：[x1]gt;[x2]gt;[x3]，即土壤—土壤恢復系數gt;土壤—土壤靜摩擦系數gt;土壤—土壤滾動摩擦因數。得靜堆積角[Y1]的回歸方程為

[Y1=27.26+4.74x1+2.85x2-0.97x3+1.59x1x2-0.005x1x3+0.86x2x3-3.15x12-1.14x22-1.03x32] （7）

在表7中回歸項[x5]、[x6]、[x5x6]、[x52]影響極顯著，其余回歸項影響不顯著，由F值大小可知，影響靜堆積角的顯著順序為：[x5]gt;[x6]gt;[x4]，即土壤—觸土鋼靜摩擦系數gt;土壤—觸土鋼滾動摩擦系數gt;土壤—觸土鋼恢復系數。得滑動摩擦角[Y2]的回歸方程為

[Y2=30.99+0.1288x4+8.93x5+1.61x6+0.095x4x5-0.0475x4x6+1.72x5x6+0.505x42-1.75x52-0.7075x62] （8）

2.4.2 交互項響應面模型分析

在Design-Expert軟件中的Analysis-Model Graphs模塊中生成3D Surface，宏觀分析交互項對堆積角和滑動摩擦角的影響趨勢（圖3、圖4）。

由圖3（a）可知，當恢復系數小于0.4時，隨著靜摩擦系數增大，等高線變化較平緩，說明此時靜摩擦系數對堆積角影響不大，當恢復系數不小于0.4時，恢復系數和靜摩擦系數響應面曲線接近橢圓且走勢相近，說明兩個系數對堆積角的影響基本相同；由圖3（b）可知，當恢復系數小于0.4時，等高線愈密集愈趨于平緩，說明此時響應面愈陡，滑動摩擦系數對堆積角影響不顯著，當恢復系數不小于0.4時，則相反；由圖3（c）可知，當靜摩擦系數小于0.3時，靜摩擦系數與滾動摩擦系數響應面曲線走勢相同，當靜摩擦系數0.3～0.5時，等高線趨于平緩，說明此時滾動摩擦系數對堆積角影響不顯著。

由圖4（a）可知，當靜摩擦系數0.3～0.7時，隨著恢復系數的變化響應面曲線大致呈平行式遞增，說明靜摩擦系數對滑動摩擦角影響更顯著；由圖4（b）同理可得，當滾動摩擦系數小于0.2時，恢復系數對堆積角的影響不顯著；由圖4（c）可知，當靜摩擦系數小于0.6時，靜摩擦系數對滑動摩擦角影響顯著，之后滾動摩擦系數對滑動摩擦角影響逐漸加大。

2.4.3 參數優化及模型驗證

利用Design-Expert軟件中Optimization-Numerical模塊，將土壤靜堆積角平均值[Y1]=30.68、滑動摩擦角平均值[Y2]=30.558分別代入式（7）、式（8）優化求解獲得最優參數組合為：[x1]=0.4，[x2]=0.8，[x3]=0.39；[x4]=0.27，[x5]=0.53，[x6]=0.16。即土壤—土壤恢復系數為0.4、土壤—土壤靜摩擦系數為0.8、土壤—土壤滾動摩擦系數為0.39；土壤—觸土鋼恢復系數為0.27、土壤—觸土鋼靜摩擦系數為0.53、土壤—觸土鋼滾動摩擦系數為0.16。最后利用上述最優參數組合依照2.3節中的仿真標定試驗設計進行5次重復仿真試驗，并按照2.2節中物理試驗平臺的搭建和求解的測量方法，得到土壤—土壤堆積角平均值為30.[47°]，土壤—觸土鋼滑動摩擦角平均值為30.[97°]，與物理堆積角、滑動摩擦角的誤差分別為0.68%，1.34%，在可接受范圍之內，試驗對比如圖5所示，進一步證明仿真試驗有效。

3 結論

1） 因土質存在復雜性和不確定性，為確定北方農田土壤經過旋耕作業后的接觸力學參數，本文采用Hertz-Mindlin接觸模型以物理試驗與仿真試驗相結合的研究方法進行離散元參數標定。搭建物理試驗平臺，通過多次堆積角試驗和滑落試驗獲得實際堆積角和滑動摩擦角。

2） 確定土壤—土壤與土壤—觸土鋼的恢復系數、靜摩擦系數、滾動摩擦系數的取值范圍后，利用Design-Expert軟件中的Box-Beheken三因素三水平正交旋轉試驗獲得二階回歸模型，并對該模型進行方差分析和響應面分析以確定其準確性和合理性。

3） 以物理堆積角和滑動摩擦角為響應值，回代二階回歸模型后獲得最優參數組合：土壤—土壤的恢復系數為0.4、靜摩擦系數為0.8、滾動摩擦系數為0.39；土壤—觸土鋼的恢復系數為0.27、靜摩擦系數為0.53、滾動摩擦系數為0.16。將兩組最優參數進行仿真驗證試驗，與物理堆積角和滑動摩擦角對比后誤差分別為0.68%和1.34%，誤差在可接受范圍之內，驗證仿真模型參數標定的真實性和可靠性。

參 考 文 獻

[ 1 ] 張瑀桐. 華北地區主要糧食作物生長水足跡及適水種植研究[D]. 北京： 中國水利水電科學研究院， 2019.

Zhang Yutong. Research on the growth water footprint and water?suitable planting of major grain crops in north China [D]. Beijing： China Institute of Water Resources and Hydropower Research， 2019.

[ 2 ] 張銳， 韓佃雷， 吉巧麗， 等. 離散元模擬中沙土參數標定方法研究[J]. 農業機械學報， 2017， 48（3）： 49-56.

Zhang Rui， Han Dianlei， Ji Qiaoli， et al. Calibration methods of sandy soil parameters in simulation of discrete element method [J]. Transaction of the Chinese Society for Agricultural Machinery， 2017， 48（3）： 49-56.

[ 3 ] 李俊偉， 佟金， 胡斌， 等. 不同含水率黏重黑土與觸土部件互作的離散元仿真參數標定[J]. 農業工程學報， 2019， 35（6）： 130-140.

Li Junwei， Tong Jin， Hu Bin， et al. Calibration of parameters of interaction between clayey black soil with different moisture content and soil?engaging component in northeast China [J]. Transaction of the Chinese Society of Agricultural Engineering， 2019， 35（6）： 130-140.

[ 4 ] 都鑫， 劉彩玲， 姜萌等. 基于離散元的包膜肥料Bonding模型參數標定[J]. 農業機械學報， 2022， 53（7）： 141-149.

Du Xin， Liu Cailing， Jiang Meng， et al. Parameter calibration of coated fertilizer Bonding model based on discrete element [J]. Transactions of the Chinese Society for Agricultural Machinery， 2022， 53（7）： 141-149.

[ 5 ] 王憲良， 胡紅， 王慶杰， 等. 基于離散元的土壤模型參數標定方法[J]. 農業機械學報， 2017， 48（12）： 78-85.

Wang Xianliang， Hu Hong， Wang Qingjie， et al. Calibration method of soil contact characteristic parameters based on DEM theory [J]. Transaction of the Chinese Society for Agricultural Machinery， 2017， 48（12）： 78-85.

[ 6 ] 馬小英， 杜伊健， 閆玉濤， 等. 觸土推板對土壤動態接觸行為的離散元分析[J]. 沈陽理工大學學報， 2022， 41（3）： 66-71.

[ 7 ] 戴飛， 張仕林， 宋學鋒， 等. 全膜雙壟溝雙幅覆膜覆土聯合作業機設計與試驗[J]. 農業機械學報， 2020， 51（5）： 108-117.

Dai Fei， Zhang Shilin， Song Xuefeng， et al. Design and test of combined operation machine for double width filming and covering soil on double ridges [J]. Transaction of the Chinese Society for Agricultural Machinery， 2020， 51（5）： 108-117.

[ 8 ] 袁娜娜. 室內環刀法測定土壤田間持水量[J]. 中國新技術新產品， 2014（9）： 184.

[ 9 ] 宋少龍， 湯智輝， 鄭炫， 等. 新疆棉田耕后土壤模型離散元參數標定[J]. 農業工程學報， 2021， 37（20）： 63-70.

Song Shaolong， Tang Zhihui， Zheng Xuan， et al. Calibration of the discrete element parameters for the soil model of cotton field after plowing in Xinjiang of China [J]. Transaction of the Chinese Society of Agricultural Engineering， 2021， 37（20）： 63-70.

[10] 戴飛， 宋學鋒， 趙武云， 等. 全膜雙壟溝覆膜土壤離散元接觸參數仿真標定[J]. 農業機械學報， 2019， 50（2）： 49-56， 77.

Dai Fei， Song Xuefeng， Zhao Wuyun， et al. Simulative calibration on contact parameters of discrete elements for covering soil on whole plastic film mulching on double ridges [J]. Transaction of the Chinese Society for Agricultural Machinery， 2019， 50（2）： 49-56， 77.

[11] Sun J， Wang Y ， et al. DEM simulation of bionic subsoilers （tillage depth gt;40 cm） with drag reduction and lower soil disturbance characteristics [J]. Advances in Engineering Software， 2018， 119： 30-37.

[12] Ucgul M， Saunders C， et al. Discrete element modelling of top soil burial using a full scale mouldboard plough under field conditions [J]. Biosystems Engineering， 2017， 160： 140-153.

[13] 孫偉， 劉小龍， 石林榕， 等. 刮板升運帶式膜上覆土裝置覆土特性[J]. 機械工程學報， 2016， 52（7）： 38-45.

Sun Wei， Liu Xiaolong， Shi Linrong， et al. Covering soil on plastic?film characteristics of scraper lifting belt mechanism [J]. Journal of Mechanical Engineering， 2016， 52（7）： 38-45.