葡萄藤防寒土與清土部件相互作用的離散元仿真參數標定

馬 帥，徐麗明，袁全春，牛 叢，曾 鑒，陳 晨，王爍爍，袁訓騰

葡萄藤防寒土與清土部件相互作用的離散元仿真參數標定

馬 帥，徐麗明※，袁全春，牛 叢，曾 鑒，陳 晨，王爍爍，袁訓騰

（中國農業大學工學院，北京 100083）

為系統地研究中國北方地區沙壤土質地的葡萄藤防寒土及其與清土機清土部件常用材料（Q235鋼、橡膠）相互作用的離散元仿真參數，以構建準確的土壤離散元仿真模型，該文選用整合延遲彈性模型（hysteretic spring contact model，HSCM）和線性粘附模型（liner cohesion model，LCM）作為土壤顆粒間的接觸模型；基于土壤堆積試驗，以土壤顆粒間恢復系數、靜摩擦系數、滾動摩擦系數和土壤粘附能量密度為因素，以土壤堆積角為指標，利用EDEM進行通用旋轉中心組合模擬試驗，采用Design-Expert軟件對試驗數據進行回歸分析，以實測的土壤堆積角作為優化目標值，獲得土壤顆粒間的最佳接觸參數組合；利用土壤屈服試驗獲得HSCM模型參數；基于斜面滑動法原理，利用傾斜板試驗臺測得土壤與Q235鋼和橡膠之間的靜摩擦系數，并以此為基礎，采用土壤滑落試驗，以滑動摩擦角為響應值，對土壤顆粒與Q235鋼和橡膠之間的恢復系數和滾動摩擦系數進行尋優，得到最優解參數組合。為驗證標定優化的離散元模型參數的準確性，采用刮土板土槽試驗和仿真試驗進行對比分析，獲得刮土板在土槽試驗和仿真試驗中的水平前進阻力分別為228.36 N和213.79 N，兩者之間的相對誤差為6.38%，表明仿真模型中土壤的物理力學特性與實際土壤基本一致，驗證了葡萄藤防寒土離散元仿真參數標定結果和研究方法準確可靠。研究結果可為基于離散元法研制適用于北方地區沙壤土質地的葡萄藤防寒土清土機提供理論基礎和技術支撐。

土壤；離散元法；清土機；葡萄藤；防寒土；參數標定

0 引 言

中國葡萄種植區域廣泛，其中北方葡萄產區由于冬季氣候寒冷干燥，需在入冬前將修剪后的葡萄藤壓倒在地面進行埋土防寒作業，等到第二年春天氣候變暖，再把防寒土清除，葡萄藤上架[1-3]。其中，清土作業工作量大、勞動強度高、機械化程度低[4-5]，為促進葡萄產業的全程機械化發展，未來加大防寒土清土機的相關研制是必然趨勢，而防寒土清土機作業時，清土部件與土壤直接接觸，對于不同類型的土壤，其作業效果差異顯著，故全面系統地研究葡萄藤防寒土離散元模型仿真參數，有助于為基于離散元法設計及優化防寒土清土機提供理論依據。

土壤離散元模型仿真參數主要包括材料本征參數、顆粒之間及顆粒與觸土部件之間的接觸參數、接觸模型參數[6-8]。由于一般應用離散元法（discrete element method，DEM）建立的土壤顆粒模型尺寸和形狀與實際土壤顆粒不同，因此，在運用離散元法對防寒土清土機與土壤相互作用仿真之前，需先根據實際土壤的質地類型，選擇合適的土壤接觸模型，并對其接觸參數和接觸模型參數進行標定[9-10]。

目前，國內外學者利用離散元法對土壤與農機部件的相互作用方面做了大量研究，對土壤等散粒體的物料特性參數標定做了大量工作。Ucgul等[11-12]通過結合Hertze-Mindlin及Hysteretic Spring接觸模型，分析了土壤間有粘結力及無粘結力時的情況，同時考慮了土壤的塑性形變，為粘結力不大的土壤選擇接觸模型提供了參考；石林榕等[13]將線性內聚力模型添加到延遲彈性模型中，建立西北旱區農田土壤模型，并對不同含水率的土壤進行了離散元仿真參數的標定，該模型不僅可以體現土壤的塑性變形，而且能體現土壤間的粘結力特性；張銳等[14]應用默認的Hertz-Mindlin 接觸模型，利用堆積角標定土壤顆粒在標準球和非標準球狀態下相互作用的參數，建立了仿真模型，該模型只針對于沙壤土使用；王憲良等[15]基于代理模型的方法，利用The Edinburgh Elasto-Plastic Cohesion Model（ECM）彈塑性接觸模型表示土壤顆粒接觸模型，并對顆粒間的靜摩擦系數和滾動摩擦系數進行了標定，該模型是考慮了土壤間彈塑性變形的一種非線性模型；李俊偉等[16]利用EDEM中的JKR模型對不同含水率黏重黑土顆粒之間的接觸模型仿真參數及其與不同觸土部件之間的接觸參數進行了標定，該模型適應于含水率較高的土壤，一般水田中應用廣泛；向偉等[17]基于土壤堆積試驗對南方黏壤土的離散元仿真參數進行標定，構建了南方黏壤土精準的離散元仿真模型；丁啟朔等[18]利用Hertz-Mindlin with Bonding模型建立了黏性水稻土的深松耕作離散元模型，該模型中顆粒間具有粘結鍵，適用于顆粒間粘結力特別大的土壤，尤其適用于混凝土和礦石等硬質顆粒。

