基于SPH算法的深松鏟破壞土壤仿真模型

張志軍,高奕玨

(常州大學(xué)環(huán)境與安全工程學(xué)院,江蘇 常州 214500)

1 引言

土壤是影響農(nóng)作物產(chǎn)量的關(guān)鍵因素。當(dāng)前,我國某些地區(qū)土壤出現(xiàn)沙化現(xiàn)象,導(dǎo)致肥力下降,農(nóng)作物減產(chǎn)。此外,長(zhǎng)期使用化學(xué)除草方式,使土壤抗藥性增強(qiáng),生態(tài)環(huán)境惡化。近年來,土壤深松方法逐漸引入到農(nóng)業(yè)生產(chǎn)中,土壤深松方法屬于一種保護(hù)性耕作方法,不會(huì)改變植被與土粒結(jié)構(gòu)。深松鏟是完成該項(xiàng)技術(shù)的主要部件,它的主要能耗是克服土壤阻力做功。由于受阻力與功耗影響較大,還需進(jìn)一步改進(jìn)優(yōu)化,同時(shí)深松作業(yè)涉及范圍廣,通過分析深松鏟切削過程對(duì)土壤造成的影響,利用數(shù)值模擬方式獲得土壤變化情況,為深松鏟優(yōu)化提供依據(jù)。

為實(shí)現(xiàn)深松鏟破壞土壤的仿真,馬躍進(jìn)[1]等人提出基于離散元法的深松鏟減阻效果仿真分析。在深松鏟頂部設(shè)置減阻降耗的凸圓刃,利用非線性黏結(jié)彈性塑性接觸模型構(gòu)建土壤仿真模型,通過插件將顆粒和深松鏟接觸作用力導(dǎo)出,完成深松鏟耕作阻力仿真。該研究可獲取耕作設(shè)備切削阻力,但沒有考慮到不同設(shè)備會(huì)產(chǎn)生不同阻力,導(dǎo)致該模型使用范圍受到限制;丁啟朔[2]等人提出基于田間攝像的深松擾動(dòng)行為和效應(yīng)研究。利用土壤耕作原位綜合測(cè)試系統(tǒng)完成深松實(shí)驗(yàn),從五個(gè)不同方位錄制深松視頻,通過視頻分析土壤緊實(shí)度與失效機(jī)理以及深松效應(yīng)影響。但是由于該方法模型的實(shí)驗(yàn)成本較高,導(dǎo)致該模型無法廣泛應(yīng)用。

為此,本文利用SPH方法構(gòu)建深松鏟破壞土壤仿真模型。SPH方法屬于拉格朗日算法[3],是解決高度非線性問題的有效方式。將此方法與有限元(Finite element method,FEM)相結(jié)合,最大程度體現(xiàn)出兩種算法優(yōu)勢(shì),更加精準(zhǔn)地實(shí)現(xiàn)深松鏟破壞土壤過程仿真。

2 深松耕作模型構(gòu)建

在深松耕作時(shí),梨底層被間隔式破壞,生成相間存在的構(gòu)造,此種構(gòu)造能夠緩解土壤退化現(xiàn)象,有助于土壤保持肥力。還能提高氣體交換能力,確保土壤充分吸收養(yǎng)分,使農(nóng)作物增產(chǎn)。深松鏟對(duì)土壤的破壞原理如圖1所示。

圖1 深松鏟破壞原理示意圖

可將整個(gè)切削過程當(dāng)作土壤受到外界作用力被劃分成不同形狀的過程。有關(guān)研究表明,入土深度不同土壤的失效方式也不同[5]。

田地可當(dāng)作沒有邊界的土壤,但受到仿真模型限制,無法建立較大的土壤模型。本文建立的土壤模型長(zhǎng)寬高分別表示為0.5m、0.3m、1.5m。結(jié)合有關(guān)約束條件仿真實(shí)際狀況,利用SPH方法分析深松鏟破壞土壤的過程。

3 基于耦合分析的深松鏟破壞土壤仿真模型

