植煙沙壤土與觸土部件相互作用的離散元仿真參數標定*

林恒矗，王琪，廖鵬，鄭文鑫，何金成，張德暉

(1. 福建農林大學機電工程學院，福州市，350002； 2. 福建省農業信息感知技術重點實驗室，福州市，350002； 3. 現代農業裝備福建省高校工程研究中心，福州市，350002)

0 引言

近年來，基于離散元法的虛擬仿真技術及專業離散元軟件EDEM在農業機械研究中得到廣泛應用，在土壤工作部件、精量播種、收獲機械等方面取得了一系列重要進展[1]。如王金武等[2]利用EDEM進行胡蘿卜收獲機松土鏟的設計與優化，其設計的松土鏟相較于鑿式松土鏟，前進阻力降低5.79%，胡蘿卜拔取力降低了20.68%；鮑玉冬等[3]進行了機械采收藍莓的數值模擬研究，通過能量變化值分析果實損傷程度，為國內機械采收藍莓的推廣提供了理論支持；宋少龍等[4]進行了分層施肥靴的參數優化，為精準對行分層施肥機的改造提供了參考；杜俊等[5]完成了對于旋埋刀輥功率分配特性的研究，得出不同工作條件下，旋埋刀輥中螺旋橫刀功率占總功率比例在60%左右；趙淑紅等[6]設計了玉米秸稈還田交互式分層深松鏟并進行了離散元仿真，得出交互式分層深松鏟破茬比率高，土壤擾動小，能耗低，滿足深松作業的要求；石志鳴等[7]通過離散元仿真試驗研究了翼型深松鏟鏟間距對土壤蓬松度、擾動系數的影響，仿真試驗與田間試驗平均誤差為11.2%、13.5%。

本研究利用萬能試驗機進行土壤堆積角試驗，而后設計顯著性影響試驗、最陡爬坡試驗、響應曲面試驗，依據試驗設計利用離散元軟件對土壤進行仿真堆積角試驗，得出仿真堆積角數值，并以目標值進行求解，確定了土壤—土壤碰撞恢復系數、靜摩擦系數及動摩擦系數等基本接觸參數，以標定的參數進行仿真堆積角試驗，對比實際堆積角試驗值，驗證標定結果可靠性，在此基礎上，以標定出的數值構建植煙土壤土槽模型，并進行土壤—觸土部件相互作用的離散元仿真試驗，得出觸土部件所受扭矩，后續進行土槽試驗，驗證土槽模型可靠性。

1 植煙土壤堆積角試驗

離散元仿真模型參數主要包括材料本征參數和基本接觸參數，本征參數包括物料密度、顆粒體幾何尺寸、物料剪切模量和泊松比等，基本接觸參數包括顆粒與顆粒、顆粒與材料之間的靜摩擦因數、滾動摩擦因數、碰撞恢復系數等。大部分參數可通過試驗測定得出，而部分參數(尤其是物理接觸參數)很難通過試驗獲取，因而許多學者提出采用實測試驗結合虛擬標定實現仿真參數的精準標定[8]。而堆積角是表征顆粒物料流動、摩擦等特性的宏觀參數，且標定的參數是否具有可信度，可以用仿真堆積角與實際堆積角對照方法驗證[1, 9]，故測量實際休止角值作為參數標定試驗的參照值[8, 10]。

1.1 材料及試驗設備

福建省煙草農業科學研究所煙草種植土壤(含水率29.9%，密度1 698 kg/m3)，島津AG-Xplus微機控制電子萬能試驗機，碳鋼圓筒(內徑55 mm，高250 mm)，碳鋼鋼板(表面光滑平整)，相機。

1.2 土壤堆積角試驗及結果

采用圓筒提升法獲得土壤顆粒堆，測量實際堆積角[9]。試驗前取鋼板水平放置于萬能試驗機平臺，將圓筒垂直放置在鋼板上并倒入土壤樣品320 g，將圓筒上端連接萬能試驗機，隨后啟動萬能試驗機，將圓筒以1 000 mm/min勻速提升，用相機拍下土壤堆積情況，試驗重復10次。

將土壤堆積圖像文件導入Matlab，對圖像進行灰度化、二值化處理，利用Canny算子進行邊緣檢測并取反，得到邊界圖像如圖1所示；將邊界圖像文件導入Origin，利用圖像數字化工具得到邊界圖像的坐標值[11-13]，對單測坐標值進行線性擬合[14]，得到散點圖與擬合直線以及線性擬合方程式，線性擬合度值R2=0.997 7，如圖2所示。對斜率取反三角函數并轉化為角度，即可得到堆積角數值；最終，得到堆積角平均值為34.43°。