本文針對中國北方地區沙壤土質地的葡萄藤防寒土（以寧夏、甘肅地區沙壤土質地的葡萄園為主），選用整合延遲彈性模型（hysteretic spring contact model，HSCM）和線性粘附模型（liner cohesion model，LCM）作為土壤顆粒間的接觸模型，采用仿真試驗與物理試驗相結合的方法，對葡萄藤防寒土離散元模型參數及其與清土機中清土部件常用材料（Q235鋼、橡膠）的接觸參數進行標定優化，并采用清土作業中常用的刮土板進行土槽試驗與仿真試驗的對比分析，對標定的優化仿真參數進行驗證，以期獲得較為準確的防寒土離散元模型參數。

1 土壤顆粒仿真模型參數標定

1.1 土壤顆粒間接觸模型的選取

不同的接觸模型使用不同的方程計算顆粒之間的接觸力和阻尼力，故針對不同的土壤質地類型，須建立不同的接觸模型。本文以中國北方地區（以寧夏、甘肅地區為主）沙壤土質地的葡萄園為主，入冬前，采用埋土機從行間取土，覆蓋在葡萄藤上，形成一個橫截面近似等腰梯形的土壟，由于經過整個冬季的沉積和雨雪，第二年春天清除時，防寒土具有一定的堅實度，平均含水率約為6%~8%，土壤顆粒之間既表現出明顯的散粒體物料特性，又有一定的粘附力，同時當土壤顆粒受到外力壓縮作用時會產生塑性形變。

在離散元軟件EDEM內嵌的接觸模型中，延遲彈性接觸模型（hysteretic spring contact model，HSCM）可用于顆粒產生塑性形變的場合；線性粘附接觸模型（liner cohesion model，LCM）可用于一般性粘結顆粒的快速計算，亦可用于具有一定含水率的顆[12]。為了真實反映沙壤土質地類型的葡萄藤防寒土宏觀特性，本文選用整合延遲彈性模型和線性粘附模型作為土壤顆粒間的接觸模型。由于假設摩擦力限制了HSCM方程中切向粒子的運動，因此將粘附力增加到法向接觸力中[19]，粘附力F/a(N)的大小為

式中為顆粒間凝聚能力密度，J/m3；A為顆粒間接觸面積，m2。A的計算公式如下：

式中r為顆粒間接觸區域的半徑，m，定義如下：

式中r為接觸顆粒等效半徑，m；E為顆粒模型等效楊氏模量，MPa；F為法向接觸力，N。法向接觸力F的計算公式如下：

式中1為加載剛度，N/m；2為卸載剛度，N/m；δ為顆粒間法向重疊量，m；0為顆粒間殘余重疊量，m。此外，顆粒間法向阻尼力的計算公式為

接觸顆粒間的總法向力為粘附力、法向接觸力和阻尼力之和，即

式中F為總法向力，N；F為法向接觸力，N；F為阻尼力，N；F為粘附力，N。

1.2 仿真模型參數標定

1.2.1 本征參數

土壤顆粒本征參數包括顆粒形狀尺寸、密度、泊松比和剪切模量，這些參數是固定的，一般從測量結果和其他參考文獻中獲取。在離散元模擬中，當涉及大量土壤作業即生成的顆粒數量非常多時，建模時生成與實際土壤顆粒尺寸大小相同的顆粒模型是不切實際的，并且若采用任意形狀顆粒進行建模，將會成倍增加仿真計算時間[20]，故為了建模的方便性和計算的快速性，本文采用離散元軟件EDEM中默認的球形顆粒，顆粒基體半徑設定為5 mm，球徑倍率分布為0.95～1.05[13]；供試土壤取自寧夏賀蘭山地區沙壤土質地的葡萄園內，其含水率在6%～8%之間，采用容積為2×10-4m3的環刀和量程為3 000 g、精度為0.01 g的電子天平測量供試土壤的密度，經10次重復測量，取其平均值1 389 kg/m3作為最終的仿真參數；泊松比和剪切模量依據文獻[13,21]分別取0.35和2.89×1010Pa。

Q235鋼與橡膠是防寒土清土機中常用的清土部件材料，根據文獻[8]，Q235鋼的密度、泊松比和剪切模量分別為7 850 kg/m3、0.28和8.2×1010Pa，橡膠的密度、泊松比和剪切模量分別為960 kg/m3、0.45和3 448 Pa。

1.2.2 接觸參數標定

土壤顆粒間接觸參數的標定基于堆積試驗，標定之前，采用漏斗法測定供試土壤的堆積角，測試裝置如圖1a所示，主要由支架、漏斗（上開口直徑為195 mm，下開口直徑為80 mm，高度為100 mm）、土壤堆積鋼板、鋼尺、數字傾角儀（分辨率：0.05°，精度：±0.2°）等組成。土壤從漏斗中緩慢滑下，在鋼板上堆積穩定后，采用鋼尺和數字傾角儀配合測量此時的土壤堆積角，試驗重復10次取平均值，最終測試結果為31.48°。

采用整合HSCM和LCM模型的土壤顆粒進行堆積仿真試驗，如圖1b所示，采用隨機分布方式在漏斗上方開口位置生成土壤顆粒，待生成的所有顆粒從漏斗中落下并停止運動時，模擬結束，利用EDEM后處理中自帶的Protractor配合Clipping工具分別從錐面的方向和方向測定土壤堆積角，取其平均值作為最終的仿真試驗結果。

