植煙沙壤土與觸土部件相互作用的離散元仿真參數(shù)標(biāo)定*

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

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

0 引言

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

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

1 植煙土壤堆積角試驗(yàn)

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

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

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

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

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

將土壤堆積圖像文件導(dǎo)入Matlab，對圖像進(jìn)行灰度化、二值化處理，利用Canny算子進(jìn)行邊緣檢測并取反，得到邊界圖像如圖1所示；將邊界圖像文件導(dǎo)入Origin，利用圖像數(shù)字化工具得到邊界圖像的坐標(biāo)值[11-13]，對單測坐標(biāo)值進(jìn)行線性擬合[14]，得到散點(diǎn)圖與擬合直線以及線性擬合方程式，線性擬合度值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)定仿真試驗(yàn)

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

福建煙草種植土壤為沙壤土，土壤顆粒之間存在一定的黏附關(guān)系，故選用Hertz-Mindlin with JKR模型，在該模型中，法向彈性接觸力的實(shí)現(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切向彈力計(jì)算公式

(5)

(6)

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

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

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

Hertz-Mindlin with JKR切向阻尼力計(jì)算公式

(7)

(8)

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

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

β——常數(shù)；

為了簡化運(yùn)算的同時(shí)保證試驗(yàn)結(jié)果的準(zhǔn)確性，土壤顆粒采用半徑為1.0 mm的球形顆粒。結(jié)合顆粒本征參數(shù)的測定及查閱相關(guān)文獻(xià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ù)，建立顆粒工廠；建立與實(shí)際圓筒、鋼板所對應(yīng)的模型。顆粒生成的初始速度設(shè)置為Z軸-2 m/s；仿真試驗(yàn)的顆粒工廠生成總質(zhì)量為0.32 kg的顆粒，顆粒生成速度為0.64 kg/s，生成時(shí)間為0.5 s，總計(jì)生成0.32 kg，顆粒穩(wěn)定后圓筒以0.166 m/s的速度勻速向上運(yùn)動(dòng)，圓筒完全移開且土壤堆積穩(wěn)定后，獲取堆積角圖像。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(9)

以實(shí)際試驗(yàn)堆積角數(shù)值為目標(biāo)值，利用Design-Expert軟件的目標(biāo)點(diǎn)求取工具求得，X2、X3、X5的值分別為0.593、0.353、0.267。以此數(shù)值再次進(jìn)行堆積角仿真試驗(yàn)，重復(fù)5次，得出堆積角平均值θ=34.403°，相對誤差δθ=1.08%，堆積角仿真結(jié)果與試驗(yàn)結(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 土壤—觸土部件相互作用離散元仿真試驗(yàn)及結(jié)果分析

3.1 觸土部件三維模型

除稈清根觸土部件由耙齒與螺旋葉片兩部分構(gòu)成，螺旋葉片由一對滿面式葉片與一對帶式葉片組成，如圖4所示。工作時(shí)，觸土部件相對于前進(jìn)方向逆時(shí)針旋轉(zhuǎn)，耙齒將土壤擊散，螺旋葉片將土壤推向兩側(cè)，從而減小作業(yè)后的土壤不平度。根據(jù)煙草根茬直徑、深度及植煙土壟寬、壟高，設(shè)計(jì)耙齒梁長度為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è)效果，并得到觸土部件所受扭矩，進(jìn)行土壤—觸土部件相互作用的離散元仿真。由Solidworks構(gòu)建觸土部件三維模型，導(dǎo)入EDEM軟件中；在EDEM軟件中建立3 000 mm×1 000 mm×700 mm大小的土槽模型，設(shè)置兩個(gè)顆粒工廠，均生成半徑為5 mm土壤顆粒；顆粒工廠1以質(zhì)量生成顆粒，共計(jì)600 kg，生成速度為300 kg/s；顆粒工廠2以質(zhì)量生成顆粒，共計(jì)146 kg，生成速度為146 kg/s，作為土壟；仿真時(shí)間步長設(shè)置為Rayleigh時(shí)間步長的10%[22]，觸土部件行進(jìn)速度為0.3 m/s，轉(zhuǎn)速為6.5 rad/s。

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

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

(a) 試驗(yàn)前

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

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

圖6 仿真試驗(yàn)扭矩值Fig. 6 Torque value of simulation test

4 土槽試驗(yàn)

4.1 試驗(yàn)條件

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

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

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

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

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

(10)

(11)

故

(12)

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

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

FP——皮帶受力，N；

r——小帶輪半徑，m；

R——大帶輪半徑，m；

η——傳動(dòng)效率。

4.3 結(jié)果分析

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

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

(a) 試驗(yàn)前

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

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

5 結(jié)論

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

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

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