基于離散元法的立式旋耕刀工作參數分析與優化*

扈偉昊，楊發展，趙國棟，林云龍，黃珂

(青島理工大學機械與汽車工程學院，山東青島，266520)

0 引言

立式旋耕與傳統臥式旋耕相比具有更深的作業深度、更大的作業幅寬和更快的作業速度，其可有效的破碎土壤、混勻前茬、降低犁底層高度、提高土壤蓄水保墑能力，具有耕后地表平整且不打亂耕層結構、一次耕作即可滿足土壤播前要求、有利于搶農時等優點[1]，這些優點使得立式旋耕在國家提倡保護性作業的背景下應用越來越多。立式旋耕作業過程中，立式旋耕刀的工作參數如前進速度、刀具轉速、耕深對自身受力和作業效果有很大影響。

國內學者關于立式旋耕方面的研究多集中于立式旋耕機結構、刀具排列及耕后土壤性狀對作物的影響等方面，而對立式旋耕工作參數的研究較少。劉芳建等就立式旋耕機的整機架構和刀片結構參數進行了分析，得到了不產生抗土現象時刀具內折彎角的臨界值；張敏等[2]發現旋耕刀72°等分角安裝方式與常規90°等分安裝方式相比運動軌跡相似，但受到的沖擊載荷更小，整機受力更平穩，使用壽命更長；聶勝委等[3]研究了立式旋耕對土壤緊實度的影響，通過試驗發現，立式旋耕能降低小麥拔節期、灌漿期10 cm的土壤緊實度，同時立式旋耕處理的小麥連續兩年產量增幅都較大，氮肥的農學效率和偏生產力均有顯著提升。離散元法(discrete element method, DEM)是一種適用于不連續介質力學行為預測的方法，其最初應用于巖石邊坡的運動分析，經過數十年的發展已廣泛應用于各個領域[4-6]。國內外學者利用離散元法對農機作業過程進行了廣泛的研究[7]。例如，熊平原等[8]用離散元法對臥式旋耕刀作業過程中所受三項阻力及扭矩進行了仿真分析，通過后期土槽試驗對比發現，仿真得到的數據與試驗得到的數據變化趨勢相同；鄧佳玉[9]使用EDEM對深松鏟作業過程的受力情況進行仿真分析，結果表明仿真值與試驗值相差在5%～15%；Ucgul等[10]通過離散元法模擬了板形犁耕作過程中的受力情況，得到了離散元法相較于解析法更加精確的結論；Makange等[11]利用EDEM軟件建立了土壤—農具間的交互模型，該模型可以精確的模擬板形犁的切削力，由此可見，眾多學者的研究都表明離散元法于農業機械作業過程研究的適用性。

本文旨在利用離散元法研究立式旋耕刀在不同前進速度、刀具轉速、耕深條件下，刀具所受阻力、扭矩及其作業功率和其碎土效果；設計正交試驗得到立式旋耕作業參數的系列預測模型，并利用所得預測模型綜合旋耕機作業功率與碎土率，優化立式旋刀的工作參數。

1 立式旋耕作業模型建立

1.1 立式旋耕刀建模

立式旋耕刀是立式旋耕的作業部件，其結構參數對作業效果及機具整體性能影響很大，本試驗采用的旋耕刀其參數的設置如表1所示。

使用SolidWorks2018軟件建立立式旋耕刀的三維模型如圖1所示，將模型文件保存為.igs格式。

表1 立式旋耕刀的參數設置Tab. 1 Parameter setting of vertical rotary tiller blade

圖1 立式旋耕刀的三維模型

1.2 土槽模型建立

1.2.1 土壤顆粒模型的建立

為探究土壤顆粒的性質，于膠東地區農田采集土樣，通過電子顯微鏡觀察發現土壤顆粒呈現不規則的球形及三角片狀，因此本試驗采用圖2所示的球體和三球組合顆粒模擬真實土壤顆粒，為減少仿真運算量，將球體顆粒和三球組合顆粒中球體半徑放大至5 mm。

圖2 球形顆粒與三球組合顆粒

1.2.2 接觸模型的選擇