1.數字傾角儀2.鋼板3.土壤4.鋼尺5.漏斗6.支架

根據1.1節選取的土壤顆粒接觸模型，在堆積角仿真試驗中，以土壤顆粒間的恢復系數、靜摩擦系數、滾動摩擦系數和LCM模型中的粘聚能量密度作為待標定參數[22]，以土壤堆積角1作為試驗指標，進行4因素通用旋轉中心組合模擬試驗，根據試驗設計原則，確定編碼系數為2.0，結合仿真預試驗結果和文獻[23-24]，確定堆積角仿真試驗因素的編碼結果如表1所示，共實施31組試驗（0水平試驗重復7次），試驗方案與結果如表2所示。

表1 土壤堆積角仿真試驗因素編碼表

應用Design-expert 8.0軟件對仿真試驗結果進行顯著性分析，結果如表3所示。由分析結果可知，在4個試驗因素中，對土壤堆積角影響顯著，對土壤堆積角影響不顯著，和對土壤堆積角影響極顯著；交互項因素中只有項對土壤堆積角影響極顯著，其余因素的交互影響均不顯著；二次方項因素中，2對土壤堆積角影響顯著，2對土壤堆積角影響極顯著，其余影響不顯著。回歸模型的值小于0.001，而失擬項的值大于0.05，說明回歸模型極其顯著且失擬不顯著，回歸有效。將不顯著項刪除后得到土壤堆積角1的回歸方程如式（7）。

表3 回歸模型顯著性分析

注：**表示極顯著（<0.01）；*表示顯著（0.01<<0.05），下同。

Note: ** Means highly significant (<0.01); * means significant (0.01<<0.05), the same below.

1=-20.75+137.62+182.41+3.23×10-3-5.93×10-3

-129.172-189.672（7）

采用Design-expert 8.0軟件Optimization-Numerical模塊，以實測的土壤堆積角31.48°為目標，對堆積角回歸模型進行尋優，得到多組優化解，經仿真驗證試驗，選取與物理試驗得到的土壤堆積角最接近的一組，即土壤間恢復系數為0.51，靜摩擦系數為0.65，滾動摩擦系數為0.06，粘聚能量密度為10 495 J/m3，此優化解下對應的堆積角仿真試驗值為31.74°，與物理試驗結果的相對誤差為0.83%，無顯著差異，說明所得的土壤顆粒間接觸參數最優值準確可靠。從圖2可以看出，采用優化參數進行仿真試驗得到的土壤堆積錐形與實際物理試驗得到的土壤錐形具有很高的相似性。

圖2 土壤堆積角仿真與物理試驗對比

1.2.3 接觸模型參數測定

對于本文選用的HSCM和LCM模型，其中LCM模型的粘聚能量密度等于土壤內聚強度，已在1.2.2節中進行了標定；HSCM模型的屈服強度通過土壤屈服試驗測得，測試裝置如圖3所示，主要由貫入圓軸（軸長150 mm、直徑15 mm，底端帶有厚度5 mm、直徑20 mm的圓盤）、100 mm×100 mm×100 mm的空心鐵盒、REGER萬能試驗機、底部支撐裝置和頂部夾持裝置等部分組成。顆粒間的阻尼系數和剛度系數采用EDEM中的默認值，分別為0.05和0.95[13,25-26]。

1.貫入圓軸 2.空心鐵盒 3.底部支撐 4.土壤 5.頂部夾持裝置 6.REGER萬能試驗機

試驗前，根據空心鐵盒容積大小，用電子天平測量1.39 kg供試土壤，裝入空心鐵盒并將頂部土壤壓平，然后放在萬能試驗機上，貫入圓軸由頂部夾持裝置固定，試驗時，設置萬能試驗機向下運行速度為200 mm/min，在夾持裝置的帶動下，貫入圓軸勻速插入土壤中，在與萬能試驗機連接的計算機上保存得到貫入壓力隨插入位移的變化值，將壓力值除以圓軸底端圓盤的水平橫截面面積，即可得到貫入應力隨插入位移的變化曲線，如圖4所示。起始階段，隨著貫入圓軸插入深度的增加，貫入應力不斷增大，當貫入圓軸的底端插入土壤中達到37 mm左右時，貫入應力開始下降，土壤發生失效，此時土壤受到的應力值為0.38 MPa，37～50 mm之間應力值趨于平緩，隨著貫入圓軸繼續下降，由于其底端不斷接近鐵盒底部，貫入圓軸本身開始受力變形，從而使應力值迅速上升。因此，HSCM模型中的屈服強度為0.38 MPa。

圖4 貫入應力與位移之間的關系圖

2 土壤與清土部件材料接觸參數的標定

針對清土機清土部件的常用材料Q235鋼和橡膠，需測量和標定土壤顆粒與其相互的接觸參數。土壤與機械部件材料之間的接觸模型為Hertz-Mindlin (no slip)，采用靜摩擦試驗測定土壤顆粒與2種材料之間的靜摩擦系數，并基于靜摩擦系數，采用土壤滑落試驗標定土壤顆粒與2種材料之間的恢復系數和滾動摩擦系數。

2.1 土壤靜摩擦試驗

2.1.1 靜摩擦試驗原理