3.1 SPH+FEM耦合分析模型

SPH的基礎(chǔ)方程為守恒方程[6],能夠有效解決大變形問題,適用于分析土壤變形破壞等非線性問題。差值理論是SPH的根本理論,利用差值函數(shù)描述全部質(zhì)點(diǎn)發(fā)生的相互作用,獲取某點(diǎn)估計(jì)值,將守恒定律變換為積分形式,完成轉(zhuǎn)換求和。

在SPH算法中,質(zhì)點(diǎn)近似函數(shù)公式如下

(1)

式中,W代表核函數(shù),h表示光滑長(zhǎng)度,該值隨時(shí)間與空間的變化而變化。

將核函數(shù)W通過輔助函數(shù)θ表示,獲取尖峰函數(shù)W(x,h)

(2)

式中,d代表空間維數(shù)。

該算法中普遍使用的光滑核表示為

(3)

式中,C為常量,取決于空間維數(shù)多少。u為光滑核系數(shù)。

h的不斷變化會(huì)對(duì)仿真精度產(chǎn)生重要影響,它可以確保存在足夠多的質(zhì)點(diǎn),使質(zhì)點(diǎn)連續(xù)有效。若h為固定值,容易出現(xiàn)數(shù)值畸變現(xiàn)象,光滑長(zhǎng)度的最佳變化范圍規(guī)定為

C1·h0

(4)

式中,h0代表原始光滑長(zhǎng)度,C1與C2表示縮放因子。且C1=C2,h值始終不變。但是SPH算法在邊緣約束處理方面還需改進(jìn),本文將SPH與FEM相結(jié)合,當(dāng)變形區(qū)域較大時(shí)利用SPH方法,如果變形區(qū)域較小則利用FEM方法,確保兩種算法發(fā)揮出最大優(yōu)勢(shì),提高仿真模型精度。

圖2代表兩種算法的耦合模型示意圖。上部分為SPH質(zhì)點(diǎn),下側(cè)區(qū)域是FEM網(wǎng)格。利用懲罰函數(shù)達(dá)到兩個(gè)區(qū)域力學(xué)參數(shù)傳遞的目的。網(wǎng)格邊界部分與SPH質(zhì)點(diǎn)相互接觸,采用點(diǎn)-面形式完成耦合。

圖2 耦合結(jié)構(gòu)示意圖

耦合算法下,粒子密度若比單元密度大,則表明仿真精度較高。為提高仿真模型精度,將0.2倍與2.0粒子間距分別當(dāng)作光滑長(zhǎng)度極小、極大值。

因此失效準(zhǔn)則[7]可定義為

(5)

式中,fn與fs分別代表法向力與剪切力,fn,fail與fs,fail分別表示法向與剪切力指數(shù)。

3.2 深松鏟破壞土壤仿真模型構(gòu)建

3.2.1 深松鏟受力分析

利用上述構(gòu)建的耦合模型分析深松鏟作業(yè)時(shí),產(chǎn)生的阻力情況。

1)鏟尖和鏟柄同時(shí)破壞土壤

深松鏟通常為平移前進(jìn),鏟尖破壞土壤時(shí)會(huì)產(chǎn)生土壤耕作阻力。因?yàn)殓P尖與土壤之間存在的深度較大,所以前進(jìn)過程中受到的力來自土壤垂直方向。鏟柄在作業(yè)時(shí)會(huì)受到土壤水平方向的耕作阻力,假設(shè)兩種阻力分別表示為k1與k2,鏟尖與鏟柄的長(zhǎng)度通過l1、l2描述,則同時(shí)生成的阻力F1的計(jì)算公式如下

F1=k1l1+k2l2

(6)

2)刀面破壞土壤

深松鏟工作時(shí),刀面對(duì)土壤產(chǎn)生擠壓作用,土壤與刀面之間的力包括法向壓力FN與摩擦阻力Ff=μFN,這兩種力在水平方向中的投影為阻力F2

F2=FNsinα+Ffcosα=FN(sinα+μcosα)

(7)

