基于SPH/FEM耦合算法的土壤切削仿真與研究

劉亞超，王秀山，馮 敏，張合虎

(河南農業大學 a.機電工程學院；b.園藝學院，鄭州 450000)

?

基于SPH/FEM耦合算法的土壤切削仿真與研究

劉亞超a，王秀山a，馮 敏b，張合虎a

(河南農業大學 a.機電工程學院；b.園藝學院，鄭州 450000)

ANSYS/LS-DYNA仿真軟件，采用SPH/FEM耦合算法模擬旋轉刀具切削土壤作業過程。在SolidWorks中建立刀具—土壤模型，用ANSYS/LS-DYNA進行前處理，修改K文件建立耦合模型，在LS-PREPOST中查看LS-DYNA971求解結果，對數據進行二次處理，分析作業過程能耗變化及應力分布，為后期優化刀具提供依據。文中土壤本構模型采用MAT_147(MAT_FHWA_SOIL)材料，對材料屈服準則進行詳細闡述，相比標準Mohr-Coulomb屈服準則,修改后的Mohr-Coulomb屈服準則引入了修正系數Ahyp，分析了Ahyp取值對屈服面的影響，得出Ahyp與土壤內聚力和內摩擦角關系；闡述SPH算法和SPH/FEM耦合算法,分別用FEM算法、SPH算法和SPH/FEM耦合算法模擬土壤切削過程。結果表明：在土壤切削前期，網格沒有發生畸變時，3種方法模擬結果相近，隨著有限元網格發生畸變，FEM算法產生了誤差，驗證了SPH/FEM耦合算法在土壤切削仿真過程中的可行性與準確性。

ANSYS/LS-DYNA；SPH/FEM耦合算法；旋轉刀具；土壤切削

0 引言

目前，我國農業生產的除草方式主要為化學除草，長期使用化學除草危害性非常大，導致雜草抗藥性增加、土壤生態環境惡化和各種災害頻發等，這與我國人民日益追求的綠色健康生活形成了突出的矛盾。因此，研究除草機器人代替傳統的人工物理除草是今后的趨勢，國內外的研究均圍繞此方面展開[1-7]。本課題組目前正在研究仿生除草機器人，如圖1所示。仿生除草機通過視頻采集系統采集根部及株間雜草分布圖像并傳輸至圖像數據處理系統，經處理后，將執行指令反饋給除草系統；除草系統根據指令驅動滑塊帶動刀具切入土壤，將雜草根部切斷，并對土壤進行松土。本文研究刀具切削土壤過程中土壤應力變化、擠壓變形、失效破壞及刀具應力分布對預測土壤結構變化、減少切削阻力、節約能耗和優化刀具設計具有重要意義。

隨著計算機技術的不斷發展，基于ANSYS/LS-DYNA的動力學仿真技術也越來越成熟。國內外利用FEM有限元法在模擬土壤切削方面已取得較大進展[8-13]。夏俊芳等[14]利用ANSYS/LS-DYNA對螺旋刀輥進行動力學仿真，所得結果與實際理論值相符。劉曉紅等[15]對振動深松鏟切削土壤進行有限元模擬分析，得出了切削過程功耗變化和等效應力情況。朱留憲[16]利用SPH算法對微耕機旋耕刀進行有限元仿真與優化。丁峻宏等[17]對自適應網格劃分、ALE及SPH等3種方法進行仿真模擬對比，結果表明：3種方法在解決土壤切削時各有特點。

圖1 除草機實驗平臺

近年來，土壤切削仿真模擬大多采用FEM有限元法，或SPH光滑粒子法，此兩種方法都存在一些缺陷：FEM有限元法在模擬大變形時容易發生網格畸變(見圖2)，需對大變形區域進行網格細化，計算時間會大大增加，效率隨之降低；SPH光滑質子法解決了FEM有限元法網格畸變的問題，但相比FEM算法計算效率大打折扣，且SPH算法邊界處理不如FEM算法方便。因此，為了發揮兩者優勢，研究人員[18-21]開始對SPH/FEM耦合算法進行研究，并進行實例仿真模擬，仿真結果與實驗結果相符。

圖2 有限元網格發生畸變