圖1 Matlab處理后的邊界圖像Fig. 1 Matlab’s processed boundary image

圖2 線性擬合圖像與方程Fig. 2 Linearly fits the image and equation

2 植煙土壤參數標定仿真試驗

2.1 仿真模型構建

福建煙草種植土壤為沙壤土，土壤顆粒之間存在一定的黏附關系，故選用Hertz-Mindlin with JKR模型，在該模型中，法向彈性接觸力的實現是基于Johnson-Kendall-Roberts理論[15-16]。

(1)

式中：FJKR——JKR法向彈性接觸力，N；

γ——相互作用參數的表面能，N/m2；

E*——等效彈性模量，Pa；

δτ——切向重疊量，m；

R*——等效接觸半徑，m。

(2)

(3)

式中：E——接觸顆粒的彈性模量，Pa；

i、j——顆粒i、j；

v——接觸顆粒的泊松比；

R——接觸顆粒的接觸半徑，m。

切向重疊量δτ與法向重疊量δn的關系

(4)

Hertz-Mindlin with JKR切向彈力計算公式

(5)

(6)

式中：Fτ——切向彈力，N；

Sτ——切向剛度，N/m；

G*——當量剪切模量，Pa。

Hertz-Mindlin with JKR切向阻尼力計算公式

(7)

(8)

式中：Fτd——切向阻尼力，N；

m*——等效質量，kg；

β——常數；

為了簡化運算的同時保證試驗結果的準確性，土壤顆粒采用半徑為1.0 mm的球形顆粒。結合顆粒本征參數的測定及查閱相關文獻[6, 9, 17-22]等，確定了仿真模型參數范圍，如表1所示。

表1 土壤顆粒堆積角仿真模型參數Tab. 1 Soil particle accumulation angle simulation model parameters

2.2 參數標定方法

在EDEM前處理界面中設置顆粒參數，建立顆粒工廠；建立與實際圓筒、鋼板所對應的模型。顆粒生成的初始速度設置為Z軸-2 m/s；仿真試驗的顆粒工廠生成總質量為0.32 kg的顆粒，顆粒生成速度為0.64 kg/s，生成時間為0.5 s，總計生成0.32 kg，顆粒穩定后圓筒以0.166 m/s的速度勻速向上運動，圓筒完全移開且土壤堆積穩定后，獲取堆積角圖像。

利用Design-Expert軟件設計二水平部分因子試驗，在EDEM軟件中展開土壤堆積角仿真試驗，篩選顯著影響堆積角的參數；針對篩選出的參數，進行最陡爬坡試驗，得到最優值所在區域；利用Design-Expert軟件設計Box-Behnken Design試驗，開展土壤堆積角仿真試驗。

2.3 因素顯著性篩選試驗及結果分析

以土壤—土壤、土壤—碳鋼之間的接觸參數為試驗因素，以堆積角仿真試驗值為評價指標，開展二水平部分因子試驗設計。試驗結束后，獲取仿真堆積角圖像，將圖像文件導入Matlab，對圖像進行灰度化、二值化處理，利用Canny算子進行邊緣檢測，對檢測后的邊界進行開運算及取反處理，將邊界圖像文件導入Origin，利用圖像數字化工具得到邊界圖像的坐標值，對單測坐標值進行線性擬合，得到邊界散點圖及其擬合直線，對擬合直線斜率取反三角函數并轉化為角度，即可得到堆積角數值，如表2所示。

對試驗結果進行顯著性分析，如表3所示。由表3可知，根據試驗結果中6個因素對于堆積角的影響程度進行排序，在土壤堆積角試驗中X2、X3與X5對堆積角影響顯著，其余因素影響較小。因此在后續試驗中，非顯著因素取中間水平，即X1取0.375、X4取0.45、X6取0.2。

表2 二水平因子試驗設計與結果Tab. 2 Design and results of the two horizontal factor tests

表3 參數顯著性分析Tab. 3 Parameter significance analysis

2.4 最陡爬坡試驗及結果分析