式中,α代表深松鏟入土角度,μ為深松鏟和土壤之間的摩擦系數(shù)[8]。

已有研究顯示,當(dāng)耕地深度在400-500mm區(qū)間內(nèi),土壤的正向壓力情況近似表示為δ0=90kN/m2,因此土壤的正向壓力[9]計(jì)算公式如下

δ=ξδ0

(8)

式中,ξ代表土壤有關(guān)系數(shù)。

假定深松鏟的工作區(qū)域面積表示為S,法向壓力FN的計(jì)算公式如下

FN=Sξδ0

(9)

則得出法相壓力FN與摩擦阻力Ff存在如下三角函數(shù)關(guān)系

Ff=FNtanφ=μFN

(10)

式中,φ代表土壤內(nèi)摩擦角。

將式(9)與式(10)代入到式(7)中,重新獲得刀面破壞土壤生成的阻力F2的計(jì)算公式

F2=Sξδ0(sinα+tanφcosα)

(11)

3.2.2 深松鏟對(duì)土壤的破壞過程

獲取深松鏟與土壤之間的摩擦力后,對(duì)深松鏟破壞土壤的過程完成仿真建模。破壞過程可分為前鏟對(duì)土壤的切削與破碎、后鏟對(duì)土壤的切削與破碎兩個(gè)階段。

1)前鏟破壞過程

前鏟的入土深度小,先與土壤發(fā)生接觸,如果鏟尖表層光滑,沒有磨損,這時(shí)鏟尖作用在無限大土體上,形成作用力P及夾角α。在此過程中,土壤受到的壓縮力σr呈現(xiàn)徑向分布,而垂直方向上表現(xiàn)為自重應(yīng)力σsz。

土壤產(chǎn)生的抗剪性能持續(xù)增大,當(dāng)出現(xiàn)最大值,鏟尖對(duì)土壤生成削切作用,整個(gè)土塊失效破壞。發(fā)生形變的土塊變換成有限體積,這時(shí)僅受到來自沒有失效土壤的作用力[10]。

前鏟逐漸前行,鏟尖使部分土壤提升。在持續(xù)作用力下導(dǎo)致上升的土壤與土體之間形成作用應(yīng)力。在不同力的相互作用下,升起的土壤表層與內(nèi)部出現(xiàn)褶皺現(xiàn)象,生成拉伸應(yīng)力。又在自重力作用下土壤中方向發(fā)生破碎。

2)后鏟破壞過程

由于前鏟先行作用于土壤,促使后鏟作業(yè)區(qū)域土體自重應(yīng)力下降。當(dāng)不存在淺層土壤影響下,后鏟可以提升土壤,確保深松深度,優(yōu)化疏松效果。

3)土壤擾動(dòng)體積

明確深松鏟前、后鏟的破壞過程后,分析土壤擾動(dòng)效果。在綜合分析深松環(huán)境下土壤的不同方向破壞差異,構(gòu)建土壤破壞模型。為便于定量研究,假設(shè)破壞滑移面的螺旋線為直線,并將土壤擾動(dòng)體積表示為V1。

因前、后鏟土壤破裂線完全吻合,則兩鏟的土壤破壞半徑表示為

Rr=d(cotα+cotβ)

(12)

Rf=df(cotα+cotβ)

(13)

式中,Rr與Rf分別代表前后鏟土壤破壞半徑,d與df分別描述前鏟與后鏟的入土深度。

前鏟和后鏟的間距計(jì)算公式如下

Xsp=Rr-Rf=(dr-df)(cotα+cotβ)

(14)

分析上述公式能夠得出,如果前鏟和后鏟的間距發(fā)生改變,則土壤抬升體積隨之改變。假設(shè)前鏟可抬升的土壤體積記為V11,兩鏟之間松動(dòng)的土壤體積為V12,后鏟作用力下提升的土壤體積設(shè)置為V13,它們的表達(dá)式分別如下所示

(15)

V12=dfXspW

(16)

(17)

前鏟與后鏟提升土壤體積V1表示為