目前，關于SPH/FEM耦合算法模擬土壤切削仿真的參考資料很少，本文基于ANSYS/LS-DYNA采用SPH/FEM耦合算法模擬旋轉刀具切削土壤作業過程，利用LS-PREPOST后處理軟件讀取LS-DYNA971計算結果，進行數據研究和二次處理[22-23]；對刀具切削土壤作業過程進行顯示動力學仿真，分析作業過程中應力分布、功耗變化，為后期優化除草機提供參考和依據。

1 算法簡介

1.1 SPH方法基本理論

傳統的有限元法和有限差分法雖然已經得到廣泛應用，但在某些方面具有局限性，尤其是在處理運動不連續時需要每一步重劃網格以保證網格線在整個求解過程中與不連續相符。但是，這樣會引起很大困難，如在求解過程中需要網格映射，會導致精度降低和計算程序復雜化，同時相應費用將會增加[24]。因此，網格劃分限制了它在高速沖擊侵徹、爆炸效應、磨蝕及材料等方面的應用。

為了克服和消除以上限制，光滑質點動力學法(SPH)應運而生。光滑質點動力學法的離散化不使用單元，而是使用固定質量的可動點，質量固定在質點的坐標系，所以SPH也算是拉格朗日型，其基本方程也是守恒方程和固體材料本構方程，標準單元和SPH節點的拉格朗日代碼非常相似[25]。由于SPH法不使用網格，沒有網格畸變問題，可以處理大變形問題，其可以簡單而精確地實現復雜的本構行為，也適用于材料在高速加載速率下斷裂這個困難問題。

SPH的基礎是插值理論[24]，各質點的相互作用借助于插值函數來描述，利用插值函數給出量場在一點的核心估計值(Kemel Estimate),將連續介質動力學的守恒定律由微分方程形式轉換為積分形式，進而轉換為求和。

SPH方法中質點近似函數為

∏h∫(x)=∫f(y)W(x-y,h)dy

(1)

其中，W是核函數(插值核)；h是光滑長度，光滑長度隨時間和空間變化。

核函數W使用輔助函數θ進行定義，有

(2)

其中，d是空間維數；W(x,h)是尖峰函數，如圖3所示。

圖3 核函數

SPH中最常用的光滑核是3次B樣條，定義為

(3)

其中，C是常量，由空間維數決定。

光滑長度h對計算效率和精度有重要影響，目的是為了保持質點鄰域內有足夠的質點，以確保質點連續變量近似有效，如果采用固定的光滑長度材料膨脹時會導致數值畸變[24]。當質點分離時，光滑長度增加；當質點匯聚時，光滑長度減少，其變化范圍在一定范圍，可表示為

C1·h0

(4)

其中，h0是初始光滑長度；C1、C2是縮放因子；如果C1=C2，h則為固定光滑長度，不隨時間和空間變化。

另外，SPH算法中的邊界處理不如有限元網格成熟，有限元網格模型可以通過對邊界節點進行約束來處理邊界，但在SPH方法中不適合。因此，對于SPH方法，邊界需要進行特別處理，通過添加關鍵字在邊界創建虛擬質點來定義邊界對稱面，如圖4所示。

對靠近邊界處的質點，通過自身映射自動創建虛擬質點，虛擬質點具有質點相同的特性，因此可以對質點產生近似對稱作用[24]。

圖4 SPH邊界處理

1.2 SPH/FEM算法

圖5為SPH/FEM的耦合示意圖。其中，左邊部分是SPH質點，右邊部分為FEM網格。兩個部分通過罰函數約束來實現力學參數傳遞；網格邊界區域和SPH質點接觸，通過接觸類型中的“點-面”固連接觸進行耦合定義[24]。

圖5 SPH/FEM的耦合示意圖

在耦合算法中，耦合界面附近的粒子密度相對單元密度越大，算法的計算精度越精確，通常粒子密度與單元密度相當或稍大時，即可獲得比較理想的計算結果；但是，耦合界面處的區域內總存在一定的應力振蕩。實驗表明：隨著光滑長度的增大，粒子區域內的應力振蕩幅度減小[19]。為了保證耦合計算精確和穩定，文中最小光滑長度取0.2倍粒子間距，最大光滑長度取2.0倍粒子間距。

固連接觸失效準則[26]為

(5)

