新疆農田粉土離散元仿真參數標定

李清超，鄭 炫，劉進寶，楊懷君，王子龍，張魯云

（1.石河子大學機械電氣工程學院，新疆石河子 832000；2.新疆農墾科學院機械裝備研究所，新疆石河子 832000）

0 引言

【研究意義】新疆地區部分土壤屬于粉土，土壤顆粒主要以粉粒為主，比重占到80%以上[1]，由于粒徑較小，再加上新疆地區覆膜滴灌技術使用使得秋季農田土壤含有一定水分，具有一定黏性，不利于犁地作業[2～3]。新疆地區粉土土壤農田不適合免耕作業，耕翻成為基礎性措施之一。目前對于新疆地區粉土土壤的犁地作業未能有準確描述，難以形成犁地翻垡作業的理論性指導，其主要的關鍵問題是缺少準確的參數模型進行仿真模擬，獲取更優的土壤翻耕離散元模型，對機械設計人員有實際意義?！厩叭搜芯窟M展】離散元法被廣泛應用在農業機械觸土部件與土壤相互作用的研究上[4～8]，耕地模型及土壤與部件模型的準確性是保證結果有效的前提。建立一個系統、準確的新疆粉土土壤模型及粉土土壤與觸土部件模型，是利用離散元分析方法對耕整地機械進行優化設計的基礎。有研究通過物理試驗與土壤仿真試驗結合的方法，實現了對土壤離散元接觸參數的標定[9～10]。Liang 等[11]以實測的土壤堆積角和土-膜的靜摩擦系數為響應值，通過B-B試驗，對土-土之間以及土-膜之間接觸模型參數做了標定。戴飛等[12]用Hertz-Mindlin 無滑動模型，通過土壤堆積角、土壤與鋼滑動摩擦角的試驗方法，對全膜雙莖溝覆膜土壤進行了標定。謝方平等[13]通過單軸密閉壓縮試驗和無側限抗壓強度試驗，標定了基于Edinburgh Elasto-Plastic Adhesion（EEPA）模型的土壤離散元參數?！颈狙芯壳腥朦c】新疆粉土土壤的待耕地，在常年滴水灌溉和耕翻作用下，其土壤耕性與其他地區不同土壤類型的耕地主要有阻力大、易與觸土部件黏附等區別，仿真模型也會存在差異，難以通用。需要技術上重新標定，以校驗用于耕地機械設計制造的依據。【擬解決的關鍵問題】以新疆地區粉土土壤的待耕地為研究對象，離散元仿真試驗與物理試驗相結合的方法，采用標定土壤與機具材料之間及土壤與土壤之間的離散元接觸模型參數，為優化農業機械觸土部件提供參考。

1 材料與方法

1.1 材料

1.1.1 土壤本構參數

土壤樣品來源于2021年秋季作物收獲后，采自新疆石河子市小麥屯村農場（44°21′31″N，85°55′59″E），前茬作物為棉花，土壤質地為粉土。采用5點取樣法取耕翻作業0 ～30 cm深度的土壤，利用環刀（100 cm2）、YH-M50002電子天平（浙江東陽英衡智能設備廠，精度0.01 g）、DHG-9070A干燥機（上海一恒科技有限公司）、ZJ型應變控制式直剪儀（南京土壤儀器廠）、BT-2001 型激光粒度分布儀（濕法），對土壤的密度、含水率、泊松比、土壤粒徑進行測量。密度、含水率、粒徑的測量按照土壤試驗標準進行；土壤泊松比的測量方法參照文獻[14]。土壤本構參數：密度1 480 kg/m3，含水率14.5%，泊松比0.37，粒徑分布＜0.001 mm、0.001 ～0.05mm和＞0.05mm分別占1.82%，85.55%和12.63%。計算泊松比。

表1 土壤本構參數Table 1 Soil constitutive parameters

1.1.2 土壤與65 Mn靜摩擦系數

土壤與65 Mn之間的靜摩擦系數是兩者之間的特有屬性，不會因形狀改變而改變。土壤與65 Mn 之間的靜摩擦系數由靜摩擦試驗臺測得。將土壤用方刀制成1 cm×1 cm×1 cm的土塊，放于65 Mn 斜面之上，旋轉推桿來逐漸提高斜面角度，通過斜面儀來記錄土塊下滑瞬間的角度，由公式（2）計算得到靜摩擦系數，重復10次，取平均值作為最終結果，得到土壤與65 Mn 之間的靜摩擦系數為0.56。圖1

