犁體曲面離散元仿真試驗與參數優化

何永強，賀俊林，楊作梅

(山西農業大學 工學院，山西 太谷 030801)

0 引言

我國是一個農業大國，耕地面積大，達1.35×108hm2[1],普遍使用鏵式犁耕作，存在犁體阻力大、耗能高的問題[2]。近年來，隨著電動拖拉機的發展與應用[3]，在保證鏵式犁耕作質量的情況下，節能更是一種必然趨勢。為減小犁體阻力，降低能耗，本文運用離散元法對犁體曲面進行仿真優化試驗。

離散元法(Distinct Element Method, DEM)是用來計算散體介質系統(干散體、濕散體或兩相懸浮體)中不連續問題力學行為的一種數值模擬方法[4]，可以通過單元間多種連接方式來體現土壤等多相介質間的不同物理關系，從而有效地模擬土壤的開裂、分離等非連續現象及工作部件與土壤顆粒模型間的相互作用力[5]。張青松等[6]利用離散元法對油菜直播機開溝犁進行了優化與試驗，犁體牽引阻力比優化前減小了11.26%。Tsuji[7]等利用離散元法模擬了推土機鏟板對土壤的開挖和推進過程，得到了土壤顆粒滾動、間歇性坍塌等運動特征。目前，對犁體曲面優化方面的研究主要運用解析法[8]、有限元[9-11]等方法，取得了一定的成果，但對利用離散元仿真試驗優化方面的研究鮮見報道。

為減小犁體耕作阻力，筆者建立了EDEM離散元土壤仿真模型，對不同犁體曲面參數的鏵式犁耕作過程進行仿真，得出犁體阻力動態變化趨勢和最佳的犁體結構參數組合，為犁體的優化提供一種離散元分析方法，且為鏵式犁的優化設計提供參考。

1 犁體曲面設計

為減小犁體阻力，以北方系列BT-30扭柱型熟地犁為例，針對其耕作阻力大、耗能高的問題，參照GB/T14225-2008，設計犁體幅寬b為300mm、最大耕深a為250mm。

1.1 導曲線的繪制

導曲線是確定水平直元線位置的引導曲線，扭柱型犁體曲面大多采用拋物線狀的導曲線。本文中導曲線的直線部分長度S根據耕深取50mm，導曲線的拋物線部分采用包絡線作圖法[12-14]繪制，如圖1所示。導曲線開度l、導曲線高度h、端點切線夾角ω分別為

l=C1b(cosΔε-sinε)

(1)

(2)

(3)

其中，l為導曲線開度；C1為常數；b為幅寬；Δε為導曲線扣垡角；ε為犁鏟安裝角；h為導曲線高度；ω為端點切線夾角。

為了提高犁體碎土能力，減小犁鏵與溝底土壤的摩擦、降低犁體損耗，一般ε的范圍為25°～30°。導曲線扣垡角Δε一般在5°～10°之間，Δε越大扣垡作用越強，可以改善碎土性能，但阻力隨之會增大。C1為常數，一般在1.0～1.7之間。取小值時，碎土作用較強，但翻垡的效果減弱且阻力增大。為了保證耕作質量，本文取C1=1.3。

圖1 導曲線

1.2 元線角變化規律

熟地型犁體曲面元線角變化規律分為兩部分，元線角θ0～θmin的變化規律一般采用直線，所以按照直線法求出各元線角即可。元線角θmin～θmax一般的變化規律為二次曲線的一部分，方程為

(4)

式中x—笛卡爾坐標系中二次曲線的橫坐標；

y—笛卡爾坐標系中二次曲線的縱坐標。

為了方便計算，現將直線部分和二次曲線部分上元線角的變化統一表達為元線角與對應元線號的函數關系，即

(5)

式中θn—元線號為n時的元線角；

θ0—初始元線角(即鏵刃角)，一般取θ0=36°～45°；

n—元線號；

nmin—最小元線號，nmin=0；

nmax—最大元線號，nmax=15；

θmin—最小元線角，θmin=40°；

θmax—最大元線角；θmax=47°；

Δz—元線之間的距離,Δz=25mm。

1.3 三維模型的建立

利用SolidWorks軟件建立犁體三維模型。首先，選擇參考平面分別繪制出鏵刃線和導曲線草圖；然后，根據元線號和與其對應的元線角繪制水平直元線，并利用“邊界-曲面”命令生成光滑的曲面； 最后，繪制出犁體曲面前視圖，將犁體曲面前視圖投影到初始光滑曲面上，得到一個封閉的空間曲線。在此基礎上，利用“裁剪-曲面”命令將犁體曲面裁剪出來，再使用“加厚”命令將犁體曲面變為三維實體，結果如圖2所示。

1.水平直元線 2.犁體實體模型 3.導曲線 4.犁體曲面前視圖 5.鏵刃線。

2 仿真模型的建立

2.1 土壤顆粒間接觸模型

土壤顆粒之間具有一定的黏結作用，故采用Hertz-Mindlin黏結接觸模型[15-16]。當顆粒在某時刻被黏結起來時，按照此模型計算可以得到顆粒之間法向和切向黏結力分別為Fn、Ft，法向和切向黏結力矩分別為Tn、Tt，每個顆粒法向和切向的角速度分別為ωn、ωt，每個顆粒法向和切向的速度分別為vn、vt。其中，Fn、Ft、Tn、Tt隨著時步δt從0開始增加。各物理量之間的關系方程為