其中，fn為法向力；fs為剪切力；fn.fail為法向失效力；fs.fail為剪切失效力；m1為法向力指數；m2為剪切力指數。

1.3 FEM、SPH、SPH/FEM耦合3種算法仿真對比

3種算法模擬土壤切削時刀具為首次切削土壤，此時最大的等效應力分別為3.047E-03、3.334E-03、3.339E-03。結果相差不大，證明耦合算法的正確性，如圖6～圖8所示。在前期研究時發現，土壤切削前期，網格沒有發生畸變時，3種方法模擬結果相近，隨著有限元網格發生畸變造成FEM算法產生了誤差。

圖6 FEM算法仿真

圖7 SPH算法仿真

圖8 SPH/FEM耦合算法

2 SPH/FEM耦合模型

2.1 作業原理

除草系統中的滑塊左右張合運動帶動左右刀具相向運動，刀片在電機的帶動下，以平行于地面的狀態高速切入土壤中，將雜草的根部切斷，同時對土壤進行松土。圖9為刀具工作示意圖；圖10為刀具在ADAMS仿真中的工作路徑。

圖9 刀具工作示意圖

圖10 刀具株間除草運動軌跡

2.2 模型建立

為了簡化模型、提高效率，本文只研究一側刀具在土壤中的工作狀況，為了減少刀具入土前的無用時間，建立模型時刀具應緊靠土壤模型但不接觸。

通常ANSYS/LS-DYNA建立復雜模型比較困難，本文先在SolidWorks進行建立0.06×0.06×0.08的刀具-土壤模型，并轉換為Parasolid文件；通過與ANSYS/LS-DYNA無縫連接，將文件導入ANSYS/LS-DYNA進行前處理。定義模型為SOLID164單元，并劃分網格，保存K文件，然后用UltraEdit對K文件進行修改建立耦合模型。圖11為SolidWorks模型，圖12為有限元網格模型，圖13為耦合模型。

圖11 SolidWorks模型

圖12 有限元網格模型

圖13 耦合模型

2.3 模型參數設定

2.3.1 刀具模型

刀具在土壤中運動時幾乎不會發生形變，定義刀具模型材料為剛性體MAT_RIGID，采用自由網格劃分法劃分網格。刀具材料參數如下：

密度/kg·m-3：7 890

彈性模量/Pa：2.1E+11

泊松比：0.3

2.3.2 土壤模型

2.3.2.1 土壤材料

土壤模型參數的選取直接決定了仿真結果的準確性，為了能真實地反映出刀具切削土壤實際情況，選取ANSYS/LS-DYNA中土壤專用材料MAT_147(MAT_FHWA_SOIL)[27-28]。這種材料考慮了土壤的密度、土粒相對密度、水密度、粘塑性、體積彈性模量、剪切模量、內摩擦角、土壤含水率及孔隙水效應等，可以很好地表現出土壤在失效后的特征，材料達到失效后可以及時刪除失效單元。

土壤的失效準則采用修正的Mohr-Coulomb(莫爾-庫倫原則)屈服準則[28]，相比標準Mohr-Coulomb屈服準則，修正的Mohr-Coulomb屈服準則引入了Drucker-Prager修正系數Ahyp，修正后的方程為

(6)

其中，σy為屈服面力；P為壓力；φ為內摩擦角；J2為第二不變應力偏量；K(θ)偏平面角函數；Ahyp為修正系數；C為內聚力。

當Ahyp為0時，修正的Mohr-Coulomb屈服準則及為標準屈服準則,如圖14為標準Mohr-Coulomb屈服準則與修正的Mohr-Coulomb屈服準則線性表示。

圖14 標準Mohr-Coulomb屈服準則與

由圖14可知：修正的 Mohr-Coulomb屈服準則在屈服面力為0時，“脊”部更為平滑，因此更接近真實失效條件。

選擇適當的Ahyp系數對仿真的穩定性很重要，研究人員[28]發現，Ahyp系數選擇與內摩擦角和內聚力可以用一個合理的關系表示，有

(7)

隨著Ahyp增大，σy隨之增大。如果Ahyp較大，則明顯偏離Mohr-Coulomb屈服準則的包絡線。為了保持Mohr-Coulomb屈服準則的包絡線，取Ahyp值為1.0E+2Pa。圖15為Ahyp影響下的屈服面。