斜面滑動法是測量靜摩擦系數的常用方法[16]，原理如圖5所示，質量為的物體放置在傾角為的斜面上，其重力分解為2個力：沿斜面向下的分力×sin和垂直于斜面的分力×cos。當斜面傾角小于滑動臨界角時，×sin小于物體與斜面間的靜摩擦力，物體保持靜止，隨著斜面傾角的增加，×sin越來越大，當大于物體滑動臨界角時，×sin＞，物體將開始沿著斜面下滑。基于圖5中物體的受力平衡分析，通過公式（8）～（11），可計算得到兩者之間的靜摩擦系數。

式中為靜摩擦系數。

注：為物體質量，kg；為斜面傾角，(°)；F為斜面對物體的支持力，N；為斜面對物體的摩擦力，N；為重力加速度，m·s2。

Note:is the mass of the object, kg;is the slope angle, (°);Fis the support force of slope on object, N;is the friction force of slope on object, N;is the acceleration of gravity, m·s2.

圖5斜面滑動法測量原理

Fig.5 Measuring principle of slope sliding method

2.1.2 土壤靜摩擦試驗

利用自制的傾斜試驗臺進行靜摩擦試驗，如圖6所示，主要由底板、傾斜板、土壤擋板、鐵絲和升降臺等組成。其中4個土壤擋板安裝在傾斜板上，形成一個無蓋土槽容器，在土槽容器中裝滿供試土壤并將其頂部土壤壓平，以消除由于表面起伏對被測清土部件材料滑動的影響，傾斜板通過鐵絲連接在升降裝置上，通過升降裝置的上下運動帶動傾斜板繞其末端轉動軸旋轉。

1.Q235鋼板 2.土壤擋板 3.底板 4.傾斜板 5.數字傾角儀 6.土壤 7.鐵絲 8.升降臺9.橡膠片

由于土壤與不同材料之間的靜摩擦系數是土壤與材料之間的特有屬性，不會因為材料形狀和大小的改變而改變，故試驗時將預先制備好的被測清土部件材料樣品（50 mm×50 mm×3 mm的Q235鋼板、50 mm×50 mm× 3 mm的橡膠片）放在土壤表面頂端，將數字傾角儀校正后放置在傾斜板上，設定升降裝置向上運行的速度為200 mm/min，在鐵絲的牽引下傾斜板開始緩慢傾斜，當被測觸土材料樣品在土壤表面開始滑動的瞬間，關停升降裝置，記下此時傾角儀的度數，每個試驗重復10次取平均值，測量結果如表4所示，依據公式（11）計算靜摩擦系數，即tan20.81°=0.38，tan25.44°=0.48，土壤與Q235鋼和橡膠之間的靜摩擦系數分別為0.38和0.48。

表4 Q235鋼板和橡膠片在土壤表面滑動瞬間的傾角

2.2 土壤滑落試驗

2.2.1 試驗方法

針對沙壤土顆粒，土壤與清土部件材料（Q235鋼和橡膠）之間的恢復系數和滾動摩擦系數不易直接測定，基于已測定的靜摩擦系數，利用土壤滑落試驗進行標定，為更精確地觀察、控制試驗狀態并測定相應的試驗結果，利用土壤顆粒在傾斜板上完全滑落時，對應滑動摩擦角的大小作為標定試驗依據值。測定裝置采用2.1.2節的自制傾斜試驗臺，當測量土壤顆粒在不同的觸土材料上滑動時，只需在傾斜板上土槽容器底部更換放置對應的觸土材料即可。經前期預試驗可知，選用土壤質量的多少與其在清土部件材料上的滑動行為及完全滑落時對應滑動摩擦角大小無關，本文選擇100 g供試土壤進行試驗。試驗時，將土壤堆放在觸土材料表面的頂部，如圖7a所示，然后設定升降裝置向上運行的速度為200 mm/min，傾斜板緩慢傾斜，當土壤從材料表面的頂部完全滑落時，如圖7b所示，關停升降裝置，利用傾角儀測定此時傾斜板傾斜的角度即為滑動摩擦角的大小，每個試驗重復10次取平均值。

圖7 土壤在橡膠上的滑落試驗

對應的土壤滑落仿真試驗如圖8所示，在EDEM軟件中利用顆粒工廠生成同等質量的土壤顆粒并堆放在觸土材料表面的頂部，如圖8a所示，然后設定傾斜板繞其底部旋轉軸開始轉動，當土壤顆粒完全滑落時，如圖8b所示，利用EDEM后處理中自帶的Protractor工具測量傾斜板的傾角，此時傾斜板的傾角即為滑動摩擦角的仿真測定結果。

圖8 土壤在橡膠上的滑落仿真試驗

2.2.2 試驗設計

以土壤顆粒與Q235鋼板和橡膠之間的恢復系數和滾動摩擦系數為試驗因素，以土壤在2種材料上完全滑落時對應的滑動摩擦角2和3為評價指標，分別進行二因素通用旋轉中心組合仿真模擬試驗。參考文獻[24,27-28]，確定仿真試驗因素水平的編碼結果如表5所示，共實施13組響應面分析試驗（零水平試驗重復5次），試驗方案與結果如表6所示。

表5 土壤滑落仿真試驗因素編碼表

表6 土壤滑落仿真試驗方案與結果

2.2.3 試驗結果分析與參數優化

應用Design-expert 8.0軟件對仿真試驗結果進行回歸顯著性分析，結果如表7所示。