由表3可知，X3、X5對堆積角的效應為正值，故設置仿真試驗參數值遞增；X2對堆積角的效應為負值，故設置仿真參數遞減。再次進行堆積角仿真試驗，記錄仿真堆積角試驗值、相對誤差，試驗參數及結果如表4所示。

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

在2個試驗因素逐漸增大、1個試驗因素逐漸減小的過程中堆積角隨之增大，與實際物理試驗的堆積角的相對誤差呈先減后增的趨勢，仿真試驗中4號試驗相對誤差最小，為3.96%，故以4號試驗中的各參數數值作為后期試驗的中心點，3號和5號試驗分別作為低水平和高水平進行后續的BBD響應曲面試驗設計。

2.5 響應曲面試驗及結果分析

根據篩選的試驗結果，取X2、X3與X5三個基本接觸參數，利用Design-Expert軟件進行三因素三水平的響應曲面試驗設計，中心水平設置5組重復，共進行17組土壤堆積角仿真試驗，試驗設計及仿真結果如表5所示。

表5 接觸參數響應曲面試驗設計與結果Tab. 5 Contact parameters respond to surface test design and results

應用Design-Expert軟件對響應曲面試驗結果進行擬合分析，選取二次全模型，模型擬合度值R2=0.922，方差分析結果如表6所示，響應曲面分析如圖3所示。由表6知，二次全模型回歸模型極其顯著(P=0.003 9<0.01)，且失擬項不顯著(P=0.097 8>0.05)，證明模型可靠，其中X2、X3和X5的平方項對堆積角的影響都顯著。堆積角仿真值θ與土壤—土壤靜摩擦系數X2、土壤—土壤動摩擦系數X3、土壤—碳鋼靜摩擦系數X5的二次多項式關系為

θ=1.488-94.872X2+195.071X3+141.233X5+28.719X2X3+156.317X2X5+481.657X3X5+30.968X22-385.320X32-791.550X52

(9)

以實際試驗堆積角數值為目標值，利用Design-Expert軟件的目標點求取工具求得，X2、X3、X5的值分別為0.593、0.353、0.267。以此數值再次進行堆積角仿真試驗，重復5次，得出堆積角平均值θ=34.403°，相對誤差δθ=1.08%，堆積角仿真結果與試驗結果接近，說明標定的參數準確可靠。

(a) X2、X3與堆積角曲面圖

(b)X2、X5與堆積角曲面圖

(c)X3、X5與堆積角曲面圖圖3 響應曲面分析圖Fig. 3 Response surface analysis diagram

3 土壤—觸土部件相互作用離散元仿真試驗及結果分析

3.1 觸土部件三維模型

除稈清根觸土部件由耙齒與螺旋葉片兩部分構成，螺旋葉片由一對滿面式葉片與一對帶式葉片組成，如圖4所示。工作時，觸土部件相對于前進方向逆時針旋轉，耙齒將土壤擊散，螺旋葉片將土壤推向兩側，從而減小作業后的土壤不平度。根據煙草根茬直徑、深度及植煙土壟寬、壟高，設計耙齒梁長度為320 mm，帶式葉片外徑為200 mm，內徑為125 mm，滿面式葉片半徑為120 mm，耙齒長度為70 mm。

(a) 觸土部件正視圖 (b) 觸土部件斜視圖圖4 觸土部件三維模型Fig. 4 3D model of soil contact parts1.帶式葉片 2.滿面式葉片 3.耙齒 4.耙齒梁固定座 5.觸土部件軸 6.耙齒梁

3.2 仿真模型

為了更直觀地了解觸土部件對土壤的作業效果，并得到觸土部件所受扭矩，進行土壤—觸土部件相互作用的離散元仿真。由Solidworks構建觸土部件三維模型，導入EDEM軟件中；在EDEM軟件中建立3 000 mm×1 000 mm×700 mm大小的土槽模型，設置兩個顆粒工廠，均生成半徑為5 mm土壤顆粒；顆粒工廠1以質量生成顆粒，共計600 kg，生成速度為300 kg/s；顆粒工廠2以質量生成顆粒，共計146 kg，生成速度為146 kg/s，作為土壟；仿真時間步長設置為Rayleigh時間步長的10%[22]，觸土部件行進速度為0.3 m/s，轉速為6.5 rad/s。

3.3 仿真結果與分析

仿真試驗前后作業效果對比如圖5所示。作業后，中間的土壟被挖掘并被螺旋葉片推向兩側。

(a) 試驗前