圖15 Ahyp影響下的屈服面

2.3.2.2 土壤參數

結合本地土壤情況并實地測量土壤參數(見圖16)，參考Soil Material Model 147[27-28]及《農業機械設計手冊》[29]，土壤的主要參數如表1所示，采用kg-m-s單位制。

圖16 土壤含水率測試

名稱單位數值土壤密度kg/m32350土粒相對密度kg/m32．65體積模量Pa2．8E+07剪切模量Pa1．8E+07內摩擦角rad1．1內聚力Pa6．2E+03含水量0．064AhypPa1．0E+2

2.4 修改K文件

1)刪除大變形區域的拉格朗日有限元網格，相應區域添加SPH粒子及相關Part信息。

2)添加SPH相關關鍵字。添加*CONTROL_SPH用于設定SPH質點計算，添加關鍵字*SECTION_SPH用來定義SPH質點屬性，質點光滑長度CSLH取值范圍在1.05～1.3；小于1是不允許的，大于1.3將會增加計算時間[24]，本文中CSLH取值1.2。

3)定義相關接觸。添加*CONTACT_ERODING_NODES_TO_SURFACE關鍵字用于定義刀具與SPH質點之間的接觸為侵蝕接觸算法;添加*CONTACT_TIED_NODES_TO_SU-RFACE_OFFSET關鍵字用于定義SPH質點與有限元網格的固連接觸。

4)處理邊界約束。考慮土壤與大地接觸的實際情況，添加*BOUNDARY_NON_REFLECTING關鍵字,定義模型中土壤有限網格前后面及底部為無反射界面。SPH無法添加無反射邊界，需添加*BOUNDARY_SPH_SYMMETRY_PLANE來定義其邊界。

5)定義速度載荷曲線。刀具移動速度和旋轉速度在模擬仿真中是恒定不變的，添加載荷曲線*DEFINE_CURV*BOUNDARY_PRESCRIBED_MOTION_RIGID關鍵字定義刀具旋轉角速度。最后保存K文件，提交到ANSYS/LS-DYNA971進行求解。

3 仿真結果分析

3.1 土壤等效應力分析

當刀具開始切削土壤時，土壤先發生形變，隨著刀具切入土壤部分的增加，變形量也隨之增加，直到達到屈服準則，土壤開始被破壞，其應力是先增加后減少最后趨于平穩。圖17為刀具作業過程中土壤應力變化情況。

(a) 運行時間0.00225s (b) 運行時間0.00325s

(c) 運行時間0.00450s (d) 運行時間0.00650s

圖17(a)為刀具剛剛切入土壤，圖17(b)為刀具部分切入土壤，圖17(c)為刀具全部切入土壤，圖17(d)為刀具從土壤中切出。從圖17中可以看出：刀具剛開始切入土壤時等效應力較大，隨著刀具切入土壤部分增多，等效應力趨于平穩。整個切削過程中，土壤受力沒有大的波動，切削過程較為平穩。由于受到土壤壓力及切削力，圖中最大應力均位于刀具內側曲面處，因此在加工刀具時應注意進行局部強化處理。

圖17(d)中，SPH整塊粒子向上翹起，說明刀具切削土壤時對土壤起到了松土的效果，達到設計要求。土壤切削過程中，應力向運動方向擴散，有利于未切削區域土壤松散，減少切削阻力。

3.2 切削能耗分析

刀具切削土壤時，刀具的動能保持不變，切削能耗主要包括為克服土壤變形和土壤運動時需要的能耗。隨著刀具切入土壤，土壤能耗逐漸增加，最后趨于平穩，如圖18所示。

將能耗曲線對時間進行求導，得出功率變化曲線，如圖19所示。

圖18 土壤能耗變化曲線

圖19 土壤功率變化曲線

當刀具切入土壤時，功率迅速增加，隨著刀具切出土壤，功率隨之迅速降低。刀具全部切入時，功率達到最大40W。由此可知，每個高速扭矩電機的輸出功率至少為40W。

4 結論

1)采用SPH/FEM耦合算法模擬刀具在繞自身中心旋轉的情況下水平運動切削土壤的作業過程，在大變形區域采用SPH算法，其余區域采用有限元網格算法。結果表明：耦合算法相比有限網格法可以很好地解決網格畸變的問題，相比SPH算法減少了計算時間，提高了工作效率。

