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

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

0 引言

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

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

1 植煙土壤堆積角試驗

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

1.1 材料及試驗設(shè)備

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

1.2 土壤堆積角試驗及結(jié)果

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

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

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

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

2 植煙土壤參數(shù)標(biāo)定仿真試驗

2.1 仿真模型構(gòu)建

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

(1)

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

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

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

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

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

(2)

(3)

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

i、j——顆粒i、j；

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

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

切向重疊量δτ與法向重疊量δn的關(guān)系

(4)

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

(5)

(6)

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

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

G*——當(dāng)量剪切模量，Pa。

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

(7)

(8)

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

m*——等效質(zhì)量，kg；

β——常數(shù)；

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

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

2.2 參數(shù)標(biāo)定方法

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

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

2.3 因素顯著性篩選試驗及結(jié)果分析

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

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

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

表3 參數(shù)顯著性分析Tab. 3 Parameter significance analysis

2.4 最陡爬坡試驗及結(jié)果分析

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

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

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

2.5 響應(yīng)曲面試驗及結(jié)果分析

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

表5 接觸參數(shù)響應(yīng)曲面試驗設(shè)計與結(jié)果Tab. 5 Contact parameters respond to surface test design and results

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

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

(9)

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

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

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

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

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

3.1 觸土部件三維模型

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

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

3.2 仿真模型

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

3.3 仿真結(jié)果與分析

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

(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 試驗條件

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

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

4.2 試驗方法及指標(biāo)

調(diào)整土槽土壤狀態(tài)使其接近于植煙沙壤土狀態(tài)，使用變頻器調(diào)整刀具轉(zhuǎn)速及土槽車行進速度，使其與仿真試驗一致，而后進行刀具試驗。作業(yè)完成后，獲取扭矩傳感器數(shù)據(jù)與仿真試驗扭矩數(shù)據(jù)進行對比分析。

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

(10)

(11)

故

(12)

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

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

FP——皮帶受力，N；

r——小帶輪半徑，m；

R——大帶輪半徑，m；

η——傳動效率。

4.3 結(jié)果分析

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

土槽試驗扭矩值與仿真試驗扭矩值對比，如圖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%，誤差在可接受范圍內(nèi)，表明仿真模型中土壤與實際土壤的物理力學(xué)特性基本一致，土壤模型可靠。

(a) 試驗前

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

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

5 結(jié)論

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

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

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