圖1 實際物理試驗Fig.1 Actual physical test

式中：為靜摩擦系數，為土塊下滑瞬間角度。

1.1.3 土壤滾動距離

用65Mn 板搭建的斜坡試驗臺，標定土壤與65 Mn 之間的接觸模型參數。土球顆粒由3D 打印的球形模具制得，半徑為5 mm。在預試驗的基礎上，取斜面的高度為30 mm，斜面長度為100 mm，水平面長度為400 mm。將土球放于斜面上端中間位置處，使其沿斜面自然滾下，至土球停止運動。用直尺測量土球沿水平方向的滾動距離，重復10 次，平均值為227.29 mm，作為最終的目標值。

1.1.4 土壤堆積角

土壤堆積角試驗裝置，主要有萬能試驗機（AI-7000LAU30，高鐵檢測儀器有限公司）、65 Mn圓筒、65 Mn板組成。圓筒[15]內徑為90 mm，高為270 mm，下方與65 Mn 板（500 mm × 500 mm）接觸，內部裝填土壤試樣1 kg。采用萬能試驗機，以0.05 m/s 的速度向上垂直升起圓筒，土壤向下滑落在65 Mn板上，待土壤靜置穩定后，用相機在垂直的方向上拍照，得到土樣堆積角的正視圖。后續使用matlab圖像識別功能，提取邊界，通過線性回歸來獲得土樣堆積角的大小。重復10 次取平均值，得到土壤的堆積角為37.64°。圖1

1.1.5 犁耕模型

選取顆粒形狀為EDEM 軟件自帶的球狀顆粒，半徑為5 mm。采用EDEM 的Geometries 功能，生成2 500 mm×3 600 mm×400 mm 的土槽，顆粒總數為3617032個，將用solidworks繪制好的犁體曲面模型保存成x_t 格式并導入到EDEM中。前犁體和后犁體設置耕寬為450 mm，耕深為300 mm，水平前進速度分別設置為8、9 和10 km/h，總仿真時間為2.3 s，其中前0.3 s 生成土壤，后2s犁體模型完成在土壤中的運動。圖2

圖2 犁體仿真耕作過程Fig.2 Plow simulation tillage process.

1.2 方法

針對常用的耕翻機具材料65 Mn，采用靜摩擦試驗測定土壤與65 Mn 材料的靜摩擦系數，運用斜坡物理試驗與斜坡EDEM仿真試驗相結合的方法來標定土壤與65 Mn鋼之間的離散元接觸模型參數。運用堆積角物理試驗與堆積角EDEM仿真試驗相結合的方法建立粉土土壤待耕地的離散元仿真模型。應用Design-Expert 軟件，根據Central Composite 試驗原理建立粉土土壤的斜坡滾動距離和堆積角與相關參數的回歸模型，并找到最優參數組合。在最優參數組合下，比較堆積角的仿真值與實際值，以及耕翻仿真試驗與耕翻田間試驗的對照，驗證模型的準確性。

1.2.1 接觸模型選取

秋季待耕期間，由于膜下滴灌的影響，粉土土壤待耕地含有一定水分，顆粒間存在明顯的黏附現象，且在耕翻時，土壤較容易黏附在觸土部件上。Hertz-Mindlin with JKR 接觸模型克服了經典的Hertz接觸模型不涉及黏結力的缺點，通過表面能來表現顆粒間的粘結力，適用于有明顯黏結和團聚的濕黏物料[16～19]。該模型的JKR 法向彈性接觸力計算公式為：

式中FJKR：JKR法向彈性接觸力，N；

E*：等效彈性模量（Pa）；

R*：等效接觸半徑（m）；

α：兩顆粒切向重疊量（m）；

γ：表面能（J/m3）；

δ：兩顆粒法向重疊量（m）。

其中，等效彈性模量和等效接觸半徑由式（5）和（6）定義：

式中E1、E2：兩顆粒的彈性模量（Pa）；

V1、V2：兩顆粒的泊松比；

