基于響應(yīng)曲面法的土壤離散元模型的參數(shù)標(biāo)定研究*

劉坤宇，蘇宏杰，李飛宇，焦巍

(1. 中國農(nóng)業(yè)科學(xué)院草原研究所/農(nóng)業(yè)農(nóng)村部草原畜牧業(yè)裝備科學(xué)觀測實驗站，呼和浩特市，010010；2. 內(nèi)蒙古農(nóng)業(yè)大學(xué)機電工程學(xué)院，呼和浩特市，010010)

0 引言

離散元方法是分析和求解一些復(fù)雜的離散系統(tǒng)的動力學(xué)問題所提出的一種數(shù)值計算方法[1]，1971年由Peter Cundall首次提出[2]，并應(yīng)用于巖土力學(xué)。土壤作為一種典型的離散系統(tǒng)，在農(nóng)業(yè)機械耕作中與土壤接觸時，土壤—機械相互作用模型的建立非常困難。而在離散元中，可以將土壤看作為由大量離散的獨立運動顆粒組成的整體。離散元法能夠從微觀的角度出發(fā)，在研究農(nóng)業(yè)機械對土壤的擾動量，土壤對農(nóng)業(yè)機械的作用力等研究時能夠更加精確，可以直觀的反應(yīng)每個顆粒所受的力的大小、位移等[3-5]。

現(xiàn)如今離散元法在農(nóng)業(yè)機械中的應(yīng)用越來越廣泛，李博，王燕[6-7]采用離散元方法，對深松鏟進行仿真，分析深松鏟在耕作中所受的阻力以及松土效果等研究，但在仿真分析過程中土壤的力學(xué)參數(shù)對仿真試驗結(jié)果影響較大。

為解決仿真土壤力學(xué)參數(shù)對仿真結(jié)果的影響[8-11]，針對內(nèi)蒙古呼和浩特地區(qū)沙壤土進行土壤標(biāo)定。采用Hertz-Mindlin with JKR接觸模型，以真實的土壤直剪切試驗，堆積試驗來獲取土壤的泊松比，以及堆積角。以堆積角為響應(yīng)值，基于響應(yīng)面優(yōu)化，標(biāo)定了土壤離散元的參數(shù)。

采用實測試驗與離散元仿真模擬相結(jié)合的方法，對研究土壤與農(nóng)機具接觸部件的影響有著至關(guān)重要的意義，如深松機，移栽機，旋耕機等農(nóng)業(yè)裝備，在微觀角度可以更直接的分析土壤與農(nóng)機具相互作用的運動規(guī)律。在離散元仿真過程中，土壤的離散元模型建立的準(zhǔn)確度，對仿真結(jié)果影響巨大，所以準(zhǔn)確地建立土壤的離散元模型有著重要的研究意義。

1 材料與方法

1.1 土壤基本參數(shù)

試驗時間：2019年10月17日。試驗地點：內(nèi)蒙古呼和浩特市耕地試驗田，采用5點取樣法，取0～50 cm 土壤，將土壤帶回實驗室，進行相關(guān)參數(shù)測量，為后續(xù)深松機，移栽機等農(nóng)業(yè)機械離散元仿真提供支持。本次試驗用烘干法測得土壤含水率為11.62%，測得容重為1.36 g/m3。土壤密度為1 165 kg/m2。使用分級篩對土壤進行顆粒分級，用電子秤稱取不同粒徑對應(yīng)質(zhì)量，計算得到顆粒粒徑的對應(yīng)質(zhì)量分?jǐn)?shù)，如表1所示。

表1 顆粒粒徑分布Tab. 1 Particle size distribution

1.2 土壤直剪試驗

試驗采用BZJ-2型應(yīng)變控制式直剪儀，如圖1所示，測試時對土壤施加4種不同壓力，分別施加100 kPa、200 kPa、300 kPa、400 kPa的垂直壓力進行剪切，每次轉(zhuǎn)動手輪一圈，記錄量力環(huán)的讀數(shù)，使得土樣在3～5 min剪破。當(dāng)土樣剪切結(jié)束后，取走砝碼、透水石，然后清掃剪切盒里面的土樣。旋轉(zhuǎn)手輪使其百分表歸零。每組試驗重復(fù)3次，記錄試驗結(jié)果。