2)FEM在處理大變形時需要格外注意網格畸變問題，SPH算法在處理大變形時不用考慮網格畸變問題，但計算效率不及FEM,在處理邊界問題上也不如其成熟，無法添加無反射邊界，需要用“BOUNDARY_SPH_SYMMETRY_PLAN”來定義邊界。

3)比較分析了FEM、SPH、SPH/FEM耦合3種算法模擬土壤切削過程，在土壤切削前期，網格沒有發生畸變時，3種方法模擬結果相近；隨著有限元網格發生畸變造成FEM算法產生了誤差，從而驗證了SPH/FEM耦合算法在土壤切削仿真過程中的可行性與準確性。

4)仿真得出刀具在切削過程中的土壤受力及應力分布情況，最大應力均位于刀具內側曲面處，為加工刀具及局部強化處理提供了依據。土壤切削過程中，應力波向運動方向擴散，有利于未切削土壤松散破碎，節約切削能耗。整個切削過程中，土壤受力沒有大的波動，切削過程比較平穩。

本實驗基于ANASYS/LS-DYNA動力學仿真軟件，采用SPH/FEM耦合算法模擬旋轉刀具切削土壤，目前處于研究階段，隨后將進行一系列系統仿真并與實驗結果進行對比，并對所建立模型進行修改完善，也為今后同類問題研究提供參考。

[1] 李江國，劉占良，張晉國，等.國內外田間機械除草技術研究現狀[J].農機化研究,2006(10)：14-16.

[2] J K Kouwebhoven.Intra-raw Mechanica Weed Control-possibilities and Problems [J]. Soil and Tillage Research,1997,41:87-104 .

[3] 張朋舉，張紋，陳樹人，等.八爪式株間機械除草裝置虛擬設計與運動仿真木[J].農業機械學報,2010,41(4)：56-59.

[4] Z Gobor，P Schulze Lammers，M Martinov. Development of a mechatronic intra-row weeding system with rotational hoeing tools: Theoretical approach and simulation[J].Computers and Electronics in Agriculture，2013，98：166-174.

[5] 陳勇,田磊,鄭加強.基于直接施藥方法的除草機器人[J].農業機械學報,2005,36(10):91-93.

[6] 梁遠，汪春，張偉，等. 3ZCS-7型復式中耕除草機的設計[J].農機化研究,2010,32(6)：21-24.

[7] 韓豹，申建英，李悅梅.3ZCF-7700型多功能中耕除草機設計與試驗[J].農業工程學報,2011,27(1)：124-129.

[8] 齊龍.基于ANSYS/LS-DYNA的松土刀切削有限元仿真[J].農機化研究，2015，37(7)：48-52.

[9] 馬愛麗，廖慶喜，田波平，等.基于ANSYS/LS-DYNA的螺旋刀具土壤切削的數值模擬[J].華中農業大學學報，2009，28(2)：248-252.

[10] 周明，張國忠，姚興林.有限元法在土壤切削過程模擬中的應用[J].農機化研究，2009，31(9)：190-192.

[11] 周明，張國忠，許綺川，等.土壤直角切削的有限元仿真[J]. 華中農業大學學報，2009，28(4)：491-494.

[12] 劉增漢，楊堅，周仕成，等.基于LS-DYNA的土壤—甘蔗—切割器系統的動力學仿真研究[J].農業裝備與車輛工程，2010(5)：29-31.

[13] 夏哲浩，姚立紅，闞江明.基于ANSYS/LS-DYNA旋轉刀具切削土壤與木材的數值模擬[J].森林工程，2016,32(1)：43-47.

[14] 夏俊芳，何小偉，余水生，等.基于ANSYS/LS-DYNA的螺旋刀輥土壤切削有限元模擬[J].農業工程學報,2013(10)：34-39.

[15] 劉曉紅，邱立春.振動深松鏟土壤切削有限元模擬分析-基于ANSYS/LS-DYNA[J].農機化研究，2017，39(1)：19-23.

[16] 朱留憲.基于SPH算法的微耕機旋耕刀有限元仿真與優化[D].重慶：西南大學，2012.

