基于SPH算法的深松鏟土壤切削過程仿真及試驗研究

曹中華，崔晉波，湛小梅，李亞麗，王圓明，余雪源，劉汶樹，宋樹民

(重慶市農業科學院，重慶 401329)

0 引言

為避免土壤耕作過程中出現的地力下降問題，進行土壤深松是迫切的需要[1]。土壤深松技術對于耕層植被和土粒的結構基本不會有影響，有利于降低能耗及成本[2]。深松鏟是實現深松技術的核心部件，深松機在田間作業的主要能耗來自于深松鏟克服土壤阻力所做的功[3]。目前，深松作業面臨的主要問題是作業時的切削阻力較大，為達到理想的耕作效果，研究深松鏟耕作過程中的切削阻力問題刻不容緩。

本文以MASCHIO PINOCCHIO系列全方位深松機的核心部件深松鏟為研究對象，該深松鏟采用鑿型鏟加雙翼鏟結構，最大耕深分別達到450、550、650mm；深松鏟通過鏟尖入土深松、鏟翼淺松側翻土壤、破土鏟對土壤進行破碎。由于工作過程中，耕深較深，土壤切削阻力較大，能耗較大，所以為了降低能耗，進行土壤切削阻力的研究必不可少。然而，目前主要是針對深松鏟的結構參數進行研究，建立深松鏟受力數學模型，對于研究深松鏟土壤切削過程的研究還比較少，不能準確描述深松鏟工作過程中的受力狀態，更不能準確地建立深松鏟土壤切削系統模型[4-8]。

傳統方法是通過田間試驗數據的理論修正來進行農業機具的優化研究，這種方法是直觀、可靠的[9-10]。

然而，田間試驗是一個復雜的系統工程，需要收到人、物、時間、天氣等條件的制約，因此試驗周期長，而且試驗結果不能有效反應機具—土壤相互作用的過程，具有一定的局限[11-14]。20世紀末以來，研究者逐漸將數值模擬方法引入到農機具與土壤相互作用的研究中，可對農機具和土壤相互作用過程進行仿真，直觀反映土壤的運動軌跡和刀具的受力大小[15-17]。因此，本文基于光滑粒子流體動力學(SPH)方法，對深松鏟土壤切削過程進行仿真，分析深松鏟在切削土壤過程中的運動特征和力學行為，并基于耕整地移動式田間動態試驗臺進行田間試驗以驗證仿真的準確性，揭示深松鏟工作過程中切削阻力的變化規律，為未來進行深松鏟的減阻耐磨特性研究和優化設計提供理論依據，也為其他耕作部件土壤切削過程的研究提供理論基礎。深松鏟主要有機械式和振動式兩種類型，本文研究機械式深松鏟。

1 深松鏟切土動力學仿真流程

深松鏟土壤切削動力學仿真流程如圖1所示。

2 深松鏟土壤切削過程仿真研究

本文以MASCHIO深松機的深松鏟為研究對象，為了使模擬更符合實際情況，運用SPH法對深松鏟土壤切削過程進行仿真分析。

2.1 SPH數值分析方法

SPH是一種基于插值理論的無網格的純拉格朗日方法，實質是將連續體離散成相互作用的粒子，由各粒子分擔連續體的質量、速度等物理量，求解時對各粒子分開求解進而獲得整體信息。其計算不需要依賴網格，因此在解決大變形問題上具有獨到優勢。深松鏟土壤切削過程是一個大變形過程，為了使仿真過程更接近實際情況，本文基于SPH法來仿真深松鏟土壤切削過程。

圖1 深松鏟切土動力學仿真流程Fig.1 The flow chart of soil cutting process by subsoiler

2.2 建立仿真模型

2.2.1 土壤SPH模型

本文所選土壤為川渝地區較為普遍的紫色土壤，參照土工試驗方法標準和相關研究對其進行參數測定[18-21]，測得土壤干密度為ρ=1 760kg/m3，內聚力c=5.2e+4Pa，內摩擦角φ=0.223rad相關參數如表1所示。

表1 土壤參數Table1 Parameters of soil

續表1

在LS-Prepost中建立SPH土壤模型，土壤大小為700mm×450mm×500mm，各方向粒子層數為70、45、50，構建完成后的土壤SPH模型如圖2所示。

圖2 土壤SPH模型Fig.2 The SPH model of soil

2.2.2 深松鏟有限元模型

深松鏟結構復雜，在ANSYS中建模困難，本文運用三維掃描儀逆向進行深松鏟點云構建，然后在CATIA中進行三維實體建模。運用Hypermesh軟件進行深松鏟網格劃分，網格單元尺寸為5mm，共452 759個單元，深松鏟有限元模型如圖3所示。深松鏟的主要參數如表2所示。

圖3 深松鏟有限元模型Fig.3 Finite element model of subsoiler表2 深松鏟參數Table 2 Parameters of subsoiler

參數名稱單位取值深松鏟前進速度m/s1深松鏟旋轉速度rad/s0