EDEM軟件中內置了Hertz-Mindlin(no slip)、Hertz-Mindlin with Bonding、Hertz-Mindlin with JKR等接觸模型，其中Hertz-Mindlin with Bonding、Hertz-Mindlin with JKR都可應用于土壤顆粒間的約束[12-13]。通過土壤質地分析試驗得到本試驗中土壤平均含水率17.94%，宜選用Hertz-Mindlin with Bonding模型模擬土壤顆粒間的粘結[8]。該模型在顆粒之間形成粘結鍵以承受一定的阻力和阻力矩，粘結鍵產生之后，顆粒上的法向力Fn、切向力Ft與法向力矩Tn、切向力矩Tt從0隨時間步按照式(1)調整。

(1)

式中：δ——時間步長；

vn——顆粒法向速度，m/s；

vt——顆粒切向速度，m/s；

Sn——顆粒法向剛度，N/m；

St——顆粒切向剛度，N/m；

ωn——顆粒法向角速度，rad/s；

ωt——顆粒切向角速度，rad/s；

A——顆粒接觸面積，m2；

J——顆粒慣性矩，m4；

RB——顆粒間粘結半徑，m。

粘結鍵所受法向力與切向力達到臨界值時斷裂，其臨界法向、切向應力如式(2)所示。

(2)

式中：Mτ——顆粒間切向力矩；

Fτ——顆粒間切向力；

Mn——顆粒間法向力矩。

1.2.3 參數設置

離散元仿真參數包括材料本征參數、材料接觸參數和接觸模型參數，其中材料本征參數包括材料泊松比、密度、剪切模量(或楊氏模量)；材料接觸參數包括恢復系數、靜摩擦因數、動摩擦因數；Hertz-Mindlin with Bonding接觸模型參數包括法向剛度、切向剛度、臨界法向應力、臨界切向應力及粘結半徑[14-15]。本試驗中立式旋耕刀材料選為65 Mn，通過查閱相關文獻[16-18]，確定離散元仿真的具體參數如表2所示。

在EDEM前處理模塊中創建材料及顆粒，并按表2設置好參數后，建立尺寸為2 000 mm×500 mm×400 mm的土槽，土槽中土壤的高度為27 cm，其中由上至下0～16 cm為耕層土壤，16～27 cm為犁底層土壤，土槽模型如圖3所示。

表2 仿真模型參數設置Tab. 2 Simulation model parameter settings

圖3 土槽模型

2 旋耕作業仿真分析

2.1 仿真過程

設定立式旋耕刀的前進速度為1.5 m/s，刀具轉速300 r/min，耕深25 cm。設定仿真總時長為1.25 s，其中0～0.25 s為刀具入土時間，0.25～1.25 s為旋耕作業時間，仿真過程如圖4所示。

(a) 0.00 s (b) 0.25 s (c) 0.50 s

(d) 0.75 s (e) 1.00 s (f) 1.25 s

2.2 旋耕作業阻力分析

定義立式旋耕刀受到阻力的方向，如圖5所示，設Fx、Fy、Fz分別為立式旋耕刀所受前進阻力、側向阻力與垂直阻力，以刀具前進方向為前進阻力Fx的正向，垂直于前進方向向內為側向阻力Fy的正向，豎直向上為垂直阻力Fz的正向。對0.44～0.66 s立式旋耕刀的阻力進行分析，作旋耕刀所受總阻力與三向阻力值隨時間的變化圖，如圖6所示。

圖5 立式旋耕刀三向阻力方向定義

圖6 旋耕刀所受阻力隨時間的變化

其中0.45～0.65 s為旋耕刀進入作業階段的第二個轉動周期，由圖5可知：旋耕刀前進阻力Fx始終為負值，這表明在整個旋耕周期內旋耕刀的前進阻力與前進方向相反；旋耕刀所受前進阻力與總阻力變化趨勢相同且同時達到峰值和谷值，這表明旋耕作業過程中前進阻力對總阻力的變化趨勢影響較大。側向阻力Fy與前進阻力Fx的變化趨勢相反，前進阻力達到峰值時，側向阻力在谷值附近，前進阻力達到谷值時，側向阻力達到峰值附近。旋耕刀垂直阻力Fz則在整個旋耕作業過程中處于相對穩定的狀態，并在300 N左右上下浮動。

在一個周期內，旋耕刀總阻力的峰值為1 374.79 N、1 460.10 N分別在0.45 s、0.56 s時刻出現，谷值為557.756 N、587.248 N，在0.51 s、0.61 s時刻出現。旋耕刀總阻力達到峰、谷值時所處位置如圖7所示。

(a) 0.45 s (b) 0.56 s

(c) 0.51 s (d) 0.61 s