[17] 丁峻宏，金先龍，郭毅之，等. 土壤切削大變形的三維數值仿真[J]. 農業機械學報， 2007， 38(4)：118-121.

[18] 張志春，強洪夫，高巍然.光滑粒子流體動力學—有限元接觸算法研究[J].高壓物理學報，2011，25(2)：97-102.

[19] 肖毅華.有限元法與光滑粒子法的耦合算法研究[D].長沙：湖南大學，2012.

[20] 吳文峰.基于SPH-FEM耦合方法的微細銑削仿真與實驗研究[D].南京：南京航空航天大學，2013.

[21] 劉春含.有限元和光滑粒子耦合方法在爆炸問題中的應用研究[D].長沙：湖南大學，2013.

[22] 尚曉江，蘇建宇.ANSYS/LS-DYNA 動力分析方法與工程實例[M].北京：中國水利水電出版社，2005：102-116.

[23] 白金澤.LS-DYNA3D理論基礎與實例分析[M].北京：科學出版社，2005：204-231.

[24] 李裕春，時黨用，趙遠. ANSYS+11.0LS-DYNA基礎理論與工程實踐[M].北京：中國水利水電出版社，2005：434-431.

[25] 宿崇，唐亮，侯俊銘，等.基于FEM與SPH耦合算法的金屬切削仿真研究[J].系統仿真學報，2009,21(16)：5002-5005.

[26] 段念，王文珊，于怡青，等.基于FEM與SPH耦合算法的單顆磨粒切削玻璃動態過程仿真[J].中國機械過程，2013,24(20)：2716-2721.

[27] Lewis B A. Manual for LS-DYNA soil material model 147(FHWA-HRT-04-095)[R]. Department of Transportation:Federal,Highway Administration,U.S.A.,2004.

[28] Reid J D,Coon B A,Lewis B A,etc. Evaluation of LS-DYNA soil material model 147(FHWA-HRT-04-095)[R].Department of Transportation: Federal Highway Administration,U.S.A.,2004.

[29] 中國農業機械化科學研究院.農業機械設計手冊(上)[K].北京：中國農業科學技術出版社，2007：237-238.Abstract ID:1003-188X(2017)07-0021-EA

The Soil Cutting Dynamics Simulation and Research Based on SPH/FEM Coupling Algorithm

Liu Yachaoa, Wang Xiushana, Fen Minb, Zhang Hehua

(a.Mechanical and Electrical Engineering College; b.College of Horticulture, Hennan Agricultural University, Zhengzhou 450000,China)

In this paper, based on ANSYS/ LS-DYNA ,The SPH/FEM coupling algorithm are adopted to simulate the rotary cutter cutting soil process,tools-soil model is established in the SolidWorks, Before processing in the ANSYS/LS-DYNA and the mod-el is set up by mod-ifeied the K file,check the LS-DYNA971 results and the datais secondary processing by th-e LS-PREPOST, analysis of process energy consumption ch-ange and stress distribution, provide the basis for later optimization tool;The soil constitutive model adopts MAT_147 (MAT_FHWA_SOIL) materials, the material yieldc-riterion are expounded in details,co-mpare with the standard mohr-coulomb, Ahyp correction coefficient is introduced of The modified mohr-coulomb，analyzed the Ahyp values influence on the yield surface, Ahyp relationship with cohesion and internal friction angle are presented also;Explain SPH algorithm and SPH/FEM coupling algori-thm, using FEM algorithm respectively, SPH algorithm and SPH/FEM coupling algorithm to simulate the soil cutting process,the results show that in the early stage of the soil cutting, gridwithout distortion occurs, three ways of the simulation results, the finite element mesh distort on caused by FEM algorithm error, thus in SPH/FEM coupling algorithm is verified the feasibility and accuracy in the process of soil cutting simulation.

ANSYS/LS-DYNA；SPH/FEM coupling algorithm；rotational cutting tool; soil cutting

2016-08-30

國家自然科學基金項目(U1204524)；河南省高等學校骨干教師資助項目(2014GGJS-034)

劉亞超(1990-),男，河南安陽人，碩士研究生，( E-mail)liuyachao0808@163.com。

王秀山(1975-)，男，鄭州人，副教授，博士，( E-mail)towxs@163.com。

S224.1

A

1003-188X(2017)07-0021-07