R1、R2：兩顆粒的接觸半徑（m）。取表面能γ為0時，此時JKR模型與Hertz接觸模型一致，即：

即使兩顆粒沒有直接接觸，該模型也會在一定距離范圍內，存在黏結力，顆粒間存在黏結力的最大間隙為：

式中，δC和αC分別為顆粒間存在粘結力的法向和切向最大間隙。

當顆粒間隙大于0，小于δC時，兩者間黏結力達到最大值，計算公式為：

JKR接觸模型中的表面能，能較接近地模擬土壤顆粒之間以及土壤與觸土部件之間地黏結力，選用Hertz-Mindlin with JKR 接觸模型來建立新疆粉土土壤類型待耕地的離散元仿真模型。

1.2.2 EDEM仿真試驗

應用離散元軟件EDEM 2018 版本建立仿真模型，所有仿真試驗的土壤顆粒形狀選為球狀，半徑設置為5 mm，生成方式為隨機分布，三維結構模型的材料均為65 Mn。離散元仿真試驗中的土壤和65 Mn 的本構參數，通過測量或查閱參考文獻[20]獲得。表2

表2 土壤與65Mn的本構參數Table 2 Constitutive parameters of soil and 65Mn

斜坡仿真試驗模型與物理試驗保持一致，將solidworks 繪制的斜坡模型保存成x_t 格式，導入到EDEM軟件中，坐標原點位置設置在斜坡頂端，中間位置處，并在此位置設置顆粒工廠，顆粒生成方式為dynamic，顆粒數量為1，數據保存間隙為0.01 s，重力加速度為9.81 m/s2。顆粒球順著斜坡滑下，到靜止時，將水平坐標導出，再減去95.39 mm（斜坡水平長度），即為顆粒的水平滾動距離。表3

表3 土-65Mn仿真試驗因素水平編碼Table 3 Factor level coding table of soil-65Mn simulation test

參考文獻[14～15，20～22]的接觸參數選取范圍，采用二次回歸正交旋轉中心組合試驗的方法，確定土壤與65 Mn接觸參數的范圍。表4

表4 土-土仿真試驗因素水平編碼表因素水平Table 4 Factor level coding table of soil-soil simulation test

將solidworks 繪制的堆積角試驗模型保存成x_t 格式，導入到EDEM 軟件中，圓筒和65 Mn 板各尺寸與物理試驗保持一致。在圓筒上方創建顆粒工廠，生成顆粒方式為dynamic，數量為1 kg，待落在圓筒中的顆粒穩定后，圓筒以0.05 m/s 的速度垂直向上運動，顆粒在重力作用下向下滑落，最終顆粒停止運動形成堆積角，測量。圖3

圖3 仿真試驗Fig.3 simulation test.

參考文獻[14～15，20～22]的接觸參數選取范圍，采用二次回歸正交旋轉中心組合試驗的方法，確定土壤接觸參數的范圍。表5

表5 斜坡仿真試驗設計及結果Table 5 Slope simulation test design and results

1.2.3 阻力對比

田間試驗在新疆石河子市小麥屯村農場（44°21′31″N，85°55′59″E）進 行，土 壤 含 水 率為15.51%，土壤堅實度為2152.8 kPa。使用的儀器有：雷沃歐M904 拖拉機、翻轉犁、遙測儀（Autobona，黑龍江省農業機械工程科學研究院制），直尺、卷尺、秒表等。量取地塊中間100 m作為測試區間，在保證相應耕作速度的前提下，通過遙測儀測得耕作阻力，并對區間內的試驗結果取平均值。

2 結果與分析

2.1 土壤與65Mn接觸參數標定

研究表明，二次回歸方程如下：

該回歸模型P＜0.000 1，失擬項P=0.132 0 ＞0.05，決定系數R2=0.983 5，校正系數Radj2=0.968 1，均＞0.9，模型具有良好的準確度，該模型可較好地表達滾動距離與土-65 Mn 恢復系數、土-65 Mn 靜摩擦系數、土-65 Mn 動摩擦系數、土-65 MnJKR表面能之間的關系。4個因素對滾動距離影響的主次順序為土-65 Mn 動摩擦系數、土-65 MnJKR 表面能、土-65Mn 恢復系數、土-65 Mn 靜摩擦系數。其中，C、D對滾動距離影響極顯著，A、B對滾動距離影響不顯著。一次交互作用均對滾動距離影響不顯著。二次交互作用中，A2、B2、D2對滾動距離影響顯著。C2對滾動距離影響極顯著。選定靜摩擦系數為0.56 的一組最優解。選擇土球在65 Mn 板上滾動距離相近的優化解。表6，表7