續表2

2.2.3 土壤—深松鏟有限元模型

在建立土壤—開溝刀具模型時，為有效避免切削土壤過程中土壤大變形問題造成的網格畸變和負體積，本文采用SPH法建立土壤模型，在LS-PREPOST中直接生成SPH粒子，建立700mm×450mm×500mm的SPH土壤模型。圖4為土壤—深松鏟的有限元模型。

圖4 土壤—深松鏟有限元模型Fig.4 The finite element model of soil-subsoiler

2.3 深松鏟切削土壤過程仿真與分析

運用LS-Prepost軟件查看LS-DYNA971求解結束后生成的d3plot文件。本文通過LS-Prepost軟件可以直觀地看到深松鏟對土壤進行深松的過程，同時可以看到深松鏟的運動情況，以及土壤的運動、破碎和變形情況，深松鏟和土壤的受力情況；而且，可以通過分析深松鏟切削土壤過程中的應變和應力云圖，了解土壤應力和應變的分布情況，可以清晰地看到應力集中的位置；通過提取切削力時程曲線分析深松鏟的切削阻力。

2.3.1 切削過程仿真

從圖4和圖5可以看出深松鏟土壤切削仿真過程的效果。圖4為深松鏟土壤切削有限元模型，也是初始狀態，深松鏟與土壤不接觸，SPH粒子不產生任何變化。圖5中，t=0.2s時刻，深松鏟開始切削土壤，鏟尖最先與土壤接觸，土壤發生破壞，深松鏟繼續切削土壤運動，鏟翼也對土壤進行切削；t=0.4s時刻，深松鏟的鏟尖和鏟翼同時進行土壤切削；t=0.6s時刻，整個深松鏟完全進入土壤，此刻也基本處于切削穩定狀態，深松鏟對土壤擾動范圍也最大；t=0.8s時刻，深松鏟土壤切削過程在設置時間內結束，深松鏟鏟尖切出土壤，深松鏟通過對土壤的抬升、扭曲和翻轉，土壤被切削部分的通透性加強，底層土壤的緊實度增加，該部分SPH粒子排列不均勻，與田間作業情況相吻合。

(a) t=0.2s (b) t=0.4s

(c) t=0.6s (d) t=0.8s圖5 深松鏟切削土壤過程Fig.5 The process of soil cutting by subsoiler

2.3.2 土壤的等效應力、應變情況分析

圖6和圖7反應了深松鏟切削土壤過程土壤的等效應力和等效應變效果。

圖6中，t=0.2s時刻深松鏟鏟尖進入土壤對土壤進行深松，此刻土壤受到最大等效應力為1.361MPa；t=0.4s時刻深松鏟鏟尖和鏟翼同時對土壤進行深松和淺松，此刻最大等效應力為1.469MPa；t=0.6s時刻整個深松鏟都進入土壤，對土壤進行深松、淺松和破土，此刻最大等效應力為3.184MPa；t=0.8s時刻深松鏟鏟尖切除土壤，切削過程結束，此刻土壤最大等效應力為3.025MPa。

通過圖7得出：土壤所受的最大等效應變和最大等效應力分布基本一致，都分布在鏟尖、鏟翼和破土鏟與土壤接觸位置。深松鏟土壤切削整個過程，SPH土壤粒子的擾動也在不斷變化，剛開始接觸時土壤粒子擾動較小。隨著深松鏟的鏟尖、鏟翼和破土鏟逐漸進入土壤， SPH土壤粒子的擾動加大， 且內側擾動明顯大于外側，說明深松鏟內側SPH土壤的應力、應變大于外側。

(a) t=0.2s

(b) t=0.4s

(c) t=0.6s

(d) t=0.8s圖6 土壤等效應力分布圖Fig.6 Equivalent stress distribution map of soil

(a) t=0.2s

(b)t=0.4s

(c)t=0.6s

(d)t=0.8s圖7 土壤等效應變分布圖Fig.7 Equivalent strain distribution map of soil

2.3.3 深松鏟切削阻力分析

深松鏟土壤切削過程中深松鏟受到土壤的反作用力即為切削阻力。深松鏟的切削阻力可以通過在LS-Prepost 中載入rcforc文件進行查看。圖8(a)為x、y、z3個方向的切削阻力；圖8(b)為切削阻力合力Resultant Force的時程曲線。

(a) x、y、z方向切削阻力曲線

(b) 深松鏟切削阻力合力曲線圖8 深松鏟切削阻力時程曲線Fig.8 Curve of tillage resistance of subsoiler

由圖8(a)可以得出：深松鏟切削阻力由x、y、z3個方向的阻力合成，值為負說明所受力為阻力。A曲線的走勢說明了在深松鏟運動方向受到最大的切削阻力；B曲線的走勢說明了在深松鏟切削運動過程中，在豎直方向也會受到較大的切削阻力；C曲線的走勢說明在深松鏟切削土壤運動過程中在兩側會受到土壤的擠壓，產生較小的土壤擠壓力，該方向的力對切削阻力的影響不大。

