何永強，賀俊林，楊作梅

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

0 引言

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

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

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

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 元線角變化規(guī)律

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

(4)

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

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

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

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

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

2 仿真模型的建立

2.1 土壤顆粒間接觸模型

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

(6)

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

Sn—法向剛度；

St—切向剛度；

t—時間。

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

(7)

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

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

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

2.2 離散元模型的建立

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

1.犁體 2.土壤

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

續(xù)表1

3 仿真試驗與結果分析

3.1 正交試驗設計

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

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

表2 試驗因素水平

3.2 試驗結果分析

3.2.1 犁體阻力動態(tài)變化規(guī)律

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

圖4 犁體阻力

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

圖5 犁體耕作過程

3.2.2 正交試驗結果與分析

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

表3 正交試驗方案與結果

表4 正交試驗方差分析

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

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

4 結論

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

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

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