(b) 試驗后圖5 仿真試驗Fig. 5 Simulated test

仿真試驗中觸土部件所受到的扭矩值在14.38～63.82 N·m間震蕩，扭矩最大值為63.82 N·m，如圖6所示。

圖6 仿真試驗扭矩值Fig. 6 Torque value of simulation test

4 土槽試驗

4.1 試驗條件

試驗地點：福建農林大學現代農業裝備福建省高校工程研究中心土槽試驗室(土槽長×寬：7.0 m×2.4 m)，試驗設備：土槽車(1.8 m×1.4 m，電機3 kW×2)，變頻器(JF1800G；MITSUBISHI，A700)，扭矩傳感器(CYT-302，精度為0.5%)，采集軟件為M400數據采集管理軟件。

將觸土部件安裝在土槽車，以變頻器調節土槽車行進速度與觸土部件轉速，設定行進速度為0.3 m/s，轉速為6.5 rad/s，傳動方式采用齒數比為25∶60的同步帶輪傳動，傳動效率取0.97。

4.2 試驗方法及指標

調整土槽土壤狀態使其接近于植煙沙壤土狀態，使用變頻器調整刀具轉速及土槽車行進速度，使其與仿真試驗一致，而后進行刀具試驗。作業完成后，獲取扭矩傳感器數據與仿真試驗扭矩數據進行對比分析。

由式(10)～式(12)計算觸土部件所受扭矩。

(10)

(11)

故

(12)

式中：Tr——傳感器所受扭矩，N·m；

T——觸土部件所受扭矩，N·m；

FP——皮帶受力，N；

r——小帶輪半徑，m；

R——大帶輪半徑，m；

η——傳動效率。

4.3 結果分析

土槽試驗前后作業效果對比如圖7所示。作業后，中間的土壟被挖掘并被螺旋葉片推向兩側。對比圖5與圖7可知，仿真試驗與土槽試驗的作業效果一致。

土槽試驗扭矩值與仿真試驗扭矩值對比，如圖12所示。土槽試驗扭矩值圖像與仿真試驗扭矩值圖像存在差異，這是由于土槽試驗土壤并非均勻分布且人工起壟在不同位置的形狀尺寸具有一定的差異性。土槽試驗扭矩值在33.40～56.49 N·m間震蕩，扭矩最大值為56.49 N·m，與仿真試驗值(63.82 N·m)誤差為11.48%，兩者扭矩曲線趨勢大致相同。對比仿真試驗中15個扭矩波峰平均值(45.08 N·m)與土槽試驗中15個扭矩波峰平均值(49.34 N·m)，誤差為9.45%，誤差在可接受范圍內，表明仿真模型中土壤與實際土壤的物理力學特性基本一致，土壤模型可靠。

(a) 試驗前

(b) 試驗后圖7 土槽試驗Fig. 7 Soil tank test

圖8 土槽試驗與仿真試驗扭矩值Fig. 8 Torque value of soil tank test and simulation test

5 結論

1) 通過顯著性仿真試驗、最陡爬坡仿真試驗、響應曲面仿真試驗，得出土壤—土壤碰撞恢復系數、土壤—碳鋼碰撞恢復系數、土壤—碳鋼動摩擦系數對堆積角影響較小，其取值分別為0.375、0.45、0.2；響應曲面仿真試驗優化求解得到：土壤—土壤靜摩擦系數、土壤—土壤動摩擦系數、土壤—碳鋼靜摩擦系數分別為0.593、0.353、0.267；利用標定的參數建立模型，進行5次堆積角仿真試驗，均值θ=34.403°，相對誤差δθ=1.08%，仿真結果與試驗結果接近，說明標定的參數準確可靠。

2) 觸土部件離散元仿真試驗與土槽試驗結果表明，觸土部件所受扭矩仿真值與實測值曲線趨勢大致相同，仿真試驗中觸土部件所受到的扭矩最大值(63.82 N·m)與實際試驗值(56.49 N·m)的誤差為11.48%。仿真試驗中15個扭矩波峰平均值(45.08 N·m)與土槽試驗中15個扭矩波峰平均值(49.34 N·m)的誤差為9.45%，說明土壤模型準確可靠。

本研究以物理試驗為參考，通過堆積角試驗和土槽試驗，確定了準確的南方植煙沙壤土離散元模型，為該類型土壤與觸土部件相互作用的動力學研究提供了技術支持。