(6)

式中RB—顆粒間的黏結半徑；

Sn—法向剛度；

St—切向剛度；

t—時間。

此模型可以阻止顆粒間切向和法向的相對運動。當法向和切向應力超過定義的臨界值時，黏結關系就會被破壞，顆粒間法向和切向應力臨界值滿足的關系式為

(7)

式中σmax—法向黏結力定義臨界值；

τmax—切向黏結力定義臨界值。

為了在數學模型上體現顆粒之間的自然接觸關系，接觸半徑應該設置大于自然中粒子實際接觸半徑。

2.2 離散元模型的建立

在EDEM軟件中定義土壤和犁體材料屬性，并設置土壤顆粒半徑；然后將犁體三維模型轉化為IGS格式導入EDEM軟件中，賦予其犁體材料屬性。設置犁體與土壤顆粒、土壤顆粒與土壤顆粒之間的相互作用的參數；然后由顆粒工廠生成土壤顆粒，設置犁體運動速度；最后，進行仿真試驗研究。仿真模型如圖3所示，參數設置如表1所示。

1.犁體 2.土壤

參數單位數值犁體密度ρkg/m37800犁體剪切模量GPa7.0×1010

續表1

3 仿真試驗與結果分析

3.1 正交試驗設計

根據水平直元線法原理和式(1)～式(5)中各參數之間的關系可知：導曲線鏵刃角、犁鏟安裝角和導曲線扣壟角不同，犁體曲面形狀就不同，可能對犁體阻力有較大影響。

為減小犁體阻力，分析了不同犁體曲面參數對犁體阻力的影響。以耕作穩定時犁體所受阻力平均值Y為試驗指標，選取鏵刃角θ0、犁鏟安裝角ε、導曲線扣垡角Δε為試驗因素，運用EDEM軟件進行L9(34)正交試驗。為了保證犁體對土壤的滑切作用，減小犁體阻力，參照文中1.1、1.2節對以上三因素進行水平選取，如表2所示。仿真試驗中，犁體作業速度設為5km/h、耕深220mm，每組試驗重復3次，取平均值。

表2 試驗因素水平

3.2 試驗結果分析

3.2.1 犁體阻力動態變化規律

犁體耕作阻力隨時間的變化曲線如圖4所示。其中，AB時間段內，鏵式犁鏵尖開始進入土壤，隨著時間變化，犁體所受阻力逐漸增加，耕作過程如圖5(a)所示；之后，犁體全部進入耕作狀態后，犁體所受阻力處于穩定階段，如BC時間段所示，耕作過程如圖5(b)所示。

圖4 犁體阻力

θ0=36°、ε=25°、Δε=5°時的犁體阻力曲線如圖4所示。

圖5 犁體耕作過程

3.2.2 正交試驗結果與分析

正交試驗結果如表3所示。利用SAS軟件對正交試驗結果進行方差統計分析，結果如表4所示。

表3 正交試驗方案與結果

表4 正交試驗方差分析

對試驗結果進行極差分析，確定各因素主次影響順序為A>C>B，其優組合為A3B1C1。模型的顯著性P值為0.004 0，決定系數R2=0.998 7，方差分析結果可靠。方差分析可知：鏵刃角對犁體阻力影響極顯著(P<0.01),犁鏟安裝角、導曲線扣垡角對犁體阻力影響顯著(P<0.05)。所以，在設計制造犁體曲面的過程中，應著重考慮鏵刃角、犁鏟安裝角和導曲線扣垡角對犁體曲面的影響，以達到減小犁體阻力的目的。

選取最優水平組合A3B1C1進行試驗，結果表明：當犁體完全進入耕作狀態時，平均阻力為2 532N，比正交試驗中9組試驗結果的平均值減小了6.36%；鏵刃角為45°時，犁體對土壤的推力方向在土壤顆粒與鏵刃的摩擦錐外，且此時會產生剪應力最大斷面，有利于犁體對土壤的滑切作用，從而減小犁體阻力；犁鏟安裝角越小，犁鏵與溝底土壤的摩擦力越小，導曲線扣垡角越小，犁翼與土垡的摩擦作用越弱。理論分析與仿真優化結果一致，說明優化結果是有效可行的。

4 結論

1)采用水平直元線法設計了犁體曲面，運用SolidWorks軟件對犁體曲面三維實體建模，且運用EDEM軟件建立了犁體—土壤離散元模型。

2)仿真試驗結果表明：犁體阻力隨犁體與土壤接觸面積的增大而增大，當犁體全部進入耕作狀態時，犁體所受阻力達到2 621～2 795N，且趨于穩定。

3)正交試驗表明：鏵刃角對犁體阻力影響極顯著(P<0.01)，犁鏟安裝角和導曲線扣垡角對犁體阻力影響顯著(P<0.05)，且犁體結構最佳參數組合為鏵刃角45°、犁鏟安裝角25°、導曲線扣垡角5°。該組合犁體阻力為2 532N，較9組試驗結果的平均值減小了6.36%。研究結果為犁體曲面優化設計提供了一種離散元分析方法，并為犁體曲面的設計提供了參考。