圖1 BZJ-2型應(yīng)變控制式直剪儀

按式(1)、式(2)計算土壤剪應(yīng)力，可得剪應(yīng)力(抗剪強度)測試結(jié)果如表2。

σ=P/S

(1)

式中：σ——剪切面法向垂直應(yīng)力，kPa；

S——土壤的剪切截面面積，m2。

根據(jù)庫倫剪切公式

τ=C+σ×tanφ

(2)

式中：τ——土壤剪應(yīng)力，kPa；

C——土體粘聚力，kPa；

σ——剪切面的法向應(yīng)力，kPa；

φ——土壤內(nèi)摩擦角，(°)。

繪制抗剪強度與垂直應(yīng)力的關(guān)系曲線如圖2所示。根據(jù)土壤的剪切強度與垂直壓力關(guān)系，獲取土壤的內(nèi)摩擦角19°，土壤的內(nèi)聚力9.06，摩擦系數(shù)0.344，根據(jù)材料力學(xué)可以求得土壤的泊松比如式(3)所示[12]。

表2 剪應(yīng)力(抗剪強度)測試結(jié)果Tab. 2 Shear stress (shear strength) test results

圖2 抗剪強度與垂直應(yīng)力的關(guān)系曲線

(3)

其中：K0=1-sinφ。

通過試驗得到內(nèi)摩擦角為19°，計算得到土壤泊松比為0.40。根據(jù)相關(guān)文獻[13-14]與樣品土壤的特性，選用剪切模量為1.2×103kPa，土壤泊松比一般為0.25～0.45之間，因此可選用實際試驗計算得到的泊松比。

1.3 土壤堆積角試驗

土壤的堆積角是農(nóng)機具—土壤相互作用過程中影響土壤應(yīng)變的重要參數(shù)，堆積角試驗選用FT-104B型堆積角測定儀，如圖3所示，測試時，將土壤倒入漏斗內(nèi)，使土壤受重力自然降落并堆積，土壤的堆積斜面與水平面的夾角就是堆積角，重復(fù)5次試驗。

圖3 FT-104B型堆積角測定儀

為防止肉眼觀測帶來的數(shù)據(jù)誤差，將堆積角利用Origin軟件讀取顆粒單側(cè)圖像，通過獲取顆粒邊界點，進行線性擬合，擬合曲線的斜率即為堆積角。具體圖像拾取與擬合如圖4所示。最終取得平均值為43.54°。

(a) 獲取顆粒邊界點

2 離散元虛擬標(biāo)定

2.1 土壤顆粒模型

在Solidworks軟件中建立堆積角試驗裝置，并將其導(dǎo)入EDEM中，試驗裝置如圖5所示。

圖5 仿真試驗裝置圖

模型中，漏斗頂部直徑為125 mm，底部直徑為25 mm。顆粒在漏斗上方，顆粒工廠內(nèi)生成顆粒，然后自由下落，經(jīng)過漏斗底部，落入接料板上。顆粒的生成方式為Dynamic，總生成顆粒重量為500 g。仿真完成后，測量顆粒堆積角度。并進行響應(yīng)面試驗。

在離散元(EDEM 2018)模擬過程中，土壤模型的準(zhǔn)確性是離散元仿真模擬的基礎(chǔ)。為提高土壤模型的準(zhǔn)確性，縮短仿真時間，采用離散元中自帶的球形顆粒代替土壤顆粒形狀，選用標(biāo)準(zhǔn)球形顆粒，顆粒半徑設(shè)置為1 mm。設(shè)置瑞利時間步長為25%。

2.2 離散元仿真基本參數(shù)

內(nèi)蒙古呼和浩特地區(qū)土壤類型為粘壤土，土壤具有一定粘性，采用EDEM中自帶Hertz-Mindlin with JKR接觸模型，顆粒之間相互吸引力用顆粒表面能表示[15]。JKR模型適用于顆粒間有明顯粘性和團聚力的物料[16]。土壤堆積試驗中，試驗裝置材料為鋼材，設(shè)置密度為7 850 kg/m3，泊松比為0.3，剪切模量為7.0×107kPa。

2.3 Plackett-Burman試驗