表6 斜坡仿真試驗方差Table 6 Variance analysis of slope simulation test

表7 土-65Mn接觸模型參數Table 7 Parameters of soil-65mn contact model

2.2 土壤與土壤接觸參數標定

研究表明，建立堆積角Y2與土-土恢復系數X1、土-土靜摩擦系數X2、土-土動摩擦系數X3、土-土JKR表面能X4之間的二階回歸模型。表8

表8 土-土參數標定試驗設計與結果Table 8 Design and results of soil-soil parameter calibration test

模型P值小于0.000 1，失擬項P值大于0.05，決定系數R2均大于0.9，具有良好的準確度，此模型可較好地表達堆積角與土-土恢復系數、土-土靜摩擦系數、土-土動摩擦系數、土-土JKR 表面能之間的關系。其中，有4個回歸項影響極顯著（P＜0.01），分別為X4、X1X4、X3X4、X22，2個回歸項對模型影響顯著（P＜0.05），分別為X2、X42。表9

表9 土-土參數標定試驗方差Table 9 Analysis of variance of soil-soil parameter calibration test

取堆積角37.64°為目標，選擇與目標值相近的優化解。選擇軟件給出的多組最優解中的第一個組解作為土-土接觸模型的參數，最終確定土-土接觸模型參數。表10

表10 土-土接觸模型參數Table 10 Parameters of soil-soil contact model

2.3 物理試驗驗證

研究表明，土壤參數標定后的滾動距離、堆積角的相對誤差分別為6.05%、1.28%。標定后堆積角的形狀與實際形狀更為接近。表11，圖4

圖4 土壤堆積角仿真試驗與物理試驗對比Fig.4 Comparison between soil accumulation angle simulation test and physical test.

表11 試驗結果對比Table 11 comparison of test results

2.4 犁體阻力試驗驗證

研究表明，仿真的滾動距離比實際滾動距離大了6.05%，堆積角小了1.28%，在耕速10 km/h情況下，犁體耕作阻力小了6.49%；土壤含水率在田間試驗時，與取土階段相比發生了輕微的變化，也是造成誤差的原因之一。所建立的新疆粉土土壤待耕地離散元仿真模型，在誤差范圍內是可以替代土壤的實際參數的，適合秋季進行耕翻的干旱地區。表12，圖5

表12 阻力對比Table 12 Resistance comparison

圖5 田間試驗Fig.5 field experiment

3 討論

研究標定的是新疆地區易于粘附的粉土土壤離散元模型，雖然與東北地區黑土的離散元模型選用的接觸模型是一致的[20]，但具體接觸參數標定結果不同，新疆農田粉土的JKR 表面能要明顯小于東北黑土；北方地區的葡萄藤防寒土所選用的基礎模型為整合延遲彈性模型和線性粘附模型，與研究的接觸模型有一定區別[21]；與新疆棉田耕后土壤離散元模型相比，研究所選的接觸模型更適合有一定粘附特性的土壤[22]。研究與黑土區玉米秸稈-土壤混料離散元模型相比，只考慮土-土之間以及土-65Mn 之間的接觸參數，并未考慮殘茬的相關作用[15]。

4 結論

4.1 利用二次回歸擬合方法獲得土-65 Mn之間最優接觸參數組合：恢復系數0.51、靜摩擦系數0.56、動摩擦系數0.08、JKR表面能4.12。

4.2 土-土之間最優接觸參數組合為：恢復系數0.57、靜摩擦系數0.65、動摩擦系數0.23、JKR表面能4.49。

4.3 對比仿真試驗和實際試驗的滾動距離被堆積角，誤差分別為：6.05%、1.25%。

4.4 犁體8、9、10 km/h 作業速度下，仿真試驗和實際試驗的阻力相對誤差為：5.69%、5.95%、2.74%，誤差在可接受范圍內，驗證了標定的新疆粉土土壤離散元接觸參數的準確性。