由分析可知，在2個試驗因素中，土壤與材料之間的滾動摩擦系數對土壤在鋼板和橡膠上的滑動摩擦角的影響極顯著，而土壤與材料之間的恢復系數對土壤在鋼板上的滑動摩擦角影響不顯著，對土壤在橡膠上的滑動摩擦角影響顯著。土壤在鋼板和橡膠上滑動摩擦角的回歸模型的值均小于0.01，而失擬項的值均大于0.05，說明回歸模型極其顯著且失擬不顯著，回歸有效。將不顯著項刪除后得到土壤在鋼板上的滑動摩擦角2、土壤在橡膠上的滑動摩擦角2的回歸方程如式（12）～（13）。

表7 回歸模型顯著性分析

2=20.23+13.13（12）

3=5.08+11.37+119.99?214.402（13）

2.2.1節的試驗測得土壤顆粒在鋼板和橡膠上完全滑落時的滑動摩擦角分別為25.30°和28.10°。通過Design-expert 8.0軟件Optimization- Numerical模塊，分別以實測的2個滑動摩擦角為目標值進行尋優，得到多組優化解，經仿真驗證試驗，選取與實測試驗得到的滑動摩擦角最接近的一組，即土壤與鋼板間的恢復系數為0.60、滾動摩擦系數為0.37，結合2.1.2節測得的靜摩擦系數0.38，此優化解下對應的土壤顆粒在鋼板上完全滑落時的滑動摩擦角為25.40°，與實測結果的相對誤差為0.40%；土壤與橡膠間的恢復系數為0.61、滾動摩擦系數為0.23，結合2.1.2節測得的靜摩擦系數0.48，此優化解下對應的土壤在橡膠上完全滑落時的滑動摩擦角為28.80°，與實測結果的相對誤差為2.49%；由此表明，標定優化后的土壤顆粒與鋼板和橡膠之間的接觸參數準確可靠。

3 刮土板刮土驗證試驗

3.1 試驗目的

標定的土壤離散元模型參數與土壤顆粒間的粘結力、內聚力、塑性形變等密切相關，選用不同的土壤模型參數建模，土壤所表現出來的宏觀性質是不同的，其中表現較明顯的是觸土部件在土壤中運動時受到的阻力大小。采用刮土板進行刮土作業時，土壤在刮土板的作用下運動，不僅可以反映土壤顆粒模型間的相互作用，而且可以反映觸土部件與土壤顆粒模型間的相互作用，故本文為驗證標定優化后參數的可靠性和準確性，選用葡萄藤防寒土清除作業中常用的刮土板進行刮土試驗，以作業時刮土板的水平前進阻力為響應值，將土槽試驗中的實測值與EDEM中的仿真值進行對比分析，利用相對誤差值判斷所構建的葡萄藤防寒土離散元模型的有效性。

3.2 刮土板刮土土槽試驗

3.2.1 試驗條件準備

試驗在中國農業大學工學院土槽中進行，為保證試驗時土槽中土壤條件與實際葡萄園內的葡萄藤防寒土條件一致，試驗前，通過人工制備土壤，包括旋耕、灑水、壓實、翻運和起壟等處理，并根據春天清土時土壟的外形尺寸，制成了橫截面近似于等腰梯形的2條土壟，每條長約22 m，土壟上底寬約25 cm，下底寬約70 cm，高約25 cm。用SC-900 型土壤緊實度儀測定土壟緊實度,從頂部到底部約在69～135 kPa之間，土壟底部緊實度大于頂部緊實度，用TDR 150型土壤含水率測定土壟的含水率約在6%～8%之間。

3.2.2 土槽試驗

本文設計的刮土板機構主要由機架、刮土板、角度調節板和裝在刮土板底部的橡膠等部分組成，通過三點懸掛系統掛接在土槽臺車上，如圖9a所示，刮土厚度通過升降系統調節，以土壟上表面為基準，設定刮土厚度為12 cm，試驗時，通過操作控制臺上的計算機，設定土槽臺車的前進速度為1 km/h，數據保存頻率為20 Hz，土槽臺車帶動刮土板前進并進行刮土作業，如圖9b所示，將由六分力測力架測得的刮土板在作業時的水平前進阻力實時存儲在計算機上。共進行2次平行刮土試驗，并將2次平行試驗測得的刮土板水平前進阻力的均值記為最終測試結果。

1.土槽 2.刮土板 3.角度調節板 4.橡膠 5.機架 6.土壟 7.六分力測力架 8.臺車 9.控制柜 10.計算機

3.3 刮土板刮土仿真試驗

運用EDEM軟件建立土槽模型，設置其基本尺寸（長′寬′高）為1 800 mm′700 mm′270 mm，采用標定優化后的土壤離散元模型參數在土槽中生成總量為220 000個土壤顆粒模型，形成1 800 mm長的土壟模型，然后在土壟模型頂部設置一個向下運動的壓板壓實土壤，使建立的土壟離散元模型的緊實度和橫截面尺寸與土槽試驗的參數一致，然后將在SolidWorks中建立的刮土板三維模型按照1:1的比例導入EDEM中[29]，保證刮土板模型與土壟模型之間的相對位置與土槽試驗時一致，試驗裝置材料為Q235鋼和橡膠，其參數見1.2.1和2.2.3節，最終建立的EDEM仿真模型如圖10a所示。