本次仿真使用EDEM2018版本，在GEMM數(shù)據(jù)庫中輸入實測土壤密度與堆積角度，獲取仿真參數(shù)的范圍值，土壤JKR表面能0～11.25 J/m2，土壤—土壤滾動摩擦0.1～0.2，土壤—土壤靜摩擦0.32～1.16，土壤—土壤恢復(fù)系數(shù)0.15～0.75。通過查閱文獻[13-14]獲得以下參數(shù)范圍，土壤—鋼材恢復(fù)系數(shù)0.2～0.5，土壤—鋼材靜摩擦0.5～1.2，土壤—鋼材滾動摩擦0.05～0.2，土壤泊松比0.2～0.5。

應(yīng)用Design Expert軟件進行Plackett-Burman試驗設(shè)計，選取上述8個真實參數(shù)，與3個虛擬參數(shù)，每個參數(shù)選取高低水平，分別編碼-1和+1表示，如表3所示，共進行12組試驗，每組試驗仿真3次，選取平均值為試驗堆積角，Plackett-Burman試驗設(shè)計與結(jié)果如表4所示。

表3 Plackett-Burman試驗參數(shù)列表Tab. 3 List of Plackett-Burman test parameters

Plackett-Burman試驗設(shè)計與結(jié)果如表4所示。運用Design Expert軟件對PB試驗結(jié)果進行方差分析，得到各參數(shù)的影響效果如表5所示。

表4 Plackett-Burman試驗設(shè)計與結(jié)果Tab. 4 Plackett-Burman test design and results

表5 Plackett-Burman試驗參數(shù)顯著性分析Tab. 5 Significance analysis of Plackett-Burman test parameters

2.4 最陡爬坡試驗

根據(jù)表5可知，X2土壤接觸模型JKR表面能(參數(shù)A)、X3土壤—土壤恢復(fù)系數(shù)(參數(shù)B)、X4土壤—土壤靜摩擦因數(shù)(參數(shù)C)對顆粒堆積角影響顯著，其余參數(shù)影響不顯著，因此針對參數(shù)A、參數(shù)B、參數(shù)C進行最陡爬坡試驗，試驗設(shè)計與結(jié)果如表6所示。

表6 最陡爬坡試驗設(shè)計與結(jié)果Tab. 6 Design and results of the steepest climbing test

由表6可得，堆積角誤差的趨勢是由大減小然后再次增大的，在3號水平時，誤差達到最小值，因此最優(yōu)值區(qū)間在3號水平附近，隨后選取3號水平為中心點，2號、4號水平分別選取為低、高水平進行后續(xù)的Box-Behnken試驗和響應(yīng)面設(shè)計。

2.5 Box-Behnken試驗和響應(yīng)面設(shè)計

2.5.1 Box-Behnken試驗

由最陡爬坡試驗結(jié)果設(shè)計Box-Behnken試驗，表7為Box-Behnken試驗設(shè)計與結(jié)果表。根據(jù)試驗結(jié)果采用Design Expert軟件建立堆積角與變量A、B、C的二階回歸方程

α=45.79+4.36A+1.37B+1.15C-0.672 5AB-

1.50AC+3.74BC-1.70A2-1.44B2+0.640 5C2

式中：α——為土壤堆積角。

表7 Box-Behnken試驗設(shè)計與結(jié)果Tab. 7 Box-Behnken test design and results

對表7進行方差分析，如表8所示。

表8 Box-Behnken試驗方差分析表Tab. 8 Box-Behnken test analysis of variance

由表8可知，土壤接觸模型JKR表面能(參數(shù)A)、土壤—土壤恢復(fù)系數(shù)(參數(shù)B)、土壤—土壤靜摩擦因數(shù)(參數(shù)C)對顆粒堆積角影響顯著，交互作用項BC對堆積角影響顯著。從單因素角度分析，各因素對堆積角的影響順序：土壤接觸模型JKR表面能(參數(shù)A)>土壤—土壤恢復(fù)系數(shù)(參數(shù)B) >土壤—土壤靜摩擦因數(shù)(參數(shù)C)。

2.5.2 回歸模型交互效應(yīng)分析