由圖8(b)可知：初始時刻，深松鏟土壤未進行接觸，切削阻力數值為0；緊接著，深松鏟鏟尖首先切入土壤，切削阻力增大，隨著深松鏟鏟翼和破土鏟逐漸切入，切削阻力逐漸增大；0.6s時刻，深松鏟全部切入土壤，切削阻力達到最大值；0.8s時刻深松鏟離開土壤，切削阻力逐漸降低。結合圖6和圖7土壤等效應力和等效應變變化情況，證明切削阻力和土壤等效應力和等效應變存在關聯性。

0.5s時刻深松鏟完全進入土壤，0.8s時刻深松鏟開始切出土壤，提取圖8中0.5～0.8s切削阻力數據，在Excel中進行計算，得出平均切削阻力為3.65kN。

3 田間試驗驗證

3.1 試驗條件

試驗在重慶市農業科學院農業機械研究所農機中試基地進行，以履帶式拖拉機作為牽引動力。試驗過程中，拖拉機前進速度為3.6km/h，耕深標定為450mm。深松鏟切削阻力測試系統采用耕整地移動式田間動態試驗臺。耕作土壤為川渝地區較為普遍的紫色土壤，其物理特性參數如表3所示。

表3 試驗區土壤物理特性參數Table 3 Soil physical parameters in the test area

該試驗臺由輕型履帶式拖拉機、田間作業機具測試系統及相應測試傳感器組成，主要適用于南方丘陵山區較平坦、開闊地塊，可用于旋耕機、深松機等農機具在不同土壤工況(不同堅實度、不同含水率)下的田間作業性能試驗，可直接測試深松鏟工作過程中x、y、z3個方向的切削阻力。該系統由服務器端平臺和遠端平臺兩大部分組成。其中，遠端平臺能夠實現傳感器數據的實時采集與處理、試驗參數配置、試驗數據文件存儲于上傳及采集模塊故障檢測等多種功能。服務器端平臺由計算機及相應測試應用軟件組成，可以實現實時查看遠端平臺測試試驗運行狀態、試驗結束后的數據文件遠程存儲等功能。

試驗之前應進行傳感器檢查、參數設置、速度標定、耕深標定和模塊檢測等事項。田間試驗過程如圖9所示。

圖9 田間試驗過程Fig.9 Field test process

其他試驗儀器包括TJSD-750堅實度儀、TZS-I水分儀、環刀。

該系統通過實時采集并處理農機具在田間作業時的工作性能參數和拖拉機常規性能參數，可用于研究分析農機具在不同作業工況下的作業性能情況及其與拖拉機能耗之間的關系，為農機具的研發、改進提供可靠的科研數據支撐。

3.2試驗內容與結果分析

本文通過耕整地移動式田間動態試驗臺進行深松鏟切削土壤田間試驗，試驗過程中深松鏟耕作深度為450mm，土壤堅實度為1 021kPa，土壤含水率為21%，土壤容積密度為1 760kg/m3。同時，運用遠端田間作業機具測試系統采集耕作過程中深松鏟切削阻力數據，提取20組數據(見表4)，求取切削阻力平均值為3.542kN，與仿真結果相比誤差為2.9%，驗證了基于SPH法進行深松鏟切削土壤過程的仿真是可行的。

表4 深松鏟切削土壤試驗前進速度-切削阻力表Table 4 Forward speed - cutting resistanceTable of soil cutting process by subsoiler

4 結論與展望

4.1 結論

1)基于SPH算法研究了深松鏟土壤切削過程仿真，建立了深松鏟-土壤有限元模型，直觀展示了深松鏟土壤切削過程，與實際情況基本吻合。

2)深松鏟土壤切削仿真表明：不同時刻、不同方向等效應力變化規律不同。切削過程中，最大等效應力為3.184MPa，土壤所受的最大等效應變和最大等效應力分布基本一致，都分布在鏟尖、鏟翼和破土鏟與土壤接觸位置。

3)在耕深為450mm、深松鏟前進速度為3.6km/h、土壤堅實度為1 021kPa、土壤含水率為21%情況下，進行深松鏟土壤切削仿真，在深松鏟完全進入土壤時求解土壤切削阻力平均值為3.65kN。

4)基于耕整地移動式田間動態試驗臺進行田間試驗，測試深松鏟切削阻力值。在耕作穩定后提取20組數據進行切削阻力平均值求解為3.542kN，與仿真結果對比誤差為3.05%，驗證了基于SPH法進行深松鏟切削土壤過程的仿真是可行的。

4.2 展望

1)通過對深松鏟的結構參數和運動參數的研究，結合數值模擬方法和田間試驗，進行深松鏟的優化設計，以達到進一步減阻降耗的目的。

2)基于流固耦合理論，采用數值模擬方法的耦合來進行深松鏟切削土壤過程的仿真分析。