由圖7(a)、圖7(b)可知，當旋耕刀總阻力為峰值時，其位于機具8刀具前進方向上的未耕區域，土壤質地堅硬，旋耕刀受到切削土壤、破壞土壤結構與土壤顆粒時的阻力；當旋耕刀總阻力為谷值時，其位于如圖7(c)、圖7(d)所示的已耕區域內，土壤質地疏松，旋耕刀只受到擾動土壤顆粒的阻力[18-19]。

2.3 旋耕刀所受扭矩分析

旋耕作業中，旋耕刀受到的扭矩是一項重要的工作指標，通過扭矩的變化情況可以判斷旋耕刀整體工作的穩定性。本試驗中，旋耕刀扭矩隨時間的變化情況如圖8所示。

圖8 旋耕刀所受扭矩隨時間變化情況

從圖8可以看出，0～0.25 s，隨著刀具入土深度的增加，其所受扭矩在波動中略有上升，0.25 s后扭矩驟升，在0.27 s達到峰值421.86 N·m。旋耕作業穩定后，旋耕刀所受扭矩在約200～310 N·m間波動，其周期與變化趨勢與旋耕總阻力相似。

3 正交試驗的設置與分析

3.1 試驗方案與結果

立式旋耕作業過程中機具的前進速度、刀具轉速及耕深影響著刀具自身受力和作業效果，為對上述三個工作參數進行分析優化，以旋耕刀所受阻力、扭矩及作業功率和碎土率為評價指標，使用Design-Expert 12軟件設計多因素正交試驗，試驗的因素水平編碼如表3所示。正交試驗設計方案及試驗結果如表4所示。

表3 仿真試驗因素水平及編碼Tab. 3 Factor levels and coding

表4 正交試驗結果Tab. 4 Orthogonal test results

表4中，旋耕刀所受阻力、扭矩取立式旋耕刀作業穩定后2個旋轉周期的平均值，碎土率取旋耕刀作業相同距離(1.2 m)時Bond鍵的破碎率，作業功率通過式(3)計算得出。

(3)

式中：P——作業功率，kW；

T——扭矩，N·m；

n——刀具轉速，r/min。

3.2 方差分析與回歸模型建立

3.2.1 旋耕刀所受阻力的回歸模型與方差分析

通過Design-Expert 12數據處分析，得到阻力F的擬合回歸方程如式(4)所示。

F=930.62+120.4A-37.84B+135.7C+

32.51AC-49.98C2

(4)

回歸模型的方差分析結果如表5所示。

表5 旋耕刀所受阻力回歸模型的方差分析Tab. 5 Variance analysis of rotary tillage resistance regression model

從表5可以看出，該回歸模型的P值小于0.000 1，表明該模型具有很高的顯著性，模型的失擬向P值與F值均大于0.05，擬合精度較高；模型的多元系數R2為0.925 8，表明模型的相關性較好；信噪比APrecision為18.767 0，遠大于4，說明該模型具有較高的可信度與準確度，通過各參數P值對比發現，前進速度A與耕深C對阻力的影響最為顯著，其次是刀具轉速B，二階耕深C2再次之，前進速度與耕深的交互項AC對阻力的影響最小。

3.2.2 旋耕刀所受扭矩的回歸模型與方差分析

通過數據分析，得到旋耕刀所受扭矩T的回歸方程如式(5)所示。

T=241.13+24.32A-5.68B+25.42C+

8.36AC-11.29C2

(5)

該回歸模型的方差分析結果如表6所示，其P值為0.000 2小于0.05，表明該模型顯著性較高；模型的失擬項P值為0.830 4大于0.05，F值為0.452 7亦大于0.05，表明模型失擬性不顯著，誤差較小，擬合程度高；多元系數R2為0.863 1，說明模型中變量的相關性較好，模型的信噪比APrecision為13.309 5，可靠性較高。通過對參數的P值、F值分析可知，各參數中耕深對扭矩的影響最為顯著，其次是前進速度和二階耕深，而刀具轉速、前進速度與耕深的交互項對扭矩的影響均不顯著(P>0.1)。

表6 旋耕刀所受扭矩回歸模型的方差分析Tab. 6 Variance analysis of rotary torque regression model

3.2.3 作業功率的回歸模型與方差分析

通過數據分析，得到作業功率P的多元擬合回歸方程如式(6)所示。

P=7.57+0.771 3A+0.566 3B+0.787 5C+

