基于響應曲面法的土壤離散元模型的參數標定研究*

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

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

0 引言

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

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

為解決仿真土壤力學參數對仿真結果的影響[8-11]，針對內蒙古呼和浩特地區沙壤土進行土壤標定。采用Hertz-Mindlin with JKR接觸模型，以真實的土壤直剪切試驗，堆積試驗來獲取土壤的泊松比，以及堆積角。以堆積角為響應值，基于響應面優化，標定了土壤離散元的參數。

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

1 材料與方法

1.1 土壤基本參數

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

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

1.2 土壤直剪試驗

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

圖1 BZJ-2型應變控制式直剪儀

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

σ=P/S

(1)

式中：σ——剪切面法向垂直應力，kPa；

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

根據庫倫剪切公式

τ=C+σ×tanφ

(2)

式中：τ——土壤剪應力，kPa；

C——土體粘聚力，kPa；

σ——剪切面的法向應力，kPa；

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

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

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

圖2 抗剪強度與垂直應力的關系曲線

(3)

其中：K0=1-sinφ。

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

1.3 土壤堆積角試驗

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

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

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

(a) 獲取顆粒邊界點

2 離散元虛擬標定

2.1 土壤顆粒模型

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

圖5 仿真試驗裝置圖

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

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

2.2 離散元仿真基本參數

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

2.3 Plackett-Burman試驗

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

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

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

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

表4 Plackett-Burman試驗設計與結果Tab. 4 Plackett-Burman test design and results

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

2.4 最陡爬坡試驗

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

表6 最陡爬坡試驗設計與結果Tab. 6 Design and results of the steepest climbing test

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

2.5 Box-Behnken試驗和響應面設計

2.5.1 Box-Behnken試驗

由最陡爬坡試驗結果設計Box-Behnken試驗，表7為Box-Behnken試驗設計與結果表。根據試驗結果采用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試驗設計與結果Tab. 7 Box-Behnken test design and results

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

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

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

2.5.2 回歸模型交互效應分析

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

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

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

(a) A和B交互作用

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

3 最優參數組合與仿真驗證

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

(a) 堆積角仿真結果

4 結論

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

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