根據(jù)Box-Behnken試驗方差分析結(jié)果可得，土壤接觸模型JKR表面能(參數(shù)A)—土壤—土壤恢復(fù)系數(shù)(參數(shù)B)的交互項、土壤接觸模型JKR表面能(參數(shù)A)—土壤—土壤靜摩擦因數(shù)(參數(shù)C)的交互項、土壤—土壤恢復(fù)系數(shù)(參數(shù)B)—土壤—土壤靜摩擦因數(shù)(參數(shù)C)的交互項對顆粒堆積角影響顯著。采用Design Expert軟件繪制3個交互作用的響應(yīng)曲面，如圖6所示。

由圖6(a)可以看出，土壤接觸模型JKR表面能(參數(shù)A)與土壤—土壤恢復(fù)系數(shù)(參數(shù)B)變化對應(yīng)的曲面坡度較大，(參數(shù)A)與(參數(shù)B)引起的堆積角變化較大，等高線曲率平緩，表明(參數(shù)A)與(參數(shù)C)交互作用不顯著。

由圖6(b)表示，土壤接觸模型JKR表面能(參數(shù)A)與土壤—土壤靜摩擦因數(shù)(參數(shù)C)對應(yīng)的曲面坡度較大表明(參數(shù)A)與(參數(shù)C)對堆積角的影響較大，圖6(b)的等高線曲率平緩，表明土壤接觸模型JKR表面能(參數(shù)A)—土壤—土壤恢復(fù)系數(shù)(參數(shù)B)交互影響不顯著。

(a) A和B交互作用

由圖6(c)可知土壤—土壤靜摩擦因數(shù)(參數(shù)C)與土壤—土壤恢復(fù)系數(shù)(參數(shù)B)對應(yīng)的曲面坡度較小表明(參數(shù)C)與(參數(shù)B)變化引起的堆積角變化較小，圖中的等高線顯示較大曲率表明(參數(shù)C)與(參數(shù)B)交互作用影響顯著。

3 最優(yōu)參數(shù)組合與仿真驗證

應(yīng)用Design Expert軟件的優(yōu)化功能，以實際測得堆積角43.54°為目標(biāo)，使得仿真結(jié)果最接近實際堆積角43.54°，所求得到若干組解，最終選取與實測堆積角最接近一組最優(yōu)解：土壤接觸模型JKR表面能為3.927 J/m2、X3土壤—土壤恢復(fù)系數(shù)為0.332、土壤—土壤靜摩擦因數(shù)為0.719。如圖7所示為實測堆積角與采用最優(yōu)參數(shù)組合的仿真對比結(jié)果。本次仿真使用最優(yōu)解獲得的參數(shù)，仿真中泊松比、剪切模量為采用文中1.2直剪試驗中得到的參數(shù)，進行3次重復(fù)模擬試驗。得到3次仿真堆積角分別為41.15°、44.20°、41.75°，取得平均值為42.36°，相對誤差2.7%。通過圖7可知，本次試驗獲得的堆積角仿真結(jié)果與真實的試驗結(jié)果在堆積角角度與堆積的形態(tài)上有很高的相似性。

(a) 堆積角仿真結(jié)果

4 結(jié)論

1) 針對內(nèi)蒙古呼和浩特地區(qū)沙壤土，通過土壤堆積角實測試驗，獲取實際堆積角度，采用實測試驗與仿真相結(jié)合的方法，使用Design Expert軟件依次設(shè)計Plackett-Burman試驗、最陡爬坡試驗與Box-Behnken試驗，篩選對堆積角影響顯著的物理參數(shù)，分析影響堆積角參數(shù)的交互作用，最終運用優(yōu)化功能，得到土壤的最優(yōu)參數(shù)組合。將最后參數(shù)輸入EDEM中仿真，得到仿真結(jié)果與實際堆積角試驗對比發(fā)現(xiàn)，堆積角角度與堆積的形態(tài)上有很高的相似性。

2) 由Plackett-Burman試驗結(jié)果可知，土壤接觸模型JKR表面能(參數(shù)A)、土壤—土壤恢復(fù)系數(shù)(參數(shù)B)、土壤—土壤靜摩擦因數(shù)(參數(shù)C)對顆粒堆積角影響顯著，其余參數(shù)影響均不顯著。由Box-Behnken試驗可知，土壤—土壤靜摩擦因數(shù)(參數(shù)C)與土壤—土壤恢復(fù)系數(shù)(參數(shù)B)對堆積角的交互作用影響為顯著。其余交互作用均布顯著。