依據土槽試驗條件，設定刮土板的刮土厚度為120 mm，水平前進速度為1 km/h，沿軸正方向。在EDEM求解器模塊對仿真時間步長、仿真時間、數據保存間隔及網格大小進行設置[30]。為保證仿真的連續性，設置其仿真時間步長為1.40×10-6s（即Rayleigh 時間步長的40%），仿真時間為7 s，網格單元尺寸設置為顆粒平均半徑的4倍，刮土板在EDEM中的仿真作業狀態如圖10b所示，待仿真作業完成后，在后處理中輸出刮土板沿軸方向受到的阻力大小。

1.土壟模型 2.機架 3.刮土板 4.橡膠

3.4 試驗結果與分析

土槽試驗中，隨著臺車的前進，刮土板開始接觸土壤，受到的水平阻力由0開始增大，以土壟起點為0點，當臺車前進7.5 m左右時，刮土板與土壤穩定接觸，阻力變化趨于平緩，如圖11a所示，隨著刮土板的前進，刮土板受到的阻力波動較大，分析原因可知，是由于人工起壟使土壟在不同位置的形狀尺寸具有一定的誤差以及土壤條件的差異造成。采用水平阻力的平均值作為最終的試驗結果，通過計算其結果為228.36 N。

分析仿真試驗可知圖11b，刮土板作業時受到的水平前進阻力的整體變化趨勢與土槽試驗中大體一致，隨著刮土板進入仿真計算區域與土壤顆粒模型開始接觸，水平阻力不斷增大，在1.5 s左右時，刮土板開始與土壤穩定接觸，阻力變化趨于平緩，5.3 s左右以后，由于被刮土板壅積在最前方的土壤開始離開仿真計算區域，水平阻力值不斷減小，水平阻力在1.5～5.3 s之間平緩變化，取其平緩階段的阻力平均值作為最終的仿真試驗結果，即為213.79 N。

圖11 刮土板水平前進阻力的變化曲線

由此可知，刮土板在土槽試驗中受到的水平前進阻力與仿真試驗結果的相對誤差為6.38%，誤差在可接受范圍內，表明了仿真模型中土壤的物理力學特性與實際土壤一致，驗證了葡萄藤防寒土離散元仿真參數標定結果和研究方法準確可靠。

4 結 論

針對中國北方地區沙壤土質地的葡萄藤防寒土，選用整合HSCM和LCM作為土壤顆粒間接觸模型，采用實測試驗與仿真試驗相結合的方法，對土壤離散元模型參數及其與清土機清土部件常用材料（Q235鋼、橡膠）的接觸參數進行系統標定，確定了最優參數組合，并通過土槽試驗和仿真試驗的對比分析進行驗證，得到如下主要結論：

1）在EDEM中，以土壤顆粒間恢復系數、靜摩擦系數、滾動摩擦系數和土壤粘附能量密度為因素，以土壤堆積角為指標，設計了通用旋轉中心組合模擬試驗，對試驗數據進行回歸分析，表明只有土壤間靜摩擦系數對土壤堆積角無影響，其他接觸參數均有顯著影響；以土壤堆積角實測試驗結果作為優化目標值，利用響應面優化方法獲得土壤間接觸參數的最優組合：恢復系數0.51、靜摩擦系數0.65、滾動摩擦系數0.06、粘聚能量密度10 495 J/m3；利用萬能試驗機，采用土壤屈服試驗測得HSCM模型參數屈服強度為0.38 MPa.

2）基于斜面滑動法原理，利用自制的傾斜試驗臺直接測定土壤與Q235鋼和橡膠之間的靜摩擦系數分別為0.38和0.48，并以此為基礎，采用土壤滑落試驗，對土壤與Q235鋼和橡膠間的恢復系數和滾動摩擦系數進行標定優化，得到土壤與Q235鋼之間的恢復系數和滾動摩擦系數分別為0.60和0.37，土壤與橡膠之間的恢復系數和滾動摩擦系數分別為0.61和0.23.

3）為驗證所標定優化的離散元模型參數的準確性，采用刮土板的土槽試驗和仿真試驗進行對比分析，獲得刮土板在土槽試驗和仿真試驗中的水平前進阻力分別為228.36和213.79 N，兩者之間的相對誤差為6.38%，誤差在可接受范圍內，表明仿真土壤模型的物理力學特性與實際土壤基本一致，驗證了葡萄藤防寒土離散元仿真參數標定結果和研究方法準確可靠。研究結果可為后期基于離散元法研制相關葡萄藤防寒土清土機提供理論基礎和技術支撐。

[1]徐麗明，李超，王文斌，等. 葡萄生產機械化技術與裝備[J]. 新疆農機化，2012(5)：24－26.

[2]袁全春，徐麗明，馬帥，等. 我國北方葡萄冬季埋土和春季清土作業的機械化現狀與思考[J]. 中外葡萄與葡萄酒，2017(6)：66－67.

[3]馬帥，徐麗明，邢潔潔，等. 葉輪旋轉式葡萄藤埋土單邊清除機研制[J]. 農業工程學報，2018，34(23)：1－10.