(18)

分別分析前后鏟對(duì)土壤的破壞過程,計(jì)算出兩鏟對(duì)土壤的破壞半徑與土壤抬升體積,完成深松鏟破壞土壤仿真模型構(gòu)建。

4 仿真設(shè)計(jì)與結(jié)果分析

為證明本文仿真模型的可行性,選取某地區(qū)的紫色土壤進(jìn)行仿真。土壤的相關(guān)參數(shù)如表1所示。

表1 仿真目標(biāo)參數(shù)表

構(gòu)建SPH土壤模型,土壤形狀為550mm×650mm×500mm。每個(gè)方向的粒子層數(shù)設(shè)置為60、40與45。深松鏟的相關(guān)參數(shù)如表2所示。

表2 深松鏟相關(guān)參數(shù)表

為利用本文方法構(gòu)建的仿真模型對(duì)深松鏟破壞土壤的輪廓狀況進(jìn)行分析,分別獲得深松操作后土壤區(qū)域內(nèi)前鏟與后鏟對(duì)土壤的宏觀破壞輪廓。繪制為橫向破壞輪廓,計(jì)算輪廓寬度。當(dāng)深松鏟入土深度分別為50mm、100mm、150mm、200mm時(shí),前后鏟對(duì)土壤的破壞輪廓分別如圖3-6所示。(圖中,深度表示深松作業(yè)完成后,溝底與未耕地表之間的垂直距離。)

由圖3-圖6可知,當(dāng)入土深度不斷變化時(shí),土壤破壞輪廓存在較大差異。隨著深度增加,破壞輪廓的平均寬度隨之?dāng)U大,土壤擾動(dòng)范圍也持續(xù)增加;前后鏟之間的土壤深度減小。結(jié)合該仿真結(jié)果,若要想改變土壤破壞體積,需在變化入土深度的同時(shí)調(diào)整前鏟與后鏟之間距離。

圖3 入土深度為50mm時(shí)土壤破壞輪廓

圖4 入土深度為100mm時(shí)土壤破壞輪廓

圖5 入土深度為150mm時(shí)土壤破壞輪廓

圖6 入土深度為200mm時(shí)土壤破壞輪廓

(19)

(20)

(21)

Uj=1-Vj

(22)

表3 不同入土深度的指標(biāo)變化表

由表3可知,當(dāng)前鏟的入土深度為100mm時(shí),標(biāo)準(zhǔn)差、變異系數(shù)最小,穩(wěn)定系數(shù)最高。這表明100mm是最佳入土深度,此時(shí)深松穩(wěn)定性較強(qiáng)。這是因?yàn)橥寥雷陨韺儆诜蔷|(zhì)的,且自重應(yīng)力會(huì)隨著入土深度的加深而提高,當(dāng)超過100mm時(shí),鏟尖容易與堅(jiān)硬石塊發(fā)生碰撞,因不能有效克服石塊的自重應(yīng)力,深松鏟會(huì)向上彈起,造成波動(dòng)較大,降低穩(wěn)定性。

經(jīng)過上述仿真,獲取了入土深度對(duì)土壤破壞輪廓的影響規(guī)律,并確定最佳入土深度,為提高深松操作性能提供理論依據(jù)。

5 結(jié)論

土壤介質(zhì)具有復(fù)雜、高度非線性等特征,本文將光滑粒子流體動(dòng)力學(xué)與有限元算法相結(jié)合,構(gòu)建深松鏟破壞土壤的仿真模型。通過該模型能夠得出隨著入土深度的加深,土壤破壞輪廓逐漸增大,同時(shí)確定了最佳入土深度。

但是本次仿真利用的土壤屬于原狀土,內(nèi)部的某些因素可能對(duì)土壤性質(zhì)造成影響,導(dǎo)致實(shí)驗(yàn)存在一定誤差。在今后研究中,需對(duì)土壤樣本做預(yù)處理,確定所有因素對(duì)指標(biāo)產(chǎn)生的影響,使獲得的結(jié)果誤差更小。