0.217 5AB+0.26AC-0.381 9C2

(6)

該回歸模型的方差分析結果如表7所示，該回歸模型的P值為0.000 3遠小于0.05，表明該模型極為顯著；對模型進行失擬項檢驗，其P值為0.783 2，F值為0.505 1均大于0.5則模型失擬性不顯著，擬合程度較高；多元系數R2為0.888 2，信噪比APrecision為12.484 0，表明該回歸模型具有較高的精確度與可靠性。通過對各參數的F、P值分析可知，立式旋耕刀的前進速度、刀具轉速、耕深對作業功率的影響都極為顯著，二階耕深顯著，各交互項則不顯著(P>0.1)。

表7 作業功率回歸模型的方差分析Tab. 7 Analysis of variance of power regression model

3.2.4 碎土率的回歸模型與方差分析

通過數據分析，得到碎土率W的多元擬合回歸方程如式(7)所示。

W=67.15-2.21A+1.74B+8.03C+

2.11BC+1.41A2+1.13C2

(7)

回歸模型的方差分析如表8所示。

表8 碎土率回歸模型的方差分析Tab. 8 Analysis of variance of regression model of broken soil rate

從表8可以看出，該回歸模型的P值小于0.000 1，表明該回歸模型極為顯著；模型的失擬項系數F值與P值皆大于0.05，表明該模型的失擬性不顯著，擬合精度較高；模型的多元系數R2為0.934 7，表明模型具有相關性，信噪比APrecision為15.608 7遠大于4，表明該模型具有較高的可信度。對參數的P值、F值分析得到各參數對碎土率影響的顯著性為耕深>前進速度>刀具轉速>轉速與耕深的交互項BC>二階前進速度A2>二階耕深C2。

3.3 各因素的響應面分析

通過Design-Expert 12軟件對數據的處理，得到工作參數對旋耕刀所受阻力、扭矩及作業功率和碎土率影響的響應面圖如圖9～圖12所示。

圖9 各因素對旋耕刀所受阻力的響應面圖

圖10 各因素對旋耕刀所受扭矩的響應面圖

圖11 各因素對旋耕刀作業功率的響應面圖

圖12 各因素對旋耕刀碎土率的響應面圖

從圖9可知，隨著前進速度、耕深的增加以及刀具轉速的減小，旋耕刀所受扭矩呈現增加的趨勢，在耕深加深的過程中，阻力的增加速率逐漸減小；較淺的耕深與較慢前進速度可得到較小的前進阻力。由圖10可知，扭矩隨各因素的變化情況與阻力相近。由圖11可知，前進速度、刀具轉速、耕深與扭矩都呈正相關；較淺的耕深與較小的前進速度可得到較小的作業功率。由圖12可知，耕深及刀具轉速的增加會帶動碎土率的增加，但增加刀具轉速對碎土率的提升效果并不明顯，而刀具的前進速度與碎土率呈負相關；較深的耕深及較快的刀具轉速會帶來較高的碎土率。

3.4 參數優化

利用Design-Expert軟件的優化模塊，通過所得預測模型對立式旋耕刀的作業功率和碎土率進行綜合優化，其約束及目標函數如式(8)所示。

(8)

得到的最佳工作參數：前進速度為1.2 m/s，刀具轉速為307.141 r/min，耕深為25 cm，此時模型預測的立式旋耕刀作業功率為7.042 kW，碎土率為80.85%，將上述參數帶入仿真試驗，得到模型預測值與仿真值的作業功率相對誤差為3.97%，碎土率相對誤差為3.45%，驗證了預測模型的可用性。

4 結論

1) 建立了立式旋耕刀作業過程的離散元模型，通過仿真試驗對立式旋耕刀作業過程中各向阻力及扭矩的變化情況進行了簡要分析。

2) 通過設計正交試驗得到立式旋耕刀所受阻力、扭矩及作業功率、碎土率隨旋耕刀前進速度、刀具轉速、耕深的變化規律，并建立相應的多元回歸預測模型。

3) 利用所得預測模型對立式旋耕作業參數進行優化，得到的結果為：前進速度1.2 m/s，刀具轉速307.141 r/min，耕深25 cm，此時立式旋耕刀的作業功率為7.042 kW，碎土率為80.85%，經過仿真驗證得到作業功率模型預測值與仿真值的相對誤差為3.97%，碎土率相對誤差為3.45%，這表明了在一定的誤差允許范圍內預測模型的可用性。