Ma Shuai, Xu Liming, Xing Jiejie, et al. Development of unilateral cleaning machine for grapevineburied by soil with rotary impeller[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(23): 1－10. (in Chinese with English abstract)

[4]周偉彬. 葡萄起藤機的設計[D]. 銀川：寧夏大學，2017.

Zhou Weibin. Design of Grape Vine Digging Machine[D]. Yinchuan: Ningxia University, 2017. (in Chinese with English abstract)

[5]郭家穩. 淺議哈密地區葡萄埋起藤全程機械化試驗推廣[J]. 農民致富之友，2017(8)：215－216.

[6]Janda A, Ooi J Y. DEM modeling of cone penetration and unconfined compression in cohesive solids[J]. Powder Technology, 2016, 293: 60－68.

[7]袁全春，徐麗明，邢潔潔，等. 機施有機肥散體顆粒離散元模型參數標定[J]. 農業工程學報，2018，34(18)：21－27.

Yuan Quanchun, Xu Liming, Xing Jiejie, et al. Parameter calibration of discrete element model of organic fertilizer particles for mechanical fertilization[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(18): 21－27. (in Chinese with English abstract)

[8]劉文政，何進，李洪文，等. 基于離散元的微型馬鈴薯仿真參數標定[J]. 農業機械學報，2018，49(5)：125－135.

Liu Wenzheng, He Jin, Li Hongwen, et al. Calibration of simulation parameters for potato minituber based on EDEM[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(5): 125－135. (in Chinese with English abstract)

[9]李博，陳軍，黃玉祥. 機械與土壤相互作用的離散元仿真研究進展[J]. 農機化研究，2015(1)：217－222.

Li Bo, Chen Jun, Huang Yuxiang. Research progress of discrete element simulation of soil-machine Interaction[J]. Journal of Agricultural Mechanization Research, 2015(1): 217－222. (in Chinese with English abstract)

[10]Milkevycha V, Munkholmb L J, Chen Y, et al. Modelling approach for soil displacement in tillage using discrete element method[J].Soil and Tillage Research, 2018, 183: 60－71.

[11]Uchul M, Fielke J M, Saunders C. 3D DEM tillage simulation: Validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil[J]. Soil and Tillage Research, 2014, 144: 220－227.

[12]Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling (DEM) of tillage: Accounting for soil cohesion and adhesion[J]. Biosystems Engineering, 2015, 129(1): 298－306.

[13]石林榕，趙武云，孫偉. 基于離散元的西北旱區農田土壤顆粒接觸模型和參數標定[J]. 農業工程學報，2017，33(21)：181－187.

Shi Linrong, Zhao Wuyun, Sun Wei. Parameter calibration of soil particles contact model of farmland soil in northwest arid region based on discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(21): 181－187. (in Chinese with English abstract)

[14]張銳，韓佃雷，吉巧麗，等. 離散元模擬中沙土參數標定方法研究[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]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 49－56. (in Chinese with English abstract)

[15]王憲良，胡紅，王慶杰，等. 基于離散元的土壤模型參數標定方法[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]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(12): 78－85. (in Chinese with English abstract)

[16]李俊偉，佟金，胡斌，等. 不同含水率黏重黑土與觸土部件互作的離散元仿真參數標定[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]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(6): 130－140. (in Chinese with English abstract)

[17]向偉，吳明亮，呂江南，等. 基于堆積試驗的黏壤土仿真物理參數標定[J]. 農業工程學報，2019，35(12)：116－123.

Xiang Wei, Wu Mingliang, Lü Jiangnan, et al. Calibration of simulation physical parameters of clay loam based on soil accumulation test[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(12): 116－123. (in Chinese with English abstract)

[18]丁啟朔，任駿，Belal E，等. 濕粘水稻土深松過程離散元分析[J]. 農業機械學報，2017，48(3)：38－48.

Ding Qishuo, Ren Jun, Belal E, et al. DEM analysis of subsoiling process in wet clayey paddy soil[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 38－48. (in Chinese with English abstract)

[19]Zeng Z, Chen Y, Zhang X. Modelling the interaction of a deep tillage tool with heterogeneous soil[J]. Computers and Electronics in Agriculture, 2017, 143: 130－138.

[20]Shmulevich I. State of the art modeling of soil-tillage interaction using discrete element method[J]. Soil and Tillage Research, 2010, 111(1): 41－53.

[21]王金武，唐漢，王金峰，等. 懸掛式水田單側修筑埂機數值模擬分析與性能優化[J]. 農業機械學報，2017，48(8)：72－80.

Wang Jinwu, Tang Han, Wang Jinfeng, et al. Numerical analysis and performance optimization experiment on hanging unilateral ridger for paddy field[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(8): 72－80. (in Chinese with English abstract)

[22]鄭侃. 深松旋耕作業次序可調式聯合作業機研究[D]. 北京：中國農業大學，2018.

Zheng Kan. Study on Subsoiling and Rotary Tillage Operation Adjustable Combine Machine[D]. Beijing: China Agricultural University, 2018. (in Chinese with English abstract)

[23]張銳，韓佃雷，吉巧麗，等. 離散元模擬中沙土參數標定方法研究[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]. Transactions of the Chinese society for Agricultural Machinery, 2017, 48(3): 49－56. (in Chinese with English abstract)

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

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]. Transactions of the Chinese society for Agricultural Machinery, 2019, 50(2): 49－56. (in Chinese with English abstract)

[25]DEM Solution Ltd. EDEM 2.5 theory reference guide[R/OL]. 2014-12-05[2015-06-07].http://www.docin.com/p-980174717.html

[26]胡國明. 顆粒系統的離散元素法分析仿真[M]. 武漢：武漢理工大學出版社，2010.

[27]趙淑紅，王加一，陳君執，等. 保護性耕作擬合曲線型深松鏟設計與試驗[J]. 農業機械學報，2018，49(2)：82－92.

Zhao Shuhong, Wang Jiayi, Chen Junzhi, et al. Design and experiment of fitting curve subsoiler of conservation tillage[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(2): 82－92. (in Chinese with English abstract)

[28]孫偉，劉小龍，石林榕，等. 刮板升運帶式膜上覆土裝置覆土特性[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. (in Chinese with English abstract)

[29]方會敏，姬長英，Farman A C，等. 基于離散元法的旋耕過程土壤運動行為分析[J]. 農業機械學報，2016，47(3)：22－28.

Fang Huimin, Ji Changying, Farman A C, et al. Analysis of soil dynamic behavior during rotary tillage based on distinct element method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(3): 22－28. (in Chinese with English abstract)

[30]Di Renzo A, Di Maio F. An improved integral non-linear model for the contact of particles in distinct element simulations[J]. Chemical Engineering Science, 2005, 60(5): 1303－1312.

Calibration of discrete element simulation parameters of grapevine antifreezing soil and its interaction with soil-cleaning components

Ma Shuai, Xu Liming※, Yuan Quanchun, Niu Cong, Zeng Jian, Chen Chen, Wang Shuoshuo, Yuan Xunteng

(100083,)

Grapevine in seasonally frozen regions needs to be warm-insulated by soil in winter with the antifreezing soil removed in spring most mechanically by a soil removal machine. The purpose of this paper is to simulate the interaction between the insulating soil (with sandy loam texture) and the soil-cleaning materials (Q235 steel, rubber) commonly used in the soil removal machine, based on the discrete element method. The simulation model was constructed based on properties of the soil by integrating the hysteretic spring contact model (HSCM) and the linear cohesion model (LCM) as the contact model between soil particles. We took soil-soil restitution coefficient, soil-soil frictional coefficient, soil-soil rolling coefficient and soil cohesion energy density as the determinants and the soil accumulation angle as an evaluation index. The 4-factor universal rotation center combination simulation test, based on the EDEM, was used to regress the relationship between the determinants and the index using the Design-Expert software. The results showed that the soil-soil frictional coefficient was the only parameter that did not have significant effects on the soil accumulation angle. The best contact-parameter combination between soil particles was obtained by taking the physically measured soil accumulation angle as the optimization objective, which gave 0.51 for the soil-soil restitution coefficient, 0.65 for the soil-soil frictional coefficient, 0.06 for the soil-soil rolling frictional coefficient, and 10 495 J/m3for the soil cohesion energy density. The associated soil accumulation angle was 31.74o, with a relative error of 0.83% compared with the physically measured results. The universal testing machine for soil yield test was used to obtain the HSCM model parameters based on the change in penetration stress with the displacement, and the resultant soil yield strength was 0.38 MPa. The static frictional coefficient between soil and Q235 steel as well as the rubber measured by the inclination test bench was 0.38 and 0.48 respectively. These data and the EDEM were used to conduct the simulation test of the soil slip, with the restitution coefficient and the rolling frictional coefficient between soil and the materials taken as the determinants and the sliding frictional angle as evaluating index. Regressing the test date with the two-factor universal rotation center combination test showed that the rolling frictional coefficient between soil and the materials had a significant effect on the sliding frictional angle between soil and the steel plate and rubber. In contrast, the restitution coefficient between the soil and the materials impacted significantly on the sliding frictional angle between the soil and the rubber but not on the soil and the steel plate. The restitution coefficient and the rolling frictional coefficient between the soil and the Q235 steel as well as the rubber were optimized by using the measured sliding frictional angle as the optimization objective. The resultant restitution coefficient and the rolling frictional coefficient were 0.60 and 0.37 respectively for the soil and the Q235 steel, and 0.61 and 0.23 respectively for the soil and the rubber. Soil bin test and simulation test of the scraper were conducted to verify the calibrated discrete element model parameters. The horizontal forward resistance of the scraper in the soil bin test and simulation test was 228.36 and 213.79 N respectively, with a relative error of 6.38%. The results presented in this paper have important implications for using discrete element method to analyze the removal of grapevine-insulating soil.

soils; discrete element method; soil-cleaning machine; grapevine; antifreezing soil; parameters calibration

馬 帥，徐麗明，袁全春，牛 叢，曾 鑒，陳 晨，王爍爍，袁訓騰. 葡萄藤防寒土與清土部件相互作用的離散元仿真參數標定[J]. 農業工程學報，2020，36(1)：40－49.doi：10.11975/j.issn.1002-6819.2020.01.005 http://www.tcsae.org

Ma Shuai, Xu Liming, Yuan Quanchun, Niu Cong, Zeng Jian, Chen Chen, Wang Shuoshuo, Yuan Xunteng. Calibration of discrete element simulation parameters of grapevine antifreezing soil and its interaction with soil-cleaning components[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(1): 40－49. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.01.005 http://www.tcsae.org

2019-09-04

2019-12-19

現代農業產業技術體系建設專項資金資助（CARS-29）

馬 帥，博士生，主要從事生物生產自動化研究。Email：13637070719@163.com

徐麗明，教授，博士生導師，主要從事生物生產自動化技術和裝備研究。Email：xlmoffice@126.com

10.11975/j.issn.1002-6819.2020.01.005

S152.9; S22

A

1002-6819(2020)-01-0040-10